Skip to content

Commit

Permalink
WIP fixing xray integrals
Browse files Browse the repository at this point in the history
  • Loading branch information
@daviesje
daviesje committed Dec 5, 2024
1 parent 2bcfeb3 commit 7bd9a18
Show file tree
Hide file tree
Showing 6 changed files with 134 additions and 116 deletions.
31 changes: 22 additions & 9 deletions src/py21cmfast/src/HaloBox.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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);
}

Expand Down Expand Up @@ -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]);
}
}
}
}
Expand Down
129 changes: 84 additions & 45 deletions src/py21cmfast/src/hmf.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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));
Expand Down Expand Up @@ -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);
Expand All @@ -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){
Expand Down Expand Up @@ -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);
}

Expand Down Expand Up @@ -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,
Expand All @@ -817,20 +855,20 @@ 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.;
}

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,
Expand All @@ -839,19 +877,19 @@ 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.;
}
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 = {
Expand All @@ -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,&params); //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,&params) / exp(lnM_cond); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta)
else
return 0.;
}
Expand All @@ -890,30 +928,30 @@ 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 = {
.growthf = growthf,
.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,&params); //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,&params) / exp(lnM_cond); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta)
else
return 0.;
}
Expand All @@ -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,
Expand All @@ -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,&params); //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,&params) / exp(lnM_cond); //NOTE: condition mass is used as if it were Lagrangian (no 1+delta)
else
return 0.;
}
Expand Down
Loading

0 comments on commit 7bd9a18

Please sign in to comment.