diff --git a/Cargo.toml b/Cargo.toml index 37fe182b..12e54cf5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -77,12 +77,13 @@ typed-builder = "0.20.0" pythonize = { version = "0.21", optional = true } snafu = { version = "0.8.3", features = ["backtrace"] } serde_dhall = "0.12" -toml = "0.8.14" + [dev-dependencies] polars = { version = "0.43.1", features = ["parquet"] } rstest = "0.22.0" pretty_env_logger = "0.5" +toml = "0.8.14" [build-dependencies] shadow-rs = "0.35.0" diff --git a/examples/03_geo_analysis/raise.rs b/examples/03_geo_analysis/raise.rs index 649caf61..d1ad84f4 100644 --- a/examples/03_geo_analysis/raise.rs +++ b/examples/03_geo_analysis/raise.rs @@ -19,7 +19,7 @@ use nyx::{ }, io::{gravity::HarmonicsMem, ExportCfg}, md::{prelude::Objective, StateParameter}, - propagators::{PropOpts, Propagator, RSSCartesianStep}, + propagators::{ErrorControl, IntegratorOptions, Propagator}, Spacecraft, }; use std::{error::Error, sync::Arc}; @@ -117,9 +117,9 @@ fn main() -> Result<(), Box> { // We specify a minimum step in the propagator because the Ruggiero control would otherwise drive this step very low. let (final_state, traj) = Propagator::rk89( sc_dynamics.clone(), - PropOpts::builder() + IntegratorOptions::builder() .min_step(10.0_f64.seconds()) - .error_ctrl(RSSCartesianStep {}) + .error_ctrl(ErrorControl::RSSCartesianStep) .build(), ) .with(sc, almanac.clone()) diff --git a/examples/03_geo_analysis/stationkeeping.rs b/examples/03_geo_analysis/stationkeeping.rs index afc6358c..eef6a45d 100644 --- a/examples/03_geo_analysis/stationkeeping.rs +++ b/examples/03_geo_analysis/stationkeeping.rs @@ -20,7 +20,7 @@ use nyx::{ io::{gravity::HarmonicsMem, ExportCfg}, mc::{MonteCarlo, MultivariateNormal, StateDispersion}, md::{prelude::Objective, StateParameter}, - propagators::{PropOpts, Propagator, RSSCartesianStep}, + propagators::{ErrorControl, IntegratorOptions, Propagator}, Spacecraft, State, }; use std::{error::Error, sync::Arc}; @@ -103,9 +103,9 @@ fn main() -> Result<(), Box> { // Build the propagator setup. let setup = Propagator::rk89( sc_dynamics.clone(), - PropOpts::builder() + IntegratorOptions::builder() .min_step(10.0_f64.seconds()) - .error_ctrl(RSSCartesianStep {}) + .error_ctrl(ErrorControl::RSSCartesianStep) .build(), ); diff --git a/src/cosmic/mod.rs b/src/cosmic/mod.rs index da23739f..4352e9b7 100644 --- a/src/cosmic/mod.rs +++ b/src/cosmic/mod.rs @@ -110,6 +110,12 @@ where fn set_value(&mut self, param: StateParameter, _val: f64) -> Result<(), StateError> { Err(StateError::Unavailable { param }) } + + /// Returns a copy of the orbit + fn orbit(&self) -> Orbit; + + /// Modifies this state's orbit + fn set_orbit(&mut self, _orbit: Orbit) {} } pub fn assert_orbit_eq_or_abs(left: &Orbit, right: &Orbit, epsilon: f64, msg: &str) { diff --git a/src/cosmic/spacecraft.rs b/src/cosmic/spacecraft.rs index eb7d0204..fcc4d16e 100644 --- a/src/cosmic/spacecraft.rs +++ b/src/cosmic/spacecraft.rs @@ -761,6 +761,14 @@ impl State for Spacecraft { fn unset_stm(&mut self) { self.stm = None; } + + fn orbit(&self) -> Orbit { + self.orbit + } + + fn set_orbit(&mut self, orbit: Orbit) { + self.orbit = orbit; + } } impl Add>> for Spacecraft { diff --git a/src/dynamics/mod.rs b/src/dynamics/mod.rs index 73f44e8a..a356a5ee 100644 --- a/src/dynamics/mod.rs +++ b/src/dynamics/mod.rs @@ -172,6 +172,7 @@ pub trait AccelModel: Send + Sync + fmt::Display { /// Stores dynamical model errors #[derive(Debug, PartialEq, Snafu)] +#[snafu(visibility(pub(crate)))] pub enum DynamicsError { /// Fuel exhausted at the provided spacecraft state #[snafu(display("fuel exhausted at {sc}"))] diff --git a/src/io/mod.rs b/src/io/mod.rs index bbb32ff9..0031377e 100644 --- a/src/io/mod.rs +++ b/src/io/mod.rs @@ -159,13 +159,13 @@ impl ExportCfg { #[derive(Debug, Snafu)] #[snafu(visibility(pub(crate)))] pub enum ConfigError { - #[snafu(display("Failed to read configuration file: {source}"))] + #[snafu(display("failed to read configuration file: {source}"))] ReadError { source: io::Error }, - #[snafu(display("Failed to parse YAML configuration file: {source}"))] + #[snafu(display("failed to parse YAML configuration file: {source}"))] ParseError { source: serde_yaml::Error }, - #[snafu(display("Invalid configuration: {msg}"))] + #[snafu(display("of invalid configuration: {msg}"))] InvalidConfig { msg: String }, } diff --git a/src/mc/montecarlo.rs b/src/mc/montecarlo.rs index 4bdb73f4..9ea01f24 100644 --- a/src/mc/montecarlo.rs +++ b/src/mc/montecarlo.rs @@ -24,7 +24,7 @@ use crate::mc::results::{PropResult, Results, Run}; use crate::mc::DispersedState; use crate::md::trajectory::Interpolatable; use crate::md::EventEvaluator; -use crate::propagators::{ErrorCtrl, Propagator}; +use crate::propagators::Propagator; #[cfg(not(target_arch = "wasm32"))] use crate::time::Unit; use crate::time::{Duration, Epoch}; @@ -89,9 +89,9 @@ where /// Generate states and propagate each independently until a specific event is found `trigger` times. #[allow(clippy::needless_lifetimes)] - pub fn run_until_nth_event<'a, D, E, F>( + pub fn run_until_nth_event( self, - prop: Propagator<'a, D, E>, + prop: Propagator, almanac: Arc, max_duration: Duration, event: &F, @@ -100,7 +100,7 @@ where ) -> Results> where D: Dynamics, - E: ErrorCtrl, + F: EventEvaluator, DefaultAllocator: Allocator<::Size> + Allocator<::Size, ::Size> @@ -113,9 +113,9 @@ where /// Generate states and propagate each independently until a specific event is found `trigger` times. #[must_use = "Monte Carlo result must be used"] #[allow(clippy::needless_lifetimes)] - pub fn resume_run_until_nth_event<'a, D, E, F>( + pub fn resume_run_until_nth_event( &self, - prop: Propagator<'a, D, E>, + prop: Propagator, almanac: Arc, skip: usize, max_duration: Duration, @@ -125,7 +125,7 @@ where ) -> Results> where D: Dynamics, - E: ErrorCtrl, + F: EventEvaluator, DefaultAllocator: Allocator<::Size> + Allocator<::Size, ::Size> @@ -188,16 +188,16 @@ where /// Generate states and propagate each independently until a specific event is found `trigger` times. #[must_use = "Monte Carlo result must be used"] #[allow(clippy::needless_lifetimes)] - pub fn run_until_epoch<'a, D, E>( + pub fn run_until_epoch( self, - prop: Propagator<'a, D, E>, + prop: Propagator, almanac: Arc, end_epoch: Epoch, num_runs: usize, ) -> Results> where D: Dynamics, - E: ErrorCtrl, + DefaultAllocator: Allocator<::Size> + Allocator<::Size, ::Size> + Allocator<::VecLength>, @@ -209,9 +209,9 @@ where /// Resumes a Monte Carlo run by skipping the first `skip` items, generating states only after that, and propagate each independently until the specified epoch. #[must_use = "Monte Carlo result must be used"] #[allow(clippy::needless_lifetimes)] - pub fn resume_run_until_epoch<'a, D, E>( + pub fn resume_run_until_epoch( &self, - prop: Propagator<'a, D, E>, + prop: Propagator, almanac: Arc, skip: usize, end_epoch: Epoch, @@ -219,7 +219,7 @@ where ) -> Results> where D: Dynamics, - E: ErrorCtrl, + DefaultAllocator: Allocator<::Size> + Allocator<::Size, ::Size> + Allocator<::VecLength>, diff --git a/src/md/mod.rs b/src/md/mod.rs index 05db997a..7045d5b7 100644 --- a/src/md/mod.rs +++ b/src/md/mod.rs @@ -25,7 +25,7 @@ use snafu::prelude::*; pub mod prelude { pub use super::{ - optimizer::*, + targeter::*, trajectory::{ExportCfg, Interpolatable, Traj}, Event, ScTraj, StateParameter, }; @@ -36,7 +36,7 @@ pub mod prelude { pub use crate::dynamics::{Dynamics, NyxError}; pub use crate::io::gravity::HarmonicsMem; pub use crate::md::objective::Objective; - pub use crate::propagators::{PropOpts, Propagator}; + pub use crate::propagators::{IntegratorOptions, Propagator}; pub use crate::time::{Duration, Epoch, TimeUnits, Unit}; pub use crate::Spacecraft; pub use crate::{State, TimeTagged}; @@ -52,7 +52,7 @@ pub use events::{Event, EventEvaluator}; pub mod objective; pub mod opti; -pub use opti::optimizer; +pub use opti::targeter; pub type ScTraj = trajectory::Traj; // pub type Ephemeris = trajectory::Traj; diff --git a/src/md/opti/mod.rs b/src/md/opti/mod.rs index 7f09afa0..1524a795 100644 --- a/src/md/opti/mod.rs +++ b/src/md/opti/mod.rs @@ -19,16 +19,16 @@ // pub mod convert_impulsive; pub mod multipleshooting; pub use multipleshooting::{ctrlnodes, multishoot}; -/// Uses a Levenberg Marquardt minimizer to solve the damped least squares problem. -// #[cfg(feature = "broken-donotuse")] -// pub mod minimize_lm; -pub mod optimizer; /// Uses a [Newton Raphson](https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization) method where the Jacobian is computed via finite differencing. pub mod raphson_finite_diff; /// Uses a [Newton Raphson](https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization) method where the Jacobian is computed via hyperdual numbers. pub mod raphson_hyperdual; pub mod solution; pub mod target_variable; +/// Uses a Levenberg Marquardt minimizer to solve the damped least squares problem. +// #[cfg(feature = "broken-donotuse")] +// pub mod minimize_lm; +pub mod targeter; #[derive(Copy, Clone, Debug, PartialEq, Eq)] pub enum DiffMethod { diff --git a/src/md/opti/multipleshooting/altitude_heuristic.rs b/src/md/opti/multipleshooting/altitude_heuristic.rs index 7f29842b..f61c419c 100644 --- a/src/md/opti/multipleshooting/altitude_heuristic.rs +++ b/src/md/opti/multipleshooting/altitude_heuristic.rs @@ -27,10 +27,9 @@ use super::{ }; use crate::errors::TargetingError; use crate::md::{prelude::*, PropSnafu}; -use crate::propagators::error_ctrl::ErrorCtrl; use crate::{Orbit, Spacecraft}; -impl<'a, E: ErrorCtrl> MultipleShooting<'a, E, Node, 3, 3> { +impl<'a> MultipleShooting<'a, Node, 3, 3> { /// Builds a multiple shooting structure assuming that the optimal trajectory is near a linear /// heuristic in geodetic altitude and direction. /// For example, if x0 has an altitude of 100 km and xf has an altitude @@ -42,7 +41,7 @@ impl<'a, E: ErrorCtrl> MultipleShooting<'a, E, Node, 3, 3> { node_count: usize, angular_velocity_deg_s: f64, body_frame: Frame, - prop: &'a Propagator<'a, SpacecraftDynamics, E>, + prop: &'a Propagator, almanac: Arc, ) -> Result { if node_count < 3 { diff --git a/src/md/opti/multipleshooting/equidistant_heuristic.rs b/src/md/opti/multipleshooting/equidistant_heuristic.rs index 976d215f..ffa6f729 100644 --- a/src/md/opti/multipleshooting/equidistant_heuristic.rs +++ b/src/md/opti/multipleshooting/equidistant_heuristic.rs @@ -21,10 +21,9 @@ use super::multishoot::MultipleShooting; pub use super::CostFunction; use crate::errors::TargetingError; use crate::md::prelude::*; -use crate::propagators::error_ctrl::ErrorCtrl; use crate::{Orbit, Spacecraft}; -impl<'a, E: ErrorCtrl> MultipleShooting<'a, E, Node, 3, 3> { +impl<'a> MultipleShooting<'a, Node, 3, 3> { /// Builds a multiple shooting structure assuming that the optimal trajectory is a straight line /// between the start and end points. The position of the nodes will be update at each iteration /// of the outer loop. @@ -33,7 +32,7 @@ impl<'a, E: ErrorCtrl> MultipleShooting<'a, E, Node, 3, 3> { x0: Spacecraft, xf: Orbit, node_count: usize, - prop: &'a Propagator<'a, SpacecraftDynamics, E>, + prop: &'a Propagator, ) -> Result { if node_count < 3 { error!("At least three nodes are needed for a multiple shooting optimization"); diff --git a/src/md/opti/multipleshooting/multishoot.rs b/src/md/opti/multipleshooting/multishoot.rs index 25669763..c31ac95b 100644 --- a/src/md/opti/multipleshooting/multishoot.rs +++ b/src/md/opti/multipleshooting/multishoot.rs @@ -22,9 +22,8 @@ pub use super::CostFunction; use super::{MultipleShootingError, TargetingSnafu}; use crate::linalg::{DMatrix, DVector, SVector}; use crate::md::opti::solution::TargeterSolution; -use crate::md::optimizer::Optimizer; +use crate::md::targeter::Targeter; use crate::md::{prelude::*, TargetingError}; -use crate::propagators::error_ctrl::ErrorCtrl; use crate::pseudo_inverse; use crate::{Orbit, Spacecraft}; @@ -39,15 +38,9 @@ pub trait MultishootNode: Copy + Into<[Objective; O]> { /// Source of implementation: "Low Thrust Optimization in Cislunar and Translunar space", 2018 Nathan Re (Parrish) /// OT: size of the objectives for each node (e.g. 3 if the objectives are X, Y, Z). /// VT: size of the variables for targeter node (e.g. 4 if the objectives are thrust direction (x,y,z) and thrust level). -pub struct MultipleShooting< - 'a, - E: ErrorCtrl, - T: MultishootNode, - const VT: usize, - const OT: usize, -> { +pub struct MultipleShooting<'a, T: MultishootNode, const VT: usize, const OT: usize> { /// The propagator setup (kind, stages, etc.) - pub prop: &'a Propagator<'a, SpacecraftDynamics, E>, + pub prop: &'a Propagator, /// List of nodes of the optimal trajectory pub targets: Vec, /// Starting point, must be a spacecraft equipped with a thruster @@ -66,9 +59,7 @@ pub struct MultipleShooting< pub all_dvs: Vec>, } -impl<'a, E: ErrorCtrl, T: MultishootNode, const VT: usize, const OT: usize> - MultipleShooting<'a, E, T, VT, OT> -{ +impl<'a, T: MultishootNode, const VT: usize, const OT: usize> MultipleShooting<'a, T, VT, OT> { /// Solve the multiple shooting problem by finding the arrangement of nodes to minimize the cost function. pub fn solve( &mut self, @@ -90,7 +81,7 @@ impl<'a, E: ErrorCtrl, T: MultishootNode, const VT: usize, const OT: usize> /* *** ** 1. Solve the delta-v differential corrector between each node ** *** */ - let tgt = Optimizer { + let tgt = Targeter { prop: self.prop, objectives: self.targets[i].into(), variables: self.variables, @@ -134,7 +125,7 @@ impl<'a, E: ErrorCtrl, T: MultishootNode, const VT: usize, const OT: usize> ** Note that because the first initial_state is x0, the i-th "initial state" ** is the initial state to reach the i-th node. ** *** */ - let inner_tgt_a = Optimizer::delta_v(self.prop, next_node); + let inner_tgt_a = Targeter::delta_v(self.prop, next_node); let inner_sol_a = inner_tgt_a .try_achieve_dual( initial_states[i], @@ -160,7 +151,7 @@ impl<'a, E: ErrorCtrl, T: MultishootNode, const VT: usize, const OT: usize> /* *** ** 2.C. Rerun the targeter from the new state at the perturbed node to the next unpertubed node ** *** */ - let inner_tgt_b = Optimizer::delta_v(self.prop, self.targets[i + 1].into()); + let inner_tgt_b = Targeter::delta_v(self.prop, self.targets[i + 1].into()); let inner_sol_b = inner_tgt_b .try_achieve_dual( inner_sol_a.achieved_state, @@ -250,7 +241,7 @@ impl<'a, E: ErrorCtrl, T: MultishootNode, const VT: usize, const OT: usize> for (i, node) in self.targets.iter().enumerate() { // Run the unpertubed targeter - let tgt = Optimizer::delta_v(self.prop, (*node).into()); + let tgt = Targeter::delta_v(self.prop, (*node).into()); let sol = tgt .try_achieve_dual( initial_states[i], @@ -287,8 +278,8 @@ impl<'a, E: ErrorCtrl, T: MultishootNode, const VT: usize, const OT: usize> } } -impl<'a, E: ErrorCtrl, T: MultishootNode, const VT: usize, const OT: usize> fmt::Display - for MultipleShooting<'a, E, T, VT, OT> +impl<'a, T: MultishootNode, const VT: usize, const OT: usize> fmt::Display + for MultipleShooting<'a, T, VT, OT> { #[allow(clippy::or_fun_call, clippy::clone_on_copy)] fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { @@ -357,15 +348,15 @@ impl, const O: usize> fmt::Display for MultipleShootingSolu impl, const O: usize> MultipleShootingSolution { /// Allows building the trajectories between different nodes /// This will rebuild the targeters and apply the solutions sequentially - pub fn build_trajectories<'a, E: ErrorCtrl>( + pub fn build_trajectories( &self, - prop: &'a Propagator<'a, SpacecraftDynamics, E>, + prop: &Propagator, almanac: Arc, ) -> Result, MultipleShootingError> { let mut trajz = Vec::with_capacity(self.nodes.len()); - for (i, node) in self.nodes.iter().enumerate() { - let (_, traj) = Optimizer::delta_v(prop, (*node).into()) + for (i, node) in self.nodes.iter().copied().enumerate() { + let (_, traj) = Targeter::delta_v(prop, node.into()) .apply_with_traj(&self.solutions[i], almanac.clone()) .context(TargetingSnafu { segment: i })?; trajz.push(traj); diff --git a/src/md/opti/raphson_finite_diff.rs b/src/md/opti/raphson_finite_diff.rs index 88b9112a..260fd045 100644 --- a/src/md/opti/raphson_finite_diff.rs +++ b/src/md/opti/raphson_finite_diff.rs @@ -16,8 +16,8 @@ along with this program. If not, see . */ -use super::optimizer::Optimizer; use super::solution::TargeterSolution; +use super::targeter::Targeter; use crate::cosmic::{AstroAlmanacSnafu, AstroPhysicsSnafu}; use crate::dynamics::guidance::{GuidanceError, LocalFrame, Mnvr}; use crate::errors::TargetingError; @@ -26,7 +26,6 @@ use crate::md::{prelude::*, AstroSnafu, GuidanceSnafu, UnderdeterminedProblemSna use crate::md::{PropSnafu, StateParameter}; pub use crate::md::{Variable, Vary}; use crate::polyfit::CommonPolynomial; -use crate::propagators::error_ctrl::ErrorCtrl; use crate::pseudo_inverse; use hifitime::TimeUnits; use rayon::prelude::*; @@ -34,7 +33,7 @@ use snafu::{ensure, ResultExt}; #[cfg(not(target_arch = "wasm32"))] use std::time::Instant; -impl<'a, E: ErrorCtrl, const V: usize, const O: usize> Optimizer<'a, E, V, O> { +impl<'a, const V: usize, const O: usize> Targeter<'a, V, O> { /// Differential correction using finite differencing #[allow(clippy::comparison_chain)] pub fn try_achieve_fd( diff --git a/src/md/opti/raphson_hyperdual.rs b/src/md/opti/raphson_hyperdual.rs index 35da9725..eda4ca32 100644 --- a/src/md/opti/raphson_hyperdual.rs +++ b/src/md/opti/raphson_hyperdual.rs @@ -25,13 +25,12 @@ use crate::linalg::{DMatrix, SVector}; use crate::md::{prelude::*, PropSnafu, UnderdeterminedProblemSnafu}; use crate::md::{AstroSnafu, StateParameter}; pub use crate::md::{Variable, Vary}; -use crate::propagators::error_ctrl::ErrorCtrl; use crate::pseudo_inverse; use crate::utils::are_eigenvalues_stable; #[cfg(not(target_arch = "wasm32"))] use std::time::Instant; -impl<'a, E: ErrorCtrl, const V: usize, const O: usize> Optimizer<'a, E, V, O> { +impl<'a, const V: usize, const O: usize> Targeter<'a, V, O> { /// Differential correction using hyperdual numbers for the objectives #[allow(clippy::comparison_chain)] pub fn try_achieve_dual( diff --git a/src/md/opti/optimizer.rs b/src/md/opti/targeter.rs similarity index 89% rename from src/md/opti/optimizer.rs rename to src/md/opti/targeter.rs index 3ef03dd7..389e4edf 100644 --- a/src/md/opti/optimizer.rs +++ b/src/md/opti/targeter.rs @@ -26,16 +26,15 @@ use crate::md::AstroSnafu; use crate::md::PropSnafu; use crate::md::StateParameter; pub use crate::md::{Variable, Vary}; -use crate::propagators::error_ctrl::ErrorCtrl; use std::fmt; use super::solution::TargeterSolution; /// An optimizer structure with V control variables and O objectives. #[derive(Clone)] -pub struct Optimizer<'a, E: ErrorCtrl, const V: usize, const O: usize> { +pub struct Targeter<'a, const V: usize, const O: usize> { /// The propagator setup (kind, stages, etc.) - pub prop: &'a Propagator<'a, SpacecraftDynamics, E>, + pub prop: &'a Propagator, /// The list of objectives of this targeter pub objectives: [Objective; O], /// An optional frame (and Cosm) to compute the objectives in. @@ -49,7 +48,7 @@ pub struct Optimizer<'a, E: ErrorCtrl, const V: usize, const O: usize> { pub iterations: usize, } -impl<'a, E: ErrorCtrl, const V: usize, const O: usize> fmt::Display for Optimizer<'a, E, V, O> { +impl<'a, const V: usize, const O: usize> fmt::Display for Targeter<'a, V, O> { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { let mut objmsg = String::from(""); for obj in &self.objectives { @@ -65,12 +64,9 @@ impl<'a, E: ErrorCtrl, const V: usize, const O: usize> fmt::Display for Optimize } } -impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 3, O> { +impl<'a, const O: usize> Targeter<'a, 3, O> { /// Create a new Targeter which will apply an impulsive delta-v correction. - pub fn delta_v( - prop: &'a Propagator<'a, SpacecraftDynamics, E>, - objectives: [Objective; O], - ) -> Self { + pub fn delta_v(prop: &'a Propagator, objectives: [Objective; O]) -> Self { Self { prop, objectives, @@ -86,10 +82,7 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 3, O> { } /// Create a new Targeter which will MOVE the position of the spacecraft at the correction epoch - pub fn delta_r( - prop: &'a Propagator<'a, SpacecraftDynamics, E>, - objectives: [Objective; O], - ) -> Self { + pub fn delta_r(prop: &'a Propagator, objectives: [Objective; O]) -> Self { Self { prop, objectives, @@ -105,10 +98,7 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 3, O> { } /// Create a new Targeter which will apply an impulsive delta-v correction on all components of the VNC frame. By default, max step is 0.5 km/s. - pub fn vnc( - prop: &'a Propagator<'a, SpacecraftDynamics, E>, - objectives: [Objective; O], - ) -> Self { + pub fn vnc(prop: &'a Propagator, objectives: [Objective; O]) -> Self { Self { prop, objectives, @@ -124,10 +114,10 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 3, O> { } } -impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 4, O> { +impl<'a, const O: usize> Targeter<'a, 4, O> { /// Create a new Targeter which will apply a continuous thrust for the whole duration of the segment pub fn thrust_dir( - prop: &'a Propagator<'a, SpacecraftDynamics, E>, + prop: &'a Propagator, objectives: [Objective; O], ) -> Self { Self { @@ -146,10 +136,10 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 4, O> { } } -impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 7, O> { +impl<'a, const O: usize> Targeter<'a, 7, O> { /// Create a new Targeter which will apply a continuous thrust for the whole duration of the segment pub fn thrust_dir_rate( - prop: &'a Propagator<'a, SpacecraftDynamics, E>, + prop: &'a Propagator, objectives: [Objective; O], ) -> Self { Self { @@ -171,10 +161,10 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 7, O> { } } -impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 10, O> { +impl<'a, const O: usize> Targeter<'a, 10, O> { /// Create a new Targeter which will apply a continuous thrust for the whole duration of the segment pub fn thrust_profile( - prop: &'a Propagator<'a, SpacecraftDynamics, E>, + prop: &'a Propagator, objectives: [Objective; O], ) -> Self { Self { @@ -199,10 +189,10 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 10, O> { } } -impl<'a, E: ErrorCtrl, const V: usize, const O: usize> Optimizer<'a, E, V, O> { +impl<'a, const V: usize, const O: usize> Targeter<'a, V, O> { /// Create a new Targeter which will apply an impulsive delta-v correction. pub fn new( - prop: &'a Propagator<'a, SpacecraftDynamics, E>, + prop: &'a Propagator, variables: [Variable; V], objectives: [Objective; O], ) -> Self { @@ -218,7 +208,7 @@ impl<'a, E: ErrorCtrl, const V: usize, const O: usize> Optimizer<'a, E, V, O> { /// Create a new Targeter which will apply an impulsive delta-v correction. pub fn in_frame( - prop: &'a Propagator<'a, SpacecraftDynamics, E>, + prop: &'a Propagator, variables: [Variable; V], objectives: [Objective; O], objective_frame: Frame, @@ -235,7 +225,7 @@ impl<'a, E: ErrorCtrl, const V: usize, const O: usize> Optimizer<'a, E, V, O> { /// Create a new Targeter which will apply an impulsive delta-v correction on the specified components of the VNC frame. pub fn vnc_with_components( - prop: &'a Propagator<'a, SpacecraftDynamics, E>, + prop: &'a Propagator, variables: [Variable; V], objectives: [Objective; O], ) -> Self { diff --git a/src/md/trajectory/interpolatable.rs b/src/md/trajectory/interpolatable.rs index 32e54641..b845cd5b 100644 --- a/src/md/trajectory/interpolatable.rs +++ b/src/md/trajectory/interpolatable.rs @@ -47,9 +47,6 @@ where /// List of state parameters that will be exported to a trajectory file in addition to the epoch (provided in this different formats). fn export_params() -> Vec; - - /// Returns the orbit - fn orbit(&self) -> &Orbit; } impl Interpolatable for Spacecraft { @@ -156,8 +153,4 @@ impl Interpolatable for Spacecraft { ] .concat() } - - fn orbit(&self) -> &Orbit { - &self.orbit - } } diff --git a/src/md/trajectory/traj.rs b/src/md/trajectory/traj.rs index 58ba2f45..01a2c302 100644 --- a/src/md/trajectory/traj.rs +++ b/src/md/trajectory/traj.rs @@ -445,8 +445,8 @@ where // Build an array of all the RIC differences let mut ric_diff = Vec::with_capacity(other_states.len()); for (ii, other_state) in other_states.iter().enumerate() { - let self_orbit = *self_states[ii].orbit(); - let other_orbit = *other_state.orbit(); + let self_orbit = self_states[ii].orbit(); + let other_orbit = other_state.orbit(); let this_ric_diff = self_orbit.ric_difference(&other_orbit).map_err(Box::new)?; diff --git a/src/od/mod.rs b/src/od/mod.rs index 821074de..c017b587 100644 --- a/src/od/mod.rs +++ b/src/od/mod.rs @@ -64,7 +64,6 @@ pub mod snc; pub type SpacecraftODProcess<'a> = self::process::ODProcess< 'a, crate::md::prelude::SpacecraftDynamics, - crate::propagators::RSSCartesianStep, msr::RangeDoppler, nalgebra::Const<3>, crate::Spacecraft, diff --git a/src/od/process/export.rs b/src/od/process/export.rs index ecdf73e1..d2e778dc 100644 --- a/src/od/process/export.rs +++ b/src/od/process/export.rs @@ -23,7 +23,6 @@ use crate::linalg::{DefaultAllocator, DimName}; use crate::md::trajectory::Interpolatable; use crate::md::StateParameter; use crate::od::estimate::*; -use crate::propagators::error_ctrl::ErrorCtrl; use crate::State; use crate::{od::*, Spacecraft}; use arrow::array::{Array, BooleanBuilder, Float64Builder, StringBuilder}; @@ -41,8 +40,8 @@ use std::path::{Path, PathBuf}; use super::ODProcess; -impl<'a, D: Dynamics, E: ErrorCtrl, Msr: Measurement, A: DimName> - ODProcess<'a, D, E, Msr, A, Spacecraft, KF> +impl<'a, D: Dynamics, Msr: Measurement, A: DimName> + ODProcess<'a, D, Msr, A, Spacecraft, KF> where D::StateType: Interpolatable + Add::Size>, Output = D::StateType>, diff --git a/src/od/process/mod.rs b/src/od/process/mod.rs index d4a44782..6d75c4a9 100644 --- a/src/od/process/mod.rs +++ b/src/od/process/mod.rs @@ -23,7 +23,6 @@ pub use crate::od::estimate::*; pub use crate::od::ground_station::*; pub use crate::od::snc::*; pub use crate::od::*; -use crate::propagators::error_ctrl::ErrorCtrl; use crate::propagators::PropInstance; pub use crate::time::{Duration, Unit}; use anise::prelude::Almanac; @@ -45,7 +44,6 @@ mod export; pub struct ODProcess< 'a, D: Dynamics, - E: ErrorCtrl, Msr: Measurement, A: DimName, S: EstimateFrom + Interpolatable, @@ -73,7 +71,7 @@ pub struct ODProcess< + Allocator::Size>, { /// PropInstance used for the estimation - pub prop: PropInstance<'a, D, E>, + pub prop: PropInstance<'a, D>, /// Kalman filter itself pub kf: K, /// Vector of estimates available after a pass @@ -91,12 +89,11 @@ pub struct ODProcess< impl< 'a, D: Dynamics, - E: ErrorCtrl, Msr: Measurement, A: DimName, S: EstimateFrom + Interpolatable, K: Filter, - > ODProcess<'a, D, E, Msr, A, S, K> + > ODProcess<'a, D, Msr, A, S, K> where D::StateType: Interpolatable + Add::Size>, Output = D::StateType>, ::VecLength>>::Buffer: Send, @@ -123,7 +120,7 @@ where { /// Initialize a new orbit determination process with an optional trigger to switch from a CKF to an EKF. pub fn new( - prop: PropInstance<'a, D, E>, + prop: PropInstance<'a, D>, kf: K, ekf_trigger: Option, resid_crit: Option, @@ -145,7 +142,7 @@ where /// Initialize a new orbit determination process with an Extended Kalman filter. The switch from a classical KF to an EKF is based on the provided trigger. pub fn ekf( - prop: PropInstance<'a, D, E>, + prop: PropInstance<'a, D>, kf: K, trigger: EkfTrigger, resid_crit: Option, @@ -286,7 +283,7 @@ where where Dev: TrackingDeviceSim, { - // TODO(now): Add ExportCfg to iterate and to process so the data can be exported as we process it. Consider a thread writing with channel for faster serialization. + // TODO: Add ExportCfg to iterate and to process so the data can be exported as we process it. Consider a thread writing with channel for faster serialization. let mut best_rms = self.rms_residual_ratios(); let mut previous_rms = best_rms; @@ -754,12 +751,11 @@ where impl< 'a, D: Dynamics, - E: ErrorCtrl, Msr: Measurement, A: DimName, S: EstimateFrom + Interpolatable, K: Filter, - > ODProcess<'a, D, E, Msr, A, S, K> + > ODProcess<'a, D, Msr, A, S, K> where D::StateType: Interpolatable + Add::Size>, Output = D::StateType>, ::VecLength>>::Buffer: Send, @@ -785,7 +781,7 @@ where + Allocator::Size>, { pub fn ckf( - prop: PropInstance<'a, D, E>, + prop: PropInstance<'a, D>, kf: K, resid_crit: Option, almanac: Arc, diff --git a/src/propagators/error_ctrl.rs b/src/propagators/error_ctrl.rs index a38e21ca..d3c0b232 100644 --- a/src/propagators/error_ctrl.rs +++ b/src/propagators/error_ctrl.rs @@ -16,43 +16,67 @@ along with this program. If not, see . */ +use serde::{Deserialize, Serialize}; + use crate::linalg::allocator::Allocator; -use crate::linalg::{DefaultAllocator, DimName, OVector, U1, U3}; +use crate::linalg::{DefaultAllocator, DimName, OVector, U3}; // This determines when to take into consideration the magnitude of the state_delta and // prevents dividing by too small of a number. const REL_ERR_THRESH: f64 = 0.1; -/// The Error Control trait manages how a propagator computes the error in the current step. -pub trait ErrorCtrl -where - Self: Copy + Send + Sync, -{ +/// The Error Control manages how a propagator computes the error in the current step. +#[derive(Copy, Clone, Debug, PartialEq, Eq, Serialize, Deserialize)] +pub enum ErrorControl { + /// An RSS state error control which effectively for the provided vector composed of two vectors of the same unit, both of size 3 (e.g. position + velocity). + RSSCartesianState, + /// An RSS step error control which effectively for the provided vector composed of two vectors of the same unit, both of size 3 (e.g. position + velocity). + RSSCartesianStep, + /// An RSS state error control: when in doubt, use this error controller, especially for high accurracy. + /// + /// Here is the warning from GMAT R2016a on this error controller: + /// > This is a more stringent error control method than [`rss_step`] that is often used as the default in other software such as STK. + /// > If you set [the] accuracy to a very small number, 1e-13 for example, and set the error control to [`rss_step`], integrator + /// > performance will be poor, for little if any improvement in the accuracy of the orbit integration. + /// > For more best practices of these integrators (which clone those in GMAT), please refer to the + /// > [GMAT reference](https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/doc/help/src/Resource_NumericalIntegrators.xml#L1292). + /// > (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3004] + RSSState, + /// An RSS step error control which effectively computes the L2 norm of the provided Vector of size 3 + /// + /// Note that this error controller should be preferably be used only with slices of a state with the same units. + /// For example, one should probably use this for position independently of using it for the velocity. + /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3045] + RSSStep, + /// A largest error control which effectively computes the largest error at each component + /// + /// This is a standard error computation algorithm, but it's arguably bad if the state's components have different units. + /// It calculates the largest local estimate of the error from the integration (`error_est`) + /// given the difference in the candidate state and the previous state (`state_delta`). + /// This error estimator is from the physical model estimator of GMAT + /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/PhysicalModel.cpp#L987] + LargestError, + /// A largest state error control + /// + /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3018] + LargestState, + + /// A largest step error control which effectively computes the L1 norm of the provided Vector of size 3 + /// + /// Note that this error controller should be preferably be used only with slices of a state with the same units. + /// For example, one should probably use this for position independently of using it for the velocity. + /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3033] + LargestStep, +} + +impl ErrorControl { /// Computes the actual error of the current step. /// /// The `error_est` is the estimated error computed from the difference in the two stages of /// of the RK propagator. The `candidate` variable is the candidate state, and `cur_state` is /// the current state. This function must return the error. - fn estimate( - error_est: &OVector, - candidate: &OVector, - cur_state: &OVector, - ) -> f64 - where - DefaultAllocator: Allocator; -} - -/// A largest error control which effectively computes the largest error at each component -/// -/// This is a standard error computation algorithm, but it's arguably bad if the state's components have different units. -/// It calculates the largest local estimate of the error from the integration (`error_est`) -/// given the difference in the candidate state and the previous state (`state_delta`). -/// This error estimator is from the physical model estimator of GMAT -/// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/PhysicalModel.cpp#L987] -#[derive(Clone, Copy)] -pub struct LargestError; -impl ErrorCtrl for LargestError { - fn estimate( + pub fn estimate( + self, error_est: &OVector, candidate: &OVector, cur_state: &OVector, @@ -60,81 +84,111 @@ impl ErrorCtrl for LargestError { where DefaultAllocator: Allocator, { - let state_delta = candidate - cur_state; - let mut max_err = 0.0; - for (i, prop_err_i) in error_est.iter().enumerate() { - let err = if state_delta[i] > REL_ERR_THRESH { - (prop_err_i / state_delta[i]).abs() - } else { - prop_err_i.abs() - }; - if err > max_err { - max_err = err; + match self { + ErrorControl::RSSCartesianState => { + if N::dim() >= 6 { + let err_radius = RSSState::estimate::( + &error_est.fixed_rows::<3>(0).into_owned(), + &candidate.fixed_rows::<3>(0).into_owned(), + &cur_state.fixed_rows::<3>(0).into_owned(), + ); + let err_velocity = RSSState::estimate::( + &error_est.fixed_rows::<3>(3).into_owned(), + &candidate.fixed_rows::<3>(3).into_owned(), + &cur_state.fixed_rows::<3>(3).into_owned(), + ); + err_radius.max(err_velocity) + } else { + RSSStep::estimate(error_est, candidate, cur_state) + } } - } - max_err - } -} + ErrorControl::RSSCartesianStep => { + if N::dim() >= 6 { + let err_radius = RSSStep::estimate::( + &error_est.fixed_rows::<3>(0).into_owned(), + &candidate.fixed_rows::<3>(0).into_owned(), + &cur_state.fixed_rows::<3>(0).into_owned(), + ); + let err_velocity = RSSStep::estimate::( + &error_est.fixed_rows::<3>(3).into_owned(), + &candidate.fixed_rows::<3>(3).into_owned(), + &cur_state.fixed_rows::<3>(3).into_owned(), + ); + err_radius.max(err_velocity) + } else { + RSSStep::estimate(error_est, candidate, cur_state) + } + } + ErrorControl::RSSState => { + let mag = 0.5 * (candidate + cur_state).norm(); + let err = error_est.norm(); + if mag > REL_ERR_THRESH { + err / mag + } else { + err + } + } + ErrorControl::RSSStep => { + let mag = (candidate - cur_state).norm(); + let err = error_est.norm(); + if mag > REL_ERR_THRESH.sqrt() { + err / mag + } else { + err + } + } + ErrorControl::LargestError => { + let state_delta = candidate - cur_state; + let mut max_err = 0.0; + for (i, prop_err_i) in error_est.iter().enumerate() { + let err = if state_delta[i] > REL_ERR_THRESH { + (prop_err_i / state_delta[i]).abs() + } else { + prop_err_i.abs() + }; + if err > max_err { + max_err = err; + } + } + max_err + } + ErrorControl::LargestState => { + let sum_state = candidate + cur_state; + let mut mag = 0.0f64; + let mut err = 0.0f64; + for i in 0..N::dim() { + mag += 0.5 * sum_state[i].abs(); + err += error_est[i].abs(); + } -/// A largest step error control which effectively computes the L1 norm of the provided Vector of size 3 -/// -/// Note that this error controller should be preferably be used only with slices of a state with the same units. -/// For example, one should probably use this for position independently of using it for the velocity. -/// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3033] -#[derive(Clone, Copy)] -pub struct LargestStep; -impl ErrorCtrl for LargestStep { - fn estimate( - error_est: &OVector, - candidate: &OVector, - cur_state: &OVector, - ) -> f64 - where - DefaultAllocator: Allocator, - { - let state_delta = candidate - cur_state; - let mut mag = 0.0f64; - let mut err = 0.0f64; - for i in 0..N::dim() { - mag += state_delta[i].abs(); - err += error_est[i].abs(); - } + if mag > REL_ERR_THRESH { + err / mag + } else { + err + } + } + ErrorControl::LargestStep => { + let state_delta = candidate - cur_state; + let mut mag = 0.0f64; + let mut err = 0.0f64; + for i in 0..N::dim() { + mag += state_delta[i].abs(); + err += error_est[i].abs(); + } - if mag > REL_ERR_THRESH { - err / mag - } else { - err + if mag > REL_ERR_THRESH { + err / mag + } else { + err + } + } } } } -/// A largest state error control -/// -/// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3018] -#[derive(Clone, Copy)] -pub struct LargestState; -impl ErrorCtrl for LargestState { - fn estimate( - error_est: &OVector, - candidate: &OVector, - cur_state: &OVector, - ) -> f64 - where - DefaultAllocator: Allocator, - { - let sum_state = candidate + cur_state; - let mut mag = 0.0f64; - let mut err = 0.0f64; - for i in 0..N::dim() { - mag += 0.5 * sum_state[i].abs(); - err += error_est[i].abs(); - } - - if mag > REL_ERR_THRESH { - err / mag - } else { - err - } +impl Default for ErrorControl { + fn default() -> Self { + Self::RSSCartesianStep } } @@ -145,8 +199,8 @@ impl ErrorCtrl for LargestState { /// (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3045] #[derive(Clone, Copy)] #[allow(clippy::upper_case_acronyms)] -pub struct RSSStep; -impl ErrorCtrl for RSSStep { +struct RSSStep; +impl RSSStep { fn estimate( error_est: &OVector, candidate: &OVector, @@ -176,8 +230,8 @@ impl ErrorCtrl for RSSStep { /// > (Source)[https://github.com/ChristopherRabotin/GMAT/blob/37201a6290e7f7b941bc98ee973a527a5857104b/src/base/forcemodel/ODEModel.cpp#L3004] #[derive(Clone, Copy)] #[allow(clippy::upper_case_acronyms)] -pub struct RSSState; -impl ErrorCtrl for RSSState { +struct RSSState; +impl RSSState { fn estimate( error_est: &OVector, candidate: &OVector, @@ -195,126 +249,3 @@ impl ErrorCtrl for RSSState { } } } - -/// An RSS state error control which effectively for the provided vector -/// composed of two vectors of the same unit, both of size 3 (e.g. position + velocity). -#[derive(Clone, Copy)] -#[allow(clippy::upper_case_acronyms)] -pub struct RSSCartesianState; -impl ErrorCtrl for RSSCartesianState { - fn estimate( - error_est: &OVector, - candidate: &OVector, - cur_state: &OVector, - ) -> f64 - where - DefaultAllocator: Allocator, - { - if N::dim() >= 6 { - let err_radius = RSSState::estimate::( - &error_est.fixed_rows::<3>(0).into_owned(), - &candidate.fixed_rows::<3>(0).into_owned(), - &cur_state.fixed_rows::<3>(0).into_owned(), - ); - let err_velocity = RSSState::estimate::( - &error_est.fixed_rows::<3>(3).into_owned(), - &candidate.fixed_rows::<3>(3).into_owned(), - &cur_state.fixed_rows::<3>(3).into_owned(), - ); - let mut remaining_err = 0.0; - for i in 6..N::dim() { - let this_err = RSSState::estimate::( - &error_est.fixed_rows::<1>(i).into_owned(), - &candidate.fixed_rows::<1>(i).into_owned(), - &cur_state.fixed_rows::<1>(i).into_owned(), - ); - if this_err > remaining_err { - remaining_err = this_err; - } - } - remaining_err.max(err_radius.max(err_velocity)) - } else { - RSSState::estimate(error_est, candidate, cur_state) - } - } -} - -/// An RSS state error control which effectively for the provided vector -/// composed of two vectors of the same unit, both of size 3 (e.g. position + velocity). -#[derive(Clone, Copy)] -#[allow(clippy::upper_case_acronyms)] -pub struct RSSCartesianStep; -impl ErrorCtrl for RSSCartesianStep { - fn estimate( - error_est: &OVector, - candidate: &OVector, - cur_state: &OVector, - ) -> f64 - where - DefaultAllocator: Allocator, - { - if N::dim() >= 6 { - let err_radius = RSSStep::estimate::( - &error_est.fixed_rows::<3>(0).into_owned(), - &candidate.fixed_rows::<3>(0).into_owned(), - &cur_state.fixed_rows::<3>(0).into_owned(), - ); - let err_velocity = RSSStep::estimate::( - &error_est.fixed_rows::<3>(3).into_owned(), - &candidate.fixed_rows::<3>(3).into_owned(), - &cur_state.fixed_rows::<3>(3).into_owned(), - ); - err_radius.max(err_velocity) - } else { - RSSStep::estimate(error_est, candidate, cur_state) - } - } -} - -/// An RSS state error control which effectively for the provided vector -/// composed of two vectors of the same unit, both of size 3 (e.g. position + velocity). -#[derive(Clone, Copy)] -#[allow(clippy::upper_case_acronyms)] -pub struct RSSCartesianStepStm; -impl ErrorCtrl for RSSCartesianStepStm { - fn estimate( - error_est: &OVector, - candidate: &OVector, - cur_state: &OVector, - ) -> f64 - where - DefaultAllocator: Allocator, - { - let err_radius = RSSStep::estimate::( - &error_est.fixed_rows::<3>(0).into_owned(), - &candidate.fixed_rows::<3>(0).into_owned(), - &cur_state.fixed_rows::<3>(0).into_owned(), - ); - let err_velocity = RSSStep::estimate::( - &error_est.fixed_rows::<3>(3).into_owned(), - &candidate.fixed_rows::<3>(3).into_owned(), - &cur_state.fixed_rows::<3>(3).into_owned(), - ); - let err_cov_radius = RSSStep::estimate::( - &OVector::::new(error_est[6], error_est[6 + 7], error_est[6 + 14]), - &OVector::::new(candidate[6], candidate[6 + 7], candidate[6 + 14]), - &OVector::::new(cur_state[6], cur_state[6 + 7], cur_state[6 + 14]), - ); - - let err_cov_velocity = RSSStep::estimate::( - &OVector::::new(error_est[6 + 21], error_est[6 + 28], error_est[6 + 35]), - &OVector::::new(candidate[6 + 21], candidate[6 + 28], candidate[6 + 35]), - &OVector::::new(cur_state[6 + 21], cur_state[6 + 28], cur_state[6 + 35]), - ); - - let errs = vec![err_radius, err_velocity, err_cov_radius, err_cov_velocity]; - let mut max_err = 0.0; - for err in errs { - if err > max_err { - max_err = err; - } - } - - max_err - } -} diff --git a/src/propagators/instance.rs b/src/propagators/instance.rs index 47777dbe..26a770bf 100644 --- a/src/propagators/instance.rs +++ b/src/propagators/instance.rs @@ -16,9 +16,8 @@ along with this program. If not, see . */ -use super::error_ctrl::ErrorCtrl; use super::{DynamicsSnafu, IntegrationDetails, PropagationError, Propagator}; -use crate::dynamics::Dynamics; +use crate::dynamics::{Dynamics, DynamicsAlmanacSnafu}; use crate::linalg::allocator::Allocator; use crate::linalg::{DefaultAllocator, OVector}; use crate::md::trajectory::{Interpolatable, Traj}; @@ -39,7 +38,7 @@ use std::time::Instant; /// A Propagator allows propagating a set of dynamics forward or backward in time. /// It is an EventTracker, without any event tracking. It includes the options, the integrator /// details of the previous step, and the set of coefficients used for the monomorphic instance. -pub struct PropInstance<'a, D: Dynamics, E: ErrorCtrl> +pub struct PropInstance<'a, D: Dynamics> where DefaultAllocator: Allocator<::Size> + Allocator<::Size, ::Size> @@ -48,7 +47,7 @@ where /// The state of this propagator instance pub state: D::StateType, /// The propagator setup (kind, stages, etc.) - pub prop: &'a Propagator<'a, D, E>, + pub prop: &'a Propagator, /// Stores the details of the previous integration step pub details: IntegrationDetails, /// Should progress reports be logged @@ -60,7 +59,7 @@ where pub(crate) k: Vec::VecLength>>, } -impl<'a, D: Dynamics, E: ErrorCtrl> PropInstance<'a, D, E> +impl<'a, D: Dynamics> PropInstance<'a, D> where DefaultAllocator: Allocator<::Size> + Allocator<::Size, ::Size> @@ -95,11 +94,6 @@ where } let stop_time = self.state.epoch() + duration; - #[cfg(not(target_arch = "wasm32"))] - let tick = Instant::now(); - #[cfg(not(target_arch = "wasm32"))] - let mut prev_tick = Instant::now(); - if self.log_progress { // Prevent the print spam for orbit determination cases info!("Propagating for {} until {}", duration, stop_time); @@ -115,6 +109,39 @@ where if backprop { self.step_size = -self.step_size; // Invert the step size } + + // Transform the state if needed + let mut original_frame = None; + if let Some(integration_frame) = self.prop.opts.integration_frame { + if integration_frame != self.state.orbit().frame { + original_frame = Some(self.state.orbit().frame); + let mut new_orbit = self + .almanac + .transform_to(self.state.orbit(), integration_frame, None) + .context(DynamicsAlmanacSnafu { + action: "transforming state into desired integration frame", + }) + .context(DynamicsSnafu)?; + // If the integration frame has parameters, we set them here. + if let Some(mu_km3_s2) = integration_frame.mu_km3_s2 { + new_orbit.frame.mu_km3_s2 = Some(mu_km3_s2); + } + // If the integration frame has parameters, we set them here. + if let Some(shape) = integration_frame.shape { + new_orbit.frame.shape = Some(shape); + } + if self.log_progress { + info!("State transformed to the integration frame {integration_frame}"); + } + self.state.set_orbit(new_orbit); + } + } + + #[cfg(not(target_arch = "wasm32"))] + let tick = Instant::now(); + #[cfg(not(target_arch = "wasm32"))] + let mut prev_tick = Instant::now(); + loop { let epoch = self.state.epoch(); if (!backprop && epoch + self.step_size > stop_time) @@ -129,6 +156,19 @@ where debug!("Done in {}", tock); } } + + // Rotate back if needed + if let Some(original_frame) = original_frame { + let new_orbit = self + .almanac + .transform_to(self.state.orbit(), original_frame, None) + .context(DynamicsAlmanacSnafu { + action: "transforming state from desired integration frame", + }) + .context(DynamicsSnafu)?; + self.state.set_orbit(new_orbit); + } + return Ok(self.state); } // Take one final step of exactly the needed duration until the stop time @@ -160,6 +200,18 @@ where } } + // Rotate back if needed + if let Some(original_frame) = original_frame { + let new_orbit = self + .almanac + .transform_to(self.state.orbit(), original_frame, None) + .context(DynamicsAlmanacSnafu { + action: "transforming state from desired integration frame", + }) + .context(DynamicsSnafu)?; + self.state.set_orbit(new_orbit); + } + return Ok(self.state); } else { #[cfg(not(target_arch = "wasm32"))] @@ -342,14 +394,14 @@ where .context(DynamicsSnafu)?; self.k[0] = ki; let mut a_idx: usize = 0; - for i in 0..(self.prop.stages - 1) { + for i in 0..(self.prop.method.stages() - 1) { // Let's compute the c_i by summing the relevant items from the list of coefficients. // \sum_{j=1}^{i-1} a_ij ∀ i ∈ [2, s] let mut ci: f64 = 0.0; // The wi stores the a_{s1} * k_1 + a_{s2} * k_2 + ... + a_{s, s-1} * k_{s-1} + let mut wi = OVector::::VecLength>::from_element(0.0); for kj in &self.k[0..i + 1] { - let a_ij = self.prop.a_coeffs[a_idx]; + let a_ij = self.prop.method.a_coeffs()[a_idx]; ci += a_ij; wi += a_ij * kj; a_idx += 1; @@ -374,9 +426,9 @@ where let mut error_est = OVector::::VecLength>::from_element(0.0); for (i, ki) in self.k.iter().enumerate() { - let b_i = self.prop.b_coeffs[i]; + let b_i = self.prop.method.b_coeffs()[i]; if !self.fixed_step { - let b_i_star = self.prop.b_coeffs[i + self.prop.stages]; + let b_i_star = self.prop.method.b_coeffs()[i + self.prop.method.stages()]; error_est += step_size * (b_i - b_i_star) * ki; } next_state += step_size * b_i * ki; @@ -388,7 +440,11 @@ where return Ok(((self.details.step), next_state)); } else { // Compute the error estimate. - self.details.error = E::estimate(&error_est, &next_state, state_vec); + self.details.error = + self.prop + .opts + .error_ctrl + .estimate(&error_est, &next_state, state_vec); if self.details.error <= self.prop.opts.tolerance || step_size <= self.prop.opts.min_step.to_seconds() || self.details.attempts >= self.prop.opts.attempts @@ -407,7 +463,7 @@ where let proposed_step = 0.9 * step_size * (self.prop.opts.tolerance / self.details.error) - .powf(1.0 / f64::from(self.prop.order)); + .powf(1.0 / f64::from(self.prop.method.order())); step_size = if proposed_step > self.prop.opts.max_step.to_seconds() { self.prop.opts.max_step.to_seconds() } else { @@ -424,7 +480,7 @@ where let proposed_step = 0.9 * step_size * (self.prop.opts.tolerance / self.details.error) - .powf(1.0 / f64::from(self.prop.order - 1)); + .powf(1.0 / f64::from(self.prop.method.order() - 1)); step_size = if proposed_step < self.prop.opts.min_step.to_seconds() { self.prop.opts.min_step.to_seconds() } else { diff --git a/src/propagators/options.rs b/src/propagators/options.rs index 1c52e1ab..8ec73d49 100644 --- a/src/propagators/options.rs +++ b/src/propagators/options.rs @@ -20,19 +20,20 @@ use std::fmt; use crate::time::{Duration, Unit}; -use super::{ErrorCtrl, RSSCartesianStep}; +use super::ErrorControl; +use anise::frames::Frame; +use serde::{Deserialize, Serialize}; use typed_builder::TypedBuilder; -/// PropOpts stores the integrator options, including the minimum and maximum step sizes, and the -/// max error size. +/// Stores the integrator options, including the minimum and maximum step sizes, and the central body to perform the integration. /// /// Note that different step sizes and max errors are only used for adaptive /// methods. To use a fixed step integrator, initialize the options using `with_fixed_step`, and /// use whichever adaptive step integrator is desired. For example, initializing an RK45 with /// fixed step options will lead to an RK4 being used instead of an RK45. -#[derive(Clone, Copy, Debug, TypedBuilder)] +#[derive(Clone, Copy, Debug, TypedBuilder, Serialize, Deserialize, PartialEq)] #[builder(doc)] -pub struct PropOpts { +pub struct IntegratorOptions { #[builder(default_code = "60.0 * Unit::Second")] pub init_step: Duration, #[builder(default_code = "0.001 * Unit::Second")] @@ -45,19 +46,24 @@ pub struct PropOpts { pub attempts: u8, #[builder(default = false)] pub fixed_step: bool, - pub error_ctrl: E, + #[builder(default)] + pub error_ctrl: ErrorControl, + /// If a frame is specified and the propagator state is in a different frame, it it changed to this frame prior to integration. + /// Note, when setting this, it's recommended to call `strip` on the Frame. + #[builder(default, setter(strip_option))] + pub integration_frame: Option, } -impl PropOpts { +impl IntegratorOptions { /// `with_adaptive_step` initializes an `PropOpts` such that the integrator is used with an /// adaptive step size. The number of attempts is currently fixed to 50 (as in GMAT). pub fn with_adaptive_step( min_step: Duration, max_step: Duration, tolerance: f64, - error_ctrl: E, + error_ctrl: ErrorControl, ) -> Self { - PropOpts { + IntegratorOptions { init_step: max_step, min_step, max_step, @@ -65,6 +71,7 @@ impl PropOpts { attempts: 50, fixed_step: false, error_ctrl, + integration_frame: None, } } @@ -72,7 +79,7 @@ impl PropOpts { min_step: f64, max_step: f64, tolerance: f64, - error_ctrl: E, + error_ctrl: ErrorControl, ) -> Self { Self::with_adaptive_step( min_step * Unit::Second, @@ -82,6 +89,41 @@ impl PropOpts { ) } + /// `with_fixed_step` initializes an `PropOpts` such that the integrator is used with a fixed + /// step size. + pub fn with_fixed_step(step: Duration) -> Self { + IntegratorOptions { + init_step: step, + min_step: step, + max_step: step, + tolerance: 0.0, + fixed_step: true, + attempts: 0, + error_ctrl: ErrorControl::RSSCartesianStep, + integration_frame: None, + } + } + + pub fn with_fixed_step_s(step: f64) -> Self { + Self::with_fixed_step(step * Unit::Second) + } + + /// Returns the default options with a specific tolerance. + #[allow(clippy::field_reassign_with_default)] + pub fn with_tolerance(tolerance: f64) -> Self { + let mut opts = Self::default(); + opts.tolerance = tolerance; + opts + } + + /// Creates a propagator with the provided max step, and sets the initial step to that value as well. + #[allow(clippy::field_reassign_with_default)] + pub fn with_max_step(max_step: Duration) -> Self { + let mut opts = Self::default(); + opts.set_max_step(max_step); + opts + } + /// Returns a string with the information about these options pub fn info(&self) -> String { format!("{self}") @@ -104,7 +146,7 @@ impl PropOpts { } } -impl fmt::Display for PropOpts { +impl fmt::Display for IntegratorOptions { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { if self.fixed_step { write!(f, "fixed step: {:e}", self.min_step,) @@ -118,86 +160,66 @@ impl fmt::Display for PropOpts { } } -impl PropOpts { - /// `with_fixed_step` initializes an `PropOpts` such that the integrator is used with a fixed - /// step size. - pub fn with_fixed_step(step: Duration) -> Self { - PropOpts { - init_step: step, - min_step: step, - max_step: step, - tolerance: 0.0, - fixed_step: true, - attempts: 0, - error_ctrl: RSSCartesianStep {}, - } - } - - pub fn with_fixed_step_s(step: f64) -> Self { - Self::with_fixed_step(step * Unit::Second) - } - - /// Returns the default options with a specific tolerance. - #[allow(clippy::field_reassign_with_default)] - pub fn with_tolerance(tolerance: f64) -> Self { - let mut opts = Self::default(); - opts.tolerance = tolerance; - opts - } - - /// Creates a propagator with the provided max step, and sets the initial step to that value as well. - #[allow(clippy::field_reassign_with_default)] - pub fn with_max_step(max_step: Duration) -> Self { - let mut opts = Self::default(); - opts.set_max_step(max_step); - opts - } -} - -impl Default for PropOpts { +impl Default for IntegratorOptions { /// `default` returns the same default options as GMAT. - fn default() -> PropOpts { - PropOpts { + fn default() -> IntegratorOptions { + IntegratorOptions { init_step: 60.0 * Unit::Second, min_step: 0.001 * Unit::Second, max_step: 2700.0 * Unit::Second, tolerance: 1e-12, attempts: 50, fixed_step: false, - error_ctrl: RSSCartesianStep {}, + error_ctrl: ErrorControl::RSSCartesianStep, + integration_frame: None, } } } -#[test] -fn test_options() { - use super::error_ctrl::RSSStep; - - let opts = PropOpts::with_fixed_step_s(1e-1); - assert_eq!(opts.min_step, 1e-1 * Unit::Second); - assert_eq!(opts.max_step, 1e-1 * Unit::Second); - assert!(opts.tolerance.abs() < f64::EPSILON); - assert!(opts.fixed_step); - - let opts = PropOpts::with_adaptive_step_s(1e-2, 10.0, 1e-12, RSSStep {}); - assert_eq!(opts.min_step, 1e-2 * Unit::Second); - assert_eq!(opts.max_step, 10.0 * Unit::Second); - assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); - assert!(!opts.fixed_step); - - let opts: PropOpts = Default::default(); - assert_eq!(opts.init_step, 60.0 * Unit::Second); - assert_eq!(opts.min_step, 0.001 * Unit::Second); - assert_eq!(opts.max_step, 2700.0 * Unit::Second); - assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); - assert_eq!(opts.attempts, 50); - assert!(!opts.fixed_step); - - let opts = PropOpts::with_max_step(1.0 * Unit::Second); - assert_eq!(opts.init_step, 1.0 * Unit::Second); - assert_eq!(opts.min_step, 0.001 * Unit::Second); - assert_eq!(opts.max_step, 1.0 * Unit::Second); - assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); - assert_eq!(opts.attempts, 50); - assert!(!opts.fixed_step); +#[cfg(test)] +mod ut_integr_opts { + use hifitime::Unit; + + use crate::propagators::{ErrorControl, IntegratorOptions}; + + #[test] + fn test_options() { + let opts = IntegratorOptions::with_fixed_step_s(1e-1); + assert_eq!(opts.min_step, 1e-1 * Unit::Second); + assert_eq!(opts.max_step, 1e-1 * Unit::Second); + assert!(opts.tolerance.abs() < f64::EPSILON); + assert!(opts.fixed_step); + + let opts = + IntegratorOptions::with_adaptive_step_s(1e-2, 10.0, 1e-12, ErrorControl::RSSStep); + assert_eq!(opts.min_step, 1e-2 * Unit::Second); + assert_eq!(opts.max_step, 10.0 * Unit::Second); + assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); + assert!(!opts.fixed_step); + + let opts: IntegratorOptions = Default::default(); + assert_eq!(opts.init_step, 60.0 * Unit::Second); + assert_eq!(opts.min_step, 0.001 * Unit::Second); + assert_eq!(opts.max_step, 2700.0 * Unit::Second); + assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); + assert_eq!(opts.attempts, 50); + assert!(!opts.fixed_step); + + let opts = IntegratorOptions::with_max_step(1.0 * Unit::Second); + assert_eq!(opts.init_step, 1.0 * Unit::Second); + assert_eq!(opts.min_step, 0.001 * Unit::Second); + assert_eq!(opts.max_step, 1.0 * Unit::Second); + assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); + assert_eq!(opts.attempts, 50); + assert!(!opts.fixed_step); + } + + #[test] + fn test_serde() { + let opts = IntegratorOptions::default(); + let serialized = toml::to_string(&opts).unwrap(); + println!("{serialized}"); + let deserd: IntegratorOptions = toml::from_str(&serialized).unwrap(); + assert_eq!(deserd, opts); + } } diff --git a/src/propagators/propagator.rs b/src/propagators/propagator.rs index 1ab9e8cc..0ed595da 100644 --- a/src/propagators/propagator.rs +++ b/src/propagators/propagator.rs @@ -20,8 +20,7 @@ use std::sync::Arc; use anise::almanac::Almanac; -use super::error_ctrl::{ErrorCtrl, RSSCartesianStep}; -use super::{Dormand78, IntegrationDetails, PropInstance, PropOpts, RK, RK89}; +use super::{IntegrationDetails, IntegratorMethod, IntegratorOptions, PropInstance}; use crate::dynamics::Dynamics; use crate::linalg::allocator::Allocator; use crate::linalg::{DefaultAllocator, OVector}; @@ -32,7 +31,7 @@ use crate::State; /// It is an EventTracker, without any event tracking. It includes the options, the integrator /// details of the previous step, and the set of coefficients used for the monomorphic instance. #[derive(Clone, Debug)] -pub struct Propagator<'a, D: Dynamics, E: ErrorCtrl> +pub struct Propagator where DefaultAllocator: Allocator<::Size> + Allocator<::Size, ::Size> @@ -40,15 +39,12 @@ where + Allocator<::VecLength>, { pub dynamics: D, // Stores the dynamics used. *Must* use this to get the latest values - pub opts: PropOpts, // Stores the integration options (tolerance, min/max step, init step, etc.) - pub(crate) order: u8, // Order of the integrator - pub(crate) stages: usize, // Number of stages, i.e. how many times the derivatives will be called - pub(crate) a_coeffs: &'a [f64], - pub(crate) b_coeffs: &'a [f64], + pub opts: IntegratorOptions, // Stores the integration options (tolerance, min/max step, init step, etc.) + pub method: IntegratorMethod, } /// The `Propagator` trait defines the functions of a propagator and of an event tracker. -impl<'a, D: Dynamics, E: ErrorCtrl> Propagator<'a, D, E> +impl Propagator where DefaultAllocator: Allocator<::Size> + Allocator<::Size, ::Size> @@ -56,14 +52,11 @@ where + Allocator<::VecLength>, { /// Each propagator must be initialized with `new` which stores propagator information. - pub fn new(dynamics: D, opts: PropOpts) -> Self { + pub fn new(dynamics: D, method: IntegratorMethod, opts: IntegratorOptions) -> Self { Self { dynamics, opts, - stages: T::STAGES, - order: T::ORDER, - a_coeffs: T::A_COEFFS, - b_coeffs: T::B_COEFFS, + method, } } @@ -82,20 +75,20 @@ where } /// An RK89 propagator (the default) with custom propagator options. - pub fn rk89(dynamics: D, opts: PropOpts) -> Self { - Self::new::(dynamics, opts) + pub fn rk89(dynamics: D, opts: IntegratorOptions) -> Self { + Self::new(dynamics, IntegratorMethod::RungeKutta89, opts) } /// A Dormand Prince 7-8 propagator with custom propagator options: it's about 20% faster than an RK98, and more stable in two body dynamics. /// WARNINGS: Dormand Prince may have issues with generating proper trajectories, leading to glitches in event finding. - pub fn dp78(dynamics: D, opts: PropOpts) -> Self { - Self::new::(dynamics, opts) + pub fn dp78(dynamics: D, opts: IntegratorOptions) -> Self { + Self::new(dynamics, IntegratorMethod::DormandPrince78, opts) } - pub fn with(&'a self, state: D::StateType, almanac: Arc) -> PropInstance<'a, D, E> { + pub fn with(&self, state: D::StateType, almanac: Arc) -> PropInstance { // Pre-allocate the k used in the propagator - let mut k = Vec::with_capacity(self.stages + 1); - for _ in 0..self.stages { + let mut k = Vec::with_capacity(self.method.stages() + 1); + for _ in 0..self.method.stages() { k.push(OVector::::VecLength>::zeros()); } PropInstance { @@ -113,24 +106,16 @@ where k, } } -} -impl<'a, D: Dynamics> Propagator<'a, D, RSSCartesianStep> -where - DefaultAllocator: Allocator<::Size> - + Allocator<::Size, ::Size> - + Allocator<::Size, ::Size> - + Allocator<::VecLength>, -{ /// Default propagator is an RK89 with the default PropOpts. pub fn default(dynamics: D) -> Self { - Self::new::(dynamics, PropOpts::default()) + Self::rk89(dynamics, IntegratorOptions::default()) } /// A default Dormand Prince 78 propagator with the default PropOpts. /// Faster and more stable than an RK89 (`default`) but seems to cause issues for event finding. /// WARNINGS: Dormand Prince may have issues with generating proper trajectories, leading to glitches in event finding. pub fn default_dp78(dynamics: D) -> Self { - Self::new::(dynamics, PropOpts::default()) + Self::dp78(dynamics, IntegratorOptions::default()) } } diff --git a/src/propagators/rk_methods/dormand.rs b/src/propagators/rk_methods/dormand.rs index b124c96d..31948135 100644 --- a/src/propagators/rk_methods/dormand.rs +++ b/src/propagators/rk_methods/dormand.rs @@ -19,7 +19,7 @@ use super::RK; /// `Dormand45` is a [Dormand-Prince integrator](https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method). -pub struct Dormand45 {} +pub(crate) struct Dormand45 {} impl RK for Dormand45 { const ORDER: u8 = 5; @@ -68,7 +68,7 @@ impl RK for Dormand45 { /// `Dormand78` is a [Dormand-Prince integrator](https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method). /// /// Coefficients taken from GMAT `src/base/propagator/PrinceDormand78.cpp`. -pub struct Dormand78 {} +pub(crate) struct Dormand78 {} impl RK for Dormand78 { const ORDER: u8 = 8; diff --git a/src/propagators/rk_methods/fehlberg.rs b/src/propagators/rk_methods/fehlberg.rs deleted file mode 100644 index 7c358bfa..00000000 --- a/src/propagators/rk_methods/fehlberg.rs +++ /dev/null @@ -1,58 +0,0 @@ -/* - Nyx, blazing fast astrodynamics - Copyright (C) 2018-onwards Christopher Rabotin - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU Affero General Public License as published - by the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Affero General Public License for more details. - - You should have received a copy of the GNU Affero General Public License - along with this program. If not, see . -*/ - -use super::RK; - -/// `Fehlberg45` is a [Runge Kutta Fehlberg integrator](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method). -pub struct Fehlberg45 {} - -impl RK for Fehlberg45 { - const ORDER: u8 = 5; - const STAGES: usize = 6; - const A_COEFFS: &'static [f64] = &[ - 1.0 / 4.0, - 3.0 / 32.0, - 9.0 / 32.0, - 1932.0 / 2197.0, - -7200.0 / 2197.0, - 7296.0 / 2197.0, - 439.0 / 216.0, - -8.0, - 3680.0 / 513.0, - -845.0 / 4104.0, - -8.0 / 27.0, - 2.0, - -3544.0 / 2565.0, - 1859.0 / 4104.0, - -11.0 / 40.0, - ]; - const B_COEFFS: &'static [f64] = &[ - 16.0 / 135.0, - 0.0, - 6656.0 / 12825.0, - 28561.0 / 56430.0, - -9.0 / 50.0, - 2.0 / 55.0, - 25.0 / 216.0, - 0.0, - 1408.0 / 2565.0, - 2197.0 / 4104.0, - -1.0 / 5.0, - 0.0, - ]; -} diff --git a/src/propagators/rk_methods/mod.rs b/src/propagators/rk_methods/mod.rs index 096474a8..f6500385 100644 --- a/src/propagators/rk_methods/mod.rs +++ b/src/propagators/rk_methods/mod.rs @@ -17,17 +17,24 @@ */ mod rk; -pub use self::rk::*; +use std::str::FromStr; + +use serde::Deserialize; +use serde::Serialize; + +use crate::io::ConfigError; + +use self::rk::*; mod dormand; -pub use self::dormand::*; -mod fehlberg; -pub use self::fehlberg::*; +use self::dormand::*; mod verner; -pub use self::verner::*; +use self::verner::*; + +use super::PropagationError; /// The `RK` trait defines a Runge Kutta integrator. #[allow(clippy::upper_case_acronyms)] -pub trait RK +trait RK where Self: Sized, { @@ -47,3 +54,135 @@ where /// Butcher table for that RK. `Self.a_coeffs().len()` must be of size (order+1)*2. const B_COEFFS: &'static [f64]; } + +/// Enum of supported integration methods, all of which are part of the Runge Kutta family of ordinary differential equation (ODE) solvers. +/// Nomenclature: X-Y means that this is an X order solver with a Y order error correction step. +#[derive(Copy, Clone, Debug, PartialEq, Eq, Serialize, Deserialize)] +pub enum IntegratorMethod { + /// Runge Kutta 8-9 is the recommended integrator for most application. + RungeKutta89, + /// `Dormand78` is a [Dormand-Prince integrator](https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method). Coefficients taken from GMAT `src/base/propagator/PrinceDormand78.cpp`. + DormandPrince78, + /// `Dormand45` is a [Dormand-Prince integrator](https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method). + DormandPrince45, + /// Runge Kutta 4 is a fixed step solver. + RungeKutta4, + /// Runge Kutta 4-5 [Cash Karp integrator](https://en.wikipedia.org/wiki/Cash%E2%80%93Karp_method). + CashKarp45, + /// Verner56 is an RK Verner integrator of order 5-6. Coefficients taken from [here (PDF)](http://people.math.sfu.ca/~jverner/classify.1992.ps). + Verner56, +} + +impl IntegratorMethod { + /// Returns the order of this integrator (as u8 because there probably isn't an order greater than 255). + /// The order is used for the adaptive step size only to compute the error between estimates. + pub const fn order(self) -> u8 { + match self { + Self::RungeKutta89 => RK89::ORDER, + Self::DormandPrince78 => Dormand78::ORDER, + Self::DormandPrince45 => Dormand45::ORDER, + Self::RungeKutta4 => RK4Fixed::ORDER, + Self::CashKarp45 => CashKarp45::ORDER, + Self::Verner56 => Verner56::ORDER, + } + } + + /// Returns the stages of this integrator, i.e. how many times the derivatives will be called + pub const fn stages(self) -> usize { + match self { + Self::RungeKutta89 => RK89::STAGES, + Self::DormandPrince78 => Dormand78::STAGES, + Self::DormandPrince45 => Dormand45::STAGES, + Self::RungeKutta4 => RK4Fixed::STAGES, + Self::CashKarp45 => CashKarp45::STAGES, + Self::Verner56 => Verner56::STAGES, + } + } + + /// Returns a pointer to a list of f64 corresponding to the A coefficients of the Butcher table for that RK. + /// This module only supports *implicit* integrators, and as such, `Self.a_coeffs().len()` must be of + /// size (order+1)*(order)/2. + /// *Warning:* this RK trait supposes that the implementation is consistent, i.e. c_i = \sum_j a_{ij}. + pub const fn a_coeffs(self) -> &'static [f64] { + match self { + Self::RungeKutta89 => RK89::A_COEFFS, + Self::DormandPrince78 => Dormand78::A_COEFFS, + Self::DormandPrince45 => Dormand45::A_COEFFS, + Self::RungeKutta4 => RK4Fixed::A_COEFFS, + Self::CashKarp45 => CashKarp45::A_COEFFS, + Self::Verner56 => Verner56::A_COEFFS, + } + } + /// Returns a pointer to a list of f64 corresponding to the b_i and b^*_i coefficients of the + /// Butcher table for that RK. `Self.a_coeffs().len()` must be of size (order+1)*2. + pub const fn b_coeffs(self) -> &'static [f64] { + match self { + Self::RungeKutta89 => RK89::B_COEFFS, + Self::DormandPrince78 => Dormand78::B_COEFFS, + Self::DormandPrince45 => Dormand45::B_COEFFS, + Self::RungeKutta4 => RK4Fixed::B_COEFFS, + Self::CashKarp45 => CashKarp45::B_COEFFS, + Self::Verner56 => Verner56::B_COEFFS, + } + } +} + +impl Default for IntegratorMethod { + fn default() -> Self { + Self::RungeKutta89 + } +} + +impl FromStr for IntegratorMethod { + type Err = PropagationError; + + fn from_str(s: &str) -> Result { + match s.to_lowercase().as_str() { + "rungekutta89" => Ok(Self::RungeKutta89), + "dormandprince78" => Ok(Self::DormandPrince78), + "dormandprince45" => Ok(Self::DormandPrince45), + "rungekutta4" => Ok(Self::RungeKutta4), + "cashkarp45" => Ok(Self::CashKarp45), + "verner56" => Ok(Self::Verner56), + _ => { + let valid = [ + "RungeKutta89", + "DormandPrince78", + "DormandPrince45", + "RungeKutta4", + "CashKarp45", + "Verner56", + ]; + let valid_msg = valid.join(","); + Err(PropagationError::PropConfigError { + source: ConfigError::InvalidConfig { + msg: format!("unknow integration method `{s}`, must be one of {valid_msg}"), + }, + }) + } + } + } +} + +#[cfg(test)] +mod ut_propagator { + use std::str::FromStr; + + use super::IntegratorMethod; + + #[test] + fn from_str_ok() { + let valid = [ + "RungeKutta89", + "DormandPrince78", + "DormandPrince45", + "RungeKutta4", + "CashKarp45", + "Verner56", + ]; + for method in valid { + assert!(IntegratorMethod::from_str(method.to_uppercase().as_str()).is_ok()); + } + assert!(IntegratorMethod::from_str("blah").is_err()); + } +} diff --git a/src/propagators/rk_methods/rk.rs b/src/propagators/rk_methods/rk.rs index 5cd3b35b..282b64cc 100644 --- a/src/propagators/rk_methods/rk.rs +++ b/src/propagators/rk_methods/rk.rs @@ -16,10 +16,10 @@ along with this program. If not, see . */ -pub use super::RK; +use super::RK; /// `CashKarp45` is a [Runge Kutta Cash Karp integrator](https://en.wikipedia.org/wiki/Cash%E2%80%93Karp_method). -pub struct CashKarp45 {} +pub(crate) struct CashKarp45 {} impl RK for CashKarp45 { const ORDER: u8 = 5; @@ -61,7 +61,7 @@ impl RK for CashKarp45 { /// /// If initialized with an `PropOpts.with_adaptive_step`, the variable step will **not** be taken into consideration. #[allow(clippy::upper_case_acronyms)] -pub struct RK4Fixed {} +pub(crate) struct RK4Fixed {} impl RK for RK4Fixed { const ORDER: u8 = 4; @@ -80,32 +80,13 @@ impl RK for RK4Fixed { ]; } -/// `RK2Fixed` is a fixed step RK4 (or midpoint method). -/// -/// If initialized with an `PropOpts.with_adaptive_step`, the variable step will **not** be taken into consideration. -#[allow(clippy::upper_case_acronyms)] -pub struct RK2Fixed {} - -impl RK for RK2Fixed { - const ORDER: u8 = 2; - const STAGES: usize = 2; - const A_COEFFS: &'static [f64] = &[2.0 / 3.0]; - const B_COEFFS: &'static [f64] = &[ - 1.0 / 4.0, - 3.0 / 4.0, - // NOTE: Duplicating the B coefficients for force the error to zero. - 1.0 / 4.0, - 3.0 / 4.0, - ]; -} - const SQRT6: f64 = 2.449_489_742_783_178; /// `RK89` is a Runge Kutta 8-9 integrator. /// /// Coefficients taken from GMAT `src/base/propagator/RungeKutta89.cpp`. #[allow(clippy::upper_case_acronyms)] -pub struct RK89 {} +pub(crate) struct RK89 {} impl RK for RK89 { const ORDER: u8 = 9; diff --git a/src/propagators/rk_methods/verner.rs b/src/propagators/rk_methods/verner.rs index 9318d93a..0adb5dbe 100644 --- a/src/propagators/rk_methods/verner.rs +++ b/src/propagators/rk_methods/verner.rs @@ -21,7 +21,7 @@ use super::RK; /// `Verner56` is an RK Verner integrator of order 5-6. /// /// Coefficients taken from [here (PDF)](http://people.math.sfu.ca/~jverner/classify.1992.ps). -pub struct Verner56 {} +pub(crate) struct Verner56 {} impl RK for Verner56 { const ORDER: u8 = 6; diff --git a/tests/cosmic/eclipse.rs b/tests/cosmic/eclipse.rs index e9cdb84c..92755df0 100644 --- a/tests/cosmic/eclipse.rs +++ b/tests/cosmic/eclipse.rs @@ -7,7 +7,7 @@ use nyx::cosmic::eclipse::EclipseLocator; use nyx::cosmic::Orbit; use nyx::dynamics::orbital::OrbitalDynamics; use nyx::dynamics::SpacecraftDynamics; -use nyx::propagators::{PropOpts, Propagator}; +use nyx::propagators::{IntegratorOptions, Propagator}; use nyx::time::{Epoch, Unit}; use std::sync::{mpsc, Arc}; use std::thread; @@ -38,7 +38,7 @@ fn leo_sun_earth_eclipses(almanac: Arc) { let almanac_c = almanac.clone(); thread::spawn(move || { let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::rk89(dynamics, PropOpts::with_fixed_step_s(60.0)); + let setup = Propagator::rk89(dynamics, IntegratorOptions::with_fixed_step_s(60.0)); setup .with(leo.into(), almanac_c) @@ -90,7 +90,7 @@ fn geo_sun_earth_eclipses(almanac: Arc) { thread::spawn(move || { let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::rk89(dynamics, PropOpts::with_fixed_step_s(60.0)); + let setup = Propagator::rk89(dynamics, IntegratorOptions::with_fixed_step_s(60.0)); setup .with(geo.into(), almanac_c) diff --git a/tests/mission_design/force_models.rs b/tests/mission_design/force_models.rs index 41fcf65c..4139c793 100644 --- a/tests/mission_design/force_models.rs +++ b/tests/mission_design/force_models.rs @@ -1,5 +1,6 @@ extern crate nyx_space as nyx; +use anise::constants::frames::MOON_J2000; use nyx::cosmic::{Orbit, Spacecraft}; use nyx::dynamics::{Drag, OrbitalDynamics, SolarPressure, SpacecraftDynamics}; use nyx::linalg::Vector6; @@ -8,6 +9,7 @@ use nyx::time::{Epoch, Unit}; use nyx::utils::rss_orbit_vec_errors; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; +use nyx_space::md::prelude::Interpolatable; use rstest::*; use std::sync::Arc; @@ -150,7 +152,7 @@ fn srp_earth_meo_ecc_inc(almanac: Arc) { let sc = Spacecraft::from_srp_defaults(orbit, dry_mass, 16.0); - let setup = Propagator::default(sc_dyn); + let setup = Propagator::default(sc_dyn.clone()); let mut prop = setup.with(sc, almanac.clone()); let final_state = prop.for_duration(prop_time).unwrap(); @@ -176,6 +178,39 @@ fn srp_earth_meo_ecc_inc(almanac: Arc) { assert!(err_r < 2e-3, "position error too large for SRP"); assert!(err_v < 1e-6, "velocity error too large for SRP"); + // Test that we can specify an integration frame in the options and that the result is the same. + // Note that we also test here that we're setting the GM and shape of the integration frame as defined + // and not just grabbing that data from the almanac. + let orbit = almanac.transform_to(orbit, MOON_J2000, None).unwrap(); + println!("{:x}", orbit); + + let sc_moon = Spacecraft::from_srp_defaults(orbit, dry_mass, 16.0); + + let mut setup_moon = Propagator::default(sc_dyn); + setup_moon.opts.integration_frame = Some(eme2k); + let mut prop = setup_moon.with(sc_moon, almanac.clone()); + let final_moon_state = prop.for_duration(prop_time).unwrap(); + assert!( + final_moon_state.frame().ephem_origin_match(MOON_J2000), + "expected a result in the Moon frame" + ); + println!("{}", final_moon_state); + + let final_earth_orbit = almanac + .transform_to(final_moon_state.orbit, EARTH_J2000, None) + .unwrap(); + + let (fw_err_r, fw_err_v) = + rss_orbit_vec_errors(&final_earth_orbit.to_cartesian_pos_vel(), &rslt); + println!( + "Error accumulated in ecc+inc MEO (with penumbras) over {} after integration frame swap: {:.6} m \t{:.6} m/s", + prop_time, + fw_err_r * 1e3, + fw_err_v * 1e3 + ); + assert!(dbg!((fw_err_r - err_r).abs() / err_r) < 1e-9); + assert!(dbg!((fw_err_v - err_v).abs() / err_v) < 1e-12); + // Compare the case with the hyperdual EOMs (computation uses another part of the code) let mut prop = setup.with(sc.with_stm(), almanac); let final_state_dual = prop.for_duration(prop_time).unwrap(); diff --git a/tests/mission_design/orbitaldyn.rs b/tests/mission_design/orbitaldyn.rs index abd0a905..3da4ad4b 100644 --- a/tests/mission_design/orbitaldyn.rs +++ b/tests/mission_design/orbitaldyn.rs @@ -8,7 +8,6 @@ use na::{Const, OMatrix}; use nyx::cosmic::{assert_orbit_eq_or_abs, Orbit}; use nyx::dynamics::{Dynamics, OrbitalDynamics, PointMasses, SpacecraftDynamics}; use nyx::linalg::Vector6; -use nyx::propagators::error_ctrl::RSSCartesianStep; use nyx::time::{Epoch, Unit}; use nyx::utils::{rss_orbit_errors, rss_orbit_vec_errors}; use nyx::State; @@ -53,9 +52,10 @@ fn energy_conservation(almanac: Arc) { eme2k, ); - let rk89_final = Propagator::new::( + let rk89_final = Propagator::new( SpacecraftDynamics::new(OrbitalDynamics::two_body()), - PropOpts::default(), + IntegratorMethod::RungeKutta89, + IntegratorOptions::default(), ) .with(Spacecraft::from(start_state), almanac.clone()) .for_duration(prop_time) @@ -74,9 +74,10 @@ fn energy_conservation(almanac: Arc) { } println!(); - let dp78_final = Propagator::new::( + let dp78_final = Propagator::new( SpacecraftDynamics::new(OrbitalDynamics::two_body()), - PropOpts::default(), + IntegratorMethod::DormandPrince78, + IntegratorOptions::default(), ) .with(start_state.into(), almanac) .for_duration(prop_time) @@ -120,7 +121,7 @@ fn val_two_body_dynamics(almanac: Arc) { ); let dynamics = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - let setup = Propagator::rk89(dynamics, PropOpts::::default()); + let setup = Propagator::rk89(dynamics, IntegratorOptions::default()); let mut prop = setup.with(state.into(), almanac); prop.for_duration(prop_time).unwrap(); @@ -209,7 +210,10 @@ fn val_halo_earth_moon_dynamics(almanac_gmat: Arc) { let bodies = vec![MOON]; let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::rk89(dynamics, PropOpts::with_fixed_step(10 * Unit::Second)); + let setup = Propagator::rk89( + dynamics, + IntegratorOptions::with_fixed_step(10 * Unit::Second), + ); let mut prop = setup.with(halo_rcvr.into(), almanac); prop.for_duration(prop_time).unwrap(); let (err_r, err_v) = rss_orbit_errors(&prop.state.orbit, &rslt); @@ -273,7 +277,7 @@ fn val_halo_earth_moon_dynamics_adaptive(almanac_gmat: Arc) { let bodies = vec![MOON]; let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::rk89(dynamics, PropOpts::default()); + let setup = Propagator::rk89(dynamics, IntegratorOptions::default()); let mut prop = setup.with(halo_rcvr.into(), almanac); prop.for_duration(prop_time).unwrap(); let (err_r, err_v) = rss_orbit_vec_errors(&prop.state.orbit.to_cartesian_pos_vel(), &rslt); @@ -339,7 +343,7 @@ fn val_llo_earth_moon_dynamics_adaptive(almanac_gmat: Arc) { let bodies = vec![MOON]; let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::rk89(dynamics, PropOpts::default()); + let setup = Propagator::rk89(dynamics, IntegratorOptions::default()); let mut prop = setup.with(llo_xmtr.into(), almanac); prop.for_duration(prop_time).unwrap(); let (err_r, err_v) = rss_orbit_vec_errors(&prop.state.orbit.to_cartesian_pos_vel(), &rslt); @@ -401,7 +405,10 @@ fn val_halo_multi_body_dynamics(almanac_gmat: Arc) { let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::rk89(dynamics, PropOpts::with_fixed_step(10 * Unit::Second)); + let setup = Propagator::rk89( + dynamics, + IntegratorOptions::with_fixed_step(10 * Unit::Second), + ); let mut prop = setup.with(halo_rcvr.into(), almanac); prop.for_duration(prop_time).unwrap(); let (err_r, err_v) = rss_orbit_vec_errors(&prop.state.orbit.to_cartesian_pos_vel(), &rslt); @@ -752,7 +759,7 @@ fn two_body_dual(almanac: Arc) { let prop_time = 2 * Unit::Minute; let step_size = 10 * Unit::Second; - let setup = Propagator::rk89(dynamics, PropOpts::with_fixed_step(step_size)); + let setup = Propagator::rk89(dynamics, IntegratorOptions::with_fixed_step(step_size)); let mut prop = setup.with(init_sc, almanac.clone()); let final_state = prop.for_duration(prop_time).unwrap(); @@ -805,7 +812,7 @@ fn multi_body_dynamics_dual(almanac: Arc) { // let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); let dynamics = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - let setup = Propagator::rk89(dynamics, PropOpts::with_fixed_step(step_size)); + let setup = Propagator::rk89(dynamics, IntegratorOptions::with_fixed_step(step_size)); let final_state = setup .with(halo_rcvr.into(), almanac.clone()) .for_duration(prop_time) @@ -915,7 +922,7 @@ fn val_earth_sph_harmonics_j2(almanac: Arc) { let dynamics = SpacecraftDynamics::new(OrbitalDynamics::from_model(harmonics)); - let prop_state = Propagator::rk89(dynamics, PropOpts::::default()) + let prop_state = Propagator::rk89(dynamics, IntegratorOptions::default()) .with(state.into(), almanac) .for_duration(1 * Unit::Day) .unwrap(); @@ -968,7 +975,7 @@ fn val_earth_sph_harmonics_12x12(almanac_gmat: Arc) { let dynamics = SpacecraftDynamics::new(OrbitalDynamics::from_model(harmonics)); - let setup = Propagator::rk89(dynamics.clone(), PropOpts::with_tolerance(1e-9)); + let setup = Propagator::rk89(dynamics.clone(), IntegratorOptions::with_tolerance(1e-9)); let prop_time = 1 * Unit::Day; let final_state = setup .with(state.into(), almanac.clone()) @@ -992,7 +999,7 @@ fn val_earth_sph_harmonics_12x12(almanac_gmat: Arc) { // We set up a new propagator with a fixed step. Without the fixed step, the error control // on the STM leads to a difference of 1.04 meters in this one day propagation. - let setup = Propagator::rk89(dynamics, PropOpts::with_fixed_step_s(30.0)); + let setup = Propagator::rk89(dynamics, IntegratorOptions::with_fixed_step_s(30.0)); let prop_time = 6 * Unit::Hour; let final_state = setup .with(state.into(), almanac.clone()) @@ -1155,7 +1162,7 @@ fn hf_prop(almanac: Arc) { harmonics, ])); - let setup = Propagator::rk89(dynamics, PropOpts::with_tolerance(1e-9)); + let setup = Propagator::rk89(dynamics, IntegratorOptions::with_tolerance(1e-9)); let rslt = setup .with(state.into(), almanac) .for_duration(30.0 * Unit::Day) @@ -1196,7 +1203,11 @@ fn val_cislunar_dynamics(almanac_gmat: Arc) { ); let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(vec![EARTH, SUN, MOON])); - let setup = Propagator::new::(dynamics, PropOpts::with_fixed_step_s(0.5)); + let setup = Propagator::new( + dynamics, + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step_s(0.5), + ); let mut prop = setup.with(state.into(), almanac); prop.for_duration(prop_time).unwrap(); diff --git a/tests/mission_design/targeter/b_plane.rs b/tests/mission_design/targeter/b_plane.rs index c1889c02..51435b57 100644 --- a/tests/mission_design/targeter/b_plane.rs +++ b/tests/mission_design/targeter/b_plane.rs @@ -4,8 +4,8 @@ use anise::constants::celestial_objects::JUPITER_BARYCENTER; use anise::constants::celestial_objects::MOON; use anise::constants::celestial_objects::SUN; use anise::constants::frames::MOON_J2000; -use nyx::md::optimizer::*; use nyx::md::prelude::*; +use nyx::md::targeter::*; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; use rstest::*; @@ -48,7 +48,7 @@ fn tgt_b_plane_earth_gravity_assist_no_propagation(almanac: Arc) { let b_plane_tgt = BPlaneTarget::from_bt_br(13135.7982982557, 5022.26511510685); - let tgt = Optimizer::delta_v(&prop, b_plane_tgt.to_objectives()); + let tgt = Targeter::delta_v(&prop, b_plane_tgt.to_objectives()); let sol = tgt .try_achieve_from(spacecraft, epoch, epoch, almanac.clone()) @@ -117,7 +117,7 @@ fn tgt_b_plane_lunar_transfer(almanac: Arc) { // GMAT truth with central differencing: 1.15740867962, -0.576350387399, 0.632247251449 // GMAT truth with forward differencing: 1.33490412071, -0.5447988683, 1.77697094604 (approach currently in Nyx) - let tgt = Optimizer::in_frame( + let tgt = Targeter::in_frame( &prop, [ Variable { @@ -209,7 +209,7 @@ fn tgt_b_plane_earth_gravity_assist_with_propagation(almanac: Arc) { let b_plane_tgt = BPlaneTarget::from_bt_br(13135.7982982557, 5022.26511510685); - let tgt = Optimizer::delta_v(&prop, b_plane_tgt.to_objectives()); + let tgt = Targeter::delta_v(&prop, b_plane_tgt.to_objectives()); let sol = tgt .try_achieve_from(prior_sc, prior_sc.epoch(), epoch, almanac.clone()) diff --git a/tests/mission_design/targeter/finite_burns.rs b/tests/mission_design/targeter/finite_burns.rs index 104099fa..8ef6f91f 100644 --- a/tests/mission_design/targeter/finite_burns.rs +++ b/tests/mission_design/targeter/finite_burns.rs @@ -4,8 +4,8 @@ use anise::constants::celestial_objects::{JUPITER_BARYCENTER, MOON, SUN}; use hifitime::TimeUnits; use nyx::dynamics::guidance::{LocalFrame, Mnvr, Thruster}; use nyx::linalg::Vector3; -use nyx::md::optimizer::*; use nyx::md::prelude::*; +use nyx::md::targeter::*; use crate::propagation::GMAT_EARTH_GM; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; @@ -52,7 +52,7 @@ fn thrust_dir_tgt_sma_aop_raan(almanac: Arc) { Objective::within_tolerance(StateParameter::RAAN, 60.000182, 1e-3), ]; - let tgt = Optimizer::thrust_dir(&setup, objectives); + let tgt = Targeter::thrust_dir(&setup, objectives); println!("{}", tgt); @@ -99,7 +99,7 @@ fn thrust_dir_rate_tgt_sma_aop_raan(almanac: Arc) { Objective::within_tolerance(StateParameter::RAAN, 60.000182, 1e-3), ]; - let tgt = Optimizer::thrust_dir_rate(&setup, objectives); + let tgt = Targeter::thrust_dir_rate(&setup, objectives); println!("{}", tgt); @@ -148,7 +148,7 @@ fn thrust_profile_tgt_sma_aop_raan(almanac: Arc) { Objective::within_tolerance(StateParameter::RAAN, 60.000182, 1e-3), ]; - let tgt = Optimizer::thrust_profile(&setup, objectives); + let tgt = Targeter::thrust_profile(&setup, objectives); println!("{}", tgt); @@ -222,7 +222,7 @@ fn val_tgt_finite_burn(almanac: Arc) { let sc_no_thrust = SpacecraftDynamics::new(orbital_dyn); let mut prop_no_thrust = Propagator::default(sc_no_thrust); prop_no_thrust.set_max_step(mnvr0.duration()); - let impulsive_tgt = Optimizer::delta_v( + let impulsive_tgt = Targeter::delta_v( &prop_no_thrust, [ Objective::within_tolerance(StateParameter::X, sc_xf_desired.orbit.radius_km.x, 1e-5), diff --git a/tests/mission_design/targeter/multi_oe.rs b/tests/mission_design/targeter/multi_oe.rs index 77caa13c..ac5b11fd 100644 --- a/tests/mission_design/targeter/multi_oe.rs +++ b/tests/mission_design/targeter/multi_oe.rs @@ -2,8 +2,8 @@ extern crate nyx_space as nyx; // use nyx::dynamics::guidance::Mnvr; use nyx::dynamics::guidance::Thruster; -use nyx::md::optimizer::*; use nyx::md::prelude::*; +use nyx::md::targeter::*; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; use rstest::*; @@ -40,7 +40,7 @@ fn tgt_c3_decl(almanac: Arc) { Objective::within_tolerance(StateParameter::C3, -5.0, 0.5), ]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); @@ -106,7 +106,7 @@ fn conv_tgt_sma_ecc(almanac: Arc) { Objective::within_tolerance(StateParameter::SMA, 8100.0, 0.1), ]; - let tgt = Optimizer::new( + let tgt = Targeter::new( &setup, [ Variable { @@ -209,7 +209,7 @@ fn tgt_hd_sma_ecc(almanac: Arc) { Objective::within_tolerance(StateParameter::SMA, 8100.0, 0.1), ]; - let tgt = Optimizer::new( + let tgt = Targeter::new( &setup, [ Variable { diff --git a/tests/mission_design/targeter/multi_oe_vnc.rs b/tests/mission_design/targeter/multi_oe_vnc.rs index d265320a..7397b602 100644 --- a/tests/mission_design/targeter/multi_oe_vnc.rs +++ b/tests/mission_design/targeter/multi_oe_vnc.rs @@ -1,7 +1,7 @@ extern crate nyx_space as nyx; -use nyx::md::optimizer::*; use nyx::md::prelude::*; +use nyx::md::targeter::*; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; use rstest::*; @@ -38,7 +38,7 @@ fn tgt_vnc_c3_decl(almanac: Arc) { Objective::within_tolerance(StateParameter::C3, -5.0, 0.5), ]; - let tgt = Optimizer::vnc(&setup, objectives); + let tgt = Targeter::vnc(&setup, objectives); println!("{}", tgt); @@ -86,7 +86,7 @@ fn tgt_vnc_sma_ecc(almanac: Arc) { Objective::within_tolerance(StateParameter::SMA, 8100.0, 0.1), ]; - let tgt = Optimizer::vnc_with_components( + let tgt = Targeter::vnc_with_components( &setup, [ Variable { diff --git a/tests/mission_design/targeter/single_oe.rs b/tests/mission_design/targeter/single_oe.rs index b7736655..33603a63 100644 --- a/tests/mission_design/targeter/single_oe.rs +++ b/tests/mission_design/targeter/single_oe.rs @@ -3,8 +3,8 @@ extern crate nyx_space as nyx; use anise::constants::celestial_objects::JUPITER_BARYCENTER; use anise::constants::celestial_objects::MOON; use anise::constants::celestial_objects::SUN; -use nyx::md::optimizer::*; use nyx::md::prelude::*; +use nyx::md::targeter::*; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; use rstest::*; @@ -43,7 +43,7 @@ fn tgt_sma_from_apo(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::SMA, xf_desired_sma)]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); @@ -103,7 +103,7 @@ fn tgt_sma_from_peri_fd(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::SMA, xf_desired_sma)]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); @@ -161,7 +161,7 @@ fn tgt_hd_sma_from_peri(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::SMA, xf_desired_sma)]; - let mut tgt = Optimizer::delta_v(&setup, objectives); + let mut tgt = Targeter::delta_v(&setup, objectives); tgt.iterations = 5; println!("{}", tgt); @@ -268,7 +268,7 @@ fn tgt_ecc_from_apo(almanac: Arc) { let xf_desired_ecc = 0.4; - let tgt = Optimizer::new( + let tgt = Targeter::new( &setup, [ Variable { @@ -342,7 +342,7 @@ fn tgt_ecc_from_peri(almanac: Arc) { let xf_desired_ecc = 0.4; - let tgt = Optimizer::new( + let tgt = Targeter::new( &setup, [ Variable { @@ -416,7 +416,7 @@ fn tgt_raan_from_apo(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::RAAN, xf_desired_raan)]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); @@ -469,7 +469,7 @@ fn tgt_raan_from_peri(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::RAAN, xf_desired_raan)]; - let tgt = Optimizer::new( + let tgt = Targeter::new( &setup, [ Variable { @@ -543,7 +543,7 @@ fn tgt_aop_from_apo(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::AoP, xf_desired_aop)]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); @@ -596,7 +596,7 @@ fn tgt_aop_from_peri(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::AoP, xf_desired_aop)]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); diff --git a/tests/orbit_determination/measurements.rs b/tests/orbit_determination/measurements.rs index 48f528c0..0317ddd6 100644 --- a/tests/orbit_determination/measurements.rs +++ b/tests/orbit_determination/measurements.rs @@ -8,6 +8,7 @@ use nyx::dynamics::SpacecraftDynamics; use nyx::od::prelude::*; use nyx::time::Epoch; use nyx::{dynamics::OrbitalDynamics, propagators::Propagator}; +use nyx_space::propagators::IntegratorMethod; use rand::SeedableRng; use rand_pcg::Pcg64Mcg; @@ -78,7 +79,6 @@ fn val_measurements_topo(almanac: Arc) { use self::nyx::cosmic::Orbit; use self::nyx::md::prelude::*; use self::nyx::od::prelude::*; - use self::nyx::propagators::RK4Fixed; use std::str::FromStr; let cislunar1 = Orbit::cartesian( @@ -119,10 +119,11 @@ fn val_measurements_topo(almanac: Arc) { // Define the propagator information. let prop_time = 12 * Unit::Hour; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); - let setup = Propagator::new::( + let setup = Propagator::new( SpacecraftDynamics::new(OrbitalDynamics::point_masses(vec![EARTH, MOON, SUN])), + IntegratorMethod::RungeKutta4, opts, ); diff --git a/tests/orbit_determination/multi_body.rs b/tests/orbit_determination/multi_body.rs index 5d7b109d..97a5ba65 100644 --- a/tests/orbit_determination/multi_body.rs +++ b/tests/orbit_determination/multi_body.rs @@ -6,12 +6,12 @@ use anise::constants::celestial_objects::SUN; use anise::constants::frames::IAU_EARTH_FRAME; use nyx::od::simulator::TrackingArcSim; use nyx::od::simulator::TrkConfig; +use nyx_space::propagators::IntegratorMethod; use self::nyx::md::prelude::*; use self::nyx::od::prelude::*; // Extra testing imports -use self::nyx::propagators::RK4Fixed; use nyx::linalg::{SMatrix, SVector}; use std::collections::BTreeMap; @@ -101,7 +101,7 @@ fn od_val_multi_body_ckf_perfect_stations( // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -111,7 +111,11 @@ fn od_val_multi_body_ckf_perfect_stations( let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let orbital_dyn = OrbitalDynamics::point_masses(bodies); // Generate the truth data. - let setup = Propagator::new::(SpacecraftDynamics::new(orbital_dyn), opts); + let setup = Propagator::new( + SpacecraftDynamics::new(orbital_dyn), + IntegratorMethod::RungeKutta4, + opts, + ); let mut prop = setup.with(initial_state.into(), almanac.clone()); let (final_truth, traj) = prop.for_duration_with_traj(prop_time).unwrap(); @@ -145,7 +149,7 @@ fn od_val_multi_body_ckf_perfect_stations( let ckf = KF::no_snc(initial_estimate); - let mut odp = ODProcess::<_, _, RangeDoppler, _, _, _>::ckf(prop_est, ckf, None, almanac); + let mut odp = ODProcess::<_, RangeDoppler, _, _, _>::ckf(prop_est, ckf, None, almanac); odp.process_arc::(&arc).unwrap(); @@ -223,7 +227,7 @@ fn multi_body_ckf_covar_map( // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -233,7 +237,11 @@ fn multi_body_ckf_covar_map( // Generate the truth data on one thread. let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let orbital_dyn = OrbitalDynamics::point_masses(bodies); - let setup = Propagator::new::(SpacecraftDynamics::new(orbital_dyn), opts); + let setup = Propagator::new( + SpacecraftDynamics::new(orbital_dyn), + IntegratorMethod::RungeKutta4, + opts, + ); let mut prop = setup.with(initial_state.into(), almanac.clone()); let (_, traj) = prop.for_duration_with_traj(prop_time).unwrap(); diff --git a/tests/orbit_determination/resid_reject.rs b/tests/orbit_determination/resid_reject.rs index ba7c8288..b6a4b076 100644 --- a/tests/orbit_determination/resid_reject.rs +++ b/tests/orbit_determination/resid_reject.rs @@ -11,7 +11,7 @@ use nyx_space::cosmic::Orbit; use nyx_space::dynamics::orbital::OrbitalDynamics; use nyx_space::md::prelude::*; use nyx_space::od::prelude::*; -use nyx_space::propagators::{PropOpts, Propagator, RK4Fixed}; +use nyx_space::propagators::{IntegratorMethod, IntegratorOptions, Propagator}; use nyx_space::time::{Epoch, TimeUnits}; use nyx_space::utils::rss_orbit_errors; use std::collections::BTreeMap; @@ -50,7 +50,7 @@ fn traj(epoch: Epoch, almanac: Arc) -> Traj { let orbital_dyn = OrbitalDynamics::point_masses(bodies); let truth_setup = Propagator::dp78( SpacecraftDynamics::new(orbital_dyn), - PropOpts::with_max_step(step_size), + IntegratorOptions::with_max_step(step_size), ); let (_, traj) = truth_setup .with(initial_state, almanac) @@ -188,7 +188,11 @@ fn od_resid_reject_inflated_snc_ckf_two_way( let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let estimator = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::new::(estimator, PropOpts::with_fixed_step(10.seconds())); + let setup = Propagator::new( + estimator, + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.seconds()), + ); let prop_est = setup.with(initial_state_dev.with_stm(), almanac.clone()); // Define the process noise to assume an unmodeled acceleration on X, Y and Z in the ECI frame @@ -302,7 +306,11 @@ fn od_resid_reject_default_ckf_two_way( // We expect the estimated orbit to be _nearly_ perfect because we've removed SATURN_BARYCENTER from the estimated trajectory let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let estimator = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::new::(estimator, PropOpts::with_fixed_step(10.seconds())); + let setup = Propagator::new( + estimator, + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.seconds()), + ); let prop_est = setup.with(initial_state_dev.with_stm(), almanac.clone()); // Define the process noise to assume an unmodeled acceleration on X, Y and Z in the ECI frame diff --git a/tests/orbit_determination/robust.rs b/tests/orbit_determination/robust.rs index 7cfcdb88..119ff775 100644 --- a/tests/orbit_determination/robust.rs +++ b/tests/orbit_determination/robust.rs @@ -9,11 +9,12 @@ use nyx::dynamics::SpacecraftDynamics; use nyx::io::ExportCfg; use nyx::md::StateParameter; use nyx::od::prelude::*; -use nyx::propagators::{PropOpts, Propagator, RK4Fixed}; +use nyx::propagators::{IntegratorOptions, Propagator}; use nyx::time::{Epoch, TimeUnits, Unit}; use nyx::utils::rss_orbit_errors; use nyx::Spacecraft; use nyx_space::mc::StateDispersion; +use nyx_space::propagators::IntegratorMethod; use polars::prelude::*; use std::collections::BTreeMap; use std::env; @@ -79,7 +80,7 @@ fn od_robust_test_ekf_realistic_one_way(almanac: Arc) { // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -119,7 +120,11 @@ fn od_robust_test_ekf_realistic_one_way(almanac: Arc) { let bodies = vec![MOON, SUN, JUPITER_BARYCENTER, SATURN_BARYCENTER]; let orbital_dyn = OrbitalDynamics::point_masses(bodies); - let truth_setup = Propagator::new::(SpacecraftDynamics::new(orbital_dyn), opts); + let truth_setup = Propagator::new( + SpacecraftDynamics::new(orbital_dyn), + IntegratorMethod::RungeKutta4, + opts, + ); let (_, traj) = truth_setup .with(initial_state, almanac.clone()) .for_duration_with_traj(prop_time) @@ -146,7 +151,7 @@ fn od_robust_test_ekf_realistic_one_way(almanac: Arc) { // We expect the estimated orbit to be _nearly_ perfect because we"ve removed SATURN_BARYCENTER from the estimated trajectory let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let estimator = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::new::(estimator, opts); + let setup = Propagator::new(estimator, IntegratorMethod::RungeKutta4, opts); let prop_est = setup.with(initial_state_dev.with_stm(), almanac.clone()); // Define the process noise to assume an unmodeled acceleration on X, Y and Z in the ECI frame diff --git a/tests/orbit_determination/spacecraft.rs b/tests/orbit_determination/spacecraft.rs index 512d2db8..e4b97313 100644 --- a/tests/orbit_determination/spacecraft.rs +++ b/tests/orbit_determination/spacecraft.rs @@ -10,11 +10,11 @@ use nyx::linalg::{SMatrix, SVector}; use nyx::md::trajectory::ExportCfg; use nyx::md::{Event, StateParameter}; use nyx::od::prelude::*; -use nyx::propagators::{PropOpts, Propagator, RK4Fixed}; +use nyx::propagators::{IntegratorOptions, Propagator}; use nyx::time::{Epoch, TimeUnits, Unit}; use nyx_space::cosmic::SrpConfig; use nyx_space::dynamics::guidance::LocalFrame; -use nyx_space::propagators::RSSCartesianStep; +use nyx_space::propagators::{ErrorControl, IntegratorMethod}; use std::collections::BTreeMap; use std::path::PathBuf; @@ -124,7 +124,7 @@ fn od_val_sc_mb_srp_reals_duals_models( // Define the propagator information. let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -142,7 +142,7 @@ fn od_val_sc_mb_srp_reals_duals_models( let sc_init_state = Spacecraft::from_srp_defaults(initial_state, dry_mass_kg, sc_area); - let setup = Propagator::new::(sc_dynamics, opts); + let setup = Propagator::new(sc_dynamics, IntegratorMethod::RungeKutta4, opts); let mut prop = setup.with(sc_init_state, almanac.clone()); let (final_truth, traj) = prop.for_duration_with_traj(prop_time).unwrap(); @@ -328,9 +328,9 @@ fn od_val_sc_srp_estimation( let setup = Propagator::dp78( sc_dynamics, - PropOpts::builder() + IntegratorOptions::builder() .max_step(1.minutes()) - .error_ctrl(RSSCartesianStep {}) + .error_ctrl(ErrorControl::RSSCartesianStep) .build(), ); let mut prop = setup.with(sc_truth, almanac.clone()); diff --git a/tests/orbit_determination/two_body.rs b/tests/orbit_determination/two_body.rs index 9ee2c56b..fed2091c 100644 --- a/tests/orbit_determination/two_body.rs +++ b/tests/orbit_determination/two_body.rs @@ -10,8 +10,9 @@ use nyx::io::ConfigRepr; use nyx::io::{gravity::*, ExportCfg}; use nyx::linalg::{SMatrix, SVector}; use nyx::od::prelude::*; -use nyx::propagators::{PropOpts, Propagator, RK4Fixed}; +use nyx::propagators::{IntegratorOptions, Propagator}; use nyx::Spacecraft; +use nyx_space::propagators::IntegratorMethod; use std::collections::BTreeMap; use std::env; use std::path::PathBuf; @@ -116,7 +117,7 @@ fn od_tb_val_ekf_fixed_step_perfect_stations( // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -125,7 +126,7 @@ fn od_tb_val_ekf_fixed_step_perfect_stations( // We're sharing both the propagator and the dynamics. let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let mut prop = setup.with(initial_state.into(), almanac.clone()); let (final_truth, traj) = prop.for_duration_with_traj(prop_time).unwrap(); @@ -245,7 +246,11 @@ fn od_tb_val_with_arc( // We're sharing both the propagator and the dynamics. let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - let setup = Propagator::new::(orbital_dyn, PropOpts::with_fixed_step_s(10.0)); + let setup = Propagator::new( + orbital_dyn, + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step_s(10.0), + ); let mut prop = setup.with(initial_state.into(), almanac.clone()); let (final_truth, traj) = prop.for_duration_with_traj(duration).unwrap(); @@ -408,7 +413,7 @@ fn od_tb_val_ckf_fixed_step_perfect_stations( // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -418,7 +423,7 @@ fn od_tb_val_ckf_fixed_step_perfect_stations( // Generate the truth data on one thread. let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let mut prop = setup.with(initial_state.into(), almanac.clone()); let (final_truth, traj) = prop.for_duration_with_traj(prop_time).unwrap(); @@ -619,7 +624,7 @@ fn od_tb_ckf_fixed_step_iteration_test( // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -627,7 +632,7 @@ fn od_tb_ckf_fixed_step_iteration_test( let initial_state = Orbit::keplerian(22000.0, 0.01, 30.0, 80.0, 40.0, 0.0, dt, eme2k); let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let mut prop = setup.with(initial_state.into(), almanac.clone()); let (final_truth, traj) = prop.for_duration_with_traj(prop_time).unwrap(); @@ -756,7 +761,7 @@ fn od_tb_ckf_fixed_step_perfect_stations_snc_covar_map( // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -764,7 +769,7 @@ fn od_tb_ckf_fixed_step_perfect_stations_snc_covar_map( let initial_state = Orbit::keplerian(22000.0, 0.01, 30.0, 80.0, 40.0, 0.0, dt, eme2k); let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let mut prop = setup.with(initial_state.into(), almanac.clone()); let (final_truth, traj) = prop.for_duration_with_traj(prop_time).unwrap(); @@ -873,9 +878,10 @@ fn od_tb_ckf_map_covar(almanac: Arc) { // We expect the estimated orbit to be perfect since we're using strictly the same dynamics, no noise on // the measurements, and the same time step. - let setup = Propagator::new::( + let setup = Propagator::new( SpacecraftDynamics::new(OrbitalDynamics::two_body()), - PropOpts::with_fixed_step(step_size), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(step_size), ); let prop_est = setup.with(Spacecraft::from(initial_state).with_stm(), almanac.clone()); let covar_radius_km = 1.0e-3; @@ -948,7 +954,7 @@ fn od_tb_val_harmonics_ckf_fixed_step_perfect( // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -959,7 +965,7 @@ fn od_tb_val_harmonics_ckf_fixed_step_perfect( let earth_sph_harm = HarmonicsMem::from_cof("data/JGM3.cof.gz", 70, 70, true).unwrap(); let harmonics = Harmonics::from_stor(iau_earth, earth_sph_harm); let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::from_model(harmonics)); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let mut prop = setup.with(initial_state.into(), almanac.clone()); let (final_truth, traj) = prop.for_duration_with_traj(prop_time).unwrap(); @@ -1056,7 +1062,7 @@ fn od_tb_ckf_fixed_step_perfect_stations_several_snc_covar_map( // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -1064,7 +1070,7 @@ fn od_tb_ckf_fixed_step_perfect_stations_several_snc_covar_map( let initial_state = Orbit::keplerian(22000.0, 0.01, 30.0, 80.0, 40.0, 0.0, dt, eme2k); let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let mut prop = setup.with(initial_state.into(), almanac.clone()); let (final_truth, traj) = prop.for_duration_with_traj(prop_time).unwrap(); diff --git a/tests/orbit_determination/xhat_dev.rs b/tests/orbit_determination/xhat_dev.rs index 9a69c06f..cb8ad9a9 100644 --- a/tests/orbit_determination/xhat_dev.rs +++ b/tests/orbit_determination/xhat_dev.rs @@ -10,9 +10,10 @@ use nyx::dynamics::SpacecraftDynamics; use nyx::io::gravity::*; use nyx::linalg::{SMatrix, SVector}; use nyx::od::prelude::*; -use nyx::propagators::{PropOpts, Propagator, RK4Fixed}; +use nyx::propagators::{IntegratorOptions, Propagator}; use nyx::utils::rss_orbit_errors; use nyx::Spacecraft; +use nyx_space::propagators::IntegratorMethod; use std::collections::BTreeMap; use std::convert::TryFrom; @@ -70,7 +71,7 @@ fn xhat_dev_test_ekf_two_body(almanac: Arc) { // Define the propagator information. let prop_time = 0.01 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -90,7 +91,7 @@ fn xhat_dev_test_ekf_two_body(almanac: Arc) { ); let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let (_, traj) = setup .with(Spacecraft::from(initial_state), almanac.clone()) .for_duration_with_traj(prop_time) @@ -298,7 +299,7 @@ fn xhat_dev_test_ekf_multi_body(almanac: Arc) { // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -319,7 +320,7 @@ fn xhat_dev_test_ekf_multi_body(almanac: Arc) { let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let (_, traj) = setup .with(initial_state.into(), almanac.clone()) @@ -468,7 +469,7 @@ fn xhat_dev_test_ekf_harmonics(almanac: Arc) { // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -499,7 +500,7 @@ fn xhat_dev_test_ekf_harmonics(almanac: Arc) { PointMasses::new(bodies), ])); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let (_, traj) = setup .with(initial_state.into(), almanac.clone()) @@ -631,7 +632,7 @@ fn xhat_dev_test_ekf_realistic(almanac: Arc) { // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -649,7 +650,7 @@ fn xhat_dev_test_ekf_realistic(almanac: Arc) { let bodies = vec![MOON, SUN, JUPITER_BARYCENTER, SATURN_BARYCENTER]; let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let truth_setup = Propagator::new::(orbital_dyn, opts); + let truth_setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let (_, traj) = truth_setup .with(initial_state.into(), almanac.clone()) @@ -666,7 +667,7 @@ fn xhat_dev_test_ekf_realistic(almanac: Arc) { // We expect the estimated orbit to be _nearly_ perfect because we've removed SATURN_BARYCENTER from the estimated trajectory let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let estimator = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::new::(estimator, opts); + let setup = Propagator::new(estimator, IntegratorMethod::RungeKutta4, opts); let prop_est = setup.with(Spacecraft::from(initial_state).with_stm(), almanac.clone()); let covar_radius_km = 1.0e2; let covar_velocity_km_s = 1.0e1; @@ -783,7 +784,7 @@ fn xhat_dev_test_ckf_smoother_multi_body(almanac: Arc) { // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -804,7 +805,7 @@ fn xhat_dev_test_ckf_smoother_multi_body(almanac: Arc) { let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let (_, traj) = setup .with(initial_state.into(), almanac.clone()) .for_duration_with_traj(prop_time) @@ -1053,7 +1054,7 @@ fn xhat_dev_test_ekf_snc_smoother_multi_body(almanac: Arc) { // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -1074,7 +1075,7 @@ fn xhat_dev_test_ekf_snc_smoother_multi_body(almanac: Arc) { let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let (_, traj) = setup .with(initial_state.into(), almanac.clone()) .for_duration_with_traj(prop_time) @@ -1324,7 +1325,7 @@ fn xhat_dev_test_ckf_iteration_multi_body(almanac: Arc) { // Define the propagator information. let prop_time = 1 * Unit::Day; let step_size = 10.0 * Unit::Second; - let opts = PropOpts::with_fixed_step(step_size); + let opts = IntegratorOptions::with_fixed_step(step_size); // Define state information. let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); @@ -1345,7 +1346,7 @@ fn xhat_dev_test_ckf_iteration_multi_body(almanac: Arc) { let bodies = vec![MOON, SUN, JUPITER_BARYCENTER]; let orbital_dyn = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies)); - let setup = Propagator::new::(orbital_dyn, opts); + let setup = Propagator::new(orbital_dyn, IntegratorMethod::RungeKutta4, opts); let (_, traj) = setup .with(initial_state.into(), almanac.clone()) .for_duration_with_traj(prop_time) diff --git a/tests/propagation/events.rs b/tests/propagation/events.rs index b48f0bad..eb8837ff 100644 --- a/tests/propagation/events.rs +++ b/tests/propagation/events.rs @@ -38,7 +38,7 @@ fn event_tracker_true_anomaly(almanac: Arc) { let events = vec![peri_event, apo_event, ta_event0, ta_event1]; let dynamics = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - let setup = Propagator::rk89(dynamics, PropOpts::with_tolerance(1e-9)); + let setup = Propagator::rk89(dynamics, IntegratorOptions::with_tolerance(1e-9)); let mut prop = setup.with(state.into(), almanac.clone()); let (_, traj) = prop.for_duration_with_traj(prop_time).unwrap(); diff --git a/tests/propagation/propagators.rs b/tests/propagation/propagators.rs index 958f45d8..4d97694a 100644 --- a/tests/propagation/propagators.rs +++ b/tests/propagation/propagators.rs @@ -5,7 +5,7 @@ use hifitime::JD_J2000; use nyx::cosmic::{assert_orbit_eq_or_abs, Orbit}; use nyx::dynamics::orbital::OrbitalDynamics; use nyx::dynamics::SpacecraftDynamics; -use nyx::propagators::error_ctrl::RSSCartesianState; +use nyx::propagators::error_ctrl::ErrorControl; use nyx::time::{Epoch, Unit}; use nyx::utils::rss_orbit_errors; use nyx::{propagators::*, Spacecraft}; @@ -40,107 +40,44 @@ fn regress_leo_day_adaptive(almanac: Arc) { )); let final_dt = dt + prop_time; - let all_rslts = vec![ - Orbit::cartesian( - -5_971.198_524_908_157, - 3_945.775_509_326_305_4, - 2_864.262_542_023_422, - 0.048_766_212_879_869_19, - -4.184_873_956_982_518, - 5.849_098_380_963_502, - final_dt, - eme2k, - ), - Orbit::cartesian( - -5_971.194_190_197_366, - 3_945.506_606_221_459_6, - 2_864.636_682_800_498_4, - 0.049_097_015_227_526_38, - -4.185_093_356_859_808, - 5.848_940_840_578_1, - final_dt, - eme2k, - ), - Orbit::cartesian( - -5_971.194_190_305_766, - 3_945.506_612_356_549_3, - 2_864.636_674_277_756_4, - 0.049_097_007_640_393_29, - -4.185_093_351_832_897, - 5.848_940_844_198_66, - final_dt, - eme2k, - ), - ]; + let rslt = Orbit::cartesian( + -5_971.194_190_197_366, + 3_945.506_606_221_459_6, + 2_864.636_682_800_498_4, + 0.049_097_015_227_526_38, + -4.185_093_356_859_808, + 5.848_940_840_578_1, + final_dt, + eme2k, + ); let dynamics = SpacecraftDynamics::new(OrbitalDynamics::two_body()); - { - let setup = Propagator::new::( - dynamics.clone(), - PropOpts::with_fixed_step(1.0 * Unit::Second), - ); - let mut prop = setup.with(init, almanac.clone()); - prop.for_duration(prop_time).unwrap(); + let setup = Propagator::new( + dynamics.clone(), + IntegratorMethod::CashKarp45, + IntegratorOptions::with_adaptive_step( + min_step, + max_step, + accuracy, + ErrorControl::RSSCartesianState, + ), + ); + let mut prop = setup.with(init, almanac.clone()); + prop.for_duration(prop_time).unwrap(); + assert_orbit_eq_or_abs( + &prop.state.orbit, + &rslt, + 1e-7, + "CashKarp45 two body prop failed", + ); + let prev_details = prop.latest_details(); + if prev_details.error > accuracy { assert_eq!( - prop.state.orbit.to_cartesian_pos_vel(), - all_rslts[0].to_cartesian_pos_vel(), - "RK2Fixed two body prop failed" - ); - let prev_details = prop.latest_details(); - if prev_details.error > accuracy { - assert_eq!( - prev_details.step, min_step, - "step size should be at its minimum because error is higher than tolerance: {:?}", - prev_details - ); - } - } - - { - let setup = Propagator::new::( - dynamics.clone(), - PropOpts::with_adaptive_step(min_step, max_step, accuracy, RSSCartesianState {}), - ); - let mut prop = setup.with(init, almanac.clone()); - prop.for_duration(prop_time).unwrap(); - assert_orbit_eq_or_abs( - &prop.state.orbit, - &all_rslts[1], - 1e-7, - "CashKarp45 two body prop failed", + prev_details.step, min_step, + "step size should be at its minimum because error is higher than tolerance: {:?}", + prev_details ); - let prev_details = prop.latest_details(); - if prev_details.error > accuracy { - assert_eq!( - prev_details.step, min_step, - "step size should be at its minimum because error is higher than tolerance: {:?}", - prev_details - ); - } - } - - { - let setup = Propagator::new::( - dynamics, - PropOpts::with_adaptive_step(min_step, max_step, accuracy, RSSCartesianState {}), - ); - let mut prop = setup.with(init, almanac.clone()); - prop.for_duration(prop_time).unwrap(); - assert_orbit_eq_or_abs( - &prop.state.orbit, - &all_rslts[2], - 1e-7, - "Fehlberg45 two body prop failed", - ); - let prev_details = prop.latest_details(); - if prev_details.error > accuracy { - assert_eq!( - prev_details.step, min_step, - "step size should be at its minimum because error is higher than tolerance: {:?}", - prev_details - ); - } } } @@ -211,9 +148,15 @@ fn gmat_val_leo_day_adaptive(almanac: Arc) { let dynamics = SpacecraftDynamics::new(OrbitalDynamics::two_body()); { - let setup = Propagator::new::( + let setup = Propagator::new( dynamics.clone(), - PropOpts::with_adaptive_step(min_step, max_step, accuracy, RSSCartesianState {}), + IntegratorMethod::DormandPrince45, + IntegratorOptions::with_adaptive_step( + min_step, + max_step, + accuracy, + ErrorControl::RSSCartesianState, + ), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -262,9 +205,15 @@ fn gmat_val_leo_day_adaptive(almanac: Arc) { } { - let setup = Propagator::new::( + let setup = Propagator::new( dynamics.clone(), - PropOpts::with_adaptive_step(min_step, max_step, accuracy, RSSCartesianState {}), + IntegratorMethod::Verner56, + IntegratorOptions::with_adaptive_step( + min_step, + max_step, + accuracy, + ErrorControl::RSSCartesianState, + ), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -291,9 +240,14 @@ fn gmat_val_leo_day_adaptive(almanac: Arc) { } { - let setup = Propagator::new::( + let setup = Propagator::dp78( dynamics.clone(), - PropOpts::with_adaptive_step(min_step, max_step, accuracy, RSSCartesianState {}), + IntegratorOptions::with_adaptive_step( + min_step, + max_step, + accuracy, + ErrorControl::RSSCartesianState, + ), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -319,9 +273,14 @@ fn gmat_val_leo_day_adaptive(almanac: Arc) { } { - let setup = Propagator::new::( + let setup = Propagator::rk89( dynamics, - PropOpts::with_adaptive_step(min_step, max_step, accuracy, RSSCartesianState {}), + IntegratorOptions::with_adaptive_step( + min_step, + max_step, + accuracy, + ErrorControl::RSSCartesianState, + ), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -418,9 +377,10 @@ fn gmat_val_leo_day_fixed(almanac: Arc) { let dynamics = SpacecraftDynamics::new(OrbitalDynamics::two_body()); { - let setup = Propagator::new::( + let setup = Propagator::new( dynamics.clone(), - PropOpts::with_fixed_step(1.0 * Unit::Second), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(1.0 * Unit::Second), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -444,9 +404,10 @@ fn gmat_val_leo_day_fixed(almanac: Arc) { } { - let setup = Propagator::new::( + let setup = Propagator::new( dynamics.clone(), - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorMethod::Verner56, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -465,9 +426,10 @@ fn gmat_val_leo_day_fixed(almanac: Arc) { } { - let setup = Propagator::new::( + let setup = Propagator::new( dynamics.clone(), - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorMethod::DormandPrince45, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -481,9 +443,10 @@ fn gmat_val_leo_day_fixed(almanac: Arc) { } { - let setup = Propagator::new::( + let setup = Propagator::new( dynamics.clone(), - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorMethod::DormandPrince78, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -497,8 +460,10 @@ fn gmat_val_leo_day_fixed(almanac: Arc) { } { - let setup = - Propagator::new::(dynamics, PropOpts::with_fixed_step(10.0 * Unit::Second)); + let setup = Propagator::rk89( + dynamics, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); assert_eq!(prop.state.orbit, all_rslts[4], "two body prop failed"); diff --git a/tests/propagation/stm.rs b/tests/propagation/stm.rs index 2cf88f41..32d8d296 100644 --- a/tests/propagation/stm.rs +++ b/tests/propagation/stm.rs @@ -32,9 +32,10 @@ fn stm_fixed_step(almanac: Arc) { .with_mu_km3_s2(GMAT_EARTH_GM); let epoch = Epoch::from_gregorian_tai_at_midnight(2020, 1, 1); - let prop = Propagator::new::( + let prop = Propagator::new( SpacecraftDynamics::new(OrbitalDynamics::two_body()), - PropOpts::with_fixed_step(10 * Unit::Second), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10 * Unit::Second), ); // First test is in mostly linear regime (low eccentricity) diff --git a/tests/propagation/stopcond.rs b/tests/propagation/stopcond.rs index 60955310..a3442517 100644 --- a/tests/propagation/stopcond.rs +++ b/tests/propagation/stopcond.rs @@ -14,11 +14,11 @@ use nyx::dynamics::guidance::{FiniteBurns, LocalFrame, Mnvr, Thruster}; use nyx::dynamics::orbital::OrbitalDynamics; use nyx::dynamics::SpacecraftDynamics; use nyx::md::{Event, StateParameter}; -use nyx::propagators::error_ctrl::RSSCartesianStep; -use nyx::propagators::{PropOpts, Propagator}; +use nyx::propagators::{IntegratorOptions, Propagator}; use nyx::time::{Epoch, TimeUnits, Unit}; use nyx::{Spacecraft, State}; +use nyx_space::propagators::ErrorControl; use rstest::*; #[fixture] @@ -163,7 +163,7 @@ fn stop_cond_nrho_apo(almanac: Arc) { let setup = Propagator::rk89( dynamics, - PropOpts::with_adaptive_step_s(1.0, 60.0, 1e-6, RSSCartesianStep {}), + IntegratorOptions::with_adaptive_step_s(1.0, 60.0, 1e-6, ErrorControl::RSSCartesianStep), ); // NOTE: Here, we will propagate for the maximum duration in the original frame diff --git a/tests/propulsion/closedloop_multi_oe_ruggiero.rs b/tests/propulsion/closedloop_multi_oe_ruggiero.rs index feae05a8..042dc8df 100644 --- a/tests/propulsion/closedloop_multi_oe_ruggiero.rs +++ b/tests/propulsion/closedloop_multi_oe_ruggiero.rs @@ -7,12 +7,13 @@ use self::nyx::cosmic::{GuidanceMode, Orbit, Spacecraft}; use self::nyx::dynamics::guidance::{Objective, Ruggiero, Thruster}; use self::nyx::dynamics::{OrbitalDynamics, SpacecraftDynamics}; use self::nyx::md::{Event, StateParameter}; -use self::nyx::propagators::{PropOpts, Propagator, RK4Fixed}; +use self::nyx::propagators::{IntegratorOptions, Propagator}; use self::nyx::time::{Epoch, Unit}; /// NOTE: Herein shows the difference between the QLaw and Ruggiero (and other control laws). /// The Ruggiero control law takes quite some longer to converge than the QLaw. use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; +use nyx_space::propagators::IntegratorMethod; use rstest::*; #[fixture] @@ -66,8 +67,11 @@ fn qlaw_as_ruggiero_case_a(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, ruggiero_ctrl); println!("[qlaw_as_ruggiero_case_a] {:x}", orbit); - let setup = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)); + let setup = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ); let mut prop = setup.with(sc_state, almanac.clone()); let (final_state, traj) = prop.for_duration_with_traj(prop_time).unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; @@ -128,11 +132,14 @@ fn qlaw_as_ruggiero_case_b(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, ruggiero_ctrl); println!("[qlaw_as_ruggiero_case_b] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[qlaw_as_ruggiero_case_b] {:x}", final_state.orbit); @@ -183,11 +190,14 @@ fn qlaw_as_ruggiero_case_c(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, ruggiero_ctrl); println!("[qlaw_as_ruggiero_case_c] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[qlaw_as_ruggiero_case_c] {:x}", final_state.orbit); @@ -241,11 +251,14 @@ fn qlaw_as_ruggiero_case_d(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, ruggiero_ctrl); println!("[qlaw_as_ruggiero_case_d] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[qlaw_as_ruggiero_case_d] {:x}", final_state.orbit); @@ -301,11 +314,14 @@ fn qlaw_as_ruggiero_case_e(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, ruggiero_ctrl); println!("[qlaw_as_ruggiero_case_e] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[qlaw_as_ruggiero_case_e] {:x}", final_state.orbit); @@ -357,8 +373,11 @@ fn qlaw_as_ruggiero_case_f(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, ruggiero_ctrl); println!("[qlaw_as_ruggiero_case_f] {:x}", orbit); - let setup = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)); + let setup = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ); let (final_state, traj) = setup .with(sc_state, almanac.clone()) .for_duration_with_traj(prop_time) @@ -417,11 +436,14 @@ fn ruggiero_iepc_2011_102(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, ruggiero_ctrl); println!("[ruggiero_iepc_2011_102] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[ruggiero_iepc_2011_102] {:x}", final_state.orbit); diff --git a/tests/propulsion/closedloop_single_oe_ruggiero.rs b/tests/propulsion/closedloop_single_oe_ruggiero.rs index 6116c134..ca3b3bf3 100644 --- a/tests/propulsion/closedloop_single_oe_ruggiero.rs +++ b/tests/propulsion/closedloop_single_oe_ruggiero.rs @@ -4,10 +4,11 @@ extern crate nyx_space as nyx; use self::nyx::cosmic::{GuidanceMode, Orbit, Spacecraft}; use self::nyx::dynamics::guidance::{Objective, Ruggiero, StateParameter, Thruster}; use self::nyx::dynamics::{OrbitalDynamics, SpacecraftDynamics}; -use self::nyx::propagators::{PropOpts, Propagator, RK4Fixed}; +use self::nyx::propagators::{IntegratorOptions, Propagator}; use self::nyx::time::{Epoch, Unit}; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; +use nyx_space::propagators::IntegratorMethod; use rstest::*; use std::sync::Arc; @@ -53,11 +54,14 @@ fn rugg_sma(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, guid_law); println!("[rugg_sma] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_sma] {:x}", final_state.orbit); @@ -103,11 +107,14 @@ fn rugg_sma_regress_threshold(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(OrbitalDynamics::two_body(), guid_law); println!("[rugg_sma_regress] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac.clone()) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac.clone()) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_sma_regress] {:x}", final_state.orbit); @@ -157,11 +164,14 @@ fn rugg_sma_decr(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, guid_law); println!("[rugg_sma_decr] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_sma_decr] {:x}", final_state.orbit); @@ -212,11 +222,14 @@ fn rugg_inc(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, guid_law); println!("[rugg_inc] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_inc] {:x}", final_state.orbit); @@ -268,11 +281,14 @@ fn rugg_inc_threshold(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, guid_law); println!("[rugg_inc_threshold] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_inc_threshold] {:x}", final_state.orbit); @@ -323,11 +339,14 @@ fn rugg_inc_decr(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, guid_law); println!("[rugg_inc_decr] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_inc_decr] {:x}", final_state.orbit); @@ -378,11 +397,14 @@ fn rugg_ecc(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, guid_law); println!("[rugg_ecc] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_ecc] {:x}", final_state.orbit); @@ -432,11 +454,14 @@ fn rugg_ecc_regress_threshold(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(OrbitalDynamics::two_body(), guid_law); println!("[rugg_ecc_regress] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac.clone()) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac.clone()) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_ecc_regress] {:x}", final_state.orbit); @@ -488,11 +513,14 @@ fn rugg_ecc_decr(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, guid_law); println!("[rugg_ecc_decr] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_ecc_decr] {:x}", final_state.orbit); @@ -545,11 +573,14 @@ fn rugg_aop(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, guid_law); println!("[rugg_aop] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_aop] {:x}", final_state.orbit); @@ -601,11 +632,14 @@ fn rugg_aop_decr(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, guid_law); println!("[rugg_aop_decr] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_aop_decr] {:x}", final_state.orbit); @@ -654,8 +688,11 @@ fn rugg_raan(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, guid_law); println!("[rugg_raan] {:x}", orbit); - let setup = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)); + let setup = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ); let mut prop = setup.with(sc_state, almanac.clone()); let (final_state, traj) = prop.for_duration_with_traj(prop_time).unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; @@ -704,11 +741,14 @@ fn rugg_raan_regress_threshold(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(OrbitalDynamics::two_body(), guid_law); println!("[rugg_raan_regress] {:x}", orbit); - let final_state = - Propagator::new::(sc.clone(), PropOpts::with_fixed_step(10.0 * Unit::Second)) - .with(sc_state, almanac.clone()) - .for_duration(prop_time) - .unwrap(); + let final_state = Propagator::new( + sc.clone(), + IntegratorMethod::RungeKutta4, + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), + ) + .with(sc_state, almanac.clone()) + .for_duration(prop_time) + .unwrap(); let fuel_usage = fuel_mass - final_state.fuel_mass_kg; println!("[rugg_raan_regress] {:x}", final_state.orbit); diff --git a/tests/propulsion/schedule.rs b/tests/propulsion/schedule.rs index 2f1b4a91..5e2aa1d8 100644 --- a/tests/propulsion/schedule.rs +++ b/tests/propulsion/schedule.rs @@ -3,7 +3,7 @@ use self::nyx::cosmic::{GuidanceMode, Orbit, Spacecraft}; use self::nyx::dynamics::guidance::{FiniteBurns, Mnvr, Thruster}; use self::nyx::dynamics::{OrbitalDynamics, SpacecraftDynamics}; use self::nyx::linalg::Vector3; -use self::nyx::propagators::{PropOpts, Propagator}; +use self::nyx::propagators::{IntegratorOptions, Propagator}; use self::nyx::time::{Epoch, Unit}; use self::nyx::utils::rss_orbit_vec_errors; use crate::propagation::GMAT_EARTH_GM; @@ -72,7 +72,7 @@ fn val_transfer_schedule_no_depl(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law_no_decr(orbital_dyn, schedule); // Setup a propagator, and propagate for that duration // NOTE: We specify the use an RK89 to match the GMAT setup. - let final_state = Propagator::rk89(sc, PropOpts::with_fixed_step(10.0 * Unit::Second)) + let final_state = Propagator::rk89(sc, IntegratorOptions::with_fixed_step(10.0 * Unit::Second)) .with(sc_state, almanac) .for_duration(prop_time) .unwrap(); @@ -163,7 +163,7 @@ fn val_transfer_schedule_depl(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, schedule); // Setup a propagator, and propagate for that duration // NOTE: We specify the use an RK89 to match the GMAT setup. - let setup = Propagator::rk89(sc, PropOpts::with_fixed_step(10.0 * Unit::Second)); + let setup = Propagator::rk89(sc, IntegratorOptions::with_fixed_step(10.0 * Unit::Second)); let final_state = setup .with(sc_state, almanac.clone()) .for_duration(prop_time) @@ -294,7 +294,7 @@ fn val_transfer_single_maneuver_depl(almanac: Arc) { let sc = SpacecraftDynamics::from_guidance_law(orbital_dyn, Arc::new(mnvr0)); // Setup a propagator, and propagate for that duration // NOTE: We specify the use an RK89 to match the GMAT setup. - let setup = Propagator::rk89(sc, PropOpts::with_fixed_step(10.0 * Unit::Second)); + let setup = Propagator::rk89(sc, IntegratorOptions::with_fixed_step(10.0 * Unit::Second)); let final_state = setup .with(sc_state, almanac.clone()) .for_duration(prop_time)