Verify dust templates

import os

os.environ[
    "OMP_NUM_THREADS"
] = "64"  # 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 dust_gnilc/

#name = "template"
#name = "Td"
name = "beta"

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

template_date = "_2023.06.06" if name in ["beta", "Td"] else ""

m = hp.read_map(f"gnilc_dust_{name}_nside4096{template_date}.fits", pol)

m2048 = hp.read_map(f"gnilc_dust_{name}_nside2048{template_date}.fits", pol)

m8192 = hp.read_map(f"gnilc_dust_{name}_nside8192{template_date}.fits", pol)

cl = hp.anafast(m, lmax=int(2.5*4096)).reshape((len(pol), -1))

cl8192 = hp.anafast(m8192, lmax=int(2*8192)).reshape((len(pol), -1))

cl2048 = hp.anafast(m2048, lmax=int(2.5*2048)).reshape((len(pol), -1))

from pathlib import Path

datadir = Path("raw")

if name == "template":
    name = "logpoltens"
    suffix="_complex64"
else:
    suffix = ""

alm_large_scale = hp.read_alm(
    datadir / f"gnilc_dust_largescale_template_{name}_alm_nside2048_lmax1024{suffix}.fits.gz",
    hdu=1,
)

cl_small_scale = hp.read_cl(
    datadir / f"gnilc_dust_small_scales_{name}_cl_lmax16384{template_date}.fits.gz"
).reshape((len(pol), -1))

if name == "logpoltens":
    name = "template"

pol_label = "TEB"

for p in pol:
    plt.figure(figsize=(12,6))
    plt.loglog(cl8192[p], label=f"map at Nside=8192")
    plt.loglog(cl[p], label="map at Nside=4096")
    plt.loglog(cl2048[p], label=f"map at Nside=2048")
    if name != "template":
        plt.loglog(hp.alm2cl(alm_large_scale.astype(complex).reshape((len(pol), -1)))[p], "--", label="large scale")
        plt.loglog(cl_small_scale[p], "--", label="small scales")
    plt.legend()
    plt.title(f"{name} maps spectra comparison {pol_label[p]}")
    plt.grid();