From 1f02cac6a9b85181d528d39bb03f3087b97898ca Mon Sep 17 00:00:00 2001 From: Connor Johnstone Date: Wed, 30 Sep 2026 16:07:54 -0400 Subject: [PATCH] Added equinoctial elements. Also a little cleanup --- benches/kepler.rs | 36 ++- benches/propagation.rs | 6 +- readme.md | 6 +- src/equinoctial.rs | 603 +++++++++++++++++++++++++++++++++++++++++ src/frames.rs | 2 +- src/kepler.rs | 10 +- src/lib.rs | 2 + src/state.rs | 3 +- 8 files changed, 654 insertions(+), 14 deletions(-) create mode 100644 src/equinoctial.rs diff --git a/benches/kepler.rs b/benches/kepler.rs index 58cae6a..1b6ae2f 100644 --- a/benches/kepler.rs +++ b/benches/kepler.rs @@ -8,6 +8,7 @@ use std::hint::black_box; use astrokit::body::Body; use astrokit::cartesian::Cartesian; +use astrokit::equinoctial::Equinoctial; use astrokit::frames::Frame; use astrokit::kepler::{self, Degeneracy, Keplerian}; use criterion::{criterion_group, criterion_main, Criterion}; @@ -121,11 +122,44 @@ fn bench_invariants(c: &mut Criterion) { group.finish(); } +/// Equinoctial against Keplerian on the same states. +/// +/// The interesting comparison is GEO: Keplerian takes its doubly-degenerate +/// branch there and equinoctial takes no branch at all. +fn bench_equinoctial(c: &mut Criterion) { + let mut group = c.benchmark_group("equinoctial"); + + let leo = leo(); + let geo = geo(); + + group.bench_function("from_cartesian/leo", |b| { + b.iter(|| Equinoctial::from_cartesian(&black_box(leo))) + }); + group.bench_function("from_cartesian/geo", |b| { + b.iter(|| Equinoctial::from_cartesian(&black_box(geo))) + }); + + let elements = Equinoctial::from_cartesian(&leo); + group.bench_function("to_cartesian", |b| { + b.iter(|| black_box(elements).to_cartesian()) + }); + group.bench_function("round_trip", |b| { + b.iter(|| Equinoctial::from_cartesian(&black_box(leo)).to_cartesian()) + }); + + let keplerian = leo.keplerian(); + group.bench_function("from_keplerian", |b| { + b.iter(|| Equinoctial::from_keplerian(&black_box(keplerian))) + }); + group.finish(); +} + criterion_group!( benches, bench_scalar_relations, bench_anomalies, bench_conversions, - bench_invariants + bench_invariants, + bench_equinoctial ); criterion_main!(benches); diff --git a/benches/propagation.rs b/benches/propagation.rs index edbda6a..a1d5b61 100644 --- a/benches/propagation.rs +++ b/benches/propagation.rs @@ -117,9 +117,9 @@ fn bench_one_period(c: &mut Criterion) { /// The cost of a state transition matrix relative to a plain propagation. /// -/// Stage 1 measured 10.6x for a 6x6 on a two-body arc. Central differences would -/// need thirteen propagations and lose about half the digits, so anything near -/// that ratio is a win; this tracks it as the force models get heavier. +/// Central differences would need thirteen propagations for a 6x6 and lose about +/// half the available digits. This tracks the ratio as the force models get +/// heavier. #[cfg(feature = "autodiff")] fn bench_state_transition_matrix(c: &mut Criterion) { use nalgebra::{Matrix6, Vector6}; diff --git a/readme.md b/readme.md index 1536b8d..da79a6f 100644 --- a/readme.md +++ b/readme.md @@ -10,9 +10,9 @@ also differentiate. The crate is still being built to a milestone plan: -- **M1** — the study-critical subset: two-body and spherical-harmonic gravity, - drag, third-body, cannonball SRP, GCRF/ITRF frames, ground stations and the - four measurement types, with a validation suite against Orekit reference +- **M1** — the core subset: two-body and spherical-harmonic gravity, drag, + third-body, cannonball SRP, GCRF/ITRF frames, ground stations and the four + measurement types, with a validation suite against Orekit reference trajectories. - **M2** — general-purpose completion: higher-degree gravity fields, more density models, box-wing SRP and eclipses, maneuvers, covariance propagation, diff --git a/src/equinoctial.rs b/src/equinoctial.rs new file mode 100644 index 0000000..3b51773 --- /dev/null +++ b/src/equinoctial.rs @@ -0,0 +1,603 @@ +//! Modified equinoctial elements. +//! +//! Unlike [`Keplerian`], these have no degenerate case on a circular or equatorial +//! orbit. The one singularity is exactly retrograde motion, `i = 180`, which this +//! module refuses rather than returning a diverging `tan(i/2)`. +//! +//! Conversions to and from Cartesian go directly, through the equinoctial basis +//! vectors, not via [`Keplerian`] — routing through the degenerate representation +//! would defeat the purpose. +//! +//! Longitudes are *true*, matching [`Keplerian`]'s use of true anomaly. +//! +//! Walker, Ireland and Owens, "A set of modified equinoctial orbit elements", +//! Celestial Mechanics 36 (1985). + +use nalgebra::Vector3; +use ordinary_diffeq::scalar::Scalar; + +use crate::cartesian::Cartesian; +use crate::frames::Frame; +use crate::kepler::{self, wrap_two_pi, Degeneracy, Keplerian}; +use crate::state::{HasPosition, HasVelocity, State}; + +/// How close to exactly retrograde an orbit may be before the prograde +/// convention is refused. +/// +/// `h` and `k` carry `tan(i/2)`, which diverges at `i = 180`. At this tolerance +/// they reach about 1e7, which is representable but already meaningless. +pub const RETROGRADE_TOLERANCE: f64 = 1e-7; + +/// Modified equinoctial orbital elements. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Equinoctial { + /// Semi-latus rectum `a(1 - e^2)`, m. + pub semi_latus_rectum: T, + /// `e cos(aop + raan)`. + pub f: T, + /// `e sin(aop + raan)`. + pub g: T, + /// `tan(i/2) cos(raan)`. + pub h: T, + /// `tan(i/2) sin(raan)`. + pub k: T, + /// True longitude `raan + aop + true anomaly`, rad. + pub true_longitude: T, + /// The frame these describe motion in. Supplies `mu`. + pub frame: Frame, +} + +impl Equinoctial { + /// Elements from a Cartesian state, without going through [`Keplerian`]. + /// + /// # Panics + /// + /// If the frame is not inertial, or if the orbit is within + /// [`RETROGRADE_TOLERANCE`] of exactly retrograde. + pub fn from_cartesian(state: &Cartesian) -> Self { + state.require_inertial("converting to equinoctial elements"); + let (position, velocity) = (state.position, state.velocity); + let mu = state.frame.mu(); + + let momentum = position.cross(&velocity); + let normal = momentum.normalize(); + + // h and k come straight from the orbit normal: for a prograde orbit the + // normal is [sin i sin raan, -sin i cos raan, cos i], and + // tan(i/2) = sin i / (1 + cos i). + let denominator = T::lift(1.0) + normal[2]; + assert!( + denominator.re() > RETROGRADE_TOLERANCE, + "equinoctial elements in the prograde convention are singular at \ + i = 180 degrees; this orbit is within {RETROGRADE_TOLERANCE:e} of it" + ); + let h = -normal[1] / denominator; + let k = normal[0] / denominator; + + // The semi-latus rectum needs no eccentricity: p = |h|^2 / mu. + let semi_latus_rectum = momentum.norm_squared() / mu; + + // f and g are the eccentricity vector's components in the equinoctial + // basis, and the true longitude is the position's angle in it. + let (f_hat, g_hat) = Self::basis(h, k); + let eccentricity_vector = (position + * (velocity.norm_squared() - T::lift(mu) / position.norm()) + - velocity * position.dot(&velocity)) + / T::lift(mu); + + Self { + semi_latus_rectum, + f: eccentricity_vector.dot(&f_hat), + g: eccentricity_vector.dot(&g_hat), + h, + k, + true_longitude: wrap_two_pi(position.dot(&g_hat).atan2(position.dot(&f_hat))), + frame: state.frame, + } + } + + /// The equinoctial basis vectors, in the frame's axes. + /// + /// `f_hat` is the direction true longitude is measured from and `g_hat` + /// completes the orbital plane. Both are unit vectors and their cross product + /// is the orbit normal. + fn basis(h: T, k: T) -> (Vector3, Vector3) { + let one = T::lift(1.0); + let scale = one / (one + h.powi(2) + k.powi(2)); + ( + Vector3::new( + (one - k.powi(2) + h.powi(2)) * scale, + h * k * 2.0 * scale, + -k * 2.0 * scale, + ), + Vector3::new( + h * k * 2.0 * scale, + (one + k.powi(2) - h.powi(2)) * scale, + h * 2.0 * scale, + ), + ) + } + + /// Checks that the basis is orthonormal and right-handed about the orbit + /// normal implied by `h` and `k`. Used by the tests. + #[cfg(test)] + fn basis_normal(h: T, k: T) -> Vector3 { + let (f_hat, g_hat) = Self::basis(h, k); + f_hat.cross(&g_hat) + } + + /// A Cartesian state, in the frame these elements carry. + pub fn to_cartesian(&self) -> Cartesian { + let one = T::lift(1.0); + let mu = self.frame.mu(); + let (sin_l, cos_l) = (self.true_longitude.sin(), self.true_longitude.cos()); + + let alpha_squared = self.h.powi(2) - self.k.powi(2); + let s_squared = one + self.h.powi(2) + self.k.powi(2); + let w = one + self.f * cos_l + self.g * sin_l; + let radius = self.semi_latus_rectum / w; + let speed_scale = (T::lift(mu) / self.semi_latus_rectum).sqrt(); + let hk2 = self.h * self.k * 2.0; + + let position = Vector3::new( + cos_l + alpha_squared * cos_l + hk2 * sin_l, + sin_l - alpha_squared * sin_l + hk2 * cos_l, + (self.h * sin_l - self.k * cos_l) * 2.0, + ) * (radius / s_squared); + + let velocity = Vector3::new( + -(sin_l + alpha_squared * sin_l - hk2 * cos_l + self.g - hk2 * self.f + + alpha_squared * self.g), + -(-cos_l + alpha_squared * cos_l + hk2 * sin_l - self.f + + hk2 * self.g + + alpha_squared * self.f), + (self.h * cos_l + self.k * sin_l + self.f * self.h + self.g * self.k) * 2.0, + ) * (speed_scale / s_squared); + + Cartesian::new(position, velocity, self.frame) + } + + /// Elements from a Keplerian set, by the defining relations. + /// + /// # Panics + /// + /// If the orbit is within [`RETROGRADE_TOLERANCE`] of exactly retrograde. + pub fn from_keplerian(elements: &Keplerian) -> Self { + assert!( + (core::f64::consts::PI - elements.inclination.re()).abs() > RETROGRADE_TOLERANCE, + "equinoctial elements in the prograde convention are singular at \ + i = 180 degrees" + ); + let longitude_of_periapsis = elements.aop + elements.raan; + let tangent = (elements.inclination * 0.5).tan(); + Self { + semi_latus_rectum: elements.semi_major_axis + * (T::lift(1.0) - elements.eccentricity.powi(2)), + f: elements.eccentricity * longitude_of_periapsis.cos(), + g: elements.eccentricity * longitude_of_periapsis.sin(), + h: tangent * elements.raan.cos(), + k: tangent * elements.raan.sin(), + true_longitude: wrap_two_pi(longitude_of_periapsis + elements.true_anomaly), + frame: elements.frame, + } + } + + /// The Keplerian set these describe. + /// + /// Where the Keplerian set is degenerate the usual substitutions apply and + /// [`Keplerian::degeneracy`] says so, exactly as + /// [`Keplerian::from_cartesian`] would report it. + pub fn to_keplerian(&self) -> Keplerian { + self.to_cartesian().keplerian() + } + + /// Semi-major axis, m. Negative on a hyperbolic orbit. + pub fn semi_major_axis(&self) -> T { + self.semi_latus_rectum / (T::lift(1.0) - self.eccentricity().powi(2)) + } + + /// Eccentricity, dimensionless. + pub fn eccentricity(&self) -> T { + (self.f.powi(2) + self.g.powi(2)).sqrt() + } + + /// Inclination, rad. + pub fn inclination(&self) -> T { + (self.h.powi(2) + self.k.powi(2)).sqrt().atan() * 2.0 + } + + /// Right ascension of the ascending node, rad. Undefined for an equatorial + /// orbit, where `h` and `k` are both zero. + pub fn raan(&self) -> T { + self.k.atan2(self.h) + } + + /// Orbital period, s. + pub fn period(&self) -> T { + kepler::period(self.semi_major_axis(), self.frame.mu()) + } + + /// Mean motion, rad/s. + pub fn mean_motion(&self) -> T { + kepler::mean_motion(self.semi_major_axis(), self.frame.mu()) + } +} + +impl Keplerian { + /// These elements as a non-singular set. + pub fn equinoctial(&self) -> Equinoctial { + Equinoctial::from_keplerian(self) + } +} + +impl Cartesian { + /// The equinoctial elements of this state. + pub fn equinoctial(&self) -> Equinoctial { + Equinoctial::from_cartesian(self) + } +} + +/// Starts a [`State`] from equinoctial elements, for any layout with a position +/// and a velocity. +pub trait FromEquinoctial: HasPosition + HasVelocity { + /// Starts a [`State`] from equinoctial elements. + fn with_equinoctial(self, elements: Equinoctial) -> State { + elements.to_cartesian().to_state(self) + } +} + +impl + HasVelocity> FromEquinoctial for L {} + +/// Whether a Keplerian set's degeneracy makes its angles untrustworthy. +/// +/// A convenience for choosing a representation: if this is true, the equinoctial +/// set is the one to work in. +pub fn prefers_equinoctial(degeneracy: Degeneracy) -> bool { + degeneracy != Degeneracy::None +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::body::Body; + use crate::frames::Orientation; + use crate::state::Orbit; + use approx::assert_relative_eq; + + const MU: f64 = Body::EARTH.mu; + + fn inclined_elliptical() -> Cartesian { + Keplerian { + semi_major_axis: 2.0e7, + eccentricity: 0.3, + inclination: 45.0_f64.to_radians(), + raan: 30.0_f64.to_radians(), + aop: 60.0_f64.to_radians(), + true_anomaly: 120.0_f64.to_radians(), + frame: Frame::EARTH_J2000, + degeneracy: Degeneracy::None, + } + .to_cartesian() + } + + fn circular_equatorial(radius: f64) -> Cartesian { + Cartesian::new( + Vector3::new(radius, 0.0, 0.0), + Vector3::new(0.0, kepler::circular_speed(radius, MU), 0.0), + Frame::EARTH_J2000, + ) + } + + fn leo() -> Cartesian { + Cartesian::new( + Vector3::new( + 4.263868426884883e6, + 5.146189057155391e6, + 1.1310208421331816e6, + ), + Vector3::new(-5923.454461876975, 4496.802639690076, 1870.3893008991558), + Frame::EARTH_J2000, + ) + } + + /// The invariant that was wrong: the basis vectors' cross product has to be + /// the orbit normal that `h` and `k` describe, not some other unit vector. + #[test] + fn test_basis_is_right_handed_about_the_orbit_normal() { + for (inclination, raan) in [ + (45.0_f64, 30.0_f64), + (0.0, 0.0), + (98.0, 200.0), + (170.0, 45.0), + ] { + let (inclination, raan) = (inclination.to_radians(), raan.to_radians()); + let tangent: f64 = (inclination / 2.0).tan(); + let (h, k) = (tangent * raan.cos(), tangent * raan.sin()); + + let expected = Vector3::new( + raan.sin() * inclination.sin(), + -raan.cos() * inclination.sin(), + inclination.cos(), + ); + assert_relative_eq!( + Equinoctial::::basis_normal(h, k), + expected, + max_relative = 1e-12, + epsilon = 1e-14 + ); + } + } + + #[test] + fn test_cartesian_round_trips() { + for state in [inclined_elliptical(), leo(), circular_equatorial(4.2164e7)] { + let back = state.equinoctial().to_cartesian(); + // An equatorial orbit's z components are both about 1e-9 m, so a purely + // relative comparison on them is meaningless. + assert_relative_eq!( + back.position, + state.position, + max_relative = 1e-11, + epsilon = 1e-6 + ); + assert_relative_eq!( + back.velocity, + state.velocity, + max_relative = 1e-11, + epsilon = 1e-9 + ); + assert_eq!(back.frame, state.frame); + } + } + + /// The direct Cartesian conversion has to agree with the route through + /// Keplerian, where Keplerian is trustworthy. + #[test] + fn test_direct_conversion_agrees_with_the_keplerian_route() { + let state = inclined_elliptical(); + let direct = state.equinoctial(); + let via_keplerian = state.keplerian().equinoctial(); + + assert_relative_eq!( + direct.semi_latus_rectum, + via_keplerian.semi_latus_rectum, + max_relative = 1e-11 + ); + assert_relative_eq!( + direct.f, + via_keplerian.f, + max_relative = 1e-11, + epsilon = 1e-12 + ); + assert_relative_eq!( + direct.g, + via_keplerian.g, + max_relative = 1e-11, + epsilon = 1e-12 + ); + assert_relative_eq!( + direct.h, + via_keplerian.h, + max_relative = 1e-11, + epsilon = 1e-12 + ); + assert_relative_eq!( + direct.k, + via_keplerian.k, + max_relative = 1e-11, + epsilon = 1e-12 + ); + assert_relative_eq!( + direct.true_longitude, + via_keplerian.true_longitude, + max_relative = 1e-11 + ); + } + + /// The point of the module: a circular equatorial orbit is a degenerate + /// Keplerian set and an ordinary equinoctial one. + #[test] + fn test_geo_is_degenerate_in_keplerian_and_not_in_equinoctial() { + let state = circular_equatorial(4.2164e7); + + assert_eq!(state.keplerian().degeneracy, Degeneracy::CircularEquatorial); + assert!(prefers_equinoctial(state.keplerian().degeneracy)); + + // Circular and equatorial is just the origin of the f-g-h-k space. + let elements = state.equinoctial(); + assert_relative_eq!(elements.f, 0.0, epsilon = 1e-15); + assert_relative_eq!(elements.g, 0.0, epsilon = 1e-15); + assert_relative_eq!(elements.h, 0.0, epsilon = 1e-15); + assert_relative_eq!(elements.k, 0.0, epsilon = 1e-15); + assert_relative_eq!(elements.semi_latus_rectum, 4.2164e7, max_relative = 1e-12); + assert_relative_eq!(elements.true_longitude, 0.0, epsilon = 1e-12); + } + + #[test] + fn test_derived_quantities_match_the_keplerian_set() { + let state = inclined_elliptical(); + let keplerian = state.keplerian(); + let equinoctial = state.equinoctial(); + + assert_relative_eq!( + equinoctial.semi_major_axis(), + keplerian.semi_major_axis, + max_relative = 1e-11 + ); + assert_relative_eq!( + equinoctial.eccentricity(), + keplerian.eccentricity, + max_relative = 1e-11 + ); + assert_relative_eq!( + equinoctial.inclination(), + keplerian.inclination, + max_relative = 1e-11 + ); + assert_relative_eq!( + equinoctial.raan(), + keplerian.raan, + max_relative = 1e-11, + epsilon = 1e-12 + ); + assert_relative_eq!( + equinoctial.period(), + keplerian.period(), + max_relative = 1e-11 + ); + assert_relative_eq!( + equinoctial.mean_motion(), + keplerian.mean_motion(), + max_relative = 1e-11 + ); + } + + #[test] + fn test_keplerian_round_trips_through_equinoctial() { + let original = inclined_elliptical().keplerian(); + let back = original.equinoctial().to_keplerian(); + + assert_relative_eq!( + back.semi_major_axis, + original.semi_major_axis, + max_relative = 1e-11 + ); + assert_relative_eq!( + back.eccentricity, + original.eccentricity, + max_relative = 1e-11 + ); + assert_relative_eq!(back.inclination, original.inclination, max_relative = 1e-11); + assert_relative_eq!(back.raan, original.raan, max_relative = 1e-11); + assert_relative_eq!( + back.aop, + original.aop, + max_relative = 1e-11, + epsilon = 1e-12 + ); + assert_relative_eq!( + back.true_anomaly, + original.true_anomaly, + max_relative = 1e-11 + ); + } + + #[test] + fn test_state_from_equinoctial() { + let elements = inclined_elliptical().equinoctial(); + let y = Orbit.with_equinoctial(elements); + let back = y.cartesian(Frame::EARTH_J2000).equinoctial(); + assert_relative_eq!(back.f, elements.f, max_relative = 1e-11, epsilon = 1e-12); + assert_relative_eq!( + back.true_longitude, + elements.true_longitude, + max_relative = 1e-11 + ); + } + + #[test] + #[should_panic(expected = "requires an inertial frame")] + fn test_a_rotating_frame_is_refused() { + let rotating = Cartesian::new( + leo().position, + leo().velocity, + Frame::EARTH_J2000.with_orientation(Orientation::BodyFixed), + ); + let _ = rotating.equinoctial(); + } + + #[test] + #[should_panic(expected = "singular at i = 180")] + fn test_exactly_retrograde_is_refused() { + let radius = 7.0e6; + // Equatorial, going the wrong way round: i = 180 exactly. + let retrograde = Cartesian::new( + Vector3::new(radius, 0.0, 0.0), + Vector3::new(0.0, -kepler::circular_speed(radius, MU), 0.0), + Frame::EARTH_J2000, + ); + let _ = retrograde.equinoctial(); + } + + /// A slightly retrograde orbit is fine — only the pole itself is singular. + #[test] + fn test_near_retrograde_is_fine() { + let radius = 7.0e6; + let speed = kepler::circular_speed(radius, MU); + // i = 170 degrees. + let inclination = 170.0_f64.to_radians(); + let state = Cartesian::new( + Vector3::new(radius, 0.0, 0.0), + Vector3::new(0.0, speed * inclination.cos(), speed * inclination.sin()), + Frame::EARTH_J2000, + ); + + let elements = state.equinoctial(); + assert_relative_eq!(elements.inclination(), inclination, max_relative = 1e-10); + let back = elements.to_cartesian(); + assert_relative_eq!(back.position, state.position, max_relative = 1e-10); + } + + /// Equinoctial elements are differentiable where Keplerian angles are not: + /// on a circular orbit, `d(aop)/d(state)` is meaningless, but `d(f)/d(state)` + /// is perfectly well behaved. + #[cfg(feature = "autodiff")] + #[test] + fn test_elements_differentiate_on_a_circular_orbit() { + use nalgebra::{SVector, Vector6}; + use num_dual::DualSVec64; + + let state = circular_equatorial(4.2164e7); + let packed = Vector6::new( + state.position[0], + state.position[1], + state.position[2], + state.velocity[0], + state.velocity[1], + state.velocity[2], + ); + + let (value, gradient): (f64, SVector) = num_dual::gradient( + |y: Vector6>| { + Cartesian::new( + y.fixed_rows::<3>(0).into_owned(), + y.fixed_rows::<3>(3).into_owned(), + Frame::EARTH_J2000, + ) + .equinoctial() + .f + }, + &packed, + ); + + assert_relative_eq!(value, 0.0, epsilon = 1e-15); + // Every partial is finite, which is the property that matters. + for partial in gradient.iter() { + assert!(partial.is_finite(), "gradient has a non-finite entry"); + } + // And it agrees with central differences. + let f_of = |y: Vector6| { + Cartesian::new( + y.fixed_rows::<3>(0).into_owned(), + y.fixed_rows::<3>(3).into_owned(), + Frame::EARTH_J2000, + ) + .equinoctial() + .f + }; + for i in 0..6 { + let h = 1e-4 * packed[i].abs().max(1.0); + let mut plus = packed; + let mut minus = packed; + plus[i] += h; + minus[i] -= h; + let difference = (f_of(plus) - f_of(minus)) / (2.0 * h); + assert_relative_eq!( + gradient[i], + difference, + max_relative = 1e-5, + epsilon = 1e-16 + ); + } + } +} diff --git a/src/frames.rs b/src/frames.rs index 102f80b..b43af95 100644 --- a/src/frames.rs +++ b/src/frames.rs @@ -38,7 +38,7 @@ pub struct Frame { } impl Frame { - /// Earth-centered, J2000-oriented. The frame the study's scenarios work in. + /// Earth-centered, J2000-oriented. pub const EARTH_J2000: Self = Self { center: Body::EARTH, orientation: Orientation::J2000, diff --git a/src/kepler.rs b/src/kepler.rs index 2a13607..26bdb19 100644 --- a/src/kepler.rs +++ b/src/kepler.rs @@ -123,7 +123,7 @@ fn angle_about(from: Vector3, to: Vector3, unit_normal: Vector3 } /// Wraps an angle into `[0, 2pi)`. -fn wrap_two_pi(angle: T) -> T { +pub(crate) fn wrap_two_pi(angle: T) -> T { let turns = (angle.re() / TAU).floor(); angle - TAU * turns } @@ -144,8 +144,8 @@ fn wrap_two_pi(angle: T) -> T { /// - **circular and equatorial** — both are undefined, so `raan` and `aop` are /// zero and `true_anomaly` carries the true longitude /// -/// GEO hits the third case and near-circular LEO the first, so these are the -/// common cases here rather than corner cases. +/// A geostationary orbit hits the third case and a near-circular one the first, +/// so these are common rather than corner cases. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Keplerian { /// Semi-major axis, m. @@ -483,8 +483,8 @@ mod tests { ); } - /// GEO is circular and equatorial, the doubly-singular case, and it is one of - /// the study's scenarios rather than a corner case. + /// A geostationary orbit is circular and equatorial: the doubly-singular case, + /// and a thoroughly ordinary orbit. #[test] fn test_geo_is_reported_as_circular_equatorial() { let radius = 4.2164e7; diff --git a/src/lib.rs b/src/lib.rs index b44f07b..1a1fce7 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -12,6 +12,7 @@ pub mod body; pub mod cartesian; pub mod dynamics; +pub mod equinoctial; pub mod frames; pub mod kepler; pub mod propagation; @@ -25,6 +26,7 @@ pub mod prelude { pub use crate::body::Body; pub use crate::cartesian::{Cartesian, FromCartesian}; pub use crate::dynamics::Dynamics::*; + pub use crate::equinoctial::{Equinoctial, FromEquinoctial}; pub use crate::frames::{Frame, Orientation}; pub use crate::kepler::{Degeneracy, Keplerian, TranslationalLayout}; pub use crate::propagation::{Accuracy, Propagation}; diff --git a/src/state.rs b/src/state.rs index 6728dac..3c015e8 100644 --- a/src/state.rs +++ b/src/state.rs @@ -756,7 +756,8 @@ mod tests { assert_eq!(speed::(&full), 13.0); } - /// The `Orbit + k estimated parameters` case the filters will need. + /// The `Orbit + k estimated parameters` case, for a filter that estimates + /// more than the trajectory. #[test] fn test_user_defined_layout_with_estimated_parameters() { define_state! {