code: def cov_pmrapmdec_to_pmllpmbb(cov_pmradec,ra,dec,degree=False,epoch=2000.0):
      
      if len(cov_pmradec.shape) == 3:
          out= sc.zeros(cov_pmradec.shape)
          ndata= out.shape[0]
          lb = radec_to_lb(ra,dec,degree=degree,epoch=epoch)
          for ii in range(ndata):
              out[ii,:,:]= cov_pmradec_to_pmllbb_single(cov_pmradec[ii,:,:],
                                                        ra[ii],dec[ii],lb[ii,1],
                                                        degree,epoch)
          return out
      else:
          l,b = radec_to_lb(ra,dec,degree=degree,epoch=epoch)
          return cov_pmradec_to_pmllbb_single(cov_pmradec,ra,dec,b,degree,epoch)
docs: NAME:
  
         cov_pmrapmdec_to_pmllpmbb
  
      PURPOSE:
  
         propagate the proper motions errors through the rotation from (ra,dec) to (l,b)
  
      INPUT:
  
         covar_pmradec - uncertainty covariance matrix of the proper motion in ra (multplied with cos(dec)) and dec [2,2] or [:,2,2]
  
         ra - right ascension
  
         dec - declination
  
         degree - if True, ra and dec 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:
  
         covar_pmllbb [2,2] or [:,2,2] [pmll here is pmll x cos(b)]
  
      HISTORY:
  
         2010-04-12 - Written - Bovy (NYU)