PROGRAM LEA406a

C EXAMPLE OF CALLING THE 'LEA406.FOR' SUBROUTINE: 
C CALCULATION OF LUNAR RECTANGULAR COORDINATES ON THE BASE OF COMPLETE LEA-406a SERIES

      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION TDB(2), P(3) 
      DIMENSION TLIM(3)

      DATA MODEL/0/ ! COMPLETE LEA-406a SERIES ARE USED      
      DATA DATE_MIN/2268932.5D0/, DATE_MAX/2634166.5D0/, STEP/36525.D0/ ! [JD]  [1500 - 2500]
C
C INITIALIZATION OF TDB ARRAY (2-WORD D.P. JULIAN DATE (TDB) AT WHICH LUNAR POSITION IS WANTED) 
C
      TDB(1)=2451545.D0
C Get the 4 values min, max, step, model from standard input 
      read(5,'(F16.8)', end=7) T
      TLIM(1) = T
      if (T.gt.DATE_MAX .or. T.lt.DATE_MIN) goto 9
      read(5,'(F16.8)', end=7) T
      TLIM(2) = T
      if (T.gt.DATE_MAX .or. T.lt.DATE_MIN) goto 9
      read(5,'(F16.8)', end=7) T
      STEP = T
      read(5, '(I1)',end=8) M
      MODEL=M
    8 if (MODEL.lt.0 .or. MODEL.gt.1) MODEL=0
      if (MODEL.eq.0) then
          write(*,'(A)') ' ====Complete LEA-406a series used'
      else
          write(*,'(A)') ' ====Simplified LEA-406 model used'
      endif
      
C
C CALCULATION OF GEOCENTRIC POSITION OF 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 
      NSTEP = 0
      TDB(2)=TLIM(1)-TDB(1)
    1 DO WHILE (TDB(1)+TDB(2).LE.TLIM(2) .and. NSTEP.lt.1000) 

       if (NSTEP.eq.0) WRITE(*,11) 'x[km]','y[km]','z[km]'
       CALL LEA406(TDB,MODEL,P)
       NSTEP = NSTEP+1
C
C PRINTING THE RESULTS
C
       WRITE (*,'(F17.8,3F17.5,F12.5)') TDB(1)+TDB(2), P

       TDB(2)=TDB(2)+STEP

      END DO 

      STOP

   7  WRITE (*,79) 'JDmin', 'JDmax', 'JSstep'
      STOP 1
   9  WRITE (*,99) DATE_MIN, DATE_MAX,T
  11  FORMAT(' JulianDate[TDB]',3A17)
  79  FORMAT(' +++Give on 3 consecutive lines the values of', 3(/A7))
  99  FORMAT(' ***Date not within JD',F9.1,' and ',F9.1,':',F9.1)
      STOP 99
      END
C
C                               RESULTS OF CALCULATION:
C
C   GEOCENTRIC RECTANGULAR COORDINATES OF THE MOON IN THE MEAN GEOEQUATOR AND EQUINOX OF J2000.0 
C   CALCULATED ON THE BASE OF COMPLETE ANALYTICAL SERIES LEA-406a (KUDRYAVTSEV 2007)
C
C   JULIAN DATE, TDB       X, KM            Y, KM            Z, KM         |   DIFFERENCE FROM DE-406, KM 
C   2268932.50000000     135292.91277    -298773.33126    -146998.83620    |   0.00310 
C   2305457.50000000    -239090.32567    -276225.69944    -158278.99337    |   0.00230 
C   2341982.50000000    -370191.65009     -70677.12970     -17271.73161    |   0.00240 
C   2378507.50000000    -222885.44533     276871.03659     157142.92190    |   0.00248 
C   2415032.50000000      55225.43107     355859.66895     154093.43547    |   0.00166 
C   2451557.50000000     364406.63531     109457.02311      10993.11763    |   0.00108 
C   2488082.50000000     354506.70853    -161535.36163     -79477.66354    |   0.00132 
C   2524607.50000000      15744.98875    -339021.57922    -117595.75588    |   0.00260 
C   2561132.50000000    -292971.78212    -264171.90202     -90864.24938    |   0.00185 
C   2597657.50000000    -362541.62407      54645.60490      -4857.32320    |   0.00085 
C   2634166.50000000     -34988.77296    -337212.38231    -111990.81334    |   0.00137