1 /*
2 * - - - - - - - - - - -
3 * g a l _ n u t 0 0 b
4 * - - - - - - - - - - -
5 *
6 * This routine is part of the General Astrodynamics Library
7 *
8 * Description:
9 *
10 * Nutation, IAU 2000B model.
11 *
12 * This routine is an independent translation of a FORTRAN routine
13 * that is part of IAU's SOFA software collection.
14 *
15 * Status:
16 *
17 * canonical model.
18 *
19 * Given:
20 *
21 * date1,date2 d TT as a 2-part Julian Date (Note 1)
22 *
23 * Returned:
24 *
25 * *dpsi,*deps d nutation, luni-solar + planetary (Note 2)
26 *
27 * Notes:
28 *
29 * 1) The TT date date1+date2 is a Julian Date, apportioned in any
30 * convenient way between the two arguments. For example,
31 * JD(TT)=2450123.7 could be expressed in any of these ways,
32 * among others:
33 *
34 * date1 date2
35 *
36 * 2450123.7 0.0 (JD method)
37 * 2451545.0 -1421.3 (J2000 method)
38 * 2400000.5 50123.2 (MJD method)
39 * 2450123.5 0.2 (date & time method)
40 *
41 * The JD method is the most natural and convenient to use in
42 * cases where the loss of several decimal digits of resolution
43 * is acceptable. The J2000 method is best matched to the way
44 * the argument is handled internally and will deliver the
45 * optimum resolution. The MJD method and the date & time methods
46 * are both good compromises between resolution and convenience.
47 *
48 * 2) The nutation components in longitude and obliquity are in radians
49 * and with respect to the equinox and ecliptic of date. The
50 * obliquity at J2000 is assumed to be the Lieske et al. (1977) value
51 * of 84381.448 arcsec. (The errors that result from using this
52 * routine with the IAU 2006 value of 84381.406 arcsec can be
53 * neglected.)
54 *
55 * The nutation model consists only of luni-solar terms, but includes
56 * also a fixed offset which compensates for certain long-period
57 * planetary terms (Note 7).
58 *
59 * 3) This routine is an implementation of the IAU 2000B abridged
60 * nutation model formally adopted by the IAU General Assembly in
61 * 2000. The routine computes the MHB_2000_SHORT luni-solar nutation
62 * series (Luzum 2001), but without the associated corrections for
63 * the precession rate adjustments and the offset between the GCRS
64 * and J2000 mean poles.
65 *
66 * 4) The full IAU 2000A (MHB2000) nutation model contains nearly 1400
67 * terms. The IAU 2000B model (McCarthy & Luzum 2003) contains only
68 * 77 terms, plus additional simplifications, yet still delivers
69 * results of 1 mas accuracy at present epochs. This combination of
70 * accuracy and size makes the IAU 2000B abridged nutation model
71 * suitable for most practical applications.
72 *
73 * The routine delivers a pole accurate to 1 mas from 1900 to 2100
74 * (usually better than 1 mas, very occasionally just outside 1 mas).
75 * The full IAU 2000A model, which is implemented in the routine
76 * gal_nut00a (q.v.), delivers considerably greater accuracy at
77 * current epochs; however, to realize this improved accuracy,
78 * corrections for the essentially unpredictable free-core-nutation
79 * (FCN) must also be included.
80 *
81 * 5) The present routine provides classical nutation. The
82 * MHB_2000_SHORT algorithm, from which it is adapted, deals also
83 * with (i) the offsets between the GCRS and mean poles and (ii) the
84 * adjustments in longitude and obliquity due to the changed
85 * precession rates. These additional functions, namely frame bias
86 * and precession adjustments, are supported by the SOFA routines
87 * gal_bi00 and gal_pr00.
88 *
89 * 6) The MHB_2000_SHORT algorithm also provides "total" nutations,
90 * comprising the arithmetic sum of the frame bias, precession
91 * adjustments, and nutation (luni-solar + planetary). These total
92 * nutations can be used in combination with an existing IAU 1976
93 * precession implementation, such as gal_pmat76, to deliver GCRS-to-
94 * true predictions of mas accuracy at current epochs. However, for
95 * symmetry with the gal_nut00a routine (q.v. for the reasons), the
96 * SOFA routines do not generate the "total nutations" directly.
97 * Should they be required, they could of course easily be generated
98 * by calling gal_bi00, gal_pr00 and the present routine and adding
99 * the results.
100 *
101 * 7) The IAU 2000B model includes "planetary bias" terms that are fixed
102 * in size but compensate for long-period nutations. The amplitudes
103 * quoted in McCarthy & Luzum (2003), namely Dpsi = -1.5835 mas and
104 * Depsilon = +1.6339 mas, are optimized for the "total nutations"
105 * method described in Note 6. The Luzum (2001) values used in this
106 * SOFA implementation, namely -0.135 mas and +0.388 mas, are
107 * optimized for the "rigorous" method, where frame bias, precession
108 * and nutation are applied separately and in that order. During the
109 * interval 1995-2050, the SOFA implementation delivers a maximum
110 * error of 1.001 mas (not including FCN).
111 *
112 * References:
113 *
114 * Lieske, J.H., Lederle, T., Fricke, W., Morando, B., "Expressions
115 * for the precession quantities based upon the IAU /1976/ system of
116 * astronomical constants", Astron.Astrophys. 58, 1-2, 1-16. (1977)
117 *
118 * Luzum, B., private communication, 2001 (Fortran code
119 * MHB_2000_SHORT)
120 *
121 * McCarthy, D.D. & Luzum, B.J., "An abridged model of the
122 * precession-nutation of the celestial pole", Cel.Mech.Dyn.Astron.
123 * 85, 37-49 (2003)
124 *
125 * Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
126 * Francou, G., Laskar, J., Astron.Astrophys. 282, 663-683 (1994)
127 *
128 * This revision:
129 *
130 * 2007 June 26 ( c version 2008 January 19 )
131 *
132 *
133 * Copyright (C) 2008 Paul C. L. Willmott. See notes at end.
134 *
135 *-----------------------------------------------------------------------
136 */
137
138 #include <math.h>
139 #include "gal_const.h"
140 #include "gal_nut00b.h"
141
142 void
gal_nut00b(double date1,double date2,double * dpsi,double * deps)143 gal_nut00b
144 (
145 double date1,
146 double date2,
147 double *dpsi,
148 double *deps
149 )
150 {
151
152 /*
153 * Miscellaneous
154 */
155
156 double t, el, elp, f, d, om, arg, dp, de, sarg, carg,
157 dpsils, depsls, dpsipl, depspl ;
158
159 int i ;
160
161 /* ----------------------------------------
162 * Tables of argument and term coefficients
163 * ----------------------------------------
164 */
165
166 /* -------------------------
167 * Luni-Solar nutation model
168 * -------------------------
169 */
170
171 /*
172 * Number of terms in the luni-solar nutation model
173 */
174
175 #define NLS 77
176
177 /*
178 * Coefficients for fundamental arguments
179 */
180
181 static const int NALS[ NLS ][5] = {
182
183 /*
184 * Luni-Solar argument multipliers:
185 *
186 * L L' F D Om
187 */
188
189 { 0, 0, 0, 0, 1, } ,
190 { 0, 0, 2, -2, 2, } ,
191 { 0, 0, 2, 0, 2, } ,
192 { 0, 0, 0, 0, 2, } ,
193 { 0, 1, 0, 0, 0, } ,
194 { 0, 1, 2, -2, 2, } ,
195 { 1, 0, 0, 0, 0, } ,
196 { 0, 0, 2, 0, 1, } ,
197 { 1, 0, 2, 0, 2, } ,
198 { 0, -1, 2, -2, 2, } ,
199 { 0, 0, 2, -2, 1, } ,
200 { -1, 0, 2, 0, 2, } ,
201 { -1, 0, 0, 2, 0, } ,
202 { 1, 0, 0, 0, 1, } ,
203 { -1, 0, 0, 0, 1, } ,
204 { -1, 0, 2, 2, 2, } ,
205 { 1, 0, 2, 0, 1, } ,
206 { -2, 0, 2, 0, 1, } ,
207 { 0, 0, 0, 2, 0, } ,
208 { 0, 0, 2, 2, 2, } ,
209 { 0, -2, 2, -2, 2, } ,
210 { -2, 0, 0, 2, 0, } ,
211 { 2, 0, 2, 0, 2, } ,
212 { 1, 0, 2, -2, 2, } ,
213 { -1, 0, 2, 0, 1, } ,
214 { 2, 0, 0, 0, 0, } ,
215 { 0, 0, 2, 0, 0, } ,
216 { 0, 1, 0, 0, 1, } ,
217 { -1, 0, 0, 2, 1, } ,
218 { 0, 2, 2, -2, 2, } ,
219 { 0, 0, -2, 2, 0, } ,
220 { 1, 0, 0, -2, 1, } ,
221 { 0, -1, 0, 0, 1, } ,
222 { -1, 0, 2, 2, 1, } ,
223 { 0, 2, 0, 0, 0, } ,
224 { 1, 0, 2, 2, 2, } ,
225 { -2, 0, 2, 0, 0, } ,
226 { 0, 1, 2, 0, 2, } ,
227 { 0, 0, 2, 2, 1, } ,
228 { 0, -1, 2, 0, 2, } ,
229 { 0, 0, 0, 2, 1, } ,
230 { 1, 0, 2, -2, 1, } ,
231 { 2, 0, 2, -2, 2, } ,
232 { -2, 0, 0, 2, 1, } ,
233 { 2, 0, 2, 0, 1, } ,
234 { 0, -1, 2, -2, 1, } ,
235 { 0, 0, 0, -2, 1, } ,
236 { -1, -1, 0, 2, 0, } ,
237 { 2, 0, 0, -2, 1, } ,
238 { 1, 0, 0, 2, 0, } ,
239 { 0, 1, 2, -2, 1, } ,
240 { 1, -1, 0, 0, 0, } ,
241 { -2, 0, 2, 0, 2, } ,
242 { 3, 0, 2, 0, 2, } ,
243 { 0, -1, 0, 2, 0, } ,
244 { 1, -1, 2, 0, 2, } ,
245 { 0, 0, 0, 1, 0, } ,
246 { -1, -1, 2, 2, 2, } ,
247 { -1, 0, 2, 0, 0, } ,
248 { 0, -1, 2, 2, 2, } ,
249 { -2, 0, 0, 0, 1, } ,
250 { 1, 1, 2, 0, 2, } ,
251 { 2, 0, 0, 0, 1, } ,
252 { -1, 1, 0, 1, 0, } ,
253 { 1, 1, 0, 0, 0, } ,
254 { 1, 0, 2, 0, 0, } ,
255 { -1, 0, 2, -2, 1, } ,
256 { 1, 0, 0, 0, 2, } ,
257 { -1, 0, 0, 1, 0, } ,
258 { 0, 0, 2, 1, 2, } ,
259 { -1, 0, 2, 4, 2, } ,
260 { -1, 1, 0, 1, 1, } ,
261 { 0, -2, 2, -2, 1, } ,
262 { 1, 0, 2, 2, 1, } ,
263 { -2, 0, 2, 2, 2, } ,
264 { -1, 0, 0, 0, 2, } ,
265 { 1, 1, 2, -2, 2, } ,
266
267 } ;
268
269 /*
270 * Longitude and obliquity coefficients
271 */
272
273 static const double CLS[ NLS ][6] = {
274
275 /*
276 * Luni-Solar nutation coefficients, unit 1e-7 arcsec:
277 * longitude (sin, t*sin, cos), obliquity (cos, t*cos, sin)
278 */
279
280 { -172064161e0, -174666e0, 33386e0, 92052331e0, 9086e0, 15377e0, } ,
281 { -13170906e0, -1675e0, -13696e0, 5730336e0, -3015e0, -4587e0, } ,
282 { -2276413e0, -234e0, 2796e0, 978459e0, -485e0, 1374e0, } ,
283 { 2074554e0, 207e0, -698e0, -897492e0, 470e0, -291e0, } ,
284 { 1475877e0, -3633e0, 11817e0, 73871e0, -184e0, -1924e0, } ,
285 { -516821e0, 1226e0, -524e0, 224386e0, -677e0, -174e0, } ,
286 { 711159e0, 73e0, -872e0, -6750e0, 0e0, 358e0, } ,
287 { -387298e0, -367e0, 380e0, 200728e0, 18e0, 318e0, } ,
288 { -301461e0, -36e0, 816e0, 129025e0, -63e0, 367e0, } ,
289 { 215829e0, -494e0, 111e0, -95929e0, 299e0, 132e0, } ,
290 { 128227e0, 137e0, 181e0, -68982e0, -9e0, 39e0, } ,
291 { 123457e0, 11e0, 19e0, -53311e0, 32e0, -4e0, } ,
292 { 156994e0, 10e0, -168e0, -1235e0, 0e0, 82e0, } ,
293 { 63110e0, 63e0, 27e0, -33228e0, 0e0, -9e0, } ,
294 { -57976e0, -63e0, -189e0, 31429e0, 0e0, -75e0, } ,
295 { -59641e0, -11e0, 149e0, 25543e0, -11e0, 66e0, } ,
296 { -51613e0, -42e0, 129e0, 26366e0, 0e0, 78e0, } ,
297 { 45893e0, 50e0, 31e0, -24236e0, -10e0, 20e0, } ,
298 { 63384e0, 11e0, -150e0, -1220e0, 0e0, 29e0, } ,
299 { -38571e0, -1e0, 158e0, 16452e0, -11e0, 68e0, } ,
300 { 32481e0, 0e0, 0e0, -13870e0, 0e0, 0e0, } ,
301 { -47722e0, 0e0, -18e0, 477e0, 0e0, -25e0, } ,
302 { -31046e0, -1e0, 131e0, 13238e0, -11e0, 59e0, } ,
303 { 28593e0, 0e0, -1e0, -12338e0, 10e0, -3e0, } ,
304 { 20441e0, 21e0, 10e0, -10758e0, 0e0, -3e0, } ,
305 { 29243e0, 0e0, -74e0, -609e0, 0e0, 13e0, } ,
306 { 25887e0, 0e0, -66e0, -550e0, 0e0, 11e0, } ,
307 { -14053e0, -25e0, 79e0, 8551e0, -2e0, -45e0, } ,
308 { 15164e0, 10e0, 11e0, -8001e0, 0e0, -1e0, } ,
309 { -15794e0, 72e0, -16e0, 6850e0, -42e0, -5e0, } ,
310 { 21783e0, 0e0, 13e0, -167e0, 0e0, 13e0, } ,
311 { -12873e0, -10e0, -37e0, 6953e0, 0e0, -14e0, } ,
312 { -12654e0, 11e0, 63e0, 6415e0, 0e0, 26e0, } ,
313 { -10204e0, 0e0, 25e0, 5222e0, 0e0, 15e0, } ,
314 { 16707e0, -85e0, -10e0, 168e0, -1e0, 10e0, } ,
315 { -7691e0, 0e0, 44e0, 3268e0, 0e0, 19e0, } ,
316 { -11024e0, 0e0, -14e0, 104e0, 0e0, 2e0, } ,
317 { 7566e0, -21e0, -11e0, -3250e0, 0e0, -5e0, } ,
318 { -6637e0, -11e0, 25e0, 3353e0, 0e0, 14e0, } ,
319 { -7141e0, 21e0, 8e0, 3070e0, 0e0, 4e0, } ,
320 { -6302e0, -11e0, 2e0, 3272e0, 0e0, 4e0, } ,
321 { 5800e0, 10e0, 2e0, -3045e0, 0e0, -1e0, } ,
322 { 6443e0, 0e0, -7e0, -2768e0, 0e0, -4e0, } ,
323 { -5774e0, -11e0, -15e0, 3041e0, 0e0, -5e0, } ,
324 { -5350e0, 0e0, 21e0, 2695e0, 0e0, 12e0, } ,
325 { -4752e0, -11e0, -3e0, 2719e0, 0e0, -3e0, } ,
326 { -4940e0, -11e0, -21e0, 2720e0, 0e0, -9e0, } ,
327 { 7350e0, 0e0, -8e0, -51e0, 0e0, 4e0, } ,
328 { 4065e0, 0e0, 6e0, -2206e0, 0e0, 1e0, } ,
329 { 6579e0, 0e0, -24e0, -199e0, 0e0, 2e0, } ,
330 { 3579e0, 0e0, 5e0, -1900e0, 0e0, 1e0, } ,
331 { 4725e0, 0e0, -6e0, -41e0, 0e0, 3e0, } ,
332 { -3075e0, 0e0, -2e0, 1313e0, 0e0, -1e0, } ,
333 { -2904e0, 0e0, 15e0, 1233e0, 0e0, 7e0, } ,
334 { 4348e0, 0e0, -10e0, -81e0, 0e0, 2e0, } ,
335 { -2878e0, 0e0, 8e0, 1232e0, 0e0, 4e0, } ,
336 { -4230e0, 0e0, 5e0, -20e0, 0e0, -2e0, } ,
337 { -2819e0, 0e0, 7e0, 1207e0, 0e0, 3e0, } ,
338 { -4056e0, 0e0, 5e0, 40e0, 0e0, -2e0, } ,
339 { -2647e0, 0e0, 11e0, 1129e0, 0e0, 5e0, } ,
340 { -2294e0, 0e0, -10e0, 1266e0, 0e0, -4e0, } ,
341 { 2481e0, 0e0, -7e0, -1062e0, 0e0, -3e0, } ,
342 { 2179e0, 0e0, -2e0, -1129e0, 0e0, -2e0, } ,
343 { 3276e0, 0e0, 1e0, -9e0, 0e0, 0e0, } ,
344 { -3389e0, 0e0, 5e0, 35e0, 0e0, -2e0, } ,
345 { 3339e0, 0e0, -13e0, -107e0, 0e0, 1e0, } ,
346 { -1987e0, 0e0, -6e0, 1073e0, 0e0, -2e0, } ,
347 { -1981e0, 0e0, 0e0, 854e0, 0e0, 0e0, } ,
348 { 4026e0, 0e0, -353e0, -553e0, 0e0, -139e0, } ,
349 { 1660e0, 0e0, -5e0, -710e0, 0e0, -2e0, } ,
350 { -1521e0, 0e0, 9e0, 647e0, 0e0, 4e0, } ,
351 { 1314e0, 0e0, 0e0, -700e0, 0e0, 0e0, } ,
352 { -1283e0, 0e0, 0e0, 672e0, 0e0, 0e0, } ,
353 { -1331e0, 0e0, 8e0, 663e0, 0e0, 4e0, } ,
354 { 1383e0, 0e0, -2e0, -594e0, 0e0, -2e0, } ,
355 { 1405e0, 0e0, 4e0, -610e0, 0e0, 2e0, } ,
356 { 1290e0, 0e0, 0e0, -556e0, 0e0, 0e0, } ,
357
358 } ;
359
360 /* ---------------------------------------
361 * Fixed offset in lieu of planetary terms (radians)
362 * ---------------------------------------
363 */
364
365 const double DPPLAN = -0.135 * GAL_MAS2R ;
366 const double DEPLAN = +0.388 * GAL_MAS2R ;
367
368 /*
369 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
370 */
371
372 /*
373 * Interval between fundamental epoch J2000.0 and given date (JC).
374 */
375
376 t = ( ( date1 - GAL_J2000 ) + date2 ) / GAL_DJC ;
377
378 /*
379 * -------------------
380 * LUNI-SOLAR NUTATION
381 * -------------------
382 */
383
384 /*
385 * Fundamental (Delaunay) arguments from Simon et al. (1994)
386 */
387
388 /*
389 * Mean anomaly of the Moon.
390 */
391
392 el = fmod ( 485868.249036 +
393 ( + 1717915923.2178 ) * t, GAL_TURNAS ) * GAL_AS2R ;
394
395 /*
396 * Mean anomaly of the Sun.
397 */
398
399 elp = fmod ( 1287104.79305 +
400 ( + 129596581.0481 ) * t, GAL_TURNAS ) * GAL_AS2R ;
401
402 /*
403 * Mean argument of the latitude of the Moon.
404 */
405
406 f = fmod ( 335779.526232 +
407 ( + 1739527262.8478 ) * t, GAL_TURNAS ) * GAL_AS2R ;
408
409 /*
410 * Mean elongation of the Moon from the Sun.
411 */
412
413 d = fmod ( 1072260.70369 +
414 ( + 1602961601.2090 ) * t, GAL_TURNAS ) * GAL_AS2R ;
415
416 /*
417 * Mean longitude of the ascending node of the Moon.
418 */
419
420 om = fmod ( 450160.398036 +
421 ( - 6962890.5431 ) * t, GAL_TURNAS ) * GAL_AS2R ;
422
423 /*
424 * Initialize the nutation values.
425 */
426
427 dp = 0.0 ;
428 de = 0.0 ;
429
430 /*
431 * Summation of luni-solar nutation series (in reverse order).
432 */
433
434 for ( i = NLS - 1; i >= 0; i-- ) {
435
436 /*
437 * Argument and functions.
438 */
439
440 arg = fmod ( ( double ) NALS[i][0] * el +
441 ( double ) NALS[i][1] * elp +
442 ( double ) NALS[i][2] * f +
443 ( double ) NALS[i][3] * d +
444 ( double ) NALS[i][4] * om, GAL_2PI ) ;
445
446 sarg = sin ( arg ) ;
447 carg = cos ( arg ) ;
448
449 /*
450 * Term.
451 */
452
453 dp += ( CLS[i][0] + CLS[i][1] * t ) * sarg + CLS[i][2] * carg ;
454 de += ( CLS[i][3] + CLS[i][4] * t ) * carg + CLS[i][5] * sarg ;
455
456 }
457
458 /*
459 * Convert from 0.1 microarcsec units to radians.
460 */
461
462 dpsils = dp * GAL_U2R ;
463 depsls = de * GAL_U2R ;
464
465 /*
466 * -----------------------------
467 * IN LIEU OF PLANETARY NUTATION
468 * -----------------------------
469 */
470
471 /*
472 * Fixed offset to correct for missing terms in truncated series.
473 */
474
475 dpsipl = DPPLAN ;
476 depspl = DEPLAN ;
477
478 /*
479 * -------
480 * RESULTS
481 * -------
482 */
483
484 /*
485 * Add luni-solar and planetary components.
486 */
487
488 *dpsi = dpsils + dpsipl ;
489 *deps = depsls + depspl ;
490
491 /*
492 * Finished.
493 */
494
495 }
496
497 /*
498 * gal - General Astrodynamics Library
499 * Copyright (C) 2008 Paul C. L. Willmott
500 *
501 * This program is free software; you can redistribute it and/or modify
502 * it under the terms of the GNU General Public License as published by
503 * the Free Software Foundation; either version 2 of the License, or
504 * (at your option) any later version.
505 *
506 * This program is distributed in the hope that it will be useful,
507 * but WITHOUT ANY WARRANTY; without even the implied warranty of
508 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
509 * GNU General Public License for more details.
510 *
511 * You should have received a copy of the GNU General Public License along
512 * with this program; if not, write to the Free Software Foundation, Inc.,
513 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
514 *
515 * Contact:
516 *
517 * Paul Willmott
518 * vp9mu@amsat.org
519 */
520