*     ---------------
*     PROGRAM TESTEPH
*     ---------------
*
*     BDL-GF9612
*
*     This program allows to compute Solar, Lunar and Planetary
*     coordinates with the Chebychev Ephemerides SLP96 read in an
*     direct access binary file SLP96 and to compare them with the
*     test values.
*
*     Uses the subroutine SLP96.
*
* --- Declarations -----------------------------------------------------
*
      implicit double precision (a-h,o-z)
      logical fexist
      character*40 fname1,fname2
*
      dimension r1(3),r2(6),d(3)
*
      parameter (prec=1.d-12,lrec=5576)
*
      keph=0
      kerr=0
*
* --- Open the direct access binary file SLP96 -------------------------
*
100   continue
*
      write (*,'(//2x,''SLP96 : TESTEPH (BDL-GF9612)''/)')
*
      write (*,'(/2x,''Name of the SLP96 binary file  : '',$)')
      read (*,'(a40)') fname1
*
      inquire (file=fname1,exist=fexist)
      if (.not.fexist) then
         write (*,'(/2x,''Error binary file'')')
         goto 100
      endif
*
      open (10,file=fname1,access='direct',recl=lrec,status='old',
     .      iostat=nerr)
*
      if (nerr.ne.0) then
         write (*,'(/2x,''Error binary file'')')
         goto 100
      endif
*
      read (10,rec=1,iostat=nerr) idf,tla,tlb
*
      if (nerr.ne.0) then
         write (*,'(/2x,''Error binary file'')')
         goto 100
      endif
*
      close (10)
*
* --- Open the test values file ----------------------------------------
*
200   continue
*
      write (*,'(/2x,''Name of the test values file   : '',$)')
      read (*,'(a40)') fname2
*
      inquire (file=fname2,exist=fexist)
      if (.not.fexist) then
         write (*,'(/2x,''Error test values file'')')
         goto 200
      endif
*
      open (20,file=fname2,status='old',iostat=nerr)
*
      if (nerr.ne.0) then
         write (*,'(/2x,''Error test values file'')')
         goto 200
      endif
*
      do i=1,13
         read (20,1000)
      enddo
*
      write (*,'(/2x,''Executing ... '')')
*
*
* --- Compute and compare the coordinates ------------------------------
*
300   continue
*
      read (20,1000,end=400) idate,tjd,ibody,icent,ipv,iframe,icoord,r1
*
      if (tjd.lt.tla) goto 300
      if (tjd.gt.tlb) goto 400
*
      keph=keph+1
      call SLP96 (tjd,ibody,icent,ipv,iframe,icoord,fname1,10,r2,ierr)
*
      if (ierr.ne.0) then
         write (*,'(/2x,''Error'',i3/)') ierr
         stop
      endif
*
      i1=(ipv-1)*3
      i2=ipv*3
*
      do j=1,3
         d(j)=r1(j)-r2(j+i1)
         if (abs(d(j)).gt.prec) then
            kerr=kerr+1
            write (*,'(/2x,''Warning : Difference > '',d8.1)') prec
            write (*,'(/2x,''Computed values : '')')
            write (*,1000) idate,tjd,ibody,icent,ipv,iframe,icoord,
     .                     (r2(i),i=i1+1,i2)
            write (*,'(/2x,''Test     values : '')')
            write (*,1000) idate,tjd,ibody,icent,ipv,iframe,icoord,r1
            write (*,'(/2x,''Press Enter '',$)')
            read (*,*)
            write (*,'(/2x,''Executing ... '')')
            goto 300
         endif
      enddo
*
      goto 300
*
* --- End of comparison ------------------------------------------------
*
400   continue
*
      write (*,'(//2x,''SLP96 : TESTEPH RESULT   prec'',d8.1)') prec
      write (*,'(/2x,''Number of computed coordinates : '',i4.4)') keph
      write (*,'(/2x,''Number of differences > prec   : '',i4.4)') kerr
      write (*,'(/2x,''Program ended''/)')
      stop
*
1000  format (i10,f10.1,2i3,3i2,3f16.12)
*
      end
*
*
*
*     ---------------------------------------------------------------
      subroutine SLP96 (tjd,ibody,icent,ipv,iframe,icoord,fname,lu,
     .                  result,ierr)
*     ---------------------------------------------------------------
*
*     BDL-GF9612
*
*     Solar, Lunar and Planetary coordinates J2000.
*
*     The coordinates are computed with the Chebychev Ephemerides SLP96
*     read in a direct access binary file SLP96.
*
* --- Input ------------------------------------------------------------
*
*     tjd       Julian Date TDB (double real).
*
*     ibody     Body Index (integer).
*
*               ibody=01 : Mercury       ibody=07 : Uranus
*               ibody=02 : Venus         ibody=08 : Neptune
*               ibody=03 : E-M Baryc.    ibody=09 : void
*               ibody=04 : Mars          ibody=10 : Moon
*               ibody=05 : Jupiter       ibody=11 : Sun
*               ibody=06 : Saturn        ibody=12 : Earth
*
*     icent     Frame Center Index (integer).
*
*               icent=00 : Barycenter of Solar System
*               icent=01 : Mercury       icent=07 : Uranus
*               icent=02 : Venus         icent=08 : Neptune
*               icent=03 : E-M Baryc.    icent=09 : void
*               icent=04 : Mars          icent=10 : Moon
*               icent=05 : Jupiter       icent=11 : Sun
*               icent=06 : Saturn        icent=12 : Earth
*
*     ipv       Position-Velocity Index (integer).
*
*               ipv=1 : Position.
*               ipv=2 : Position and Velocity.
*
*     iframe    Frame Index (integer).
*
*               iframe=1 : Equinox and Equator  J2000 (FK5)
*               iframe=2 : Equinox and Ecliptic J2000 (dynamical)
*
*     icoord    Coordinates Index (integer).
*
*               icoord=1 : Rectangular Coordinates.
*               icoord=2 : Spherical Coordinates.
*
*     fname     Name of the file SLP96 (character*40).
*
*     lu        Logical Unit of file SLP96 (integer).
*
* --- Output -----------------------------------------------------------
*
*     result    Results table (double real).
*
*               Rectangular Coordinates (icoord=1) :
*
*               Position :
*               result(1) : Equatorial or Ecliptic Component X (au).
*               result(2) : Equatorial or Ecliptic Component Y (au).
*               result(3) : Equatorial or Ecliptic Component Z (au).
*
*               Velocity :
*               result(4) : Equatorial or Ecliptic Component X (au/d).
*               result(5) : Equatorial or Ecliptic Component Y (au/d).
*               result(6) : Equatorial or Ecliptic Component Z (au/d).
*
*               Spherical Coordinates (icoord=2) :
*
*               Position :
*               result(1) : Right Ascension or Longitude (rad).
*               result(2) : Declination or Latitude (rad).
*               result(3) : Geometric Distance (au).
*
*               Velocity :
*               result(4) : Right Ascension or Longitude (rad/d).
*               result(5) : Declination or Latitude (rad/d).
*               result(6) : Geometric Distance (au/d).
*
*     ierr      Error Index (integer).
*
*               ierr=0  : No error.
*               ierr=10 : File error.
*               ierr=11 : Date error (tjd).
*               ierr=12 : Body error (ibody).
*               ierr=13 : Frame Center error (icent).
*               ierr=14 : Position-Velocity error (ipv).
*               ierr=15 : Frame error (iframe).
*               ierr=16 : Coordinates error (icoord).
*
* --- Declarations -----------------------------------------------------
*
      implicit double precision (a-h,o-z)
*
      character*40 fname
*
      dimension result(6)
      dimension r(6,3),v(6),w(6)
      dimension ic(3),ic0(3),tjd0(3),ipv0(3)
      dimension coef(696),loc(11,3)
      dimension pc(16),vc(16)
      dimension eq2000(3,3)
*
* --- Initializations --------------------------------------------------
*
      data icall/1/,id0/0/
      data ic0/3*0/,tjd0/3*0.d0/,ipv0/3*0/
      data ncp0/0/,ncv0/0/
      data pc/1.d0,2.d0,14*0.d0/
      data vc/0.d0,1.d0,14*0.d0/
      data pc2/2.d0/,vc3/100.d0/
*
* --- Ratio Earth-Moon mass --------------------------------------------
*
      data em/81.3007d0/
*
* --- Rotation matrix (Equator FK5) ------------------------------------
*
      data eq2000/
     .+1.000000000000d0,-0.000000479966d0, 0.000000000000d0,
     .+0.000000440360d0,+0.917482137087d0,+0.397776982902d0,
     .-0.000000190919d0,-0.397776982902d0,+0.917482137087d0/
*
* --- First call -------------------------------------------------------
*
      if (icall.eq.1) then
*
* ------ Read Chebychev Ephemeris parameters ---------------------------
*
         ierr=10
         nbytes=232
         open (lu,file=fname,access='direct',status='old',
     .         recl=nbytes,iostat=nerr)
         if (nerr.ne.0) return
         read (lu,rec=1,iostat=nerr)
     .         idf,tla,tlb,delta,ncoef,nblock,nbytes,loc
         if (nerr.ne.0) return
         close (lu)
         if (idf.ne.96) return
*
* ------ Constants -----------------------------------------------------
*
         pi=3.141592653589793d0
         dpi=pi*2.d0
         f1=-1.d0/(1.d0+em)
         f2=em/(1.d0+em)
*
* ------ Open the direct access binary file SLP96 ----------------------
*
         open (lu,file=fname,access='direct',status='old',
     .         recl=nbytes,iostat=nerr)
         if (nerr.ne.0) return
*
         icall=0
*
      endif
*
* --- Check the parameters ---------------------------------------------
*
      ierr=0
      if (icoord.lt.1.or.icoord.gt.2) ierr=16
      if (iframe.lt.1.or.iframe.gt.2) ierr=15
      if (ipv.lt.1.or.ipv.gt.2) ierr=14
      if (icent.lt.0.or.icent.eq.9.or.icent.gt.12) ierr=13
      if (ibody.lt.1.or.ibody.eq.9.or.ibody.gt.12) ierr=12
      if (tjd.lt.tla.or.tjd.gt.tlb) ierr=11
      if (ierr.ne.0) return
      ncd=3*ipv
      ncase=0
      if (ibody.eq.12.or.icent.eq.12) f=f1
      if (ibody.eq.10.or.icent.eq.10) f=f2
*
* --- Read the Chebychev's coefficients --------------------------------
*
      id=dint((tjd-tla)/delta)+2
      if (id.gt.nblock) id=nblock
      if (id.ne.id0) then
         read (lu,rec=id,iostat=nerr) t0,(coef(i),i=1,ncoef)
         if (nerr.ne.0) then
            ierr=10
            return
         endif
         id0=id
      endif
*
* --- Case : Earth and Moon --------------------------------------------
*
      if ((icent.eq.12.and.ibody.eq.10).or.
     .    (icent.eq.10.and.ibody.eq.12)) then
         ncase=1
         nsub=1
         ic(1)=10
      endif
*
* --- Case : Barycenter of Solar System --------------------------------
*
      if (ncase.eq.0.and.icent.eq.0) then
         ncase=2
         if (ibody.eq.11) then
            nsub=1
            ic(1)=11
         else
            if (ibody.eq.12.or.ibody.eq.10) then
               nsub=3
               ic(1)=11
               ic(2)=3
               ic(3)=10
            else
               nsub=2
               ic(1)=11
               ic(2)=ibody
            endif
         endif
      endif
*
* --- Case : Sun -------------------------------------------------------
*
      if (ncase.eq.0.and.(icent.eq.11.or.ibody.eq.11)) then
         ncase=3
         if (ibody.eq.12.or.icent.eq.12.or.
     .       ibody.eq.10.or.icent.eq.10) then
            nsub=2
            ic(1)=3
            ic(2)=10
         else
            nsub=1
            if (icent.eq.11) ic(1)=ibody
            if (ibody.eq.11) ic(1)=icent
         endif
      endif
*
* --- Case : Earth or Moon ---------------------------------------------
*
      if (ncase.eq.0.and.
     .   (ibody.eq.12.or.icent.eq.12.or.
     .    ibody.eq.10.or.icent.eq.10)) then
         ncase=4
         if (ibody.eq.3.or.icent.eq.3) then
            nsub=1
            ic(1)=10
         else
            nsub=3
            if (icent.eq.12.or.icent.eq.10) ic(1)=ibody
            if (ibody.eq.12.or.ibody.eq.10) ic(1)=icent
            ic(2)=3
            ic(3)=10
         endif
      endif
*
* --- Other cases ------------------------------------------------------
*
      if (ncase.eq.0) then
         ncase=5
         nsub=2
         ic(1)=ibody
         ic(2)=icent
      endif
*
* --- Substitution of time ---------------------------------------------
*
      do n=1,nsub
*
         if (ic(n).ne.ic0(n).or.
     .       tjd.ne.tjd0(n).or.
     .       ipv.gt.ipv0(n)) then

            ic0(n)=ic(n)
            tjd0(n)=tjd
            ipv0(n)=ipv
*
            k=ic(n)
            iad=loc(k,1)
            dlt=loc(k,2)
            ncf=loc(k,3)
*
            ni=dint((tjd-t0)/dlt)+1
            if (tjd.eq.t0+delta) ni=ni-1
            l=iad+3*ncf*(ni-1)
            tz=t0+(ni-1)*dlt
*
* --------- Compute the position ---------------------------------------
*
            cp=2.d0*(tjd-tz)/dlt-1.d0
            pc(2)=cp
*
            if (pc(2).ne.pc2.or.ncf.gt.ncp0) then
               ncp0=ncf
               pc2=pc(2)
               dt=cp+cp
               do i=3,ncf
                  pc(i)=dt*pc(i-1)-pc(i-2)
               enddo
            endif
*
            do i=1,3
               r(i,n)=0.d0
               do j=1,ncf
                  jp=ncf-j+1
                  jt=l+ncf*(i-1)+jp-1
                  r(i,n)=r(i,n)+pc(jp)*coef(jt)
               enddo
            enddo
*
* --------- Compute the velovity ---------------------------------------
*
            if (ipv.eq.2) then
*
               cv=2.d0/dlt
               vc(3)=dt+dt
*
               if (vc(3).ne.vc3.or.ncf.gt.ncv0) then
                  ncv0=ncf
                  vc3=vc(3)
                  do i=4,ncf
                     vc(i)=dt*vc(i-1)+2.d0*pc(i-1)-vc(i-2)
                  enddo
               endif
*
               do i=1,3
                  r(i+3,n)=0.d0
                  do j=2,ncf
                     jp=ncf-j+2
                     jt=l+ncf*(i-1)+jp-1
                     r(i+3,n)=r(i+3,n)+vc(jp)*coef(jt)
                  enddo
                  r(i+3,n)=r(i+3,n)*cv
               enddo
*
            endif
*
         endif
*
      enddo
*
* --- Select the case --------------------------------------------------
*
      select case (ncase)
*
* ------ Case : Earth and Moon -----------------------------------------
*
         case(1)
            is=+1
            if (icent.eq.10) is=-1
            do i=1,ncd
               v(i)=r(i,1)*is
            enddo
*
* ------ Case : Barycenter of Solar System -----------------------------
*
         case(2)
            do i=1,ncd
               v(i)=r(i,1)
               if (ibody.ne.11) then
                  v(i)=v(i)+r(i,2)
                  if (ibody.eq.12.or.ibody.eq.10) then
                     v(i)=v(i)+f*r(i,3)
                  endif
               endif
            enddo
*
* ------ Case : Sun ----------------------------------------------------
*
         case(3)
            is=+1
            if (ibody.eq.11) is=-1
            do i=1,ncd
               v(i)=r(i,1)*is
               if (ibody.eq.12.or.icent.eq.12.or.
     .             ibody.eq.10.or.icent.eq.10) then
                  v(i)=v(i)+f*r(i,2)*is
               endif
            enddo
*
* ------ Case : Earth or Moon ------------------------------------------
*
         case(4)
            is=+1
            if (ibody.eq.12.or.ibody.eq.10) is=-1
            if (ibody.eq.3.or.icent.eq.3) then
               do i=1,ncd
                  v(i)=-f*r(i,1)*is
               enddo
            else
               do i=1,ncd
                  v(i)=(r(i,1)-r(i,2)-f*r(i,3))*is
               enddo
            endif
*
* ------ Other cases ---------------------------------------------------
*
         case (5)
            do i=1,ncd
               v(i)=r(i,1)-r(i,2)
            enddo
*
      endselect
*
* --- Equatorial coordinates -------------------------------------------
*
      if (iframe.eq.1) then
         do l=1,ipv
            k=3*(l-1)
            do i=1,3
               w(i+k)=0.d0
               do j=1,3
                  w(i+k)=w(i+k)+eq2000(i,j)*v(j+k)
               enddo
            enddo
         enddo
      endif
*
* --- Ecliptic coordinates ---------------------------------------------
*
      if (iframe.eq.2) then
         do i=1,ncd
            w(i)=v(i)
         enddo
      endif
*
* --- Rectangular coordinates ------------------------------------------
*
      if (icoord.eq.1) then
         do i=1,ncd
            result(i)=w(i)
         enddo
         return
      endif
*
* --- Spherical coordinates --------------------------------------------
*
      if (icoord.eq.2) then
         p2=w(1)*w(1)+w(2)*w(2)
         p=sqrt(p2)
         d2=p2+w(3)*w(3)
         d=sqrt(d2)
         x=atan2(w(2),w(1))
         if (x.lt.0.0) x=x+dpi
         result(1)=x
         result(2)=atan(w(3)/p)
         result(3)=d
         if (ipv.eq.2) then
            result(4)=(w(1)*w(5)-w(2)*w(4))/p2
            result(5)=(w(6)*p2-w(3)*(w(1)*w(4)+w(2)*w(5)))/p/d2
            result(6)=(w(1)*w(4)+w(2)*w(5)+w(3)*w(6))/d
         endif
      endif
      return
*
      end

