Merge pull request #340 from litebird/hwp_angle

Redefinition of hwp_angle

CHANGELOG.md

@@ -1,5 +1,7 @@
 # HEAD
 
+-   **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)
+
 -   Module for including nonlinearity in the simulations [#331](https://github.com/litebird/litebird_sim/pull/331)
 
 -   Improve the documentation of the binner and the destriper [#333](https://github.com/litebird/litebird_sim/pull/333)

docs/source/non_linearity.rst

@@ -160,9 +160,9 @@ If the 2f amplitude is not read from the IMo, one has to specify :math:`A_2` usi
     sim.add_2f(component="tod_2f")
 
 
-Refer to the :ref:`gd-api-reference` for the full list of non-linearity simulation parameters.
+Refer to the :ref:`nl-api-reference` for the full list of non-linearity simulation parameters.
 
-.. _gd-api-reference:
+.. _nl-api-reference:
 
 API reference
 -------------

docs/source/scanning.rst

@@ -168,15 +168,15 @@ similar to what is going to be used for LiteBIRD:
 
 .. testoutput::
 
-  Shape: (1, 600, 3)
+  Shape: (600, 3)
   Pointings:
-  [[[ 2.182 -0.    -1.571]
-    [ 2.182 -0.006 -1.576]
-    [ 2.182 -0.012 -1.582]
-    ...
-    [ 0.089 -2.967 -1.738]
-    [ 0.088 -3.021 -1.687]
-    [ 0.087 -3.075 -1.635]]]
+  [[ 2.182 -0.    -1.571]
+   [ 2.182 -0.006 -1.576]
+   [ 2.182 -0.012 -1.582]
+   ...
+   [ 0.089 -2.967 -1.738]
+   [ 0.088 -3.021 -1.687]
+   [ 0.087 -3.075 -1.635]]
 
 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:
@@ -201,12 +201,27 @@ for now just keep in mind the overall shape of the code:
    the example above, this is done by the function
    :func:`.get_pointings`.
 
-3. The method :meth:`.Observation.get_pointings` returns a ``(D, N, 3)``
-   matrix, where D represents the detector index, N the index of the sample
-   and the three final columns contain the colatitude :math:`\theta`, 
-   the longitude :math:`\phi`, and the orientation angle :math:`\psi`, 
-   all expressed in radians. These angles are expressed in the Ecliptic
-   Coordinate System, where the Equator is aligned with the Ecliptic Plane of
+3. The method :meth:`.Observation.get_pointings` returns an array with
+   either 2 or 3 fields depending on the argument passed:
+
+   - if an integer is passed, this is interpreted as the index of the
+     detector in the observation, and a ``(N, 3)`` matrix is returned 
+     where the first column contains the colatitude :math:`\theta`, 
+     the second column the longitude :math:`\phi`, and the third column
+     the orientation angle :math:`\psi`, all expressed in radians.
+
+   - if a list containing indices is passed, this is interpreted as
+     a list of detectors in the observation for which we want to compute
+     the pointing. It returns a ``(D, N, 3)`` matrix where D represents 
+     the detector index, N the index of the sample and the three final
+     columns are the same described in the first case.
+
+   - if the string "all" is passed then a ``(D, N, 3)`` matrix is returned
+     containig the pointing information for all the detectors in the 
+     observation.
+
+   These angles are expressed in the Ecliptic Coordinate
+   System, where the Equator is aligned with the Ecliptic Plane of
    the Solar System.
 
 
@@ -361,7 +376,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 o]rientation
+rate, the direction of the detector on the sky and its orientation
 angle can be computed via a call to :meth:`.Observation.get_pointings`.
 The calculation works as follows:
 
@@ -679,6 +694,11 @@ few lines of code:
 
 The following code implements our mock scanning strategy::
 
+   import litebird_sim as lbs
+   from litebird_sim import RotQuaternion
+   import astropy
+   from typing import Union
+
    class SimpleScanningStrategy(lbs.ScanningStrategy):
        def generate_spin2ecl_quaternions(
            self,
@@ -733,7 +753,7 @@ The following code implements our mock scanning strategy::
            # "RotQuaternion"
            return lbs.RotQuaternion(
                start_time=start_time,
-               pointing_freq_hz=1.0 / delta_time_s,
+               sampling_rate_hz=1.0 / delta_time_s,
                quats=spin2ecliptic_quats,
            )
 
@@ -770,8 +790,10 @@ computing one quaternion every minute, we compute one quaternion every
        name="foo",
        sampling_rate_hz=1.0 / ((1.0 * u.day).to("s").value),
    )
+
    (obs,) = sim.create_observations(detectors=[det])
-   pointings = lbs.get_pointings(obs, sim.spin2ecliptic_quats, np.array([det.quat]))
+   lbs.prepare_pointings(obs, lbs.InstrumentInfo(), sim.spin2ecliptic_quats)
+   pointings, _ = obs.get_pointings("all")
 
    m = np.zeros(healpy.nside2npix(64))
    pixidx = healpy.ang2pix(64, pointings[0, :, 0], pointings[0, :, 1])
@@ -810,7 +832,7 @@ of a descendant of the class :class:`.HWP` to the method
     )
 
     sim.set_instrument(
-        instr = lbs.InstrumentInfo(
+        instrument = lbs.InstrumentInfo(
             boresight_rotangle_rad=0.0,
             spin_boresight_angle_rad=0.872_664_625_997_164_8,
             spin_rotangle_rad=3.141_592_653_589_793,
@@ -830,9 +852,23 @@ of a descendant of the class :class:`.HWP` to the method
 
     sim.prepare_pointings()
 
+    pointing, hwp_angle = sim.observations[0].get_pointings()
 
 This example uses the :class:`.IdealHWP`, which represents an ideal
 spinning HWP.
+The method :func:`.get_pointings` returns both pointing and hwp angle
+on the fly.
+As shown before a similar syntax allow to precompute the pointing and the 
+hwp angle::
+
+    sim.prepare_pointings()
+    sim.precompute_pointings()
+
+    sim.observations[0].pointing_matrix.shape
+    sim.observations[0].hwp_angle.shape
+
+This fills the fields `pointing_matrix` and `hwp_angle` for all the
+observations in the current simulation. 
 
 
 Observing point sources in the sky

litebird_sim/dipole.py

@@ -225,7 +225,7 @@ def add_dipole(
         if type(pointings) is np.ndarray:
             theta_phi_det = pointings[detector_idx, :, :]
         else:
-            theta_phi_det = pointings(detector_idx)[0][0, :, 0:2]
+            theta_phi_det = pointings(detector_idx)[0][:, 0:2]
 
         add_dipole_for_one_detector(
             tod_det=tod[detector_idx],

litebird_sim/hwp.py

@@ -130,7 +130,7 @@ def _get_ideal_hwp_angle(
     for sample_idx in range(output_buffer.size):
         angle = (
             start_angle_rad
-            + (start_time_s + delta_time_s * sample_idx) * 2 * ang_speed_radpsec
+            + (start_time_s + delta_time_s * sample_idx) * ang_speed_radpsec
         ) % (2 * np.pi)
 
         output_buffer[sample_idx] = angle
@@ -139,7 +139,14 @@ def _get_ideal_hwp_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
+    shape = pointing_buffer.shape
+    if len(shape) == 3:
+        detectors, samples, _ = shape
+    elif len(shape) == 2:
+        detectors, samples = shape
+    else:
+        ValueError("Wrong shape for the pointing_buffer")
+
     hwp_angles = np.empty(samples, dtype=pointing_buffer.dtype)
     _get_ideal_hwp_angle(
         output_buffer=hwp_angles,
@@ -148,9 +155,12 @@ def _add_ideal_hwp_angle(
         start_angle_rad=start_angle_rad,
         ang_speed_radpsec=ang_speed_radpsec,
     )
-
-    for det_idx in range(detectors):
-        pointing_buffer[det_idx, :, 2] += hwp_angles
+    if len(shape) == 3:
+        for det_idx in range(detectors):
+            pointing_buffer[det_idx, :, 2] += 2 * hwp_angles
+    elif len(shape) == 2:
+        for det_idx in range(detectors):
+            pointing_buffer[det_idx, :] += 2 * hwp_angles
 
 
 class IdealHWP(HWP):
@@ -202,7 +212,7 @@ def apply_hwp_to_pointings(
         bore2ecl_quaternions_inout: npt.NDArray,  # Boresight→Ecliptic quaternions at 19 Hz
         hwp_angle_out: npt.NDArray,
     ) -> None:
-        # We do not touch `bore2ecl_quaternions_inout`, as an ideal HWP does not
+        # We do not touch `bore2ecl_quaternions_input`, as an ideal HWP does not
         # alter the (θ, φ) direction of the boresight nor the orientation ψ
         self.get_hwp_angle(
             output_buffer=hwp_angle_out,

litebird_sim/hwp_diff_emiss.py

@@ -20,7 +20,7 @@
 # here we calculate 2f directly.
 @njit
 def compute_2f_for_one_sample(angle_rad, amplitude_k, monopole_k):
-    return amplitude_k * np.cos(angle_rad) + monopole_k
+    return amplitude_k * np.cos(2 * angle_rad) + monopole_k
 
 
 @njit(parallel=True)

litebird_sim/hwp_sys/hwp_sys.py

@@ -1148,7 +1148,7 @@ def fill_tod(
                 If ``observations`` is a list, ``pointings`` must be a list of same length.
 
 
-        hwp_angle (optional): `2ωt`, hwp rotation angles (radians). If ``pointings`` is passed,
+        hwp_angle (optional): `ωt`, hwp rotation angles (radians). If ``pointings`` is passed,
                 ``hwp_angle`` must be passed as well, otherwise both must be ``None``.
                 If not passed, it is computed on the fly (generated by :func:`lbs.get_pointings`
                 per detector).
@@ -1256,7 +1256,6 @@ def fill_tod(
                         cur_point, cur_hwp_angle = cur_obs.get_pointings(
                             detector_idx=idet, pointings_dtype=dtype_pointings
                         )
-                        cur_point = cur_point.reshape(-1, 3)
                     else:
                         cur_point = ptg_list[idx_obs][idet, :, :]
                         cur_hwp_angle = hwp_angle_list[idx_obs]
@@ -1297,7 +1296,7 @@ def fill_tod(
                         mUQ=self.mUQ,
                         mQU=self.mQU,
                         pixel_ind=pix,
-                        theta=cur_hwp_angle / 2,  # hwp angle returns 2ωt
+                        theta=cur_hwp_angle,
                         psi=psi,
                         xi=xi,
                         maps=self.maps,
@@ -1315,7 +1314,7 @@ def fill_tod(
                         mUQ=self.mUQ,
                         mQU=self.mQU,
                         pixel_ind=pix,
-                        theta=cur_hwp_angle / 2,  # hwp angle returns 2ωt
+                        theta=cur_hwp_angle,
                         psi=psi,
                         xi=xi,
                         maps=self.maps,
@@ -1340,7 +1339,7 @@ def fill_tod(
                                 mUQs=self.mUQs,
                                 mQUs=self.mQUs,
                                 pixel_ind=pix,
-                                theta=cur_hwp_angle / 2,  # hwp angle returns 2ωt
+                                theta=cur_hwp_angle,
                                 psi=psi,
                                 xi=xi,
                             )
@@ -1359,19 +1358,15 @@ def fill_tod(
                                 mUQs=self.mUQs,
                                 mQUs=self.mQUs,
                                 pixel_ind=pix,
-                                theta=cur_hwp_angle / 2,  # hwp angle returns 2ωt
+                                theta=cur_hwp_angle,
                                 psi=psi,
                                 xi=xi,
                             )
 
                     else:
                         # in this case factor 4 included here
-                        ca = np.cos(
-                            2 * cur_point[:, 2] + 4 * cur_hwp_angle / 2
-                        )  # hwp angle returns 2ωt
-                        sa = np.sin(
-                            2 * cur_point[:, 2] + 4 * cur_hwp_angle / 2
-                        )  # hwp angle returns 2ωt
+                        ca = np.cos(2 * cur_point[:, 2] + 4 * cur_hwp_angle)
+                        sa = np.sin(2 * cur_point[:, 2] + 4 * cur_hwp_angle)
 
                         self.atd[pix, 0] += tod * 0.5
                         self.atd[pix, 1] += tod * ca * 0.5

litebird_sim/mapmaking/common.py

@@ -205,15 +205,14 @@ def _compute_pixel_indices(
             curr_pointings_det = pointings[idet, :, :]
         else:
             curr_pointings_det, hwp_angle = pointings(idet)
-            curr_pointings_det = curr_pointings_det.reshape(-1, 3)
 
         if hwp_angle is None:
             hwp_angle = 0
 
         if output_coordinate_system == CoordinateSystem.Galactic:
             curr_pointings_det = rotate_coordinates_e2g(curr_pointings_det)
 
-        polang_all[idet] = curr_pointings_det[:, 2] + hwp_angle
+        polang_all[idet] = curr_pointings_det[:, 2] + 2 * hwp_angle
 
         pixidx_all[idet] = hpx.ang2pix(curr_pointings_det[:, :2])
 

litebird_sim/observations.py

@@ -680,20 +680,30 @@ def get_pointings(
 
             # All the detectors are included
             pointings, hwp_angle = cur_obs.get_pointings("all")
+            # This returns ``(N_det, N_samples, 3)`` array
 
             # Only the first two detectors are included
             pointings, hwp_angle = cur_obs.get_pointings([0, 1])
+            # This returns ``(2, N_samples, 3)`` array
 
             # Only the first detector is used
             pointings, hwp_angle = cur_obs.get_pointings(0)
+            # This returns ``(N_samples, 3)`` array
+
+            # NB if a list of indices is passed an array of dimension ``(N_det, N_samples, 3)``
+            # is always returned
+            # For example this
+            pointings, hwp_angle = cur_obs.get_pointings([0])
+            # returns ``(1, N_samples, 3)`` array
 
         The return value is a pair containing (1) the pointing matrix and (2) the
         HWP angle. The pointing matrix is a NumPy array with shape ``(N_det, N_samples, 3)``,
         where ``N_det` is the number of detectors and ``N_samples`` is the number of
         samples in the TOD (the field ``Observation.n_samples``). The last dimension
         spans the three angles θ (colatitude, in radians), φ (longitude, in radians),
         and ψ (orientation angle, in radians). *Important*: if you ask for just *one*
-        detector, the shape of the pointing matrix will always be ``(N_samples, 3)``.
+        detector passing the index of the detector, the shape of the pointing matrix
+        will always be ``(N_samples, 3)``.
         The HWP angle is always a vector with shape ``(N_samples,)``, as it does
         not depend on the list of detectors.
 
@@ -712,7 +722,7 @@ def get_pointings(
                 (detector_idx >= 0) and (detector_idx < self.n_detectors)
             ), f"Invalid detector index {detector_idx}, it must be a number between 0 and {self.n_detectors - 1}"
 
-            pointing, hwp = self.pointing_provider.get_pointings(
+            return self.pointing_provider.get_pointings(
                 detector_quat=self.quat[detector_idx],
                 start_time=self.start_time,
                 start_time_global=self.start_time_global,
@@ -722,7 +732,6 @@ def get_pointings(
                 hwp_buffer=hwp_buffer,
                 pointings_dtype=pointings_dtype,
             )
-            return pointing.reshape(1, -1, 3), hwp
 
         # More complex case: we need all the detectors
         if isinstance(detector_idx, str):

litebird_sim/scan_map.py

@@ -63,7 +63,6 @@ def scan_map(
             curr_pointings_det = pointings[detector_idx, :, :]
         else:
             curr_pointings_det, hwp_angle = pointings(detector_idx)
-            curr_pointings_det = curr_pointings_det.reshape(-1, 3)
 
         if hwp_angle is None:
             hwp_angle = 0
@@ -87,7 +86,7 @@ def scan_map(
                 input_T=maps_det[0, pixel_ind_det],
                 input_Q=maps_det[1, pixel_ind_det],
                 input_U=maps_det[2, pixel_ind_det],
-                pol_angle_det=curr_pointings_det[:, 2] + hwp_angle,
+                pol_angle_det=curr_pointings_det[:, 2] + 2 * hwp_angle,
             )
 
         elif interpolation == "linear":
@@ -102,7 +101,7 @@ def scan_map(
                 input_U=hp.get_interp_val(
                     maps_det[2, :], curr_pointings_det[:, 0], curr_pointings_det[:, 1]
                 ),
-                pol_angle_det=curr_pointings_det[:, 2] + hwp_angle,
+                pol_angle_det=curr_pointings_det[:, 2] + 2 * hwp_angle,
             )
 
         else: