*
*   This program substitutes the time in the variable p of the 
*   precession-nutation series SMART97 (dynamical system).
*
*   Computation for 21 dates from the Julian date 2415020.0 with
*   a 3652.5 day step.
*
*   The series p is read in the file "figdyn97.p"
*
* --- Declarations ---------------------------------------------------
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*50 fich
*
      parameter(t0=2415020.d0,ndat=20,step=3652.5d0)
      parameter(fich='figdyn97.p')
*
      dimension abc(3,5000),itpw(5000),ifi(12),jfi(12,5000),sc(2,5000)
*
      OPEN(10,file=fich,status='old')
*
* --- Loading of the series in the arrays itpw(5000),abc(3,5000) -----
*                                    and jfi(12,5000),sc(2,5000) -----
*
      nbterm=0
11    continue
      read(10,1000,end=12) ia,ib,ic,k,i,(ifi(j),j=1,12),s,c,aa,bb,cc
      nbterm=nbterm+1
      itpw(nbterm)=k
      abc(1,nbterm)=aa
      abc(2,nbterm)=bb
      abc(3,nbterm)=cc
      do jj=1,12
        jfi(jj,nbterm)=ifi(jj)
      enddo
      sc(1,nbterm)=s
      sc(2,nbterm)=c
      go to 11
12    continue
      if(nbterm.gt.5000) stop 5000
*
      write(*,1001)
      write(*,1002) 
      do j=1,10
        write(*,1003) ia,ib,ic,j,(jfi(jj,j),jj=1,12),
     .        sc(1,j),itpw(j),sc(2,j),itpw(j)
      enddo
*
      write(*,1004)
      do j=1,10
        if(abc(3,j).lt.0.d0) then
          ccc=-abc(3,j)
          write(*,1005) ia,ib,ic,j,itpw(j),(abc(l,j),l=1,2),ccc
        else
          write(*,1006) ia,ib,ic,j,itpw(j),(abc(l,j),l=1,3)
        endif
      enddo
*
* --- Substitution of the time in the series -------------------------
*
      write(*,1007)
      do it=0,ndat
        tjd=t0+step*it
        call smart97(tjd,abc,itpw,nbterm,result,ierr)
        if(ierr.ne.0) then
          write (*,1008) tjd
          stop
        endif
        write(*,1009) tjd,result
      enddo
*
      stop
*
* --- Formats -------------------------------------------------------------
*
1000  format(3i1,i2,i5,12i3,f15.4,f24.4,f24.4,f13.10,f18.8)
1001  format(/,2x,'PROGRAM EXAMPLE SMART97',
     .       //,'Printing of the first 10 terms of the',
     .          ' series p in microarcseconds, first form',/)
1002  format(27x,'ARGUMENT',24x,'SIN',26x,'COS')
1003  format(3i1,i5,2x,3(1x,4i3),f15.4,'*t**',i1,f24.4,'*t**',i1)
1004  format(//,'Printing of the first 10 terms of the',
     .          ' series p in microarcseconds, second form',/)
1005  format(3i1,i5,'  t**',i2,' *',f24.4,' * cos(',f13.10,
     .    ' -',f17.8,' *t)')
1006  format(3i1,i5,'  t**',i2,' *',f24.4,' * cos(',f13.10,
     .    ' +',f17.8,' *t)')
1007  format(//,'SMART97, Axe de figure, dynamical reference frame', 
     .              ', variable p',//,9x,'Julian date',9x,
     .              'result (microarcsecond)',/)
1008  format(/,1x,'date JD',f11.2,' out of the interval 1900-2100')
1009  format(7x,f14.5,f27.4)
*
      END
*                                                            
      subroutine SMART97(tjd,abc,itpw,nbterm,result,ierr)
*
* --- Object --------------------------------------------------------------
*
* substitution of the time in the precession-nutation series SMART97
*
* --- Input ---------------------------------------------------------------
*
*     tjd         : Julian date TT (real)
*     abc(3,5000) : coefficients a, b, c of the series (real)
*     itpw(5000)  : degree of time in the series (integer)
*     nbterm      : number of terms of the series (integer)
*
* --- Output --------------------------------------------------------------
* 
*     result      : result in microarcseconds (real)
*     ierr        : error index (integer)
*                   ierr = 0 : no error,
*                   ierr = 1 : error, date out of the interval 1900-2100.
*
* --- Declarations --------------------------------------------------------
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*
      dimension abc(3,5000),itpw(5000)
*
* --- Check up the date ---------------------------------------------------
*
      tt=(tjd-2451545.d0)/365250.d0
*
      if(abs(tt).ge.0.101d0) then
        ierr=1
        return
      endif
*
* --- Substitution of time ------------------------------------------------
*
      result=0.d0
      tpw=1.d0
*
      do nb=1,nbterm
        if(itpw(nb).ge.1) tpw=tt**itpw(nb)
        result=result+tpw*abc(1,nb)*cos(abc(2,nb)+abc(3,nb)*tt)
      enddo
*
      ierr=0
      return
*
      END
