1 /*
2 * - - - - - - - -
3 * g a l _ s 0 6
4 * - - - - - - - -
5 *
6 * This routine is part of the General Astrodynamics Library
7 *
8 * Description:
9 *
10 * The CIO locator s, positioning the Celestial Intermediate Origin on
11 * the equator of the Celestial Intermediate Pole, given the CIP's X,Y
12 * coordinates. Compatible with IAU 2006/2000A precession-nutation.
13 *
14 * This routine is an independent translation of a FORTRAN routine
15 * that is part of IAU's SOFA software collection.
16 *
17 * Status:
18 *
19 * canonical model.
20 *
21 * Given:
22 *
23 * date1,date2 d TT as a 2-part Julian Date (Note 1)
24 * x,y d CIP coordinates (Note 3)
25 *
26 * Returned:
27 *
28 * gal_s06 d the CIO locator s in radians (Note 2)
29 *
30 * Notes:
31 *
32 * 1) The TT date date1+date2 is a Julian Date, apportioned in any
33 * convenient way between the two arguments. For example,
34 * JD(TT)=2450123.7 could be expressed in any of these ways,
35 * among others:
36 *
37 * date1 date2
38 *
39 * 2450123.7 0.0 (JD method)
40 * 2451545.0 -1421.3 (J2000 method)
41 * 2400000.5 50123.2 (MJD method)
42 * 2450123.5 0.2 (date & time method)
43 *
44 * The JD method is the most natural and convenient to use in
45 * cases where the loss of several decimal digits of resolution
46 * is acceptable. The J2000 method is best matched to the way
47 * the argument is handled internally and will deliver the
48 * optimum resolution. The MJD method and the date & time methods
49 * are both good compromises between resolution and convenience.
50 *
51 * 2) The CIO locator s is the difference between the right ascensions
52 * of the same point in two systems: the two systems are the GCRS
53 * and the CIP,CIO, and the point is the ascending node of the
54 * CIP equator. The quantity s remains below 0.1 arcsecond
55 * throughout 1900-2100.
56 *
57 * 3) The series used to compute s is in fact for s+xy/2, where x and y
58 * are the x and y components of the CIP unit vector; this series is
59 * more compact than a direct series for s would be. This routine
60 * requires X,Y to be supplied by the caller, who is responsible for
61 * providing values that are consistent with the supplied date.
62 *
63 * 4) The model is consistent with the "P03" precession (Capitaine et
64 * al. 2003), adopted by IAU 2006 Resolution 1, 2006, and the
65 * IAU 2000A nutation (with P03 adjustments).
66 *
67 * Called:
68 *
69 * gal_fal03 mean anomaly of the Moon
70 * gal_falp03 mean anomaly of the Sun
71 * gal_faf03 mean argument of the latitude of the Moon
72 * gal_fad03 mean elongation of the Moon from the Sun
73 * gal_faom03 mean longitude of the Moon's ascending node
74 * gal_fave03 mean longitude of Venus
75 * gal_fae03 mean longitude of Earth
76 * gal_fapa03 general accumulated precession in longitude
77 *
78 * References:
79 *
80 * Capitaine, N., Wallace, P.T. & Chapront, J., 2003, Astron.
81 * Astrophys. 432, 355
82 *
83 * McCarthy, D.D., Petit, G. (eds.) 2004, IERS Conventions (2003),
84 * IERS Technical Note No. 32, BKG
85 *
86 * This revision:
87 *
88 * 2007 May 29 ( c version 2008 January 19 )
89 *
90 *
91 * Copyright (C) 2008 Paul C. L. Willmott. See notes at end.
92 *
93 *-----------------------------------------------------------------------
94 */
95
96 #include <math.h>
97 #include "gal_const.h"
98 #include "gal_s06.h"
99 #include "gal_fal03.h"
100 #include "gal_falp03.h"
101 #include "gal_faf03.h"
102 #include "gal_fad03.h"
103 #include "gal_faom03.h"
104 #include "gal_fave03.h"
105 #include "gal_fae03.h"
106 #include "gal_fapa03.h"
107
108 double
gal_s06(double date1,double date2,double x,double y)109 gal_s06
110 (
111 double date1,
112 double date2,
113 double x,
114 double y
115 )
116 {
117
118 /*
119 * Time since J2000, in Julian centuries
120 */
121
122 double t ;
123
124 /*
125 * Miscellaneous
126 */
127
128 int i, j ;
129 double a, s0, s1, s2, s3, s4, s5 ;
130
131 /*
132 * Fundamental arguments
133 */
134
135 double fa[8] ;
136
137 /*
138 * ---------------------
139 * The series for s+XY/2
140 * ---------------------
141 */
142
143 /*
144 * Number of terms in the series
145 */
146
147 #define NSP 6
148 #define NS0 33
149 #define NS1 3
150 #define NS2 25
151 #define NS3 4
152 #define NS4 1
153
154 /*
155 * Polynomial coefficients
156 */
157
158 static const double SP[NSP] = {
159 94e-6 ,
160 3808.65e-6 ,
161 -122.68e-6 ,
162 -72574.11e-6 ,
163 27.98e-6 ,
164 15.62e-6 ,
165 } ;
166
167 /*
168 * Coefficients of l,l',F,D,Om,LVe,LE,pA
169 * Argument coefficients for t^0
170 */
171
172 static const int KS0[NS0][8] = {
173
174 { 0, 0, 0, 0, 1, 0, 0, 0, } ,
175 { 0, 0, 0, 0, 2, 0, 0, 0, } ,
176 { 0, 0, 2, -2, 3, 0, 0, 0, } ,
177 { 0, 0, 2, -2, 1, 0, 0, 0, } ,
178 { 0, 0, 2, -2, 2, 0, 0, 0, } ,
179 { 0, 0, 2, 0, 3, 0, 0, 0, } ,
180 { 0, 0, 2, 0, 1, 0, 0, 0, } ,
181 { 0, 0, 0, 0, 3, 0, 0, 0, } ,
182 { 0, 1, 0, 0, 1, 0, 0, 0, } ,
183 { 0, 1, 0, 0, -1, 0, 0, 0, } ,
184 { 1, 0, 0, 0, -1, 0, 0, 0, } ,
185 { 1, 0, 0, 0, 1, 0, 0, 0, } ,
186 { 0, 1, 2, -2, 3, 0, 0, 0, } ,
187 { 0, 1, 2, -2, 1, 0, 0, 0, } ,
188 { 0, 0, 4, -4, 4, 0, 0, 0, } ,
189 { 0, 0, 1, -1, 1, -8, 12, 0, } ,
190 { 0, 0, 2, 0, 0, 0, 0, 0, } ,
191 { 0, 0, 2, 0, 2, 0, 0, 0, } ,
192 { 1, 0, 2, 0, 3, 0, 0, 0, } ,
193 { 1, 0, 2, 0, 1, 0, 0, 0, } ,
194 { 0, 0, 2, -2, 0, 0, 0, 0, } ,
195 { 0, 1, -2, 2, -3, 0, 0, 0, } ,
196 { 0, 1, -2, 2, -1, 0, 0, 0, } ,
197 { 0, 0, 0, 0, 0, 8,-13, -1, } ,
198 { 0, 0, 0, 2, 0, 0, 0, 0, } ,
199 { 2, 0, -2, 0, -1, 0, 0, 0, } ,
200 { 0, 1, 2, -2, 2, 0, 0, 0, } ,
201 { 1, 0, 0, -2, 1, 0, 0, 0, } ,
202 { 1, 0, 0, -2, -1, 0, 0, 0, } ,
203 { 0, 0, 4, -2, 4, 0, 0, 0, } ,
204 { 0, 0, 2, -2, 4, 0, 0, 0, } ,
205 { 1, 0, -2, 0, -3, 0, 0, 0, } ,
206 { 1, 0, -2, 0, -1, 0, 0, 0, } ,
207
208 } ;
209
210 /*
211 * Argument coefficients for t^1
212 */
213
214 static const int KS1[NS1][8] = {
215
216 { 0, 0, 0, 0, 2, 0, 0, 0, } ,
217 { 0, 0, 0, 0, 1, 0, 0, 0, } ,
218 { 0, 0, 2, -2, 3, 0, 0, 0, } ,
219
220 } ;
221
222 /*
223 * Argument coefficients for t^2
224 */
225
226 static const int KS2[NS2][8] = {
227
228 { 0, 0, 0, 0, 1, 0, 0, 0, } ,
229 { 0, 0, 2, -2, 2, 0, 0, 0, } ,
230 { 0, 0, 2, 0, 2, 0, 0, 0, } ,
231 { 0, 0, 0, 0, 2, 0, 0, 0, } ,
232 { 0, 1, 0, 0, 0, 0, 0, 0, } ,
233 { 1, 0, 0, 0, 0, 0, 0, 0, } ,
234 { 0, 1, 2, -2, 2, 0, 0, 0, } ,
235 { 0, 0, 2, 0, 1, 0, 0, 0, } ,
236 { 1, 0, 2, 0, 2, 0, 0, 0, } ,
237 { 0, 1, -2, 2, -2, 0, 0, 0, } ,
238 { 1, 0, 0, -2, 0, 0, 0, 0, } ,
239 { 0, 0, 2, -2, 1, 0, 0, 0, } ,
240 { 1, 0, -2, 0, -2, 0, 0, 0, } ,
241 { 0, 0, 0, 2, 0, 0, 0, 0, } ,
242 { 1, 0, 0, 0, 1, 0, 0, 0, } ,
243 { 1, 0, -2, -2, -2, 0, 0, 0, } ,
244 { 1, 0, 0, 0, -1, 0, 0, 0, } ,
245 { 1, 0, 2, 0, 1, 0, 0, 0, } ,
246 { 2, 0, 0, -2, 0, 0, 0, 0, } ,
247 { 2, 0, -2, 0, -1, 0, 0, 0, } ,
248 { 0, 0, 2, 2, 2, 0, 0, 0, } ,
249 { 2, 0, 2, 0, 2, 0, 0, 0, } ,
250 { 2, 0, 0, 0, 0, 0, 0, 0, } ,
251 { 1, 0, 2, -2, 2, 0, 0, 0, } ,
252 { 0, 0, 2, 0, 0, 0, 0, 0, } ,
253
254 } ;
255
256 /*
257 * Argument coefficients for t^3
258 */
259
260 static const int KS3[NS3][8] = {
261
262 { 0, 0, 0, 0, 1, 0, 0, 0, } ,
263 { 0, 0, 2, -2, 2, 0, 0, 0, } ,
264 { 0, 0, 2, 0, 2, 0, 0, 0, } ,
265 { 0, 0, 0, 0, 2, 0, 0, 0, } ,
266
267 } ;
268
269 /*
270 * Argument coefficients for t^4
271 */
272
273 static const int KS4[NS4][8] = {
274
275 { 0, 0, 0, 0, 1, 0, 0, 0, } ,
276
277 } ;
278
279 /*
280 * Sine and cosine coefficients
281 */
282
283 /*
284 * Sine ane cosine coefficients for t^0
285 */
286
287 static const double SS0[NS0][2] = {
288
289 { -2640.73e-6, +0.39e-6, } ,
290 { -63.53e-6, +0.02e-6, } ,
291 { -11.75e-6, -0.01e-6, } ,
292 { -11.21e-6, -0.01e-6, } ,
293 { +4.57e-6, 0.00e-6, } ,
294 { -2.02e-6, 0.00e-6, } ,
295 { -1.98e-6, 0.00e-6, } ,
296 { +1.72e-6, 0.00e-6, } ,
297 { +1.41e-6, +0.01e-6, } ,
298 { +1.26e-6, +0.01e-6, } ,
299 { +0.63e-6, 0.00e-6, } ,
300 { +0.63e-6, 0.00e-6, } ,
301 { -0.46e-6, 0.00e-6, } ,
302 { -0.45e-6, 0.00e-6, } ,
303 { -0.36e-6, 0.00e-6, } ,
304 { +0.24e-6, +0.12e-6, } ,
305 { -0.32e-6, 0.00e-6, } ,
306 { -0.28e-6, 0.00e-6, } ,
307 { -0.27e-6, 0.00e-6, } ,
308 { -0.26e-6, 0.00e-6, } ,
309 { +0.21e-6, 0.00e-6, } ,
310 { -0.19e-6, 0.00e-6, } ,
311 { -0.18e-6, 0.00e-6, } ,
312 { +0.10e-6, -0.05e-6, } ,
313 { -0.15e-6, 0.00e-6, } ,
314 { +0.14e-6, 0.00e-6, } ,
315 { +0.14e-6, 0.00e-6, } ,
316 { -0.14e-6, 0.00e-6, } ,
317 { -0.14e-6, 0.00e-6, } ,
318 { -0.13e-6, 0.00e-6, } ,
319 { +0.11e-6, 0.00e-6, } ,
320 { -0.11e-6, 0.00e-6, } ,
321 { -0.11e-6, 0.00e-6, } ,
322
323 } ;
324
325 /*
326 * Sine ane cosine coefficients for t^1
327 */
328
329 static const double SS1[NS1][2] = {
330
331 { -0.07e-6, +3.57e-6, } ,
332 { +1.73e-6, -0.03e-6, } ,
333 { 0.00e-6, +0.48e-6, } ,
334
335 } ;
336
337 /*
338 * Sine ane cosine coefficients for t^2
339 */
340
341 static const double SS2[NS2][2] = {
342
343 { +743.52e-6, -0.17e-6, } ,
344 { 56.91e-6, +0.06e-6, } ,
345 { +9.84e-6, -0.01e-6, } ,
346 { -8.85e-6, +0.01e-6, } ,
347 { -6.38e-6, -0.05e-6, } ,
348 { -3.07e-6, 0.00e-6, } ,
349 { +2.23e-6, 0.00e-6, } ,
350 { +1.67e-6, 0.00e-6, } ,
351 { +1.30e-6, 0.00e-6, } ,
352 { +0.93e-6, 0.00e-6, } ,
353 { +0.68e-6, 0.00e-6, } ,
354 { -0.55e-6, 0.00e-6, } ,
355 { +0.53e-6, 0.00e-6, } ,
356 { -0.27e-6, 0.00e-6, } ,
357 { -0.27e-6, 0.00e-6, } ,
358 { -0.26e-6, 0.00e-6, } ,
359 { -0.25e-6, 0.00e-6, } ,
360 { +0.22e-6, 0.00e-6, } ,
361 { -0.21e-6, 0.00e-6, } ,
362 { +0.20e-6, 0.00e-6, } ,
363 { +0.17e-6, 0.00e-6, } ,
364 { +0.13e-6, 0.00e-6, } ,
365 { -0.13e-6, 0.00e-6, } ,
366 { -0.12e-6, 0.00e-6, } ,
367 { -0.11e-6, 0.00e-6, } ,
368
369 } ;
370
371 /*
372 * Sine ane cosine coefficients for t^3
373 */
374
375 static const double SS3[NS3][2] = {
376
377 { +0.30e-6, -23.42e-6, } ,
378 { -0.03e-6, -1.46e-6, } ,
379 { -0.01e-6, -0.25e-6, } ,
380 { 0.00e-6, +0.23e-6, } ,
381
382 } ;
383
384 /*
385 * Sine ane cosine coefficients for t^4
386 */
387
388 static const double SS4[NS4][2] = {
389
390 { -0.26e-6, -0.01e-6, } ,
391
392 } ;
393
394 /*
395 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
396 */
397
398 /*
399 * Interval between fundamental epoch J2000.0 and current date (JC).
400 */
401
402 t = ( ( date1 - GAL_J2000 ) + date2 ) / GAL_DJC ;
403
404 /*
405 * Fundamental Arguments (from IERS Conventions 2003)
406 */
407
408 /*
409 * Mean anomaly of the Moon.
410 */
411
412 fa[0] = gal_fal03 ( t ) ;
413
414 /*
415 * Mean anomaly of the Sun.
416 */
417
418 fa[1] = gal_falp03 ( t ) ;
419
420 /*
421 * Mean longitude of the Moon minus that of the ascending node.
422 */
423
424 fa[2] = gal_faf03 ( t ) ;
425
426 /*
427 * Mean elongation of the Moon from the Sun.
428 */
429
430 fa[3] = gal_fad03 ( t ) ;
431
432 /*
433 * Mean longitude of the ascending node of the Moon.
434 */
435
436 fa[4] = gal_faom03 ( t ) ;
437
438 /*
439 * Mean longitude of Venus.
440 */
441
442 fa[5] = gal_fave03 ( t ) ;
443
444 /*
445 * Mean longitude of Earth.
446 */
447
448 fa[6] = gal_fae03 ( t ) ;
449
450 /*
451 * General precession in longitude.
452 */
453
454 fa[7] = gal_fapa03 ( t ) ;
455
456 /*
457 * Evaluate S.
458 */
459
460 s0 = SP[0] ;
461 s1 = SP[1] ;
462 s2 = SP[2] ;
463 s3 = SP[3] ;
464 s4 = SP[4] ;
465 s5 = SP[5] ;
466
467 for ( i = NS0 - 1; i >= 0; i-- ) {
468 a = 0.0 ;
469 for ( j = 0; j < 8; j++ ) {
470 a += ( double )(KS0[i][j]) * fa[j] ;
471 }
472 s0 += ( SS0[i][0] * sin ( a ) + SS0[i][1] * cos ( a ) ) ;
473 }
474
475 for ( i = NS1 - 1; i >= 0; i-- ) {
476 a = 0.0 ;
477 for ( j = 0; j < 8; j++ ) {
478 a += ( double )(KS1[i][j]) * fa[j] ;
479 }
480 s1 += ( SS1[i][0] * sin ( a ) + SS1[i][1] * cos ( a ) ) ;
481 }
482
483 for ( i = NS2 - 1; i >= 0; i-- ) {
484 a = 0.0 ;
485 for ( j = 0; j < 8; j++ ) {
486 a += ( double )(KS2[i][j]) * fa[j] ;
487 }
488 s2 += ( SS2[i][0] * sin ( a ) + SS2[i][1] * cos ( a ) ) ;
489 }
490
491 for ( i = NS3 - 1; i >= 0; i-- ) {
492 a = 0.0 ;
493 for ( j = 0; j < 8; j++ ) {
494 a += ( double )(KS3[i][j]) * fa[j] ;
495 }
496 s3 += ( SS3[i][0] * sin ( a ) + SS3[i][1] * cos ( a ) ) ;
497 }
498
499 for ( i = NS4 - 1; i >= 0; i-- ) {
500 a = 0.0 ;
501 for ( j = 0; j < 8; j++ ) {
502 a += ( double )(KS4[i][j]) * fa[j] ;
503 }
504 s4 += ( SS4[i][0] * sin ( a ) + SS4[i][1] * cos ( a ) ) ;
505 }
506
507 return ( s0 +
508 ( s1 +
509 ( s2 +
510 ( s3 +
511 ( s4 +
512 s5 * t ) * t ) * t ) * t ) * t ) * GAL_AS2R - x * y / 2.0 ;
513
514 /*
515 * Finished.
516 */
517
518 }
519
520 /*
521 * gal - General Astrodynamics Library
522 * Copyright (C) 2008 Paul C. L. Willmott
523 *
524 * This program is free software; you can redistribute it and/or modify
525 * it under the terms of the GNU General Public License as published by
526 * the Free Software Foundation; either version 2 of the License, or
527 * (at your option) any later version.
528 *
529 * This program is distributed in the hope that it will be useful,
530 * but WITHOUT ANY WARRANTY; without even the implied warranty of
531 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
532 * GNU General Public License for more details.
533 *
534 * You should have received a copy of the GNU General Public License along
535 * with this program; if not, write to the Free Software Foundation, Inc.,
536 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
537 *
538 * Contact:
539 *
540 * Paul Willmott
541 * vp9mu@amsat.org
542 */
543