Packages/MoonriseKit/Sources/SkyEngine/Geomagnetic.swift

 1import Foundation
 2
 3/// Geomagnetic latitude from a centred-dipole model of Earth's field.
 4public enum Geomagnetic {
 5    /// North geomagnetic pole, IGRF-14 for 2025 (approximately 80.8° N, 72.8° W).
 6    public static let poleLatitude = 80.8
 7    public static let poleLongitude = -72.8
 8
 9    /// Geomagnetic latitude in degrees for a geographic position.
10    public static func latitude(_ observer: Observer) -> Double {
11        let r = Double.pi / 180
12        let s = sin(observer.latitude * r) * sin(poleLatitude * r)
13            + cos(observer.latitude * r) * cos(poleLatitude * r) * cos((observer.longitude - poleLongitude) * r)
14        return asin(max(-1, min(1, s))) / r
15    }
16
17    /// The backend's 5° band, named by its equatorward edge: 50 covers 50–55° N, -55 covers 50–55° S.
18    /// Must match `bandFor` in server/src/aurora.ts.
19    public static func band(_ observer: Observer) -> Int {
20        let m = latitude(observer)
21        let b = Int((abs(m) / 5).rounded(.down)) * 5
22        return m < 0 ? -(b + 5) : b
23    }
24
25    /// Geomagnetic latitude of the auroral oval's equatorward edge for a Kp value.
26    /// Must match `ovalEdge` in server/src/aurora.ts.
27    public static func ovalEdge(kp: Double) -> Double {
28        66.5 - 2.04 * kp
29    }
30
31    /// Lowest geomagnetic latitude (absolute) where aurora may be seen low on the poleward horizon.
32    /// Must match `visibleFrom` in server/src/aurora.ts.
33    public static func visibleFrom(kp: Double) -> Double {
34        ovalEdge(kp: kp) - 6
35    }
36
37    /// Kp at which aurora may be seen on the horizon from here; 0 under or poleward of the quiet oval.
38    public static func kpNeeded(_ observer: Observer) -> Double {
39        kpNeeded(geomagneticLatitude: latitude(observer))
40    }
41
42    /// `kpNeeded(_:)` for a geomagnetic latitude in degrees.
43    public static func kpNeeded(geomagneticLatitude: Double) -> Double {
44        max(0, (visibleFrom(kp: 0) - abs(geomagneticLatitude)) / 2.04)
45    }
46}