diff --git a/src/py21cmfast/outputs.py b/src/py21cmfast/outputs.py index 621ffa34..8ccb933e 100644 --- a/src/py21cmfast/outputs.py +++ b/src/py21cmfast/outputs.py @@ -544,6 +544,8 @@ def _get_box_structures(self) -> dict[str, dict | tuple[int]]: "Nion_box": shape, "z_re_box": shape, "Nrec_box": shape, + "nHI_box": shape, + "C_box": shape, "temp_kinetic_all_gas": shape, "Fcoll": filter_shape, } @@ -581,6 +583,8 @@ def get_required_input_arrays(self, input_box: _BaseOutputStruct) -> list[str]: if self.flag_options.INHOMO_RECO: required += [ "Nrec_box", + "nHI_box", + "C_box", ] if ( self.flag_options.USE_MASS_DEPENDENT_ZETA @@ -1119,6 +1123,7 @@ def __init__( random_seed, lightcones, node_redshifts=None, + node_redshifts_adjusted=None, mean_f_colls=None, mean_f_coll_MINIs=None, global_quantities=None, @@ -1135,6 +1140,7 @@ def __init__( self.astro_params = astro_params self.flag_options = flag_options self.node_redshifts = node_redshifts + self.node_redshifts_adjusted = node_redshifts_adjusted self.cache_files = cache_files self.log10_mturnovers = log10_mturnovers self.log10_mturnovers_mini = log10_mturnovers_mini @@ -1146,7 +1152,7 @@ def __init__( mean_f_colls * 10**astro_params.F_STAR10 * 10**astro_params.F_ESC10 - * ((1.0 + self.node_redshifts) / 8) ** astro_params.BETA_ESC + * ((1.0 + self.node_redshifts_adjusted) / 8) ** astro_params.BETA_ESC * self.global_params["Pop2_ion"] ) self.Nion_mcg = ( @@ -1249,6 +1255,9 @@ def _write_particulars(self, fname): global_q[k] = v f["node_redshifts"] = self.node_redshifts + f["node_redshifts_adjusted"] = self.node_redshifts_adjusted + f["Nion_acg"] = self.Nion_acg + f["Nion_mcg"] = self.Nion_mcg @classmethod def _read_inputs(cls, fname): diff --git a/src/py21cmfast/src/21cmFAST.h b/src/py21cmfast/src/21cmFAST.h index a29200fb..a0c1bb5a 100644 --- a/src/py21cmfast/src/21cmFAST.h +++ b/src/py21cmfast/src/21cmFAST.h @@ -126,16 +126,17 @@ struct TsBox{ }; struct IonizedBox{ + double redshift_PC; double mean_f_coll; double mean_f_coll_MINI; - double mean_f_coll_PC; - double mean_f_coll_MINI_PC; double log10_Mturnover_ave; double log10_Mturnover_MINI_ave; float *xH_box; float *Gamma12_box; float *MFP_box; float *z_re_box; + float *C_box; + float *nHI_box; float *Nrec_box; float *Nion_box; float *temp_kinetic_all_gas; diff --git a/src/py21cmfast/src/IonisationBox.c b/src/py21cmfast/src/IonisationBox.c index 42a6b5be..24f6a03f 100644 --- a/src/py21cmfast/src/IonisationBox.c +++ b/src/py21cmfast/src/IonisationBox.c @@ -117,7 +117,7 @@ LOG_SUPER_DEBUG("defined parameters"); pixel_volume = pow(user_params->BOX_LEN/((float)(user_params->HII_DIM)), 3); - double F_ESC10_zterm, prev_F_ESC10_zterm; + double F_ESC10_zterm, prev_F_ESC10_zterm; // For recombinations if(flag_options->INHOMO_RECO) { @@ -181,6 +181,7 @@ LOG_DEBUG("original prev_redshift=%f, updated prev_redshift=%f delta-z = %f", st // Throw(ParameterError); Throw(PhotonConsError); } + box->redshift_PC = redshift; } if(flag_options->USE_MASS_DEPENDENT_ZETA) { @@ -192,7 +193,7 @@ LOG_DEBUG("original prev_redshift=%f, updated prev_redshift=%f delta-z = %f", st if(flag_options->PHOTON_CONS) { stored_F_ESC10_zterm = pow((1.+stored_redshift)/8., astro_params->BETA_ESC); stored_prev_F_ESC10_zterm = pow((1.+stored_prev_redshift)/8., astro_params->BETA_ESC); - } + } } else { ION_EFF_FACTOR = astro_params->HII_EFF_FACTOR; @@ -380,10 +381,10 @@ LOG_DEBUG("first redshift, do some initialization"); R_BUBBLE_MAX = 25.483241248322766 / cosmo_params->hlittle; else R_BUBBLE_MAX = 112 / cosmo_params->hlittle * pow( (1.+redshift) / 5. , -4.4); - if (flag_options->USE_MINI_HALOS){ - R_BUBBLE_MAX_MAX = astro_params->R_BUBBLE_MAX; - LOG_DEBUG("set max radius at this redshift as %f vs all %f", R_BUBBLE_MAX, R_BUBBLE_MAX_MAX); - } + if (flag_options->USE_MINI_HALOS){ + R_BUBBLE_MAX_MAX = astro_params->R_BUBBLE_MAX; + LOG_DEBUG("set max radius at this redshift as %f vs all %f", R_BUBBLE_MAX, R_BUBBLE_MAX_MAX); + } } else R_BUBBLE_MAX = astro_params->R_BUBBLE_MAX; @@ -399,10 +400,10 @@ LOG_DEBUG("first redshift, do some initialization"); // really painful to get the length... counter = 1; R=fmax(global_params.R_BUBBLE_MIN, (cell_length_factor*user_params->BOX_LEN/(float)user_params->HII_DIM)); - if (flag_options->EVOLVING_R_BUBBLE_MAX) - R_BUBBLE_MAX_tmp = R_BUBBLE_MAX_MAX; - else - R_BUBBLE_MAX_tmp = R_BUBBLE_MAX; + if (flag_options->EVOLVING_R_BUBBLE_MAX) + R_BUBBLE_MAX_tmp = R_BUBBLE_MAX_MAX; + else + R_BUBBLE_MAX_tmp = R_BUBBLE_MAX; while ((R - fmin(R_BUBBLE_MAX_tmp, L_FACTOR*user_params->BOX_LEN)) <= FRACT_FLOAT_ERR ){ if(R >= fmin(R_BUBBLE_MAX_tmp, L_FACTOR*user_params->BOX_LEN)) { stored_R = R/(global_params.DELTA_R_HII_FACTOR); @@ -415,10 +416,6 @@ LOG_DEBUG("first redshift, do some initialization"); previous_ionize_box->Fcoll_MINI = (float *) calloc(HII_TOT_NUM_PIXELS*counter, sizeof(float)); previous_ionize_box->mean_f_coll = 0.0; previous_ionize_box->mean_f_coll_MINI = 0.0; - if(flag_options->PHOTON_CONS) { - previous_ionize_box->mean_f_coll_PC = 0.0; - previous_ionize_box->mean_f_coll_MINI_PC = 0.0; - } #pragma omp parallel shared(prev_deltax_unfiltered) private(i,j,k) num_threads(user_params->N_THREADS) { @@ -524,10 +521,10 @@ LOG_SUPER_DEBUG("average turnover masses are %.2f and %.2f for ACGs and MCGs", b Mlim_Fstar = Mass_limit_bisection(M_MIN, 1e16, astro_params->ALPHA_STAR, astro_params->F_STAR10); Mlim_Fesc = Mass_limit_bisection(M_MIN, 1e16, astro_params->ALPHA_ESC, astro_params->F_ESC10* F_ESC10_zterm); prev_Mlim_Fesc = Mass_limit_bisection(M_MIN, 1e16, astro_params->ALPHA_ESC, astro_params->F_ESC10* prev_F_ESC10_zterm); - if(flag_options->PHOTON_CONS){ + if(flag_options->PHOTON_CONS){ stored_Mlim_Fesc = Mass_limit_bisection(M_MIN, 1e16, astro_params->ALPHA_ESC, astro_params->F_ESC10* stored_F_ESC10_zterm); stored_prev_Mlim_Fesc = Mass_limit_bisection(M_MIN, 1e16, astro_params->ALPHA_ESC, astro_params->F_ESC10* stored_prev_F_ESC10_zterm); - } + } } else { @@ -605,39 +602,20 @@ LOG_SUPER_DEBUG("sigma table has been initialised"); // Determine the normalisation for the excursion set algorithm if (flag_options->USE_MASS_DEPENDENT_ZETA) { if (flag_options->USE_MINI_HALOS){ - if (previous_ionize_box->mean_f_coll * prev_ION_EFF_FACTOR < 1e-4){ + if (previous_ionize_box->mean_f_coll * prev_ION_EFF_FACTOR < 1e-4) box->mean_f_coll = Nion_General(redshift,M_MIN,Mturnover,astro_params->ALPHA_STAR,astro_params->ALPHA_ESC, astro_params->F_STAR10,astro_params->F_ESC10* F_ESC10_zterm,Mlim_Fstar,Mlim_Fesc); - if(flag_options->PHOTON_CONS) { - box->mean_f_coll_PC = Nion_General(stored_redshift,M_MIN,Mturnover,astro_params->ALPHA_STAR,astro_params->ALPHA_ESC, - astro_params->F_STAR10,astro_params->F_ESC10* stored_F_ESC10_zterm,Mlim_Fstar,stored_Mlim_Fesc); - } - } - else{ + else box->mean_f_coll = previous_ionize_box->mean_f_coll + \ Nion_General(redshift,M_MIN,Mturnover,astro_params->ALPHA_STAR,astro_params->ALPHA_ESC, astro_params->F_STAR10,astro_params->F_ESC10 * F_ESC10_zterm,Mlim_Fstar,Mlim_Fesc) - \ Nion_General(prev_redshift,M_MIN,Mturnover,astro_params->ALPHA_STAR,astro_params->ALPHA_ESC, astro_params->F_STAR10,astro_params->F_ESC10 * F_ESC10_zterm,Mlim_Fstar,Mlim_Fesc);// similar to Mturnover, F_ESC10_zterm should not be prev_F_ESC10_zterm - if(flag_options->PHOTON_CONS) { - box->mean_f_coll_PC = previous_ionize_box->mean_f_coll_PC + \ - Nion_General(stored_redshift,M_MIN,Mturnover,astro_params->ALPHA_STAR,astro_params->ALPHA_ESC, - astro_params->F_STAR10,astro_params->F_ESC10 *stored_F_ESC10_zterm,Mlim_Fstar,stored_Mlim_Fesc) - \ - Nion_General(stored_prev_redshift,M_MIN,Mturnover,astro_params->ALPHA_STAR,astro_params->ALPHA_ESC, - astro_params->F_STAR10,astro_params->F_ESC10 *stored_F_ESC10_zterm,Mlim_Fstar,stored_Mlim_Fesc); - } - } - if (previous_ionize_box->mean_f_coll_MINI * ION_EFF_FACTOR_MINI < 1e-4){ + if (previous_ionize_box->mean_f_coll_MINI * ION_EFF_FACTOR_MINI < 1e-4) box->mean_f_coll_MINI = Nion_General_MINI(redshift,M_MIN,Mturnover_MINI,Mcrit_atom, astro_params->ALPHA_STAR_MINI,astro_params->ALPHA_ESC,astro_params->F_STAR7_MINI, astro_params->F_ESC7_MINI,Mlim_Fstar_MINI,Mlim_Fesc_MINI); - if(flag_options->PHOTON_CONS) { - box->mean_f_coll_MINI_PC = Nion_General_MINI(stored_redshift,M_MIN,Mturnover_MINI,Mcrit_atom, - astro_params->ALPHA_STAR_MINI,astro_params->ALPHA_ESC,astro_params->F_STAR7_MINI, - astro_params->F_ESC7_MINI,Mlim_Fstar_MINI,Mlim_Fesc_MINI); - } - } - else{ + else box->mean_f_coll_MINI = previous_ionize_box->mean_f_coll_MINI + \ Nion_General_MINI(redshift,M_MIN,Mturnover_MINI,Mcrit_atom,astro_params->ALPHA_STAR_MINI, astro_params->ALPHA_ESC,astro_params->F_STAR7_MINI,astro_params->F_ESC7_MINI @@ -645,16 +623,6 @@ LOG_SUPER_DEBUG("sigma table has been initialised"); Nion_General_MINI(prev_redshift,M_MIN,Mturnover_MINI,Mcrit_atom,astro_params->ALPHA_STAR_MINI, astro_params->ALPHA_ESC,astro_params->F_STAR7_MINI,astro_params->F_ESC7_MINI, Mlim_Fstar_MINI,Mlim_Fesc_MINI); - if(flag_options->PHOTON_CONS) { - box->mean_f_coll_MINI_PC = previous_ionize_box->mean_f_coll_MINI_PC + \ - Nion_General_MINI(stored_redshift,M_MIN,Mturnover_MINI,Mcrit_atom,astro_params->ALPHA_STAR_MINI, - astro_params->ALPHA_ESC,astro_params->F_STAR7_MINI,astro_params->F_ESC7_MINI - ,Mlim_Fstar_MINI,Mlim_Fesc_MINI) - \ - Nion_General_MINI(stored_prev_redshift,M_MIN,Mturnover_MINI,Mcrit_atom,astro_params->ALPHA_STAR_MINI, - astro_params->ALPHA_ESC,astro_params->F_STAR7_MINI,astro_params->F_ESC7_MINI, - Mlim_Fstar_MINI,Mlim_Fesc_MINI); - } - } f_coll_min = Nion_General(global_params.Z_HEAT_MAX,M_MIN,Mturnover,astro_params->ALPHA_STAR, astro_params->ALPHA_ESC,astro_params->F_STAR10,astro_params->F_ESC10*F_ESC10_zterm,Mlim_Fstar,Mlim_Fesc); f_coll_min_MINI = Nion_General_MINI(global_params.Z_HEAT_MAX,M_MIN,Mturnover_MINI,Mcrit_atom, @@ -665,20 +633,12 @@ LOG_SUPER_DEBUG("sigma table has been initialised"); box->mean_f_coll = Nion_General(redshift,M_MIN,Mturnover,astro_params->ALPHA_STAR,astro_params->ALPHA_ESC, astro_params->F_STAR10,astro_params->F_ESC10*F_ESC10_zterm,Mlim_Fstar,Mlim_Fesc); box->mean_f_coll_MINI = 0.; - if(flag_options->PHOTON_CONS) { - box->mean_f_coll_PC = Nion_General(stored_redshift,M_MIN,Mturnover,astro_params->ALPHA_STAR,astro_params->ALPHA_ESC, - astro_params->F_STAR10,astro_params->F_ESC10*stored_F_ESC10_zterm,Mlim_Fstar,stored_Mlim_Fesc); - box->mean_f_coll_MINI_PC = 0.; - } f_coll_min = Nion_General(global_params.Z_HEAT_MAX,M_MIN,Mturnover,astro_params->ALPHA_STAR,astro_params->ALPHA_ESC, astro_params->F_STAR10,astro_params->F_ESC10*F_ESC10_zterm,Mlim_Fstar,Mlim_Fesc); } } else { box->mean_f_coll = FgtrM_General(redshift, M_MIN); - if(flag_options->PHOTON_CONS) { - box->mean_f_coll_PC = FgtrM_General(stored_redshift, M_MIN); - } } if(isfinite(box->mean_f_coll)==0) { @@ -819,10 +779,10 @@ LOG_SUPER_DEBUG("excursion set normalisation, mean_f_coll_MINI: %e", box->mean_f // loop through the filter radii (in Mpc) R=fmax(global_params.R_BUBBLE_MIN, (cell_length_factor*user_params->BOX_LEN/(float)user_params->HII_DIM)); - if ((flag_options->EVOLVING_R_BUBBLE_MAX) && (flag_options->USE_MINI_HALOS)) - R_BUBBLE_MAX_tmp = R_BUBBLE_MAX_MAX; - else - R_BUBBLE_MAX_tmp = R_BUBBLE_MAX; + if ((flag_options->EVOLVING_R_BUBBLE_MAX) && (flag_options->USE_MINI_HALOS)) + R_BUBBLE_MAX_tmp = R_BUBBLE_MAX_MAX; + else + R_BUBBLE_MAX_tmp = R_BUBBLE_MAX; while ((R - fmin(R_BUBBLE_MAX_tmp, L_FACTOR * user_params->BOX_LEN)) <= FRACT_FLOAT_ERR) { R *= global_params.DELTA_R_HII_FACTOR; @@ -1414,10 +1374,10 @@ LOG_ULTRA_DEBUG("while loop for until RtoM(R)=%f reaches M_MIN=%f", RtoM(R), M_M // check if fully ionized! if ( (f_coll * ION_EFF_FACTOR + f_coll_MINI * ION_EFF_FACTOR_MINI> (1. - xHII_from_xrays)*(1.0+rec)) ){ //IONIZED!! - // This is for the evolving Rmax case in minihalo runs - // We need to make sure the filtering radii are consitent across all snapshots - // But we do not need, at higher redshift, to assign reionization for radius larger than the real Rmax - if (~((flag_options->EVOLVING_R_BUBBLE_MAX) && (flag_options->USE_MINI_HALOS) && (R > R_BUBBLE_MAX))){ + // This is for the evolving Rmax case in minihalo runs + // We need to make sure the filtering radii are consitent across all snapshots + // But we do not need, at higher redshift, to assign reionization for radius larger than the real Rmax + if (~((flag_options->EVOLVING_R_BUBBLE_MAX) && (flag_options->USE_MINI_HALOS) && (R > R_BUBBLE_MAX))){ // if this is the first crossing of the ionization barrier for this cell (largest R), record the gamma // this assumes photon-starved growth of HII regions... breaks down post EoR if (flag_options->INHOMO_RECO && (box->xH_box[HII_R_INDEX(x,y,z)] > FRACT_FLOAT_ERR) ){ @@ -1439,7 +1399,7 @@ LOG_ULTRA_DEBUG("while loop for until RtoM(R)=%f reaches M_MIN=%f", RtoM(R), M_M if (global_params.FIND_BUBBLE_ALGORITHM == 1) // sphere method update_in_sphere(box->xH_box, user_params->HII_DIM, HII_D_PARA, R/(user_params->BOX_LEN), \ x/(user_params->HII_DIM+0.0), y/(user_params->HII_DIM+0.0), z/(HII_D_PARA+0.0)); - } + } } // end ionized // If not fully ionized, then assign partial ionizations else if (LAST_FILTER_STEP && (box->xH_box[HII_R_INDEX(x, y, z)] > TINY)) { @@ -1590,6 +1550,8 @@ LOG_ULTRA_DEBUG("while loop for until RtoM(R)=%f reaches M_MIN=%f", RtoM(R), M_M dNrec = splined_recombination_rate(z_eff - 1., box->Gamma12_box[HII_R_INDEX(x, y, z)]) * fabs_dtdz * ZSTEP * (1. - box->xH_box[HII_R_INDEX(x, y, z)]); + box->C_box[HII_R_INDEX(x, y, z)] = splined_clumping_factor(z_eff - 1., box->Gamma12_box[HII_R_INDEX(x, y, z)]); + box->nHI_box[HII_R_INDEX(x, y, z)] = splined_residual_neutral_hydrogen(z_eff - 1., box->Gamma12_box[HII_R_INDEX(x, y, z)]); if (isfinite(dNrec) == 0) { something_finite_or_infinite = 1; diff --git a/src/py21cmfast/src/recombinations.c b/src/py21cmfast/src/recombinations.c index 198356ff..bf42db71 100755 --- a/src/py21cmfast/src/recombinations.c +++ b/src/py21cmfast/src/recombinations.c @@ -20,18 +20,25 @@ static gsl_spline *beta_spline; #define RR_lnGamma_min (double) (-10) // min ln gamma12 used #define RR_DEL_lnGamma (float) (0.1) static double RR_table[RR_Z_NPTS][RR_lnGamma_NPTS], lnGamma_values[RR_lnGamma_NPTS]; -static gsl_interp_accel *RR_acc[RR_Z_NPTS]; -static gsl_spline *RR_spline[RR_Z_NPTS]; +static double CF_table[RR_Z_NPTS][RR_lnGamma_NPTS], RNH_table[RR_Z_NPTS][RR_lnGamma_NPTS]; +static gsl_interp_accel *RR_acc[RR_Z_NPTS], *CF_acc[RR_Z_NPTS], *RNH_acc[RR_Z_NPTS]; +static gsl_spline *RR_spline[RR_Z_NPTS], *CF_spline[RR_Z_NPTS], *RNH_spline[RR_Z_NPTS]; /*** FUNCTION PROTOTYPES ***/ double splined_recombination_rate(double z_eff, double gamma12_bg); // assumes T=1e4 and case B +double splined_clumping_factor(double z_eff, double gamma12_bg); // assumes T=1e4 and case B +double splined_residual_neutral_hydrogen(double z_eff, double gamma12_bg); // assumes T=1e4 and case B double recombination_rate(double z_eff, double gamma12_bg, double T4, int usecaseB); +double clumping_factor(double z_eff, double gamma12_bg, double T4, int usecaseB); +double residual_neutral_hydrogen(double z_eff, double gamma12_bg, double T4, int usecaseB); void init_MHR(); /*initializes the lookup table for the PDF density integral in MHR00 model at redshift z*/ void free_MHR(); /* deallocates the gsl structures from init_MHR */ double Gamma_SS(double Gamma_bg, double Delta, double T_4, double z);//ionization rate w. self shielding double MHR_rr (double del, void *params); +double MHR_cf (double del, void *params); +double MHR_rnh (double del, void *params); double A_MHR(double z); /*returns the A parameter in MHR00model*/ double C_MHR(double z); /*returns the C parameter in MHR00model*/ double beta_MHR(double z); /*returns the beta parameter in MHR00model*/ @@ -69,6 +76,52 @@ double splined_recombination_rate(double z_eff, double gamma12_bg){ return gsl_spline_eval(RR_spline[z_ct], lnGamma, RR_acc[z_ct]); } +double splined_clumping_factor(double z_eff, double gamma12_bg){ + int z_ct = (int) (z_eff / RR_DEL_Z + 0.5); // round to nearest int + double lnGamma = log(gamma12_bg); + + // check out of bounds + if ( z_ct < 0 ){ // out of array bounds + z_ct = 0; + } + else if (z_ct >= RR_Z_NPTS){ + z_ct = RR_Z_NPTS-1; + } + + if (lnGamma < RR_lnGamma_min){ + return 0; + } + else if (lnGamma >= (RR_lnGamma_min + RR_DEL_lnGamma * ( RR_lnGamma_NPTS - 1 )) ){ + LOG_WARNING("splined_clumping_factor: Gamma12 of %g is outside of interpolation array", gamma12_bg); + lnGamma = RR_lnGamma_min + RR_DEL_lnGamma * ( RR_lnGamma_NPTS - 1 ) - FRACT_FLOAT_ERR; + } + + return gsl_spline_eval(CF_spline[z_ct], lnGamma, CF_acc[z_ct]); +} + +double splined_residual_neutral_hydrogen(double z_eff, double gamma12_bg){ + int z_ct = (int) (z_eff / RR_DEL_Z + 0.5); // round to nearest int + double lnGamma = log(gamma12_bg); + + // check out of bounds + if ( z_ct < 0 ){ // out of array bounds + z_ct = 0; + } + else if (z_ct >= RR_Z_NPTS){ + z_ct = RR_Z_NPTS-1; + } + + if (lnGamma < RR_lnGamma_min){ + return 0; + } + else if (lnGamma >= (RR_lnGamma_min + RR_DEL_lnGamma * ( RR_lnGamma_NPTS - 1 )) ){ + LOG_WARNING("splined_residual_neutral_hydrogen: Gamma12 of %g is outside of interpolation array", gamma12_bg); + lnGamma = RR_lnGamma_min + RR_DEL_lnGamma * ( RR_lnGamma_NPTS - 1 ) - FRACT_FLOAT_ERR; + } + + return gsl_spline_eval(RNH_spline[z_ct], lnGamma, RNH_acc[z_ct]); +} + void init_MHR(){ int z_ct, gamma_ct; float z, gamma; @@ -88,12 +141,20 @@ void init_MHR(){ lnGamma_values[gamma_ct] = RR_lnGamma_min + gamma_ct*RR_DEL_lnGamma; // ln of Gamma12 gamma = exp(lnGamma_values[gamma_ct]); RR_table[z_ct][gamma_ct] = recombination_rate(z, gamma, 1, 1); // CHANGE THIS TO INCLUDE TEMPERATURE + CF_table[z_ct][gamma_ct] = clumping_factor(z, gamma, 1, 1); // CHANGE THIS TO INCLUDE TEMPERATURE + RNH_table[z_ct][gamma_ct] = residual_neutral_hydrogen(z, gamma, 1, 1); // CHANGE THIS TO INCLUDE TEMPERATURE } // set up the spline in gamma RR_acc[z_ct] = gsl_interp_accel_alloc(); + CF_acc[z_ct] = gsl_interp_accel_alloc(); + RNH_acc[z_ct] = gsl_interp_accel_alloc(); RR_spline[z_ct] = gsl_spline_alloc (gsl_interp_cspline, RR_lnGamma_NPTS); gsl_spline_init(RR_spline[z_ct], lnGamma_values, RR_table[z_ct], RR_lnGamma_NPTS); + CF_spline[z_ct] = gsl_spline_alloc (gsl_interp_cspline, RR_lnGamma_NPTS); + gsl_spline_init(CF_spline[z_ct], lnGamma_values, CF_table[z_ct], RR_lnGamma_NPTS); + RNH_spline[z_ct] = gsl_spline_alloc (gsl_interp_cspline, RR_lnGamma_NPTS); + gsl_spline_init(RNH_spline[z_ct], lnGamma_values, RNH_table[z_ct], RR_lnGamma_NPTS); } // go to next redshift @@ -111,6 +172,10 @@ void free_MHR(){ for (z_ct=0; z_ct < RR_Z_NPTS; z_ct++){ gsl_spline_free (RR_spline[z_ct]); gsl_interp_accel_free(RR_acc[z_ct]); + gsl_spline_free (CF_spline[z_ct]); + gsl_interp_accel_free(CF_acc[z_ct]); + gsl_spline_free (RNH_spline[z_ct]); + gsl_interp_accel_free(RNH_acc[z_ct]); } return; @@ -147,6 +212,33 @@ double MHR_rr (double lnD, void *params){ return 1e15*n_H * PDelta * alpha * x_e * x_e * del * del;//note extra D since we are integrating over lnD } +double MHR_cf (double lnD, void *params){ + double del=exp(lnD); + RR_par p = *(RR_par *) params; + double z = p.z; + double gamma = Gamma_SS(p.gamma12_bg, del, p.T4, z); + double n_H = p.avenH*del; + double x_e = 1.0 - neutral_fraction(n_H, p.T4, gamma, p.usecaseB); + double PDelta; + + PDelta = p.A * exp( - 0.5*pow((pow(del,-2.0/3.0)- p.C_0 ) / ((2.0*7.61/(3.0*(1.0+z)))), 2)) * pow(del, p.beta); + + return del * PDelta * x_e * x_e * del * del;//note extra D since we are integrating over lnD +} + +double MHR_rnh (double lnD, void *params){ + double del=exp(lnD); + RR_par p = *(RR_par *) params; + double z = p.z; + double gamma = Gamma_SS(p.gamma12_bg, del, p.T4, z); + double n_H = p.avenH*del; + double x_HI = neutral_fraction(n_H, p.T4, gamma, p.usecaseB); + double PDelta; + + PDelta = p.A * exp( - 0.5*pow((pow(del,-2.0/3.0)- p.C_0 ) / ((2.0*7.61/(3.0*(1.0+z)))), 2)) * pow(del, p.beta); + + return 1e6 * del * PDelta * x_HI * del;//note extra D since we are integrating over lnD +} // returns the recombination rate per baryon (1/1e15s), integrated over the MHR density PDF, // given an ionizing background of gamma12_bg @@ -184,6 +276,70 @@ double recombination_rate(double z, double gamma12_bg, double T4, int usecaseB){ return result; } +double clumping_factor(double z, double gamma12_bg, double T4, int usecaseB){ + double result, error, lower_limit, upper_limit, A, C_0, beta, avenH; + gsl_function F; + double rel_tol = 0.01; //<- relative tolerance + gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000); + RR_par p = {z, gamma12_bg, T4, A_MHR(z), C_MHR(z), beta_MHR(z), No*pow( 1+z, 3), usecaseB}; + + F.function = &MHR_cf; + F.params=&p; + lower_limit = log(0.01); + upper_limit = log(200); + + int status; + + gsl_set_error_handler_off(); + + status = gsl_integration_qag (&F, lower_limit, upper_limit, 0, rel_tol, + 1000, GSL_INTEG_GAUSS61, w, &result, &error); + + if(status!=0) { + LOG_ERROR("gsl integration error occured!"); + LOG_ERROR("(function argument): lower_limit=%e upper_limit=%e rel_tol=%e result=%e error=%e",lower_limit,upper_limit,rel_tol,result,error); + LOG_ERROR("data: z=%e gamma12_bg=%e T4=%e A_MHR(z)=%e",z,gamma12_bg,T4,A_MHR(z)); + LOG_ERROR("data: C_MHR(z)=%e beta_MHR(z)=%e nH=%e usecaseB=%d\n",C_MHR(z),beta_MHR(z),No*pow( 1+z, 3),usecaseB); + GSL_ERROR(status); + } + + gsl_integration_workspace_free (w); + + return result; +} + +double residual_neutral_hydrogen(double z, double gamma12_bg, double T4, int usecaseB){ + double result, error, lower_limit, upper_limit, A, C_0, beta, avenH; + gsl_function F; + double rel_tol = 0.01; //<- relative tolerance + gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000); + RR_par p = {z, gamma12_bg, T4, A_MHR(z), C_MHR(z), beta_MHR(z), No*pow( 1+z, 3), usecaseB}; + + F.function = &MHR_rnh; + F.params=&p; + lower_limit = log(0.01); + upper_limit = log(200); + + int status; + + gsl_set_error_handler_off(); + + status = gsl_integration_qag (&F, lower_limit, upper_limit, 0, rel_tol, + 1000, GSL_INTEG_GAUSS61, w, &result, &error); + + if(status!=0) { + LOG_ERROR("gsl integration error occured!"); + LOG_ERROR("(function argument): lower_limit=%e upper_limit=%e rel_tol=%e result=%e error=%e",lower_limit,upper_limit,rel_tol,result,error); + LOG_ERROR("data: z=%e gamma12_bg=%e T4=%e A_MHR(z)=%e",z,gamma12_bg,T4,A_MHR(z)); + LOG_ERROR("data: C_MHR(z)=%e beta_MHR(z)=%e nH=%e usecaseB=%d\n",C_MHR(z),beta_MHR(z),No*pow( 1+z, 3),usecaseB); + GSL_ERROR(status); + } + + gsl_integration_workspace_free (w); + + return result; +} + double aux_function(double del, void *params){ double result; double z = *(double *) params; diff --git a/src/py21cmfast/wrapper.py b/src/py21cmfast/wrapper.py index 5716f77c..abc3d930 100644 --- a/src/py21cmfast/wrapper.py +++ b/src/py21cmfast/wrapper.py @@ -2793,6 +2793,7 @@ def run_lightcone( global_q = {quantity: np.zeros(len(scrollz)) for quantity in global_quantities} mean_f_colls = np.zeros(len(scrollz)) mean_f_coll_MINIs = np.zeros(len(scrollz)) + node_redshifts_adjusted = np.zeros(len(scrollz)) pf = None perturb_files = [] @@ -2863,16 +2864,11 @@ def run_lightcone( write=write, # quick hack for running MultiNest cleanup=(cleanup and iz == (len(scrollz) - 1)), ) - mean_f_colls[iz] = ( - ib2.mean_f_coll_PC if flag_options.PHOTON_CONS else ib2.mean_f_coll - ) - mean_f_coll_MINIs[iz] = ( - ib2.mean_f_coll_MINI_PC - if ib2.flag_options.PHOTON_CONS - else ib2.mean_f_coll_MINI - ) + mean_f_colls[iz] = ib2.mean_f_coll + mean_f_coll_MINIs[iz] = ib2.mean_f_coll_MINI log10_mturnovers[iz] = ib2.log10_Mturnover_ave log10_mturnovers_mini[iz] = ib2.log10_Mturnover_MINI_ave + node_redshifts_adjusted[iz] = ib2.redshift_PC bt2 = brightness_temperature( ionized_box=ib2, @@ -2927,7 +2923,7 @@ def run_lightcone( # Save mean/global quantities for quantity in global_quantities: - if quantity == "Nion_box": + if quantity in ["Nion_box", "C_box", "nHI_box"]: global_q[quantity][iz] = np.ma.masked_equal( getattr(outs[_fld_names[quantity]][1], quantity), 0 ).mean() @@ -2998,6 +2994,7 @@ def run_lightcone( init_box.random_seed, lc, node_redshifts=scrollz, + node_redshifts_adjusted=node_redshifts_adjusted, global_quantities=global_q, mean_f_colls=mean_f_colls, mean_f_coll_MINIs=mean_f_coll_MINIs,