program spline_pgplot_main


      implicit none
      integer max_tabu, ngrids_max, n_tabu, i, choice
      integer more_points, ngrids, pgopen
      parameter (max_tabu = 100, ngrids_max = 2000)
      real xa(max_tabu), ya(max_tabu), y2a(max_tabu), x, y, delta_x
      real x_plot(ngrids_max), y_plot(ngrids_max), yp1, ypn
      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 slopes at two ends. (If > 10E30 then natural spline is set.)
      write(*,*) 'Please type in slopes at bothend (> 10E30 -> y2a=0.'
      read (*,*) yp1, ypn

c Call subroutine spline to solve for y2a array.
      call spline(xa, ya, n_tabu, yp1, ypn, y2a)

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.
      if (choice.eq.0) then
1002    write(*,*) 'Type in your x:'
        read(*,*) x
        call splint(xa, ya, y2a, n_tabu, x, y)
        write(*,*) 'Interpolated value is ', y
        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.
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)
        do i=1, ngrids
          x_plot(i) = x_min + (i-1)*delta_x
          call splint(xa,ya,y2a,n_tabu,x_plot(i),y_plot(i))
        end do

      if ( pgopen('?') .le. 0 ) stop
      call pgpap(8.0,0.75)
      call pgenv(0.0, 4.0, 0.0, 20.0, 0, 1)
      call pgline(ngrids, x_plot, y_plot)
      call pgpt(n_tabu, xa, ya, 9)
      call pgclos

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

      end if






      end