diff --git a/doc/changes.rst b/doc/changes.rst index b75817df..be1a039d 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -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) ------------------ diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index bae9e80b..63efe732 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -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']) @@ -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 @@ -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] @@ -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 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: diff --git a/py/fastspecfit/data/DLAcoeff.txt b/py/fastspecfit/data/DLAcoeff.txt new file mode 100644 index 00000000..ba9ff068 --- /dev/null +++ b/py/fastspecfit/data/DLAcoeff.txt @@ -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 diff --git a/py/fastspecfit/data/LAFcoeff.txt b/py/fastspecfit/data/LAFcoeff.txt new file mode 100644 index 00000000..3023deee --- /dev/null +++ b/py/fastspecfit/data/LAFcoeff.txt @@ -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 diff --git a/py/fastspecfit/util.py b/py/fastspecfit/util.py index 2a5eccb8..ee08a6f3 100644 --- a/py/fastspecfit/util.py +++ b/py/fastspecfit/util.py @@ -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.