Astrodynamics coordinate frame transforms
chore: stage linter changes across packages
+22
-116
+13
-53
lib/coordinate.ml
+13
-53
lib/coordinate.ml
···
5
5
6
6
(** Astrodynamics coordinate frame transforms. *)
7
7
8
-
(* ------------------------------------------------------------------ *)
9
-
(* Types *)
10
-
(* ------------------------------------------------------------------ *)
11
-
12
-
type vec3 = { x : float; y : float; z : float }
13
8
type geodetic = { lat : float; lon : float; alt : float }
14
9
15
-
(* ------------------------------------------------------------------ *)
16
-
(* Constants (WGS-84) *)
17
-
(* ------------------------------------------------------------------ *)
10
+
(* ── Constants (WGS-84) ───────────────────────────────────────────── *)
18
11
19
12
let earth_radius = 6378.137
20
13
let earth_flattening = 1. /. 298.257223563
···
23
16
let earth_rotation_rate = 7.2921151467e-5
24
17
let twopi = 2. *. Float.pi
25
18
26
-
(* ------------------------------------------------------------------ *)
27
-
(* Utility *)
28
-
(* ------------------------------------------------------------------ *)
29
-
30
-
let vec3 x y z = { x; y; z }
31
-
let vec3_length v = sqrt ((v.x *. v.x) +. (v.y *. v.y) +. (v.z *. v.z))
19
+
(* ── Utility ───────────────────────────────────────────────────────── *)
32
20
33
21
let normalize_angle a =
34
22
let a = Float.rem a twopi in
35
23
if a < 0. then a +. twopi else a
36
24
37
-
(* ------------------------------------------------------------------ *)
38
-
(* Time conversions *)
39
-
(* ------------------------------------------------------------------ *)
25
+
(* ── Time ──────────────────────────────────────────────────────────── *)
40
26
41
27
let julian_date_of_unix t = (t /. 86400.) +. 2440587.5
42
28
let unix_of_julian_date jd = (jd -. 2440587.5) *. 86400.
···
45
31
let jd = julian_date_of_unix unix_t in
46
32
(jd -. 2451545.0) /. 36525.0
47
33
48
-
(** GMST using the IAU 1982 model.
49
-
50
-
Reference: Vallado, Eq. 3-47. GMST = 67310.54841 + (876600h +
51
-
8640184.812866)T
52
-
+ 0.093104 T^2 - 6.2e-6 T^3 where T is Julian centuries from J2000.0. *)
53
34
let gmst_of_unix unix_time =
54
35
let t = julian_centuries unix_time in
55
36
let gmst_sec =
···
62
43
63
44
let gmst ptime = gmst_of_unix (Ptime.to_float_s ptime)
64
45
65
-
(* ------------------------------------------------------------------ *)
66
-
(* Frame transforms *)
67
-
(* ------------------------------------------------------------------ *)
46
+
(* ── Frame transforms ──────────────────────────────────────────────── *)
68
47
69
-
(** Z-axis rotation by angle theta. Rz(θ) =
70
-
[[cos θ, sin θ, 0], [-sin θ, cos θ, 0], [0, 0, 1]] *)
71
-
let rotate_z theta v =
48
+
let rotate_z theta (v : Vec3.t) =
72
49
let c = cos theta and s = sin theta in
73
-
{ x = (c *. v.x) +. (s *. v.y); y = (-.s *. v.x) +. (c *. v.y); z = v.z }
50
+
Vec3.v ((c *. v.x) +. (s *. v.y)) ((-.s *. v.x) +. (c *. v.y)) v.z
74
51
75
52
let teme_to_ecef ~gmst v = rotate_z gmst v
76
53
let ecef_to_teme ~gmst v = rotate_z (-.gmst) v
77
54
let j2000_to_ecef ~gmst v = rotate_z gmst v
78
55
let ecef_to_j2000 ~gmst v = rotate_z (-.gmst) v
79
56
80
-
(* ------------------------------------------------------------------ *)
81
-
(* Geodetic conversions (WGS-84 ellipsoid) *)
82
-
(* ------------------------------------------------------------------ *)
57
+
(* ── Geodetic conversions ──────────────────────────────────────────── *)
83
58
84
-
(** Bowring's iterative method for ECEF → geodetic.
85
-
86
-
Reference: Vallado Algorithm 12, Montenbruck & Gill §5.3. Converges in 2-3
87
-
iterations for typical positions. *)
88
-
let ecef_to_geodetic v =
59
+
let ecef_to_geodetic (v : Vec3.t) =
89
60
let a = earth_radius in
90
61
let f = earth_flattening in
91
62
let e2 = (2. *. f) -. (f *. f) in
92
-
(* first eccentricity squared *)
93
63
let p = sqrt ((v.x *. v.x) +. (v.y *. v.y)) in
94
64
let lon_rad = atan2 v.y v.x in
95
-
(* Initial latitude estimate *)
96
65
let lat_rad = ref (atan2 v.z p) in
97
-
(* Iterate Bowring's method *)
98
66
for _ = 1 to 5 do
99
67
let sin_lat = sin !lat_rad in
100
68
let n = a /. sqrt (1. -. (e2 *. sin_lat *. sin_lat)) in
···
108
76
in
109
77
{ lat = !lat_rad *. 180. /. Float.pi; lon = lon_rad *. 180. /. Float.pi; alt }
110
78
111
-
(** Geodetic → ECEF.
112
-
113
-
Reference: Vallado Algorithm 11. *)
114
79
let geodetic_to_ecef g =
115
80
let a = earth_radius in
116
81
let f = earth_flattening in
···
120
85
let sin_lat = sin lat_rad in
121
86
let cos_lat = cos lat_rad in
122
87
let n = a /. sqrt (1. -. (e2 *. sin_lat *. sin_lat)) in
123
-
{
124
-
x = (n +. g.alt) *. cos_lat *. cos lon_rad;
125
-
y = (n +. g.alt) *. cos_lat *. sin lon_rad;
126
-
z = ((n *. (1. -. e2)) +. g.alt) *. sin_lat;
127
-
}
88
+
Vec3.v
89
+
((n +. g.alt) *. cos_lat *. cos lon_rad)
90
+
((n +. g.alt) *. cos_lat *. sin lon_rad)
91
+
(((n *. (1. -. e2)) +. g.alt) *. sin_lat)
128
92
129
93
let teme_to_geodetic ~gmst v = ecef_to_geodetic (teme_to_ecef ~gmst v)
130
94
131
-
(* ------------------------------------------------------------------ *)
132
-
(* Pretty-printing *)
133
-
(* ------------------------------------------------------------------ *)
134
-
135
-
let pp_vec3 ppf v = Fmt.pf ppf "(%.3f, %.3f, %.3f)" v.x v.y v.z
95
+
(* ── Pretty-printing ───────────────────────────────────────────────── *)
136
96
137
97
let pp_geodetic ppf g =
138
98
Fmt.pf ppf "%.4f°%s %.4f°%s %.3f km" (abs_float g.lat)
+8
-62
lib/coordinate.mli
+8
-62
lib/coordinate.mli
···
15
15
16
16
(** {1 Types} *)
17
17
18
-
type vec3 = { x : float; y : float; z : float }
19
-
(** Cartesian position or velocity vector (km or km/s). *)
20
-
21
18
type geodetic = { lat : float; lon : float; alt : float }
22
-
(** Geodetic coordinates: latitude (deg), longitude (deg), altitude (km).
23
-
Latitude: [-90, +90], north positive. Longitude: [-180, +180], east
24
-
positive. *)
19
+
(** Geodetic coordinates: latitude (deg), longitude (deg), altitude (km). *)
25
20
26
21
(** {1 Constants} *)
27
22
28
23
val earth_radius : float
29
-
(** [earth_radius] is Earth's mean equatorial radius (6378.137 km, WGS-84). *)
30
-
31
24
val earth_flattening : float
32
-
(** [earth_flattening] is Earth's flattening factor (1/298.257223563, WGS-84).
33
-
*)
34
-
35
25
val earth_mu : float
36
-
(** [earth_mu] is Earth's gravitational parameter (398600.4418 km^3/s^2). *)
37
-
38
26
val earth_j2 : float
39
-
(** [earth_j2] is Earth's second zonal harmonic (0.00108263). *)
40
-
41
27
val earth_rotation_rate : float
42
-
(** [earth_rotation_rate] is Earth's rotation rate (7.2921151467e-5 rad/s). *)
43
28
44
29
(** {1 Time} *)
45
30
46
31
val gmst_of_unix : float -> float
47
-
(** [gmst_of_unix t] is Greenwich Mean Sidereal Time (radians) at Unix timestamp
48
-
[t]. Uses the IAU 1982 model (Vallado Eq. 3-47). *)
49
-
50
32
val gmst : Ptime.t -> float
51
-
(** [gmst t] is GMST (radians) at time [t]. *)
52
-
53
33
val julian_date_of_unix : float -> float
54
-
(** [julian_date_of_unix t] is the Julian Date for Unix timestamp [t]. *)
55
-
56
34
val unix_of_julian_date : float -> float
57
-
(** [unix_of_julian_date jd] is the Unix timestamp for Julian Date [jd]. *)
58
-
59
35
val julian_centuries : float -> float
60
-
(** [julian_centuries unix_t] is Julian centuries from J2000.0 epoch. *)
61
36
62
37
(** {1 Frame transforms} *)
63
38
64
-
val teme_to_ecef : gmst:float -> vec3 -> vec3
65
-
(** [teme_to_ecef ~gmst v] rotates vector [v] from TEME to ECEF frame using GMST
66
-
angle. This is a simple Z-axis rotation. (Vallado Algorithm 24). *)
67
-
68
-
val ecef_to_teme : gmst:float -> vec3 -> vec3
69
-
(** [ecef_to_teme ~gmst v] is the inverse of {!teme_to_ecef}. *)
70
-
71
-
val j2000_to_ecef : gmst:float -> vec3 -> vec3
72
-
(** [j2000_to_ecef ~gmst v] rotates from J2000/EME2000 to ECEF. At our precision
73
-
level (~arcsecond), J2000≈TEME, so this uses the same GMST rotation. *)
74
-
75
-
val ecef_to_j2000 : gmst:float -> vec3 -> vec3
76
-
(** [ecef_to_j2000 ~gmst v] is the inverse of {!j2000_to_ecef}. *)
39
+
val teme_to_ecef : gmst:float -> Vec3.t -> Vec3.t
40
+
val ecef_to_teme : gmst:float -> Vec3.t -> Vec3.t
41
+
val j2000_to_ecef : gmst:float -> Vec3.t -> Vec3.t
42
+
val ecef_to_j2000 : gmst:float -> Vec3.t -> Vec3.t
77
43
78
44
(** {1 Geodetic conversions} *)
79
45
80
-
val ecef_to_geodetic : vec3 -> geodetic
81
-
(** [ecef_to_geodetic v] converts ECEF position (km) to geodetic coordinates
82
-
using Bowring's iterative method on the WGS-84 ellipsoid. (Vallado Algorithm
83
-
12). *)
84
-
85
-
val geodetic_to_ecef : geodetic -> vec3
86
-
(** [geodetic_to_ecef g] converts geodetic coordinates to ECEF (km). (Vallado
87
-
Algorithm 11). *)
88
-
89
-
val teme_to_geodetic : gmst:float -> vec3 -> geodetic
90
-
(** [teme_to_geodetic ~gmst v] converts TEME position to geodetic. Convenience
91
-
for [ecef_to_geodetic (teme_to_ecef ~gmst v)]. *)
46
+
val ecef_to_geodetic : Vec3.t -> geodetic
47
+
val geodetic_to_ecef : geodetic -> Vec3.t
48
+
val teme_to_geodetic : gmst:float -> Vec3.t -> geodetic
92
49
93
50
(** {1 Utility} *)
94
51
95
-
val vec3 : float -> float -> float -> vec3
96
-
(** [vec3 x y z] is [{x; y; z}]. *)
97
-
98
-
val vec3_length : vec3 -> float
99
-
(** [vec3_length v] is the Euclidean length of [v]. *)
100
-
101
52
val normalize_angle : float -> float
102
-
(** [normalize_angle a] normalizes angle [a] to \[0, 2pi). *)
103
53
104
54
(** {1 Pretty-printing} *)
105
55
106
-
val pp_vec3 : vec3 Fmt.t
107
-
(** [pp_vec3] pretty-prints a vector. *)
108
-
109
56
val pp_geodetic : geodetic Fmt.t
110
-
(** [pp_geodetic] pretty-prints geodetic coordinates. *)