Adding observational effects in your model power spectrum

You should first check why is it important to have the correct survey dimension. Then, hopefully now you are already familiar with creating a survey volume and corresponding k-modes:

import numpy as np
import matplotlib.pyplot as plt
from meer21cm import util
from meer21cm.plot import plot_map
from meer21cm import PowerSpectrum

z_min = 0.5
z_max = 1.5
# frequency resolution in Hz, use your own number
freq_resol = 0.5e6
num_ch = int((util.redshift_to_freq(z_min)-util.redshift_to_freq(z_max))/freq_resol)
# obtain the frequency channels
nu = np.linspace(0,num_ch-1,num_ch) * freq_resol + util.redshift_to_freq(z_max)

ra_range = [-45, 80]
dec_range = [-70, 5]
ang_resol = 1.0 # in deg, similar to beam size here for simplicity
num_pix_x = 102 # one extra on each side as buffer
num_pix_y = 102 # one extra on each side as buffer
wcs, num_pix_x, num_pix_y = util.create_wcs_with_range(
    ra_range,
    dec_range,
    [ang_resol,ang_resol],
)

ps = PowerSpectrum(
    wproj=wcs,
    num_pix_x=num_pix_x,
    num_pix_y=num_pix_y,
    nu=nu,
    ra_range=ra_range,
    dec_range=dec_range,
    # downscale the resolution along line-of-sight
    downres_factor_radial = 1.5,
    # downscale the resolution on the transverse plane
    downres_factor_transverse = 1.2,
    omega_hi = 5e-4,
    tracer_bias_1 = 1.5,
    mean_amp_1 = 'average_hi_temp'
)

ps.get_enclosing_box()

/Users/zhaotingchen/miniconda3/envs/meer21cm/lib/python3.10/site-packages/halomod/halo_model.py:32: UserWarning: Warning: Some Halo-Exclusion models have significant speedup when using Numba
  from .halo_exclusion import Exclusion, NoExclusion

Now we have the realistic k-modes. In principle, you can then use your own modelling to generate the 3D power spectrum, by inputting ps.kmode and ps.mumode. For now, let us just use the simple linearly biased HI power spectrum with a kaiser RSD effect:

ps.k1dbins = np.geomspace(2e-3, 0.1, 11)
ps.kperpbins = np.linspace(0, 0.066, 12)
ps.kparabins = np.linspace(0, 0.1, 21)
pscy,_ = ps.get_cy_power(ps.auto_power_tracer_1_model)
ps1d,keff,_ = ps.get_1d_power(ps.auto_power_tracer_1_model)

plt.pcolormesh(ps.kperpbins, ps.kparabins, np.log10(pscy).T)
plt.colorbar()
plt.xlabel(r'k$_\perp$')
plt.ylabel(r'k$_\parallel$')

Text(0, 0.5, 'k$_\\parallel$')

plt.plot(keff, ps1d)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'k')
plt.ylabel(r'P(k)')

Text(0, 0.5, 'P(k)')

Beam attenuation

The first thing you can add is the beam. This is taken care of by adding a Gaussian attenuation, with the beam width (sigma, not FWHM) specified by ps.sigma_beam_ch:

ps.sigma_beam_ch = 0.5 # in deg
ps.sigma_beam_in_mpc

np.float64(27.242727548225005)

You can put in a frequency dependent beam as well by passing an array to ps.sigma_beam_ch. It matter only in mock simulation. At modelling level the only thing that matters is the ps.sigma_beam_in_mpc (which of course is a simplifcation that may be inaccurate). Then you can retrieve the beam attenuation term:

ps.beam_attenuation().shape

(87, 119, 378)

You can then simply multiply this term squared (in auto power spectrum) to your model. If you use the simple model provided in meer21cm, then this is automatically applied when you set sigma_beam_ch and include_beam;

ps.include_beam = [True, False]# first element is HI

ps1d_beam,keff,_ = ps.get_1d_power(ps.auto_power_tracer_1_model)
# this is the same as manually applying average temperature and beam
ps1d_beam2,keff,_ = ps.get_1d_power(
    ps.average_hi_temp**2 * ps.auto_power_tracer_1_model_noobs*ps.beam_attenuation()**2
)

plt.plot(keff, ps1d)
plt.plot(keff, ps1d_beam)
plt.plot(keff, ps1d_beam2,ls=':')

plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'k')
plt.ylabel(r'P(k)')

Text(0, 0.5, 'P(k)')

Map resolution

HI intensity maps are averaged into sky pixels and frequency channels. This leads to a window function effect in the power spectrum. In meer21cm, you can retrieve this effect by first setting the sampling resolution of the map cube:

ps.sampling_resol = [
    ps.pix_resol_in_mpc,
    ps.pix_resol_in_mpc,
    ps.los_resol_in_mpc,
]

ps.map_sampling().shape

(87, 119, 378)

Again, you can multiply this term squared to the power spectrum. In meer21cm, this is automatically included if you turn on ps.include_sky_sampling:

ps.include_sky_sampling = [True, False]

ps1d_beam_map,keff,_ = ps.get_1d_power(ps.auto_power_tracer_1_model)
# this is the same as manually applying average temperature and beam
ps1d_beam_map2,keff,_ = ps.get_1d_power(
    ps.average_hi_temp**2
    * ps.auto_power_tracer_1_model_noobs
    * ps.beam_attenuation()**2
    * ps.map_sampling()**2
)

plt.plot(keff, ps1d)
plt.plot(keff, ps1d_beam_map)
plt.plot(keff, ps1d_beam_map2,ls=':')

plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'k')
plt.ylabel(r'P(k)')

Text(0, 0.5, 'P(k)')

Gridding compensation

HI intensity maps, as well as galaxy catalogues, are gridded onto the rectangular grid. This creates an additional windowing effect that attenuates the power spectrum. You can add in this effect by having set the box dimensions using ps.get_enclosing_box(), specify the mass assignment scheme, and then turn on ps.compensate:

ps.compensate = [True, True]
ps.grid_scheme = 'cic' # [nnb, cic, tsc, pcs]

ps.gridding_compensation().shape

(87, 119, 378)

Again, you can multiply this term squared to your power spectrum:

ps1d_beam_map_comp,keff,_ = ps.get_1d_power(ps.auto_power_tracer_1_model)
# this is the same as manually applying average temperature and beam
ps1d_beam_map_comp2,keff,_ = ps.get_1d_power(
    ps.average_hi_temp**2
    * ps.auto_power_tracer_1_model_noobs
    * ps.beam_attenuation()**2
    * ps.map_sampling()**2
    * ps.gridding_compensation()**2
)

plt.plot(keff, ps1d)
plt.plot(keff, ps1d_beam_map_comp)
plt.plot(keff, ps1d_beam_map_comp2,ls=':')

plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'k')
plt.ylabel(r'P(k)')

Text(0, 0.5, 'P(k)')

Weights convolution

Finally, the gridded box is not uniformly weighted when performing Fourier transformation. The weighting effectively creates a convolution in Fourier space with the power spectrum.

In meer21cm, you can easily retrieve the box grid sampling by performing a gridding (note this is just to invoke the calculation, you don’t need any data in):

ps.grid_data_to_field();

fig, axes = plt.subplots(1,3,figsize=(15,5))
for i, ax in enumerate(axes):
    ax.imshow(ps.counts_in_box.mean(i),aspect='auto')

The weighting is then the counts of pixels in the box grids. In meer21cm you can set the associated weights to toggle the convolution, which under the hood uses get_modelpk_conv:

from meer21cm.power import get_modelpk_conv

ps.weights_1 = ps.counts_in_box.astype('float')

ps1d_beam_map_comp_w,keff,_ = ps.get_1d_power(ps.auto_power_tracer_1_model)
# this is the same as manually applying average temperature and beam
ps1d_beam_map_comp_w2,keff,_ = ps.get_1d_power(
    get_modelpk_conv(ps.average_hi_temp**2
    * ps.auto_power_tracer_1_model_noobs
    * ps.beam_attenuation()**2
    * ps.map_sampling()**2
    * ps.gridding_compensation()**2,
    ps.weights_1,
    ps.weights_1,
    )
)

plt.plot(keff, ps1d)
plt.plot(keff, ps1d_beam_map_comp_w)
plt.plot(keff, ps1d_beam_map_comp_w2,ls=':')

plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'k')
plt.ylabel(r'P(k)')

Text(0, 0.5, 'P(k)')

Note that the weights convolution is normalised in amplitude. To visualise its effect you should check mode-mixing:

ps.weights_1 = None
pscy_noweights,_ = ps.get_cy_power(ps.auto_power_tracer_1_model)
ps.weights_1 = ps.counts_in_box.astype('float')
pscy_weights,_ = ps.get_cy_power(ps.auto_power_tracer_1_model)

np.log10(pscy_noweights).min(),np.log10(pscy_noweights).max(),

(np.float64(-6.6116079471293245), np.float64(-2.52761634643415))

fig, axes = plt.subplots(1,2,figsize=(12,5))
ax = axes[0]
ax.pcolormesh(ps.kperpbins, ps.kparabins, np.log10(pscy_noweights).T, vmin=-6.7,vmax=-2.5)
ax = axes[1]
ax.pcolormesh(ps.kperpbins, ps.kparabins, np.log10(pscy_weights).T, vmin=-6.7,vmax=-2.5)

<matplotlib.collections.QuadMesh at 0x139f212a0>