1 | /* nutation (in IAU (1980) expression) and abberation; stern
|
---|
2 | * on an HP PA processor, this reproduces the Almanac nutation values
|
---|
3 | * (given to 0.001") EXACTLY over 750 days (1995 and 1996)
|
---|
4 | */
|
---|
5 | #include <stdio.h>
|
---|
6 | #include <math.h>
|
---|
7 |
|
---|
8 | #include "astro.h"
|
---|
9 |
|
---|
10 | #define NUT_SCALE 1e4
|
---|
11 | #define NUT_SERIES 106
|
---|
12 | #define NUT_MAXMUL 4
|
---|
13 | #define SECPERCIRC (3600.*360.)
|
---|
14 |
|
---|
15 | /* Delaunay arguments, in arc seconds; they differ slightly from ELP82B */
|
---|
16 | static double delaunay[5][4] = {
|
---|
17 | {485866.733, 1717915922.633, 31.310, 0.064}, /* M', moon mean anom */
|
---|
18 | {1287099.804, 129596581.224, -0.577, -0.012}, /* M, sun mean anom */
|
---|
19 | {335778.877, 1739527263.137, -13.257, 0.011}, /* F, moon arg lat */
|
---|
20 | {1072261.307, 1602961601.328, -6.891, 0.019}, /* D, elong moon sun */
|
---|
21 | {450160.280, -6962890.539, 7.455, 0.008}, /* Om, moon l asc node */
|
---|
22 | };
|
---|
23 |
|
---|
24 | /* multipliers for Delaunay arguments */
|
---|
25 | static short multarg[NUT_SERIES][5] = {
|
---|
26 | /* bounds: -2..3, -2..2, -2/0/2/4, -4..4, 0..2 */
|
---|
27 | {0, 0, 0, 0, 1},
|
---|
28 | {0, 0, 0, 0, 2},
|
---|
29 | {-2, 0, 2, 0, 1},
|
---|
30 | {2, 0, -2, 0, 0},
|
---|
31 | {-2, 0, 2, 0, 2},
|
---|
32 | {1, -1, 0, -1, 0},
|
---|
33 | {0, -2, 2, -2, 1},
|
---|
34 | {2, 0, -2, 0, 1},
|
---|
35 | {0, 0, 2, -2, 2},
|
---|
36 | {0, 1, 0, 0, 0},
|
---|
37 | {0, 1, 2, -2, 2},
|
---|
38 | {0, -1, 2, -2, 2},
|
---|
39 | {0, 0, 2, -2, 1},
|
---|
40 | {2, 0, 0, -2, 0},
|
---|
41 | {0, 0, 2, -2, 0},
|
---|
42 | {0, 2, 0, 0, 0},
|
---|
43 | {0, 1, 0, 0, 1},
|
---|
44 | {0, 2, 2, -2, 2},
|
---|
45 | {0, -1, 0, 0, 1},
|
---|
46 | {-2, 0, 0, 2, 1},
|
---|
47 | {0, -1, 2, -2, 1},
|
---|
48 | {2, 0, 0, -2, 1},
|
---|
49 | {0, 1, 2, -2, 1},
|
---|
50 | {1, 0, 0, -1, 0},
|
---|
51 | {2, 1, 0, -2, 0},
|
---|
52 | {0, 0, -2, 2, 1},
|
---|
53 | {0, 1, -2, 2, 0},
|
---|
54 | {0, 1, 0, 0, 2},
|
---|
55 | {-1, 0, 0, 1, 1},
|
---|
56 | {0, 1, 2, -2, 0},
|
---|
57 | {0, 0, 2, 0, 2},
|
---|
58 | {1, 0, 0, 0, 0},
|
---|
59 | {0, 0, 2, 0, 1},
|
---|
60 | {1, 0, 2, 0, 2},
|
---|
61 | {1, 0, 0, -2, 0},
|
---|
62 | {-1, 0, 2, 0, 2},
|
---|
63 | {0, 0, 0, 2, 0},
|
---|
64 | {1, 0, 0, 0, 1},
|
---|
65 | {-1, 0, 0, 0, 1},
|
---|
66 | {-1, 0, 2, 2, 2},
|
---|
67 | {1, 0, 2, 0, 1},
|
---|
68 | {0, 0, 2, 2, 2},
|
---|
69 | {2, 0, 0, 0, 0},
|
---|
70 | {1, 0, 2, -2, 2},
|
---|
71 | {2, 0, 2, 0, 2},
|
---|
72 | {0, 0, 2, 0, 0},
|
---|
73 | {-1, 0, 2, 0, 1},
|
---|
74 | {-1, 0, 0, 2, 1},
|
---|
75 | {1, 0, 0, -2, 1},
|
---|
76 | {-1, 0, 2, 2, 1},
|
---|
77 | {1, 1, 0, -2, 0},
|
---|
78 | {0, 1, 2, 0, 2},
|
---|
79 | {0, -1, 2, 0, 2},
|
---|
80 | {1, 0, 2, 2, 2},
|
---|
81 | {1, 0, 0, 2, 0},
|
---|
82 | {2, 0, 2, -2, 2},
|
---|
83 | {0, 0, 0, 2, 1},
|
---|
84 | {0, 0, 2, 2, 1},
|
---|
85 | {1, 0, 2, -2, 1},
|
---|
86 | {0, 0, 0, -2, 1},
|
---|
87 | {1, -1, 0, 0, 0},
|
---|
88 | {2, 0, 2, 0, 1},
|
---|
89 | {0, 1, 0, -2, 0},
|
---|
90 | {1, 0, -2, 0, 0},
|
---|
91 | {0, 0, 0, 1, 0},
|
---|
92 | {1, 1, 0, 0, 0},
|
---|
93 | {1, 0, 2, 0, 0},
|
---|
94 | {1, -1, 2, 0, 2},
|
---|
95 | {-1, -1, 2, 2, 2},
|
---|
96 | {-2, 0, 0, 0, 1},
|
---|
97 | {3, 0, 2, 0, 2},
|
---|
98 | {0, -1, 2, 2, 2},
|
---|
99 | {1, 1, 2, 0, 2},
|
---|
100 | {-1, 0, 2, -2, 1},
|
---|
101 | {2, 0, 0, 0, 1},
|
---|
102 | {1, 0, 0, 0, 2},
|
---|
103 | {3, 0, 0, 0, 0},
|
---|
104 | {0, 0, 2, 1, 2},
|
---|
105 | {-1, 0, 0, 0, 2},
|
---|
106 | {1, 0, 0, -4, 0},
|
---|
107 | {-2, 0, 2, 2, 2},
|
---|
108 | {-1, 0, 2, 4, 2},
|
---|
109 | {2, 0, 0, -4, 0},
|
---|
110 | {1, 1, 2, -2, 2},
|
---|
111 | {1, 0, 2, 2, 1},
|
---|
112 | {-2, 0, 2, 4, 2},
|
---|
113 | {-1, 0, 4, 0, 2},
|
---|
114 | {1, -1, 0, -2, 0},
|
---|
115 | {2, 0, 2, -2, 1},
|
---|
116 | {2, 0, 2, 2, 2},
|
---|
117 | {1, 0, 0, 2, 1},
|
---|
118 | {0, 0, 4, -2, 2},
|
---|
119 | {3, 0, 2, -2, 2},
|
---|
120 | {1, 0, 2, -2, 0},
|
---|
121 | {0, 1, 2, 0, 1},
|
---|
122 | {-1, -1, 0, 2, 1},
|
---|
123 | {0, 0, -2, 0, 1},
|
---|
124 | {0, 0, 2, -1, 2},
|
---|
125 | {0, 1, 0, 2, 0},
|
---|
126 | {1, 0, -2, -2, 0},
|
---|
127 | {0, -1, 2, 0, 1},
|
---|
128 | {1, 1, 0, -2, 1},
|
---|
129 | {1, 0, -2, 2, 0},
|
---|
130 | {2, 0, 0, 2, 0},
|
---|
131 | {0, 0, 2, 4, 2},
|
---|
132 | {0, 1, 0, 1, 0}
|
---|
133 | };
|
---|
134 |
|
---|
135 | /* amplitudes which have secular terms; in 1/NUT_SCALE arc seconds
|
---|
136 | * {index, constant dPSI, T/10 in dPSI, constant in dEPS, T/10 in dEPS}
|
---|
137 | */
|
---|
138 | static long ampsecul[][5] = {
|
---|
139 | {0 ,-171996 ,-1742 ,92025 ,89},
|
---|
140 | {1 ,2062 ,2 ,-895 ,5},
|
---|
141 | {8 ,-13187 ,-16 ,5736 ,-31},
|
---|
142 | {9 ,1426 ,-34 ,54 ,-1},
|
---|
143 | {10 ,-517 ,12 ,224 ,-6},
|
---|
144 | {11 ,217 ,-5 ,-95 ,3},
|
---|
145 | {12 ,129 ,1 ,-70 ,0},
|
---|
146 | {15 ,17 ,-1 ,0 ,0},
|
---|
147 | {17 ,-16 ,1 ,7 ,0},
|
---|
148 | {30 ,-2274 ,-2 ,977 ,-5},
|
---|
149 | {31 ,712 ,1 ,-7 ,0},
|
---|
150 | {32 ,-386 ,-4 ,200 ,0},
|
---|
151 | {33 ,-301 ,0 ,129 ,-1},
|
---|
152 | {37 ,63 ,1 ,-33 ,0},
|
---|
153 | {38 ,-58 ,-1 ,32 ,0},
|
---|
154 | /* termination */ { -1, }
|
---|
155 | };
|
---|
156 |
|
---|
157 | /* amplitudes which only have constant terms; same unit as above
|
---|
158 | * {dPSI, dEPS}
|
---|
159 | * indexes which are already in ampsecul[][] are zeroed
|
---|
160 | */
|
---|
161 | static short ampconst[NUT_SERIES][2] = {
|
---|
162 | {0,0},
|
---|
163 | {0,0},
|
---|
164 | {46,-24},
|
---|
165 | {11,0},
|
---|
166 | {-3,1},
|
---|
167 | {-3,0},
|
---|
168 | {-2,1},
|
---|
169 | {1,0},
|
---|
170 | {0,0},
|
---|
171 | {0,0},
|
---|
172 | {0,0},
|
---|
173 | {0,0},
|
---|
174 | {0,0},
|
---|
175 | {48,1},
|
---|
176 | {-22,0},
|
---|
177 | {0,0},
|
---|
178 | {-15,9},
|
---|
179 | {0,0},
|
---|
180 | {-12,6},
|
---|
181 | {-6,3},
|
---|
182 | {-5,3},
|
---|
183 | {4,-2},
|
---|
184 | {4,-2},
|
---|
185 | {-4,0},
|
---|
186 | {1,0},
|
---|
187 | {1,0},
|
---|
188 | {-1,0},
|
---|
189 | {1,0},
|
---|
190 | {1,0},
|
---|
191 | {-1,0},
|
---|
192 | {0,0},
|
---|
193 | {0,0},
|
---|
194 | {0,0},
|
---|
195 | {0,0},
|
---|
196 | {-158,-1},
|
---|
197 | {123,-53},
|
---|
198 | {63,-2},
|
---|
199 | {0,0},
|
---|
200 | {0,0},
|
---|
201 | {-59,26},
|
---|
202 | {-51,27},
|
---|
203 | {-38,16},
|
---|
204 | {29,-1},
|
---|
205 | {29,-12},
|
---|
206 | {-31,13},
|
---|
207 | {26,-1},
|
---|
208 | {21,-10},
|
---|
209 | {16,-8},
|
---|
210 | {-13,7},
|
---|
211 | {-10,5},
|
---|
212 | {-7,0},
|
---|
213 | {7,-3},
|
---|
214 | {-7,3},
|
---|
215 | {-8,3},
|
---|
216 | {6,0},
|
---|
217 | {6,-3},
|
---|
218 | {-6,3},
|
---|
219 | {-7,3},
|
---|
220 | {6,-3},
|
---|
221 | {-5,3},
|
---|
222 | {5,0},
|
---|
223 | {-5,3},
|
---|
224 | {-4,0},
|
---|
225 | {4,0},
|
---|
226 | {-4,0},
|
---|
227 | {-3,0},
|
---|
228 | {3,0},
|
---|
229 | {-3,1},
|
---|
230 | {-3,1},
|
---|
231 | {-2,1},
|
---|
232 | {-3,1},
|
---|
233 | {-3,1},
|
---|
234 | {2,-1},
|
---|
235 | {-2,1},
|
---|
236 | {2,-1},
|
---|
237 | {-2,1},
|
---|
238 | {2,0},
|
---|
239 | {2,-1},
|
---|
240 | {1,-1},
|
---|
241 | {-1,0},
|
---|
242 | {1,-1},
|
---|
243 | {-2,1},
|
---|
244 | {-1,0},
|
---|
245 | {1,-1},
|
---|
246 | {-1,1},
|
---|
247 | {-1,1},
|
---|
248 | {1,0},
|
---|
249 | {1,0},
|
---|
250 | {1,-1},
|
---|
251 | {-1,0},
|
---|
252 | {-1,0},
|
---|
253 | {1,0},
|
---|
254 | {1,0},
|
---|
255 | {-1,0},
|
---|
256 | {1,0},
|
---|
257 | {1,0},
|
---|
258 | {-1,0},
|
---|
259 | {-1,0},
|
---|
260 | {-1,0},
|
---|
261 | {-1,0},
|
---|
262 | {-1,0},
|
---|
263 | {-1,0},
|
---|
264 | {-1,0},
|
---|
265 | {1,0},
|
---|
266 | {-1,0},
|
---|
267 | {1,0}
|
---|
268 | };
|
---|
269 |
|
---|
270 | /* given the modified JD, mj, find the nutation in obliquity, *deps, and
|
---|
271 | * the nutation in longitude, *dpsi, each in radians.
|
---|
272 | */
|
---|
273 | void
|
---|
274 | nutation (
|
---|
275 | double mj,
|
---|
276 | double *deps, /* on input: precision parameter in arc seconds */
|
---|
277 | double *dpsi)
|
---|
278 | {
|
---|
279 | static double lastmj = -10000, lastdeps, lastdpsi;
|
---|
280 | double T, T2, T3, T10; /* jul cent since J2000 */
|
---|
281 | double prec; /* series precis in arc sec */
|
---|
282 | int i, isecul; /* index in term table */
|
---|
283 | static double delcache[5][2*NUT_MAXMUL+1];
|
---|
284 | /* cache for multiples of delaunay args
|
---|
285 | * [M',M,F,D,Om][-min*x, .. , 0, .., max*x]
|
---|
286 | * make static to have unfilled fields cleared on init
|
---|
287 | */
|
---|
288 |
|
---|
289 | if (mj == lastmj) {
|
---|
290 | *deps = lastdeps;
|
---|
291 | *dpsi = lastdpsi;
|
---|
292 | return;
|
---|
293 | }
|
---|
294 |
|
---|
295 | prec = 0.0;
|
---|
296 |
|
---|
297 | #if 0 /* this is if deps should contain a precision value */
|
---|
298 | prec =* deps;
|
---|
299 | if (prec < 0.0 || prec > 1.0) /* accept only sane value */
|
---|
300 | prec = 1.0;
|
---|
301 | #endif
|
---|
302 |
|
---|
303 | /* augment for abundance of small terms */
|
---|
304 | prec *= NUT_SCALE/10;
|
---|
305 |
|
---|
306 | T = (mj - J2000)/36525.;
|
---|
307 | T2 = T * T;
|
---|
308 | T3 = T2 * T;
|
---|
309 | T10 = T/10.;
|
---|
310 |
|
---|
311 | /* calculate delaunay args and place in cache */
|
---|
312 | for (i = 0; i < 5; ++i) {
|
---|
313 | double x;
|
---|
314 | short j;
|
---|
315 |
|
---|
316 | x = delaunay[i][0] +
|
---|
317 | delaunay[i][1] * T +
|
---|
318 | delaunay[i][2] * T2 +
|
---|
319 | delaunay[i][3] * T3;
|
---|
320 |
|
---|
321 | /* convert to radians */
|
---|
322 | x /= SECPERCIRC;
|
---|
323 | x -= floor(x);
|
---|
324 | x *= 2.*PI;
|
---|
325 |
|
---|
326 | /* fill cache table */
|
---|
327 | for (j = 0; j <= 2*NUT_MAXMUL; ++j)
|
---|
328 | delcache[i][j] = (j - NUT_MAXMUL) * x;
|
---|
329 | }
|
---|
330 |
|
---|
331 | /* find dpsi and deps */
|
---|
332 | lastdpsi = lastdeps = 0.;
|
---|
333 | for (i = isecul = 0; i < NUT_SERIES ; ++i) {
|
---|
334 | double arg = 0., ampsin, ampcos;
|
---|
335 | short j;
|
---|
336 |
|
---|
337 | if (ampconst[i][0] || ampconst[i][1]) {
|
---|
338 | /* take non-secular terms from simple array */
|
---|
339 | ampsin = ampconst[i][0];
|
---|
340 | ampcos = ampconst[i][1];
|
---|
341 | } else {
|
---|
342 | /* secular terms from different array */
|
---|
343 | ampsin = ampsecul[isecul][1] + ampsecul[isecul][2] * T10;
|
---|
344 | ampcos = ampsecul[isecul][3] + ampsecul[isecul][4] * T10;
|
---|
345 | ++isecul;
|
---|
346 | }
|
---|
347 |
|
---|
348 | for (j = 0; j < 5; ++j)
|
---|
349 | arg += delcache[j][NUT_MAXMUL + multarg[i][j]];
|
---|
350 |
|
---|
351 | if (fabs(ampsin) >= prec)
|
---|
352 | lastdpsi += ampsin * sin(arg);
|
---|
353 |
|
---|
354 | if (fabs(ampcos) >= prec)
|
---|
355 | lastdeps += ampcos * cos(arg);
|
---|
356 |
|
---|
357 | }
|
---|
358 |
|
---|
359 | /* convert to radians.
|
---|
360 | */
|
---|
361 | lastdpsi = degrad(lastdpsi/3600./NUT_SCALE);
|
---|
362 | lastdeps = degrad(lastdeps/3600./NUT_SCALE);
|
---|
363 |
|
---|
364 | lastmj = mj;
|
---|
365 | *deps = lastdeps;
|
---|
366 | *dpsi = lastdpsi;
|
---|
367 | }
|
---|
368 |
|
---|
369 | /* given the modified JD, mj, correct, IN PLACE, the right ascension *ra
|
---|
370 | * and declination *dec (both in radians) for nutation.
|
---|
371 | */
|
---|
372 | void
|
---|
373 | nut_eq (double mj, double *ra, double *dec)
|
---|
374 | {
|
---|
375 | static double lastmj = -10000;
|
---|
376 | static double a[3][3]; /* rotation matrix */
|
---|
377 | double xold, yold, zold, x, y, z;
|
---|
378 |
|
---|
379 | if (mj != lastmj) {
|
---|
380 | double epsilon, dpsi, deps;
|
---|
381 | double se, ce, sp, cp, sede, cede;
|
---|
382 |
|
---|
383 | obliquity(mj, &epsilon);
|
---|
384 | nutation(mj, &deps, &dpsi);
|
---|
385 |
|
---|
386 | /* the rotation matrix a applies the nutation correction to
|
---|
387 | * a vector of equatoreal coordinates Xeq to Xeq' by 3 subsequent
|
---|
388 | * rotations: R1 - from equatoreal to ecliptic system by
|
---|
389 | * rotation of angle epsilon about x, R2 - rotate ecliptic
|
---|
390 | * system by -dpsi about its z, R3 - from ecliptic to equatoreal
|
---|
391 | * by rotation of angle -(epsilon + deps)
|
---|
392 | *
|
---|
393 | * Xeq' = A * Xeq = R3 * R2 * R1 * Xeq
|
---|
394 | *
|
---|
395 | * [ 1 0 0 ]
|
---|
396 | * R1 = [ 0 cos(eps) sin(eps) ]
|
---|
397 | * [ 0 - sin(eps) cos(eps) ]
|
---|
398 | *
|
---|
399 | * [ cos(dpsi) - sin(dpsi) 0 ]
|
---|
400 | * R2 = [ sin(dpsi) cos(dpsi) 0 ]
|
---|
401 | * [ 0 0 1 ]
|
---|
402 | *
|
---|
403 | * [ 1 0 0 ]
|
---|
404 | * R3 = [ 0 cos(eps + deps) - sin(eps + deps) ]
|
---|
405 | * [ 0 sin(eps + deps) cos(eps + deps) ]
|
---|
406 | *
|
---|
407 | * for efficiency, here is a explicitely:
|
---|
408 | */
|
---|
409 |
|
---|
410 | se = sin(epsilon);
|
---|
411 | ce = cos(epsilon);
|
---|
412 | sp = sin(dpsi);
|
---|
413 | cp = cos(dpsi);
|
---|
414 | sede = sin(epsilon + deps);
|
---|
415 | cede = cos(epsilon + deps);
|
---|
416 |
|
---|
417 | a[0][0] = cp;
|
---|
418 | a[0][1] = -sp*ce;
|
---|
419 | a[0][2] = -sp*se;
|
---|
420 |
|
---|
421 | a[1][0] = cede*sp;
|
---|
422 | a[1][1] = cede*cp*ce+sede*se;
|
---|
423 | a[1][2] = cede*cp*se-sede*ce;
|
---|
424 |
|
---|
425 | a[2][0] = sede*sp;
|
---|
426 | a[2][1] = sede*cp*ce-cede*se;
|
---|
427 | a[2][2] = sede*cp*se+cede*ce;
|
---|
428 |
|
---|
429 | lastmj = mj;
|
---|
430 | }
|
---|
431 |
|
---|
432 | sphcart(*ra, *dec, 1.0, &xold, &yold, &zold);
|
---|
433 | x = a[0][0] * xold + a[0][1] * yold + a[0][2] * zold;
|
---|
434 | y = a[1][0] * xold + a[1][1] * yold + a[1][2] * zold;
|
---|
435 | z = a[2][0] * xold + a[2][1] * yold + a[2][2] * zold;
|
---|
436 | cartsph(x, y, z, ra, dec, &zold); /* radius should be 1.0 */
|
---|
437 | if (*ra < 0.) *ra += 2.*PI; /* make positive for display */
|
---|
438 | }
|
---|
439 |
|
---|
440 | /* For RCS Only -- Do Not Edit */
|
---|
441 | static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: nutation.c,v $ $Date: 2008-03-25 17:45:17 $ $Revision: 1.7 $ $Name: not supported by cvs2svn $"};
|
---|