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        DeltaT.installed
 75        // Days since J2000 (2000-01-01 12:00 UTC).
 76        return Astronomy_TimeFromDays((date.timeIntervalSince1970 - 946_728_000) / 86_400)
 77    }
 78
 79    static func check(_ status: astro_status_t) throws {
 80        guard status == ASTRO_SUCCESS else { throw SkyEngineError.astronomy(status: Int(status.rawValue)) }
 81    }
 82
 83    public static func equatorial(_ body: Body, at date: Date, from observer: Observer) throws -> Equatorial {
 84        lock.lock()
 85        defer { lock.unlock() }
 86        var t = time(date)
 87        let eq = Astronomy_Equator(body.astro, &t, observer.astro, EQUATOR_OF_DATE, ABERRATION)
 88        try check(eq.status)
 89        return Equatorial(rightAscension: eq.ra, declination: eq.dec, distance: eq.dist)
 90    }
 91
 92    public static func horizontal(_ body: Body, at date: Date, from observer: Observer, refraction: Bool = false) throws -> Horizontal {
 93        let eq = try equatorial(body, at: date, from: observer)
 94        return horizontal(rightAscension: eq.rightAscension, declination: eq.declination, at: date, from: observer, refraction: refraction)
 95    }
 96
 97    /// For fixed objects such as stars: RA in hours, declination in degrees, both of date.
 98    public static func horizontal(rightAscension: Double, declination: Double, at date: Date, from observer: Observer, refraction: Bool = false) -> Horizontal {
 99        lock.lock()
100        defer { lock.unlock() }
101        var t = time(date)
102        let hor = Astronomy_Horizon(&t, observer.astro, rightAscension, declination, refraction ? REFRACTION_NORMAL : REFRACTION_NONE)
103        return Horizontal(altitude: hor.altitude, azimuth: hor.azimuth)
104    }
105
106    /// Horizontal positions for many fixed objects at once, from J2000 RA and declination in degrees.
107    /// Applies precession and nutation; no aberration or refraction.
108    /// 2000-01-01 12:00 TT.
109    public static let j2000 = Date(timeIntervalSince1970: 946_727_935.816)
110
111    /// Right ascension and declination on the equator and equinox of `date` (precession and nutation applied), degrees.
112    public static func equatorOfDate(j2000 c: (ra: Double, dec: Double), at date: Date) -> (ra: Double, dec: Double) {
113        lock.lock()
114        defer { lock.unlock() }
115        var t = time(date)
116        let rotation = Astronomy_Rotation_EQJ_EQD(&t)
117        let ra = c.ra * .pi / 180, dec = c.dec * .pi / 180
118        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))
119        var r = atan2(v.y, v.x) * 180 / .pi
120        if r < 0 { r += 360 }
121        return (r, asin(v.z) * 180 / .pi)
122    }
123
124    public static func horizontal(j2000 coordinates: [(ra: Double, dec: Double)], at date: Date, from observer: Observer) -> [Horizontal] {
125        lock.lock()
126        defer { lock.unlock() }
127        var t = time(date)
128        let rotation = Astronomy_Rotation_EQJ_HOR(&t, observer.astro)
129        return coordinates.map { c in
130            let ra = c.ra * .pi / 180, dec = c.dec * .pi / 180
131            let v = astro_vector_t(status: ASTRO_SUCCESS, x: cos(dec) * cos(ra), y: cos(dec) * sin(ra), z: sin(dec), t: t)
132            let s = Astronomy_HorizonFromVector(Astronomy_RotateVector(rotation, v), REFRACTION_NONE)
133            return Horizontal(altitude: s.lat, azimuth: s.lon)
134        }
135    }
136
137    /// Heliocentric position in AU, J2000 equatorial axes.
138    public static func heliocentric(_ body: Body, at date: Date) throws -> SIMD3<Double> {
139        lock.lock()
140        defer { lock.unlock() }
141        let v = Astronomy_HelioVector(body.astro, time(date))
142        try check(v.status)
143        return SIMD3(v.x, v.y, v.z)
144    }
145
146    /// Heliocentric position in AU on J2000 ecliptic axes (x toward the equinox, z toward the ecliptic north pole).
147    public static func heliocentricEcliptic(_ body: Body, at date: Date) throws -> SIMD3<Double> {
148        ecliptic(try heliocentric(body, at: date))
149    }
150
151    /// Earth's heliocentric position in AU on J2000 ecliptic axes.
152    public static func heliocentricEarth(at date: Date) throws -> SIMD3<Double> {
153        lock.lock()
154        defer { lock.unlock() }
155        let v = Astronomy_HelioVector(BODY_EARTH, time(date))
156        try check(v.status)
157        return ecliptic(SIMD3(v.x, v.y, v.z))
158    }
159
160    /// The Moon's geocentric position in km on J2000 ecliptic axes (as JPL Horizons gives with REF_PLANE ECLIPTIC).
161    public static func geocentricMoonEcliptic(at date: Date) -> SIMD3<Double> {
162        lock.lock()
163        defer { lock.unlock() }
164        let v = Astronomy_GeoMoon(time(date))
165        return ecliptic(SIMD3(v.x, v.y, v.z)) * 1.495978707e8
166    }
167
168    /// J2000 equatorial to J2000 ecliptic axes.
169    static func ecliptic(_ v: SIMD3<Double>) -> SIMD3<Double> {
170        let e = 23.4392911 * Double.pi / 180
171        return SIMD3(v.x, v.y * cos(e) + v.z * sin(e), -v.y * sin(e) + v.z * cos(e))
172    }
173
174    public static func illumination(_ body: Body, at date: Date) throws -> Illumination {
175        lock.lock()
176        defer { lock.unlock() }
177        let i = Astronomy_Illumination(body.astro, time(date))
178        try check(i.status)
179        return Illumination(magnitude: i.mag, phaseFraction: i.phase_fraction)
180    }
181
182    static func date(_ t: astro_time_t) -> Date {
183        Date(timeIntervalSince1970: 946_728_000 + t.ut * 86_400)
184    }
185
186    /// Next rise and next set within `limitDays` after `date`; nil if none in that window.
187    public static func riseSet(_ body: Body, after date: Date, from observer: Observer, limitDays: Double = 1) -> (rise: Date?, set: Date?) {
188        lock.lock()
189        defer { lock.unlock() }
190        let t = time(date)
191        let rise = Astronomy_SearchRiseSetEx(body.astro, observer.astro, DIRECTION_RISE, t, limitDays, 0)
192        let set = Astronomy_SearchRiseSetEx(body.astro, observer.astro, DIRECTION_SET, t, limitDays, 0)
193        return (rise.status == ASTRO_SUCCESS ? Self.date(rise.time) : nil,
194                set.status == ASTRO_SUCCESS ? Self.date(set.time) : nil)
195    }
196
197    public enum MoonQuarter: Int, Sendable {
198        case new, firstQuarter, full, lastQuarter
199    }
200
201    /// The first lunar quarter after `date`.
202    public static func nextMoonQuarter(after date: Date) throws -> (quarter: MoonQuarter, date: Date) {
203        lock.lock()
204        defer { lock.unlock() }
205        let q = Astronomy_SearchMoonQuarter(time(date))
206        try check(q.status)
207        return (MoonQuarter(rawValue: Int(q.quarter))!, Self.date(q.time))
208    }
209
210    /// Moon's ecliptic longitude east of the sun in degrees: 0 new, 90 first quarter, 180 full, 270 last quarter.
211    public static func moonPhase(at date: Date) throws -> Double {
212        lock.lock()
213        defer { lock.unlock() }
214        let a = Astronomy_MoonPhase(time(date))
215        try check(a.status)
216        return a.angle
217    }
218}