!>
!> @file bvalue.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>
!>
function bvalue ( t, bcoef, n, k, x, jderiv )

!*************************************************************************
!
!! BVALUE evaluates a derivative of a spline from its B-spline representation.
!
!  Discussion:
!
!    The spline is taken to be continuous from the right.
!
!    The nontrivial knot interval (T(I),T(I+1)) containing X is
!    located with the aid of INTERV.  The K B-spline coefficients
!    of F relevant for this interval are then obtained from BCOEF,
!    or are taken to be zero if not explicitly available, and are
!    then differenced JDERIV times to obtain the B-spline
!    coefficients of (D**JDERIV)F relevant for that interval.
!
!    Precisely, with J = JDERIV, we have from X.(12) of the text that:
!
!      (D**J)F = sum ( BCOEF(.,J)*B(.,K-J,T) )
!
!    where
!                      / BCOEF(.),                    ,  J == 0
!                     /
!       BCOEF(.,J) = / BCOEF(.,J-1) - BCOEF(.-1,J-1)
!                   / -----------------------------,  0 < J
!                  /    (T(.+K-J) - T(.))/(K-J)
!
!    Then, we use repeatedly the fact that
!
!      sum ( A(.)*B(.,M,T)(X) ) = sum ( A(.,X)*B(.,M-1,T)(X) )
!
!    with
!                   (X - T(.))*A(.) + (T(.+M-1) - X)*A(.-1)
!      A(.,X) =   ---------------------------------------
!                   (X - T(.))      + (T(.+M-1) - X)
!
!    to write (D**J)F(X) eventually as a linear combination of
!    B-splines of order 1, and the coefficient for B(I,1,T)(X)
!    must then be the desired number (D**J)F(X).
!    See x.(17)-(19) of text.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) T(N+K), the knot sequence.  T is assumed
!    to be nondecreasing.
!
!    Input, real ( kind = 8 ) BCOEF(N), B-spline coefficient sequence.
!
!    Input, integer N, the length of BCOEF.
!
!    Input, integer K, the order of the spline.
!
!    Input, real ( kind = 8 ) X, the point at which to evaluate.
!
!    Input, integer JDERIV, the order of the derivative to
!    be evaluated.  JDERIV is assumed to be zero or positive.
!
!    Output, real ( kind = 8 ) BVALUE, the value of the (JDERIV)-th
!    derivative of the spline at X.
!
  implicit none

  integer k
  integer n

  real ( kind = 8 ) aj(k)
  real ( kind = 8 ) bcoef(n)
  real ( kind = 8 ) bvalue
  real ( kind = 8 ) dl(k)
  real ( kind = 8 ) dr(k)
  integer i
  integer ilo
  integer j
  integer jc
  integer jcmax
  integer jcmin
  integer jderiv
  integer jj
  integer mflag
  real ( kind = 8 ) t(n+k)
  real ( kind = 8 ) x

  bvalue = 0.0D+00

  if ( k <= jderiv ) then
    return
  end if
!
!  Find I so that 1 <= i < n+k and t(i) < t(i+1) and t(i) <= x < t(i+1).
!
!  If no such i can be found, X lies
!  outside the support of  the spline F and  bvalue=0.
!  (the asymmetry in this choice of  i  makes F rightcontinuous)
!
  call interv ( t, n+k, x, i, mflag )

  if ( mflag /= 0 ) then
    return
  end if
!
!  If K=1 (and jderiv = 0), bvalue = bcoef(i).
!
  if ( k <= 1 ) then
    bvalue = bcoef(i)
    return
  end if
!
!  Store the K b-spline coefficients relevant for the knot interval
!  (T(i),T(i+1)) in aj(1),...,aj(k) and compute dl(j)=x-t(i+1-j),
!  dr(j)=T(i+j)-x, j=1,...,k-1 . set any of the aj not obtainable
!  from input to zero.  Set any T's not obtainable equal to T(1) or
!  to T(n+k) appropriately.
!
  jcmin = 1

  if ( k <= i ) then

    do j = 1, k-1
      dl(j) = x-t(i+1-j)
    end do

  else

    jcmin = 1-(i-k)

    do j = 1, i
      dl(j) = x-t(i+1-j)
    end do

    do j = i, k-1
      aj(k-j) = 0.0D+00
      dl(j) = dl(i)
    end do

  end if

  jcmax = k

  if ( i <= n ) then
    go to 90
  end if

  jcmax = k + n - i
  do j = 1, k+n-i
    dr(j) = t(i+j)-x
  end do

  do j = k+n-i, k-1
    aj(j+1) = 0.0D+00
    dr(j) = dr(k+n-i)
  end do

  go to 110

   90 continue

  do j = 1, k-1
    dr(j) = t(i+j)-x
  end do

  110 continue

  do jc = jcmin, jcmax
    aj(jc) = bcoef(i-k+jc)
  end do
!
!  Difference the coefficients JDERIV times.
!
  do j = 1, jderiv

    ilo = k-j
    do jj = 1, k-j
      aj(jj) = ((aj(jj+1)-aj(jj))/(dl(ilo)+dr(jj))) * real ( k - j, kind = 8 )
      ilo = ilo-1
    end do

  end do
!
!  Compute value at X in (t(i),t(i+1)) of jderiv-th derivative,
!  given its relevant b-spline coefficients in aj(1),...,aj(k-jderiv).
!
  do j = jderiv+1, k-1
    ilo = k-j
    do jj = 1, k-j
      aj(jj) = ( aj(jj+1) * dl(ilo) + aj(jj) * dr(jj) ) &
        / ( dl(ilo) + dr(jj) )
      ilo = ilo-1
    end do
  end do

  bvalue = aj(1)

  return
end