code: def rvir(self,H=70.,Om=0.3,overdens=200.,wrtcrit=False,ro=None,vo=None,
               use_physical=False): # use_physical necessary bc of pop=False, does nothing inside
          
          if ro is None: ro= self._ro
          if vo is None: vo= self._vo
          if wrtcrit:
              od= overdens/bovy_conversion.dens_in_criticaldens(vo,ro,H=H)
          else:
              od= overdens/bovy_conversion.dens_in_meanmatterdens(vo,ro,
                                                                  H=H,Om=Om)
          dc= 12.*self.dens(self.a,0.,use_physical=False)/od
          x= optimize.brentq(lambda y: (numpy.log(1.+y)-y/(1.+y))/y**3.-1./dc,
                             0.01,100.)
          return x*self.a
docs: NAME:
  
             rvir
  
          PURPOSE:
  
             calculate the virial radius for this density distribution
  
          INPUT:
  
             H= (default: 70) Hubble constant in km/s/Mpc
             
             Om= (default: 0.3) Omega matter
         
             overdens= (200) overdensity which defines the virial radius
  
             wrtcrit= (False) if True, the overdensity is wrt the critical density rather than the mean matter density
  
             ro= distance scale in kpc or as Quantity (default: object-wide, which if not set is 8 kpc))
  
             vo= velocity scale in km/s or as Quantity (default: object-wide, which if not set is 220 km/s))
  
          OUTPUT:
          
             virial radius
          
          HISTORY:
  
             2014-01-29 - Written - Bovy (IAS)