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

output_dir = Path("production-data") / "dust_gnilc"

datadir = output_dir / "raw"

output_nside = 2048

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

Large scales

alm_log_pol_tens_large_scale = hp.read_alm(
    datadir
    / "gnilc_dust_largescale_template_logpoltens_alm_nside2048_lmax1024_complex64.fits.gz",
    hdu=(1, 2, 3),
)

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

map_log_pol_tens_large_scale

Small scales modulation

modulate_alm = {
    k: hp.read_alm(datadir / f"gnilc_dust_{k}_modulation_alms_lmax768.fits.gz").astype(
        np.complex128
    )
    for k in ["temperature", "polarization"]
}

Small scales

cl_small_scale = hp.read_cl(
    datadir / "gnilc_dust_small_scales_logpoltens_cl_lmax16384.fits.gz"
)

synalm_lmax = 8192 * 2  # it needs to be the same for all output nside
# synalm_lmax = 1024
np.random.seed(8192)

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))
    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
)

map_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

map_log_pol_tens_small_scale

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 += map_log_pol_tens_small_scale

del map_log_pol_tens_small_scale

from pysm3.utils import log_pol_tens_to_map

output_map = log_pol_tens_to_map(map_log_pol_tens)

output_map

del map_log_pol_tens

Galactic plane fix

galplane_fix = hp.read_map(datadir / "gnilc_dust_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)

output_map

Color correction

Planck 353 GHz color correction https://github.com/galsci/pysm/issues/99

output_map *= 0.911

output_map

hp.write_map(
    output_dir / f"gnilc_dust_template_nside{output_nside}.fits",
    output_map,
    dtype=np.float32,
    overwrite=True,
    column_units = "uK_RJ",
    extra_header = [("lmax", output_lmax), ("ref_freq", "353 GHz")]
)