source: trunk/FACT++/sofa/src/plan94.c@ 18355

Last change on this file since 18355 was 18346, checked in by tbretz, 11 years ago
File size: 22.6 KB
Line 
1#include "sofa.h"
2
3int iauPlan94(double date1, double date2, int np, double pv[2][3])
4/*
5** - - - - - - - - - -
6** i a u P l a n 9 4
7** - - - - - - - - - -
8**
9** This function is part of the International Astronomical Union's
10** SOFA (Standards Of Fundamental Astronomy) software collection.
11**
12** Status: support function.
13**
14** Approximate heliocentric position and velocity of a nominated major
15** planet: Mercury, Venus, EMB, Mars, Jupiter, Saturn, Uranus or
16** Neptune (but not the Earth itself).
17**
18** Given:
19** date1 double TDB date part A (Note 1)
20** date2 double TDB date part B (Note 1)
21** np int planet (1=Mercury, 2=Venus, 3=EMB, 4=Mars,
22** 5=Jupiter, 6=Saturn, 7=Uranus, 8=Neptune)
23**
24** Returned (argument):
25** pv double[2][3] planet p,v (heliocentric, J2000.0, AU,AU/d)
26**
27** Returned (function value):
28** int status: -1 = illegal NP (outside 1-8)
29** 0 = OK
30** +1 = warning: year outside 1000-3000
31** +2 = warning: failed to converge
32**
33** Notes:
34**
35** 1) The date date1+date2 is in the TDB time scale (in practice TT can
36** be used) and is a Julian Date, apportioned in any convenient way
37** between the two arguments. For example, JD(TDB)=2450123.7 could
38** be expressed in any of these ways, among others:
39**
40** date1 date2
41**
42** 2450123.7 0.0 (JD method)
43** 2451545.0 -1421.3 (J2000 method)
44** 2400000.5 50123.2 (MJD method)
45** 2450123.5 0.2 (date & time method)
46**
47** The JD method is the most natural and convenient to use in cases
48** where the loss of several decimal digits of resolution is
49** acceptable. The J2000 method is best matched to the way the
50** argument is handled internally and will deliver the optimum
51** resolution. The MJD method and the date & time methods are both
52** good compromises between resolution and convenience. The limited
53** accuracy of the present algorithm is such that any of the methods
54** is satisfactory.
55**
56** 2) If an np value outside the range 1-8 is supplied, an error status
57** (function value -1) is returned and the pv vector set to zeroes.
58**
59** 3) For np=3 the result is for the Earth-Moon Barycenter. To obtain
60** the heliocentric position and velocity of the Earth, use instead
61** the SOFA function iauEpv00.
62**
63** 4) On successful return, the array pv contains the following:
64**
65** pv[0][0] x }
66** pv[0][1] y } heliocentric position, AU
67** pv[0][2] z }
68**
69** pv[1][0] xdot }
70** pv[1][1] ydot } heliocentric velocity, AU/d
71** pv[1][2] zdot }
72**
73** The reference frame is equatorial and is with respect to the
74** mean equator and equinox of epoch J2000.0.
75**
76** 5) The algorithm is due to J.L. Simon, P. Bretagnon, J. Chapront,
77** M. Chapront-Touze, G. Francou and J. Laskar (Bureau des
78** Longitudes, Paris, France). From comparisons with JPL
79** ephemeris DE102, they quote the following maximum errors
80** over the interval 1800-2050:
81**
82** L (arcsec) B (arcsec) R (km)
83**
84** Mercury 4 1 300
85** Venus 5 1 800
86** EMB 6 1 1000
87** Mars 17 1 7700
88** Jupiter 71 5 76000
89** Saturn 81 13 267000
90** Uranus 86 7 712000
91** Neptune 11 1 253000
92**
93** Over the interval 1000-3000, they report that the accuracy is no
94** worse than 1.5 times that over 1800-2050. Outside 1000-3000 the
95** accuracy declines.
96**
97** Comparisons of the present function with the JPL DE200 ephemeris
98** give the following RMS errors over the interval 1960-2025:
99**
100** position (km) velocity (m/s)
101**
102** Mercury 334 0.437
103** Venus 1060 0.855
104** EMB 2010 0.815
105** Mars 7690 1.98
106** Jupiter 71700 7.70
107** Saturn 199000 19.4
108** Uranus 564000 16.4
109** Neptune 158000 14.4
110**
111** Comparisons against DE200 over the interval 1800-2100 gave the
112** following maximum absolute differences. (The results using
113** DE406 were essentially the same.)
114**
115** L (arcsec) B (arcsec) R (km) Rdot (m/s)
116**
117** Mercury 7 1 500 0.7
118** Venus 7 1 1100 0.9
119** EMB 9 1 1300 1.0
120** Mars 26 1 9000 2.5
121** Jupiter 78 6 82000 8.2
122** Saturn 87 14 263000 24.6
123** Uranus 86 7 661000 27.4
124** Neptune 11 2 248000 21.4
125**
126** 6) The present SOFA re-implementation of the original Simon et al.
127** Fortran code differs from the original in the following respects:
128**
129** * C instead of Fortran.
130**
131** * The date is supplied in two parts.
132**
133** * The result is returned only in equatorial Cartesian form;
134** the ecliptic longitude, latitude and radius vector are not
135** returned.
136**
137** * The result is in the J2000.0 equatorial frame, not ecliptic.
138**
139** * More is done in-line: there are fewer calls to subroutines.
140**
141** * Different error/warning status values are used.
142**
143** * A different Kepler's-equation-solver is used (avoiding
144** use of double precision complex).
145**
146** * Polynomials in t are nested to minimize rounding errors.
147**
148** * Explicit double constants are used to avoid mixed-mode
149** expressions.
150**
151** None of the above changes affects the result significantly.
152**
153** 7) The returned status indicates the most serious condition
154** encountered during execution of the function. Illegal np is
155** considered the most serious, overriding failure to converge,
156** which in turn takes precedence over the remote date warning.
157**
158** Called:
159** iauAnp normalize angle into range 0 to 2pi
160**
161** Reference: Simon, J.L, Bretagnon, P., Chapront, J.,
162** Chapront-Touze, M., Francou, G., and Laskar, J.,
163** Astron. Astrophys. 282, 663 (1994).
164**
165** This revision: 2013 June 18
166**
167** SOFA release 2015-02-09
168**
169** Copyright (C) 2015 IAU SOFA Board. See notes at end.
170*/
171{
172/* Gaussian constant */
173 static const double GK = 0.017202098950;
174
175/* Sin and cos of J2000.0 mean obliquity (IAU 1976) */
176 static const double SINEPS = 0.3977771559319137;
177 static const double COSEPS = 0.9174820620691818;
178
179/* Maximum number of iterations allowed to solve Kepler's equation */
180 static const int KMAX = 10;
181
182 int jstat, i, k;
183 double t, da, dl, de, dp, di, dom, dmu, arga, argl, am,
184 ae, dae, ae2, at, r, v, si2, xq, xp, tl, xsw,
185 xcw, xm2, xf, ci2, xms, xmc, xpxq2, x, y, z;
186
187/* Planetary inverse masses */
188 static const double amas[] = { 6023600.0, /* Mercury */
189 408523.5, /* Venus */
190 328900.5, /* EMB */
191 3098710.0, /* Mars */
192 1047.355, /* Jupiter */
193 3498.5, /* Saturn */
194 22869.0, /* Uranus */
195 19314.0 }; /* Neptune */
196
197/*
198** Tables giving the mean Keplerian elements, limited to t^2 terms:
199**
200** a semi-major axis (AU)
201** dlm mean longitude (degree and arcsecond)
202** e eccentricity
203** pi longitude of the perihelion (degree and arcsecond)
204** dinc inclination (degree and arcsecond)
205** omega longitude of the ascending node (degree and arcsecond)
206*/
207
208 static const double a[][3] = {
209 { 0.3870983098, 0.0, 0.0 }, /* Mercury */
210 { 0.7233298200, 0.0, 0.0 }, /* Venus */
211 { 1.0000010178, 0.0, 0.0 }, /* EMB */
212 { 1.5236793419, 3e-10, 0.0 }, /* Mars */
213 { 5.2026032092, 19132e-10, -39e-10 }, /* Jupiter */
214 { 9.5549091915, -0.0000213896, 444e-10 }, /* Saturn */
215 { 19.2184460618, -3716e-10, 979e-10 }, /* Uranus */
216 { 30.1103868694, -16635e-10, 686e-10 } /* Neptune */
217 };
218
219 static const double dlm[][3] = {
220 { 252.25090552, 5381016286.88982, -1.92789 },
221 { 181.97980085, 2106641364.33548, 0.59381 },
222 { 100.46645683, 1295977422.83429, -2.04411 },
223 { 355.43299958, 689050774.93988, 0.94264 },
224 { 34.35151874, 109256603.77991, -30.60378 },
225 { 50.07744430, 43996098.55732, 75.61614 },
226 { 314.05500511, 15424811.93933, -1.75083 },
227 { 304.34866548, 7865503.20744, 0.21103 }
228 };
229
230 static const double e[][3] = {
231 { 0.2056317526, 0.0002040653, -28349e-10 },
232 { 0.0067719164, -0.0004776521, 98127e-10 },
233 { 0.0167086342, -0.0004203654, -0.0000126734 },
234 { 0.0934006477, 0.0009048438, -80641e-10 },
235 { 0.0484979255, 0.0016322542, -0.0000471366 },
236 { 0.0555481426, -0.0034664062, -0.0000643639 },
237 { 0.0463812221, -0.0002729293, 0.0000078913 },
238 { 0.0094557470, 0.0000603263, 0.0 }
239 };
240
241 static const double pi[][3] = {
242 { 77.45611904, 5719.11590, -4.83016 },
243 { 131.56370300, 175.48640, -498.48184 },
244 { 102.93734808, 11612.35290, 53.27577 },
245 { 336.06023395, 15980.45908, -62.32800 },
246 { 14.33120687, 7758.75163, 259.95938 },
247 { 93.05723748, 20395.49439, 190.25952 },
248 { 173.00529106, 3215.56238, -34.09288 },
249 { 48.12027554, 1050.71912, 27.39717 }
250 };
251
252 static const double dinc[][3] = {
253 { 7.00498625, -214.25629, 0.28977 },
254 { 3.39466189, -30.84437, -11.67836 },
255 { 0.0, 469.97289, -3.35053 },
256 { 1.84972648, -293.31722, -8.11830 },
257 { 1.30326698, -71.55890, 11.95297 },
258 { 2.48887878, 91.85195, -17.66225 },
259 { 0.77319689, -60.72723, 1.25759 },
260 { 1.76995259, 8.12333, 0.08135 }
261 };
262
263 static const double omega[][3] = {
264 { 48.33089304, -4515.21727, -31.79892 },
265 { 76.67992019, -10008.48154, -51.32614 },
266 { 174.87317577, -8679.27034, 15.34191 },
267 { 49.55809321, -10620.90088, -230.57416 },
268 { 100.46440702, 6362.03561, 326.52178 },
269 { 113.66550252, -9240.19942, -66.23743 },
270 { 74.00595701, 2669.15033, 145.93964 },
271 { 131.78405702, -221.94322, -0.78728 }
272 };
273
274/* Tables for trigonometric terms to be added to the mean elements of */
275/* the semi-major axes */
276
277 static const double kp[][9] = {
278 { 69613, 75645, 88306, 59899, 15746, 71087, 142173, 3086, 0 },
279 { 21863, 32794, 26934, 10931, 26250, 43725, 53867, 28939, 0 },
280 { 16002, 21863, 32004, 10931, 14529, 16368, 15318, 32794, 0 },
281 { 6345, 7818, 15636, 7077, 8184, 14163, 1107, 4872, 0 },
282 { 1760, 1454, 1167, 880, 287, 2640, 19, 2047, 1454 },
283 { 574, 0, 880, 287, 19, 1760, 1167, 306, 574 },
284 { 204, 0, 177, 1265, 4, 385, 200, 208, 204 },
285 { 0, 102, 106, 4, 98, 1367, 487, 204, 0 }
286 };
287
288 static const double ca[][9] = {
289 { 4, -13, 11, -9, -9, -3, -1, 4, 0 },
290 { -156, 59, -42, 6, 19, -20, -10, -12, 0 },
291 { 64, -152, 62, -8, 32, -41, 19, -11, 0 },
292 { 124, 621, -145, 208, 54, -57, 30, 15, 0 },
293 { -23437, -2634, 6601, 6259, -1507,-1821, 2620, -2115, -1489 },
294 { 62911,-119919, 79336,17814,-24241,12068, 8306, -4893, 8902 },
295 { 389061,-262125,-44088, 8387,-22976,-2093, -615, -9720, 6633 },
296 { -412235,-157046,-31430,37817, -9740, -13, -7449, 9644, 0 }
297 };
298
299 static const double sa[][9] = {
300 { -29, -1, 9, 6, -6, 5, 4, 0, 0 },
301 { -48, -125, -26, -37, 18, -13, -20, -2, 0 },
302 { -150, -46, 68, 54, 14, 24, -28, 22, 0 },
303 { -621, 532, -694, -20, 192, -94, 71, -73, 0 },
304 { -14614,-19828, -5869, 1881, -4372, -2255, 782, 930, 913 },
305 { 139737, 0, 24667, 51123, -5102, 7429, -4095, -1976, -9566 },
306 { -138081, 0, 37205,-49039,-41901,-33872,-27037,-12474, 18797 },
307 { 0, 28492,133236, 69654, 52322,-49577,-26430, -3593, 0 }
308 };
309
310/* Tables giving the trigonometric terms to be added to the mean */
311/* elements of the mean longitudes */
312
313 static const double kq[][10] = {
314 { 3086,15746,69613,59899,75645,88306, 12661, 2658, 0, 0 },
315 { 21863,32794,10931, 73, 4387,26934, 1473, 2157, 0, 0 },
316 { 10,16002,21863,10931, 1473,32004, 4387, 73, 0, 0 },
317 { 10, 6345, 7818, 1107,15636, 7077, 8184, 532, 10, 0 },
318 { 19, 1760, 1454, 287, 1167, 880, 574, 2640, 19, 1454 },
319 { 19, 574, 287, 306, 1760, 12, 31, 38, 19, 574 },
320 { 4, 204, 177, 8, 31, 200, 1265, 102, 4, 204 },
321 { 4, 102, 106, 8, 98, 1367, 487, 204, 4, 102 }
322 };
323
324 static const double cl[][10] = {
325 { 21, -95, -157, 41, -5, 42, 23, 30, 0, 0 },
326 { -160, -313, -235, 60, -74, -76, -27, 34, 0, 0 },
327 { -325, -322, -79, 232, -52, 97, 55, -41, 0, 0 },
328 { 2268, -979, 802, 602, -668, -33, 345, 201, -55, 0 },
329 { 7610, -4997,-7689,-5841,-2617, 1115,-748,-607, 6074, 354 },
330 { -18549, 30125,20012, -730, 824, 23,1289,-352, -14767, -2062 },
331 { -135245,-14594, 4197,-4030,-5630,-2898,2540,-306, 2939, 1986 },
332 { 89948, 2103, 8963, 2695, 3682, 1648, 866,-154, -1963, -283 }
333 };
334
335 static const double sl[][10] = {
336 { -342, 136, -23, 62, 66, -52, -33, 17, 0, 0 },
337 { 524, -149, -35, 117, 151, 122, -71, -62, 0, 0 },
338 { -105, -137, 258, 35, -116, -88,-112, -80, 0, 0 },
339 { 854, -205, -936, -240, 140, -341, -97, -232, 536, 0 },
340 { -56980, 8016, 1012, 1448,-3024,-3710, 318, 503, 3767, 577 },
341 { 138606,-13478,-4964, 1441,-1319,-1482, 427, 1236, -9167, -1918 },
342 { 71234,-41116, 5334,-4935,-1848, 66, 434, -1748, 3780, -701 },
343 { -47645, 11647, 2166, 3194, 679, 0,-244, -419, -2531, 48 }
344 };
345
346/*--------------------------------------------------------------------*/
347
348/* Validate the planet number. */
349 if ((np < 1) || (np > 8)) {
350 jstat = -1;
351
352 /* Reset the result in case of failure. */
353 for (k = 0; k < 2; k++) {
354 for (i = 0; i < 3; i++) {
355 pv[k][i] = 0.0;
356 }
357 }
358
359 } else {
360
361 /* Decrement the planet number to start at zero. */
362 np--;
363
364 /* Time: Julian millennia since J2000.0. */
365 t = ((date1 - DJ00) + date2) / DJM;
366
367 /* OK status unless remote date. */
368 jstat = fabs(t) <= 1.0 ? 0 : 1;
369
370 /* Compute the mean elements. */
371 da = a[np][0] +
372 (a[np][1] +
373 a[np][2] * t) * t;
374 dl = (3600.0 * dlm[np][0] +
375 (dlm[np][1] +
376 dlm[np][2] * t) * t) * DAS2R;
377 de = e[np][0] +
378 ( e[np][1] +
379 e[np][2] * t) * t;
380 dp = iauAnpm((3600.0 * pi[np][0] +
381 (pi[np][1] +
382 pi[np][2] * t) * t) * DAS2R);
383 di = (3600.0 * dinc[np][0] +
384 (dinc[np][1] +
385 dinc[np][2] * t) * t) * DAS2R;
386 dom = iauAnpm((3600.0 * omega[np][0] +
387 (omega[np][1] +
388 omega[np][2] * t) * t) * DAS2R);
389
390 /* Apply the trigonometric terms. */
391 dmu = 0.35953620 * t;
392 for (k = 0; k < 8; k++) {
393 arga = kp[np][k] * dmu;
394 argl = kq[np][k] * dmu;
395 da += (ca[np][k] * cos(arga) +
396 sa[np][k] * sin(arga)) * 1e-7;
397 dl += (cl[np][k] * cos(argl) +
398 sl[np][k] * sin(argl)) * 1e-7;
399 }
400 arga = kp[np][8] * dmu;
401 da += t * (ca[np][8] * cos(arga) +
402 sa[np][8] * sin(arga)) * 1e-7;
403 for (k = 8; k < 10; k++) {
404 argl = kq[np][k] * dmu;
405 dl += t * (cl[np][k] * cos(argl) +
406 sl[np][k] * sin(argl)) * 1e-7;
407 }
408 dl = fmod(dl, D2PI);
409
410 /* Iterative soln. of Kepler's equation to get eccentric anomaly. */
411 am = dl - dp;
412 ae = am + de * sin(am);
413 k = 0;
414 dae = 1.0;
415 while (k < KMAX && fabs(dae) > 1e-12) {
416 dae = (am - ae + de * sin(ae)) / (1.0 - de * cos(ae));
417 ae += dae;
418 k++;
419 if (k == KMAX-1) jstat = 2;
420 }
421
422 /* True anomaly. */
423 ae2 = ae / 2.0;
424 at = 2.0 * atan2(sqrt((1.0 + de) / (1.0 - de)) * sin(ae2),
425 cos(ae2));
426
427 /* Distance (AU) and speed (radians per day). */
428 r = da * (1.0 - de * cos(ae));
429 v = GK * sqrt((1.0 + 1.0 / amas[np]) / (da * da * da));
430
431 si2 = sin(di / 2.0);
432 xq = si2 * cos(dom);
433 xp = si2 * sin(dom);
434 tl = at + dp;
435 xsw = sin(tl);
436 xcw = cos(tl);
437 xm2 = 2.0 * (xp * xcw - xq * xsw);
438 xf = da / sqrt(1 - de * de);
439 ci2 = cos(di / 2.0);
440 xms = (de * sin(dp) + xsw) * xf;
441 xmc = (de * cos(dp) + xcw) * xf;
442 xpxq2 = 2 * xp * xq;
443
444 /* Position (J2000.0 ecliptic x,y,z in AU). */
445 x = r * (xcw - xm2 * xp);
446 y = r * (xsw + xm2 * xq);
447 z = r * (-xm2 * ci2);
448
449 /* Rotate to equatorial. */
450 pv[0][0] = x;
451 pv[0][1] = y * COSEPS - z * SINEPS;
452 pv[0][2] = y * SINEPS + z * COSEPS;
453
454 /* Velocity (J2000.0 ecliptic xdot,ydot,zdot in AU/d). */
455 x = v * (( -1.0 + 2.0 * xp * xp) * xms + xpxq2 * xmc);
456 y = v * (( 1.0 - 2.0 * xq * xq) * xmc - xpxq2 * xms);
457 z = v * (2.0 * ci2 * (xp * xms + xq * xmc));
458
459 /* Rotate to equatorial. */
460 pv[1][0] = x;
461 pv[1][1] = y * COSEPS - z * SINEPS;
462 pv[1][2] = y * SINEPS + z * COSEPS;
463
464 }
465
466/* Return the status. */
467 return jstat;
468
469/*----------------------------------------------------------------------
470**
471** Copyright (C) 2015
472** Standards Of Fundamental Astronomy Board
473** of the International Astronomical Union.
474**
475** =====================
476** SOFA Software License
477** =====================
478**
479** NOTICE TO USER:
480**
481** BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
482** CONDITIONS WHICH APPLY TO ITS USE.
483**
484** 1. The Software is owned by the IAU SOFA Board ("SOFA").
485**
486** 2. Permission is granted to anyone to use the SOFA software for any
487** purpose, including commercial applications, free of charge and
488** without payment of royalties, subject to the conditions and
489** restrictions listed below.
490**
491** 3. You (the user) may copy and distribute SOFA source code to others,
492** and use and adapt its code and algorithms in your own software,
493** on a world-wide, royalty-free basis. That portion of your
494** distribution that does not consist of intact and unchanged copies
495** of SOFA source code files is a "derived work" that must comply
496** with the following requirements:
497**
498** a) Your work shall be marked or carry a statement that it
499** (i) uses routines and computations derived by you from
500** software provided by SOFA under license to you; and
501** (ii) does not itself constitute software provided by and/or
502** endorsed by SOFA.
503**
504** b) The source code of your derived work must contain descriptions
505** of how the derived work is based upon, contains and/or differs
506** from the original SOFA software.
507**
508** c) The names of all routines in your derived work shall not
509** include the prefix "iau" or "sofa" or trivial modifications
510** thereof such as changes of case.
511**
512** d) The origin of the SOFA components of your derived work must
513** not be misrepresented; you must not claim that you wrote the
514** original software, nor file a patent application for SOFA
515** software or algorithms embedded in the SOFA software.
516**
517** e) These requirements must be reproduced intact in any source
518** distribution and shall apply to anyone to whom you have
519** granted a further right to modify the source code of your
520** derived work.
521**
522** Note that, as originally distributed, the SOFA software is
523** intended to be a definitive implementation of the IAU standards,
524** and consequently third-party modifications are discouraged. All
525** variations, no matter how minor, must be explicitly marked as
526** such, as explained above.
527**
528** 4. You shall not cause the SOFA software to be brought into
529** disrepute, either by misuse, or use for inappropriate tasks, or
530** by inappropriate modification.
531**
532** 5. The SOFA software is provided "as is" and SOFA makes no warranty
533** as to its use or performance. SOFA does not and cannot warrant
534** the performance or results which the user may obtain by using the
535** SOFA software. SOFA makes no warranties, express or implied, as
536** to non-infringement of third party rights, merchantability, or
537** fitness for any particular purpose. In no event will SOFA be
538** liable to the user for any consequential, incidental, or special
539** damages, including any lost profits or lost savings, even if a
540** SOFA representative has been advised of such damages, or for any
541** claim by any third party.
542**
543** 6. The provision of any version of the SOFA software under the terms
544** and conditions specified herein does not imply that future
545** versions will also be made available under the same terms and
546** conditions.
547*
548** In any published work or commercial product which uses the SOFA
549** software directly, acknowledgement (see www.iausofa.org) is
550** appreciated.
551**
552** Correspondence concerning SOFA software should be addressed as
553** follows:
554**
555** By email: sofa@ukho.gov.uk
556** By post: IAU SOFA Center
557** HM Nautical Almanac Office
558** UK Hydrographic Office
559** Admiralty Way, Taunton
560** Somerset, TA1 2DN
561** United Kingdom
562**
563**--------------------------------------------------------------------*/
564}
Note: See TracBrowser for help on using the repository browser.