From 8de46e3ffa0133d4836091d989880c5a5752c977 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 29 Apr 2026 16:11:17 +0000 Subject: [PATCH 1/2] Initial plan From d2b1e7847180e12dea7d768676c2295450aaab06 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 29 Apr 2026 16:12:34 +0000 Subject: [PATCH 2/2] Replace np.trapz with np.trapezoid and require numpy>=2.0 Agent-Logs-Url: https://github.com/ZeusCosmo/Zeus21/sessions/125bc937-b800-45fc-b854-32e56c5f75df Co-authored-by: JulianBMunoz <22434409+JulianBMunoz@users.noreply.github.com> --- requirements.txt | 2 +- setup.py | 2 +- tests/test_inputs.py | 8 ++++---- zeus21/UVLFs.py | 4 ++-- zeus21/correlations.py | 4 ++-- zeus21/cosmology.py | 8 ++++---- zeus21/sfrd.py | 38 +++++++++++++++++++------------------- zeus21/xrays.py | 4 ++-- 8 files changed, 35 insertions(+), 35 deletions(-) diff --git a/requirements.txt b/requirements.txt index 54f6d49..bbc242b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ pytest -numpy +numpy>=2.0 scipy mcfit numexpr 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_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..ec08a94 100644 --- a/zeus21/UVLFs.py +++ b/zeus21/UVLFs.py @@ -82,7 +82,7 @@ def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwi 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) + UVLF_filtered = np.trapezoid(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1) if(Astro_Parameters.USE_POPIII==False): return UVLF_filtered @@ -98,7 +98,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/sfrd.py b/zeus21/sfrd.py index 1bcd979..6819cc5 100644 --- a/zeus21/sfrd.py +++ b/zeus21/sfrd.py @@ -138,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 @@ -153,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: @@ -179,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 @@ -241,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: @@ -256,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) @@ -376,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 @@ -500,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 @@ -658,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 @@ -708,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 @@ -716,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) @@ -798,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): @@ -932,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