diff --git a/src/py21cmfast/plotting.py b/src/py21cmfast/plotting.py index 08f26bb5..4eae3517 100644 --- a/src/py21cmfast/plotting.py +++ b/src/py21cmfast/plotting.py @@ -99,7 +99,7 @@ def _imshow_slice( if not rotate: slc = slc.T - if cmap == "EoR": + if cmap == "EoR" and "norm" not in imshow_kw: imshow_kw["norm"] = colors.Normalize(vmin=-150, vmax=30) norm_kw = {k: imshow_kw.pop(k) for k in ["vmin", "vmax"] if k in imshow_kw} diff --git a/src/py21cmfast/src/HaloField.c b/src/py21cmfast/src/HaloField.c index daf2587f..ee482c2e 100644 --- a/src/py21cmfast/src/HaloField.c +++ b/src/py21cmfast/src/HaloField.c @@ -180,27 +180,6 @@ int ComputeHaloField(float redshift_desc, float redshift, UserParams *user_param // Sheth-Tormen mass function (as of right now, We do not even reproduce EPS results) delta_crit = growth_factor*sheth_delc_dexm(Deltac/growth_factor, sigma_z0(M)); - // if(global_params.DELTA_CRIT_MODE == 1){ - // //This algorithm does not use the sheth tormen OR Jenkins parameters, - // // rather it uses a reduced barrier to correct for some part of this algorithm, - // // which would otherwise result in ~3x halos at redshift 6 (including with EPS). - // // Once I Figure out the cause of the discrepancy I can adjust again but for now - // // we will use the parameters that roughly give the ST mass function. - // // delta_crit = growth_factor*sheth_delc(Deltac/growth_factor, sigma_z0(M)); - // if(user_params->HMF==1) { - // // use sheth tormen correction - // delta_crit = growth_factor*sheth_delc(Deltac/growth_factor, sigma_z0(M)); - // } - // else if(user_params->HMF==4) { - // // use Delos 2023 flat barrier - // delta_crit = DELTAC_DELOS; - // } - // else if(user_params->HMF!=0){ - // LOG_WARNING("Halo Finder: You have selected DELTA_CRIT_MODE==1 with HMF %d which does not have a barrier - // , using EPS deltacrit = 1.68",user_params->HMF); - // } - // } - // first let's check if virialized halos of this size are rare enough // that we don't have to worry about them (let's define 7 sigma away, as in Mesinger et al 05) if ((sigma_z0(M)*growth_factor*7.) < delta_crit){ diff --git a/src/py21cmfast/src/InputParameters.c b/src/py21cmfast/src/InputParameters.c index 935550e1..e2a5a42b 100644 --- a/src/py21cmfast/src/InputParameters.c +++ b/src/py21cmfast/src/InputParameters.c @@ -18,30 +18,35 @@ CosmoParams *cosmo_params_global; AstroParams *astro_params_global; FlagOptions *flag_options_global; +// These need to be removed, I am commenting below roughly where/when +// they are used, and my tentative decision on what to do with them GlobalParams global_params = { - .ALPHA_UVB = 5.0, - .EVOLVE_DENSITY_LINEARLY = 0, - .SMOOTH_EVOLVED_DENSITY_FIELD = 0, - .R_smooth_density = 0.2, - .HII_ROUND_ERR = 1e-5, - .FIND_BUBBLE_ALGORITHM = 2, - .N_POISSON = 5, - .T_USE_VELOCITIES = 1, - .MAX_DVDR = 0.2, - .DELTA_R_HII_FACTOR = 1.1, - .DELTA_R_FACTOR = 1.1, - .HII_FILTER = 0, - .INITIAL_REDSHIFT = 300., - .R_OVERLAP_FACTOR = 1., - .DELTA_CRIT_MODE = 1, - .HALO_FILTER = 0, - .OPTIMIZE = 0, - .OPTIMIZE_MIN_MASS = 1e11, - - - .CRIT_DENS_TRANSITION = 1.2, - .MIN_DENSITY_LOW_LIMIT = 9e-8, + .ALPHA_UVB = 5.0, // IonisationBox, for Gamma (move to AstroParams) + .EVOLVE_DENSITY_LINEARLY = 0, //PerturnbedField, uses linear growth instead of 1/2LPT (move to UserParams) + .SMOOTH_EVOLVED_DENSITY_FIELD = 0, //PerturbedField, smooth field after perturb (delete) + .R_smooth_density = 0.2, //With above, radius for smoothing (delete) + .HII_ROUND_ERR = 1e-5, // IonisationBox, Nion threshold to run the ES (move to define) + .FIND_BUBBLE_ALGORITHM = 2, // IonisationBox central pixel vs sphere (move to UserParams) + .MAX_DVDR = 0.2, //BrightnessTemp, maximum velocity gradient for RSDS (move to AstroParams?) + .DELTA_R_HII_FACTOR = 1.1, //IonisationBox, radius factor for reion ES (move to AstroParams) + .DELTA_R_FACTOR = 1.1, //HaloField, radius factor for DexM halo ES (move to UserParams) + .HII_FILTER = 0, //IonisationBox, filter for reion ES (move to AstroParams) + .HALO_FILTER = 0, //HaloField, filter for halo ES (move to UserParams) + .HEAT_FILTER = 0, //SpinTemperature, filter for Ts model (move to AstroParams) + .FILTER = 0, //Filter for calculating sigma(M) and R <---> M conversions Unlike others this is 0: Tophat, 1: Gaussian + + .INITIAL_REDSHIFT = 300., //PerturbField, initial redshift for moving the particles (move to UserParams) + + .OPTIMIZE = 0, //HaloField, excluding zones around existing halos for speed (delete?) + .R_OVERLAP_FACTOR = 1., //HaloField, exclude cells within X*radius for optimisation (delete?) + .OPTIMIZE_MIN_MASS = 1e11, //don't optimize below this mass (delete?) + + .CRIT_DENS_TRANSITION = 1.2, //Gauss-Legendre does poorly at high delta, switch to GSL-QAG here (define?) + .MIN_DENSITY_LOW_LIMIT = 9e-8, //lower limit for extrapolated density in ionised temperatures (define) + + //all the z-photoncons parameters + //TODO: walk through the model again, it's confusing .RecombPhotonCons = 0, .PhotonConsStart = 0.995, .PhotonConsEnd = 0.3, @@ -49,41 +54,41 @@ GlobalParams global_params = { .PhotonConsEndCalibz = 3.5, .PhotonConsSmoothing = 1, - .HEAT_FILTER = 0, - .CLUMPING_FACTOR = 2., - .R_XLy_MAX = 500., - .NUM_FILTER_STEPS_FOR_Ts = 40, - .TK_at_Z_HEAT_MAX = -1, - .XION_at_Z_HEAT_MAX = -1, - .Pop = 2, - .Pop2_ion = 5000, - .Pop3_ion = 44021, + .CLUMPING_FACTOR = 2., //SpinTemperature, for recombinations in neutral IGM + .R_XLy_MAX = 500., //SpinTemperature, max radius for ES (Move to AstroParams) + .NUM_FILTER_STEPS_FOR_Ts = 40, //SpinTemperature, self-explanatory (move to AstroParams) - .NU_X_BAND_MAX = 2000.0, - .NU_X_MAX = 10000.0, + .Pop = 2, //Only used in lyman alpha spectral calculations, when minihalos are off (delete) + .Pop2_ion = 5000, //Move to AstroParams + .Pop3_ion = 44021, // Move to AstroParams - .NBINS_LF = 100, + .NU_X_BAND_MAX = 2000.0, //Used in SFR -> Lx conversion factor (define) + .NU_X_MAX = 10000.0, //Limit for nu integrals (define) - .P_CUTOFF = 0, - .M_WDM = 2, - .g_x = 1.5, + .P_CUTOFF = 0, //Used in EH power spectrum and with WDM in setting minimum mass, seems inconsistent (delete?) + .M_WDM = 2, //Used in setting minimum mass if P_CUTOFF is True. (delete?) + .g_x = 1.5, //Used in setting minimum mass if P_CUTOFF is True. (delete?) + + //Move These to CosmoParams .OMn = 0.0, .OMk = 0.0, .OMr = 8.6e-5, .OMtot = 1.0, .Y_He = 0.245, .wl = -1.0, - .SHETH_b = 0.15, //In the literature this is 0.485 (RM08) or 0.5 (SMT01) or 0.34 (Barkana+01) Master 21cmFAST currently has 0.15 - .SHETH_c = 0.05, //In the literature this is 0.615 (RM08) or 0.6 (SMT01) or 0.81 (Barkana+01) Master 21cmFAST currently has 0.05 - .Zreion_HeII = 3.0, - .FILTER = 0, - .R_BUBBLE_MIN = 0.620350491, - .M_MIN_INTEGRAL = 1e5, - .M_MAX_INTEGRAL = 1e16, - .T_RE = 2e4, + //Adjustments to SMT Barrier so that DexM fits the mass function (Move to UserParams?) + .SHETH_b = 0.15, //In the literature this is 0.485 (RM08) or 0.5 (SMT01) or 0.34 (Barkana+01) + .SHETH_c = 0.05, //In the literature this is 0.615 (RM08) or 0.6 (SMT01) or 0.81 (Barkana+01) + + .Zreion_HeII = 3.0, //used in tau computation, move to tau arguments + .R_BUBBLE_MIN = 0.620350491, //minimum bubble radius in cMpc, used as a ceiling for the last radius (define?) + .M_MIN_INTEGRAL = 1e5, //minimum mass for the integrals of the mass function (define?) + .M_MAX_INTEGRAL = 1e16, //maximum mass for the integrals of the mass function (define?) + + .T_RE = 2e4, //reionisation temperature, move to AstroParams - .VAVG=25.86, + .VAVG=25.86, //Move to AstroParams, when FlagOptions.FIX_VCB_AVG is true, this is the value - .USE_ADIABATIC_FLUCTUATIONS = 1, + .USE_ADIABATIC_FLUCTUATIONS = 1, //For first Ts Box, (fix to True) }; diff --git a/src/py21cmfast/src/Stochasticity.c b/src/py21cmfast/src/Stochasticity.c index 2441f37d..0c3a8527 100644 --- a/src/py21cmfast/src/Stochasticity.c +++ b/src/py21cmfast/src/Stochasticity.c @@ -460,8 +460,8 @@ int stoc_partition_sample(struct HaloSamplingConstants * hs_constants, gsl_rng * double nu_fudge_factor = user_params_global->HALOMASS_CORRECTION; //set initial amount - // M_remaining = M_cond; // full condition - M_remaining = exp_M; //subtract unresolved mass + M_remaining = M_cond; // full condition + // M_remaining = exp_M; //subtract unresolved mass lnM_remaining = log(M_remaining); double nu_min; diff --git a/src/py21cmfast/src/_inputparams_wrapper.h b/src/py21cmfast/src/_inputparams_wrapper.h index b9fa648a..3d7f9413 100644 --- a/src/py21cmfast/src/_inputparams_wrapper.h +++ b/src/py21cmfast/src/_inputparams_wrapper.h @@ -128,15 +128,12 @@ typedef struct GlobalParams{ float R_smooth_density; float HII_ROUND_ERR; int FIND_BUBBLE_ALGORITHM; - int N_POISSON; - int T_USE_VELOCITIES; float MAX_DVDR; float DELTA_R_HII_FACTOR; float DELTA_R_FACTOR; int HII_FILTER; float INITIAL_REDSHIFT; float R_OVERLAP_FACTOR; - int DELTA_CRIT_MODE; int HALO_FILTER; int OPTIMIZE; float OPTIMIZE_MIN_MASS; @@ -156,8 +153,7 @@ typedef struct GlobalParams{ double CLUMPING_FACTOR; float R_XLy_MAX; int NUM_FILTER_STEPS_FOR_Ts; - double TK_at_Z_HEAT_MAX; - double XION_at_Z_HEAT_MAX; + int Pop; float Pop2_ion; float Pop3_ion; @@ -165,8 +161,6 @@ typedef struct GlobalParams{ float NU_X_BAND_MAX; float NU_X_MAX; - int NBINS_LF; - int P_CUTOFF; float M_WDM; float g_x;