!>
!> @file ppvalu.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine ppvalu ( break, coef, l, k, x, jderiv, value )

!*******************************************************************************
!
!! PPVALU evaluates a piecewise polynomial function or its derivative.
!
!  Discussion:
!
!    PPVALU calculates the value at X of the JDERIV-th derivative of
!    the piecewise polynomial function F from its piecewise
!    polynomial representation.
!
!    The interval index I, appropriate for X, is found through a
!    call to INTERV.  The formula for the JDERIV-th derivative
!    of F is then evaluated by nested multiplication.
!
!    The J-th derivative of F is given by:
!
!      (d**j)f(x) =
!        coef(j+1,i) + h * (
!        coef(j+2,i) + h * (
!        ...
!        coef(k-1,i) + h * (
!        coef(k,i) / (k-j-1) ) / (k-j-2) ... ) / 2 ) / 1
!
!    with
!
!      H=X-BREAK(I)
!
!    and
!
!      i = max( 1 , max( j ,  break(j) <= x , 1 <= j <= l ) ).
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) BREAK(L+1), real COEF(*), integer L, for
!    piecewise polynomial representation of the function F to
!    be evaluated.
!
!    Input, integer K, the order of the polynomial pieces
!    that make up the function F.  The usual value for
!    K is 4, signifying a piecewise cubic polynomial.
!
!    Input, real ( kind = 8 ) X, the point at which to evaluate F or
!    of its derivatives.
!
!    Input, integer JDERIV, the order of the derivative to be
!    evaluated.  If JDERIV is 0, then F itself is evaluated,
!    which is actually the most common case.  It is assumed
!    that JDERIV is zero or positive.
!
!    Output, real ( kind = 8 ) VALUE, the value of the JDERIV-th
!    derivative of F at X.
!
  implicit none

  integer k
  integer l

  real ( kind = 8 ) break(l+1)
  real ( kind = 8 ) coef(k,l)
  real ( kind = 8 ) fmmjdr
  real ( kind = 8 ) h
  integer i
  integer jderiv
  integer m
  integer ndummy
  real ( kind = 8 ) value
  real ( kind = 8 ) x

  value = 0.0D+00

  fmmjdr = k - jderiv
!
!  Derivatives of order K or higher are identically zero.
!
  if ( k <= jderiv ) then
    return
  end if
!
!  Find the index I of the largest breakpoint to the left of X.
!
  call interv ( break, l+1, x, i, ndummy )
!
!  Evaluate the JDERIV-th derivative of the I-th polynomial piece at X.
!
  h = x - break(i)
  m = k

  do

    value = ( value / fmmjdr ) * h + coef(m,i)
    m = m - 1
    fmmjdr = fmmjdr - 1.0D+00

    if ( fmmjdr <= 0.0D+00 ) then
      exit
    end if

  end do

  return
end