Two-body Keplerian orbit propagation
Replace Printf/Format with Fmt, rename test_cara
Style: Fmt.pr/str/pf across collision, spacedata, stix, ocm.
Rename test_cara → test_collision_properties.
+66
-66
+1
-1
fuzz/dune
+1
-1
fuzz/dune
+7
-7
fuzz/fuzz_kepler.ml
+7
-7
fuzz/fuzz_kepler.ml
···
26
26
Float.is_finite px && Float.is_finite py && Float.is_finite pz
27
27
&& Float.is_finite vx && Float.is_finite vy && Float.is_finite vz
28
28
then begin
29
-
let pos = Kepler.Vec3.v px py pz in
30
-
let vel = Kepler.Vec3.v vx vy vz in
29
+
let pos = Vec3.v px py pz in
30
+
let vel = Vec3.v vx vy vz in
31
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
···
45
45
Int64.float_of_bits (Bytes.get_int64_le b (i * 8))
46
46
in
47
47
let x = get 0 and y = get 1 and z = get 2 and s = get 3 in
48
-
let v = Kepler.Vec3.v x y z in
49
-
let _ = Kepler.Vec3.length v in
50
-
let _ = Kepler.Vec3.scale s v in
51
-
let _ = Kepler.Vec3.add v v in
52
-
let _ = Fmt.str "%a" Kepler.Vec3.pp v in
48
+
let v = Vec3.v x y z in
49
+
let _ = Vec3.length v in
50
+
let _ = Vec3.scale s v in
51
+
let _ = Vec3.add v v in
52
+
let _ = Fmt.str "%a" Vec3.pp v in
53
53
()
54
54
55
55
let suite =
+1
-1
test/dune
+1
-1
test/dune
+15
-15
test/test_analytic.ml
+15
-15
test/test_analytic.ml
···
45
45
46
46
(** ISS-like circular orbit: analytic and RK4 should agree. *)
47
47
let test_iss_cross_check () =
48
-
let pos = Kepler.Vec3.v 6778. 0. 0. in
49
-
let vel = Kepler.Vec3.v 0. (sqrt (398600.4418 /. 6778.)) 0. in
48
+
let pos = Vec3.v 6778. 0. 0. in
49
+
let vel = Vec3.v 0. (sqrt (398600.4418 /. 6778.)) 0. in
50
50
let dt = 2700. in
51
51
(* 45 minutes *)
52
52
let p_rk4 = Kepler.Propagate.at ~pos ~vel ~dt in
···
59
59
60
60
(** Molniya HEO: analytic and RK4 should agree within tolerance. *)
61
61
let test_molniya_cross_check () =
62
-
let pos = Kepler.Vec3.v 6878. 0. 0. in
63
-
let vel = Kepler.Vec3.v 0. 10.015 0. in
62
+
let pos = Vec3.v 6878. 0. 0. in
63
+
let vel = Vec3.v 0. 10.015 0. in
64
64
let dt = 5400. in
65
65
(* 90 minutes *)
66
66
let p_rk4 = Kepler.Propagate.at ~pos ~vel ~dt in
···
77
77
78
78
(** After one full period, satellite returns to start. *)
79
79
let test_full_period_return () =
80
-
let pos = Kepler.Vec3.v 6778. 0. 0. in
81
-
let vel = Kepler.Vec3.v 0. (sqrt (398600.4418 /. 6778.)) 0. in
80
+
let pos = Vec3.v 6778. 0. 0. in
81
+
let vel = Vec3.v 0. (sqrt (398600.4418 /. 6778.)) 0. in
82
82
let period = 2. *. Float.pi *. sqrt ((6778. ** 3.) /. 398600.4418) in
83
83
let p = Kepler.Analytic.at ~pos ~vel ~dt:period in
84
84
check_float "return x" 2. 6778. p.x;
···
86
86
87
87
(** At dt=0, position is unchanged. *)
88
88
let test_zero_dt () =
89
-
let pos = Kepler.Vec3.v 7000. 1000. 3000. in
90
-
let vel = Kepler.Vec3.v (-1.) 6.5 2. in
89
+
let pos = Vec3.v 7000. 1000. 3000. in
90
+
let vel = Vec3.v (-1.) 6.5 2. in
91
91
let p = Kepler.Analytic.at ~pos ~vel ~dt:0. in
92
92
check_float "dt=0 x" 1e-6 pos.x p.x;
93
93
check_float "dt=0 y" 1e-6 pos.y p.y;
···
95
95
96
96
(** Precomputed elements give same result as direct. *)
97
97
let test_precomputed () =
98
-
let pos = Kepler.Vec3.v 6778. 0. 0. in
99
-
let vel = Kepler.Vec3.v 0. (sqrt (398600.4418 /. 6778.)) 0. in
98
+
let pos = Vec3.v 6778. 0. 0. in
99
+
let vel = Vec3.v 0. (sqrt (398600.4418 /. 6778.)) 0. in
100
100
let el = Kepler.Analytic.precompute ~pos ~vel in
101
101
let p1 = Kepler.Analytic.at ~pos ~vel ~dt:1800. in
102
102
let p2 = Kepler.Analytic.at_precomputed el ~dt:1800. in
···
106
106
107
107
(** Arc should produce correct number of points. *)
108
108
let test_arc_length () =
109
-
let pos = Kepler.Vec3.v 7000. 0. 0. in
110
-
let vel = Kepler.Vec3.v 0. 7.5 0. in
109
+
let pos = Vec3.v 7000. 0. 0. in
110
+
let vel = Vec3.v 0. 7.5 0. in
111
111
let arc = Kepler.Analytic.arc ~pos ~vel ~duration_s:3600. ~num_points:42 in
112
112
Alcotest.(check int) "arc length" 42 (Array.length arc)
113
113
114
114
(** Radius should be approximately constant for circular orbit. *)
115
115
let test_circular_radius () =
116
-
let pos = Kepler.Vec3.v 6778. 0. 0. in
117
-
let vel = Kepler.Vec3.v 0. (sqrt (398600.4418 /. 6778.)) 0. in
116
+
let pos = Vec3.v 6778. 0. 0. in
117
+
let vel = Vec3.v 0. (sqrt (398600.4418 /. 6778.)) 0. in
118
118
let arc = Kepler.Analytic.arc ~pos ~vel ~duration_s:5400. ~num_points:100 in
119
119
Array.iter
120
120
(fun p ->
121
-
let r = Kepler.Vec3.length p in
121
+
let r = Vec3.length p in
122
122
Alcotest.(check (float 2.)) (Fmt.str "r=%.0f" r) 6778. r)
123
123
arc
124
124
+33
-33
test/test_propagate.ml
+33
-33
test/test_propagate.ml
···
18
18
(* --- Energy conservation --- *)
19
19
20
20
(** Specific orbital energy: E = v²/2 - μ/r. Must be conserved. *)
21
-
let specific_energy (pos : Kepler.Vec3.t) (vel : Kepler.Vec3.t) =
22
-
let r = Kepler.Vec3.length pos in
23
-
let v = Kepler.Vec3.length vel in
21
+
let specific_energy (pos : Vec3.t) (vel : Vec3.t) =
22
+
let r = Vec3.length pos in
23
+
let v = Vec3.length vel in
24
24
(v *. v /. 2.) -. (Kepler.Propagate.mu /. r)
25
25
26
26
(* --- Test vector 1: ISS-like circular LEO orbit ---
27
27
From Vallado Example 2-4 (adapted):
28
28
Circular orbit at ~408 km altitude.
29
29
r = 6778 km, v = 7.669 km/s (circular velocity). *)
30
-
let iss_pos = Kepler.Vec3.v 6778.0 0.0 0.0
31
-
let iss_vel = Kepler.Vec3.v 0.0 7.669 0.0
30
+
let iss_pos = Vec3.v 6778.0 0.0 0.0
31
+
let iss_vel = Vec3.v 0.0 7.669 0.0
32
32
33
33
let test_iss_energy_conservation () =
34
34
let e0 = specific_energy iss_pos iss_vel in
···
37
37
~num_points:100
38
38
in
39
39
(* Check that the last point has similar orbital radius (circular orbit) *)
40
-
let r_last = Kepler.Vec3.length arc.(99) in
41
-
let r_first = Kepler.Vec3.length arc.(0) in
40
+
let r_last = Vec3.length arc.(99) in
41
+
let r_first = Vec3.length arc.(0) in
42
42
(* For a circular orbit, all points should be at ~same radius *)
43
43
check_float "radius conservation" 50. r_first r_last;
44
44
(* Energy at first/last points approximately matches *)
···
49
49
Perigee: ~500 km (r_p = 6878 km), Apogee: ~39,873 km (r_a = 46244 km)
50
50
At perigee: v_p = sqrt(μ * (2/r_p - 1/a)) where a = (r_p + r_a)/2
51
51
a = 26561 km, v_p = 10.015 km/s *)
52
-
let molniya_pos = Kepler.Vec3.v 6878.0 0.0 0.0
53
-
let molniya_vel = Kepler.Vec3.v 0.0 10.015 0.0
52
+
let molniya_pos = Vec3.v 6878.0 0.0 0.0
53
+
let molniya_vel = Vec3.v 0.0 10.015 0.0
54
54
55
55
let test_molniya_energy_conservation () =
56
56
let e0 = specific_energy molniya_pos molniya_vel in
···
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 (fun acc p -> Float.max acc (Kepler.Vec3.length p)) 0. arc
66
+
Array.fold_left (fun acc p -> Float.max acc (Vec3.length p)) 0. arc
67
67
in
68
68
(* Should reach at least 30000 km (partway to apogee) *)
69
69
Alcotest.(check (float 1000.)) "molniya max radius" 37000. max_r;
···
73
73
From Curtis Example 2.12:
74
74
GEO radius: 42164 km, circular velocity: 3.075 km/s
75
75
Period: 86164 s (sidereal day). *)
76
-
let geo_pos = Kepler.Vec3.v 42164.0 0.0 0.0
77
-
let geo_vel = Kepler.Vec3.v 0.0 3.075 0.0
76
+
let geo_pos = Vec3.v 42164.0 0.0 0.0
77
+
let geo_vel = Vec3.v 0.0 3.075 0.0
78
78
79
79
let test_geo_period () =
80
80
(* Propagate one full period *)
···
84
84
in
85
85
(* After one full period, should return near starting point *)
86
86
let last = arc.(359) in
87
-
let r_last = Kepler.Vec3.length last in
87
+
let r_last = Vec3.length last in
88
88
check_float "GEO radius after 1 period" 200. 42164. r_last
89
89
90
90
(* --- Test vector 4: Polar orbit ---
91
91
Inclined orbit: position in XZ plane, velocity in Y direction.
92
92
Sun-synchronous ~800 km: r = 7178 km, v = 7.452 km/s *)
93
-
let polar_pos = Kepler.Vec3.v 7178.0 0.0 0.0
94
-
let polar_vel = Kepler.Vec3.v 0.0 0.0 7.452
93
+
let polar_pos = Vec3.v 7178.0 0.0 0.0
94
+
let polar_vel = Vec3.v 0.0 0.0 7.452
95
95
96
96
let test_polar_orbit () =
97
97
let arc =
···
103
103
(* Actually with velocity in Z, the orbit plane is XZ, so Y≈0 *)
104
104
let max_y =
105
105
Array.fold_left
106
-
(fun acc (p : Kepler.Vec3.t) -> Float.max acc (abs_float p.y))
106
+
(fun acc (p : Vec3.t) -> Float.max acc (abs_float p.y))
107
107
0. arc
108
108
in
109
109
Alcotest.(check bool) "stays in XZ plane" true (max_y < 1.0)
···
112
112
Same as ISS but negative velocity: clockwise orbit.
113
113
r = 6778 km, v = -7.669 km/s *)
114
114
let test_retrograde () =
115
-
let pos = Kepler.Vec3.v 6778.0 0.0 0.0 in
116
-
let vel = Kepler.Vec3.v 0.0 (-7.669) 0.0 in
115
+
let pos = Vec3.v 6778.0 0.0 0.0 in
116
+
let vel = Vec3.v 0.0 (-7.669) 0.0 in
117
117
let arc = Kepler.Propagate.arc ~pos ~vel ~duration_s:5400. ~num_points:100 in
118
118
(* Should orbit in opposite direction. After quarter period (~1350s),
119
119
should be at ~(0, -6778, 0) instead of (0, 6778, 0) for prograde *)
···
124
124
From BMW: escape velocity at 6678 km = sqrt(2μ/r) = 10.926 km/s
125
125
With v > v_escape, object should fly away indefinitely. *)
126
126
let test_escape () =
127
-
let pos = Kepler.Vec3.v 6678.0 0.0 0.0 in
128
-
let vel = Kepler.Vec3.v 0.0 12.0 0.0 in
127
+
let pos = Vec3.v 6678.0 0.0 0.0 in
128
+
let vel = Vec3.v 0.0 12.0 0.0 in
129
129
(* > escape velocity *)
130
130
let e = specific_energy pos vel in
131
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
-
let r_last = Kepler.Vec3.length arc.(99) in
135
-
let r_mid = Kepler.Vec3.length arc.(75) in
134
+
let r_last = Vec3.length arc.(99) in
135
+
let r_mid = Vec3.length arc.(75) in
136
136
Alcotest.(check bool) "escaping: monotonic" true (r_last > r_mid);
137
137
Alcotest.(check bool) "escaping: r_mid far from Earth" true (r_mid > 20000.)
138
138
···
142
142
This is a standard test case from Vallado for coordinate transforms
143
143
and orbit determination. We verify energy conservation. *)
144
144
let test_vallado_ex1 () =
145
-
let pos = Kepler.Vec3.v 6524.834 6862.875 6448.296 in
146
-
let vel = Kepler.Vec3.v 4.901327 5.533756 (-1.976341) in
145
+
let pos = Vec3.v 6524.834 6862.875 6448.296 in
146
+
let vel = Vec3.v 4.901327 5.533756 (-1.976341) in
147
147
let e0 = specific_energy pos vel in
148
148
Alcotest.(check bool) "bound orbit" true (e0 < 0.);
149
149
let arc = Kepler.Propagate.arc ~pos ~vel ~duration_s:7200. ~num_points:100 in
150
150
(* Orbit should stay within reasonable bounds *)
151
151
let max_r =
152
-
Array.fold_left (fun acc p -> Float.max acc (Kepler.Vec3.length p)) 0. arc
152
+
Array.fold_left (fun acc p -> Float.max acc (Vec3.length p)) 0. arc
153
153
in
154
154
(* High velocity → large apogee; just verify finite positions *)
155
155
Alcotest.(check bool) "max radius finite" true (Float.is_finite max_r)
···
159
159
let test_near_parabolic () =
160
160
let r = 6678. in
161
161
let v_esc = sqrt (2. *. Kepler.Propagate.mu /. r) in
162
-
let pos = Kepler.Vec3.v r 0. 0. in
163
-
let vel = Kepler.Vec3.v 0. (v_esc -. 0.001) 0. in
162
+
let pos = Vec3.v r 0. 0. in
163
+
let vel = Vec3.v 0. (v_esc -. 0.001) 0. in
164
164
let e = specific_energy pos vel in
165
165
Alcotest.(check bool) "near-zero energy" true (abs_float e < 0.1);
166
166
let arc = Kepler.Propagate.arc ~pos ~vel ~duration_s:36000. ~num_points:50 in
167
167
(* Should still produce valid positions *)
168
168
Array.iter
169
-
(fun (p : Kepler.Vec3.t) ->
170
-
let r = Kepler.Vec3.length p in
169
+
(fun (p : Vec3.t) ->
170
+
let r = Vec3.length p in
171
171
Alcotest.(check bool) "above surface" true (r > 6000.))
172
172
arc
173
173
···
175
175
Verify that propagation is time-reversible:
176
176
propagate forward then backward should return to start. *)
177
177
let test_time_reversibility () =
178
-
let pos = Kepler.Vec3.v 8000.0 0.0 6000.0 in
179
-
let vel = Kepler.Vec3.v 0.0 (-5.0) 6.0 in
178
+
let pos = Vec3.v 8000.0 0.0 6000.0 in
179
+
let vel = Vec3.v 0.0 (-5.0) 6.0 in
180
180
(* Forward 3600s *)
181
181
let fwd = Kepler.Propagate.arc ~pos ~vel ~duration_s:7200. ~num_points:100 in
182
182
(* The midpoint (index 50) should be close to the initial position *)
···
189
189
190
190
(* --- Test vector 10: Array length --- *)
191
191
let test_arc_length () =
192
-
let pos = Kepler.Vec3.v 7000.0 0.0 0.0 in
193
-
let vel = Kepler.Vec3.v 0.0 7.5 0.0 in
192
+
let pos = Vec3.v 7000.0 0.0 0.0 in
193
+
let vel = Vec3.v 0.0 7.5 0.0 in
194
194
let arc = Kepler.Propagate.arc ~pos ~vel ~duration_s:3600. ~num_points:42 in
195
195
Alcotest.(check int) "arc length" 42 (Array.length arc)
196
196
+9
-9
test/test_vec3.ml
+9
-9
test/test_vec3.ml
···
3
3
let check_float msg expected actual =
4
4
Alcotest.(check (float eps)) msg expected actual
5
5
6
-
let check_vec3 msg (ex, ey, ez) (v : Kepler.Vec3.t) =
6
+
let check_vec3 msg (ex, ey, ez) (v : Vec3.t) =
7
7
check_float (msg ^ ".x") ex v.x;
8
8
check_float (msg ^ ".y") ey v.y;
9
9
check_float (msg ^ ".z") ez v.z
10
10
11
11
let test_add () =
12
-
let a = Kepler.Vec3.v 1. 2. 3. in
13
-
let b = Kepler.Vec3.v 4. 5. 6. in
14
-
check_vec3 "add" (5., 7., 9.) (Kepler.Vec3.add a b)
12
+
let a = Vec3.v 1. 2. 3. in
13
+
let b = Vec3.v 4. 5. 6. in
14
+
check_vec3 "add" (5., 7., 9.) (Vec3.add a b)
15
15
16
16
let test_scale () =
17
-
let v = Kepler.Vec3.v 1. 2. 3. in
18
-
check_vec3 "scale" (2., 4., 6.) (Kepler.Vec3.scale 2. v)
17
+
let v = Vec3.v 1. 2. 3. in
18
+
check_vec3 "scale" (2., 4., 6.) (Vec3.scale 2. v)
19
19
20
20
let test_length () =
21
-
let v = Kepler.Vec3.v 3. 4. 0. in
22
-
check_float "length" 5.0 (Kepler.Vec3.length v)
21
+
let v = Vec3.v 3. 4. 0. in
22
+
check_float "length" 5.0 (Vec3.length v)
23
23
24
-
let test_zero () = check_vec3 "zero" (0., 0., 0.) Kepler.Vec3.zero
24
+
let test_zero () = check_vec3 "zero" (0., 0., 0.) Vec3.zero
25
25
26
26
let suite =
27
27
( "vec3",