More cleanup of pointing functions

CMB-S4 · Oct 17, 2023 · 2dca0a0 · 2dca0a0
1 parent d77404a
commit 2dca0a0
Show file tree

Hide file tree

Showing 6 changed files with 181 additions and 415 deletions.
diff --git a/maps/include/maps/pointing.h b/maps/include/maps/pointing.h
@@ -5,40 +5,49 @@
 
 #include <maps/G3SkyMap.h>
 
-// XXX: Define what the following functions do!
-
+// Compute a vector of detector pointing quaternions from a vector of
+// boresight transform quaternions
 G3VectorQuat
 get_detector_pointing_quats(double x_offset, double y_offset,
-    const G3VectorQuat & trans_quat, MapCoordReference coord_sys);
+    const G3VectorQuat &trans_quat, MapCoordReference coord_sys);
 
+// Compute a vector of detector pointing pixels from a vector of
+// boresight transform quaternions
 std::vector<size_t>
 get_detector_pointing_pixels(double x_offset, double y_offset,
-    const G3VectorQuat & trans_quat, G3SkyMapConstPtr skymap);
-
-void get_detector_pointing( 
-	double x_offset, double y_offset,
-	const G3VectorQuat & trans_quat,
-	MapCoordReference coord_sys,
-	std::vector<double> & alpha, std::vector<double> & delta);
-
-void get_detector_rotation( 
-	double x_offset, double y_offset,
-	const G3VectorQuat & trans_quat,
-	std::vector<double> & rot);
-
-quat get_transform_quat(double as_0, double ds_0,
-			double ae_0, double de_0,
-			double as_1, double ds_1,
-			double ae_1, double de_1);
+    const G3VectorQuat &trans_quat, G3SkyMapConstPtr skymap);
+
+// Compute vectors of detector pointing angles from a vector of
+// boresight transform quaternions
+void
+get_detector_pointing(double x_offset, double y_offset,
+    const G3VectorQuat &trans_quat, MapCoordReference coord_sys,
+    std::vector<double> &alpha, std::vector<double> &delta);
+
+// Compute a vector of detector pointing rotation angles from a vector of
+// boresight transform quaternions
+std::vector<double>
+get_detector_rotation(double x_offset, double y_offset,
+    const G3VectorQuat &trans_quat);
+
+// Compute a vector quaternion that is the boresight rotated
+// by the given x and y offsets.
 quat offsets_to_quat(double x_offset, double y_offset);
-quat get_fk5_j2000_to_gal_quat();
-void quat_to_ang(quat q, double & alpha, double & delta);
+
+// Conversion functions between sky coordinates and vector quaternions
+void quat_to_ang(quat q, double &alpha, double &delta);
 quat ang_to_quat(double alpha, double delta);
-quat coord_quat_to_delta_hat(quat q);
+
+// Compute the angular separation between two vector quaternions
+double quat_ang_sep(quat q0, quat q1);
+
+// Compute a rotation quaternion that would rotate the boresight vector
+// to point in the given sky direction.
 quat get_origin_rotator(double alpha, double delta);
-double get_rot_ang(quat start_q, quat end_q, quat trans);
 
-double angular_distance(double alpha0, double delta0, double alpha1, double delta1);
+// Compute the quaternion for rotating FK5 J2000 coordinates to
+// Galactic J2000 coordinates
+quat get_fk5_j2000_to_gal_quat();
 
 #endif
 
diff --git a/maps/python/coordsysmodules.py b/maps/python/coordsysmodules.py
@@ -1,9 +1,8 @@
 from spt3g import core
 from spt3g.core import G3TimestreamQuat
-from spt3g.maps.azel import convert_azel_to_radec
 from spt3g.maps import MapCoordReference
-from spt3g.maps import create_det_az_el_trans, create_lazy_det_ra_dec_trans
-from spt3g.maps import create_det_ra_dec_trans, convert_ra_dec_trans_to_gal
+from spt3g.maps import get_origin_rotator_timestream, get_boresight_rotator_timestream
+from spt3g.maps import get_fk5_j2000_to_gal_quat
 
 
 __all__ = [
@@ -68,16 +67,18 @@ def FillCoordTransRotations(
         return
 
     trans = G3TimestreamQuat()
-    trans.start = frame[bs_az_key].start
-    trans.stop = frame[bs_az_key].stop
     if end_coord_sys == MapCoordReference.Local:
-        create_det_az_el_trans(frame[bs_az_key], frame[bs_el_key], trans)
+        trans = get_origin_rotator_timestream(
+            frame[bs_az_key], frame[bs_el_key], end_coord_sys
+        )
     else:
         if do_bad_transform:
             core.log_debug("You are doing the old calculation for pointing")
-            create_lazy_det_ra_dec_trans(frame[bs_ra_key], frame[bs_dec_key], trans)
+            trans = get_origin_rotator_timestream(
+                frame[bs_ra_key], frame[bs_dec_key], end_coord_sys
+            )
         else:
-            create_det_ra_dec_trans(
+            trans = get_boresight_rotator_timestream(
                 frame[bs_az_key],
                 frame[bs_el_key],
                 frame[bs_ra_key],
@@ -86,7 +87,6 @@ def FillCoordTransRotations(
                 frame[offset_el_key],
                 frame[offset_ra_key],
                 frame[offset_dec_key],
-                trans,
             )
     frame[transform_store_key] = trans
 
@@ -106,10 +106,7 @@ def EquatorialToGalacticTransRotations(
 
     if frame.type != core.G3FrameType.Scan:
         return
-    gal_trans = G3TimestreamQuat()
-    gal_trans.start = frame[eq_trans_key].start
-    gal_trans.stop = frame[eq_trans_key].stop
-    convert_ra_dec_trans_to_gal(frame[eq_trans_key], gal_trans)
+    gal_trans = get_fk5_j2000_to_gal_quat() * frame[eq_trans_key]
     frame[out_key] = gal_trans
 
 
@@ -129,9 +126,8 @@ def AddLocalTransRotations(
 
     if frame.type != core.G3FrameType.Scan:
         return
-    local_trans = G3TimestreamQuat()
-    local_trans.start = frame[az_key].start
-    local_trans.stop = frame[az_key].stop
-    create_det_az_el_trans(frame[az_key], frame[el_key], local_trans)
+    local_trans = get_origin_rotator_timestream(
+        frame[az_key], frame[el_key], MapCoordReference.Local
+    )
     frame[out_key] = local_trans
 
diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx
@@ -426,18 +426,23 @@ G3SkyMap::QueryAlphaEllipse(quat q, double a, double b) const
 	double rmin = a > b ? b : a;
 	double sd = q.R_component_4();
 	double cd = sqrt((1 - sd) * (1 + sd));
+
+	// focus distance from center
 	double da = ACOS(COS(rmaj) / COS(rmin)) / cd;
 
+	// focus locations
 	quat qda = get_origin_rotator(da, 0);
 	quat ql = qda * q * ~qda;
 	quat qr = ~qda * q * qda;
 
+	// narrow search to pixels within the major disc
 	auto disc = QueryDisc(q, rmaj);
 
+	// narrow further to locus of points within ellipse
 	std::vector<size_t> pixels;
 	for (auto i: disc) {
 		quat qp = PixelToQuat(i);
-		double d = ACOS(dot3(qp, ql)) + ACOS(dot3(qp, qr));
+		double d = quat_ang_sep(qp, ql) + quat_ang_sep(qp, qr);
 		if (d < 2 * rmaj)
 			pixels.push_back(i);
 	}

diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx
@@ -286,7 +286,7 @@ void ReprojMap(G3SkyMapConstPtr in_map, G3SkyMapPtr out_map, int rebin, bool int
 		rotate = true;
 		q_rot = get_fk5_j2000_to_gal_quat();
 		if (in_map->coord_ref == MapCoordReference::Equatorial)
-			q_rot = 1. / q_rot;
+			q_rot = ~q_rot;
 	} else if (in_map->coord_ref != out_map->coord_ref) {
 		log_fatal("Cannot convert input coord_ref %d to output coord_ref %d",
 		    in_map->coord_ref, out_map->coord_ref);
@@ -312,7 +312,7 @@ void ReprojMap(G3SkyMapConstPtr in_map, G3SkyMapPtr out_map, int rebin, bool int
 			double val = 0;
 			auto quats = out_map->GetRebinQuats(i, rebin);
 			if (rotate)
-				quats = q_rot * quats / q_rot;
+				quats = q_rot * quats * ~q_rot;
 			if (interp)
 				for (size_t j = 0; j < quats.size(); j++)
 					val += in_map->GetInterpValue(quats[j]);
@@ -329,7 +329,7 @@ void ReprojMap(G3SkyMapConstPtr in_map, G3SkyMapPtr out_map, int rebin, bool int
 			double val = 0;
 			quat q = out_map->PixelToQuat(i);
 			if (rotate)
-				q = q_rot * q / q_rot;
+				q = q_rot * q * ~q_rot;
 			if (interp)
 				val = in_map->GetInterpValue(q);
 			else