Foreground Simulation

In meer21cm, we provide several basic foreground emission templates. They are meant for visual checks and templates for simulation, and do not incorporate effects such as polarization (and leakages). Nevertheless, by interfacing with pysm and GDSM, the templates should be of reasonable accuracy, ideal for studying foreground removal strategies and running forecasts.

The ForegroundSimulation class handles the simulation. The minimum input is simply hp_side, which is the HEALPix nside that gives the angular resolution of the foreground sky map:

[1]:

import healpy as hp
from meer21cm.fg import ForegroundSimulation
import numpy as np
from meer21cm import PowerSpectrum
from meer21cm.plot import plot_map
from meer21cm.telescope import dish_beam_sigma

[2]:

fgsim = ForegroundSimulation(
    hp_nside=128,
    backend="haslam", # default is a Haslam 408 MHz map
    sp_indx_for_haslam_backend=-2.0, # default is -2.0
)

The default is just to interpolate the Haslam map to any given frequency:

[3]:

fg_map = fgsim.fg_wcs_cube(np.array([1e9])) # 1 GHz

[4]:

hp.mollview(np.log10(fg_map[0]), title="Foreground at 1 GHz")

The default coordinate system is C (Celestial). You can also use G (Galactic):

[5]:

fgsim.coord_system = "G"
fg_map = fgsim.fg_wcs_cube(np.array([1e9])) # 1 GHz
hp.mollview(np.log10(fg_map[0]), title="Foreground at 1 GHz in Galactic coordinates")

For more complicated simulations, you can either use GDSM or PySM3:

[6]:

fgsim.backend = "gdsm"
fg_map = fgsim.fg_wcs_cube(np.array([1e9])) # 1 GHz
hp.mollview(np.log10(fg_map[0]), title="Foreground at 1 GHz from GDSM")

fgsim.backend = "pysm"
fg_map = fgsim.fg_wcs_cube(np.array([1e9])) # 1 GHz
hp.mollview(np.log10(fg_map[0]), title="Foreground at 1 GHz from PySM3")

For pysm, there is a choice of which component to include, by setting pysm_preset_strings. See pysm documentation for more information.

Say only a simple synchrotron map:

[7]:

fgsim.pysm_preset_strings = ["s1"]
fg_map = fgsim.fg_wcs_cube(np.array([1e9])) # 1 GHz
hp.mollview(np.log10(fg_map[0]), title="Foreground at 1 GHz from PySM3, only s1")

If you already have a wcs survey area specified, you can pass it into fgsim and it will automatically project the healpix map onto the wcs grids:

[8]:

ps = PowerSpectrum(
    survey='meerklass_2021',
    band='L',
)

[9]:

fgsim.coord_system = 'C'
fgsim.wproj = ps.wproj
fgsim.num_pix_x = ps.num_pix_x
fgsim.num_pix_y = ps.num_pix_y
fg_map = fgsim.fg_wcs_cube(ps.nu) # 1 GHz

The result is the foreground map cube at each frequency:

[10]:

fg_map.shape

[10]:

(133, 73, 252)

[11]:

plot_map(fg_map,ps.wproj,)

If you have a beam model you can convolve it with the beam:

[12]:

ps.sigma_beam_ch = dish_beam_sigma(13.5, ps.nu,)

[13]:

fg_w_beam,_ = ps.convolve_data(
    ps.beam_image,
    data=fg_map,
    weights=ps.w_HI,
    assign_to_self=False, # do not overwrite ps.data and ps.w_HI
)

[14]:

plot_map(fg_w_beam,ps.wproj,)