!>
!> @file bsplpp.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/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine bsplpp ( t, bcoef, n, k, scrtch, break, coef, l )

!*************************************************************************
!
!! BSPLPP converts from B-spline to piecewise polynomial form.
!
!  Discussion:
!
!    The B-spline representation of a spline is ( T, BCOEF, N, K ),
!    while the piecewise polynomial representation is
!    ( BREAK, COEF, L, K ).
!
!    For each breakpoint interval, the K relevant B-spline coefficients
!    of the spline are found and then differenced repeatedly to get the
!    B-spline coefficients of all the derivatives of the spline on that
!    interval.
!
!    The spline and its first K-1 derivatives are then evaluated at the
!    left end point of that interval, using BSPLVB repeatedly to obtain
!    the values of all B-splines of the appropriate order at that point.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) T(N+K), the knot sequence.
!
!    Input, real ( kind = 8 ) BCOEF(N), the B spline coefficient sequence.
!
!    Input, integer N, the number of B spline coefficients.
!
!    Input, integer K, the order of the spline.
!
!    Work array, real ( kind = 8 ) SCRTCH(K,K).
!
!    Output, real ( kind = 8 ) BREAK(L+1), the piecewise polynomial breakpoint
!    sequence.  BREAK contains the distinct points in the
!    sequence T(K),...,T(N+1)
!
!    Output, real ( kind = 8 ) COEF(K,N), with COEF(I,J) = (I-1)st derivative
!    of the spline at BREAK(J) from the right.
!
!    Output, integer L, the number of polynomial pieces which
!    make up the spline in the interval (T(K),T(N+1)).
!
  implicit none

  integer k
  integer l
  integer n

  real ( kind = 8 ) bcoef(n)
  real ( kind = 8 ) biatx(k)
  real ( kind = 8 ) break(*)
  real ( kind = 8 ) coef(k,n)
  real ( kind = 8 ) diff
  integer i
  integer j
  integer jp1
  integer left
  integer lsofar
  real ( kind = 8 ) scrtch(k,k)
  real ( kind = 8 ) sum1
  real ( kind = 8 ) t(n+k)

  lsofar = 0
  break(1) = t(k)

  do left = k, n
!
!  Find the next nontrivial knot interval.
!
    if ( t(left+1) == t(left) ) then
      cycle
    end if

    lsofar = lsofar + 1
    break(lsofar+1) = t(left+1)

    if ( k <= 1 ) then
      coef(1,lsofar) = bcoef(left)
      cycle
    end if
!
!  Store the K B-spline coefficients relevant to current knot
!  interval in SCRTCH(*,1).
!
    do i = 1, k
      scrtch(i,1) = bcoef(left-k+i)
    end do
!
!  For j=1,...,k-1, compute the  k-j  b-spline coefficients relevant to
!  current knot interval for the j-th derivative by differencing
!  those for the (j-1)st derivative, and store in scrtch(.,j+1) .
!
    do jp1 = 2, k
      j = jp1-1
      do i = 1, k-j
        diff = t(left+i)-t(left+i-(k-j))
        if ( 0.0D+00 < diff ) then
          scrtch(i,jp1)=((scrtch(i+1,j)-scrtch(i,j)) / diff ) &
            * real ( k - j, kind = 8 )
        end if
      end do
    end do
!
!  For J=0, ..., K-1, find the values at T(left)  of the  j+1
!  B-splines of order J+1 whose support contains the current
!  knot interval from those of order J (in  biatx ), then comb-
!  ine with the B-spline coefficients (in scrtch(.,k-j) ) found earlier
!  to compute the (k-j-1)st derivative at  t(left)  of the given
!  spline.
!
    call bsplvb ( t, 1, 1, t(left), left, biatx )

    coef(k,lsofar) = scrtch(1,k)

    do jp1 = 2, k

      call bsplvb ( t, jp1, 2, t(left), left, biatx )

      sum1 = 0.0D+00
      do i = 1, jp1
        sum1 = sum1 + biatx(i) * scrtch(i,k+1-jp1)
      end do

      coef(k+1-jp1,lsofar) = sum1

    end do

  end do

  l = lsofar

  return
end