diff --git a/src/py21cmfast/drivers/param_config.py b/src/py21cmfast/drivers/param_config.py index a09b56cd2..cb4ddc305 100644 --- a/src/py21cmfast/drivers/param_config.py +++ b/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") diff --git a/src/py21cmfast/src/HaloBox.c b/src/py21cmfast/src/HaloBox.c index a9b226f13..52335c278 100644 --- a/src/py21cmfast/src/HaloBox.c +++ b/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 diff --git a/src/py21cmfast/src/hmf.c b/src/py21cmfast/src/hmf.c index 8af2e27ed..7f33aa08c 100644 --- a/src/py21cmfast/src/hmf.c +++ b/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)) diff --git a/src/py21cmfast/src/hmf.h b/src/py21cmfast/src/hmf.h index e053134e1..033fcac7c 100644 --- a/src/py21cmfast/src/hmf.h +++ b/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); diff --git a/src/py21cmfast/src/interp_tables.c b/src/py21cmfast/src/interp_tables.c index 8f2b54bb6..d2413935e 100644 --- a/src/py21cmfast/src/interp_tables.c +++ b/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){ diff --git a/src/py21cmfast/src/interp_tables.h b/src/py21cmfast/src/interp_tables.h index 87e8fc4c6..ba06dcbe0 100644 --- a/src/py21cmfast/src/interp_tables.h +++ b/src/py21cmfast/src/interp_tables.h @@ -33,7 +33,12 @@ double EvaluateNion_Conditional(double delta, double log10Mturn, double growthf, double Mlim_Fstar, double Mlim_Fesc, bool prev); double EvaluateNion_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_Fesc, bool prev); - +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, double l_x, double l_x_mini, double t_h, double t_star, + int method, int method_mini, bool minihalos); +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); void initialise_SFRD_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); diff --git a/src/py21cmfast/src/scaling_relations.c b/src/py21cmfast/src/scaling_relations.c index 00e075b52..93810a971 100644 --- a/src/py21cmfast/src/scaling_relations.c +++ b/src/py21cmfast/src/scaling_relations.c @@ -4,6 +4,7 @@ #include #include #include +#include #include "cexcept.h" #include "exceptions.h" #include "logger.h" @@ -11,6 +12,9 @@ #include "Constants.h" #include "InputParameters.h" #include "cosmology.h" +#include "photoncons.h" +#include "hmf.h" +#include "thermochem.h" #include "scaling_relations.h" @@ -28,7 +32,8 @@ void set_scaling_constants(double redshift, AstroParams *astro_params, FlagOptio consts->alpha_upper = astro_params->UPPER_STELLAR_TURNOVER_INDEX; consts->pivot_upper = astro_params->UPPER_STELLAR_TURNOVER_MASS; - consts->upper_pivot_ratio = pow(astro_params->UPPER_STELLAR_TURNOVER_MASS/1e10,astro_params->ALPHA_STAR); + consts->upper_pivot_ratio = pow(consts->pivot_upper/1e10,consts->alpha_star) + \ + pow(consts->pivot_upper/1e10,consts->alpha_upper); consts->fstar_7 = astro_params->F_STAR7_MINI; consts->alpha_star_mini = astro_params->ALPHA_STAR_MINI; @@ -81,6 +86,14 @@ by keeping them in log. Since they are called within integrals they don't use th ScalingConstants. Instead pulling from the GSL integration parameter structs */ +double scaling_single_PL(double M, double alpha, double pivot){ + return pow(M/pivot,alpha); +} + +double log_scaling_single_PL(double lnM, double alpha, double ln_pivot){ + return alpha*(lnM - ln_pivot); +} + //single power-law with provided limit for scaling == 1, returns f(M)/f(pivot) == norm //limit is provided ahead of time for efficiency double scaling_PL_limit(double M, double norm, double alpha, double pivot, double limit){ @@ -88,7 +101,7 @@ double scaling_PL_limit(double M, double norm, double alpha, double pivot, doubl return 1/norm; //if alpha is zero, this returns 1 as expected (note strict inequalities above) - return scaling_single_PL(M,alpha_pivot); + return scaling_single_PL(M,alpha,pivot); } //log version for possible acceleration @@ -100,7 +113,7 @@ double log_scaling_PL_limit(double lnM, double ln_norm, double alpha, double ln_ return log_scaling_single_PL(lnM,alpha,ln_pivot); } -//concave-down double power-law, we pass pivot_ratio == (pivot_hi/pivot_lo)^alpha_lo +//concave-down double power-law, we pass pivot_ratio == denominator(M=pivot_lo) // to save pow() calls. Due to the double power-law we gain little from a log verison // returns f(M)/f(M==pivot_lo) double scaling_double_PL(double M, double alpha_lo, double pivot_ratio, @@ -109,15 +122,6 @@ double scaling_double_PL(double M, double alpha_lo, double pivot_ratio, return pivot_ratio / (pow(M/pivot_hi,-alpha_lo) + pow(M/pivot_hi,-alpha_hi)); } -double scaling_single_PL(double M, double alpha, double pivot){ - return pow(M/pivot,alpha); -} - -double log_scaling_single_PL(double lnM, double alpha, double ln_pivot){ - return alpha*(lnM - ln_pivot); -} - - /* Scaling relations used in the halo sampler. Quantities calculated for a single halo mass and are summed onto grids. @@ -152,11 +156,11 @@ double lx_on_sfr_Lehmer(double metallicity){ //double power-law in Z with the low-metallicity PL fixed as constant double lx_on_sfr_doublePL(double metallicity, double lx_constant){ - double hiz_index = -0.64; //power-law index of LX/SFR at high-z - double loz_index = 0.; - double z_pivot = 0.05; // Z at which LX/SFR == lx_constant + double hi_z_index = -0.64; //power-law index of LX/SFR at high-z + double lo_z_index = 0.; + double z_pivot = 0.05; // Z at which LX/SFR == lx_constant / 2 - return lx_constant*scaling_double_PL(metallicity,loz_index.,1.,z_index,z_pivot); + return lx_constant*scaling_double_PL(metallicity,lo_z_index,1.,hi_z_index,z_pivot); } //first order power law Lx with cross-term (e.g Kaur+22) @@ -223,11 +227,10 @@ void get_halo_stellarmass(double halo_mass, double mturn_acg, double mturn_mcg, //We don't want an upturn even with a negative ALPHA_STAR if(flag_options_global->USE_UPPER_STELLAR_TURNOVER && (f_a > fu_a)){ - fstar_mean = scaling_double_PL(halo_mass,) - fstar_mean = consts->upper_pivot_ratio / (pow(halo_mass/fu_p,-f_a)+pow(halo_mass/fu_p,-fu_a)); + fstar_mean = scaling_double_PL(halo_mass,f_a,consts->upper_pivot_ratio,fu_a,fu_p); } else{ - fstar_mean = pow(halo_mass/1e10,f_a); //PL term + fstar_mean = scaling_single_PL(halo_mass,consts->alpha_star,1e10); //PL term } f_sample = f_10 * fstar_mean * exp(-mturn_acg/halo_mass + star_rng*sigma_star - stoc_adjustment_term); //1e10 normalisation of stellar mass if(f_sample > 1.) f_sample = 1.;