Generate synchrotron spectral index map

import numpy as np
import healpy as hp
from pathlib import Path

datadir=Path("data/")
output_dir = Path("production-data/synch")
output_dir_raw = output_dir / "raw"

output_nside = 2048

output_lmax = int(min(2.5 * output_nside, 8192 * 2))

Large scales

alm_large_scale = hp.read_alm(
    output_dir_raw / f"synch_largescale_beta_alm_nside512_lmax768.fits.gz",
    hdu=1,
)

map_large_scale = hp.alm2map(alm_large_scale.astype(np.complex128), nside=output_nside)

Small scales modulation

modulate_alm = hp.read_alm(
    output_dir_raw / f"synch_amplitude_modulation_alms_lmax768.fits.gz"
).astype(np.complex128)

Small scales

cl_small_scale = hp.read_cl(
    output_dir_raw / f"synch_small_scales_beta_cl_lmax16384.fits.gz"
)

synalm_lmax = 8192 * 2  # for reproducibility
# synalm_lmax = output_lmax
seed = 444
np.random.seed(seed)

alm_small_scale = hp.synalm(
    cl_small_scale,
    lmax=synalm_lmax,
    new=True,
)

alm_small_scale = hp.almxfl(alm_small_scale, np.ones(output_lmax+1))
map_small_scale = hp.alm2map(alm_small_scale, nside=output_nside)
map_small_scale *= hp.alm2map(modulate_alm, output_nside)
assert np.isnan(map_small_scale).sum() == 0

Combine scales

• Combine small and large scale maps

• Write output map

output_map = map_large_scale + map_small_scale - 3.1

hp.write_map(output_dir / f"synch_beta_nside{output_nside}.fits", output_map, dtype=np.float32, overwrite=True)

from pysm3.utils import add_metadata

add_metadata([output_dir / f"synch_beta_nside{output_nside}.fits"], coord="G", unit="", ref_freq="23 GHz", ell_max=str(output_lmax))

hp.mollview(map_large_scale, title="Large scale")
hp.mollview(map_small_scale, title="Small scale")
hp.mollview(output_map, title="Total")