Scanning strategy

The LiteBIRD Simulation Framework provides a set of tools to simulate the orbit of the spacecraft and to compute the directions where each detector is looking at the sky, as a function of time. The time stream of directions and orientations of each detector is usually called pointing information, and we’ll consistently use this jargon in the documentation.

Note that this chapter only deals with the direction some detector is looking at, but the actual position/velocity of the spacecraft is not needed to do this calculation. The framework provide other facilities to compute this information, and they are described in Dipole anisotropy.

This chapter provides an in-depth explanation about how to use the facilities provided by the framework to compute the pointing information for any detector in one of the focal planes.

The spacecraft’s motion

In the case of a space mission like LiteBIRD, the motion of the spacecraft and its structure decide where each detector is looking at each time. The following video shows qualitatively which kind of motion is simulated by our framework:

(The animation was created using POV-Ray, and the files used to create it are available on GitHub .)

You can see that there are two rotations in the animation: the spacecraft spins quickly around its spin axis (grey axis), but this axis does not stand still: it performs a precession around the blue axis, which represents the Sun-Earth direction. (You should imagine the Earth on the left, and the Sun on the very far left.)

Note that the detectors are not necessarily aligned with the spin axis; in fact, the animation shows the actual direction of observation for two different detectors as two red and green lines: you can see that they are looking quite at two opposite sides of the spin axis. Every detector looks along its own direction, but detectors belonging to the same instrument (e.g., LFT) look not far away from each other; it is customary to express their pointing directions relative to an «average» direction, called the boresight direction, which is the main optical axis of the instrument. Since in LiteBIRD there are three instruments (LFT, MFT, HFT), there should be three boresight directions; however, MFT and HFT share the same telescope, and thus it’s customary to show only one boresight for both. This is the true meaning of the red and green axes in the video above: the red axis represents the «average» direction where LFT detectors are looking at, and the green axis is the same for MFT/HFT.

The animation does not show a third rotation happening, which is the revolution of the spacecraft around the Sun, taking one year to complete. (Including it in the video would have been useless, as it is really slow when compared with the spin and precession!). This means that the motion of the spacecraft is the composition of three rotations:

1. Rotation of the spacecraft around the spin axis (grey line);

2. Rotation (precession) of the spin axis around the Sun-Earth axis (blue line);

3. Yearly rotation of the Sun-Earth axis around the Sun.

If you think about it, you will realize that the kinematics of this motion can be fully described by the following quantities:

• The angle between the spin axis and the boresight direction(s), usually called β;

• The angle between the spin axis and the Sun-Earth axis, usually called α.

• The speed of the rotation of the boresight direction around the spin axis;

• The speed of the precession around the Sun-Earth axis, which is usually slower than the rotation speed;

They are sketched in the following diagram:

These parameters define the so-called scanning strategy, i.e., the way the instruments observe the sky during the mission’s lifetime. The LiteBIRD Simulation Framework provides all the tools necessary to simulate the composition of these rotations, and it can produce pointing information from the synthetic description of the scanning strategy. Here is a full example, using the scanning strategy proposed for CORE ([3]), which is qualitatively similar to what is going to be used for LiteBIRD:

import litebird_sim as lbs
import astropy.units as u
import numpy as np

sim = lbs.Simulation(
    start_time=0,
    duration_s=60.0,
    description="Simple simulation",
    random_seed=12345,
)

# We now simulate the motion of the spacecraft over a time span
# of one minute (specified in the `duration_s` parameter above).
sim.set_scanning_strategy(
    scanning_strategy=lbs.SpinningScanningStrategy(
        spin_sun_angle_rad=np.deg2rad(30), # CORE-specific parameter
        spin_rate_hz=0.5 / 60.0,     # Ditto
        # We use astropy to convert the period (4 days) in
        # seconds
        precession_rate_hz=1.0 / (4 * u.day).to("s").value,
    )
)

# Here we specify the β angle of the focal plane of
# the instrument
sim.set_instrument(
    lbs.InstrumentInfo(
        name="core",
        spin_boresight_angle_rad=np.deg2rad(65),
    ),
)

# The motion of the spacecraft is now encoded in a set of quaternions,
# in the field `sim.spin2ecliptic_quats`. We use it to produce the
# pointing information for a fake boresight detector `det`, belonging
# to the instrument `core` (unlike LiteBIRD, CORE had only one focal
# plane and one instrument)
det = lbs.DetectorInfo(name="foo", sampling_rate_hz=10)

# By default, `create_observations` creates just *one* observation
obs, = sim.create_observations(detectors=[det])

# Compute the pointings at the same sampling frequency as the
# TOD (10 Hz, see the variable `det` above)
sim.compute_pointings()

pointings = obs.pointings

print("Shape:", pointings.shape)
print("Pointings:")
print(np.array_str(pointings, precision=3))

Shape: (1, 600, 2)
Pointings:
[[[ 2.182  0.   ]
  [ 2.182 -0.006]
  [ 2.182 -0.012]
  ...
  [ 0.089 -2.967]
  [ 0.088 -3.021]
  [ 0.087 -3.075]]]

All the details in this code are explained in the next sections, so for now just keep in mind the overall shape of the code:

1. Once the duration of the simulation (one minute in the example above) is set, we call the method Simulation.set_scanning_strategy(), which forces the framework to compute how the orientation of the spacecraft with respect to the sky sphere evolves with time. This method produces a set of quaternions, which encode the result of the composition of all the rotations (spin, precession, revolution around the Sun) described above; these quaternions are saved in the spin2ecliptic_quats field of the sim class.

2. When the simulation code needs to determine where a detector is pointing to (the detector det in our example), the quaternions are used to retrieve (1) the coordinates on the Sky sphere, and (2) the orientation angle (ψ). Both quantities are computed in the Ecliptic reference frame using the sampling rate of the detector, which in our example is 10 Hz (i.e., ten samples per second). In the example above, this is done by the function get_pointings().

3. The function get_pointings() returns a (N, 3) matrix, where the first column contains the colatitude \(\theta\), the second column the longitude \(\phi\), and the third column the orientation angle \(\psi\), all expressed in radians.

Computing the spacecraft’s orientation

To compute where a detector is looking at the sky sphere, there is a number of transformations that need to be carried out:

We start from the detector’s reference frame, which assumes that the main beam of the radiation pattern is aligned with the z axis, and that the beam of the detector is oriented using the x axis as the reference axis. (In other words, the x axis provides a reference frame for asymmetric beams.)

The next reference frame is the boresight, and to convert from the detector’s reference frame to the boresight there is a rotation, which is encoded in a quaternion that is saved in the IMO.

Next, we move from the reference frame of the boresight to that of the spacecraft. The information about the placement of the boresight with respect to the spin axis is encoded in the class InstrumentInfo. After this transformation, the spin axis is aligned with the z axis.

The next transformation goes from the spacecraft’s to the Ecliptic reference frame; the Ecliptic is on the xy plane, and the z axis points towards the Ecliptic North Pole. In this case, the framework provides two ways to compute the transformation:

1. The revolution of the Earth around the Sun is modelled using a plain circular motion, and the starting position is arbitrary; this mode is triggered whenever the time is tracked using floating-point numbers (i.e., the parameter start_time in the constructor of Simulation is a float).

2. The proper motion of the Earth around the Sun is computed using ephemeridis tables; the calculation is much slower, but the result enables to properly simulate the observation of moving objects in the Solar System, like planets or asteroids. In this case, the parameter start_time must be an instance of the class astropy.time.Time. In the example above, we would enable the computation of proper Earth’s motion with the following small change:

   import astropy.time
   
   sim = lbs.Simulation(
       # Turn on full computation of the Earth's orbit around the Sun
       start_time=astropy.time.Time("2020-01-01"),
       duration_s=60.0,
       description="Simple simulation",
       random_seed=12345,
   )
   

You should compute the proper motion of the Earth around the Sun only if you absolutely need to, as the computation can be 10÷100 times slower.

From quaternions to detector pointings

To compute the pointing information for a detector, the quaternions computed through the call to Simulation.set_scanning_strategy() are not enough, as they only tell how to convert a vector from the spin axis reference frame to the Ecliptic reference frame. We need two more quaternions that tell how to convert from the reference frame of the detector to that of the spin axis:

1. The first quaternion describes how the detector reference frame (with the z axis aligned with the main axis of the radiation pattern) can be converted to the reference frame of the focal plane (with the z axis aligned with the boresight). This information is included in the IMO and is properly initialized if you call DetectorInfo.from_imo(). If you do not specify any quaternion, the constructor for DetectorInfo will assume that the detector is looking at the boresight, and it will thus use the quaternion \((0 0 0 1)\); this is the case of the simple example we presented above.

2. The second quaternion describes how to convert the reference frame of the focal plane (with the z axis aligned with the boresight) to the reference frame of the spacecraft (where the z axis is aligned with its spin axis). This quaternion is stored in the field bore2spin_quat of the class InstrumentInfo, and it is initialized when you call the method Simulation.set_scanning_strategy().

The LiteBIRD Simulation Framework recomputes the orientation of the spacecraft with a regular spacing in time (the default is one minute). However, pointings need to be known at the same sampling frequency used by the detector, which is usually much higher than the frequency used to compute the quaternions (in our example above, the sampling frequency of detector det was 10 Hz, but the sampling frequency of the quaternions was 1/60 Hz). Since the framework uses quaternions to encode the orientation of the spacecraft, oversampling them to the sampling frequency of the detector is just a matter of applying a spherical linear interpolation (abbr. slerp), according to the following figure:

To be sure to include an additional quaternion after the last sample, like in the figure above, the framework provides the static method ScanningStrategy.optimal_num_of_quaternions(), which computes how many quaternions need to be calculated to cover some time span with a given interval between quaternions. For instance, if our simulation lasts 100 s and we want one quaternion every minute, then the expected number of quaternions to be computed is 3: one at \(t = 0\), one at \(t = 60\,\mathrm{s}\), and one at \(t = 120\,\mathrm{s}\), so that the latter two can be interpolated for the samples in the range \(60\,\mathrm{s} \leq t \leq 100\,\mathrm{s}\). Here is how the function works:

print(lbs.ScanningStrategy.optimal_num_of_quaternions(
    time_span_s=100,
    delta_time_s=60,
))

3

When using MPI, the relatively small size in memory of the quaternions (the thick black lines in the figure) enables the framework to keep a duplicate of the list in all the MPI processes. This is unlike what happens with the data in TODs (the thin gray lines), which are split in several blocks inside the Observation class.

Note

Slerp assumes a rotation with constant angular speed and axis between consecutive quaternions, and thus it only approximates the true composition of all the rotations (spin, precession, revolution around the Sun) that we have discussed above. However, don’t forget that the real spacecraft will follow a scanning strategy that will be more complex than the one described by our geometrical model, because of many non-idealities that are unavoidable in a spacecraft. The approximation of the «slerp» operation is thus unlikely to be relevant.

Once all the quaternions have been computed at the proper sampling rate, the direction of the detector on the sky and its orientation angle are computed as follows:

• The direction is the vector \(\vec d = R \hat e_z\), where \(R\) is the overall rotation from the detector’s reference frame to the Ecliptic reference frame, and \(\hat e_z = (0, 0, 1)\) is the one-length vector aligned with the z axis.

• The orientation angle is given by the angle between the North direction passing through the vector \(\vec d\) (i.e., along the meridian of \(\vec d\)) and the vector \(\vec p = R \hat e_x\), where \(R\) is the same as above and \(\hat e_x = (1, 0, 0)\), as shown in the following figure (note that \(\hat e_x\) has been drawn twice because the one in the upper side represents the orientation direction in the detector’s reference frame):

  

The purpose of the method Simulation.compute_pointings(), used in the example at the beginning of this chapter, is to call get_det2ecl_quaternions() to compute the quaternions at the same sampling frequency as the scientific datastream, and then to apply the two definitions above to compute the direction and the orientation angle.

How the boresight is specified

As LiteBIRD includes three focal planes and two telescopes, the specification of the boresight requires some care. In [3] and [4], the boresight direction is encoded using just one number, the angle between the boresight and the spin axis. However, both papers deal with spacecrafts hosting only one focal plane.

The full orientation of the boresight direction is specified using three angles:

1. The ψ angle encodes the rotation of the focal plane with respect to the boresight direction itself, and it is ideally 0°;

2. The angle between the boresight direction and the spin axis is usually notated with the symbol β (among the three, this is the most important number: it’s 65° for CORE, 69° for PICO);

3. Finally, the boresight can be rotated by an angle φ around the spin axis: this is important only when you have more than one focal plane. For LiteBIRD, \(\phi_\text{LFT} - \phi_\text{MHFT} \approx 180^\circ\).

Interpretation of pointings

With «pointing», we refer to two different concept:

1. The direction where the detector is looking at;

2. The orientation of the detector while it’s looking at the sky.

The direction can be encoded either as a one-length vector (x, y, z) or as a couple of angles; the LiteBIRD simulation framework adopts the second option to save memory, and it encodes directions using the colatitude (i.e., 90° minus the latitude) and the longitude, commonly indicated with the letters θ (colatitude) and φ (longitude).

The orientation of the detector (second point above) can be expressed either as a vector that is tangent to the sky sphere, or as an angle calculated with respect to the meridian/parallel going through the point the detector is looking at. Again, to reduce memory usage, our framework only encodes the angle.

The method Simulation.compute_pointings() stores the pointings of the \(n_d\) detectors kept in the field pointings of the Observation; they are laid out in memory as a \((n_d, N, 2)\) matrix, where \(N\) is the number of samples in the timeline, and the last dimension holds the colatitude and longitude (in radians). The orientation angle is kept in Observation.psi. Let’s visualize the position of these pointings on a Healpix map:

import healpy, numpy as np
import matplotlib.pylab as plt

nside = 64
pixidx = healpy.ang2pix(nside, pointings[0, :, 0], pointings[0, :, 1])
m = np.zeros(healpy.nside2npix(nside))
m[pixidx] = 1
healpy.mollview(m)

Custom scanning strategies

In this section we explain how scanning strategies different from the nominal, «spinning» strategy can be modelled. You will need to understand the functions provided by the framework to deal with quaternions.

The framework uses a right-handed coordinate system, like the one shown in figure:

where the grey arrows indicate the verse of positive rotations (they follow the usual right-hand rule: point your thumb along the axis, and the other fingers will point towards the positive direction of the rotation).

A few words about quaternions

To describe a rotation in 3D space, there are several choices: Euler angles, rotation matrices, quaternions, etc. Each of these systems has its own share of advantages and disadvantages: for instance, rotation matrices are handy when you have a vector and want to rotate it, as it’s just a matter of doing a matrix-vector multiplication. Quaternions are more complicated on this regard, but they offer a mathematical operation called slerp (shorthand for spherical linear interpolation), which is incredibly useful and not available with other representations, like rotation matrices. We assume that the reader knows what quaternion are and their mathematical properties; if you are not, be sure to read the book Visualizing quaternions, by Andrew J. Hanson (Elsevier, 2006, ISBN-0080474772) and the provocative essay by Marc ten Bosch, Let’s remove Quaternions from every 3D engine.

The LiteBIRD simulation framework models quaternions using the convention \((v_x, v_y, v_z, w)\); be aware that some textbooks use the order \((w, v_x, v_y, v_z)\). As the framework uses quaternions only to model rotations, they all must obey the relation \(v_x^2 + v_y^2 + v_z^2 + w^2 = 1\) (normalized quaternions), which is a property satisfied by rotation quaternions.

Python functions for quaternions

The LiteBIRD Simulation Framework provides three functions, quat_rotation_x(), quat_rotation_y(), quat_rotation_z() to compute simple rotation quaternions; they return the normalized quaternion representing a rotation by an angle \(\theta\) around one of the three axis x, y, and z:

import litebird_sim as lbs
import numpy as np

def print_quaternion(q):
    print("{:.3f} {:.3f} {:.3f} {:.3f}".format(*q))

print("Rotation by π/3 around x:")
print_quaternion(lbs.quat_rotation_x(theta_rad=np.pi/3))
print("Rotation by π/3 around y:")
print_quaternion(lbs.quat_rotation_y(theta_rad=np.pi/3))
print("Rotation by π/3 around z:")
print_quaternion(lbs.quat_rotation_z(theta_rad=np.pi/3))

Rotation by π/3 around x:
0.500 0.000 0.000 0.866
Rotation by π/3 around y:
0.000 0.500 0.000 0.866
Rotation by π/3 around z:
0.000 0.000 0.500 0.866

There are two functions that implement in-place multiplication of quaternions: quat_right_multiply() performs the calculation \(r \leftarrow r \times q\), and quat_left_multiply() performs the calculation \(r \leftarrow q \times r\) (where \(\leftarrow\) indicates the assignment operation):

quat = np.array(lbs.quat_rotation_x(np.pi / 3))
lbs.quat_left_multiply(quat, *lbs.quat_rotation_z(np.pi / 2))
print("Rotation by π/3 around x and then by π/2 around z:")
print_quaternion(quat)

Rotation by π/3 around x and then by π/2 around z:
0.354 0.354 0.612 0.612

Note the syntax for quat_left_multiply(): you are supposed to pass the four components of the quaternion \(q\) as separate arguments, and thus we need to prepend the call to lbs.quat_rotation_z with * to expand the result (a 4-element tuple) into the four parameters required by quat_left_multiply(). The reason for this weird syntax is efficiency, as this kind of function call can be easily optimized by Numba (which is used extensively in the code).

Finally, the framework provides the function rotate_vector(), which applies the rotation described by a normalized quaternion to a vector. There are faster versions in rotate_x_vector(), rotate_y_vector(), and rotate_z_vector() that rotate the three basis vectors (1, 0, 0), (0, 1, 0), and (0, 0, 1). The functions all_rotate_vectors(), all_rotate_x_vectors(), all_rotate_y_vectors(), and all_rotate_z_vectors() can be applied to arrays of quaternions and vectors.

A simple scanning strategy

We are now ready to discuss how to implement other types of scanning strategies. There are plenty of reasons why one would like to go beyond the class SpinningScanningStrategy:

1. You want to study the effect of non-idealities, like second-order effects caused by contractions/dilations in the mechanical structure of the telescope that make the angle between the spin axis and the boresight vary with time.

2. You are thinking about how to make dedicated observations of some celestial source (e.g., the Crab Nebula) for the purpose of calibrating the instruments.

To define a new scanning strategy, we define a descendeant of the ScanningStrategy class, an Abstract Base Class (ABC); the only method that must be defined is ScanningStrategy.generate_spin2ecl_quaternions(), which takes as inputs the start time, the length of the simulation, and the time interval to be used between consecutive quaternions. The method must return an instance of the Spin2EclipticQuaternions, containing the computed sequence of quaternions.

We’ll code here a very simple scanning strategy, which does not involve anything fancy: the spacecraft will just spin around the Sun-Earth axis, and the boresight direction will be along the same spin axis. Thus, the boresight detector is going to see only the points along the Ecliptic plane. This scanning strategy is scientifically useless, but it’s simple enough to be implemented in a few lines of code:

1. The transformation from boresight to the spin axis reference frame is the identity;

2. There is no precession of the spin axis; therefore, the latter stays on the Ecliptic axis;

3. The only rotation is caused by the revolution of the Sun-Earth axis around the Sun, which is implemented as a rotation on the xy plane, i.e., around the z axis.

The following code implements our mock scanning strategy:

class SimpleScanningStrategy(lbs.ScanningStrategy):
    def generate_spin2ecl_quaternions(
        self, start_time, time_span_s, delta_time_s,
    ):
        # Compute how many quaternions are needed to cover
        # the time interval specified by "time_span_s"
        num_of_quaternions = (
            lbs.ScanningStrategy.optimal_num_of_quaternions(
                time_span_s=time_span_s,
                delta_time_s=delta_time_s,
            )
        )

        # Make room for the quaternions
        spin2ecliptic_quats = np.empty((num_of_quaternions, 4))

        # We compute the times when the quaternions need to be
        # calculated. Note that ScanningStrategy returns two
        # arrays ("time" and "time_s"), but we neglect the second
        # because we don't need it in this very simple case
        (time, _) = lbs.ScanningStrategy.get_times(
            start_time=start_time,
            delta_time_s=delta_time_s,
            num_of_quaternions=num_of_quaternions,
        )

        # Compute the angle on the Ecliptic plane between the x
        # axis and the Sun-Earth axis, possibly using AstroPy
        sun_earth_angles_rad = (
            lbs.calculate_sun_earth_angles_rad(time)
        )

        # This code is *not* optimized: in a real-world case,
        # you'll probably want to use Numba instead of the
        # following "for" loop
        for i in range(num_of_quaternions):
            # Rotate by 90° around the y axis (move the boresight
            # to the Ecliptic xy plane)
            spin2ecliptic_quats[i, :] = lbs.quat_rotation_y(np.pi / 2)

            # Simulate the revolution of the spacecraft around
            # the Sun using the angles computed above
            lbs.quat_left_multiply(
                spin2ecliptic_quats[i, :],
                *lbs.quat_rotation_z(sun_earth_angles_rad[i]),
            )

        # Return the quaternions wrapped in an instance of
        # "Spin2EclipticQuaternions"
        return lbs.Spin2EclipticQuaternions(
            start_time=start_time,
            pointing_freq_hz=1.0 / delta_time_s,
            quats=spin2ecliptic_quats,
        )

To test the class SimpleScanningStrategy, we write some code very similar to the example presented at the beginning of this section. However, we cannot run the simulation for just one hour, as it would be not enough to see any change in the pointing direction: the only things that changes as time passes is the position of the Earth on the Ecliptic plane, and it takes 365 days to do one revolution. Therefore, we extend the length of the simulation to 365 days. Of course, there is no need to use an high sampling frequency in our example, so we use here just one sample per day; for the same reason, instead of computing one quaternion every minute, we compute one quaternion every 30 days:

import astropy.units as u
import healpy
import numpy as np

sim = lbs.Simulation(
    start_time=0,
    duration_s=(365 * u.day).to("s").value,
    description="Simple simulation",
    random_seed=12345,
)

sim.set_scanning_strategy(
    scanning_strategy=SimpleScanningStrategy(),
    delta_time_s=(30 * u.day).to("s").value
)

det = lbs.DetectorInfo(
    name="foo",
    sampling_rate_hz=1.0 / ((1.0 * u.day).to("s").value),
)
(obs,) = sim.create_observations(detectors=[det])
pointings = lbs.get_pointings(obs, sim.spin2ecliptic_quats, np.array([det.quat]))

m = np.zeros(healpy.nside2npix(64))
pixidx = healpy.ang2pix(64, pointings[:, 0], pointings[:, 1])
m[pixidx] = 1
healpy.mollview(m)

Here is the result: we’re indeed scanning the Ecliptic plane!

Half Wave Plate

The rotation of the polarization angle induced by a HWP can be included in the returned pointing information by passing an instance of a descendant of the class HWP to the method Simulation.set_hwp(). Here is an example:

import litebird_sim as lbs

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

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=sim.start_time,
    ),
    delta_time_s=60,
)

sim.set_instrument(
    instr = 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,
    )
)

sim.set_hwp(
    lbs.IdealHWP(ang_speed_radpsec=4.084_070_449_666_731),
)

det = lbs.DetectorInfo(
    name="Boresight_detector",
    sampling_rate_hz=1.0,
    quat=[0.0, 0.0, 0.0, 1.0],
)
obs, = sim.create_observations(detectors=[det])

sim.compute_pointings()

This example uses the IdealHWP, which represents an ideal spinning HWP.

Observing point sources in the sky

It is useful to simulate the observation of point sources in the sky, both for a scientific purpose or for instrument calibration. For instance, an important task in the calibration of a CMB space experiment is the estimation of the radiation pattern \(\gamma(\theta, \phi)\) for each detector (sometimes \(\gamma\) is called the beam function). This task can be done through the observation of a bright point source, like one of the outer planets (Mars, Jupiter, Saturn, etc.): assuming that the source is really pointlike and neglecting every other emission from the sky, the response measured by a detector is proportional to the radiation pattern \(\gamma(\theta, \phi)\), where the angles \(\theta, \phi\) identify the position of the planet in the reference frame of the detector, i.e., where \(\theta = 0\) is the direction of the main beam axis.

The functions described in this chapter can be used to analyze how detectors are going to observe point sources in the sky, properly taking into account proper motions of the sources (this applies to Solar System objects, like planets and comets). The library provides the functions get_ecl2det_quaternions(), which has the same syntax as get_pointings() but returns a matrix with shape (N, 4) containing the N quaternions that transform from the Ecliptic reference frame to the detector’s. Thus, this method can be used to estimate how far from the main beam axis a celestial object is, and its orientation with respect to the orientation of the detector.

Here we show a simple example; the first part is identical to the examples shown above (using the same scanning strategy as CORE’s), but here we employ AstroPy to compute the Ecliptic coordinates of Jupiter during the simulation and convert them in the reference frame of the boresight detector using get_ecl2det_quaternions():

import numpy as np
import litebird_sim as lbs
import astropy.time, astropy.units as u
from astropy.coordinates import (
    ICRS,
    get_body_barycentric,
    BarycentricMeanEcliptic,
    solar_system_ephemeris,
)

sim = lbs.Simulation(
    # We use AstroPy times here!
    start_time=astropy.time.Time("2020-01-01T00:00:00"),
    duration_s=60.0,
    description="Simple simulation",
    random_seed=12345,
)

sim.set_scanning_strategy(
    scanning_strategy=lbs.SpinningScanningStrategy(
        spin_sun_angle_rad=np.deg2rad(30),
        spin_rate_hz=0.5 / 60.0,
        precession_rate_hz=1.0 / (4 * u.day).to("s").value,
    )
)

sim.set_instrument(
    lbs.InstrumentInfo(
        name="core",
        spin_boresight_angle_rad=np.deg2rad(65),
    ),
)

det = lbs.DetectorInfo(name="foo", sampling_rate_hz=10)
obs, = sim.create_observations(detectors=[det])

#################################################################
# Here begins the juicy part

solar_system_ephemeris.set("builtin")

# The variable "icrs_pos" contains the x,y,z coordinates of Jupiter
# in the ICRS reference frame for each sample time in the
# observation.
icrs_pos = get_body_barycentric(
    "jupiter",
    obs.get_times(astropy_times=True),
)
# Convert the ICRS r.f. into the barycentric mean Ecliptic r.f.,
# which is the reference frame used by the LiteBIRD simulation
# framework
ecl_vec = (ICRS(icrs_pos)
    .transform_to(BarycentricMeanEcliptic())
    .cartesian
    .get_xyz()
    .value
)

# The variable ecl_vec is a 3×N matrix containing the vectors.
# We normalize them so that each has length one (using the L_2
# norm, hence ord=2)
ecl_vec /= np.linalg.norm(ecl_vec, axis=0, ord=2)

# Convert the matrix to a N×3 shape
ecl_vec = ecl_vec.transpose()

# Calculate the quaternions that convert the Ecliptic
# reference system into the detector's reference system
quats = lbs.get_ecl2det_quaternions(
    obs,
    sim.spin2ecliptic_quats,
    bore2spin_quat=sim.instrument.bore2spin_quat,
    detector_quats=[det.quat],
)

# Make room for the xyz vectors in the detector's reference frame
det_vec = np.empty_like(ecl_vec)

# Do the rotation!
lbs.all_rotate_vectors(det_vec, quats[0], ecl_vec)

print(det_vec)

[[ 0.57053937  0.07219102 -0.81809124]
 [ 0.57038372  0.06957116 -0.8184267 ]
 [ 0.57023386  0.0669494  -0.81874973]
 ...
 [ 0.99293109 -0.0800506   0.08763421]
 [ 0.99310516 -0.07743726  0.08800916]
 [ 0.99327345 -0.07482179  0.08837171]]

Again, the vectors printed by this script are in the reference frame of the detector, where the vector [0 0 1] indicates the main axis of the detector. We can inspect how close Jupiter moves to the main beam axis during the simulation if we convert the set of (x, y, z) vectors into the angles \(\theta\) (colatitude) and \(\phi\) (longitude), as the colatitude is simply the angular distance from the main beam axis (\(\theta = 0\)):

import healpy
theta, phi = healpy.vec2ang(det_vec)

import matplotlib.pylab as plt

times = obs.get_times()
plt.plot(times - times[0], np.rad2deg(theta))
plt.xlabel("Time [s]")
plt.ylabel("Angular separation [deg]")

We see that Jupiter is ~10° away from the beam axis after ~30 seconds since the start of the simulation.

API reference

class litebird_sim.scanning.ScanningStrategy

Bases: ABC

A class that simulate a scanning strategy

This is an abstract base class; you should probably use SpinningScanningStrategy, unless you are interested in simulating other kinds of scanning strategies. If this is the case, refer to the documentation.

abstract generate_spin2ecl_quaternions(start_time: float | Time, time_span_s: float, delta_time_s: float) → Spin2EclipticQuaternions

Generate the quaternions for spin-axis-to-Ecliptic rotations

This method simulates the scanning strategy of the spacecraft assuming that the mission begins at some time start_time and lasts for time_span_s seconds. The purpose of the function is to compute the orientation of the spacecraft once every delta_time_s seconds for the whole duration of the mission; the orientation is expressed as a quaternion that encodes the rotation from the reference frame of the spacecraft’s spin axis (aligned with the y axis) to the reference frame of the Ecliptic Coordinate System.

The function returns a Spin2EclipticQuaternions object that fully covers the time interval between start_time and start_time + time_span_s: this means that an additional quaternion after the time t_end = start_time + time_span_s might be appended.

Parameters:

• start_time (Union[float, astropy.time.Time]) – start time of the simulation. If it is a floating-point number, it is arbitrary and can usually be set to 0.0; otherwise, it must be a astropy.time.Time object, and in this case a more precise computation of the orientation of the spacecraft is used. Depending on the duration of the simulation, the second case can be a few orders of magnitude slower: it should be used only when the simulation needs to track the position of moving objects (e.g., planets).

• time_span_s (float) – interval of time that needs to be simulated, in seconds. These seconds are added to start_time, and their meaning depends on its type: if start_time is a float, you should consider the duration as a sidereal time, but if it’s a astropy.time.Time time, time_span_s is expressed as a Solar time.

• delta_time_s (float) – for efficiency, quaternions are not sampled at the same sample rate as the scientific data, but at a much lower rate. The default should be good for all the most common cases, but you can tune it with this parameter.

static get_times(start_time: float | Time, delta_time_s: float, num_of_quaternions: int)

Return a vector of equally-separated times

Depending on the type of the parameter start_time (either a float or a astropy.time.Time instance), return a vector of times that mark the instant when a quaternion must be computed by the class.

The class returns a 2-element tuple, containing (1) the time expressed using the same type as start_time (either float or astropy.time.Time), and (2) a vector containing the time measured in seconds. The latter is useful when your scanning strategy depends on the time for the computation of angles and rotation, e.g., if you need to compute \(2\pi\nu t\).

static optimal_num_of_quaternions(time_span_s: float, delta_time_s: float) → int

Return the number of quaternions to compute

Given a time span and a time interval between consecutive quaternions, this static method computes how many quaternions are needed to properly cover the time span.

class litebird_sim.scanning.Spin2EclipticQuaternions(start_time: float | Time, pointing_freq_hz: float, quats)

Bases: object

A matrix of quaternions sampled uniformly over time

This class is used to hold quaternions that represent the transformation from the reference frame of the LiteBIRD spin axis to the Ecliptic reference frame.

The class has the following members:

• start_time is either a floating-point number or an astropy.time.Time object.

• pointing_freq_hz is the sampling frequency of the quaternions, in Hertz

• quats is a NumPy array of shape (N × 4), containing the N (normalized) quaternions

get_detector_quats(detector_quat, bore2spin_quat, time0: float | Time, sampling_rate_hz: float, nsamples: int)

Return detector-to-Ecliptic quaternions

This method combines the spin-axis-to-Ecliptic quaternions in self.quat with two additional rotations (detector_quat, bore2spin_quat), representing the transformation from the reference frame of a detector to the boresight reference frame and the transformation from the boresight to the spin reference frame. The result is a quaternion that directly transforms the reference frame of the detector to Ecliptic coordinates.

Usually, the parameter detector_quat is read from the IMO, and the parameter bore2spin_quat is calculated through the class InstrumentInfo, which has the field bore2spin_quat. If all you have is the angle β (in radians) between the boresight and the spin axis, just pass quat_rotation_y(β) here.

As this kind of quaternion is used to compute the pointings of a detector, which are used in map-making, it applies a «slerp» operation on the quaternion, oversampling them to the sampling frequency of the detector, expressed through the parameter sampling_rate_hz.

The parameters time0 and nsamples specify which is the time interval that needs to be covered by the quaternions computed by this method. The type of the parameter time0 must match that of self.start_time.

nbytes()

Return the number of bytes allocated for the quaternions

class litebird_sim.scanning.SpinningScanningStrategy(spin_sun_angle_rad, precession_rate_hz, spin_rate_hz, start_time=<Time object: scale='tdb' format='iso' value=2027-01-01 00:00:00.000>)

Bases: ScanningStrategy

A class containing the parameters of the sky scanning strategy

This class is used to hold together the parameters that define the nominal scanning strategy of the LiteBIRD spacecraft. It’s a simple scanning strategy that closely matches the ones proposed for other CMB experiments from space like CORE and Pico: a spinning motion of the spacecraft around some axis, composed with a precession motion around the Sun-Earth-spacecraft axis (assuming that the spacecraft flies around the L_2 point of the Sun-Earth system).

The constructor accepts the following parameters:

• spin_sun_angle_rad: angle between the spin axis and the Sun-LiteBIRD direction (floating-point number, in radians)

• precession_rate_hz: the period of the precession rotation (floating-point number, in minutes)

• spin_rate_hz: the number of rotations per minute (RPM) around the spin axis (floating-point number)

• start_time: an astropy.time.Time object representing the start of the observation. It’s currently unused, but it is meant to represent the time when the rotation starts (i.e., the angle ωt is zero).

These fields are available once the object has been initialized.

You can create an instance of this class using the class method from_imo(), which reads the parameters from the IMO.

all_spin_to_ecliptic(result_matrix, sun_earth_angles_rad, time_vector_s)

static from_imo(imo: Imo, url: str | UUID)

Read the definition of the scanning strategy from the IMO

This function returns a SpinningScanningStrategy object containing the set of parameters that define the scanning strategy of the spacecraft, i.e., the way it observes the sky during the nominal mission.

Parameters:

• imo (Imo) – an instance of the Imo class

• url (str or UUID) – a reference to the data file containing the definition of the scanning strategy. It can be either a string like /releases/v1.0/satellite/scanning_parameters/ or a UUID.

Example:

imo = Imo()
sstr = SpinningScanningStrategy.from_imo(
    imo=imo,
    url="/releases/v1.0/satellite/scanning_parameters/",
)
print(sstr)

generate_spin2ecl_quaternions(start_time: float | Time, time_span_s: float, delta_time_s: float) → Spin2EclipticQuaternions

Generate the quaternions for spin-axis-to-Ecliptic rotations

This method simulates the scanning strategy of the spacecraft assuming that the mission begins at some time start_time and lasts for time_span_s seconds. The purpose of the function is to compute the orientation of the spacecraft once every delta_time_s seconds for the whole duration of the mission; the orientation is expressed as a quaternion that encodes the rotation from the reference frame of the spacecraft’s spin axis (aligned with the y axis) to the reference frame of the Ecliptic Coordinate System.

The function returns a Spin2EclipticQuaternions object that fully covers the time interval between start_time and start_time + time_span_s: this means that an additional quaternion after the time t_end = start_time + time_span_s might be appended.

Parameters:

• start_time (Union[float, astropy.time.Time]) – start time of the simulation. If it is a floating-point number, it is arbitrary and can usually be set to 0.0; otherwise, it must be a astropy.time.Time object, and in this case a more precise computation of the orientation of the spacecraft is used. Depending on the duration of the simulation, the second case can be a few orders of magnitude slower: it should be used only when the simulation needs to track the position of moving objects (e.g., planets).

• time_span_s (float) – interval of time that needs to be simulated, in seconds. These seconds are added to start_time, and their meaning depends on its type: if start_time is a float, you should consider the duration as a sidereal time, but if it’s a astropy.time.Time time, time_span_s is expressed as a Solar time.

• delta_time_s (float) – for efficiency, quaternions are not sampled at the same sample rate as the scientific data, but at a much lower rate. The default should be good for all the most common cases, but you can tune it with this parameter.

litebird_sim.scanning.all_compute_pointing_and_orientation(result_matrix, quat_matrix)

Repeatedly apply compute_pointing_and_orientation()

Prototype:

all_compute_pointing_and_orientation(
    result_matrix: numpy.array[D, N, 3],
    quat_matrix: numpy.array[N, D, 4],
)

Assuming that result_matrix is a (D, N, 3) matrix and quat_matrix a (N, D, 4) matrix, iterate over all the N samples and D detectors and apply compute_pointing_and_orientation() to every item.

litebird_sim.scanning.all_spin_to_ecliptic(result_matrix, sun_earth_angles_rad, spin_sun_angle_rad, precession_rate_hz, spin_rate_hz, time_vector_s)

Apply spin_to_ecliptic() to each row of a matrix

Prototype:

all_spin_to_ecliptic(
    result_matrix: numpy.array[N, 4],
    sun_earth_angle_rad: float,
    spin_sun_angle_rad: float,
    precession_rate_hz: float,
    spin_rate_hz: float,
    time_vector_s: numpy.array[N],
)

This function extends spin_to_ecliptic() to work with the vector of times time_vector_s; all the other parameters must still be float as in spin_to_ecliptic; the variable result_matrix must be a matrix of shape (len(time_vector_s), 4).

litebird_sim.scanning.calculate_sun_earth_angles_rad(time_vector)

Compute the angle between the x axis and the Earth

This function computes the angle on the plane of the Ecliptic (assuming to be the xy plane) between the Sun-Earth direction and the x axis. Depending on the type of the parameter time_vector, the result is computed differently:

• If time_vector is a astropy.time.Time object, the angle is computed using the Barycentric Mean Ecliptic reference frame and the Ephemerides tables provided by AstroPy (slow but accurate)

• Otherwise, time_vector is assumed to be a NumPy array of floats, and a simple circular motion with constant angular velocity is assumed. The angular velocity is YEARLY_OMEGA_SPIN_HZ, which is equal to \(2π/T\), with T being the average duration of one year in seconds, and it is assumed that at time t = 0 the angle is zero.

litebird_sim.scanning.compute_pointing_and_orientation(result, quaternion)

Store in “result” the pointing direction and polarization angle.

Prototype:

compute_pointing_and_orientation(
    result: numpy.array[3],
    quaternion: numpy.array[4],
)

The function assumes that quaternion encodes a rotation which transforms the z axis into the direction of a beam in the sky, i.e., it assumes that the beam points towards z in its own reference frame and that quaternion transforms the reference frame to celestial coordinates.

The variable result is used to save the result of the computation, and it should be a 3-element NumPy array. On exit, its values will be:

• result[0]: the colatitude of the sky direction, in radians

• result[1]: the longitude of the sky direction, in radians

• result[2]: the orientation angle (assuming that in the beam reference frame points towards x), measured with respect to the North and East directions in the celestial sphere

This function does not support broadcasting; use all_compute_pointing_and_orientation() if you need to transform several quaternions at once.

Example:

import numpy as np
result = np.empty(3)
compute_pointing_and_orientation(result, np.array([
    0.0, np.sqrt(2) / 2, 0.0, np.sqrt(2) / 2,
])

litebird_sim.scanning.get_det2ecl_quaternions(obs, spin2ecliptic_quats: ~litebird_sim.scanning.Spin2EclipticQuaternions, detector_quats, bore2spin_quat, quaternion_buffer=None, dtype=<class 'numpy.float64'>)

Return the detector-to-Ecliptic quaternions

This function returns a (D, N, 4) tensor containing the quaternions that convert a vector in detector’s coordinates into the frame of reference of the Ecliptic. The number of quaternions is equal to the number of samples hold in this observation, obs.n_samples. Given that the z axis in the frame of reference of a detector points along the main beam axis, this means that if you use these quaternions to rotate the vector z = [0, 0, 1], you will end up with the sequence of vectors pointing towards the points in the sky (in Ecliptic coordinates) that are observed by the detector. This is a low-level function; you should usually call the function get_pointings(), which wraps this function to compute both the pointing direction and the polarization angle. See also the method get_ecl2det_quaternions(), which mirrors this one. If you plan to call this function repeatedly, you can save some running time by pre-allocating the buffer used to hold the quaternions with the parameter quaternion_buffer. This must be a NumPy floating-point array whose shape can be computed using get_quaternion_buffer_shape(). If you pass quaternion_buffer, the return value will be a pointer to this buffer.

litebird_sim.scanning.get_ecl2det_quaternions(obs, spin2ecliptic_quats: ~litebird_sim.scanning.Spin2EclipticQuaternions, detector_quats, bore2spin_quat, quaternion_buffer=None, dtype=<class 'numpy.float64'>)

Return the Ecliptic-to-detector quaternions

This function returns a (D, N, 4) matrix containing the N quaternions of all the D detectors that convert a vector in Ecliptic coordinates into the frame of reference of the detector itself. The number of quaternions is equal to the number of samples hold in this observation. This function is useful when you want to simulate how a point source is observed by the detector’s beam: if you know the Ecliptic coordinates of the point sources, you can easily derive the location of the source with respect to the reference frame of the detector’s beam.

litebird_sim.scanning.get_quaternion_buffer_shape(obs, num_of_detectors=None)

Return the shape of the buffer used to hold detector quaternions.

This function can be used to pre-allocate the buffer used by get_det2ecl_quaternions() and get_ecl2det_quaternions() to save the quaternions representing the change of the orientation of the detectors with time.

Here is a typical use:

import numpy as np
import litebird_sim as lbs
obs = lbs.Observation(...)
bufshape = get_quaternion_buffer_shape(obs, n_detectors)
quaternions = np.empty(bufshape, dtype=np.float64)
quats = get_det2ecl_quaternions(
    ...,
    quaternion_buffer=quaternions,
)

litebird_sim.scanning.orientation_angle(theta_rad, phi_rad, ordir)

Compute the orientation of a detector at a given point on the sky

Prototype:

orientation_angle(
    theta_rad: float,
    phi_rad: float,
    ordir: numpy.array[3],
)

This function returns the orientation (in radians) with respect to the North Pole of the celestial sphere for the point at coordinates theta_rad (colatitude, in radians) and phi_rad (longitude, in radians), assuming that ordir is a 3-element NumPy array representing a normalized vector which departs from the point on the celestial sphere and is aligned with the orientation direction.

litebird_sim.scanning.spin_to_ecliptic(result, sun_earth_angle_rad, spin_sun_angle_rad, precession_rate_hz, spin_rate_hz, time_s)

Compute a quaternion with the spin-axis-to-Ecliptic rotation

Prototype:

spin_to_ecliptic(
    result: numpy.array[4],
    sun_earth_angle_rad: float,
    spin_sun_angle_rad: float,
    precession_rate_hz: float,
    spin_rate_hz: float,
    time_s: float,
)

This function computes the (normalized) quaternion that encodes the rotation which transforms the frame of reference of the spacecraft’s spin axis into the Ecliptic frame of reference. The result is saved in the parameter result, which must be a 4-element NumPy array; the order of the elements of the quaternion is (vx, vy, vz, w).

The function computes the quaternion as the following sequence of rotations:

1. A rotation around the z axis by the angle \(2π ν t\), with ν being the parameter spin_rate_hz and t the parameter time_s (this rotation accounts for the rotation of the spacecraft around the spin axis)

2. A rotation around the y axis by the angle \(π/2 - \alpha\), with ɑ being the parameter spin_sun_angle_rad (this accounts for the inclination of the spin axis with respect to the Ecliptic plane)

3. A rotation around the x axis by the angle \(2π ν t\), with ν being the parameter precession_rate_hz and t the parameter time_s (this rotation accounts for the rotation of the spin axis because of the precessional motion)

4. A rotation around the z axis by the angle sun_earth_angle_rad (this accounts for the yearly revolution of the spacecraft around the Sun)

Parameters:

• sun_earth_angle_rad (float) – Angle between the x axis and the Sun-Earth direction on the xy Ecliptic plane (in radians)

• spin_sun_angle_rad (float) – Angle between the spin axis of the spacecraft and the Sun-Earth direction (in radians); this angle is sometimes called ɑ

• precession_rate_hz (float) – The frequency of rotations around the precession axis (in rotations/sec)

• spin_rate_hz (float) – The frequency of rotations around the spin axis (in rotations/sec)

• time_s (float) – the time when to compute the quaternion

litebird_sim.pointings.apply_hwp_to_obs(obs, hwp: HWP, pointing_matrix)

Modify a pointing matrix to consider the effect of a HWP

This function modifies the variable pointing_matrix (a D×N×3 matrix, with D the number of detectors and N the number of samples) so that the orientation angle considers the behavior of the half-wave plate in hwp.

litebird_sim.pointings.get_pointing_buffer_shape(obs: Observation)

Return the shape of the pointing matrix for a given observation.

This function can be used to determine the size to be passed to NumPy methods that allocate a new matrix in memory, such as numpy.empty and numpy.zeros.

litebird_sim.pointings.get_pointings(obs, spin2ecliptic_quats: ~litebird_sim.scanning.Spin2EclipticQuaternions, bore2spin_quat, detector_quats=None, quaternion_buffer=None, pointing_buffer=None, dtype_pointing=<class 'numpy.float32'>, hwp: ~litebird_sim.hwp.HWP | None = None, store_pointings_in_obs=True)

Return the time stream of pointings for the detector

Given a Spin2EclipticQuaternions and a quaternion representing the transformation from the reference frame of a detector to the boresight reference frame, compute a set of pointings for the detector that encompasses the time span covered by observation obs (i.e., starting from obs.start_time and including obs.n_samples pointings). The parameter spin2ecliptic_quats can be easily retrieved by the field spin2ecliptic_quats in a object of Simulation object, once the method Simulation.set_scanning_strategy() is called.

The parameter bore2spin_quat is calculated through the class Instrument, which has the field bore2spin_quat. If all you have is the angle β between the boresight and the spin axis, just pass quat_rotation_y(β) here.

The parameter detector_quats is optional. By default is None, in this case, if you passed an array of DetectorInfo objects to the method Simulation.create_observations() through the parameter detectors, get_pointings will use the detector quaternions from the same DetectorInfo objects. Otherwise it can contain a stack of detector quaternions. For example, it can be:

• The stack of the field quat of an instance of the class

  DetectorInfo

• If all you want to do is a simulation using a boresight

  direction, you can pass the value np.array([[0., 0., 0., 1.]]), which represents the null rotation.

If HWP is not None, this specifies the HWP to use for the computation of proper polarization angles.

Warning: if hwp is not None, the code adds the α angle of the HWP to the orientation angle ψ, which is generally not correct! This is going to be fixed in the next release of the LiteBIRD Simulation Framework.

The return value is a (D x N × 3) tensor: the colatitude (in radians) is stored in column 0 (e.g., result[:, :, 0]), the longitude (ditto) in column 1, and the polarization angle (ditto) in column 2. You can extract the three vectors using the following idiom:

pointings = obs.get_pointings(...)

# Extract the colatitude (theta), longitude (phi), and
# polarization angle (psi) from pointings
theta, phi, psi = [pointings[:, :, i] for i in (0, 1, 2)]

If you plan to call this function repeatedly, you can save some running time by pre-allocating the buffer used to hold the pointings and the quaternions with the parameters pointing_buffer and quaternion_buffer. Both must be a NumPy floating-point array whose shape can be computed using get_quaternion_buffer_shape() and get_pointing_buffer_shape(). If you use these parameters, the return value will be a pointer to the pointing_buffer.

litebird_sim.pointings.get_pointings_for_observations(obs: ~litebird_sim.observations.Observation | ~typing.List[~litebird_sim.observations.Observation], spin2ecliptic_quats: ~litebird_sim.scanning.Spin2EclipticQuaternions, bore2spin_quat, hwp: ~litebird_sim.hwp.HWP | None = None, store_pointings_in_obs=True, dtype_pointing=<class 'numpy.float32'>)

Obtain pointings for a list of observations

This is a wrapper around the get_pointings() function that computes pointing information for a list of observations and returns a list of pointings. If a single observation is passed then a single pointing array is returned, and, practically, this function only calls get_pointings().

class litebird_sim.hwp.HWP

Bases: object

Abstract class that represents a generic HWP

Being an abstract class, you should never instantiate it. It is used to signal the type of parameters to some functions (e.g., get_pointings()).

If you need to use a HWP object, you should better use derived classes like IdealHWP.

add_hwp_angle(pointing_buffer, start_time_s: float, delta_time_s: float) → None

Modify pointings so that they include the effect of the HWP

This method must be redefined in derived classes. The parameter pointing_buffer must be a D×N×3 matrix representing the three angles (colatitude, longitude, orientation) for D detectors and N measurements. The function only alters orientation and returns nothing.

The parameters start_time_s and delta_time_s specify the time of the first sample in pointings and must be floating-point values; this means that you should already have converted any AstroPy time to a plain scalar before calling this method.

Warning: this changes the interpretation of the ψ variable in the pointings! You should better use get_hwp_angle() and keep ψ and the α angle of the HWP separate. This method is going to be deprecated in a future release of the LiteBIRD Simulation Framework.

get_hwp_angle(output_buffer, start_time_s: float, delta_time_s: float) → None

Calculate the rotation angle of the HWP

This method must be redefined in derived classes. The parameter start_time_s specifies the time of the first sample in pointings and must be a floating-point value; this means that you should already have converted any AstroPy time to a plain scalar before calling this method. The parameter delta_time_s is the inverse of the sampling frequency and must be expressed in seconds.

The result will be saved in output_buffer, which must have already been allocated with the appropriate number of samples.

class litebird_sim.hwp.IdealHWP(ang_speed_radpsec: float, start_angle_rad=0.0)

Bases: HWP

A ideal half-wave plate that spins regularly

This class represents a perfect HWP that spins with constant angular velocity. The constructor accepts the angular speed, expressed in rad/sec, and the start angle (in radians). The latter should be referred to the first time sample in the simulation, i.e., the earliest sample simulated in any of the MPI processes used for the simulation.

Given a polarization angle \(\psi\), this class turns it into \(\psi + \psi_\text{hwp,0} + 2 \omega_\text{hwp} t\), where \(\psi_\text{hwp,0}\) is the start angle specified in the constructor and \(\omega_\text{hwp}\) is the angular speed of the HWP.

add_hwp_angle(pointing_buffer, start_time_s: float, delta_time_s: float) → None

Modify pointings so that they include the effect of the HWP

This method must be redefined in derived classes. The parameter pointing_buffer must be a D×N×3 matrix representing the three angles (colatitude, longitude, orientation) for D detectors and N measurements. The function only alters orientation and returns nothing.

The parameters start_time_s and delta_time_s specify the time of the first sample in pointings and must be floating-point values; this means that you should already have converted any AstroPy time to a plain scalar before calling this method.

Warning: this changes the interpretation of the ψ variable in the pointings! You should better use get_hwp_angle() and keep ψ and the α angle of the HWP separate. This method is going to be deprecated in a future release of the LiteBIRD Simulation Framework.

get_hwp_angle(output_buffer, start_time_s: float, delta_time_s: float) → None

Calculate the rotation angle of the HWP

This method must be redefined in derived classes. The parameter start_time_s specifies the time of the first sample in pointings and must be a floating-point value; this means that you should already have converted any AstroPy time to a plain scalar before calling this method. The parameter delta_time_s is the inverse of the sampling frequency and must be expressed in seconds.

The result will be saved in output_buffer, which must have already been allocated with the appropriate number of samples.

litebird_sim.quaternions.all_rotate_vectors(result_matrix, quat_matrix, vec_matrix)

Rotate a set of vectors using quaternions

Prototype:

all_rotate_vectors(
    result_matrix: numpy.array[N, 3],
    quat_matrix: numpy.array[N, 4],
    vec_matrix: numpy.array[N, 3],
)

Assuming that result_matrix and vec_matrix are two NumPy arrays with shape (N, 3) and quat_matrix with shape (N, 4), apply rotate_vector() to each row.

litebird_sim.quaternions.all_rotate_x_vectors(result_matrix, quat_matrix)

Rotate the vector [1, 0, 0] using quaternions

Prototype:

all_rotate_x_vectors(
    result_matrix: numpy.array[N, 3],
    quat_matrix: numpy.array[N, 4],
)

Assuming that result_matrix is a NumPy array with shape (N, 3) and quat_matrix with shape (N, 4), apply rotate_x_vector() to each row.

litebird_sim.quaternions.all_rotate_y_vectors(result_matrix, quat_matrix)

Rotate the vector [0, 1, 0] using quaternions

Prototype:

all_rotate_y_vectors(
    result_matrix: numpy.array[N, 3],
    quat_matrix: numpy.array[N, 4],
)

Assuming that result_matrix is a NumPy array with shape (N, 3) and quat_matrix with shape (N, 4), apply rotate_y_vector() to each row.

litebird_sim.quaternions.all_rotate_z_vectors(result_matrix, quat_matrix)

Rotate the vector [0, 0, 1] using quaternions

Prototype:

all_rotate_z_vectors(
    result_matrix: numpy.array[N, 3],
    quat_matrix: numpy.array[N, 4],
)

Assuming that result_matrix is a NumPy array with shape (N, 3) and quat_matrix with shape (N, 4), apply rotate_z_vector() to each row.

litebird_sim.quaternions.quat_left_multiply(result, other_v1, other_v2, other_v3, other_w)

Perform a multiplication between two quaternions

Prototype:

quat_left_multiply(
    result: numpy.array[3],
    other_v1: float,
    other_v2: float,
    other_v3: float,
    other_w: float,
)

This function implements the computation \(r = q \times r\); see also quat_right_multiply() for the computation \(r = r\times q\).

It’s easy to use NumPy quaternions for q as well:

import numpy as np
r = np.array([1.0, 2.0, 3.0, 4.0])
q = np.array([0.1, 0.2, 0.3, 0.4])
quat_right_multiply(r, *q)  # Unpack "q"
print("Result:", r)

litebird_sim.quaternions.quat_right_multiply(result, other_v1, other_v2, other_v3, other_w)

Perform a multiplication between two quaternions

Prototype:

quat_right_multiply(
    result: numpy.array[3],
    other_v1: float,
    other_v2: float,
    other_v3: float,
    other_w: float,
)

This function implements the computation \(r = r imes q\), where r is the parameter result (a 3-element NumPy array) and q is the set of parameters other_v1, other_v2, other_v3, other_w. The reason why the elements of quaternion q are passed one by one is efficiency: in this way, the caller does not have to allocate a numpy.array for simple quaternions (like the ones returned by quat_rotation_x(), quat_rotation_y(), quat_rotation_z()).

It’s easy to use NumPy quaternions for q as well:

import numpy as np
r = np.array([1.0, 2.0, 3.0, 4.0])
q = np.array([0.1, 0.2, 0.3, 0.4])
quat_right_multiply(r, *q)  # Unpack "q"
print("Result:", r)

See also quat_left_multiply() for the computation \(r = q imes r\).

litebird_sim.quaternions.quat_rotation_x(theta_rad)

Return a quaternion representing a rotation around the x axis

Prototype:

quat_rotation_x(theta_rad: float) -> Tuple[float, float, float, float]

The angle theta_rad must be expressed in radians. The return value is the quaternion, using the order (v_x, v_y, v_z, w); it is returned as a 4-element tuple.

The fact that the result is a tuple instead of a NumPy array is because of speed: it helps in preventing unnecessary allocations in performance-critical code.

See also quat_rotation_y() and quat_rotation_z()

litebird_sim.quaternions.quat_rotation_y(theta_rad)

Return a quaternion representing a rotation around the y axis

Prototype:

quat_rotation_y(theta_rad: float) -> Tuple[float, float, float, float]

See also quat_rotation_x() and quat_rotation_z()

litebird_sim.quaternions.quat_rotation_z(theta_rad)

Return a quaternion representing a rotation around the y axis

Prototype:

quat_rotation_z(theta_rad: float) -> Tuple[float, float, float, float]

See also quat_rotation_x() and quat_rotation_y()

litebird_sim.quaternions.rotate_vector(result, vx, vy, vz, w, vect)

Rotate a vector using a quaternion

Prototype:

rotate_vector(
    result: numpy.array[3],
    vx: float,
    vy: float,
    vz: float,
    w: float,
    vect: numpy.array[3],
)

Applies a rotation, encoded through the quaternion vx, vy, vz, vw, to the vector vect (a 3-element NumPy array), storing the result in result (again a 3-element array).

Note: do not pass the same variable to vect and result!

The formula to rotate a vector v by a quaternion (q_v, w) is the following: \(v' = v + 2q_v ⨯ (q_v ⨯ v + w v)\), where q_v is the vector (vx, vy, vz).

litebird_sim.quaternions.rotate_x_vector(result, vx, vy, vz, w)

Rotate the x vector using the quaternion (vx, vy, vz, w)

Prototype:

rotate_x_vector(
    result: numpy.array[3],
    vx: float,
    vy: float,
    vz: float,
    w: float,
)

This function is equivalent to rotate_vector(result, vx, vy, vz, w, [1, 0, 0]), but it’s faster.

litebird_sim.quaternions.rotate_y_vector(result, vx, vy, vz, w)

Rotate the x vector using the quaternion (vx, vy, vz, w)

Prototype:

rotate_y_vector(
    result: numpy.array[3],
    vx: float,
    vy: float,
    vz: float,
    w: float,
)

This function is equivalent to rotate_vector(result, vx, vy, vz, w, [0, 1, 0]), but it’s faster.

litebird_sim.quaternions.rotate_z_vector(result, vx, vy, vz, w)

Rotate the x vector using the quaternion (vx, vy, vz, w)

Prototype:

rotate_z_vector(
    result: numpy.array[3],
    vx: float,
    vy: float,
    vz: float,
    w: float,
)

This function is equivalent to rotate_vector(result, vx, vy, vz, w, [0, 0, 1]), but it’s faster.