nyx_space/tools/lambert/
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 crate::cosmic::AstroError;
20use crate::errors::LambertError;
21use crate::linalg::Vector3;
22use std::f64::consts::PI;
23
24mod godding;
25mod izzo;
26
27use anise::errors::PhysicsError;
28use anise::prelude::Orbit;
29pub use godding::gooding;
30pub use izzo::izzo;
31
32const TAU: f64 = 2.0 * PI;
33const LAMBERT_EPSILON: f64 = 1e-4; // General epsilon
34const LAMBERT_EPSILON_TIME: f64 = 1e-4; // Time epsilon
35const LAMBERT_EPSILON_RAD: f64 = (5e-5 / 180.0) * PI; // 0.00005 degrees
36/// Maximum number of iterations allowed in the Lambert problem solver.
37/// This is a safety measure to prevent infinite loops in case a solution cannot be found.
38const MAX_ITERATIONS: usize = 1000;
39
40/// Define the transfer kind for a Lambert
41pub enum TransferKind {
42    Auto,
43    ShortWay,
44    LongWay,
45    NRevs(u8),
46}
47
48impl TransferKind {
49    /// Calculate the direction multiplier based on the transfer kind.
50    ///
51    /// # Arguments
52    ///
53    /// * `r_final` - The final radius vector.
54    /// * `r_init` - The initial radius vector.
55    ///
56    /// # Returns
57    ///
58    /// * `Result<f64, NyxError>` - The direction multiplier or an error if the transfer kind is not supported.
59    fn direction_of_motion(
60        self,
61        r_final: &Vector3<f64>,
62        r_init: &Vector3<f64>,
63    ) -> Result<f64, LambertError> {
64        match self {
65            TransferKind::Auto => {
66                let mut dnu = r_final[1].atan2(r_final[0]) - r_init[1].atan2(r_final[1]);
67                if dnu > TAU {
68                    dnu -= TAU;
69                } else if dnu < 0.0 {
70                    dnu += TAU;
71                }
72
73                if dnu > std::f64::consts::PI {
74                    Ok(-1.0)
75                } else {
76                    Ok(1.0)
77                }
78            }
79            TransferKind::ShortWay => Ok(1.0),
80            TransferKind::LongWay => Ok(-1.0),
81            _ => Err(LambertError::MultiRevNotSupported),
82        }
83    }
84}
85
86#[derive(Copy, Clone, Debug)]
87pub struct LambertInput {
88    pub initial_state: Orbit,
89    pub final_state: Orbit,
90}
91
92impl LambertInput {
93    pub fn from_planetary_states(
94        initial_state: Orbit,
95        final_state: Orbit,
96    ) -> Result<Self, AstroError> {
97        if final_state.frame != initial_state.frame {
98            return Err(AstroError::AstroPhysics {
99                source: PhysicsError::FrameMismatch {
100                    action: "Lambert solver requires both states to be in the same frame",
101                    frame1: final_state.frame.into(),
102                    frame2: initial_state.frame.into(),
103                },
104            });
105        }
106        // Ensure that the GM is set
107        initial_state
108            .frame
109            .mu_km3_s2()
110            .map_err(|e| AstroError::AstroPhysics { source: e })?;
111
112        Ok(Self {
113            initial_state,
114            final_state,
115        })
116    }
117
118    /// Return the gravitational parameter of this Lambert problem
119    pub fn mu_km2_s3(&self) -> f64 {
120        self.initial_state.frame.mu_km3_s2().unwrap()
121    }
122}
123
124/// A solution to the Lambert problem.
125#[derive(Copy, Clone, Debug)]
126pub struct LambertSolution {
127    pub v_init_km_s: Vector3<f64>,
128    pub v_final_km_s: Vector3<f64>,
129    /// Turn angle ONLY computed with Godding method
130    pub phi_rad: f64,
131    pub input: LambertInput,
132}
133
134impl LambertSolution {
135    /// Return the v infinity vector at departure, in km/s
136    pub fn v_inf_depart_km_s(&self) -> Vector3<f64> {
137        self.input.initial_state.velocity_km_s - self.v_init_km_s
138    }
139
140    /// Return v infinity vector at arrival, in km/s
141    pub fn v_inf_arrive_km_s(&self) -> Vector3<f64> {
142        self.input.final_state.velocity_km_s - self.v_final_km_s
143    }
144
145    /// Return the transfer orbit computed from setting the departure velocity to the initial state.
146    pub fn transfer_orbit(mut self) -> Orbit {
147        self.input.initial_state.velocity_km_s = self.v_init_km_s;
148        self.input.initial_state
149    }
150
151    /// Return the arrival orbit computed from setting the arrival velocity to the final state.
152    pub fn arrival_orbit(mut self) -> Orbit {
153        self.input.final_state.velocity_km_s = self.v_final_km_s;
154        self.input.final_state
155    }
156
157    /// Return the declination of the departure v infinity, in degrees
158    pub fn v_inf_depart_declination_deg(&self) -> f64 {
159        let v_inf_km_s = self.v_inf_depart_km_s();
160        (v_inf_km_s.z / v_inf_km_s.norm()).asin().to_degrees()
161    }
162
163    /// Returns the c3 computed as the departure v infinity norm squared
164    pub fn c3_km2_s2(&self) -> f64 {
165        self.v_inf_depart_km_s().norm_squared()
166    }
167}
nyx_space/tools/lambert/mod.rs

nyx_space/tools/lambert/
mod.rs
