From b84251e64b03719dede287c3d686b3512cc5eac4 Mon Sep 17 00:00:00 2001 From: yonboyage <59982772+yonboyage@users.noreply.github.com> Date: Wed, 8 Nov 2023 17:36:51 -0600 Subject: [PATCH] Lots of changes --- zeus21/UVLFs.py | 121 ++++++++++++++++++++++++++++++++++++++++++++ zeus21/__init__.py | 1 + zeus21/constants.py | 12 ++++- zeus21/cosmology.py | 45 ++++++++++++++++ zeus21/inputs.py | 24 +++++++-- zeus21/sfrd.py | 57 +++++++++++++++++---- 6 files changed, 244 insertions(+), 16 deletions(-) create mode 100644 zeus21/UVLFs.py diff --git a/zeus21/UVLFs.py b/zeus21/UVLFs.py new file mode 100644 index 0000000..03d6601 --- /dev/null +++ b/zeus21/UVLFs.py @@ -0,0 +1,121 @@ +""" + +Compute UVLFs given our SFR and HMF models. + +Author: Julian B. Muñoz +UT Austin - June 2023 + +""" + +from . import cosmology +from . import constants +from .sfrd import SFR +from .cosmology import bias_Tinker + +import numpy as np +from scipy.special import erf + + + + + + +def MUV_of_SFR(SFRtab, kappaUV): + 'returns MUV, uses SFR. Dust added later in loglike.' + #convert SFR to MUVs + LUVtab = SFRtab/kappaUV + MUVtab = 51.63 - 2.5 * np.log10(LUVtab) #AB magnitude + return MUVtab + + +#and combine to get UVLF: +def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwidth, MUVcenters, MUVwidths, DUST_FLAG=True, RETURNBIAS = False): + 'Binned UVLF in units of 1/Mpc^3/mag, for bins at with a Gaussian width zwidth, centered at MUV centers with tophat width MUVwidths. z width only in HMF since that varies the most rapidly. If flag RETURNBIAS set to true it returns number-avgd bias instead of UVLF, still have to divide by UVLF' + + if(constants.NZ_TOINT>1): + DZ_TOINT = np.linspace(-np.sqrt(constants.NZ_TOINT/3.),np.sqrt(constants.NZ_TOINT/3.),constants.NZ_TOINT) #in sigmas around zcenter + else: + DZ_TOINT = np.array([0.0]) + WEIGHTS_TOINT = np.exp(-DZ_TOINT**2/2.)/np.sum(np.exp(-DZ_TOINT**2/2.)) #assumed Gaussian in z, fair + + + + + SFRlist = SFR(Astro_Parameters,Cosmo_Parameters,HMF_interpolator,zcenter) + sigmaUV = Astro_Parameters.sigmaUV + + if (constants.FLAG_RENORMALIZE_LUV == True): #lower the LUV (or SFR) to recover the true avg, not log-avg + SFRlist/= np.exp((np.log(10)/2.5*sigmaUV)**2/2.0) + + MUVbarlist = MUV_of_SFR(SFRlist, Astro_Parameters._kappaUV) #avg for each Mh + MUVbarlist = np.fmin(MUVbarlist,constants._MAGMAX) + + + + if(RETURNBIAS==True): # weight by bias + biasM = np.array([bias_Tinker(Cosmo_Parameters, HMF_interpolator.sigma_int(HMF_interpolator.Mhtab,zcenter+dz*zwidth)) for dz in DZ_TOINT]) + else: # do not weight by bias + biasM = np.ones_like(WEIGHTS_TOINT) + + + HMFtab = np.array([HMF_interpolator.HMF_int(HMF_interpolator.Mhtab,zcenter+dz*zwidth) for dz in DZ_TOINT]) + HMFcurr = np.sum(WEIGHTS_TOINT * HMFtab.T * biasM.T,axis=1) + + #cannot directly 'dust' the theory since the properties of the IRX-beta relation are calibrated on observed MUV. Recursion instead: + currMUV = MUVbarlist + if(DUST_FLAG==True): + currMUV2 = np.ones_like(currMUV) + while(np.sum(np.abs((currMUV2-currMUV)/currMUV)) > 0.02): + currMUV = MUVbarlist + AUV(Astro_Parameters,zcenter,currMUV) + currMUV2 = currMUV + + + MUVcuthi = MUVcenters + MUVwidths/2. + MUVcutlo = MUVcenters - MUVwidths/2. + + xhi = np.subtract.outer(MUVcuthi , currMUV)/(np.sqrt(2) * sigmaUV) + xlo = np.subtract.outer(MUVcutlo, currMUV )/(np.sqrt(2) * sigmaUV) + weights = (erf(xhi) - erf(xlo)).T/(2.0 * MUVwidths) + + UVLF_filtered = np.trapz(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1) + + return UVLF_filtered + + + + + +#####Here the dust attenuation +def AUV(Astro_Parameters, z, MUV, HIGH_Z_DUST = True, _zmaxdata=8.0): + 'Average attenuation A as a function of OBSERVED z and magnitude. If using on theory iterate until convergence. HIGH_Z_DUST is whether to do dust at higher z than 0 or set to 0. Fix at \beta(z=8) result if so' + + betacurr = beta(z,MUV) + + C0, C1 = Astro_Parameters.C0dust, Astro_Parameters.C1dust + + sigmabeta = 0.34 #from Bouwens 2014 + + Auv = C0 + 0.2*np.log(10)*sigmabeta**2 * C1**2 + C1 * betacurr + Auv=Auv.T + if not (HIGH_Z_DUST): + Auv*=np.heaviside(_zmaxdata - z,0.5) + Auv=Auv.T + return np.fmax(Auv, 0.0) + +def beta(z, MUV): + 'Color as a function of redshift and mag, interpolated from Bouwens 2013-14 data.' + + zdatbeta = [2.5,3.8,5.0,5.9,7.0,8.0] + betaMUVatM0 = [-1.7,-1.85,-1.91,-2.00,-2.05,-2.13] + dbeta_dMUV = [-0.20,-0.11,-0.14,-0.20,-0.20,-0.15] + + _MUV0 = -19.5 + _c = -2.33 + + betaM0 = np.interp(z, zdatbeta, betaMUVatM0, left=betaMUVatM0[0], right=betaMUVatM0[-1]) + dbetaM0 = (MUV - _MUV0).T * np.interp(z, zdatbeta, dbeta_dMUV, left=dbeta_dMUV[0], right=dbeta_dMUV[-1]) + + sol1 = (betaM0-_c) * np.exp(dbetaM0/(betaM0-_c))+_c #for MUV > MUV0 + sol2 = dbetaM0 + betaM0 #for MUV < MUV0 + + return sol1.T * np.heaviside(MUV - _MUV0, 0.5) + sol2.T * np.heaviside(_MUV0 - MUV, 0.5) diff --git a/zeus21/__init__.py b/zeus21/__init__.py index e1b81fd..957357d 100644 --- a/zeus21/__init__.py +++ b/zeus21/__init__.py @@ -4,6 +4,7 @@ from .correlations import * from .sfrd import get_T21_coefficients from .xrays import Xray_class +from .UVLFs import UVLF_binned import warnings warnings.filterwarnings("ignore", category=UserWarning) #to silence unnecessary warning in mcfit diff --git a/zeus21/constants.py b/zeus21/constants.py index 986564a..ac80aa4 100644 --- a/zeus21/constants.py +++ b/zeus21/constants.py @@ -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 @@ -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 @@ -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 diff --git a/zeus21/cosmology.py b/zeus21/cosmology.py index 75dcae5..94feb2f 100644 --- a/zeus21/cosmology.py +++ b/zeus21/cosmology.py @@ -12,7 +12,24 @@ from scipy.interpolate import RegularGridInterpolator from . import constants +from .inputs import Cosmo_Parameters, Cosmo_Parameters_Input +from .correlations import Correlations +def cosmo_wrapper(Cosmo_Parameters_Input): + """ + Wrapper function for all the cosmology. It takes Cosmo_Parameters_Input and returns: + Cosmo_Parameters, Class_Cosmo, Correlations, HMF_interpolator + """ + + ClassCosmo = Class() + ClassCosmo.compute() + + ClassyCosmo = runclass(Cosmo_Parameters_Input) + CosmoParams = Cosmo_Parameters(Cosmo_Parameters_Input, ClassyCosmo) + CorrFClass = Correlations(CosmoParams, ClassyCosmo) + HMFintclass = HMF_interpolator(CosmoParams,ClassyCosmo) + + return CosmoParams, ClassyCosmo, CorrFClass, HMFintclass @@ -240,6 +257,34 @@ def T021(Cosmo_Parameters, z): return 34 * pow((1+z)/16.,0.5) * (Cosmo_Parameters.omegab/0.022) * pow(Cosmo_Parameters.omegam/0.14,-0.5) +#UNUSED bias, just for reference +def bias_ST(Cosmo_Parameters, sigmaM): + # from https://arxiv.org/pdf/1007.4201.pdf Table 1 + a_ST = Cosmo_Parameters.a_ST + p_ST = Cosmo_Parameters.p_ST + delta_crit_ST = Cosmo_Parameters.delta_crit_ST + nu = delta_crit_ST/sigmaM + nutilde = np.sqrt(a_ST) * nu + + return 1.0 + (nutilde**2 - 1.0 + 2. * p_ST/(1.0 + nutilde**(2. * p_ST) ) )/delta_crit_ST + +def bias_Tinker(Cosmo_Parameters, sigmaM): + #from https://arxiv.org/pdf/1001.3162.pdf, Delta=200 + delta_crit_ST = Cosmo_Parameters.delta_crit_ST + nu = delta_crit_ST/sigmaM + + #Tinker fit + _Deltahalo = 200; + _yhalo = np.log10(_Deltahalo) + _Abias = 1.0 + 0.24 * _yhalo * np.exp(-(4.0/_yhalo)**4.) + _abias = 0.44*_yhalo-0.88 + _Bbias = 0.183 + _bbias = 1.5 + _Cbias = 0.019 + 0.107 * _yhalo + 0.19 * np.exp(-(4.0/_yhalo)**4.) + _cbias = 2.4 + + return 1.0 - _Abias*(nu**_abias/(nu**_abias + delta_crit_ST**_abias)) + _Bbias * nu**_bbias + _Cbias * nu**_cbias + #UNUSED: # def interp2Dlinear_only_y(arrayxy, arrayz, x, y): # "2D interpolator where the x axis is assumed to be an array identical to the trained x. That is, an array of 1D linear interpolators. arrayxy is [x,y]. arrayz is result. x is the x input (=arrayxy[0]), and y the y input. Returns z result (array)" diff --git a/zeus21/inputs.py b/zeus21/inputs.py index 5a76998..6899134 100644 --- a/zeus21/inputs.py +++ b/zeus21/inputs.py @@ -17,7 +17,7 @@ class Cosmo_Parameters_Input: "Class to pass the 6 LCDM parameters as input" - def __init__(self, omegab= 0.0223828, omegac = 0.1201075, h_fid = 0.67810, As = 2.100549e-09, ns = 0.9660499, tau_fid = 0.05430842, kmax_CLASS = 500., zmax_CLASS = 50., Flag_emulate_21cmfast = False): + def __init__(self, omegab= 0.0223828, omegac = 0.1201075, h_fid = 0.67810, As = 2.100549e-09, ns = 0.9660499, tau_fid = 0.05430842, kmax_CLASS = 500., zmax_CLASS = 50.,zmin_CLASS = 5., Flag_emulate_21cmfast = False): self.omegab = omegab self.omegac = omegac @@ -29,6 +29,7 @@ def __init__(self, omegab= 0.0223828, omegac = 0.1201075, h_fid = 0.67810, As = #other params for CLASS self.kmax_CLASS = kmax_CLASS self.zmax_CLASS = zmax_CLASS + self.zmin_CLASS = zmin_CLASS #and whether to emulate 21cmFAST self.Flag_emulate_21cmfast = Flag_emulate_21cmfast #whether to emulate 21cmFAST in HMF, LyA, and X-ray opacity calculations @@ -49,7 +50,7 @@ def __init__(self, CosmoParams_input, ClassCosmo): #other params in the input self.kmax_CLASS = CosmoParams_input.kmax_CLASS self.zmax_CLASS = CosmoParams_input.zmax_CLASS - self.zmin_CLASS = 4.0 #when to start the HMF calcs., not an input strictly + self.zmin_CLASS = CosmoParams_input.zmin_CLASS #when to start the HMF calcs., not an input strictly self.Flag_emulate_21cmfast = CosmoParams_input.Flag_emulate_21cmfast #whether to emulate 21cmFAST in HMF, LyA, and X-ray opacity calculations #derived params @@ -76,8 +77,12 @@ def __init__(self, CosmoParams_input, ClassCosmo): self._ztabinchi = np.linspace(0.0, 100. , 10000) #cheap so do a lot - self._chitab = ClassCosmo.z_of_r(self._ztabinchi)[0] + # self._chitab = ClassCosmo.z_of_r(self._ztabinchi)[0] + # self.zfofRint = interp1d(self._chitab, self._ztabinchi) + self._chitab, self._Hztab = ClassCosmo.z_of_r(self._ztabinchi) #chi and dchi/dz self.zfofRint = interp1d(self._chitab, self._ztabinchi) + self.chiofzint = interp1d(self._ztabinchi,self._chitab) + self.Hofzint = interp1d(self._ztabinchi,self._Hztab) _thermo = ClassCosmo.get_thermodynamics() self.Tadiabaticint = interp1d(_thermo['z'], _thermo['Tb [K]']) @@ -107,7 +112,7 @@ def __init__(self, CosmoParams_input, ClassCosmo): #HMF-related constants if(self.Flag_emulate_21cmfast == False): #standard, best fit ST from Schneider+ - self.a_ST = 0.85 #0.85 to fit 1805.00021, 0.707 for original ST + self.a_ST = 0.707 #OG ST fit, or 0.85 to fit 1805.00021 self.p_ST = 0.3 self.Amp_ST = 0.3222 self.delta_crit_ST = 1.686 @@ -126,7 +131,9 @@ def __init__(self, CosmoParams_input, ClassCosmo): class Astro_Parameters: "Class to pass the astro parameters as input" - def __init__(self, Cosmo_Parameters, astromodel = 0, epsstar = 0.1, alphastar = 0.5, betastar = -0.5, Mc = 3e11, fesc10 = 0.1, alphaesc = 0.0, L40_xray = 3.0, E0_xray = 500., alpha_xray = -1.0, Emax_xray_norm=2000, Nalpha_lyA = 9690, Mturn_fixed = None, FLAG_MTURN_SHARP= False, accretion_model = 1): + def __init__(self, Cosmo_Parameters, astromodel = 0, epsstar = 0.1, dlog10epsstardz = 0.0, alphastar = 0.5, betastar = -0.5, Mc = 3e11, fesc10 = 0.1, alphaesc = 0.0, \ + L40_xray = 3.0, E0_xray = 500., alpha_xray = -1.0, Emax_xray_norm=2000, Nalpha_lyA = 9690, Mturn_fixed = None, FLAG_MTURN_SHARP= False, \ + accretion_model = 0, sigmaUV=0.5, C0dust = 4.43, C1dust = 1.99): if(Cosmo_Parameters.Flag_emulate_21cmfast==True and astromodel == 0): @@ -138,9 +145,12 @@ def __init__(self, Cosmo_Parameters, astromodel = 0, epsstar = 0.1, alphastar = #SFR(Mh) parameters: self.epsstar = epsstar #epsilon_* = f* at Mc + self.dlog10epsstardz = dlog10epsstardz #dlog10epsilon/dz + self._zpivot = 8.0 #fixed, at which z we evaluate eps and dlogeps/dz self.alphastar = alphastar #powerlaw index for lower masses self.betastar = betastar #powerlaw index for higher masses, only for model 0 self.Mc = Mc # mass at which the power law cuts, only for model 0 + self.sigmaUV = sigmaUV #stochasticity (gaussian rms) in the halo-galaxy connection P(MUV | Mh) - TODO: only used in UVLF not sfrd if self.astromodel == 0: #GALUMI-like self.accretion_model = accretion_model #0 = exponential, 1= EPS #choose the accretion model. Default = EPS @@ -191,6 +201,10 @@ def __init__(self, Cosmo_Parameters, astromodel = 0, epsstar = 0.1, alphastar = self.Mturn_fixed = Mturn_fixed self.FLAG_MTURN_SHARP = FLAG_MTURN_SHARP #whether to do sharp cut at Mturn_fixed or regular exponential cutoff. Only active if FLAG_MTURN_FIXED and turned on by hand. + #dust parameters for UVLFs: + self.C0dust, self.C1dust = C0dust, C1dust #4.43, 1.99 is Meurer99; 4.54, 2.07 is Overzier01 + self._kappaUV = 1.15e-28 #SFR/LUV, value from Madau+Dickinson14, fully degenerate with epsilon + diff --git a/zeus21/sfrd.py b/zeus21/sfrd.py index 7790cc2..178a519 100644 --- a/zeus21/sfrd.py +++ b/zeus21/sfrd.py @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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