function [ yval, ypval, yppval ] = spline_cubic_val ( n, t, y, ypp, tval )

%*****************************************************************************80
%
%% SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point.
%
%  Discussion:
%
%    SPLINE_CUBIC_SET must have already been called to define the
%    values of YPP.
%
%    For any point T in the interval T(IVAL), T(IVAL+1), the form of
%    the spline is
%
%      SPL(T) = A
%             + B * ( T - T(IVAL) )
%             + C * ( T - T(IVAL) )**2
%             + D * ( T - T(IVAL) )**3
%
%    Here:
%      A = Y(IVAL)
%      B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
%        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
%      C = YPP(IVAL) / 2
%      D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl de Boor,
%    A Practical Guide to Splines,
%    Springer Verlag, 1978.
%
%  Parameters:
%
%    Input, integer N, the number of data values.
%
%    Input, real T(N), the knot values.
%
%    Input, real Y(N), the data values at the knots.
%
%    Input, real YPP(N), the second derivatives of the spline at the knots.
%
%    Input, real TVAL, a point, typically between T(1) and T(N), at
%    which the spline is to be evalulated.  If TVAL lies outside
%    this range, extrapolation is used.
%
%    Output, real YVAL, YPVAL, YPPVAL, the value of the spline, and
%    its first two derivatives at TVAL.
%

%
%  Determine the interval [T(LEFT), T(RIGHT)] that contains TVAL.
%  Values below T(1) or above T(N) use extrapolation.
%
  [ left, right ] = r8vec_bracket ( n, t, tval );
%
%  Evaluate the polynomial.
%
  dt = tval - t(left);
  h = t(right) - t(left);

  yval = y(left) ...
       + dt * ( ( y(right) - y(left) ) / h ...
              - ( ypp(right) / 6.0 + ypp(left) / 3.0 ) * h ...
       + dt * ( 0.5 * ypp(left) ...
       + dt * ( ( ypp(right) - ypp(left) ) / ( 6.0 * h ) ) ) );

  ypval = ( y(right) - y(left) ) / h ...
       - ( ypp(right) / 6.0 + ypp(left) / 3.0 ) * h ...
       + dt * ( ypp(left) ...
       + dt * ( 0.5 * ( ypp(right) - ypp(left) ) / h ) );

  yppval = ypp(left) + dt * ( ypp(right) - ypp(left) ) / h;

  return
end