diff --git a/requirements.txt b/requirements.txt index 54f6d49..9b133ad 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ pytest -numpy +numpy>=2.0 scipy mcfit numexpr @@ -7,4 +7,5 @@ astropy powerbox pyfftw sphinx -myst_parser \ No newline at end of file +myst_parser +tqdm diff --git a/setup.py b/setup.py index 898653c..7a64876 100755 --- a/setup.py +++ b/setup.py @@ -14,7 +14,7 @@ packages=['zeus21'], long_description=open('README.md').read(), install_requires=[ - "numpy", + "numpy>=2.0", "scipy", "mcfit", "classy", diff --git a/tests/test_UVLFs.py b/tests/test_UVLFs.py index e63322b..7c34c47 100644 --- a/tests/test_UVLFs.py +++ b/tests/test_UVLFs.py @@ -138,4 +138,67 @@ def test_UVLF_binned(): assert uvlf_nodust.shape == (3,) # Without dust, we expect different values than with dust - assert not np.array_equal(uvlf, uvlf_nodust) \ No newline at end of file + assert not np.array_equal(uvlf, uvlf_nodust) + + +def test_UVLF_binned_with_min_t_formation(): + """Test that min_t_formation_Myr produces finite outputs and suppresses the bright end. + + When sigmaUV is large, scatter can push small halos into unphysically bright bins. + Setting min_t_formation_Myr places a physical upper limit on each halo's SFR based on + its maximum stellar mass (all baryons converted to stars) and the minimum formation time. + This should suppress the very bright end of the UVLF without affecting the faint end. + """ + UserParams = zeus21.User_Parameters() + CosmoParams_input = zeus21.Cosmo_Parameters_Input(kmax_CLASS=10., zmax_CLASS=20.) + ClassyCosmo = zeus21.runclass(CosmoParams_input) + CosmoParams = zeus21.Cosmo_Parameters(UserParams, CosmoParams_input, ClassyCosmo) + HMFintclass = zeus21.HMF_interpolator(UserParams, CosmoParams, ClassyCosmo) + + # Use a large sigmaUV to create unphysical scatter into the bright end + large_sigmaUV = 2.0 + min_t_Myr = 10.0 + + # AstroParams with the physicality cutoff applied + AstroParams_cut = zeus21.Astro_Parameters( + UserParams, CosmoParams, + sigmaUV=large_sigmaUV, + min_t_formation_Myr=min_t_Myr + ) + + # AstroParams without the cutoff (default None) + AstroParams_nocut = zeus21.Astro_Parameters( + UserParams, CosmoParams, + sigmaUV=large_sigmaUV + ) + + z_center = 6.0 + z_width = 0.5 + # Include a very bright bin (-25) where small-halo scatter is cut off, + # a typical bin (-20), and a faint bin (-15) that should be unaffected + MUV_centers = np.array([-25.0, -20.0, -15.0]) + MUV_widths = np.full_like(MUV_centers, 1.0) + + uvlf_cut = UVLF_binned( + AstroParams_cut, CosmoParams, HMFintclass, + z_center, z_width, MUV_centers, MUV_widths, + DUST_FLAG=False, RETURNBIAS=False + ) + uvlf_nocut = UVLF_binned( + AstroParams_nocut, CosmoParams, HMFintclass, + z_center, z_width, MUV_centers, MUV_widths, + DUST_FLAG=False, RETURNBIAS=False + ) + + # Output must be finite (no NaNs or Infs) with the cutoff applied + assert np.all(np.isfinite(uvlf_cut)), "UVLF with min_t_formation_Myr cutoff contains NaN or Inf values" + + # All values must be non-negative + assert np.all(uvlf_cut >= 0.0), "UVLF with min_t_formation_Myr cutoff contains negative values" + + # The cutoff should suppress the very bright end: small halos that could not + # physically produce MUV=-25 galaxies (min_MUV~-18.7 for 1e8 Msun with t_min=10 Myr) + # no longer contribute via scatter, so the bright-end UVLF should be lower + assert uvlf_cut[0] < uvlf_nocut[0], ( + "min_t_formation_Myr cutoff should suppress the very bright end (MUV=-25) of the UVLF" + ) \ No newline at end of file diff --git a/tests/test_inputs.py b/tests/test_inputs.py index 4a3102d..606f143 100644 --- a/tests/test_inputs.py +++ b/tests/test_inputs.py @@ -96,23 +96,23 @@ def test_inputs(): #test Pop II Xray SED Energylisttest = np.logspace(2,np.log10(AstroParams.Emax_xray_norm),100) SEDXtab_test = AstroParams.SED_XRAY(Energylisttest, 2) #same in both models - normalization_XraySED = np.trapz(Energylisttest * SEDXtab_test,Energylisttest) + normalization_XraySED = np.trapezoid(Energylisttest * SEDXtab_test,Energylisttest) assert( normalization_XraySED == pytest.approx(1.0, 0.05) ) #5% is enough here #test Pop III Xray SED SEDXtab_test = AstroParams.SED_XRAY(Energylisttest, 3) #same in both models - normalization_XraySED = np.trapz(Energylisttest * SEDXtab_test,Energylisttest) + normalization_XraySED = np.trapezoid(Energylisttest * SEDXtab_test,Energylisttest) assert( normalization_XraySED == pytest.approx(1.0, 0.05) ) #5% is enough here #test Pop II LyA SED nulisttest = np.linspace(zeus21.constants.freqLyA, zeus21.constants.freqLyCont, 100) SEDLtab_test = AstroParams.SED_LyA(nulisttest, 2) #same in both models - normalization_LyASED = np.trapz(SEDLtab_test,nulisttest) + normalization_LyASED = np.trapezoid(SEDLtab_test,nulisttest) assert( normalization_LyASED == pytest.approx(1.0, 0.05) ) #5% is enough here #test Pop III LyA SED nulisttest = np.linspace(zeus21.constants.freqLyA, zeus21.constants.freqLyCont, 100) SEDLtab_test = AstroParams.SED_LyA(nulisttest, 3) #same in both models - normalization_LyASED = np.trapz(SEDLtab_test,nulisttest) + normalization_LyASED = np.trapezoid(SEDLtab_test,nulisttest) assert( normalization_LyASED == pytest.approx(1.0, 0.05) ) #5% is enough here diff --git a/zeus21/UVLFs.py b/zeus21/UVLFs.py index b654c53..7fd1d22 100644 --- a/zeus21/UVLFs.py +++ b/zeus21/UVLFs.py @@ -80,9 +80,26 @@ def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwi 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) + + if (Astro_Parameters.min_t_formation_Myr == None): + min_MUV = -100.0 # essentially no cutoff, since the scatter is large at low masses and can cause numerical issues if we try to integrate over unphysically bright galaxies there. This is just a numerical cutoff, not a physical one, and the exact value doesn't matter much since the scatter is large there anyway. + else: + Mstarmax = HMF_interpolator.Mhtab * Cosmo_Parameters.OmegaB /Cosmo_Parameters.OmegaM #max stellar mass in each halo, if all baryons turned to stars + _tmaxSFR = Astro_Parameters.min_t_formation_Myr * 1e6 #arbitrary timescale to determine max SFR in yrs + SFRmax = Mstarmax / (_tmaxSFR) + min_MUV = MUV_of_SFR(SFRmax, Astro_Parameters._kappaUV) #min MUV in each halo, if all baryons turned to stars at max SFR for 10 Myr. This is a very rough cutoff to avoid unphysically small MUVs (bright galaxies) at low masses, which can cause numerical issues since the scatter is large there. It's not a physical cutoff, just a numerical one. The exact value doesn't matter much since the scatter is large there anyway, but it prevents the code from trying to integrate over unphysically bright galaxies in low-mass halos. + x_min = (min_MUV - currMUV)/(np.sqrt(2) * sigmaUV) + xhi_cut = np.fmax(xhi, x_min) + xlo_cut = np.fmax(xlo, x_min) + + weights_unnormalized = (erf(xhi_cut) - erf(xlo_cut)).T/(2.0 * MUVwidths) + weights = weights_unnormalized/ (0.5*(1-erf(x_min)+1e-6))[:,None] # Renormalize distributions based on the portion cut off by min_MUV + + ### Standard as usual, no cuts: + # weights = (erf(xhi) - erf(xlo)).T/(2.0 * MUVwidths) #comment to myself, this 2 in denominator is correct here, nothing to do with the MUVwidths/2 a few lines above - UVLF_filtered = np.trapz(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1) + UVLF_filtered = np.trapezoid(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1) + if(Astro_Parameters.USE_POPIII==False): return UVLF_filtered @@ -98,7 +115,7 @@ def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwi xlo = np.subtract.outer(MUVcutlo, MUVbarlist_III)/(np.sqrt(2) * sigmaUV) weights = (erf(xhi) - erf(xlo)).T/(2.0 * MUVwidths) - UVLF_filtered_III = np.trapz(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1) + UVLF_filtered_III = np.trapezoid(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1) return UVLF_filtered, UVLF_filtered_III diff --git a/zeus21/correlations.py b/zeus21/correlations.py index d6ab1c2..c79d72c 100644 --- a/zeus21/correlations.py +++ b/zeus21/correlations.py @@ -382,7 +382,7 @@ def __init__(self, User_Parameters, Cosmo_Parameters, Astro_Parameters, ClassCos #the z>zmax part of the integral we do aside. Assume Tk=Tadiabatic from CLASS. _zlisthighz_ = np.linspace(T21_coefficients.zintegral[-1], 99., 100) #beyond z=100 need to explictly tell CLASS to save growth _dgrowthhighz_ = cosmology.dgrowth_dz(Cosmo_Parameters, _zlisthighz_) - _hizintegral = np.trapz(cosmology.Tadiabatic(Cosmo_Parameters,_zlisthighz_) + _hizintegral = np.trapezoid(cosmology.Tadiabatic(Cosmo_Parameters,_zlisthighz_) /(1+_zlisthighz_)**2 * _dgrowthhighz_, _zlisthighz_) self._betaTad_ = -2./3. * _factor_adi_/self._lingrowthd * (np.cumsum(_integrand_adi[::-1])[::-1] + _hizintegral) #units of Tk_avg. Internal sum goes from high to low z (backwards), minus sign accounts for it properly so it's positive. @@ -1261,7 +1261,7 @@ def get_Pk_from_xi(self, rsinput, xiinput): # # Probdtab = np.exp(-dtab**2/sigmaRRsq/2.0) # -# norm = np.trapz(NionEPS * Probdtab, dtab) +# norm = np.trapezoid(NionEPS * Probdtab, dtab) # NionEPS/=norm # # bindex = min(range(len(NionEPS)), key=lambda i: abs(NionEPS[i]-_invQbar)) diff --git a/zeus21/cosmology.py b/zeus21/cosmology.py index c2500b7..c468dea 100644 --- a/zeus21/cosmology.py +++ b/zeus21/cosmology.py @@ -75,13 +75,13 @@ def runclass(CosmologyIn): theta_b = velTransFunc['t_b'] theta_c = velTransFunc['t_cdm'] - sigma_vcb = np.sqrt(np.trapz(CosmologyIn.As * (kVel/0.05)**(CosmologyIn.ns-1) /kVel * (theta_b - theta_c)**2/kVel**2, kVel)) * constants.c_kms + sigma_vcb = np.sqrt(np.trapezoid(CosmologyIn.As * (kVel/0.05)**(CosmologyIn.ns-1) /kVel * (theta_b - theta_c)**2/kVel**2, kVel)) * constants.c_kms ClassCosmo.pars['sigma_vcb'] = sigma_vcb ###HAC: now computing average velocity assuming a Maxwell-Boltzmann distribution of velocities velArr = np.geomspace(0.01, constants.c_kms, 1000) #in km/s vavgIntegrand = (3 / (2 * np.pi * sigma_vcb**2))**(3/2) * 4 * np.pi * velArr**2 * np.exp(-3 * velArr**2 / (2 * sigma_vcb**2)) - ClassCosmo.pars['v_avg'] = np.trapz(vavgIntegrand * velArr, velArr) + ClassCosmo.pars['v_avg'] = np.trapezoid(vavgIntegrand * velArr, velArr) ###HAC: Computing Vcb Power Spectrum ClassCosmo.pars['k_vcb'] = kVel @@ -99,8 +99,8 @@ def runclass(CosmologyIn): j0bessel = lambda x: np.sin(x)/x j2bessel = lambda x: (3 / x**2 - 1) * np.sin(x)/x - 3*np.cos(x)/x**2 - psi0 = 1 / 3 / (sigma_vcb/constants.c_kms)**2 * np.trapz(kVelIntp**2 / 2 / np.pi**2 * p_vcb_intp(np.log(kVelIntp)) * j0bessel(kVelIntp * np.transpose([rVelIntp])), kVelIntp, axis = 1) - psi2 = -2 / 3 / (sigma_vcb/constants.c_kms)**2 * np.trapz(kVelIntp**2 / 2 / np.pi**2 * p_vcb_intp(np.log(kVelIntp)) * j2bessel(kVelIntp * np.transpose([rVelIntp])), kVelIntp, axis = 1) + psi0 = 1 / 3 / (sigma_vcb/constants.c_kms)**2 * np.trapezoid(kVelIntp**2 / 2 / np.pi**2 * p_vcb_intp(np.log(kVelIntp)) * j0bessel(kVelIntp * np.transpose([rVelIntp])), kVelIntp, axis = 1) + psi2 = -2 / 3 / (sigma_vcb/constants.c_kms)**2 * np.trapezoid(kVelIntp**2 / 2 / np.pi**2 * p_vcb_intp(np.log(kVelIntp)) * j2bessel(kVelIntp * np.transpose([rVelIntp])), kVelIntp, axis = 1) k_eta, P_eta = mcfit.xi2P(rVelIntp, l=0, lowring = True)((6 * psi0**2 + 3 * psi2**2), extrap = False) diff --git a/zeus21/inputs.py b/zeus21/inputs.py index e838903..65184d7 100644 --- a/zeus21/inputs.py +++ b/zeus21/inputs.py @@ -270,7 +270,8 @@ def __init__(self, UserParams, Cosmo_Parameters, A_vcb = 1.0, beta_vcb = 1.8, - quadratic_SFRD_lognormal = False # Sarah Libanore, use second order in lognormal + quadratic_SFRD_lognormal = False, # Sarah Libanore, use second order in lognormal + min_t_formation_Myr = None ): @@ -412,7 +413,13 @@ def __init__(self, UserParams, Cosmo_Parameters, if Cosmo_Parameters.Flag_emulate_21cmfast: print('Quadratic SFRD not yet implemented when Flag_emulate_21cmfast = True; the code will use quadratic_SFRD_lognormal = False') self.quadratic_SFRD_lognormal = False - + + if min_t_formation_Myr is not None: + if (not np.isscalar(min_t_formation_Myr) + or not np.isfinite(min_t_formation_Myr) + or min_t_formation_Myr <= 0): + raise ValueError("min_t_formation_Myr must be None or a strictly positive finite number.") + self.min_t_formation_Myr = min_t_formation_Myr #Minimum formation time of galaxies in Myr for UVLF, sets a minimum M*dot = M*/t_formation with fstar = 1 def SED_XRAY(self, En, pop = 0): #pop set to zero as default, but it must be set to either 2 or 3 diff --git a/zeus21/sfrd.py b/zeus21/sfrd.py index 4471784..6819cc5 100644 --- a/zeus21/sfrd.py +++ b/zeus21/sfrd.py @@ -8,6 +8,9 @@ Edited by Hector Afonso G. Cruz JHU - July 2024 +Edited by Emily Bregou +UT Austin - October 2025 + Edited by Sarah Libanore BGU - July 2025 @@ -135,14 +138,14 @@ def __init__(self, User_Parameters, Cosmo_Parameters, ClassCosmo, Astro_Paramete zSFRD, mArray = np.meshgrid(zSFRDflat, HMF_interpolator.Mhtab, indexing = 'ij', sparse = True) J21LW_interp = interpolate.interp1d(zSFRDflat, np.zeros_like(zSFRDflat), kind = 'linear', bounds_error = False, fill_value = 0,) #no LW background. Controls only Mmol() function, NOT the individual Pop II and III LW background - SFRD_II_avg = np.trapz(SFRD_II_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mArray, zSFRD, zSFRD), HMF_interpolator.logtabMh, axis = 1) #never changes with J_LW + SFRD_II_avg = np.trapezoid(SFRD_II_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mArray, zSFRD, zSFRD), HMF_interpolator.logtabMh, axis = 1) #never changes with J_LW SFRD_II_interp = interpolate.interp1d(zSFRDflat, SFRD_II_avg, kind = 'cubic', bounds_error = False, fill_value = 0,) J21LW_II = 1e21 * J_LW(Astro_Parameters, Cosmo_Parameters, SFRD_II_avg, zSFRDflat, 2) #this never changes; only Pop III Quanties change self.J_21_LW_II = interpolate.interp1d(zSFRDflat, J21LW_II, kind = 'cubic')(self.zintegral) #different from J21LW_interp if Astro_Parameters.USE_POPIII == True: - SFRD_III_Iter_Matrix = [np.trapz(SFRD_III_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mArray, J21LW_interp, zSFRD, zSFRD, ClassCosmo.pars['v_avg']), HMF_interpolator.logtabMh, axis = 1)] #changes with each iteration + SFRD_III_Iter_Matrix = [np.trapezoid(SFRD_III_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mArray, J21LW_interp, zSFRD, zSFRD, ClassCosmo.pars['v_avg']), HMF_interpolator.logtabMh, axis = 1)] #changes with each iteration errorTolerance = 0.001 # 0.1 percent accuracy recur_iterate_Flag = True @@ -150,7 +153,7 @@ def __init__(self, User_Parameters, Cosmo_Parameters, ClassCosmo, Astro_Paramete J21LW_III_iter = 1e21 * J_LW(Astro_Parameters, Cosmo_Parameters, SFRD_III_Iter_Matrix[-1], zSFRDflat, 3) J21LW_interp = interpolate.interp1d(zSFRDflat, J21LW_II + J21LW_III_iter, kind = 'linear', fill_value = 0, bounds_error = False) - SFRD_III_avg_n = np.trapz(SFRD_III_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mArray, J21LW_interp, zSFRD, zSFRD, ClassCosmo.pars['v_avg']), HMF_interpolator.logtabMh, axis = 1) + SFRD_III_avg_n = np.trapezoid(SFRD_III_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mArray, J21LW_interp, zSFRD, zSFRD, ClassCosmo.pars['v_avg']), HMF_interpolator.logtabMh, axis = 1) SFRD_III_Iter_Matrix.append(SFRD_III_avg_n) if max(SFRD_III_Iter_Matrix[-1]/SFRD_III_Iter_Matrix[-2]) < 1.0 + errorTolerance and min(SFRD_III_Iter_Matrix[-1]/SFRD_III_Iter_Matrix[-2]) > 1.0 - errorTolerance: @@ -176,8 +179,8 @@ def __init__(self, User_Parameters, Cosmo_Parameters, ClassCosmo, Astro_Paramete zpTable, tempTable, mTable = np.meshgrid(self.zintegral, self.Rtabsmoo, HMF_interpolator.Mhtab, indexing = 'ij', sparse = True) zppTable = self.zGreaterMatrix.reshape((len(self.zintegral), len(self.Rtabsmoo), 1)) - self.SFRDbar2D_II = np.trapz(SFRD_II_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mTable, zppTable, zpTable), HMF_interpolator.logtabMh, axis = 2) - self.SFRDbar2D_III = np.trapz(SFRD_III_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mTable, J21LW_interp, zppTable, zpTable, ClassCosmo.pars['v_avg']), HMF_interpolator.logtabMh, axis = 2) + self.SFRDbar2D_II = np.trapezoid(SFRD_II_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mTable, zppTable, zpTable), HMF_interpolator.logtabMh, axis = 2) + self.SFRDbar2D_III = np.trapezoid(SFRD_III_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mTable, J21LW_interp, zppTable, zpTable, ClassCosmo.pars['v_avg']), HMF_interpolator.logtabMh, axis = 2) self.SFRDbar2D_II[np.isnan(self.SFRDbar2D_II)] = 0.0 self.SFRDbar2D_III[np.isnan(self.SFRDbar2D_III)] = 0.0 @@ -238,13 +241,13 @@ def __init__(self, User_Parameters, Cosmo_Parameters, ClassCosmo, Astro_Paramete ######## # Compute SFRD quantities - SFRD_II_dR = np.trapz(integrand_II, HMF_interpolator.logtabMh, axis = 2) + SFRD_II_dR = np.trapezoid(integrand_II, HMF_interpolator.logtabMh, axis = 2) # SarahLibanore: to compute reionization - niondot_II_dR = np.trapz(integrand_II*fesctab_II[None, None, :, None], HMF_interpolator.logtabMh, axis = 2) + niondot_II_dR = np.trapezoid(integrand_II*fesctab_II[None, None, :, None], HMF_interpolator.logtabMh, axis = 2) # SarahLibanore: compute quantities in Lagrangian space to get gamma in Lagrangian space - SFRD_II_dR_Lag = np.trapz(integrand_II_Lag, HMF_interpolator.logtabMh, axis = 2) - niondot_II_dR_Lag = np.trapz(integrand_II_Lag*fesctab_II[None, None, :, None], HMF_interpolator.logtabMh, axis = 2) + SFRD_II_dR_Lag = np.trapezoid(integrand_II_Lag, HMF_interpolator.logtabMh, axis = 2) + niondot_II_dR_Lag = np.trapezoid(integrand_II_Lag*fesctab_II[None, None, :, None], HMF_interpolator.logtabMh, axis = 2) ### if Astro_Parameters.USE_POPIII == True: @@ -253,9 +256,9 @@ def __init__(self, User_Parameters, Cosmo_Parameters, ClassCosmo, Astro_Paramete elif(Cosmo_Parameters.Flag_emulate_21cmfast==True): integrand_III = PS_HMF_corr * SFR_III(Astro_Parameters, Cosmo_Parameters, ClassCosmo, HMF_interpolator, mArray, J21LW_interp, zGreaterArray, zGreaterArray, ClassCosmo.pars['v_avg']) * mArray - SFRD_III_dR = np.trapz(integrand_III, HMF_interpolator.logtabMh, axis = 2) + SFRD_III_dR = np.trapezoid(integrand_III, HMF_interpolator.logtabMh, axis = 2) # SarahLibanore: reionization - niondot_III_dR = np.trapz(integrand_III*fesctab_III[None, None, :, None], HMF_interpolator.logtabMh, axis = 2) + niondot_III_dR = np.trapezoid(integrand_III*fesctab_III[None, None, :, None], HMF_interpolator.logtabMh, axis = 2) else: SFRD_III_dR = np.zeros_like(SFRD_II_dR) @@ -373,7 +376,7 @@ def __init__(self, User_Parameters, Cosmo_Parameters, ClassCosmo, Astro_Paramete else: print("ERROR: Need to set FLAG_EMULATE_21CMFAST at True or False in the self.gamma_index2D calculation.") - SFRD_III_dR_V = np.trapz(integrand_III, HMF_interpolator.logtabMh, axis = 2) + SFRD_III_dR_V = np.trapezoid(integrand_III, HMF_interpolator.logtabMh, axis = 2) SFRDIII_Ratio = SFRD_III_dR_V / SFRD_III_dR_V[:,:,len(vAvg_array)//2].reshape((len(self.zintegral), len(self.Rtabsmoo), 1)) SFRDIII_Ratio[np.isnan(SFRDIII_Ratio)] = 0.0 @@ -497,7 +500,7 @@ def __init__(self, User_Parameters, Cosmo_Parameters, ClassCosmo, Astro_Paramete opticalDepthIntegrand = 1 / cosmology.HubinvMpc(Cosmo_Parameters, zPPCube) / (1+zPPCube) * sigmatot * cosmology.n_H(Cosmo_Parameters, zPPCube) * constants.Mpctocm #this uses atom fractions of 1 for HI and x_He for HeI # opticalDepthIntegrand = 1 / cosmology.HubinvMpc(Cosmo_Parameters, zPPCube) / (1+zPPCube) * sigmatot * cosmology.n_baryon(Cosmo_Parameters, zPPCube) * constants.Mpctocm - tauCube = np.trapz(opticalDepthIntegrand, zPPCube, axis = 3) + tauCube = np.trapezoid(opticalDepthIntegrand, zPPCube, axis = 3) indextautoolarge = np.array(tauCube>=Xrays.TAUMAX) tauCube[indextautoolarge] = Xrays.TAUMAX @@ -655,8 +658,8 @@ def __init__(self, User_Parameters, Cosmo_Parameters, ClassCosmo, Astro_Paramete integrand_II_table = SFRD_II_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mArray, zArray, zArray) integrand_III_table = SFRD_III_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, mArray, J21LW_interp, zArray, zArray, ClassCosmo.pars['v_avg']) - self.niondot_avg_II = Astro_Parameters.N_ion_perbaryon_II/cosmology.rho_baryon(Cosmo_Parameters,0.) * np.trapz(integrand_II_table * fesctab_II, HMF_interpolator.logtabMh, axis = 1) - self.niondot_avg_III = Astro_Parameters.N_ion_perbaryon_II/cosmology.rho_baryon(Cosmo_Parameters,0.) * np.trapz(integrand_III_table * fesctab_III, HMF_interpolator.logtabMh, axis = 1) + self.niondot_avg_II = Astro_Parameters.N_ion_perbaryon_II/cosmology.rho_baryon(Cosmo_Parameters,0.) * np.trapezoid(integrand_II_table * fesctab_II, HMF_interpolator.logtabMh, axis = 1) + self.niondot_avg_III = Astro_Parameters.N_ion_perbaryon_II/cosmology.rho_baryon(Cosmo_Parameters,0.) * np.trapezoid(integrand_III_table * fesctab_III, HMF_interpolator.logtabMh, axis = 1) self.niondot_avg = self.niondot_avg_II + self.niondot_avg_III if(Cosmo_Parameters.Flag_emulate_21cmfast==False): #regular calculation, integrating over time and accounting for recombinations in the exponent @@ -705,7 +708,7 @@ def tau_reio(Cosmo_Parameters, T21_coefficients): _nelistlowz = cosmology.n_H(Cosmo_Parameters,_zlistlowz)*(1.0 + Cosmo_Parameters.x_He + Cosmo_Parameters.x_He * np.heaviside(constants.zHeIIreio - _zlistlowz,0.5)) # _nelistlowz = cosmology.n_baryon(Cosmo_Parameters,_zlistlowz)*(Cosmo_Parameters.f_H + Cosmo_Parameters.f_He + Cosmo_Parameters.f_He * np.heaviside(constants.zHeIIreio - _zlistlowz,0.5)) _distlistlowz = 1.0/cosmology.HubinvMpc(Cosmo_Parameters,_zlistlowz)/(1+_zlistlowz) - _lowzint = constants.sigmaT * np.trapz(_nelistlowz*_distlistlowz,_zlistlowz) * constants.Mpctocm + _lowzint = constants.sigmaT * np.trapezoid(_nelistlowz*_distlistlowz,_zlistlowz) * constants.Mpctocm _zlisthiz = T21_coefficients.zintegral @@ -713,7 +716,7 @@ def tau_reio(Cosmo_Parameters, T21_coefficients): # _nelistlhiz = cosmology.n_baryon(Cosmo_Parameters,_zlisthiz) * (1.0 - T21_coefficients.xHI_avg) _distlisthiz = 1.0/cosmology.HubinvMpc(Cosmo_Parameters,_zlisthiz)/(1+_zlisthiz) - _hizint = constants.sigmaT * np.trapz(_nelistlhiz*_distlisthiz,_zlisthiz) * constants.Mpctocm + _hizint = constants.sigmaT * np.trapezoid(_nelistlhiz*_distlisthiz,_zlisthiz) * constants.Mpctocm return(_lowzint + _hizint) @@ -751,10 +754,10 @@ def Mmol(Astro_Parameters, Cosmo_Parameters, J21LW_interp, z, vCB): def fstarofz(Astro_Parameters, Cosmo_Parameters, z, Mhlist): epsstar_ofz = Astro_Parameters.epsstar * 10**(Astro_Parameters.dlog10epsstardz * (z-Astro_Parameters._zpivot) ) if Cosmo_Parameters.Flag_emulate_21cmfast == False: - return 2.0 * Cosmo_Parameters.OmegaB/Cosmo_Parameters.OmegaM * epsstar_ofz\ - /(pow(Mhlist/Astro_Parameters.Mc,- Astro_Parameters.alphastar) + pow(Mhlist/Astro_Parameters.Mc,- Astro_Parameters.betastar) ) + return Cosmo_Parameters.OmegaB/Cosmo_Parameters.OmegaM * np.clip(2.0 * epsstar_ofz\ + /(pow(Mhlist/Astro_Parameters.Mc,- Astro_Parameters.alphastar) + pow(Mhlist/Astro_Parameters.Mc,- Astro_Parameters.betastar) ), 0, 1) elif Cosmo_Parameters.Flag_emulate_21cmfast == True: - return Cosmo_Parameters.OmegaB/Cosmo_Parameters.OmegaM * epsstar_ofz /(pow(Mhlist/Astro_Parameters.Mc,- Astro_Parameters.alphastar)) + return Cosmo_Parameters.OmegaB/Cosmo_Parameters.OmegaM * np.clip(epsstar_ofz /(pow(Mhlist/Astro_Parameters.Mc,- Astro_Parameters.alphastar)), 0, 1) ###HAC: Added fstar for PopIII @@ -795,7 +798,7 @@ def J_LW(Astro_Parameters, Cosmo_Parameters, sfrdIter, z, pop): # integrandLW *= (1+z)**3 / cosmology.Hubinvyr(Cosmo_Parameters,zIntMatrix) / (1+zIntMatrix) #HAC: delete this and comment above back in!!! integrandLW *= Nlw * Elw / massProton / deltaNulw integrandLW = integrandLW * sfrdIterMatrix_LW * (1 /u.Mpc**2).to(1/u.cm**2).value #broadcasting doesn't like augmented assignment operations (like *=) for some reason - return np.trapz(integrandLW, x = zIntMatrix, axis = 0) + return np.trapezoid(integrandLW, x = zIntMatrix, axis = 0) def SFRD_II_integrand(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, massVector, z, z2): @@ -929,7 +932,7 @@ def dSFRDIII_dJ(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, J21LW_inte integrand_III = HMF_curr * SFRtab_currIII * HMF_interpolator.Mhtab integrand_III *= Astro_Parameters.A_LW * Astro_Parameters.beta_LW * J21LW_interp(z)**(Astro_Parameters.beta_LW - 1) integrand_III *= -1 * Mmol_vcb(Astro_Parameters, Cosmo_Parameters, z, Cosmo_Parameters.vcb_avg)/ HMF_interpolator.Mhtab - return np.trapz(integrand_III, HMF_interpolator.logtabMh) + return np.trapezoid(integrand_III, HMF_interpolator.logtabMh) def fesc_II(Astro_Parameters, Mh): diff --git a/zeus21/xrays.py b/zeus21/xrays.py index 6666014..78cdaa1 100644 --- a/zeus21/xrays.py +++ b/zeus21/xrays.py @@ -43,7 +43,7 @@ def optical_depth(self, User_Parameters, Cosmo_Parameters, En,z,zp): # divided by factor of H(z')(1+z') because of variable of integration change from proper distance to redshift integrand = 1.0/HubinvMpc(Cosmo_Parameters, zinttau)/(1+zinttau) * sigmatot * n_H(Cosmo_Parameters, zinttau) * constants.Mpctocm # integrand = 1.0/HubinvMpc(Cosmo_Parameters, zinttau)/(1+zinttau) * sigmatot * n_baryon(Cosmo_Parameters, zinttau) * constants.Mpctocm - taulist = np.trapz(integrand, zinttau, axis=1) + taulist = np.trapezoid(integrand, zinttau, axis=1) #OLD: kept for reference only. # taulist = 1.0*np.zeros_like(Envec) @@ -55,7 +55,7 @@ def optical_depth(self, User_Parameters, Cosmo_Parameters, En,z,zp): # # integrand = 1.0/HubinvMpc(Cosmo_Parameters, zinttau)/(1+zinttau) * sigmatot * n_baryon(Cosmo_Parameters, zinttau) * constants.Mpctocm # - # taulist[iE] = np.trapz(integrand, zinttau) + # taulist[iE] = np.trapezoid(integrand, zinttau) indextautoolarge = np.array(taulist>=self.TAUMAX) taulist [indextautoolarge] = self.TAUMAX