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
|
use crate::solar_system::{Angle, BodyId, Kilometers, OrbitalBody, Percentage, SystemId};
use crate::timeman::{Second};
pub struct StaticOrbit
{
parent: BodyId,
system: SystemId,
eccentricity: Percentage,
inclination: Angle,
long_asc_node: Angle,
long_periapsis: Angle,
mean_long: Angle,
semi_major_axis: Kilometers,
}
pub trait StaticOrbiter
{
fn orbit(&self) -> Option<&StaticOrbit>;
fn sgp(&self) -> f64;
}
impl StaticOrbit
{
pub fn new(
parent: BodyId,
system: SystemId,
eccentricity: Percentage,
inclination: Angle,
long_asc_node: Angle,
long_periapsis: Angle,
mean_long: Angle,
semi_major_axis: Kilometers)
-> Self {
Self {
parent,
system,
eccentricity,
inclination,
long_asc_node,
long_periapsis,
mean_long,
semi_major_axis
}
}
pub fn new_circular(
parent: &OrbitalBody,
radius: Kilometers)
-> Self
{
Self {
parent: parent.id(),
system: parent.system(),
eccentricity: 1e-6,
inclination: 0.0,
long_asc_node: 0.0,
long_periapsis: 0.0,
mean_long: 0.0,
semi_major_axis: radius
}
}
pub fn period<T: StaticOrbiter>(
&self,
this_body: &T)
-> Second
{
((self.semi_major_axis.powf(3.0) / this_body.sgp()).sqrt() * std::f64::consts::TAU) as Second
}
pub fn parent(&self) -> BodyId { self.parent }
pub fn system(&self) -> SystemId { self.system }
pub fn sma(&self) -> Kilometers { self.semi_major_axis }
pub fn calculate_position_at<T: StaticOrbiter>(
&self,
this_body: &T,
time: Second)
-> cgmath::Vector3<f64>
{
let eccentricity = self.eccentricity;
let long_asc_node = self.long_asc_node;
let arg_periaps = self.long_periapsis - long_asc_node;
let sma_cubed = self.semi_major_axis.powf(3.0);
let mean_motion = (this_body.sgp() / sma_cubed).sqrt();
let mean_anomaly_epoch = self.mean_long - self.long_periapsis;
let mean_anomaly = mean_anomaly_epoch + (mean_motion * time as f64);
//Find the eccentric anomaly via newton's method
let mut eccentric_anomaly = mean_anomaly;
for _ in 0..100 {
let (e_sin, e_cos) = eccentric_anomaly.sin_cos();
let new_anomaly = eccentric_anomaly -
(
(eccentric_anomaly - eccentricity * e_sin - mean_anomaly) /
(1.0 - eccentricity * e_cos)
);
let diff = eccentric_anomaly - new_anomaly;
eccentric_anomaly = new_anomaly;
if diff.abs() < 1e-6 {
break;
}
}
let true_anomaly = 2.0 * (
((1.0 + eccentricity) / (1.0 - eccentricity)).sqrt()
* (eccentric_anomaly / 2.0).tan()
).atan();
let vw = true_anomaly + arg_periaps;
let (vw_sin, vw_cos) = vw.sin_cos();
let (omega_sin, omega_cos) = long_asc_node.sin_cos();
let (i_sin, i_cos) = self.inclination.sin_cos();
let r = self.semi_major_axis * (1.0 - eccentricity * eccentric_anomaly.cos());
let x = r * (omega_cos * vw_cos - omega_sin * vw_sin * i_cos);
let y = r * i_sin * vw_sin;
let z = r * (omega_sin * vw_cos + omega_cos * vw_sin * i_cos);
cgmath::Vector3::<f64>::new(x, y, z)
}
}
|
