Skip to main content

nyx_space/md/trajectory/
mod.rs

1/*
2    Nyx, blazing fast astrodynamics
3    Copyright (C) 2018-onwards Christopher Rabotin <christopher.rabotin@gmail.com>
4
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU Affero General Public License as published
7    by the Free Software Foundation, either version 3 of the License, or
8    (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU Affero General Public License for more details.
14
15    You should have received a copy of the GNU Affero General Public License
16    along with this program.  If not, see <https://www.gnu.org/licenses/>.
17*/
18
19use anise::{
20    errors::PhysicsError,
21    math::{cartesian::CartesianState, interpolation::InterpolationError},
22};
23use snafu::prelude::*;
24
25mod interpolatable;
26mod sc_traj;
27mod traj;
28mod traj_it;
29
30pub(crate) use interpolatable::INTERPOLATION_SAMPLES;
31pub use interpolatable::Interpolatable;
32pub use traj::Traj;
33
34pub use crate::io::ExportCfg;
35
36use super::StateParameter;
37use crate::time::{Duration, Epoch};
38
39#[derive(Clone, PartialEq, Debug, Snafu)]
40pub enum TrajError {
41    #[snafu(display("Event {event} not found between {start} and {end}"))]
42    EventNotFound {
43        start: Epoch,
44        end: Epoch,
45        event: String,
46    },
47    #[snafu(display("No interpolation data at {epoch}"))]
48    NoInterpolationData {
49        epoch: Epoch,
50    },
51    #[snafu(display("Failed to create trajectory: {msg}"))]
52    CreationError {
53        msg: String,
54    },
55    #[snafu(display(
56        "Probable bug: Requested epoch {req_epoch}, corresponding to an offset of {req_dur} in a spline of duration {spline_dur}"
57    ))]
58    OutOfSpline {
59        req_epoch: Epoch,
60        req_dur: Duration,
61        spline_dur: Duration,
62    },
63    #[snafu(display("Interpolation failed: {source}"))]
64    Interpolation {
65        source: InterpolationError,
66    },
67    TrajPhysics {
68        source: PhysicsError,
69    },
70    TrajGeneric {
71        err: String,
72    },
73}
74
75/// Smooth the RIC differences using an in-line median filter.
76/// This avoids allocations and operates directly on the Cartesian components.
77fn smooth_state_diff_in_place(ric_diff: &mut [CartesianState], window_size: usize) {
78    assert!(
79        window_size % 2 == 1,
80        "Window size must be odd for proper median calculation"
81    );
82    let half_window = window_size / 2;
83
84    // Temporary buffer to store sorted values for median calculation
85    let mut temp_buffer = vec![0.0; window_size];
86
87    // Iterate over each state in the array
88    for i in 0..ric_diff.len() {
89        let start = i.saturating_sub(half_window);
90        let end = (i + half_window + 1).min(ric_diff.len());
91
92        // Smooth each component independently
93        for component in 0..6 {
94            // Fill the temporary buffer with values from the current window
95            for (j, idx) in (start..end).enumerate() {
96                temp_buffer[j] = match component {
97                    0 => ric_diff[idx].radius_km.x,
98                    1 => ric_diff[idx].radius_km.y,
99                    2 => ric_diff[idx].radius_km.z,
100                    3 => ric_diff[idx].velocity_km_s.x,
101                    4 => ric_diff[idx].velocity_km_s.y,
102                    5 => ric_diff[idx].velocity_km_s.z,
103                    _ => unreachable!(),
104                };
105            }
106
107            // Sort the buffer to find the median
108            temp_buffer[..end - start].sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
109
110            // Replace the current value with the median
111            let median = temp_buffer[(end - start) / 2];
112            match component {
113                0 => ric_diff[i].radius_km.x = median,
114                1 => ric_diff[i].radius_km.y = median,
115                2 => ric_diff[i].radius_km.z = median,
116                3 => ric_diff[i].velocity_km_s.x = median,
117                4 => ric_diff[i].velocity_km_s.y = median,
118                5 => ric_diff[i].velocity_km_s.z = median,
119                _ => unreachable!(),
120            }
121        }
122    }
123}