Time Ordered Simulations

The LiteBIRD simulation framework intends to offer the capability to simulate many effects in detector timestreams. This page will detail these simulated effects as they become available. The timestreams are stored in the tod field of Observation objects. These arrays have shape (n_detectors, n_samples), where detectors are indexed in the same order as the Observation.detector_global_info array. All the time streams inside the framework use Kelvin as the default unit.

Filling TOD with signal

The framework provides scan_map(), a routine which scans an input map accordingly to the scanning strategy and fills the detector timestreams. 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])

# Compute the pointing information
sim.compute_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 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 obs 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 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.

Adding Noise

The ability to add noise to your detector timestreams is supported through the function add_noise_to_observations() and the low-level versions add_noise(), add_white_noise(), and add_one_over_f_noise().

Here is a short example that shows how to add noise:

import litebird_sim as lbs
import numpy as np

# Create a simulation lasting 100 seconds
sim = lbs.Simulation(
    base_path='./output',
    start_time=0,
    duration_s=100,
    random_seed=12345,
)

# Create a detector object
det = lbs.DetectorInfo(
  net_ukrts=100,
  sampling_rate_hz=10
)

obs = sim.create_observations(detectors=[det])

# Here we add white noise using the detector
# noise parameters from the `det` object.
# We use the random number generator provided
# by `sim`, which is initialized with the
# seed we passed to the Simulation constructor
# to ensure repeatability.
lbs.noise.add_noise_to_observations(obs, 'white', random=sim.random)

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

5.65263e-04
-1.54522e-04
3.42276e-04
-1.42274e-04
-1.71110e-04
8.72188e-05
-1.23400e-04
-6.99311e-05
6.58389e-05
5.51306e-04

Note that we pass sim.random as the number generator to use. This is a member variable that is initialized by the constructor of the class Simulation, and it is safe to be used with multiple MPI processes as it ensures that each process has its own random number generator with a different seed. You can also pass another random number generator, as long as it has the normal method. More information on the generation of random numbers can be found in Random numbers in litebird_sim.

To add white noise using a custom white noise sigma, in µK, we can call the low level function directly:

import litebird_sim as lbs

sim = lbs.Simulation(
    base_path='./output',
    start_time=0,
    duration_s=100,
    random_seed=12345,
)

det = lbs.DetectorInfo(
  net_ukrts=100,
  sampling_rate_hz=10,
)

obs = sim.create_observations(detectors=[det])

custom_sigma_uk = 1234
lbs.noise.add_white_noise(obs[0].tod[0], custom_sigma_uk, random=sim.random)

We can also add 1/f noise using a very similar call to the above:

import litebird_sim as lbs

sim = lbs.Simulation(
    base_path='./output',
    start_time=0,
    duration_s=100,
    random_seed=12345,
)

det = lbs.DetectorInfo(
  net_ukrts=100,
  sampling_rate_hz=10,
  alpha=1,
  fknee_mhz=10
)

obs = sim.create_observations(detectors=[det])

# Here we add 1/f noise using the detector noise
# parameters from the detector object
lbs.noise.add_noise_to_observations(obs, 'one_over_f', random=sim.random)

Again, to generate noise with custom parameters, we can either use the low-level function or edit the Observation object to contain the desired noise parameters.

import litebird_sim as lbs
import numpy as np

sim = lbs.Simulation(
    base_path='./output',
    start_time=0,
    duration_s=100,
    random_seed=12345,
)

det = lbs.DetectorInfo(
  net_ukrts=100,
  sampling_rate_hz=10,
  alpha=1,
  fknee_mhz=10,
  fmin_hz=0.001,
)

obs = sim.create_observations(detectors=[det])

custom_sigma_uk = 1234
custom_fknee_mhz = 12.34
custom_alpha = 1.234
custom_fmin_hz = 0.0123

# Option 1: we call the low-level function directly
lbs.noise.add_one_over_f_noise(
    obs[0].tod[0],
    custom_fknee_mhz,
    custom_fmin_hz,
    custom_alpha,
    custom_sigma_uk,
    obs[0].sampling_rate_hz,
    sim.random,
)

# Option 2: we change the values in `obs`
obs[0].fknee_mhz[0] = custom_fknee_mhz
obs[0].fmin_hz[0] = custom_fmin_hz
obs[0].alpha[0] = custom_alpha
obs[0].net_ukrts[0] = (
    custom_sigma_uk / np.sqrt(obs[0].sampling_rate_hz)
)

lbs.noise.add_noise_to_observations(obs, 'one_over_f', random=sim.random)

Warning

Be sure you understand the difference between the noise level in a timestream and the noise level in a map. Although of course the latter depends on the former, the conversion depends on several factors.

A common mistake is use the mission time divided by the number of pixels in the map in a call to add_white_noise(). This is wrong, as the noise level per pixel depends on the overall integration time, which is always less than the mission time because of cosmic ray loss, repointing maneuvers, etc. These effects reduce the number of samples in the timeline that can be used to estimate the map, but they do not affect the noise of the timeline.

Methods of the Simulation class

The class Simulation provides two simple functions that fill with sky signal and noise all the observations of a given simulation. The function Simulation.fill_tods() takes a map and scans it, while the function Simulation.add_noise() adds noise to the timelines. Thanks to these functions 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.compute_pointings()

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

sim.fill_tods(sky_signal)

sim.add_noise(noise_type='white', random=sim.random)

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)

Bolometric equation

litebird_sim.scan_map.scan_map(tod, pointings, pol_angle, maps: Dict[str, ndarray], input_names, 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)

litebird_sim.scan_map.scan_map_in_observations(obs: Observation | List[Observation], maps: Dict[str, ndarray], pointings: ndarray | List[ndarray] | 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 obs 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 `obs.sky_tod`
scan_map_in_observations(sim.observations, …, component="sky_tod")

litebird_sim.noise.add_noise(tod, noise_type: str, sampling_rate_hz: float, net_ukrts, fknee_mhz, fmin_hz, alpha, random, scale=1.0)

Add noise (white or 1/f) to a 2D array of floating-point values

This function sums an array of random number following a white noise model with an optional 1/f component to data, which is assumed to be a D×N array containing the TOD for D detectors, each containing N samples.

The parameter noisetype must either be white or one_over_f; in the latter case, the noise will contain a 1/f part and a white noise part.

Be sure not to include cosmic ray loss, repointing maneuvers, etc., in the value passed as net_ukrts, as these affect the integration time but not the white noise per sample.

The parameter scale can be used to introduce measurement unit conversions when appropriate. Default units: [K].

The parameter random must be specified and must be a random number generator that implements the normal method. You should typically use the random field of a Simulation object for this.

The parameters net_ukrts, fknee_mhz, fmin_hz, alpha, and scale can either be scalars or arrays; in the latter case, their size must be the same as tod.shape[0], which is the number of detectors in the TOD.

litebird_sim.noise.add_noise_to_observations(obs: Observation | List[Observation], noise_type: str, random: Generator, scale: float = 1.0, component: str = 'tod')

Add noise of the defined type to the observations in obs

This class provides an interface to the low-level function add_noise(). The parameter obs can either be one Observation instance or a list of observations, which are typically taken from the field observations of a Simulation object. Unlike add_noise(), it is not needed to pass the noise parameters here, as they are taken from the characteristics of the detectors saved in obs. The parameter random must be specified and must be a random number generator that implements the normal method. You should typically use the random field of a Simulation object for this.

By default, the noise 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 noise alone cur_obs.noise_tod = np.zeros_like(cur_obs.tod)

# Ask add_noise_to_observations to store the noise # in obs.noise_tod add_noise_to_observations(sim.observations, …, component=”noise_tod”)

See add_noise() for more information.

litebird_sim.noise.add_one_over_f_noise(data, fknee_mhz: float, fmin_hz: float, alpha: float, sigma: float, sampling_rate_hz: float, random)

Adds a 1/f noise timestream with the given f knee and alpha to data To be called from add_noise_to_observations

Parameters:

• data – 1-D numpy array

• fknee_mhz – knee frequency in mHz

• fmin_hz – kmin frequency for high pass in Hz

• alpha – low frequency spectral tilt

• sigma – the white noise level per sample. Be sure not to include cosmic ray loss, repointing maneuvers, etc., as these affect the integration time but not the white noise per sample.

• sampling_rate_hz – the sampling frequency of the data

• random – a random number generator that implements the normal method. You should typically use the random field of a Simulation object for this. It must be specified

litebird_sim.noise.add_white_noise(data, sigma: float, random)

Adds white noise with the given sigma to the array data.

To be called from add_noise_to_observations.

Parameters:

• data – 1-D numpy array

• sigma – the white noise level per sample. Be sure not to include cosmic ray loss, repointing maneuvers, etc., as these affect the integration time but not the white noise per sample.

• random – a random number generator that implements the normal method. You should typically use the random field of a Simulation object for this. It must be specified

litebird_sim.noise.build_one_over_f_model(ft, freqs, fknee_mhz, fmin_hz, alpha, sigma)

litebird_sim.noise.nearest_pow2(data)

returns the next largest power of 2 that will encompass the full data set

data: 1-D numpy array

litebird_sim.noise.rescale_noise(net_ukrts: float, sampling_rate_hz: float, scale: float)