v2: Satellite passes !24
9 files changed, +517 −1
Layout: unified · split
App/Sources/SkyModel.swift +20
| @@ -51,6 +51,9 @@ final class SkyModel: NSObject, CLLocationManagerDelegate { | |||
| 51 | private(set) var lastBuild: Duration = .zero | 51 | private(set) var lastBuild: Duration = .zero |
| 52 | private(set) var tonight: Tonight? | 52 | private(set) var tonight: Tonight? |
| 53 | private(set) var spaceWeather: SpaceWeather? | 53 | private(set) var spaceWeather: SpaceWeather? |
| 54 | private(set) var passes: [SatellitePass] = [] | ||
| 55 | private var elements: [OrbitalElements] = [] | ||
| 56 | private var elementsFetched: Date? | ||
| 54 | private var weatherFetched: Date? | 57 | private var weatherFetched: Date? |
| 55 | /// Hours from now shown by the time scrubber. | 58 | /// Hours from now shown by the time scrubber. |
| 56 | var offsetHours = 0.0 { didSet { refresh() } } | 59 | var offsetHours = 0.0 { didSet { refresh() } } |
| @@ -109,6 +112,23 @@ final class SkyModel: NSObject, CLLocationManagerDelegate { | |||
| 109 | } | 112 | } |
| 110 | } | 113 | } |
| 111 | 114 | ||
| 115 | private struct SatelliteList: Decodable { var satellites: [OrbitalElements] } | ||
| 116 | |||
| 117 | /// Orbital elements from the backend (refreshed hourly), then passes for the next three days. | ||
| 118 | func refreshPasses() async { | ||
| 119 | if elementsFetched.map({ Date.now.timeIntervalSince($0) > 3600 }) ?? true, let apiBase { | ||
| 120 | if let (data, _) = try? await URLSession.shared.data(from: apiBase.appending(path: "v1/satellites")), | ||
| 121 | let list = try? JSONDecoder().decode(SatelliteList.self, from: data), !list.satellites.isEmpty { | ||
| 122 | elements = list.satellites | ||
| 123 | elementsFetched = .now | ||
| 124 | } | ||
| 125 | } | ||
| 126 | let elements = elements, observer = observer | ||
| 127 | passes = await Task.detached { | ||
| 128 | elements.flatMap { SkyEngine.passes(of: $0, from: .now, hours: 72, observer: observer) }.sorted { $0.rise < $1.rise } | ||
| 129 | }.value | ||
| 130 | } | ||
| 131 | |||
| 112 | /// SWPC asks for no more than once a minute; this checks every 10. | 132 | /// SWPC asks for no more than once a minute; this checks every 10. |
| 113 | func refreshSpaceWeather() async { | 133 | func refreshSpaceWeather() async { |
| 114 | if let weatherFetched, Date.now.timeIntervalSince(weatherFetched) < 600 { return } | 134 | if let weatherFetched, Date.now.timeIntervalSince(weatherFetched) < 600 { return } |
App/Sources/TodaySheet.swift +28
| @@ -16,6 +16,7 @@ struct TodaySheet: View { | |||
| 16 | section("Planets tonight") { planets(t) } | 16 | section("Planets tonight") { planets(t) } |
| 17 | section("Sun") { sun(t) } | 17 | section("Sun") { sun(t) } |
| 18 | } | 18 | } |
| 19 | section("Satellite passes") { satellites } | ||
| 19 | section("Space weather") { weather } | 20 | section("Space weather") { weather } |
| 20 | Text("Space weather data: NOAA Space Weather Prediction Center.") | 21 | Text("Space weather data: NOAA Space Weather Prediction Center.") |
| 21 | .font(.caption2).dimmed(0.5).padding(.top, 8) | 22 | .font(.caption2).dimmed(0.5).padding(.top, 8) |
| @@ -25,6 +26,7 @@ struct TodaySheet: View { | |||
| 25 | } | 26 | } |
| 26 | .foregroundStyle(ink) | 27 | .foregroundStyle(ink) |
| 27 | .presentationBackground(background) | 28 | .presentationBackground(background) |
| 29 | .task { await model.refreshPasses() } | ||
| 28 | } | 30 | } |
| 29 | 31 | ||
| 30 | private func section(_ title: String, @ViewBuilder content: () -> some View) -> some View { | 32 | private func section(_ title: String, @ViewBuilder content: () -> some View) -> some View { |
| @@ -77,6 +79,32 @@ struct TodaySheet: View { | |||
| 77 | ]) | 79 | ]) |
| 78 | } | 80 | } |
| 79 | 81 | ||
| 82 | @ViewBuilder | ||
| 83 | private var satellites: some View { | ||
| 84 | let shown = model.passes.prefix(8) | ||
| 85 | if shown.isEmpty { | ||
| 86 | Text("No passes above 10° in the next three days.").font(.subheadline).dimmed(0.6) | ||
| 87 | } else { | ||
| 88 | VStack(alignment: .leading, spacing: 10) { | ||
| 89 | ForEach(Array(shown)) { p in | ||
| 90 | VStack(alignment: .leading, spacing: 2) { | ||
| 91 | HStack { | ||
| 92 | Text(p.name.replacingOccurrences(of: " (ZARYA)", with: "")).font(.subheadline.weight(.medium)) | ||
| 93 | if p.visible { | ||
| 94 | Text("Visible").font(.caption2.weight(.medium)).padding(.horizontal, 5).padding(.vertical, 1) | ||
| 95 | .overlay(RoundedRectangle(cornerRadius: 3).stroke(lineWidth: 1)) | ||
| 96 | } | ||
| 97 | } | ||
| 98 | Text(p.rise.formatted(.dateTime.weekday(.abbreviated).hour().minute()) + String(format: " · rises %@, highest %.0f°, sets %@", compass(p.riseAzimuth), p.maxAltitude, compass(p.setAzimuth))) | ||
| 99 | .font(.footnote.monospacedDigit()).dimmed(0.8) | ||
| 100 | } | ||
| 101 | .accessibilityElement(children: .combine) | ||
| 102 | } | ||
| 103 | Text("Visible means sunlit while your sky is dark. Orbits: CelesTrak.").font(.caption2).dimmed(0.6) | ||
| 104 | } | ||
| 105 | } | ||
| 106 | } | ||
| 107 | |||
| 80 | @ViewBuilder | 108 | @ViewBuilder |
| 81 | private var weather: some View { | 109 | private var weather: some View { |
| 82 | if let w = model.spaceWeather { | 110 | if let w = model.spaceWeather { |
Packages/MoonriseKit/Sources/SkyEngine/Passes.swift added +103
| @@ -0,0 +1,103 @@ | |||
| 1 | import CAstronomy | ||
| 2 | import Foundation | ||
| 3 | |||
| 4 | public 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 | public var id: String { "\(name)-\(Int(rise.timeIntervalSince1970))" } | ||
| 15 | } | ||
| 16 | |||
| 17 | extension SkyEngine { | ||
| 18 | /// Greenwich mean sidereal time in radians (IAU 1982, as used with SGP4's TEME frame). | ||
| 19 | static func gmst(_ date: Date) -> Double { | ||
| 20 | let jd = date.timeIntervalSince1970 / 86_400 + 2_440_587.5 | ||
| 21 | let t = (jd - 2_451_545) / 36_525 | ||
| 22 | var s = -6.2e-6 * t * t * t + 0.093104 * t * t + (876_600 * 3600 + 8_640_184.812866) * t + 67_310.54841 | ||
| 23 | s = fmod(s * Double.pi / 180 / 240, 2 * Double.pi) | ||
| 24 | return s < 0 ? s + 2 * Double.pi : s | ||
| 25 | } | ||
| 26 | |||
| 27 | /// Observer position in an Earth-fixed frame (km), WGS-84. | ||
| 28 | static func observerECEF(_ o: Observer) -> SIMD3<Double> { | ||
| 29 | let a = 6378.137, f = 1 / 298.257223563 | ||
| 30 | let e2 = f * (2 - f) | ||
| 31 | let lat = o.latitude * .pi / 180, lon = o.longitude * .pi / 180, h = o.height / 1000 | ||
| 32 | let n = a / (1 - e2 * sin(lat) * sin(lat)).squareRoot() | ||
| 33 | return SIMD3((n + h) * cos(lat) * cos(lon), (n + h) * cos(lat) * sin(lon), (n * (1 - e2) + h) * sin(lat)) | ||
| 34 | } | ||
| 35 | |||
| 36 | /// Altitude and azimuth of a TEME position as seen by the observer. | ||
| 37 | static func lookAngles(teme r: SIMD3<Double>, at date: Date, from o: Observer) -> Horizontal { | ||
| 38 | let g = gmst(date) | ||
| 39 | let ecef = SIMD3(r.x * cos(g) + r.y * sin(g), -r.x * sin(g) + r.y * cos(g), r.z) | ||
| 40 | let d = ecef - observerECEF(o) | ||
| 41 | let lat = o.latitude * .pi / 180, lon = o.longitude * .pi / 180 | ||
| 42 | let south = sin(lat) * cos(lon) * d.x + sin(lat) * sin(lon) * d.y - cos(lat) * d.z | ||
| 43 | let east = -sin(lon) * d.x + cos(lon) * d.y | ||
| 44 | let up = cos(lat) * cos(lon) * d.x + cos(lat) * sin(lon) * d.y + sin(lat) * d.z | ||
| 45 | let range = (d * d).sum().squareRoot() | ||
| 46 | var az = atan2(east, -south) * 180 / .pi | ||
| 47 | if az < 0 { az += 360 } | ||
| 48 | return Horizontal(altitude: asin(up / range) * 180 / .pi, azimuth: az) | ||
| 49 | } | ||
| 50 | |||
| 51 | /// Unit vector to the Sun on equator-of-date axes, close enough to TEME for Earth's shadow. | ||
| 52 | static func sunDirection(_ date: Date) -> SIMD3<Double> { | ||
| 53 | var t = time(date) | ||
| 54 | let v = Astronomy_GeoVector(BODY_SUN, t, NO_ABERRATION) | ||
| 55 | let rot = Astronomy_Rotation_EQJ_EQD(&t) | ||
| 56 | let e = Astronomy_RotateVector(rot, v) | ||
| 57 | let s = SIMD3(e.x, e.y, e.z) | ||
| 58 | return s / (s * s).sum().squareRoot() | ||
| 59 | } | ||
| 60 | |||
| 61 | /// Cylindrical Earth-shadow test. | ||
| 62 | static func sunlit(_ r: SIMD3<Double>, sun: SIMD3<Double>) -> Bool { | ||
| 63 | let along = (r * sun).sum() | ||
| 64 | if along > 0 { return true } | ||
| 65 | let perp = r - along * sun | ||
| 66 | return (perp * perp).sum().squareRoot() > 6378.135 | ||
| 67 | } | ||
| 68 | |||
| 69 | /// Passes above `minAltitude` in the next `hours`, sampled every 20 s. | ||
| 70 | public static func passes(of elements: OrbitalElements, from start: Date, hours: Double = 72, observer: Observer, minAltitude: Double = 10) -> [SatellitePass] { | ||
| 71 | guard let sgp = try? SGP4(elements) else { return [] } | ||
| 72 | lock.lock() | ||
| 73 | defer { lock.unlock() } | ||
| 74 | var out: [SatellitePass] = [] | ||
| 75 | var current: (rise: Date, riseAz: Double, peak: Date, maxAlt: Double, visible: Bool)? | ||
| 76 | var lastAz = 0.0 | ||
| 77 | let step = 20.0 | ||
| 78 | var t = start | ||
| 79 | let end = start.addingTimeInterval(hours * 3600) | ||
| 80 | while t <= end { | ||
| 81 | guard let r = sgp.position(at: t) else { break } | ||
| 82 | let h = lookAngles(teme: r, at: t, from: observer) | ||
| 83 | if h.altitude >= 0 { | ||
| 84 | if current == nil { current = (t, h.azimuth, t, h.altitude, false) } | ||
| 85 | if h.altitude > current!.maxAlt { current!.maxAlt = h.altitude; current!.peak = t } | ||
| 86 | if !current!.visible, h.altitude >= minAltitude { | ||
| 87 | var tt = time(t) | ||
| 88 | let sun = Astronomy_Equator(BODY_SUN, &tt, observer.astro, EQUATOR_OF_DATE, ABERRATION) | ||
| 89 | let sunAlt = Astronomy_Horizon(&tt, observer.astro, sun.ra, sun.dec, REFRACTION_NORMAL).altitude | ||
| 90 | if sunAlt < -6, sunlit(r, sun: sunDirection(t)) { current!.visible = true } | ||
| 91 | } | ||
| 92 | lastAz = h.azimuth | ||
| 93 | } else if let c = current { | ||
| 94 | if c.maxAlt >= minAltitude { | ||
| 95 | 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)) | ||
| 96 | } | ||
| 97 | current = nil | ||
| 98 | } | ||
| 99 | t = t.addingTimeInterval(step) | ||
| 100 | } | ||
| 101 | return out | ||
| 102 | } | ||
| 103 | } | ||
Packages/MoonriseKit/Sources/SkyEngine/SGP4.swift added +232
| @@ -0,0 +1,232 @@ | |||
| 1 | import Foundation | ||
| 2 | |||
| 3 | /// General perturbations orbital elements in CelesTrak's OMM JSON form. | ||
| 4 | public struct OrbitalElements: Codable, Sendable, Hashable { | ||
| 5 | public var OBJECT_NAME: String | ||
| 6 | public var NORAD_CAT_ID: Int | ||
| 7 | /// UTC, ISO 8601 without a zone, e.g. "2026-10-18T03:28:15.076416". | ||
| 8 | public var EPOCH: String | ||
| 9 | /// Revolutions per day. | ||
| 10 | public var MEAN_MOTION: Double | ||
| 11 | public var ECCENTRICITY: Double | ||
| 12 | /// Degrees. | ||
| 13 | public var INCLINATION: Double | ||
| 14 | public var RA_OF_ASC_NODE: Double | ||
| 15 | public var ARG_OF_PERICENTER: Double | ||
| 16 | public var MEAN_ANOMALY: Double | ||
| 17 | public var BSTAR: Double | ||
| 18 | |||
| 19 | public init(name: String, norad: Int, epoch: String, meanMotion: Double, eccentricity: Double, inclination: Double, raan: Double, argPerigee: Double, meanAnomaly: Double, bstar: Double) { | ||
| 20 | OBJECT_NAME = name; NORAD_CAT_ID = norad; EPOCH = epoch; MEAN_MOTION = meanMotion; ECCENTRICITY = eccentricity | ||
| 21 | INCLINATION = inclination; RA_OF_ASC_NODE = raan; ARG_OF_PERICENTER = argPerigee; MEAN_ANOMALY = meanAnomaly; BSTAR = bstar | ||
| 22 | } | ||
| 23 | |||
| 24 | public var epochDate: Date? { | ||
| 25 | let parts = EPOCH.split(separator: ".") | ||
| 26 | guard let whole = try? Date(String(parts[0]) + "Z", strategy: .iso8601) else { return nil } | ||
| 27 | let frac = parts.count > 1 ? Double("0." + parts[1]) ?? 0 : 0 | ||
| 28 | return whole.addingTimeInterval(frac) | ||
| 29 | } | ||
| 30 | } | ||
| 31 | |||
| 32 | /// Near-Earth SGP4 (Vallado et al. 2006, "Revisiting Spacetrack Report #3"), WGS-72 constants. | ||
| 33 | /// Deep-space orbits (period of 225 minutes or more) are not supported. | ||
| 34 | public struct SGP4: Sendable { | ||
| 35 | // WGS-72 | ||
| 36 | static let mu = 398600.8 | ||
| 37 | static let re = 6378.135 | ||
| 38 | static let xke = 60.0 / ((re * re * re) / mu).squareRoot() | ||
| 39 | static let j2 = 0.001082616 | ||
| 40 | static let j3 = -0.00000253881 | ||
| 41 | static let j4 = -0.00000165597 | ||
| 42 | static let j3oj2 = j3 / j2 | ||
| 43 | static let twoPi = 2 * Double.pi | ||
| 44 | |||
| 45 | public let epoch: Date | ||
| 46 | let ecco, inclo, nodeo, argpo, mo, bstar: Double | ||
| 47 | let noUnkozai, ao, con41, x1mth2, x7thm1, eta, cc1, cc4, cc5, d2, d3, d4, t2cof, t3cof, t4cof, t5cof: Double | ||
| 48 | let mdot, argpdot, nodedot, nodecf, omgcof, xmcof, xlcof, aycof, delmo, sinmao, cosio, sinio: Double | ||
| 49 | let isimp: Bool | ||
| 50 | |||
| 51 | public enum Failure: Error { case deepSpace, badEpoch } | ||
| 52 | |||
| 53 | public init(_ el: OrbitalElements) throws { | ||
| 54 | guard let epoch = el.epochDate else { throw Failure.badEpoch } | ||
| 55 | self.epoch = epoch | ||
| 56 | let r = Double.pi / 180 | ||
| 57 | ecco = el.ECCENTRICITY | ||
| 58 | inclo = el.INCLINATION * r | ||
| 59 | nodeo = el.RA_OF_ASC_NODE * r | ||
| 60 | argpo = el.ARG_OF_PERICENTER * r | ||
| 61 | mo = el.MEAN_ANOMALY * r | ||
| 62 | bstar = el.BSTAR | ||
| 63 | let noKozai = el.MEAN_MOTION * Self.twoPi / 1440 | ||
| 64 | let x2o3 = 2.0 / 3.0 | ||
| 65 | |||
| 66 | // initl | ||
| 67 | let eccsq = ecco * ecco, omeosq = 1 - eccsq, rteosq = omeosq.squareRoot() | ||
| 68 | cosio = cos(inclo) | ||
| 69 | let cosio2 = cosio * cosio | ||
| 70 | let ak = pow(Self.xke / noKozai, x2o3) | ||
| 71 | let d1 = 0.75 * Self.j2 * (3 * cosio2 - 1) / (rteosq * omeosq) | ||
| 72 | var del = d1 / (ak * ak) | ||
| 73 | let adel = ak * (1 - del * del - del * (1.0 / 3.0 + 134 * del * del / 81)) | ||
| 74 | del = d1 / (adel * adel) | ||
| 75 | noUnkozai = noKozai / (1 + del) | ||
| 76 | guard Self.twoPi / noUnkozai < 225 else { throw Failure.deepSpace } | ||
| 77 | ao = pow(Self.xke / noUnkozai, x2o3) | ||
| 78 | sinio = sin(inclo) | ||
| 79 | let po = ao * omeosq | ||
| 80 | let con42 = 1 - 5 * cosio2 | ||
| 81 | con41 = -con42 - cosio2 - cosio2 | ||
| 82 | let posq = po * po | ||
| 83 | let rp = ao * (1 - ecco) | ||
| 84 | |||
| 85 | // sgp4init | ||
| 86 | let ss = 78 / Self.re + 1 | ||
| 87 | let qzms2t = pow((120 - 78) / Self.re, 4) | ||
| 88 | isimp = rp < 220 / Self.re + 1 | ||
| 89 | var sfour = ss, qzms24 = qzms2t | ||
| 90 | let perige = (rp - 1) * Self.re | ||
| 91 | if perige < 156 { | ||
| 92 | sfour = perige < 98 ? 20 : perige - 78 | ||
| 93 | qzms24 = pow((120 - sfour) / Self.re, 4) | ||
| 94 | sfour = sfour / Self.re + 1 | ||
| 95 | } | ||
| 96 | let pinvsq = 1 / posq | ||
| 97 | let tsi = 1 / (ao - sfour) | ||
| 98 | eta = ao * ecco * tsi | ||
| 99 | let etasq = eta * eta, eeta = ecco * eta | ||
| 100 | let psisq = abs(1 - etasq) | ||
| 101 | let coef = qzms24 * pow(tsi, 4) | ||
| 102 | let coef1 = coef / pow(psisq, 3.5) | ||
| 103 | let cc2 = coef1 * noUnkozai * (ao * (1 + 1.5 * etasq + eeta * (4 + etasq)) + 0.375 * Self.j2 * tsi / psisq * con41 * (8 + 3 * etasq * (8 + etasq))) | ||
| 104 | cc1 = bstar * cc2 | ||
| 105 | let cc3 = ecco > 1e-4 ? -2 * coef * tsi * Self.j3oj2 * noUnkozai * sinio / ecco : 0 | ||
| 106 | x1mth2 = 1 - cosio2 | ||
| 107 | cc4 = 2 * noUnkozai * coef1 * ao * omeosq * (eta * (2 + 0.5 * etasq) + ecco * (0.5 + 2 * etasq) | ||
| 108 | - Self.j2 * tsi / (ao * psisq) * (-3 * con41 * (1 - 2 * eeta + etasq * (1.5 - 0.5 * eeta)) + 0.75 * x1mth2 * (2 * etasq - eeta * (1 + etasq)) * cos(2 * argpo))) | ||
| 109 | cc5 = 2 * coef1 * ao * omeosq * (1 + 2.75 * (etasq + eeta) + eeta * etasq) | ||
| 110 | let cosio4 = cosio2 * cosio2 | ||
| 111 | let temp1 = 1.5 * Self.j2 * pinvsq * noUnkozai | ||
| 112 | let temp2 = 0.5 * temp1 * Self.j2 * pinvsq | ||
| 113 | let temp3 = -0.46875 * Self.j4 * pinvsq * pinvsq * noUnkozai | ||
| 114 | mdot = noUnkozai + 0.5 * temp1 * rteosq * con41 + 0.0625 * temp2 * rteosq * (13 - 78 * cosio2 + 137 * cosio4) | ||
| 115 | argpdot = -0.5 * temp1 * con42 + 0.0625 * temp2 * (7 - 114 * cosio2 + 395 * cosio4) + temp3 * (3 - 36 * cosio2 + 49 * cosio4) | ||
| 116 | let xhdot1 = -temp1 * cosio | ||
| 117 | nodedot = xhdot1 + (0.5 * temp2 * (4 - 19 * cosio2) + 2 * temp3 * (3 - 7 * cosio2)) * cosio | ||
| 118 | omgcof = bstar * cc3 * cos(argpo) | ||
| 119 | xmcof = ecco > 1e-4 ? -x2o3 * coef * bstar / eeta : 0 | ||
| 120 | nodecf = 3.5 * omeosq * xhdot1 * cc1 | ||
| 121 | t2cof = 1.5 * cc1 | ||
| 122 | xlcof = -0.25 * Self.j3oj2 * sinio * (3 + 5 * cosio) / (abs(cosio + 1) > 1.5e-12 ? 1 + cosio : 1.5e-12) | ||
| 123 | aycof = -0.5 * Self.j3oj2 * sinio | ||
| 124 | delmo = pow(1 + eta * cos(mo), 3) | ||
| 125 | sinmao = sin(mo) | ||
| 126 | x7thm1 = 7 * cosio2 - 1 | ||
| 127 | if !isimp { | ||
| 128 | let cc1sq = cc1 * cc1 | ||
| 129 | d2 = 4 * ao * tsi * cc1sq | ||
| 130 | let temp = d2 * tsi * cc1 / 3 | ||
| 131 | d3 = (17 * ao + sfour) * temp | ||
| 132 | d4 = 0.5 * temp * ao * tsi * (221 * ao + 31 * sfour) * cc1 | ||
| 133 | t3cof = d2 + 2 * cc1sq | ||
| 134 | t4cof = 0.25 * (3 * d3 + cc1 * (12 * d2 + 10 * cc1sq)) | ||
| 135 | t5cof = 0.2 * (3 * d4 + 12 * cc1 * d3 + 6 * d2 * d2 + 15 * cc1sq * (2 * d2 + cc1sq)) | ||
| 136 | } else { | ||
| 137 | d2 = 0; d3 = 0; d4 = 0; t3cof = 0; t4cof = 0; t5cof = 0 | ||
| 138 | } | ||
| 139 | } | ||
| 140 | |||
| 141 | /// Position (km) and velocity (km/s) in the TEME frame, `minutes` after epoch. | ||
| 142 | public func propagate(minutes t: Double) -> (r: SIMD3<Double>, v: SIMD3<Double>)? { | ||
| 143 | let x2o3 = 2.0 / 3.0 | ||
| 144 | let xmdf = mo + mdot * t | ||
| 145 | let argpdf = argpo + argpdot * t | ||
| 146 | let nodedf = nodeo + nodedot * t | ||
| 147 | var argpm = argpdf, mm = xmdf | ||
| 148 | let t2 = t * t | ||
| 149 | var nodem = nodedf + nodecf * t2 | ||
| 150 | var tempa = 1 - cc1 * t | ||
| 151 | var tempe = bstar * cc4 * t | ||
| 152 | var templ = t2cof * t2 | ||
| 153 | if !isimp { | ||
| 154 | let delomg = omgcof * t | ||
| 155 | let delm = xmcof * (pow(1 + eta * cos(xmdf), 3) - delmo) | ||
| 156 | let temp = delomg + delm | ||
| 157 | mm = xmdf + temp | ||
| 158 | argpm = argpdf - temp | ||
| 159 | let t3 = t2 * t, t4 = t3 * t | ||
| 160 | tempa = tempa - d2 * t2 - d3 * t3 - d4 * t4 | ||
| 161 | tempe = tempe + bstar * cc5 * (sin(mm) - sinmao) | ||
| 162 | templ = templ + t3cof * t3 + t4 * (t4cof + t * t5cof) | ||
| 163 | } | ||
| 164 | let am = pow(Self.xke / noUnkozai, x2o3) * tempa * tempa | ||
| 165 | let nm = Self.xke / pow(am, 1.5) | ||
| 166 | var em = ecco - tempe | ||
| 167 | guard em < 1, am >= 0.95 else { return nil } | ||
| 168 | if em < 1e-6 { em = 1e-6 } | ||
| 169 | mm += noUnkozai * templ | ||
| 170 | var xlm = mm + argpm + nodem | ||
| 171 | nodem = fmod(nodem, Self.twoPi) | ||
| 172 | argpm = fmod(argpm, Self.twoPi) | ||
| 173 | xlm = fmod(xlm, Self.twoPi) | ||
| 174 | mm = fmod(xlm - argpm - nodem, Self.twoPi) | ||
| 175 | |||
| 176 | // Long-period periodics | ||
| 177 | let axnl = em * cos(argpm) | ||
| 178 | var temp = 1 / (am * (1 - em * em)) | ||
| 179 | let aynl = em * sin(argpm) + temp * aycof | ||
| 180 | let xl = mm + argpm + nodem + temp * xlcof * axnl | ||
| 181 | |||
| 182 | // Kepler's equation | ||
| 183 | let u = fmod(xl - nodem, Self.twoPi) | ||
| 184 | var eo1 = u, tem5 = 9999.9, ktr = 1 | ||
| 185 | var sineo1 = 0.0, coseo1 = 0.0 | ||
| 186 | while abs(tem5) >= 1e-12 && ktr <= 10 { | ||
| 187 | sineo1 = sin(eo1); coseo1 = cos(eo1) | ||
| 188 | tem5 = 1 - coseo1 * axnl - sineo1 * aynl | ||
| 189 | tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5 | ||
| 190 | if abs(tem5) >= 0.95 { tem5 = tem5 > 0 ? 0.95 : -0.95 } | ||
| 191 | eo1 += tem5 | ||
| 192 | ktr += 1 | ||
| 193 | } | ||
| 194 | |||
| 195 | // Short-period periodics | ||
| 196 | let ecose = axnl * coseo1 + aynl * sineo1 | ||
| 197 | let esine = axnl * sineo1 - aynl * coseo1 | ||
| 198 | let el2 = axnl * axnl + aynl * aynl | ||
| 199 | let pl = am * (1 - el2) | ||
| 200 | guard pl >= 0 else { return nil } | ||
| 201 | let rl = am * (1 - ecose) | ||
| 202 | let rdotl = am.squareRoot() * esine / rl | ||
| 203 | let rvdotl = pl.squareRoot() / rl | ||
| 204 | let betal = (1 - el2).squareRoot() | ||
| 205 | temp = esine / (1 + betal) | ||
| 206 | let sinu = am / rl * (sineo1 - aynl - axnl * temp) | ||
| 207 | let cosu = am / rl * (coseo1 - axnl + aynl * temp) | ||
| 208 | var su = atan2(sinu, cosu) | ||
| 209 | let sin2u = (cosu + cosu) * sinu | ||
| 210 | let cos2u = 1 - 2 * sinu * sinu | ||
| 211 | temp = 1 / pl | ||
| 212 | let temp1 = 0.5 * Self.j2 * temp | ||
| 213 | let temp2 = temp1 * temp | ||
| 214 | let mrt = rl * (1 - 1.5 * temp2 * betal * con41) + 0.5 * temp1 * x1mth2 * cos2u | ||
| 215 | su -= 0.25 * temp2 * x7thm1 * sin2u | ||
| 216 | let xnode = nodem + 1.5 * temp2 * cosio * sin2u | ||
| 217 | let xinc = inclo + 1.5 * temp2 * cosio * sinio * cos2u | ||
| 218 | let mvt = rdotl - nm * temp1 * x1mth2 * sin2u / Self.xke | ||
| 219 | let rvdot = rvdotl + nm * temp1 * (x1mth2 * cos2u + 1.5 * con41) / Self.xke | ||
| 220 | |||
| 221 | let sinsu = sin(su), cossu = cos(su), snod = sin(xnode), cnod = cos(xnode), sini = sin(xinc), cosi = cos(xinc) | ||
| 222 | let xmx = -snod * cosi, xmy = cnod * cosi | ||
| 223 | let uv = SIMD3(xmx * sinsu + cnod * cossu, xmy * sinsu + snod * cossu, sini * sinsu) | ||
| 224 | let vv = SIMD3(xmx * cossu - cnod * sinsu, xmy * cossu - snod * sinsu, sini * cossu) | ||
| 225 | let vkmpersec = Self.re * Self.xke / 60 | ||
| 226 | return (mrt * uv * Self.re, (mvt * uv + rvdot * vv) * vkmpersec) | ||
| 227 | } | ||
| 228 | |||
| 229 | public func position(at date: Date) -> SIMD3<Double>? { | ||
| 230 | propagate(minutes: date.timeIntervalSince(epoch) / 60)?.r | ||
| 231 | } | ||
| 232 | } | ||
Packages/MoonriseKit/Tests/SkyEngineTests/SGP4Tests.swift added +62
| @@ -0,0 +1,62 @@ | |||
| 1 | import Foundation | ||
| 2 | import Testing | ||
| 3 | @testable import SkyEngine | ||
| 4 | |||
| 5 | // Vallado's SGP4 verification case for satellite 00005 (SGP4-VER.TLE / tcppver.out): | ||
| 6 | // 1 00005U 58002B 00179.78495062 .00000023 00000-0 28098-4 0 4753 | ||
| 7 | // 2 00005 34.2682 348.7242 1859667 331.7664 19.3264 10.82419157413667 | ||
| 8 | // Epoch 2000 day 179.78495062 = 2000-06-27T18:50:19.733568Z. | ||
| 9 | let vanguard = OrbitalElements(name: "00005", norad: 5, epoch: "2000-06-27T18:50:19.733568", | ||
| 10 | meanMotion: 10.82419157, eccentricity: 0.1859667, inclination: 34.2682, | ||
| 11 | raan: 348.7242, argPerigee: 331.7664, meanAnomaly: 19.3264, bstar: 0.28098e-4) | ||
| 12 | |||
| 13 | func dist(_ a: SIMD3<Double>, _ b: SIMD3<Double>) -> Double { let d = a - b; return (d * d).sum().squareRoot() } | ||
| 14 | |||
| 15 | @Test(arguments: [ | ||
| 16 | (0.0, SIMD3(7022.46529266, -1400.08296755, 0.03995155), SIMD3(1.893841015, 6.405893759, 4.534807250)), | ||
| 17 | (360.0, SIMD3(-7154.03120202, -3783.17682504, -3536.19412294), SIMD3(4.741887409, -4.151817765, -2.093935425)), | ||
| 18 | (720.0, SIMD3(-7134.59340119, 6531.68641334, 3260.27186483), SIMD3(-4.113793027, -2.911922039, -2.557327851)), | ||
| 19 | ]) | ||
| 20 | func matchesValladoVerification(_ t: Double, _ r: SIMD3<Double>, _ v: SIMD3<Double>) throws { | ||
| 21 | let s = try SGP4(vanguard) | ||
| 22 | let out = try #require(s.propagate(minutes: t)) | ||
| 23 | #expect(dist(out.r, r) < 0.01, "t=\(t): r off by \(dist(out.r, r)) km") | ||
| 24 | #expect(dist(out.v, v) < 1e-5, "t=\(t): v off by \(dist(out.v, v)) km/s") | ||
| 25 | } | ||
| 26 | |||
| 27 | @Test func epochParsesFractionalSeconds() throws { | ||
| 28 | let d = try #require(vanguard.epochDate) | ||
| 29 | #expect(abs(d.timeIntervalSince1970 - 962_131_819.733568) < 1e-3) | ||
| 30 | } | ||
| 31 | |||
| 32 | @Test func issIsInLowOrbit() throws { | ||
| 33 | // CelesTrak GP for the ISS, 2026-10-18. | ||
| 34 | let iss = OrbitalElements(name: "ISS (ZARYA)", norad: 25544, epoch: "2026-10-18T03:28:15.076416", | ||
| 35 | meanMotion: 15.48918166, eccentricity: 0.000678, inclination: 51.6328, | ||
| 36 | raan: 51.6259, argPerigee: 277.7856, meanAnomaly: 82.2363, bstar: 0.00014456) | ||
| 37 | let s = try SGP4(iss) | ||
| 38 | let r = try #require(s.propagate(minutes: 60)?.r) | ||
| 39 | let altitude = (r * r).sum().squareRoot() - 6378.135 | ||
| 40 | #expect(altitude > 380 && altitude < 450, "\(altitude) km") | ||
| 41 | } | ||
| 42 | |||
| 43 | @Test func gmstAtJ2000() { | ||
| 44 | // GMST at 2000-01-01 12:00 UT1 is 280.46061837 degrees. | ||
| 45 | let g = SkyEngine.gmst(Date(timeIntervalSince1970: 946_728_000)) * 180 / .pi | ||
| 46 | #expect(abs(g - 280.46061837) < 1e-4) | ||
| 47 | } | ||
| 48 | |||
| 49 | @Test func issPassesOverChicago() throws { | ||
| 50 | let iss = OrbitalElements(name: "ISS (ZARYA)", norad: 25544, epoch: "2026-10-18T03:28:15.076416", | ||
| 51 | meanMotion: 15.48918166, eccentricity: 0.000678, inclination: 51.6328, | ||
| 52 | raan: 51.6259, argPerigee: 277.7856, meanAnomaly: 82.2363, bstar: 0.00014456) | ||
| 53 | let start = try #require(iss.epochDate) | ||
| 54 | let passes = SkyEngine.passes(of: iss, from: start, hours: 72, observer: chicago) | ||
| 55 | // The ISS (51.6° inclination) passes over 42° N several times a day. | ||
| 56 | #expect(passes.count >= 6, "\(passes.count) passes") | ||
| 57 | for p in passes { | ||
| 58 | #expect(p.maxAltitude >= 10 && p.maxAltitude <= 90) | ||
| 59 | let minutes = p.set.timeIntervalSince(p.rise) / 60 | ||
| 60 | #expect(minutes > 2 && minutes < 12, "\(minutes) min") | ||
| 61 | } | ||
| 62 | } | ||
server/src/celestrak.ts added +37
| @@ -0,0 +1,37 @@ | |||
| 1 | // CelesTrak GP data. CelesTrak updates every 2 hours and blocks IPs that fetch more | ||
| 2 | // often or keep requesting after errors, so poll every 2 hours and stop on any non-200. | ||
| 3 | |||
| 4 | /** ISS, Tiangong (CSS Tianhe) and Hubble. */ | ||
| 5 | export const SATELLITES = [25544, 48274, 20580]; | ||
| 6 | |||
| 7 | export interface OMM { | ||
| 8 | OBJECT_NAME: string; | ||
| 9 | NORAD_CAT_ID: number; | ||
| 10 | EPOCH: string; | ||
| 11 | MEAN_MOTION: number; | ||
| 12 | ECCENTRICITY: number; | ||
| 13 | INCLINATION: number; | ||
| 14 | RA_OF_ASC_NODE: number; | ||
| 15 | ARG_OF_PERICENTER: number; | ||
| 16 | MEAN_ANOMALY: number; | ||
| 17 | BSTAR: number; | ||
| 18 | } | ||
| 19 | |||
| 20 | const KEYS: (keyof OMM)[] = ["OBJECT_NAME", "NORAD_CAT_ID", "EPOCH", "MEAN_MOTION", "ECCENTRICITY", "INCLINATION", "RA_OF_ASC_NODE", "ARG_OF_PERICENTER", "MEAN_ANOMALY", "BSTAR"]; | ||
| 21 | |||
| 22 | export function pick(rows: Record<string, unknown>[]): OMM[] { | ||
| 23 | return rows.map((r) => Object.fromEntries(KEYS.map((k) => [k, r[k]])) as unknown as OMM); | ||
| 24 | } | ||
| 25 | |||
| 26 | /** Fetches each satellite in turn; throws on the first non-200 so nothing more is requested. */ | ||
| 27 | export async function fetchElements(): Promise<OMM[]> { | ||
| 28 | const out: OMM[] = []; | ||
| 29 | for (const id of SATELLITES) { | ||
| 30 | const res = await fetch(`https://celestrak.org/NORAD/elements/gp.php?CATNR=${id}&FORMAT=JSON`, { | ||
| 31 | headers: { "user-agent": "moonrise-server (sh.krz.moonrise)" }, | ||
| 32 | }); | ||
| 33 | if (res.status !== 200) throw new Error(`CelesTrak ${res.status} for ${id}`); | ||
| 34 | out.push(...pick((await res.json()) as Record<string, unknown>[])); | ||
| 35 | } | ||
| 36 | return out; | ||
| 37 | } | ||
server/src/index.ts +18
| @@ -1,5 +1,6 @@ | |||
| 1 | import { send, type ApnsConfig } from "./apns"; | 1 | import { send, type ApnsConfig } from "./apns"; |
| 2 | import { actionableKp, alertBands, eventKey, isNight } from "./aurora"; | 2 | import { actionableKp, alertBands, eventKey, isNight } from "./aurora"; |
| 3 | import { fetchElements } from "./celestrak"; | ||
| 3 | import { fetchUpcoming, type Launch } from "./ll2"; | 4 | import { fetchUpcoming, type Launch } from "./ll2"; |
| 4 | import { PRIVACY_HTML } from "./privacy"; | 5 | import { PRIVACY_HTML } from "./privacy"; |
| 5 | import { fetchSpaceWeather, next24h, type SpaceWeather } from "./swpc"; | 6 | import { fetchSpaceWeather, next24h, type SpaceWeather } from "./swpc"; |
| @@ -23,6 +24,7 @@ interface DeviceRequest { | |||
| 23 | } | 24 | } |
| 24 | 25 | ||
| 25 | const LAUNCH_CRON = "*/10 * * * *"; | 26 | const LAUNCH_CRON = "*/10 * * * *"; |
| 27 | const SATELLITE_CRON = "7 */2 * * *"; | ||
| 26 | 28 | ||
| 27 | const json = (body: unknown, status = 200) => Response.json(body, { status }); | 29 | const json = (body: unknown, status = 200) => Response.json(body, { status }); |
| 28 | 30 | ||
| @@ -58,6 +60,10 @@ export default { | |||
| 58 | return new Response(PRIVACY_HTML, { headers: { "content-type": "text/html; charset=utf-8", "cache-control": "public, max-age=3600" } }); | 60 | return new Response(PRIVACY_HTML, { headers: { "content-type": "text/html; charset=utf-8", "cache-control": "public, max-age=3600" } }); |
| 59 | if (request.method === "GET" && url.pathname === "/v1/today") return today(env, now); | 61 | if (request.method === "GET" && url.pathname === "/v1/today") return today(env, now); |
| 60 | if (request.method === "GET" && url.pathname === "/v1/launches") return launches(env, now); | 62 | if (request.method === "GET" && url.pathname === "/v1/launches") return launches(env, now); |
| 63 | if (request.method === "GET" && url.pathname === "/v1/satellites") { | ||
| 64 | const row = await env.DB.prepare("SELECT fetched_at, body FROM feed_snapshots WHERE feed = 'celestrak' AND ok = 1").first<{ fetched_at: string; body: string }>(); | ||
| 65 | return json({ fetchedAt: row?.fetched_at ?? null, satellites: row ? JSON.parse(row.body) : [] }); | ||
| 66 | } | ||
| 61 | 67 | ||
| 62 | if (request.method === "POST" && url.pathname === "/v1/devices") { | 68 | if (request.method === "POST" && url.pathname === "/v1/devices") { |
| 63 | const d = (await request.json().catch(() => ({}))) as Partial<DeviceRequest>; | 69 | const d = (await request.json().catch(() => ({}))) as Partial<DeviceRequest>; |
| @@ -85,6 +91,7 @@ export default { | |||
| 85 | const now = new Date(controller.scheduledTime); | 91 | const now = new Date(controller.scheduledTime); |
| 86 | // Launches every 10 minutes: 6 of the 15 anonymous Launch Library 2 requests per hour. | 92 | // Launches every 10 minutes: 6 of the 15 anonymous Launch Library 2 requests per hour. |
| 87 | if (controller.cron === LAUNCH_CRON) ctx.waitUntil(pollLaunches(env, now)); | 93 | if (controller.cron === LAUNCH_CRON) ctx.waitUntil(pollLaunches(env, now)); |
| 94 | else if (controller.cron === SATELLITE_CRON) ctx.waitUntil(pollSatellites(env, now)); | ||
| 88 | else ctx.waitUntil(pollSpaceWeather(env, now)); | 95 | else ctx.waitUntil(pollSpaceWeather(env, now)); |
| 89 | }, | 96 | }, |
| 90 | } satisfies ExportedHandler<Env>; | 97 | } satisfies ExportedHandler<Env>; |
| @@ -142,6 +149,17 @@ async function pollLaunches(env: Env, now: Date): Promise<void> { | |||
| 142 | ]); | 149 | ]); |
| 143 | } | 150 | } |
| 144 | 151 | ||
| 152 | async function pollSatellites(env: Env, now: Date): Promise<void> { | ||
| 153 | try { | ||
| 154 | const elements = await fetchElements(); | ||
| 155 | await env.DB.prepare("INSERT INTO feed_snapshots (feed, fetched_at, ok, body) VALUES ('celestrak', ?1, 1, ?2) ON CONFLICT (feed) DO UPDATE SET fetched_at = ?1, ok = 1, body = ?2") | ||
| 156 | .bind(now.toISOString(), JSON.stringify(elements)).run(); | ||
| 157 | } catch (e) { | ||
| 158 | console.error("celestrak", String(e)); | ||
| 159 | await env.DB.prepare("INSERT INTO feed_snapshots (feed, fetched_at, ok) VALUES ('celestrak:error', ?1, 0) ON CONFLICT (feed) DO UPDATE SET fetched_at = ?1").bind(now.toISOString()).run(); | ||
| 160 | } | ||
| 161 | } | ||
| 162 | |||
| 145 | async function pollSpaceWeather(env: Env, now: Date): Promise<void> { | 163 | async function pollSpaceWeather(env: Env, now: Date): Promise<void> { |
| 146 | let w: SpaceWeather; | 164 | let w: SpaceWeather; |
| 147 | try { | 165 | try { |
server/test/celestrak.test.ts added +16
| @@ -0,0 +1,16 @@ | |||
| 1 | import assert from "node:assert/strict"; | ||
| 2 | import { test } from "node:test"; | ||
| 3 | import { pick } from "../src/celestrak.ts"; | ||
| 4 | |||
| 5 | test("keeps only the SGP4 fields", () => { | ||
| 6 | // CelesTrak GP JSON for the ISS, 2026-10-18. | ||
| 7 | const row = { | ||
| 8 | OBJECT_NAME: "ISS (ZARYA)", OBJECT_ID: "1998-067A", EPOCH: "2026-10-18T03:28:15.076416", MEAN_MOTION: 15.48918166, | ||
| 9 | ECCENTRICITY: 0.000678, INCLINATION: 51.6328, RA_OF_ASC_NODE: 51.6259, ARG_OF_PERICENTER: 277.7856, MEAN_ANOMALY: 82.2363, | ||
| 10 | EPHEMERIS_TYPE: 0, CLASSIFICATION_TYPE: "U", NORAD_CAT_ID: 25544, ELEMENT_SET_NO: 999, REV_AT_EPOCH: 59081, BSTAR: 0.00014456, | ||
| 11 | MEAN_MOTION_DOT: 7.503e-05, MEAN_MOTION_DDOT: 0, | ||
| 12 | }; | ||
| 13 | const [o] = pick([row]); | ||
| 14 | assert.deepEqual(Object.keys(o).sort(), ["ARG_OF_PERICENTER", "BSTAR", "ECCENTRICITY", "EPOCH", "INCLINATION", "MEAN_ANOMALY", "MEAN_MOTION", "NORAD_CAT_ID", "OBJECT_NAME", "RA_OF_ASC_NODE"]); | ||
| 15 | assert.equal(o.NORAD_CAT_ID, 25544); | ||
| 16 | }); | ||
server/wrangler.jsonc +1 −1
| @@ -2,7 +2,7 @@ | |||
| 2 | "name": "moonrise-server", | 2 | "name": "moonrise-server", |
| 3 | "main": "src/index.ts", | 3 | "main": "src/index.ts", |
| 4 | "compatibility_date": "2026-10-01", | 4 | "compatibility_date": "2026-10-01", |
| 5 | "triggers": { "crons": ["*/5 * * * *", "*/10 * * * *"] }, | 5 | "triggers": { "crons": ["*/5 * * * *", "*/10 * * * *", "7 */2 * * *"] }, |
| 6 | "d1_databases": [ | 6 | "d1_databases": [ |
| 7 | { | 7 | { |
| 8 | "binding": "DB", | 8 | "binding": "DB", |