Collision probability computation for conjunction assessment
Wire RS into tm-sync, implement convolutional codec, add interop vectors
- tm-sync Reed_solomon: now uses ocaml-reed-solomon with CCSDS
interleaving (I=1..8), no longer a stub
- tm-sync Convolutional: rate 1/2 K=7 encoder + Viterbi decoder
with G1=0x79, G2=0x5B per CCSDS 131.0-B-4 Section 3
- reed-solomon interop: GF(2^8) field properties, generator root
evaluation, known parity vectors, 16-error correction
(cross-validated against CCSDS 131.0-B-4 Annex F)
- tm-sync interop: full 255-byte PN sequence from bg2bhc/gr-lilacsat,
convolutional impulse response matching generator polynomials
+7
-4
+2
-2
lib/collision.ml
+2
-2
lib/collision.ml
···
444
444
where c = y1, and a,b from the two other points.
445
445
Minimum at t1 - b/(2a). *)
446
446
let h0 = t0 -. t1 and h2 = t2 -. t1 in
447
-
let denom = (h0 *. h2 *. (h2 -. h0)) in
447
+
let denom = h0 *. h2 *. (h2 -. h0) in
448
448
if Float.abs denom < 1e-30 then
449
449
let d = sqrt !min_d2 in
450
450
let dv = vec3_norm (vec3_sub vel1.(i) vel2.(i)) in
···
454
454
let a = ((h2 *. (y0 -. y1)) -. (h0 *. (y2 -. y1))) /. denom in
455
455
(* b = (h0*h0*(y2-y1) - h2*h2*(y0-y1)) / (h0*h2*(h2-h0)) *)
456
456
let b =
457
-
(((h0 *. h0) *. (y2 -. y1)) -. ((h2 *. h2) *. (y0 -. y1))) /. denom
457
+
((h0 *. h0 *. (y2 -. y1)) -. (h2 *. h2 *. (y0 -. y1))) /. denom
458
458
in
459
459
if a <= 0.0 then
460
460
(* Concave — no minimum, use step-based *)
+5
-2
test/interop/gmat/test.ml
+5
-2
test/interop/gmat/test.ml
···
18
18
let v_sub (a : Vec3.t) (b : Vec3.t) : Vec3.t =
19
19
{ x = a.x -. b.x; y = a.y -. b.y; z = a.z -. b.z }
20
20
21
-
let v_norm (v : Vec3.t) = Float.sqrt ((v.x *. v.x) +. (v.y *. v.y) +. (v.z *. v.z))
21
+
let v_norm (v : Vec3.t) =
22
+
Float.sqrt ((v.x *. v.x) +. (v.y *. v.y) +. (v.z *. v.z))
22
23
23
24
(* Source: GMAT R2026a, conjunction_leo.script.
24
25
Sat1: ISS-like (51.6 deg inc, 400 km alt).
···
66
67
min_dist := d;
67
68
min_i := i)
68
69
done;
69
-
let dv = v_norm (v_sub (to_v sv1.(!min_i).Odm.vel) (to_v sv2.(!min_i).Odm.vel)) in
70
+
let dv =
71
+
v_norm (v_sub (to_v sv1.(!min_i).Odm.vel) (to_v sv2.(!min_i).Odm.vel))
72
+
in
70
73
(* Crossing orbits at ~90 deg: relative velocity ~sqrt(2) * 7.5 ~= 10.6 km/s *)
71
74
Alcotest.(check bool) "high relative velocity" true (dv > 5.0);
72
75
Printf.printf " Relative velocity at TCA: %.3f km/s\n" dv