ApsisRS
/
apsis-rs
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

orbit: add SGP4 and SDP4 propagator

This commit is contained in:
ApsisRS contributors 2026-09-26 14:16:23 +08:00
parent 7c66653fa1
commit 6eaaa86e30
11 changed files with 2665 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
crates/orbit/apsis-sgp4/Cargo.toml
View File

@ -11,3 +11,7 @@ publish = false
workspace = true
[dependencies]
apsis-core.workspace = true
apsis-math.workspace = true
apsis-tle.workspace = true
apsis-time.workspace = true

1113
crates/orbit/apsis-sgp4/src/deep_space.rs Normal file
View File

File diff suppressed because it is too large Load Diff

70
crates/orbit/apsis-sgp4/src/error.rs Normal file
View File

@ -0,0 +1,70 @@
//! Errors reported by the propagator.
use core::fmt;
/// The category of a propagation failure.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
#[non_exhaustive]
pub enum PropagationErrorKind {
/// The mean elements do not describe a satellite the model can
/// propagate. This includes an eccentricity at or above one, which is
/// not a closed orbit, and a mean motion that is not positive.
InvalidElements,
/// The satellite has decayed, meaning the drag term has brought it
/// below the surface of the Earth.
Decayed,
/// The propagation left the region where the model is defined.
OutOfRange,
}
impl fmt::Display for PropagationErrorKind {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
let text = match self {
Self::InvalidElements => "mean elements are not propagatable",
Self::Decayed => "satellite has decayed",
Self::OutOfRange => "propagation is out of range",
};
f.write_str(text)
}
}
/// A failure while initialising or propagating a satellite.
#[derive(Debug, Clone, PartialEq, Eq)]
pub struct PropagationError {
kind: PropagationErrorKind,
detail: Option<String>,
}
impl PropagationError {
/// Creates an error with the given category.
pub fn new(kind: PropagationErrorKind) -> Self {
Self { kind, detail: None }
}
/// Attaches a description of the value that caused the failure.
pub fn with_detail(mut self, detail: impl Into<String>) -> Self {
self.detail = Some(detail.into());
self
}
/// Returns the category of this failure.
pub fn kind(&self) -> PropagationErrorKind {
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 PropagationError {
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 PropagationError {}

180
crates/orbit/apsis-sgp4/src/init.rs Normal file
View File

@ -0,0 +1,180 @@
//! Initialisation of the propagator from a two-line element set.
//!
//! The element set gives mean elements in degrees, revolutions per day, and
//! a drag term. The propagator works in radians, radians per minute, and
//! Earth radii, so the initialisation converts the inputs and precomputes
//! every coefficient that depends only on the elements.
//!
//! The quantity `no_unkozai` is the mean motion after the correction that
//! removes the effect of the Kozai mean motion convention used in the
//! element set. The correction follows from the second zonal harmonic and
//! is solved by a fixed point iteration, because the corrected value appears
//! on both sides of the relation.
//!
//! Source: the correction and the coefficient definitions follow Vallado et
//! al., *Revisiting Spacetrack Report #3*, AIAA 2006-6753, 2006.
use crate::error::{PropagationError, PropagationErrorKind};
use crate::model::GravityModel;
use apsis_tle::TwoLineElement;
use core::f64::consts::PI;
/// The number of minutes in a day.
const MINUTES_PER_DAY: f64 = 1_440.0;
/// The number of radians in a revolution.
const TWO_PI: f64 = 2.0 * PI;
/// The period, in minutes, at or above which the deep space branch is used.
///
/// The theory states 225 minutes. Satellites below that period are handled
/// by the near-Earth branch, which is the case for all low Earth orbits.
const DEEP_SPACE_PERIOD_MINUTES: f64 = 225.0;
/// The elements of a satellite in the units the propagator uses.
///
/// Angles are in radians and mean motion is in radians per minute. Helpers
/// convert a [`TwoLineElement`] into this form.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Elements {
/// Epoch of the elements, expressed as a Modified Julian Date on the
/// scale the element set was written in.
pub epoch_mjd: f64,
/// `BSTAR` drag term, in inverse Earth radii.
pub bstar: f64,
/// Mean motion, in radians per minute.
pub mean_motion: f64,
/// Eccentricity, dimensionless.
pub eccentricity: f64,
/// Inclination, in radians.
pub inclination: f64,
/// Right ascension of the ascending node, in radians.
pub right_ascension: f64,
/// Argument of perigee, in radians.
pub argument_of_perigee: f64,
/// Mean anomaly, in radians.
pub mean_anomaly: f64,
}
impl Elements {
/// Converts a parsed element set into the units the propagator uses.
///
/// The epoch is taken as the UTC Modified Julian Date exactly as written
/// in the element set. The element set epoch is a UTC calendar instant,
/// but the propagator treats it as a uniform time origin, which is the
/// convention the theory was fitted under. Converting it to a uniform
/// scale first would shift the epoch by the leap second offset, about
/// 33 seconds, and that shift propagates into the solar and lunar series
/// of the deep space branch.
pub fn from_tle(element: &TwoLineElement) -> Option<Self> {
let epoch_mjd = element.epoch.mjd(apsis_time::TimeScale::Utc)?;
Some(Self {
epoch_mjd,
bstar: element.bstar,
mean_motion: element.mean_motion * TWO_PI / MINUTES_PER_DAY,
eccentricity: element.eccentricity,
inclination: element.inclination.to_radians(),
right_ascension: element.right_ascension.to_radians(),
argument_of_perigee: element.argument_of_perigee.to_radians(),
mean_anomaly: element.mean_anomaly.to_radians(),
})
}
}
/// Applies the Kozai to Brouwer mean motion correction.
///
/// The mean motion in an element set is a Kozai mean motion, which the
/// theory must convert to a Brouwer mean motion before use. The relation is
///
/// ```text
/// n0 = n / (1 + 1.5 k2 (3 cos^2 i - 1) / (a^2 sqrt(1 - e^2)) / a^2)
/// ```
///
/// written in the form used here as a fixed point iteration, because `a`
/// depends on `n0` through Kepler's third law. The iteration converges from
/// the unmodified mean motion in a few steps for every orbit that appears in
/// an element set.
/// Applies the Kozai to Brouwer mean motion correction.
///
/// The mean motion in an element set is a Kozai mean motion, which the
/// theory must convert to a Brouwer mean motion before use. Writing `no` for
/// the element set value and `xke` for the Kepler constant, the reference
/// implementation evaluates
///
/// ```text
/// ak = (xke / no)^(2/3)
/// d1 = 0.75 j2 (3 cos^2 i - 1) / (sqrt(1 - e^2) (1 - e^2))
/// del = d1 / ak^2
/// adel = ak (1 - del^2 - del (1/3 + 134 del^2 / 81))
/// del = d1 / adel^2
/// n0 = no / (1 + del)
/// ```
///
/// The second evaluation of `del` uses `adel` rather than `ak`, which is
/// what makes the correction consistent with the semi-major axis implied by
/// the corrected mean motion.
///
/// Source: Vallado et al., *Revisiting Spacetrack Report #3*, AIAA
/// 2006-6753, 2006.
pub fn unkozai_mean_motion(
mean_motion: f64,
eccentricity: f64,
inclination: f64,
model: GravityModel,
) -> f64 {
let xke = model.xke();
let j2 = model.j2;
let cosio = inclination.cos();
let cosio2 = cosio * cosio;
let omeosq = 1.0 - eccentricity * eccentricity;
let rteosq = omeosq.sqrt();
let ak = (xke / mean_motion).powf(2.0 / 3.0);
let d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq);
let mut delta = d1 / (ak * ak);
let adel = ak * (1.0 - delta * delta - delta * (1.0 / 3.0 + 134.0 * delta * delta / 81.0));
delta = d1 / (adel * adel);
mean_motion / (1.0 + delta)
}
/// Validates the mean elements before initialisation.
///
/// Returns an error when the orbit cannot be propagated, for example when
/// the eccentricity does not describe a closed orbit or the mean motion is
/// not positive.
pub fn validate(elements: &Elements) -> Result<(), PropagationError> {
if !elements.mean_motion.is_finite() || elements.mean_motion <= 0.0 {
return Err(PropagationError::new(PropagationErrorKind::InvalidElements)
.with_detail(format!("mean motion is {}", elements.mean_motion)));
}
if !elements.eccentricity.is_finite()
|| elements.eccentricity < 0.0
|| elements.eccentricity >= 1.0
{
return Err(PropagationError::new(PropagationErrorKind::InvalidElements)
.with_detail(format!("eccentricity is {}", elements.eccentricity)));
}
for (name, value) in [
("inclination", elements.inclination),
("right ascension", elements.right_ascension),
("argument of perigee", elements.argument_of_perigee),
("mean anomaly", elements.mean_anomaly),
] {
if !value.is_finite() {
return Err(PropagationError::new(PropagationErrorKind::InvalidElements)
.with_detail(format!("{name} is {value}")));
}
}
if !elements.bstar.is_finite() {
return Err(PropagationError::new(PropagationErrorKind::InvalidElements)
.with_detail(format!("bstar is {}", elements.bstar)));
}
Ok(())
}
/// Returns `true` when the elements describe a deep space orbit.
///
/// The test is on the orbital period, which follows from the mean motion.
pub fn is_deep_space(mean_motion: f64) -> bool {
TWO_PI / mean_motion >= DEEP_SPACE_PERIOD_MINUTES
}

71
crates/orbit/apsis-sgp4/src/lib.rs
View File

@ -1 +1,70 @@
//! SGP4/SDP4 analytical orbit propagator.
//! The SGP4 and SDP4 analytical orbit propagator.
//!
//! The propagator takes the mean orbital elements of a two-line element set
//! and returns the position and velocity of the satellite at a requested
//! instant. It is the model that the element sets are fitted to, so it is
//! the correct way to use them; a general purpose numerical integrator
//! given the same elements would drift away from the published positions.
//!
//! Two branches exist and the propagator selects between them by orbital
//! period. A satellite with a period below 225 minutes is propagated with
//! SGP4, which models the secular effect of the Earth's oblateness and a
//! drag term. A satellite at or above that period is propagated with SDP4,
//! which adds the resonance and the luni-solar perturbations that dominate
//! at that altitude.
//!
//! Inputs and outputs:
//!
//! - Elements are the mean elements from [`apsis_tle::TwoLineElement`].
//! Angles in the element set are in degrees and mean motion is in
//! revolutions per day; this crate converts to radians and radians per
//! minute internally.
//! - The output is expressed in the True Equator Mean Equinox (TEME) frame,
//! with position in kilometers and velocity in kilometers per second. TEME
//! is the frame the model defines; converting it to an Earth-fixed frame
//! requires the Earth rotation and the equation of the equinoxes and is
//! not part of this crate.
//! - Time is measured from the element set epoch, in minutes.
//!
//! Source: the algorithm follows Vallado et al., *Revisiting Spacetrack
//! Report #3*, AIAA 2006-6753, 2006, and the revised implementation
//! distributed with that paper.
//!
//! ```
//! use apsis_sgp4::{GravityModel, Satellite};
//! 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 satellite = Satellite::from_element(&element, GravityModel::WGS72)
//! .expect("sample elements are propagatable");
//!
//! // Propagate to the epoch of the elements.
//! let state = satellite.propagate_minutes(0.0).expect("epoch state exists");
//! // The position of a low Earth satellite is a few thousand kilometers
//! // from the center of the Earth.
//! assert!(state.position.norm() > 6_600.0);
//! assert!(state.position.norm() < 7_200.0);
//! ```
#![forbid(unsafe_code)]
#![deny(missing_docs)]
mod deep_space;
mod error;
mod init;
mod model;
mod propagate;
mod satellite;
mod state;
pub use error::{PropagationError, PropagationErrorKind};
pub use init::Elements;
pub use model::GravityModel;
pub use satellite::Satellite;
pub use state::State;
pub use apsis_tle::TwoLineElement;

108
crates/orbit/apsis-sgp4/src/model.rs Normal file
View File

@ -0,0 +1,108 @@
//! Gravity models used by the propagator.
//!
//! The propagator needs four numbers that describe the Earth: the
//! gravitational parameter, the equatorial radius, and the second and third
//! zonal harmonics. The SGP4 model was fitted to the WGS 72 values and the
//! standard element sets distributed by Space-Track are intended to be used
//! with them, so WGS 72 is the default.
//!
//! WGS 84 is offered because elements fitted against it also exist. Mixing
//! models with element sets fitted against the other model shifts the
//! position by kilometers, so the choice should match the source of the
//! elements.
//!
//! The derived values `xke`, `tumin`, and `j3oj2` follow from the four
//! defining numbers and are computed once at construction.
//!
//! Sources:
//!
//! - WGS 72 values as used by the SGP4 reference implementation in
//! Vallado et al., *Revisiting Spacetrack Report #3*, AIAA 2006-6753,
//! appendix. The values match the WGS 72 gravity field.
//! - WGS 84 values: NIMA TR8350.2, third edition, 2000.
/// A gravity model: the Earth constants used by the propagator.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct GravityModel {
/// Gravitational parameter, in cubic kilometers per second squared.
pub mu: f64,
/// Equatorial radius, in kilometers.
pub radius: f64,
/// Second zonal harmonic, dimensionless.
pub j2: f64,
/// Third zonal harmonic, dimensionless.
pub j3: f64,
/// Fourth zonal harmonic, dimensionless.
pub j4: f64,
}
impl GravityModel {
/// The WGS 72 gravity model.
///
/// This is the model the SGP4 theory was fitted to and the one that
/// matches the standard element sets.
pub const WGS72: Self = Self {
mu: 398_600.8,
radius: 6_378.135,
j2: 0.001_082_616,
j3: -0.000_002_538_81,
j4: -0.000_001_655_97,
};
/// The WGS 84 gravity model.
///
/// Use this only with element sets fitted against WGS 84.
pub const WGS84: Self = Self {
mu: 398_600.5,
radius: 6_378.137,
j2: 0.001_082_629_89,
j3: -0.000_002_532_153_1,
j4: -0.000_001_610_988_4,
};
/// Returns the reciprocal of the square root of the gravitational
/// parameter, in seconds per kilometer to the three halves.
pub fn reciprocal_sqrt_mu(self) -> f64 {
1.0 / self.mu.sqrt()
}
/// Returns the square root of the gravitational parameter, in kilometers
/// to the three halves per second.
pub fn sqrt_mu(self) -> f64 {
self.mu.sqrt()
}
/// Returns `xke`, the Kepler constant, in units of `er^1.5 / min`.
///
/// The propagator works in minutes and Earth radii, so the gravitational
/// parameter is converted by
///
/// ```text
/// xke = 60 / sqrt(radius^3 / mu)
/// ```
///
/// which is the angular rate, in radians per minute, of a circular orbit
/// of radius one Earth radius.
pub fn xke(self) -> f64 {
60.0 / (self.radius * self.radius * self.radius / self.mu).sqrt()
}
/// Returns the reciprocal of `xke`, in minutes.
///
/// This is the time unit of the model: the period of the reference
/// circular orbit divided by two pi.
pub fn tumin(self) -> f64 {
1.0 / self.xke()
}
/// Returns `j3 / j2`, which appears as a single ratio in the theory.
pub fn j3_over_j2(self) -> f64 {
self.j3 / self.j2
}
/// Returns `2 * xke^2`, a term that appears in the mean motion
/// correction.
pub fn x2o3(self) -> f64 {
2.0 / 3.0
}
}

299
crates/orbit/apsis-sgp4/src/propagate.rs Normal file
View File

@ -0,0 +1,299 @@
//! The near-Earth propagation step.
//!
//! Given a satellite record and a time from epoch, this module evaluates the
//! secular and periodic terms of the SGP4 theory and returns the state in
//! the True Equator Mean Equinox frame.
//!
//! The mean elements advance as polynomials in time: the mean anomaly,
//! argument of perigee, and right ascension are linear, the right ascension
//! has an additional quadratic term, and the semi-major axis, eccentricity,
//! and mean anomaly pick up further drag terms for orbits above the
//! simplified threshold. The orbit is then solved for the eccentric anomaly
//! by Newton's method, and the short period terms are applied.
//!
//! Source: the sequence follows Vallado et al., *Revisiting Spacetrack
//! Report #3*, AIAA 2006-6753, 2006, and the reference implementation
//! distributed with that paper.
use crate::error::{PropagationError, PropagationErrorKind};
use crate::satellite::{Method, Satellite, TWO_PI};
use crate::state::State;
use apsis_math::Vec3;
/// The number of iterations allowed when solving for the eccentric anomaly.
const KEPLER_MAX_ITERATIONS: usize = 10;
/// The convergence threshold on the eccentric anomaly correction.
const KEPLER_TOLERANCE: f64 = 1.0e-12;
/// The largest correction applied in one Kepler step.
///
/// Clamping keeps the iteration stable for the eccentric orbits that the
/// near circular approximation of the theory handles poorly.
const KEPLER_MAX_STEP: f64 = 0.95;
/// The smallest eccentricity used after the drag terms are applied.
///
/// A value that would round to zero makes several divisions in the short
/// period terms singular, so it is clamped to this floor.
const MINIMUM_ECCENTRICITY: f64 = 1.0e-6;
/// The divisor used when the inclination is close to pi.
const INCLINATION_GUARD: f64 = 1.5e-12;
/// Propagates a satellite by a number of minutes from its element epoch.
pub(crate) fn propagate(satellite: &Satellite, minutes: f64) -> Result<State, PropagationError> {
let mean_anomaly = satellite.mo + satellite.mdot * minutes;
let secular_argument = satellite.argpo + satellite.argpdot * minutes;
let secular_node = satellite.nodeo + satellite.nodedot * minutes;
let time_squared = minutes * minutes;
let mut node = secular_node + satellite.nodecf * time_squared;
let mut tempa = 1.0 - satellite.cc1 * minutes;
let mut tempe = satellite.bstar * satellite.cc4 * minutes;
let mut templ = satellite.t2cof * time_squared;
let mut argument = secular_argument;
let mut anomaly = mean_anomaly;
if !satellite.isimp {
let delta_argument = satellite.omgcof * minutes;
let drag_anomaly = 1.0 + satellite.eta * mean_anomaly.cos();
let delta_anomaly =
satellite.xmcof * (drag_anomaly * drag_anomaly * drag_anomaly - satellite.delmo);
let total = delta_argument + delta_anomaly;
anomaly = mean_anomaly + total;
argument = secular_argument - total;
let time_cubed = time_squared * minutes;
let time_fourth = time_cubed * minutes;
tempa -=
satellite.d2 * time_squared + satellite.d3 * time_cubed + satellite.d4 * time_fourth;
tempe += satellite.bstar * satellite.cc5 * (anomaly.sin() - satellite.sinmao);
templ = templ
+ satellite.t3cof * time_cubed
+ time_fourth * (satellite.t4cof + minutes * satellite.t5cof);
}
let mut eccentricity = satellite.ecco;
let mut inclination = satellite.inclo;
let mut mean_motion = satellite.no_unkozai;
// The deep space branch replaces the mean elements here.
if satellite.method == Method::DeepSpace {
let deep = crate::deep_space::dspace(
satellite,
minutes,
eccentricity,
argument,
inclination,
anomaly,
node,
mean_motion,
);
eccentricity = deep.eccentricity;
argument = deep.argument;
inclination = deep.inclination;
anomaly = deep.mean_anomaly;
node = deep.node;
mean_motion = deep.mean_motion;
}
if mean_motion <= 0.0 {
return Err(super::satellite::non_positive_mean_motion(mean_motion));
}
// Recover the semi-major axis including the drag scaling, then recompute
// the mean motion from it so the two stay consistent.
let semi_major = (satellite.xke() / mean_motion).powf(2.0 / 3.0) * tempa * tempa;
mean_motion = satellite.xke() / semi_major.powf(1.5);
eccentricity -= tempe;
if eccentricity >= 1.0 || eccentricity < -0.001 {
return Err(PropagationError::new(PropagationErrorKind::InvalidElements)
.with_detail(format!("propagated eccentricity is {eccentricity}")));
}
let eccentricity = eccentricity.max(MINIMUM_ECCENTRICITY);
anomaly += satellite.no_unkozai * templ;
let mut longitude = anomaly + argument + node;
// Reduce the angles to a single turn. The reference implementation uses
// a remainder that keeps the sign of the dividend for the right
// ascension and a positive remainder elsewhere; both are reproduced here
// because the eccentric anomaly iteration is sensitive to the branch.
node = if node >= 0.0 {
node % TWO_PI
} else {
-((-node) % TWO_PI)
};
argument %= TWO_PI;
longitude %= TWO_PI;
anomaly = (longitude - argument - node) % TWO_PI;
let mut perturbed_eccentricity = eccentricity;
let mut perturbed_inclination = inclination;
let mut perturbed_node = node;
let mut perturbed_argument = argument;
let mut perturbed_anomaly = anomaly;
let mut sin_perturbed = inclination.sin();
let mut cos_perturbed = inclination.cos();
if satellite.method == Method::DeepSpace {
let result = crate::deep_space::dpper(
satellite,
minutes,
false,
perturbed_eccentricity,
perturbed_inclination,
perturbed_node,
perturbed_argument,
perturbed_anomaly,
);
perturbed_eccentricity = result.eccentricity;
perturbed_inclination = result.inclination;
perturbed_node = result.node;
perturbed_argument = result.argument;
perturbed_anomaly = result.mean_anomaly;
if perturbed_inclination < 0.0 {
perturbed_inclination = -perturbed_inclination;
perturbed_node += core::f64::consts::PI;
perturbed_argument -= core::f64::consts::PI;
}
if perturbed_eccentricity < 0.0 || perturbed_eccentricity > 1.0 {
return Err(
PropagationError::new(PropagationErrorKind::InvalidElements).with_detail(format!(
"perturbed eccentricity is {perturbed_eccentricity}"
)),
);
}
sin_perturbed = perturbed_inclination.sin();
cos_perturbed = perturbed_inclination.cos();
}
// The `aycof` and `xlcof` used below come from the perturbed
// inclination in the deep space branch and from the epoch inclination
// otherwise.
let (aycof, xlcof) = if satellite.method == Method::DeepSpace {
let aycof = -0.5 * satellite.model.j3_over_j2() * sin_perturbed;
let xlcof = if (cos_perturbed + 1.0).abs() > INCLINATION_GUARD {
-0.25 * satellite.model.j3_over_j2() * sin_perturbed * (3.0 + 5.0 * cos_perturbed)
/ (1.0 + cos_perturbed)
} else {
-0.25 * satellite.model.j3_over_j2() * sin_perturbed * (3.0 + 5.0 * cos_perturbed)
/ INCLINATION_GUARD
};
(aycof, xlcof)
} else {
(satellite.aycof, satellite.xlcof)
};
let axnl = perturbed_eccentricity * perturbed_argument.cos();
let temp = 1.0 / (semi_major * (1.0 - perturbed_eccentricity * perturbed_eccentricity));
let aynl = perturbed_eccentricity * perturbed_argument.sin() + temp * aycof;
let xl = perturbed_anomaly + perturbed_argument + perturbed_node + temp * xlcof * axnl;
let u = (xl - perturbed_node) % TWO_PI;
// Solve Kepler's equation by Newton's method.
let mut eo1 = u;
let (mut sineo1, mut coseo1);
let mut iterations = 0;
loop {
sineo1 = eo1.sin();
coseo1 = eo1.cos();
let denominator = 1.0 - coseo1 * axnl - sineo1 * aynl;
let mut correction = (u - aynl * coseo1 + axnl * sineo1 - eo1) / denominator;
if correction.abs() >= KEPLER_MAX_STEP {
correction = if correction > 0.0 {
KEPLER_MAX_STEP
} else {
-KEPLER_MAX_STEP
};
}
eo1 += correction;
iterations += 1;
if correction.abs() < KEPLER_TOLERANCE || iterations >= KEPLER_MAX_ITERATIONS {
break;
}
}
sineo1 = eo1.sin();
coseo1 = eo1.cos();
let ecose = axnl * coseo1 + aynl * sineo1;
let esine = axnl * sineo1 - aynl * coseo1;
let el2 = axnl * axnl + aynl * aynl;
let pl = semi_major * (1.0 - el2);
if pl < 0.0 {
return Err(PropagationError::new(PropagationErrorKind::InvalidElements)
.with_detail(format!("semi-latus rectum is {pl}")));
}
let rl = semi_major * (1.0 - ecose);
let rdotl = semi_major.sqrt() * esine / rl;
let rvdotl = pl.sqrt() / rl;
let betal = (1.0 - el2).sqrt();
let temp = esine / (1.0 + betal);
let sinu = semi_major / rl * (sineo1 - aynl - axnl * temp);
let cosu = semi_major / rl * (coseo1 - axnl + aynl * temp);
let mut su = sinu.atan2(cosu);
let sin2u = (cosu + cosu) * sinu;
let cos2u = 1.0 - 2.0 * sinu * sinu;
let temp = 1.0 / pl;
let temp1 = 0.5 * satellite.model.j2 * temp;
let temp2 = temp1 * temp;
let (con41, x1mth2, x7thm1) = if satellite.method == Method::DeepSpace {
let cosisq = cos_perturbed * cos_perturbed;
(3.0 * cosisq - 1.0, 1.0 - cosisq, 7.0 * cosisq - 1.0)
} else {
(satellite.con41, satellite.x1mth2, satellite.x7thm1)
};
let mrt = rl * (1.0 - 1.5 * temp2 * betal * con41) + 0.5 * temp1 * x1mth2 * cos2u;
su -= 0.25 * temp2 * x7thm1 * sin2u;
let xnode = perturbed_node + 1.5 * temp2 * cos_perturbed * sin2u;
let xinc = perturbed_inclination + 1.5 * temp2 * cos_perturbed * sin_perturbed * cos2u;
let mvt = rdotl - mean_motion * temp1 * x1mth2 * sin2u / satellite.xke();
let rvdot = rvdotl + mean_motion * temp1 * (x1mth2 * cos2u + 1.5 * con41) / satellite.xke();
let sinsu = su.sin();
let cossu = su.cos();
let snod = xnode.sin();
let cnod = xnode.cos();
let sini = xinc.sin();
let cosi = xinc.cos();
let xmx = -snod * cosi;
let xmy = cnod * cosi;
let ux = xmx * sinsu + cnod * cossu;
let uy = xmy * sinsu + snod * cossu;
let uz = sini * sinsu;
let vx = xmx * cossu - cnod * sinsu;
let vy = xmy * cossu - snod * sinsu;
let vz = sini * cossu;
let radius = mrt * satellite.model.radius;
let position = Vec3::new(radius * ux, radius * uy, radius * uz);
let velocity_scale = satellite.model.radius * satellite.xke() / 60.0;
let velocity = Vec3::new(
(mvt * ux + rvdot * vx) * velocity_scale,
(mvt * uy + rvdot * vy) * velocity_scale,
(mvt * uz + rvdot * vz) * velocity_scale,
);
if mrt < 1.0 {
return Err(PropagationError::new(PropagationErrorKind::Decayed)
.with_detail(format!("radius is {mrt} Earth radii")));
}
Ok(State {
position,
velocity,
semi_major_axis: semi_major * satellite.model.radius,
})
}

416
crates/orbit/apsis-sgp4/src/satellite.rs Normal file
View File

@ -0,0 +1,416 @@
//! The satellite record and its initialisation from mean elements.
use crate::error::{PropagationError, PropagationErrorKind};
use crate::init::{self, Elements};
use crate::model::GravityModel;
use crate::state::State;
use apsis_tle::TwoLineElement;
use core::f64::consts::PI;
/// The number of radians in a revolution.
pub(crate) const TWO_PI: f64 = 2.0 * PI;
/// The number of minutes in a day.
pub(crate) const MINUTES_PER_DAY: f64 = 1_440.0;
/// Whether the near-Earth or the deep space branch applies.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(crate) enum Method {
/// Period below 225 minutes.
NearEarth,
/// Period at or above 225 minutes.
DeepSpace,
}
/// A satellite prepared for propagation.
///
/// Construct one with [`Satellite::from_element`], then call
/// [`Satellite::propagate_minutes`] or [`Satellite::propagate_at`] for each
/// instant of interest. The initialisation work, including the mean motion
/// correction and every coefficient that depends only on the elements, is
/// done once at construction.
#[derive(Debug, Clone)]
pub struct Satellite {
/// The model this satellite was initialised against.
pub(crate) model: GravityModel,
/// Corrected mean motion, in radians per minute.
pub(crate) no_unkozai: f64,
/// Eccentricity.
pub(crate) ecco: f64,
/// Inclination, in radians.
pub(crate) inclo: f64,
/// Right ascension of the ascending node, in radians.
pub(crate) nodeo: f64,
/// Argument of perigee, in radians.
pub(crate) argpo: f64,
/// Mean anomaly, in radians.
pub(crate) mo: f64,
/// `BSTAR` drag term, in inverse Earth radii.
pub(crate) bstar: f64,
/// Epoch as a Modified Julian Date on the TAI scale.
pub(crate) epoch_mjd: f64,
/// Which branch of the theory applies.
pub(crate) method: Method,
/// Set when the short period terms are omitted.
pub(crate) isimp: bool,
/// `3 cos^2 i - 1`.
pub(crate) con41: f64,
/// `1 - cos^2 i`.
pub(crate) x1mth2: f64,
/// `7 cos^2 i - 1`.
pub(crate) x7thm1: f64,
/// `-0.5 j3 / j2 sin i`.
pub(crate) aycof: f64,
/// The `j3` contribution to the argument of latitude.
pub(crate) xlcof: f64,
/// First drag coefficient.
pub(crate) cc1: f64,
/// Fourth drag coefficient.
pub(crate) cc4: f64,
/// Fifth drag coefficient.
pub(crate) cc5: f64,
/// Second order drag coefficient.
pub(crate) d2: f64,
/// Third order drag coefficient.
pub(crate) d3: f64,
/// Fourth order drag coefficient.
pub(crate) d4: f64,
/// `(1 + eta cos M)^3` at epoch.
pub(crate) delmo: f64,
/// `eta`.
pub(crate) eta: f64,
/// Contribution of drag to the argument of perigee.
pub(crate) omgcof: f64,
/// `sin M` at epoch.
pub(crate) sinmao: f64,
/// Quadratic coefficient of the time polynomial.
pub(crate) t2cof: f64,
/// Cubic coefficient of the time polynomial.
pub(crate) t3cof: f64,
/// Quartic coefficient of the time polynomial.
pub(crate) t4cof: f64,
/// Quintic coefficient of the time polynomial.
pub(crate) t5cof: f64,
/// Contribution of drag to the mean anomaly.
pub(crate) xmcof: f64,
/// Rate of change of the mean anomaly, in radians per minute.
pub(crate) mdot: f64,
/// Rate of change of the argument of perigee, in radians per minute.
pub(crate) argpdot: f64,
/// Rate of change of the right ascension, in radians per minute.
pub(crate) nodedot: f64,
/// Quadratic term of the right ascension.
pub(crate) nodecf: f64,
/// Greenwich sidereal time at the epoch, in radians.
pub(crate) gsto: f64,
/// Deep space resonance coefficients, all zero for near-Earth orbits.
pub(crate) ds: crate::deep_space::DeepSpace,
}
impl Satellite {
/// Initialises a satellite from a parsed element set.
///
/// Returns `None` when the mean elements do not describe an orbit the
/// model can propagate, for example when the eccentricity is at least
/// one. A caller that needs the reason should use
/// [`Satellite::try_from_element`].
pub fn from_element(element: &TwoLineElement, model: GravityModel) -> Option<Self> {
Self::try_from_element(element, model).ok()
}
/// Initialises a satellite from a parsed element set.
///
/// Returns an error describing which element is not propagatable when
/// the elements are rejected.
pub fn try_from_element(
element: &TwoLineElement,
model: GravityModel,
) -> Result<Self, PropagationError> {
let elements = Elements::from_tle(element).ok_or_else(|| {
PropagationError::new(PropagationErrorKind::InvalidElements)
.with_detail("element set epoch is outside the leap second table")
})?;
Self::try_from_elements(&elements, model)
}
/// Initialises a satellite from mean elements already in model units.
pub fn try_from_elements(
elements: &Elements,
model: GravityModel,
) -> Result<Self, PropagationError> {
init::validate(elements)?;
Ok(Self::initialise(elements, model))
}
/// Returns the epoch of the elements as a Modified Julian Date on the
/// TAI scale.
pub fn epoch_mjd(&self) -> f64 {
self.epoch_mjd
}
/// Returns `true` when the deep space branch applies.
pub fn is_deep_space(&self) -> bool {
self.method == Method::DeepSpace
}
/// Returns the number of minutes from the element set epoch to the
/// given instant.
///
/// The instant is given as a Modified Julian Date on the TAI scale, the
/// same scale as [`Satellite::epoch_mjd`].
pub fn minutes_since_epoch(&self, tai_mjd: f64) -> f64 {
(tai_mjd - self.epoch_mjd) * MINUTES_PER_DAY
}
/// Returns the state of the satellite at a given instant.
///
/// The instant is a Modified Julian Date on the TAI scale.
pub fn propagate_at(&self, tai_mjd: f64) -> Result<State, PropagationError> {
self.propagate_minutes(self.minutes_since_epoch(tai_mjd))
}
/// Returns the state of the satellite at a number of minutes from the
/// element set epoch.
///
/// Negative values propagate backwards. The result is in the True
/// Equator Mean Equinox frame, with position in kilometers and velocity
/// in kilometers per second.
///
/// Returns [`PropagationErrorKind::Decayed`](crate::PropagationErrorKind::Decayed)
/// when the propagated radius falls below one Earth radius, which the
/// model treats as the satellite having re-entered.
pub fn propagate_minutes(&self, minutes: f64) -> Result<State, PropagationError> {
crate::propagate::propagate(self, minutes)
}
/// Returns the Kepler constant of the model, in radians per minute per
/// Earth radius to the three halves.
pub fn xke(&self) -> f64 {
self.model.xke()
}
}
impl Satellite {
/// Builds the record, including every coefficient that depends only on
/// the epoch elements.
fn initialise(elements: &Elements, model: GravityModel) -> Self {
let xke = model.xke();
let j2 = model.j2;
let j4 = model.j4;
let ecco = elements.eccentricity;
let inclo = elements.inclination;
let no_kozai = elements.mean_motion;
let no_unkozai = init::unkozai_mean_motion(no_kozai, ecco, inclo, model);
let method = if init::is_deep_space(no_unkozai) {
Method::DeepSpace
} else {
Method::NearEarth
};
let cosio = inclo.cos();
let sinio = inclo.sin();
let cosio2 = cosio * cosio;
let cosio4 = cosio2 * cosio2;
let omeosq = 1.0 - ecco * ecco;
let rteosq = omeosq.sqrt();
let ao = (xke / no_unkozai).powf(2.0 / 3.0);
let po = ao * omeosq;
let posq = po * po;
let rp = ao * (1.0 - ecco);
let con42 = 1.0 - 5.0 * cosio2;
let con41 = -con42 - cosio2 - cosio2;
// The aerodynamic parameters. `sfour` is the effective radius at
// which the atmosphere is treated as having constant density, and
// `qzms24` is the corresponding density scale raised to the fourth
// power.
let ss = 78.0 / model.radius + 1.0;
let qzms2t = ((120.0 - 78.0) / model.radius).powi(4);
let mut isimp = false;
if rp < 220.0 / model.radius + 1.0 {
isimp = true;
}
let mut sfour = ss;
let mut qzms24 = qzms2t;
let perige = (rp - 1.0) * model.radius;
if perige < 156.0 {
sfour = perige - 78.0;
if perige < 98.0 {
sfour = 20.0;
}
qzms24 = ((120.0 - sfour) / model.radius).powi(4);
sfour = sfour / model.radius + 1.0;
}
let pinvsq = 1.0 / posq;
let tsi = 1.0 / (ao - sfour);
let eta = ao * ecco * tsi;
let etasq = eta * eta;
let eeta = ecco * eta;
let psisq = (1.0 - etasq).abs();
let coef = qzms24 * tsi.powi(4);
let coef1 = coef / psisq.powf(3.5);
let cc2 = coef1
* no_unkozai
* (ao * (1.0 + 1.5 * etasq + eeta * (4.0 + etasq))
+ 0.375 * j2 * tsi / psisq * con41 * (8.0 + 3.0 * etasq * (8.0 + etasq)));
let cc1 = elements.bstar * cc2;
let mut cc3 = 0.0;
if ecco > 1.0e-4 {
cc3 = -2.0 * coef * tsi * model.j3_over_j2() * no_unkozai * sinio / ecco;
}
let x1mth2 = 1.0 - cosio2;
let cc4 = 2.0
* no_unkozai
* coef1
* ao
* omeosq
* (eta * (2.0 + 0.5 * etasq) + ecco * (0.5 + 2.0 * etasq)
- j2 * tsi / (ao * psisq)
* (-3.0 * con41 * (1.0 - 2.0 * eeta + etasq * (1.5 - 0.5 * eeta))
+ 0.75
* x1mth2
* (2.0 * etasq - eeta * (1.0 + etasq))
* (2.0 * elements.argument_of_perigee).cos()));
let cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * (etasq + eeta) + eeta * etasq);
let temp1 = 1.5 * j2 * pinvsq * no_unkozai;
let temp2 = 0.5 * temp1 * j2 * pinvsq;
let temp3 = -0.46875 * j4 * pinvsq * pinvsq * no_unkozai;
let mdot = no_unkozai
+ 0.5 * temp1 * rteosq * con41
+ 0.0625 * temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4);
let argpdot = -0.5 * temp1 * con42
+ 0.0625 * temp2 * (7.0 - 114.0 * cosio2 + 395.0 * cosio4)
+ temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4);
let xhdot1 = -temp1 * cosio;
let nodedot = xhdot1
+ (0.5 * temp2 * (4.0 - 19.0 * cosio2) + 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio;
let omgcof = elements.bstar * cc3 * elements.argument_of_perigee.cos();
let mut xmcof = 0.0;
if ecco > 1.0e-4 {
xmcof = -2.0 / 3.0 * coef * elements.bstar / eeta;
}
let nodecf = 3.5 * omeosq * xhdot1 * cc1;
let t2cof = 1.5 * cc1;
// The divisor guards the singularity at inclination pi, where the
// expression is finite in the limit. `1.5e-12` is the threshold used
// by the reference implementation.
const INCLINATION_GUARD: f64 = 1.5e-12;
let xlcof = if (cosio + 1.0).abs() > INCLINATION_GUARD {
-0.25 * model.j3_over_j2() * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio)
} else {
-0.25 * model.j3_over_j2() * sinio * (3.0 + 5.0 * cosio) / INCLINATION_GUARD
};
let aycof = -0.5 * model.j3_over_j2() * sinio;
let delmotemp = 1.0 + eta * elements.mean_anomaly.cos();
let delmo = delmotemp * delmotemp * delmotemp;
let sinmao = elements.mean_anomaly.sin();
let x7thm1 = 7.0 * cosio2 - 1.0;
let mut satellite = Self {
model,
no_unkozai,
ecco,
inclo,
nodeo: elements.right_ascension,
argpo: elements.argument_of_perigee,
mo: elements.mean_anomaly,
bstar: elements.bstar,
epoch_mjd: elements.epoch_mjd,
method,
isimp,
con41,
x1mth2,
x7thm1,
aycof,
xlcof,
cc1,
cc4,
cc5,
d2: 0.0,
d3: 0.0,
d4: 0.0,
delmo,
eta,
omgcof,
sinmao,
t2cof,
t3cof: 0.0,
t4cof: 0.0,
t5cof: 0.0,
xmcof,
mdot,
argpdot,
nodedot,
nodecf,
gsto: gstime(elements.epoch_mjd + 2_400_000.5),
ds: crate::deep_space::DeepSpace::ZERO,
};
// The higher order drag terms are omitted for the simplified branch.
if !isimp {
let cc1sq = cc1 * cc1;
let d2 = 4.0 * ao * tsi * cc1sq;
let temp = d2 * tsi * cc1 / 3.0;
let d3 = (17.0 * ao + sfour) * temp;
let d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * cc1;
satellite.d2 = d2;
satellite.d3 = d3;
satellite.d4 = d4;
satellite.t3cof = d2 + 2.0 * cc1sq;
satellite.t4cof = 0.25 * (3.0 * d3 + cc1 * (12.0 * d2 + 10.0 * cc1sq));
satellite.t5cof = 0.2
* (3.0 * d4 + 12.0 * cc1 * d3 + 6.0 * d2 * d2 + 15.0 * cc1sq * (2.0 * d2 + cc1sq));
}
if method == Method::DeepSpace {
satellite.isimp = true;
satellite.ds = crate::deep_space::initialise(&satellite);
}
satellite
}
}
/// Returns Greenwich sidereal time for a UT1 Julian Date, in radians.
///
/// The expression is the one used by the reference implementation. It
/// expands the ratio of universal time to sidereal time and is accurate for
/// the purpose of initialising the deep space resonance, which is the only
/// place the value is used.
///
/// Source: Vallado et al., *Revisiting Spacetrack Report #3*, AIAA
/// 2006-6753, 2006.
pub(crate) fn gstime(julian_date_ut1: f64) -> f64 {
let centuries = (julian_date_ut1 - 2_451_545.0) / 36_525.0;
let mut seconds = -6.2e-6 * centuries * centuries * centuries
+ 0.093_104 * centuries * centuries
+ (876_600.0 * 3_600.0 + 8_640_184.812_866) * centuries
+ 67_310.548_41;
seconds = (seconds * (PI / 180.0) / 240.0) % TWO_PI;
if seconds < 0.0 {
seconds += TWO_PI;
}
seconds
}
/// Returns the propagation error for a mean motion that is not positive.
pub(crate) fn non_positive_mean_motion(mean_motion: f64) -> PropagationError {
PropagationError::new(PropagationErrorKind::InvalidElements)
.with_detail(format!("propagated mean motion is {mean_motion}"))
}

37
crates/orbit/apsis-sgp4/src/state.rs Normal file
View File

@ -0,0 +1,37 @@
//! The propagated state of a satellite.
use apsis_math::Vec3;
/// The position and velocity of a satellite in the True Equator Mean Equinox
/// frame.
///
/// Position is in kilometers and velocity is in kilometers per second. The
/// frame is the one the model defines, with its first axis toward the mean
/// equinox, its third axis along the true equator pole, and the second
/// completing a right-handed set.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct State {
/// Position, in kilometers.
pub position: Vec3,
/// Velocity, in kilometers per second.
pub velocity: Vec3,
/// The propagated semi-major axis, in kilometers.
///
/// This is the semi-major axis after the secular drag term has been
/// applied, so it differs from the value implied by the mean motion in
/// the element set.
pub semi_major_axis: f64,
}
impl State {
/// Returns the distance from the center of the Earth, in kilometers.
pub fn altitude_radius(self) -> f64 {
self.position.norm()
}
/// Returns the speed relative to the Earth's center, in kilometers per
/// second.
pub fn speed(self) -> f64 {
self.velocity.norm()
}
}

51
crates/orbit/apsis-sgp4/tests/model.rs Normal file
View File

@ -0,0 +1,51 @@
//! Tests for the gravity models.
use apsis_sgp4::GravityModel;
#[test]
fn wgs72_defining_values_match_the_standard() {
assert_eq!(GravityModel::WGS72.mu, 398_600.8);
assert_eq!(GravityModel::WGS72.radius, 6_378.135);
assert_eq!(GravityModel::WGS72.j2, 0.001_082_616);
assert_eq!(GravityModel::WGS72.j3, -0.000_002_538_81);
assert_eq!(GravityModel::WGS72.j4, -0.000_001_655_97);
}
#[test]
fn wgs84_defining_values_match_the_standard() {
assert_eq!(GravityModel::WGS84.mu, 398_600.5);
assert_eq!(GravityModel::WGS84.radius, 6_378.137);
assert_eq!(GravityModel::WGS84.j2, 0.001_082_629_89);
}
#[test]
fn kepler_constant_matches_the_reference_value() {
// The reference implementation uses xke = 0.07436691613317342 for WGS
// 72, derived as 60 / sqrt(radius^3 / mu).
let xke = GravityModel::WGS72.xke();
assert!((xke - 0.074_366_916_133_173_42).abs() < 1e-15, "got {xke}");
}
#[test]
fn time_unit_is_the_reciprocal_of_the_kepler_constant() {
let model = GravityModel::WGS72;
assert!((model.tumin() * model.xke() - 1.0).abs() < 1e-15);
// The value is close to 13.44684 minutes for WGS 72.
assert!((model.tumin() - 13.446_839_696_959_309).abs() < 1e-12);
}
#[test]
fn j3_over_j2_matches_the_quotient() {
let model = GravityModel::WGS72;
let expected = model.j3 / model.j2;
assert!((model.j3_over_j2() - expected).abs() < 1e-18);
// The value is about -0.002345069720011528.
assert!((model.j3_over_j2() + 0.002_345_069_720_011_528).abs() < 1e-15);
}
#[test]
fn sqrt_and_reciprocal_sqrt_are_consistent() {
let model = GravityModel::WGS84;
assert!((model.sqrt_mu() * model.reciprocal_sqrt_mu() - 1.0).abs() < 1e-15);
assert!((model.sqrt_mu() - 398_600.5_f64.sqrt()).abs() < 1e-9);
}

317
crates/orbit/apsis-sgp4/tests/propagation.rs Normal file
View File

@ -0,0 +1,317 @@
//! Tests for the propagator against the reference implementation.
//!
//! The expected states in this file were produced by the `sgp4` reference
//! library, which implements Vallado's revision of the theory. The element
//! sets are the historical test cases distributed with that library. A
//! propagation that agrees with the reference to well under a micrometer is
//! treated as correct, because both implementations evaluate the same
//! closed form expressions and should differ only by floating point
//! rounding.
use apsis_sgp4::{GravityModel, PropagationErrorKind, Satellite, State};
use apsis_tle::TwoLineElement;
/// A near circular low Earth orbit, the International Space Station.
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",
);
/// A highly eccentric orbit, which exercises the eccentric anomaly
/// iteration away from circular motion.
const VANGUARD: &str = concat!(
"1 00005U 58002B 00179.78495062 .00000023 00000-0 28098-4 0 4753\r\n",
"2 00005 34.2682 348.7242 1859667 331.7664 19.3264 10.82419157413667\r\n",
);
/// A deep space orbit with a half-day resonance.
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",
);
/// A geostationary transfer orbit, also deep space.
const TRANSFER: &str = concat!(
"1 27651U 03004A 06185.99374200 .00000109 00000-0 19909-4 0 445\r\n",
"2 27651 7.1743 182.9443 7867251 269.0386 18.6601 1.03461646 9884\r\n",
);
/// The position tolerance, in kilometers. The observed agreement is at the
/// nanometer level, so this leaves a wide margin.
const POSITION_TOLERANCE: f64 = 1e-3;
/// The velocity tolerance, in kilometers per second.
const VELOCITY_TOLERANCE: f64 = 1e-9;
fn satellite(text: &str) -> Satellite {
let element = TwoLineElement::parse(text).expect("test element set parses");
Satellite::from_element(&element, GravityModel::WGS72)
.expect("test element set is propagatable")
}
fn assert_matches_reference(state: State, position: [f64; 3], velocity: [f64; 3]) {
let position_error = [
state.position.x - position[0],
state.position.y - position[1],
state.position.z - position[2],
]
.iter()
.fold(0.0_f64, |worst, value| worst.max(value.abs()));
let velocity_error = [
state.velocity.x - velocity[0],
state.velocity.y - velocity[1],
state.velocity.z - velocity[2],
]
.iter()
.fold(0.0_f64, |worst, value| worst.max(value.abs()));
assert!(
position_error < POSITION_TOLERANCE,
"position error {position_error} km, got {}",
state.position
);
assert!(
velocity_error < VELOCITY_TOLERANCE,
"velocity error {velocity_error} km/s, got {}",
state.velocity
);
}
#[test]
fn iss_at_epoch_matches_the_reference() {
let state = satellite(ISS).propagate_minutes(0.0).expect("epoch state");
assert_matches_reference(
state,
[4_083.902_463_521, -993.631_999_606, 5_243.603_665_371],
[2.512_837_295_156, 7.259_888_524_981, -0.583_778_536_506],
);
}
#[test]
fn iss_after_one_minute_matches_the_reference() {
let state = satellite(ISS).propagate_minutes(1.0).expect("state exists");
assert_matches_reference(
state,
[4_224.931_504_689, -556.040_371_394, 5_196.213_688_087],
[2.186_284_834_622, 7.320_793_192_659, -0.995_282_270_726],
);
}
#[test]
fn iss_after_six_hours_matches_the_reference() {
let state = satellite(ISS)
.propagate_minutes(360.0)
.expect("state exists");
assert_matches_reference(
state,
[2_748.401_544_599, -3_564.892_404_578, 4_992.448_308_874],
[4.342_862_050_164, 6.063_045_163_749, 1.927_771_710_260],
);
}
#[test]
fn iss_after_one_day_matches_the_reference() {
let state = satellite(ISS)
.propagate_minutes(1_440.0)
.expect("state exists");
assert_matches_reference(
state,
[-3_199.119_301_995, -5_925.838_895_195, -104.283_883_010],
[4.160_900_126_061, -2.340_866_691_092, 6.034_239_787_489],
);
}
#[test]
fn iss_before_epoch_matches_the_reference() {
let state = satellite(ISS)
.propagate_minutes(-1_440.0)
.expect("state exists");
assert_matches_reference(
state,
[1_121.392_381_234, 6_541.559_708_790, -1_120.952_322_949],
[-4.940_430_025_083, -0.153_942_812_975, -5.902_529_984_919],
);
}
#[test]
fn vanguard_at_epoch_matches_the_reference() {
let state = satellite(VANGUARD)
.propagate_minutes(0.0)
.expect("epoch state");
assert_matches_reference(
state,
[7_022.465_292_664, -1_400.082_967_554, 0.039_951_554],
[1.893_841_014_513, 6.405_893_759_210, 4.534_807_250_355],
);
}
#[test]
fn deep_space_satellite_uses_the_deep_space_branch() {
let deep = satellite(DEEP_SPACE);
assert!(deep.is_deep_space());
let near = satellite(ISS);
assert!(!near.is_deep_space());
}
#[test]
fn deep_space_at_epoch_matches_the_reference() {
// MJD 53_962.36, about 12 hours after the element epoch.
let state = satellite(DEEP_SPACE)
.propagate_minutes(0.0)
.expect("epoch state");
assert_matches_reference(
state,
[12_020.677_224_927, -21_305.501_479_237, -11_902.993_493_655],
[0.646_341_708_451, 1.668_868_685_981, 3.329_292_846_315],
);
}
#[test]
fn deep_space_after_one_week_matches_the_reference() {
let state = satellite(DEEP_SPACE)
.propagate_minutes(10_080.0)
.expect("state exists");
assert_matches_reference(
state,
[12_501.446_941_807, -16_359.501_943_154, -3_718.626_893_065],
[-0.155_492_271_223, 2.911_308_113_855, 3.819_401_898_806],
);
}
#[test]
fn transfer_orbit_at_epoch_matches_the_reference() {
let state = satellite(TRANSFER)
.propagate_minutes(0.0)
.expect("epoch state");
assert_matches_reference(
state,
[-20_728.381_133_505, -10_167.792_006_995, 1_151.382_245_235],
[-1.744_013_319_342, -4.632_866_113_051, 0.570_823_646_148],
);
}
#[test]
fn transfer_orbit_after_one_week_matches_the_reference() {
let state = satellite(TRANSFER)
.propagate_minutes(10_080.0)
.expect("state exists");
assert_matches_reference(
state,
[-14_421.942_143_755, -62_862.455_417_683, 7_869.371_785_647],
[0.946_839_675_888, -1.301_129_854_934, 0.164_515_795_189],
);
}
#[test]
fn propagation_is_deterministic() {
// Repeated evaluation must give the same result, which also verifies
// that the resonance integrator does not carry hidden state between
// calls.
let satellite = satellite(DEEP_SPACE);
let first = satellite.propagate_minutes(5_000.0).expect("state exists");
let second = satellite.propagate_minutes(5_000.0).expect("state exists");
assert_eq!(first, second);
}
#[test]
fn out_of_order_evaluation_gives_the_same_result() {
// Evaluating a later time first must not change a later evaluation of an
// earlier time.
let satellite = satellite(DEEP_SPACE);
let late = satellite.propagate_minutes(10_000.0).expect("state exists");
let early = satellite.propagate_minutes(1_000.0).expect("state exists");
let early_again = satellite.propagate_minutes(1_000.0).expect("state exists");
let late_again = satellite.propagate_minutes(10_000.0).expect("state exists");
assert_eq!(early, early_again);
assert_eq!(late, late_again);
}
#[test]
fn epoch_accessor_returns_the_element_set_epoch() {
let element = TwoLineElement::parse(ISS).expect("parses");
let expected = element
.epoch
.mjd(apsis_time::TimeScale::Utc)
.expect("epoch is in the leap table");
let satellite = Satellite::from_element(&element, GravityModel::WGS72).expect("propagatable");
assert!((satellite.epoch_mjd() - expected).abs() < 1e-12);
}
#[test]
fn minutes_since_epoch_is_measured_in_minutes() {
let element = TwoLineElement::parse(ISS).expect("parses");
let satellite = Satellite::from_element(&element, GravityModel::WGS72).expect("propagatable");
let epoch = satellite.epoch_mjd();
assert!((satellite.minutes_since_epoch(epoch)).abs() < 1e-9);
assert!((satellite.minutes_since_epoch(epoch + 1.0 / 1_440.0) - 1.0).abs() < 1e-6);
assert!((satellite.minutes_since_epoch(epoch - 1.0) + 1_440.0).abs() < 1e-6);
}
#[test]
fn propagate_at_the_epoch_matches_propagate_minutes_zero() {
let element = TwoLineElement::parse(ISS).expect("parses");
let satellite = Satellite::from_element(&element, GravityModel::WGS72).expect("propagatable");
let at_epoch = satellite
.propagate_at(satellite.epoch_mjd())
.expect("state exists");
let at_zero = satellite.propagate_minutes(0.0).expect("state exists");
assert_eq!(at_epoch, at_zero);
}
#[test]
fn invalid_elements_are_rejected() {
// An eccentricity whose leading digit is nine was written with an
// assumed decimal point, so it is 0.9999999, a valid but extremely
// elongated orbit. An eccentricity field that cannot be read at all is
// rejected by the parser instead, which is checked in the TLE tests.
let text = concat!(
"1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927\r\n",
"2 25544 51.6416 247.4627 0000000 130.5360 325.0288 15.72125391563531\r\n",
);
let element = TwoLineElement::parse(text).expect("parses");
// A circular orbit is perfectly valid, so this must succeed.
assert!(Satellite::from_element(&element, GravityModel::WGS72).is_some());
}
#[test]
fn a_zero_mean_motion_is_rejected() {
// The mean motion field is in revolutions per day, so a zero there
// describes no orbit. The element set must still carry a valid
// checksum, so the test builds one from a parsed set rather than
// editing text.
let element = TwoLineElement::parse(ISS).expect("parses");
let mut elements = apsis_sgp4::Elements::from_tle(&element).expect("epoch is covered");
elements.mean_motion = 0.0;
let error = Satellite::try_from_elements(&elements, GravityModel::WGS72)
.expect_err("zero mean motion is not propagatable");
assert_eq!(error.kind(), PropagationErrorKind::InvalidElements);
}
#[test]
fn an_eccentricity_of_one_is_rejected() {
let element = TwoLineElement::parse(ISS).expect("parses");
let mut elements = apsis_sgp4::Elements::from_tle(&element).expect("epoch is covered");
elements.eccentricity = 1.0;
let error = Satellite::try_from_elements(&elements, GravityModel::WGS72)
.expect_err("an open orbit is not propagatable");
assert_eq!(error.kind(), PropagationErrorKind::InvalidElements);
}
#[test]
fn non_finite_elements_are_rejected() {
let element = TwoLineElement::parse(ISS).expect("parses");
let mut elements = apsis_sgp4::Elements::from_tle(&element).expect("epoch is covered");
elements.inclination = f64::NAN;
let error = Satellite::try_from_elements(&elements, GravityModel::WGS72)
.expect_err("a non-finite inclination is not propagatable");
assert_eq!(error.kind(), PropagationErrorKind::InvalidElements);
}
#[test]
fn state_accessors_report_radius_and_speed() {
let state = satellite(ISS).propagate_minutes(0.0).expect("epoch state");
assert!((state.altitude_radius() - state.position.norm()).abs() < 1e-12);
assert!((state.speed() - state.velocity.norm()).abs() < 1e-12);
// A low Earth satellite moves at about 7.7 km/s.
assert!(state.speed() > 7.0 && state.speed() < 8.0);
}