Use quaternions for mapmaking

maps/include/maps/pointing.h

@@ -7,6 +7,10 @@
 
 // XXX: Define what the following functions do!
 
+G3VectorQuat
+get_detector_pointing_quats(double x_offset, double y_offset,
+    const G3VectorQuat & trans_quat, MapCoordReference coord_sys);
+
 void get_detector_pointing( 
 	double x_offset, double y_offset,
 	const G3VectorQuat & trans_quat,

maps/src/HitsBinner.cxx

@@ -256,11 +256,10 @@ HitsBinner::BinHits(const BolometerProperties &bp, const G3VectorQuat &pointing,
 {
 	// Get per-detector pointing timestream
 
-	std::vector<double> alpha, delta;
-	get_detector_pointing(bp.x_offset, bp.y_offset, pointing, H->coord_ref,
-	    alpha, delta);
+	auto detquats = get_detector_pointing_quats(bp.x_offset, bp.y_offset,
+	    pointing, H->coord_ref);
 
-	auto detpointing = H->AnglesToPixels(alpha, delta);
+	auto detpointing = H->QuatsToPixels(detquats);
 
 	for (size_t i = 0; i < detpointing.size(); i++) {
 		(*H)[detpointing[i]] += 1.0;

maps/src/MapBinner.cxx

@@ -372,11 +372,10 @@ MapBinner::BinTimestream(const G3Timestream &det, double weight,
 {
 	// Get per-detector pointing timestream
 
-	std::vector<double> alpha, delta;
-	get_detector_pointing(bp.x_offset, bp.y_offset, pointing, T->coord_ref,
-	    alpha, delta);
+	auto detquats = get_detector_pointing_quats(bp.x_offset, bp.y_offset,
+	    pointing, T->coord_ref);
 
-	auto detpointing = T->AnglesToPixels(alpha, delta);
+	auto detpointing = T->QuatsToPixels(detquats);
 
 	if (Q) {
 		// And polarization coupling

maps/src/MapMockObserver.cxx

@@ -202,12 +202,11 @@ MapMockObserver::Process(G3FramePtr frame, std::deque<G3FramePtr> &out)
 		det.stop = pointing->stop;
 
 		// Get per-detector pointing timestream
-		std::vector<double> alpha, delta;
-		get_detector_pointing(bp.x_offset, bp.y_offset, *pointing,
-		    T_->coord_ref, alpha, delta);
+		auto detquats = get_detector_pointing_quats(bp.x_offset, bp.y_offset,
+		    *pointing, T_->coord_ref);
 		std::vector<size_t> detpointing;
 		if (!interp_)
-			detpointing = T_->AnglesToPixels(alpha, delta);
+			detpointing = T_->QuatsToPixels(detquats);
 
 		if (Q_) {
 			StokesVector pcoupling;
@@ -217,7 +216,7 @@ MapMockObserver::Process(G3FramePtr frame, std::deque<G3FramePtr> &out)
 				if (interp_) {
 					std::vector<size_t> pixels;
 					std::vector<double> weights;
-					T_->GetInterpPixelsWeights(alpha[i], delta[i], pixels, weights);
+					T_->GetInterpPixelsWeights(detquats[i], pixels, weights);
 					det[i] =
 					    T_->GetInterpPrecalc(pixels, weights) * pcoupling.t +
 					    Q_->GetInterpPrecalc(pixels, weights) * pcoupling.q +
@@ -232,7 +231,7 @@ MapMockObserver::Process(G3FramePtr frame, std::deque<G3FramePtr> &out)
 		} else {
 			for (size_t i = 0; i < det.size(); i++) {
 				if (interp_)
-					det[i] = T_->GetInterpValue(alpha[i], delta[i]);
+					det[i] = T_->GetInterpValue(detquats[i]);
 				else
 					det[i] = T_->at(detpointing[i]);
 			}

maps/src/MapTODMasker.cxx

@@ -105,10 +105,9 @@ MapTODMasker::Process(G3FramePtr frame, std::deque<G3FramePtr> &out)
 		const BolometerProperties &bp = boloprops_->at(dets[i]);
 
 		// Get per-detector pointing timestream
-		std::vector<double> alpha, delta;
-		get_detector_pointing(bp.x_offset, bp.y_offset, *pointing,
-		    skymap->coord_ref, alpha, delta);
-		auto detpointing = skymap->AnglesToPixels(alpha, delta);
+		auto detquats = get_detector_pointing_quats(bp.x_offset, bp.y_offset,
+		    *pointing, skymap->coord_ref);
+		auto detpointing = skymap->QuatsToPixels(detquats);
 
 		det.resize(pointing->size());
 		for (size_t j = 0; j < det.size(); j++)

maps/src/MapTODPointing.cxx

@@ -152,15 +152,14 @@ MapTODPointing::Process(G3FramePtr frame, std::deque<G3FramePtr> &out)
 	#pragma omp parallel for
 #endif
 	for (size_t i = 0; i < dets.size(); i++) {
-		std::vector<int> &det = output->at(dets[i]);
+		auto &det = output->at(dets[i]);
 		const BolometerProperties &bp = boloprops_->at(dets[i]);
 
 		// Get per-detector pointing timestream
-		std::vector<double> alpha, delta;
-		get_detector_pointing(bp.x_offset, bp.y_offset, *pointing,
-		    stub_->coord_ref, alpha, delta);
+		auto detquats = get_detector_pointing_quats(bp.x_offset, bp.y_offset,
+		    *pointing, stub_->coord_ref);
 
-		auto detpointing = stub_->AnglesToPixels(alpha, delta);
+		auto detpointing = stub_->QuatsToPixels(detquats);
 
 		for (size_t j = 0; j < sz; j++)
 			det[j] = detpointing[j];

maps/src/SingleDetectorBoresightBinner.cxx

@@ -9,6 +9,7 @@
 #include <G3Data.h>
 #include <G3Map.h>
 #include <maps/G3SkyMap.h>
+#include <maps/pointing.h>
 
 class SingleDetectorBoresightBinner : public G3Module {
 public:
@@ -183,7 +184,8 @@ SingleDetectorBoresightBinner::Process(G3FramePtr frame,
 	g3_assert(timestreams->NSamples() == pointing->size());
 
 	// Conjugate pointing rotation with boresight vector
-	auto pointing_vec = (*pointing)*quat(0,1,0,0)/(*pointing);
+	auto pointing_vec = get_detector_pointing_quats(0, 0,
+	    *pointing, template_->coord_ref);
 	auto pixels = template_->QuatsToPixels(pointing_vec);
 
 	for (size_t i = 0; i < pixels.size(); i++)

maps/src/SingleDetectorMapBinner.cxx

@@ -156,14 +156,14 @@ SingleDetectorMapBinner::Process(G3FramePtr frame,
 		}
 
 		// Get per-detector pointing timestream
-        	std::vector<double> alpha, delta;
 		auto bp = boloprops_->find(det);
 		if (bp == boloprops_->end())
 			log_fatal("Missing bolometer properties for %s",
 			    det.c_str());
-		get_detector_pointing(bp->second.x_offset, bp->second.y_offset,
-		    *pointing, m->coord_ref, alpha, delta);
-		auto pixels = m->AnglesToPixels(alpha, delta);
+		auto detquats = get_detector_pointing_quats(
+		    bp->second.x_offset, bp->second.y_offset,
+		    *pointing, m->coord_ref);
+		auto pixels = m->QuatsToPixels(detquats);
 
 		g3_assert(ts->size() == pixels.size());
 		for (size_t j = 0; j < ts->size(); j++) {

maps/src/pointing.cxx

@@ -361,6 +361,28 @@ convert_ra_dec_trans_to_gal(const G3VectorQuat &radec_trans,
 		gal_trans[i] = gt*radec_trans[i];
 }
 
+G3VectorQuat
+get_detector_pointing_quats(double x_offset, double y_offset,
+    const G3VectorQuat &trans_quat, MapCoordReference coord_sys)
+{
+	quat q_off = offsets_to_quat(x_offset, y_offset);
+	size_t nsamp = trans_quat.size();
+	G3VectorQuat det_quats(nsamp, quat(0, 1, 0, 0));
+
+	for (size_t i = 0; i < nsamp; i++)
+		det_quats[i] = trans_quat[i] * q_off * ~trans_quat[i];
+
+	if (coord_sys == Local) {
+		for (size_t i = 0; i < nsamp; i++) {
+			const quat &q = det_quats[i];
+			det_quats[i] = quat(q.R_component_1(), q.R_component_2(),
+			    q.R_component_3(), -q.R_component_4());
+		}
+	}
+
+	return det_quats;
+}
+
 void
 get_detector_pointing(double x_offset, double y_offset,
     const G3VectorQuat &trans_quat, MapCoordReference coord_sys,