Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Summer + fall work #10

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Summer + fall work #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b84251e
Lots of changes
yonboyage Nov 8, 2023
515f0e2
Merge branch 'JulianBMunoz:main' into main
yonboyage Nov 8, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 11 additions & 1 deletion zeus21/constants.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
######################



precisionboost = 1.0 #makes integrals take more points for boost in precision. Baseline = 1.0
FLAG_FORCE_LINEAR_CF = 0 #0 to do standard calculation, 1 to force linearization of correlation function

Expand Down Expand Up @@ -54,7 +55,7 @@
EN_ION_HeI = 24.59 #eV
sigma0norm = 1e-18 #cm^2 normalization of xray cross sections

ZMAX_INTEGRAL = 35.0 #at which z we start the integrals. We take 35. as fiducial since there is not much SFRD. Make sure to test if you have exotic cosmology/astrophysics
ZMAX_INTEGRAL = 50.0 #35.0 #at which z we start the integrals. We take 35. as fiducial since there is not much SFRD. Make sure to test if you have exotic cosmology/astrophysics

#LyA related constants
wavelengthLyC = 91.1753 ##lyman continuum in nm
Expand All @@ -75,6 +76,15 @@
Tstar_21 = 0.0682 #T* in K for the 21cm transition
A10_21 = 2.85e-15 #1/s, Einstein 10 coeff for HI

############################################################################################

Tk_HH = [1, 2, 4, 6, 8, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 500, 700, 1000, 2000, 3000, 5000, 7000, 10000]
k_HH = [1.38e-13, 1.43e-13, 2.71e-13, 6.60e-13, 1.47e-12, 2.88e-12, 9.10e-12, 1.78e-11, 2.73e-11, 3.67e-11, 5.38e-11, 6.86e-11, 8.14e-11, 9.25e-11, 1.02e-10, 1.11e-10, 1.19e-10, 1.75e-10, 2.09e-10, 2.56e-10, 2.91e-10, 3.31e-10, 4.27e-10, 4.97e-10, 6.03e-10, 6.87e-10, 7.87e-10]
Tk_He = [1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000]
k_He = [2.39e-10, 3.37e-10, 5.30e-10, 7.46e-10, 1.05e-9, 1.63e-9, 2.26e-9, 3.11e-9, 4.59e-9, 5.92e-9, 7.15e-9, 7.71e-9, 8.17e-9, 8.32e-9, 8.37e-9, 8.29e-9, 8.11e-9]

############################################################################################


#whether to renormalize the C2 oefficients (appendix in 2302.08506)
C2_RENORMALIZATION_FLAG = 1 - FLAG_FORCE_LINEAR_CF
Expand Down
57 changes: 47 additions & 10 deletions zeus21/sfrd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,11 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
self.SFRD_avg = np.zeros_like(self.zintegral)
self.SFRDbar2D = np.zeros((self.Nzintegral, Cosmo_Parameters.NRs)) #SFR at z=zprime when averaged over a radius R (so up to a higher z)
self.gamma_index2D = np.zeros_like(self.SFRDbar2D) #index of SFR ~ exp(\gamma delta)
self.gamma2_index2D = np.zeros_like(self.SFRDbar2D) #quadratic fit gamma like ~ exp(\gamma \delta + \gamma_2 \delta^2)

self.niondot_avg = np.zeros_like(self.zintegral) #\dot nion at each z (int d(SFRD)/dM *fesc(M) dM)/rhobaryon
self.gamma_Niondot_index2D = np.zeros_like(self.SFRDbar2D) #index of SFR ~ exp(\gamma delta)
self.gamma_Niondot_index2D = np.zeros_like(self.SFRDbar2D) #index of SFR ~ exp(\gamma \delta)
self.gamma2_Niondot_index2D = np.zeros_like(self.SFRDbar2D) #index of SFR ~ exp(\gamma \delta + \gamma_2 \delta^2)

self.ztabRsmoo = np.zeros_like(self.SFRDbar2D) #z's that correspond to each Radius R around each zp

Expand All @@ -64,8 +67,8 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola


#and EPS factors
Nsigmad = 0.2 #how many sigmas we explore
Nds = 2 #how many deltas
Nsigmad = .2 #how many sigmas we explore (.2 originally)
Nds = 3 #how many deltas (2 originally)
deltatab_norm = np.linspace(-Nsigmad,Nsigmad,Nds)


Expand Down Expand Up @@ -197,6 +200,7 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
deltatab = deltatab_norm * sigmaR

SFRD_delta = np.zeros_like(deltatab)
Niondot_delta = np.zeros_like(deltatab)

nu0 = Cosmo_Parameters.delta_crit_ST/sigmacurr #the EPS delta = 0 result. Note 1/sigmacurr and not sigmamod
nu0[indextoobig]=1.0 #set to 1 to avoid under/overflows, we don't sum over those masses since they're too big
Expand All @@ -220,10 +224,38 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola

integrand_EPS *= (1.0 + dd) #to convert from Lagrangian to Eulerian
SFRD_delta[idelta] = np.trapz(integrand_EPS, HMF_interpolator.logtabMh)



self.gamma_index2D[izp,iR] = np.log(SFRD_delta[-1]/SFRD_delta[0])/(deltatab[-1] - deltatab[0]) #notice its der wrt delta, not delta/sigma or growth
Niondot_delta[idelta] = np.trapz(integrand_EPS * fesctab, HMF_interpolator.logtabMh) #ionizing emissivity (Niondot) is basically SFRD but with an escape fraction that varies with halo mass


midpoint = len(SFRD_delta)//2

self.gamma_index2D[izp,iR] = np.log(SFRD_delta[midpoint+1]/SFRD_delta[midpoint-1])/(deltatab[midpoint+1] - deltatab[midpoint-1]) #notice its der wrt delta, not delta/sigma or growth

self.gamma_Niondot_index2D[izp, iR] = np.log(Niondot_delta[midpoint+1]/Niondot_delta[midpoint-1])/(deltatab[midpoint+1] - deltatab[midpoint-1])


'''if (izp == 35)*(iR == 16): #z~10 and R~10Mpc
print(izp)
print(iR)
print(zp)
print(RR)
self.testSFRD = SFRD_delta
self.testdeltatab = deltatab
self.testNiondot = Niondot_delta
self.testfesctab = fesctab'''



der1 = np.log(SFRD_delta[midpoint]/SFRD_delta[midpoint-1])/(deltatab[midpoint] - deltatab[midpoint-1]) #ln(y2/y1)/(x2-x1)
der2 = np.log(SFRD_delta[midpoint+1]/SFRD_delta[midpoint])/(deltatab[midpoint+1] - deltatab[midpoint]) #ln(y3/y2)/(x3-x2)

self.gamma2_index2D[izp,iR] = (der2 - der1)/(deltatab[midpoint+1] - deltatab[midpoint-1]) #second derivative: (der2-der1)/((x3-x1)/2)

der1_N = np.log(Niondot_delta[midpoint]/Niondot_delta[midpoint-1])/(deltatab[midpoint] - deltatab[midpoint-1])
der2_N = np.log(Niondot_delta[midpoint+1]/Niondot_delta[midpoint])/(deltatab[midpoint+1] - deltatab[midpoint])

self.gamma2_Niondot_index2D[izp, iR] = (der2_N - der1_N)/(deltatab[midpoint+1] - deltatab[midpoint-1])



self.coeff2LyAzpRR[izp] = self.Rtabsmoo * self.dlogRR * self.SFRDbar2D[izp,:] * LyAintegral/ constants.yrTos/constants.Mpctocm**2
Expand Down Expand Up @@ -305,6 +337,11 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
self._coeff_Ja_xa_0 = 1.66e11/(1+self.zintegral) #They use a fixed (and slightly ~10% off) value.
else:
self._coeff_Ja_xa_0 = 8.0*np.pi*(constants.wavelengthLyA/1e7)**2 * constants.widthLyA * constants.Tstar_21/(9.0*constants.A10_21*self.T_CMB) #units of (cm^2 s Hz sr), convert from Ja to xa. should give 1.81e11/(1+self.zintegral) for Tcmb_0=2.725 K

#collision coefficient fh(1-x_e)
self.xc_HH = Cosmo_Parameters.f_H * (1.0 - self.xe_avg) * cosmology.n_baryon(Cosmo_Parameters, self.zintegral) * np.interp(self.Tk_avg, constants.Tk_HH, constants.k_HH) / constants.A10_21 * constants.Tstar_21 / cosmology.Tcmb(ClassCosmo, self.zintegral)
self.xc_He = self.xe_avg * cosmology.n_baryon(Cosmo_Parameters, self.zintegral) * np.interp(self.Tk_avg, constants.Tk_He, constants.k_He) / constants.A10_21 * constants.Tstar_21 / cosmology.Tcmb(ClassCosmo, self.zintegral) #xe
self.xc_avg = self.xc_HH + self.xc_He


if(constants.FLAG_WF_ITERATIVE==True): #iteratively find Tcolor and Ts. Could initialize one to zero, but this should converge faster
Expand All @@ -314,7 +351,7 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
self.invTcol_avg = 1.0 / self.Tk_avg
self.coeff_Ja_xa = self._coeff_Ja_xa_0 * Salpha_exp(self.zintegral, self.Tk_avg, self.xe_avg)
self.xa_avg = self.coeff_Ja_xa * self.Jalpha_avg
self._invTs_avg = (1.0/self.T_CMB+self.xa_avg*self.invTcol_avg)/(1+self.xa_avg)
self._invTs_avg = (1.0/self.T_CMB+self.xa_avg*self.invTcol_avg+self.xc_avg*self.invTk_avg)/(1+self.xa_avg+self.xc_avg)

_invTs_tryfirst = self._invTs_avg #so the loop below does not trigger

Expand All @@ -332,7 +369,7 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
self.invTcol_avg = 1.0/self.Tk_avg + constants.gcolorfactorHirata * 1.0/self.Tk_avg * (_invTs_tryfirst - 1.0/self.Tk_avg)

#and finally Ts^-1
self._invTs_avg = (1.0/self.T_CMB+self.xa_avg * self.invTcol_avg)/(1+self.xa_avg)
self._invTs_avg = (1.0/self.T_CMB+self.xa_avg * self.invTcol_avg + self.xc_avg * 1.0/self.Tk_avg)/(1+self.xa_avg+self.xc_avg)



Expand Down Expand Up @@ -376,7 +413,7 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola


#and finally, get the signal
self.T21avg = cosmology.T021(Cosmo_Parameters,self.zintegral) * self.xa_avg/(1.0 + self.xa_avg) * (1.0 - self.T_CMB * self.invTcol_avg) * self.xHI_avg
self.T21avg = cosmology.T021(Cosmo_Parameters,self.zintegral) * (self.xa_avg + self.xc_avg)/(1.0 + self.xa_avg + self.xc_avg) * (1.0 - self.T_CMB * self.invTcol_avg) * self.xHI_avg



Expand Down