Astrodynamics coordinate frame transforms
fix(ocaml-requests): update tests and fuzz for cstruct→Bytes migration
Test files still referenced Cstruct.t where the API now uses bytes.
Fixed all H2 frame, HPACK, client, and connection tests.
Fixed fuzz test. 330 tests pass.
+112
-94
+27
-2
coordinate.opam
+27
-2
coordinate.opam
···
1
1
# This file is generated by dune, edit dune-project instead
2
2
opam-version: "2.0"
3
-
name: "coordinate"
4
3
synopsis: "Astrodynamics coordinate frame transforms"
5
-
depends: []
4
+
maintainer: ["Thomas Gazagnaire <thomas@gazagnaire.org>"]
5
+
authors: ["Thomas Gazagnaire <thomas@gazagnaire.org>"]
6
+
license: "ISC"
7
+
homepage: "https://tangled.org/gazagnaire.org/ocaml-coordinate"
8
+
bug-reports: "https://tangled.org/gazagnaire.org/ocaml-coordinate/issues"
9
+
depends: [
10
+
"dune" {>= "3.21"}
11
+
"fmt"
12
+
"ptime"
13
+
"odoc" {with-doc}
14
+
]
15
+
build: [
16
+
["dune" "subst"] {dev}
17
+
[
18
+
"dune"
19
+
"build"
20
+
"-p"
21
+
name
22
+
"-j"
23
+
jobs
24
+
"@install"
25
+
"@runtest" {with-test}
26
+
"@doc" {with-doc}
27
+
]
28
+
]
29
+
dev-repo: "git+https://tangled.org/gazagnaire.org/ocaml-coordinate"
30
+
x-maintenance-intent: ["(latest)"]
+13
dune-project
+13
dune-project
···
2
2
(name coordinate)
3
3
(source (tangled gazagnaire.org/ocaml-coordinate))
4
4
(formatting (enabled_for ocaml))
5
+
6
+
(generate_opam_files true)
7
+
8
+
(license ISC)
9
+
(authors "Thomas Gazagnaire <thomas@gazagnaire.org>")
10
+
(maintainers "Thomas Gazagnaire <thomas@gazagnaire.org>")
11
+
12
+
(package
13
+
(name coordinate)
14
+
(synopsis "Astrodynamics coordinate frame transforms")
15
+
(depends
16
+
fmt
17
+
ptime))
+5
-5
fuzz/fuzz_coordinate.ml
+5
-5
fuzz/fuzz_coordinate.ml
···
14
14
let b = Bytes.of_string buf in
15
15
let get i = Int64.float_of_bits (Bytes.get_int64_le b (i * 8)) in
16
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
17
+
if
18
+
Float.is_finite x && Float.is_finite y && Float.is_finite z
19
+
&& Float.is_finite gmst
19
20
then begin
20
21
let v = Coordinate.vec3 x y z in
21
22
let ecef = Coordinate.teme_to_ecef ~gmst v in
···
40
41
let lat = Float.rem (get 0) 90. in
41
42
let lon = Float.rem (get 1) 180. in
42
43
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
44
+
if Float.is_finite lat && Float.is_finite lon && Float.is_finite alt then begin
45
45
let g = { Coordinate.lat; lon; alt } in
46
46
let ecef = Coordinate.geodetic_to_ecef g in
47
47
let g' = Coordinate.ecef_to_geodetic ecef in
···
60
60
let g = Coordinate.gmst_of_unix t in
61
61
guard (Float.is_finite g);
62
62
guard (g >= 0.);
63
-
guard (g < 2. *. Float.pi +. 0.001)
63
+
guard (g < (2. *. Float.pi) +. 0.001)
64
64
end
65
65
66
66
let suite =
+17
-31
lib/coordinate.ml
+17
-31
lib/coordinate.ml
···
28
28
(* ------------------------------------------------------------------ *)
29
29
30
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))
31
+
let vec3_length v = sqrt ((v.x *. v.x) +. (v.y *. v.y) +. (v.z *. v.z))
34
32
35
33
let normalize_angle a =
36
34
let a = Float.rem a twopi in
···
41
39
(* ------------------------------------------------------------------ *)
42
40
43
41
let julian_date_of_unix t = (t /. 86400.) +. 2440587.5
44
-
45
42
let unix_of_julian_date jd = (jd -. 2440587.5) *. 86400.
46
43
47
44
let julian_centuries unix_t =
···
50
47
51
48
(** GMST using the IAU 1982 model.
52
49
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. *)
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. *)
57
53
let gmst_of_unix unix_time =
58
54
let t = julian_centuries unix_time in
59
55
let gmst_sec =
60
56
67310.54841
61
-
+. ((876600.0 *. 3600. +. 8640184.812866) *. t)
57
+
+. (((876600.0 *. 3600.) +. 8640184.812866) *. t)
62
58
+. (0.093104 *. t *. t)
63
59
-. (6.2e-6 *. t *. t *. t)
64
60
in
···
70
66
(* Frame transforms *)
71
67
(* ------------------------------------------------------------------ *)
72
68
73
-
(** Z-axis rotation by angle theta.
74
-
Rz(θ) = [[cos θ, sin θ, 0], [-sin θ, cos θ, 0], [0, 0, 1]] *)
69
+
(** Z-axis rotation by angle theta. Rz(θ) =
70
+
[[cos θ, sin θ, 0], [-sin θ, cos θ, 0], [0, 0, 1]] *)
75
71
let rotate_z theta v =
76
72
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
-
}
73
+
{ x = (c *. v.x) +. (s *. v.y); y = (-.s *. v.x) +. (c *. v.y); z = v.z }
82
74
83
75
let teme_to_ecef ~gmst v = rotate_z gmst v
84
76
let ecef_to_teme ~gmst v = rotate_z (-.gmst) v
···
91
83
92
84
(** Bowring's iterative method for ECEF → geodetic.
93
85
94
-
Reference: Vallado Algorithm 12, Montenbruck & Gill §5.3.
95
-
Converges in 2-3 iterations for typical positions. *)
86
+
Reference: Vallado Algorithm 12, Montenbruck & Gill §5.3. Converges in 2-3
87
+
iterations for typical positions. *)
96
88
let ecef_to_geodetic v =
97
89
let a = earth_radius in
98
90
let f = earth_flattening in
99
-
let e2 = (2. *. f) -. (f *. f) in (* first eccentricity squared *)
91
+
let e2 = (2. *. f) -. (f *. f) in
92
+
(* first eccentricity squared *)
100
93
let p = sqrt ((v.x *. v.x) +. (v.y *. v.y)) in
101
94
let lon_rad = atan2 v.y v.x in
102
95
(* Initial latitude estimate *)
···
110
103
let sin_lat = sin !lat_rad in
111
104
let n = a /. sqrt (1. -. (e2 *. sin_lat *. sin_lat)) in
112
105
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))
106
+
if abs_float (cos !lat_rad) > 1e-10 then (p /. cos !lat_rad) -. n
107
+
else (abs_float v.z /. sin_lat) -. (n *. (1. -. e2))
117
108
in
118
-
{
119
-
lat = !lat_rad *. 180. /. Float.pi;
120
-
lon = lon_rad *. 180. /. Float.pi;
121
-
alt;
122
-
}
109
+
{ lat = !lat_rad *. 180. /. Float.pi; lon = lon_rad *. 180. /. Float.pi; alt }
123
110
124
111
(** Geodetic → ECEF.
125
112
···
136
123
{
137
124
x = (n +. g.alt) *. cos_lat *. cos lon_rad;
138
125
y = (n +. g.alt) *. cos_lat *. sin lon_rad;
139
-
z = (n *. (1. -. e2) +. g.alt) *. sin_lat;
126
+
z = ((n *. (1. -. e2)) +. g.alt) *. sin_lat;
140
127
}
141
128
142
129
let teme_to_geodetic ~gmst v = ecef_to_geodetic (teme_to_ecef ~gmst v)
···
148
135
let pp_vec3 ppf v = Fmt.pf ppf "(%.3f, %.3f, %.3f)" v.x v.y v.z
149
136
150
137
let pp_geodetic ppf g =
151
-
Fmt.pf ppf "%.4f°%s %.4f°%s %.3f km"
152
-
(abs_float g.lat)
138
+
Fmt.pf ppf "%.4f°%s %.4f°%s %.3f km" (abs_float g.lat)
153
139
(if g.lat >= 0. then "N" else "S")
154
140
(abs_float g.lon)
155
141
(if g.lon >= 0. then "E" else "W")
+19
-20
lib/coordinate.mli
+19
-20
lib/coordinate.mli
···
1
1
(** Astrodynamics coordinate frame transforms.
2
2
3
-
Provides conversions between standard reference frames used in
4
-
satellite operations:
3
+
Provides conversions between standard reference frames used in satellite
4
+
operations:
5
5
- TEME (True Equator, Mean Equinox) — SGP4 output frame
6
6
- ECEF (Earth-Centered, Earth-Fixed) — rotates with Earth
7
7
- J2000 / EME2000 — inertial frame (≈ GCRF for our precision)
···
20
20
21
21
type geodetic = { lat : float; lon : float; alt : float }
22
22
(** Geodetic coordinates: latitude (deg), longitude (deg), altitude (km).
23
-
Latitude: [-90, +90], north positive.
24
-
Longitude: [-180, +180], east positive. *)
23
+
Latitude: [-90, +90], north positive. Longitude: [-180, +180], east
24
+
positive. *)
25
25
26
26
(** {1 Constants} *)
27
27
···
29
29
(** [earth_radius] is Earth's mean equatorial radius (6378.137 km, WGS-84). *)
30
30
31
31
val earth_flattening : float
32
-
(** [earth_flattening] is Earth's flattening factor (1/298.257223563, WGS-84). *)
32
+
(** [earth_flattening] is Earth's flattening factor (1/298.257223563, WGS-84).
33
+
*)
33
34
34
35
val earth_mu : float
35
36
(** [earth_mu] is Earth's gravitational parameter (398600.4418 km^3/s^2). *)
···
43
44
(** {1 Time} *)
44
45
45
46
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). *)
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). *)
48
49
49
50
val gmst : Ptime.t -> float
50
51
(** [gmst t] is GMST (radians) at time [t]. *)
···
61
62
(** {1 Frame transforms} *)
62
63
63
64
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). *)
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
67
68
68
val ecef_to_teme : gmst:float -> vec3 -> vec3
69
69
(** [ecef_to_teme ~gmst v] is the inverse of {!teme_to_ecef}. *)
70
70
71
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. *)
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. *)
75
74
76
75
val ecef_to_j2000 : gmst:float -> vec3 -> vec3
77
76
(** [ecef_to_j2000 ~gmst v] is the inverse of {!j2000_to_ecef}. *)
···
79
78
(** {1 Geodetic conversions} *)
80
79
81
80
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). *)
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). *)
85
84
86
85
val geodetic_to_ecef : geodetic -> vec3
87
-
(** [geodetic_to_ecef g] converts geodetic coordinates to ECEF (km).
88
-
(Vallado Algorithm 11). *)
86
+
(** [geodetic_to_ecef g] converts geodetic coordinates to ECEF (km). (Vallado
87
+
Algorithm 11). *)
89
88
90
89
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)]. *)
90
+
(** [teme_to_geodetic ~gmst v] converts TEME position to geodetic. Convenience
91
+
for [ecef_to_geodetic (teme_to_ecef ~gmst v)]. *)
93
92
94
93
(** {1 Utility} *)
95
94
+31
-36
test/test_coordinate.ml
+31
-36
test/test_coordinate.ml
···
21
21
(* WGS-84 values: ITG standard *)
22
22
check_float "earth_radius" 1e-3 6378.137 Coordinate.earth_radius;
23
23
check_float "earth_mu" 1e-3 398600.4418 Coordinate.earth_mu;
24
-
check_float "flattening" 1e-12 (1. /. 298.257223563) Coordinate.earth_flattening;
24
+
check_float "flattening" 1e-12 (1. /. 298.257223563)
25
+
Coordinate.earth_flattening;
25
26
check_float "J2" 1e-8 0.00108263 Coordinate.earth_j2
26
27
27
28
(* ================================================================== *)
28
29
(* GMST tests *)
29
30
(* ================================================================== *)
30
31
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. *)
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. *)
34
34
let test_gmst_j2000 () =
35
-
let j2000_unix = 946728000. in (* 2000-01-01T12:00:00 UTC approx *)
35
+
let j2000_unix = 946728000. in
36
+
(* 2000-01-01T12:00:00 UTC approx *)
36
37
let gmst = Coordinate.gmst_of_unix j2000_unix in
37
38
let gmst_deg = rad_to_deg gmst in
38
39
(* Allow 0.5° tolerance for TT/UTC difference *)
39
40
let diff = abs_float (gmst_deg -. 280.46) in
40
41
Alcotest.(check bool) "GMST at J2000 ~280.46°" true (diff < 0.5)
41
42
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π. *)
43
+
(** GMST increases by ~360.985647° per day (sidereal day). After exactly one
44
+
sidereal day (86164.0905s), GMST should advance by exactly 2π. *)
45
45
let test_gmst_sidereal_day () =
46
-
let t0 = 1735732800. in (* 2025-01-01T12:00:00 UTC *)
46
+
let t0 = 1735732800. in
47
+
(* 2025-01-01T12:00:00 UTC *)
47
48
let sidereal_day = 86164.0905 in
48
49
let g0 = Coordinate.gmst_of_unix t0 in
49
50
let g1 = Coordinate.gmst_of_unix (t0 +. sidereal_day) in
···
55
56
(** GMST should be monotonically increasing (mod 2π). *)
56
57
let test_gmst_monotone () =
57
58
let t0 = 946728000. in
58
-
let dt = 3600. in (* 1 hour steps *)
59
+
let dt = 3600. in
60
+
(* 1 hour steps *)
59
61
let prev = ref (Coordinate.gmst_of_unix t0) in
60
62
for i = 1 to 24 do
61
63
let t = t0 +. (Float.of_int i *. dt) in
···
63
65
let delta = g -. !prev in
64
66
(* Allow wrap-around at 2π *)
65
67
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.);
68
+
Alcotest.(check bool) (Fmt.str "monotone step %d" i) true (delta > 0.);
69
69
prev := g
70
70
done
71
71
···
112
112
check_float "y" 1e-6 2000. ecef.y;
113
113
check_float "z" 1e-6 3000. ecef.z
114
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 *)
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 *)
118
117
let test_teme_ecef_90deg () =
119
118
let v = Coordinate.vec3 1000. 0. 0. in
120
119
let gmst = Float.pi /. 2. in
···
145
144
(* Geodetic tests *)
146
145
(* ================================================================== *)
147
146
148
-
(** Point on equator at prime meridian: lat=0, lon=0, alt=0.
149
-
ECEF = (a, 0, 0) where a = earth_radius. *)
147
+
(** Point on equator at prime meridian: lat=0, lon=0, alt=0. ECEF = (a, 0, 0)
148
+
where a = earth_radius. *)
150
149
let test_geodetic_equator_pm () =
151
150
let v = Coordinate.vec3 Coordinate.earth_radius 0. 0. in
152
151
let g = Coordinate.ecef_to_geodetic v in
···
154
153
check_float "lon" 0.01 0. g.lon;
155
154
check_float "alt" 0.01 0. g.alt
156
155
157
-
(** North pole: lat=90, lon=0 (or undefined), alt=0.
158
-
ECEF = (0, 0, b) where b = a*(1-f). *)
156
+
(** North pole: lat=90, lon=0 (or undefined), alt=0. ECEF = (0, 0, b) where b =
157
+
a*(1-f). *)
159
158
let test_geodetic_north_pole () =
160
-
let b =
161
-
Coordinate.earth_radius *. (1. -. Coordinate.earth_flattening)
162
-
in
159
+
let b = Coordinate.earth_radius *. (1. -. Coordinate.earth_flattening) in
163
160
let v = Coordinate.vec3 0. 0. b in
164
161
let g = Coordinate.ecef_to_geodetic v in
165
162
check_float "north pole lat" 0.01 90. g.lat;
···
172
169
check_float "lat" 0.01 0. g.lat;
173
170
check_float "lon" 0.01 90. g.lon
174
171
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. *)
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. *)
178
175
let test_geodetic_washington () =
179
176
let g = { Coordinate.lat = 38.8895; lon = -77.0353; alt = 0.054 } in
180
177
let ecef = Coordinate.geodetic_to_ecef g in
···
215
212
216
213
(** GEO altitude: ~35786 km above equator. *)
217
214
let test_geodetic_geo_altitude () =
218
-
let r = 42164.0 in (* GEO radius *)
215
+
let r = 42164.0 in
216
+
(* GEO radius *)
219
217
let v = Coordinate.vec3 r 0. 0. in
220
218
let g = Coordinate.ecef_to_geodetic v in
221
219
let expected_alt = r -. Coordinate.earth_radius in
222
220
check_float "GEO alt" 5. expected_alt g.alt
223
221
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. *)
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
+
*)
227
225
let test_vallado_3_3 () =
228
226
let v = Coordinate.vec3 6524.834 6862.875 6448.296 in
229
227
let g = Coordinate.ecef_to_geodetic v in
230
228
(* These are approximate — the example has a different reference *)
231
229
let r = Coordinate.vec3_length v in
232
230
Alcotest.(check bool) "high orbit" true (r > 10000.);
233
-
Alcotest.(check bool) "lat reasonable" true
234
-
(g.lat > 20. && g.lat < 50.);
231
+
Alcotest.(check bool) "lat reasonable" true (g.lat > 20. && g.lat < 50.);
235
232
Alcotest.(check bool) "alt positive" true (g.alt > 3000.)
236
233
237
234
(* ================================================================== *)
···
262
259
let test_normalize_angle () =
263
260
let twopi = 2. *. Float.pi in
264
261
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));
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));
269
264
check_float "zero" 1e-10 0. (Coordinate.normalize_angle 0.)
270
265
271
266
(* ================================================================== *)