SUBROUTINE NUTATE( T, DPSI, DEPS, DDPSI, DDEPS ) C C GIVEN DYNAMIC TIME T (JULIAN MILLENIA SINCE J2000), THIS ROUTINE C CALCULATES THE NUTATION IN LONGITUDE, DPSI, AND IN THE OBLIQUITY C OF THE ECLIPTIC, DEPS, AS WELL AS THE ASSOCIATED VELOCITY CORREC- C TIONS, DDPSI AND DDEPS. DPSI AND DEPS ARE RETURNED IN RADIANS, C AND DDPSI AND DDEPS IN RADIANS/JULIAN MILLENIA. C C REFERENCE: P. K. SEIDELMANN, 1980 IAU THEORY OF NUTATION, C CELESTIAL MECHANICS V. 27 (1982), PP. 79-106 C C B. KNAPP, 1999-06-11, 2000-09-01, 2001-04-24 (i*4) C C C RCS DATA C C $Header$ C C $Log$ C C IMPLICIT NONE C C INPUT REAL*8 T C C OUTPUT REAL*8 DPSI, DEPS, DDPSI, DDEPS C C CONSTANTS REAL*8 D2R, S2R PARAMETER (D2R = 3.14159265358979324D0/180.D0) PARAMETER (S2R = 1.D-4/3600.D0*D2R) C (CONVERT UNITS OF 10000 ARC SECONDS TO RADIANS) C C LOCAL VARIABLES REAL*8 M, E, F, D, O, DM, DE, DF, DD, DO, L, DL, & SL, CL, PSI_K, EPS_K INTEGER*4 K C C NUTATION COEFFICIENTS (CONSTANTS) INTEGER*4 NTERMS PARAMETER (NTERMS=106) C REAL*8 CM(NTERMS) / & 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, -1, 1, -1, -1, & 0, 1, -2, 2, 0, 2, 1, 2, 0, -1, 0, -1, 0, 0, 0, 1, & 0, -1, 2, 1, 0, 0, 0, 1, -2, 2, 0, 1, 1, 2, 0, 0, & 1, 2, 0, 0, 1, 0, 1, 3, 1, -1, 0, -2, 1, 1, 1, 0, & 1, -2, -1, 2, 0, -1, 1, 3, -2, 1, 2, -1, 1, -2, 0, 0, & 0, 2, -1, 0, 1, 2, 1, 0, 3, 1, 0, -1, 0, 0, 0, 1, & 0, 1, 2, -1, 2, 1, 1, 2, 0, 0 & / REAL*8 CE(NTERMS) / & 0, 0, 0, 0, 1, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 0, & -1, 0, 0, 0, 1, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, & -1, 0, 1, 1, 0, 0, 0, 0, -1, -1, -1, 0, -1, 0, 1, 0, & 0, 0, 0, 0, -2, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, & 0, 1, 0, 1, 0, 0, -1, 0, 0, 0, 1, -1, 0, 0, 1, 0, & -1, 1, 0, 0, 0, 0, 0, 0, 0, 1 & / REAL*8 CF(NTERMS) / & 0, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 2, 2, 0, 0, 2, & 0, 2, 2, 0, 2, 2, 2, 0, 2, 2, 2, 0, 0, 2, 0, 0, & 0, 2, -2, 2, 2, 2, 2, 0, 0, 2, 0, 2, 0, 2, 0, 2, & 0, 0, 2, 0, -2, 0, 0, 2, 2, 2, 2, 2, 0, 2, 0, 2, & 0, 0, 2, 0, 2, 2, 2, 0, 2, 2, 2, 0, 2, 2, 0, -2, & -2, 0, 0, 2, 0, 2, 0, 4, 2, 2, 2, 0, -2, 2, 0, -2, & 2, 0, -2, 4, 0, 0, -2, 0, 2, 0 & / REAL*8 CD(NTERMS) / & 0, -2, 0, 0, 0, 0, -2, 0, 0, -2, -2, -2, 0, 0, 0, 2, & 2, 0, 0, -2, 2, 0, -2, 0, 0, 0, -2, 2, 0, -2, 0, -2, & 0, 2, 0, 2, 0, 0, 2, -2, 2, -2, 2, -2, 2, 0, -2, -2, & 0, -2, -2, -2, 0, 1, -1, 0, 0, 2, 2, 0, -1, 0, 0, 1, & 0, 0, 4, 0, -2, -2, 0, 0, 2, 2, -2, 0, -2, 4, 0, 2, & 2, -2, 1, -2, 2, 2, -2, -2, -2, -2, 0, 2, 0, -1, 2, -2, & 0, -2, 0, 0, -4, -4, 2, 2, 4, 1 & / REAL*8 CO(NTERMS) / & 1, 2, 2, 2, 0, 0, 2, 1, 2, 2, 0, 1, 2, 1, 1, 2, & 0, 1, 1, 0, 2, 2, 2, 0, 0, 1, 0, 1, 1, 2, 0, 1, & 1, 1, 0, 2, 2, 2, 1, 0, 1, 2, 1, 1, 0, 1, 1, 1, & 0, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 2, 0, 0, 0, 2, & 2, 1, 2, 1, 1, 1, 2, 0, 2, 1, 1, 2, 2, 2, 2, 0, & 1, 0, 1, 0, 1, 2, 0, 2, 2, 0, 1, 1, 1, 2, 0, 0, & 1, 1, 1, 2, 0, 0, 0, 0, 2, 0 & / REAL*8 PSI0(NTERMS) / & -171996, -13187, -2274, 2062, 1426, 712, -517, -386, & -301, 217, -158, 129, 123, 63, -58, -59, & 63, -51, 46, 48, -38, -31, 29, 29, & 26, 21, -22, 16, -15, -16, 17, -13, & -12, -10, 11, -8, 7, -7, -7, -7, & -6, 6, -6, 6, 6, -5, -5, -5, & 5, 4, 4, -4, 4, -4, -4, -3, & -3, -3, -3, -3, -3, 3, -3, 2, & -2, -2, -2, 2, -2, -2, 2, 2, & 1, -1, 1, 1, 1, -1, 1, -1, & 1, 1, 1, -1, -1, -1, 1, 1, & 1, -1, 1, 1, -1, -1, -1, -1, & -1, -1, 1, 1, -1, -1, -1, 1, & -1, 1 & / REAL*8 PSI1(NTERMS) / & -1742, -16, -2, 2, -34, 1, 12, -4, & 0, -5, 0, 1, 0, 1, 1, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 1, -1, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0 & / REAL*8 EPS0(NTERMS) / & 92025, 5736, 977, -895, 54, -7, 224, 200, & 129, -95, -1, -70, -53, -33, 32, 26, & -2, 27, -24, 1, 16, 13, -12, -1, & -1, -10, 0, -8, 9, 7, 0, 7, & 6, 5, 0, 3, -3, 3, 3, 0, & 3, -3, 3, -3, 0, 3, 3, 3, & 0, -2, -2, 0, 0, 0, 0, 1, & 1, 1, 1, 1, 0, 0, 0, -1, & 1, 1, 1, -1, 1, 1, -1, 0, & -1, 1, -1, -1, -1, 1, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0 & / REAL*8 EPS1(NTERMS) / & 89, -31, -5, 5, -1, 0, -6, 0, & -1, 3, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0 & / C C POLYNOMIALS (IN DEGREES) FOR THE FUNDAMENTAL ARGUMENTS (CONSTANTS) INTEGER*4 PDEG PARAMETER (PDEG=3) C C MEAN ANOMALY OF THE MOON REAL*8 PM(0:PDEG) / 1.349629813888889D+2, 4.771988673980556D+6, & 8.697222222222222D-1, 1.777777777777778D-2 / C C MEAN ANOMALY OF THE SUN (EARTH) REAL*8 PE(0:PDEG) / 3.575277233333333D+2, 3.599905034000000D+5, & -1.602777777777778D-2, -3.333333333333333D-3 / C C MOON'S ARGUMENT OF LATITUDE REAL*8 PF(0:PDEG) / 9.327191027777778D+1, 4.832020175380556D+6, & -3.682500000000000D-1, 3.055555555555556D-3 / C C MEAN ELONGATION OF THE MOON FROM THE SUN REAL*8 PD(0:PDEG) / 2.978503630555556D+2, 4.452671114800000D+6, & -1.914166666666667D-1, 5.277777777777778D-3 / C C LONGITUDE OF THE ASCENDING NODE OF THE MOON'S MEAN ORBIT ON THE C ECLIPTIC, MEASURED FROM THE MEAN EQUINOX OF THE DATE REAL*8 PO(0:PDEG) / 1.250445222222222D+2, -1.934136260833333D+4, & 2.070833333333333D-1, 2.222222222222222D-3 / C SAVE CM, CE, CF, CD, CO, PSI0, PSI1, EPS0, EPS1, & PM, PE, PF, PD, PO C C EVALUATE THE FUNDAMENTAL ARGUMENTS M = PM(0)+T*(PM(1)+T*( PM(2)+T* PM(3))) DM = PM(1)+T*(2.D0*PM(2)+T*3.D0*PM(3)) E = PE(0)+T*(PE(1)+T*( PE(2)+T* PE(3))) DE = PE(1)+T*(2.D0*PE(2)+T*3.D0*PE(3)) F = PF(0)+T*(PF(1)+T*( PF(2)+T* PF(3))) DF = PF(1)+T*(2.D0*PF(2)+T*3.D0*PF(3)) D = PD(0)+T*(PD(1)+T*( PD(2)+T* PD(3))) DD = PD(1)+T*(2.D0*PD(2)+T*3.D0*PD(3)) O = PO(0)+T*(PO(1)+T*( PO(2)+T* PO(3))) DO = PO(1)+T*(2.D0*PO(2)+T*3.D0*PO(3)) C C COMPUTE THE SUM OF ALL TERMS DPSI = 0.D0 DEPS = 0.D0 DDPSI = 0.D0 DDEPS = 0.D0 DO 1 K=NTERMS,1,-1 L = MOD( CM(K)* M + CE(K)* E + CF(K)* F + CD(K)* D + CO(K)* O, & 360.D0 )*D2R DL = ( CM(K)*DM + CE(K)*DE + CF(K)*DF + CD(K)*DD + CO(K)*DO & )*D2R SL = SIN( L ) CL = COS( L ) PSI_K = PSI0(K) + PSI1(K)*T EPS_K = EPS0(K) + EPS1(K)*T DPSI = DPSI + PSI_K*SL DEPS = DEPS + EPS_K*CL DDPSI = DDPSI + PSI_K*CL*DL + SL*PSI1(K) DDEPS = DDEPS - EPS_K*SL*DL + CL*EPS1(K) 1 CONTINUE C DPSI = DPSI*S2R DEPS = DEPS*S2R DDPSI = DDPSI*S2R DDEPS = DDEPS*S2R C RETURN END