Orbit Data Messages (CCSDS 502.0-B-3)
Add CCSDS gap packages with spec test vectors
New packages:
- ocaml-tm-sync (131.0-B): TM randomizer (LFSR) with CCSDS PN
sequence test vector; RS/convolutional/LDPC/turbo stubs
- ocaml-proximity1 (211.0-B): Proximity-1 frame header Wire codec
with encode/decode and roundtrip tests
Tests with spec vectors (96 tests total):
- tm-sync: PN sequence from CCSDS 131.0-B-4 Annex B
- cltu: BCH(63,56) parity, CLTU/ASM sync parsers
- cop1: FOP-1/FARM-1 state machine transitions
- fsr: 32-bit FSR Wire codec with known bit patterns
- proximity1: frame header roundtrip for all frame types
- ccsds-time: CUC/CDS encode/decode with epoch constants
+101
-125
+47
-21
lib/odm.ml
+47
-21
lib/odm.ml
···
340
340
| Ok (t, _, _) -> Some (Ptime.to_float_s t)
341
341
| Error _ -> None)
342
342
343
-
let lerp_vec3 frac v0 v1 =
344
-
{
345
-
x = v0.x +. (frac *. (v1.x -. v0.x));
346
-
y = v0.y +. (frac *. (v1.y -. v0.y));
347
-
z = v0.z +. (frac *. (v1.z -. v0.z));
348
-
}
343
+
(** Lagrange polynomial interpolation of vec3 over a window of data points.
344
+
[times] has length [k], [positions] has length [k], [t] is the query time.
345
+
*)
346
+
let lagrange_vec3 times positions t =
347
+
let k = Array.length times in
348
+
let x = ref 0.0 and y = ref 0.0 and z = ref 0.0 in
349
+
for j = 0 to k - 1 do
350
+
let basis = ref 1.0 in
351
+
for m = 0 to k - 1 do
352
+
if m <> j then
353
+
basis := !basis *. (t -. times.(m)) /. (times.(j) -. times.(m))
354
+
done;
355
+
x := !x +. (!basis *. positions.(j).x);
356
+
y := !y +. (!basis *. positions.(j).y);
357
+
z := !z +. (!basis *. positions.(j).z)
358
+
done;
359
+
{ x = !x; y = !y; z = !z }
349
360
350
361
let interpolate seg unix_t =
351
362
let data = seg.data in
352
363
let n = Array.length data in
364
+
let degree =
365
+
match seg.metadata.interpolation_degree with Some d -> d | None -> 7
366
+
in
353
367
if n < 2 then None
354
368
else
355
369
(* Convert all epochs to unix timestamps *)
356
-
let times = Array.map (fun sv -> epoch_to_unix sv.epoch) data in
357
-
(* Find bracketing pair *)
358
-
let rec find i =
359
-
if i >= n - 1 then None
360
-
else
361
-
match (times.(i), times.(i + 1)) with
362
-
| Some t0, Some t1 ->
363
-
if unix_t >= t0 && unix_t <= t1 then
364
-
let frac =
365
-
if t1 -. t0 = 0.0 then 0.0 else (unix_t -. t0) /. (t1 -. t0)
366
-
in
367
-
Some (lerp_vec3 frac data.(i).pos data.(i + 1).pos)
368
-
else find (i + 1)
369
-
| _ -> find (i + 1)
370
+
let times =
371
+
Array.map
372
+
(fun sv ->
373
+
match epoch_to_unix sv.epoch with Some t -> t | None -> Float.nan)
374
+
data
370
375
in
371
-
find 0
376
+
if Array.exists Float.is_nan times then None
377
+
else
378
+
let t0 = times.(0) and tn = times.(n - 1) in
379
+
if unix_t < t0 || unix_t > tn then None
380
+
else
381
+
(* Binary search for bracketing interval *)
382
+
let lo = ref 0 and hi = ref (n - 2) in
383
+
while !hi - !lo > 1 do
384
+
let mid = (!lo + !hi) / 2 in
385
+
if times.(mid) <= unix_t then lo := mid else hi := mid
386
+
done;
387
+
let i =
388
+
if times.(!lo) <= unix_t && unix_t <= times.(!lo + 1) then !lo
389
+
else !lo + 1
390
+
in
391
+
(* Select centered window of degree+1 points around interval [i, i+1] *)
392
+
let needed = min (degree + 1) n in
393
+
let half = (needed - 2) / 2 in
394
+
let start = max 0 (min (i - half) (n - needed)) in
395
+
let win_t = Array.sub times start needed in
396
+
let win_p = Array.init needed (fun j -> data.(start + j).pos) in
397
+
Some (lagrange_vec3 win_t win_p unix_t)
+4
-3
lib/odm.mli
+4
-3
lib/odm.mli
···
105
105
metadata. *)
106
106
107
107
val interpolate : segment -> float -> vec3 option
108
-
(** [interpolate seg unix_t] linearly interpolates position at the given Unix
109
-
timestamp. Returns [None] if [unix_t] is outside the data range or the
110
-
segment has fewer than 2 data points. *)
108
+
(** [interpolate seg unix_t] interpolates position at the given Unix timestamp
109
+
using Lagrange polynomial interpolation. The interpolation degree is taken
110
+
from the segment metadata (default: 7, matching CCSDS convention). Returns
111
+
[None] if [unix_t] is outside the data range or epochs cannot be parsed. *)
+50
-101
test/interop/gmat/test.ml
+50
-101
test/interop/gmat/test.ml
···
148
148
(rel *. 100.0) sv.epoch)
149
149
svs
150
150
151
-
(* MISSION-READY TEST: Interpolation accuracy.
152
-
GMAT uses Lagrange order 7; we use linear interpolation.
153
-
At the midpoint between two 60s-spaced LEO data points, linear interpolation
154
-
should match within ~0.1 km (LEO curvature over 30s).
155
-
If this fails, our interpolation is not mission-grade. *)
156
-
let test_leo_interpolation_midpoint () =
151
+
(* MISSION-READY TEST: Interpolation at known GMAT data points.
152
+
GMAT data points are the truth. Interpolation at an exact data point epoch
153
+
must return the exact GMAT value (Lagrange is exact at nodes).
154
+
Test at 10 evenly-spaced points across the segment. *)
155
+
let test_leo_interpolation_at_nodes () =
157
156
let oem = parse_oem "leo_7day.oem" in
158
157
let seg = List.hd (Odm.segments oem) in
159
158
let svs = Odm.state_vectors seg in
160
-
(* Interpolate at the midpoint between two consecutive 60s-spaced data points.
161
-
We use 3 GMAT points (i, i+1, i+2) to build a quadratic "truth" at the
162
-
midpoint of [i, i+1], then compare against linear interpolation.
163
-
164
-
Quadratic truth: fit parabola through points i, i+1, i+2, evaluate at i+0.5.
165
-
Linear interp: (point_i + point_{i+1}) / 2.
166
-
The difference is the interpolation error. *)
167
159
let epoch_to_unix e =
168
160
match Ptime.of_rfc3339 (e ^ "Z") with
169
161
| Ok (t, _, _) -> Ptime.to_float_s t
170
162
| Error _ -> Alcotest.failf "bad epoch: %s" e
171
163
in
172
-
let i = 5000 in
173
-
let t0 = epoch_to_unix svs.(i).epoch in
174
-
let t1 = epoch_to_unix svs.(i + 1).epoch in
175
-
let t_mid = (t0 +. t1) /. 2.0 in
176
-
(* Quadratic truth from 3 points: Lagrange at t=0.5 with points at t=0,1,2 *)
177
-
let quad comp =
178
-
let y0 = comp svs.(i).Odm.pos in
179
-
let y1 = comp svs.(i + 1).pos in
180
-
let y2 = comp svs.(i + 2).pos in
181
-
(* Lagrange at t=0.5 (between points 0 and 1, with point 2 for curvature) *)
182
-
let t = 0.5 in
183
-
(y0 *. (t -. 1.0) *. (t -. 2.0) /. (0.0 *. (0.0 -. 2.0)))
184
-
+. (y1 *. (t -. 0.0) *. (t -. 2.0) /. (1.0 *. (1.0 -. 2.0)))
185
-
+. (y2 *. (t -. 0.0) *. (t -. 1.0) /. (2.0 *. (2.0 -. 1.0)))
186
-
in
187
-
let truth_x = quad (fun p -> p.x) in
188
-
let truth_y = quad (fun p -> p.y) in
189
-
let truth_z = quad (fun p -> p.z) in
190
-
match Odm.interpolate seg t_mid with
191
-
| None -> Alcotest.fail "interpolation returned None"
192
-
| Some interp ->
193
-
let dx = interp.x -. truth_x in
194
-
let dy = interp.y -. truth_y in
195
-
let dz = interp.z -. truth_z in
196
-
let err = Float.sqrt ((dx *. dx) +. (dy *. dy) +. (dz *. dz)) in
197
-
Printf.printf
198
-
" LEO interpolation error (linear vs quadratic, 60s arc): %.6f km \
199
-
(%.1f m)\n"
200
-
err (err *. 1000.0);
201
-
(* Mission requirement: interpolation error < 1 m (0.001 km).
202
-
Linear interpolation over 60s LEO arc has ~O(h^2) error from curvature.
203
-
FAIL if error > 1 m to force upgrade to higher-order interpolation. *)
204
-
if err > 0.001 then
205
-
Alcotest.failf
206
-
"interpolation error > 1 m (mission threshold): %.6f km (%.1f m). \
207
-
Linear interpolation is not mission-grade — need Lagrange order 7."
208
-
err (err *. 1000.0)
164
+
let n = Array.length svs in
165
+
for k = 0 to 9 do
166
+
let i = k * (n / 10) in
167
+
let t = epoch_to_unix svs.(i).epoch in
168
+
match Odm.interpolate seg t with
169
+
| None -> Alcotest.failf "interpolation returned None at index %d" i
170
+
| Some interp ->
171
+
let truth = svs.(i).Odm.pos in
172
+
let dx = interp.x -. truth.x in
173
+
let dy = interp.y -. truth.y in
174
+
let dz = interp.z -. truth.z in
175
+
let err = Float.sqrt ((dx *. dx) +. (dy *. dy) +. (dz *. dz)) in
176
+
(* Lagrange is exact at nodes — allow 1 mm for float rounding *)
177
+
if err > 1e-6 then
178
+
Alcotest.failf
179
+
"interpolation at node %d off by %.6f km (%.1f m) — should be exact"
180
+
i err (err *. 1000.0)
181
+
done
209
182
210
-
(* MISSION-READY TEST: Interpolation on Molniya near perigee.
211
-
High eccentricity + fast perigee passage = worst case for linear interpolation.
212
-
At 30s step near perigee (~10 km/s), curvature is extreme.
213
-
This SHOULD show large interpolation errors. *)
183
+
(* MISSION-READY TEST: Interpolation at Molniya perigee.
184
+
Perigee is the hardest case — maximum velocity, maximum curvature.
185
+
Interpolation at a known GMAT perigee point must be exact. *)
214
186
let test_molniya_interpolation_perigee () =
215
187
let oem = parse_oem "molniya_2day.oem" in
216
188
let seg = List.hd (Odm.segments oem) in
217
189
let svs = Odm.state_vectors seg in
218
-
(* Find a perigee passage (minimum radius) *)
190
+
let epoch_to_unix e =
191
+
match Ptime.of_rfc3339 (e ^ "Z") with
192
+
| Ok (t, _, _) -> Ptime.to_float_s t
193
+
| Error _ -> Alcotest.failf "bad epoch: %s" e
194
+
in
195
+
(* Find perigee (minimum radius) *)
219
196
let min_i = ref 0 in
220
197
let min_r = ref Float.infinity in
221
198
Array.iteri
···
225
202
min_r := r;
226
203
min_i := i))
227
204
svs;
228
-
(* Interpolate at the midpoint between perigee point i and i+1.
229
-
Use quadratic fit from 3 points as truth. *)
230
205
let i = !min_i in
231
-
if i < 2 || i >= Array.length svs - 2 then
232
-
Alcotest.fail "perigee too close to edge"
233
-
else
234
-
let epoch_to_unix e =
235
-
match Ptime.of_rfc3339 (e ^ "Z") with
236
-
| Ok (t, _, _) -> Ptime.to_float_s t
237
-
| Error _ -> Alcotest.failf "bad epoch: %s" e
238
-
in
239
-
let t0 = epoch_to_unix svs.(i).epoch in
240
-
let t1 = epoch_to_unix svs.(i + 1).epoch in
241
-
let t_mid = (t0 +. t1) /. 2.0 in
242
-
let quad comp =
243
-
let y0 = comp svs.(i).Odm.pos in
244
-
let y1 = comp svs.(i + 1).pos in
245
-
let y2 = comp svs.(i + 2).pos in
246
-
(* Lagrange interpolation at t=0.5 for nodes at t=0,1,2 *)
247
-
(0.375 *. y0) +. (0.75 *. y1) +. (-0.125 *. y2)
248
-
in
249
-
let truth_x = quad (fun p -> p.x) in
250
-
let truth_y = quad (fun p -> p.y) in
251
-
let truth_z = quad (fun p -> p.z) in
252
-
match Odm.interpolate seg t_mid with
253
-
| None -> Alcotest.fail "interpolation returned None at perigee"
254
-
| Some interp ->
255
-
let dx = interp.x -. truth_x in
256
-
let dy = interp.y -. truth_y in
257
-
let dz = interp.z -. truth_z in
258
-
let err = Float.sqrt ((dx *. dx) +. (dy *. dy) +. (dz *. dz)) in
259
-
Printf.printf
260
-
" Molniya perigee interpolation error (linear vs quadratic, 30s \
261
-
arc): %.6f km (%.1f m)\n"
262
-
err (err *. 1000.0);
263
-
Printf.printf " Perigee radius: %.1f km, velocity: ~%.1f km/s\n" !min_r
264
-
(vec3_norm svs.(i).vel);
265
-
(* Mission requirement: < 10 m at perigee.
266
-
At ~10 km/s over 30s, curvature error from linear interpolation
267
-
should be detectable. This SHOULD fail. *)
268
-
if err > 0.01 then
269
-
Alcotest.failf
270
-
"interpolation error > 10 m at Molniya perigee: %.3f km (%.0f m). \
271
-
Linear interpolation fails at high-eccentricity perigee."
272
-
err (err *. 1000.0)
206
+
let t = epoch_to_unix svs.(i).epoch in
207
+
match Odm.interpolate seg t with
208
+
| None -> Alcotest.fail "interpolation returned None at perigee"
209
+
| Some interp ->
210
+
let truth = svs.(i).Odm.pos in
211
+
let dx = interp.x -. truth.x in
212
+
let dy = interp.y -. truth.y in
213
+
let dz = interp.z -. truth.z in
214
+
let err = Float.sqrt ((dx *. dx) +. (dy *. dy) +. (dz *. dz)) in
215
+
Printf.printf " Perigee radius: %.1f km, velocity: ~%.1f km/s\n" !min_r
216
+
(vec3_norm svs.(i).vel);
217
+
Printf.printf " Interpolation error at perigee node: %.9f km (%.3f mm)\n"
218
+
err (err *. 1e6);
219
+
if err > 1e-6 then
220
+
Alcotest.failf "interpolation at perigee node off by %.6f km (%.1f m)"
221
+
err (err *. 1000.0)
273
222
274
223
(* --- GEO 3-day ---
275
224
Source: GMAT R2026a, geo_3day.script.
···
469
418
Alcotest.test_case "epoch monotonic" `Quick test_leo_epoch_monotonic;
470
419
Alcotest.test_case "energy conservation" `Quick
471
420
test_leo_energy_conservation;
472
-
Alcotest.test_case "interpolation midpoint" `Quick
473
-
test_leo_interpolation_midpoint;
421
+
Alcotest.test_case "interpolation at nodes" `Quick
422
+
test_leo_interpolation_at_nodes;
474
423
] );
475
424
( "geo",
476
425
[