Packages/MoonriseKit/Tests/SkyEngineTests/HorizonsTests.swift

 1import Foundation
 2import Testing
 3@testable import SkyEngine
 4
 5// Reference values from JPL Horizons (ssd.jpl.nasa.gov/api/horizons.api), fetched 2026-10-17.
 6// Observer table: site -87.6298 E, 41.8781 N, 0 km; 2026-10-17 02:00 UT; quantities 2 (apparent
 7// RA/Dec, degrees) and 4 (airless apparent azimuth/elevation).
 8// Vector table: CENTER 500@10, REF_PLANE FRAME, 2026-10-17 00:00 UT, AU.
 9
10let chicago = Observer(latitude: 41.8781, longitude: -87.6298)
11let observed = Date(timeIntervalSince1970: 1_792_202_400) // 2026-10-17T02:00:00Z
12let vectorEpoch = Date(timeIntervalSince1970: 1_792_195_200) // 2026-10-17T00:00:00Z
13
14struct HorizonsRow {
15    var body: Body
16    var raDegrees: Double
17    var dec: Double
18    var azimuth: Double
19    var altitude: Double
20}
21
22let observerRows: [HorizonsRow] = [
23    HorizonsRow(body: .moon, raDegrees: 276.63979, dec: -28.01156, azimuth: 223.832359, altitude: 5.580273),
24    HorizonsRow(body: .mars, raDegrees: 133.60867, dec: 18.78232, azimuth: 15.406574, altitude: -27.902306),
25    HorizonsRow(body: .jupiter, raDegrees: 144.84484, dec: 14.69461, azimuth: 3.616329, altitude: -33.353920),
26    HorizonsRow(body: .saturn, raDegrees: 10.55474, dec: 1.59612, azimuth: 124.822854, altitude: 34.510943),
27]
28
29/// One arcminute, Astronomy Engine's stated accuracy.
30let arcminute = 1.0 / 60
31
32func angularDifference(_ a: Double, _ b: Double) -> Double {
33    abs(remainder(a - b, 360))
34}
35
36@Test(arguments: observerRows)
37func apparentEquatorialMatchesHorizons(_ row: HorizonsRow) throws {
38    let eq = try SkyEngine.equatorial(row.body, at: observed, from: chicago)
39    let raError = angularDifference(eq.rightAscension * 15, row.raDegrees) * cos(row.dec * .pi / 180)
40    #expect(raError < arcminute, "\(row.body) RA off by \(raError * 60)'")
41    #expect(abs(eq.declination - row.dec) < arcminute, "\(row.body) Dec off by \(abs(eq.declination - row.dec) * 60)'")
42}
43
44@Test(arguments: observerRows)
45func horizontalMatchesHorizons(_ row: HorizonsRow) throws {
46    let hor = try SkyEngine.horizontal(row.body, at: observed, from: chicago)
47    let azError = angularDifference(hor.azimuth, row.azimuth) * cos(row.altitude * .pi / 180)
48    #expect(azError < arcminute, "\(row.body) azimuth off by \(azError * 60)'")
49    #expect(abs(hor.altitude - row.altitude) < arcminute, "\(row.body) altitude off by \(abs(hor.altitude - row.altitude) * 60)'")
50}
51
52@Test(arguments: [
53    (Body.jupiter, SIMD3(-3.581981680856141, 3.573182031399209, 1.618755638647523)),
54    (Body.saturn, SIMD3(9.237054585455835, 1.860186986025053, 0.3705349659085789)),
55])
56func heliocentricMatchesHorizons(_ body: Body, _ expected: SIMD3<Double>) throws {
57    let v = try SkyEngine.heliocentric(body, at: vectorEpoch)
58    let d = v - expected
59    let error = (d * d).sum().squareRoot()
60    let angle = error / (expected * expected).sum().squareRoot() * 180 / .pi
61    #expect(angle < arcminute, "\(body) off by \(error) AU (\(angle * 3600)\")")
62}
63
64@Test func moonPhaseIsWaxing() throws {
65    let phase = try SkyEngine.moonPhase(at: observed)
66    #expect(phase > 0 && phase < 180)
67}