21cmfast / 21cmFAST

Official repository for 21cmFAST: a code for generating fast simulations of the cosmological 21cm signal
MIT License
58 stars 37 forks source link

Fixing some issues in the HaloBox struct and Mass Function Integrals. #379

Closed daviesje closed 5 months ago

daviesje commented 6 months ago

This solves issue (#371) by improving the hardcoded double power law in the stellar to halo mass relation stellar mass, so that it uses a flag option USE_UPPER_STELLAR_TURNOVER to switch it on or off, and two global params UPPER_STELLAR_TURNOVER_MASS and UPPER_STELLAR_TURNOVER_INDEX to define the power law segment. This is currently only implmented in the halo model since it will only effect galaxy-based observables (the effect on the radiative backgrounds is minimal).

I have also fixed some issues with the function set_fixed_grids, a function used to estimate the average quantities of halo properties in each cell. This is used either when we do not want stochastic halo samples but we still want to use the halo-based radiative backgrounds (USE_HALO_FIELD=True FIXED_HALO_GRIDS=True), or when we want to add a population of halos below the minimum sampler mass (AVG_BELOW_SAMPLER=True). Since these mass function integrals are now used over a wider variety of upper/lower mass limits we now have to be more careful about how we set limits around the critical density, and so here I assume that if the cell is above critical density, the entire cell mass is within a single halo with mass exactly equal to the cell's Lagrangian mass (without (1+delta) since that is how the CMFs are defined).

This does involve a rather clumsy way of checking the mass limits (e.g line 1555 in ps.c if(M_cond*(1-FRACT_FLOAT_ERR) <= exp(lnM2))) since we only want to return the halo if it is within the mass integral. If the upper limit of the integral is defined exactly at the condition mass, floating point errors can result from taking the exponent and we end up on the wrong side of the conditional if we don't include the tolerance. If you think it's necessary, I can create a more robust way of checking the mass limit but I think this works well enough.