code: def estimateBIsochrone(pot,R,z,phi=None):
      
      if pot is None: #pragma: no cover
          raise IOError("pot= needs to be set to a Potential instance or list thereof")
      if isinstance(R,nu.ndarray):
          if phi is None: phi= [None for r in R]
          bs= nu.array([estimateBIsochrone(pot,R[ii],z[ii],phi=phi[ii],
                                           use_physical=False)
                        for ii in range(len(R))])
          return nu.array([nu.amin(bs[True^nu.isnan(bs)]),
                           nu.median(bs[True^nu.isnan(bs)]),
                           nu.amax(bs[True^nu.isnan(bs)])])
      else:
          r2= R**2.+z**2
          r= math.sqrt(r2)
          dlvcdlr= dvcircdR(pot,r,phi=phi,use_physical=False)/vcirc(pot,r,phi=phi,use_physical=False)*r
          try:
              b= optimize.brentq(lambda x: dlvcdlr-(x/math.sqrt(r2+x**2.)-0.5*r2/(r2+x**2.)),
                                 0.01,100.)
          except: #pragma: no cover
              b= nu.nan
          return b
docs: NAME:
  
         estimateBIsochrone
  
      PURPOSE:
  
         Estimate a good value for the scale of the isochrone potential by matching the slope of the rotation curve
  
      INPUT:
  
         pot- Potential instance or list thereof
  
         R,z - coordinates (if these are arrays, the median estimated delta is returned, i.e., if this is an orbit)
  
         phi= (None) azimuth to use for non-axisymmetric potentials (array if R and z are arrays)
  
      OUTPUT:
  
         b if 1 R,Z given
  
         bmin,bmedian,bmax if multiple R given       
  
      HISTORY:
  
         2013-09-12 - Written - Bovy (IAS)
  
         2016-02-20 - Changed input order to allow physical conversions - Bovy (UofT)
  
         2016-06-28 - Added phi= keyword for non-axisymmetric potential - Bovy (UofT)