! nuderiv_test.f90  non uniformly spaced derivative coefficients test 

program nuderiv_test
  implicit none
  integer ::  order, npoint, point
  double precision, dimension(5) :: x 
  data x / 1.0, 2.0, 4.0, 7.0, 8.5 /
  double precision, dimension(5) :: y
  double precision, dimension(5) :: yp
  double precision, dimension(5) :: ypp
  double precision, dimension(5) :: yppp
  double precision, dimension(5) :: ypppp
  double precision, dimension(5) :: c
  double precision :: nuyp ! any order 
  integer :: i, j

  interface
    subroutine nuderiv(order, npoint, point, X, C)
      integer, intent(in) :: order           ! derivative order
      integer, intent(in) :: npoint          ! number of points
      integer, intent(in) :: point           ! where derivative computed
      double precision, dimension(npoint), intent(in) :: X ! X input values
      double precision, dimension(npoint), intent(inout) :: C ! coef of Y's
    end subroutine nuderiv
    function yf4(x) result(value) ! test function 4 
      double precision, intent(in) :: x
      double precision :: value
    end function yf4
    function ypf4(x) result(value) ! derivative of test function 4
      double precision, intent(in) :: x
      double precision :: value
    end function ypf4
    function yppf4(x) result(value) ! second derivative of test function 4 
      double precision, intent(in) :: x
      double precision :: value
    end function yppf4
    function ypppf4(x) result(value) ! third derivative of test function 4 
      double precision, intent(in) :: x
      double precision :: value
    end function ypppf4
    function yf5(x) result(value) ! test function 5 
      double precision, intent(in) :: x
      double precision :: value
    end function yf5
    function ypf5(x) result(value) ! derivative of test function 5
      double precision, intent(in) :: x
      double precision :: value
    end function ypf5
    function yppf5(x) result(value) ! second derivative of test function 5 
      double precision, intent(in) :: x
      double precision :: value
    end function yppf5
    function ypppf5(x) result(value) ! third derivative of test function 5 
      double precision, intent(in) :: x
      double precision :: value
    end function ypppf5
    function yppppf5(x) result(value) ! fourth derivative of test function 5 
      double precision, intent(in) :: x
      double precision :: value
    end function yppppf5
  end interface

  print *, 'nuderiv_test.f90 running'
  npoint = 5
  print *, 'testing  npoint=',npoint
  do i=1,npoint
    y(i) = yf4(x(i))
    yp(i) = ypf4(x(i))
    ypp(i) = yppf4(x(i))
    yppp(i) = ypppf4(x(i))
    print *, 'x(',i,')=',x(i),', y=',y(i),', yp=',yp(i),', ypp=',ypp(i), &
            ', yppp=',yppp(i)
  end do ! i
  do order=1,npoint-1
    do point=1,npoint
      print *, 'testing order=',order,', npoint=',npoint,', point=',point
      call nuderiv(order, npoint, point, x, c)
      do i=1,npoint
         print *, 'nu c(',i,')=',c(i)
      end do ! i
      nuyp = 0.0
      do i=1,npoint
        nuyp = nuyp + c(i)*y(i)
      end do ! i
      if(order .eq. 1) then
        print *, 'order=',order,', npoint=',npoint,', point=',point, &
                 ', yp=',yp(point),', nuyp=',nuyp,', err=',yp(point)-nuyp
      else if(order .eq. 2) then
        print *, 'order=',order,', npoint=',npoint,', point=',point, &
                 ', ypp=',ypp(point),', nuyp=',nuyp,', err=',ypp(point)-nuyp
      else if(order .eq. 3) then
        print *, 'order=',order,', npoint=',npoint,', point=',point, &
                 ', yppp=',yppp(point),', nuyp=',nuyp,', err=',yppp(point)-nuyp
      end if
    end do ! point
  end do ! order

  npoint = 5
  print *, ' testing  npoint=',npoint
  do i=i,npoint
    y(i) = yf5(x(i))
    yp(i) = ypf5(x(i))
    ypp(i) = yppf5(x(i))
    yppp(i) = ypppf5(x(i))
    ypppp(i) = yppppf5(x(i))
    print *, 'x(%d)=',x(i),', y=',y(i),', yp=',yp(i),', ypp=',ypp(i), &
             ', yppp=',yppp(i),', ypppp=',ypppp(i)
  end do ! 
  do order=1,npoint-1
    do point=1,npoint
      print *, 'testing order=',order,', npoint=',npoint,', point=',point
      call nuderiv(order, npoint, point, x, c)
      do i=1,npoint
        print *, 'nu c(',i,')=',c(i)
      end do ! i
      nuyp = 0.0D0
      do i=1,npoint
        nuyp = nuyp + c(i)*y(i)
      end do ! i
      if(order .eq. 1) then
        print *, 'order=',order,', npoint=',npoint,', point=',point, &
                 ', yp=',yp(point),', nuyp=',nuyp,', err=',yp(point)-nuyp
      else if(order .eq. 2) then
        print *, 'order=',order,', npoint=',npoint,', point=',point, &
                 ', ypp=',ypp(point),', nuyp=',nuyp,', err=',ypp(point)-nuyp
      else if(order .eq. 3) then
        print *, 'order=',order,', npoint=',npoint,', point=',point, &
                 ', yppp=',yppp(point),', nuyp=',nuyp,', err=',yppp(point)-nuyp
      else if(order .eq. 4) then
        print *, 'order=',order,', npoint=',npoint,', point=',point, &
                 ', ypppp=',ypppp(point),', nuyp=',nuyp, &
                 ', err=',ypppp(point)-nuyp
      end if
    end do ! point
  end do ! order
end program nuderiv_test


function yf4(x) result(value) ! test function 4 
  implicit none
  double precision, intent(in) :: x
  double precision :: value
  value = 3.0 -0.5* x + 0.2*x*x - 0.1*x*x*x 
end function yf4

function ypf4(x) result(value) ! derivative of test function 4
  implicit none
  double precision, intent(in) :: x
  double precision :: value
  value = -0.5 + 0.4*x - 0.3*x*x 
end function ypf4

function yppf4(x) result(value) ! second derivative of test function 4 
  implicit none
  double precision, intent(in) :: x
  double precision :: value
  value = 0.4 - 0.6*x 
end function yppf4

function ypppf4(x) result(value) ! third derivative of test function 4 
  implicit none
  double precision, intent(in) :: x
  double precision :: value
  value = -0.6 
end function ypppf4

function yf5(x) result(value) ! test function 5 
  implicit none
  double precision, intent(in) :: x
  double precision :: value
  value = 3.0 -0.5* x + 0.1*x*x - 0.05*x*x*x + 0.002*x*x*x*x 
end function yf5

function ypf5(x) result(value) ! derivative of test function 5
  implicit none
  double precision, intent(in) :: x
  double precision :: value
  value = -0.5 + 0.2*x - 0.15*x*x + 0.008*x*x*x 
end function ypf5

function yppf5(x) result(value) ! second derivative of test function 5 
  implicit none
  double precision, intent(in) :: x
  double precision :: value
  value = 0.2 - 0.3*x + 0.024*x*x 
end function yppf5

function ypppf5(x) result(value) ! third derivative of test function 5 
  implicit none
  double precision, intent(in) :: x
  double precision :: value
  value = -0.3 + 0.048*x 
end function ypppf5

function yppppf5(x) result(value) ! fourth derivative of test function 5 
  implicit none
  double precision, intent(in) :: x
  double precision :: value
  value = 0.048 
end function yppppf5