Merge 10b8287 into 9dd77fe

py/desiutil/dust.py

@@ -14,7 +14,7 @@
 from astropy.io.fits import getdata
 from astropy.coordinates import SkyCoord
 from astropy import units as u
-from .log import get_logger
+from desiutil.log import get_logger
 log = get_logger()
 
 def extinction_total_to_selective_ratio(band, photsys) :
@@ -181,6 +181,162 @@ def ext_ccm(wave, Rv=3.1):
 
     return A
 
+def ext_fitzpatrick(wave, R_V=3.1, avglmc=False, lmc2=False,
+                    x0=None,gamma=None,
+                    c1=None,c2=None,c3=None,c4=None) :
+    """
+    Return extinction curve from Fitzpatrick (1999).
+
+    Args:
+        wave : 1D array of vacuum wavelength [Angstroms]
+
+    Returns:
+        1D array of A(lambda)/A(V)
+
+    OPTIONAL INPUT KEYWORDS
+      R_V - scalar specifying the ratio of total to selective extinction
+               R(V) = A(V) / E(B - V).  If not specified, then R = 3.1
+               Extreme values of R(V) range from 2.3 to 5.3
+
+      avglmc - if set, then the default fit parameters c1,c2,c3,c4,gamma,x0
+             are set to the average values determined for reddening in the
+             general Large Magellanic Cloud (LMC) field by Misselt et al.
+            (1999, ApJ, 515, 128)
+      lmc2 - if set, then the fit parameters are set to the values determined
+             for the LMC2 field (including 30 Dor) by Misselt et al.
+             Note that neither /AVGLMC or /LMC2 will alter the default value
+             of R_V which is poorly known for the LMC.
+
+      The following five input keyword parameters allow the user to customize
+      the adopted extinction curve.  For example, see Clayton et al. (2003,
+      ApJ, 588, 871) for examples of these parameters in different interstellar
+      environments.
+
+      x0 - Centroid of 2200 A bump in microns (default = 4.596)
+      gamma - Width of 2200 A bump in microns (default = 0.99)
+      c3 - Strength of the 2200 A bump (default = 3.23)
+      c4 - FUV curvature (default = 0.41)
+      c2 - Slope of the linear UV extinction component
+           (default = -0.824 + 4.717 / R)
+      c1 - Intercept of the linear UV extinction component
+           (default = 2.030 - 3.007 * c2)
+
+    NOTES:
+       (1) The following comparisons between the FM curve and that of Cardelli,
+           Clayton, & Mathis (1989), (see ccm_unred.pro):
+           (a) - In the UV, the FM and CCM curves are similar for R < 4.0, but
+                 diverge for larger R
+           (b) - In the optical region, the FM more closely matches the
+                 monochromatic extinction, especially near the R band.
+       (2)  Many sightlines with peculiar ultraviolet interstellar extinction
+               can be represented with the FM curve, if the proper value of
+               R(V) is supplied.
+    REQUIRED MODULES:
+       scipy, numpy
+    REVISION HISTORY:
+       Written   W. Landsman        Raytheon  STX   October, 1998
+       Based on FMRCurve by E. Fitzpatrick (Villanova)
+       Added /LMC2 and /AVGLMC keywords,  W. Landsman   August 2000
+       Added ExtCurve keyword, J. Wm. Parker   August 2000
+       Assume since V5.4 use COMPLEMENT to WHERE  W. Landsman April 2006
+       Ported to Python, C. Theissen August 2012
+    """
+
+    # copied by J. Guy from E. Schlafly fm_unred function
+
+    from scipy.interpolate import CubicSpline
+
+    x = 10000. / np.array(wave)  # Convert to inverse microns
+    curve = np.zeros(x.shape)
+
+    if lmc2 :
+        if x0 == None: x0 = 4.626
+        if gamma == None: gamma =  1.05
+        if c4 == None: c4 = 0.42
+        if c3 == None: c3 = 1.92
+        if c2 == None: c2 = 1.31
+        if c1 == None: c1 = -2.16
+    elif avglmc :
+        if x0 == None: x0 = 4.596
+        if gamma == None: gamma = 0.91
+        if c4 == None: c4 = 0.64
+        if c3 == None: c3 =  2.73
+        if c2 == None: c2 = 1.11
+        if c1 == None: c1 = -1.28
+    else:
+        if x0 == None: x0 = 4.596
+        if gamma == None: gamma = 0.99
+        if c4 == None: c4 = 0.41
+        if c3 == None: c3 =  3.23
+        if c2 == None: c2 = -0.824 + 4.717 / R_V
+        if c1 == None: c1 = 2.030 - 3.007 * c2
+
+    # Compute UV portion of A(lambda)/E(B-V) curve using FM fitting function and
+    # R-dependent coefficients
+
+    xcutuv = 10000.0 / 2700.0
+    xspluv = 10000.0 / np.array([2700.0, 2600.0])
+
+    iuv = x >= xcutuv
+    iuv_comp = ~iuv
+
+    if len(x[iuv]) > 0: xuv = np.concatenate( (xspluv, x[iuv]) )
+    else: xuv = xspluv.copy()
+
+    yuv = c1  + c2 * xuv
+    yuv = yuv + c3 * xuv**2 / ( ( xuv**2 - x0**2 )**2 + ( xuv * gamma )**2 )
+
+    filter1 = xuv.copy()
+    filter1[xuv <= 5.9] = 5.9
+
+    yuv = yuv + c4 * ( 0.5392 * ( filter1 - 5.9 )**2 + 0.05644 * ( filter1 - 5.9 )**3 )
+    yuv = yuv + R_V
+    yspluv = yuv[0:2].copy()                  # save spline points
+
+    if len(x[iuv]) > 0: curve[iuv] = yuv[2:len(yuv)]      # remove spline points
+
+    # Compute optical portion of A(lambda)/E(B-V) curve
+    # using cubic spline anchored in UV, optical, and IR
+
+    xsplopir = np.concatenate(([0], 10000.0 / np.array([26500.0, 12200.0, 6000.0, 5470.0, 4670.0, 4110.0])))
+    ysplir   = np.array([0.0, 0.26469, 0.82925]) * R_V / 3.1
+    ysplop   = [np.polyval(np.array([2.13572e-04, 1.00270, -4.22809e-01]), R_V ),
+                np.polyval(np.array([-7.35778e-05, 1.00216, -5.13540e-02]), R_V ),
+                np.polyval(np.array([-3.32598e-05, 1.00184, 7.00127e-01]), R_V ),
+                np.polyval(np.array([-4.45636e-05, 7.97809e-04, -5.46959e-03, 1.01707, 1.19456] ), R_V ) ]
+
+    ysplopir = np.concatenate( (ysplir, ysplop) )
+
+    if len(iuv_comp) > 0:
+        cubic = CubicSpline(np.concatenate( (xsplopir,xspluv) ),
+                            np.concatenate( (ysplopir,yspluv) ), bc_type='natural')
+        curve[iuv_comp] = cubic( x[iuv_comp] )
+
+    return curve/R_V
+
+# https://ui.adsabs.harvard.edu/abs/2011ApJ...737..103S
+ebv_sfd_scaling = 0.86 #
+
+def dust_transmission(wave,ebv_sfd) :
+    """
+    Return the dust transmission.
+    The routine does internally multiply ebv_sfd by dust.ebv_sfd_scaling.
+    The Fitzpatrick dust extinction law is used, assuming R_V=3.1.
+
+    Args:
+        wave : 1D array of vacuum wavelength [Angstroms]
+        ebv_sfd : E(B-V) as derived from the map of Schlegel, Finkbeiner and Davis (1998)
+
+    Returns:
+        1D array of dust transmission (between 0 and 1)
+    """
+    Rv = 3.1
+
+    extinction = ext_fitzpatrick(np.atleast_1d(wave),R_V=Rv)
+    if np.isscalar(wave) :
+        extinction=float(extinction[0])
+
+    return 10**(-Rv*extinction*ebv_sfd_scaling * ebv_sfd/2.5)
 
 # The SFDMap and _Hemisphere classes and the _bilinear_interpolate and ebv
 # functions below were copied on Nov/20/2016 from
@@ -473,6 +629,8 @@ def __repr__(self):
                         self.fnames['south'], self.scaling))
 
 
+
+
 def ebv(*args, **kwargs):
     """Convenience function, equivalent to ``SFDMap().ebv(*args)``.
     """
@@ -482,3 +640,84 @@ def ebv(*args, **kwargs):
                south=kwargs.get('south', "SFD_dust_4096_sgp.fits"),
                scaling=kwargs.get('scaling', 1.))
     return m.ebv(*args, **kwargs)
+
+def main():
+    import matplotlib.pyplot as plt
+
+    plt.figure("reddening")
+
+    wave=np.linspace(3000,11000,1000)
+    rv=3.1
+    plt.plot(wave,ebv_sfd_scaling*rv*ext_odonnell(wave, Rv=rv),label="{} x Odonnell (1994)".format(ebv_sfd_scaling))
+    plt.plot(wave,ebv_sfd_scaling*rv*ext_ccm(wave, Rv=rv),label="{} x CCM (1989)".format(ebv_sfd_scaling))
+    plt.plot(wave,ebv_sfd_scaling*rv*ext_fitzpatrick(wave, R_V=rv),label="{} x Fitzpatrick (1999)".format(ebv_sfd_scaling))
+
+
+    # load filters and compute extinctions here
+    try :
+        import speclite.filters
+
+
+        sed=(wave/6000.)**-2.2 # dummy star with a spectrum close to a 7000K star (for wavelength>4000A)
+
+        fluxunits = 1e-17 * u.erg / u.s / u.cm**2 / u.Angstrom
+        photsys_survey={"N":"BASS+MZLS","S":"DECALS"}
+        count=0
+
+
+
+        for photsys in ["N","S"] :
+            for band in ["G","R","Z"] :
+                count += 1
+                color="C{}".format(count%10)
+
+                if photsys=="N" :
+                    if band.upper() in ["G","R"] :
+                        filtername="BASS-{}".format(band.lower())
+                    else :
+                        filtername="MzLS-z"
+                elif photsys=="S" :
+                    filtername="decam2014-{}".format(band.lower())
+
+                R=extinction_total_to_selective_ratio(band, photsys)
+
+                filter_response=speclite.filters.load_filter(filtername)
+
+                rv  = 3.1
+                ebv = 0.1
+                mag_no_extinction   = filter_response.get_ab_magnitude(sed*fluxunits,wave)
+                mag_with_extinction = filter_response.get_ab_magnitude((sed*10**(-0.4*rv*ebv*ext_fitzpatrick(wave, R_V=rv)))*fluxunits,wave)
+
+                f=np.interp(wave,filter_response.wavelength,filter_response.response)
+                fwave = np.sum(sed*f*wave)/np.sum(sed*f)
+                Rbis=(mag_with_extinction-mag_no_extinction)/ebv
+
+                if count==1 :
+                    label1="with extinction_total_to_selective_ratio"
+                    label2="with {} x ext_fitzpatrick".format(ebv_sfd_scaling)
+                else :
+                    label1=None
+                    label2=None
+
+                plt.plot(fwave,R,"x",color=color,label=label1)
+                plt.plot(fwave,ebv_sfd_scaling*Rbis,"o",color=color,label=label2)
+
+                f=np.interp(wave,filter_response.wavelength,filter_response.response)
+                ii=(f>0.01)
+                plt.plot(wave[ii],f[ii],"--",color=color,label=filtername)
+                print("R_{}({}:{})= {:4.3f} =? {}x{:4.3f} = {:4.3f} , delta = {:5.4f}".format(band,photsys,photsys_survey[photsys],R,ebv_sfd_scaling,Rbis,ebv_sfd_scaling*Rbis,R-ebv_sfd_scaling*Rbis))
+                #print("R_ter=",-2.5*np.log10(dust_transmission(fwave,ebv_sfd=ebv))/ebv)
+
+
+    except ImportError as e :
+        print(e)
+        print("Cannot import speclite, so no broadband comparison")
+
+    plt.legend(title="For Rv={}".format(rv))
+    plt.xlabel("wavelength (A)")
+    plt.ylabel("A(wavelength)/E(B-V)")
+    plt.grid()
+    plt.show()
+
+if __name__ == '__main__':
+    main()