Generate GNILC spectral index and temperature maps

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
name = "beta"

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

Galactic mask

galactic_mask = (
    hp.ud_grade(
        hp.read_map(datadir / "HFI_Mask_GalPlane-apo2_2048_R2.00_GAL080_noapo.fits.gz"),
        output_nside,
    )
    == 1
)

Large scales

alm_large_scale = hp.read_alm(
    datadir / f"gnilc_dust_largescale_template_{name}_alm_nside2048_lmax1024.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(
    datadir / f"gnilc_dust_temperature_modulation_alms_lmax768.fits.gz"
).astype(np.complex128)

Small scales

cl_small_scale = hp.read_cl(
    datadir / f"gnilc_dust_small_scales_{name}_cl_lmax16384_2023.06.06.fits.gz"
)

len(cl_small_scale)

cl_small_scale

import matplotlib.pyplot as plt
plt.loglog(cl_small_scale)

synalm_lmax = 8192 * 2  # for reproducibility
# synalm_lmax = 1024
seed = 777 if name == "beta" else 888
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(int(2.5 * output_nside+1)))
map_small_scale = hp.alm2map(alm_small_scale, nside=output_nside)
assert np.isnan(map_small_scale).sum() == 0

map_small_scale = hp.alm2map(alm_small_scale, nside=output_nside)

hp.mollview(map_small_scale)

map_small_scale *= hp.alm2map(modulate_alm, output_nside)

hp.mollview(map_small_scale)

hp.mollview(galactic_mask)

Combine scales

• Combine small and large scale maps

• Transform from logpoltens to IQU

• Write output map

map_total = map_large_scale + map_small_scale

galplane_fix_mask = hp.read_map(datadir / "gnilc_dust_galplane.fits.gz", 3)

ext_number = {"beta":0, "Td":1}

galplane_fix = hp.read_map(datadir / "gnilc_dust_beta_Td_galplane.fits.gz", ext_number[name])

map_total *= hp.ud_grade(galplane_fix_mask, output_nside)
map_total += hp.ud_grade(galplane_fix * (1 - galplane_fix_mask), output_nside)

hp.mollview(map_total - map_large_scale - map_small_scale, title="Effect of galactic plane fix")

hp.write_map(
    output_dir / f"gnilc_dust_{name}_nside{output_nside}_2023.06.06.fits",
    map_total,
    coord="G",
    column_units="" if name == "beta" else "K",
    extra_header = [("lmax", output_lmax), ("ref_freq", "353 GHz")],
    dtype=np.float32,
    overwrite=True,
)

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