Add option in the detector pointing operator to deflect the pointing … (

#776)

* Add option in the detector pointing operator to deflect the pointing based on the HWP angle

* Make toast the first import

* Add HWP phase and unit test

src/toast/ops/pointing_detector/pointing_detector.py

@@ -1,15 +1,16 @@
-# Copyright (c) 2015-2020 by the parties listed in the AUTHORS file.
+# Copyright (c) 2015-2024 by the parties listed in the AUTHORS file.
 # All rights reserved.  Use of this source code is governed by
 # a BSD-style license that can be found in the LICENSE file.
 
+import astropy.units as u
 import numpy as np
 import traitlets
 
 from ... import qarray as qa
 from ...accelerator import ImplementationType
 from ...observation import default_values as defaults
 from ...timing import function_timer
-from ...traits import Bool, Int, Unicode, UseEnum, trait_docs
+from ...traits import Bool, Int, Quantity, Unicode, UseEnum, trait_docs
 from ...utils import Logger
 from ..operator import Operator
 from .kernels import pointing_detector
@@ -46,6 +47,21 @@ class PointingDetectorSimple(Operator):
         defaults.boresight_radec, help="Observation shared key for boresight"
     )
 
+    hwp_angle = Unicode(
+        defaults.hwp_angle, allow_none=True, help="Observation shared key for HWP angle"
+    )
+
+    hwp_angle_offset = Quantity(
+        0 * u.deg, help="HWP angle offset to apply when constructing deflection"
+    )
+
+    hwp_deflection_radius = Quantity(
+        None,
+        allow_none=True,
+        help="If non-zero, nominal detector pointing will be deflected in a circular "
+        "pattern according to HWP phase.",
+    )
+
     quats = Unicode(
         defaults.quats,
         allow_none=True,
@@ -226,6 +242,20 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                 use_accel=use_accel,
             )
 
+            # Optionally apply HWP deflection
+            if self.hwp_deflection_radius is not None and \
+               self.hwp_deflection_radius.value != 0:
+                xaxis, yaxis, zaxis = np.eye(3)
+                deflection = qa.rotation(
+                    xaxis, self.hwp_deflection_radius.to_value(u.radian)
+                )
+                hwp_angle = ob.shared[self.hwp_angle].data
+                hwp_angle += self.hwp_angle_offset.to_value(u.rad)
+                hwp_quat = qa.rotation(zaxis, hwp_angle)
+                det_quats = ob.detdata[self.quats].data
+                for idet in range(len(dets)):
+                    det_quats[idet] = qa.mult(det_quats[idet], qa.mult(hwp_quat, deflection))
+
         return
 
     def _finalize(self, data, **kwargs):

src/toast/tests/CMakeLists.txt

@@ -37,6 +37,7 @@ install(FILES
     ops_mapmaker_solve.py
     ops_mapmaker.py
     covariance.py
+    ops_pointing_detector.py
     ops_pointing_healpix.py
     ops_pointing_wcs.py
     ops_memory_counter.py

src/toast/tests/ops_pointing_detector.py

@@ -0,0 +1,88 @@
+# Copyright (c) 2024 by the parties listed in the AUTHORS file.
+# All rights reserved.  Use of this source code is governed by
+# a BSD-style license that can be found in the LICENSE file.
+
+import os
+
+import astropy.units as u
+import healpy as hp
+import numpy as np
+
+from .. import ops as ops
+from .. import qarray as qa
+from .._libtoast import pixels_healpix, stokes_weights_IQU
+from ..accelerator import ImplementationType, accel_enabled
+from ..intervals import IntervalList, interval_dtype
+from ..observation import default_values as defaults
+from ._helpers import (
+    close_data, create_outdir, create_satellite_data, create_ground_data
+)
+from .mpi import MPITestCase
+
+
+class PointingDetectorTest(MPITestCase):
+    def setUp(self):
+        fixture_name = os.path.splitext(os.path.basename(__file__))[0]
+        self.outdir = create_outdir(self.comm, fixture_name)
+
+    def test_detector_pointing_simple(self):
+        # Create a fake satellite data set for testing
+        data = create_satellite_data(self.comm)
+
+        detpointing = ops.PointingDetectorSimple()
+        detpointing.apply(data)
+
+        # Also make a copy using a python codepath
+        detpointing.kernel_implementation = ImplementationType.NUMPY
+        detpointing.quats = "pyquat"
+        detpointing.apply(data)
+
+        for ob in data.obs:
+            np.testing.assert_array_equal(
+                ob.detdata[defaults.quats], ob.detdata["pyquat"]
+            )
+
+        rank = 0
+        if self.comm is not None:
+            rank = self.comm.rank
+
+        handle = None
+        if rank == 0:
+            handle = open(os.path.join(self.outdir, "out_test_detpointing_simple_info"), "w")
+        data.info(handle=handle)
+        if rank == 0:
+            handle.close()
+
+        close_data(data)
+
+    def test_detector_pointing_hwp_deflect(self):
+        # Create fake observing of a small patch
+        data = create_ground_data(self.comm)
+
+        # Regular pointing
+        detpointing1 = ops.PointingDetectorSimple()
+        detpointing1.apply(data)
+
+        # Pointing with deflection
+        detpointing2 = ops.PointingDetectorSimple(
+            hwp_angle=defaults.hwp_angle,
+            hwp_angle_offset=15 * u.deg,
+            hwp_deflection_radius=1.0 * u.deg,
+            quats="deflected",
+        )
+        detpointing2.apply(data)
+
+        # Compare
+        for ob in data.obs:
+            quats1 = ob.detdata[defaults.quats]
+            quats2 = ob.detdata["deflected"]
+            ndet = quats1.shape[0]
+            for idet in range(ndet):
+                theta1, phi1, psi1 = qa.to_iso_angles(quats1[idet])
+                theta2, phi2, psi2 = qa.to_iso_angles(quats2[idet])
+                rms_theta = np.degrees(np.std(theta1 - theta2))
+                rms_phi = np.degrees(np.std(phi1 - phi2))
+                self.assertTrue(rms_theta > .5 and rms_theta < 1.5)
+                self.assertTrue(rms_phi > .5 and rms_phi < 1.5)
+
+        close_data(data)

src/toast/tests/runner.py

@@ -45,6 +45,7 @@
 from . import ops_noise_estim as test_ops_noise_estim
 from . import ops_perturbhwp as test_ops_perturbhwp
 from . import ops_pixels_healpix as test_ops_pixels_healpix
+from . import ops_pointing_detector as test_ops_pointing_detector
 from . import ops_pointing_healpix as test_ops_pointing_healpix
 from . import ops_pointing_wcs as test_ops_pointing_wcs
 from . import ops_polyfilter as test_ops_polyfilter
@@ -181,6 +182,7 @@ def test(name=None, verbosity=2):
         suite.addTest(loader.loadTestsFromModule(test_ops_sim_satellite))
         suite.addTest(loader.loadTestsFromModule(test_ops_sim_ground))
         suite.addTest(loader.loadTestsFromModule(test_ops_memory_counter))
+        suite.addTest(loader.loadTestsFromModule(test_ops_pointing_detector))
         suite.addTest(loader.loadTestsFromModule(test_ops_pointing_healpix))
         suite.addTest(loader.loadTestsFromModule(test_ops_pointing_wcs))
         suite.addTest(loader.loadTestsFromModule(test_ops_sim_tod_noise))

workflows/toast_sim_ground.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python3
 
-# Copyright (c) 2015-2021 by the parties listed in the AUTHORS file.
+# Copyright (c) 2015-2024 by the parties listed in the AUTHORS file.
 # All rights reserved.  Use of this source code is governed by
 # a BSD-style license that can be found in the LICENSE file.
 
@@ -25,6 +25,8 @@
 
 """
 
+import toast
+
 import argparse
 import datetime
 import os