How to use the galpy.util.bovy_coords function in galpy

MJ= (object-wide default) absolute magnitude in J or an array of samples of absolute  magnitudes in J for the tracer population
           JK= (object-wide default) J-Ks color or an array of samples of the J-Ks color 
        OUTPUT
           effective selection fraction
        HISTORY:
           2017-01-18 - Written - Bovy (UofT/CCA)
        """
        if isinstance(dist,(int,float)):
            dist= numpy.array([dist])
        elif isinstance(dist,list):
            dist= numpy.array(dist)
        MJ, JK= self._parse_mj_jk(MJ,JK)
        distmod= 5.*numpy.log10(dist)+10.
        # Extract the distribution of A_J and A_J-A_Ks at this distance 
        # from the dust map, use twice the radius of a pixel for this
        lcen, bcen= bovy_coords.radec_to_lb(ra,dec,degree=True)
        pixarea, aj= self._dmap3d.dust_vals_disk(\
            lcen,bcen,dist,healpy.nside2resol(_BASE_NSIDE)/numpy.pi*180.)
        totarea= numpy.sum(pixarea)
        ejk= aj*(1.-1./2.5) # Assume AJ/AK = 2.5
        distmod= numpy.tile(distmod,(aj.shape[0],1))
        pixarea= numpy.tile(pixarea,(len(dist),1)).T
        out= numpy.zeros_like(dist)
        for mj,jk in zip(MJ,JK):
            tj= mj+distmod+aj
            tjk= jk+ejk
            out+= numpy.sum(pixarea*self._tgasSel(tj,tjk,ra,dec),axis=0)
        if not self._maxd is None:
            out[dist > self._maxd]= 0.
        return out/totarea/len(MJ)

self._dec_cen_4vol= 90.-180./numpy.pi*theta
            if linearDist:
                dists= numpy.linspace(0.001,10.,ndists)
                dms= 5.*numpy.log10(dists)+10.
                self._deltadm_4vol= dists[1]-dists[0]
            else:
                dms= numpy.linspace(0.,18.,ndists)
                self._deltadm_4vol= (dms[1]-dms[0])*numpy.log(10.)/5.
            self._dists_4vol= 10.**(0.2*dms-2.)
            self._tiled_dists3_4vol= numpy.tile(\
                self._dists_4vol**(3.-linearDist),(len(self._ra_cen_4vol),1))
            self._tiled_ra_cen_4vol= numpy.tile(self._ra_cen_4vol,
                                                 (len(self._dists_4vol),1)).T
            self._tiled_dec_cen_4vol= numpy.tile(self._dec_cen_4vol,
                                                 (len(self._dists_4vol),1)).T
            lb= bovy_coords.radec_to_lb(phi,numpy.pi/2.-theta)
            l= numpy.tile(lb[:,0],(len(self._dists_4vol),1)).T.flatten()
            b= numpy.tile(lb[:,1],(len(self._dists_4vol),1)).T.flatten()
            XYZ_4vol= \
                bovy_coords.lbd_to_XYZ(l,b,
                   numpy.tile(self._dists_4vol,
                              (len(self._ra_cen_4vol),1)).flatten())
            self._X_4vol= numpy.reshape(XYZ_4vol[:,0],(len(self._ra_cen_4vol),
                                                       len(self._dists_4vol)))
            self._Y_4vol= numpy.reshape(XYZ_4vol[:,1],(len(self._ra_cen_4vol),
                                                       len(self._dists_4vol)))
            self._Z_4vol= numpy.reshape(XYZ_4vol[:,2],(len(self._ra_cen_4vol),
                                                       len(self._dists_4vol)))
        # Cache effective-selection function
        MJ, JK= self._parse_mj_jk(MJ,JK)
        new_hash= hashlib.md5(numpy.array([MJ,JK])).hexdigest()
        if not hasattr(self,'_vol_MJ_hash') or new_hash != self._vol_MJ_hash \