From 5581882da200e2414fa3336719de78d1547f76a9 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 12 Oct 2023 11:07:21 -0500 Subject: [PATCH] Use quaternions for mapmaking --- maps/include/maps/pointing.h | 4 ++++ maps/src/HitsBinner.cxx | 7 +++---- maps/src/MapBinner.cxx | 7 +++---- maps/src/MapMockObserver.cxx | 11 +++++------ maps/src/MapTODMasker.cxx | 7 +++---- maps/src/MapTODPointing.cxx | 9 ++++----- maps/src/SingleDetectorBoresightBinner.cxx | 4 +++- maps/src/SingleDetectorMapBinner.cxx | 8 ++++---- maps/src/pointing.cxx | 22 ++++++++++++++++++++++ 9 files changed, 51 insertions(+), 28 deletions(-) diff --git a/maps/include/maps/pointing.h b/maps/include/maps/pointing.h index 647ad4c5..ecb47d48 100644 --- a/maps/include/maps/pointing.h +++ b/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, diff --git a/maps/src/HitsBinner.cxx b/maps/src/HitsBinner.cxx index ae6533dd..3f3a0bef 100644 --- a/maps/src/HitsBinner.cxx +++ b/maps/src/HitsBinner.cxx @@ -256,11 +256,10 @@ HitsBinner::BinHits(const BolometerProperties &bp, const G3VectorQuat &pointing, { // Get per-detector pointing timestream - std::vector 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; diff --git a/maps/src/MapBinner.cxx b/maps/src/MapBinner.cxx index d47665d7..5623730c 100644 --- a/maps/src/MapBinner.cxx +++ b/maps/src/MapBinner.cxx @@ -372,11 +372,10 @@ MapBinner::BinTimestream(const G3Timestream &det, double weight, { // Get per-detector pointing timestream - std::vector 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 diff --git a/maps/src/MapMockObserver.cxx b/maps/src/MapMockObserver.cxx index ba81012d..befc928f 100644 --- a/maps/src/MapMockObserver.cxx +++ b/maps/src/MapMockObserver.cxx @@ -202,12 +202,11 @@ MapMockObserver::Process(G3FramePtr frame, std::deque &out) det.stop = pointing->stop; // Get per-detector pointing timestream - std::vector 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 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 &out) if (interp_) { std::vector pixels; std::vector 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 &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]); } diff --git a/maps/src/MapTODMasker.cxx b/maps/src/MapTODMasker.cxx index 0c956ca6..0212fbc1 100644 --- a/maps/src/MapTODMasker.cxx +++ b/maps/src/MapTODMasker.cxx @@ -105,10 +105,9 @@ MapTODMasker::Process(G3FramePtr frame, std::deque &out) const BolometerProperties &bp = boloprops_->at(dets[i]); // Get per-detector pointing timestream - std::vector 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++) diff --git a/maps/src/MapTODPointing.cxx b/maps/src/MapTODPointing.cxx index 5dd11208..dd243174 100644 --- a/maps/src/MapTODPointing.cxx +++ b/maps/src/MapTODPointing.cxx @@ -152,15 +152,14 @@ MapTODPointing::Process(G3FramePtr frame, std::deque &out) #pragma omp parallel for #endif for (size_t i = 0; i < dets.size(); i++) { - std::vector &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 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]; diff --git a/maps/src/SingleDetectorBoresightBinner.cxx b/maps/src/SingleDetectorBoresightBinner.cxx index 6c9b5d4f..2170eddc 100644 --- a/maps/src/SingleDetectorBoresightBinner.cxx +++ b/maps/src/SingleDetectorBoresightBinner.cxx @@ -9,6 +9,7 @@ #include #include #include +#include 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++) diff --git a/maps/src/SingleDetectorMapBinner.cxx b/maps/src/SingleDetectorMapBinner.cxx index 93f85933..c485b0ff 100644 --- a/maps/src/SingleDetectorMapBinner.cxx +++ b/maps/src/SingleDetectorMapBinner.cxx @@ -156,14 +156,14 @@ SingleDetectorMapBinner::Process(G3FramePtr frame, } // Get per-detector pointing timestream - std::vector 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++) { diff --git a/maps/src/pointing.cxx b/maps/src/pointing.cxx index 24b9b832..15b7ab5d 100644 --- a/maps/src/pointing.cxx +++ b/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,