Packages/MoonriseKit/Sources/SkyEngine/Kepler.swift

  1import CAstronomy
  2import Foundation
  3
  4/// A two-body orbit around the Sun from osculating elements (J2000 ecliptic), for small bodies the
  5/// engine doesn't cover. Planetary perturbations are ignored, so positions drift away from the epoch:
  6/// fine for drawing, not for observing.
  7public struct KeplerOrbit: Sendable, Identifiable {
  8    public var id: String
  9    public var name: String
 10    /// Library page, if there is one.
 11    public var page: String?
 12    public var a: Double, e: Double, i: Double, node: Double, peri: Double
 13    /// Time of perihelion (TDB Julian date). Propagating from it rather than from the epoch's mean
 14    /// anomaly keeps comets right near perihelion, where they are worth drawing.
 15    public var perihelionJD: Double
 16    public var isComet: Bool
 17
 18    /// Orbital period in days.
 19    public var periodDays: Double { 365.256898 * pow(a, 1.5) }
 20
 21    /// Heliocentric position in AU on J2000 ecliptic axes.
 22    public func position(at date: Date) -> SIMD3<Double> {
 23        let jd = date.timeIntervalSince1970 / 86_400 + 2_440_587.5
 24        let n = 360 / periodDays
 25        let M = ((n * (jd - perihelionJD)).truncatingRemainder(dividingBy: 360) + 360).truncatingRemainder(dividingBy: 360) * .pi / 180
 26        var E = e < 0.8 ? M : .pi
 27        for _ in 0..<50 {
 28            let d = (E - e * sin(E) - M) / (1 - e * cos(E))
 29            E -= d
 30            if abs(d) < 1e-12 { break }
 31        }
 32        let x = a * (cos(E) - e), y = a * (1 - e * e).squareRoot() * sin(E)
 33        let r = Double.pi / 180
 34        let (cO, sO) = (cos(node * r), sin(node * r)), (cw, sw) = (cos(peri * r), sin(peri * r)), (ci, si) = (cos(i * r), sin(i * r))
 35        return SIMD3(
 36            (cO * cw - sO * sw * ci) * x + (-cO * sw - sO * cw * ci) * y,
 37            (sO * cw + cO * sw * ci) * x + (-sO * sw + cO * cw * ci) * y,
 38            (sw * si) * x + (cw * si) * y
 39        )
 40    }
 41
 42    /// One full orbit, sampled evenly in eccentric anomaly so long, thin comet orbits stay smooth near the Sun.
 43    public func orbitPath(steps: Int = 240) -> [SIMD3<Double>] {
 44        let r = Double.pi / 180
 45        let (cO, sO) = (cos(node * r), sin(node * r)), (cw, sw) = (cos(peri * r), sin(peri * r)), (ci, si) = (cos(i * r), sin(i * r))
 46        return (0..<steps).map { k in
 47            let E = Double(k) / Double(steps) * 2 * .pi
 48            let x = a * (cos(E) - e), y = a * (1 - e * e).squareRoot() * sin(E)
 49            return SIMD3(
 50                (cO * cw - sO * sw * ci) * x + (-cO * sw - sO * cw * ci) * y,
 51                (sO * cw + cO * sw * ci) * x + (-sO * sw + cO * cw * ci) * y,
 52                (sw * si) * x + (cw * si) * y
 53            )
 54        }
 55    }
 56
 57    /// The next perihelion after `date`.
 58    public func nextPerihelion(after date: Date) -> Date {
 59        let jd = date.timeIntervalSince1970 / 86_400 + 2_440_587.5
 60        let k = ((jd - perihelionJD) / periodDays).rounded(.up)
 61        return Date(timeIntervalSince1970: (perihelionJD + k * periodDays - 2_440_587.5) * 86_400)
 62    }
 63
 64    /// JPL Small-Body Database osculating elements (full precision), fetched 2026-10-19.
 65    public static let dwarfPlanets: [KeplerOrbit] = [
 66        KeplerOrbit(id: "ceres", name: "Ceres", page: "ceres", a: 2.765552595034094, e: 0.07969229514816586, i: 10.58802780183462, node: 80.24862682043221, peri: 73.29421453021587, perihelionJD: 2461599.841466614, isComet: false),
 67        KeplerOrbit(id: "haumea", name: "Haumea", page: "haumea", a: 43.06029023650952, e: 0.1944430148898797, i: 28.20847393040364, node: 121.7860561329425, peri: 240.6905472508661, perihelionJD: 2500416.599615283, isComet: false),
 68        KeplerOrbit(id: "makemake", name: "Makemake", page: "makemake", a: 45.57093317300052, e: 0.1588889953992523, i: 29.02785603743067, node: 79.2948338209406, peri: 297.0922733397207, perihelionJD: 2408158.6940986165, isComet: false),
 69        KeplerOrbit(id: "eris", name: "Eris", page: "eris", a: 67.93394687853566, e: 0.4382385347971672, i: 43.9258279471791, node: 36.00477044417249, peri: 150.7949235840312, perihelionJD: 2545407.716847013, isComet: false),
 70    ]
 71
 72    public static let comets: [KeplerOrbit] = [
 73        KeplerOrbit(id: "halley", name: "Halley's Comet", page: nil, a: 17.92863504856923, e: 0.9679359956953211, i: 162.1905300439129, node: 59.09894720612437, peri: 112.2414314637764, perihelionJD: 2446469.9736161465, isComet: true),
 74        KeplerOrbit(id: "encke", name: "Comet Encke", page: nil, a: 2.218624763741819, e: 0.846940430888432, i: 11.33675914978793, node: 334.0171400057891, peri: 187.2881998732444, perihelionJD: 2460240.0285326503, isComet: true),
 75        KeplerOrbit(id: "67p", name: "67P/Churyumov–Gerasimenko", page: "rosetta", a: 3.462249490129549, e: 0.6409081308996354, i: 7.040294937543767, node: 50.13557377155012, peri: 12.79824970228189, perihelionJD: 2457247.588657812, isComet: true),
 76        KeplerOrbit(id: "swift-tuttle", name: "Comet Swift–Tuttle", page: nil, a: 26.0920694978266, e: 0.963225755046038, i: 113.453816997171, node: 139.3811920815948, peri: 152.9821676305871, perihelionJD: 2448968.499784556, isComet: true),
 77        KeplerOrbit(id: "tempel-tuttle", name: "Comet Tempel–Tuttle", page: nil, a: 10.3383382297577, e: 0.905552720972412, i: 162.486575379434, node: 235.270989149082, peri: 172.5002736828059, perihelionJD: 2450872.5977341123, isComet: true),
 78        KeplerOrbit(id: "phaethon", name: "3200 Phaethon", page: nil, a: 1.271464620762422, e: 0.8896722844001859, i: 22.31052728977475, node: 265.0988060326684, peri: 322.3001684954223, perihelionJD: 2461285.6164380596, isComet: true),
 79        KeplerOrbit(id: "tuttle", name: "Comet 8P/Tuttle", page: nil, a: 5.700291361749948, e: 0.8197876162331501, i: 54.98196882510133, node: 270.3421615706648, peri: 207.5115175085094, perihelionJD: 2454492.5180154694, isComet: true),
 80    ]
 81}
 82
 83extension MeteorShower {
 84    /// The comet (or asteroid) whose debris Earth passes through, by KeplerOrbit id.
 85    public var parent: String? {
 86        switch name {
 87        case "Perseids": "swift-tuttle"
 88        case "Leonids": "tempel-tuttle"
 89        case "Orionids", "Eta Aquariids": "halley"
 90        case "Geminids": "phaethon"
 91        case "Ursids": "tuttle"
 92        default: nil
 93        }
 94    }
 95
 96    /// Days from `date` to the nearest peak (negative once it has passed), in the given calendar.
 97    public func daysToPeak(from date: Date, calendar: Calendar = .current) -> Double {
 98        let year = calendar.component(.year, from: date)
 99        return (year - 1...year + 1).compactMap { y in
100            calendar.date(from: DateComponents(year: y, month: month, day: day, hour: 12)).map { $0.timeIntervalSince(date) / 86_400 }
101        }
102        .min { abs($0) < abs($1) } ?? .infinity
103    }
104}
105
106extension SkyEngine {
107    /// Io, Europa, Ganymede and Callisto relative to Jupiter, AU on J2000 ecliptic axes.
108    public static func galileanMoons(at date: Date) -> [SIMD3<Double>] {
109        lock.lock()
110        defer { lock.unlock() }
111        let m = Astronomy_JupiterMoons(time(date))
112        return [m.io, m.europa, m.ganymede, m.callisto].map { ecliptic(SIMD3($0.x, $0.y, $0.z)) }
113    }
114}