Two-body Keplerian orbit propagation
Deduplicate math: remove pi/twopi aliases, add Coordinate.deg_to_rad
- Remove `let pi = Float.pi` aliases (collision, geometry, sphere,
spacedata) — use Float.pi directly
- Remove `let twopi` aliases (coordinate, kepler/analytic) — inline
- Add Coordinate.deg_to_rad and rad_to_deg, replace inline conversions
in coordinate, heatmap, geometry
- mu already exported as Coordinate.earth_mu and Propagate.mu
+8
-9
+8
-9
lib/analytic.ml
+8
-9
lib/analytic.ml
···
15
15
- Curtis, "Orbital Mechanics for Engineering Students", Ch. 4 *)
16
16
17
17
let mu = Propagate.mu
18
-
let twopi = 2. *. Float.pi
19
18
20
19
(** {1 Orbital elements} *)
21
20
···
80
79
else
81
80
let cos_raan = n_vec.x /. n_mag in
82
81
let raan = acos (Float.max (-1.) (Float.min 1. cos_raan)) in
83
-
if n_vec.y < 0. then twopi -. raan else raan
82
+
if n_vec.y < 0. then (2. *. Float.pi) -. raan else raan
84
83
in
85
84
(* Argument of periapsis: cos(ω) = n·e / (|n||e|) *)
86
85
let argp =
···
88
87
else
89
88
let cos_argp = dot n_vec e_vec /. (n_mag *. e) in
90
89
let argp = acos (Float.max (-1.) (Float.min 1. cos_argp)) in
91
-
if e_vec.z < 0. then twopi -. argp else argp
90
+
if e_vec.z < 0. then (2. *. Float.pi) -. argp else argp
92
91
in
93
92
(* True anomaly: cos(ν) = e·r / (|e||r|) *)
94
93
let nu =
···
96
95
if n_mag < 1e-10 then
97
96
(* Equatorial circular: use true longitude from X axis *)
98
97
let l = atan2 pos.y pos.x in
99
-
if l < 0. then l +. twopi else l
98
+
if l < 0. then l +. (2. *. Float.pi) else l
100
99
else
101
100
(* Circular inclined: use argument of latitude *)
102
101
let cos_u = dot n_vec pos /. (n_mag *. r_mag) in
103
102
let u = acos (Float.max (-1.) (Float.min 1. cos_u)) in
104
-
if pos.z < 0. then twopi -. u else u
103
+
if pos.z < 0. then (2. *. Float.pi) -. u else u
105
104
else
106
105
let cos_nu = dot e_vec pos /. (e *. r_mag) in
107
106
let nu = acos (Float.max (-1.) (Float.min 1. cos_nu)) in
108
-
if rdotv < 0. then twopi -. nu else nu
107
+
if rdotv < 0. then (2. *. Float.pi) -. nu else nu
109
108
in
110
109
(* Mean motion *)
111
110
let n =
···
120
119
(** Solve M = E - e·sin(E) for E given M and e. Newton-Raphson iteration,
121
120
converges in 3-5 steps for e < 0.97. Vallado Algorithm 2. *)
122
121
let solve_kepler ~mean_anomaly ~eccentricity =
123
-
let m = Float.rem mean_anomaly twopi in
124
-
let m = if m < 0. then m +. twopi else m in
122
+
let m = Float.rem mean_anomaly (2. *. Float.pi) in
123
+
let m = if m < 0. then m +. (2. *. Float.pi) else m in
125
124
let e = eccentricity in
126
125
(* Initial guess: E₀ = M + e·sin(M) for low eccentricity *)
127
126
let ea = ref (m +. (e *. sin m)) in
···
175
174
176
175
(** {1 Element queries} *)
177
176
178
-
let period el = if el.a <= 0. then infinity (* hyperbolic *) else twopi /. el.n
177
+
let period el = if el.a <= 0. then infinity (* hyperbolic *) else (2. *. Float.pi) /. el.n
179
178
let eccentricity el = el.e
180
179
let semi_major_axis el = el.a
181
180
let inclination el = el.i