source: trunk/FACT++/erfa/src/plan94.c@ 18403

Last change on this file since 18403 was 18348, checked in by tbretz, 9 years ago
File size: 21.2 KB
Line 
1#include "erfa.h"
2
3int eraPlan94(double date1, double date2, int np, double pv[2][3])
4/*
5** - - - - - - - - - -
6** e r a P l a n 9 4
7** - - - - - - - - - -
8**
9** Approximate heliocentric position and velocity of a nominated major
10** planet: Mercury, Venus, EMB, Mars, Jupiter, Saturn, Uranus or
11** Neptune (but not the Earth itself).
12**
13** Given:
14** date1 double TDB date part A (Note 1)
15** date2 double TDB date part B (Note 1)
16** np int planet (1=Mercury, 2=Venus, 3=EMB, 4=Mars,
17** 5=Jupiter, 6=Saturn, 7=Uranus, 8=Neptune)
18**
19** Returned (argument):
20** pv double[2][3] planet p,v (heliocentric, J2000.0, AU,AU/d)
21**
22** Returned (function value):
23** int status: -1 = illegal NP (outside 1-8)
24** 0 = OK
25** +1 = warning: year outside 1000-3000
26** +2 = warning: failed to converge
27**
28** Notes:
29**
30** 1) The date date1+date2 is in the TDB time scale (in practice TT can
31** be used) and is a Julian Date, apportioned in any convenient way
32** between the two arguments. For example, JD(TDB)=2450123.7 could
33** be expressed in any of these ways, among others:
34**
35** date1 date2
36**
37** 2450123.7 0.0 (JD method)
38** 2451545.0 -1421.3 (J2000 method)
39** 2400000.5 50123.2 (MJD method)
40** 2450123.5 0.2 (date & time method)
41**
42** The JD method is the most natural and convenient to use in cases
43** where the loss of several decimal digits of resolution is
44** acceptable. The J2000 method is best matched to the way the
45** argument is handled internally and will deliver the optimum
46** resolution. The MJD method and the date & time methods are both
47** good compromises between resolution and convenience. The limited
48** accuracy of the present algorithm is such that any of the methods
49** is satisfactory.
50**
51** 2) If an np value outside the range 1-8 is supplied, an error status
52** (function value -1) is returned and the pv vector set to zeroes.
53**
54** 3) For np=3 the result is for the Earth-Moon Barycenter. To obtain
55** the heliocentric position and velocity of the Earth, use instead
56** the ERFA function eraEpv00.
57**
58** 4) On successful return, the array pv contains the following:
59**
60** pv[0][0] x }
61** pv[0][1] y } heliocentric position, AU
62** pv[0][2] z }
63**
64** pv[1][0] xdot }
65** pv[1][1] ydot } heliocentric velocity, AU/d
66** pv[1][2] zdot }
67**
68** The reference frame is equatorial and is with respect to the
69** mean equator and equinox of epoch J2000.0.
70**
71** 5) The algorithm is due to J.L. Simon, P. Bretagnon, J. Chapront,
72** M. Chapront-Touze, G. Francou and J. Laskar (Bureau des
73** Longitudes, Paris, France). From comparisons with JPL
74** ephemeris DE102, they quote the following maximum errors
75** over the interval 1800-2050:
76**
77** L (arcsec) B (arcsec) R (km)
78**
79** Mercury 4 1 300
80** Venus 5 1 800
81** EMB 6 1 1000
82** Mars 17 1 7700
83** Jupiter 71 5 76000
84** Saturn 81 13 267000
85** Uranus 86 7 712000
86** Neptune 11 1 253000
87**
88** Over the interval 1000-3000, they report that the accuracy is no
89** worse than 1.5 times that over 1800-2050. Outside 1000-3000 the
90** accuracy declines.
91**
92** Comparisons of the present function with the JPL DE200 ephemeris
93** give the following RMS errors over the interval 1960-2025:
94**
95** position (km) velocity (m/s)
96**
97** Mercury 334 0.437
98** Venus 1060 0.855
99** EMB 2010 0.815
100** Mars 7690 1.98
101** Jupiter 71700 7.70
102** Saturn 199000 19.4
103** Uranus 564000 16.4
104** Neptune 158000 14.4
105**
106** Comparisons against DE200 over the interval 1800-2100 gave the
107** following maximum absolute differences. (The results using
108** DE406 were essentially the same.)
109**
110** L (arcsec) B (arcsec) R (km) Rdot (m/s)
111**
112** Mercury 7 1 500 0.7
113** Venus 7 1 1100 0.9
114** EMB 9 1 1300 1.0
115** Mars 26 1 9000 2.5
116** Jupiter 78 6 82000 8.2
117** Saturn 87 14 263000 24.6
118** Uranus 86 7 661000 27.4
119** Neptune 11 2 248000 21.4
120**
121** 6) The present ERFA re-implementation of the original Simon et al.
122** Fortran code differs from the original in the following respects:
123**
124** * C instead of Fortran.
125**
126** * The date is supplied in two parts.
127**
128** * The result is returned only in equatorial Cartesian form;
129** the ecliptic longitude, latitude and radius vector are not
130** returned.
131**
132** * The result is in the J2000.0 equatorial frame, not ecliptic.
133**
134** * More is done in-line: there are fewer calls to subroutines.
135**
136** * Different error/warning status values are used.
137**
138** * A different Kepler's-equation-solver is used (avoiding
139** use of double precision complex).
140**
141** * Polynomials in t are nested to minimize rounding errors.
142**
143** * Explicit double constants are used to avoid mixed-mode
144** expressions.
145**
146** None of the above changes affects the result significantly.
147**
148** 7) The returned status indicates the most serious condition
149** encountered during execution of the function. Illegal np is
150** considered the most serious, overriding failure to converge,
151** which in turn takes precedence over the remote date warning.
152**
153** Called:
154** eraAnp normalize angle into range 0 to 2pi
155**
156** Reference: Simon, J.L, Bretagnon, P., Chapront, J.,
157** Chapront-Touze, M., Francou, G., and Laskar, J.,
158** Astron. Astrophys. 282, 663 (1994).
159**
160** Copyright (C) 2013-2015, NumFOCUS Foundation.
161** Derived, with permission, from the SOFA library. See notes at end of file.
162*/
163{
164/* Gaussian constant */
165 static const double GK = 0.017202098950;
166
167/* Sin and cos of J2000.0 mean obliquity (IAU 1976) */
168 static const double SINEPS = 0.3977771559319137;
169 static const double COSEPS = 0.9174820620691818;
170
171/* Maximum number of iterations allowed to solve Kepler's equation */
172 static const int KMAX = 10;
173
174 int jstat, i, k;
175 double t, da, dl, de, dp, di, dom, dmu, arga, argl, am,
176 ae, dae, ae2, at, r, v, si2, xq, xp, tl, xsw,
177 xcw, xm2, xf, ci2, xms, xmc, xpxq2, x, y, z;
178
179/* Planetary inverse masses */
180 static const double amas[] = { 6023600.0, /* Mercury */
181 408523.5, /* Venus */
182 328900.5, /* EMB */
183 3098710.0, /* Mars */
184 1047.355, /* Jupiter */
185 3498.5, /* Saturn */
186 22869.0, /* Uranus */
187 19314.0 }; /* Neptune */
188
189/*
190** Tables giving the mean Keplerian elements, limited to t^2 terms:
191**
192** a semi-major axis (AU)
193** dlm mean longitude (degree and arcsecond)
194** e eccentricity
195** pi longitude of the perihelion (degree and arcsecond)
196** dinc inclination (degree and arcsecond)
197** omega longitude of the ascending node (degree and arcsecond)
198*/
199
200 static const double a[][3] = {
201 { 0.3870983098, 0.0, 0.0 }, /* Mercury */
202 { 0.7233298200, 0.0, 0.0 }, /* Venus */
203 { 1.0000010178, 0.0, 0.0 }, /* EMB */
204 { 1.5236793419, 3e-10, 0.0 }, /* Mars */
205 { 5.2026032092, 19132e-10, -39e-10 }, /* Jupiter */
206 { 9.5549091915, -0.0000213896, 444e-10 }, /* Saturn */
207 { 19.2184460618, -3716e-10, 979e-10 }, /* Uranus */
208 { 30.1103868694, -16635e-10, 686e-10 } /* Neptune */
209 };
210
211 static const double dlm[][3] = {
212 { 252.25090552, 5381016286.88982, -1.92789 },
213 { 181.97980085, 2106641364.33548, 0.59381 },
214 { 100.46645683, 1295977422.83429, -2.04411 },
215 { 355.43299958, 689050774.93988, 0.94264 },
216 { 34.35151874, 109256603.77991, -30.60378 },
217 { 50.07744430, 43996098.55732, 75.61614 },
218 { 314.05500511, 15424811.93933, -1.75083 },
219 { 304.34866548, 7865503.20744, 0.21103 }
220 };
221
222 static const double e[][3] = {
223 { 0.2056317526, 0.0002040653, -28349e-10 },
224 { 0.0067719164, -0.0004776521, 98127e-10 },
225 { 0.0167086342, -0.0004203654, -0.0000126734 },
226 { 0.0934006477, 0.0009048438, -80641e-10 },
227 { 0.0484979255, 0.0016322542, -0.0000471366 },
228 { 0.0555481426, -0.0034664062, -0.0000643639 },
229 { 0.0463812221, -0.0002729293, 0.0000078913 },
230 { 0.0094557470, 0.0000603263, 0.0 }
231 };
232
233 static const double pi[][3] = {
234 { 77.45611904, 5719.11590, -4.83016 },
235 { 131.56370300, 175.48640, -498.48184 },
236 { 102.93734808, 11612.35290, 53.27577 },
237 { 336.06023395, 15980.45908, -62.32800 },
238 { 14.33120687, 7758.75163, 259.95938 },
239 { 93.05723748, 20395.49439, 190.25952 },
240 { 173.00529106, 3215.56238, -34.09288 },
241 { 48.12027554, 1050.71912, 27.39717 }
242 };
243
244 static const double dinc[][3] = {
245 { 7.00498625, -214.25629, 0.28977 },
246 { 3.39466189, -30.84437, -11.67836 },
247 { 0.0, 469.97289, -3.35053 },
248 { 1.84972648, -293.31722, -8.11830 },
249 { 1.30326698, -71.55890, 11.95297 },
250 { 2.48887878, 91.85195, -17.66225 },
251 { 0.77319689, -60.72723, 1.25759 },
252 { 1.76995259, 8.12333, 0.08135 }
253 };
254
255 static const double omega[][3] = {
256 { 48.33089304, -4515.21727, -31.79892 },
257 { 76.67992019, -10008.48154, -51.32614 },
258 { 174.87317577, -8679.27034, 15.34191 },
259 { 49.55809321, -10620.90088, -230.57416 },
260 { 100.46440702, 6362.03561, 326.52178 },
261 { 113.66550252, -9240.19942, -66.23743 },
262 { 74.00595701, 2669.15033, 145.93964 },
263 { 131.78405702, -221.94322, -0.78728 }
264 };
265
266/* Tables for trigonometric terms to be added to the mean elements of */
267/* the semi-major axes */
268
269 static const double kp[][9] = {
270 { 69613, 75645, 88306, 59899, 15746, 71087, 142173, 3086, 0 },
271 { 21863, 32794, 26934, 10931, 26250, 43725, 53867, 28939, 0 },
272 { 16002, 21863, 32004, 10931, 14529, 16368, 15318, 32794, 0 },
273 { 6345, 7818, 15636, 7077, 8184, 14163, 1107, 4872, 0 },
274 { 1760, 1454, 1167, 880, 287, 2640, 19, 2047, 1454 },
275 { 574, 0, 880, 287, 19, 1760, 1167, 306, 574 },
276 { 204, 0, 177, 1265, 4, 385, 200, 208, 204 },
277 { 0, 102, 106, 4, 98, 1367, 487, 204, 0 }
278 };
279
280 static const double ca[][9] = {
281 { 4, -13, 11, -9, -9, -3, -1, 4, 0 },
282 { -156, 59, -42, 6, 19, -20, -10, -12, 0 },
283 { 64, -152, 62, -8, 32, -41, 19, -11, 0 },
284 { 124, 621, -145, 208, 54, -57, 30, 15, 0 },
285 { -23437, -2634, 6601, 6259, -1507,-1821, 2620, -2115, -1489 },
286 { 62911,-119919, 79336,17814,-24241,12068, 8306, -4893, 8902 },
287 { 389061,-262125,-44088, 8387,-22976,-2093, -615, -9720, 6633 },
288 { -412235,-157046,-31430,37817, -9740, -13, -7449, 9644, 0 }
289 };
290
291 static const double sa[][9] = {
292 { -29, -1, 9, 6, -6, 5, 4, 0, 0 },
293 { -48, -125, -26, -37, 18, -13, -20, -2, 0 },
294 { -150, -46, 68, 54, 14, 24, -28, 22, 0 },
295 { -621, 532, -694, -20, 192, -94, 71, -73, 0 },
296 { -14614,-19828, -5869, 1881, -4372, -2255, 782, 930, 913 },
297 { 139737, 0, 24667, 51123, -5102, 7429, -4095, -1976, -9566 },
298 { -138081, 0, 37205,-49039,-41901,-33872,-27037,-12474, 18797 },
299 { 0, 28492,133236, 69654, 52322,-49577,-26430, -3593, 0 }
300 };
301
302/* Tables giving the trigonometric terms to be added to the mean */
303/* elements of the mean longitudes */
304
305 static const double kq[][10] = {
306 { 3086,15746,69613,59899,75645,88306, 12661, 2658, 0, 0 },
307 { 21863,32794,10931, 73, 4387,26934, 1473, 2157, 0, 0 },
308 { 10,16002,21863,10931, 1473,32004, 4387, 73, 0, 0 },
309 { 10, 6345, 7818, 1107,15636, 7077, 8184, 532, 10, 0 },
310 { 19, 1760, 1454, 287, 1167, 880, 574, 2640, 19, 1454 },
311 { 19, 574, 287, 306, 1760, 12, 31, 38, 19, 574 },
312 { 4, 204, 177, 8, 31, 200, 1265, 102, 4, 204 },
313 { 4, 102, 106, 8, 98, 1367, 487, 204, 4, 102 }
314 };
315
316 static const double cl[][10] = {
317 { 21, -95, -157, 41, -5, 42, 23, 30, 0, 0 },
318 { -160, -313, -235, 60, -74, -76, -27, 34, 0, 0 },
319 { -325, -322, -79, 232, -52, 97, 55, -41, 0, 0 },
320 { 2268, -979, 802, 602, -668, -33, 345, 201, -55, 0 },
321 { 7610, -4997,-7689,-5841,-2617, 1115,-748,-607, 6074, 354 },
322 { -18549, 30125,20012, -730, 824, 23,1289,-352, -14767, -2062 },
323 { -135245,-14594, 4197,-4030,-5630,-2898,2540,-306, 2939, 1986 },
324 { 89948, 2103, 8963, 2695, 3682, 1648, 866,-154, -1963, -283 }
325 };
326
327 static const double sl[][10] = {
328 { -342, 136, -23, 62, 66, -52, -33, 17, 0, 0 },
329 { 524, -149, -35, 117, 151, 122, -71, -62, 0, 0 },
330 { -105, -137, 258, 35, -116, -88,-112, -80, 0, 0 },
331 { 854, -205, -936, -240, 140, -341, -97, -232, 536, 0 },
332 { -56980, 8016, 1012, 1448,-3024,-3710, 318, 503, 3767, 577 },
333 { 138606,-13478,-4964, 1441,-1319,-1482, 427, 1236, -9167, -1918 },
334 { 71234,-41116, 5334,-4935,-1848, 66, 434, -1748, 3780, -701 },
335 { -47645, 11647, 2166, 3194, 679, 0,-244, -419, -2531, 48 }
336 };
337
338/*--------------------------------------------------------------------*/
339
340/* Validate the planet number. */
341 if ((np < 1) || (np > 8)) {
342 jstat = -1;
343
344 /* Reset the result in case of failure. */
345 for (k = 0; k < 2; k++) {
346 for (i = 0; i < 3; i++) {
347 pv[k][i] = 0.0;
348 }
349 }
350
351 } else {
352
353 /* Decrement the planet number to start at zero. */
354 np--;
355
356 /* Time: Julian millennia since J2000.0. */
357 t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJM;
358
359 /* OK status unless remote date. */
360 jstat = fabs(t) <= 1.0 ? 0 : 1;
361
362 /* Compute the mean elements. */
363 da = a[np][0] +
364 (a[np][1] +
365 a[np][2] * t) * t;
366 dl = (3600.0 * dlm[np][0] +
367 (dlm[np][1] +
368 dlm[np][2] * t) * t) * ERFA_DAS2R;
369 de = e[np][0] +
370 ( e[np][1] +
371 e[np][2] * t) * t;
372 dp = eraAnpm((3600.0 * pi[np][0] +
373 (pi[np][1] +
374 pi[np][2] * t) * t) * ERFA_DAS2R);
375 di = (3600.0 * dinc[np][0] +
376 (dinc[np][1] +
377 dinc[np][2] * t) * t) * ERFA_DAS2R;
378 dom = eraAnpm((3600.0 * omega[np][0] +
379 (omega[np][1] +
380 omega[np][2] * t) * t) * ERFA_DAS2R);
381
382 /* Apply the trigonometric terms. */
383 dmu = 0.35953620 * t;
384 for (k = 0; k < 8; k++) {
385 arga = kp[np][k] * dmu;
386 argl = kq[np][k] * dmu;
387 da += (ca[np][k] * cos(arga) +
388 sa[np][k] * sin(arga)) * 1e-7;
389 dl += (cl[np][k] * cos(argl) +
390 sl[np][k] * sin(argl)) * 1e-7;
391 }
392 arga = kp[np][8] * dmu;
393 da += t * (ca[np][8] * cos(arga) +
394 sa[np][8] * sin(arga)) * 1e-7;
395 for (k = 8; k < 10; k++) {
396 argl = kq[np][k] * dmu;
397 dl += t * (cl[np][k] * cos(argl) +
398 sl[np][k] * sin(argl)) * 1e-7;
399 }
400 dl = fmod(dl, ERFA_D2PI);
401
402 /* Iterative soln. of Kepler's equation to get eccentric anomaly. */
403 am = dl - dp;
404 ae = am + de * sin(am);
405 k = 0;
406 dae = 1.0;
407 while (k < KMAX && fabs(dae) > 1e-12) {
408 dae = (am - ae + de * sin(ae)) / (1.0 - de * cos(ae));
409 ae += dae;
410 k++;
411 if (k == KMAX-1) jstat = 2;
412 }
413
414 /* True anomaly. */
415 ae2 = ae / 2.0;
416 at = 2.0 * atan2(sqrt((1.0 + de) / (1.0 - de)) * sin(ae2),
417 cos(ae2));
418
419 /* Distance (AU) and speed (radians per day). */
420 r = da * (1.0 - de * cos(ae));
421 v = GK * sqrt((1.0 + 1.0 / amas[np]) / (da * da * da));
422
423 si2 = sin(di / 2.0);
424 xq = si2 * cos(dom);
425 xp = si2 * sin(dom);
426 tl = at + dp;
427 xsw = sin(tl);
428 xcw = cos(tl);
429 xm2 = 2.0 * (xp * xcw - xq * xsw);
430 xf = da / sqrt(1 - de * de);
431 ci2 = cos(di / 2.0);
432 xms = (de * sin(dp) + xsw) * xf;
433 xmc = (de * cos(dp) + xcw) * xf;
434 xpxq2 = 2 * xp * xq;
435
436 /* Position (J2000.0 ecliptic x,y,z in AU). */
437 x = r * (xcw - xm2 * xp);
438 y = r * (xsw + xm2 * xq);
439 z = r * (-xm2 * ci2);
440
441 /* Rotate to equatorial. */
442 pv[0][0] = x;
443 pv[0][1] = y * COSEPS - z * SINEPS;
444 pv[0][2] = y * SINEPS + z * COSEPS;
445
446 /* Velocity (J2000.0 ecliptic xdot,ydot,zdot in AU/d). */
447 x = v * (( -1.0 + 2.0 * xp * xp) * xms + xpxq2 * xmc);
448 y = v * (( 1.0 - 2.0 * xq * xq) * xmc - xpxq2 * xms);
449 z = v * (2.0 * ci2 * (xp * xms + xq * xmc));
450
451 /* Rotate to equatorial. */
452 pv[1][0] = x;
453 pv[1][1] = y * COSEPS - z * SINEPS;
454 pv[1][2] = y * SINEPS + z * COSEPS;
455
456 }
457
458/* Return the status. */
459 return jstat;
460
461}
462/*----------------------------------------------------------------------
463**
464**
465** Copyright (C) 2013-2015, NumFOCUS Foundation.
466** All rights reserved.
467**
468** This library is derived, with permission, from the International
469** Astronomical Union's "Standards of Fundamental Astronomy" library,
470** available from http://www.iausofa.org.
471**
472** The ERFA version is intended to retain identical functionality to
473** the SOFA library, but made distinct through different function and
474** file names, as set out in the SOFA license conditions. The SOFA
475** original has a role as a reference standard for the IAU and IERS,
476** and consequently redistribution is permitted only in its unaltered
477** state. The ERFA version is not subject to this restriction and
478** therefore can be included in distributions which do not support the
479** concept of "read only" software.
480**
481** Although the intent is to replicate the SOFA API (other than
482** replacement of prefix names) and results (with the exception of
483** bugs; any that are discovered will be fixed), SOFA is not
484** responsible for any errors found in this version of the library.
485**
486** If you wish to acknowledge the SOFA heritage, please acknowledge
487** that you are using a library derived from SOFA, rather than SOFA
488** itself.
489**
490**
491** TERMS AND CONDITIONS
492**
493** Redistribution and use in source and binary forms, with or without
494** modification, are permitted provided that the following conditions
495** are met:
496**
497** 1 Redistributions of source code must retain the above copyright
498** notice, this list of conditions and the following disclaimer.
499**
500** 2 Redistributions in binary form must reproduce the above copyright
501** notice, this list of conditions and the following disclaimer in
502** the documentation and/or other materials provided with the
503** distribution.
504**
505** 3 Neither the name of the Standards Of Fundamental Astronomy Board,
506** the International Astronomical Union nor the names of its
507** contributors may be used to endorse or promote products derived
508** from this software without specific prior written permission.
509**
510** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
511** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
512** LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
513** FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
514** COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
515** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
516** BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
517** LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
518** CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
519** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
520** ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
521** POSSIBILITY OF SUCH DAMAGE.
522**
523*/
Note: See TracBrowser for help on using the repository browser.