solarsystems 0.0.1

N-body solar system engine — gravitational dynamics, orbital mechanics, perturbations, event detection, and full celestial orchestration
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202

use crate::Vec3;
use sciforge::hub::domain::common::constants::G;
use std::f64::consts::PI;

#[derive(Clone, Debug)]
pub struct OrbitalElements {
    pub semi_major_axis: f64,
    pub eccentricity: f64,
    pub inclination: f64,
    pub longitude_ascending_node: f64,
    pub argument_periapsis: f64,
    pub true_anomaly: f64,
}

impl OrbitalElements {
    pub fn from_state_vectors(pos: Vec3, vel: Vec3, mu: f64) -> Self {
        let r = pos.magnitude();
        let v2 = vel.magnitude_squared();

        let h = pos.cross(&vel);
        let h_mag = h.magnitude();

        let n = Vec3::new(-h.y, h.x, 0.0);
        let n_mag = n.magnitude();

        let e_vec = vel.cross(&h) / mu - pos.normalize();
        let ecc = e_vec.magnitude();

        let energy = 0.5 * v2 - mu / r;
        let sma = if (1.0 - ecc).abs() > 1e-10 {
            -mu / (2.0 * energy)
        } else {
            h_mag * h_mag / mu
        };

        let inc = (h.z / h_mag).clamp(-1.0, 1.0).acos();

        let raan = if n_mag > 1e-10 {
            let val = (n.x / n_mag).clamp(-1.0, 1.0).acos();
            if n.y >= 0.0 { val } else { 2.0 * PI - val }
        } else {
            0.0
        };

        let argp = if n_mag > 1e-10 && ecc > 1e-10 {
            let val = (n.dot(&e_vec) / (n_mag * ecc)).clamp(-1.0, 1.0).acos();
            if e_vec.z >= 0.0 { val } else { 2.0 * PI - val }
        } else if ecc > 1e-10 {
            let val = e_vec.y.atan2(e_vec.x);
            if val < 0.0 { val + 2.0 * PI } else { val }
        } else {
            0.0
        };

        let ta = if ecc > 1e-10 {
            let val = (e_vec.dot(&pos) / (ecc * r)).clamp(-1.0, 1.0).acos();
            if pos.dot(&vel) >= 0.0 {
                val
            } else {
                2.0 * PI - val
            }
        } else {
            0.0
        };

        Self {
            semi_major_axis: sma,
            eccentricity: ecc,
            inclination: inc,
            longitude_ascending_node: raan,
            argument_periapsis: argp,
            true_anomaly: ta,
        }
    }

    pub fn to_state_vectors(&self, mu: f64) -> (Vec3, Vec3) {
        let p = self.semi_major_axis * (1.0 - self.eccentricity * self.eccentricity);
        let r_mag = p / (1.0 + self.eccentricity * self.true_anomaly.cos());

        let r_pqw = Vec3::new(
            r_mag * self.true_anomaly.cos(),
            r_mag * self.true_anomaly.sin(),
            0.0,
        );
        let v_pqw = Vec3::new(
            -(mu / p).sqrt() * self.true_anomaly.sin(),
            (mu / p).sqrt() * (self.eccentricity + self.true_anomaly.cos()),
            0.0,
        );

        let pos = pqw_to_ijk(
            r_pqw,
            self.longitude_ascending_node,
            self.argument_periapsis,
            self.inclination,
        );
        let vel = pqw_to_ijk(
            v_pqw,
            self.longitude_ascending_node,
            self.argument_periapsis,
            self.inclination,
        );
        (pos, vel)
    }

    pub fn orbital_period(&self, mu: f64) -> f64 {
        2.0 * PI * (self.semi_major_axis.powi(3) / mu).sqrt()
    }

    pub fn mean_motion(&self, mu: f64) -> f64 {
        (mu / self.semi_major_axis.powi(3)).sqrt()
    }

    pub fn periapsis(&self) -> f64 {
        self.semi_major_axis * (1.0 - self.eccentricity)
    }

    pub fn apoapsis(&self) -> f64 {
        self.semi_major_axis * (1.0 + self.eccentricity)
    }

    pub fn semi_latus_rectum(&self) -> f64 {
        self.semi_major_axis * (1.0 - self.eccentricity * self.eccentricity)
    }

    pub fn radius_at_true_anomaly(&self, nu: f64) -> f64 {
        self.semi_latus_rectum() / (1.0 + self.eccentricity * nu.cos())
    }

    pub fn velocity_at_radius(&self, r: f64, mu: f64) -> f64 {
        (mu * (2.0 / r - 1.0 / self.semi_major_axis)).abs().sqrt()
    }

    pub fn eccentric_anomaly(&self) -> f64 {
        let e = self.eccentricity;
        let nu = self.true_anomaly;
        ((1.0 - e).sqrt() * (nu / 2.0).sin()).atan2((1.0 + e).sqrt() * (nu / 2.0).cos()) * 2.0
    }

    pub fn mean_anomaly(&self) -> f64 {
        let ea = self.eccentric_anomaly();
        ea - self.eccentricity * ea.sin()
    }

    pub fn true_anomaly_from_mean(mean_anomaly: f64, eccentricity: f64) -> f64 {
        let ea = solve_kepler(mean_anomaly, eccentricity);
        let e = eccentricity;
        let ta =
            2.0 * ((1.0 + e).sqrt() * (ea / 2.0).sin()).atan2((1.0 - e).sqrt() * (ea / 2.0).cos());
        ta.rem_euclid(2.0 * PI)
    }

    pub fn advance(&mut self, dt: f64, mu: f64) {
        let n = self.mean_motion(mu);
        let ma = (self.mean_anomaly() + n * dt).rem_euclid(2.0 * PI);
        self.true_anomaly = Self::true_anomaly_from_mean(ma, self.eccentricity);
    }

    pub fn from_body_state(pos: Vec3, vel: Vec3, central_mass: f64) -> Self {
        Self::from_state_vectors(pos, vel, G * central_mass)
    }
}

pub fn solve_kepler(mean_anomaly: f64, eccentricity: f64) -> f64 {
    let mut ea = mean_anomaly;
    for _ in 0..50 {
        let delta = ea - eccentricity * ea.sin() - mean_anomaly;
        let denom = 1.0 - eccentricity * ea.cos();
        if denom.abs() < 1e-15 {
            break;
        }
        ea -= delta / denom;
        if delta.abs() < 1e-14 {
            break;
        }
    }
    ea
}

fn pqw_to_ijk(v: Vec3, raan: f64, argp: f64, inc: f64) -> Vec3 {
    let (so, co) = raan.sin_cos();
    let (sw, cw) = argp.sin_cos();
    let (si, ci) = inc.sin_cos();
    Vec3::new(
        v.x * (co * cw - so * sw * ci) - v.y * (co * sw + so * cw * ci),
        v.x * (so * cw + co * sw * ci) - v.y * (so * sw - co * cw * ci),
        v.x * (sw * si) + v.y * (cw * si),
    )
}

pub fn hohmann_transfer_dv(r1: f64, r2: f64, mu: f64) -> (f64, f64) {
    let v_circ1 = (mu / r1).sqrt();
    let v_circ2 = (mu / r2).sqrt();
    let a_transfer = (r1 + r2) / 2.0;
    let v_periapsis = (mu * (2.0 / r1 - 1.0 / a_transfer)).sqrt();
    let v_apoapsis = (mu * (2.0 / r2 - 1.0 / a_transfer)).sqrt();
    ((v_periapsis - v_circ1).abs(), (v_circ2 - v_apoapsis).abs())
}

pub fn synodic_period(t1: f64, t2: f64) -> f64 {
    (t1 * t2 / (t1 - t2)).abs()
}