Packages/MoonriseKit/Sources/SkyEngine/Passes.swift

  1import CAstronomy
  2import Foundation
  3
  4public struct SatellitePass: Sendable, Identifiable, Equatable {
  5    public var name: String
  6    public var rise: Date
  7    public var set: Date
  8    public var peak: Date
  9    public var maxAltitude: Double
 10    public var riseAzimuth: Double
 11    public var setAzimuth: Double
 12    /// Sunlit while the observer's sky is dark: can be seen with the eye.
 13    public var visible: Bool
 14    /// ISS, Tiangong, Hubble or a Starlink train, as opposed to the rest of CelesTrak's visual group.
 15    public var featured = false
 16    public var id: String { "\(name)-\(Int(rise.timeIntervalSince1970))" }
 17
 18    /// "ISS" for "ISS (ZARYA)", "Tiangong" for "CSS (TIANHE)", "SL-8 rocket body" for "SL-8 R/B".
 19    public var displayName: String {
 20        switch name {
 21        case "ISS (ZARYA)": "ISS"
 22        case "CSS (TIANHE)": "Tiangong"
 23        case "HST": "Hubble"
 24        default: name.replacingOccurrences(of: " R/B", with: " rocket body")
 25        }
 26    }
 27
 28    /// When a Lock Screen activity for this pass starts: 10 minutes before rise, or `now` if that has passed.
 29    /// Nil once the pass has set, or when rise is more than a day away and the elements will have drifted.
 30    public func activityStart(now: Date) -> Date? {
 31        guard set > now, rise.timeIntervalSince(now) <= 86_400 else { return nil }
 32        return max(now, rise.addingTimeInterval(-600))
 33    }
 34
 35    /// "Moon up, 80% lit" for a visible pass with the Moon above the horizon at its peak and lit past `moonlightLimit`, when
 36    /// moonlight can hide a faint satellite. Nil for the ISS and Tiangong, which are bright enough that it does not matter.
 37    public func moonNote(from observer: Observer) -> String? {
 38        guard visible, name != "ISS (ZARYA)", name != "CSS (TIANHE)",
 39              let moon = try? SkyEngine.horizontal(.moon, at: peak, from: observer, refraction: true), moon.altitude > 0,
 40              let lit = try? SkyEngine.illumination(.moon, at: peak).phaseFraction, lit >= SkyEngine.moonlightLimit else { return nil }
 41        return "Moon up, \(Int((lit * 100).rounded()))% lit"
 42    }
 43}
 44
 45extension OrbitalElements {
 46    public static let featuredIDs: Set<Int> = [25544, 48274, 20580]
 47
 48    public var isFeatured: Bool { Self.featuredIDs.contains(NORAD_CAT_ID) || OBJECT_NAME.hasPrefix("Starlink train") }
 49}
 50
 51extension SkyEngine {
 52    /// Eight-point compass direction for an azimuth in degrees.
 53    public static func compass(_ azimuth: Double) -> String {
 54        ["N", "NE", "E", "SE", "S", "SW", "W", "NW"][Int(((azimuth + 22.5) / 45).rounded(.down)) % 8]
 55    }
 56
 57    /// Greenwich mean sidereal time in radians (IAU 1982, as used with SGP4's TEME frame).
 58    static func gmst(_ date: Date) -> Double {
 59        let jd = date.timeIntervalSince1970 / 86_400 + 2_440_587.5
 60        let t = (jd - 2_451_545) / 36_525
 61        var s = -6.2e-6 * t * t * t + 0.093104 * t * t + (876_600 * 3600 + 8_640_184.812866) * t + 67_310.54841
 62        s = fmod(s * Double.pi / 180 / 240, 2 * Double.pi)
 63        return s < 0 ? s + 2 * Double.pi : s
 64    }
 65
 66    /// Observer position in an Earth-fixed frame (km), WGS-84.
 67    static func observerECEF(_ o: Observer) -> SIMD3<Double> {
 68        let a = 6378.137, f = 1 / 298.257223563
 69        let e2 = f * (2 - f)
 70        let lat = o.latitude * .pi / 180, lon = o.longitude * .pi / 180, h = o.height / 1000
 71        let n = a / (1 - e2 * sin(lat) * sin(lat)).squareRoot()
 72        return SIMD3((n + h) * cos(lat) * cos(lon), (n + h) * cos(lat) * sin(lon), (n * (1 - e2) + h) * sin(lat))
 73    }
 74
 75    /// Altitude and azimuth of a TEME position as seen by the observer.
 76    static func lookAngles(teme r: SIMD3<Double>, at date: Date, from o: Observer) -> Horizontal {
 77        let g = gmst(date)
 78        let ecef = SIMD3(r.x * cos(g) + r.y * sin(g), -r.x * sin(g) + r.y * cos(g), r.z)
 79        let d = ecef - observerECEF(o)
 80        let lat = o.latitude * .pi / 180, lon = o.longitude * .pi / 180
 81        let south = sin(lat) * cos(lon) * d.x + sin(lat) * sin(lon) * d.y - cos(lat) * d.z
 82        let east = -sin(lon) * d.x + cos(lon) * d.y
 83        let up = cos(lat) * cos(lon) * d.x + cos(lat) * sin(lon) * d.y + sin(lat) * d.z
 84        let range = (d * d).sum().squareRoot()
 85        var az = atan2(east, -south) * 180 / .pi
 86        if az < 0 { az += 360 }
 87        return Horizontal(altitude: asin(up / range) * 180 / .pi, azimuth: az)
 88    }
 89
 90    /// Unit vector to the Sun on equator-of-date axes, close enough to TEME for Earth's shadow.
 91    static func sunDirection(_ date: Date) -> SIMD3<Double> {
 92        var t = time(date)
 93        let v = Astronomy_GeoVector(BODY_SUN, t, NO_ABERRATION)
 94        let rot = Astronomy_Rotation_EQJ_EQD(&t)
 95        let e = Astronomy_RotateVector(rot, v)
 96        let s = SIMD3(e.x, e.y, e.z)
 97        return s / (s * s).sum().squareRoot()
 98    }
 99
100    /// Cylindrical Earth-shadow test.
101    static func sunlit(_ r: SIMD3<Double>, sun: SIMD3<Double>) -> Bool {
102        let along = (r * sun).sum()
103        if along > 0 { return true }
104        let perp = r - along * sun
105        return (perp * perp).sum().squareRoot() > 6378.135
106    }
107
108    /// Passes above `minAltitude` in the next `hours`, sampled every 20 s.
109    public static func passes(of elements: OrbitalElements, from start: Date, hours: Double = 72, observer: Observer, minAltitude: Double = 10) -> [SatellitePass] {
110        guard let sgp = try? SGP4(elements) else { return [] }
111        lock.lock()
112        defer { lock.unlock() }
113        var out: [SatellitePass] = []
114        var current: (rise: Date, riseAz: Double, peak: Date, maxAlt: Double, visible: Bool)?
115        var lastAz = 0.0
116        let step = 20.0
117        var t = start
118        let end = start.addingTimeInterval(hours * 3600)
119        while t <= end {
120            guard let r = sgp.position(at: t) else { break }
121            let h = lookAngles(teme: r, at: t, from: observer)
122            if h.altitude >= 0 {
123                if current == nil { current = (t, h.azimuth, t, h.altitude, false) }
124                if h.altitude > current!.maxAlt { current!.maxAlt = h.altitude; current!.peak = t }
125                if !current!.visible, h.altitude >= minAltitude {
126                    var tt = time(t)
127                    let sun = Astronomy_Equator(BODY_SUN, &tt, observer.astro, EQUATOR_OF_DATE, ABERRATION)
128                    let sunAlt = Astronomy_Horizon(&tt, observer.astro, sun.ra, sun.dec, REFRACTION_NORMAL).altitude
129                    if sunAlt < -6, sunlit(r, sun: sunDirection(t)) { current!.visible = true }
130                }
131                lastAz = h.azimuth
132            } else if let c = current {
133                if c.maxAlt >= minAltitude {
134                    out.append(SatellitePass(name: elements.OBJECT_NAME, rise: c.rise, set: t, peak: c.peak, maxAltitude: c.maxAlt, riseAzimuth: c.riseAz, setAzimuth: lastAz, visible: c.visible, featured: elements.isFeatured))
135                }
136                current = nil
137            }
138            t = t.addingTimeInterval(step)
139        }
140        return out
141    }
142
143    /// Visible passes of any of `elements`, earliest first, each set searched for `hours` from the later of `start`
144    /// and its epoch: a set published before launch is epoched at deployment and means nothing earlier.
145    public static func visiblePasses(of elements: [OrbitalElements], from start: Date, hours: Double = 72, observer: Observer) -> [SatellitePass] {
146        elements.flatMap { e in
147            passes(of: e, from: max(start, e.epochDate ?? start), hours: hours, observer: observer).filter(\.visible)
148        }
149        .sorted { $0.rise < $1.rise }
150    }
151}