code: def local_fehdist(feh):
      
      fehdist= 0.8/0.15*np.exp(-0.5*(feh-0.016)**2./0.15**2.)\
          +0.2/0.22*np.exp(-0.5*(feh+0.15)**2./0.22**2.)
  
      return fehdist
docs: feh PDF based on local SDSS distribution
      
      From Jo Bovy:
      https://github.com/jobovy/apogee/blob/master/apogee/util/__init__.py#L3
      2D gaussian fit based on Casagrande (2011)