!>
!> @file spline_hermite_val.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 spline_hermite_val ( ndata, tdata, c, tval, sval )

!*************************************************************************
!
!! SPLINE_HERMITE_VAL evaluates a piecewise cubic Hermite interpolant.
!
!  Discussion:
!
!    SPLINE_HERMITE_SET must be called first, to set up the
!    spline data from the raw function and derivative data.
!
!  Modified:
!
!    06 April 1999
!
!  Reference:
!
!    Conte and de Boor,
!    Algorithm PCUBIC,
!    Elementary Numerical Analysis,
!    1973, page 234.
!
!  Parameters:
!
!    Input, integer NDATA, the number of data points.
!    NDATA is assumed to be at least 2.
!
!    Input, real ( kind = 8 ) TDATA(NDATA), the abscissas of the data points.
!    The entries of TDATA are assumed to be strictly increasing.
!
!    Input, real ( kind = 8 ) C(4,NDATA), contains the data computed by
!    SPLINE_HERMITE_SET.
!
!    Input, real ( kind = 8 ) TVAL, the point where the interpolant is to
!    be evaluated.
!
!    Output, real ( kind = 8 ) SVAL, the value of the interpolant at TVAL.
!
  implicit none

  integer ndata

  real ( kind = 8 ) c(4,ndata)
  real ( kind = 8 ) dt
  integer i
  integer j
  real ( kind = 8 ) sval
  real ( kind = 8 ) tdata(ndata)
  real ( kind = 8 ) tval
!
!  Find the interval J = [ TDATA(J), TDATA(J+1) ] that contains
!  or is nearest to TVAL.
!
  j = ndata - 1

  do i = 1, ndata-2

    if ( tval < tdata(i+1) ) then
      j = i
      exit
    end if

  end do
!
!  Evaluate the cubic polynomial.
!
  dt = tval - tdata(j)

  sval = c(1,j) + dt * ( c(2,j) + dt * ( c(3,j) + dt * c(4,j) ) )

  return
end