!  Simple tests for cspline.f90.  Uses programs from lineqn.f90.
!  Source file:  csptst.f90

!  Test N-point spline

      subroutine testspline(N)
      integer N

      INCLUDE 'lineqn.hp'
      INCLUDE 'cspline.hI'

      real,    dimension(1:Nx) :: X, Xs, Ys, Y2, Z, U, W, R
      integer, dimension(1:Nx) :: I			! Coordinates
      real h, xa, xb, da, db, f
      integer ia, ib

      I = [1:Nx]
      print *, "Cubic splines using ln(x) in the interval [xa:xb]."
      print *, "End conditions are i?-th derivative at x? = d?.?"
      print *, "Specify xa <= 0 to exit, f < 0 to try new data."
      do while ( .true. )
	print '(" xa ia da > "$)'
	read *, xa, ia, da
	if ( xa <= 0.0 ) go to 99		! Exit
	print '(" xb ib db > "$)'
	read *, xb, ib, db
	h  = (xb - xa)/(N - 1)
	Xs = xa + (I - 1)*h
	Ys = log(Xs)
	Y2 = splined2(Xs, Ys, N, ia, da, ib, db)
	print *, "2nd derivatives:"
	print *, "Range: ", minval(Y2), maxval(Y2)
	print 21
21	FORMAT ("   f   min error   max error  max rel err")
        do while ( .true.)
	  print '(" f > "$)'
	  read *, f
	  if ( f < 0.0 ) go to 88		! To read new data
	  X = Xs + h*f
	  Z = csplint(X, Xs, Ys, Y2, N)
	  U = log(X)
	  W = U - Z				! Error
	  R = abs(W/U)				! Rel. error
	  where ( U == 0.0 ) R = abs(Z)
	  print 22, f, minval(W(1:N-1)), maxval(W(1:N-1)), &
      			 maxval(R(1:N-1))
 22	  FORMAT(X, F4.2, 3(X, E11.4))
	end do
 88	continue				! with new data
      end do
 99   continue					! For exit
      end subroutine testspline


!  Read problem size and test.

      program tspline
      INCLUDE 'lineqn.hp'
      integer N

      print *, "N = problem size, <= ", Nx, ".  Specify N <= 0 to quit."
      do while ( .true. )
	print '(" N > "$)'
	read *, N
	if ( N <= 0 ) go to 99			! Exit
	call testspline(N)
      end do
 99   continue					! For exit
      end program tspline
