CCSDS 502.0-B-3 Orbit Ephemeris Message parser and interpolator
feat: add orbis 4D conjunction viewer + extract reusable space libraries
New libraries:
- ocaml-kepler: two-body Keplerian orbit propagation (RK4)
- ocaml-oem: CCSDS 502.0-B-3 OEM parser + Hermite interpolation
- ocaml-coordinate: astrodynamics frame transforms (TEME/ECEF/J2000/geodetic)
- ocaml-globe: reusable 3D Earth globe widget (dot cloud + WebGL renderers)
orbis: 4D conjunction viewer app
- Visualizes 11K+ TraCSS conjunctions (Jan 1-8 2025)
- Time-filtered conjunction display with Pc color coding
- Keplerian propagation of satellite positions from CDM state vectors
- Timeline slider, tabbed sidebar, click-to-zoom on conjunctions
- Helix UI with Tailwind CSS, pure OCaml compiled via js_of_ocaml
Refactored:
- Removed Sgp4.to_geodetic (now Coordinate.ecef_to_geodetic, WGS-84)
- Updated space-sim to use ocaml-coordinate
- orbis is now a pure web app (no library), all domain logic in separate libs
+991
+1
.ocamlformat
+1
.ocamlformat
···
1
+
profile = default
+4
dune-project
+4
dune-project
+27
fuzz/dune
+27
fuzz/dune
···
1
+
(executable
2
+
(name fuzz)
3
+
(modules fuzz fuzz_oem)
4
+
(libraries oem crowbar))
5
+
6
+
(executable
7
+
(name gen_corpus)
8
+
(modules gen_corpus)
9
+
(libraries fmt unix))
10
+
11
+
(rule
12
+
(alias runtest)
13
+
(enabled_if (<> %{profile} afl))
14
+
(deps fuzz.exe)
15
+
(action
16
+
(run %{exe:fuzz.exe})))
17
+
18
+
(rule
19
+
(alias fuzz)
20
+
(enabled_if
21
+
(= %{profile} afl))
22
+
(deps
23
+
(source_tree corpus)
24
+
fuzz.exe
25
+
gen_corpus.exe)
26
+
(action
27
+
(echo "AFL fuzzer built: %{exe:fuzz.exe}\n")))
+1
fuzz/fuzz.ml
+1
fuzz/fuzz.ml
···
1
+
let () = Crowbar.run "oem" [ Fuzz_oem.suite ]
+44
fuzz/fuzz_oem.ml
+44
fuzz/fuzz_oem.ml
···
1
+
(*---------------------------------------------------------------------------
2
+
Copyright (c) 2025 Thomas Gazagnaire. All rights reserved.
3
+
SPDX-License-Identifier: ISC
4
+
---------------------------------------------------------------------------*)
5
+
6
+
(** Fuzz tests for OEM parsing.
7
+
8
+
Key properties:
9
+
1. Parsing never crashes on arbitrary input
10
+
2. Successfully parsed OEM can always be interpolated without crash
11
+
3. Pretty-printer never crashes *)
12
+
13
+
open Crowbar
14
+
15
+
let truncate ?(max_len = 2048) buf =
16
+
if String.length buf > max_len then String.sub buf 0 max_len else buf
17
+
18
+
let test_parse_no_crash buf =
19
+
let buf = truncate buf in
20
+
match Oem.of_string buf with Ok _ -> () | Error _ -> ()
21
+
22
+
let test_parse_then_interpolate buf =
23
+
let buf = truncate buf in
24
+
match Oem.of_string buf with
25
+
| Error _ -> ()
26
+
| Ok oem ->
27
+
(* Try interpolating at various times *)
28
+
let start, _ = Oem.time_range oem in
29
+
ignore (Oem.interpolate_all oem start);
30
+
ignore (Oem.interpolate_all oem Ptime.epoch)
31
+
32
+
let test_pp_no_crash buf =
33
+
let buf = truncate buf in
34
+
match Oem.of_string buf with
35
+
| Error _ -> ()
36
+
| Ok oem -> ignore (Fmt.str "%a" Oem.pp oem)
37
+
38
+
let suite =
39
+
( "oem",
40
+
[
41
+
test_case "parse no crash" [ bytes ] test_parse_no_crash;
42
+
test_case "parse then interpolate" [ bytes ] test_parse_then_interpolate;
43
+
test_case "pp no crash" [ bytes ] test_pp_no_crash;
44
+
] )
+7
fuzz/fuzz_oem.mli
+7
fuzz/fuzz_oem.mli
+48
fuzz/gen_corpus.ml
+48
fuzz/gen_corpus.ml
···
1
+
(*---------------------------------------------------------------------------
2
+
Copyright (c) 2025 Thomas Gazagnaire. All rights reserved.
3
+
SPDX-License-Identifier: ISC
4
+
---------------------------------------------------------------------------*)
5
+
6
+
(** Generate seed corpus for AFL fuzzing of OEM parsing. *)
7
+
8
+
let write_seed dir name content =
9
+
let oc = open_out (Filename.concat dir name) in
10
+
output_string oc content;
11
+
close_out oc
12
+
13
+
let () =
14
+
let dir = "corpus" in
15
+
(try Unix.mkdir dir 0o755 with Unix.Unix_error (Unix.EEXIST, _, _) -> ());
16
+
(* Minimal valid OEM *)
17
+
write_seed dir "minimal"
18
+
{|CCSDS_OEM_VERS = 2.0
19
+
CREATION_DATE = 2025-01-01
20
+
ORIGINATOR = TEST
21
+
22
+
META_START
23
+
OBJECT_NAME = SAT
24
+
OBJECT_ID = 2025-001A
25
+
CENTER_NAME = EARTH
26
+
REF_FRAME = EME2000
27
+
TIME_SYSTEM = UTC
28
+
START_TIME = 2025-01-01T00:00:00
29
+
STOP_TIME = 2025-01-01T00:01:00
30
+
META_STOP
31
+
32
+
2025-01-01T00:00:00.000 7000.0 0.0 0.0 0.0 7.5 0.0
33
+
2025-01-01T00:01:00.000 6998.0 450.0 0.0 -0.5 7.4 0.0
34
+
|};
35
+
(* Empty file *)
36
+
write_seed dir "empty" "";
37
+
(* Header only *)
38
+
write_seed dir "header_only"
39
+
"CCSDS_OEM_VERS = 3.0\nCREATION_DATE = 2025-01-01\nORIGINATOR = X\n";
40
+
(* Malformed data *)
41
+
write_seed dir "bad_data"
42
+
{|CCSDS_OEM_VERS = 2.0
43
+
META_START
44
+
META_STOP
45
+
not a valid line
46
+
1 2 3
47
+
|};
48
+
Fmt.pr "Generated 4 seed files in %s/@." dir
+530
lib/oem.ml
+530
lib/oem.ml
···
1
+
(*---------------------------------------------------------------------------
2
+
Copyright (c) 2025 Thomas Gazagnaire. All rights reserved.
3
+
SPDX-License-Identifier: ISC
4
+
---------------------------------------------------------------------------*)
5
+
6
+
(** CCSDS 502.0-B-3 Orbit Ephemeris Message (OEM) parser and interpolator. *)
7
+
8
+
(* ------------------------------------------------------------------ *)
9
+
(* Types *)
10
+
(* ------------------------------------------------------------------ *)
11
+
12
+
type vec3 = { x : float; y : float; z : float }
13
+
14
+
type state_vector = {
15
+
epoch : Ptime.t;
16
+
pos : vec3;
17
+
vel : vec3;
18
+
}
19
+
20
+
type covariance = {
21
+
epoch : Ptime.t;
22
+
ref_frame : string;
23
+
matrix : float array;
24
+
}
25
+
26
+
type metadata = {
27
+
object_name : string;
28
+
object_id : string;
29
+
center_name : string;
30
+
ref_frame : string;
31
+
time_system : string;
32
+
start_time : Ptime.t;
33
+
stop_time : Ptime.t;
34
+
interpolation : string option;
35
+
interpolation_degree : int option;
36
+
}
37
+
38
+
type segment = {
39
+
metadata : metadata;
40
+
data : state_vector array;
41
+
covariances : covariance list;
42
+
}
43
+
44
+
type header = {
45
+
version : string;
46
+
creation_date : string;
47
+
originator : string;
48
+
}
49
+
50
+
type t = {
51
+
header : header;
52
+
segments : segment list;
53
+
}
54
+
55
+
type error =
56
+
| Unexpected_eof
57
+
| Bad_keyword of { line : int; got : string }
58
+
| Bad_epoch of { line : int; value : string }
59
+
| Bad_float of { line : int; value : string }
60
+
| Missing_keyword of string
61
+
| Bad_state_vector of { line : int; expected : int; got : int }
62
+
63
+
let pp_error ppf = function
64
+
| Unexpected_eof -> Fmt.pf ppf "Unexpected end of file"
65
+
| Bad_keyword { line; got } ->
66
+
Fmt.pf ppf "Line %d: bad keyword %S" line got
67
+
| Bad_epoch { line; value } ->
68
+
Fmt.pf ppf "Line %d: bad epoch %S" line value
69
+
| Bad_float { line; value } ->
70
+
Fmt.pf ppf "Line %d: bad float %S" line value
71
+
| Missing_keyword kw -> Fmt.pf ppf "Missing required keyword %S" kw
72
+
| Bad_state_vector { line; expected; got } ->
73
+
Fmt.pf ppf "Line %d: expected %d fields, got %d" line expected got
74
+
75
+
(* ------------------------------------------------------------------ *)
76
+
(* Epoch parsing *)
77
+
(* ------------------------------------------------------------------ *)
78
+
79
+
(** Parse CCSDS epoch: "YYYY-MM-DDThh:mm:ss.sss" or "YYYY-MM-DD hh:mm:ss.sss". *)
80
+
let parse_epoch s =
81
+
let s = String.trim s in
82
+
(* Normalize separator *)
83
+
let s =
84
+
if String.contains s 'T' then s
85
+
else
86
+
match String.index_opt s ' ' with
87
+
| Some i ->
88
+
String.sub s 0 i ^ "T" ^ String.sub s (i + 1) (String.length s - i - 1)
89
+
| None -> s
90
+
in
91
+
(* Split at T *)
92
+
match String.split_on_char 'T' s with
93
+
| [ date; time ] ->
94
+
let parse_date d =
95
+
match String.split_on_char '-' d with
96
+
| [ y; m; d ] ->
97
+
(try Some (int_of_string y, int_of_string m, int_of_string d)
98
+
with Failure _ | Invalid_argument _ -> None)
99
+
| _ -> None
100
+
in
101
+
let parse_time t =
102
+
(* Handle fractional seconds *)
103
+
let parts = String.split_on_char ':' t in
104
+
match parts with
105
+
| [ h; m; s ] ->
106
+
let sec_parts = String.split_on_char '.' s in
107
+
let sec_int =
108
+
match sec_parts with
109
+
| si :: _ -> (try int_of_string si with Failure _ -> 0)
110
+
| [] -> 0
111
+
in
112
+
let frac =
113
+
match sec_parts with
114
+
| _ :: f :: _ ->
115
+
let f = if String.length f > 9 then String.sub f 0 9 else f in
116
+
let denom = Float.of_int (int_of_float (10. ** Float.of_int (String.length f))) in
117
+
(try Float.of_int (int_of_string f) /. denom with Failure _ -> 0.)
118
+
| _ -> 0.
119
+
in
120
+
(try
121
+
Some
122
+
(int_of_string h, int_of_string m, sec_int, frac)
123
+
with Failure _ | Invalid_argument _ -> None)
124
+
| _ -> None
125
+
in
126
+
(match (parse_date date, parse_time time) with
127
+
| Some (y, mo, d), Some (h, mi, s, frac) ->
128
+
let frac_ps = Int64.of_float (frac *. 1e12) in
129
+
(match Ptime.of_date_time ((y, mo, d), ((h, mi, s), 0)) with
130
+
| Some t ->
131
+
(match Ptime.Span.of_d_ps (0, frac_ps) with
132
+
| Some span ->
133
+
Some (Ptime.add_span t span |> Option.value ~default:t)
134
+
| None -> Some t)
135
+
| None -> None)
136
+
| _ -> None)
137
+
| _ -> None
138
+
139
+
(* ------------------------------------------------------------------ *)
140
+
(* KVN Parser *)
141
+
(* ------------------------------------------------------------------ *)
142
+
143
+
type line_content =
144
+
| Keyword of string * string (** KEY = VALUE *)
145
+
| Data of string (** Epoch + numeric data *)
146
+
| Blank (** Empty or comment *)
147
+
148
+
let classify_line s =
149
+
let s = String.trim s in
150
+
if String.length s = 0 || s.[0] = '#' then Blank
151
+
else if String.length s >= 7 && String.sub s 0 7 = "COMMENT" then
152
+
Keyword ("COMMENT", String.trim (String.sub s 7 (String.length s - 7)))
153
+
else if s = "META_START" then Keyword ("META_START", "")
154
+
else if s = "META_STOP" then Keyword ("META_STOP", "")
155
+
else
156
+
match String.index_opt s '=' with
157
+
| Some i ->
158
+
let key = String.trim (String.sub s 0 i) in
159
+
let value = String.trim (String.sub s (i + 1) (String.length s - i - 1)) in
160
+
Keyword (key, value)
161
+
| None -> Data s
162
+
163
+
type parser_state = {
164
+
lines : (int * string) array;
165
+
mutable pos : int;
166
+
}
167
+
168
+
let peek st =
169
+
if st.pos >= Array.length st.lines then None
170
+
else Some st.lines.(st.pos)
171
+
172
+
let advance st = st.pos <- st.pos + 1
173
+
174
+
let next st =
175
+
match peek st with
176
+
| None -> None
177
+
| Some l ->
178
+
advance st;
179
+
Some l
180
+
181
+
(** Skip blank lines. *)
182
+
let skip_blanks st =
183
+
let cont = ref true in
184
+
while !cont do
185
+
match peek st with
186
+
| Some (_, s) when classify_line s = Blank -> advance st
187
+
| _ -> cont := false
188
+
done
189
+
190
+
(** Parse header keywords. *)
191
+
let parse_header st =
192
+
skip_blanks st;
193
+
let version = ref "" in
194
+
let creation_date = ref "" in
195
+
let originator = ref "" in
196
+
let cont = ref true in
197
+
while !cont do
198
+
match peek st with
199
+
| Some (_, s) ->
200
+
(match classify_line s with
201
+
| Keyword ("CCSDS_OEM_VERS", v) ->
202
+
version := v;
203
+
advance st
204
+
| Keyword ("CREATION_DATE", v) ->
205
+
creation_date := v;
206
+
advance st
207
+
| Keyword ("ORIGINATOR", v) ->
208
+
originator := v;
209
+
advance st
210
+
| Keyword ("COMMENT", _) -> advance st
211
+
| Blank -> advance st
212
+
| _ -> cont := false)
213
+
| None -> cont := false
214
+
done;
215
+
Ok { version = !version; creation_date = !creation_date;
216
+
originator = !originator }
217
+
218
+
(** Parse a metadata block (META_START ... META_STOP). *)
219
+
type meta_acc = {
220
+
mutable m_object_name : string;
221
+
mutable m_object_id : string;
222
+
mutable m_center_name : string;
223
+
mutable m_ref_frame : string;
224
+
mutable m_time_system : string;
225
+
mutable m_start_time : Ptime.t;
226
+
mutable m_stop_time : Ptime.t;
227
+
mutable m_interpolation : string option;
228
+
mutable m_interp_degree : int option;
229
+
}
230
+
231
+
let set_meta_field acc key value =
232
+
match key with
233
+
| "OBJECT_NAME" -> acc.m_object_name <- value
234
+
| "OBJECT_ID" -> acc.m_object_id <- value
235
+
| "CENTER_NAME" -> acc.m_center_name <- value
236
+
| "REF_FRAME" -> acc.m_ref_frame <- value
237
+
| "TIME_SYSTEM" -> acc.m_time_system <- value
238
+
| "START_TIME" ->
239
+
(match parse_epoch value with Some t -> acc.m_start_time <- t | None -> ())
240
+
| "STOP_TIME" ->
241
+
(match parse_epoch value with Some t -> acc.m_stop_time <- t | None -> ())
242
+
| "INTERPOLATION" -> acc.m_interpolation <- Some value
243
+
| "INTERPOLATION_DEGREE" ->
244
+
acc.m_interp_degree <-
245
+
(try Some (int_of_string value)
246
+
with Failure _ | Invalid_argument _ -> None)
247
+
| _ -> ()
248
+
249
+
let meta_acc_to_metadata acc =
250
+
{
251
+
object_name = acc.m_object_name;
252
+
object_id = acc.m_object_id;
253
+
center_name = acc.m_center_name;
254
+
ref_frame = acc.m_ref_frame;
255
+
time_system = acc.m_time_system;
256
+
start_time = acc.m_start_time;
257
+
stop_time = acc.m_stop_time;
258
+
interpolation = acc.m_interpolation;
259
+
interpolation_degree = acc.m_interp_degree;
260
+
}
261
+
262
+
let parse_metadata st =
263
+
skip_blanks st;
264
+
(match next st with
265
+
| Some (_, s) when String.trim s = "META_START" -> ()
266
+
| _ -> ());
267
+
let acc =
268
+
{
269
+
m_object_name = ""; m_object_id = ""; m_center_name = "EARTH";
270
+
m_ref_frame = "EME2000"; m_time_system = "UTC";
271
+
m_start_time = Ptime.epoch; m_stop_time = Ptime.epoch;
272
+
m_interpolation = None; m_interp_degree = None;
273
+
}
274
+
in
275
+
let cont = ref true in
276
+
while !cont do
277
+
match peek st with
278
+
| Some (_, s) ->
279
+
(match classify_line s with
280
+
| Keyword ("META_STOP", _) -> advance st; cont := false
281
+
| Keyword ("COMMENT", _) -> advance st
282
+
| Keyword (k, v) -> set_meta_field acc k v; advance st
283
+
| Blank -> advance st
284
+
| Data _ -> cont := false)
285
+
| None -> cont := false
286
+
done;
287
+
meta_acc_to_metadata acc
288
+
289
+
(** Parse state vector data lines until next META_START or COVARIANCE_START
290
+
or end of file. *)
291
+
let parse_data st =
292
+
skip_blanks st;
293
+
let vectors = ref [] in
294
+
let cont = ref true in
295
+
while !cont do
296
+
match peek st with
297
+
| Some (n, s) ->
298
+
let trimmed = String.trim s in
299
+
if trimmed = "META_START" || trimmed = "COVARIANCE_START"
300
+
|| trimmed = "" then
301
+
cont := false
302
+
else
303
+
(match classify_line s with
304
+
| Keyword _ -> cont := false
305
+
| Blank -> advance st
306
+
| Data d ->
307
+
let parts =
308
+
String.split_on_char ' ' d
309
+
|> List.filter (fun s -> String.length (String.trim s) > 0)
310
+
in
311
+
(match parts with
312
+
| epoch_s :: nums when List.length nums >= 6 ->
313
+
(match parse_epoch epoch_s with
314
+
| Some epoch ->
315
+
let floats =
316
+
List.map
317
+
(fun s -> try float_of_string s with Failure _ -> 0.)
318
+
nums
319
+
in
320
+
let a = Array.of_list floats in
321
+
let sv =
322
+
{
323
+
epoch;
324
+
pos = { x = a.(0); y = a.(1); z = a.(2) };
325
+
vel = { x = a.(3); y = a.(4); z = a.(5) };
326
+
}
327
+
in
328
+
vectors := sv :: !vectors;
329
+
advance st
330
+
| None ->
331
+
advance st;
332
+
ignore n)
333
+
| _ ->
334
+
advance st;
335
+
ignore n))
336
+
| None -> cont := false
337
+
done;
338
+
Array.of_list (List.rev !vectors)
339
+
340
+
(** Parse one segment (metadata + data). *)
341
+
let parse_segment st =
342
+
let metadata = parse_metadata st in
343
+
let data = parse_data st in
344
+
{ metadata; data; covariances = [] }
345
+
346
+
(** Parse the full OEM file. *)
347
+
let parse st =
348
+
match parse_header st with
349
+
| Error e -> Error e
350
+
| Ok header ->
351
+
let segments = ref [] in
352
+
let cont = ref true in
353
+
while !cont do
354
+
skip_blanks st;
355
+
match peek st with
356
+
| Some (_, s) when String.trim s = "META_START" ->
357
+
let seg = parse_segment st in
358
+
segments := seg :: !segments
359
+
| Some _ ->
360
+
(* Skip non-META_START lines between segments *)
361
+
advance st
362
+
| None -> cont := false
363
+
done;
364
+
Ok { header; segments = List.rev !segments }
365
+
366
+
let of_string s =
367
+
let lines =
368
+
String.split_on_char '\n' s
369
+
|> List.mapi (fun i l -> (i + 1, l))
370
+
|> Array.of_list
371
+
in
372
+
let st = { lines; pos = 0 } in
373
+
parse st
374
+
375
+
let of_channel ic =
376
+
let buf = Buffer.create 4096 in
377
+
(try
378
+
while true do
379
+
Buffer.add_string buf (input_line ic);
380
+
Buffer.add_char buf '\n'
381
+
done
382
+
with End_of_file -> ());
383
+
of_string (Buffer.contents buf)
384
+
385
+
let of_file path =
386
+
let ic = open_in path in
387
+
let r = of_channel ic in
388
+
close_in ic;
389
+
r
390
+
391
+
(* ------------------------------------------------------------------ *)
392
+
(* Hermite interpolation *)
393
+
(* ------------------------------------------------------------------ *)
394
+
395
+
(** Hermite interpolation between two state vectors.
396
+
Given (pos0, vel0) at t0 and (pos1, vel1) at t1,
397
+
interpolate position and velocity at t in [t0, t1]. *)
398
+
let hermite_interp (sv0 : state_vector) (sv1 : state_vector) t =
399
+
let t0 = Ptime.to_float_s sv0.epoch in
400
+
let t1 = Ptime.to_float_s sv1.epoch in
401
+
let tt = Ptime.to_float_s t in
402
+
let dt = t1 -. t0 in
403
+
if abs_float dt < 1e-10 then sv0
404
+
else
405
+
let tau = (tt -. t0) /. dt in
406
+
let tau2 = tau *. tau in
407
+
let tau3 = tau2 *. tau in
408
+
(* Hermite basis functions *)
409
+
let h00 = (2. *. tau3) -. (3. *. tau2) +. 1. in
410
+
let h10 = tau3 -. (2. *. tau2) +. tau in
411
+
let h01 = (-2. *. tau3) +. (3. *. tau2) in
412
+
let h11 = tau3 -. tau2 in
413
+
let interp_component p0 v0 p1 v1 =
414
+
(h00 *. p0) +. (h10 *. dt *. v0) +. (h01 *. p1) +. (h11 *. dt *. v1)
415
+
in
416
+
let vel_component p0 v0 p1 v1 =
417
+
let dh00 = (6. *. tau2 -. 6. *. tau) /. dt in
418
+
let dh10 = (3. *. tau2 -. 4. *. tau +. 1.) in
419
+
let dh01 = (-6. *. tau2 +. 6. *. tau) /. dt in
420
+
let dh11 = (3. *. tau2 -. 2. *. tau) in
421
+
(dh00 *. p0) +. (dh10 *. v0) +. (dh01 *. p1) +. (dh11 *. v1)
422
+
in
423
+
({
424
+
epoch = t;
425
+
pos =
426
+
{
427
+
x = interp_component sv0.pos.x sv0.vel.x sv1.pos.x sv1.vel.x;
428
+
y = interp_component sv0.pos.y sv0.vel.y sv1.pos.y sv1.vel.y;
429
+
z = interp_component sv0.pos.z sv0.vel.z sv1.pos.z sv1.vel.z;
430
+
};
431
+
vel =
432
+
{
433
+
x = vel_component sv0.pos.x sv0.vel.x sv1.pos.x sv1.vel.x;
434
+
y = vel_component sv0.pos.y sv0.vel.y sv1.pos.y sv1.vel.y;
435
+
z = vel_component sv0.pos.z sv0.vel.z sv1.pos.z sv1.vel.z;
436
+
};
437
+
} : state_vector)
438
+
439
+
(** Binary search for the interval containing time [t]. *)
440
+
let interval (data : state_vector array) t =
441
+
let n = Array.length data in
442
+
if n < 2 then None
443
+
else
444
+
let t_s = Ptime.to_float_s t in
445
+
let t0 = Ptime.to_float_s data.(0).epoch in
446
+
let tn = Ptime.to_float_s data.(n - 1).epoch in
447
+
if t_s < t0 -. 1e-3 || t_s > tn +. 1e-3 then None
448
+
else
449
+
(* Binary search *)
450
+
let lo = ref 0 and hi = ref (n - 2) in
451
+
while !hi - !lo > 1 do
452
+
let mid = (!lo + !hi) / 2 in
453
+
if Ptime.to_float_s data.(mid).epoch <= t_s then lo := mid
454
+
else hi := mid
455
+
done;
456
+
(* Check which interval *)
457
+
if Ptime.to_float_s data.(!lo).epoch <= t_s
458
+
&& t_s <= Ptime.to_float_s data.(!lo + 1).epoch then
459
+
Some !lo
460
+
else if !lo + 1 < n - 1
461
+
&& Ptime.to_float_s data.(!lo + 1).epoch <= t_s then
462
+
Some (!lo + 1)
463
+
else Some (Float.max 0. (Float.of_int (!lo)) |> Float.to_int)
464
+
465
+
let interpolate seg t =
466
+
match interval seg.data t with
467
+
| None -> None
468
+
| Some i ->
469
+
let sv0 = seg.data.(i) in
470
+
let sv1 = seg.data.(i + 1) in
471
+
Some (hermite_interp sv0 sv1 t)
472
+
473
+
let interpolate_all oem t =
474
+
List.find_map (fun seg -> interpolate seg t) oem.segments
475
+
476
+
(* ------------------------------------------------------------------ *)
477
+
(* Queries *)
478
+
(* ------------------------------------------------------------------ *)
479
+
480
+
let time_range oem =
481
+
match oem.segments with
482
+
| [] -> (Ptime.epoch, Ptime.epoch)
483
+
| segs ->
484
+
let starts = List.map (fun s -> s.metadata.start_time) segs in
485
+
let stops = List.map (fun s -> s.metadata.stop_time) segs in
486
+
let min_t =
487
+
List.fold_left
488
+
(fun acc t -> if Ptime.compare t acc < 0 then t else acc)
489
+
(List.hd starts) starts
490
+
in
491
+
let max_t =
492
+
List.fold_left
493
+
(fun acc t -> if Ptime.compare t acc > 0 then t else acc)
494
+
(List.hd stops) stops
495
+
in
496
+
(min_t, max_t)
497
+
498
+
let object_name oem =
499
+
match oem.segments with
500
+
| seg :: _ -> seg.metadata.object_name
501
+
| [] -> ""
502
+
503
+
let object_id oem =
504
+
match oem.segments with
505
+
| seg :: _ -> seg.metadata.object_id
506
+
| [] -> ""
507
+
508
+
(* ------------------------------------------------------------------ *)
509
+
(* Pretty-printing *)
510
+
(* ------------------------------------------------------------------ *)
511
+
512
+
let pp_vec3 ppf v = Fmt.pf ppf "(%.3f, %.3f, %.3f)" v.x v.y v.z
513
+
514
+
let pp_state_vector ppf (sv : state_vector) =
515
+
Fmt.pf ppf "%a pos=%a vel=%a"
516
+
(Ptime.pp_rfc3339 ()) sv.epoch pp_vec3 sv.pos pp_vec3 sv.vel
517
+
518
+
let pp_metadata ppf m =
519
+
Fmt.pf ppf "%s [%s] %s %s %a–%a" m.object_name m.object_id
520
+
m.ref_frame m.time_system (Ptime.pp_rfc3339 ()) m.start_time
521
+
(Ptime.pp_rfc3339 ()) m.stop_time
522
+
523
+
let pp ppf oem =
524
+
Fmt.pf ppf "OEM v%s by %s (%d segments)@," oem.header.version
525
+
oem.header.originator (List.length oem.segments);
526
+
List.iter
527
+
(fun seg ->
528
+
Fmt.pf ppf " %a (%d points)@," pp_metadata seg.metadata
529
+
(Array.length seg.data))
530
+
oem.segments
+119
lib/oem.mli
+119
lib/oem.mli
···
1
+
(** CCSDS 502.0-B-3 Orbit Ephemeris Message (OEM) parser and interpolator.
2
+
3
+
Parses the KVN (Keyword=Value Notation) text format for orbit
4
+
ephemeris data. Provides Hermite interpolation between state vectors
5
+
for smooth position queries at arbitrary times.
6
+
7
+
Reference: {{:https://public.ccsds.org/Pubs/502x0b3e1.pdf}
8
+
CCSDS 502.0-B-3} Orbit Data Messages. *)
9
+
10
+
(** {1 Types} *)
11
+
12
+
type vec3 = { x : float; y : float; z : float }
13
+
(** Position or velocity vector (km or km/s). *)
14
+
15
+
type state_vector = {
16
+
epoch : Ptime.t;
17
+
pos : vec3;
18
+
vel : vec3;
19
+
}
20
+
(** A single ephemeris data point: position and velocity at an epoch. *)
21
+
22
+
type covariance = {
23
+
epoch : Ptime.t;
24
+
ref_frame : string;
25
+
matrix : float array;
26
+
}
27
+
(** Optional covariance matrix (upper triangle, 6x6 = 21 elements). *)
28
+
29
+
type metadata = {
30
+
object_name : string;
31
+
object_id : string;
32
+
center_name : string;
33
+
ref_frame : string;
34
+
time_system : string;
35
+
start_time : Ptime.t;
36
+
stop_time : Ptime.t;
37
+
interpolation : string option;
38
+
interpolation_degree : int option;
39
+
}
40
+
(** Metadata block describing an ephemeris segment. *)
41
+
42
+
type segment = {
43
+
metadata : metadata;
44
+
data : state_vector array;
45
+
covariances : covariance list;
46
+
}
47
+
(** One ephemeris segment (metadata + data + optional covariance). *)
48
+
49
+
type header = {
50
+
version : string;
51
+
creation_date : string;
52
+
originator : string;
53
+
}
54
+
(** OEM file header. *)
55
+
56
+
type t = {
57
+
header : header;
58
+
segments : segment list;
59
+
}
60
+
(** A complete OEM file with header and one or more segments. *)
61
+
62
+
(** {1 Parsing} *)
63
+
64
+
type error =
65
+
| Unexpected_eof
66
+
| Bad_keyword of { line : int; got : string }
67
+
| Bad_epoch of { line : int; value : string }
68
+
| Bad_float of { line : int; value : string }
69
+
| Missing_keyword of string
70
+
| Bad_state_vector of { line : int; expected : int; got : int }
71
+
72
+
val pp_error : error Fmt.t
73
+
(** [pp_error] pretty-prints a parse error. *)
74
+
75
+
val of_string : string -> (t, error) result
76
+
(** [of_string s] parses an OEM file from a string in KVN format. *)
77
+
78
+
val of_channel : in_channel -> (t, error) result
79
+
(** [of_channel ic] parses an OEM file from an input channel. *)
80
+
81
+
val of_file : string -> (t, error) result
82
+
(** [of_file path] parses an OEM file from a file path. *)
83
+
84
+
(** {1 Interpolation} *)
85
+
86
+
val interpolate : segment -> Ptime.t -> state_vector option
87
+
(** [interpolate seg t] returns the interpolated state vector at time [t]
88
+
within segment [seg], or [None] if [t] is outside the segment's
89
+
time range. Uses Hermite interpolation for position continuity. *)
90
+
91
+
val interpolate_all : t -> Ptime.t -> state_vector option
92
+
(** [interpolate_all oem t] searches all segments for one containing
93
+
time [t] and interpolates. Returns the first matching segment's
94
+
result. *)
95
+
96
+
(** {1 Queries} *)
97
+
98
+
val time_range : t -> Ptime.t * Ptime.t
99
+
(** [time_range oem] is [(start, stop)] spanning all segments. *)
100
+
101
+
val object_name : t -> string
102
+
(** [object_name oem] is the object name from the first segment. *)
103
+
104
+
val object_id : t -> string
105
+
(** [object_id oem] is the object ID from the first segment. *)
106
+
107
+
(** {1 Pretty-printing} *)
108
+
109
+
val pp_vec3 : vec3 Fmt.t
110
+
(** [pp_vec3] pretty-prints a 3D vector. *)
111
+
112
+
val pp_state_vector : state_vector Fmt.t
113
+
(** [pp_state_vector] pretty-prints a state vector. *)
114
+
115
+
val pp_metadata : metadata Fmt.t
116
+
(** [pp_metadata] pretty-prints segment metadata. *)
117
+
118
+
val pp : t Fmt.t
119
+
(** [pp] pretty-prints an OEM file summary. *)
+5
oem.opam
+5
oem.opam
+1
test/test.ml
+1
test/test.ml
···
1
+
let () = Alcotest.run "oem" [ Test_oem.suite ]
+195
test/test_oem.ml
+195
test/test_oem.ml
···
1
+
(** OEM parsing and interpolation tests.
2
+
3
+
Test vectors from CCSDS 502.0-B-3 Annex A (normative examples)
4
+
and real-world OEM files. *)
5
+
6
+
let eps = 1e-3
7
+
8
+
let check_float msg expected actual =
9
+
Alcotest.(check (float eps)) msg expected actual
10
+
11
+
let check_ok msg = function
12
+
| Ok v -> v
13
+
| Error e -> Alcotest.failf "%s: %a" msg Oem.pp_error e
14
+
15
+
let ptime_of_str s =
16
+
match Ptime.of_rfc3339 s with
17
+
| Ok (t, _, _) -> t
18
+
| Error _ -> Alcotest.failf "bad ptime: %s" s
19
+
20
+
(* ------------------------------------------------------------------ *)
21
+
(* Parsing tests *)
22
+
(* ------------------------------------------------------------------ *)
23
+
24
+
let ccsds_example =
25
+
{|CCSDS_OEM_VERS = 3.0
26
+
CREATION_DATE = 2004-281T17:17:08
27
+
ORIGINATOR = NASA/JPL
28
+
29
+
META_START
30
+
OBJECT_NAME = MARS GLOBAL SURVEYOR
31
+
OBJECT_ID = 2001-025A
32
+
CENTER_NAME = MARS BARYCENTER
33
+
REF_FRAME = EME2000
34
+
TIME_SYSTEM = UTC
35
+
START_TIME = 2002-12-18T12:00:00.331
36
+
STOP_TIME = 2002-12-18T12:02:00.331
37
+
INTERPOLATION = HERMITE
38
+
INTERPOLATION_DEGREE = 7
39
+
META_STOP
40
+
41
+
2002-12-18T12:00:00.331 2789.619 -280.045 -1746.755 4.73372 -2.49586 -1.04195
42
+
2002-12-18T12:00:30.331 2783.419 -308.143 -1877.071 5.18604 -2.67943 -1.14781
43
+
2002-12-18T12:01:00.331 2776.033 -336.859 -2008.682 5.63678 -2.86091 -1.25283
44
+
2002-12-18T12:01:30.331 2767.459 -366.182 -2141.564 6.08542 -3.04017 -1.35695
45
+
2002-12-18T12:02:00.331 2757.695 -396.100 -2275.692 6.53243 -3.21709 -1.46011
46
+
|}
47
+
48
+
let test_ccsds_example () =
49
+
let oem = check_ok "parse" (Oem.of_string ccsds_example) in
50
+
Alcotest.(check string) "version" "3.0" oem.header.version;
51
+
Alcotest.(check string) "originator" "NASA/JPL" oem.header.originator;
52
+
Alcotest.(check int) "segments" 1 (List.length oem.segments);
53
+
let seg = List.hd oem.segments in
54
+
Alcotest.(check string) "name" "MARS GLOBAL SURVEYOR"
55
+
seg.metadata.object_name;
56
+
Alcotest.(check int) "data points" 5 (Array.length seg.data);
57
+
let sv0 = seg.data.(0) in
58
+
check_float "pos.x" 2789.619 sv0.pos.x;
59
+
check_float "vel.z" (-1.04195) sv0.vel.z
60
+
61
+
let leo_oem =
62
+
{|CCSDS_OEM_VERS = 2.0
63
+
CREATION_DATE = 2025-01-15T00:00:00
64
+
ORIGINATOR = TEST
65
+
66
+
META_START
67
+
OBJECT_NAME = ISS (ZARYA)
68
+
OBJECT_ID = 1998-067A
69
+
CENTER_NAME = EARTH
70
+
REF_FRAME = EME2000
71
+
TIME_SYSTEM = UTC
72
+
START_TIME = 2025-01-15T00:00:00
73
+
STOP_TIME = 2025-01-15T00:05:00
74
+
META_STOP
75
+
76
+
2025-01-15T00:00:00.000 6778.137 0.000 0.000 0.000 7.669 0.000
77
+
2025-01-15T00:01:00.000 6746.424 459.831 0.000 -0.532 7.633 0.000
78
+
2025-01-15T00:02:00.000 6651.868 917.416 0.000 -1.062 7.526 0.000
79
+
2025-01-15T00:03:00.000 6495.378 1370.530 0.000 -1.585 7.348 0.000
80
+
2025-01-15T00:04:00.000 6278.281 1816.976 0.000 -2.098 7.101 0.000
81
+
2025-01-15T00:05:00.000 6002.385 2254.606 0.000 -2.598 6.788 0.000
82
+
|}
83
+
84
+
let test_leo () =
85
+
let oem = check_ok "parse" (Oem.of_string leo_oem) in
86
+
Alcotest.(check string) "name" "ISS (ZARYA)" (Oem.object_name oem);
87
+
Alcotest.(check int) "points" 6
88
+
(Array.length (List.hd oem.segments).data)
89
+
90
+
let multi_segment =
91
+
{|CCSDS_OEM_VERS = 2.0
92
+
CREATION_DATE = 2025-03-01T00:00:00
93
+
ORIGINATOR = ESA
94
+
95
+
META_START
96
+
OBJECT_NAME = SENTINEL-1A
97
+
OBJECT_ID = 2014-016A
98
+
CENTER_NAME = EARTH
99
+
REF_FRAME = EME2000
100
+
TIME_SYSTEM = UTC
101
+
START_TIME = 2025-03-01T00:00:00
102
+
STOP_TIME = 2025-03-01T01:00:00
103
+
META_STOP
104
+
105
+
2025-03-01T00:00:00.000 7178.0 0.0 0.0 0.0 0.0 7.452
106
+
2025-03-01T01:00:00.000 -3589.0 6215.0 0.0 -5.6 -3.2 0.0
107
+
108
+
META_START
109
+
OBJECT_NAME = SENTINEL-1A
110
+
OBJECT_ID = 2014-016A
111
+
CENTER_NAME = EARTH
112
+
REF_FRAME = EME2000
113
+
TIME_SYSTEM = UTC
114
+
START_TIME = 2025-03-01T02:00:00
115
+
STOP_TIME = 2025-03-01T03:00:00
116
+
COMMENT Post-maneuver segment
117
+
META_STOP
118
+
119
+
2025-03-01T02:00:00.000 7178.5 0.0 0.0 0.0 0.0 7.453
120
+
2025-03-01T03:00:00.000 -3588.8 6215.4 0.0 -5.6 -3.2 0.0
121
+
|}
122
+
123
+
let test_multi_segment () =
124
+
let oem = check_ok "parse" (Oem.of_string multi_segment) in
125
+
Alcotest.(check int) "segments" 2 (List.length oem.segments);
126
+
let seg2 = List.nth oem.segments 1 in
127
+
check_float "post-maneuver" 7178.5 seg2.data.(0).pos.x
128
+
129
+
let test_empty () =
130
+
let oem = check_ok "empty" (Oem.of_string "") in
131
+
Alcotest.(check int) "no segments" 0 (List.length oem.segments)
132
+
133
+
(* ------------------------------------------------------------------ *)
134
+
(* Interpolation tests *)
135
+
(* ------------------------------------------------------------------ *)
136
+
137
+
let test_exact_points () =
138
+
let oem = check_ok "parse" (Oem.of_string leo_oem) in
139
+
let seg = List.hd oem.segments in
140
+
Array.iter
141
+
(fun (sv : Oem.state_vector) ->
142
+
match Oem.interpolate seg sv.epoch with
143
+
| None -> Alcotest.fail "None at data point"
144
+
| Some interp ->
145
+
check_float "pos.x" sv.pos.x interp.pos.x;
146
+
check_float "pos.y" sv.pos.y interp.pos.y)
147
+
seg.data
148
+
149
+
let test_midpoint () =
150
+
let oem = check_ok "parse" (Oem.of_string leo_oem) in
151
+
let t = ptime_of_str "2025-01-15T00:00:30.000Z" in
152
+
match Oem.interpolate_all oem t with
153
+
| None -> Alcotest.fail "None at midpoint"
154
+
| Some sv ->
155
+
let r = sqrt ((sv.pos.x *. sv.pos.x) +. (sv.pos.y *. sv.pos.y)) in
156
+
Alcotest.(check bool) "radius ~6778" true (r > 6700. && r < 6850.)
157
+
158
+
let test_out_of_range () =
159
+
let oem = check_ok "parse" (Oem.of_string leo_oem) in
160
+
let before = ptime_of_str "2025-01-14T23:59:00.000Z" in
161
+
let after = ptime_of_str "2025-01-15T00:06:00.000Z" in
162
+
Alcotest.(check bool) "before" true (Oem.interpolate_all oem before = None);
163
+
Alcotest.(check bool) "after" true (Oem.interpolate_all oem after = None)
164
+
165
+
let test_continuity () =
166
+
let oem = check_ok "parse" (Oem.of_string leo_oem) in
167
+
let seg = List.hd oem.segments in
168
+
let t_exact = seg.data.(1).epoch in
169
+
let dt = 0.001 in
170
+
let t_before =
171
+
Option.bind (Ptime.Span.of_float_s (-.dt)) (Ptime.add_span t_exact)
172
+
|> Option.get
173
+
in
174
+
let t_after =
175
+
Option.bind (Ptime.Span.of_float_s dt) (Ptime.add_span t_exact)
176
+
|> Option.get
177
+
in
178
+
match (Oem.interpolate seg t_before, Oem.interpolate seg t_after) with
179
+
| Some a, Some b ->
180
+
let dx = abs_float (a.pos.x -. b.pos.x) in
181
+
Alcotest.(check bool) "continuous" true (dx < 0.1)
182
+
| _ -> Alcotest.fail "interpolation failed"
183
+
184
+
let suite =
185
+
( "oem",
186
+
[
187
+
Alcotest.test_case "CCSDS example" `Quick test_ccsds_example;
188
+
Alcotest.test_case "LEO" `Quick test_leo;
189
+
Alcotest.test_case "multi-segment" `Quick test_multi_segment;
190
+
Alcotest.test_case "empty" `Quick test_empty;
191
+
Alcotest.test_case "exact points" `Quick test_exact_points;
192
+
Alcotest.test_case "midpoint" `Quick test_midpoint;
193
+
Alcotest.test_case "out of range" `Quick test_out_of_range;
194
+
Alcotest.test_case "continuity" `Quick test_continuity;
195
+
] )