! test_laphi.f90  family of Lagrange functions and derivatives
!             L_n,j(x) = product i=1,n i/=j (x-x_i)/(x_j-x_i)
!             for points x_1, x_2, ,. , x_n
!  generate L_n,j(x) for n=2,6  j=1,n and derivatives to 4th
!  should be same results as produced by test_laphi.adb


program test_laphi
  use laphi
  implicit none
  integer, parameter :: NN = 6
  double precision :: h, hf
  double precision, dimension(NN) :: xg
  double precision :: x, y, yp, ypp, yppp, ypppp
  double precision :: xmin = 0.0 ! xg(0)
  double precision :: xmax = 1.0 ! xg(n)
  integer :: n, j, k

  print *, 'test_laphi.f90 running on 0<x<1 for N=',NN
  print *, ' L_order,point  order is highest power of x, point 1,n'
  print *, '    L_n,j(x) = product i=1,n i/=j (x-x_i)/(x_j-x_i)'
  print *, '    for points x_1, x_2, ,. , x_n'
  print *, '    generate L_n,j(x) for n=2,6  j=1,n and derivatives to 4th'

  do n=2,NN  ! points go 1 to n
    h = (xmax-xmin)/(n-1)
    do j=1,n 
      xg(j) = xmin+(j-1)*h
      print *, 'xg(',j,')=',xg(j)
    end do
    hf = h/4.0
    do j=1,n 
      do k=1,4*(n-1)+1 
        x = xmin+(k-1)*hf

        ! evaluate function, n,j,x,xg,denom
        y = phi(x, j, n, xg)

        ! evaluate first derivative of function, n,j,x,xg,denom
        yp = phip(x, j, n, xg)

        ! evaluate second derivative of function, n,j,x,xg,denom
        ypp = phipp(x, j, n, xg)

        ! evaluate third derivative of function, n,j,x,xg,denom
        yppp = phippp(x, j, n, xg)

        ! evaluate fourth derivative of function, n,j,x,xg,denom
        ypppp = phipppp(x, j, n, xg)

        print "('L_',I1,',',I1,'(',f6.3,')=',f7.3,', yp=',f6.2,', ypp=',f7.2,', yppp=',f7.1,', ypppp=',f7.1)",  &
               n, j, x, y, yp, ypp, yppp, ypppp
      end do ! on k,x
      print *, ' '
    end do ! on j
    print *, ' '
  end do ! on n
  print *, 'end test_laphi.f90'
end program test_laphi