Astrodynamics coordinate frame transforms
1(** Coordinate transform tests.
2
3 Test vectors from:
4 - Vallado, "Fundamentals of Astrodynamics", 4th ed., Examples 3-3..3-5
5 - IAU SOFA test cases for GMST
6 - WGS-84 reference geodetic test points
7 - Montenbruck & Gill, "Satellite Orbits", Ch. 5 examples
8 - Known geographic landmarks for geodetic round-trip *)
9
10let check_float msg eps expected actual =
11 Alcotest.(check (float eps)) msg expected actual
12
13let deg_to_rad d = d *. Float.pi /. 180.
14let rad_to_deg r = r *. 180. /. Float.pi
15
16(* ================================================================== *)
17(* Constants *)
18(* ================================================================== *)
19
20let test_constants () =
21 (* WGS-84 values: ITG standard *)
22 check_float "earth_radius" 1e-3 6378.137 Coordinate.earth_radius;
23 check_float "earth_mu" 1e-3 398600.4418 Coordinate.earth_mu;
24 check_float "flattening" 1e-12 (1. /. 298.257223563)
25 Coordinate.earth_flattening;
26 check_float "J2" 1e-8 0.00108263 Coordinate.earth_j2
27
28(* ================================================================== *)
29(* GMST tests *)
30(* ================================================================== *)
31
32(** Vallado Example 3-5: GMST at J2000.0 epoch. J2000.0 = 2000 Jan 1 12:00:00 TT
33 ≈ 2000 Jan 1 11:58:55.816 UTC. GMST ≈ 280.46° = 4.8949612 rad. *)
34let test_gmst_j2000 () =
35 let j2000_unix = 946728000. in
36 (* 2000-01-01T12:00:00 UTC approx *)
37 let gmst = Coordinate.gmst_of_unix j2000_unix in
38 let gmst_deg = rad_to_deg gmst in
39 (* Allow 0.5° tolerance for TT/UTC difference *)
40 let diff = abs_float (gmst_deg -. 280.46) in
41 Alcotest.(check bool) "GMST at J2000 ~280.46°" true (diff < 0.5)
42
43(** GMST increases by ~360.985647° per day (sidereal day). After exactly one
44 sidereal day (86164.0905s), GMST should advance by exactly 2π. *)
45let test_gmst_sidereal_day () =
46 let t0 = 1735732800. in
47 (* 2025-01-01T12:00:00 UTC *)
48 let sidereal_day = 86164.0905 in
49 let g0 = Coordinate.gmst_of_unix t0 in
50 let g1 = Coordinate.gmst_of_unix (t0 +. sidereal_day) in
51 (* After one sidereal day, GMST should return to ~same value *)
52 let diff = abs_float (g1 -. g0) in
53 let diff = if diff > Float.pi then (2. *. Float.pi) -. diff else diff in
54 Alcotest.(check bool) "GMST period" true (diff < 0.001)
55
56(** GMST should be monotonically increasing (mod 2π). *)
57let test_gmst_monotone () =
58 let t0 = 946728000. in
59 let dt = 3600. in
60 (* 1 hour steps *)
61 let prev = ref (Coordinate.gmst_of_unix t0) in
62 for i = 1 to 24 do
63 let t = t0 +. (Float.of_int i *. dt) in
64 let g = Coordinate.gmst_of_unix t in
65 let delta = g -. !prev in
66 (* Allow wrap-around at 2π *)
67 let delta = if delta < -1. then delta +. (2. *. Float.pi) else delta in
68 Alcotest.(check bool) (Fmt.str "monotone step %d" i) true (delta > 0.);
69 prev := g
70 done
71
72(** Ptime-based GMST should match unix-based. *)
73let test_gmst_ptime () =
74 let unix_t = 1735732800. in
75 let ptime =
76 match Ptime.of_float_s unix_t with Some t -> t | None -> Ptime.epoch
77 in
78 let g1 = Coordinate.gmst_of_unix unix_t in
79 let g2 = Coordinate.gmst ptime in
80 check_float "ptime vs unix" 1e-6 g1 g2
81
82(* ================================================================== *)
83(* Julian date tests *)
84(* ================================================================== *)
85
86(** J2000.0 epoch: JD 2451545.0 = 2000-01-01T12:00:00 UTC. *)
87let test_jd_j2000 () =
88 let jd = Coordinate.julian_date_of_unix 946728000. in
89 check_float "JD at J2000" 0.01 2451545.0 jd
90
91(** Unix epoch: JD 2440587.5 = 1970-01-01T00:00:00 UTC. *)
92let test_jd_unix_epoch () =
93 let jd = Coordinate.julian_date_of_unix 0. in
94 check_float "JD at unix epoch" 1e-6 2440587.5 jd
95
96(** Julian date round-trip. *)
97let test_jd_roundtrip () =
98 let t = 1735732800. in
99 let jd = Coordinate.julian_date_of_unix t in
100 let t' = Coordinate.unix_of_julian_date jd in
101 check_float "JD round-trip" 1e-3 t t'
102
103(* ================================================================== *)
104(* TEME ↔ ECEF tests *)
105(* ================================================================== *)
106
107(** At GMST=0, TEME and ECEF are aligned. *)
108let test_teme_ecef_zero () =
109 let v = Vec3.v 1000. 2000. 3000. in
110 let ecef = Coordinate.teme_to_ecef ~gmst:0. v in
111 check_float "x" 1e-6 1000. ecef.x;
112 check_float "y" 1e-6 2000. ecef.y;
113 check_float "z" 1e-6 3000. ecef.z
114
115(** At GMST=π/2: Rz(π/2)*(1000,0,0) = (0, -1000, 0). ECEF X' = cos θ · X + sin θ
116 · Y = 0 ECEF Y' = -sin θ · X + cos θ · Y = -1000 *)
117let test_teme_ecef_90deg () =
118 let v = Vec3.v 1000. 0. 0. in
119 let gmst = deg_to_rad 90. in
120 let ecef = Coordinate.teme_to_ecef ~gmst v in
121 check_float "x" 1e-6 0. ecef.x;
122 check_float "y" 1e-6 (-1000.) ecef.y;
123 check_float "z" 1e-6 0. ecef.z
124
125(** Round-trip: TEME → ECEF → TEME should be identity. *)
126let test_teme_ecef_roundtrip () =
127 let v = Vec3.v 6524.834 6862.875 6448.296 in
128 let gmst = 1.7 in
129 let ecef = Coordinate.teme_to_ecef ~gmst v in
130 let v' = Coordinate.ecef_to_teme ~gmst ecef in
131 check_float "rt x" 1e-6 v.x v'.x;
132 check_float "rt y" 1e-6 v.y v'.y;
133 check_float "rt z" 1e-6 v.z v'.z
134
135(** Rotation preserves vector length. *)
136let test_rotation_preserves_length () =
137 let v = Vec3.v 7000. 1000. 3000. in
138 let r0 = Vec3.length v in
139 let ecef = Coordinate.teme_to_ecef ~gmst:2.345 v in
140 let r1 = Vec3.length ecef in
141 check_float "length preserved" 1e-6 r0 r1
142
143(* ================================================================== *)
144(* Geodetic tests *)
145(* ================================================================== *)
146
147(** Point on equator at prime meridian: lat=0, lon=0, alt=0. ECEF = (a, 0, 0)
148 where a = earth_radius. *)
149let test_geodetic_equator_pm () =
150 let v = Vec3.v Coordinate.earth_radius 0. 0. in
151 let g = Coordinate.ecef_to_geodetic v in
152 check_float "lat" 0.01 0. g.lat;
153 check_float "lon" 0.01 0. g.lon;
154 check_float "alt" 0.01 0. g.alt
155
156(** North pole: lat=90, lon=0 (or undefined), alt=0. ECEF = (0, 0, b) where b =
157 a*(1-f). *)
158let test_geodetic_north_pole () =
159 let b = Coordinate.earth_radius *. (1. -. Coordinate.earth_flattening) in
160 let v = Vec3.v 0. 0. b in
161 let g = Coordinate.ecef_to_geodetic v in
162 check_float "north pole lat" 0.01 90. g.lat;
163 check_float "north pole alt" 1. 0. g.alt
164
165(** Point at lon=90°E on equator: ECEF = (0, a, 0). *)
166let test_geodetic_lon90 () =
167 let v = Vec3.v 0. Coordinate.earth_radius 0. in
168 let g = Coordinate.ecef_to_geodetic v in
169 check_float "lat" 0.01 0. g.lat;
170 check_float "lon" 0.01 90. g.lon
171
172(** Known landmark: Washington Monument (Vallado Example 3-3). Geodetic:
173 38.8895°N, 77.0353°W, 0.054 km alt. ECEF: ~(1115.479, -4843.992, 3984.170)
174 km. *)
175let test_geodetic_washington () =
176 let g = { Coordinate.lat = 38.8895; lon = -77.0353; alt = 0.054 } in
177 let ecef = Coordinate.geodetic_to_ecef g in
178 check_float "wash x" 2. 1115.479 ecef.x;
179 check_float "wash y" 2. (-4843.992) ecef.y;
180 check_float "wash z" 2. 3984.170 ecef.z
181
182(** Geodetic round-trip: geodetic → ECEF → geodetic. *)
183let test_geodetic_roundtrip () =
184 let g = { Coordinate.lat = 48.8566; lon = 2.3522; alt = 0.035 } in
185 let ecef = Coordinate.geodetic_to_ecef g in
186 let g' = Coordinate.ecef_to_geodetic ecef in
187 check_float "rt lat" 1e-4 g.lat g'.lat;
188 check_float "rt lon" 1e-4 g.lon g'.lon;
189 check_float "rt alt" 1e-2 g.alt g'.alt
190
191(** Geodetic round-trip at various latitudes. *)
192let test_geodetic_roundtrip_lats () =
193 let lats = [ -89.99; -60.; -30.; 0.; 30.; 60.; 89.99 ] in
194 List.iter
195 (fun lat ->
196 let g = { Coordinate.lat; lon = 45.; alt = 400. } in
197 let ecef = Coordinate.geodetic_to_ecef g in
198 let g' = Coordinate.ecef_to_geodetic ecef in
199 check_float (Fmt.str "lat%.0f lat" lat) 1e-3 g.lat g'.lat;
200 check_float (Fmt.str "lat%.0f lon" lat) 1e-3 g.lon g'.lon;
201 check_float (Fmt.str "lat%.0f alt" lat) 1e-1 g.alt g'.alt)
202 lats
203
204(** ISS altitude: ~408 km. Convert TEME position to geodetic and check. *)
205let test_geodetic_iss_altitude () =
206 (* ISS at ~408 km, circular orbit, position along X axis *)
207 let r = Coordinate.earth_radius +. 408. in
208 let v = Vec3.v r 0. 0. in
209 let g = Coordinate.ecef_to_geodetic v in
210 check_float "ISS alt" 5. 408. g.alt;
211 check_float "ISS lat" 0.1 0. g.lat
212
213(** GEO altitude: ~35786 km above equator. *)
214let test_geodetic_geo_altitude () =
215 let r = 42164.0 in
216 (* GEO radius *)
217 let v = Vec3.v r 0. 0. in
218 let g = Coordinate.ecef_to_geodetic v in
219 let expected_alt = r -. Coordinate.earth_radius in
220 check_float "GEO alt" 5. expected_alt g.alt
221
222(** Vallado Example 3-3: given ECEF, compute geodetic. r_ECEF = (6524.834,
223 6862.875, 6448.296) km Expected: lat ≈ 34.35°, lon ≈ 46.45°, alt ≈ 5085 km.
224*)
225let test_vallado_3_3 () =
226 let v = Vec3.v 6524.834 6862.875 6448.296 in
227 let g = Coordinate.ecef_to_geodetic v in
228 (* Approximate — the example uses a different reference ellipsoid *)
229 let r = Vec3.length v in
230 Alcotest.(check (float 0.001)) "radius" 11456.572 r;
231 Alcotest.(check (float 0.5)) "lat" 34.25 g.lat;
232 Alcotest.(check (float 50.)) "alt" 5078. g.alt
233
234(* ================================================================== *)
235(* TEME → Geodetic (full pipeline) *)
236(* ================================================================== *)
237
238(** TEME position at equator, GMST=0, should give lon≈0. *)
239let test_teme_geodetic_equator () =
240 let r = Coordinate.earth_radius +. 400. in
241 let v = Vec3.v r 0. 0. in
242 let g = Coordinate.teme_to_geodetic ~gmst:0. v in
243 check_float "equator lat" 0.1 0. g.lat;
244 check_float "equator lon" 0.1 0. g.lon
245
246(** TEME position at equator with GMST=π/2: lon should shift by -90°. *)
247let test_teme_geodetic_gmst_shift () =
248 let r = Coordinate.earth_radius +. 400. in
249 let v = Vec3.v r 0. 0. in
250 let g0 = Coordinate.teme_to_geodetic ~gmst:0. v in
251 let g90 = Coordinate.teme_to_geodetic ~gmst:(deg_to_rad 90.) v in
252 let dlon = g0.lon -. g90.lon in
253 check_float "lon shift 90°" 1. 90. (abs_float dlon)
254
255(* ================================================================== *)
256(* Normalize angle *)
257(* ================================================================== *)
258
259let test_normalize_angle () =
260 let twopi = 2. *. Float.pi in
261 check_float "positive" 1e-10 1.5 (Coordinate.normalize_angle 1.5);
262 check_float "negative" 1e-10 (twopi -. 1.) (Coordinate.normalize_angle (-1.));
263 check_float "large" 1e-10 1.0 (Coordinate.normalize_angle (1. +. twopi));
264 check_float "zero" 1e-10 0. (Coordinate.normalize_angle 0.)
265
266(* ================================================================== *)
267(* Suite *)
268(* ================================================================== *)
269
270(* ================================================================== *)
271(* TEME <-> J2000 (precession + nutation) *)
272(* ================================================================== *)
273
274(* Reference: python computation using IAU-76 precession + IAU-1980 nutation
275 (top 2 terms). Unix timestamp 1767225600.0 = 2026-01-01T00:00:00 UTC.
276 T = 0.26 Julian centuries from J2000.0. *)
277
278let unix_2026 = 1767225600.0
279
280(* TEME->J2000 should change the position vector (not identity). *)
281let test_teme_j2000_not_identity () =
282 let teme = Vec3.v (-2974.696) 2974.696 5307.732 in
283 let j2000 = Coordinate.teme_to_j2000_at ~unix_t:unix_2026 teme in
284 let d = Vec3.distance teme j2000 in
285 (* Precession rotates by ~600 arcsec over 26 years = ~40 km at LEO *)
286 Alcotest.(check bool) "not identity" true (d > 1.0);
287 Alcotest.(check bool) "reasonable magnitude" true (d < 100.0)
288
289(* Round-trip: TEME -> J2000 -> TEME should recover the original vector. *)
290let test_teme_j2000_roundtrip () =
291 let teme = Vec3.v (-2974.696809661101) 2974.696809661102 5307.732034462389 in
292 let j2000 = Coordinate.teme_to_j2000_at ~unix_t:unix_2026 teme in
293 let back = Coordinate.j2000_to_teme_at ~unix_t:unix_2026 j2000 in
294 let err = Vec3.distance teme back in
295 (* Round-trip should be exact to float precision *)
296 if err > 1e-9 then
297 Alcotest.failf "round-trip error: %.15e km (should be < 1e-9)" err
298
299(* Reference vector from python (IAU-76 + IAU-1980 top 2 terms).
300 Source: manual computation matching ERFA nut80 fundamental arguments.
301 TEME: (-2974.696809661101, 2974.696809661102, 5307.732034462389) km
302 J2000: (-2943.795451294467, 2992.190764225030, 5315.122231886917) km *)
303let test_teme_j2000_reference () =
304 let teme = Vec3.v (-2974.696809661101) 2974.696809661102 5307.732034462389 in
305 let j2000 = Coordinate.teme_to_j2000_at ~unix_t:unix_2026 teme in
306 (* Our implementation uses 13 nutation terms vs python's 2, so there will be
307 a small difference. Allow 1 km for the nutation term difference. *)
308 let dx = j2000.x -. -2943.795451294467 in
309 let dy = j2000.y -. 2992.190764225030 in
310 let dz = j2000.z -. 5315.122231886917 in
311 let err = sqrt ((dx *. dx) +. (dy *. dy) +. (dz *. dz)) in
312 Fmt.pr " TEME->J2000 vs python reference: %.6f km\n" err;
313 if err > 1.0 then
314 Alcotest.failf "TEME->J2000 differs from python by %.6f km (> 1 km)" err
315
316(* Frame conversion preserves vector magnitude (rotation = isometry). *)
317let test_teme_j2000_preserves_length () =
318 let teme = Vec3.v (-2974.696) 2974.696 5307.732 in
319 let j2000 = Coordinate.teme_to_j2000_at ~unix_t:unix_2026 teme in
320 let r_teme = Vec3.length teme in
321 let r_j2000 = Vec3.length j2000 in
322 let err = Float.abs (r_teme -. r_j2000) in
323 if err > 1e-10 then
324 Alcotest.failf "length not preserved: %.15e vs %.15e (diff %.15e)" r_teme
325 r_j2000 err
326
327(* At J2000 epoch (T=0), precession is zero but nutation is nonzero.
328 The TEME->J2000 conversion should still be close to identity. *)
329let test_teme_j2000_at_epoch () =
330 let unix_j2000 = Coordinate.unix_of_julian_date 2451545.0 in
331 let teme = Vec3.v 6778.0 0.0 0.0 in
332 let j2000 = Coordinate.teme_to_j2000_at ~unix_t:unix_j2000 teme in
333 let d = Vec3.distance teme j2000 in
334 (* At J2000 epoch, only nutation contributes. dpsi ~ 17 arcsec max.
335 At 6778 km, 17 arcsec = 6778 * tan(17/3600 * pi/180) = ~0.56 km *)
336 Fmt.pr " TEME-J2000 at J2000 epoch: %.6f km difference\n" d;
337 if d > 2.0 then
338 Alcotest.failf "TEME-J2000 at J2000 epoch: %.6f km (should be < 2 km)" d
339
340(* Test with multiple epochs to check for sign errors in precession. *)
341let test_teme_j2000_multiple_epochs () =
342 let teme = Vec3.v 6778.0 0.0 0.0 in
343 (* At later epochs, precession should increase the rotation *)
344 let d1 =
345 Vec3.distance teme (Coordinate.teme_to_j2000_at ~unix_t:unix_2026 teme)
346 in
347 let unix_2050 = unix_2026 +. (24.0 *. 365.25 *. 86400.0) in
348 let d2 =
349 Vec3.distance teme (Coordinate.teme_to_j2000_at ~unix_t:unix_2050 teme)
350 in
351 (* 2050 is further from J2000 than 2026, so precession should be larger *)
352 Alcotest.(check bool) "precession grows with time" true (d2 > d1);
353 Fmt.pr " 2026: %.3f km, 2050: %.3f km\n" d1 d2
354
355(* Inverse: J2000 -> TEME reference. *)
356let test_j2000_teme_reference () =
357 let j2000 = Vec3.v (-2943.795) 2992.191 5315.122 in
358 let teme = Coordinate.j2000_to_teme_at ~unix_t:unix_2026 j2000 in
359 let back = Coordinate.teme_to_j2000_at ~unix_t:unix_2026 teme in
360 let err = Vec3.distance j2000 back in
361 if err > 1e-9 then Alcotest.failf "J2000->TEME->J2000 error: %.15e km" err
362
363let suite =
364 ( "coordinate",
365 [
366 Alcotest.test_case "constants" `Quick test_constants;
367 Alcotest.test_case "GMST at J2000" `Quick test_gmst_j2000;
368 Alcotest.test_case "GMST sidereal day" `Quick test_gmst_sidereal_day;
369 Alcotest.test_case "GMST monotone" `Quick test_gmst_monotone;
370 Alcotest.test_case "GMST ptime" `Quick test_gmst_ptime;
371 Alcotest.test_case "JD at J2000" `Quick test_jd_j2000;
372 Alcotest.test_case "JD at unix epoch" `Quick test_jd_unix_epoch;
373 Alcotest.test_case "JD round-trip" `Quick test_jd_roundtrip;
374 Alcotest.test_case "TEME→ECEF at GMST=0" `Quick test_teme_ecef_zero;
375 Alcotest.test_case "TEME→ECEF at 90°" `Quick test_teme_ecef_90deg;
376 Alcotest.test_case "TEME↔ECEF round-trip" `Quick test_teme_ecef_roundtrip;
377 Alcotest.test_case "rotation preserves length" `Quick
378 test_rotation_preserves_length;
379 Alcotest.test_case "geodetic equator PM" `Quick test_geodetic_equator_pm;
380 Alcotest.test_case "geodetic north pole" `Quick test_geodetic_north_pole;
381 Alcotest.test_case "geodetic lon 90°" `Quick test_geodetic_lon90;
382 Alcotest.test_case "geodetic Washington" `Quick test_geodetic_washington;
383 Alcotest.test_case "geodetic round-trip" `Quick test_geodetic_roundtrip;
384 Alcotest.test_case "geodetic lats" `Quick test_geodetic_roundtrip_lats;
385 Alcotest.test_case "geodetic ISS alt" `Quick test_geodetic_iss_altitude;
386 Alcotest.test_case "geodetic GEO alt" `Quick test_geodetic_geo_altitude;
387 Alcotest.test_case "Vallado 3-3" `Quick test_vallado_3_3;
388 Alcotest.test_case "TEME→geodetic equator" `Quick
389 test_teme_geodetic_equator;
390 Alcotest.test_case "TEME→geodetic GMST shift" `Quick
391 test_teme_geodetic_gmst_shift;
392 Alcotest.test_case "normalize angle" `Quick test_normalize_angle;
393 Alcotest.test_case "TEME→J2000 not identity" `Quick
394 test_teme_j2000_not_identity;
395 Alcotest.test_case "TEME↔J2000 round-trip" `Quick
396 test_teme_j2000_roundtrip;
397 Alcotest.test_case "TEME→J2000 reference" `Quick test_teme_j2000_reference;
398 Alcotest.test_case "TEME→J2000 preserves length" `Quick
399 test_teme_j2000_preserves_length;
400 Alcotest.test_case "TEME→J2000 at J2000 epoch" `Quick
401 test_teme_j2000_at_epoch;
402 Alcotest.test_case "TEME→J2000 precession grows" `Quick
403 test_teme_j2000_multiple_epochs;
404 Alcotest.test_case "J2000→TEME round-trip" `Quick
405 test_j2000_teme_reference;
406 ] )