diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index f3fc8974c..47912f904 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -9,6 +9,7 @@ repos: hooks: - id: trailing-whitespace - id: check-added-large-files + args: ['--maxkb=2000'] - id: check-ast - id: check-json - id: check-merge-conflict diff --git a/docs/acknowledge.rst b/docs/acknowledge.rst index e81d3ad92..ee99c6457 100644 --- a/docs/acknowledge.rst +++ b/docs/acknowledge.rst @@ -15,11 +15,24 @@ If you use ``21cmFAST v3+`` in your research please cite both of: In addition, the following papers introduce various features into ``21cmFAST``. If you use these features, please cite the relevant papers. +Lyman alpha multiple scattering + + Flitter, J., Munoz, J. B., Mesinger, A., + "Semi-analytical approach to Lyman-alpha multiple-scattering in 21-cm signal simulations", + eprint arXiv:2601.14360, 2026. https://doi.org/10.48550/arXiv.2601.14360. + Discrete Halo Sampler / version 4: Davies, J. E., Mesinger, A., Murray, S. G., "Efficient simulation of discrete galaxy populations and associated radiation fields during the first billion years", - eprint arXiv:2504.17254, 2025. https://doi.org/10.48550/arXiv.2504.17254 + Astronomy and Astrophysics, vol. 701, A. 236, 2025. https://doi.org/10.1051/0004-6361/202554951. + +Photon non-conservation correction: + + Park, J., Greig, B., Mesinger, A., + "Calibrating excursion set reionization models to approximately conserve ionizing photons", + Monthly Notices of the Royal Astronomical Society, vol. 517, no. 1, + pp 192-200, 2022 https://doi.org/10.1093/mnras/stac2756. Mini-halos: @@ -42,6 +55,20 @@ Mass-dependent ionizing efficiency: functions and the 21-cm signal”, Monthly Notices of the Royal Astronomical Society, vol. 484, no. 1, pp. 933–949, 2019. https://doi.org/10.1093/mnras/stz032. +Lightcone and redshift space distortions: + + Greig, B., Mesinger, A., + "21CMMC with a 3D light-cone: the impact of the co-evolution approximation on the astrophysics of reionisation and cosmic dawn.", + Monthly Notices of the Royal Astronomical Society, vol. 477, no. 3, + pp 3217-3229, 2018. https://doi.org/10.1093/mnras/sty796. + +Inhomogeneous recombination: + + Sobacchi, E., Mesinger, A., + "Inhomogeneous recombinations during cosmic reionization", + Monthly Notices of the Royal Astronomical Society, vol. 440, no. 2, + pp 1662-1673, 2014. https://doi.org/10.1093/mnras/stu377. + If you are unsure which modules are used within your simulations, we provide a handy function to print out which works to refer ``py21cmfast.utils.show_references``, which accepts a single instance of the ``InputParameters`` class and shows which papers are relevant for your simulation. diff --git a/docs/conf.py b/docs/conf.py index d5e89b760..3ca470299 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -84,7 +84,7 @@ def check_module(self): 'undoc-members', 'show-inheritance', 'show-module-summary', - 'special-members', + #'special-members', 'imported-members', 'inherited-members', ] @@ -94,7 +94,8 @@ def check_module(self): autoapi_member_order = 'groupwise' autoapi_own_page_level = 'class' autoapi_keep_files = True -autodoc_typehints = 'description' +autodoc_typehints = 'both' +autoapi_template_dir = 'templates/autoapi' autosummary_generate = False numpydoc_show_class_members = False @@ -159,11 +160,3 @@ def check_module(self): "templates", "**.ipynb_checkpoints", ] - -def setup(app): - """Setup function for Sphinx.""" - print("SETTING UP THE CLASS BASED DECORATOR DOCUMENTER") - app.add_autodocumenter(ClassDecoratedDocumenter) - return { - 'parallel_read_safe': True - } \ No newline at end of file diff --git a/docs/images/papers/multiple_scattering/coevals_z11.png b/docs/images/papers/multiple_scattering/coevals_z11.png new file mode 100644 index 000000000..cdfc6c1f7 Binary files /dev/null and b/docs/images/papers/multiple_scattering/coevals_z11.png differ diff --git a/docs/images/papers/multiple_scattering/coevals_z20.png b/docs/images/papers/multiple_scattering/coevals_z20.png new file mode 100644 index 000000000..f16c6183e Binary files /dev/null and b/docs/images/papers/multiple_scattering/coevals_z20.png differ diff --git a/docs/index.rst b/docs/index.rst index b77ad26f0..d1a84b500 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -8,6 +8,7 @@ quickstart tutorials performance + models faqs/index updates_from_v3 acknowledge @@ -15,6 +16,8 @@ .. toctree:: :caption: API Reference :hidden: + :maxdepth: 1 + autoapi/py21cmfast/index .. toctree:: diff --git a/docs/models.rst b/docs/models.rst new file mode 100644 index 000000000..45741a15d --- /dev/null +++ b/docs/models.rst @@ -0,0 +1,84 @@ +Description of Physical Models +============================== + +``21cmFAST`` supports many physical models and effects that can often be toggled on or +off interchangeably (though some depend on others). The +:class:`~py21cmfast.MatterOptions` and +:class:`~py21cmfast.AstroOptions` classes contain all of the flags +that control which models to include in the simulation. + + +Below we provide a brief explanation on how some of the flags modify the output of +``21cmFAST``. + +Models of the Matter Field +-------------------------- + +The parameters and flags of :class:`~py21cmfast.MatterOptions` +are used to control how cosmological matter fields (e.g. densities, velocities and halo +properties) are evaluated in the simulation. + +Source Model +~~~~~~~~~~~~ + +.. note:: Set with the ``SOURCE_MODEL`` parameter of :class:`~py21cmfast.MatterOptions`. + +To be filled. + +Astrophysical Models +-------------------- + +The parameters and flags of :class:`~py21cmfast.AstroOptions` +are used to control how astrophysical quantities (e.g. star formation rate, UV and +ionizing radiation) are evaluated in the simulation. + +It is important to stress that the generation of +:class:`~py21cmfast.PerturbedField` and +:class:`~py21cmfast.HaloCatalog` objects do not depend on these +parameters (nor on :class:`~py21cmfast.AstroParams`). Therefore, if +the cache contains :class:`~py21cmfast.PerturbedField` and +:class:`~py21cmfast.HaloCatalog` objects that had been previously +generated with a different set of :class:`~py21cmfast.AstroOptions`, these objects will +be loaded from the cache, instead of being re-evaluated. +This architecture allows one to quickly simulate different astrophysical models, given a +cosmological model. + +Spin Temperature Fluctuations +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +To be filled. + +Minihalos +~~~~~~~~~ + +To be filled. + +Multiple Scattering of Lyman Alpha Photons +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. note:: Toggled with the ``LYA_MULTIPLE_SCATTERING`` flag of + :class:`~py21cmfast.AstroOptions` (default: ``False``). + +The physical effect that enables the absorption feature in the 21-cm signal during +cosmic dawn is the strong coupling between the spin temperature and the gas kinetic +temperature. This coupling is obtained through Lyman alpha radiation that comes from the +first stars and is absorbed by the IGM. Since the cross section for the interaction +between photons near the Lyman alpha resonance frequency and HI atoms in the IGM has a +non-negligible width, the Lyman alpha photons that are absorbed by the IGM had not +traveled in straight lines, but rather had scattered along their path. This means that +the effective Lyman alpha emissivity that the IGM "sees" becomes more local, thereby +increasing the contrast in ``J_alpha`` (Lyman alpha flux) maps in the simulation, as can +be seen below. In the context of the 21-cm signal, this effect becomes negligible at +sufficiently low redshifts (normally below :math:`z<15`), once the spin temperature +completely follows gas kinetic temperature. For more information on this effect in the +simulation, see `Flitter, Munoz and Mesinger 2026 `_. + +.. image:: ./images/papers/multiple_scattering/coevals_z20.png + :width: 800px + :align: center + :alt: the effect of the Lyman alpha multiple scattering effect in the simulation, at z=20. Figure is taken from Flitter, Munoz and Mesinger 2026 (arxiv: 2601.14360). + +.. image:: ./images/papers/multiple_scattering/coevals_z11.png + :width: 800px + :align: center + :alt: the effect of the Lyman alpha multiple scattering effect in the simulation, at z=20. Figure is taken from Flitter, Munoz and Mesinger 2026 (arxiv: 2601.14360). diff --git a/docs/updates_from_v3.rst b/docs/updates_from_v3.rst index 72fb3b111..9a08154fe 100644 --- a/docs/updates_from_v3.rst +++ b/docs/updates_from_v3.rst @@ -1,6 +1,6 @@ -============================== -Changes from 21cmFAST v3 to v4 -============================== +===================== +Changes from v3 to v4 +===================== This page outlines the main differences between the usage of ``21cmFAST`` versions 3 and 4. Use this documentation as a reference for understanding and adapting to the updates in diff --git a/setup.py b/setup.py index 462968819..7dd103ac7 100644 --- a/setup.py +++ b/setup.py @@ -24,7 +24,7 @@ def _read(name: str): test_req = [ "clang-format", "clang-tidy", - "hmf", + "mpmath", "pre-commit", "pytest>=5.0", "pytest-cov", diff --git a/src/py21cmfast/__init__.py b/src/py21cmfast/__init__.py index af1be6d5a..9c1c606f6 100644 --- a/src/py21cmfast/__init__.py +++ b/src/py21cmfast/__init__.py @@ -12,7 +12,7 @@ __version__ = "unknown" __all__ = [ - "DATA_PATH", + "_DATA_PATH", "AngularLightconer", "AstroOptions", "AstroParams", @@ -73,7 +73,7 @@ from . import lightconers, plotting, wrapper from ._cfg import config -from ._data import DATA_PATH +from ._data import _DATA_PATH from ._logging import configure_logging from ._templates import create_params_from_template, list_templates, write_template from .drivers.coeval import Coeval, generate_coeval, run_coeval diff --git a/src/py21cmfast/_cfg.py b/src/py21cmfast/_cfg.py index 302b8e323..47d2d2564 100644 --- a/src/py21cmfast/_cfg.py +++ b/src/py21cmfast/_cfg.py @@ -7,7 +7,7 @@ from typing import ClassVar from . import yaml -from ._data import DATA_PATH +from ._data import _DATA_PATH from .c_21cmfast import ffi, lib from .wrapper.structs import StructInstanceWrapper @@ -23,7 +23,7 @@ class Config(dict): _defaults: ClassVar = { "direc": "~/21cmFAST-cache", "ignore_R_BUBBLE_MAX_error": False, - "external_table_path": DATA_PATH, + "external_table_path": _DATA_PATH, "HALO_CATALOG_MEM_FACTOR": 1.2, "EXTRA_HALOBOX_FIELDS": False, "safe_read": True, @@ -106,3 +106,4 @@ def load(cls, file_name: str | Path): # On import, load the default config config = Config() +"""A singleton Config object used to manage global configuration options.""" diff --git a/src/py21cmfast/_data/__init__.py b/src/py21cmfast/_data/__init__.py index ae4da01b0..52b916822 100644 --- a/src/py21cmfast/_data/__init__.py +++ b/src/py21cmfast/_data/__init__.py @@ -1,3 +1,3 @@ from pathlib import Path -DATA_PATH = Path(__file__).parent +_DATA_PATH = Path(__file__).parent diff --git a/src/py21cmfast/drivers/coeval.py b/src/py21cmfast/drivers/coeval.py index 31186a0c0..9a9e1bcb5 100644 --- a/src/py21cmfast/drivers/coeval.py +++ b/src/py21cmfast/drivers/coeval.py @@ -756,6 +756,7 @@ def _redshift_loop_generator( this_xraysource = sf.compute_xray_source_field( redshift=z, hboxes=[*hbox_arr, this_halobox], + previous_ionize_box=getattr(prev_coeval, "ionized_box", None), write=write.xray_source_box, **kw, ) diff --git a/src/py21cmfast/drivers/single_field.py b/src/py21cmfast/drivers/single_field.py index 29b6755d3..1c7f443a2 100644 --- a/src/py21cmfast/drivers/single_field.py +++ b/src/py21cmfast/drivers/single_field.py @@ -9,6 +9,7 @@ import warnings import numpy as np +from astropy import constants from astropy import units as un from astropy.cosmology import z_at_value @@ -462,6 +463,7 @@ def compute_xray_source_field( initial_conditions: InitialConditions, hboxes: list[HaloBox], redshift: float, + previous_ionize_box: IonizedBox | None = None, ) -> XraySourceBox: r""" Compute filtered grid of SFR for use in spin temperature calculation. @@ -478,6 +480,9 @@ def compute_xray_source_field( The initial conditions of the run. The user and cosmo params hboxes: Sequence of :class:`~HaloBox` instances This contains the list of Halobox instances which are used to create this source field + previous_ionize_box: :class:`IonizedBox` or None + An ionized box at higher redshift. This is only used if `LYA_MULTIPLE_SCATTERING` is true. + Returns ------- @@ -528,6 +533,34 @@ def compute_xray_source_field( # inner and outer redshifts (following the C code) zpp_avg = zpp_edges - np.diff(np.insert(zpp_edges, 0, redshift)) / 2 + # Compute the comoving diffusion scale in the case of Lyman alpha multiple scattering + if inputs.astro_options.LYA_MULTIPLE_SCATTERING: + # TODO: In principle, the diffusion scale varies locally but for simplicty, we consider the global ionization value. + # See https://github.com/21cmfast/21cmFAST/issues/606. + if previous_ionize_box is None: + x_HI = 1.0 + else: + x_HI = previous_ionize_box.neutral_fraction.value.mean() + A_alpha = 6.25e8 * un.Hz + nu_Lya = 2.46606727e15 * un.Hz + n_H_z0 = ( + (1.0 - inputs.cosmo_params.Y_He) + * inputs.cosmo_params.cosmo.critical_density(0) + * inputs.cosmo_params.OMb + / constants.m_p + ) + # Eq. (24) in arxiv: 2601.14360 + R_star = 3.0 * constants.c**4 * A_alpha**2 * n_H_z0 * x_HI * (1.0 + redshift) + R_star /= ( + 32.0 + * np.pi**3 + * nu_Lya**4 + * inputs.cosmo_params.cosmo.H0**2 + * inputs.cosmo_params.OMm + ) + else: + R_star = 0.0 * un.Mpc + interp_fields = ["halo_sfr", "halo_xray"] if inputs.astro_options.USE_MINI_HALOS: interp_fields += ["halo_sfr_mini", "log10_Mcrit_MCG_ave"] @@ -544,6 +577,9 @@ def compute_xray_source_field( if inputs.astro_options.USE_MINI_HALOS: box.filtered_sfr_mini.value[i] = 0 box.mean_log10_Mcrit_LW.value[i] = inputs.astro_params.M_TURN # minimum + if inputs.astro_options.LYA_MULTIPLE_SCATTERING: + box.filtered_sfr_lw.value[i] = 0 + box.filtered_sfr_mini_lw.value[i] = 0 logger.debug(f"ignoring Radius {i} which is above Z_HEAT_MAX") continue @@ -563,6 +599,9 @@ def compute_xray_source_field( if inputs.astro_options.USE_MINI_HALOS: box.filtered_sfr_mini.value[i] = 0 box.mean_log10_Mcrit_LW.value[i] = hbox_interp.log10_Mcrit_MCG_ave + if inputs.astro_options.LYA_MULTIPLE_SCATTERING: + box.filtered_sfr_lw.value[i] = 0 + box.filtered_sfr_mini_lw.value[i] = 0 logger.debug(f"ignoring Radius {i} due to no stars") continue @@ -571,6 +610,7 @@ def compute_xray_source_field( R_inner=R_inner, R_outer=R_outer, R_ct=i, + R_star=R_star.to("Mpc").value, allow_already_computed=True, ) diff --git a/src/py21cmfast/io/h5.py b/src/py21cmfast/io/h5.py index 7cd1b39c3..695c66f50 100644 --- a/src/py21cmfast/io/h5.py +++ b/src/py21cmfast/io/h5.py @@ -205,7 +205,7 @@ def write_outputs_to_group( new = array.written_to_disk(H5Backend(group.file.filename, f"{group.name}/{k}")) setattr(output, k, new) - for k in output.struct.primitive_fields: + for k in output._struct.primitive_fields: try: group.attrs[k] = getattr(output, k) except TypeError as e: diff --git a/src/py21cmfast/plotting.py b/src/py21cmfast/plotting.py index 62dafd912..deee7e4ef 100644 --- a/src/py21cmfast/plotting.py +++ b/src/py21cmfast/plotting.py @@ -170,7 +170,7 @@ def coeval_sliceplot( """ if kind is None: if isinstance(struct, outputs.OutputStruct): - kind = struct.struct.fieldnames[0] + kind = struct._struct.fieldnames[0] elif isinstance(struct, Coeval): kind = "brightness_temp" diff --git a/src/py21cmfast/src/HaloCatalog.c b/src/py21cmfast/src/HaloCatalog.c index 052365dc8..ab7490be9 100644 --- a/src/py21cmfast/src/HaloCatalog.c +++ b/src/py21cmfast/src/HaloCatalog.c @@ -207,7 +207,7 @@ int ComputeHaloCatalog(float redshift_desc, float redshift, InitialConditions *b // now filter the box on scale R // 0 = top hat in real space, 1 = top hat in k space - filter_box(density_field, box_dim, matter_options_global->HALO_FILTER, R, 0.); + filter_box(density_field, box_dim, matter_options_global->HALO_FILTER, R, 0., 0.); // do the FFT to get delta_m box dft_c2r_cube(matter_options_global->USE_FFTW_WISDOM, grid_dim, z_dim, @@ -369,7 +369,7 @@ int ComputeHaloCatalog(float redshift_desc, float redshift, InitialConditions *b filter_box(density_field, box_dim, 0, physconst.l_factor * simulation_options_global->BOX_LEN / (simulation_options_global->HII_DIM + 0.0), - 0.); + 0., 0.); } dft_c2r_cube(matter_options_global->USE_FFTW_WISDOM, simulation_options_global->DIM, D_PARA, simulation_options_global->N_THREADS, density_field); diff --git a/src/py21cmfast/src/InitialConditions.c b/src/py21cmfast/src/InitialConditions.c index ce372ae1e..e50df8a7d 100644 --- a/src/py21cmfast/src/InitialConditions.c +++ b/src/py21cmfast/src/InitialConditions.c @@ -198,7 +198,7 @@ void compute_relative_velocities(fftwf_complex *box, fftwf_complex *box_saved, f filter_box(box, hi_dim, 0, physconst.l_factor * simulation_options_global->BOX_LEN / (simulation_options_global->HII_DIM + 0.0), - 0.); + 0., 0.); } // fft each velocity component back to real space @@ -332,7 +332,7 @@ void compute_velocity_fields(fftwf_complex *box, fftwf_complex *box_saved, float filter_box(box, hi_dim, 0, physconst.l_factor * simulation_options_global->BOX_LEN / (simulation_options_global->HII_DIM + 0.0), - 0.); + 0., 0.); } } @@ -515,7 +515,7 @@ void compute_velocity_fields_2LPT(fftwf_complex *box, fftwf_complex *box_saved, filter_box(box, hi_dim, 0, physconst.l_factor * simulation_options_global->BOX_LEN / (simulation_options_global->HII_DIM + 0.0), - 0.); + 0., 0.); } } @@ -703,7 +703,7 @@ int ComputeInitialConditions(unsigned long long random_seed, InitialConditions * filter_box(HIRES_box, hi_dim, 0, physconst.l_factor * simulation_options_global->BOX_LEN / (simulation_options_global->HII_DIM + 0.0), - 0.); + 0., 0.); } // FFT back to real space diff --git a/src/py21cmfast/src/IonisationBox.c b/src/py21cmfast/src/IonisationBox.c index 29297767a..da98f5c52 100644 --- a/src/py21cmfast/src/IonisationBox.c +++ b/src/py21cmfast/src/IonisationBox.c @@ -599,25 +599,27 @@ void copy_filter_transform(struct FilteredGrids *fg_struct, struct IonBoxConstan if (rspec.R_index > 0) { double R = rspec.R; - filter_box(fg_struct->deltax_filtered, box_dim, consts->hii_filter, R, 0.); + filter_box(fg_struct->deltax_filtered, box_dim, consts->hii_filter, R, 0., 0.); if (astro_options_global->USE_TS_FLUCT) { - filter_box(fg_struct->xe_filtered, box_dim, consts->hii_filter, R, 0.); + filter_box(fg_struct->xe_filtered, box_dim, consts->hii_filter, R, 0., 0.); } if (consts->filter_recombinations) { - filter_box(fg_struct->N_rec_filtered, box_dim, consts->hii_filter, R, 0.); + filter_box(fg_struct->N_rec_filtered, box_dim, consts->hii_filter, R, 0., 0.); } if (consts->lagrangian_source_grids) { int filter_hf = astro_options_global->USE_EXP_FILTER ? 3 : consts->hii_filter; - filter_box(fg_struct->stars_filtered, box_dim, filter_hf, R, consts->mfp_meandens); + filter_box(fg_struct->stars_filtered, box_dim, filter_hf, R, consts->mfp_meandens, 0.); if (astro_options_global->INHOMO_RECO) { - filter_box(fg_struct->sfr_filtered, box_dim, filter_hf, R, consts->mfp_meandens); + filter_box(fg_struct->sfr_filtered, box_dim, filter_hf, R, consts->mfp_meandens, + 0.); } } else { if (astro_options_global->USE_MINI_HALOS) { - filter_box(fg_struct->prev_deltax_filtered, box_dim, consts->hii_filter, R, 0.); + filter_box(fg_struct->prev_deltax_filtered, box_dim, consts->hii_filter, R, 0., 0.); filter_box(fg_struct->log10_Mturnover_MINI_filtered, box_dim, consts->hii_filter, R, + 0., 0.); + filter_box(fg_struct->log10_Mturnover_filtered, box_dim, consts->hii_filter, R, 0., 0.); - filter_box(fg_struct->log10_Mturnover_filtered, box_dim, consts->hii_filter, R, 0.); } } } diff --git a/src/py21cmfast/src/PerturbedField.c b/src/py21cmfast/src/PerturbedField.c index 8157c2b6c..6a1f613ca 100644 --- a/src/py21cmfast/src/PerturbedField.c +++ b/src/py21cmfast/src/PerturbedField.c @@ -152,7 +152,7 @@ void assign_to_lowres_grid(fftwf_complex *hires_grid, fftwf_complex *lowres_grid // Now filter the box filter_box(hires_grid, hi_dim, 0, - physconst.l_factor * simulation_options_global->BOX_LEN / (lo_dim[0] + 0.0), 0.); + physconst.l_factor * simulation_options_global->BOX_LEN / (lo_dim[0] + 0.0), 0., 0.); // FFT back to real space dft_c2r_cube(matter_options_global->USE_FFTW_WISDOM, hi_dim[0], hi_dim[2], @@ -222,7 +222,7 @@ void smooth_and_clip_density(fftwf_complex *lowres_grid, fftwf_complex *density_ simulation_options_global->DENSITY_SMOOTH_RADIUS * simulation_options_global->BOX_LEN / (float)simulation_options_global->HII_DIM, - 0.); + 0., 0.); } LOG_SUPER_DEBUG("delta_k after smoothing: "); @@ -354,7 +354,7 @@ void compute_perturbed_velocities(unsigned short axis, double redshift, filter_box(velocity_fft_grid, box_dim, 0, physconst.l_factor * simulation_options_global->BOX_LEN / (simulation_options_global->HII_DIM + 0.0), - 0.); + 0., 0.); } dft_c2r_cube(matter_options_global->USE_FFTW_WISDOM, box_dim[0], box_dim[2], diff --git a/src/py21cmfast/src/SpinTemperatureBox.c b/src/py21cmfast/src/SpinTemperatureBox.c index b350ce482..5cfeed32a 100644 --- a/src/py21cmfast/src/SpinTemperatureBox.c +++ b/src/py21cmfast/src/SpinTemperatureBox.c @@ -594,7 +594,7 @@ void fill_Rbox_table(float **result, fftwf_complex *unfiltered_box, double *R_ar // don't filter on cell size if (R > physconst.l_factor * (simulation_options_global->BOX_LEN / simulation_options_global->HII_DIM)) { - filter_box(box, box_dim, astro_options_global->HEAT_FILTER, R, 0.); + filter_box(box, box_dim, astro_options_global->HEAT_FILTER, R, 0., 0.); } LOG_ULTRA_DEBUG("db2 %d", R_ct); @@ -648,7 +648,7 @@ void fill_Rbox_table(float **result, fftwf_complex *unfiltered_box, double *R_ar // (better memory locality) // TODO: filter speed tests void one_annular_filter(float *input_box, float *output_box, double R_inner, double R_outer, - double *u_avg, double *f_avg) { + double R_star, int filter_type, double *u_avg, double *f_avg) { int i, j, k; unsigned long long int ct; double unfiltered_avg = 0; @@ -702,7 +702,7 @@ void one_annular_filter(float *input_box, float *output_box, double R_inner, dou // Don't filter on the cell scale if (R_inner > 0) { - filter_box(filtered_box, box_dim, 4, R_inner, R_outer); + filter_box(filtered_box, box_dim, filter_type, R_inner, R_outer, R_star); } // now fft back to real space @@ -748,11 +748,12 @@ void one_annular_filter(float *input_box, float *output_box, double R_inner, dou // fill a box[R_ct][box_ct] array for use in TS by filtering on different scales and storing results // Similar to fill_Rbox_table but called using different redshifts for each scale -int UpdateXraySourceBox(HaloBox *halobox, double R_inner, double R_outer, int R_ct, +int UpdateXraySourceBox(HaloBox *halobox, double R_inner, double R_outer, int R_ct, double R_star, XraySourceBox *source_box) { - int status; + int status, filter_type; Try { // the indexing needs these + filter_type = astro_options_global->LYA_MULTIPLE_SCATTERING ? 5 : 4; // only print once, since this is called for every R if (R_ct == 0) LOG_DEBUG("starting XraySourceBox"); @@ -761,19 +762,27 @@ int UpdateXraySourceBox(HaloBox *halobox, double R_inner, double R_outer, int R_ double xray_avg, fxray_avg; one_annular_filter(halobox->halo_sfr, &(source_box->filtered_sfr[R_ct * HII_TOT_NUM_PIXELS]), R_inner, R_outer, - &sfr_avg, &fsfr_avg); + R_star, filter_type, &sfr_avg, &fsfr_avg); one_annular_filter(halobox->halo_xray, &(source_box->filtered_xray[R_ct * HII_TOT_NUM_PIXELS]), R_inner, - R_outer, &xray_avg, &fxray_avg); - + R_outer, R_star, 4, &xray_avg, &fxray_avg); source_box->mean_sfr[R_ct] = fsfr_avg; if (astro_options_global->USE_MINI_HALOS) { one_annular_filter(halobox->halo_sfr_mini, &(source_box->filtered_sfr_mini[R_ct * HII_TOT_NUM_PIXELS]), R_inner, - R_outer, &sfr_avg_mini, &fsfr_avg_mini); - + R_outer, R_star, filter_type, &sfr_avg_mini, &fsfr_avg_mini); source_box->mean_sfr_mini[R_ct] = fsfr_avg_mini; source_box->mean_log10_Mcrit_LW[R_ct] = halobox->log10_Mcrit_MCG_ave; + // In case of multiple scattering and mini-halos, we need to filter the SFRD fields + // again for the the LW feedback, as these photons travel in straight lines + if (astro_options_global->LYA_MULTIPLE_SCATTERING) { + one_annular_filter(halobox->halo_sfr, + &(source_box->filtered_sfr_lw[R_ct * HII_TOT_NUM_PIXELS]), + R_inner, R_outer, R_star, 4, &sfr_avg, &fsfr_avg); + one_annular_filter(halobox->halo_sfr_mini, + &(source_box->filtered_sfr_mini_lw[R_ct * HII_TOT_NUM_PIXELS]), + R_inner, R_outer, R_star, 4, &sfr_avg_mini, &fsfr_avg_mini); + } } if (R_ct == astro_params_global->N_STEP_TS - 1) LOG_DEBUG("finished XraySourceBox"); @@ -891,7 +900,11 @@ void init_first_Ts(TsBox *box, float *dens, float z, float zp, double *x_e_ave, xe = xion_RECFAST(zp, 0); TK = T_RECFAST(zp, 0); - cT_ad = cT_approx(zp); + if (astro_options_global->USE_ADIABATIC_FLUCTUATIONS) { + cT_ad = cT_approx(zp); + } else { + cT_ad = 0.; + } growth_factor_zp = dicke(zp); inverse_growth_factor_z = 1 / dicke(z); @@ -1653,6 +1666,7 @@ void ts_main(float redshift, float prev_redshift, float perturbed_field_redshift int xidx; double ival, sfr_term, xray_sfr; double sfr_term_mini = 0; + double sfr_term_lw, sfr_term_mini_lw; #pragma omp for for (box_ct = 0; box_ct < HII_TOT_NUM_PIXELS; box_ct++) { // sum each R contribution together @@ -1668,6 +1682,14 @@ void ts_main(float redshift, float prev_redshift, float perturbed_field_redshift // Minihalos and s->yr conversion are already included here xray_sfr = source_box->filtered_xray[R_ct * HII_TOT_NUM_PIXELS + box_ct] * z_edge_factor * xray_R_factor * 1e38; + if (astro_options_global->USE_MINI_HALOS && + astro_options_global->LYA_MULTIPLE_SCATTERING) { + sfr_term_lw = + source_box->filtered_sfr_lw[R_ct * HII_TOT_NUM_PIXELS + box_ct] * + z_edge_factor; + } else { + sfr_term_lw = sfr_term; + } } else { // NOTE: for SOURCE_MODEL==0 F_STAR10 is still used for constant // stellar fraction @@ -1681,15 +1703,25 @@ void ts_main(float redshift, float prev_redshift, float perturbed_field_redshift sfr_term_mini = source_box->filtered_sfr_mini[R_ct * HII_TOT_NUM_PIXELS + box_ct] * z_edge_factor; + if (astro_options_global->LYA_MULTIPLE_SCATTERING) { + sfr_term_mini_lw = + source_box + ->filtered_sfr_mini_lw[R_ct * HII_TOT_NUM_PIXELS + box_ct] * + z_edge_factor; + } else { + sfr_term_mini_lw = sfr_term_mini; + } } else { sfr_term_mini = del_fcoll_Rct_MINI[box_ct] * z_edge_factor * avg_fix_term_MINI * astro_params_global->F_STAR7_MINI; xray_sfr += sfr_term_mini * astro_params_global->L_X_MINI * xray_R_factor * physconst.s_per_yr; + sfr_term_lw = sfr_term; + sfr_term_mini_lw = sfr_term_mini; } dstarlyLW_dt_box[box_ct] += - sfr_term * dstarlyLW_dt_prefactor[R_ct] + - sfr_term_mini * dstarlyLW_dt_prefactor_MINI[R_ct]; + sfr_term_lw * dstarlyLW_dt_prefactor[R_ct] + + sfr_term_mini_lw * dstarlyLW_dt_prefactor_MINI[R_ct]; } xidx = m_xHII_low_box[box_ct]; ival = inverse_val_box[box_ct]; diff --git a/src/py21cmfast/src/SpinTemperatureBox.h b/src/py21cmfast/src/SpinTemperatureBox.h index 2c607ff38..807232d77 100644 --- a/src/py21cmfast/src/SpinTemperatureBox.h +++ b/src/py21cmfast/src/SpinTemperatureBox.h @@ -8,7 +8,7 @@ int ComputeTsBox(float redshift, float prev_redshift, float perturbed_field_reds PerturbedField *perturbed_field, XraySourceBox *source_box, TsBox *previous_spin_temp, InitialConditions *ini_boxes, TsBox *this_spin_temp); -int UpdateXraySourceBox(HaloBox *halobox, double R_inner, double R_outer, int R_ct, +int UpdateXraySourceBox(HaloBox *halobox, double R_inner, double R_outer, int R_ct, double R_star, XraySourceBox *source_box); #endif diff --git a/src/py21cmfast/src/_functionprototypes_wrapper.h b/src/py21cmfast/src/_functionprototypes_wrapper.h index 59a0984d3..47e77fede 100644 --- a/src/py21cmfast/src/_functionprototypes_wrapper.h +++ b/src/py21cmfast/src/_functionprototypes_wrapper.h @@ -31,7 +31,7 @@ int ComputeBrightnessTemp(float redshift, TsBox *spin_temp, IonizedBox *ionized_ int ComputeHaloBox(double redshift, InitialConditions *ini_boxes, HaloCatalog *halos, TsBox *previous_spin_temp, IonizedBox *previous_ionize_box, HaloBox *grids); -int UpdateXraySourceBox(HaloBox *halobox, double R_inner, double R_outer, int R_ct, +int UpdateXraySourceBox(HaloBox *halobox, double R_inner, double R_outer, int R_ct, double R_star, XraySourceBox *source_box); /*--------------------------*/ @@ -118,7 +118,13 @@ int single_test_sample(unsigned long long int seed, int n_condition, float *cond int test_halo_props(double redshift, float *vcb_grid, float *J21_LW_grid, float *z_re_grid, float *Gamma12_ion_grid, int n_halos, float *halo_masses, float *halo_coords, float *star_rng, float *sfr_rnd, float *xray_rng, float *halo_props_out); -int test_filter(float *input_box, double R, double R_param, int filter_flag, double *result); +int test_filter(float *input_box, double R, double R_param, double R_star, int filter_flag, + double *result); + +// test functions for multiple scattering +double compute_mu_for_multiple_scattering(double x_em); +double compute_eta_for_multiple_scattering(double x_em); +double hyper_2F3(double kR, double alpha, double beta); /* Functions required to access cosmology & mass functions directly */ double dicke(double z); diff --git a/src/py21cmfast/src/_inputparams_wrapper.h b/src/py21cmfast/src/_inputparams_wrapper.h index 0f0410ae0..e13e0589e 100644 --- a/src/py21cmfast/src/_inputparams_wrapper.h +++ b/src/py21cmfast/src/_inputparams_wrapper.h @@ -143,6 +143,8 @@ typedef struct AstroOptions { bool FIX_VCB_AVG; bool USE_EXP_FILTER; bool CELL_RECOMB; + bool LYA_MULTIPLE_SCATTERING; + bool USE_ADIABATIC_FLUCTUATIONS; int PHOTON_CONS_TYPE; bool USE_UPPER_STELLAR_TURNOVER; bool HALO_SCALING_RELATIONS_MEDIAN; diff --git a/src/py21cmfast/src/_outputstructs_wrapper.h b/src/py21cmfast/src/_outputstructs_wrapper.h index 1a220fd78..52b040bb2 100644 --- a/src/py21cmfast/src/_outputstructs_wrapper.h +++ b/src/py21cmfast/src/_outputstructs_wrapper.h @@ -68,6 +68,8 @@ typedef struct XraySourceBox { float *filtered_sfr; float *filtered_xray; float *filtered_sfr_mini; + float *filtered_sfr_lw; + float *filtered_sfr_mini_lw; double *mean_log10_Mcrit_LW; double *mean_sfr; diff --git a/src/py21cmfast/src/filtering.c b/src/py21cmfast/src/filtering.c index 0f2029d25..a879d1407 100644 --- a/src/py21cmfast/src/filtering.c +++ b/src/py21cmfast/src/filtering.c @@ -1,6 +1,7 @@ #include #include +#include #include #include #include @@ -102,7 +103,7 @@ double exp_mfp_filter(double k, double R, double mfp, double exp_term) { return f; } -double spherical_shell_filter(double k, double R_outer, double R_inner) { +double spherical_shell_filter(double k, double R_inner, double R_outer) { double kR_inner = k * R_inner; double kR_outer = k * R_outer; @@ -115,18 +116,213 @@ double spherical_shell_filter(double k, double R_outer, double R_inner) { (sin(kR_outer) - cos(kR_outer) * kR_outer - sin(kR_inner) + cos(kR_inner) * kR_inner); } -void filter_box(fftwf_complex *box, int box_dim[3], int filter_type, float R, float R_param) { +struct multiple_scattering_params { + double alpha_outer; + double beta_outer; + double alpha_inner; + double beta_inner; +}; + +double compute_mu_for_multiple_scattering(double x_em) { + double mu; + double zeta_em = log10(x_em); + // Eq. (29) in arxiv: 2601.14360 + if (x_em > 30) { + mu = 1. - 1.0478 * pow(x_em, -0.7266); + } else if (x_em > 3.) { + mu = -0.104 * pow(zeta_em, 5) + 0.4867 * pow(zeta_em, 4) - 0.8217 * pow(zeta_em, 3) + + 0.4889 * zeta_em * zeta_em + 0.264 * zeta_em + 0.518; + } else if (x_em > 0.2) { + mu = -0.0285 * pow(zeta_em, 5) + 0.087 * pow(zeta_em, 4) - 0.1205 * pow(zeta_em, 3) - + 0.0456 * zeta_em * zeta_em + 0.3787 * zeta_em + 0.5285; + } else { + mu = 0.3982 * pow(x_em, 0.1592); + } + return mu; +} + +double compute_eta_for_multiple_scattering(double x_em) { + double eta; + double zeta_em = log10(x_em); + // Eq. (30) in arxiv: 2601.14360 + if (x_em > 20.) { + eta = 1. - 2.804 * pow(x_em, -1.242); + } else if (x_em > 3.) { + eta = 2.17 * pow(zeta_em, 5) - 8.832 * pow(zeta_em, 4) + 13.579 * pow(zeta_em, 3) - + 10.04 * zeta_em * zeta_em + 4.166 * zeta_em - 0.17; + } else if (x_em > 0.2) { + eta = 0.352 * pow(zeta_em, 5) - 0.0516 * pow(zeta_em, 4) - 0.293 * pow(zeta_em, 3) + + 0.342 * zeta_em * zeta_em + 0.582 * zeta_em + 0.266; + } else { + eta = 0.4453 * pow(x_em, 1.296); + } + return eta; +} + +void initialize_alphas_and_betas_for_multiple_scattering( + double R_inner, double R_outer, double R_star, struct multiple_scattering_params *consts) { + if (R_star == 0.) { + // R_star == 0 could happen after reionization since R_star is proportional to x_HI + consts->alpha_inner = 1.; + consts->alpha_outer = 1.; + consts->beta_inner = 1.; + consts->beta_outer = 0.; + } else { + // Eq. (25) in arxiv: 2601.14360 + double x_em_inner = R_inner / R_star; + double x_em_outer = R_outer / R_star; + + double mu_inner = compute_mu_for_multiple_scattering(x_em_inner); + double eta_inner = compute_eta_for_multiple_scattering(x_em_inner); + double mu_outer = compute_mu_for_multiple_scattering(x_em_outer); + double eta_outer = compute_eta_for_multiple_scattering(x_em_outer); + + // Eq. (28) in arxiv: 2601.14360 (mu = alpha/(alpha+beta), eta = alpha/(alpha+beta^2)) + consts->alpha_inner = (1. / eta_inner - 1.) / pow(1. / mu_inner - 1., 2); + consts->beta_inner = (1. / eta_inner - 1.) / (1. / mu_inner - 1.); + consts->alpha_outer = (1. / eta_outer - 1.) / pow(1. / mu_outer - 1., 2); + consts->beta_outer = (1. / eta_outer - 1.) / (1. / mu_outer - 1.); + } +} + +double asymptotic_2F3(double kR, double alpha, double beta) { + // Approximation to the hypergeometric function 2F3(a1,a2;b1,b2,b3;-z) for z >> 1, + // as given in + // https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F3/06/02/03/01/0003/. In + // our case, z=-kR^2/4 and the parameters of the hypergeometric function are given below + // (see also Eq. D8 in arxiv: 2601.14360) + double a1 = (2. + alpha) / 2.; + double a2 = (3. + alpha) / 2.; + double b1 = 5. / 2.; + double b2 = (2. + alpha + beta) / 2.; + double b3 = (3. + alpha + beta) / 2.; + + double gamma_a1 = tgamma(a1); + double gamma_a2 = tgamma(a2); + double gamma_b1 = 3. / 4.; // Actually Gamma(5/2) is 3/4*sqrt(pi), but we absorbed the sqrt(pi) + // in the other terms below + double gamma_b2 = tgamma(b2); + double gamma_b3 = tgamma(b3); + + double gamma_b2_over_a1, gamma_b3_over_a2, decay_term1, decay_term2, F_asymp; + // If alpha >> beta, the gamma function becomes very large and overflow can happen if they are + // naively computed. In this limit, we use the following asymptotic approximation, which is + // based on Stirling's formula + if (a1 < 20.) { + gamma_b2_over_a1 = gamma_b2 / gamma_a1; + gamma_b3_over_a2 = gamma_b3 / gamma_a2; + } else { + double y = beta / 2; // b2-a1 = b3-a2 + gamma_b2_over_a1 = pow(a1, y) * exp((a1 + y - 0.5) * (y / a1 - y * y / (2. * a1 * a1) + + y * y * y / (3. * a1 * a1 * a1)) - + y); + gamma_b3_over_a2 = pow(a2, y) * exp((a2 + y - 0.5) * (y / a2 - y * y / (2. * a2 * a1) + + y * y * y / (3. * a2 * a2 * a2)) - + y); + } + // If alpha is very big the following decaying terms can be completely neglected + if (alpha < 10.) { + // There are some Gamma function whose argument could be close to a pole (non-positive + // integers), but they are at the denominator, so they behave nicely and we actually need + // the reciprocal gamma function + double gamma_b1_minus_a1_inv = gsl_sf_gammainv(b1 - a1); // 1/Gamma((3-alpha)/2) + double gamma_b1_minus_a2_inv = gsl_sf_gammainv(b1 - a2); // 1/Gamma((2-alpha)/2) + double gamma_b2_minus_a2_inv = gsl_sf_gammainv(b2 - a2); // 1/Gamma((beta-1)/2) + decay_term1 = M_PI * gamma_a1 * gamma_b1_minus_a1_inv / tgamma(b2 - a1) / tgamma(b3 - a1) / + pow(kR / 2., alpha + 2.); // gamma(a2-a1) = sqrt(pi) + decay_term2 = -2. * M_PI * gamma_a2 * gamma_b1_minus_a2_inv * gamma_b2_minus_a2_inv / + tgamma(b3 - a2) / pow(kR / 2., alpha + 3.); // gamma(a1-a2) = -2*sqrt(pi) + } else { + decay_term1 = 0.; + decay_term2 = 0.; + } + + F_asymp = (cos(kR - M_PI * (2. + beta) / 2.) - + (1. + (alpha - 1.) * beta) / kR * sin(kR - M_PI * (2. + beta) / 2.)) / + pow(kR / 2, beta + 2); + F_asymp += decay_term1 + decay_term2; + F_asymp *= gamma_b1 * gamma_b2_over_a1 * gamma_b3_over_a2; + + return F_asymp; +} + +// Implementation of 2F3((alpha+2)/2, (alpha+3)/2 ; (alpha+beta+2)/2, (alpha+beta+3)/2 ; -kR^2 /4)) +// This is Eq. (32) in arxiv: 2601.14360 +double hyper_2F3(double kR, double alpha, double beta) { + if (beta == 0.) { + // beta=0 during reionization. + // In this case we return the straight-line window function (no multiple scattering if there + // are no HI atoms!) + return 3.0 / (pow(kR, 3)) * (sin(kR) - cos(kR) * kR); + } + // For a small argument, we compute the hypergeometric function through power-law expansion + // This is Eq. (D7) in arxiv: 2601.14360 + if (kR < 30.) { + int n; + int max_terms = 1000; + double sum = 0.; + double term = 1.; + for (n = 1; n < max_terms; n++) { + sum += term; + term *= -1. / (1. + beta / (alpha + 2. * n)) / (1. + beta / (alpha + 1 + 2. * n)) * kR * + kR / (2. * n) / (2. * n + 3.); + if (fabs(term) < fabs(sum) * 1e-4) { + break; + } + } + return sum; + } + // For large arguments, the above sum becomes numerically unstable, and we use instead + // asymptotic approximation + else { + double F_ms, F_sl; + F_ms = asymptotic_2F3(kR, alpha, beta); + F_sl = 3.0 / (pow(kR, 3)) * (sin(kR) - cos(kR) * kR); + /* At large arguments, the hypergeometric function (multiple scattering window function) + should be below the straight-line window function. However, for large alpha values the + asymptotic approximation is not adequate at 30 < kR < 100 and the asymptotic formula + diverges at this range. We could of course increase the threshold for "big" argument, but + then the above power-law expansion diverges. I therefore use this rule of thumb, which is + necessary only for large alpha values (where the hypergeometric function approaches the + straight-line window function) and intermediate kR values. */ + if (fabs(F_ms) < fabs(F_sl)) { + return F_ms; + } else { + return F_sl; + } + } +} + +double multiple_scattering_filter(double k, double R_inner, double R_outer, + struct multiple_scattering_params *consts) { + double kR_inner = k * R_inner; + double kR_outer = k * R_outer; + double W; + // Eq. (11) in arxiv: 2601.14360 + W = pow(R_outer, 3.) * hyper_2F3(kR_outer, consts->alpha_outer, consts->beta_outer) - + pow(R_inner, 3.) * hyper_2F3(kR_inner, consts->alpha_inner, consts->beta_inner); + W /= pow(R_outer, 3.) - pow(R_inner, 3.); + return W; +} + +void filter_box(fftwf_complex *box, int box_dim[3], int filter_type, float R, float R_param, + float R_star) { double delta_k[3] = { 2.0 * M_PI / simulation_options_global->BOX_LEN, 2.0 * M_PI / simulation_options_global->BOX_LEN, 2.0 * M_PI / (simulation_options_global->BOX_LEN * simulation_options_global->NON_CUBIC_FACTOR)}; + struct multiple_scattering_params consts_for_ms; + // setup constants if needed double R_const; if (filter_type == 3) { R_const = exp(-R / R_param); } + if (filter_type == 5) { + initialize_alphas_and_betas_for_multiple_scattering(R, R_param, R_star, &consts_for_ms); + } // loop through k-box #pragma omp parallel num_threads(simulation_options_global->N_THREADS) @@ -181,6 +377,9 @@ void filter_box(fftwf_complex *box, int box_dim[3], int filter_type, float R, fl box[grid_index] *= exp_mfp_filter(sqrt(k_mag_sq), R, R_param, R_const); } else if (filter_type == 4) { // spherical shell, R_param == inner radius box[grid_index] *= spherical_shell_filter(sqrt(k_mag_sq), R, R_param); + } else if (filter_type == 5) { // multiple scattering window function + box[grid_index] *= + multiple_scattering_filter(sqrt(k_mag_sq), R, R_param, &consts_for_ms); } else { if ((n_x == 0) && (n_y == 0) && (n_z == 0)) LOG_WARNING("Filter type %i is undefined. Box is unfiltered.", @@ -195,7 +394,8 @@ void filter_box(fftwf_complex *box, int box_dim[3], int filter_type, float R, fl } // Test function to filter a box without computing a whole output box -int test_filter(float *input_box, double R, double R_param, int filter_flag, double *result) { +int test_filter(float *input_box, double R, double R_param, double R_star, int filter_flag, + double *result) { int i, j, k; unsigned long long int ii, jj; int box_dim[3] = { @@ -225,7 +425,7 @@ int test_filter(float *input_box, double R, double R_param, int filter_flag, dou memcpy(box_filtered, box_unfiltered, sizeof(fftwf_complex) * HII_KSPACE_NUM_PIXELS); - filter_box(box_filtered, box_dim, filter_flag, R, R_param); + filter_box(box_filtered, box_dim, filter_flag, R, R_param, R_star); dft_c2r_cube(matter_options_global->USE_FFTW_WISDOM, simulation_options_global->HII_DIM, HII_D_PARA, simulation_options_global->N_THREADS, box_filtered); diff --git a/src/py21cmfast/src/filtering.h b/src/py21cmfast/src/filtering.h index 48ab03ed3..14782ba57 100644 --- a/src/py21cmfast/src/filtering.h +++ b/src/py21cmfast/src/filtering.h @@ -5,8 +5,10 @@ #include "InputParameters.h" -void filter_box(fftwf_complex *box, int box_dim[3], int filter_type, float R, float R_param); -int test_filter(float *input_box, double R, double R_param, int filter_flag, double *result); +void filter_box(fftwf_complex *box, int box_dim[3], int filter_type, float R, float R_param, + float R_star); +int test_filter(float *input_box, double R, double R_param, double R_star, int filter_flag, + double *result); double filter_function(double k, int filter_type); double dwdm_filter(double k, double R, int filter_type); diff --git a/src/py21cmfast/utils.py b/src/py21cmfast/utils.py index 5feddd4c7..c356f10e2 100644 --- a/src/py21cmfast/utils.py +++ b/src/py21cmfast/utils.py @@ -42,7 +42,7 @@ def recursive_difference( return out -def show_references(inputs: InputParameters, print_to_stdout=True): +def show_references(inputs: InputParameters, lightcone=True, print_to_stdout=True): """Print out all the relevant references for a simulation based on input parameters.""" ref_string = ( "The main reference for 21cmFAST v3+:\n" @@ -58,6 +58,29 @@ def show_references(inputs: InputParameters, print_to_stdout=True): "https://ui.adsabs.harvard.edu/link_gateway/2011MNRAS.411..955M/doi:10.1111/j.1365-2966.2010.17731.x\n\n" ) + if inputs.astro_options.INHOMO_RECO: + ref_string += ( + "Inhomogeneous recombination model introduced in:\n" + "=============================================\n" + "Sobacchi, E., Mesinger, A.,\n" + "“Inhomogeneous recombinations during cosmic reionization”,\n" + "Monthly Notices of the Royal Astronomical Society,\n" + "vol. 440, no. 2, pp. 1662-1673, 2014.\n" + "https://doi.org/10.1093/mnras/stu377.\n\n" + ) + + if lightcone: + ref_string += ( + "Lightcone and redshift space distortions first introduced in:\n" + "=============================================\n" + "Greig, B., Mesinger, A.,\n" + "“21CMMC with a 3D light-cone: the impact of the co-evolution approximation on \n" + "the astrophysics of reionisation and cosmic dawn.”,\n" + "Monthly Notices of the Royal Astronomical Society,\n" + "vol. 477, no. 3, 3217-3229, 2018.\n" + "https://doi.org/10.1093/mnras/sty796.\n\n" + ) + if inputs.matter_options.SOURCE_MODEL in [ "E-INTEGRAL", "L-INTEGRAL", @@ -81,12 +104,27 @@ def show_references(inputs: InputParameters, print_to_stdout=True): "“A tale of two sites - I. Inferring the properties of minihalo-hosted galaxies from\n" "current observations”, Monthly Notices of the Royal Astronomical Society, vol. 495,\n" "no. 1, pp. 123–140, 2020. https://doi.org/10.1093/mnras/staa1131.\n\n" + ) + + if inputs.matter_options.USE_RELATIVE_VELOCITIES: + ref_string += ( "With improvements and DM-baryon relative velocity model in:\n" "===========================================================\n" "Muñoz, J.B., Qin, Y., Mesinger, A., Murray, S., Greig, B., and Mason, C.,\n" '"The Impact of the First Galaxies on Cosmic Dawn and Reionization",\n' "Monthly Notices of the Royal Astronomical Society, vol. 511, no. 3,\n" - "pp 3657-3681, 2022 https://doi.org/10.1093/mnras/stac185\n" + "pp 3657-3681, 2022 https://doi.org/10.1093/mnras/stac185\n\n" + ) + + if inputs.astro_options.PHOTON_CONS_TYPE != "no-photoncons": + ref_string += ( + "Photon conservation correction introduced in:\n" + "=============================================\n" + "Park, J., Greig, B., Mesinger, A.,\n" + "“Calibrating excursion set reionization models to approximately conserve ionizing photons”,\n" + "Monthly Notices of the Royal Astronomical Society,\n" + "vol. 517, no. 1, pp. 192-200, 2022.\n" + "https://doi.org/10.1093/mnras/stac2756.\n\n" ) if inputs.matter_options.lagrangian_source_grid: @@ -96,7 +134,17 @@ def show_references(inputs: InputParameters, print_to_stdout=True): "Davies, J. E., Mesinger, A., Murray, S. G.,\n" '"Efficient simulation of discrete galaxy populations and associated radiation\n' 'fields during the first billion years",\n' - "eprint arXiv:2504.17254, 2025. https://doi.org/10.48550/arXiv.2504.17254\n\n" + "Astronomy and Astrophysics, vol. 701, A. 236, 2025.\n" + "https://doi.org/10.1051/0004-6361/202554951.\n\n" + ) + + if inputs.astro_options.LYA_MULTIPLE_SCATTERING: + ref_string += ( + "Lyman alpha multiple scattering effect introduced in:\n" + "========================================================\n" + "Flitter, J., Munoz, J. B., Mesinger, A.,\n" + '"Semi-analytical approach to Lyman-alpha multiple-scattering in 21-cm signal simulations",\n' + "eprint arXiv:2601.14360, 2026. https://doi.org/10.48550/arXiv.2601.14360.\n\n" ) if print_to_stdout: diff --git a/src/py21cmfast/wrapper/cfuncs.py b/src/py21cmfast/wrapper/cfuncs.py index 29eaf0090..7a452609a 100644 --- a/src/py21cmfast/wrapper/cfuncs.py +++ b/src/py21cmfast/wrapper/cfuncs.py @@ -27,12 +27,12 @@ def broadcast_input_struct(inputs: InputParameters): """Broadcast the parameters to the C library.""" lib.Broadcast_struct_global_all( - inputs.simulation_options.cstruct, - inputs.matter_options.cstruct, - inputs.cosmo_params.cstruct, - inputs.astro_params.cstruct, - inputs.astro_options.cstruct, - inputs.cosmo_tables.cstruct, + inputs.simulation_options._cstruct, + inputs.matter_options._cstruct, + inputs.cosmo_params._cstruct, + inputs.astro_params._cstruct, + inputs.astro_options._cstruct, + inputs.cosmo_tables._cstruct, ) diff --git a/src/py21cmfast/wrapper/inputs.py b/src/py21cmfast/wrapper/inputs.py index ad1ee0350..053cb16b9 100644 --- a/src/py21cmfast/wrapper/inputs.py +++ b/src/py21cmfast/wrapper/inputs.py @@ -158,17 +158,17 @@ def new(cls, x: dict | Self | None = None, **kwargs): Examples -------- - >>> up = SimulationOptions({'HII_DIM': 250}) - >>> up.HII_DIM + >>> params = SimulationOptions.new({'HII_DIM': 250}) + >>> params.HII_DIM 250 - >>> up = SimulationOptions(up) - >>> up.HII_DIM + >>> params = SimulationOptions.new(params) + >>> params.HII_DIM 250 - >>> up = SimulationOptions() - >>> up.HII_DIM + >>> params = SimulationOptions.new() + >>> params.HII_DIM 200 - >>> up = SimulationOptions(HII_DIM=256) - >>> up.HII_DIM + >>> params = SimulationOptions.new(HII_DIM=256) + >>> params.HII_DIM 256 """ if isinstance(x, dict): @@ -187,15 +187,15 @@ def __init_subclass__(cls) -> None: InputStruct._subclasses[cls.__name__] = cls @cached_property - def struct(self) -> StructWrapper: + def _struct(self) -> StructWrapper: """The python-wrapped struct associated with this input object.""" return StructWrapper(self.__class__.__name__) @cached_property - def cstruct(self) -> StructWrapper: + def _cstruct(self) -> StructWrapper: """The object pointing to the memory accessed by C-code for this struct.""" cdict = self.cdict - for k in self.struct.fieldnames: + for k in self._struct.fieldnames: val = cdict[k] # TODO: is this really required here? (I don't think the wrapper can satisfy this condition) @@ -203,9 +203,9 @@ def cstruct(self) -> StructWrapper: # If it is a string, need to convert it to C string ourselves. val = self.ffi.new("char[]", val.encode()) - setattr(self.struct.cstruct, k, val) + setattr(self._struct.cstruct, k, val) - return self.struct.cstruct + return self._struct.cstruct def clone(self, **kwargs): """Make a fresh copy of the instance with arbitrary parameters updated.""" @@ -245,7 +245,7 @@ def cdict(self) -> dict: } out = {} - for k in self.struct.fieldnames: + for k in self._struct.fieldnames: val = getattr(self, k) att = attrs.fields_dict(self.__class__).get(k, None) # we assume properties (as opposed to attributes) are already converted @@ -261,8 +261,8 @@ def __str__(self) -> str: return f"""{self.__class__.__name__}:{params} """ @classmethod - def from_subdict(cls, dct, safe=True): - """Construct an instance of a parameter structure from a dictionary.""" + def from_subdict(cls, dct: dict[str, Any], safe: bool = True) -> Self: + """Construct an instance of this class from a dictionary.""" fieldnames = [ field.name for field in attrs.fields(cls) @@ -317,7 +317,7 @@ class Table1D: ) @cached_property - def cstruct(self): + def _cstruct(self): """Cached pointer to the memory of the object in C.""" ctab = ffi.new("Table1D *") ctab.size = self.size @@ -365,21 +365,21 @@ def new(cls, x: dict | Self | None = None, **kwargs): ) @cached_property - def struct(self): + def _struct(self): """The python-wrapped struct associated with this input object.""" return StructWrapper("CosmoTables") @cached_property - def cstruct(self) -> StructWrapper: + def _cstruct(self) -> StructWrapper: """The object pointing to the memory accessed by C-code for this struct.""" - for k in self.struct.fieldnames: + for k in self._struct.fieldnames: val = getattr(self, k) if isinstance(val, Table1D): - setattr(self.struct.cstruct, k, val.cstruct) - elif val is not None: - setattr(self.struct.cstruct, k, val) + setattr(self._struct.cstruct, k, val._cstruct) + elif isinstance(val, (float | bool)) or val is not None: + setattr(self._struct.cstruct, k, val) - return self.struct.cstruct + return self._struct.cstruct def clone(self, **kwargs): """Make a fresh copy of the instance with arbitrary parameters updated.""" @@ -388,30 +388,59 @@ def clone(self, **kwargs): @attrs.define(frozen=True, kw_only=True) class CosmoParams(InputStruct): - """ - Cosmological parameters (with defaults) which translates to a C struct. - - To see default values for each parameter, use ``CosmoParams._defaults_``. - All parameters passed in the constructor are also saved as instance attributes which should - be considered read-only. This is true of all input-parameter classes. + r""" + Cosmological parameters. - Default parameters are based on Plank18, https://arxiv.org/pdf/1807.06209.pdf, - Table 2, last column. [TT,TE,EE+lowE+lensing+BAO] + Default parameters are based on `Plank18 `_, + Table 2, last column [TT,TE,EE+lowE+lensing+BAO]. - Parameters + Attributes ---------- - SIGMA_8 : float, optional - RMS mass variance (power spectrum normalisation). - hlittle : float, optional - The hubble parameter, H_0/100. - OMm : float, optional - Omega matter. - OMb : float, optional - Omega baryon, the baryon component. - POWER_INDEX : float, optional - Spectral index of the power spectrum. - A_s: float, optional - Amplitude of primordial curvature power spectrum, at k_pivot = 0.05 Mpc^-1. + hlittle : float + The hubble parameter, :math:`H_0/100`. + OMm : float + The fractional matter density, :math:`\Omega_m`. This is the sum of the baryon + and CDM densities. + OMb : float + The fraction baryon density, :math:`\Omega_b`. + POWER_INDEX : float + Spectral index of the power spectrum, :math:`n_s` + SIGMA_8 : float + RMS mass variance of the matter field in spheres of 8 Mpc/h (power spectrum + normalisation). If not given explicitly but A_s is given, it is auto-calculated + via A_s and the other cosmological parameters. Only one of SIGMA_8 and A_s can + be given, otherwise an error is raised. If neither is given, it is set to the + default value of 0.8102. + A_s : float + Amplitude of primordial curvature power spectrum, at + :math:`k_{pivot} = 0.05 {\rm Mpc}^{-1}`. + If not given explicitly and SIGMA_8 is given, it is auto-calculated via SIGMA_8 + and the other cosmological parameters, using CLASS. Only one of SIGMA_8 and A_s + can be given, otherwise an error is raised. If neither is given, it is set to + the default value of 2.105e-9. + base_cosmo : :class:`astropy.cosmology.FLRW`, optional + The base astropy cosmology object to use for this cosmology. + You can set this to something other than LambdaCDM if you want to use a + different cosmology. Note that the cosmological parameters in ``base_cosmo`` + are not used at all---only the cosmography. You almost certainly should not + set this to anything other than the default, unless you know what you are doing. + OMn : float + Omega neutrino, the neutrino component. Note that this is not used in any of the + calculations in the code, but is just provided for completeness. + OMk : float + Omega curvature, the curvature component. Note that this is not used in any of + the calculations in the code, but is just provided for completeness. + OMr : float + Omega radiation, the radiation component. Note that this is not used in any of + the calculations in the code, but is just provided for completeness. + OMtot : float + Omega total, the total density. Note that this is not used in any of the + calculations and is assumed to be unity (i.e. a flat universe), but is just + provided for completeness. + Y_He : float + Helium mass fraction. + wl : float + Dark energy equation of state parameter. """ _DEFAULT_SIGMA_8: ClassVar[float] = 0.8102 @@ -420,23 +449,28 @@ class CosmoParams(InputStruct): _base_cosmo: Annotated[FLRW, Parameter(show=False, parse=False)] = field( default=Planck18, validator=validators.instance_of(FLRW), eq=False, repr=False ) + _SIGMA_8: float = field( default=None, converter=attrs.converters.optional(float), validator=validators.optional(validators.gt(0)), ) + hlittle: float = field( default=Planck18.h, converter=float, validator=validators.gt(0) ) OMm: float = field( default=Planck18.Om0, converter=float, validator=validators.gt(0) ) + OMb: float = field( default=Planck18.Ob0, converter=float, validator=validators.gt(0) ) + POWER_INDEX: float = field( default=0.9665, converter=float, validator=validators.gt(0) ) + _A_s: float = field( default=None, converter=attrs.converters.optional(float), @@ -444,12 +478,17 @@ class CosmoParams(InputStruct): ) OMn: float = field(default=0.0, converter=float, validator=validators.ge(0)) + OMk: float = field(default=0.0, converter=float, validator=validators.ge(0)) + OMr: float = field(default=8.6e-5, converter=float, validator=validators.ge(0)) + OMtot: float = field( default=1.0, converter=float, validator=validators.ge(0) ) # TODO: force this to be the sum of the others + Y_He: float = field(default=0.24, converter=float, validator=validators.ge(0)) + wl: float = field(default=-1.0, converter=float) # TODO: Combined validation via Astropy? @@ -555,97 +594,113 @@ def asdict(self) -> dict: @attrs.define(frozen=True, kw_only=True) class MatterOptions(InputStruct): - """ + r""" Structure containing options which affect the matter field (ICs, perturbedfield, halos). - Parameters + Attributes ---------- - HMF: str, optional + HMF Determines which halo mass function to be used for the normalisation of the collapsed fraction (default Sheth-Tormen). Should be one of the following codes: - PS (Press-Schechter) - ST (Sheth-Tormen) - Watson (Watson FOF) - Watson-z (Watson FOF-z) - Delos (Delos+23) - USE_RELATIVE_VELOCITIES: int, optional + + * PS (Press-Schechter) + * ST (Sheth-Tormen) + * Watson (Watson FOF) + * Watson-z (Watson FOF-z) + * Delos (Delos+23) + + USE_RELATIVE_VELOCITIES Flag to decide whether to use relative velocities. If True, POWER_SPECTRUM is automatically set to 5. Default False. - POWER_SPECTRUM: str, optional + POWER_SPECTRUM Determines which power spectrum to use, default EH (unless `USE_RELATIVE_VELOCITIES` is True). Use the following codes: - EH : Eisenstein & Hu 1999 - BBKS: Bardeen et al. 1986 - EFSTATHIOU: Efstathiou et al. 1992 - PEEBLES: Peebles 1980 - WHITE: White 1985 - CLASS: Uses fits from the CLASS code - PERTURB_ON_HIGH_RES : bool, optional + + * EH : Eisenstein & Hu 1999 + * BBKS: Bardeen et al. 1986 + * EFSTATHIOU: Efstathiou et al. 1992 + * PEEBLES: Peebles 1980 + * WHITE: White 1985 + * CLASS: Uses fits from the CLASS code + PERTURB_ON_HIGH_RES Whether to perform the Zel'Dovich or 2LPT perturbation on the low or high resolution grid. - USE_FFTW_WISDOM : bool, optional + USE_FFTW_WISDOM Whether or not to use stored FFTW_WISDOMs for improving performance of FFTs - USE_INTERPOLATION_TABLES: str, optional + USE_INTERPOLATION_TABLES Defines the interpolation tables used in the code. Default is 'hmf-interpolation'. There are three levels available: - 'no-interpolation': No interpolation tables are used. - 'sigma-interpolation': Interpolation tables are used for sigma(M) only. - 'hmf-interpolation': Interpolation tables are used for sigma(M) the halo mass function. - PERTURB_ALGORITHM: str, optional - Whether to use second-order Lagrangian perturbation theory (2LPT), Zel'dovich (ZELDOVICH), - or linear evolution (LINEAR). + + * 'no-interpolation': No interpolation tables are used. + * 'sigma-interpolation': Interpolation tables are used for sigma(M) only. + * 'hmf-interpolation': Interpolation tables are used for sigma(M) the halo mass + function. + PERTURB_ALGORITHM + Whether to use second-order Lagrangian perturbation theory (2LPT), Zel'dovich + (ZELDOVICH), or linear evolution (LINEAR). Set this to 2LPT if the density field or the halo positions are extrapolated to low redshifts. The current implementation is very naive and adds a factor ~6 to the memory requirements. Reference: Scoccimarro R., 1998, MNRAS, 299, 1097-1118 Appendix D. - MINIMIZE_MEMORY: bool, optional + MINIMIZE_MEMORY If set, the code will run in a mode that minimizes memory usage, at the expense of some CPU/disk-IO. Good for large boxes / small computers. - SAMPLE_METHOD: str, optional + SAMPLE_METHOD The sampling method to use in the halo sampler when calculating progenitor populations: - MASS-LIMITED : Mass-limited CMF sampling, where samples are drawn until the expected mass is reached - NUMBER-LIMITED : Number-limited CMF sampling, where we select a number of halos from the Poisson distribution - and then sample the CMF that many times - PARTITION : Sheth et al 1999 Partition sampling, where the EPS collapsed fraction is sampled (gaussian tail) - and then the condition is updated using the conservation of mass. - BINARY-SPLIT : Parkinsson et al 2008 Binary split model as in DarkForest (Qiu et al 2021) where the EPS merger rate - is sampled on small internal timesteps such that only binary splits can occur. - NOTE: The initial sampling from the density grid will ALWAYS use number-limited sampling (method 1) - FILTER : string, optional + + * MASS-LIMITED : Mass-limited CMF sampling, where samples are drawn until the + expected mass is reached + * NUMBER-LIMITED : Number-limited CMF sampling, where we select a number of + halos from the Poisson distribution and then sample the CMF that many times + * PARTITION : Sheth et al 1999 Partition sampling, where the EPS collapsed + fraction is sampled (gaussian tail) and then the condition is updated using + the conservation of mass. + * BINARY-SPLIT : Parkinson et al 2008 Binary split model as in DarkForest + (Qiu et al 2021) where the EPS merger rate is sampled on small internal + timesteps such that only binary splits can occur. + + .. note:: The initial sampling from the density grid will ALWAYS use NUMBER-LIMITED + sampling (method 1) + FILTER Filter to use for sigma (matter field variance) and radius to mass conversions. available options are: `spherical-tophat` and `gaussian` - HALO_FILTER : string, optional + HALO_FILTER Filter to use for the DexM halo finder. available options are: `spherical-tophat`, `sharp-k` and `gaussian` - SMOOTH_EVOLVED_DENSITY_FIELD: bool, optional + SMOOTH_EVOLVED_DENSITY_FIELD Smooth the evolved density field after perturbation. - DEXM_OPTMIZE: bool, optional - Use a faster version of the DexM halo finder which excludes halos from forming within a certain distance of larger halos. - KEEP_3D_VELOCITIES: bool, optional + DEXM_OPTMIZE + Use a faster version of the DexM halo finder which excludes halos from forming + within a certain distance of larger halos. + KEEP_3D_VELOCITIES Whether to keep the 3D velocities in the ICs. If False, only the z velocity is kept. - SOURCE_MODEL: str, optional + SOURCE_MODEL The source model to use in the simulation. Options are: - E-INTEGRAL : The traditional excursion-set formalism, where source properties are - defined on the Eulerian grid after 2LPT in regions of filter scale R (see the X_FILTER options for filter shapes). - This integrates over the CHMF using the smoothed density grids, then multiplies the result. - by (1 + delta) to get the source properties in each cell. - CONST-ION-EFF: Similar to E-INTEGRAL, but ionizing efficiency is constant and does not depend on the halo mass - (see Mesinger+ 2010). - L-INTEGRAL : Analagous to the 'ESF-L' model described in Trac+22, where source properties - are defined on the Lagrangian (IC) grid by integrating the CHMF prior to the IGM physics - and then mapping properties to the Eulerian grid using 2LPT. - DEXM-ESF : The DexM excursion-set formalism, where discrete halo catalogues are generated - on the Lagrangian (IC) grid using an excursion-set halo finder. Source properties - are defined on the Lagrangian grid and then mapped to the Eulerian grid using 2LPT. - This model utilised the 'L-INTEGRAL' method for halos below the DexM mass resolution, - which is the mass of the high-resolution (DIM^3) cells. - CHMF-SAMPLER : The CHMF sampler, where discrete halo catalogues are generated by sampling - the CHMF on the IC grid, between the low-resolution (HII_DIM^3) cell mass and a minimum - mass defined by the user (SAMPLER_MIN_MASS). This model uses the 'L-INTEGRAL' method for - halos below the SAMPLER_MIN_MASS, and the 'DEXM-ESF' method for halos above the HII_DIM - cell mass. + + * ``E-INTEGRAL``: The traditional excursion-set formalism, where source properties + are defined on the Eulerian grid after 2LPT in regions of filter scale R + (see the X_FILTER options for filter shapes). This integrates over the CHMF + using the smoothed density grids, then multiplies the result by (1 + delta) + to get the source properties in each cell. + * ``CONST-ION-EFF``: Similar to E-INTEGRAL, but ionizing efficiency is constant and + does not depend on the halo mass (see Mesinger+ 2010). + * ``L-INTEGRAL``: Analagous to the 'ESF-L' model described in Trac+22, where source + properties are defined on the Lagrangian (IC) grid by integrating the CHMF + prior to the IGM physics and then mapping properties to the Eulerian grid + using 2LPT. + * ``DEXM-ESF``: The DexM excursion-set formalism, where discrete halo catalogues + are generated on the Lagrangian (IC) grid using an excursion-set halo + finder. Source properties are defined on the Lagrangian grid and then mapped + to the Eulerian grid using 2LPT. This model utilised the 'L-INTEGRAL' method + for halos below the DexM mass resolution, which is the mass of the + high-resolution (DIM^3) cells. + * ``CHMF-SAMPLER``: The CHMF sampler, where discrete halo catalogues are generated + by sampling the CHMF on the IC grid, between the low-resolution (HII_DIM^3) + cell mass and a minimum mass defined by the user (SAMPLER_MIN_MASS). This + model uses the 'L-INTEGRAL' method for halos below the SAMPLER_MIN_MASS, and + the 'DEXM-ESF' method for halos above the HII_DIM cell mass. """ HMF: Literal["PS", "ST", "WATSON", "WATSON-Z", "DELOS"] = choice_field(default="ST") @@ -736,7 +791,7 @@ class SimulationOptions(InputStruct): """ Structure containing broad simulation options. - Parameters + Attributes ---------- HII_DIM : int, optional Number of cells for the low-res box (after smoothing the high-resolution matter @@ -764,18 +819,20 @@ class SimulationOptions(InputStruct): defines this setting. Specifying both will result in an error. By default, the BOX_LEN will be calculated as 1.5 * HII_DIM. NON_CUBIC_FACTOR : float, optional - Factor which allows the creation of non-cubic boxes. It will shorten/lengthen the line - of sight dimension of all boxes. NON_CUBIC_FACTOR * DIM/HII_DIM must result in an integer. + Factor which allows the creation of non-cubic boxes. It will shorten/lengthen + the line of sight dimension of all boxes. ``NON_CUBIC_FACTOR * DIM/HII_DIM`` + must result in an integer. N_THREADS : int, optional Sets the number of processors (threads) to be used for performing 21cmFAST. Default 1. SAMPLER_MIN_MASS: float, optional - The minimum mass to sample in the halo sampler when SOURCE_MODEL is "CHMF-SAMPLER", - decreasing this can drastically increase both compute time and memory usage. + The minimum mass to sample in the halo sampler when ``SOURCE_MODEL`` is + "CHMF-SAMPLER", decreasing this can drastically increase both compute time and + memory usage. SAMPLER_BUFFER_FACTOR: float, optional - The arrays for the halo sampler will have size of SAMPLER_BUFFER_FACTOR multiplied by the expected - number of halos in the box. Ideally this should be close to unity but one may wish to increase it to - test alternative scenarios + The arrays for the halo sampler will have size of SAMPLER_BUFFER_FACTOR + multiplied by the expected number of halos in the box. Ideally this should be + close to unity but one may wish to increase it to test alternative scenarios. N_COND_INTERP: int, optional The number of condition bins in the inverse CMF tables. N_PROB_INTERP: int, optional @@ -784,13 +841,14 @@ class SimulationOptions(InputStruct): The minimum log-probability of the inverse CMF tables. HALOMASS_CORRECTION: float, optional This provides a corrective factor to the mass-limited (SAMPLE_METHOD==0) sampling, which multiplies the - expected mass from a condition by this number. The default value of 0.9 is calibrated to the mass-limited - sampling on a timestep of ZPRIME_STEP_FACTOR=1.02. + expected mass from a condition by this number. The default value of 0.9 is + calibrated to the mass-limited sampling on a timestep of ZPRIME_STEP_FACTOR=1.02. If ZPRIME_STEP_FACTOR is increased, this value should be set closer to 1. - This factor is also used in the partition (SAMPLE_METHOD==2) sampler, dividing nu(M) of each sample drawn. + This factor is also used in the partition (SAMPLE_METHOD==2) sampler, dividing + nu(M) of each sample drawn. PARKINSON_G0: float, optional - Only used when SAMPLE_METHOD==3, sets the normalisation of the correction to the extended press-schecter - used in Parkinson et al. 2008. + Only used when SAMPLE_METHOD==3, sets the normalisation of the correction to the + extended press-schecter used in Parkinson et al. 2008. PARKINSON_y1: float, optional Only used when SAMPLE_METHOD==3, sets the index of the sigma power-law term of the correction to the extended Press-Schechter mass function used in Parkinson et al. 2008. @@ -808,29 +866,38 @@ class SimulationOptions(InputStruct): INITIAL_REDSHIFT : float, optional Initial redshift used to perturb field from DELTA_R_FACTOR: float, optional - The factor by which to decrease the size of the filter in DexM when creating halo catalogues. + The factor by which to decrease the size of the filter in DexM when creating + halo catalogues. DENSITY_SMOOTH_RADIUS: float, optional The radius of the smoothing kernel in Mpc. DEXM_OPTIMIZE_MINMASS: float, optional - The minimum mass of a halo for which to use the DexM optimization if DEXM_OPTIMIZE is True. + The minimum mass of a halo for which to use the DexM optimization if + DEXM_OPTIMIZE is True. DEXM_R_OVERLAP: float, optional - The factor by which to multiply the halo radius to determine the distance within which smaller halos are excluded. + The factor by which to multiply the halo radius to determine the distance within + which smaller halos are excluded. CORR_STAR : float, optional Self-correlation length used for updating halo properties. To model the - correlation in the SHMR between timesteps, we sample from a conditional bivariate gaussian - with correlation factor given by exp(-dz/CORR_STAR). This value is placed in SimulationOptions - since it is used in the halo sampler, and not in the ionization routines. + correlation in the SHMR between timesteps, we sample from a conditional + bivariate gaussian with correlation factor given by exp(-dz/CORR_STAR). This + value is placed in SimulationOptions since it is used in the halo sampler, and + not in the ionization routines. CORR_SFR : float, optional - Self-correlation length used for updating star formation rate, see "CORR_STAR" for details. + Self-correlation length used for updating star formation rate, see "CORR_STAR" + for details. CORR_LX : float, optional - Self-correlation length used for updating xray luminosity, see "CORR_STAR" for details. + Self-correlation length used for updating xray luminosity, see "CORR_STAR" for + details. K_MAX_FOR_CLASS: float, optional - Maximum wavenumber to run CLASS, in 1/Mpc. Becomes relevant only if matter_options.POWER_SPECTRUM = "CLASS". + Maximum wavenumber to run CLASS, in 1/Mpc. Becomes relevant only if + ``matter_options.POWER_SPECTRUM = "CLASS"``. MIN_XE_FOR_FCOLL_IN_TAUX: float, optional - Minimum global x_e value for which the collapsed fraction (f_coll) is evaluated in the tau_X integral (X-ray optical depth). - When x_e is above this threshold value, it is assumed that f_coll=0, in order to speed up the calculations. - For now, this parameter becomes relevant only when run_global_evolution is called, as it controls the runtime of this - function (higher values reduce the runtime, in expense of degraded precision). + Minimum global x_e value for which the collapsed fraction (f_coll) is evaluated + in the tau_X integral (X-ray optical depth). When x_e is above this threshold + value, it is assumed that f_coll=0, in order to speed up the calculations. + For now, this parameter becomes relevant only when run_global_evolution is + called, as it controls the runtime of this function (higher values reduce the + runtime, in expense of degraded precision). """ _DEFAULT_HIRES_TO_LOWRES_FACTOR: ClassVar[float] = 3 @@ -1015,9 +1082,11 @@ class AstroOptions(InputStruct): USE_X_RAY_HEATING : bool, optional Whether to include X-ray heating (useful for debugging). USE_CMB_HEATING : bool, optional - Whether to include CMB heating. (cf Eq.4 of Meiksin 2021, arxiv.org/abs/2105.14516) + Whether to include CMB heating. (cf Eq.4 of Meiksin 2021, + https://arxiv.org/abs/2105.14516) USE_LYA_HEATING : bool, optional - Whether to use Lyman-alpha heating. (cf Sec. 3 of Reis+2021, doi.org/10.1093/mnras/stab2089) + Whether to use Lyman-alpha heating. (cf Sec. 3 of Reis+2021, + https://doi.org/10.1093/mnras/stab2089) INHOMO_RECO : bool, optional Whether to perform inhomogeneous recombinations. Increases the computation time. @@ -1026,55 +1095,87 @@ class AstroOptions(InputStruct): Dramatically increases the computation time. M_MIN_in_Mass : bool, optional Whether the minimum halo mass (for ionization) is defined by - mass or virial temperature. Only has an effect when SOURCE_MODEL == 'CONST-ION-EFF' + mass or virial temperature. Only has an effect when + ``SOURCE_MODEL == 'CONST-ION-EFF'`` PHOTON_CONS_TYPE : str, optional Whether to perform a small correction to account for the inherent photon non-conservation. This can be one of three types of correction: - no-photoncons: No photon cosnervation correction, - z-photoncons: Photon conservation correction by adjusting the redshift of the N_ion source field (Park+22) - alpha-photoncons: Adjustment to the escape fraction power-law slope, based on fiducial results in Park+22, This runs a - series of global xH evolutions and one calibration simulation to find the adjustment as a function of xH - f-photoncons: Adjustment to the escape fraction normalisation, runs one calibration simulation to find the - adjustment as a function of xH where f'/f = xH_global/xH_calibration + * no-photoncons: No photon cosnervation correction, + * z-photoncons: Photon conservation correction by adjusting the redshift of the + N_ion source field (Park+22) + * alpha-photoncons: Adjustment to the escape fraction power-law slope, based on + fiducial results in Park+22, This runs a series of global xH evolutions and + one calibration simulation to find the adjustment as a function of xH + * f-photoncons: Adjustment to the escape fraction normalisation, runs one + calibration simulation to find the adjustment as a function of xH where + ``f'/f = xH_global/xH_calibration`` FIX_VCB_AVG: bool, optional - Determines whether to use a fixed vcb=VAVG (*regardless* of USE_RELATIVE_VELOCITIES). It includes the average effect of velocities but not its fluctuations. See Muñoz+21 (2110.13919). + Determines whether to use a fixed vcb=VAVG (*regardless* of + USE_RELATIVE_VELOCITIES). It includes the average effect of velocities but not + its fluctuations. See `Muñoz+21 `_) USE_EXP_FILTER: bool, optional - Use the exponential filter (MFP-epsilon(r) from Davies & Furlanetto 2021) when calculating ionising emissivity fields - NOTE: this does not affect other field filters, and should probably be used with HII_FILTER==0 (real-space top-hat) - CELL_RECOMB: bool, optional - An alternate way of counting recombinations based on the local cell rather than the filter region. - This is part of the perspective shift (see Davies & Furlanetto 2021) from counting photons/atoms in a sphere and flagging a central - pixel to counting photons which we expect to reach the central pixel, and taking the ratio of atoms in the pixel. - This flag simply turns off the filtering of N_rec grids, and takes the recombinations in the central cell. + Use the exponential filter (MFP-epsilon(r) from Davies & Furlanetto 2021) when + calculating ionising emissivity fields. + + .. warning: this does not affect other field filters, and should probably be + used with HII_FILTER==0 (real-space top-hat). + + CELL_RECOMB : bool, optional + An alternate way of counting recombinations based on the local cell rather than + the filter region. This is part of the perspective shift (see Davies & + Furlanetto 2021) from counting photons/atoms in a sphere and flagging a central + pixel to counting photons which we expect to reach the central pixel, and taking + the ratio of atoms in the pixel. This flag simply turns off the filtering of + N_rec grids, and takes the recombinations in the central cell. + LYA_MULTIPLE_SCATTERING: bool, optional + If True, multiple scattering window function is used for the computation of + Lyman alpha photons. If False, the straight-line window function is used (see + more info in `2601.14360 `_). + This feature can be turned on only when the source model is defined on the + Lagrangian grid, otherwise an error is raised. + USE_ADIABATIC_FLUCTUATIONS: bool, optional + Whether to apply adiabatic fluctuations to the initial temperature box, see + Munoz 2023. If set to False, the initial temperature box is completely + homogeneous. Default is True. USE_UPPER_STELLAR_TURNOVER: bool, optional - Whether to use an additional powerlaw in stellar mass fraction at high halo mass. The pivot mass scale and power-law index are - controlled by two parameters, UPPER_STELLAR_TURNOVER_MASS and UPPER_STELLAR_TURNOVER_INDEX respectively. + Whether to use an additional powerlaw in stellar mass fraction at high halo + mass. The pivot mass scale and power-law index are controlled by two parameters, + UPPER_STELLAR_TURNOVER_MASS and UPPER_STELLAR_TURNOVER_INDEX respectively. This is currently only implemented using the discrete halo model, and has no effect otherwise. HALO_SCALING_RELATIONS_MEDIAN: bool, optional - If True, halo scaling relation parameters (F_STAR10,t_STAR etc...) define the median of their conditional distributions - If False, they describe the mean. - This becomes important when using non-symmetric dristributions such as the log-normal - HII_FILTER : string + If True, halo scaling relation parameters (F_STAR10,t_STAR etc...) define the + median of their conditional distributions. If False, they describe the mean. + This becomes important when using non-symmetric dristributions such as the log-normal. + HII_FILTER : str Filter for the halo or density field used to generate ionization field Available options are: 'spherical-tophat', 'sharp-k', and 'gaussian' - HEAT_FILTER : int + HEAT_FILTER : str Filter for the halo or density field used to generate the spin-temperature field Available options are: 'spherical-tophat', 'sharp-k', and 'gaussian' IONISE_ENTIRE_SPHERE: bool, optional - If True, ionises the entire sphere on the filter scale when an ionised region is found - in the excursion set. + If True, ionises the entire sphere on the filter scale when an ionised region is + found in the excursion set. INTEGRATION_METHOD_ATOMIC: str, optional - The integration method to use for conditional MF integrals of atomic halos in the grids: - NOTE: global integrals will use GSL QAG adaptive integration - 'GSL-QAG': GSL QAG adaptive integration, - 'GAUSS-LEGENDRE': Gauss-Legendre integration, previously forced in the interpolation tables, - 'GAMMA-APPROX': Approximate integration, assuming sharp cutoffs and a triple power-law for sigma(M) based on EPS + The integration method to use for conditional MF integrals of atomic halos in + the grids: + + * 'GSL-QAG': GSL QAG adaptive integration, + * 'GAUSS-LEGENDRE': Gauss-Legendre integration, previously forced in the + interpolation tables, + * 'GAMMA-APPROX': Approximate integration, assuming sharp cutoffs and a triple + power-law for sigma(M) based on EPS + + .. note:: Global integrals will use GSL QAG adaptive integration INTEGRATION_METHOD_MINI: str, optional - The integration method to use for conditional MF integrals of minihalos in the grids: - 'GSL-QAG': GSL QAG adaptive integration, - 'GAUSS-LEGENDRE': Gauss-Legendre integration, previously forced in the interpolation tables, - 'GAMMA-APPROX': Approximate integration, assuming sharp cutoffs and a triple power-law for sigma(M) based on EPS + The integration method to use for conditional MF integrals of minihalos in the + grids: + + * 'GSL-QAG': GSL QAG adaptive integration, + * 'GAUSS-LEGENDRE': Gauss-Legendre integration, previously forced in the + interpolation tables, + * 'GAMMA-APPROX': Approximate integration, assuming sharp cutoffs and a triple + power-law for sigma(M) based on EPS """ USE_MINI_HALOS: bool = field(default=False, converter=bool) @@ -1086,6 +1187,8 @@ class AstroOptions(InputStruct): FIX_VCB_AVG: bool = field(default=False, converter=bool) USE_EXP_FILTER: bool = field(default=True, converter=bool) CELL_RECOMB: bool = field(default=True, converter=bool) + LYA_MULTIPLE_SCATTERING = field(default=False, converter=bool) + USE_ADIABATIC_FLUCTUATIONS = field(default=True, converter=bool) PHOTON_CONS_TYPE: Literal[ "no-photoncons", "z-photoncons", "alpha-photoncons", "f-photoncons" ] = choice_field( @@ -1143,29 +1246,29 @@ class AstroParams(InputStruct): """ Astrophysical parameters. - NB: All Mean scaling relations are defined in log-space, such that the lines they produce - give exp(), this means that increasing the lognormal scatter in these relations - will increase the but not + .. attention: All Mean scaling relations are defined in log-space, such that the + lines they produce give ``exp(mean(log(property)))``. This means that increasing + the lognormal scatter in these relations will increase the ``mean(property)`` but + not ``mean(log(property))`` Parameters ---------- INHOMO_RECO : bool, optional - Whether inhomogeneous recombinations are being calculated. This is not a part of the - astro parameters structure, but is required by this class to set some default behaviour. + Whether inhomogeneous recombinations are being calculated. This is not a part of + the astro parameters structure, but is required by this class to set some + default behaviour. HII_EFF_FACTOR : float, optional The ionizing efficiency of high-z galaxies (zeta, from Eq. 2 of Greig+2015). Higher values tend to speed up reionization. F_STAR10 : float, optional The fraction of galactic gas in stars for 10^10 solar mass haloes. - Only used in the "new" parameterization, - This is used along with `F_ESC10` to determine `HII_EFF_FACTOR` (which - is then unused). See Eq. 11 of Greig+2018 and Sec 2.1 of Park+2018. + Only used if ``MASS_DEPENDENT_ZETA`` is True in :class:`AstroOptions`. + See Eq. 11 of Greig+2018 and Sec 2.1 of Park+2018. Given in log10 units. F_STAR7_MINI : float, optional The fraction of galactic gas in stars for 10^7 solar mass minihaloes. Only used in the "minihalo" parameterization, i.e. when `USE_MINI_HALOS` is set to True - (in :class:`AstroOptions`). If so, this is used along with `F_ESC7_MINI` to - determine `HII_EFF_FACTOR_MINI` (which is then unused). See Eq. 8 of Qin+2020. + (in :class:`AstroOptions`).. See Eq. 8 of Qin+2020. If the MCG scaling relations are not provided explicitly, we extend the ACG ones by default. Given in log10 units. ALPHA_STAR : float, optional @@ -1179,20 +1282,21 @@ class AstroParams(InputStruct): Lognormal scatter (dex) of the halo mass to stellar mass relation. Uniform across all masses and redshifts. SIGMA_SFR_LIM : float, optional - Lognormal scatter (dex) of the stellar mass to SFR relation above a stellar mass of 1e10 solar. + Lognormal scatter (dex) of the stellar mass to SFR relation above a stellar mass + of 1e10 solar. SIGMA_SFR_INDEX : float, optional index of the power-law between SFMS scatter and stellar mass below 1e10 solar. F_ESC10 : float, optional The "escape fraction", i.e. the fraction of ionizing photons escaping into the - IGM, for 10^10 solar mass haloes. Only used in the "new" parameterization. - This is used along with `F_STAR10` to determine `HII_EFF_FACTOR` (which + IGM, for 10^10 solar mass haloes. Only used if ``MASS_DEPENDENT_ZETA`` is True + in :class:`AstroOptions`. This is used along with `F_STAR10` to determine + ``HII_EFF_FACTOR`` (which is then unused). See Eq. 11 of Greig+2018 and Sec 2.1 of Park+2018. F_ESC7_MINI: float, optional The "escape fraction for minihalos", i.e. the fraction of ionizing photons escaping into the IGM, for 10^7 solar mass minihaloes. Only used in the "minihalo" parameterization, i.e. when `USE_MINI_HALOS` is set to True (in - :class:`AstroOptions`). If so, this is used along with `F_ESC7_MINI` to determine - `HII_EFF_FACTOR_MINI` (which is then unused). See Eq. 17 of Qin+2020. If the MCG + :class:`AstroOptions`). See Eq. 17 of Qin+2020. If the MCG scaling relations are not provided explicitly, we extend the ACG ones by default. Given in log10 units. ALPHA_ESC : float, optional @@ -1204,18 +1308,18 @@ class AstroParams(InputStruct): See Sec 2.1 of Park+2018. R_BUBBLE_MAX : float, optional Mean free path in Mpc of ionizing photons within ionizing regions (Sec. 2.1.2 of - Greig+2015). Default is 50 if `INHOMO_RECO` is True, or 15.0 if not. + Greig+2015). Default is 50 if ``INHOMO_RECO`` is True, or 15.0 if not. ION_Tvir_MIN : float, optional Minimum virial temperature of star-forming haloes (Sec 2.1.3 of Greig+2015). Given in log10 units. L_X : float, optional The specific X-ray luminosity per unit star formation escaping host galaxies. - Cf. Eq. 6 of Greig+2018. Given in log10 units. For the double power-law used in the Halo Model - This gives the low-z limite. + Cf. Eq. 6 of Greig+2018. Given in log10 units. For the double power-law used in + the Halo Model. This gives the low-z limit. L_X_MINI: float, optional The specific X-ray luminosity per unit star formation escaping host galaxies for minihalos. Cf. Eq. 23 of Qin+2020. Given in log10 units. For the double - power-law used in the Halo Model. This gives the low-z limite. If the MCG + power-law used in the Halo Model. This gives the low-z limit. If the MCG scaling relations are not provided explicitly, we extend the ACG ones by default. NU_X_THRESH : float, optional X-ray energy threshold for self-absorption by host galaxies (in eV). Also called @@ -1234,30 +1338,40 @@ class AstroParams(InputStruct): Fractional characteristic time-scale (fraction of hubble time) defining the star-formation rate of galaxies. See Sec 2.1, Eq. 3 of Park+2018. A_LW, BETA_LW: float, optional - Impact of the LW feedback on Mturn for minihaloes. Default is 22.8685 and 0.47 following Machacek+01, respectively. Latest simulations suggest 2.0 and 0.6. See Sec 2 of Muñoz+21 (2110.13919). + Impact of the LW feedback on Mturn for minihaloes. Default is 22.8685 and 0.47 + following Machacek+01, respectively. Latest simulations suggest 2.0 and 0.6. + See Sec 2 of Muñoz+21 (2110.13919). A_VCB, BETA_VCB: float, optional - Impact of the DM-baryon relative velocities on Mturn for minihaloes. Default is 1.0 and 1.8, and agrees between different sims. See Sec 2 of Muñoz+21 (2110.13919). - UPPER_STELLAR_TURNOVER_MASS: - The pivot mass associated with the optional upper mass power-law of the stellar-halo mass relation - (see AstroOptions.USE_UPPER_STELLAR_TURNOVER) - UPPER_STELLAR_TURNOVER_INDEX: - The power-law index associated with the optional upper mass power-law of the stellar-halo mass relation - (see AstroOptions.USE_UPPER_STELLAR_TURNOVER) + Impact of the DM-baryon relative velocities on Mturn for minihaloes. Default is + 1.0 and 1.8, and agrees between different sims. See Sec 2 of Muñoz+21 (2110.13919). + UPPER_STELLAR_TURNOVER_MASS : float + The pivot mass associated with the optional upper mass power-law of the + stellar-halo mass relation (see ``USE_UPPER_STELLAR_TURNOVER`` + in :class:`AstroOptions`) + UPPER_STELLAR_TURNOVER_INDEX : float + The power-law index associated with the optional upper mass power-law of the + stellar-halo mass relation (see ``USE_UPPER_STELLAR_TURNOVER`` + in :class:`AstroOptions`) SIGMA_LX: float, optional - Lognormal scatter (dex) of the Xray luminosity relation (a function of stellar mass, star formation rate and redshift). - This scatter is uniform across all halo properties and redshifts. + Lognormal scatter (dex) of the Xray luminosity relation (a function of stellar + mass, star formation rate and redshift). This scatter is uniform across all + halo properties and redshifts. FIXED_VAVG : float, optional - The fixed value of the average velocity used when AstroOptions.FIX_VCB_AVG is set to True. + The fixed value of the average velocity used when AstroOptions.FIX_VCB_AVG is + set to True. POP2_ION: float, optional Number of ionizing photons per baryon produced by Pop II stars. POP3_ION: float, optional Number of ionizing photons per baryon produced by Pop III stars. CLUMPING_FACTOR: float, optional - Clumping factor of the IGM used ONLY in the x-ray partial ionisations (not the reionsiation model). Default is 2.0. + Clumping factor of the IGM used ONLY in the x-ray partial ionisations (not the + reionsiation model). Default is 2.0. ALPHA_UVB: float, optional - The power-law index of the UVB spectrum. Used for Gamma12 in the recombination model + The power-law index of the UVB spectrum. Used for Gamma12 in the recombination + model DELTA_R_HII_FACTOR: float, optional - The factor by which to decrease the size of the HII filter when calculating the HII regions. + The factor by which to decrease the size of the HII filter when calculating the + HII regions. R_BUBBLE_MIN: float, optional Minimum size of ionized regions in Mpc. Default is 0.620350491. MAX_DVDR: float, optional @@ -1443,7 +1557,7 @@ def get_logspaced_redshifts( min_redshift: float, z_step_factor: float, max_redshift: float, -) -> tuple[float]: +) -> tuple[float, ...]: """Compute a sequence of redshifts to evolve over that are log-spaced.""" redshifts = [min_redshift] while redshifts[-1] < max_redshift: @@ -1452,7 +1566,7 @@ def get_logspaced_redshifts( return tuple(redshifts[::-1]) -def _node_redshifts_converter(value) -> tuple[float] | None: +def _node_redshifts_converter(value) -> tuple[float, ...] | None: if value is None or len(value) == 0: return () if hasattr(value, "__len__"): @@ -1470,27 +1584,27 @@ class InputParameters: Parameters ---------- - random_seed + random_seed : float The seed that determines the realization produced by a 21cmFAST run. - node_redshifts + node_redshifts : list[float] The redshifts at which coeval boxes will be computed. By default, empty if no evolution is required, and logarithmically spaced in (1+z) between z=5.5 and SimulationOptions.Z_HEAT_MAX if evolution is required. - cosmo_params + cosmo_params : :class:`~CosmoParams` Cosmological parameters of a 21cmFAST run. - simulation_options + simulation_options : :class:`~SimulationOptions` Parameters controlling the simulation as a whole, e.g. the box size and dimensionality. - matter_options + matter_options : :class:`~MatterOptions` Parameters controlling the matter field generated by 21cmFAST. - astro_options + astro_options : :class:`~AstroOptions` Options for which physical processes to include in the simulation. - astro_params + astro_params : :class:`~AstroParams` Astrophysical parameter values. """ - random_seed = _field(converter=int) + random_seed: int = _field(converter=int) cosmo_params: CosmoParams = input_param_field(CosmoParams) matter_options: MatterOptions = input_param_field(MatterOptions) simulation_options: SimulationOptions = input_param_field(SimulationOptions) @@ -1629,6 +1743,10 @@ def _astro_options_validator(self, att, val): raise ValueError( f"USE_EXP_FILTER is not compatible with SOURCE_MODEL == {self.matter_options.SOURCE_MODEL}" ) + if val.LYA_MULTIPLE_SCATTERING: + raise ValueError( + f"LYA_MULTIPLE_SCATTERING is not compatible with SOURCE_MODEL == {self.matter_options.SOURCE_MODEL}" + ) if ( not self.matter_options.has_discrete_halos and val.USE_UPPER_STELLAR_TURNOVER @@ -1748,7 +1866,7 @@ def is_compatible_with(self, other: Self) -> bool: for key in self.merge_keys() ) - def evolve_input_structs(self, **kwargs): + def evolve_input_structs(self, **kwargs) -> Self: """Return an altered clone of the `InputParameters` structs. Unlike clone(), this function takes fields from the constituent `InputStruct` classes @@ -1802,7 +1920,7 @@ def from_template( random_seed: int, node_redshifts: tuple[float] | None = None, **kwargs, - ): + ) -> Self: """Construct full InputParameters instance from native or TOML file template. Parameters @@ -1835,11 +1953,11 @@ def from_template( dct.pop("cosmo_tables") return cls(**dct, **cls_kw) - def clone(self, **kwargs): + def clone(self, **kwargs) -> Self: """Generate a copy of the InputParameter structure with specified changes.""" return evolve(self, **kwargs) - def __repr__(self): + def __repr__(self) -> str: """ Return a string representation of the structure. diff --git a/src/py21cmfast/wrapper/outputs.py b/src/py21cmfast/wrapper/outputs.py index af4d44765..81e74bd20 100644 --- a/src/py21cmfast/wrapper/outputs.py +++ b/src/py21cmfast/wrapper/outputs.py @@ -154,14 +154,14 @@ def arrays(self) -> dict[str, Array]: return {k: x for k, x in me.items() if isinstance(x, Array)} @cached_property - def struct(self) -> StructWrapper: + def _struct(self) -> StructWrapper: """The python-wrapped struct associated with this input object.""" return StructWrapper(self._name) @cached_property - def cstruct(self) -> StructWrapper: + def _cstruct(self) -> StructWrapper: """The object pointing to the memory accessed by C-code for this struct.""" - return self.struct.cstruct + return self._struct.cstruct def _init_arrays(self): for k, array in self.arrays.items(): @@ -193,11 +193,11 @@ def push_to_backend(self): # to unnecessarily load things in. We leave it to the user to ensure that all # required arrays are loaded into memory before calling this function. if array.state.initialized: - self.struct.expose_to_c(array, name) + self._struct.expose_to_c(array, name) - for k in self.struct.primitive_fields: + for k in self._struct.primitive_fields: if getattr(self, k) is not None: - setattr(self.cstruct, k, getattr(self, k)) + setattr(self._cstruct, k, getattr(self, k)) def pull_from_backend(self): """Sync the current state of the object with the underlying C-struct. @@ -205,8 +205,8 @@ def pull_from_backend(self): This will pull any primitives calculated in the backend to the python object. Arrays are passed in as pointers, and do not need to be copied back. """ - for k in self.struct.primitive_fields: - setattr(self, k, getattr(self.cstruct, k)) + for k in self._struct.primitive_fields: + setattr(self, k, getattr(self._cstruct, k)) def get(self, ary: str | Array): """If possible, load an array from disk, storing it and returning the underlying array.""" @@ -309,7 +309,7 @@ def _remove_array(self, k: str, *, force=False): return if state.c_has_active_memory: - lib.free(getattr(self.cstruct, k)) + lib.free(getattr(self._cstruct, k)) setattr(self, k, array.without_value()) @@ -407,7 +407,7 @@ def summarize(self, indent: int = 0) -> str: # print primitive fields out += "".join( f"{indent} {fieldname:>25}: {getattr(self, fieldname, 'non-existent')}\n" - for fieldname in self.struct.primitive_fields + for fieldname in self._struct.primitive_fields ) return out @@ -454,7 +454,7 @@ def _compute(self, allow_already_computed: bool = False, *args): # Construct the args. All StructWrapper objects need to actually pass their # underlying cstruct, rather than themselves. inputs = [ - arg.cstruct if isinstance(arg, OutputStruct | InputStruct) else arg + arg._cstruct if isinstance(arg, OutputStruct | InputStruct) else arg for arg in args ] # Sync the python/C memory @@ -471,7 +471,7 @@ def _compute(self, allow_already_computed: bool = False, *args): # Perform the C computation try: - exitcode = self._c_compute_function(*inputs, self.cstruct) + exitcode = self._c_compute_function(*inputs, self._cstruct) except TypeError as e: logger.error(f"Arguments to {self._c_compute_function.__name__}: {inputs}") raise e @@ -1146,6 +1146,8 @@ class XraySourceBox(OutputStructZ): filtered_sfr = _arrayfield() filtered_sfr_mini = _arrayfield(optional=True) filtered_xray = _arrayfield() + filtered_sfr_lw = _arrayfield(optional=True) + filtered_sfr_mini_lw = _arrayfield(optional=True) mean_sfr = _arrayfield() mean_sfr_mini = _arrayfield(optional=True) mean_log10_Mcrit_LW = _arrayfield(optional=True) @@ -1190,6 +1192,9 @@ def new(cls, inputs: InputParameters, redshift: float, **kw) -> Self: out["mean_log10_Mcrit_LW"] = Array( (inputs.astro_params.N_STEP_TS,), dtype=np.float64 ) + if inputs.astro_options.LYA_MULTIPLE_SCATTERING: + out["filtered_sfr_lw"] = Array(shape, dtype=np.float32) + out["filtered_sfr_mini_lw"] = Array(shape, dtype=np.float32) return cls( inputs=inputs, @@ -1216,6 +1221,7 @@ def compute( R_inner, R_outer, R_ct, + R_star, allow_already_computed: bool = False, ): """Compute the function.""" @@ -1225,6 +1231,7 @@ def compute( R_inner, R_outer, R_ct, + R_star, ) @@ -1270,6 +1277,7 @@ def new(cls, inputs: InputParameters, redshift: float, **kw) -> Self: } if inputs.astro_options.USE_MINI_HALOS: out["J_21_LW"] = Array(shape, dtype=np.float32) + return cls(inputs=inputs, redshift=redshift, **out, **kw) @cached_property diff --git a/tests/io/test_caching.py b/tests/io/test_caching.py index 855de4f08..94fcb3f0e 100644 --- a/tests/io/test_caching.py +++ b/tests/io/test_caching.py @@ -38,7 +38,7 @@ def create_full_run_cache(cachedir: Path) -> caching.RunCache: setattr(o, k, v.with_value(v.value)) # Mock the primitive fields as well... - for fld in o.struct.primitive_fields: + for fld in o._struct.primitive_fields: setattr(o, fld, 0.0) h5.write_output_to_hdf5(o, fname) diff --git a/tests/produce_integration_test_data.py b/tests/produce_integration_test_data.py index b7337b646..575ee3028 100644 --- a/tests/produce_integration_test_data.py +++ b/tests/produce_integration_test_data.py @@ -241,6 +241,30 @@ }, ], "fftw_wisdom": [18, {"USE_FFTW_WISDOM": True}], + "multiple_scattering": [ + 18, + { + "LYA_MULTIPLE_SCATTERING": True, + "SOURCE_MODEL": "L-INTEGRAL", + "USE_TS_FLUCT": True, + }, + ], + "multiple_scattering_mini": [ + 18, + { + "LYA_MULTIPLE_SCATTERING": True, + "SOURCE_MODEL": "L-INTEGRAL", + "USE_TS_FLUCT": True, + "USE_MINI_HALOS": True, + "INHOMO_RECO": True, + "N_THREADS": 4, + "USE_RELATIVE_VELOCITIES": True, + "POWER_SPECTRUM": "CLASS", + "K_MAX_FOR_CLASS": 1.0, + "R_BUBBLE_MAX": 50.0, + "M_TURN": 5.0, + }, + ], } if len(set(OPTIONS_TESTRUNS.keys())) != len(list(OPTIONS_TESTRUNS.keys())): diff --git a/tests/test_data/power_spectra_multiple_scattering.h5 b/tests/test_data/power_spectra_multiple_scattering.h5 new file mode 100644 index 000000000..b1e20513f Binary files /dev/null and b/tests/test_data/power_spectra_multiple_scattering.h5 differ diff --git a/tests/test_data/power_spectra_multiple_scattering_mini.h5 b/tests/test_data/power_spectra_multiple_scattering_mini.h5 new file mode 100644 index 000000000..34b7eee20 Binary files /dev/null and b/tests/test_data/power_spectra_multiple_scattering_mini.h5 differ diff --git a/tests/test_data_exists.py b/tests/test_data_exists.py index 82c578213..17c481a2a 100644 --- a/tests/test_data_exists.py +++ b/tests/test_data_exists.py @@ -2,20 +2,20 @@ import pytest -from py21cmfast import DATA_PATH +from py21cmfast import _DATA_PATH @pytest.fixture(scope="module") def datafiles(): - return list(DATA_PATH.glob("*.dat")) + return list(_DATA_PATH.glob("*.dat")) def test_exists(datafiles): for fl in datafiles: assert fl.exists() - assert (DATA_PATH / "x_int_tables").exists() - assert (DATA_PATH / "x_int_tables").is_dir() + assert (_DATA_PATH / "x_int_tables").exists() + assert (_DATA_PATH / "x_int_tables").is_dir() def test_readable(datafiles): diff --git a/tests/test_filtering.py b/tests/test_filtering.py index d6242675a..d8720b7e4 100644 --- a/tests/test_filtering.py +++ b/tests/test_filtering.py @@ -1,6 +1,7 @@ """Test filtering of the density field.""" import matplotlib as mpl +import mpmath import numpy as np import pytest from matplotlib.colors import Normalize @@ -53,8 +54,8 @@ def get_expected_output_centre(r_in, R_filter, R_param, filter_flag): return (R_ratio < 1) * np.exp(-r_in / R_param) / exp_vol elif filter_flag == 4: # output is spherical shell - exp_vol = 4 / 3 * np.pi * (R_filter**3 - R_param**3) - return (R_ratio < 1) * (R_param <= r_in) / exp_vol + exp_vol = 4 / 3 * np.pi * (R_param**3 - R_filter**3) + return (R_ratio > 1) * (R_param >= r_in) / exp_vol # return binned quantities @@ -100,7 +101,7 @@ def test_filters(filter_flag, R, plt): if filter_flag == 3: R_param = 20 elif filter_flag == 4: - R_param = max(R - 4 * (up.BOX_LEN / up.HII_DIM), 0) + R_param = R + 4 * (up.BOX_LEN / up.HII_DIM) else: R_param = 0 @@ -109,6 +110,7 @@ def test_filters(filter_flag, R, plt): ffi.cast("float *", input_box_centre.ctypes.data), R, R_param, + 0.0, filter_flag, ffi.cast("double *", output_box_centre.ctypes.data), ) @@ -298,3 +300,79 @@ def filter_plot( axs[idx, 3].grid() axs[idx, 3].set_xlabel("expected cell value") axs[idx, 3].set_ylabel("filtered cell value") + + +@pytest.mark.parametrize("R_inner", R_PARAM_LIST) +@pytest.mark.parametrize("n", [2, 4, 6, 8]) +@pytest.mark.parametrize("R_star", [1e-6, 5, 10, 20]) +def test_MS_filter(R_inner, n, R_star): + """Test the multiple scattering filter.""" + opts = prd.get_all_options_struct(redshift=10.0) + inputs = opts["inputs"] + up = inputs.simulation_options + + # testing a random input box + rng = np.random.default_rng(12345) + input_box = rng.random((up.HII_DIM,) * 3) + output_box_SL = np.zeros((up.HII_DIM,) * 3, dtype="f8") + output_box_MS = np.zeros((up.HII_DIM,) * 3, dtype="f8") + + R_outer = n * R_inner + + broadcast_input_struct(inputs) + lib.test_filter( + ffi.cast("float *", input_box.ctypes.data), + R_inner, + R_outer, + R_star, + 4, + ffi.cast("double *", output_box_SL.ctypes.data), + ) + + broadcast_input_struct(inputs) + lib.test_filter( + ffi.cast("float *", input_box.ctypes.data), + R_inner, + R_outer, + R_star, + 5, + ffi.cast("double *", output_box_MS.ctypes.data), + ) + + if R_star < 1: + # Test that MS filter returns the same output as the SL filter in the SL limit (R_star -> 0). + np.testing.assert_allclose(output_box_SL, output_box_MS, atol=1e-4) + else: + # Test that the SL and MS filters yield different outputs if R_star is not too small + assert not np.allclose(output_box_SL, output_box_MS, atol=1e-4) + + +@pytest.mark.parametrize("x_em", [0.0, 0.1, 0.5, 1.0, 5.0, 10.0, 50.0, 100.0, 500.0]) +def test_hyper_2F3(x_em): + """Test the implementation of the hypergeometric function in the C code.""" + mu = lib.compute_mu_for_multiple_scattering(x_em) + eta = lib.compute_eta_for_multiple_scattering(x_em) + + # Eq. (28) in arxiv: 2601.14360 (mu = alpha/(alpha+beta), eta = alpha/(alpha+beta^2)) + if mu == 0.0 and eta == 0.0: + alpha = np.inf + beta = 0.0 + else: + alpha = (1.0 / eta - 1.0) / pow(1.0 / mu - 1.0, 2) + beta = (1.0 / eta - 1.0) / (1.0 / mu - 1.0) + + kR_array = np.logspace(-1, 3, 100) + hyper_python = np.array( + [ + float( + mpmath.hyper( + [(2.0 + alpha) / 2.0, (3.0 + alpha) / 2.0], + [5.0 / 2.0, (2.0 + alpha + beta) / 2.0, (3.0 + alpha + beta) / 2.0], + -0.25 * kR**2, + ) + ) + for kR in kR_array + ] + ) + hyper_C = np.array([lib.hyper_2F3(kR, alpha, beta) for kR in kR_array]) + np.testing.assert_allclose(hyper_python, hyper_C, rtol=0.0, atol=2e-3) diff --git a/tests/test_input_structs.py b/tests/test_input_structs.py index 798e2d7da..9db1c7203 100644 --- a/tests/test_input_structs.py +++ b/tests/test_input_structs.py @@ -98,7 +98,7 @@ def test_pickle(self): assert self.cosmo == c4 # Make sure the c data gets loaded fine. - assert c4.cstruct.hlittle == self.cosmo.cstruct.hlittle + assert c4._cstruct.hlittle == self.cosmo._cstruct.hlittle def test_asdict(self): """Test the asdict() method works.""" @@ -435,6 +435,20 @@ class TestInputParameters: "astro_options": AstroOptions(PHOTON_CONS_TYPE="z-photoncons"), }, ), + ( + ValueError, + "LYA_MULTIPLE_SCATTERING is not compatible with SOURCE_MODEL == E-INTEGRAL", + { + "matter_options": MatterOptions( + SOURCE_MODEL="E-INTEGRAL", + ), + "astro_options": AstroOptions( + LYA_MULTIPLE_SCATTERING=True, + USE_EXP_FILTER=False, + USE_UPPER_STELLAR_TURNOVER=False, + ), + }, + ), ( ValueError, "USE_EXP_FILTER is not compatible with SOURCE_MODEL == E-INTEGRAL", diff --git a/tests/test_singlefield.py b/tests/test_singlefield.py index 9fd36ca06..2bf4043df 100644 --- a/tests/test_singlefield.py +++ b/tests/test_singlefield.py @@ -20,6 +20,7 @@ PerturbedField, TsBox, ) +from py21cmfast.wrapper.arrays import Array @pytest.fixture(scope="module") @@ -541,3 +542,54 @@ def test_bad_input_structs(default_input_struct_ts): halobox=hb, previous_ionized_box=ib_p, ) + + +@pytest.mark.parametrize("lya_multiple_scattering", [False, True]) +@pytest.mark.parametrize("use_mini_halos", [False, True]) +def test_xray_source_field_with_zero_sfr( + default_input_struct_ts, redshift, use_mini_halos, lya_multiple_scattering +): + """Test compute_xray_source_field with zero sfr boxes.""" + inputs = default_input_struct_ts.evolve_input_structs( + USE_MINI_HALOS=use_mini_halos, + INHOMO_RECO=True, + LYA_MULTIPLE_SCATTERING=lya_multiple_scattering, + SOURCE_MODEL="L-INTEGRAL", + ) + ics = p21c.InitialConditions.new(inputs=inputs) + + hbox1 = HaloBox.new(redshift=redshift + 1, inputs=inputs) + hbox2 = HaloBox.new(redshift=redshift, inputs=inputs) + + # This is needed because the input arrays must be in a computed state. + fields = ["halo_sfr", "halo_xray"] + if use_mini_halos: + fields += ["halo_sfr_mini", "log10_Mcrit_MCG_ave"] + shape = hbox1.halo_sfr.shape + array = ( + Array(shape=shape, dtype=np.float32) + .initialize() + .with_value(val=np.zeros(shape)) + ) + for hbox in [hbox1, hbox2]: + for name in fields: + setattr(hbox, name, array.computed()) + if use_mini_halos: + hbox.log10_Mcrit_MCG_ave = 5.0 + + xraysource = p21c.compute_xray_source_field( + initial_conditions=ics, + hboxes=[hbox1, hbox2], + redshift=redshift, + ) + + output_fields = ["filtered_sfr", "filtered_xray"] + if use_mini_halos: + output_fields += [ + "filtered_sfr_mini", + ] + if lya_multiple_scattering: + output_fields += ["filtered_sfr_lw", "filtered_sfr_mini_lw"] + + for field in output_fields: + assert np.all(getattr(xraysource, field).value == 0.0) diff --git a/tests/test_utils.py b/tests/test_utils.py index bcb176e32..a8cd0145c 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -9,34 +9,80 @@ def test_ref_printing(): inputs = InputParameters.from_template("latest", random_seed=1234) - ref_str = show_references(inputs, print_to_stdout=False) - - assert "2011MNRAS.411..955M" in ref_str - assert "10.21105/joss.02582" in ref_str - assert "10.1093/mnras/stz032" in ref_str - assert "10.1093/mnras/staa1131" not in ref_str - assert "10.1093/mnras/stac185" not in ref_str - assert "10.48550/arXiv.2504.17254" not in ref_str - - inputs = InputParameters.from_template("mini-dhalos", random_seed=1234) - ref_str = show_references(inputs, print_to_stdout=False) - - assert "2011MNRAS.411..955M" in ref_str - assert "10.21105/joss.02582" in ref_str - assert "10.1093/mnras/stz032" in ref_str - assert "10.1093/mnras/staa1131" in ref_str - assert "10.1093/mnras/stac185" in ref_str - assert "10.48550/arXiv.2504.17254" in ref_str + ref_str = show_references(inputs, lightcone=True, print_to_stdout=False) + + assert "2011MNRAS.411..955M" in ref_str # 21cmFAST first paper + assert "10.21105/joss.02582" in ref_str # v3 (wrapper) + assert "10.1093/mnras/stu377" in ref_str # INHOMO_RECO + assert "10.1093/mnras/sty796" in ref_str # LIGHTCONE + RSD + assert "10.1093/mnras/stz032" in ref_str # USE_MASS_DEPENDENT_ZETA + assert "10.1093/mnras/staa1131" not in ref_str # USE_MINI_HALOS + assert "10.1093/mnras/stac185" not in ref_str # USE_RELATIVE_VELOCITIES + assert "10.1093/mnras/stac2756" not in ref_str # PHOTON_CONS + assert "10.1051/0004-6361/202554951" not in ref_str # LAGRANGIAN_SOURCE_MODEL + assert "10.48550/arXiv.2601.14360" not in ref_str # LYA_MULTIPLE_SCATTERING + + inputs = InputParameters.from_template( + "mini-dhalos", random_seed=1234, K_MAX_FOR_CLASS=1.0 + ) + ref_str = show_references(inputs, lightcone=True, print_to_stdout=False) + + assert "2011MNRAS.411..955M" in ref_str # 21cmFAST first paper + assert "10.21105/joss.02582" in ref_str # v3 (wrapper) + assert "10.1093/mnras/stu377" in ref_str # INHOMO_RECO + assert "10.1093/mnras/sty796" in ref_str # LIGHTCONE + RSD + assert "10.1093/mnras/stz032" in ref_str # USE_MASS_DEPENDENT_ZETA + assert "10.1093/mnras/staa1131" in ref_str # USE_MINI_HALOS + assert "10.1093/mnras/stac185" in ref_str # USE_RELATIVE_VELOCITIES + assert "10.1093/mnras/stac2756" not in ref_str # PHOTON_CONS + assert "10.1051/0004-6361/202554951" in ref_str # LAGRANGIAN_SOURCE_MODEL + assert "10.48550/arXiv.2601.14360" not in ref_str # LYA_MULTIPLE_SCATTERING inputs = InputParameters.from_template("const-zeta", random_seed=1234) - ref_str = show_references(inputs, print_to_stdout=True) - - assert "2011MNRAS.411..955M" in ref_str - assert "10.21105/joss.02582" in ref_str - assert "10.1093/mnras/stz032" not in ref_str - assert "10.1093/mnras/staa1131" not in ref_str - assert "10.1093/mnras/stac185" not in ref_str - assert "10.48550/arXiv.2504.17254" not in ref_str + ref_str = show_references(inputs, lightcone=True, print_to_stdout=True) + + assert "2011MNRAS.411..955M" in ref_str # 21cmFAST first paper + assert "10.21105/joss.02582" in ref_str # v3 (wrapper) + assert "10.1093/mnras/stu377" not in ref_str # INHOMO_RECO + assert "10.1093/mnras/sty796" in ref_str # LIGHTCONE + RSD + assert "10.1093/mnras/stz032" not in ref_str # USE_MASS_DEPENDENT_ZETA + assert "10.1093/mnras/staa1131" not in ref_str # USE_MINI_HALOS + assert "10.1093/mnras/stac185" not in ref_str # USE_RELATIVE_VELOCITIES + assert "10.1093/mnras/stac2756" not in ref_str # PHOTON_CONS + assert "10.1051/0004-6361/202554951" not in ref_str # LAGRANGIAN_SOURCE_MODEL + assert "10.48550/arXiv.2601.14360" not in ref_str # LYA_MULTIPLE_SCATTERING + + inputs = InputParameters.from_template( + "fixed-halos", random_seed=1234, LYA_MULTIPLE_SCATTERING=True + ) + ref_str = show_references(inputs, lightcone=True, print_to_stdout=True) + + assert "2011MNRAS.411..955M" in ref_str # 21cmFAST first paper + assert "10.21105/joss.02582" in ref_str # v3 (wrapper) + assert "10.1093/mnras/stu377" in ref_str # INHOMO_RECO + assert "10.1093/mnras/sty796" in ref_str # LIGHTCONE + RSD + assert "10.1093/mnras/stz032" in ref_str # USE_MASS_DEPENDENT_ZETA + assert "10.1093/mnras/staa1131" not in ref_str # USE_MINI_HALOS + assert "10.1093/mnras/stac185" not in ref_str # USE_RELATIVE_VELOCITIES + assert "10.1093/mnras/stac2756" not in ref_str # PHOTON_CONS + assert "10.1051/0004-6361/202554951" in ref_str # LAGRANGIAN_SOURCE_MODEL + assert "10.48550/arXiv.2601.14360" in ref_str # LYA_MULTIPLE_SCATTERING + + inputs = InputParameters.from_template( + "latest", random_seed=1234, PHOTON_CONS_TYPE="z-photoncons" + ) + ref_str = show_references(inputs, lightcone=False, print_to_stdout=True) + + assert "2011MNRAS.411..955M" in ref_str # 21cmFAST first paper + assert "10.21105/joss.02582" in ref_str # v3 (wrapper) + assert "10.1093/mnras/stu377" in ref_str # INHOMO_RECO + assert "10.1093/mnras/sty796" not in ref_str # LIGHTCONE + RSD + assert "10.1093/mnras/stz032" in ref_str # USE_MASS_DEPENDENT_ZETA + assert "10.1093/mnras/staa1131" not in ref_str # USE_MINI_HALOS + assert "10.1093/mnras/stac185" not in ref_str # USE_RELATIVE_VELOCITIES + assert "10.1093/mnras/stac2756" in ref_str # PHOTON_CONS + assert "10.1051/0004-6361/202554951" not in ref_str # LAGRANGIAN_SOURCE_MODEL + assert "10.48550/arXiv.2601.14360" not in ref_str # LYA_MULTIPLE_SCATTERING class TestRecursiveDifference: