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/mc/montecarlo.rs b/src/mc/montecarlo.rs index f582d8f7..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( + pub fn run_until_nth_event( self, - prop: Propagator, + 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( + pub fn resume_run_until_nth_event( &self, - prop: Propagator, + 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( + pub fn run_until_epoch( self, - prop: Propagator, + 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( + pub fn resume_run_until_epoch( &self, - prop: Propagator, + 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..ec7a6377 100644 --- a/src/md/mod.rs +++ b/src/md/mod.rs @@ -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}; diff --git a/src/md/opti/multipleshooting/altitude_heuristic.rs b/src/md/opti/multipleshooting/altitude_heuristic.rs index 768b2543..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, + 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 6f6abbda..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, + 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 8c92b907..579bf3e7 100644 --- a/src/md/opti/multipleshooting/multishoot.rs +++ b/src/md/opti/multipleshooting/multishoot.rs @@ -24,7 +24,6 @@ use crate::linalg::{DMatrix, DVector, SVector}; use crate::md::opti::solution::TargeterSolution; use crate::md::optimizer::Optimizer; 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, + 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, @@ -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,9 +348,9 @@ 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( + pub fn build_trajectories( &self, - prop: &Propagator, + prop: &Propagator, almanac: Arc, ) -> Result, MultipleShootingError> { let mut trajz = Vec::with_capacity(self.nodes.len()); diff --git a/src/md/opti/optimizer.rs b/src/md/opti/optimizer.rs index 856747d1..dda79528 100644 --- a/src/md/opti/optimizer.rs +++ b/src/md/opti/optimizer.rs @@ -26,16 +26,17 @@ 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; +// TODO(now): rename to Differential Controller + /// 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 Optimizer<'a, const V: usize, const O: usize> { /// The propagator setup (kind, stages, etc.) - pub prop: &'a Propagator, + 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 +50,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 Optimizer<'a, V, O> { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { let mut objmsg = String::from(""); for obj in &self.objectives { @@ -65,12 +66,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> Optimizer<'a, 3, O> { /// Create a new Targeter which will apply an impulsive delta-v correction. - pub fn delta_v( - prop: &'a Propagator, - objectives: [Objective; O], - ) -> Self { + pub fn delta_v(prop: &'a Propagator, objectives: [Objective; O]) -> Self { Self { prop, objectives, @@ -86,10 +84,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, - objectives: [Objective; O], - ) -> Self { + pub fn delta_r(prop: &'a Propagator, objectives: [Objective; O]) -> Self { Self { prop, objectives, @@ -105,7 +100,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, objectives: [Objective; O]) -> Self { + pub fn vnc(prop: &'a Propagator, objectives: [Objective; O]) -> Self { Self { prop, objectives, @@ -121,10 +116,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> Optimizer<'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, + prop: &'a Propagator, objectives: [Objective; O], ) -> Self { Self { @@ -143,10 +138,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> Optimizer<'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, + prop: &'a Propagator, objectives: [Objective; O], ) -> Self { Self { @@ -168,10 +163,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> Optimizer<'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, + prop: &'a Propagator, objectives: [Objective; O], ) -> Self { Self { @@ -196,10 +191,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> Optimizer<'a, V, O> { /// Create a new Targeter which will apply an impulsive delta-v correction. pub fn new( - prop: &'a Propagator, + prop: &'a Propagator, variables: [Variable; V], objectives: [Objective; O], ) -> Self { @@ -215,7 +210,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, + prop: &'a Propagator, variables: [Variable; V], objectives: [Objective; O], objective_frame: Frame, @@ -232,7 +227,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, + prop: &'a Propagator, variables: [Variable; V], objectives: [Objective; O], ) -> Self { diff --git a/src/md/opti/raphson_finite_diff.rs b/src/md/opti/raphson_finite_diff.rs index 88b9112a..5d3f3104 100644 --- a/src/md/opti/raphson_finite_diff.rs +++ b/src/md/opti/raphson_finite_diff.rs @@ -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> Optimizer<'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..7d4ffa30 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> Optimizer<'a, V, O> { /// Differential correction using hyperdual numbers for the objectives #[allow(clippy::comparison_chain)] pub fn try_achieve_dual( 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..0fc2efde 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,80 +84,104 @@ 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 - } -} - -/// 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::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 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 + } + } + 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 + } + } } } } @@ -145,8 +193,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 +224,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 +243,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 c1be7125..7214a69d 100644 --- a/src/propagators/instance.rs +++ b/src/propagators/instance.rs @@ -16,7 +16,6 @@ along with this program. If not, see . */ -use super::error_ctrl::ErrorCtrl; use super::{DynamicsSnafu, IntegrationDetails, PropagationError, Propagator}; use crate::dynamics::Dynamics; use crate::linalg::allocator::Allocator; @@ -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, + 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> @@ -388,7 +387,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 diff --git a/src/propagators/options.rs b/src/propagators/options.rs index 1c52e1ab..2d8e59c8 100644 --- a/src/propagators/options.rs +++ b/src/propagators/options.rs @@ -20,7 +20,7 @@ use std::fmt; use crate::time::{Duration, Unit}; -use super::{ErrorCtrl, RSSCartesianStep}; +use super::ErrorControl; use typed_builder::TypedBuilder; /// PropOpts stores the integrator options, including the minimum and maximum step sizes, and the @@ -32,7 +32,7 @@ use typed_builder::TypedBuilder; /// fixed step options will lead to an RK4 being used instead of an RK45. #[derive(Clone, Copy, Debug, TypedBuilder)] #[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 +45,20 @@ pub struct PropOpts { pub attempts: u8, #[builder(default = false)] pub fixed_step: bool, - pub error_ctrl: E, + pub error_ctrl: ErrorControl, + // TODO(now): Add an optional central body that will swap central bodies on init } -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, @@ -72,7 +73,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 +83,40 @@ 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, + } + } + + 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 +139,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,74 +153,36 @@ 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, } } } #[test] fn test_options() { - use super::error_ctrl::RSSStep; - - let opts = PropOpts::with_fixed_step_s(1e-1); + 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 = PropOpts::with_adaptive_step_s(1e-2, 10.0, 1e-12, RSSStep {}); + 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: PropOpts = Default::default(); + 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); @@ -193,7 +190,7 @@ fn test_options() { assert_eq!(opts.attempts, 50); assert!(!opts.fixed_step); - let opts = PropOpts::with_max_step(1.0 * Unit::Second); + 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); diff --git a/src/propagators/propagator.rs b/src/propagators/propagator.rs index 831a06a8..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::{IntegrationDetails, IntegratorMethod, PropInstance, PropOpts}; +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 +pub struct Propagator where DefaultAllocator: Allocator<::Size> + Allocator<::Size, ::Size> @@ -40,12 +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 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 Propagator +impl Propagator where DefaultAllocator: Allocator<::Size> + Allocator<::Size, ::Size> @@ -53,7 +52,7 @@ where + Allocator<::VecLength>, { /// Each propagator must be initialized with `new` which stores propagator information. - pub fn new(dynamics: D, method: IntegratorMethod, opts: PropOpts) -> Self { + pub fn new(dynamics: D, method: IntegratorMethod, opts: IntegratorOptions) -> Self { Self { dynamics, opts, @@ -76,17 +75,17 @@ where } /// An RK89 propagator (the default) with custom propagator options. - pub fn rk89(dynamics: D, opts: PropOpts) -> Self { + 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 { + pub fn dp78(dynamics: D, opts: IntegratorOptions) -> Self { Self::new(dynamics, IntegratorMethod::DormandPrince78, opts) } - pub fn with(&self, state: D::StateType, almanac: Arc) -> PropInstance { + 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.method.stages() + 1); for _ in 0..self.method.stages() { @@ -107,24 +106,16 @@ where k, } } -} -impl Propagator -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::rk89(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::dp78(dynamics, PropOpts::default()) + Self::dp78(dynamics, IntegratorOptions::default()) } } diff --git a/src/propagators/rk_methods/mod.rs b/src/propagators/rk_methods/mod.rs index 15fc0c8e..dea49f02 100644 --- a/src/propagators/rk_methods/mod.rs +++ b/src/propagators/rk_methods/mod.rs @@ -19,6 +19,9 @@ mod rk; use std::str::FromStr; +use serde::Deserialize; +use serde::Serialize; + use crate::io::ConfigError; pub use self::rk::*; @@ -31,7 +34,7 @@ use super::PropagationError; /// The `RK` trait defines a Runge Kutta integrator. #[allow(clippy::upper_case_acronyms)] -pub trait RK +trait RK where Self: Sized, { @@ -54,7 +57,7 @@ where /// 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)] +#[derive(Copy, Clone, Debug, PartialEq, Eq, Serialize, Deserialize)] pub enum IntegratorMethod { /// Runge Kutta 8-9 is the recommended integrator for most application. RungeKutta89, diff --git a/src/propagators/rk_methods/rk.rs b/src/propagators/rk_methods/rk.rs index 731f8746..fe78c203 100644 --- a/src/propagators/rk_methods/rk.rs +++ b/src/propagators/rk_methods/rk.rs @@ -16,7 +16,7 @@ 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 {} 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/orbitaldyn.rs b/tests/mission_design/orbitaldyn.rs index 1ad7e7af..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; @@ -56,7 +55,7 @@ fn energy_conservation(almanac: Arc) { let rk89_final = Propagator::new( SpacecraftDynamics::new(OrbitalDynamics::two_body()), IntegratorMethod::RungeKutta89, - PropOpts::default(), + IntegratorOptions::default(), ) .with(Spacecraft::from(start_state), almanac.clone()) .for_duration(prop_time) @@ -78,7 +77,7 @@ fn energy_conservation(almanac: Arc) { let dp78_final = Propagator::new( SpacecraftDynamics::new(OrbitalDynamics::two_body()), IntegratorMethod::DormandPrince78, - PropOpts::default(), + IntegratorOptions::default(), ) .with(start_state.into(), almanac) .for_duration(prop_time) @@ -122,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(); @@ -211,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); @@ -275,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); @@ -341,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); @@ -403,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); @@ -754,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(); @@ -807,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) @@ -917,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(); @@ -970,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()) @@ -994,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()) @@ -1157,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) @@ -1201,7 +1206,7 @@ fn val_cislunar_dynamics(almanac_gmat: Arc) { let setup = Propagator::new( dynamics, IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step_s(0.5), + 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/orbit_determination/measurements.rs b/tests/orbit_determination/measurements.rs index 5ecdd84b..0317ddd6 100644 --- a/tests/orbit_determination/measurements.rs +++ b/tests/orbit_determination/measurements.rs @@ -119,7 +119,7 @@ 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( SpacecraftDynamics::new(OrbitalDynamics::point_masses(vec![EARTH, MOON, SUN])), diff --git a/tests/orbit_determination/multi_body.rs b/tests/orbit_determination/multi_body.rs index b594145d..97a5ba65 100644 --- a/tests/orbit_determination/multi_body.rs +++ b/tests/orbit_determination/multi_body.rs @@ -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(); @@ -149,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(); @@ -227,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(); diff --git a/tests/orbit_determination/resid_reject.rs b/tests/orbit_determination/resid_reject.rs index bb0cc99d..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::{IntegratorMethod, PropOpts, Propagator}; +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) @@ -191,7 +191,7 @@ fn od_resid_reject_inflated_snc_ckf_two_way( let setup = Propagator::new( estimator, IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.seconds()), + IntegratorOptions::with_fixed_step(10.seconds()), ); let prop_est = setup.with(initial_state_dev.with_stm(), almanac.clone()); @@ -309,7 +309,7 @@ fn od_resid_reject_default_ckf_two_way( let setup = Propagator::new( estimator, IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.seconds()), + IntegratorOptions::with_fixed_step(10.seconds()), ); let prop_est = setup.with(initial_state_dev.with_stm(), almanac.clone()); diff --git a/tests/orbit_determination/robust.rs b/tests/orbit_determination/robust.rs index 7bbf515d..119ff775 100644 --- a/tests/orbit_determination/robust.rs +++ b/tests/orbit_determination/robust.rs @@ -9,7 +9,7 @@ use nyx::dynamics::SpacecraftDynamics; use nyx::io::ExportCfg; use nyx::md::StateParameter; use nyx::od::prelude::*; -use nyx::propagators::{PropOpts, Propagator}; +use nyx::propagators::{IntegratorOptions, Propagator}; use nyx::time::{Epoch, TimeUnits, Unit}; use nyx::utils::rss_orbit_errors; use nyx::Spacecraft; @@ -80,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(); diff --git a/tests/orbit_determination/spacecraft.rs b/tests/orbit_determination/spacecraft.rs index 3305b6e1..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}; +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::{IntegratorMethod, 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(); @@ -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 4ee51f1f..fed2091c 100644 --- a/tests/orbit_determination/two_body.rs +++ b/tests/orbit_determination/two_body.rs @@ -10,7 +10,7 @@ use nyx::io::ConfigRepr; use nyx::io::{gravity::*, ExportCfg}; use nyx::linalg::{SMatrix, SVector}; use nyx::od::prelude::*; -use nyx::propagators::{PropOpts, Propagator}; +use nyx::propagators::{IntegratorOptions, Propagator}; use nyx::Spacecraft; use nyx_space::propagators::IntegratorMethod; use std::collections::BTreeMap; @@ -117,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(); @@ -249,7 +249,7 @@ fn od_tb_val_with_arc( let setup = Propagator::new( orbital_dyn, IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step_s(10.0), + IntegratorOptions::with_fixed_step_s(10.0), ); let mut prop = setup.with(initial_state.into(), almanac.clone()); @@ -413,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(); @@ -624,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(); @@ -761,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(); @@ -881,7 +881,7 @@ fn od_tb_ckf_map_covar(almanac: Arc) { let setup = Propagator::new( SpacecraftDynamics::new(OrbitalDynamics::two_body()), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(step_size), + 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; @@ -954,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(); @@ -1062,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(); diff --git a/tests/orbit_determination/xhat_dev.rs b/tests/orbit_determination/xhat_dev.rs index 56d37e04..cb8ad9a9 100644 --- a/tests/orbit_determination/xhat_dev.rs +++ b/tests/orbit_determination/xhat_dev.rs @@ -10,7 +10,7 @@ use nyx::dynamics::SpacecraftDynamics; use nyx::io::gravity::*; use nyx::linalg::{SMatrix, SVector}; use nyx::od::prelude::*; -use nyx::propagators::{PropOpts, Propagator}; +use nyx::propagators::{IntegratorOptions, Propagator}; use nyx::utils::rss_orbit_errors; use nyx::Spacecraft; use nyx_space::propagators::IntegratorMethod; @@ -71,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(); @@ -299,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(); @@ -469,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(); @@ -632,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(); @@ -784,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(); @@ -1054,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(); @@ -1325,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(); 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 c7fe4e83..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}; @@ -56,7 +56,12 @@ fn regress_leo_day_adaptive(almanac: Arc) { let setup = Propagator::new( dynamics.clone(), IntegratorMethod::CashKarp45, - 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(); @@ -146,7 +151,12 @@ fn gmat_val_leo_day_adaptive(almanac: Arc) { let setup = Propagator::new( dynamics.clone(), IntegratorMethod::DormandPrince45, - 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(); @@ -198,7 +208,12 @@ fn gmat_val_leo_day_adaptive(almanac: Arc) { let setup = Propagator::new( dynamics.clone(), IntegratorMethod::Verner56, - 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(); @@ -227,7 +242,12 @@ fn gmat_val_leo_day_adaptive(almanac: Arc) { { 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(); @@ -255,7 +275,12 @@ fn gmat_val_leo_day_adaptive(almanac: Arc) { { 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(); @@ -355,7 +380,7 @@ fn gmat_val_leo_day_fixed(almanac: Arc) { let setup = Propagator::new( dynamics.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(1.0 * Unit::Second), + IntegratorOptions::with_fixed_step(1.0 * Unit::Second), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -382,7 +407,7 @@ fn gmat_val_leo_day_fixed(almanac: Arc) { let setup = Propagator::new( dynamics.clone(), IntegratorMethod::Verner56, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -404,7 +429,7 @@ fn gmat_val_leo_day_fixed(almanac: Arc) { let setup = Propagator::new( dynamics.clone(), IntegratorMethod::DormandPrince45, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -421,7 +446,7 @@ fn gmat_val_leo_day_fixed(almanac: Arc) { let setup = Propagator::new( dynamics.clone(), IntegratorMethod::DormandPrince78, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ); let mut prop = setup.with(init, almanac.clone()); prop.for_duration(prop_time).unwrap(); @@ -435,7 +460,10 @@ fn gmat_val_leo_day_fixed(almanac: Arc) { } { - let setup = Propagator::rk89(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 848a0b66..32d8d296 100644 --- a/tests/propagation/stm.rs +++ b/tests/propagation/stm.rs @@ -35,7 +35,7 @@ fn stm_fixed_step(almanac: Arc) { let prop = Propagator::new( SpacecraftDynamics::new(OrbitalDynamics::two_body()), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10 * Unit::Second), + 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 19c55e5b..042dc8df 100644 --- a/tests/propulsion/closedloop_multi_oe_ruggiero.rs +++ b/tests/propulsion/closedloop_multi_oe_ruggiero.rs @@ -7,7 +7,7 @@ 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}; +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). @@ -70,7 +70,7 @@ fn qlaw_as_ruggiero_case_a(almanac: Arc) { let setup = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + 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(); @@ -135,7 +135,7 @@ fn qlaw_as_ruggiero_case_b(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -193,7 +193,7 @@ fn qlaw_as_ruggiero_case_c(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -254,7 +254,7 @@ fn qlaw_as_ruggiero_case_d(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -317,7 +317,7 @@ fn qlaw_as_ruggiero_case_e(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -376,7 +376,7 @@ fn qlaw_as_ruggiero_case_f(almanac: Arc) { let setup = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ); let (final_state, traj) = setup .with(sc_state, almanac.clone()) @@ -439,7 +439,7 @@ fn ruggiero_iepc_2011_102(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) diff --git a/tests/propulsion/closedloop_single_oe_ruggiero.rs b/tests/propulsion/closedloop_single_oe_ruggiero.rs index c4da0f73..ca3b3bf3 100644 --- a/tests/propulsion/closedloop_single_oe_ruggiero.rs +++ b/tests/propulsion/closedloop_single_oe_ruggiero.rs @@ -4,7 +4,7 @@ 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}; +use self::nyx::propagators::{IntegratorOptions, Propagator}; use self::nyx::time::{Epoch, Unit}; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; @@ -57,7 +57,7 @@ fn rugg_sma(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -110,7 +110,7 @@ fn rugg_sma_regress_threshold(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac.clone()) .for_duration(prop_time) @@ -167,7 +167,7 @@ fn rugg_sma_decr(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -225,7 +225,7 @@ fn rugg_inc(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -284,7 +284,7 @@ fn rugg_inc_threshold(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -342,7 +342,7 @@ fn rugg_inc_decr(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -400,7 +400,7 @@ fn rugg_ecc(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -457,7 +457,7 @@ fn rugg_ecc_regress_threshold(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac.clone()) .for_duration(prop_time) @@ -516,7 +516,7 @@ fn rugg_ecc_decr(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -576,7 +576,7 @@ fn rugg_aop(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -635,7 +635,7 @@ fn rugg_aop_decr(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac) .for_duration(prop_time) @@ -691,7 +691,7 @@ fn rugg_raan(almanac: Arc) { let setup = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + 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(); @@ -744,7 +744,7 @@ fn rugg_raan_regress_threshold(almanac: Arc) { let final_state = Propagator::new( sc.clone(), IntegratorMethod::RungeKutta4, - PropOpts::with_fixed_step(10.0 * Unit::Second), + IntegratorOptions::with_fixed_step(10.0 * Unit::Second), ) .with(sc_state, almanac.clone()) .for_duration(prop_time) 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)