From 37633b63103ffd961a85211577853e051972aa60 Mon Sep 17 00:00:00 2001 From: Emily Bregou Date: Wed, 11 Feb 2026 14:20:19 -0600 Subject: [PATCH 1/4] added halo mass depenent min_MUV to P(MUV|Mh) functions for calculating UVLFs & bias --- zeus21/UVLFs.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/zeus21/UVLFs.py b/zeus21/UVLFs.py index b654c53..b353c9b 100644 --- a/zeus21/UVLFs.py +++ b/zeus21/UVLFs.py @@ -8,8 +8,8 @@ Edited by Hector Afonso G. Cruz JHU - July 2024 -Bug fix by Emily Bregou -UT Austin - June 2025 +Edited by Emily Bregou +UT Austin - February 2026 """ from . import cosmology @@ -35,7 +35,7 @@ def MUV_of_SFR(SFRtab, kappaUV): #and combine to get UVLF: -def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwidth, MUVcenters, MUVwidths, DUST_FLAG=True, RETURNBIAS = False): +def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwidth, MUVcenters, MUVwidths, min_MUV = None, DUST_FLAG=True, RETURNBIAS = False, RETURNWEIGHTS = 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): @@ -57,7 +57,7 @@ def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwi MUVbarlist = np.fmin(MUVbarlist,constants._MAGMAX) - if(RETURNBIAS==True): # weight by bias + 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) @@ -80,8 +80,21 @@ 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) - + + # Cut distributions based on min_MUV (user-input, halo mass depenent): + if min_MUV is None: + min_MUV = np.full_like(-100, HMF_interpolator.Mhtab) + 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)))[:,None] # Renormalize distributions based on the portion cut off by min_MUV + + if RETURNWEIGHTS: + return weights + + UVLF_filtered = np.trapz(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1) if(Astro_Parameters.USE_POPIII==False): From 9eabc5dbb52d7ef29a05c015efb72f601617433d Mon Sep 17 00:00:00 2001 From: Emily Bregou Date: Wed, 11 Feb 2026 15:51:32 -0600 Subject: [PATCH 2/4] fixed bug --- zeus21/UVLFs.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/zeus21/UVLFs.py b/zeus21/UVLFs.py index b353c9b..f7905f9 100644 --- a/zeus21/UVLFs.py +++ b/zeus21/UVLFs.py @@ -35,7 +35,7 @@ def MUV_of_SFR(SFRtab, kappaUV): #and combine to get UVLF: -def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwidth, MUVcenters, MUVwidths, min_MUV = None, DUST_FLAG=True, RETURNBIAS = False, RETURNWEIGHTS = False): +def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwidth, MUVcenters, MUVwidths, minMUV = None, DUST_FLAG=True, RETURNBIAS = False, RETURNWEIGHTS = 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): @@ -81,15 +81,16 @@ 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) - # Cut distributions based on min_MUV (user-input, halo mass depenent): - if min_MUV is None: - min_MUV = np.full_like(-100, HMF_interpolator.Mhtab) - x_min = (min_MUV - currMUV)/(np.sqrt(2) * sigmaUV) + # Cut distributions based on minMUV (user-input, halo mass depenent) + #MUVmin can be set to be MUV_of_SFR(max_SFR, Astro_Parameters._kappaUV), with max_SFR = Mstar/min_t & Mstar = fb*Mh + if minMUV is None: + minMUV = np.full_like(HMF_interpolator.Mhtab, -100) + x_min = (minMUV - 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)))[:,None] # Renormalize distributions based on the portion cut off by min_MUV + weights = weights_unnormalized/ (0.5*(1-erf(x_min)))[:,None] # Renormalize distributions based on the portion cut off by minMUV if RETURNWEIGHTS: return weights From 76eedce7d267616bd467d13ee69e4fea4fe6dc29 Mon Sep 17 00:00:00 2001 From: Emily Bregou Date: Mon, 13 Apr 2026 15:58:21 -0500 Subject: [PATCH 3/4] Implemented Rodriguez-Puebla 16 accretion --- zeus21/inputs.py | 10 +++++++--- zeus21/sfrd.py | 20 ++++++++++++++++---- 2 files changed, 23 insertions(+), 7 deletions(-) diff --git a/zeus21/inputs.py b/zeus21/inputs.py index 4ceadc7..d705c71 100644 --- a/zeus21/inputs.py +++ b/zeus21/inputs.py @@ -7,6 +7,9 @@ Edited by Hector Afonso G. Cruz JHU - July 2024 + +Edited by Emily Bregou +UT Austin - March 2026 """ from . import constants @@ -220,8 +223,9 @@ class Astro_Parameters: def __init__(self, UserParams, Cosmo_Parameters, astromodel = 0, - accretion_model = 0, - + accretion_model = 'Exp', # Options are exponential (Exp), extended Press-Schechter (EPS), + # or a fitting function from Nbody simulations (RP16) + alphastar = 0.5, betastar = -0.5, epsstar = 0.1, @@ -324,7 +328,7 @@ def __init__(self, UserParams, Cosmo_Parameters, self.fstarmax = 1.0 #where we cap it if self.astromodel == 0: #GALUMI-like - self.accretion_model = accretion_model #0 = exponential, 1= EPS #choose the accretion model. Default = EPS + self.accretion_model = accretion_model #choose the accretion model. Default = Exp. Exp = exponential, EPS = extended Press-Schechter, Yung = Yung+24 fitting function elif self.astromodel == 1: #21cmfast-like, ignores Mc and beta and has a t* later in SFR() self.tstar = 0.5 self.fstar10 = self.epsstar diff --git a/zeus21/sfrd.py b/zeus21/sfrd.py index b36c4f9..9f206d5 100644 --- a/zeus21/sfrd.py +++ b/zeus21/sfrd.py @@ -9,7 +9,7 @@ JHU - July 2024 Edited by Emily Bregou -UT Austin - October 2025 +UT Austin - March 2026 """ from . import cosmology @@ -792,10 +792,11 @@ def dMh_dt(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, massVector, z): Mh = massVector if(Astro_Parameters.astromodel == False): #GALLUMI-like - if(Astro_Parameters.accretion_model == False): #exponential accretion + if(Astro_Parameters.accretion_model == 'Exp'): #exponential accretion dMhdz = massVector * constants.ALPHA_accretion_exponential + Mhdot = dMhdz*cosmology.Hubinvyr(Cosmo_Parameters,z)*(1.0+z) - elif(Astro_Parameters.accretion_model == True): #EPS accretion + elif(Astro_Parameters.accretion_model == 'EPS'): #EPS accretion Mh2 = Mh * constants.EPSQ_accretion indexMh2low = Mh2 < Mh.flatten()[0] @@ -809,10 +810,21 @@ def dMh_dt(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, massVector, z): dzgrow = z*0.01 dgrowthdz = (cosmology.growth(Cosmo_Parameters,z+dzgrow) - cosmology.growth(Cosmo_Parameters,z-dzgrow))/(2.0 * dzgrow) dMhdz = - Mh * np.sqrt(2/np.pi)/np.sqrt(sigmaMh2**2 - sigmaMh**2) *dgrowthdz/growth * Cosmo_Parameters.delta_crit_ST + + Mhdot = dMhdz*cosmology.Hubinvyr(Cosmo_Parameters,z)*(1.0+z) + + elif(Astro_Parameters.accretion_model == 'RP16'): # Fitting function to Rodríguez-Puebla+16 N-body simulations (eq. 11, dynamically + # averaged parameters from table 2) + a = (1+z)**-1 + beta = 10**(2.73-(1.828*a)+(0.654*a**2)) + alpha = 1 + (0.329*a) - (0.206*a**2) + + # factors of h are accounted for to give units of M_sun/year for halo masses in units of M_sun: + Mhdot = beta * (Mh/1e12)**alpha * cosmology.Hub(Cosmo_Parameters, z) / (100*Cosmo_Parameters.h_fid) else: print("ERROR! Have to choose an accretion model in Astro_Parameters (accretion_model)") - Mhdot = dMhdz*cosmology.Hubinvyr(Cosmo_Parameters,z)*(1.0+z) + return Mhdot elif(Astro_Parameters.astromodel == True): #21cmfast-like From 7b4ecce12e3e3c279c89d1e0f475d61aa47eb5ed Mon Sep 17 00:00:00 2001 From: Emily Bregou Date: Mon, 13 Apr 2026 16:03:48 -0500 Subject: [PATCH 4/4] minor formatting changes --- zeus21/sfrd.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/zeus21/sfrd.py b/zeus21/sfrd.py index 9f206d5..7a374b6 100644 --- a/zeus21/sfrd.py +++ b/zeus21/sfrd.py @@ -9,7 +9,7 @@ JHU - July 2024 Edited by Emily Bregou -UT Austin - March 2026 +UT Austin - April 2026 """ from . import cosmology @@ -792,6 +792,7 @@ def dMh_dt(Astro_Parameters, Cosmo_Parameters, HMF_interpolator, massVector, z): Mh = massVector if(Astro_Parameters.astromodel == False): #GALLUMI-like + if(Astro_Parameters.accretion_model == 'Exp'): #exponential accretion dMhdz = massVector * constants.ALPHA_accretion_exponential Mhdot = dMhdz*cosmology.Hubinvyr(Cosmo_Parameters,z)*(1.0+z)