From 2ac98368c518e176bffa258c35fc3dbd9ba36758 Mon Sep 17 00:00:00 2001
From: Maurizio Tomasi <ziotom78@gmail.com>
Date: Fri, 8 Mar 2024 11:12:55 +0100
Subject: [PATCH 1/5] =?UTF-8?q?Substitute=20=E2=80=9Cpolarization=20angle?=
 =?UTF-8?q?=E2=80=9D=20with=20=E2=80=9Corientation=E2=80=9D=20where=20appr?=
 =?UTF-8?q?opriate?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 benchmarks/pointing_generation.py             |  8 ++-
 ...irection.svg => orientation-direction.svg} |  0
 docs/source/scanning.rst                      | 24 +++----
 docs/source/singularity_shell.cast            |  4 +-
 litebird_sim/__init__.py                      |  8 +--
 litebird_sim/hwp.py                           |  4 +-
 litebird_sim/pointings.py                     |  6 +-
 litebird_sim/scanning.py                      | 62 +++++++++----------
 litebird_sim/toast_destriper/__init__.py      | 14 +++--
 test/test_scanning.py                         | 59 +++++++++---------
 10 files changed, 100 insertions(+), 89 deletions(-)
 rename docs/source/images/{polarization-direction.svg => orientation-direction.svg} (100%)

diff --git a/benchmarks/pointing_generation.py b/benchmarks/pointing_generation.py
index fc2d7fe8..b4a34696 100644
--- a/benchmarks/pointing_generation.py
+++ b/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
+    ),
 )
diff --git a/docs/source/images/polarization-direction.svg b/docs/source/images/orientation-direction.svg
similarity index 100%
rename from docs/source/images/polarization-direction.svg
rename to docs/source/images/orientation-direction.svg
diff --git a/docs/source/scanning.rst b/docs/source/scanning.rst
index 3be5e1f0..a4fffe53 100644
--- a/docs/source/scanning.rst
+++ b/docs/source/scanning.rst
@@ -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,7 +351,7 @@ 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
@@ -357,22 +359,22 @@ angle are computed as follows:
   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
diff --git a/docs/source/singularity_shell.cast b/docs/source/singularity_shell.cast
index 53701bd6..e1aba8d1 100644
--- a/docs/source/singularity_shell.cast
+++ b/docs/source/singularity_shell.cast
@@ -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 "]
diff --git a/litebird_sim/__init__.py b/litebird_sim/__init__.py
index d58cd802..09c5e99b 100644
--- a/litebird_sim/__init__.py
+++ b/litebird_sim/__init__.py
@@ -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,
@@ -215,8 +215,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",
diff --git a/litebird_sim/hwp.py b/litebird_sim/hwp.py
index 775f4373..e8fe8b3f 100644
--- a/litebird_sim/hwp.py
+++ b/litebird_sim/hwp.py
@@ -22,8 +22,8 @@ def add_hwp_angle(self, pointing_buffer, start_time_s: float, delta_time_s: floa
 
         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
diff --git a/litebird_sim/pointings.py b/litebird_sim/pointings.py
index fd6ca65a..15d0ce5e 100644
--- a/litebird_sim/pointings.py
+++ b/litebird_sim/pointings.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`.
     """
 
@@ -152,7 +152,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)
             ),
diff --git a/litebird_sim/scanning.py b/litebird_sim/scanning.py
index 6e958a98..cb9aa7b6 100644
--- a/litebird_sim/scanning.py
+++ b/litebird_sim/scanning.py
@@ -37,30 +37,30 @@ def _clip_sincos(x):
 
 
 @njit
-def polarization_angle(theta_rad, phi_rad, poldir):
-    """Compute the polarization angle at a given point on the sky
+def orientation_angle(theta_rad, phi_rad, ordir):
+    """Compute the orientation of a detector at a given point on the sky
 
     Prototype::
 
-        polarization_angle(
+        orientation_angle(
             theta_rad: float,
             phi_rad: float,
-            poldir: numpy.array[3],
+            ordir: 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
+    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 `poldir` is a 3-element
+    (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
-    polarization direction.
+    orientation direction.
 
     """
-    # We want here to associate a polarization angle with a specific
-    # direction in the sky P = (θ, ϕ) and a polarization direction,
+    # We want here to associate a orientation with a specific
+    # direction in the sky P = (θ, ϕ) and a orientation direction,
     # which is a vector of length one starting from P. To compute the
-    # polarization angle with respect to a fixed frame on the
+    # orientation angle with respect to a fixed frame on the
     # celestial sphere, we need first to derive the two vectors
     # pointing towards North and East.
     #
@@ -104,28 +104,28 @@ def polarization_angle(theta_rad, phi_rad, poldir):
     #   North = [-cos(θ) * cos(ϕ), -cos(θ) * sin(ϕ), sin(θ)]
     #   East  = [-sin(ϕ), cos(ϕ), 0]
     #
-    # To compute the polarization angle, we're just looking at the dot
-    # product between "poldir" and these two directions. We use
+    # To compute the orientation angle, we're just looking at the dot
+    # product between "ordir" and these two directions. We use
     # `clip_sincos` to prevent problems from values that are slightly
     # outside the allowed range [-1,1] because of numerical roundoff
     # errors.
 
-    cos_psi = _clip_sincos(-np.sin(phi_rad) * poldir[0] + np.cos(phi_rad) * poldir[1])
+    cos_psi = _clip_sincos(-np.sin(phi_rad) * ordir[0] + np.cos(phi_rad) * ordir[1])
     sin_psi = _clip_sincos(
-        (-np.cos(theta_rad) * np.cos(phi_rad) * poldir[0])
-        + (-np.cos(theta_rad) * np.sin(phi_rad) * poldir[1])
-        + (np.sin(theta_rad) * poldir[2])
+        (-np.cos(theta_rad) * np.cos(phi_rad) * ordir[0])
+        + (-np.cos(theta_rad) * np.sin(phi_rad) * ordir[1])
+        + (np.sin(theta_rad) * ordir[2])
     )
     return np.arctan2(sin_psi, cos_psi)
 
 
 @njit
-def compute_pointing_and_polangle(result, quaternion):
+def compute_pointing_and_orientation(result, quaternion):
     """Store in "result" the pointing direction and polarization angle.
 
     Prototype::
 
-        compute_pointing_and_polangle(
+        compute_pointing_and_orientation(
             result: numpy.array[3],
             quaternion: numpy.array[4],
         )
@@ -144,19 +144,19 @@ def compute_pointing_and_polangle(result, quaternion):
 
     - ``result[1]``: the longitude of the sky direction, in radians
 
-    - ``result[2]``: the polarization angle (assuming that in the beam
+    - ``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
-    :func:`all_compute_pointing_and_polangle` if you need to
+    :func:`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_polangle(result, np.array([
+        compute_pointing_and_orientation(result, np.array([
             0.0, np.sqrt(2) / 2, 0.0, np.sqrt(2) / 2,
         ])
 
@@ -176,30 +176,30 @@ def compute_pointing_and_polangle(result, quaternion):
     rotate_x_vector(result, vx, vy, vz, w)
 
     # Compute the polarization angle
-    pol_angle = polarization_angle(
-        theta_rad=theta_pointing, phi_rad=phi_pointing, poldir=result
+    orientation = orientation_angle(
+        theta_rad=theta_pointing, phi_rad=phi_pointing, ordir=result
     )
 
     # Finally, set "result" to the true result of the computation
     result[0] = theta_pointing
     result[1] = phi_pointing
-    result[2] = pol_angle
+    result[2] = orientation
 
 
 @njit(parallel=True)
-def all_compute_pointing_and_polangle(result_matrix, quat_matrix):
-    """Repeatedly apply :func:`compute_pointing_and_polangle`
+def all_compute_pointing_and_orientation(result_matrix, quat_matrix):
+    """Repeatedly apply :func:`compute_pointing_and_orientation`
 
     Prototype::
 
-        all_compute_pointing_and_polangle(
+        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
-    :func:`compute_pointing_and_polangle` to every item.
+    :func:`compute_pointing_and_orientation` to every item.
 
     """
 
@@ -212,7 +212,7 @@ def all_compute_pointing_and_polangle(result_matrix, quat_matrix):
 
     for det_idx in range(n_dets):
         for sample_idx in prange(n_samples):
-            compute_pointing_and_polangle(
+            compute_pointing_and_orientation(
                 result_matrix[det_idx, sample_idx, :],
                 quat_matrix[sample_idx, det_idx, :],
             )
@@ -895,7 +895,7 @@ def _precompile():
     )
 
     result = np.empty((1, 10, 3))
-    all_compute_pointing_and_polangle(
+    all_compute_pointing_and_orientation(
         result_matrix=result,
         quat_matrix=np.random.rand(result.shape[1], result.shape[0], 4),
     )
diff --git a/litebird_sim/toast_destriper/__init__.py b/litebird_sim/toast_destriper/__init__.py
index 9ffd58d6..e8a62a92 100644
--- a/litebird_sim/toast_destriper/__init__.py
+++ b/litebird_sim/toast_destriper/__init__.py
@@ -71,14 +71,14 @@ def __init__(
                 self.keydict[f"signal_{det}"] = np.array(obs_tod[i], dtype=np.float64)
 
             theta_phi = curpnt[:, 0:2]
-            polangle = curpnt[:, 2]
+            orientation = curpnt[:, 2]
 
             if coordinates == CoordinateSystem.Galactic:
-                theta_phi, polangle = rotate_coordinates_e2g(
-                    pointings_ecl=theta_phi, pol_angle_ecl=polangle
+                theta_phi, orientation = rotate_coordinates_e2g(
+                    pointings_ecl=theta_phi, pol_angle_ecl=orientation
                 )
             elif coordinates == CoordinateSystem.Ecliptic:
-                pass  # Do nothing, "theta_phi" and "polangle" are ok
+                pass  # Do nothing, "theta_phi" and "orientation" are ok
             else:
                 assert ValueError(
                     "unable to handle coordinate system {coordinates} in `destripe`"
@@ -88,7 +88,11 @@ def __init__(
 
             if polarization:
                 weights = np.stack(
-                    (np.ones(nsamples), np.cos(2 * polangle), np.sin(2 * polangle))
+                    (
+                        np.ones(nsamples),
+                        np.cos(2 * orientation),
+                        np.sin(2 * orientation),
+                    )
                 ).transpose()
             else:
                 weights = np.ones(nsamples).reshape((-1, 1))
diff --git a/test/test_scanning.py b/test/test_scanning.py
index 94ae3b15..be4bfa2c 100644
--- a/test/test_scanning.py
+++ b/test/test_scanning.py
@@ -8,17 +8,17 @@
 from litebird_sim import IdealHWP
 
 
-def test_compute_pointing_and_polangle():
+def test_compute_pointing_and_orientation():
     quat = np.array(lbs.quat_rotation_y(np.pi / 2))
     result = np.empty(3)
-    lbs.compute_pointing_and_polangle(result, quat)
+    lbs.compute_pointing_and_orientation(result, quat)
     assert np.allclose(result, [np.pi / 2, 0.0, -np.pi / 2])
 
     # We stay along the same pointing, but we're rotating the detector
-    # by 90°, so the polarization angle is the only number that
+    # by 90°, so the orientation angle is the only number that
     # changes
     lbs.quat_left_multiply(quat, *lbs.quat_rotation_x(np.pi / 4))
-    lbs.compute_pointing_and_polangle(result, quat)
+    lbs.compute_pointing_and_orientation(result, quat)
     assert np.allclose(result, [np.pi / 2, 0.0, -np.pi / 4])
 
 
@@ -91,19 +91,19 @@ def test_simulation_pointings_still():
     assert np.allclose(np.arctan2(boresight[1], boresight[0]), 2 * np.pi / 365.25)
 
     # Now redo the calculation using get_pointings
-    pointings_and_polangle = lbs.get_pointings(
+    pointings_and_orientation = lbs.get_pointings(
         obs,
         spin2ecliptic_quats=sim.spin2ecliptic_quats,
         bore2spin_quat=instr.bore2spin_quat,
         detector_quats=np.array([[0.0, 0.0, 0.0, 1.0]]),
     )
 
-    colatitude = pointings_and_polangle[..., 0]
-    longitude = pointings_and_polangle[..., 1]
-    polangle = pointings_and_polangle[..., 2]
+    colatitude = pointings_and_orientation[..., 0]
+    longitude = pointings_and_orientation[..., 1]
+    orientation = pointings_and_orientation[..., 2]
 
     assert np.allclose(colatitude, np.pi / 2), colatitude
-    assert np.allclose(np.abs(polangle), np.pi / 2), polangle
+    assert np.allclose(np.abs(orientation), np.pi / 2), orientation
 
     # The longitude should have changed by a fraction 23 hours /
     # 365.25 days of a complete circle (we have 24 samples, from t = 0
@@ -141,26 +141,30 @@ def test_simulation_two_detectors():
 
     instr = lbs.InstrumentInfo(spin_boresight_angle_rad=0.0)
 
-    pointings_and_polangle = lbs.get_pointings(
+    pointings_and_orientation = lbs.get_pointings(
         obs,
         spin2ecliptic_quats=sim.spin2ecliptic_quats,
         bore2spin_quat=instr.bore2spin_quat,
         detector_quats=quaternions,
     )
 
-    assert pointings_and_polangle.shape == (2, 24, 3)
+    assert pointings_and_orientation.shape == (2, 24, 3)
 
-    assert np.allclose(pointings_and_polangle[0, :, 0], pointings_and_polangle[1, :, 0])
-    assert np.allclose(pointings_and_polangle[0, :, 1], pointings_and_polangle[1, :, 1])
+    assert np.allclose(
+        pointings_and_orientation[0, :, 0], pointings_and_orientation[1, :, 0]
+    )
+    assert np.allclose(
+        pointings_and_orientation[0, :, 1], pointings_and_orientation[1, :, 1]
+    )
 
     # The ψ angle should differ by 45°
     assert np.allclose(
-        np.abs(pointings_and_polangle[0, :, 2] - pointings_and_polangle[1, :, 2]),
+        np.abs(pointings_and_orientation[0, :, 2] - pointings_and_orientation[1, :, 2]),
         np.pi / 2,
     )
 
 
-def test_simulation_pointings_polangle(tmp_path):
+def test_simulation_pointings_orientation(tmp_path):
     sim = lbs.Simulation(
         base_path=tmp_path / "simulation_dir",
         start_time=0.0,
@@ -182,18 +186,17 @@ def test_simulation_pointings_polangle(tmp_path):
 
     instr = lbs.InstrumentInfo(spin_boresight_angle_rad=0.0)
 
-    pointings_and_polangle = lbs.get_pointings(
+    pointings_and_orientation = lbs.get_pointings(
         obs,
         spin2ecliptic_quats=sim.spin2ecliptic_quats,
         bore2spin_quat=instr.bore2spin_quat,
         detector_quats=np.array([[0.0, 0.0, 0.0, 1.0]]),
     )
-    polangle = pointings_and_polangle[..., 2]
+    orientation = pointings_and_orientation[..., 2]
 
-    # Check that the polarization angle scans every value between -π
-    # and +π
-    assert np.allclose(np.max(polangle), np.pi, atol=0.01)
-    assert np.allclose(np.min(polangle), -np.pi, atol=0.01)
+    # Check that the orientation scans every value in [-π, +π]
+    assert np.allclose(np.max(orientation), np.pi, atol=0.01)
+    assert np.allclose(np.min(orientation), -np.pi, atol=0.01)
 
     # Simulate the generation of a report
     sim.flush()
@@ -221,13 +224,13 @@ def test_simulation_pointings_spinning(tmp_path):
 
     instr = lbs.InstrumentInfo(spin_boresight_angle_rad=np.deg2rad(15.0))
 
-    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]]),
         bore2spin_quat=instr.bore2spin_quat,
     )
-    colatitude = pointings_and_polangle[..., 0]
+    colatitude = pointings_and_orientation[..., 0]
 
     reference_spin2ecliptic_file = Path(__file__).parent / "reference_spin2ecl.txt.gz"
     reference = np.loadtxt(reference_spin2ecliptic_file)
@@ -235,7 +238,7 @@ def test_simulation_pointings_spinning(tmp_path):
 
     reference_pointings_file = Path(__file__).parent / "reference_pointings.txt.gz"
     reference = np.loadtxt(reference_pointings_file)
-    assert np.allclose(pointings_and_polangle[0, :, :], reference)
+    assert np.allclose(pointings_and_orientation[0, :, :], reference)
 
     # Check that the colatitude does not depart more than ±15° from
     # the Ecliptic
@@ -266,7 +269,7 @@ def test_simulation_pointings_mjd(tmp_path):
     instr = lbs.InstrumentInfo(spin_boresight_angle_rad=np.deg2rad(20.0))
 
     for idx, obs in enumerate(sim.observations):
-        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]]),
@@ -275,7 +278,7 @@ def test_simulation_pointings_mjd(tmp_path):
 
         filename = Path(__file__).parent / f"reference_obs_pointings{idx:03d}.npy"
         reference = np.load(filename, allow_pickle=False)
-        assert np.allclose(pointings_and_polangle, reference)
+        assert np.allclose(pointings_and_orientation, reference)
 
 
 def test_simulation_pointings_hwp_mjd(tmp_path):
@@ -299,7 +302,7 @@ def test_simulation_pointings_hwp_mjd(tmp_path):
     instr = lbs.InstrumentInfo(spin_boresight_angle_rad=np.deg2rad(20.0))
 
     for idx, obs in enumerate(sim.observations):
-        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]]),
@@ -309,7 +312,7 @@ def test_simulation_pointings_hwp_mjd(tmp_path):
 
         filename = Path(__file__).parent / f"reference_obs_pointings_hwp{idx:03d}.npy"
         reference = np.load(filename, allow_pickle=False)
-        assert np.allclose(pointings_and_polangle, reference)
+        assert np.allclose(pointings_and_orientation, reference)
 
 
 def test_scanning_quaternions(tmp_path):

From aa26336a38e0cd9ac6ddc1192a77d5bf67a6b139 Mon Sep 17 00:00:00 2001
From: Maurizio Tomasi <ziotom78@gmail.com>
Date: Fri, 8 Mar 2024 11:13:14 +0100
Subject: [PATCH 2/5] Fix warnings produced by Sphinx

---
 litebird_sim/hwp_sys/hwp_sys.py | 27 ++++++++++++++++-----------
 1 file changed, 16 insertions(+), 11 deletions(-)

diff --git a/litebird_sim/hwp_sys/hwp_sys.py b/litebird_sim/hwp_sys/hwp_sys.py
index 1d1f12ed..add365ae 100644
--- a/litebird_sim/hwp_sys/hwp_sys.py
+++ b/litebird_sim/hwp_sys/hwp_sys.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:

From 2980e76ce0f11f654975b9e7ce7b993ce2e4ab1e Mon Sep 17 00:00:00 2001
From: Maurizio Tomasi <ziotom78@gmail.com>
Date: Mon, 11 Mar 2024 17:21:20 +0100
Subject: [PATCH 3/5] Add a warning in get_pointings if a HWP is being used.

---
 litebird_sim/pointings.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/litebird_sim/pointings.py b/litebird_sim/pointings.py
index 15d0ce5e..082231eb 100644
--- a/litebird_sim/pointings.py
+++ b/litebird_sim/pointings.py
@@ -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

From 0bff1e52a9ea89caf90b6b267760ff82c64b00e1 Mon Sep 17 00:00:00 2001
From: Maurizio Tomasi <ziotom78@gmail.com>
Date: Mon, 11 Mar 2024 17:22:16 +0100
Subject: [PATCH 4/5] Implement a new API for the HWP

---
 litebird_sim/hwp.py | 72 +++++++++++++++++++++++++++++++++++++++------
 1 file changed, 63 insertions(+), 9 deletions(-)

diff --git a/litebird_sim/hwp.py b/litebird_sim/hwp.py
index e8fe8b3f..ba163f95 100644
--- a/litebird_sim/hwp.py
+++ b/litebird_sim/hwp.py
@@ -16,7 +16,28 @@ 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
 
@@ -29,6 +50,11 @@ def add_hwp_angle(self, pointing_buffer, start_time_s: float, delta_time_s: floa
         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,

From db5ae895605172b4e81a67f46ea701052bf7243d Mon Sep 17 00:00:00 2001
From: Maurizio Tomasi <ziotom78@gmail.com>
Date: Fri, 15 Mar 2024 12:18:05 +0100
Subject: [PATCH 5/5] [skip ci] Update CHANGELOG

---
 CHANGELOG.md | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6aa2e125..4d02eb4a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,8 +1,9 @@
 # HEAD
 
--   Fixing bug in mbs to pass general bandpass to mbs [#271](https://github.com/litebird/litebird_sim/pull/271) 
+-   **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)  
+-   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)