Program Interp
c===================================================================
c Driver program for interpolation with fintr.f and spline.f/seval/f
c September 2006
c The program does following
c 1. read inital data for preparing an interpolated data table
c 2. generate xi() and yi() table for interpolation later
c 3. does spline interpolation
c 4. does polinomial interpolation (1st, 3rd, 5th and 7th order)
c-------------------------------------------------------------------
c n=24 - max nuber of points for the initial table
c nxi  - real number of points for the initial table
c nout - number of output points (the final table)
c key  = 0 equidistant mesh for interpolation (xmax-xmin)/(nxi-1)
c      = 1 reading xi(i) from the input stream
c xmin - left point for the interval
c xmax - right point
c f(w) - the function for interpolation
c
c files:
c int_i.dat  - initial information for interpolation
c int_o1.dat - generated data table forinterpolation
c int_o2.dat - results of interpolation
c===================================================================
      IMPLICIT REAL*8 (A-H,O-Z),INTEGER * 4(I-N)
      Parameter (n=24)
      Dimension xi(n), yi(n), b(n), c(n), d(n)

c      f(w) = sin(w**2)
c      f(w) = exp(sin(w))
c      f(w) = 0.200/( (w-3.2)**2 + 0.04)
       f(w) = 1.0/(1.0+25*x*x)
c      f(w) = cos(3.0*w)/(0.4+(w-2.0)**2)


c open files
      open (unit=1, file='inter_i1.dat')
      open (unit=2, file='inter_o1.dat')
      open (unit=3, file='inter_o2.dat')
c
c reading data
c
      read (1,1000) nxi, nout, key
      write(*,1000) nxi, nout, key
      if (nxi.gt.n) then
        write (*,2001)
        stop
      end if
2001  format(' too many point for the table')
      read (1,1001) xmin, xmax
      write(*,1001) xmin, xmax
c
c generation or reading x_i array
c
      if(key.eq.0) then
          do i=1,nxi
             xi(i)=xmin+real(i-1)*(xmax-xmin)/real(nxi-1)
          end do
          else
             read (1,1001) (xi(i),i=1,n)
      end if
      write(*,1001) (xi(i),i=1,nxi)
c
c generating the table of f(x_i) for x_i
c
      do i=1,nxi
         x = xi(i)
         yi(i) = f(x)
         write(2,1002) x, yi(i)
      end do
c
c preparation
c
      call spline(nxi,xi,yi,b,c,d)
      step = (xmax-xmin)/real(nout-1)
c
c interpolation
c
      write (3,1004)
      do i=1,nout
        x = xmin+step*real(i-1)
        y = f(x)
        y1 = fintr (x,yi,xi,nxi,2)
        if (nxi.ge.4) y3 = fintr (x,yi,xi,nxi,4)
        if (nxi.ge.6) y5 = fintr (x,yi,xi,nxi,6)
        if (nxi.ge.8) y7 = fintr (x,yi,xi,nxi,8)
        ys = seval (nxi,x,xi,yi,b,c,d)
        write (3,1002) x, y, y1, y3, y5, y7, ys
      end do

1000   format (3i5)
1001   format (12f5.2)
1002   format (f5.2,6f9.4)
1003   format (6f9.4)
1004   format ('  x     f',8x,'p1       p3       p5       p7    spline')
      stop
        end