Bug fix in FlatSkyMap.query_disc method

The method for determining the box of valid pixels that contain the disc to
query used an incorrect calculation for offsetting the disc center to sweep out
a circle.  This PR fixes this calculation by using a slightly less efficient but
reliable method for computing the offset quaternion.

Author: arahlin

maps/src/FlatSkyProjection.cxx

@@ -602,8 +602,9 @@ FlatSkyProjection::QueryDisc(const Quat &q, double radius) const
 {
 	static const size_t npts = 72;
 	double dd = -2.0 * radius / sqrt(2.0);
-	auto qd = get_origin_rotator(0, dd);
-	auto p = qd * q * ~qd;
+	double alpha, delta;
+	quat_to_ang(q, alpha, delta);
+	auto p = ang_to_quat(alpha, delta + (((90*deg - delta) < -dd) ? dd : -dd));
 	double pva = q.b();
 	double pvb = q.c();
 	double pvc = q.d();

maps/src/G3SkyMap.cxx

@@ -436,9 +436,10 @@ G3SkyMap::QueryAlphaEllipse(const Quat &q, double a, double b) const
 	double da = ACOS(COS(rmaj) / COS(rmin)) / cd;
 
 	// focus locations
-	auto qda = get_origin_rotator(da, 0);
-	auto ql = qda * q * ~qda;
-	auto qr = ~qda * q * qda;
+	double alpha, delta;
+	quat_to_ang(q, alpha, delta);
+	auto ql = ang_to_quat(alpha - da, delta);
+	auto qr = ang_to_quat(alpha + da, delta);
 
 	// narrow search to pixels within the major disc
 	auto disc = QueryDisc(q, rmaj);

maps/tests/flatsky_maps.py

@@ -109,27 +109,35 @@
 ).ravel()
 radius = 5 * x.res
 
-avec = []
-dvec = []
-rvec = []
-masked = set()
-
-for a in numpy.linspace(numpy.min(alpha), numpy.max(alpha), 20):
-    for d in numpy.linspace(numpy.min(delta), numpy.max(delta), 20):
-        avec.append(a)
-        dvec.append(d)
-        rvec.append(radius)
-        pix1 = x.query_disc(a, d, radius)
-        pix2 = numpy.where(
-            angles.separation(
-                SkyCoord(
-                    a / core.G3Units.deg * u.degree,
-                    d / core.G3Units.deg * u.degree,
-                )
-            ).deg * core.G3Units.deg < radius
-        )[0]
-        masked |= set(pix2)
-        assert not (set(pix1) ^ set(pix2))
-
-mask = maps.make_point_source_mask(x, numpy.asarray(avec), numpy.asarray(dvec), numpy.asarray(rvec))
-assert not (set(mask.nonzero()) ^ masked)
+for alpha1 in numpy.linspace(0, 360, 13) * core.G3Units.deg:
+    for delta1 in numpy.linspace(-90, 90, 5) * core.G3Units.deg:
+        avec = []
+        dvec = []
+        rvec = []
+        masked = set()
+
+        x = FlatSkyMap(50, 50, core.G3Units.arcmin, proj=MapProjection.ProjZEA, alpha_center=alpha1, delta_center=delta1)
+        alpha, delta = maps.get_ra_dec_map(x)
+        angles = SkyCoord(
+            numpy.asarray(alpha) / core.G3Units.deg * u.degree,
+            numpy.asarray(delta) / core.G3Units.deg * u.degree,
+        ).ravel()
+        for a in numpy.linspace(numpy.min(alpha), numpy.max(alpha), 20):
+            for d in numpy.linspace(numpy.min(delta), numpy.max(delta), 20):
+                avec.append(a)
+                dvec.append(d)
+                rvec.append(radius)
+                pix1 = x.query_disc(a, d, radius)
+                pix2 = numpy.where(
+                    angles.separation(
+                        SkyCoord(
+                            a / core.G3Units.deg * u.degree,
+                            d / core.G3Units.deg * u.degree,
+                        )
+                    ).deg * core.G3Units.deg < radius
+                )[0]
+                masked |= set(pix2)
+                assert not (set(pix1) ^ set(pix2))
+
+        mask = maps.make_point_source_mask(x, numpy.asarray(avec), numpy.asarray(dvec), numpy.asarray(rvec))
+        assert not (set(mask.nonzero()) ^ masked)