Verify synchrotron templates

import os

os.environ[
    "OMP_NUM_THREADS"
] = "200"  # for jupyter.nersc.gov otherwise the notebook only uses 2 cores

import pysm3
import healpy as hp
import numpy as np
import matplotlib.pyplot as plt

cd /global/cfs/cdirs/cmb/www/pysm-data

cd synch/

ls

name = "template"
#name = "curvature"
#name = "beta"

pol = (0,1,2) if name == "template" else (0,)

template_date = "_2023.02.25" if name in ["template"] else ""

nsides = [2048, 4096, 8192]

m = {}
for nside in nsides:
    m[nside] = hp.read_map(f"synch_{name}_nside{nside}{template_date}.fits", pol)

cl = {}
for nside in nsides:
    cl[nside] = hp.anafast(m[nside], lmax=int(min(2.5*nside,2*8192)))
    if cl[nside].ndim == 1:
        cl[nside] = cl[nside].reshape((1, -1))
    print(nside, cl[nside].shape)

from pathlib import Path

datadir = Path("raw")

if name == "template":
    largescale_filename = datadir / "synch_largescale_template_logpoltens_alm_lmax128_2023.02.24.fits.gz"
elif name == "curvature":
    largescale_filename = datadir / "synch_curvature_alm_nside8192_lmax16384.fits"
else:
    largescale_filename = datadir / f"synch_largescale_{name}_alm_nside512_lmax768.fits.gz"

alm_large_scale = hp.read_alm(
    largescale_filename,
    hdu=1,
)

if name == "template":
    smallscale_filename = datadir / "synch_small_scales_cl_lmax16384_2023.02.24.fits.gz"
else:
    smallscale_filename = datadir / f"synch_small_scales_{name}_cl_lmax16384.fits.gz"

if name != "curvature":
    cl_small_scale = hp.read_cl(
        smallscale_filename
    )
    cl_small_scale = cl_small_scale.reshape((len(pol) if name != "template" else 4, -1))

pol_label = "TEB"

for p in pol:
    plt.figure(figsize=(12,6))
    for nside in reversed(nsides):
        plt.loglog(cl[nside][p], label=f"map at Nside={nside}")
    if name != "template":
        plt.loglog(hp.alm2cl(alm_large_scale.astype(complex).reshape((len(pol), -1)))[p], "--", alpha=.5, label="large scale")
        if name != "curvature":
            plt.loglog(cl_small_scale[p], "--", label="small scales")
    plt.title(f"{name} maps spectra comparison {pol_label[p]}")
    plt.axvline([3150], linestyle="--", color="black", label="ell = 3150", alpha=.5)
    plt.legend()

    plt.grid();

modulation_alm = hp.read_alm(datadir / "synch_amplitude_modulation_alms_lmax768.fits.gz")

modulation_map = hp.alm2map(modulation_alm.astype(complex), nside=2048)

hp.mollview(modulation_map, title=f"Modulation map, mean: {modulation_map.mean():.2f}", unit="dimensionless")