RP16 halo accretion

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, minMUV = None, DUST_FLAG=True, RETURNBIAS = False, RETURNWEIGHTS = False):
     'Binned UVLF in units of 1/Mpc^3/mag, for bins at <zcenter> 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,22 @@ 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 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 minMUV
+
+    if RETURNWEIGHTS:
+        return weights
+
+
     UVLF_filtered = np.trapz(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1)
 
     if(Astro_Parameters.USE_POPIII==False):

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

zeus21/sfrd.py

@@ -9,7 +9,7 @@
 JHU - July 2024
 
 Edited by Emily Bregou
-UT Austin - October 2025
+UT Austin - April 2026
 """
 
 from . import cosmology
@@ -792,10 +792,12 @@ 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 +811,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