      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
