def galcenrect_to_XYZ(X,Y,Z,Xsun=1.,Zsun=0.,_extra_rot=True):
    """
    NAME:

       galcenrect_to_XYZ

    PURPOSE:

       transform rectangular Galactocentric to XYZ coordinates (wrt Sun) coordinates

    INPUT:

       X, Y, Z - Galactocentric rectangular coordinates

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

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

    OUTPUT:

       (X, Y, Z)

    HISTORY:

       2011-02-23 - 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 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([[X,X,X],[Y,Y,Y],[Z,Z,Z]]).T).sum(-1)
                 +nu.array([dgc,zero,zero]).T)
    else:
        out= nu.dot(nu.array([[-costheta,0.,-sintheta],
                              [0.,1.,0.],
                              [-nu.sign(Xsun)*sintheta,0.,
                                nu.sign(Xsun)*costheta]]),
                    nu.array([X,Y,Z])).T+nu.array([dgc,0.,0.])
    if _extra_rot:
        return nu.dot(galcen_extra_rot.T,out.T).T    
    else:
        return out