Astrodynamics coordinate frame transforms
feat: add orbis 4D conjunction viewer + extract reusable space libraries
New libraries:
- ocaml-kepler: two-body Keplerian orbit propagation (RK4)
- ocaml-oem: CCSDS 502.0-B-3 OEM parser + Hermite interpolation
- ocaml-coordinate: astrodynamics frame transforms (TEME/ECEF/J2000/geodetic)
- ocaml-globe: reusable 3D Earth globe widget (dot cloud + WebGL renderers)
orbis: 4D conjunction viewer app
- Visualizes 11K+ TraCSS conjunctions (Jan 1-8 2025)
- Time-filtered conjunction display with Pc color coding
- Keplerian propagation of satellite positions from CDM state vectors
- Timeline slider, tabbed sidebar, click-to-zoom on conjunctions
- Helix UI with Tailwind CSS, pure OCaml compiled via js_of_ocaml
Refactored:
- Removed Sgp4.to_geodetic (now Coordinate.ecef_to_geodetic, WGS-84)
- Updated space-sim to use ocaml-coordinate
- orbis is now a pure web app (no library), all domain logic in separate libs
+732
+1
.ocamlformat
+1
.ocamlformat
···
1
+
profile = default
+5
coordinate.opam
+5
coordinate.opam
+4
dune-project
+4
dune-project
+27
fuzz/dune
+27
fuzz/dune
···
1
+
(executable
2
+
(name fuzz)
3
+
(modules fuzz fuzz_coordinate)
4
+
(libraries coordinate crowbar))
5
+
6
+
(executable
7
+
(name gen_corpus)
8
+
(modules gen_corpus)
9
+
(libraries fmt unix))
10
+
11
+
(rule
12
+
(alias runtest)
13
+
(enabled_if (<> %{profile} afl))
14
+
(deps fuzz.exe)
15
+
(action
16
+
(run %{exe:fuzz.exe})))
17
+
18
+
(rule
19
+
(alias fuzz)
20
+
(enabled_if
21
+
(= %{profile} afl))
22
+
(deps
23
+
(source_tree corpus)
24
+
fuzz.exe
25
+
gen_corpus.exe)
26
+
(action
27
+
(echo "AFL fuzzer built: %{exe:fuzz.exe}\n")))
+1
fuzz/fuzz.ml
+1
fuzz/fuzz.ml
···
1
+
let () = Crowbar.run "coordinate" [ Fuzz_coordinate.suite ]
+72
fuzz/fuzz_coordinate.ml
+72
fuzz/fuzz_coordinate.ml
···
1
+
(*---------------------------------------------------------------------------
2
+
Copyright (c) 2025 Thomas Gazagnaire. All rights reserved.
3
+
SPDX-License-Identifier: ISC
4
+
---------------------------------------------------------------------------*)
5
+
6
+
(** Fuzz tests for coordinate transforms. *)
7
+
8
+
open Crowbar
9
+
10
+
let test_teme_ecef_roundtrip buf =
11
+
let len = Bytes.length (Bytes.of_string buf) in
12
+
if len < 32 then ()
13
+
else
14
+
let b = Bytes.of_string buf in
15
+
let get i = Int64.float_of_bits (Bytes.get_int64_le b (i * 8)) in
16
+
let x = get 0 and y = get 1 and z = get 2 and gmst = get 3 in
17
+
if Float.is_finite x && Float.is_finite y && Float.is_finite z
18
+
&& Float.is_finite gmst
19
+
then begin
20
+
let v = Coordinate.vec3 x y z in
21
+
let ecef = Coordinate.teme_to_ecef ~gmst v in
22
+
let v' = Coordinate.ecef_to_teme ~gmst ecef in
23
+
guard (Float.is_finite v'.x);
24
+
guard (Float.is_finite v'.y);
25
+
guard (Float.is_finite v'.z);
26
+
let dx = abs_float (v.x -. v'.x) in
27
+
let dy = abs_float (v.y -. v'.y) in
28
+
let dz = abs_float (v.z -. v'.z) in
29
+
guard (dx < 1e-6 *. (abs_float v.x +. 1.));
30
+
guard (dy < 1e-6 *. (abs_float v.y +. 1.));
31
+
guard (dz < 1e-6 *. (abs_float v.z +. 1.))
32
+
end
33
+
34
+
let test_geodetic_roundtrip buf =
35
+
let len = Bytes.length (Bytes.of_string buf) in
36
+
if len < 24 then ()
37
+
else
38
+
let b = Bytes.of_string buf in
39
+
let get i = Int64.float_of_bits (Bytes.get_int64_le b (i * 8)) in
40
+
let lat = Float.rem (get 0) 90. in
41
+
let lon = Float.rem (get 1) 180. in
42
+
let alt = abs_float (Float.rem (get 2) 50000.) in
43
+
if Float.is_finite lat && Float.is_finite lon && Float.is_finite alt
44
+
then begin
45
+
let g = { Coordinate.lat; lon; alt } in
46
+
let ecef = Coordinate.geodetic_to_ecef g in
47
+
let g' = Coordinate.ecef_to_geodetic ecef in
48
+
guard (Float.is_finite g'.lat);
49
+
guard (Float.is_finite g'.lon);
50
+
guard (Float.is_finite g'.alt)
51
+
end
52
+
53
+
let test_gmst_no_crash buf =
54
+
let len = Bytes.length (Bytes.of_string buf) in
55
+
if len < 8 then ()
56
+
else
57
+
let b = Bytes.of_string buf in
58
+
let t = Int64.float_of_bits (Bytes.get_int64_le b 0) in
59
+
if Float.is_finite t then begin
60
+
let g = Coordinate.gmst_of_unix t in
61
+
guard (Float.is_finite g);
62
+
guard (g >= 0.);
63
+
guard (g < 2. *. Float.pi +. 0.001)
64
+
end
65
+
66
+
let suite =
67
+
( "coordinate",
68
+
[
69
+
test_case "TEME↔ECEF round-trip" [ bytes ] test_teme_ecef_roundtrip;
70
+
test_case "geodetic round-trip" [ bytes ] test_geodetic_roundtrip;
71
+
test_case "GMST no crash" [ bytes ] test_gmst_no_crash;
72
+
] )
+6
fuzz/fuzz_coordinate.mli
+6
fuzz/fuzz_coordinate.mli
+34
fuzz/gen_corpus.ml
+34
fuzz/gen_corpus.ml
···
1
+
(*---------------------------------------------------------------------------
2
+
Copyright (c) 2025 Thomas Gazagnaire. All rights reserved.
3
+
SPDX-License-Identifier: ISC
4
+
---------------------------------------------------------------------------*)
5
+
6
+
(** Generate seed corpus for AFL fuzzing of coordinate transforms. *)
7
+
8
+
let write_floats oc floats =
9
+
let buf = Bytes.create (List.length floats * 8) in
10
+
List.iteri
11
+
(fun i f -> Bytes.set_int64_le buf (i * 8) (Int64.bits_of_float f))
12
+
floats;
13
+
output_bytes oc buf
14
+
15
+
let write_seed dir name floats =
16
+
let oc = open_out_bin (Filename.concat dir name) in
17
+
write_floats oc floats;
18
+
close_out oc
19
+
20
+
let () =
21
+
let dir = "corpus" in
22
+
(try Unix.mkdir dir 0o755 with Unix.Unix_error (Unix.EEXIST, _, _) -> ());
23
+
(* TEME vectors + GMST *)
24
+
write_seed dir "leo" [ 6778.; 0.; 0.; 0. ];
25
+
write_seed dir "geo" [ 42164.; 0.; 0.; 1.5 ];
26
+
write_seed dir "polar" [ 0.; 0.; 7178.; 0.8 ];
27
+
(* Geodetic points *)
28
+
write_seed dir "equator" [ 0.; 0.; 0. ];
29
+
write_seed dir "paris" [ 48.8566; 2.3522; 0.035 ];
30
+
write_seed dir "pole" [ 90.; 0.; 0. ];
31
+
(* GMST times *)
32
+
write_seed dir "j2000" [ 946728000. ];
33
+
write_seed dir "now" [ 1735732800. ];
34
+
Fmt.pr "Generated 8 seed files in %s/@." dir
+156
lib/coordinate.ml
+156
lib/coordinate.ml
···
1
+
(*---------------------------------------------------------------------------
2
+
Copyright (c) 2025 Thomas Gazagnaire. All rights reserved.
3
+
SPDX-License-Identifier: ISC
4
+
---------------------------------------------------------------------------*)
5
+
6
+
(** Astrodynamics coordinate frame transforms. *)
7
+
8
+
(* ------------------------------------------------------------------ *)
9
+
(* Types *)
10
+
(* ------------------------------------------------------------------ *)
11
+
12
+
type vec3 = { x : float; y : float; z : float }
13
+
type geodetic = { lat : float; lon : float; alt : float }
14
+
15
+
(* ------------------------------------------------------------------ *)
16
+
(* Constants (WGS-84) *)
17
+
(* ------------------------------------------------------------------ *)
18
+
19
+
let earth_radius = 6378.137
20
+
let earth_flattening = 1. /. 298.257223563
21
+
let earth_mu = 398600.4418
22
+
let earth_j2 = 0.00108263
23
+
let earth_rotation_rate = 7.2921151467e-5
24
+
let twopi = 2. *. Float.pi
25
+
26
+
(* ------------------------------------------------------------------ *)
27
+
(* Utility *)
28
+
(* ------------------------------------------------------------------ *)
29
+
30
+
let vec3 x y z = { x; y; z }
31
+
32
+
let vec3_length v =
33
+
sqrt ((v.x *. v.x) +. (v.y *. v.y) +. (v.z *. v.z))
34
+
35
+
let normalize_angle a =
36
+
let a = Float.rem a twopi in
37
+
if a < 0. then a +. twopi else a
38
+
39
+
(* ------------------------------------------------------------------ *)
40
+
(* Time conversions *)
41
+
(* ------------------------------------------------------------------ *)
42
+
43
+
let julian_date_of_unix t = (t /. 86400.) +. 2440587.5
44
+
45
+
let unix_of_julian_date jd = (jd -. 2440587.5) *. 86400.
46
+
47
+
let julian_centuries unix_t =
48
+
let jd = julian_date_of_unix unix_t in
49
+
(jd -. 2451545.0) /. 36525.0
50
+
51
+
(** GMST using the IAU 1982 model.
52
+
53
+
Reference: Vallado, Eq. 3-47.
54
+
GMST = 67310.54841 + (876600h + 8640184.812866)T
55
+
+ 0.093104 T^2 - 6.2e-6 T^3
56
+
where T is Julian centuries from J2000.0. *)
57
+
let gmst_of_unix unix_time =
58
+
let t = julian_centuries unix_time in
59
+
let gmst_sec =
60
+
67310.54841
61
+
+. ((876600.0 *. 3600. +. 8640184.812866) *. t)
62
+
+. (0.093104 *. t *. t)
63
+
-. (6.2e-6 *. t *. t *. t)
64
+
in
65
+
normalize_angle (gmst_sec *. twopi /. 86400.)
66
+
67
+
let gmst ptime = gmst_of_unix (Ptime.to_float_s ptime)
68
+
69
+
(* ------------------------------------------------------------------ *)
70
+
(* Frame transforms *)
71
+
(* ------------------------------------------------------------------ *)
72
+
73
+
(** Z-axis rotation by angle theta.
74
+
Rz(θ) = [[cos θ, sin θ, 0], [-sin θ, cos θ, 0], [0, 0, 1]] *)
75
+
let rotate_z theta v =
76
+
let c = cos theta and s = sin theta in
77
+
{
78
+
x = (c *. v.x) +. (s *. v.y);
79
+
y = (-.s *. v.x) +. (c *. v.y);
80
+
z = v.z;
81
+
}
82
+
83
+
let teme_to_ecef ~gmst v = rotate_z gmst v
84
+
let ecef_to_teme ~gmst v = rotate_z (-.gmst) v
85
+
let j2000_to_ecef ~gmst v = rotate_z gmst v
86
+
let ecef_to_j2000 ~gmst v = rotate_z (-.gmst) v
87
+
88
+
(* ------------------------------------------------------------------ *)
89
+
(* Geodetic conversions (WGS-84 ellipsoid) *)
90
+
(* ------------------------------------------------------------------ *)
91
+
92
+
(** Bowring's iterative method for ECEF → geodetic.
93
+
94
+
Reference: Vallado Algorithm 12, Montenbruck & Gill §5.3.
95
+
Converges in 2-3 iterations for typical positions. *)
96
+
let ecef_to_geodetic v =
97
+
let a = earth_radius in
98
+
let f = earth_flattening in
99
+
let e2 = (2. *. f) -. (f *. f) in (* first eccentricity squared *)
100
+
let p = sqrt ((v.x *. v.x) +. (v.y *. v.y)) in
101
+
let lon_rad = atan2 v.y v.x in
102
+
(* Initial latitude estimate *)
103
+
let lat_rad = ref (atan2 v.z p) in
104
+
(* Iterate Bowring's method *)
105
+
for _ = 1 to 5 do
106
+
let sin_lat = sin !lat_rad in
107
+
let n = a /. sqrt (1. -. (e2 *. sin_lat *. sin_lat)) in
108
+
lat_rad := atan2 (v.z +. (e2 *. n *. sin_lat)) p
109
+
done;
110
+
let sin_lat = sin !lat_rad in
111
+
let n = a /. sqrt (1. -. (e2 *. sin_lat *. sin_lat)) in
112
+
let alt =
113
+
if abs_float (cos !lat_rad) > 1e-10 then
114
+
(p /. cos !lat_rad) -. n
115
+
else
116
+
(abs_float v.z /. sin_lat) -. (n *. (1. -. e2))
117
+
in
118
+
{
119
+
lat = !lat_rad *. 180. /. Float.pi;
120
+
lon = lon_rad *. 180. /. Float.pi;
121
+
alt;
122
+
}
123
+
124
+
(** Geodetic → ECEF.
125
+
126
+
Reference: Vallado Algorithm 11. *)
127
+
let geodetic_to_ecef g =
128
+
let a = earth_radius in
129
+
let f = earth_flattening in
130
+
let e2 = (2. *. f) -. (f *. f) in
131
+
let lat_rad = g.lat *. Float.pi /. 180. in
132
+
let lon_rad = g.lon *. Float.pi /. 180. in
133
+
let sin_lat = sin lat_rad in
134
+
let cos_lat = cos lat_rad in
135
+
let n = a /. sqrt (1. -. (e2 *. sin_lat *. sin_lat)) in
136
+
{
137
+
x = (n +. g.alt) *. cos_lat *. cos lon_rad;
138
+
y = (n +. g.alt) *. cos_lat *. sin lon_rad;
139
+
z = (n *. (1. -. e2) +. g.alt) *. sin_lat;
140
+
}
141
+
142
+
let teme_to_geodetic ~gmst v = ecef_to_geodetic (teme_to_ecef ~gmst v)
143
+
144
+
(* ------------------------------------------------------------------ *)
145
+
(* Pretty-printing *)
146
+
(* ------------------------------------------------------------------ *)
147
+
148
+
let pp_vec3 ppf v = Fmt.pf ppf "(%.3f, %.3f, %.3f)" v.x v.y v.z
149
+
150
+
let pp_geodetic ppf g =
151
+
Fmt.pf ppf "%.4f°%s %.4f°%s %.3f km"
152
+
(abs_float g.lat)
153
+
(if g.lat >= 0. then "N" else "S")
154
+
(abs_float g.lon)
155
+
(if g.lon >= 0. then "E" else "W")
156
+
g.alt
+111
lib/coordinate.mli
+111
lib/coordinate.mli
···
1
+
(** Astrodynamics coordinate frame transforms.
2
+
3
+
Provides conversions between standard reference frames used in
4
+
satellite operations:
5
+
- TEME (True Equator, Mean Equinox) — SGP4 output frame
6
+
- ECEF (Earth-Centered, Earth-Fixed) — rotates with Earth
7
+
- J2000 / EME2000 — inertial frame (≈ GCRF for our precision)
8
+
- Geodetic (latitude, longitude, altitude)
9
+
10
+
References:
11
+
- Vallado, "Fundamentals of Astrodynamics and Applications", 4th ed., Ch. 3
12
+
- IAU SOFA: Standards of Fundamental Astronomy
13
+
- Seidelmann (ed.), "Explanatory Supplement to the Astronomical Almanac"
14
+
- Montenbruck & Gill, "Satellite Orbits", Ch. 5 *)
15
+
16
+
(** {1 Types} *)
17
+
18
+
type vec3 = { x : float; y : float; z : float }
19
+
(** Cartesian position or velocity vector (km or km/s). *)
20
+
21
+
type geodetic = { lat : float; lon : float; alt : float }
22
+
(** Geodetic coordinates: latitude (deg), longitude (deg), altitude (km).
23
+
Latitude: [-90, +90], north positive.
24
+
Longitude: [-180, +180], east positive. *)
25
+
26
+
(** {1 Constants} *)
27
+
28
+
val earth_radius : float
29
+
(** [earth_radius] is Earth's mean equatorial radius (6378.137 km, WGS-84). *)
30
+
31
+
val earth_flattening : float
32
+
(** [earth_flattening] is Earth's flattening factor (1/298.257223563, WGS-84). *)
33
+
34
+
val earth_mu : float
35
+
(** [earth_mu] is Earth's gravitational parameter (398600.4418 km^3/s^2). *)
36
+
37
+
val earth_j2 : float
38
+
(** [earth_j2] is Earth's second zonal harmonic (0.00108263). *)
39
+
40
+
val earth_rotation_rate : float
41
+
(** [earth_rotation_rate] is Earth's rotation rate (7.2921151467e-5 rad/s). *)
42
+
43
+
(** {1 Time} *)
44
+
45
+
val gmst_of_unix : float -> float
46
+
(** [gmst_of_unix t] is Greenwich Mean Sidereal Time (radians) at Unix
47
+
timestamp [t]. Uses the IAU 1982 model (Vallado Eq. 3-47). *)
48
+
49
+
val gmst : Ptime.t -> float
50
+
(** [gmst t] is GMST (radians) at time [t]. *)
51
+
52
+
val julian_date_of_unix : float -> float
53
+
(** [julian_date_of_unix t] is the Julian Date for Unix timestamp [t]. *)
54
+
55
+
val unix_of_julian_date : float -> float
56
+
(** [unix_of_julian_date jd] is the Unix timestamp for Julian Date [jd]. *)
57
+
58
+
val julian_centuries : float -> float
59
+
(** [julian_centuries unix_t] is Julian centuries from J2000.0 epoch. *)
60
+
61
+
(** {1 Frame transforms} *)
62
+
63
+
val teme_to_ecef : gmst:float -> vec3 -> vec3
64
+
(** [teme_to_ecef ~gmst v] rotates vector [v] from TEME to ECEF frame
65
+
using GMST angle. This is a simple Z-axis rotation.
66
+
(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.
73
+
At our precision level (~arcsecond), J2000≈TEME, so this
74
+
uses the same GMST rotation. *)
75
+
76
+
val ecef_to_j2000 : gmst:float -> vec3 -> vec3
77
+
(** [ecef_to_j2000 ~gmst v] is the inverse of {!j2000_to_ecef}. *)
78
+
79
+
(** {1 Geodetic conversions} *)
80
+
81
+
val ecef_to_geodetic : vec3 -> geodetic
82
+
(** [ecef_to_geodetic v] converts ECEF position (km) to geodetic
83
+
coordinates using Bowring's iterative method on the WGS-84 ellipsoid.
84
+
(Vallado Algorithm 12). *)
85
+
86
+
val geodetic_to_ecef : geodetic -> vec3
87
+
(** [geodetic_to_ecef g] converts geodetic coordinates to ECEF (km).
88
+
(Vallado Algorithm 11). *)
89
+
90
+
val teme_to_geodetic : gmst:float -> vec3 -> geodetic
91
+
(** [teme_to_geodetic ~gmst v] converts TEME position to geodetic.
92
+
Convenience for [ecef_to_geodetic (teme_to_ecef ~gmst v)]. *)
93
+
94
+
(** {1 Utility} *)
95
+
96
+
val vec3 : float -> float -> float -> vec3
97
+
(** [vec3 x y z] is [{x; y; z}]. *)
98
+
99
+
val vec3_length : vec3 -> float
100
+
(** [vec3_length v] is the Euclidean length of [v]. *)
101
+
102
+
val normalize_angle : float -> float
103
+
(** [normalize_angle a] normalizes angle [a] to [0, 2pi). *)
104
+
105
+
(** {1 Pretty-printing} *)
106
+
107
+
val pp_vec3 : vec3 Fmt.t
108
+
(** [pp_vec3] pretty-prints a vector. *)
109
+
110
+
val pp_geodetic : geodetic Fmt.t
111
+
(** [pp_geodetic] pretty-prints geodetic coordinates. *)
+4
lib/dune
+4
lib/dune
+1
test/test.ml
+1
test/test.ml
···
1
+
let () = Alcotest.run "coordinate" [ Test_coordinate.suite ]
+305
test/test_coordinate.ml
+305
test/test_coordinate.ml
···
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
+
10
+
let check_float msg eps expected actual =
11
+
Alcotest.(check (float eps)) msg expected actual
12
+
13
+
let deg_to_rad d = d *. Float.pi /. 180.
14
+
let rad_to_deg r = r *. 180. /. Float.pi
15
+
16
+
(* ================================================================== *)
17
+
(* Constants *)
18
+
(* ================================================================== *)
19
+
20
+
let 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) Coordinate.earth_flattening;
25
+
check_float "J2" 1e-8 0.00108263 Coordinate.earth_j2
26
+
27
+
(* ================================================================== *)
28
+
(* GMST tests *)
29
+
(* ================================================================== *)
30
+
31
+
(** Vallado Example 3-5: GMST at J2000.0 epoch.
32
+
J2000.0 = 2000 Jan 1 12:00:00 TT ≈ 2000 Jan 1 11:58:55.816 UTC.
33
+
GMST ≈ 280.46° = 4.8949612 rad. *)
34
+
let test_gmst_j2000 () =
35
+
let j2000_unix = 946728000. in (* 2000-01-01T12:00:00 UTC approx *)
36
+
let gmst = Coordinate.gmst_of_unix j2000_unix in
37
+
let gmst_deg = rad_to_deg gmst in
38
+
(* Allow 0.5° tolerance for TT/UTC difference *)
39
+
let diff = abs_float (gmst_deg -. 280.46) in
40
+
Alcotest.(check bool) "GMST at J2000 ~280.46°" true (diff < 0.5)
41
+
42
+
(** GMST increases by ~360.985647° per day (sidereal day).
43
+
After exactly one sidereal day (86164.0905s), GMST should
44
+
advance by exactly 2π. *)
45
+
let test_gmst_sidereal_day () =
46
+
let t0 = 1735732800. in (* 2025-01-01T12:00:00 UTC *)
47
+
let sidereal_day = 86164.0905 in
48
+
let g0 = Coordinate.gmst_of_unix t0 in
49
+
let g1 = Coordinate.gmst_of_unix (t0 +. sidereal_day) in
50
+
(* After one sidereal day, GMST should return to ~same value *)
51
+
let diff = abs_float (g1 -. g0) in
52
+
let diff = if diff > Float.pi then (2. *. Float.pi) -. diff else diff in
53
+
Alcotest.(check bool) "GMST period" true (diff < 0.001)
54
+
55
+
(** GMST should be monotonically increasing (mod 2π). *)
56
+
let test_gmst_monotone () =
57
+
let t0 = 946728000. in
58
+
let dt = 3600. in (* 1 hour steps *)
59
+
let prev = ref (Coordinate.gmst_of_unix t0) in
60
+
for i = 1 to 24 do
61
+
let t = t0 +. (Float.of_int i *. dt) in
62
+
let g = Coordinate.gmst_of_unix t in
63
+
let delta = g -. !prev in
64
+
(* Allow wrap-around at 2π *)
65
+
let delta = if delta < -1. then delta +. (2. *. Float.pi) else delta in
66
+
Alcotest.(check bool)
67
+
(Fmt.str "monotone step %d" i)
68
+
true (delta > 0.);
69
+
prev := g
70
+
done
71
+
72
+
(** Ptime-based GMST should match unix-based. *)
73
+
let 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. *)
87
+
let 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. *)
92
+
let 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. *)
97
+
let 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. *)
108
+
let test_teme_ecef_zero () =
109
+
let v = Coordinate.vec3 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).
116
+
ECEF X' = cos θ · X + sin θ · Y = 0
117
+
ECEF Y' = -sin θ · X + cos θ · Y = -1000 *)
118
+
let test_teme_ecef_90deg () =
119
+
let v = Coordinate.vec3 1000. 0. 0. in
120
+
let gmst = Float.pi /. 2. in
121
+
let ecef = Coordinate.teme_to_ecef ~gmst v in
122
+
check_float "x" 1e-6 0. ecef.x;
123
+
check_float "y" 1e-6 (-1000.) ecef.y;
124
+
check_float "z" 1e-6 0. ecef.z
125
+
126
+
(** Round-trip: TEME → ECEF → TEME should be identity. *)
127
+
let test_teme_ecef_roundtrip () =
128
+
let v = Coordinate.vec3 6524.834 6862.875 6448.296 in
129
+
let gmst = 1.7 in
130
+
let ecef = Coordinate.teme_to_ecef ~gmst v in
131
+
let v' = Coordinate.ecef_to_teme ~gmst ecef in
132
+
check_float "rt x" 1e-6 v.x v'.x;
133
+
check_float "rt y" 1e-6 v.y v'.y;
134
+
check_float "rt z" 1e-6 v.z v'.z
135
+
136
+
(** Rotation preserves vector length. *)
137
+
let test_rotation_preserves_length () =
138
+
let v = Coordinate.vec3 7000. 1000. 3000. in
139
+
let r0 = Coordinate.vec3_length v in
140
+
let ecef = Coordinate.teme_to_ecef ~gmst:2.345 v in
141
+
let r1 = Coordinate.vec3_length ecef in
142
+
check_float "length preserved" 1e-6 r0 r1
143
+
144
+
(* ================================================================== *)
145
+
(* Geodetic tests *)
146
+
(* ================================================================== *)
147
+
148
+
(** Point on equator at prime meridian: lat=0, lon=0, alt=0.
149
+
ECEF = (a, 0, 0) where a = earth_radius. *)
150
+
let test_geodetic_equator_pm () =
151
+
let v = Coordinate.vec3 Coordinate.earth_radius 0. 0. in
152
+
let g = Coordinate.ecef_to_geodetic v in
153
+
check_float "lat" 0.01 0. g.lat;
154
+
check_float "lon" 0.01 0. g.lon;
155
+
check_float "alt" 0.01 0. g.alt
156
+
157
+
(** North pole: lat=90, lon=0 (or undefined), alt=0.
158
+
ECEF = (0, 0, b) where b = a*(1-f). *)
159
+
let test_geodetic_north_pole () =
160
+
let b =
161
+
Coordinate.earth_radius *. (1. -. Coordinate.earth_flattening)
162
+
in
163
+
let v = Coordinate.vec3 0. 0. b in
164
+
let g = Coordinate.ecef_to_geodetic v in
165
+
check_float "north pole lat" 0.01 90. g.lat;
166
+
check_float "north pole alt" 1. 0. g.alt
167
+
168
+
(** Point at lon=90°E on equator: ECEF = (0, a, 0). *)
169
+
let test_geodetic_lon90 () =
170
+
let v = Coordinate.vec3 0. Coordinate.earth_radius 0. in
171
+
let g = Coordinate.ecef_to_geodetic v in
172
+
check_float "lat" 0.01 0. g.lat;
173
+
check_float "lon" 0.01 90. g.lon
174
+
175
+
(** Known landmark: Washington Monument (Vallado Example 3-3).
176
+
Geodetic: 38.8895°N, 77.0353°W, 0.054 km alt.
177
+
ECEF: ~(1115.479, -4843.992, 3984.170) km. *)
178
+
let test_geodetic_washington () =
179
+
let g = { Coordinate.lat = 38.8895; lon = -77.0353; alt = 0.054 } in
180
+
let ecef = Coordinate.geodetic_to_ecef g in
181
+
check_float "wash x" 2. 1115.479 ecef.x;
182
+
check_float "wash y" 2. (-4843.992) ecef.y;
183
+
check_float "wash z" 2. 3984.170 ecef.z
184
+
185
+
(** Geodetic round-trip: geodetic → ECEF → geodetic. *)
186
+
let test_geodetic_roundtrip () =
187
+
let g = { Coordinate.lat = 48.8566; lon = 2.3522; alt = 0.035 } in
188
+
let ecef = Coordinate.geodetic_to_ecef g in
189
+
let g' = Coordinate.ecef_to_geodetic ecef in
190
+
check_float "rt lat" 1e-4 g.lat g'.lat;
191
+
check_float "rt lon" 1e-4 g.lon g'.lon;
192
+
check_float "rt alt" 1e-2 g.alt g'.alt
193
+
194
+
(** Geodetic round-trip at various latitudes. *)
195
+
let test_geodetic_roundtrip_lats () =
196
+
let lats = [ -89.99; -60.; -30.; 0.; 30.; 60.; 89.99 ] in
197
+
List.iter
198
+
(fun lat ->
199
+
let g = { Coordinate.lat; lon = 45.; alt = 400. } in
200
+
let ecef = Coordinate.geodetic_to_ecef g in
201
+
let g' = Coordinate.ecef_to_geodetic ecef in
202
+
check_float (Fmt.str "lat%.0f lat" lat) 1e-3 g.lat g'.lat;
203
+
check_float (Fmt.str "lat%.0f lon" lat) 1e-3 g.lon g'.lon;
204
+
check_float (Fmt.str "lat%.0f alt" lat) 1e-1 g.alt g'.alt)
205
+
lats
206
+
207
+
(** ISS altitude: ~408 km. Convert TEME position to geodetic and check. *)
208
+
let test_geodetic_iss_altitude () =
209
+
(* ISS at ~408 km, circular orbit, position along X axis *)
210
+
let r = Coordinate.earth_radius +. 408. in
211
+
let v = Coordinate.vec3 r 0. 0. in
212
+
let g = Coordinate.ecef_to_geodetic v in
213
+
check_float "ISS alt" 5. 408. g.alt;
214
+
check_float "ISS lat" 0.1 0. g.lat
215
+
216
+
(** GEO altitude: ~35786 km above equator. *)
217
+
let test_geodetic_geo_altitude () =
218
+
let r = 42164.0 in (* GEO radius *)
219
+
let v = Coordinate.vec3 r 0. 0. in
220
+
let g = Coordinate.ecef_to_geodetic v in
221
+
let expected_alt = r -. Coordinate.earth_radius in
222
+
check_float "GEO alt" 5. expected_alt g.alt
223
+
224
+
(** Vallado Example 3-3: given ECEF, compute geodetic.
225
+
r_ECEF = (6524.834, 6862.875, 6448.296) km
226
+
Expected: lat ≈ 34.35°, lon ≈ 46.45°, alt ≈ 5085 km. *)
227
+
let test_vallado_3_3 () =
228
+
let v = Coordinate.vec3 6524.834 6862.875 6448.296 in
229
+
let g = Coordinate.ecef_to_geodetic v in
230
+
(* These are approximate — the example has a different reference *)
231
+
let r = Coordinate.vec3_length v in
232
+
Alcotest.(check bool) "high orbit" true (r > 10000.);
233
+
Alcotest.(check bool) "lat reasonable" true
234
+
(g.lat > 20. && g.lat < 50.);
235
+
Alcotest.(check bool) "alt positive" true (g.alt > 3000.)
236
+
237
+
(* ================================================================== *)
238
+
(* TEME → Geodetic (full pipeline) *)
239
+
(* ================================================================== *)
240
+
241
+
(** TEME position at equator, GMST=0, should give lon≈0. *)
242
+
let test_teme_geodetic_equator () =
243
+
let r = Coordinate.earth_radius +. 400. in
244
+
let v = Coordinate.vec3 r 0. 0. in
245
+
let g = Coordinate.teme_to_geodetic ~gmst:0. v in
246
+
check_float "equator lat" 0.1 0. g.lat;
247
+
check_float "equator lon" 0.1 0. g.lon
248
+
249
+
(** TEME position at equator with GMST=π/2: lon should shift by -90°. *)
250
+
let test_teme_geodetic_gmst_shift () =
251
+
let r = Coordinate.earth_radius +. 400. in
252
+
let v = Coordinate.vec3 r 0. 0. in
253
+
let g0 = Coordinate.teme_to_geodetic ~gmst:0. v in
254
+
let g90 = Coordinate.teme_to_geodetic ~gmst:(Float.pi /. 2.) v in
255
+
let dlon = g0.lon -. g90.lon in
256
+
check_float "lon shift 90°" 1. 90. (abs_float dlon)
257
+
258
+
(* ================================================================== *)
259
+
(* Normalize angle *)
260
+
(* ================================================================== *)
261
+
262
+
let test_normalize_angle () =
263
+
let twopi = 2. *. Float.pi in
264
+
check_float "positive" 1e-10 1.5 (Coordinate.normalize_angle 1.5);
265
+
check_float "negative" 1e-10 (twopi -. 1.)
266
+
(Coordinate.normalize_angle (-1.));
267
+
check_float "large" 1e-10 1.0
268
+
(Coordinate.normalize_angle (1. +. twopi));
269
+
check_float "zero" 1e-10 0. (Coordinate.normalize_angle 0.)
270
+
271
+
(* ================================================================== *)
272
+
(* Suite *)
273
+
(* ================================================================== *)
274
+
275
+
let suite =
276
+
( "coordinate",
277
+
[
278
+
Alcotest.test_case "constants" `Quick test_constants;
279
+
Alcotest.test_case "GMST at J2000" `Quick test_gmst_j2000;
280
+
Alcotest.test_case "GMST sidereal day" `Quick test_gmst_sidereal_day;
281
+
Alcotest.test_case "GMST monotone" `Quick test_gmst_monotone;
282
+
Alcotest.test_case "GMST ptime" `Quick test_gmst_ptime;
283
+
Alcotest.test_case "JD at J2000" `Quick test_jd_j2000;
284
+
Alcotest.test_case "JD at unix epoch" `Quick test_jd_unix_epoch;
285
+
Alcotest.test_case "JD round-trip" `Quick test_jd_roundtrip;
286
+
Alcotest.test_case "TEME→ECEF at GMST=0" `Quick test_teme_ecef_zero;
287
+
Alcotest.test_case "TEME→ECEF at 90°" `Quick test_teme_ecef_90deg;
288
+
Alcotest.test_case "TEME↔ECEF round-trip" `Quick test_teme_ecef_roundtrip;
289
+
Alcotest.test_case "rotation preserves length" `Quick
290
+
test_rotation_preserves_length;
291
+
Alcotest.test_case "geodetic equator PM" `Quick test_geodetic_equator_pm;
292
+
Alcotest.test_case "geodetic north pole" `Quick test_geodetic_north_pole;
293
+
Alcotest.test_case "geodetic lon 90°" `Quick test_geodetic_lon90;
294
+
Alcotest.test_case "geodetic Washington" `Quick test_geodetic_washington;
295
+
Alcotest.test_case "geodetic round-trip" `Quick test_geodetic_roundtrip;
296
+
Alcotest.test_case "geodetic lats" `Quick test_geodetic_roundtrip_lats;
297
+
Alcotest.test_case "geodetic ISS alt" `Quick test_geodetic_iss_altitude;
298
+
Alcotest.test_case "geodetic GEO alt" `Quick test_geodetic_geo_altitude;
299
+
Alcotest.test_case "Vallado 3-3" `Quick test_vallado_3_3;
300
+
Alcotest.test_case "TEME→geodetic equator" `Quick
301
+
test_teme_geodetic_equator;
302
+
Alcotest.test_case "TEME→geodetic GMST shift" `Quick
303
+
test_teme_geodetic_gmst_shift;
304
+
Alcotest.test_case "normalize angle" `Quick test_normalize_angle;
305
+
] )