meer21cm package

meer21cm.cosmology module

This module contains the class for storing the cosmological model used for calculation.

The class CosmologyParameters is the class for storing the cosmological parameters, and settings for computing matter power spectrum.

The class CosmologyCalculator is the base class for storing the cosmological model used for calculation. It is typically used as a base class for other classes that inherit from it, and not used directly.

Note that, there are always two sets of cosmological parameters defined in the class: - the fiducial cosmology, which is the cosmology that is used to transform sky coordinates to comoving coordinates. - the true cosmology, which is the cosmology that is used to compute the model power spectra.

The Alcock–Paczynski effect is then always included in the model power spectrum calculation as well as converting the field k-modes to model k-modes (see meer21cm.power.PowerSpectrum).

class meer21cm.cosmology.CosmologyCalculator(backend: str = 'camb', omega_hi: ndarray | float = 0.0005, fiducial_cosmology: str | dict = {'As': np.float64(2.105209331337507e-09), 'Neff': 3.046, 'h': np.float64(0.6766), 'neutrino_mass': np.float64(0.06), 'ns': 0.9665, 'omega_baryon': 0.04897, 'omega_cold': 0.30966, 'tau': 0.0561, 'w0': -1.0, 'wa': 0.0}, true_cosmology: str | dict | None = None, ps_type: str = 'linear', kmin: float = 0.001, kmax: float = 3.0, num_kpoints: int = 200, cold: bool = True, **params)[source]

Bases: Specification

The class for storing the cosmological model used for calculation.

The underlying cosmological model is defined via astropy.cosmology.w0waCDM with all the background properties calculated via astropy.

The matter density fluctuation is calculated using camb or baccoemu based on the input backend.

Parameters:

backend (str, default "camb") – The backend to use for computing the matter power spectrum. Either “camb” or “bacco”.
omega_hi (float or np.ndarray, default 5e-4) – The HI density as a function of redshift, over the critical density of the Universe at z=0. If a float is provided, it will be used as the HI density at all redshifts. If an array is provided, it will be used as the HI density at each frequency channel.
fiducial_cosmology (dict, default Planck18) – The fiducial cosmology parameters. Fiducial cosmology should be used to perform the power spectrum estimation, such as transforming sky coordinates to comoving coordinates.
true_cosmology (dict, default Planck18) – The true cosmology parameters. True cosmology should be varied during parameter inference.
ps_type (str, default "linear") – The type of the matter power spectrum.
kmin (float, default 1e-3) – The minimum k in Mpc^-1 for calculating matter power. k below kmin will be extrapolated.
kmax (float, default 3.0) – The maximum k in Mpc^-1 for calculating matter power. k above kmax will be extrapolated.
num_kpoints (int, default 200) – The number of k points to compute the interpolation of the matter power spectrum.
cold (bool, default True) – Whether to use the cold matter power spectrum. If True, the matter power spectrum refers to CDM+Baryon. If False, the matter power spectrum refers to all matter including massive neutrinos.
**params (dict) – Additional parameters to be passed to the base class meer21cm.dataanalysis.Specification

Attributes:

As: The amplitude of the primordial power spectrum for the true (fitted) cosmology.
W_HI: A binary window for whether a pixel has been sampled
alpha_AP: The anisotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel} / \alpha_{\perp}\).
alpha_iso: The isotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel}^{1/3} \alpha_{\perp}^{2/3}\).
alpha_parallel: The line of sight Alcock–Paczynski effect parameter.
alpha_perp
average_hi_temp: The average HI brightness temperature in Kelvin at central redshift self.z.
backend: Which backend to use for computing the matter power.
batch_number: Number of sequential batches used by gridding routines.
beam_image: Returns the beam image projected onto the sky map for the input beam model.
beam_model: The name of the beam function.
beam_type: The beam type that can be either be isotropic or anisotropic.
beam_unit: The unit of input beam size parameter sigma
beam_window_ch: Per-channel spherical-harmonic beam window \(B_\ell\) for HEALPix hp.smoothing workflows.
ch_id_gal: The channel id (0-indexed) of each galaxy in the catalogue for cross-correlation.
cold: If True (recommended), the matter power spectrum is the CDM ps.
cosmo: A shortcut to the astropy.cosmology.Cosmology object for the true cosmology.
cospar_fiducial
cospar_true
counts: The number of hits per pixel for the map data
data: The map data
dec_gal: The declination of each galaxy in the catalogue for cross-correlation.
dec_map: The declination of each pixel in the map.
dvdf: velocity resolution per unit frequency on average, in km/s/Hz
dvdf_ch: velocity resolution per unit frequency in each channel, in km/s/Hz
expfactor: The expansion factor
f_growth_fiducial: The growth factor at self.expfactor for the fiducial cosmology.
f_growth_true: The growth factor at self.expfactor for the true (fitted) cosmology.
fiducial_cosmology
freq_gal: The 21cm line frequency for each galaxy in Hz.
freq_resol: frequency resolution in Hz
h: The Hubble constant for the true (fitted) cosmology.
hp_nside: HEALPix \(N_{side}\) for HEALPix-backed specifications.
kmax: The maximum k in Mpc^-1 for calculating matter power.
kmin: The minimum k in Mpc^-1 for calculating matter power.
los_resol_in_mpc: effective frequency resolution in Mpc for the fiducial cosmology.
map_has_sampling: A binary window for whether a pixel has been sampled
map_unit_type: The type of the map unit.
matter_power_spectrum_fnc: Interpolation function for the real-space isotropic matter power spectrum.
maximum_sampling_channel: Returns the index of the frequency channel with the maximum sampling on the sky map.
neutrino_mass: The neutrino mass for the true (fitted) cosmology.
ns: The primordial power spectrum spectral index for the true (fitted) cosmology.
nu: The input frequencies of the survey
num_kpoints: The number of k points for calculating matter power.
num_pix_x: The number of pixels along the first map axis.
num_pix_y: The number of pixels along the second map axis.
omega_baryon: The baryon matter density parameter for the true (fitted) cosmology.
omega_cold: The cold matter density parameter for the true (fitted) cosmology.
omega_hi: The HI density as a function of redshift, over the critical density of the Universe at z=0, defined at each frequency channel so that self.omega_hi corresponds to self.z_ch.
omega_hi_z_func: Interpolate the input self.omega_hi at each frequency channel self.z_ch to a function of redshift.
omega_hi_z_mean: The mean HI density at the central redshift self.z, over the critical density of the Universe at z=0.
pix_resol: angular resolution of the map pixel in deg
pix_resol_in_mpc: angular resolution of the map pixel in Mpc for the fiducial cosmology.
pixel_area: angular area of the map pixel in deg^2
pixel_id: Sparse HEALPix pixel indices (RING, nest=False).
precision: Floating precision flag.
ps_type: linear or nonlinear for the matter power.
ra_gal: The right ascension of each galaxy in the catalogue for cross-correlation.
ra_map: The right ascension of each pixel in the map.
real_dtype: Active real floating dtype controlled by precision.
sigma_beam_ch: The input beam size parameter sigma for each channel.
sigma_beam_ch_in_mpc: The input beam size parameter in Mpc in each channel.
sigma_beam_in_mpc: The channel averaged beam size in Mpc
sound_horizon_drag_fiducial: The sound horizon at drag epoch for the fiducial cosmology.
sound_horizon_drag_true: The sound horizon at drag epoch for the true (fitted) cosmology.
survey_volume: Total survey volume in Mpc^3.
true_cosmology
vel_resol: velocity resolution on average in km/s
vel_resol_ch: velocity resolution of each channel in km/s
w0: The dark energy equation of state parameter for the true (fitted) cosmology.
w_HI: The weights per map pixel.
wa: The dark energy equation of state parameter for the true (fitted) cosmology.
weights_map_pixel: The weights per map pixel.
wproj: The WCS projection object for the map geometry.
z: The effective centre redshift of the frequency range
z_as_func_of_comov_dist: Returns a function that returns the redshift for input comoving distance.
z_ch: The redshift of each frequency channel
z_gal: The redshifts of each galaxy in the catalogue for cross-correlation.

Methods

`check_is_map_noiselike_using_pca`(A_mat[, ...])	Use the source mixing matrix from eigendecomposition of the covariance matrix, project out the map data with more and more modes, and check if the variance of the residual map behaves like white noise.
`clean_cache`(attr)	Set the attributes to None.
`convolve_data`([kernel, data, weights, ...])	Convolve map data.
`create_white_noise_map`(sigma_N[, counts, ...])	Create a white noise map with the given standard deviation.
`deltav_to_deltar`(delta_v)	Convert a velocity interval delta_v to a comoving distance interval delta_r.
`deltaz_to_deltar`(delta_z)	Convert a redshift interval delta_z to a comoving distance interval delta_r.
`generate_full_healpix_map`([data, fill_value])	Generate a full healpix map from the data.
`get_alpha_parallel`()	Calculate the line of sight Alcock–Paczynski effect parameter.
`get_alpha_perp`()	Calculate the transverse Alcock–Paczynski effect parameter.
`get_beam_image`([wproj, num_pix_x, ...])	Calculate the beam image projected onto the sky map for the input beam model.
`get_beam_window_ch`([ch_sel, cache, lmax, ...])	Build per-channel harmonic beam windows for HEALPix smoothing.
`get_cospar`(cosmology)	Generate a `CosmologyParameters` object from the input cosmology.
`get_jackknife_patches`(ra_patch_num, ...[, ...])	Split the map into roughly equal patches.
`get_matter_power_spectrum`()	Calculate the matter power spectrum, interpolate it, and save it into the class attribute matter_power_spectrum_fnc.
`get_sound_horizon_drag`(Omh2, Obh2, Onuh2)	Calculate the sound horizon at drag epoch for a set of cosmological parameters.
`get_weights_none_to_one`(attr_name)	Get the weights, and if it is None, convert it to 1.0 of size of kmode.
`get_z_as_func_of_comov_dist`()	Calculate an array of comoving distances with redshifts, and construct a function that returns the redshift for input comoving distance.
`read_from_fits`()	Read in a FITS file for the map data and hit counts.
`read_from_pickle`()	Read in a pickle file for cross-correlation and save the data into the class attributes.
`read_gal_cat`([ra_col, dec_col, z_col, trim])	Read in a galaxy catalogue for cross-correlation and save the data into the class attributes.
`set_radecnu_bounds_from_map`()	Set `ra_range`, `dec_range`, `nu_min`, and `nu_max` from the loaded `_ra_map`, `_dec_map`, and `nu` (tight RA interval, declination min/max, frequency min/max).
`trim_gal_to_range`()	Trim the galaxy catalogue to the specified range.
`trim_map_to_range`()	Trim the map to the specified range.

property As: The amplitude of the primordial power spectrum for the true (fitted) cosmology. Note it does not update fiducial cosmology for power spectrum estimation, and only used for model fitting.

property alpha_AP: The anisotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel} / \alpha_{\perp}\).

property alpha_iso: The isotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel}^{1/3} \alpha_{\perp}^{2/3}\).

property alpha_parallel: The line of sight Alcock–Paczynski effect parameter.

property alpha_perp

property average_hi_temp: The average HI brightness temperature in Kelvin at central redshift self.z. Calculation is based on the true (fitted) cosmology.

property backend: Which backend to use for computing the matter power. Either camb or bacco.

property cold: If True (recommended), the matter power spectrum is the CDM ps. If False, it will include massive neutrino for bacco and total matter for camb.

property cosmo

A shortcut to the astropy.cosmology.Cosmology object for the true cosmology.

Should only be used when true and fiducial cosmology are the same, and returns a warning if not.

If you set cosmo to a new value, it will set true_cosmology and fiducial_cosmology to the same value. If true_cosmology and fiducial_cosmology are different, an error will be raised.

property cospar_fiducial

property cospar_true

deltav_to_deltar(delta_v)[source]

Convert a velocity interval delta_v to a comoving distance interval delta_r. The conversion is based on the true (fitted) cosmology.

Parameters:: delta_v (float.) – The velocity interval in km/s.
Returns:: delta_r – The comoving distance interval in Mpc.
Return type:: float.

deltaz_to_deltar(delta_z)[source]

Convert a redshift interval delta_z to a comoving distance interval delta_r. The conversion is based on the true (fitted) cosmology.

Note that, the usual redshift error defined in galaxy survey is usually delta_z / (1+z).

Parameters:: delta_z (float.) – The redshift interval.
Returns:: delta_r – The comoving distance interval in Mpc.
Return type:: float.

property expfactor: The expansion factor

property f_growth_fiducial: The growth factor at self.expfactor for the fiducial cosmology.

property f_growth_true: The growth factor at self.expfactor for the true (fitted) cosmology.

property fiducial_cosmology

get_alpha_parallel()[source]: Calculate the line of sight Alcock–Paczynski effect parameter.

get_alpha_perp()[source]: Calculate the transverse Alcock–Paczynski effect parameter.

get_cospar(cosmology: dict)[source]

Generate a CosmologyParameters object from the input cosmology.

Parameters:: cosmology (dict) – The cosmology parameters.
Returns:: cospar – The cosmology parameters object.
Return type:: CosmologyParameters

get_matter_power_spectrum()[source]

Calculate the matter power spectrum, interpolate it, and save it into the class attribute matter_power_spectrum_fnc.

Note that, Alcock–Paczynski effect is included here in the power spectrum, and therefore for all model power spectra that depends on the matter power spectrum.

get_sound_horizon_drag(Omh2: float, Obh2: float, Onuh2: float)[source]

Calculate the sound horizon at drag epoch for a set of cosmological parameters.

Uses the fitting formula from [1411.1074](https://arxiv.org/abs/1411.1074)

get_z_as_func_of_comov_dist()[source]: Calculate an array of comoving distances with redshifts, and construct a function that returns the redshift for input comoving distance. The function is saved into the class attribute z_as_func_of_comov_dist. The comoving distance is calculated for the fiducial cosmology.

property h: The Hubble constant for the true (fitted) cosmology.

property kmax: The maximum k in Mpc^-1 for calculating matter power. k above kmax will be extrapolated.

property kmin: The minimum k in Mpc^-1 for calculating matter power. k below kmin will be extrapolated.

property los_resol_in_mpc: effective frequency resolution in Mpc for the fiducial cosmology.

property matter_power_spectrum_fnc: Interpolation function for the real-space isotropic matter power spectrum. The matter power spectrum is calculated for the true (fitted) cosmology.

property neutrino_mass: The neutrino mass for the true (fitted) cosmology.

property ns: The primordial power spectrum spectral index for the true (fitted) cosmology.

property num_kpoints: The number of k points for calculating matter power.

property omega_baryon: The baryon matter density parameter for the true (fitted) cosmology.

property omega_cold: The cold matter density parameter for the true (fitted) cosmology.

property omega_hi: The HI density as a function of redshift, over the critical density of the Universe at z=0, defined at each frequency channel so that self.omega_hi corresponds to self.z_ch. Interpolation will be used to get the HI density at other redshifts. If a float is provided, it will be used as the HI density at all redshifts.

property omega_hi_z_func: Interpolate the input self.omega_hi at each frequency channel self.z_ch to a function of redshift.

property omega_hi_z_mean: The mean HI density at the central redshift self.z, over the critical density of the Universe at z=0. Interpolated from the input self.omega_hi at each frequency channel self.z_ch.

property pix_resol_in_mpc: angular resolution of the map pixel in Mpc for the fiducial cosmology.

property ps_type: linear or nonlinear for the matter power.

property sigma_beam_ch_in_mpc: The input beam size parameter in Mpc in each channel. The comoving distance is calculated for the true (fitted) cosmology.

property sound_horizon_drag_fiducial: The sound horizon at drag epoch for the fiducial cosmology.

property sound_horizon_drag_true: The sound horizon at drag epoch for the true (fitted) cosmology.

property survey_volume

Total survey volume in Mpc^3. The volume is calculated for the fiducial cosmology.

Note that, the sampling along the sky map is assumed to be the same for all frequency channels, and the code by default uses the maximum sampling channel to calculate the area. This is desired, as the survey lightcone can contain holes inside, which is considered part of the volume.

Parameters:: i (int, default None) – The index of the frequency channel to calculate the survey volume. Default is None, which uses the maximum sampling channel.

property true_cosmology

property w0: The dark energy equation of state parameter for the true (fitted) cosmology.

property wa: The dark energy equation of state parameter for the true (fitted) cosmology.

property z_as_func_of_comov_dist: Returns a function that returns the redshift for input comoving distance. The comoving distance is calculated for the fiducial cosmology.

class meer21cm.cosmology.CosmologyParameters(ps_type='linear', kmin=0.001, kmax=3.0, num_kpoints=200, expfactor=1.0, cold=True, omega_de=None, tau=0.0561, Neff=3.046, omega_cold=0.30966, As=np.float64(2.105209331337507e-09), omega_baryon=0.04897, ns=0.9665, h=np.float64(0.6766), neutrino_mass=np.float64(0.06), w0=-1.0, wa=0.0)[source]

Bases: object

The class for storing cosmological parameters, and settings for computing matter power spectrum. The naming of the input arguments for cosmological parameters follow baccoemu . It either uses camb or baccoemu to compute the matter power spectrum.

Note that everything is not in h unit unless explicitly specified in name (of course except sigma_8 which follows the definition of 8 Mpc/h).

Further note that, baccoemu is trained on CLASS . Therefore, in the usual range of parameters in the LCDM, you should see the <1% difference between these two backends as differences between the Boltzmann solver codes (although this is not well tested on our end). Use it with precaution if you want to do precision cosmology type of forecasts and sims.

Parameters:

ps_type (str, default "linear") – The type of the matter power spectrum.
kmin (float, default 1e-3) – The minimum k in Mpc^-1 for calculating matter power. k below kmin will be extrapolated.
kmax (float, default 3.0) – The maximum k in Mpc^-1 for calculating matter power. k above kmax will be extrapolated.
omega_cold (float, default astropy.cosmology.Planck18.Om0) – The density fraction of CDM+Baryon at z=0.
As (float, default astropy.cosmology.Planck18.As) – The amplitude of the initial power spectrum.
omega_baryon (float, default astropy.cosmology.Planck18.Ob0) – The density fraction of baryons at z=0.
ns (float, default astropy.cosmology.Planck18.meta["n"]) – The spectral index of the initial power spectrum.
h (float, default astropy.cosmology.Planck18.h) – The Hubble parameter over 100km/s/Mpc.
neutrino_mass (float, default astropy.cosmology.Planck18.m_nu.sum().value) – The sum of the neutrino mass in eV.
w0 (float, default -1.0) – The dark energy equation of state at a=1 (z=0).
wa (float, default 0.0) – The redshift-dependent part of the dark energy equation of state. \(w(a) = w_0 + w_a (1 - a)\).
expfactor (float, default 1.0) – The expansion factor which is calculated as \(a = 1 / (1 + z)\).
cold (bool, default True) – Whether to use the cold matter power spectrum. If True, the matter power spectrum refers to CDM+Baryon. If False, the matter power spectrum refers to all matter including massive neutrinos.
num_kpoints (int, default 200) – The number of k points to compute the interpolation of the matter power spectrum.
omega_de (float, default None) – The density fraction of dark energy at z=0. If None, it will be calculated using camb from the rest of the input parameters.
tau (float, default 0.0561) – The optical depth of the reionization. Note that it does not affect the matter and tracer power spectrum calculation.
Neff (float, default 3.046) – The effective number of neutrino species. Note that it does not affect the matter and tracer power spectrum calculation.

Attributes:

omega_de: The dark energy density fraction at z=0.

Methods

`get_bacco_pars`()	Generate a dictionary that can be used as input for the `bacco` emulator.
`get_camb_pars`()	Generate a `camb.model.CAMBparams` set by the input cosmology.
`get_derived_Ode`()	Use camb to calculate the Ode0 given input parameters.
`get_matter_power_spectrum_bacco`()	Emulate the CDM power spectrum using bacco.
`get_matter_power_spectrum_camb`()	Compute the CDM power spectrum using camb.
`set_astropy_cosmo`([name])	Generate a `astropy.cosmology.w0waCDM` set by the input cosmology.

get_bacco_pars()[source]: Generate a dictionary that can be used as input for the bacco emulator. Currently only support non-baryonic matter power.

get_camb_pars()[source]: Generate a camb.model.CAMBparams set by the input cosmology.

get_derived_Ode()[source]: Use camb to calculate the Ode0 given input parameters.

get_matter_power_spectrum_bacco()[source]: Emulate the CDM power spectrum using bacco.

get_matter_power_spectrum_camb()[source]: Compute the CDM power spectrum using camb.

property omega_de: The dark energy density fraction at z=0.

set_astropy_cosmo(name='new')[source]: Generate a astropy.cosmology.w0waCDM set by the input cosmology. Note that the dark energy density Ode0 is a derived quantity which is calculated using camb if not put in.

meer21cm.cosmology.get_cosmo_dict(cosmo: str)[source]: Get the cosmology dictionary from the input cosmology.

meer21cm.cosmology.get_ns_from_astropy(x)

meer21cm.dataanalysis module

This module contains the base class for reading and visualizing the map data cube.

Note that, the defined class, Specification, is the base class for reading and visualizing the map data cube. It is typically used as a base class for other classes that inherit from it, and not used directly.

class meer21cm.dataanalysis.Specification(nu=None, wproj=None, num_pix_x=None, num_pix_y=None, map_has_sampling=None, sigma_beam_ch=None, beam_unit=Unit('deg'), map_unit=Unit('K'), map_file=None, counts_file=None, pickle_file=None, los_axis=-1, nu_min=None, nu_max=None, filter_map_los=True, gal_file=None, weighting='counts', ra_range=(0, 360), dec_range=(-90, 90), beam_model='gaussian', data=None, weights_map_pixel=None, counts=None, survey='', band='', z_interp_max=6.0, soft_filter_los=True, filter_los_threshold=None, data_column='map', counts_column='hit', freq_column='freq', wcs_column='wcs', auto_set_radecnu_bounds=True, precision=True, batch_number=1, skymap=None, hp_nside=None, healpix_pixel_id=None, b_ell_l_max=8192, **kwparams)[source]

Bases: object

Base class for reading and visualizing the map data cube.

Parameters:

nu (np.ndarray, default None) – The frequencies of the survey in Hz.
wproj (astropy.wcs.WCS, default None) – The WCS object for the map.
num_pix_x (int, default None) – The number of pixels in the first axis of the map data (WCS only).
num_pix_y (int, default None) – The number of pixels in the second axis of the map data (WCS only).
hp_nside (int, default None) – HEALPix \(N_{side}\). Implies sparse (n_pix, n_chan) layout via HealpixSkyMap. Incompatible with predefined survey/band maps. Mutually exclusive with passing skymap.
healpix_pixel_id (array-like, default None) – Optional explicit sparse pixel indices at hp_nside (RING, nest=False). If omitted, pixels are derived from ra_range and dec_range.
map_has_sampling (np.ndarray, default None) – A binary window for whether a pixel has been sampled.
sigma_beam_ch (np.ndarray, default None) – The beam size parameter for each frequency channel.
beam_unit (astropy.units.Unit, default astropy.units.deg) – The unit of the beam size parameter.
map_unit (astropy.units.Unit, default astropy.units.K) – The unit of the map data.
map_file (str, default None) – The file path of the map data. Supports automatic reading of the MeerKLASS L-band data. For UHF data use pickle_file for the file path of the pickle file.
counts_file (str, default None) – The file path of the hit counts data. Supports automatic reading of the MeerKLASS L-band data. For UHF data use pickle_file for the file path of the pickle file.
pickle_file (str, default None) – The file path of the pickle file. Supports automatic reading of the MeerKLASS UHF data.
los_axis (int, default -1) – The axis of the map data that corresponds to the line of sight. Warning: Tranposing the data to align the los axis is not properly taken care of in the code. If your map los axis is not the last axis, it is recommended to manually transpose the data so that the los axis is the last axis.
nu_min (float, default None,) – The minimum frequency of the survey in Hz. Data below this frequency will be clipped.
nu_max (float, default None,) – The maximum frequency of the survey in Hz. Data above this frequency will be clipped.
filter_map_los (bool, default True) – Whether to filter the map data along the line of sight. See meer21cm.io.filter_incomplete_los()
gal_file (str, default None,) – The file path of the galaxy catalogue.
weighting (str, default "counts") – The weighting scheme for the map data.
ra_range (tuple, default (0, 360)) – The range of the right ascension of the map data in degrees. Data outside this range will be masked.
dec_range (tuple, default (-90, 90)) – The range of the declination of the map data in degrees. Data outside this range will be masked.
beam_model (str, default "gaussian") – The shape of the beam.
data (np.ndarray, default None) – The map data.
weights_map_pixel (np.ndarray, default None) – The weights per map pixel.
counts (np.ndarray, default None) – The number of hits per pixel for the map data.
survey (str, default "") – The survey name.
band (str, default "") – The band of the survey.
z_interp_max (float, default 6.0) – The maximum redshift to interpolate the redshift as a function of comoving distance. See meer21cm.dataanalysis.Specification.get_z_as_func_of_comov_dist().
soft_filter_los (bool, default True) – If filter_map_los is True, whether to use a soft criterion. If False, any line of sight that is not 100% sampled will be removed. If True, the maximum sampling fraction of the map cube is calculated and used as the criterion. See meer21cm.io.filter_incomplete_los().
filter_los_threshold (float, default None) – If given, instead of filtering out incomplete los by checking the maximum sampling fraction along the los, a fixed threshold is used to filter out incomplete los. See meer21cm.io.filter_incomplete_los().
data_column (str, default "map") – The column name of the map data.
counts_column (str, default "hit") – The column name of the number of sampling for each pixel.
freq_column (str, default "freq") – The column name of the frequencies of each channel in the data.
wcs_column (str, default "wcs") – The column name of the astropy.wcs.WCS object for the map.
auto_set_radecnu_bounds (bool, default True) – If True, read_from_fits() and read_from_pickle() call set_radecnu_bounds_from_map() after loading so ra_range, dec_range, nu_min, and nu_max match the loaded grid and channels.
precision (bool, default True) – Floating precision selector for core numeric arrays. If True, use double precision (np.float64); if False, use single precision (np.float32).
batch_number (int, default 1) – Number of sequential batches used by various routines. A value of 1 means no batching.
skymap (SkyMap, default None) – Injected angular geometry (WcsSkyMap or HealpixSkyMap). Mutually exclusive with hp_nside.
b_ell_l_max (int, default 8192) – The maximum multipole for the beam window function for healpix map smoothing.

Attributes:

W_HI: A binary window for whether a pixel has been sampled
batch_number: Number of sequential batches used by gridding routines.
beam_image: Returns the beam image projected onto the sky map for the input beam model.
beam_model: The name of the beam function.
beam_type: The beam type that can be either be isotropic or anisotropic.
beam_unit: The unit of input beam size parameter sigma
beam_window_ch: Per-channel spherical-harmonic beam window \(B_\ell\) for HEALPix hp.smoothing workflows.
ch_id_gal: The channel id (0-indexed) of each galaxy in the catalogue for cross-correlation.
counts: The number of hits per pixel for the map data
data: The map data
dec_gal: The declination of each galaxy in the catalogue for cross-correlation.
dec_map: The declination of each pixel in the map.
dvdf: velocity resolution per unit frequency on average, in km/s/Hz
dvdf_ch: velocity resolution per unit frequency in each channel, in km/s/Hz
freq_gal: The 21cm line frequency for each galaxy in Hz.
freq_resol: frequency resolution in Hz
hp_nside: HEALPix \(N_{side}\) for HEALPix-backed specifications.
map_has_sampling: A binary window for whether a pixel has been sampled
map_unit_type: The type of the map unit.
maximum_sampling_channel: Returns the index of the frequency channel with the maximum sampling on the sky map.
nu: The input frequencies of the survey
num_pix_x: The number of pixels along the first map axis.
num_pix_y: The number of pixels along the second map axis.
pix_resol: angular resolution of the map pixel in deg
pixel_area: angular area of the map pixel in deg^2
pixel_id: Sparse HEALPix pixel indices (RING, nest=False).
precision: Floating precision flag.
ra_gal: The right ascension of each galaxy in the catalogue for cross-correlation.
ra_map: The right ascension of each pixel in the map.
real_dtype: Active real floating dtype controlled by precision.
sigma_beam_ch: The input beam size parameter sigma for each channel.
sigma_beam_in_mpc: The channel averaged beam size in Mpc
vel_resol: velocity resolution on average in km/s
vel_resol_ch: velocity resolution of each channel in km/s
w_HI: The weights per map pixel.
weights_map_pixel: The weights per map pixel.
wproj: The WCS projection object for the map geometry.
z: The effective centre redshift of the frequency range
z_ch: The redshift of each frequency channel
z_gal: The redshifts of each galaxy in the catalogue for cross-correlation.

Methods

`check_is_map_noiselike_using_pca`(A_mat[, ...])	Use the source mixing matrix from eigendecomposition of the covariance matrix, project out the map data with more and more modes, and check if the variance of the residual map behaves like white noise.
`clean_cache`(attr)	Set the attributes to None.
`convolve_data`([kernel, data, weights, ...])	Convolve map data.
`create_white_noise_map`(sigma_N[, counts, ...])	Create a white noise map with the given standard deviation.
`generate_full_healpix_map`([data, fill_value])	Generate a full healpix map from the data.
`get_beam_image`([wproj, num_pix_x, ...])	Calculate the beam image projected onto the sky map for the input beam model.
`get_beam_window_ch`([ch_sel, cache, lmax, ...])	Build per-channel harmonic beam windows for HEALPix smoothing.
`get_jackknife_patches`(ra_patch_num, ...[, ...])	Split the map into roughly equal patches.
`get_weights_none_to_one`(attr_name)	Get the weights, and if it is None, convert it to 1.0 of size of kmode.
`read_from_fits`()	Read in a FITS file for the map data and hit counts.
`read_from_pickle`()	Read in a pickle file for cross-correlation and save the data into the class attributes.
`read_gal_cat`([ra_col, dec_col, z_col, trim])	Read in a galaxy catalogue for cross-correlation and save the data into the class attributes.
`set_radecnu_bounds_from_map`()	Set `ra_range`, `dec_range`, `nu_min`, and `nu_max` from the loaded `_ra_map`, `_dec_map`, and `nu` (tight RA interval, declination min/max, frequency min/max).
`trim_gal_to_range`()	Trim the galaxy catalogue to the specified range.
`trim_map_to_range`()	Trim the map to the specified range.

property W_HI: A binary window for whether a pixel has been sampled

property batch_number: Number of sequential batches used by gridding routines.

property beam_image: Returns the beam image projected onto the sky map for the input beam model.

property beam_model: The name of the beam function.

property beam_type: The beam type that can be either be isotropic or anisotropic.

property beam_unit: The unit of input beam size parameter sigma

property beam_window_ch: Per-channel spherical-harmonic beam window \(B_\ell\) for HEALPix hp.smoothing workflows. Populated lazily via get_beam_window_ch().

property ch_id_gal: The channel id (0-indexed) of each galaxy in the catalogue for cross-correlation. Galaxies out of the frequency range will be given len(self.nu) as indices.

check_is_map_noiselike_using_pca(A_mat, data=None, sigma_N=1.0)[source]

Use the source mixing matrix from eigendecomposition of the covariance matrix, project out the map data with more and more modes, and check if the variance of the residual map behaves like white noise.

You can use meer21cm.util.pca_clean() to retrieve the source mixing matrix: .. code-block:: python

N_fg = 15 # check 15 modes removed res_map, A_mat = pca_clean(ps.data, N_fg, weights=ps.W_HI, return_A=True) res, noise = ps.check_is_map_noiselike_using_pca(A_mat) plt.plot(res / noise)

If the residual map is noise-like, the plot should decrease and eventually reach a plateau.

If you know the expected std of the map (per hit), you can pass it to sigma_N to scale the noise variance, and the plateau should be close to 1.

Note that, the input data should be the mean-centered data. You can use meer21cm.util.mean_center_signal() to mean-center the data if needed.

Parameters:

A_mat (array.) – The source mixing matrix.
data (array, default None.) – The data to be projected out. If None, the class attribute self.data will be used.
sigma_N (float.)

clean_cache(attr)[source]: Set the attributes to None. This is used to clear the cache of the attributes.

convolve_data(kernel=None, data=None, weights=None, assign_to_self=True)[source]

Convolve map data.

WCS maps use raster-space weighted_convolution() with a beam image kernel (e.g. from beam_image()).

HEALPix maps ignore raster kernel and use harmonic weighted_smoothing_healpix() with per-channel beam windows from get_beam_window_ch() — pass kernel=None.

Parameters:

kernel (np.ndarray or None) – Raster beam cube (..., nx, ny, n_ch) aligned with LOS on wcs. Ignored when kernel is None on HEALPix (then harmonic smoothing is applied). On WCS maps, passing kernel=None raises ValueError (provide e.g. self.beam_image).
data (np.ndarray, default None) – Map to convolve; default self.data.
weights (np.ndarray, default None) – Per-pixel weights (e.g. self.w_HI); required semantics match data.
assign_to_self (bool, default True) – Assign results to self.data and self.w_HI.

Returns:

data (np.ndarray)
weights (np.ndarray)

property counts: The number of hits per pixel for the map data

create_white_noise_map(sigma_N, counts=None, seed=None, inf_to_zero=True)[source]

Create a white noise map with the given standard deviation. The sigma in each pixel is then scaled by the counts 1/sqrt(counts).

If you want to generate multiple random catalogues, you need to set a different seed manually for each catalogue.

If you want to use different noise level per pixel, you can either pass a 3D array of sigma_N, or a single number and a 3D array of counts. You can usually pass self.counts as the counts array, but do check the counts are set up correctly by plot_map(self.counts, self.wproj).

Finally, note that the noise map is not masked by the survey selection function. You can mask the noise map manually by noise_map *= self.W_HI.

Parameters:

sigma_N (float or array.) – The standard deviation of the white noise.
counts (array, default None.) – The counts in each pixel. If None, the counts will be one across the cube.
seed (int, default None.) – If none, the seed is pulled from OS as described by the numpy documentation.
inf_to_zero (bool, default True.) – If True, the inf values in the noise map will be set to zero.

Returns:

noise_map – The white noise map.

Return type:

array.

property data: The map data

property dec_gal: The declination of each galaxy in the catalogue for cross-correlation.

property dec_map: The declination of each pixel in the map.

property dvdf: velocity resolution per unit frequency on average, in km/s/Hz

property dvdf_ch: velocity resolution per unit frequency in each channel, in km/s/Hz

property freq_gal: The 21cm line frequency for each galaxy in Hz.

property freq_resol: frequency resolution in Hz

generate_full_healpix_map(data=None, fill_value=nan)[source]

Generate a full healpix map from the data. Pixels not included in the data will be filled with the fill_value (default is np.nan).

Parameters:

data (array, default None.) – The data to be included in the full healpix map. If None, the class attribute self.data will be used.
fill_value (float, default np.nan.) – The value to fill the pixels not included in the data.

Returns:

full_healpix_map – The full healpix map.

Return type:

array.

get_beam_image(wproj=None, num_pix_x=None, num_pix_y=None, cache=True, ch_sel=None)[source]

Calculate the beam image projected onto the sky map for the input beam model.

Parameters:

wproj (astropy.wcs.WCS, default None) – The WCS object for the map. Default uses self.wproj.
num_pix_x (int, default None) – The number of pixels in the first axis of the map data. Default uses self.num_pix_x.
num_pix_y (int, default None) – The number of pixels in the second axis of the map data. Default uses self.num_pix_y.
cache (bool, default True) – Whether to cache the beam image. Default is True. If True, the beam image will be cached and returned directly if it is already computed. If False, the beam image will be computed and returned. The cache is saved into the class attribute beam_image.
ch_sel (array-like, default None) – Optional channel selection. If provided, returns beam image only for selected channels. Caching is only applied when ch_sel is None.

get_beam_window_ch(ch_sel=None, cache=True, lmax=None, hp_nside=None)[source]

Build per-channel harmonic beam windows for HEALPix smoothing.

Parameters:

ch_sel (array-like of int, optional) – Channel indices. If omitted, uses all len(self.nu) channels.
cache (bool, default True) – Cache the full-channel result in _beam_window_ch when ch_sel spans all channels.
lmax (int, optional) – Maximum multipole (inclusive). Defaults to min(3 * hp_nside - 1, 8192).
hp_nside (int, optional) – HEALPix \(N_{\rm side}\) used to set the default lmax. Defaults to self.hp_nside.

Returns:

Beam window rows \(B_\ell\) per selected channel.

Return type:

ndarray of shape (len(ch_sel), lmax + 1)

get_jackknife_patches(ra_patch_num, dec_patch_num, nu_patch_num, ra_range=None, dec_range=None, nu_range=None)[source]

Split the map into roughly equal patches. Each patch can then be masked, which can be used for jackknife resampling for covariance estimation. Note that the masks=True is where the pixels should be masked. So for example, if you want to exclude a patch, the correct survey window is then self.W_HI * (1-mask_arr[i]) and the weights self.w_HI * (1-mask_arr[i]).

If you want to examine the patch splits, you can visualise the mask array by using meer21cm.plot.visualise_patch_split().

Parameters:

ra_patch_num (int) – The number of patche grids in the right ascension direction.
dec_patch_num (int) – The number of patche grids in the declination direction.
nu_patch_num (int) – The number of patche grids in the frequency direction.
ra_range (tuple, default None) – The range of the right ascension of the map data in degrees. Default uses self.ra_range.
dec_range (tuple, default None) – The range of the declination of the map data in degrees. Default uses self.dec_range.
nu_range (tuple, default None) – The range of the frequency of the map data in Hz. Default uses [self.nu.min() - self.freq_resol/2, self.nu.max() + self.freq_resol/2].

get_weights_none_to_one(attr_name)[source]: Get the weights, and if it is None, convert it to 1.0 of size of kmode. Only used for power spectrum calculation. Defined here for inheritance.

property hp_nside: HEALPix \(N_{side}\) for HEALPix-backed specifications.

property map_has_sampling: A binary window for whether a pixel has been sampled

property map_unit_type: The type of the map unit. If the map unit is temperature, return “T”. If the map unit is flux density, return “F”. If the map unit is not temperature or flux density, raise an error.

property maximum_sampling_channel: Returns the index of the frequency channel with the maximum sampling on the sky map.

property nu: The input frequencies of the survey

property num_pix_x: The number of pixels along the first map axis.

property num_pix_y: The number of pixels along the second map axis.

property pix_resol: angular resolution of the map pixel in deg

property pixel_area: angular area of the map pixel in deg^2

property pixel_id: Sparse HEALPix pixel indices (RING, nest=False).

property precision: Floating precision flag. True for float64, False for float32.

property ra_gal: The right ascension of each galaxy in the catalogue for cross-correlation.

property ra_map: The right ascension of each pixel in the map.

read_from_fits()[source]

Read in a FITS file for the map data and hit counts. The FITS file need to follow the format of the MeerKLASS L-band data. See meer21cm.io.read_map() for more details.

After reading the data, the map data and hit counts are filtered along the frequency direction (see meer21cm.io.filter_incomplete_los()), and trimmed to the specified range (see meer21cm.dataanalysis.Specification.trim_map_to_range()). The weights per pixel are set to the hit counts if self.weighting is “counts”, or set to 1 if self.weighting is “uniform”.

read_from_pickle()[source]: Read in a pickle file for cross-correlation and save the data into the class attributes. See meer21cm.io.read_pickle() for more details.

read_gal_cat(ra_col='RA', dec_col='DEC', z_col='Z', trim=True)[source]

Read in a galaxy catalogue for cross-correlation and save the data into the class attributes. The data is read from the gal_file, which has to be a FITS file.

Parameters:

ra_col (str, default "RA") – The column name of the right ascension in the galaxy catalogue.
dec_col (str, default "DEC") – The column name of the declination in the galaxy catalogue.
z_col (str, default "Z") – The column name of the redshift in the galaxy catalogue.
trim (bool, default True) – Whether to trim the galaxy catalogue to the ra,dec,z range of the map. See meer21cm.dataanalysis.Specification.trim_gal_to_range().

property real_dtype: Active real floating dtype controlled by precision.

set_radecnu_bounds_from_map()[source]: Set ra_range, dec_range, nu_min, and nu_max from the loaded _ra_map, _dec_map, and nu (tight RA interval, declination min/max, frequency min/max). Only consider unmaksed pixels.

property sigma_beam_ch: The input beam size parameter sigma for each channel. If one number is provided, it will be used for all channels.

property sigma_beam_in_mpc: The channel averaged beam size in Mpc

trim_gal_to_range()[source]

Trim the galaxy catalogue to the specified range. The galaxy catalogue outside the ra-dec-z range will be removed.

Note that, a small buffer corresponding to half of the frequency channel bandwidth is added to the redshift range.

trim_map_to_range()[source]: Trim the map to the specified range. The map data and counts outside the range will be set to zero. The map_has_sampling and weights_map_pixel will be set to False outside the range.

property vel_resol: velocity resolution on average in km/s

property vel_resol_ch: velocity resolution of each channel in km/s

property w_HI: The weights per map pixel.

property weights_map_pixel: The weights per map pixel.

property wproj: The WCS projection object for the map geometry.

property z: The effective centre redshift of the frequency range

property z_ch: The redshift of each frequency channel

property z_gal: The redshifts of each galaxy in the catalogue for cross-correlation.

meer21cm.fg module

class meer21cm.fg.ForegroundSimulation(hp_nside=256, wproj=None, num_pix_x=None, num_pix_y=None, backend='haslam', pysm_preset_strings=['d1', 's1', 'f1', 'a1', 'c1'], sp_indx_for_haslam_backend=-2.0, coord_system='C')[source]

Bases: object

Foreground simulation class. All outputs are in units of K_RJ (Kelvin).

Parameters:

hp_nside (int) – HEALPix nside.
wproj (tuple) – WCS projection parameters.
num_pix_x (int) – Number of pixels in x direction.
num_pix_y (int) – Number of pixels in y direction.
backend (str, default 'haslam') – Backend to use for foreground simulation. Options: ‘gdsm’ (Global Sky Model), ‘pysm’ (PySM3), ‘haslam’ (Haslam 408 MHz map).
pysm_preset_strings (list, default ['d1', 's1', 'f1', 'a1', 'c1']) – List of PySM3 preset strings for included components. Default: [‘d1’, ‘s1’, ‘f1’, ‘a1’, ‘c1’].
sp_indx_for_haslam_backend (float, default -2.0) – Index of the spectral index for the Haslam 408 MHz map. Only used for ‘haslam’ backend.
coord_system (str, default 'C') – Coordinate system to use for the foreground simulation. Options: ‘G’ (Galactic), ‘C’ (Celestial), ‘E’ (Ecliptic).

Methods

`fg_wcs_cube`(freq)	Generate WCS cube of foregrounds.
`healpix_gen_gdsm`(freq)	Generate HEALPix map of foregrounds using Global Sky Model.
`healpix_gen_haslam`(freq)	Generate HEALPix map of foregrounds using Haslam 408 MHz map.
`healpix_gen_pysm`(freq)	Generate HEALPix map of foregrounds using PySM3.

fg_wcs_cube(freq)[source]

Generate WCS cube of foregrounds. The function will call the appropriate healpix_gen_* function based on the backend, and then project the map to the WCS coordinates if provided.

Parameters:: freq (float or array-like) – Frequency in Hz.
Returns:: out_map – WCS cube of foregrounds in units of K_RJ.
Return type:: numpy.ndarray

healpix_gen_gdsm(freq)[source]

Generate HEALPix map of foregrounds using Global Sky Model.

Parameters:: freq (float or array-like) – Frequency in Hz.
Returns:: cube – HEALPix map of foregrounds in units of K_RJ.
Return type:: numpy.ndarray

healpix_gen_haslam(freq)[source]

Generate HEALPix map of foregrounds using Haslam 408 MHz map.

Parameters:: freq (float or array-like) – Frequency in Hz.
Returns:: cube – HEALPix map of foregrounds in units of K_RJ.
Return type:: numpy.ndarray

healpix_gen_pysm(freq)[source]

Generate HEALPix map of foregrounds using PySM3.

Parameters:: freq (float or array-like) – Frequency in Hz.
Returns:: cube – HEALPix map of foregrounds in units of K_RJ.
Return type:: numpy.ndarray

meer21cm.grid module

meer21cm.grid.compensate_grid_window_effects(field_in_real_space, grid_scheme='nnb')[source]: Apply correction to cancel the windowing effects from discretization of fields into grids.

meer21cm.grid.find_rotation_matrix(vec)[source]

find the rotation matrix to rotate the input vector to (0,0,1).

Note that in 3D space, the rotation is not unique. For simplicity, this function first finds the rotational matrix so that the vector (x,y,z) is first rotated to \((\sqrt{x^2+y^2},0,z)\), and then find another matrix to rotate the vector to (0,0,1).

Parameters:: vec (numpy array.) – The input unit vector
Returns:: rot_mat – The rotational matrix so that rot_mat @ vec is np.array([0,0,1]).
Return type:: numpy array.

meer21cm.grid.fourier_window_for_assignment(num_mesh, window='nnb')[source]

Calculate the effective window function in Fourier space from mass assignment scheme that sample continueous fields to discrete grids.

The window function can be written as [1]

\[W(k_x,k_y,k_z) = \Bigg({\rm sinc}\bigg(\frac{k_x H_x}{2}\bigg) {\rm sinc}\bigg(\frac{k_y H_y}{2}\bigg) {\rm sinc}\bigg(\frac{k_z H_z}{2}\bigg)\Bigg)^p,\]

where \(k_{x,y,z}\) is the wavenumber of the grid in Fourier space and \(H_{x,y,z}\) is the length of the grid in real space. \(p\) is the power index related to the mass assignment scheme, and is equal to [1,2,3,4] for [nnb,cic,tsc,pcs]

Parameters:

num_mesh (list) – The number of grids on each side
window (str, default "nnb".) – The mass assignment scheme

Returns:

window_in_fourier – The window function in Fourier space

Return type:

numpy array

References

[1]

Sefusatti, E. et al., “Accurate Estimators of Correlation Functions in Fourier Space”, https://ui.adsabs.harvard.edu/abs/2016MNRAS.460.3624S.

meer21cm.grid.interlace_two_fields(real_field_1, real_field_2, shift)[source]

Interlacing two fields, where one is the shifted version of the other.

Parameters:

real_field_1 (array-like.) – The first field for interlacing
real_field_2 (array-like.) – The second field for interlacing, which should be a shifted version of the first.
shift (float.) – The shift of the field when performing Fourier transform, in the unit of cell size.

Returns:

interlaced_field – The interlaced field.

Return type:

array-like.

meer21cm.grid.minimum_enclosing_box_of_lightcone(ra_arr, dec_arr, freq, cosmo=FlatLambdaCDM(name='Planck18', H0=<Quantity 67.66 km / (Mpc s)>, Om0=0.30966, Tcmb0=<Quantity 2.7255 K>, Neff=3.046, m_nu=<Quantity [0., 0., 0.06] eV>, Ob0=0.04897), ang_unit='deg', tile=True, return_coord=False, buffkick=0.0, rot_mat=None)[source]

This functions finds a rotational axis to rotate the sky vectors of input coordinates so that the (crude) mean of the coordinates is at (0,0,1), and then finds the enclosing cuboid box for the coordinates. The box is not really minimum but should be quite optimal.

The function also returns a rotational matrix for rotating the coordinates in the cuboid back to the sky positions. For any point in the box pos = np.array([x,y,z]), you can find its RA and Dec by performing the rotation

vec = inv_rot @ pos
vec /= np.sqrt(np.sum(vec**2))
ra_pos, dec_pos = hp.vec2ang(vec,lonlat=True)

Parameters:

ra_arr (numpy array.) – The RA of the coordinates
dec_arr (numpy array.) – The Dec of the coordinates
freq (numpy array.) – The frequencies of the coordinates.
cosmo (astropy.cosmology.Cosmology object, default astropy.cosmology.Planck18.) – The input cosmology for converting frequencies to los length.
ang_unit (str or astropy.units.Unit) – The unit of the input angular coordinates.
tile (bool, default True.) – Whether to tile the input cooridnates so that the output is a meshgrid of input angular coordinates and frequencies.
return_coord (bool, default True.) – If True, also returns the corrosponding (x,y,z) coordinate of the input coordinates.
buffkick (float, default 0.0.) – The box is extended by buffkick on each end of each dimension.
rot_mat (numpy array, default None.) – If specified, override the rotation matrix calculated from the mean cooridnate.

Returns:

x_min (float.) – The origin of the box along x-axis.
y_min (float.) – The origin of the box along y-axis.
z_min (float.) – The origin of the box along z-axis.
x_len (float.) – The length of the box along x-axis.
y_len (float.) – The length of the box along y-axis.
z_len (float.) – The length of the box along z-axis.
inv_rot (numpy array.) – The rotational matrix to rotate the box back to the sky positions.
pos_arr (``numpy’’ array.) – Only returns if ``return_coord = True’’. The Cartesian coordinates of the input ra and dec.

meer21cm.grid.particle_to_mesh_distance(particle_pos, box_len, box_ndim)[source]

Calculate the distance between particles and the nearest mesh center. The distance is in the unit of the cell size. For particles outside the box, the nearest mesh center is the one on the boundary.

Parameters:

particle_pos (array.) – The coordinates of the particles.
box_len (array.) – The length of the box on each side
box_ndim (array.) – The number of grids on each side

Returns:

dist (array.) – The distance between the particles and the nearest mesh center.
indx_grid (array.) – The index of the nearest mesh center.

meer21cm.grid.project_function(s_arr, grid_scheme='nnb')[source]

Return the weighting function for the given mass assignment scheme and input 1D distance. The distance is in the unit of the cell size.

Parameters:

s_arr (float array.) – The distance in the unit of the cell size.
grid_scheme (str, default 'nnb') – The mass assignment scheme.

Returns:

weight_arr – The weighting function.

Return type:

float array.

meer21cm.grid.project_particle_to_regular_grid(particle_pos, box_len, box_ndim, grid_scheme='nnb', particle_mass=None, particle_weights=None, average=True, shift=0.0, compensate=False)[source]

Project particles into a regular grid with a certain mass assignment scheme.

Parameters:

particle_pos (array.) – The coordinates of the particles.
box_len (array.) – The length of the box on each side
box_ndim (array.) – The number of grids on each side
grid_scheme (str, default 'nnb') – The mass assignment scheme.
particle_mass (array, default None.) – The mass of each particle.
particle_weights (array, default None.) – The weights of each particle.
average (bool, default True.) – The grid values are weighted averages of the particles if True and weighted sums of the particles if False.
shift (float, default 0.0.) – Shift the position of the particles by the same amount in all directions, in the unit of cell size.

Returns:

mesh_mass (array.) – The mass of each grid.
mesh_weights (array.) – The weights of each grid.
mesh_counts (array.) – The effective number of particles in each grid.

meer21cm.grid.rotation_matrix_to_radec0(ra, dec)[source]: Find the rotation matrix to rotate the input point at (ra, dec) to (0, 0), by first rotating to (0, dec) and then to (0, 0).

meer21cm.grid.shot_noise_correction_from_gridding(box_ndim, grid_scheme)[source]

Calculate the multiplicative correction from gridding to the shot noise. Support ‘nnb’, ‘cic’ and ‘tsc’. The correction is taken from Jing (2005), astro-ph/0409240.

Parameters:

box_ndim (array.) – The number of grids on each side.
grid_scheme (str.) – The mass assignment scheme.

Returns:

shot_noise_correction – The multiplicative correction from gridding to the shot noise.

Return type:

array.

meer21cm.grid.sky_partition_for_radecrange(ra_range, dec_range, nside_out=128, nside_in=1024, dec_pad=0)[source]

Find a partition of the sky, so that each patch can be rotated to cover the specified RA and Dec range.

Parameters:

ra_range (array_like) – The range of RA to cover.
dec_range (array_like) – The range of Dec to cover.
nside_out (int, default 128) – The HEALPix NSIDE of the output map pixel id.
nside_in (int, default 1024) – The HEALPix NSIDE of the map pixel id for inner calculation.
dec_pad (int, default 0) – The number of extra rows to pad in Dec. Increasing this number will result in patches overlapping with each other.

Returns:

pix_id_for_patch_i (list) – The list of pixel id for each patch.
rot_mat_for_patch_i (list) – The list of rotation matrix for each patch, to rotate the patch back to cover the range.

meer21cm.inference module

This module contains the sampler classes for performing model fitting.

Similar to meer21cm.transfer, sampler class takes in the attributes of a meer21cm.power.PowerSpectrum instance, and perform the model fitting that includes the observational effects, window functions etc defined in the PS instance. The difference is that, the sampler class takes the dictionary of attributes instead of the PS itself, so you should always use meer21cm.inference.extract_model_fitting_inputs() to extract the attributes from the PS instance, and use the dictionary to initialize the sampler class. For example:

>>> from meer21cm.inference import extract_model_fitting_inputs, SamplerEmcee
>>> from meer21cm.power import PowerSpectrum
>>> # the ps instance is usually something you have already defined and used to read data
>>> # or it can be a MockSimulation instance you have used to generate the mock
>>> ps = PowerSpectrum(band="L", survey="meerklass_2021", sigma_v_1=100, tracer_bias_1=1.5, tracer_bias_2=2.0)
>>> ps_dict = extract_model_fitting_inputs(ps)
>>> sampler = SamplerEmcee(ps_dict, ...)

similarily, you can also use meer21cm.inference.SamplerNautilus to perform model fitting.

Caution:

data_vector, data_covariance and observables must be in the same order. For example, if observables is [“cross”, “hi-auto”], then the first half of data_vector should be the cross-power and second half should be the hi-auto power. Similarly, the [:half, :half] of data_covariance should be the cross-power covariance and the [half:, half:] should be the hi-auto power covariance (the rest of the matrix should be the cross-correlation between cross and auto powers).
params_name must be in the same order as the params_prior.

class meer21cm.inference.SamplerBase(ps_dict: dict, data_vector: ndarray, data_covariance: ndarray, params_name: list[str], params_prior: list[tuple[str, float, float]], observables: list[str] = ['cross'], save: bool = False, save_filename: str | None = None, save_model_blobs: bool = False, do_hartlap_correction: bool = False, num_mocks: int | None = None, do_percival_correction: bool = False)[source]

Bases: object

Base class for all samplers.

Parameters:

ps_dict (dict) – A dictionary of model fitting inputs.
data_vector (np.ndarray) – The data vector.
data_covariance (np.ndarray) – The data covariance matrix.
params_name (list[str]) – The names of the parameters. Must match the attribute names of the meer21cm.power.PowerSpectrum instance.
params_prior (list[tuple[str, float, float]]) – The prior distribution of the parameters. The first element of the tuple is the distribution type, which can be “uniform” or “gaussian”. The second and third elements are the parameters of the distribution. For uniform distribution, the second and third elements are the lower and upper bounds of the distribution. For gaussian distribution, the second and third elements are the mean and standard deviation of the distribution.
observables (list[str], default ["cross"]) – The observables to fit. Must be a list of “cross”, “hi-auto”, or “gg-auto”.
save (bool, default False) – Whether to save the fitting results while running the sampler.
save_filename (str | None, default None) – The filename to save the fitting results when save is True. If not provided, the results will not be saved.
save_model_blobs (bool, default False) – Whether to save the model blobs. The model blobs are the model vectors listed in the observables argument.
do_hartlap_correction (bool, default False) – Whether to apply the Hartlap correction to the data covariance matrix.
num_mocks (int, default None) – The number of mock data vectors to use for the Hartlap correction. If not provided, the Hartlap correction will not be applied.
do_percival_correction (bool, default False) – Whether to apply the Percival correction to the data covariance matrix.

Attributes:

data_covariance: The data covariance matrix.
do_hartlap_correction: Whether to apply the Hartlap correction to the data covariance matrix.
do_percival_correction: Whether to apply the Percival correction to the data covariance matrix.
hartlap_factor: The Hartlap factor.
inverse_covariance: The matrix inverse of the data covariance.
num_mocks: The number of mock data vectors to use for the Hartlap and Percival corrections.
percival_factor: The Percival factor.

Methods

`compute_log_likelihood`(params)	Calculate the log likelihood from the input parameter values.
`get_model_instance`(params)	Get the model instance from the input `ps_dict` and the updated parameters `params`.
`get_model_vector`(params)	Calculate the model vector from the input parameter values.
`validate_input`()	Validate the input parameters.

compute_log_likelihood(params: list[float] | ndarray) → float | tuple[float, ndarray][source]

Calculate the log likelihood from the input parameter values.

Parameters:: params (list[float] | np.ndarray) – The updated values of the parameters defined in the params_name argument.
Returns:: log_likelihood – The log likelihood.
Return type:: float

property data_covariance: ndarray: The data covariance matrix.

property do_hartlap_correction: bool: Whether to apply the Hartlap correction to the data covariance matrix.

property do_percival_correction: bool: Whether to apply the Percival correction to the data covariance matrix.

get_model_instance(params: list[float] | ndarray) → PowerSpectrum[source]

Get the model instance from the input ps_dict and the updated parameters params.

Parameters:: params (list[float] | np.ndarray) – The updated values of the parameters defined in the params_name argument.
Returns:: model_instance – The model PS instance to calculate the model vector.
Return type:: meer21cm.power.PowerSpectrum

get_model_vector(params: list[float] | ndarray) → ndarray[source]

Calculate the model vector from the input parameter values.

Parameters:: params (list[float] | np.ndarray) – The updated values of the parameters defined in the params_name argument.
Returns:: model_vector – The model vector. The shape is (len(observables), len(k1dbins) - 1).
Return type:: np.ndarray

property hartlap_factor: float

The Hartlap factor. The Hartlap correction is given by [1] :

\[\hat{\Sigma^{-1}} = \frac{N_m - N_d - 2}{N_m - 1} \Sigma^{-1}\]

where \(N_m\) is the number of mock realisations for calculating the covariance, \(N_d\) is the number of data points, and \(\Sigma\) is the mock-based covariance matrix.

References

[1]

Hartlap, J., et al. 2007, A&A, astro-ph/0608064

property inverse_covariance: ndarray: The matrix inverse of the data covariance.

property num_mocks: int | None: The number of mock data vectors to use for the Hartlap and Percival corrections.

property percival_factor: float

The Percival factor. The Percival correction is given by [1] :

\[\hat{\Sigma^{-1}} = \frac{1 + B(N_d - N_p)}{1 + A + B(N_p + 1)} \Sigma^{-1}\]

\[A = \frac{2}{(N_m - N_d - 1) (N_m - N_d - 4)}\]

\[B = \frac{N_m - N_d - 2}{(N_m - N_d - 1)(N_m - N_d - 4)}\]

where \(N_d\) is the number of data points, \(N_p\) is the number of parameters, \(N_m\) is the number of mock realisations for calculating the covariance, and \(\Sigma\) is the mock-based covariance matrix.

References

[1]

Percival, W. J. et al. 2014, MNRAS, 1312.4841

validate_input()[source]: Validate the input parameters.

class meer21cm.inference.SamplerEmcee(ps_dict: dict, data_vector: ndarray, data_covariance: ndarray, params_name: list[str], params_prior: list[tuple[str, float, float]], nwalkers: int, nsteps: int, nthreads: int, observables: list[str] = ['cross'], mp_backend: str = 'multiprocessing', save: bool = False, save_filename: str | None = None, save_model_blobs: bool = False, do_hartlap_correction: bool = False, num_mocks: int | None = None, do_percival_correction: bool = False)[source]

Bases: SamplerBase

A sampler class to perform model fitting using the emcee sampler.

Parameters:

ps_dict (dict) – A dictionary of model fitting inputs.
data_vector (np.ndarray) – The data vector.
data_covariance (np.ndarray) – The data covariance matrix.
params_name (list[str]) – The names of the parameters. Must match the attribute names of the meer21cm.power.PowerSpectrum instance.
params_prior (list[tuple[str, float, float]]) – The prior distribution of the parameters. The first element of the tuple is the distribution type, which can be “uniform” or “gaussian”. The second and third elements are the parameters of the distribution. For uniform distribution, the second and third elements are the lower and upper bounds of the distribution. For gaussian distribution, the second and third elements are the mean and standard deviation of the distribution.
nwalkers (int) – The number of random walkers.
nsteps (int) – The number of sampling steps.
nthreads (int) – The number of parallel threads.
observables (list[str], default ["cross"]) – The observables to fit. Must be a list of “cross”, “hi-auto”, or “gg-auto”.
mp_backend (str, default "multiprocessing") – The backend to use for the multiprocessing. Can be “multiprocessing” or “mpi”.
save (bool, default False) – Whether to save the fitting results while running the sampler.
save_filename (str | None, default None) – The filename to save the fitting results when save is True. If not provided, the results will not be saved.
save_model_blobs (bool, default False) – Whether to save the model blobs. The model blobs are the model vectors listed in the observables argument.
do_hartlap_correction (bool, default False) – Whether to apply the Hartlap correction to the data covariance matrix.
num_mocks (int, default None) – The number of mock data vectors to use for the Hartlap correction. If not provided, the Hartlap correction will not be applied.
do_percival_correction (bool, default False) – Whether to apply the Percival correction to the data covariance matrix.

Attributes:

data_covariance: The data covariance matrix.
do_hartlap_correction: Whether to apply the Hartlap correction to the data covariance matrix.
do_percival_correction: Whether to apply the Percival correction to the data covariance matrix.
hartlap_factor: The Hartlap factor.
inverse_covariance: The matrix inverse of the data covariance.
ndim: The number of parameters to fit.
num_mocks: The number of mock data vectors to use for the Hartlap and Percival corrections.
percival_factor: The Percival factor.

Methods

`compute_log_likelihood`(params)	Calculate the log likelihood from the input parameter values.
`get_backend`()	Get the backend of the sampler.
`get_blobs`([sampler])	Get the blobs from the sampler.
`get_log_prob`([sampler])	Get the log probability from the sampler.
`get_model_instance`(params)	Get the model instance from the input `ps_dict` and the updated parameters `params`.
`get_model_vector`(params)	Calculate the model vector from the input parameter values.
`get_points`([sampler])	Get the chain from the sampler.
`log_likelihood`(params_values)	Calculate the log likelihood plus the log prior from the input parameter values.
`log_prior`(params_values)	Calculate the log prior from the input parameter values.
`log_prior_gaussian`(value, mean, sigma)	Calculate the log prior for a gaussian distribution.
`log_prior_uniform`(value, low, high)	Calculate the log prior for a flat (uniform) distribution.
`run`(resume[, progress, start_coord, run_sampler])	Run the sampler.
`validate_input`()	Validate the input parameters.

get_backend() → HDFBackend[source]: Get the backend of the sampler.

get_blobs(sampler: EnsembleSampler | None = None) → ndarray[source]

Get the blobs from the sampler.

Parameters:: sampler (emcee.EnsembleSampler | None, default None) – The sampler object. If not provided, the sampler will use the backend to get the blobs (only works if save is True).
Returns:: blobs – The model vectorblobs.
Return type:: np.ndarray

get_log_prob(sampler: EnsembleSampler | None = None) → ndarray[source]

Get the log probability from the sampler.

Parameters:: sampler (emcee.EnsembleSampler | None, default None) – The sampler object. If not provided, the sampler will use the backend to get the log probability (only works if save is True).
Returns:: log_prob – The log probability. The shape is (nwalkers, nsteps).
Return type:: np.ndarray

get_points(sampler: EnsembleSampler | None = None) → ndarray[source]

Get the chain from the sampler.

Parameters:: sampler (emcee.EnsembleSampler | None, default None) – The sampler object. If not provided, the sampler will use the backend to get the chain (only works if save is True).
Returns:: chain – The chain. The shape is (nwalkers, nsteps, ndim).
Return type:: np.ndarray

log_likelihood(params_values)[source]

Calculate the log likelihood plus the log prior from the input parameter values.

Parameters:

params_values (list[float] | np.ndarray) – The updated values of the parameters defined in the params_name argument.

Returns:

log_likelihood (float) – The log likelihood plus the log prior.
model_vector (np.ndarray, optional) – The model vector. Only returned when save_model_blobs is True.

log_prior(params_values)[source]

Calculate the log prior from the input parameter values.

Parameters:: params_values (list[float] | np.ndarray) – The updated values of the parameters defined in the params_name argument.
Returns:: log_prior – The log prior.
Return type:: float

log_prior_gaussian(value, mean, sigma)[source]

Calculate the log prior for a gaussian distribution.

Parameters:

value (float) – The value to calculate the log prior for.
mean (float) – The mean of the distribution.
sigma (float) – The standard deviation of the distribution.

Returns:

log_prior – The log prior.

Return type:

float

log_prior_uniform(value, low, high)[source]

Calculate the log prior for a flat (uniform) distribution.

Parameters:

value (float) – The value to calculate the log prior for.
low (float) – The lower bound of the distribution.
high (float) – The upper bound of the distribution.

Returns:

log_prior – The log prior.

Return type:

float

property ndim: The number of parameters to fit.

run(resume: bool, progress: bool = True, start_coord: ndarray | None = None, run_sampler: bool = True)[source]

Run the sampler.

Parameters:

resume (bool) – Whether to resume the sampler from the last saved state. If resume is False and the save file already exists, the sampler will reset the save file.
progress (bool, default True) – Whether to show the progress bar.
start_coord (np.ndarray | None, default None) – The initial coordinates of the walkers. If not provided, the sampler will randomly initialize the walkers.
run_sampler (bool, default True) – Whether to run the sampler. If False, the sampler will only initialize the walkers and return the sampler object.

Returns:

sampler – The sampler object.

Return type:

emcee.EnsembleSampler

class meer21cm.inference.SamplerNautilus(ps_dict: dict, data_vector: ndarray, data_covariance: ndarray, params_name: list[str], params_prior: list[tuple[str, float, float]], observables: list[str] = ['cross'], n_live_points: int = 2000, f_live: float = 0.01, n_shell: int = 1, n_eff: int = 10000, save: bool = False, save_filename: str | None = None, mp_backend: str = 'multiprocessing', nthreads: int = 1, timeout: float = inf, save_model_blobs: bool = False, do_hartlap_correction: bool = False, num_mocks: int | None = None, do_percival_correction: bool = False)[source]

Bases: SamplerBase

A sampler class to perform model fitting using the nautilus sampler.

Parameters:

ps_dict (dict) – A dictionary of model fitting inputs.
data_vector (np.ndarray) – The data vector.
data_covariance (np.ndarray) – The data covariance matrix.
params_name (list[str]) – The names of the parameters. Must match the attribute names of the meer21cm.power.PowerSpectrum instance.
params_prior (list[tuple[str, float, float]]) – The prior distribution of the parameters. The first element of the tuple is the distribution type, which can be “uniform” or “gaussian”. The second and third elements are the parameters of the distribution. For uniform distribution, the second and third elements are the lower and upper bounds of the distribution. For gaussian distribution, the second and third elements are the mean and standard deviation of the distribution.
observables (list[str], default ["cross"]) – The observables to fit. Must be a list of “cross”, “hi-auto”, or “gg-auto”.
n_live_points (int, default 2000) – The number of live points.
f_live (float, default 0.01) – The fraction of live points to keep.
n_shell (int, default 1) – The number of shells to use.
n_eff (int, default 10000) – The effective number of samples.
save (bool, default False) – Whether to save the fitting results while running the sampler.
save_filename (str | None, default None) – The filename to save the fitting results when save is True. If not provided, the results will not be saved.
mp_backend (str, default "multiprocessing") – The backend to use for the multiprocessing. Can be “multiprocessing” or “mpi”.
nthreads (int, default 1) – The number of parallel threads.
timeout (float, default np.inf) – The timeout for the sampler. If the sampler runs longer than the timeout, it will be terminated.
save_model_blobs (bool, default False) – Whether to save the model blobs. The model blobs are the model vectors listed in the observables argument.
hartlap_correction (bool, default False) – Whether to apply the Hartlap correction to the data covariance matrix.
num_mocks (int, default None) – The number of mock data vectors to use for the Hartlap correction. If not provided, the Hartlap correction will not be applied.
percival_correction (bool, default False) – Whether to apply the Percival correction to the data covariance matrix.

Attributes:

data_covariance: The data covariance matrix.
do_hartlap_correction: Whether to apply the Hartlap correction to the data covariance matrix.
do_percival_correction: Whether to apply the Percival correction to the data covariance matrix.
hartlap_factor: The Hartlap factor.
inverse_covariance: The matrix inverse of the data covariance.
num_mocks: The number of mock data vectors to use for the Hartlap and Percival corrections.
percival_factor: The Percival factor.

Methods

`compute_log_likelihood`(params)	Calculate the log likelihood from the input parameter values.
`get_model_instance`(params)	Get the model instance from the input `ps_dict` and the updated parameters `params`.
`get_model_vector`(params)	Calculate the model vector from the input parameter values.
`get_nautilus_prior`()	Generate the `nautilus.Prior` objects from the input parameters.
`get_posterior`([sampler])	Get the posterior from the sampler.
`run`(resume[, progress, run_sampler])	Run the sampler.
`validate_input`()	Validate the input parameters.

get_nautilus_prior()[source]: Generate the nautilus.Prior objects from the input parameters.

get_posterior(sampler: Sampler | None = None)[source]

Get the posterior from the sampler.

Parameters:: sampler (nautilus.Sampler | None, default None) – The sampler object. If not provided, the sampler will use the backend to get the posterior (only works if save is True).
Returns:: posterior – The posterior. The first element is the points. The second element is the log weights. The third element is the log likelihoods. If save_model_blobs is True, the fourth element is the model blobs.
Return type:: tuple

run(resume: bool, progress: bool = True, run_sampler: bool = True)[source]

Run the sampler.

Parameters:

resume (bool) – Whether to resume the sampler from the last saved state.
progress (bool, default True) – Whether to show the progress bar.
run_sampler (bool, default True) – Whether to run the sampler. If False, the sampler will only initialize the walkers and return the sampler object.

Returns:

sampler – The sampler object.

Return type:

nautilus.Sampler

meer21cm.inference.extract_model_fitting_inputs(ps: PowerSpectrum)[source]

Extract a dictionary of model fitting inputs from a meer21cm.power.PowerSpectrum instance.

Parameters:: ps (meer21cm.power.PowerSpectrum) – A meer21cm.power.PowerSpectrum instance.
Returns:: ps_dict – A dictionary of model fitting inputs.
Return type:: dict

meer21cm.io module

Module for reading and pre-processing MeerKLASS maps.

meer21cm.io.cal_freq(ch_id, band='L', nu_min=None, delta_nu=None)[source]

returns the centre of the frequency channel for channel id ch_id of the meerkat telescope.

Parameters:

ch_id (int.) – The channel id.
band (str, default 'L'.) – Frequency band, can either be ‘L’ or ‘UHF’. Retrieves default MeerKAT setting. If nu_min and delta_nu are passed, the default settings are overridden.
nu_min (float, default 856.0*1e6 Hz.) – The lower end of the frequency range.
delta_nu (float, default 0.208984375*1e6 Hz.) – The channel bandwidth.

Returns:

freq – The frequency of the channel.

Return type:

float.

meer21cm.io.filter_incomplete_los(map_intensity, map_has_sampling, map_weight, map_pix_counts, los_axis=-1, soft_mask=False, threshold_instead_of_filter=None)[source]

Filter the map so that along the line-of-sight, only pixels that has sampling at every channel gets selected.

If soft_mask is True, instead of filtering out incomplete los, the filtering is applied by checking the maximum sampling fraction along the los, and the pixels with less than the maximum sampling fraction are masked.

If threshold_instead_of_filter is given, instead of filtering out incomplete los, the filtering is applied by checking the maximum sampling fraction along the los, and the pixels with less than the maximum sampling fraction are masked.

Parameters:

map_intensity (array.) – The input map.
map_has_sampling (boolean array.) – Whether the pixel has measurement.
map_weight (array.) – the pixel weights.
map_pix_counts (array.) – The channel bandwidth.
los_axis (int, default -1.) – which axis is the los.
soft_mask (boolean, default False.) – whether to apply soft masking.
threshold_instead_of_filter (float, default None.) – if given, instead of filtering out incomplete los, the filtering is applied by checking the sampling fraction along the los, and pixels with less than the threshold are masked.

Returns:

map_intensity (array.) – map after pixels that have incomplete los sampling are masked.
map_has_sampling (array.) – sampling after pixels that have incomplete los sampling are masked.
map_weight (array.) – weights after pixels that have incomplete los sampling are masked.
map_pix_counts (array.) – counts after pixels that have incomplete los sampling are masked.

meer21cm.io.read_map(map_file, counts_file=None, nu_min=-inf, nu_max=inf, ch_start=1, los_axis=-1, band='L')[source]

Read fits files of MeerKLASS 4k L-band data into arrays.

Parameters:

map_file (str.) – The input map file.
counts_file (str, default None.) – The input pixel counts file.
nu_min (float, default -np.inf.) – The lower end of frequency cut.
nu_max (float, default np.inf.) – The higher end of freuqency cut.
ch_start (int, default 1.) – The starting channel of the data.
los_axis (int, default -1.) – which axis is the los.

Returns:

map_data (array.) – The map data.
counts (array.) – The number of sampling for each pixel. If no counts_file is specified, it is the same as map_has_sampling.
map_has_sampling (boolean array.) – Whether the pixels are samplied.
ra (array.) – The RA coordinates of each pixel
dec (array.) – The Dec coordinates of each pixel
nu (array.) – The frequencies of each channel in the data.
wproj (astropy.wcs.WCS object.) – The two-dimensional wcs object for the map.

meer21cm.io.read_pickle(pickle_file, nu_min=-inf, nu_max=inf, los_axis=-1, data_column='map', counts_column='hit', freq_column='freq', wcs_column='wcs')[source]

Read pickle file of MeerKLASS UHF-band data into arrays. The file format requires the following keys: - map: the map data. - hit: the number of sampling for each pixel. - freq: the frequencies of each channel in the data in MHz. meer21cm then converts it to Hz. - wcs: the astropy.wcs.WCS object for the map.

Parameters:

pickle_file (str.) – The input pickle file.
nu_min (float, default -np.inf.) – The lower end of frequency cut. Channels below this frequency will be thrown away.
nu_max (float, default np.inf.) – The higher end of freuqency cut. Channels above this frequency will be thrown away.
los_axis (int, default -1.) – which axis is the los.
data_column (str, default "map".) – The column name of the map data.
counts_column (str, default "hit".) – The column name of the number of sampling for each pixel.
freq_column (str, default "freq".) – The column name of the frequencies of each channel in the data.
wcs_column (str, default "wcs".) – The column name of the astropy.wcs.WCS object for the map.

Returns:

map_data – The map data.

Return type:

array.

meer21cm.mock module

This module contains the classes for generating mock intensity mapping data, galaxy catalogues for cross-correlation, and HI galaxy catalogues for HI emission line observations.

class meer21cm.mock.HIGalaxySimulation(no_vel=True, tf_slope=None, tf_zero=None, halo_model=None, hi_mass_from='hod', himf_pars=None, num_ch_ext_on_each_side=5, **params)[source]

Bases: MockSimulation

The class for generating mock HI galaxy catalogues for HI emission line observations.

Parameters:

no_vel (bool, default True) – If True, HI sources will have no velocity width, and the total flux will be allocated to the central frequency channel.
tf_slope (float, default None) – The slope of the Tully-Fisher relation when no_vel is False.
tf_zero (float, default None) – The intercept of the Tully-Fisher relation when no_vel is False.
halo_model (halomod.TracerHaloModel, default None) – The halo model to be used. Only used if hi_mass_from is “hod”.
hi_mass_from (str, default "hod") – The method for calculating the HI mass of the mock tracers. Can be “hod” or “himf”.
himf_pars (dict, default None) – The parameters for the HI mass function. Only used if hi_mass_from is “himf”.
num_ch_ext_on_each_side (int, default 5) – The number of frequency channels to be extended on each side of the central frequency channel when generating the HI profile. A good rule of thumb is 200km/s over self.vel_resol.
**params (dict) – Additional parameters to be passed to the base class meer21cm.mock.MockSimulation.

Attributes:

As: The amplitude of the primordial power spectrum for the true (fitted) cosmology.
W_HI: A binary window for whether a pixel has been sampled
alpha_AP: The anisotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel} / \alpha_{\perp}\).
alpha_iso: The isotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel}^{1/3} \alpha_{\perp}^{2/3}\).
alpha_parallel: The line of sight Alcock–Paczynski effect parameter.
alpha_perp
auto_power_3d_1: The 3D power spectrum of the first tracer.
auto_power_3d_2: The 3D power spectrum of the second tracer.
auto_power_matter_model: The model matter power spectrum with RSD effects.
auto_power_matter_model_r: The model matter power spectrum in real space (without RSD).
auto_power_tracer_1_model: The 3D model power spectrum for the first tracer.
auto_power_tracer_1_model_noobs: The model power spectrum for the first tracer without observational effects.
auto_power_tracer_2_model: The 3D model power spectrum for the second tracer.
auto_power_tracer_2_model_noobs: The model power spectrum for the second tracer without observational effects.
average_hi_temp: The average HI brightness temperature in Kelvin at central redshift self.z.
average_model_hi_temp: Calculate the average HI brightness temperature in the map cube, taking care of redshift evolution and map sampling.
backend: Which backend to use for computing the matter power.
batch_number: Number of sequential batches used by gridding routines.
beam_image: Returns the beam image projected onto the sky map for the input beam model.
beam_model: The name of the beam function.
beam_type: The beam type that can be either be isotropic or anisotropic.
beam_unit: The unit of input beam size parameter sigma
beam_window_ch: Per-channel spherical-harmonic beam window \(B_\ell\) for HEALPix hp.smoothing workflows.
box_buffkick: The buffer kick for the box on each side when gridding.
box_len: The length of all sides of the box in Mpc.
box_ndim: The number of grids along each side of the enclosing box.
box_origin: The coordinate of the origin of the box in Mpc.
box_resol: The grid length of each side of the enclosing box in Mpc.
box_voxel_redshift: The redshift of each voxel in the rectangular box.
ch_id_gal: The channel id (0-indexed) of each galaxy in the catalogue for cross-correlation.
cold: If True (recommended), the matter power spectrum is the CDM ps.
compensate: Whether the gridded fields are compensated according to the mass assignment scheme.
cosmo: A shortcut to the astropy.cosmology.Cosmology object for the true cosmology.
cospar_fiducial
cospar_true
counts: The number of hits per pixel for the map data
counts_in_box: The counts of the map cube voxels in the rectangular box.
cross_coeff: The cross-correlation coefficient between the two tracers.
cross_power_3d: The 3D cross power spectrum between the two tracers.
cross_power_tracer_model: The 3D model cross power spectrum between the two tracers.
cross_power_tracer_model_noobs: The model power spectrum for the cross correlation between the two tracers without observational effects.
data: The map data
dec_gal: The declination of each galaxy in the catalogue for cross-correlation.
dec_map: The declination of each pixel in the map.
dec_mock_tracer: The Dec coordinate of mock galaxies on the sky
discrete_base_field: Which tracer (1 or 2) field to sample the discrete tracer positions.
discrete_source_dndz: The redshift kernel of the discrete tracer sources.
downres_factor_radial: The down-sampling factor for the radial direction of the rectangular box for gridding.
downres_factor_transverse: The down-sampling factor for the transverse direction of the rectangular box for gridding.
dvdf: velocity resolution per unit frequency on average, in km/s/Hz
dvdf_ch: velocity resolution per unit frequency in each channel, in km/s/Hz
expfactor: The expansion factor
f_growth_fiducial: The growth factor at self.expfactor for the fiducial cosmology.
f_growth_true: The growth factor at self.expfactor for the true (fitted) cosmology.
fiducial_cosmology
field_1: The density field of the first tracer.
field_2: The density field of the second tracer.
flat_sky: Whether to use flat sky approximation.
flat_sky_padding: Pad the rectangular box in the flat sky approximation.
fog_profile: The shape of the finger-of-god profile to be used in the model calculation.
fourier_field_1: The Fourier transform of the density field of the first tracer.
fourier_field_2: The Fourier transform of the density field of the second tracer.
freq_gal: The 21cm line frequency for each galaxy in Hz.
freq_resol: frequency resolution in Hz
h: The Hubble constant for the true (fitted) cosmology.
halo_mass_mock_tracer: The halo mass of the mock tracers in log10 M_sun/h.
halo_model: A halomod.TracerHaloModel object for storing the halo model.
hi_mass_from: Methods for calculating HI mass.
hi_mass_mock_tracer: The HI mass of the mock tracers in log10 M_sun.
hi_profile_mock_tracer: The emission line profiles of the mock tracers in Jansky
highres_sim: Optional integer upsampling factor for HIGalaxySimulation.propagate_hi_profile_to_map().
hp_nside: HEALPix \(N_{side}\) for HEALPix-backed specifications.
include_beam: Whether the beam attenuation is included in the model calculation.
interlace_shift: The length in the unit of grid cell size for shifting the gridded field for interlacing.
k_mode: The fiducial (observed) mode of the 3D k-vector.
k_nyquist: The Nyquist frequency of the 3D box along each axis.
k_para: The fiducial parallel k-mode of the 3D box.
k_perp: The fiducial perpendicular k-vector of the 3D box.
k_vec: The 3D k-vector of the box.
kaiser_rsd: Whether RSD is included in the simulation and model calculation.
kmax: The maximum k in Mpc^-1 for calculating matter power.
kmin: The minimum k in Mpc^-1 for calculating matter power.
kmode: The input kmode for the model calculation.
los_resol_in_mpc: effective frequency resolution in Mpc for the fiducial cosmology.
map_has_sampling: A binary window for whether a pixel has been sampled
map_unit_type: The type of the map unit.
matter_power_spectrum_fnc: Interpolation function for the real-space isotropic matter power spectrum.
maximum_sampling_channel: Returns the index of the frequency channel with the maximum sampling on the sky map.
mean_center_1: Whether field_1 needs to be mean centered
mean_center_2: Whether field_2 needs to be mean centered
mock_amp_1: The overall amplitude of the mock tracer field 1 in each grid.
mock_amp_2: The overall amplitude of the mock tracer field 2 in each grid.
mock_inside_range: Whether the mock galaxies are inside the survey area and frequency range
mock_kaiser_field_k_matter: The Kaiser rsd effect correction for the mock matter field in k-space.
mock_kaiser_field_k_tracer_1: The Kaiser rsd effect correction for the mock tracer field 1 in k-space.
mock_kaiser_field_k_tracer_2: The Kaiser rsd effect correction for the mock tracer field 2 in k-space.
mock_matter_field: The simulated dark matter density field in redshift space.
mock_matter_field_r: The simulated dark matter density field in real space.
mock_tracer_field_1: The simulated tracer field 1 in redshift space with unit if mock_amp_1 or mean_amp_1 is given.
mock_tracer_field_1_r: The simulated tracer field 1 unitsless density contrast in real space.
mock_tracer_field_2: The simulated tracer field 2 in redshift space with unit if mock_amp_2 or mean_amp_2 is given.
mock_tracer_field_2_r: The simulated tracer field 2 unitsless density contrast in real space.
mock_tracer_position_in_box: The simulated tracer positions in the box in real space.
mock_tracer_position_in_radecz: The simulated tracer positions projected onto the grid.
mock_velocity_u_matter: The normalised peculiar velocity field in real space, defined as
mock_velocity_u_tracer_1: The normalised peculiar velocity field used for the first tracer.
mock_velocity_u_tracer_2: The normalised peculiar velocity field used for the second tracer.
model_hi_temp_in_box: Based on the redshift dependence of Omega_HI, calculate the average HI brightness temperature for each grid in the rectangular box.
mu_mode: The fiducial (observed) mu values of each k-mode so that \(k_\parallel = k imes \mu\).
mumode: The mu values of each 3D k-mode.
neutrino_mass: The neutrino mass for the true (fitted) cosmology.
no_vel: If True, HI sources will have no velocity width.
ns: The primordial power spectrum spectral index for the true (fitted) cosmology.
nu: The input frequencies of the survey
num_discrete_source: The total number of discrete tracer sources.
num_kpoints: The number of k points for calculating matter power.
num_particle_per_pixel: The number of random sampling particles for each sky map pixel.
num_pix_x: The number of pixels along the first map axis.
num_pix_y: The number of pixels along the second map axis.
omega_baryon: The baryon matter density parameter for the true (fitted) cosmology.
omega_cold: The cold matter density parameter for the true (fitted) cosmology.
omega_hi: The HI density as a function of redshift, over the critical density of the Universe at z=0, defined at each frequency channel so that self.omega_hi corresponds to self.z_ch.
omega_hi_z_func: Interpolate the input self.omega_hi at each frequency channel self.z_ch to a function of redshift.
omega_hi_z_mean: The mean HI density at the central redshift self.z, over the critical density of the Universe at z=0.
parallel_plane: Whether the mock field is generated in the parallel-plane limit.
pix_coor_in_box: The cartesian coordinate of the pixels in Mpc, shifted so that the origin is the origin of the enclosing box.
pix_coor_in_cartesian: The cartesian coordinate of the pixels in Mpc.
pix_resol: angular resolution of the map pixel in deg
pix_resol_in_mpc: angular resolution of the map pixel in Mpc for the fiducial cosmology.
pixel_area: angular area of the map pixel in deg^2
pixel_id: Sparse HEALPix pixel indices (RING, nest=False).
precision: Floating precision flag.
ps_type: linear or nonlinear for the matter power.
ra_gal: The right ascension of each galaxy in the catalogue for cross-correlation.
ra_map: The right ascension of each pixel in the map.
ra_mock_tracer: The RA coordinate of mock galaxies on the sky
real_dtype: Active real floating dtype controlled by precision.
renorm_ps_1: The renormalization factor of the power spectrum of the first tracer.
renorm_ps_2: The renormalization factor of the power spectrum of the second tracer.
renorm_ps_cross: The renormalization factor of the cross power spectrum.
rot_mat_sky_to_box: The rotational matrix from spheircal cooridnate to regular box.
rsd_from_field: If True, the kaiser rsd effect is generated at field level by calculating the corresponding peculiar velocity field.
sampling_resol: The sampling resolution corresponding to the map-making/gridding of the density field.
seed: Seed value for RNG calls throughout the instance.
sigma_beam_ch: The input beam size parameter sigma for each channel.
sigma_beam_ch_in_mpc: The input beam size parameter in Mpc in each channel.
sigma_beam_in_mpc: The channel averaged beam size in Mpc
sigma_v_1: The velocity dispersion of the first tracer in km/s.
sigma_v_2: The velocity dispersion of the second tracer in km/s.
sigma_z_1: The redshift error of the first tracer.
sigma_z_2: The redshift error of the second tracer.
sound_horizon_drag_fiducial: The sound horizon at drag epoch for the fiducial cosmology.
sound_horizon_drag_true: The sound horizon at drag epoch for the true (fitted) cosmology.
survey_volume: Total survey volume in Mpc^3.
tf_slope: Slope of Tully-Fisher relation.
tf_zero: The intercept of the T-F relation.
tot_num_source_in_box: The total number of mock sources in the box needed to achieve self.num_discrete_source number of sources in the survey volume.
tracer_bias_1: The linear bias of the first tracer.
tracer_bias_2: The linear bias of the second tracer.
true_cosmology
unitless_1: Whether field_1 needs to be divided by its mean
unitless_2: Whether field_2 needs to be divided by its mean
vel_resol: velocity resolution on average in km/s
vel_resol_ch: velocity resolution of each channel in km/s
w0: The dark energy equation of state parameter for the true (fitted) cosmology.
w_HI: The weights per map pixel.
wa: The dark energy equation of state parameter for the true (fitted) cosmology.
weights_1: The weights of the first tracer in the rectangular grid.
weights_2: The weights of the second tracer in the rectangular grid.
weights_field_1: The weights of the first tracer in the density field.
weights_field_2: The weights of the second tracer in the density field.
weights_grid_1: The weights of the first tracer in the rectangular grid.
weights_grid_2: The weights of the second tracer in the rectangular grid.
weights_map_pixel: The weights per map pixel.
wproj: The WCS projection object for the map geometry.
x_mode: The mode of the 3D x-vector.
x_vec: The 3D x-vector of the box.
z: The effective centre redshift of the frequency range
z_as_func_of_comov_dist: Returns a function that returns the redshift for input comoving distance.
z_ch: The redshift of each frequency channel
z_gal: The redshifts of each galaxy in the catalogue for cross-correlation.
z_mock_tracer: The redshift of mock galaxies

Methods

`apply_taper_to_field`(field[, taper_func, axis])	Apply a taper to the field, by multiplying the taper function to the corresponding weights of the field.
`beam_attenuation`()	The beam attenuation factor.
`cal_rsd_power`(power_in_real_space, beta1, ...)	Calculate the redshift space power spectrum.
`check_is_map_noiselike_using_pca`(A_mat[, ...])	Use the source mixing matrix from eigendecomposition of the covariance matrix, project out the map data with more and more modes, and check if the variance of the residual map behaves like white noise.
`clean_cache`(attr)	Set the attributes to None.
`convolve_data`([kernel, data, weights, ...])	Convolve map data.
`create_white_noise_map`(sigma_N[, counts, ...])	Create a white noise map with the given standard deviation.
`deltav_to_deltar`(delta_v)	Convert a velocity interval delta_v to a comoving distance interval delta_r.
`deltaz_to_deltar`(delta_z)	Convert a redshift interval delta_z to a comoving distance interval delta_r.
`fog_gaussian`(sigma_r[, kmode, mumode])	The Gaussian finger-of-god profile.
`fog_lorentz`(sigma_r[, kmode, mumode])	The Lorentzian finger-of-god profile.
`fog_term`(sigma_r[, kmode, mumode])	The FoG term for the model calculation.
`gen_random_poisson_galaxy`([sel, num_g_rand, ...])	Generate a random galaxy catalogue from the map cube following the Poisson distribution.
`generate_full_healpix_map`([data, fill_value])	Generate a full healpix map from the data.
`get_1d_power`(power3d[, k1dbins, k1dweights, ...])	Bin the 3D power spectrum into 1D power spectrum.
`get_alpha_parallel`()	Calculate the line of sight Alcock–Paczynski effect parameter.
`get_alpha_perp`()	Calculate the transverse Alcock–Paczynski effect parameter.
`get_beam_image`([wproj, num_pix_x, ...])	Calculate the beam image projected onto the sky map for the input beam model.
`get_beam_window_ch`([ch_sel, cache, lmax, ...])	Build per-channel harmonic beam windows for HEALPix smoothing.
`get_cospar`(cosmology)	Generate a `CosmologyParameters` object from the input cosmology.
`get_counts_in_box`([partial_sel])	Grid the counts of the map cube voxels into the rectangular box, and return the effective counts per rectangular grid.
`get_cy_power`(power3d[, kperpbins, ...])	Bin the 3D power spectrum into cylindrical k_perp-k_para power spectrum.
`get_enclosing_box`([rot_mat])	invoke to calculate the box dimensions for enclosing all the map pixels.
`get_fourier_field_1`()	Calculate the Fourier transform of the density field of the first tracer.
`get_fourier_field_2`()	Calculate the Fourier transform of the density field of the second tracer.
`get_halo_mass_mock_tracer`()	Calculate the halo mass of the mock tracers.
`get_hi_mass_mock_tracer_hod`()	Calculate the HI mass of the mock tracers using the halo occupation distribution based on the HOD model stored in `self.halo_model`.
`get_hi_profile_mock_tracer`()	Calculate the emission line profiles of the mock tracers.
`get_jackknife_patches`(ra_patch_num, ...[, ...])	Split the map into roughly equal patches.
`get_matter_power_spectrum`()	Calculate the matter power spectrum, interpolate it, and save it into the class attribute matter_power_spectrum_fnc.
`get_mock_field_from_power`(power_array)	Generate a mock field that follows the input matter power spectrum `power_array`.
`get_mock_field_in_redshift_space`(delta_x, ...)	Generate the mock field in redshift space.
`get_mock_field_r`([bias])	Generate a mock field in real space that follows the input matter power spectrum `self.matter_power_spectrum_fnc` * bias**2.
`get_mock_kaiser_field_k`(field)	Generate the Kaiser rsd effect for the mock matter field in k-space.
`get_mock_matter_field`()	Generate the mock matter field in redshift space.
`get_mock_matter_field_r`()	Generate the mock matter field in real space.
`get_mock_tracer_field`(tracer_i)	Generate the mock tracer field in redshift space.
`get_mock_tracer_field_r`(tracer_i)	Generate the mock tracer field in real space.
`get_mock_tracer_position_in_box`(tracer_i[, ...])	Function to retrieve the tracer positions in redshift space.
`get_mock_tracer_position_in_radecz`()	Project the mock tracer positions in the rectangular box to the sky coordinates and redshifts.
`get_mock_velocity_u_field`(mock_field)	Generate the normalised peculiar velocity field in real space
`get_model_matter_power`()	Calculate the model matter power spectrum.
`get_model_matter_power_r`()	Calculate the model matter power spectrum in real space (without RSD).
`get_model_power_cross`()	Calculate the model cross power spectrum between the two tracers.
`get_model_power_i`(i)	Calculate the model power spectrum for the i-th tracer.
`get_model_power_noobs_cross`()	Calculate the model cross power spectrum between the two tracers without observational effects.
`get_model_power_noobs_i`(i)	Calculate the model power spectrum for the i-th tracer without observational effects.
`get_n_bar_correction`()	Calculate the number density correction for the galaxy catalogue.
`get_sound_horizon_drag`(Omh2, Obh2, Onuh2)	Calculate the sound horizon at drag epoch for a set of cosmological parameters.
`get_weights_none_to_one`(attr_name)	Get the weights, and if it is None, convert it to 1.0 of size of kmode.
`get_z_as_func_of_comov_dist`()	Calculate an array of comoving distances with redshifts, and construct a function that returns the redshift for input comoving distance.
`grid_data_to_field`([flat_sky, partial_sel])	Grid the stored data map to a rectangular grid field.
`grid_field_to_sky_map`(field[, average, ...])	Grid a field in the rectangular box onto the sky.
`grid_gal_to_field`([radecfreq, flat_sky])	Grid the galaxy catalogue to a rectangular grid field.
`gridding_compensation`()	The sampling window function to be compensated for the gridding mass assignment scheme.
`map_sampling`([sampling_resol, p])	The sampling window function from the map cube to be convolved with model power spectrum.
`propagate_field_k_to_model`()	Use field k-modes for the model, taking into account the Alcock–Paczynski effect.
`propagate_hi_profile_to_map`([...])	Project the `hi_profile_mock_tracer` onto sky maps and convolve with the beam (if `beam`).
`propagate_mock_field_to_data`(field[, beam, ...])	Grid the mock tracer field onto the sky map cube.
`propagate_mock_tracer_to_gal_cat`([trim])	Propagate the mock tracer positions to the galaxy data catalogue.
`ra_dec_z_for_coord_in_box`(pos_in_box)	Convert the coordinates in the box to ra, dec, z, and also return the comoving distance to the observer for each point.
`read_from_fits`()	Read in a FITS file for the map data and hit counts.
`read_from_pickle`()	Read in a pickle file for cross-correlation and save the data into the class attributes.
`read_gal_cat`([ra_col, dec_col, z_col, trim])	Read in a galaxy catalogue for cross-correlation and save the data into the class attributes.
`set_corr_type`(corr_type, tracer_indx)	A utility function to help decide whether a tracer field needs to be mean centred, renormalised by its mean, and shot noise removed.
`set_radecnu_bounds_from_map`()	Set `ra_range`, `dec_range`, `nu_min`, and `nu_max` from the loaded `_ra_map`, `_dec_map`, and `nu` (tight RA interval, declination min/max, frequency min/max).
`trim_gal_to_range`()	Trim the galaxy catalogue to the specified range.
`trim_map_to_range`()	Trim the map to the specified range.
`use_flat_sky_box`([flat_sky_padding])	Use flat sky approximation to calculate the box dimensions.

get_halo_mass_mock_tracer()[source]

Calculate the halo mass of the mock tracers.

In this case, the discrete tracer is considered to be the halos instead of galaxies, so corroespondingly you should set self.num_discrete_source to the number of halos, and self.tracer_bias_2 to the effective bias of the halo number density field.

The sampling is then done by using the halo mass function stored in self.halo_model. See hmf package for more details about the halo mass function.

get_hi_mass_mock_tracer_hod()[source]

Calculate the HI mass of the mock tracers using the halo occupation distribution based on the HOD model stored in self.halo_model.

See halomod package for more details about the HOD model.

get_hi_profile_mock_tracer()[source]

Calculate the emission line profiles of the mock tracers.

The profiles are calculated using the HI mass of the mock tracers, and the Tully-Fisher relation stored in self.tf_slope and self.tf_zero.

property halo_mass_mock_tracer: The halo mass of the mock tracers in log10 M_sun/h.

property halo_model: A halomod.TracerHaloModel object for storing the halo model.

property hi_mass_from: Methods for calculating HI mass. Can either be ‘hod’ or ‘himf’.

property hi_mass_mock_tracer: The HI mass of the mock tracers in log10 M_sun.

property hi_profile_mock_tracer: The emission line profiles of the mock tracers in Jansky

property no_vel: If True, HI sources will have no velocity width.

propagate_hi_profile_to_map(return_highres=False, beam=True, beam_image=None)[source]: Project the hi_profile_mock_tracer onto sky maps and convolve with the beam (if beam). If return_highres, the returned map is the higher resolution map specified by highres_sim. If not, the map will be downsampled to the original resolution.

property tf_slope: Slope of Tully-Fisher relation. See meer21cm.util.tully_fisher().

property tf_zero: The intercept of the T-F relation. See meer21cm.util.tully_fisher().

class meer21cm.mock.MockSimulation(density='lognormal', num_discrete_source=100, discrete_base_field=2, highres_sim=None, parallel_plane=True, rsd_from_field=False, discrete_source_dndz=<function ones_like>, mock_amp_1=None, mock_amp_2=None, **params)[source]

Bases: PowerSpectrum

The class for generating mock intensity mapping data cube and galaxy catalogues for cross-correlation.

Parameters:

density (str, default "lognormal") – The density distribution of the mock field. Can be “lognormal” or “gaussian”.
num_discrete_source (int, default 100) – The number of discrete tracer sources.
discrete_base_field (int, default 2) – The tracer (1 or 2) field to sample the discrete tracer positions.
highres_sim (int, default None) – For WCS maps, this factor is deprecated in propagate_mock_field_to_data() (a warning is issued and native resolution is used). It still controls the optional high-resolution path in HIGalaxySimulation.propagate_hi_profile_to_map().
parallel_plane (bool, default True) – Whether to simulate the mock field in the parallel-plane limit (i.e. k_z = k_parallel). Only used if rsd_from_field is True, otherwise parallel-plane is hard-coded to True.
rsd_from_field (bool, default False) – Whether to generate the RSD effect at field level. If False, the RSD effect is generated at power spectrum level.
discrete_source_dndz (function, default np.ones_like) – The redshift distribution of the discrete tracer sources. Must be a function of redshift. Note that the overall number of discrete sources is given by self.num_discrete_source, so only the shape of the dndz is used, not the normalization.
**params (dict) – Additional parameters to be passed to the base class meer21cm.power.PowerSpectrum.

Attributes:

As: The amplitude of the primordial power spectrum for the true (fitted) cosmology.
W_HI: A binary window for whether a pixel has been sampled
alpha_AP: The anisotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel} / \alpha_{\perp}\).
alpha_iso: The isotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel}^{1/3} \alpha_{\perp}^{2/3}\).
alpha_parallel: The line of sight Alcock–Paczynski effect parameter.
alpha_perp
auto_power_3d_1: The 3D power spectrum of the first tracer.
auto_power_3d_2: The 3D power spectrum of the second tracer.
auto_power_matter_model: The model matter power spectrum with RSD effects.
auto_power_matter_model_r: The model matter power spectrum in real space (without RSD).
auto_power_tracer_1_model: The 3D model power spectrum for the first tracer.
auto_power_tracer_1_model_noobs: The model power spectrum for the first tracer without observational effects.
auto_power_tracer_2_model: The 3D model power spectrum for the second tracer.
auto_power_tracer_2_model_noobs: The model power spectrum for the second tracer without observational effects.
average_hi_temp: The average HI brightness temperature in Kelvin at central redshift self.z.
average_model_hi_temp: Calculate the average HI brightness temperature in the map cube, taking care of redshift evolution and map sampling.
backend: Which backend to use for computing the matter power.
batch_number: Number of sequential batches used by gridding routines.
beam_image: Returns the beam image projected onto the sky map for the input beam model.
beam_model: The name of the beam function.
beam_type: The beam type that can be either be isotropic or anisotropic.
beam_unit: The unit of input beam size parameter sigma
beam_window_ch: Per-channel spherical-harmonic beam window \(B_\ell\) for HEALPix hp.smoothing workflows.
box_buffkick: The buffer kick for the box on each side when gridding.
box_len: The length of all sides of the box in Mpc.
box_ndim: The number of grids along each side of the enclosing box.
box_origin: The coordinate of the origin of the box in Mpc.
box_resol: The grid length of each side of the enclosing box in Mpc.
box_voxel_redshift: The redshift of each voxel in the rectangular box.
ch_id_gal: The channel id (0-indexed) of each galaxy in the catalogue for cross-correlation.
cold: If True (recommended), the matter power spectrum is the CDM ps.
compensate: Whether the gridded fields are compensated according to the mass assignment scheme.
cosmo: A shortcut to the astropy.cosmology.Cosmology object for the true cosmology.
cospar_fiducial
cospar_true
counts: The number of hits per pixel for the map data
counts_in_box: The counts of the map cube voxels in the rectangular box.
cross_coeff: The cross-correlation coefficient between the two tracers.
cross_power_3d: The 3D cross power spectrum between the two tracers.
cross_power_tracer_model: The 3D model cross power spectrum between the two tracers.
cross_power_tracer_model_noobs: The model power spectrum for the cross correlation between the two tracers without observational effects.
data: The map data
dec_gal: The declination of each galaxy in the catalogue for cross-correlation.
dec_map: The declination of each pixel in the map.
dec_mock_tracer: The Dec coordinate of mock galaxies on the sky
discrete_base_field: Which tracer (1 or 2) field to sample the discrete tracer positions.
discrete_source_dndz: The redshift kernel of the discrete tracer sources.
downres_factor_radial: The down-sampling factor for the radial direction of the rectangular box for gridding.
downres_factor_transverse: The down-sampling factor for the transverse direction of the rectangular box for gridding.
dvdf: velocity resolution per unit frequency on average, in km/s/Hz
dvdf_ch: velocity resolution per unit frequency in each channel, in km/s/Hz
expfactor: The expansion factor
f_growth_fiducial: The growth factor at self.expfactor for the fiducial cosmology.
f_growth_true: The growth factor at self.expfactor for the true (fitted) cosmology.
fiducial_cosmology
field_1: The density field of the first tracer.
field_2: The density field of the second tracer.
flat_sky: Whether to use flat sky approximation.
flat_sky_padding: Pad the rectangular box in the flat sky approximation.
fog_profile: The shape of the finger-of-god profile to be used in the model calculation.
fourier_field_1: The Fourier transform of the density field of the first tracer.
fourier_field_2: The Fourier transform of the density field of the second tracer.
freq_gal: The 21cm line frequency for each galaxy in Hz.
freq_resol: frequency resolution in Hz
h: The Hubble constant for the true (fitted) cosmology.
highres_sim: Optional integer upsampling factor for HIGalaxySimulation.propagate_hi_profile_to_map().
hp_nside: HEALPix \(N_{side}\) for HEALPix-backed specifications.
include_beam: Whether the beam attenuation is included in the model calculation.
interlace_shift: The length in the unit of grid cell size for shifting the gridded field for interlacing.
k_mode: The fiducial (observed) mode of the 3D k-vector.
k_nyquist: The Nyquist frequency of the 3D box along each axis.
k_para: The fiducial parallel k-mode of the 3D box.
k_perp: The fiducial perpendicular k-vector of the 3D box.
k_vec: The 3D k-vector of the box.
kaiser_rsd: Whether RSD is included in the simulation and model calculation.
kmax: The maximum k in Mpc^-1 for calculating matter power.
kmin: The minimum k in Mpc^-1 for calculating matter power.
kmode: The input kmode for the model calculation.
los_resol_in_mpc: effective frequency resolution in Mpc for the fiducial cosmology.
map_has_sampling: A binary window for whether a pixel has been sampled
map_unit_type: The type of the map unit.
matter_power_spectrum_fnc: Interpolation function for the real-space isotropic matter power spectrum.
maximum_sampling_channel: Returns the index of the frequency channel with the maximum sampling on the sky map.
mean_center_1: Whether field_1 needs to be mean centered
mean_center_2: Whether field_2 needs to be mean centered
mock_amp_1: The overall amplitude of the mock tracer field 1 in each grid.
mock_amp_2: The overall amplitude of the mock tracer field 2 in each grid.
mock_inside_range: Whether the mock galaxies are inside the survey area and frequency range
mock_kaiser_field_k_matter: The Kaiser rsd effect correction for the mock matter field in k-space.
mock_kaiser_field_k_tracer_1: The Kaiser rsd effect correction for the mock tracer field 1 in k-space.
mock_kaiser_field_k_tracer_2: The Kaiser rsd effect correction for the mock tracer field 2 in k-space.
mock_matter_field: The simulated dark matter density field in redshift space.
mock_matter_field_r: The simulated dark matter density field in real space.
mock_tracer_field_1: The simulated tracer field 1 in redshift space with unit if mock_amp_1 or mean_amp_1 is given.
mock_tracer_field_1_r: The simulated tracer field 1 unitsless density contrast in real space.
mock_tracer_field_2: The simulated tracer field 2 in redshift space with unit if mock_amp_2 or mean_amp_2 is given.
mock_tracer_field_2_r: The simulated tracer field 2 unitsless density contrast in real space.
mock_tracer_position_in_box: The simulated tracer positions in the box in real space.
mock_tracer_position_in_radecz: The simulated tracer positions projected onto the grid.
mock_velocity_u_matter: The normalised peculiar velocity field in real space, defined as
mock_velocity_u_tracer_1: The normalised peculiar velocity field used for the first tracer.
mock_velocity_u_tracer_2: The normalised peculiar velocity field used for the second tracer.
model_hi_temp_in_box: Based on the redshift dependence of Omega_HI, calculate the average HI brightness temperature for each grid in the rectangular box.
mu_mode: The fiducial (observed) mu values of each k-mode so that \(k_\parallel = k imes \mu\).
mumode: The mu values of each 3D k-mode.
neutrino_mass: The neutrino mass for the true (fitted) cosmology.
ns: The primordial power spectrum spectral index for the true (fitted) cosmology.
nu: The input frequencies of the survey
num_discrete_source: The total number of discrete tracer sources.
num_kpoints: The number of k points for calculating matter power.
num_particle_per_pixel: The number of random sampling particles for each sky map pixel.
num_pix_x: The number of pixels along the first map axis.
num_pix_y: The number of pixels along the second map axis.
omega_baryon: The baryon matter density parameter for the true (fitted) cosmology.
omega_cold: The cold matter density parameter for the true (fitted) cosmology.
omega_hi: The HI density as a function of redshift, over the critical density of the Universe at z=0, defined at each frequency channel so that self.omega_hi corresponds to self.z_ch.
omega_hi_z_func: Interpolate the input self.omega_hi at each frequency channel self.z_ch to a function of redshift.
omega_hi_z_mean: The mean HI density at the central redshift self.z, over the critical density of the Universe at z=0.
parallel_plane: Whether the mock field is generated in the parallel-plane limit.
pix_coor_in_box: The cartesian coordinate of the pixels in Mpc, shifted so that the origin is the origin of the enclosing box.
pix_coor_in_cartesian: The cartesian coordinate of the pixels in Mpc.
pix_resol: angular resolution of the map pixel in deg
pix_resol_in_mpc: angular resolution of the map pixel in Mpc for the fiducial cosmology.
pixel_area: angular area of the map pixel in deg^2
pixel_id: Sparse HEALPix pixel indices (RING, nest=False).
precision: Floating precision flag.
ps_type: linear or nonlinear for the matter power.
ra_gal: The right ascension of each galaxy in the catalogue for cross-correlation.
ra_map: The right ascension of each pixel in the map.
ra_mock_tracer: The RA coordinate of mock galaxies on the sky
real_dtype: Active real floating dtype controlled by precision.
renorm_ps_1: The renormalization factor of the power spectrum of the first tracer.
renorm_ps_2: The renormalization factor of the power spectrum of the second tracer.
renorm_ps_cross: The renormalization factor of the cross power spectrum.
rot_mat_sky_to_box: The rotational matrix from spheircal cooridnate to regular box.
rsd_from_field: If True, the kaiser rsd effect is generated at field level by calculating the corresponding peculiar velocity field.
sampling_resol: The sampling resolution corresponding to the map-making/gridding of the density field.
seed: Seed value for RNG calls throughout the instance.
sigma_beam_ch: The input beam size parameter sigma for each channel.
sigma_beam_ch_in_mpc: The input beam size parameter in Mpc in each channel.
sigma_beam_in_mpc: The channel averaged beam size in Mpc
sigma_v_1: The velocity dispersion of the first tracer in km/s.
sigma_v_2: The velocity dispersion of the second tracer in km/s.
sigma_z_1: The redshift error of the first tracer.
sigma_z_2: The redshift error of the second tracer.
sound_horizon_drag_fiducial: The sound horizon at drag epoch for the fiducial cosmology.
sound_horizon_drag_true: The sound horizon at drag epoch for the true (fitted) cosmology.
survey_volume: Total survey volume in Mpc^3.
tot_num_source_in_box: The total number of mock sources in the box needed to achieve self.num_discrete_source number of sources in the survey volume.
tracer_bias_1: The linear bias of the first tracer.
tracer_bias_2: The linear bias of the second tracer.
true_cosmology
unitless_1: Whether field_1 needs to be divided by its mean
unitless_2: Whether field_2 needs to be divided by its mean
vel_resol: velocity resolution on average in km/s
vel_resol_ch: velocity resolution of each channel in km/s
w0: The dark energy equation of state parameter for the true (fitted) cosmology.
w_HI: The weights per map pixel.
wa: The dark energy equation of state parameter for the true (fitted) cosmology.
weights_1: The weights of the first tracer in the rectangular grid.
weights_2: The weights of the second tracer in the rectangular grid.
weights_field_1: The weights of the first tracer in the density field.
weights_field_2: The weights of the second tracer in the density field.
weights_grid_1: The weights of the first tracer in the rectangular grid.
weights_grid_2: The weights of the second tracer in the rectangular grid.
weights_map_pixel: The weights per map pixel.
wproj: The WCS projection object for the map geometry.
x_mode: The mode of the 3D x-vector.
x_vec: The 3D x-vector of the box.
z: The effective centre redshift of the frequency range
z_as_func_of_comov_dist: Returns a function that returns the redshift for input comoving distance.
z_ch: The redshift of each frequency channel
z_gal: The redshifts of each galaxy in the catalogue for cross-correlation.
z_mock_tracer: The redshift of mock galaxies

Methods

`apply_taper_to_field`(field[, taper_func, axis])	Apply a taper to the field, by multiplying the taper function to the corresponding weights of the field.
`beam_attenuation`()	The beam attenuation factor.
`cal_rsd_power`(power_in_real_space, beta1, ...)	Calculate the redshift space power spectrum.
`check_is_map_noiselike_using_pca`(A_mat[, ...])	Use the source mixing matrix from eigendecomposition of the covariance matrix, project out the map data with more and more modes, and check if the variance of the residual map behaves like white noise.
`clean_cache`(attr)	Set the attributes to None.
`convolve_data`([kernel, data, weights, ...])	Convolve map data.
`create_white_noise_map`(sigma_N[, counts, ...])	Create a white noise map with the given standard deviation.
`deltav_to_deltar`(delta_v)	Convert a velocity interval delta_v to a comoving distance interval delta_r.
`deltaz_to_deltar`(delta_z)	Convert a redshift interval delta_z to a comoving distance interval delta_r.
`fog_gaussian`(sigma_r[, kmode, mumode])	The Gaussian finger-of-god profile.
`fog_lorentz`(sigma_r[, kmode, mumode])	The Lorentzian finger-of-god profile.
`fog_term`(sigma_r[, kmode, mumode])	The FoG term for the model calculation.
`gen_random_poisson_galaxy`([sel, num_g_rand, ...])	Generate a random galaxy catalogue from the map cube following the Poisson distribution.
`generate_full_healpix_map`([data, fill_value])	Generate a full healpix map from the data.
`get_1d_power`(power3d[, k1dbins, k1dweights, ...])	Bin the 3D power spectrum into 1D power spectrum.
`get_alpha_parallel`()	Calculate the line of sight Alcock–Paczynski effect parameter.
`get_alpha_perp`()	Calculate the transverse Alcock–Paczynski effect parameter.
`get_beam_image`([wproj, num_pix_x, ...])	Calculate the beam image projected onto the sky map for the input beam model.
`get_beam_window_ch`([ch_sel, cache, lmax, ...])	Build per-channel harmonic beam windows for HEALPix smoothing.
`get_cospar`(cosmology)	Generate a `CosmologyParameters` object from the input cosmology.
`get_counts_in_box`([partial_sel])	Grid the counts of the map cube voxels into the rectangular box, and return the effective counts per rectangular grid.
`get_cy_power`(power3d[, kperpbins, ...])	Bin the 3D power spectrum into cylindrical k_perp-k_para power spectrum.
`get_enclosing_box`([rot_mat])	invoke to calculate the box dimensions for enclosing all the map pixels.
`get_fourier_field_1`()	Calculate the Fourier transform of the density field of the first tracer.
`get_fourier_field_2`()	Calculate the Fourier transform of the density field of the second tracer.
`get_jackknife_patches`(ra_patch_num, ...[, ...])	Split the map into roughly equal patches.
`get_matter_power_spectrum`()	Calculate the matter power spectrum, interpolate it, and save it into the class attribute matter_power_spectrum_fnc.
`get_mock_field_from_power`(power_array)	Generate a mock field that follows the input matter power spectrum `power_array`.
`get_mock_field_in_redshift_space`(delta_x, ...)	Generate the mock field in redshift space.
`get_mock_field_r`([bias])	Generate a mock field in real space that follows the input matter power spectrum `self.matter_power_spectrum_fnc` * bias**2.
`get_mock_kaiser_field_k`(field)	Generate the Kaiser rsd effect for the mock matter field in k-space.
`get_mock_matter_field`()	Generate the mock matter field in redshift space.
`get_mock_matter_field_r`()	Generate the mock matter field in real space.
`get_mock_tracer_field`(tracer_i)	Generate the mock tracer field in redshift space.
`get_mock_tracer_field_r`(tracer_i)	Generate the mock tracer field in real space.
`get_mock_tracer_position_in_box`(tracer_i[, ...])	Function to retrieve the tracer positions in redshift space.
`get_mock_tracer_position_in_radecz`()	Project the mock tracer positions in the rectangular box to the sky coordinates and redshifts.
`get_mock_velocity_u_field`(mock_field)	Generate the normalised peculiar velocity field in real space
`get_model_matter_power`()	Calculate the model matter power spectrum.
`get_model_matter_power_r`()	Calculate the model matter power spectrum in real space (without RSD).
`get_model_power_cross`()	Calculate the model cross power spectrum between the two tracers.
`get_model_power_i`(i)	Calculate the model power spectrum for the i-th tracer.
`get_model_power_noobs_cross`()	Calculate the model cross power spectrum between the two tracers without observational effects.
`get_model_power_noobs_i`(i)	Calculate the model power spectrum for the i-th tracer without observational effects.
`get_n_bar_correction`()	Calculate the number density correction for the galaxy catalogue.
`get_sound_horizon_drag`(Omh2, Obh2, Onuh2)	Calculate the sound horizon at drag epoch for a set of cosmological parameters.
`get_weights_none_to_one`(attr_name)	Get the weights, and if it is None, convert it to 1.0 of size of kmode.
`get_z_as_func_of_comov_dist`()	Calculate an array of comoving distances with redshifts, and construct a function that returns the redshift for input comoving distance.
`grid_data_to_field`([flat_sky, partial_sel])	Grid the stored data map to a rectangular grid field.
`grid_field_to_sky_map`(field[, average, ...])	Grid a field in the rectangular box onto the sky.
`grid_gal_to_field`([radecfreq, flat_sky])	Grid the galaxy catalogue to a rectangular grid field.
`gridding_compensation`()	The sampling window function to be compensated for the gridding mass assignment scheme.
`map_sampling`([sampling_resol, p])	The sampling window function from the map cube to be convolved with model power spectrum.
`propagate_field_k_to_model`()	Use field k-modes for the model, taking into account the Alcock–Paczynski effect.
`propagate_mock_field_to_data`(field[, beam, ...])	Grid the mock tracer field onto the sky map cube.
`propagate_mock_tracer_to_gal_cat`([trim])	Propagate the mock tracer positions to the galaxy data catalogue.
`ra_dec_z_for_coord_in_box`(pos_in_box)	Convert the coordinates in the box to ra, dec, z, and also return the comoving distance to the observer for each point.
`read_from_fits`()	Read in a FITS file for the map data and hit counts.
`read_from_pickle`()	Read in a pickle file for cross-correlation and save the data into the class attributes.
`read_gal_cat`([ra_col, dec_col, z_col, trim])	Read in a galaxy catalogue for cross-correlation and save the data into the class attributes.
`set_corr_type`(corr_type, tracer_indx)	A utility function to help decide whether a tracer field needs to be mean centred, renormalised by its mean, and shot noise removed.
`set_radecnu_bounds_from_map`()	Set `ra_range`, `dec_range`, `nu_min`, and `nu_max` from the loaded `_ra_map`, `_dec_map`, and `nu` (tight RA interval, declination min/max, frequency min/max).
`trim_gal_to_range`()	Trim the galaxy catalogue to the specified range.
`trim_map_to_range`()	Trim the map to the specified range.
`use_flat_sky_box`([flat_sky_padding])	Use flat sky approximation to calculate the box dimensions.

property dec_mock_tracer: The Dec coordinate of mock galaxies on the sky

property discrete_base_field: Which tracer (1 or 2) field to sample the discrete tracer positions.

property discrete_source_dndz: The redshift kernel of the discrete tracer sources. Must be a function of redshift. ** Only the shape of the dndz is used, not the normalization. ** The overall number of discrete sources is given by self.num_discrete_source. Must be in the unit of per volume instead of unitless , i.e. the integral of the dndz over redshift must be a number density instead of number of sources.

get_mock_field_from_power(power_array)[source]: Generate a mock field that follows the input matter power spectrum power_array.

get_mock_field_in_redshift_space(delta_x, field, sigma_v)[source]

Generate the mock field in redshift space.

Parameters:

delta_x (np.ndarray) – The mock field in real space.
field (str) – The field to be simulated. Can be “matter” or “tracer_1” or “tracer_2”.
sigma_v (float) – The velocity dispersion of the tracer in km/s.

Returns:

delta_x – The mock field in redshift space.

Return type:

np.ndarray

get_mock_field_r(bias=1)[source]: Generate a mock field in real space that follows the input matter power spectrum self.matter_power_spectrum_fnc * bias**2.

get_mock_kaiser_field_k(field)[source]

Generate the Kaiser rsd effect for the mock matter field in k-space.

In the parrallel-plane limit, the Kaiser rsd effect is given by:

\[\delta_{\rm rsd} = f \mu^2 \delta_k\]

where \(f\) is the growth rate, \(\mu\) is the cosine of the angle between the wave vector and the line of sight, and \(\delta_k\) is the Fourier transform of the real-space matter density field.

This function returns \(\delta_{\rm rsd}\) in k-space so that for any mock tracer field, the Kaiser effect can be applied by adding \(\delta_{\rm rsd}\) to the real-space tracer field in k-space.

get_mock_matter_field()[source]: Generate the mock matter field in redshift space.

get_mock_matter_field_r()[source]: Generate the mock matter field in real space.

get_mock_tracer_field(tracer_i)[source]

Generate the mock tracer field in redshift space.

Parameters:: tracer_i (int) – The index of the tracer. Can be 1 or 2.
Returns:: delta_x – The mock tracer field in redshift space.
Return type:: np.ndarray

get_mock_tracer_field_r(tracer_i)[source]

Generate the mock tracer field in real space.

Parameters:: tracer_i (int) – The index of the tracer. Can be 1 or 2.
Returns:: delta_x – The mock tracer field in real space.
Return type:: np.ndarray

get_mock_tracer_position_in_box(tracer_i, density_field=None)[source]

Function to retrieve the tracer positions in redshift space. Modified from powerbox.

Parameters:

tracer_i (int) – The index of the tracer. Can be 1 or 2.
density_field (np.ndarray, default None) – The density field of the tracer. If None, the simulated mock tracer field will be used.

Returns:

tracer_positions – The tracer positions in the rectangular box.

Return type:

np.ndarray

get_mock_tracer_position_in_radecz()[source]: Project the mock tracer positions in the rectangular box to the sky coordinates and redshifts.

get_mock_velocity_u_field(mock_field)[source]

Generate the normalised peculiar velocity field in real space

Parameters:: mock_field (np.ndarray) – The mock field in real space.
Returns:: u_r – The normalised peculiar velocity field in real space.
Return type:: np.ndarray

property highres_sim

Optional integer upsampling factor for HIGalaxySimulation.propagate_hi_profile_to_map().

For WCS, propagate_mock_field_to_data() ignores this value and emits a DeprecationWarning when it is not None.

property mock_amp_1: The overall amplitude of the mock tracer field 1 in each grid. This can be T_HI(z) at each grid for example, which is multiplied to the simulated density field to get the final mock field 1 self.mock_tracer_field_1. If not given, default is to use self.mean_amp_1 as the amplitude.

property mock_amp_2: The overall amplitude of the mock tracer field 2 in each grid. This can be T_HI(z) at each grid for example, which is multiplied to the simulated density field to get the final mock field 2 self.mock_tracer_field_2. If not given, default is to use self.mean_amp_2 as the amplitude.

property mock_inside_range: Whether the mock galaxies are inside the survey area and frequency range

property mock_kaiser_field_k_matter: The Kaiser rsd effect correction for the mock matter field in k-space.

property mock_kaiser_field_k_tracer_1

The Kaiser rsd effect correction for the mock tracer field 1 in k-space.

See the docstring of get_mock_kaiser_field_k() for more details.

property mock_kaiser_field_k_tracer_2

The Kaiser rsd effect correction for the mock tracer field 2 in k-space.

See the docstring of get_mock_kaiser_field_k() for more details.

property mock_matter_field: The simulated dark matter density field in redshift space. If self.kaiser_rsd is False, it is simply set to the real space field.

property mock_matter_field_r: The simulated dark matter density field in real space.

property mock_tracer_field_1: The simulated tracer field 1 in redshift space with unit if mock_amp_1 or mean_amp_1 is given.

property mock_tracer_field_1_r: The simulated tracer field 1 unitsless density contrast in real space.

property mock_tracer_field_2: The simulated tracer field 2 in redshift space with unit if mock_amp_2 or mean_amp_2 is given.

property mock_tracer_field_2_r: The simulated tracer field 2 unitsless density contrast in real space.

property mock_tracer_position_in_box: The simulated tracer positions in the box in real space.

property mock_tracer_position_in_radecz

The simulated tracer positions projected onto the grid. The tracers outside the binary selection window map_has_sampling or outside the frequency range can be trimmed off. Furthermore, excess tracers are also excluded. The selection function is stored at inside_range.

Returns a list in the order of (ra,dec,z,inside_range).

property mock_velocity_u_matter

The normalised peculiar velocity field in real space, defined as

\[u_i = -\frac{v_i}{\mathcal{H} f}\]

where \(v_i\) is the peculiar velocity, \(\mathcal{H}\) is the conformal Hubble parameter, and \(f\) is the growth rate.

property mock_velocity_u_tracer_1

The normalised peculiar velocity field used for the first tracer. While the peculiar velocity field is only dependent on the matter field, note that if you simulate lognormal tracer fields, the two fields are not exactly correlated. One simple way to think about it is that field_1/bias_1 is not equal to field_2/bias_2, because both fields range from -1 to +max, and similarly the mock matter field is also not field_1/bias_1 or field_2/bias_2. This is an intrinsic weakness of the lognormal mock that you should be aware of.

As a result, for each tracer the velocity field should be recalculated based on field_i/tracer_bias_i.

property mock_velocity_u_tracer_2

The normalised peculiar velocity field used for the second tracer.

See the docstring of mock_velocity_u_tracer_1() for more details.

property num_discrete_source: The total number of discrete tracer sources. Note that the final mock catalogue is not exactly this number, due to Poisson sampling errors.

property parallel_plane: Whether the mock field is generated in the parallel-plane limit.

propagate_mock_field_to_data(field, beam=True, highres_sim=None, average=True)[source]

Grid the mock tracer field onto the sky map cube.

Parameters:

field (np.ndarray) – The mock tracer field to be gridded.
beam (bool, default True) – Whether to convolve the beam model to the sky map.
highres_sim (int, default None) – Deprecated for WCS. If this argument or self.highres_sim is not None while the sky map uses WCS, a DeprecationWarning is emitted and the factor is ignored; gridding always uses the native map resolution.
average (bool, default True) – If True, then the map pixel value will be the average of the rectangular box values (i.e. the mock field is in temperature units). If False, then the map pixel value will be the sum of the rectangular box values (i.e. the mock field is in flux density units).

Returns:

The output sky map cube.

Return type:

np.ndarray

propagate_mock_tracer_to_gal_cat(trim=True)[source]

Propagate the mock tracer positions to the galaxy data catalogue.

If trim, only tracers inside the ra-dec-z range will be propagated.

Parameters:: trim (bool, default True) – If True, only the mock tracer positions inside the ra-dec-z range will be propagated.

property ra_mock_tracer: The RA coordinate of mock galaxies on the sky

property rsd_from_field: If True, the kaiser rsd effect is generated at field level by calculating the corresponding peculiar velocity field. This allows the lognormal mock to go beyond the parallel-plane limit. If False, the kaiser rsd effect is generated at power spectrum level assuming parallel-plane, and then the field is generated from the anistropic power spectrum.

property tot_num_source_in_box: The total number of mock sources in the box needed to achieve self.num_discrete_source number of sources in the survey volume. Only used internally, no physical meaning. Note that if you change the simulation settings such as self.num_discrete_source, self.discrete_source_dndz, self.z_ch, self.W_HI, etc, this property will be automatically updated but the mock catalog is not.

property z_mock_tracer: The redshift of mock galaxies

meer21cm.mock.generate_colored_noise(x_size, x_len, power_spectrum, seed=None)[source]

Generate random 1D gaussian fluctuations following a specific spectrum. This is similar to colorednoise package for generating colored noise. It is simply wrapping the generate_gaussian_field function under the hood. Note that the Fourier convention used should be consistent with np.fft, and the power spectrum is dimensionless.

Parameters:

x_size (int) – The number of sampling.
x_len (float) – The total length of the sampling.
power_spectrum (array) – The power spectrum of the random noise in Fourier space.
seed (int, default None) – The seed number for random generator for sampling. If None, a random seed is used.

Returns:

rand_arr – The random noise.

Return type:

float array.

meer21cm.mock.generate_gaussian_field(box_ndim, box_len, power_spectrum, seed, ps_has_volume=True)[source]

Generate a Gaussian field with the given power spectrum. If ps_has_volume is True, the power spectrum is assumed to have the volume unit, usually in the unit of \((Mpc)^{-3}\).

The length unit of the box is arbitrary, as long as the unit in box_len and power_spectrum are consistent.

Parameters:

box_ndim (array) – The number of grid points in each dimension.
box_len (array) – The length of the box in each dimension.
power_spectrum (array) – The input power spectrum.
seed (int) – The seed for the random number generator.
ps_has_volume (bool, default True) – If True, the power spectrum is assumed to have the volume unit.

Returns:

delta_x – The generated Gaussian field.

Return type:

array

meer21cm.mock.generate_lognormal_field(box_ndim, box_len, power_spectrum, seed, ps_has_volume=True)[source]

Generate a lognormal field with the given power spectrum. If ps_has_volume is True, the power spectrum is assumed to have the volume unit, usually in the unit of \((Mpc)^{-3}\).

The length unit of the box is arbitrary, as long as the unit in box_len and power_spectrum are consistent.

Parameters:

box_ndim (array) – The number of grid points in each dimension.
box_len (array) – The length of the box in each dimension.
power_spectrum (array) – The input power spectrum.
seed (int) – The seed for the random number generator.
ps_has_volume (bool, default True) – If True, the power spectrum is assumed to have the volume unit.

Returns:

delta_x – The generated lognormal field.

Return type:

array

meer21cm.mock.hi_mass_to_flux_profile(loghimass, z_g, nu, tf_slope=None, tf_zero=None, cosmo=FlatLambdaCDM(name='Planck18', H0=<Quantity 67.66 km / (Mpc s)>, Om0=0.30966, Tcmb0=<Quantity 2.7255 K>, Neff=3.046, m_nu=<Quantity [0., 0., 0.06] eV>, Ob0=0.04897), seed=None, num_ch_ext_on_each_side=5, internal_step=1001, no_vel=True)[source]

Convert HI mass to emission line profile.

The relation between mass and flux (flux density integrated) of a HI source can be written as [1]

\[M_{\rm HI} = \frac{16\pi m_H}{3 h f_{21} A_{10}}\, D_L^2\, S,\]

where \(M_{\rm HI}\) is the HI mass, \(m_H\) is the mass of neutral hydrogen atom, \(h\) is the Planck constant, \(f_{21}\) is the rest frequency of 21cm line, \(A_{10}\) is the spontaneous emission rate of 21cm line, \(D_L\) is the luminosity distance of the source, and \(S\) is the flux.

If no_vel is set to False, the w50 parameters of the HI sources are calculated using a Tully-Fisher relation. Random inclinations and busy function parameters are assigned to the sources. The busy functions are then used as the emission profiles. The profiles are gridded into frequency channels. The gridded profile, hifluxd_ch has the shape of (ch_offset,num_source), with the zeroth axis corresponding to the channel offset, and varies from (-N,-N+1,…,0,…,N-1,N).

Parameters:

loghimass (array) – The HI mass of sources in log10 solmar mass (no h).
z_g (array) – The redshifts of the sources
nu (array) – The frequency channels in Hz.
tf_slope (float, default None.) – The slope of the T-F relation. See meer21cm.util.tully_fisher().
tf_zero (float, default None.) – The intercept of the T-F relation. See meer21cm.util.tully_fisher().
cosmo (optional, default Planck18.) – The cosmology used.
seed (optional, default None.) – The seed number for rng.
num_ch_ext_on_each_side (optional, default 5.) – Internal parameter, no need to change.
internal_step (optional, default 1001.) – Internal parameter, decrease for less accuracy. Can be lower when velocity resolution is low.
no_vel (optional, default True.) – If True, source will have no emission line profile and just a delta peak.

Returns:

hifluxd_ch – The flux density of each source in Jy.

Return type:

array

References

[1]

Meyer et al., “Tracing HI Beyond the Local Universe”, https://arxiv.org/abs/1705.04210

meer21cm.power module

This module handles computation of power spectrum from gridded fields and its corresponding model power spectrum from theory.

The class ModelPowerSpectrum is the class for computing the model power spectrum of an LSS tracer field.

The class FieldPowerSpectrum is the class for computing the power spectrum of a gridded field from LSS data.

The class PowerSpectrum coherently combines the two classes above, and provides an interface for gridding the intensity mapping data as well as the galaxy catalogue to perform power spectrum estimation and for auto-correlation and cross-correlation.

class meer21cm.power.FieldPowerSpectrum(field_1, box_len, weights_1=None, mean_center_1=False, unitless_1=False, field_2=None, weights_2=None, mean_center_2=False, unitless_2=False, **params)[source]

Bases: Specification

The class for computing the power spectrum of a gridded field from LSS data.

Parameters:

field_1 (np.ndarray.) – The density field of the first tracer.
field_2 (np.ndarray, default None) – The density field of the second tracer. If None, calculation of the second tracer and the cross-correlation will be skipped.
box_len (list of 3 floats.) – The length of the box along each axis.
weights_1 (np.ndarray, default None) – The weights of the first tracer. Default is uniform weights.
mean_center_1 (bool, default False) – Whether to mean-center the first tracer field.
unitless_1 (bool, default False) – Whether to divide the first tracer field by its mean.
weights_2 (np.ndarray, default None) – The weights of the second tracer. Default is uniform weights.
mean_center_2 (bool, default False) – Whether to mean-center the second tracer field.
unitless_2 (bool, default False) – Whether to divide the second tracer field by its mean.
**params (dict) – Additional parameters to be passed to the base class meer21cm.dataanalysis.Specification.

Attributes:

W_HI: A binary window for whether a pixel has been sampled
auto_power_3d_1: The 3D power spectrum of the first tracer.
auto_power_3d_2: The 3D power spectrum of the second tracer.
batch_number: Number of sequential batches used by gridding routines.
beam_image: Returns the beam image projected onto the sky map for the input beam model.
beam_model: The name of the beam function.
beam_type: The beam type that can be either be isotropic or anisotropic.
beam_unit: The unit of input beam size parameter sigma
beam_window_ch: Per-channel spherical-harmonic beam window \(B_\ell\) for HEALPix hp.smoothing workflows.
box_len: The length of all sides of the box in Mpc.
box_ndim: The number of grids along each side of the enclosing box.
box_resol: The grid length of each side of the enclosing box in Mpc.
ch_id_gal: The channel id (0-indexed) of each galaxy in the catalogue for cross-correlation.
counts: The number of hits per pixel for the map data
cross_power_3d: The 3D cross power spectrum between the two tracers.
data: The map data
dec_gal: The declination of each galaxy in the catalogue for cross-correlation.
dec_map: The declination of each pixel in the map.
dvdf: velocity resolution per unit frequency on average, in km/s/Hz
dvdf_ch: velocity resolution per unit frequency in each channel, in km/s/Hz
field_1: The density field of the first tracer.
field_2: The density field of the second tracer.
fourier_field_1: The Fourier transform of the density field of the first tracer.
fourier_field_2: The Fourier transform of the density field of the second tracer.
freq_gal: The 21cm line frequency for each galaxy in Hz.
freq_resol: frequency resolution in Hz
hp_nside: HEALPix \(N_{side}\) for HEALPix-backed specifications.
k_mode: The fiducial (observed) mode of the 3D k-vector.
k_nyquist: The Nyquist frequency of the 3D box along each axis.
k_para: The fiducial parallel k-mode of the 3D box.
k_perp: The fiducial perpendicular k-vector of the 3D box.
k_vec: The 3D k-vector of the box.
map_has_sampling: A binary window for whether a pixel has been sampled
map_unit_type: The type of the map unit.
maximum_sampling_channel: Returns the index of the frequency channel with the maximum sampling on the sky map.
mean_center_1: Whether field_1 needs to be mean centered
mean_center_2: Whether field_2 needs to be mean centered
mu_mode: The fiducial (observed) mu values of each k-mode so that \(k_\parallel = k imes \mu\).
nu: The input frequencies of the survey
num_pix_x: The number of pixels along the first map axis.
num_pix_y: The number of pixels along the second map axis.
pix_resol: angular resolution of the map pixel in deg
pixel_area: angular area of the map pixel in deg^2
pixel_id: Sparse HEALPix pixel indices (RING, nest=False).
precision: Floating precision flag.
ra_gal: The right ascension of each galaxy in the catalogue for cross-correlation.
ra_map: The right ascension of each pixel in the map.
real_dtype: Active real floating dtype controlled by precision.
renorm_ps_1: The renormalization factor of the power spectrum of the first tracer.
renorm_ps_2: The renormalization factor of the power spectrum of the second tracer.
renorm_ps_cross: The renormalization factor of the cross power spectrum.
sigma_beam_ch: The input beam size parameter sigma for each channel.
sigma_beam_in_mpc: The channel averaged beam size in Mpc
unitless_1: Whether field_1 needs to be divided by its mean
unitless_2: Whether field_2 needs to be divided by its mean
vel_resol: velocity resolution on average in km/s
vel_resol_ch: velocity resolution of each channel in km/s
w_HI: The weights per map pixel.
weights_map_pixel: The weights per map pixel.
wproj: The WCS projection object for the map geometry.
x_mode: The mode of the 3D x-vector.
x_vec: The 3D x-vector of the box.
z: The effective centre redshift of the frequency range
z_ch: The redshift of each frequency channel
z_gal: The redshifts of each galaxy in the catalogue for cross-correlation.

Methods

`check_is_map_noiselike_using_pca`(A_mat[, ...])	Use the source mixing matrix from eigendecomposition of the covariance matrix, project out the map data with more and more modes, and check if the variance of the residual map behaves like white noise.
`clean_cache`(attr)	Set the attributes to None.
`convolve_data`([kernel, data, weights, ...])	Convolve map data.
`create_white_noise_map`(sigma_N[, counts, ...])	Create a white noise map with the given standard deviation.
`generate_full_healpix_map`([data, fill_value])	Generate a full healpix map from the data.
`get_beam_image`([wproj, num_pix_x, ...])	Calculate the beam image projected onto the sky map for the input beam model.
`get_beam_window_ch`([ch_sel, cache, lmax, ...])	Build per-channel harmonic beam windows for HEALPix smoothing.
`get_fourier_field_1`()	Calculate the Fourier transform of the density field of the first tracer.
`get_fourier_field_2`()	Calculate the Fourier transform of the density field of the second tracer.
`get_jackknife_patches`(ra_patch_num, ...[, ...])	Split the map into roughly equal patches.
`get_weights_none_to_one`(attr_name)	Get the weights, and if it is None, convert it to 1.0 of size of kmode.
`read_from_fits`()	Read in a FITS file for the map data and hit counts.
`read_from_pickle`()	Read in a pickle file for cross-correlation and save the data into the class attributes.
`read_gal_cat`([ra_col, dec_col, z_col, trim])	Read in a galaxy catalogue for cross-correlation and save the data into the class attributes.
`set_corr_type`(corr_type, tracer_indx)	A utility function to help decide whether a tracer field needs to be mean centred, renormalised by its mean, and shot noise removed.
`set_radecnu_bounds_from_map`()	Set `ra_range`, `dec_range`, `nu_min`, and `nu_max` from the loaded `_ra_map`, `_dec_map`, and `nu` (tight RA interval, declination min/max, frequency min/max).
`trim_gal_to_range`()	Trim the galaxy catalogue to the specified range.
`trim_map_to_range`()	Trim the map to the specified range.

property auto_power_3d_1: The 3D power spectrum of the first tracer.

property auto_power_3d_2: The 3D power spectrum of the second tracer.

property box_len: The length of all sides of the box in Mpc.

property box_ndim: The number of grids along each side of the enclosing box. To ensure even sampling of +k and -k modes, the number of grids along every axis needs to be odd.

property box_resol: The grid length of each side of the enclosing box in Mpc.

property cross_power_3d: The 3D cross power spectrum between the two tracers.

property field_1: The density field of the first tracer.

property field_2: The density field of the second tracer.

property fourier_field_1: The Fourier transform of the density field of the first tracer.

property fourier_field_2: The Fourier transform of the density field of the second tracer.

get_fourier_field_1()[source]: Calculate the Fourier transform of the density field of the first tracer.

get_fourier_field_2()[source]: Calculate the Fourier transform of the density field of the second tracer.

property k_mode: The fiducial (observed) mode of the 3D k-vector.

property k_nyquist: The Nyquist frequency of the 3D box along each axis.

property k_para: The fiducial parallel k-mode of the 3D box.

property k_perp: The fiducial perpendicular k-vector of the 3D box.

property k_vec: The 3D k-vector of the box.

property mean_center_1: Whether field_1 needs to be mean centered

property mean_center_2: Whether field_2 needs to be mean centered

property mu_mode: The fiducial (observed) mu values of each k-mode so that \(k_\parallel = k imes \mu\).

property renorm_ps_1: The renormalization factor of the power spectrum of the first tracer.

property renorm_ps_2: The renormalization factor of the power spectrum of the second tracer.

property renorm_ps_cross: The renormalization factor of the cross power spectrum.

set_corr_type(corr_type, tracer_indx)[source]

A utility function to help decide whether a tracer field needs to be mean centred, renormalised by its mean, and shot noise removed. Currently only two types are supported, “Gal” and “HI” (case-insensitive). If the tracer is galaxy (number counts), the auto power spectrum is mean centred, renormalised, and then shot noise removed. If HI, none of the above will be performed.

Parameters:

corr_type (str) – The tracer type.
tracer_indx (int) – Either 1 or 2.

property unitless_1: Whether field_1 needs to be divided by its mean

property unitless_2: Whether field_2 needs to be divided by its mean

property x_mode: The mode of the 3D x-vector.

property x_vec: The 3D x-vector of the box.

class meer21cm.power.ModelPowerSpectrum(kmode=None, mumode=None, tracer_bias_1=1.0, sigma_v_1=0.0, tracer_bias_2=None, sigma_v_2=0.0, include_beam=[True, False], fog_profile='lorentz', cross_coeff=1.0, weights_field_1=None, weights_field_2=None, weights_grid_1=None, weights_grid_2=None, mean_amp_1=1.0, mean_amp_2=1.0, sampling_resol=None, include_sky_sampling=[True, False], compensate=[True, True], kaiser_rsd=True, sigma_z_1=0.0, sigma_z_2=0.0, **params)[source]

Bases: CosmologyCalculator

The class for computing the model power spectrum of an LSS tracer field.

Parameters:

kmode (np.ndarray, default None) – The true mode of k in Mpc-1.
mumode (np.ndarray, default None) – The true mu values of each k-mode so that \(k_\parallel = k \times \mu\).
tracer_bias_1 (float, default 1.0) – The linear bias of the first tracer.
sigma_v_1 (float, default 0.0) – The velocity dispersion of the first tracer in km/s.
tracer_bias_2 (float, default None) – The linear bias of the second tracer.
sigma_v_2 (float, default 0.0) – The velocity dispersion of the second tracer in km/s.
include_beam (list, default [True, False]) – Whether to include the beam attenuation in the model calculation. Must be a list of two booleans, the first for the first tracer and the second for the second tracer.
fog_profile (str, default "lorentz") – The shape of the finger-of-god profile to be used in the model calculation. Either “lorentz” or “gaussian”.
cross_coeff (float, default 1.0) – The cross-correlation coefficient between the two tracers.
weights_field_1 (np.ndarray, default None) – The field-level weights of the first tracer in the density field.
weights_field_2 (np.ndarray, default None) – The field-level weights of the second tracer in the density field.
weights_grid_1 (np.ndarray, default None) – The grid-level weights of the first tracer in the density field.
weights_grid_2 (np.ndarray, default None) – The grid-level weights of the second tracer in the density field.
mean_amp_1 (float, default 1.0) – The mean amplitude of the first tracer. Can be used to rescale the power spectrum, for example by the average brightness temperature.
mean_amp_2 (float, default 1.0) – The mean amplitude of the second tracer. Can be used to rescale the power spectrum, for example by the average brightness temperature.
sampling_resol (list, default None) – The sampling resolution of the field in Mpc. If sampling_resol is “auto”, the sampling resolution will be set to the pixel size of the map.
include_sky_sampling (list, default [True, False]) – Whether to include the sky sampling in the model calculation. If just a boolean is provided, it will be used for both tracers.
compensate (list, default [True, True]) – Whether the gridded fields are compensated according to the mass assignment scheme. Note that the compensation is applied to the model power spectrum, and not to the gridded data fields.
kaiser_rsd (bool, default True) – Whether to include the RSD effect in the model calculation and mock simulation.
sigma_z_1 (float, default 0.0) – The redshift error of the first tracer.
sigma_z_2 (float, default 0.0) – The redshift error of the second tracer.
**params (dict) – Additional parameters to be passed to the base class meer21cm.cosmology.CosmologyCalculator.

Attributes:

As: The amplitude of the primordial power spectrum for the true (fitted) cosmology.
W_HI: A binary window for whether a pixel has been sampled
alpha_AP: The anisotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel} / \alpha_{\perp}\).
alpha_iso: The isotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel}^{1/3} \alpha_{\perp}^{2/3}\).
alpha_parallel: The line of sight Alcock–Paczynski effect parameter.
alpha_perp
auto_power_matter_model: The model matter power spectrum with RSD effects.
auto_power_matter_model_r: The model matter power spectrum in real space (without RSD).
auto_power_tracer_1_model: The 3D model power spectrum for the first tracer.
auto_power_tracer_1_model_noobs: The model power spectrum for the first tracer without observational effects.
auto_power_tracer_2_model: The 3D model power spectrum for the second tracer.
auto_power_tracer_2_model_noobs: The model power spectrum for the second tracer without observational effects.
average_hi_temp: The average HI brightness temperature in Kelvin at central redshift self.z.
backend: Which backend to use for computing the matter power.
batch_number: Number of sequential batches used by gridding routines.
beam_image: Returns the beam image projected onto the sky map for the input beam model.
beam_model: The name of the beam function.
beam_type: The beam type that can be either be isotropic or anisotropic.
beam_unit: The unit of input beam size parameter sigma
beam_window_ch: Per-channel spherical-harmonic beam window \(B_\ell\) for HEALPix hp.smoothing workflows.
ch_id_gal: The channel id (0-indexed) of each galaxy in the catalogue for cross-correlation.
cold: If True (recommended), the matter power spectrum is the CDM ps.
compensate: Whether the gridded fields are compensated according to the mass assignment scheme.
cosmo: A shortcut to the astropy.cosmology.Cosmology object for the true cosmology.
cospar_fiducial
cospar_true
counts: The number of hits per pixel for the map data
cross_coeff: The cross-correlation coefficient between the two tracers.
cross_power_tracer_model: The 3D model cross power spectrum between the two tracers.
cross_power_tracer_model_noobs: The model power spectrum for the cross correlation between the two tracers without observational effects.
data: The map data
dec_gal: The declination of each galaxy in the catalogue for cross-correlation.
dec_map: The declination of each pixel in the map.
dvdf: velocity resolution per unit frequency on average, in km/s/Hz
dvdf_ch: velocity resolution per unit frequency in each channel, in km/s/Hz
expfactor: The expansion factor
f_growth_fiducial: The growth factor at self.expfactor for the fiducial cosmology.
f_growth_true: The growth factor at self.expfactor for the true (fitted) cosmology.
fiducial_cosmology
fog_profile: The shape of the finger-of-god profile to be used in the model calculation.
freq_gal: The 21cm line frequency for each galaxy in Hz.
freq_resol: frequency resolution in Hz
h: The Hubble constant for the true (fitted) cosmology.
hp_nside: HEALPix \(N_{side}\) for HEALPix-backed specifications.
include_beam: Whether the beam attenuation is included in the model calculation.
kaiser_rsd: Whether RSD is included in the simulation and model calculation.
kmax: The maximum k in Mpc^-1 for calculating matter power.
kmin: The minimum k in Mpc^-1 for calculating matter power.
kmode: The input kmode for the model calculation.
los_resol_in_mpc: effective frequency resolution in Mpc for the fiducial cosmology.
map_has_sampling: A binary window for whether a pixel has been sampled
map_unit_type: The type of the map unit.
matter_power_spectrum_fnc: Interpolation function for the real-space isotropic matter power spectrum.
maximum_sampling_channel: Returns the index of the frequency channel with the maximum sampling on the sky map.
mumode: The mu values of each 3D k-mode.
neutrino_mass: The neutrino mass for the true (fitted) cosmology.
ns: The primordial power spectrum spectral index for the true (fitted) cosmology.
nu: The input frequencies of the survey
num_kpoints: The number of k points for calculating matter power.
num_pix_x: The number of pixels along the first map axis.
num_pix_y: The number of pixels along the second map axis.
omega_baryon: The baryon matter density parameter for the true (fitted) cosmology.
omega_cold: The cold matter density parameter for the true (fitted) cosmology.
omega_hi: The HI density as a function of redshift, over the critical density of the Universe at z=0, defined at each frequency channel so that self.omega_hi corresponds to self.z_ch.
omega_hi_z_func: Interpolate the input self.omega_hi at each frequency channel self.z_ch to a function of redshift.
omega_hi_z_mean: The mean HI density at the central redshift self.z, over the critical density of the Universe at z=0.
pix_resol: angular resolution of the map pixel in deg
pix_resol_in_mpc: angular resolution of the map pixel in Mpc for the fiducial cosmology.
pixel_area: angular area of the map pixel in deg^2
pixel_id: Sparse HEALPix pixel indices (RING, nest=False).
precision: Floating precision flag.
ps_type: linear or nonlinear for the matter power.
ra_gal: The right ascension of each galaxy in the catalogue for cross-correlation.
ra_map: The right ascension of each pixel in the map.
real_dtype: Active real floating dtype controlled by precision.
sampling_resol: The sampling resolution corresponding to the map-making/gridding of the density field.
sigma_beam_ch: The input beam size parameter sigma for each channel.
sigma_beam_ch_in_mpc: The input beam size parameter in Mpc in each channel.
sigma_beam_in_mpc: The channel averaged beam size in Mpc
sigma_v_1: The velocity dispersion of the first tracer in km/s.
sigma_v_2: The velocity dispersion of the second tracer in km/s.
sigma_z_1: The redshift error of the first tracer.
sigma_z_2: The redshift error of the second tracer.
sound_horizon_drag_fiducial: The sound horizon at drag epoch for the fiducial cosmology.
sound_horizon_drag_true: The sound horizon at drag epoch for the true (fitted) cosmology.
survey_volume: Total survey volume in Mpc^3.
tracer_bias_1: The linear bias of the first tracer.
tracer_bias_2: The linear bias of the second tracer.
true_cosmology
vel_resol: velocity resolution on average in km/s
vel_resol_ch: velocity resolution of each channel in km/s
w0: The dark energy equation of state parameter for the true (fitted) cosmology.
w_HI: The weights per map pixel.
wa: The dark energy equation of state parameter for the true (fitted) cosmology.
weights_1: The weights of the first tracer in the rectangular grid.
weights_2: The weights of the second tracer in the rectangular grid.
weights_field_1: The weights of the first tracer in the density field.
weights_field_2: The weights of the second tracer in the density field.
weights_grid_1: The weights of the first tracer in the rectangular grid.
weights_grid_2: The weights of the second tracer in the rectangular grid.
weights_map_pixel: The weights per map pixel.
wproj: The WCS projection object for the map geometry.
z: The effective centre redshift of the frequency range
z_as_func_of_comov_dist: Returns a function that returns the redshift for input comoving distance.
z_ch: The redshift of each frequency channel
z_gal: The redshifts of each galaxy in the catalogue for cross-correlation.

Methods

`beam_attenuation`()	The beam attenuation factor.
`cal_rsd_power`(power_in_real_space, beta1, ...)	Calculate the redshift space power spectrum.
`check_is_map_noiselike_using_pca`(A_mat[, ...])	Use the source mixing matrix from eigendecomposition of the covariance matrix, project out the map data with more and more modes, and check if the variance of the residual map behaves like white noise.
`clean_cache`(attr)	Set the attributes to None.
`convolve_data`([kernel, data, weights, ...])	Convolve map data.
`create_white_noise_map`(sigma_N[, counts, ...])	Create a white noise map with the given standard deviation.
`deltav_to_deltar`(delta_v)	Convert a velocity interval delta_v to a comoving distance interval delta_r.
`deltaz_to_deltar`(delta_z)	Convert a redshift interval delta_z to a comoving distance interval delta_r.
`fog_gaussian`(sigma_r[, kmode, mumode])	The Gaussian finger-of-god profile.
`fog_lorentz`(sigma_r[, kmode, mumode])	The Lorentzian finger-of-god profile.
`fog_term`(sigma_r[, kmode, mumode])	The FoG term for the model calculation.
`generate_full_healpix_map`([data, fill_value])	Generate a full healpix map from the data.
`get_alpha_parallel`()	Calculate the line of sight Alcock–Paczynski effect parameter.
`get_alpha_perp`()	Calculate the transverse Alcock–Paczynski effect parameter.
`get_beam_image`([wproj, num_pix_x, ...])	Calculate the beam image projected onto the sky map for the input beam model.
`get_beam_window_ch`([ch_sel, cache, lmax, ...])	Build per-channel harmonic beam windows for HEALPix smoothing.
`get_cospar`(cosmology)	Generate a `CosmologyParameters` object from the input cosmology.
`get_jackknife_patches`(ra_patch_num, ...[, ...])	Split the map into roughly equal patches.
`get_matter_power_spectrum`()	Calculate the matter power spectrum, interpolate it, and save it into the class attribute matter_power_spectrum_fnc.
`get_model_matter_power`()	Calculate the model matter power spectrum.
`get_model_matter_power_r`()	Calculate the model matter power spectrum in real space (without RSD).
`get_model_power_cross`()	Calculate the model cross power spectrum between the two tracers.
`get_model_power_i`(i)	Calculate the model power spectrum for the i-th tracer.
`get_model_power_noobs_cross`()	Calculate the model cross power spectrum between the two tracers without observational effects.
`get_model_power_noobs_i`(i)	Calculate the model power spectrum for the i-th tracer without observational effects.
`get_sound_horizon_drag`(Omh2, Obh2, Onuh2)	Calculate the sound horizon at drag epoch for a set of cosmological parameters.
`get_weights_none_to_one`(attr_name)	Get the weights, and if it is None, convert it to 1.0 of size of kmode.
`get_z_as_func_of_comov_dist`()	Calculate an array of comoving distances with redshifts, and construct a function that returns the redshift for input comoving distance.
`gridding_compensation`()	The sampling window function to be compensated for the gridding mass assignment scheme.
`map_sampling`()	The sampling window function from the map cube to be convolved with data.
`read_from_fits`()	Read in a FITS file for the map data and hit counts.
`read_from_pickle`()	Read in a pickle file for cross-correlation and save the data into the class attributes.
`read_gal_cat`([ra_col, dec_col, z_col, trim])	Read in a galaxy catalogue for cross-correlation and save the data into the class attributes.
`set_radecnu_bounds_from_map`()	Set `ra_range`, `dec_range`, `nu_min`, and `nu_max` from the loaded `_ra_map`, `_dec_map`, and `nu` (tight RA interval, declination min/max, frequency min/max).
`trim_gal_to_range`()	Trim the galaxy catalogue to the specified range.
`trim_map_to_range`()	Trim the map to the specified range.

property auto_power_matter_model: The model matter power spectrum with RSD effects. The 3D k-modes corrospond to the input kmode and mumode.

property auto_power_matter_model_r: The model matter power spectrum in real space (without RSD).

property auto_power_tracer_1_model: The 3D model power spectrum for the first tracer. The 3D k-modes corrospond to the input kmode and mumode. Unlike noobs power, the mean amplitude is applied.

property auto_power_tracer_1_model_noobs: The model power spectrum for the first tracer without observational effects. Note that the power is in units of volume, so the mean amplitude is not applied.

property auto_power_tracer_2_model: The 3D model power spectrum for the second tracer. The 3D k-modes corrospond to the input kmode and mumode. Unlike noobs power, the mean amplitude is applied.

property auto_power_tracer_2_model_noobs: The model power spectrum for the second tracer without observational effects. Note that the power is in units of volume, so the mean amplitude is not applied.

beam_attenuation()[source]: The beam attenuation factor.

cal_rsd_power(power_in_real_space, beta1, sigma_v_1, sigma_z_1, beta2=None, sigma_v_2=None, sigma_z_2=None, r=1.0, mumode=None)[source]

Calculate the redshift space power spectrum. If properties of the second tracer are not set, they will be set to the same as the first tracer so the result will be the auto power spectrum.

Parameters:

power_in_real_space (np.ndarray) – The power spectrum in real space.
beta1 (float) – The growth rate over bias of the first tracer.
sigma_v_1 (float) – The velocity dispersion of the first tracer.
sigma_z_1 (float) – The redshift dispersion of the first tracer.
beta2 (float, default None) – The growth rate over bias of the second tracer.
sigma_v_2 (float, default None) – The velocity dispersion of the second tracer.
sigma_z_2 (float, default None) – The redshift dispersion of the second tracer.
r (float, default 1.0) – The correlation coefficient between the two tracers.
mumode (np.ndarray, default None) – The mu values of each 3D k-mode.

Returns:

power_in_redshift_space – The power spectrum in redshift space.

Return type:

np.ndarray

property compensate: Whether the gridded fields are compensated according to the mass assignment scheme. Note that the compensation is applied to the model power spectrum, and not to the gridded data fields.

property cross_coeff: The cross-correlation coefficient between the two tracers.

property cross_power_tracer_model: The 3D model cross power spectrum between the two tracers. The 3D k-modes corrospond to the input kmode and mumode. Unlike noobs power, the mean amplitude is applied.

property cross_power_tracer_model_noobs: The model power spectrum for the cross correlation between the two tracers without observational effects. Note that the power is in units of volume, so the mean amplitude is not applied.

fog_gaussian(sigma_r, kmode=None, mumode=None)[source]

The Gaussian finger-of-god profile.

\[{\rm FoG} = {\rm exp}(-(\sigma_r k_\parallel/H)^2/2)\]

Note the power spectrum has FoG squared with the two FoG terms that can be different for two tracers.

Parameters:

sigma_r (float.) – The velocity dispersion in terms of the comoving distance in Mpc.
kmode (float, None.) – The mode of 3D k in Mpc-1. If None, self.kmode will be used.
mumode (float, None.) – The mu values of each 3D k-mode. In None, self.mumode will be used.

Returns:

fog – The FoG term.

Return type:

float.

fog_lorentz(sigma_r, kmode=None, mumode=None)[source]

The Lorentzian finger-of-god profile.

\[{\rm FoG} = \sqrt{1/(1+(\sigma_r k_\parallel/H)^2)}\]

Note the power spectrum has FoG squared with the two FoG terms that can be different for two tracers.

Parameters:

sigma_r (float.) – The velocity dispersion in terms of the comoving distance in Mpc.
kmode (float, None.) – The mode of 3D k in Mpc-1. If None, self.kmode will be used.
mumode (float, None.) – The mu values of each 3D k-mode. In None, self.mumode will be used.

Returns:

fog – The FoG term.

Return type:

float.

property fog_profile: The shape of the finger-of-god profile to be used in the model calculation. Either “lorentz” or “gaussian”.

fog_term(sigma_r, kmode=None, mumode=None)[source]

The FoG term for the model calculation. It reads the profile type from the attribute fog_profile.

Parameters:

sigma_r (float.) – The velocity dispersion in terms of the comoving distance in Mpc.
kmode (float, None.) – The mode of 3D k in Mpc-1. If None, self.kmode will be used.
mumode (float, None.) – The mu values of each 3D k-mode. In None, self.mumode will be used.

Returns:

fog – The FoG term.

Return type:

float.

get_model_matter_power()[source]: Calculate the model matter power spectrum. The attribute f”_auto_power_matter_model” will be set by the output.

get_model_matter_power_r()[source]: Calculate the model matter power spectrum in real space (without RSD). The attribute f”_auto_power_matter_model_r” will be set by the output.

get_model_power_cross()[source]: Calculate the model cross power spectrum between the two tracers. The attribute f”_cross_power_tracer_model” will be set by the output.

get_model_power_i(i)[source]

Calculate the model power spectrum for the i-th tracer. The attribute f”_auto_power_tracer_{i}_model” will be set by the output.

Parameters:: i (int) – The index of the tracer.
Returns:: auto_power_model – The model power spectrum for the i-th tracer.
Return type:: np.ndarray

get_model_power_noobs_cross()[source]: Calculate the model cross power spectrum between the two tracers without observational effects. The attribute f”_cross_power_tracer_model_noobs” will be set by the output.

get_model_power_noobs_i(i)[source]: Calculate the model power spectrum for the i-th tracer without observational effects. The attribute f”_auto_power_tracer_{i}_model_noobs” will be set by the output.

gridding_compensation()[source]: The sampling window function to be compensated for the gridding mass assignment scheme. Note that the window can only be calculated in Cartesian grids, so it is not used in ModelPowerSpectrum and only in PowerSpectrum.

property include_beam: Whether the beam attenuation is included in the model calculation. Must be a list of two booleans, the first for the first tracer and the second for the second tracer. If just a boolean is provided, it will be used for both tracers.

property kaiser_rsd: Whether RSD is included in the simulation and model calculation. If True, uses the linear Kaiser effect and the FoG profile to compute the model power spectrum.

property kmode: The input kmode for the model calculation.

map_sampling()[source]: The sampling window function from the map cube to be convolved with data. Note that the window can only be calculated in Cartesian grids, so it is not used in ModelPowerSpectrum and only in PowerSpectrum.

property mumode: The mu values of each 3D k-mode.

property sampling_resol: The sampling resolution corresponding to the map-making/gridding of the density field.

property sigma_v_1: The velocity dispersion of the first tracer in km/s.

property sigma_v_2: The velocity dispersion of the second tracer in km/s.

property sigma_z_1: The redshift error of the first tracer.

property sigma_z_2: The redshift error of the second tracer.

property tracer_bias_1: The linear bias of the first tracer.

property tracer_bias_2: The linear bias of the second tracer.

property weights_1: The weights of the first tracer in the rectangular grid.

property weights_2: The weights of the second tracer in the rectangular grid.

property weights_field_1: The weights of the first tracer in the density field.

property weights_field_2: The weights of the second tracer in the density field.

property weights_grid_1: The weights of the first tracer in the rectangular grid.

property weights_grid_2: The weights of the second tracer in the rectangular grid.

class meer21cm.power.PowerSpectrum(field_1=None, box_len=None, weights_field_1=None, weights_grid_1=None, mean_center_1=False, unitless_1=False, field_2=None, weights_field_2=None, weights_grid_2=None, mean_center_2=False, unitless_2=False, k1dbins=None, kmode=None, mumode=None, tracer_bias_1=1.0, sigma_v_1=0.0, tracer_bias_2=None, sigma_v_2=0.0, include_beam=[True, False], fog_profile='lorentz', cross_coeff=1.0, model_k_from_field=True, mean_amp_1=1.0, mean_amp_2=1.0, sampling_resol=None, include_sky_sampling=[True, False], downres_factor_transverse=1.2, downres_factor_radial=2.0, init_box_from_map_data=False, box_buffkick=5, compensate=[False, False], taper_func=<function blackmanharris>, kaiser_rsd=True, grid_scheme='nnb', interlace_shift=0.0, num_particle_per_pixel=1, seed=None, kperpbins=None, kparabins=None, flat_sky=False, flat_sky_padding=[0, 0, 0], k1dweights=None, **params)[source]

Bases: FieldPowerSpectrum, ModelPowerSpectrum

The class for coherently combining the FieldPowerSpectrum and ModelPowerSpectrum classes, and providing an interface for gridding the intensity mapping data as well as the galaxy catalogue to perform power spectrum estimation and for auto-correlation and cross-correlation.

Note that, while you can manually set most of the attributes such as the box length, the density field, the weights, etc., the class is mainly used for first gridding the intensity mapping data and the galaxy catalogue to a rectangular grid field, which then set these attributes automatically. The usual input parameters are those of meer21cm.dataanalysis.Specification.

For usage, check the tutorial notebooks in the examples and cookbook sections.

Parameters:

field_1 (np.ndarray, default None) – The density field of the first tracer.
box_len (list of 3 floats.) – The length of the box along each axis.
weights_field_1 (np.ndarray, default None) – The field-level weights of the first tracer. Default is uniform weights.
weights_grid_1 (np.ndarray, default None) – The grid-level weights of the first tracer. Default is uniform weights.
mean_center_1 (bool, default False) – Whether to mean-center the first tracer field.
unitless_1 (bool, default False) – Whether to divide the first tracer field by its mean.
field_2 (np.ndarray, default None) – The density field of the second tracer. If None, calculation of the second tracer and the cross-correlation will be skipped.
weights_field_2 (np.ndarray, default None) – The field-level weights of the second tracer. Default is uniform weights.
weights_grid_2 (np.ndarray, default None) – The grid-level weights of the second tracer. Default is uniform weights.
mean_center_2 (bool, default False) – Whether to mean-center the second tracer field.
unitless_2 (bool, default False) – Whether to divide the second tracer field by its mean.
k1dbins (list of floats, default None) – The bin edges of k in Mpc-1 for the 1D power spectrum.
kmode (float, default None) – The mode of 3D k in Mpc-1 for the model calculation.
mumode (float, default None) – The mu mode of each 3D k-mode for the model calculation.
tracer_bias_1 (float, default 1.0) – The linear bias of the first tracer.
sigma_v_1 (float, default 0.0) – The velocity dispersion of the first tracer in km/s.
tracer_bias_2 (float, default None) – The linear bias of the second tracer.
sigma_v_2 (float, default 0.0) – The velocity dispersion of the second tracer in km/s.
include_beam (list, default [True, False]) – Whether to include the beam attenuation in the model calculation.
fog_profile (str, default "lorentz") – The shape of the finger-of-god profile to be used in the model calculation.
cross_coeff (float, default 1.0) – The cross-correlation coefficient between the two tracers.
model_k_from_field (bool, default True) – Whether to calculate the model k-mode self.kmode from the field k-mode self.k_mode.
mean_amp_1 (float, default 1.0) – The mean amplitude of the first tracer.
mean_amp_2 (float, default 1.0) – The mean amplitude of the second tracer.
sampling_resol (list, default None) – The sampling resolution of the field in Mpc. If sampling_resol is “auto”, the sampling resolution will be set to the pixel size of the map.
include_sky_sampling (list, default [True, False]) – Whether to include the sky sampling in the model calculation. If just a boolean is provided, it will be used for both tracers.
downres_factor_transverse (float, default None) – The down-sampling factor for the transverse direction of the rectangular box for gridding.
downres_factor_radial (float, default None) – The down-sampling factor for the radial direction of the rectangular box for gridding.
init_box_from_map_data (bool, default False) – If True, the box dimensions will be calculated from the input data cube during initialization. You can always manually set the box dimensions later by self.get_enclosing_box().
box_buffkick (float, default 5) – The buffer kick for the box on each side when gridding. In the unit of Mpc.
compensate (list, default [False, False]) – Whether the gridded fields are compensated according to the mass assignment scheme. Note that the compensation is applied to the model power spectrum, and not to the gridded data fields.
taper_func (function, default windows.blackmanharris) – The taper function to be applied to the gridded field. Note that the taper function is not automatically applied, and is only used when calling PowerSpectrum.apply_taper_to_field().
kaiser_rsd (bool, default True) – Whether to include the RSD effect in the model calculation and mock simulation.
grid_scheme (str, default "nnb") – The grid scheme to be used for gridding. Can be “nnb”, “cic”, “tsc” or “pcs”.
interlace_shift (float, default 0.0) – The shift in the unit of grid cell size for interlacing.
num_particle_per_pixel (int, default 1) – The number of random sampling particles for each sky map pixel.
seed (int, default None) – The seed for the random number generator.
kperpbins (list of floats, default None) – The bin edges of k_perp in Mpc-1 for the cylindrical average power spectrum.
kparabins (list of floats, default None) – The bin edges of k_para in Mpc-1 for the cylindrical average power spectrum.
flat_sky (bool, default False) – Whether to use the flat sky approximation.
flat_sky_padding (list, default [0, 0, 0]) – The padding for the flat sky box.
**params (dict) – Additional parameters to be passed to the base class meer21cm.cosmology.CosmologyCalculator.

Attributes:

As: The amplitude of the primordial power spectrum for the true (fitted) cosmology.
W_HI: A binary window for whether a pixel has been sampled
alpha_AP: The anisotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel} / \alpha_{\perp}\).
alpha_iso: The isotropic Alcock–Paczynski effect parameter, which is \(\alpha_{\parallel}^{1/3} \alpha_{\perp}^{2/3}\).
alpha_parallel: The line of sight Alcock–Paczynski effect parameter.
alpha_perp
auto_power_3d_1: The 3D power spectrum of the first tracer.
auto_power_3d_2: The 3D power spectrum of the second tracer.
auto_power_matter_model: The model matter power spectrum with RSD effects.
auto_power_matter_model_r: The model matter power spectrum in real space (without RSD).
auto_power_tracer_1_model: The 3D model power spectrum for the first tracer.
auto_power_tracer_1_model_noobs: The model power spectrum for the first tracer without observational effects.
auto_power_tracer_2_model: The 3D model power spectrum for the second tracer.
auto_power_tracer_2_model_noobs: The model power spectrum for the second tracer without observational effects.
average_hi_temp: The average HI brightness temperature in Kelvin at central redshift self.z.
average_model_hi_temp: Calculate the average HI brightness temperature in the map cube, taking care of redshift evolution and map sampling.
backend: Which backend to use for computing the matter power.
batch_number: Number of sequential batches used by gridding routines.
beam_image: Returns the beam image projected onto the sky map for the input beam model.
beam_model: The name of the beam function.
beam_type: The beam type that can be either be isotropic or anisotropic.
beam_unit: The unit of input beam size parameter sigma
beam_window_ch: Per-channel spherical-harmonic beam window \(B_\ell\) for HEALPix hp.smoothing workflows.
box_buffkick: The buffer kick for the box on each side when gridding.
box_len: The length of all sides of the box in Mpc.
box_ndim: The number of grids along each side of the enclosing box.
box_origin: The coordinate of the origin of the box in Mpc.
box_resol: The grid length of each side of the enclosing box in Mpc.
box_voxel_redshift: The redshift of each voxel in the rectangular box.
ch_id_gal: The channel id (0-indexed) of each galaxy in the catalogue for cross-correlation.
cold: If True (recommended), the matter power spectrum is the CDM ps.
compensate: Whether the gridded fields are compensated according to the mass assignment scheme.
cosmo: A shortcut to the astropy.cosmology.Cosmology object for the true cosmology.
cospar_fiducial
cospar_true
counts: The number of hits per pixel for the map data
counts_in_box: The counts of the map cube voxels in the rectangular box.
cross_coeff: The cross-correlation coefficient between the two tracers.
cross_power_3d: The 3D cross power spectrum between the two tracers.
cross_power_tracer_model: The 3D model cross power spectrum between the two tracers.
cross_power_tracer_model_noobs: The model power spectrum for the cross correlation between the two tracers without observational effects.
data: The map data
dec_gal: The declination of each galaxy in the catalogue for cross-correlation.
dec_map: The declination of each pixel in the map.
downres_factor_radial: The down-sampling factor for the radial direction of the rectangular box for gridding.
downres_factor_transverse: The down-sampling factor for the transverse direction of the rectangular box for gridding.
dvdf: velocity resolution per unit frequency on average, in km/s/Hz
dvdf_ch: velocity resolution per unit frequency in each channel, in km/s/Hz
expfactor: The expansion factor
f_growth_fiducial: The growth factor at self.expfactor for the fiducial cosmology.
f_growth_true: The growth factor at self.expfactor for the true (fitted) cosmology.
fiducial_cosmology
field_1: The density field of the first tracer.
field_2: The density field of the second tracer.
flat_sky: Whether to use flat sky approximation.
flat_sky_padding: Pad the rectangular box in the flat sky approximation.
fog_profile: The shape of the finger-of-god profile to be used in the model calculation.
fourier_field_1: The Fourier transform of the density field of the first tracer.
fourier_field_2: The Fourier transform of the density field of the second tracer.
freq_gal: The 21cm line frequency for each galaxy in Hz.
freq_resol: frequency resolution in Hz
h: The Hubble constant for the true (fitted) cosmology.
hp_nside: HEALPix \(N_{side}\) for HEALPix-backed specifications.
include_beam: Whether the beam attenuation is included in the model calculation.
interlace_shift: The length in the unit of grid cell size for shifting the gridded field for interlacing.
k_mode: The fiducial (observed) mode of the 3D k-vector.
k_nyquist: The Nyquist frequency of the 3D box along each axis.
k_para: The fiducial parallel k-mode of the 3D box.
k_perp: The fiducial perpendicular k-vector of the 3D box.
k_vec: The 3D k-vector of the box.
kaiser_rsd: Whether RSD is included in the simulation and model calculation.
kmax: The maximum k in Mpc^-1 for calculating matter power.
kmin: The minimum k in Mpc^-1 for calculating matter power.
kmode: The input kmode for the model calculation.
los_resol_in_mpc: effective frequency resolution in Mpc for the fiducial cosmology.
map_has_sampling: A binary window for whether a pixel has been sampled
map_unit_type: The type of the map unit.
matter_power_spectrum_fnc: Interpolation function for the real-space isotropic matter power spectrum.
maximum_sampling_channel: Returns the index of the frequency channel with the maximum sampling on the sky map.
mean_center_1: Whether field_1 needs to be mean centered
mean_center_2: Whether field_2 needs to be mean centered
model_hi_temp_in_box: Based on the redshift dependence of Omega_HI, calculate the average HI brightness temperature for each grid in the rectangular box.
mu_mode: The fiducial (observed) mu values of each k-mode so that \(k_\parallel = k imes \mu\).
mumode: The mu values of each 3D k-mode.
neutrino_mass: The neutrino mass for the true (fitted) cosmology.
ns: The primordial power spectrum spectral index for the true (fitted) cosmology.
nu: The input frequencies of the survey
num_kpoints: The number of k points for calculating matter power.
num_particle_per_pixel: The number of random sampling particles for each sky map pixel.
num_pix_x: The number of pixels along the first map axis.
num_pix_y: The number of pixels along the second map axis.
omega_baryon: The baryon matter density parameter for the true (fitted) cosmology.
omega_cold: The cold matter density parameter for the true (fitted) cosmology.
omega_hi: The HI density as a function of redshift, over the critical density of the Universe at z=0, defined at each frequency channel so that self.omega_hi corresponds to self.z_ch.
omega_hi_z_func: Interpolate the input self.omega_hi at each frequency channel self.z_ch to a function of redshift.
omega_hi_z_mean: The mean HI density at the central redshift self.z, over the critical density of the Universe at z=0.
pix_coor_in_box: The cartesian coordinate of the pixels in Mpc, shifted so that the origin is the origin of the enclosing box.
pix_coor_in_cartesian: The cartesian coordinate of the pixels in Mpc.
pix_resol: angular resolution of the map pixel in deg
pix_resol_in_mpc: angular resolution of the map pixel in Mpc for the fiducial cosmology.
pixel_area: angular area of the map pixel in deg^2
pixel_id: Sparse HEALPix pixel indices (RING, nest=False).
precision: Floating precision flag.
ps_type: linear or nonlinear for the matter power.
ra_gal: The right ascension of each galaxy in the catalogue for cross-correlation.
ra_map: The right ascension of each pixel in the map.
real_dtype: Active real floating dtype controlled by precision.
renorm_ps_1: The renormalization factor of the power spectrum of the first tracer.
renorm_ps_2: The renormalization factor of the power spectrum of the second tracer.
renorm_ps_cross: The renormalization factor of the cross power spectrum.
rot_mat_sky_to_box: The rotational matrix from spheircal cooridnate to regular box.
sampling_resol: The sampling resolution corresponding to the map-making/gridding of the density field.
seed: Seed value for RNG calls throughout the instance.
sigma_beam_ch: The input beam size parameter sigma for each channel.
sigma_beam_ch_in_mpc: The input beam size parameter in Mpc in each channel.
sigma_beam_in_mpc: The channel averaged beam size in Mpc
sigma_v_1: The velocity dispersion of the first tracer in km/s.
sigma_v_2: The velocity dispersion of the second tracer in km/s.
sigma_z_1: The redshift error of the first tracer.
sigma_z_2: The redshift error of the second tracer.
sound_horizon_drag_fiducial: The sound horizon at drag epoch for the fiducial cosmology.
sound_horizon_drag_true: The sound horizon at drag epoch for the true (fitted) cosmology.
survey_volume: Total survey volume in Mpc^3.
tracer_bias_1: The linear bias of the first tracer.
tracer_bias_2: The linear bias of the second tracer.
true_cosmology
unitless_1: Whether field_1 needs to be divided by its mean
unitless_2: Whether field_2 needs to be divided by its mean
vel_resol: velocity resolution on average in km/s
vel_resol_ch: velocity resolution of each channel in km/s
w0: The dark energy equation of state parameter for the true (fitted) cosmology.
w_HI: The weights per map pixel.
wa: The dark energy equation of state parameter for the true (fitted) cosmology.
weights_1: The weights of the first tracer in the rectangular grid.
weights_2: The weights of the second tracer in the rectangular grid.
weights_field_1: The weights of the first tracer in the density field.
weights_field_2: The weights of the second tracer in the density field.
weights_grid_1: The weights of the first tracer in the rectangular grid.
weights_grid_2: The weights of the second tracer in the rectangular grid.
weights_map_pixel: The weights per map pixel.
wproj: The WCS projection object for the map geometry.
x_mode: The mode of the 3D x-vector.
x_vec: The 3D x-vector of the box.
z: The effective centre redshift of the frequency range
z_as_func_of_comov_dist: Returns a function that returns the redshift for input comoving distance.
z_ch: The redshift of each frequency channel
z_gal: The redshifts of each galaxy in the catalogue for cross-correlation.

Methods

`apply_taper_to_field`(field[, taper_func, axis])	Apply a taper to the field, by multiplying the taper function to the corresponding weights of the field.
`beam_attenuation`()	The beam attenuation factor.
`cal_rsd_power`(power_in_real_space, beta1, ...)	Calculate the redshift space power spectrum.
`check_is_map_noiselike_using_pca`(A_mat[, ...])	Use the source mixing matrix from eigendecomposition of the covariance matrix, project out the map data with more and more modes, and check if the variance of the residual map behaves like white noise.
`clean_cache`(attr)	Set the attributes to None.
`convolve_data`([kernel, data, weights, ...])	Convolve map data.
`create_white_noise_map`(sigma_N[, counts, ...])	Create a white noise map with the given standard deviation.
`deltav_to_deltar`(delta_v)	Convert a velocity interval delta_v to a comoving distance interval delta_r.
`deltaz_to_deltar`(delta_z)	Convert a redshift interval delta_z to a comoving distance interval delta_r.
`fog_gaussian`(sigma_r[, kmode, mumode])	The Gaussian finger-of-god profile.
`fog_lorentz`(sigma_r[, kmode, mumode])	The Lorentzian finger-of-god profile.
`fog_term`(sigma_r[, kmode, mumode])	The FoG term for the model calculation.
`gen_random_poisson_galaxy`([sel, num_g_rand, ...])	Generate a random galaxy catalogue from the map cube following the Poisson distribution.
`generate_full_healpix_map`([data, fill_value])	Generate a full healpix map from the data.
`get_1d_power`(power3d[, k1dbins, k1dweights, ...])	Bin the 3D power spectrum into 1D power spectrum.
`get_alpha_parallel`()	Calculate the line of sight Alcock–Paczynski effect parameter.
`get_alpha_perp`()	Calculate the transverse Alcock–Paczynski effect parameter.
`get_beam_image`([wproj, num_pix_x, ...])	Calculate the beam image projected onto the sky map for the input beam model.
`get_beam_window_ch`([ch_sel, cache, lmax, ...])	Build per-channel harmonic beam windows for HEALPix smoothing.
`get_cospar`(cosmology)	Generate a `CosmologyParameters` object from the input cosmology.
`get_counts_in_box`([partial_sel])	Grid the counts of the map cube voxels into the rectangular box, and return the effective counts per rectangular grid.
`get_cy_power`(power3d[, kperpbins, ...])	Bin the 3D power spectrum into cylindrical k_perp-k_para power spectrum.
`get_enclosing_box`([rot_mat])	invoke to calculate the box dimensions for enclosing all the map pixels.
`get_fourier_field_1`()	Calculate the Fourier transform of the density field of the first tracer.
`get_fourier_field_2`()	Calculate the Fourier transform of the density field of the second tracer.
`get_jackknife_patches`(ra_patch_num, ...[, ...])	Split the map into roughly equal patches.
`get_matter_power_spectrum`()	Calculate the matter power spectrum, interpolate it, and save it into the class attribute matter_power_spectrum_fnc.
`get_model_matter_power`()	Calculate the model matter power spectrum.
`get_model_matter_power_r`()	Calculate the model matter power spectrum in real space (without RSD).
`get_model_power_cross`()	Calculate the model cross power spectrum between the two tracers.
`get_model_power_i`(i)	Calculate the model power spectrum for the i-th tracer.
`get_model_power_noobs_cross`()	Calculate the model cross power spectrum between the two tracers without observational effects.
`get_model_power_noobs_i`(i)	Calculate the model power spectrum for the i-th tracer without observational effects.
`get_n_bar_correction`()	Calculate the number density correction for the galaxy catalogue.
`get_sound_horizon_drag`(Omh2, Obh2, Onuh2)	Calculate the sound horizon at drag epoch for a set of cosmological parameters.
`get_weights_none_to_one`(attr_name)	Get the weights, and if it is None, convert it to 1.0 of size of kmode.
`get_z_as_func_of_comov_dist`()	Calculate an array of comoving distances with redshifts, and construct a function that returns the redshift for input comoving distance.
`grid_data_to_field`([flat_sky, partial_sel])	Grid the stored data map to a rectangular grid field.
`grid_field_to_sky_map`(field[, average, ...])	Grid a field in the rectangular box onto the sky.
`grid_gal_to_field`([radecfreq, flat_sky])	Grid the galaxy catalogue to a rectangular grid field.
`gridding_compensation`()	The sampling window function to be compensated for the gridding mass assignment scheme.
`map_sampling`([sampling_resol, p])	The sampling window function from the map cube to be convolved with model power spectrum.
`propagate_field_k_to_model`()	Use field k-modes for the model, taking into account the Alcock–Paczynski effect.
`ra_dec_z_for_coord_in_box`(pos_in_box)	Convert the coordinates in the box to ra, dec, z, and also return the comoving distance to the observer for each point.
`read_from_fits`()	Read in a FITS file for the map data and hit counts.
`read_from_pickle`()	Read in a pickle file for cross-correlation and save the data into the class attributes.
`read_gal_cat`([ra_col, dec_col, z_col, trim])	Read in a galaxy catalogue for cross-correlation and save the data into the class attributes.
`set_corr_type`(corr_type, tracer_indx)	A utility function to help decide whether a tracer field needs to be mean centred, renormalised by its mean, and shot noise removed.
`set_radecnu_bounds_from_map`()	Set `ra_range`, `dec_range`, `nu_min`, and `nu_max` from the loaded `_ra_map`, `_dec_map`, and `nu` (tight RA interval, declination min/max, frequency min/max).
`trim_gal_to_range`()	Trim the galaxy catalogue to the specified range.
`trim_map_to_range`()	Trim the map to the specified range.
`use_flat_sky_box`([flat_sky_padding])	Use flat sky approximation to calculate the box dimensions.

apply_taper_to_field(field, taper_func=None, axis=[2])[source]

Apply a taper to the field, by multiplying the taper function to the corresponding weights of the field.

Parameters:

field (int.) – The index of the field to be tapered, either 1 or 2.
taper_func (function, default None.) – The taper function. Default uses the stored self.taper_func.
axis (list, default [2,].) – The axis to apply the taper to. Default is the z-axis which is approximately the los.

property average_model_hi_temp: Calculate the average HI brightness temperature in the map cube, taking care of redshift evolution and map sampling. Calculation is based on the true (fitted) cosmology.

property box_buffkick: The buffer kick for the box on each side when gridding. In the unit of Mpc.

property box_origin: The coordinate of the origin of the box in Mpc. See meer21cm.grid.minimum_enclosing_box_of_lightcone() for definition.

property box_voxel_redshift: The redshift of each voxel in the rectangular box.

property counts_in_box: The counts of the map cube voxels in the rectangular box.

property downres_factor_radial: The down-sampling factor for the radial direction of the rectangular box for gridding. The box resolution is then multiplied by this factor from the resolution of the frequency channel specified by los_resol_in_mpc. For example, if los_resol_in_mpc is 0.1 Mpc, and downres_factor_radial is 2.0, the box resolution will be 0.2 Mpc.

property downres_factor_transverse: The down-sampling factor for the transverse direction of the rectangular box for gridding. The box resolution is then multiplied by this factor from the resolution of the sky map pixel specified by pix_resol_in_mpc. For example, if pix_resol_in_mpc is 0.1 Mpc, and downres_factor_transverse is 2.0, the box resolution will be 0.2 Mpc.

property flat_sky: Whether to use flat sky approximation. If True, no proper projection and sky rotation is considered. Instead, the sky map cube is assumed to be a rectangular grid with equal voxel size specified by pix_resol_in_mpc and los_resol_in_mpc.

property flat_sky_padding

Pad the rectangular box in the flat sky approximation.

The input should be a list of 3 integers, corresponding to number of padding cells along each dimension in both directions. For example, [1,1,1] will pad 2x2x2 cells.

gen_random_poisson_galaxy(sel=None, num_g_rand=None, seed=None)[source]

Generate a random galaxy catalogue from the map cube following the Poisson distribution. The generation of the sample does not use the instance seed if not explicitly passed and will use a random one otherwise. If you want to generate multiple random catalogues, you need to set a different seed manually for each catalogue.

Parameters:

sel (array, default None) – The selection function of the galaxy catalogue. If None, use the class attribute self.W_HI.
num_g_rand (int, default None) – The number of galaxies to generate. Default uses the number of galaxies stored in the data in self.ra_gal.
seed (int, default None) – The seed for the random number generator.

Returns:

ra_rand (np.ndarray.) – The ra of the random galaxies.
dec_rand (np.ndarray.) – The dec of the random galaxies.
freq_rand (np.ndarray.) – The frequency of the random galaxies. The redshift of the random galaxies can be calculated by meer21cm.util.redshift_to_freq(z_rand).

get_1d_power(power3d, k1dbins=None, k1dweights=None, k_xyz_min=None, k_xyz_max=None, k_perppara_min=None, k_perppara_max=None, multipole_ell=0, mu_model=None)[source]

Bin the 3D power spectrum into 1D power spectrum. If the input power3d is a string, it is assumed to be an attribute of the class, for example auto_power_3d_1. Also see meer21cm.power.bin_3d_to_1d() for more details.

By default the 1D power spectrum is calculated for the monopole. Passing multipole_ell will calculate the 1D power spectrum for the given multipole.

Parameters:

power3d (np.ndarray or str) – The 3D power spectrum.
k1dbins (np.ndarray, default None) – The bins for the 1D power spectrum. Default is the same as the class attribute.
k1dweights (np.ndarray, default None) – The weights for the 3D power spectrum. Default is equal weights for every k-mode.
k_xyz_min (list of size 3, default None) – The minimum k-mode for the 1D power spectrum in x, y, z directions.
k_xyz_max (list of size 3, default None) – The maximum k-mode for the 1D power spectrum in x, y, z directions.
k_perppara_min (list of size 2, default None) – The minimum k_perp and k_para for the 1D power spectrum.
k_perppara_max (list of size 2, default None) – The maximum k_perp and k_para for the 1D power spectrum.
multipole_ell (int, default 0) – The multipole order for the 1D power spectrum. By default the 1D power spectrum is calculated for the monopole.
mu_model (np.ndarray, default None) – The mu-modes for the legendre polynomial. If None, use the class attribute mumode.

Returns:

power1d (np.ndarray) – The 1D power spectrum.
k1deff (np.ndarray) – The effective k-mode for each bin.
nmodes (np.ndarray) – The number of modes in each bin.

get_counts_in_box(partial_sel=None)[source]

Grid the counts of the map cube voxels into the rectangular box, and return the effective counts per rectangular grid.

Parameters:: partial_sel (array, default None) – An additional selection function of the data on top of W_HI. Allows hacking for batch processing.
Returns:: counts_in_grids – The counts of the map cube voxels in the rectangular box.
Return type:: array.

get_cy_power(power3d, kperpbins=None, kparabins=None, kcyweights=None, multipole_ell=0, mu_model=None)[source]

Bin the 3D power spectrum into cylindrical k_perp-k_para power spectrum. If the input power3d is a string, it is assumed to be an attribute of the class, for example auto_power_3d_1. Also see meer21cm.power.bin_3d_to_cy() for more details.

Passing multipole_ell will calculate the cylindrical power spectrum multiplied by the Legendre polynomial.

Parameters:

power3d (np.ndarray or str) – The 3D power spectrum.
kperpbins (np.ndarray, default None) – The k_perp bins for the cylindrical ps. Default is the same as the class attribute.
kparabins (np.ndarray, default None) – The k_para bins for the cylindrical ps. Default is the same as the class attribute.
kcyweights (np.ndarray, default None) – The weights for the 3D power spectrum. Default is equal weights for every k-mode.
multipole_ell (int, default 0) – The multipole order for the cylindrical power spectrum. By default the cylindrical power spectrum is calculated for the monopole.
mu_model (np.ndarray, default None) – The mu-modes for the legendre polynomial. If None, use the class attribute mumode.

Returns:

powercy (np.ndarray) – The cylindrical power spectrum.
weightscy (np.ndarray) – The weights for the cylindrical k-modes.

get_enclosing_box(rot_mat=None)[source]

invoke to calculate the box dimensions for enclosing all the map pixels.

Parameters:: rot_mat (np.ndarray, default None) – The rotational matrix from the sky to the box. If None, calculates the suitable rotation matrix automatically.

get_n_bar_correction()[source]: Calculate the number density correction for the galaxy catalogue.

grid_data_to_field(flat_sky=None, partial_sel=None)[source]

Grid the stored data map to a rectangular grid field.

If flat_sky is True, no gridding is performed. Instead, the map cube dimensions are taken to be a rectangular grid, with the grid length corresponding to the pixel resolution on x/y and los frequency resolution as z.

If flat_sky is False, the data is gridded onto a regular grid using the input grid scheme and performing the proper curved sky projection.

The gridded field is stored as field_1 and the weights are stored as weights_1.

Parameters:

flat_sky (bool, default None) – If True, use flat sky approximation.
partial_sel (array, default None) – An additional selection function of the data on top of W_HI. Allows hacking for batch processing.

grid_field_to_sky_map(field, average=True, mask=True, wproj=None, num_pix_x=None, num_pix_y=None, los_sel=None)[source]

Grid a field in the rectangular box onto the sky.

Parameters:

field (array.) – The field in the box to be gridded.
average (bool, default True.) – Whether the field grids are averaged or summed into sky pixels.
mask (bool, default True.) – If True, the sky map is then masked by the survey selection function.
wproj (astropy.wcs.WCS object, default None.) – WCS only. The WCS for the output sky map. Default uses self.wproj.
num_pix_x (int, default None.) – WCS only. Number of pixels along the first map axis. Default self.num_pix_x.
num_pix_y (int, default None.) – WCS only. Number of pixels along the second map axis. Default self.num_pix_y.
los_sel (array-like, default None.) – Optional selector of line-of-sight (last-axis) indices represented by field. If None, field is assumed to contain all LOS slices with shape self.box_ndim. If provided, field must have shape (self.box_ndim[0], self.box_ndim[1], len(los_sel)) and the projected output still uses the full LOS map axis (self.nu.size), so multiple chunked calls can be merged by accumulating mass/count outputs.

Returns:

map_bin (array.) – The output sky map (WCS (nx, ny, n_ch) or HEALPix (n_pix, n_ch)).
count_bin (array.) – Per-cell accumulation used for averaging (WCS) or voxel counts (HEALPix).

grid_gal_to_field(radecfreq=None, flat_sky=None)[source]

Grid the galaxy catalogue to a rectangular grid field.

If flat_sky is True, no gridding is performed. Instead, the map cube dimensions are taken to be a rectangular grid, with the grid length corresponding to the pixel resolution on x/y and los frequency resolution as z.

gridding_compensation()[source]: The sampling window function to be compensated for the gridding mass assignment scheme.

property interlace_shift: The length in the unit of grid cell size for shifting the gridded field for interlacing. 0 corresponds to no interlacing.

map_sampling(sampling_resol=None, p=1)[source]

The sampling window function from the map cube to be convolved with model power spectrum. This should correspond to the resolution of map-making on the sky and the frequency channel, not the resolution of the gridded field.

Parameters:

sampling_resol (list, default None) – The sampling resolution of the field in Mpc. Default is the class attribute sampling_resol.
p (int, default 1) – The index of assignment scheme.

Returns:

B_sampling – The sampling window function in 3D k-space.

Return type:

np.ndarray.

property model_hi_temp_in_box: Based on the redshift dependence of Omega_HI, calculate the average HI brightness temperature for each grid in the rectangular box. This can be used in a way that the weighted average of the model_hi_temp_in_box is used as the average t_bar in the model power spectrum (by passing it to mean_amp_1), whereas the 3D model_hi_temp_in_box is used as the field weight to account for the effect of the redshift evolution of Omega_HI in the power spectrum. Calculation is based on the true (fitted) cosmology.

property num_particle_per_pixel: The number of random sampling particles for each sky map pixel.

property pix_coor_in_box: The cartesian coordinate of the pixels in Mpc, shifted so that the origin is the origin of the enclosing box.

property pix_coor_in_cartesian: The cartesian coordinate of the pixels in Mpc.

propagate_field_k_to_model()[source]: Use field k-modes for the model, taking into account the Alcock–Paczynski effect.

\[ \begin{align}\begin{aligned}k_\perp^\text{fiducial} = k_\perp^\text{true} \times \alpha_\perp\\k_\parallel^\text{fiducial} = k_\parallel^\text{true} \times \alpha_\parallel\end{aligned}\end{align} \]

ra_dec_z_for_coord_in_box(pos_in_box)[source]

Convert the coordinates in the box to ra, dec, z, and also return the comoving distance to the observer for each point.

Parameters:

pos_in_box (array.) – The coordinates in the box.

Returns:

pos_ra (array.) – The ra of the points.
pos_dec (array.) – The dec of the points.
pos_z (array.) – The redshift of the points.
pos_comov_dist (array.) – The comoving distance to the observer for each point.

property rot_mat_sky_to_box

The rotational matrix from spheircal cooridnate to regular box.

See meer21cm.grid.minimum_enclosing_box_of_lightcone() for definition.

property seed: Seed value for RNG calls throughout the instance.

use_flat_sky_box(flat_sky_padding=None)[source]

Use flat sky approximation to calculate the box dimensions.

Parameters:: flat_sky_padding (list, default None) – The padding for the flat sky box. If None, use the class attribute flat_sky_padding.

meer21cm.power.bin_3d_to_1d(ps3d, kfield, k1dedges, weights=None, error=False, vectorize=False)[source]

Bin a 3d distribution, e.g. power spectrum \(P_{3D}(\vec{k})\), into 1D average.

Note that, the distribution is unraveled to a 1D array, so essentially an array of any dimension would do, as long as ps3d, kfield, and weights have the same size.

The mean of the 1D average is calculated as

\[\hat{P}_{\rm 1D}^{i} = \big(\sum_j P_{\rm 3D}^{ j} w_{ j} \big)/\big(\sum_j w_{ j}\big),\]

where j loops over all the modess that fall into the \(i^{\rm th}\) bin and \(w_{ j}\) is the weights.

If error is set to True, a sampling error is calculated and returned so that

\[(\Delta P_{\rm 1D}^{\rm i})^2 = \big(\sum_j (P_{\rm 3D}^{\rm j}-\hat{P}_{\rm 1D}^{\rm i})^2 w_{\rm j}^2 \big) \Big/ \big(\sum_j w_{\rm j}\big)^2.\]

Parameters:

ps3d (np.ndarray) – The 3D distribution to be binned.
kfield (np.ndarray) – The k-field of the 3D distribution.
k1dedges (np.ndarray) – The bin edges for the 1D power spectrum.
weights (np.ndarray, default None) – The weights for each 3D k-mode of the power spectrum.
error (bool, default False) – Whether to calculate the sampling error.
vectorize (bool, default False) – Whether to vectorize the calculation, assuming the first axis is independent realisations.

Returns:

ps1d (np.ndarray) – The 1D power spectrum.
ps1derr (np.ndarray) – The sampling error for the 1D power spectrum. Returned only if error is True.
k1deff (np.ndarray) – The effective k-mode for each bin.
nmodes (np.ndarray) – The number of modes in each bin.

meer21cm.power.bin_3d_to_cy(ps3d, kperp_i, kperpedges, weights=None, average=True, vectorize=False)[source]

Function to bin a 3D distribution (e.g. power spectrum) into cylindrical average. The arrays are unravelled to 2D before binning by keeping the last axis. The 2D array is then binned along the first axis.

The output is flipped so that the first axis is the original last axis. Therefore, to bin a 3D power spectrum to a cylindrical average, one can simply run bin_3d_to_cy twice (see PowerSpectrum.get_cy_power).

Parameters:

ps3d (array.) – The 3D distribution to be binned.
kperp_i (array.) – The perpendicular k-mode corresponding to the first axis.
kperpedges (array.) – The bin edges for the perpendicular k-mode.
weights (array, None.) – The weights for each 3D k-mode of the power spectrum.
average (bool, default True.) – If True, calculate the weighted average of the power spectrum in each bin. Else, calculate the weighted sum.
vectorize (bool, default False) – Whether to vectorize the calculation, assuming the first axis is independent realisations.

Returns:

pscy – The cylindrical average of the 3D distribution.

Return type:

np.ndarray

meer21cm.power.gaussian_beam_attenuation(k_perp, beam_sigma_in_mpc)[source]: The beam attenuation term to be multiplied to model power spectrum assuming a Gaussian beam.

meer21cm.power.get_fourier_density(real_field, weights=None, mean_center=False, unitless=False, norm='forward')[source]

Perform Fourier transform of a density field in real space. Note that this is deliberately written in a way that is not dimension specific. It can be used to calculate power spectrum of arbitrary dimension.

Note that, the field is multiplied by the weights and then Fourier-transformed, and is not weight normalised.

Parameters:

real_field (np.ndarray) – The real-space field.
weights (np.ndarray, default None) – The weights of the field.
mean_center (bool, default False) – Whether to mean center the field.
unitless (bool, default False) – Whether to make the field unitless.
norm (str, default "forward") – The normalization of the Fourier transform. Naming is the same as np.fft.

Returns:

fourier_field – The Fourier transform of the field.

Return type:

np.ndarray

meer21cm.power.get_gaussian_noise_floor(sigma_n, box_ndim, box_volume=1.0, counts=None)[source]

Calculate the Gaussian noise floor for a given field.

Parameters:

sigma_n (float) – The standard deviation of the noise before being averaged down by the sampling.
box_ndim (tuple) – The number of grids of the box along each direction.
box_volume (float, default 1.0) – The volume of the box.
counts (np.ndarray, default None) – The number of sampling in the box. If None, it is set to be 1.0.

Returns:

noise_floor – The noise floor.

Return type:

float

meer21cm.power.get_k_vector(box_ndim, box_resol)[source]

Get the wavenumber vector along each direction for a given box.

Parameters:

box_ndim (int) – The number of dimensions of the box.
box_resol (float) – The resolution of the box.

Returns:

kvecarr – The wavenumber vector along each direction.

Return type:

tuple

meer21cm.power.get_modelpk_conv(psmod, weights1_in_real=None, weights2=None, renorm=True)[source]

Convolve a model power spectrum with real-space weights.

Parameters:

psmod (np.ndarray) – The model power spectrum.
weights1_in_real (np.ndarray, default None) – The real-space weights for the first field. Default is None, which means no weights.
weights2 (np.ndarray, default None) – The real-space weights for the second field. Default is None, which means no weights.
renorm (bool, default True) – Whether to renormalize the power spectrum.

Returns:

power_conv – The convolved power spectrum.

Return type:

np.ndarray

meer21cm.power.get_power_spectrum(fourier_field, box_len, weights=None, field_2=None, weights_2=None, renorm=True)[source]

Calculate the power spectrum for one/two given Fourier fields.

Parameters:

fourier_field (np.ndarray) – The Fourier field of the first tracer.
box_len (tuple) – The length of the box along each direction.
weights (np.ndarray, default None) – The weights of the first tracer in real space.
field_2 (np.ndarray, default None) – The Fourier field of the second tracer. If None, it is set to be the same as the first field.
weights_2 (np.ndarray, default None) – The weights of the second tracer in real space. Must be provided if field_2 is provided.
renorm (bool, default True) – Whether to renormalize the power spectrum by the weights.

Returns:

power – The power spectrum.

Return type:

np.ndarray

meer21cm.power.get_renormed_field(real_field, weights=None, mean_center=False, unitless=False)[source]

Mean center the field and renormalise it by dividing the mean.

Parameters:

real_field (np.ndarray) – The real-space field.
weights (np.ndarray, default None) – The weights of the field.
mean_center (bool, default False) – Whether to mean center the field.
unitless (bool, default False) – Whether to make the field unitless.

Returns:

field – The renormalized field.

Return type:

np.ndarray

meer21cm.power.get_shot_noise(real_field, box_len, weights=None)[source]

Calculate the shot noise of a field.

Parameters:

real_field (np.ndarray) – The real-space field.
box_len (tuple) – The length of the box along each direction.
weights (np.ndarray, default None) – The weights of the field.

Returns:

shot_noise – The shot noise of the field.

Return type:

float

meer21cm.power.get_shot_noise_galaxy(gal_count, box_len, weights_grid=None, weights_field=None)[source]: Calculate the shot noise of a galaxy number count field.

meer21cm.power.get_vec_mode(vecarr)[source]

Calculate the mode of the n-dimensional vectors on the grids

Parameters:: vecarr (tuple) – The vectors.
Returns:: mode – The mode of the vectors.
Return type:: np.ndarray

meer21cm.power.get_x_vector(box_ndim, box_resol)[source]

Get the position vector along each direction for a given box.

Parameters:

box_ndim (int) – The number of dimensions of the box.
box_resol (float) – The resolution of the box.

Returns:

xvecarr – The position vector along each direction.

Return type:

tuple

meer21cm.power.power_weights_renorm(weights1_in_real=None, weights2=None)[source]

Calculate the renormalization coefficient based on the weights on the density field when calculating power spectrum. The renormalization is defined as

\[\frac{{N_{\rm grid}}} {\sum_{i} w_1(x_i) w_2(x_i)},\]

where \(N_{\rm grid}\) is the number of grids in the box and \(i\) loops over all the grids.

Note that this renormaliszation corresponds to the diagonal renormalisation matrix that does not change the window function convolution, but only renormalises the sum of each row of the window function matrix. See Chen (2025) [1] for more details.

Parameters:

weights1_in_real (array, default None.) – The weights of the density field in real space. Must be in the shape of the regular grid field.
weights2 (array, default None.) – If cross-correlation, the weights for the second field.

Returns:

weights_norm – The renormalization coefficient.

Return type:

float.

References

[1] Chen, Z., 2025, “A quadratic estimator view of the transfer function correction in intensity mapping surveys”, https://ui.adsabs.harvard.edu/abs/2025MNRAS.542L…1C/abstract.

meer21cm.power.step_window_attenuation(k_dir, step_size_in_mpc, p=1)[source]: The beam attenuation term to be multiplied to model power spectrum assuming a Gaussian beam.

meer21cm.stack module

meer21cm.stack.stack(sp, stack_angular_num_nearby_pix=10, symmetrize=False)[source]

Calculate a stacked 3D cubelet using the intensity maps and source positions stored in sp. The stacked signal in the cubelet is calculated as:

\[S(\Delta\alpha,\Delta\phi,\Delta\nu) = \frac{\sum_i S(\alpha_i{+}\Delta\alpha,\,\phi_i{+}\Delta\phi,\,\nu_i{+}\Delta\nu) w_i}{\sum_i w_i},\]

where i loops over each source stored in the (sp.ra_gal, sp.dec_gal, sp.z_gal) positions. The pixels each source falls into are at (\(\alpha_i\), \(\phi_i\), \(\nu_i\)). \(S\) is the intensity map stored in sp.data. \(w\) is the weights stored in sp.weights_map_pixel and \(w_i\) is the weights at \((\alpha_i{+}\Delta\alpha,\,\phi_i{+}\Delta\phi,\,\nu_i{+}\Delta\nu)\).

The cubelet extends towards the entire frequency range of the map so \(\delta\nu\) is sampled at [\(-N_{\rm ch} \delta\nu\),…,0,…, \(N_{\rm ch} \delta\nu\)]. The angular sampling of the cubelet corresponds to the map pixels, and the size of the angular plane is set by stack_angular_num_nearby_pix. Note that stack_angular_num_nearby_pix is the number of pixels each side of the centre so the size of the angular plane is (2 * stack_angular_num_nearby_pix + 1)**2.

If symmetrize, a mirroring of the individual cubelets is performed along \(\Delta\nu=0\). This corresponds to the 180deg rotation along the spectral axis described in Sinigaglia et al. (2022) [1] and is the only symmetry that single-dish IM stacking is sensitive to.

Parameters:

sp (meer21cm.Specification object.) – The data used for stacking.
stack_angular_num_nearby_pix (optional, default 10.) – The number of map pixels sampled on each side relative to the source centre.
symmetrize (optional, default False.) – Whether to symmetrize the stacking.

Returns:

stack_3D_map (array.) – The averaged cubelet for the stacking.
stack_3D_weight (array.) – The weights in each voxel in the averaged cubelet.

References

[1]

Sinigaglia, F. et al., “Optimizing spectral stacking for 21-cm observations of galaxies: accuracy assessment and symmetrized stacking”, https://ui.adsabs.harvard.edu/abs/2022MNRAS.514.4205S.

meer21cm.stack.sum_3d_stack(stack_3D_map, vel_ch_avg=5, ang_sum_dist=3.0)[source]

Collapse a stacked cubelet into stacked image and stacked spectrum.

Note that for stacked image, vel_ch_avg is the number of channels that go into the summation on each side of the centre channel so that the total number of channels that are summed is (2 * vel_ch_avg + 1).

For stacked spectrum, the angular pixels that go into the summation are determined by the distance to the center pixel. Note that the distance is in cell length not physical angular unit.

Parameters:

stack_3D_map (array.) – The stacked cubelet.
vel_ch_avg (optional, default 5.) – How many channels on each side of the center to sum into stacked image.
ang_sum_dist (optional, default 3.0.) – The distance within which the angular pixels are summed to stacked spectrum

Returns:

angular_stack_map (array.) – The stacked image.
spectral_stack_map (array.) – The stacked spectrum.

meer21cm.telescope module

meer21cm.telescope.cmb_temperature(nu, tcmb0=np.float64(2.7255))[source]

Calculate the background CMB temperature at given frequencies.

Parameters:

nu (float.) – The observing frequency in Hz.
tcmb0 (float, default Planck18.Tcmb0.value.) – The background CMB temperature at z=0.

Returns:

tcmb – The CMB temperature at given frequencies in Kelvin.

Return type:

float.

meer21cm.telescope.cos_beam(sigma)[source]

Returns a cosine-tapered beam function [1]

\[B(\theta) = \bigg[ \frac{\cos \big( 1.189 \pi \theta / \theta_b \big)} {1-4\big( 1.189 \theta / \theta_b \big)} \bigg]^2\]

for given input parameter \(\sigma\), the FWHM is set to \(\theta_b = 2\sqrt{2{\rm log}2 \sigma}\).

Parameters:: sigma (float.) – The width of the beam profile.
Returns:: beam_func – The beam function.
Return type:: function.

References

[1]

Mauch et al., “The 1.28 GHz MeerKAT DEEP2 Image”, https://arxiv.org/abs/1912.06212

meer21cm.telescope.dish_beam_sigma(dish_diameter, nu, gamma=1.0, ang_unit=Unit('deg'))[source]

Calculate the beam size of a dish telescope assuming

\[\theta_{\rm FWHM} = \gamma \frac{\lambda}{D},\]

where \(\theta_{\rm FWHM}\) is the FWHM of the beam, \(\gamma\) is the aperture efficiency, \(\lambda\) is the observing wavelength, and D is the dish diameter.

The sigma of the Gaussian beam is then \(\sigma = \theta_{\rm FWHM}/ 2\sqrt{2 {\rm ln}2}\).

Parameters:

dish_diameter (float.) – The diameter of the dish in metre.
nu (float.) – The observing frequency in Hz.
gamma (float, default 1.0.) – The aperture efficiency.
ang_unit (str or astropy.units.Unit, default deg.) – The unit of the output.

Returns:

beam_sigma – The sigma of the beam.

Return type:

float.

meer21cm.telescope.galaxy_temperature(nu, tgal_408MHz=25, sp_indx=-2.75)[source]

The temperature template of the Milky Way.

Note that, for an accurate T_sky, you can instead use meer21cm.fg.ForegroundSimulation to generate the foregrounds.

Parameters:

nu (float.) – The observing frequency in Hz.
tgal_408MHz (float.) – The average galaxy temperature at 408MHz in Kelvin.
sp_indx (float.) – The spectral index to extrapolate it to input frequencies.

Returns:

Tgal – The galaxy temperature at given frequencies in Kelvin.

Return type:

float.

meer21cm.telescope.gaussian_beam(sigma)[source]

Returns a Gaussian beam function

\[B(\theta) = {\rm exp}[-\frac{\theta^2}{2\sigma^2}]\]

when the beam width \(\sigma\) is specified.

Parameters:: sigma (float.) – The width of the gaussian beam profile.
Returns:: beam_func – The beam function.
Return type:: function.

meer21cm.telescope.gaussian_beam_window(sigma_rad, lmax)[source]

Closed-form Gaussian beam window function in spherical-harmonic space.

\[B_\ell = \exp\left(-\frac{\ell(\ell+1)\sigma^2}{2}\right),\]

where \(\sigma\) is the Gaussian beam dispersion in radians.

Parameters:

sigma_rad (float) – Gaussian beam dispersion in radians.
lmax (int) – Maximum multipole; the returned array has length lmax + 1.

Returns:

b_ell

Return type:

ndarray of shape (lmax + 1,)

meer21cm.telescope.get_beam_xy(wproj, xdim, ydim)[source]: Get the x and y angular coordinates of the given wcs.

meer21cm.telescope.isotropic_beam_profile(xdim, ydim, wproj, beam_func, ang_unit=Unit('deg'))[source]

Generate an isotropic image of the beam for given wproj and beam_func. The image can later be used to convolve or deconvolve with intensity maps.

Parameters:

xdim (int.) – The number of pixels in the first axis.
ydim (int.) – The number of pixels in the second axis.
wproj (astropy.wcs.WCS object.) – The two-dimensional wcs object for the map.
beam_func (function.) – The beam function.
ang_unit (str or astropy.units.Unit.) – The unit of input values for beam_func.

Returns:

beam_image – The image of the beam.

Return type:

float array.

meer21cm.telescope.isotropic_beam_window(beam_func, sigma, lmax, theta_grid=None)[source]

Numerical beam window function for an isotropic beam via healpy.sphtfunc.beam2bl.

The beam profile beam_func(theta) is sampled on a fine theta grid covering \([0, \theta_{\rm max}]\) with \(\theta_{\rm max} \approx 8\sigma\) (truncated at \(\pi\)), and the resulting \(B(\theta)\) is converted to \(B_\ell\) with hp.sphtfunc.beam2bl.

Parameters:

beam_func (callable) – Maps theta (radians or the ang_unit chosen by the caller’s beam factory) to beam amplitude. Must accept vector inputs.
sigma (float) – Beam dispersion in the same angular units as beam_func expects. Used to pick the default theta_grid when none is provided.
lmax (int) – Maximum multipole.
theta_grid (ndarray, optional) – Explicit theta sample points (same units as expected by beam_func). Must start at 0 and be monotonically increasing.

Returns:

b_ell

Return type:

ndarray of shape (lmax + 1,)

meer21cm.telescope.kat_beam(nu, wproj, xdim, ydim, band='L')[source]

Returns a beam model from the katbeam model, which is a simplification of the model reported in Asad et al. [1]. The katbeam implementation here still needs validation. Use it with caution, especially if you want correct orientation of the beam.

References

[1]

Asad et al., “Primary beam effects of radio astronomy antennas – II. Modelling the MeerKAT L-band beam”, https://arxiv.org/abs/1904.07155

meer21cm.telescope.receiver_temperature_meerkat(nu)[source]

The receiver temperature of MeerKAT.

Parameters:: nu (float.) – The observing frequency in Hz.
Returns:: Trx – The receiver temperature at given frequencies in Kelvin.
Return type:: float.

meer21cm.telescope.weighted_convolution(signal, kernel, weights, kernel_renorm=True, los_axis=-1)[source]

Perform weighted convolution of signal. The weighted convolution of the signal is defined as

\[\tilde{s} = [(s \cdot w) * b]/[w * b],\]

where \(s\) is the signal, \(w\) is the weight, \(b\) is the convolution kernel and \(w\) denotes convolution.

The convolution also creates new weights for the output signal so that

\[\tilde{w} = [w * b]^2 / [w * b^2]\]

Parameters:

signal (float.) – The input signal to be convolved
kernel (float.) – The convolution kernel
weights (float.) – The weights for the signal
kernel_renorm (boolean, default True.) – Whether to renormalise the kernel so that the sum of the kernel is one. Should be set to True for temperature and False for flux density.
los_axis (int, default -1.) – which axis is the los.

Returns:

conv_signal (float.) – The convolved signal.
conv_weights (float.) – The convolved weights.

meer21cm.telescope.weighted_smoothing_healpix(data_pix, weights_pix, beam_window_ch, hp_nside, pixel_id, kernel_renorm=True, nside_out=None, pixel_id_out=None)[source]

Weighted HEALPix smoothing, mirroring weighted_convolution() semantics.

For each channel \(c\) with beam window \(B^{(c)}_\ell\), this returns

\[\tilde s^{(c)} = \frac{\mathcal{S}[s^{(c)} w^{(c)} \, B^{(c)}_\ell]} {\mathcal{S}[w^{(c)} \, B^{(c)}_\ell]},\qquad \tilde w^{(c)} = \frac{\bigl(\mathcal{S}[w^{(c)} B^{(c)}_\ell]\bigr)^2} {\mathcal{S}[w^{(c)} (B^{(c)}_\ell)^2]},\]

where \(\mathcal{S}\) denotes harmonic-space smoothing via hp.smoothing(..., beam_window=B_\ell) on a full-sphere scratch buffer (allocated per channel, not kept).

Parameters:

data_pix ((n_pix, n_ch) ndarray)
weights_pix ((n_pix, n_ch) ndarray)
beam_window_ch ((n_ch, lmax+1) ndarray or (lmax+1,) ndarray) – Per-channel beam window. A 1D array is broadcast over all channels.
hp_nside (int) – HEALPix \(N_{\rm side}\) that pixel_id refers to.
pixel_id (ndarray of int) – HEALPix indices at hp_nside, matching rows of data_pix / weights_pix, used to paint the high-resolution sphere before smoothing.
pixel_id_out (ndarray of int, optional) – If nside_out is below hp_nside, indices at nside_out where outputs are sampled (survey footprint). Required in that case. When no downgrade, defaults to pixel_id.
kernel_renorm (bool, default True) – Kept for API symmetry with weighted_convolution(). B_\ell is already normalised by hp.smoothing (B_0 = 1 for Gaussian windows), so this flag is informational only.
nside_out (int, optional) – If set to a HEALPix \(N_{\\rm side}\) strictly less than hp_nside, harmonic smoothing is carried out at hp_nside on full spheres, maps are hp.ud_grade down to nside_out, then values are sampled at pixel_id_out.

Returns:

conv_data_pix ((n_pix_out, n_ch) ndarray)
conv_weights_pix ((n_pix_out, n_ch) ndarray)

meer21cm.transfer module

This module contains the transfer function analysis class and related functions.

In data analysis, you should have already used meer21cm.power.PowerSpectrum to calculate the power spectrum of the data, or meer21cm.mock.MockSimulation to generate the mock data for power spectrum estimation. In that case, you can pass the meer21cm.power.PowerSpectrum or meer21cm.mock.MockSimulation instance as an input to meer21cm.transfer.TransferFunction. The transfer function class then takes into account of the settings of the power spectrum instance, and perform mock simulations to estimate the numerical transfer function, keeping the consistency of the gridding, beam, and other settings.

Finally note that, the transfer function calculation supports parallelisation using MPI, by using mpi4py.futures.MPIPoolExecutor. The default installation of meer21cm does not include mpi4py, and you need to configure MPI first (usually with openmpi or mpich), and then install mpi4py manually.

class meer21cm.transfer.TransferFunction(ps: PowerSpectrum, N_fg: int, R_mat: ndarray | None = None, uncleaned_data: ndarray | None = None, highres_sim: int | None = None, upres_transverse: float = 4, upres_radial: float = 4, mean_center_map: bool = True, pca_map_weights: ndarray | None = None, parallel_plane: bool = True, rsd_from_field: bool = False, discrete_source_dndz: list | Callable = <function ones_like>, pool: str = 'multiprocessing', num_process: int | None = None, unmask_during_mock: bool = False)[source]

Bases: object

A class to perform transfer function analysis.

By default, the transfer function is calculated by generating mock HI signal, and perform PCA cleaning by injecting the mock HI signal into the data (unless R_mat is provided, in which case no PCA is performed and R_mat is used to clean the injected data). The transfer function is then calculated by dividing the 1D power spectrum of the cleaned data by the 1D power spectrum of the mock data. A minimum example running 10 realizations:

>>> tf = TransferFunction(ps, N_fg=3)
>>> tf1d_arr = tf.run(range(10), type="auto") # gives you the 1D transfer function with 10 realizations

If type is "auto", the numerator is the cross-correlation between cleaned HI mock and the original HI mock (and therefore the result is based on HI “auto” power), and the denominator is the auto power of the original HI mock. If type is "cross", the numerator is the cross-correlation between cleaned HI mock and galaxy mock, and the denominator is the cross power of the original HI mock and galaxy mock.

Additionally, if type is "null", instead of running transfer function calculation, no HI mock is generated, and only mock galaxy data is used to cross-correlate with the map data as a null test.

Parameters:

ps (meer21cm.power.PowerSpectrum) – The power spectrum instance.
N_fg (int) – The number of eigenmodes for PCA cleaning.
R_mat (np.ndarray, default None) – The PCA matrix. If not provided, it will be calculated from the data. If provided, no map injection is performed and the cleaning uses this matrix.
uncleaned_data (np.ndarray, default None) – The uncleaned data, which the mock HI signal is to be injected into. If not provided, the mock HI signal is injected into ps.data. Note that this can be wrong if ps.data is already overridden by the residual map. In that case, you should provide the original map.
highres_sim (int, default Nones) – The mock field will be gridded to a high resolution sky map, with the resolution specified by this integer times the pixel resolution, and then down-sampled to the pixel resolution.
upres_transverse (float, default 4) – The mock field is first generated in a rectangular box. This is the up-sampling factor of the box comparing to the pixel resolution in the transverse direction.
upres_radial (float, default 4) – The up-sampling factor in the radial direction comparing to the frequency resolution.
mean_center_map (bool, default True) – Whether to mean center the map data for PCA cleaning.
pca_map_weights (np.ndarray, default None) – If mean centering is performed, this is the weights for the mean calculation. Also used to weight the data for covariance calculation for PCA eigendecomposition.
parallel_plane (bool, default True) – Whether to simulate the mock field in the parallel-plane limit (i.e. k_z = k_parallel). Only used if rsd_from_field is True, otherwise parallel-plane is hard-coded to True.
rsd_from_field (bool, default False) – Whether to generate the RSD effect at field level. If False, the RSD effect is generated at power spectrum level.
discrete_source_dndz (list or function, default np.ones_like) – The redshift distribution of the discrete tracer sources. Must be a function of redshift. Note that the overall number of discrete sources is set to ps.ra_gal.size or ps.num_discrete_source, so only the shape of the dndz is used, not the normalization. If a list is provided, the first element is the redshift array, and the second element is the dndz array, and a 1D linear interpolation is used to get the dndz function.
pool (str, default "multiprocessing") – The pool to use for parallelisation. Can be “multiprocessing” or “mpi”.
num_process (int, default None) – The number of processes to use for parallelisation. If not provided, the number of processes is set to the number of cores available.
unmask_during_mock (bool, default False) – Whether to unmask the wcs during the mock simulation. If True, the survey area will be unmasked during the mock simulation, and then masked again after mock HI is generated. This is useful for mitigating edge effects and when box buffer is small. For a large patch with a wide z-range, it is recommended to set this to True, otherwise the mock HI map may contain nan and cause errors.

Methods

`get_arg_list_for_parallel_auto`(seed_list[, ...])	Generate a list of arguments for parallelisation of the auto-correlation runs (i.e. the transfer function is calculated between cleaned mock HI and original mock HI).
`get_arg_list_for_parallel_cross`(seed_list[, ...])	Generate a list of arguments for parallelisation of the cross-correlation runs (i.e. the transfer function is calculated between mock HI and mock galaxy).
`get_arg_list_for_parallel_null`(seed_list[, ...])	Generate a list of arguments for parallelisation of the null test runs.
`get_mock_instance_attr_dict`(seed)	Generate the attribute dictionary for the mock instance.
`run`(seed_list[, type, return_power_3d, ...])	Run the transfer function calculation.

get_arg_list_for_parallel_auto(seed_list, return_power_3d=False, return_power_1d=False)[source]

Generate a list of arguments for parallelisation of the auto-correlation runs (i.e. the transfer function is calculated between cleaned mock HI and original mock HI). This list is then used for pool.starmap.

Parameters:

seed_list (list) – The list of seeds for the mock instances.
return_power_3d (bool, default False) – Whether to return the 3D power spectrum of the auto-correlation.
return_power_1d (bool, default False) – Whether to return the 1D power spectrum of the auto-correlation.

Returns:

arg_list – The list of arguments for the parallelisation.

Return type:

list

get_arg_list_for_parallel_cross(seed_list, return_power_3d=False, return_power_1d=False)[source]

Generate a list of arguments for parallelisation of the cross-correlation runs (i.e. the transfer function is calculated between mock HI and mock galaxy). This list is then used for pool.starmap.

Parameters:

seed_list (list) – The list of seeds for the mock instances.
return_power_3d (bool, default False) – Whether to return the 3D power spectrum of the cross-correlation.
return_power_1d (bool, default False) – Whether to return the 1D power spectrum of the cross-correlation.

Returns:

arg_list – The list of arguments for the parallelisation.

Return type:

list

get_arg_list_for_parallel_null(seed_list, return_power_3d=False, return_power_1d=False)[source]

Generate a list of arguments for parallelisation of the null test runs. This list is then used for pool.starmap.

Parameters:

seed_list (list) – The list of seeds for the mock instances.
return_power_3d (bool, default False) – Whether to return the 3D power spectrum of the null test.
return_power_1d (bool, default False) – Whether to return the 1D power spectrum of the null test. Note that for null test, the result itself is the 1D power, so this is a dummy argument to keep consistency with cross and auto runs.

Returns:

arg_list – The list of arguments for the parallelisation.

Return type:

list

get_mock_instance_attr_dict(seed)[source]

Generate the attribute dictionary for the mock instance. It reads the attributes from the input power spectrum instance, and then sets a few mock-specific attributes specified in the input parameters of meer21cm.transfer.TransferFunction. The seed attribute is given at the input of this function.

Parameters:: seed (int) – The seed for the mock instance.
Returns:: attr_dict – The attribute dictionary for the mock instance.
Return type:: dict

run(seed_list, type='cross', return_power_3d=False, return_power_1d=False)[source]

Run the transfer function calculation.

Note that, run automatically uses a parallel pool to loop over the seed_list. If you believe the parallel behaviour is not as expected, you can manually extract the argument list and map the function yourself. For example: .. code-block:: python

>>> tf = TransferFunction(ps, N_fg=3)
>>> tf1d_arr = tf.run(range(10), type="auto")

is the same as: .. code-block:: python

>>> tf = TransferFunction(ps, N_fg=3)
>>> arg_list = tf.get_arg_list_for_parallel_auto(range(10))
>>> tf1d_arr = []
>>> with Pool(tf.num_process) as pool:
>>>     for result_i in pool.starmap(run_tf_calculation_auto, arg_list):
>>>         tf1d_arr.append(result_i)

Parameters:

seed_list (list) – The list of seeds for the mock instances.
type (str, default "cross") – The type of transfer function calculation to run. Can be “cross”, “auto”, or “null”.
return_power_3d (bool, default False) – Whether to return the 3D power spectrum of the transfer function.
return_power_1d (bool, default False) – Whether to return the 1D power spectrum of the transfer function.

Returns:

results_arr – The list of results for the transfer function calculation, with each element being a sublist for the result of one seed. The first element of each sublist is the transfer function. If return_power_3d is True, the second element is the 3D power spectrum of the original mock data, and the third element is the 3D power spectrum of the cleaned mock data. If return_power_1d is True, the next element is the 1D power spectrum of the original mock data, and the next element after that is the 1D power spectrum of the cleaned mock data.

Return type:

list

meer21cm.transfer.analytic_transfer_function(clean_mat_1, clean_mat_2=None)[source]

Calculate the analytic transfer function of a clean matrix. See Chen 2025 [1] for derivations.

For a foreground cleaning matrix \(R_{ab}\), the residual data vector for power spectrum estimation is \(r_{ija} = \sum_b R_{ab} m_{ijb}\), where i,j are the pixel indices, and a,b are the frequency indices.

Under the flat sky approximation, the signal loss, as well as the mode-mixing, is along the line-of-sight (k_para) direction. The unnormalised window function matrix is

\[H_{ab} = |\tilde{R}^1_{ab} (\tilde{R}^2_{ab})^*|_{\rm Re}\]

The corresponding signal loss is \(\sum_b H_{ab}\), and the analytical transfer function is the inverse of the signal loss.

After normalisation, the window function matrix is

\[W = {\rm diag}\Big(\sum_b H_{ab}\Big)^{-1} H\]

Parameters:

clean_mat_1 (np.ndarray) – The clean matrix that applies to the data vector.
clean_mat_2 (np.ndarray, optional) – The clean matrix that applies to the second data vector for cross-correlation. If not provided, it is assumed to be the same as clean_mat_1 for auto-correlation.

Returns:

transfer_func (np.ndarray) – The analytical transfer function.
Wab (np.ndarray) – The normalised window function matrix.

References

[1]

Chen, Z.,, “A quadratic estimator view of the transfer function correction in intensity mapping surveys”, https://ui.adsabs.harvard.edu/abs/2025MNRAS.542L…1C/abstract.

meer21cm.transfer.fft_matrix(mat, norm='backward')[source]

Perform the Fourier transform of a matrix.

\[\tilde{M} = \mathcal{F} M \mathcal{F}^{-1}\]

where \(\mathcal{F}\) is the Fourier transform matrix. See also meer21cm.util.dft_matrix().

Parameters:

mat (np.ndarray) – The matrix to be transformed.
norm (str, default "backward") – The normalization of the Fourier transform.

Returns:

mat_fft – The Fourier transformed matrix.

Return type:

np.ndarray

meer21cm.transfer.get_pca_matrix(map, N_fg, weights, mean_center_map)[source]

Get the PCA matrix that operates on the map data.

The eigendecomposition of the map data gives the eigenvectors \(v_{i}\), which can be used to form the source mixing matrix A so that

\[A = \{v_1, v_2, ..., v_{N_fg}\}\]

The PCA matrix is then

\[R = I - A A^T\]

The residual data vector after cleaning is then

\[r_{ija} = \sum_{b} R_{ab} m_{ijb}\]

where \(R_{ab}\) is the PCA matrix, and \(m_{ijb}\) is the uncleaned map data.

Parameters:

map (np.ndarray) – The map data.
N_fg (int) – The number of foreground components.
weights (np.ndarray) – The weights of the map data.
mean_center_map (bool) – Whether to mean center the map data.

meer21cm.transfer.run_null_test(mock_attr_dict, hi_map_rg, downres_factor_radial, downres_factor_transverse, weights_field_1, weights_grid_1, weights_field_2, weights_grid_2, k_sel_3d_to_1d, return_power_3d=False, unmask_during_mock=False, ra_range=None, dec_range=None, discrete_source_dndz=None)[source]

Run null test realisations by calculating the 1D cross-power spectrum of the mock galaxy x data map.

Parameters:

mock_attr_dict (dict) – The attribute dictionary to initialize the mock instance.
hi_map_rg (np.ndarray) – The HI data map.
downres_factor_radial (float) – The downres factor for the radial direction for gridding the map data.
downres_factor_transverse (float) – The downres factor for the transverse direction for gridding the map data.
weights_field_1 (np.ndarray) – The field-level weights for the field 1.
weights_grid_1 (np.ndarray) – The grid-level weights for the field 1.
weights_field_2 (np.ndarray) – The field-level weights for the field 2.
weights_grid_2 (np.ndarray) – The grid-level weights for the field 2.
k_sel_3d_to_1d (np.ndarray) – The weights for averaging the 3D power spectrum k-modes to 1D power spectrum.
return_power_3d (bool, default False) – Whether to return the 3D power spectrum of the mock results.
unmask_during_mock (bool, default False) – Whether to unmask the wcs during the mock simulation. If True, the survey area will be unmasked during the mock simulation, and then masked again after mock HI is generated. This is useful for mitigating edge effects and when box buffer is small. For a large patch with a wide z-range, it is recommended to set this to True, otherwise the mock HI map may contain nan and cause errors.
ra_range (tuple, default None) – The range of the RA in the sky map. Map and galaxy will be trimmed to this range.
dec_range (tuple, default None) – The range of the Dec in the sky map. Map and galaxy will be trimmed to this range.
discrete_source_dndz (list or function, default None) – The redshift distribution of the discrete tracer sources. Must be a function of redshift. Note that the overall number of discrete sources is set to ps.ra_gal.size or ps.num_discrete_source, so only the shape of the dndz is used, not the normalization. If a list is provided, the first element is the redshift array, and the second element is the dndz array, and a 1D linear interpolation is used to get the dndz function.

Returns:

result – The list of results, including the 1D cross-power spectrum of the mock galaxy x data map. If return_power_3d is True, the second element is the 3D power spectrum of the mock results.

Return type:

list

meer21cm.transfer.run_tf_calculation_auto(mock_attr_dict, downres_factor_radial, downres_factor_transverse, weights_field_1, weights_grid_1, k_sel_3d_to_1d, N_fg, pca_map_weights, mean_center_map, R_mat=None, uncleaned_data=0.0, return_power_3d=False, return_power_1d=False, unmask_during_mock=False, ra_range=None, dec_range=None)[source]

Run the transfer function calculation by calculating the ratio between the 1D power spectrum of cleaned mock HI data x original mock HI data to the 1D power spectrum of the original mock HI data x original mock HI data.

Parameters:

mock_attr_dict (dict) – The attribute dictionary to initialize the mock instance.
downres_factor_radial (float) – The downres factor for the radial direction for gridding the map data.
downres_factor_transverse (float) – The downres factor for the transverse direction for gridding the map data.
weights_field_1 (np.ndarray) – The field-level weights for the field 1.
weights_grid_1 (np.ndarray) – The grid-level weights for the field 1.
k_sel_3d_to_1d (np.ndarray) – The weights for averaging the 3D power spectrum k-modes to 1D power spectrum.
N_fg (int) – The number of foreground components to clean.
pca_map_weights (np.ndarray) – The weights for the mean centering and covariance matrix calculation during PCA.
mean_center_map (bool) – Whether to mean-center the map before PCA cleaning.
R_mat (np.ndarray, default None) – The PCA matrix to clean the map. If not provided, it will be calculated by injecting the original mock HI data into the data map. If provided, it will be used to clean the map, and no injection + PCA is performed.
uncleaned_data (np.ndarray, default 0.0) – The uncleaned data map.
return_power_3d (bool, default False) – Whether to return the 3D power spectrum of the mock results, including the uncleaned and cleaned data.
return_power_1d (bool, default False) – Whether to return the 1D power spectrum of the mock results, including the uncleaned and cleaned data.
unmask_during_mock (bool, default False) – Whether to unmask the wcs during the mock simulation. If True, the survey area will be unmasked during the mock simulation, and then masked again after mock HI is generated. This is useful for mitigating edge effects and when box buffer is small. For a large patch with a wide z-range, it is recommended to set this to True, otherwise the mock HI map may contain nan and cause errors.
ra_range (tuple, default None) – The range of the RA in the sky map. Map and galaxy will be trimmed to this range.
dec_range (tuple, default None) – The range of the Dec in the sky map. Map and galaxy will be trimmed to this range.

Returns:

result – The list of results, including the transfer function, the 3D power spectrum of the uncleaned and cleaned data, and the 1D power spectrum of the uncleaned and cleaned data.

Return type:

list

meer21cm.transfer.run_tf_calculation_cross(mock_attr_dict, downres_factor_radial, downres_factor_transverse, weights_field_1, weights_grid_1, weights_field_2, weights_grid_2, k_sel_3d_to_1d, N_fg, pca_map_weights, mean_center_map, R_mat=None, uncleaned_data=0.0, return_power_3d=False, return_power_1d=False, unmask_during_mock=False, ra_range=None, dec_range=None, discrete_source_dndz=None)[source]

Run the transfer function calculation by calculating the ratio of the 1D cross-power spectrum of the cleaned mock HI data x mock galaxy data to the 1D cross-power spectrum of the original mock HI data x mock galaxy data.

Parameters:

mock_attr_dict (dict) – The attribute dictionary to initialize the mock instance.
downres_factor_radial (float) – The downres factor for the radial direction for gridding the map data.
downres_factor_transverse (float) – The downres factor for the transverse direction for gridding the map data.
weights_field_1 (np.ndarray) – The field-level weights for the field 1.
weights_grid_1 (np.ndarray) – The grid-level weights for the field 1.
weights_field_2 (np.ndarray) – The field-level weights for the field 2.
weights_grid_2 (np.ndarray) – The grid-level weights for the field 2.
k_sel_3d_to_1d (np.ndarray) – The weights for averaging the 3D power spectrum k-modes to 1D power spectrum.
N_fg (int) – The number of foreground components to clean.
pca_map_weights (np.ndarray) – The weights for the mean centering and covariance matrix calculation during PCA.
mean_center_map (bool) – Whether to mean-center the map before PCA cleaning.
R_mat (np.ndarray, default None) – The PCA matrix to clean the map. If not provided, it will be calculated by injecting the original mock HI data into the data map. If provided, it will be used to clean the map, and no injection + PCA is performed.
uncleaned_data (np.ndarray, default 0.0) – The uncleaned data map.
return_power_3d (bool, default False) – Whether to return the 3D power spectrum of the mock results, including the uncleaned and cleaned data.
return_power_1d (bool, default False) – Whether to return the 1D power spectrum of the mock results, including the uncleaned and cleaned data.
unmask_during_mock (bool, default False) – Whether to unmask the wcs during the mock simulation. If True, the survey area will be unmasked during the mock simulation, and then masked again after mock HI is generated. This is useful for mitigating edge effects and when box buffer is small. For a large patch with a wide z-range, it is recommended to set this to True, otherwise the mock HI map may contain nan and cause errors.
ra_range (tuple, default None) – The range of the RA in the sky map. Map and galaxy will be trimmed to this range.
dec_range (tuple, default None) – The range of the Dec in the sky map. Map and galaxy will be trimmed to this range.
discrete_source_dndz (list or function, default None) – The redshift distribution of the discrete tracer sources. Must be a function of redshift. Note that the overall number of discrete sources is set to ps.ra_gal.size or ps.num_discrete_source, so only the shape of the dndz is used, not the normalization. If a list is provided, the first element is the redshift array, and the second element is the dndz array, and a 1D linear interpolation is used to get the dndz function.

Returns:

result – The list of results, including the transfer function, the 3D power spectrum of the uncleaned and cleaned data, and the 1D power spectrum of the uncleaned and cleaned data.

Return type:

list

meer21cm.util module

class meer21cm.util.HiddenPrints[source]

Bases: object

A context manager that suppresses print statements. This is useful for suppressing print statements in the middle of a calculation.

class meer21cm.util.Obuljen18(**model_parameters)[source]

Bases: HODBulk

A halomod.hod.HODBulk object for the HI-halo mass relation reported in Obuljen et al. (2018) [1].

Attributes:

mmin: Defines a reasonable minimum mass to set for this HOD to converge when integrated.

Methods

`central_occupation`(m)	The occupation function of the central component.
`cs_pairs`(m)	The average amount of the tracer coupled with itself in haloes of mass m, <T_s T_c>.
`get_models`()	Get a dictionary of all implemented models for this component.
`nc`(m)	Density of Central Tracer.
`ns`(m)	Density of Satellite Tracer.
`satellite_occupation`(m)	The occupation function of the satellite (or profile-dependent) component.
`sigma_central`(m)	The standard deviation of the central tracer amount in haloes of mass m.
`sigma_satellite`(m)	The standard deviation of the satellite tracer amount in haloes of mass m.
`ss_pairs`(m)	The average amount of the tracer coupled with itself in haloes of mass m, <T_s T_s>.
`total_occupation`(m)	The total (average) occupation of the halo.
`total_pair_function`(m)	The total weight of the occupation paired with itself.
`unit_conversion`(cosmo, z)	A factor to convert the total occupation to a desired unit.

sigma_satellite(m)[source]: The standard deviation of the satellite tracer amount in haloes of mass m.

meer21cm.util.angle_in_range(alpha, lower, upper)[source]: Determines whether the angle alpha is in a range. All inputs in degree. Stolen from [this url](https://stackoverflow.com/questions/66799475/how-to-elegantly-find-if-an-angle-is-between-a-range).

meer21cm.util.busy_function_simple(xarr, par_a, par_b, par_c, width)[source]

The simplified busy function that assumes mirror symmetry around x=0 [1].

\[B_2(x) = \frac{a}{2} \times ({\rm erf}[b(w^2-x^2)]+1) \times (cx^2+1)\]

Parameters:

xarr (float array.) – the input x values
par_a (float.) – amplitude parameter
par_b (float.) – b parameter that controls the sharpness of the double peaks
par_c (float.) – c parameter that controls the height of the double peaks
width (float.) – the width of the profile

Returns:

b2x – the busy function values at xarr

Return type:

float array.

References

[1]

Westmeier, T. et al., “The busy function: a new analytic function for describing the integrated 21-cm spectral profile of galaxies”, https://ui.adsabs.harvard.edu/abs/arXiv:1311.5308 .

meer21cm.util.cal_himf(x, mmin, cosmo, mmax=11, integrate_step=500)[source]

Calculate the integrated quantity related to the HIMF.

Parameters:

x (list of float.) – need to be [phi_s,log10(m_s),alpha_s]
mmin (float.) – The minimum mass to integrate from in log10
cosmo (an astropy.cosmology.Cosmology object.) – The cosmology object to calculate critical density
mmax (Optional float, default 11.) – The maximum mass to integrate to in log10.
integrate_step (Optional int, default 500.) – The number of steps to integrate the HIMF.

Returns:

nhi (float.) – The number density of HI galaxies, in the units of phi_s * dex
omegahi (float.) – The HI density over the critical density of the present day (assuming the recommended units for x are used)
psn (float.) – The shot noise in the units of Mpc:sup:3 (assuming the recommended units for x are used)

meer21cm.util.center_to_edges(arr)[source]: Extend a linear spaced monotonic array so that the original array is the middle point of the output array.

meer21cm.util.check_unit_equiv(u1, u2)[source]

Check if two units are equivelant

Parameters:

u1 (astropy.units.Unit object.) – The first input unit.
u2 (astropy.units.Unit object.) – The second input unit.

Returns:

result – Whether they are the same unit.

Return type:

bool.

meer21cm.util.coeff_hi_density_to_temp(z=0, cosmo=FlatLambdaCDM(name='Planck18', H0=<Quantity 67.66 km / (Mpc s)>, Om0=0.30966, Tcmb0=<Quantity 2.7255 K>, Neff=3.046, m_nu=<Quantity [0., 0., 0.06] eV>, Ob0=0.04897))[source]

The conversion coefficient \(C_{\rm HI}\) so that

\[\bar{T}_{\rm HI} = C_{\rm HI} \rho_{\rm HI}\]

Parameters:

z (float, defulat 0.0.) – The redshift
cosmo (cosmology, default Planck18.) – The cosmology used

Returns:

C_HI – The coefficient

Return type:

quantity.

meer21cm.util.convert_hpmap_in_jy_to_temp(hp_map, freq)[source]

meer21cm.util.create_udres_wproj(wproj, udres_scale)[source]

Take an input wcs and create a wcs object that upgrades/downgrades the resolution. The ratio between the output resolution and input is determined by udres_scale.

Parameters:

wproj (astropy.wcs.WCS object.) – The input wcs.
udres_scale (float) – The ratio between input and output resolution.

Returns:

w – The output wcs.

Return type:

astropy.wcs.WCS object.

meer21cm.util.create_wcs(ra_cr, dec_cr, ngrid, resol, crpix=None, ctype=['RA---ZEA', 'DEC--ZEA'])[source]

Create a wcs object that can be used to map a index array to sky cooridnates, for a given central coordinate and grid dimensions.

Parameters:

ra_cr (float.) – The coordinate of the critical pixel in RA in degree.
dec_cr (float.) – The coordinate of the critical pixel in Dec in degree.
ngrid (array-like of 2 elements.) – The number of pixels on the two axis.
resol (array-like of 2 elements.) – The resolution along the two axis in degree.
crpix (array-like of 2 elements, default None.) – The index of the critical pixel. Default is the center of the array.
ctype (list of two strings, default ['RA---ZEA', 'DEC--ZEA'].) – The projection type.

Returns:

w – The output wcs.

Return type:

astropy.wcs.WCS object.

meer21cm.util.create_wcs_with_range(ra_range, dec_range, resol=[0.3, 0.3], buffer=[1.2, 1.2], ctype=['RA---ZEA', 'DEC--ZEA'])[source]

Create a wcs object that can be used to map a index array to sky cooridnates, based on an approximate survey area specified by ra_range and dec_range.

A naive number of pixels will be calculated for both directions based on the difference between the lower and upper limit of the range, divided by the input resolution. It is then multiplied by the buffer because usually you need more than that due to curved sky.

Parameters:

ra_range (array-like of 2 elements.) – The lower and upper limit of the RA range in degree.
dec_range (array-like of 2 elements.) – The lower and upper limit of the Dec range in degree.
resol (array-like of 2 elements, default [0.3,0.3].) – The resolution along the two axis in degree.
buffer (array-like of 2 elements, default [1.2,1.2].) – The proportional increase to the range.
ctype (list of two strings, default ['RA---ZEA', 'DEC--ZEA'].) – The projection type.

Returns:

w (astropy.wcs.WCS object.) – The output wcs.
num_pix_x (int.) – Number of pixels on the first axis.
num_pix_y (int.) – Number of pixels on the second axis.

meer21cm.util.cumu_nhi_from_himf(m, mmin, x)[source]

The integrated source number density from HIMF.

Parameters:

m (float array.) – The higher end of integration in log10
mmin (float.) – The minimum mass to integrate from in log10
x (list of float.) – need to be [phi_s,log10(m_s),alpha_s]

Returns:

nhi – The integrated number density of HI galaxies, in the units of phi_s * dex

Return type:

float array.

meer21cm.util.dft_matrix(N, norm='backward')[source]

Generate the DFT matrix for a given N. The default is the backward normalization, which is the same as np.fft.fft.

Parameters:

N (int) – The size of the matrix.
norm (str, default 'backward') – The normalization of the DFT matrix.

Returns:

dft_mat – The DFT matrix.

Return type:

np.ndarray

meer21cm.util.find_block_id(filename)[source]

Find the MeerKAT data block id from the filename. MeerKAT data blocks are identified by a 10-digit number. Use vfind_id to find the block id for an array of filenames.

Parameters:: filename (str) – The filename to find the block id from.
Returns:: block_id – The block id.
Return type:: str

meer21cm.util.find_ch_id(nu_inp, nu_ch)[source]

For which channel does the input frequency fall into. The channel ids are zero-indexed. Input frequencies outside the frequency range are assigned \(N_{\rm ch}\) (note the last channel id is \(N_{\rm ch}-1\))

Parameters:

nu_inp (array.) – The input frequencies
nu_ch (array.) – Must be monotonically increasing. The centre frequencies of each channel.

Returns:

which_ch – The ch ids.

Return type:

int array.

meer21cm.util.find_indx_for_subarr(subarr, arr)[source]

Find the indices of the elements of an array in another array.

Parameters:

subarr (numpy array.) – The sub-array to search for. Elements can be repeated.
arr (numpy array.) – The larger array that contains all elements of subarr.

Returns:

indices – The position of the elements of subarr in arr.

Return type:

int array.

meer21cm.util.find_property_with_tags(obj)[source]: Retrieve a dictionary for all the properties of a class that has tags. The keys of the dictionary are the property names and the values are the tags of each property.

meer21cm.util.freq_to_redshift(freq)[source]

Convert frequency of 21cm emission to redshift

Parameters:: freq (float)
Returns:: redshift
Return type:: float

meer21cm.util.get_ang_between_coord(ra1, dec1, ra2, dec2, unit='deg')[source]

Calculate the angle between two points on the sphere.

Parameters:

ra1 (float array) – The RA of the first point
dec1 (float array) – The Dec of the first point
ra2 (float array) – The RA of the second point
dec2 (float array) – The Dec of the second point
unit (str, default 'deg'.)

Returns:

result – The angle in the specified unit.

Return type:

float array.

meer21cm.util.get_default_args(func)[source]

meer21cm.util.get_nd_slicer(ndim=3)[source]

Get a list of slice objects that can be used to slice an ndarray. For example, if you have a list of k-vectors, then you can this to match the vectors for array index propagation.

Parameters:: ndim (int, default 3.) – The number of dimensions of the array to slice.
Returns:: out
Return type:: list of slice objects.

meer21cm.util.get_wcs_coor(wcs, xx, yy, ang_unit='deg')[source]

Retrieve RA and Dec coordinates of pixels in the WCS.

Parameters:

wproj (astropy.wcs.WCS object.) – The two-dimensional wcs object for the map.
xx (array) – The first indx of the pixel.
yy (array) – The second indx of the pixel.

Returns:

ra (array.) – The RA of the pixels.
dec (array.) – The Dec of the pixels.

meer21cm.util.healpix_to_wcs(hp_map, xx, yy, wcs)[source]: Project the healpix map onto a 2-D wcs grid. Map has to be in temperature unit.

meer21cm.util.himf(m, phi_s, m_s, alpha_s)[source]

Analytical HIMF function (or any other Schechter function).

\[\phi = {\rm ln}(10) \times \phi_* \times \bigg(\frac{m}{m_*} \bigg)^{\alpha_*+1}\times e^{-m/m_*}\]

While the units are arbitrary, it is recommended that phi_s is in the unit of Mpc:sup:-3`dex:sup:-1`, m_s is in the unit of log10 M_sun, and alpha_s has no unit.

Parameters:

m (float array.) – mass in log10.
phi_s (float.) – HIMF amplitude
m_s (float.) – knee mass in log10.
alpha_s (float.) – slope

Returns:

out – The HIMF values at m

Return type:

float array.

meer21cm.util.himf_pars_jones18(h_70)[source]

The HIMF parameters measured in Jones+18 [1].

Parameters:

h_70 (float.) – The Hubble parameter over 70 km/s/Mpc.

Returns:

phi_star (float.) – The amplitude of HIMF in Mpc-3 dex-1.
m_star (float.) – The knee mass of HIMF in log10 solar mass.
alpha (float.) – The slope of HIMF.

References

[1]

Jones, M. et al., “The ALFALFA HI mass function: A dichotomy in the low-mass slope and a locally suppressed ‘knee’ mass”, https://ui.adsabs.harvard.edu/abs/arXiv:1802.00053.

meer21cm.util.hod_obuljen18(logmh, m0h=9.52, mminh=11.27, alpha=0.44, cosmo=FlatLambdaCDM(name='Planck18', H0=<Quantity 67.66 km / (Mpc s)>, Om0=0.30966, Tcmb0=<Quantity 2.7255 K>, Neff=3.046, m_nu=<Quantity [0., 0., 0.06] eV>, Ob0=0.04897), input_has_h=True, output_has_h=False)[source]

HI-halo mass relation reported in Obuljen et al. (2018) [1]. The default settings assume the mass is in M_sun/h, and returns mass without h unit. Note that, if you want input without h unit, you must change the m0h and mminh manually as well. The HI-halo mass relation follows

\[M_{\rm HI} (M_h) = M_0 (M_h/M_{\rm min})^\alpha \, {\rm exp}[-M_{\rm min}/M_h],\]

Parameters:

logmh (float.) – input halo mass in log10
m0h (optional, default 9.52.) – The \(M_0\) parameter in the HOD in log10
mminh (optional, default 11.27.) – The \(M_{min}\) parameter in the HOD in log10
alpha (optional, default 0.44.) – The \(\alpha\) parameter in the HOD.
cosmo (optional, default Planck18.) – The default cosmology, used only for h unit.
input_has_h (optional, default True.) – Whether the input mass has h unit.
output_has_h (optional, default False.) – Whether the output mass has h unit.

Return type:

The HI mass for the input halo mass.

References

[1]

Obuljen, A. et al., “The H I content of dark matter haloes at z ≈ 0 from ALFALFA “, https://ui.adsabs.harvard.edu/abs/2019MNRAS.486.5124O.

meer21cm.util.inv_dft_matrix(N, norm='backward')[source]

Generate the inverse DFT matrix for a given N. The default is the backward normalization, which is the same as np.fft.fft.

Parameters:

N (int) – The size of the matrix.
norm (str, default "backward") – The normalization of the DFT matrix.

Returns:

inv_dft_mat – The inverse DFT matrix.

Return type:

np.ndarray

meer21cm.util.jy_to_kelvin(val, omega, freq)[source]

convert Jy/beam to brightness temperature in Kelvin.

Parameters:

val (numpy array.) – The input value(s) in Jy/beam or Jy/pix
omega (float.) – beam or pixel area in Steradian.
freq (float.) – the frequency for conversion in Hz.

Returns:

result – The brightness temperature in Kelvin.

Return type:

float array.

meer21cm.util.legendre_polynomial_with_factor(ell: int, return_coeff: bool = True)[source]

Return the Legendre polynomial with the given ell, with the (2*ell+1) factor.

Parameters:

ell (int) – The Legendre polynomial order.
return_coeff (bool, default True) – Whether to return the coefficients or the polynomial or the callable function.

Returns:

coeff (np.ndarray) – The coefficients of the Legendre polynomial.
poly (np.poly1d) – The Legendre polynomial function.

meer21cm.util.mean_center_signal(signal, weights=None, los_axis=-1)[source]

Mean-center the signal along the LOS (spectral) axis.

Supports 3D WCS (nx, ny, n_ch) and 2D HEALPix (n_pix, n_ch) layouts; los_axis picks which axis is frequency (default -1).

Parameters:

signal (array.) – The input signal
weights (array.) – The weights of each element in the signal.
los_axis (int, default -1.) – Index of the spectral axis (-1 is the last axis; e.g. (n_pix, n_ch)).

Returns:

centered_signal – The mean-centered signal.

Return type:

array.

meer21cm.util.omega_hi_to_average_temp(omega_hi, z=0, cosmo=FlatLambdaCDM(name='Planck18', H0=<Quantity 67.66 km / (Mpc s)>, Om0=0.30966, Tcmb0=<Quantity 2.7255 K>, Neff=3.046, m_nu=<Quantity [0., 0., 0.06] eV>, Ob0=0.04897))[source]

The average HI brightness temperature from given HI density.

Parameters:

omega_hi (float) – The HI density relative to the z=0 critical density.
z (float, defulat 0.0.) – The redshift
cosmo (cosmology, default Planck18.) – The cosmology used

Returns:

t_bar – The average HI temperature in Kelvin

Return type:

float.

meer21cm.util.pca_clean(signal, N_fg, weights=None, return_analysis=False, mean_center=False, los_axis=-1, return_A=False, mean_center_weights=None, ignore_nan=False, covariance=None)[source]

Performs PCA cleaning of map data (3D WCS cube or 2D HEALPix (n_pix, n_ch)). If mean_center is set to True, then the input signal is first mean-centered so that

\[\vec{d}^{\rm mc} = \vec{d} - \langle \vec{d} \rangle,\]

where the mean of the data vector at a specific channel i is

\[\langle \vec{d} \rangle_i = \frac{\sum_a w_{ia} d_{ia}}{\sum_a w_{ia}},\]

where a loops over each sampling in channel i and \(w\) is the weight of each element.

A covariance matrix is then calculated on the mean-centered data

\[C_{ij} = \frac{\sum_{a,b }w_{ia} d^{\rm mc}_{ia} w_{jb} d^{\rm mc}_{jb}}{\sum_{a,b }w_{ia}w_{jb}},\]

over which the eigenvalue decomposition is performed. See Section 4.3. of MeerKLASS Collaboration 2024 [1] for more details.

Note that, in the rigorous definition, the data vector should be zero-meaned, and the weights used to calculate the mean and the covariance should be the same. However, in practice, many people don’t remove the mean of the data. Some use one type of weights (often just binary masks) for mean calculation, and then use different weights for covariance calculation. While it is not encouraged, that flexibility is allowed in the function, by setting a different weight mean_center_weights.

If there are frequency gaps in the data, you can set ignore_nan=True to ignore these channels, and treat the rest of the data as a continuous spectrum.

Parameters:

signal (array) – Map cube: either 3D WCS (nx, ny, n_ch) or 2D HEALPix (n_pix, n_ch) with the spectral axis given by los_axis (default last axis for both).
N_fg (int.) – Number of modes to be removed.
weights (array, default None) – The weights of each element in the signal. Default will set uniform weights for each element.
return_analysis (bool, default False) – If True, instead of residual maps the function will return eigenanalysis quantities.
mean_center (bool, default False) – Whether to mean-center the input data vector
los_axis (int, default -1.) – Index of the spectral axis (-1 is the last axis; e.g. (n_pix, n_ch)).
return_A (bool, default False.) – Whether to return the mixing matrix A.
mean_center_weights (array, default None.) – The weights of each element for mean calculation. Default follows the weights argument.
ignore_nan (bool, default False.) – Whether to ignore NaN values in the data. If True, channels with NaN values are skipped.
covariance (array, default None.) – The covariance matrix of the input signal. If not provided, it will be calculated from the input signal.

Returns:

residual (array.) – The residual after PCA cleaning.
A (array, if return_A=True.) – The mixing matrix.
covariance (array, if return_analysis=True.) – The covariance matrix of the input signal.
eignumb (array, if return_analysis=True.) – The number indexing of the eigenvalues starting from 1.
eigenval (array, if return_analysis=True.) – The eigenvalues of the covariance matrix.
V (array, if return_analysis=True.) – The eigenvectors of the covariance matrix.

References

[1]

MeerKLASS collab, “MeerKLASS L-band deep-field intensity maps: entering the HI dominated regime”, https://arxiv.org/abs/2407.21626.

meer21cm.util.ra_array_crosses_zero(ra_deg, atol=1e-09)[source]

Return True if the smallest RA arc containing all values passes through 0° / 360°.

Longitudes are taken modulo 360°. For example, an array mixing values near 350° and near 10° returns True (the short arc bridges the branch cut). An array confined to a single segment that does not wrap, such as 10°–30°, returns False.

Parameters:

ra_deg (array_like) – Right ascensions in degrees.
atol (float) – Numerical tolerance in degrees (for comparing arc lengths).

Return type:

bool

meer21cm.util.radec_to_indx(ra_arr, dec_arr, wproj, to_int=True, ang_unit='deg')[source]

meer21cm.util.random_sample_indx(tot_len, num_sub_sample, seed=None)[source]

Generate a random sub-sample indices.

Parameters:

tot_len (int.) – The total number of data points to sample from.
num_sub_sample (int.) – Number of sub-samples.
seed (int, default None.) – The seed for the random number generator.

Returns:

sub_indx – The sub-sample indices.

Return type:

array.

meer21cm.util.read_healpix_fits(file)[source]

meer21cm.util.real_dtype_from_array(array, default=<class 'numpy.float64'>)[source]

Return a real floating dtype compatible with an input array.

Parameters:

array (np.ndarray) – The input array.
default (np.dtype, default np.float64.) – The default dtype to return if the input array is not complex or floating.

Returns:

dtype – The real floating dtype compatible with the input array.

Return type:

np.dtype.

meer21cm.util.rebin_spectrum(spectrum, rebin_width=3, mode='avg')[source]

meer21cm.util.redshift_to_freq(redshift)[source]

Convert redshift to frequency

Parameters:: redshift (float)
Returns:: freq
Return type:: float

meer21cm.util.sample_from_dist(func, xmin, xmax, size=1, cdf=False, seed=None)[source]

Sample from custom distribution.

Parameters:

func (distribution function.) – The probability distribution function (cdf=False) or the cumulative distribution function (cdf=True).
xmin (float.) – The minimum value to sample from
xmax (float.) – The maximum value to sample from
size (int or list of int, default 1.) – The size of the sample array
cdf (bool, default False.) – Wheter PDF or CDF is used.
seed (int, default None.) – Seed for the random number generator. Fix for reproducible samples.

Returns:

sample – The random sample following the input distribution function.

Return type:

float array.

meer21cm.util.sample_map_from_highres(map_highres, ra_map, dec_map, wproj_lowres, num_pix_x, num_pix_y, average=True)[source]

grid a highres map into lowres.

Parameters:

map_highres (array) – The input high-res map.
ra_map (array) – The RA coordinates of the input map pixels.
dec_map (array) – The Dec coordinates of the input map pixels.
wproj_lowres (astropy.wcs.WCS object.) – The wcs object for the low-res map.
num_pix_x (int) – The number of pixels along the first axis for the low-res map.
num_pix_y (int) – The number of pixels along the second axis for the low-res map.
average (bool, default True.) – Whether the sampling is an average of the high-res pixels (True) or the sum (False).

meer21cm.util.super_sample_array(arr_in, super_factor)[source]

Super-sample an array along each direction by a factor

Parameters:

arr_in (array.) – The input array to sample from.
super_factor (list of int.) – super sample factor along each axis.

Returns:

arr_out – The super-sampled array.

Return type:

array.

meer21cm.util.tagging(*tags)[source]

A decorator that does one simple thing: Adding tags to functions.

For example, you can add tags when defining this function

from meer21cm.util import tagging
@tagging('test')
def foo(x):
   return x
print(foo.tags)

You will find that foo.tags is ('test',). meer21cm uses this function to keep track of parameter dependencies of functions in classes.

meer21cm.util.tightest_ra_interval(ra_deg)[source]

Smallest RA interval (lower, upper) in degrees that contains all values, using the same wrap convention as angle_in_range() (lower > upper when the arc crosses 0°).

Parameters:: ra_deg (array_like) – Right ascensions in degrees.
Returns:: lower, upper – Bounds in degrees. If all values coincide, both are that longitude. If the points span the full circle, returns (0.0, 360.0).
Return type:: float
Raises:: ValueError – If ra_deg is empty.

meer21cm.util.tully_fisher(xarr, slope, zero_point, inv=False)[source]

Tully-Fisher relation.

Note that, regardless of inv, the slope and zero_point always refer to the Tully-Fisher relation and not the inverse. For example, zero_point is always in the unit of log10 mass.

Parameters:

xarr (float array.) – input velocity if inv=False and mass if inv=True.
slope (float.) – the slope of Tully-Fisher relation.
zero_point (float.) – the intercept of Tully-Fisher relation
inv (bool, default False.) – if True, calculate velocity based on input mass.

Returns:

out – The output mass if inv=False and velocity if inv=True.

Return type:

float array.

meer21cm.util.weighted_covariance(signal, weights, renorm=True)[source]

Calculate the weighted covariance matrix of the signal.

The signal is assumed to be mean-centered, with the shape of (num_ch, num_pix).

Parameters:

signal (array.) – The input signal
weights (array.) – The weights of each element in the signal.
renorm (bool, default True.) – Whether to renormalize the covariance matrix over the sum of weights.

Returns:

covariance – The covariance matrix of the input signal.

Return type:

array.

meer21cm.util.which_ra_range_is_tighter(range_a, range_b, atol=1e-09)[source]

Compare two RA intervals given as (lower, upper) in degrees, using the same wrap convention as angle_in_range() (lower may be greater than upper).

Exactly one range must enclose the other (as sets of longitudes on the circle), or the two ranges must be equivalent; otherwise a ValueError is raised.

Parameters:

range_a (tuple of float) – (ra_lower_deg, ra_upper_deg).
range_b (tuple of float) – (ra_lower_deg, ra_upper_deg).
atol (float) – Numerical tolerance in degrees (full-sky detection and span comparisons).

Returns:

-1 if range_a is strictly tighter (proper subset of range_b), 1 if range_b is strictly tighter, 0 if the two ranges describe the same coverage.

Return type:

int

Raises:

ValueError – If neither range encloses the other (including partial overlap without nesting).