From b089b206f522c1bcbe35d872f011d91c27e48bd3 Mon Sep 17 00:00:00 2001 From: Connor Johnstone Date: Wed, 30 Sep 2026 15:38:04 -0400 Subject: [PATCH] Added full cartesian support with (the beginning of) frames --- benches/kepler.rs | 37 +++--- benches/propagation.rs | 8 +- examples/two_body.rs | 62 +++++---- src/body.rs | 128 ++++++++++++++++++ src/cartesian.rs | 225 ++++++++++++++++++++++++++++++++ src/dynamics/mod.rs | 44 ++++--- src/dynamics/models/two_body.rs | 31 +++-- src/frames.rs | 136 +++++++++++++++++++ src/kepler.rs | 173 ++++++++++++------------ src/lib.rs | 6 + src/propagation.rs | 48 ++++++- 11 files changed, 731 insertions(+), 167 deletions(-) create mode 100644 src/body.rs create mode 100644 src/cartesian.rs create mode 100644 src/frames.rs diff --git a/benches/kepler.rs b/benches/kepler.rs index e50d9f6..58cae6a 100644 --- a/benches/kepler.rs +++ b/benches/kepler.rs @@ -6,31 +6,36 @@ use std::hint::black_box; +use astrokit::body::Body; +use astrokit::cartesian::Cartesian; +use astrokit::frames::Frame; use astrokit::kepler::{self, Degeneracy, Keplerian}; use criterion::{criterion_group, criterion_main, Criterion}; use nalgebra::Vector3; -const MU: f64 = 3.986004415e14; +const MU: f64 = Body::EARTH.mu; /// Inclined, near-circular LEO. Just outside the degeneracy tolerance, so it /// takes the general branch. -fn leo() -> (Vector3, Vector3) { - ( +fn leo() -> Cartesian { + Cartesian::new( Vector3::new( 4.263868426884883e6, 5.146189057155391e6, 1.1310208421331816e6, ), Vector3::new(-5923.454461876975, 4496.802639690076, 1870.3893008991558), + Frame::EARTH_J2000, ) } /// GEO: circular and equatorial, so it takes the doubly-degenerate branch. -fn geo() -> (Vector3, Vector3) { +fn geo() -> Cartesian { let radius = 4.2164e7; - ( + Cartesian::new( Vector3::new(radius, 0.0, 0.0), Vector3::new(0.0, kepler::circular_speed(radius, MU), 0.0), + Frame::EARTH_J2000, ) } @@ -73,14 +78,14 @@ fn bench_anomalies(c: &mut Criterion) { fn bench_conversions(c: &mut Criterion) { let mut group = c.benchmark_group("kepler/conversions"); - let (r, v) = leo(); + let leo = leo(); group.bench_function("from_cartesian/leo", |b| { - b.iter(|| Keplerian::from_cartesian(black_box(r), black_box(v), black_box(MU))) + b.iter(|| Keplerian::from_cartesian(&black_box(leo))) }); - let (r_geo, v_geo) = geo(); + let geo = geo(); group.bench_function("from_cartesian/geo", |b| { - b.iter(|| Keplerian::from_cartesian(black_box(r_geo), black_box(v_geo), black_box(MU))) + b.iter(|| Keplerian::from_cartesian(&black_box(geo))) }); let elements = Keplerian { @@ -90,30 +95,28 @@ fn bench_conversions(c: &mut Criterion) { raan: 0.5, aop: 1.0, true_anomaly: 2.0, + frame: Frame::EARTH_J2000, degeneracy: Degeneracy::None, }; group.bench_function("to_cartesian", |b| { - b.iter(|| black_box(elements).to_cartesian(black_box(MU))) + b.iter(|| black_box(elements).to_cartesian()) }); group.bench_function("round_trip", |b| { - b.iter(|| { - let elements = Keplerian::from_cartesian(black_box(r), black_box(v), MU); - black_box(elements.to_cartesian(MU)) - }) + b.iter(|| Keplerian::from_cartesian(&black_box(leo)).to_cartesian()) }); group.finish(); } fn bench_invariants(c: &mut Criterion) { let mut group = c.benchmark_group("kepler/invariants"); - let (r, v) = leo(); + let leo = leo(); group.bench_function("specific_energy", |b| { - b.iter(|| kepler::specific_energy(black_box(r), black_box(v), black_box(MU))) + b.iter(|| black_box(leo).specific_energy()) }); group.bench_function("specific_angular_momentum", |b| { - b.iter(|| kepler::specific_angular_momentum(black_box(r), black_box(v))) + b.iter(|| black_box(leo).specific_angular_momentum()) }); group.finish(); } diff --git a/benches/propagation.rs b/benches/propagation.rs index c81dfe9..edbda6a 100644 --- a/benches/propagation.rs +++ b/benches/propagation.rs @@ -36,7 +36,7 @@ fn velocity() -> Vector3 { } fn two_body() -> Vec { - vec![TwoBody { mu: MU }] + vec![TwoBody { body: Body::EARTH }] } /// One derivative evaluation, per layout. The layout should cost nothing. @@ -45,7 +45,7 @@ fn bench_derivative(c: &mut Criterion) { let assembly = Assembly::new( dynamics::build::(&two_body()), - Context::new(epoch()), + Context::new(epoch(), Frame::EARTH_J2000), ); let y = Orbit.with_position(position()).with_velocity(velocity()); group.bench_function("two_body/Orbit", |b| { @@ -54,7 +54,7 @@ fn bench_derivative(c: &mut Criterion) { let assembly = Assembly::new( dynamics::build::(&two_body()), - Context::new(epoch()), + Context::new(epoch(), Frame::EARTH_J2000), ); let y = OrbitMass .with_position(position()) @@ -66,7 +66,7 @@ fn bench_derivative(c: &mut Criterion) { let assembly = Assembly::new( dynamics::build::(&two_body()), - Context::new(epoch()), + Context::new(epoch(), Frame::EARTH_J2000), ); let y = Full .with_position(position()) diff --git a/examples/two_body.rs b/examples/two_body.rs index cc7ec25..e80d062 100644 --- a/examples/two_body.rs +++ b/examples/two_body.rs @@ -2,46 +2,38 @@ //! //! Run with `cargo run --example two_body`. -use astrokit::kepler; use astrokit::prelude::*; use nalgebra::Vector3; -/// Earth's gravitational parameter, m^3/s^2. -const MU: f64 = 3.986004415e14; - -/// Semi-major axis of the test orbit, m. -const SMA: f64 = 6.7781363e6; - -/// Earth's equatorial radius, m. Comes from a body once bodies land. -const EARTH_RADIUS: f64 = 6.378137e6; - fn main() { let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0); - let period = Duration::from_seconds(kepler::period(SMA, MU)); + let frame = Frame::EARTH_J2000; - let y0 = Orbit - .with_position(Vector3::new( + let initial = Cartesian::new( + Vector3::new( 4.263868426884883e6, 5.146189057155391e6, 1.1310208421331816e6, - )) - .with_velocity(Vector3::new( - -5923.454461876975, - 4496.802639690076, - 1870.3893008991558, - )); + ), + Vector3::new(-5923.454461876975, 4496.802639690076, 1870.3893008991558), + frame, + ); - let forces = vec![TwoBody { mu: MU }]; + let elements = initial.keplerian(); + let period = Duration::from_seconds(initial.period()); + let forces = vec![TwoBody { body: Body::EARTH }]; - let mut propagation = Propagation::new(epoch, y0, &forces) + let mut propagation = Propagation::from_cartesian(epoch, initial, &forces) .duration(period) .accuracy(Accuracy::High); let solution = propagation.run(); - let end = Orbit.state(*solution.states.last().unwrap()); + let end = Orbit + .state(*solution.states.last().unwrap()) + .cartesian(frame); - let elements = y0.keplerian(MU); println!("epoch: {epoch}"); + println!("body: {}", initial.frame.center.name()); println!("period: {period}"); println!( "elements: a = {:.3} km, e = {:.6}, i = {:.3} deg ({:?})", @@ -50,24 +42,30 @@ fn main() { elements.inclination.to_degrees(), elements.degeneracy ); + println!( + "altitude: {:.3} km", + initial.spherical_altitude() / 1e3 + ); println!("steps: {}", solution.states.len() - 1); println!( "position error: {:.3e} m", - (end.position() - y0.position()).norm() + (end.position - initial.position).norm() ); println!( "velocity error: {:.3e} m/s", - (end.velocity() - y0.velocity()).norm() + (end.velocity - initial.velocity).norm() + ); + println!( + "energy drift: {:.3e} J/kg", + end.specific_energy() - initial.specific_energy() ); // Dense output: the state at any time in the arc, not only at a step. - let half = Orbit.state(solution.interpolate(period.to_seconds() / 2.0)); + let half = Orbit + .state(solution.interpolate(period.to_seconds() / 2.0)) + .cartesian(frame); println!( - "altitude at half a period: {:.3} km", - (half.position().norm() - EARTH_RADIUS) / 1e3 - ); - println!( - "specific energy drift: {:.3e} J/kg", - end.specific_energy(MU) - y0.specific_energy(MU) + "half a period: {:.3} km altitude", + half.spherical_altitude() / 1e3 ); } diff --git a/src/body.rs b/src/body.rs new file mode 100644 index 0000000..83d44fa --- /dev/null +++ b/src/body.rs @@ -0,0 +1,128 @@ +//! Celestial bodies. + +use serde::{Deserialize, Serialize}; + +/// A celestial body: the gravitational parameter and shape that force models and +/// element conversions need. +/// +/// Carried by a [`Frame`](crate::frames::Frame) and by the dynamics that need a +/// central body, so that no signature takes a bare `mu`. A propagation about the +/// Moon is then a matter of saying so rather than of working around an assumption. +#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)] +pub struct Body { + /// NAIF integer id, for ANISE ephemeris and orientation lookups. + pub naif_id: i32, + /// Gravitational parameter, m^3/s^2. + pub mu: f64, + /// Equatorial radius, m. + pub equatorial_radius: f64, + /// Flattening, dimensionless. Zero for a sphere. + pub flattening: f64, +} + +impl Body { + /// Earth. + /// + /// `mu` is EGM96's. **It must agree with whatever gravity field a + /// propagation uses**: the point-mass term and the spherical harmonics are two + /// halves of one expansion, and taking `mu` from one model and the harmonics + /// from another is an inconsistency. Once a gravity + /// field is loaded it should supply its own `mu` and override this. + /// + /// Radius and flattening are WGS84. + pub const EARTH: Self = Self { + naif_id: 399, + mu: 3.986004415e14, + equatorial_radius: 6378137.0, + flattening: 1.0 / 298.257223563, + }; + + /// The Moon. `mu` and radius from DE440; the figure is very nearly spherical. + pub const MOON: Self = Self { + naif_id: 301, + mu: 4.902800118e12, + equatorial_radius: 1737400.0, + flattening: 0.0012, + }; + + /// The Sun, treated as spherical. + pub const SUN: Self = Self { + naif_id: 10, + mu: 1.32712440018e20, + equatorial_radius: 6.957e8, + flattening: 0.0, + }; + + /// A spherical body, for a test case or a body whose figure does not matter. + pub const fn spherical(naif_id: i32, mu: f64, radius: f64) -> Self { + Self { + naif_id, + mu, + equatorial_radius: radius, + flattening: 0.0, + } + } + + /// Polar radius, m. + pub fn polar_radius(&self) -> f64 { + self.equatorial_radius * (1.0 - self.flattening) + } + + /// A name for the recognized bodies, or a rendering of the NAIF id. + /// + /// Diagnostics only. Provenance records the fields, which is the thing that + /// actually pins a run down. + pub fn name(&self) -> String { + match self.naif_id { + 10 => "Sun".to_string(), + 199 => "Mercury".to_string(), + 299 => "Venus".to_string(), + 301 => "Moon".to_string(), + 399 => "Earth".to_string(), + 499 => "Mars".to_string(), + 599 => "Jupiter".to_string(), + 699 => "Saturn".to_string(), + 799 => "Uranus".to_string(), + 899 => "Neptune".to_string(), + id => format!("NAIF {id}"), + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_earth_matches_wgs84() { + assert_eq!(Body::EARTH.equatorial_radius, 6378137.0); + // WGS84's polar radius, to the millimeter. + assert_relative_eq!(Body::EARTH.polar_radius(), 6356752.314245, epsilon = 1e-3); + assert_eq!(Body::EARTH.name(), "Earth"); + } + + #[test] + fn test_spherical_body_has_equal_radii() { + let body = Body::spherical(9999, 1.0e14, 1.0e6); + assert_eq!(body.polar_radius(), body.equatorial_radius); + assert_eq!(body.name(), "NAIF 9999"); + } + + /// Bodies go into a run's provenance, so they have to round-trip. + #[test] + fn test_body_round_trips_through_json() { + let json = serde_json::to_string(&Body::EARTH).unwrap(); + assert_eq!(serde_json::from_str::(&json).unwrap(), Body::EARTH); + } + + /// The gravitational parameters should agree with ANISE's, which come from + /// JPL kernels. ANISE works in km^3/s^2. + #[test] + fn test_gravitational_parameters_are_the_expected_magnitude() { + // Earth: EGM96 against the DE440 value, which differ in the 9th digit. + assert_relative_eq!(Body::EARTH.mu * 1e-9, 398600.4415, max_relative = 1e-10); + assert_relative_eq!(Body::MOON.mu * 1e-9, 4902.800118, max_relative = 1e-10); + assert_relative_eq!(Body::SUN.mu * 1e-9, 132712440018.0, max_relative = 1e-10); + } +} diff --git a/src/cartesian.rs b/src/cartesian.rs new file mode 100644 index 0000000..bc0264f --- /dev/null +++ b/src/cartesian.rs @@ -0,0 +1,225 @@ +//! A translational state as position, velocity and frame. +//! +//! This is the interchange representation. An integrated +//! [`State`] is layout-shaped and may carry mass, attitude or +//! estimated parameters; a `Cartesian` is the translational part of it together +//! with the frame it is expressed in, which is what every element conversion needs +//! and what no bare pair of vectors can supply. + +use nalgebra::Vector3; +use ordinary_diffeq::scalar::Scalar; + +use crate::frames::Frame; +use crate::kepler::{self, Keplerian}; +use crate::state::{HasPosition, HasVelocity, State}; + +/// A position and velocity in a reference frame. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Cartesian { + /// Position, m. + pub position: Vector3, + /// Velocity, m/s. + pub velocity: Vector3, + /// The frame both are expressed in. + pub frame: Frame, +} + +impl Cartesian { + /// A state from a position and velocity in `frame`. + pub fn new(position: Vector3, velocity: Vector3, frame: Frame) -> Self { + Self { + position, + velocity, + frame, + } + } + + /// Distance from the frame's origin, m. + pub fn radius(&self) -> T { + self.position.norm() + } + + /// Speed, m/s. + pub fn speed(&self) -> T { + self.velocity.norm() + } + + /// Height above the central body's equatorial radius, m. + /// + /// Spherical, not geodetic: it subtracts the equatorial radius rather than + /// solving for the distance to the reference ellipsoid. For the Earth the two + /// differ by up to about 21 km at the poles. Geodetic altitude arrives with + /// the frame transformations, which is also where it belongs, since it is only + /// meaningful in a body-fixed frame. + pub fn spherical_altitude(&self) -> T { + self.radius() - self.frame.center.equatorial_radius + } + + /// Specific orbital energy, J/kg. + pub fn specific_energy(&self) -> T { + kepler::specific_energy(self.position, self.velocity, self.frame.mu()) + } + + /// Specific angular momentum, m^2/s. + pub fn specific_angular_momentum(&self) -> Vector3 { + kepler::specific_angular_momentum(self.position, self.velocity) + } + + /// The Keplerian elements of this state. + /// + /// # Panics + /// + /// If the frame is not inertial. Orbital elements describe inertial motion, and + /// a velocity in a rotating frame is not the inertial velocity — converting one + /// anyway would return a plausible-looking wrong answer. + pub fn keplerian(&self) -> Keplerian { + Keplerian::from_cartesian(self) + } + + /// Orbital period, s. + pub fn period(&self) -> T { + kepler::period(self.keplerian().semi_major_axis, self.frame.mu()) + } + + /// This state in a layout, with everything the layout also carries left zero. + pub fn to_state(&self, layout: L) -> State + where + L: HasPosition + HasVelocity, + { + layout + .with_position(self.position) + .with_velocity(self.velocity) + } + + /// Asserts the frame is inertial, for the operations that require it. + pub(crate) fn require_inertial(&self, what: &str) { + assert!( + self.frame.is_inertial(), + "{what} requires an inertial frame, but this state is in {:?} about {}", + self.frame.orientation, + self.frame.center.name() + ); + } +} + +impl + HasVelocity> State { + /// The translational part of this state, in `frame`. + /// + /// A layout carries no frame of its own — a propagation holds that — so it has + /// to be supplied here. + pub fn cartesian(&self, frame: Frame) -> Cartesian { + Cartesian::new(self.position(), self.velocity(), frame) + } +} + +/// Layouts that can be filled from a [`Cartesian`]. +pub trait FromCartesian: HasPosition + HasVelocity { + /// Starts a [`State`] from a Cartesian state, discarding its frame. + fn with_cartesian(self, state: &Cartesian) -> State { + state.to_state(self) + } +} + +impl + HasVelocity> FromCartesian for L {} + +#[cfg(test)] +mod tests { + use super::*; + use crate::body::Body; + use crate::frames::Orientation; + use crate::state::{Full, Orbit, OrbitMass}; + use approx::assert_relative_eq; + + fn leo() -> Cartesian { + Cartesian::new( + Vector3::new( + 4.263868426884883e6, + 5.146189057155391e6, + 1.1310208421331816e6, + ), + Vector3::new(-5923.454461876975, 4496.802639690076, 1870.3893008991558), + Frame::EARTH_J2000, + ) + } + + #[test] + fn test_radius_speed_and_altitude() { + let state = leo(); + assert_relative_eq!(state.radius(), 6.778135e6, max_relative = 1e-5); + assert_relative_eq!(state.speed(), 7668.5, max_relative = 1e-3); + // A 400 km circular orbit, measured off the equatorial radius. + assert_relative_eq!(state.spherical_altitude(), 399998.0, max_relative = 1e-3); + } + + #[test] + fn test_invariants_read_mu_from_the_frame() { + let state = leo(); + assert_relative_eq!( + state.specific_energy(), + kepler::specific_energy(state.position, state.velocity, Body::EARTH.mu), + max_relative = 1e-15 + ); + assert_eq!( + state.specific_angular_momentum(), + state.position.cross(&state.velocity) + ); + } + + /// The same Cartesian state about a different body gives different elements, + /// because the frame supplies `mu`. This is the whole reason for the type. + #[test] + fn test_the_frame_changes_the_elements() { + let about_earth = leo(); + let about_moon = Cartesian::new( + about_earth.position, + about_earth.velocity, + Frame::MOON_J2000, + ); + + let earth = about_earth.keplerian(); + let moon = about_moon.keplerian(); + + // Same geometry, different central mass, so a different semi-major axis. + assert_relative_eq!(earth.inclination, moon.inclination, max_relative = 1e-12); + assert!(earth.semi_major_axis != moon.semi_major_axis); + // Far too fast for the Moon to hold on to: hyperbolic, so sma is negative. + assert!(moon.semi_major_axis < 0.0); + } + + #[test] + fn test_round_trips_through_a_layout() { + let state = leo(); + + for_each_layout(&state); + + fn for_each_layout(state: &Cartesian) { + let on_orbit = state.to_state(Orbit); + assert_eq!(on_orbit.position(), state.position); + assert_eq!(on_orbit.cartesian(state.frame), *state); + + let on_mass = Orbit.with_cartesian(state); + assert_eq!(on_mass.velocity(), state.velocity); + + // A wider layout takes the translational part and leaves the rest zero. + let on_full = state.to_state(Full); + assert_eq!(on_full.position(), state.position); + assert_eq!(on_full.mass(), 0.0); + assert_eq!(on_full.cartesian(state.frame), *state); + + let on_orbit_mass = OrbitMass.with_cartesian(state).with_mass(500.0); + assert_eq!(on_orbit_mass.mass(), 500.0); + assert_eq!(on_orbit_mass.cartesian(state.frame), *state); + } + } + + #[test] + #[should_panic(expected = "requires an inertial frame")] + fn test_elements_from_a_rotating_frame_are_refused() { + let rotating = Cartesian::new( + leo().position, + leo().velocity, + Frame::EARTH_J2000.with_orientation(Orientation::BodyFixed), + ); + let _ = rotating.keplerian(); + } +} diff --git a/src/dynamics/mod.rs b/src/dynamics/mod.rs index c92fa74..fd68eaa 100644 --- a/src/dynamics/mod.rs +++ b/src/dynamics/mod.rs @@ -16,7 +16,7 @@ //! use astrokit::dynamics::{build, Dynamics}; //! use astrokit::state::Orbit; //! -//! let spec = vec![Dynamics::TwoBody { mu: 3.986004415e14 }]; +//! let spec = vec![Dynamics::TwoBody { body: astrokit::body::Body::EARTH }]; //! //! // The same spec, at two scalars. //! let plain = build::(&spec); @@ -29,8 +29,9 @@ //! hidden behind a mode flag. //! //! A side benefit: because a model is constructed at the propagated scalar, it -//! stores its constants already converted. [`models::TwoBody`] holds `mu` as a -//! `T`, so `T::lift` happens once at build time and never in the derivative. +//! stores its constants already converted. [`models::TwoBody`] holds its body's +//! `mu` as a `T`, so `T::lift` happens once at build time and never in the +//! derivative. //! //! # Force or acceleration //! @@ -47,6 +48,8 @@ use ordinary_diffeq::scalar::Scalar; use serde::{Deserialize, Serialize}; use std::ops::{Add, AddAssign}; +use crate::body::Body; +use crate::frames::Frame; use crate::state::{HasPosition, HasVelocity, Layout, State}; pub mod models; @@ -124,17 +127,20 @@ impl Add for Contribution { } /// What a [`DynamicsModel`] needs that is not in the state vector. -#[derive(Clone, Copy, Debug, PartialEq, Eq)] +#[derive(Clone, Copy, Debug, PartialEq)] pub struct Context { /// The epoch the propagation started from. The `t` a model receives is /// seconds after this. pub epoch: Epoch, + /// The frame the state is expressed in, and so the central body a model can + /// assume without being told. + pub frame: Frame, } impl Context { - /// A context anchored at `epoch`. - pub fn new(epoch: Epoch) -> Self { - Self { epoch } + /// A context anchored at `epoch`, in `frame`. + pub fn new(epoch: Epoch, frame: Frame) -> Self { + Self { epoch, frame } } /// The epoch `t` seconds after the start of the propagation. @@ -179,10 +185,13 @@ pub trait DynamicsModel> { #[non_exhaustive] #[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)] pub enum Dynamics { - /// Point-mass gravity of a central body. + /// Point-mass gravity of a body. + /// + /// Named explicitly rather than taken from the propagation's frame, because an + /// n-body composition needs several and only one of them is the center. TwoBody { - /// Gravitational parameter, m^3/s^2. - mu: f64, + /// The attracting body. + body: Body, }, } @@ -194,7 +203,7 @@ impl Dynamics { L: HasPosition + HasVelocity + 'static, { match *self { - Self::TwoBody { mu } => Box::new(models::TwoBody::new(mu)), + Self::TwoBody { body } => Box::new(models::TwoBody::new(body)), } } } @@ -246,7 +255,7 @@ mod tests { #[test] fn test_spec_round_trips_through_json() { - let spec = vec![Dynamics::TwoBody { mu: 3.986004415e14 }]; + let spec = vec![Dynamics::TwoBody { body: Body::EARTH }]; let json = serde_json::to_string(&spec).unwrap(); let back: Vec = serde_json::from_str(&json).unwrap(); assert_eq!(spec, back); @@ -256,7 +265,7 @@ mod tests { /// built at two different scalars. #[test] fn test_one_spec_builds_at_two_scalars() { - let spec = vec![Dynamics::TwoBody { mu: 3.986004415e14 }]; + let spec = vec![Dynamics::TwoBody { body: Body::EARTH }]; let plain = build::(&spec); assert_eq!(plain.len(), 1); @@ -272,7 +281,7 @@ mod tests { #[test] fn test_context_advances_an_epoch() { let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0); - let context = Context::new(epoch); + let context = Context::new(epoch, Frame::EARTH_J2000); assert_eq!(context.epoch_at(0.0), epoch); assert_eq!( context.epoch_at(3600.0), @@ -283,11 +292,14 @@ mod tests { /// A model bounded on a capability works for any layout that has it. #[test] fn test_two_body_acceleration_is_layout_agnostic() { - let context = Context::new(Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0)); + let context = Context::new( + Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0), + Frame::EARTH_J2000, + ); let mu = 3.986004415e14; let r = Vector3::new(7.0e6, 0.0, 0.0); - let model = models::TwoBody::::new(mu); + let model = models::TwoBody::::new(Body::EARTH); let on_orbit = Orbit.with_position(r); let on_full = Full.with_position(r).with_mass(500.0); diff --git a/src/dynamics/models/two_body.rs b/src/dynamics/models/two_body.rs index 72ad0c9..c103d2c 100644 --- a/src/dynamics/models/two_body.rs +++ b/src/dynamics/models/two_body.rs @@ -3,22 +3,31 @@ use ordinary_diffeq::scalar::Scalar; use super::super::{Context, Contribution, DynamicsModel}; +use crate::body::Body; use crate::state::{HasPosition, State}; -/// Point-mass gravity of a single central body, at the origin of the state's -/// frame. +/// Point-mass gravity of a single body, at the origin of the state's frame. /// -/// `mu` is stored in the propagated scalar, converted once at construction, so -/// the derivative does no lifting. +/// The body's `mu` is stored in the propagated scalar, converted once at +/// construction, so the derivative does no lifting. #[derive(Clone, Copy, Debug, PartialEq)] pub struct TwoBody { + body: Body, mu: T, } impl TwoBody { - /// Point-mass gravity with gravitational parameter `mu`, in m^3/s^2. - pub fn new(mu: f64) -> Self { - Self { mu: T::lift(mu) } + /// Point-mass gravity of `body`. + pub fn new(body: Body) -> Self { + Self { + body, + mu: T::lift(body.mu), + } + } + + /// The attracting body. + pub fn body(&self) -> Body { + self.body } /// The gravitational parameter, m^3/s^2. @@ -42,6 +51,7 @@ impl> DynamicsModel for Tw #[cfg(test)] mod tests { use super::*; + use crate::frames::Frame; use crate::state::Orbit; use approx::assert_relative_eq; use hifitime::Epoch; @@ -51,8 +61,11 @@ mod tests { #[test] fn test_acceleration_points_at_the_origin() { - let context = Context::new(Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0)); - let model = TwoBody::::new(MU); + let context = Context::new( + Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0), + Frame::EARTH_J2000, + ); + let model = TwoBody::::new(Body::EARTH); let y = Orbit.with_position(Vector3::new(3.0e6, 4.0e6, 0.0)); let c = model.contribution(0.0, &y, &context); diff --git a/src/frames.rs b/src/frames.rs new file mode 100644 index 0000000..102f80b --- /dev/null +++ b/src/frames.rs @@ -0,0 +1,136 @@ +//! Reference frames: where a state's origin is and how its axes are oriented. +//! +//! This module carries the frame *type*. Transformations between orientations — +//! J2000 to ITRF, topocentric, geodetic site coordinates — are a separate piece of +//! work and will land here on top of ANISE. + +use serde::{Deserialize, Serialize}; + +use crate::body::Body; + +/// How a frame's axes are oriented. +/// +/// Provisional and deliberately coarse. The distinctions that matter for +/// sub-arcsecond work — ICRF against EME2000, ITRF against a body-fixed IAU frame +/// — arrive with the transformations, which is the point at which getting them +/// wrong would show up. Marked `#[non_exhaustive]` so refining it is not a +/// breaking change. +#[non_exhaustive] +#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)] +pub enum Orientation { + /// Inertial, aligned with the J2000 mean equator and equinox. + J2000, + /// Rotating with the central body. + BodyFixed, +} + +/// A reference frame: a central body and an orientation. +/// +/// Held by every [`Cartesian`](crate::cartesian::Cartesian) and every element set, +/// which is what lets conversions read `mu` from the state rather than take it as +/// an argument. +#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)] +pub struct Frame { + /// The body at the origin. + pub center: Body, + /// How the axes are oriented. + pub orientation: Orientation, +} + +impl Frame { + /// Earth-centered, J2000-oriented. The frame the study's scenarios work in. + pub const EARTH_J2000: Self = Self { + center: Body::EARTH, + orientation: Orientation::J2000, + }; + + /// Earth-centered, rotating with the Earth. + pub const EARTH_FIXED: Self = Self { + center: Body::EARTH, + orientation: Orientation::BodyFixed, + }; + + /// Moon-centered, J2000-oriented. + pub const MOON_J2000: Self = Self { + center: Body::MOON, + orientation: Orientation::J2000, + }; + + /// A J2000-oriented frame centered on `center`. + pub const fn inertial(center: Body) -> Self { + Self { + center, + orientation: Orientation::J2000, + } + } + + /// A frame rotating with `center`. + pub const fn body_fixed(center: Body) -> Self { + Self { + center, + orientation: Orientation::BodyFixed, + } + } + + /// The central body's gravitational parameter, m^3/s^2. + pub fn mu(&self) -> f64 { + self.center.mu + } + + /// Whether orbital motion in this frame is unaccelerated by the frame itself. + /// + /// Element conversions and the two-body relations are only meaningful in an + /// inertial frame; in a rotating one the velocity is not the inertial velocity. + pub fn is_inertial(&self) -> bool { + matches!(self.orientation, Orientation::J2000) + } + + /// The same frame with a different orientation. + pub fn with_orientation(&self, orientation: Orientation) -> Self { + Self { + orientation, + ..*self + } + } +} + +impl Default for Frame { + fn default() -> Self { + Self::EARTH_J2000 + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_frame_carries_the_central_bodys_mu() { + assert_eq!(Frame::EARTH_J2000.mu(), Body::EARTH.mu); + assert_eq!(Frame::MOON_J2000.mu(), Body::MOON.mu); + } + + #[test] + fn test_only_inertial_frames_report_inertial() { + assert!(Frame::EARTH_J2000.is_inertial()); + assert!(!Frame::EARTH_FIXED.is_inertial()); + assert!(Frame::inertial(Body::SUN).is_inertial()); + assert!(!Frame::body_fixed(Body::SUN).is_inertial()); + } + + #[test] + fn test_orientation_can_be_swapped_keeping_the_center() { + let fixed = Frame::EARTH_J2000.with_orientation(Orientation::BodyFixed); + assert_eq!(fixed, Frame::EARTH_FIXED); + assert_eq!(fixed.center, Body::EARTH); + } + + #[test] + fn test_frame_round_trips_through_json() { + let json = serde_json::to_string(&Frame::MOON_J2000).unwrap(); + assert_eq!( + serde_json::from_str::(&json).unwrap(), + Frame::MOON_J2000 + ); + } +} diff --git a/src/kepler.rs b/src/kepler.rs index b5282b5..2a13607 100644 --- a/src/kepler.rs +++ b/src/kepler.rs @@ -15,6 +15,8 @@ use core::f64::consts::TAU; use nalgebra::Vector3; use ordinary_diffeq::scalar::Scalar; +use crate::cartesian::Cartesian; +use crate::frames::Frame; use crate::state::{HasPosition, HasVelocity, State}; /// How close to zero an eccentricity or inclination has to be before @@ -158,6 +160,9 @@ pub struct Keplerian { pub aop: T, /// True anomaly, rad. pub true_anomaly: T, + /// The frame these describe motion in. Supplies `mu`, so no method here takes + /// one. + pub frame: Frame, /// Which elements, if any, are not what they are named. pub degeneracy: Degeneracy, } @@ -176,11 +181,17 @@ pub enum Degeneracy { } impl Keplerian { - /// Elements from an inertial position and velocity. + /// Elements from a Cartesian state. /// - /// Vallado's RV2COE, with the three singular cases handled as - /// [`Keplerian`] documents. - pub fn from_cartesian(position: Vector3, velocity: Vector3, mu: f64) -> Self { + /// Vallado's RV2COE, with the three singular cases handled as [`Keplerian`] + /// documents and the angles measured by `atan2` rather than `acos`. + /// + /// # Panics + /// + /// If the state's frame is not inertial. + pub fn from_cartesian(state: &Cartesian) -> Self { + state.require_inertial("converting to Keplerian elements"); + let (position, velocity, mu) = (state.position, state.velocity, state.frame.mu()); let zero = T::lift(0.0); let r = position.norm(); @@ -265,17 +276,19 @@ impl Keplerian { raan, aop, true_anomaly, + frame: state.frame, degeneracy, } } - /// An inertial position and velocity, in meters and meters per second. + /// A Cartesian state, in the frame these elements carry. /// - /// Inverts [`Keplerian::from_cartesian`] for a non-degenerate set. A - /// degenerate set round-trips too, because the substituted angles are zero - /// where they are undefined and the perifocal construction then puts the - /// argument of latitude or true longitude in the right place. - pub fn to_cartesian(&self, mu: f64) -> (Vector3, Vector3) { + /// Inverts [`Keplerian::from_cartesian`] for a non-degenerate set. A degenerate + /// set round-trips too, because the substituted angles are zero where they are + /// undefined and the perifocal construction then puts the argument of latitude + /// or true longitude in the right place. + pub fn to_cartesian(&self) -> Cartesian { + let mu = self.frame.mu(); let one = T::lift(1.0); let semi_latus = self.semi_major_axis * (one - self.eccentricity.powi(2)); let radius = semi_latus / (one + self.eccentricity * self.true_anomaly.cos()); @@ -293,9 +306,10 @@ impl Keplerian { T::lift(0.0), ); - ( + Cartesian::new( self.perifocal_to_inertial(in_plane), self.perifocal_to_inertial(velocity_in_plane), + self.frame, ) } @@ -319,13 +333,13 @@ impl Keplerian { } /// Orbital period, in seconds. - pub fn period(&self, mu: f64) -> T { - period(self.semi_major_axis, mu) + pub fn period(&self) -> T { + period(self.semi_major_axis, self.frame.mu()) } /// Mean motion, in radians per second. - pub fn mean_motion(&self, mu: f64) -> T { - mean_motion(self.semi_major_axis, mu) + pub fn mean_motion(&self) -> T { + mean_motion(self.semi_major_axis, self.frame.mu()) } /// Eccentric anomaly, in radians. @@ -355,31 +369,20 @@ impl Keplerian { /// these without doing anything. pub trait TranslationalLayout: HasPosition + HasVelocity { /// Starts a [`State`] from Keplerian elements. - fn with_keplerian(self, elements: Keplerian, mu: f64) -> State { - let (position, velocity) = elements.to_cartesian(mu); - self.with_position(position).with_velocity(velocity) + fn with_keplerian(self, elements: Keplerian) -> State { + elements.to_cartesian().to_state(self) } } impl + HasVelocity> TranslationalLayout for L {} impl + HasVelocity> State { - /// The Keplerian elements of this state, about a body of parameter `mu`. - pub fn keplerian(&self, mu: f64) -> Keplerian { - Keplerian::from_cartesian(self.position(), self.velocity(), mu) + /// The Keplerian elements of this state, in `frame`. + pub fn keplerian(&self, frame: Frame) -> Keplerian { + self.cartesian(frame).keplerian() } - /// Orbital period, in seconds. - pub fn period(&self, mu: f64) -> T { - period(self.keplerian(mu).semi_major_axis, mu) - } - - /// Specific orbital energy, in J/kg. - pub fn specific_energy(&self, mu: f64) -> T { - specific_energy(self.position(), self.velocity(), mu) - } - - /// Specific angular momentum, in m^2/s. + /// Specific angular momentum, in m^2/s. Needs no frame. pub fn specific_angular_momentum(&self) -> Vector3 { specific_angular_momentum(self.position(), self.velocity()) } @@ -388,21 +391,33 @@ impl + HasVelocity> State (Vector3, Vector3) { - ( + fn leo() -> Cartesian { + Cartesian::new( Vector3::new( 4.263868426884883e6, 5.146189057155391e6, 1.1310208421331816e6, ), Vector3::new(-5923.454461876975, 4496.802639690076, 1870.3893008991558), + Frame::EARTH_J2000, + ) + } + + /// A circular orbit of the given radius in the equatorial plane. + fn circular_equatorial(radius: f64) -> Cartesian { + Cartesian::new( + Vector3::new(radius, 0.0, 0.0), + Vector3::new(0.0, circular_speed(radius, MU), 0.0), + Frame::EARTH_J2000, ) } @@ -418,23 +433,18 @@ mod tests { #[test] fn test_circular_speed_gives_a_closed_orbit() { let radius = 7.0e6; - let speed = circular_speed(radius, MU); - let elements = Keplerian::from_cartesian( - Vector3::new(radius, 0.0, 0.0), - Vector3::new(0.0, speed, 0.0), - MU, - ); + let elements = circular_equatorial(radius).keplerian(); assert_relative_eq!(elements.semi_major_axis, radius, max_relative = 1e-12); assert!(elements.eccentricity < 1e-15); } #[test] fn test_cartesian_round_trips_through_elements() { - let (r, v) = leo(); - let elements = Keplerian::from_cartesian(r, v, MU); - let (r2, v2) = elements.to_cartesian(MU); - assert_relative_eq!(r2, r, max_relative = 1e-12); - assert_relative_eq!(v2, v, max_relative = 1e-12); + let state = leo(); + let back = state.keplerian().to_cartesian(); + assert_relative_eq!(back.position, state.position, max_relative = 1e-12); + assert_relative_eq!(back.velocity, state.velocity, max_relative = 1e-12); + assert_eq!(back.frame, state.frame); } /// An inclined elliptical orbit, where no element is degenerate. @@ -447,10 +457,10 @@ mod tests { 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, }; - let (r, v) = elements.to_cartesian(MU); - let back = Keplerian::from_cartesian(r, v, MU); + let back = elements.to_cartesian().keplerian(); assert_eq!(back.degeneracy, Degeneracy::None); assert_relative_eq!( @@ -478,11 +488,7 @@ mod tests { #[test] fn test_geo_is_reported_as_circular_equatorial() { let radius = 4.2164e7; - let elements = Keplerian::from_cartesian( - Vector3::new(radius, 0.0, 0.0), - Vector3::new(0.0, circular_speed(radius, MU), 0.0), - MU, - ); + let elements = circular_equatorial(radius).keplerian(); assert_eq!(elements.degeneracy, Degeneracy::CircularEquatorial); assert_relative_eq!(elements.semi_major_axis, radius, max_relative = 1e-12); @@ -492,11 +498,12 @@ mod tests { assert_relative_eq!(elements.true_anomaly, 0.0, epsilon = 1e-12); // A quarter turn later the true longitude is 90 degrees. - let quarter = Keplerian::from_cartesian( + let quarter = Cartesian::new( Vector3::new(0.0, radius, 0.0), Vector3::new(-circular_speed(radius, MU), 0.0, 0.0), - MU, - ); + Frame::EARTH_J2000, + ) + .keplerian(); assert_relative_eq!(quarter.true_anomaly, PI / 2.0, max_relative = 1e-12); } @@ -506,11 +513,12 @@ mod tests { let radius = 7.0e6; let speed = circular_speed(radius, MU); // Ascending node along +x, so at +x the argument of latitude is zero. - let elements = Keplerian::from_cartesian( + let elements = Cartesian::new( Vector3::new(radius, 0.0, 0.0), Vector3::new(0.0, speed * 0.5, speed * 3.0_f64.sqrt() / 2.0), - MU, - ); + Frame::EARTH_J2000, + ) + .keplerian(); assert_eq!(elements.degeneracy, Degeneracy::Circular); assert_relative_eq!( @@ -554,6 +562,7 @@ mod tests { raan: 0.0, aop: 0.0, true_anomaly: 0.0, + frame: Frame::EARTH_J2000, degeneracy: Degeneracy::None, }; assert_relative_eq!(elements.periapsis(), 1.4e7, max_relative = 1e-12); @@ -565,21 +574,19 @@ mod tests { #[test] fn test_state_accessors() { - let (r, v) = leo(); - let y = Orbit.with_position(r).with_velocity(v); + let state = leo(); + let y = Orbit.with_cartesian(&state); - let elements = y.keplerian(MU); + let elements = y.keplerian(state.frame); assert_relative_eq!( - y.period(MU), + elements.period(), period(elements.semi_major_axis, MU), max_relative = 1e-12 ); - assert_relative_eq!( - y.specific_energy(MU), - specific_energy(r, v, MU), - max_relative = 1e-12 + assert_eq!( + y.specific_angular_momentum(), + state.position.cross(&state.velocity) ); - assert_eq!(y.specific_angular_momentum(), r.cross(&v)); } #[test] @@ -591,10 +598,11 @@ mod tests { raan: 0.3, aop: 0.6, true_anomaly: 1.2, + frame: Frame::EARTH_J2000, degeneracy: Degeneracy::None, }; - let y = Orbit.with_keplerian(elements, MU); - let back = y.keplerian(MU); + let y = Orbit.with_keplerian(elements); + let back = y.keplerian(Frame::EARTH_J2000); assert_relative_eq!( back.semi_major_axis, elements.semi_major_axis, @@ -616,8 +624,9 @@ mod tests { use anise::prelude::Orbit as AniseOrbit; use hifitime::Epoch; - let (r, v) = leo(); - let ours = Keplerian::from_cartesian(r, v, MU); + let state = leo(); + let (r, v) = (state.position, state.velocity); + let ours = state.keplerian(); let frame = EARTH_J2000.with_mu_km3_s2(MU * 1e-9); let theirs = AniseOrbit::new( @@ -662,7 +671,7 @@ mod tests { max_relative = 1e-9 ); assert_relative_eq!( - period(ours.semi_major_axis, MU), + ours.period(), theirs.period().unwrap().to_seconds(), max_relative = 1e-9 ); @@ -676,7 +685,8 @@ mod tests { use nalgebra::{SVector, Vector6}; use num_dual::DualSVec64; - let (r, v) = leo(); + let leo = leo(); + let (r, v) = (leo.position, leo.velocity); let state = Vector6::new(r[0], r[1], r[2], v[0], v[1], v[2]); // d(sma)/d(state), by automatic differentiation. @@ -684,25 +694,24 @@ mod tests { |y: Vector6>| { let position = y.fixed_rows::<3>(0).into_owned(); let velocity = y.fixed_rows::<3>(3).into_owned(); - Keplerian::from_cartesian(position, velocity, MU).semi_major_axis + Cartesian::new(position, velocity, Frame::EARTH_J2000) + .keplerian() + .semi_major_axis }, &state, ); - assert_relative_eq!( - sma, - Keplerian::from_cartesian(r, v, MU).semi_major_axis, - max_relative = 1e-12 - ); + assert_relative_eq!(sma, leo.keplerian().semi_major_axis, max_relative = 1e-12); // Cross-check the gradient against central differences, which agree only // to about the square root of machine epsilon. let sma_of = |y: Vector6| { - Keplerian::from_cartesian( + Cartesian::new( y.fixed_rows::<3>(0).into_owned(), y.fixed_rows::<3>(3).into_owned(), - MU, + Frame::EARTH_J2000, ) + .keplerian() .semi_major_axis }; for i in 0..6 { diff --git a/src/lib.rs b/src/lib.rs index 78f623b..b44f07b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -9,7 +9,10 @@ #![cfg_attr(feature = "autodiff", doc = include_str!("../readme.md"))] #![warn(missing_docs)] +pub mod body; +pub mod cartesian; pub mod dynamics; +pub mod frames; pub mod kepler; pub mod propagation; pub mod state; @@ -19,7 +22,10 @@ pub mod prelude { pub use crate::dynamics::{Context, Contribution, Dynamics, DynamicsModel}; // The variants, so a script writes `TwoBody { mu }` rather than // `Dynamics::TwoBody { mu }`. Note this puts every model's name in scope. + pub use crate::body::Body; + pub use crate::cartesian::{Cartesian, FromCartesian}; pub use crate::dynamics::Dynamics::*; + pub use crate::frames::{Frame, Orientation}; pub use crate::kepler::{Degeneracy, Keplerian, TranslationalLayout}; pub use crate::propagation::{Accuracy, Propagation}; pub use crate::state::{ diff --git a/src/propagation.rs b/src/propagation.rs index 36d5dd2..ac674e3 100644 --- a/src/propagation.rs +++ b/src/propagation.rs @@ -6,8 +6,10 @@ use ordinary_diffeq::integrator::bs3::BS3; use ordinary_diffeq::prelude::*; use serde::{Deserialize, Serialize}; +use crate::cartesian::Cartesian; use crate::dynamics::{self, Context, Contribution, Dynamics, DynamicsModel}; -use crate::state::{HasPosition, HasVelocity, Layout, State}; +use crate::frames::Frame; +use crate::state::{HasPosition, HasVelocity, Layout, Orbit, State}; /// How precisely the equations are solved. /// @@ -189,6 +191,7 @@ pub struct Propagation< > { epoch: Epoch, initial_state: State, + frame: Frame, duration: Option, specs: Vec, integrator: S, @@ -210,6 +213,7 @@ where Self { epoch, initial_state, + frame: Frame::default(), duration: None, specs: dynamics.to_vec(), integrator: Vern7::new().a_tol(1e-10).r_tol(1e-10), @@ -219,12 +223,40 @@ where } } +impl Propagation> +where + T: Scalar + 'static, + Vern7<6>: Integrator<6>, +{ + /// A propagation of a translational state, in that state's own frame. + /// + /// The shorthand for the common case: a [`Cartesian`] already carries its + /// frame, so neither the layout nor the frame has to be named. + /// + /// Use [`Propagation::new`] with a wider layout when the state also carries + /// mass, attitude or estimated parameters. + pub fn from_cartesian(epoch: Epoch, state: Cartesian, dynamics: &[Dynamics]) -> Self { + let frame = state.frame; + Self::new(epoch, state.to_state(Orbit), dynamics).frame(frame) + } +} + impl Propagation where T: Scalar + 'static, L: HasPosition + HasVelocity + 'static, S: Tolerances + Copy, { + /// The frame the state is expressed in. + /// + /// Defaults to [`Frame::EARTH_J2000`]. The default lives here rather than in + /// the force models, so a propagation about another body is a matter of saying + /// so rather than of overriding an assumption in several places. + pub fn frame(mut self, frame: Frame) -> Self { + self.frame = frame; + self + } + /// How long to propagate for. pub fn duration(mut self, duration: Duration) -> Self { self.duration = Some(duration); @@ -282,6 +314,7 @@ where Propagation { epoch: self.epoch, initial_state: self.initial_state, + frame: self.frame, duration: self.duration, specs: self.specs, integrator, @@ -312,7 +345,7 @@ where let assembly = Assembly { models: dynamics::build::(&self.specs), - context: Context::new(self.epoch), + context: Context::new(self.epoch, self.frame), }; let ode = ODE::new( @@ -331,6 +364,7 @@ where #[cfg(test)] mod tests { use super::*; + use crate::body::Body; use crate::state::{Full, Orbit, OrbitMass}; use approx::assert_relative_eq; use nalgebra::{Vector3, Vector6}; @@ -344,7 +378,7 @@ mod tests { } fn two_body() -> Vec { - vec![Dynamics::TwoBody { mu: MU }] + vec![Dynamics::TwoBody { body: Body::EARTH }] } fn leo() -> State { @@ -565,7 +599,7 @@ mod tests { vec![Box::new(ConstantForce { force: Vector3::new(10.0, 0.0, 0.0), })], - Context::new(epoch()), + Context::new(epoch(), Frame::EARTH_J2000), ); let y = OrbitMass @@ -587,9 +621,9 @@ mod tests { Box::new(ConstantForce { force: Vector3::new(10.0, 0.0, 0.0), }), - Dynamics::TwoBody { mu: MU }.build(), + Dynamics::TwoBody { body: Body::EARTH }.build(), ], - Context::new(epoch()), + Context::new(epoch(), Frame::EARTH_J2000), ); let y = OrbitMass @@ -612,7 +646,7 @@ mod tests { vec![Box::new(ConstantForce { force: Vector3::new(10.0, 0.0, 0.0), })], - Context::new(epoch()), + Context::new(epoch(), Frame::EARTH_J2000), ); let _ = state_derivative(0.0, leo().vector, &assembly); }