!>
!> @file cubslo.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 cubslo ( tau, c, n )

!*******************************************************************************
!
!! CUBSLO solves for slopes defining a cubic spline.
!
!  Discussion:
!
!    A tridiagonal linear system for the unknown slopes S(I) of
!    F at TAU(I), I=1,..., N, is generated and then solved by Gauss
!    elimination, with S(I) ending up in C(2,I), for all I.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) TAU(N), the abscissas or X values of
!    the data points.  The entries of TAU are assumed to be
!    strictly increasing.
!
!    Input, integer N, the number of data points.  N is
!    assumed to be at least 2.
!
!    Input/output, real ( kind = 8 ) C(4,N).
!    On input, C(1,I) contains the function value at TAU(I),
!    for I = 1 to N.
!    C(2,1) contains the slope at TAU(1) and C(2,N) contains
!    the slope at TAU(N).
!    On output, the intermediate slopes at TAU(I) have been
!    stored in C(2,I), for I = 2 to N-1.
!
  implicit none

  integer n

  real ( kind = 8 ) c(4,n)
  real ( kind = 8 ) g
  integer i
  integer ibcbeg
  integer ibcend
  real ( kind = 8 ) tau(n)
!
!  Set up the right hand side of the linear system.
!  C(2,1) and C(2,N) are presumably already set.
!
  do i = 2, n-1
    c(2,i) = 3.0D+00 * ( &
      ( tau(i) - tau(i-1) ) * ( c(1,i+1) - c(1,i) ) / ( tau(i+1) - tau(i) ) + &
      ( tau(i+1) - tau(i) ) * ( c(1,i) - c(1,i-1) ) / ( tau(i) - tau(i-1) ) )
  end do
!
!  Set the diagonal coefficients.
!
  c(4,1) = 1.0D+00
  do i = 2, n-1
    c(4,i) = 2.0D+00 * ( tau(i+1) - tau(i-1) )
  end do
  c(4,n) = 1.0D+00
!
!  Set the off-diagonal coefficients.
!
  c(3,1) = 0.0D+00
  do i = 2, n
    c(3,i) = tau(i) - tau(i-1)
  end do
!
!  Forward elimination.
!
  do i = 2, n-1
    g = -c(3,i+1) / c(4,i-1)
    c(4,i) = c(4,i) + g * c(3,i-1)
    c(2,i) = c(2,i) + g * c(2,i-1)
  end do
!
!  Back substitution for the interior slopes.
!
  do i = n-1, 2, -1
    c(2,i) = ( c(2,i) - c(3,i) * c(2,i+1) ) / c(4,i)
  end do

  return
end