gazagnaire.org / ocaml-coordinate
Astrodynamics coordinate frame transforms
0
fork

Configure Feed

Select the types of activity you want to include in your feed.

Add TEME/J2000 conversion tests for ocaml-coordinate

7 new tests validated against python reference computation:
- Not-identity (precession moves vector ~40 km at LEO over 26 years)
- Round-trip TEME->J2000->TEME (< 1e-9 km)
- Reference vector comparison (< 1 km vs python IAU-76/IAU-1980)
- Length preservation (rotation is isometry)
- At J2000 epoch (only nutation contributes, < 2 km)
- Precession grows with time (2050 > 2026)
- J2000->TEME->J2000 inverse round-trip

Thomas Gazagnaire 9f1e76f3 66fcce5e

+108
+108
test/test_coordinate.ml
··· 267 267 (* Suite *) 268 268 (* ================================================================== *) 269 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 + 278 + let unix_2026 = 1767225600.0 279 + 280 + (* TEME->J2000 should change the position vector (not identity). *) 281 + let 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. *) 290 + let 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 *) 303 + let 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 + Printf.printf " 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). *) 317 + let 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. *) 329 + let 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 + Printf.printf " 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. *) 341 + let 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 + Printf.printf " 2026: %.3f km, 2050: %.3f km\n" d1 d2 354 + 355 + (* Inverse: J2000 -> TEME reference. *) 356 + let 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 362 + Alcotest.failf "J2000->TEME->J2000 error: %.15e km" err 363 + 270 364 let suite = 271 365 ( "coordinate", 272 366 [ ··· 297 391 Alcotest.test_case "TEME→geodetic GMST shift" `Quick 298 392 test_teme_geodetic_gmst_shift; 299 393 Alcotest.test_case "normalize angle" `Quick test_normalize_angle; 394 + Alcotest.test_case "TEME→J2000 not identity" `Quick 395 + test_teme_j2000_not_identity; 396 + Alcotest.test_case "TEME↔J2000 round-trip" `Quick 397 + test_teme_j2000_roundtrip; 398 + Alcotest.test_case "TEME→J2000 reference" `Quick 399 + test_teme_j2000_reference; 400 + Alcotest.test_case "TEME→J2000 preserves length" `Quick 401 + test_teme_j2000_preserves_length; 402 + Alcotest.test_case "TEME→J2000 at J2000 epoch" `Quick 403 + test_teme_j2000_at_epoch; 404 + Alcotest.test_case "TEME→J2000 precession grows" `Quick 405 + test_teme_j2000_multiple_epochs; 406 + Alcotest.test_case "J2000→TEME round-trip" `Quick 407 + test_j2000_teme_reference; 300 408 ] )