1 #include <cstdlib>
2 using namespace std;
3
4 #include "body.h"
5 #include "xpUtil.h"
6
7 /*
8 GUST86 ephemeris is described in Laskar & Jacobson,
9 Astron. Astrophys. 188, 212-224 (1987)
10
11 Much of this code is translated from the GUST86 FORTRAN code which
12 is at ftp://ftp.bdl.fr/pub/ephem/satel/gust86
13 */
14
15 static double
solveKepler(const double L,const double K,const double H)16 solveKepler(const double L, const double K, const double H)
17 {
18 if (L == 0) return(0);
19
20 double F = L;
21 double E;
22
23 double F0 = L;
24 double E0 = fabs(L);
25
26 const double eps = 1e-16;
27 for (int i = 0; i < 20; i++)
28 {
29 const double SF = sin(F0);
30 const double CF = cos(F0);
31 const double FF0 = F0 - K*SF + H*CF - L;
32 const double FPF0 = 1 - K*CF - H*SF;
33 double SDIR = FF0/FPF0;
34
35 double denom = 1;
36 while (1)
37 {
38 F = F0 - SDIR / denom;
39 E = fabs(F-F0);
40 if (E <= E0) break;
41 denom *= 2;
42 }
43
44 if (denom == 1 && E <= eps && FF0 <= eps) return(F);
45
46 F0 = F;
47 E0 = E;
48 }
49 return(F);
50 }
51
52 static void
calcRectangular(const double N,const double L,const double K,const double H,const double Q,const double P,const double GMS,double & X,double & Y,double & Z)53 calcRectangular(const double N, const double L, const double K,
54 const double H, const double Q, const double P,
55 const double GMS,
56 double &X, double &Y, double &Z)
57 {
58 // Calculate the semi-major axis
59 const double A = pow(GMS/(N*N), 1./3.) / AU_to_km;
60
61 const double PHI = sqrt(1 - K*K - H*H);
62 const double PSI = 1/(1+PHI);
63
64 const double RKI = sqrt(1 - Q*Q - P*P);
65
66 const double F = solveKepler(L, K, H);
67
68 const double SF = sin(F);
69 const double CF = cos(F);
70
71 const double RLMF = -K*SF + H*CF;
72
73 double rot[3][2];
74 rot[0][0] = 1 - 2*P*P;
75 rot[0][1] = 2*P*Q;
76 rot[1][0] = 2*P*Q;
77 rot[1][1] = 1 - 2*Q*Q;
78 rot[2][0] = -2*P*RKI;
79 rot[2][1] = 2*Q*RKI;
80
81 double TX[2];
82 TX[0] = A*(CF - PSI * H * RLMF - K);
83 TX[1] = A*(SF + PSI * K * RLMF - H);
84
85 X = rot[0][0] * TX[0] + rot[0][1] * TX[1];
86 Y = rot[1][0] * TX[0] + rot[1][1] * TX[1];
87 Z = rot[2][0] * TX[0] + rot[2][1] * TX[1];
88 }
89
90 // convert UME50* coordinates to EME50
91 static void
UranicentricToGeocentricEquatorial(double & X,double & Y,double & Z)92 UranicentricToGeocentricEquatorial(double &X, double &Y, double &Z)
93 {
94 const double alpha0 = 76.6067 * deg_to_rad;
95 const double delta0 = 15.0322 * deg_to_rad;
96
97 const double sa = sin(alpha0);
98 const double sd = sin(delta0);
99 const double ca = cos(alpha0);
100 const double cd = cos(delta0);
101
102 const double oldX = X;
103 const double oldY = Y;
104 const double oldZ = Z;
105
106 X = sa * oldX + ca * sd * oldY + ca * cd * oldZ;
107 Y = -ca * oldX + sa * sd * oldY + sa * cd * oldZ;
108 Z = -cd * oldY + sd * oldZ;
109 }
110
111 void
urasat(const double jd,const body b,double & X,double & Y,double & Z)112 urasat(const double jd, const body b,
113 double &X, double &Y, double &Z)
114 {
115 const double t = jd - 2444239.5;
116 const double tcen = t/365.25;
117
118 const double N1 = fmod(4.445190550 * t - 0.238051, TWO_PI);
119 const double N2 = fmod(2.492952519 * t + 3.098046, TWO_PI);
120 const double N3 = fmod(1.516148111 * t + 2.285402, TWO_PI);
121 const double N4 = fmod(0.721718509 * t + 0.856359, TWO_PI);
122 const double N5 = fmod(0.466692120 * t - 0.915592, TWO_PI);
123
124 const double E1 = (20.082 * deg_to_rad * tcen + 0.611392);
125 const double E2 = ( 6.217 * deg_to_rad * tcen + 2.408974);
126 const double E3 = ( 2.865 * deg_to_rad * tcen + 2.067774);
127 const double E4 = ( 2.078 * deg_to_rad * tcen + 0.735131);
128 const double E5 = ( 0.386 * deg_to_rad * tcen + 0.426767);
129
130 const double I1 = (-20.309 * deg_to_rad * tcen + 5.702313);
131 const double I2 = ( -6.288 * deg_to_rad * tcen + 0.395757);
132 const double I3 = ( -2.836 * deg_to_rad * tcen + 0.589326);
133 const double I4 = ( -1.843 * deg_to_rad * tcen + 1.746237);
134 const double I5 = ( -0.259 * deg_to_rad * tcen + 4.206896);
135
136 const double GM1 = 4.4;
137 const double GM2 = 86.1;
138 const double GM3 = 84.0;
139 const double GM4 = 230.0;
140 const double GM5 = 200.0;
141 const double GMU = 5794554.5 - (GM1 + GM2 + GM3 + GM4 + GM5);
142
143 double N = 0, L = 0, K = 0, H = 0, Q = 0, P = 0, GMS = 0;
144
145 switch (b)
146 {
147 case MIRANDA:
148 N = (4443522.67
149 - 34.92 * cos(N1 - 3*N2 + 2*N3)
150 + 8.47 * cos(2*N1 - 6*N2 + 4*N3)
151 + 1.31 * cos(3*N1 - 9*N2 + 6*N3)
152 - 52.28 * cos(N1 - N2)
153 -136.65 * cos(2*N1 - 2*N2)) * 1e-6;
154
155 L = (-238051.58
156 + 4445190.55 * t
157 + 25472.17 * sin(N1 - 3*N2 + 2*N3)
158 - 3088.31 * sin(2*N1 - 6*N2 + 4*N3)
159 - 318.10 * sin(3*N1 - 9*N2 + 6*N3)
160 - 37.49 * sin(4*N1 - 12*N2 + 8*N3)
161 - 57.85 * sin(N1 - N2)
162 - 62.32 * sin(2*N1 - 2*N2)
163 - 27.95 * sin(3*N1 - 3*N2)) * 1e-6;
164
165 K = (1312.38 * cos(E1)
166 + 71.81 * cos(E2)
167 + 69.77 * cos(E3)
168 + 6.75 * cos(E4)
169 + 6.27 * cos(E5)
170 - 123.31 * cos(-N1 + 2*N2)
171 + 39.52 * cos(-2*N1 + 3*N2)
172 + 194.10 * cos(N1)) * 1e-6;
173
174 H = (1312.38 * sin(E1)
175 + 71.81 * sin(E2)
176 + 69.77 * sin(E3)
177 + 6.75 * sin(E4)
178 + 6.27 * sin(E5)
179 - 123.31 * sin(-N1 + 2*N2)
180 + 39.52 * sin(-2*N1 + 3*N2)
181 + 194.10 * sin(N1)) * 1e-6;
182
183 Q = (37871.71 * cos(I1)
184 + 27.01 * cos(I2)
185 + 30.76 * cos(I3)
186 + 12.18 * cos(I4)
187 + 5.37 * cos(I5)) * 1e-6;
188
189 P = (37871.71 * sin(I1)
190 + 27.01 * sin(I2)
191 + 30.76 * sin(I3)
192 + 12.18 * sin(I4)
193 + 5.37 * sin(I5)) * 1e-6;
194
195 GMS = GMU + GM1;
196 break;
197 case ARIEL:
198 N = (2492542.57
199 + 2.55 * cos(N1 - 3*N2 + 2*N3)
200 - 42.16 * cos(N2 - N3)
201 - 102.56 * cos(2*N2 - 2*N3)) * 1e-6;
202
203 L = (3098046.41
204 + 2492952.52 * t
205 - 1860.50 * sin(N1 - 3*N2 + 2*N3)
206 + 219.99 * sin(2*N1 - 6*N2 + 4*N3)
207 + 23.10 * sin(3*N1 - 9*N2 + 6*N3)
208 + 4.30 * sin(4*N1 - 12*N2 + 8*N3)
209 - 90.11 * sin(N2 - N3)
210 - 91.07 * sin(2*(N2 - N3))
211 - 42.75 * sin(3*(N2 - N3))
212 - 16.49 * sin(2*(N2 - N4))) * 1e-6;
213
214 K = (- 3.35 * cos(E1)
215 + 1187.63 * cos(E2)
216 + 861.59 * cos(E3)
217 + 71.50 * cos(E4)
218 + 55.59 * cos(E5)
219 - 84.60 * cos(-N2 + 2*N3)
220 + 91.81 * cos(-2*N2 + 3*N3)
221 + 20.03 * cos(-N2 + 2*N4)
222 + 89.77 * cos(N2)) * 1e-6;
223
224 H = (- 3.35 * sin(E1)
225 + 1187.63 * sin(E2)
226 + 861.59 * sin(E3)
227 + 71.50 * sin(E4)
228 + 55.59 * sin(E5)
229 - 84.60 * sin(-N2 + 2*N3)
230 + 91.81 * sin(-2*N2 + 3*N3)
231 + 20.03 * sin(-N2 + 2*N4)
232 + 89.77 * sin(N2)) * 1e-6;
233
234 Q = (- 121.75 * cos(I1)
235 + 358.25 * cos(I2)
236 + 290.08 * cos(I3)
237 + 97.78 * cos(I4)
238 + 33.97 * cos(I5)) * 1e-6;
239
240 P = (- 121.75 * sin(I1)
241 + 358.25 * sin(I2)
242 + 290.08 * sin(I3)
243 + 97.78 * sin(I4)
244 + 33.97 * sin(I5)) * 1e-6;
245
246 GMS = GMU + GM2;
247 break;
248 case UMBRIEL:
249 N = (1515954.90
250 + 9.74 * cos(N3 - 2*N4 + E3)
251 - 106.00 * cos(N2 - N3)
252 + 54.16 * cos(2*(N2 - N3))
253 - 23.59 * cos(N3 - N4)
254 - 70.70 * cos(2*(N3 - N4))
255 - 36.28 * cos(3*(N3 - N4))) * 1e-6;
256
257 L = (2285401.69
258 + 1516148.11 * t
259 + 660.57 * sin( N1 - 3*N2 + 2*N3)
260 - 76.51 * sin(2*N1 - 6*N2 + 4*N3)
261 - 8.96 * sin(3*N1 - 9*N2 + 6*N3)
262 - 2.53 * sin(4*N1 - 12*N2 + 8*N3)
263 - 52.91 * sin(N3 - 4*N4 + 3*N5)
264 - 7.34 * sin(N3 - 2*N4 + E5)
265 - 1.83 * sin(N3 - 2*N4 + E4)
266 + 147.91 * sin(N3 - 2*N4 + E3)
267 - 7.77 * sin(N3 - 2*N4 + E2)
268 + 97.76 * sin(N2 - N3)
269 + 73.13 * sin(2*(N2 - N3))
270 + 34.71 * sin(3*(N2 - N3))
271 + 18.89 * sin(4*(N2 - N3))
272 - 67.89 * sin(N3 - N4)
273 - 82.86 * sin(2*(N3 - N4))
274 - 33.81 * sin(3*(N3 - N4))
275 - 15.79 * sin(4*(N3 - N4))
276 - 10.21 * sin(N3 - N5)
277 - 17.08 * sin(2*(N3 - N5))) * 1e-6;
278
279 K = (- 0.21 * cos(E1)
280 - 227.95 * cos(E2)
281 + 3904.69 * cos(E3)
282 + 309.17 * cos(E4)
283 + 221.92 * cos(E5)
284 + 29.34 * cos(N2)
285 + 26.20 * cos(N3)
286 + 51.19 * cos(-N2+2*N3)
287 - 103.86 * cos(-2*N2+3*N3)
288 - 27.16 * cos(-3*N2+4*N3)
289 - 16.22 * cos(N4)
290 + 549.23 * cos(-N3 + 2*N4)
291 + 34.70 * cos(-2*N3 + 3*N4)
292 + 12.81 * cos(-3*N3 + 4*N4)
293 + 21.81 * cos(-N3 + 2*N5)
294 + 46.25 * cos(N3)) * 1e-6;
295
296 H = (- 0.21 * sin(E1)
297 - 227.95 * sin(E2)
298 + 3904.69 * sin(E3)
299 + 309.17 * sin(E4)
300 + 221.92 * sin(E5)
301 + 29.34 * sin(N2)
302 + 26.20 * sin(N3)
303 + 51.19 * sin(-N2+2*N3)
304 - 103.86 * sin(-2*N2+3*N3)
305 - 27.16 * sin(-3*N2+4*N3)
306 - 16.22 * sin(N4)
307 + 549.23 * sin(-N3 + 2*N4)
308 + 34.70 * sin(-2*N3 + 3*N4)
309 + 12.81 * sin(-3*N3 + 4*N4)
310 + 21.81 * sin(-N3 + 2*N5)
311 + 46.25 * sin(N3)) * 1e-6;
312
313 Q = (- 10.86 * cos(I1)
314 - 81.51 * cos(I2)
315 + 1113.36 * cos(I3)
316 + 350.14 * cos(I4)
317 + 106.50 * cos(I5)) * 1e-6;
318
319 P = (- 10.86 * sin(I1)
320 - 81.51 * sin(I2)
321 + 1113.36 * sin(I3)
322 + 350.14 * sin(I4)
323 + 106.50 * sin(I5)) * 1e-6;
324
325 GMS = GMU + GM3;
326 break;
327 case TITANIA:
328 N = (721663.16
329 - 2.64 * cos(N3 - 2*N4 + E3)
330 - 2.16 * cos(2*N4 - 3*N5 + E5)
331 + 6.45 * cos(2*N4 - 3*N5 + E4)
332 - 1.11 * cos(2*N4 - 3*N5 + E3)
333 - 62.23 * cos(N2 - N4)
334 - 56.13 * cos(N3 - N4)
335 - 39.94 * cos(N4 - N5)
336 - 91.85 * cos(2*(N4 - N5))
337 - 58.31 * cos(3*(N4 - N5))
338 - 38.60 * cos(4*(N4 - N5))
339 - 26.18 * cos(5*(N4 - N5))
340 - 18.06 * cos(6*(N4 - N5))) * 1e-6;
341
342 L = (856358.79
343 + 721718.51 * t
344 + 20.61 * sin(N3 - 4*N4 + 3*N5)
345 - 2.07 * sin(N3 - 2*N4 + E5)
346 - 2.88 * sin(N3 - 2*N4 + E4)
347 - 40.79 * sin(N3 - 2*N4 + E3)
348 + 2.11 * sin(N3 - 2*N4 + E2)
349 - 51.83 * sin(2*N4 - 3*N5 + E5)
350 + 159.87 * sin(2*N4 - 3*N5 + E4)
351 - 35.05 * sin(2*N4 - 3*N5 + E3)
352 - 1.56 * sin(3*N4 - 4*N5 + E5)
353 + 40.54 * sin(N2 - N4)
354 + 46.17 * sin(N3 - N4)
355 - 317.76 * sin(N4 - N5)
356 - 305.59 * sin(2*(N4 - N5))
357 - 148.36 * sin(3*(N4 - N5))
358 - 82.92 * sin(4*(N4 - N5))
359 - 49.98 * sin(5*(N4 - N5))
360 - 31.56 * sin(6*(N4 - N5))
361 - 20.56 * sin(7*(N4 - N5))
362 - 13.69 * sin(8*(N4 - N5))) * 1e-6;
363
364 K = (- 0.02 * cos(E1)
365 - 1.29 * cos(E2)
366 - 324.51 * cos(E3)
367 + 932.81 * cos(E4)
368 + 1120.89 * cos(E5)
369 + 33.86 * cos(N2)
370 + 17.46 * cos(N4)
371 + 16.58 * cos(-N2 + 2*N4)
372 + 28.89 * cos(N3)
373 - 35.86 * cos(-N3 + 2*N4)
374 - 17.86 * cos(N4)
375 - 32.10 * cos(N5)
376 - 177.83 * cos(-N4 + 2*N5)
377 + 793.43 * cos(-2*N4 + 3*N5)
378 + 99.48 * cos(-3*N4 + 4*N5)
379 + 44.83 * cos(-4*N4 + 5*N5)
380 + 25.13 * cos(-5*N4 + 6*N5)
381 + 15.43 * cos(-6*N4 + 7*N5)) * 1e-6;
382
383 H = (- 0.02 * sin(E1)
384 - 1.29 * sin(E2)
385 - 324.51 * sin(E3)
386 + 932.81 * sin(E4)
387 + 1120.89 * sin(E5)
388 + 33.86 * sin(N2)
389 + 17.46 * sin(N4)
390 + 16.58 * sin(-N2 + 2*N4)
391 + 28.89 * sin(N3)
392 - 35.86 * sin(-N3 + 2*N4)
393 - 17.86 * sin(N4)
394 - 32.10 * sin(N5)
395 - 177.83 * sin(-N4 + 2*N5)
396 + 793.43 * sin(-2*N4 + 3*N5)
397 + 99.48 * sin(-3*N4 + 4*N5)
398 + 44.83 * sin(-4*N4 + 5*N5)
399 + 25.13 * sin(-5*N4 + 6*N5)
400 + 15.43 * sin(-6*N4 + 7*N5)) * 1e-6;
401
402 Q = (- 1.43 * cos(I1)
403 - 1.06 * cos(I2)
404 - 140.13 * cos(I3)
405 + 685.72 * cos(I4)
406 + 378.32 * cos(I5)) * 1e-6;
407
408 P = (- 1.43 * sin(I1)
409 - 1.06 * sin(I2)
410 - 140.13 * sin(I3)
411 + 685.72 * sin(I4)
412 + 378.32 * sin(I5)) * 1e-6;
413
414 GMS = GMU + GM4;
415 break;
416 case OBERON:
417 N = (466580.54
418 + 2.08 * cos(2*N4 - 3*N5 + E5)
419 - 6.22 * cos(2*N4 - 3*N5 + E4)
420 + 1.07 * cos(2*N4 - 3*N5 + E3)
421 - 43.10 * cos(N2 - N5)
422 - 38.94 * cos(N3 - N5)
423 - 80.11 * cos(N4 - N5)
424 + 59.06 * cos(2*(N4 - N5))
425 + 37.49 * cos(3*(N4 - N5))
426 + 24.82 * cos(4*(N4 - N5))
427 + 16.84 * cos(5*(N4 - N5))) * 1e-6;
428
429 L = (-915591.80
430 + 466692.12 * t
431 - 7.82 * sin(N3 - 4*N4 + 3*N5)
432 + 51.29 * sin(2*N4 - 3*N5 + E5)
433 - 158.24 * sin(2*N4 - 3*N5 + E4)
434 + 34.51 * sin(2*N4 - 3*N5 + E3)
435 + 47.51 * sin(N2 - N5)
436 + 38.96 * sin(N3 - N5)
437 + 359.73 * sin(N4 - N5)
438 + 282.78 * sin(2*(N4 - N5))
439 + 138.60 * sin(3*(N4 - N5))
440 + 78.03 * sin(4*(N4 - N5))
441 + 47.29 * sin(5*(N4 - N5))
442 + 30.00 * sin(6*(N4 - N5))
443 + 19.62 * sin(7*(N4 - N5))
444 + 13.11 * sin(8*(N4 - N5))) * 1e-6;
445
446 K = ( 0 * cos(E1)
447 - 0.35 * cos(E2)
448 + 74.53 * cos(E3)
449 - 758.68 * cos(E4)
450 + 1397.34 * cos(E5)
451 + 39.00 * cos(N2)
452 + 17.66 * cos(-N2 + 2*N5)
453 + 32.42 * cos(N3)
454 + 79.75 * cos(N4)
455 + 75.66 * cos(N5)
456 + 134.04 * cos(-N4 + 2*N5)
457 - 987.26 * cos(-2*N4 + 3*N5)
458 - 126.09 * cos(-3*N4 + 4*N5)
459 - 57.42 * cos(-4*N4 + 5*N5)
460 - 32.41 * cos(-5*N4 + 6*N5)
461 - 19.99 * cos(-6*N4 + 7*N5)
462 - 12.94 * cos(-7*N4 + 8*N5)) * 1e-6;
463
464 H = (0 * sin(E1)
465 - 0.35 * sin(E2)
466 + 74.53 * sin(E3)
467 - 758.68 * sin(E4)
468 + 1397.34 * sin(E5)
469 + 39.00 * sin(N2)
470 + 17.66 * sin(-N2 + 2*N5)
471 + 32.42 * sin(N3)
472 + 79.75 * sin(N4)
473 + 75.66 * sin(N5)
474 + 134.04 * sin(-N4 + 2*N5)
475 - 987.26 * sin(-2*N4 + 3*N5)
476 - 126.09 * sin(-3*N4 + 4*N5)
477 - 57.42 * sin(-4*N4 + 5*N5)
478 - 32.41 * sin(-5*N4 + 6*N5)
479 - 19.99 * sin(-6*N4 + 7*N5)
480 - 12.94 * sin(-7*N4 + 8*N5)) * 1e-6;
481
482 Q = (- 0.44 * cos(I1)
483 - 0.31 * cos(I2)
484 + 36.89 * cos(I3)
485 - 596.33 * cos(I4)
486 + 451.69 * cos(I5)) * 1e-6;
487
488 P = (- 0.44 * sin(I1)
489 - 0.31 * sin(I2)
490 + 36.89 * sin(I3)
491 - 596.33 * sin(I4)
492 + 451.69 * sin(I5)) * 1e-6;
493
494 GMS = GMU + GM5;
495 break;
496 default:
497 xpExit("Unknown Uranus satellite\n", __FILE__, __LINE__);
498 }
499
500 X = 0;
501 Y = 0;
502 Z = 0;
503
504 N /= 86400;
505 L = fmod(L, TWO_PI);
506
507 calcRectangular(N, L, K, H, Q, P, GMS, X, Y, Z);
508
509 UranicentricToGeocentricEquatorial(X, Y, Z);
510
511 // precess to J2000
512 precessB1950J2000(X, Y, Z);
513 }
514