Collision probability computation for conjunction assessment
Upgrade to ocamlformat 0.29.0; fix csvt/sexpt streaming; reformat
- Update .ocamlformat to 0.29.0 across all 591 files
- csvt: reuse single Buffer.t for field reads (no alloc per field)
- sexpt: Obj members decoded from stream into Dict, typed Variant GADT
- Reformat all source files for 0.29.0
+39
-51
+34
-46
lib/collision.ml
+34
-46
lib/collision.ml
···
18
18
hbr : float;
19
19
}
20
20
21
-
(* {1 Vector/Matrix operations} *)
22
-
23
-
let vec3_sub a b = { x = a.x -. b.x; y = a.y -. b.y; z = a.z -. b.z }
24
-
let vec3_dot a b = (a.x *. b.x) +. (a.y *. b.y) +. (a.z *. b.z)
25
-
let vec3_norm v = sqrt (vec3_dot v v)
26
-
27
-
let vec3_cross a b =
28
-
{
29
-
x = (a.y *. b.z) -. (a.z *. b.y);
30
-
y = (a.z *. b.x) -. (a.x *. b.z);
31
-
z = (a.x *. b.y) -. (a.y *. b.x);
32
-
}
33
-
34
-
let vec3_scale s v = { x = s *. v.x; y = s *. v.y; z = s *. v.z }
35
-
36
-
let vec3_normalize v =
37
-
let n = vec3_norm v in
38
-
if n < 1e-15 then v else vec3_scale (1.0 /. n) v
39
-
40
21
(* {1 RTN to ECI covariance rotation}
41
22
42
23
The CDM gives covariance in RTN (Radial-Transverse-Normal) frame.
···
48
29
T = N × R (transverse / in-track) *)
49
30
50
31
let rtn_basis (state : Cdm.state_vector) =
51
-
let r_hat = vec3_normalize state.pos in
52
-
let h = vec3_cross state.pos state.vel in
53
-
let n_hat = vec3_normalize h in
54
-
let t_hat = vec3_cross n_hat r_hat in
32
+
let r_hat = Vec3.normalize state.pos in
33
+
let h = Vec3.cross state.pos state.vel in
34
+
let n_hat = Vec3.normalize h in
35
+
let t_hat = Vec3.cross n_hat r_hat in
55
36
(r_hat, t_hat, n_hat)
56
37
57
38
(* Rotate 3x3 symmetric RTN covariance to ECI.
···
126
107
z_hat = relative velocity direction (normal to plane).
127
108
x_hat is chosen in the miss direction projected onto the plane. *)
128
109
let conjunction_plane_basis ~z_hat ~miss =
129
-
let miss_proj = vec3_sub miss (vec3_scale (vec3_dot miss z_hat) z_hat) in
130
-
let miss_proj_norm = vec3_norm miss_proj in
110
+
let miss_proj = Vec3.sub miss (Vec3.scale (Vec3.dot miss z_hat) z_hat) in
111
+
let miss_proj_norm = Vec3.length miss_proj in
131
112
let x_hat =
132
-
if miss_proj_norm > 1e-12 then vec3_scale (1.0 /. miss_proj_norm) miss_proj
113
+
if miss_proj_norm > 1e-12 then Vec3.scale (1.0 /. miss_proj_norm) miss_proj
133
114
else
134
115
let trial =
135
116
if Float.abs z_hat.x < 0.9 then { x = 1.0; y = 0.0; z = 0.0 }
136
117
else { x = 0.0; y = 1.0; z = 0.0 }
137
118
in
138
-
let perp = vec3_sub trial (vec3_scale (vec3_dot trial z_hat) z_hat) in
139
-
vec3_normalize perp
119
+
let perp = Vec3.sub trial (Vec3.scale (Vec3.dot trial z_hat) z_hat) in
120
+
Vec3.normalize perp
140
121
in
141
-
let y_hat = vec3_cross z_hat x_hat in
122
+
let y_hat = Vec3.cross z_hat x_hat in
142
123
(x_hat, y_hat)
143
124
144
125
(* Project a 3x3 ECI covariance onto two basis vectors, returning a 2x2. *)
···
151
132
z = (c.s13 *. ej.x) +. (c.s23 *. ej.y) +. (c.s33 *. ej.z);
152
133
}
153
134
in
154
-
vec3_dot ei cx
135
+
Vec3.dot ei cx
155
136
in
156
137
(proj x_hat x_hat, proj x_hat y_hat, proj y_hat y_hat)
157
138
···
179
160
(sigma_x, sigma_y, miss_x, miss_y)
180
161
181
162
let encounter_of_cdm ~hbr (cdm : Cdm.t) : encounter =
182
-
let v_rel = vec3_sub cdm.obj2.state.vel cdm.obj1.state.vel in
183
-
let z_hat = vec3_normalize v_rel in
184
-
let miss = vec3_sub cdm.obj2.state.pos cdm.obj1.state.pos in
163
+
let v_rel = Vec3.sub cdm.obj2.state.vel cdm.obj1.state.vel in
164
+
let z_hat = Vec3.normalize v_rel in
165
+
let miss = Vec3.sub cdm.obj2.state.pos cdm.obj1.state.pos in
185
166
let x_hat, y_hat = conjunction_plane_basis ~z_hat ~miss in
186
167
let c1_eci = rtn_to_eci cdm.obj1.cov cdm.obj1.state in
187
168
let c2_eci = rtn_to_eci cdm.obj2.cov cdm.obj2.state in
188
169
let c = sym3_add c1_eci c2_eci in
189
170
let cp11, cp12, cp22 = project_covariance c x_hat y_hat in
190
-
let mx = vec3_dot miss x_hat in
191
-
let my = vec3_dot miss y_hat in
171
+
let mx = Vec3.dot miss x_hat in
172
+
let my = Vec3.dot miss y_hat in
192
173
let sigma_x, sigma_y, miss_x, miss_y =
193
174
diagonalize_and_project cp11 cp12 cp22 mx my
194
175
in
···
258
239
259
240
where f(r,θ) = (r·cos(θ)-mx)²/σx² + (r·sin(θ)-my)²/σy² *)
260
241
242
+
(* Pre-compute quadrature nodes once at module init *)
243
+
let gl_nodes_64, gl_weights_64 = gauss_legendre_nodes 64
244
+
261
245
let pc_foster ?(n = 64) (enc : encounter) =
262
246
let sx = enc.sigma_x in
263
247
let sy = enc.sigma_y in
···
270
254
let sx2 = sx *. sx in
271
255
let sy2 = sy *. sy in
272
256
(* Quadrature over θ ∈ [0, 2π] and r ∈ [0, hbr] *)
273
-
let nodes_t, weights_t = gauss_legendre_nodes n in
274
-
let nodes_r, weights_r = gauss_legendre_nodes n in
257
+
let nodes_t, weights_t =
258
+
if n = 64 then (gl_nodes_64, gl_weights_64) else gauss_legendre_nodes n
259
+
in
260
+
let nodes_r, weights_r =
261
+
if n = 64 then (gl_nodes_64, gl_weights_64) else gauss_legendre_nodes n
262
+
in
275
263
let sum = ref 0.0 in
276
264
for i = 0 to n - 1 do
277
265
(* Map [-1,1] → [0, 2π] *)
···
427
415
if i = 0 || i = n - 1 then
428
416
(* Minimum at edge — cannot refine *)
429
417
let d = sqrt !min_d2 in
430
-
let dv = vec3_norm (vec3_sub vel1.(i) vel2.(i)) in
418
+
let dv = Vec3.length (Vec3.sub vel1.(i) vel2.(i)) in
431
419
Some { time = times.(i); miss_distance = d; relative_velocity = dv }
432
420
else
433
421
(* Step 2: quadratic fit on d^2(t) using 3 points around minimum.
···
447
435
let denom = h0 *. h2 *. (h2 -. h0) in
448
436
if Float.abs denom < 1e-30 then
449
437
let d = sqrt !min_d2 in
450
-
let dv = vec3_norm (vec3_sub vel1.(i) vel2.(i)) in
438
+
let dv = Vec3.length (Vec3.sub vel1.(i) vel2.(i)) in
451
439
Some { time = t1; miss_distance = d; relative_velocity = dv }
452
440
else
453
441
(* a = (h2*(y0-y1) - h0*(y2-y1)) / (h0*h2*(h2-h0)) *)
···
459
447
if a <= 0.0 then
460
448
(* Concave — no minimum, use step-based *)
461
449
let d = sqrt !min_d2 in
462
-
let dv = vec3_norm (vec3_sub vel1.(i) vel2.(i)) in
450
+
let dv = Vec3.length (Vec3.sub vel1.(i) vel2.(i)) in
463
451
Some { time = t1; miss_distance = d; relative_velocity = dv }
464
452
else
465
453
let dt = -.b /. (2.0 *. a) in
···
475
463
(* Relative velocity: interpolate linearly *)
476
464
let frac = (t_star -. t0) /. (t2 -. t0) in
477
465
let frac = Float.max 0.0 (Float.min 1.0 frac) in
478
-
let dv_lo = vec3_norm (vec3_sub vel1.(i - 1) vel2.(i - 1)) in
479
-
let dv_hi = vec3_norm (vec3_sub vel1.(i + 1) vel2.(i + 1)) in
466
+
let dv_lo = Vec3.length (Vec3.sub vel1.(i - 1) vel2.(i - 1)) in
467
+
let dv_hi = Vec3.length (Vec3.sub vel1.(i + 1) vel2.(i + 1)) in
480
468
let dv = dv_lo +. (frac *. (dv_hi -. dv_lo)) in
481
469
Some { time = t_star; miss_distance = d_star; relative_velocity = dv }
482
470
···
498
486
let assess ?(hbr = 0.015) (cdm : Cdm.t) =
499
487
let enc = encounter_of_cdm ~hbr cdm in
500
488
let miss =
501
-
let d = vec3_sub cdm.obj2.state.pos cdm.obj1.state.pos in
502
-
vec3_norm d *. 1000.0
489
+
let d = Vec3.sub cdm.obj2.state.pos cdm.obj1.state.pos in
490
+
Vec3.length d *. 1000.0
503
491
in
504
492
let vrel =
505
-
let dv = vec3_sub cdm.obj2.state.vel cdm.obj1.state.vel in
506
-
vec3_norm dv
493
+
let dv = Vec3.sub cdm.obj2.state.vel cdm.obj1.state.vel in
494
+
Vec3.length dv
507
495
in
508
496
{
509
497
pc = pc_foster enc;
+4
-4
lib/collision.mli
+4
-4
lib/collision.mli
···
3
3
{2 Quick start}
4
4
5
5
{[
6
-
(* Assess a CDM in one line *)
7
-
let a = Collision.assess cdm in
8
-
Printf.printf "Pc = %e, miss = %.0f m, risk = %s\n" a.pc a.miss_distance
9
-
(Collision.risk_level a)
6
+
(* Assess a CDM in one line *)
7
+
let a = Collision.assess cdm in
8
+
Printf.printf "Pc = %e, miss = %.0f m, risk = %s\n" a.pc a.miss_distance
9
+
(Collision.risk_level a)
10
10
]}
11
11
12
12
{b References}: