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