Skip to content

Commit

Permalink
Merge pull request #1 from yasada96/update_igm
Browse files Browse the repository at this point in the history
Update igm
  • Loading branch information
yasada96 authored Oct 23, 2024
2 parents e524d3f + 64d7d75 commit 32e5210
Show file tree
Hide file tree
Showing 5 changed files with 336 additions and 19 deletions.
Binary file added .DS_Store
Binary file not shown.
7 changes: 6 additions & 1 deletion eazy/data/zphot.param.default
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ SMOOTH_SIGMA 100. # Gaussian sigma (in Angstroms) to smoot

## Templates
TEMPLATES_FILE templates/fsps_full/tweak_fsps_QSF_12_v3.param # Template definition file
TEMPLATE_COMBOS a # Template combination options:
TEMPLATE_COMBOS a # Template combination options:
# 1 : one template at a time
# 2 : two templates, read allowed combinations from TEMPLATES_FILE
# -2 : two templates, all permutations
Expand All @@ -22,6 +22,11 @@ TEMP_ERR_A2 0.20 # Template error amplitude
SYS_ERR 0.01 # Systematic flux error (% of flux)
APPLY_IGM y # Apply Madau 1995 IGM absorption
IGM_SCALE_TAU 1.0 # Scale factor times Inoue14 IGM tau
ADD_CGM y # Add Asada24 CGM damping wing absorption
SIGMOID_PARAM1 3.4835 # Sigmoid func parameter (A) for the N_HI-z relation in Asada24
SIGMOID_PARAM2 1.2581 # Sigmoid func parameter (a) for the N_HI-z relation in Asada24
SIGMOID_PARAM3 18.249 # Sigmoid func parameter (C) for the N_HI-z relation in Asada24


SCALE_2175_BUMP 0.00 # Scaling of 2175A bump. Values 0.13 (0.27) absorb ~10 (20) % at peak.
TEMPLATE_SMOOTH 0.0 # Velocity smoothing (km/s) for templates, < 0 for no smoothing
Expand Down
275 changes: 272 additions & 3 deletions eazy/igm.py
Original file line number Diff line number Diff line change
@@ -1,9 +1,243 @@
import os
import numpy as np
#from math import pow as _pow

from . import __file__ as filepath

__all__ = ["Inoue14"]
__all__ = ["Asada24","Inoue14"]

class Asada24(object):
def __init__(self, sigmoid_params=(3.48347968, 1.25809685, 18.24922789), scale_tau=1., add_cgm=True):
"""
Compute IGM+CGM transmission from Asada et al. 2024, in prep.
The IGM model is from Inoue+ (2014).
Parameters
----------
sigmoid_params : 3-tuple of float
Parameters that controll the redshift evolution of the CGM HI gas column density.
The defaul values are from Asada et al. (2024).
scale_tau : float
Scalar multiplied to tau_igm
add_cgm : bool
Add the additional LyA damping absorption at z>6 as described in Asada+24.
If False, the transmission will be identical to Inoue+ 2014
"""
self._load_data()
self.sigmoid_params = sigmoid_params
self.scale_tau = scale_tau
self.add_cgm = add_cgm


def tau_cgm(self, N_HI, lam, z):
"""
CGM optical depth given by Totani+06, eqn (1)
Parameters
----------
N_HI : float
HI column density [cm-2]
lam : 1D array
wavelength array in the observed frame [AA]
z : float
Redshift of the source
Returns
-------
1D array of tau_CGM
"""

lam_rest = lam/(1+z)
nu_rest = 3e18/lam_rest

tau = np.zeros(len(lam))

for i, nu in enumerate(nu_rest):
tau[i] = N_HI * sigma_a(nu)

return tau

def lgNHI_z(self, z):
"""
HI column density as a function of redshift, calibrated in Asada+ 2024.
Only valid at z>=6
Parameters
----------
z : float
Redshift of the source
Returns
-------
log10(HI column density [cm-2])
"""

lgN_HI = sigmoid(z, self.sigmoid_params[0], self.sigmoid_params[1], self.sigmoid_params[2])

return lgN_HI

### IGM part -- identical to Inoue+14
def _load_data(self):
path = os.path.join(os.path.dirname(filepath),'data')
#print path

LAF_file = os.path.join(path, 'LAFcoeff.txt')
DLA_file = os.path.join(path, '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 * _pow(1.0+zS, 2) - 0.07661 * _pow(1.0+zS, 2.3) * _pow(lobs[x0]/lamL, (-3e-1)) - 0.1347 * _pow(lobs[x0]/lamL, 2)
else:
x1 = lobs >= lamL*(1.+z1DLA)

tLCDLA_value[x0 & x1] = 0.04696 * _pow(1.0+zS, 3) - 0.01779 * _pow(1.0+zS, 3.3) * _pow(lobs[x0 & x1]/lamL, (-3e-1)) - 0.02916 * _pow(lobs[x0 & x1]/lamL, 3)
tLCDLA_value[x0 & ~x1] =0.6340 + 0.04696 * _pow(1.0+zS, 3) - 0.01779 * _pow(1.0+zS, 3.3) * _pow(lobs[x0 & ~x1]/lamL, (-3e-1)) - 0.1347 * _pow(lobs[x0 & ~x1]/lamL, 2) - 0.2905 * _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 * (_pow(lobs[x0]/lamL, 1.2) - _pow(1.0+zS, -9e-1) * _pow(lobs[x0]/lamL, 2.1))
elif zS < z2LAF:
x1 = lobs >= lamL*(1+z1LAF)
tLCLAF_value[x0 & x1] = 2.545e-2 * (_pow(1.0+zS, 1.6) * _pow(lobs[x0 & x1]/lamL, 2.1) - _pow(lobs[x0 & x1]/lamL, 3.7))
tLCLAF_value[x0 & ~x1] = 2.545e-2 * _pow(1.0+zS, 1.6) * _pow(lobs[x0 & ~x1]/lamL, 2.1) + 0.3248 * _pow(lobs[x0 & ~x1]/lamL, 1.2) - 0.2496 * _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 * (_pow(1.0+zS, 3.4) * _pow(lobs[x0 & x1]/lamL, 2.1) - _pow(lobs[x0 & x1]/lamL, 5.5))
tLCLAF_value[x0 & x2] = 5.221e-4 * _pow(1.0+zS, 3.4) * _pow(lobs[x0 & x2]/lamL, 2.1) + 0.2182 * _pow(lobs[x0 & x2]/lamL, 2.1) - 2.545e-2 * _pow(lobs[x0 & x2]/lamL, 3.7)
tLCLAF_value[x0 & x3] = 5.221e-4 * _pow(1.0+zS, 3.4) * _pow(lobs[x0 & x3]/lamL, 2.1) + 0.3248 * _pow(lobs[x0 & x3]/lamL, 1.2) - 3.140e-2 * _pow(lobs[x0 & x3]/lamL, 2.1)

return tLCLAF_value

def full_IGM(self, z, lobs):
"""Get full IGM+CGM absorption
Parameters
----------
z : float
Redshift to evaluate IGM absorption
lobs : array-like
Observed-frame wavelength(s) in Angstroms.
Returns
-------
abs : array-like
IGM+CGM transmission factor
"""

tau_LS = self.tLSLAF(z, lobs) + self.tLSDLA(z, lobs)
tau_LC = self.tLCLAF(z, lobs) + self.tLCDLA(z, lobs)
tau_clip = 0.

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

if self.add_cgm:
if (z<6):
tau_C = np.zeros(len(lobs))
else:
NHI = 10**(self.lgNHI_z(z))
tau_C = self.tau_cgm(NHI, lobs, z)
cgmz = np.exp(-tau_C)
else:
cgmz = np.ones(len(lobs))

return igmz * cgmz


class Inoue14(object):
def __init__(self, scale_tau=1.):
Expand All @@ -13,8 +247,8 @@ def __init__(self, scale_tau=1.):
Parameters
----------
scale_tau : float
Parameter multiplied to the IGM :math:`\tau` values (exponential
in the linear absorption fraction).
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}`.
"""
self._load_data()
Expand Down Expand Up @@ -192,3 +426,38 @@ def _pow(a, b):
"""
return a**b

def sigmoid(x,A,a,c):
"""
Sigmoid function centered at x=6
"""

return A/(1+np.exp(-a*(x-6))) + c


def sigma_a(nu_rest):
"""
Lyα absorption cross section for the restframe frequency \nu given by Totani+06, eqn (1)
for CGM damping wing calculation
Parameters
----------
nu : float
rest-frame frequency [Hz]
Returns
-------
Lyα absorption cross section at the restframe frequency \nu_rest [cm2]
"""

Lam_a = 6.255486e8 ## Hz
nu_lya = 2.46607e15 ## Hz

C = 6.9029528e22 ## Constant factor [Ang2/s2]

sig = C *(nu_rest/nu_lya)**4 / (4*np.pi**2 * (nu_rest - nu_lya)**2 + Lam_a**2*(nu_rest/nu_lya)**6/4)

s = sig * 1e-16 ## convert AA-2 to cm-2

return s


17 changes: 13 additions & 4 deletions eazy/photoz.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@

from . import utils

IGM_OBJECT = igm_module.Inoue14()
IGM_OBJECT = igm_module.Asada24()

__all__ = ["PhotoZ", "TemplateGrid", "template_lsq", "fit_by_redshift"]

Expand Down Expand Up @@ -272,6 +272,15 @@ def __init__(self, param_file=None, translate_file=None, zeropoint_file=None, lo

if 'IGM_SCALE_TAU' in self.param.params:
IGM_OBJECT.scale_tau = self.param['IGM_SCALE_TAU']

if self.param['ADD_CGM'] in utils.TRUE_VALUES:
IGM_OBJECT.add_cgm = True
else:
IGM_OBJECT.add_cgm = False

sigmoid_params = (self.param.params['SIGMOID_PARAM1'], self.param.params['SIGMOID_PARAM2'], self.param.params['SIGMOID_PARAM3'])
IGM_OBJECT.sigmoid_params = sigmoid_params


### Read templates
kws = dict(templates_file=self.param['TEMPLATES_FILE'],
Expand Down Expand Up @@ -2673,7 +2682,7 @@ def show_fit(self, id, id_is_idx=False, zshow=None, show_fnu=0, get_spec=False,

if self.tempfilt.add_igm:
igmz = templ.wave*0.+1
lyman = templ.wave < 1300
lyman = templ.wave < 2000
igmz[lyman] = IGM_OBJECT.full_IGM(z, templz[lyman])
else:
igmz = 1.
Expand Down Expand Up @@ -2832,7 +2841,7 @@ def show_fit(self, id, id_is_idx=False, zshow=None, show_fnu=0, get_spec=False,
templzi = templ.wave*(1+zi)
if self.tempfilt.add_igm:
igmz = templ.wave*0.+1
lyman = templ.wave < 1300
lyman = templ.wave < 2000
igmz[lyman] = IGM_OBJECT.full_IGM(zi, templzi[lyman])
else:
igmz = 1.
Expand Down Expand Up @@ -5688,7 +5697,7 @@ def _integrate_tempfilt(itemp, templ, zgrid, RES, f_numbers, add_igm, galactic_e
# IGM absorption
if add_igm:
igmz = templ.wave*0.+1
lyman = templ.wave < 1300
lyman = templ.wave < 2000
igmz[lyman] = igm.full_IGM(zgrid[iz], lz[lyman])
else:
igmz = 1.
Expand Down
Loading

0 comments on commit 32e5210

Please sign in to comment.