Two-body Keplerian orbit propagation
fix(ocaml-requests): update tests and fuzz for cstruct→Bytes migration
Test files still referenced Cstruct.t where the API now uses bytes.
Fixed all H2 frame, HPACK, client, and connection tests.
Fixed fuzz test. 330 tests pass.
+72
-60
+12
dune-project
+12
dune-project
···
2
2
(name kepler)
3
3
(source (tangled gazagnaire.org/ocaml-kepler))
4
4
(formatting (enabled_for ocaml))
5
+
6
+
(generate_opam_files true)
7
+
8
+
(license ISC)
9
+
(authors "Thomas Gazagnaire <thomas@gazagnaire.org>")
10
+
(maintainers "Thomas Gazagnaire <thomas@gazagnaire.org>")
11
+
12
+
(package
13
+
(name kepler)
14
+
(synopsis "Two-body Keplerian orbit propagation")
15
+
(depends
16
+
fmt))
+7
-7
fuzz/fuzz_kepler.ml
+7
-7
fuzz/fuzz_kepler.ml
···
5
5
6
6
(** Fuzz tests for Keplerian propagation.
7
7
8
-
Key properties tested:
9
-
1. Propagation never crashes (no unhandled exceptions)
10
-
2. No NaN/infinity in output positions for valid inputs
11
-
3. Vec3 operations are safe on all floats *)
8
+
Key properties tested: 1. Propagation never crashes (no unhandled
9
+
exceptions) 2. No NaN/infinity in output positions for valid inputs 3. Vec3
10
+
operations are safe on all floats *)
12
11
13
12
open Alcobar
14
13
···
23
22
in
24
23
let px = get 0 and py = get 1 and pz = get 2 in
25
24
let vx = get 3 and vy = get 4 and vz = get 5 in
26
-
if Float.is_finite px && Float.is_finite py && Float.is_finite pz
27
-
&& Float.is_finite vx && Float.is_finite vy && Float.is_finite vz
25
+
if
26
+
Float.is_finite px && Float.is_finite py && Float.is_finite pz
27
+
&& Float.is_finite vx && Float.is_finite vy && Float.is_finite vz
28
28
then begin
29
29
let pos = Kepler.Vec3.v px py pz in
30
30
let vel = Kepler.Vec3.v vx vy vz in
31
-
let dur = abs_float (Float.of_int n) *. 10. +. 1. in
31
+
let dur = (abs_float (Float.of_int n) *. 10.) +. 1. in
32
32
let _arc =
33
33
Kepler.Propagate.arc ~pos ~vel ~duration_s:dur ~num_points:10
34
34
in
+26
-2
kepler.opam
+26
-2
kepler.opam
···
1
1
# This file is generated by dune, edit dune-project instead
2
2
opam-version: "2.0"
3
-
name: "kepler"
4
3
synopsis: "Two-body Keplerian orbit propagation"
5
-
depends: []
4
+
maintainer: ["Thomas Gazagnaire <thomas@gazagnaire.org>"]
5
+
authors: ["Thomas Gazagnaire <thomas@gazagnaire.org>"]
6
+
license: "ISC"
7
+
homepage: "https://tangled.org/gazagnaire.org/ocaml-kepler"
8
+
bug-reports: "https://tangled.org/gazagnaire.org/ocaml-kepler/issues"
9
+
depends: [
10
+
"dune" {>= "3.21"}
11
+
"fmt"
12
+
"odoc" {with-doc}
13
+
]
14
+
build: [
15
+
["dune" "subst"] {dev}
16
+
[
17
+
"dune"
18
+
"build"
19
+
"-p"
20
+
name
21
+
"-j"
22
+
jobs
23
+
"@install"
24
+
"@runtest" {with-test}
25
+
"@doc" {with-doc}
26
+
]
27
+
]
28
+
dev-repo: "git+https://tangled.org/gazagnaire.org/ocaml-kepler"
29
+
x-maintenance-intent: ["(latest)"]
+4
-6
lib/propagate.ml
+4
-6
lib/propagate.ml
···
5
5
6
6
(** Two-body (Keplerian) orbit propagation via RK4 integration.
7
7
8
-
Given position and velocity at epoch, propagates forward/backward
9
-
using numerical integration. Suitable for short arcs (~1 orbit). *)
8
+
Given position and velocity at epoch, propagates forward/backward using
9
+
numerical integration. Suitable for short arcs (~1 orbit). *)
10
10
11
11
let mu = 398600.4418
12
12
···
29
29
let a4 = gravity p4 in
30
30
let new_pos =
31
31
add pos
32
-
(scale (dt /. 6.)
33
-
(add vel (add (scale 2. v2) (add (scale 2. v3) v4))))
32
+
(scale (dt /. 6.) (add vel (add (scale 2. v2) (add (scale 2. v3) v4))))
34
33
in
35
34
let new_vel =
36
35
add vel
37
-
(scale (dt /. 6.)
38
-
(add a1 (add (scale 2. a2) (add (scale 2. a3) a4))))
36
+
(scale (dt /. 6.) (add a1 (add (scale 2. a2) (add (scale 2. a3) a4))))
39
37
in
40
38
(new_pos, new_vel)
41
39
+6
-10
lib/propagate.mli
+6
-10
lib/propagate.mli
···
4
4
(** [mu] is Earth's gravitational parameter (km^3/s^2). *)
5
5
6
6
val at : pos:Vec3.t -> vel:Vec3.t -> dt:float -> Vec3.t
7
-
(** [at ~pos ~vel ~dt] propagates a state vector by [dt] seconds
8
-
and returns the new position. Uses adaptive step count. *)
7
+
(** [at ~pos ~vel ~dt] propagates a state vector by [dt] seconds and returns the
8
+
new position. Uses adaptive step count. *)
9
9
10
10
val arc :
11
-
pos:Vec3.t ->
12
-
vel:Vec3.t ->
13
-
duration_s:float ->
14
-
num_points:int ->
15
-
Vec3.t array
16
-
(** [arc ~pos ~vel ~duration_s ~num_points] propagates a state vector
17
-
over [duration_s] seconds centered on epoch, returning [num_points]
18
-
positions in the same frame as the input. *)
11
+
pos:Vec3.t -> vel:Vec3.t -> duration_s:float -> num_points:int -> Vec3.t array
12
+
(** [arc ~pos ~vel ~duration_s ~num_points] propagates a state vector over
13
+
[duration_s] seconds centered on epoch, returning [num_points] positions in
14
+
the same frame as the input. *)
+16
-33
test/test_propagate.ml
+16
-33
test/test_propagate.ml
···
6
6
- Curtis, "Orbital Mechanics for Engineering Students"
7
7
- STK validation data (circular, elliptical, polar orbits) *)
8
8
9
-
let eps_km = 1.0
10
9
(** Position accuracy ~1 km after short propagation (RK4 vs analytical). *)
10
+
let eps_km = 1.0
11
11
12
-
let eps_energy = 1e-4
13
12
(** Specific energy conservation tolerance. *)
13
+
let eps_energy = 1e-4
14
14
15
15
let check_float msg eps expected actual =
16
16
Alcotest.(check (float eps)) msg expected actual
···
56
56
let e0 = specific_energy molniya_pos molniya_vel in
57
57
(* Propagate half an orbit (~6 hours = 21600s) *)
58
58
let arc =
59
-
Kepler.Propagate.arc ~pos:molniya_pos ~vel:molniya_vel
60
-
~duration_s:21600. ~num_points:200
59
+
Kepler.Propagate.arc ~pos:molniya_pos ~vel:molniya_vel ~duration_s:21600.
60
+
~num_points:200
61
61
in
62
62
(* Energy at various points should be conserved *)
63
63
(* The arc array gives positions; we need to recompute vel for energy.
64
64
Instead, verify that the orbit reaches near-apogee altitude. *)
65
65
let max_r =
66
-
Array.fold_left
67
-
(fun acc p -> Float.max acc (Kepler.Vec3.length p))
68
-
0. arc
66
+
Array.fold_left (fun acc p -> Float.max acc (Kepler.Vec3.length p)) 0. arc
69
67
in
70
68
(* Should reach at least 30000 km (partway to apogee) *)
71
69
Alcotest.(check bool) "reaches high altitude" true (max_r > 25000.);
···
116
114
let test_retrograde () =
117
115
let pos = Kepler.Vec3.v 6778.0 0.0 0.0 in
118
116
let vel = Kepler.Vec3.v 0.0 (-7.669) 0.0 in
119
-
let arc =
120
-
Kepler.Propagate.arc ~pos ~vel ~duration_s:5400. ~num_points:100
121
-
in
117
+
let arc = Kepler.Propagate.arc ~pos ~vel ~duration_s:5400. ~num_points:100 in
122
118
(* Should orbit in opposite direction. After quarter period (~1350s),
123
119
should be at ~(0, -6778, 0) instead of (0, 6778, 0) for prograde *)
124
120
let quarter = arc.(75) in
···
129
125
With v > v_escape, object should fly away indefinitely. *)
130
126
let test_escape () =
131
127
let pos = Kepler.Vec3.v 6678.0 0.0 0.0 in
132
-
let vel = Kepler.Vec3.v 0.0 12.0 0.0 in (* > escape velocity *)
128
+
let vel = Kepler.Vec3.v 0.0 12.0 0.0 in
129
+
(* > escape velocity *)
133
130
let e = specific_energy pos vel in
134
131
Alcotest.(check bool) "positive energy (hyperbolic)" true (e > 0.);
135
-
let arc =
136
-
Kepler.Propagate.arc ~pos ~vel ~duration_s:36000. ~num_points:100
137
-
in
132
+
let arc = Kepler.Propagate.arc ~pos ~vel ~duration_s:36000. ~num_points:100 in
138
133
(* Distance should monotonically increase in forward direction *)
139
134
let r_last = Kepler.Vec3.length arc.(99) in
140
135
let r_mid = Kepler.Vec3.length arc.(75) in
···
150
145
let vel = Kepler.Vec3.v 4.901327 5.533756 (-1.976341) in
151
146
let e0 = specific_energy pos vel in
152
147
Alcotest.(check bool) "bound orbit" true (e0 < 0.);
153
-
let arc =
154
-
Kepler.Propagate.arc ~pos ~vel ~duration_s:7200. ~num_points:100
155
-
in
148
+
let arc = Kepler.Propagate.arc ~pos ~vel ~duration_s:7200. ~num_points:100 in
156
149
(* Orbit should stay within reasonable bounds *)
157
150
let max_r =
158
-
Array.fold_left
159
-
(fun acc p -> Float.max acc (Kepler.Vec3.length p))
160
-
0. arc
151
+
Array.fold_left (fun acc p -> Float.max acc (Kepler.Vec3.length p)) 0. arc
161
152
in
162
153
(* High velocity → large apogee; just verify finite positions *)
163
154
Alcotest.(check bool) "max radius finite" true (Float.is_finite max_r)
···
171
162
let vel = Kepler.Vec3.v 0. (v_esc -. 0.001) 0. in
172
163
let e = specific_energy pos vel in
173
164
Alcotest.(check bool) "near-zero energy" true (abs_float e < 0.1);
174
-
let arc =
175
-
Kepler.Propagate.arc ~pos ~vel ~duration_s:36000. ~num_points:50
176
-
in
165
+
let arc = Kepler.Propagate.arc ~pos ~vel ~duration_s:36000. ~num_points:50 in
177
166
(* Should still produce valid positions *)
178
167
Array.iter
179
168
(fun (p : Kepler.Vec3.t) ->
···
188
177
let pos = Kepler.Vec3.v 8000.0 0.0 6000.0 in
189
178
let vel = Kepler.Vec3.v 0.0 (-5.0) 6.0 in
190
179
(* Forward 3600s *)
191
-
let fwd =
192
-
Kepler.Propagate.arc ~pos ~vel ~duration_s:7200. ~num_points:100
193
-
in
180
+
let fwd = Kepler.Propagate.arc ~pos ~vel ~duration_s:7200. ~num_points:100 in
194
181
(* The midpoint (index 50) should be close to the initial position *)
195
182
let mid = fwd.(50) in
196
183
let dx = abs_float (mid.x -. pos.x) in
···
203
190
let test_arc_length () =
204
191
let pos = Kepler.Vec3.v 7000.0 0.0 0.0 in
205
192
let vel = Kepler.Vec3.v 0.0 7.5 0.0 in
206
-
let arc =
207
-
Kepler.Propagate.arc ~pos ~vel ~duration_s:3600. ~num_points:42
208
-
in
193
+
let arc = Kepler.Propagate.arc ~pos ~vel ~duration_s:3600. ~num_points:42 in
209
194
Alcotest.(check int) "arc length" 42 (Array.length arc)
210
195
211
196
let suite =
···
213
198
[
214
199
Alcotest.test_case "ISS energy conservation" `Quick
215
200
test_iss_energy_conservation;
216
-
Alcotest.test_case "Molniya HEO" `Quick
217
-
test_molniya_energy_conservation;
201
+
Alcotest.test_case "Molniya HEO" `Quick test_molniya_energy_conservation;
218
202
Alcotest.test_case "GEO period" `Quick test_geo_period;
219
203
Alcotest.test_case "polar orbit" `Quick test_polar_orbit;
220
204
Alcotest.test_case "retrograde" `Quick test_retrograde;
221
205
Alcotest.test_case "escape trajectory" `Quick test_escape;
222
206
Alcotest.test_case "Vallado ex1" `Quick test_vallado_ex1;
223
207
Alcotest.test_case "near-parabolic" `Quick test_near_parabolic;
224
-
Alcotest.test_case "time reversibility" `Quick
225
-
test_time_reversibility;
208
+
Alcotest.test_case "time reversibility" `Quick test_time_reversibility;
226
209
Alcotest.test_case "arc length" `Quick test_arc_length;
227
210
] )
+1
-2
test/test_vec3.ml
+1
-2
test/test_vec3.ml