diff --git a/Cargo.toml b/Cargo.toml index eda5c98..56bdb73 100644 --- a/Cargo.toml +++ b/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" diff --git a/src/kete/rust/lib.rs b/src/kete/rust/lib.rs index b8729f4..da8a8ec 100644 --- a/src/kete/rust/lib.rs +++ b/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)?)?; diff --git a/src/kete/rust/propagation.rs b/src/kete/rust/propagation.rs index a09e29d..b9a26e5 100644 --- a/src/kete/rust/propagation.rs +++ b/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, + jd: PyTime, + include_asteroids: bool, + non_gravs: Option>>, + suppress_errors: bool, + suppress_impact_errors: bool, +) -> PyResult { + let (states, was_vec): (Vec<_>, bool) = states.into(); + let non_gravs: Option>> = 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 = 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::::into(State::::new_nan( + desig, jd, center, + )) + .change_frame(frame)); + } + match propagation::propagate_picard_n_body_spk(state, jd, include_asteroids, model) + { + Ok(state) => Ok(Into::::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::::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::::into(State::::new_nan( + desig, time, center, + )) + .change_frame(frame)) + } else { + Ok( + Into::::into(State::::new_nan( + desig, jd, center, + )) + .change_frame(frame), + ) + } + } + } + }) + .collect::>>()?; + res.append(&mut proc_chunk); + } + + maybe_vec_to_pyobj(py, res, was_vec) +} diff --git a/src/kete/rust/time.rs b/src/kete/rust/time.rs index 12d6869..5677033 100644 --- a/src/kete/rust/time.rs +++ b/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::::from_iso(&(s.to_owned() + "+00:00")) { - return Ok(PyTime(t.tdb())); - } + if !s.contains('+') + && let Ok(t) = Time::::from_iso(&(s.to_owned() + "+00:00")) + { + return Ok(PyTime(t.tdb())); } Ok(PyTime(Time::::from_iso(s)?.tdb())) } diff --git a/src/kete_core/data/masses.tsv b/src/kete_core/data/masses.tsv index 2166262..561093a 100644 --- a/src/kete_core/data/masses.tsv +++ b/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 diff --git a/src/kete_core/src/errors.rs b/src/kete_core/src/errors.rs index e697ea9..41d4cb3 100644 --- a/src/kete_core/src/errors.rs +++ b/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 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::(s), diff --git a/src/kete_core/src/flux/reflected.rs b/src/kete_core/src/flux/reflected.rs index af7d5bd..435869c 100644 --- a/src/kete_core/src/flux/reflected.rs +++ b/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)); } diff --git a/src/kete_core/src/flux/sun.rs b/src/kete_core/src/flux/sun.rs index f906735..13900ef 100644 --- a/src/kete_core/src/flux/sun.rs +++ b/src/kete_core/src/flux/sun.rs @@ -106,6 +106,7 @@ pub fn solar_flux_black_body(dist: f64, wavelength: f64) -> f64 { /// /// /// 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], diff --git a/src/kete_core/src/fov/fov_like.rs b/src/kete_core/src/fov/fov_like.rs index f286787..687c764 100644 --- a/src/kete_core/src/fov/fov_like.rs +++ b/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()?; diff --git a/src/kete_core/src/frames/definitions.rs b/src/kete_core/src/frames/definitions.rs index a179634..55ec0a1 100644 --- a/src/kete_core/src/frames/definitions.rs +++ b/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 ))) diff --git a/src/kete_core/src/io/bytes.rs b/src/kete_core/src/io/bytes.rs index 4bb7bda..7e3978a 100644 --- a/src/kete_core/src/io/bytes.rs +++ b/src/kete_core/src/io/bytes.rs @@ -57,7 +57,7 @@ pub(crate) fn bytes_to_f64(bytes: &[u8], little_endian: bool) -> KeteResult /// Change a collection of bytes into a vector of f64s. pub(crate) fn bytes_to_f64_vec(bytes: &[u8], little_endian: bool) -> KeteResult> { 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> { 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) diff --git a/src/kete_core/src/propagation/acceleration.rs b/src/kete_core/src/propagation/acceleration.rs index f14d2e7..e8d1176 100644 --- a/src/kete_core/src/propagation/acceleration.rs +++ b/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> { let mut accel = Vector3::::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, + meta: &mut AccelSPKMeta<'_>, + exact_eval: bool, +) -> KeteResult> { + let pos: Vector3 = state_vec.fixed_rows::<3>(0).into(); + let vel: Vector3 = state_vec.fixed_rows::<3>(3).into(); + let accel = spk_accel(time, &pos, &vel, meta, exact_eval)?; + let mut res = SVector::::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, _obj_vel: &Vector3, diff --git a/src/kete_core/src/propagation/kepler.rs b/src/kete_core/src/propagation/kepler.rs index d0e1b27..accb714 100644 --- a/src/kete_core/src/propagation/kepler.rs +++ b/src/kete_core/src/propagation/kepler.rs @@ -82,7 +82,7 @@ pub fn compute_eccentric_anomaly(ecc: f64, mean_anom: f64, peri_dist: f64) -> Ke } ecc => { // Hyperbolic - let f = |ecc_anom: f64| (ecc * ecc_anom.sinh() - ecc_anom - mean_anom); + let f = |ecc_anom: f64| ecc * ecc_anom.sinh() - ecc_anom - mean_anom; let d = |ecc_anom: f64| -1.0 + ecc * ecc_anom.cosh(); let guess = (mean_anom / ecc).asinh(); Ok(newton_raphson(f, d, guess, 1e-11)?) diff --git a/src/kete_core/src/propagation/mod.rs b/src/kete_core/src/propagation/mod.rs index 4c5e637..ede93fe 100644 --- a/src/kete_core/src/propagation/mod.rs +++ b/src/kete_core/src/propagation/mod.rs @@ -38,25 +38,28 @@ use crate::frames::Equatorial; use crate::prelude::{Desig, KeteResult}; use crate::spice::LOADED_SPK; use crate::state::State; -use nalgebra::{DVector, Vector3}; +use nalgebra::{DVector, SMatrix, SVector, Vector3}; mod acceleration; mod kepler; mod nongrav; +mod picard; mod radau; mod runge_kutta; mod state_transition; +mod util; // expose the public methods in spk to the outside world. pub use acceleration::{ AccelSPKMeta, AccelVecMeta, CentralAccelMeta, accel_grad, central_accel, central_accel_grad, - spk_accel, vec_accel, + spk_accel, spk_accel_first_order, vec_accel, }; pub use kepler::{ PARABOLIC_ECC_LIMIT, analytic_2_body, compute_eccentric_anomaly, compute_true_anomaly, eccentric_anomaly_from_true, moid, propagate_two_body, }; pub use nongrav::NonGravModel; +pub use picard::{PC15, PC25, PicardIntegrator, PicardStep, dumb_picard_init}; pub use radau::RadauIntegrator; pub use runge_kutta::RK45Integrator; pub use state_transition::compute_state_transition; @@ -141,6 +144,85 @@ pub fn propagate_n_body_spk( Ok(new_state) } +/// Initialization function for the picard integrator which initializes the state +/// using two body mechanics. +fn picard_two_body_init( + times: &[f64; N], + init_pos: &SVector, +) -> SMatrix { + let pos: Vector3 = init_pos.fixed_rows::<3>(0).into(); + let vel: Vector3 = init_pos.fixed_rows::<3>(3).into(); + let t0 = times[0]; + + let mut res: SMatrix = SMatrix::zeros(); + res.fixed_rows_mut::<3>(0).set_column(0, &pos); + res.fixed_rows_mut::<3>(3).set_column(0, &vel); + + for (idx, t) in times.iter().enumerate().skip(1) { + let dt = t - t0; + let (p, v) = analytic_2_body(dt, &pos, &vel, None).unwrap(); + + res.fixed_rows_mut::<3>(0).set_column(idx, &p); + res.fixed_rows_mut::<3>(3).set_column(idx, &v); + } + res +} + +/// Propagate an object using full N-Body physics with the Radau 15th order integrator. +pub fn propagate_picard_n_body_spk( + mut state: State, + jd_final: f64, + include_extended: bool, + non_grav_model: Option, +) -> KeteResult> { + let center = state.center_id; + let spk = &LOADED_SPK.try_read().unwrap(); + spk.try_change_center(&mut state, 0)?; + + let mass_list = { + if include_extended { + &GravParams::selected_masses() + } else { + &GravParams::planets() + } + }; + + let mut metadata = AccelSPKMeta { + close_approach: None, + non_grav_model, + massive_obj: mass_list, + }; + + let integrator = &PC15; + + let mut state_vec = SVector::::zeros(); + state_vec + .fixed_rows_mut::<3>(0) + .set_column(0, &state.pos.into()); + state_vec + .fixed_rows_mut::<3>(3) + .set_column(0, &state.vel.into()); + + let final_state_vec = { + integrator.integrate( + &spk_accel_first_order, + &picard_two_body_init, + state_vec, + state.jd, + jd_final, + 1.0, + &mut metadata, + )? + }; + + let pos: Vector3 = final_state_vec.fixed_rows::<3>(0).into(); + let vel: Vector3 = final_state_vec.fixed_rows::<3>(3).into(); + + let mut new_state = State::new(state.desig.to_owned(), jd_final, pos.into(), vel.into(), 0); + spk.try_change_center(&mut new_state, center)?; + Ok(new_state) +} + /// Propagate an object using two body mechanics. /// This is a brute force way to solve the kepler equations of motion as it uses Radau /// as an integrator. diff --git a/src/kete_core/src/propagation/picard.rs b/src/kete_core/src/propagation/picard.rs new file mode 100644 index 0000000..46f1e88 --- /dev/null +++ b/src/kete_core/src/propagation/picard.rs @@ -0,0 +1,481 @@ +/// # Picard-Chebyshev Numerical Integrator +/// +/// First Order ODE Integrator. +/// +/// This is based on a whole slew of papers, though this implementation leans heavily +/// on the algorithm described in: +/// +/// "Surfing Chaotic Perturbations in Interplanetary Multi-Flyby Trajectories: +/// Augmented Picard-Chebyshev Integration for Parallel and GPU Computing +/// Architectures", 2022, `` +/// +/// Though that paper appears to have a typo in its `A` matrix definition, so some of +/// the thesis written by Darin Koblick "Parallel High-Precision Orbit Propagation +/// Using The Modified Picard-Chebyshev Method", 2012 +/// was used to correct the matrix definition. +/// +/// This is not a GPU implementation, though some small simplifications were applied to +/// the mathematics. This implementation is designed to reduce the total number of +/// SPICE kernel calls to a minimum, as it was found that the Radau integration time +/// was dominated by these queries. +/// +/// Since this integrator fits Chebyshev polynomials at the same time that it performs +/// the integration, the integrator is designed to record the polynomial coefficients. +/// These are stored in the [`PicardStep`], which exposes a function allowing the +/// user to query the state of the system at any point between the start of the +/// integration and the end. +/// +// BSD 3-Clause License +// +// Copyright (c) 2025, Dar Dahlen +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// 1. Redistributions of source code must retain the above copyright notice, this +// list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// 3. Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +use std::f64::consts::PI; + +use itertools::Itertools; +use nalgebra::{SMatrix, SVector}; + +use crate::errors::{Error, KeteResult}; +use crate::propagation::util::FirstOrderODE; + +/// Initialization function which fills the integrators initial guess. +type PicardInitFunc<'a, const DIM: usize, const N: usize> = + &'a dyn Fn(&[f64; N], &SVector) -> SMatrix; + +/// Initialization function for the picard integrator where the first value is repeated +pub fn dumb_picard_init( + _times: &[f64; N], + initial_pos: &SVector, +) -> SMatrix { + let mut init: SMatrix = SMatrix::zeros(); + for idx in 0..N { + init.set_column(idx, initial_pos); + } + init +} + +// ______ +// ___.--------'------`---------.____ +// _.---'----------------------------------`---.__ +// .'___=]=========================================== +// ,-----------------------..__/.' >--.______ _______.---' +// ]====================<==||(__) .' `------' +// `-----------------------`' ----.___--/ +// / /---' `/ +// /_______(______________________/ Boldly going where many numerical +// `-------------.--------------.' integrators have gone before. +// \________|_.-' +// +/// Picard-Chebyshev integrator for solving first order ODEs. +/// +/// This integrator fits Chebyshev polynomials to the provided first order ODEs. +/// These polynomials can be retrieved and queried at any point between the start and +/// end point of the integration. This is designed so that an object may have its orbit +/// integrated any length of time, and there is a very fast lookup for the object's +/// state at any point along the integration. +/// +/// This integrator has a set of matrices which are pre-computed and stored. +/// It is recommended to use one of the pre-build integrators [`PC15`] or [`PC25`] to +/// avoid having to reconstruct these matrices every time the integrator is used. +/// +/// The mathematics for this integrator can be found in a number of papers, but +/// `` is a good place to start. +/// +/// The primary entry-point for the integrator is the [`PicardIntegrator::integrate`] +/// function. See that function for more details on its use. +/// +#[derive(Debug, Clone)] +pub struct PicardIntegrator { + // These parameters names are taken directly from the paper cited above. + // SA represents the matrix product of S and A. + c: SMatrix, + + a: SMatrix, + + sa: SVector, + + tau: [f64; N], +} + +/// Preconstructed [`PicardIntegrator`] for a 14th order polynomial +pub static PC15: std::sync::LazyLock> = + std::sync::LazyLock::new(PicardIntegrator::default); + +/// Preconstructed [`PicardIntegrator`] for a 24th order polynomial +pub static PC25: std::sync::LazyLock> = + std::sync::LazyLock::new(PicardIntegrator::default); + +/// A single stop of the Picard Integrator +#[derive(Debug, Clone)] +pub struct PicardStep { + /// Chebyshev polynomial coefficients of the first type + b: SMatrix, + + /// Final result + pub y: SMatrix, + + /// Start time + pub t0: f64, + + /// End time + pub t1: f64, +} + +impl<'a, const N: usize, const NM1: usize> PicardIntegrator { + /// Integrate from the start point to the end. + /// + /// # Arguments + /// + /// * `func` - First order ODE function which accepts time, state of the system, a + /// metadata class, and a bool value to indicate if the provided state is + /// believed to be an exact evaluated position. + /// * `init_func` - A function which can initialize the starting guess for the + /// integrator. If this is unavailable, then [`dumb_picard_init`] can be used as + /// a basic initializer. + /// * `t0` - Start time of the integration. + /// * `t1` - End time of the integration. + /// * `step_size` - Initial step size of the integration, a rough guess is all that + /// is required, as the integrator will find its own valid step size. This speeds + /// up the initial search for the valid step size. + /// * `metadata` - An initial packet of metadata which will be passed through to + /// the `func` during evaluation, this may include information such as masses, + /// non-gravitational forces, record keeping of close encounters, etc. This + /// metadata is mutated in place during integration. + /// + /// # Returns + /// + /// Returns the expected state of the system at the final time step `t1`. + /// + pub fn integrate( + &self, + func: FirstOrderODE<'a, MType, DIM>, + init_func: PicardInitFunc<'a, DIM, N>, + initial_pos: SVector, + t0: f64, + t1: f64, + step_size: f64, + metadata: &mut MType, + ) -> KeteResult> { + let mut cur_t0 = t0; + + let mut cur_stepsize = step_size; + let mut cur_t1 = t0 + cur_stepsize; + if (t0 - t1).abs() < step_size.abs() { + cur_stepsize = t1 - t0; + cur_t1 = t1; + } + + let mut times = self.evaluation_times(cur_t0, cur_t1); + let mut initial_pos = init_func(×, &initial_pos); + + let mut has_failed = false; + loop { + let step = self.step(func, cur_t0, cur_t1, initial_pos, metadata); + match step { + Ok(res) => { + if cur_t1 == t1 { + return Ok(res.y.column(N - 1).into()); + } + times = self.evaluation_times(cur_t0, cur_t1); + initial_pos = init_func(×, &res.y.column(N - 1).into()); + + cur_t0 = cur_t1; + cur_t1 += cur_stepsize; + if (cur_stepsize.is_sign_positive() && cur_t1 > t1) + || (cur_stepsize.is_sign_negative() && cur_t1 < t1) + { + cur_t1 = t1; + } + if !has_failed { + cur_stepsize *= 1.5; + } + } + Err(Error::Impact(idx, jd)) => { + return Err(Error::Impact(idx, jd)); + } + Err(_) => { + has_failed = true; + cur_stepsize *= 0.7; + cur_t1 = cur_t0 + cur_stepsize; + } + } + if cur_stepsize.abs() < 1e-10 { + return Err(Error::Convergence(format!( + "Failed to converge, step sizes became too small. {}", + cur_stepsize.abs() + ))); + } + } + } + + /// Single fallible step of the integrator + /// + /// This returns a single [`PicardStep`] which may be used to evaluate the state + /// of the system at any point between the start time `t0` and end `t1`. + /// + fn step( + &self, + func: FirstOrderODE<'a, MType, DIM>, + t0: f64, + t1: f64, + mut cur_pos: SMatrix, + metadata: &mut MType, + ) -> KeteResult> { + let mut b: SMatrix = SMatrix::zeros(); + let mut f: SMatrix = SMatrix::zeros(); + let times = self.evaluation_times(t0, t1); + let w2 = (t1 - t0) / 2.0; + let mut error: f64 = 100.0; + let mut last_error = 50.0; + + // when the answer converges to within tolerance, do one last evaluation and + // trigger the metadata exact evaluation, allowing the `func` to update the + // metadata with exact values. + let mut final_iteration = false; + + for _ in 0..150 { + for ((&time, y), mut f_row) in times + .iter() + .zip(cur_pos.column_iter()) + .zip(f.column_iter_mut()) + { + // during the final iteration, let `func` know that we should update + // the metadata + f_row.set_column(0, &(w2 * func(time, &y.into(), metadata, final_iteration)?)); + } + + b.set_column(0, &(f * self.sa + 2.0 * cur_pos.column(0))); + + b.fixed_view_mut::(0, 1) + .iter_mut() + .zip((f * self.a).into_iter()) + .for_each(|(x, &n)| *x = n); + + let last_pos = cur_pos; + cur_pos = b * self.c; + + // diverging, quit quickly + if error > 1000.0 && last_error > error { + return Err(Error::Convergence( + "Integrator solution is diverging.".into(), + )); + } + + // if we are on the final evaluation, return the answer. + if final_iteration { + return Ok(PicardStep { + b, + y: cur_pos, + t0, + t1, + }); + } + + // if we have converged, then trigger the final iteration + if error < 2.0 * f64::EPSILON { + final_iteration = true; + } + + last_error = error; + error = (last_pos - cur_pos).abs().max(); + } + + Err(Error::Convergence("Integrator failed to converge.".into())) + } + + /// Given the start and stop times of a desired integration, return the actual + /// times where the function will be evaluated starting at t0 and ending at t1. + fn evaluation_times(&self, t0: f64, t1: f64) -> [f64; N] { + let w2 = (t1 - t0) / 2.0; + let w1 = (t1 + t0) / 2.0; + let mut times = self.tau; + times.iter_mut().for_each(|x| *x = w2 * (*x) + w1); + times + } +} + +impl PicardStep { + /// Evaluate the fitted integration solution at the specified time. + /// This will fail if the requested time is outside of the integration bounds. + pub fn evaluate(&self, t: f64) -> KeteResult<[f64; DIM]> { + let w1 = (self.t0 + self.t1) * 0.5; + let w2 = (self.t1 - self.t0) * 0.5; + let tau_time = ((t - w1) * w2).acos(); + if tau_time.is_nan() { + return Err(Error::ExceedsLimits( + "Queried time it outside of the fitted time span".into(), + )); + } + Ok(chebyshev_eval(tau_time, &self.b.transpose().into())) + } +} + +impl Default for PicardIntegrator { + fn default() -> Self { + // NM1 must be equal to N-1, I wish I could do this with const generics + // but that would require switching to nightly rust and I wont do that. + // So this compile time check is done instead. + const { assert!(N - 1 == NM1, "NM1 must be 1 less than N.") }; + + // Combining the cos(k*arccos(-cos(...))) simplifies to this equation + // This is used throughout the integrator's construction + fn cheb_t(k: f64, j: f64, n: f64) -> f64 { + (k * (j + n - 1.0) * PI / (n - 1.0)).cos() + } + + let n = N as f64; + + // Construct the C matrix + let mut c = SMatrix::::zeros(); + // note that calling cheb_t with k=0 evaluates to 1, and in the paper this + // gets divided by 2, resulting in 0.5. + c.fill_column(0, 0.5); + + for j in 0..N { + for k in 1..N { + c[(j, k)] = cheb_t(k as f64, j as f64, n); + } + } + + // Construct the A matrix from `` + // Note that the paper appears to have typos in the indexing, I suspect + // there should be a k in several places where there is currently a j. + // This implementation is a blend of that paper, and the original math + // that the paper cites. + let mut a = SMatrix::::zeros(); + for k in 0..N - 2 { + for j in 0..N { + a[(k, j)] = (cheb_t(k as f64, j as f64, n) - cheb_t((k + 2) as f64, j as f64, n)) + / (((N - 1) * (k + 1)) as f64); + } + } + + for j in 0..N { + a[(N - 2, j)] = cheb_t(n - 2.0, j as f64, n) / ((N - 1).pow(2) as f64); + } + + a[(N - 2, 0)] /= 2.0; + a[(N - 2, N - 1)] /= 2.0; + + // Now we construct the SA Vector, in the original paper they keep S as a + // separate vector and multiply it against A over and over, whereas we are just + // doing it once and saving it. + let s: SVector = + SVector::from_iterator((0..N).map(|idx| 2.0 * (-1_f64).powi(idx as i32))); + + let sa: SVector = a.transpose() * s; + + let tau: [f64; N] = (0..N) + .map(|j| -(j as f64 * PI / (n - 1.0)).cos()) + .collect_vec() + .try_into() + .unwrap(); + + Self { + a: a.transpose(), + sa, + c: c.transpose(), + tau, + } + } +} + +/// Evaluate chebyshev polynomials of the first type. +#[inline(always)] +fn chebyshev_eval(x: f64, coef: &[[f64; R]; C]) -> [f64; C] { + const { assert!(R > 2, "Dimension R must be greater than 2.") }; + const { assert!(C > 0, "Dimension C must be greater than 0.") }; + + let mut val = [0.0; C]; + + for idx in 0..C { + val[idx] = coef[idx][0] + coef[idx][1] * x; + } + + let mut second_t = 1.0; + let mut last_t = x; + let mut next_t; + + let x2: f64 = 2.0 * x; + for idy in 2..R { + next_t = x2 * last_t - second_t; + for idx in 0..C { + val[idx] += coef[idx][idy] * next_t; + } + second_t = last_t; + last_t = next_t; + } + + val +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn exponential_decay_test() { + // integrate a first order ode of the form f' = -c f'(0) = k + // where c = [1, 0.5, 0.1] and k = [1, 2, 5] + // This has the analytic solution of f=k * exp(-c*t) + + // define the f' function + fn func( + _: f64, + vals: &SVector, + _: &mut (), + _: bool, + ) -> KeteResult> { + let v: SVector = [1.0, 0.5, 0.1].into(); + Ok(-vals.component_mul(&v)) + } + + // integrate with the initial conditions + let p = &PC15; + let t1 = 5.0; + let res = p + .integrate( + &func, + &dumb_picard_init, + [1.0, 2.0, 5.0].into(), + 0.0, + t1, + 0.1, + &mut (), + ) + .unwrap(); + + // test against the analytic solution + assert!((res[0] - (1.0 * (-t1).exp())).abs() < 1e-15); + assert!((res[1] - (2.0 * (-0.5 * t1).exp())).abs() < 1e-15); + assert!( + (res[2] - (5.0 * (-0.1 * t1).exp())).abs() < 1e-14, + "{}", + (res[2] - (5.0 * (-0.1 * t1).exp())).abs() + ); + } +} diff --git a/src/kete_core/src/propagation/radau.rs b/src/kete_core/src/propagation/radau.rs index ac1a150..2d40df1 100644 --- a/src/kete_core/src/propagation/radau.rs +++ b/src/kete_core/src/propagation/radau.rs @@ -31,21 +31,12 @@ // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. use crate::errors::Error; use crate::prelude::KeteResult; +use crate::propagation::util::SecondOrderODE; use itertools::izip; use nalgebra::allocator::Allocator; use nalgebra::*; use nalgebra::{DefaultAllocator, Dim, OMatrix, OVector, RowSVector, SMatrix, SVector, U1, U7}; -/// Function will be of the form y'' = F(t, y, y', metadata) -/// This is the second-order general solver (class II in the Everhart paper). -type RadauFunc<'a, MType, D> = &'a dyn Fn( - f64, - &OVector, - &OVector, - &mut MType, - bool, -) -> KeteResult>; - /// Integrator will return a result of this type. type RadauResult = KeteResult<(OVector, OVector, MType)>; @@ -118,7 +109,7 @@ pub struct RadauIntegrator<'a, MType, D: Dim> where DefaultAllocator: Allocator + Allocator, { - func: RadauFunc<'a, MType, D>, + func: SecondOrderODE<'a, MType, D>, metadata: MType, final_time: f64, @@ -142,7 +133,7 @@ where DefaultAllocator: Allocator + Allocator, { fn new( - func: RadauFunc<'a, MType, D>, + func: SecondOrderODE<'a, MType, D>, state_init: OVector, state_der_init: OVector, time_init: f64, @@ -183,7 +174,7 @@ where /// Integrate the functions from the initial time to the final time. pub fn integrate( - func: RadauFunc<'a, MType, D>, + func: SecondOrderODE<'a, MType, D>, state_init: OVector, state_der_init: OVector, time_init: f64, @@ -226,7 +217,7 @@ where } Err(error) => match error { Error::Impact(_, _) => Err(error)?, - Error::DAFLimits(_) => Err(error)?, + Error::ExceedsLimits(_) => Err(error)?, _ => { step_failures += 1; next_step_size *= 0.7; diff --git a/src/kete_core/src/propagation/state_transition.rs b/src/kete_core/src/propagation/state_transition.rs index f7156fe..f2b59eb 100644 --- a/src/kete_core/src/propagation/state_transition.rs +++ b/src/kete_core/src/propagation/state_transition.rs @@ -30,7 +30,9 @@ use crate::constants::GMS_SQRT; use crate::frames::Equatorial; use crate::prelude::{KeteResult, State}; -use crate::propagation::{CentralAccelMeta, RK45Integrator, central_accel, central_accel_grad}; +use crate::propagation::{ + CentralAccelMeta, PC15, central_accel, central_accel_grad, dumb_picard_init, +}; use nalgebra::{Const, Matrix6, SVector, U1, U6, Vector3}; fn stm_ivp_eqn( @@ -78,13 +80,13 @@ fn stm_ivp_eqn( /// Compute a state transition matrix assuming only 2-body mechanics. /// -/// This uses a Runge-Kutta 4/5 algorithm. +/// This uses the Picard-Chebyshev integrator [`PC15`]. pub fn compute_state_transition( state: &mut State, jd: f64, central_mass: f64, ) -> ([[f64; 3]; 2], Matrix6) { - let meta = CentralAccelMeta { + let mut meta = CentralAccelMeta { mass_scaling: central_mass, ..Default::default() }; @@ -102,18 +104,21 @@ pub fn compute_state_transition( initial_state .fixed_rows_mut::<3>(3) .set_column(0, &(Vector3::from(state.vel) / GMS_SQRT)); - let rad = RK45Integrator::integrate( - &stm_ivp_eqn, - initial_state, - state.jd * GMS_SQRT, - jd * GMS_SQRT, - meta, - 1e-12, - ) - .unwrap(); + + let integrator = &PC15; + let rad = integrator + .integrate( + &stm_ivp_eqn, + &dumb_picard_init, + initial_state, + state.jd * GMS_SQRT, + jd * GMS_SQRT, + 1.0, + &mut meta, + ) + .unwrap(); let vec_reshape = rad - .0 .fixed_rows::<36>(6) .into_owned() .reshape_generic(U6, U6) @@ -134,8 +139,8 @@ pub fn compute_state_transition( Matrix6::::from_diagonal(&[1.0, 1.0, 1.0, GMS_SQRT, GMS_SQRT, GMS_SQRT].into()); ( [ - rad.0.fixed_rows::<3>(0).into(), - (rad.0.fixed_rows::<3>(3) * GMS_SQRT).into(), + rad.fixed_rows::<3>(0).into(), + (rad.fixed_rows::<3>(3) * GMS_SQRT).into(), ], scaling_a * vec_reshape * scaling_b, ) diff --git a/src/kete_core/src/propagation/util.rs b/src/kete_core/src/propagation/util.rs new file mode 100644 index 0000000..a70189d --- /dev/null +++ b/src/kete_core/src/propagation/util.rs @@ -0,0 +1,18 @@ +use nalgebra::{OVector, SVector}; + +use crate::errors::KeteResult; + +/// Function will be of the form y' = F(time, y, metadata, bool) +/// This is the first-order general IVP solver. +pub(crate) type FirstOrderODE<'a, MType, const DIM: usize> = + &'a dyn Fn(f64, &SVector, &mut MType, bool) -> KeteResult>; + +/// Function will be of the form y'' = F(time, y, y', metadata, bool) +/// This is the input for a second-order general IVP solver. +pub(crate) type SecondOrderODE<'a, MType, D> = &'a dyn Fn( + f64, + &OVector, + &OVector, + &mut MType, + bool, +) -> KeteResult>; diff --git a/src/kete_core/src/spice/ck.rs b/src/kete_core/src/spice/ck.rs index ced371c..aca7d35 100644 --- a/src/kete_core/src/spice/ck.rs +++ b/src/kete_core/src/spice/ck.rs @@ -81,7 +81,8 @@ impl CkCollection { /// Get the closest record to the given JD for the specified instrument ID. /// /// # Errors - /// [`Error::DAFLimits`] if the instrument ID does not have a record for the target JD. + /// [`Error::ExceedsLimits`] if the instrument ID does not have a record for the + /// target JD. /// pub fn try_get_frame( &self, @@ -100,7 +101,7 @@ impl CkCollection { } } - Err(Error::DAFLimits(format!( + Err(Error::ExceedsLimits(format!( "Instrument ({instrument_id}) does not have an CK record for the target JD." )))? } diff --git a/src/kete_core/src/spice/ck_segments.rs b/src/kete_core/src/spice/ck_segments.rs index b3d6021..6881125 100644 --- a/src/kete_core/src/spice/ck_segments.rs +++ b/src/kete_core/src/spice/ck_segments.rs @@ -49,7 +49,7 @@ impl CkSegment { ) -> KeteResult<(Time, NonInertialFrame)> { let arr_ref: &CkArray = self.into(); if arr_ref.instrument_id != instrument_id { - return Err(Error::DAFLimits(format!( + return Err(Error::ExceedsLimits(format!( "Instrument ID is not present in this record. {}", arr_ref.instrument_id ))); @@ -138,7 +138,7 @@ impl CkSegmentType2 { ) -> KeteResult<(Time, NonInertialFrame)> { let sclk = LOADED_SCLK .try_read() - .map_err(|_| Error::DAFLimits("Failed to read SCLK data.".into()))?; + .map_err(|_| Error::ExceedsLimits("Failed to read SCLK data.".into()))?; let tick = sclk.try_time_to_tick(self.array.naif_id, time)?; // get the time of the last record and its index @@ -165,7 +165,7 @@ impl CkSegmentType2 { let dt = tick - record_time; if dt < 0.0 { - return Err(Error::DAFLimits(format!( + return Err(Error::ExceedsLimits(format!( "Requested time {record_idx} is before the start of the segment." ))); } @@ -204,13 +204,13 @@ impl TryFrom for CkSegmentType2 { dir_size = n_records / 100; if array_len != (n_records * 10 + dir_size) { - return Err(Error::DAFLimits( + return Err(Error::ExceedsLimits( "CK File is not formatted correctly, directory size of segments appear incorrect." .into(), )); } if n_records == 0 { - return Err(Error::DAFLimits( + return Err(Error::ExceedsLimits( "CK File does not contain any records.".into(), )); } @@ -348,7 +348,7 @@ impl CkSegmentType3 { ) -> KeteResult<(Time, UnitQuaternion, Option<[f64; 3]>)> { let sclk = LOADED_SCLK .try_read() - .map_err(|_| Error::DAFLimits("Failed to read SCLK data.".into()))?; + .map_err(|_| Error::ExceedsLimits("Failed to read SCLK data.".into()))?; let tick = sclk.try_time_to_tick(self.array.naif_id, time)?; // If there is only one record, return it immediately. @@ -427,12 +427,12 @@ impl TryFrom for CkSegmentType3 { let n_intervals = array.daf[array.daf.len() - 2] as usize; if n_records == 0 { - return Err(Error::DAFLimits( + return Err(Error::ExceedsLimits( "CK File does not contain any records.".into(), )); } if n_intervals == 0 { - return Err(Error::DAFLimits( + return Err(Error::ExceedsLimits( "CK File does not contain any intervals of records.".into(), )); } @@ -455,7 +455,7 @@ impl TryFrom for CkSegmentType3 { expected_size += time_dir_size + interval_dir_size; if expected_size != array.daf.len() { - return Err(Error::DAFLimits( + return Err(Error::ExceedsLimits( "CK File not formatted correctly. Number of records found in file don't match expected." .into(), )); diff --git a/src/kete_core/src/spice/pck.rs b/src/kete_core/src/spice/pck.rs index eab92a4..9b93465 100644 --- a/src/kete_core/src/spice/pck.rs +++ b/src/kete_core/src/spice/pck.rs @@ -82,7 +82,7 @@ impl PckCollection { } } - Err(Error::DAFLimits(format!( + Err(Error::ExceedsLimits(format!( "Object ({id}) does not have an PCK record for the target JD." )))? } @@ -124,10 +124,10 @@ impl PckCollection { let path = entry.path(); if path.is_file() { let filename = path.to_str().unwrap(); - if filename.to_lowercase().ends_with(".bpc") { - if let Err(err) = self.load_file(filename) { - eprintln!("Failed to load PCK file {filename}: {err}"); - } + if filename.to_lowercase().ends_with(".bpc") + && let Err(err) = self.load_file(filename) + { + eprintln!("Failed to load PCK file {filename}: {err}"); } } }); diff --git a/src/kete_core/src/spice/pck_segments.rs b/src/kete_core/src/spice/pck_segments.rs index ffd405c..daff0bc 100644 --- a/src/kete_core/src/spice/pck_segments.rs +++ b/src/kete_core/src/spice/pck_segments.rs @@ -95,7 +95,7 @@ impl PckSegment { let arr_ref: &PckArray = self.into(); if center_id != arr_ref.frame_id { - Err(Error::DAFLimits( + Err(Error::ExceedsLimits( "Center ID is not present in this record.".into(), ))?; } @@ -103,7 +103,9 @@ impl PckSegment { let jds = jd_to_spice_jd(jd); if jds < arr_ref.jds_start || jds > arr_ref.jds_end { - Err(Error::DAFLimits("JD is not present in this record.".into()))?; + Err(Error::ExceedsLimits( + "JD is not present in this record.".into(), + ))?; } if arr_ref.reference_frame_id != 17 { Err(Error::ValueError(format!( diff --git a/src/kete_core/src/spice/sclk.rs b/src/kete_core/src/spice/sclk.rs index b94cbe3..10c4e07 100644 --- a/src/kete_core/src/spice/sclk.rs +++ b/src/kete_core/src/spice/sclk.rs @@ -194,10 +194,10 @@ impl SclkCollection { let path = entry.path(); if path.is_file() { let filename = path.to_str().unwrap(); - if filename.to_lowercase().ends_with(".tsc") { - if let Err(err) = self.load_file(filename) { - eprintln!("Failed to load SCLK file {filename}: {err}"); - } + if filename.to_lowercase().ends_with(".tsc") + && let Err(err) = self.load_file(filename) + { + eprintln!("Failed to load SCLK file {filename}: {err}"); } } }); diff --git a/src/kete_core/src/spice/spk.rs b/src/kete_core/src/spice/spk.rs index 8b30965..60d0b1e 100644 --- a/src/kete_core/src/spice/spk.rs +++ b/src/kete_core/src/spice/spk.rs @@ -107,7 +107,7 @@ impl SpkCollection { } } } - Err(Error::DAFLimits(format!( + Err(Error::ExceedsLimits(format!( "Object ({id}) does not have an SPK record for the target JD." ))) } @@ -269,7 +269,7 @@ impl SpkCollection { if let Some((v, _)) = result { Ok(v.iter().skip(1).map(|x| x.1).collect()) } else { - Err(Error::DAFLimits(format!( + Err(Error::ExceedsLimits(format!( "SPK files are missing information to be able to map from obj {start} to obj {goal}" ))) } @@ -396,10 +396,10 @@ impl SpkCollection { let path = entry.path(); if path.is_file() { let filename = path.to_str().unwrap(); - if filename.to_lowercase().ends_with(".bsp") { - if let Err(err) = self.load_file(filename) { - eprintln!("Failed to load SPK file {filename}: {err}"); - } + if filename.to_lowercase().ends_with(".bsp") + && let Err(err) = self.load_file(filename) + { + eprintln!("Failed to load SPK file {filename}: {err}"); } } }); diff --git a/src/kete_core/src/spice/spk_segments.rs b/src/kete_core/src/spice/spk_segments.rs index 90287e2..8ba5e18 100644 --- a/src/kete_core/src/spice/spk_segments.rs +++ b/src/kete_core/src/spice/spk_segments.rs @@ -117,7 +117,7 @@ impl SpkSegment { // this is faster than calling contains, probably because the || instead of && if jds < arr_ref.jds_start || jds > arr_ref.jds_end { - return Err(Error::DAFLimits( + return Err(Error::ExceedsLimits( "JD is not present in this record.".to_string(), )); }