Collision probability computation for conjunction assessment
Simplify public APIs for odm, collision, and cam
odm: Add Odm.position_at — queries position across all segments in one
call. No need to manually iterate segments or convert epochs.
collision: Add Collision.assess — one-call CDM risk assessment returning
Pc, miss distance, relative velocity, and risk level (Critical/High/
Watch/Low enum). encounter/pc_foster/pc_chan/pc_max remain available for
custom workflows.
cam: Add Cam.avoid — one-call CDM avoidance. Takes a CDM and time-to-burn,
returns the minimum maneuver to reduce Pc below threshold (default 1e-5).
Extracts conjunction geometry from CDM automatically.
+99
-48
+40
lib/collision.ml
+40
lib/collision.ml
···
486
486
let enc = encounter_of_cdm ~hbr cdm in
487
487
pc_foster enc
488
488
489
+
(* {1 High-level assessment} *)
490
+
491
+
type assessment = {
492
+
pc : float;
493
+
pc_max : float;
494
+
miss_distance : float;
495
+
relative_velocity : float;
496
+
}
497
+
498
+
let assess ?(hbr = 0.015) (cdm : Cdm.t) =
499
+
let enc = encounter_of_cdm ~hbr cdm in
500
+
let miss =
501
+
let d = vec3_sub cdm.obj2.state.pos cdm.obj1.state.pos in
502
+
vec3_norm d *. 1000.0
503
+
in
504
+
let vrel =
505
+
let dv = vec3_sub cdm.obj2.state.vel cdm.obj1.state.vel in
506
+
vec3_norm dv
507
+
in
508
+
{ pc = pc_foster enc; pc_max = pc_max enc; miss_distance = miss;
509
+
relative_velocity = vrel }
510
+
511
+
type risk = Critical | High | Watch | Low
512
+
513
+
let risk_level a =
514
+
if a.pc > 1e-4 then Critical
515
+
else if a.pc > 1e-5 then High
516
+
else if a.pc > 1e-7 then Watch
517
+
else Low
518
+
519
+
let pp_risk ppf = function
520
+
| Critical -> Fmt.pf ppf "critical"
521
+
| High -> Fmt.pf ppf "high"
522
+
| Watch -> Fmt.pf ppf "watch"
523
+
| Low -> Fmt.pf ppf "low"
524
+
525
+
let pp_assessment ppf a =
526
+
Fmt.pf ppf "@[<v>Pc: %.2e (%a)@,Pc_max: %.2e@,miss: %.0f m@,vrel: %.3f km/s@]"
527
+
a.pc pp_risk (risk_level a) a.pc_max a.miss_distance a.relative_velocity
528
+
489
529
(* {1 Pretty-printing} *)
490
530
491
531
let pp_sym2 ppf s =
+59
-48
lib/collision.mli
+59
-48
lib/collision.mli
···
1
1
(** Collision probability computation for conjunction assessment.
2
2
3
-
Computes the probability of collision (Pc) between two space objects given
4
-
their state vectors, position covariances, and hard-body radii.
3
+
{2 Quick start}
5
4
6
-
Implements the standard 2D Foster method (NASA CARA), the Alfano maximum Pc
7
-
bound, and the Chan series expansion.
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"
9
+
a.pc a.miss_distance (Collision.risk_level a)
10
+
]}
8
11
9
12
{b References}:
10
13
- Foster, J.L., Estes, H.S., "A Parametric Analysis of Orbital Debris
···
13
16
- Alfano, S., "A Numerical Implementation of Spherical Object Collision
14
17
Probability", JOTA 2005 *)
15
18
16
-
(** {1 Types} *)
19
+
(** {1 Assessment} *)
17
20
18
-
type vec3 = Vec3.t = { x : float; y : float; z : float }
19
-
20
-
type encounter = {
21
-
miss_x : float; (** Miss distance x-component in conjunction plane (km). *)
22
-
miss_y : float; (** Miss distance y-component in conjunction plane (km). *)
23
-
sigma_x : float; (** 1-sigma uncertainty along x-axis (km). *)
24
-
sigma_y : float; (** 1-sigma uncertainty along y-axis (km). *)
25
-
hbr : float; (** Combined hard-body radius (km). *)
21
+
type assessment = {
22
+
pc : float; (** Probability of collision (0-1). *)
23
+
pc_max : float; (** Upper bound on Pc. *)
24
+
miss_distance : float; (** Miss distance in meters. *)
25
+
relative_velocity : float; (** Relative velocity in km/s. *)
26
26
}
27
-
(** Conjunction geometry projected onto the conjunction plane. *)
27
+
(** Summary of conjunction risk. *)
28
28
29
-
(** {1 Conjunction Plane Projection} *)
29
+
val assess : ?hbr:float -> Cdm.t -> assessment
30
+
(** [assess cdm] computes all risk metrics from a CDM.
30
31
31
-
val encounter_of_cdm : hbr:float -> Cdm.t -> encounter
32
-
(** [encounter_of_cdm ~hbr cdm] projects a CDM's conjunction geometry onto the
33
-
conjunction plane (perpendicular to relative velocity at TCA).
32
+
[hbr] is the combined hard-body radius in km (default: 0.015 km = 15 m,
33
+
typical for two LEO spacecraft). *)
34
34
35
-
The combined covariance C = C1 + C2 is projected onto the conjunction plane
36
-
and diagonalized. [hbr] is the combined hard-body radius in km (typically
37
-
0.005-0.020 km for LEO satellites). *)
35
+
type risk = Critical | High | Watch | Low
36
+
(** Risk level based on Pc thresholds (NASA CARA standard):
37
+
- [Critical]: Pc > 1e-4 — maneuver required
38
+
- [High]: Pc > 1e-5 — active monitoring, maneuver likely
39
+
- [Watch]: Pc > 1e-7 — monitor conjunction evolution
40
+
- [Low]: Pc <= 1e-7 — no action needed *)
38
41
39
-
(** {1 Probability of Collision} *)
42
+
val risk_level : assessment -> risk
43
+
(** [risk_level a] classifies the conjunction risk. *)
40
44
41
-
val pc_foster : ?n:int -> encounter -> float
42
-
(** [pc_foster ?n enc] computes the 2D probability of collision using the Foster
43
-
method (numerical integration via Gauss-Legendre quadrature).
45
+
val pp_risk : risk Fmt.t
46
+
(** Pretty-print a risk level. *)
44
47
45
-
[n] is the number of quadrature points per dimension (default 64). This is
46
-
the standard method used by NASA CARA and 18th SDS. *)
47
-
48
-
val pc_chan : ?terms:int -> encounter -> float
49
-
(** [pc_chan ?terms enc] computes Pc using the Chan series expansion. [terms] is
50
-
the number of series terms (default 50). *)
51
-
52
-
val pc_max : encounter -> float
53
-
(** [pc_max enc] returns an upper bound on Pc: [Pc_max = hbr² / (2·σx·σy)].
54
-
55
-
This bounds the integral by the peak density times the disk area. If Pc_max
56
-
is below threshold, the conjunction can be dismissed without computing the
57
-
exact Pc. *)
58
-
59
-
val pc : hbr:float -> Cdm.t -> float
60
-
(** [pc ~hbr cdm] is a convenience function that projects [cdm] onto the
61
-
conjunction plane and computes Pc using the Foster method. *)
48
+
val pp_assessment : assessment Fmt.t
49
+
(** Pretty-print an assessment. *)
62
50
63
51
(** {1 TCA refinement} *)
64
52
···
69
57
}
70
58
(** Time of Closest Approach result. *)
71
59
60
+
type vec3 = Vec3.t = { x : float; y : float; z : float }
61
+
72
62
val find_tca :
73
63
times:float array ->
74
64
pos1:vec3 array ->
···
76
66
vel1:vec3 array ->
77
67
vel2:vec3 array ->
78
68
tca option
79
-
(** [find_tca ~times ~pos1 ~pos2 ~vel1 ~vel2] finds the time of closest approach
80
-
between two objects using quadratic refinement. Arrays must be the same
81
-
length (synchronized time steps). Returns [None] if fewer than 3 points. The
82
-
refined TCA is found by fitting a parabola to d²(t) around the step-based
83
-
minimum. *)
69
+
(** [find_tca ~times ~pos1 ~pos2 ~vel1 ~vel2] finds the time of closest
70
+
approach between two objects using quadratic refinement. Arrays must be the
71
+
same length (synchronized time steps). Returns [None] if fewer than 3
72
+
points. *)
73
+
74
+
(** {1 Encounter geometry} *)
75
+
76
+
type encounter = {
77
+
miss_x : float; (** Miss distance x-component in conjunction plane (km). *)
78
+
miss_y : float; (** Miss distance y-component in conjunction plane (km). *)
79
+
sigma_x : float; (** 1-sigma uncertainty along x-axis (km). *)
80
+
sigma_y : float; (** 1-sigma uncertainty along y-axis (km). *)
81
+
hbr : float; (** Combined hard-body radius (km). *)
82
+
}
83
+
(** Conjunction geometry projected onto the conjunction plane. Most users
84
+
should use {!assess} instead of constructing this directly. *)
85
+
86
+
val encounter_of_cdm : hbr:float -> Cdm.t -> encounter
87
+
(** Project a CDM onto the conjunction plane. *)
88
+
89
+
val pc_foster : ?n:int -> encounter -> float
90
+
(** Foster 2D Pc via Gauss-Legendre quadrature (NASA CARA standard). *)
84
91
85
-
(** {1 Pretty-printing} *)
92
+
val pc_chan : ?terms:int -> encounter -> float
93
+
(** Chan series expansion Pc. *)
94
+
95
+
val pc_max : encounter -> float
96
+
(** Upper bound: [hbr^2 / (2 * sigma_x * sigma_y)]. *)
86
97
87
98
val pp_encounter : encounter Fmt.t
88
99
(** Pretty-print encounter geometry. *)