SUBROUTINE LEA406(TDB,MODEL,P) C THIS SUBROUTINE READS THE LEA-406 LUNAR ANALYTICAL SERIES AND C GIVES THE GEOCENTRIC POSITION OF THE MOON REFERENCED TO C THE EARTH MEAN EQUATOR AND EQUINOX OF J2000. C C CALLING SEQUENCE PARAMETERS: C C TDB = INPUT 2-WORD D.P. JULIAN DATE (TDB) AT WHICH LUNAR POSITION IS WANTED. C C MODEL = INPUT INTEGER DEFINING THE MODEL OF LEA-406 ANALYTICAL SERIES C TO BE USED. THE FOLLOWING OPTIONS ARE AVAILABLE: C MODEL = 0 - THE COMPLETE ANALYTICAL SERIES LEA-406a ARE USED C (VALID OVER 1500 - 2500 OR FROM JD 2268932.5 TO JD 2634166.5). C MODEL = 1 - THE SIMPLIFIED ANALYTICAL SERIES LEA-406b ARE USED C (VALID OVER 3000BC - 3000AD OR FROM JD 625673.5 TO JD 2816787.5). C C P = OUTPUT 3-WORD D.P. ARRAY CONTAINING GEOCENTRIC POSITION OF C THE MOON REFERENCED TO THE EARTH MEAN EQUATOR AND EQUINOX OF J2000 C (THE REFERENCE FRAME OF NUMERICAL LUNAR EPHEMERIS LE-405/406). C THE UNITS ARE KM. C C EXAMPLE OF USE: C C ... C C TDB(1)=2451545.D0 C TDB(2)=1.D0 C C MODEL=0 C C CALL LEA406(TDB, MODEL, P) C C ... C C SO THE LUNAR POSITION (3-WORD D.P. ARRAY P) WILL BE CALCULATED AT EPOCH TDB(1)+TDB(2)=JD 2451546.D0 C ON THE BASE OF COMPLETE (MODEL=0) LEA-406 SERIES. C C VERSION 1.1 (AUGUST 9, 2007) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (N_REC=11200) DIMENSION TDB(2),TDB_0(2),P(3),P0(3),P1(3),PRM(3,3),IND(14), * A_R(N_REC), AT_R(N_REC), ATT_R(N_REC), * C_R(N_REC), CT_R(N_REC), CTT_R(N_REC), * F0_R(N_REC), F1_R(N_REC), F2_R(N_REC), F3_R(N_REC), F4_R(N_REC), * A_V(N_REC), AT_V(N_REC), ATT_V(N_REC), * C_V(N_REC), CT_V(N_REC), CTT_V(N_REC), * F0_V(N_REC), F1_V(N_REC), F2_V(N_REC), F3_V(N_REC), F4_V(N_REC), * A_U(N_REC), AT_U(N_REC), ATT_U(N_REC), * C_U(N_REC), CT_U(N_REC), CTT_U(N_REC), * F0_U(N_REC), F1_U(N_REC), F2_U(N_REC), F3_U(N_REC), F4_U(N_REC) LOGICAL LOADDATA COMMON /BLOCK_DATA/ PI, RADEG, RASEC, FR(0:4,14), FRM(0:4) COMMON /SAVE_DP/ TDB_0, FRM0, * A_R, AT_R, ATT_R, C_R, CT_R, CTT_R, F0_R, F1_R, F2_R, F3_R, F4_R, * A_V, AT_V, ATT_V, C_V, CT_V, CTT_V, F0_V, F1_V, F2_V, F3_V, F4_V, * A_U, AT_U, ATT_U, C_U, CT_U, CTT_U, F0_U, F1_U, F2_U, F3_U, F4_U COMMON /SAVE_INT/ MODELCUR, I_R, I_V, I_U COMMON /SAVE_LOG/ LOADDATA DATA TDB_0/2451545.D0, 0.D0/ ! INITIAL EPOCH (J2000.0) DATA MODELCUR/-1/ C C CHECK IF THE MODEL AND INPUT EPOCH ARE CONSISTENT C IF (MODEL.EQ.0) THEN IF ((TDB(1)+TDB(2)).LT.2268932.5.OR. * (TDB(1)+TDB(2)).GT.2634166.5D0) THEN WRITE (*,*) ' THE INPUT EPOCH IS OUTSIDE THE TIME INTERVAL WH *ERE THE COMPLETE LEA-406a SERIES ARE VALID. STOP!' STOP END IF ELSEIF (MODEL.EQ.1) THEN IF ((TDB(1)+TDB(2)).LT.625673.5D0.OR. * (TDB(1)+TDB(2)).GT.2816787.5D0) THEN WRITE (*,*) ' THE INPUT EPOCH IS OUTSIDE THE TIME INTERVAL WH *ERE THE SIMPLIFIED LEA-406b SERIES ARE VALID. STOP!' STOP END IF ELSE WRITE (*,*) ' MODEL=', MODEL, 'IS NOT SUPPORTED.' WRITE (*,*) ' SELECT EITHER MODEL=0 OR MODEL=1. STOP!' STOP END IF C C CALCULATION OF TIME DIFFERENCE BETWEEN THE INPUT EPOCH AND INITIAL EPOCH C DTC=((TDB(1)-TDB_0(1))+TDB(2))/36525.D0 ! TIME DIFEFRENCE IN TDB CENTURIES DTM=DTC/10.D0 ! TIME DIFEFRENCE IN TDB MILLENNIUMS C C CHECK IF THE MODEL IS CHANGED AND THE DATA ARRAYS SHOULD BE RELOADED C IF (MODEL.NE.MODELCUR) THEN LOADDATA=.TRUE. MODELCUR=MODEL ELSE LOADDATA=.FALSE. END IF C C SUBSEQUENT READING THE DATA ARRAYS IF THE MODEL IS CHANGED. C IF (LOADDATA.AND.MODEL.EQ.0) THEN ! COMPLETE ANALYTICAL SERIES LEA-406a ARE USED C C READING COMPLETE SERIES FOR COORDINATE R (LUNAR DISTANCE) C OPEN (1,FILE='table6.dat',STATUS='OLD') DO WHILE (.TRUE.) READ (1,'(I6,2X,5I3,1X,8I3,1X,I3,F16.7,2F11.6,3F19.12)',END=6) * I,IND(12),IND(11),IND(13),IND(10),IND(9),(IND(J),J=1,8),IND(14) * ,A_R(I),AT_R(I),ATT_R(I),C_R(I),CT_R(I),CTT_R(I) C_R(I)=C_R(I)*RADEG CT_R(I)=CT_R(I)*RADEG CTT_R(I)=CTT_R(I)*RADEG C C CALCULATING THE ARGUMENTS [RAD] C F0_R(I)=0.D0 F1_R(I)=0.D0 F2_R(I)=0.D0 F3_R(I)=0.D0 F4_R(I)=0.D0 DO K=1,14 F0_R(I)=F0_R(I)+IND(K)*FR(0,K) F1_R(I)=F1_R(I)+IND(K)*FR(1,K) F2_R(I)=F2_R(I)+IND(K)*FR(2,K) F3_R(I)=F3_R(I)+IND(K)*FR(3,K) F4_R(I)=F4_R(I)+IND(K)*FR(4,K) END DO F0_R(I)=F0_R(I)*RADEG F1_R(I)=F1_R(I)*RASEC F2_R(I)=F2_R(I)*RASEC F3_R(I)=F3_R(I)*RASEC F4_R(I)=F4_R(I)*RASEC END DO 6 I_R=I CLOSE (1) C C READING COMPLETE SERIES FOR COORDINATE V (LUNAR ECLIPTIC LONGITUDE) C OPEN (1,FILE='table7.dat',STATUS='OLD') DO WHILE (.TRUE.) READ (1,'(I6,2X,5I3,1X,8I3,1X,I3,F16.7,2F11.6,3F19.12)',END=7) * I,IND(12),IND(11),IND(13),IND(10),IND(9),(IND(J),J=1,8),IND(14) * ,A_V(I),AT_V(I),ATT_V(I),C_V(I),CT_V(I),CTT_V(I) C_V(I)=C_V(I)*RADEG CT_V(I)=CT_V(I)*RADEG CTT_V(I)=CTT_V(I)*RADEG C C CALCULATING THE ARGUMENTS [RAD] C F0_V(I)=0.D0 F1_V(I)=0.D0 F2_V(I)=0.D0 F3_V(I)=0.D0 F4_V(I)=0.D0 DO K=1,14 F0_V(I)=F0_V(I)+IND(K)*FR(0,K) F1_V(I)=F1_V(I)+IND(K)*FR(1,K) F2_V(I)=F2_V(I)+IND(K)*FR(2,K) F3_V(I)=F3_V(I)+IND(K)*FR(3,K) F4_V(I)=F4_V(I)+IND(K)*FR(4,K) END DO F0_V(I)=F0_V(I)*RADEG F1_V(I)=F1_V(I)*RASEC F2_V(I)=F2_V(I)*RASEC F3_V(I)=F3_V(I)*RASEC F4_V(I)=F4_V(I)*RASEC FRM0=FRM(0)*3600.D0 END DO 7 I_V=I CLOSE (1) C C READING COMPLETE SERIES FOR COORDINATE U (LUNAR ECLIPTIC LATITUDE) C OPEN (1,FILE='table8.dat',STATUS='OLD') DO WHILE (.TRUE.) READ (1,'(I6,2X,5I3,1X,8I3,1X,I3,F16.7,2F11.6,3F19.12)',END=8) * I,IND(12),IND(11),IND(13),IND(10),IND(9),(IND(J),J=1,8),IND(14) * ,A_U(I),AT_U(I),ATT_U(I),C_U(I),CT_U(I),CTT_U(I) C_U(I)=C_U(I)*RADEG CT_U(I)=CT_U(I)*RADEG CTT_U(I)=CTT_U(I)*RADEG C C CALCULATING THE ARGUMENTS [RAD] C F0_U(I)=0.D0 F1_U(I)=0.D0 F2_U(I)=0.D0 F3_U(I)=0.D0 F4_U(I)=0.D0 DO K=1,14 F0_U(I)=F0_U(I)+IND(K)*FR(0,K) F1_U(I)=F1_U(I)+IND(K)*FR(1,K) F2_U(I)=F2_U(I)+IND(K)*FR(2,K) F3_U(I)=F3_U(I)+IND(K)*FR(3,K) F4_U(I)=F4_U(I)+IND(K)*FR(4,K) END DO F0_U(I)=F0_U(I)*RADEG F1_U(I)=F1_U(I)*RASEC F2_U(I)=F2_U(I)*RASEC F3_U(I)=F3_U(I)*RASEC F4_U(I)=F4_U(I)*RASEC END DO 8 I_U=I CLOSE (1) END IF IF (LOADDATA.AND.MODEL.EQ.1) THEN ! SIMPLIFIED ANALYTICAL SERIES LEA-406b ARE USED C C READING SIMPLIFIED SERIES FOR COORDINATE R (LUNAR DISTANCE) C OPEN (1,FILE='table9.dat',STATUS='OLD') DO WHILE (.TRUE.) READ (1,'(I6,2X,5I3,1X,8I3,1X,I3,F16.7,2F11.6,3F19.12)',END=9) * I,IND(12),IND(11),IND(13),IND(10),IND(9),(IND(J),J=1,8),IND(14) * ,A_R(I),AT_R(I),ATT_R(I),C_R(I),CT_R(I),CTT_R(I) C_R(I)=C_R(I)*RADEG CT_R(I)=CT_R(I)*RADEG CTT_R(I)=CTT_R(I)*RADEG C C CALCULATING THE ARGUMENTS [RAD] C F0_R(I)=0.D0 F1_R(I)=0.D0 F2_R(I)=0.D0 F3_R(I)=0.D0 F4_R(I)=0.D0 DO K=1,14 F0_R(I)=F0_R(I)+IND(K)*FR(0,K) F1_R(I)=F1_R(I)+IND(K)*FR(1,K) F2_R(I)=F2_R(I)+IND(K)*FR(2,K) F3_R(I)=F3_R(I)+IND(K)*FR(3,K) F4_R(I)=F4_R(I)+IND(K)*FR(4,K) END DO F0_R(I)=F0_R(I)*RADEG F1_R(I)=F1_R(I)*RASEC F2_R(I)=F2_R(I)*RASEC F3_R(I)=F3_R(I)*RASEC F4_R(I)=F4_R(I)*RASEC END DO 9 I_R=I CLOSE (1) C C READING SIMPLIFIED SERIES FOR COORDINATE V (LUNAR ECLIPTIC LONGITUDE) C OPEN (1,FILE='table10.dat',STATUS='OLD') DO WHILE (.TRUE.) READ (1,'(I6,2X,5I3,1X,8I3,1X,I3,F16.7,2F11.6,3F19.12)',END=10) * I,IND(12),IND(11),IND(13),IND(10),IND(9),(IND(J),J=1,8),IND(14) * ,A_V(I),AT_V(I),ATT_V(I),C_V(I),CT_V(I),CTT_V(I) C_V(I)=C_V(I)*RADEG CT_V(I)=CT_V(I)*RADEG CTT_V(I)=CTT_V(I)*RADEG C C CALCULATING THE ARGUMENTS [RAD] C F0_V(I)=0.D0 F1_V(I)=0.D0 F2_V(I)=0.D0 F3_V(I)=0.D0 F4_V(I)=0.D0 DO K=1,14 F0_V(I)=F0_V(I)+IND(K)*FR(0,K) F1_V(I)=F1_V(I)+IND(K)*FR(1,K) F2_V(I)=F2_V(I)+IND(K)*FR(2,K) F3_V(I)=F3_V(I)+IND(K)*FR(3,K) F4_V(I)=F4_V(I)+IND(K)*FR(4,K) END DO F0_V(I)=F0_V(I)*RADEG F1_V(I)=F1_V(I)*RASEC F2_V(I)=F2_V(I)*RASEC F3_V(I)=F3_V(I)*RASEC F4_V(I)=F4_V(I)*RASEC FRM0=FRM(0)*3600.D0 END DO 10 I_V=I CLOSE (1) C C READING SIMPLIFIED SERIES FOR COORDINATE U (LUNAR ECLIPTIC LATITUDE) C OPEN (1,FILE='table11.dat',STATUS='OLD') DO WHILE (.TRUE.) READ (1,'(I6,2X,5I3,1X,8I3,1X,I3,F16.7,2F11.6,3F19.12)',END=11) * I,IND(12),IND(11),IND(13),IND(10),IND(9),(IND(J),J=1,8),IND(14) * ,A_U(I),AT_U(I),ATT_U(I),C_U(I),CT_U(I),CTT_U(I) C_U(I)=C_U(I)*RADEG CT_U(I)=CT_U(I)*RADEG CTT_U(I)=CTT_U(I)*RADEG C C CALCULATING THE ARGUMENTS [RAD] C F0_U(I)=0.D0 F1_U(I)=0.D0 F2_U(I)=0.D0 F3_U(I)=0.D0 F4_U(I)=0.D0 DO K=1,14 F0_U(I)=F0_U(I)+IND(K)*FR(0,K) F1_U(I)=F1_U(I)+IND(K)*FR(1,K) F2_U(I)=F2_U(I)+IND(K)*FR(2,K) F3_U(I)=F3_U(I)+IND(K)*FR(3,K) F4_U(I)=F4_U(I)+IND(K)*FR(4,K) END DO F0_U(I)=F0_U(I)*RADEG F1_U(I)=F1_U(I)*RASEC F2_U(I)=F2_U(I)*RASEC F3_U(I)=F3_U(I)*RASEC F4_U(I)=F4_U(I)*RASEC END DO 11 I_U=I CLOSE (1) END IF C C CALCULATION OF THE COORDINATE R (LUNAR DISTANCE) AT THE INPUT EPOCH [KM] C R=0.D0 DO I=1,I_R ARG=(((F4_R(I)*DTC+F3_R(I))*DTC+F2_R(I))*DTC+F1_R(I))*DTC+F0_R(I) R=R+A_R(I)*DCOS(ARG+C_R(I)) IF (AT_R(I).GT.0.D0) R=R+AT_R(I)*DTM*DCOS(ARG+CT_R(I)) IF (ATT_R(I).GT.0.D0) R=R+ATT_R(I)*DTM*DTM*DCOS(ARG+CTT_R(I)) END DO C C CALCULATION OF THE COORDINATE V (LUNAR ECLIPTIC LONGITUDE) AT THE INPUT EPOCH [RAD] C V=FRM0+(((FRM(4)*DTC+FRM(3))*DTC+FRM(2))*DTC+FRM(1))*DTC DO I=1,I_V ARG=(((F4_V(I)*DTC+F3_V(I))*DTC+F2_V(I))*DTC+F1_V(I))*DTC+F0_V(I) V=V+A_V(I)*DSIN(ARG+C_V(I)) IF (AT_V(I).GT.0.D0) V=V+AT_V(I)*DTM*DSIN(ARG+CT_V(I)) IF (ATT_V(I).GT.0.D0) V=V+ATT_V(I)*DTM*DTM*DSIN(ARG+CTT_V(I)) END DO V=DMOD(V*RASEC,PI) C C CALCULATION OF THE COORDINATE U (LUNAR ECLIPTIC LATITUDE) AT THE INPUT EPOCH [RAD] C U=0.D0 DO I=1,I_U ARG=(((F4_U(I)*DTC+F3_U(I))*DTC+F2_U(I))*DTC+F1_U(I))*DTC+F0_U(I) U=U+A_U(I)*DSIN(ARG+C_U(I)) IF (AT_U(I).GT.0.D0) U=U+AT_U(I)*DTM*DSIN(ARG+CT_U(I)) IF (ATT_U(I).GT.0.D0) U=U+ATT_U(I)*DTM*DTM*DSIN(ARG+CTT_U(I)) END DO U=DMOD(U*RASEC,PI) C C CALCULATION OF GEOCENTRIC POSITION OF THE MOON RELATIVE TO THE ECLIPTIC AND EQUINOX OF EPOCH C P0(1)=R*DCOS(U)*DCOS(V) P0(2)=R*DCOS(U)*DSIN(V) P0(3)=R*DSIN(U) C C CALCULATION OF GEOCENTRIC POSITION OF THE MOON RELATIVE TO THE GEOEQUATOR AND EQUINOX OF EPOCH C EPSD=-EPS(DTM) ! THE OBLIQUITY AT EPOCH (SIMON et al. 1994A&A..282..663) WITH OPPOSITE SIGN CE=DCOS(EPSD) SE=DSIN(EPSD) P1(1)=P0(1) P1(2)=CE*P0(2)+SE*P0(3) P1(3)=-SE*P0(2)+CE*P0(3) C C CALCULATION OF GEOCENTRIC POSITION OF THE MOON RELATIVE TO THE GEOEQUATOR AND EQUINOX OF J2000 C CALL PREC(TDB_0,TDB,PRM) ! CALCULATION OF PRECESSION MATRIX (SIMON et al. 1994A&A..282..663) CALL VM(P1,PRM,P) END BLOCK DATA CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C BLOCK DATA ARRAY FOR *LEA-406* C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMPLICIT REAL*8(A-H,O-Z) COMMON /BLOCK_DATA/ PI, RADEG, RASEC, FR(0:4,14), FRM(0:4) DATA PI/6.283185307179586D0/, * RADEG/1.7453292519943296D-2/, RASEC/4.8481368110953599D-6/ C C SOURCE FOR FREQUENCIES GIVEN BELOW: SIMON ET AL. 1994A&A..282..663 ("1992" VALUES") C THE UNITS ARE DEG, ARCSEC/CENTURY, ARCSEC/CENTURY^2, ARCSEC/CENTURY^3, ARCSEC/CENTURY^4 C DATA (FR(I,1),I=0,4) /252.25090552D0, ! FREQUENCIES FOR L_MERCURY (J2000) * 538101628.688982D0, -0.0192789D0, 0.00000639D0, 0.D0/ DATA (FR(I,2),I=0,4) /181.97980085D0, ! FREQUENCIES FOR L_VENUS (J2000) * 210664136.433548D0, 0.0059381D0, -0.00000627D0, 0.D0/ DATA (FR(I,3),I=0,4) /100.46645683D0, ! FREQUENCIES FOR L_EARTH (J2000) * 129597742.283429D0, -0.0204411D0, -.00000523D0, 0.D0/ DATA (FR(I,4),I=0,4) /355.43299958D0, ! FREQUENCIES FOR L_MARS (J2000) * 68905077.493988D0, 0.0094264D0, -.00001043D0, 0.D0/ DATA (FR(I,5),I=0,4) / 34.35151874D0, ! FREQUENCIES FOR L_JUPITER (J2000) * 10925660.377991D0, -0.3060378D0, 0.00005706D0, .000004667D0/ DATA (FR(I,6),I=0,4) / 50.07744430D0, ! FREQUENCIES FOR L_SATURN (J2000) * 4399609.855732D0, 0.7561614D0, -.00016618D0, -.000011484D0/ DATA (FR(I,7),I=0,4) /314.05500511D0, ! FREQUENCIES FOR L_URANUS (J2000) * 1542481.193933D0, -0.0175083D0, 0.00002156D0, 0.D0/ DATA (FR(I,8),I=0,4) /304.34866548D0, ! FREQUENCIES FOR L_NEPTUNE (J2000) * 786550.320744D0, 0.0021103D0, -0.00000895D0, 0.D0/ DATA (FR(I,9),I=0,4) /125.04455501D0, !FREQUENCIES FOR OMEGA_LUNE (J2000) * -6967919.3631D0, 6.3602D0, .007625D0, -.00003586D0/ DATA (FR(I,10),I=0,4) /297.85019547D0, ! FREQUENCIES FOR DELAUNAY ARG. D * 1602961601.2090D0, -6.3706D0, 0.006593D0,-0.00003169D0/ DATA (FR(I,11),I=0,4) /357.52910918D0, ! FREQUENCIES FOR DELAUNAY ARG. l' * 129596581.0481D0, -0.5532D0, 0.000136D0,-0.00001149D0/ DATA (FR(I,12),I=0,4) /134.96340251D0, ! FREQUENCIES FOR DELAUNAY ARG. l * 1717915923.2178D0, 31.8792D0, 0.051635D0,-0.00024470D0/ DATA (FR(I,13),I=0,4) / 93.27209062D0, ! FREQUENCIES FOR DELAUNAY ARG. F * 1739527262.8478D0,-12.7512D0,-0.001037D0, 0.00000417D0/ DATA (FR(I,14),I=0,4) / 0.D0, ! FREQUENCIES FOR GENERAL PRECISION Pa * 5028.82D0, 1.112022D0, 0.0000773D0, -0.00002353D0/ DATA (FRM(I),I=0,4) /218.31664563D0, ! FREQUENCIES FOR L_MOON (Williams 1991 VALUE FOR Pa IS USED) * 1732564372.30470D0,-5.2790D0, 0.006665D0,-0.00005522D0/ END DOUBLE PRECISION FUNCTION EPS(DTM) C FUNCTION RETURNS THE OBLIQUITY OF ECLIPTIC AT EPOCH C DIFFERING FROM J2000.0 (JD2451545) BY VALUE DTM. C C DATA SOURCE: SIMON ET AL. 1994A&A..282..66 C C CALLING SEQUENCE PARAMETERS: C C DTM = INPUT D.P. TIME DIFFERENCE FROM THE EPOCH J2000.0 [THOUSANDS OF YEARS]. C C EPS = OUTPUT (FUNCTION) D.P. VALUE OF THE OBLIQUITY [RAD] IMPLICIT REAL*8 (A-H,O-Z) PARAMETER *(EPS0=(23.D0+26.D0/60.D0+21.412D0/3600.D0)*1.7453292519943296D-2) COMMON /BLOCK_DATA/ PI, RADEG, RASEC, FR(0:4,14), FRM(0:4) EPS=EPS0+((((-.0025D0*DTM-.0051D0)*DTM+1.9989D0)*DTM-.0152D0)*DTM * -468.0927D0)*DTM*RASEC RETURN END SUBROUTINE PREC(TDB_0,TDB_1,PM) C SUBROUTINE CALCULATES THE PRECESSION MATRIX DEFINING POSITION OF THE C MEAN EQUINOX AND GEOEQUATOR OF EPOCH TDB_1 RELATIVE TO THE MEAN EQUINOX C AND EQUATOR OF INITIAL EPOCH TDB_0. C C DATA SOURCE: SIMON ET AL. 1994A&A..282..663 C C CALLING SEQUENCE PARAMETERS: C C TDB_0 = INPUT 2-WORD D.P. TDB DATE OF INITIAL (FIXED) EPOCH. C C TDB_1 = INPUT 2-WORD D.P. TDB DATE OF ANOTHER EPOCH. C C PM = OUTPUT (3x3) D.P. PRECESSION MATRIX SO THAT: C C R_1 = PM x R_0 C WHERE C R_0 IS A 3-WORD COORDINATE VECTOR OF AN ARBITRARY POINT IN THE REFERENCE C FRAME DEFINED BY THE MEAN EQUINOX AND GEOEQUATOR OF EPOCH TDB_0, C AND C R_1 IS THE 3-WORD COORDINATE VECTOR OF THE SAME POINT IN THE REFERENCE C FRAME DEFINED BY THE MEAN EQUINOX AND GEOEQUATOR OF EPOCH TDB_1. IMPLICIT REAL*8 (A-H,O-Z) DIMENSION TDB_0(2),TDB_1(2),PM(3,3) COMMON /BLOCK_DATA/ PI, RADEG, RASEC, FR(0:4,14), FRM(0:4) C C CALCULATION OF TIME DIFFERENCES C T=((TDB_0(1)-2451545.D0)+TDB_0(2))/365250.D0 T2=T*T T3=T2*T T4=T3*T T5=T4*T D=((TDB_1(1)-TDB_0(1))+(TDB_1(2)-TDB_0(2)))/365250.D0 D2=D*D D3=D2*D D4=D3*D D5=D4*D D6=D5*D C C CALCULATION OF PRECESSION QUANTITIES C DZA=(23060.9097D0+139.7495D0*T-.0038D0*T2 * -0.5918D0*T3-0.0037D0*T4+0.0007D0*T5)*D ZA=(DZA+(109.5270D0+.2446D0*T-1.3913D0*T2-.0134D0*T3 * +.0026D0*T4)*D2 * +(18.2667D0-1.1400D0*T-0.0173D0*T2+0.0044D0*T3)*D3 * +(-0.2821D0-0.0093D0*T+0.0032D0*T2)*D4 * +(-0.0301D0+0.0006D0*T)*D5-0.0001D0*D6)*RASEC DZA=(DZA+(30.2226D0-.2523D0*T-.3840D0*T2-.0014D0*T3 * +.0007D0*T4)*D2 * +(18.0183D0-.1326D0*T+0.0006D0*T2+0.0005D0*T3)*D3 * +(-0.0583D0-0.0001D0*T+0.0007D0*T2)*D4 * -0.0285D0*D5-0.0002D0*D6)*RASEC TA=((20042.0207D0-85.3131D0*T-.2111D0*T2+0.3642D0*T3 * +.0008D0*T4-.0005D0*T5)*D * +(-42.6566D0-.2111D0*T+.5463D0*T2+.0017D0*T3 * -.0012D0*T4)*D2 * +(-41.8238D0+.0359D0*T+0.0027D0*T2-0.0001D0*T3)*D3 * +(-0.0731D0+0.0019D0*T+0.0009D0*T2)*D4 * +(-0.0127D0+0.0011D0*T)*D5+0.0004D0*D6)*RASEC C C CALCULATION OF PRECESSION MATRIX C CDZ=DCOS(DZA) SDZ=DSIN(DZA) CZ=DCOS(ZA) SZ=DSIN(ZA) CT=DCOS(TA) ST=DSIN(TA) PM(1,1)=-SDZ*SZ+CDZ*CZ*CT PM(1,2)=-CDZ*SZ-SDZ*CZ*CT PM(1,3)=-CZ*ST PM(2,1)=SDZ*CZ+CDZ*SZ*CT PM(2,2)=CDZ*CZ-SDZ*SZ*CT PM(2,3)=-SZ*ST PM(3,1)=CDZ*ST PM(3,2)=-SDZ*ST PM(3,3)=CT RETURN END SUBROUTINE VM(A,B,C) C SUBROUTINE MULTIPLIES VECTOR AND MATRIX C C CALLING SEQUENCE PARAMETERS: C C A = INPUT D.P. VECTOR (3) C B = INPUT D.P. MATRIX (3*3) C C C = OUTPUT D.P. VECTOR (3) DOUBLE PRECISION A(3),B(3,3),C(3) C(1)=A(1)*B(1,1)+A(2)*B(2,1)+A(3)*B(3,1) C(2)=A(1)*B(1,2)+A(2)*B(2,2)+A(3)*B(3,2) C(3)=A(1)*B(1,3)+A(2)*B(2,3)+A(3)*B(3,3) RETURN END