gridded intensity mapping power
Toolkit for regridding line intensity maps onto a Cartesian grid on which fast Fourier transforms can be run for analysing n-point clustering statistics (such as the power spectrum) in Fourier space. The pipeline provides two novel features:
- a validated Monte-Carlo style resampling from map (R.A., Dec. Frequency) space onto a Cartesian (x, y, z in [Mpc/h]) grid,
- option of higher-order interpolation of the sampling particles (using
pmesh
), which mitigates aliasing effects.
The algorithms in this package are discussed in detail in the accompanying paper (https://arxiv.org/pdf/2312.07289.pdf), where validation on a suite of simulations is presented.
- Clone this repo
git clone https://github.com/stevecunnington/gridimp
. - Change directory to the cloned gridimp repo.
cd gridimp
- Activate a python environment with the required packages (see full list of dependencies in the
environments.yml
file). The simplest way to create and activate a new anaconda environment with the necessary packages is with:
conda env create -f environment.yml
conda activate gridimp
- Install the repo with:
python setup.py install
- Validate the package installation by running the test demo script:
python examples/toy.py
Generate a fast cubic lognormal mock over n=128^3 cells with size l=500^3(Mpc/h)^3:
from gridimp import cosmo
from gridimp import mock
cosmo.SetCosmology(z=1) # set cosmology
Pmod = cosmo.GetModelPk(z=1) # obtain matter power spectrum
l = 500 # length of box down one side [Mpc/h]
n = 128 # number of cells down one side
dims = [l,l,l,n,n,n,0,0,0] # lengths,cells,origins for each dimension
delta_0 = mock.Generate(Pmod,dims) # generate LIM mock (default b=1,T=1)
By defining input simulated survey parameters, a lognormal mock can be generated onto a Cartesian grid calculated to fully enclose the survey footprint.
### Define survey parameters:
ramin,ramax = 10,30 # [deg]
decmin,decmax = 10,30 # [deg]
numin,numax = 900,1110 # [MHz]
nnu = 60 # number of frequency channels
nside = 128 # HEALPix resolution: determines pixel size
n0 = 256 # n0^3 will be number of cells for input grid cube
### Initialise healpix and get sky coordinates for voxels covering survey:
from gridimp import grid
ra,dec,nu,dims_0 = grid.init_healpix(nside,ramin,ramax,decmin,decmax,numin,numax,nnu,n0)
delta_0 = mock.Generate(Pmod,dims_0) # generate LIM mock (default b=1,T=1)
# Create healpy sky maps for each frequency channel:
m = grid.lightcone_healpy(delta_0,dims_0,ra,dec,nu,nside)
n = 64 # n^3 will be number of cells for out FFT grid
Np = 5 # number of sampling particles per map voxel for regridding
window = 'ngp' # mass assignment function
compensate = True # correct for interpolaton effects at field level
interlace = False # interlace FFT fields
# Regrid to FFT:
from astropy import units as u
from gridimp import line
ra_p,dec_p,nu_p,pixvals = grid.SkyPixelParticles_healpy(ra,dec,nu,nside,m,Np=Np)
xp,yp,zp = grid.SkyCoordtoCartesian(ra_p.to(u.deg).value,dec_p.to(u.deg).value,line.nu21cm_to_z(nu_p),ramean_arr=ra.to(u.deg).value,decmean_arr=dec.to(u.deg).value,doTile=False)
dims_fft = grid.comoving_dims(ra,dec,nu,nside,(n,n,n))
delta_G,W,counts = grid.mesh(xp,yp,zp,pixvals,dims_fft,window,compensate,interlace)
examples/toy.py
has a walkthrough of simulation generation -> Healpix sky mapping -> Cartesian regridding -> power spectrum estimation.
scripts
folder has the original in-depth code used to generate the main results in the accompanying paper (https://arxiv.org/pdf/2312.07289.pdf).
If you use any of the code in gridimp, we kindly ask that you cite the accompanying paper (https://arxiv.org/pdf/2312.07289.pdf).
Feel free to get in contact with any comments, suggestions, or bug reports at:
Steve Cunnington [[email protected]],
alternatively, please open a New issue from the Issues tab.