use std::sync::{Arc, Mutex};
use nalgebra::{DMatrix, DVector, Vector3, Vector6};
use crate::constants::{GM_EARTH, R_EARTH};
use crate::earth_models::{density_harris_priester, density_nrlmsise00};
use crate::frames::{earth_rotation, rotation_eci_to_ecef};
use crate::integrators::traits::DIntegrator;
use crate::math::SMatrix3;
use crate::math::interpolation::{
CovarianceInterpolationConfig, CovarianceInterpolationMethod, InterpolationConfig,
InterpolationMethod, interpolate_hermite_cubic_dvector6,
};
use crate::math::jacobian::DNumericalJacobian;
use crate::math::jacobian::DifferenceMethod;
use crate::math::sensitivity::DNumericalSensitivity;
use crate::orbit_dynamics::{
GravityModel, accel_drag, accel_gravity_spherical_harmonics_with_workspace,
accel_point_mass_gravity, accel_relativity, accel_solar_radiation_pressure, accel_third_body,
eclipse_conical, eclipse_cylindrical, get_global_gravity_model, sun_position,
};
use crate::propagators::{
AtmosphericModel, EclipseModel, ForceModelConfig, FrameTransformationModel,
GravityConfiguration, GravityModelSource,
};
use crate::relative_motion::rotation_eci_to_rtn;
use crate::time::Epoch;
use crate::traits::{OrbitFrame, OrbitRepresentation};
use crate::trajectories::DOrbitTrajectory;
use crate::trajectories::traits::{
InterpolatableTrajectory, STMStorage, SensitivityStorage, Trajectory,
};
use crate::utils::errors::BraheError;
use crate::utils::identifiable::Identifiable;
use crate::utils::state_providers::{
DCovarianceProvider, DOrbitCovarianceProvider, DOrbitStateProvider, DStateProvider,
};
use crate::{AngleFormat, accel_earth_zonal_gravity, state_eci_to_koe};
use super::TrajectoryMode;
use super::traits::DStatePropagator;
use crate::events::{DDetectedEvent, DEventDetector, EventAction, EventType, dscan_for_event};
use crate::integrators::traits::{DControlInput, DStateDynamics};
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[allow(dead_code)]
enum PropagationMode {
StateOnly,
WithSTM,
WithSensitivity,
WithSTMAndSensitivity,
}
#[derive(Debug)]
enum EventProcessingResult {
NoEvents,
Restart { epoch: Epoch, state: DVector<f64> },
Terminal { epoch: Epoch, state: DVector<f64> },
}
type SharedDynamics =
Arc<dyn Fn(f64, &DVector<f64>, Option<&DVector<f64>>) -> DVector<f64> + Send + Sync>;
const ROTATION_CACHE_CAPACITY: usize = 20;
struct RotationCache {
entries: [Option<(f64, SMatrix3)>; ROTATION_CACHE_CAPACITY],
next: usize,
model: FrameTransformationModel,
}
impl RotationCache {
fn new(model: FrameTransformationModel) -> Self {
Self {
entries: [None; ROTATION_CACHE_CAPACITY],
next: 0,
model,
}
}
fn get_or_compute(&mut self, t: f64, epoch: Epoch) -> SMatrix3 {
for entry in &self.entries {
if let Some((cached_t, r)) = entry
&& *cached_t == t
{
return *r;
}
}
let r = match self.model {
FrameTransformationModel::FullEarthRotation => rotation_eci_to_ecef(epoch),
FrameTransformationModel::EarthRotationOnly => earth_rotation(epoch),
};
self.entries[self.next] = Some((t, r));
self.next = (self.next + 1) % ROTATION_CACHE_CAPACITY;
r
}
}
struct DynamicsWorkspace {
rotation_cache: RotationCache,
sh_v: DMatrix<f64>,
sh_w: DMatrix<f64>,
}
impl DynamicsWorkspace {
fn new(frame_transform: FrameTransformationModel) -> Self {
Self {
rotation_cache: RotationCache::new(frame_transform),
sh_v: DMatrix::<f64>::zeros(0, 0),
sh_w: DMatrix::<f64>::zeros(0, 0),
}
}
}
pub struct DNumericalOrbitPropagator {
epoch_initial: Epoch,
epoch_current: Epoch,
t_rel: f64,
integrator: Box<dyn DIntegrator>,
dt: f64,
dt_next: f64,
x_initial: DVector<f64>,
x_curr: DVector<f64>,
params: DVector<f64>,
state_dim: usize,
propagation_mode: PropagationMode,
stm: Option<DMatrix<f64>>,
sensitivity: Option<DMatrix<f64>>,
store_stm_history: bool,
store_sensitivity_history: bool,
initial_covariance: Option<DMatrix<f64>>,
current_covariance: Option<DMatrix<f64>>,
trajectory: DOrbitTrajectory,
trajectory_mode: TrajectoryMode,
interpolation_method: InterpolationMethod,
covariance_interpolation_method: CovarianceInterpolationMethod,
store_accelerations: bool,
shared_dynamics: SharedDynamics,
event_detectors: Vec<Box<dyn DEventDetector>>,
event_log: Vec<DDetectedEvent>,
terminated: bool,
pub name: Option<String>,
pub id: Option<u64>,
pub uuid: Option<uuid::Uuid>,
}
pub struct Set;
pub struct Unset;
struct BuilderData {
epoch: Option<Epoch>,
state: Option<DVector<f64>>,
force_config: Option<ForceModelConfig>,
propagation_config: Option<super::NumericalPropagationConfig>,
params: Option<DVector<f64>>,
additional_dynamics: Option<DStateDynamics>,
control_input: DControlInput,
initial_covariance: Option<DMatrix<f64>>,
}
pub struct DNumericalOrbitPropagatorBuilder<E, S, F> {
data: BuilderData,
_phantom: std::marker::PhantomData<(E, S, F)>,
}
impl<S, F> DNumericalOrbitPropagatorBuilder<Unset, S, F> {
pub fn epoch(mut self, epoch: Epoch) -> DNumericalOrbitPropagatorBuilder<Set, S, F> {
self.data.epoch = Some(epoch);
DNumericalOrbitPropagatorBuilder {
data: self.data,
_phantom: std::marker::PhantomData,
}
}
}
impl<E, F> DNumericalOrbitPropagatorBuilder<E, Unset, F> {
pub fn state(mut self, state: DVector<f64>) -> DNumericalOrbitPropagatorBuilder<E, Set, F> {
self.data.state = Some(state);
DNumericalOrbitPropagatorBuilder {
data: self.data,
_phantom: std::marker::PhantomData,
}
}
}
impl<E, S> DNumericalOrbitPropagatorBuilder<E, S, Unset> {
pub fn force_config(
mut self,
config: ForceModelConfig,
) -> DNumericalOrbitPropagatorBuilder<E, S, Set> {
self.data.force_config = Some(config);
DNumericalOrbitPropagatorBuilder {
data: self.data,
_phantom: std::marker::PhantomData,
}
}
}
impl<E, S, F> DNumericalOrbitPropagatorBuilder<E, S, F> {
pub fn propagation_config(mut self, config: super::NumericalPropagationConfig) -> Self {
self.data.propagation_config = Some(config);
self
}
pub fn params(mut self, params: DVector<f64>) -> Self {
self.data.params = Some(params);
self
}
pub fn additional_dynamics(mut self, dynamics: DStateDynamics) -> Self {
self.data.additional_dynamics = Some(dynamics);
self
}
pub fn control_input(
mut self,
control: Box<
dyn Fn(f64, &DVector<f64>, Option<&DVector<f64>>) -> DVector<f64> + Send + Sync,
>,
) -> Self {
self.data.control_input = Some(control);
self
}
pub fn initial_covariance(mut self, covariance: DMatrix<f64>) -> Self {
self.data.initial_covariance = Some(covariance);
self
}
}
impl DNumericalOrbitPropagatorBuilder<Set, Set, Set> {
pub fn build(self) -> Result<DNumericalOrbitPropagator, BraheError> {
DNumericalOrbitPropagator::new(
self.data.epoch.unwrap(),
self.data.state.unwrap(),
self.data.propagation_config.unwrap_or_default(),
self.data.force_config.unwrap(),
self.data.params,
self.data.additional_dynamics,
self.data.control_input,
self.data.initial_covariance,
)
}
}
impl DNumericalOrbitPropagator {
#[allow(clippy::too_many_arguments)]
pub fn new(
epoch: Epoch,
state: DVector<f64>,
propagation_config: super::NumericalPropagationConfig,
force_config: ForceModelConfig,
params: Option<DVector<f64>>,
additional_dynamics: Option<DStateDynamics>,
control_input: DControlInput,
initial_covariance: Option<DMatrix<f64>>,
) -> Result<Self, BraheError> {
propagation_config.validate()?;
force_config.validate_params(params.as_ref())?;
let params = params.unwrap_or_else(|| DVector::zeros(0));
let state_eci = state;
let state_dim = state_eci.len();
let enable_stm = propagation_config.variational.enable_stm || initial_covariance.is_some();
let enable_sensitivity = propagation_config.variational.enable_sensitivity;
if enable_sensitivity && params.is_empty() {
return Err(BraheError::PropagatorError(
"Sensitivity propagation requires params to be provided".to_string(),
));
}
let gravity_model = match &force_config.gravity {
GravityConfiguration::SphericalHarmonic {
source: GravityModelSource::ModelType(model_type),
degree,
order,
} => {
let mut arc = GravityModel::shared(model_type)?;
if *degree < arc.n_max || *order < arc.m_max {
Arc::make_mut(&mut arc).set_max_degree_order(*degree, *order)?;
}
Some(arc)
}
_ => None,
};
let shared_dynamics = Self::build_shared_dynamics(
epoch,
force_config.clone(),
params.clone(),
additional_dynamics,
gravity_model.clone(),
);
let dynamics = Self::wrap_for_integrator(Arc::clone(&shared_dynamics));
let initial_dt = propagation_config.integrator.initial_step.unwrap_or(60.0);
let jacobian_provider = if enable_stm || enable_sensitivity {
Some(Self::build_jacobian_provider(
Arc::clone(&shared_dynamics),
propagation_config.variational.jacobian_method,
))
} else {
None
};
let sensitivity_provider = if enable_sensitivity {
Some(Self::build_sensitivity_provider(
Arc::clone(&shared_dynamics),
propagation_config.variational.sensitivity_method,
))
} else {
None
};
let integrator = crate::integrators::create_dintegrator(
propagation_config.method,
state_dim,
dynamics,
jacobian_provider,
sensitivity_provider,
control_input,
propagation_config.integrator,
);
let mut trajectory = DOrbitTrajectory::new(
state_dim,
OrbitFrame::ECI,
OrbitRepresentation::Cartesian,
None,
);
trajectory.set_interpolation_method(propagation_config.interpolation_method);
if propagation_config.variational.store_stm_history {
trajectory.enable_stm_storage();
}
if propagation_config.variational.store_sensitivity_history && !params.is_empty() {
trajectory.enable_sensitivity_storage(params.len());
}
if propagation_config.store_accelerations {
trajectory.enable_acceleration_storage(3); }
let initial_stm = if propagation_config.variational.store_stm_history {
Some(DMatrix::identity(state_dim, state_dim))
} else {
None
};
let initial_sensitivity =
if propagation_config.variational.store_sensitivity_history && !params.is_empty() {
Some(DMatrix::zeros(state_dim, params.len()))
} else {
None
};
let initial_acceleration = if propagation_config.store_accelerations {
let dx = shared_dynamics(0.0, &state_eci, Some(¶ms));
Some(DVector::from_vec(vec![dx[3], dx[4], dx[5]]))
} else {
None
};
trajectory.add_full(
epoch,
state_eci.clone(),
initial_covariance.clone(),
initial_stm,
initial_sensitivity,
initial_acceleration,
);
let propagation_mode = match (enable_stm, enable_sensitivity) {
(false, false) => PropagationMode::StateOnly,
(true, false) => PropagationMode::WithSTM,
(false, true) => PropagationMode::WithSensitivity,
(true, true) => PropagationMode::WithSTMAndSensitivity,
};
let stm = if enable_stm {
Some(DMatrix::identity(state_dim, state_dim))
} else {
None
};
let sensitivity = if enable_sensitivity {
Some(DMatrix::zeros(state_dim, params.len()))
} else {
None
};
let current_covariance = initial_covariance.clone();
let mut result = Ok(Self {
epoch_initial: epoch,
epoch_current: epoch,
t_rel: 0.0,
integrator,
dt: initial_dt,
dt_next: initial_dt,
x_initial: state_eci.clone(),
x_curr: state_eci,
params,
state_dim,
propagation_mode,
stm,
sensitivity,
store_stm_history: propagation_config.variational.store_stm_history,
store_sensitivity_history: propagation_config.variational.store_sensitivity_history,
initial_covariance,
current_covariance,
trajectory,
trajectory_mode: TrajectoryMode::AllSteps,
interpolation_method: propagation_config.interpolation_method,
covariance_interpolation_method: CovarianceInterpolationMethod::TwoWasserstein,
store_accelerations: propagation_config.store_accelerations,
shared_dynamics,
event_detectors: Vec::new(),
event_log: Vec::new(),
terminated: false,
name: None,
id: None,
uuid: None,
});
if let Ok(ref mut prop) = result {
let uuid = uuid::Uuid::now_v7();
prop.uuid = Some(uuid);
prop.trajectory.uuid = Some(uuid);
}
result
}
pub fn from_gp_record(
record: &crate::types::GPRecord,
force_config: ForceModelConfig,
propagation_config: Option<super::NumericalPropagationConfig>,
params: Option<DVector<f64>>,
) -> Result<Self, BraheError> {
use crate::propagators::traits::SStatePropagator;
let sgp = crate::propagators::SGPPropagator::from_gp_record(record, 60.0)?;
let epoch = sgp.epoch;
let state = DVector::from_column_slice(sgp.initial_state().as_slice());
Self::new(
epoch,
state,
propagation_config.unwrap_or_default(),
force_config,
params,
None,
None,
None,
)
}
pub fn builder() -> DNumericalOrbitPropagatorBuilder<Unset, Unset, Unset> {
DNumericalOrbitPropagatorBuilder {
data: BuilderData {
epoch: None,
state: None,
force_config: None,
propagation_config: None,
params: None,
additional_dynamics: None,
control_input: None,
initial_covariance: None,
},
_phantom: std::marker::PhantomData,
}
}
fn compute_acceleration(&self, epoch: Epoch, state: &DVector<f64>) -> Option<DVector<f64>> {
if !self.store_accelerations {
return None;
}
let t_rel = epoch - self.epoch_initial;
let dx = (self.shared_dynamics)(t_rel, state, Some(&self.params));
Some(DVector::from_vec(vec![dx[3], dx[4], dx[5]]))
}
fn process_initial_events(&mut self) {
if self.event_detectors.is_empty() || self.terminated {
return;
}
let epoch = self.current_epoch();
let state = self.x_curr.clone();
let params_ref = if self.params.is_empty() {
None
} else {
Some(&self.params)
};
let mut events_to_process = Vec::new();
for (idx, detector) in self.event_detectors.iter().enumerate() {
if detector.is_processed() {
continue;
}
let value = detector.evaluate(epoch, &state, params_ref);
let target = detector.target_value();
if (value - target).abs() <= detector.value_tolerance() {
let event = DDetectedEvent {
window_open: epoch,
window_close: epoch,
entry_state: state.clone(),
exit_state: state.clone(),
value,
name: detector.name().to_string(),
action: detector.action(),
event_type: EventType::Instantaneous,
detector_index: idx,
};
events_to_process.push((idx, event));
}
}
if !events_to_process.is_empty() {
for &(idx, _) in &events_to_process {
self.event_detectors[idx].mark_processed();
}
match self.process_events_smart(events_to_process) {
EventProcessingResult::NoEvents => {}
EventProcessingResult::Restart { epoch, state } => {
self.t_rel = epoch - self.epoch_initial;
self.epoch_current = epoch;
self.x_curr = state;
self.integrator.reset_cache();
}
EventProcessingResult::Terminal {
epoch: term_epoch,
state: term_state,
} => {
self.epoch_current = term_epoch;
self.x_curr = term_state;
}
}
}
}
fn scan_all_events(
&self,
epoch_prev: Epoch,
epoch_new: Epoch,
x_prev: &DVector<f64>,
x_new: &DVector<f64>,
) -> Vec<(usize, DDetectedEvent)> {
let mut events = Vec::new();
let t_prev_rel = epoch_prev - self.epoch_initial;
let t_new_rel = epoch_new - self.epoch_initial;
let state_fn = |epoch: Epoch| -> DVector<f64> {
let t_rel = epoch - self.epoch_initial;
interpolate_hermite_cubic_dvector6(t_prev_rel, t_new_rel, x_prev, x_new, t_rel)
};
for (idx, detector) in self.event_detectors.iter().enumerate() {
if let Some(event) = dscan_for_event(
detector.as_ref(),
idx,
&state_fn,
epoch_prev,
epoch_new,
x_prev,
x_new,
Some(&self.params),
) {
events.push((idx, event));
}
}
events
}
fn refine_event_states(
&self,
events: &mut [(usize, DDetectedEvent)],
epoch_prev: Epoch,
x_prev: &DVector<f64>,
) {
if events.is_empty() {
return;
}
let t_prev_rel = epoch_prev - self.epoch_initial;
let params_ref = if self.params.is_empty() {
None
} else {
Some(&self.params)
};
self.integrator.reset_cache();
for (_, event) in events.iter_mut() {
let dt_to_event = event.window_open - epoch_prev;
if dt_to_event.abs() <= 1e-12 {
event.entry_state = x_prev.clone();
event.exit_state = x_prev.clone();
continue;
}
let max_sub_dt = 10.0;
let n_steps = ((dt_to_event.abs() / max_sub_dt).ceil() as usize).max(1);
let sub_dt = dt_to_event / n_steps as f64;
let mut t = t_prev_rel;
let mut x = x_prev.clone();
for _ in 0..n_steps {
let result = self.integrator.step(t, x, params_ref, Some(sub_dt));
x = result.state;
t += result.dt_used;
}
event.entry_state = x.clone();
event.exit_state = x;
}
self.integrator.reset_cache();
}
fn process_events_smart(
&mut self,
detected_events: Vec<(usize, DDetectedEvent)>,
) -> EventProcessingResult {
if detected_events.is_empty() {
return EventProcessingResult::NoEvents;
}
let mut sorted_events = detected_events;
sorted_events.sort_by(|(_, a), (_, b)| a.window_open.partial_cmp(&b.window_open).unwrap());
let first_callback_idx = sorted_events
.iter()
.position(|(det_idx, _)| self.event_detectors[*det_idx].callback().is_some());
match first_callback_idx {
None => {
let mut terminal_event = None;
for (det_idx, event) in sorted_events {
self.event_log.push(event.clone());
self.event_detectors[det_idx].mark_processed();
if !matches!(self.trajectory_mode, TrajectoryMode::Disabled) {
self.trajectory
.add(event.window_open, event.entry_state.clone());
}
let detector = &self.event_detectors[det_idx];
if detector.action() == EventAction::Stop {
terminal_event = Some(event);
break; }
}
if let Some(term_event) = terminal_event {
self.terminated = true;
return EventProcessingResult::Terminal {
epoch: term_event.window_open,
state: term_event.entry_state.clone(),
};
}
EventProcessingResult::NoEvents }
Some(callback_idx) => {
for (det_idx, event) in sorted_events.iter().take(callback_idx) {
self.event_log.push(event.clone());
self.event_detectors[*det_idx].mark_processed();
if !matches!(self.trajectory_mode, TrajectoryMode::Disabled) {
self.trajectory
.add(event.window_open, event.entry_state.clone());
}
}
let (det_idx, callback_event) = &sorted_events[callback_idx];
let detector = &self.event_detectors[*det_idx];
self.event_log.push(callback_event.clone());
if !matches!(self.trajectory_mode, TrajectoryMode::Disabled) {
let acc = self.compute_acceleration(
callback_event.window_open,
&callback_event.entry_state,
);
self.trajectory.add_full(
callback_event.window_open,
callback_event.entry_state.clone(),
None,
None,
None,
acc,
);
}
if let Some(callback) = detector.callback() {
let (new_state, new_params, action) = callback(
callback_event.window_open,
&callback_event.entry_state,
Some(&self.params),
);
let y_after = if let Some(y_new) = new_state {
if !matches!(self.trajectory_mode, TrajectoryMode::Disabled) {
let acc = self.compute_acceleration(callback_event.window_open, &y_new);
self.trajectory.add_full(
callback_event.window_open,
y_new.clone(),
None,
None,
None,
acc,
);
}
y_new
} else {
callback_event.entry_state.clone()
};
if let Some(p_new) = new_params {
self.params = p_new;
}
detector.mark_processed();
if action == EventAction::Stop {
self.terminated = true;
return EventProcessingResult::Terminal {
epoch: callback_event.window_open,
state: y_after,
};
}
return EventProcessingResult::Restart {
epoch: callback_event.window_open,
state: y_after,
};
}
unreachable!("Callback event must have callback");
}
}
}
fn build_shared_dynamics(
epoch_initial: Epoch,
force_config: ForceModelConfig,
params: DVector<f64>,
additional_dynamics: Option<DStateDynamics>,
gravity_model: Option<Arc<GravityModel>>,
) -> SharedDynamics {
let workspace = Arc::new(Mutex::new(DynamicsWorkspace::new(
force_config.frame_transform.clone(),
)));
Arc::new(
move |t: f64,
state: &DVector<f64>,
params_opt: Option<&DVector<f64>>|
-> DVector<f64> {
let mut ws = workspace.lock().unwrap();
let DynamicsWorkspace {
rotation_cache,
sh_v,
sh_w,
} = &mut *ws;
let epoch = epoch_initial + t;
let r_i2b = rotation_cache.get_or_compute(t, epoch);
let dx_orbital = Self::compute_dynamics(
t,
state,
epoch_initial,
&force_config,
params_opt.or(Some(¶ms)),
gravity_model.as_ref(),
r_i2b,
sh_v,
sh_w,
);
drop(ws);
let mut dx = DVector::zeros(state.len());
dx[0] = dx_orbital[0];
dx[1] = dx_orbital[1];
dx[2] = dx_orbital[2];
dx[3] = dx_orbital[3];
dx[4] = dx_orbital[4];
dx[5] = dx_orbital[5];
if let Some(ref add_dyn) = additional_dynamics
&& state.len() > 6
{
dx += add_dyn(t, state, params_opt);
}
dx
},
)
}
fn wrap_for_integrator(shared: SharedDynamics) -> DStateDynamics {
Box::new(
move |t: f64, state: &DVector<f64>, params: Option<&DVector<f64>>| -> DVector<f64> {
shared(t, state, params)
},
)
}
fn build_jacobian_provider(
shared_dynamics: SharedDynamics,
method: DifferenceMethod,
) -> Box<dyn crate::math::jacobian::DJacobianProvider> {
let dynamics_for_jacobian = Box::new(
move |t: f64, state: &DVector<f64>, params: Option<&DVector<f64>>| -> DVector<f64> {
shared_dynamics(t, state, params)
},
);
match method {
DifferenceMethod::Forward => {
Box::new(DNumericalJacobian::forward(dynamics_for_jacobian))
}
DifferenceMethod::Central => {
Box::new(DNumericalJacobian::central(dynamics_for_jacobian))
}
DifferenceMethod::Backward => {
Box::new(DNumericalJacobian::backward(dynamics_for_jacobian))
}
}
}
fn build_sensitivity_provider(
shared_dynamics: SharedDynamics,
method: DifferenceMethod,
) -> Box<dyn crate::math::sensitivity::DSensitivityProvider> {
let dynamics_with_params = Box::new(
move |t: f64, state: &DVector<f64>, params: &DVector<f64>| -> DVector<f64> {
shared_dynamics(t, state, Some(params))
},
);
match method {
DifferenceMethod::Forward => {
Box::new(DNumericalSensitivity::forward(dynamics_with_params))
}
DifferenceMethod::Central => {
Box::new(DNumericalSensitivity::central(dynamics_with_params))
}
DifferenceMethod::Backward => {
Box::new(DNumericalSensitivity::backward(dynamics_with_params))
}
}
}
fn compute_dynamics(
t: f64,
state: &DVector<f64>,
epoch_initial: Epoch,
force_config: &ForceModelConfig,
params_opt: Option<&DVector<f64>>,
gravity_model: Option<&Arc<GravityModel>>,
r_i2b: SMatrix3,
sh_v: &mut DMatrix<f64>,
sh_w: &mut DMatrix<f64>,
) -> Vector6<f64> {
let epoch = epoch_initial + t;
let r = Vector3::new(state[0], state[1], state[2]);
let v = Vector3::new(state[3], state[4], state[5]);
let x_eci = Vector6::new(state[0], state[1], state[2], state[3], state[4], state[5]);
let mut a_total = Vector3::zeros();
match &force_config.gravity {
GravityConfiguration::PointMass => {
a_total += accel_point_mass_gravity(r, Vector3::zeros(), GM_EARTH);
}
GravityConfiguration::EarthZonal { degree } => {
let r_ecef = r_i2b * r;
let a_ecef = accel_earth_zonal_gravity(r_ecef, degree.into());
a_total += r_i2b.transpose() * a_ecef;
}
GravityConfiguration::SphericalHarmonic {
source,
degree,
order,
} => {
match source {
GravityModelSource::Global => {
let global_model: std::sync::RwLockReadGuard<'_, Box<GravityModel>> =
get_global_gravity_model();
a_total += accel_gravity_spherical_harmonics_with_workspace(
r,
r_i2b,
&global_model,
*degree,
*order,
sh_v,
sh_w,
);
}
GravityModelSource::ModelType(_) => {
if let Some(model) = gravity_model {
a_total += accel_gravity_spherical_harmonics_with_workspace(
r,
r_i2b,
model.as_ref(),
*degree,
*order,
sh_v,
sh_w,
);
}
}
}
}
}
if let Some(drag_config) = &force_config.drag {
let mass = force_config
.mass
.as_ref()
.expect("Mass must be configured for drag calculation")
.get_value(params_opt);
let drag_area = drag_config.area.get_value(params_opt);
let cd = drag_config.cd.get_value(params_opt);
let density = match &drag_config.model {
AtmosphericModel::HarrisPriester => {
let r_sun = sun_position(epoch);
density_harris_priester(r, r_sun)
}
AtmosphericModel::NRLMSISE00 => {
let r_ecef = r_i2b * r;
density_nrlmsise00(&epoch, r_ecef).unwrap_or(0.0)
}
AtmosphericModel::Exponential {
scale_height,
rho0,
h0,
} => {
let altitude = r.norm() - R_EARTH;
let h_diff = altitude - h0;
rho0 * (-h_diff / scale_height).exp()
}
};
a_total += accel_drag(x_eci, density, mass, drag_area, cd, r_i2b);
}
if let Some(srp_config) = &force_config.srp {
let mass = force_config
.mass
.as_ref()
.expect("Mass must be configured for SRP calculation")
.get_value(params_opt);
let srp_area = srp_config.area.get_value(params_opt);
let cr = srp_config.cr.get_value(params_opt);
let r_sun = sun_position(epoch);
let mut a_srp = accel_solar_radiation_pressure(r, r_sun, mass, cr, srp_area, 4.56e-6);
let eclipse_factor = match srp_config.eclipse_model {
EclipseModel::None => 1.0,
EclipseModel::Cylindrical => eclipse_cylindrical(r, r_sun),
EclipseModel::Conical => eclipse_conical(r, r_sun),
};
a_srp *= eclipse_factor;
a_total += a_srp;
}
if let Some(tb_config) = &force_config.third_body {
for body in &tb_config.bodies {
a_total += accel_third_body(body.clone(), tb_config.ephemeris_source, epoch, r);
}
}
if force_config.relativity {
a_total += accel_relativity(x_eci);
}
Vector6::new(v[0], v[1], v[2], a_total[0], a_total[1], a_total[2])
}
pub fn set_trajectory_mode(&mut self, mode: TrajectoryMode) {
self.trajectory_mode = mode;
}
pub fn trajectory_mode(&self) -> TrajectoryMode {
self.trajectory_mode
}
pub fn set_output_step(&mut self, _step: f64) {
}
pub fn disable_output_step(&mut self) {
}
pub fn current_state_ref(&self) -> &DVector<f64> {
&self.x_curr
}
pub fn current_params(&self) -> &DVector<f64> {
&self.params
}
pub fn trajectory(&self) -> &DOrbitTrajectory {
&self.trajectory
}
pub fn stm(&self) -> Option<&DMatrix<f64>> {
self.stm.as_ref()
}
pub fn sensitivity(&self) -> Option<&DMatrix<f64>> {
self.sensitivity.as_ref()
}
pub fn stm_at_idx(&self, index: usize) -> Option<&DMatrix<f64>> {
self.trajectory.stm_at_idx(index)
}
pub fn sensitivity_at_idx(&self, index: usize) -> Option<&DMatrix<f64>> {
self.trajectory.sensitivity_at_idx(index)
}
pub fn stm_at(&self, epoch: Epoch) -> Option<DMatrix<f64>> {
self.trajectory.stm_at(epoch)
}
pub fn sensitivity_at(&self, epoch: Epoch) -> Option<DMatrix<f64>> {
self.trajectory.sensitivity_at(epoch)
}
pub fn reinitialize(
&mut self,
epoch: Epoch,
state: DVector<f64>,
covariance: Option<DMatrix<f64>>,
) {
self.t_rel = epoch - self.epoch_initial;
self.epoch_current = epoch;
self.x_curr = state;
match self.propagation_mode {
PropagationMode::StateOnly => {
}
PropagationMode::WithSTM => {
self.stm = Some(DMatrix::identity(self.state_dim, self.state_dim));
}
PropagationMode::WithSensitivity => {
let param_dim = self.params.len();
self.sensitivity = Some(DMatrix::zeros(self.state_dim, param_dim));
}
PropagationMode::WithSTMAndSensitivity => {
self.stm = Some(DMatrix::identity(self.state_dim, self.state_dim));
let param_dim = self.params.len();
self.sensitivity = Some(DMatrix::zeros(self.state_dim, param_dim));
}
}
self.integrator.reset_cache();
self.initial_covariance = covariance.clone();
self.current_covariance = covariance;
self.terminated = false;
}
pub fn current_covariance(&self) -> Option<&DMatrix<f64>> {
self.current_covariance.as_ref()
}
pub fn has_stm(&self) -> bool {
matches!(
self.propagation_mode,
PropagationMode::WithSTM | PropagationMode::WithSTMAndSensitivity
)
}
pub fn terminated(&self) -> bool {
self.terminated
}
pub fn add_event_detector(&mut self, detector: Box<dyn DEventDetector>) {
self.event_detectors.push(detector);
}
pub fn take_event_detectors(&mut self) -> Vec<Box<dyn DEventDetector>> {
std::mem::take(&mut self.event_detectors)
}
pub fn set_event_detectors(&mut self, detectors: Vec<Box<dyn DEventDetector>>) {
self.event_detectors = detectors;
}
pub fn take_event_log(&mut self) -> Vec<DDetectedEvent> {
std::mem::take(&mut self.event_log)
}
pub fn set_event_log(&mut self, log: Vec<DDetectedEvent>) {
self.event_log = log;
}
pub fn set_terminated(&mut self, terminated: bool) {
self.terminated = terminated;
}
pub fn is_terminated(&self) -> bool {
self.terminated
}
pub fn event_log(&self) -> &[DDetectedEvent] {
&self.event_log
}
pub fn events_by_name(&self, name: &str) -> Vec<&DDetectedEvent> {
self.query_events().by_name_contains(name).collect()
}
pub fn latest_event(&self) -> Option<&DDetectedEvent> {
self.event_log.last()
}
pub fn events_in_range(&self, start: Epoch, end: Epoch) -> Vec<&DDetectedEvent> {
self.query_events().in_time_range(start, end).collect()
}
pub fn query_events(
&self,
) -> crate::events::EventQuery<'_, std::slice::Iter<'_, DDetectedEvent>> {
crate::events::EventQuery::new(self.event_log.iter())
}
pub fn events_by_detector_index(&self, index: usize) -> Vec<&DDetectedEvent> {
self.query_events().by_detector_index(index).collect()
}
pub fn events_by_detector_index_in_range(
&self,
index: usize,
start: Epoch,
end: Epoch,
) -> Vec<&DDetectedEvent> {
self.query_events()
.by_detector_index(index)
.in_time_range(start, end)
.collect()
}
pub fn events_by_name_in_range(
&self,
name: &str,
start: Epoch,
end: Epoch,
) -> Vec<&DDetectedEvent> {
self.query_events()
.by_name_contains(name)
.in_time_range(start, end)
.collect()
}
pub fn clear_events(&mut self) {
self.event_log.clear();
}
pub fn reset_termination(&mut self) {
self.terminated = false;
}
fn step_once(&mut self) {
let epoch_prev = self.current_epoch();
let x_prev = self.x_curr.clone();
let dt_requested = self.dt_next;
let params_ref = if self.params.is_empty() {
None
} else {
Some(&self.params)
};
match self.propagation_mode {
PropagationMode::StateOnly => {
let result = self.integrator.step(
self.t_rel,
self.x_curr.clone(),
params_ref,
Some(dt_requested),
);
self.x_curr = result.state;
self.dt = result.dt_used;
self.dt_next = result.dt_next;
}
PropagationMode::WithSTM => {
let phi = self.stm.take().unwrap();
let result = self.integrator.step_with_varmat(
self.t_rel,
self.x_curr.clone(),
params_ref,
phi,
Some(dt_requested),
);
self.x_curr = result.state;
self.stm = result.phi;
self.dt = result.dt_used;
self.dt_next = result.dt_next;
if let Some(ref p0) = self.initial_covariance
&& let Some(ref phi_result) = self.stm
{
let p_new = phi_result * p0 * phi_result.transpose();
self.current_covariance = Some(p_new);
}
}
PropagationMode::WithSensitivity => {
let sens = self.sensitivity.take().unwrap();
let result = self.integrator.step_with_sensmat(
self.t_rel,
self.x_curr.clone(),
sens,
&self.params,
Some(dt_requested),
);
self.x_curr = result.state;
self.sensitivity = result.sens;
self.dt = result.dt_used;
self.dt_next = result.dt_next;
}
PropagationMode::WithSTMAndSensitivity => {
let phi = self.stm.take().unwrap();
let sens = self.sensitivity.take().unwrap();
let result = self.integrator.step_with_varmat_sensmat(
self.t_rel,
self.x_curr.clone(),
phi,
sens,
&self.params,
Some(dt_requested),
);
self.x_curr = result.state;
self.stm = result.phi.clone();
self.sensitivity = result.sens;
self.dt = result.dt_used;
self.dt_next = result.dt_next;
if let Some(ref p0) = self.initial_covariance
&& let Some(ref phi_result) = self.stm
{
let p_new = phi_result * p0 * phi_result.transpose();
self.current_covariance = Some(p_new);
}
}
}
self.t_rel += self.dt;
self.epoch_current = self.epoch_initial + self.t_rel;
let epoch_new = self.epoch_current;
let mut detected_events =
self.scan_all_events(epoch_prev, epoch_new, &x_prev, &self.x_curr);
self.refine_event_states(&mut detected_events, epoch_prev, &x_prev);
match self.process_events_smart(detected_events) {
EventProcessingResult::NoEvents => {
if self.should_store_state() {
let cov = self.current_covariance.clone();
let stm_to_store = if self.store_stm_history {
self.stm.clone()
} else {
None
};
let sens_to_store = if self.store_sensitivity_history {
self.sensitivity.clone()
} else {
None
};
let acc = self.compute_acceleration(epoch_new, &self.x_curr);
self.trajectory.add_full(
epoch_new,
self.x_curr.clone(),
cov,
stm_to_store,
sens_to_store,
acc,
);
}
}
EventProcessingResult::Restart { epoch, state } => {
self.t_rel = epoch - self.epoch_initial;
self.epoch_current = epoch;
self.x_curr = state;
self.integrator.reset_cache();
}
EventProcessingResult::Terminal {
epoch: term_epoch,
state: term_state,
} => {
self.epoch_current = term_epoch;
self.x_curr = term_state.clone();
if self.should_store_state() {
let cov = self.current_covariance.clone();
let stm_to_store = if self.store_stm_history {
self.stm.clone()
} else {
None
};
let sens_to_store = if self.store_sensitivity_history {
self.sensitivity.clone()
} else {
None
};
let acc = self.compute_acceleration(term_epoch, &term_state);
self.trajectory.add_full(
term_epoch,
term_state,
cov,
stm_to_store,
sens_to_store,
acc,
);
}
}
}
}
fn should_store_state(&self) -> bool {
match self.trajectory_mode {
TrajectoryMode::OutputStepsOnly => true,
TrajectoryMode::AllSteps => true,
TrajectoryMode::Disabled => false,
}
}
}
impl super::traits::DStatePropagator for DNumericalOrbitPropagator {
fn step_by(&mut self, step_size: f64) {
let target_t = self.t_rel + step_size;
self.process_initial_events();
let is_forward = step_size >= 0.0;
while !self.terminated {
if is_forward {
if self.t_rel >= target_t {
break;
}
} else if self.t_rel <= target_t {
break;
}
let remaining = target_t - self.t_rel;
if remaining.abs() <= 1e-12 {
break;
}
let dt_max = if is_forward {
remaining.min(self.dt_next.abs())
} else {
-remaining.abs().min(self.dt_next.abs())
};
let saved_dt_next = self.dt_next;
self.dt_next = dt_max;
self.step_once();
if self.dt_next.abs() > saved_dt_next.abs() {
self.dt_next = saved_dt_next;
}
}
}
fn propagate_to(&mut self, target_epoch: Epoch) {
let target_rel = target_epoch - self.epoch_initial;
self.process_initial_events();
let is_forward = target_rel >= self.t_rel;
while !self.terminated {
if is_forward {
if self.t_rel >= target_rel {
break;
}
} else if self.t_rel <= target_rel {
break;
}
let remaining = target_rel - self.t_rel;
if remaining.abs() <= 1e-12 {
break;
}
let dt_max = if is_forward {
remaining.min(self.dt_next.abs())
} else {
-remaining.abs().min(self.dt_next.abs())
};
let saved_dt_next = self.dt_next;
self.dt_next = dt_max;
self.step_once();
if self.dt_next.abs() >= saved_dt_next.abs() {
self.dt_next = saved_dt_next;
}
}
}
fn current_epoch(&self) -> Epoch {
self.epoch_current
}
fn current_state(&self) -> DVector<f64> {
self.x_curr.clone()
}
fn initial_epoch(&self) -> Epoch {
self.epoch_initial
}
fn initial_state(&self) -> DVector<f64> {
self.x_initial.clone()
}
fn state_dim(&self) -> usize {
self.state_dim
}
fn step_size(&self) -> f64 {
self.dt
}
fn set_step_size(&mut self, step_size: f64) {
self.dt = step_size;
self.dt_next = step_size;
}
fn reset(&mut self) {
self.t_rel = 0.0;
self.epoch_current = self.epoch_initial;
self.x_curr = self.x_initial.clone();
let initial_dt = self.dt; self.dt_next = initial_dt;
match self.propagation_mode {
PropagationMode::StateOnly => {
self.stm = None;
self.sensitivity = None;
}
PropagationMode::WithSTM => {
self.stm = Some(DMatrix::identity(self.state_dim, self.state_dim));
self.sensitivity = None;
}
PropagationMode::WithSensitivity => {
self.stm = None;
let param_dim = self.params.len();
self.sensitivity = Some(DMatrix::zeros(self.state_dim, param_dim));
}
PropagationMode::WithSTMAndSensitivity => {
self.stm = Some(DMatrix::identity(self.state_dim, self.state_dim));
let param_dim = self.params.len();
self.sensitivity = Some(DMatrix::zeros(self.state_dim, param_dim));
}
}
self.trajectory = DOrbitTrajectory::new(
self.state_dim,
OrbitFrame::ECI,
OrbitRepresentation::Cartesian,
None,
);
self.event_log.clear();
self.terminated = false;
for detector in &self.event_detectors {
detector.reset_processed();
}
}
fn set_eviction_policy_max_size(&mut self, max_size: usize) -> Result<(), BraheError> {
self.trajectory.set_eviction_policy_max_size(max_size)
}
fn set_eviction_policy_max_age(&mut self, max_age: f64) -> Result<(), BraheError> {
self.trajectory.set_eviction_policy_max_age(max_age)
}
}
impl InterpolationConfig for DNumericalOrbitPropagator {
fn with_interpolation_method(mut self, method: InterpolationMethod) -> Self {
self.interpolation_method = method;
self.trajectory.set_interpolation_method(method);
self
}
fn set_interpolation_method(&mut self, method: InterpolationMethod) {
self.interpolation_method = method;
self.trajectory.set_interpolation_method(method);
}
fn get_interpolation_method(&self) -> InterpolationMethod {
self.interpolation_method
}
}
impl CovarianceInterpolationConfig for DNumericalOrbitPropagator {
fn with_covariance_interpolation_method(
mut self,
method: CovarianceInterpolationMethod,
) -> Self {
self.covariance_interpolation_method = method;
self.trajectory.set_covariance_interpolation_method(method);
self
}
fn set_covariance_interpolation_method(&mut self, method: CovarianceInterpolationMethod) {
self.covariance_interpolation_method = method;
self.trajectory.set_covariance_interpolation_method(method);
}
fn get_covariance_interpolation_method(&self) -> CovarianceInterpolationMethod {
self.covariance_interpolation_method
}
}
impl DStateProvider for DNumericalOrbitPropagator {
fn state(&self, epoch: Epoch) -> Result<DVector<f64>, BraheError> {
if let Ok(state) = self.trajectory.interpolate(&epoch) {
return Ok(state);
}
if (self.current_epoch() - epoch).abs() < 1e-9 {
return Ok(self.x_curr.clone());
}
let start = self.trajectory.start_epoch().unwrap_or(self.epoch_initial);
let end = self.current_epoch();
Err(BraheError::OutOfBoundsError(format!(
"Cannot get state at epoch {}: outside propagator time range [{}, {}]. \
Call step_by() or propagate_to() to advance the propagator first.",
epoch, start, end
)))
}
fn state_dim(&self) -> usize {
self.state_dim
}
}
impl DOrbitStateProvider for DNumericalOrbitPropagator {
fn state_eci(&self, epoch: Epoch) -> Result<Vector6<f64>, BraheError> {
if let Ok(state) = self.trajectory.interpolate(&epoch) {
return Ok(state.fixed_rows::<6>(0).into());
}
if (self.current_epoch() - epoch).abs() < 1e-9 {
return Ok(self.x_curr.fixed_rows::<6>(0).into());
}
let start = self.trajectory.start_epoch().unwrap_or(self.epoch_initial);
let end = self.current_epoch();
Err(BraheError::OutOfBoundsError(format!(
"Cannot get state at epoch {}: outside propagator time range [{}, {}]. \
Call step_by() or propagate_to() to advance the propagator first.",
epoch, start, end
)))
}
fn state_ecef(&self, epoch: Epoch) -> Result<Vector6<f64>, BraheError> {
let eci_state = self.state_eci(epoch)?;
Ok(crate::frames::state_eci_to_ecef(epoch, eci_state))
}
fn state_gcrf(&self, epoch: Epoch) -> Result<Vector6<f64>, BraheError> {
self.state_eci(epoch)
}
fn state_itrf(&self, epoch: Epoch) -> Result<Vector6<f64>, BraheError> {
let gcrf_state = self.state_gcrf(epoch)?;
Ok(crate::frames::state_gcrf_to_itrf(epoch, gcrf_state))
}
fn state_eme2000(&self, epoch: Epoch) -> Result<Vector6<f64>, BraheError> {
let gcrf_state = self.state_gcrf(epoch)?;
Ok(crate::frames::state_gcrf_to_eme2000(gcrf_state))
}
fn state_koe_osc(
&self,
epoch: Epoch,
angle_format: AngleFormat,
) -> Result<Vector6<f64>, BraheError> {
let eci_state = self.state_eci(epoch)?;
Ok(state_eci_to_koe(eci_state, angle_format))
}
}
impl DCovarianceProvider for DNumericalOrbitPropagator {
fn covariance(&self, epoch: Epoch) -> Result<DMatrix<f64>, BraheError> {
if self.current_covariance.is_none() {
return Err(BraheError::InitializationError(
"Covariance not available: covariance tracking was not enabled for this propagator"
.to_string(),
));
}
if let Some(start) = self.trajectory.start_epoch()
&& epoch < start
{
return Err(BraheError::OutOfBoundsError(format!(
"Cannot get covariance at epoch {}: before trajectory start {}. \
Call step_by() or propagate_to() to advance the propagator first.",
epoch, start
)));
}
if let Some(end) = self.trajectory.end_epoch()
&& epoch > end
{
return Err(BraheError::OutOfBoundsError(format!(
"Cannot get covariance at epoch {}: after trajectory end {}. \
Call step_by() or propagate_to() to advance the propagator first.",
epoch, end
)));
}
self.trajectory.covariance_at(epoch).ok_or_else(|| {
BraheError::OutOfBoundsError(format!(
"Cannot get covariance at epoch {}: no covariance data available at this epoch",
epoch
))
})
}
fn covariance_dim(&self) -> usize {
self.state_dim
}
}
impl DOrbitCovarianceProvider for DNumericalOrbitPropagator {
fn covariance_eci(&self, epoch: Epoch) -> Result<DMatrix<f64>, BraheError> {
DCovarianceProvider::covariance(self, epoch)
}
fn covariance_gcrf(&self, epoch: Epoch) -> Result<DMatrix<f64>, BraheError> {
DCovarianceProvider::covariance(self, epoch)
}
fn covariance_rtn(&self, epoch: Epoch) -> Result<DMatrix<f64>, BraheError> {
let cov_eci = DCovarianceProvider::covariance(self, epoch)?;
let state_eci = self.state_eci(epoch)?;
let rot_eci_to_rtn = rotation_eci_to_rtn(state_eci);
let r = state_eci.fixed_rows::<3>(0);
let v = state_eci.fixed_rows::<3>(3);
let f_dot = (r.cross(&v)).norm() / (r.norm().powi(2));
let omega = nalgebra::Vector3::new(0.0, 0.0, f_dot);
let omega_skew = nalgebra::SMatrix::<f64, 3, 3>::new(
0.0, -omega[2], omega[1], omega[2], 0.0, -omega[0], -omega[1], omega[0], 0.0,
);
let j21 = -omega_skew * rot_eci_to_rtn;
let dim = cov_eci.nrows();
let mut jacobian = DMatrix::<f64>::zeros(dim, dim);
for i in 0..3 {
for j in 0..3 {
jacobian[(i, j)] = rot_eci_to_rtn[(i, j)];
jacobian[(3 + i, 3 + j)] = rot_eci_to_rtn[(i, j)];
}
}
for i in 3..6 {
for j in 0..3 {
jacobian[(i, j)] = j21[(i - 3, j)];
}
}
for i in 6..dim {
jacobian[(i, i)] = 1.0;
}
let cov_rtn = &jacobian * &cov_eci * jacobian.transpose();
Ok(cov_rtn)
}
}
impl Identifiable for DNumericalOrbitPropagator {
fn with_name(mut self, name: &str) -> Self {
self.name = Some(name.to_string());
self
}
fn with_uuid(mut self, uuid: uuid::Uuid) -> Self {
self.uuid = Some(uuid);
self
}
fn with_new_uuid(mut self) -> Self {
self.uuid = Some(uuid::Uuid::now_v7());
self
}
fn with_id(mut self, id: u64) -> Self {
self.id = Some(id);
self
}
fn with_identity(
mut self,
name: Option<&str>,
uuid: Option<uuid::Uuid>,
id: Option<u64>,
) -> Self {
self.name = name.map(|s| s.to_string());
self.uuid = uuid;
self.id = id;
self
}
fn set_identity(&mut self, name: Option<&str>, uuid: Option<uuid::Uuid>, id: Option<u64>) {
self.name = name.map(|s| s.to_string());
self.uuid = uuid;
self.id = id;
}
fn set_id(&mut self, id: Option<u64>) {
self.id = id;
}
fn set_name(&mut self, name: Option<&str>) {
self.name = name.map(|s| s.to_string());
}
fn generate_uuid(&mut self) {
self.uuid = Some(uuid::Uuid::now_v7());
}
fn get_id(&self) -> Option<u64> {
self.id
}
fn get_name(&self) -> Option<&str> {
self.name.as_deref()
}
fn get_uuid(&self) -> Option<uuid::Uuid> {
self.uuid
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::constants::units::AngleFormat;
use crate::coordinates::position_ecef_to_geodetic;
use crate::eop::{EOPExtrapolation, FileEOPProvider, set_global_eop_provider};
use crate::events::{DAltitudeEvent, DTimeEvent, EventDirection};
use crate::frames::position_eci_to_ecef;
use crate::propagators::NumericalPropagationConfig;
use crate::propagators::force_model_config::{
AtmosphericModel, DragConfiguration, EphemerisSource, GravityConfiguration,
ParameterSource, ThirdBody,
};
use crate::propagators::traits::DStatePropagator;
use crate::time::TimeSystem;
use crate::{orbital_period, state_koe_to_eci};
fn setup_global_test_eop() {
let eop = FileEOPProvider::from_default_standard(true, EOPExtrapolation::Hold).unwrap();
set_global_eop_provider(eop);
}
fn setup_global_test_space_weather() {
crate::utils::testing::setup_global_test_space_weather();
}
fn default_test_params() -> DVector<f64> {
DVector::from_vec(vec![
1000.0, 10.0, 2.2, 10.0, 1.3, ])
}
fn test_force_config_with_params() -> ForceModelConfig {
ForceModelConfig {
gravity: GravityConfiguration::PointMass,
drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::ParameterIndex(1),
cd: ParameterSource::ParameterIndex(2),
}),
srp: None,
third_body: None,
relativity: false,
mass: Some(ParameterSource::ParameterIndex(0)),
frame_transform: FrameTransformationModel::default(),
}
}
#[test]
fn test_dnumericalorbitpropagator_construction_default() {
setup_global_test_eop();
setup_global_test_space_weather();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0, 0.0,
7500.0,
0.0, ]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::default(),
Some(default_test_params()),
None,
None,
None,
);
assert!(prop.is_ok());
}
fn jgm3_force_config(degree: usize, order: usize) -> ForceModelConfig {
ForceModelConfig {
gravity: GravityConfiguration::SphericalHarmonic {
source: GravityModelSource::ModelType(
crate::orbit_dynamics::GravityModelType::JGM3,
),
degree,
order,
},
drag: None,
srp: None,
third_body: None,
relativity: false,
mass: None,
frame_transform: FrameTransformationModel::default(),
}
}
fn jgm3_initial_state() -> DVector<f64> {
DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0])
}
#[test]
#[serial_test::serial]
fn test_propagator_uses_gravity_model_cache_without_truncation() {
setup_global_test_eop();
crate::orbit_dynamics::clear_gravity_model_cache();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let prop = DNumericalOrbitPropagator::new(
epoch,
jgm3_initial_state(),
NumericalPropagationConfig::default(),
jgm3_force_config(70, 70),
None,
None,
None,
None,
)
.unwrap();
let from_cache = crate::orbit_dynamics::GravityModel::shared(
&crate::orbit_dynamics::GravityModelType::JGM3,
)
.unwrap();
assert_eq!(
std::sync::Arc::strong_count(&from_cache),
3,
"no-truncation propagator should share the cached Arc \
(cache + propagator + test binding = 3)"
);
drop(prop);
assert_eq!(
std::sync::Arc::strong_count(&from_cache),
2,
"dropping the propagator should release its Arc clone, \
leaving cache + this test binding"
);
}
#[test]
#[serial_test::serial]
fn test_propagator_clones_gravity_model_when_truncating() {
setup_global_test_eop();
crate::orbit_dynamics::clear_gravity_model_cache();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let _prop = DNumericalOrbitPropagator::new(
epoch,
jgm3_initial_state(),
NumericalPropagationConfig::default(),
jgm3_force_config(10, 10),
None,
None,
None,
None,
)
.unwrap();
let from_cache = crate::orbit_dynamics::GravityModel::shared(
&crate::orbit_dynamics::GravityModelType::JGM3,
)
.unwrap();
assert_eq!(
std::sync::Arc::strong_count(&from_cache),
2,
"truncating propagator should own a separate Arc, \
not share the cached one"
);
assert_eq!(from_cache.n_max, 70);
assert_eq!(from_cache.m_max, 70);
}
#[test]
fn test_rotation_cache_hit_returns_same_matrix() {
setup_global_test_eop();
let mut cache = RotationCache::new(FrameTransformationModel::FullEarthRotation);
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let first = cache.get_or_compute(0.0, epoch);
let second = cache.get_or_compute(0.0, epoch);
assert_eq!(
first, second,
"cache hit on the same t must return the bit-identical matrix"
);
}
#[test]
fn test_rotation_cache_matches_uncached_call() {
setup_global_test_eop();
let mut cache = RotationCache::new(FrameTransformationModel::FullEarthRotation);
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let cached = cache.get_or_compute(0.0, epoch);
let fresh = rotation_eci_to_ecef(epoch);
assert_eq!(
cached, fresh,
"cached rotation must match a fresh call to rotation_eci_to_ecef"
);
}
#[test]
fn test_rotation_cache_distinct_epochs_distinct_entries() {
setup_global_test_eop();
let mut cache = RotationCache::new(FrameTransformationModel::FullEarthRotation);
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let r_at_zero = cache.get_or_compute(0.0, epoch);
let r_at_3600 = cache.get_or_compute(3600.0, epoch + 3600.0);
assert_ne!(
r_at_zero, r_at_3600,
"rotations one hour apart should differ (Earth has rotated ~15 deg)"
);
assert_eq!(cache.get_or_compute(0.0, epoch), r_at_zero);
assert_eq!(cache.get_or_compute(3600.0, epoch + 3600.0), r_at_3600);
}
#[test]
fn test_rotation_cache_earth_rotation_only_path() {
setup_global_test_eop();
let mut cache = RotationCache::new(FrameTransformationModel::EarthRotationOnly);
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let cached = cache.get_or_compute(0.0, epoch);
let fresh = earth_rotation(epoch);
assert_eq!(
cached, fresh,
"EarthRotationOnly cache must dispatch to earth_rotation, not the full chain"
);
}
#[test]
fn test_rotation_cache_rk4_stage_pattern_hits() {
setup_global_test_eop();
let mut cache = RotationCache::new(FrameTransformationModel::FullEarthRotation);
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let h = 30.0;
let _stage0 = cache.get_or_compute(0.0, epoch);
let stage1 = cache.get_or_compute(h / 2.0, epoch + h / 2.0);
let stage2 = cache.get_or_compute(h / 2.0, epoch + h / 2.0);
assert_eq!(
stage1, stage2,
"RK4 stage 1 and stage 2 (both at t + h/2) must return the same cached matrix"
);
let stage3 = cache.get_or_compute(h, epoch + h);
let step2_stage0 = cache.get_or_compute(h, epoch + h);
assert_eq!(
stage3, step2_stage0,
"cross-step boundary (stage 3 of step N == stage 0 of step N+1) must be a cache hit"
);
}
#[test]
fn test_rotation_cache_capacity_evicts_oldest() {
setup_global_test_eop();
let mut cache = RotationCache::new(FrameTransformationModel::FullEarthRotation);
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let inserts: Vec<(f64, SMatrix3)> = (0..=ROTATION_CACHE_CAPACITY)
.map(|i| {
let t = i as f64;
let r = cache.get_or_compute(t, epoch + t);
(t, r)
})
.collect();
for &(t, expected) in &inserts[1..] {
assert_eq!(
cache.get_or_compute(t, epoch + t),
expected,
"entry at t={t} should still be cached"
);
}
let recomputed = cache.get_or_compute(0.0, epoch);
assert_eq!(
recomputed, inserts[0].1,
"post-eviction recompute must still match the original value"
);
}
#[test]
fn test_force_model_configuration_construction_variants() {
let _ = ForceModelConfig::default();
let _ = ForceModelConfig::high_fidelity();
let _ = ForceModelConfig::earth_gravity();
let _ = ForceModelConfig::leo_default();
let _ = ForceModelConfig::geo_default();
}
#[test]
fn test_dnumericalorbitpropagator_dstate_propagator_step_by() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0, 0.0,
7500.0,
0.0, ]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
fn use_trait<P: DStatePropagator>(prop: &mut P) {
let initial_epoch = prop.initial_epoch();
let _initial_state = prop.initial_state();
let _state_dim = prop.state_dim();
prop.step_by(60.0);
let current_epoch = prop.current_epoch();
let _current_state = prop.current_state();
assert!((current_epoch - initial_epoch - 60.0).abs() < 1.0);
}
use_trait(&mut prop);
}
#[test]
fn test_dnumericalorbitpropagator_dstatepropagator_propagate_to_forward() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0, 0.0,
7500.0,
0.0, ]);
let initial_pos = state[0];
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let target_time = 3600.0; let target_epoch = epoch + target_time;
prop.propagate_to(target_epoch);
assert!(
(prop.current_epoch() - target_epoch).abs() < 0.1,
"Should reach target epoch within 0.1s"
);
let current_state = prop.current_state();
assert_ne!(
current_state[0], initial_pos,
"Position should have changed"
);
}
#[test]
fn test_dnumericalorbitpropagator_dstatepropagator_step_by_backward() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0, 0.0,
7500.0,
0.0, ]);
let initial_pos = state[0];
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let step_size = -1800.0;
prop.step_by(step_size);
let time_diff: f64 = prop.current_epoch() - epoch;
assert!(
time_diff < -1790.0 && time_diff > -1810.0,
"Should have stepped backward ~1800s, got: {}",
time_diff
);
assert!(
prop.current_epoch() < epoch,
"Current epoch should be before initial epoch"
);
let current_state = prop.current_state();
assert_ne!(
current_state[0], initial_pos,
"Position should have changed during backward propagation"
);
}
#[test]
fn test_dnumericalorbitpropagator_dstatepropagator_propagate_to_backward() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0, 0.0,
7500.0,
0.0, ]);
let initial_pos = state[0];
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let target_time = -3600.0; let target_epoch = epoch + target_time;
prop.propagate_to(target_epoch);
assert!(
(prop.current_epoch() - target_epoch).abs() < 0.1,
"Should reach past epoch within 0.1s, got diff: {}",
prop.current_epoch() - target_epoch
);
assert!(
prop.current_epoch() < epoch,
"Current epoch should be before initial epoch"
);
let current_state = prop.current_state();
assert_ne!(
current_state[0], initial_pos,
"Position should have changed during backward propagation"
);
}
#[test]
fn test_dnumericalorbitpropagator_dstatepropagator_propagate_steps() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0, 0.0,
7500.0,
0.0, ]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let step_size = 60.0;
let num_steps = 10;
prop.propagate_steps(num_steps);
let elapsed = prop.current_epoch() - epoch;
let expected_min = step_size * (num_steps as f64) * 0.5; let expected_max = step_size * (num_steps as f64) * 2.0;
assert!(
elapsed > expected_min && elapsed < expected_max,
"Should have propagated approximately {} steps, elapsed: {}, expected range: [{}, {}]",
num_steps,
elapsed,
expected_min,
expected_max
);
let current_state = prop.current_state();
assert_ne!(current_state[0], state[0], "Position should have changed");
assert!(
prop.trajectory().len() > 0,
"Trajectory should have at least one entry"
);
}
#[test]
fn test_dnumericalorbitpropagator_dstatepropagator_reset() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0, 0.0,
7500.0,
0.0, ]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.propagate_to(epoch + 1800.0);
let state_after = prop.current_state();
assert_ne!(state_after[0], state[0], "State should have changed");
assert!(prop.trajectory().len() > 0, "Should have trajectory");
prop.reset();
assert_eq!(
prop.current_epoch(),
prop.initial_epoch(),
"Should be at initial epoch after reset"
);
let state_after_reset = prop.current_state();
for i in 0..6 {
assert!(
(state_after_reset[i] - state[i]).abs() < 1e-10,
"State should match initial state after reset, index {}: {} vs {}",
i,
state_after_reset[i],
state[i]
);
}
assert_eq!(
prop.trajectory().len(),
0,
"Trajectory should be empty after reset"
);
let target = epoch + 900.0;
prop.propagate_to(target);
let time_diff: f64 = prop.current_epoch() - target;
assert!(
time_diff.abs() < 0.1,
"Should propagate correctly after reset"
);
}
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_state_eci() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let initial_pos = state[0];
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let eci_state = prop.state_eci(query_epoch).unwrap();
assert_eq!(eci_state.len(), 6);
assert_ne!(eci_state[0], initial_pos);
let pos_mag = (eci_state[0].powi(2) + eci_state[1].powi(2) + eci_state[2].powi(2)).sqrt();
assert!(pos_mag > R_EARTH && pos_mag < R_EARTH + 1000e3);
}
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_state_ecef() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let ecef_state = prop.state_ecef(query_epoch).unwrap();
assert_eq!(ecef_state.len(), 6);
let pos_mag =
(ecef_state[0].powi(2) + ecef_state[1].powi(2) + ecef_state[2].powi(2)).sqrt();
assert!(pos_mag > R_EARTH && pos_mag < R_EARTH + 1000e3);
}
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_state_gcrf() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let gcrf_state = prop.state_gcrf(query_epoch).unwrap();
assert_eq!(gcrf_state.len(), 6);
let pos_mag =
(gcrf_state[0].powi(2) + gcrf_state[1].powi(2) + gcrf_state[2].powi(2)).sqrt();
assert!(pos_mag > R_EARTH && pos_mag < R_EARTH + 1000e3);
}
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_state_itrf() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let itrf_state = prop.state_itrf(query_epoch).unwrap();
assert_eq!(itrf_state.len(), 6);
let pos_mag =
(itrf_state[0].powi(2) + itrf_state[1].powi(2) + itrf_state[2].powi(2)).sqrt();
assert!(pos_mag > R_EARTH && pos_mag < R_EARTH + 1000e3);
}
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_state_eme2000() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let eme2000_state = prop.state_eme2000(query_epoch).unwrap();
assert_eq!(eme2000_state.len(), 6);
let pos_mag =
(eme2000_state[0].powi(2) + eme2000_state[1].powi(2) + eme2000_state[2].powi(2)).sqrt();
assert!(pos_mag > R_EARTH && pos_mag < R_EARTH + 1000e3);
}
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_osculating_elements_radians() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let elements = prop
.state_koe_osc(query_epoch, AngleFormat::Degrees)
.unwrap();
assert_eq!(elements.len(), 6);
assert!(elements[0] > R_EARTH + 300e3 && elements[0] < R_EARTH + 700e3);
assert!(elements[1] < 0.1);
assert!(elements[2].abs() < 0.1); }
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_osculating_elements_degrees() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let elements = prop
.state_koe_osc(query_epoch, AngleFormat::Degrees)
.unwrap();
assert_eq!(elements.len(), 6);
assert!(elements[0] > R_EARTH + 300e3 && elements[0] < R_EARTH + 700e3);
assert!(elements[1] < 0.1);
assert!(elements[2].abs() < 10.0); }
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_frame_conversion_roundtrip() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let eci_state = prop.state_eci(query_epoch).unwrap();
let ecef_state = prop.state_ecef(query_epoch).unwrap();
let eci_from_ecef = crate::frames::state_ecef_to_eci(query_epoch, ecef_state);
for i in 0..6 {
let diff = (eci_state[i] - eci_from_ecef[i]).abs();
let tolerance = if i < 3 { 1e-6 } else { 1e-9 }; assert!(
diff < tolerance,
"ECI ↔ ECEF round-trip failed at index {}: diff = {} m or m/s",
i,
diff
);
}
}
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_representation_conversion_roundtrip() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let cartesian = prop.state_eci(query_epoch).unwrap();
let keplerian = prop
.state_koe_osc(query_epoch, AngleFormat::Degrees)
.unwrap();
let cartesian_from_keplerian = state_koe_to_eci(keplerian, AngleFormat::Degrees);
for i in 0..6 {
let diff = (cartesian[i] - cartesian_from_keplerian[i]).abs();
let tolerance = if i < 3 { 1.0 } else { 1e-3 }; assert!(
diff < tolerance,
"Cartesian ↔ Keplerian round-trip failed at index {}: diff = {} m or m/s",
i,
diff
);
}
}
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_extended_state_preservation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0,
0.0,
7500.0,
0.0,
42.0,
99.0, ]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let eci_state = prop.state_eci(query_epoch).unwrap();
assert_eq!(eci_state.len(), 6, "ECI state should be 6D");
let full_state = prop.trajectory.interpolate(&query_epoch).unwrap();
assert_eq!(full_state.len(), 8, "Full state should preserve 8D");
assert_eq!(full_state.len(), 8);
assert_eq!(
full_state[6], state[6],
"Extended state index 6 should be preserved"
);
assert_eq!(
full_state[7], state[7],
"Extended state index 7 should be preserved"
);
}
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_angle_format_handling() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let elements_rad = prop
.state_koe_osc(query_epoch, AngleFormat::Radians)
.unwrap();
let elements_deg = prop
.state_koe_osc(query_epoch, AngleFormat::Degrees)
.unwrap();
assert!(
(elements_rad[0] - elements_deg[0]).abs() < 1e-6,
"Semi-major axis should match"
);
assert!(
(elements_rad[1] - elements_deg[1]).abs() < 1e-9,
"Eccentricity should match"
);
use crate::constants::math::RAD2DEG;
for i in 2..6 {
let expected_deg = elements_rad[i] * RAD2DEG;
let diff = (elements_deg[i] - expected_deg).abs();
assert!(
diff < 1e-9,
"Angle at index {} should differ by RAD2DEG factor: {} deg vs {} deg (diff: {})",
i,
elements_deg[i],
expected_deg,
diff
);
}
}
#[test]
fn test_dnumericalorbitpropagator_dorbitstateprovider_interpolation_accuracy() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 450.0; let interpolated_state = prop.state_eci(query_epoch).unwrap();
let pos_mag = (interpolated_state[0].powi(2)
+ interpolated_state[1].powi(2)
+ interpolated_state[2].powi(2))
.sqrt();
assert!(
pos_mag > R_EARTH && pos_mag < R_EARTH + 1000e3,
"Interpolated position should be in LEO range"
);
let initial_state = prop.initial_state();
let pos_diff = (interpolated_state[0] - initial_state[0]).abs();
assert!(
pos_diff > 1.0,
"Interpolated state should differ from initial"
);
}
#[test]
fn test_dnumericalorbitpropagator_dorbitcovarianceprovider_covariance_eci() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut initial_cov = DMatrix::zeros(6, 6);
for i in 0..3 {
initial_cov[(i, i)] = 100.0; }
for i in 3..6 {
initial_cov[(i, i)] = 1.0; }
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov.clone()),
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let cov_eci = prop.covariance_eci(query_epoch).unwrap();
assert_eq!(cov_eci.nrows(), 6);
assert_eq!(cov_eci.ncols(), 6);
assert!(
cov_eci[(0, 0)] > initial_cov[(0, 0)],
"Position variance should grow"
);
}
#[test]
fn test_dnumericalorbitpropagator_dorbitcovarianceprovider_covariance_gcrf() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let initial_cov = DMatrix::identity(6, 6) * 100.0;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov),
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let cov_gcrf = prop.covariance_gcrf(query_epoch).unwrap();
assert_eq!(cov_gcrf.nrows(), 6);
assert_eq!(cov_gcrf.ncols(), 6);
let cov_eci = prop.covariance_eci(query_epoch).unwrap();
for i in 0..6 {
for j in 0..6 {
let diff = (cov_gcrf[(i, j)] - cov_eci[(i, j)]).abs();
assert!(diff < 1e-6, "GCRF and ECI covariances should match closely");
}
}
}
#[test]
fn test_dnumericalorbitpropagator_dorbitcovarianceprovider_covariance_rtn() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let initial_cov = DMatrix::identity(6, 6) * 100.0;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov),
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let cov_rtn = prop.covariance_rtn(query_epoch).unwrap();
assert_eq!(cov_rtn.nrows(), 6);
assert_eq!(cov_rtn.ncols(), 6);
let cov_eci = prop.covariance_eci(query_epoch).unwrap();
let mut has_difference = false;
for i in 0..6 {
for j in 0..6 {
if (cov_rtn[(i, j)] - cov_eci[(i, j)]).abs() > 1.0 {
has_difference = true;
break;
}
}
}
assert!(has_difference, "RTN covariance should differ from ECI");
for i in 0..6 {
for j in 0..6 {
let diff = (cov_rtn[(i, j)] - cov_rtn[(j, i)]).abs();
assert!(
diff < 1e-6,
"RTN covariance should be symmetric at ({},{}): diff = {}",
i,
j,
diff
);
}
}
}
#[test]
fn test_dnumericalorbitpropagator_dorbitcovarianceprovider_interpolation_accuracy() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut initial_cov = DMatrix::zeros(6, 6);
for i in 0..3 {
initial_cov[(i, i)] = 100.0; }
for i in 3..6 {
initial_cov[(i, i)] = 1.0; }
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov.clone()),
)
.unwrap();
for _ in 0..10 {
prop.step_by(180.0); }
let query_epoch1 = epoch + 360.0; let query_epoch2 = epoch + 900.0; let query_epoch3 = epoch + 1440.0;
let cov1 = prop.covariance_eci(query_epoch1).unwrap();
let cov2 = prop.covariance_eci(query_epoch2).unwrap();
let cov3 = prop.covariance_eci(query_epoch3).unwrap();
assert_eq!(cov1.nrows(), 6);
assert_eq!(cov2.nrows(), 6);
assert_eq!(cov3.nrows(), 6);
let var1 = cov1[(0, 0)];
let var2 = cov2[(0, 0)];
let var3 = cov3[(0, 0)];
assert!(
var1 > initial_cov[(0, 0)],
"First covariance should be larger than initial"
);
assert!(
var2 > var1,
"Second covariance should be larger than first: {} > {}",
var2,
var1
);
assert!(
var3 > var2,
"Third covariance should be larger than second: {} > {}",
var3,
var2
);
}
#[test]
#[cfg_attr(not(feature = "manual"), ignore)] fn test_dnumericalorbitpropagator_dorbitcovarianceprovider_positive_definiteness() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut initial_cov = DMatrix::zeros(6, 6);
for i in 0..6 {
initial_cov[(i, i)] = 100.0; }
for i in 0..5 {
initial_cov[(i, i + 1)] = 10.0;
initial_cov[(i + 1, i)] = 10.0;
}
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov),
)
.unwrap();
prop.step_by(1800.0);
let query_epoch = epoch + 900.0;
let cov_eci = prop.covariance_eci(query_epoch).unwrap();
let cov_gcrf = prop.covariance_gcrf(query_epoch).unwrap();
let cov_rtn = prop.covariance_rtn(query_epoch).unwrap();
for i in 0..6 {
assert!(
cov_eci[(i, i)] > 0.0,
"ECI diagonal element {} should be positive",
i
);
assert!(
cov_gcrf[(i, i)] > 0.0,
"GCRF diagonal element {} should be positive",
i
);
assert!(
cov_rtn[(i, i)] > 0.0,
"RTN diagonal element {} should be positive",
i
);
}
for i in 0..6 {
for j in 0..6 {
assert!(
(cov_eci[(i, j)] - cov_eci[(j, i)]).abs() < 1e-9,
"ECI covariance should be symmetric"
);
assert!(
(cov_rtn[(i, j)] - cov_rtn[(j, i)]).abs() < 1e-6,
"RTN covariance should be symmetric"
);
}
}
}
#[test]
fn test_dnumericalorbitpropagator_dorbitcovarianceprovider_error_handling() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop_no_cov = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None, )
.unwrap();
prop_no_cov.step_by(1800.0);
let query_epoch = epoch + 900.0;
let result = prop_no_cov.covariance_eci(query_epoch);
assert!(result.is_err(), "Should error when covariance not enabled");
let initial_cov = DMatrix::identity(6, 6) * 100.0;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov),
)
.unwrap();
prop.step_by(1800.0);
let before_epoch = epoch - 1000.0; let result = prop.covariance_eci(before_epoch);
assert!(
result.is_err(),
"Should error for epoch before trajectory start"
);
}
#[test]
fn test_dnumericalorbitpropagator_rejects_hermite_quintic_without_accelerations() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let config = NumericalPropagationConfig::default()
.with_interpolation_method(InterpolationMethod::HermiteQuintic);
let result = DNumericalOrbitPropagator::new(
epoch,
state,
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
);
let err = result
.err()
.expect("expected propagator construction to fail");
let msg = err.to_string();
assert!(
msg.contains("HermiteQuintic interpolation requires store_accelerations = true"),
"unexpected error message: {msg}"
);
assert!(
msg.contains("with_store_accelerations(true)"),
"error should mention the fix: {msg}"
);
}
#[test]
fn test_dnumericalorbitpropagator_interpolationconfig_with_method() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap()
.with_interpolation_method(InterpolationMethod::Linear);
assert_eq!(prop.get_interpolation_method(), InterpolationMethod::Linear);
}
#[test]
fn test_dnumericalorbitpropagator_interpolationconfig_set_and_get() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let initial_method = prop.get_interpolation_method();
assert_eq!(initial_method, InterpolationMethod::Linear);
prop.set_interpolation_method(InterpolationMethod::HermiteCubic);
assert_eq!(
prop.get_interpolation_method(),
InterpolationMethod::HermiteCubic
);
}
#[test]
fn test_dnumericalorbitpropagator_covarianceinterpolationconfig_with_method() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let initial_cov = DMatrix::identity(6, 6) * 100.0;
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov),
)
.unwrap()
.with_covariance_interpolation_method(CovarianceInterpolationMethod::MatrixSquareRoot);
assert_eq!(
prop.get_covariance_interpolation_method(),
CovarianceInterpolationMethod::MatrixSquareRoot
);
}
#[test]
fn test_dnumericalorbitpropagator_covarianceinterpolationconfig_set_and_get() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let initial_cov = DMatrix::identity(6, 6) * 100.0;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov),
)
.unwrap();
let initial_method = prop.get_covariance_interpolation_method();
assert_eq!(
initial_method,
CovarianceInterpolationMethod::TwoWasserstein
);
prop.set_covariance_interpolation_method(CovarianceInterpolationMethod::MatrixSquareRoot);
assert_eq!(
prop.get_covariance_interpolation_method(),
CovarianceInterpolationMethod::MatrixSquareRoot
);
}
#[test]
fn test_dnumericalorbitpropagator_identifiable_with_name() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap()
.with_name("TestSat");
assert_eq!(prop.get_name(), Some("TestSat"));
}
#[test]
fn test_dnumericalorbitpropagator_identifiable_with_id() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap()
.with_id(12345);
assert_eq!(prop.get_id(), Some(12345));
}
#[test]
fn test_dnumericalorbitpropagator_identifiable_with_uuid() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let test_uuid = uuid::Uuid::now_v7();
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap()
.with_uuid(test_uuid);
assert_eq!(prop.get_uuid(), Some(test_uuid));
}
#[test]
fn test_dnumericalorbitpropagator_identifiable_with_new_uuid() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap()
.with_new_uuid();
assert!(prop.get_uuid().is_some());
}
#[test]
fn test_dnumericalorbitpropagator_identifiable_setters() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert_eq!(prop.get_name(), None);
assert_eq!(prop.get_id(), None);
assert!(prop.get_uuid().is_some());
prop.set_name(Some("UpdatedSat"));
assert_eq!(prop.get_name(), Some("UpdatedSat"));
prop.set_id(Some(99999));
assert_eq!(prop.get_id(), Some(99999));
prop.generate_uuid();
assert!(prop.get_uuid().is_some());
prop.set_name(None);
assert_eq!(prop.get_name(), None);
}
#[test]
fn test_dnumericalorbitpropagator_identifiable_with_identity() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let test_uuid = uuid::Uuid::now_v7();
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap()
.with_identity(Some("CompleteSat"), Some(test_uuid), Some(42));
assert_eq!(prop.get_name(), Some("CompleteSat"));
assert_eq!(prop.get_uuid(), Some(test_uuid));
assert_eq!(prop.get_id(), Some(42));
}
#[test]
fn test_dnumericalorbitpropagator_identifiable_set_identity() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let test_uuid = uuid::Uuid::now_v7();
prop.set_identity(Some("ModifiedSat"), Some(test_uuid), Some(777));
assert_eq!(prop.get_name(), Some("ModifiedSat"));
assert_eq!(prop.get_uuid(), Some(test_uuid));
assert_eq!(prop.get_id(), Some(777));
prop.set_identity(None, None, None);
assert_eq!(prop.get_name(), None);
assert_eq!(prop.get_uuid(), None);
assert_eq!(prop.get_id(), None);
}
#[test]
fn test_dnumericalorbitpropagator_identifiable_persistence_through_propagation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let test_uuid = uuid::Uuid::now_v7();
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap()
.with_name("PersistentSat")
.with_id(123)
.with_uuid(test_uuid);
assert_eq!(prop.get_name(), Some("PersistentSat"));
assert_eq!(prop.get_id(), Some(123));
assert_eq!(prop.get_uuid(), Some(test_uuid));
prop.step_by(1800.0);
assert_eq!(prop.get_name(), Some("PersistentSat"));
assert_eq!(prop.get_id(), Some(123));
assert_eq!(prop.get_uuid(), Some(test_uuid));
prop.reset();
assert_eq!(prop.get_name(), Some("PersistentSat"));
assert_eq!(prop.get_id(), Some(123));
assert_eq!(prop.get_uuid(), Some(test_uuid));
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_api_methods() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert_eq!(prop.event_log().len(), 0);
assert!(!prop.terminated());
assert!(prop.latest_event().is_none());
let time_event = DTimeEvent::new(epoch + 3600.0, "Test Event");
prop.add_event_detector(Box::new(time_event));
assert_eq!(prop.event_log().len(), 0);
prop.terminated = true;
assert!(prop.terminated());
prop.reset_termination();
assert!(!prop.terminated());
prop.clear_events();
assert_eq!(prop.event_log().len(), 0);
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_time_event() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let event_time = epoch + 1800.0;
let time_event = DTimeEvent::new(event_time, "30 Minute Mark");
prop.add_event_detector(Box::new(time_event));
prop.propagate_to(epoch + 3600.0);
let events = prop.event_log();
assert_eq!(events.len(), 1);
assert_eq!(events[0].name, "30 Minute Mark");
let event_epoch = events[0].window_open;
assert!((event_epoch - event_time).abs() < 0.1);
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_altitude_event() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let a = R_EARTH + 500e3; let e = 0.02; let i = 0.0;
let raan = 0.0;
let argp = 0.0;
let ta = 0.0;
let oe = DVector::from_vec(vec![a, e, i, raan, argp, ta]);
let state = state_koe_to_eci(
Vector6::from_column_slice(oe.as_slice()),
AngleFormat::Radians,
);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_column_slice(state.as_slice()),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let alt_event = DAltitudeEvent::new(450e3, "Low Alt", EventDirection::Any);
prop.add_event_detector(Box::new(alt_event));
let period = 2.0 * orbital_period(a);
prop.propagate_to(epoch + period);
let events = prop.event_log();
assert!(
!events.is_empty(),
"Expected at least 1 altitude crossing (should be 2), got {}. Period: {} s",
events.len(),
period
);
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_no_altitude_events() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let a = R_EARTH + 500e3; let e = 0.00; let i = 0.0;
let raan = 0.0;
let argp = 0.0;
let ta = 0.0;
let oe = DVector::from_vec(vec![a, e, i, raan, argp, ta]);
let state = state_koe_to_eci(
Vector6::from_column_slice(oe.as_slice()),
AngleFormat::Radians,
);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_column_slice(state.as_slice()),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let alt_event = DAltitudeEvent::new(450e3, "Low Alt", EventDirection::Any);
prop.add_event_detector(Box::new(alt_event));
let period = 2.0 * orbital_period(a);
prop.propagate_to(epoch + period);
let events = prop.event_log();
assert!(
events.is_empty(),
"Expected at least 0 altitude crossings, got {}. Period: {} s",
events.len(),
period
);
}
#[test]
fn test_dnumericalorbitpropagator_value_event_matches_altitude_event() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let a = R_EARTH + 500e3;
let e = 0.02;
let i = 97.8f64;
let raan = 0.0;
let argp = 0.0;
let mean_anomaly = 0.0;
let oe = DVector::from_vec(vec![a, e, i, raan, argp, mean_anomaly]);
let state = state_koe_to_eci(
Vector6::from_column_slice(oe.as_slice()),
AngleFormat::Degrees,
);
let mut prop_builtin = DNumericalOrbitPropagator::new(
epoch,
DVector::from_column_slice(state.as_slice()),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let mut prop_manual = DNumericalOrbitPropagator::new(
epoch,
DVector::from_column_slice(state.as_slice()),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let alt_event_builtin = DAltitudeEvent::new(450e3, "Altitude", EventDirection::Any);
prop_builtin.add_event_detector(Box::new(alt_event_builtin));
use crate::events::DValueEvent;
use nalgebra::Vector3;
let alt_event_manual = DValueEvent::new(
"ManualAltitude",
|t: Epoch, state: &DVector<f64>, _params: Option<&DVector<f64>>| {
let r_eci = Vector3::new(state[0], state[1], state[2]);
let r_ecef = position_eci_to_ecef(t, r_eci);
let geodetic = position_ecef_to_geodetic(r_ecef, AngleFormat::Radians);
geodetic[2]
},
450e3,
EventDirection::Any,
);
prop_manual.add_event_detector(Box::new(alt_event_manual));
let period = 2.0 * orbital_period(a);
prop_builtin.propagate_to(epoch + period);
prop_manual.propagate_to(epoch + period);
let events_builtin = prop_builtin.events_by_name("Altitude");
let events_manual = prop_manual.events_by_name("ManualAltitude");
assert_eq!(
events_builtin.len(),
events_manual.len(),
"Built-in and manual value events should detect same number of crossings"
);
assert!(
!events_builtin.is_empty(),
"Should detect at least one altitude crossing"
);
for (builtin, manual) in events_builtin.iter().zip(events_manual.iter()) {
let time_diff = (builtin.window_open - manual.window_open).abs();
assert!(
time_diff < 0.1,
"Event times should match within 0.1s, got diff: {:.6}s",
time_diff
);
}
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_callback_state_mutation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let maneuver_time = epoch + 1800.0;
let delta_v = 10.0;
let maneuver = DTimeEvent::new(maneuver_time, "Maneuver").with_callback(Box::new(
move |_t, state, _params| {
let mut new_state = state.clone();
new_state[3] += delta_v; (Some(new_state), None, EventAction::Continue)
},
));
prop.add_event_detector(Box::new(maneuver));
prop.propagate_to(epoch + 3600.0);
assert_eq!(prop.event_log().len(), 1);
assert!(!prop.terminated());
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_callback_parameter_mutation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
test_force_config_with_params(),
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
let initial_mass = prop.current_params()[0];
let event_time = epoch + 1800.0;
let new_mass = 500.0;
let mass_update = DTimeEvent::new(event_time, "Mass Update").with_callback(Box::new(
move |_t, _state, params| {
let mut new_params = params.unwrap().clone();
new_params[0] = new_mass; (None, Some(new_params), EventAction::Continue)
},
));
prop.add_event_detector(Box::new(mass_update));
prop.propagate_to(epoch + 3600.0);
let final_mass = prop.current_params()[0];
assert_ne!(initial_mass, final_mass);
assert_eq!(final_mass, new_mass);
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_terminal_event() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let terminal_time = epoch + 1800.0;
let terminal_event = DTimeEvent::new(terminal_time, "Terminal").set_terminal();
prop.add_event_detector(Box::new(terminal_event));
prop.propagate_to(epoch + 3600.0);
let events = prop.event_log();
println!("Events detected: {}", events.len());
if !events.is_empty() {
println!(
"Event: {} at {}",
events[0].name,
events[0].window_open - epoch
);
}
assert!(
prop.terminated(),
"Propagator should be terminated but isn't. Events detected: {}",
events.len()
);
let current = prop.current_epoch();
assert!(
(current - terminal_time).abs() < 10.0,
"Current epoch {} is not close to terminal time {}",
current - epoch,
terminal_time - epoch
); }
#[test]
fn test_dnumericalorbitpropagator_event_detection_multiple_no_callbacks() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 1200.0, "Event 1")));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 1800.0, "Event 2")));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 2400.0, "Event 3")));
prop.propagate_to(epoch + 3600.0);
let events = prop.event_log();
assert_eq!(events.len(), 3);
assert_eq!(events[0].name, "Event 1");
assert_eq!(events[1].name, "Event 2");
assert_eq!(events[2].name, "Event 3");
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_smart_processing() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 600.0, "Event 1 - No CB")));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 1200.0, "Event 2 - No CB")));
let callback_event = DTimeEvent::new(epoch + 1800.0, "Event 3 - WITH CB").with_callback(
Box::new(|_t, _state, _params| (None, None, EventAction::Continue)),
);
prop.add_event_detector(Box::new(callback_event));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 2400.0, "Event 4 - No CB")));
prop.propagate_to(epoch + 3600.0);
let events = prop.event_log();
assert_eq!(events.len(), 4);
assert_eq!(events[0].name, "Event 1 - No CB");
assert_eq!(events[1].name, "Event 2 - No CB");
assert_eq!(events[2].name, "Event 3 - WITH CB");
assert_eq!(events[3].name, "Event 4 - No CB");
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_at_initial_epoch() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let event_time = epoch + 1.0; prop.add_event_detector(Box::new(DTimeEvent::new(event_time, "Near-Initial Event")));
prop.step_by(1800.0);
let events = prop.event_log();
assert_eq!(
events.len(),
1,
"Event near initial epoch should be detected"
);
assert_eq!(events[0].name, "Near-Initial Event");
assert!((events[0].window_open - event_time).abs() < 0.1);
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_at_final_epoch() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let target_epoch = epoch + 1800.0;
prop.add_event_detector(Box::new(DTimeEvent::new(target_epoch, "Final Epoch Event")));
prop.propagate_to(target_epoch);
let events = prop.event_log();
assert_eq!(events.len(), 1, "Event at final epoch should be detected");
assert_eq!(events[0].name, "Final Epoch Event");
assert!((events[0].window_open - target_epoch).abs() < 1e-6);
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_simultaneous_events() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let event_time = epoch + 900.0;
prop.add_event_detector(Box::new(DTimeEvent::new(event_time, "Event A")));
prop.add_event_detector(Box::new(DTimeEvent::new(event_time, "Event B")));
prop.add_event_detector(Box::new(DTimeEvent::new(event_time, "Event C")));
prop.step_by(1800.0);
let events = prop.event_log();
assert_eq!(
events.len(),
3,
"All simultaneous events should be detected"
);
for event in events {
assert!((event.window_open - event_time).abs() < 1e-6);
}
let names: Vec<&str> = events.iter().map(|e| e.name.as_str()).collect();
assert!(names.contains(&"Event A"));
assert!(names.contains(&"Event B"));
assert!(names.contains(&"Event C"));
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_rapid_crossings() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
for i in 1..=10 {
let event_time = epoch + (i as f64) * 1.0;
prop.add_event_detector(Box::new(DTimeEvent::new(
event_time,
format!("Rapid Event {}", i),
)));
}
prop.step_by(360.0);
let events = prop.event_log();
assert_eq!(
events.len(),
10,
"All rapid crossing events should be detected"
);
for i in 0..9 {
assert!(
events[i].window_open < events[i + 1].window_open,
"Events should be in chronological order"
);
}
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_with_backward_propagation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 900.0, "Forward Event")));
prop.step_by(1800.0);
assert_eq!(prop.event_log().len(), 1);
let forward_final_epoch = prop.current_epoch();
prop.step_by(-900.0);
assert!(
prop.current_epoch() < forward_final_epoch,
"Backward propagation should move time backwards"
);
assert!(
!prop.event_log().is_empty(),
"Event log should have at least the original event"
);
assert!(
prop.event_log().iter().any(|e| e.name == "Forward Event"),
"Forward Event should be in event log"
);
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_log_persistence_across_reset_termination() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let terminal_event = DTimeEvent::new(epoch + 900.0, "Terminal Event").set_terminal();
prop.add_event_detector(Box::new(terminal_event));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 1800.0, "Post-Terminal")));
prop.step_by(3600.0);
assert_eq!(prop.event_log().len(), 1);
assert!(prop.terminated());
prop.reset_termination();
assert!(!prop.terminated());
assert_eq!(
prop.event_log().len(),
1,
"Event log should persist after reset_termination"
);
prop.step_by(1800.0);
let events = prop.event_log();
assert_eq!(
events.len(),
2,
"Should detect post-terminal event after reset"
);
assert_eq!(events[0].name, "Terminal Event");
assert_eq!(events[1].name, "Post-Terminal");
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_clear_vs_remove() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 600.0, "Event 1")));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 1200.0, "Event 2")));
prop.step_by(900.0);
assert_eq!(prop.event_log().len(), 1);
prop.clear_events();
assert_eq!(
prop.event_log().len(),
0,
"Event log should be empty after clear"
);
prop.step_by(900.0);
assert_eq!(
prop.event_log().len(),
1,
"Should detect new event after clearing log"
);
assert_eq!(prop.event_log()[0].name, "Event 2");
prop.reset();
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 600.0, "Reset Event")));
prop.step_by(900.0);
let events = prop.event_log();
assert!(!events.is_empty(), "Should detect events after reset");
let reset_event_found = events.iter().any(|e| e.name == "Reset Event");
assert!(reset_event_found, "Should find reset event");
let event1_found = events.iter().any(|e| e.name == "Event 1");
assert!(event1_found, "Should find Event 1 from before reset");
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_multiple_callbacks_same_step() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let a = R_EARTH + 500e3; let e = 0.02; let i = 0.0;
let raan = 0.0;
let argp = 0.0;
let ta = 0.0;
let oe = DVector::from_vec(vec![a, e, i, raan, argp, ta]);
let state = state_koe_to_eci(
Vector6::from_column_slice(oe.as_slice()),
AngleFormat::Radians,
);
use crate::integrators::IntegratorConfig;
use crate::propagators::IntegratorMethod;
let prop_config = NumericalPropagationConfig {
method: IntegratorMethod::RK4,
integrator: IntegratorConfig::fixed_step(900.0),
..Default::default()
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_column_slice(state.as_slice()),
prop_config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let callback1_executed = Arc::new(std::sync::atomic::AtomicBool::new(false));
let callback2_executed = Arc::new(std::sync::atomic::AtomicBool::new(false));
let cb1_flag = callback1_executed.clone();
let cb2_flag = callback2_executed.clone();
use crate::events::EventAction;
let event1 = DAltitudeEvent::new(410e3, "Event 1 - 410 km", EventDirection::Any)
.with_callback(Box::new(move |_t, state, _params| {
cb1_flag.store(true, std::sync::atomic::Ordering::SeqCst);
let mut new_state = state.clone();
new_state[3] += 5.0; (Some(new_state), None, EventAction::Continue)
}));
let event2 = DAltitudeEvent::new(380e3, "Event 2 - 380 km", EventDirection::Any)
.with_callback(Box::new(move |_t, state, _params| {
cb2_flag.store(true, std::sync::atomic::Ordering::SeqCst);
let mut new_state = state.clone();
new_state[4] += 10.0; (Some(new_state), None, EventAction::Continue)
}));
prop.add_event_detector(Box::new(event1));
prop.add_event_detector(Box::new(event2));
let period = orbital_period(a);
prop.propagate_to(epoch + period / 4.0);
assert!(
callback2_executed.load(std::sync::atomic::Ordering::SeqCst),
"Event 2 callback (380 km, first chronological) should execute"
);
assert!(
callback1_executed.load(std::sync::atomic::Ordering::SeqCst),
"Event 1 callback (410 km, second chronological) should also execute"
);
let events = prop.event_log();
assert!(
events.len() >= 2,
"At least 2 events should be in the log, got {}",
events.len()
);
let has_380km = events.iter().any(|e| e.name == "Event 2 - 380 km");
let has_410km = events.iter().any(|e| e.name == "Event 1 - 410 km");
assert!(has_380km, "380 km event should be in the log");
assert!(has_410km, "410 km event should be in the log");
}
#[test]
fn test_dnumericalorbitpropagator_event_detection_time_events_no_infinite_loop() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let callback1_count = Arc::new(std::sync::atomic::AtomicUsize::new(0));
let callback2_count = Arc::new(std::sync::atomic::AtomicUsize::new(0));
let cb1_count = callback1_count.clone();
let cb2_count = callback2_count.clone();
use crate::events::EventAction;
let event1 = DTimeEvent::new(epoch + 1800.0, "Event 1").with_callback(Box::new(
move |_t, _state, _params| {
cb1_count.fetch_add(1, std::sync::atomic::Ordering::SeqCst);
(None, None, EventAction::Continue)
},
));
let event2 = DTimeEvent::new(epoch + 1800.1, "Event 2").with_callback(Box::new(
move |_t, _state, _params| {
cb2_count.fetch_add(1, std::sync::atomic::Ordering::SeqCst);
(None, None, EventAction::Continue)
},
));
prop.add_event_detector(Box::new(event1));
prop.add_event_detector(Box::new(event2));
prop.propagate_to(epoch + 3600.0);
let cb1_executions = callback1_count.load(std::sync::atomic::Ordering::SeqCst);
let cb2_executions = callback2_count.load(std::sync::atomic::Ordering::SeqCst);
assert_eq!(
cb1_executions, 1,
"Callback 1 should execute exactly once, got {}",
cb1_executions
);
assert_eq!(
cb2_executions, 1,
"Callback 2 should execute exactly once, got {}",
cb2_executions
);
let events = prop.event_log();
assert_eq!(events.len(), 2, "Both events should be logged");
}
#[test]
fn test_dnumericalorbitpropagator_event_callback_accuracy_at_different_times() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = Vector6::new(
R_EARTH + 700e3,
0.001,
97.8_f64.to_radians(),
0.26,
0.52,
0.79,
);
let state = state_koe_to_eci(oe, AngleFormat::Radians);
let x_init = DVector::from_column_slice(state.as_slice());
let force_config = ForceModelConfig::earth_gravity();
let prop_config = NumericalPropagationConfig::default();
let target_epoch = epoch + 120.0;
let dv = DVector::from_vec(vec![0.0, 0.0, 0.0, 0.1, -0.05, 0.02]);
let offsets = [0.0, 0.0001, 1.0, 10.0, 30.0, 60.0];
for &offset in &offsets {
let baseline_state = if offset <= 1e-12 {
let mut x_burned = x_init.clone();
for i in 3..6 {
x_burned[i] += dv[i];
}
let mut p = DNumericalOrbitPropagator::new(
epoch,
x_burned,
prop_config.clone(),
force_config.clone(),
None,
None,
None,
None,
)
.unwrap();
p.propagate_to(target_epoch);
p.current_state().clone()
} else {
let mut p1 = DNumericalOrbitPropagator::new(
epoch,
x_init.clone(),
prop_config.clone(),
force_config.clone(),
None,
None,
None,
None,
)
.unwrap();
p1.propagate_to(epoch + offset);
let mut state_at_burn = p1.current_state().clone();
for i in 3..6 {
state_at_burn[i] += dv[i];
}
let mut p2 = DNumericalOrbitPropagator::new(
epoch + offset,
state_at_burn,
prop_config.clone(),
force_config.clone(),
None,
None,
None,
None,
)
.unwrap();
p2.propagate_to(target_epoch);
p2.current_state().clone()
};
let dv_clone = dv.clone();
let mut prop_te = DNumericalOrbitPropagator::new(
epoch,
x_init.clone(),
prop_config.clone(),
force_config.clone(),
None,
None,
None,
None,
)
.unwrap();
use crate::events::EventAction;
let event = DTimeEvent::new(epoch + offset, format!("Burn_T+{offset}")).with_callback(
Box::new(move |_t, state, _params| {
let mut new_state = state.clone();
for i in 3..6 {
new_state[i] += dv_clone[i];
}
(Some(new_state), None, EventAction::Continue)
}),
);
prop_te.add_event_detector(Box::new(event));
prop_te.propagate_to(target_epoch);
let te_state = prop_te.current_state();
let pos_err = (te_state.rows(0, 3) - baseline_state.rows(0, 3)).norm();
let vel_err = (te_state.rows(3, 3) - baseline_state.rows(3, 3)).norm();
assert!(
pos_err < 1e-2,
"T+{offset}s: position error {pos_err:.2e} m exceeds 10 mm threshold"
);
assert!(
vel_err < 1e-5,
"T+{offset}s: velocity error {vel_err:.2e} m/s exceeds 10 μm/s threshold"
);
let events = prop_te.event_log();
assert_eq!(
events.len(),
1,
"T+{offset}s: expected 1 event, got {}",
events.len()
);
}
}
#[test]
fn test_dnumericalorbitpropagator_event_callback_multi_impulse_accuracy() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = Vector6::new(
R_EARTH + 700e3,
0.001,
97.8_f64.to_radians(),
0.26,
0.52,
0.79,
);
let state = state_koe_to_eci(oe, AngleFormat::Radians);
let x_init = DVector::from_column_slice(state.as_slice());
let force_config = ForceModelConfig::earth_gravity();
let prop_config = NumericalPropagationConfig::default();
let dvs = [
DVector::from_vec(vec![0.0, 0.0, 0.0, 0.1, -0.05, 0.02]),
DVector::from_vec(vec![0.0, 0.0, 0.0, -0.03, 0.08, -0.01]),
DVector::from_vec(vec![0.0, 0.0, 0.0, 0.02, -0.02, 0.005]),
];
let burn_times = [0.0001, 60.0, 120.0]; let target_epoch = epoch + 180.0;
let mut x_state = x_init.clone();
let mut cur_ep = epoch;
for (dv, &t_burn) in dvs.iter().zip(&burn_times) {
let burn_ep = epoch + t_burn;
if burn_ep > cur_ep {
let mut p = DNumericalOrbitPropagator::new(
cur_ep,
x_state.clone(),
prop_config.clone(),
force_config.clone(),
None,
None,
None,
None,
)
.unwrap();
p.propagate_to(burn_ep);
x_state = p.current_state().clone();
cur_ep = burn_ep;
}
for i in 3..6 {
x_state[i] += dv[i];
}
}
let mut p_final = DNumericalOrbitPropagator::new(
cur_ep,
x_state,
prop_config.clone(),
force_config.clone(),
None,
None,
None,
None,
)
.unwrap();
p_final.propagate_to(target_epoch);
let baseline_state = p_final.current_state().clone();
use crate::events::EventAction;
let mut prop_te = DNumericalOrbitPropagator::new(
epoch,
x_init.clone(),
prop_config.clone(),
force_config.clone(),
None,
None,
None,
None,
)
.unwrap();
for (k, (dv, &t_burn)) in dvs.iter().zip(&burn_times).enumerate() {
let dv_clone = dv.clone();
let event = DTimeEvent::new(epoch + t_burn, format!("Burn_{k}")).with_callback(
Box::new(move |_t, state, _params| {
let mut new_state = state.clone();
for i in 3..6 {
new_state[i] += dv_clone[i];
}
(Some(new_state), None, EventAction::Continue)
}),
);
prop_te.add_event_detector(Box::new(event));
}
prop_te.propagate_to(target_epoch);
let te_state = prop_te.current_state();
let pos_err = (te_state.rows(0, 3) - baseline_state.rows(0, 3)).norm();
let vel_err = (te_state.rows(3, 3) - baseline_state.rows(3, 3)).norm();
assert!(
pos_err < 1e-2,
"Multi-impulse position error {pos_err:.2e} m exceeds 10 mm threshold"
);
assert!(
vel_err < 1e-5,
"Multi-impulse velocity error {vel_err:.2e} m/s exceeds 10 μm/s threshold"
);
let events = prop_te.event_log();
assert_eq!(events.len(), 3, "Expected 3 events, got {}", events.len());
}
#[test]
fn test_dnumericalorbitpropagator_continuous_control_via_control_input() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let a = R_EARTH + 500e3;
let e = 0.0;
let i = 0.0;
let raan = 0.0;
let argp = 0.0;
let ta = 0.0;
let oe = DVector::from_vec(vec![a, e, i, raan, argp, ta]);
let state = state_koe_to_eci(
Vector6::from_column_slice(oe.as_slice()),
AngleFormat::Radians,
);
let state = DVector::from_column_slice(state.as_slice());
let control_fn: crate::integrators::traits::DControlInput =
Some(Box::new(|_t, state, _params| {
let v = Vector3::new(state[3], state[4], state[5]);
let v_mag = v.norm();
let a_control = if v_mag > 1e-6 {
v * (0.0005 / v_mag) } else {
Vector3::zeros()
};
let mut dx = DVector::zeros(state.len());
dx[3] = a_control[0]; dx[4] = a_control[1]; dx[5] = a_control[2]; dx
}));
let mut prop_ref = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None, None,
)
.unwrap();
let mut prop_test = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
control_fn, None,
)
.unwrap();
let period = orbital_period(a);
prop_ref.propagate_to(epoch + period);
prop_test.propagate_to(epoch + period);
let state_ref = prop_ref.current_state();
let state_test = prop_test.current_state();
let oe_ref = state_eci_to_koe(
Vector6::from_column_slice(state_ref.as_slice()),
AngleFormat::Radians,
);
let oe_test = state_eci_to_koe(
Vector6::from_column_slice(state_test.as_slice()),
AngleFormat::Radians,
);
let a_ref = oe_ref[0];
let a_test = oe_test[0];
let delta_a = a_test - a_ref;
assert!(
delta_a > 4000.0,
"Semi-major axis should increase by > 4 km with tangential thrust, got {:.2} km",
delta_a / 1000.0
);
let expected_delta_a = 5100.0; let tolerance = 0.20; assert!(
(delta_a - expected_delta_a).abs() < expected_delta_a * tolerance,
"Semi-major axis increase should be ~{:.2} km ± 20%, got {:.2} km",
expected_delta_a / 1000.0,
delta_a / 1000.0
);
}
#[test]
fn test_dnumericalorbitpropagator_accuracy_vs_keplerian() {
use crate::propagators::KeplerianPropagator;
use crate::propagators::traits::SStatePropagator;
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.001, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut num_prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let mut kep_prop = KeplerianPropagator::new(
epoch,
state,
OrbitFrame::ECI,
OrbitRepresentation::Cartesian,
None,
60.0,
);
let orbital_period = orbital_period(oe[0]);
num_prop.step_by(orbital_period);
kep_prop.step_by(orbital_period);
let num_final = num_prop.current_state();
let kep_final = kep_prop.current_state();
let pos_error = (num_final.fixed_rows::<3>(0) - kep_final.fixed_rows::<3>(0)).norm();
assert!(
pos_error < 10.0,
"Position error vs Keplerian should be < 10 m, got {:.3} m",
pos_error
);
let vel_error = (num_final.fixed_rows::<3>(3) - kep_final.fixed_rows::<3>(3)).norm();
assert!(
vel_error < 0.01,
"Velocity error vs Keplerian should be < 10 mm/s, got {:.6} m/s",
vel_error
);
}
#[test]
fn test_dnumericalorbitpropagator_accuracy_orbital_period() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let a = R_EARTH + 500e3;
let oe = nalgebra::Vector6::new(a, 0.01, 45.0_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let t_theory = orbital_period(a);
let n_orbits = 5.0;
prop.step_by(n_orbits * t_theory);
let t_actual = (prop.current_epoch() - epoch) / n_orbits;
let period_error = ((t_actual - t_theory) / t_theory).abs();
assert!(
period_error < 1e-3,
"Orbital period error should be < 0.1%, got {:.6}%",
period_error * 100.0
);
}
#[test]
fn test_dnumericalorbitpropagator_accuracy_adaptive_step_behavior() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut config_tight = NumericalPropagationConfig::default();
config_tight.integrator.abs_tol = 1e-12;
config_tight.integrator.rel_tol = 1e-12;
let mut config_loose = NumericalPropagationConfig::default();
config_loose.integrator.abs_tol = 1e-8;
config_loose.integrator.rel_tol = 1e-8;
let mut prop_tight = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
config_tight,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let mut prop_loose = DNumericalOrbitPropagator::new(
epoch,
state,
config_loose,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop_tight.step_by(1800.0);
prop_loose.step_by(1800.0);
let state_tight = prop_tight.current_state();
let state_loose = prop_loose.current_state();
let difference = (state_tight.fixed_rows::<3>(0) - state_loose.fixed_rows::<3>(0)).norm();
assert!(
difference > 0.0 && difference < 1.0,
"Tolerance should affect accuracy, diff = {:.6} m",
difference
);
}
#[test]
fn test_dnumericalorbitpropagator_accuracy_energy_conservation_point_mass() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.001, 45.0_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let r0 = state.fixed_rows::<3>(0).norm();
let v0 = state.fixed_rows::<3>(3).norm();
let e0 = 0.5 * v0 * v0 - GM_EARTH / r0;
let orbital_period = orbital_period(oe[0]);
prop.step_by(10.0 * orbital_period);
let final_state = prop.current_state();
let r1 = final_state.fixed_rows::<3>(0).norm();
let v1 = final_state.fixed_rows::<3>(3).norm();
let e1 = 0.5 * v1 * v1 - GM_EARTH / r1;
let rel_energy_error = ((e1 - e0) / e0).abs();
assert!(
rel_energy_error < 1e-6,
"Energy conservation error should be < 1e-6 (10 orbits), got {:.3e}",
rel_energy_error
);
}
#[test]
fn test_dnumericalorbitpropagator_accuracy_angular_momentum_conservation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 55.0_f64, 30.0_f64, 45.0_f64, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let r0 = state.fixed_rows::<3>(0);
let v0 = state.fixed_rows::<3>(3);
let h0 = r0.cross(&v0);
let orbital_period = orbital_period(oe[0]);
prop.step_by(10.0 * orbital_period);
let final_state = prop.current_state();
let r1 = final_state.fixed_rows::<3>(0);
let v1 = final_state.fixed_rows::<3>(3);
let h1 = r1.cross(&v1);
let h_error = (h1 - h0).norm();
let h0_mag = h0.norm();
let rel_h_error = h_error / h0_mag;
assert!(
rel_h_error < 5e-6,
"Angular momentum conservation error should be < 5e-6 (10 orbits), got {:.3e}",
rel_h_error
);
}
#[test]
#[cfg_attr(not(feature = "ci"), ignore)]
fn test_dnumericalorbitpropagator_accuracy_energy_drift_long_term() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.001, 45.0_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let r0 = state.fixed_rows::<3>(0).norm();
let v0 = state.fixed_rows::<3>(3).norm();
let e0 = 0.5 * v0 * v0 - GM_EARTH / r0;
let orbital_period = orbital_period(oe[0]);
prop.step_by(100.0 * orbital_period);
let final_state = prop.current_state();
let r1 = final_state.fixed_rows::<3>(0).norm();
let v1 = final_state.fixed_rows::<3>(3).norm();
let e1 = 0.5 * v1 * v1 - GM_EARTH / r1;
let rel_energy_drift = ((e1 - e0) / e0).abs();
assert!(
rel_energy_drift < 1e-5,
"Energy drift over 100 orbits should be < 1e-5, got {:.3e}",
rel_energy_drift
);
}
#[test]
#[cfg_attr(not(feature = "ci"), ignore)]
fn test_dnumericalorbitpropagator_accuracy_orbital_stability_long_term() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe_initial =
nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 55.0_f64, 30.0_f64, 45.0_f64, 0.0);
let state = state_koe_to_eci(oe_initial, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe_initial[0]);
prop.step_by(100.0 * orbital_period);
let final_state = prop.current_state();
let oe_final =
state_eci_to_koe(final_state.fixed_rows::<6>(0).into(), AngleFormat::Radians);
let a_drift = (oe_final[0] - oe_initial[0]).abs();
assert!(
a_drift < 10.0,
"Semi-major axis drift over 100 orbits should be < 10 m, got {:.1} m",
a_drift
);
let e_drift = (oe_final[1] - oe_initial[1]).abs();
assert!(
e_drift < 0.001,
"Eccentricity drift over 100 orbits should be < 0.001, got {:.6}",
e_drift
);
}
#[test]
fn test_dnumericalorbitpropagator_accuracy_leo_regime() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(
R_EARTH + 400e3, 0.001,
51.6_f64, 0.0,
0.0,
0.0,
);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(3.0 * 86400.0);
assert!(
prop.current_epoch() > epoch,
"Should successfully propagate LEO orbit for 3 days"
);
let final_state = prop.current_state();
let r_final = final_state.fixed_rows::<3>(0).norm();
assert!(
r_final > R_EARTH + 300e3 && r_final < R_EARTH + 500e3,
"LEO orbit should remain stable, altitude = {:.1} km",
(r_final - R_EARTH) / 1000.0
);
}
#[test]
#[cfg_attr(not(feature = "ci"), ignore)]
fn test_dnumericalorbitpropagator_accuracy_geo_regime() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(
R_EARTH + 35786e3, 0.001,
0.1_f64, 0.0,
0.0,
0.0,
);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(7.0 * 86400.0);
assert!(
prop.current_epoch() > epoch,
"Should successfully propagate GEO orbit for 7 days"
);
let final_state = prop.current_state();
let r_final = final_state.fixed_rows::<3>(0).norm();
let altitude_km = (r_final - R_EARTH) / 1000.0;
assert!(
(altitude_km - 35786.0).abs() < 100.0,
"GEO altitude should remain near 35786 km, got {:.1} km",
altitude_km
);
}
#[test]
fn test_dnumericalorbitpropagator_accuracy_heo_regime() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(
26554e3, 0.72, 63.4_f64, 0.0, 270.0_f64, 0.0,
);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe[0]);
prop.step_by(2.0 * orbital_period);
assert!(
prop.current_epoch() > epoch,
"Should successfully propagate HEO orbit for 2 periods"
);
let final_state = prop.current_state();
let oe_final =
state_eci_to_koe(final_state.fixed_rows::<6>(0).into(), AngleFormat::Radians);
let e_error = (oe_final[1] - oe[1]).abs();
assert!(
e_error < 0.01,
"HEO eccentricity should be preserved, error = {:.6}",
e_error
);
}
#[test]
fn test_dnumericalorbitpropagator_accuracy_near_circular_stability() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(
R_EARTH + 600e3,
0.0001, 45.0_f64,
0.0,
0.0,
0.0,
);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe[0]);
prop.step_by(50.0 * orbital_period);
let final_state = prop.current_state();
let oe_final =
state_eci_to_koe(final_state.fixed_rows::<6>(0).into(), AngleFormat::Radians);
assert!(
oe_final[1] < 0.001,
"Orbit should remain nearly circular, e_final = {:.6}",
oe_final[1]
);
let r_final = final_state.fixed_rows::<3>(0).norm();
let r_variation = (r_final - oe[0]).abs();
assert!(
r_variation < 5000.0,
"Radius variation should be small for circular orbit, got {:.1} m",
r_variation
);
}
#[test]
fn test_dnumericalorbitpropagator_edge_case_high_eccentricity() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(
R_EARTH + 15000e3, 0.95, 45.0_f64,
0.0,
0.0,
0.0,
);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe[0]);
prop.step_by(orbital_period);
assert!(
prop.current_epoch() > epoch,
"Should successfully propagate high eccentricity orbit"
);
let final_state = prop.current_state();
let oe_final =
state_eci_to_koe(final_state.fixed_rows::<6>(0).into(), AngleFormat::Radians);
assert!(
(oe_final[1] - oe[1]).abs() < 0.01,
"High eccentricity should be preserved, e_error = {:.6}",
(oe_final[1] - oe[1]).abs()
);
}
#[test]
fn test_dnumericalorbitpropagator_edge_case_equatorial_orbit() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(
R_EARTH + 500e3,
0.01,
0.0, 0.0,
0.0,
0.0,
);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe[0]);
prop.step_by(10.0 * orbital_period);
assert!(
prop.current_epoch() > epoch,
"Should successfully propagate equatorial orbit"
);
}
#[test]
fn test_dnumericalorbitpropagator_edge_case_polar_orbit() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(
R_EARTH + 800e3,
0.001,
90.0_f64, 0.0,
0.0,
0.0,
);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe[0]);
prop.step_by(10.0 * orbital_period);
let final_state = prop.current_state();
let oe_final =
state_eci_to_koe(final_state.fixed_rows::<6>(0).into(), AngleFormat::Degrees);
assert!(
(oe_final[2] - oe[2]).abs() < 0.01,
"Polar inclination should be preserved, i_error = {:.6} deg",
(oe_final[2] - oe[2]).abs()
);
}
#[test]
fn test_dnumericalorbitpropagator_edge_case_very_short_step() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let initial_state = prop.current_state().clone();
prop.step_by(0.01);
assert!(
(prop.current_epoch() - epoch) > 0.0,
"Should handle very short time steps"
);
let final_state = prop.current_state();
let state_diff = (final_state - initial_state).norm();
assert!(
state_diff > 0.0 && state_diff < 100.0,
"State should change slightly, diff = {:.6} m",
state_diff
);
}
#[test]
fn test_dnumericalorbitpropagator_edge_case_propagate_to_same_epoch() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.propagate_to(epoch);
let final_state = prop.current_state();
let state_diff = (final_state - state).norm();
assert!(
state_diff < 1e-10,
"State should remain unchanged when propagating to same epoch"
);
}
#[test]
fn test_dnumericalorbitpropagator_edge_case_backward_then_forward() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(-900.0);
let backward_epoch = prop.current_epoch();
assert!(backward_epoch < epoch, "Should propagate backward");
prop.propagate_to(epoch);
let final_state = prop.current_state();
let state_diff = (final_state - state).norm();
assert!(
state_diff < 100.0, "Should return close to original state after backward/forward, diff = {:.3} m",
state_diff
);
}
#[test]
fn test_dnumericalorbitpropagator_edge_case_single_step_propagation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.step_by(1.0);
assert_eq!(
(prop.current_epoch() - epoch),
1.0,
"Should take exactly one second step"
);
}
#[test]
fn test_dnumericalorbitpropagator_construction_with_custom_params() {
setup_global_test_eop();
setup_global_test_space_weather();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let force_config = ForceModelConfig::default();
let params = DVector::from_vec(vec![500.0, 5.0, 2.0, 8.0, 1.5]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config,
Some(params),
None,
None,
None,
);
assert!(prop.is_ok());
let prop = prop.unwrap();
assert_eq!(DStatePropagator::state_dim(&prop), 6);
assert_eq!(prop.initial_epoch(), epoch);
}
#[test]
fn test_dnumericalorbitpropagator_construction_ecef_frame() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let eci_state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::new(
epoch,
eci_state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
);
assert!(prop.is_ok());
let prop = prop.unwrap();
let internal_eci = prop.current_state();
assert!((internal_eci[0] - eci_state[0]).abs() < 1.0);
assert!((internal_eci[1] - eci_state[1]).abs() < 1.0);
assert!((internal_eci[2] - eci_state[2]).abs() < 1.0);
}
#[test]
fn test_dnumericalorbitpropagator_construction_extended_state() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let extended_state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0, 0.0,
7500.0,
0.0, 100.0, 200.0, ]);
let prop = DNumericalOrbitPropagator::new(
epoch,
extended_state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
);
assert!(prop.is_ok());
let prop = prop.unwrap();
assert_eq!(DStatePropagator::state_dim(&prop), 8);
let current = prop.current_state();
assert_eq!(current.len(), 8);
assert_eq!(current[6], 100.0);
assert_eq!(current[7], 200.0);
}
#[test]
fn test_dnumericalorbitpropagator_construction_with_additional_dynamics() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let extended_state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0,
0.0,
7500.0,
0.0,
1000.0, ]);
let additional_dynamics: DStateDynamics = Box::new(|_t, state, _params| {
let mut dx = DVector::zeros(state.len());
dx[6] = -0.1; dx
});
let prop = DNumericalOrbitPropagator::new(
epoch,
extended_state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
Some(additional_dynamics),
None,
None,
);
assert!(prop.is_ok());
let mut prop = prop.unwrap();
assert_eq!(DStatePropagator::state_dim(&prop), 7);
let initial_mass = prop.current_state()[6];
prop.step_by(10.0);
let final_mass = prop.current_state()[6];
assert!((final_mass - (initial_mass - 1.0)).abs() < 1e-3);
}
#[test]
fn test_dnumericalorbitpropagator_trajectory_stores_additional_states() {
use approx::assert_abs_diff_eq;
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let extended_state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0,
0.0,
7500.0,
0.0,
1000.0, 100.0, ]);
let additional_dynamics: DStateDynamics = Box::new(|_t, state, _params| {
let mut dx = DVector::zeros(state.len());
dx[6] = -0.1; dx[7] = -0.05; dx
});
let config = NumericalPropagationConfig {
interpolation_method: InterpolationMethod::Linear,
..Default::default()
};
let prop = DNumericalOrbitPropagator::new(
epoch,
extended_state.clone(),
config,
ForceModelConfig::earth_gravity(),
None,
Some(additional_dynamics),
None,
None,
);
assert!(prop.is_ok());
let mut prop = prop.unwrap();
prop.set_trajectory_mode(TrajectoryMode::AllSteps);
assert_eq!(DStatePropagator::state_dim(&prop), 8);
assert_eq!(prop.trajectory().dimension(), 8);
assert_eq!(prop.trajectory().additional_dimension(), 2);
prop.step_by(10.0);
let traj = prop.trajectory();
assert!(traj.len() >= 1);
let final_epoch = prop.current_epoch();
let orbital_state = prop.state_eci(final_epoch).unwrap();
assert_eq!(orbital_state.len(), 6);
let final_idx = traj.len() - 1;
let (_final_epoch, final_full_state) = traj.get(final_idx).unwrap();
assert_eq!(final_full_state.len(), 8);
let final_mass = final_full_state[6];
let final_fuel = final_full_state[7];
assert_abs_diff_eq!(final_mass, 999.0, epsilon = 0.2); assert_abs_diff_eq!(final_fuel, 99.5, epsilon = 0.2);
if traj.len() >= 2 {
let start = traj.start_epoch().unwrap();
let end = traj.end_epoch().unwrap();
let mid_epoch = start + (end - start) / 2.0;
let mid_state = traj.interpolate(&mid_epoch).unwrap();
assert_eq!(mid_state.len(), 8);
assert!(mid_state[6] < 1000.0); assert!(mid_state[7] < 100.0); }
}
#[test]
fn test_dnumericalorbitpropagator_construction_multiple_integrators() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
use crate::propagators::IntegratorMethod;
let prop_dp54 = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
);
assert!(prop_dp54.is_ok());
let prop_rkf45 = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::with_method(IntegratorMethod::RKF45),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
);
assert!(prop_rkf45.is_ok());
let prop_rkn = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::with_method(IntegratorMethod::RKN1210),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
);
assert!(prop_rkn.is_ok());
let prop_rk4 = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::with_method(IntegratorMethod::RK4),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
);
assert!(prop_rk4.is_ok());
}
#[test]
fn test_dnumericalorbitpropagator_current_epoch() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert_eq!(prop.current_epoch(), epoch);
prop.step_by(100.0);
let new_epoch = prop.current_epoch();
assert!((new_epoch - epoch - 100.0).abs() < 1.0);
}
#[test]
fn test_dnumericalorbitpropagator_current_state() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let eci_state = prop.current_state();
assert_eq!(eci_state.len(), 6);
for i in 0..6 {
assert!((eci_state[i] - state[i]).abs() < 1e-6);
}
}
#[test]
fn test_dnumericalorbitpropagator_current_params() {
setup_global_test_eop();
setup_global_test_space_weather();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let params = prop.current_params();
assert_eq!(params.len(), 0);
let custom_params = DVector::from_vec(vec![1000.0, 10.0, 2.2, 10.0, 1.3]);
let prop_with_params = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::default(),
Some(custom_params),
None,
None,
None,
)
.unwrap();
let params = prop_with_params.current_params();
assert_eq!(params.len(), 5);
assert_eq!(params[0], 1000.0); assert_eq!(params[1], 10.0); assert_eq!(params[2], 2.2); assert_eq!(params[3], 10.0); assert_eq!(params[4], 1.3); }
#[test]
fn test_dnumericalorbitpropagator_initial_epoch() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 12, 30, 45.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert_eq!(prop.initial_epoch(), epoch);
prop.step_by(1000.0);
assert_eq!(prop.initial_epoch(), epoch);
}
#[test]
fn test_dnumericalorbitpropagator_initial_state() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let initial = prop.initial_state();
for i in 0..6 {
assert!((initial[i] - state[i]).abs() < 1e-6);
}
prop.step_by(1000.0);
let initial_after = prop.initial_state();
for i in 0..6 {
assert!((initial_after[i] - state[i]).abs() < 1e-6);
}
}
#[test]
fn test_dnumericalorbitpropagator_state_dim() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state_6d = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop_6d = DNumericalOrbitPropagator::new(
epoch,
state_6d,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert_eq!(DStatePropagator::state_dim(&prop_6d), 6);
let state_8d = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0,
0.0,
7500.0,
0.0,
100.0,
200.0,
]);
let prop_8d = DNumericalOrbitPropagator::new(
epoch,
state_8d,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert_eq!(DStatePropagator::state_dim(&prop_8d), 8);
}
#[test]
fn test_dnumericalorbitpropagator_trajectory_access() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let traj = prop.trajectory();
assert_eq!(traj.len(), 1);
prop.step_by(100.0);
prop.step_by(100.0);
let traj_after = prop.trajectory();
assert!(traj_after.len() > 1);
}
#[test]
fn test_dnumericalorbitpropagator_stm_access() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop_no_stm = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert!(prop_no_stm.stm().is_none());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert!(prop.stm().is_some());
let stm = prop.stm().unwrap();
assert_eq!(stm.nrows(), 6);
assert_eq!(stm.ncols(), 6);
}
#[test]
fn test_dnumericalorbitpropagator_sensitivity_access() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop_no_sens = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert!(prop_no_sens.sensitivity().is_none());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true; let params = DVector::from_vec(vec![1000.0, 10.0, 2.2, 10.0, 1.3]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
config,
ForceModelConfig::earth_gravity(),
Some(params),
None,
None,
None,
)
.unwrap();
assert!(prop.sensitivity().is_some());
let sens = prop.sensitivity().unwrap();
assert_eq!(sens.nrows(), 6);
assert_eq!(sens.ncols(), 5);
}
#[test]
fn test_dnumericalorbitpropagator_terminated_flag() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert!(!prop.terminated());
let terminal_event = DTimeEvent::new(epoch + 100.0, "Terminal").set_terminal();
prop.add_event_detector(Box::new(terminal_event));
prop.propagate_to(epoch + 200.0);
assert!(prop.terminated());
}
#[test]
fn test_dnumericalorbitpropagator_set_trajectory_mode() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let force_config = ForceModelConfig {
gravity: GravityConfiguration::PointMass,
drag: None,
srp: None,
third_body: None,
relativity: false,
mass: None,
frame_transform: FrameTransformationModel::default(),
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
None,
None,
None,
None,
)
.unwrap();
prop.set_trajectory_mode(TrajectoryMode::AllSteps);
assert_eq!(prop.trajectory_mode(), TrajectoryMode::AllSteps);
prop.step_by(100.0);
prop.step_by(100.0);
assert!(prop.trajectory().len() > 0);
prop.reset();
prop.set_trajectory_mode(TrajectoryMode::Disabled);
assert_eq!(prop.trajectory_mode(), TrajectoryMode::Disabled);
prop.step_by(100.0);
prop.step_by(100.0);
assert_eq!(prop.trajectory().len(), 0);
prop.reset();
prop.set_trajectory_mode(TrajectoryMode::OutputStepsOnly);
assert_eq!(prop.trajectory_mode(), TrajectoryMode::OutputStepsOnly);
}
#[test]
fn test_dnumericalorbitpropagator_trajectory_mode_getter() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert_eq!(prop.trajectory_mode(), TrajectoryMode::AllSteps);
}
#[test]
fn test_dnumericalorbitpropagator_set_step_size() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.set_step_size(30.0);
assert_eq!(prop.step_size(), 30.0);
prop.set_step_size(120.0);
assert_eq!(prop.step_size(), 120.0);
}
#[test]
fn test_dnumericalorbitpropagator_step_size_getter() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let step_size = prop.step_size();
assert!(step_size > 0.0);
}
#[test]
fn test_dnumericalorbitpropagator_set_eviction_policy_max_size() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let result = prop.set_eviction_policy_max_size(10);
assert!(result.is_ok());
for _ in 0..20 {
prop.step_by(10.0);
}
let traj_len = prop.trajectory().len();
assert!(
traj_len <= 10,
"Trajectory length {} exceeds max size 10",
traj_len
);
}
#[test]
fn test_dnumericalorbitpropagator_set_eviction_policy_max_age() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let result = prop.set_eviction_policy_max_age(100.0);
assert!(result.is_ok());
for _ in 0..20 {
prop.step_by(10.0);
}
let traj = prop.trajectory();
if traj.len() > 0 {
let epochs = &traj.epochs;
let current_time = prop.current_epoch() - prop.initial_epoch();
let oldest_time = epochs[0] - prop.initial_epoch();
assert!(
current_time - oldest_time <= 100.0 + 10.0,
"Oldest state age {} exceeds max age 100s",
current_time - oldest_time
);
}
}
#[test]
fn test_dnumericalorbitpropagator_force_gravity_point_mass() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let force_config = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: None,
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
let initial_energy = compute_orbital_energy(&state);
prop.step_by(5400.0);
let final_state = prop.current_state();
let final_energy = compute_orbital_energy(&final_state);
assert!(
(final_energy - initial_energy).abs() / initial_energy.abs() < 1e-6,
"Energy not conserved: {} vs {}",
initial_energy,
final_energy
);
}
fn compute_orbital_energy(state: &DVector<f64>) -> f64 {
let r = (state[0].powi(2) + state[1].powi(2) + state[2].powi(2)).sqrt();
let v = (state[3].powi(2) + state[4].powi(2) + state[5].powi(2)).sqrt();
0.5 * v.powi(2) - GM_EARTH / r
}
#[test]
fn test_dnumericalorbitpropagator_force_gravity_spherical_harmonic() {
use crate::orbit_dynamics::gravity::GravityModelType;
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let force_config = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::SphericalHarmonic {
source: GravityModelSource::ModelType(GravityModelType::EGM2008_360),
degree: 4,
order: 4,
},
drag: None,
srp: None,
third_body: None,
relativity: false,
};
let prop_result = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
);
assert!(
prop_result.is_ok(),
"Spherical harmonic config should construct successfully"
);
}
#[test]
fn test_dnumericalorbitpropagator_force_gravity_j2_perturbation() {
use crate::orbit_dynamics::gravity::GravityModelType;
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(
R_EARTH + 700e3, 0.001, 97.8_f64, 0.0, 0.0, 0.0, );
let state_initial = state_koe_to_eci(oe, AngleFormat::Degrees);
let force_config = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::SphericalHarmonic {
source: GravityModelSource::ModelType(GravityModelType::EGM2008_360),
degree: 2,
order: 0,
},
drag: None,
srp: None,
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state_initial.as_slice().to_vec()),
NumericalPropagationConfig::default(),
force_config,
None,
None,
None,
None,
)
.unwrap();
prop.step_by(86400.0);
let state_final = prop.current_state();
let oe_final =
state_eci_to_koe(state_final.fixed_rows::<6>(0).into(), AngleFormat::Degrees);
let delta_raan = oe_final[3] - oe[3];
assert!(
delta_raan.abs() > 1e-6,
"J2 should cause RAAN drift, but delta = {} rad",
delta_raan
);
let delta_a = oe_final[0] - oe[0];
assert!(
delta_a.abs() < 10000.0,
"Semi-major axis drift should be reasonable, delta = {} m",
delta_a
);
}
#[test]
fn test_dnumericalorbitpropagator_force_gravity_degree_order_convergence() {
use crate::orbit_dynamics::gravity::GravityModelType;
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let configs = vec![(2, 0), (4, 4), (8, 8)];
let mut final_states = Vec::new();
for (degree, order) in configs {
let force_config = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::SphericalHarmonic {
source: GravityModelSource::ModelType(GravityModelType::EGM2008_360),
degree,
order,
},
drag: None,
srp: None,
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config,
None,
None,
None,
None,
)
.unwrap();
prop.step_by(5400.0);
let state_vec: Vector6<f64> = prop.current_state().fixed_rows::<6>(0).into();
final_states.push(state_vec);
}
let diff_2_4: f64 = (final_states[0] - final_states[1]).norm();
let diff_4_8: f64 = (final_states[1] - final_states[2]).norm();
assert!(
diff_4_8 < diff_2_4,
"Higher degree/order should converge: diff(2x0,4x4)={:.3}m > diff(4x4,8x8)={:.3}m",
diff_2_4,
diff_4_8
);
assert!(diff_2_4 > 10.0, "J2 vs 4x4 should differ significantly");
assert!(diff_4_8 < diff_2_4, "Higher orders should show convergence");
}
#[test]
fn test_dnumericalorbitpropagator_force_gravity_global_vs_modeltype() {
use crate::orbit_dynamics::gravity::GravityModelType;
use crate::utils::testing::setup_global_test_gravity_model;
setup_global_test_eop();
setup_global_test_gravity_model();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let force_config_modeltype = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::SphericalHarmonic {
source: GravityModelSource::ModelType(GravityModelType::EGM2008_360),
degree: 4,
order: 4,
},
drag: None,
srp: None,
third_body: None,
relativity: false,
};
let mut prop_modeltype = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config_modeltype,
None,
None,
None,
None,
)
.unwrap();
let force_config_global = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::SphericalHarmonic {
source: GravityModelSource::Global,
degree: 4,
order: 4,
},
drag: None,
srp: None,
third_body: None,
relativity: false,
};
let prop_global_result = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config_global,
None,
None,
None,
None,
);
let _ = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
prop_modeltype.step_by(1800.0);
}));
match prop_global_result {
Ok(mut prop_global) => {
let _ = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
prop_global.step_by(1800.0);
}));
}
Err(_) => {
}
}
}
#[test]
fn test_dnumericalorbitpropagator_force_drag_harris_priester() {
use crate::propagators::force_model_config::{DragConfiguration, ParameterSource};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 300e3,
0.0,
0.0, 0.0,
7700.0,
0.0,
]);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)),
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::Value(10.0),
cd: ParameterSource::Value(2.2),
}),
srp: None,
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
let initial_energy = compute_orbital_energy(&state);
prop.step_by(600.0);
let final_state = prop.current_state();
let final_energy = compute_orbital_energy(&final_state);
assert!(
final_energy < initial_energy,
"Drag should decrease orbital energy"
);
}
#[test]
fn test_dnumericalorbitpropagator_force_drag_exponential() {
use crate::propagators::force_model_config::{DragConfiguration, ParameterSource};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 300e3, 0.0, 0.0, 0.0, 7700.0, 0.0]);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)),
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: Some(DragConfiguration {
model: AtmosphericModel::Exponential {
scale_height: 60000.0, rho0: 3.614e-13, h0: 300000.0, },
area: ParameterSource::Value(10.0),
cd: ParameterSource::Value(2.2),
}),
srp: None,
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
let initial_energy = compute_orbital_energy(&state);
prop.step_by(600.0);
let final_energy = compute_orbital_energy(&prop.current_state());
assert!(final_energy < initial_energy);
}
#[test]
fn test_dnumericalorbitpropagator_force_drag_nrlmsise00() {
use crate::propagators::force_model_config::{DragConfiguration, ParameterSource};
setup_global_test_eop();
setup_global_test_space_weather();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 400e3, 0.0, 0.0, 0.0, 7700.0, 0.0]);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::Value(1000.0)),
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: Some(DragConfiguration {
model: AtmosphericModel::NRLMSISE00,
area: ParameterSource::Value(10.0),
cd: ParameterSource::Value(2.2),
}),
srp: None,
third_body: None,
relativity: false,
};
let prop_result = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
);
assert!(
prop_result.is_ok(),
"NRLMSISE00 config should construct successfully"
);
if let Ok(mut prop) = prop_result {
let _ = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
prop.step_by(600.0);
}));
}
}
#[test]
fn test_dnumericalorbitpropagator_force_drag_magnitude_direction() {
use crate::propagators::force_model_config::{DragConfiguration, ParameterSource};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 300e3,
0.0,
0.0, 0.0,
7700.0,
0.0, ]);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)),
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::Value(10.0),
cd: ParameterSource::Value(2.2),
}),
srp: None,
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop.step_by(60.0);
let final_state = prop.current_state();
let initial_vy = state[4];
let final_vy = final_state[4];
assert!(
final_vy.abs() < initial_vy.abs(),
"Drag should reduce velocity magnitude: |{}| should be < |{}|",
final_vy,
initial_vy
);
}
#[test]
fn test_dnumericalorbitpropagator_force_drag_orbital_decay() {
use crate::propagators::force_model_config::{DragConfiguration, ParameterSource};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe_initial = nalgebra::Vector6::new(
R_EARTH + 250e3, 0.001, 45.0_f64,
0.0,
0.0,
0.0,
);
let state = state_koe_to_eci(oe_initial, AngleFormat::Degrees);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)),
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::Value(10.0),
cd: ParameterSource::Value(2.2),
}),
srp: None,
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
DVector::from_vec(state.as_slice().to_vec()),
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe_initial[0]);
prop.step_by(10.0 * orbital_period);
let final_state = prop.current_state();
let oe_final =
state_eci_to_koe(final_state.fixed_rows::<6>(0).into(), AngleFormat::Degrees);
let delta_a = oe_final[0] - oe_initial[0];
assert!(
delta_a < -1000.0,
"Semi-major axis should decay by at least 1 km over 10 orbits at 250 km: delta = {} m",
delta_a
);
let altitude_initial = oe_initial[0] - R_EARTH;
let altitude_final = oe_final[0] - R_EARTH;
assert!(
altitude_final < altitude_initial,
"Altitude should decrease: {} km -> {} km",
altitude_initial / 1000.0,
altitude_final / 1000.0
);
}
#[test]
fn test_dnumericalorbitpropagator_force_srp_no_eclipse() {
use crate::propagators::force_model_config::{
ParameterSource, SolarRadiationPressureConfiguration,
};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 20000e3,
0.0,
0.0, 0.0,
4000.0,
0.0,
]);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)),
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: Some(SolarRadiationPressureConfiguration {
area: ParameterSource::Value(10.0),
cr: ParameterSource::Value(1.3),
eclipse_model: EclipseModel::None,
}),
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop.step_by(3600.0);
assert!(prop.current_epoch() > epoch);
}
#[test]
fn test_dnumericalorbitpropagator_force_srp_cylindrical_eclipse() {
use crate::propagators::force_model_config::{
ParameterSource, SolarRadiationPressureConfiguration,
};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)),
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: Some(SolarRadiationPressureConfiguration {
area: ParameterSource::Value(10.0),
cr: ParameterSource::Value(1.3),
eclipse_model: EclipseModel::Cylindrical,
}),
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop.step_by(5400.0);
assert!(prop.current_epoch() > epoch);
}
#[test]
fn test_dnumericalorbitpropagator_force_srp_conical_eclipse() {
use crate::propagators::force_model_config::{
ParameterSource, SolarRadiationPressureConfiguration,
};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)),
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: Some(SolarRadiationPressureConfiguration {
area: ParameterSource::Value(10.0),
cr: ParameterSource::Value(1.3),
eclipse_model: EclipseModel::Conical,
}),
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop.step_by(5400.0);
assert!(prop.current_epoch() > epoch);
}
#[test]
fn test_dnumericalorbitpropagator_force_srp_eclipse_transition() {
use crate::propagators::force_model_config::{
ParameterSource, SolarRadiationPressureConfiguration,
};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 3, 20, 0, 0, 0.0, 0.0, TimeSystem::UTC); let state = DVector::from_vec(vec![R_EARTH + 800e3, 0.0, 0.0, 0.0, 7450.0, 0.0]);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)),
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: Some(SolarRadiationPressureConfiguration {
area: ParameterSource::Value(10.0),
cr: ParameterSource::Value(1.3),
eclipse_model: EclipseModel::Conical, }),
third_body: None,
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(R_EARTH + 800e3);
prop.step_by(orbital_period);
assert!(
prop.current_epoch() > epoch,
"Should successfully handle eclipse transitions"
);
assert!(
(prop.current_epoch() - epoch) > (orbital_period - 60.0),
"Should complete nearly full orbit"
);
}
#[test]
fn test_dnumericalorbitpropagator_force_third_body_sun() {
use crate::propagators::force_model_config::ThirdBodyConfiguration;
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 35786e3,
0.0,
0.0, 0.0,
3075.0,
0.0,
]);
let force_config = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: None,
third_body: Some(ThirdBodyConfiguration {
ephemeris_source: EphemerisSource::LowPrecision,
bodies: vec![ThirdBody::Sun],
}),
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop.step_by(86400.0); assert!(prop.current_epoch() > epoch);
}
#[test]
fn test_dnumericalorbitpropagator_force_third_body_moon() {
use crate::propagators::force_model_config::ThirdBodyConfiguration;
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 35786e3, 0.0, 0.0, 0.0, 3075.0, 0.0]);
let force_config = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: None,
third_body: Some(ThirdBodyConfiguration {
ephemeris_source: EphemerisSource::LowPrecision,
bodies: vec![ThirdBody::Moon],
}),
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop.step_by(86400.0);
assert!(prop.current_epoch() > epoch);
}
#[test]
fn test_dnumericalorbitpropagator_force_third_body_planets_de() {
use crate::propagators::force_model_config::ThirdBodyConfiguration;
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 35786e3, 0.0, 0.0, 0.0, 3075.0, 0.0]);
let force_config = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: None,
third_body: Some(ThirdBodyConfiguration {
ephemeris_source: EphemerisSource::DE440s,
bodies: vec![ThirdBody::Sun, ThirdBody::Moon, ThirdBody::Jupiter],
}),
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop.step_by(86400.0);
assert!(prop.current_epoch() > epoch);
}
#[test]
fn test_dnumericalorbitpropagator_force_third_body_jupiter() {
use crate::propagators::force_model_config::ThirdBodyConfiguration;
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 35786e3, 0.0, 0.0, 0.0, 3075.0, 0.0]);
let force_config = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: None,
third_body: Some(ThirdBodyConfiguration {
ephemeris_source: EphemerisSource::DE440s,
bodies: vec![ThirdBody::Jupiter],
}),
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop.step_by(86400.0);
assert!(
prop.current_epoch() > epoch,
"Should successfully propagate with Jupiter perturbation"
);
}
#[test]
fn test_dnumericalorbitpropagator_force_third_body_mars() {
use crate::propagators::force_model_config::ThirdBodyConfiguration;
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 35786e3, 0.0, 0.0, 0.0, 3075.0, 0.0]);
let force_config = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: None,
third_body: Some(ThirdBodyConfiguration {
ephemeris_source: EphemerisSource::DE440s,
bodies: vec![ThirdBody::Mars],
}),
relativity: false,
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop.step_by(86400.0);
assert!(
prop.current_epoch() > epoch,
"Should successfully propagate with Mars perturbation"
);
}
#[test]
fn test_dnumericalorbitpropagator_force_relativity() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let force_config = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: None,
srp: None,
third_body: None,
relativity: true, };
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop.step_by(5400.0);
assert!(prop.current_epoch() > epoch);
}
#[test]
fn test_dnumericalorbitpropagator_force_combined_leo() {
setup_global_test_eop();
setup_global_test_space_weather();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let force_config = ForceModelConfig::leo_default();
let prop_result = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
);
assert!(
prop_result.is_ok(),
"LEO default config should construct successfully"
);
}
#[test]
fn test_dnumericalorbitpropagator_force_combined_geo() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 35786e3, 0.0, 0.0, 0.0, 3075.0, 0.0]);
let force_config = ForceModelConfig::geo_default();
let prop_result = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
);
assert!(
prop_result.is_ok(),
"GEO default config should construct successfully"
);
}
#[test]
fn test_dnumericalorbitpropagator_force_combined_high_fidelity() {
use crate::propagators::force_model_config::{
DragConfiguration, ParameterSource, SolarRadiationPressureConfiguration,
ThirdBodyConfiguration,
};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 800e3, 0.0, 0.0, 0.0, 7450.0, 0.0]);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::Value(1000.0)), frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass, drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::Value(5.0),
cd: ParameterSource::Value(2.2),
}),
srp: Some(SolarRadiationPressureConfiguration {
area: ParameterSource::Value(5.0),
cr: ParameterSource::Value(1.3),
eclipse_model: EclipseModel::Conical,
}),
third_body: Some(ThirdBodyConfiguration {
ephemeris_source: EphemerisSource::LowPrecision,
bodies: vec![ThirdBody::Sun, ThirdBody::Moon],
}),
relativity: true, };
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
let initial_energy = {
let r = prop.current_state().fixed_rows::<3>(0).norm();
let v = prop.current_state().fixed_rows::<3>(3).norm();
0.5 * v * v - GM_EARTH / r
};
prop.step_by(3600.0);
let final_energy = {
let r = prop.current_state().fixed_rows::<3>(0).norm();
let v = prop.current_state().fixed_rows::<3>(3).norm();
0.5 * v * v - GM_EARTH / r
};
assert!(
final_energy < initial_energy,
"Energy should decrease with drag at 800 km altitude"
);
assert!(
prop.current_epoch() > epoch,
"Should successfully propagate with all force models"
);
}
#[test]
fn test_dnumericalorbitpropagator_force_parameter_sensitivity() {
use crate::propagators::force_model_config::{DragConfiguration, ParameterSource};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 300e3, 0.0, 0.0, 0.0, 7700.0, 0.0]);
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)),
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::PointMass,
drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::Value(10.0),
cd: ParameterSource::ParameterIndex(2), }),
srp: None,
third_body: None,
relativity: false,
};
let mut prop1 = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
force_config.clone(),
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
prop1.step_by(600.0);
let final_state1 = prop1.current_state().clone();
let mut prop2 = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
force_config,
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
let params = prop2.current_params();
assert_eq!(params[2], 2.2);
prop2.step_by(600.0);
let final_state2 = prop2.current_state();
let position_diff = (final_state1[0] - final_state2[0]).abs()
+ (final_state1[1] - final_state2[1]).abs()
+ (final_state1[2] - final_state2[2]).abs();
assert!(
position_diff < 100.0,
"Position difference too large: {}",
position_diff
);
}
#[test]
fn test_dnumericalorbitpropagator_propagation_mode_state_only() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert!(prop.stm().is_none());
assert!(prop.sensitivity().is_none());
}
#[test]
fn test_dnumericalorbitpropagator_propagation_mode_with_stm() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert!(prop.stm().is_some());
assert!(prop.sensitivity().is_none());
prop.step_by(100.0);
let stm = prop.stm().unwrap();
assert_eq!(stm.nrows(), 6);
assert_eq!(stm.ncols(), 6);
let identity = DMatrix::<f64>::identity(6, 6);
let mut differs = false;
for i in 0..6 {
for j in 0..6 {
let diff: f64 = stm[(i, j)] - identity[(i, j)];
if diff.abs() > 1e-6 {
differs = true;
break;
}
}
}
assert!(differs, "STM should have evolved from identity");
}
#[test]
fn test_dnumericalorbitpropagator_propagation_mode_with_sensitivity() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
config,
test_force_config_with_params(),
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
assert!(prop.stm().is_none());
assert!(prop.sensitivity().is_some());
prop.step_by(100.0);
let sens = prop.sensitivity().unwrap();
assert_eq!(sens.nrows(), 6);
assert_eq!(sens.ncols(), 5);
}
#[test]
fn test_dnumericalorbitpropagator_propagation_mode_with_stm_and_sensitivity() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
config,
test_force_config_with_params(),
Some(default_test_params()),
None,
None,
None,
)
.unwrap();
assert!(prop.stm().is_some());
assert!(prop.sensitivity().is_some());
prop.step_by(100.0);
assert_eq!(prop.stm().unwrap().nrows(), 6);
assert_eq!(prop.sensitivity().unwrap().nrows(), 6);
}
#[test]
fn test_dnumericalorbitpropagator_config_stm() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert!(prop.stm().is_some());
}
#[test]
fn test_dnumericalorbitpropagator_config_sensitivity() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true; let params = DVector::from_vec(vec![1000.0, 10.0, 2.2, 10.0, 1.3]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
config,
ForceModelConfig::earth_gravity(),
Some(params),
None,
None,
None,
)
.unwrap();
assert!(prop.sensitivity().is_some());
}
#[test]
fn test_dnumericalorbitpropagator_config_stm_and_sensitivity() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true; let params = DVector::from_vec(vec![1000.0, 10.0, 2.2, 10.0, 1.3]);
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
config,
ForceModelConfig::earth_gravity(),
Some(params),
None,
None,
None,
)
.unwrap();
assert!(prop.stm().is_some());
assert!(prop.sensitivity().is_some());
}
#[test]
fn test_covariance_propagation_initialization() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut initial_cov = DMatrix::zeros(6, 6);
for i in 0..3 {
initial_cov[(i, i)] = 100.0; }
for i in 3..6 {
initial_cov[(i, i)] = 1.0; }
let prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov.clone()),
)
.unwrap();
assert!(prop.stm().is_some());
assert!(prop.initial_covariance.is_some());
assert!(prop.current_covariance.is_some());
let stored_cov = prop.current_covariance.as_ref().unwrap();
for i in 0..6 {
for j in 0..6 {
assert!((stored_cov[(i, j)] - initial_cov[(i, j)]).abs() < 1e-10);
}
}
}
#[test]
fn test_covariance_propagates_with_stm() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let initial_cov = DMatrix::identity(6, 6) * 10.0;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov.clone()),
)
.unwrap();
prop.propagate_to(epoch + 600.0);
let final_cov = prop.current_covariance.as_ref().unwrap();
let initial_pos_variance = initial_cov[(0, 0)];
let final_pos_variance = final_cov[(0, 0)];
assert!(
final_pos_variance > initial_pos_variance,
"Position variance should grow: initial={}, final={}",
initial_pos_variance,
final_pos_variance
);
}
#[test]
fn test_covariance_stored_in_trajectory() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let initial_cov = DMatrix::identity(6, 6) * 100.0;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(initial_cov),
)
.unwrap();
prop.propagate_steps(5);
assert!(prop.trajectory().covariances.is_some());
let covs = prop.trajectory().covariances.as_ref().unwrap();
assert!(
!covs.is_empty(),
"Trajectory should contain covariance matrices"
);
}
#[test]
fn test_dnumericalorbitpropagator_stm_identity_initial_condition() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let stm = prop.stm().unwrap();
for i in 0..6 {
for j in 0..6 {
let expected = if i == j { 1.0 } else { 0.0 };
assert!(
(stm[(i, j)] - expected).abs() < 1e-10,
"STM[{},{}] = {}, expected {}",
i,
j,
stm[(i, j)],
expected
);
}
}
}
#[test]
fn test_dnumericalorbitpropagator_stm_determinant_preservation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let det_initial = prop.stm().unwrap().determinant();
assert!((det_initial - 1.0).abs() < 1e-10);
let orbital_period = orbital_period(oe[0]);
prop.step_by(orbital_period);
let det_final = prop.stm().unwrap().determinant();
assert!(
(det_final - 1.0).abs() < 1e-6,
"STM determinant should be 1 for Hamiltonian system, got {}",
det_final
);
}
#[test]
fn test_dnumericalorbitpropagator_stm_composition_property() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let mut prop1 = DNumericalOrbitPropagator::new(
epoch,
state_dvec.clone(),
config.clone(),
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let t1 = epoch + 600.0;
prop1.propagate_to(t1);
let phi_t1_t0 = prop1.stm().unwrap().clone();
let state_t1 = prop1.current_state();
let mut prop2 = DNumericalOrbitPropagator::new(
t1,
state_t1,
config.clone(),
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let t2 = t1 + 600.0;
prop2.propagate_to(t2);
let phi_t2_t1 = prop2.stm().unwrap().clone();
let mut prop3 = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop3.propagate_to(t2);
let phi_t2_t0 = prop3.stm().unwrap().clone();
let phi_composed = phi_t2_t1 * phi_t1_t0;
for i in 0..6 {
for j in 0..6 {
let diff = (phi_t2_t0[(i, j)] - phi_composed[(i, j)]).abs();
assert!(
diff < 1e-6,
"STM composition failed at [{},{}]: direct={}, composed={}, diff={}",
i,
j,
phi_t2_t0[(i, j)],
phi_composed[(i, j)],
diff
);
}
}
}
#[test]
fn test_dnumericalorbitpropagator_stm_vs_direct_perturbation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let mut prop_nominal = DNumericalOrbitPropagator::new(
epoch,
state_dvec.clone(),
config.clone(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop_nominal.step_by(100.0);
let stm = prop_nominal.stm().unwrap().clone();
let state_nominal = prop_nominal.current_state();
let delta_x0 = DVector::from_vec(vec![1.0, 0.0, 0.0, 0.0, 0.0, 0.0]); let state_perturbed = state_dvec.clone() + delta_x0.clone();
let mut prop_perturbed = DNumericalOrbitPropagator::new(
epoch,
state_perturbed,
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop_perturbed.step_by(100.0);
let state_perturbed_final = prop_perturbed.current_state();
let delta_x_direct = state_perturbed_final - state_nominal.clone();
let delta_x_stm = stm * delta_x0;
let error = (delta_x_direct - delta_x_stm).norm();
assert!(
error < 0.01,
"STM vs direct perturbation error = {} m, expected < 0.01 m",
error
);
}
#[test]
fn test_dnumericalorbitpropagator_stm_at_methods() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
config.variational.store_stm_history = true;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.propagate_steps(5);
let traj = prop.trajectory();
assert!(traj.len() > 0, "Trajectory should have stored states");
for idx in 0..traj.len() {
let stm_at_idx = prop.stm_at_idx(idx);
assert!(
stm_at_idx.is_some(),
"Should be able to retrieve STM at index {}",
idx
);
let stm = stm_at_idx.unwrap();
assert_eq!(stm.nrows(), 6);
assert_eq!(stm.ncols(), 6);
}
let test_epoch = traj.epochs[traj.len() / 2]; let stm_at_epoch = prop.stm_at(test_epoch);
assert!(
stm_at_epoch.is_some(),
"Should retrieve STM at stored epoch"
);
let stm_by_idx = prop.stm_at_idx(traj.len() / 2).unwrap();
let stm_by_epoch = stm_at_epoch.unwrap();
let diff = (stm_by_idx - stm_by_epoch).norm();
assert!(
diff < 1e-12,
"STM retrieval by epoch and index should match"
);
}
#[test]
fn test_dnumericalorbitpropagator_stm_with_different_jacobian_methods() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let mut config_forward = NumericalPropagationConfig::default();
config_forward.variational.enable_stm = true;
config_forward.variational.jacobian_method = DifferenceMethod::Forward;
let mut prop_forward = DNumericalOrbitPropagator::new(
epoch,
state_dvec.clone(),
config_forward,
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop_forward.step_by(100.0);
let stm_forward = prop_forward.stm().unwrap().clone();
let mut config_central = NumericalPropagationConfig::default();
config_central.variational.enable_stm = true;
config_central.variational.jacobian_method = DifferenceMethod::Central;
let mut prop_central = DNumericalOrbitPropagator::new(
epoch,
state_dvec.clone(),
config_central,
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop_central.step_by(100.0);
let stm_central = prop_central.stm().unwrap().clone();
let mut config_backward = NumericalPropagationConfig::default();
config_backward.variational.enable_stm = true;
config_backward.variational.jacobian_method = DifferenceMethod::Backward;
let mut prop_backward = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config_backward,
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop_backward.step_by(100.0);
let stm_backward = prop_backward.stm().unwrap().clone();
for i in 0..6 {
for j in 0..6 {
let avg = (stm_forward[(i, j)].abs()
+ stm_central[(i, j)].abs()
+ stm_backward[(i, j)].abs())
/ 3.0;
if avg > 1e-6 {
let diff_fc = (stm_forward[(i, j)] - stm_central[(i, j)]).abs();
let diff_cb = (stm_central[(i, j)] - stm_backward[(i, j)]).abs();
assert!(
diff_fc / avg < 0.01,
"Forward vs Central mismatch at [{},{}]",
i,
j
);
assert!(
diff_cb / avg < 0.01,
"Central vs Backward mismatch at [{},{}]",
i,
j
);
}
}
}
}
#[test]
fn test_dnumericalorbitpropagator_stm_eigenvalue_analysis() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe[0]);
prop.step_by(orbital_period / 2.0);
let stm = prop.stm().unwrap().clone();
let eigenvalues = stm.complex_eigenvalues();
for (i, eigenvalue) in eigenvalues.iter().enumerate() {
let magnitude = (eigenvalue.re * eigenvalue.re + eigenvalue.im * eigenvalue.im).sqrt();
assert!(
(magnitude - 1.0).abs() < 0.1,
"Eigenvalue {} has magnitude {}, expected ≈ 1 (neutral stability)",
i,
magnitude
);
}
}
#[test]
fn test_dnumericalorbitpropagator_stm_with_different_force_models() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let mut prop_gravity = DNumericalOrbitPropagator::new(
epoch,
state_dvec.clone(),
config.clone(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop_gravity.step_by(300.0);
let stm_gravity = prop_gravity.stm().unwrap();
let det_gravity = stm_gravity.determinant();
assert!((det_gravity.abs() - 1.0).abs() < 1e-6);
use crate::orbit_dynamics::gravity::GravityModelType;
let force_config_j2 = ForceModelConfig {
mass: None,
frame_transform: FrameTransformationModel::default(),
gravity: GravityConfiguration::SphericalHarmonic {
source: GravityModelSource::ModelType(GravityModelType::EGM2008_360),
degree: 2,
order: 0,
},
drag: None,
srp: None,
third_body: None,
relativity: false,
};
let mut prop_j2 = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
force_config_j2,
None,
None,
None,
None,
)
.unwrap();
prop_j2.step_by(300.0);
let stm_j2 = prop_j2.stm().unwrap();
let stm_diff = (stm_gravity - stm_j2).norm();
assert!(
stm_diff > 1e-6,
"STM should differ between gravity models, diff = {}",
stm_diff
);
}
#[test]
#[cfg_attr(not(feature = "manual"), ignore)] fn test_dnumericalorbitpropagator_stm_accuracy_degradation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
ForceModelConfig::two_body_gravity(),
None,
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe[0]);
prop.step_by(orbital_period);
let det_1_orbit = prop.stm().unwrap().determinant();
prop.step_by(9.0 * orbital_period);
let det_10_orbits = prop.stm().unwrap().determinant();
assert!(
(det_1_orbit - 1.0).abs() < 1e-6,
"STM determinant after 1 orbit = {}",
det_1_orbit
);
assert!(
(det_10_orbits - 1.0).abs() < 1e-4,
"STM determinant after 10 orbits = {} (some degradation expected)",
det_10_orbits
);
}
#[test]
fn test_dnumericalorbitpropagator_stm_interpolation_accuracy() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
config.variational.store_stm_history = true;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec.clone(),
config.clone(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
for _ in 0..10 {
prop.step_by(60.0);
}
let mid_epoch = epoch + 300.0; let stm_interpolated = prop.stm_at(mid_epoch);
assert!(stm_interpolated.is_some());
let mut prop_exact = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop_exact.propagate_to(mid_epoch);
let stm_exact = prop_exact.stm().unwrap();
let stm_interp = stm_interpolated.unwrap();
let error = (stm_interp - stm_exact).norm();
let stm_norm = stm_exact.norm();
assert!(
error / stm_norm < 0.01,
"STM interpolation relative error = {}, expected < 1%",
error / stm_norm
);
}
#[test]
fn test_dnumericalorbitpropagator_sensitivity_vs_finite_difference() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let params_nominal = DVector::from_vec(vec![500.0, 2.2, 1.5, 10.0, 0.3]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true;
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)), frame_transform: FrameTransformationModel::default(),
drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::Value(10.0),
cd: ParameterSource::Value(2.2),
}),
srp: None,
third_body: None,
..Default::default()
};
let mut prop_sens = DNumericalOrbitPropagator::new(
epoch,
state_dvec.clone(),
config,
force_config.clone(),
Some(params_nominal.clone()),
None,
None,
None,
)
.unwrap();
prop_sens.step_by(300.0); let sensitivity = prop_sens.sensitivity().unwrap().clone();
let state_nominal = prop_sens.current_state();
let delta_p = 1.0; let mut params_perturbed = params_nominal.clone();
params_perturbed[0] += delta_p;
let mut prop_perturbed = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
NumericalPropagationConfig::default(),
force_config,
Some(params_perturbed),
None,
None,
None,
)
.unwrap();
prop_perturbed.step_by(300.0);
let state_perturbed = prop_perturbed.current_state();
let sensitivity_fd = (state_perturbed - state_nominal) / delta_p;
for i in 0..6 {
let rel_error = if sensitivity_fd[i].abs() > 1e-6 {
((sensitivity[(i, 0)] - sensitivity_fd[i]) / sensitivity_fd[i]).abs()
} else {
(sensitivity[(i, 0)] - sensitivity_fd[i]).abs()
};
assert!(
rel_error < 0.10, "Sensitivity[{}] mismatch: variational={}, FD={}, rel_error={}",
i,
sensitivity[(i, 0)],
sensitivity_fd[i],
rel_error
);
}
}
#[test]
fn test_dnumericalorbitpropagator_sensitivity_mass_physical_reasonableness() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(
R_EARTH + 400e3, 0.01,
97.8_f64,
0.0,
0.0,
0.0,
);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let params = DVector::from_vec(vec![500.0, 2.2, 1.5, 10.0, 0.3]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true;
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)), frame_transform: FrameTransformationModel::default(),
drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::Value(10.0),
cd: ParameterSource::Value(2.2),
}),
srp: None,
third_body: None,
..Default::default()
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
force_config,
Some(params),
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe[0]);
prop.step_by(orbital_period);
let sensitivity = prop.sensitivity().unwrap();
let dv_dmass = sensitivity.fixed_rows::<3>(3);
let dv_dmass_mag = dv_dmass.norm();
assert!(
dv_dmass_mag > 1e-6,
"Mass sensitivity should be non-negligible for drag, got {}",
dv_dmass_mag
);
}
#[test]
fn test_dnumericalorbitpropagator_sensitivity_drag_coefficient() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 400e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let params = DVector::from_vec(vec![500.0, 2.2, 1.5, 10.0, 0.3]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true;
let force_config = ForceModelConfig {
mass: Some(ParameterSource::Value(500.0)), frame_transform: FrameTransformationModel::default(),
drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::Value(10.0),
cd: ParameterSource::ParameterIndex(1), }),
srp: None,
third_body: None,
..Default::default()
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
force_config,
Some(params),
None,
None,
None,
)
.unwrap();
let orbital_period = orbital_period(oe[0]);
prop.step_by(orbital_period);
let sensitivity = prop.sensitivity().unwrap();
let sensitivity_norm = sensitivity.norm();
assert!(
sensitivity_norm > 1e-6,
"Cd sensitivity should be non-negligible, got {}",
sensitivity_norm
);
}
#[test]
#[cfg_attr(not(feature = "manual"), ignore)] fn test_dnumericalorbitpropagator_sensitivity_srp_coefficient() {
use crate::propagators::force_model_config::{
ParameterSource, SolarRadiationPressureConfiguration, ThirdBodyConfiguration,
};
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(
R_EARTH + 35786e3, 0.01,
0.1_f64,
0.0,
0.0,
0.0,
);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let params = DVector::from_vec(vec![1000.0, 2.2, 1.5, 20.0, 0.3]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true;
let force_config = ForceModelConfig {
mass: Some(ParameterSource::Value(1000.0)),
frame_transform: FrameTransformationModel::default(),
drag: None,
srp: Some(SolarRadiationPressureConfiguration {
area: ParameterSource::Value(20.0),
cr: ParameterSource::ParameterIndex(2), eclipse_model: EclipseModel::None,
}),
third_body: Some(ThirdBodyConfiguration {
ephemeris_source: EphemerisSource::LowPrecision,
bodies: vec![ThirdBody::Sun],
}),
..Default::default()
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
force_config,
Some(params),
None,
None,
None,
)
.unwrap();
prop.step_by(86400.0);
let sensitivity = prop.sensitivity().unwrap();
let sensitivity_norm = sensitivity.norm();
assert!(
sensitivity_norm > 1e-6,
"Cr sensitivity should be non-negligible at GEO, got {}",
sensitivity_norm
);
}
#[test]
fn test_dnumericalorbitpropagator_sensitivity_zero_for_unused_parameters() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let params = DVector::from_vec(vec![500.0, 2.2, 1.5]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true;
let force_config = ForceModelConfig::earth_gravity();
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
force_config,
Some(params),
None,
None,
None,
)
.unwrap();
prop.step_by(300.0);
let sensitivity = prop.sensitivity().unwrap();
let sensitivity_norm = sensitivity.norm();
assert!(
sensitivity_norm < 1e-10,
"Sensitivity for unused parameter should be ~0, got {}",
sensitivity_norm
);
}
#[test]
fn test_dnumericalorbitpropagator_sensitivity_storage_in_trajectory() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let params = DVector::from_vec(vec![500.0, 2.2]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true;
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)), frame_transform: FrameTransformationModel::default(),
drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::Value(10.0),
cd: ParameterSource::ParameterIndex(1),
}),
srp: None,
third_body: None,
..Default::default()
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
force_config,
Some(params),
None,
None,
None,
)
.unwrap();
prop.propagate_steps(5);
assert!(
prop.trajectory().sensitivities.is_some(),
"Trajectory should contain sensitivity matrices"
);
let sensitivities = prop.trajectory().sensitivities.as_ref().unwrap();
assert!(
!sensitivities.is_empty(),
"Sensitivity storage should not be empty"
);
}
#[test]
fn test_dnumericalorbitpropagator_sensitivity_at_methods() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let oe = nalgebra::Vector6::new(R_EARTH + 500e3, 0.01, 97.8_f64, 0.0, 0.0, 0.0);
let state = state_koe_to_eci(oe, AngleFormat::Degrees);
let state_dvec = DVector::from_vec(state.as_slice().to_vec());
let params = DVector::from_vec(vec![500.0]);
let mut config = NumericalPropagationConfig::default();
config.variational.enable_sensitivity = true;
config.variational.store_sensitivity_history = true;
let force_config = ForceModelConfig {
mass: Some(ParameterSource::ParameterIndex(0)), frame_transform: FrameTransformationModel::default(),
drag: Some(DragConfiguration {
model: AtmosphericModel::HarrisPriester,
area: ParameterSource::Value(10.0),
cd: ParameterSource::Value(2.2),
}),
srp: None,
third_body: None,
..Default::default()
};
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
config,
force_config,
Some(params),
None,
None,
None,
)
.unwrap();
prop.propagate_steps(5);
let traj = prop.trajectory();
assert!(traj.len() > 0, "Trajectory should have stored states");
for idx in 0..traj.len() {
let sens_at_idx = prop.sensitivity_at_idx(idx);
assert!(
sens_at_idx.is_some(),
"Should be able to retrieve sensitivity at index {}",
idx
);
let sens = sens_at_idx.unwrap();
assert_eq!(sens.nrows(), 6);
assert_eq!(sens.ncols(), 1); }
let test_epoch = traj.epochs[traj.len() / 2]; let sens_at_epoch = prop.sensitivity_at(test_epoch);
assert!(
sens_at_epoch.is_some(),
"Should retrieve sensitivity at stored epoch"
);
let sens_by_idx = prop.sensitivity_at_idx(traj.len() / 2).unwrap();
let sens_by_epoch = sens_at_epoch.unwrap();
let diff = (sens_by_idx - sens_by_epoch).norm();
assert!(
diff < 1e-12,
"Sensitivity retrieval by epoch and index should match"
);
let sens_final = prop.sensitivity_at_idx(traj.len() - 1).unwrap();
assert!(
sens_final.norm() > 1e-6,
"Final sensitivity should be non-zero"
);
}
#[test]
fn test_dnumericalorbitpropagator_covariance_stm_formula_verification() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7.6e3, 0.0]);
let mut p0 = DMatrix::from_element(6, 6, 0.0);
for i in 0..3 {
p0[(i, i)] = 100.0 * 100.0;
p0[(i + 3, i + 3)] = 0.1 * 0.1;
}
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(p0.clone()),
)
.unwrap();
prop.step_by(300.0);
let phi = prop.stm().unwrap();
let p_propagated = prop.current_covariance.as_ref().unwrap();
let p_computed = phi * &p0 * phi.transpose();
let diff = (p_propagated - p_computed).norm();
let rel_error = diff / p_propagated.norm();
assert!(
rel_error < 1e-10,
"Covariance should match P(t) = Φ·P₀·Φᵀ formula, rel_error: {}",
rel_error
);
}
#[test]
fn test_dnumericalorbitpropagator_covariance_monte_carlo_validation() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, TimeSystem::UTC);
let state_nominal = vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7.6e3, 0.0];
let mut p0 = DMatrix::from_element(6, 6, 0.0);
for i in 0..3 {
p0[(i, i)] = 100.0 * 100.0; p0[(i + 3, i + 3)] = 0.1 * 0.1; }
let state_dvec = DVector::from_vec(state_nominal.clone());
let mut prop_cov = DNumericalOrbitPropagator::new(
epoch,
state_dvec,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(p0.clone()),
)
.unwrap();
let dt = 300.0; prop_cov.step_by(dt);
let p_propagated = prop_cov.current_covariance.as_ref().unwrap();
use rand::SeedableRng;
use rv::dist::Gaussian;
use rv::prelude::Sampleable;
let n_samples = 100;
let mut final_states = Vec::new();
let mut rng = rand::rngs::StdRng::seed_from_u64(12345);
for _ in 0..n_samples {
let mut state_perturbed = state_nominal.clone();
for i in 0..6 {
let std_dev = p0[(i, i)].sqrt();
let gaussian = Gaussian::new(0.0, std_dev).unwrap();
let sample: f64 = gaussian.draw(&mut rng);
state_perturbed[i] += sample;
}
let state_dvec_pert = DVector::from_vec(state_perturbed);
let mut prop_mc = DNumericalOrbitPropagator::new(
epoch,
state_dvec_pert,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop_mc.step_by(dt);
final_states.push(prop_mc.current_state_ref().clone());
}
let mut mean_state = DVector::from_element(6, 0.0);
for state in &final_states {
mean_state += state;
}
mean_state /= n_samples as f64;
let mut p_mc = DMatrix::from_element(6, 6, 0.0);
for state in &final_states {
let delta = state - &mean_state;
p_mc += &delta * delta.transpose();
}
p_mc /= (n_samples - 1) as f64;
for i in 0..6 {
let rel_error = (p_propagated[(i, i)] - p_mc[(i, i)]).abs() / p_propagated[(i, i)];
assert!(
rel_error < 0.5, "Variance {} should match Monte Carlo, rel_error: {}",
i,
rel_error
);
}
}
#[test]
fn test_dnumericalorbitpropagator_events_by_detector_index() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 1800.0, "Event A"))); prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 2700.0, "Event B"))); prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 3600.0, "Event C")));
prop.propagate_to(epoch + 7200.0);
let events_0 = prop.events_by_detector_index(0);
assert_eq!(events_0.len(), 1);
assert_eq!(events_0[0].name, "Event A");
assert_eq!(events_0[0].detector_index, 0);
let events_1 = prop.events_by_detector_index(1);
assert_eq!(events_1.len(), 1);
assert_eq!(events_1[0].name, "Event B");
assert_eq!(events_1[0].detector_index, 1);
let events_2 = prop.events_by_detector_index(2);
assert_eq!(events_2.len(), 1);
assert_eq!(events_2[0].name, "Event C");
assert_eq!(events_2[0].detector_index, 2);
let events_99 = prop.events_by_detector_index(99);
assert_eq!(events_99.len(), 0);
}
#[test]
fn test_dnumericalorbitpropagator_events_combined_filters() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 1000.0, "Early Event")));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 2000.0, "Mid Event")));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 3000.0, "Late Event")));
prop.propagate_to(epoch + 4000.0);
let events_in_range = prop.events_by_detector_index_in_range(0, epoch, epoch + 1500.0);
assert_eq!(events_in_range.len(), 1);
assert_eq!(events_in_range[0].name, "Early Event");
let events_by_name = prop.events_by_name_in_range("Event", epoch, epoch + 2500.0);
assert_eq!(events_by_name.len(), 2);
let events_filtered: Vec<_> = prop
.query_events()
.by_detector_index(1)
.in_time_range(epoch, epoch + 2500.0)
.collect();
assert_eq!(events_filtered.len(), 1);
assert_eq!(events_filtered[0].name, "Mid Event");
}
#[test]
fn test_dnumericalorbitpropagator_query_with_iterator_methods() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
for i in 0..5 {
prop.add_event_detector(Box::new(DTimeEvent::new(
epoch + (i as f64 + 1.0) * 1000.0,
format!("Event {}", i),
)));
}
prop.propagate_to(epoch + 6000.0);
let count = prop.query_events().by_name_contains("Event").count();
assert_eq!(count, 5);
let first_two: Vec<_> = prop
.query_events()
.by_name_contains("Event")
.take(2)
.collect();
assert_eq!(first_two.len(), 2);
let has_event = prop.query_events().by_detector_index(0).next().is_some();
assert!(has_event);
let found = prop.query_events().find(|e| e.name.contains("Event 3"));
assert!(found.is_some());
assert_eq!(found.unwrap().name, "Event 3");
}
#[test]
fn test_dnumericalorbitpropagator_query_edge_cases() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let empty_events = prop.events_by_detector_index(0);
assert_eq!(empty_events.len(), 0);
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 10000.0, "Never Fires")));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 1800.0, "Fires")));
prop.propagate_to(epoch + 3600.0);
let events_0 = prop.events_by_detector_index(0);
assert_eq!(events_0.len(), 0);
let events_1 = prop.events_by_detector_index(1);
assert_eq!(events_1.len(), 1);
assert_eq!(events_1[0].name, "Fires");
assert_eq!(events_1[0].detector_index, 1);
let no_match = prop.events_by_name("NonExistent");
assert_eq!(no_match.len(), 0);
let no_events_in_range = prop.events_in_range(epoch + 5000.0, epoch + 6000.0);
assert_eq!(no_events_in_range.len(), 0);
}
#[test]
fn test_dnumericalorbitpropagator_existing_methods_unchanged() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 1000.0, "Altitude Event")));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 2000.0, "Range Event")));
prop.add_event_detector(Box::new(DTimeEvent::new(epoch + 3000.0, "Altitude Check")));
prop.propagate_to(epoch + 4000.0);
let altitude_events = prop.events_by_name("Altitude");
assert_eq!(altitude_events.len(), 2);
let events_in_range = prop.events_in_range(epoch, epoch + 2500.0);
assert_eq!(events_in_range.len(), 2);
let all_events = prop.event_log();
assert_eq!(all_events.len(), 3);
let latest = prop.latest_event();
assert!(latest.is_some());
assert_eq!(latest.unwrap().name, "Altitude Check");
}
#[test]
fn test_dnumericalorbitpropagator_take_event_detectors() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let detector1 = DTimeEvent::new(epoch + 1000.0, "Event1");
let detector2 = DTimeEvent::new(epoch + 2000.0, "Event2");
prop.add_event_detector(Box::new(detector1));
prop.add_event_detector(Box::new(detector2));
let taken = prop.take_event_detectors();
assert_eq!(taken.len(), 2);
prop.propagate_to(epoch + 3000.0);
assert!(prop.event_log().is_empty());
}
#[test]
fn test_dnumericalorbitpropagator_set_event_detectors() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let detectors: Vec<Box<dyn crate::events::DEventDetector>> = vec![
Box::new(DTimeEvent::new(epoch + 1000.0, "Event1")),
Box::new(DTimeEvent::new(epoch + 2000.0, "Event2")),
];
prop.set_event_detectors(detectors);
prop.propagate_to(epoch + 3000.0);
assert_eq!(prop.event_log().len(), 2);
}
#[test]
fn test_dnumericalorbitpropagator_take_event_log() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let detector = DTimeEvent::new(epoch + 1000.0, "TestEvent");
prop.add_event_detector(Box::new(detector));
prop.propagate_to(epoch + 1500.0);
assert_eq!(prop.event_log().len(), 1);
let taken_log = prop.take_event_log();
assert_eq!(taken_log.len(), 1);
assert_eq!(taken_log[0].name, "TestEvent");
assert!(prop.event_log().is_empty());
}
#[test]
fn test_dnumericalorbitpropagator_set_terminated_is_terminated() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
assert!(!prop.is_terminated());
prop.set_terminated(true);
assert!(prop.is_terminated());
prop.set_terminated(false);
assert!(!prop.is_terminated());
}
#[test]
fn test_dnumericalorbitpropagator_event_detector_roundtrip() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state,
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let detector = DTimeEvent::new(epoch + 1500.0, "RoundtripEvent");
prop.add_event_detector(Box::new(detector));
let taken = prop.take_event_detectors();
assert_eq!(taken.len(), 1);
prop.propagate_to(epoch + 1000.0);
assert!(prop.event_log().is_empty());
prop.set_event_detectors(taken);
prop.propagate_to(epoch + 2000.0);
assert_eq!(prop.event_log().len(), 1);
assert!(prop.event_log()[0].name.contains("RoundtripEvent"));
}
fn iss_gp_record() -> crate::types::GPRecord {
let json = r#"{
"OBJECT_NAME": "ISS (ZARYA)",
"OBJECT_ID": "1998-067A",
"EPOCH": "2024-01-15T12:00:00.000000",
"MEAN_MOTION": 15.50000000,
"ECCENTRICITY": 0.00010000,
"INCLINATION": 51.6400,
"RA_OF_ASC_NODE": 200.0000,
"ARG_OF_PERICENTER": 100.0000,
"MEAN_ANOMALY": 260.0000,
"EPHEMERIS_TYPE": 0,
"CLASSIFICATION_TYPE": "U",
"NORAD_CAT_ID": 25544,
"ELEMENT_SET_NO": 999,
"REV_AT_EPOCH": 45000,
"BSTAR": 0.00034100,
"MEAN_MOTION_DOT": 0.00001000,
"MEAN_MOTION_DDOT": 0.00000000
}"#;
serde_json::from_str(json).unwrap()
}
#[test]
fn test_from_gp_record_two_body() {
setup_global_test_eop();
let record = iss_gp_record();
let prop = DNumericalOrbitPropagator::from_gp_record(
&record,
ForceModelConfig::two_body_gravity(),
None,
None,
);
assert!(prop.is_ok());
let prop = prop.unwrap();
let state = prop.current_state();
let pos_mag = (state[0] * state[0] + state[1] * state[1] + state[2] * state[2]).sqrt();
assert!(
pos_mag > 6000e3 && pos_mag < 7000e3,
"Position magnitude: {}",
pos_mag
);
}
#[test]
fn test_from_gp_record_earth_gravity() {
setup_global_test_eop();
let record = iss_gp_record();
let prop = DNumericalOrbitPropagator::from_gp_record(
&record,
ForceModelConfig::earth_gravity(),
None,
None,
);
assert!(prop.is_ok());
let prop = prop.unwrap();
let state = prop.current_state();
let pos_mag = (state[0] * state[0] + state[1] * state[1] + state[2] * state[2]).sqrt();
assert!(
pos_mag > 6000e3 && pos_mag < 7000e3,
"Position magnitude: {}",
pos_mag
);
}
#[test]
fn test_dnumericalorbitpropagator_reinitialize_epoch_and_state() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
6878.0e3, 0.0, 0.0, 0.0, 7612.0, 0.0, ]);
let config = NumericalPropagationConfig::default();
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop.propagate_to(epoch + 60.0);
let state_at_60 = prop.current_state();
assert_ne!(
state_at_60, state,
"State should have changed after propagation"
);
let new_state = DVector::from_vec(vec![
6900.0e3, 100.0e3, 0.0, 0.0, 7600.0, 0.0,
]);
let new_epoch = epoch + 60.0;
prop.reinitialize(new_epoch, new_state.clone(), None);
assert_eq!(prop.current_epoch(), new_epoch);
assert_eq!(prop.current_state(), new_state);
prop.propagate_to(new_epoch + 60.0);
let state_at_120 = prop.current_state();
assert_ne!(
state_at_120, new_state,
"State should change after propagation from reinitialize"
);
let expected_epoch = epoch + 120.0;
let epoch_diff: f64 = prop.current_epoch() - expected_epoch;
assert!(epoch_diff.abs() < 1e-6);
}
#[test]
fn test_dnumericalorbitpropagator_reinitialize_stm_reset() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![6878.0e3, 0.0, 0.0, 0.0, 7612.0, 0.0]);
let p0 = DMatrix::from_diagonal(&DVector::from_vec(vec![1e6, 1e6, 1e6, 1e2, 1e2, 1e2]));
let mut config = NumericalPropagationConfig::default();
config.variational.enable_stm = true;
let mut prop = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
Some(p0),
)
.unwrap();
let stm_initial = prop.stm().unwrap().clone();
assert_eq!(stm_initial, DMatrix::identity(6, 6));
prop.propagate_to(epoch + 60.0);
let stm_at_60 = prop.stm().unwrap().clone();
assert_ne!(
stm_at_60,
DMatrix::identity(6, 6),
"STM should not be identity after propagation"
);
let new_state = prop.current_state();
prop.reinitialize(epoch + 60.0, new_state, None);
let stm_after_reinit = prop.stm().unwrap().clone();
assert_eq!(
stm_after_reinit,
DMatrix::identity(6, 6),
"STM should be identity after reinitialize"
);
prop.propagate_to(epoch + 120.0);
let stm_at_120 = prop.stm().unwrap().clone();
assert_ne!(
stm_at_120,
DMatrix::identity(6, 6),
"STM should change after propagation from reinitialize"
);
}
#[test]
fn test_dnumericalorbitpropagator_reinitialize_preserves_dynamics() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![6878.0e3, 0.0, 0.0, 0.0, 7612.0, 0.0]);
let config = NumericalPropagationConfig::default();
let mut prop1 = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
config.clone(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop1.propagate_to(epoch + 60.0);
let state_at_60 = prop1.current_state();
prop1.reinitialize(epoch + 60.0, state_at_60.clone(), None);
prop1.propagate_to(epoch + 120.0);
let result_reinit = prop1.current_state();
let mut prop2 = DNumericalOrbitPropagator::new(
epoch + 60.0,
state_at_60,
config,
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
prop2.propagate_to(epoch + 120.0);
let result_fresh = prop2.current_state();
for i in 0..6 {
assert!(
(result_reinit[i] - result_fresh[i]).abs() < 1.0, "Component {}: reinit={}, fresh={}, diff={}",
i,
result_reinit[i],
result_fresh[i],
(result_reinit[i] - result_fresh[i]).abs()
);
}
}
#[test]
fn test_dnumericalorbitpropagator_builder_minimal() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::builder()
.epoch(epoch)
.state(state)
.force_config(ForceModelConfig::earth_gravity())
.build();
assert!(prop.is_ok());
}
#[test]
fn test_dnumericalorbitpropagator_builder_with_propagation_config() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::builder()
.epoch(epoch)
.state(state)
.force_config(ForceModelConfig::earth_gravity())
.propagation_config(NumericalPropagationConfig::high_precision())
.build();
assert!(prop.is_ok());
}
#[test]
fn test_dnumericalorbitpropagator_builder_with_params() {
setup_global_test_eop();
setup_global_test_space_weather();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let prop = DNumericalOrbitPropagator::builder()
.epoch(epoch)
.state(state)
.force_config(ForceModelConfig::leo_default())
.params(default_test_params())
.build();
assert!(prop.is_ok());
}
#[test]
fn test_dnumericalorbitpropagator_builder_with_covariance() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let p0 = DMatrix::<f64>::identity(6, 6) * 1e6;
let prop = DNumericalOrbitPropagator::builder()
.epoch(epoch)
.state(state)
.force_config(ForceModelConfig::earth_gravity())
.initial_covariance(p0)
.build();
assert!(prop.is_ok());
assert!(prop.unwrap().current_covariance().is_some());
}
#[test]
fn test_dnumericalorbitpropagator_builder_full() {
setup_global_test_eop();
setup_global_test_space_weather();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let p0 = DMatrix::<f64>::identity(6, 6) * 1e6;
let prop = DNumericalOrbitPropagator::builder()
.epoch(epoch)
.state(state)
.force_config(ForceModelConfig::leo_default())
.propagation_config(NumericalPropagationConfig::default())
.params(default_test_params())
.initial_covariance(p0)
.build();
assert!(prop.is_ok());
let prop = prop.unwrap();
assert!(prop.current_covariance().is_some());
assert_eq!(DStatePropagator::state_dim(&prop), 6);
}
#[test]
fn test_dnumericalorbitpropagator_builder_with_additional_dynamics_and_control_input() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![
R_EARTH + 500e3,
0.0,
0.0,
0.0,
7500.0,
0.0,
1000.0, ]);
let additional_dynamics: DStateDynamics = Box::new(|_t, state, _params| {
let mut dx = DVector::zeros(state.len());
dx[6] = -0.1; dx
});
let control: Box<
dyn Fn(f64, &DVector<f64>, Option<&DVector<f64>>) -> DVector<f64> + Send + Sync,
> = Box::new(|_t, state, _params| {
let v = Vector3::new(state[3], state[4], state[5]);
let v_mag = v.norm();
let mut dx = DVector::zeros(state.len());
if v_mag > 1e-6 {
let a = v * (0.0001 / v_mag);
dx[3] = a[0];
dx[4] = a[1];
dx[5] = a[2];
}
dx
});
let prop = DNumericalOrbitPropagator::builder()
.epoch(epoch)
.state(state)
.force_config(ForceModelConfig::earth_gravity())
.additional_dynamics(additional_dynamics)
.control_input(control)
.build();
assert!(prop.is_ok());
let mut prop = prop.unwrap();
assert_eq!(DStatePropagator::state_dim(&prop), 7);
let initial_mass = prop.current_state()[6];
prop.step_by(10.0);
assert!((prop.current_state()[6] - (initial_mass - 1.0)).abs() < 1e-3);
}
#[test]
fn test_dnumericalorbitpropagator_builder_matches_new() {
setup_global_test_eop();
let epoch = Epoch::from_datetime(2024, 1, 1, 0, 0, 0.0, 0.0, TimeSystem::UTC);
let state = DVector::from_vec(vec![R_EARTH + 500e3, 0.0, 0.0, 0.0, 7500.0, 0.0]);
let via_new = DNumericalOrbitPropagator::new(
epoch,
state.clone(),
NumericalPropagationConfig::default(),
ForceModelConfig::earth_gravity(),
None,
None,
None,
None,
)
.unwrap();
let via_builder = DNumericalOrbitPropagator::builder()
.epoch(epoch)
.state(state)
.force_config(ForceModelConfig::earth_gravity())
.build()
.unwrap();
assert_eq!(via_new.initial_epoch(), via_builder.initial_epoch());
assert_eq!(
DStatePropagator::state_dim(&via_new),
DStatePropagator::state_dim(&via_builder)
);
}
}