Skip to content

Commit

Permalink
Merge pull request #136 from desihub/igm-attenuation
Browse files Browse the repository at this point in the history 
update IGM attenuation model
  • Loading branch information
@moustakas
moustakas authored Jul 25, 2023
2 parents d746b76 + de33adf commit a33514c
Show file tree
Hide file tree
Showing 5 changed files with 293 additions and 36 deletions.
2 changes: 2 additions & 0 deletions doc/changes.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,14 @@ Change Log
* Bug fix of reversed tied flux ratio of [OII]7320,7330 doublet [`PR #120`_].
* Do not constrain the SPS age by default [`PR #132`_].
* Bug fix of emission-line subtracted Dn(4000) measurement [`PR #135`_].
* Update IGM attenuation coefficients [`PR #136`_].

.. _`PR #115`: https://github.com/desihub/fastspecfit/pull/115
.. _`PR #116`: https://github.com/desihub/fastspecfit/pull/116
.. _`PR #120`: https://github.com/desihub/fastspecfit/pull/120
.. _`PR #132`: https://github.com/desihub/fastspecfit/pull/132
.. _`PR #135`: https://github.com/desihub/fastspecfit/pull/135
.. _`PR #136`: https://github.com/desihub/fastspecfit/pull/136

2.1.2 (2023-04-01)
------------------
Expand Down
239 changes: 213 additions & 26 deletions py/fastspecfit/continuum.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -531,6 +531,8 @@ def __init__(self, nophoto=False, load_filters=True):
"""Class to load filters, dust, and filter- and dust-related methods.
"""
super(Filters, self).__init__()

from speclite import filters

self.bands = np.array(['g', 'r', 'z', 'W1', 'W2', 'W3', 'W4'])
Expand Down Expand Up @@ -811,7 +813,215 @@ def _integrate(wave, flux, ivar, w1, w2):

return dn4000, dn4000_ivar

class ContinuumTools(Filters):
class Inoue14(object):
def __init__(self, scale_tau=1.):
"""
IGM absorption from Inoue et al. (2014)
Parameters
----------
scale_tau : float
Parameter multiplied to the IGM :math:`\tau` values (exponential
in the linear absorption fraction).
I.e., :math:`f_\mathrm{igm} = e^{-\mathrm{scale\_tau} \tau}`.
Copyright (c) 2016-2022 Gabriel Brammer
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""
super(Inoue14, self).__init__()

self._load_data()
self.scale_tau = scale_tau

@staticmethod
def _pow(a, b):
"""C-like power, a**b
"""
return a**b

def _load_data(self):
from importlib import resources
LAF_file = resources.files('fastspecfit').joinpath('data/LAFcoeff.txt')
DLA_file = resources.files('fastspecfit').joinpath('data/DLAcoeff.txt')

data = np.loadtxt(LAF_file, unpack=True)
ix, lam, ALAF1, ALAF2, ALAF3 = data
self.lam = lam[:,np.newaxis]
self.ALAF1 = ALAF1[:,np.newaxis]
self.ALAF2 = ALAF2[:,np.newaxis]
self.ALAF3 = ALAF3[:,np.newaxis]

data = np.loadtxt(DLA_file, unpack=True)
ix, lam, ADLA1, ADLA2 = data
self.ADLA1 = ADLA1[:,np.newaxis]
self.ADLA2 = ADLA2[:,np.newaxis]

return True


@property
def NA(self):
"""
Number of Lyman-series lines
"""
return self.lam.shape[0]


def tLSLAF(self, zS, lobs):
"""
Lyman series, Lyman-alpha forest
"""
z1LAF = 1.2
z2LAF = 4.7

l2 = self.lam #[:, np.newaxis]
tLSLAF_value = np.zeros_like(lobs*l2).T

x0 = (lobs < l2*(1+zS))
x1 = x0 & (lobs < l2*(1+z1LAF))
x2 = x0 & ((lobs >= l2*(1+z1LAF)) & (lobs < l2*(1+z2LAF)))
x3 = x0 & (lobs >= l2*(1+z2LAF))

tLSLAF_value = np.zeros_like(lobs*l2)
tLSLAF_value[x1] += ((self.ALAF1/l2**1.2)*lobs**1.2)[x1]
tLSLAF_value[x2] += ((self.ALAF2/l2**3.7)*lobs**3.7)[x2]
tLSLAF_value[x3] += ((self.ALAF3/l2**5.5)*lobs**5.5)[x3]

return tLSLAF_value.sum(axis=0)


def tLSDLA(self, zS, lobs):
"""
Lyman Series, DLA
"""
z1DLA = 2.0

l2 = self.lam #[:, np.newaxis]
tLSDLA_value = np.zeros_like(lobs*l2)

x0 = (lobs < l2*(1+zS)) & (lobs < l2*(1.+z1DLA))
x1 = (lobs < l2*(1+zS)) & ~(lobs < l2*(1.+z1DLA))

tLSDLA_value[x0] += ((self.ADLA1/l2**2)*lobs**2)[x0]
tLSDLA_value[x1] += ((self.ADLA2/l2**3)*lobs**3)[x1]

return tLSDLA_value.sum(axis=0)


def tLCDLA(self, zS, lobs):
"""
Lyman continuum, DLA
"""
z1DLA = 2.0
lamL = 911.8

tLCDLA_value = np.zeros_like(lobs)

x0 = lobs < lamL*(1.+zS)
if zS < z1DLA:
tLCDLA_value[x0] = 0.2113 * self._pow(1.0+zS, 2) - 0.07661 * self._pow(1.0+zS, 2.3) * self._pow(lobs[x0]/lamL, (-3e-1)) - 0.1347 * self._pow(lobs[x0]/lamL, 2)
else:
x1 = lobs >= lamL*(1.+z1DLA)

tLCDLA_value[x0 & x1] = 0.04696 * self._pow(1.0+zS, 3) - 0.01779 * self._pow(1.0+zS, 3.3) * self._pow(lobs[x0 & x1]/lamL, (-3e-1)) - 0.02916 * self._pow(lobs[x0 & x1]/lamL, 3)
tLCDLA_value[x0 & ~x1] =0.6340 + 0.04696 * self._pow(1.0+zS, 3) - 0.01779 * self._pow(1.0+zS, 3.3) * self._pow(lobs[x0 & ~x1]/lamL, (-3e-1)) - 0.1347 * self._pow(lobs[x0 & ~x1]/lamL, 2) - 0.2905 * self._pow(lobs[x0 & ~x1]/lamL, (-3e-1))

return tLCDLA_value


def tLCLAF(self, zS, lobs):
"""
Lyman continuum, LAF
"""
z1LAF = 1.2
z2LAF = 4.7
lamL = 911.8

tLCLAF_value = np.zeros_like(lobs)

x0 = lobs < lamL*(1.+zS)

if zS < z1LAF:
tLCLAF_value[x0] = 0.3248 * (self._pow(lobs[x0]/lamL, 1.2) - self._pow(1.0+zS, -9e-1) * self._pow(lobs[x0]/lamL, 2.1))
elif zS < z2LAF:
x1 = lobs >= lamL*(1+z1LAF)
tLCLAF_value[x0 & x1] = 2.545e-2 * (self._pow(1.0+zS, 1.6) * self._pow(lobs[x0 & x1]/lamL, 2.1) - self._pow(lobs[x0 & x1]/lamL, 3.7))
tLCLAF_value[x0 & ~x1] = 2.545e-2 * self._pow(1.0+zS, 1.6) * self._pow(lobs[x0 & ~x1]/lamL, 2.1) + 0.3248 * self._pow(lobs[x0 & ~x1]/lamL, 1.2) - 0.2496 * self._pow(lobs[x0 & ~x1]/lamL, 2.1)
else:
x1 = lobs > lamL*(1.+z2LAF)
x2 = (lobs >= lamL*(1.+z1LAF)) & (lobs < lamL*(1.+z2LAF))
x3 = lobs < lamL*(1.+z1LAF)

tLCLAF_value[x0 & x1] = 5.221e-4 * (self._pow(1.0+zS, 3.4) * self._pow(lobs[x0 & x1]/lamL, 2.1) - self._pow(lobs[x0 & x1]/lamL, 5.5))
tLCLAF_value[x0 & x2] = 5.221e-4 * self._pow(1.0+zS, 3.4) * self._pow(lobs[x0 & x2]/lamL, 2.1) + 0.2182 * self._pow(lobs[x0 & x2]/lamL, 2.1) - 2.545e-2 * self._pow(lobs[x0 & x2]/lamL, 3.7)
tLCLAF_value[x0 & x3] = 5.221e-4 * self._pow(1.0+zS, 3.4) * self._pow(lobs[x0 & x3]/lamL, 2.1) + 0.3248 * self._pow(lobs[x0 & x3]/lamL, 1.2) - 3.140e-2 * self._pow(lobs[x0 & x3]/lamL, 2.1)

return tLCLAF_value


def full_IGM(self, z, lobs):
"""Get full Inoue IGM absorption
Parameters
----------
z : float
Redshift to evaluate IGM absorption
lobs : array
Observed-frame wavelength(s) in Angstroms.
Returns
-------
abs : array
IGM absorption
"""
tau_LS = self.tLSLAF(z, lobs) + self.tLSDLA(z, lobs)
tau_LC = self.tLCLAF(z, lobs) + self.tLCDLA(z, lobs)

### Upturn at short wavelengths, low-z
#k = 1./100
#l0 = 600-6/k
#clip = lobs/(1+z) < 600.
#tau_clip = 100*(1-1./(1+np.exp(-k*(lobs/(1+z)-l0))))
tau_clip = 0.

return np.exp(-self.scale_tau*(tau_LC + tau_LS + tau_clip))


def build_grid(self, zgrid, lrest):
"""Build a spline interpolation object for fast IGM models
Returns: self.interpolate
"""

from scipy.interpolate import CubicSpline
igm_grid = np.zeros((len(zgrid), len(lrest)))
for iz in range(len(zgrid)):
igm_grid[iz,:] = self.full_IGM(zgrid[iz], lrest*(1+zgrid[iz]))

self.interpolate = CubicSpline(zgrid, igm_grid)


class ContinuumTools(Filters, Inoue14):
"""Tools for dealing with stellar continua.
Parameters
Expand All @@ -835,7 +1045,7 @@ class ContinuumTools(Filters):
def __init__(self, nophoto=False, continuum_pixkms=25.0, pixkms_wavesplit=1e4):

super(ContinuumTools, self).__init__(nophoto=nophoto)

from fastspecfit.emlines import read_emlines

self.massnorm = 1e10 # stellar mass normalization factor [Msun]
Expand Down Expand Up @@ -1075,29 +1285,6 @@ def build_linemask(wave, flux, ivar, redshift=0.0, nsig=7.0, linetable=None,

return linemask_dict

@staticmethod
def transmission_Lyman(zObj, lObs):
"""Calculate the transmitted flux fraction from the Lyman series
This returns the transmitted flux fraction:
1 -> everything is transmitted (medium is transparent)
0 -> nothing is transmitted (medium is opaque)
Args:
zObj (float): Redshift of object
lObs (array of float): wavelength grid
Returns:
array of float: transmitted flux fraction
"""
from fastspecfit.util import Lyman_series

lRF = lObs/(1.+zObj)
T = np.ones(lObs.size)
for l in list(Lyman_series.keys()):
w = lRF<Lyman_series[l]['line']
zpix = lObs[w]/Lyman_series[l]['line']-1.
tauEff = Lyman_series[l]['A']*(1.+zpix)**Lyman_series[l]['B']
T[w] *= np.exp(-tauEff)
return T

@staticmethod
def smooth_and_resample(templateflux, templatewave, specwave=None, specres=None):
Expand Down Expand Up @@ -1216,7 +1403,7 @@ def templates2data(self, _templateflux, _templatewave, redshift=0.0, dluminosity
# stellar mass.
if redshift > 0:
ztemplatewave = templatewave * (1.0 + redshift)
T = self.transmission_Lyman(redshift, ztemplatewave)
T = self.full_IGM(redshift, ztemplatewave)
T *= FLUXNORM * self.massnorm * (10.0 / (1e6 * dluminosity))**2 / (1.0 + redshift)
ztemplateflux = templateflux * T[:, np.newaxis]
else:
Expand Down
39 changes: 39 additions & 0 deletions py/fastspecfit/data/DLAcoeff.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
2 1215.670 1.61698E-04 5.38995E-05
3 1025.720 1.54539E-04 5.15129E-05
4 972.537 1.49767E-04 4.99222E-05
5 949.743 1.46031E-04 4.86769E-05
6 937.803 1.42893E-04 4.76312E-05
7 930.748 1.40159E-04 4.67196E-05
8 926.226 1.37714E-04 4.59048E-05
9 923.150 1.35495E-04 4.51650E-05
10 920.963 1.33452E-04 4.44841E-05
11 919.352 1.31561E-04 4.38536E-05
12 918.129 1.29785E-04 4.32617E-05
13 917.181 1.28117E-04 4.27056E-05
14 916.429 1.26540E-04 4.21799E-05
15 915.824 1.25041E-04 4.16804E-05
16 915.329 1.23614E-04 4.12046E-05
17 914.919 1.22248E-04 4.07494E-05
18 914.576 1.20938E-04 4.03127E-05
19 914.286 1.19681E-04 3.98938E-05
20 914.039 1.18469E-04 3.94896E-05
21 913.826 1.17298E-04 3.90995E-05
22 913.641 1.16167E-04 3.87225E-05
23 913.480 1.15071E-04 3.83572E-05
24 913.339 1.14011E-04 3.80037E-05
25 913.215 1.12983E-04 3.76609E-05
26 913.104 1.11972E-04 3.73241E-05
27 913.006 1.11002E-04 3.70005E-05
28 912.918 1.10051E-04 3.66836E-05
29 912.839 1.09125E-04 3.63749E-05
30 912.768 1.08220E-04 3.60734E-05
31 912.703 1.07337E-04 3.57789E-05
32 912.645 1.06473E-04 3.54909E-05
33 912.592 1.05629E-04 3.52096E-05
34 912.543 1.04802E-04 3.49340E-05
35 912.499 1.03991E-04 3.46636E-05
36 912.458 1.03198E-04 3.43994E-05
37 912.420 1.02420E-04 3.41402E-05
38 912.385 1.01657E-04 3.38856E-05
39 912.353 1.00908E-04 3.36359E-05
40 912.324 1.00168E-04 3.33895E-05
39 changes: 39 additions & 0 deletions py/fastspecfit/data/LAFcoeff.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
2 1215.670 1.68976E-02 2.35379E-03 1.02611E-04
3 1025.720 4.69229E-03 6.53625E-04 2.84940E-05
4 972.537 2.23898E-03 3.11884E-04 1.35962E-05
5 949.743 1.31901E-03 1.83735E-04 8.00974E-06
6 937.803 8.70656E-04 1.21280E-04 5.28707E-06
7 930.748 6.17843E-04 8.60640E-05 3.75186E-06
8 926.226 4.60924E-04 6.42055E-05 2.79897E-06
9 923.150 3.56887E-04 4.97135E-05 2.16720E-06
10 920.963 2.84278E-04 3.95992E-05 1.72628E-06
11 919.352 2.31771E-04 3.22851E-05 1.40743E-06
12 918.129 1.92348E-04 2.67936E-05 1.16804E-06
13 917.181 1.62155E-04 2.25878E-05 9.84689E-07
14 916.429 1.38498E-04 1.92925E-05 8.41033E-07
15 915.824 1.19611E-04 1.66615E-05 7.26340E-07
16 915.329 1.04314E-04 1.45306E-05 6.33446E-07
17 914.919 9.17397E-05 1.27791E-05 5.57091E-07
18 914.576 8.12784E-05 1.13219E-05 4.93564E-07
19 914.286 7.25069E-05 1.01000E-05 4.40299E-07
20 914.039 6.50549E-05 9.06198E-06 3.95047E-07
21 913.826 5.86816E-05 8.17421E-06 3.56345E-07
22 913.641 5.31918E-05 7.40949E-06 3.23008E-07
23 913.480 4.84261E-05 6.74563E-06 2.94068E-07
24 913.339 4.42740E-05 6.16726E-06 2.68854E-07
25 913.215 4.06311E-05 5.65981E-06 2.46733E-07
26 913.104 3.73821E-05 5.20723E-06 2.27003E-07
27 913.006 3.45377E-05 4.81102E-06 2.09731E-07
28 912.918 3.19891E-05 4.45601E-06 1.94255E-07
29 912.839 2.97110E-05 4.13867E-06 1.80421E-07
30 912.768 2.76635E-05 3.85346E-06 1.67987E-07
31 912.703 2.58178E-05 3.59636E-06 1.56779E-07
32 912.645 2.41479E-05 3.36374E-06 1.46638E-07
33 912.592 2.26347E-05 3.15296E-06 1.37450E-07
34 912.543 2.12567E-05 2.96100E-06 1.29081E-07
35 912.499 1.99967E-05 2.78549E-06 1.21430E-07
36 912.458 1.88476E-05 2.62543E-06 1.14452E-07
37 912.420 1.77928E-05 2.47850E-06 1.08047E-07
38 912.385 1.68222E-05 2.34330E-06 1.02153E-07
39 912.353 1.59286E-05 2.21882E-06 9.67268E-08
40 912.324 1.50996E-05 2.10334E-06 9.16925E-08
10 changes: 0 additions & 10 deletions py/fastspecfit/util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,16 +17,6 @@
except:
C_LIGHT = 299792.458 # [km/s]

# Lyman-alpha from eqn 5 of Calura et al. 2012 (Arxiv: 1201.5121)
# Other from eqn 1.1 of Irsic et al. 2013 , (Arxiv: 1307.3403)
Lyman_series = {
'Lya' : { 'line':1215.67, 'A':0.0023, 'B':3.64, 'var_evol':3.8 },
'Lyb' : { 'line':1025.72, 'A':0.0023/5.2615, 'B':3.64, 'var_evol':3.8 },
'Ly3' : { 'line':972.537, 'A':0.0023/14.356, 'B':3.64, 'var_evol':3.8 },
'Ly4' : { 'line':949.7431, 'A':0.0023/29.85984, 'B':3.64, 'var_evol':3.8 },
'Ly5' : { 'line':937.8035, 'A':0.0023/53.36202, 'B':3.64, 'var_evol':3.8 },
}

def ivar2var(ivar, clip=1e-3, sigma=False, allmasked_ok=False):
"""Safely convert an inverse variance to a variance. Note that we clip at 1e-3
by default, not zero.
Expand Down

0 comments on commit a33514c

Please sign in to comment.