Skip to content

Commit

Permalink
Merge pull request #41 from mtakahiro/igm_attenuation
Browse files Browse the repository at this point in the history 
Igm attenuation
  • Loading branch information
@mtakahiro
mtakahiro authored Dec 9, 2022
2 parents d6792ae + dbabffd commit 126413f
Show file tree
Hide file tree
Showing 9 changed files with 380 additions and 140 deletions.
56 changes: 37 additions & 19 deletions gsf/fitting.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -637,8 +637,15 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set:
except:
self.force_agefix = False

print('\n')
# SFH prior;
# try:
# self.norder_sfh_prior = int(inputs['norder_sfh_prior'])
# self.f_prior_sfh = True
# except:
# self.norder_sfh_prior = None
# self.f_prior_sfh = False

print('\n')

def get_lines(self, LW0):
fLW = np.zeros(len(LW0), dtype='int')
Expand Down Expand Up @@ -886,7 +893,7 @@ def residual_z(pars,z):

def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
zlimu=None, snlim=0, priors=None, f_bb_zfit=True, f_line_check=False,
f_norm=True, f_lambda=False, zmax=20, fit_photometry=False, f_exclude_negative=False):
f_norm=True, f_lambda=False, zmax=20, include_bb=False, f_exclude_negative=False):
'''
Purpose
-------
Expand All @@ -909,8 +916,8 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
Minimum redshift uncertainty.
snlim : float
SN limit for data points. Those below the number will be cut from the fit.
f_bb_zfit : bool
Redshift fitting if only BB data. If False, returns nothing.
include_bb : bool
Turn this True if Redshift fitting is requred for only BB data.
f_line_check : bool
If True, line masking.
priors : dict, optional
Expand Down Expand Up @@ -942,17 +949,11 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
sn = self.dict['fy'] / self.dict['ey']

# Only spec data?
if fit_photometry:
if include_bb:
con_cz = ()#(sn>snlim)
else:
con_cz = (self.dict['NR']<NRbb_lim) #& (sn>snlim)

if len(self.dict['fy'][con_cz])==0:
if f_bb_zfit:
con_cz = ()#(sn>snlim)
else:
return 'y'

if f_exclude_negative:
con_cz &= (sn>snlim)

Expand Down Expand Up @@ -1004,7 +1005,10 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
self.p_prior = prior_s

# Plot;
if self.fzvis==1:
if len(self.dict['fy'][con_cz])==0:
return 'y'

if self.fzvis==1 and (include_bb | len(fy_cz)>0):
import matplotlib as mpl
mpl.use('TkAgg')
data_model = np.zeros((len(x_cz),4),'float')
Expand All @@ -1027,8 +1031,13 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
print('############################')
print('Start MCMC for redshift fit')
print('############################')
res_cz, fitc_cz = check_redshift(fy_cz, ey_cz, x_cz, fm_tmp, xm_tmp/(1+self.zgal), self.zgal, self.z_prior, self.p_prior, \
NR_cz, zliml, zlimu, self.nmc_cz, self.nwalk_cz, include_photometry=False)
res_cz, fitc_cz = check_redshift(
fy_cz, ey_cz, x_cz, fm_tmp, xm_tmp/(1+self.zgal),
self.zgal, self.z_prior, self.p_prior,
NR_cz, zliml, zlimu, self.nmc_cz, self.nwalk_cz,
include_photometry=include_photometry_zfit
)

z_cz = np.percentile(res_cz.flatchain['z'], [16,50,84])
scl_cz0 = np.percentile(res_cz.flatchain['Cz0'], [16,50,84])
scl_cz1 = np.percentile(res_cz.flatchain['Cz1'], [16,50,84])
Expand Down Expand Up @@ -1361,24 +1370,30 @@ def set_param(self):
self.Aini = 1

if self.SFH_FORM==-99:
self.age_vary = []
if len(self.age) != len(self.aamin):
for aa in range(len(self.age)):
if aa not in self.aamin:
fit_params.add('A'+str(aa), value=self.Amin, vary=False)
self.ndim -= 1
self.age_vary.append(False)
else:
fit_params.add('A'+str(aa), value=self.Aini, min=self.Amin, max=self.Amax)
self.age_vary.append(True)
else:
for aa in range(len(self.age)):
if self.age[aa] == 99:
fit_params.add('A'+str(aa), value=self.Amin, vary=False)
self.ndim -= 1
self.age_vary.append(False)
elif self.age[aa]>agemax and not self.force_agefix:
print('At this redshift, A%d is beyond the age of universe and not being used.'%(aa))
fit_params.add('A'+str(aa), value=self.Amin, vary=False)
self.ndim -= 1
self.age_vary.append(False)
else:
fit_params.add('A'+str(aa), value=self.Aini, min=self.Amin, max=self.Amax)
self.age_vary.append(True)

else:
for aa in range(self.npeak):
Expand Down Expand Up @@ -1581,7 +1596,7 @@ def get_shuffle(self, out, nshuf=3.0, amp=1e-4):
def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0,
f_move:bool=False, verbose:bool=False, skip_fitz:bool=False, out=None, f_plot_accept:bool=True,
f_shuffle:bool=True, amp_shuffle=1e-2, check_converge:bool=True, Zini=None, f_plot_chain:bool=True,
f_chind:bool=True, ncpu:int=0, f_prior_sfh:bool=False):
f_chind:bool=True, ncpu:int=0, f_prior_sfh:bool=False, norder_sfh_prior:int=3):
'''
Main module of this script.
Expand Down Expand Up @@ -1671,7 +1686,7 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0,
if skip_fitz:
flag_z = 'y'
else:
flag_z = self.fit_redshift(xm_tmp, fm_tmp, delzz=0.001, fit_photometry=False)
flag_z = self.fit_redshift(xm_tmp, fm_tmp, delzz=0.001, include_bb=False)

#################################################
# Gor for mcmc phase
Expand Down Expand Up @@ -1752,7 +1767,8 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0,
sampler = zeus.EnsembleSampler(self.nwalk, self.ndim, class_post.lnprob_emcee, \
args=[out.params, self.dict['fy'], self.dict['ey'], self.dict['wht2'], self.f_dust], \
moves=moves, maxiter=1e6,\
kwargs={'f_val':True, 'out':out, 'lnpreject':-np.inf, 'f_chind':f_chind, 'f_prior_sfh':f_prior_sfh},\
kwargs={'f_val':True, 'out':out, 'lnpreject':-np.inf, 'f_chind':f_chind,
'f_prior_sfh':f_prior_sfh, 'norder_sfh_prior':norder_sfh_prior},\
)
# Run MCMC
nburn = int(self.nmc/10)
Expand All @@ -1772,7 +1788,8 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0,
sampler = zeus.EnsembleSampler(self.nwalk, self.ndim, class_post.lnprob_emcee, \
args=[out.params, self.dict['fy'], self.dict['ey'], self.dict['wht2'], self.f_dust], \
moves=moves, maxiter=1e4,\
kwargs={'f_val':True, 'out':out, 'lnpreject':-np.inf, 'f_prior_sfh':f_prior_sfh},\
kwargs={'f_val':True, 'out':out, 'lnpreject':-np.inf,
'f_prior_sfh':f_prior_sfh, 'norder_sfh_prior':norder_sfh_prior},\
)

else:
Expand All @@ -1781,7 +1798,8 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0,
sampler = emcee.EnsembleSampler(self.nwalk, self.ndim, class_post.lnprob_emcee, \
args=(out.params, self.dict['fy'], self.dict['ey'], self.dict['wht2'], self.f_dust),\
#moves=moves,\
kwargs={'f_val': True, 'out': out, 'lnpreject':-np.inf, 'f_prior_sfh':f_prior_sfh},\
kwargs={'f_val': True, 'out': out, 'lnpreject':-np.inf,
'f_prior_sfh':f_prior_sfh, 'norder_sfh_prior':norder_sfh_prior},\
)

if check_converge:
Expand Down
9 changes: 6 additions & 3 deletions gsf/function_class.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@

from .function import *
from .basic_func import Basic
from .function_igm import dijkstra_igm_abs


class Func:
Expand Down Expand Up @@ -173,7 +174,7 @@ def open_spec_neb_fits(self, fall:int = 0, orig:bool = False):
return lib


def open_spec_fits_dir(self, nage:int, nz:int, kk, Av00:float, zgal:float, A00:float):
def open_spec_fits_dir(self, nage:int, nz:int, kk, Av00:float, zgal:float, A00:float, f_IGM=True):
'''
Load template in obs range.
But for weird template.
Expand Down Expand Up @@ -430,7 +431,8 @@ def tmp03_neb(self, A00, Av00, logU, nmodel, Z, zgal, lib, f_apply_dust=True, EB
return A00 * yyd_sort, xxd_sort


def tmp04(self, par, f_Alog:bool=True, nprec:int=1, pp:int = 0, f_val:bool=False, lib_all:bool=False, f_nrd:bool=False):
def tmp04(self, par, f_Alog:bool=True, nprec:int=1, pp:int = 0, f_val:bool=False, lib_all:bool=False, f_nrd:bool=False,
f_IGM=True):
'''
Makes model template with a given param set.
Also dust attenuation.
Expand Down Expand Up @@ -509,6 +511,7 @@ def tmp04(self, par, f_Alog:bool=True, nprec:int=1, pp:int = 0, f_val:bool=False

self.MB.logMtmp = np.log10(Mtot)

# @@@ Filter convolution may need to happpen here
if round(zmc,nprec) != round(self.MB.zgal,nprec):
fint = interpolate.interp1d(xx, yy, kind='nearest', fill_value="extrapolate")
xx_s = xx / (1+self.MB.zgal) * (1+zmc)
Expand All @@ -532,8 +535,8 @@ def tmp04(self, par, f_Alog:bool=True, nprec:int=1, pp:int = 0, f_val:bool=False
yyd, xxd, nrd = dust_kc(xx/(1.+zmc), yy, Av00, nr, Rv=4.05, gamma=-0.2)
else:
yyd, xxd, nrd = dust_calz(xx/(1.+zmc), yy, Av00, nr)
xxd *= (1.+zmc)

xxd *= (1.+zmc)
if self.dust_model != 0:
# This may be needed when not calzetti model
nrd_yyd = np.zeros((len(nrd),3), dtype=float)
Expand Down
152 changes: 149 additions & 3 deletions gsf/function_igm.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,153 @@
import numpy as np
import sys
from scipy.integrate import simps
from scipy import integrate


def get_XI(z, zout=6, zin=8):
'''
Very simplified model.
'''
zs = np.linspace(zout, zin, 100)
if z < zout:
XI = 0
elif z > zin:
XI = 1.
else:
XI = (z-zout) / (zin-zout) * 1.0

return XI


def get_dtdz(z, zs, dtdzs):
'''
'''
iix = np.argmin(np.abs(zs[:-1]-z))
return dtdzs[iix]


def masongronke_igm_abs(xtmp, ytmp, zin, cosmo=None, xLL=1216., c=3e18, ckms=3e5, zobs=6, xLLL=1400):
'''
Purpose
-------
Parameters
----------
xtmp : float array
Rest-frame wavelength, in AA.
ytmp : float array
flux, in f_lambda.
zin : target redshift of IGM application
Returns
-------
IGM attenuated flux.
'''
if cosmo == None:
from astropy.cosmology import WMAP9 as cosmo

tau = np.zeros(len(xtmp), dtype=float)
xobs = xtmp * (1.*zin)
ytmp_abs = np.zeros(len(ytmp), float)
zs = np.linspace(zobs,zin,1000)

xtmp_obs = xtmp * (1+zin)
x_HI = 0.8
T = 1e0 #1e4 # K
sigma_0 = 5.9e-14 * (T / 1e4)**(-1/2) # cm2
a_V = 4.7e-4 * (T / 1e4)**(-1/2) #
nu_a = 2.46e15 #Hz
k_B = 1.380649e-23 / (1e3)**2 #m2 kg s-2 K-1 * (km/m)**2 = km2 kg/s2/K
m_p = 1.67262192e-27 #kilograms
delta_nu_d = nu_a * np.sqrt(2 * k_B * T / m_p / ckms**2)

dtdzs = (cosmo.age(zs)[0:-1].value - cosmo.age(zs)[1:].value) * 1e9 * 365.25 * 3600 * 24 / np.diff(zs) # s

# xLL = 1390
for ii in range(len(xtmp_obs)):
if xtmp[ii] < xLL:
tau[ii] = 100
elif xtmp[ii] < xLLL:
nu = c / xtmp_obs[ii] # Hz
x = (nu - nu_a) / delta_nu_d
phi_x = get_H(x,a_V) #

# tau[ii] = integrate.quad(lambda z: ckms * get_dtdz(z, zs, dtdzs) * x_HI * sigma_0 * phi_x, zobs, zin)[0] * get_column(zin, cosmo)
tau[ii] = integrate.quad(lambda z: ckms * get_dtdz(z, zs, dtdzs) * x_HI * (1.88e-7 * (1+z)**3) * sigma_0 * phi_x, zobs, zin)[0]
print(tau[ii], xtmp_obs[ii])
else:
tau[ii] = 1e-9

# R_b1 = 0.0 # Mpc
# x_D = 0.8 # neutral fraction
#NH = get_column(zin, cosmo)
ytmp_abs = ytmp * np.exp(-tau)
# con = ()
# ytmp_abs[con] = ytmp[con] * np.exp(-tau[con])

return ytmp_abs


def dijkstra_igm_abs(xtmp, ytmp, zin, cosmo=None, xLL=1216., ckms=3e5,
R_b1=1.0, delta_v_0=600, alpha_x=1.0):
'''
Purpose
-------
Apply IMG-attenuation of Dijikstra (2014).
https://www.cambridge.org/core/services/aop-cambridge-core/content/view/S1323358014000332
Parameters
----------
xtmp : float array
Rest-frame wavelength, in AA.
ytmp : float array
flux, in f_lambda.
zin :
target redshift of IGM application
R_b1 : float
Bubble size, in Mpc
Returns
-------
IGM attenuated flux.
'''
import scipy.interpolate as interpolate
if cosmo == None:
from astropy.cosmology import WMAP9 as cosmo

tau = np.zeros(len(xtmp), dtype=float)
xobs = xtmp * (1.*zin)
ytmp_abs = np.zeros(len(ytmp), float)

xtmp_obs = xtmp * (1+zin)

x_HI = get_XI(zin) # neutral fraction
x_D = alpha_x * x_HI # x_D is not clear..
delta_lam = (xtmp - xLL) * (zin + 1)
delta_lam_fine = (np.linspace(900,2000,1000) - xLL) * (zin + 1)

delta_v = ckms * delta_lam_fine / (xLL * (1.+zin))
delta_v_b1 = delta_v
if R_b1>0:
delta_v_b1 += cosmo.H(zin).value * R_b1 / (1.+zin) # km / (Mpc s) * Mpc

tau_fine = 2.3 * x_D * (delta_v_b1/delta_v_0)**(-1) * ((1+zin)/10)**(3/2)
con_tau = (tau_fine < 0) | (delta_v_b1 == 0)
tau_fine[con_tau] = 100

fint = interpolate.interp1d(delta_lam_fine, tau_fine, kind='nearest', fill_value="extrapolate")
tau = fint(delta_lam)

# import matplotlib.pyplot as plt
# plt.close()
# plt.plot(xtmp, tau[:] )
# plt.xlim(1200,1500)
# plt.show()
# hoge

ytmp_abs = ytmp * np.exp(-tau)

return ytmp_abs


def madau_igm_abs(xtmp, ytmp, zin, cosmo=None, xLL=1216.):
Expand All @@ -25,10 +172,9 @@ def madau_igm_abs(xtmp, ytmp, zin, cosmo=None, xLL=1216.):
if cosmo == None:
from astropy.cosmology import WMAP9 as cosmo

tau = np.zeros(len(xtmp), dtype='float')
#xLL = 912.
tau = np.zeros(len(xtmp), dtype=float)
xobs = xtmp * (1.*zin)
ytmp_abs = np.zeros(len(ytmp), 'float')
ytmp_abs = np.zeros(len(ytmp), float)

NH = get_column(zin, cosmo)
tau = (NH/1.6e17) * (xtmp/xLL)**(3.)
Expand Down
Loading

0 comments on commit 126413f

Please sign in to comment.