source: Sophya/trunk/SophyaExt/XephemAstroLib/nutation.c@ 3586

Last change on this file since 3586 was 3477, checked in by cmv, 18 years ago

mise a jour Xephem 3.7.3 , cmv 25/03/2008

File size: 9.9 KB
Line 
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 */
16static 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 */
25static 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 */
138static 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 */
161static 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 */
273void
274nutation (
275double mj,
276double *deps, /* on input: precision parameter in arc seconds */
277double *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 */
372void
373nut_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 */
441static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: nutation.c,v $ $Date: 2008-03-25 17:45:17 $ $Revision: 1.7 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.