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

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

use crate::Vec3;

pub struct TransferResult {
    pub delta_v1: f64,
    pub delta_v2: f64,
    pub total_delta_v: f64,
    pub transfer_time: f64,
    pub semi_major_axis: f64,
}

pub fn hohmann_transfer(r1: f64, r2: f64, mu: f64) -> TransferResult {
    let a_t = (r1 + r2) / 2.0;
    let v_circ1 = (mu / r1).sqrt();
    let v_circ2 = (mu / r2).sqrt();
    let v_peri = (mu * (2.0 / r1 - 1.0 / a_t)).sqrt();
    let v_apo = (mu * (2.0 / r2 - 1.0 / a_t)).sqrt();
    let dv1 = (v_peri - v_circ1).abs();
    let dv2 = (v_circ2 - v_apo).abs();
    let tof = PI * (a_t.powi(3) / mu).sqrt();
    TransferResult {
        delta_v1: dv1,
        delta_v2: dv2,
        total_delta_v: dv1 + dv2,
        transfer_time: tof,
        semi_major_axis: a_t,
    }
}

pub struct BiEllipticResult {
    pub delta_v1: f64,
    pub delta_v2: f64,
    pub delta_v3: f64,
    pub total_delta_v: f64,
    pub transfer_time: f64,
}

pub fn bielliptic_transfer(r1: f64, r2: f64, r_intermediate: f64, mu: f64) -> BiEllipticResult {
    let a1 = (r1 + r_intermediate) / 2.0;
    let a2 = (r2 + r_intermediate) / 2.0;

    let v_circ1 = (mu / r1).sqrt();
    let v_circ2 = (mu / r2).sqrt();

    let v1_depart = (mu * (2.0 / r1 - 1.0 / a1)).sqrt();
    let v1_arrive = (mu * (2.0 / r_intermediate - 1.0 / a1)).sqrt();
    let v2_depart = (mu * (2.0 / r_intermediate - 1.0 / a2)).sqrt();
    let v2_arrive = (mu * (2.0 / r2 - 1.0 / a2)).sqrt();

    let dv1 = (v1_depart - v_circ1).abs();
    let dv2 = (v2_depart - v1_arrive).abs();
    let dv3 = (v_circ2 - v2_arrive).abs();

    let tof = PI * (a1.powi(3) / mu).sqrt() + PI * (a2.powi(3) / mu).sqrt();

    BiEllipticResult {
        delta_v1: dv1,
        delta_v2: dv2,
        delta_v3: dv3,
        total_delta_v: dv1 + dv2 + dv3,
        transfer_time: tof,
    }
}

pub fn hohmann_vs_bielliptic_threshold(r1: f64, mu: f64) -> f64 {
    11.94 * r1 * (mu / (mu + 1.0)).powf(0.0)
}

pub fn synodic_period(t1: f64, t2: f64) -> f64 {
    let diff = (1.0 / t1 - 1.0 / t2).abs();
    if diff < 1e-15 {
        f64::INFINITY
    } else {
        1.0 / diff
    }
}

pub fn launch_window_phase_angle(r1: f64, r2: f64, mu: f64) -> f64 {
    let a_t = (r1 + r2) / 2.0;
    let tof = PI * (a_t.powi(3) / mu).sqrt();
    let omega2 = (mu / r2.powi(3)).sqrt();
    PI - omega2 * tof
}

pub fn launch_window_next(phase_current: f64, phase_required: f64, synodic: f64) -> f64 {
    let diff = (phase_required - phase_current) % (2.0 * PI);
    let diff = if diff < 0.0 { diff + 2.0 * PI } else { diff };
    diff / (2.0 * PI) * synodic
}

pub fn characteristic_energy(v_infinity: f64) -> f64 {
    v_infinity * v_infinity
}

pub fn v_infinity_from_c3(c3: f64) -> f64 {
    c3.abs().sqrt()
}

pub fn hyperbolic_excess_velocity(v_departure: f64, v_circular: f64) -> f64 {
    (v_departure * v_departure - v_circular * v_circular)
        .abs()
        .sqrt()
}

pub fn gravity_assist_delta_v(v_infinity_in: f64, flyby_mass: f64, periapsis_radius: f64) -> f64 {
    let mu = G * flyby_mass;
    let deflection = 2.0 * (mu / (periapsis_radius * v_infinity_in * v_infinity_in + mu)).asin();
    2.0 * v_infinity_in * (deflection / 2.0).sin()
}

pub fn gravity_assist_max_deflection(v_infinity: f64, flyby_mass: f64, min_periapsis: f64) -> f64 {
    let mu = G * flyby_mass;
    2.0 * (mu / (min_periapsis * v_infinity * v_infinity + mu)).asin()
}

pub fn oberth_effect_dv(v_periapsis: f64, delta_v_applied: f64) -> f64 {
    let e_before = 0.5 * v_periapsis * v_periapsis;
    let v_after = v_periapsis + delta_v_applied;
    let e_after = 0.5 * v_after * v_after;
    (2.0 * (e_after - e_before)).abs().sqrt()
}

pub fn plane_change_dv(v: f64, angle: f64) -> f64 {
    2.0 * v * (angle / 2.0).sin()
}

pub fn combined_plane_change_dv(v1: f64, v2: f64, angle: f64) -> f64 {
    (v1 * v1 + v2 * v2 - 2.0 * v1 * v2 * angle.cos()).sqrt()
}

pub fn lambert_tof_estimate(r1: Vec3, r2: Vec3, mu: f64, prograde: bool) -> f64 {
    let r1_mag = r1.magnitude();
    let r2_mag = r2.magnitude();
    let cos_dnu = r1.dot(&r2) / (r1_mag * r2_mag);
    let sin_dnu = if prograde {
        r1.cross(&r2).z.signum() * (1.0 - cos_dnu * cos_dnu).sqrt()
    } else {
        -r1.cross(&r2).z.signum() * (1.0 - cos_dnu * cos_dnu).sqrt()
    };
    let dnu = sin_dnu.atan2(cos_dnu);
    let dnu = if dnu < 0.0 { dnu + 2.0 * PI } else { dnu };

    let k = r1_mag * r2_mag * (1.0 - cos_dnu);
    let l = r1_mag + r2_mag;
    let m = r1_mag * r2_mag * (1.0 + cos_dnu);

    let p_min = k / (l + (2.0 * m).sqrt());
    let a = m * k * p_min / (2.0 * m - l * l * p_min);
    let f = 1.0 - r2_mag / p_min * (1.0 - cos_dnu);
    let g = r1_mag * r2_mag * sin_dnu / (mu * p_min).sqrt();

    let tof_parabolic = g.abs();
    let tof_elliptic = PI * (a.abs().powi(3) / mu).sqrt();
    let correction = dnu.sin() * f.abs() * 1e-15;
    tof_parabolic.max(tof_elliptic) + correction
}

pub fn vis_viva(mu: f64, r: f64, a: f64) -> f64 {
    (mu * (2.0 / r - 1.0 / a)).abs().sqrt()
}

pub fn circular_velocity(mu: f64, r: f64) -> f64 {
    (mu / r).sqrt()
}

pub fn escape_velocity(mu: f64, r: f64) -> f64 {
    (2.0 * mu / r).sqrt()
}

pub fn flight_path_angle(r: f64, v: f64, h: f64) -> f64 {
    let cos_fpa = h / (r * v);
    cos_fpa.clamp(-1.0, 1.0).acos()
}

pub fn delta_v_budget(maneuvers: &[f64]) -> f64 {
    maneuvers.iter().map(|dv| dv.abs()).sum()
}

pub fn tsiolkovsky_mass_ratio(delta_v: f64, isp: f64) -> f64 {
    let g0 = 9.80665;
    (delta_v / (isp * g0)).exp()
}

pub fn burn_time(delta_v: f64, thrust: f64, initial_mass: f64) -> f64 {
    initial_mass * delta_v / thrust
}