diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index a2695678..20bb97f3 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -43,6 +43,7 @@ jobs: wget -O data/de430.bsp http://public-data.nyxspace.com/anise/de430.bsp wget -O data/de440s.bsp http://public-data.nyxspace.com/anise/de440s.bsp wget -O data/de440.bsp http://public-data.nyxspace.com/anise/de440.bsp + wget -O data/de440_type3.bsp http://public-data.nyxspace.com/anise/de440_type3.bsp wget -O data/pck08.pca http://public-data.nyxspace.com/anise/v0.4/pck08.pca wget -O data/pck11.pca http://public-data.nyxspace.com/anise/v0.4/pck11.pca wget -O data/moon_fk.epa http://public-data.nyxspace.com/anise/v0.4/moon_fk.epa @@ -112,6 +113,7 @@ jobs: wget -O data/de430.bsp http://public-data.nyxspace.com/anise/de430.bsp wget -O data/de440s.bsp http://public-data.nyxspace.com/anise/de440s.bsp wget -O data/de440.bsp http://public-data.nyxspace.com/anise/de440.bsp + wget -O data/de440_type3.bsp http://public-data.nyxspace.com/anise/de440_type3.bsp wget -O data/pck08.pca http://public-data.nyxspace.com/anise/v0.4/pck08.pca wget -O data/pck11.pca http://public-data.nyxspace.com/anise/v0.4/pck11.pca wget -O data/gmat-hermite.bsp http://public-data.nyxspace.com/anise/ci/gmat-hermite.bsp @@ -182,6 +184,7 @@ jobs: wget -O data/de430.bsp http://public-data.nyxspace.com/anise/de430.bsp wget -O data/de440s.bsp http://public-data.nyxspace.com/anise/de440s.bsp wget -O data/de440.bsp http://public-data.nyxspace.com/anise/de440.bsp + wget -O data/de440_type3.bsp http://public-data.nyxspace.com/anise/de440_type3.bsp wget -O data/pck08.pca http://public-data.nyxspace.com/anise/v0.4/pck08.pca wget -O data/pck11.pca http://public-data.nyxspace.com/anise/v0.4/pck11.pca wget -O data/gmat-hermite.bsp http://public-data.nyxspace.com/anise/ci/gmat-hermite.bsp @@ -214,6 +217,7 @@ jobs: cargo llvm-cov test --no-report validate_bpc_to_iau_rotations -- --nocapture --ignored cargo llvm-cov test --no-report validate_jplde_de440s_no_aberration --features spkezr_validation -- --nocapture --ignored cargo llvm-cov test --no-report validate_jplde_de440s_aberration_lt --features spkezr_validation -- --nocapture --ignored + cargo llvm-cov test --no-report validate_jplde_de440_type3_no_aberration --features spkezr_validation -- --nocapture --ignored cargo llvm-cov test --no-report validate_hermite_type13_from_gmat --features spkezr_validation -- --nocapture --ignored cargo llvm-cov test --no-report validate_lagrange_type9_with_varying_segment_sizes --features spkezr_validation -- --nocapture --ignored cargo llvm-cov test --no-report ut_embed --features embed_ephem diff --git a/anise/src/ephemerides/translate_to_parent.rs b/anise/src/ephemerides/translate_to_parent.rs index 710f7d50..20c529f4 100644 --- a/anise/src/ephemerides/translate_to_parent.rs +++ b/anise/src/ephemerides/translate_to_parent.rs @@ -17,7 +17,9 @@ use crate::ephemerides::EphemInterpolationSnafu; use crate::hifitime::Epoch; use crate::math::cartesian::CartesianState; use crate::math::Vector3; -use crate::naif::daf::datatypes::{HermiteSetType13, LagrangeSetType9, Type2ChebyshevSet}; +use crate::naif::daf::datatypes::{ + HermiteSetType13, LagrangeSetType9, Type2ChebyshevSet, Type3ChebyshevSet, +}; use crate::naif::daf::{DAFError, DafDataType, NAIFDataSet, NAIFSummaryRecord}; use crate::prelude::Frame; @@ -64,6 +66,16 @@ impl Almanac { data.evaluate(epoch, summary) .context(EphemInterpolationSnafu)? } + DafDataType::Type3ChebyshevSextuplet => { + let data = + spk_data + .nth_data::(idx_in_spk) + .context(SPKSnafu { + action: "fetching data for interpolation", + })?; + data.evaluate(epoch, summary) + .context(EphemInterpolationSnafu)? + } DafDataType::Type9LagrangeUnequalStep => { let data = spk_data .nth_data::(idx_in_spk) diff --git a/anise/src/math/interpolation/chebyshev.rs b/anise/src/math/interpolation/chebyshev.rs index a7bbeb1e..11e0b433 100644 --- a/anise/src/math/interpolation/chebyshev.rs +++ b/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 { + // 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) +} diff --git a/anise/src/math/interpolation/mod.rs b/anise/src/math/interpolation/mod.rs index 61d7977b..fabdee67 100644 --- a/anise/src/math/interpolation/mod.rs +++ b/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; diff --git a/anise/src/naif/daf/datatypes/chebyshev.rs b/anise/src/naif/daf/datatypes/chebyshev.rs index 441eaa69..32911b9c 100644 --- a/anise/src/naif/daf/datatypes/chebyshev.rs +++ b/anise/src/naif/daf/datatypes/chebyshev.rs @@ -271,51 +271,6 @@ impl<'a> NAIFDataRecord<'a> for Type2ChebyshevRecord<'a> { } } -#[derive(PartialEq)] -pub struct Type3ChebyshevRecord<'a> { - pub midpoint: Epoch, - pub radius: Duration, - pub x_coeffs: &'a [f64], - pub y_coeffs: &'a [f64], - pub z_coeffs: &'a [f64], - pub vx_coeffs: &'a [f64], - pub vy_coeffs: &'a [f64], - pub vz_coeffs: &'a [f64], -} - -impl<'a> fmt::Display for Type3ChebyshevRecord<'a> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!( - f, - "start: {}\tend: {}\nx: {:?}\ny: {:?}\nz: {:?}\nvx: {:?}\nvy: {:?}\nvz: {:?}", - self.midpoint - self.radius, - self.midpoint + self.radius, - self.x_coeffs, - self.y_coeffs, - self.z_coeffs, - self.vx_coeffs, - self.vy_coeffs, - self.vz_coeffs - ) - } -} - -impl<'a> NAIFDataRecord<'a> for Type3ChebyshevRecord<'a> { - fn from_slice_f64(slice: &'a [f64]) -> Self { - let num_coeffs = (slice.len() - 2) / 6; - Self { - midpoint: Epoch::from_et_seconds(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..], - } - } -} - #[cfg(test)] mod chebyshev_ut { use crate::{ diff --git a/anise/src/naif/daf/datatypes/chebyshev3.rs b/anise/src/naif/daf/datatypes/chebyshev3.rs new file mode 100644 index 00000000..e0f0082c --- /dev/null +++ b/anise/src/naif/daf/datatypes/chebyshev3.rs @@ -0,0 +1,380 @@ +/* + * ANISE Toolkit + * Copyright (C) 2021-onward Christopher Rabotin et al. (cf. AUTHORS.md) + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at https://mozilla.org/MPL/2.0/. + * + * Documentation: https://nyxspace.com/ + */ + +use core::fmt; +use hifitime::{Duration, Epoch, TimeUnits}; +use snafu::{ensure, ResultExt}; + +use crate::{ + errors::{DecodingError, IntegrityError, TooFewDoublesSnafu}, + math::{ + interpolation::{chebyshev_eval_poly, InterpDecodingSnafu, InterpolationError}, + Vector3, + }, + naif::daf::{NAIFDataRecord, NAIFDataSet, NAIFSummaryRecord}, +}; + +#[derive(PartialEq)] +pub struct Type3ChebyshevSet<'a> { + pub init_epoch: Epoch, + pub interval_length: Duration, + pub rsize: usize, + pub num_records: usize, + pub record_data: &'a [f64], +} + +impl<'a> Type3ChebyshevSet<'a> { + pub fn degree(&self) -> usize { + (self.rsize - 2) / 6 - 1 + } + + fn spline_idx( + &self, + epoch: Epoch, + summary: &S, + ) -> Result { + if epoch < summary.start_epoch() - 1_i64.nanoseconds() + || epoch > summary.end_epoch() + 1_i64.nanoseconds() + { + // No need to go any further. + return Err(InterpolationError::NoInterpolationData { + req: epoch, + start: summary.start_epoch(), + end: summary.end_epoch(), + }); + } + + let window_duration_s = self.interval_length.to_seconds(); + + let ephem_start_delta_s = epoch.to_et_seconds() - summary.start_epoch_et_s(); + + Ok(((ephem_start_delta_s / window_duration_s) as usize + 1).min(self.num_records)) + } +} + +impl<'a> fmt::Display for Type3ChebyshevSet<'a> { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!( + f, + "start: {:E}\tlength: {}\trsize: {}\tnum_records: {}\tlen data: {}", + self.init_epoch, + self.interval_length, + self.rsize, + self.num_records, + self.record_data.len() + ) + } +} + +impl<'a> NAIFDataSet<'a> for Type3ChebyshevSet<'a> { + type StateKind = (Vector3, Vector3); + type RecordKind = Type3ChebyshevRecord<'a>; + const DATASET_NAME: &'static str = "Chebyshev Type 3"; + + fn from_f64_slice(slice: &'a [f64]) -> Result { + ensure!( + slice.len() >= 5, + TooFewDoublesSnafu { + dataset: Self::DATASET_NAME, + need: 5_usize, + got: slice.len() + } + ); + // For this kind of record, the data is stored at the very end of the dataset + let seconds_since_j2000 = slice[slice.len() - 4]; + if !seconds_since_j2000.is_finite() { + return Err(DecodingError::Integrity { + source: IntegrityError::SubNormal { + dataset: Self::DATASET_NAME, + variable: "seconds since J2000 ET", + }, + }); + } + + let start_epoch = Epoch::from_et_seconds(seconds_since_j2000); + + let interval_length_s = slice[slice.len() - 3]; + if !interval_length_s.is_finite() { + return Err(DecodingError::Integrity { + source: IntegrityError::SubNormal { + dataset: Self::DATASET_NAME, + variable: "interval length in seconds", + }, + }); + } else if interval_length_s <= 0.0 { + return Err(DecodingError::Integrity { + source: IntegrityError::InvalidValue { + dataset: Self::DATASET_NAME, + variable: "interval length in seconds", + value: interval_length_s, + reason: "must be strictly greater than zero", + }, + }); + } + + let interval_length = interval_length_s.seconds(); + let rsize = slice[slice.len() - 2] as usize; + let num_records = slice[slice.len() - 1] as usize; + + Ok(Self { + init_epoch: start_epoch, + interval_length, + rsize, + num_records, + record_data: &slice[0..slice.len() - 4], + }) + } + + fn nth_record(&self, n: usize) -> Result { + Ok(Self::RecordKind::from_slice_f64( + self.record_data + .get(n * self.rsize..(n + 1) * self.rsize) + .ok_or(DecodingError::InaccessibleBytes { + start: n * self.rsize, + end: (n + 1) * self.rsize, + size: self.record_data.len(), + })?, + )) + } + + fn evaluate( + &self, + epoch: Epoch, + summary: &S, + ) -> Result<(Vector3, Vector3), InterpolationError> { + let spline_idx = self.spline_idx(epoch, summary)?; + + let window_duration_s = self.interval_length.to_seconds(); + let radius_s = window_duration_s / 2.0; + + let record = self + .nth_record(spline_idx - 1) + .context(InterpDecodingSnafu)?; + + let normalized_time = (epoch.to_et_seconds() - record.midpoint_et_s) / radius_s; + + let mut state = Vector3::zeros(); + let mut rate = Vector3::zeros(); + + for (cno, coeffs) in [record.x_coeffs, record.y_coeffs, record.z_coeffs] + .iter() + .enumerate() + { + 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_poly(normalized_time, coeffs, epoch, self.degree())?; + rate[cno] = val; + } + + Ok((state, rate)) + } + + fn check_integrity(&self) -> Result<(), IntegrityError> { + // Verify that none of the data is invalid once when we load it. + for val in self.record_data { + if !val.is_finite() { + return Err(IntegrityError::SubNormal { + dataset: Self::DATASET_NAME, + variable: "one of the record data", + }); + } + } + + Ok(()) + } + + fn truncate( + mut self, + summary: &S, + new_start: Option, + new_end: Option, + ) -> Result { + let start_idx = if let Some(start) = new_start { + self.spline_idx(start, summary)? - 1 + } else { + 0 + }; + + let end_idx = if let Some(end) = new_end { + self.spline_idx(end, summary)? + } else { + self.num_records - 1 + }; + + self.record_data = &self.record_data[start_idx * self.rsize..(end_idx + 1) * self.rsize]; + self.num_records = (self.record_data.len() / self.rsize) - 1; + self.init_epoch = self.nth_record(0).unwrap().midpoint_epoch() - 0.5 * self.interval_length; + + Ok(self) + } + + /// Builds the DAF array representing a Chebyshev Type 2 interpolation set. + fn to_f64_daf_vec(&self) -> Result, InterpolationError> { + let mut data = self.record_data.to_vec(); + data.push(self.init_epoch.to_et_seconds()); + data.push(self.interval_length.to_seconds()); + data.push(self.rsize as f64); + data.push(self.num_records as f64); + + Ok(data) + } +} + +#[derive(Debug, PartialEq)] +pub struct Type3ChebyshevRecord<'a> { + pub midpoint_et_s: f64, + pub radius: Duration, + pub x_coeffs: &'a [f64], + pub y_coeffs: &'a [f64], + pub z_coeffs: &'a [f64], + pub vx_coeffs: &'a [f64], + pub vy_coeffs: &'a [f64], + pub vz_coeffs: &'a [f64], +} + +impl<'a> Type3ChebyshevRecord<'a> { + pub fn midpoint_epoch(&self) -> Epoch { + Epoch::from_et_seconds(self.midpoint_et_s) + } +} + +impl<'a> fmt::Display for Type3ChebyshevRecord<'a> { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!( + f, + "start: {}\tend: {}\nx: {:?}\ny: {:?}\nz: {:?}\nvx: {:?}\nvy: {:?}\nvz: {:?}", + self.midpoint_epoch() - self.radius, + self.midpoint_epoch() + self.radius, + self.x_coeffs, + self.y_coeffs, + self.z_coeffs, + self.vx_coeffs, + self.vy_coeffs, + self.vz_coeffs + ) + } +} + +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..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..], + } + } +} + +#[cfg(test)] +mod chebyshev_ut { + use crate::{ + errors::{DecodingError, IntegrityError}, + naif::daf::NAIFDataSet, + }; + + use super::Type3ChebyshevSet; + + #[test] + fn too_small() { + if Type3ChebyshevSet::from_f64_slice(&[0.1, 0.2, 0.3, 0.4]) + != Err(DecodingError::TooFewDoubles { + dataset: "Chebyshev Type 3", + got: 4, + need: 5, + }) + { + panic!("test failure"); + } + } + + #[test] + fn subnormal() { + match Type3ChebyshevSet::from_f64_slice(&[0.0, f64::INFINITY, 0.0, 0.0, 0.0]) { + Ok(_) => panic!("test failed on invalid init_epoch"), + Err(e) => { + assert_eq!( + e, + DecodingError::Integrity { + source: IntegrityError::SubNormal { + dataset: "Chebyshev Type 3", + variable: "seconds since J2000 ET", + }, + } + ); + } + } + + match Type3ChebyshevSet::from_f64_slice(&[0.0, 0.0, f64::INFINITY, 0.0, 0.0]) { + Ok(_) => panic!("test failed on invalid interval_length"), + Err(e) => { + assert_eq!( + e, + DecodingError::Integrity { + source: IntegrityError::SubNormal { + dataset: "Chebyshev Type 3", + variable: "interval length in seconds", + }, + } + ); + } + } + + match Type3ChebyshevSet::from_f64_slice(&[0.0, 0.0, -1e-16, 0.0, 0.0]) { + Ok(_) => panic!("test failed on invalid interval_length"), + Err(e) => { + assert_eq!( + e, + DecodingError::Integrity { + source: IntegrityError::InvalidValue { + dataset: "Chebyshev Type 3", + variable: "interval length in seconds", + value: -1e-16, + reason: "must be strictly greater than zero" + }, + } + ); + } + } + + // Load a slice whose metadata is OK but the record data is not + let dataset = + Type3ChebyshevSet::from_f64_slice(&[f64::INFINITY, 0.0, 2e-16, 0.0, 0.0]).unwrap(); + match dataset.check_integrity() { + Ok(_) => panic!("test failed on invalid interval_length"), + Err(e) => { + assert_eq!( + e, + IntegrityError::SubNormal { + dataset: "Chebyshev Type 3", + variable: "one of the record data", + }, + ); + } + } + } +} diff --git a/anise/src/naif/daf/datatypes/mod.rs b/anise/src/naif/daf/datatypes/mod.rs index dd2cedc9..eb442183 100644 --- a/anise/src/naif/daf/datatypes/mod.rs +++ b/anise/src/naif/daf/datatypes/mod.rs @@ -9,10 +9,12 @@ */ pub mod chebyshev; +pub mod chebyshev3; pub mod hermite; pub mod lagrange; pub mod posvel; pub use chebyshev::*; +pub use chebyshev3::*; pub use hermite::*; pub use lagrange::*; diff --git a/anise/tests/almanac/mod.rs b/anise/tests/almanac/mod.rs index dbc03dca..08975601 100644 --- a/anise/tests/almanac/mod.rs +++ b/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}, }; @@ -77,3 +77,31 @@ fn test_state_transformation() { assert_eq!(orig_state, from_state_itrf93_to_eme2k); } + +#[test] +fn test_type3_state_transformation() { + // Load BSP and BPC + let ctx = Almanac::default(); + + let spk = SPK::load("../data/de440_type3.bsp").unwrap(); + + let almanac = ctx + .with_spk(spk) + .unwrap() + .load("../data/pck08.pca") + .unwrap(); + + let epoch = Epoch::from_str("2021-10-29 12:34:56 TDB").unwrap(); + + let to_parent = almanac.translate_to_parent(EARTH_J2000, epoch).unwrap(); + + println!("{to_parent}"); + + // Ensure that we can query the type 3 chebyshev DE440 file + + let state = almanac + .translate(EARTH_J2000, SUN_J2000, epoch, None) + .expect("type 3 chebyshev could not be queried"); + + println!("{state:x}"); +} diff --git a/anise/tests/ephemerides/validation/mod.rs b/anise/tests/ephemerides/validation/mod.rs index 1e362fec..000dd8dc 100644 --- a/anise/tests/ephemerides/validation/mod.rs +++ b/anise/tests/ephemerides/validation/mod.rs @@ -9,6 +9,7 @@ */ mod type02_chebyshev_jpl_de; +mod type03_chebyshev_jpl_de; mod type09_lagrange; mod type13_hermite; diff --git a/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs b/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs new file mode 100644 index 00000000..34e77f6d --- /dev/null +++ b/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs @@ -0,0 +1,34 @@ +/* + * ANISE Toolkit + * Copyright (C) 2021-onward Christopher Rabotin et al. (cf. AUTHORS.md) + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at https://mozilla.org/MPL/2.0/. + * + * Documentation: https://nyxspace.com/ + */ + +use super::{compare::*, validate::Validation}; + +#[ignore = "Requires Rust SPICE -- must be executed serially"] +#[test] +fn validate_jplde_de440_type3_no_aberration() { + let file_name = "spk-type3-validation-de440".to_string(); + let comparator = CompareEphem::new( + vec!["../data/de440_type3.bsp".to_string()], + file_name.clone(), + 1_000, + None, + ); + + let err_count = comparator.run(); + + assert_eq!(err_count, 0, "None of the queries should fail!"); + + let validator = Validation { + file_name, + ..Default::default() + }; + + validator.validate(); +} diff --git a/data/de440_type3.bsp b/data/de440_type3.bsp new file mode 100644 index 00000000..473880e1 --- /dev/null +++ b/data/de440_type3.bsp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7e655129d94609bd7d22e792d8e7b2772e314cc42f03ed04f2ea46e0fa75b993 +size 360009728