Skip to content

Commit

Permalink
Merge 8b6af6d into 9dd77fe
Browse files Browse the repository at this point in the history
  • Loading branch information
@julienguy
julienguy committed Apr 17, 2021
2 parents 9dd77fe + 8b6af6d commit 6812ad5
Showing 1 changed file with 279 additions and 16 deletions.
295 changes: 279 additions & 16 deletions py/desiutil/dust.py
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@

# Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*- coding: utf-8 -*-
"""
Expand All @@ -14,10 +15,10 @@
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) :
def extinction_total_to_selective_ratio(band, photsys, match_legacy_surveys=False) :
"""Return the linear coefficient R_B = A(B)/E(B-V) where
A(B) = -2.5*log10(transmission in B band), for band B in 'G','R' or 'Z',
the optical bands of the legacy surveys. photsys = 'N' or 'S' specifies
Expand All @@ -30,23 +31,44 @@ def extinction_total_to_selective_ratio(band, photsys) :
Returns:
scalar, total extinction A(band) = -2.5*log10(transmission(band))
"""
if match_legacy_surveys :
# Based on the fit from the columns MW_TRANSMISSION_X and EBV
# for the DR8 target catalogs and propagated in fibermaps
# R_X = -2.5*log10(MW_TRANSMISSION_X) / EBV
# It is the same value for the N and S surveys in DR8 and DR9 catalogs.
R={"G_N":3.2140,
"R_N":2.1650,
"Z_N":1.2110,
"G_S":3.2140,
"R_S":2.1650,
"Z_S":1.2110,
}
else :
# From https://desi.lbl.gov/trac/wiki/ImagingStandardBandpass
# DECam u 3881.6 3.994
# DECam g 4830.8 3.212
# DECam r 6409.0 2.164
# DECam i 7787.5 1.591
# DECam z 9142.7 1.211
# DECam Y 9854.5 1.063
# BASS g 4772.1 3.258
# BASS r 6383.6 2.176
# MzLS z 9185.1 1.199
# Consistent with the synthetic magnitudes and function dust_transmission

R={"G_N":3.258,
"R_N":2.176,
"Z_N":1.199,
"G_S":3.212,
"R_S":2.164,
"Z_S":1.211,
}

# Based on the fit from the columns MW_TRANSMISSION_X and EBV
# 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.

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}
assert(band.upper() in ["G","R","Z"])
assert(photsys.upper() in ["N","S"])
return R["{}_{}".format(band.upper(),photsys.upper())]

def mwdust_transmission(ebv, band, photsys):
def mwdust_transmission(ebv, band, photsys, match_legacy_surveys=False):
"""Convert SFD E(B-V) value to dust transmission 0-1 for band and photsys
Args:
Expand All @@ -61,7 +83,7 @@ def mwdust_transmission(ebv, band, photsys):
However, `ebv` can be an array with a str `photsys`.
"""
if isinstance(photsys, str):
r_band = extinction_total_to_selective_ratio(band, photsys)
r_band = extinction_total_to_selective_ratio(band, photsys, match_legacy_surveys=match_legacy_surveys)
a_band = r_band * ebv
transmission = 10**(-a_band / 2.5)
return transmission
Expand All @@ -76,7 +98,7 @@ def mwdust_transmission(ebv, band, photsys):
transmission = np.zeros(len(ebv))
for p in np.unique(photsys):
ii = (photsys == p)
r_band = extinction_total_to_selective_ratio(band, p)
r_band = extinction_total_to_selective_ratio(band, p, match_legacy_surveys=match_legacy_surveys)
a_band = r_band * ebv[ii]
transmission[ii] = 10**(-a_band / 2.5)

Expand Down Expand Up @@ -181,6 +203,170 @@ 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
Expand Down Expand Up @@ -473,6 +659,8 @@ def __repr__(self):
self.fnames['south'], self.scaling))




def ebv(*args, **kwargs):
"""Convenience function, equivalent to ``SFDMap().ebv(*args)``.
"""
Expand All @@ -482,3 +670,78 @@ 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"
else :
label1=None

plt.plot(fwave,R,"x",color=color,label=label1)

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()

0 comments on commit 6812ad5

Please sign in to comment.