Two-body Keplerian orbit propagation
Fix merlint E005/E010: split long functions and flatten deep nesting
Refactored 21 functions across 17 files to bring them under merlint
thresholds. Pure structural refactoring — no behavior changes.
Packages: irmin, merlint, monopam, ocaml-btree, ocaml-cdm, ocaml-git,
ocaml-globe, ocaml-kepler, ocaml-osv, ocaml-sigstore, ocaml-stix, prune.
+27
-22
+27
-22
lib/analytic.ml
+27
-22
lib/analytic.ml
···
32
32
33
33
(** {1 State vector ↔ orbital elements} *)
34
34
35
+
(** Clamp [x] to [-1, 1] before [acos]. *)
36
+
let clamp_acos x = acos (Float.max (-1.) (Float.min 1. x))
37
+
38
+
(** Compute true anomaly from eccentricity vector, node vector, position, and
39
+
radial-velocity sign. Handles circular/equatorial special cases. *)
40
+
let compute_true_anomaly ~(e_vec : Vec3.t) ~e ~(n_vec : Vec3.t) ~n_mag
41
+
~(pos : Vec3.t) ~r_mag ~rdotv =
42
+
let open Vec3 in
43
+
if e < 1e-10 then
44
+
if n_mag < 1e-10 then
45
+
(* Equatorial circular: use true longitude from X axis *)
46
+
let l = atan2 pos.y pos.x in
47
+
if l < 0. then l +. (2. *. Float.pi) else l
48
+
else
49
+
(* Circular inclined: use argument of latitude *)
50
+
let cos_u = dot n_vec pos /. (n_mag *. r_mag) in
51
+
let u = clamp_acos cos_u in
52
+
if pos.z < 0. then (2. *. Float.pi) -. u else u
53
+
else
54
+
let cos_nu = dot e_vec pos /. (e *. r_mag) in
55
+
let nu = clamp_acos cos_nu in
56
+
if rdotv < 0. then (2. *. Float.pi) -. nu else nu
57
+
35
58
(** Convert state vector to classical orbital elements. Vallado Algorithm 9. *)
36
59
let elements_of_state (pos : Vec3.t) (vel : Vec3.t) =
37
60
let open Vec3 in
···
72
95
else -.mu /. (2. *. energy)
73
96
in
74
97
(* Inclination: cos(i) = hz / |h| *)
75
-
let i = acos (Float.max (-1.) (Float.min 1. (h.z /. h_mag))) in
98
+
let i = clamp_acos (h.z /. h_mag) in
76
99
(* RAAN: cos(Ω) = nx / |n| *)
77
100
let raan =
78
101
if n_mag < 1e-10 then 0.
79
102
else
80
-
let cos_raan = n_vec.x /. n_mag in
81
-
let raan = acos (Float.max (-1.) (Float.min 1. cos_raan)) in
103
+
let raan = clamp_acos (n_vec.x /. n_mag) in
82
104
if n_vec.y < 0. then (2. *. Float.pi) -. raan else raan
83
105
in
84
106
(* Argument of periapsis: cos(ω) = n·e / (|n||e|) *)
85
107
let argp =
86
108
if n_mag < 1e-10 || e < 1e-10 then 0.
87
109
else
88
-
let cos_argp = dot n_vec e_vec /. (n_mag *. e) in
89
-
let argp = acos (Float.max (-1.) (Float.min 1. cos_argp)) in
110
+
let argp = clamp_acos (dot n_vec e_vec /. (n_mag *. e)) in
90
111
if e_vec.z < 0. then (2. *. Float.pi) -. argp else argp
91
112
in
92
-
(* True anomaly: cos(ν) = e·r / (|e||r|) *)
93
-
let nu =
94
-
if e < 1e-10 then
95
-
if n_mag < 1e-10 then
96
-
(* Equatorial circular: use true longitude from X axis *)
97
-
let l = atan2 pos.y pos.x in
98
-
if l < 0. then l +. (2. *. Float.pi) else l
99
-
else
100
-
(* Circular inclined: use argument of latitude *)
101
-
let cos_u = dot n_vec pos /. (n_mag *. r_mag) in
102
-
let u = acos (Float.max (-1.) (Float.min 1. cos_u)) in
103
-
if pos.z < 0. then (2. *. Float.pi) -. u else u
104
-
else
105
-
let cos_nu = dot e_vec pos /. (e *. r_mag) in
106
-
let nu = acos (Float.max (-1.) (Float.min 1. cos_nu)) in
107
-
if rdotv < 0. then (2. *. Float.pi) -. nu else nu
108
-
in
113
+
let nu = compute_true_anomaly ~e_vec ~e ~n_vec ~n_mag ~pos ~r_mag ~rdotv in
109
114
(* Mean motion *)
110
115
let n =
111
116
if a > 0. then sqrt (mu /. (a *. a *. a))