[1]:

import os

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

Preprocess HASLAM map as it is done for Synchrotron Amplitude template

[2]:

from pathlib import Path
import healpy as hp
import matplotlib.pyplot as plt
import numpy as np

# import pymaster as nmt
from astropy.io import fits

%matplotlib inline

[3]:

plt.style.use("seaborn-talk")

[4]:

import pysm3 as pysm
import pysm3.units as u

nside = 512
lmax = 3 * nside
comp = "IQU"
components = list(enumerate(comp))
components

[4]:

[(0, 'I'), (1, 'Q'), (2, 'U')]

[5]:

datadir = Path("data")
imapfile = datadir / "haslam408_dsds_Remazeilles2014.fits"
qumapfile = datadir / "wmap_band_iqumap_r9_9yr_K_v5.fits"
if not imapfile.exists():
    !wget -O $imapfile   https://lambda.gsfc.nasa.gov/data/foregrounds/haslam_2014/haslam408_dsds_Remazeilles2014.fits
if not qumapfile.exists():
    !wget -O $qumapfile   https://lambda.gsfc.nasa.gov/data/map/dr5/skymaps/9yr/raw/wmap_band_iqumap_r9_9yr_K_v5.fits
imap = hp.read_map(imapfile)
qumap = hp.read_map(qumapfile, field=[1, 2])

[6]:

rescaling_factor = (23.0 / 0.408) ** -3.1

imap *= rescaling_factor
imap <<= u.K_RJ
imap = imap.to(u.uK_RJ)
qumap <<= u.mK_RJ
qumap = qumap.to("uK_RJ")
IQU = np.array([imap, qumap[0], qumap[1]])
FWHM_SMOOTHING = 2  # deg
IQU = hp.smoothing(IQU, fwhm=np.radians(FWHM_SMOOTHING))
IQU <<= u.uK_RJ
plt.figure(figsize=(20, 5))
for i_pol, pol in components:
    hp.mollview(
        IQU[i_pol],
        title="Haslam/WMAP K synch " + pol,
        sub=131 + i_pol,
        unit=IQU.unit,
        min=-300,
        max=300,
    )
../_images/preprocess-templates_synchrotron_curvature_6_0.png 
[7]:

from pysm3.utils import map_to_log_pol_tens

iqu = map_to_log_pol_tens(IQU.value)

[8]:

datadir = Path("data")

[9]:

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

[10]:

output_nside = 2048

[11]:

output_lmax = min(int(2.5 * output_nside), 8192*2)
output_ell = np.arange(output_lmax + 1)
output_cl_norm = output_ell * (output_ell + 1) / np.pi / 2
output_cl_norm[0] = 1

[12]:

print(
    f"{np.isnan(iqu[0]).sum()  } pixels out of { iqu[0].size} are NaNs in Log Pol Tens maps "
)
for i in range(3):
    iqu[i, np.isnan(iqu[i])] = np.nanmedian(iqu[i])
assert np.isnan(iqu).sum() == 0
print(
    f"{np.isnan(iqu[0]).sum()  } pixels out of { iqu[0].size} are NaNs in Log Pol Tens maps "
)

527 pixels out of 3145728 are NaNs in Log Pol Tens maps
0 pixels out of 3145728 are NaNs in Log Pol Tens maps

[13]:

for i_pol, pol in components:
    hp.mollview(iqu[i_pol], title="Log Pol tensor " + pol, sub=131 + i_pol)
../_images/preprocess-templates_synchrotron_curvature_13_0.png

Template for \(C_s\) is the \(i\) map smoothed at 5deg

The modulation amplitude is exactly the same as the one adopted for \(\beta_s\) and \(intensity \) small scales

[14]:

ismooth = hp.smoothing(iqu[0], fwhm=np.radians(5))
minmax = lambda m, a, b: a + (b - a) * (m - m.min()) / (m.max() - m.min())

modulate_amp = (ismooth) * 1.0

mskmd = ismooth > 5.5
b1 = 1.5
a = 0.5

mskmd = ismooth > 5

modulate_amp[mskmd] = minmax(ismooth[mskmd], 1.1, 2)
modulate_amp[~mskmd] = minmax(ismooth[~mskmd], 0.1, 1.1)

Tweak the values for MinMax rescaling so that we match Kogut 2012 measurements within ARCADE patch

[15]:

Cmap = minmax(ismooth, a=-0.041, b=-0.076)

Load ARCADE patch

make wget link

[16]:

arcadefile = datadir / "arc2_3150_v19.fits"
if not arcadefile.exists():
    !wget -O $arcadefile   https://lambda.gsfc.nasa.gov/data/suborbital/ARCADE/arc2_3150_v19.fits

arcade = hp.read_map(arcadefile)
mask = arcade != 0
mask_arcade = hp.ud_grade(mask, nside_out=hp.get_nside(Cmap))
hp.mollview(mask_arcade)
../_images/preprocess-templates_synchrotron_curvature_19_0.png 
[17]:

print(
    f" Kogut 2012 measured Cs = -0.052\pm 0.005, we get within ARCADE patch:\n"
    + f"{Cmap[mask_arcade].mean():.4f}\pm { Cmap[mask_arcade] .std():.4f}"
)

 Kogut 2012 measured Cs = -0.052\pm 0.005, we get within ARCADE patch:
-0.0517\pm 0.0054

[18]:

hp.gnomview(
    (Cmap * mask_arcade) / mask_arcade,
    title="Matched synchrotron curvature - ARCADE",
    rot=[66, 0],
    reso=10,
    xsize=780,
    cmap="jet",
    min=-0.1,
    max=0,
)

hp.mollview(Cmap, title="synchrotron curvature", cmap="jet", min=-0.1, max=0)
hp.graticule()

/tmp/ipykernel_17482/3364598603.py:2: RuntimeWarning: invalid value encountered in true_divide
  (Cmap * mask_arcade) / mask_arcade,
../_images/preprocess-templates_synchrotron_curvature_21_1.png
../_images/preprocess-templates_synchrotron_curvature_21_2.png

Inject Modulated small scales to \(C_s\)

[19]:

ell = np.arange(lmax + 1)
cl_spass_fit = 10 ** (-0.82655497) * ell ** (-2.60860208)
cl_spass_fit[0] = 0
ell_fit_low = 10
ell_fit_high = 36

/tmp/ipykernel_17482/1537673065.py:2: RuntimeWarning: divide by zero encountered in power
  cl_spass_fit = 10 ** (-0.82655497) * ell ** (-2.60860208)

[20]:

def run_anafast(m, lmax):
    clanaf = hp.anafast(m, lmax=lmax)
    cl = {}
    cl["TT"] = clanaf
    ell = np.arange(lmax + 1)

    cl_norm = ell * (ell + 1) / np.pi / 2
    cl_norm[0] = 1
    return ell, cl_norm, cl


def model(ell, A, gamma):
    out = A * ell ** gamma
    out[:2] = 0
    return out

[21]:

ell, cl_norm, cl_c = run_anafast(Cmap, lmax=lmax)

[22]:

gamma_spass = -2.60860208
A0 = cl_c["TT"][ell_fit_high] / ell_fit_high ** (gamma_spass)

if output_nside < 4096:

    plt.loglog(ell, cl_c["TT"], label="Anafast $C_\ell$")
    plt.plot(
        ell[ell_fit_high:],
        model(ell[ell_fit_high:], A0, gamma_spass),
        label="power law SPASS ",
    )

    plt.axvline(
        ell_fit_low, linestyle="--", color="black", label="$ \ell={} $".format(ell_fit_low)
    )
    plt.axvline(
        ell_fit_high, linestyle="--", color="gray", label="$ \ell={} $".format(ell_fit_high)
    )
    plt.grid()
    plt.title(r"$C_s  $" + f" power spectrum ")
    plt.legend()
    plt.ylabel("$ C_\ell [\mu K_{RJ}^2]$")
    plt.xlabel(("$\ell$"))
    plt.xlim(2, lmax)
../_images/preprocess-templates_synchrotron_curvature_26_0.png 
[23]:

from pysm3.utils import add_metadata, sigmoid

[24]:

sig_func = sigmoid(output_ell, x0=ell_fit_high, width=20)
np.random.seed(777)
alm_Cs = hp.map2alm(Cmap, lmax=lmax, use_pixel_weights=True)
Cs_ls_alm = hp.almxfl(alm_Cs, np.sqrt(1.0 - sig_func))
Cs_ss_alm = hp.synalm(model(output_ell, A0, gamma_spass), lmax=output_lmax)
Cs_ss_alm = hp.almxfl(Cs_ss_alm, np.sqrt(sig_func))
Cs_ss = hp.alm2map(Cs_ss_alm, lmax=output_lmax, nside=output_nside)
Cs_ss *= hp.ud_grade(modulate_amp, output_nside)

Cs_ls = hp.alm2map(Cs_ls_alm, nside=output_nside)
Cs_out = Cs_ss + Cs_ls
ell, _, clcs_out = run_anafast(Cs_out, lmax=lmax)
ell, _, clls = run_anafast(Cs_ls, lmax=lmax)
ell, _, clss = run_anafast(Cs_ss, lmax=lmax)

/tmp/ipykernel_17482/1065277110.py:13: RuntimeWarning: divide by zero encountered in power
  out = A * ell ** gamma

[25]:

output_map_alm = hp.map2alm(Cs_out, lmax=output_lmax)

[26]:

if output_nside > 4096:

    hp.write_alm(
        output_dir / "raw" / f"synch_curvature_alm_nside{output_nside}_lmax{output_lmax}.fits",
        output_map_alm,
        out_dtype = np.complex64,
        overwrite=True,
    )
    add_metadata([output_dir / "raw" / f"synch_curvature_alm_nside{output_nside}_lmax{output_lmax}.fits"], coord="G", unit="", ref_freq="23 GHz")

[27]:

if output_nside < 4096:

    for ii in range(1):
        plt.loglog(
            ell, cl_c["TT"], alpha=0.5, color="C%d" % ii, label=r"$ C_s$ w/o small scales"
        )
        plt.loglog(
            ell,
            clcs_out["TT"],
            label=r"$  C_s$ w/ small scales ",
        )
        # plt.loglog(ell,  clls[pol], label=r"$  \beta_s$ PySM3 ",  )
        # plt.loglog(ell,  clss[pol], label=r"$  \beta_s$ PySM3 ",  )

        plt.grid(True)

        plt.legend()

        plt.ylabel("$ C_\ell $")
        plt.xlabel(("$\ell$"))
        # plt.xlim(10, 55 )
        plt.xlim(2, lmax)
    plt.ylim(1e-15, 1e-4)
../_images/preprocess-templates_synchrotron_curvature_31_0.png 
[28]:

mask_arcade_high_nside = hp.ud_grade(mask_arcade, output_nside).astype(np.bool)

/tmp/ipykernel_17482/3137707393.py:1: DeprecationWarning: `np.bool` is a deprecated alias for the builtin `bool`. To silence this warning, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  mask_arcade_high_nside = hp.ud_grade(mask_arcade, output_nside).astype(np.bool)

[29]:

if output_nside < 4096:

    hp.gnomview(
        (Cmap * mask_arcade) / mask_arcade,
        title="Cs w/o small scales - ARCADE",
        rot=[66, 0],
        reso=10,
        sub=121,
        xsize=780,
        cmap="jet",
        min=-0.1,
        max=0,
    )
    hp.gnomview(
        (Cs_out * mask_arcade_high_nside) / mask_arcade_high_nside,
        title="Cs w/ small scales - ARCADE",
        rot=[66, 0],
        reso=10,
        sub=122,
        xsize=780,
        cmap="jet",
        min=-0.1,
        max=0,
    )
    plt.figure()
    hp.mollview(Cmap, title="Cs w/o small scales", cmap="jet", min=-0.1, max=0, sub=121)


    hp.mollview(Cs_out, title="Cs w/ small scales", cmap="jet", min=-0.1, max=0, sub=122)
    hp.graticule()

/tmp/ipykernel_17482/3357961219.py:4: RuntimeWarning: invalid value encountered in true_divide
  (Cmap * mask_arcade) / mask_arcade,
/tmp/ipykernel_17482/3357961219.py:15: RuntimeWarning: invalid value encountered in true_divide
  (Cs_out * mask_arcade_high_nside) / mask_arcade_high_nside,
../_images/preprocess-templates_synchrotron_curvature_33_1.png
../_images/preprocess-templates_synchrotron_curvature_33_2.png 
[30]:

print(
    f" Kogut 2012 measured Cs = -0.052\pm 0.005, we get within ARCADE patch:\n"
    + f"{Cs_out[mask_arcade_high_nside].mean().mean():.4f}\pm { Cs_out[mask_arcade_high_nside]  .std():.4f}"
)

 Kogut 2012 measured Cs = -0.052\pm 0.005, we get within ARCADE patch:
-0.0517\pm 0.0054

[ ]: