Scanning a map to fill a TOD

The framework provides scan_map(), a routine which scans an input map accordingly to the scanning strategy and fills the detector timestreams. It implements three possible algebras:

• No HWP:

(1)\[d_t = T + \gamma \left(Q\cos(2\theta + 2\psi_t)+U\sin(2\theta + 2\psi_t)\right),\]

where \(\theta\) is the polarization angle of the detecotor, \(\psi_t\) is the orientation of the telescope at the time \(t\), and \(\gamma\) is the polarization efficiency.

• Ideal HWP:

(2)\[d_t = T + \gamma \left(Q\cos(4\alpha_t - 2\theta + 2\psi_t)+U\sin(4\alpha_t - 2\theta + 2\psi_t)\right),\]

where \(\alpha_t\) is the HWP angle at the time \(t\).

• Generic HWP:

(3)\[d_t = (1,0,0,0) \times M_{\rm pol} \times R(\theta) \times R^T(\alpha_t) \times M_{\rm HWP} \times R(\alpha_t) \times R(\psi_t) \times \vec{S}\]

where

• \(M_{\rm pol}\) is mueller matrix of the polarimeter;

• \(M_{\rm HWP}\) is mueller matrix of the HWP;

• \(R\) is rotation matrix;

• \(\vec{S}\) is the Stokes vector \((T, Q, U, V)^T\).

Supported Sky Inputs

The routine supports different types of sky descriptions:

1. HealpixMap: Real-space maps sampled using nearest-neighbor lookup via ang2pix.

2. SphericalHarmonics: Harmonic space input, interpolated on-the-fly using the ducc0 backend.

3. Dictionary: A mapping where keys are detector or channel names and values are any of the above objects.

Implementation Example

You can fill with signal an existing TOD by using the function scan_map_in_observations(), as the following example shows:

import litebird_sim as lbs
import numpy as np

start_time_s = 0
time_span_s = 1

nside = 256
npix = 12 * nside * nside

# Create a simulation
sim = lbs.Simulation(
    base_path="./output",
    start_time=start_time_s,
    duration_s=time_span_s,
    random_seed=12345,
)

# Define the scanning strategy
sim.set_scanning_strategy(
    lbs.SpinningScanningStrategy(
        spin_sun_angle_rad=0.785_398_163_397_448_3,
        precession_rate_hz=8.664_850_513_998_931e-05,
        spin_rate_hz=0.000_833_333_333_333_333_4,
        start_time=start_time_s,
    ),
    delta_time_s=7200,
)

sim.set_instrument(
    lbs.InstrumentInfo(
        boresight_rotangle_rad=0.0,
        spin_boresight_angle_rad=0.872_664_625_997_164_8,
        spin_rotangle_rad=3.141_592_653_589_793,
    ),
)

# Create a detector object
det = lbs.DetectorInfo(
    name="Detector",
    sampling_rate_hz=10,
    quat=[0.0, 0.0, 0.0, 1.0],
)

# Initialize the observation
(obs,) = sim.create_observations(detectors=[det])

# Prepare the quaternions used to compute the pointings
sim.prepare_pointings()

# Create a map to scan (in realistic simulations,
# use the input_sky module provided by litebird_sim)
values_tuple = (np.ones(npix), np.ones(npix), np.ones(npix))
maps = lbs.HealpixMap(values=values_tuple, nside=nside, coordinates=lbs.CoordinateSystem.Ecliptic)

# Here scan the map and fill tod
lbs.scan_map_in_observations(
    obs,
    maps,
)

for i in range(obs.n_samples):
    # + 0. removes leading minus from negative zero
    value = np.round(obs.tod[0][i], 5) + 0.
    print(f"{value:.5f}")

0.00000
-0.00075
-0.00151
-0.00226
-0.00301
-0.00376
-0.00451
-0.00526
-0.00601
-0.00676

The code automatically selects the fastest algebra based on the provided HWP.

The pointing information can be included in the observation or passed through pointings. If both observations and pointings are provided, they must be coherent, so either a single Observation and a single numpy array, or same lenght list of Observations and numpy arrays. The effect of a possible HWP is included in the pointing information, see Scanning strategy, the polarization angle of the detectors is taken from the corresponding attributes included in the observations. The same applies to the polarization efficiency.

The function automatically detects the coordinate system of the input map (e.g., Galactic, Ecliptic) and rotates the pointings accordingly.

The low-level function, scan_map(), allows for direct manipulation of TOD arrays and pointing matrices for advanced users.

Methods of the Simulation class

The class Simulation provides the function Simulation.fill_tods(), which takes a map and scans it. Using this with Simulation.add_noise(), the generation of a simulation becomes quite transparent:

import litebird_sim as lbs
from astropy.time import Time
import numpy as np

start_time = 0
time_span_s = 1000.0
sampling_hz = 10.0
nside = 128

sim = lbs.Simulation(
    start_time=start_time,
    duration_s=time_span_s,
    random_seed=12345,
)

# We pick a simple scanning strategy where the spin axis is aligned
# with the Sun-Earth axis, and the spacecraft spins once every minute
sim.set_scanning_strategy(
    lbs.SpinningScanningStrategy(
        spin_sun_angle_rad=np.deg2rad(0),
        precession_rate_hz=0,
        spin_rate_hz=1 / 60,
        start_time=start_time,
    ),
    delta_time_s=5.0,
 )

# We simulate an instrument whose boresight is perpendicular to
# the spin axis.
sim.set_instrument(
    lbs.InstrumentInfo(
        boresight_rotangle_rad=0.0,
        spin_boresight_angle_rad=np.deg2rad(90),
        spin_rotangle_rad=np.deg2rad(75),
    )
)

# A simple detector looking along the boresight direction
det = lbs.DetectorInfo(
    name="Boresight_detector",
    sampling_rate_hz=sampling_hz,
    bandcenter_ghz=100.0,
    net_ukrts=50.0,
)

sim.create_observations(detectors=det)

sim.prepare_pointings()

values_tuple = (np.ones(12*nside*nside)*1e-4, np.ones(12*nside*nside)*1e-4, np.ones(12*nside*nside)*1e-4)
sky_signal = lbs.HealpixMap(values=values_tuple, nside=nside, coordinates=lbs.CoordinateSystem.Galactic)

sim.fill_tods(sky_signal)

sim.add_noise(noise_type='white')

for i in range(5):
    print(f"{sim.observations[0].tod[0][i]:.5e}")

6.26181e-05
3.14752e-04
2.56129e-04
2.21924e-04
1.07646e-04

The input sky can be generated using the function Simulation.get_sky() which produces sky maps for all the detectors of a given observation based on an instance of SkyGenerationParams. These methods are MPI-compatible, distributing inputs based on the job’s detector configuration without requiring broadcast operations.

API reference

litebird_sim.scan_map.compute_signal_for_one_sample(T, Q, U, co, si, gamma)

Bolometric equation

litebird_sim.scan_map.compute_signal_generic_hwp_for_one_sample(Stokes, Vpol, Rhwp, Mhwp, Rtel)

Bolometric equation for generic HWP Mueller matrix (1,0,0,0) x Mpol x Rpol x Rhwp^T x Mhwp x Rhwp x Rtel x Stokes

litebird_sim.scan_map.rot_matrix(mat, angle)

litebird_sim.scan_map.scan_map(tod, pointings, maps: ~litebird_sim.maps_and_harmonics.HealpixMap | dict[str, ~litebird_sim.maps_and_harmonics.HealpixMap | ~litebird_sim.input_sky.SkyGenerationParams] | ~litebird_sim.maps_and_harmonics.SphericalHarmonics | dict[str, ~litebird_sim.maps_and_harmonics.SphericalHarmonics | ~litebird_sim.input_sky.SkyGenerationParams], pol_angle_detectors: ~numpy.ndarray | None = None, pol_eff_detectors: ~numpy.ndarray | None = None, hwp: ~litebird_sim.hwp.HWP | None = None, hwp_angle: ~numpy.ndarray | None = None, mueller_hwp: ~numpy.ndarray | None = None, input_names: str | None = None, interpolation: str | None = '', pointings_dtype=<class 'numpy.float64'>, nthreads: int = 0)

Scan a sky map and fill time-ordered data (TOD) based on detector observations.

This function modifies the values in tod by adding the contribution of the bolometric equation given a T/Q/U description of the sky. The pointings argument must be a DxNx2 matrix containing the pointing information, where D is the number of detectors for the current observation and N is the size of the tod array.

Supported sky descriptions

• HealpixMap

• SphericalHarmonics

• dict[str, HealpixMap]

• dict[str, SphericalHarmonics]

In the dictionary case, the key selected for each detector is taken from input_names[detector_idx].

The sky coordinate system is taken from the coordinates attribute of the selected object (either HealpixMap or SphericalHarmonics). If coordinates is None, Galactic coordinates are assumed and a warning is issued.

param tod:

Time-ordered data (TOD) array of shape (n_detectors, n_samples) that will be filled with the simulated sky signal.

type tod:

np.ndarray

param pointings:

Pointing information for each detector. If an array, it should have shape (n_detectors, n_samples, 2 or 3), where the first two entries are (theta, phi) in radians, and the optional third entry is the telescope orientation. If a callable, it should return pointing data when passed a detector index.

type pointings:

np.ndarray or callable

param maps:

dict[str, SphericalHarmonics] Sky model to be scanned. In dictionary form, the keys must match the entries of input_names.

type maps:

HealpixMap | SphericalHarmonics | dict[str, HealpixMap] |

param pol_angle_detectors:

Polarization angles of detectors in radians. If None, all angles are set to zero.

type pol_angle_detectors:

np.ndarray or None, default=None

param pol_eff_detectors:

Polarization efficiency of detectors. If None, all detectors have unit efficiency.

type pol_eff_detectors:

np.ndarray or None, default=None

param hwp:

Half-wave plate (HWP) model. If None, HWP effects are ignored unless the Observation object contains HWP data.

type hwp:

HWP, optional

param hwp_angle:

Half-wave plate angles of an external HWP object.

type hwp_angle:

np.ndarray or None, default=None

param mueller_hwp:

Mueller matrices for a non-ideal HWP.

type mueller_hwp:

np.ndarray or None, default=None

param input_names:

Per-detector keys to select entries in maps when maps is a dictionary.

type input_names:

array-like of str or None, default=None

param interpolation:

Currently unused here for the harmonic case; real-space maps are sampled by nearest-neighbour (Healpix ang2pix).

type interpolation:

str or None, default=””

param pointings_dtype:

Data type for pointings generated on the fly.

type pointings_dtype:

dtype, optional

param nthreads:

Number of threads to use for convolution. If set to 0, all available CPU cores will be used.

type nthreads:

int, default=0

raises TypeError:

If maps is None at this level, or has an unsupported type.

raises AssertionError:

If tod and pointings shapes are inconsistent.

Notes

• The function modifies tod in place by adding the scanned sky signal.

• If mueller_hwp is provided, a full HWP Mueller matrix transformation is applied.

• Polarization angles are corrected based on telescope orientation and HWP effects.

• This function is crucial for simulating realistic observations in CMB and astrophysical experiments.

litebird_sim.scan_map.scan_map_for_one_detector(tod_det, input_T, input_Q, input_U, pol_angle_det, pol_eff_det)

litebird_sim.scan_map.scan_map_generic_hwp_for_one_detector(tod_det, input_T, input_Q, input_U, orientation_telescope, pol_angle_det, pol_eff_det, hwp_angle, mueller_hwp)

litebird_sim.scan_map.scan_map_in_observations(observations: ~litebird_sim.observations.Observation | list[~litebird_sim.observations.Observation], maps: ~litebird_sim.maps_and_harmonics.HealpixMap | dict[str, ~litebird_sim.maps_and_harmonics.HealpixMap | ~litebird_sim.input_sky.SkyGenerationParams] | ~litebird_sim.maps_and_harmonics.SphericalHarmonics | dict[str, ~litebird_sim.maps_and_harmonics.SphericalHarmonics | ~litebird_sim.input_sky.SkyGenerationParams] | None = None, pointings: ~numpy.ndarray | list[~numpy.ndarray] | None = None, hwp: ~litebird_sim.hwp.HWP | None = None, component: str = 'tod', pointings_dtype: type[~typing.Any] | ~numpy.dtype[~typing.Any] | ~numpy._typing._dtype_like._HasDType[~numpy.dtype[~typing.Any]] | ~numpy._typing._dtype_like._HasNumPyDType[~numpy.dtype[~typing.Any]] | tuple[~typing.Any, ~typing.Any] | list[~typing.Any] | ~numpy._typing._dtype_like._DTypeDict | str = <class 'numpy.float64'>, save_tod: bool = True, apply_non_linearity: bool = False, add_2f_hwpss: bool = False, mueller_phases: dict | None = None, integrate_in_band: bool = False, nthreads: int | None = None)

Scan a sky map and fill time-ordered data (TOD) for a set of observations.

This is a wrapper around the scan_map() function (and, for the harmonic-expansion case, fill_tod()) that applies to the TOD stored in observations and the pointings stored in pointings. The two types can either be a single Observation instance and a NumPy matrix, or a list of observations and a list of NumPy matrices; in the latter case, they must have the same number of elements.

The field maps can be:

• a single HealpixMap instance, used for all detectors;

• a single SphericalHarmonics instance, used for all detectors;

• a dictionary mapping detector or channel names (strings) to HealpixMap or SphericalHarmonics objects.

In the harmonic case (when the underlying routines see SphericalHarmonics objects), the code performs interpolation using the ducc0 spherical harmonics backend. For real-space inputs (HealpixMap), no interpolation in harmonic space is performed.

By default, the signal is added to Observation.tod. If you want to add it to some other field of the Observation class, use component:

for cur_obs in sim.observations:
    # Allocate a new TOD for the sky signal alone
    cur_obs.sky_tod = np.zeros_like(cur_obs.tod)

# Ask `scan_map_in_observations` to store the sky signal
# in `observations.sky_tod`
scan_map_in_observations(sim.observations, ..., component="sky_tod")

Parameters:

• observations (Observation or list[Observation]) – One or more Observation objects containing detector names, pointings, and TOD data, to which the computed sky signal will be added.

• maps (HealpixMap, SphericalHarmonics, dict[str, HealpixMap] or dict[str, SphericalHarmonics], optional) – Sky description. If a dictionary, keys should match detector or channel names, and values should be HealpixMap or SphericalHarmonics instances. If a single object is provided, the same sky is used for all detectors. If None, the function attempts to read sky from the first observation.

• pointings (np.ndarray or list[np.ndarray], optional) – Pointing matrices associated with the observations. If None, the function extracts pointing information from the Observation objects.

• hwp (HWP, optional) – Half-wave plate (HWP) model. If None, HWP effects are ignored unless the Observation object itself contains HWP data.

• component (str, default="tod") – The TOD component in the Observation object where the computed signal will be stored.

• pointings_dtype (dtype, optional) – Data type for pointings generated on the fly. If the pointing is passed or already precomputed this parameter is ineffective. Default is np.float64.

• apply_non_linearity (bool, optional) – (For the harmonics expansion case) applies the coupling of the non-linearity systematics with hwp_sys.

• add_2f_hwpss (bool, optional) – (For the harmonics expansion case) adds the 2f HWP synchronous signal to the TOD.

• mueller_phases (dict or np.ndarray, optional) – (For the harmonics expansion case) the non-ideal phases for the Mueller matrix elements. When None is given, temporary values from Patanchon et al. [2021] are used.

• integrate_in_band (bool, default=False) – Whether to integrate the signal over the detector’s frequency band. Only implemented for the Jones formalism and explicit harmonics expansion case.

• nthreads (int, default=None) – Number of threads to use in the convolution. If None, the function reads from the OMP_NUM_THREADS environment variable.

Raises:

• ValueError – If a dictionary maps does not contain the required detector or channel keys.

• AssertionError – If the number of observations and pointings do not match. If tod and pointings shapes are inconsistent.

Notes

• This function modifies observations in place by adding the computed sky signal to the specified component field.

• If pointings is None, the function attempts to extract them from Observation objects. If the pointing is generated on the fly, pointings_dtype specifies its type.

• If an HWP model is provided, the function computes HWP angles and applies the corresponding Mueller matrices.

• This function supports both single observations and lists of observations, handling each one separately.

litebird_sim.scan_map.vec_polarimeter(angle, gamma)

litebird_sim.scan_map.vec_stokes(stokes, T, Q, U)