From 573ad1c6cbeee159afd0732ad0e311497ea6bf39 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Tue, 17 Sep 2024 23:56:56 -0600 Subject: [PATCH 1/4] Initial support for Chebyshev type 3 --- anise/src/ephemerides/translate_to_parent.rs | 14 +- anise/src/naif/daf/datatypes/chebyshev.rs | 207 ++++++++++++++++++ anise/tests/almanac/mod.rs | 53 +++++ anise/tests/ephemerides/validation/mod.rs | 1 + .../validation/type03_chebyshev_jpl_de.rs | 35 +++ data/de440_type3.bsp | 3 + 6 files changed, 312 insertions(+), 1 deletion(-) create mode 100644 anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs create mode 100644 data/de440_type3.bsp 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/naif/daf/datatypes/chebyshev.rs b/anise/src/naif/daf/datatypes/chebyshev.rs index 441eaa69..216bb885 100644 --- a/anise/src/naif/daf/datatypes/chebyshev.rs +++ b/anise/src/naif/daf/datatypes/chebyshev.rs @@ -316,6 +316,213 @@ impl<'a> NAIFDataRecord<'a> for Type3ChebyshevRecord<'a> { } } +#[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) / 3 - 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 = Type2ChebyshevRecord<'a>; + const DATASET_NAME: &'static str = "Chebyshev Type 2"; + + 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; + + // Now, build the X, Y, Z data from the record data. + 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, deriv) = + chebyshev_eval(normalized_time, coeffs, radius_s, epoch, self.degree())?; + state[cno] = val; + rate[cno] = deriv; + } + + 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) + } +} + #[cfg(test)] mod chebyshev_ut { use crate::{ diff --git a/anise/tests/almanac/mod.rs b/anise/tests/almanac/mod.rs index dbc03dca..cd8a01db 100644 --- a/anise/tests/almanac/mod.rs +++ b/anise/tests/almanac/mod.rs @@ -77,3 +77,56 @@ 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 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 + ); + + println!("{orig_state:x}"); + println!("{state_itrf93:X}"); + + // 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(); + + println!("{from_state_itrf93_to_eme2k}"); + println!("{from_state_itrf93_to_eme2k:x}"); + + assert_eq!(orig_state, from_state_itrf93_to_eme2k); +} 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..a2e90b92 --- /dev/null +++ b/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs @@ -0,0 +1,35 @@ +/* + * 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}; +use anise::prelude::Aberration; + +#[ignore = "Requires Rust SPICE -- must be executed serially"] +#[test] +fn validate_jplde_de440_full() { + 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 From e8e2a84e47d74bcbdc84a3e3b7977be152645eca Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Wed, 18 Sep 2024 00:15:51 -0600 Subject: [PATCH 2/4] Private ChebyshevSet trait reduces code duplication --- .github/workflows/rust.yml | 4 + anise/src/naif/daf/datatypes/chebyshev.rs | 412 +++++++----------- .../validation/type03_chebyshev_jpl_de.rs | 2 +- 3 files changed, 173 insertions(+), 245 deletions(-) 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/naif/daf/datatypes/chebyshev.rs b/anise/src/naif/daf/datatypes/chebyshev.rs index 216bb885..015a17c5 100644 --- a/anise/src/naif/daf/datatypes/chebyshev.rs +++ b/anise/src/naif/daf/datatypes/chebyshev.rs @@ -30,9 +30,64 @@ pub struct Type2ChebyshevSet<'a> { pub record_data: &'a [f64], } -impl<'a> Type2ChebyshevSet<'a> { - pub fn degree(&self) -> usize { - (self.rsize - 2) / 3 - 1 +impl<'a> fmt::Display for Type2ChebyshevSet<'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() + ) + } +} + +#[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> 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() + ) + } +} + +trait ChebyshevSet<'a>: PartialEq + fmt::Display { + fn new( + init_epoch: Epoch, + interval_length: Duration, + rsize: usize, + num_records: usize, + record_data: &'a [f64], + ) -> Self; + + fn init_epoch(&self) -> Epoch; + fn interval_length(&self) -> Duration; + fn rsize(&self) -> usize; + fn num_records(&self) -> usize; + fn record_data(&self) -> &'a [f64]; + + fn set_init_epoch(&mut self, init_epoch: Epoch); + fn set_num_records(&mut self, num_records: usize); + fn set_record_data(&mut self, record_data: &'a [f64]); + + fn degree(&self) -> usize { + (self.rsize() - 2) / 3 - 1 } fn spline_idx( @@ -51,29 +106,101 @@ impl<'a> Type2ChebyshevSet<'a> { }); } - let window_duration_s = self.interval_length.to_seconds(); + 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)) + Ok(((ephem_start_delta_s / window_duration_s) as usize + 1).min(self.num_records())) } } -impl<'a> fmt::Display for Type2ChebyshevSet<'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> ChebyshevSet<'a> for Type2ChebyshevSet<'a> { + fn new( + init_epoch: Epoch, + interval_length: Duration, + rsize: usize, + num_records: usize, + record_data: &'a [f64], + ) -> Self { + Self { + init_epoch, + interval_length, + rsize, + num_records, + record_data, + } + } + + fn init_epoch(&self) -> Epoch { + self.init_epoch + } + fn interval_length(&self) -> Duration { + self.interval_length + } + fn rsize(&self) -> usize { + self.rsize + } + fn num_records(&self) -> usize { + self.num_records + } + fn record_data(&self) -> &'a [f64] { + self.record_data + } + + fn set_init_epoch(&mut self, init_epoch: Epoch) { + self.init_epoch = init_epoch; + } + fn set_num_records(&mut self, num_records: usize) { + self.num_records = num_records; + } + fn set_record_data(&mut self, record_data: &'a [f64]) { + self.record_data = record_data; + } +} +impl<'a> ChebyshevSet<'a> for Type3ChebyshevSet<'a> { + fn new( + init_epoch: Epoch, + interval_length: Duration, + rsize: usize, + num_records: usize, + record_data: &'a [f64], + ) -> Self { + Self { + init_epoch, + interval_length, + rsize, + num_records, + record_data, + } + } + + fn init_epoch(&self) -> Epoch { + self.init_epoch + } + fn interval_length(&self) -> Duration { + self.interval_length + } + fn rsize(&self) -> usize { + self.rsize + } + fn num_records(&self) -> usize { + self.num_records + } + fn record_data(&self) -> &'a [f64] { + self.record_data + } + fn set_init_epoch(&mut self, init_epoch: Epoch) { + self.init_epoch = init_epoch; + } + fn set_num_records(&mut self, num_records: usize) { + self.num_records = num_records; + } + fn set_record_data(&mut self, record_data: &'a [f64]) { + self.record_data = record_data; } } -impl<'a> NAIFDataSet<'a> for Type2ChebyshevSet<'a> { +impl<'a, T: ChebyshevSet<'a>> NAIFDataSet<'a> for T { type StateKind = (Vector3, Vector3); type RecordKind = Type2ChebyshevRecord<'a>; const DATASET_NAME: &'static str = "Chebyshev Type 2"; @@ -123,23 +250,23 @@ impl<'a> NAIFDataSet<'a> for Type2ChebyshevSet<'a> { let rsize = slice[slice.len() - 2] as usize; let num_records = slice[slice.len() - 1] as usize; - Ok(Self { - init_epoch: start_epoch, + Ok(Self::new( + start_epoch, interval_length, rsize, num_records, - record_data: &slice[0..slice.len() - 4], - }) + &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) + 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(), + start: n * self.rsize(), + end: (n + 1) * self.rsize(), + size: self.record_data().len(), })?, )) } @@ -151,7 +278,7 @@ impl<'a> NAIFDataSet<'a> for Type2ChebyshevSet<'a> { ) -> Result<(Vector3, Vector3), InterpolationError> { let spline_idx = self.spline_idx(epoch, summary)?; - let window_duration_s = self.interval_length.to_seconds(); + 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. @@ -179,7 +306,7 @@ impl<'a> NAIFDataSet<'a> for Type2ChebyshevSet<'a> { 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 { + for val in self.record_data() { if !val.is_finite() { return Err(IntegrityError::SubNormal { dataset: Self::DATASET_NAME, @@ -206,23 +333,27 @@ impl<'a> NAIFDataSet<'a> for Type2ChebyshevSet<'a> { let end_idx = if let Some(end) = new_end { self.spline_idx(end, summary)? } else { - self.num_records - 1 + 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; + self.set_record_data( + &self.record_data()[start_idx * self.rsize()..(end_idx + 1) * self.rsize()], + ); + self.set_num_records((self.record_data().len() / self.rsize()) - 1); + self.set_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); + 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) } @@ -316,213 +447,6 @@ impl<'a> NAIFDataRecord<'a> for Type3ChebyshevRecord<'a> { } } -#[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) / 3 - 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 = Type2ChebyshevRecord<'a>; - const DATASET_NAME: &'static str = "Chebyshev Type 2"; - - 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; - - // Now, build the X, Y, Z data from the record data. - 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, deriv) = - chebyshev_eval(normalized_time, coeffs, radius_s, epoch, self.degree())?; - state[cno] = val; - rate[cno] = deriv; - } - - 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) - } -} - #[cfg(test)] mod chebyshev_ut { use crate::{ diff --git a/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs b/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs index a2e90b92..705d3312 100644 --- a/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs +++ b/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs @@ -13,7 +13,7 @@ use anise::prelude::Aberration; #[ignore = "Requires Rust SPICE -- must be executed serially"] #[test] -fn validate_jplde_de440_full() { +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()], From 0c84f72f9cd296d9e1c1833c98f060b987c8eb6d Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 22 Sep 2024 20:56:42 -0600 Subject: [PATCH 3/4] Revert changes tot ype2 cheby, copied over relevant type3 cheby The abstraction barely reduced code duplication but greatly increased complexity --- anise/src/naif/daf/datatypes/chebyshev.rs | 250 ++------------ anise/src/naif/daf/datatypes/chebyshev3.rs | 375 +++++++++++++++++++++ anise/src/naif/daf/datatypes/mod.rs | 2 + 3 files changed, 414 insertions(+), 213 deletions(-) create mode 100644 anise/src/naif/daf/datatypes/chebyshev3.rs diff --git a/anise/src/naif/daf/datatypes/chebyshev.rs b/anise/src/naif/daf/datatypes/chebyshev.rs index 015a17c5..32911b9c 100644 --- a/anise/src/naif/daf/datatypes/chebyshev.rs +++ b/anise/src/naif/daf/datatypes/chebyshev.rs @@ -30,64 +30,9 @@ pub struct Type2ChebyshevSet<'a> { pub record_data: &'a [f64], } -impl<'a> fmt::Display for Type2ChebyshevSet<'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() - ) - } -} - -#[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> 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() - ) - } -} - -trait ChebyshevSet<'a>: PartialEq + fmt::Display { - fn new( - init_epoch: Epoch, - interval_length: Duration, - rsize: usize, - num_records: usize, - record_data: &'a [f64], - ) -> Self; - - fn init_epoch(&self) -> Epoch; - fn interval_length(&self) -> Duration; - fn rsize(&self) -> usize; - fn num_records(&self) -> usize; - fn record_data(&self) -> &'a [f64]; - - fn set_init_epoch(&mut self, init_epoch: Epoch); - fn set_num_records(&mut self, num_records: usize); - fn set_record_data(&mut self, record_data: &'a [f64]); - - fn degree(&self) -> usize { - (self.rsize() - 2) / 3 - 1 +impl<'a> Type2ChebyshevSet<'a> { + pub fn degree(&self) -> usize { + (self.rsize - 2) / 3 - 1 } fn spline_idx( @@ -106,101 +51,29 @@ trait ChebyshevSet<'a>: PartialEq + fmt::Display { }); } - let window_duration_s = self.interval_length().to_seconds(); + 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())) + Ok(((ephem_start_delta_s / window_duration_s) as usize + 1).min(self.num_records)) } } -impl<'a> ChebyshevSet<'a> for Type2ChebyshevSet<'a> { - fn new( - init_epoch: Epoch, - interval_length: Duration, - rsize: usize, - num_records: usize, - record_data: &'a [f64], - ) -> Self { - Self { - init_epoch, - interval_length, - rsize, - num_records, - record_data, - } - } - - fn init_epoch(&self) -> Epoch { - self.init_epoch - } - fn interval_length(&self) -> Duration { - self.interval_length - } - fn rsize(&self) -> usize { - self.rsize - } - fn num_records(&self) -> usize { - self.num_records - } - fn record_data(&self) -> &'a [f64] { - self.record_data - } - - fn set_init_epoch(&mut self, init_epoch: Epoch) { - self.init_epoch = init_epoch; - } - fn set_num_records(&mut self, num_records: usize) { - self.num_records = num_records; - } - fn set_record_data(&mut self, record_data: &'a [f64]) { - self.record_data = record_data; - } -} -impl<'a> ChebyshevSet<'a> for Type3ChebyshevSet<'a> { - fn new( - init_epoch: Epoch, - interval_length: Duration, - rsize: usize, - num_records: usize, - record_data: &'a [f64], - ) -> Self { - Self { - init_epoch, - interval_length, - rsize, - num_records, - record_data, - } - } - - fn init_epoch(&self) -> Epoch { - self.init_epoch - } - fn interval_length(&self) -> Duration { - self.interval_length - } - fn rsize(&self) -> usize { - self.rsize - } - fn num_records(&self) -> usize { - self.num_records - } - fn record_data(&self) -> &'a [f64] { - self.record_data - } - fn set_init_epoch(&mut self, init_epoch: Epoch) { - self.init_epoch = init_epoch; - } - fn set_num_records(&mut self, num_records: usize) { - self.num_records = num_records; - } - fn set_record_data(&mut self, record_data: &'a [f64]) { - self.record_data = record_data; +impl<'a> fmt::Display for Type2ChebyshevSet<'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, T: ChebyshevSet<'a>> NAIFDataSet<'a> for T { +impl<'a> NAIFDataSet<'a> for Type2ChebyshevSet<'a> { type StateKind = (Vector3, Vector3); type RecordKind = Type2ChebyshevRecord<'a>; const DATASET_NAME: &'static str = "Chebyshev Type 2"; @@ -250,23 +123,23 @@ impl<'a, T: ChebyshevSet<'a>> NAIFDataSet<'a> for T { let rsize = slice[slice.len() - 2] as usize; let num_records = slice[slice.len() - 1] as usize; - Ok(Self::new( - start_epoch, + Ok(Self { + init_epoch: start_epoch, interval_length, rsize, num_records, - &slice[0..slice.len() - 4], - )) + 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()) + 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(), + start: n * self.rsize, + end: (n + 1) * self.rsize, + size: self.record_data.len(), })?, )) } @@ -278,7 +151,7 @@ impl<'a, T: ChebyshevSet<'a>> NAIFDataSet<'a> for T { ) -> Result<(Vector3, Vector3), InterpolationError> { let spline_idx = self.spline_idx(epoch, summary)?; - let window_duration_s = self.interval_length().to_seconds(); + 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. @@ -306,7 +179,7 @@ impl<'a, T: ChebyshevSet<'a>> NAIFDataSet<'a> for T { 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() { + for val in self.record_data { if !val.is_finite() { return Err(IntegrityError::SubNormal { dataset: Self::DATASET_NAME, @@ -333,27 +206,23 @@ impl<'a, T: ChebyshevSet<'a>> NAIFDataSet<'a> for T { let end_idx = if let Some(end) = new_end { self.spline_idx(end, summary)? } else { - self.num_records() - 1 + self.num_records - 1 }; - self.set_record_data( - &self.record_data()[start_idx * self.rsize()..(end_idx + 1) * self.rsize()], - ); - self.set_num_records((self.record_data().len() / self.rsize()) - 1); - self.set_init_epoch( - self.nth_record(0).unwrap().midpoint_epoch() - 0.5 * self.interval_length(), - ); + 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); + 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) } @@ -402,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..bdfa9ee2 --- /dev/null +++ b/anise/src/naif/daf/datatypes/chebyshev3.rs @@ -0,0 +1,375 @@ +/* + * 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, 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) / 3 - 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; + + // Now, build the X, Y, Z data from the record data. + 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(normalized_time, coeffs, radius_s, 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())?; + 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(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; + 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..], + } + } +} + +#[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::*; From 91bdc93d5183cb9aee661024ea7f1cd52c0bc6be Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Mon, 23 Sep 2024 20:54:51 -0600 Subject: [PATCH 4/4] Fix Chebyshev type 3 --- anise/src/math/interpolation/chebyshev.rs | 30 ++++++++++++++ anise/src/math/interpolation/mod.rs | 2 +- anise/src/naif/daf/datatypes/chebyshev3.rs | 29 +++++++------ anise/tests/almanac/mod.rs | 41 ++++--------------- .../validation/type03_chebyshev_jpl_de.rs | 1 - 5 files changed, 56 insertions(+), 47 deletions(-) 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/chebyshev3.rs b/anise/src/naif/daf/datatypes/chebyshev3.rs index bdfa9ee2..e0f0082c 100644 --- a/anise/src/naif/daf/datatypes/chebyshev3.rs +++ b/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( @@ -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,7 +167,7 @@ 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; } @@ -176,7 +175,7 @@ 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())?; 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..], } } } diff --git a/anise/tests/almanac/mod.rs b/anise/tests/almanac/mod.rs index cd8a01db..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}, }; @@ -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}"); } diff --git a/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs b/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs index 705d3312..34e77f6d 100644 --- a/anise/tests/ephemerides/validation/type03_chebyshev_jpl_de.rs +++ b/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]