code: def rl(Pot,lz):
      
      Pot= flatten(Pot)
      if _APY_LOADED and isinstance(lz,units.Quantity):
          if hasattr(Pot,'_ro'):
              lz= lz.to(units.km/units.s*units.kpc).value/Pot._vo/Pot._ro
          elif hasattr(Pot[0],'_ro'):
              lz= lz.to(units.km/units.s*units.kpc).value/Pot[0]._vo/Pot[0]._ro
      #Find interval
      rstart= _rlFindStart(math.fabs(lz),#assumes vo=1.
                           math.fabs(lz),
                           Pot)
      try:
          return optimize.brentq(_rlfunc,10.**-5.,rstart,
                                 args=(math.fabs(lz),
                                       Pot),
                                 maxiter=200,disp=False)
      except ValueError: #Probably lz small and starting lz to great
          rlower= _rlFindStart(10.**-5.,
                               math.fabs(lz),
                               Pot,lower=True)
          return optimize.brentq(_rlfunc,rlower,rstart,
                                 args=(math.fabs(lz),
                                       Pot))
docs: NAME:
  
         rl
  
      PURPOSE:
  
         calculate the radius of a circular orbit of Lz
  
      INPUT:
  
         Pot - Potential instance or list thereof
  
         lz - Angular momentum (can be Quantity)
  
      OUTPUT:
  
         radius
  
      HISTORY:
  
         2012-07-30 - Written - Bovy (IAS@MPIA)
  
      NOTE:
  
         seems to take about ~0.5 ms for a Miyamoto-Nagai potential; 
         ~0.75 ms for a MWPotential