diff --git a/src/py21cmfast/src/HaloBox.c b/src/py21cmfast/src/HaloBox.c index c38312bb..1aad5b5d 100644 --- a/src/py21cmfast/src/HaloBox.c +++ b/src/py21cmfast/src/HaloBox.c @@ -75,9 +75,13 @@ void set_halo_properties(double halo_mass, double M_turn_a, double M_turn_m, } // no rng for escape fraction yet - fesc = fmin(consts->fesc_10 * pow(halo_mass / 1e10, consts->alpha_esc), 1); + // add redshift dependence to fesc + fesc = fmin(consts->fesc_10 * pow(halo_mass / 1e10, consts->alpha_esc) * consts->zesc_power_law, + 1); if (astro_options_global->USE_MINI_HALOS) - fesc_mini = fmin(consts->fesc_7 * pow(halo_mass / 1e7, consts->alpha_esc), 1); + fesc_mini = fmin( + consts->fesc_7 * pow(halo_mass / 1e7, consts->alpha_esc) * consts->zesc_power_law_mini, + 1); n_ion_sample = stellar_mass * consts->pop2_ion * fesc + stellar_mass_mini * consts->pop3_ion * fesc_mini; @@ -111,10 +115,14 @@ int get_uhmf_averages(double M_min, double M_max, double M_turn_a, double M_turn double prefactor_sfr = prefactor_stars / consts->t_star / t_h; double prefactor_sfr_mini = prefactor_stars_mini / consts->t_star / t_h; - double prefactor_nion = prefactor_stars * consts->fesc_10 * consts->pop2_ion; - double prefactor_nion_mini = prefactor_stars_mini * consts->fesc_7 * consts->pop3_ion; - double prefactor_wsfr = prefactor_sfr * consts->fesc_10 * consts->pop2_ion; - double prefactor_wsfr_mini = prefactor_sfr_mini * consts->fesc_7 * consts->pop3_ion; + double prefactor_nion = + prefactor_stars * consts->fesc_10 * consts->pop2_ion * consts->zesc_power_law; + double prefactor_nion_mini = + prefactor_stars_mini * consts->fesc_7 * consts->pop3_ion * consts->zesc_power_law_mini; + double prefactor_wsfr = + prefactor_sfr * consts->fesc_10 * consts->pop2_ion * consts->zesc_power_law; + double prefactor_wsfr_mini = + prefactor_sfr_mini * consts->fesc_7 * consts->pop3_ion * consts->zesc_power_law_mini; double mass_intgrl; double intgrl_fesc_weighted, intgrl_stars_only; diff --git a/src/py21cmfast/src/IonisationBox.c b/src/py21cmfast/src/IonisationBox.c index aee309f7..d5f1b4d1 100644 --- a/src/py21cmfast/src/IonisationBox.c +++ b/src/py21cmfast/src/IonisationBox.c @@ -161,8 +161,9 @@ void set_ionbox_constants(double redshift, double prev_redshift, struct IonBoxCo consts->T_re = astro_params_global->T_RE; if (matter_options_global->SOURCE_MODEL > 0) { - consts->ion_eff_factor_gl = sc.pop2_ion * sc.fstar_10 * sc.fesc_10; - consts->ion_eff_factor_mini_gl = sc.pop3_ion * sc.fstar_7 * sc.fesc_7; + consts->ion_eff_factor_gl = sc.pop2_ion * sc.fstar_10 * sc.fesc_10 * sc.zesc_power_law; + consts->ion_eff_factor_mini_gl = + sc.pop3_ion * sc.fstar_7 * sc.fesc_7 * sc.zesc_power_law_mini; } else { consts->ion_eff_factor_gl = astro_params_global->HII_EFF_FACTOR; consts->ion_eff_factor_mini_gl = 0.; diff --git a/src/py21cmfast/src/_inputparams_wrapper.h b/src/py21cmfast/src/_inputparams_wrapper.h index a9a6d93f..2a86cab4 100644 --- a/src/py21cmfast/src/_inputparams_wrapper.h +++ b/src/py21cmfast/src/_inputparams_wrapper.h @@ -97,6 +97,8 @@ typedef struct AstroParams { // Escape Fraction float F_ESC10; float ALPHA_ESC; + float BETA_ESC; + float BETA_ESC_MINI; float F_ESC7_MINI; float T_RE; diff --git a/src/py21cmfast/src/debugging.c b/src/py21cmfast/src/debugging.c index 31e072d4..27f2ca39 100644 --- a/src/py21cmfast/src/debugging.c +++ b/src/py21cmfast/src/debugging.c @@ -126,15 +126,17 @@ void writeAstroParams(AstroParams *p) { " ALPHA_STAR=%8.3f\n" " F_ESC10=%8.3f\n" " ALPHA_ESC=%8.3f\n" + " BETA_ESC=%8.3f\n" + " BETA_ESC_MINI=%8.3f\n" " t_STAR=%8.3f\n" " L_X=%10.3e\n" " NU_X_THRESH=%8.3f\n" " X_RAY_SPEC_INDEX=%8.3f,\n" " UPPER_STELLAR_TURNOVER_MASS=%10.3e\n" " UPPER_STELLAR_TURNOVER_INDEX=%8.3e\n", - p->M_TURN, p->R_BUBBLE_MAX, p->F_STAR10, p->ALPHA_STAR, p->F_ESC10, p->ALPHA_ESC, p->t_STAR, - p->L_X, p->NU_X_THRESH, p->X_RAY_SPEC_INDEX, p->UPPER_STELLAR_TURNOVER_MASS, - p->UPPER_STELLAR_TURNOVER_INDEX); + p->M_TURN, p->R_BUBBLE_MAX, p->F_STAR10, p->ALPHA_STAR, p->F_ESC10, p->ALPHA_ESC, + p->BETA_ESC, p->BETA_ESC_MINI, p->t_STAR, p->L_X, p->NU_X_THRESH, p->X_RAY_SPEC_INDEX, + p->UPPER_STELLAR_TURNOVER_MASS, p->UPPER_STELLAR_TURNOVER_INDEX); LOG_INFO( "\n HaloCatalog AstroParams:\n" " SIGMA_STAR=%8.3f\n" diff --git a/src/py21cmfast/src/hmf.c b/src/py21cmfast/src/hmf.c index 4d8391aa..70e7f832 100644 --- a/src/py21cmfast/src/hmf.c +++ b/src/py21cmfast/src/hmf.c @@ -1038,7 +1038,6 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double lnM_co .delta = delta2, .gamma_type = -3, }; - if (lnM1 >= lnM_cond) return 0.; // return 1 halo at the condition mass if delta is exceeded and the condition is within the // integral limits diff --git a/src/py21cmfast/src/scaling_relations.c b/src/py21cmfast/src/scaling_relations.c index 93f5a4d3..380aff38 100644 --- a/src/py21cmfast/src/scaling_relations.c +++ b/src/py21cmfast/src/scaling_relations.c @@ -23,7 +23,8 @@ void print_sc_consts(ScalingConstants *c) { LOG_DEBUG("SHMR: f10 %.2e a %.2e f7 %.2e a_mini %.2e sigma %.2e", c->fstar_10, c->alpha_star, c->fstar_7, c->alpha_star_mini, c->sigma_star); LOG_DEBUG("Upper: a_upper %.2e pivot %.2e", c->alpha_upper, c->pivot_upper); - LOG_DEBUG("FESC: f10 %.2e a %.2e f7 %.2e", c->fesc_10, c->alpha_esc, c->fesc_7); + LOG_DEBUG("FESC: f10 %.2e a %.2e b %.2e f7 %.2e zesc_pl %.2e", c->fesc_10, c->alpha_esc, + c->beta_esc, c->beta_esc_mini, c->fesc_7, c->zesc_power_law, c->zesc_power_law_mini); LOG_DEBUG("SSFR: t* %.2e th %.8e sigma %.2e idx %.2e", c->t_star, c->t_h, c->sigma_sfr_lim, c->sigma_sfr_idx); LOG_DEBUG("Turnovers (nofb) ACG %.2e MCG %.2e Upper %.2e", c->mturn_a_nofb, c->mturn_m_nofb, @@ -63,9 +64,15 @@ void set_scaling_constants(double redshift, ScalingConstants *consts, bool use_p consts->sigma_xray = astro_params_global->SIGMA_LX; consts->alpha_esc = astro_params_global->ALPHA_ESC; + consts->beta_esc = astro_params_global->BETA_ESC; + consts->beta_esc_mini = astro_params_global->BETA_ESC_MINI; consts->fesc_10 = astro_params_global->F_ESC10; consts->fesc_7 = astro_params_global->F_ESC7_MINI; + // Calculate the redshift-dependent power-law term for escape fraction once per redshift + consts->zesc_power_law = pow((1 + redshift) / 8.0, consts->beta_esc); + consts->zesc_power_law_mini = pow((1 + redshift) / 8.0, consts->beta_esc_mini); + if (use_photoncons) { if (astro_options_global->PHOTON_CONS_TYPE == 2) consts->alpha_esc = get_fesc_fit(redshift); @@ -108,6 +115,9 @@ ScalingConstants evolve_scaling_constants_sfr(ScalingConstants *sc) { sc_sfrd.fesc_10 = 1.; sc_sfrd.fesc_7 = 1.; sc_sfrd.alpha_esc = 0.; + sc_sfrd.beta_esc = 0.; + sc_sfrd.beta_esc_mini = 0.; + sc_sfrd.zesc_power_law = 1.; sc_sfrd.Mlim_Fesc = 0.; sc_sfrd.Mlim_Fesc_mini = 0.; @@ -149,6 +159,10 @@ ScalingConstants evolve_scaling_constants_to_redshift(double redshift, ScalingCo sc_z.mturn_m_nofb = lyman_werner_threshold(redshift, 0., sc_z.vcb_norel); } + // Update the redshift-dependent power-law term for escape fraction + sc_z.zesc_power_law = pow((1 + redshift) / 8.0, sc_z.beta_esc); + sc_z.zesc_power_law_mini = pow((1 + redshift) / 8.0, sc_z.beta_esc_mini); + return sc_z; } diff --git a/src/py21cmfast/src/scaling_relations.h b/src/py21cmfast/src/scaling_relations.h index 0743ae98..7eddb029 100644 --- a/src/py21cmfast/src/scaling_relations.h +++ b/src/py21cmfast/src/scaling_relations.h @@ -36,10 +36,15 @@ typedef struct ScalingConstants { double fesc_10; double alpha_esc; + double beta_esc; double fesc_7; + double beta_esc_mini; double pop2_ion; double pop3_ion; + double zesc_power_law; + double zesc_power_law_mini; + double vcb_norel; double acg_thresh; double mturn_a_nofb; diff --git a/src/py21cmfast/wrapper/inputs.py b/src/py21cmfast/wrapper/inputs.py index cb2780e3..a010b8a1 100644 --- a/src/py21cmfast/wrapper/inputs.py +++ b/src/py21cmfast/wrapper/inputs.py @@ -1103,11 +1103,11 @@ def _USE_MINI_HALOS_vald(self, att, val): @PHOTON_CONS_TYPE.validator def _PHOTON_CONS_TYPE_vld(self, att, val): - """Raise an error if using PHOTON_CONS_TYPE='z_photoncons' and USE_MINI_HALOS is True.""" + """Raise an error if using PHOTON_CONS_TYPE='z-photoncons' and USE_MINI_HALOS is True.""" if self.USE_MINI_HALOS and val == "z-photoncons": raise ValueError( "USE_MINI_HALOS is not compatible with the redshift-based" - " photon conservation corrections (PHOTON_CONS_TYPE=='z_photoncons')! " + " photon conservation corrections (PHOTON_CONS_TYPE=='z-photoncons')! " ) @USE_EXP_FILTER.validator @@ -1182,6 +1182,12 @@ class AstroParams(InputStruct): ALPHA_ESC : float, optional Power-law index of escape fraction as a function of halo mass. See Sec 2.1 of Park+2018. + BETA_ESC : float, optional + Power-law index of escape fraction as a function of redshift. See Eq. 2 of + Qin+2025. + BETA_ESC_MINI : float, optional + Power-law index of escape fraction as a function of redshift for minihalos. + If the scaling relations are not provided explicitly, we extend the ACG ones by default. M_TURN : float, optional Turnover mass (in log10 solar mass units) for quenching of star formation in halos, due to SNe or photo-heating feedback, or inefficient gas accretion. @@ -1276,6 +1282,11 @@ class AstroParams(InputStruct): default=-0.5, converter=float, ) + BETA_ESC: float = field( + default=0.0, + converter=float, + ) + BETA_ESC_MINI: float = field(converter=float) F_ESC7_MINI: float = field( default=-2.0, converter=float, @@ -1372,6 +1383,10 @@ def _F_STAR7_MINI_default(self): def _ALPHA_STAR_MINI_default(self): return self.ALPHA_STAR + @BETA_ESC_MINI.default + def _BETA_ESC_MINI_default(self): + return self.BETA_ESC + @L_X_MINI.default def _L_X_MINI_default(self): return self.L_X @@ -1430,6 +1445,8 @@ def get_logspaced_redshifts( ) -> tuple[float]: """Compute a sequence of redshifts to evolve over that are log-spaced.""" redshifts = [min_redshift] + if z_step_factor <= 1.0: + return (min_redshift,) while redshifts[-1] < max_redshift: redshifts.append((redshifts[-1] + 1.0) * z_step_factor - 1.0) @@ -1651,6 +1668,18 @@ def _astro_params_validator(self, att, val): "update of M_TURN", stacklevel=2, ) + if ( + self.astro_options.USE_MASS_DEPENDENT_ZETA + and val.BETA_ESC != 0 + and self.astro_options.PHOTON_CONS_TYPE + not in ["no-photoncons", "alpha-photoncons"] + ): + warnings.warn( + f"You have set BETA_ESC != 0 but PHOTON_CONS_TYPE is {self.astro_options.PHOTON_CONS_TYPE}. " + "This changes the escape fraction so it is not consistent with the manual setting of scaling." + "Set PHOTON_CONS_TYPE to 'no-photoncons' or 'alpha-photoncons' if you want the scaling to be exact.", + stacklevel=2, + ) if ( self.astro_options.HII_FILTER == "sharp-k" diff --git a/tests/test_input_structs.py b/tests/test_input_structs.py index 3f94d298..c00a9486 100644 --- a/tests/test_input_structs.py +++ b/tests/test_input_structs.py @@ -368,6 +368,22 @@ def test_dim_setting_serialization( new = new.evolve_input_structs(**evolved) self.check_attributes_dim(new.simulation_options, expected) + with pytest.warns( + UserWarning, + match="You have set BETA_ESC != 0 but PHOTON_CONS_TYPE is f-photoncons.", + ): + InputParameters( + cosmo_params=CosmoParams(), + astro_params=AstroParams(BETA_ESC=0.1), + simulation_options=SimulationOptions(), + matter_options=MatterOptions(), + astro_options=AstroOptions( + PHOTON_CONS_TYPE="f-photoncons", + USE_MASS_DEPENDENT_ZETA=True, + ), + random_seed=1, + ) + class TestMatterOptions: """Tests of the MatterOptions class."""