From da86f2f2ede444a6834dd1f4dd2a077c6b66aade Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 14:21:55 +0800 Subject: [PATCH] orbit: add layer facade that combines propagation and elements --- Cargo.lock | 8 ++ crates/orbit/apsis-orbit/Cargo.toml | 6 + crates/orbit/apsis-orbit/src/lib.rs | 53 ++++++++ crates/orbit/apsis-orbit/src/orbit.rs | 151 +++++++++++++++++++++ crates/orbit/apsis-orbit/tests/orbit.rs | 172 ++++++++++++++++++++++++ 5 files changed, 390 insertions(+) create mode 100644 crates/orbit/apsis-orbit/src/orbit.rs create mode 100644 crates/orbit/apsis-orbit/tests/orbit.rs diff --git a/Cargo.lock b/Cargo.lock index 78b9d75..274b747 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -121,6 +121,14 @@ version = "0.1.0" [[package]] name = "apsis-orbit" version = "0.1.0" +dependencies = [ + "apsis-core", + "apsis-kepler", + "apsis-math", + "apsis-sgp4", + "apsis-time", + "apsis-tle", +] [[package]] name = "apsis-output" diff --git a/crates/orbit/apsis-orbit/Cargo.toml b/crates/orbit/apsis-orbit/Cargo.toml index a221472..070365e 100644 --- a/crates/orbit/apsis-orbit/Cargo.toml +++ b/crates/orbit/apsis-orbit/Cargo.toml @@ -11,3 +11,9 @@ publish = false workspace = true [dependencies] +apsis-core.workspace = true +apsis-math.workspace = true +apsis-time.workspace = true +apsis-tle.workspace = true +apsis-sgp4.workspace = true +apsis-kepler.workspace = true diff --git a/crates/orbit/apsis-orbit/src/lib.rs b/crates/orbit/apsis-orbit/src/lib.rs index 1496f63..decaa4c 100644 --- a/crates/orbit/apsis-orbit/src/lib.rs +++ b/crates/orbit/apsis-orbit/src/lib.rs @@ -1 +1,54 @@ //! Orbital state facade that integrates the orbit layer. +//! +//! This crate is the single entry point to the orbit layer. It re-exports +//! the element parsing, the analytical propagator, and the two-body model, +//! and adds one type that ties them together. +//! +//! [`Orbit`] holds a satellite prepared for propagation together with the +//! gravity model it was initialised against, and offers the conversions +//! between the two ways of describing an orbit that the layer supports: +//! +//! - Mean elements propagated by SGP4, which is the model a two-line +//! element set is written for and the only way to reproduce the published +//! positions of a catalogued satellite. +//! - Instantaneous Keplerian elements, which describe the osculating +//! two-body orbit at one instant and are useful for inspection and for +//! planning a maneuver. +//! +//! The two descriptions are not interchangeable. Propagating mean elements +//! with the two-body model, or osculating elements with SGP4, gives a +//! different and wrong answer; the conversions here are explicit about which +//! direction they go. +//! +//! ``` +//! use apsis_orbit::{GravityModel, Orbit}; +//! use apsis_tle::TwoLineElement; +//! +//! let text = concat!( +//! "1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927\r\n", +//! "2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537\r\n", +//! ); +//! let element = TwoLineElement::parse(text).expect("sample is well formed"); +//! let orbit = Orbit::from_element(&element, GravityModel::WGS72) +//! .expect("sample elements are propagatable"); +//! +//! let state = orbit.propagate_minutes(60.0).expect("state exists"); +//! assert!(state.speed() > 7.0); +//! +//! // The osculating elements at that instant. +//! let osculating = orbit.osculating_elements(60.0).expect("state exists"); +//! assert!(osculating.eccentricity < 0.01); +//! ``` + +#![forbid(unsafe_code)] +#![deny(missing_docs)] + +mod orbit; + +pub use orbit::Orbit; + +pub use apsis_kepler::{GravityParameter, KeplerianElements}; +pub use apsis_sgp4::{ + Elements, GravityModel, PropagationError, PropagationErrorKind, Satellite, State, +}; +pub use apsis_tle::{Classification, ElementError, ElementErrorKind, TwoLineElement}; diff --git a/crates/orbit/apsis-orbit/src/orbit.rs b/crates/orbit/apsis-orbit/src/orbit.rs new file mode 100644 index 0000000..b915e82 --- /dev/null +++ b/crates/orbit/apsis-orbit/src/orbit.rs @@ -0,0 +1,151 @@ +//! The combined orbit type. + +use apsis_kepler::{GravityParameter, KeplerianElements}; +use apsis_math::Vec3; +use apsis_sgp4::{Elements, GravityModel, PropagationError, Satellite, State}; +use apsis_tle::TwoLineElement; + +/// A satellite orbit, held as mean elements ready for propagation. +/// +/// The type pairs a prepared SGP4 satellite with the gravity model it was +/// initialised against, because the osculating element conversion needs the +/// same gravitational parameter to be consistent with the propagation. +#[derive(Debug, Clone)] +pub struct Orbit { + satellite: Satellite, + parameter: GravityParameter, +} + +impl Orbit { + /// Prepares an orbit from a parsed element set and a gravity model. + /// + /// The gravity model must match the one the element set was fitted + /// against, which for a catalogued satellite means + /// [`GravityModel::WGS72`]. Returns `None` when the elements are not + /// propagatable. + pub fn from_element(element: &TwoLineElement, model: GravityModel) -> Option { + Self::try_from_element(element, model).ok() + } + + /// Prepares an orbit from a parsed element set and a gravity model. + /// + /// Returns an error describing which element is not propagatable. + pub fn try_from_element( + element: &TwoLineElement, + model: GravityModel, + ) -> Result { + let satellite = Satellite::try_from_element(element, model)?; + Ok(Self { + satellite, + parameter: parameter_for(model), + }) + } + + /// Prepares an orbit from mean elements already in model units. + pub fn from_elements( + elements: &Elements, + model: GravityModel, + ) -> Result { + let satellite = Satellite::try_from_elements(elements, model)?; + Ok(Self { + satellite, + parameter: parameter_for(model), + }) + } + + /// Returns the underlying propagator. + pub fn satellite(&self) -> &Satellite { + &self.satellite + } + + /// Returns the gravitational parameter paired with the gravity model. + /// + /// The value is the one used for the osculating element conversion, so + /// that the conversion is consistent with the propagation. + pub fn gravity_parameter(&self) -> GravityParameter { + self.parameter + } + + /// Returns the epoch of the elements as a Modified Julian Date on the + /// UTC scale. + /// + /// The epoch is the value written in the element set, which is the + /// origin the propagator measures time from. + pub fn epoch_mjd(&self) -> f64 { + self.satellite.epoch_mjd() + } + + /// Returns `true` when the propagator uses the deep space branch. + pub fn is_deep_space(&self) -> bool { + self.satellite.is_deep_space() + } + + /// Returns the state at a number of minutes from the element epoch. + /// + /// Position is in kilometers and velocity in kilometers per second, in + /// the True Equator Mean Equinox frame. + pub fn propagate_minutes(&self, minutes: f64) -> Result { + self.satellite.propagate_minutes(minutes) + } + + /// Returns the state at a Modified Julian Date on the UTC scale. + pub fn propagate_at(&self, utc_mjd: f64) -> Result { + self.propagate_minutes((utc_mjd - self.satellite.epoch_mjd()) * 1_440.0) + } + + /// Returns the osculating Keplerian elements at a number of minutes from + /// the element epoch. + /// + /// The osculating elements describe the two-body orbit that matches the + /// position and velocity at that instant. They are not the mean elements + /// of the element set, and they change along the orbit because SGP4 + /// includes perturbations that the two-body model does not. + /// + /// Returns an error when the state cannot be propagated, or when the + /// instantaneous orbit is not elliptical, which can happen for a state + /// the propagation has pushed onto an escape trajectory. + pub fn osculating_elements(&self, minutes: f64) -> Result { + let state = self.propagate_minutes(minutes)?; + // The two crates define their own state type, so the components are + // transferred explicitly rather than relying on a shared type. + let kepler_state = apsis_kepler::State { + position: state.position, + velocity: state.velocity, + }; + KeplerianElements::from_state(&kepler_state, self.parameter).ok_or_else(|| { + PropagationError::new(apsis_sgp4::PropagationErrorKind::OutOfRange) + .with_detail("the instantaneous orbit is not elliptical") + }) + } + + /// Returns the magnitude of the state vector. + /// + /// This is a convenience for the common case of wanting only the + /// distance from the center of the Earth. + pub fn radius_at(&self, minutes: f64) -> Result { + self.propagate_minutes(minutes) + .map(|state| state.position.norm()) + } + + /// Returns the position at a number of minutes from the element epoch. + pub fn position_at(&self, minutes: f64) -> Result { + self.propagate_minutes(minutes).map(|state| state.position) + } +} + +/// Returns the gravitational parameter that matches a gravity model. +/// +/// The two models differ slightly in their value of the gravitational +/// parameter, and using the wrong one shifts the osculating semi-major axis +/// by a small amount. The values are taken from the model definitions rather +/// than duplicated, so a change to a model propagates here. +fn parameter_for(model: GravityModel) -> GravityParameter { + // The two parameters differ by 0.3 km^3/s^2, about one part in a + // million. The comparison is on the model's defining value, which is + // exact for the constants defined in this workspace. + if model.mu == GravityModel::WGS84.mu { + GravityParameter::WGS84 + } else { + GravityParameter::WGS72 + } +} diff --git a/crates/orbit/apsis-orbit/tests/orbit.rs b/crates/orbit/apsis-orbit/tests/orbit.rs new file mode 100644 index 0000000..06fc2bb --- /dev/null +++ b/crates/orbit/apsis-orbit/tests/orbit.rs @@ -0,0 +1,172 @@ +//! Tests for the combined orbit type. + +use apsis_orbit::{GravityModel, Orbit}; +use apsis_tle::TwoLineElement; + +const ISS: &str = concat!( + "1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927\r\n", + "2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537\r\n", +); + +const DEEP_SPACE: &str = concat!( + "1 09880U 77021A 06174.85710247 .00000090 00000-0 10000-3 0 7322\r\n", + "2 09880 64.1526 313.1296 6546925 103.4683 305.3980 2.00754276 13047\r\n", +); + +fn orbit(text: &str) -> Orbit { + let element = TwoLineElement::parse(text).expect("test element set parses"); + Orbit::from_element(&element, GravityModel::WGS72).expect("test element set is propagatable") +} + +#[test] +fn construction_succeeds_for_a_valid_element_set() { + let orbit = orbit(ISS); + assert!(!orbit.is_deep_space()); + assert_eq!( + orbit.gravity_parameter(), + apsis_orbit::GravityParameter::WGS72 + ); +} + +#[test] +fn deep_space_orbits_are_reported() { + assert!(orbit(DEEP_SPACE).is_deep_space()); +} + +#[test] +fn propagation_matches_the_underlying_propagator() { + let orbit = orbit(ISS); + let via_orbit = orbit.propagate_minutes(360.0).expect("state exists"); + let via_satellite = orbit + .satellite() + .propagate_minutes(360.0) + .expect("state exists"); + assert_eq!(via_orbit, via_satellite); +} + +#[test] +fn position_and_radius_accessors_agree_with_the_state() { + let orbit = orbit(ISS); + let state = orbit.propagate_minutes(120.0).expect("state exists"); + assert_eq!( + orbit.position_at(120.0).expect("state exists"), + state.position + ); + assert!((orbit.radius_at(120.0).expect("state exists") - state.position.norm()).abs() < 1e-12); +} + +#[test] +fn propagate_at_uses_the_element_epoch() { + let orbit = orbit(ISS); + let epoch = orbit.epoch_mjd(); + let at_epoch = orbit.propagate_at(epoch).expect("state exists"); + let at_zero = orbit.propagate_minutes(0.0).expect("state exists"); + assert!((at_epoch.position - at_zero.position).norm() < 1e-6); + + // One day later is 1440 minutes. + let next_day = orbit.propagate_at(epoch + 1.0).expect("state exists"); + let one_day = orbit.propagate_minutes(1440.0).expect("state exists"); + assert!((next_day.position - one_day.position).norm() < 1e-3); +} + +#[test] +fn osculating_elements_describe_a_low_earth_orbit() { + let orbit = orbit(ISS); + let elements = orbit.osculating_elements(0.0).expect("state exists"); + // These values were checked against an independent computation from the + // same state vector. The inclination agrees with the catalogued value + // closely, because the oblateness has little effect on it, while the + // eccentricity differs from the catalogued 0.0006703. That difference + // is expected: the element set stores mean elements and the osculating + // eccentricity includes the short period variation the two-body model + // does not reproduce. + assert!( + (elements.semi_major_axis - 6_725.542).abs() < 0.01, + "semi-major axis {}", + elements.semi_major_axis + ); + assert!( + (elements.eccentricity - 0.000_832_1).abs() < 1e-6, + "eccentricity {}", + elements.eccentricity + ); + assert!( + (elements.inclination - 51.6217_f64.to_radians()).abs() < 1e-4, + "inclination {}", + elements.inclination.to_degrees() + ); +} + +#[test] +fn osculating_elements_reproduce_the_state() { + let orbit = orbit(ISS); + let state = orbit.propagate_minutes(45.0).expect("state exists"); + let elements = orbit.osculating_elements(45.0).expect("state exists"); + let rebuilt = elements + .to_state(orbit.gravity_parameter()) + .expect("elements are valid"); + assert!((rebuilt.position - state.position).norm() < 1e-6); + assert!((rebuilt.velocity - state.velocity).norm() < 1e-9); +} + +#[test] +fn osculating_elements_change_along_the_orbit() { + // The osculating elements are not the mean elements, so they vary + // because the perturbations SGP4 models are absent from the two-body + // reconstruction. + let orbit = orbit(ISS); + let start = orbit.osculating_elements(0.0).expect("state exists"); + let later = orbit.osculating_elements(720.0).expect("state exists"); + assert!((start.eccentricity - later.eccentricity).abs() > 1e-7); +} + +#[test] +fn osculating_elements_work_for_a_deep_space_orbit() { + let orbit = orbit(DEEP_SPACE); + let elements = orbit.osculating_elements(0.0).expect("state exists"); + // The catalogued mean values are an eccentricity of 0.6546925 and an + // inclination of 64.1526 degrees. The osculating values are close but + // not identical for the reason given in the low Earth orbit test. + assert!( + (elements.eccentricity - 0.655_871_3).abs() < 1e-6, + "eccentricity {}", + elements.eccentricity + ); + assert!( + (elements.inclination - 64.1696_f64.to_radians()).abs() < 1e-4, + "inclination {}", + elements.inclination.to_degrees() + ); + assert!( + (elements.semi_major_axis - 26_543.982).abs() < 0.01, + "semi-major axis {}", + elements.semi_major_axis + ); +} + +#[test] +fn the_gravity_parameter_matches_the_model() { + let element = TwoLineElement::parse(ISS).expect("parses"); + let wgs72 = Orbit::from_element(&element, GravityModel::WGS72).expect("propagatable"); + assert_eq!( + wgs72.gravity_parameter(), + apsis_orbit::GravityParameter::WGS72 + ); + + let wgs84 = Orbit::from_element(&element, GravityModel::WGS84).expect("propagatable"); + assert_eq!( + wgs84.gravity_parameter(), + apsis_orbit::GravityParameter::WGS84 + ); +} + +#[test] +fn the_facade_reexports_the_layer() { + // The facade must make the element parser and the propagator types + // reachable without a separate dependency. + let element = apsis_orbit::TwoLineElement::parse(ISS).expect("parses"); + assert_eq!(element.satellite_number, 25544); + let satellite = + apsis_orbit::Satellite::from_element(&element, GravityModel::WGS72).expect("propagatable"); + assert!(satellite.propagate_minutes(0.0).is_ok()); +}