From 6097fdeddfe7aaf769fe67c80be679a9596230a4 Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Thu, 10 Sep 2026 18:05:45 -0600 Subject: [PATCH 1/6] Add FMA dispatching to speed up gravity calcs by 5x --- src/earthgravity.rs | 60 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 59 insertions(+), 1 deletion(-) diff --git a/src/earthgravity.rs b/src/earthgravity.rs index 58c8be01..9946423f 100644 --- a/src/earthgravity.rs +++ b/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( &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::(pos, max_order) }; + } + } + self.accel_and_partials_t_inner::(pos, max_order) + } + + #[cfg(target_arch = "x86_64")] + #[target_feature(enable = "fma")] + unsafe fn accel_and_partials_t_fma( + &self, + pos: &Vector3, + max_order: usize, + ) -> (Vector3, Matrix3) { + self.accel_and_partials_t_inner::(pos, max_order) + } + + #[inline(always)] + fn accel_and_partials_t_inner( + &self, + pos: &Vector3, + max_order: usize, ) -> (Vector3, Matrix3) { let (v, w) = self.compute_legendre::(pos); let accel = self.accel_from_legendre_t::(&v, &w, max_order); @@ -448,12 +478,38 @@ impl Gravity { pos: &Vector3, max_order: usize, ) -> Vector3 { - let (v, w) = self.compute_legendre::(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::(pos, max_order) }; + } + } + self.accel_t_inner::(pos, max_order) + } + #[cfg(target_arch = "x86_64")] + #[target_feature(enable = "fma")] + unsafe fn accel_t_fma( + &self, + pos: &Vector3, + max_order: usize, + ) -> Vector3 { + self.accel_t_inner::(pos, max_order) + } + + #[inline(always)] + fn accel_t_inner( + &self, + pos: &Vector3, + max_order: usize, + ) -> Vector3 { + let (v, w) = self.compute_legendre::(pos); self.accel_from_legendre_t::(&v, &w, max_order) } // Equations 7.65 to 7.69 in Montenbruck & Gill + #[inline(always)] fn partials_from_legendre_t( &self, v: &Legendre, @@ -556,6 +612,7 @@ impl Gravity { } /// See Equation 3.33 in Montenbruck & Gill + #[inline(always)] fn accel_from_legendre_t( &self, v: &Legendre, @@ -601,6 +658,7 @@ impl Gravity { numeris::vector![ax, ay, az] * self.gravity_constant / self.radius / self.radius } + #[inline(always)] fn compute_legendre(&self, pos: &Vector3) -> (Legendre, Legendre) { let rsq = pos.norm_squared(); let scale = self.radius / rsq; From 487cc69749e27afb3b72f9ac0b3d672f1bf05c64 Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Thu, 10 Sep 2026 18:22:17 -0600 Subject: [PATCH 2/6] Compute globe7 and glob7s once per step --- src/nrlmsise.rs | 42 ++++++++++++++++++++++++++++-------------- 1 file changed, 28 insertions(+), 14 deletions(-) diff --git a/src/nrlmsise.rs b/src/nrlmsise.rs index dedd3f2d..28657511 100644 --- a/src/nrlmsise.rs +++ b/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 { From 1767c9dd0f4959b8514f77913d9f84ec5f126998 Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Thu, 10 Sep 2026 18:22:59 -0600 Subject: [PATCH 3/6] Switch to a single Bowring geodetic conversion step --- src/itrfcoord.rs | 28 +++++++++++----------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/src/itrfcoord.rs b/src/itrfcoord.rs index 8aacd66e..02c09aed 100644 --- a/src/itrfcoord.rs +++ b/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(); From 6c2c808ceceeab5f4436e646039adfb373a44b92 Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Thu, 10 Sep 2026 18:29:11 -0600 Subject: [PATCH 4/6] Reduce space weather reads --- src/nrlmsise.rs | 14 +++++++++----- src/spaceweather.rs | 8 ++++++-- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/nrlmsise.rs b/src/nrlmsise.rs index 28657511..7c840be8 100644 --- a/src/nrlmsise.rs +++ b/src/nrlmsise.rs @@ -5026,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) { diff --git a/src/spaceweather.rs b/src/spaceweather.rs index c04f5a01..f1b22ecc 100644 --- a/src/spaceweather.rs +++ b/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(tm: &T) -> Result { 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". From abe2f0ae0e2ab99b6b29e348f16a51a6eb29c5e2 Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Thu, 10 Sep 2026 19:19:22 -0600 Subject: [PATCH 5/6] Speed up import --- python/satkit/__init__.py | 5 +---- python/src/lib.rs | 1 + 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/python/satkit/__init__.py b/python/satkit/__init__.py index 8f222bea..916b23f1 100644 --- a/python/satkit/__init__.py +++ b/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. diff --git a/python/src/lib.rs b/python/src/lib.rs index d11ef29d..bd7ae476 100644 --- a/python/src/lib.rs +++ b/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::()?; m.add_class::()?; m.add_class::()?; From fdc5b2011475caaf527ef329942621fe0e40cc4e Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Thu, 10 Sep 2026 19:20:23 -0600 Subject: [PATCH 6/6] Skip parsing unneeded gravity data --- src/earthgravity.rs | 38 +++++++++++++++----------------------- 1 file changed, 15 insertions(+), 23 deletions(-) diff --git a/src/earthgravity.rs b/src/earthgravity.rs index 9946423f..0135556d 100644 --- a/src/earthgravity.rs +++ b/src/earthgravity.rs @@ -753,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), @@ -795,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()?; } }