Change the convention for the polarization angle

docs/source/timeordered.rst

@@ -92,17 +92,16 @@ shows:
 
 .. testoutput::
 
-    0.00000
-    -0.00075
-    -0.00151
-    -0.00226
-    -0.00301
-    -0.00376
-    -0.00451
-    -0.00526
-    -0.00601
-    -0.00676   
-
+    2.00000
+    2.00075
+    2.00151
+    2.00226
+    2.00301
+    2.00376
+    2.00451
+    2.00526
+    2.00601
+    2.00676
 
 The input maps to scan can be either included in a dictionary with the name of
 the channel or the name of the dectector as keyword (the routines described in 
@@ -341,12 +340,11 @@ transparent:
 
 .. testoutput::
 
-    4.14241e-04
-    5.46700e-05
-    3.03378e-04
-    6.13975e-05
-    4.72613e-05
-
+    3.51023e-04
+    -9.19179e-06
+    2.38899e-04
+    -3.67181e-06
+    -1.83718e-05
 
 API reference
 -------------

litebird_sim/scanning.py

@@ -106,18 +106,23 @@ def polarization_angle(theta_rad, phi_rad, poldir):
     #   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
-    # `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])
-    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])
+    # product between "poldir" and these two directions.
+
+    sin_theta = np.sin(theta_rad)
+    rx = sin_theta * np.cos(phi_rad)
+    ry = sin_theta * np.sin(phi_rad)
+    rz = np.cos(theta_rad)
+
+    # This formula assumes that ψ=0 is aligned along the
+    # "East" vector indicated above
+    by = poldir[0] * ry - poldir[1] * rx
+    bx = (
+        poldir[0] * (-rz * rx)
+        + poldir[1] * (-rz * ry)
+        + poldir[2] * (rx**2 + ry**2)
     )
-    return np.arctan2(sin_psi, cos_psi)
+
+    return np.arctan2(by, bx)
 
 
 @njit

test/test_scanning.py

@@ -12,14 +12,14 @@ def test_compute_pointing_and_polangle():
     quat = np.array(lbs.quat_rotation_y(np.pi / 2))
     result = np.empty(3)
     lbs.compute_pointing_and_polangle(result, quat)
-    assert np.allclose(result, [np.pi / 2, 0.0, -np.pi / 2])
+    assert np.allclose(result, [np.pi / 2, 0.0, np.pi])
 
     # We stay along the same pointing, but we're rotating the detector
     # by 90°, so the polarization 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)
-    assert np.allclose(result, [np.pi / 2, 0.0, -np.pi / 4])
+    assert np.allclose(result, [np.pi / 2, 0.0, -3 * np.pi / 4])
 
 
 def test_spin_to_ecliptic():
@@ -103,7 +103,7 @@ def test_simulation_pointings_still():
     polangle = pointings_and_polangle[..., 2]
 
     assert np.allclose(colatitude, np.pi / 2), colatitude
-    assert np.allclose(np.abs(polangle), np.pi / 2), polangle
+    assert np.allclose(np.abs(polangle), np.pi), polangle
 
     # The longitude should have changed by a fraction 23 hours /
     # 365.25 days of a complete circle (we have 24 samples, from t = 0
@@ -155,7 +155,8 @@ def test_simulation_two_detectors():
 
     # 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_polangle[0, :, 2] - pointings_and_polangle[1, :, 2])
+        % np.pi,
         np.pi / 2,
     )
 
@@ -224,10 +225,12 @@ def test_simulation_pointings_spinning(tmp_path):
     colatitude = pointings_and_polangle[..., 0]
 
     reference_spin2ecliptic_file = Path(__file__).parent / "reference_spin2ecl.txt.gz"
+    # np.savetxt(reference_spin2ecliptic_file, sim.spin2ecliptic_quats.quats)
     reference = np.loadtxt(reference_spin2ecliptic_file)
     assert np.allclose(sim.spin2ecliptic_quats.quats, reference)
 
     reference_pointings_file = Path(__file__).parent / "reference_pointings.txt.gz"
+    # np.savetxt(reference_pointings_file, pointings_and_polangle[0, :, :])
     reference = np.loadtxt(reference_pointings_file)
     assert np.allclose(pointings_and_polangle[0, :, :], reference)
 
@@ -267,6 +270,7 @@ def test_simulation_pointings_mjd(tmp_path):
         )
 
         filename = Path(__file__).parent / f"reference_obs_pointings{idx:03d}.npy"
+        # np.save(filename, pointings_and_polangle)
         reference = np.load(filename, allow_pickle=False)
         assert np.allclose(pointings_and_polangle, reference)
 
@@ -300,6 +304,7 @@ def test_simulation_pointings_hwp_mjd(tmp_path):
         )
 
         filename = Path(__file__).parent / f"reference_obs_pointings_hwp{idx:03d}.npy"
+        # np.save(filename, pointings_and_polangle)
         reference = np.load(filename, allow_pickle=False)
         assert np.allclose(pointings_and_polangle, reference)