1use std::f32::consts::{PI, TAU};
2
3const G: f32 = 6.67384e-11;
4
5#[inline]
6pub fn calculate_position_at_time(
7 semi_major_axis: f32,
8 eccentricity: f32,
9 argument_of_periapsis: f32,
10 initial_mean_anomaly: f32,
11 parent_mass: f32,
12 time: f32,
13) -> (f32, f32, f32) {
14 let period = calculate_period(semi_major_axis, parent_mass);
15 let mean_motion = calculate_mean_motion(period);
16 let mean_anomaly = calculate_mean_anomaly(mean_motion, initial_mean_anomaly, time);
17 let eccentric_anomaly = calculate_eccentric_anomaly(eccentricity, mean_anomaly);
18 let true_anomaly = calculate_true_anomaly(eccentricity, eccentric_anomaly);
19 let heliocentric_distance = calculate_heliocentric_distance(semi_major_axis, eccentricity, true_anomaly);
20 calculate_position(true_anomaly, heliocentric_distance, argument_of_periapsis, mean_anomaly)
21}
22
23#[inline]
24pub fn calculate_period(semi_major_axis: f32, parent_mass: f32) -> f32 {
25 TAU * (semi_major_axis.powi(3) / (G * parent_mass)).sqrt()
26}
27
28#[inline]
29pub fn calculate_mean_motion(period: f32) -> f32 {
30 TAU / period
31}
32
33#[inline]
34pub fn calculate_mean_anomaly(mean_motion: f32, initial_mean_anomaly: f32, time: f32) -> f32 {
35 (initial_mean_anomaly + mean_motion * time).rem_euclid(TAU)
36}
37
38#[inline]
39pub fn calculate_initial_mean_anomaly(mean_anomaly: f32, period: f32, time: f32) -> f32 {
40 let mean_motion = calculate_mean_motion(period);
41 (mean_anomaly - mean_motion * time).rem_euclid(TAU)
42}
43
44#[inline]
45pub fn calculate_eccentric_anomaly(eccentricity: f32, mean_anomaly: f32) -> f32 {
46 let e = eccentricity;
47 let ma = mean_anomaly;
48 let mut ea = ma;
49 for _i in 0..5 {
51 ea = ea - (ea - e * ea.sin() - ma) / (1.0 - e * ea.cos());
52 }
53 ea
54}
55
56#[inline]
57pub fn calculate_true_anomaly(eccentricity: f32, eccentric_anomaly: f32) -> f32 {
58 let e = eccentricity;
59 let e_a = eccentric_anomaly;
60 2.0 * (((1.0 + e) / (1.0 - e) * ((e_a / 2.0).tan()).powi(2)).sqrt()).atan()
61}
62
63#[inline]
64pub fn calculate_heliocentric_distance(semi_major_axis: f32, eccentricity: f32, true_anomaly: f32) -> f32 {
65 let semilatus_rectum = semi_major_axis * (1.0 - eccentricity.powi(2));
66 semilatus_rectum / (1.0 + eccentricity * true_anomaly.cos())
67}
68
69#[inline]
70pub fn calculate_position(
71 true_anomaly: f32,
72 heliocentric_distance: f32,
73 argument_of_periapsis: f32,
74 mean_anomaly: f32,
75) -> (f32, f32, f32) {
76 let zmod = if (mean_anomaly % TAU) < PI { -1.0 } else { 1.0 };
77
78 let x = heliocentric_distance * true_anomaly.cos();
79 let z = heliocentric_distance * true_anomaly.sin() * zmod;
80
81 let rotated_x = x * argument_of_periapsis.cos() - z * argument_of_periapsis.sin();
82 let rotated_z = x * argument_of_periapsis.sin() + z * argument_of_periapsis.cos();
83
84 (rotated_x, 0.0, rotated_z)
85}