Simulated Sky Signal in time domain

In this lesson we will use the TOAST Operator OpSimPySM to create timestreams for an instrument given a sky model.

In [ ]:
# Load common tools for all lessons
import sys
sys.path.insert(0, "..")
from lesson_tools import (
    check_nersc,
    fake_focalplane
)
In [ ]:
import toast
import healpy as hp
import numpy as np
In [ ]:
env = toast.Environment.get()                         
env.set_log_level("DEBUG")

Scanning strategy

Before being able to scan a map into a timestream we need to define a scanning strategy and get pointing information for each channel.

We use the same satellite scanning used in lesson 2 about scanning strategies, see the 02_Simulated_Scan_Strategies/simscan_satellite.ipynb for more details.

In [ ]:
focal_plane = fake_focalplane()
In [ ]:
focal_plane.keys()
In [ ]:
focal_plane["0A"]["fwhm_arcmin"]
In [ ]:
# Scan parameters

alpha = 50.0      # precession opening angle, degrees
beta = 45.0       # spin opening angle, degrees
p_alpha = 25.0    # precession period, minutes
p_beta = 1.25     # spin period, minutes
samplerate = 0.5  # sample rate, Hz
hwprpm = 5.0      # HWP rotation in RPM
nside = 64        # Healpix NSIDE

# We will use one observation per day, with no gaps in between, and
# run for one year.

obs_samples = int(24 * 3600.0 * samplerate) - 1
nobs = 366

# Slew the precession axis so that it completes one circle

deg_per_day = 360.0 / nobs
In [ ]:
from toast.todmap import TODSatellite, slew_precession_axis
In [ ]:
detquat = {ch: focal_plane[ch]["quat"] for ch in focal_plane}
In [ ]:
# Create distributed data

comm = toast.Comm()
data = toast.Data(comm)

# Append observations

for ob in range(nobs):
    obsname = "{:03d}".format(ob)
    obsfirst = ob * (obs_samples + 1)
    obsstart = 24 * 3600.0
    tod = TODSatellite(
        comm.comm_group, 
        detquat, 
        obs_samples, 
        firstsamp=obsfirst,
        firsttime=obsstart,
        rate=samplerate,
        spinperiod=p_beta,
        spinangle=beta,
        precperiod=p_alpha,
        precangle=alpha,
        coord="E",
        hwprpm=hwprpm
    )
    qprec = np.empty(4 * tod.local_samples[1], dtype=np.float64).reshape((-1, 4))
    slew_precession_axis(
        qprec,
        firstsamp=obsfirst,
        samplerate=samplerate,
        degday=deg_per_day,
    )
    tod.set_prec_axis(qprec=qprec)
    obs = dict()
    obs["tod"] = tod
    data.obs.append(obs)
In [ ]:
from toast.todmap import (
    get_submaps_nested,
    OpPointingHpix,
    OpAccumDiag
)
from toast.map import (
    DistPixels
)

# Make a simple pointing matrix

pointing = OpPointingHpix(nside=nside, nest=True, mode="IQU")
pointing.exec(data)

# Compute the locally hit pixels

localpix, localsm, subnpix = get_submaps_nested(data, nside)

# Construct a distributed map to store the hit map

npix = 12 * nside**2

hits = DistPixels(
    comm=data.comm.comm_world,
    size=npix,
    nnz=1,
    dtype=np.int64,
    submap=subnpix,
    local=localsm,
)
hits.data.fill(0)

# Accumulate the hit map locally

build_hits = OpAccumDiag(hits=hits)
build_hits.exec(data)

# Reduce the map across processes (a No-op in this case)

hits.allreduce()
In [ ]:
%matplotlib inline
hp.mollview(hits.data.flatten(), nest=True)

Define PySM parameters and instrument bandpasses

Then we define the sky model parameters, choosing the desired set of PySM models and then we specify the band center and the bandwidth for a top-hat bandpass. Currently top-hat bandpasses are the only type supported by the operator, in the future we will implement arbitrary bandpasses.

Then bandpass parameters can be added directly to the focal_plane dictionary:

In [ ]:
for ch in focal_plane:
    focal_plane[ch]["bandcenter_ghz"] = 70                                                                                              
    focal_plane[ch]["bandwidth_ghz"] = 10
    focal_plane[ch]["fwhm"] = focal_plane[ch]["fwhm_arcmin"] # to be fixed in TOAST
In [ ]:
pysm_sky_config = ["s1", "f1", "a1", "d1"]

Run the OpSimPySM operator

The OpSimPySM operator:

  • Creates top-hat bandpasses arrays (frequency axis and weights) as expected by PySM
  • Loops by channel and for each:
    • Creates a PySMSky object just with 1 channel at a time
    • Executes PySMSky to evaluate the sky models and bandpass-integrate
    • Calls PySM to perform distributed smoothing with libsharp
    • Gathers the map on the first MPI process
    • Applies coordinate transformation if necessary (not currently implemented in libsharp)
    • Use the DistMap object to communicate to each process the part of the sky they observe
    • Calls OpSimScan to rescan the map to a timeline
In [ ]:
from toast.todmap import OpSimPySM
OpSimPySM?
In [ ]:
opsim_pysm = OpSimPySM(
    comm=None,
    pysm_model=pysm_sky_config,
    nside=nside,
    apply_beam=True,
    debug=True,
    focalplanes=[focal_plane],
    subnpix=subnpix,
    localsm=localsm
)
In [ ]:
opsim_pysm.exec(data)

Plot output timelines

In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
In [ ]:
tod = data.obs[0]['tod']
In [ ]:
pix = tod.cache.reference("pixels_0A")
In [ ]:
import toast.qarray as qa
theta, phi, pa = qa.to_angles(tod.read_pntg(detector="0A"))
In [ ]:
pix
In [ ]:
num = 10000
plt.figure(figsize=(7, 5))
plt.plot(np.degrees(theta[:num]), tod.cache.reference("signal_0A")[:num], ".")
plt.xlabel("$Colatitude [deg]$")
plt.ylabel("$Signal [ \mu K_{RJ} ]$");

Bin the output to a map

In [ ]:
from numba import njit
In [ ]:
@njit
def just_make_me_a_map(output_map, signals):
    """Temperature only binner
    
    Bins a list of (pix, signal) tuples into an output map,
    it does not support polarization, so it just averages it out.
    
    Parameters
    ----------
    output_map : np.array
        already zeroed output map
    signals : numba.typed.List of (np.array[int64] pix, np.array[np.double] signal)
    
    Returns
    -------
    hits : np.array[np.int64]
        hitmap
    """
    hits = np.zeros(len(output_map), dtype=np.int64)
    for pix, signal in signals:
        for p,s in zip(pix, signal):
            output_map[p] += s
            hits[p] += 1
    output_map[hits != 0] /= hits[hits != 0]
    return hits
In [ ]:
from numba.typed import List

signals = List()
for obs in data.obs:
    for ch in focal_plane:
        signals.append((obs["tod"].cache.reference("pixels_%s" % ch), obs["tod"].cache.reference("signal_%s" % ch)))
In [ ]:
output_map = np.zeros(npix, dtype=np.double)
h = just_make_me_a_map(output_map, signals)
In [ ]:
hp.mollview(h, title="hitmap", nest=True)
In [ ]:
hp.mollview(output_map, nest=True, min=0, max=1e-3, cmap="coolwarm")
In [ ]:
hp.gnomview(output_map, rot=(0,0), xsize=5000, ysize=2000, cmap="coolwarm", nest=True, min=0, max=1e-2)

Custom sky components

  • pysm_component_objects: pass custom PySM component objects, see for example the WebSkyCIB model in the so_pysm_models repository, it provides a Cosmic Infrared Background computed from
In [ ]: