Pointing systematics HWP class

CHANGELOG.md

@@ -1,4 +1,5 @@
 # HEAD
+-   **Breaking change**: `PointingSys()` now requires `Observation` as an argument. And several functions to add pointing systematics are merged into `left_multiply_syst_quats()`.
 
 -   **Breaking change**: Redefinition (ωt instead of 2ωt) of the hwp angle returned by `get_pointings()`. Change in the pointing returned by `Observation.get_pointings()`, now behaving as it was before this [commit](https://github.com/litebird/litebird_sim/pull/319/commits/b3bc3bb2049c152cc183d6cfc68f4598f5b93ec0). Documentation updated accordingly. [#340](https://github.com/litebird/litebird_sim/pull/340)
 

README.md

@@ -119,13 +119,16 @@ We can visualize detectors in the focal plane by:
 python -m litebird_sim.plot_fp
 ```
 This software loads the IMo which is installed in the machine you are using. 
+
 As the conversation unfolds, an interactive Matplotlib window will appear. 
 Detectors corresponding to the specified channels are represented as blue dots. 
+
 Clicking on a dot reveals the `DetectorInfo` for that detector in real time, highlighted with a red star. 
+
 Additionally, if you agree during the conversation to generate a detector file,
 a list of starred detectors will be saved into a text file at the designated location after you closed the plot.
 
-![plot_fp_usage](https://private-user-images.githubusercontent.com/83496454/388083441-93876c66-8b61-40c0-b23f-79a6998e5e33.gif?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MzIxMDQ2MDcsIm5iZiI6MTczMjEwNDMwNywicGF0aCI6Ii84MzQ5NjQ1NC8zODgwODM0NDEtOTM4NzZjNjYtOGI2MS00MGMwLWIyM2YtNzlhNjk5OGU1ZTMzLmdpZj9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDExMjAlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQxMTIwVDEyMDUwN1omWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPTZlOTk2M2M3NDgyZjUyY2M4NjFiOGE3ZjliM2RhMWU1YWU2NjFkNmM1ZGQxYzAzM2Y3ZDVmOGI3ZTdhMTAwNjEmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0In0.MoAArSYcXhQgVYaer_O72XHIamO_Ex6B5ITuyVgB8_g)
+![plot_fp_usage](https://private-user-images.githubusercontent.com/83496454/388083441-93876c66-8b61-40c0-b23f-79a6998e5e33.gif?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MzIzNjAzOTAsIm5iZiI6MTczMjM2MDA5MCwicGF0aCI6Ii84MzQ5NjQ1NC8zODgwODM0NDEtOTM4NzZjNjYtOGI2MS00MGMwLWIyM2YtNzlhNjk5OGU1ZTMzLmdpZj9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDExMjMlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQxMTIzVDExMDgxMFomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPWJkMmM3NTYzMzc4MTU2NGU3ODE2YWMzZTFlZDVjZmEzODJhNWMwMjhhMDY4MzUxMWU3OGYzZmZiYWFiZDcwNDQmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0In0.zG2QVKA9VeL5j79jTjiPNQ9kH8uyxVSeCzd3GBj3QVk)
 
 ## Roadmap
 

docs/source/part4.rst

@@ -7,3 +7,4 @@ Systematic effects
    gaindrifts.rst
    non_linearity.rst
    hwp_sys.rst
+   pointing_sys.rst

docs/source/plot_fp.rst

@@ -1,7 +1,7 @@
 .. _plot_fp:
 
 Focal plane visualization
-===========
+=========================
 
 We can visualize detectors in the focal plane by:
 

docs/source/pointing_sys.rst

@@ -0,0 +1,243 @@
+.. _pointing_sys:
+
+Pointing systematics
+====================
+The pointing systematics causes some disturbances in the pointing direction of the detectors. Due to the systematics, we will have a signal from the sky where we are not pointing, and the obtained TOD will be perturbed by the systematics.
+
+These TODs with systematics will be used for map-making, though the pointing matrix in the map-maker will be created with the expected pointings, i.e., pointings without systematics.
+
+In order to simulate the pointing systematics, we multiply the systematic quaternions, which describe perturbations, by the expected quaternions.
+For example, to add a 5-degree offset around the :math:`x`-axis to the positional quaternion :math:`q_z=(0,0,1,0)` indicating the :math:`z`-axis, the following operation is available:
+
+.. code-block:: python
+
+    import litebird_sim as lbs
+    import numpy as np
+
+    def print_quaternion(q):
+        print("{:.3f} {:.3f} {:.3f} {:.3f}".format(*q))
+
+    det = lbs.DetectorInfo(
+        name="000_000_003_QA_040_T",
+        sampling_rate_hz=1.0,
+    )
+
+    # Positional quaternion indicating z-axis
+    q_z = lbs.RotQuaternion(
+        quats=np.array([0.0, 0.0, 1.0, 0.0]),
+    )
+
+    # Systematic quaternion indicating 5-degree rotation around x-axis
+    q_sys = lbs.RotQuaternion(
+        quats=np.array(
+            lbs.quat_rotation_x(np.deg2rad(5.0))
+        ),
+    )
+
+    # The systematic quaternions are multiplied by the expected quaternions in-place.
+    lbs.left_multiply_syst_quats(
+        q_z,
+        q_sys,
+        det,
+        start_time=0.0,
+        sampling_rate_hz=det.sampling_rate_hz,
+    )
+
+    print("Rotation by 5 deg. around x-axis:")
+    print_quaternion(q_z.quats[0])
+
+.. testoutput::
+
+    Rotation by 5 deg. around x-axis:
+    -0.000 -0.044 0.999 -0.000
+
+The function :func:`.left_multiply_syst_quats` multiplies the systematic quaternions by the expected quaternions in-place.
+
+By using this function, the :class:`.PointingSys` class is implemented to inject the systematic quaternions into the quaternions in a specific coordinate.
+
+:class:`.PointingSys` class
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+The :class:`.PointingSys` class is used to inject the systematic quaternions into the quaternions in a specific coordinate.
+
+**Note that every coordinate is defined by a left-handed coordinate system, and projected pointings by HEALPix are defined as a view from the spacecraft into the sky.**
+
+The supported coordinates are:
+
+* :class:`.FocalplaneCoord`: The coordinate of the focal plane with the :math:`z`-axis along the boresight.
+    * :meth:`.FocalplaneCoord.add_offset`: Add the offset to the detectors in the focal plane.
+    * :meth:`.FocalplaneCoord.add_disturb`: Add the time-dependent disturbance to the detectors in the focal plane.
+* :class:`.SpacecraftCoord`: The coordinate of the spacecraft with the :math:`z`-axis along its spin axis.
+    * :meth:`.SpacecraftCoord.add_offset`: Add the offset to the entire spacecraft.
+    * :meth:`.SpacecraftCoord.add_disturb`: Add the time-dependent disturbance to the entire spacecraft.
+* :class:`.HWPCoord`: Same coordinates as the focal plane but treats the pointing systematics due to the HWP rotation.
+    * :meth:`.HWPCoord.add_hwp_rot_disturb`: Add the rotational disturbance synchronized with the HWP rotation.
+
+These classes are accessible through the :attr:`.PointingSys.focalplane`, :attr:`.PointingSys.spacecraft`, and :attr:`.PointingSys.hwp` attributes.
+
+How to inject the pointing systematics
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+It is possible to operate :class:`.FocalplaneCoord`, :class:`.SpacecraftCoord`, and :class:`.HWPCoord` directly, but it is recommended to use the :class:`.PointingSys` class to inject the pointing systematics. The following is an example of injecting vibrations into the focal plane.
+
+.. code-block:: python
+
+    import numpy as np
+    import litebird_sim as lbs
+    from pathlib import Path
+    import tomlkit
+    import healpy as hp
+
+    def get_hitmap(nside, pointings):
+        npix = hp.nside2npix(nside)
+        ipix = hp.ang2pix(nside, pointings[:, :, 0], pointings[:, :, 1])
+        hitmap, _ = np.histogram(ipix, bins=np.arange(npix + 1))
+        return hitmap
+
+    start_time = 0.0
+    sampling_hz = 1.0
+    telescope = "LFT"
+    dets = []
+
+    # Load the mock focal plane configuration
+    path_of_toml = (
+        Path(lbs.__file__).resolve().parent.parent
+        / "test"
+        / "pointing_sys_reference"
+        / "mock_focalplane.toml"
+    )
+    with open(path_of_toml, "r", encoding="utf-8") as toml_file:
+        toml_data = tomlkit.parse(toml_file.read())
+        for i in range(len(toml_data[telescope])):
+            dets.append(lbs.DetectorInfo.from_dict(toml_data[telescope][f"det_{i:03}"]))
+
+    # Create a simulation object
+    sim = lbs.Simulation(
+        base_path="pntsys_example",
+        start_time=start_time,
+        duration_s=1000.0,
+        random_seed=12345,
+    )
+
+    # Set the scanning strategy
+    sim.set_scanning_strategy(
+        scanning_strategy=lbs.SpinningScanningStrategy(
+            spin_sun_angle_rad=np.deg2rad(45.0),
+            spin_rate_hz=0.05 / 60.0,
+            precession_rate_hz=1.0 / (3.2 * 60 * 60),
+        ),
+        delta_time_s=1.0 / sampling_hz,
+    )
+
+    # Set the instrument
+    sim.set_instrument(
+        lbs.InstrumentInfo(
+            name="mock_LiteBIRD",
+            spin_boresight_angle_rad=np.deg2rad(50.0),
+        ),
+    )
+
+    # Set the HWP
+    sim.set_hwp(lbs.IdealHWP(sim.instrument.hwp_rpm * 2 * np.pi / 60))
+
+    # Create the observations
+    (obs,) = sim.create_observations(
+        detectors=dets,
+        split_list_over_processes=False,
+    )
+
+    # Initialize the pointing systematics object
+    pntsys = lbs.PointingSys(sim, obs, dets)
+    nquats = obs.n_samples + 1
+    axes = ["x", "y"]
+    noise_sigma_deg = 1.0
+
+    for ax in axes:
+        # Prepare the noise to add vibration to the focal plane
+        noise_rad_array = np.zeros(nquats)
+        lbs.add_white_noise(
+            noise_rad_array, sigma=np.deg2rad(noise_sigma_deg), random=sim.random
+        )
+        # Add the vibration to the pointings
+        pntsys.focalplane.add_disturb(noise_rad_array, ax)
+
+    lbs.prepare_pointings(
+        obs,
+        sim.instrument,
+        sim.spin2ecliptic_quats,
+        sim.hwp,
+    )
+
+    pointings, hwp_angle = obs.get_pointings("all")
+    nside = 64
+    hitmap = get_hitmap(nside, pointings)
+    hp.mollview(hitmap, title="Hit-map with focal plane vibration")
+
+.. image:: images/hitmap_with_vibration.png
+
+In this code snippet, the pointing systematics are injected into the focal plane by adding white noise with a standard deviation of 1 degree to the focal plane. The hit-map is created with the injected pointing systematics.
+
+Tips for MPI distribution
+~~~~~~~~~~~~~~~~~~~~~~~~~
+In the case we consider time-dependent pointing systematics, we may have to increase the number of quaternions by changing ``delta_time_s`` in the :meth:`.Simulation.set_scanning_strategy` method. For example, it is changed by considering the nominal sampling rate of LiteBIRD, i.e., 19 Hz:
+
+.. code-block:: python
+
+    sim.set_scanning_strategy(
+        scanning_strategy=lbs.SpinningScanningStrategy(
+            spin_sun_angle_rad=np.deg2rad(45.0),
+            spin_rate_hz=0.05 / 60.0,
+            precession_rate_hz=1.0 / (3.2 * 60 * 60),
+        ),
+        delta_time_s=1.0 / 19.0,
+    )
+
+However, after this execution, all of :attr:`.Simulation.spin2ecliptic_quats` from start to end in the observation are calculated and stored in memory. And this is not distributed by MPI, so that even when we use MPI, every process will calculate :attr:`.Simulation.spin2ecliptic_quats` from start to end of the observation. This causes a memory issue.
+
+To avoid this issue, a simple trick to distribute the quaternion calculations is:
+
+.. code-block:: python
+
+    # After preparing:
+    # `Simulation` with MPI_COMM_WORLD,
+    # list of `DetectorInfo`,
+    # `Imo` object...
+
+    # Create observation
+    (obs,) = sim.create_observations(
+        detectors=dets,
+        n_blocks_det=1,
+        n_blocks_time=size,
+        split_list_over_processes=False
+    )
+
+    pntsys = lbs.PointingSys(sim, obs, dets)
+
+    # Define scanning strategy parameters from IMO
+    scanning_strategy = lbs.SpinningScanningStrategy.from_imo(
+        url= f"/releases/vPTEP/satellite/scanning_parameters/",
+        imo=imo,
+    )
+
+    # Calculate quaternions for each process
+    sim.spin2ecliptic_quats = scanning_strategy.generate_spin2ecl_quaternions(
+        start_time=obs.start_time,
+        time_span_s=obs.n_samples/obs.sampling_rate_hz,
+        delta_time_s=1.0/19.0,
+    )
+
+    # After injecting the pointing systematics...
+
+    lbs.prepare_pointings(
+        obs,
+        sim.instrument,
+        sim.spin2ecliptic_quats,
+        sim.hwp,
+    )
+
+Especially, this trick is needed when the time-dependent pointing systematics' spectrum has higher frequency than the quaternion's sampling rate given by ``delta_time_s``.
+
+API Reference
+-------------
+.. automodule:: litebird_sim.pointing_sys
+    :members:
+    :undoc-members:
+    :show-inheritance:

litebird_sim/__init__.py

@@ -89,14 +89,12 @@
     prepare_pointings,
     precompute_pointings,
 )
-from .pointing_sys.pointing_sys import (
+from .pointing_sys import (
     get_detector_orientation,
-    left_multiply_offset2det,
-    left_multiply_disturb2det,
-    left_multiply_offset2quat,
-    left_multiply_disturb2quat,
+    left_multiply_syst_quats,
     FocalplaneCoord,
     SpacecraftCoord,
+    HWPCoord,
     PointingSys,
 )
 from .profiler import TimeProfiler, profile_list_to_speedscope
@@ -319,12 +317,10 @@ def destripe_with_toast2(*args, **kwargs):
     "apply_gaindrift_to_observations",
     # pointing_sys.py
     "get_detector_orientation",
-    "left_multiply_offset2det",
-    "left_multiply_disturb2det",
-    "left_multiply_offset2quat",
-    "left_multiply_disturb2quat",
+    "left_multiply_syst_quats",
     "FocalplaneCoord",
     "SpacecraftCoord",
+    "HWPCoord",
     "PointingSys",
     # hwp_diff_emiss.py
     "add_2f",