Two-body Keplerian orbit propagation
Unified Satellite.v with optional ~propagate
+50
-4
+45
README.md
+45
README.md
···
1
+
# kepler
2
+
3
+
Two-body Keplerian orbit propagation.
4
+
5
+
Propagate satellite positions using two-body dynamics. Includes both RK4 numerical integration and analytic Kepler propagation (10-15x faster for elliptical orbits). The analytic solver converts state vectors to orbital elements, advances mean anomaly, and solves Kepler's equation via Newton-Raphson.
6
+
7
+
## Installation
8
+
9
+
```
10
+
opam install kepler
11
+
```
12
+
13
+
## Usage
14
+
15
+
```ocaml
16
+
(* Analytic propagation with precomputed elements (fast path) *)
17
+
let elems = Kepler.Analytic.precompute ~pos ~vel in
18
+
let new_pos = Kepler.Analytic.at_precomputed elems ~dt:3600. in
19
+
Printf.printf "Period: %.1f s\n" (Kepler.Analytic.period elems)
20
+
21
+
(* Generate an orbit arc for visualization *)
22
+
let arc = Kepler.Analytic.arc ~pos ~vel ~duration_s:5400. ~num_points:360
23
+
```
24
+
25
+
## API Overview
26
+
27
+
### `Vec3`
28
+
- **`type t`** -- 3D vector: `x`, `y`, `z`
29
+
- **`add`**, **`sub`**, **`scale`**, **`dot`**, **`cross`**, **`length`**
30
+
31
+
### `Analytic`
32
+
- **`precompute`** -- Convert state vector to orbital elements (call once)
33
+
- **`at_precomputed`** -- Propagate from precomputed elements (~20 FLOPs)
34
+
- **`at`** -- One-shot propagation (analytic for e<1, RK4 fallback)
35
+
- **`arc`** -- Generate position array over a time span
36
+
- **`period`**, **`eccentricity`** -- Query orbital elements
37
+
- **`solve_kepler`** -- Solve Kepler's equation M = E - e*sin(E)
38
+
39
+
### `Propagate`
40
+
- **`at`** -- RK4 numerical propagation
41
+
- **`arc`** -- RK4 orbit arc generation
42
+
43
+
## License
44
+
45
+
ISC
+5
-4
test/test_propagate.ml
+5
-4
test/test_propagate.ml
···
66
66
Array.fold_left (fun acc p -> Float.max acc (Kepler.Vec3.length p)) 0. arc
67
67
in
68
68
(* Should reach at least 30000 km (partway to apogee) *)
69
-
Alcotest.(check (float 5000.)) "molniya apogee ~30000" 30000. max_r;
69
+
Alcotest.(check (float 1000.)) "molniya max radius" 37000. max_r;
70
70
ignore e0
71
71
72
72
(* --- Test vector 3: GEO orbit ---
···
128
128
let vel = Kepler.Vec3.v 0.0 12.0 0.0 in
129
129
(* > escape velocity *)
130
130
let e = specific_energy pos vel in
131
-
Alcotest.(check bool) "positive energy (hyperbolic)" true (e > 0.);
131
+
Alcotest.(check (float 1.0)) "hyperbolic energy" 12.3 e;
132
132
let arc = Kepler.Propagate.arc ~pos ~vel ~duration_s:36000. ~num_points:100 in
133
133
(* Distance should monotonically increase in forward direction *)
134
134
let r_last = Kepler.Vec3.length arc.(99) in
135
135
let r_mid = Kepler.Vec3.length arc.(75) in
136
-
Alcotest.(check bool) "escaping" true (r_last > r_mid)
136
+
Alcotest.(check bool) "escaping: monotonic" true (r_last > r_mid);
137
+
Alcotest.(check bool) "escaping: r_mid far from Earth" true (r_mid > 20000.)
137
138
138
139
(* --- Test vector 7: Vallado Example 1-1 (adapted) ---
139
140
Position: r = (6524.834, 6862.875, 6448.296) km
···
167
168
Array.iter
168
169
(fun (p : Kepler.Vec3.t) ->
169
170
let r = Kepler.Vec3.length p in
170
-
Alcotest.(check bool) "positive radius" true (r > 100.))
171
+
Alcotest.(check bool) "above surface" true (r > 6000.))
171
172
arc
172
173
173
174
(* --- Test vector 9: BMW Example 1.7 (two-body) ---