nyx_space/mc/
helpers.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 super::{Distribution, Normal, Rng, Uniform};
20use crate::linalg::Vector3;
21use crate::NyxError;
22
23/// Returns a unit vector from a normal distribution.
24/// Implements the Sphere Point Picking method: https://mathworld.wolfram.com/SpherePointPicking.html
25pub fn unit_vector_from_seed<R: Rng>(rng: &mut R) -> Vector3<f64> {
26    let distr = Uniform::new_inclusive(0.0, 1.0);
27    let u = distr.sample(rng);
28    let v = distr.sample(rng);
29    let theta = std::f64::consts::TAU * u;
30    let phi = (2.0 * v - 1.0).acos();
31    Vector3::new(theta.cos() * phi.sin(), theta.sin() * phi.sin(), phi.cos())
32}
33
34/// Apply a Normal pointing error to the provided delta-v vector (must be in km/s).
35/// Requires a 3-sigma error percentage (e.g. 0.05 is a 5% pointing error) and a pseudo random number generator
36/// Returns the delta-v vector with the applying pointing error
37pub fn dv_pointing_error<R: Rng>(
38    cur_pointing: &Vector3<f64>,
39    dv: Vector3<f64>,
40    error_prct3s: f64,
41    rng: &mut R,
42) -> Result<Vector3<f64>, NyxError> {
43    if !(0.0..1.0).contains(&error_prct3s) {
44        return Err(NyxError::MonteCarlo {
45            msg: format!("Pointing error percentage must be between 0 and 1, got {error_prct3s}"),
46        });
47    }
48
49    let dv_mag = dv.norm();
50    if dv_mag < f64::EPSILON {
51        return Err(NyxError::MonteCarlo {
52            msg: format!("Delta-v vector is nil, cannot apply a pointing error: {dv}"),
53        });
54    }
55
56    let dv_hat = dv / dv_mag;
57    let cur_pointing_mag = cur_pointing.norm();
58    let cur_angle = (cur_pointing.dot(&dv_hat) / cur_pointing_mag).acos();
59
60    let new_angle = Normal::new(cur_angle, error_prct3s / 3.0)
61        .unwrap()
62        .sample(rng);
63    // Return initial delta-v with pointing error
64    Ok(dv_hat * new_angle.cos() * dv_mag)
65}
66
67/// Returns the provided Delta V vector with a magnitude and pointing error
68pub fn dv_execution_error<R: Rng>(
69    cur_pointing: &Vector3<f64>,
70    dv: Vector3<f64>,
71    pointing_3s: f64,
72    mag_3s: f64,
73    rng: &mut R,
74) -> Result<Vector3<f64>, NyxError> {
75    let dv = dv_pointing_error(cur_pointing, dv, pointing_3s, rng)?;
76
77    let new_mag = Normal::new(dv.norm(), mag_3s / 3.0).unwrap().sample(rng);
78
79    Ok(new_mag * (dv / dv.norm()))
80}
81
82#[allow(clippy::float_equality_without_abs)]
83#[test]
84fn test_dv_mag_fixed() {
85    use super::thread_rng;
86    use crate::time::Epoch;
87    use anise::constants::frames::EARTH_J2000;
88    use anise::prelude::Orbit;
89
90    let orbit = Orbit::new(
91        -2436.45,
92        -2436.45,
93        6891.037,
94        5.088_611,
95        -5.088_611,
96        0.0,
97        Epoch::from_gregorian_tai_at_noon(2021, 3, 24),
98        EARTH_J2000,
99    );
100
101    let dv_mag_distr = Normal::new(5e-3, 5e-4).unwrap();
102
103    for _ in 0..=1000 {
104        let dv_mag = dv_mag_distr.sample(&mut thread_rng());
105        let dv_point = unit_vector_from_seed(&mut thread_rng());
106        let dv = dv_point * dv_mag;
107        let dv_w_err = dv_pointing_error(&orbit.velocity_km_s, dv, 0.1, &mut thread_rng()).unwrap();
108        assert!(
109            (dv_w_err.norm() - dv_mag) < f64::EPSILON,
110            "{:.1e}",
111            (dv_w_err.norm() - dv_mag)
112        );
113    }
114}
nyx_space/mc/helpers.rs

nyx_space/mc/
helpers.rs
