orbit: add layer facade that combines propagation and elements

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"

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

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};

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> {
+        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<Self, PropagationError> {
+        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<Self, PropagationError> {
+        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<State, PropagationError> {
+        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<State, PropagationError> {
+        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<KeplerianElements, PropagationError> {
+        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<f64, PropagationError> {
+        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<Vec3, PropagationError> {
+        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
+    }
+}

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());
+}