diff --git a/src/astrometry.ts b/src/astrometry.ts index dcfe966..6243db9 100644 --- a/src/astrometry.ts +++ b/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 + ) + ) } /*****************************************************************************************************************/ diff --git a/tests/astrometry.spec.ts b/tests/astrometry.spec.ts index 99cfe64..08262a5 100644 --- a/tests/astrometry.spec.ts +++ b/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) + } + }) }) /*****************************************************************************************************************/