Quick start on mock simulations using meer21cm

In this notebook, we give a brief overview of using meer21cm to generate mock simulations of HI temperature field and galaxy sample, and show how to use the field and galaxy sample for power spectrum estimation.

meer21cm has two key features:

1. The package is easy to interact. A lot of the functionalities are wrapped under the hood.

2. A smart cache system is used, so that computationally heavy functions are only run once. However, if the relevant parameters are updated, the cache is automatically cleared so that it can be updated.

from meer21cm import MockSimulation
import numpy as np
import matplotlib.pyplot as plt
from meer21cm.plot import plot_map
import time
from meer21cm.power import get_shot_noise_galaxy
from meer21cm.grid import shot_noise_correction_from_gridding

Generation of mock signal maps:

In order to generate a mock signal and use it for power spectrum estimation, three steps must be performed:

1. Based on the survey area and redshift range of choice, generate a regular grid box with relatively fine resolution. Use the parameters to generate either a density field or a discrete source catalogue in the box.

2. Project the generated tracer onto the sky in (RA,Dec,z). For density field, it is averaged onto the map-level pixels. Noise, foregrounds, beam and other effects can be then added.

3. Use the map or the source catalogue, and grid the data into a regular grid box with relatively coarse reolustion. Perform power spectrum estimation.

Note that, each of the step has many nuances and examples to discuss these steps in details are being prepared. For now, let us just go through the whole thing without exhausting all the details.

We can construct a mock simulation by simply:

# fixing seed for reproducable results
mock = MockSimulation(
    seed=12345,
    # use pre-defined survey dimensions
    survey="meerklass_2021",
    band="L",
)

The mock object contains all sorts of configurations and settings for the simulations, and the default corresponds to the MeerKLASS L-band deep-field. For example, the frequency range and resolution are:

# in Hz
mock.nu.min(),mock.nu.max(),np.diff(mock.nu).mean()

(np.float64(971150390.625), np.float64(1023605468.75), np.float64(208984.375))

You can always update attributes of mock on the fly. For example, now I want to limit my survey area:

raminMK,ramaxMK = 334,357
decminMK,decmaxMK = -35,-26.5
ra_range = (raminMK,ramaxMK)
dec_range = (decminMK,decmaxMK)

mock.ra_range = ra_range
mock.dec_range = dec_range
mock.trim_map_to_range()

plot_map(mock.data,mock.wproj,W=mock.W_HI)

Here, W_HI is the binary window where 1 suggests that the pixels are sampled. wproj is an astropy WCS object for storing the survey area.

Now we can simply generate a mock signal map of HI overdensity:

# this is to tell mock to use **relatively fine** grid
mock.downres_factor_radial = 1/2.0
mock.downres_factor_transverse = 1/2.0
mock.data = mock.propagate_mock_field_to_data(mock.mock_tracer_field_1)

plot_map(mock.data,mock.wproj,W=mock.W_HI)

Note that, here I have already asked mock to do three different things under the hood:

• Calculate box dimensions for the density simulation

• Simulate a lognormal realization of the tracer density

• Project the tracer density field onto the map

While this is all done behine the scenes, you can also asked for it to be done separately:

# calculate box dimensions
mock.get_enclosing_box()
# generate tracer field
mock.get_mock_tracer_field(1)
# then grid
mock.data = mock.grid_field_to_sky_map(mock.mock_tracer_field_1)[0]

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[39], line 6
      4 mock.get_mock_tracer_field(1)
      5 # then grid
----> 6 mock.data = mock.grid_field_to_sky_map(mock.mock_tracer_field_1)[0]

File ~/Desktop/work/meer21cm/src/meer21cm/power.py:3274, in PowerSpectrum.grid_field_to_sky_map(self, field, average, mask, wproj, num_pix_x, num_pix_y)
   3272 pos_xyz = np.array(np.meshgrid(*self.x_vec, indexing="ij")).reshape((3, -1)).T
   3273 pos_ra, pos_dec, pos_z, _ = self.ra_dec_z_for_coord_in_box(pos_xyz)
-> 3274 pos_indx_1, pos_indx_2 = radec_to_indx(pos_ra, pos_dec, wproj, to_int=False)
   3275 pos_indx_z = redshift_to_freq(pos_z) - self.nu.min()
   3276 pos_indx_array = np.array([pos_indx_1, pos_indx_2, pos_indx_z]).T

File ~/Desktop/work/meer21cm/src/meer21cm/util.py:808, in radec_to_indx(ra_arr, dec_arr, wproj, to_int, ang_unit)
    806 def radec_to_indx(ra_arr, dec_arr, wproj, to_int=True, ang_unit="deg"):
    807     coor = SkyCoord(ra_arr, dec_arr, unit=ang_unit)
--> 808     indx_1, indx_2 = wproj.world_to_pixel(coor)
    809     if to_int:
    810         indx_1 = np.round(indx_1).astype("int")

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcsapi/high_level_api.py:321, in HighLevelWCSMixin.world_to_pixel(self, *world_objects)
    316 world_values = high_level_objects_to_values(
    317     *world_objects, low_level_wcs=self.low_level_wcs
    318 )
    320 # Finally we convert to pixel coordinates
--> 321 pixel_values = self.low_level_wcs.world_to_pixel_values(*world_values)
    323 return pixel_values

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcsapi/fitswcs.py:344, in FITSWCSAPIMixin.world_to_pixel_values(self, *world_arrays)
    341 from astropy.wcs.wcs import NoConvergence
    343 try:
--> 344     pixel = self.all_world2pix(*world_arrays, 0)
    345 except NoConvergence as e:
    346     warnings.warn(str(e))

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcs.py:2110, in WCS.all_world2pix(self, tolerance, maxiter, adaptive, detect_divergence, quiet, *args, **kwargs)
   2107 if self.wcs is None:
   2108     raise ValueError("No basic WCS settings were created.")
-> 2110 return self._array_converter(
   2111     lambda *args, **kwargs: self._all_world2pix(
   2112         *args,
   2113         tolerance=tolerance,
   2114         maxiter=maxiter,
   2115         adaptive=adaptive,
   2116         detect_divergence=detect_divergence,
   2117         quiet=quiet,
   2118     ),
   2119     "input",
   2120     *args,
   2121     **kwargs,
   2122 )

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcs.py:1546, in WCS._array_converter(self, func, sky, ra_dec_order, *args)
   1540     except Exception:
   1541         raise TypeError(
   1542             "When providing more than two arguments, they must be "
   1543             "a 1-D array for each axis, followed by an origin."
   1544         )
-> 1546     return _return_list_of_arrays(axes, origin)
   1548 raise TypeError(
   1549     f"WCS projection has {self.naxis} dimensions, so expected 2 (an Nx{self.naxis} array "
   1550     f"and the origin argument) or {self.naxis + 1} arguments (the position in each "
   1551     f"dimension, and the origin argument). Instead, {len(args)} arguments were "
   1552     "given."
   1553 )

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcs.py:1496, in WCS._array_converter.<locals>._return_list_of_arrays(axes, origin)
   1494 if ra_dec_order and sky == "input":
   1495     xy = self._denormalize_sky(xy)
-> 1496 output = func(xy, origin)
   1497 if ra_dec_order and sky == "output":
   1498     output = self._normalize_sky(output)

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcs.py:2111, in WCS.all_world2pix.<locals>.<lambda>(*args, **kwargs)
   2107 if self.wcs is None:
   2108     raise ValueError("No basic WCS settings were created.")
   2110 return self._array_converter(
-> 2111     lambda *args, **kwargs: self._all_world2pix(
   2112         *args,
   2113         tolerance=tolerance,
   2114         maxiter=maxiter,
   2115         adaptive=adaptive,
   2116         detect_divergence=detect_divergence,
   2117         quiet=quiet,
   2118     ),
   2119     "input",
   2120     *args,
   2121     **kwargs,
   2122 )

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcs.py:1902, in WCS._all_world2pix(self, world, origin, tolerance, maxiter, adaptive, detect_divergence, quiet)
   1692 def _all_world2pix(
   1693     self, world, origin, tolerance, maxiter, adaptive, detect_divergence, quiet
   1694 ):
   (...)
   1900
   1901     # initial approximation (linear WCS based only)
-> 1902     pix0 = self.wcs_world2pix(world, origin)
   1904     # Check that an iterative solution is required at all
   1905     # (when any of the non-CD-matrix-based corrections are
   1906     # present). If not required return the initial
   1907     # approximation (pix0).
   1908     if not self.has_distortion:
   1909         # No non-WCS corrections detected so
   1910         # simply return initial approximation:

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcs.py:2451, in WCS.wcs_world2pix(self, *args, **kwargs)
   2449 if self.wcs is None:
   2450     raise ValueError("No basic WCS settings were created.")
-> 2451 return self._array_converter(
   2452     lambda xy, o: self.wcs.s2p(xy, o)["pixcrd"], "input", *args, **kwargs
   2453 )

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcs.py:1532, in WCS._array_converter(self, func, sky, ra_dec_order, *args)
   1530     if xy.shape == () or len(xy.shape) == 1:
   1531         return _return_list_of_arrays([xy], origin)
-> 1532     return _return_single_array(xy, origin)
   1534 elif len(args) == self.naxis + 1:
   1535     axes = args[:-1]

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcs.py:1515, in WCS._array_converter.<locals>._return_single_array(xy, origin)
   1513 if ra_dec_order and sky == "input":
   1514     xy = self._denormalize_sky(xy)
-> 1515 result = func(xy, origin)
   1516 if ra_dec_order and sky == "output":
   1517     result = self._normalize_sky(result)

File ~/miniconda3/envs/meer21cm/lib/python3.10/site-packages/astropy/wcs/wcs.py:2452, in WCS.wcs_world2pix.<locals>.<lambda>(xy, o)
   2449 if self.wcs is None:
   2450     raise ValueError("No basic WCS settings were created.")
   2451 return self._array_converter(
-> 2452     lambda xy, o: self.wcs.s2p(xy, o)["pixcrd"], "input", *args, **kwargs
   2453 )

KeyboardInterrupt:

Note that here, I choose the first tracer mock.mock_tracer_field_1 to grid onto the sky. When it is invoked, mock.mock_tracer_field_1 automatically gets calculated under the hood. When you want to use it again, it is already cached:

s_t = time.time()
mock.mock_tracer_field_1
print(time.time()-s_t)

0.0075490474700927734

But if you modify related quantities, for example the underlying cosmology, it will be updated automatically:

print('Cosmology was:', mock.cosmo)
mock.cosmo = 'WMAP1'
print('Cosmology now:', mock.cosmo)
s_t = time.time()
mock.mock_tracer_field_1
print(time.time()-s_t)

Cosmology was: w0waCDM(name="new", H0=67.66 km / (Mpc s), Om0=0.30966, Ode0=0.6888523996875395, Tcmb0=2.7255 K, Neff=3.046, m_nu=[0.   0.   0.06] eV, Ob0=0.04897, w0=-1.0, wa=0.0)
Cosmology now: w0waCDM(name="new", H0=72.0 km / (Mpc s), Om0=0.257, Ode0=0.7429192958579401, Tcmb0=2.7255 K, Neff=3.046, m_nu=[0. 0. 0.] eV, Ob0=0.0436, w0=-1.0, wa=0.0)
1.202225685119629

Note that here the original ‘new’ cosmology corresponds to Planck18. Note that, meer21cm assumes w0waCDM, and the input parameter set for the cosmology is:

from meer21cm.cosmology import fiducial_dict
print(fiducial_dict)

{'tau': 0.0561, 'Neff': 3.046, 'omega_cold': 0.30966, 'As': np.float64(2.105209331337507e-09), 'omega_baryon': 0.04897, 'ns': 0.9665, 'h': np.float64(0.6766), 'neutrino_mass': np.float64(0.06), 'w0': -1.0, 'wa': 0.0}

You can also manually set any of the above parameters and pass it back:

cosmology = fiducial_dict.copy()
cosmology['h'] = 0.7
cosmology['omega_cold'] = 0.25
cosmology['omega_baryon'] = 0.02
cosmology['neutrino_mass'] = 0.01
cosmology['w0']=-0.95
cosmology['wa']=-0.05
mock.cosmo = cosmology

s_t = time.time()
mock.mock_tracer_field_1
print(time.time()-s_t)

1.685748815536499

mock.data = mock.propagate_mock_field_to_data(mock.mock_tracer_field_1)

plot_map(mock.data,mock.wproj,W=mock.W_HI)

Note the small difference in colour bar, you can see indeed it has been updated. If you change irrelevant quantity, for example the bias for tracer 2, it does not clear the cache:

# change tracer 2
mock.tracer_bias_2 = 1.0
# still use tracer 1
mock.data = mock.propagate_mock_field_to_data(mock.mock_tracer_field_1)

plot_map(mock.data,mock.wproj,W=mock.W_HI)

The overdensity by default does not have unit. However, you can specify a number, or even a function as the mean amplitude mock.mean_amp_1. Currently, HI temperature is supported:

mock.omegahi = 5.4e-4
# in Kelvin
print(mock.average_hi_temp)
# you can simply tell mock to use this as the amplitude
mock.mean_amp_1 = 'average_hi_temp'

0.00010919695114461744

Again, the stored overdensity will be multiplied by the amplitude automatically. Just regrid it onto the sky:

mock.data = mock.propagate_mock_field_to_data(mock.mock_tracer_field_1)

plot_map(mock.data,mock.wproj,W=mock.W_HI)

Now your map is in Kelvin unit.

Generate of mock galaxy catalogues

We may also want to generate some mock galaxies for cross-correlation. It is again done almost automatically:

# base the galaxy catalogue on tracer 2
mock.discrete_base_field = 2
# specify galaxy bias
mock.tracer_bias_2 = 1.0
# specify how many galaxies you want
mock.num_discrete_source = 100000
mock.propagate_mock_tracer_to_gal_cat()

Once again, many things are done under the hood.

First, a discrete sample of tracers are generated in the box.

mock.get_mock_tracer_position_in_box(mock.discrete_base_field)

Then, it is rotated back to (RA,Dec,z) coordinates.

Note that the number of source positions generated will be larger than mock.num_discrete_source. This is because the box will be larger than the survey lightcone to ensure the lightcone will be fully covered.

plot_map(mock.W_HI,mock.wproj,W=mock.W_HI)
ax = plt.gca()
plt.scatter(
    mock.ra_mock_tracer,
    mock.dec_mock_tracer,
    transform=ax.get_transform('world'),s=1,label='mock tracers',alpha=0.2
)

<matplotlib.collections.PathCollection at 0x1217fee00>

And finally, excess sources will be excluded and the rest stored as the galaxy catalogue:

mock.propagate_mock_tracer_to_gal_cat()

plot_map(mock.data,mock.wproj,W=mock.W_HI)
ax = plt.gca()
plt.scatter(
    mock.ra_gal,
    mock.dec_gal,
    transform=ax.get_transform('world'),s=1,label='mock gals'
)

<matplotlib.collections.PathCollection at 0x122babe20>

You can then use the signal map and the galaxy catalogue for power spectrum estimation, and compare it with model power spectrum.

# internal parameter for compensating gridding effects
# will have more description in another tutorial
mock.compensate= True
# this is to tell mock to use **relatively coarse** resolution
mock.downres_factor_radial = 4.0
mock.downres_factor_transverse = 1.2
mock.grid_scheme = 'cic'
# this is used for shot noise calculation later
himap_rg,_,_ = mock.grid_data_to_field()
# project the galaxy number density to a regular grid
galcount_rg,_,_ = mock.grid_gal_to_field()

By default, sky map is stored as field 1 and galaxy catalogue is stored as field 2. You can retrieve both the field-level power spectrum and the model power very easily:

# 1d k-bins in Mpc-1
k1dedges = np.linspace(0.05, 0.5, 10)
# pass it to mock
mock.k1dbins = k1dedges
# get data
mock.field_1 = himap_rg
mock.field_2 = galcount_rg
mock.weights_grid_1 =  mock.counts_in_box.astype('float')
mock.weights_grid_2 = None
mock.weights_field_1 = None
mock.weights_field_2 =  ((mock.counts_in_box) > 0.0).astype('float')
pdata_1d_gg,keff,nmodes = mock.get_1d_power(
    'auto_power_3d_2',
)
pdata_1d_hi,keff_hi,nmodes_hi = mock.get_1d_power(
    'auto_power_3d_1',
)
pdata_1d_c,keff_c,nmodes_c = mock.get_1d_power(
    'cross_power_3d',
)
shot_noise_g = get_shot_noise_galaxy(galcount_rg,mock.box_len, mock.weights_field_2)
shot_noise_3d = shot_noise_correction_from_gridding(
    mock.box_ndim,
    mock.grid_scheme,
) * shot_noise_g
psn_1d,_,_ = mock.get_1d_power(shot_noise_3d)
# get model
pmod_1d_gg,_,_ = mock.get_1d_power(
    'auto_power_tracer_2_model'
)
pmod_1d_hi,_,_ = mock.get_1d_power(
    'auto_power_tracer_1_model'
)
pmod_1d_cross,_,_ = mock.get_1d_power(
    'cross_power_tracer_model'
)

# a naive error bar
perror_1d_gg = (pdata_1d_gg)/np.sqrt(nmodes)
perror_1d_hi = (pdata_1d_hi)/np.sqrt(nmodes_hi)
perror_1d_cross = np.sqrt(pdata_1d_c**2 + pdata_1d_gg*pdata_1d_hi)/np.sqrt(nmodes_c)

plt.rcParams.update({'font.size':15})
plt.errorbar(keff,pdata_1d_gg-psn_1d,yerr=perror_1d_gg)

plt.plot(keff,pmod_1d_gg)

plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'$\rm P_g(k)\,[Mpc^3]$')
plt.xlabel(r'$\rm k\,[Mpc^{-1}]$')

Text(0.5, 0, '$\\rm k\\,[Mpc^{-1}]$')

plt.rcParams.update({'font.size':15})
plt.errorbar(keff_hi,pdata_1d_hi,yerr=perror_1d_hi)

plt.plot(keff_hi,pmod_1d_hi)

plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'$\rm P_{HI}(k)\,[K^2Mpc^3]$')
plt.xlabel(r'$\rm k\,[Mpc^{-1}]$')

Text(0.5, 0, '$\\rm k\\,[Mpc^{-1}]$')

plt.rcParams.update({'font.size':15})
plt.errorbar(keff_c,pdata_1d_c,yerr=perror_1d_cross)

plt.plot(keff,pmod_1d_cross)

plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'$\rm P_{cross}(k)\,[KMpc^3]$')
plt.xlabel(r'$\rm k\,[Mpc^{-1}]$')

Text(0.5, 0, '$\\rm k\\,[Mpc^{-1}]$')

Note that, the model power spectrum is automatically set to include all the effects from potentially beam convolution, weights convolution, sampling effects etc. The weighting and impacts on the power spectrum can be found in the dedicated cookbook

Internally, mock uses a fixed seed number to generate the random fields. If not specified, it will use a random integer from 0 to 2^32:

mock_temp = MockSimulation()
mock_temp.seed

2268282742

Therefore, you can easily run independent realizations of the signal map for various purposes. Let us run a bunch of signal-only realizations to see if our power spectrum estimation matches the input, and also to calculate the signal covariance:

def run_one_realization(seed=None):
    mock = MockSimulation(
        survey="meerklass_2021",band="L",
        ra_range=ra_range,dec_range=dec_range,mean_amp_1='average_hi_temp',
        tracer_bias_2=1.0,tracer_bias_1=1.0,
        downres_factor_radial=0.5,downres_factor_transverse=0.5,
        k1dbins = k1dedges,num_discrete_source=100000,compensate=True,
        seed=seed,kmax=5.0,grid_scheme='cic',
    )
    mock.data = mock.propagate_mock_field_to_data(mock.mock_tracer_field_1)
    mock.propagate_mock_tracer_to_gal_cat()
    mock.downres_factor_radial=4.0
    mock.downres_factor_transverse=1.2
    himap_rg,_,_ = mock.grid_data_to_field()
    galcount_rg,_,_ = mock.grid_gal_to_field()
    mock.sampling_resol = np.array([
        mock.pix_resol_in_mpc,
        mock.pix_resol_in_mpc,
        mock.los_resol_in_mpc,
    ])
    mock.compensate = True
    mock.include_sky_sampling = [True, False]
    mock.field_1 = himap_rg
    mock.field_2 = galcount_rg
    mock.weights_grid_1 = (mock.counts_in_box).astype('float')
    mock.weights_grid_2 = None
    mock.weights_field_1 = None
    mock.weights_field_2 = ((mock.counts_in_box) > 0.0).astype('float')
    shot_noise_g = get_shot_noise_galaxy(galcount_rg,mock.box_len, mock.weights_field_2)
    shot_noise_3d = shot_noise_correction_from_gridding(
        mock.box_ndim,
        mock.grid_scheme,
    ) * shot_noise_g
    # get model
    pmod_1d_gg,_,_ = mock.get_1d_power(
        'auto_power_tracer_2_model'
    )
    pmod_1d_hi,_,_ = mock.get_1d_power(
        'auto_power_tracer_1_model'
    )
    pmod_1d_cross,_,_ = mock.get_1d_power(
        'cross_power_tracer_model'
    )
    # get data
    pdata_1d_gg,keff,nmodes = mock.get_1d_power(
        'auto_power_3d_2',
    )
    pdata_1d_hi,keff_hi,nmodes_hi = mock.get_1d_power(
        'auto_power_3d_1',
    )
    pdata_1d_c,keff_c,nmodes_c = mock.get_1d_power(
        'cross_power_3d',
    )
    psn_1d,_,_ = mock.get_1d_power(shot_noise_3d)
    return [[pmod_1d_gg,pmod_1d_hi,pmod_1d_cross,pdata_1d_gg-psn_1d,pdata_1d_hi,pdata_1d_c],]

power_arr = []
s_t = time.time()
for i in range(10):
    power_arr += run_one_realization(i)
t_used = time.time()-s_t

print('Time used in seconds: ',time.time()-s_t)

Time used in seconds:  92.46432900428772

power_arr = np.array(power_arr)

i = 0
plt.rcParams.update({'font.size':15})
plt.errorbar(keff,power_arr.mean(0)[i+3],yerr=power_arr.std(0)[i+3])

plt.plot(keff,power_arr.mean(0)[i])

plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'$\rm P_g(k)\,[Mpc^3]$')
plt.xlabel(r'$\rm k\,[Mpc^{-1}]$')

Text(0.5, 0, '$\\rm k\\,[Mpc^{-1}]$')

i = 1
plt.rcParams.update({'font.size':15})
plt.errorbar(keff,power_arr.mean(0)[i+3],yerr=power_arr.std(0)[i+3])

plt.plot(keff,power_arr.mean(0)[i])

plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'$\rm P_{HI}(k)\,[K^2Mpc^3]$')
plt.xlabel(r'$\rm k\,[Mpc^{-1}]$')

Text(0.5, 0, '$\\rm k\\,[Mpc^{-1}]$')

i = 2
plt.rcParams.update({'font.size':15})
plt.errorbar(keff,power_arr.mean(0)[i+3],yerr=power_arr.std(0)[i+3])

plt.plot(keff,power_arr.mean(0)[i])

plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'$\rm P_{cross}(k)\,[KMpc^3]$')
plt.xlabel(r'$\rm k\,[Mpc^{-1}]$')

Text(0.5, 0, '$\\rm k\\,[Mpc^{-1}]$')

You can also visualise the correlation matrices of the signals (illustration only, covariance needs much more than just 10 realizations to converge)

power_cov = power_arr-power_arr.mean(0)[None]
power_cov = (power_cov[:,:,:,None]*power_cov[:,:,None,:]).mean(0)

fig,axes= plt.subplots(1,3)
titlearr = ['galaxy','HI','cross']
for i in range(3):
    cov_arr = power_cov[i+3]
    cov_diag = np.sqrt((np.diagonal(cov_arr)))
    cov_arr /= cov_diag[:,None] * cov_diag[None,:]
    axes[i].imshow(cov_arr,vmin=0,vmax=1,cmap='bwr')
    axes[i].set_title(titlearr[i])

Note that, the pixel resolution in the default setting is:

mock.pix_resol_in_mpc, mock.pix_resol # in Mpc and in deg

(np.float64(8.832036679403524), np.float64(0.3))

Which corresponds to the scale of:

np.pi/mock.pix_resol_in_mpc

np.float64(0.3557042127005707)

Scales close to the map resolution will inevitably be harder to measure. Note that in reality, MeerKAT has a ~1 deg beam. We can add in the beam effects quite easily (see separate tutorials on beam (TODO) for more details):

from meer21cm.telescope import dish_beam_sigma

D_dish = 13.5 #in m
mock.sigma_beam_ch = dish_beam_sigma(D_dish,mock.nu)

The default beam is a simple Gaussian:

plot_map(mock.beam_image,mock.wproj)

You can convolve this beam with the mock signal map:

mock.convolve_data(mock.beam_image)

plot_map(mock.data,mock.wproj,W=mock.W_HI)

Then redo the gridding and power spectrum estimation. The model power spectrum will also be automatically updated to include the beam attenuation effects. Let’s redo the realizations to see the differences:

def run_one_realization_beam(seed=None):
    mock = MockSimulation(
        survey="meerklass_2021",band="L",
        ra_range=ra_range,dec_range=dec_range,mean_amp_1='average_hi_temp',
        tracer_bias_2=1.0,tracer_bias_1=1.0,
        downres_factor_radial=0.5,downres_factor_transverse=0.5,
        k1dbins = k1dedges,num_discrete_source=100000,compensate=True,
        seed=seed,kmax=5.0,grid_scheme='cic',
    )
    mock.sigma_beam_ch = dish_beam_sigma(D_dish,mock.nu)
    mock.data = mock.propagate_mock_field_to_data(mock.mock_tracer_field_1)
    mock.propagate_mock_tracer_to_gal_cat()
    mock.downres_factor_radial=4.0
    mock.downres_factor_transverse=1.2
    himap_rg,_,_ = mock.grid_data_to_field()
    galcount_rg,_,_ = mock.grid_gal_to_field()
    mock.sampling_resol = np.array([
        mock.pix_resol_in_mpc,
        mock.pix_resol_in_mpc,
        mock.los_resol_in_mpc,
    ])
    mock.compensate = True
    mock.include_sky_sampling = [True, False]
    mock.field_1 = himap_rg
    mock.field_2 = galcount_rg
    mock.weights_grid_1 = (mock.counts_in_box).astype('float')
    mock.weights_grid_2 = None
    mock.weights_field_1 = None
    mock.weights_field_2 = ((mock.counts_in_box) > 0.0).astype('float')
    shot_noise_g = get_shot_noise_galaxy(galcount_rg,mock.box_len, mock.weights_field_2)
    shot_noise_3d = shot_noise_correction_from_gridding(
        mock.box_ndim,
        mock.grid_scheme,
    ) * shot_noise_g
    # get model
    pmod_1d_gg,_,_ = mock.get_1d_power(
        'auto_power_tracer_2_model'
    )
    pmod_1d_hi,_,_ = mock.get_1d_power(
        'auto_power_tracer_1_model'
    )
    pmod_1d_cross,_,_ = mock.get_1d_power(
        'cross_power_tracer_model'
    )
    # get data
    pdata_1d_gg,keff,nmodes = mock.get_1d_power(
        'auto_power_3d_2',
    )
    pdata_1d_hi,keff_hi,nmodes_hi = mock.get_1d_power(
        'auto_power_3d_1',
    )
    pdata_1d_c,keff_c,nmodes_c = mock.get_1d_power(
        'cross_power_3d',
    )
    psn_1d,_,_ = mock.get_1d_power(shot_noise_3d)
    return [[pmod_1d_gg,pmod_1d_hi,pmod_1d_cross,pdata_1d_gg-psn_1d,pdata_1d_hi,pdata_1d_c],]

%%capture
power_arr_beam = []
s_t = time.time()
for i in range(10,20):
    power_arr_beam += run_one_realization_beam(i)
t_used = time.time()-s_t

print('Time used in seconds: ',time.time()-s_t)

Time used in seconds:  95.62058305740356

power_arr_beam = np.array(power_arr_beam)

Galaxy is not affected by beam:

i = 0
plt.rcParams.update({'font.size':15})
plt.errorbar(keff,power_arr_beam.mean(0)[i+3],yerr=power_arr_beam.std(0)[i+3])

plt.plot(keff,power_arr_beam.mean(0)[i])
plt.plot(keff,power_arr.mean(0)[i],color='black',ls='--',alpha=0.5)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'$\rm P_g(k)\,[Mpc^3]$')
plt.xlabel(r'$\rm k\,[Mpc^{-1}]$')

Text(0.5, 0, '$\\rm k\\,[Mpc^{-1}]$')

HI is:

i = 1
plt.rcParams.update({'font.size':15})
plt.errorbar(keff,power_arr_beam.mean(0)[i+3],yerr=power_arr_beam.std(0)[i+3],label='With Beam')
plt.errorbar(keff,power_arr.mean(0)[i+3],yerr=power_arr.std(0)[i+3],label='Without Beam',
             color='black',ls='--',alpha=0.5)

plt.plot(keff,power_arr_beam.mean(0)[i])
plt.plot(keff,power_arr.mean(0)[i],color='black',ls='--',alpha=0.5)

plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'$\rm P_{HI}(k)\,[K^2Mpc^3]$')
plt.xlabel(r'$\rm k\,[Mpc^{-1}]$')
plt.legend()

<matplotlib.legend.Legend at 0x121e7ab00>

i = 2
plt.rcParams.update({'font.size':15})
plt.errorbar(keff,power_arr_beam.mean(0)[i+3],yerr=power_arr_beam.std(0)[i+3],label='With Beam')
plt.errorbar(keff,power_arr.mean(0)[i+3],yerr=power_arr.std(0)[i+3],label='Without Beam',
             color='black',ls='--',alpha=0.5)

plt.plot(keff,power_arr_beam.mean(0)[i])
plt.plot(keff,power_arr.mean(0)[i],color='black',ls='--',alpha=0.5)

plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'$\rm P_{cross}(k)\,[KMpc^3]$')
plt.xlabel(r'$\rm k\,[Mpc^{-1}]$')
plt.legend()

<matplotlib.legend.Legend at 0x121dad1e0>