From 2c8002f65b21ea22b837b863bb89b8c2229f7dda Mon Sep 17 00:00:00 2001 From: keskitalo Date: Tue, 6 Aug 2024 12:22:11 -0700 Subject: [PATCH 1/3] Add option in the detector pointing operator to deflect the pointing based on the HWP angle --- .../pointing_detector/pointing_detector.py | 28 +++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/src/toast/ops/pointing_detector/pointing_detector.py b/src/toast/ops/pointing_detector/pointing_detector.py index 91a1407a2..918132aa6 100644 --- a/src/toast/ops/pointing_detector/pointing_detector.py +++ b/src/toast/ops/pointing_detector/pointing_detector.py @@ -1,7 +1,8 @@ -# 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 @@ -9,7 +10,7 @@ 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,17 @@ 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_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 +238,18 @@ 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_quat = qa.rotation(zaxis, hwp_angle) + det_quats = ob.detdata[self.quats].data + det_quats[:] = qa.mult(det_quats, qa.mult(hwp_quat, deflection)) + return def _finalize(self, data, **kwargs): From f4710fe5e09bc5a12a4c5ec8fc7227433f336f78 Mon Sep 17 00:00:00 2001 From: Reijo Keskitalo Date: Fri, 9 Aug 2024 09:54:53 -0700 Subject: [PATCH 2/3] Make toast the first import --- workflows/toast_sim_ground.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/workflows/toast_sim_ground.py b/workflows/toast_sim_ground.py index 7a29ed0d8..240a36c99 100644 --- a/workflows/toast_sim_ground.py +++ b/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 From 522b284ecea462d4c221238b842e47fbd9b33a9a Mon Sep 17 00:00:00 2001 From: keskitalo Date: Thu, 15 Aug 2024 11:53:06 -0700 Subject: [PATCH 3/3] Add HWP phase and unit test --- .../pointing_detector/pointing_detector.py | 28 +++--- src/toast/tests/CMakeLists.txt | 1 + src/toast/tests/ops_pointing_detector.py | 88 +++++++++++++++++++ src/toast/tests/runner.py | 2 + 4 files changed, 108 insertions(+), 11 deletions(-) create mode 100644 src/toast/tests/ops_pointing_detector.py diff --git a/src/toast/ops/pointing_detector/pointing_detector.py b/src/toast/ops/pointing_detector/pointing_detector.py index 918132aa6..29712a607 100644 --- a/src/toast/ops/pointing_detector/pointing_detector.py +++ b/src/toast/ops/pointing_detector/pointing_detector.py @@ -51,6 +51,10 @@ class PointingDetectorSimple(Operator): 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, @@ -238,17 +242,19 @@ 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_quat = qa.rotation(zaxis, hwp_angle) - det_quats = ob.detdata[self.quats].data - det_quats[:] = qa.mult(det_quats, qa.mult(hwp_quat, deflection)) + # 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 diff --git a/src/toast/tests/CMakeLists.txt b/src/toast/tests/CMakeLists.txt index 9aa3745b2..b61e7f47e 100644 --- a/src/toast/tests/CMakeLists.txt +++ b/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 diff --git a/src/toast/tests/ops_pointing_detector.py b/src/toast/tests/ops_pointing_detector.py new file mode 100644 index 000000000..4a1243165 --- /dev/null +++ b/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) diff --git a/src/toast/tests/runner.py b/src/toast/tests/runner.py index a1139e3ea..8bc091d62 100644 --- a/src/toast/tests/runner.py +++ b/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))