Merge pull request #40 from ZeusCosmo/implement_max_SFR_physical

Implement max sfr physical

Author: JulianBMunoz

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)
+    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"
+    )

zeus21/UVLFs.py

@@ -80,10 +80,27 @@ 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.trapezoid(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1)
 
+
     if(Astro_Parameters.USE_POPIII==False):
         return UVLF_filtered
     else:

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