Generate synchrotron amplitude 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))

comp = "synch"
sub = "template"

Large scales

largescale_filename = [f"{comp}_largescale_{sub}_logpoltens_alm_lmax128_2023.02.24.fits.gz"]

assert len(largescale_filename) == 1

largescale_filename = largescale_filename[0]

largescale_filename

alm_log_pol_tens_large_scale = hp.read_alm(largescale_filename, hdu=(1,2,3))

map_log_pol_tens_large_scale = hp.alm2map(alm_log_pol_tens_large_scale.astype(np.complex128), nside=output_nside)

del alm_log_pol_tens_large_scale

Small scales modulation

modulate_alm = { k:hp.read_alm(output_dir_raw/f"synch_{k}_modulation_alms_lmax64_2023.02.24.fits.gz").astype(np.complex128) for k in ["temperature","polarization"] }

Small scales

cl_small_scale = hp.read_cl(output_dir_raw / "synch_small_scales_cl_lmax16384_2023.02.24.fits.gz")

synalm_lmax = 8192*2 # for reproducibility
# synalm_lmax = 512

np.random.seed(555)

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

alm_log_pol_tens_small_scale = [hp.almxfl(each, np.ones(output_lmax+1)) for each in alm_log_pol_tens_small_scale]
map_log_pol_tens_small_scale = hp.alm2map(alm_log_pol_tens_small_scale, nside=output_nside)
del alm_log_pol_tens_small_scale
map_log_pol_tens_small_scale[0] *= hp.alm2map(modulate_alm["temperature"], output_nside)
map_log_pol_tens_small_scale[1:] *= hp.alm2map(modulate_alm["polarization"], output_nside)
assert np.isnan(map_log_pol_tens_small_scale).sum() == 0

Combine scales

• Combine small and large scale maps

• Transform from logpoltens to IQU

• Write output map

map_log_pol_tens = map_log_pol_tens_large_scale + map_log_pol_tens_small_scale

del map_log_pol_tens_large_scale, map_log_pol_tens_small_scale

from pysm3.utils import log_pol_tens_to_map, add_metadata

output_map = log_pol_tens_to_map(map_log_pol_tens)

Galactic plane fix

galplane_fix = hp.read_map(output_dir_raw / "synch_galplane.fits.gz", (0, 1, 2, 3))

output_map *= hp.ud_grade(galplane_fix[3], output_nside)
output_map += hp.ud_grade(galplane_fix[:3] * (1 - galplane_fix[3]), output_nside)

if output_nside < 4096:
    hp.mollview((output_map - log_pol_tens_to_map(map_log_pol_tens))[0])
    hp.mollview(output_map[0])

Write outputs

from datetime import date
today = date.today()
version = today.strftime("%Y.%m.%d")
version

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

add_metadata([output_dir / f"synch_template_nside{output_nside}_{version}.fits"], coord="G", unit="uK_RJ", ref_freq="23 GHz")