diff --git a/maps/include/maps/pointing.h b/maps/include/maps/pointing.h index 412bd74b..87b27eca 100644 --- a/maps/include/maps/pointing.h +++ b/maps/include/maps/pointing.h @@ -5,40 +5,49 @@ #include -// 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 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 & alpha, std::vector & delta); - -void get_detector_rotation( - double x_offset, double y_offset, - const G3VectorQuat & trans_quat, - std::vector & 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 &alpha, std::vector &delta); + +// Compute a vector of detector pointing rotation angles from a vector of +// boresight transform quaternions +std::vector +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 index 415af5b2..36291e2f 100644 --- 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 index 4598721e..1a4c9191 100644 --- 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 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 index a04d013b..cd397788 100644 --- 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]); diff --git a/maps/src/pointing.cxx b/maps/src/pointing.cxx index 2031f157..888eda5f 100644 --- a/maps/src/pointing.cxx +++ b/maps/src/pointing.cxx @@ -21,8 +21,12 @@ const double PI = constants::pi(); #define ASIN asin #define ATAN2 atan2 -//#define CHECK_QUAT_INVERSE - +#define QNORM(q) \ + { \ + double n = dot3(q, q); \ + if (fabs(n - 1.0) > 1e-6) \ + q /= sqrt(n); \ + } /* * Quaternions cannot represent parity flips. Since celestial coordinates @@ -43,21 +47,12 @@ project_on_plane(quat plane_normal, quat point) quat out_q(point); //ensure unit vec - plane_normal /= sqrt(dot3(plane_normal,plane_normal)); + QNORM(plane_normal); out_q -= plane_normal * dot3(plane_normal, point); + QNORM(out_q); return out_q; } -static bool -sloppy_eq(quat a, quat b, double slop = 1e-6) -{ - // Fuzzy quaternion equality comparison - return ((fabs(a.R_component_1() - (b.R_component_1())) < slop ) && - (fabs(a.R_component_2() - (b.R_component_2())) < slop ) && - (fabs(a.R_component_3() - (b.R_component_3())) < slop ) && - (fabs(a.R_component_4() - (b.R_component_4())) < slop )); -} - quat ang_to_quat(double alpha, double delta) { @@ -71,10 +66,7 @@ ang_to_quat(double alpha, double delta) void quat_to_ang(quat q, double &alpha, double &delta) { - double d = dot3(q,q); - if (fabs(d - 1.0) > 1e-6){ - q /= sqrt(d); - } + QNORM(q); delta = ASIN(q.R_component_4()) * G3Units::rad; alpha = ATAN2(q.R_component_3(), q.R_component_2())*G3Units::rad; } @@ -88,51 +80,53 @@ py_quat_to_ang(quat q) return boost::python::make_tuple(a, d); } -quat +double +quat_ang_sep(quat a, quat b) +{ + QNORM(a); + QNORM(b); + + double d = dot3(a, b); + if (d > 1) + return 0; + if (d < -1) + return PI * G3Units::rad; + return acos(d) * G3Units::rad; +} + +static quat coord_quat_to_delta_hat(quat q) { // computes the delta hat vector for a given point on the unit sphere // specified by q // // (The delta hat is equal to -alpha hat) - - q /= sqrt(dot3(q,q)); + QNORM(q); double st = sqrt(1 - (q.R_component_4()*q.R_component_4())); quat u= quat(0, -1 * (q.R_component_2() * q.R_component_4())/st, -1 * (q.R_component_3() * q.R_component_4())/st, st); - u /= sqrt(dot3(u,u)); + QNORM(u); return u; } -double -get_rot_ang(quat start_q, quat end_q, quat trans) +static double +get_rot_ang(quat start_q, quat trans) { // delta is the physicist spherical coordinates delta // Computes delta hat for the start q applies trans to it // and then computes the angle between that and end_q's delta hat. - quat t = trans * coord_quat_to_delta_hat(start_q) / trans; + quat t = trans * coord_quat_to_delta_hat(start_q) * ~trans; + quat end_q = trans * start_q * ~trans; quat t_p = coord_quat_to_delta_hat(end_q); - - t /= sqrt(dot3(t,t)); - t_p /= sqrt(dot3(t_p,t_p)); - double d = dot3(t,t_p); double sf = (dot3(end_q, cross3(t, t_p)) < 0) ? -1 : 1; - if (d > 1) { - g3_assert(d < 1.01); - return 0; - } else if (d < -1) { - g3_assert(d > -1.01); - return PI * G3Units::rad; - } else { - return sf * acos(d) * G3Units::rad; - } + return sf * quat_ang_sep(t, t_p); } -quat +static 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) { @@ -156,40 +150,24 @@ get_transform_quat(double as_0, double ds_0, double ae_0, double de_0, quat tquat = cross3(asds_0, aede_0); double mag = sqrt(dot3(tquat, tquat)); - double cos_ang = dot3(asds_0, aede_0); - double ang; - if (cos_ang < -1) - ang = PI * G3Units::rad; - else if (cos_ang > 1) - ang = 0; - else - ang = acos(cos_ang); + double ang = quat_ang_sep(asds_0, aede_0); tquat *= sin(ang/2.0) / mag; tquat += quat(cos(ang/2.0),0,0,0); // trans_asds_1 and aede_1 should now be the same up to a rotation // around aede_0 - quat trans_asds_1 = tquat * asds_1 / tquat; + quat trans_asds_1 = tquat * asds_1 * ~tquat; // Project them on to a plane and find the angle between the two vectors // using (ae_0, de_0) as the normal since we are rotating around that // vector. quat p_asds1 = project_on_plane(aede_0, trans_asds_1); quat p_aede1 = project_on_plane(aede_0, aede_1); - p_asds1 /= sqrt(dot3(p_asds1,p_asds1)); - p_aede1 /= sqrt(dot3(p_aede1,p_aede1)); - - double cos_rot_ang = dot3(p_asds1, p_aede1); - double rot_ang; - if (cos_rot_ang < -1) - rot_ang = PI * G3Units::rad; - else if (cos_rot_ang > 1) - rot_ang = 0; - else - rot_ang = acos(cos_rot_ang); + + double rot_ang = quat_ang_sep(p_asds1, p_aede1); double sf = (dot3(aede_0, cross3(p_asds1, p_aede1)) < 0) ? -1 : 1; rot_ang *= sf; - + double sin_rot_ang_ov_2 = sin(rot_ang/2.0); quat rot_quat = quat(cos(rot_ang/2.0), sin_rot_ang_ov_2 * aede_0.R_component_2(), @@ -200,57 +178,15 @@ get_transform_quat(double as_0, double ds_0, double ae_0, double de_0, return final_trans; } -static std::vector -test_trans(double az_0, double el_0, double ra_0, double dec_0, - double az_1, double el_1, double ra_1, double dec_1, - double az_t, double el_t) -{ - // computes the transform from the first 4 variables - // returns that transform to the last 2 variables - - double ra_t, dec_t; - quat q = get_transform_quat(az_0, -el_0, - ra_0, dec_0, - az_1, -el_1, - ra_1, dec_1); - quat azel = ang_to_quat(az_t, -el_t); - quat rad = q * azel / q; - quat_to_ang(rad, ra_t, dec_t); - std::vector r(2,0); - r[0] = ra_t; - r[1] = dec_t; - return r; -} - -static std::vector -test_gal_trans(double az_0, double el_0, double ra_0, double dec_0, - double az_1, double el_1, double ra_1, double dec_1, - double az_t, double el_t) -{ - double l_t, b_t; - quat q = get_fk5_j2000_to_gal_quat() *get_transform_quat(az_0, -el_0, - ra_0, dec_0, - az_1, -el_1, - ra_1, dec_1); - quat azel = ang_to_quat(az_t, -el_t); - quat rad = q * azel / q; - quat_to_ang(rad, l_t, b_t); - std::vector r(2,0); - r[0] = l_t; - r[1] = b_t; - return r; -} - -static double -test_gal_trans_rot(double ra, double dec) +quat +get_fk5_j2000_to_gal_quat() { - quat start = ang_to_quat(ra,dec); - quat trans = get_fk5_j2000_to_gal_quat(); - - return get_rot_ang(start, trans*start/trans, trans); + // returns the quaternion that rotates FK5 J2000 to galactic J2000 coordinates + // return get_transform_quat(0,0, 1.6814025470759737, -1.050488399695429, + // 0,-0.7853981633974483, 5.750520098164818, -1.2109809382060603); + return quat(0.4889475076,-0.483210684,0.1962537583,0.699229742); } - quat offsets_to_quat(double x_offset, double y_offset) { @@ -272,56 +208,32 @@ get_origin_rotator(double alpha, double delta) quat(cos(delta/2.0), 0, -sin(delta/2.0), 0)); } -static void -print_fk5_j2000_to_gal_quat() -{ - // uhh, so, this code was a super lazy way to get the quaternion - // that takes fk5 j2000 to galactic j2000 - std::cout << std::setprecision(10) << std::endl; - std::cout << get_transform_quat( - 0,0, 1.6814025470759737, -1.050488399695429, - 0,-0.7853981633974483, 5.750520098164818, -1.2109809382060603) - << std::endl; -} - -quat -get_fk5_j2000_to_gal_quat() -{ - // returns the quaternion that rotates fk5j2000 to galactic J2000 - // coordinates - return quat(0.4889475076,-0.483210684,0.1962537583,0.699229742); -} - -static void -create_det_az_el_trans(const G3Timestream &az, const G3Timestream &el, - G3VectorQuat &trans_quats) // XXX: switch to G3TimestreamQuat? +static G3TimestreamQuat +get_origin_rotator_timestream(const G3Timestream &alpha, const G3Timestream &delta, + MapCoordReference coord_sys) { // Creates the transform that takes (1,0,0) to az, -el // for why it's -el see the comment at the top of this document - g3_assert(az.size() == el.size()); - trans_quats = G3VectorQuat(az.size(), quat(1,0,0,0)); - for (size_t i = 0; i < az.size(); i++) - trans_quats[i] = get_origin_rotator(az[i], -el[i]); -} + g3_assert(alpha.size() == delta.size()); + G3TimestreamQuat trans_quats(alpha.size(), quat(1,0,0,0)); + trans_quats.start = alpha.start; + trans_quats.stop = alpha.stop; + if (coord_sys == Local) + for (size_t i = 0; i < alpha.size(); i++) + trans_quats[i] = get_origin_rotator(alpha[i], -delta[i]); + else + for (size_t i = 0; i < alpha.size(); i++) + trans_quats[i] = get_origin_rotator(alpha[i], delta[i]); -static void -create_lazy_det_ra_dec_trans(const G3Timestream &ra, const G3Timestream &dec, - G3VectorQuat &trans_quats) -{ - // Creates the transform that takes (1,0,0) to ra,dec - g3_assert(ra.size() == dec.size()); - trans_quats = G3VectorQuat(ra.size(), quat(1,0,0,0)); - for (size_t i = 0; i < ra.size(); i++) - trans_quats[i] = get_origin_rotator(ra[i], dec[i]); + return trans_quats; } -static void -create_det_ra_dec_trans(const G3Timestream &az_0, const G3Timestream &el_0, +static G3TimestreamQuat +get_boresight_rotator_timestream(const G3Timestream &az_0, const G3Timestream &el_0, const G3Timestream &ra_0, const G3Timestream &dec_0, const G3Timestream &az_1, const G3Timestream &el_1, - const G3Timestream &ra_1, const G3Timestream &dec_1, - G3VectorQuat & trans_quats) + const G3Timestream &ra_1, const G3Timestream &dec_1) { // Computes the transform that takes (1,0,0) to the point (ra_0, dec_0) // and properly handles rotation about the (ra_0, dec_0) point with the @@ -336,7 +248,9 @@ create_det_ra_dec_trans(const G3Timestream &az_0, const G3Timestream &el_0, g3_assert(az_0.size() == dec_1.size()); g3_assert(az_0.size() == ra_0.size()); g3_assert(az_0.size() == ra_1.size()); - trans_quats = G3VectorQuat(ra_0.size(), quat(1,0,0,0)); + G3TimestreamQuat trans_quats(ra_0.size(), quat(1,0,0,0)); + trans_quats.start = az_0.start; + trans_quats.stop = az_0.stop; for (size_t i = 0; i < ra_0.size(); i++) { trans_quats[i] = get_transform_quat( @@ -346,19 +260,8 @@ create_det_ra_dec_trans(const G3Timestream &az_0, const G3Timestream &el_0, ra_1[i], dec_1[i] )*get_origin_rotator(az_0[i], -el_0[i]); } -} - -static void -convert_ra_dec_trans_to_gal(const G3VectorQuat &radec_trans, - G3VectorQuat &gal_trans) -{ - // Converts a rotation from (1,0,0) to fk5 j2000 into a rotation that - // takes (1,0,0) to our galactic (l,b) - gal_trans = G3VectorQuat(radec_trans.size(), quat(1,0,0,0)); - quat gt = get_fk5_j2000_to_gal_quat(); - for (size_t i = 0; i < radec_trans.size(); i++) - gal_trans[i] = gt*radec_trans[i]; + return trans_quats; } G3VectorQuat @@ -432,24 +335,9 @@ get_detector_pointing(double x_offset, double y_offset, } for (size_t i = 0; i < alpha.size(); i++) { - //using boost inverse - //quat q=trans_quat[i]*det_pos/trans_quat[i]; - //uses an inverse that assumes we are on the unit sphere - const quat & t = trans_quat[i]; - quat q=trans_quat[i]*det_pos * quat( t.R_component_1(), - -t.R_component_2(), -t.R_component_3(), -t.R_component_4()); - + quat q = trans_quat[i] * det_pos * ~trans_quat[i]; quat_to_ang(q, alpha[i], delta[i]); - - #ifdef CHECK_QUAT_INVERSE - double a,d; - quat u = trans_quat[i]*det_pos/trans_quat[i]; - quat_to_ang(u, a, d); - if( fabs(a - alpha[i]) > 1e-5 || fabs(d - delta[i]) > 1e-5){ - log_fatal("Failed trans %lf %lf %lf %lf\n", a, alpha[i], d, delta[i]); - } - #endif } if (coord_sys == Local) { for (size_t i = 0; i < delta.size(); i++) @@ -458,145 +346,52 @@ get_detector_pointing(double x_offset, double y_offset, } -void +std::vector get_detector_rotation(double x_offset, double y_offset, - const G3VectorQuat &trans_quat, std::vector &rot) + const G3VectorQuat &trans_quat) { - // Computes the polarization angle rotation that occurs under the - // transform trans_quat and stores it in rot. - - rot = std::vector(trans_quat.size(), 0); + // Computes the polarization angle rotation that occurs under the + // transform trans_quat and stores it in rot. + std::vector rot(trans_quat.size(), 0); quat det_pos = offsets_to_quat(x_offset, y_offset); - for (size_t i = 0; i < rot.size(); i++) { - quat q = trans_quat[i]*det_pos/trans_quat[i]; - rot[i] = get_rot_ang(det_pos, q, trans_quat[i]); - } -} - -static std::vector -convert_celestial_offsets_to_local_offsets(G3VectorQuat trans_vec, - double x_offset_celest, double y_offset_celest) -{ - double alpha, delta; - double x_offset_local, y_offset_local; - - // written at a time when we only have python bindings for the vector quat - g3_assert( trans_vec.size() == 1); - quat trans = trans_vec[0]; - - // First we compute the celestial position of boresight. - quat bs_cel = trans * quat(0,1,0,0) / trans; - quat_to_ang(bs_cel, alpha, delta); //alpha, delta are celestial boresight - - // Next, with our celestial position offsets we compute the position - // of our detector in celestial coordinates with the celestial offsets. - // because the offset is in celestial coordinates we do that with the transform - // that takes (1,0,0) in celestial coordinates to the boresight while ignoring - //boresight rotation - - quat celestial_trans = get_origin_rotator(alpha, delta); - - quat det_cel = celestial_trans * offsets_to_quat(x_offset_celest, y_offset_celest) / - celestial_trans; - - // We then transform those positions to local coordinates using the full transform - quat inv_trans = boost::math::conj(trans); - quat det_local = inv_trans * det_cel / inv_trans; - - #ifdef CHECK_QUAT_INVERSE - //just to check I'm not insane can be dropped in the future - g3_assert(sloppy_eq(inv_trans * bs_cel / inv_trans, quat(0,1,0,0))); - #endif - - // we read off the position of the detector. - quat_to_ang(det_local, x_offset_local, y_offset_local); - - // C++11 is magic. - return {x_offset_local, -y_offset_local}; -} - -static double -angle_d2(double ra_0, double dec_0, double ra_1, double dec_1) -{ - double delta_ra_abs = fabs(ra_0 - ra_1); - double delta_dec = dec_0 - dec_1; - if (delta_ra_abs > PI){ - delta_ra_abs = 2 * PI - delta_ra_abs; - } - delta_ra_abs /= cos(dec_0); - return delta_ra_abs*delta_ra_abs + delta_dec*delta_dec; - -} - -static G3VectorQuat -get_closest_transform(double ra, double dec, - const G3VectorDouble &ras, const G3VectorDouble &decs, - const G3VectorQuat &trans) -{ - // Returns quaternion transformation nearest the requested ra/dec - - double dist = angle_d2(ras[0]/G3Units::rad, decs[0]/G3Units::rad, - ra/G3Units::rad, dec/G3Units::rad); - size_t ind = 0; - for (size_t i = 0; i < ras.size(); i++) { - double td = angle_d2(ras[i]/G3Units::rad, decs[i]/G3Units::rad, - ra/G3Units::rad, dec/G3Units::rad); - if (td < dist) { - dist = td; - ind = i; - } - } - G3VectorQuat v; - v.push_back(trans[ind]); - return v; -} - -double -angular_distance(double alpha0, double delta0, double alpha1, double delta1) -{ - alpha0 /= G3Units::rad; - delta0 /= G3Units::rad; - alpha1 /= G3Units::rad; - delta1 /= G3Units::rad; - - double sa = sin((alpha0 - alpha1) / 2.0); - double sd = sin((delta0 - delta1) / 2.0); - double v = sd * sd + cos(delta0) * cos(delta1) * sa * sa; + for (size_t i = 0; i < rot.size(); i++) + rot[i] = get_rot_ang(det_pos, trans_quat[i]); - return 2 * asin(sqrt(v)) * G3Units::rad; + return rot; } PYBINDINGS("maps") { using namespace boost::python; - def("test_trans_", test_trans); - def("test_gal_trans_", test_gal_trans); - def("test_gal_trans_rot_", test_gal_trans_rot); - def("print_fk5_j2000_to_gal_quat_", print_fk5_j2000_to_gal_quat); + // for testing def("c_quat_to_ang_", py_quat_to_ang); def("c_ang_to_quat_", ang_to_quat); - def("get_origin_rotator", get_origin_rotator, (arg("alpha"), arg("delta")), - "Compute the transformation quaternion that would rotate the " - "vector (1, 0, 0) to point in the given direction."); - def("offsets_to_quat", offsets_to_quat, (arg("x"), arg("y")), - "Returns the vector quaternion (0,1,0,0) rotated by the given " - "x and y offsets. Equivalent to ``t * quat(0,1,0,0) / t``, where " - "``t = get_origin_rotator(x, -y)``"); - def("create_det_az_el_trans", create_det_az_el_trans, - "Construct a quaternion vector from timestreams of detector " - "azimuth and elevation. Equivalent to ``R_z(az) * R_y(-el)``."); - def("create_lazy_det_ra_dec_trans", create_lazy_det_ra_dec_trans, - "Construct a quaternion vector from timestreams of detector " - "RA and declination. Equivalent to ``R_z(ra) * R_y(dec)``"); - def("create_det_ra_dec_trans", create_det_ra_dec_trans, - "Construct a quaternion vector from timestreams of detector " - "coordinates. Computes the transformation from local " - "(az_0, el_0) coordinates to celestial (ra_0, dec_0), " - "accounting for rotation about the boresight by including " - "the second set of points."); - def("convert_ra_dec_trans_to_gal", convert_ra_dec_trans_to_gal, - "Rotate a vector of quaternions from Equatorial to Galactic " - "coordinates"); + def("get_fk5_j2000_to_gal_quat", get_fk5_j2000_to_gal_quat, + "Return the rotation quaternion to rotate from equatorial to " + "galactic coordinates."); + def("get_origin_rotator", get_origin_rotator, (arg("alpha"), arg("delta")), + "Compute the transformation quaternion that would rotate the " + "vector (1, 0, 0) to point in the given direction."); + def("offsets_to_quat", offsets_to_quat, (arg("x"), arg("y")), + "Returns the vector quaternion (0,1,0,0) rotated by the given " + "x and y offsets. Equivalent to ``t * quat(0,1,0,0) / t``, where " + "``t = get_origin_rotator(x, -y)``"); + def("get_transform_quat", get_transform_quat, + "Computes a rotation that will take (as_0,ds_0) to (ae_0, de_0) and " + "(as_1, ds_1) to (ae_1, de_1)"); + def("get_rot_ang", get_rot_ang, (arg("start_q"), arg("trans")), + "Computes the boresight rotation of the vector start_q when rotated " + "by trans."); + def("get_origin_rotator_timestream", get_origin_rotator_timestream, + (arg("alpha"), arg("delta"), arg("coord_sys")), + "Construct a transform quaternion timestream from timestreams of sky " + "coordinates. Equivalent to ``R_z(alpha) * R_y(delta)``."); + def("get_boresight_rotator_timestream", get_boresight_rotator_timestream, + "Construct a transform quaternion timestream from timestreams of " + "local and equatorial boresight pointing coordinates. Computes the " + "transform from local (az_0, el_0) coordinates to equatorial " + "(ra_0, dec_0), accounting for rotation about the boresight by " + "including the second set of points."); } diff --git a/maps/tests/transtest.py b/maps/tests/transtest.py index aef5c99c..ea66f117 100755 --- a/maps/tests/transtest.py +++ b/maps/tests/transtest.py @@ -29,75 +29,36 @@ import astropy.coordinates from astropy import units as u c = astropy.coordinates.FK5(ra=np.asarray(o_ra_0)*u.rad, dec=np.asarray(o_dec_0)*u.rad) -g = c.transform_to(astropy.coordinates.Galactic) +g = c.transform_to(astropy.coordinates.Galactic()) l_test = g.l.radian b_test = g.b.radian -def wrap_ring(a): - np.asarray(a)[np.where(a > np.pi)] -= 2*np.pi -def sloppy_eq(f,g, eps = 1e-5): - return np.abs(f-g) < eps +def sloppy_eq(f, g, eps=1e-5): + np.testing.assert_allclose(f, g, rtol=0, atol=eps) +def apply_trans(alpha, delta, trans): + q_off = maps.ang_to_quat(alpha, delta) + ra, dec = maps.quat_to_ang(trans * q_off / trans) + if ra < 0: + ra += 2 * np.pi + return ra, dec for i in range(n_samps): - t_ra_0, t_dec_0 = maps.test_trans_( az_0[i], el_0[i], ra_0[i], dec_0[i], - az_1[i], el_1[i], ra_1[i], dec_1[i], - az_0[i], el_0[i]) - if t_ra_0 <0: - t_ra_0 += 2*np.pi - if not sloppy_eq(t_ra_0, ra_0[i]): - print(t_ra_0, ra_0[i]) - print(az_0[i], el_0[i], ra_0[i], dec_0[i], - az_1[i], el_1[i], ra_1[i], dec_1[i], - az_0[i], el_0[i]) - assert(0) - if not sloppy_eq(t_dec_0, dec_0[i]): - print(t_dec_0, dec_0[i]) - print(az_0[i], el_0[i], ra_0[i], dec_0[i], - az_1[i], el_1[i], ra_1[i], dec_1[i], - az_0[i], el_0[i]) - assert(0) - - t_ra_1, t_dec_1 = maps.test_trans_( az_0[i], el_0[i], ra_0[i], dec_0[i], - az_1[i], el_1[i], ra_1[i], dec_1[i], - az_1[i], el_1[i]) - if t_ra_1 <0: - t_ra_1 += 2*np.pi - - if not sloppy_eq(t_ra_1, ra_1[i]): - print(t_ra_1, ra_1[i]) - print(az_0[i], el_0[i], ra_0[i], dec_0[i], - az_1[i], el_1[i], ra_1[i], dec_1[i]) - assert(0) - if not sloppy_eq(t_dec_1, dec_1[i]): - print(t_dec_1, dec_1[i]) - print(az_0[i], el_0[i], ra_0[i], dec_0[i], - az_1[i], el_1[i], ra_1[i], dec_1[i]) - assert(0) - - t_ra_o, t_dec_o = maps.test_trans_( az_0[i], el_0[i], ra_0[i], dec_0[i], - az_1[i], el_1[i], ra_1[i], dec_1[i], - o_az_0[i], o_el_0[i]) - if t_ra_o <0: - t_ra_o += 2*np.pi - if not sloppy_eq(t_ra_o, o_ra_0[i]): - print(t_ra_o, o_ra_0[i]) - assert(0) - if not sloppy_eq(t_dec_o, o_dec_0[i]): - print(t_dec_o, o_dec_0[i]) - assert(0) - - - - t_l_o, t_b_o = maps.test_gal_trans_( az_0[i], el_0[i], ra_0[i], dec_0[i], - az_1[i], el_1[i], ra_1[i], dec_1[i], - o_az_0[i], o_el_0[i]) - - if t_l_o < 0: - t_l_o += 2*np.pi - if not sloppy_eq(t_l_o, l_test[i]): - assert(0) - if not sloppy_eq(t_b_o, b_test[i]): - assert(0) + q_trans = maps.get_transform_quat( + az_0[i], -el_0[i], ra_0[i], dec_0[i], + az_1[i], -el_1[i], ra_1[i], dec_1[i] + ); + t_ra_0, t_dec_0 = apply_trans(az_0[i], -el_0[i], q_trans) + sloppy_eq([t_ra_0, t_dec_0], [ra_0[i], dec_0[i]]) + + t_ra_1, t_dec_1 = apply_trans(az_1[i], -el_1[i], q_trans) + sloppy_eq([t_ra_1, t_dec_1], [ra_1[i], dec_1[i]]) + + t_ra_o, t_dec_o = apply_trans(o_az_0[i], -o_el_0[i], q_trans) + sloppy_eq([t_ra_o, t_dec_o], [o_ra_0[i], o_dec_0[i]]) + + q_trans_gal = maps.get_fk5_j2000_to_gal_quat() * q_trans + t_l_o, t_b_o = apply_trans(o_az_0[i], -o_el_0[i], q_trans_gal) + sloppy_eq([t_l_o, t_b_o], [l_test[i], b_test[i]])