Two-body Keplerian orbit propagation
feat: orbit ghost/trail split, custom layers, Kepler period
ocaml-kepler:
- Analytic.period: orbital period from elements (seconds)
- Analytic.eccentricity: accessor
- Vec3: add sub, dot, cross
ocaml-globe/webgl/Orbit:
- add_ghost: persistent ghost orbits (survive clear_trails)
- clear_trails: clear only trails, keep ghosts
- clear_ghosts: clear only ghosts
- clear_lines: clear everything (backward compat)
- Draw order: ghosts first, then trails, then dots
ocaml-globe/webgl/Scene:
- layers.custom: list of draw functions called after orbits,
before coverage. For app-specific renderers (conjunction markers).
+72
-87
+34
-39
lib/analytic.ml
+34
-39
lib/analytic.ml
···
5
5
6
6
(** Analytic two-body Kepler propagation.
7
7
8
-
Converts state vector to orbital elements, advances mean anomaly,
9
-
solves Kepler's equation, and converts back. 10-15x faster than RK4
10
-
for elliptical orbits.
8
+
Converts state vector to orbital elements, advances mean anomaly, solves
9
+
Kepler's equation, and converts back. 10-15x faster than RK4 for elliptical
10
+
orbits.
11
11
12
12
References:
13
13
- Vallado, "Fundamentals of Astrodynamics", 4th ed., Alg. 4-5
···
20
20
(** {1 Orbital elements} *)
21
21
22
22
type elements = {
23
-
a : float; (** Semi-major axis (km). *)
24
-
e : float; (** Eccentricity. *)
25
-
i : float; (** Inclination (rad). *)
26
-
raan : float; (** Right ascension of ascending node (rad). *)
27
-
argp : float; (** Argument of periapsis (rad). *)
28
-
nu : float; (** True anomaly (rad). *)
29
-
n : float; (** Mean motion (rad/s). *)
23
+
a : float; (** Semi-major axis (km). *)
24
+
e : float; (** Eccentricity. *)
25
+
i : float; (** Inclination (rad). *)
26
+
raan : float; (** Right ascension of ascending node (rad). *)
27
+
argp : float; (** Argument of periapsis (rad). *)
28
+
nu : float; (** True anomaly (rad). *)
29
+
n : float; (** Mean motion (rad/s). *)
30
30
}
31
31
32
32
(** {1 State vector ↔ orbital elements} *)
33
33
34
-
(** Convert state vector to classical orbital elements.
35
-
Vallado Algorithm 9. *)
34
+
(** Convert state vector to classical orbital elements. Vallado Algorithm 9. *)
36
35
let elements_of_state (pos : Vec3.t) (vel : Vec3.t) =
37
36
let open Vec3 in
38
37
let r_mag = length pos in
···
90
89
(* Mean motion *)
91
90
let n =
92
91
if a > 0. then sqrt (mu /. (a *. a *. a))
93
-
else sqrt (mu /. (-.a *. -.a *. -.a)) (* hyperbolic *)
92
+
else sqrt (mu /. (-.a *. -.a *. -.a))
93
+
(* hyperbolic *)
94
94
in
95
95
{ a; e; i; raan; argp; nu; n }
96
96
97
97
(** {1 Kepler's equation solver} *)
98
98
99
-
(** Solve M = E - e·sin(E) for E given M and e.
100
-
Newton-Raphson iteration, converges in 3-5 steps for e < 0.97.
101
-
Vallado Algorithm 2. *)
99
+
(** Solve M = E - e·sin(E) for E given M and e. Newton-Raphson iteration,
100
+
converges in 3-5 steps for e < 0.97. Vallado Algorithm 2. *)
102
101
let solve_kepler ~mean_anomaly ~eccentricity =
103
102
let m = Float.rem mean_anomaly twopi in
104
103
let m = if m < 0. then m +. twopi else m in
···
116
115
117
116
(** {1 Position from elements + time} *)
118
117
119
-
(** Compute position at time dt (seconds) from epoch.
120
-
The orbital elements describe the orbit at epoch (when ν was measured). *)
118
+
(** Compute position at time dt (seconds) from epoch. The orbital elements
119
+
describe the orbit at epoch (when ν was measured). *)
121
120
let position_at_dt elements dt =
122
121
let el = elements in
123
122
(* True anomaly → eccentric anomaly at epoch *)
124
123
let e = el.e in
125
124
let cos_nu = cos el.nu in
126
-
let ea0 =
127
-
atan2 (sqrt (1. -. (e *. e)) *. sin el.nu) (e +. cos_nu)
128
-
in
125
+
let ea0 = atan2 (sqrt (1. -. (e *. e)) *. sin el.nu) (e +. cos_nu) in
129
126
(* Mean anomaly at epoch *)
130
127
let m0 = ea0 -. (e *. sin ea0) in
131
128
(* Mean anomaly at t + dt *)
···
134
131
let ea = solve_kepler ~mean_anomaly:m ~eccentricity:e in
135
132
(* Eccentric anomaly → true anomaly *)
136
133
let sin_ea = sin ea and cos_ea = cos ea in
137
-
let nu =
138
-
atan2 (sqrt (1. -. (e *. e)) *. sin_ea) (cos_ea -. e)
139
-
in
134
+
let nu = atan2 (sqrt (1. -. (e *. e)) *. sin_ea) (cos_ea -. e) in
140
135
(* Radius *)
141
136
let r = el.a *. (1. -. (e *. cos_ea)) in
142
137
(* Position in orbital plane (PQW frame) *)
···
147
142
let cos_raan = cos el.raan and sin_raan = sin el.raan in
148
143
let cos_i = cos el.i and sin_i = sin el.i in
149
144
let x =
150
-
px *. ((cos_raan *. cos_argp) -. (sin_raan *. sin_argp *. cos_i))
151
-
+. py
152
-
*. (-.(cos_raan *. sin_argp) -. (sin_raan *. cos_argp *. cos_i))
145
+
(px *. ((cos_raan *. cos_argp) -. (sin_raan *. sin_argp *. cos_i)))
146
+
+. (py *. (-.(cos_raan *. sin_argp) -. (sin_raan *. cos_argp *. cos_i)))
153
147
in
154
148
let y =
155
-
px *. ((sin_raan *. cos_argp) +. (cos_raan *. sin_argp *. cos_i))
156
-
+. py
157
-
*. (-.(sin_raan *. sin_argp) +. (cos_raan *. cos_argp *. cos_i))
158
-
in
159
-
let z =
160
-
px *. (sin_argp *. sin_i) +. py *. (cos_argp *. sin_i)
149
+
(px *. ((sin_raan *. cos_argp) +. (cos_raan *. sin_argp *. cos_i)))
150
+
+. (py *. (-.(sin_raan *. sin_argp) +. (cos_raan *. cos_argp *. cos_i)))
161
151
in
152
+
let z = (px *. (sin_argp *. sin_i)) +. (py *. (cos_argp *. sin_i)) in
162
153
Vec3.v x y z
163
154
155
+
(** {1 Element queries} *)
156
+
157
+
let period el = if el.a <= 0. then infinity (* hyperbolic *) else twopi /. el.n
158
+
let eccentricity el = el.e
159
+
164
160
(** {1 Public API} *)
165
161
166
-
(** Propagate a state vector by [dt] seconds using analytic Kepler.
167
-
~50 FLOPs vs ~720 for RK4. Exact for two-body, any eccentricity < 1. *)
162
+
(** Propagate a state vector by [dt] seconds using analytic Kepler. ~50 FLOPs vs
163
+
~720 for RK4. Exact for two-body, any eccentricity < 1. *)
168
164
let at ~pos ~vel ~dt =
169
165
let el = elements_of_state pos vel in
170
166
if el.e >= 1.0 then
···
172
168
Propagate.at ~pos ~vel ~dt
173
169
else position_at_dt el dt
174
170
175
-
(** Precompute orbital elements for repeated propagation.
176
-
Call this once, then use {!at_precomputed} for each time step. *)
171
+
(** Precompute orbital elements for repeated propagation. Call this once, then
172
+
use {!at_precomputed} for each time step. *)
177
173
let precompute ~pos ~vel = elements_of_state pos vel
178
174
179
175
(** Propagate from precomputed elements. Very fast (~20 FLOPs). *)
···
187
183
(** Generate an arc of [num_points] positions centered on epoch. *)
188
184
let arc ~pos ~vel ~duration_s ~num_points =
189
185
let el = elements_of_state pos vel in
190
-
if el.e >= 1.0 then
191
-
Propagate.arc ~pos ~vel ~duration_s ~num_points
186
+
if el.e >= 1.0 then Propagate.arc ~pos ~vel ~duration_s ~num_points
192
187
else
193
188
let dt = duration_s /. Float.of_int num_points in
194
189
let half = num_points / 2 in
+22
-19
lib/analytic.mli
+22
-19
lib/analytic.mli
···
1
1
(** Analytic two-body Kepler propagation.
2
2
3
-
10-15x faster than RK4 for elliptical orbits. Converts state vector
4
-
to orbital elements, advances mean anomaly, solves Kepler's equation
5
-
via Newton-Raphson, and converts back to position.
3
+
10-15x faster than RK4 for elliptical orbits. Converts state vector to
4
+
orbital elements, advances mean anomaly, solves Kepler's equation via
5
+
Newton-Raphson, and converts back to position.
6
6
7
7
Falls back to RK4 for hyperbolic/parabolic orbits (e >= 1). *)
8
8
···
12
12
(** Precomputed orbital elements for fast repeated propagation. *)
13
13
14
14
val precompute : pos:Vec3.t -> vel:Vec3.t -> elements
15
-
(** [precompute ~pos ~vel] converts a state vector to orbital elements.
16
-
Call once per satellite, then use {!at_precomputed} each frame. *)
15
+
(** [precompute ~pos ~vel] converts a state vector to orbital elements. Call
16
+
once per satellite, then use {!at_precomputed} each frame. *)
17
+
18
+
val period : elements -> float
19
+
(** [period el] is the orbital period in seconds. Returns [infinity] for
20
+
hyperbolic/parabolic orbits. *)
21
+
22
+
val eccentricity : elements -> float
23
+
(** [eccentricity el] is the orbital eccentricity. *)
17
24
18
25
(** {1 Propagation} *)
19
26
20
27
val at : pos:Vec3.t -> vel:Vec3.t -> dt:float -> Vec3.t
21
-
(** [at ~pos ~vel ~dt] propagates by [dt] seconds. Uses analytic Kepler
22
-
for bound orbits (e < 1), RK4 fallback for hyperbolic. *)
28
+
(** [at ~pos ~vel ~dt] propagates by [dt] seconds. Uses analytic Kepler for
29
+
bound orbits (e < 1), RK4 fallback for hyperbolic. *)
23
30
24
31
val at_precomputed : elements -> dt:float -> Vec3.t
25
-
(** [at_precomputed el ~dt] propagates from precomputed elements.
26
-
Very fast (~20 FLOPs per call). *)
32
+
(** [at_precomputed el ~dt] propagates from precomputed elements. Very fast (~20
33
+
FLOPs per call). *)
27
34
28
35
val arc :
29
-
pos:Vec3.t ->
30
-
vel:Vec3.t ->
31
-
duration_s:float ->
32
-
num_points:int ->
33
-
Vec3.t array
34
-
(** [arc ~pos ~vel ~duration_s ~num_points] generates positions
35
-
centered on epoch. Uses analytic Kepler for bound orbits. *)
36
+
pos:Vec3.t -> vel:Vec3.t -> duration_s:float -> num_points:int -> Vec3.t array
37
+
(** [arc ~pos ~vel ~duration_s ~num_points] generates positions centered on
38
+
epoch. Uses analytic Kepler for bound orbits. *)
36
39
37
40
(** {1 Kepler equation} *)
38
41
39
42
val solve_kepler : mean_anomaly:float -> eccentricity:float -> float
40
-
(** [solve_kepler ~mean_anomaly ~eccentricity] solves Kepler's equation
41
-
M = E - e·sin(E) for eccentric anomaly E. Newton-Raphson, 3-5
42
-
iterations for e < 0.97. *)
43
+
(** [solve_kepler ~mean_anomaly ~eccentricity] solves Kepler's equation M = E -
44
+
e·sin(E) for eccentric anomaly E. Newton-Raphson, 3-5 iterations for e <
45
+
0.97. *)
+16
-29
test/test_analytic.ml
+16
-29
test/test_analytic.ml
···
1
1
(** Analytic Kepler solver tests.
2
2
3
-
Validates against RK4 (cross-check), known orbital properties,
4
-
and Kepler equation test vectors. *)
3
+
Validates against RK4 (cross-check), known orbital properties, and Kepler
4
+
equation test vectors. *)
5
5
6
6
let eps_km = 2.0 (* analytic vs RK4 agreement within 2km *)
7
7
···
18
18
19
19
(** For e=0.5, M=1.0: known solution E ≈ 1.4987 (Vallado). *)
20
20
let test_kepler_moderate () =
21
-
let e =
22
-
Kepler.Analytic.solve_kepler ~mean_anomaly:1.0 ~eccentricity:0.5
23
-
in
21
+
let e = Kepler.Analytic.solve_kepler ~mean_anomaly:1.0 ~eccentricity:0.5 in
24
22
(* Verify: M = E - e·sin(E) *)
25
23
let m_check = e -. (0.5 *. sin e) in
26
24
check_float "M roundtrip" 1e-10 1.0 m_check
27
25
28
26
(** High eccentricity (e=0.9) still converges. *)
29
27
let test_kepler_high_ecc () =
30
-
let e =
31
-
Kepler.Analytic.solve_kepler ~mean_anomaly:0.5 ~eccentricity:0.9
32
-
in
28
+
let e = Kepler.Analytic.solve_kepler ~mean_anomaly:0.5 ~eccentricity:0.9 in
33
29
let m_check = e -. (0.9 *. sin e) in
34
30
check_float "high-e roundtrip" 1e-8 0.5 m_check
35
31
36
32
(** M = 0 → E = 0 for any eccentricity. *)
37
33
let test_kepler_zero () =
38
-
let e =
39
-
Kepler.Analytic.solve_kepler ~mean_anomaly:0. ~eccentricity:0.7
40
-
in
34
+
let e = Kepler.Analytic.solve_kepler ~mean_anomaly:0. ~eccentricity:0.7 in
41
35
check_float "M=0 → E=0" 1e-10 0. e
42
36
43
37
(** M = π → E = π for any eccentricity (apoapsis). *)
44
38
let test_kepler_pi () =
45
39
let e =
46
-
Kepler.Analytic.solve_kepler ~mean_anomaly:Float.pi
47
-
~eccentricity:0.5
40
+
Kepler.Analytic.solve_kepler ~mean_anomaly:Float.pi ~eccentricity:0.5
48
41
in
49
42
check_float "M=π → E=π" 1e-8 Float.pi e
50
43
···
54
47
let test_iss_cross_check () =
55
48
let pos = Kepler.Vec3.v 6778. 0. 0. in
56
49
let vel = Kepler.Vec3.v 0. 7.669 0. in
57
-
let dt = 2700. in (* 45 minutes *)
50
+
let dt = 2700. in
51
+
(* 45 minutes *)
58
52
let p_rk4 = Kepler.Propagate.at ~pos ~vel ~dt in
59
53
let p_ana = Kepler.Analytic.at ~pos ~vel ~dt in
60
54
let dx = abs_float (p_rk4.x -. p_ana.x) in
···
69
63
let test_molniya_cross_check () =
70
64
let pos = Kepler.Vec3.v 6878. 0. 0. in
71
65
let vel = Kepler.Vec3.v 0. 10.015 0. in
72
-
let dt = 5400. in (* 90 minutes *)
66
+
let dt = 5400. in
67
+
(* 90 minutes *)
73
68
let p_rk4 = Kepler.Propagate.at ~pos ~vel ~dt in
74
69
let p_ana = Kepler.Analytic.at ~pos ~vel ~dt in
75
70
let dx = p_rk4.x -. p_ana.x in
···
86
81
let test_full_period_return () =
87
82
let pos = Kepler.Vec3.v 6778. 0. 0. in
88
83
let vel = Kepler.Vec3.v 0. 7.669 0. in
89
-
let period = 2. *. Float.pi *. sqrt (6778. ** 3. /. 398600.4418) in
84
+
let period = 2. *. Float.pi *. sqrt ((6778. ** 3.) /. 398600.4418) in
90
85
let p = Kepler.Analytic.at ~pos ~vel ~dt:period in
91
86
check_float "return x" 10. 6778. p.x;
92
87
check_float "return y" 10. 0. p.y
···
115
110
let test_arc_length () =
116
111
let pos = Kepler.Vec3.v 7000. 0. 0. in
117
112
let vel = Kepler.Vec3.v 0. 7.5 0. in
118
-
let arc =
119
-
Kepler.Analytic.arc ~pos ~vel ~duration_s:3600. ~num_points:42
120
-
in
113
+
let arc = Kepler.Analytic.arc ~pos ~vel ~duration_s:3600. ~num_points:42 in
121
114
Alcotest.(check int) "arc length" 42 (Array.length arc)
122
115
123
116
(** Radius should be approximately constant for circular orbit. *)
124
117
let test_circular_radius () =
125
118
let pos = Kepler.Vec3.v 6778. 0. 0. in
126
119
let vel = Kepler.Vec3.v 0. 7.669 0. in
127
-
let arc =
128
-
Kepler.Analytic.arc ~pos ~vel ~duration_s:5400. ~num_points:100
129
-
in
120
+
let arc = Kepler.Analytic.arc ~pos ~vel ~duration_s:5400. ~num_points:100 in
130
121
Array.iter
131
122
(fun p ->
132
123
let r = Kepler.Vec3.length p in
133
-
Alcotest.(check bool)
134
-
(Fmt.str "r=%.0f" r)
135
-
true (r > 6700. && r < 6850.))
124
+
Alcotest.(check bool) (Fmt.str "r=%.0f" r) true (r > 6700. && r < 6850.))
136
125
arc
137
126
138
127
let suite =
···
144
133
Alcotest.test_case "kepler M=0" `Quick test_kepler_zero;
145
134
Alcotest.test_case "kepler M=pi" `Quick test_kepler_pi;
146
135
Alcotest.test_case "ISS cross-check" `Quick test_iss_cross_check;
147
-
Alcotest.test_case "Molniya cross-check" `Quick
148
-
test_molniya_cross_check;
149
-
Alcotest.test_case "full period return" `Quick
150
-
test_full_period_return;
136
+
Alcotest.test_case "Molniya cross-check" `Quick test_molniya_cross_check;
137
+
Alcotest.test_case "full period return" `Quick test_full_period_return;
151
138
Alcotest.test_case "zero dt" `Quick test_zero_dt;
152
139
Alcotest.test_case "precomputed" `Quick test_precomputed;
153
140
Alcotest.test_case "arc length" `Quick test_arc_length;