keplerian_orbits/
lib.rs

1#![allow(non_snake_case)]
2
3use std::f32::consts::PI;
4
5use glam::{Mat3, Quat, Vec3};
6
7pub struct Orbit {
8    pub eccentricity: f32,
9    pub semi_major_axis: f32,
10    pub inclination: f32,
11    pub longitude_of_ascending_node: f32,
12    pub argument_of_periapsis: f32,
13}
14
15pub struct StateVectors {
16    pub position: Vec3,
17    pub velocity: Vec3,
18}
19
20impl StateVectors {
21    pub fn new(position: Vec3, velocity: Vec3) -> Self {
22        Self { position, velocity }
23    }
24}
25
26/// Gravitational constant
27const G: f32 = 6.67430e-11;
28
29impl Orbit {
30    pub fn state_vectors_to_orbit(state_vectors: StateVectors, central_body_mass: f32) -> Self {
31        // Compute the magnitude of the position vector and velocity vector
32        let r_mag = state_vectors.position.length();
33        let v_mag = state_vectors.velocity.length();
34
35        // Compute the specific angular momentum vector
36        let h = state_vectors.position.cross(state_vectors.velocity);
37        let h_mag = h.length();
38
39        // Compute the eccentricity vector
40        let e = (state_vectors.velocity.cross(h) / central_body_mass)
41            - (state_vectors.position / r_mag);
42        let e_mag = e.length();
43
44        // Compute the semi-major axis
45        let a = 1.0 / (2.0 / r_mag - v_mag * v_mag / (central_body_mass * central_body_mass));
46
47        // Compute the inclination
48        let i = h.z / h_mag;
49
50        // Compute the longitude of ascending node
51        let n = Vec3::new(-h.y, h.x, 0.0).normalize();
52        let o = n.cross(Vec3::Z);
53        let sign = o.z.signum();
54        let cos_omega = n.dot(Vec3::X) / o.length();
55        let sin_omega = sign * o.y / o.length();
56        let omega = sin_omega.atan2(cos_omega);
57
58        // Compute the argument of periapsis
59        let cos_w = e.dot(n) / (e_mag * n.length());
60        let sin_w = e.dot(o) / (e_mag * o.length());
61        let w = sin_w.atan2(cos_w);
62
63        // Create a new Orbit object with the computed Keplerian elements
64        Orbit {
65            eccentricity: e_mag,
66            semi_major_axis: a,
67            inclination: i.asin(),
68            longitude_of_ascending_node: omega,
69            argument_of_periapsis: w,
70        }
71    }
72
73    /// https://en.wikipedia.org/wiki/Standard_gravitational_parameter
74    pub fn standard_gravitational_parameter(mass: f32) -> f32 {
75        G * mass
76    }
77
78    /// https://en.wikipedia.org/wiki/Orbital_period
79    pub fn period(&self, mass: f32) -> f32 {
80        let sm_cubed = self.semi_major_axis * self.semi_major_axis * self.semi_major_axis;
81
82        2.0 * PI * (sm_cubed / Self::standard_gravitational_parameter(mass)).sqrt()
83    }
84
85    /// https://en.wikipedia.org/wiki/Mean_motion
86    pub fn mean_motion(&self, mass: f32) -> f32 {
87        2.0 * PI / self.period(mass)
88    }
89
90    /// https://en.wikipedia.org/wiki/Mean_anomaly
91    pub fn mean_anomaly(&self, mass: f32, time: f32) -> f32 {
92        self.mean_motion(mass) * time
93    }
94
95    /// https://en.wikipedia.org/wiki/Eccentric_anomaly
96    ///
97    /// Eccentric Anomaly (EA) is given by the equation:
98    /// M = EA - e*sin(EA)
99    /// where
100    /// M is the mean anomaly
101    /// e is the eccentricity
102    ///
103    /// To estimate the Eccentric Anomaly we use Newton's method
104    /// https://en.wikipedia.org/wiki/Newton%27s_method
105    pub fn estimate_eccentric_anomaly(&self, mass: f32, time: f32, tolerance: f32) -> f32 {
106        let mean_anomaly = self.mean_anomaly(mass, time);
107        let mut eccentric_anomaly = mean_anomaly;
108
109        let mut error = 1.0;
110
111        while error > tolerance {
112            // f(E) = E - e*sin(E) - M
113            let fe =
114                eccentric_anomaly - (self.eccentricity * eccentric_anomaly.sin()) - mean_anomaly;
115
116            // f'(E) = 1 - e*cos(E)
117            let fe_prime = 1.0 - (self.eccentricity * eccentric_anomaly.cos());
118
119            let next_eccentric_anomaly = eccentric_anomaly - (fe / fe_prime);
120
121            error = (next_eccentric_anomaly - eccentric_anomaly).abs();
122            eccentric_anomaly = next_eccentric_anomaly;
123        }
124
125        eccentric_anomaly
126    }
127
128    pub fn state_vectors_at_epoch(&self, mass: f32, time: f32, tolerance: f32) -> StateVectors {
129        let Ω = self.longitude_of_ascending_node;
130        let ω = self.argument_of_periapsis;
131        let i = self.inclination;
132        let a = self.semi_major_axis;
133        let e = self.eccentricity;
134        let μ = Self::standard_gravitational_parameter(mass);
135
136        let E = self.estimate_eccentric_anomaly(mass, time, tolerance);
137
138        // Specific angular momentum
139        let h = f32::sqrt(μ * a * (1.0 - e.powi(2)));
140
141        // Perifocal x and y
142        let x = a * (f32::cos(E) - e);
143        let y = a * f32::sqrt(1.0 - e.powi(2)) * f32::sin(E);
144
145        // Perifocal velocity
146        let vx = -f32::sin(E) * h / (a * (1.0 - e * f32::cos(E)));
147        let vy = f32::sqrt(1.0 - e.powi(2)) * f32::cos(E) * h / (a * (1.0 - e * f32::cos(E)));
148
149        // Vectors
150        let position_perifocal = Vec3::new(x, 0.0, y);
151        let velocity_perifocal = Vec3::new(vx, 0.0, vy);
152
153        // Compute rotation matrices to transform perifocal frame to ecliptic frame
154        let i = Mat3::from_axis_angle(Vec3::X, i);
155        let o = Mat3::from_axis_angle(Vec3::Z, Ω);
156        let w = Mat3::from_axis_angle(Vec3::Z, ω);
157
158        // Compute rotation quaternion
159        let q = Quat::from_mat3(&(o * i * w));
160
161        let position = q.mul_vec3(position_perifocal);
162        let velocity = q.mul_vec3(velocity_perifocal);
163
164        StateVectors { position, velocity }
165    }
166}
keplerian_orbits/lib.rs

keplerian_orbits/
lib.rs
