Commit be791d60 authored by Guang's avatar Guang

Updating IGM:

1. remove "weight"
2. allow for short wavelength (<912 A)
3. remove ad hoc "X-ray wavelength"
parent fbaa7b1b
......@@ -22,12 +22,11 @@ from collections import OrderedDict
import numpy as np
from scipy.constants import parsec
from scipy.misc import factorial
from scipy.special import factorial
from astropy.cosmology import WMAP7 as cosmology
from . import SedModule
def igm_transmission(wavelength, redshift):
"""Intergalactic transmission (Meiksin, 2006)
......@@ -51,7 +50,6 @@ def igm_transmission(wavelength, redshift):
gamma = 0.2788 # Gamma(0.5,1) i.e., Gamma(2-beta,1) with beta = 1.5
n0 = 0.25
lambda_limit = 91.2 # Lyman limit in nm
lam_xray_max = 10**0.5 # Maximum X-ray wavelength in nm
lambda_n = np.empty(n_transitions_max)
z_n = np.empty((n_transitions_max, len(wavelength)))
......@@ -92,11 +90,12 @@ def igm_transmission(wavelength, redshift):
(float(n) * (float(n*n - 1.))))
for n in range(2, n_transitions_max):
w = np.where(z_n[n, :] >= redshift)
# If z_n>=redshift or z_n<0, the photon cannot be absorbed by Lyman n->1
w = np.where( (z_n[n, :] >= redshift) | (z_n[n, :] < 0) )
tau_n[n, w] = 0.
z_l = wavelength / lambda_limit - 1.
w = np.where(z_l < redshift)
w = np.where( z_l < redshift )
tau_l_igm = np.zeros_like(wavelength)
tau_l_igm[w] = (0.805 * np.power(1. + z_l[w], 3) *
......@@ -120,23 +119,18 @@ def igm_transmission(wavelength, redshift):
tau_l_lls = np.zeros_like(wavelength)
tau_l_lls[w] = n0 * ((term1 - term2) * term3 - term4)
tau_taun = np.sum(tau_n[2:n_transitions_max, :], axis=0)
lambda_min_igm = (1+redshift)*70.
w = np.where(wavelength < lambda_min_igm)
# Reset for short wavelength (z_l<0)
w = np.where( z_l<0 )
# Get the normalization factor at z_l=0
tau_norm_l_igm = np.interp(0, z_l, tau_l_igm)
tau_norm_l_lls = np.interp(0, z_l, tau_l_lls)
# Calculate tau_l_igm & tau_l_lls, assuming cross section ~ lambda^2.75 (O'Meara et al. 2013)
tau_l_igm[w] = tau_norm_l_igm * (z_l[w]+1)**2.75
tau_l_lls[w] = tau_norm_l_lls * (z_l[w]+1)**2.75
weight = np.ones_like(wavelength)
weight[w] = np.power(wavelength[w]/lambda_min_igm, 2.)
# Another weight using erf function can be used.
# However, you would need to add: from scipy.special import erf
# weight[w] = 0.5*(1.+erf(0.05*(wavelength[w]-lambda_min_igm)))
igm_transmission = np.exp(-tau_taun-tau_l_igm-tau_l_lls) * weight
tau_taun = np.sum(tau_n[2:n_transitions_max, :], axis=0)
# Reset the transmission for X-ray (E>0.4 keV) to 1
igm_transmission[np.where( \
(wavelength/(1+redshift)) <= lam_xray_max \
)] = 1.
igm_transmission = np.exp(-tau_taun-tau_l_igm-tau_l_lls)
return igm_transmission
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment