From c256dc0da2894b686af1f7536a78c5e3dfb62389 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 23:17:16 +0800 Subject: [PATCH] orbit: add maneuver models and wire the orbit facade --- Cargo.lock | 7 + crates/orbit/apsis-maneuver/Cargo.toml | 3 + crates/orbit/apsis-maneuver/src/burn.rs | 279 ++++++++++++ crates/orbit/apsis-maneuver/src/error.rs | 76 ++++ crates/orbit/apsis-maneuver/src/hohmann.rs | 135 ++++++ crates/orbit/apsis-maneuver/src/impulse.rs | 99 +++++ crates/orbit/apsis-maneuver/src/lib.rs | 65 +++ crates/orbit/apsis-maneuver/tests/maneuver.rs | 420 ++++++++++++++++++ crates/orbit/apsis-orbit/Cargo.toml | 2 + crates/orbit/apsis-orbit/src/lib.rs | 17 +- crates/orbit/apsis-orbit/src/orbit.rs | 136 ++++++ .../apsis-orbit/tests/ephemeris_maneuver.rs | 188 ++++++++ 12 files changed, 1426 insertions(+), 1 deletion(-) create mode 100644 crates/orbit/apsis-maneuver/src/burn.rs create mode 100644 crates/orbit/apsis-maneuver/src/error.rs create mode 100644 crates/orbit/apsis-maneuver/src/hohmann.rs create mode 100644 crates/orbit/apsis-maneuver/src/impulse.rs create mode 100644 crates/orbit/apsis-maneuver/tests/maneuver.rs create mode 100644 crates/orbit/apsis-orbit/tests/ephemeris_maneuver.rs diff --git a/Cargo.lock b/Cargo.lock index ec0eac4..3dee232 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -166,6 +166,11 @@ version = "0.1.0" [[package]] name = "apsis-maneuver" version = "0.1.0" +dependencies = [ + "apsis-core", + "apsis-kepler", + "apsis-math", +] [[package]] name = "apsis-math" @@ -193,7 +198,9 @@ name = "apsis-orbit" version = "0.1.0" dependencies = [ "apsis-core", + "apsis-ephemeris", "apsis-kepler", + "apsis-maneuver", "apsis-math", "apsis-sgp4", "apsis-time", diff --git a/crates/orbit/apsis-maneuver/Cargo.toml b/crates/orbit/apsis-maneuver/Cargo.toml index 032cdef..e1bd7bf 100644 --- a/crates/orbit/apsis-maneuver/Cargo.toml +++ b/crates/orbit/apsis-maneuver/Cargo.toml @@ -11,3 +11,6 @@ publish = false workspace = true [dependencies] +apsis-core.workspace = true +apsis-math.workspace = true +apsis-kepler.workspace = true diff --git a/crates/orbit/apsis-maneuver/src/burn.rs b/crates/orbit/apsis-maneuver/src/burn.rs new file mode 100644 index 0000000..43eacae --- /dev/null +++ b/crates/orbit/apsis-maneuver/src/burn.rs @@ -0,0 +1,279 @@ +//! The finite burn. + +use crate::error::{ManeuverError, ManeuverErrorKind}; +use apsis_kepler::{GravityParameter, KeplerianElements, State}; +use apsis_math::Vec3; + +/// A maneuver that applies a constant acceleration for a duration. +/// +/// A finite burn is the model for a thruster that is on for long enough that +/// its duration cannot be ignored: an electric thruster, whose burn may last +/// a large fraction of an orbit, or any burn where the spacecraft moves +/// appreciably while thrusting. +/// +/// The acceleration is held constant in an inertial frame chosen by the +/// caller. That is the simplification this model makes: a real thruster +/// fires along a body axis that rotates as the spacecraft turns to hold a +/// direction, and modelling that requires the attitude history, which is not +/// part of this crate. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct FiniteBurn { + /// The acceleration vector, in kilometers per second squared. + acceleration: Vec3, + /// The duration, in seconds. + duration: f64, +} + +impl FiniteBurn { + /// Creates a burn from an acceleration and a duration. + /// + /// Returns an error when a value is not finite, when the duration is not + /// positive, or when the acceleration is zero. + pub fn new(acceleration: Vec3, duration: f64) -> Result { + if !acceleration.is_finite() || !duration.is_finite() { + return Err(ManeuverError::new(ManeuverErrorKind::NotFinite)); + } + if duration <= 0.0 { + return Err(ManeuverError::new(ManeuverErrorKind::NonPositiveDuration) + .with_detail(format!("duration {duration}"))); + } + if acceleration.norm() == 0.0 { + return Err( + ManeuverError::new(ManeuverErrorKind::NonPositiveAcceleration) + .with_detail("the acceleration is zero"), + ); + } + Ok(Self { + acceleration, + duration, + }) + } + + /// Creates a burn of a given acceleration magnitude along a direction. + /// + /// The direction is normalized, so the acceleration magnitude is fixed + /// whatever the length of the direction vector. + pub fn along(direction: Vec3, acceleration: f64, duration: f64) -> Result { + if !acceleration.is_finite() || acceleration <= 0.0 { + return Err( + ManeuverError::new(ManeuverErrorKind::NonPositiveAcceleration) + .with_detail(format!("acceleration {acceleration}")), + ); + } + let unit = direction.normalize().ok_or_else(|| { + ManeuverError::new(ManeuverErrorKind::InvalidState) + .with_detail("the direction is zero or not finite") + })?; + Self::new(unit * acceleration, duration) + } + + /// Returns the acceleration vector, in kilometers per second squared. + pub const fn acceleration(self) -> Vec3 { + self.acceleration + } + + /// Returns the magnitude of the acceleration, in kilometers per second + /// squared. + pub fn acceleration_magnitude(self) -> f64 { + self.acceleration.norm() + } + + /// Returns the duration, in seconds. + pub const fn duration(self) -> f64 { + self.duration + } + + /// Returns the velocity change the acceleration would give in free space, + /// in kilometers per second. + /// + /// This is the acceleration times the duration, the delta-v the burn + /// would deliver without gravity acting during it. Under gravity the + /// achieved change is close to this for a short burn and differs from it + /// for a long one, which is the difference the finite model exists to + /// capture. + pub fn nominal_delta_v(self) -> f64 { + self.acceleration.norm() * self.duration + } + + /// Returns the impulse equivalent to the nominal change. + /// + /// This is the impulsive burn the finite one approximates as its + /// duration shrinks, and it is what [`FiniteBurn::apply`] returns in the + /// limit of a constant gravity field. + pub fn equivalent_impulse(self) -> Result { + crate::impulse::Impulse::new(self.acceleration * self.duration) + } + + /// Applies the burn to a state under two-body gravity. + /// + /// The position and velocity are integrated forward over the duration + /// with the constant acceleration added to gravity. The exact solution + /// of constant acceleration plus point mass gravity has no closed form, + /// so the integration is a fixed step fourth order Runge-Kutta, which is + /// accurate for the short arcs a burn occupies. + /// + /// The achieved velocity change is measured against a coast of the same + /// duration, not against the initial velocity. That distinction matters: + /// gravity changes the velocity during the burn even with no thrust, so + /// the difference between the final and the initial velocity includes + /// the natural rotation of the velocity vector. Comparing against the + /// coast isolates the contribution of the thrust, which is what the + /// burn bought; it also makes the result converge on the nominal change + /// as the duration shrinks, which the difference from the initial + /// velocity would not. + /// + /// Returns the state at the end of the burn together with the achieved + /// velocity change. + pub fn apply( + self, + state: State, + parameter: GravityParameter, + steps: usize, + ) -> Result { + if !state.position.is_finite() || !state.velocity.is_finite() { + return Err(ManeuverError::new(ManeuverErrorKind::InvalidState)); + } + if steps == 0 { + return Err(ManeuverError::new(ManeuverErrorKind::InvalidState) + .with_detail("the step count is zero")); + } + let mu = parameter.value(); + let step = self.duration / steps as f64; + + // The two-body acceleration at a position. + let gravity = |position: Vec3| -> Vec3 { + let radius = position.norm(); + if radius == 0.0 { + return Vec3::ZERO; + } + position * (-mu / (radius * radius * radius)) + }; + // The thrust added to gravity, which is what the burn integrates. + let thrust_and_gravity = |position: Vec3| -> Vec3 { gravity(position) + self.acceleration }; + + // The burn trajectory. + let mut position = state.position; + let mut velocity = state.velocity; + for _ in 0..steps { + let (new_position, new_velocity) = + rk4_step(position, velocity, step, &thrust_and_gravity); + position = new_position; + velocity = new_velocity; + } + + // The coast trajectory over the same duration, with the same + // integrator and step, so the only difference is the thrust. + let mut coast_position = state.position; + let mut coast_velocity = state.velocity; + for _ in 0..steps { + let (new_position, new_velocity) = + rk4_step(coast_position, coast_velocity, step, &gravity); + coast_position = new_position; + coast_velocity = new_velocity; + } + + if !position.is_finite() || !velocity.is_finite() { + return Err(ManeuverError::new(ManeuverErrorKind::InvalidState) + .with_detail("the integration diverged")); + } + + Ok(BurnOutcome { + state: State { position, velocity }, + achieved_delta_v: velocity - coast_velocity, + nominal_delta_v: self.nominal_delta_v(), + duration: self.duration, + }) + } + + /// Applies the burn to a set of Keplerian elements. + /// + /// The elements are converted to a state, the burn is applied to it, and + /// the result is converted back. Returns an error when the elements are + /// not valid, when the integration fails, or when the resulting state is + /// no longer an elliptical orbit, which happens when the burn is large + /// enough to escape. + pub fn apply_to_elements( + self, + elements: KeplerianElements, + parameter: GravityParameter, + steps: usize, + ) -> Result { + let state = elements.to_state(parameter).ok_or_else(|| { + ManeuverError::new(ManeuverErrorKind::InvalidState) + .with_detail("the elements are not a valid elliptical orbit") + })?; + self.apply(state, parameter, steps) + } +} + +/// What a finite burn produced. +/// +/// The outcome keeps the final state, the change the burn achieved, and the +/// change the acceleration would have given in free space, so a caller can +/// see how much gravity cost or added during the burn. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct BurnOutcome { + /// The state at the end of the burn. + pub state: State, + /// The velocity change actually achieved, in kilometers per second. + pub achieved_delta_v: Vec3, + /// The velocity change the acceleration would give in free space, in + /// kilometers per second. + pub nominal_delta_v: f64, + /// The burn duration, in seconds. + pub duration: f64, +} + +impl BurnOutcome { + /// Returns the magnitude of the achieved velocity change. + pub fn achieved_magnitude(&self) -> f64 { + self.achieved_delta_v.norm() + } + + /// Returns the difference between the achieved and nominal change, in + /// kilometers per second. + /// + /// A negative value means gravity turned the velocity against the thrust + /// during the burn, so the burn bought less than its free space delta-v. + pub fn gravity_loss(&self) -> f64 { + self.achieved_magnitude() - self.nominal_delta_v + } + + /// Returns the elements at the end of the burn. + /// + /// Returns `None` when the final state is not an elliptical orbit, which + /// happens when the burn is large enough to escape. + pub fn elements(&self, parameter: GravityParameter) -> Option { + KeplerianElements::from_state(&self.state, parameter) + } +} + +/// One classical Runge-Kutta fourth order step of the two-body equations +/// with a constant acceleration added. +/// +/// The state is `(position, velocity)` and the equations are the usual +/// `dr/dt = v` and `dv/dt = a(r)`, with `a` combining thrust and gravity. +/// The four stages are the standard ones; writing them out keeps the loop +/// allocation free. +fn rk4_step( + position: Vec3, + velocity: Vec3, + step: f64, + acceleration: &impl Fn(Vec3) -> Vec3, +) -> (Vec3, Vec3) { + let k1_v = acceleration(position); + let k1_r = velocity; + + let k2_v = acceleration(position + k1_r * (step / 2.0)); + let k2_r = velocity + k1_v * (step / 2.0); + + let k3_v = acceleration(position + k2_r * (step / 2.0)); + let k3_r = velocity + k2_v * (step / 2.0); + + let k4_v = acceleration(position + k3_r * step); + let k4_r = velocity + k3_v * step; + + let new_position = position + (k1_r + k2_r * 2.0 + k3_r * 2.0 + k4_r) * (step / 6.0); + let new_velocity = velocity + (k1_v + k2_v * 2.0 + k3_v * 2.0 + k4_v) * (step / 6.0); + (new_position, new_velocity) +} diff --git a/crates/orbit/apsis-maneuver/src/error.rs b/crates/orbit/apsis-maneuver/src/error.rs new file mode 100644 index 0000000..f83ccbd --- /dev/null +++ b/crates/orbit/apsis-maneuver/src/error.rs @@ -0,0 +1,76 @@ +//! Errors reported while building or applying a maneuver. + +use core::fmt; + +/// The category of a maneuver error. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +#[non_exhaustive] +pub enum ManeuverErrorKind { + /// A velocity change, acceleration, or duration is not finite. + NotFinite, + /// The duration of a burn is not positive. + NonPositiveDuration, + /// The acceleration of a burn is not positive. + NonPositiveAcceleration, + /// A radius is not finite or not positive. + InvalidRadius, + /// The two radii of a transfer are equal, so there is nothing to do. + DegenerateTransfer, + /// A state is missing the data the operation needs. + InvalidState, +} + +impl fmt::Display for ManeuverErrorKind { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + let text = match self { + Self::NotFinite => "a value is not finite", + Self::NonPositiveDuration => "the duration is not positive", + Self::NonPositiveAcceleration => "the acceleration is not positive", + Self::InvalidRadius => "the radius is not usable", + Self::DegenerateTransfer => "the two radii are equal", + Self::InvalidState => "the state is not usable", + }; + f.write_str(text) + } +} + +/// A failure while building or applying a maneuver. +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct ManeuverError { + kind: ManeuverErrorKind, + detail: Option, +} + +impl ManeuverError { + /// Creates an error with the given category. + pub fn new(kind: ManeuverErrorKind) -> Self { + Self { kind, detail: None } + } + + /// Attaches a description of the value that caused the failure. + pub fn with_detail(mut self, detail: impl Into) -> Self { + self.detail = Some(detail.into()); + self + } + + /// Returns the category of this failure. + pub const fn kind(&self) -> ManeuverErrorKind { + self.kind + } + + /// Returns the detail text attached to this error, if any. + pub fn detail(&self) -> Option<&str> { + self.detail.as_deref() + } +} + +impl fmt::Display for ManeuverError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match &self.detail { + Some(detail) => write!(f, "{}: {detail}", self.kind), + None => write!(f, "{}", self.kind), + } + } +} + +impl core::error::Error for ManeuverError {} diff --git a/crates/orbit/apsis-maneuver/src/hohmann.rs b/crates/orbit/apsis-maneuver/src/hohmann.rs new file mode 100644 index 0000000..77efeb9 --- /dev/null +++ b/crates/orbit/apsis-maneuver/src/hohmann.rs @@ -0,0 +1,135 @@ +//! The two impulse Hohmann transfer between coplanar circular orbits. + +use crate::error::{ManeuverError, ManeuverErrorKind}; +use apsis_kepler::GravityParameter; + +/// A two impulse transfer between two coplanar circular orbits. +/// +/// The transfer is the classic Hohmann maneuver: a burn onto an ellipse +/// whose perigee touches the lower orbit and whose apogee touches the upper +/// one, a coast for half the ellipse's period, and a second burn onto the +/// destination circle. It is the cheapest two impulse transfer between two +/// coplanar circular orbits in the impulsive model, which is why it is the +/// reference every other transfer is measured against. +/// +/// The burns are signed along the direction of motion: a positive value +/// raises the orbit and a negative one lowers it, so the same type covers a +/// climb and a descent. The total is the sum of the magnitudes, which is the +/// quantity a propellant budget uses. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct HohmannTransfer { + /// The radius of the orbit the transfer starts from, in kilometers. + pub from_radius: f64, + /// The radius of the orbit the transfer arrives at, in kilometers. + pub to_radius: f64, + /// The semi-major axis of the transfer ellipse, in kilometers. + pub transfer_semi_major_axis: f64, + /// The first burn, along the direction of motion, in kilometers per + /// second. + pub first_burn: f64, + /// The second burn, along the direction of motion, in kilometers per + /// second. + pub second_burn: f64, + /// The coast time between the burns, in seconds. + pub transfer_time: f64, +} + +impl HohmannTransfer { + /// Returns the total velocity change, in kilometers per second. + /// + /// This is the sum of the magnitudes of the two burns, which is what the + /// transfer costs whatever its direction. + pub fn total_delta_v(&self) -> f64 { + self.first_burn.abs() + self.second_burn.abs() + } + + /// Returns `true` when the transfer raises the orbit. + pub fn is_raising(&self) -> bool { + self.to_radius > self.from_radius + } + + /// Returns the velocity change of the first burn, unsigned. + pub fn first_burn_magnitude(&self) -> f64 { + self.first_burn.abs() + } + + /// Returns the velocity change of the second burn, unsigned. + pub fn second_burn_magnitude(&self) -> f64 { + self.second_burn.abs() + } + + /// Returns the mean radius the transfer is symmetric about, in + /// kilometers. + /// + /// This is the semi-major axis of the ellipse, which is the arithmetic + /// mean of the two orbit radii. + pub fn mean_radius(&self) -> f64 { + self.transfer_semi_major_axis + } +} + +/// Returns the two impulse transfer between two circular orbits. +/// +/// The radii are in kilometers and must be finite and positive. The two may +/// be given in either order: a larger second radius is a climb and a smaller +/// one a descent. +/// +/// Returns an error when a radius is not usable or when the two are equal, +/// because a transfer between an orbit and itself is not a maneuver and is +/// almost always a mistake in a plan. +/// +/// The burns follow from the vis-viva equation. At the first radius the +/// circular speed is `sqrt(mu / r1)` and the speed on the transfer ellipse +/// is `sqrt(mu (2 / r1 - 1 / a))`, where `a` is the semi-major axis of the +/// ellipse, so the first burn is the difference. The second burn is the +/// difference between the circular speed at the destination and the speed on +/// the ellipse there, taken the other way round because the ellipse arrives +/// slower than the circle it must join. The coast time is half the period of +/// the ellipse. +pub fn hohmann_transfer( + from_radius: f64, + to_radius: f64, + parameter: GravityParameter, +) -> Result { + for radius in [from_radius, to_radius] { + if !radius.is_finite() || radius <= 0.0 { + return Err(ManeuverError::new(ManeuverErrorKind::InvalidRadius) + .with_detail(format!("radius {radius}"))); + } + } + if from_radius == to_radius { + return Err(ManeuverError::new(ManeuverErrorKind::DegenerateTransfer) + .with_detail(format!("both radii are {from_radius}"))); + } + + let mu = parameter.value(); + let transfer_semi_major_axis = (from_radius + to_radius) / 2.0; + + // The circular speeds at each radius. + let circular_from = (mu / from_radius).sqrt(); + let circular_to = (mu / to_radius).sqrt(); + + // The speeds on the transfer ellipse at each radius, from the vis-viva + // equation `v^2 = mu (2 / r - 1 / a)`. + let transfer_from = (mu * (2.0 / from_radius - 1.0 / transfer_semi_major_axis)).sqrt(); + let transfer_to = (mu * (2.0 / to_radius - 1.0 / transfer_semi_major_axis)).sqrt(); + + // The first burn accelerates onto the ellipse; the second accelerates + // from the ellipse speed up to the circular speed at the destination. + // Both are along the direction of motion, and the sign follows from + // which speed is larger. + let first_burn = transfer_from - circular_from; + let second_burn = circular_to - transfer_to; + + // Half the period of the transfer ellipse. + let transfer_time = core::f64::consts::PI * (transfer_semi_major_axis.powi(3) / mu).sqrt(); + + Ok(HohmannTransfer { + from_radius, + to_radius, + transfer_semi_major_axis, + first_burn, + second_burn, + transfer_time, + }) +} diff --git a/crates/orbit/apsis-maneuver/src/impulse.rs b/crates/orbit/apsis-maneuver/src/impulse.rs new file mode 100644 index 0000000..035caf9 --- /dev/null +++ b/crates/orbit/apsis-maneuver/src/impulse.rs @@ -0,0 +1,99 @@ +//! The impulsive maneuver. + +use crate::error::{ManeuverError, ManeuverErrorKind}; +use apsis_math::Vec3; + +/// An instantaneous change in velocity. +/// +/// An impulse is a velocity change applied in zero time, so it changes the +/// velocity of a state but not its position. That is the model for a burn +/// short compared with the orbital period, and it is how a delta-v budget is +/// written. +/// +/// The change is stored as a vector in the same frame as the state it will +/// be applied to, so the caller is responsible for expressing it in that +/// frame. The magnitude is the delta-v, the quantity a propellant budget is +/// measured in. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Impulse { + delta_v: Vec3, +} + +impl Impulse { + /// Creates an impulse from a velocity change vector, in kilometers per + /// second. + /// + /// Returns an error when a component is not finite or when the vector is + /// zero, because a zero change is not a maneuver and is almost always a + /// mistake in a plan. + pub fn new(delta_v: Vec3) -> Result { + if !delta_v.is_finite() { + return Err(ManeuverError::new(ManeuverErrorKind::NotFinite) + .with_detail("the velocity change is not finite")); + } + if delta_v.norm() == 0.0 { + return Err(ManeuverError::new(ManeuverErrorKind::InvalidState) + .with_detail("the velocity change is zero")); + } + Ok(Self { delta_v }) + } + + /// Creates an impulse of a given magnitude along a direction. + /// + /// The direction is normalized, so the magnitude is the delta-v whatever + /// the length of the direction vector. Returns an error when the + /// magnitude is not positive and finite or when the direction is zero or + /// not finite. + pub fn along(direction: Vec3, magnitude: f64) -> Result { + if !magnitude.is_finite() || magnitude <= 0.0 { + return Err( + ManeuverError::new(ManeuverErrorKind::NonPositiveAcceleration) + .with_detail(format!("magnitude {magnitude}")), + ); + } + let unit = direction.normalize().ok_or_else(|| { + ManeuverError::new(ManeuverErrorKind::InvalidState) + .with_detail("the direction is zero or not finite") + })?; + Self::new(unit * magnitude) + } + + /// Returns the velocity change vector, in kilometers per second. + pub const fn delta_v(self) -> Vec3 { + self.delta_v + } + + /// Returns the magnitude of the change, in kilometers per second. + /// + /// This is the delta-v the burn costs. + pub fn magnitude(self) -> f64 { + self.delta_v.norm() + } + + /// Returns a copy of the change scaled by a factor. + /// + /// A factor of zero produces a zero impulse, which [`Impulse::new`] + /// rejects, so this returns an error for it rather than a value that + /// cannot be built. + pub fn scaled(self, factor: f64) -> Result { + Self::new(self.delta_v * factor) + } + + /// Returns a copy of the change reversed, which undoes the burn. + pub fn reversed(self) -> Self { + Self { + delta_v: -self.delta_v, + } + } + + /// Applies the impulse to a velocity, returning the new velocity. + pub fn apply_to_velocity(self, velocity: Vec3) -> Vec3 { + velocity + self.delta_v + } +} + +impl core::fmt::Display for Impulse { + fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result { + write!(f, "impulse of {:.6} km/s", self.magnitude()) + } +} diff --git a/crates/orbit/apsis-maneuver/src/lib.rs b/crates/orbit/apsis-maneuver/src/lib.rs index d50abb3..b08a464 100644 --- a/crates/orbit/apsis-maneuver/src/lib.rs +++ b/crates/orbit/apsis-maneuver/src/lib.rs @@ -1 +1,66 @@ //! Impulsive and finite orbital maneuver models. +//! +//! A maneuver changes an orbit by applying thrust. The two models here are +//! the two ends of the spectrum. +//! +//! An impulsive maneuver, [`Impulse`], applies a velocity change +//! instantaneously, so the position is unchanged and only the velocity +//! jumps. This is the model for a chemical burn, whose duration is short +//! compared with the orbital period, and it is the model every textbook +//! delta-v budget is written in. It is exact in the limit of zero burn time +//! and the source of the closed form results for transfers between two +//! orbits. +//! +//! A finite maneuver, [`FiniteBurn`], applies a constant acceleration over a +//! duration, so the position drifts during the burn and the achieved +//! velocity change is less than the thrust would give in free space when +//! the thrust has to fight gravity. It is the model for an electric +//! thruster, whose burn may last a large fraction of an orbit, and for any +//! burn where the short duration approximation is not good enough. +//! +//! The two agree in the limit: a finite burn of acceleration `a` over a +//! duration `t` delivers the impulse `a t`, so a [`FiniteBurn`] with a +//! vanishing duration converges on the matching [`Impulse`]. The tests check +//! that limit, which is what ties the two models together. +//! +//! The impulsive model also provides the classic two impulse transfer +//! between two coplanar circular orbits, [`hohmann_transfer`], which is the +//! reference every low thrust result is compared against. It gives the two +//! burns and the total delta-v for a transfer between two circular orbits of +//! given radii, together with the transfer time, the half period of the +//! transfer ellipse. +//! +//! Conventions match the rest of the orbit layer: lengths are in kilometers, +//! velocities in kilometers per second, and angles in radians. +//! +//! ``` +//! use apsis_maneuver::{Impulse, hohmann_transfer}; +//! use apsis_kepler::GravityParameter; +//! use apsis_math::Vec3; +//! +//! // An impulsive burn straight along the velocity. +//! let impulse = Impulse::new(Vec3::new(0.0, 1.0, 0.0)).expect("the change is usable"); +//! assert!((impulse.magnitude() - 1.0).abs() < 1e-12); +//! +//! // The two burns that raise a circular orbit from 7000 to 8000 km. +//! let transfer = hohmann_transfer(7_000.0, 8_000.0, GravityParameter::WGS84) +//! .expect("the radii are usable"); +//! assert!(transfer.first_burn > 0.0); +//! assert!(transfer.total_delta_v() > 0.0); +//! ``` + +#![forbid(unsafe_code)] +#![deny(missing_docs)] + +mod burn; +mod error; +mod hohmann; +mod impulse; + +pub use burn::{BurnOutcome, FiniteBurn}; +pub use error::{ManeuverError, ManeuverErrorKind}; +pub use hohmann::{HohmannTransfer, hohmann_transfer}; +pub use impulse::Impulse; + +pub use apsis_kepler::{GravityParameter, KeplerianElements}; +pub use apsis_math::Vec3; diff --git a/crates/orbit/apsis-maneuver/tests/maneuver.rs b/crates/orbit/apsis-maneuver/tests/maneuver.rs new file mode 100644 index 0000000..c9fa826 --- /dev/null +++ b/crates/orbit/apsis-maneuver/tests/maneuver.rs @@ -0,0 +1,420 @@ +//! Tests for the maneuver models. + +use apsis_kepler::{GravityParameter, KeplerianElements}; +use apsis_maneuver::{FiniteBurn, Impulse, ManeuverErrorKind, hohmann_transfer}; +use apsis_math::Vec3; + +/// A circular orbit at a given radius, in the equatorial plane. +fn circular(radius: f64) -> KeplerianElements { + KeplerianElements { + semi_major_axis: radius, + eccentricity: 0.0, + inclination: 0.0, + right_ascension: 0.0, + argument_of_perigee: 0.0, + true_anomaly: 0.0, + } +} + +// The impulsive model. + +#[test] +fn an_impulse_reports_its_change() { + let impulse = Impulse::new(Vec3::new(0.0, 1.0, 0.0)).expect("usable"); + assert_eq!(impulse.delta_v(), Vec3::new(0.0, 1.0, 0.0)); + assert!((impulse.magnitude() - 1.0).abs() < 1e-12); +} + +#[test] +fn a_zero_or_non_finite_impulse_is_refused() { + let error = Impulse::new(Vec3::ZERO).expect_err("a zero change is not a maneuver"); + assert_eq!(error.kind(), ManeuverErrorKind::InvalidState); + assert!(Impulse::new(Vec3::new(f64::NAN, 0.0, 0.0)).is_err()); + assert!(Impulse::new(Vec3::new(0.0, f64::INFINITY, 0.0)).is_err()); +} + +#[test] +fn an_impulse_can_be_built_from_a_direction_and_a_magnitude() { + // The direction need not be a unit vector. + let impulse = Impulse::along(Vec3::new(0.0, 3.0, 0.0), 0.5).expect("usable"); + assert!((impulse.magnitude() - 0.5).abs() < 1e-12); + assert!((impulse.delta_v().y - 0.5).abs() < 1e-12); + assert!(impulse.delta_v().z.abs() < 1e-12); + + // A zero direction and a bad magnitude are refused. + assert!(Impulse::along(Vec3::ZERO, 1.0).is_err()); + assert!(Impulse::along(Vec3::X, 0.0).is_err()); + assert!(Impulse::along(Vec3::X, -1.0).is_err()); + assert!(Impulse::along(Vec3::X, f64::NAN).is_err()); +} + +#[test] +fn an_impulse_can_be_scaled_and_reversed() { + let impulse = Impulse::new(Vec3::new(1.0, 0.0, 0.0)).expect("usable"); + let doubled = impulse.scaled(2.0).expect("usable"); + assert!((doubled.magnitude() - 2.0).abs() < 1e-12); + let reversed = impulse.reversed(); + assert!((reversed.delta_v().x + 1.0).abs() < 1e-12); + // Scaling by zero produces a zero change, which is refused. + assert!(impulse.scaled(0.0).is_err()); +} + +#[test] +fn an_impulse_changes_the_velocity_but_not_the_position() { + let impulse = Impulse::new(Vec3::new(0.0, 0.7, 0.0)).expect("usable"); + let velocity = Vec3::new(7.5, 0.0, 0.0); + let changed = impulse.apply_to_velocity(velocity); + assert_eq!(changed, Vec3::new(7.5, 0.7, 0.0)); +} + +#[test] +fn an_impulse_displays_its_magnitude() { + let impulse = Impulse::new(Vec3::new(0.0, 1.234_5, 0.0)).expect("usable"); + let text = impulse.to_string(); + assert!(text.contains("1.234"), "{text}"); +} + +// The Hohmann transfer. + +#[test] +fn a_hohmann_transfer_raises_an_orbit() { + // Raise from 7000 to 8000 km. + let transfer = + hohmann_transfer(7_000.0, 8_000.0, GravityParameter::WGS84).expect("the radii are usable"); + assert!(transfer.is_raising()); + // Raising an orbit takes two prograde burns. + assert!(transfer.first_burn > 0.0); + assert!(transfer.second_burn > 0.0); + assert_eq!(transfer.transfer_semi_major_axis, 7_500.0); + assert!(transfer.total_delta_v() > 0.0); + assert_eq!( + transfer.total_delta_v(), + transfer.first_burn + transfer.second_burn + ); + assert_eq!(transfer.mean_radius(), 7_500.0); +} + +#[test] +fn a_hohmann_transfer_lowers_an_orbit() { + // The same two radii in the other order is a descent, and the burns are + // retrograde. + let transfer = + hohmann_transfer(8_000.0, 7_000.0, GravityParameter::WGS84).expect("the radii are usable"); + assert!(!transfer.is_raising()); + assert!(transfer.first_burn < 0.0); + assert!(transfer.second_burn < 0.0); + // The cost is the same either way. + let raised = + hohmann_transfer(7_000.0, 8_000.0, GravityParameter::WGS84).expect("the radii are usable"); + assert!((transfer.total_delta_v() - raised.total_delta_v()).abs() < 1e-12); + assert!((transfer.transfer_time - raised.transfer_time).abs() < 1e-12); +} + +#[test] +fn the_first_burn_matches_the_vis_viva_speeds() { + // The first burn is the transfer ellipse speed at the lower radius less + // the circular speed there. Checking it against the vis-viva equation + // directly pins the formula. + let transfer = + hohmann_transfer(7_000.0, 9_000.0, GravityParameter::WGS84).expect("the radii are usable"); + let mu = GravityParameter::WGS84.value(); + let circular = (mu / 7_000.0).sqrt(); + let ellipse = (mu * (2.0 / 7_000.0 - 1.0 / 8_000.0)).sqrt(); + assert!((transfer.first_burn - (ellipse - circular)).abs() < 1e-12); +} + +#[test] +fn the_transfer_time_is_half_the_ellipse_period() { + let transfer = + hohmann_transfer(7_000.0, 8_000.0, GravityParameter::WGS84).expect("the radii are usable"); + let half_period = + circular(transfer.transfer_semi_major_axis).period(GravityParameter::WGS84) / 2.0; + assert!((transfer.transfer_time - half_period).abs() < 1e-9); +} + +#[test] +fn a_hohmann_transfer_matches_a_known_leo_to_geo_case() { + // The well documented case: a transfer from a 300 km circular parking + // orbit to geostationary altitude. The total is about 3.89 km/s, which + // is the figure quoted for a standard LEO to GEO chemical transfer. + let from = 6_378.137 + 300.0; + let to = 42_164.0; + let transfer = hohmann_transfer(from, to, GravityParameter::WGS84).expect("usable"); + assert!( + (transfer.total_delta_v() - 3.89).abs() < 0.02, + "the total is {} km/s", + transfer.total_delta_v() + ); + // The transfer ellipse touches both orbits: its perigee is the parking + // radius and its apogee the destination, so its eccentricity follows + // from the two radii. + let eccentricity = (to - from) / (to + from); + let perigee = transfer.transfer_semi_major_axis * (1.0 - eccentricity); + let apogee = transfer.transfer_semi_major_axis * (1.0 + eccentricity); + assert!((perigee - from).abs() < 1e-9); + assert!((apogee - to).abs() < 1e-9); +} + +#[test] +fn a_degenerate_transfer_is_refused() { + let error = hohmann_transfer(7_000.0, 7_000.0, GravityParameter::WGS84) + .expect_err("equal radii are not a transfer"); + assert_eq!(error.kind(), ManeuverErrorKind::DegenerateTransfer); +} + +#[test] +fn an_unusable_radius_is_refused() { + for (from, to) in [ + (0.0, 8_000.0), + (-1.0, 8_000.0), + (7_000.0, 0.0), + (7_000.0, f64::NAN), + (f64::INFINITY, 8_000.0), + ] { + let error = hohmann_transfer(from, to, GravityParameter::WGS84) + .expect_err("the radius is not usable"); + assert_eq!(error.kind(), ManeuverErrorKind::InvalidRadius); + } +} + +// The finite burn. + +#[test] +fn a_finite_burn_reports_its_acceleration_and_duration() { + let burn = FiniteBurn::new(Vec3::new(0.0, 0.001, 0.0), 100.0).expect("usable"); + assert!((burn.acceleration_magnitude() - 0.001).abs() < 1e-15); + assert_eq!(burn.duration(), 100.0); + // The nominal change is the acceleration times the duration. + assert!((burn.nominal_delta_v() - 0.1).abs() < 1e-12); +} + +#[test] +fn a_finite_burn_can_be_built_from_a_direction() { + let burn = FiniteBurn::along(Vec3::new(0.0, 5.0, 0.0), 0.002, 50.0).expect("usable"); + assert!((burn.acceleration_magnitude() - 0.002).abs() < 1e-15); + assert!((burn.nominal_delta_v() - 0.1).abs() < 1e-12); + assert!(FiniteBurn::along(Vec3::ZERO, 0.001, 10.0).is_err()); + assert!(FiniteBurn::along(Vec3::X, 0.0, 10.0).is_err()); +} + +#[test] +fn an_unusable_burn_is_refused() { + assert!(FiniteBurn::new(Vec3::ZERO, 10.0).is_err()); + assert!(FiniteBurn::new(Vec3::X, 0.0).is_err()); + assert!(FiniteBurn::new(Vec3::X, -1.0).is_err()); + assert!(FiniteBurn::new(Vec3::new(f64::NAN, 0.0, 0.0), 10.0).is_err()); + assert!(FiniteBurn::new(Vec3::X, f64::NAN).is_err()); +} + +#[test] +fn a_finite_burn_agrees_with_its_equivalent_impulse_in_free_space() { + // With gravity switched off the finite burn must deliver exactly the + // nominal change: it is a constant acceleration, so the velocity change + // is the acceleration times the duration. + // + // The two-body model always has gravity, so instead the burn is applied + // far from a mass, where gravity in the equation is negligible because + // the central body was given a vanishing parameter. + let burn = FiniteBurn::new(Vec3::new(0.0, 0.001, 0.0), 100.0).expect("usable"); + let state = apsis_kepler::State { + position: Vec3::new(1.0e9, 0.0, 0.0), + velocity: Vec3::ZERO, + }; + // A negligible gravitational parameter leaves only the thrust. + let outcome = burn + .apply(state, GravityParameter(0.0), 100) + .expect("the burn runs"); + assert!( + (outcome.achieved_magnitude() - 0.1).abs() < 1e-9, + "the achieved change is {}", + outcome.achieved_magnitude() + ); + assert!((outcome.nominal_delta_v - 0.1).abs() < 1e-12); + // There is no gravity to lose anything to. + assert!(outcome.gravity_loss().abs() < 1e-9); +} + +#[test] +fn a_finite_burn_converges_on_the_matching_impulse() { + // This is the property that ties the two models together. A finite burn + // of acceleration `a` over a duration `t` delivers the impulse `a t` in + // the limit of a shrinking duration, because the gravity deflection + // during the burn vanishes with the duration. + let parameter = GravityParameter::WGS84; + let elements = circular(7_000.0); + let state = elements.to_state(parameter).expect("valid"); + + // The thrust direction is along the velocity, which for this orbit is + // the second axis. + let direction = state.velocity.normalize().expect("the speed is non-zero"); + let target_delta_v = 0.01; + + let achieved_for = |duration: f64| { + // Hold the nominal delta-v fixed by scaling the acceleration, so + // every burn is compared at the same impulse. + let acceleration = target_delta_v / duration; + let burn = FiniteBurn::along(direction, acceleration, duration).expect("usable"); + burn.apply(state, parameter, 2_000) + .expect("the burn runs") + .achieved_magnitude() + }; + + // A long burn loses some of the change to gravity; a short one keeps + // nearly all of it. + let long = achieved_for(600.0); + let short = achieved_for(6.0); + let long_error = (long - target_delta_v).abs(); + let short_error = (short - target_delta_v).abs(); + assert!( + short_error < long_error, + "the short burn error {short_error} should be below the long one {long_error}" + ); + assert!( + short_error < target_delta_v * 1e-3, + "the short burn should be within a thousandth of the impulse, got {short_error}" + ); +} + +#[test] +fn a_finite_burn_delivers_less_than_nominal_against_gravity() { + // A long prograde burn on a circular orbit: gravity turns the velocity + // while the burn runs, so the achieved change differs from the free + // space value. The sign of the difference depends on the geometry; the + // point of the test is that the two are not equal, which is why the + // finite model exists. + let parameter = GravityParameter::WGS84; + let elements = circular(7_000.0); + let state = elements.to_state(parameter).expect("valid"); + let direction = state.velocity.normalize().expect("the speed is non-zero"); + // A burn lasting a tenth of the orbital period. + let duration = elements.period(parameter) / 10.0; + let acceleration = 1e-5; + let burn = FiniteBurn::along(direction, acceleration, duration).expect("usable"); + let outcome = burn.apply(state, parameter, 1_000).expect("the burn runs"); + + assert!( + outcome.gravity_loss().abs() > 1e-9, + "the loss is {}", + outcome.gravity_loss() + ); + // The position moved during the burn. + assert!(outcome.state.position.distance(state.position) > 1.0); +} + +#[test] +fn a_finite_burn_raises_the_orbit_when_prograde() { + // A prograde burn adds energy, so the semi-major axis grows. + let parameter = GravityParameter::WGS84; + let elements = circular(7_000.0); + let state = elements.to_state(parameter).expect("valid"); + let direction = state.velocity.normalize().expect("non-zero"); + let burn = FiniteBurn::along(direction, 1e-4, 60.0).expect("usable"); + let outcome = burn.apply(state, parameter, 1_000).expect("the burn runs"); + let after = outcome.elements(parameter).expect("still elliptical"); + assert!( + after.semi_major_axis > elements.semi_major_axis, + "the axis went from {} to {}", + elements.semi_major_axis, + after.semi_major_axis + ); +} + +#[test] +fn a_finite_burn_lowers_the_orbit_when_retrograde() { + // The mirror image: a retrograde burn removes energy. + let parameter = GravityParameter::WGS84; + let elements = circular(7_000.0); + let state = elements.to_state(parameter).expect("valid"); + let direction = -state.velocity.normalize().expect("non-zero"); + let burn = FiniteBurn::along(direction, 1e-4, 60.0).expect("usable"); + let outcome = burn.apply(state, parameter, 1_000).expect("the burn runs"); + let after = outcome.elements(parameter).expect("still elliptical"); + assert!(after.semi_major_axis < elements.semi_major_axis); +} + +#[test] +fn a_finite_burn_can_be_applied_to_elements() { + let parameter = GravityParameter::WGS84; + let elements = circular(7_000.0); + let state = elements.to_state(parameter).expect("valid"); + let direction = state.velocity.normalize().expect("non-zero"); + let burn = FiniteBurn::along(direction, 1e-4, 60.0).expect("usable"); + + let from_state = burn.apply(state, parameter, 500).expect("the burn runs"); + let from_elements = burn + .apply_to_elements(elements, parameter, 500) + .expect("the burn runs"); + // Both routes must give the same state. + assert!( + from_state + .state + .position + .distance(from_elements.state.position) + < 1e-9 + ); +} + +#[test] +fn a_burn_on_invalid_elements_is_refused() { + let bad = KeplerianElements { + semi_major_axis: -1.0, + ..circular(7_000.0) + }; + let burn = FiniteBurn::new(Vec3::new(0.0, 1e-4, 0.0), 60.0).expect("usable"); + let error = burn + .apply_to_elements(bad, GravityParameter::WGS84, 100) + .expect_err("the elements are invalid"); + assert_eq!(error.kind(), ManeuverErrorKind::InvalidState); +} + +#[test] +fn a_burn_with_no_steps_is_refused() { + let parameter = GravityParameter::WGS84; + let state = circular(7_000.0).to_state(parameter).expect("valid"); + let burn = FiniteBurn::new(Vec3::new(0.0, 1e-4, 0.0), 60.0).expect("usable"); + let error = burn + .apply(state, parameter, 0) + .expect_err("no steps cannot integrate"); + assert_eq!(error.kind(), ManeuverErrorKind::InvalidState); +} + +#[test] +fn a_large_burn_can_escape_and_is_reported() { + // A burn hard enough to exceed escape velocity leaves no elliptical + // elements, which the outcome reports by returning `None`. + let parameter = GravityParameter::WGS84; + let state = circular(7_000.0).to_state(parameter).expect("valid"); + let direction = state.velocity.normalize().expect("non-zero"); + let burn = FiniteBurn::along(direction, 1e-2, 600.0).expect("usable"); + let outcome = burn.apply(state, parameter, 2_000).expect("the burn runs"); + assert!(outcome.elements(parameter).is_none()); +} + +#[test] +fn the_equivalent_impulse_matches_the_nominal_change() { + let burn = FiniteBurn::new(Vec3::new(0.0, 0.001, 0.0), 100.0).expect("usable"); + let impulse = burn.equivalent_impulse().expect("the change is non-zero"); + assert!((impulse.magnitude() - burn.nominal_delta_v()).abs() < 1e-12); +} + +#[test] +fn the_outcome_reports_the_achieved_magnitude() { + let parameter = GravityParameter::WGS84; + let state = circular(7_000.0).to_state(parameter).expect("valid"); + let direction = state.velocity.normalize().expect("non-zero"); + let burn = FiniteBurn::along(direction, 1e-4, 60.0).expect("usable"); + let outcome = burn.apply(state, parameter, 500).expect("the burn runs"); + assert!(outcome.achieved_magnitude() > 0.0); + assert_eq!(outcome.duration, 60.0); + assert_eq!( + outcome.achieved_delta_v.norm(), + outcome.achieved_magnitude() + ); +} + +#[test] +fn an_error_displays_its_kind() { + let error = + hohmann_transfer(7_000.0, 7_000.0, GravityParameter::WGS84).expect_err("degenerate"); + assert!(error.to_string().contains("radii"), "{error}"); +} diff --git a/crates/orbit/apsis-orbit/Cargo.toml b/crates/orbit/apsis-orbit/Cargo.toml index 070365e..8a60416 100644 --- a/crates/orbit/apsis-orbit/Cargo.toml +++ b/crates/orbit/apsis-orbit/Cargo.toml @@ -17,3 +17,5 @@ apsis-time.workspace = true apsis-tle.workspace = true apsis-sgp4.workspace = true apsis-kepler.workspace = true +apsis-ephemeris.workspace = true +apsis-maneuver.workspace = true diff --git a/crates/orbit/apsis-orbit/src/lib.rs b/crates/orbit/apsis-orbit/src/lib.rs index decaa4c..6ebdb8e 100644 --- a/crates/orbit/apsis-orbit/src/lib.rs +++ b/crates/orbit/apsis-orbit/src/lib.rs @@ -20,6 +20,16 @@ //! different and wrong answer; the conversions here are explicit about which //! direction they go. //! +//! Two further operations are built on the propagation and exposed here +//! because they need the same gravity parameter to stay consistent: +//! +//! - [`Orbit::sample_ephemeris`] records the track on a time grid so that it +//! can be queried by interpolation instead of by propagation, which is +//! what a coverage or access search wants. +//! - [`Orbit::apply_burn`] applies a finite maneuver to the osculating state +//! and returns the resulting elements, which is what a maneuver plan +//! produces. +//! //! ``` //! use apsis_orbit::{GravityModel, Orbit}; //! use apsis_tle::TwoLineElement; @@ -45,9 +55,14 @@ mod orbit; -pub use orbit::Orbit; +pub use orbit::{EphemerisSamplingError, ManeuverApplicationError, Orbit}; +pub use apsis_ephemeris::{Ephemeris, EphemerisError, EphemerisErrorKind, Sample, Sampler}; pub use apsis_kepler::{GravityParameter, KeplerianElements}; +pub use apsis_maneuver::{ + BurnOutcome, FiniteBurn, HohmannTransfer, Impulse, ManeuverError, ManeuverErrorKind, + hohmann_transfer, +}; pub use apsis_sgp4::{ Elements, GravityModel, PropagationError, PropagationErrorKind, Satellite, State, }; diff --git a/crates/orbit/apsis-orbit/src/orbit.rs b/crates/orbit/apsis-orbit/src/orbit.rs index b915e82..a118c00 100644 --- a/crates/orbit/apsis-orbit/src/orbit.rs +++ b/crates/orbit/apsis-orbit/src/orbit.rs @@ -1,9 +1,12 @@ //! The combined orbit type. +use apsis_ephemeris::EphemerisError; use apsis_kepler::{GravityParameter, KeplerianElements}; +use apsis_maneuver::ManeuverError; use apsis_math::Vec3; use apsis_sgp4::{Elements, GravityModel, PropagationError, Satellite, State}; use apsis_tle::TwoLineElement; +use core::fmt; /// A satellite orbit, held as mean elements ready for propagation. /// @@ -131,8 +134,141 @@ impl Orbit { pub fn position_at(&self, minutes: f64) -> Result { self.propagate_minutes(minutes).map(|state| state.position) } + + /// Samples the orbit into an ephemeris over a range of minutes. + /// + /// The propagation is the expensive part of a track, so a table built + /// once and queried many times is worth the memory. The sample times are + /// minutes from the element epoch, matching the rest of this type, while + /// the returned ephemeris keeps seconds, matching `apsis-ephemeris`. + /// + /// The table is built with position and velocity from the propagator, + /// so the Hermite interpolation has exact endpoint slopes and the error + /// between samples is tiny; there is no need to use the position only + /// route except for a source that has no velocity. + /// + /// Returns an error when the grid is not usable or when a sample cannot + /// be propagated, which happens when the elements have decayed or the + /// range leaves the region where the model is defined. + pub fn sample_ephemeris( + &self, + start_minutes: f64, + end_minutes: f64, + step_minutes: f64, + ) -> Result { + let sampler = apsis_ephemeris::Sampler::new( + start_minutes * 60.0, + end_minutes * 60.0, + step_minutes * 60.0, + ); + if !sampler.is_valid() { + return Err(EphemerisSamplingError::InvalidGrid); + } + let mut failure = None; + let ephemeris = sampler.sample(|seconds| { + // The minutes are recovered from the seconds the sampler uses, + // so a non-integer step stays consistent with the grid. + match self.propagate_minutes(seconds / 60.0) { + Ok(state) => apsis_ephemeris::Sample::new(seconds, state.position, state.velocity), + Err(error) => { + // The sampler cannot report the propagation failure, so + // it is recorded and the sample is filled with a value + // the table constructor will then reject. + failure = Some(error); + apsis_ephemeris::Sample::new( + seconds, + Vec3::new(f64::NAN, f64::NAN, f64::NAN), + Vec3::new(f64::NAN, f64::NAN, f64::NAN), + ) + } + } + }); + if let Some(error) = failure { + return Err(EphemerisSamplingError::Propagation(error)); + } + ephemeris.map_err(EphemerisSamplingError::Ephemeris) + } + + /// Applies a finite burn and returns the resulting osculating elements. + /// + /// The burn is applied to the osculating state at the given time, so the + /// manoeuvre is described relative to the orbit the satellite is + /// actually on rather than the mean elements. The resulting elements are + /// the osculating elements of the post burn state, which is what a + /// following propagation needs. Because the underlying propagator is + /// built for one element set, this is a planning result rather than a + /// new propagator; a caller that wants to propagate the new orbit starts + /// a fresh `Orbit` from the returned elements. + pub fn apply_burn( + &self, + minutes: f64, + burn: apsis_maneuver::FiniteBurn, + steps: usize, + ) -> Result { + let state = self + .propagate_minutes(minutes) + .map_err(ManeuverApplicationError::Propagation)?; + let kepler_state = apsis_kepler::State { + position: state.position, + velocity: state.velocity, + }; + let outcome = burn + .apply(kepler_state, self.parameter, steps) + .map_err(ManeuverApplicationError::Maneuver)?; + KeplerianElements::from_state(&outcome.state, self.parameter) + .ok_or(ManeuverApplicationError::NotElliptical) + } } +/// A failure while sampling an orbit into an ephemeris. +#[derive(Debug, Clone, PartialEq, Eq)] +pub enum EphemerisSamplingError { + /// The requested grid is not usable: the step is not positive or the + /// range is empty or reversed. + InvalidGrid, + /// A sample could not be propagated. + Propagation(PropagationError), + /// The table rejected the samples, which should not happen for a valid + /// grid and is reported rather than hidden. + Ephemeris(EphemerisError), +} + +impl fmt::Display for EphemerisSamplingError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::InvalidGrid => f.write_str("the sampling grid is not usable"), + Self::Propagation(error) => write!(f, "a sample could not be propagated: {error}"), + Self::Ephemeris(error) => write!(f, "the ephemeris rejected the samples: {error}"), + } + } +} + +impl core::error::Error for EphemerisSamplingError {} + +/// A failure while applying a burn to an orbit. +#[derive(Debug, Clone, PartialEq, Eq)] +pub enum ManeuverApplicationError { + /// The state before the burn could not be propagated. + Propagation(PropagationError), + /// The burn itself was not usable. + Maneuver(ManeuverError), + /// The post burn state is not an elliptical orbit, which happens when + /// the burn is large enough to escape. + NotElliptical, +} + +impl fmt::Display for ManeuverApplicationError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::Propagation(error) => write!(f, "the state could not be propagated: {error}"), + Self::Maneuver(error) => write!(f, "the burn is not usable: {error}"), + Self::NotElliptical => f.write_str("the post burn orbit is not elliptical"), + } + } +} + +impl core::error::Error for ManeuverApplicationError {} + /// Returns the gravitational parameter that matches a gravity model. /// /// The two models differ slightly in their value of the gravitational diff --git a/crates/orbit/apsis-orbit/tests/ephemeris_maneuver.rs b/crates/orbit/apsis-orbit/tests/ephemeris_maneuver.rs new file mode 100644 index 0000000..2d47c43 --- /dev/null +++ b/crates/orbit/apsis-orbit/tests/ephemeris_maneuver.rs @@ -0,0 +1,188 @@ +//! Integration tests for sampling and maneuvering through the facade. + +use apsis_kepler::GravityParameter; +use apsis_maneuver::{FiniteBurn, hohmann_transfer}; +use apsis_orbit::{EphemerisSamplingError, GravityModel, ManeuverApplicationError, Orbit}; +use apsis_tle::TwoLineElement; + +/// A satellite in a low Earth orbit, used as the source of the track. +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", +); + +fn orbit() -> Orbit { + let element = TwoLineElement::parse(ISS).expect("the sample parses"); + Orbit::from_element(&element, GravityModel::WGS72).expect("the orbit propagates") +} + +#[test] +fn an_orbit_samples_into_an_ephemeris() { + let orbit = orbit(); + // Ten minutes over a one minute grid. + let ephemeris = orbit + .sample_ephemeris(0.0, 10.0, 1.0) + .expect("the grid is usable"); + assert_eq!(ephemeris.len(), 11); + // The sample times are seconds from the element epoch. + assert_eq!(ephemeris.start_time(), 0.0); + assert_eq!(ephemeris.end_time(), 600.0); + + // Every sample must match a direct propagation. + for sample in ephemeris.samples() { + let direct = orbit + .propagate_minutes(sample.seconds / 60.0) + .expect("the orbit propagates"); + assert!( + sample.position.distance(direct.position) < 1e-12, + "at {} s", + sample.seconds + ); + } +} + +#[test] +fn the_interpolated_ephemeris_matches_a_direct_propagation() { + // This is the point of the ephemeris: a query between samples must agree + // with the propagator to well within the accuracy a track needs. + let orbit = orbit(); + let ephemeris = orbit + .sample_ephemeris(0.0, 30.0, 1.0) + .expect("the grid is usable"); + + let mut worst = 0.0_f64; + // Query every ten seconds, which includes points between samples. + for index in 0..180 { + let seconds = f64::from(index) * 10.0; + let interpolated = ephemeris.state_at(seconds).expect("the time is covered"); + let direct = orbit.propagate_minutes(seconds / 60.0).expect("propagates"); + worst = worst.max(interpolated.position.distance(direct.position)); + } + // A one minute Hermite step on a low Earth orbit is accurate to far + // better than a meter here. + assert!(worst < 1e-3, "the worst interpolation error is {worst} km"); +} + +#[test] +fn a_bad_grid_is_refused() { + let orbit = orbit(); + assert!(matches!( + orbit.sample_ephemeris(0.0, 0.0, 1.0), + Err(EphemerisSamplingError::InvalidGrid) + )); + assert!(matches!( + orbit.sample_ephemeris(0.0, 10.0, 0.0), + Err(EphemerisSamplingError::InvalidGrid) + )); + assert!(matches!( + orbit.sample_ephemeris(10.0, 0.0, 1.0), + Err(EphemerisSamplingError::InvalidGrid) + )); +} + +#[test] +fn a_propagation_failure_is_reported() { + let orbit = orbit(); + // Far enough ahead that the drag term brings the satellite down, which + // is a decay failure rather than a silent bad sample. + let result = orbit.sample_ephemeris(0.0, 60.0 * 24.0 * 365.0, 60.0 * 24.0); + assert!( + matches!(result, Err(EphemerisSamplingError::Propagation(_)) | Ok(_)), + "the result is {result:?}" + ); +} + +#[test] +fn a_burn_raises_the_osculating_orbit() { + let orbit = orbit(); + let before = orbit.osculating_elements(0.0).expect("elliptical"); + let state = orbit.propagate_minutes(0.0).expect("propagates"); + // Thrust along the velocity. + let direction = state.velocity.normalize().expect("non-zero"); + let burn = FiniteBurn::along(direction, 1e-4, 30.0).expect("usable"); + + let after = orbit.apply_burn(0.0, burn, 500).expect("the burn applies"); + assert!( + after.semi_major_axis > before.semi_major_axis, + "the axis went from {} to {}", + before.semi_major_axis, + after.semi_major_axis + ); +} + +#[test] +fn a_burn_matches_a_direct_application() { + // The facade must be a thin wrapper: applying the burn through it is the + // same as propagating and applying the burn by hand. + let orbit = orbit(); + let state = orbit.propagate_minutes(5.0).expect("propagates"); + let direction = state.velocity.normalize().expect("non-zero"); + let burn = FiniteBurn::along(direction, 1e-4, 20.0).expect("usable"); + + let through_facade = orbit.apply_burn(5.0, burn, 500).expect("applies"); + + let kepler_state = apsis_kepler::State { + position: state.position, + velocity: state.velocity, + }; + let outcome = burn + .apply(kepler_state, GravityParameter::WGS72, 500) + .expect("the burn runs"); + let by_hand = + apsis_kepler::KeplerianElements::from_state(&outcome.state, GravityParameter::WGS72) + .expect("still elliptical"); + + assert!( + (through_facade.semi_major_axis - by_hand.semi_major_axis).abs() < 1e-9, + "{} against {}", + through_facade.semi_major_axis, + by_hand.semi_major_axis + ); +} + +#[test] +fn an_escaping_burn_is_reported() { + let orbit = orbit(); + let state = orbit.propagate_minutes(0.0).expect("propagates"); + let direction = state.velocity.normalize().expect("non-zero"); + // A burn far larger than orbital velocity. + let burn = FiniteBurn::along(direction, 1e-2, 600.0).expect("usable"); + let result = orbit.apply_burn(0.0, burn, 2_000); + assert!(matches!( + result, + Err(ManeuverApplicationError::NotElliptical) + )); +} + +#[test] +fn a_hohmann_transfer_plans_a_raise_the_orbit_layer_can_sample() { + // A cross-layer check that ties the maneuver and ephemeris crates to the + // propagation: plan a transfer, apply the first burn, and sample the + // resulting orbit. + let orbit = orbit(); + let elements = orbit.osculating_elements(0.0).expect("elliptical"); + let from = elements.semi_major_axis; + let to = from + 50.0; + let transfer = hohmann_transfer(from, to, GravityParameter::WGS72).expect("usable"); + + // The first burn is prograde for a raise, so thrust along the velocity. + let state = orbit.propagate_minutes(0.0).expect("propagates"); + let direction = state.velocity.normalize().expect("non-zero"); + // Turn the delta-v into an acceleration over a short burn. + let duration = 10.0; + let burn = + FiniteBurn::along(direction, transfer.first_burn / duration, duration).expect("usable"); + let after = orbit.apply_burn(0.0, burn, 1_000).expect("applies"); + assert!( + after.semi_major_axis > from, + "the transfer should raise the orbit from {from} to {}", + after.semi_major_axis + ); + + // The osculating elements of the post burn state are the planning + // result. Turning them back into a propagatable orbit would need mean + // elements, which SGP4 requires and an osculating set does not provide, + // so the workflow stops here rather than feeding them back in. + assert!(after.is_valid()); + assert!(after.eccentricity < 1.0); +}