Fix Chebyshev type 3

anise/src/math/interpolation/chebyshev.rs

@@ -57,3 +57,33 @@ pub fn chebyshev_eval(
     let deriv = (w[0] + normalized_time * dw[0] - dw[1]) / spline_radius_s;
     Ok((val, deriv))
 }
+
+/// Attempts to evaluate a Chebyshev polynomial given the coefficients, returning only the value
+///
+/// # Notes
+/// 1. At this point, the splines are expected to be in Chebyshev format and no verification is done.
+pub fn chebyshev_eval_poly(
+    normalized_time: f64,
+    spline_coeffs: &[f64],
+    eval_epoch: Epoch,
+    degree: usize,
+) -> Result<f64, InterpolationError> {
+    // Workspace array
+    let mut w = [0.0_f64; 3];
+
+    for j in (2..=degree + 1).rev() {
+        w[2] = w[1];
+        w[1] = w[0];
+        w[0] = (spline_coeffs
+            .get(j - 1)
+            .ok_or(InterpolationError::MissingInterpolationData { epoch: eval_epoch })?)
+            + (2.0 * normalized_time * w[1] - w[2]);
+    }
+
+    let val = (spline_coeffs
+        .first()
+        .ok_or(InterpolationError::MissingInterpolationData { epoch: eval_epoch })?)
+        + (normalized_time * w[0]);
+
+    Ok(val)
+}

anise/src/math/interpolation/mod.rs

@@ -12,7 +12,7 @@ mod chebyshev;
 mod hermite;
 mod lagrange;
 
-pub use chebyshev::chebyshev_eval;
+pub use chebyshev::{chebyshev_eval, chebyshev_eval_poly};
 pub use hermite::hermite_eval;
 use hifitime::Epoch;
 pub use lagrange::lagrange_eval;

anise/src/naif/daf/datatypes/chebyshev3.rs

@@ -15,7 +15,7 @@ use snafu::{ensure, ResultExt};
 use crate::{
     errors::{DecodingError, IntegrityError, TooFewDoublesSnafu},
     math::{
-        interpolation::{chebyshev_eval, InterpDecodingSnafu, InterpolationError},
+        interpolation::{chebyshev_eval_poly, InterpDecodingSnafu, InterpolationError},
         Vector3,
     },
     naif::daf::{NAIFDataRecord, NAIFDataSet, NAIFSummaryRecord},
@@ -32,7 +32,7 @@ pub struct Type3ChebyshevSet<'a> {
 
 impl<'a> Type3ChebyshevSet<'a> {
     pub fn degree(&self) -> usize {
-        (self.rsize - 2) / 3 - 1
+        (self.rsize - 2) / 6 - 1
     }
 
     fn spline_idx<S: NAIFSummaryRecord>(
@@ -154,7 +154,6 @@ impl<'a> NAIFDataSet<'a> for Type3ChebyshevSet<'a> {
         let window_duration_s = self.interval_length.to_seconds();
         let radius_s = window_duration_s / 2.0;
 
-        // Now, build the X, Y, Z data from the record data.
         let record = self
             .nth_record(spline_idx - 1)
             .context(InterpDecodingSnafu)?;
@@ -168,15 +167,15 @@ impl<'a> NAIFDataSet<'a> for Type3ChebyshevSet<'a> {
             .iter()
             .enumerate()
         {
-            let (val, _) = chebyshev_eval(normalized_time, coeffs, radius_s, epoch, self.degree())?;
+            let val = chebyshev_eval_poly(normalized_time, coeffs, epoch, self.degree())?;
             state[cno] = val;
         }
 
         for (cno, coeffs) in [record.vx_coeffs, record.vy_coeffs, record.vz_coeffs]
             .iter()
             .enumerate()
         {
-            let (val, _) = chebyshev_eval(normalized_time, coeffs, radius_s, epoch, self.degree())?;
+            let val = chebyshev_eval_poly(normalized_time, coeffs, epoch, self.degree())?;
             rate[cno] = val;
         }
 
@@ -234,7 +233,7 @@ impl<'a> NAIFDataSet<'a> for Type3ChebyshevSet<'a> {
     }
 }
 
-#[derive(PartialEq)]
+#[derive(Debug, PartialEq)]
 pub struct Type3ChebyshevRecord<'a> {
     pub midpoint_et_s: f64,
     pub radius: Duration,
@@ -272,15 +271,21 @@ impl<'a> fmt::Display for Type3ChebyshevRecord<'a> {
 impl<'a> NAIFDataRecord<'a> for Type3ChebyshevRecord<'a> {
     fn from_slice_f64(slice: &'a [f64]) -> Self {
         let num_coeffs = (slice.len() - 2) / 6;
+
+        let end_x_idx = num_coeffs + 2;
+        let end_y_idx = 2 * num_coeffs + 2;
+        let end_z_idx = 3 * num_coeffs + 2;
+        let end_vx_idx = 4 * num_coeffs + 2;
+        let end_vy_idx = 5 * num_coeffs + 2;
         Self {
             midpoint_et_s: slice[0],
             radius: slice[1].seconds(),
-            x_coeffs: &slice[2..num_coeffs],
-            y_coeffs: &slice[2 + num_coeffs..num_coeffs * 2],
-            z_coeffs: &slice[2 + num_coeffs * 2..num_coeffs * 3],
-            vx_coeffs: &slice[2 + num_coeffs * 3..num_coeffs * 4],
-            vy_coeffs: &slice[2 + num_coeffs * 4..num_coeffs * 5],
-            vz_coeffs: &slice[2 + num_coeffs * 5..],
+            x_coeffs: &slice[2..end_x_idx],
+            y_coeffs: &slice[end_x_idx..end_y_idx],
+            z_coeffs: &slice[end_y_idx..end_z_idx],
+            vx_coeffs: &slice[end_z_idx..end_vx_idx],
+            vy_coeffs: &slice[end_vx_idx..end_vy_idx],
+            vz_coeffs: &slice[end_vy_idx..],
         }
     }
 }

anise/tests/almanac/mod.rs

@@ -1,6 +1,6 @@
 // Start by creating the ANISE planetary data
 use anise::{
-    constants::frames::{EARTH_ITRF93, EARTH_J2000},
+    constants::frames::{EARTH_ITRF93, EARTH_J2000, SUN_J2000},
     naif::kpl::parser::convert_tpc,
     prelude::{Aberration, Almanac, Orbit, BPC, SPK},
 };
@@ -84,49 +84,24 @@ fn test_type3_state_transformation() {
     let ctx = Almanac::default();
 
     let spk = SPK::load("../data/de440_type3.bsp").unwrap();
-    let bpc = BPC::load("../data/earth_latest_high_prec.bpc").unwrap();
 
     let almanac = ctx
         .with_spk(spk)
         .unwrap()
-        .with_bpc(bpc)
-        .unwrap()
         .load("../data/pck08.pca")
         .unwrap();
 
-    // Let's build an orbit
-    // Start by grabbing a copy of the frame.
-    let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
-    // Define an epoch
     let epoch = Epoch::from_str("2021-10-29 12:34:56 TDB").unwrap();
 
-    let orig_state = Orbit::keplerian(
-        8_191.93, 1e-6, 12.85, 306.614, 314.19, 99.887_7, epoch, eme2k,
-    );
-
-    // Transform that into another frame.
-    let state_itrf93 = almanac
-        .transform_to(orig_state, EARTH_ITRF93, Aberration::NONE)
-        .unwrap();
-
-    // Check that the frame is correctly set.
-    assert_eq!(state_itrf93.frame.ephemeris_id, EARTH_ITRF93.ephemeris_id);
-    assert_eq!(
-        state_itrf93.frame.orientation_id,
-        EARTH_ITRF93.orientation_id
-    );
+    let to_parent = almanac.translate_to_parent(EARTH_J2000, epoch).unwrap();
 
-    println!("{orig_state:x}");
-    println!("{state_itrf93:X}");
+    println!("{to_parent}");
 
-    // Convert back.
-    // Note that the Aberration correction constants are actually options!
-    let from_state_itrf93_to_eme2k = almanac
-        .transform_to(state_itrf93, EARTH_J2000, None)
-        .unwrap();
+    // Ensure that we can query the type 3 chebyshev DE440 file
 
-    println!("{from_state_itrf93_to_eme2k}");
-    println!("{from_state_itrf93_to_eme2k:x}");
+    let state = almanac
+        .translate(EARTH_J2000, SUN_J2000, epoch, None)
+        .expect("type 3 chebyshev could not be queried");
 
-    assert_eq!(orig_state, from_state_itrf93_to_eme2k);
+    println!("{state:x}");
 }

anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs

@@ -9,7 +9,6 @@
  */
 
 use super::{compare::*, validate::Validation};
-use anise::prelude::Aberration;
 
 #[ignore = "Requires Rust SPICE -- must be executed serially"]
 #[test]