matchRaDec: fix imprecision at small distance

The formula used for converting unit sphere distance^2 to an Angle
isn't precise for small distance. Using an alternative formula that
(unfortunately) includes a sqrt is more precise.

Bug identified by Dominique Boutigny.

src/table/Match.cc

@@ -137,8 +137,8 @@ double toUnitSphereDistanceSquared(afw::geom::Angle theta) noexcept {
  * @returns the angle between the two vectors
  */
 Angle fromUnitSphereDistanceSquared(double d2) noexcept {
-    return (std::acos(1. - d2 / 2.)) * afw::geom::radians;
-    // == 2.0 * asin(0.5 * sqrt(d2))
+    // acos(1 - 0.5*d2) doesn't require sqrt but isn't as precise for small d2
+    return 2.0*std::asin(0.5*std::sqrt(d2))*afw::geom::radians;
 }
 
 }  // anonymous

tests/test_sourceMatch.py

@@ -350,6 +350,41 @@ def assertEqualFloat(self, value1, value2):
         self.assertTrue(value1 == value2 or
                         (np.isnan(value1) and np.isnan(value2)))
 
+    def testDistancePrecision(self):
+        """Test for precision of the calculated distance
+
+        Check that the distance produced by matchRaDec is the same
+        as the distance produced from calculating the separation
+        between the matched coordinates.
+
+        Based on DM-13891.
+        """
+        rng = np.random.RandomState(12345)  # I have the same combination on my luggage
+        radius = 0.5*afwGeom.arcseconds
+        num = 1000
+        raKey = afwTable.SourceTable.getCoordKey().getRa()
+        decKey = afwTable.SourceTable.getCoordKey().getDec()
+        for ii in range(num):
+            src1 = self.ss1.addNew()
+            src1.setId(ii)
+            src1.set(afwTable.SourceTable.getCoordKey().getRa(),
+                     (10 + 0.001*ii) * afwGeom.degrees)
+            src1.set(afwTable.SourceTable.getCoordKey().getDec(),
+                     (10 + 0.001*ii) * afwGeom.degrees)
+
+            src2 = self.ss2.addNew()
+            src2.setId(2*num + ii)
+            coord = src1.getCoord().offset(rng.uniform(high=360)*afwGeom.degrees,
+                                           rng.uniform(high=radius.asArcseconds())*afwGeom.arcseconds)
+            src2.set(afwTable.SourceTable.getCoordKey().getRa(), coord.getLongitude())
+            src2.set(afwTable.SourceTable.getCoordKey().getDec(), coord.getLatitude())
+
+        matches = afwTable.matchRaDec(self.ss1, self.ss2, radius)
+        dist1 = np.array([(mm.distance*afwGeom.radians).asArcseconds() for mm in matches])
+        dist2 = np.array([mm.first.getCoord().separation(mm.second.getCoord()).asArcseconds() for mm in matches])
+        diff = dist1 - dist2
+        self.assertLess(diff.std(), 1.0e-10)  # I get 4e-12
+
 
 class TestMemory(lsst.utils.tests.MemoryTestCase):
     pass