/* Nutation in longitude and obliquity * computed at Julian date J. * * References: * "Summary of 1980 IAU Theory of Nutation (Final Report of the * IAU Working Group on Nutation)", P. K. Seidelmann et al., in * Transactions of the IAU Vol. XVIII A, Reports on Astronomy, * P. A. Wayman, ed.; D. Reidel Pub. Co., 1982. * * "Nutation and the Earth's Rotation", * I.A.U. Symposium No. 78, May, 1977, page 256. * I.A.U., 1980. * * Woolard, E.W., "A redevelopment of the theory of nutation", * The Astronomical Journal, 58, 1-3 (1953). * * This program implements all of the 1980 IAU nutation series. * Results checked at 100 points against the 1986 AA; all agreed. * * * - S. L. Moshier, November 1987 */ /* The answers are posted here by nutlo(): */ double jdnut = -1.0; /* time to which the nutation applies */ double nutl = 0.0; /* nutation in longitude (radians) */ double nuto = 0.0; /* nutation in obliquity (radians) */ extern double eps, coseps, sineps, DTR; /* Each term in the expansion has a trigonometric * argument given by * W = i*MM + j*MS + k*FF + l*DD + m*OM * where the variables are defined below. * The nutation in longitude is a sum of terms of the * form (a + bT) * sin(W). The terms for nutation in obliquity * are of the form (c + dT) * cos(W). The coefficients * are arranged in the tabulation as follows: * * Coefficient: * i j k l m a b c d * 0, 0, 0, 0, 1, -171996, -1742, 92025, 89, * The first line of the table, above, is done separately * since two of the values do not fit into 16 bit integers. * The values a and c are arc seconds times 10000. b and d * are arc seconds per Julian century times 100000. i through m * are integers. See the program for interpretation of MM, MS, * etc., which are mean orbital elements of the Sun and Moon. * * If terms with coefficient less than X are omitted, the peak * errors will be: * * omit error, omit error, * a < longitude c < obliquity * .0005" .0100" .0008" .0094" * .0046 .0492 .0095 .0481 * .0123 .0880 .0224 .0905 * .0386 .1808 .0895 .1129 */ short nt[105*9] = { 0, 0, 0, 0, 2, 2062, 2,-895, 5, -2, 0, 2, 0, 1, 46, 0,-24, 0, 2, 0,-2, 0, 0, 11, 0, 0, 0, -2, 0, 2, 0, 2,-3, 0, 1, 0, 1,-1, 0,-1, 0,-3, 0, 0, 0, 0,-2, 2,-2, 1,-2, 0, 1, 0, 2, 0,-2, 0, 1, 1, 0, 0, 0, 0, 0, 2,-2, 2,-13187,-16, 5736,-31, 0, 1, 0, 0, 0, 1426,-34, 54,-1, 0, 1, 2,-2, 2,-517, 12, 224,-6, 0,-1, 2,-2, 2, 217,-5,-95, 3, 0, 0, 2,-2, 1, 129, 1,-70, 0, 2, 0, 0,-2, 0, 48, 0, 1, 0, 0, 0, 2,-2, 0,-22, 0, 0, 0, 0, 2, 0, 0, 0, 17,-1, 0, 0, 0, 1, 0, 0, 1,-15, 0, 9, 0, 0, 2, 2,-2, 2,-16, 1, 7, 0, 0,-1, 0, 0, 1,-12, 0, 6, 0, -2, 0, 0, 2, 1,-6, 0, 3, 0, 0,-1, 2,-2, 1,-5, 0, 3, 0, 2, 0, 0,-2, 1, 4, 0,-2, 0, 0, 1, 2,-2, 1, 4, 0,-2, 0, 1, 0, 0,-1, 0,-4, 0, 0, 0, 2, 1, 0,-2, 0, 1, 0, 0, 0, 0, 0,-2, 2, 1, 1, 0, 0, 0, 0, 1,-2, 2, 0,-1, 0, 0, 0, 0, 1, 0, 0, 2, 1, 0, 0, 0, -1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 2,-2, 0,-1, 0, 0, 0, 0, 0, 2, 0, 2,-2274,-2, 977,-5, 1, 0, 0, 0, 0, 712, 1,-7, 0, 0, 0, 2, 0, 1,-386,-4, 200, 0, 1, 0, 2, 0, 2,-301, 0, 129,-1, 1, 0, 0,-2, 0,-158, 0,-1, 0, -1, 0, 2, 0, 2, 123, 0,-53, 0, 0, 0, 0, 2, 0, 63, 0,-2, 0, 1, 0, 0, 0, 1, 63, 1,-33, 0, -1, 0, 0, 0, 1,-58,-1, 32, 0, -1, 0, 2, 2, 2,-59, 0, 26, 0, 1, 0, 2, 0, 1,-51, 0, 27, 0, 0, 0, 2, 2, 2,-38, 0, 16, 0, 2, 0, 0, 0, 0, 29, 0,-1, 0, 1, 0, 2,-2, 2, 29, 0,-12, 0, 2, 0, 2, 0, 2,-31, 0, 13, 0, 0, 0, 2, 0, 0, 26, 0,-1, 0, -1, 0, 2, 0, 1, 21, 0,-10, 0, -1, 0, 0, 2, 1, 16, 0,-8, 0, 1, 0, 0,-2, 1,-13, 0, 7, 0, -1, 0, 2, 2, 1,-10, 0, 5, 0, 1, 1, 0,-2, 0,-7, 0, 0, 0, 0, 1, 2, 0, 2, 7, 0,-3, 0, 0,-1, 2, 0, 2,-7, 0, 3, 0, 1, 0, 2, 2, 2,-8, 0, 3, 0, 1, 0, 0, 2, 0, 6, 0, 0, 0, 2, 0, 2,-2, 2, 6, 0,-3, 0, 0, 0, 0, 2, 1,-6, 0, 3, 0, 0, 0, 2, 2, 1,-7, 0, 3, 0, 1, 0, 2,-2, 1, 6, 0,-3, 0, 0, 0, 0,-2, 1,-5, 0, 3, 0, 1,-1, 0, 0, 0, 5, 0, 0, 0, 2, 0, 2, 0, 1,-5, 0, 3, 0, 0, 1, 0,-2, 0,-4, 0, 0, 0, 1, 0,-2, 0, 0, 4, 0, 0, 0, 0, 0, 0, 1, 0,-4, 0, 0, 0, 1, 1, 0, 0, 0,-3, 0, 0, 0, 1, 0, 2, 0, 0, 3, 0, 0, 0, 1,-1, 2, 0, 2,-3, 0, 1, 0, -1,-1, 2, 2, 2,-3, 0, 1, 0, -2, 0, 0, 0, 1,-2, 0, 1, 0, 3, 0, 2, 0, 2,-3, 0, 1, 0, 0,-1, 2, 2, 2,-3, 0, 1, 0, 1, 1, 2, 0, 2, 2, 0,-1, 0, -1, 0, 2,-2, 1,-2, 0, 1, 0, 2, 0, 0, 0, 1, 2, 0,-1, 0, 1, 0, 0, 0, 2,-2, 0, 1, 0, 3, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 2, 1, 2, 2, 0,-1, 0, -1, 0, 0, 0, 2, 1, 0,-1, 0, 1, 0, 0,-4, 0,-1, 0, 0, 0, -2, 0, 2, 2, 2, 1, 0,-1, 0, -1, 0, 2, 4, 2,-2, 0, 1, 0, 2, 0, 0,-4, 0,-1, 0, 0, 0, 1, 1, 2,-2, 2, 1, 0,-1, 0, 1, 0, 2, 2, 1,-1, 0, 1, 0, -2, 0, 2, 4, 2,-1, 0, 1, 0, -1, 0, 4, 0, 2, 1, 0, 0, 0, 1,-1, 0,-2, 0, 1, 0, 0, 0, 2, 0, 2,-2, 1, 1, 0,-1, 0, 2, 0, 2, 2, 2,-1, 0, 0, 0, 1, 0, 0, 2, 1,-1, 0, 0, 0, 0, 0, 4,-2, 2, 1, 0, 0, 0, 3, 0, 2,-2, 2, 1, 0, 0, 0, 1, 0, 2,-2, 0,-1, 0, 0, 0, 0, 1, 2, 0, 1, 1, 0, 0, 0, -1,-1, 0, 2, 1, 1, 0, 0, 0, 0, 0,-2, 0, 1,-1, 0, 0, 0, 0, 0, 2,-1, 2,-1, 0, 0, 0, 0, 1, 0, 2, 0,-1, 0, 0, 0, 1, 0,-2,-2, 0,-1, 0, 0, 0, 0,-1, 2, 0, 1,-1, 0, 0, 0, 1, 1, 0,-2, 1,-1, 0, 0, 0, 1, 0,-2, 2, 0,-1, 0, 0, 0, 2, 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 2, 4, 2,-1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, }; /* arrays to hold sines and cosines of multiple angles */ double ss[5][8] = {0.0}; double cc[5][8] = {0.0}; extern double ss[5][8], cc[5][8]; nutlo(J) double J; { double f, g, w, T; double MM, MS, FF, DD, OM; double cu, su, cv, sv, cw, sw; double C, D; int i, j, k, k1, m; short *p; double sin(), cos(), mod360(); if( jdnut == J ) return(0); jdnut = J; /* Julian centuries from 2000 January 1.5, * barycentric dynamical time */ T = (J-2451545.0)/36525.0; /* Fundamental arguments in the FK5 reference system. /* The coefficients, originally given to 0.001", * are converted here to degrees. */ /* longitude of the mean ascending node of the lunar orbit * on the ecliptic, measured from the mean equinox of date */ OM = ((2.22e-6*T + 0.00207833)*T - 1934.1362608)*T + 125.0445222; OM = DTR * mod360(OM); /* mean longitude of the Sun minus the * mean longitude of the Sun's perigee */ MS = ((-3.33e-6*T - 0.0001603)*T + 35999.0503400)*T + 357.5277233; MS = DTR * mod360(MS); /* mean longitude of the Moon minus the * mean longitude of the Moon's perigee */ MM = ((1.778e-5*T + 0.0086972)*T + 477198.8673981)*T + 134.9629814; MM = DTR * mod360(MM); /* mean longitude of the Moon minus the * mean longitude of the Moon's node */ FF = ((3.056e-6*T - 0.00368250)*T + 483202.0175381)*T + 93.2719103; FF = DTR * mod360(FF); /* mean elongation of the Moon from the Sun. */ DD = (( 5.278e-6*T - 0.0019142)*T + 445267.1114800)*T + 297.8503631; DD = DTR * mod360(DD); /* Calculate sin( i*MM ), etc. for needed multiple angles */ sscc( 0, MM, 3 ); sscc( 1, MS, 2 ); sscc( 2, FF, 4 ); sscc( 3, DD, 4 ); sscc( 4, OM, 2 ); /* first terms, not in table: */ C = (-0.01742*T - 17.1996)*ss[4][0]; /* sin(OM) */ D = ( 0.00089*T + 9.2025)*cc[4][0]; /* cos(OM) */ p = &nt[0]; /* point to start of table */ for( i=0; i<105; i++ ) { /* argument of sine and cosine */ k1 = 0; for( m=0; m<5; m++ ) { j = *p++; if( j ) { k = j; if( j < 0 ) k = -k; su = ss[m][k-1]; /* sin(k*angle) */ if( j < 0 ) su = -su; cu = cc[m][k-1]; if( k1 == 0 ) { /* set first angle */ sv = su; cv = cu; k1 = 1; } else { /* combine angles */ sw = su*cv + cu*sv; cv = cu*cv - su*sv; sv = sw; } } } /* longitude coefficient */ f = *p++ * 0.0001; if( (k = *p++) != 0 ) f += 0.00001 * T * k; /* obliquity coefficient */ g = *p++ * 0.0001; if( (k = *p++) != 0 ) g += 0.00001 * T * k; /* accumulate the terms */ C += f * sv; D += g * cv; } /* printf( "nutation: in longitude %.3f\", in obliquity %.3f\"\n", C, D ); */ /* Save answers, expressed in radians */ nutl = DTR * C / 3600.0; nuto = DTR * D / 3600.0; return(1); } /* Nutation -- AA page B20 * using nutation in longitude and obliquity from nutlo() * and obliquity of the ecliptic from epsiln() * both calculated for Julian date J. * * p[] = equatorial rectangular position vector of object for * mean ecliptic and equinox of date. */ nutate( J, p ) double J; double p[]; { double ce, se, cl, sl, sino, f; double dp[3], p1[3]; int i; double sin(), cos(); nutlo(J); /* be sure we calculated nutl and nuto */ epsiln(J); /* and also the obliquity of date */ f = eps + nuto; ce = cos( f ); se = sin( f ); sino = sin(nuto); cl = cos( nutl ); sl = sin( nutl ); /* Apply adjustment * to equatorial rectangular coordinates of object. * * This is a composite of three rotations: rotate about x axis * to ecliptic of date; rotate about new z axis by the nutation * in longitude; rotate about new x axis back to equator of date * plus nutation in obliquity. */ p1[0] = cl*p[0] - sl*coseps*p[1] - sl*sineps*p[2]; p1[1] = sl*ce*p[0] + ( cl*coseps*ce + sineps*se )*p[1] - ( sino + (1.0-cl)*sineps*se )*p[2]; p1[2] = sl*se*p[0] + ( sino + (cl-1.0)*se*coseps )*p[1] + ( cl*sineps*se + coseps*ce )*p[2]; for( i=0; i<3; i++ ) dp[i] = p1[i] - p[i]; showcor( "nutation", p, dp ); for( i=0; i<3; i++ ) p[i] = p1[i]; } /* Prepare lookup table of sin and cos ( i*Lj ) * for required multiple angles */ sscc( k, arg, n ) int k; double arg; int n; { double cu, su, cv, sv, s; int i; double sin(), cos(); su = sin(arg); cu = cos(arg); ss[k][0] = su; /* sin(L) */ cc[k][0] = cu; /* cos(L) */ sv = 2.0*su*cu; cv = cu*cu - su*su; ss[k][1] = sv; /* sin(2L) */ cc[k][1] = cv; for( i=2; i