fix: resolve the angular separation of coincident coordinates in astrometry module in @observerly/astrometry

src/astrometry.ts

@@ -49,21 +49,34 @@ import {
  *
  */
 export const getAngularSeparation = (A: SphericalCoordinate, B: SphericalCoordinate): number => {
-  // Calculate the angular separation between A and B (in degrees):
-  let θ =
-    degrees(
-      Math.acos(
-        Math.sin(radians(A.θ)) * Math.sin(radians(B.θ)) +
-          Math.cos(radians(A.θ)) * Math.cos(radians(B.θ)) * Math.cos(radians(A.φ - B.φ))
-      )
-    ) % 360
-
-  // Correct for negative angles:
-  if (θ < 0) {
-    θ += 360
+  // The unit vector of A, in the frame the two coordinates are both given in:
+  const a = {
+    x: Math.cos(radians(A.θ)) * Math.cos(radians(A.φ)),
+    y: Math.cos(radians(A.θ)) * Math.sin(radians(A.φ)),
+    z: Math.sin(radians(A.θ))
   }
 
-  return θ
+  // The unit vector of B, in the same frame:
+  const b = {
+    x: Math.cos(radians(B.θ)) * Math.cos(radians(B.φ)),
+    y: Math.cos(radians(B.θ)) * Math.sin(radians(B.φ)),
+    z: Math.sin(radians(B.θ))
+  }
+
+  // The separation is taken from the magnitude of the cross product of the two unit vectors, e.g.,
+  // its sine, against their dot product, e.g., its cosine, and not as the arc cosine of the dot
+  // product alone, which is out of domain for two coincident coordinates, e.g., the rounding of
+  // the dot product carries it just beyond one and the separation is not a number, and which is
+  // ill-conditioned for a separation near zero, e.g., a milliarcsecond resolves as none at all.
+  //
+  // N.B. The separation is the angle between two directions, and so it is bounded to [0°, 180°]
+  // by construction, and is not normalised against a whole turn:
+  return degrees(
+    Math.atan2(
+      Math.hypot(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x),
+      a.x * b.x + a.y * b.y + a.z * b.z
+    )
+  )
 }
 
 /*****************************************************************************************************************/

tests/astrometry.spec.ts

@@ -68,7 +68,7 @@ describe('getAngularSeparation', () => {
         φ: spica.ra
       }
     )
-    expect(θ).toBe(32.79290589269233)
+    expect(θ).toBe(32.79290589269235)
   })
 
   it('should return the greatest possible angular separation for two objects at their antipodes', () => {
@@ -98,6 +98,36 @@ describe('getAngularSeparation', () => {
     )
     expect(θ).toBe(180)
   })
+
+  it('should return no separation for a coordinate compared with itself', () => {
+    // The dot product of a unit vector with itself rounds to just beyond one, and so the arc
+    // cosine of it alone is out of domain, e.g., the separation is not a number:
+    for (let θ = -90; θ <= 90; θ += 0.25) {
+      for (let φ = 0; φ < 360; φ += 7) {
+        expect(getAngularSeparation({ θ, φ }, { θ, φ })).toBe(0)
+      }
+    }
+  })
+
+  it('should return the greatest separation for two coordinates at their antipodes, wherever they lie', () => {
+    for (let θ = -90; θ <= 90; θ += 0.5) {
+      for (let φ = 0; φ < 360; φ += 7) {
+        expect(getAngularSeparation({ θ, φ }, { θ: -θ, φ: φ + 180 })).toBeCloseTo(180, 9)
+      }
+    }
+  })
+
+  it('should resolve a separation that is small against the coordinates themselves', () => {
+    // The separation is ill-conditioned towards zero where it is taken as the arc cosine of the
+    // dot product, e.g., a separation of a milliarcsecond is resolved as none at all:
+    const target = { θ: 7.4070639, φ: 88.7929583 }
+
+    for (const arcsec of [1e-3, 1e-2, 1e-1, 1, 10]) {
+      const separation = getAngularSeparation(target, { ...target, θ: target.θ + arcsec / 3600 })
+
+      expect(separation * 3600).toBeCloseTo(arcsec, 9)
+    }
+  })
 })
 
 /*****************************************************************************************************************/