Add Picard-Chebyshev numerical integrator

Cargo.toml

@@ -18,7 +18,7 @@ default-members = ["src/kete_core"]
 [workspace.package]
 version = "2.1.5"
 edition = "2024"
-rust-version = "1.85"
+rust-version = "1.90"
 license = "BSD-3-Clause"
 license-file = "LICENSE"
 repository = "https://github.com/dahlend/kete"

src/kete/rust/lib.rs

@@ -112,6 +112,7 @@ fn _core(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(propagation::propagation_n_body_spk_py, m)?)?;
     m.add_function(wrap_pyfunction!(propagation::propagation_n_body_py, m)?)?;
     m.add_function(wrap_pyfunction!(propagation::moid_py, m)?)?;
+    m.add_function(wrap_pyfunction!(propagation::picard, m)?)?;
 
     m.add_function(wrap_pyfunction!(fovs::fov_checks_py, m)?)?;
     m.add_function(wrap_pyfunction!(fovs::fov_spk_checks_py, m)?)?;

src/kete/rust/propagation.rs

@@ -284,3 +284,103 @@ pub fn propagation_n_body_py(
     }
     Ok((final_states, final_planets))
 }
+
+/// Debugging for picard
+#[pyfunction]
+#[pyo3(name = "picard", signature = (states, jd, include_asteroids=false,
+    non_gravs=None, suppress_errors=true, suppress_impact_errors=false))]
+pub fn picard(
+    py: Python<'_>,
+    states: MaybeVec<PyState>,
+    jd: PyTime,
+    include_asteroids: bool,
+    non_gravs: Option<MaybeVec<Option<PyNonGravModel>>>,
+    suppress_errors: bool,
+    suppress_impact_errors: bool,
+) -> PyResult<PyObject> {
+    let (states, was_vec): (Vec<_>, bool) = states.into();
+    let non_gravs: Option<Vec<Option<PyNonGravModel>>> = non_gravs.map(|x| x.into());
+    let non_gravs = non_gravs.unwrap_or(vec![None; states.len()]);
+
+    if states.len() != non_gravs.len() {
+        Err(Error::ValueError(
+            "non_gravs must be the same length as states.".into(),
+        ))?;
+    }
+
+    let mut res: Vec<PyState> = Vec::new();
+    let jd = jd.jd();
+
+    // propagation is broken into chunks, every time a chunk is completed
+    // python is checked for signals. This allows keyboard interrupts to be caught
+    // and the process interrupted.
+
+    for chunk in states
+        .into_iter()
+        .zip(non_gravs.into_iter())
+        .collect_vec()
+        .chunks(1000)
+    {
+        py.check_signals()?;
+
+        let mut proc_chunk = chunk
+            .to_owned()
+            .into_par_iter()
+            .with_min_len(5)
+            .map(|(state, model)| {
+                let model = model.map(|x| x.0);
+                let center = state.center_id();
+                let frame = state.frame;
+                let state = state.raw;
+                let desig = state.desig.clone();
+
+                // if the input has a NAN in it, skip the propagation entirely and return
+                // the nans.
+                if !state.is_finite() {
+                    if !suppress_errors {
+                        Err(Error::ValueError("Input state contains NaNs.".into()))?;
+                    };
+                    return Ok(Into::<PyState>::into(State::<Ecliptic>::new_nan(
+                        desig, jd, center,
+                    ))
+                    .change_frame(frame));
+                }
+                match propagation::propagate_picard_n_body_spk(state, jd, include_asteroids, model)
+                {
+                    Ok(state) => Ok(Into::<PyState>::into(state).change_frame(frame)),
+                    Err(er) => {
+                        if !suppress_errors {
+                            Err(er)?
+                        } else if let Error::Impact(id, time) = er {
+                            if !suppress_impact_errors {
+                                eprintln!(
+                                    "Impact detected between ({}) <-> {} at time {} ({})",
+                                    desig,
+                                    spice::try_name_from_id(id).unwrap_or(id.to_string()),
+                                    time,
+                                    Time::<TDB>::new(time).utc().to_iso().unwrap()
+                                );
+                            }
+                            // if we get an impact, we return a state with NaNs
+                            // but put the impact time into the new state.
+                            Ok(Into::<PyState>::into(State::<Ecliptic>::new_nan(
+                                desig, time, center,
+                            ))
+                            .change_frame(frame))
+                        } else {
+                            Ok(
+                                Into::<PyState>::into(State::<Ecliptic>::new_nan(
+                                    desig, jd, center,
+                                ))
+                                .change_frame(frame),
+                            )
+                        }
+                    }
+                }
+            })
+            .collect::<PyResult<Vec<_>>>()?;
+        res.append(&mut proc_chunk);
+    }
+
+    maybe_vec_to_pyobj(py, res, was_vec)
+}

src/kete/rust/time.rs

@@ -128,10 +128,10 @@ impl PyTime {
         // attempt to make life easier for the user by checking if they are missing
         // the timezone information. If they are, append it and return. Otherwise
         // let the conversion fail as it normally would.
-        if !s.contains('+') {
-            if let Ok(t) = Time::<UTC>::from_iso(&(s.to_owned() + "+00:00")) {
-                return Ok(PyTime(t.tdb()));
-            }
+        if !s.contains('+')
+            && let Ok(t) = Time::<UTC>::from_iso(&(s.to_owned() + "+00:00"))
+        {
+            return Ok(PyTime(t.tdb()));
         }
         Ok(PyTime(Time::<UTC>::from_iso(s)?.tdb()))
     }

src/kete_core/data/masses.tsv

@@ -10,17 +10,17 @@
 6        0.00028588567002459455 0.0004028667
 7        4.36624961322212e-05   0.0001708514
 8        5.15138377265457e-05   0.0001655371
-1000012  4.989735701e-18        2.27e-8
-20000001 4.7191659919436e-10    6.276827e-6
-20000002 1.0259143964958e-10    3.415824e-6
+1000012  4.989735701e-18
+20000001 4.7191659919436e-10
+20000002 1.0259143964958e-10
 20000003 1.447167e-11
-20000004 1.3026452498655e-10    3.512082e-6
+20000004 1.3026452498655e-10
 20000005 1.086874e-12
 20000006 4.874281e-12
 20000007 8.589039e-12
 20000008 1.939855e-12
 20000009 4.896769e-12
-20000010 4.3953391300849e-11    2.894426e-6
+20000010 4.3953391300849e-11
 20000011 3.466748e-12
 20000012 1.109906e-12
 20000013 3.848924e-12
@@ -316,7 +316,7 @@
 20000694 1.723815e-13
 20000696 2.289096e-13
 20000702 5.653793e-12
-20000704 1.911e-11            2.219283e-6
+20000704 1.911e-11
 20000705 1.321904e-12
 20000709 4.032293e-13
 20000712 4.258305e-13

src/kete_core/src/errors.rs

@@ -48,8 +48,8 @@ pub enum Error {
     /// Input or variable exceeded expected or allowed bounds.
     ValueError(String),
 
-    /// Querying an SPK file failed due to it missing the requisite data.
-    DAFLimits(String),
+    /// Attempting to query outside of data limits.
+    ExceedsLimits(String),
 
     /// Attempting to load or convert to/from an Frame of reference which is not known.
     UnknownFrame(i32),
@@ -66,7 +66,10 @@ impl error::Error for Error {}
 impl fmt::Display for Error {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         match self {
-            Self::Convergence(s) | Self::ValueError(s) | Self::DAFLimits(s) | Self::IOError(s) => {
+            Self::Convergence(s)
+            | Self::ValueError(s)
+            | Self::ExceedsLimits(s)
+            | Self::IOError(s) => {
                 write!(f, "{s}")
             }
             Self::UnknownFrame(_) => {
@@ -87,7 +90,7 @@ impl From<Error> for PyErr {
     fn from(err: Error) -> Self {
         match err {
             Error::IOError(s)
-            | Error::DAFLimits(s)
+            | Error::ExceedsLimits(s)
             | Error::ValueError(s)
             | Error::Convergence(s) => Self::new::<exceptions::PyValueError, _>(s),
 

src/kete_core/src/flux/reflected.rs

@@ -277,12 +277,12 @@ impl HGParams {
         if let Some(albedo) = vis_albedo {
             // h is defined and albedo is defined, meaning diameter may be calculated.
             let expected_diam = c_hg / albedo.sqrt() * (10_f64).powf(-0.2 * h_mag);
-            if let Some(diam) = diam {
-                if (expected_diam - diam).abs() > 1e-8 {
-                    Err(Error::ValueError(format!(
-                        "Provided diameter doesn't match with computed diameter. {expected_diam} != {diam}"
-                    )))?;
-                }
+            if let Some(diam) = diam
+                && (expected_diam - diam).abs() > 1e-8
+            {
+                Err(Error::ValueError(format!(
+                    "Provided diameter doesn't match with computed diameter. {expected_diam} != {diam}"
+                )))?;
             }
             return Ok((h_mag, Some(albedo), Some(expected_diam), c_hg));
         }

src/kete_core/src/flux/sun.rs

@@ -106,6 +106,7 @@ pub fn solar_flux_black_body(dist: f64, wavelength: f64) -> f64 {
 /// <https://www.nrel.gov/grid/solar-resource/spectra-astm-e490.html>
 ///
 /// First column is wavelength in um, second column is Watts / m^2 / micron at 1 AU.
+#[allow(clippy::approx_constant, reason = "These are not constants")]
 const E490_DATA: &[[f64; 2]] = &[
     [0.1195, 0.0619],
     [0.1205, 0.5614],

src/kete_core/src/fov/fov_like.rs

@@ -182,10 +182,10 @@ pub trait FovLike: Sync + Sized {
 
                 if (state.jd - obs_state.jd).abs() < dt_limit {
                     let (_, contains, _) = self.check_linear(state);
-                    if let Contains::Outside(dist) = contains {
-                        if dist > max_dist {
-                            return None;
-                        }
+                    if let Contains::Outside(dist) = contains
+                        && dist > max_dist
+                    {
+                        return None;
                     }
                     let (idx, contains, state) = self.check_two_body(state).ok()?;
                     match contains {
@@ -194,10 +194,10 @@ pub trait FovLike: Sync + Sized {
                     }
                 } else {
                     let (_, contains, _) = self.check_two_body(state).ok()?;
-                    if let Contains::Outside(dist) = contains {
-                        if dist > max_dist {
-                            return None;
-                        }
+                    if let Contains::Outside(dist) = contains
+                        && dist > max_dist
+                    {
+                        return None;
                     }
                     let (idx, contains, state) =
                         self.check_n_body(state, include_asteroids).ok()?;

src/kete_core/src/frames/definitions.rs

@@ -226,13 +226,13 @@ impl NonInertialFrame {
                     ));
                 }
             }
-            Err(Error::DAFLimits(format!(
+            Err(Error::ExceedsLimits(format!(
                 "Reference frame ID {} not found in CK data.",
                 self.reference_frame_id
             )))
         } else {
             // Unsupported frame
-            Err(Error::DAFLimits(format!(
+            Err(Error::ExceedsLimits(format!(
                 "Reference frame ID {} is not supported.",
                 self.reference_frame_id
             )))

src/kete_core/src/io/bytes.rs

@@ -57,7 +57,7 @@ pub(crate) fn bytes_to_f64(bytes: &[u8], little_endian: bool) -> KeteResult<f64>
 /// Change a collection of bytes into a vector of f64s.
 pub(crate) fn bytes_to_f64_vec(bytes: &[u8], little_endian: bool) -> KeteResult<Box<[f64]>> {
     let byte_len = bytes.len();
-    if byte_len % 8 != 0 {
+    if !byte_len.is_multiple_of(8) {
         Err(Error::IOError("File is not correctly formatted".into()))?;
     }
     let res: Box<[f64]> = (0..byte_len / 8)
@@ -69,7 +69,7 @@ pub(crate) fn bytes_to_f64_vec(bytes: &[u8], little_endian: bool) -> KeteResult<
 /// Change a collection of bytes into a vector of i32s.
 pub(crate) fn bytes_to_i32_vec(bytes: &[u8], little_endian: bool) -> KeteResult<Box<[i32]>> {
     let byte_len = bytes.len();
-    if byte_len % 4 != 0 {
+    if !byte_len.is_multiple_of(4) {
         Err(Error::IOError("File is not correctly formatted".into()))?;
     }
     let res: Box<[i32]> = (0..byte_len / 4)

src/kete_core/src/propagation/acceleration.rs

@@ -53,7 +53,7 @@ use crate::prelude::KeteResult;
 use crate::spice::LOADED_SPK;
 use crate::{constants, errors::Error, propagation::nongrav::NonGravModel};
 use nalgebra::allocator::Allocator;
-use nalgebra::{DefaultAllocator, Dim, Matrix3, OVector, U1, U2, Vector3};
+use nalgebra::{DefaultAllocator, Dim, Matrix3, OVector, SVector, U1, U2, Vector3};
 use std::ops::AddAssign;
 
 /// Metadata object used by the [`central_accel`] function below.
@@ -158,12 +158,11 @@ pub fn spk_accel(
 ) -> KeteResult<Vector3<f64>> {
     let mut accel = Vector3::<f64>::zeros();
 
-    if exact_eval {
-        if let Some(close_approach) = meta.close_approach.as_mut() {
-            if close_approach.2 == 0.0 {
-                *close_approach = (-1, 1000000.0, 0.0);
-            }
-        }
+    if exact_eval
+        && let Some(close_approach) = meta.close_approach.as_mut()
+        && close_approach.2 == 0.0
+    {
+        *close_approach = (-1, 1000000.0, 0.0);
     }
 
     let spk = &LOADED_SPK.try_read().unwrap();
@@ -177,10 +176,10 @@ pub fn spk_accel(
 
         if exact_eval {
             let r = rel_pos.norm();
-            if let Some(close_approach) = meta.close_approach.as_mut() {
-                if close_approach.2 >= r {
-                    *close_approach = (id, r, time);
-                }
+            if let Some(close_approach) = meta.close_approach.as_mut()
+                && close_approach.2 >= r
+            {
+                *close_approach = (id, r, time);
             }
 
             if r as f32 <= radius {
@@ -190,15 +189,35 @@ pub fn spk_accel(
         grav_params.add_acceleration(&mut accel, &rel_pos, &rel_vel);
 
         // If the center is the sun, add non-gravitational forces
-        if grav_params.naif_id == 10 {
-            if let Some(non_grav) = &meta.non_grav_model {
-                non_grav.add_acceleration(&mut accel, &rel_pos, &rel_vel);
-            }
+        if grav_params.naif_id == 10
+            && let Some(non_grav) = &meta.non_grav_model
+        {
+            non_grav.add_acceleration(&mut accel, &rel_pos, &rel_vel);
         }
     }
     Ok(accel)
 }
 
+/// Convert the second order ODE acceleration function into a first order.
+/// This allows the second order ODE to be used with the picard integrator.
+///
+/// The `state_vec` is made up of concatenated position and velocity vectors.
+/// Otherwise this is just a thin wrapper over the [`spk_accel`] function.
+pub fn spk_accel_first_order(
+    time: f64,
+    state_vec: &SVector<f64, 6>,
+    meta: &mut AccelSPKMeta<'_>,
+    exact_eval: bool,
+) -> KeteResult<SVector<f64, 6>> {
+    let pos: Vector3<f64> = state_vec.fixed_rows::<3>(0).into();
+    let vel: Vector3<f64> = state_vec.fixed_rows::<3>(3).into();
+    let accel = spk_accel(time, &pos, &vel, meta, exact_eval)?;
+    let mut res = SVector::<f64, 6>::zeros();
+    res.fixed_rows_mut::<3>(0).set_column(0, &vel);
+    res.fixed_rows_mut::<3>(3).set_column(0, &accel);
+    Ok(res)
+}
+
 /// Metadata for the [`vec_accel`] function defined below.
 #[derive(Debug)]
 pub struct AccelVecMeta<'a> {
@@ -269,11 +288,12 @@ where
             grav_params.add_acceleration(&mut accel_working, &rel_pos, &rel_vel);
 
             // If the center is the sun, add non-gravitational forces
-            if (grav_params.naif_id == 10) & (idx > n_massive) {
-                if let Some(non_grav) = &meta.non_gravs[idx - n_massive] {
-                    non_grav.add_acceleration(&mut accel_working, &rel_pos, &rel_vel);
-                }
+            if (grav_params.naif_id == 10) & (idx > n_massive)
+                && let Some(non_grav) = &meta.non_gravs[idx - n_massive]
+            {
+                non_grav.add_acceleration(&mut accel_working, &rel_pos, &rel_vel);
             }
+
             accel.fixed_rows_mut::<3>(idx * 3).add_assign(accel_working);
         }
     }
@@ -296,7 +316,7 @@ pub fn central_accel_grad(
 
 /// Calculate the Jacobian for the [`central_accel`] function.
 ///
-/// This enables the computation of the STM.
+/// This enables the computation of the two body STM.
 pub fn accel_grad(
     obj_pos: &Vector3<f64>,
     _obj_vel: &Vector3<f64>,