code: def estimate_hz(self,R,z,dz=10.**-8.,**kwargs):
          
          if z == 0.:
              zs= [z,z+dz]
          else:
              zs= [z-dz/2.,z+dz/2.]
          sf= numpy.array([self.density(R,zz,use_physical=False,
                                        **kwargs) for zz in zs])
          lsf= numpy.log(sf)
          return -dz/(lsf[1]-lsf[0])
docs: NAME:
  
             estimate_hz
  
          PURPOSE:
  
             estimate the exponential scale height at R
  
          INPUT:
  
             R - Galactocentric radius (can be Quantity)
  
             dz - z range to use (can be Quantity)
  
             density kwargs
  
          OUTPUT:
  
             estimated hz
  
          HISTORY:
  
             2012-08-30 - Written - Bovy (IAS)
  
             2013-01-28 - Re-written - Bovy