Packages/MoonriseKit/Sources/SkyEngine/SkyEngine.swift

  1import CAstronomy
  2import Foundation
  3
  4/// A place on Earth. Latitude and longitude in degrees (east positive), height in metres.
  5public struct Observer: Sendable, Equatable {
  6    public var latitude: Double
  7    public var longitude: Double
  8    public var height: Double
  9
 10    public init(latitude: Double, longitude: Double, height: Double = 0) {
 11        self.latitude = latitude
 12        self.longitude = longitude
 13        self.height = height
 14    }
 15
 16    var astro: astro_observer_t { Astronomy_MakeObserver(latitude, longitude, height) }
 17}
 18
 19public enum Body: String, CaseIterable, Sendable {
 20    case sun, moon, mercury, venus, mars, jupiter, saturn, uranus, neptune, pluto
 21
 22    public static let planets: [Body] = [.mercury, .venus, .mars, .jupiter, .saturn, .uranus, .neptune]
 23
 24    var astro: astro_body_t {
 25        switch self {
 26        case .sun: BODY_SUN
 27        case .moon: BODY_MOON
 28        case .mercury: BODY_MERCURY
 29        case .venus: BODY_VENUS
 30        case .mars: BODY_MARS
 31        case .jupiter: BODY_JUPITER
 32        case .saturn: BODY_SATURN
 33        case .uranus: BODY_URANUS
 34        case .neptune: BODY_NEPTUNE
 35        case .pluto: BODY_PLUTO
 36        }
 37    }
 38}
 39
 40/// Apparent right ascension (hours) and declination (degrees) of date, distance in AU.
 41public struct Equatorial: Sendable, Equatable {
 42    public var rightAscension: Double
 43    public var declination: Double
 44    public var distance: Double
 45}
 46
 47/// Altitude above the horizon and azimuth from north through east, both in degrees.
 48public struct Horizontal: Sendable, Equatable {
 49    public var altitude: Double
 50    public var azimuth: Double
 51
 52    public init(altitude: Double, azimuth: Double) {
 53        self.altitude = altitude
 54        self.azimuth = azimuth
 55    }
 56}
 57
 58public struct Illumination: Sendable, Equatable {
 59    public var magnitude: Double
 60    /// Fraction of the disc that is lit, 0...1.
 61    public var phaseFraction: Double
 62}
 63
 64public enum SkyEngineError: Error, Equatable {
 65    case astronomy(status: Int)
 66}
 67
 68/// Positions and phases from Astronomy Engine. Its C code keeps internal caches,
 69/// so every call into it holds one process-wide lock.
 70public enum SkyEngine {
 71    static let lock = NSRecursiveLock()
 72
 73    static func time(_ date: Date) -> astro_time_t {
 74        // Days since J2000 (2000-01-01 12:00 UTC).
 75        Astronomy_TimeFromDays((date.timeIntervalSince1970 - 946_728_000) / 86_400)
 76    }
 77
 78    static func check(_ status: astro_status_t) throws {
 79        guard status == ASTRO_SUCCESS else { throw SkyEngineError.astronomy(status: Int(status.rawValue)) }
 80    }
 81
 82    public static func equatorial(_ body: Body, at date: Date, from observer: Observer) throws -> Equatorial {
 83        lock.lock()
 84        defer { lock.unlock() }
 85        var t = time(date)
 86        let eq = Astronomy_Equator(body.astro, &t, observer.astro, EQUATOR_OF_DATE, ABERRATION)
 87        try check(eq.status)
 88        return Equatorial(rightAscension: eq.ra, declination: eq.dec, distance: eq.dist)
 89    }
 90
 91    public static func horizontal(_ body: Body, at date: Date, from observer: Observer, refraction: Bool = false) throws -> Horizontal {
 92        let eq = try equatorial(body, at: date, from: observer)
 93        return horizontal(rightAscension: eq.rightAscension, declination: eq.declination, at: date, from: observer, refraction: refraction)
 94    }
 95
 96    /// For fixed objects such as stars: RA in hours, declination in degrees, both of date.
 97    public static func horizontal(rightAscension: Double, declination: Double, at date: Date, from observer: Observer, refraction: Bool = false) -> Horizontal {
 98        lock.lock()
 99        defer { lock.unlock() }
100        var t = time(date)
101        let hor = Astronomy_Horizon(&t, observer.astro, rightAscension, declination, refraction ? REFRACTION_NORMAL : REFRACTION_NONE)
102        return Horizontal(altitude: hor.altitude, azimuth: hor.azimuth)
103    }
104
105    /// Horizontal positions for many fixed objects at once, from J2000 RA and declination in degrees.
106    /// Applies precession and nutation; no aberration or refraction.
107    /// 2000-01-01 12:00 TT.
108    public static let j2000 = Date(timeIntervalSince1970: 946_727_935.816)
109
110    /// Right ascension and declination on the equator and equinox of `date` (precession and nutation applied), degrees.
111    public static func equatorOfDate(j2000 c: (ra: Double, dec: Double), at date: Date) -> (ra: Double, dec: Double) {
112        lock.lock()
113        defer { lock.unlock() }
114        var t = time(date)
115        let rotation = Astronomy_Rotation_EQJ_EQD(&t)
116        let ra = c.ra * .pi / 180, dec = c.dec * .pi / 180
117        let v = Astronomy_RotateVector(rotation, astro_vector_t(status: ASTRO_SUCCESS, x: cos(dec) * cos(ra), y: cos(dec) * sin(ra), z: sin(dec), t: t))
118        var r = atan2(v.y, v.x) * 180 / .pi
119        if r < 0 { r += 360 }
120        return (r, asin(v.z) * 180 / .pi)
121    }
122
123    public static func horizontal(j2000 coordinates: [(ra: Double, dec: Double)], at date: Date, from observer: Observer) -> [Horizontal] {
124        lock.lock()
125        defer { lock.unlock() }
126        var t = time(date)
127        let rotation = Astronomy_Rotation_EQJ_HOR(&t, observer.astro)
128        return coordinates.map { c in
129            let ra = c.ra * .pi / 180, dec = c.dec * .pi / 180
130            let v = astro_vector_t(status: ASTRO_SUCCESS, x: cos(dec) * cos(ra), y: cos(dec) * sin(ra), z: sin(dec), t: t)
131            let s = Astronomy_HorizonFromVector(Astronomy_RotateVector(rotation, v), REFRACTION_NONE)
132            return Horizontal(altitude: s.lat, azimuth: s.lon)
133        }
134    }
135
136    /// Heliocentric position in AU, J2000 equatorial axes.
137    public static func heliocentric(_ body: Body, at date: Date) throws -> SIMD3<Double> {
138        lock.lock()
139        defer { lock.unlock() }
140        let v = Astronomy_HelioVector(body.astro, time(date))
141        try check(v.status)
142        return SIMD3(v.x, v.y, v.z)
143    }
144
145    /// Heliocentric position in AU on J2000 ecliptic axes (x toward the equinox, z toward the ecliptic north pole).
146    public static func heliocentricEcliptic(_ body: Body, at date: Date) throws -> SIMD3<Double> {
147        ecliptic(try heliocentric(body, at: date))
148    }
149
150    /// Earth's heliocentric position in AU on J2000 ecliptic axes.
151    public static func heliocentricEarth(at date: Date) throws -> SIMD3<Double> {
152        lock.lock()
153        defer { lock.unlock() }
154        let v = Astronomy_HelioVector(BODY_EARTH, time(date))
155        try check(v.status)
156        return ecliptic(SIMD3(v.x, v.y, v.z))
157    }
158
159    /// The Moon's geocentric position in km on J2000 ecliptic axes (as JPL Horizons gives with REF_PLANE ECLIPTIC).
160    public static func geocentricMoonEcliptic(at date: Date) -> SIMD3<Double> {
161        lock.lock()
162        defer { lock.unlock() }
163        let v = Astronomy_GeoMoon(time(date))
164        return ecliptic(SIMD3(v.x, v.y, v.z)) * 1.495978707e8
165    }
166
167    /// J2000 equatorial to J2000 ecliptic axes.
168    static func ecliptic(_ v: SIMD3<Double>) -> SIMD3<Double> {
169        let e = 23.4392911 * Double.pi / 180
170        return SIMD3(v.x, v.y * cos(e) + v.z * sin(e), -v.y * sin(e) + v.z * cos(e))
171    }
172
173    public static func illumination(_ body: Body, at date: Date) throws -> Illumination {
174        lock.lock()
175        defer { lock.unlock() }
176        let i = Astronomy_Illumination(body.astro, time(date))
177        try check(i.status)
178        return Illumination(magnitude: i.mag, phaseFraction: i.phase_fraction)
179    }
180
181    static func date(_ t: astro_time_t) -> Date {
182        Date(timeIntervalSince1970: 946_728_000 + t.ut * 86_400)
183    }
184
185    /// Next rise and next set within `limitDays` after `date`; nil if none in that window.
186    public static func riseSet(_ body: Body, after date: Date, from observer: Observer, limitDays: Double = 1) -> (rise: Date?, set: Date?) {
187        lock.lock()
188        defer { lock.unlock() }
189        let t = time(date)
190        let rise = Astronomy_SearchRiseSetEx(body.astro, observer.astro, DIRECTION_RISE, t, limitDays, 0)
191        let set = Astronomy_SearchRiseSetEx(body.astro, observer.astro, DIRECTION_SET, t, limitDays, 0)
192        return (rise.status == ASTRO_SUCCESS ? Self.date(rise.time) : nil,
193                set.status == ASTRO_SUCCESS ? Self.date(set.time) : nil)
194    }
195
196    public enum MoonQuarter: Int, Sendable {
197        case new, firstQuarter, full, lastQuarter
198    }
199
200    /// The first lunar quarter after `date`.
201    public static func nextMoonQuarter(after date: Date) throws -> (quarter: MoonQuarter, date: Date) {
202        lock.lock()
203        defer { lock.unlock() }
204        let q = Astronomy_SearchMoonQuarter(time(date))
205        try check(q.status)
206        return (MoonQuarter(rawValue: Int(q.quarter))!, Self.date(q.time))
207    }
208
209    /// Moon's ecliptic longitude east of the sun in degrees: 0 new, 90 first quarter, 180 full, 270 last quarter.
210    public static func moonPhase(at date: Date) throws -> Double {
211        lock.lock()
212        defer { lock.unlock() }
213        let a = Astronomy_MoonPhase(time(date))
214        try check(a.status)
215        return a.angle
216    }
217}