Merge a0bb1fa 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) :
@@ -35,13 +35,15 @@ def extinction_total_to_selective_ratio(band, photsys) :
     # for the DR8 target catalogs and propagated in fibermaps
     # R_X = -2.5*log10(MW_TRANSMISSION_X) / EBV
     # excerpt from https://www.legacysurvey.org/dr8/catalogs/#galactic-extinction-coefficients :  Eddie Schlafly has computed the extinction coefficients for the DECam filters through airmass=1.3, computed for a 7000K source spectrum as was done in the Appendix of Schlafly & Finkbeiner (2011). These coefficients are A / E(B-V) = 3.995, 3.214, 2.165, 1.592, 1.211, 1.064 (note that these are slightly different from the coefficients in Schlafly & Finkbeiner 2011). The coefficients are multiplied by the SFD98 E(B-V) values at the coordinates of each object to derive the g, r and z mw_transmission values in the Legacy Surveys catalogs. The coefficients at different airmasses only change by a small amount, with the largest effect in g -band where the coefficient would be 3.219 at airmass=1 and 3.202 at airmass=2. We calculate Galactic extinction for BASS and MzLS as if they are on the DECam filter system.
-
+    # It is the same value for the N and S surveys in DR8 and DR9 catalogs.
+    # So the values for the N survey (BASS+MzLS) are not very accurate but we do not change them here for consistency with the Legacy Survey data.
     R={"G_N":3.2140,
        "R_N":2.1650,
        "Z_N":1.2110,
-       "G_S":3.2829,
-       "R_S":2.1999,
-       "Z_S":1.2150}
+       "G_S":3.2140,
+       "R_S":2.1650,
+       "Z_S":1.2110,
+       }
     assert(band.upper() in ["G","R","Z"])
     assert(photsys.upper() in ["N","S"])
     return R["{}_{}".format(band.upper(),photsys.upper())]
@@ -181,6 +183,165 @@ 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
+
+# based on the work from https://ui.adsabs.harvard.edu/abs/2011ApJ...737..103S
+# note from Eddie: I recommend applying the SF11 calibration in the following way:
+# A(lambda, F99) = A(lambda, F99, rv)/A(1 micron, F99, rv) * SFD_EBV * 1.029.
+# That's a definition that only uses monochromatic extinctions,
+# so a lot of ambiguity in what extinction means goes away.
+# That doesn't look like the 0.86 rescaling that we quote in abstract,
+# but it's convenient because it uses only monochromatic extinctions.
+
+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) / ext_fitzpatrick(np.array([10000.]),R_V=Rv) * ebv_sfd * 1.029
+    if np.isscalar(wave) :
+        extinction=float(extinction[0])
+    return 10**(-extinction/2.5)
 
 # The SFDMap and _Hemisphere classes and the _bilinear_interpolate and ebv
 # functions below were copied on Nov/20/2016 from
@@ -473,6 +634,8 @@ def __repr__(self):
                         self.fnames['south'], self.scaling))
 
 
+
+
 def ebv(*args, **kwargs):
     """Convenience function, equivalent to ``SFDMap().ebv(*args)``.
     """
@@ -482,3 +645,81 @@ 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
+    ebv_sfd = 1.
+    trans = dust_transmission(wave,ebv_sfd)
+    ext   = -2.5*np.log10(trans)
+    plt.plot(wave,ext,label="-2.5*log10(dust_transmission)")
+
+    # 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*dust_transmission(wave,ebv)*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"
+                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} =? {:4.3f}, delta = {:5.4f}".format(band,photsys,photsys_survey[photsys],R,Rbis,R-Rbis))
+
+
+    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()

py/desiutil/test/test_dust.py

@@ -132,10 +132,10 @@ def test_total_to_selective(self):
                 self.assertEqual(rb1, rb2)
 
         #- North and South should be different
-        for band in ['G', 'R', 'Z']:
-            rb1 = dust.extinction_total_to_selective_ratio(band, 'N')
-            rb2 = dust.extinction_total_to_selective_ratio(band, 'S')
-            self.assertNotEqual(rb1, rb2)
+        #for band in ['G', 'R', 'Z']:
+        #    rb1 = dust.extinction_total_to_selective_ratio(band, 'N')
+        #    rb2 = dust.extinction_total_to_selective_ratio(band, 'S')
+        #    self.assertNotEqual(rb1, rb2)
 
         #- B is not a supported band (G,R,Z)
         with self.assertRaises(AssertionError):
@@ -164,8 +164,8 @@ def test_mwdust_transmission(self):
         self.assertEqual(len(tn), len(ebv))
         self.assertEqual(len(ts), len(ebv))
         #- N vs. S should be different where ebv>0
-        ii = (ebv>0)
-        self.assertTrue(np.all(tn[ii] != ts[ii]))
+        #ii = (ebv>0)
+        #self.assertTrue(np.all(tn[ii] != ts[ii]))
 
         #- array photsys must have ebv array of same length
         with self.assertRaises(ValueError):