From 45c90118bd03d2042bbd9e5edd87e1dae87fda2d Mon Sep 17 00:00:00 2001 From: Julian Munoz Date: Wed, 29 Apr 2026 12:05:50 -0500 Subject: [PATCH 1/6] Added max SFR for UVLFs In case sigmaUV is huge. Basically M*/t for some t parameter min_t_formation_Myr, by default set to None --- zeus21/UVLFs.py | 23 +++++++++++++++++++++-- zeus21/inputs.py | 4 +++- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/zeus21/UVLFs.py b/zeus21/UVLFs.py index ec08a94..baff27b 100644 --- a/zeus21/UVLFs.py +++ b/zeus21/UVLFs.py @@ -80,9 +80,28 @@ 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 = np.zeros_like(weights_unnormalized) + + 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) + UVLF_filtered = np.trapz(weights.T * HMFcurr, HMF_interpolator.Mhtab, axis=-1) + if(Astro_Parameters.USE_POPIII==False): return UVLF_filtered diff --git a/zeus21/inputs.py b/zeus21/inputs.py index e838903..1fbca4a 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 ): @@ -413,6 +414,7 @@ def __init__(self, UserParams, Cosmo_Parameters, print('Quadratic SFRD not yet implemented when Flag_emulate_21cmfast = True; the code will use quadratic_SFRD_lognormal = False') self.quadratic_SFRD_lognormal = False + 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 From 84df2d1e8d0968c0106f8b3ea3065a6a73471f91 Mon Sep 17 00:00:00 2001 From: Julian Munoz Date: Wed, 29 Apr 2026 12:06:29 -0500 Subject: [PATCH 2/6] Update UVLFs.py comment out old weights --- zeus21/UVLFs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/zeus21/UVLFs.py b/zeus21/UVLFs.py index baff27b..1f706fb 100644 --- a/zeus21/UVLFs.py +++ b/zeus21/UVLFs.py @@ -98,7 +98,7 @@ def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwi 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 + # 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) From 1d37a33feeb175919ec8e1f329937fabb3bdebfa Mon Sep 17 00:00:00 2001 From: Julian Munoz Date: Wed, 29 Apr 2026 12:06:57 -0500 Subject: [PATCH 3/6] Update UVLFs.py and fix trapz to trapezoid... --- zeus21/UVLFs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/zeus21/UVLFs.py b/zeus21/UVLFs.py index 1f706fb..0bd43ef 100644 --- a/zeus21/UVLFs.py +++ b/zeus21/UVLFs.py @@ -100,7 +100,7 @@ def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwi ### 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): From a22b5f8d0629ab9bb5f70f35ff72aae53308a6f9 Mon Sep 17 00:00:00 2001 From: Julian Munoz Date: Wed, 29 Apr 2026 12:16:31 -0500 Subject: [PATCH 4/6] Update zeus21/UVLFs.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- zeus21/UVLFs.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/zeus21/UVLFs.py b/zeus21/UVLFs.py index 0bd43ef..7fd1d22 100644 --- a/zeus21/UVLFs.py +++ b/zeus21/UVLFs.py @@ -93,8 +93,6 @@ def UVLF_binned(Astro_Parameters,Cosmo_Parameters,HMF_interpolator, zcenter, zwi xlo_cut = np.fmax(xlo, x_min) weights_unnormalized = (erf(xhi_cut) - erf(xlo_cut)).T/(2.0 * MUVwidths) - weights = np.zeros_like(weights_unnormalized) - 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: From 33d2d9c70fca936be019d0aaabb1be5e6c1b36dd Mon Sep 17 00:00:00 2001 From: Julian Munoz Date: Wed, 29 Apr 2026 12:17:45 -0500 Subject: [PATCH 5/6] Update zeus21/inputs.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- zeus21/inputs.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/zeus21/inputs.py b/zeus21/inputs.py index 1fbca4a..65184d7 100644 --- a/zeus21/inputs.py +++ b/zeus21/inputs.py @@ -413,7 +413,12 @@ 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 From 3e643acc5e3af851f949429d05bdd07883396eee Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 29 Apr 2026 17:23:25 +0000 Subject: [PATCH 6/6] Add test for min_t_formation_Myr in test_UVLFs.py Agent-Logs-Url: https://github.com/ZeusCosmo/Zeus21/sessions/2c5ae1f8-5c98-4298-b426-afd7b6757fb8 Co-authored-by: JulianBMunoz <22434409+JulianBMunoz@users.noreply.github.com> --- tests/test_UVLFs.py | 65 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 64 insertions(+), 1 deletion(-) 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