!
! --- PROGRAM VSOP2013_compute -----------------------------------------
!
!     COMPUTATION OF PLANETARY COORDINATES VSOP2013 
!     WITH THE CHEBYSHEV EPHEMERIDES VSOP2013.
!
!     THE USER HAS TO CHOICE THE VALUES OF THE PARAMETERS
!     NAMEF : NAME OF THE SEQUENTIAL FILE (character*40)
!     IRECL : INDEX OF RECORD LENGH IN THE DIRECT ACCESS FILE (integer).
!     TZERO : JULIAN DATE OF THE TIME ORIGIN (real*8)
!
!     *** CAUTION ***
!     The length of the records of the direct acces file is expressed 
!     in unit which depends on the Fortran compiler: 
!     7840 in bytes or 
!     3920 in words of 4 bytes or
!     1960 in words of 2 bytes.
!
! --- Declarations -----------------------------------------------------
!
      implicit none
!
      integer                      :: n,ipla,ndat,maxpla,nul,ierr,irecl
      real*8                       :: TZERO,step,jd,time
      character*12                 :: message
      character*40                 :: namef
!
      real*8,        dimension(6)  :: r
      character*10,  dimension(9)  :: planet
!
      real*8,        parameter     :: jd2000=2451545.d0
      character*40,  parameter     :: namef0='VSOP2013.p2000.bin'
!
      data planet/
     &     'MERCURY','VENUS','EARTH-MOON','MARS','JUPITER',
     &     'SATURN','URANUS','NEPTUNE','PLUTO'/
!
! --- Parameters -------------------------------------------------------
!
      NAMEF='VSOP2013.m4000.bin'  ! [-4500 -3500]
      NAMEF='VSOP2013.m2000.bin'  ! [-3500 -1500]
      NAMEF='VSOP2013.m1000.bin'  ! [-1500  0]
      NAMEF='VSOP2013.p1000.bin'  ! [ 0  1500]
      NAMEF='VSOP2013.p2000.bin'  ! [+1500 +3000]
      NAMEF='VSOP2013.p4000.bin'  ! [+3000 +4500]
!
      IRECL=4  !!! *** record length in words of 4 bytes.
      IRECL=2  !!! *** record length in words of 2 bytes.
      IRECL=1  !!! *** record length in bytes.
!
      TZERO=77432.5d0     ! -4500/01/01 0h
      TZERO=625307.5d0    ! -3000/01/01 0h
      TZERO=1173182.5d0   ! -1500/01/01 0h
      TZERO=1721057.5d0   !  0000/01/01 0h
      TZERO=2268932.5d0   !  1500/01/01 0h
      TZERO=2816818.5d0   !  3000/02/01 0h
!
! --- Parameters -------------------------------------------------------
!
      ndat=5
      step=136798.d0
!
      open (10,file='VSOP2013_compute.txt')
      nul=20
!
      maxpla=9
      if (NAMEF.ne.namef0) maxpla=8
!
! --- Computation ------------------------------------------------------
!
      do ipla=1,maxpla
         do n=1,ndat
            jd=TZERO+(n-1)*step
            time=jd-jd2000
            call EPHVSOP2013 (time,ipla,NAMEF,nul,IRECL,r,ierr)
            if (ierr.ne.0) then
               if (ierr.eq.10) message='Error File'
               if (ierr.eq.11) message='Error Time'
               if (ierr.eq.12) message='Error Planet'
               if (ierr.eq.13) message='Error irecl index'
               write (10,1001) planet(ipla),jd,message
               stop
            endif
            write (10,1002) planet(ipla),jd,r
         enddo
      enddo
!
      stop
!
! --- Formats ----------------------------------------------------------
!
1001  format (/2x,a10,2x,'JD',f9.1,2x,a12)
1002  format (2x,a10,2x,'JD',f9.1,
     &        2x,'X :',f16.12,' ua  ',
     &        2x,'Y :',f16.12,' ua  ',
     &        2x,'Z :',f16.12,' ua  '/
     &        27x,'X'':',f16.12,' ua/d',
     &        2x,'Y'':',f16.12,' ua/d',
     &        2x,'Z'':',f16.12,' ua/d')
!
      end
!
!
!
      subroutine EPHVSOP2013 (djb,ipla,namef,nul,irecl,r,ierr)
!-----------------------------------------------------------------------
!
!     Reference : GF-JLS-1212
!
! --- Object -----------------------------------------------------------
!
!     Heliocentric rectangular coordinates of planets.
!     Computation using a file of Chebyshev Ephemeris VSOP2013.
!     Frame : Dynamical Ecliptic and Equinox J2000.
!
! --- Input -----------------------------------------------------------
!
!     djb     : Time TDB in days from the origin J2000 (real*8).
!     ipla    : Index of the planet (integer).
!               ipla=1 Mercury     ipla=5 Jupiter
!               ipla=2 Venus       ipla=6 Saturn
!               ipla=3 E-M Baryc.  ipla=7 Uranus
!               ipla=4 Mars        ipla=8 Neptune
!                                  ipla=9 Pluto
!
!     namef   : Name of the File of Chebyshev Ephemeris (character*40).
!     nul     : Unit of the File of Chebyshev Ephemeris (integer).
!     irecl   : record lengh index (integer).
!               irecl=1 record size in bytes.
!               irecl=2 record size in words of 2 bytes.
!               irecl=4 record size in words of 4 bytes.
!
! --- Output -----------------------------------------------------------
!
!     r(6)    : Table of the results (real*8).
!               r(1:3) : position  XYZ (au).
!               r(4:6) : velocity  XYZ (au/j).
!
!     ierr    : Error Index (integer).
!               ierr=00 no error.
!               ierr=10 error in the file VSOP2013.
!               ierr=11 error in the parameter djb.
!               ierr=12 error in the parameter ipla.
!               ierr=13 error in the parameter irecl.
!
! --- Declarations -----------------------------------------------------
!
      implicit none
!
      logical                            :: fexist
!
      real*8                             :: djb
      integer                            :: ipla,nul,irecl,ierr
      character*40                       :: namef
      real*8,       dimension(6)         :: r
!
      integer                            :: i,ik,j,jp,jt,n
      integer                            :: nerr,iad,ncf,nsi
      integer                            :: iloc,irec,nul0,ir0
      real*8                             :: dj0,dj1,dj2,x
      character*40                       :: namef0
!
      real*8,       dimension(20)        :: tn
      real*8,       dimension(9,3)       :: resul
!
      real*8,       parameter            :: t2000=2451545.d0
!
! --- Parameters of the File of Chebyshev Ephemeris VSOP2013 -----------
!
      integer                            :: idf,nintv,ncoef,lrec
      integer,      dimension(3,9)       :: loc
      real*8                             :: t1,t2,delta,delta2
      real*8,       dimension(978)       :: coef
      integer,      parameter            :: krecl=7840
!
! --- Initializations --------------------------------------------------
!
      data namef0 /'***'/
      data ir0  /0/
      data nul0 /0/
!
      save
!
! --- Reading the Header -----------------------------------------------
!
      if (namef.ne.namef0.or.nul.ne.nul0) then
         ierr=13
         if (irecl.ne.1.and.irecl.ne.2.and.irecl.ne.4) return
         lrec=krecl/irecl
         ierr=10
         inquire (file=namef,exist=fexist)
         if (.not.fexist) return
         open (nul,file=namef,access='direct',recl=lrec,iostat=nerr)
         if (nerr.ne.0) return
         read (nul,rec=1,iostat=nerr)  idf,t1,t2,delta,nintv,ncoef,loc
         if (nerr.ne.0) return
         namef0=namef
         nul0=nul
      endif
!
! --- Error in the parameters ------------------------------------------
!
      ierr=11
      if (djb.lt.(t1-t2000).or.djb.ge.(t2-t2000)) return
      ierr=12
      if (ipla.lt.1.or.ipla.gt.9) return
      ierr=10
!
! --- Reading the Chebychev coefficients -------------------------------
!
      irec=(djb-(t1-t2000))/delta+2
      if (djb.eq.(t2-t2000)) irec=nintv+1
      if (irec.ne.ir0) then
         read (nul,rec=irec,iostat=nerr) dj1,dj2,coef   
         if (nerr.ne.0) return
         ir0=irec
      endif
!
! --- Computation : Poistions & Velocity -------------------------------
!
      iad=loc(1,ipla)
      ncf=loc(2,ipla)
      nsi=loc(3,ipla)
      delta2=delta/nsi
!
      ik=(djb-(dj1-t2000))/delta2
      if (ik.eq.nsi) ik=ik-1
      iloc=iad+6*ncf*ik
      dj0=dj1+ik*delta2-t2000
!
      x=2.d0*(djb-dj0)/delta2-1.d0
      tn(1)=1.d0
      tn(2)=x
      do i=3,ncf
         tn(i)=2.d0*x*tn(i-1)-tn(i-2)
      enddo
!
      do i=1,6
         r(i)=0.d0
         do j=1,ncf
            jp=ncf-j+1
            jt=iloc+ncf*(i-1)+jp-1
            r(i)=r(i)+tn(jp)*coef(jt)
         enddo
      enddo
!
! --- Exit -------------------------------------------------------------
!
      ierr=0
      return
!
      end subroutine EPHVSOP2013
