source: Sophya/trunk/SophyaExt/XephemAstroLib/vsop87.c@ 2937

Last change on this file since 2937 was 2818, checked in by cmv, 20 years ago

Update de Xephem 3.7 cmv 21/08/2005

File size: 6.2 KB
Line 
1/* VSOP87 planetary theory
2 *
3 * currently uses version VSOP87D:
4 * heliocentric spherical, mean ecliptic of date.
5 *
6 * calculation of rates (daily changes) is optional;
7 * see header file for the necessary #define's
8 *
9 * rough orientation on calculation time, miliseconds
10 * on an HP 715/75, all planets Mercury to Neptune, prec=0.0:
11 *
12 * terms with rates without rates
13 * 3598 11 7.1
14 * 31577 51 44
15 *
16 * with secular terms for JD 2232395.0 19/12/1399 0h TDB:
17 *
18 * FULL PRECISION code (31577 terms), milliseconds
19 * prec terms rates no rates
20 * 1e-8 15086 62 36
21 * 1e-7 10105 44 25
22 * 1e-6 3725 20 13
23 * 1e-5 1324 11 7.8
24 * 1e-4 443 7.0 6.0
25 * 1e-3 139 6.0 5.0
26 *
27 * REDUCED PRECISION code (3598 terms), milliseconds
28 * prec terms rates no rates
29 * 1e-7 2463 9.9 5.5
30 * 1e-6 1939 8.0 4.5
31 * 1e-5 1131 4.9 2.9
32 * 1e-4 443 2.2 1.5
33 * 1e-3 139 1.0 0.9
34 */
35
36#include <math.h>
37
38#include "astro.h"
39#include "vsop87.h"
40
41#define VSOP_A1000 365250.0 /* days per millenium */
42#define VSOP_MAXALPHA 5 /* max degree of time */
43
44/******************************************************************
45 * adapted from BdL FORTRAN Code; stern
46 *
47 * Reference : Bureau des Longitudes - PBGF9502
48 *
49 * Object : calculate a VSOP87 position for a given time.
50 *
51 * Input :
52 *
53 * mj modified julian date, counted from J1900.0
54 * time scale : dynamical time TDB.
55 *
56 * obj object number as in astro.h, NB: not for pluto
57 *
58 * prec relative precision
59 *
60 * if prec is equal to 0 then the precision is the precision
61 * p0 of the complete solution VSOP87.
62 * Mercury p0 = 0.6 10**-8
63 * Venus p0 = 2.5 10**-8
64 * Earth p0 = 2.5 10**-8
65 * Mars p0 = 10.0 10**-8
66 * Jupiter p0 = 35.0 10**-8
67 * Saturn p0 = 70.0 10**-8
68 * Uranus p0 = 8.0 10**-8
69 * Neptune p0 = 42.0 10**-8
70 *
71 * if prec is not equal to 0, let us say in between p0 and
72 * 10**-3, the precision is :
73 * for the positions :
74 * - prec*a0 au for the distances.
75 * - prec rad for the other variables.
76 * for the velocities :
77 * - prec*a0 au/day for the distances.
78 * - prec rad/day for the other variables.
79 * a0 is the semi-major axis of the body.
80 *
81 * Output :
82 *
83 * ret[6] array of the results (double).
84 *
85 * for spherical coordinates :
86 * 1: longitude (rd)
87 * 2: latitude (rd)
88 * 3: radius (au)
89 * #if VSOP_GETRATE:
90 * 4: longitude velocity (rad/day)
91 * 5: latitude velocity (rad/day)
92 * 6: radius velocity (au/day)
93 *
94 * return: error index (int)
95 * 0: no error.
96 * 2: object out of range [MERCURY .. NEPTUNE, SUN]
97 * 3: precision out of range [0.0 .. 1e-3]
98 ******************************************************************/
99int
100vsop87 (double mj, int obj, double prec, double *ret)
101{
102 static double (*vx_map[])[3] = { /* data tables */
103 vx_mercury, vx_venus, vx_mars, vx_jupiter,
104 vx_saturn, vx_uranus, vx_neptune, 0, vx_earth,
105 };
106 static int (*vn_map[])[3] = { /* indexes */
107 vn_mercury, vn_venus, vn_mars, vn_jupiter,
108 vn_saturn, vn_uranus, vn_neptune, 0, vn_earth,
109 };
110 static double a0[] = { /* semimajor axes; for precision ctrl only */
111 0.39, 0.72, 1.5, 5.2, 9.6, 19.2, 30.1, 39.5, 1.0,
112 };
113 double (*vx_obj)[3] = vx_map[obj]; /* VSOP87 data and indexes */
114 int (*vn_obj)[3] = vn_map[obj];
115
116 double t[VSOP_MAXALPHA+1]; /* powers of time */
117 double t_abs[VSOP_MAXALPHA+1]; /* powers of abs(time) */
118 double q; /* aux for precision control */
119 int i, cooidx, alpha; /* misc indexes */
120
121 if (obj == PLUTO || obj > SUN)
122 return (2);
123
124 if (prec < 0.0 || prec > 1e-3)
125 return(3);
126
127 /* zero result array */
128 for (i = 0; i < 6; ++i) ret[i] = 0.0;
129
130 /* time and its powers */
131 t[0] = 1.0;
132 t[1] = (mj - J2000)/VSOP_A1000;
133 for (i = 2; i <= VSOP_MAXALPHA; ++i) t[i] = t[i-1] * t[1];
134 t_abs[0] = 1.0;
135 for (i = 1; i <= VSOP_MAXALPHA; ++i) t_abs[i] = fabs(t[i]);
136
137 /* precision control */
138 q = -log10(prec + 1e-35) - 2; /* decades below 1e-2 */
139 q = VSOP_ASCALE * prec / 10.0 / q; /* reduce threshold progressively
140 * for higher precision */
141
142 /* do the term summation; first the spatial dimensions */
143 for (cooidx = 0; cooidx < 3; ++cooidx) {
144
145 /* then the powers of time */
146 for (alpha = 0; vn_obj[alpha+1][cooidx] ; ++alpha) {
147 double p, term, termdot;
148
149 /* precision threshold */
150 p= alpha ? q/(t_abs[alpha] + alpha*t_abs[alpha-1]*1e-4 + 1e-35) : q;
151#if VSOP_SPHERICAL
152 if (cooidx == 2) /* scale by semimajor axis for radius */
153#endif
154 p *= a0[obj];
155
156 term = termdot = 0.0;
157 for (i = vn_obj[alpha][cooidx]; i < vn_obj[alpha+1][cooidx]; ++i) {
158 double a, b, c, arg;
159
160 a = vx_obj[i][0];
161 if (a < p) continue; /* ignore small terms */
162
163 b = vx_obj[i][1];
164 c = vx_obj[i][2];
165
166 arg = b + c * t[1];
167 term += a * cos(arg);
168#if VSOP_GETRATE
169 termdot += -c * a * sin(arg);
170#endif
171 }
172
173 ret[cooidx] += t[alpha] * term;
174#if VSOP_GETRATE
175 ret[cooidx + 3] += t[alpha] * termdot +
176 ((alpha > 0) ? alpha * t[alpha - 1] * term : 0.0);
177#endif
178 } /* alpha */
179 } /* cooidx */
180
181 for (i = 0; i < 6; ++i) ret[i] /= VSOP_ASCALE;
182
183#if VSOP_SPHERICAL
184 /* reduce longitude to 0..2pi */
185 ret[0] -= floor(ret[0]/(2.*PI)) * (2.*PI);
186#endif
187
188#if VSOP_GETRATE
189 /* convert millenium rate to day rate */
190 for (i = 3; i < 6; ++i) ret[i] /= VSOP_A1000;
191#endif
192
193#if VSOP_SPHERICAL
194 /* reduction from dynamical equinox of VSOP87 to FK5;
195 */
196 if (prec < 5e-7) { /* 5e-7 rad = 0.1 arc seconds */
197 double L1, c1, s1;
198 L1 = ret[0] - degrad(13.97 * t[1] - 0.031 * t[2]);
199 c1 = cos(L1); s1 = sin(L1);
200 ret[0] += degrad(-0.09033 + 0.03916 * (c1 + s1) * tan(ret[1]))/3600.0;
201 ret[1] += degrad(0.03916 * (c1 - s1))/3600.0;
202 }
203#endif
204
205 return (0);
206}
207
208/* For RCS Only -- Do Not Edit */
209static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: vsop87.c,v $ $Date: 2005-08-21 10:02:40 $ $Revision: 1.6 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.