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}, traj valid from {start} to {end}"))]
48    NoInterpolationData {
49        epoch: Epoch,
50        start: Epoch,
51        end: Epoch,
52    },
53    #[snafu(display("cannot interpolate empty trajectory at {epoch}"))]
54    EmptyTrajectory {
55        epoch: Epoch,
56    },
57    #[snafu(display("Failed to create trajectory: {msg}"))]
58    CreationError {
59        msg: String,
60    },
61    #[snafu(display(
62        "Probable bug: Requested epoch {req_epoch}, corresponding to an offset of {req_dur} in a spline of duration {spline_dur}"
63    ))]
64    OutOfSpline {
65        req_epoch: Epoch,
66        req_dur: Duration,
67        spline_dur: Duration,
68    },
69    #[snafu(display("Interpolation failed: {source}"))]
70    Interpolation {
71        source: InterpolationError,
72    },
73    TrajPhysics {
74        source: PhysicsError,
75    },
76    TrajGeneric {
77        err: String,
78    },
79}
80
81/// Smooth the RIC differences using an in-line median filter.
82/// This avoids allocations and operates directly on the Cartesian components.
83fn smooth_state_diff_in_place(ric_diff: &mut [CartesianState], window_size: usize) {
84    assert!(
85        window_size % 2 == 1,
86        "Window size must be odd for proper median calculation"
87    );
88    let half_window = window_size / 2;
89
90    // Temporary buffer to store sorted values for median calculation
91    let mut temp_buffer = vec![0.0; window_size];
92
93    // Iterate over each state in the array
94    for i in 0..ric_diff.len() {
95        let start = i.saturating_sub(half_window);
96        let end = (i + half_window + 1).min(ric_diff.len());
97
98        // Smooth each component independently
99        for component in 0..6 {
100            // Fill the temporary buffer with values from the current window
101            for (j, idx) in (start..end).enumerate() {
102                temp_buffer[j] = match component {
103                    0 => ric_diff[idx].radius_km.x,
104                    1 => ric_diff[idx].radius_km.y,
105                    2 => ric_diff[idx].radius_km.z,
106                    3 => ric_diff[idx].velocity_km_s.x,
107                    4 => ric_diff[idx].velocity_km_s.y,
108                    5 => ric_diff[idx].velocity_km_s.z,
109                    _ => unreachable!(),
110                };
111            }
112
113            // Sort the buffer to find the median
114            temp_buffer[..end - start].sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
115
116            // Replace the current value with the median
117            let median = temp_buffer[(end - start) / 2];
118            match component {
119                0 => ric_diff[i].radius_km.x = median,
120                1 => ric_diff[i].radius_km.y = median,
121                2 => ric_diff[i].radius_km.z = median,
122                3 => ric_diff[i].velocity_km_s.x = median,
123                4 => ric_diff[i].velocity_km_s.y = median,
124                5 => ric_diff[i].velocity_km_s.z = median,
125                _ => unreachable!(),
126            }
127        }
128    }
129}