compiles

src/py21cmfast/drivers/param_config.py

@@ -129,6 +129,12 @@ def _flag_options_validator(self, att, val):
                     "Simultaneously."
                 )
                 raise NotImplementedError(msg)
+            if val.USE_HALO_FIELD and "GAMMA-APPROX" in (
+                self.user_params.INTEGRATION_METHOD_MINI,
+                self.user_params.INTEGRATION_METHOD_ATOMIC,
+            ):
+                msg = "the USE_HALO_FIELD mode uses more complex scaling relations which are not compatible with 'GAMMA-APPROX' integration"
+                raise ValueError(msg)
 
         if val.USE_EXP_FILTER and not val.USE_HALO_FIELD:
             warnings.warn("USE_EXP_FILTER has no effect unless USE_HALO_FIELD is true")

src/py21cmfast/src/HaloBox.c

@@ -18,7 +18,6 @@
 #include "interp_tables.h"
 #include "thermochem.h"
 #include "hmf.h"
-#include "photoncons.h"
 #include "scaling_relations.h"
 
 #include "HaloBox.h"
@@ -107,8 +106,8 @@ int get_box_averages(double M_min, double M_max, double M_turn_a, double M_turn_
     double prefactor_wsfr_mini = prefactor_sfr_mini * consts->fesc_7;
 
     double mass_intgrl;
-    double intgrl_fesc_weighted, intgrl_stars_only, integral_xray;
-    double intgrl_fesc_weighted_mini=0., intgrl_stars_only_mini=0.;
+    double intgrl_fesc_weighted, intgrl_stars_only;
+    double intgrl_fesc_weighted_mini=0., intgrl_stars_only_mini=0., integral_xray=0.;
 
     //NOTE: we use the atomic method for all halo mass/count here
     mass_intgrl = Fcoll_General(consts->redshift,lnMmin,lnMmax);
@@ -128,7 +127,8 @@ int get_box_averages(double M_min, double M_max, double M_turn_a, double M_turn_
     if(flag_options_global->USE_TS_FLUCT){
         integral_xray = Xray_General(consts->redshift, lnMmin, lnMmax, M_turn_m, M_turn_a,
                                             consts->alpha_star, consts->alpha_star_mini, consts->fstar_10,
-                                            consts->fesc_7, consts->Mlim_Fstar, consts->Mlim_Fstar_mini);
+                                            consts->fesc_7, consts->l_x, consts->l_x_mini, consts->t_h, consts->t_star,
+                                            consts->Mlim_Fstar, consts->Mlim_Fstar_mini);
     }
 
     averages_out->halo_mass = mass_intgrl * prefactor_mass;
@@ -301,8 +301,8 @@ int set_fixed_grids(double M_min, double M_max, InitialConditions *ini_boxes,
         double dens;
         double l10_mturn_a,l10_mturn_m;
         double mass_intgrl, h_count;
-        double intgrl_fesc_weighted, intgrl_stars_only, integral_xray;
-        double intgrl_fesc_weighted_mini=0., intgrl_stars_only_mini=0.;
+        double intgrl_fesc_weighted, intgrl_stars_only;
+        double intgrl_fesc_weighted_mini=0., intgrl_stars_only_mini=0., integral_xray=0;
         double dens_fac;
 
 #pragma omp for reduction(+:hm_sum,sm_sum,sm_sum_mini,sfr_sum,sfr_sum_mini,xray_sum,nion_sum,wsfr_sum)
@@ -335,7 +335,7 @@ int set_fixed_grids(double M_min, double M_max, InitialConditions *ini_boxes,
                 //MAKE A NEW TABLEdouble delta, double log10Mturn_m, double growthf, double M_min, double M_max, double M_cond, double sigma_max,
                                   //   double Mturn_a, double Mlim_Fstar, double Mlim_Fstar_MINI
                 integral_xray = EvaluateXray_Conditional(dens,l10_mturn_m,growth_z,M_min,M_max,M_cell,sigma_cell,
-                                                            l10_mturn_a,consts->Mlim_Fstar,consts->Mlim_Fstar_mini);
+                                                            l10_mturn_a,consts->t_h,consts->Mlim_Fstar,consts->Mlim_Fstar_mini);
             }
 
             grids->count[i] = (int)(h_count * M_cell * dens_fac); //NOTE: truncated

src/py21cmfast/src/hmf.c

@@ -13,6 +13,7 @@
 #include "InputParameters.h"
 #include "cosmology.h"
 #include "interp_tables.h"
+#include "scaling_relations.h"
 
 #include "hmf.h"
 
@@ -305,12 +306,12 @@ double nion_fraction_mini(double lnM, void *param_struct){
 double xray_fraction_doublePL(double lnM, void *param_struct){
     struct parameters_gsl_MF_integrals p = *(struct parameters_gsl_MF_integrals *)param_struct;
     double M = exp(lnM);
-    double Fstar = scaling_PL_limit(M,p.f_star_norm,p.alpha_star,p.Mlim_star,1e10) * exp(-p.Mturn/M);
+    double Fstar = p.f_star_norm * scaling_PL_limit(M,p.f_star_norm,p.alpha_star,p.Mlim_star,1e10) * exp(-p.Mturn_upper/M);
 
     //using the escape fraction variables for minihalos
     double Fstar_mini = 0.;
     if(flag_options_global->USE_MINI_HALOS)
-        Fstar_mini = scaling_PL_limit(M,p.f_esc_norm,p.alpha_esc,p.Mlim_esc,1e7) * exp(-p.Mturn/M - M/p.Mturn_upper);
+        Fstar_mini = p.f_esc_norm * scaling_PL_limit(M,p.f_esc_norm,p.alpha_esc,p.Mlim_esc,1e7) * exp(-p.Mturn/M - M/p.Mturn_upper);
 
     double sfr = M*Fstar/(p.t_star*p.t_h);
     double sfr_mini = M*Fstar_mini/(p.t_star*p.t_h);
@@ -545,7 +546,7 @@ double Fcollapprox_condition(double numin, double nucondition, double beta){
 //  and takes advantage of the fact that Gamma_inc(x,min) = integral_min^inf (t^(x-1)exp(-t)) dt which is satisfied for the HMF when the
 //  above approximations are made
 //Originally written by JBM within the GL integration before it was separated here and generalised to the other integrals
-double MFIntegral_Approx(double lnM_lo, double lnM_hi, struct parameters_gsl_MF_integrals params, double (*integrand)(double,void*)){
+double MFIntegral_Approx(double lnM_lo, double lnM_hi, struct parameters_gsl_MF_integrals params){
     //variables used in the calculation
     double delta,sigma_c;
     double index_base;
@@ -555,9 +556,16 @@ double MFIntegral_Approx(double lnM_lo, double lnM_hi, struct parameters_gsl_MF_
         LOG_ERROR("Ensure parameter input specifically to this function has HMF==0");
         Throw(TableGenerationError);
     }
+    if(abs(params.gamma_type) > 4 || abs(params.gamma_type)){
+        LOG_ERROR("Approximate Fcoll only works for single power-law scaling relations");
+        LOG_ERROR("These include the following General/Conditional integration functions");
+        LOG_ERROR("Nhalo, Fcoll, Nion, Nion_MINI");
+        LOG_ERROR("Something has gone wrong in the backend such that the 'Gamma-Approx'");
+        LOG_ERROR("integration method was used on a more complex scaling relation");
+        Throw(TableGenerationError);
+    }
     double growthf = params.growthf;
-    FIXALLTHETYPESTUFF
-    if(type < 0){
+    if(params.gamma_type < 0){
         //we are a conditional mf
         delta = params.delta;
         sigma_c = params.sigma_cond;
@@ -580,10 +588,10 @@ double MFIntegral_Approx(double lnM_lo, double lnM_hi, struct parameters_gsl_MF_
     //The below limit setting is done simply so that variables which do not conern particular integrals
     //      can be left undefined, rather than explicitly set to some value (0 or 1e20)
     //Mass and number integrals set the lower cutoff to the integral limit
-    if(fabs(type) >= 3 && lnMturn_l > lnM_lo_limit)
+    if(abs(params.gamma_type) >= 3 && lnMturn_l > lnM_lo_limit)
         lnM_lo_limit = lnMturn_l;
     //non-minihalo integrals set the upper cutoff to the integral limit
-    if(fabs(type) == 4 && lnMturn_u < lnM_hi_limit)
+    if(abs(params.gamma_type) == 4 && lnMturn_u < lnM_hi_limit)
         lnM_hi_limit = lnMturn_u;
 
     //it is possible for the lower turnover (LW crit or reion feedback)
@@ -593,10 +601,10 @@ double MFIntegral_Approx(double lnM_lo, double lnM_hi, struct parameters_gsl_MF_
     }
 
     //n_ion or MINI
-    if(fabs(type) >= 3)
+    if(abs(params.gamma_type) >= 3)
         index_base = params.alpha_star + params.alpha_esc;
     //fcoll
-    else if(fabs(type)==2)
+    else if(abs(params.gamma_type)==2)
         index_base = 0.;
     //nhalo
     else
@@ -633,7 +641,7 @@ double MFIntegral_Approx(double lnM_lo, double lnM_hi, struct parameters_gsl_MF_
     //NOTES: For speed the minihalos ignore the condition mass limit (assuming nu_hi_limit(tilde) < nu_condition (no tilde))
     //    and never get into the high mass power law (nu_hi_limit < nu_pivot1 (both tilde))
     //ACGs ignore the upper mass limit (no upper turnover), both assume the condition is above the highest pivot
-    if(fabs(type) == 4){
+    if(abs(params.gamma_type) == 4){
       // re-written for further speedups
       if (nu_hi_limit <= nu_pivot2){ //if both are below pivot2 don't bother adding and subtracting the high contribution
         fcoll += (Fcollapprox(nu_lo_limit,beta3))*pow(nu_pivot2_umf,-beta3);
@@ -685,7 +693,7 @@ double IntegratedNdM(double lnM_lo, double lnM_hi, struct parameters_gsl_MF_inte
     if(method==1)
         return IntegratedNdM_GL(lnM_lo, lnM_hi, params, integrand);
     if(method==2)
-        return MFIntegral_Approx(lnM_lo, lnM_hi, params, integrand);
+        return MFIntegral_Approx(lnM_lo, lnM_hi, params);
 
     LOG_ERROR("Invalid integration method %d",method);
     Throw(ValueError);
@@ -721,18 +729,20 @@ double FgtrM_wsigma(double z, double sig){
 
 double Nhalo_General(double z, double lnM_min, double lnM_max){
     struct parameters_gsl_MF_integrals integral_params = {
-                .redshift = z,
-                .growthf = dicke(z),
-                .HMF = user_params_global->HMF,
+        .redshift = z,
+        .growthf = dicke(z),
+        .HMF = user_params_global->HMF,
+        .gamma_type=1,
     };
     return IntegratedNdM(lnM_min, lnM_max, integral_params, &u_mf_integrand, 0);
 }
 
 double Fcoll_General(double z, double lnM_min, double lnM_max){
     struct parameters_gsl_MF_integrals integral_params = {
-                .redshift = z,
-                .growthf = dicke(z),
-                .HMF = user_params_global->HMF,
+        .redshift = z,
+        .growthf = dicke(z),
+        .HMF = user_params_global->HMF,
+        .gamma_type=2,
     };
     return IntegratedNdM(lnM_min, lnM_max, integral_params, &u_fcoll_integrand, 0);
 }
@@ -750,6 +760,7 @@ double Nion_General(double z, double lnM_Min, double lnM_Max, double MassTurnove
         .Mlim_star = Mlim_Fstar,
         .Mlim_esc = Mlim_Fesc,
         .HMF = user_params_global->HMF,
+        .gamma_type=3,
     };
     return IntegratedNdM(lnM_Min,lnM_Max,params,&u_nion_integrand,0);
 }
@@ -768,6 +779,7 @@ double Nion_General_MINI(double z, double lnM_Min, double lnM_Max, double MassTu
         .Mlim_star = Mlim_Fstar,
         .Mlim_esc = Mlim_Fesc,
         .HMF = user_params_global->HMF,
+        .gamma_type=4,
     };
     return IntegratedNdM(lnM_Min,lnM_Max,params,&u_nion_integrand_mini,0);
 }
@@ -791,6 +803,7 @@ double Xray_General(double z, double lnM_Min, double lnM_Max, double MassTurnove
         .l_x_norm_mini = l_x_mini,
         .t_h = t_h,
         .t_star = t_star,
+        .gamma_type=5,
     };
     return IntegratedNdM(lnM_Min,lnM_Max,params,&u_xray_integrand,0);
 }
@@ -801,6 +814,7 @@ double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double M_cond
         .HMF = user_params_global->HMF,
         .sigma_cond = sigma,
         .delta = delta,
+        .gamma_type=-1,
     };
 
     if(delta <= -1. || lnM1 >= log(M_cond))
@@ -822,6 +836,7 @@ double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double M_cond
         .HMF = user_params_global->HMF,
         .sigma_cond = sigma,
         .delta = delta,
+        .gamma_type=-2,
     };
 
     if(delta <= -1. || lnM1 >= log(M_cond))
@@ -852,6 +867,7 @@ double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M
         .HMF = user_params_global->HMF,
         .sigma_cond = sigma2,
         .delta = delta2,
+        .gamma_type=-4,
     };
 
     if(delta2 <= -1. || lnM1 >= log(M_cond))
@@ -889,6 +905,7 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond
         .HMF = user_params_global->HMF,
         .sigma_cond = sigma2,
         .delta = delta2,
+        .gamma_type=-3,
     };
 
     if(delta2 <= -1. || lnM1 >= log(M_cond))
@@ -909,12 +926,14 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond
     return IntegratedNdM(lnM1,lnM2,params,&c_nion_integrand_mini,method);
 }
 
-double Xray_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond, double sigma2, double delta2, double MassTurnover,
+double Xray_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond, double sigma2, double delta2,
+                         double MassTurnover, double MassTurnover_upper,
                         double Alpha_star, double Alpha_star_mini, double Fstar10, double Fstar7, double Mlim_Fstar,
                         double Mlim_Fstar_mini, double l_x, double l_x_mini, double t_h, double t_star, int method){
     struct parameters_gsl_MF_integrals params = {
         .growthf = growthf,
         .Mturn = MassTurnover,
+        .Mturn_upper = MassTurnover_upper,
         .alpha_star = Alpha_star,
         .alpha_esc = Alpha_star_mini,
         .f_star_norm = Fstar10,
@@ -926,6 +945,7 @@ double Xray_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond
         .l_x_norm_mini = l_x_mini,
         .t_h = t_h,
         .t_star = t_star,
+        .gamma_type=-5,
     };
 
     if(delta2 <= -1. || lnM1 >= log(M_cond))

src/py21cmfast/src/hmf.h

@@ -37,6 +37,9 @@ struct parameters_gsl_MF_integrals{
     double l_x_norm_mini;
     double t_h;
     double t_star;
+
+    //needed for FAST_FCOLL gamma approximations
+    int gamma_type;
 };
 
 /* HMF Integrals */
@@ -49,10 +52,17 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond
 double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M_cond, double sigma2, double delta2, double MassTurnover,
                             double MassTurnover_upper, double Alpha_star, double Alpha_esc, double Fstar7,
                             double Fesc7, double Mlim_Fstar, double Mlim_Fesc, int method);
+double Xray_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond, double sigma2, double delta2,
+                         double MassTurnover, double MassTurnover_upper,
+                        double Alpha_star, double Alpha_star_mini, double Fstar10, double Fstar7, double Mlim_Fstar,
+                        double Mlim_Fstar_mini, double l_x, double l_x_mini, double t_h, double t_star, int method);
 double Nion_General(double z, double lnM_Min, double lnM_Max, double MassTurnover, double Alpha_star, double Alpha_esc, double Fstar10,
                      double Fesc10, double Mlim_Fstar, double Mlim_Fesc);
 double Nion_General_MINI(double z, double lnM_Min, double lnM_Max, double MassTurnover, double MassTurnover_upper, double Alpha_star,
                          double Alpha_esc, double Fstar7_MINI, double Fesc7_MINI, double Mlim_Fstar, double Mlim_Fesc);
+double Xray_General(double z, double lnM_Min, double lnM_Max, double MassTurnover, double MassTurnover_upper, double Alpha_star,
+                     double Alpha_star_mini, double Fstar10, double Fstar7, double l_x, double l_x_mini, double t_h,
+                     double t_star, double Mlim_Fstar, double Mlim_Fstar_mini);
 double Nhalo_General(double z, double lnM_min, double lnM_max);
 double Fcoll_General(double z, double lnM_min, double lnM_max);
 double unconditional_mf(double growthf, double lnM, double z, int HMF);

src/py21cmfast/src/interp_tables.c

@@ -485,7 +485,8 @@ void initialise_SFRD_Conditional_table(double min_density, double max_density, d
 //This function initialises one table, for table Rx arrays I will call this function in a loop
 void initialise_Xray_Conditional_table(double min_density, double max_density, double growthf,
                                     float Mcrit_atom, double Mmin, double Mmax, double Mcond, float Alpha_star, float Alpha_star_mini,
-                                    float Fstar10, float Fstar7_MINI, int method, int method_mini, bool minihalos){
+                                    float Fstar10, float Fstar7_MINI, double l_x, double l_x_mini, double t_h, double t_star,
+                                    int method, int method_mini, bool minihalos){
     float Mlim_Fstar,sigma2;
     float Mlim_Fstar_MINI=0.;
     int i,k;
@@ -534,7 +535,9 @@ void initialise_Xray_Conditional_table(double min_density, double max_density, d
             curr_dens = min_density + (float)i/((float)NDELTA-1.)*(max_density - min_density);
             if(!minihalos){
                 Xray_conditional_table_1D.y_arr[i] = log(Xray_ConditionalM(growthf,lnMmin,lnMmax,Mcond,sigma2,curr_dens,
-                                                Mcrit_atom,Alpha_star,0.,Fstar10,0.,Mlim_Fstar,0.,user_params_global->INTEGRATION_METHOD_ATOMIC));
+                                            0.,Mcrit_atom,Alpha_star,0.,Fstar10,0.,Mlim_Fstar,0.,
+                                            l_x, l_x_mini, t_h, t_star,
+                                            user_params_global->INTEGRATION_METHOD_ATOMIC));
 
                 if(Xray_conditional_table_1D.y_arr[i] < -50.)
                     Xray_conditional_table_1D.y_arr[i] = -50.;
@@ -545,6 +548,7 @@ void initialise_Xray_Conditional_table(double min_density, double max_density, d
                 Xray_conditional_table_2D.z_arr[i][k] = log(Xray_ConditionalM(growthf,lnMmin,lnMmax,Mcond,sigma2,
                                             curr_dens,MassTurnover[k],Mcrit_atom,
                                             Alpha_star,Alpha_star_mini,Fstar10,Fstar7_MINI,Mlim_Fstar,Mlim_Fstar_MINI,
+                                            l_x, l_x_mini, t_h, t_star,
                                             user_params_global->INTEGRATION_METHOD_MINI)); //Using mini integration method for both
 
                 if(Xray_conditional_table_2D.z_arr[i][k] < -50.)
@@ -855,7 +859,7 @@ void free_conditional_tables(){
     free_RGTable2D_f(&Nion_conditional_table_MINI);
     free_RGTable2D_f(&Nion_conditional_table_prev);
     free_RGTable2D_f(&Nion_conditional_table_MINI_prev);
-    free_RGTable2D_f(&Xray_conditional_table_1D);
+    free_RGTable1D_f(&Xray_conditional_table_1D);
     free_RGTable2D_f(&Xray_conditional_table_2D);
 }
 
@@ -866,7 +870,7 @@ void free_global_tables(){
     free_RGTable2D(&Nion_z_table_MINI);
     free_RGTable1D(&fcoll_z_table);
     free_RGTable1D(&Xray_z_table_1D);
-    free_RGTable1D(&Xray_z_table_2D);
+    free_RGTable2D(&Xray_z_table_2D);
 }
 
 //JD: moving the interp table evaluations here since some of them are needed in nu_tau_one
@@ -992,17 +996,19 @@ double EvaluateNion_Conditional_MINI(double delta, double log10Mturn_m, double g
                                 Mlim_Fesc,user_params_global->INTEGRATION_METHOD_MINI);
 }
 
-double EvaluateXray_Conditional_MINI(double delta, double log10Mturn_m, double growthf, double M_min, double M_max, double M_cond, double sigma_max,
-                                     double Mturn_a, double Mlim_Fstar, double Mlim_Fstar_MINI){
+double EvaluateXray_Conditional(double delta, double log10Mturn_m, double growthf, double M_min, double M_max, double M_cond, double sigma_max,
+                                     double Mturn_a, double t_h, double Mlim_Fstar, double Mlim_Fstar_MINI){
     if(user_params_global->USE_INTERPOLATION_TABLES){
         if(flag_options_global->USE_MINI_HALOS)
             return exp(EvaluateRGTable2D_f(delta,log10Mturn_m,&Xray_conditional_table_2D));
         return exp(EvaluateRGTable1D_f(delta,&Xray_conditional_table_1D));
     }
 
-    return Xray_ConditionalM(growthf,log(M_min),log(M_max),M_cond,sigma_max,delta,pow(10,log10Mturn_m),Mturn_a,astro_params_global->ALPHA_STAR,
-                                astro_params_global->ALPHA_STAR_MINI,astro_params_global->F_STAR10,astro_params_global->F_STAR7_MINI,
-                                Mlim_Fstar, Mlim_Fstar_MINI, user_params_global->INTEGRATION_METHOD_MINI);
+    return Xray_ConditionalM(growthf, log(M_min), log(M_max), M_cond, sigma_max, delta, pow(10,log10Mturn_m),
+                                Mturn_a, astro_params_global->ALPHA_STAR, astro_params_global->ALPHA_STAR_MINI,
+                                astro_params_global->F_STAR10, astro_params_global->F_STAR7_MINI,
+                                Mlim_Fstar, Mlim_Fstar_MINI, astro_params_global->L_X, astro_params_global->L_X_MINI,t_h,
+                                astro_params_global->t_STAR, user_params_global->INTEGRATION_METHOD_MINI);
 }
 
 double EvaluateFcoll_delta(double delta, double growthf, double sigma_min, double sigma_max){