source: trunk/MagicSoft/slalib/planet.c@ 2572

Last change on this file since 2572 was 731, checked in by tbretz, 25 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 29.1 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaPlanet ( double date, int np, double pv[6], int *jstat )
4/*
5** - - - - - - - - - -
6** s l a P l a n e t
7** - - - - - - - - - -
8**
9** Approximate heliocentric position and velocity of a specified
10** major planet.
11**
12** Given:
13** date double TDB (loosely ET) as a Modified Julian Date
14** (JD-2400000.5)
15** np int planet (1=Mercury, 2=Venus, 3=EMB, ...
16** ... 9=Pluto)
17**
18** Returned:
19** pv double[6] heliocentric x,y,z,xdot,ydot,zdot, J2000
20** equatorial triad (AU,AU/s)
21**
22** *jstat int status: +1 = warning: date outside 1000-3000
23** *jstat int status: 0 = OK
24** -1 = illegal NP (outside 1-9)
25** -2 = solution didn't converge
26**
27** Called: slaPlanel
28**
29** Notes
30**
31** 1 The epoch, date, is in the TDB timescale and is a Modified
32** Julian Date (JD-2400000.5).
33**
34** 2 The reference frame is equatorial and is with respect to the
35** mean equinox and ecliptic of epoch J2000.
36**
37** 3 If an np value outside the range 1-9 is supplied, an error
38** status (jstat = -1) is returned and the pv vector set to zeroes.
39**
40** 4 The algorithm for obtaining the mean elements of the planets
41** from Mercury to Neptune is due to J.L. Simon, P. Bretagnon,
42** J. Chapront, M. Chapront-Touze, G. Francou and J. Laskar
43** (Bureau des Longitudes, Paris). The (completely different)
44** algorithm for calculating the ecliptic coordinates of Pluto
45** is by Meeus.
46**
47** 5 Comparisons of the present routine with the JPL DE200 ephemeris
48** give the following RMS errors over the interval 1960-2025:
49**
50** position (km) speed (metre/sec)
51**
52** Mercury 334 0.437
53** Venus 1060 0.855
54** EMB 2010 0.815
55** Mars 7690 1.98
56** Jupiter 71700 7.70
57** Saturn 199000 19.4
58** Uranus 564000 16.4
59** Neptune 158000 14.4
60**
61** From comparisons with DE102, Simon et al quote the following
62** longitude accuracies over the interval 1800-2200:
63**
64** Mercury 4"
65** Venus 5"
66** EMB 6"
67** Mars 17"
68** Jupiter 71"
69** Saturn 81"
70** Uranus 86"
71** Neptune 11"
72**
73** In the case of Pluto, Meeus quotes an accuracy of 0.6 arcsec
74** in longitude and 0.2 arcsec in latitude for the period
75** 1885-2099.
76**
77** For all except Pluto, over the period 1000-3000 the accuracy
78** is better than 1.5 times that over 1800-2200. Outside the
79** period 1000-3000 the accuracy declines. For Pluto the
80** accuracy declines rapidly outside the period 1885-2099.
81** Outside these ranges (1885-2099 for Pluto, 1000-3000 for
82** the rest) a "date out of range" warning status (JSTAT=+1)
83** is returned.
84**
85** 6 The algorithms for (i) Mercury through Neptune and (ii) Pluto
86** are completely independent. In the Mercury through Neptune
87** case, the present SLALIB C implementation follows the original
88** Simon et al Fortran code closely, and delivers essentially
89** the same results. The changes are these:
90**
91** * The date is supplied as a Modified Julian Date rather
92** than a Julian Date (MJD = JD - 2400000.5).
93**
94** * The result is returned only in equatorial Cartesian form;
95** the ecliptic longitude, latitude and radius vector are not
96** returned.
97**
98** * The velocity is in AU per second, not AU per day.
99**
100** * Different error/warning status values are used.
101**
102** * Kepler's equation is not solved inline.
103**
104** * Polynomials in T are nested to minimize rounding errors.
105**
106** * Explicit double-precision constants are used to avoid
107** mixed-mode expressions.
108**
109** 7 For np=3 the result is for the Earth-Moon Barycentre. To
110** obtain the heliocentric position and velocity of the Earth,
111** either use the SLALIB routine slaEvp or use slaDmoon and
112** subtract 0.012150581 times the geocentric Moon vector from
113** the EMB vector produced by the present routine. (The Moon
114** vector should be precessed to J2000 first, but this can
115** be omitted for modern epochs without introducing significant
116** inaccuracy.)
117**
118** References: Simon et al., Astron. Astrophys. 282, 663 (1994).
119** Meeus, Astronomical Algorithms, Willmann-Bell (1991).
120**
121** Defined in slamac.h: D2PI, DAS2R, DD2R, dmod
122**
123** Last revision: 27 May 1997
124**
125** Copyright P.T.Wallace. All rights reserved.
126*/
127
128/* Gaussian gravitational constant (exact) */
129#define GCON 0.01720209895
130
131/* Canonical days to seconds */
132#define CD2S ( GCON / 86400.0 )
133
134/* Seconds per Julian century */
135#define SPC ( 36525.0 * 86400.0 )
136
137/* Sin and cos of J2000 mean obliquity (IAU 1976) */
138#define SE 0.3977771559319137
139#define CE 0.9174820620691818
140
141{
142 int ip, i, j;
143 double t, da, de, dpe, di, dom, dmu, arga, argl, dm,
144 dj, ds, dp, wlbr[3], wlbrd[3],
145 wj, ws, wp, al, ald, sal, cal,
146 ac, bc, dl, dld, db, dbd, dr, drd,
147 sl, cl, sb, cb, slcb, clcb, x, y, z, xd, yd, zd;
148
149/*
150** -----------------------
151** Mercury through Neptune
152** -----------------------
153*/
154
155/* Planetary inverse masses */
156 static double amas[] = {
157 6023600.0,
158 408523.5,
159 328900.5,
160 3098710.0,
161 1047.355,
162 3498.5,
163 22869.0,
164 19314.0
165 };
166
167/*
168** Tables giving the mean Keplerian elements, limited to T^2 terms:
169**
170** a semi-major axis (AU)
171** dlm mean longitude (degree and arcsecond)
172** e eccentricity
173** pi longitude of the perihelion (degree and arcsecond)
174** dinc inclination (degree and arcsecond)
175** omega longitude of the ascending node (degree and arcsecond)
176*/
177 static double a[8][3] = {
178 { 0.3870983098, 0.0, 0.0 },
179 { 0.7233298200, 0.0, 0.0 },
180 { 1.0000010178, 0.0, 0.0 },
181 { 1.5236793419, 3e-10, 0.0 },
182 { 5.2026032092, 19132e-10, -39e-10 },
183 { 9.5549091915, -0.0000213896, 444e-10 },
184 { 19.2184460618, -3716e-10, 979e-10 },
185 { 30.1103868694, -16635e-10, 686e-10 }
186 };
187 static double dlm[8][3] = {
188 { 252.25090552, 5381016286.88982, -1.92789 },
189 { 181.97980085, 2106641364.33548, 0.59381 },
190 { 100.46645683, 1295977422.83429, -2.04411 },
191 { 355.43299958, 689050774.93988, 0.94264 },
192 { 34.35151874, 109256603.77991, -30.60378 },
193 { 50.07744430, 43996098.55732, 75.61614 },
194 { 314.05500511, 15424811.93933, -1.75083 },
195 { 304.34866548, 7865503.20744, 0.21103 }
196 };
197 static double e[8][3] = {
198 { 0.2056317526, 0.0002040653, -28349e-10 },
199 { 0.0067719164, -0.0004776521, 98127e-10 },
200 { 0.0167086342, -0.0004203654, -0.0000126734 },
201 { 0.0934006477, 0.0009048438, -80641e-10 },
202 { 0.0484979255, 0.0016322542, -0.0000471366 },
203 { 0.0555481426, -0.0034664062, -0.0000643639 },
204 { 0.0463812221, -0.0002729293, 0.0000078913 },
205 { 0.0094557470, 0.0000603263, 0.0 }
206 };
207 static double pi[8][3] = {
208 { 77.45611904, 5719.11590, -4.83016 },
209 { 131.56370300, 175.48640, -498.48184 },
210 { 102.93734808, 11612.35290, 53.27577 },
211 { 336.06023395, 15980.45908, -62.32800 },
212 { 14.33120687, 7758.75163, 259.95938 },
213 { 93.05723748, 20395.49439, 190.25952 },
214 { 173.00529106, 3215.56238, -34.09288 },
215 { 48.12027554, 1050.71912, 27.39717 }
216 };
217 static double dinc[8][3] = {
218 { 7.00498625, -214.25629, 0.28977 },
219 { 3.39466189, -30.84437, -11.67836 },
220 { 0.0, 469.97289, -3.35053 },
221 { 1.84972648, -293.31722, -8.11830 },
222 { 1.30326698, -71.55890, 11.95297 },
223 { 2.48887878, 91.85195, -17.66225 },
224 { 0.77319689, -60.72723, 1.25759 },
225 { 1.76995259, 8.12333, 0.08135 }
226 };
227 static double omega[8][3] = {
228 { 48.33089304, -4515.21727, -31.79892 },
229 { 76.67992019, -10008.48154, -51.32614 },
230 { 174.87317577, -8679.27034, 15.34191 },
231 { 49.55809321, -10620.90088, -230.57416 },
232 { 100.46440702, 6362.03561, 326.52178 },
233 { 113.66550252, -9240.19942, -66.23743 },
234 { 74.00595701, 2669.15033, 145.93964 },
235 { 131.78405702, -221.94322, -0.78728 }
236 };
237
238/*
239** Tables for trigonometric terms to be added to the mean elements
240** of the semi-major axes.
241*/
242 static double dkp[8][9] = {
243 { 69613.0, 75645.0, 88306.0, 59899.0, 15746.0, 71087.0,
244 142173.0, 3086.0, 0.0 },
245 { 21863.0, 32794.0, 26934.0, 10931.0, 26250.0, 43725.0,
246 53867.0, 28939.0, 0.0 },
247 { 16002.0, 21863.0, 32004.0, 10931.0, 14529.0, 16368.0,
248 15318.0, 32794.0, 0.0 },
249 { 6345.0, 7818.0, 15636.0, 7077.0, 8184.0, 14163.0,
250 1107.0, 4872.0, 0.0 },
251 { 1760.0, 1454.0, 1167.0, 880.0, 287.0, 2640.0,
252 19.0, 2047.0, 1454.0 },
253 { 574.0, 0.0, 880.0, 287.0, 19.0, 1760.0,
254 1167.0, 306.0, 574.0 },
255 { 204.0, 0.0, 177.0, 1265.0, 4.0, 385.0,
256 200.0, 208.0, 204.0 },
257 { 0.0, 102.0, 106.0, 4.0, 98.0, 1367.0,
258 487.0, 204.0, 0.0 }
259 };
260 static double ca[8][9] = {
261 { 4.0, -13.0, 11.0, -9.0, -9.0, -3.0,
262 -1.0, 4.0, 0.0 },
263 { -156.0, 59.0, -42.0, 6.0, 19.0, -20.0,
264 -10.0, -12.0, 0.0 },
265 { 64.0, -152.0, 62.0, -8.0, 32.0, -41.0,
266 19.0, -11.0, 0.0 },
267 { 124.0, 621.0, -145.0, 208.0, 54.0, -57.0,
268 30.0, 15.0, 0.0 },
269 { -23437.0, -2634.0, 6601.0, 6259.0, -1507.0, -1821.0,
270 2620.0, -2115.0,-1489.0 },
271 { 62911.0,-119919.0, 79336.0, 17814.0,-24241.0, 12068.0,
272 8306.0, -4893.0, 8902.0 },
273 { 389061.0,-262125.0,-44088.0, 8387.0,-22976.0, -2093.0,
274 -615.0, -9720.0, 6633.0 },
275 { -412235.0,-157046.0,-31430.0, 37817.0, -9740.0, -13.0,
276 -7449.0, 9644.0, 0.0 }
277 };
278 static double sa[8][9] = {
279 { -29.0, -1.0, 9.0, 6.0, -6.0, 5.0,
280 4.0, 0.0, 0.0 },
281 { -48.0, -125.0, -26.0, -37.0, 18.0, -13.0,
282 -20.0, -2.0, 0.0 },
283 { -150.0, -46.0, 68.0, 54.0, 14.0, 24.0,
284 -28.0, 22.0, 0.0 },
285 { -621.0, 532.0, -694.0, -20.0, 192.0, -94.0,
286 71.0, -73.0, 0.0 },
287 { -14614.0,-19828.0, -5869.0, 1881.0, -4372.0, -2255.0,
288 782.0, 930.0, 913.0 },
289 { 139737.0, 0.0, 24667.0, 51123.0, -5102.0, 7429.0,
290 -4095.0, -1976.0,-9566.0 },
291 { -138081.0, 0.0, 37205.0,-49039.0,-41901.0,-33872.0,
292 -27037.0,-12474.0,18797.0 },
293 { 0.0, 28492.0,133236.0, 69654.0, 52322.0,-49577.0,
294 -26430.0, -3593.0, 0.0 }
295 };
296
297/*
298** Tables giving the trigonometric terms to be added to the mean
299** elements of the mean longitudes.
300*/
301 static double dkq[8][10] = {
302 { 3086.0, 15746.0, 69613.0, 59899.0, 75645.0,
303 88306.0, 12661.0, 2658.0, 0.0, 0.0 },
304 { 21863.0, 32794.0, 10931.0, 73.0, 4387.0,
305 26934.0, 1473.0, 2157.0, 0.0, 0.0 },
306 { 10.0, 16002.0, 21863.0, 10931.0, 1473.0,
307 32004.0, 4387.0, 73.0, 0.0, 0.0 },
308 { 10.0, 6345.0, 7818.0, 1107.0, 15636.0,
309 7077.0, 8184.0, 532.0, 10.0, 0.0 },
310 { 19.0, 1760.0, 1454.0, 287.0, 1167.0,
311 880.0, 574.0, 2640.0, 19.0,1454.0 },
312 { 19.0, 574.0, 287.0, 306.0, 1760.0,
313 12.0, 31.0, 38.0, 19.0, 574.0 },
314 { 4.0, 204.0, 177.0, 8.0, 31.0,
315 200.0, 1265.0, 102.0, 4.0, 204.0 },
316 { 4.0, 102.0, 106.0, 8.0, 98.0,
317 1367.0, 487.0, 204.0, 4.0, 102.0 }
318 };
319 static double clo[8][10] = {
320 { 21.0, -95.0, -157.0, 41.0, -5.0,
321 42.0, 23.0, 30.0, 0.0, 0.0 },
322 { -160.0, -313.0, -235.0, 60.0, -74.0,
323 -76.0, -27.0, 34.0, 0.0, 0.0 },
324 { -325.0, -322.0, -79.0, 232.0, -52.0,
325 97.0, 55.0, -41.0, 0.0, 0.0 },
326 { 2268.0, -979.0, 802.0, 602.0, -668.0,
327 -33.0, 345.0, 201.0, -55.0, 0.0 },
328 { 7610.0, -4997.0, -7689.0, -5841.0, -2617.0,
329 1115.0, -748.0, -607.0, 6074.0, 354.0 },
330 { -18549.0, 30125.0, 20012.0, -730.0, 824.0,
331 23.0, 1289.0, -352.0,-14767.0,-2062.0 },
332 { -135245.0, -14594.0, 4197.0, -4030.0, -5630.0,
333 -2898.0, 2540.0, -306.0, 2939.0, 1986.0 },
334 { 89948.0, 2103.0, 8963.0, 2695.0, 3682.0,
335 1648.0, 866.0, -154.0, -1963.0, -283.0 }
336 };
337 static double slo[8][10] = {
338 { -342.0, 136.0, -23.0, 62.0, 66.0,
339 -52.0, -33.0, 17.0, 0.0, 0.0 },
340 { 524.0, -149.0, -35.0, 117.0, 151.0,
341 122.0, -71.0, -62.0, 0.0, 0.0 },
342 { -105.0, -137.0, 258.0, 35.0, -116.0,
343 -88.0, -112.0, -80.0, 0.0, 0.0 },
344 { 854.0, -205.0, -936.0, -240.0, 140.0,
345 -341.0, -97.0, -232.0, 536.0, 0.0 },
346 { -56980.0, 8016.0, 1012.0, 1448.0, -3024.0,
347 -3710.0, 318.0, 503.0, 3767.0, 577.0 },
348 { 138606.0, -13478.0, -4964.0, 1441.0, -1319.0,
349 -1482.0, 427.0, 1236.0, -9167.0, -1918.0 },
350 { 71234.0, -41116.0, 5334.0, -4935.0, -1848.0,
351 66.0, 434.0, -1748.0, 3780.0, -701.0 },
352 { -47645.0, 11647.0, 2166.0, 3194.0, 679.0,
353 0.0, -244.0, -419.0, -2531.0, 48.0 }
354 };
355
356/*
357** -----
358** Pluto
359** -----
360*/
361
362/*
363** Coefficients for fundamental arguments: mean longitudes (degrees)
364** and mean rate of change of longitude (degrees per Julian century)
365** for Jupiter, Saturn and Pluto
366*/
367 static double dj0 = 34.35, djd = 3034.9057,
368 ds0 = 50.08, dsd = 1222.1138,
369 dp0 = 238.96, dpd = 144.9600;
370
371/* Coefficients for latitude, longitude, radius vector */
372 static double dl0 = 238.956785, dld0 = 144.96,
373 db0 = -3.908202,
374 dr0 = 40.7247248;
375
376/*
377** Coefficients for periodic terms (Meeus's Table 36.A)
378*/
379 struct ab {
380 double a; /* sine component */
381 double b; /* cosine component */
382 };
383 struct tm {
384 int ij; /* Jupiter contribution to argument */
385 int is; /* Saturn contribution to argument */
386 int ip; /* Pluto contribution to argument */
387 struct ab dlbr[3]; /* longitude (degrees),
388 latitude (degrees),
389 radius vector (AU) */
390 };
391 static struct tm term[] = {
392
393 /* 1 */ { 0, 0, 1, { { -19798886e-6, 19848454e-6 },
394 { -5453098e-6, -14974876e-6 },
395 { 66867334e-7, 68955876e-7 } } },
396 /* 2 */ { 0, 0, 2, { { 897499e-6, -4955707e-6 },
397 { 3527363e-6, 1672673e-6 },
398 { -11826086e-7, -333765e-7 } } },
399 /* 3 */ { 0, 0, 3, { { 610820e-6, 1210521e-6 },
400 { -1050939e-6, 327763e-6 },
401 { 1593657e-7, -1439953e-7 } } },
402 /* 4 */ { 0, 0, 4, { { -341639e-6, -189719e-6 },
403 { 178691e-6, -291925e-6 },
404 { -18948e-7, 482443e-7 } } },
405 /* 5 */ { 0, 0, 5, { { 129027e-6, -34863e-6 },
406 { 18763e-6, 100448e-6 },
407 { -66634e-7, -85576e-7 } } },
408 /* 6 */ { 0, 0, 6, { { -38215e-6, 31061e-6 },
409 { -30594e-6, -25838e-6 },
410 { 30841e-7, -5765e-7 } } },
411 /* 7 */ { 0, 1, -1, { { 20349e-6, -9886e-6 },
412 { 4965e-6, 11263e-6 },
413 { -6140e-7, 22254e-7 } } },
414 /* 8 */ { 0, 1, 0, { { -4045e-6, -4904e-6 },
415 { 310e-6, -132e-6 },
416 { 4434e-7, 4443e-7 } } },
417 /* 9 */ { 0, 1, 1, { { -5885e-6, -3238e-6 },
418 { 2036e-6, -947e-6 },
419 { -1518e-7, 641e-7 } } },
420 /* 10 */ { 0, 1, 2, { { -3812e-6, 3011e-6 },
421 { -2e-6, -674e-6 },
422 { -5e-7, 792e-7 } } },
423 /* 11 */ { 0, 1, 3, { { -601e-6, 3468e-6 },
424 { -329e-6, -563e-6 },
425 { 518e-7, 518e-7 } } },
426 /* 12 */ { 0, 2, -2, { { 1237e-6, 463e-6 },
427 { -64e-6, 39e-6 },
428 { -13e-7, -221e-7 } } },
429 /* 13 */ { 0, 2, -1, { { 1086e-6, -911e-6 },
430 { -94e-6, 210e-6 },
431 { 837e-7, -494e-7 } } },
432 /* 14 */ { 0, 2, 0, { { 595e-6, -1229e-6 },
433 { -8e-6, -160e-6 },
434 { -281e-7, 616e-7 } } },
435 /* 15 */ { 1, -1, 0, { { 2484e-6, -485e-6 },
436 { -177e-6, 259e-6 },
437 { 260e-7, -395e-7 } } },
438 /* 16 */ { 1, -1, 1, { { 839e-6, -1414e-6 },
439 { 17e-6, 234e-6 },
440 { -191e-7, -396e-7 } } },
441 /* 17 */ { 1, 0, -3, { { -964e-6, 1059e-6 },
442 { 582e-6, -285e-6 },
443 { -3218e-7, 370e-7 } } },
444 /* 18 */ { 1, 0, -2, { { -2303e-6, -1038e-6 },
445 { -298e-6, 692e-6 },
446 { 8019e-7, -7869e-7 } } },
447 /* 19 */ { 1, 0, -1, { { 7049e-6, 747e-6 },
448 { 157e-6, 201e-6 },
449 { 105e-7, 45637e-7 } } },
450 /* 20 */ { 1, 0, 0, { { 1179e-6, -358e-6 },
451 { 304e-6, 825e-6 },
452 { 8623e-7, 8444e-7 } } },
453 /* 21 */ { 1, 0, 1, { { 393e-6, -63e-6 },
454 { -124e-6, -29e-6 },
455 { -896e-7, -801e-7 } } },
456 /* 22 */ { 1, 0, 2, { { 111e-6, -268e-6 },
457 { 15e-6, 8e-6 },
458 { 208e-7, -122e-7 } } },
459 /* 23 */ { 1, 0, 3, { { -52e-6, -154e-6 },
460 { 7e-6, 15e-6 },
461 { -133e-7, 65e-7 } } },
462 /* 24 */ { 1, 0, 4, { { -78e-6, -30e-6 },
463 { 2e-6, 2e-6 },
464 { -16e-7, 1e-7 } } },
465 /* 25 */ { 1, 1, -3, { { -34e-6, -26e-6 },
466 { 4e-6, 2e-6 },
467 { -22e-7, 7e-7 } } },
468 /* 26 */ { 1, 1, -2, { { -43e-6, 1e-6 },
469 { 3e-6, 0e-6 },
470 { -8e-7, 16e-7 } } },
471 /* 27 */ { 1, 1, -1, { { -15e-6, 21e-6 },
472 { 1e-6, -1e-6 },
473 { 2e-7, 9e-7 } } },
474 /* 28 */ { 1, 1, 0, { { -1e-6, 15e-6 },
475 { 0e-6, -2e-6 },
476 { 12e-7, 5e-7 } } },
477 /* 29 */ { 1, 1, 1, { { 4e-6, 7e-6 },
478 { 1e-6, 0e-6 },
479 { 1e-7, -3e-7 } } },
480 /* 30 */ { 1, 1, 3, { { 1e-6, 5e-6 },
481 { 1e-6, -1e-6 },
482 { 1e-7, 0e-7 } } },
483 /* 31 */ { 2, 0, -6, { { 8e-6, 3e-6 },
484 { -2e-6, -3e-6 },
485 { 9e-7, 5e-7 } } },
486 /* 32 */ { 2, 0, -5, { { -3e-6, 6e-6 },
487 { 1e-6, 2e-6 },
488 { 2e-7, -1e-7 } } },
489 /* 33 */ { 2, 0, -4, { { 6e-6, -13e-6 },
490 { -8e-6, 2e-6 },
491 { 14e-7, 10e-7 } } },
492 /* 34 */ { 2, 0, -3, { { 10e-6, 22e-6 },
493 { 10e-6, -7e-6 },
494 { -65e-7, 12e-7 } } },
495 /* 35 */ { 2, 0, -2, { { -57e-6, -32e-6 },
496 { 0e-6, 21e-6 },
497 { 126e-7, -233e-7 } } },
498 /* 36 */ { 2, 0, -1, { { 157e-6, -46e-6 },
499 { 8e-6, 5e-6 },
500 { 270e-7, 1068e-7 } } },
501 /* 37 */ { 2, 0, 0, { { 12e-6, -18e-6 },
502 { 13e-6, 16e-6 },
503 { 254e-7, 155e-7 } } },
504 /* 38 */ { 2, 0, 1, { { -4e-6, 8e-6 },
505 { -2e-6, -3e-6 },
506 { -26e-7, -2e-7 } } },
507 /* 39 */ { 2, 0, 2, { { -5e-6, 0e-6 },
508 { 0e-6, 0e-6 },
509 { 7e-7, 0e-7 } } },
510 /* 40 */ { 2, 0, 3, { { 3e-6, 4e-6 },
511 { 0e-6, 1e-6 },
512 { -11e-7, 4e-7 } } },
513 /* 41 */ { 3, 0, -2, { { -1e-6, -1e-6 },
514 { 0e-6, 1e-6 },
515 { 4e-7, -14e-7 } } },
516 /* 42 */ { 3, 0, -1, { { 6e-6, -3e-6 },
517 { 0e-6, 0e-6 },
518 { 18e-7, 35e-7 } } },
519 /* 43 */ { 3, 0, 0, { { -1e-6, -2e-6 },
520 { 0e-6, 1e-6 },
521 { 13e-7, 3e-7 } } } };
522
523
524
525/* Validate the planet number. */
526 if ( np < 1 || np > 9 ) {
527 *jstat = -1;
528 for ( i = 0; i <= 5; i++ ) pv[i] = 0.0;
529 return;
530 } else {
531 ip = np - 1;
532 }
533
534/* Separate algorithms for Pluto and the rest. */
535 if ( np != 9 ) {
536
537 /* ----------------------- */
538 /* Mercury through Neptune */
539 /* ----------------------- */
540
541 /* Time: Julian millennia since J2000. */
542 t = ( date - 51544.5 ) / 365250.0;
543
544 /* OK status unless remote epoch. */
545 *jstat = ( fabs ( t ) <= 1.0 ) ? 0 : 1;
546
547 /* Compute the mean elements. */
548 da = a[ip][0] + ( a[ip][1] + a[ip][2] * t ) * t;
549 dl = ( 3600.0 * dlm[ip][0] + ( dlm[ip][1] + dlm[ip][2] * t ) * t )
550 * DAS2R;
551 de = e[ip][0] + ( e[ip][1] + e[ip][2] * t ) * t;
552 dpe = dmod ( ( 3600.0 * pi[ip][0] + ( pi[ip][1] + pi[ip][2] * t ) * t )
553 * DAS2R,D2PI );
554 di = ( 3600.0 * dinc[ip][0] + ( dinc[ip][1] + dinc[ip][2] * t ) * t )
555 * DAS2R;
556 dom = dmod( ( 3600.0 * omega[ip][0] + ( omega[ip][1]
557 + omega[ip][2] * t ) * t ) * DAS2R, D2PI );
558
559 /* Apply the trigonometric terms. */
560 dmu = 0.35953620 * t;
561 for ( j = 0; j <= 7; j++ ) {
562 arga = dkp[ip][j] * dmu;
563 argl = dkq[ip][j] * dmu;
564 da += ( ca[ip][j] * cos ( arga ) +
565 sa[ip][j] * sin ( arga ) ) * 1e-7;
566 dl += ( clo[ip][j] * cos ( argl ) +
567 slo[ip][j] * sin ( argl ) ) * 1e-7;
568 }
569 arga = dkp[ip][8] * dmu;
570 da += t * ( ca[ip][8] * cos ( arga ) +
571 sa[ip][8] * sin ( arga ) ) * 1e-7;
572 for ( j = 8; j <= 9; j++ ) {
573 argl = dkq[ip][j] * dmu;
574 dl += t * ( clo[ip][j] * cos ( argl ) +
575 slo[ip][j] * sin ( argl ) ) * 1e-7;
576 }
577 dl = dmod ( dl, D2PI );
578
579 /* Daily motion. */
580 dm = GCON * sqrt ( ( 1.0 + 1.0 / amas[ip] ) / ( da * da * da ) );
581
582 /* Make the prediction. */
583 slaPlanel ( date, 1, date, di, dom, dpe, da, de, dl, dm, pv, &j );
584 if ( j < 0 ) *jstat = -2;
585
586
587 } else {
588
589 /* ----- */
590 /* Pluto */
591 /* ----- */
592
593 /* Time: Julian centuries since J2000. */
594 t = ( date - 51544.5 ) / 36525.0;
595
596 /* OK status unless remote epoch. */
597 *jstat = t >= -1.15 && t <= 1.0 ? 0 : -1;
598
599 /* Fundamental arguments (radians). */
600 dj = ( dj0 + djd * t ) * DD2R;
601 ds = ( ds0 + dsd * t ) * DD2R;
602 dp = ( dp0 + dpd * t ) * DD2R;
603
604 /* Initialize coefficients and derivatives. */
605 for ( i = 0; i < 3; i++ ) {
606 wlbr[i] = 0.0;
607 wlbrd[i] = 0.0;
608 }
609
610 /* Term by term through Meeus Table 36.A. */
611 for ( j = 0; j < ( sizeof term / sizeof term[0] ); j++ ) {
612
613 /* Argument and derivative (radians, radians per century). */
614 wj = (double) ( term[j].ij );
615 ws = (double) ( term[j].is );
616 wp = (double) ( term[j].ip );
617 al = wj * dj + ws * ds + wp * dp;
618 ald = ( wj * djd + ws * dsd + wp * dpd ) * DD2R;
619
620 /* Functions of argument. */
621 sal = sin ( al );
622 cal = cos ( al );
623
624 /* Periodic terms in longitude, latitude, radius vector. */
625 for ( i = 0; i < 3; i++ ) {
626
627 /* A and B coefficients (deg, AU). */
628 ac = term[j].dlbr[i].a;
629 bc = term[j].dlbr[i].b;
630
631 /* Periodic terms (deg, AU, deg/Jc, AU/Jc). */
632 wlbr[i] = wlbr[i] + ac * sal + bc * cal;
633 wlbrd[i] = wlbrd[i] + ( ac * cal - bc * sal ) * ald;
634 }
635 }
636
637 /* Heliocentric longitude and derivative (radians, radians/sec). */
638 dl = ( dl0 + dld0 * t + wlbr[0] ) * DD2R;
639 dld = ( dld0 + wlbrd[0] ) * DD2R / SPC;
640
641 /* Heliocentric latitude and derivative (radians, radians/sec). */
642 db = ( db0 + wlbr[1] ) * DD2R;
643 dbd = wlbrd[1] * DD2R / SPC;
644
645 /* Heliocentric radius vector and derivative (AU, AU/sec). */
646 dr = dr0 + wlbr[2];
647 drd = wlbrd[2] / SPC;
648
649 /* Functions of latitude, longitude, radius vector. */
650 sl = sin ( dl );
651 cl = cos ( dl );
652 sb = sin ( db );
653 cb = cos ( db );
654 slcb = sl * cb;
655 clcb = cl * cb;
656
657 /* Heliocentric vector and derivative, J2000 ecliptic and equinox. */
658 x = dr * clcb;
659 y = dr * slcb;
660 z = dr * sb;
661 xd = drd * clcb - dr * ( cl * sb * dbd + slcb * dld );
662 yd = drd * slcb + dr * ( - sl * sb * dbd + clcb * dld );
663 zd = drd * sb + dr * cb * dbd;
664
665 /* Transform to J2000 equator and equinox. */
666 pv[0] = x;
667 pv[1] = y * CE - z * SE;
668 pv[2] = y * SE + z * CE;
669 pv[3] = xd;
670 pv[4] = yd * CE - zd * SE;
671 pv[5] = yd * SE + zd * CE;
672 }
673}
Note: See TracBrowser for help on using the repository browser.