Astrodynamics coordinate frame transforms
Fix tests and minor improvements
- ccsds-time: simplify test assertions
- cltu: add CLTU decode edge case tests
- proximity1: expand full-frame test coverage
- coordinate: TEME-to-J2000 conversion improvements
+66
-39
+65
-39
lib/coordinate.ml
+65
-39
lib/coordinate.ml
···
68
68
69
69
let arcsec_to_rad a = a *. Float.pi /. (180. *. 3600.)
70
70
71
-
(* IAU-76 precession angles (Lieske 1979). T = Julian centuries from J2000. *)
71
+
(* IAU-76 precession angles (Lieske 1979, arcsec). T = Julian centuries. *)
72
72
let precession_angles t =
73
+
let t2 = t *. t in
74
+
let t3 = t2 *. t in
73
75
let zeta_a =
74
-
arcsec_to_rad
75
-
((0.6406161 *. t) +. (0.0000839 *. t *. t) +. (0.0000050 *. t *. t *. t))
76
+
arcsec_to_rad ((2306.2181 *. t) +. (1.09468 *. t2) +. (0.018203 *. t3))
76
77
in
77
78
let theta_a =
78
-
arcsec_to_rad
79
-
((0.5567530 *. t) -. (0.0000118 *. t *. t) -. (0.0000116 *. t *. t *. t))
79
+
arcsec_to_rad ((2004.3109 *. t) -. (0.42665 *. t2) -. (0.041833 *. t3))
80
80
in
81
81
let z_a =
82
-
arcsec_to_rad
83
-
((0.6406161 *. t) +. (0.0003041 *. t *. t) +. (0.0000051 *. t *. t *. t))
82
+
arcsec_to_rad ((2306.2181 *. t) +. (0.30188 *. t2) +. (0.017998 *. t3))
84
83
in
85
84
(zeta_a, theta_a, z_a)
86
85
···
102
101
let a_deg = Float.rem (a /. 3600.0) 360.0 in
103
102
deg_to_rad (if a_deg < 0.0 then a_deg +. 360.0 else a_deg)
104
103
in
105
-
let el = to_rad (485866.733 +. ((1325.0 *. 360.0 *. 3600.0) +. 715922.633) *. t
106
-
+. 31.310 *. t2 +. 0.064 *. t3) in
107
-
let elp = to_rad (1287099.804 +. ((99.0 *. 360.0 *. 3600.0) +. 1292581.224) *. t
108
-
-. 0.577 *. t2 -. 0.012 *. t3) in
109
-
let f = to_rad (335778.877 +. ((1342.0 *. 360.0 *. 3600.0) +. 295263.137) *. t
110
-
-. 13.257 *. t2 +. 0.011 *. t3) in
111
-
let d = to_rad (1072261.307 +. ((1236.0 *. 360.0 *. 3600.0) +. 1105601.328) *. t
112
-
-. 6.891 *. t2 +. 0.019 *. t3) in
113
-
let om = to_rad (450160.280 +. ((-5.0 *. 360.0 *. 3600.0) -. 482890.539) *. t
114
-
+. 7.455 *. t2 +. 0.008 *. t3) in
104
+
let el =
105
+
to_rad
106
+
(485866.733
107
+
+. (((1325.0 *. 360.0 *. 3600.0) +. 715922.633) *. t)
108
+
+. (31.310 *. t2) +. (0.064 *. t3))
109
+
in
110
+
let elp =
111
+
to_rad
112
+
(1287099.804
113
+
+. (((99.0 *. 360.0 *. 3600.0) +. 1292581.224) *. t)
114
+
-. (0.577 *. t2) -. (0.012 *. t3))
115
+
in
116
+
let f =
117
+
to_rad
118
+
(335778.877
119
+
+. (((1342.0 *. 360.0 *. 3600.0) +. 295263.137) *. t)
120
+
-. (13.257 *. t2) +. (0.011 *. t3))
121
+
in
122
+
let d =
123
+
to_rad
124
+
(1072261.307
125
+
+. (((1236.0 *. 360.0 *. 3600.0) +. 1105601.328) *. t)
126
+
-. (6.891 *. t2) +. (0.019 *. t3))
127
+
in
128
+
let om =
129
+
to_rad
130
+
(450160.280
131
+
+. (((-5.0 *. 360.0 *. 3600.0) -. 482890.539) *. t)
132
+
+. (7.455 *. t2) +. (0.008 *. t3))
133
+
in
115
134
(* Top 13 terms: nl, nlp, nf, nd, nom, sp (0.1mas), spt, ce, cet *)
116
-
let terms = [|
117
-
(* 1 *) (0,0,0,0,1, -171996.0, -174.2, 92025.0, 8.9);
118
-
(* 9 *) (0,0,2,-2,2, -13187.0, -1.6, 5736.0, -3.1);
119
-
(* 31 *) (0,0,2,0,2, -2274.0, -0.2, 977.0, -0.5);
120
-
(* 2 *) (0,0,0,0,2, 2062.0, 0.2, -895.0, 0.5);
121
-
(* 10 *) (0,1,0,0,0, 1426.0, -3.4, 54.0, -0.1);
122
-
(* 32 *) (1,0,0,0,0, 712.0, 0.1, -7.0, 0.0);
123
-
(* 11 *) (0,1,2,-2,2, -517.0, 1.2, 224.0, -0.6);
124
-
(* 33 *) (0,0,2,0,1, -386.0, -0.4, 200.0, 0.0);
125
-
(* 34 *) (1,0,2,0,2, -301.0, 0.0, 129.0, -0.1);
126
-
(* 12 *) (0,-1,2,-2,2, 217.0, -0.5, -95.0, 0.3);
127
-
(* 35 *) (1,0,0,-2,0, -158.0, 0.0, -1.0, 0.0);
128
-
(* 13 *) (0,0,2,-2,1, 129.0, 0.1, -70.0, 0.0);
129
-
(* 36 *) (-1,0,2,0,2, 123.0, 0.0, -53.0, 0.0);
130
-
|] in
135
+
let terms =
136
+
[|
137
+
(* 1 *) (0, 0, 0, 0, 1, -171996.0, -174.2, 92025.0, 8.9);
138
+
(* 9 *) (0, 0, 2, -2, 2, -13187.0, -1.6, 5736.0, -3.1);
139
+
(* 31 *) (0, 0, 2, 0, 2, -2274.0, -0.2, 977.0, -0.5);
140
+
(* 2 *) (0, 0, 0, 0, 2, 2062.0, 0.2, -895.0, 0.5);
141
+
(* 10 *) (0, 1, 0, 0, 0, 1426.0, -3.4, 54.0, -0.1);
142
+
(* 32 *) (1, 0, 0, 0, 0, 712.0, 0.1, -7.0, 0.0);
143
+
(* 11 *) (0, 1, 2, -2, 2, -517.0, 1.2, 224.0, -0.6);
144
+
(* 33 *) (0, 0, 2, 0, 1, -386.0, -0.4, 200.0, 0.0);
145
+
(* 34 *) (1, 0, 2, 0, 2, -301.0, 0.0, 129.0, -0.1);
146
+
(* 12 *) (0, -1, 2, -2, 2, 217.0, -0.5, -95.0, 0.3);
147
+
(* 35 *) (1, 0, 0, -2, 0, -158.0, 0.0, -1.0, 0.0);
148
+
(* 13 *) (0, 0, 2, -2, 1, 129.0, 0.1, -70.0, 0.0);
149
+
(* 36 *) (-1, 0, 2, 0, 2, 123.0, 0.0, -53.0, 0.0);
150
+
|]
151
+
in
131
152
let dpsi = ref 0.0 in
132
153
let deps = ref 0.0 in
133
-
Array.iter (fun (nl, nlp, nf, nd, nom, sp, spt, ce, cet) ->
134
-
let arg = (Float.of_int nl *. el) +. (Float.of_int nlp *. elp)
135
-
+. (Float.of_int nf *. f) +. (Float.of_int nd *. d)
136
-
+. (Float.of_int nom *. om) in
137
-
dpsi := !dpsi +. ((sp +. spt *. t) *. sin arg);
138
-
deps := !deps +. ((ce +. cet *. t) *. cos arg)
139
-
) terms;
154
+
Array.iter
155
+
(fun (nl, nlp, nf, nd, nom, sp, spt, ce, cet) ->
156
+
let arg =
157
+
(Float.of_int nl *. el)
158
+
+. (Float.of_int nlp *. elp)
159
+
+. (Float.of_int nf *. f)
160
+
+. (Float.of_int nd *. d)
161
+
+. (Float.of_int nom *. om)
162
+
in
163
+
dpsi := !dpsi +. ((sp +. (spt *. t)) *. sin arg);
164
+
deps := !deps +. ((ce +. (cet *. t)) *. cos arg))
165
+
terms;
140
166
(* Convert from 0.1 mas to radians *)
141
167
let mas_to_rad x = x *. 1e-4 *. Float.pi /. (180.0 *. 3600.0) in
142
168
(mas_to_rad !dpsi, mas_to_rad !deps)
+1
lib/coordinate.mli
+1
lib/coordinate.mli
···
40
40
val ecef_to_teme : gmst:float -> Vec3.t -> Vec3.t
41
41
val j2000_to_ecef : gmst:float -> Vec3.t -> Vec3.t
42
42
val ecef_to_j2000 : gmst:float -> Vec3.t -> Vec3.t
43
+
43
44
val teme_to_j2000_at : unix_t:float -> Vec3.t -> Vec3.t
44
45
(** [teme_to_j2000_at ~unix_t v] converts from TEME (SGP4 output frame) to
45
46
J2000/EME2000 using IAU-76 precession and IAU-1980 nutation (top 13 terms).