diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index 87980330..98ba532c 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -55,9 +55,13 @@ jobs: run: sh dev-env-setup.sh && cd .. # Return to root - name: Test debug + env: + LAGRANGE_BSP: ${{ secrets.LAGRANGE_BSP }} run: cargo test --workspace --exclude anise-gui --exclude anise-py - name: Test release + env: + LAGRANGE_BSP: ${{ secrets.LAGRANGE_BSP }} run: cargo test --release --workspace --exclude anise-gui --exclude anise-py lints: @@ -85,7 +89,7 @@ jobs: steps: - name: Checkout sources uses: actions/checkout@v4 - + - name: Download data run: | wget -O data/de421.bsp http://public-data.nyxspace.com/anise/de421.bsp @@ -115,12 +119,17 @@ jobs: - name: Rust-SPICE hermite validation run: RUST_BACKTRACE=1 cargo test validate_hermite_type13_ --features spkezr_validation --release --workspace --exclude anise-gui --exclude anise-py -- --nocapture --include-ignored --test-threads 1 - + + - name: Rust-SPICE Lagrange validation + env: + LAGRANGE_BSP: ${{ secrets.LAGRANGE_BSP }} + run: RUST_BACKTRACE=1 cargo test validate_lagrange_type9_with_varying_segment_sizes --features spkezr_validation --release --workspace --exclude anise-gui --exclude anise-py -- --nocapture --include-ignored --test-threads 1 + - name: Rust-SPICE PCK validation run: RUST_BACKTRACE=1 cargo test validate_iau_rotation_to_parent --release --workspace --exclude anise-gui --exclude anise-py -- --nocapture --ignored - + - name: Rust-SPICE BPC validation - run: | + run: | RUST_BACKTRACE=1 cargo test validate_bpc_ --release --workspace --exclude anise-gui --exclude anise-py -- --nocapture --include-ignored --test-threads 1 RUST_BACKTRACE=1 cargo test de440s_translation_verif_venus2emb --release --workspace --exclude anise-gui --exclude anise-py -- --nocapture --include-ignored --test-threads 1 @@ -148,7 +157,7 @@ jobs: steps: - name: Checkout uses: actions/checkout@v4 - + - name: Download data run: | wget -O data/de421.bsp http://public-data.nyxspace.com/anise/de421.bsp @@ -173,6 +182,9 @@ jobs: uses: taiki-e/install-action@cargo-llvm-cov - name: Generate coverage report + env: + LAGRANGE_BSP: ${{ secrets.LAGRANGE_BSP }} + RUSTFLAGS: --cfg __ui_tests run: | cd anise # Prevent the workspace flag cargo llvm-cov clean --workspace @@ -183,10 +195,9 @@ jobs: 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_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 cargo llvm-cov report --lcov > ../lcov.txt - env: - RUSTFLAGS: --cfg __ui_tests - name: Upload coverage report uses: codecov/codecov-action@v3 @@ -217,4 +228,3 @@ jobs: cd anise # Jump into the package cargo login $TOKEN cargo publish - \ No newline at end of file diff --git a/anise/Cargo.toml b/anise/Cargo.toml index 0bfa07f6..db91b18c 100644 --- a/anise/Cargo.toml +++ b/anise/Cargo.toml @@ -37,6 +37,7 @@ rust-embed = { version = "8.4.0", features = [ "interpolate-folder-path", "include-exclude", ], optional = true } +regex = {version = "1.10.5" , optional = true} [dev-dependencies] rust-spice = "0.7.6" @@ -53,7 +54,7 @@ default = ["metaload"] # Enabling this flag significantly increases compilation times due to Arrow and Polars. spkezr_validation = [] python = ["pyo3", "pyo3-log"] -metaload = ["url", "reqwest/blocking", "platform-dirs"] +metaload = ["url", "reqwest/blocking", "platform-dirs", "regex"] embed_ephem = ["rust-embed"] [[bench]] diff --git a/anise/src/almanac/metaload/metafile.rs b/anise/src/almanac/metaload/metafile.rs index a111be33..1d3ab6de 100644 --- a/anise/src/almanac/metaload/metafile.rs +++ b/anise/src/almanac/metaload/metafile.rs @@ -10,9 +10,10 @@ use log::{debug, info}; use platform_dirs::AppDirs; - +use regex::Regex; use serde_derive::{Deserialize, Serialize}; use serde_dhall::StaticType; +use std::env; use std::fs::{create_dir_all, remove_file, File}; use std::io::Write; use std::path::Path; @@ -32,6 +33,11 @@ use crate::prelude::InputOutputError; use super::MetaAlmanacError; +/// MetaFile allows downloading a remote file from a URL (http, https only), and interpolation of paths in environment variable using the Dhall syntax `env:MY_ENV_VAR`. +/// +/// The data is stored in the user's local temp directory (i.e. `~/.local/share/nyx-space/anise/` on Linux and `AppData/Local/nyx-space/anise/` on Windows). +/// Prior to loading a remote resource, if the local resource exists, its CRC32 will be computed: if it matches the CRC32 of this instance of MetaFile, +/// then the file will not be downloaded a second time. #[cfg_attr(feature = "python", pyclass)] #[cfg_attr(feature = "python", pyo3(module = "anise"))] #[cfg_attr(feature = "python", pyo3(get_all, set_all))] @@ -53,6 +59,8 @@ impl MetaFile { } pub(crate) fn _process(&mut self) -> Result<(), MetaAlmanacError> { + // First, parse environment variables if any. + self.uri = replace_env_vars(&self.uri); match Url::parse(&self.uri) { Err(e) => { debug!("parsing {} caused {e} -- assuming local path", self.uri); @@ -274,6 +282,15 @@ impl MetaFile { } } +fn replace_env_vars(input: &str) -> String { + let re = Regex::new(r"env:([A-Z_][A-Z0-9_]*)").unwrap(); + re.replace_all(input, |caps: ®ex::Captures| { + let var_name = &caps[1]; + env::var(var_name).unwrap_or_else(|_| format!("env:{}", var_name)) + }) + .to_string() +} + #[cfg(test)] mod ut_metafile { use super::MetaFile; @@ -308,4 +325,25 @@ mod ut_metafile { assert!(unix_rel_path._process().is_ok()); assert_eq!(unix_rel_path.uri, "../Users/me/meta.dhall".to_string()); } + + #[test] + fn test_metafile_regex() { + use std::env; + let mut user_path = MetaFile { + uri: "env:USER/.cargo/env".to_string(), + crc32: None, + }; + user_path._process().unwrap(); + assert_eq!(user_path.uri, env::var("USER").unwrap() + "/.cargo/env"); + + let mut unknown_path = MetaFile { + uri: "env:BLAH_BLAH_NO_EXIST/.cargo/env".to_string(), + crc32: None, + }; + unknown_path._process().unwrap(); + assert_eq!( + unknown_path.uri, + "env:BLAH_BLAH_NO_EXIST/.cargo/env".to_string() + ); + } } diff --git a/anise/src/math/interpolation/hermite.rs b/anise/src/math/interpolation/hermite.rs index 4da845c3..7bcd3a28 100644 --- a/anise/src/math/interpolation/hermite.rs +++ b/anise/src/math/interpolation/hermite.rs @@ -16,7 +16,6 @@ 2. The relevant comments (including authors) from hrmint are kept. 3. The tests are not part of the original SPICE code. 4. The transliteration in itself justifies the change of license from unrestricted to MPL. - */ /* hrmint.f -- translated by f2c (version 19980913). You must link the resulting object file with the libraries: @@ -25,37 +24,37 @@ /* $ Restrictions */ -/* None. */ +/* None. */ /* $ Literature_References */ -/* [1] "Numerical Recipes---The Art of Scientific Computing" by */ -/* William H. Press, Brian P. Flannery, Saul A. Teukolsky, */ -/* William T. Vetterling (see sections 3.0 and 3.1). */ +/* [1] "Numerical Recipes---The Art of Scientific Computing" by */ +/* William H. Press, Brian P. Flannery, Saul A. Teukolsky, */ +/* William T. Vetterling (see sections 3.0 and 3.1). */ -/* [2] "Elementary Numerical Analysis---An Algorithmic Approach" */ -/* by S. D. Conte and Carl de Boor. See p. 64. */ +/* [2] "Elementary Numerical Analysis---An Algorithmic Approach" */ +/* by S. D. Conte and Carl de Boor. See p. 64. */ /* $ Author_and_Institution */ -/* N.J. Bachman (JPL) */ +/* N.J. Bachman (JPL) */ /* $ Version */ /* - SPICELIB Version 1.2.1, 28-JAN-2014 (NJB) */ -/* Fixed a few comment typos. */ +/* Fixed a few comment typos. */ /* - SPICELIB Version 1.2.0, 01-FEB-2002 (NJB) (EDW) */ -/* Bug fix: declarations of local variables XI and XIJ */ -/* were changed from DOUBLE PRECISION to INTEGER. */ -/* Note: bug had no effect on behavior of this routine. */ +/* Bug fix: declarations of local variables XI and XIJ */ +/* were changed from DOUBLE PRECISION to INTEGER. */ +/* Note: bug had no effect on behavior of this routine. */ /* - SPICELIB Version 1.1.0, 28-DEC-2001 (NJB) */ -/* Blanks following final newline were truncated to */ -/* suppress compilation warnings on the SGI-N32 platform. */ +/* Blanks following final newline were truncated to */ +/* suppress compilation warnings on the SGI-N32 platform. */ /* - SPICELIB Version 1.0.0, 01-MAR-2000 (NJB) */ @@ -93,27 +92,27 @@ pub fn hermite_eval( // At this point, we know that the lengths of items is correct, so we can directly address them without worry for overflowing the array. - let work: &mut [f64] = &mut [0.0; 256]; + let work: &mut [f64] = &mut [0.0; 8 * MAX_SAMPLES]; let n: usize = xs.len(); - /* Copy the input array into WORK. After this, the first column */ - /* of WORK represents the first column of our triangular */ - /* interpolation table. */ + /* Copy the input array into WORK. After this, the first column */ + /* of WORK represents the first column of our triangular */ + /* interpolation table. */ for i in 0..n { work[2 * i] = ys[i]; work[2 * i + 1] = ydots[i]; } - /* Compute the second column of the interpolation table: this */ - /* consists of the N-1 values obtained by evaluating the */ - /* first-degree interpolants at X. We'll also evaluate the */ - /* derivatives of these interpolants at X and save the results in */ - /* the second column of WORK. Because the derivative computations */ - /* depend on the function computations from the previous column in */ - /* the interpolation table, and because the function interpolation */ - /* overwrites the previous column of interpolated function values, */ - /* we must evaluate the derivatives first. */ + /* Compute the second column of the interpolation table: this */ + /* consists of the N-1 values obtained by evaluating the */ + /* first-degree interpolants at X. We'll also evaluate the */ + /* derivatives of these interpolants at X and save the results in */ + /* the second column of WORK. Because the derivative computations */ + /* depend on the function computations from the previous column in */ + /* the interpolation table, and because the function interpolation */ + /* overwrites the previous column of interpolated function values, */ + /* we must evaluate the derivatives first. */ for i in 1..=n - 1 { let c1 = xs[i] - x_eval; @@ -128,51 +127,50 @@ pub fn hermite_eval( }); } - /* The second column of WORK contains interpolated derivative */ - /* values. */ + /* The second column of WORK contains interpolated derivative */ + /* values. */ - /* The odd-indexed interpolated derivatives are simply the input */ - /* derivatives. */ + /* The odd-indexed interpolated derivatives are simply the input */ + /* derivatives. */ let prev = 2 * i - 1; let curr = 2 * i; work[prev + 2 * n - 1] = work[prev]; - /* The even-indexed interpolated derivatives are the slopes of */ - /* the linear interpolating polynomials for adjacent input */ - /* abscissa/ordinate pairs. */ + /* The even-indexed interpolated derivatives are the slopes of */ + /* the linear interpolating polynomials for adjacent input */ + /* abscissa/ordinate pairs. */ work[prev + 2 * n] = (work[curr] - work[prev - 1]) / denom; - /* The first column of WORK contains interpolated function values. */ - /* The odd-indexed entries are the linear Taylor polynomials, */ - /* for each input abscissa value, evaluated at X. */ + /* The first column of WORK contains interpolated function values. */ + /* The odd-indexed entries are the linear Taylor polynomials, */ + /* for each input abscissa value, evaluated at X. */ let temp = work[prev] * (x_eval - xs[i - 1]) + work[prev - 1]; work[prev] = (c1 * work[prev - 1] + c2 * work[curr]) / denom; work[prev - 1] = temp; } - /* The last column entries were not computed by the preceding loop; */ - /* compute them now. */ + /* The last column entries were not computed by the preceding loop; */ + /* compute them now. */ work[4 * n - 2] = work[(2 * n) - 1]; work[2 * (n - 1)] += work[(2 * n) - 1] * (x_eval - xs[n - 1]); - /* Compute columns 3 through 2*N of the table. */ + /* Compute columns 3 through 2*N of the table. */ for j in 2..=(2 * n) - 1 { for i in 1..=(2 * n) - j { - /* In the theoretical construction of the interpolation table, - */ - /* there are 2*N abscissa values, since each input abcissa */ - /* value occurs with multiplicity two. In this theoretical */ - /* construction, the Jth column of the interpolation table */ - /* contains results of evaluating interpolants that span J+1 */ - /* consecutive abscissa values. The indices XI and XIJ below */ - /* are used to pick the correct abscissa values out of the */ - /* physical XVALS array, in which the abscissa values are not */ - /* repeated. */ + /* In the theoretical construction of the interpolation table, */ + /* there are 2*N abscissa values, since each input abcissa */ + /* value occurs with multiplicity two. In this theoretical */ + /* construction, the Jth column of the interpolation table */ + /* contains results of evaluating interpolants that span J+1 */ + /* consecutive abscissa values. The indices XI and XIJ below */ + /* are used to pick the correct abscissa values out of the */ + /* physical XVALS array, in which the abscissa values are not */ + /* repeated. */ let xi = (i + 1) / 2; let xij = (i + j + 1) / 2; @@ -181,36 +179,35 @@ pub fn hermite_eval( let denom = xs[xij - 1] - xs[xi - 1]; if denom.abs() < f64::EPSILON { return Err(InterpolationError::InterpMath { - source:MathError::DivisionByZero { - action: - "hermite data contains likely duplicate abcissa, remove duplicate states", - }}); + source: MathError::DivisionByZero { + action: "hermite data contains duplicate states", + }, + }); } - /* Compute the interpolated derivative at X for the Ith */ - /* interpolant. This is the derivative with respect to X of */ - /* the expression for the interpolated function value, which */ - /* is the second expression below. This derivative computation - */ - /* is done first because it relies on the interpolated */ - /* function values from the previous column of the */ - /* interpolation table. */ + /* Compute the interpolated derivative at X for the Ith */ + /* interpolant. This is the derivative with respect to X of */ + /* the expression for the interpolated function value, which */ + /* is the second expression below. This derivative computation */ + /* is done first because it relies on the interpolated */ + /* function values from the previous column of the */ + /* interpolation table. */ - /* The derivative expression here corresponds to equation */ - /* 2.35 on page 64 in reference [2]. */ + /* The derivative expression here corresponds to equation */ + /* 2.35 on page 64 in reference [2]. */ work[i + 2 * n - 1] = (c1 * work[i + 2 * n - 1] + c2 * work[i + 2 * n] + (work[i] - work[i - 1])) / denom; - /* Compute the interpolated function value at X for the Ith */ - /* interpolant. */ + /* Compute the interpolated function value at X for the Ith */ + /* interpolant. */ work[i - 1] = (c1 * work[i - 1] + c2 * work[i]) / denom; } } - /* Our interpolated function value is sitting in WORK(1,1) at this */ - /* point. The interpolated derivative is located in WORK(1,2). */ + /* Our interpolated function value is sitting in WORK(1,1) at this */ + /* point. The interpolated derivative is located in WORK(1,2). */ let f = work[0]; let df = work[2 * n]; diff --git a/anise/src/math/interpolation/lagrange.rs b/anise/src/math/interpolation/lagrange.rs new file mode 100644 index 00000000..c2524e75 --- /dev/null +++ b/anise/src/math/interpolation/lagrange.rs @@ -0,0 +1,102 @@ +/* + * 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 crate::errors::MathError; + +use super::{InterpolationError, MAX_SAMPLES}; + +pub fn lagrange_eval( + xs: &[f64], + ys: &[f64], + x_eval: f64, +) -> Result<(f64, f64), InterpolationError> { + if xs.len() != ys.len() { + return Err(InterpolationError::CorruptedData { + what: "lengths of abscissas (xs), ordinates (ys), and first derivatives (ydots) differ", + }); + } else if xs.is_empty() { + return Err(InterpolationError::CorruptedData { + what: "list of abscissas (xs) is empty", + }); + } + + // At this point, we know that the lengths of items is correct, so we can directly address them without worry for overflowing the array. + + let work: &mut [f64] = &mut [0.0; MAX_SAMPLES]; + let dwork: &mut [f64] = &mut [0.0; MAX_SAMPLES]; + for (ii, y) in ys.iter().enumerate() { + work[ii] = *y; + } + + let n = xs.len(); + + for j in 1..n { + for i in 0..(n - j) { + let xi = xs[i]; + let xij = xs[i + j]; + + let denom = xi - xij; + if denom.abs() < f64::EPSILON { + return Err(InterpolationError::InterpMath { + source: MathError::DivisionByZero { + action: "lagrange data contains duplicate states", + }, + }); + } + + let work_i = work[i]; + let work_ip1 = work[i + 1]; + + work[i] = ((x_eval - xij) * work_i + (xi - x_eval) * work_ip1) / denom; + + let deriv = (work_i - work_ip1) / denom; + dwork[i] = ((x_eval - xij) * dwork[i] + (xi - x_eval) * dwork[i + 1]) / denom + deriv; + } + } + + let f = work[0]; + let df = dwork[0]; + Ok((f, df)) +} + +#[test] +fn lagrange_spice_docs_example() { + let ts = [-1.0, 0.0, 3.0, 5.0]; + let yvals = [-2.0, -7.0, -8.0, 26.0]; + let dyvals = [ + -4.6333333333333334, + -4.9833333333333333, + 7.76666666666666666, + 27.7666666666666662, + ]; + + // Check that we can interpolate the values exactly. + for (i, t) in ts.iter().enumerate() { + let (eval, deval) = lagrange_eval(&ts, &yvals, *t).unwrap(); + let eval_err = (eval - yvals[i]).abs(); + assert!(eval_err < f64::EPSILON, "f(x) error is {eval_err:e}"); + + let deval_err = (deval - dyvals[i]).abs(); + assert!( + deval_err < f64::EPSILON, + "#{i}: f'(x) error is {deval_err:e}" + ); + } + + // Check the interpolation from the SPICE documentation + let (x, dx) = lagrange_eval(&ts, &yvals, 2.0).unwrap(); + + // WARNING: The documentation data is wrong! Evaluation at 2.0 returns -12.2999... in spiceypy. + let expected_x = -12.299999999999999; + let expected_dx = 1.2166666666666666; + dbg!(x, dx); + assert!((x - expected_x).abs() < f64::EPSILON, "X error"); + assert!((dx - expected_dx).abs() < f64::EPSILON, "dX error"); +} diff --git a/anise/src/math/interpolation/mod.rs b/anise/src/math/interpolation/mod.rs index 2fa871ca..61d7977b 100644 --- a/anise/src/math/interpolation/mod.rs +++ b/anise/src/math/interpolation/mod.rs @@ -10,10 +10,12 @@ mod chebyshev; mod hermite; +mod lagrange; pub use chebyshev::chebyshev_eval; pub use hermite::hermite_eval; use hifitime::Epoch; +pub use lagrange::lagrange_eval; use snafu::Snafu; use crate::errors::{DecodingError, MathError}; @@ -50,4 +52,8 @@ pub enum InterpolationError { kind: &'static str, op: &'static str, }, + #[snafu(display( + "{dataset} is not yet supported -- https://github.com/nyx-space/anise/issues/{issue}" + ))] + UnimplementedType { issue: u32, dataset: &'static str }, } diff --git a/anise/src/naif/daf/datatypes/hermite.rs b/anise/src/naif/daf/datatypes/hermite.rs index 450c59dd..2d451e34 100644 --- a/anise/src/naif/daf/datatypes/hermite.rs +++ b/anise/src/naif/daf/datatypes/hermite.rs @@ -115,7 +115,10 @@ impl<'a> NAIFDataSet<'a> for HermiteSetType12<'a> { _epoch: Epoch, _: &S, ) -> Result { - todo!("https://github.com/anise-toolkit/anise.rs/issues/14") + Err(InterpolationError::UnimplementedType { + dataset: Self::DATASET_NAME, + issue: 14, + }) } fn check_integrity(&self) -> Result<(), IntegrityError> { diff --git a/anise/src/naif/daf/datatypes/lagrange.rs b/anise/src/naif/daf/datatypes/lagrange.rs index f1110535..c57554fd 100644 --- a/anise/src/naif/daf/datatypes/lagrange.rs +++ b/anise/src/naif/daf/datatypes/lagrange.rs @@ -10,11 +10,15 @@ use core::fmt; use hifitime::{Duration, Epoch, TimeUnits}; -use snafu::ensure; +use snafu::{ensure, ResultExt}; use crate::{ errors::{DecodingError, IntegrityError, TooFewDoublesSnafu}, - math::{cartesian::CartesianState, interpolation::InterpolationError, Vector3}, + math::{ + cartesian::CartesianState, + interpolation::{lagrange_eval, InterpDecodingSnafu, InterpolationError, MAX_SAMPLES}, + Vector3, + }, naif::daf::{NAIFDataRecord, NAIFDataSet, NAIFRecord, NAIFSummaryRecord}, DBL_SIZE, }; @@ -112,7 +116,10 @@ impl<'a> NAIFDataSet<'a> for LagrangeSetType8<'a> { _epoch: Epoch, _: &S, ) -> Result { - todo!("https://github.com/anise-toolkit/anise.rs/issues/12") + Err(InterpolationError::UnimplementedType { + dataset: Self::DATASET_NAME, + issue: 12, + }) } fn check_integrity(&self) -> Result<(), IntegrityError> { @@ -202,10 +209,115 @@ impl<'a> NAIFDataSet<'a> for LagrangeSetType9<'a> { fn evaluate( &self, - _epoch: Epoch, + epoch: Epoch, _: &S, ) -> Result { - todo!("https://github.com/anise-toolkit/anise.rs/issues/13") + // Start by doing a binary search on the epoch registry to limit the search space in the total number of epochs. + // TODO: use the epoch registry to reduce the search space + // Check that we even have interpolation data for that time + if epoch.to_et_seconds() + 1e-9 < self.epoch_data[0] + || epoch.to_et_seconds() - 1e-9 > *self.epoch_data.last().unwrap() + { + return Err(InterpolationError::NoInterpolationData { + req: epoch, + start: Epoch::from_et_seconds(self.epoch_data[0]), + end: Epoch::from_et_seconds(*self.epoch_data.last().unwrap()), + }); + } + // Now, perform a binary search on the epochs themselves. + match self.epoch_data.binary_search_by(|epoch_et| { + epoch_et + .partial_cmp(&epoch.to_et_seconds()) + .expect("epochs in Hermite data is now NaN or infinite but was not before") + }) { + Ok(idx) => { + // Oh wow, this state actually exists, no interpolation needed! + Ok(self + .nth_record(idx) + .context(InterpDecodingSnafu)? + .to_pos_vel()) + } + Err(idx) => { + // We didn't find it, so let's build an interpolation here. + let group_size = self.degree + 1; + let num_left = group_size / 2; + + // Ensure that we aren't fetching out of the window + let mut first_idx = idx.saturating_sub(num_left); + let last_idx = self.num_records.min(first_idx + group_size); + + // Check that we have enough samples + if last_idx == self.num_records { + first_idx = last_idx - 2 * num_left; + } + + // Statically allocated arrays of the maximum number of samples + let mut epochs = [0.0; MAX_SAMPLES]; + let mut xs = [0.0; MAX_SAMPLES]; + let mut ys = [0.0; MAX_SAMPLES]; + let mut zs = [0.0; MAX_SAMPLES]; + let mut vxs = [0.0; MAX_SAMPLES]; + let mut vys = [0.0; MAX_SAMPLES]; + let mut vzs = [0.0; MAX_SAMPLES]; + + for (cno, idx) in (first_idx..last_idx).enumerate() { + let record = self.nth_record(idx).context(InterpDecodingSnafu)?; + xs[cno] = record.x_km; + ys[cno] = record.y_km; + zs[cno] = record.z_km; + vxs[cno] = record.vx_km_s; + vys[cno] = record.vy_km_s; + vzs[cno] = record.vz_km_s; + epochs[cno] = self.epoch_data[idx]; + } + + // TODO: Build a container that uses the underlying data and provides an index into it. + + // Build the interpolation polynomials making sure to limit the slices to exactly the number of items we actually used + // The other ones are zeros, which would cause the interpolation function to fail. + let (x_km, _) = lagrange_eval( + &epochs[..group_size], + &xs[..group_size], + epoch.to_et_seconds(), + )?; + + let (y_km, _) = lagrange_eval( + &epochs[..group_size], + &ys[..group_size], + epoch.to_et_seconds(), + )?; + + let (z_km, _) = lagrange_eval( + &epochs[..group_size], + &zs[..group_size], + epoch.to_et_seconds(), + )?; + + let (vx_km_s, _) = lagrange_eval( + &epochs[..group_size], + &vxs[..group_size], + epoch.to_et_seconds(), + )?; + + let (vy_km_s, _) = lagrange_eval( + &epochs[..group_size], + &vys[..group_size], + epoch.to_et_seconds(), + )?; + + let (vz_km_s, _) = lagrange_eval( + &epochs[..group_size], + &vzs[..group_size], + epoch.to_et_seconds(), + )?; + + // And build the result + let pos_km = Vector3::new(x_km, y_km, z_km); + let vel_km_s = Vector3::new(vx_km_s, vy_km_s, vz_km_s); + + Ok((pos_km, vel_km_s)) + } + } } fn check_integrity(&self) -> Result<(), IntegrityError> { diff --git a/anise/tests/ephemerides/translation.rs b/anise/tests/ephemerides/translation.rs index afe75c88..48190454 100644 --- a/anise/tests/ephemerides/translation.rs +++ b/anise/tests/ephemerides/translation.rs @@ -299,19 +299,6 @@ fn de438s_translation_verif_emb2moon() { fn spk_hermite_type13_verif() { let _ = pretty_env_logger::try_init().is_err(); - // "Load" the file via a memory map (avoids allocations) - // let path = "../data/de440s.bsp"; - // let buf = file2heap!(path).unwrap(); - // let spk = SPK::parse(buf).unwrap(); - - // let buf = file2heap!("../data/gmat-hermite.bsp").unwrap(); - // let spacecraft = SPK::parse(buf).unwrap(); - - // let ctx = Almanac::from_spk(spk) - // .unwrap() - // .with_spk(spacecraft) - // .unwrap(); - let ctx = Almanac::default() .load("../data/de440s.bsp") .and_then(|ctx| ctx.load("../data/gmat-hermite.bsp")) @@ -633,3 +620,89 @@ fn de440s_translation_verif_aberrations() { ); } } + +#[cfg(feature = "metaload")] +#[test] +fn type9_lagrange_query() { + use anise::almanac::metaload::MetaFile; + use anise::constants::frames::EARTH_J2000; + use anise::prelude::Frame; + + let lagrange_meta = MetaFile { + uri: "http://public-data.nyxspace.com/anise/ci/env:LAGRANGE_BSP".to_string(), + crc32: None, + }; + + let almanac = Almanac::default() + .load_from_metafile(lagrange_meta) + .unwrap(); + + let obj_id = -10; + let obj_frame = Frame::from_ephem_j2000(obj_id); + + let (start, end) = almanac.spk_domain(obj_id).unwrap(); + + // Query near the start, but somewhere where the state is not exactly defined + let state = almanac + .translate( + obj_frame, + EARTH_J2000, + start + 1.159_f64.seconds(), + Aberration::NONE, + ) + .unwrap(); + + let expected_pos_km = Vector3::new(-7338.44373643, 3159.7629953, 760.74472775); + let expected_vel_km_s = Vector3::new(-7.21781188, -5.26834555, -4.12581558); + + dbg!(state.radius_km - expected_pos_km); + dbg!(state.velocity_km_s - expected_vel_km_s); + + assert!( + relative_eq!(state.radius_km, expected_pos_km, epsilon = 5e-6), + "got {} but want {} => err = {:.3e} km", + state.radius_km, + expected_pos_km, + (state.radius_km - expected_pos_km).norm() + ); + + assert!( + relative_eq!(state.velocity_km_s, expected_vel_km_s, epsilon = 5e-9), + "got {} but want {} => err = {:.3e} km/s", + state.velocity_km_s, + expected_vel_km_s, + (state.velocity_km_s - expected_vel_km_s).norm() + ); + + // Query near the end, but not in the registry either + let state = almanac + .translate( + obj_frame, + EARTH_J2000, + end - 1159.56_f64.seconds(), + Aberration::NONE, + ) + .unwrap(); + + let expected_pos_km = Vector3::new(10106.15561792, -166810.06321791, -95547.93140678); + let expected_vel_km_s = Vector3::new(0.43375448, -1.07193959, -0.55979951); + + dbg!(state.radius_km - expected_pos_km); + dbg!(state.velocity_km_s - expected_vel_km_s); + + assert!( + relative_eq!(state.radius_km, expected_pos_km, epsilon = 5e-6), + "got {} but want {} => err = {:.3e} km", + state.radius_km, + expected_pos_km, + (state.radius_km - expected_pos_km).norm() + ); + + assert!( + relative_eq!(state.velocity_km_s, expected_vel_km_s, epsilon = 5e-9), + "got {} but want {} => err = {:.3e} km/s", + state.velocity_km_s, + expected_vel_km_s, + (state.velocity_km_s - expected_vel_km_s).norm() + ); +} diff --git a/anise/tests/ephemerides/validation/compare.rs b/anise/tests/ephemerides/validation/compare.rs index 00afeba4..558be03e 100644 --- a/anise/tests/ephemerides/validation/compare.rs +++ b/anise/tests/ephemerides/validation/compare.rs @@ -68,7 +68,6 @@ impl EphemValData { /// An ephemeris comparison tool that writes the differences between ephemerides to a Parquet file. pub struct CompareEphem { pub input_file_names: Vec, - pub output_file_name: String, pub num_queries_per_pair: usize, pub dry_run: bool, pub aberration: Option, @@ -109,7 +108,6 @@ impl CompareEphem { let writer = ArrowWriter::try_new(file, Arc::new(schema), Some(props)).unwrap(); Self { - output_file_name, input_file_names, num_queries_per_pair, aberration, @@ -313,8 +311,8 @@ impl CompareEphem { 0 => (data.spice_val_x_km, data.anise_val_x_km), 1 => (data.spice_val_y_km, data.anise_val_y_km), 2 => (data.spice_val_z_km, data.anise_val_z_km), - 3 => (data.spice_val_vy_km_s, data.anise_val_vy_km_s), - 4 => (data.spice_val_vz_km_s, data.anise_val_vz_km_s), + 3 => (data.spice_val_vx_km_s, data.anise_val_vx_km_s), + 4 => (data.spice_val_vy_km_s, data.anise_val_vy_km_s), 5 => (data.spice_val_vz_km_s, data.anise_val_vz_km_s), _ => unreachable!(), }; diff --git a/anise/tests/ephemerides/validation/mod.rs b/anise/tests/ephemerides/validation/mod.rs index b99437b6..1e362fec 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 type09_lagrange; mod type13_hermite; mod compare; diff --git a/anise/tests/ephemerides/validation/type09_lagrange.rs b/anise/tests/ephemerides/validation/type09_lagrange.rs new file mode 100644 index 00000000..e843b2df --- /dev/null +++ b/anise/tests/ephemerides/validation/type09_lagrange.rs @@ -0,0 +1,38 @@ +/* + * 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::almanac::metaload::MetaFile; + +#[ignore = "Requires Rust SPICE -- must be executed serially"] +#[test] +fn validate_lagrange_type9_with_varying_segment_sizes() { + let mut lagrange_meta = MetaFile { + uri: "http://public-data.nyxspace.com/anise/ci/env:LAGRANGE_BSP".to_string(), + crc32: None, + }; + lagrange_meta.process().unwrap(); + + let file_name = "spk-type9-validation-variable-seg-size".to_string(); + let comparator = CompareEphem::new(vec![lagrange_meta.uri], file_name.clone(), 10_000, None); + + let err_count = comparator.run(); + + assert_eq!(err_count, 0, "None of the queries should fail!"); + + let validator = Validation { + file_name, + max_q75_err: 5e-9, + max_q99_err: 2e-7, + max_abs_err: 0.05, + }; + + validator.validate(); +} diff --git a/anise/tests/ephemerides/validation/type13_hermite.rs b/anise/tests/ephemerides/validation/type13_hermite.rs index a9304847..8d089a41 100644 --- a/anise/tests/ephemerides/validation/type13_hermite.rs +++ b/anise/tests/ephemerides/validation/type13_hermite.rs @@ -36,7 +36,7 @@ fn validate_hermite_type13_from_gmat() { #[ignore = "Requires Rust SPICE -- must be executed serially"] #[test] fn validate_hermite_type13_with_varying_segment_sizes() { - // ISSUE: This file is corrupt. I messed up the rewrite. + // ISSUE: This file is corrupt, cf. https://github.com/nyx-space/anise/issues/262 let file_name = "spk-type13-validation-variable-seg-size".to_string(); let comparator = CompareEphem::new( vec!["../data/variable-seg-size-hermite.bsp".to_string()],