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.

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 boresight directions; however, MFT and HFT share the same telescope, and thus it’s customary to show only one boresight. 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:

_images/litebird-scanning-strategy.svg

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 ([1]), 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="Sample simulation",
)

# We now simulate the motion of the spacecraft over a time span
# of one minute (specified in the `duration_s` parameter above).
sim.generate_spin2ecl_quaternions(
    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,
    )
)
instr = 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)
obs, = sim.create_observations(detectors=[det])
pointings = obs.get_pointings(
    sim.spin2ecliptic_quats,
    detector_quats=[det.quat],
    bore2spin_quat=instr.bore2spin_quat,
)

print("Shape:", pointings.shape)
print("Pointings:")
print(pointings)

Shape: (1, 600, 3)
Pointings:
[[[ 2.18166156  0.         -1.57079633]
  [ 2.18164984 -0.00579129 -1.57633021]
  [ 2.18161186 -0.01158231 -1.58186382]
  ...
  [ 0.08854811 -2.96748668 -1.73767999]
  [ 0.08788054 -3.0208052  -1.68678085]
  [ 0.08745441 -3.07477554 -1.63522983]]]

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 hour in the example above), we call the method Simulation.generate_spin2ecl_quaternions(), 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 polarization 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 method Observation.get_pointings().

  3. The method Observation.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 polarization 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 transformation that need to be carried out:

_images/coordinate-systems.svg

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 detector is sensitive to the polarization along the x axis.

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="Sample simulation",
)

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

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:

_images/quaternion-sampling.svg

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 the 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 instances of the Observation class.

Note

Slerp assumes a rotation with constant angular speed 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 polarization angle are computed as follows:

  • The direction is simply 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 polarization angle is given by the angle between the North direction passing through the vector \(\vec d\) (i.e., along the parallel 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:

_images/polarization-direction.svg

The purpose of the method Observation.get_pointings(), used in the example at the beginning of this chapter, is to call Observation.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 polarization angle.

How the boresight is specified

As LiteBIRD includes three focal planes and two telescopes, the specification of the boresight requires some care. In [1] and [2], 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

We saw in the example above that the pointings for detector det were saved in a matrix of shape N × 3, but the meaning of those numbers are probably unclear. 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.

Thus, the method Observation.get_pointings() returns a N×3 matrix whose columns contain the colatitude (in radians), longitude (ditto), and polarization angle with respect to the North and South poles of the sky. 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)
_images/scanning-strategy-example.png

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:

_images/right-handed-coordinates.svg

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="Sample simulation",
)

sim.generate_spin2ecl_quaternions(
    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 = obs.get_pointings(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!

_images/simple-scanning-strategy.png

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 class Observation provides the method 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 the polarization axis.

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 Observation.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="Sample simulation",
)
sim.generate_spin2ecl_quaternions(
    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,
    )
)
instr = 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 = obs.get_ecl2det_quaternions(
    sim.spin2ecliptic_quats,
    detector_quats=[det.quat],
    bore2spin_quat=instr.bore2spin_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.15605956 -0.96302723  0.21959958]
 [-0.15393029 -0.96231749  0.22416628]
 [-0.15180263 -0.96158138  0.22872956]
 ...
 [-0.26224282  0.96497767 -0.00684089]
 [-0.2643768   0.96435188 -0.01141765]
 [-0.26650937  0.96369966 -0.01599139]]

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]")
_images/jupiter-angular-distance.svg

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

Bibliography

1(1,2)

J. Delabrouille and The CORE collaboration. Exploring cosmic origins with CORE: Survey requirements and mission design. \jcap , 2018(4):014, April 2018. arXiv:1706.04516, doi:10.1088/1475-7516/2018/04/014.

2

Shaul Hanany and The PICO collaboration. PICO: Probe of Inflation and Cosmic Origins. arXiv e-prints, pages arXiv:1902.10541, February 2019. arXiv:1902.10541.

API reference

class litebird_sim.scanning.ScanningStrategy

Bases: abc.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: Union[float, astropy.time.core.Time], time_span_s: float, delta_time_s: float) → litebird_sim.scanning.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: Union[float, astropy.time.core.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: Union[float, astropy.time.core.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: Union[float, astropy.time.core.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: litebird_sim.scanning.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: litebird_sim.imo.imo.Imo, url: Union[str, uuid.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: Union[float, astropy.time.core.Time], time_span_s: float, delta_time_s: float) → litebird_sim.scanning.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_polangle

Repeatedly apply compute_pointing_and_polangle()

Prototype:

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

Assuming that result_matrix is a (N×3) matrix and quat_matrix a (N×4) matrix, iterate over all the N rows and apply compute_pointing_and_polangle() to every row.

litebird_sim.scanning.all_spin_to_ecliptic

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_polangle

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

Prototype:

compute_pointing_and_polangle(
    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 polarization 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_polangle() if you need to transform several quaternions at once.

Example:

import numpy as np
result = np.empty(3)
compute_pointing_and_polangle(result, np.array([
    0.0, np.sqrt(2) / 2, 0.0, np.sqrt(2) / 2,
])
litebird_sim.scanning.polarization_angle

Compute the polarization angle at a given point on the sky

Prototype:

polarization_angle(
    theta_rad: float,
    phi_rad: float,
    poldir: numpy.array[3],
)

This function returns the polarization angle (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 poldir is a 3-element NumPy array representing a normalized vector which departs from the point on the celestial sphere and is aligned with the polarization direction.

litebird_sim.scanning.spin_to_ecliptic

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.quaternions.all_rotate_vectors

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

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

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

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

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

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

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

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

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

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

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

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.