Disambiguate between “polarization angle” and “orientation” (#305)

* Substitute “polarization angle” with “orientation” where appropriate

* Fix warnings produced by Sphinx

* Add a warning in get_pointings if a HWP is being used.

* Implement a new API for the HWP

* [skip ci] Update CHANGELOG

Author: ziotom78

CHANGELOG.md

@@ -1,14 +1,18 @@
 # HEAD
 
--   Add data splits in time and detector space to binned maps [#291](https://github.com/litebird/litebird_sim/pull/291)
+-   **Breaking change**: Disambiguate between “polarization angle” and “orientation” [#305](https://github.com/litebird/litebird_sim/pull/305). A few functions have been renamed as a consequence of this change; however, they are low-level functions that are used internally (`compute_pointing_and_polangle`, `all_compute_pointing_and_polangle`, `polarization_angle`), so external codes should be unaffected by this PR.
+
+-   **Breaking change**: Reworking of the IO, `write_observations` and `read_observations` are now part of the class simulation [#293](https://github.com/litebird/litebird_sim/pull/293)  
 
 -   Mbs optionally returns alms instead of maps [#306](https://github.com/litebird/litebird_sim/pull/306)
 
 -   Include the possibility to pass components to fill_tods,  add_dipole and add_noise [#302](https://github.com/litebird/litebird_sim/issues/302)
 
--   Fixing bug in mbs to pass general bandpass to mbs [#271](https://github.com/litebird/litebird_sim/pull/271) 
+-   Add data splits in time and detector space to binned maps [#291](https://github.com/litebird/litebird_sim/pull/291)
 
--   **Breaking change**: Reworking of the IO, `write_observations` and `read_observations` are now part of the class simulation [#293](https://github.com/litebird/litebird_sim/pull/293)  
+-   Add support for partial multithreading using Numba [#276](https://github.com/litebird/litebird_sim/pull/276)
+
+-   Fixing bug in mbs to pass general bandpass to mbs [#271](https://github.com/litebird/litebird_sim/pull/271) 
 
 -   Support for numpy.float128 made optional, this fixes importing issue on ARM architectures [#286](https://github.com/litebird/litebird_sim/pull/286)
 
@@ -20,8 +24,6 @@
 
 -   New module to simulate HWP systematics [#232](https://github.com/litebird/litebird_sim/pull/232)
 
--   Add support for partial multithreading using Numba [#276](https://github.com/litebird/litebird_sim/pull/276)
-
 # Version 0.11.0
 
 -   **Breaking change**: Change the interface to the binner, implement a new destriper, and make the dependency on TOAST optional [#260](https://github.com/litebird/litebird_sim/pull/260)

benchmarks/pointing_generation.py

@@ -61,7 +61,7 @@
 
 # Compute the pointings by running a "slerp" operation
 start = time.perf_counter_ns()
-pointings_and_polangle = lbs.get_pointings(
+pointings_and_orientation = lbs.get_pointings(
     obs,
     spin2ecliptic_quats=sim.spin2ecliptic_quats,
     detector_quats=np.array([[0.0, 0.0, 0.0, 1.0]]),
@@ -71,7 +71,9 @@
 elapsed_time = (stop - start) * 1.0e-9
 
 print("Elapsed time for get_pointings: {} s".format((stop - start) * 1e-9))
-print("Shape of the pointings: ", pointings_and_polangle.shape)
+print("Shape of the pointings: ", pointings_and_orientation.shape)
 print(
-    "Speed: {:.1e} pointings/s".format(pointings_and_polangle.shape[1] / elapsed_time),
+    "Speed: {:.1e} pointings/s".format(
+        pointings_and_orientation.shape[1] / elapsed_time
+    ),
 )

docs/source/images/polarization-direction.svg renamed to docs/source/images/orientation-direction.svg

@@ -192,7 +192,7 @@ for now just keep in mind the overall shape of the code:
 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
+   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
@@ -201,7 +201,7 @@ for now just keep in mind the overall shape of the code:
 3. The function :func:`.get_pointings` returns a ``(N, 3)``
    matrix, where the first column contains the colatitude
    :math:`\theta`, the second column the longitude :math:`\phi`, and
-   the third column the polarization angle :math:`\psi`, all expressed
+   the third column the orientation angle :math:`\psi`, all expressed
    in radians.
 
 
@@ -215,7 +215,9 @@ 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 detector is sensitive to the polarization along the `x` axis.
+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
@@ -349,30 +351,30 @@ split in several blocks inside the :class:`.Observation` class.
    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
+rate, the direction of the detector on the sky and its orientation
 angle are computed as follows:
 
 - The direction is the vector :math:`\vec d = R \hat e_z`, where
   :math:`R` is the overall rotation from the detector's reference
   frame to the Ecliptic reference frame, and :math:`\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
+- The orientation angle is given by the angle between the North
   direction passing through the vector :math:`\vec d` (i.e., along the
   meridian of :math:`\vec d`) and the vector :math:`\vec p = R \hat
   e_x`, where :math:`R` is the same as above and :math:`\hat e_x = (1,
   0, 0)`, as shown in the following figure (note that :math:`\hat e_x`
   has been drawn twice because the one in the upper side represents
-  the polarization direction in the detector's reference frame):
+  the orientation direction in the detector's reference frame):
 
-  .. image:: images/polarization-direction.svg
+  .. image:: images/orientation-direction.svg
 
 The purpose of the method :meth:`.Simulation.compute_pointings`, used
 in the example at the beginning of this chapter, is to call
 :func:`.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.
+orientation angle.
 
 
 How the boresight is specified
@@ -437,7 +439,7 @@ of the :math:`n_d` detectors kept in the field ``pointings`` of the
 :class:`.Observation`; they are laid out in memory as a :math:`(n_d,
 N, 2)` matrix, where :math:`N` is the number of samples in the
 timeline, and the last dimension holds the colatitude and longitude
-(in radians). The polarization angle is kept in ``Observation.psi``.
+(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
@@ -719,7 +721,7 @@ 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
+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 :class:`.HWP` to the method
 :meth:`.Simulation.set_hwp`. Here is an example::
@@ -796,7 +798,7 @@ provides the functions
 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.
+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

docs/source/scanning.rst

@@ -202,7 +202,7 @@
 [37.248091, "o", "\r\n"]
 [37.248202, "o", "test/test_quaternions.py::test_collective_quick_rotations "]
 [37.426887, "o", "\u001b[32mPASSED\u001b[0m\u001b[32m   [ 66%]\u001b[0m"]
-[37.42769, "o", "\r\ntest/test_scanning.py::test_compute_pointing_and_polangle "]
+[37.42769, "o", "\r\ntest/test_scanning.py::test_compute_pointing_and_orientation "]
 [37.610971, "o", "\u001b[32mPASSED\u001b[0m\u001b[32m   [ 68%]\u001b[0m"]
 [37.611738, "o", "\r\ntest/test_scanning.py::test_spin_to_ecliptic "]
 [37.612417, "o", "\u001b[32mPASSED\u001b[0m\u001b[32m                [ 70%]\u001b[0m"]
@@ -211,7 +211,7 @@
 [37.61575, "o", "\u001b[32m  [ 72%]\u001b[0m"]
 [37.616256, "o", "\r\ntest/test_scanning.py::test_simulation_pointings_still "]
 [38.069883, "o", "\u001b[32mPASSED\u001b[0m\u001b[32m      [ 75%]\u001b[0m"]
-[38.07046, "o", "\r\ntest/test_scanning.py::test_simulation_pointings_polangle "]
+[38.07046, "o", "\r\ntest/test_scanning.py::test_simulation_pointings_orientation "]
 [38.810859, "o", "\u001b[32mPASSED\u001b[0m"]
 [38.810898, "o", "\u001b[32m   [ 77%]\u001b[0m"]
 [38.811539, "o", "\r\ntest/test_scanning.py::test_simulation_pointings_spinning "]

docs/source/singularity_shell.cast

@@ -64,8 +64,8 @@
     all_rotate_z_vectors,
 )
 from .scanning import (
-    compute_pointing_and_polangle,
-    all_compute_pointing_and_polangle,
+    compute_pointing_and_orientation,
+    all_compute_pointing_and_orientation,
     spin_to_ecliptic,
     all_spin_to_ecliptic,
     calculate_sun_earth_angles_rad,
@@ -216,8 +216,8 @@ def destripe_with_toast2(*args, **kwargs):
     "all_rotate_y_vectors",
     "all_rotate_z_vectors",
     # scanning.py
-    "compute_pointing_and_polangle",
-    "all_compute_pointing_and_polangle",
+    "compute_pointing_and_orientation",
+    "all_compute_pointing_and_orientation",
     "spin_to_ecliptic",
     "all_spin_to_ecliptic",
     "calculate_sun_earth_angles_rad",

litebird_sim/__init__.py

@@ -16,19 +16,45 @@ class HWP:
     classes like :class:`.IdealHWP`.
     """
 
-    def add_hwp_angle(self, pointing_buffer, start_time_s: float, delta_time_s: float):
+    def get_hwp_angle(
+        self, 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.
+        """
+        raise NotImplementedError(
+            "You should not use the HWP class in your code, use IdealHWP instead"
+        )
+
+    def add_hwp_angle(
+        self, 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, polangle)`` for D detectors and N measurements.
-        The function only alters ``polangle`` and returns nothing.
+        ``(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 :meth:`.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.
         """
         raise NotImplementedError(
             "You should not use the HWP class in your code, use IdealHWP instead"
@@ -41,18 +67,33 @@ def __str__(self):
 
 
 @njit
+def _get_ideal_hwp_angle(
+    output_buffer, start_time_s, delta_time_s, start_angle_rad, ang_speed_radpsec
+):
+    for sample_idx in range(output_buffer.size):
+        angle = (
+            start_angle_rad
+            + (start_time_s + delta_time_s * sample_idx) * 2 * ang_speed_radpsec
+        ) % (2 * np.pi)
+
+        output_buffer[sample_idx] = angle
+
+
 def _add_ideal_hwp_angle(
     pointing_buffer, start_time_s, delta_time_s, start_angle_rad, ang_speed_radpsec
 ):
     detectors, samples, _ = pointing_buffer.shape
-    for det_idx in range(detectors):
-        for sample_idx in range(samples):
-            angle = (
-                start_angle_rad
-                + (start_time_s + delta_time_s * sample_idx) * 2 * ang_speed_radpsec
-            ) % (2 * np.pi)
+    hwp_angles = np.empty(samples, dtype=pointing_buffer.dtype)
+    _get_ideal_hwp_angle(
+        output_buffer=hwp_angles,
+        start_time_s=start_time_s,
+        delta_time_s=delta_time_s,
+        start_angle_rad=start_angle_rad,
+        ang_speed_radpsec=ang_speed_radpsec,
+    )
 
-            pointing_buffer[det_idx, sample_idx, 2] += angle
+    for det_idx in range(detectors):
+        pointing_buffer[det_idx, :, 2] += hwp_angles
 
 
 class IdealHWP(HWP):
@@ -75,7 +116,20 @@ def __init__(self, ang_speed_radpsec: float, start_angle_rad=0.0):
         self.ang_speed_radpsec = ang_speed_radpsec
         self.start_angle_rad = start_angle_rad
 
-    def add_hwp_angle(self, pointing_buffer, start_time_s: float, delta_time_s: float):
+    def get_hwp_angle(
+        self, output_buffer, start_time_s: float, delta_time_s: float
+    ) -> None:
+        _get_ideal_hwp_angle(
+            output_buffer=output_buffer,
+            start_time_s=start_time_s,
+            delta_time_s=delta_time_s,
+            start_angle_rad=self.start_angle_rad,
+            ang_speed_radpsec=self.ang_speed_radpsec,
+        )
+
+    def add_hwp_angle(
+        self, pointing_buffer, start_time_s: float, delta_time_s: float
+    ) -> None:
         _add_ideal_hwp_angle(
             pointing_buffer=pointing_buffer,
             start_time_s=start_time_s,

litebird_sim/hwp.py

@@ -1128,19 +1128,24 @@ def fill_tod(
     ):
         r"""It fills tod and/or :math:`A^T A` and :math:`A^T d` for the
         "on the fly" map production
+
         Args:
-            obs class:`Observations`: container for tod. If the tod is
-                 not required, obs.tod can be not allocated
-                 i.e. in lbs.Observation allocate_tod=False.
-            pointings (float): pointing for each sample and detector
-                 generated by func:lbs.get_pointings
-                 Optional if already allocated in obs.
+
+        obs (:class:`Observation`): container for tod. If the tod is
+                 not required, you can avoid allocating ``obs.tod``
+                 i.e. in ``lbs.Observation`` use ``allocate_tod=False``.
+
+        pointings (float): pointing for each sample and detector
+                 generated by :func:`lbs.get_pointings`
+                 Optional if already allocated in ``obs``.
                  When generating pointing information, set the variable
-                 `hwp` to None since the hwp rotation angle is added to
-                 the polarization angle within the `fill_tod` function.
-            hwp_radpsec (float): hwp rotation speed in radiants per second
-            save_tod (bool): if True, the `tod` is saved in `obs.tod`, if False,
-                 `tod` gets deleted
+                 ``hwp`` to None since the hwp rotation angle is added to
+                 the orientation angle within the ``fill_tod`` function.
+
+        hwp_radpsec (float): hwp rotation speed in radiants per second
+
+        save_tod (bool): if True, ``tod`` is saved in ``obs.tod``; if False,
+                 ``tod`` gets deleted
         """
 
         if pointings is None:

litebird_sim/hwp_sys/hwp_sys.py

@@ -14,7 +14,7 @@
 
 from .scanning import (
     get_det2ecl_quaternions,
-    all_compute_pointing_and_polangle,
+    all_compute_pointing_and_orientation,
 )
 
 from .coordinates import CoordinateSystem
@@ -25,7 +25,7 @@ def apply_hwp_to_obs(obs, hwp: HWP, pointing_matrix):
 
     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
-    polarization angle considers the behavior of the half-wave plate in
+    orientation angle considers the behavior of the half-wave plate in
     `hwp`.
     """
 
@@ -99,6 +99,11 @@ def get_pointings(
     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
@@ -152,7 +157,7 @@ def get_pointings(
         )
 
         # Compute the pointing direction for each sample
-        all_compute_pointing_and_polangle(
+        all_compute_pointing_and_orientation(
             result_matrix=pointing_buffer[idx, :, :].reshape(
                 (1, pointing_buffer.shape[1], 3)
             ),