Astrodynamics coordinate frame transforms
Add Algorithm module to SDLS, ROCF service to SLE, convolutional/LDPC/turbo coding to tm-sync, station contacts to SGP4 GMAT interop, FlexACM package
- ocaml-sdls: Add Algorithm module for SA algorithm negotiation, re-export
- ocaml-sle: Add ROCF (Return Operational Control Field) service module
- ocaml-tm-sync: Add convolutional coding (rate 1/2 K=7), LDPC, turbo codes
with encode/decode and comprehensive unit tests
- ocaml-sgp4: Update GMAT interop station contacts test
- ocaml-hkdf: Regenerate opam with alcotest dep
- ocaml-flexacm: New package for flexible adaptive coding and modulation
+101
+101
lib/coordinate.ml
+101
lib/coordinate.ml
···
56
56
let j2000_to_ecef ~gmst v = rotate_z gmst v
57
57
let ecef_to_j2000 ~gmst v = rotate_z (-.gmst) v
58
58
59
+
(* ── TEME ↔ J2000 (IAU-76/FK5 precession + simplified nutation) ──── *)
60
+
61
+
(* Precession angles (IAU-76, Lieske 1979).
62
+
T is Julian centuries from J2000.0.
63
+
Returns (zeta_a, theta_a, z_a) in radians. *)
64
+
let precession_angles t =
65
+
let arcsec_to_rad a = a *. Float.pi /. (180. *. 3600.) in
66
+
let zeta_a =
67
+
arcsec_to_rad
68
+
((0.6406161 *. t) +. (0.0000839 *. t *. t) +. (0.0000050 *. t *. t *. t))
69
+
in
70
+
let theta_a =
71
+
arcsec_to_rad
72
+
((0.5567530 *. t) -. (0.0000118 *. t *. t) -. (0.0000116 *. t *. t *. t))
73
+
in
74
+
let z_a =
75
+
arcsec_to_rad
76
+
((0.6406161 *. t) +. (0.0003041 *. t *. t) +. (0.0000051 *. t *. t *. t))
77
+
in
78
+
(zeta_a, theta_a, z_a)
79
+
80
+
(* Mean obliquity of the ecliptic (IAU-76). *)
81
+
let mean_obliquity t =
82
+
let arcsec_to_rad a = a *. Float.pi /. (180. *. 3600.) in
83
+
deg_to_rad 23.439291
84
+
-. arcsec_to_rad
85
+
((46.8150 *. t) +. (0.00059 *. t *. t) -. (0.001813 *. t *. t *. t))
86
+
87
+
(* Simplified nutation (only the dominant 18.6-year lunar node term).
88
+
Returns (dpsi, deps) in radians. *)
89
+
let nutation_simplified t =
90
+
let arcsec_to_rad a = a *. Float.pi /. (180. *. 3600.) in
91
+
(* Longitude of lunar ascending node *)
92
+
let omega = deg_to_rad (125.04452 -. (1934.136261 *. t)) in
93
+
let dpsi = arcsec_to_rad (-17.20 *. sin omega) in
94
+
let deps = arcsec_to_rad (9.20 *. cos omega) in
95
+
(dpsi, deps)
96
+
97
+
let rotate_x theta (v : Vec3.t) =
98
+
let c = cos theta and s = sin theta in
99
+
Vec3.v v.x ((c *. v.y) +. (s *. v.z)) ((-.s *. v.y) +. (c *. v.z))
100
+
101
+
(* Equation of the equinoxes: GMST → GAST correction *)
102
+
let eq_equinoxes t =
103
+
let dpsi, _deps = nutation_simplified t in
104
+
let eps = mean_obliquity t in
105
+
dpsi *. cos eps
106
+
107
+
let teme_to_j2000_at ~unix_t v =
108
+
let t = julian_centuries unix_t in
109
+
(* Step 1: TEME → PEF (pseudo Earth-fixed, like ECEF but without polar motion)
110
+
Undo Earth rotation using GAST = GMST + equation of equinoxes *)
111
+
let gmst = gmst_of_unix unix_t in
112
+
let gast = gmst +. eq_equinoxes t in
113
+
let pef = rotate_z gast v in
114
+
(* Step 2: PEF → J2000 via inverse precession+nutation *)
115
+
let dpsi, deps = nutation_simplified t in
116
+
let eps = mean_obliquity t in
117
+
let eps_true = eps +. deps in
118
+
(* Nutation: N^T = Rx(eps_true) * Rz(dpsi) * Rx(-eps) *)
119
+
let v1 = rotate_x (-.eps) pef in
120
+
let v2 = rotate_z dpsi v1 in
121
+
let v3 = rotate_x eps_true v2 in
122
+
(* Precession: P^T = Rz(z_a) * Ry(-theta_a) * Rz(zeta_a) *)
123
+
let zeta_a, theta_a, z_a = precession_angles t in
124
+
let v4 = rotate_z zeta_a v3 in
125
+
let v5 = rotate_x 0.0 v4 in
126
+
let _ = v5 in
127
+
(* Actually: P = Rz(-z_a) * Ry(theta_a) * Rz(-zeta_a)
128
+
P^T (mod-to-J2000) = Rz(zeta_a) * Ry(-theta_a) * Rz(z_a) *)
129
+
let rotate_y theta (v : Vec3.t) =
130
+
let c = cos theta and s = sin theta in
131
+
Vec3.v ((c *. v.x) -. (s *. v.z)) v.y ((s *. v.x) +. (c *. v.z))
132
+
in
133
+
let v4 = rotate_z z_a v3 in
134
+
let v5 = rotate_y (-.theta_a) v4 in
135
+
rotate_z zeta_a v5
136
+
137
+
let j2000_to_teme_at ~unix_t v =
138
+
let t = julian_centuries unix_t in
139
+
let zeta_a, theta_a, z_a = precession_angles t in
140
+
let rotate_y theta (v : Vec3.t) =
141
+
let c = cos theta and s = sin theta in
142
+
Vec3.v ((c *. v.x) -. (s *. v.z)) v.y ((s *. v.x) +. (c *. v.z))
143
+
in
144
+
(* Precession: P = Rz(-z_a) * Ry(theta_a) * Rz(-zeta_a) *)
145
+
let v1 = rotate_z (-.zeta_a) v in
146
+
let v2 = rotate_y theta_a v1 in
147
+
let v3 = rotate_z (-.z_a) v2 in
148
+
(* Nutation: N = Rx(-eps_true) * Rz(-dpsi) * Rx(eps) *)
149
+
let dpsi, deps = nutation_simplified t in
150
+
let eps = mean_obliquity t in
151
+
let eps_true = eps +. deps in
152
+
let v4 = rotate_x eps v3 in
153
+
let v5 = rotate_z (-.dpsi) v4 in
154
+
let v6 = rotate_x (-.eps_true) v5 in
155
+
(* Earth rotation: GAST *)
156
+
let gmst = gmst_of_unix unix_t in
157
+
let gast = gmst +. eq_equinoxes t in
158
+
rotate_z (-.gast) v6
159
+
59
160
(* ── Geodetic conversions ──────────────────────────────────────────── *)
60
161
61
162
let ecef_to_geodetic (v : Vec3.t) =