From 9ff6896df9a622befc6f56793b12fa54de594e0f Mon Sep 17 00:00:00 2001
From: aashrita <aashrita4@gmail.com>
Date: Mon, 13 Jul 2026 20:41:14 -0700
Subject: [PATCH] Add rotations for focal plane and boresight in planet
 mapmaker. Allows for HWP deflection to be applied cleanly.

---
 sotodlib/coords/helpers.py | 35 ++++++++++++++++++++++++++---------
 sotodlib/coords/planets.py | 11 ++++++++++-
 2 files changed, 36 insertions(+), 10 deletions(-)

diff --git a/sotodlib/coords/helpers.py b/sotodlib/coords/helpers.py
index 7dd069616..4f89fe38c 100644
--- a/sotodlib/coords/helpers.py
+++ b/sotodlib/coords/helpers.py
@@ -739,6 +739,30 @@ def to_g3(self):
             return so3g.proj.quat.G3VectorQuat(temp)
         raise ValueError("Can only convert 1- or 2-d arrays to G3.")
 
+def get_deflection_quat(aman, wobble_meta):
+    '''
+    Derives the deflection quaternion.
+
+    Parameters
+    ----------
+    aman : AxisManager
+
+    wobble_meta : AxisManager
+
+    '''
+    # the amp and phase are the same for a given wafer, so we can take any of them, in this case for detector index 0
+    # !!!!! this won't work for mixing more than one wafer.
+    # the metadata has amplitudes in arcmin, and phases in radians
+    amp = wobble_meta.amp[0]/60.*np.pi/180.0
+    phase = wobble_meta.phase[0]
+
+    dxi = amp * np.cos(aman.hwp_angle - phase)
+    deta = -amp * np.sin(aman.hwp_angle - phase)
+    deflq = so3g.proj.quat.rotation_xieta(xi=dxi, eta=deta)
+
+    return deflq
+
+
 def get_deflected_sightline(aman, wobble_meta, site='so', weather='typical'):
     """
     Constructs a deflected CelestialSightLine using HWP-synchronous
@@ -776,15 +800,8 @@ def get_deflected_sightline(aman, wobble_meta, site='so', weather='typical'):
 
     if len(wafer_slots) != 1 or len(bands) != 1:
         raise ValueError("Detectors span multiple wafer_slots or bands.")
-    # the amp and phase are the same for a given wafer, so we can take any of them, in this case for detector index 0
-    # !!!!! this won't work for mixing more than one wafer.
-    # the metadata has amplitudes in arcmin, and phases in radians
-    amp = wobble_meta.amp[0]/60.*np.pi/180.0
-    phase = wobble_meta.phase[0]
-
-    dxi = amp * np.cos(aman.hwp_angle - phase)
-    deta = -amp * np.sin(aman.hwp_angle - phase)
-    deflq = so3g.proj.quat.rotation_xieta(xi=dxi, eta=deta)
+    
+    deflq = get_deflection_quat(aman, wobble_meta)
 
     sight = so3g.proj.CelestialSightLine.az_el(
         aman.timestamps,
diff --git a/sotodlib/coords/planets.py b/sotodlib/coords/planets.py
index 5b2e4ad83..33802bd30 100644
--- a/sotodlib/coords/planets.py
+++ b/sotodlib/coords/planets.py
@@ -762,7 +762,8 @@ def make_map(tod, center_on=None, scan_coords=True, thread_algo=False,
              filename=None, source_flags=None, cuts=None,
              data_splits=None,
              low_pass=None, n_modes=10,
-             eigentol=1e-3, info={}):
+             eigentol=1e-3, info={},
+             rot_bs=None, rot_fp=None):
     """Make a compact source map from the TOD.  Specify filename to write
     things to disk; this should be a format string, for example
     '{obs_id}_{map}.fits', where 'map' will be given values of
@@ -804,6 +805,14 @@ class MmTimer(coords.Timer):
                                  cuts=cuts,
                                  threads=thread_algo,
                                  wcs_kernel=wcsk)
+        
+        # apply an optional boresight rotation or focal plane rotation
+        # e.g. HWP deflection is rot_fp
+        if rot_bs is not None:
+            P.sight.Q = rot_bs * P.sight.Q
+        if rot_fp is not None:
+            P.sight.Q = P.sight.Q * rot_fp
+            
     with MmTimer('get_proj_threads'):
         P._get_proj_threads()