code: def estimate_hr(self,R,z=0.,dR=10.**-8.,**kwargs):
          
          Rs= [R-dR/2.,R+dR/2.]
          if z is None:
              sf= numpy.array([self.surfacemass_z(r,use_physical=False,
                                                  **kwargs) for r in Rs])
          else:
              sf= numpy.array([self.density(r,z,use_physical=False,
                                            **kwargs) for r in Rs])
          lsf= numpy.log(sf)
          return -dR/(lsf[1]-lsf[0])
docs: NAME:
  
             estimate_hr
  
          PURPOSE:
  
             estimate the exponential scale length at R
  
          INPUT:
  
             R - Galactocentric radius (can be Quantity)
  
             z= height (default: 0 pc) (can be Quantity)
  
             dR- range in R to use (can be Quantity)
  
             density kwargs
  
          OUTPUT:
  
             estimated hR
  
          HISTORY:
  
             2012-09-11 - Written - Bovy (IAS)
  
             2013-01-28 - Re-written - Bovy