def galcenrect_to_vxvyvz(vXg,vYg,vZg,vsun=[0.,1.,0.],Xsun=1.,Zsun=0.,
                         _extra_rot=True):
    """
    NAME:

       galcenrect_to_vxvyvz

    PURPOSE:

       transform rectangular Galactocentric coordinates to XYZ coordinates (wrt Sun) for velocities

    INPUT:

       vXg - Galactocentric x-velocity

       vYg - Galactocentric y-velocity

       vZg - Galactocentric z-velocity

       vsun - velocity of the sun in the GC frame ndarray[3] (can be array of same length as vXg; shape [3,N])

       Xsun - cylindrical distance to the GC (can be array of same length as vXg)
       
       Zsun - Sun's height above the midplane (can be array of same length as vXg)

       _extra_rot= (True) if True, perform an extra tiny rotation to align the Galactocentric coordinate frame with astropy's definition

    OUTPUT:

       [:,3]= vx, vy, vz

    HISTORY:

       2011-02-24 - Written - Bovy (NYU)

       2016-05-12 - Edited to properly take into account the Sun's vertical position; dropped Ysun keyword - Bovy (UofT)

       2017-10-24 - Allowed vsun/Xsun/Zsun to be arrays - Bovy (UofT)

       2018-04-18 - Tweaked to be consistent with astropy's Galactocentric frame - Bovy (UofT)

    """
    dgc= nu.sqrt(Xsun**2.+Zsun**2.)
    costheta, sintheta= Xsun/dgc, Zsun/dgc
    if isinstance(Xsun,nu.ndarray):
        zero= nu.zeros(len(Xsun))
        one= nu.ones(len(Xsun))
        Carr= nu.rollaxis(nu.array([[-costheta,zero,-sintheta],
                                    [zero,one,zero],
                                    [-nu.sign(Xsun)*sintheta,zero,
                                      nu.sign(Xsun)*costheta]]),2)
        out= ((Carr
               *nu.array([[vXg-vsun[0],vXg-vsun[0],vXg-vsun[0]],
                          [vYg-vsun[1],vYg-vsun[1],vYg-vsun[1]],
                          [vZg-vsun[2],vZg-vsun[2],vZg-vsun[2]]]).T).sum(-1))
    else:
        out= nu.dot(nu.array([[-costheta,0.,-sintheta],
                              [0.,1.,0.],
                              [-nu.sign(Xsun)*sintheta,0.,
                                nu.sign(Xsun)*costheta]]),
                    nu.array([vXg-vsun[0],vYg-vsun[1],vZg-vsun[2]])).T
    if _extra_rot:
        return nu.dot(galcen_extra_rot.T,out.T).T
    else:
        return out
NAME:

       galcenrect_to_vxvyvz

    PURPOSE:

       transform rectangular Galactocentric coordinates to XYZ coordinates (wrt Sun) for velocities

    INPUT:

       vXg - Galactocentric x-velocity

       vYg - Galactocentric y-velocity

       vZg - Galactocentric z-velocity

       vsun - velocity of the sun in the GC frame ndarray[3] (can be array of same length as vXg; shape [3,N])

       Xsun - cylindrical distance to the GC (can be array of same length as vXg)
       
       Zsun - Sun's height above the midplane (can be array of same length as vXg)

       _extra_rot= (True) if True, perform an extra tiny rotation to align the Galactocentric coordinate frame with astropy's definition

    OUTPUT:

       [:,3]= vx, vy, vz

    HISTORY:

       2011-02-24 - Written - Bovy (NYU)

       2016-05-12 - Edited to properly take into account the Sun's vertical position; dropped Ysun keyword - Bovy (UofT)

       2017-10-24 - Allowed vsun/Xsun/Zsun to be arrays - Bovy (UofT)

       2018-04-18 - Tweaked to be consistent with astropy's Galactocentric frame - Bovy (UofT)