From 7bd9a181f669505fb6de0f517b2d3e60e15cdaf1 Mon Sep 17 00:00:00 2001 From: James Davies Date: Thu, 5 Dec 2024 18:48:16 +0100 Subject: [PATCH] WIP fixing xray integrals --- src/py21cmfast/src/HaloBox.c | 31 +++++-- src/py21cmfast/src/hmf.c | 129 +++++++++++++++++++---------- src/py21cmfast/src/hmf.h | 45 ++-------- src/py21cmfast/src/interp_tables.c | 33 ++++---- src/py21cmfast/src/interp_tables.h | 4 +- src/py21cmfast/wrapper/cfuncs.py | 8 +- 6 files changed, 134 insertions(+), 116 deletions(-) diff --git a/src/py21cmfast/src/HaloBox.c b/src/py21cmfast/src/HaloBox.c index 52335c278..b655d937a 100644 --- a/src/py21cmfast/src/HaloBox.c +++ b/src/py21cmfast/src/HaloBox.c @@ -202,10 +202,10 @@ int set_fixed_grids(double M_min, double M_max, InitialConditions *ini_boxes, //find grid limits for tables double min_density = 0.; double max_density = 0.; - double min_log10_mturn_a = 0.; - double min_log10_mturn_m = 0.; - double max_log10_mturn_a = 0.; - double max_log10_mturn_m = 0.; + double min_log10_mturn_a = log10(global_params.M_MAX_INTEGRAL); + double min_log10_mturn_m = log10(global_params.M_MAX_INTEGRAL); + double max_log10_mturn_a = log10(astro_params_global->M_TURN); + double max_log10_mturn_m = log10(astro_params_global->M_TURN); float *mturn_a_grid = calloc(HII_TOT_NUM_PIXELS,sizeof(float)); float *mturn_m_grid = calloc(HII_TOT_NUM_PIXELS,sizeof(float)); #pragma omp parallel num_threads(user_params_global->N_THREADS) @@ -293,6 +293,17 @@ int set_fixed_grids(double M_min, double M_max, InitialConditions *ini_boxes, flag_options_global->USE_MINI_HALOS, false); initialise_dNdM_tables(min_density, max_density, lnMmin, lnMmax, growth_z, lnMcell, false); + if(flag_options_global->USE_TS_FLUCT){ + initialise_Xray_Conditional_table( + min_density, max_density, consts->redshift, growth_z, consts->mturn_a_nofb, M_min, M_max, M_cell, + consts->alpha_star, consts->alpha_star_mini, + consts->fstar_10, consts->fstar_7, consts->l_x, consts->l_x_mini, consts->t_h, + consts->t_star, + user_params_global->INTEGRATION_METHOD_ATOMIC, + user_params_global->INTEGRATION_METHOD_MINI, + flag_options_global->USE_MINI_HALOS + ); + } } #pragma omp parallel num_threads(user_params_global->N_THREADS) @@ -334,7 +345,7 @@ int set_fixed_grids(double M_min, double M_max, InitialConditions *ini_boxes, if(flag_options_global->USE_TS_FLUCT){ //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, + integral_xray = EvaluateXray_Conditional(dens,l10_mturn_m,consts->redshift,growth_z,M_min,M_max,M_cell,sigma_cell, l10_mturn_a,consts->t_h,consts->Mlim_Fstar,consts->Mlim_Fstar_mini); } @@ -842,10 +853,12 @@ int test_halo_props(double redshift, UserParams *user_params, CosmoParams *cosmo halo_props_out[12*i_halo + 10] = M_turn_r; halo_props_out[12*i_halo + 11] = out_props.metallicity; - LOG_ULTRA_DEBUG("HM %.2e SM %.2e SF %.2e NI %.2e LX %.2e",out_props.halo_mass,out_props.stellar_mass,out_props.halo_sfr,out_props.n_ion,out_props.halo_xray); - LOG_ULTRA_DEBUG("MINI: SM %.2e SF %.2e WSF %.2e Mta %.2e Mtm %.2e Mtr %.2e",out_props.stellar_mass_mini,out_props.sfr_mini,out_props.fescweighted_sfr, - out_props.m_turn_acg,out_props.m_turn_mcg,out_props.m_turn_reion); - LOG_ULTRA_DEBUG("RNG: STAR %.2e SFR %.2e XRAY %.2e",in_props[0],in_props[1],in_props[2]); + if (i_halo < 10){ + LOG_ULTRA_DEBUG("HM %.2e SM %.2e SF %.2e NI %.2e LX %.2e",out_props.halo_mass,out_props.stellar_mass,out_props.halo_sfr,out_props.n_ion,out_props.halo_xray); + LOG_ULTRA_DEBUG("MINI: SM %.2e SF %.2e WSF %.2e",out_props.stellar_mass_mini,out_props.sfr_mini,out_props.fescweighted_sfr); + LOG_ULTRA_DEBUG("Mturns ACG %.2e MCG %.2e Reion %.2e",out_props.m_turn_acg,out_props.m_turn_mcg,out_props.m_turn_reion); + LOG_ULTRA_DEBUG("RNG: STAR %.2e SFR %.2e XRAY %.2e",in_props[0],in_props[1],in_props[2]); + } } } } diff --git a/src/py21cmfast/src/hmf.c b/src/py21cmfast/src/hmf.c index 7f33aa08c..5277c9504 100644 --- a/src/py21cmfast/src/hmf.c +++ b/src/py21cmfast/src/hmf.c @@ -29,6 +29,42 @@ static float xi_GL[NGL_INT+1], wi_GL[NGL_INT+1]; static float GL_limit[2] = {0}; + +//Parameters used for gsl integral on the mass function +struct parameters_gsl_MF_integrals{ + //parameters for all MF integrals + double redshift; + double growthf; + int HMF; + + //Conditional parameters + double sigma_cond; + double delta; + + //SFR additions + double Mturn; + double f_star_norm; + double alpha_star; + double Mlim_star; + + //Nion additions + double f_esc_norm; + double alpha_esc; + double Mlim_esc; + + //Minihalo additions + double Mturn_upper; + + //X-ray additions + double l_x_norm; + double l_x_norm_mini; + double t_h; + double t_star; + + //needed for FAST_FCOLL gamma approximations + int gamma_type; +}; + /* sheth correction to delta crit */ double sheth_delc_dexm(double del, double sig){ return sqrt(SHETH_a)*del*(1. + global_params.SHETH_b*pow(sig*sig/(SHETH_a*del*del), global_params.SHETH_c)); @@ -285,16 +321,16 @@ double dNdlnM_WatsonFOF_z(double z, double growthf, double lnM){ //scaling relation for M_halo --> n_ion used in integrands double nion_fraction(double lnM, void *param_struct){ struct parameters_gsl_MF_integrals p = *(struct parameters_gsl_MF_integrals *)param_struct; - double Fstar = log_scaling_PL_limit(lnM,p.f_star_norm,p.alpha_star,p.Mlim_star,10*LN10); - double Fesc = log_scaling_PL_limit(lnM,p.f_esc_norm,p.alpha_esc,p.Mlim_esc,10*LN10); + double Fstar = log_scaling_PL_limit(lnM,p.f_star_norm,p.alpha_star,10*LN10,p.Mlim_star); + double Fesc = log_scaling_PL_limit(lnM,p.f_esc_norm,p.alpha_esc,10*LN10,p.Mlim_esc); return exp(Fstar + Fesc - p.Mturn/exp(lnM) + lnM); } double nion_fraction_mini(double lnM, void *param_struct){ struct parameters_gsl_MF_integrals p = *(struct parameters_gsl_MF_integrals *)param_struct; - double Fstar = log_scaling_PL_limit(lnM,p.f_star_norm,p.alpha_star,p.Mlim_star,7*LN10); - double Fesc = log_scaling_PL_limit(lnM,p.f_esc_norm,p.alpha_esc,p.Mlim_esc,7*LN10); + double Fstar = log_scaling_PL_limit(lnM,p.f_star_norm,p.alpha_star,7*LN10,p.Mlim_star); + double Fesc = log_scaling_PL_limit(lnM,p.f_esc_norm,p.alpha_esc,7*LN10,p.Mlim_esc); double M = exp(lnM); return exp(Fstar + Fesc - M/p.Mturn_upper - p.Mturn/M + lnM); @@ -306,22 +342,24 @@ 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 = p.f_star_norm * scaling_PL_limit(M,p.f_star_norm,p.alpha_star,p.Mlim_star,1e10) * exp(-p.Mturn_upper/M); + double Fstar = exp(log_scaling_PL_limit(lnM,p.f_star_norm,p.alpha_star,10*LN10,p.Mlim_star) - p.Mturn_upper/M + p.f_star_norm); //using the escape fraction variables for minihalos double Fstar_mini = 0.; if(flag_options_global->USE_MINI_HALOS) - 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); + Fstar_mini = exp(log_scaling_PL_limit(lnM,p.f_esc_norm,p.alpha_esc,7*LN10,p.Mlim_esc) - p.Mturn/M - M/p.Mturn_upper + p.f_esc_norm); - double sfr = M*Fstar/(p.t_star*p.t_h); - double sfr_mini = M*Fstar_mini/(p.t_star*p.t_h); + double stars = M*Fstar*cosmo_params_global->OMb/cosmo_params_global->OMm; + double stars_mini = M*Fstar_mini*cosmo_params_global->OMb/cosmo_params_global->OMm; + double sfr = stars/(p.t_star*p.t_h); + double sfr_mini = stars_mini/(p.t_star*p.t_h); double metallicity; - get_halo_metallicity(sfr+sfr_mini,M*(Fstar+Fstar_mini),p.redshift,&metallicity); + get_halo_metallicity(sfr+sfr_mini,stars+stars_mini,p.redshift,&metallicity); double l_x = get_lx_on_sfr(sfr,metallicity,p.l_x_norm); double l_x_mini = get_lx_on_sfr(sfr_mini,metallicity,p.l_x_norm_mini); - return sfr*l_x + sfr_mini*l_x_mini; + return SperYR*(sfr*l_x + sfr_mini*l_x_mini); } double conditional_mf(double growthf, double lnM, double delta, double sigma, int HMF){ @@ -447,7 +485,7 @@ double IntegratedNdM_QAG(double lnM_lo, double lnM_hi, struct parameters_gsl_MF_ LOG_ERROR("sigma=%.3e delta=%.3e",params.sigma_cond,params.delta); LOG_ERROR("Mturn_lo=%.3e f*=%.3e a*=%.3e Mlim*=%.3e",params.Mturn,params.f_star_norm,params.alpha_star,params.Mlim_star); LOG_ERROR("f_escn=%.3e a_esc=%.3e Mlim_esc=%.3e",params.f_esc_norm,params.alpha_esc,params.Mlim_esc); - LOG_ERROR("Mturn_hi %.3e",params.Mturn_upper); + LOG_ERROR("Mturn_hi %.3e gamma_type %d",params.Mturn_upper,params.gamma_type); CATCH_GSL_ERROR(status); } @@ -808,7 +846,7 @@ double Xray_General(double z, double lnM_Min, double lnM_Max, double MassTurnove return IntegratedNdM(lnM_Min,lnM_Max,params,&u_xray_integrand,0); } -double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double M_cond, double sigma, double delta, int method){ +double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double lnM_cond, double sigma, double delta, int method){ struct parameters_gsl_MF_integrals params = { .growthf = growthf, .HMF = user_params_global->HMF, @@ -817,12 +855,12 @@ double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double M_cond .gamma_type=-1, }; - if(delta <= -1. || lnM1 >= log(M_cond)) + if(delta <= -1. || lnM1 >= lnM_cond) return 0.; //return 1 halo AT THE CONDITION MASS if delta is exceeded if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)){ - if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) - return 1./M_cond; + if(lnM_cond*(1-FRACT_FLOAT_ERR) <= lnM2) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) + return 1./lnM_cond; else return 0.; } @@ -830,7 +868,7 @@ double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double M_cond return IntegratedNdM(lnM1,lnM2,params,&c_mf_integrand, method); } -double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double M_cond, double sigma, double delta, int method){ +double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double lnM_cond, double sigma, double delta, int method){ struct parameters_gsl_MF_integrals params = { .growthf = growthf, .HMF = user_params_global->HMF, @@ -839,11 +877,11 @@ double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double M_cond .gamma_type=-2, }; - if(delta <= -1. || lnM1 >= log(M_cond)) + if(delta <= -1. || lnM1 >= lnM_cond) return 0.; //return 100% of mass AT THE CONDITION MASS if delta is exceeded if(delta > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma,growthf)){ - if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) + if(lnM_cond*(1-FRACT_FLOAT_ERR) <= lnM2) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) return 1.; else return 0.; @@ -851,7 +889,7 @@ double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double M_cond return IntegratedNdM(lnM1,lnM2,params,&c_fcoll_integrand, method); } -double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M_cond, double sigma2, double delta2, double MassTurnover, +double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double lnM_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){ struct parameters_gsl_MF_integrals params = { @@ -860,24 +898,24 @@ double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M .Mturn_upper = MassTurnover_upper, .alpha_star = Alpha_star, .alpha_esc = Alpha_esc, - .f_star_norm = Fstar7, - .f_esc_norm = Fesc7, - .Mlim_star = Mlim_Fstar, - .Mlim_esc = Mlim_Fesc, + .f_star_norm = log(Fstar7), + .f_esc_norm = log(Fesc7), + .Mlim_star = log(Mlim_Fstar), + .Mlim_esc = log(Mlim_Fesc), .HMF = user_params_global->HMF, .sigma_cond = sigma2, .delta = delta2, .gamma_type=-4, }; - if(delta2 <= -1. || lnM1 >= log(M_cond)) + if(delta2 <= -1. || lnM1 >= lnM_cond) return 0.; //return 1 halo at the condition mass if delta is exceeded //NOTE: this will almost always be zero, due to the upper turover, // however this replaces an integral so it won't be slow if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)){ - if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) - return nion_fraction_mini(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) + if(lnM_cond*(1-FRACT_FLOAT_ERR) <= lnM2) //this limit is not ideal, but covers floating point errors when we set lnM2 == log(M_cond) + return nion_fraction_mini(lnM_cond,¶ms) / exp(lnM_cond); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) else return 0.; } @@ -890,7 +928,7 @@ double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M return IntegratedNdM(lnM1,lnM2,params,&c_nion_integrand_mini,method); } -double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond, double sigma2, double delta2, double MassTurnover, +double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double lnM_cond, double sigma2, double delta2, double MassTurnover, double Alpha_star, double Alpha_esc, double Fstar10, double Fesc10, double Mlim_Fstar, double Mlim_Fesc, int method){ struct parameters_gsl_MF_integrals params = { @@ -898,22 +936,22 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond .Mturn = MassTurnover, .alpha_star = Alpha_star, .alpha_esc = Alpha_esc, - .f_star_norm = Fstar10, - .f_esc_norm = Fesc10, - .Mlim_star = Mlim_Fstar, - .Mlim_esc = Mlim_Fesc, + .f_star_norm = log(Fstar10), + .f_esc_norm = log(Fesc10), + .Mlim_star = log(Mlim_Fstar), + .Mlim_esc = log(Mlim_Fesc), .HMF = user_params_global->HMF, .sigma_cond = sigma2, .delta = delta2, .gamma_type=-3, }; - if(delta2 <= -1. || lnM1 >= log(M_cond)) + if(delta2 <= -1. || lnM1 >= lnM_cond) return 0.; - //return 1 halo at the condition mass if delta is exceeded + //return 1 halo at the condition mass if delta is exceeded and the condition is within the integral limits if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)){ - if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) - return nion_fraction(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) + if(lnM_cond*(1-FRACT_FLOAT_ERR) <= lnM2) + return nion_fraction(lnM_cond,¶ms) / exp(lnM_cond); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) else return 0.; } @@ -922,24 +960,25 @@ double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond //NOTE: it's possible we may want to use another default if(params.HMF != 0 && params.HMF != 1 && params.HMF != 4) params.HMF = 0; - - return IntegratedNdM(lnM1,lnM2,params,&c_nion_integrand_mini,method); + return IntegratedNdM(lnM1,lnM2,params,&c_nion_integrand,method); } -double Xray_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond, double sigma2, double delta2, +double Xray_ConditionalM(double redshift, double growthf, double lnM1, double lnM2, double lnM_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){ + //re-using escape fraction for minihalo parameters struct parameters_gsl_MF_integrals params = { + .redshift = redshift, .growthf = growthf, .Mturn = MassTurnover, .Mturn_upper = MassTurnover_upper, .alpha_star = Alpha_star, .alpha_esc = Alpha_star_mini, - .f_star_norm = Fstar10, - .f_esc_norm = Fstar7, - .Mlim_star = Mlim_Fstar, - .Mlim_esc = Mlim_Fstar_mini, + .f_star_norm = log(Fstar10), + .f_esc_norm = log(Fstar7), + .Mlim_star = log(Mlim_Fstar), + .Mlim_esc = log(Mlim_Fstar_mini), .HMF = user_params_global->HMF, .l_x_norm = l_x, .l_x_norm_mini = l_x_mini, @@ -948,12 +987,12 @@ double Xray_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond .gamma_type=-5, }; - if(delta2 <= -1. || lnM1 >= log(M_cond)) + if(delta2 <= -1. || lnM1 >= lnM_cond) return 0.; //return 1 halo at the condition mass if delta is exceeded if(delta2 > MAX_DELTAC_FRAC*get_delta_crit(params.HMF,sigma2,growthf)){ - if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2)) - return xray_fraction_doublePL(M_cond,¶ms); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) + if(lnM_cond*(1-FRACT_FLOAT_ERR) <= lnM2) + return xray_fraction_doublePL(lnM_cond,¶ms) / exp(lnM_cond); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta) else return 0.; } diff --git a/src/py21cmfast/src/hmf.h b/src/py21cmfast/src/hmf.h index 033fcac7c..c7afe5283 100644 --- a/src/py21cmfast/src/hmf.h +++ b/src/py21cmfast/src/hmf.h @@ -7,52 +7,17 @@ #define MAX_DELTAC_FRAC (float)0.99 //max delta/deltac for the mass function integrals #define DELTA_MIN -1 //minimum delta for Lagrangian mass function integrals -//Parameters used for gsl integral on the mass function -struct parameters_gsl_MF_integrals{ - //parameters for all MF integrals - double redshift; - double growthf; - int HMF; - - //Conditional parameters - double sigma_cond; - double delta; - - //SFR additions - double Mturn; - double f_star_norm; - double alpha_star; - double Mlim_star; - - //Nion additions - double f_esc_norm; - double alpha_esc; - double Mlim_esc; - - //Minihalo additions - double Mturn_upper; - - //X-ray additions - double l_x_norm; - double l_x_norm_mini; - double t_h; - double t_star; - - //needed for FAST_FCOLL gamma approximations - int gamma_type; -}; - /* HMF Integrals */ void initialise_GL(float lnM_Min, float lnM_Max); -double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double M_cond, double sigma, double delta, int method); -double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double M_cond, double sigma, double delta, int method); -double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double M_cond, double sigma2, double delta2, double MassTurnover, +double Nhalo_Conditional(double growthf, double lnM1, double lnM2, double lnM_cond, double sigma, double delta, int method); +double Mcoll_Conditional(double growthf, double lnM1, double lnM2, double lnM_cond, double sigma, double delta, int method); +double Nion_ConditionalM(double growthf, double lnM1, double lnM2, double lnM_cond, double sigma2, double delta2, double MassTurnover, double Alpha_star, double Alpha_esc, double Fstar10, double Fesc10, double Mlim_Fstar, double Mlim_Fesc, int method); -double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double M_cond, double sigma2, double delta2, double MassTurnover, +double Nion_ConditionalM_MINI(double growthf, double lnM1, double lnM2, double lnM_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 Xray_ConditionalM(double redshift, double growthf, double lnM1, double lnM2, double lnM_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); diff --git a/src/py21cmfast/src/interp_tables.c b/src/py21cmfast/src/interp_tables.c index d2413935e..7383faf29 100644 --- a/src/py21cmfast/src/interp_tables.c +++ b/src/py21cmfast/src/interp_tables.c @@ -293,6 +293,7 @@ void initialise_Nion_Conditional_spline(float z, float Mcrit_atom, float min_den RGTable2D_f *table_2d, *table_mini; LOG_SUPER_DEBUG("Initialising Nion conditional table at mass %.2e from delta %.2e to %.2e",Mcond,min_density,max_density); + LOG_SUPER_DEBUG("l10Mturns ACG %.2e %.2e MCG %.2e %.2e",log10Mturn_min,log10Mturn_max,log10Mturn_min_MINI,log10Mturn_max_MINI); growthf = dicke(z); double lnMmin = log(Mmin); @@ -351,7 +352,7 @@ void initialise_Nion_Conditional_spline(float z, float Mcrit_atom, float min_den for(i=0;iINTEGRATION_METHOD_ATOMIC)); if(Nion_conditional_table1D.y_arr[i] < -40.) @@ -360,14 +361,14 @@ void initialise_Nion_Conditional_spline(float z, float Mcrit_atom, float min_den continue; } for (j=0; jz_arr[i][j] = log(Nion_ConditionalM(growthf,lnMmin,lnMmax,Mcond,sigma2, + table_2d->z_arr[i][j] = log(Nion_ConditionalM(growthf,lnMmin,lnMmax,log(Mcond),sigma2, overdense_table[i],mturns[j],Alpha_star,Alpha_esc, Fstar10,Fesc10,Mlim_Fstar,Mlim_Fesc,user_params_global->INTEGRATION_METHOD_ATOMIC)); if(table_2d->z_arr[i][j] < -40.) table_2d->z_arr[i][j] = -40.; - table_mini->z_arr[i][j] = log(Nion_ConditionalM_MINI(growthf,lnMmin,lnMmax,Mcond,sigma2,overdense_table[i], + table_mini->z_arr[i][j] = log(Nion_ConditionalM_MINI(growthf,lnMmin,lnMmax,log(Mcond),sigma2,overdense_table[i], mturns_MINI[j],Mcrit_atom,Alpha_star_mini,Alpha_esc,Fstar7_MINI,Fesc7_MINI, Mlim_Fstar_MINI,Mlim_Fesc_MINI,user_params_global->INTEGRATION_METHOD_MINI)); @@ -448,7 +449,7 @@ void initialise_SFRD_Conditional_table(double min_density, double max_density, d #pragma omp for for (i=0; iINTEGRATION_METHOD_ATOMIC)); if(SFRD_conditional_table.y_arr[i] < -50.) @@ -457,7 +458,7 @@ void initialise_SFRD_Conditional_table(double min_density, double max_density, d if(!minihalos) continue; for (k=0; kINTEGRATION_METHOD_MINI)); @@ -483,7 +484,7 @@ 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, +void initialise_Xray_Conditional_table(double min_density, double max_density, double redshift, 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){ @@ -534,7 +535,7 @@ void initialise_Xray_Conditional_table(double min_density, double max_density, d for (i=0; iINTEGRATION_METHOD_ATOMIC)); @@ -545,7 +546,7 @@ void initialise_Xray_Conditional_table(double min_density, double max_density, d } for (k=0; kUSE_INTERPOLATION_TABLES){ return exp(EvaluateRGTable1D_f(delta,&SFRD_conditional_table)); } - return Nion_ConditionalM(growthf,log(M_min),log(M_max),M_cond,sigma_max,delta,Mturn_a, + return Nion_ConditionalM(growthf,log(M_min),log(M_max),log(M_cond),sigma_max,delta,Mturn_a, astro_params_global->ALPHA_STAR,0.,astro_params_global->F_STAR10,1.,Mlim_Fstar,0., user_params_global->INTEGRATION_METHOD_ATOMIC); } @@ -966,7 +967,7 @@ double EvaluateSFRD_Conditional_MINI(double delta, double log10Mturn_m, double g return exp(EvaluateRGTable2D_f(delta,log10Mturn_m,&SFRD_conditional_table_MINI)); } - return Nion_ConditionalM_MINI(growthf,log(M_min),log(M_max),M_cond,sigma_max,delta,pow(10,log10Mturn_m),Mturn_a,astro_params_global->ALPHA_STAR_MINI, + return Nion_ConditionalM_MINI(growthf,log(M_min),log(M_max),log(M_cond),sigma_max,delta,pow(10,log10Mturn_m),Mturn_a,astro_params_global->ALPHA_STAR_MINI, 0.,astro_params_global->F_STAR7_MINI,1.,Mlim_Fstar, 0., user_params_global->INTEGRATION_METHOD_MINI); } @@ -979,7 +980,7 @@ double EvaluateNion_Conditional(double delta, double log10Mturn, double growthf, return exp(EvaluateRGTable1D_f(delta, &Nion_conditional_table1D)); } - return Nion_ConditionalM(growthf,log(M_min),log(M_max),M_cond,sigma_max,delta,pow(10,log10Mturn), + return Nion_ConditionalM(growthf,log(M_min),log(M_max),log(M_cond),sigma_max,delta,pow(10,log10Mturn), astro_params_global->ALPHA_STAR,astro_params_global->ALPHA_ESC,astro_params_global->F_STAR10, astro_params_global->F_ESC10,Mlim_Fstar,Mlim_Fesc,user_params_global->INTEGRATION_METHOD_ATOMIC); } @@ -991,12 +992,12 @@ double EvaluateNion_Conditional_MINI(double delta, double log10Mturn_m, double g return exp(EvaluateRGTable2D_f(delta,log10Mturn_m,table)); } - return Nion_ConditionalM_MINI(growthf,log(M_min),log(M_max),M_cond,sigma_max,delta,pow(10,log10Mturn_m),Mturn_a,astro_params_global->ALPHA_STAR_MINI, + return Nion_ConditionalM_MINI(growthf,log(M_min),log(M_max),log(M_cond),sigma_max,delta,pow(10,log10Mturn_m),Mturn_a,astro_params_global->ALPHA_STAR_MINI, astro_params_global->ALPHA_ESC,astro_params_global->F_STAR7_MINI,astro_params_global->F_ESC7_MINI,Mlim_Fstar, Mlim_Fesc,user_params_global->INTEGRATION_METHOD_MINI); } -double EvaluateXray_Conditional(double delta, double log10Mturn_m, double growthf, double M_min, double M_max, double M_cond, double sigma_max, +double EvaluateXray_Conditional(double delta, double log10Mturn_m, double redshift, 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) @@ -1004,7 +1005,7 @@ double EvaluateXray_Conditional(double delta, double log10Mturn_m, double growth 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), + return Xray_ConditionalM(redshift, growthf, log(M_min), log(M_max), log(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, @@ -1028,13 +1029,13 @@ double EvaluatedFcolldz(double delta, double redshift, double sigma_min, double double EvaluateNhalo(double condition, double growthf, double lnMmin, double lnMmax, double M_cond, double sigma, double delta){ if(user_params_global->USE_INTERPOLATION_TABLES) return EvaluateRGTable1D(condition,&Nhalo_table); - return Nhalo_Conditional(growthf, lnMmin, lnMmax, M_cond, sigma, delta, 0); + return Nhalo_Conditional(growthf, lnMmin, lnMmax, log(M_cond), sigma, delta, 0); } double EvaluateMcoll(double condition, double growthf, double lnMmin, double lnMmax, double M_cond, double sigma, double delta){ if(user_params_global->USE_INTERPOLATION_TABLES) return EvaluateRGTable1D(condition,&Mcoll_table); - return Mcoll_Conditional(growthf, lnMmin, lnMmax, M_cond, sigma, delta, 0); + return Mcoll_Conditional(growthf, lnMmin, lnMmax, log(M_cond), sigma, delta, 0); } //extrapolation function for log-probability based tables diff --git a/src/py21cmfast/src/interp_tables.h b/src/py21cmfast/src/interp_tables.h index ba06dcbe0..5cbea872a 100644 --- a/src/py21cmfast/src/interp_tables.h +++ b/src/py21cmfast/src/interp_tables.h @@ -33,11 +33,11 @@ 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, +void initialise_Xray_Conditional_table(double min_density, double max_density, double redshift, 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 EvaluateXray_Conditional(double delta, double log10Mturn_m, double redshift, 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, diff --git a/src/py21cmfast/wrapper/cfuncs.py b/src/py21cmfast/wrapper/cfuncs.py index de703eab6..ce2c740f8 100644 --- a/src/py21cmfast/wrapper/cfuncs.py +++ b/src/py21cmfast/wrapper/cfuncs.py @@ -936,7 +936,7 @@ def evaluate_SFRD_cond( growthf, np.log(M_min), np.log(M_max), - cond_mass, + np.log(cond_mass), sigma_cond, densities, mcrit_atom, @@ -953,7 +953,7 @@ def evaluate_SFRD_cond( growthf, np.log(M_min), np.log(M_max), - cond_mass, + np.log(cond_mass), sigma_cond, densities[:, None], 10 ** l10mturns[None, :], @@ -1083,7 +1083,7 @@ def evaluate_Nion_cond( growthf, np.log(M_min), np.log(M_max), - cond_mass, + np.log(cond_mass), sigma_cond, densities[:, None] if flag_options.USE_MINI_HALOS else densities, 10 ** l10mturns[None, :] if flag_options.USE_MINI_HALOS else mcrit_atom, @@ -1100,7 +1100,7 @@ def evaluate_Nion_cond( growthf, np.log(M_min), np.log(M_max), - cond_mass, + np.log(cond_mass), sigma_cond, densities[:, None], 10 ** l10mturns[None, :],