Main modules

Angular coordinates

Methods to convert data between physical (cMpc) coordinates and observational (angular-frequency) coordinates.

tools21cm.angular_coordinates.angular_slice_to_physical(input_slice, z, slice_size_deg, output_cell_size, output_size_mpc, order=0, prefilter=True)

Interpolate a slice in angular coordinates to physical

Parameters:

• input_slice (numpy array) – the 2D slice in observational coordinates

• z (float) – the redshift of the input slice

• slice_size_deg (float) – the size of the input slice in deg

• output_cell_size (float) – the output cell size in cMpc

• output_size_mpc (float) – the output size in mpc

Returns:

(physical_slice, size_mpc)

tools21cm.angular_coordinates.bin_lightcone_in_frequency(lightcone, z_low, box_size_mpc, dnu)

Bin a lightcone in frequency bins.

Parameters:

• lightcone (numpy array) – the lightcone in length units

• z_low (float) – the lowest redshift of the lightcone

• box_size_mpc (float) – the side of the lightcone in Mpc

• dnu (float) – the width of the frequency bins in MHz

Returns:

• The lightcone, binned in frequencies with high frequencies first

• The frequencies along the line of sight in MHz

tools21cm.angular_coordinates.bin_lightcone_in_mpc(lightcone, frequencies, cell_size_mpc)

Bin a lightcone in Mpc slices along the LoS

tools21cm.angular_coordinates.observational_lightcone_to_physical(observational_lightcone, input_freqs, input_dtheta, verbose=True, order=2)

Interpolate a lightcone volume measured in observational (angle/frequency) units into physical (length) units. The output resolution will be set to the coarest one, as determined either by the angular or the frequency resolution. The lightcone must have the LoS as the last index, with frequencies decreasing along the LoS.

Parameters:

• observational_lightcone (numpy array) – the input lightcone volume

• input_freqs (numpy array) – the frequency in MHz of each slice along the line of sight of the input

• input_dheta (float) – the angular size of a cell in arcmin

• verbose (bool) – show progress bar

• order (int) – The order of the spline interpolation, default is 2. The order has to be in the range 0-5. Use order=0 for ionization fraction data.

Returns:

• The output volume

• The redshifts along the LoS of the output

• The output cell size in Mpc

tools21cm.angular_coordinates.padding_lightcone(lc, padded_n=None, mode='wrap', verbose=True)

Pad lightcone in the field of view direction.

Parameters:

• lc (numpy array) – the light-cone in physical coordinates with third axis as the light-of-sight.

• padded_n (int) – number of cells to pad. The default value is half of the number of cells.

• mode (str) – mode used for padding (see the documentation of numpy.pad).

• verbose (bool) – show progress bar.

Returns:

padded lightcone

tools21cm.angular_coordinates.physical_lightcone_to_observational(physical_lightcone, input_z_low, output_dnu, output_dtheta, input_box_size_mpc=None, verbose=True, order=2, mode='pad')

Interpolate a lightcone volume from physical (length) units to observational (angle/frequency) units.

Parameters:

• physical_lightcone (numpy.ndarray) – The input lightcone volume in physical coordinates with shape (N, N, M) where N is the number of cells per side in angular direction and M is the number of cells along the line-of-sight

• input_z_low (float or array-like) – The lowest redshift of the input lightcone. If array-like, should contain redshifts for each slice along the line-of-sight

• output_dnu (float) – The frequency resolution of the output volume in MHz

• output_dtheta (float) – The angular resolution of the output in arcmin

• input_box_size_mpc (float, optional) – The size of the input field of view in Mpc. If None (default), uses the value from conv.LB

• verbose (bool, optional) – Whether to show progress messages (default: True)

• order (int, optional) – The order of the spline interpolation (0-5). Default is 2. Use order=0 for ionization fraction data.

• mode (str, optional) – How to handle the field of view at different redshifts:

  • ‘pad’ (default): Fix FoV at lowest redshift (largest angular size) and pad higher redshift slices

  • ‘crop’: Fix FoV at highest redshift (smallest angular size) and crop lower redshift slices

  • ‘full’, ‘extend’: Same as ‘pad’

  • ‘valid’: Same as ‘crop’

Returns:

tuple –

output_volumenumpy.ndarray

The transformed lightcone in observational coordinates with shape (N_theta, N_theta, N_nu)

output_freqsnumpy.ndarray

Array of output frequencies in MHz corresponding to the slices along the line-of-sight

tools21cm.angular_coordinates.physical_slice_to_angular(input_slice, z, slice_size_mpc, fov_deg, dtheta, order=0)

Interpolate a slice in physical coordinates to angular coordinates.

Parameters:

• input_slice (numpy array) – the 2D slice in physical coordinates

• z (float) – the redshift of the input slice

• slice_size_Mpc (float) – the size of the input slice in cMpc

• fov_deg (float) – the field-of-view in degrees. The output will be padded to match this size

• dtheta (float) – the target resolution in arcmin

Returns:

(angular_slice, size_deg)

tools21cm.angular_coordinates.resample_slice(input_slice, n_output_cells, order=0, prefilter=True)

Resample a 2D slice to new dimensions.

Parameters:

• input_slice (numpy array) – the input slice

• n_output_cells (int) – the number of output cells

Returns:

output slice

Useful Conversion Functions

tools21cm.astro_tools.FoV_cMpc2(FoV, redshift, cosmo=None)

Convert field of view from square degrees to comoving megaparsecs squared (cMpc^2).

tools21cm.astro_tools.Muv_to_L1500(M_AB)

Convert absolute UV magnitude to UV luminosity at 1500 Å (erg/s/Hz).

tools21cm.astro_tools.Muv_to_muv(Muv, z, cosmo=None)

Convert absolute magnitude to apparent magnitude in the UV band.

tools21cm.astro_tools.mag_to_flux(mag)

Convert apparent magnitude to flux density in CGS units (erg/s/cm^2/Hz).

tools21cm.astro_tools.muv_to_Muv(muv, z, cosmo=None)

Convert apparent magnitude to absolute magnitude in the UV band.

Bubble Statistics

Methods to calcluate the sizes of the regions of interest and estimate the size distributions.

tools21cm.bubble_stats.averaged_mfp(data, xth=0.5, boxsize=None, iterations=10000000, rays_per_point=10, verbose=True, upper_lim=False, bins=None, r_min=None, r_max=None, n_jobs=-1)

Determines the sizes using the averaged Mean-Free-Path (MFP) approach.

Parameters:

• input (ndarray) – 2D/3D array of ionization fraction/brightness temperature.

• xth (float) – The threshold value (Default: 0.5).

• boxsize (float) – The boxsize in cMpc can be given (Default: conv.LB).

• iterations (int) – Number of iterations (Default: 10_000_000).

• iterations (int) – Number of iterations (Default: 10).

• verbose (bool) – It prints the progress of the program (Default: True).

• upper_lim (bool) – It decides if the threshold is the upper limit or the lower limit (Default: False).

• bins (int) – Give number of bins or an array of sizes to re-bin into (Default: None).

• r_min (float) – Minimum size after rebinning (Default: None).

• r_max (float) – Maximum size after rebinning (Default: None).

Returns:

• r (ndarray) – sizes of the regions

• dn (ndarray) – probability of finding the corresponding size

tools21cm.bubble_stats.dist_from_volumes(sizes, resolution=1.0, bins=100, null_factor=None)

Volume distribution and effective radius distribution.

Parameters:

• sizes (list) – List of volumes in pixel**3 units.

• resolution (float) – Distance between two pixels in cMpc (Default: 1).

• bins (int) – Number of bins of volumes in log space.

Returns:

• v (ndarray) – volumes of the regions

• Pv (ndarray) – probability of finding the corresponding volume

• r (ndarray) – effective radii of the regions

• Pr (ndarray) – probability of finding the corresponding effective radius

tools21cm.bubble_stats.fof(data, xth=0.5, connectivity=1)

Determines the sizes using the friends-of-friends approach. It assumes the length of the grid as the linking length.

Parameters:

• data (ndarray) – The array containing the input data

• xth (float) – The threshold value (Default: 0.5)

Returns:

• map (ndarray) – array with each regions of interest label

• sizes (list) – all the sizes

tools21cm.bubble_stats.granulometry_bsd(data, xth=0.5, boxsize=None, verbose=True, upper_lim=False, sampling=2, log_bins=None, n_jobs=1)

Determined the sizes using the Granulometry (Gran) approach. It is based on Kakiichi et al. (2017)

Parameters:

• input (ndarray) – 2D/3D array of ionization fraction/brightness temperature.

• xth (float) – The threshold value (Default: 0.5).

• boxsize (float) – The boxsize in cMpc can be given (Default: conv.LB).

• verbose (bool) – It prints the progress of the program (Default: True).

• upper_lim (bool) – It decides if the threshold is the upper limit or the lower limit (Default: False).

• sampling (int) – Give the resolution of the radii in the pixel units (Default: 2).

• n_jobs (int) – Give the number of CPUs.

Returns:

• r (ndarray) – sizes of the regions

• dn (ndarray) – probability of finding the corresponding size

tools21cm.bubble_stats.mfp(data, xth=0.5, boxsize=None, iterations=10000000, verbose=True, upper_lim=False, bins=None, r_min=None, r_max=None)

Determines the sizes using the Mean-Free-Path (MFP) approach.

Parameters:

• input (ndarray) – 2D/3D array of ionization fraction/brightness temperature.

• xth (float) – The threshold value (Default: 0.5).

• boxsize (float) – The boxsize in cMpc can be given (Default: conv.LB).

• iterations (int) – Number of iterations (Default: 10_000_000).

• verbose (bool) – It prints the progress of the program (Default: True).

• upper_lim (bool) – It decides if the threshold is the upper limit or the lower limit (Default: False).

• bins (int) – Give number of bins or an array of sizes to re-bin into (Default: None).

• r_min (float) – Minimum size after rebinning (Default: None).

• r_max (float) – Maximum size after rebinning (Default: None).

Returns:

• r (ndarray) – sizes of the regions

• dn (ndarray) – probability of finding the corresponding size

tools21cm.bubble_stats.mfp_from_point(data, point, xth=0.5, boxsize=None, iterations=10000000, verbose=True, upper_lim=False, bins=None, r_min=None, r_max=None)

Determines the sizes using the Mean-Free-Path (MFP) approach.

Parameters:

• input (ndarray) – 2D/3D array of ionization fraction/brightness temperature.

• point (ndarray or list) – [x, y, z] of the point.

• xth (float) – The threshold value (Default: 0.5).

• boxsize (float) – The boxsize in cMpc can be given (Default: conv.LB).

• iterations (float) – Number of iterations (Default: 1e7).

• verbose (bool) – It prints the progress of the program (Default: True).

• upper_lim (bool) – It decides if the threshold is the upper limit or the lower limit (Default: False).

• bins (int) – Give number of bins or an array of sizes to re-bin into (Default: None).

• r_min (float) – Minimum size after rebinning (Default: None).

• r_max (float) – Maximum size after rebinning (Default: None).

Returns:

• r (ndarray) – sizes of the regions

• dn (ndarray) – probability of finding the corresponding size

tools21cm.bubble_stats.spa(data, xth=0.95, boxsize=None, nscales=20, upper_lim=False, binning='log', verbose=True, use_fftconvolve=True)

Determines the sizes using the Spherical-Averege (SPA) approach.

Parameters:

• input (ndarray) – 3D array of ionization fraction.

• xth (float) – The threshold value (Default: 0.5).

• boxsize (float) – The boxsize in cMpc can be given (Default: conv.LB).

• nscales (int) – The number of different radii to consider (Default: 20).

• upper_lim (bool) – It decides if the threshold is the upper limit or the lower limit (Default: True).

• verbose (bool) – Print verbose (Default: True).

• use_fftconvolve (bool) – Use FFT-based convolution (Default: True).

Returns:

• r (ndarray) – sizes of the regions

• dn (ndarray) – probability of finding the corresponding size

tools21cm.bubble_stats.watershed_bsd(data, xth=0.5, boxsize=None, verbose=True, upper_lim=False, h=None, n_jobs=1)

Still needs testing.

Determined the sizes using the Watershed approach. It is based on Lin et al. (2016)

Parameters:

• input (ndarray) – 2D/3D array of ionization fraction/brightness temperature.

• xth (float) – The threshold value (Default: 0.5).

• boxsize (float) – The boxsize in cMpc can be given (Default: conv.LB).

• verbose (bool) – It prints the progress of the program (Default: True).

• upper_lim (bool) – It decides if the threshold is the upper limit or the lower limit (Default: False).

• h (float) – Parameter for performing h-transform to smooth out local peaks.

• n_jobs (int) – Give the number of CPUs.

Returns:

• r (ndarray) – sizes of the regions

• dn (ndarray) – probability of finding the corresponding size

Cosmological calculators

Various functions for calculating some cosmological stuff

tools21cm.cosmo.M_to_Tvir(M, z)

Convert mass to virial temperature.

Parameters:

• M (float or array) – The mass(es) in solar mass unit.

• z (float) – the redshift.

Returns:

Virial temperature in K.

tools21cm.cosmo.Tvir_to_M(Tvir, z)

Convert virial temperature to mass.

Parameters:

• Tvir (float or array) – The virial temperature(s) in K.

• z (float) – the redshift.

Returns:

Mass in solar mass unit.

tools21cm.cosmo.angular_size(dl, z)

Calculate the angular size of an object at a given redshift.

Parameters:

• dl (float or array) – the physical size in kpc

• z (float or array) – the redshift of the object

Returns:

The angluar size in arcseconds

tools21cm.cosmo.angular_size_comoving(cMpc, z)

Calculate the angular size in degrees of an object with a given comoving size.

Parameters:

• cMpc (float or array) – the size in comoving Mpc

• z (float or array) – the redshift of the object

Returns:

The angular size in degrees

tools21cm.cosmo.c_to_p(z_to_cdist, z)

Convert comoving distance to proper distance

Parameters:

• z_to_cdist (float or array) – The comoving distance

• z (float) – the redshift

Returns:

Proper distance

tools21cm.cosmo.cdist_to_z(dist)

Calculate the redshift correspoding to the given comoving distance.

Parameters:

dist (float or array) – the distance in comoving Mpc

Returns:

redshift corresponding to the distance.

Note

Uses a precalculated table for interpolation. Only valid for 0 <= z < 100

tools21cm.cosmo.deg_to_cdist(deg, z)

Calculate the size in cMpc of an object with given angular diameter.

Parameters:

• deg (float or array) – the size in degrees

• z (float or array) – the redshift

Returns:

The size in cMpc

tools21cm.cosmo.hubble_parameter(z)

It calculates the Hubble parameter at any redshift in flat LCDM cosmology.

tools21cm.cosmo.luminosity_distance(z, k=0)

Calculate the luminosity distance for a given redshift.

Parameters:

• z (float or array) – the redshift(s)

• k (float) – the curvature constant.

Returns:

The luminosity distance in Mpc

tools21cm.cosmo.nu_to_cdist(nu21)

Calculate the comoving distance to a given 21 cm frequency

Parameters:

nu21 (float or array) – redshifted 21 cm frequency in MHz

Returns:

Comoving distance in Mpc

tools21cm.cosmo.nu_to_wavel(nu)

Convert frequency to wavelength

Parameters:

nu (float or array) – the frequency in MHz

Returns:

The wavelength in meters

tools21cm.cosmo.nu_to_z(nu21)

Convert 21 cm frequency in MHz to redshift

Parameters:

nu21 (float or array) – redshifted 21 cm frequency in MHz

Returns:

Redshift

tools21cm.cosmo.p_to_c(pdist, z)

Convert proper distance to comoving distance

Parameters:

• pdist (float or array) – The proper distance

• z (float) – the redshift

Returns:

Comoving distance

tools21cm.cosmo.z_to_age(z)

Calculate the age of the universe

Parameters:

z (float or array) – redshift

Returns:

Comoving distance in Gyr

tools21cm.cosmo.z_to_cdist(z)

Calculate the comoving distance

Parameters:

z (float or array) – redshift

Returns:

Comoving distance in Mpc

tools21cm.cosmo.z_to_nu(z)

Get the 21 cm frequency that corresponds to redshift z

Parameters:

z (float or array) – redshift

Returns:

redshifted 21 cm frequency in MHz

Foreground model

Methods to simulate and analyse the foreground signal for 21 cm signal.

tools21cm.foreground_model.extragalactic_pointsource_fg(z, ncells, boxsize, rseed=False, S_max=100)

@ Ghara et al. (2017)

Parameters:

• z (float) – Redshift observed with 21-cm.

• ncells (int) – Number of cells on each axis.

• boxsize (float) – Size of the FOV in Mpc.

• S_max (float) – Maximum flux of the point source to model in muJy (Default: 100).

Returns:

A 2D numpy array of brightness temperature in mK.

tools21cm.foreground_model.galactic_synch_fg(z, ncells, boxsize, rseed=False)

@ Ghara et al. (2017)

Parameters:

• z (float) – Redshift observed with 21-cm.

• ncells (int) – Number of cells on each axis.

• boxsize (float) – Size of the FOV in Mpc.

• rseed (int) – random seed to have the same realisation (Default: False).

Returns:

A 2D numpy array of brightness temperature in mK.

tools21cm.foreground_model.rolling_wedge_removal_lightcone(lightcone, redshifts, cell_size=None, chunk_length=None, OMm=None, buffer_threshold=1e-10)

Rolling over the lightcone and removing the wedge for every frequency channel.

Parameters:

• lightcone (array) – The lightcone, of shape (n_cells, n_cells, n_redshfts).

• redshifts (array) – The redshifts of every frequency channel.

• cell_size (float) – Resolution of the lightcone voxels. It is assumed that the lightcone is made of cubical voxels of volume cell_size^3.

• chunk_length (int) – Length of the box in z-direction used for wedge removal. Defaults to n_cells.

• OMm (float) – Omega matter.

• buffer_threshold (float) – Threshold which defines a wedge buffer. Buffer is then calculated as k-value for which Blackman taper power is below threshold.

Returns:

Cleaned lightcone.

Identifying regions

Methods to identify regions of interest in images.

tools21cm.identify_regions.bubbles_from_fixed_threshold(data, threshold=0, upper_lim=True)

@ Giri at al. (2018a)

It is a method to identify regions of interest in noisy images. The method uses a fixed threshold.

Parameters:

• data (ndarray) – The brightness temperature or ionization fraction cube.

• threshold (float) – The fixed threshold value (Default: 0).

• upper_lim (bool) – This decides which mode in the PDF is to be identified. ‘True’ identifies ionized regions in brightness temperature, while ‘False’ identifies in the xfrac data (Default: True).

Returns:

Binary cube where pixels identified as region of interest are the True.

tools21cm.identify_regions.bubbles_from_kmeans(data, upper_lim=True, n_jobs=1, n_clusters=3)

@ Giri at al. (2018a)

It is a method to identify regions of interest in noisy images. The method finds the optimal threshold using the 1D PDF of the image. It gives similar results compared to the Otsu’s method.

Parameters:

• data (ndarray) – The brightness temperature/xfrac cube.

• upper_lim (bool) – This decides which mode in the PDF is to be identified. ‘True’ identifies ionized regions in brightness temperature, while ‘False’ identifies in the xfrac data (Default: True).

• n_jobs (int) – Number of cores to use (Default: 1).

• n_cluster (int) – Number of clusters found in the PDF (Default: 3).

Returns:

Binary cube where pixels identified as region of interest are the True.

tools21cm.identify_regions.bubbles_from_slic(data, n_segments=5000, bins='knuth', verbose=True)

@ Giri at al. (2018b) It is a method to identify regions of interest in noisy images. The method is an implementation of the superpixel based approach, called SLIC (Achanta et al. 2010), used in the field of computer vision.

Parameters:

• data (ndarray) – The brightness temperature cube.

• n_segments (int) – Number of superpixels (Default: 2000).

• bins (int or str) – Number of bins for the PDF used for stitching. ‘blocks’, ‘knuth’, ‘scotts’, ‘freedman’ rules to determine bins automatically can also be choosen. (Default: ‘knuth’).

Returns:

Binary cube where pixels identified as region of interest are the True.

tools21cm.identify_regions.bubbles_from_triangle2D(data, upper_lim=True, nbins=256, loc_axis=2)

Gazagnes et al. (2021) https://arxiv.org/abs/2011.08260 Giri at al. (2018b) https://arxiv.org/abs/1801.06550 Zack et al. (1977) https://journals.sagepub.com/doi/10.1177/25.7.70454

It is a method to identify regions of interest in noisy images. The method uses a triangle threshold proposed in Zack et al. (1977) for 2D images. Here we use the technique proposed in Gazagnes et al. (2021) and determine threshold from 2D slices and determine take the median value.

Parameters:

• data (ndarray) – The brightness temperature or ionization fraction cube.

• upper_lim (bool) – This decides which mode in the PDF is to be identified. ‘True’ identifies ionized regions in brightness temperature, while ‘False’ identifies in the xfrac data (Default: True).

• nbins (int) – The number of bins used for histogram (Default: 256).

Returns:

Binary cube where pixels identified as region of interest are the True.

tools21cm.identify_regions.bubbles_from_triangle3D(data, upper_lim=True, nbins=256)

Gazagnes et al. (2021) https://arxiv.org/abs/2011.08260 Giri at al. (2018b) https://arxiv.org/abs/1801.06550 Zack et al. (1977) https://journals.sagepub.com/doi/10.1177/25.7.70454

It is a method to identify regions of interest in noisy images. The method uses a triangle threshold proposed in Zack et al. (1977) for 2D images. Here we use the same method for 3D data as proposed in Giri at al. (2018b) and Gazagnes et al. (2021).

Parameters:

• data (ndarray) – The brightness temperature or ionization fraction cube.

• upper_lim (bool) – This decides which mode in the PDF is to be identified. ‘True’ identifies ionized regions in brightness temperature, while ‘False’ identifies in the xfrac data (Default: True).

• nbins (int) – The number of bins used for histogram (Default: 256).

Returns:

Binary cube where pixels identified as region of interest are the True.

Lightcone

Methods to construct lightcones.

tools21cm.lightcone.get_lightcone_subvolume(lightcone, redshifts, central_z, depth_mhz=None, depth_mpc=None, odd_num_cells=True, subtract_mean=True, fov_Mpc=None, verbose=True)

Extract a subvolume from a lightcone, at a given central redshift, and with a given depth. The depth can be specified in Mpc or MHz. You must give exactly one of these parameters.

Parameters:

• ligthcone (numpy array) – the lightcone

• redshifts (numpy array) – the redshifts along the LOS

• central_z (float) – the central redshift of the subvolume

• depth_mhz (float) – the depth of the subvolume in MHz

• depth_mpc (float) – the depth of the subvolume in Mpc

• odd_num_cells (bool) – if true, the depth of the box will always be an odd number of cells. This avoids problems with power spectrum calculations.

• subtract_mean (bool) – if true, subtract the mean of the signal (Default: True)

• fov_Mpc (float) – the FoV size in Mpc

• verbose (bool) – turning the verbosity.

Returns:

Tuple with (subvolume, dims, sub_redshifts) where dims is a tuple with the dimensions of the subvolume in Mpc

tools21cm.lightcone.make_lightcone(filenames, z_low=None, z_high=None, depth_mhz=None, file_redshifts=None, cbin_bits=32, cbin_order='c', los_axis=0, raw_density=False, interpolation='linear', reading_function=None, box_length_mpc=None)

Generate a lightcone from xfrac, density, or dT data cubes.

tools21cm.lightcone.make_velocity_lightcone(vel_filenames, dens_filenames, z_low=None, z_high=None, file_redshifts=None, los_axis=0, vel_reading_function=None, dens_reading_function=None)

Generate a lightcone from velocity data. Since velocity files contain momentum rather than actual velocity, filenames for both velocity and density must be specified.

tools21cm.lightcone.redshifts_at_equal_comoving_distance(z_low, z_high, box_grid_n=256, box_length_mpc=None)

Make a frequency axis vector with equal spacing in co-moving LOS coordinates. The comoving distance between each frequency will be the same as the cell size of the box.

Parameters:

• z_low (float) – The lower redshift

• z_high (float) – The upper redhisft

• box_grid_n = 256 (int) – the number of slices in an input box

• box_length_mpc (float) – the size of the box in cMpc. If None,

• set to conv.LB

Returns:

numpy array containing the redshifts

tools21cm.lightcone.redshifts_at_equal_frequncy_width(z_low, z_high, dnu)

Make a redshift vector with equally space in the frequency domain.

Parameters:

• z_low (float) – The lower redshift

• z_high (float) – The upper redhisft

• dnu (float) – The frequency channel width in MHz

Returns:

numpy array containing the redshifts

N-Body Tools

class tools21cm.nbody_file.NbodyFile(filespath=None, z=None, node=None, alt_format=False, pid_file=True)

A CubeP3M Nbody file.

Use the read_from_file method to load a Nbody file, or pass the filespath to the constructor.

Some useful attributes of this class are:

nrpart (int): total number of particles a (float): the scale factor of the file z (float): the redshift of the file

Methods

get_npart
read_from_file

class tools21cm.nbody_file.NbodyParticle

A simple struct to hold info about a single halo

tools21cm.nbody_file.halo_list_to_grid(mass, pos_xyz, box_dim, n_grid)

Put the list of haloes on a grid.

Parameters:

• mass (ndarray) – Mass of haloes.

• pos_xyz (ndarray) – Position of the haloes.

• box_dim (float) – Length of simulation in each direction.

• n_grid (int) – Number of grids along each direction.

Returns:

• binned_mhalo (ndarray) – The haloes put on grids of shape (nGrid, nGrid, nGrid).

• bin_edges (ndarray) – Edges of the bins used to create the grids.

• bin_num (ndarray) – Number of bins.

tools21cm.nbody_file.read_pkdgrav_density_grid(filename, box_dim, nGrid=None)

Loads the density field from a file and computes the density contrast.

Parameters:

file (str) – Path to the pkdgrav density field file.

Returns:

delta (ndarray) – The density contrast delta, defined as (rho_m / rho_mean - 1). It is a 3-D mesh grid of size (nGrid, nGrid, nGrid).

tools21cm.nbody_file.read_pkdgrav_fof_halo_catalogue(filename, box_dim, nGrid=None, hlittle=0.67)

Loads the fof data created by the FOF halo finder implemented in the pkdgrav framework.

Parameters:

• filename (str) – Path to the fof data file.

• box_dim (float) – Length of the simulation volume used while running the pkdgrav code in Mpc.

Returns:

fof_data (ndarray) – An array of halo positions and mass [mass, pos_x, pos_y, pos_z].

Point statistics

This file contains various useful statistical methods

tools21cm.statistics.apply_func_along_los(signal, func, los_axis)

Apply a function, such as np.var() or np.mean(), along the line-of-sight axis of a signal on a per-slice basis.

Parameters:

• signal (ndarray) – the signal

• func (callable) – the function to apply

• los_axis (int) – the line-of-sight axis

Returns:

An array of length signal.shape[los_axis]

Example

Calculate the variance of a lightcone along the line-of-sight:

>>> lightcone = t2c.read_cbin('my_lightcone.cbin')
>>> dT_var = t2c.apply_func_along_los(lightcone, np.var, 2)

tools21cm.statistics.kurtosis(x)

Calculate the kurtosis of an array. It uses the definition given in Ross et al. (2017).

Parameters:

x (ndarray) – The array containing the input data

Returns:

The kurtosis.

tools21cm.statistics.mass_weighted_mean_xi(xi, rho)

Calculate the mass-weighted mean ionization fraction.

Parameters:

• xi (ndarray) – the ionized fraction

• rho (ndarray) – the density (arbitrary units)

Returns:

The mean mass-weighted ionized fraction.

tools21cm.statistics.signal_overdensity(signal, los_axis)

Divide by the mean of the signal along the los axis and subtract one.

Parameters:

• signal (ndarray) – the signal to subtract the mean from

• los_axis (int) – the line-of-sight axis

Returns:

The signal with the mean subtracted

TODO:vectorize

tools21cm.statistics.skewness(x)

Calculate the skewness of an array. Note that IDL calculates the skewness in a slightly different way than Python. This routine uses the IDL definition.

Parameters:

x (ndarray) – The array containing the input data.

Returns:

The skewness.

tools21cm.statistics.subtract_mean_signal(signal, los_axis=2)

Subtract the mean of the signal along the los axis.

Parameters:

• signal (ndarray) – the signal to subtract the mean from

• los_axis (int) – the line-of-sight axis (Default: 2)

Returns:

The signal with the mean subtracted

TODO:vectorize

Power spectrum

Contains functions to estimate various two point statistics.

tools21cm.power_spectrum.anisotropy_ratio_r_mu(input_array, los_axis=0, mu_cut=0.5, mubins=20, kbins=10, box_dims=None, weights=None, exclude_zero_modes=True, return_n_modes=False, absolute_mus=True, k_limit=None, **kwargs)

Calculate the anisotropy ratio r_mu(k).

Parameters:

• input_array (numpy array) – the data array

• los_axis = 0 (integer) – the line-of-sight axis

• mu_cut = 0.5 (float) – the value of mu to split the P(k,mu) data.

• mubins = 20 (integer) – the number of mu bins

• kbins = 10 (integer or array-like) – The number of bins, or a list containing the bin edges. If an integer is given, the bins are logarithmically spaced.

• box_dims = None (float or array-like) – the dimensions of the box in Mpc. If this is None, the current box volume is used along all dimensions. If it is a float, this is taken as the box length along all dimensions. If it is an array-like, the elements are taken as the box length along each axis.

• return_n_modes = False (bool) – if true, also return the number of modes in each bin

• exlude_zero_modes = True (bool) – if true, modes with any components of k equal to zero will be excluded.

• absolute_mus = True (boolean) – if true, use the absolute values of mu, range [0,1]. If false, use the range [-1,1]

• k_limit = None – Upper k limit for the bins. Options: None — use the grid maximum (default, no restriction) ‘nyquist’ — cap at k_Nyquist = π·min(N_i/L_i) ‘aliasing’ — cap at k_Nyquist/2 (alias-free limit)

Returns:

A tuple with (rmu, kbins).

tools21cm.power_spectrum.cross_power_spectrum_1d(input_array1_nd, input_array2_nd, kbins=100, box_dims=None, return_n_modes=False, binning='log', k_limit=None, fill_empty_bins=True, **kwargs)

Calculate the spherically averaged cross power spectrum of two arrays and return it as a one-dimensional array.

Parameters:

• input_array1_nd (numpy array) – the first data array

• input_array2_nd (numpy array) – the second data array

• kbins = 100 (integer or array-like) – number of bins or explicit bin edges.

• box_dims = None (float or array-like) – box dimensions in Mpc.

• return_n_modes = False (bool) – if True, also return n_modes array.

• binning = ‘log’ – binning scheme (‘log’, ‘linear’, ‘adaptive’, ‘kde’).

• k_limit = None – upper k limit (None, ‘nyquist’, or ‘aliasing’).

• fill_empty_bins = True – fill zero-mode bins via 1D nearest-neighbour.

Returns:

(Pk, bins) or (Pk, bins, n_modes) if return_n_modes is True.

tools21cm.power_spectrum.cross_power_spectrum_mu(input_array1, input_array2, los_axis=0, mubins=20, kbins=10, box_dims=None, weights=None, exclude_zero_modes=True, return_n_modes=False, absolute_mus=True, k_limit=None, **kwargs)

Calculate the cross power spectrum and bin it in mu=cos(theta) and k.

Parameters:

• input_array1 (numpy array) – the first data array

• input_array2 (numpy array) – the second data array

• los_axis = 0 (integer) – the line-of-sight axis

• mubins = 20 (integer) – the number of mu bins

• kbins = 10 (integer or array-like) – The number of bins, or a list containing the bin edges. If an integer is given, the bins are logarithmically spaced.

• box_dims = None (float or array-like) – the dimensions of the box in Mpc. If this is None, the current box volume is used along all dimensions. If it is a float, this is taken as the box length along all dimensions. If it is an array-like, the elements are taken as the box length along each axis.

• return_n_modes = False (bool) – if true, also return the number of modes in each bin

• exlude_zero_modes = True (bool) – if true, modes with any components of k equal to zero will be excluded.

• absolute_mus = True (boolean) – if true, use the absolute values of mu, range [0,1]. If false, use the range [-1,1]

Returns:

A tuple with (Pk, mubins, kbins), where Pk is an array with the cross power spectrum of dimensions (n_mubins x n_kbins), mubins is an array with the mu bin centers and kbins is an array with the k bin centers.

tools21cm.power_spectrum.cross_power_spectrum_nd(input_array1, input_array2, box_dims, **kwargs)

Calculate the cross power spectrum two arrays and return it as an n-dimensional array.

Parameters:

• input_array1 (numpy array) – the first array to calculate the power spectrum of. Can be of any dimensions.

• input_array2 (numpy array) – the second array. Must have same dimensions as input_array1.

• box_dims = None (float or array-like) – the dimensions of the box in Mpc. If this is None, the current box volume is used along all dimensions. If it is a float, this is taken as the box length along all dimensions. If it is an array-like, the elements are taken as the box length along each axis.

Returns:

The cross power spectrum in the same dimensions as the input arrays.

tools21cm.power_spectrum.dimensionless_ps(data, kbins=100, box_dims=None, binning='log', factor=10, k_limit=None)

Dimensionless power spectrum is P(k)*k^3/(2pi^2)

Parameters:

• data (ndarray) – The numpy data whose power spectrum is to be determined.

• kbins (int) – Number of bins for in the k-space (Default: 100).

• box_dims (float) – The size of the box in Mpc (Default: Takes the value from the set_sim_constants).

• binning (str) – The type of binning to be used for the k-space (Default: ‘log’).

• factor (int) – The factor multiplied to the given kbins to smooth the spectrum from (Default: 10).

• k_limit (None, ‘nyquist’, or ‘aliasing’) – Upper k limit passed to power_spectrum_1d (Default: None).

Returns:

(Delta^2, ks)

tools21cm.power_spectrum.horizon_wedge_equation(z, fov_deg=90.0)

Horizon wedge boundary in the cylindrical power spectrum. Equation 14 in arXiv:2201.10798.

Parameters:

• z (float) – redshift

• fov_deg (float) – field of view in degrees (default 90)

Returns:

A lambda function k_parallel(k_perpendicular).

tools21cm.power_spectrum.mu_binning(powerspectrum, los_axis=0, mubins=20, kbins=10, box_dims=None, weights=None, exclude_zero_modes=True, binning='log', absolute_mus=True, kmax=None)

This function is for internal use only.

tools21cm.power_spectrum.plot_2d_power(ps, xlabel='$k_\\perp$', ylabel='$k_\\parallel$', ps_label='$P(k_\\perp,k_\\parallel)$', fig=None, ax=None, plotting_scale=None, draw_wedge=None, **kwargs)

Plot the 2D cylindrical power spectrum.

Parameters:

• ps (tuple or dict) – (Pk, kper_bins, kpar_bins) or {‘Pk’: Pk, ‘kper’: kper_bins, ‘kpar’: kpar_bins}

• xlabel, ylabel, ps_label – axis and colorbar labels

• fig, ax – matplotlib figure and axes (created if None)

• plotting_scale (dict) – axis scales, e.g. {‘x’: ‘log’, ‘y’: ‘log’, ‘z’: ‘log’}

• draw_wedge (dict or None) – wedge params e.g. {‘z’: 9.0, ‘fov_deg’: 90.0, ‘ls’: ‘–’, ‘color’: ‘k’}

Returns:

matplotlib figure

tools21cm.power_spectrum.power_spectrum_1d(input_array_nd, kbins=100, box_dims=None, return_n_modes=False, binning='log', window=None, k_limit=None, fill_empty_bins=True, **kwargs)

Calculate the spherically averaged power spectrum of an array and return it as a one-dimensional array.

Parameters:

• input_array_nd (numpy array) – the data array

• kbins = 100 (integer or array-like) – number of bins or explicit bin edges.

• box_dims = None (float or array-like) – box dimensions in Mpc.

• return_n_modes = False (bool) – if True, also return n_modes array.

• binning = ‘log’ – binning scheme. Options: ‘log’ — uniform log-spaced bins (default) ‘linear’ — uniform linear-spaced bins ‘adaptive’ — edges at geometric midpoints between discrete

  k-modes, centers at geometric means; no empty bins

  ‘adaptive_linear’ — same but arithmetic midpoints/means (linear-k) ‘kde’ — Gaussian kernel in log-k space (bandwidth = 1 bin

  spacing); avoids empty bins, mild cross-bin mixing

  ‘kde_linear’ — Gaussian kernel in linear-k space

• window = None – taper (‘blackmanharris’, ‘tukey’, or array).

• k_limit = None – upper k limit (None, ‘nyquist’, or ‘aliasing’).

• fill_empty_bins = True – fill zero-mode bins via 1D nearest-neighbour after tophat/adaptive binning. No effect for kde variants.

Returns:

(Pk, bins) or (Pk, bins, n_modes) if return_n_modes is True.

tools21cm.power_spectrum.power_spectrum_2d(input_array, kbins=10, binning='log', box_dims=348.5714285714286, return_modes=False, nu_axis=2, window=None, k_limit=None, fill_empty_bins=True, **kwargs)

Calculate the power spectrum and bin it in kper and kpar.

Parameters:

• input_array (numpy array) – the data array

• nu_axis = 2 (integer) – the line-of-sight axis

• kbins = 10 (integer or [n_kper, n_kpar]) – number of bins per direction.

• box_dims = 244/.7 (float or array-like) – box dimensions in Mpc.

• return_modes = False (bool) – if True, also return n_modes array.

• binning = ‘log’ – binning scheme. Options: ‘log’ — uniform log-spaced bins (default) ‘linear’ — uniform linear-spaced bins ‘adaptive’ — edges at geometric midpoints between discrete

  k-modes, centers at geometric means; no empty bins

  ‘adaptive_linear’ — same but arithmetic midpoints/means (linear-k) ‘kde’ — Gaussian kernel in log-k space (bandwidth = 1 bin

  spacing); avoids empty bins, mild cross-bin mixing

  ‘kde_linear’ — Gaussian kernel in linear-k space

• window = None – taper (‘blackmanharris’, ‘tukey’, or array).

• k_limit = None – upper k limit (None, ‘nyquist’, or ‘aliasing’).

• fill_empty_bins = True – after tophat/adaptive binning, fill any remaining NaN (zero-mode) bins via 2D nearest-neighbour in index space. Set False to keep NaN. Ignored for kde variants.

Returns:

(Pk, kper_mid, kpar_mid) or (Pk, kper_mid, kpar_mid, n_modes).

tools21cm.power_spectrum.power_spectrum_mu(input_array, los_axis=0, mubins=20, kbins=10, box_dims=None, weights=None, exclude_zero_modes=True, return_n_modes=False, absolute_mus=True, k_limit=None, **kwargs)

Calculate the power spectrum and bin it in mu=cos(theta) and k.

Parameters:

• input_array (numpy array) – the data array

• los_axis = 0 (integer) – the line-of-sight axis

• mubins = 20 (integer) – the number of mu bins

• kbins = 10 (integer or array-like) – The number of bins, or a list containing the bin edges. If an integer is given, the bins are logarithmically spaced.

• box_dims = None (float or array-like) – the dimensions of the box in Mpc. If this is None, the current box volume is used along all dimensions. If it is a float, this is taken as the box length along all dimensions. If it is an array-like, the elements are taken as the box length along each axis.

• return_n_modes = False (bool) – if true, also return the number of modes in each bin

• exlude_zero_modes = True (bool) – if true, modes with any components of k equal to zero will be excluded.

• absolute_mus = True (boolean) – if true, use the absolute values of mu, range [0,1]. If false, use the range [-1,1]

• k_limit = None – Upper k limit for the bins. Options: None — use the grid maximum (default, no restriction) ‘nyquist’ — cap at k_Nyquist = π·min(N_i/L_i) ‘aliasing’ — cap at k_Nyquist/2 (alias-free limit)

Returns:

A tuple with (Pk, mubins, kbins), where Pk is an array with the power spectrum of dimensions (n_mubins x n_kbins), mubins is an array with the mu bin centers and kbins is an array with the k bin centers.

tools21cm.power_spectrum.power_spectrum_nd(input_array, box_dims=None, verbose=False, **kwargs)

Calculate the power spectrum of input_array and return it as an n-dimensional array.

Parameters:

• input_array (numpy array) – the array to calculate the power spectrum of. Can be of any dimensions.

• box_dims = None (float or array-like) – the dimensions of the box in Mpc. If this is None, the current box volume is used along all dimensions. If it is a float, this is taken as the box length along all dimensions. If it is an array-like, the elements are taken as the box length along each axis.

Returns:

The power spectrum in the same dimensions as the input array.

tools21cm.power_spectrum.radial_average(input_array, box_dims, kbins=10, binning='log', kmax=None, fill_empty_bins=True, **kwargs)

Radially average data. Mostly for internal use.

Parameters:

• input_array (numpy array) – the data array

• box_dims (float or array-like) – box dimensions in Mpc.

• kbins = 10 (integer or array-like) – number of bins or explicit bin edges.

• binning = ‘log’ – binning scheme (‘log’, ‘linear’, ‘adaptive’, ‘adaptive_linear’, ‘kde’, ‘kde_linear’).

• kmax = None (float) – upper k bin edge (from _resolve_k_limit).

• fill_empty_bins = True – fill zero-mode bins via 1D nearest-neighbour.

Returns:

(data, bin_centers, n_modes)

Primary beam

Methods to simulate the primary beam of radio telescope

tools21cm.primary_beam.circular_beam(ncells, z, D_station=40.0, beam_func='Gaussian', boxsize=None, l_null=None)

Generates a 2D circular beam pattern.

tools21cm.primary_beam.primary_beam(array, z, nu_axis=2, beam_func='Gaussian', boxsize=None, D_station=40.0, l_null=None, verbose=True)

arrayndarray

The array of brightness temperature.

zfloat

redshift. With one value, the function will assume the array to be coeval.

In case of light-cone, provide size of z should be equivalent to the length of frequency axis of the array.

nu_axisint

The frequency axis of the array (Default: 2)

beam_func: str

The type of function to model the primary beam. The options are ‘gaussian’, ‘bessel’, ‘sigmoid’ and ‘step’. Default: ‘gaussian’

boxsizefloat

Size of the box in physical units (cMpc). Default: From set simulation constants.

D_station: float

Diameter of the station in metres. Default: 40.

l_null: float

Length (in Mpc) of the first null. Default: None,

which will use the diameter of the station to estimate this value.

tools21cm.primary_beam.primary_beam_null(z, D_station=40.0)

Calculates the physical diameter of the first null of an Airy disk projected on the sky at redshift z.

Reading simuation outputs

Methods to read simulations outputs of few popular codes.

tools21cm.read_files.coeval_21cm_c2ray(xfrac_dir, dens_dir, z, interpolation='linear', mean_subtract=False)

Parameters:

• xfrac_dir (str) – Give the path that contains the xfrac-files.

• dens_dir (str) – Give the path that contains the density-files.

• z (float) – Redshift.

• interpolation (str) – This is used when the xfrac cube at that redshift is not available.

Returns:

numpy array of brightness temperature in mK.

tools21cm.read_files.coeval_dens_c2ray(dens_dir, z)

Parameters:

• dens_dir (str) – Give the path that contains the density-files.

• z (float) – Redshift.

Returns:

numpy array

tools21cm.read_files.coeval_xfrac_c2ray(xfrac_dir, z, interpolation='linear')

Parameters:

• xfrac_dir (str) – Give the path that contains the xfrac-files.

• z (float) – Redshift.

• interpolation (str) – This is used when the coveal cube at that redshift is not available.

Returns:

numpy array

tools21cm.read_files.get_package_resource_path(package_name, resource_path)

Get the path to a resource file within a package, with automatic fallback for older Python versions.

Parameters:

• package_name (str) – Name of the package (e.g., ‘tools21cm’)

• resource_path (str) – Path to the resource relative to the package root (e.g., ‘input_data/central_geographic_position.txt’)

Returns:

str or Path – Path to the resource file that can be used with file reading functions

Examples

>>> path = get_package_resource_path('tools21cm', 'input_data/central_geographic_position.txt')
>>> data = np.loadtxt(path)

tools21cm.read_files.read_21cmfast_files(filename)

Parameters:

filename (str) – Give the filename of the 21cmFAST output files.

Returns:

numpy array

tools21cm.read_files.read_c2ray_files(filename, file_type='xfrac', density_file=None)

Parameters:

• filename (str) – Give the filename of the C2Ray output files.

• file_type (str) –

  The file of file being read has to be mentioned. The options are

  ‘xfrac’, ‘dens’, ‘vel’. (Default: ‘xfrac’)

• density_file (str) – This is necessary if file_type==’vel’. (Default: None)

Returns:

numpy array

tools21cm.read_files.read_dictionary_data(filename, format=None)

Reads a dictionary from various file formats. Supported formats: pickle, json, yaml, hdf5, npz, netcdf (xarray), csv, xlsx

tools21cm.read_files.read_grizzly_files(filename)

Parameters:

filename (str) – Give the filename of the GRIZZLY xfrac files.

Returns:

numpy array

tools21cm.read_files.write_dictionary_data(data_dict, filename, format=None)

Writes a dictionary to various file formats. Supports: pickle, json, yaml, hdf5, npz, netcdf (xarray), csv, xlsx

Redshift space distortion

tools21cm.pv_mpm.get_distorted_dt(dT, kms, redsh, los_axis=0, velocity_axis=0, num_particles=10, periodic=True)

Apply peculiar velocity distortions to a differential temperature box, using the Mesh-Particle-Mesh method, as described in http://arxiv.org/abs/1303.5627

Parameters:

• * dT (numpy array) – the differential temperature box

• * kms (numpy array) – velocity in km/s, array of dimensions (3,mx,my,mz) where (mx,my,mz) is dimensions of dT

• * redsh (float) – the redshift

• * los_axis = 0 (int) – the line-of-sight axis of the output volume (must be 0, 1 or 2)

• * velocity_axis = 0 (int) – the index that indicates los velocity

• * num_particles = 10 (int) – the number of particles to use per cell A higher number gives better accuracy, but worse performance.

• * periodic = True (bool) – whether or not to apply periodic boundary conditions along the line-of-sight. If you are making a lightcone volume, this should be False.

Returns:

The redshift space box as a numpy array with same dimensions as dT.

Example

Read a density file, a velocity file and an xfrac file, calculate the brightness temperature, and convert it to redshift space.

>>> vfile = c2t.VelocityFile('/path/to/data/8.515v_all.dat')
>>> dfile = c2t.DensityFile('/path/to/data/8.515n_all.dat')
>>> xfile = c2t.XfracFile('/path/to/data/xfrac3d_8.515.bin')
>>> dT = c2t.calc_dt(xfile, dfile)
>>> kms = vfile.get_kms_from_density(dfile)
>>> dT_zspace = get_distorted_dt(dT, kms, dfile.z, los_axis = 0)

Note

At the moment, it is a requirement that dimensions perpendicular to the line-of-sight are equal. For example, if the box dimensions are (mx, my, mz) and the line-of-sight is along the z axis, then mx has to be equal to my.

Note

If dT is a lightcone volume, los_axis is not necessarily the same as velocity_axis. The lightcone volume methods in c2raytools all give output volumes that have the line-of-sight as the last index, regardless of the line-of-sight axis. For these volumes, you should always use los_axis=2 and set velocity_axis equal to whatever was used when producing the real-space lightcones.

tools21cm.pv_mpm.make_pv_box(dT_filename, vel_filename, dens_filename, z, los=0, num_particles=10)

Convenience method to read files and make a distorted box.

Parameters:

• * dT_filename (string) – the name of the dT file

• * vel_filename (string) – the name of the velocity file

• * dens_filename (string) – the name of the density file

• * z (float) – the redshift

• * los (integer) – the line-of-sight axis

• * num_particles (integer) – the number of particles to pass to get_distorted_dt

Returns:

The redshift space box

Smoothing

Methods to smooth or reduce resolution of the data to reduce noise.

tools21cm.smoothing.convolve_uvmap_coeval(cube, z, box_size_mpc=False, max_baseline=2.0, ratio=1.0, nu_axis=2, verbose=True, filename=None, total_int_time=6.0, int_time=10.0, declination=-30.0, uv_map=None)

This smooths the coeval cube by Gaussian in angular direction and by tophat along the third axis.

Parameters:

• coeval_cube (numpy array) – The data cube that is to be smoothed.

• z (float) – The redshift of the coeval cube.

• box_size_mpc (float) – The box size in Mpc. Default value is determined from the box size set for the simulation (set_sim_constants)

• max_baseline (float) – The maximun baseline of the telescope in km. Default value is set as 2 km (SKA core).

• ratio (int) – It is the ratio of smoothing scale in frequency direction and the angular direction (Default value: 1).

• nu_axis (int) – Frequency axis

Returns:

Smoothed_coeval_cube

tools21cm.smoothing.gauss_kernel(size, sigma=1.0, fwhm=None)

Generate a normalized gaussian kernel, defined as exp(-(x^2 + y^2)/(2sigma^2)).

Parameters:

• size (int) – Width of output array in pixels.

• sigma = 1.0 (float) – The sigma parameter for the Gaussian.

• fwhm = None (float or None) – The full width at half maximum. If this parameter is given, it overrides sigma.

Returns:

numpy array with the Gaussian. The dimensions will be size x size or size x sizey depending on whether sizey is set. The Gaussian is normalized so that its integral is 1.

tools21cm.smoothing.gauss_kernel_3d(size, sigma=1.0, fwhm=None)

Generate a normalized gaussian kernel, defined as exp(-(x^2 + y^2 + z^2)/(2sigma^2)).

Parameters:

• size (int) – Width of output array in pixels.

• sigma = 1.0 (float) – The sigma parameter for the Gaussian.

• fwhm = None (float or None) – The full width at half maximum. If this parameter is given, it overrides sigma.

Returns:

numpy array with the Gaussian. The dimensions will be size x size x size. The Gaussian is normalized so that its integral is 1.

tools21cm.smoothing.get_beam_w(baseline, z)

Calculate the width of the beam for an interferometer with a given maximum baseline. It is assumed that observations are done at lambda = 21*(1+z) cm

Parameters:

• baseline (float) – the maximum baseline in meters

• z (float) – the redshift

Returns:

The beam width in arcminutes

tools21cm.smoothing.interpolate2d(input_array, x, y, order=0)

Same as interpolate2d but for 2D data

Parameters:

• input_array (numpy array) – the array to interpolate

• x (numpy array) – the output coordinates along the x axis expressed as (fractional) indices

• y (numpy array) – the output coordinates along the y axis expressed as (fractional) indices

• order (int) – the order of the spline interpolation. Default is 0 (linear interpolation). Setting order=1 gives the same results as IDL’s interpolate function

Returns:

Interpolated array with shape (nx, ny), where nx and ny are the lengths of the arrays x and y respectively.

tools21cm.smoothing.interpolate3d(input_array, x, y, z, order=0)

This function is a recreation of IDL’s interpolate routine. It takes an input array, and interpolates it to a new size, which can be irregularly spaced.

Parameters:

• input_array (numpy array) – the array to interpolate

• x (numpy array) – the output coordinates along the x axis expressed as (fractional) indices

• y (numpy array) – the output coordinates along the y axis expressed as (fractional) indices

• z (numpy array) – the output coordinates along the z axis expressed as (fractional) indices

• order (int) – the order of the spline interpolation. Default is 0 (linear interpolation). Setting order=1 gives the same behaviour as IDL’s interpolate function with default parameters.

Returns:

Interpolated array with shape (nx, ny, nz), where nx, ny and nz are the lengths of the arrays x, y and z respectively.

tools21cm.smoothing.lanczos_kernel(size, kernel_width)

Generate a 2D Lanczos kernel.

Parameters:

• size (int) – the size of the array

• kernel_width (int) – the width of the kernel

Returns:

The kernel as a (size,size) array

tools21cm.smoothing.smooth_coeval(cube, z, box_size_mpc=False, max_baseline=2.0, ratio=1.0, nu_axis=2, verbose=True)

This smooths the coeval cube by Gaussian in angular direction and by tophat along the third axis.

Parameters:

• coeval_cube (numpy array) – The data cube that is to be smoothed.

• z (float) – The redshift of the coeval cube.

• box_size_mpc (float) – The box size in Mpc. Default value is determined from the box size set for the simulation (set_sim_constants)

• max_baseline (float) – The maximun baseline of the telescope in km. Default value is set as 2 km (SKA core).

• ratio (int) – It is the ratio of smoothing scale in frequency direction and the angular direction (Default value: 1).

• nu_axis (int) – Frequency axis

Returns:

Smoothed_coeval_cube

tools21cm.smoothing.smooth_coeval_gauss(cube, fwhm, nu_axis)

This smooths the slices parallel to the given axis of the cube by Gaussian filter.

Parameters:

• cube (numpy array) – The data cube that is to be smoothed.

• fwhm (float) – The fwhm of the Gaussian filter.

• nu_axis (int) – Frequency axis

Returns:

Smoothed_cube

tools21cm.smoothing.smooth_coeval_tophat(cube, width, nu_axis, verbose=True)

This smooths the slices perpendicular to the given axis of the cube by tophat filter.

Parameters:

• cube (numpy array) – The data cube that is to be smoothed.

• width (float) – The width of the tophat filter.

• nu_axis (int) – Frequency axis

Returns:

Smoothed_cube

tools21cm.smoothing.smooth_gauss(input_array, sigma=1.0, fwhm=None)

Smooth the input array with a Gaussian kernel specified either by sigma (standard deviation of the Gaussian function) or FWHM (Full Width Half Maximum). The latter is more appropriate when considering the resolution of a telescope.

Parameters:

• input_array (numpy array) – the array to smooth

• sigma=1.0 (float) – the width of the kernel (variance)

• fwhm = None (float or None) – The full width at half maximum. If this parameter is given, it overrides sigma.

Returns:

The smoothed array. A numpy array with the same dimensions as the input.

tools21cm.smoothing.smooth_lanczos(input_array, kernel_width)

Smooth the input array with a Lanczos kernel.

Parameters:

• input_array (numpy array) – the array to smooth

• kernel_width (int) – the width of the kernel in cells

Returns:

The smoothed array. A numpy array with the same dimensions as the input.

tools21cm.smoothing.smooth_lightcone(lightcone, z_array, box_size_mpc=False, max_baseline=2.0, ratio=1.0, verbose=True)

This smooths in both angular and frequency direction assuming both to be smoothed by same scale.

Parameters:

• lightcone (numpy array) – The lightcone that is to be smoothed.

• z_array (float) – The lowest value of the redshift in the lightcone or the whole redshift array.

• box_size_mpc (float) – The box size in Mpc. Default value is determined from the box size set for the simulation (set_sim_constants)

• max_baseline (float) – The maximun baseline of the telescope in km. Default value is set as 2 km (SKA core).

• ratio (int) – It is the ratio of smoothing scale in frequency direction and the angular direction (Default value: 1).

Returns:

(Smoothed_lightcone, redshifts)

tools21cm.smoothing.smooth_lightcone_gauss(lightcone, fwhm, nu_axis=2)

This smooths the slices perpendicular to the third axis of the lightcone by tophat filter.

Parameters:

• lightcone (numpy array) – The lightcone that is to be smoothed.

• fwhm (numpy array) – The fwhm values of the Gaussian filter at each slice along frequency axis.

• nu_axis (int) – frequency axis

Returns:

Smoothed_lightcone

tools21cm.smoothing.smooth_lightcone_tophat(lightcone, redshifts, dz, verbose=True)

This smooths the slices perpendicular to the third axis of the lightcone by tophat filter.

Parameters:

• lightcone (numpy array) – The lightcone that is to be smoothed.

• redshifts (numpy array) – The redshift of each slice along the third axis.

• dz (float) – redshift width

Returns:

Smoothed_lightcone

tools21cm.smoothing.smooth_lightcone_uv_threshold(lightcone, uv_box, threshold=0.0)

Smoothing (de-noising) the lightcone by removing noisy UV cells below threshold.

Beware, if maximal baseline is to be set, one should clean the uv_box beforehand.

Parameters:

• lightcone (numpy array) – The lightcone to be smoothed.

• uv_box (numpy array) – The UV box.

• threshold (float) – The threshold below which UV cells are removed.

Returns:

Smoothed lightcone.

tools21cm.smoothing.smooth_line(y, window=3, kind='tophat')

Smooths a 1D array using a specified window type.

tools21cm.smoothing.smooth_tophat(input_array, tophat_width)

Smooth the input array with a square tophat kernel.

Parameters:

• input_array (numpy array) – the array to smooth

• tophat_width (int) – the width of the kernel in cells

Returns:

The smoothed array. A numpy array with the same dimensions as the input.

tools21cm.smoothing.smooth_with_kernel(input_array, kernel)

Smooth the input array with an arbitrary kernel.

Parameters:

• input_array (numpy array) – the array to smooth

• kernel (numpy array) – the smoothing kernel. Must be the same size as the input array

Returns:

The smoothed array. A numpy array with the same dimensions as the input.

tools21cm.smoothing.tophat_kernel(size, tophat_width)

Generate a square tophat kernel

Parameters:

• size (int) – the size of the array

• tophat_width (int) – the size of the tophat kernel

Returns:

The kernel as a (size,size) array

tools21cm.smoothing.tophat_kernel_3d(size, tophat_width, shape='cube')

Generate a 3-dimensional tophat kernel with the specified size

Parameters:

• size (integer or list-like) – the size of the tophat kernel along each dimension.

• tophat_width (int) – the size of the tophat kernel

• shape (string) – “cube”: cubic tophat; “sphere”: spherical tophat

Returns:

The normalized kernel

Radio telescope

tools21cm.radio_telescope_layout.antenna_positions_to_baselines(antxyz)

Calculates all possible baseline vectors from a set of antenna positions.

Parameters:

antxyz (np.ndarray or astropy.units.Quantity) – An array (N_ant, 3) of antenna positions.

Returns:

np.ndarray or astropy.units.Quantity – An array (N_ant * (N_ant - 1), 3) of baseline vectors.

tools21cm.radio_telescope_layout.cartesian_to_geographic_coordinate_system(antxyz)

Converts Cartesian coordinates (x, y, z) back to geographic coordinates (longitude, latitude, height above Earth’s surface), using Astropy for unit handling.

Parameters:

antxyz (np.ndarray) – A 2D NumPy array of shape (N, 3), where: - Column 0 is the x-coordinate in meters. - Column 1 is the y-coordinate in meters. - Column 2 is the z-coordinate in meters.

Returns:

antll (np.ndarray) – A 2D NumPy array of shape (N, 3) containing the geographic coordinates: - Column 0: Longitude (λ) in degrees. - Column 1: Latitude (φ) in degrees. - Column 2: Height (h) in meters above the Earth’s surface.

Notes

• The Earth’s radius (mean radius) is assumed to be 6,371 km.

• Astropy is used to ensure unit consistency.

tools21cm.radio_telescope_layout.earth_rotation_effect(Nbase, slice_nums, int_time, total_int_time, declination=-30.0)

Applies Earth’s rotation to baseline vectors for given time slices.

This function simulates the change in projected baseline coordinates (u,v,w) over the course of an observation due to the Earth’s rotation.

Parameters:

• Nbase (np.ndarray) – Array of initial (u,v,w) baseline coordinates.

• slice_nums (int or np.ndarray) – The index or indices of the time slice(s) to compute.

• int_time (float) – Integration time per slice in seconds.

• total_int_time (float) – Total duration of the observation in hours.

• declination (float, optional) – The declination of the pointing center in degrees. Default is -30.

Returns:

new_Nbase (np.ndarray) – The rotated (u,v,w) coordinates. The shape will be (N_baselines, 3) if slice_nums is a single integer, or (N_slices, N_baselines, 3) if slice_nums is an array.

tools21cm.radio_telescope_layout.ecef_to_xyz(ecef, lon_deg, origin_ecef)

Converts ECEF coordinates to a local interferometric XYZ frame.

The XYZ frame is defined with X pointing to local East, Y to local North, and Z aligned with the ECEF Z-axis at the array center.

Parameters:

• ecef (np.ndarray) – Antenna positions (N, 3) in ECEF coordinates (meters).

• lon_deg (float) – Longitude of the array center in degrees.

• origin_ecef (np.ndarray) – ECEF position (3,) of the array center (meters).

Returns:

np.ndarray – Antenna positions (N, 3) in the local XYZ frame (meters).

tools21cm.radio_telescope_layout.enu_to_ecef(enu, lat_deg, lon_deg, origin_ecef)

Converts local East-North-Up (ENU) coordinates to ECEF.

Parameters:

• enu (np.ndarray) – Positions (N, 3) in the local ENU frame (meters).

• lat_deg (float) – Latitude of the local frame’s origin in degrees.

• lon_deg (float) – Longitude of the local frame’s origin in degrees.

• origin_ecef (np.ndarray) – ECEF position (3,) of the local frame’s origin (meters).

Returns:

np.ndarray – Positions (N, 3) in the ECEF frame (meters).

tools21cm.radio_telescope_layout.from_antenna_config(antxyz, z, nu=None)

Calculates baseline vectors (u,v,w) from antenna XYZ positions.

The function converts physical antenna separations into baseline coordinates in units of wavelength.

Parameters:

• antxyz (astropy.units.Quantity or np.ndarray) – The radio telescope antenna configuration in meters.

• z (float) – Redhsift of the slice observed, used to calculate frequency.

• nu (float, optional) – The frequency of observation in MHz. If None, it is calculated from z.

Returns:

• Nbase (np.ndarray) – Numpy array of shape (N_baselines, 3) containing the (u,v,w) values.

• N_ant (int) – The total number of antennas.

tools21cm.radio_telescope_layout.geographic_to_cartesian_coordinate_system(antll)

Converts geographic coordinates (longitude, latitude) into Cartesian coordinates (x, y, z) assuming a spherical Earth with a fixed radius, using Astropy for unit handling.

Parameters:

antll (np.ndarray) – A 2D NumPy array of shape (N, 2), where: - Column 0 is longitude in degrees (λ). - Column 1 is latitude in degrees (φ).

Returns:

antxyz (np.ndarray) – A 2D NumPy array of shape (N, 3) containing the Cartesian coordinates: - Column 0: x-coordinate in meters. - Column 1: y-coordinate in meters. - Column 2: z-coordinate in meters.

Notes

• The Earth’s radius (mean radius) is assumed to be 6,371 km.

• Astropy is used to ensure unit consistency.

tools21cm.radio_telescope_layout.get_SKA_Low_layout(subarray_type='AA4', unit='ecef', verbose=True)

Loads SKA Low antenna layout data from pre-defined configuration files.

Parameters:

• subarray_type (str, optional) – The name of the SKA Low subarray, e.g., “AA4”, “AA2”, “AASTAR”.

• unit (str, optional) – The desired output coordinate system: ‘ecef’, ‘enu’, or ‘xyz’.

• verbose (bool, optional) – If True, prints information about the loaded layout.

Returns:

astropy.units.Quantity – An (N_ant, 3) array of antenna positions in the specified unit.

tools21cm.radio_telescope_layout.get_latest_SKA_Low_layout(subarray_type='AA*', custom_stations=None, add_stations=None, exclude_stations=None, external_telescopes=None)

Retrieves the latest SKA Low layout using the official ska-ost-array-config tool.

This function acts as a wrapper for the SKAO library, providing a convenient way to access up-to-date and standardized array configurations.

Parameters:

• subarray_type (str, optional) – Name of the subarray configuration.

• custom_stations (list, optional) – A list of station names to create a custom subarray.

• add_stations (list, optional) – A list of station names to add to the specified subarray.

• exclude_stations (list, optional) – A list of station names to exclude from the specified subarray.

• external_telescopes (list, optional) – A list of external telescope configurations to include.

Returns:

ska_ost_array_config.array_config.LowSubArray or None – An object representing the SKA Low subarray configuration, or None if the required ska_ost_array_config library is not installed.

tools21cm.radio_telescope_layout.get_uv_daily_observation(ncells, z, antxyz=None, total_int_time=4.0, int_time=10.0, boxsize=None, declination=-30.0, include_mirror_baselines=False, verbose=True)

Simulates a full observation to generate a UV coverage map.

This function combines baseline calculation, Earth rotation simulation, and UV track gridding to produce a final 2D UV coverage map, representing the density of measurements in the Fourier plane.

Parameters:

• ncells (int) – Number of cells in each dimension of the output grid.

• z (float) – Redshift of the observed slice.

• antxyz (np.ndarray, optional) – The radio telescope antenna configuration. If None, a default is loaded.

• total_int_time (float) – Total observation time per day in hours.

• int_time (float) – Integration time per observation in seconds.

• boxsize (float, optional) – Comoving size of the observed sky area in Mpc.

• declination (float) – Declination of the telescope pointing in degrees.

• include_mirror_baselines (bool, optional) – If True, includes mirrored baselines for hermitian conjugation.

• verbose (bool, optional) – If True, prints progress information.

Returns:

• uv_map (np.ndarray) – 2D array (ncells, ncells) of baseline counts per pixel.

• N_ant (int) – The total number of antennas.

tools21cm.radio_telescope_layout.grid_uv_tracks(Nbase, z, ncells, boxsize=None, include_mirror_baselines=False)

Grids baseline (u,v) tracks onto a 2D map for a single time snapshot.

This function takes a set of (u,v,w) coordinates, projects them onto a 2D grid corresponding to the field of view, and counts how many baselines fall into each grid cell.

Parameters:

• Nbase (np.ndarray) – Array containing (u,v,w) values for the antenna baselines.

• z (float) – Redshift of the slice observed.

• ncells (int) – Number of cells in each dimension of the output grid.

• boxsize (float, optional) – Comoving size of the sky observed. Uses a default if None.

• include_mirror_baselines (bool, optional) – If True, also grids the hermitian conjugate baselines at (-u, -v).

Returns:

uv_map (np.ndarray) – A 2D array (ncells, ncells) of baseline counts per pixel, shifted so the zero-frequency is at the center.

tools21cm.radio_telescope_layout.subarray_type_to_antxyz(subarray_type, verbose=True)

Resolves antenna configuration from various input types.

Parameters:

• subarray_type (str or np.ndarray) – If string, loads a pre-defined SKA-Low layout (e.g., “AA4”). If numpy array, it’s used directly as the antenna positions. Otherwise, loads a default SKA1-Low configuration.

• verbose (bool, optional) – If True, enables verbose output during layout loading.

Returns:

• antxyz (astropy.units.Quantity) – A 2D array of antenna (x, y, z) positions with units of meters.

• N_ant (int) – The total number of antennas in the configuration.

tools21cm.radio_telescope_sensitivity.apply_uv_response(array, uv_map)

Parameters:

• array (A complex 2d array of signal in the uv field.)

• uv_map (Numpy array containing the number of baselines observing each pixel.)

Returns:

new_array (It is the ‘array’ after degrading the resoltion with the baseline configuration.)

tools21cm.radio_telescope_sensitivity.get_SEFD(nu_obs, T_sys=None, sefd_data=None, nu_data=None, D_station=40.0, ep_aperture=None, verbose=False)

Calculates the System Equivalent Flux Density (SEFD) for a radio antenna.

This function can operate in two modes: 1. Interpolation: If sefd_data and nu_data are provided, it

interpolates from the given data to find the SEFD at nu_obs.

2. Direct Calculation: If sefd_data is not provided, it calculates

   the SEFD from fundamental parameters like system temperature (T_sys), station diameter (D_station), and aperture efficiency (ep_aperture).

Parameters:

• nu_obs (float or np.ndarray) – Observing frequency or frequencies in MHz.

• T_sys (float, callable, or None, optional) – System temperature in Kelvin. Can be a single float value or a callable function of frequency in MHz. If None, a default model approximating the SKA-Low sky temperature plus a receiver temperature is used. Default is None.

• sefd_data (np.ndarray or None, optional) – An array of known SEFD values for interpolation. If this is provided, nu_data must also be given. Default is None.

• nu_data (np.ndarray or None, optional) – An array of frequencies in MHz corresponding to sefd_data. Default is None.

• D_station (float, optional) – Diameter of the antenna station in meters. Default is 40.0.

• ep_aperture (float, callable, or None, optional) – Aperture efficiency. Can be a single float value or a callable function of frequency in MHz. If None, a default frequency-dependent model is used. Default is None.

• verbose (bool, optional) – If True, enables progress bars and informational messages.

Returns:

sefd (float or np.ndarray) – The calculated System Equivalent Flux Density in Janskys (Jy). Returns a float or array matching the shape of nu_obs.

tools21cm.radio_telescope_sensitivity.jansky_2_kelvin(array, z, boxsize=None, ncells=None)

Parameters:

• array (ndarray) – Numpy array containing the values in jansky.

• z (float) – Redshift

• boxsize (float) – The comoving size of the sky observed. Default: It is determined from the simulation constants set.

• ncells (int) – The number of grid cells. Default: None

Returns:

A numpy array with values in mK.

tools21cm.radio_telescope_sensitivity.kelvin_2_jansky(array, z, boxsize=None, ncells=None)

Parameters:

• array (ndarray) – Numpy array containing the values in mK.

• z (float) – Redshift

• boxsize (float) – The comoving size of the sky observed. Default: It is determined from the simulation constants set.

Returns:

A numpy array with values in Jy.

tools21cm.radio_telescope_sensitivity.kelvin_jansky_conversion(ncells, z, boxsize=None)

Parameters:

• ncells (int) – Number of cells/pixels in the image.

• z (float) – Redshift

• boxsize (float) – The comoving size of the sky observed. Default: It is determined from the simulation constants set.

Returns:

The conversion factor multiplied to values in kelvin to get values in jansky.

tools21cm.radio_telescope_sensitivity.sigma_noise_radio(z, depth_mhz, obs_time, int_time, uv_map=None, N_ant=512, verbose=True, T_sys=None, sefd_data=None, nu_data=None, D_station=40.0, ep_aperture=None)

Calculate the rms of the noise added by radio interferometers.

Parameters:

zfloat

Redshift of the slice observed.

depth_mhzfloat

The bandwidth of the observation (in MHz).

obs_timefloat

The total hours of observation time.

int_timefloat

The integration time.

uv_mapndarray

ncells x ncells numpy array containing the number of baselines observing each pixel.

N_antfloat, optional

Number of antennas in SKA. Default is 564.

verbosebool, optional

If True, print detailed information. Default is True.

T_sys : float, callable, or None, optional System temperature in Kelvin. Can be a single float value or a callable function of frequency in MHz. If None, a default model approximating the SKA-Low sky temperature plus a receiver temperature is used. Default is None.

sefd_datanp.ndarray or None, optional

An array of known SEFD values for interpolation. If this is provided, nu_data must also be given. Default is None.

nu_datanp.ndarray or None, optional

An array of frequencies in MHz corresponding to sefd_data. Default is None.

D_stationfloat, optional

Diameter of the antenna station in meters. Default is 40.0.

ep_aperturefloat, callable, or None, optional

Aperture efficiency. Can be a single float value or a callable function of frequency in MHz. If None, a default frequency-dependent model is used. Default is None.

Returns:

sigmafloat

The rms of the noise in the image produced by SKA for uniformly distributed antennas (in µJy).

rms_noisefloat

The rms of the noise due to the antenna positions in the uv field (in µJy).

Methods:

• simulate the radio telescope observation strategy

• simulate telescope noise

tools21cm.radio_telescope_noise.apply_uv_response_noise(noise, uv_map, boxsize=None, uv_weighting='natural', uv_map_min=0.01, suppress_sharp_features_uv_map=False)

(Internal) Applies the instrumental UV response to a noise array.

This function weights the noise in the Fourier domain according to the UV coverage and the chosen weighting scheme.

Parameters:

• noise (np.ndarray) – The complex noise array in the Fourier domain.

• uv_map (np.ndarray) – The 2D UV coverage map.

• boxsize (float, optional) – The comoving size of the box, needed if smoothing is applied.

• uv_weighting ({‘natural’, ‘uniform’}, optional) – The weighting scheme.

• uv_map_min (float, optional) – Threshold below which UV cells are considered un-sampled and zeroed out.

• suppress_sharp_features_uv_map (bool or str, optional) – Method to smooth the UV map before applying weights.

Returns:

np.ndarray – The noise array after applying the UV response.

tools21cm.radio_telescope_noise.apply_uv_response_noise_briggs(noise, uv_map, robust=2.0, epsilon=1e-06)

Applies Briggs weighting to a noise array.

Parameters:

• noise (np.ndarray) – 2D array of base noise values (assumed uniform before weighting).

• uv_map (np.ndarray) – 2D array of natural weights (hit counts).

• robust (float, optional) – Briggs robustness parameter (-2 for uniform to +2 for near-natural).

• epsilon (float, optional) – A small value to prevent division by zero.

Returns:

np.ndarray – Noise map after applying Briggs weighting.

tools21cm.radio_telescope_noise.apply_uv_response_on_coeval(array, z, subarray_type='AA4', boxsize=None, total_int_time=6.0, int_time=10.0, declination=-30.0, uv_map=None, N_ant=None)

Applies the UV response slice-by-slice to a 3D coeval data cube.

Parameters:

• array (np.ndarray) – The input 3D data cube (nx, ny, nz).

• z, subarray_type, boxsize, etc. (various, optional) – Parameters passed to get_uv_map if uv_map is not provided.

Returns:

np.ndarray – The observed 3D data cube after applying the UV mask to each slice.

tools21cm.radio_telescope_noise.apply_uv_response_on_image(array, uv_map)

Simulates the effect of incomplete UV coverage on a true sky image.

This function mimics the observation process by transforming an image to the Fourier (UV) domain, setting all un-sampled UV cells to zero, and transforming back to the image domain.

Parameters:

• array (np.ndarray) – The 2D input image array.

• uv_map (np.ndarray) – The 2D UV coverage map, where non-zero values indicate sampled cells.

Returns:

np.ndarray – The resulting “dirty” image after applying the UV mask.

tools21cm.radio_telescope_noise.gauss_kernel_3d(size, sigma=1.0, fwhm=None)

Generates a 3D normalized Gaussian kernel.

Parameters:

• size (int) – Width of the output array in pixels.

• sigma (float, optional) – The standard deviation of the Gaussian.

• fwhm (float, optional) – The Full Width at Half Maximum. Overrides sigma if provided.

Returns:

np.ndarray – A 3D numpy array with the Gaussian kernel, normalized to sum to 1.

tools21cm.radio_telescope_noise.get_uv_map(ncells, z, subarray_type='AA4', total_int_time=6.0, int_time=10.0, boxsize=None, declination=-30.0, include_mirror_baselines=True, verbose=True, max_baseline=None)

Generates a 2D UV coverage map for a given observation.

This function simulates a daily observation to calculate the UV coverage, optionally applying a cut based on maximum baseline length.

Parameters:

• ncells (int) – Number of cells in the grid.

• z (float) – Redshift of observation.

• subarray_type (str, optional) – The telescope layout configuration name.

• total_int_time, int_time, declination (float, optional) – Observation parameters.

• boxsize (float, optional) – Comoving size of the simulation box in Mpc.

• include_mirror_baselines (bool, optional) – If True, includes mirrored baselines.

• verbose (bool, optional) – If True, enables verbose output.

• max_baseline (float, optional) – The maximum baseline length in km. If provided, UV tracks beyond the corresponding k-mode are cut.

Returns:

• uv_map (np.ndarray) – The 2D array of gridded UV coverage.

• N_ant (int) – The number of antennas.

tools21cm.radio_telescope_noise.get_uv_map_lightcone(ncells, zs, subarray_type='AA4', total_int_time=6.0, int_time=10.0, boxsize=None, declination=-30.0, save_uvmap=None, n_jobs=4, verbose=True, checkpoint=16)

Generates or loads a lightcone of UV coverage maps, one for each redshift.

This function handles the creation of UV maps across a range of redshifts, using parallel processing for efficiency and providing an option to cache the results to a file to speed up subsequent runs.

Parameters:

• ncells (int) – The number of grid cells.

• zs (np.ndarray or list) – An array or list of redshift values for the lightcone.

• subarray_type, total_int_time, etc. (various, optional) – Observational parameters passed to get_uv_map.

• save_uvmap (str, optional) – File path to save or load the UV map dictionary. If the file exists, maps are loaded; otherwise, they are generated and saved.

• n_jobs (int, optional) – Number of CPUs for parallel generation of UV maps.

• verbose (bool, optional) – If True, enables progress bars and informational messages.

• checkpoint (int, optional) – If provided, the number of redshifts to process before saving the results to save_uvmap. This is useful for long runs to prevent data loss.

Returns:

dict – A dictionary containing the UV maps for each redshift, keyed by the redshift value formatted to three decimal places, along with metadata about the simulation parameters.

tools21cm.radio_telescope_noise.ifft2_wrap(nn1)

Performs a 2D inverse FFT with wrapping to handle boundaries.

This is useful for avoiding artifacts when the data is not perfectly periodic.

Parameters:

nn1 (np.ndarray) – The 2D array in Fourier space.

Returns:

np.ndarray – The central part of the inverse transformed array.

tools21cm.radio_telescope_noise.make_uv_map_lightcone(ncells, zs, subarray_type='AA4', total_int_time=6.0, int_time=10.0, boxsize=None, declination=-30.0, verbose=True)

Creates a 3D cube of UV maps, one for each redshift slice.

Parameters:

• ncells (int) – Number of cells in the grid.

• zs (np.ndarray) – An array of redshift values for the lightcone.

• subarray_type, total_int_time, int_time, declination, boxsize (various, optional) – Observational parameters passed to get_uv_map.

• verbose (bool, optional) – If True, enables verbose output.

Returns:

• uv_lc (np.ndarray) – A 3D array (ncells, ncells, n_redshifts) containing the UV maps.

• N_ant (int) – The number of antennas.

tools21cm.radio_telescope_noise.max_baseline_to_max_k(redshift, max_baseline)

Converts a maximum physical baseline length to a maximum k-mode.

Parameters:

• redshift (float) – The redshift of observation.

• max_baseline (float or astropy.units.Quantity) – The maximum baseline length of the telescope, assumed in km if no units are given.

Returns:

astropy.units.Quantity – The corresponding maximum k-mode (spatial frequency).

tools21cm.radio_telescope_noise.noise_coeval_power_spectrum_1d(ncells, z, depth_mhz, obs_time=1000, subarray_type='AA4', kbins=10, boxsize=None, binning='log', return_n_modes=False, total_int_time=6.0, int_time=10.0, declination=-30.0, uv_map=None, N_ant=None, uv_weighting='natural', T_sys=None, sefd_data=None, nu_data=None, D_station=40.0, ep_aperture=None, fft_wrap=False, verbose=True)

Computes the 1D power spectrum of instrumental noise for a coeval observation.

This function simulates the instrumental noise properties in the Fourier domain based on the telescope configuration and observation strategy, and then computes the radially averaged 1D power spectrum.

Parameters:

ncellsint

The number of grid cells in each spatial dimension.

zfloat

The redshift of the observation.

depth_mhzfloat

The observational bandwidth in MHz.

obs_timefloat, optional

The total on-source observation time in hours.

subarray_typestr, optional

The telescope layout configuration name.

kbinsint, optional

The number of bins for radial averaging of the power spectrum.

boxsizefloat, optional

The comoving size of the simulation box in Mpc.

binning{‘log’, ‘linear’}, optional

The type of binning to use for k-modes.

return_n_modesbool, optional

If True, also returns the number of modes in each k-bin.

total_int_timefloat, optional

Total observation time per day in hours for UV coverage simulation.

int_timefloat, optional

Integration time in seconds for UV coverage simulation.

declinationfloat, optional

The pointing declination in degrees.

uv_mapnp.ndarray, optional

A pre-computed UV coverage map. If None, it will be simulated.

N_antint, optional

The number of antennas. If None, it will be determined.

uv_weighting{‘natural’, ‘uniform’}, optional

The weighting scheme to apply in the UV plane.

T_sys : float, callable, or None, optional System temperature in Kelvin. Can be a single float value or a callable function of frequency in MHz. If None, a default model approximating the SKA-Low sky temperature plus a receiver temperature is used. Default is None.

sefd_datanp.ndarray or None, optional

An array of known SEFD values for interpolation. If this is provided, nu_data must also be given. Default is None.

nu_datanp.ndarray or None, optional

An array of frequencies in MHz corresponding to sefd_data. Default is None.

D_stationfloat, optional

Diameter of the antenna station in meters. Default is 40.0.

ep_aperturefloat, callable, or None, optional

Aperture efficiency. Can be a single float value or a callable function of frequency in MHz. If None, a default frequency-dependent model is used. Default is None. verbose : bool, optional

If True, enables verbose output.

Returns:

pnnp.ndarray

The 1D noise power spectrum values.

knnp.ndarray

The central k-mode values for each bin.

n_modesnp.ndarray, optional

The number of modes in each k-bin (if return_n_modes is True).

tools21cm.radio_telescope_noise.noise_cube_coeval(ncells, z, depth_mhz=None, obs_time=1000, subarray_type='AA4', boxsize=None, total_int_time=6.0, int_time=10.0, declination=-30.0, uv_map=None, N_ant=None, uv_weighting='natural', verbose=True, fft_wrap=False, T_sys=None, sefd_data=None, nu_data=None, D_station=40.0, ep_aperture=None, suppress_sharp_features_uv_map=False)

Generates a 3D coeval cube of instrumental noise.

This function simulates a noise cube where the noise properties are constant along the line-of-sight (frequency) axis.

Parameters:

• ncells, z, obs_time, etc. (various) – Parameters passed to noise_map for each slice.

• depth_mhz (float, optional) – The bandwidth per channel. If None, it is calculated from boxsize.

Returns:

np.ndarray – A 3D cube (ncells, ncells, ncells) of instrumental noise in mK.

tools21cm.radio_telescope_noise.noise_cube_lightcone(ncells, z, obs_time=1000, subarray_type='AA4', boxsize=None, save_uvmap=None, total_int_time=6.0, int_time=10.0, declination=-30.0, N_ant=None, uv_weighting='natural', fft_wrap=False, verbose=True, n_jobs=4, checkpoint=64, T_sys=None, sefd_data=None, nu_data=None, D_station=40.0, ep_aperture=None, suppress_sharp_features_uv_map=False)

Generates a 3D lightcone of instrumental noise around a central redshift.

This function is ideal for single redshift studies with a narrow bandwidth, where each slice along the line-of-sight corresponds to a slightly different redshift and frequency.

Parameters:

• ncells (int) – The grid size.

• z (float) – The central redshift of the lightcone.

• obs_time, subarray_type, etc. (various) – Observational parameters.

• save_uvmap (str, optional) – File path to save or load pre-computed UV maps to speed up re-runs.

• n_jobs, checkpoint (int, optional) – Parameters for parallel processing of UV map generation.

Returns:

np.ndarray – A 3D (ncells, ncells, ncells) lightcone of instrumental noise in mK.

tools21cm.radio_telescope_noise.noise_lightcone(ncells, zs, obs_time=1000, subarray_type='AA4', boxsize=None, save_uvmap=None, total_int_time=6.0, int_time=10.0, declination=-30.0, uv_weighting='natural', fft_wrap=False, verbose=True, n_jobs=4, checkpoint=16, T_sys=None, sefd_data=None, nu_data=None, D_station=40.0, ep_aperture=None, suppress_sharp_features_uv_map=False)

Generates a 3D lightcone of instrumental noise over a list of redshifts.

Each slice along the line-of-sight corresponds to a different redshift, and noise properties are calculated accordingly. This function first calls get_uv_map_lightcone to efficiently generate or load all required UV maps.

Parameters:

• ncells (int) – The number of grid cells along each spatial dimension.

• zs (np.ndarray or list) – An array of redshift values defining the slices of the lightcone.

• obs_time (float, optional) – Total observation time in hours. Default is 1000.

• subarray_type (str, optional) – The type of telescope subarray (e.g., “AA4”). Default is “AA4”.

• boxsize (float, optional) – The comoving size of the simulation box in Mpc. If None, a default cosmological value is used.

• save_uvmap (str, optional) – File path to save or load the UV coverage maps. Caching these maps can significantly speed up subsequent runs.

• total_int_time (float, optional) – Total integration time in hours used for generating the UV coverage. Default is 6.

• int_time (float, optional) – The integration time for a single visibility measurement in seconds. Default is 10.

• declination (float, optional) – The pointing declination of the telescope in degrees. Default is -30.

• uv_weighting (str, optional) – The UV weighting scheme to use (e.g., ‘natural’, ‘uniform’). Default is ‘natural’.

• fft_wrap (bool, optional) – If True, use pyfft_wrap for FFTs, which can be faster. Default is False.

• verbose (bool, optional) – If True, print progress bars and status messages. Default is True.

• n_jobs (int, optional) – The number of CPU cores to use for parallel generation of UV maps. Default is 4.

• checkpoint (int, optional) – Number of redshift slices to process before saving the UV maps to a checkpoint file. Useful for long runs. Default is 16.

• sefd_data (str, optional) – Path to a file containing System Equivalent Flux Density (SEFD) data.

• nu_data (str, optional) – Path to a file containing frequencies corresponding to the SEFD data.

• suppress_sharp_features_uv_map (bool, optional) – If True, applies a suppression filter to the UV map to mitigate sharp features. Default is False.

Returns:

np.ndarray – A 3D array of shape (ncells, ncells, len(zs)) representing the noise lightcone in units of mK.

See also

get_uv_map_lightcone

Generates the underlying UV coverage maps.

noise_map

Generates a 2D noise map for a single redshift.

tools21cm.radio_telescope_noise.noise_map(ncells, z, depth_mhz, obs_time=1000, subarray_type='AA4', boxsize=None, total_int_time=6.0, int_time=10.0, declination=-30.0, uv_map=None, N_ant=None, uv_weighting='natural', T_sys=None, sefd_data=None, nu_data=None, D_station=40.0, ep_aperture=None, fft_wrap=False, verbose=True, suppress_sharp_features_uv_map=False)

Creates a 2D map of instrumental noise for a single observation slice.

This function simulates the thermal noise from a radio interferometer, taking into account the UV coverage and weighting scheme, and produces a noise map in the image domain.

Parameters:

ncellsint

The number of grid cells in each spatial dimension.

zfloat

The redshift of the observation.

depth_mhzfloat

The observational bandwidth in MHz.

obs_timefloat, optional

The total on-source observation time in hours.

subarray_typestr, optional

The telescope layout configuration name.

boxsizefloat, optional

The comoving size of the simulation box in Mpc.

total_int_time, int_time, declinationfloat, optional

Parameters for get_uv_map if uv_map is not provided.

uv_map, N_antvarious, optional

Pre-computed UV map and number of antennas.

uv_weighting{‘natural’, ‘uniform’}, optional

The weighting scheme to apply.

T_sys : float, callable, or None, optional System temperature in Kelvin. Can be a single float value or a callable function of frequency in MHz. If None, a default model approximating the SKA-Low sky temperature plus a receiver temperature is used. Default is None.

sefd_datanp.ndarray or None, optional

An array of known SEFD values for interpolation. If this is provided, nu_data must also be given. Default is None.

nu_datanp.ndarray or None, optional

An array of frequencies in MHz corresponding to sefd_data. Default is None.

D_stationfloat, optional

Diameter of the antenna station in meters. Default is 40.0.

ep_aperturefloat, callable, or None, optional

Aperture efficiency. Can be a single float value or a callable function of frequency in MHz. If None, a default frequency-dependent model is used. Default is None. fft_wrap : bool, optional

If True, use a wrapped FFT to handle periodic boundary conditions.

verbosebool, optional

If True, enables verbose output.

suppress_sharp_features_uv_mapbool or str, optional

If a string, specifies a method to smooth the UV map to reduce artifacts.

Returns:

np.ndarray

A 2D array representing the noise map in the image domain (in muJy).

tools21cm.radio_telescope_noise.signal_window(ncells, method, ndim=1, extra_param=30)

Generates a 1D, 2D, or 3D window function for signal processing.

Parameters:

• ncells (int) – The size of the window along each dimension.

• method (str) – The name of the window function to use (e.g., ‘blackmanharris’, ‘gaussian’).

• ndim (int, optional) – The number of dimensions for the window (1, 2, or 3).

• extra_param (float, optional) – An additional parameter required by some window functions, like the standard deviation for ‘gaussian’ or beta for ‘kaiser’.

Returns:

np.ndarray or None – The generated window array, or None if the method is not implemented or the dimension is not supported.

tools21cm.radio_telescope_noise.smooth_gauss_3d(array, fwhm)

Smooths a 3D array using a Gaussian kernel via FFT convolution.

Parameters:

• array (np.ndarray) – The 3D input array to be smoothed.

• fwhm (float) – The Full Width at Half Maximum of the Gaussian kernel in pixels.

Returns:

np.ndarray – The smoothed 3D array.

Temperature

Methods to estimate the brightness temperature.

tools21cm.temperature.calc_dt(xfrac, dens, z=-1)

Calculate the differential brightness temperature assuming T_s >> T_CMB

Parameters:

• xfrac (XfracFile object, string or numpy array) – the ionization fraction

• dens (DensityFile object, string or numpy array) – density in cgs units

• z = -1 (float) – The redshift (if < 0 this will be figured out from the files)

Returns:

The differential brightness temperature (in mK) as a numpy array with the same dimensions as xfrac.

tools21cm.temperature.calc_dt_full(xfrac, dens, temp, z=-1, correct=True)

Calculate the differential brightness temperature assuming only that Lyman alpha is fully coupled so T_s = T_k

(NOT T_s >> T_CMB)

Parameters:

• xfrac (XfracFile object, string or numpy array) – the ionization fraction

• dens (DensityFile object, string or numpy array) – density in cgs units

• temp (TemperFile object, string or numpy array) – the temperature in K

• z = -1 (float) – The redshift (if < 0 this will be figured out from the files)

• correct = True (bool) – if true include a correction for partially ionized cells.

Returns:

The differential brightness temperature (in mK) as a numpy array with the same dimensions as xfrac.

tools21cm.temperature.calc_dt_full_lightcone(xfrac, temp, dens, lowest_z, los_axis=2, correct=True)

Calculate the differential brightness temperature assuming only that Lyman alpha is fully coupled so T_s = T_k

(NOT T_s >> T_CMB) for lightcone data. UNTESTED

Parameters:

• xfrac (string or numpy array) – the name of the ionization fraction file (must be cbin), or the xfrac lightcone data

• temp (string or numpy array) – the name of the temperature file (must be cbin), or the temp lightcone data

• dens (string or numpy array) – the name of the density file (must be cbin), or the density data

• lowest_z (float) – the lowest redshift of the lightcone volume

• los_axis = 2 (int) – the line-of-sight axis

• correct = True (bool) – if true include a correction for partially ionized cells.

Returns:

The differential brightness temperature (in mK) as a numpy array with the same dimensions as xfrac.

tools21cm.temperature.calc_dt_halo(mhalo, box_dim, z, Y_p=0.249, cosmo=None)

Calculate the differential brightness temperature assuming T_s >> T_CMB

Parameters:

• mhalo (numpy array) – the ionization fraction

• box_dim (float) – Length of the simulation box in each direction.

• z (float) – The redshift.

• Y_p = 0.249 (float) – Helium abundance mass factor.

• cosmo – astropy.cosmology.Cosmology, optional

Cosmology object. If None, assumes Planck18 cosmology.

Returns:

The differential brightness temperature (in mK) as a numpy array with the same dimensions as xfrac.

tools21cm.temperature.calc_dt_lightcone(xfrac, dens, lowest_z, los_axis=2)

Calculate the differential brightness temperature assuming T_s >> T_CMB for lightcone data.

Parameters:

• xfrac (string or numpy array) – the name of the ionization fraction file (must be cbin), or the xfrac lightcone data

• dens (string or numpy array) – the name of the density file (must be cbin), or the density data

• lowest_z (float) – the lowest redshift of the lightcone volume

• los_axis = 2 (int) – the line-of-sight axis

Returns:

The differential brightness temperature (in mK) as a numpy array with the same dimensions as xfrac.

tools21cm.temperature.mean_dt(z)

Get the mean dT at redshift z

Parameters:

z (float or numpy array) – the redshift

Returns:

dT (float or numpy array) the mean brightness temperature in mK

Thompson Optical Depth

tools21cm.tau.tau(ionfractions, redshifts, num_points=50)

Calculate the optical depth to Thomson scattering.

Parameters:

• ionfractions (ndarray) – An array containing the ionized fraction at various points along the line-of-sight.

• redshifts (ndarray) – An array containing the redshifts corresponding to the same points as ionfractions. Must be the same length as ionfractions.

• num_points (int, optional) – The number of initial points used for the integration to account for high-redshift contributions where ionization is assumed to be negligible. Default is 50.

Returns:

tuple – A tuple containing:

tau_0ndarray

Optical depth at each redshift. The length is equal to the length of redshifts + num_points.

tau_zndarray

The corresponding redshift array. The length is equal to the length of redshifts + num_points.

Notes

• The Universe is assumed to be fully ionized at the lowest redshift supplied.

• To get the total optical depth, refer to the last value in tau_0.

tools21cm.tau.tau_map(ionfractions, redshifts=None, num_points=50, reading_function=None)

Calculate the optical depth to Thomson scattering for lightcones or maps.

Parameters:

• ionfractions (ndarray or dict) – Ionized fraction data across various points along the line-of-sight.

• redshifts (ndarray, optional) – Array of redshift values corresponding to the ionfractions data. Must be the same length as ionfractions if ionfractions is an ndarray. If ionfractions is a dict and redshifts is None, redshifts are inferred from the keys of the dictionary.

• num_points (int, optional) – Number of initial points used for the integration to account for high-redshift contributions where ionization is assumed to be negligible. Default is 50.

• reading_function (callable, optional) – Custom function to read and process ionization fraction data when ionfractions is provided in a format that requires specialized reading. The function should take a filename or data identifier as input and return an ndarray of ionization fractions.

Returns:

tuple – A tuple containing:

tau_0ndarray

Optical depth values at each spatial position and redshift. The shape is (N_x, N_y, N_z + num_points), where N_x and N_y are spatial dimensions, and N_z is the number of redshift slices in output_z.

tau_zndarray

Array of redshift values corresponding to each slice in tau_0. The length is N_z + num_points.

Notes

• The Universe is assumed to be fully ionized at the lowest redshift supplied.

• The reading_function is useful when working with custom data formats.

Topology

tools21cm.topology.EulerCharacteristic(data, thres=0.5, neighbors=6, speed_up='cython', verbose=True, n_jobs=None)

Parameters:

• data (ndarray) – The data cube containing the structure.

• thres (float) – The threshold to create the binary field from the data (Default: 0.5) Ignore this parameter if data is already a binary field.

• neighbors (int) – Define the connectivity to the neighbors (Default: 6).

• speed_up (str) – Method used to speed up calculation (Default: cython). Options are cython, numba, torch and numpy. The calculation with torch uses GPUs if available on the device.

• verbose (bool) – If True, verbose is printed.

Returns:

Euler characteristics value.

tools21cm.topology.betti0(data, thres=0.5, neighbors=6)

Parameters:

• data (ndarray) – The data cube containing the structure.

• thres (float) – The threshold to create the binary field from the data (Default: 0.5) Ignore this parameter if data is already a binary field.

• neighbors (int) – Define the connectivity to the neighbors (Default: 6).

Returns:

Betti 0

tools21cm.topology.betti1(data, thres=0.5, neighbors=6, b0=None, b2=None, chi=None, speed_up='cython', verbose=True, n_jobs=None)

Parameters:

• data (ndarray) – The data cube containing the structure.

• thres (float) – The threshold to create the binary field from the data (Default: 0.5) Ignore this parameter if data is already a binary field.

• neighbors (int) – Define the connectivity to the neighbors (Default: 6).

• speed_up (str) – Method used to speed up calculation (Default: cython).

• verbose (bool) – If True, verbose is printed.

Returns:

Betti 1

tools21cm.topology.betti2(data, thres=0.5, neighbors=6)

Parameters:

• data (ndarray) – The data cube containing the structure.

• thres (float) – The threshold to create the binary field from the data (Default: 0.5) Ignore this parameter if data is already a binary field.

• neighbors (int) – Define the connectivity to the neighbors (Default: 6).

Returns:

Betti 2

tools21cm.topology.genus(data, xth=0.5)

@ Chen & Rong (2010)

It is an implementation of Gauss-Bonnet theorem in digital spaces.

Parameters:

• data (ndarray) – The input data set which contains the structures.

• xth (float) – The threshold to put on the data set (Default: 0.5).

Returns:

Genus