precess.c 8.51 KB
Newer Older
Jeroen Vreeken's avatar
Jeroen Vreeken committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
/* Precession of the equinox and ecliptic
 * from epoch Julian date J to or from J2000.0
 *
 * Program by Steve Moshier.  */

#define WILLIAMS 1
/* James G. Williams, "Contributions to the Earth's obliquity rate,
   precession, and nutation,"  Astron. J. 108, 711-724 (1994)  */

#define SIMON 0
/* J. L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touze', G. Francou,
   and J. Laskar, "Numerical Expressions for precession formulae and
   mean elements for the Moon and the planets," Astronomy and Astrophysics
   282, 663-683 (1994)  */

#define IAU 0
/* IAU Coefficients are from:
 * J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
 * "Expressions for the Precession Quantities Based upon the IAU
 * (1976) System of Astronomical Constants,"  Astronomy and
 * Astrophysics 58, 1-16 (1977).
 */

#define LASKAR 0
/* Newer formulas that cover a much longer time span are from:
 * J. Laskar, "Secular terms of classical planetary theories
 * using the results of general theory," Astronomy and Astrophysics
 * 157, 59070 (1986).
 *
 * See also:
 * P. Bretagnon and G. Francou, "Planetary theories in rectangular
 * and spherical variables. VSOP87 solutions," Astronomy and
 * Astrophysics 202, 309-315 (1988).
 *
 * Laskar's expansions are said by Bretagnon and Francou
 * to have "a precision of about 1" over 10000 years before
 * and after J2000.0 in so far as the precession constants p^0_A
 * and epsilon^0_A are perfectly known."
 *
 * Bretagnon and Francou's expansions for the node and inclination
 * of the ecliptic were derived from Laskar's data but were truncated
 * after the term in T**6. I have recomputed these expansions from
 * Laskar's data, retaining powers up to T**10 in the result.
 *
 * The following table indicates the differences between the result
 * of the IAU formula and Laskar's formula using four different test
 * vectors, checking at J2000 plus and minus the indicated number
 * of years.
 *
 *   Years       Arc
 * from J2000  Seconds
 * ----------  -------
 *        0	  0
 *      100	.006	
 *      200     .006
 *      500     .015
 *     1000     .28
 *     2000    6.4
 *     3000   38.
 *    10000 9400.
 */

#define DOUBLE double
#if __STDC__
double cos (double);
double sin (double);
int epsiln (double);
#else
double cos(), sin();
extern int epsiln();
#endif
#define COS cos
#define SIN sin
extern DOUBLE J2000; /* = 2451545.0, 2000 January 1.5 */
extern DOUBLE STR; /* = 4.8481368110953599359e-6 radians per arc second */
extern DOUBLE coseps, sineps; /* see epsiln.c */

/* In WILLIAMS and SIMON, Laskar's terms of order higher than t^4
   have been retained, because Simon et al mention that the solution
   is the same except for the lower order terms.  */
#if WILLIAMS
static DOUBLE pAcof[] = {
#if 1
  /* Corrections to Williams (1994) introduced in DE403.  */
 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
 -0.235316, 0.076, 110.5414, 50287.91959
#else
 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
 -0.235316, 0.076, 110.5407, 50287.70000
#endif
 };
#endif
#if SIMON
/* Precession coefficients from Simon et al: */
static DOUBLE pAcof[] = {
 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
 -0.235316, 0.07732, 111.2022, 50288.200 };
#endif
#if LASKAR
/* Precession coefficients taken from Laskar's paper: */
static DOUBLE pAcof[] = {
 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
 -0.235316, 0.07732, 111.1971, 50290.966 };
#endif
#if WILLIAMS
/* Pi from Williams' 1994 paper, in radians.  No change in DE403.  */
static DOUBLE nodecof[] = {
6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 1.9e-10, 
-3.54e-9, -1.8103e-7,  1.26e-7,  7.436169e-5,
-0.04207794833,  3.052115282424};
/* pi from Williams' 1994 paper, in radians.  No change in DE403.  */
static DOUBLE inclcof[] = {
1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11, 
-5.4000441e-11, 1.32115526e-9, -6.012e-7, -1.62442e-5,
 0.00227850649, 0.0 };
#endif
#if SIMON
static DOUBLE nodecof[] = {
6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 1.9e-10, 
-3.54e-9, -1.8103e-7, 2.579e-8, 7.4379679e-5,
-0.0420782900, 3.0521126906};

static DOUBLE inclcof[] = {
1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11, 
-5.4000441e-11, 1.32115526e-9, -5.99908e-7, -1.624383e-5,
 0.002278492868, 0.0 };
#endif
#if LASKAR
/* Node and inclination of the earth's orbit computed from
 * Laskar's data as done in Bretagnon and Francou's paper.
 * Units are radians.
 */
static DOUBLE nodecof[] = {
6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 6.3190131e-10, 
-3.48388152e-9, -1.813065896e-7, 2.75036225e-8, 7.4394531426e-5,
-0.042078604317, 3.052112654975 };

static DOUBLE inclcof[] = {
1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11, 
-5.4000441e-11, 1.32115526e-9, -5.998737027e-7, -1.6242797091e-5,
 0.002278495537, 0.0 };
#endif


/* Subroutine arguments:
 *
 * R = rectangular equatorial coordinate vector to be precessed.
 *     The result is written back into the input vector.
 * J = Julian date
 * direction =
 *      Precess from J to J2000: direction = 1
 *      Precess from J2000 to J: direction = -1
 * Note that if you want to precess from J1 to J2, you would
 * first go from J1 to J2000, then call the program again
 * to go from J2000 to J2.
 */

int precess( R, J, direction )
DOUBLE R[], J;
int direction;
{
DOUBLE A, B, T, pA, W, z;
DOUBLE x[3];
DOUBLE *p;
int i;
#if IAU
DOUBLE sinth, costh, sinZ, cosZ, sinz, cosz, Z, TH;
#endif

if( J == J2000 )
	return(0);
/* Each precession angle is specified by a polynomial in
 * T = Julian centuries from J2000.0.  See AA page B18.
 */
T = (J - J2000)/36525.0;

#if IAU
/* Use IAU formula only for a few centuries, if at all.  */
if( FABS(T) > Two )
	goto laskar;

Z =  (( 0.017998*T + 0.30188)*T + 2306.2181)*T*STR;
z =  (( 0.018203*T + 1.09468)*T + 2306.2181)*T*STR;
TH = ((-0.041833*T - 0.42665)*T + 2004.3109)*T*STR;

sinth = SIN(TH);
costh = COS(TH);
sinZ = SIN(Z);
cosZ = COS(Z);
sinz = SIN(z);
cosz = COS(z);
A = cosZ*costh;
B = sinZ*costh;

if( direction < 0 )
	{ /* From J2000.0 to J */
	x[0] =    (A*cosz - sinZ*sinz)*R[0]
	        - (B*cosz + cosZ*sinz)*R[1]
	                  - sinth*cosz*R[2];

	x[1] =    (A*sinz + sinZ*cosz)*R[0]
	        - (B*sinz - cosZ*cosz)*R[1]
	                  - sinth*sinz*R[2];

	x[2] =              cosZ*sinth*R[0]
	                  - sinZ*sinth*R[1]
	                       + costh*R[2];
	}
else
	{ /* From J to J2000.0 */
	x[0] =    (A*cosz - sinZ*sinz)*R[0]
	        + (A*sinz + sinZ*cosz)*R[1]
	                  + cosZ*sinth*R[2];

	x[1] =   -(B*cosz + cosZ*sinz)*R[0]
	        - (B*sinz - cosZ*cosz)*R[1]
	                  - sinZ*sinth*R[2];

	x[2] =             -sinth*cosz*R[0]
	                  - sinth*sinz*R[1]
	                       + costh*R[2];
	}	
goto done;

laskar:
#endif /* IAU */

/* Implementation by elementary rotations using Laskar's expansions.
 * First rotate about the x axis from the initial equator
 * to the ecliptic. (The input is equatorial.)
 */
if( direction == 1 )
	epsiln( J ); /* To J2000 */
else
	epsiln( J2000 ); /* From J2000 */
x[0] = R[0];
z = coseps*R[1] + sineps*R[2];
x[2] = -sineps*R[1] + coseps*R[2];
x[1] = z;

/* Precession in longitude
 */
T /= 10.0; /* thousands of years */
p = pAcof;
pA = *p++;
for( i=0; i<9; i++ )
	pA = pA * T + *p++;
pA *= STR * T;

/* Node of the moving ecliptic on the J2000 ecliptic.
 */
p = nodecof;
W = *p++;
for( i=0; i<10; i++ )
	W = W * T + *p++;

/* Rotate about z axis to the node.
 */
if( direction == 1 )
	z = W + pA;
else
	z = W;
B = COS(z);
A = SIN(z);
z = B * x[0] + A * x[1];
x[1] = -A * x[0] + B * x[1];
x[0] = z;

/* Rotate about new x axis by the inclination of the moving
 * ecliptic on the J2000 ecliptic.
 */
p = inclcof;
z = *p++;
for( i=0; i<10; i++ )
	z = z * T + *p++;
if( direction == 1 )
	z = -z;
B = COS(z);
A = SIN(z);
z = B * x[1] + A * x[2];
x[2] = -A * x[1] + B * x[2];
x[1] = z;

/* Rotate about new z axis back from the node.
 */
if( direction == 1 )
	z = -W;
else
	z = -W - pA;
B = COS(z);
A = SIN(z);
z = B * x[0] + A * x[1];
x[1] = -A * x[0] + B * x[1];
x[0] = z;

/* Rotate about x axis to final equator.
 */
if( direction == 1 )
	epsiln( J2000 );
else
	epsiln( J );
z = coseps * x[1] - sineps * x[2];
x[2] = sineps * x[1] + coseps * x[2];
x[1] = z;

#if IAU
done:
#endif

for( i=0; i<3; i++ )
	R[i] = x[i];
return(0);
}