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.

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

hwp_radpsec = np.pi / 8
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 MBS module provided by litebird_sim)
maps = np.ones((3, npix))
in_map = {"Detector": maps, "Coordinates": lbs.CoordinateSystem.Ecliptic}

# Here scan the map and fill tod
lbs.scan_map_in_observations(
    obs,
    in_map,
    input_map_in_galactic = False,
)

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 input maps to scan can be either included in a dictionary with the name of the channel or the name of the dectector as keyword (the routines described in Synthetic sky maps already provied the inputs in the correct format), or a numpy array with shape (3, n_pixels).

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. If the input map is ecliptic coordinates set input_map_in_galactic to False. 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 routine provides an on-the-fly interpolation of the input maps. This option is available through the argument interpolation which specifies the type of TOD interpolation (”” for no interpolation, “linear” for linear interpolation). Default: no interpolation.

The low level function, scan_map(), allows a more refined handling of the inputs.

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()

sky_signal = np.ones((3,12*nside*nside))*1e-4

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}")

4.14241e-04
5.46700e-05
3.03378e-04
6.13975e-05
4.72613e-05

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: Dict[str, ndarray], pol_angle_detectors: None | ndarray = None, pol_eff_detectors: None | ndarray = None, hwp_angle: None | ndarray = None, mueller_hwp: None | ndarray = None, input_names: str | None = None, input_map_in_galactic: bool = True, interpolation: str | None = '')#

Scan a map filling time-ordered data

This function modifies the values in tod by adding the contribution of the bolometric equation given a list of TQU maps maps. The pointings argument must be a DxNx2 matrix containing the pointing information, where D is the number of detector for the current observation and N is the size of the tod array. pol_angle is the array of size DxN containing the polarization angle in radiants. input_names is an array containing the keywords that allow to select the proper input in maps for each detector in the TOD. If input_map_in_galactic is set to False the input map is assumed in ecliptic coordinates, default galactic. The interpolation argument specifies the type of TOD interpolation (”” for no interpolation, “linear” for linear interpolation)

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: Observation | List[Observation], maps: Dict[str, ndarray], pointings: ndarray | List[ndarray] | None = None, hwp: HWP | None = None, input_map_in_galactic: bool = True, component: str = 'tod', interpolation: str | None = '')#

Scan a map filling time-ordered data

This is a wrapper around the scan_map() function that applies to the TOD stored in observations and the pointings stored in pointings. The two types can either bed a 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 must either be a dictionary associating the name of each detector with a (3, NPIX) array containing the three I/Q/U maps or a plain (3, NPIX) array. In the latter case, the I/Q/U maps will be used for all the detectors.

The coordinate system is usually specified using the key Coordinates in the dictionary passed to the maps argument, and it must be an instance of the class CoordinateSystem. If you are using a plain NumPy array instead of a dictionary for maps, you should specify whether to use Ecliptic or Galactic coordinates through the parameter input_map_in_galactic. If maps["Coordinates"] is present, it must be consistent with the value for input_map_in_galactic; if not, the code prints a warning and uses the former.

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 `add_noise_to_observations` to store the noise
# in `observations.sky_tod`
scan_map_in_observations(sim.observations, …, component="sky_tod")
litebird_sim.scan_map.vec_polarimeter(angle, gamma)#
litebird_sim.scan_map.vec_stokes(stokes, T, Q, U)#