Packages/MoonriseKit/Sources/SkyEngine/EclipsePath.swift

  1import CAstronomy
  2import Foundation
  3
  4/// Ground track of the Moon's umbra (total) or antumbra (annular) for one solar eclipse,
  5/// from the Sun–Moon shadow axis intersected with the Earth's ellipsoid.
  6public struct EclipsePath: Sendable {
  7    public struct Point: Sendable, Equatable {
  8        public var date: Date
  9        public var center: (lat: Double, lon: Double)
 10        public var north: (lat: Double, lon: Double)?
 11        public var south: (lat: Double, lon: Double)?
 12        /// Width of the path at this point, km.
 13        public var widthKm: Double
 14
 15        public static func == (a: Point, b: Point) -> Bool { a.date == b.date }
 16    }
 17
 18    public var kind: String
 19    public var peak: Date
 20    /// Greatest eclipse.
 21    public var greatest: (lat: Double, lon: Double)
 22    /// Duration of totality or annularity at greatest eclipse, seconds.
 23    public var centralDuration: TimeInterval?
 24    public var points: [Point]
 25
 26    static let sunRadiusKm = 695_700.0
 27    static let moonRadiusKm = 1_737.4
 28    static let earthA = 6_378.137
 29    static let earthF = 1 / 298.257223563
 30    static let kmPerAU = 1.495978707e8
 31
 32    /// The path of the first total or annular eclipse whose peak is within a day of `date`.
 33    public static func compute(near date: Date, step: TimeInterval = 60) -> EclipsePath? {
 34        SkyEngine.lock.lock()
 35        defer { SkyEngine.lock.unlock() }
 36        let g = Astronomy_SearchGlobalSolarEclipse(SkyEngine.time(date.addingTimeInterval(-86_400)))
 37        guard g.status == ASTRO_SUCCESS, g.kind == ECLIPSE_TOTAL || g.kind == ECLIPSE_ANNULAR else { return nil }
 38        let peak = SkyEngine.date(g.peak)
 39        guard abs(peak.timeIntervalSince(date)) < 86_400 else { return nil }
 40
 41        var points: [Point] = []
 42        var t = peak.addingTimeInterval(-4 * 3600)
 43        while t <= peak.addingTimeInterval(4 * 3600) {
 44            if let p = point(at: t, step: step) { points.append(p) }
 45            t = t.addingTimeInterval(step)
 46        }
 47
 48        let at = Observer(latitude: g.latitude, longitude: g.longitude)
 49        let local = Astronomy_SearchLocalSolarEclipse(SkyEngine.time(peak.addingTimeInterval(-6 * 3600)), at.astro)
 50        var duration: TimeInterval?
 51        if local.status == ASTRO_SUCCESS, local.kind == g.kind {
 52            duration = SkyEngine.date(local.total_end.time).timeIntervalSince(SkyEngine.date(local.total_begin.time))
 53        }
 54        return EclipsePath(kind: SkyEngine.eclipseKind(g.kind), peak: peak, greatest: (g.latitude, g.longitude), centralDuration: duration, points: points)
 55    }
 56
 57    /// Centre and limits of the path at `date`, or nil if the shadow axis misses the Earth.
 58    static func point(at date: Date, step: TimeInterval) -> Point? {
 59        guard let (center, n, r) = axis(at: date), let next = axis(at: date.addingTimeInterval(step)) else { return nil }
 60        guard let c = center.hit, let c2 = next.0.hit else { return nil }
 61        // Unit vector across the track: perpendicular to the shadow axis and the ground motion.
 62        let across = normalize(cross(n, c2 - c))
 63        let offset = across * (abs(r) / kmPerAU)
 64        let north = Ray(origin: center.origin + offset, direction: n).hit
 65        let south = Ray(origin: center.origin - offset, direction: n).hit
 66        let t = SkyEngine.time(date)
 67        func geo(_ v: SIMD3<Double>) -> (lat: Double, lon: Double) {
 68            var vec = astro_vector_t(status: ASTRO_SUCCESS, x: v.x, y: v.y, z: v.z, t: t)
 69            let o = Astronomy_VectorObserver(&vec, EQUATOR_OF_DATE)
 70            return (o.latitude, o.longitude)
 71        }
 72        var a = geo(north ?? c), b = geo(south ?? c)
 73        if a.lat < b.lat { swap(&a, &b) }
 74        let width = (north != nil && south != nil) ? length(north! - south!) * kmPerAU : 2 * abs(r)
 75        return Point(date: date, center: geo(c), north: north == nil ? nil : a, south: south == nil ? nil : b, widthKm: width)
 76    }
 77
 78    struct Ray {
 79        var origin: SIMD3<Double>
 80        var direction: SIMD3<Double>
 81
 82        /// First intersection with the Earth's ellipsoid, AU, equator-of-date axes.
 83        var hit: SIMD3<Double>? {
 84            let k = 1 / (1 - EclipsePath.earthF)
 85            let o = SIMD3(origin.x, origin.y, origin.z * k), d = SIMD3(direction.x, direction.y, direction.z * k)
 86            let a = EclipsePath.earthA / EclipsePath.kmPerAU
 87            let A = dot(d, d), B = 2 * dot(o, d), C = dot(o, o) - a * a
 88            let disc = B * B - 4 * A * C
 89            guard disc >= 0 else { return nil }
 90            let s = (-B - disc.squareRoot()) / (2 * A)
 91            guard s > 0 else { return nil }
 92            let p = o + d * s
 93            return SIMD3(p.x, p.y, p.z / k)
 94        }
 95    }
 96
 97    /// Shadow axis from the Moon away from the Sun (equator of date), its direction,
 98    /// and the umbra radius where it meets the Earth (km; negative for the antumbra).
 99    static func axis(at date: Date) -> (Ray, SIMD3<Double>, Double)? {
100        var t = SkyEngine.time(date)
101        let rot = Astronomy_Rotation_EQJ_EQD(&t)
102        let sun = Astronomy_RotateVector(rot, Astronomy_GeoVector(BODY_SUN, t, NO_ABERRATION))
103        let moon = Astronomy_RotateVector(rot, Astronomy_GeoMoon(t))
104        let s = SIMD3(sun.x, sun.y, sun.z), m = SIMD3(moon.x, moon.y, moon.z)
105        let n = normalize(m - s)
106        let ray = Ray(origin: m, direction: n)
107        guard let hit = ray.hit else { return nil }
108        let D = length(m - s) * kmPerAU
109        let L = length(hit - m) * kmPerAU
110        let r = moonRadiusKm - L * (sunRadiusKm - moonRadiusKm) / D
111        return (ray, n, r)
112    }
113}
114
115private func dot(_ a: SIMD3<Double>, _ b: SIMD3<Double>) -> Double { (a * b).sum() }
116private func length(_ a: SIMD3<Double>) -> Double { dot(a, a).squareRoot() }
117private func normalize(_ a: SIMD3<Double>) -> SIMD3<Double> { a / length(a) }
118private func cross(_ a: SIMD3<Double>, _ b: SIMD3<Double>) -> SIMD3<Double> {
119    SIMD3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x)
120}
121
122/// A solar eclipse as seen from one place.
123public struct LocalEclipse: Sendable {
124    public var kind: String
125    public var obscuration: Double
126    public var partialBegin: Date
127    public var centralBegin: Date?
128    public var peak: Date
129    public var centralEnd: Date?
130    public var partialEnd: Date
131    /// Sun's altitude at the peak, degrees.
132    public var peakAltitude: Double
133
134    /// The observer's view of the solar eclipse peaking within six hours of `date`, if any part is above the horizon.
135    public static func at(_ observer: Observer, near date: Date) -> LocalEclipse? {
136        SkyEngine.lock.lock()
137        defer { SkyEngine.lock.unlock() }
138        let l = Astronomy_SearchLocalSolarEclipse(SkyEngine.time(date.addingTimeInterval(-6 * 3600)), observer.astro)
139        guard l.status == ASTRO_SUCCESS else { return nil }
140        let peak = SkyEngine.date(l.peak.time)
141        guard abs(peak.timeIntervalSince(date)) < 6 * 3600 else { return nil }
142        guard l.partial_begin.altitude > 0 || l.peak.altitude > 0 || l.partial_end.altitude > 0 else { return nil }
143        let central = l.kind == ECLIPSE_TOTAL || l.kind == ECLIPSE_ANNULAR
144        return LocalEclipse(
145            kind: SkyEngine.eclipseKind(l.kind), obscuration: l.obscuration,
146            partialBegin: SkyEngine.date(l.partial_begin.time),
147            centralBegin: central ? SkyEngine.date(l.total_begin.time) : nil,
148            peak: peak,
149            centralEnd: central ? SkyEngine.date(l.total_end.time) : nil,
150            partialEnd: SkyEngine.date(l.partial_end.time),
151            peakAltitude: l.peak.altitude
152        )
153    }
154}