Implement origin change

Cargo.toml

@@ -24,17 +24,16 @@ lox-time = { path = "crates/lox-time", version = "0.1.0-alpha.0" }
 csv = "1.3.0"
 divan = "0.1.14"
 dyn-clone = "1.0.17"
-fast-float = "0.2.0"
 fast_polynomial = "0.1.0"
 float_eq = "1.0.1"
 glam = "0.28.0"
 itertools = "0.13.0"
 libm = "0.2.8"
 nom = "7.1.3"
 num = "0.4.1"
-numpy = "0.21.0"
+numpy = "0.22.1"
 proptest = "1.4.0"
-pyo3 = "0.21.1"
+pyo3 = "0.22.6"
 quick-xml = { version = "0.31.0", features = ["serde", "serialize"] }
 regex = "1.10.4"
 rstest = "0.21.0"

crates/lox-ephem/Cargo.toml

@@ -9,5 +9,12 @@ authors.workspace = true
 repository.workspace = true
 
 [dependencies]
-nom.workspace = true
 lox-math.workspace = true
+
+nom.workspace = true
+pyo3 = { workspace = true, optional = true }
+numpy = { workspace = true, optional = true }
+thiserror.workspace = true
+
+[features]
+python = ["dep:pyo3", "dep:numpy"]

crates/lox-ephem/src/lib.rs

@@ -1 +1,77 @@
+use lox_math::types::julian_dates::Epoch;
+
+#[cfg(feature = "python")]
+pub mod python;
 pub mod spk;
+
+pub(crate) type Position = (f64, f64, f64);
+pub(crate) type Velocity = (f64, f64, f64);
+pub(crate) type Body = i32;
+
+pub trait Ephemeris {
+    type Error: std::error::Error;
+
+    fn position(&self, epoch: Epoch, origin: Body, target: Body) -> Result<Position, Self::Error>;
+    fn velocity(&self, epoch: Epoch, origin: Body, target: Body) -> Result<Velocity, Self::Error>;
+    fn state(
+        &self,
+        epoch: Epoch,
+        origin: Body,
+        target: Body,
+    ) -> Result<(Position, Velocity), Self::Error>;
+}
+
+fn ancestors(id: i32) -> Vec<i32> {
+    let mut ancestors = vec![id];
+    let mut current = id;
+    while current != 0 {
+        current /= 100;
+        ancestors.push(current);
+    }
+    ancestors
+}
+
+pub fn path_from_ids(origin: i32, target: i32) -> Vec<i32> {
+    let ancestors_origin = ancestors(origin);
+    let ancestors_target = ancestors(target);
+    let n = ancestors_target.len();
+    let mut path = ancestors_origin;
+
+    ancestors_target
+        .into_iter()
+        .take(n - 1)
+        .rev()
+        .for_each(|id| path.push(id));
+
+    if *path.first().unwrap() != 0 && *path.last().unwrap() != 0 {
+        let idx = path.iter().position(|&id| id == 0).unwrap();
+        if path[idx - 1] == path[idx + 1] {
+            let common_ancestor = vec![path[idx - 1]];
+            path.splice((idx - 1)..=(idx + 1), common_ancestor);
+        }
+    }
+
+    path
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_ancestors() {
+        assert_eq!(ancestors(0), vec![0]);
+        assert_eq!(ancestors(3), vec![3, 0]);
+        assert_eq!(ancestors(399), vec![399, 3, 0]);
+    }
+
+    #[test]
+    fn test_path_from_ids() {
+        assert_eq!(path_from_ids(399, 499), [399, 3, 0, 4, 499]);
+        assert_eq!(path_from_ids(399, 0), [399, 3, 0]);
+        assert_eq!(path_from_ids(0, 399), [0, 3, 399]);
+        assert_eq!(path_from_ids(399, 3), [399, 3]);
+        assert_eq!(path_from_ids(3, 399), [3, 399]);
+        assert_eq!(path_from_ids(399, 301), [399, 3, 301]);
+    }
+}

crates/lox-ephem/src/python.rs

@@ -0,0 +1,22 @@
+use pyo3::{exceptions::PyValueError, pyclass, pymethods, PyErr, PyResult};
+
+use crate::spk::parser::{parse_daf_spk, DafSpkError, Spk};
+
+impl From<DafSpkError> for PyErr {
+    fn from(err: DafSpkError) -> Self {
+        PyValueError::new_err(err.to_string())
+    }
+}
+
+#[pyclass(name = "SPK", module = "lox_space", frozen)]
+pub struct PySpk(pub Spk);
+
+#[pymethods]
+impl PySpk {
+    #[new]
+    fn new(path: &str) -> PyResult<Self> {
+        let data = std::fs::read(path)?;
+        let spk = parse_daf_spk(&data)?;
+        Ok(PySpk(spk))
+    }
+}

crates/lox-ephem/src/spk/api.rs

@@ -10,11 +10,9 @@ use std::collections::HashMap;
 
 use lox_math::types::julian_dates::Epoch;
 
-use super::parser::{DafSpkError, Spk, SpkSegment, SpkType2Array, SpkType2Coefficients};
+use crate::{Body, Ephemeris, Position, Velocity};
 
-type Position = (f64, f64, f64);
-type Velocity = (f64, f64, f64);
-type Body = i32;
+use super::parser::{DafSpkError, Spk, SpkSegment, SpkType2Array, SpkType2Coefficients};
 
 impl Spk {
     fn find_segment(
@@ -105,13 +103,12 @@ impl Spk {
 
         Ok((coefficients, record))
     }
+}
 
-    pub fn position(
-        &self,
-        epoch: Epoch,
-        origin: Body,
-        target: Body,
-    ) -> Result<Position, DafSpkError> {
+impl Ephemeris for Spk {
+    type Error = DafSpkError;
+
+    fn position(&self, epoch: Epoch, origin: Body, target: Body) -> Result<Position, DafSpkError> {
         let (segment, sign) = self.find_segment(origin, target)?;
 
         if epoch < segment.initial_epoch || epoch > segment.final_epoch {
@@ -141,12 +138,7 @@ impl Spk {
         Ok((x, y, z))
     }
 
-    pub fn velocity(
-        &self,
-        epoch: Epoch,
-        origin: Body,
-        target: Body,
-    ) -> Result<Velocity, DafSpkError> {
+    fn velocity(&self, epoch: Epoch, origin: Body, target: Body) -> Result<Velocity, DafSpkError> {
         let (segment, sign) = self.find_segment(origin, target)?;
 
         if epoch < segment.initial_epoch || epoch > segment.final_epoch {
@@ -197,7 +189,7 @@ impl Spk {
         Ok((x, y, z))
     }
 
-    pub fn state(
+    fn state(
         &self,
         epoch: Epoch,
         origin: Body,

crates/lox-ephem/src/spk/parser.rs

@@ -12,6 +12,7 @@ use std::iter::zip;
 use nom::bytes::complete as nb;
 use nom::error::ErrorKind;
 use nom::number::complete as nn;
+use thiserror::Error;
 
 type BodyId = i32;
 
@@ -47,17 +48,19 @@ pub struct DafSummary {
     pub final_address: usize,
 }
 
-#[derive(Debug, PartialEq)]
+#[derive(Debug, PartialEq, Error)]
 pub enum DafSpkError {
-    // The data type integer value does not match the ones in the spec
+    #[error("the data type integer value does not match the ones in the spec")]
     InvalidSpkSegmentDataType,
-    // The number of DAF components does not match the SPK specification
+    #[error("the number of DAF components does not match the SPK specification")]
     UnexpectedNumberOfComponents,
+    #[error("unable to parse")]
     UnableToParse,
+    #[error("unsupported SPK type {data_type}")]
     UnsupportedSpkArrayType { data_type: i32 },
-    // Unable to find the segment for a given center body and target body
+    #[error("unable to find the segment for a given center body and target body")]
     UnableToFindMatchingSegment,
-    // Unable to find record for a given date
+    #[error("unable to find record for a given date")]
     UnableToFindMatchingRecord,
 }
 

crates/lox-io/Cargo.toml

@@ -9,19 +9,17 @@ authors.workspace = true
 repository.workspace = true
 
 [dependencies]
+lox-derive.workspace = true
 lox-math.workspace = true
 
 csv.workspace = true
 nom.workspace = true
-serde.workspace = true
-thiserror.workspace = true
-
 quick-xml.workspace = true
-serde_json.workspace = true
-serde-aux.workspace = true
 regex.workspace = true
-fast-float.workspace = true
-lox-derive.workspace = true
+serde-aux.workspace = true
+serde.workspace = true
+serde_json.workspace = true
+thiserror.workspace = true
 
 [dev-dependencies]
 rstest.workspace = true

crates/lox-io/src/ndm.rs

@@ -13,10 +13,10 @@
 //! input that they accept. Some relaxations:
 //!
 //! - The KVN floating point numbers defined in the specification can have only
-//! one character in the integer part of the number, but we accept any regular
-//! float number.
+//!   one character in the integer part of the number, but we accept any regular
+//!   float number.
 //! - The KVN strings are defined in the specification as being either only
-//! lower-case or only upper-case, but we accept any combination of cases.
+//!   lower-case or only upper-case, but we accept any combination of cases.
 //!
 //! The XML deserializer does not perform any validation on the schema types
 //! defined (e.g. non-positive double, lat-long, angle). The validation only

crates/lox-io/src/ndm/kvn/parser.rs

@@ -317,7 +317,9 @@ fn parse_kvn_covariance_matrix_line<'a, T: Iterator<Item = &'a str> + ?Sized>(
     let result: Result<Vec<f64>, _> = next_line
         .split_whitespace()
         .map(|matrix_element| {
-            fast_float::parse(matrix_element.trim())
+            matrix_element
+                .trim()
+                .parse::<f64>()
                 .map_err(|_| KvnCovarianceMatrixParserErr::InvalidFormat { input: next_line })
         })
         .collect();
@@ -558,8 +560,9 @@ pub fn parse_kvn_numeric_line(
     let value = captures.name("value").unwrap().as_str();
     let unit = captures.name("unit").map(|x| x.as_str().to_string());
 
-    let value =
-        fast_float::parse(value).map_err(|_| KvnNumberParserErr::InvalidFormat { input })?;
+    let value = value
+        .parse::<f64>()
+        .map_err(|_| KvnNumberParserErr::InvalidFormat { input })?;
 
     Ok(KvnValue { value, unit })
 }

crates/lox-orbits/Cargo.toml

@@ -10,6 +10,7 @@ repository.workspace = true
 
 [dependencies]
 lox-bodies.workspace = true
+lox-ephem.workspace = true
 lox-time.workspace = true
 lox-math.workspace = true
 

crates/lox-orbits/src/analysis.rs

@@ -130,15 +130,7 @@ mod tests {
         let actual: Vec<Radians> = gs
             .times()
             .iter()
-            .map(|t| {
-                elevation(
-                    t.clone(),
-                    &frame,
-                    &gs,
-                    &sc,
-                    &NoOpFrameTransformationProvider,
-                )
-            })
+            .map(|t| elevation(*t, &frame, &gs, &sc, &NoOpFrameTransformationProvider))
             .collect();
         for (actual, expected) in actual.iter().zip(expected.iter()) {
             assert_close!(actual, expected, 1e-1);

crates/lox-orbits/src/propagators/semi_analytical.rs

@@ -178,7 +178,7 @@ mod tests {
         let s0 = k0.to_cartesian();
         let t1 = time + k0.orbital_period();
 
-        let propagator = Vallado::new(s0.clone());
+        let propagator = Vallado::new(s0);
         let s1 = propagator
             .propagate(t1)
             .expect("propagator should converge");
@@ -222,7 +222,7 @@ mod tests {
         let period = k0.orbital_period();
         let t_end = period.to_decimal_seconds().ceil() as i64;
         let steps = TimeDelta::range(0..=t_end).map(|dt| time + dt);
-        let trajectory = Vallado::new(s0.clone()).propagate_all(steps).unwrap();
+        let trajectory = Vallado::new(s0).propagate_all(steps).unwrap();
         let s1 = trajectory.interpolate(period);
         let k1 = s1.to_keplerian();