diff --git a/.gitattributes b/.gitattributes index 292d92e5..9ef06974 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,7 +1,3 @@ data/*.bsp filter=lfs diff=lfs merge=lfs -text data/*.bpc filter=lfs diff=lfs merge=lfs -text -data/de421.anise filter=lfs diff=lfs merge=lfs -text -data/de430.anise filter=lfs diff=lfs merge=lfs -text -data/de438s.anise filter=lfs diff=lfs merge=lfs -text -data/de440.anise filter=lfs diff=lfs merge=lfs -text data/*.pca filter=lfs diff=lfs merge=lfs -text \ No newline at end of file diff --git a/.github/workflows/benchmarks.yml b/.github/workflows/benchmarks.yml index 894b354c..0731e735 100644 --- a/.github/workflows/benchmarks.yml +++ b/.github/workflows/benchmarks.yml @@ -16,8 +16,17 @@ jobs: steps: - name: Checkout sources uses: actions/checkout@v3 - with: - lfs: true + + - name: Download data + run: | + wget -O data/de421.bsp http://public-data.nyxspace.com/anise/de421.bsp + 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/pck08.pca http://public-data.nyxspace.com/anise/pck08.pca + wget -O data/gmat-hermite.bsp http://public-data.nyxspace.com/anise/ci/gmat-hermite.bsp + wget -O data/variable-seg-size-hermite.bsp http://public-data.nyxspace.com/anise/ci/variable-seg-size-hermite.bsp + wget -O data/earth_latest_high_prec.bpc http://public-data.nyxspace.com/anise/ci/earth_latest_high_prec-2023-09-08.bpc - name: Install stable toolchain uses: actions-rs/toolchain@v1 diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml index 0ad072ad..4812c5b1 100644 --- a/.github/workflows/python.yml +++ b/.github/workflows/python.yml @@ -26,8 +26,14 @@ jobs: target: [x86_64, x86, aarch64, armv7, s390x, ppc64le] steps: - uses: actions/checkout@v3 - with: - lfs: true + + - name: Download data + run: | + wget -O data/de421.bsp http://public-data.nyxspace.com/anise/de421.bsp + 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/pck08.pca http://public-data.nyxspace.com/anise/pck08.pca - uses: actions/setup-python@v4 with: @@ -172,8 +178,14 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v3 - with: - lfs: true + + - name: Download data + run: | + wget -O data/de421.bsp http://public-data.nyxspace.com/anise/de421.bsp + 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/pck08.pca http://public-data.nyxspace.com/anise/pck08.pca - name: Build sdist uses: PyO3/maturin-action@v1 diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index cbe620cf..8cbf0f9b 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -40,8 +40,17 @@ jobs: steps: - name: Checkout sources uses: actions/checkout@v2 - with: - lfs: true + + - name: Download data + run: | + wget -O data/de421.bsp http://public-data.nyxspace.com/anise/de421.bsp + 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/pck08.pca http://public-data.nyxspace.com/anise/pck08.pca + wget -O data/gmat-hermite.bsp http://public-data.nyxspace.com/anise/ci/gmat-hermite.bsp + wget -O data/variable-seg-size-hermite.bsp http://public-data.nyxspace.com/anise/ci/variable-seg-size-hermite.bsp + wget -O data/earth_latest_high_prec.bpc http://public-data.nyxspace.com/anise/ci/earth_latest_high_prec-2023-09-08.bpc - name: Install stable toolchain uses: actions-rs/toolchain@v1 @@ -93,8 +102,17 @@ jobs: steps: - name: Checkout sources uses: actions/checkout@v2 - with: - lfs: true + + - name: Download data + run: | + wget -O data/de421.bsp http://public-data.nyxspace.com/anise/de421.bsp + 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/pck08.pca http://public-data.nyxspace.com/anise/pck08.pca + wget -O data/gmat-hermite.bsp http://public-data.nyxspace.com/anise/ci/gmat-hermite.bsp + wget -O data/variable-seg-size-hermite.bsp http://public-data.nyxspace.com/anise/ci/variable-seg-size-hermite.bsp + wget -O data/earth_latest_high_prec.bpc http://public-data.nyxspace.com/anise/ci/earth_latest_high_prec-2023-09-08.bpc - name: Install stable toolchain uses: actions-rs/toolchain@v1 @@ -150,8 +168,17 @@ jobs: steps: - name: Checkout uses: actions/checkout@v3 - with: - lfs: true + + - name: Download data + run: | + wget -O data/de421.bsp http://public-data.nyxspace.com/anise/de421.bsp + 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/pck08.pca http://public-data.nyxspace.com/anise/pck08.pca + wget -O data/gmat-hermite.bsp http://public-data.nyxspace.com/anise/ci/gmat-hermite.bsp + wget -O data/variable-seg-size-hermite.bsp http://public-data.nyxspace.com/anise/ci/variable-seg-size-hermite.bsp + wget -O data/earth_latest_high_prec.bpc http://public-data.nyxspace.com/anise/ci/earth_latest_high_prec-2023-09-08.bpc - name: Install stable toolchain uses: actions-rs/toolchain@v1 @@ -169,13 +196,15 @@ jobs: - name: Generate coverage report run: | + cd anise # Prevent the workspace flag cargo llvm-cov clean --workspace - cargo llvm-cov test --workspace --exclude anise-gui --exclude anise-py --no-report -- --test-threads=1 - cargo llvm-cov test --workspace --exclude anise-gui --exclude anise-py --no-report --tests -- compile_fail - cargo llvm-cov test --workspace --exclude anise-gui --exclude anise-py --no-report validate_iau_rotation_to_parent -- --nocapture --ignored - cargo llvm-cov test --workspace --exclude anise-gui --exclude anise-py --no-report validate_bpc_to_iau_rotations -- --nocapture --ignored - cargo llvm-cov test --workspace --exclude anise-gui --exclude anise-py --no-report validate_jplde_de440s --features spkezr_validation -- --nocapture --ignored - cargo llvm-cov test --workspace --exclude anise-gui --exclude anise-py --no-report validate_hermite_type13_from_gmat --features spkezr_validation -- --nocapture --ignored + cargo llvm-cov test --no-report -- --test-threads=1 + cargo llvm-cov test --no-report --tests -- compile_fail + cargo llvm-cov test --no-report validate_iau_rotation_to_parent -- --nocapture --ignored + 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_hermite_type13_from_gmat --features spkezr_validation -- --nocapture --ignored cargo llvm-cov report --lcov > lcov.txt env: RUSTFLAGS: --cfg __ui_tests diff --git a/Cargo.toml b/Cargo.toml index a410226f..be830bc5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,7 +3,7 @@ resolver = "2" members = ["anise", "anise-cli", "anise-gui", "anise-py"] [workspace.package] -version = "0.1.0" +version = "0.2.0" edition = "2021" authors = ["Christopher Rabotin "] description = "ANISE provides a toolkit and files for Attitude, Navigation, Instrument, Spacecraft, and Ephemeris data. It's a modern replacement of NAIF SPICE file." @@ -26,7 +26,7 @@ exclude = [ ] [workspace.dependencies] -hifitime = { version = "3.8.6", default-features = true } +hifitime = { version = "3.9.0", default-features = true } memmap2 = "=0.9.3" crc32fast = "=1.3.2" der = { version = "0.7.8", features = ["derive", "alloc", "real"] } @@ -45,7 +45,7 @@ rstest = "0.18.2" pyo3 = { version = "0.20.0", features = ["multiple-pymethods"] } pyo3-log = "0.9.0" -anise = { version = "0.1.0", path = "anise", default-features = false } +anise = { version = "0.2.0", path = "anise", default-features = false } [profile.bench] debug = true diff --git a/README.md b/README.md index 00200ead..283e607a 100644 --- a/README.md +++ b/README.md @@ -99,7 +99,7 @@ let orig_state = Orbit::keplerian( // Transform that orbit into another frame. let state_itrf93 = almanac - .transform_to(orig_state, EARTH_ITRF93, Aberration::NotSet) + .transform_to(orig_state, EARTH_ITRF93, None) .unwrap(); // The `:x` prints this orbit's Keplerian elements @@ -109,7 +109,7 @@ println!("{state_itrf93:X}"); // Convert back let from_state_itrf93_to_eme2k = almanac - .transform_to(state_itrf93, EARTH_J2000, Aberration::NotSet) + .transform_to(state_itrf93, EARTH_J2000, Aberration::NONE) .unwrap(); println!("{from_state_itrf93_to_eme2k}"); @@ -181,11 +181,11 @@ let ctx = Almanac::from_spk(spk).unwrap(); let epoch = Epoch::from_str("2020-11-15 12:34:56.789 TDB").unwrap(); let state = ctx - .translate_from_to( - VENUS_J2000, - EARTH_MOON_BARYCENTER_J2000, + .translate( + VENUS_J2000, // Target + EARTH_MOON_BARYCENTER_J2000, // Observer epoch, - Aberration::NotSet, + None, ) .unwrap(); @@ -242,7 +242,7 @@ def test_state_transformation(): assert orig_state.tlong_deg() == 0.6916999999999689 state_itrf93 = ctx.transform_to( - orig_state, Frames.EARTH_ITRF93, Aberration.NotSet + orig_state, Frames.EARTH_ITRF93, None ) print(orig_state) @@ -254,7 +254,7 @@ def test_state_transformation(): # Convert back from_state_itrf93_to_eme2k = ctx.transform_to( - state_itrf93, Frames.EARTH_J2000, Aberration.NotSet + state_itrf93, Frames.EARTH_J2000, None ) print(from_state_itrf93_to_eme2k) diff --git a/anise-py/src/lib.rs b/anise-py/src/lib.rs index 076f12ae..062f59dd 100644 --- a/anise-py/src/lib.rs +++ b/anise-py/src/lib.rs @@ -8,7 +8,7 @@ * Documentation: https://nyxspace.com/ */ -use ::anise::almanac::meta::MetaAlmanac; +use ::anise::almanac::meta::{MetaAlmanac, MetaFile}; use ::anise::almanac::Almanac; use ::anise::astro::Aberration; use hifitime::leap_seconds::{LatestLeapSeconds, LeapSecondsFile}; @@ -28,6 +28,7 @@ fn anise(py: Python, m: &PyModule) -> PyResult<()> { m.add_class::()?; m.add_class::()?; m.add_class::()?; + m.add_class::()?; Ok(()) } diff --git a/anise-py/tests/test_almanac.py b/anise-py/tests/test_almanac.py index 056e8ea6..1f43e49a 100644 --- a/anise-py/tests/test_almanac.py +++ b/anise-py/tests/test_almanac.py @@ -1,24 +1,40 @@ from pathlib import Path import pickle -from anise import Aberration, Almanac, MetaAlmanac +from anise import Almanac, MetaAlmanac from anise.astro import * from anise.astro.constants import Frames from anise.time import Epoch +from os import environ + def test_state_transformation(): """ This is the Python equivalent to anise/tests/almanac/mod.rs + but the data is loaded from the remote servers """ - data_path = Path(__file__).parent.joinpath("..", "..", "data") - # Must ensure that the path is a string - ctx = Almanac(str(data_path.joinpath("de440s.bsp"))) - # Let's add another file here -- note that the Almanac will load into a NEW variable, so we must overwrite it! - # This prevents memory leaks (yes, I promise) - ctx = ctx.load(str(data_path.joinpath("pck08.pca"))).load( - str(data_path.joinpath("earth_latest_high_prec.bpc")) - ) + + if environ.get("CI", False): + # # Load from meta kernel to not use Git LFS quota + # data_path = Path(__file__).parent.joinpath("..", "..", "data", "default_meta.dhall") + # meta = MetaAlmanac(str(data_path)) + # print(meta) + # # Process the files to be loaded + # ctx = meta.process() + print("I don't know where the files are in the Python CI") + + return + else: + data_path = Path(__file__).parent.joinpath("..", "..", "data") + # Must ensure that the path is a string + ctx = Almanac(str(data_path.joinpath("de440s.bsp"))) + # Let's add another file here -- note that the Almanac will load into a NEW variable, so we must overwrite it! + # This prevents memory leaks (yes, I promise) + ctx = ctx.load(str(data_path.joinpath("pck08.pca"))).load( + str(data_path.joinpath("earth_latest_high_prec.bpc")) + ) + eme2k = ctx.frame_info(Frames.EME2000) assert eme2k.mu_km3_s2() == 398600.435436096 assert eme2k.shape.polar_radius_km == 6356.75 @@ -43,7 +59,10 @@ def test_state_transformation(): assert abs(orig_state.raan_deg() - 306.614) < 1e-10 assert abs(orig_state.tlong_deg() - 0.6916999999999689) < 1e-10 - state_itrf93 = ctx.transform_to(orig_state, Frames.EARTH_ITRF93, Aberration.NotSet) + # In Python, we can set the aberration to None + aberration = None + + state_itrf93 = ctx.transform_to(orig_state, Frames.EARTH_ITRF93, aberration) print(orig_state) print(state_itrf93) @@ -54,7 +73,7 @@ def test_state_transformation(): # Convert back from_state_itrf93_to_eme2k = ctx.transform_to( - state_itrf93, Frames.EARTH_J2000, Aberration.NotSet + state_itrf93, Frames.EARTH_J2000, aberration ) print(from_state_itrf93_to_eme2k) @@ -102,6 +121,7 @@ def test_meta_load(): assert eme2k.shape.polar_radius_km == 6356.75 assert abs(eme2k.shape.flattening() - 0.0033536422844278) < 2e-16 + def test_exports(): for cls in [Frame, Ellipsoid, Orbit]: print(f"{cls} OK") diff --git a/anise/benches/crit_jpl_ephemerides.rs b/anise/benches/crit_jpl_ephemerides.rs index 368ae7e2..7612b755 100644 --- a/anise/benches/crit_jpl_ephemerides.rs +++ b/anise/benches/crit_jpl_ephemerides.rs @@ -22,7 +22,7 @@ fn benchmark_spice_single_hop_type2_cheby(time_it: TimeSeries) { fn benchmark_anise_single_hop_type2_cheby(ctx: &Almanac, time_it: TimeSeries) { for epoch in time_it { black_box( - ctx.translate_from_to_geometric(EARTH_J2000, LUNA_J2000, epoch) + ctx.translate_geometric(EARTH_J2000, LUNA_J2000, epoch) .unwrap(), ); } diff --git a/anise/benches/crit_spacecraft_ephemeris.rs b/anise/benches/crit_spacecraft_ephemeris.rs index 8c6d7898..3ed0148e 100644 --- a/anise/benches/crit_spacecraft_ephemeris.rs +++ b/anise/benches/crit_spacecraft_ephemeris.rs @@ -25,7 +25,7 @@ fn benchmark_anise_single_hop_type13_hermite(ctx: &Almanac, time_it: TimeSeries) let my_sc_j2k = Frame::from_ephem_j2000(-10000001); for epoch in time_it { black_box( - ctx.translate_from_to_geometric(my_sc_j2k, EARTH_J2000, epoch) + ctx.translate_geometric(my_sc_j2k, EARTH_J2000, epoch) .unwrap(), ); } diff --git a/anise/benches/iai_jpl_ephemerides.rs b/anise/benches/iai_jpl_ephemerides.rs index 79f8171c..07a49700 100644 --- a/anise/benches/iai_jpl_ephemerides.rs +++ b/anise/benches/iai_jpl_ephemerides.rs @@ -43,7 +43,7 @@ fn benchmark_anise_single_hop_type2_cheby() { for epoch in time_it { black_box( - ctx.translate_from_to_geometric(EARTH_J2000, LUNA_J2000, epoch) + ctx.translate_geometric(EARTH_J2000, LUNA_J2000, epoch) .unwrap(), ); } diff --git a/anise/benches/iai_spacecraft_ephemeris.rs b/anise/benches/iai_spacecraft_ephemeris.rs index b87758b7..3b43e458 100644 --- a/anise/benches/iai_spacecraft_ephemeris.rs +++ b/anise/benches/iai_spacecraft_ephemeris.rs @@ -37,7 +37,7 @@ fn benchmark_anise_single_hop_type13_hermite() { let my_sc_j2k = Frame::from_ephem_j2000(-10000001); black_box( - ctx.translate_from_to_geometric(my_sc_j2k, EARTH_J2000, epoch) + ctx.translate_geometric(my_sc_j2k, EARTH_J2000, epoch) .unwrap(), ); } diff --git a/anise/src/almanac/meta.rs b/anise/src/almanac/meta.rs index 1fe95b33..bbb79067 100644 --- a/anise/src/almanac/meta.rs +++ b/anise/src/almanac/meta.rs @@ -17,6 +17,7 @@ use snafu::prelude::*; use std::fs::{create_dir_all, File}; use std::io::Write; use std::path::Path; +use std::time::Duration; use url::Url; #[cfg(feature = "python")] @@ -65,12 +66,13 @@ pub enum MetaAlmanacError { #[derive(Clone, Debug, Serialize, Deserialize, PartialEq)] #[cfg_attr(feature = "python", pyclass)] #[cfg_attr(feature = "python", pyo3(module = "anise"))] +#[cfg_attr(feature = "python", pyo3(get_all, set_all))] pub struct MetaAlmanac { - files: Vec, + pub files: Vec, } impl MetaAlmanac { - /// Loads the provided path as a Dhall file and processes each file. + /// Loads the provided path as a Dhall configuration file and processes each file. pub fn new(path: String) -> Result { match serde_dhall::from_file(&path).parse::() { Err(e) => Err(MetaAlmanacError::ParseDhall { @@ -84,11 +86,14 @@ impl MetaAlmanac { #[cfg_attr(feature = "python", pymethods)] impl MetaAlmanac { - /// Loads the provided path as a Dhall file and processes each file. + /// Loads the provided path as a Dhall file. If no path is provided, creates an empty MetaAlmanac that can store MetaFiles. #[cfg(feature = "python")] #[new] - pub fn py_new(path: String) -> Result { - Self::new(path) + pub fn py_new(maybe_path: Option) -> Result { + match maybe_path { + Some(path) => Self::new(path), + None => Ok(Self { files: Vec::new() }), + } } /// Fetch all of the data and return a loaded Almanac @@ -182,14 +187,50 @@ impl Default for MetaAlmanac { } } +#[cfg_attr(feature = "python", pyclass)] +#[cfg_attr(feature = "python", pyo3(module = "anise"))] +#[cfg_attr(feature = "python", pyo3(get_all, set_all))] #[derive(Clone, Debug, Default, Serialize, Deserialize, PartialEq, StaticType)] pub struct MetaFile { + /// URI of this meta file pub uri: String, /// Optionally specify the CRC32 of this file, which will be checked prior to loading. pub crc32: Option, } +#[cfg_attr(feature = "python", pymethods)] impl MetaFile { + /// Builds a new MetaFile from the provided URI and optionally its CRC32 checksum. + #[cfg(feature = "python")] + #[new] + pub fn py_new(uri: String, crc32: Option) -> Self { + Self { uri, crc32 } + } + + #[cfg(feature = "python")] + fn __str__(&self) -> String { + format!("{self:?}") + } + + #[cfg(feature = "python")] + fn __repr__(&self) -> String { + format!("{self:?} (@{self:p})") + } + + #[cfg(feature = "python")] + fn __richcmp__(&self, other: &Self, op: CompareOp) -> Result { + match op { + CompareOp::Eq => Ok(self == other), + CompareOp::Ne => Ok(self != other), + _ => Err(PyErr::new::(format!( + "{op:?} not available" + ))), + } + } + + /// Processes this MetaFile by downloading it if it's a URL. + /// + /// This function modified `self` and changes the URI to be the path to the downloaded file. fn process(&mut self) -> Result<(), MetaAlmanacError> { match Url::parse(&self.uri) { Err(e) => { @@ -240,7 +281,13 @@ impl MetaFile { } // At this stage, either the dest path does not exist, or the CRC32 check failed. - match reqwest::blocking::get(url.clone()) { + let client = reqwest::blocking::Client::builder() + .connect_timeout(Duration::from_secs(30)) + .timeout(Duration::from_secs(30)) + .build() + .unwrap(); + + match client.get(url.clone()).send() { Ok(resp) => { if resp.status().is_success() { // Downloaded the file, let's store it locally. @@ -279,6 +326,7 @@ impl MetaFile { } } } else { + println!("err"); let err = resp.error_for_status().unwrap(); Err(MetaAlmanacError::FetchError { status: err.status(), diff --git a/anise/src/almanac/transform.rs b/anise/src/almanac/transform.rs index 63a52373..697f9176 100644 --- a/anise/src/almanac/transform.rs +++ b/anise/src/almanac/transform.rs @@ -35,11 +35,11 @@ impl Almanac { from_frame: Frame, to_frame: Frame, epoch: Epoch, - ab_corr: Aberration, + ab_corr: Option, ) -> Result { // Translate let state = self - .translate_from_to(from_frame, to_frame, epoch, ab_corr) + .translate(from_frame, to_frame, epoch, ab_corr) .with_context(|_| EphemerisSnafu { action: "transform from/to", })?; @@ -65,7 +65,7 @@ impl Almanac { &self, state: CartesianState, to_frame: Frame, - ab_corr: Aberration, + ab_corr: Option, ) -> Result { let state = self .translate_to(state, to_frame, ab_corr) @@ -96,7 +96,7 @@ impl Almanac { object: NaifId, observer: Frame, epoch: Epoch, - ab_corr: Aberration, + ab_corr: Option, ) -> Result { self.transform_from_to(Frame::from_ephem_j2000(object), observer, epoch, ab_corr) } @@ -114,7 +114,7 @@ impl Almanac { from_frame: Frame, to_frame: Frame, epoch: Epoch, - ab_corr: Aberration, + ab_corr: Option, distance_unit: LengthUnit, time_unit: TimeUnit, ) -> Result { diff --git a/anise/src/astro/aberration.rs b/anise/src/astro/aberration.rs new file mode 100644 index 00000000..b46871a4 --- /dev/null +++ b/anise/src/astro/aberration.rs @@ -0,0 +1,304 @@ +/* + * ANISE Toolkit + * Copyright (C) 2021-2023 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::{ + constants::SPEED_OF_LIGHT_KM_S, + errors::{AberrationSnafu, VelocitySnafu}, + math::{rotate_vector, Vector3}, +}; +use core::f64::EPSILON; +use core::fmt; + +#[cfg(feature = "python")] +use pyo3::prelude::*; +use snafu::ensure; + +use super::PhysicsResult; +use crate::errors::PhysicsError; + +/// Represents the aberration correction options in ANISE. +/// +/// In space science and engineering, accurately pointing instruments (like optical cameras or radio antennas) at a target is crucial. This task is complicated by the finite speed of light, necessitating corrections for the apparent position of the target. +/// +/// This structure holds parameters for aberration corrections applied to a target's position or state vector. These corrections account for the difference between the target's geometric (true) position and its apparent position as observed. +/// +/// # Rule of tumb +/// In most Earth orbits, one does _not_ need to provide any aberration corrections. Light time to the target is less than one second (the Moon is about one second away). +/// In near Earth orbits, e.g. inner solar system, preliminary analysis can benefit from enabling unconverged light time correction. Stellar aberration is probably not required. +/// For deep space missions, preliminary analysis would likely require both light time correction and stellar aberration. Mission planning and operations will definitely need converged light-time calculations. +/// +/// For more details, . +/// +/// # SPICE Validation +/// +/// The validation test `validate_jplde_de440s_aberration_lt` checks 101,000 pairs of ephemeris computations and shows that the unconverged Light Time computation matches the SPICE computations almost all the time. +/// More specifically, the 99th percentile of error is less than 5 meters, the 75th percentile is less than one meter, and the median error is less than 2 millimeters. +#[derive(Copy, Clone, Default, PartialEq, Eq)] +#[cfg_attr(feature = "python", pyclass)] +#[cfg_attr(feature = "python", pyo3(module = "anise"))] +#[cfg_attr(feature = "python", pyo3(get_all, set_all))] +pub struct Aberration { + /// Indicates whether the light time calculations should be iterated upon (more precise but three times as many CPU cycles). + pub converged: bool, + /// Flag to denote if stellar aberration correction is applied. Stellar aberration is due to the motion of the observer (caused by Earth's orbit around the Sun). + pub stellar: bool, + /// Specifies whether in reception or transmission mode. True for 'transmit' mode, indicating the correction is applied to the transmitted signal from the observer to the target. False for 'receive' mode, for signals received from the target. + pub transmit_mode: bool, +} + +impl Aberration { + /// Disables aberration corrections, e.g. all translations are geometric only (typical use case). + pub const NONE: Option = None; + /// Unconverged light time correction in reception mode without stellar aberration (e.g. a ground station targeting a spacecraft near the Moon) + pub const LT: Option = Some(Self { + converged: false, + stellar: false, + transmit_mode: false, + }); + /// Unconverged light time correction in reception mode with stellar aberration + pub const LT_S: Option = Some(Self { + converged: false, + stellar: true, + transmit_mode: false, + }); + /// Converged light time correction in reception mode without stellar aberration + pub const CN: Option = Some(Self { + converged: true, + stellar: false, + transmit_mode: false, + }); + /// Converged light time correction in reception mode with stellar aberration + pub const CN_S: Option = Some(Self { + converged: true, + stellar: true, + transmit_mode: false, + }); + /// Unconverged light time correction in transmission mode without stellar aberration (e.g. a Moon orbiter contacting a ground station) + pub const XLT: Option = Some(Self { + converged: false, + stellar: false, + transmit_mode: true, + }); + /// Unconverged light time correction in transmission mode with stellar aberration + pub const XLT_S: Option = Some(Self { + converged: false, + stellar: true, + transmit_mode: true, + }); + /// Converged light time correction in transmission mode without stellar aberration + pub const XCN: Option = Some(Self { + converged: true, + stellar: false, + transmit_mode: true, + }); + /// Converged light time correction in transmission mode with stellar aberration + pub const XCN_S: Option = Some(Self { + converged: true, + stellar: true, + transmit_mode: true, + }); + + /// Initializes a new Aberration structure from one of the following (SPICE compatibility): + /// + `NONE`: No correction + /// + `LT`: unconverged light time, no stellar aberration, reception mode + /// + `LT+S`: unconverged light time, with stellar aberration, reception mode + /// + `CN`: converged light time, no stellar aberration, reception mode + /// + `CN+S`: converged light time, with stellar aberration, reception mode + /// + `XLT`: unconverged light time, no stellar aberration, transmission mode + /// + `XLT+S`: unconverged light time, with stellar aberration, transmission mode + /// + `XCN`: converged light time, no stellar aberration, transmission mode + /// + `XCN+S`: converged light time, with stellar aberration, transmission mode + pub fn new(flag: &str) -> PhysicsResult> { + match flag.trim() { + "NONE" => Ok(Self::NONE), + "LT" => Ok(Self::LT), + "LT+S" => Ok(Self::LT_S), + "CN" => Ok(Self::CN), + "CN+S" => Ok(Self::CN_S), + "XLT" => Ok(Self::XLT), + "XLT+S" => Ok(Self::XLT_S), + "XCN" => Ok(Self::XCN), + "XCN+S" => Ok(Self::XCN_S), + _ => Err(PhysicsError::AberrationError { + action: "unknown aberration configuration name", + }), + } + } +} + +#[cfg(feature = "python")] +#[pymethods] +impl Aberration { + /// Initializes a new Aberration structure from one of the following (SPICE compatibility): + /// + `NONE`: No correction + /// + `LT`: unconverged light time, no stellar aberration, reception mode + /// + `LT+S`: unconverged light time, with stellar aberration, reception mode + /// + `CN`: converged light time, no stellar aberration, reception mode + /// + `CN+S`: converged light time, with stellar aberration, reception mode + /// + `XLT`: unconverged light time, no stellar aberration, transmission mode + /// + `XLT+S`: unconverged light time, with stellar aberration, transmission mode + /// + `XCN`: converged light time, no stellar aberration, transmission mode + /// + `XCN+S`: converged light time, with stellar aberration, transmission mode + #[new] + fn py_new(name: String) -> PhysicsResult { + match Self::new(&name)? { + Some(ab_corr) => Ok(ab_corr), + None => Err(PhysicsError::AberrationError { + action: "just uses `None` in Python instead", + }), + } + } + + fn __eq__(&self, other: &Self) -> bool { + self == other + } + + fn __str__(&self) -> String { + format!("{self}") + } + + fn __repr__(&self) -> String { + format!("{self:?} (@{self:p})") + } +} + +impl fmt::Debug for Aberration { + /// Prints this configuration as the SPICE name. + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + if self.transmit_mode { + write!(f, "X")?; + } + if self.converged { + write!(f, "CN")?; + } else { + write!(f, "LT")?; + } + if self.stellar { + write!(f, "+S")?; + } + Ok(()) + } +} + +impl fmt::Display for Aberration { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + if self.converged { + write!(f, "converged ")?; + } else { + write!(f, "unconverged ")?; + } + write!(f, "light-time ")?; + if self.stellar { + write!(f, "and stellar aberration")?; + } else { + write!(f, "aberration")?; + } + if self.transmit_mode { + write!(f, " in transmit mode")?; + } + Ok(()) + } +} + +/// Returns the provided target [Orbit] with respect to any observer corrected for steller aberration. +/// +/// # Arguments +/// +/// + `target_pos_km`: the position of a target object with respect to the observer in kilometers +/// + `obs_wrt_ssb_vel_km_s`: the velocity of the observer with respect to the Solar System Barycenter in kilometers per second +/// + `ab_corr`: the [Aberration] correction +/// +/// # Errors +/// +/// This function will return an error in the following cases: +/// 1. the aberration is not set to include stellar corrections; +/// 1. the `target` is moving faster than the speed of light. +/// +/// # Algorithm +/// Source: this algorithm and documentation were rewritten from NAIF's [`stelab`](https://github.com/nasa/kepler-pipeline/blob/f58b21df2c82969d8bd3e26a269bd7f5b9a770e1/source-code/matlab/fc/cspice-src-i686/cspice/stelab.c#L13) function: +/// +/// Let r be the vector from the observer to the object, and v be the velocity of the observer with respect to the Solar System barycenter. +/// Let w be the angle between them. The aberration angle phi is given by +/// +/// `sin(phi) = v sin(w) / c` +/// +/// Let h be the vector given by the cross product +/// +/// `h = r X v` +/// +/// Rotate r by phi radians about h to obtain the apparent position of the object. +/// +/// +pub fn stellar_aberration( + target_pos_km: Vector3, + obs_wrt_ssb_vel_km_s: Vector3, + ab_corr: Aberration, +) -> PhysicsResult { + ensure!( + ab_corr.stellar, + AberrationSnafu { + action: "stellar correction not available for this aberration" + } + ); + + // Obtain the negative of the observer's velocity. This velocity, combined + // with the target's position, will yield the inverse of the usual stellar + // aberration correction, which is exactly what we seek. + + let obs_velocity_km_s = if ab_corr.transmit_mode { + -obs_wrt_ssb_vel_km_s + } else { + obs_wrt_ssb_vel_km_s + }; + + // Get the velocity vector scaled with respect to the speed of light (v/c) + let vbyc = obs_velocity_km_s / SPEED_OF_LIGHT_KM_S; + + ensure!( + vbyc.dot(&vbyc) < 1.0, + VelocitySnafu { + action: "observer moving at or faster than light, cannot compute stellar aberration" + } + ); + + // Get a unit vector that points in the direction of the object + let u = target_pos_km.normalize(); + + // Compute u_obj x (v/c) + let h = u.cross(&vbyc); + + // Correct for stellar aberration + let mut app_target_pos_km = target_pos_km; + let sin_phi = h.norm(); + if sin_phi > EPSILON { + let phi = sin_phi.asin(); + app_target_pos_km = rotate_vector(&target_pos_km, &h, phi); + } + + Ok(app_target_pos_km) +} + +#[cfg(test)] +mod ut_aberration { + #[test] + fn test_display() { + use super::Aberration; + + assert_eq!(format!("{:?}", Aberration::LT.unwrap()), "LT"); + assert_eq!(format!("{:?}", Aberration::LT_S.unwrap()), "LT+S"); + assert_eq!(format!("{:?}", Aberration::CN.unwrap()), "CN"); + assert_eq!(format!("{:?}", Aberration::CN_S.unwrap()), "CN+S"); + + assert_eq!(format!("{:?}", Aberration::XLT.unwrap()), "XLT"); + assert_eq!(format!("{:?}", Aberration::XLT_S.unwrap()), "XLT+S"); + assert_eq!(format!("{:?}", Aberration::XCN.unwrap()), "XCN"); + assert_eq!(format!("{:?}", Aberration::XCN_S.unwrap()), "XCN+S"); + } +} diff --git a/anise/src/astro/mod.rs b/anise/src/astro/mod.rs index e808c7c4..8ffeb996 100644 --- a/anise/src/astro/mod.rs +++ b/anise/src/astro/mod.rs @@ -10,17 +10,8 @@ use crate::errors::PhysicsError; -#[cfg(feature = "python")] -use pyo3::prelude::*; - -/// Defines the aberration corrections to the state of the target body to account for one-way light time and stellar aberration. -/// **WARNING:** This enum is a placeholder until [https://github.com/anise-toolkit/anise.rs/issues/26] is implemented. -#[derive(Copy, Clone, Debug, PartialEq, Eq)] -#[cfg_attr(feature = "python", pyclass)] -#[cfg_attr(feature = "python", pyo3(module = "anise."))] -pub enum Aberration { - NotSet, -} +pub(crate) mod aberration; +pub use aberration::Aberration; pub mod orbit; pub mod orbit_geodetic; diff --git a/anise/src/astro/orbit.rs b/anise/src/astro/orbit.rs index 2d4c9404..7cb1631b 100644 --- a/anise/src/astro/orbit.rs +++ b/anise/src/astro/orbit.rs @@ -12,7 +12,7 @@ use super::PhysicsResult; use crate::{ errors::{ HyperbolicTrueAnomalySnafu, InfiniteValueSnafu, ParabolicEccentricitySnafu, - ParabolicSemiParamSnafu, PhysicsError, RadiusSnafu, + ParabolicSemiParamSnafu, PhysicsError, RadiusSnafu, VelocitySnafu, }, math::{ angles::{between_0_360, between_pm_180}, @@ -251,7 +251,7 @@ impl CartesianState { ); ensure!( self.vmag_km_s() > EPSILON, - RadiusSnafu { + VelocitySnafu { action: "cannot compute orbital momentum vector with zero velocity" } ); diff --git a/anise/src/constants.rs b/anise/src/constants.rs index e5dee64b..f55f8d87 100644 --- a/anise/src/constants.rs +++ b/anise/src/constants.rs @@ -8,6 +8,9 @@ * Documentation: https://nyxspace.com/ */ +/// Speed of light in kilometers per second (km/s) +pub const SPEED_OF_LIGHT_KM_S: f64 = 299_792.458; + pub mod celestial_objects { use crate::NaifId; diff --git a/anise/src/ephemerides/mod.rs b/anise/src/ephemerides/mod.rs index b1dda0ef..ed5c5b45 100644 --- a/anise/src/ephemerides/mod.rs +++ b/anise/src/ephemerides/mod.rs @@ -45,8 +45,9 @@ pub enum EphemerisError { #[snafu(backtrace)] source: DAFError, }, - #[snafu(display("during an ephemeris operation: {source}"))] + #[snafu(display("when {action} for ephemeris {source}"))] EphemerisPhysics { + action: &'static str, #[snafu(backtrace)] source: PhysicsError, }, diff --git a/anise/src/ephemerides/translate_to_parent.rs b/anise/src/ephemerides/translate_to_parent.rs index b4260d71..d43a07f8 100644 --- a/anise/src/ephemerides/translate_to_parent.rs +++ b/anise/src/ephemerides/translate_to_parent.rs @@ -13,7 +13,6 @@ use snafu::ResultExt; use super::{EphemerisError, SPKSnafu}; use crate::almanac::Almanac; -use crate::astro::Aberration; use crate::ephemerides::EphemInterpolationSnafu; use crate::hifitime::Epoch; use crate::math::cartesian::CartesianState; @@ -35,7 +34,6 @@ impl Almanac { &self, source: Frame, epoch: Epoch, - _ab_corr: Aberration, ) -> Result<(Vector3, Vector3, Frame), EphemerisError> { // First, let's find the SPK summary for this frame. let (summary, spk_no, idx_in_spk) = @@ -51,6 +49,7 @@ impl Almanac { .ok_or(EphemerisError::Unreachable)?; // Now let's simply evaluate the data + let (pos_km, vel_km_s) = match summary.data_type()? { DafDataType::Type2ChebyshevTriplet => { let data = spk_data @@ -93,14 +92,13 @@ impl Almanac { Ok((pos_km, vel_km_s, new_frame)) } + /// Performs the GEOMETRIC translation to the parent. Use translate_from_to for aberration. pub fn translate_to_parent( &self, source: Frame, epoch: Epoch, - ab_corr: Aberration, ) -> Result { - let (radius_km, velocity_km_s, frame) = - self.translation_parts_to_parent(source, epoch, ab_corr)?; + let (radius_km, velocity_km_s, frame) = self.translation_parts_to_parent(source, epoch)?; Ok(CartesianState { radius_km, diff --git a/anise/src/ephemerides/translations.rs b/anise/src/ephemerides/translations.rs index 93907555..ec21fcae 100644 --- a/anise/src/ephemerides/translations.rs +++ b/anise/src/ephemerides/translations.rs @@ -13,7 +13,10 @@ use snafu::ResultExt; use super::EphemerisError; use super::EphemerisPhysicsSnafu; use crate::almanac::Almanac; +use crate::astro::aberration::stellar_aberration; use crate::astro::Aberration; +use crate::constants::frames::SSB_J2000; +use crate::constants::SPEED_OF_LIGHT_KM_S; use crate::hifitime::Epoch; use crate::math::cartesian::CartesianState; use crate::math::units::*; @@ -28,106 +31,172 @@ use pyo3::prelude::*; #[cfg_attr(feature = "python", pymethods)] impl Almanac { - /// Returns the Cartesian state needed to translate the `from_frame` to the `to_frame`. + /// Returns the Cartesian state of the target frame as seen from the observer frame at the provided epoch, and optionally given the aberration correction. + /// + /// # SPICE Compatibility + /// This function is the SPICE equivalent of spkezr: `spkezr(TARGET_ID, EPOCH_TDB_S, ORIENTATION_ID, ABERRATION, OBSERVER_ID)` + /// In ANISE, the TARGET_ID and ORIENTATION are provided in the first argument (TARGET_FRAME), as that frame includes BOTH + /// the target ID and the orientation of that target. The EPOCH_TDB_S is the epoch in the TDB time system, which is computed + /// in ANISE using Hifitime. THe ABERRATION is computed by providing the optional Aberration flag. Finally, the OBSERVER + /// argument is replaced by OBSERVER_FRAME: if the OBSERVER_FRAME argument has the same orientation as the TARGET_FRAME, then this call + /// will return exactly the same data as the spkerz SPICE call. /// /// # Warning - /// This function only performs the translation and no rotation whatsoever. Use the `transform_from_to` function instead to include rotations. + /// This function only performs the translation and no rotation whatsoever. Use the `transform` function instead to include rotations. /// /// # Note /// This function performs a recursion of no more than twice the [MAX_TREE_DEPTH]. - pub fn translate_from_to( + pub fn translate( &self, - from_frame: Frame, - to_frame: Frame, + target_frame: Frame, + observer_frame: Frame, epoch: Epoch, - ab_corr: Aberration, + ab_corr: Option, ) -> Result { - let mut to_frame: Frame = to_frame; - - // If there is no frame info, the user hasn't loaded this frame, but might still want to compute a translation. - if let Ok(to_frame_info) = self.frame_from_uid(to_frame) { - // User has loaded the planetary data for this frame, so let's use that as the to_frame. - to_frame = to_frame_info; - } - - if from_frame == to_frame { + if observer_frame == target_frame { // Both frames match, return this frame's hash (i.e. no need to go higher up). - return Ok(CartesianState::zero(from_frame)); + return Ok(CartesianState::zero(observer_frame)); } - let (node_count, path, common_node) = - self.common_ephemeris_path(from_frame, to_frame, epoch)?; - - // The fwrd variables are the states from the `from frame` to the common node - let (mut pos_fwrd, mut vel_fwrd, mut frame_fwrd) = - if from_frame.ephem_origin_id_match(common_node) { - (Vector3::zeros(), Vector3::zeros(), from_frame) - } else { - self.translation_parts_to_parent(from_frame, epoch, ab_corr)? - }; - - // The bwrd variables are the states from the `to frame` back to the common node - let (mut pos_bwrd, mut vel_bwrd, mut frame_bwrd) = - if to_frame.ephem_origin_id_match(common_node) { - (Vector3::zeros(), Vector3::zeros(), to_frame) - } else { - self.translation_parts_to_parent(to_frame, epoch, ab_corr)? - }; - - for cur_node_id in path.iter().take(node_count) { - if !frame_fwrd.ephem_origin_id_match(common_node) { - let (cur_pos_fwrd, cur_vel_fwrd, cur_frame_fwrd) = - self.translation_parts_to_parent(frame_fwrd, epoch, ab_corr)?; - - pos_fwrd += cur_pos_fwrd; - vel_fwrd += cur_vel_fwrd; - frame_fwrd = cur_frame_fwrd; - } + let mut obs_frame: Frame = observer_frame; - if !frame_bwrd.ephem_origin_id_match(common_node) { - let (cur_pos_bwrd, cur_vel_bwrd, cur_frame_bwrd) = - self.translation_parts_to_parent(frame_bwrd, epoch, ab_corr)?; + // If there is no frame info, the user hasn't loaded this frame, but might still want to compute a translation. + if let Ok(obs_frame_info) = self.frame_from_uid(obs_frame) { + // User has loaded the planetary data for this frame, so let's use that as the to_frame. + obs_frame = obs_frame_info; + } - pos_bwrd += cur_pos_bwrd; - vel_bwrd += cur_vel_bwrd; - frame_bwrd = cur_frame_bwrd; + match ab_corr { + None => { + let (node_count, path, common_node) = + self.common_ephemeris_path(observer_frame, target_frame, epoch)?; + + // The fwrd variables are the states from the `from frame` to the common node + let (mut pos_fwrd, mut vel_fwrd, mut frame_fwrd) = + if observer_frame.ephem_origin_id_match(common_node) { + (Vector3::zeros(), Vector3::zeros(), observer_frame) + } else { + self.translation_parts_to_parent(observer_frame, epoch)? + }; + + // The bwrd variables are the states from the `to frame` back to the common node + let (mut pos_bwrd, mut vel_bwrd, mut frame_bwrd) = + if target_frame.ephem_origin_id_match(common_node) { + (Vector3::zeros(), Vector3::zeros(), target_frame) + } else { + self.translation_parts_to_parent(target_frame, epoch)? + }; + + for cur_node_id in path.iter().take(node_count) { + if !frame_fwrd.ephem_origin_id_match(common_node) { + let (cur_pos_fwrd, cur_vel_fwrd, cur_frame_fwrd) = + self.translation_parts_to_parent(frame_fwrd, epoch)?; + + pos_fwrd += cur_pos_fwrd; + vel_fwrd += cur_vel_fwrd; + frame_fwrd = cur_frame_fwrd; + } + + if !frame_bwrd.ephem_origin_id_match(common_node) { + let (cur_pos_bwrd, cur_vel_bwrd, cur_frame_bwrd) = + self.translation_parts_to_parent(frame_bwrd, epoch)?; + + pos_bwrd += cur_pos_bwrd; + vel_bwrd += cur_vel_bwrd; + frame_bwrd = cur_frame_bwrd; + } + + // We know this exist, so we can safely unwrap it + if cur_node_id.unwrap() == common_node { + break; + } + } + + Ok(CartesianState { + radius_km: pos_bwrd - pos_fwrd, + velocity_km_s: vel_bwrd - vel_fwrd, + epoch, + frame: obs_frame.with_orient(target_frame.orientation_id), + }) } - - // We know this exist, so we can safely unwrap it - if cur_node_id.unwrap() == common_node { - break; + Some(ab_corr) => { + // This is a rewrite of NAIF SPICE's `spkapo` + + // Find the geometric position of the observer body with respect to the solar system barycenter. + let obs_ssb = self.translate(observer_frame, SSB_J2000, epoch, None)?; + let obs_ssb_pos_km = obs_ssb.radius_km; + let obs_ssb_vel_km_s = obs_ssb.velocity_km_s; + + // Find the geometric position of the target body with respect to the solar system barycenter. + let tgt_ssb = self.translate(target_frame, SSB_J2000, epoch, None)?; + let tgt_ssb_pos_km = tgt_ssb.radius_km; + let tgt_ssb_vel_km_s = tgt_ssb.velocity_km_s; + + // Subtract the position of the observer to get the relative position. + let mut rel_pos_km = tgt_ssb_pos_km - obs_ssb_pos_km; + // NOTE: We never correct the velocity, so the geometric velocity is what we're seeking. + let mut rel_vel_km_s = tgt_ssb_vel_km_s - obs_ssb_vel_km_s; + + // Use this to compute the one-way light time in seconds. + let mut one_way_lt_s = rel_pos_km.norm() / SPEED_OF_LIGHT_KM_S; + + // To correct for light time, find the position of the target body at the current epoch + // minus the one-way light time. Note that the observer remains where he is. + + let num_it = if ab_corr.converged { 3 } else { 1 }; + let lt_sign = if ab_corr.transmit_mode { 1.0 } else { -1.0 }; + + for _ in 0..num_it { + let epoch_lt = epoch + lt_sign * one_way_lt_s * TimeUnit::Second; + let tgt_ssb = self.translate(target_frame, SSB_J2000, epoch_lt, None)?; + let tgt_ssb_pos_km = tgt_ssb.radius_km; + let tgt_ssb_vel_km_s = tgt_ssb.velocity_km_s; + + rel_pos_km = tgt_ssb_pos_km - obs_ssb_pos_km; + rel_vel_km_s = tgt_ssb_vel_km_s - obs_ssb_vel_km_s; + one_way_lt_s = rel_pos_km.norm() / SPEED_OF_LIGHT_KM_S; + } + + // If stellar aberration correction is requested, perform it now. + if ab_corr.stellar { + // Modifications based on transmission versus reception case is done in the function directly. + rel_pos_km = stellar_aberration(rel_pos_km, obs_ssb_vel_km_s, ab_corr) + .with_context(|_| EphemerisPhysicsSnafu { + action: "computing stellar aberration", + })?; + } + + Ok(CartesianState { + radius_km: rel_pos_km, + velocity_km_s: rel_vel_km_s, + epoch, + frame: observer_frame.with_orient(target_frame.orientation_id), + }) } } - - Ok(CartesianState { - radius_km: pos_fwrd - pos_bwrd, - velocity_km_s: vel_fwrd - vel_bwrd, - epoch, - frame: to_frame.with_orient(from_frame.orientation_id), - }) } /// Returns the geometric position vector, velocity vector, and acceleration vector needed to translate the `from_frame` to the `to_frame`, where the distance is in km, the velocity in km/s, and the acceleration in km/s^2. - pub fn translate_from_to_geometric( + pub fn translate_geometric( &self, - from_frame: Frame, - to_frame: Frame, + target_frame: Frame, + observer_frame: Frame, epoch: Epoch, ) -> Result { - self.translate_from_to(from_frame, to_frame, epoch, Aberration::NotSet) + self.translate(target_frame, observer_frame, epoch, Aberration::NONE) } - /// Translates the provided Cartesian state into the requested frame + /// Translates the provided Cartesian state into the requested observer frame /// /// **WARNING:** This function only performs the translation and no rotation _whatsoever_. Use the [transform_to] function instead to include rotations. #[allow(clippy::too_many_arguments)] pub fn translate_to( &self, state: CartesianState, - to_frame: Frame, - ab_corr: Aberration, + observer_frame: Frame, + ab_corr: Option, ) -> Result { - let frame_state = self.translate_from_to(state.frame, to_frame, state.epoch, ab_corr)?; + let frame_state = self.translate(state.frame, observer_frame, state.epoch, ab_corr)?; Ok(state.add_unchecked(&frame_state)) } @@ -143,14 +212,14 @@ impl Almanac { position: Vector3, velocity: Vector3, from_frame: Frame, - to_frame: Frame, + observer_frame: Frame, epoch: Epoch, - ab_corr: Aberration, + ab_corr: Option, distance_unit: LengthUnit, time_unit: TimeUnit, ) -> Result { // Compute the frame translation - let frame_state = self.translate_from_to(from_frame, to_frame, epoch, ab_corr)?; + let frame_state = self.translate(from_frame, observer_frame, epoch, ab_corr)?; let dist_unit_factor = LengthUnit::Kilometer.from_meters() * distance_unit.to_meters(); let time_unit_factor = time_unit.in_seconds(); @@ -162,6 +231,8 @@ impl Almanac { frame: from_frame, }; - (input_state + frame_state).with_context(|_| EphemerisPhysicsSnafu {}) + (input_state + frame_state).with_context(|_| EphemerisPhysicsSnafu { + action: "translating states (likely a bug!)", + }) } } diff --git a/anise/src/errors.rs b/anise/src/errors.rs index fc8d665e..6c3865c7 100644 --- a/anise/src/errors.rs +++ b/anise/src/errors.rs @@ -186,10 +186,12 @@ pub enum PhysicsError { InfiniteValue { action: &'static str }, #[snafu(display("{source}"))] AppliedMath { source: MathError }, - #[snafu(display("invalid radius {action}"))] + #[snafu(display("invalid radius: {action}"))] RadiusError { action: &'static str }, - #[snafu(display("invalid velocity {action}"))] + #[snafu(display("invalid velocity: {action}"))] VelocityError { action: &'static str }, + #[snafu(display("invalid aberration: {action}"))] + AberrationError { action: &'static str }, } impl From for InputOutputError { diff --git a/anise/src/math/cartesian.rs b/anise/src/math/cartesian.rs index d3c1537d..cfd67a14 100644 --- a/anise/src/math/cartesian.rs +++ b/anise/src/math/cartesian.rs @@ -8,15 +8,16 @@ * Documentation: https://nyxspace.com/ */ -use super::{perpv, Vector3}; +use super::{perp_vector, Vector3}; use crate::{ astro::PhysicsResult, + constants::SPEED_OF_LIGHT_KM_S, errors::{EpochMismatchSnafu, FrameMismatchSnafu, PhysicsError}, prelude::Frame, }; use core::fmt; use core::ops::{Add, Neg, Sub}; -use hifitime::Epoch; +use hifitime::{Duration, Epoch, TimeUnits}; use nalgebra::Vector6; use snafu::ensure; @@ -141,7 +142,7 @@ impl CartesianState { /// Returns the unit vector in the direction of the state velocity pub fn v_hat(&self) -> Vector3 { - perpv(&self.velocity_km_s, &self.r_hat()) / self.rmag_km() + perp_vector(&self.velocity_km_s, &self.r_hat()) / self.rmag_km() } /// Adds the other state to this state WITHOUT checking if the frames match. @@ -179,7 +180,7 @@ impl CartesianState { } /// Returns the distance in kilometers between this state and another state, if both frame match (epoch does not need to match). - pub fn distance_to(&self, other: &Self) -> PhysicsResult { + pub fn distance_to_km(&self, other: &Self) -> PhysicsResult { ensure!( self.frame == other.frame, FrameMismatchSnafu { @@ -203,6 +204,11 @@ impl CartesianState { && (self.velocity_km_s.z - other.velocity_km_s.z).abs() < velocity_tol_km_s && self.frame == other.frame } + + /// Returns the light time duration between this object and the origin of its reference frame. + pub fn light_time(&self) -> Duration { + (self.radius_km.norm() / SPEED_OF_LIGHT_KM_S).seconds() + } } impl Add for CartesianState { @@ -323,7 +329,7 @@ impl fmt::LowerExp for CartesianState { mod cartesian_state_ut { use std::f64::EPSILON; - use hifitime::{Epoch, TimeUnits}; + use hifitime::{Duration, Epoch, TimeUnits}; use crate::constants::frames::{EARTH_J2000, VENUS_J2000}; use crate::errors::PhysicsError; @@ -388,7 +394,7 @@ mod cartesian_state_ut { let s1 = CartesianState::new(10.0, 20.0, 30.0, 1.0, 2.0, 2.0, e, frame); let s2 = CartesianState::new(10.0, 20.0, 30.0, 1.0, 2.0, 2.0, e, frame); - assert!(s1.distance_to(&s2).unwrap().abs() < EPSILON); + assert!(s1.distance_to_km(&s2).unwrap().abs() < EPSILON); let as_vec6 = Vector6::new(10.0, 20.0, 30.0, 1.0, 2.0, 2.0); assert_eq!(s1.to_cartesian_pos_vel(), as_vec6); @@ -410,5 +416,7 @@ mod cartesian_state_ut { let s = CartesianState::zero_at_epoch(e, frame); assert!(s.hmag().is_err()); + + assert_eq!(s.light_time(), Duration::ZERO); } } diff --git a/anise/src/math/mod.rs b/anise/src/math/mod.rs index 1270ce56..b2c84933 100644 --- a/anise/src/math/mod.rs +++ b/anise/src/math/mod.rs @@ -24,12 +24,14 @@ pub mod rotation; pub mod units; /// Returns the projection of a onto b -pub fn projv(a: &Vector3, b: &Vector3) -> Vector3 { +/// Converted from NAIF SPICE's `projv` +pub fn project_vector(a: &Vector3, b: &Vector3) -> Vector3 { b * a.dot(b) / b.dot(b) } /// Returns the components of vector a orthogonal to b -pub fn perpv(a: &Vector3, b: &Vector3) -> Vector3 { +/// Converted from NAIF SPICE's `prepv` +pub fn perp_vector(a: &Vector3, b: &Vector3) -> Vector3 { let big_a = a[0].abs().max(a[1].abs().max(a[2].abs())); let big_b = b[0].abs().max(b[1].abs().max(b[2].abs())); if big_a < f64::EPSILON { @@ -39,7 +41,29 @@ pub fn perpv(a: &Vector3, b: &Vector3) -> Vector3 { } else { let a_scl = a / big_a; let b_scl = b / big_b; - let v = projv(&a_scl, &b_scl); + let v = project_vector(&a_scl, &b_scl); big_a * (a_scl - v) } } + +/// Rotate the vector a around the provided axis by angle theta. +pub fn rotate_vector(a: &Vector3, axis: &Vector3, theta_rad: f64) -> Vector3 { + let k_hat = axis.normalize(); + a.scale(theta_rad.cos()) + + k_hat.cross(a).scale(theta_rad.sin()) + + k_hat.scale(k_hat.dot(a) * (1.0 - theta_rad.cos())) +} + +#[cfg(test)] +mod math_ut { + use super::{rotate_vector, Vector3}; + #[test] + fn test_rotate_vector() { + use approx::assert_abs_diff_eq; + let a = Vector3::new(1.0, 0.0, 0.0); + let axis = Vector3::new(0.0, 0.0, 1.0); + let theta_rad = std::f64::consts::PI / 2.0; + let result = rotate_vector(&a, &axis, theta_rad); + assert_abs_diff_eq!(result, Vector3::new(0.0, 1.0, 0.0), epsilon = 1e-7); + } +} diff --git a/anise/tests/almanac/mod.rs b/anise/tests/almanac/mod.rs index 43e51bdf..285584f3 100644 --- a/anise/tests/almanac/mod.rs +++ b/anise/tests/almanac/mod.rs @@ -53,15 +53,16 @@ fn test_state_transformation() { // Transform that into another frame. let state_itrf93 = almanac - .transform_to(orig_state, EARTH_ITRF93, Aberration::NotSet) + .transform_to(orig_state, EARTH_ITRF93, Aberration::NONE) .unwrap(); println!("{orig_state:x}"); println!("{state_itrf93:X}"); - // Convert back + // Convert back. + // Note that the Aberration correction constants are actually options! let from_state_itrf93_to_eme2k = almanac - .transform_to(state_itrf93, EARTH_J2000, Aberration::NotSet) + .transform_to(state_itrf93, EARTH_J2000, None) .unwrap(); println!("{from_state_itrf93_to_eme2k}"); diff --git a/anise/tests/ephemerides/parent_translation_verif.rs b/anise/tests/ephemerides/parent_translation_verif.rs index cd7670a9..7e688766 100644 --- a/anise/tests/ephemerides/parent_translation_verif.rs +++ b/anise/tests/ephemerides/parent_translation_verif.rs @@ -23,7 +23,7 @@ fn invalid_load_from_static() { } #[test] -fn de438s_parent_translation_verif() { +fn de440s_parent_translation_verif() { let _ = pretty_env_logger::try_init(); let bytes = file2heap!("../data/de440s.bsp").unwrap(); @@ -44,9 +44,7 @@ fn de438s_parent_translation_verif() { ['9.5205530594596043e+07', '-4.6160758818180226e+07', '-2.6779476581501361e+07', '1.6612048969243794e+01', '2.8272067093941200e+01', '1.1668575714409423e+01'] */ - let state = ctx - .translate_to_parent(VENUS_J2000, epoch, Aberration::NotSet) - .unwrap(); + let state = ctx.translate_to_parent(VENUS_J2000, epoch).unwrap(); let pos_km = state.radius_km; let vel_km_s = state.velocity_km_s; @@ -63,10 +61,6 @@ fn de438s_parent_translation_verif() { 1.166_857_571_440_942_3e1, ); - // We expect exactly the same output as SPICE to machine precision. - assert!((pos_km - pos_expct_km).norm() < EPSILON); - assert!((vel_km_s - vel_expct_km_s).norm() < EPSILON); - // We expect exactly the same output as SPICE to machine precision. assert!( (pos_km - pos_expct_km).norm() < EPSILON, diff --git a/anise/tests/ephemerides/transform.rs b/anise/tests/ephemerides/transform.rs index 6397c4e1..f8af9e80 100644 --- a/anise/tests/ephemerides/transform.rs +++ b/anise/tests/ephemerides/transform.rs @@ -42,7 +42,7 @@ fn de440s_transform_verif_venus2emb() { let epoch = Epoch::from_gregorian_utc_at_midnight(2020, 2, 7); let state = almanac - .transform_from_to(VENUS_J2000, EARTH_ITRF93, epoch, Aberration::NotSet) + .transform_from_to(VENUS_J2000, EARTH_ITRF93, epoch, Aberration::NONE) .unwrap(); let (spice_state, _) = spice::spkezr("VENUS", epoch.to_et_seconds(), "ITRF93", "NONE", "EARTH"); diff --git a/anise/tests/ephemerides/translation.rs b/anise/tests/ephemerides/translation.rs index 1c89ad27..fa8c5f27 100644 --- a/anise/tests/ephemerides/translation.rs +++ b/anise/tests/ephemerides/translation.rs @@ -19,6 +19,8 @@ use anise::prelude::*; const POSITION_EPSILON_KM: f64 = 2e-8; // Corresponds to an error of 5e-6 meters per second, or 5.0 micrometers per second const VELOCITY_EPSILON_KM_S: f64 = 5e-9; +// Light time velocity error is too large! Cf. https://github.com/nyx-space/anise/issues/157 +const ABERRATION_VELOCITY_EPSILON_KM_S: f64 = 1e-4; #[test] fn de440s_translation_verif_venus2emb() { @@ -47,11 +49,11 @@ fn de440s_translation_verif_venus2emb() { */ let state = ctx - .translate_from_to( + .translate( VENUS_J2000, EARTH_MOON_BARYCENTER_J2000, epoch, - Aberration::NotSet, + Aberration::NONE, ) .unwrap(); @@ -84,7 +86,7 @@ fn de440s_translation_verif_venus2emb() { // Test the opposite translation let state = ctx - .translate_from_to_geometric(EARTH_MOON_BARYCENTER_J2000, VENUS_J2000, epoch) + .translate_geometric(EARTH_MOON_BARYCENTER_J2000, VENUS_J2000, epoch) .unwrap(); // We expect exactly the same output as SPICE to machine precision. @@ -134,7 +136,7 @@ fn de438s_translation_verif_venus2luna() { */ let state = ctx - .translate_from_to(VENUS_J2000, LUNA_J2000, epoch, Aberration::NotSet) + .translate(VENUS_J2000, LUNA_J2000, epoch, Aberration::NONE) .unwrap(); let pos_expct_km = Vector3::new( @@ -170,7 +172,7 @@ fn de438s_translation_verif_venus2luna() { // Test the opposite translation let state = ctx - .translate_from_to_geometric(LUNA_J2000, VENUS_J2000, epoch) + .translate_geometric(LUNA_J2000, VENUS_J2000, epoch) .unwrap(); // We expect exactly the same output as SPICE to machine precision. @@ -223,11 +225,11 @@ fn de438s_translation_verif_emb2luna() { */ let state = ctx - .translate_from_to( + .translate( EARTH_MOON_BARYCENTER_J2000, LUNA_J2000, epoch, - Aberration::NotSet, + Aberration::NONE, ) .unwrap(); @@ -267,11 +269,11 @@ fn de438s_translation_verif_emb2luna() { // Try the opposite let state = ctx - .translate_from_to( + .translate( LUNA_J2000, EARTH_MOON_BARYCENTER_J2000, epoch, - Aberration::NotSet, + Aberration::NONE, ) .unwrap(); @@ -322,7 +324,7 @@ fn spk_hermite_type13_verif() { let my_sc_j2k = Frame::from_ephem_j2000(-10000001); let state = ctx - .translate_from_to_geometric(my_sc_j2k, EARTH_J2000, epoch) + .translate_geometric(my_sc_j2k, EARTH_J2000, epoch) .unwrap(); println!("{state:?}"); @@ -381,7 +383,7 @@ fn multithread_query() { let epochs: Vec = time_it.collect(); epochs.into_par_iter().for_each(|epoch| { let state = ctx - .translate_from_to_geometric(LUNA_J2000, EARTH_MOON_BARYCENTER_J2000, epoch) + .translate_geometric(LUNA_J2000, EARTH_MOON_BARYCENTER_J2000, epoch) .unwrap(); println!("{state:?}"); }); @@ -417,11 +419,11 @@ fn hermite_query() { // Query in the middle to the parent, since we don't have anything else loaded. let state = ctx - .translate_from_to( + .translate( summary.target_frame(), summary.center_frame(), summary.start_epoch() + summary_duration * 0.5, - Aberration::NotSet, + Aberration::NONE, ) .unwrap(); @@ -430,11 +432,11 @@ fn hermite_query() { // Fetch the state at the start of this spline to make sure we don't glitch. assert!(ctx - .translate_from_to( + .translate( summary.target_frame(), summary.center_frame(), summary.start_epoch(), - Aberration::NotSet, + Aberration::NONE, ) .is_ok()); @@ -450,3 +452,189 @@ fn hermite_query() { // ) // .is_ok()); } + +/// This tests that the rotation from Moon to Earth matches SPICE with different aberration corrections. +/// We test Moon->Earth Moon Barycenter (instead of Venus->SSB as above) because there is no stellar correction possible +/// when the parent is the solar system barycenter. +#[test] +fn de440s_translation_verif_aberrations() { + let _ = pretty_env_logger::try_init(); + + let ctx = Almanac::new("../data/de440s.bsp").unwrap(); + + let epoch = Epoch::from_gregorian_utc_at_midnight(2002, 2, 7); + + /* + Python code: + >>> import spiceypy as sp + >>> sp.furnsh('data/de440s.bsp') + >>> et = 66312064.18493876 + >>> ['{:.16e}'.format(x) for x in sp.spkez(301, et, "J2000", "LT", 3)[0]]: + ['-8.1551741540104151e+04', + '-3.4544933489888906e+05', + '-1.4438031089871377e+05', + '9.6070843890026225e-01', + '-2.0357817054602378e-01', + '-1.8380326019667059e-01'] + */ + + struct AberrationCase { + correction: Option, + pos_expct_km: Vector3, + vel_expct_km_s: Vector3, + } + + let cases = [ + AberrationCase { + correction: Aberration::LT, + pos_expct_km: Vector3::new( + -8.1551741540104151e+04, + -3.4544933489888906e+05, + -1.4438031089871377e+05, + ), + vel_expct_km_s: Vector3::new( + 9.6070843890026225e-01, + -2.0357817054602378e-01, + -1.8380326019667059e-01, + ), + }, + AberrationCase { + correction: Aberration::LT_S, + pos_expct_km: Vector3::new( + -8.1570721849324545e+04, + -3.4544537500374130e+05, + -1.4437906334030110e+05, + ), + vel_expct_km_s: Vector3::new( + 9.6061748706693784e-01, + -2.0361038608395909e-01, + -1.8380826287127400e-01, + ), + }, + AberrationCase { + correction: Aberration::CN, + pos_expct_km: Vector3::new( + -8.1551743705525994e+04, + -3.4544933719548583e+05, + -1.4438031190508604e+05, + ), + vel_expct_km_s: Vector3::new( + 9.6070843946986884e-01, + -2.0357817069716688e-01, + -1.8380326026637128e-01, + ), + }, + AberrationCase { + correction: Aberration::CN_S, + pos_expct_km: Vector3::new( + -8.1570724014738982e+04, + -3.4544537730026408e+05, + -1.4437906434664151e+05, + ), + vel_expct_km_s: Vector3::new( + 9.6061748763649357e-01, + -2.0361038623448113e-01, + -1.8380826294069577e-01, + ), + }, + AberrationCase { + correction: Aberration::XLT, + pos_expct_km: Vector3::new( + -8.1601439447537065e+04, + -3.4550204350015521e+05, + -1.4440340782643855e+05, + ), + vel_expct_km_s: Vector3::new( + 9.6071525662101465e-01, + -2.0358827342129260e-01, + -1.8380776693460277e-01, + ), + }, + AberrationCase { + correction: Aberration::XLT_S, + pos_expct_km: Vector3::new( + -8.1582459098574356e+04, + -3.4550600420432026e+05, + -1.4440465574480488e+05, + ), + vel_expct_km_s: Vector3::new( + 9.6080620884495171e-01, + -2.0355606455727215e-01, + -1.8380276724235226e-01, + ), + }, + AberrationCase { + correction: Aberration::XCN, + pos_expct_km: Vector3::new( + -8.1601441613525152e+04, + -3.4550204579737782e+05, + -1.4440340883307904e+05, + ), + vel_expct_km_s: Vector3::new( + 9.6071525719129625e-01, + -2.0358827357191700e-01, + -1.8380776700407786e-01, + ), + }, + AberrationCase { + correction: Aberration::XCN_S, + pos_expct_km: Vector3::new( + -8.1582461264569836e+04, + -3.4550600650161679e+05, + -1.4440465675147722e+05, + ), + vel_expct_km_s: Vector3::new( + 9.6080620941528405e-01, + -2.0355606470851764e-01, + -1.8380276731210626e-01, + ), + }, + ]; + + for (cno, case) in cases.iter().enumerate() { + let state = ctx + .translate( + LUNA_J2000, + EARTH_MOON_BARYCENTER_J2000, + epoch, + case.correction, + ) + .unwrap(); + + let pos_km = state.radius_km; + let vel_km_s = state.velocity_km_s; + + println!("{state}"); + + // We expect exactly the same output as SPICE to machine precision. + assert!( + relative_eq!(pos_km, case.pos_expct_km, epsilon = EPSILON), + "got {} but want {} with {} (#{cno}) => err = {:.3e} km", + pos_km, + case.pos_expct_km, + case.correction.unwrap(), + (pos_km - case.pos_expct_km).norm() + ); + + assert!( + relative_eq!( + vel_km_s, + case.vel_expct_km_s, + epsilon = ABERRATION_VELOCITY_EPSILON_KM_S + ), + "got {} but want {} with {} (#{cno}) => err = {:.3e} km/s", + vel_km_s, + case.vel_expct_km_s, + case.correction.unwrap(), + (vel_km_s - case.vel_expct_km_s).norm() + ); + + println!( + "got {} but want {} with {} (#{cno}) => err = {:.3e} km/s", + vel_km_s, + case.vel_expct_km_s, + case.correction.unwrap(), + (vel_km_s - case.vel_expct_km_s).norm() + ); + } +} diff --git a/anise/tests/ephemerides/validation/compare.rs b/anise/tests/ephemerides/validation/compare.rs index a78cb730..04ed7988 100644 --- a/anise/tests/ephemerides/validation/compare.rs +++ b/anise/tests/ephemerides/validation/compare.rs @@ -71,6 +71,7 @@ pub struct CompareEphem { pub output_file_name: String, pub num_queries_per_pair: usize, pub dry_run: bool, + pub aberration: Option, pub writer: ArrowWriter, pub batch_src_frame: Vec, pub batch_dst_frame: Vec, @@ -86,6 +87,7 @@ impl CompareEphem { input_file_names: Vec, output_file_name: String, num_queries_per_pair: usize, + aberration: Option, ) -> Self { let _ = pretty_env_logger::try_init(); @@ -110,6 +112,7 @@ impl CompareEphem { output_file_name, input_file_names, num_queries_per_pair, + aberration, writer, dry_run: false, batch_src_frame: Vec::new(), @@ -137,6 +140,14 @@ impl CompareEphem { spice::furnsh(path); } + // If there is a light time correction, start after the epoch because the light time correction + // will cause us to seek out of the definition bounds. + + let bound_offset = match self.aberration { + None => 0.0_f64.seconds(), + Some(_) => 36.0_f64.hours(), + }; + // Build the pairs of the SPICE and ANISE queries at the same time as we create those instances. let mut pairs: HashMap<(i32, i32), (Frame, Frame, Epoch, Epoch)> = HashMap::new(); @@ -173,9 +184,9 @@ impl CompareEphem { let to_frame = Frame::from_ephem_j2000(ephem2.target_id); // Query the ephemeris data for a bunch of different times. - let start_epoch = ephem1.start_epoch().max(ephem2.start_epoch()); + let start_epoch = ephem1.start_epoch().max(ephem2.start_epoch()) + bound_offset; - let end_epoch = ephem1.end_epoch().min(ephem2.end_epoch()); + let end_epoch = ephem1.end_epoch().min(ephem2.end_epoch()) - bound_offset; pairs.insert(key, (from_frame, to_frame, start_epoch, end_epoch)); } @@ -228,7 +239,7 @@ impl CompareEphem { } for epoch in time_it { - let data = match ctx.translate_from_to_geometric(*from_frame, *to_frame, epoch) { + let data = match ctx.translate(*from_frame, *to_frame, epoch, self.aberration) { Ok(state) => { // Find the SPICE names let targ = @@ -246,8 +257,18 @@ impl CompareEphem { }; // Perform the same query in SPICE - let (spice_state, _) = - spice::spkezr(&targ, epoch.to_et_seconds(), "J2000", "NONE", &obs); + let spice_ab_corr = match self.aberration { + None => "NONE".to_string(), + Some(corr) => format!("{corr:?}"), + }; + + let (spice_state, _) = spice::spkezr( + &targ, + epoch.to_et_seconds(), + "J2000", + &spice_ab_corr, + &obs, + ); EphemValData { src_frame: format!("{from_frame:e}"), diff --git a/anise/tests/ephemerides/validation/type02_chebyshev_jpl_de.rs b/anise/tests/ephemerides/validation/type02_chebyshev_jpl_de.rs index 9bf60d57..0d89df4e 100644 --- a/anise/tests/ephemerides/validation/type02_chebyshev_jpl_de.rs +++ b/anise/tests/ephemerides/validation/type02_chebyshev_jpl_de.rs @@ -9,6 +9,7 @@ */ use super::{compare::*, validate::Validation}; +use anise::prelude::Aberration; #[ignore = "Requires Rust SPICE -- must be executed serially"] #[test] @@ -18,6 +19,7 @@ fn validate_jplde_de440_full() { vec!["../data/de440.bsp".to_string()], file_name.clone(), 1_000, + None, ); let err_count = comparator.run(); @@ -34,12 +36,13 @@ fn validate_jplde_de440_full() { #[ignore = "Requires Rust SPICE -- must be executed serially"] #[test] -fn validate_jplde_de440s() { +fn validate_jplde_de440s_no_aberration() { let output_file_name = "spk-type2-validation-de440s".to_string(); let comparator = CompareEphem::new( vec!["../data/de440s.bsp".to_string()], output_file_name.clone(), 1_000, + None, ); let err_count = comparator.run(); @@ -53,3 +56,29 @@ fn validate_jplde_de440s() { validator.validate(); } + +#[ignore = "Requires Rust SPICE -- must be executed serially"] +#[test] +fn validate_jplde_de440s_aberration_lt() { + let output_file_name = "spk-type2-validation-de440s-lt-aberration".to_string(); + let comparator = CompareEphem::new( + vec!["../data/de440s.bsp".to_string()], + output_file_name.clone(), + 1_000, + Aberration::LT, + ); + + let err_count = comparator.run(); + + assert_eq!(err_count, 10, "A few are expected to fail"); + + let validator = Validation { + file_name: output_file_name, + max_q75_err: 1e-3, + max_q99_err: 5e-3, + max_abs_err: 0.09, + ..Default::default() + }; + + validator.validate(); +} diff --git a/anise/tests/ephemerides/validation/type13_hermite.rs b/anise/tests/ephemerides/validation/type13_hermite.rs index ed27787d..ca7f17ce 100644 --- a/anise/tests/ephemerides/validation/type13_hermite.rs +++ b/anise/tests/ephemerides/validation/type13_hermite.rs @@ -18,6 +18,7 @@ fn validate_hermite_type13_from_gmat() { vec!["../data/gmat-hermite.bsp".to_string()], file_name.clone(), 10_000, + None, ); let err_count = comparator.run(); @@ -41,6 +42,7 @@ fn validate_hermite_type13_with_varying_segment_sizes() { vec!["../data/variable-seg-size-hermite.bsp".to_string()], file_name.clone(), 10_000, + None, ); let err_count = comparator.run();