Speed up propagation step time by 3x

python/satkit/__init__.py

@@ -1,8 +1,5 @@
-from importlib.metadata import version
-
-__version__ = version("satkit")
-
 from .satkit import *  # type: ignore
+from .satkit import __version__
 
 # The core data (IERS nutation tables, gravity models to degree 70) is
 # compiled into the extension, so satkit works with no data directory at all.

python/src/lib.rs

@@ -176,6 +176,7 @@ fn frametransform(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
 
 #[pymodule]
 pub fn satkit(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
+    m.add("__version__", env!("CARGO_PKG_VERSION"))?;
     m.add_class::<PyInstant>()?;
     m.add_class::<PyDuration>()?;
     m.add_class::<pyinstant::PyTimeScale>()?;

src/earthgravity.rs

@@ -432,10 +432,40 @@ impl Gravity {
         )
     }
 
+    // On baseline x86-64, `f64::mul_add` is a call into the `fma` runtime
+    // function per term, so the kernels are also compiled with the `fma`
+    // feature and dispatched at runtime.
+
     fn accel_and_partials_t<const N: usize, const NP4: usize>(
         &self,
         pos: &Vector3,
         max_order: usize,
+    ) -> (Vector3, Matrix3) {
+        #[cfg(target_arch = "x86_64")]
+        {
+            if std::arch::is_x86_feature_detected!("fma") {
+                // SAFETY: the `fma` feature was detected on this CPU.
+                return unsafe { self.accel_and_partials_t_fma::<N, NP4>(pos, max_order) };
+            }
+        }
+        self.accel_and_partials_t_inner::<N, NP4>(pos, max_order)
+    }
+
+    #[cfg(target_arch = "x86_64")]
+    #[target_feature(enable = "fma")]
+    unsafe fn accel_and_partials_t_fma<const N: usize, const NP4: usize>(
+        &self,
+        pos: &Vector3,
+        max_order: usize,
+    ) -> (Vector3, Matrix3) {
+        self.accel_and_partials_t_inner::<N, NP4>(pos, max_order)
+    }
+
+    #[inline(always)]
+    fn accel_and_partials_t_inner<const N: usize, const NP4: usize>(
+        &self,
+        pos: &Vector3,
+        max_order: usize,
     ) -> (Vector3, Matrix3) {
         let (v, w) = self.compute_legendre::<NP4>(pos);
         let accel = self.accel_from_legendre_t::<N, NP4>(&v, &w, max_order);
@@ -448,12 +478,38 @@ impl Gravity {
         pos: &Vector3,
         max_order: usize,
     ) -> Vector3 {
-        let (v, w) = self.compute_legendre::<NP4>(pos);
+        #[cfg(target_arch = "x86_64")]
+        {
+            if std::arch::is_x86_feature_detected!("fma") {
+                // SAFETY: the `fma` feature was detected on this CPU.
+                return unsafe { self.accel_t_fma::<N, NP4>(pos, max_order) };
+            }
+        }
+        self.accel_t_inner::<N, NP4>(pos, max_order)
+    }
 
+    #[cfg(target_arch = "x86_64")]
+    #[target_feature(enable = "fma")]
+    unsafe fn accel_t_fma<const N: usize, const NP4: usize>(
+        &self,
+        pos: &Vector3,
+        max_order: usize,
+    ) -> Vector3 {
+        self.accel_t_inner::<N, NP4>(pos, max_order)
+    }
+
+    #[inline(always)]
+    fn accel_t_inner<const N: usize, const NP4: usize>(
+        &self,
+        pos: &Vector3,
+        max_order: usize,
+    ) -> Vector3 {
+        let (v, w) = self.compute_legendre::<NP4>(pos);
         self.accel_from_legendre_t::<N, NP4>(&v, &w, max_order)
     }
 
     // Equations 7.65 to 7.69 in Montenbruck & Gill
+    #[inline(always)]
     fn partials_from_legendre_t<const N: usize, const NP4: usize>(
         &self,
         v: &Legendre<NP4>,
@@ -556,6 +612,7 @@ impl Gravity {
     }
 
     /// See Equation 3.33 in Montenbruck & Gill
+    #[inline(always)]
     fn accel_from_legendre_t<const N: usize, const NP4: usize>(
         &self,
         v: &Legendre<NP4>,
@@ -601,6 +658,7 @@ impl Gravity {
         numeris::vector![ax, ay, az] * self.gravity_constant / self.radius / self.radius
     }
 
+    #[inline(always)]
     fn compute_legendre<const NP4: usize>(&self, pos: &Vector3) -> (Legendre<NP4>, Legendre<NP4>) {
         let rsq = pos.norm_squared();
         let scale = self.radius / rsq;
@@ -695,14 +753,11 @@ impl Gravity {
         let mut gravity_constant: f64 = 0.0;
         let mut radius: f64 = 0.0;
         let mut max_degree: usize = 0;
-        let mut header_cnt = 0;
 
-        let lines: Vec<&str> = text.lines().collect();
+        let mut lines = text.lines();
 
         // Read header lines
-        for line in &lines {
-            header_cnt += 1;
-
+        for line in lines.by_ref() {
             let s: Vec<&str> = line.split_whitespace().collect();
             // Check for the header terminator before the two-token guard: the
             // ICGEM spec allows a bare "end_of_head" line (no ==== filler),
@@ -737,34 +792,29 @@ impl Gravity {
         let table_dim = (max_degree + 1).min(MAX_COEFF_DIM);
         let mut cs: CoeffTable = CoeffTable::zeros(table_dim, table_dim);
 
-        for line in &lines[header_cnt..] {
-            let s: Vec<&str> = line.split_whitespace().collect();
-            // Need at least keyword, degree, order, and the C coefficient
-            // (index 3); the S coefficient (index 4) is required only when m > 0.
-            if s.len() < 4 {
-                return Err(Error::InvalidLine((*line).to_string()));
-            }
-
-            let n: usize = s[1].parse()?;
-            let m: usize = s[2].parse()?;
+        for line in lines {
+            let invalid = || Error::InvalidLine(line.to_string());
+            // Need at least keyword, degree, order, and the C coefficient;
+            // the S coefficient is required only when m > 0. The tokens are
+            // read one at a time so the lines beyond the stored degree, most
+            // of a full-resolution file, cost only two integer parses.
+            let mut s = line.split_whitespace().skip(1);
+            let n: usize = s.next().ok_or_else(invalid)?.parse()?;
+            let m: usize = s.next().ok_or_else(invalid)?.parse()?;
             // The gfc format requires order <= degree; a violating line would
             // index outside the triangular layout below (panicking for large m,
             // silently aliasing another coefficient for moderate m).
             if m > n {
-                return Err(Error::InvalidLine((*line).to_string()));
+                return Err(invalid());
             }
+            let c = s.next().ok_or_else(invalid)?;
             // Skip coefficients beyond the stored/evaluated degree.
             if n >= table_dim {
                 continue;
             }
-            let v1: f64 = s[3].parse()?;
-            cs[(n, m)] = v1;
+            cs[(n, m)] = c.parse()?;
             if m > 0 {
-                if s.len() < 5 {
-                    return Err(Error::InvalidLine((*line).to_string()));
-                }
-                let v2: f64 = s[4].parse()?;
-                cs[(m - 1, n)] = v2;
+                cs[(m - 1, n)] = s.next().ok_or_else(invalid)?.parse()?;
             }
         }
 

src/itrfcoord.rs

@@ -278,25 +278,19 @@ impl ITRFCoord {
         const E2: f64 = 1.0 - (1.0 - WGS84_F) * (1.0 - WGS84_F);
         const EP2: f64 = E2 / (1.0 - E2);
 
+        // One refinement of the reduced latitude reaches double precision
+        // from 20 km below the surface to beyond GEO.
         let rho = self.itrf[0].hypot(self.itrf[1]);
-        let mut beta: f64 = f64::atan2(self.itrf[2], (1.0 - WGS84_F) * rho);
-        let mut sinbeta: f64 = beta.sin();
-        let mut cosbeta: f64 = beta.cos();
-        let mut phi: f64 = f64::atan2(
-            (B * EP2).mul_add(sinbeta.powi(3), self.itrf[2]),
-            (WGS84_A * E2).mul_add(-cosbeta.powi(3), rho),
+        let beta: f64 = f64::atan2(self.itrf[2], (1.0 - WGS84_F) * rho);
+        let phi: f64 = f64::atan2(
+            (B * EP2).mul_add(beta.sin().powi(3), self.itrf[2]),
+            (WGS84_A * E2).mul_add(-beta.cos().powi(3), rho),
+        );
+        let beta: f64 = f64::atan2((1.0 - WGS84_F) * phi.sin(), phi.cos());
+        let phi: f64 = f64::atan2(
+            (B * EP2).mul_add(beta.sin().powi(3), self.itrf[2]),
+            (WGS84_A * E2).mul_add(-beta.cos().powi(3), rho),
         );
-        let mut betanew: f64 = f64::atan2((1.0 - WGS84_F) * phi.sin(), phi.cos());
-        for _x in 0..5 {
-            beta = betanew;
-            sinbeta = beta.sin();
-            cosbeta = beta.cos();
-            phi = f64::atan2(
-                (B * EP2).mul_add(sinbeta.powi(3), self.itrf[2]),
-                (WGS84_A * E2).mul_add(-cosbeta.powi(3), rho),
-            );
-            betanew = f64::atan2((1.0 - WGS84_F) * phi.sin(), phi.cos());
-        }
         let lat: f64 = phi;
         let lon: f64 = f64::atan2(self.itrf[1], self.itrf[0]);
         let sinphi: f64 = phi.sin();

src/nrlmsise.rs

@@ -3388,6 +3388,8 @@ struct NrlmsiseState {
     s2tloc: f64,
     s3tloc: f64,
     c3tloc: f64,
+    clong: f64,
+    slong: f64,
     apdf: f64,
     apt: [f64; 4],
 }
@@ -3419,6 +3421,8 @@ impl NrlmsiseState {
             s2tloc: 0.0,
             s3tloc: 0.0,
             c3tloc: 0.0,
+            clong: 0.0,
+            slong: 0.0,
             apdf: 0.0,
             apt: [0.0; 4],
         }
@@ -3814,17 +3818,11 @@ fn sg0(ex: f64, p: &[f64], ap: &[f64]) -> f64 {
         / sumex(ex)
 }
 
-fn globe7(
-    p: &[f64],
-    input: &NrlmsiseInput,
-    flags: &NrlmsiseFlags,
-    state: &mut NrlmsiseState,
-) -> f64 {
-    let sr = 7.2722E-5_f64;
+/// Terms of `globe7` and `glob7s` that depend only on the input, computed
+/// once per density evaluation rather than in each of their calls.
+fn globe_setup(input: &NrlmsiseInput, flags: &NrlmsiseFlags, state: &mut NrlmsiseState) {
     let dgtr = 1.74533E-2_f64;
-    let dr = 1.72142E-2_f64;
     let hr = 0.2618_f64;
-    let mut t = [0.0_f64; 15];
     let tloc = input.lst;
     let c = (input.g_lat * dgtr).sin();
     let s = (input.g_lat * dgtr).cos();
@@ -3861,6 +3859,22 @@ fn globe7(
         state.s3tloc = (3.0 * hr * tloc).sin();
         state.c3tloc = (3.0 * hr * tloc).cos();
     }
+    state.clong = (dgtr * input.g_long).cos();
+    state.slong = (dgtr * input.g_long).sin();
+}
+
+fn globe7(
+    p: &[f64],
+    input: &NrlmsiseInput,
+    flags: &NrlmsiseFlags,
+    state: &mut NrlmsiseState,
+) -> f64 {
+    let sr = 7.2722E-5_f64;
+    let dgtr = 1.74533E-2_f64;
+    let dr = 1.72142E-2_f64;
+    let hr = 0.2618_f64;
+    let mut t = [0.0_f64; 15];
+    let tloc = input.lst;
     let cd32 = (dr * (input.doy as f64 - p[31])).cos();
     let cd18 = (2.0 * dr * (input.doy as f64 - p[17])).cos();
     let cd14 = (dr * (input.doy as f64 - p[13])).cos();
@@ -3975,7 +3989,7 @@ fn globe7(
                             + p[110] * state.plg[1][3]
                             + p[111] * state.plg[1][5])
                         * cd14)
-                    * (dgtr * input.g_long).cos()
+                    * state.clong
                     + (p[90] * state.plg[1][2]
                         + p[91] * state.plg[1][4]
                         + p[92] * state.plg[1][6]
@@ -3987,7 +4001,7 @@ fn globe7(
                                 + p[113] * state.plg[1][3]
                                 + p[114] * state.plg[1][5])
                             * cd14)
-                        * (dgtr * input.g_long).sin());
+                        * state.slong);
         }
         if flags.sw[12] != 0.0 {
             t[11] = (1.0 + p[95] * state.plg[0][1])
@@ -4061,7 +4075,6 @@ fn globe7(
 
 fn glob7s(p: &[f64], input: &NrlmsiseInput, flags: &NrlmsiseFlags, state: &NrlmsiseState) -> f64 {
     let dr = 1.72142E-2_f64;
-    let dgtr = 1.74533E-2_f64;
     let _hr = 0.2618_f64;
     let mut t = [0.0_f64; 14];
     let cd32 = (dr * (input.doy as f64 - p[31])).cos();
@@ -4115,14 +4128,14 @@ fn glob7s(p: &[f64], input: &NrlmsiseInput, flags: &NrlmsiseFlags, state: &Nrlms
                 + p[74] * state.plg[1][1]
                 + p[75] * state.plg[1][3]
                 + p[76] * state.plg[1][5])
-                * (dgtr * input.g_long).cos()
+                * state.clong
                 + (p[90] * state.plg[1][2]
                     + p[91] * state.plg[1][4]
                     + p[92] * state.plg[1][6]
                     + p[77] * state.plg[1][1]
                     + p[78] * state.plg[1][3]
                     + p[79] * state.plg[1][5])
-                    * (dgtr * input.g_long).sin());
+                    * state.slong);
     }
     let mut tt = 0.0_f64;
     for i in 0..14 {
@@ -4715,6 +4728,7 @@ fn gtd7(
     let zn2 = [72.5, 55.0, 45.0, 32.5_f64];
     let zmix = 62.5_f64;
     tselec(flags);
+    globe_setup(input, flags, state);
     let xlat = if flags.sw[2] == 0.0 {
         45.0
     } else {
@@ -5012,20 +5026,24 @@ pub fn nrlmsise(
                     _ if r.f10p7_obs_c81 >= 0.0 => r.f10p7_obs_c81,
                     _ => r.f10p7_obs,
                 };
-                match today.map(|t| t.ap_avg) {
+                match today.as_ref().map(|t| t.ap_avg) {
                     Some(a) if a >= 0 => ap = a as f64,
                     _ if r.ap_avg >= 0 => ap = r.ap_avg as f64,
                     _ => {}
                 }
                 // Record for exactly the UTC day `n` days back (spaceweather::get
                 // returns the most recent *prior* record for a missing day, which
-                // must not masquerade as that day's 3-hourly values).
+                // must not masquerade as that day's 3-hourly values). The
+                // current and previous day reuse the records already fetched.
                 if let Ok(day0) = Instant::from_date(year, mon, day) {
                     ap_a = ap_history(sec_of_day, |n| {
                         let d = day0 - Duration::from_days(n as f64);
-                        spaceweather::get(&d)
-                            .ok()
-                            .filter(|rec| (rec.date - d).as_days().abs() < 0.5)
+                        let record = match n {
+                            0 => today.clone(),
+                            1 => Some(r.clone()),
+                            _ => spaceweather::get(&d).ok(),
+                        };
+                        record.filter(|rec| (rec.date - d).as_days().abs() < 0.5)
                     });
                 }
             } else if let Some(predicted) = solar_cycle_forecast::get_predicted_f107(&time) {

src/spaceweather.rs

@@ -277,8 +277,12 @@ fn ensure_default_loaded() {
 /// * Space weather is updated daily in a file: SW-All.csv
 pub fn get<T: TimeLike>(tm: &T) -> Result<SpaceWeatherRecord> {
     let tm = tm.as_instant();
-    ensure_default_loaded();
-    let guard = SPACE_WEATHER.read();
+    let mut guard = SPACE_WEATHER.read();
+    if guard.is_none() {
+        drop(guard);
+        ensure_default_loaded();
+        guard = SPACE_WEATHER.read();
+    }
     let sw = guard.as_ref().ok_or(Error::NoRecordForDate)?;
     // Guard empty data (e.g. a header-only CSV) so the indexing below can't
     // panic; treat it the same as "not loaded".