def pmllpmbb_to_pmrapmdec(pmll,pmbb,l,b,degree=False,epoch=2000.0):
    """
    NAME:

       pmllpmbb_to_pmrapmdec

    PURPOSE:

       rotate proper motions in (l,b) into proper motions in (ra,dec)

    INPUT:

       pmll - proper motion in l (multplied with cos(b)) [mas/yr]

       pmbb - proper motion in b [mas/yr]

       l - Galactic longitude

       b - Galactic lattitude

       degree - if True, l and b are given in degrees (default=False)

       epoch - epoch of ra,dec (right now only 2000.0 and 1950.0 are supported when not using astropy's transformations internally; when internally using astropy's coordinate transformations, epoch can be None for ICRS, 'JXXXX' for FK5, and 'BXXXX' for FK4)

    OUTPUT:

       (pmra x cos(dec),pmdec), for vector inputs [:,2]

    HISTORY:

       2010-04-07 - Written - Bovy (NYU)

       2014-06-14 - Re-written w/ numpy functions for speed and w/ decorators for beauty - Bovy (IAS)

    """
    theta,dec_ngp,ra_ngp= get_epoch_angles(epoch)
    #Whether to use degrees and scalar input is handled by decorators
    radec = lb_to_radec(l,b,degree=False,epoch=epoch)
    ra= radec[:,0]
    dec= radec[:,1]
    dec[dec == dec_ngp]+= 10.**-16 #deal w/ pole.
    sindec_ngp= nu.sin(dec_ngp)
    cosdec_ngp= nu.cos(dec_ngp)
    sindec= nu.sin(dec)
    cosdec= nu.cos(dec)
    sinrarangp= nu.sin(ra-ra_ngp)
    cosrarangp= nu.cos(ra-ra_ngp)
    #These were replaced by Poleski (2013)'s equivalent form that is better at the poles
    #cosphi= (sindec_ngp-sindec*sinb)/cosdec/cosb
    #sinphi= sinrarangp*cosdec_ngp/cosb
    cosphi= sindec_ngp*cosdec-cosdec_ngp*sindec*cosrarangp
    sinphi= sinrarangp*cosdec_ngp
    norm= nu.sqrt(cosphi**2.+sinphi**2.)
    cosphi/= norm
    sinphi/= norm
    return (nu.array([[cosphi,sinphi],[-sinphi,cosphi]]).T\
                *nu.array([[pmll,pmll],[pmbb,pmbb]]).T).sum(-1)