source: trunk/MagicSoft/slalib/nutc.c

Last change on this file was 731, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 12.6 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaNutc ( double date, double *dpsi, double *deps, double *eps0 )
4/*
5** - - - - - - - -
6** s l a N u t c
7** - - - - - - - -
8**
9** Nutation: longitude & obliquity components and
10** mean obliquity (IAU 1980 theory).
11**
12** (double precision)
13**
14** References:
15** Final report of the IAU working group on nutation,
16** chairman P.K.Seidelmann, 1980.
17** Kaplan,G.H., 1981, USNO circular No. 163, pa3-6.
18**
19** Given:
20** date double TDB (loosely ET) as Modified Julian Date
21** (JD-2400000.5)
22**
23** Returned:
24** *dpsi,*deps double nutation in longitude,obliquity
25** *eps0 double mean obliquity
26**
27** Called: slaDrange
28**
29** Defined in slamac.h: DAS2R, dmod
30**
31** Last revision: 19 March 1996
32**
33** Copyright P.T.Wallace. All rights reserved.
34*/
35
36#define T2AS 1296000.0 /* Turns to arc seconds */
37#define U2R 0.4848136811095359949e-9 /* Units of 0.0001 arcsec to radians */
38
39{
40 double t, el, el2, el3, elp, elp2,
41 f, f2, f4,
42 d, d2, d4,
43 om, om2,
44 dp, de, a;
45
46/* Interval between basic epoch J2000.0 and current epoch (JC) */
47 t = ( date - 51544.5 ) / 36525.0;
48
49/* Fundamental arguments in the FK5 reference system */
50
51/* Mean longitude of the Moon minus mean longitude of the Moon's perigee */
52 el = slaDrange ( DAS2R * dmod ( 485866.733 + ( 1325.0 * T2AS + 715922.633
53 + ( 31.310 + 0.064 * t ) * t ) * t , T2AS ) );
54
55/* Mean longitude of the Sun minus mean longitude of the Sun's perigee */
56 elp = slaDrange ( DAS2R * dmod ( 1287099.804 + ( 99.0 * T2AS + 1292581.224
57 + ( -0.577 - 0.012 * t ) * t ) * t, T2AS ) );
58
59/* Mean longitude of the Moon minus mean longitude of the Moon's node */
60 f = slaDrange ( DAS2R * dmod ( 335778.877 + ( 1342.0 * T2AS + 295263.137
61 + ( -13.257 + 0.011 * t ) * t ) * t, T2AS ) );
62
63/* Mean elongation of the Moon from the Sun */
64 d = slaDrange ( DAS2R * dmod ( 1072261.307 + ( 1236.0 * T2AS + 1105601.328
65 + ( -6.891 + 0.019 * t ) * t ) * t, T2AS ) );
66
67/* Longitude of the mean ascending node of the lunar orbit on the
68 ecliptic, measured from the mean equinox of date */
69 om = slaDrange ( DAS2R * dmod ( 450160.280 + ( -5.0 * T2AS - 482890.539
70 + ( 7.455 + 0.008 * t ) * t ) * t, T2AS ) );
71
72/* Multiples of arguments */
73 el2 = el + el;
74 el3 = el2 + el;
75 elp2 = elp + elp;
76 f2 = f + f;
77 f4 = f2 + f2;
78 d2 = d + d;
79 d4 = d2 + d2;
80 om2 = om + om;
81
82/* Series for the nutation */
83 dp = 0.0;
84 de = 0.0;
85
86 dp += sin ( elp + d ); /* 106 */
87
88 dp -= sin ( f2 + d4 + om2 ); /* 105 */
89
90 dp += sin ( el2 + d2 ); /* 104 */
91
92 dp -= sin ( el - f2 + d2 ); /* 103 */
93
94 dp -= sin ( el + elp - d2 + om ); /* 102 */
95
96 dp -= sin ( - elp + f2 + om ); /* 101 */
97
98 dp -= sin ( el - f2 - d2 ); /* 100 */
99
100 dp -= sin ( elp + d2 ); /* 99 */
101
102 dp -= sin ( f2 - d + om2 ); /* 98 */
103
104 dp -= sin ( - f2 + om ); /* 97 */
105
106 dp += sin ( - el - elp + d2 + om ); /* 96 */
107
108 dp += sin ( elp + f2 + om ); /* 95 */
109
110 dp -= sin ( el + f2 - d2 ); /* 94 */
111
112 dp += sin ( el3 + f2 - d2 + om2 ); /* 93 */
113
114 dp += sin ( f4 - d2 + om2 ); /* 92 */
115
116 dp -= sin ( el + d2 + om ); /* 91 */
117
118 dp -= sin ( el2 + f2 + d2 + om2 ); /* 90 */
119
120 a = el2 + f2 - d2 + om; /* 89 */
121 dp += sin ( a );
122 de -= cos ( a );
123
124 dp += sin ( el - elp - d2 ); /* 88 */
125
126 dp += sin ( - el + f4 + om2 ); /* 87 */
127
128 a = - el2 + f2 + d4 + om2; /* 86 */
129 dp -= sin ( a );
130 de += cos ( a );
131
132 a = el + f2 + d2 + om; /* 85 */
133 dp -= sin ( a );
134 de += cos ( a );
135
136 a = el + elp + f2 - d2 + om2; /* 84 */
137 dp += sin ( a );
138 de -= cos ( a );
139
140 dp -= sin ( el2 - d4 ); /* 83 */
141
142 a = - el + f2 + d4 + om2; /* 82 */
143 dp -= 2.0 * sin ( a );
144 de += cos ( a );
145
146 a = - el2 + f2 + d2 + om2; /* 81 */
147 dp += sin ( a );
148 de = de - cos ( a );
149
150 dp -= sin ( el - d4 ); /* 80 */
151
152 a = - el + om2; /* 79 */
153 dp += sin ( a );
154 de = de - cos ( a );
155
156 a = f2 + d + om2; /* 78 */
157 dp += 2.0 * sin ( a );
158 de = de - cos ( a );
159
160 dp += 2.0 * sin ( el3 ); /* 77 */
161
162 a = el + om2; /* 76 */
163 dp -= 2.0 * sin ( a );
164 de += cos ( a );
165
166 a = el2 + om; /* 75 */
167 dp += 2.0 * sin ( a );
168 de -= cos ( a );
169
170 a = - el + f2 - d2 + om; /* 74 */
171 dp -= 2.0 * sin ( a );
172 de += cos ( a );
173
174 a = el + elp + f2 + om2; /* 73 */
175 dp += 2.0 * sin ( a );
176 de = de - cos ( a );
177
178 a = - elp + f2 + d2 + om2; /* 72 */
179 dp -= 3.0 * sin ( a );
180 de += cos ( a );
181
182 a = el3 + f2 + om2; /* 71 */
183 dp -= 3.0 * sin ( a );
184 de += cos ( a );
185
186 a = - el2 + om; /* 70 */
187 dp -= 2.0 * sin ( a );
188 de += cos ( a );
189
190 a = - el - elp + f2 + d2 + om2; /* 69 */
191 dp -= 3.0 * sin ( a );
192 de += cos ( a );
193
194 a = el - elp + f2 + om2; /* 68 */
195 dp -= 3.0 * sin ( a );
196 de += cos ( a );
197
198 dp += 3.0 * sin ( el + f2 ); /* 67 */
199
200 dp -= 3.0 * sin ( el + elp ); /* 66 */
201
202 dp -= 4.0 * sin ( d ); /* 65 */
203
204 dp += 4.0 * sin ( el - f2 ); /* 64 */
205
206 dp -= 4.0 * sin ( elp - d2 ); /* 63 */
207
208 a = el2 + f2 + om; /* 62 */
209 dp -= 5.0 * sin ( a );
210 de += 3.0 * cos ( a );
211
212 dp += 5.0 * sin ( el - elp ); /* 61 */
213
214 a = - d2 + om; /* 60 */
215 dp -= 5.0 * sin ( a );
216 de += 3.0 * cos ( a );
217
218 a = el + f2 - d2 + om; /* 59 */
219 dp += 6.0 * sin ( a );
220 de -= 3.0 * cos ( a );
221
222 a = f2 + d2 + om; /* 58 */
223 dp -= 7.0 * sin ( a );
224 de += 3.0 * cos ( a );
225
226 a = d2 + om; /* 57 */
227 dp -= 6.0 * sin ( a );
228 de += 3.0 * cos ( a );
229
230 a = el2 + f2 - d2 + om2; /* 56 */
231 dp += 6.0 * sin ( a );
232 de -= 3.0 * cos ( a );
233
234 dp += 6.0 * sin ( el + d2); /* 55 */
235
236 a = el + f2 + d2 + om2; /* 54 */
237 dp -= 8.0 * sin ( a );
238 de += 3.0 * cos ( a );
239
240 a = - elp + f2 + om2; /* 53 */
241 dp -= 7.0 * sin ( a );
242 de += 3.0 * cos ( a );
243
244 a = elp + f2 + om2; /* 52 */
245 dp += 7.0 * sin ( a );
246 de -= 3.0 * cos ( a );
247
248 dp -= 7.0 * sin ( el + elp - d2 ); /* 51 */
249
250 a = - el + f2 + d2 + om; /* 50 */
251 dp -= 10.0 * sin ( a );
252 de += 5.0 * cos ( a );
253
254 a = el - d2 + om; /* 49 */
255 dp -= 13.0 * sin ( a );
256 de += 7.0 * cos ( a );
257
258 a = - el + d2 + om; /* 48 */
259 dp += 16.0 * sin ( a );
260 de -= 8.0 * cos ( a );
261
262 a = - el + f2 + om; /* 47 */
263 dp += 21.0 * sin ( a );
264 de -= 10.0 * cos ( a );
265
266 dp += 26.0 * sin ( f2 ); /* 46 */
267 de -= cos( f2 );
268
269 a = el2 + f2 + om2; /* 45 */
270 dp -= 31.0 * sin ( a );
271 de += 13.0 * cos ( a );
272
273 a = el + f2 - d2 + om2; /* 44 */
274 dp += 29.0 * sin ( a );
275 de -= 12.0 * cos ( a );
276
277 dp += 29.0 * sin ( el2 ); /* 43 */
278 de -= cos( el2 );
279
280 a = f2 + d2 + om2; /* 42 */
281 dp -= 38.0 * sin ( a );
282 de += 16.0 * cos ( a );
283
284 a = el + f2 + om; /* 41 */
285 dp -= 51.0 * sin ( a );
286 de += 27.0 * cos ( a );
287
288 a = - el + f2 + d2 + om2; /* 40 */
289 dp -= 59.0 * sin ( a );
290 de += 26.0 * cos ( a );
291
292 a = - el + om; /* 39 */
293 dp += ( - 58.0 - 0.1 * t ) * sin ( a );
294 de += 32.0 * cos ( a );
295
296 a = el + om; /* 38 */
297 dp += ( 63.0 + 0.1 * t ) * sin ( a );
298 de -= 33.0 * cos ( a );
299
300 dp += 63.0 * sin ( d2 ); /* 37 */
301 de -= 2.0 * cos( d2 );
302
303 a = - el + f2 + om2; /* 36 */
304 dp += 123.0 * sin ( a );
305 de -= 53.0 * cos ( a );
306
307 a = el - d2; /* 35 */
308 dp -= 158.0 * sin ( a );
309 de -= cos ( a );
310
311 a = el + f2 + om2; /* 34 */
312 dp -= 301.0 * sin ( a );
313 de += ( 129.0 - 0.1 * t ) * cos ( a );
314
315 a = f2 + om; /* 33 */
316 dp += ( - 386.0 - 0.4 * t ) * sin ( a );
317 de += 200.0 * cos ( a );
318
319 dp += ( 712.0 + 0.1 * t ) * sin ( el ); /* 32 */
320 de -= 7.0 * cos( el );
321
322 a = f2 + om2; /* 31 */
323 dp += ( -2274.0 - 0.2 * t ) * sin ( a );
324 de += ( 977.0 - 0.5 * t ) * cos ( a );
325
326 dp -= sin ( elp + f2 - d2 ); /* 30 */
327
328 dp += sin ( - el + d + om ); /* 29 */
329
330 dp += sin ( elp + om2 ); /* 28 */
331
332 dp -= sin ( elp - f2 + d2 ); /* 27 */
333
334 dp += sin ( - f2 + d2 + om ); /* 26 */
335
336 dp += sin ( el2 + elp - d2 ); /* 25 */
337
338 dp -= 4.0 * sin ( el - d ); /* 24 */
339
340 a = elp + f2 - d2 + om; /* 23 */
341 dp += 4.0 * sin ( a );
342 de -= 2.0 * cos ( a );
343
344 a = el2 - d2 + om; /* 22 */
345 dp += 4.0 * sin ( a );
346 de -= 2.0 * cos ( a );
347
348 a = - elp + f2 - d2 + om; /* 21 */
349 dp -= 5.0 * sin ( a );
350 de += 3.0 * cos ( a );
351
352 a = - el2 + d2 + om; /* 20 */
353 dp -= 6.0 * sin ( a );
354 de += 3.0 * cos ( a );
355
356 a = - elp + om; /* 19 */
357 dp -= 12.0 * sin ( a );
358 de += 6.0 * cos ( a );
359
360 a = elp2 + f2 - d2 + om2; /* 18 */
361 dp += ( - 16.0 + 0.1 * t) * sin ( a );
362 de += 7.0 * cos ( a );
363
364 a = elp + om; /* 17 */
365 dp -= 15.0 * sin ( a );
366 de += 9.0 * cos ( a );
367
368 dp += ( 17.0 - 0.1 * t ) * sin ( elp2 ); /* 16 */
369
370 dp -= 22.0 * sin ( f2 - d2 ); /* 15 */
371
372 a = el2 - d2; /* 14 */
373 dp += 48.0 * sin ( a );
374 de += cos ( a );
375
376 a = f2 - d2 + om; /* 13 */
377 dp += ( 129.0 + 0.1 * t ) * sin ( a );
378 de -= 70.0 * cos ( a );
379
380 a = - elp + f2 - d2 + om2; /* 12 */
381 dp += ( 217.0 - 0.5 * t ) * sin ( a );
382 de += ( -95.0 + 0.3 * t ) * cos ( a );
383
384 a = elp + f2 - d2 + om2; /* 11 */
385 dp += ( - 517.0 + 1.2 * t ) * sin ( a );
386 de += ( 224.0 - 0.6 * t ) * cos ( a );
387
388 dp += ( 1426.0 - 3.4 * t ) * sin ( elp ); /* 10 */
389 de += ( 54.0 - 0.1 * t) * cos ( elp );
390
391 a = f2 - d2 + om2; /* 9 */
392 dp += ( - 13187.0 - 1.6 * t ) * sin ( a );
393 de += ( 5736.0 - 3.1 * t ) * cos ( a );
394
395 dp += sin ( el2 - f2 + om ); /* 8 */
396
397 a = - elp2 + f2 - d2 + om; /* 7 */
398 dp -= 2.0 * sin ( a );
399 de += cos ( a );
400
401 dp -= 3.0 * sin ( el - elp - d ); /* 6 */
402
403 a = - el2 + f2 + om2; /* 5 */
404 dp -= 3.0 * sin ( a );
405 de += cos ( a );
406
407 dp += 11.0 * sin ( el2 - f2 ); /* 4 */
408
409 a = - el2 + f2 + om; /* 3 */
410 dp += 46.0 * sin ( a );
411 de -= 24.0 * cos ( a );
412
413 dp += ( 2062.0 + 0.2 * t ) * sin ( om2 ); /* 2 */
414 de += ( - 895.0 + 0.5 * t ) * cos ( om2 );
415
416 dp += ( - 171996.0 - 174.2 * t) * sin ( om ); /* 1 */
417 de += ( 92025.0 + 8.9 * t ) * cos ( om );
418
419/* Convert results to radians */
420 *dpsi = dp * U2R;
421 *deps = de * U2R;
422
423/* Mean obliquity */
424 *eps0 = DAS2R * ( 84381.448 +
425 ( - 46.8150 +
426 ( - 0.00059 + 0.001813 * t ) * t ) * t );
427}
Note: See TracBrowser for help on using the repository browser.