program polint_main

c In subroutine polint, array declaration is used for xa(n) and ya(n),
c for this we need to reclare larger or equal size xa(n) and ya(n) in
c this main program. As for NMAX and c(NAMX) and d(NAMX) in that routine,
c they are taken care locally in that routine alreay with specific size
c (parameter) given when declaring array, therefore we don't need to deal
c with that.

      implicit none
      integer max_tabu, ngrids_max, n_tabu, m, i, j, k, choice
      integer k1, k2, more_points, ngrids
      parameter (max_tabu = 100, ngrids_max = 2000)
      real xa(max_tabu), ya(max_tabu), x, y, dy, delta_x
      real x_plot(ngrids_max), y_plot(ngrids_max),dy_plot(ngrids_max)
      real x_min, x_max 
      character*40 filename1

c Ask filename and the number of tabulated points in the file.
      write(*,*) 'Type in the file name of tabulated values (x,y)'
      read (*,*) filename1
      write(*,*) 'How many points are there in the tabulated set ?'
      read (*,*) n_tabu
      if (n_tabu.gt.max_tabu) stop 'Too many tabulated values,
     &please recompile the main program with a larger max_tabu.'
      
c Open the file and read in the tabulated values, and write to screen
c for double check.      
      open (unit = 20, file = filename1)

      do i = 1, n_tabu
        read (20,*) xa(i), ya(i)
        write(*,*) xa(i), ya(i)
      end do

      close (20)

c Ask the order of polynomial to be usedfor interpolation, not higher
c than 10.
1000  continue
      write(*,*) 'How many points do you want to form polynomial ?'
      write(*,*) 'For eample, 2 points mean linear interpolation.'
      read (*,*) m

      if (m.gt.10) then
        write(*,*) 'Waring : Neumarial Recipe default max is 10.'
        write(*,*) 'Please decress your value or set a bigger NMAX'
        write(*,*) 'in the subroutine polint.'
        write(*,*) 'Your choice (0 : exit) or (1: input again).'
        read (*,*) choice
        if (choice.eq.0) stop 'Terminaed by user.'
        if (choice.eq.1) goto 1000
        goto 1000
      endif

c One can choose to type in x one by one or to generte a continuous set
c for graph plotting.
      write(*,*) 'What do you want to do ?'
      write(*,*)'0:input x, find y ; 1:generate a set of points (x,y)'
      read (*,*) choice

c This part get x one by one, using "locate" to find the right interval.
      if (choice.eq.0) then
1002    write(*,*) 'Type in your x:'
        read(*,*) x
        call locate (xa, n_tabu, x, j)
        k = min( max(j-(m-1)/2,1) , n_tabu+1-m )
        call polint(xa(k), ya(k), m, x, y, dy)
        write(*,*) 'Interpolated value is ', y, ', error is ', dy
        write(*,*) 'One more point (yes=1;no=0) ?'
        read (*,*) more_points
        if (more_points.eq.1) goto 1002
      endif  

c This part generate a set of x, using "hunt" to find the right interval.
c The resulting (x,y) is written out into a file named by user.
      if (choice.eq.1) then

        write(*,*) 'Type your range x_initial and x-final:'
        read (*,*) x_min, x_max
        write(*,*) 'How many points do you want to plot ?'
        read (*,*) ngrids
        if(ngrids.le.1) stop 'Need at least two grids.'

        delta_x = ( x_max - x_min ) / (ngrids - 1)
        j = 1
        do i=1, ngrids
          x = x_min + (i-1)*delta_x
          call hunt(xa, n_tabu, x, j)
          k = min( max(j-(m-1)/2,1) , n_tabu+1-m )
          x_plot(i) = x
          write(*,*) 'x is ',x_plot(i),' j is ',j,' k is',k
          call polint(xa(k),ya(k),m,x_plot(i),y_plot(i),dy_plot(i))
        end do

        write(*,*) 'Plrease give a name to your output file:'
        read (*,*) filename1
        open (unit=30, file=filename1)
        do i=1,ngrids
          write(30,*) x_plot(i), y_plot(i)
        end do
        close(30)

        write(*,*) 'Info: Error estimate will be written into file
     &error.dat'
        open (unit=40, file='error.dat')
        do i=1,ngrids
          write(40,*) x_plot(i), dy_plot(i)
        end do
        close(40)

      end if

      end