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
 36extension OrbitalElements {
 37    public static let featuredIDs: Set<Int> = [25544, 48274, 20580]
 38
 39    public var isFeatured: Bool { Self.featuredIDs.contains(NORAD_CAT_ID) || OBJECT_NAME.hasPrefix("Starlink train") }
 40}
 41
 42extension SkyEngine {
 43    /// Eight-point compass direction for an azimuth in degrees.
 44    public static func compass(_ azimuth: Double) -> String {
 45        ["N", "NE", "E", "SE", "S", "SW", "W", "NW"][Int(((azimuth + 22.5) / 45).rounded(.down)) % 8]
 46    }
 47
 48    /// Greenwich mean sidereal time in radians (IAU 1982, as used with SGP4's TEME frame).
 49    static func gmst(_ date: Date) -> Double {
 50        let jd = date.timeIntervalSince1970 / 86_400 + 2_440_587.5
 51        let t = (jd - 2_451_545) / 36_525
 52        var s = -6.2e-6 * t * t * t + 0.093104 * t * t + (876_600 * 3600 + 8_640_184.812866) * t + 67_310.54841
 53        s = fmod(s * Double.pi / 180 / 240, 2 * Double.pi)
 54        return s < 0 ? s + 2 * Double.pi : s
 55    }
 56
 57    /// Observer position in an Earth-fixed frame (km), WGS-84.
 58    static func observerECEF(_ o: Observer) -> SIMD3<Double> {
 59        let a = 6378.137, f = 1 / 298.257223563
 60        let e2 = f * (2 - f)
 61        let lat = o.latitude * .pi / 180, lon = o.longitude * .pi / 180, h = o.height / 1000
 62        let n = a / (1 - e2 * sin(lat) * sin(lat)).squareRoot()
 63        return SIMD3((n + h) * cos(lat) * cos(lon), (n + h) * cos(lat) * sin(lon), (n * (1 - e2) + h) * sin(lat))
 64    }
 65
 66    /// Altitude and azimuth of a TEME position as seen by the observer.
 67    static func lookAngles(teme r: SIMD3<Double>, at date: Date, from o: Observer) -> Horizontal {
 68        let g = gmst(date)
 69        let ecef = SIMD3(r.x * cos(g) + r.y * sin(g), -r.x * sin(g) + r.y * cos(g), r.z)
 70        let d = ecef - observerECEF(o)
 71        let lat = o.latitude * .pi / 180, lon = o.longitude * .pi / 180
 72        let south = sin(lat) * cos(lon) * d.x + sin(lat) * sin(lon) * d.y - cos(lat) * d.z
 73        let east = -sin(lon) * d.x + cos(lon) * d.y
 74        let up = cos(lat) * cos(lon) * d.x + cos(lat) * sin(lon) * d.y + sin(lat) * d.z
 75        let range = (d * d).sum().squareRoot()
 76        var az = atan2(east, -south) * 180 / .pi
 77        if az < 0 { az += 360 }
 78        return Horizontal(altitude: asin(up / range) * 180 / .pi, azimuth: az)
 79    }
 80
 81    /// Unit vector to the Sun on equator-of-date axes, close enough to TEME for Earth's shadow.
 82    static func sunDirection(_ date: Date) -> SIMD3<Double> {
 83        var t = time(date)
 84        let v = Astronomy_GeoVector(BODY_SUN, t, NO_ABERRATION)
 85        let rot = Astronomy_Rotation_EQJ_EQD(&t)
 86        let e = Astronomy_RotateVector(rot, v)
 87        let s = SIMD3(e.x, e.y, e.z)
 88        return s / (s * s).sum().squareRoot()
 89    }
 90
 91    /// Cylindrical Earth-shadow test.
 92    static func sunlit(_ r: SIMD3<Double>, sun: SIMD3<Double>) -> Bool {
 93        let along = (r * sun).sum()
 94        if along > 0 { return true }
 95        let perp = r - along * sun
 96        return (perp * perp).sum().squareRoot() > 6378.135
 97    }
 98
 99    /// Passes above `minAltitude` in the next `hours`, sampled every 20 s.
100    public static func passes(of elements: OrbitalElements, from start: Date, hours: Double = 72, observer: Observer, minAltitude: Double = 10) -> [SatellitePass] {
101        guard let sgp = try? SGP4(elements) else { return [] }
102        lock.lock()
103        defer { lock.unlock() }
104        var out: [SatellitePass] = []
105        var current: (rise: Date, riseAz: Double, peak: Date, maxAlt: Double, visible: Bool)?
106        var lastAz = 0.0
107        let step = 20.0
108        var t = start
109        let end = start.addingTimeInterval(hours * 3600)
110        while t <= end {
111            guard let r = sgp.position(at: t) else { break }
112            let h = lookAngles(teme: r, at: t, from: observer)
113            if h.altitude >= 0 {
114                if current == nil { current = (t, h.azimuth, t, h.altitude, false) }
115                if h.altitude > current!.maxAlt { current!.maxAlt = h.altitude; current!.peak = t }
116                if !current!.visible, h.altitude >= minAltitude {
117                    var tt = time(t)
118                    let sun = Astronomy_Equator(BODY_SUN, &tt, observer.astro, EQUATOR_OF_DATE, ABERRATION)
119                    let sunAlt = Astronomy_Horizon(&tt, observer.astro, sun.ra, sun.dec, REFRACTION_NORMAL).altitude
120                    if sunAlt < -6, sunlit(r, sun: sunDirection(t)) { current!.visible = true }
121                }
122                lastAz = h.azimuth
123            } else if let c = current {
124                if c.maxAlt >= minAltitude {
125                    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))
126                }
127                current = nil
128            }
129            t = t.addingTimeInterval(step)
130        }
131        return out
132    }
133
134    /// Visible passes of any of `elements`, earliest first, each set searched for `hours` from the later of `start`
135    /// and its epoch: a set published before launch is epoched at deployment and means nothing earlier.
136    public static func visiblePasses(of elements: [OrbitalElements], from start: Date, hours: Double = 72, observer: Observer) -> [SatellitePass] {
137        elements.flatMap { e in
138            passes(of: e, from: max(start, e.epochDate ?? start), hours: hours, observer: observer).filter(\.visible)
139        }
140        .sorted { $0.rise < $1.rise }
141    }
142}