mo_spline.f90

!> \file mo_spline.f90

!> \brief Module which defines and evaluates spline functions.

!> \details SPLINE is a FORTRAN90 library which defines and evaluates spline functions.
!>
!>          These spline functions are typically used to:
!>
!>          interpolate data exactly at a set of points;
!>          approximate data at many points, or over an interval.
!>          The most common use of this software is for situations where a set of (X,Y) data points
!>          is known, and it is desired to determine a smooth function which passes exactly through
!>          those points, and which can be evaluated everywhere. Thus, it is possible to get a formula
!>          that allows you to "connect the dots".
!>
!>          Of course, you could could just connect the dots with straight lines, but that would look ugly,
!>          and if there really is some function that explains your data, you'd expect it to curve around
!>          rather than make sudden angular turns. The functions in SPLINE offer a variety of choices for
!>          slinky curves that will make pleasing interpolants of your data.
!>
!>          There are a variety of types of approximation curves available, including:
!>              least squares polynomials,
!>              divided difference polynomials,
!>              piecewise polynomials,
!>              B splines,
!>              Bernstein splines,
!>              beta splines,
!>              Bezier splines,
!>              Hermite splines,
!>              Overhauser (or Catmull-Rom) splines.
!>
!>          Also included are a set of routines that return the local "basis matrix", which allows the
!>          evaluation of the spline in terms of local function data.

!> \authors John Burkardt, Matthias Cuntz (JAMS)
!> \date Jul 2013
module mo_spline

  ! This module is part of the JAMS Fortran library.

  ! Adapted to the CHS library by Matthias Cuntz, July 2013

  ! License
  ! -------
  ! This file is part of the JAMS Fortran package, distributed under the MIT License.
  !
  ! Copyright (c) 1999-2013, John Burkhardt, Matthias Cuntz - mc (at) macu (dot) de
  !
  ! Permission is hereby granted, free of charge, to any person obtaining a copy
  ! of this software and associated documentation files (the "Software"), to deal
  ! in the Software without restriction, including without limitation the rights
  ! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  ! copies of the Software, and to permit persons to whom the Software is
  ! furnished to do so, subject to the following conditions:
  !
  ! The above copyright notice and this permission notice shall be included in all
  ! copies or substantial portions of the Software.
  !
  ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  ! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  ! SOFTWARE.

  use mo_kind, only: dp, i4
  use mo_utils, only: eq, le, ne

  implicit none

contains

  subroutine basis_function_b_val ( tdata, tval, yval )

    !*****************************************************************************80
    !
    !! BASIS_FUNCTION_B_VAL evaluates the B spline basis function.
    !
    !  Discussion:
    !
    !    The B spline basis function is a piecewise cubic which
    !    has the properties that:
    !
    !    * it equals 2/3 at TDATA(3), 1/6 at TDATA(2) and TDATA(4);
    !    * it is 0 for TVAL <= TDATA(1) or TDATA(5) <= TVAL;
    !    * it is strictly increasing from TDATA(1) to TDATA(3),
    !      and strictly decreasing from TDATA(3) to TDATA(5);
    !    * the function and its first two derivatives are continuous
    !      at each node TDATA(I).
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    06 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Alan Davies, Philip Samuels,
    !    An Introduction to Computational Geometry for Curves and Surfaces,
    !    Clarendon Press, 1996,
    !    ISBN: 0-19-851478-6,
    !    LC: QA448.D38.
    !
    !  Parameters:
    !
    !    Input, real(dp) TDATA(5), the nodes associated with the
    !    basis function.  The entries of TDATA are assumed to be distinct
    !    and increasing.
    !
    !    Input, real(dp) TVAL, a point at which the B spline basis
    !    function is to be evaluated.
    !
    !    Output, real(dp) YVAL, the value of the function at TVAL.
    !
    implicit none

    integer(i4), parameter :: ndata = 5

    integer(i4) left
    integer(i4) right
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) u
    real(dp) yval

    if ( le(tval,tdata(1)) .or. le(tdata(ndata),tval) ) then
       yval = 0.0_dp
       return
    end if
    !
    !  Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] containing TVAL.
    !
    call r8vec_bracket ( ndata, tdata, tval, left, right )
    !
    !  U is the normalized coordinate of TVAL in this interval.
    !
    u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) )
    !
    !  Now evaluate the function.
    !
    if ( tval < tdata(2) ) then
       yval = u**3 / 6.0_dp
    else if ( tval < tdata(3) ) then
       yval = ( ( (    - 3.0_dp   &
            * u + 3.0_dp ) &
            * u + 3.0_dp ) &
            * u + 1.0_dp ) / 6.0_dp
    else if ( tval < tdata(4) ) then
       yval = ( ( (    + 3.0_dp   &
            * u - 6.0_dp ) &
            * u + 0.0_dp ) &
            * u + 4.0_dp ) / 6.0_dp
    else if ( tval < tdata(5) ) then
       yval = ( 1.0_dp - u )**3 / 6.0_dp
    end if

    return
  end subroutine basis_function_b_val
  subroutine basis_function_beta_val ( beta1, beta2, tdata, tval, yval )

    !*****************************************************************************80
    !
    !! BASIS_FUNCTION_BETA_VAL evaluates the beta spline basis function.
    !
    !  Discussion:
    !
    !    With BETA1 = 1 and BETA2 = 0, the beta spline basis function
    !    equals the B spline basis function.
    !
    !    With BETA1 large, and BETA2 = 0, the beta spline basis function
    !    skews to the right, that is, its maximum increases, and occurs
    !    to the right of the center.
    !
    !    With BETA1 = 1 and BETA2 large, the beta spline becomes more like
    !    a linear basis function; that is, its value in the outer two intervals
    !    goes to zero, and its behavior in the inner two intervals approaches
    !    a piecewise linear function that is 0 at the second node, 1 at the
    !    third, and 0 at the fourth.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    09 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Alan Davies, Philip Samuels,
    !    An Introduction to Computational Geometry for Curves and Surfaces,
    !    Clarendon Press, 1996,
    !    ISBN: 0-19-851478-6,
    !    LC: QA448.D38.
    !
    !  Parameters:
    !
    !    Input, real(dp) BETA1, the skew or bias parameter.
    !    BETA1 = 1 for no skew or bias.
    !
    !    Input, real(dp) BETA2, the tension parameter.
    !    BETA2 = 0 for no tension.
    !
    !    Input, real(dp) TDATA(5), the nodes associated with the
    !    basis function.  The entries of TDATA are assumed to be distinct
    !    and increasing.
    !
    !    Input, real(dp) TVAL, a point at which the B spline
    !    basis function is to be evaluated.
    !
    !    Output, real(dp) YVAL, the value of the function at TVAL.
    !
    implicit none

    integer(i4), parameter :: ndata = 5

    real(dp) a
    real(dp) b
    real(dp) beta1
    real(dp) beta2
    real(dp) c
    real(dp) d
    integer(i4) left
    integer(i4) right
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) u
    real(dp) yval

    if ( le(tval,tdata(1)) .or. le(tdata(ndata),tval) ) then
       yval = 0.0_dp
       return
    end if
    !
    !  Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] containing TVAL.
    !
    call r8vec_bracket ( ndata, tdata, tval, left, right )
    !
    !  U is the normalized coordinate of TVAL in this interval.
    !
    u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) )
    !
    !  Now evaluate the function.
    !
    if ( tval < tdata(2) ) then

       yval = 2.0_dp * u**3

    else if ( tval < tdata(3) ) then

       a = beta2 + 4.0_dp * beta1 + 4.0_dp * beta1 * beta1 &
            + 6.0_dp * ( 1.0_dp - beta1 * beta1 ) &
            - 3.0_dp * ( 2.0_dp + beta2 + 2.0_dp * beta1 ) &
            + 2.0_dp * ( 1.0_dp + beta2 + beta1 + beta1 * beta1 )

       b = - 6.0_dp * ( 1.0_dp - beta1 * beta1 ) &
            + 6.0_dp * ( 2.0_dp + beta2 + 2.0_dp * beta1 ) &
            - 6.0_dp * ( 1.0_dp + beta2 + beta1 + beta1 * beta1 )

       c = - 3.0_dp * ( 2.0_dp + beta2 + 2.0_dp * beta1 ) &
            + 6.0_dp * ( 1.0_dp + beta2 + beta1 + beta1 * beta1 )

       d = - 2.0_dp * ( 1.0_dp + beta2 + beta1 + beta1 * beta1 )

       yval = ( ( d * u + c ) * u + b ) * u + a

    else if ( tval < tdata(4) ) then

       a = beta2 + 4.0_dp * beta1 + 4.0_dp * beta1 * beta1

       b = - 6.0_dp * beta1 * ( 1.0_dp - beta1 * beta1 )

       c = - 3.0_dp * ( beta2 + 2.0_dp * beta1**2 + 2.0_dp * beta1**3 )

       d = 2.0_dp * ( beta2 + beta1 + beta1**2 + beta1**3 )

       yval = ( ( d * u + c ) * u + b ) * u + a

    else if ( tval < tdata(5) ) then

       yval = 2.0_dp * beta1**3 * ( 1.0_dp - u )**3

    end if

    yval = yval / ( 2.0_dp + beta2 + 4.0_dp * beta1 + 4.0_dp * beta1**2 &
         + 2.0_dp * beta1**3 )

    return
  end subroutine basis_function_beta_val
  subroutine basis_matrix_b_uni ( mbasis )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_B_UNI sets up the uniform B spline basis matrix.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    17 June 2006
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    James Foley, Andries vanDam, Steven Feiner, John Hughes,
    !    Computer Graphics, Principles and Practice,
    !    Second Edition,
    !    Addison Wesley, 1995,
    !    ISBN: 0201848406,
    !    LC: T385.C5735.
    !
    !  Parameters:
    !
    !    Output, real(dp) MBASIS(4,4), the basis matrix.
    !
    implicit none

    real(dp) mbasis(4,4)
    !
    !  In the following statement, the matrix appears as though it
    !  has been transposed.
    !
    mbasis(1:4,1:4) = real ( reshape ( &
         (/ -1,  3, -3, 1,                &
         3, -6,  0, 4,                &
         -3,  3,  3, 1,                &
         1,  0,  0, 0 /),             &
         (/ 4, 4 /) ), kind=dp ) / 6.0_dp

    return
  end subroutine basis_matrix_b_uni
  subroutine basis_matrix_beta_uni ( beta1, beta2, mbasis )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_BETA_UNI sets up the uniform beta spline basis matrix.
    !
    !  Discussion:
    !
    !    If BETA1 = 1 and BETA2 = 0, then the beta spline reduces to
    !    the B spline.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    12 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    James Foley, Andries vanDam, Steven Feiner, John Hughes,
    !    Computer Graphics, Principles and Practice,
    !    Second Edition,
    !    Addison Wesley, 1995,
    !    ISBN: 0201848406,
    !    LC: T385.C5735.
    !
    !  Parameters:
    !
    !    Input, real(dp) BETA1, the skew or bias parameter.
    !    BETA1 = 1 for no skew or bias.
    !
    !    Input, real(dp) BETA2, the tension parameter.
    !    BETA2 = 0 for no tension.
    !
    !    Output, real(dp) MBASIS(4,4), the basis matrix.
    !
    implicit none

    real(dp) beta1
    real(dp) beta2
    real(dp) delta
    real(dp) mbasis(4,4)

    mbasis(1,1) = - 2.0_dp * beta1 * beta1 * beta1
    mbasis(1,2) =   2.0_dp * beta2 &
         + 2.0 * beta1 * ( beta1 * beta1 + beta1 + 1.0_dp )
    mbasis(1,3) = - 2.0_dp * ( beta2 + beta1 * beta1 + beta1 + 1.0_dp )
    mbasis(1,4) =   2.0_dp

    mbasis(2,1) =   6.0_dp * beta1 * beta1 * beta1
    mbasis(2,2) = - 3.0_dp * beta2 &
         - 6.0_dp * beta1 * beta1 * ( beta1 + 1.0_dp )
    mbasis(2,3) =   3.0_dp * beta2 + 6.0_dp * beta1 * beta1
    mbasis(2,4) =   0.0_dp

    mbasis(3,1) = - 6.0_dp * beta1 * beta1 * beta1
    mbasis(3,2) =   6.0_dp * beta1 * ( beta1 - 1.0_dp ) * ( beta1 + 1.0_dp )
    mbasis(3,3) =   6.0_dp * beta1
    mbasis(3,4) =   0.0_dp

    mbasis(4,1) =   2.0_dp * beta1 * beta1 * beta1
    mbasis(4,2) =   4.0_dp * beta1 * ( beta1 + 1.0_dp ) + beta2
    mbasis(4,3) =   2.0_dp
    mbasis(4,4) =   0.0_dp

    delta = ( ( 2.0_dp   &
         * beta1 + 4.0_dp ) &
         * beta1 + 4.0_dp ) &
         * beta1 + 2.0_dp + beta2

    mbasis(1:4,1:4) = mbasis(1:4,1:4) / delta

    return
  end subroutine basis_matrix_beta_uni
  subroutine basis_matrix_bezier ( mbasis )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_BEZIER sets up the cubic Bezier spline basis matrix.
    !
    !  Discussion:
    !
    !    This basis matrix assumes that the data points are stored as
    !    ( P1, P2, P3, P4 ).  P1 is the function value at T = 0, while
    !    P2 is used to approximate the derivative at T = 0 by
    !    dP/dt = 3 * ( P2 - P1 ).  Similarly, P4 is the function value
    !    at T = 1, and P3 is used to approximate the derivative at T = 1
    !    by dP/dT = 3 * ( P4 - P3 ).
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    05 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    James Foley, Andries vanDam, Steven Feiner, John Hughes,
    !    Computer Graphics, Principles and Practice,
    !    Second Edition,
    !    Addison Wesley, 1995,
    !    ISBN: 0201848406,
    !    LC: T385.C5735.
    !
    !  Parameters:
    !
    !    Output, real(dp) MBASIS(4,4), the basis matrix.
    !
    implicit none

    real(dp) mbasis(4,4)

    mbasis(1,1) = -1.0_dp
    mbasis(1,2) =  3.0_dp
    mbasis(1,3) = -3.0_dp
    mbasis(1,4) =  1.0_dp

    mbasis(2,1) =  3.0_dp
    mbasis(2,2) = -6.0_dp
    mbasis(2,3) =  3.0_dp
    mbasis(2,4) =  0.0_dp

    mbasis(3,1) = -3.0_dp
    mbasis(3,2) =  3.0_dp
    mbasis(3,3) =  0.0_dp
    mbasis(3,4) =  0.0_dp

    mbasis(4,1) =  1.0_dp
    mbasis(4,2) =  0.0_dp
    mbasis(4,3) =  0.0_dp
    mbasis(4,4) =  0.0_dp

    return
  end subroutine basis_matrix_bezier
  subroutine basis_matrix_hermite ( mbasis )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_HERMITE sets up the Hermite spline basis matrix.
    !
    !  Discussion:
    !
    !    This basis matrix assumes that the data points are stored as
    !    ( P1, P2, P1', P2' ), with P1 and P1' being the data value and
    !    the derivative dP/dT at T = 0, while P2 and P2' apply at T = 1.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    06 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    James Foley, Andries vanDam, Steven Feiner, John Hughes,
    !    Computer Graphics, Principles and Practice,
    !    Second Edition,
    !    Addison Wesley, 1995,
    !    ISBN: 0201848406,
    !    LC: T385.C5735.
    !
    !  Parameters:
    !
    !    Output, real(dp) MBASIS(4,4), the basis matrix.
    !
    implicit none

    real(dp) mbasis(4,4)

    mbasis(1,1) =  2.0_dp
    mbasis(1,2) = -2.0_dp
    mbasis(1,3) =  1.0_dp
    mbasis(1,4) =  1.0_dp

    mbasis(2,1) = -3.0_dp
    mbasis(2,2) =  3.0_dp
    mbasis(2,3) = -2.0_dp
    mbasis(2,4) = -1.0_dp

    mbasis(3,1) =  0.0_dp
    mbasis(3,2) =  0.0_dp
    mbasis(3,3) =  1.0_dp
    mbasis(3,4) =  0.0_dp

    mbasis(4,1) =  1.0_dp
    mbasis(4,2) =  0.0_dp
    mbasis(4,3) =  0.0_dp
    mbasis(4,4) =  0.0_dp

    return
  end subroutine basis_matrix_hermite
  subroutine basis_matrix_overhauser_nonuni ( alpha, beta, mbasis )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_OVERHAUSER_NONUNI: nonuniform Overhauser spline basis matrix.
    !
    !  Discussion:
    !
    !    This basis matrix assumes that the data points P1, P2, P3 and
    !    P4 are not uniformly spaced in T, and that P2 corresponds to T = 0,
    !    and P3 to T = 1.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    05 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, real(dp) ALPHA, BETA.
    !    ALPHA = || P2 - P1 || / ( || P3 - P2 || + || P2 - P1 || )
    !    BETA  = || P3 - P2 || / ( || P4 - P3 || + || P3 - P2 || ).
    !
    !    Output, real(dp) MBASIS(4,4), the basis matrix.
    !
    implicit none

    real(dp) alpha
    real(dp) beta
    real(dp) mbasis(4,4)

    mbasis(1,1) = - ( 1.0_dp - alpha ) * ( 1.0_dp - alpha ) / alpha
    mbasis(1,2) =   beta + ( 1.0_dp - alpha ) / alpha
    mbasis(1,3) =   alpha - 1.0_dp / ( 1.0_dp - beta )
    mbasis(1,4) =   beta * beta / ( 1.0_dp - beta )

    mbasis(2,1) =   2.0_dp * ( 1.0_dp - alpha ) * ( 1.0_dp - alpha ) / alpha
    mbasis(2,2) = ( - 2.0_dp * ( 1.0_dp - alpha ) - alpha * beta ) / alpha
    mbasis(2,3) = ( 2.0_dp * ( 1.0_dp - alpha ) &
         - beta * ( 1.0_dp - 2.0_dp * alpha ) ) / ( 1.0_dp - beta )
    mbasis(2,4) = - beta * beta / ( 1.0_dp - beta )

    mbasis(3,1) = - ( 1.0_dp - alpha ) * ( 1.0_dp - alpha ) / alpha
    mbasis(3,2) =   ( 1.0_dp - 2.0_dp * alpha ) / alpha
    mbasis(3,3) =   alpha
    mbasis(3,4) =   0.0_dp

    mbasis(4,1) =   0.0_dp
    mbasis(4,2) =   1.0_dp
    mbasis(4,3) =   0.0_dp
    mbasis(4,4) =   0.0_dp

    return
  end subroutine basis_matrix_overhauser_nonuni
  subroutine basis_matrix_overhauser_nul ( alpha, mbasis )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_OVERHAUSER_NUL sets the nonuniform left Overhauser basis matrix.
    !
    !  Discussion:
    !
    !    This basis matrix assumes that the data points P1, P2, and
    !    P3 are not uniformly spaced in T, and that P1 corresponds to T = 0,
    !    and P2 to T = 1.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    27 August 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, real(dp) ALPHA.
    !    ALPHA = || P2 - P1 || / ( || P3 - P2 || + || P2 - P1 || )
    !
    !    Output, real(dp) MBASIS(3,3), the basis matrix.
    !
    implicit none

    real(dp) alpha
    real(dp) mbasis(3,3)

    mbasis(1,1) =   1.0_dp / alpha
    mbasis(1,2) = - 1.0_dp / ( alpha * ( 1.0_dp - alpha ) )
    mbasis(1,3) =   1.0_dp / ( 1.0_dp - alpha )

    mbasis(2,1) = - ( 1.0_dp + alpha ) / alpha
    mbasis(2,2) =   1.0_dp / ( alpha * ( 1.0_dp - alpha ) )
    mbasis(2,3) = - alpha / ( 1.0_dp - alpha )

    mbasis(3,1) =   1.0_dp
    mbasis(3,2) =   0.0_dp
    mbasis(3,3) =   0.0_dp

    return
  end subroutine basis_matrix_overhauser_nul
  subroutine basis_matrix_overhauser_nur ( beta, mbasis )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_OVERHAUSER_NUR: the nonuniform right Overhauser basis matrix.
    !
    !  Discussion:
    !
    !    This basis matrix assumes that the data points PN-2, PN-1, and
    !    PN are not uniformly spaced in T, and that PN-1 corresponds to T = 0,
    !    and PN to T = 1.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    27 August 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, real(dp) BETA.
    !    BETA = || P(N) - P(N-1) || 
    !         / ( || P(N) - P(N-1) || + || P(N-1) - P(N-2) || )
    !
    !    Output, real(dp) MBASIS(3,3), the basis matrix.
    !
    implicit none

    real(dp) beta
    real(dp) mbasis(3,3)

    mbasis(1,1) =   1.0_dp / beta
    mbasis(1,2) = - 1.0_dp / ( beta * ( 1.0_dp - beta ) )
    mbasis(1,3) =   1.0_dp / ( 1.0_dp - beta )

    mbasis(2,1) = - ( 1.0_dp + beta ) / beta
    mbasis(2,2) =   1.0_dp / ( beta * ( 1.0_dp - beta ) )
    mbasis(2,3) = - beta / ( 1.0_dp - beta )

    mbasis(3,1) =   1.0_dp
    mbasis(3,2) =   0.0_dp
    mbasis(3,3) =   0.0_dp

    return
  end subroutine basis_matrix_overhauser_nur
  subroutine basis_matrix_overhauser_uni ( mbasis )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_OVERHAUSER_UNI sets the uniform Overhauser spline basis matrix.
    !
    !  Discussion:
    !
    !    This basis matrix assumes that the data points P1, P2, P3 and
    !    P4 are uniformly spaced in T, and that P2 corresponds to T = 0,
    !    and P3 to T = 1.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    05 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    James Foley, Andries vanDam, Steven Feiner, John Hughes,
    !    Computer Graphics, Principles and Practice,
    !    Second Edition,
    !    Addison Wesley, 1995,
    !    ISBN: 0201848406,
    !    LC: T385.C5735.
    !
    !  Parameters:
    !
    !    Output, real(dp) MBASIS(4,4), the basis matrix.
    !
    implicit none

    real(dp) mbasis(4,4)

    mbasis(1,1) = - 1.0_dp / 2.0_dp
    mbasis(1,2) =   3.0_dp / 2.0_dp
    mbasis(1,3) = - 3.0_dp / 2.0_dp
    mbasis(1,4) =   1.0_dp / 2.0_dp

    mbasis(2,1) =   2.0_dp / 2.0_dp
    mbasis(2,2) = - 5.0_dp / 2.0_dp
    mbasis(2,3) =   4.0_dp / 2.0_dp
    mbasis(2,4) = - 1.0_dp / 2.0_dp

    mbasis(3,1) = - 1.0_dp / 2.0_dp
    mbasis(3,2) =   0.0_dp
    mbasis(3,3) =   1.0_dp / 2.0_dp
    mbasis(3,4) =   0.0_dp

    mbasis(4,1) =   0.0_dp
    mbasis(4,2) =   2.0_dp / 2.0_dp
    mbasis(4,3) =   0.0_dp
    mbasis(4,4) =   0.0_dp

    return
  end subroutine basis_matrix_overhauser_uni
  subroutine basis_matrix_overhauser_uni_l ( mbasis )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_OVERHAUSER_UNI_L sets the left uniform Overhauser basis matrix.
    !
    !  Discussion:
    !
    !    This basis matrix assumes that the data points P1, P2, and P3
    !    are not uniformly spaced in T, and that P1 corresponds to T = 0,
    !    and P2 to T = 1.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    05 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Output, real(dp) MBASIS(3,3), the basis matrix.
    !
    implicit none

    real(dp) mbasis(3,3)

    mbasis(1,1) =   2.0_dp
    mbasis(1,2) = - 4.0_dp
    mbasis(1,3) =   2.0_dp

    mbasis(2,1) = - 3.0_dp
    mbasis(2,2) =   4.0_dp
    mbasis(2,3) = - 1.0_dp

    mbasis(3,1) =   1.0_dp
    mbasis(3,2) =   0.0_dp
    mbasis(3,3) =   0.0_dp

    return
  end subroutine basis_matrix_overhauser_uni_l
  subroutine basis_matrix_overhauser_uni_r ( mbasis )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_OVERHAUSER_UNI_R sets the right uniform Overhauser basis matrix.
    !
    !  Discussion:
    !
    !    This basis matrix assumes that the data points P(N-2), P(N-1),
    !    and P(N) are uniformly spaced in T, and that P(N-1) corresponds to
    !    T = 0, and P(N) to T = 1.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    05 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Output, real(dp) MBASIS(3,3), the basis matrix.
    !
    implicit none

    real(dp) mbasis(3,3)

    mbasis(1,1) =   2.0_dp
    mbasis(1,2) = - 4.0_dp
    mbasis(1,3) =   2.0_dp

    mbasis(2,1) = - 3.0_dp
    mbasis(2,2) =   4.0_dp
    mbasis(2,3) = - 1.0_dp

    mbasis(3,1) =   1.0_dp
    mbasis(3,2) =   0.0_dp
    mbasis(3,3) =   0.0_dp

    return
  end subroutine basis_matrix_overhauser_uni_r
  subroutine basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval, yval )

    !*****************************************************************************80
    !
    !! BASIS_MATRIX_TMP computes Q = T * MBASIS * P
    !
    !  Discussion:
    !
    !    YDATA is a vector of data values, most frequently the values of some
    !    function sampled at uniformly spaced points.  MBASIS is the basis
    !    matrix for a particular kind of spline.  T is a vector of the
    !    powers of the normalized difference between TVAL and the left
    !    endpoint of the interval.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    10 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) LEFT, indicats that TVAL is in the interval
    !    [ TDATA(LEFT), TDATA(LEFT+1) ], or that this is the "nearest"
    !    interval to TVAL.
    !    For TVAL < TDATA(1), use LEFT = 1.
    !    For TDATA(NDATA) < TVAL, use LEFT = NDATA - 1.
    !
    !    Input, integer(i4) N, the order of the basis matrix.
    !
    !    Input, real(dp) MBASIS(N,N), the basis matrix.
    !
    !    Input, integer(i4) NDATA, the dimension of the vectors TDATA
    !    and YDATA.
    !
    !    Input, real(dp) TDATA(NDATA), the abscissa values.  This routine
    !    assumes that the TDATA values are uniformly spaced, with an
    !    increment of 1.0.
    !
    !    Input, real(dp) YDATA(NDATA), the data values to be
    !    interpolated or approximated.
    !
    !    Input, real(dp) TVAL, the value of T at which the spline is to be
    !    evaluated.
    !
    !    Output, real(dp) YVAL, the value of the spline at TVAL.
    !
    implicit none

    integer(i4), parameter :: maxn = 4
    integer(i4) n
    integer(i4) ndata

    real(dp) arg
    integer(i4) first
    integer(i4) i
    integer(i4) j
    integer(i4) left
    real(dp) mbasis(n,n)
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) tvec(maxn)
    real(dp) ydata(ndata)
    real(dp) yval

    if ( left == 1 ) then
       arg = 0.5_dp * ( tval - tdata(left) )
       first = left
    else if ( left < ndata - 1 ) then
       arg = tval - tdata(left)
       first = left - 1
    else if ( left == ndata - 1 ) then
       arg = 0.5_dp * ( 1.0_dp + tval - tdata(left) )
       first = left - 1
    end if
    !
    !  TVEC(I) = ARG^(N-I).
    !
    tvec(n) = 1.0_dp
    do i = n-1, 1, -1
       tvec(i) = arg * tvec(i+1)
    end do

    yval = 0.0_dp
    do j = 1, n
       yval = yval + dot_product ( tvec(1:n), mbasis(1:n,j) ) &
            * ydata(first - 1 + j)
    end do

    return
  end subroutine basis_matrix_tmp
  subroutine bc_val ( n, t, xcon, ycon, xval, yval )

    !*****************************************************************************80
    !
    !! BC_VAL evaluates a parameterized N-th degree Bezier curve in 2D.
    !
    !  Discussion:
    !
    !    BC_VAL(T) is the value of a vector function of the form
    !
    !      BC_VAL(T) = ( X(T), Y(T) )
    !
    !    where
    !
    !      X(T) = sum ( 0 <= I <= N ) XCON(I) * BERN(I,N)(T)
    !      Y(T) = sum ( 0 <= I <= N ) YCON(I) * BERN(I,N)(T)
    !
    !    BERN(I,N)(T) is the I-th Bernstein polynomial of order N
    !    defined on the interval [0,1],
    !
    !    XCON(0:N) and YCON(0:N) are the coordinates of N+1 "control points".
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    12 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    David Kahaner, Cleve Moler, Steven Nash,
    !    Numerical Methods and Software,
    !    Prentice Hall, 1989,
    !    ISBN: 0-13-627258-4,
    !    LC: TA345.K34.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the degree of the Bezier curve.
    !    N must be at least 0.
    !
    !    Input, real(dp) T, the point at which the Bezier curve should
    !    be evaluated.  The best results are obtained within the interval
    !    [0,1] but T may be anywhere.
    !
    !    Input, real(dp) XCON(0:N), YCON(0:N), the X and Y coordinates
    !    of the control points.  The Bezier curve will pass through
    !    the points ( XCON(0), YCON(0) ) and ( XCON(N), YCON(N) ), but
    !    generally NOT through the other control points.
    !
    !    Output, real(dp) XVAL, YVAL, the X and Y coordinates of the point
    !    on the Bezier curve corresponding to the given T value.
    !
    implicit none

    integer(i4) n

    real(dp) bval(0:n)
    real(dp) t
    real(dp) xcon(0:n)
    real(dp) xval
    real(dp) ycon(0:n)
    real(dp) yval

    call bp01 ( n, t, bval )

    xval = dot_product ( xcon(0:n), bval(0:n) )
    yval = dot_product ( ycon(0:n), bval(0:n) )

    return
  end subroutine bc_val
  function bez_val ( n, x, a, b, y )

    !*****************************************************************************80
    !
    !! BEZ_VAL evaluates an N-th degree Bezier function at a point.
    !
    !  Discussion:
    !
    !    The Bezier function has the form:
    !
    !      BEZ(X) = sum ( 0 <= I <= N ) Y(I) * BERN(N,I)( (X-A)/(B-A) )
    !
    !    BERN(N,I)(X) is the I-th Bernstein polynomial of order N
    !    defined on the interval [0,1],
    !
    !    Y(0:N) is a set of coefficients,
    !
    !    and if, for I = 0 to N, we define the N+1 points
    !
    !      X(I) = ( (N-I)*A + I*B) / N,
    !
    !    equally spaced in [A,B], the pairs ( X(I), Y(I) ) can be regarded as
    !    "control points".
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    12 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    David Kahaner, Cleve Moler, Steven Nash,
    !    Numerical Methods and Software,
    !    Prentice Hall, 1989,
    !    ISBN: 0-13-627258-4,
    !    LC: TA345.K34.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the degree of the Bezier function.
    !    N must be at least 0.
    !
    !    Input, real(dp) X, the point at which the Bezier function should
    !    be evaluated.  The best results are obtained within the interval
    !    [A,B] but X may be anywhere.
    !
    !    Input, real(dp) A, B, the interval over which the Bezier function
    !    has been defined.  This is the interval in which the control
    !    points have been set up.  Note BEZ(A) = Y(0) and BEZ(B) = Y(N),
    !    although BEZ will not, in general pass through the other
    !    control points.  A and B must not be equal.
    !
    !    Input, real(dp) Y(0:N), a set of data defining the Y coordinates
    !    of the control points.
    !
    !    Output, real(dp) BEZ_VAL, the value of the Bezier function at X.
    !
    implicit none

    integer(i4) n

    real(dp) a
    real(dp) b
    real(dp) bez_val
    real(dp) bval(0:n)
    real(dp) x
    real(dp) x01
    real(dp) y(0:n)

    if ( eq(b - a, 0.0_dp) ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'BEZ_VAL - Fatal error!'
       write ( *, '(a,g14.6)' ) '  Null interval, A = B = ', a
       stop
    end if
    !
    !  X01 lies in [0,1], in the same relative position as X in [A,B].
    !
    x01 = ( x - a ) / ( b - a )

    call bp01 ( n, x01, bval )

    bez_val = dot_product ( y(0:n), bval(0:n) )

    return
  end function bez_val
  subroutine bp01 ( n, x, bern )

    !*****************************************************************************80
    !
    !! BP01 evaluates the Bernstein basis polynomials for [0,1] at a point.
    !
    !  Discussion:
    !
    !    For any N greater than or equal to 0, there is a set of N+1 Bernstein
    !    basis polynomials, each of degree N, which form a basis for
    !    all polynomials of degree N on [0,1].
    !
    !    BERN(N,I,X) = [N!/(I!*(N-I)!)] * (1-X)^(N-I) * X^I
    !
    !    N is the degree;
    !
    !    0 <= I <= N indicates which of the N+1 basis polynomials
    !    of degree N to choose;
    !
    !    X is a point in [0,1] at which to evaluate the basis polynomial.
    !
    !  First values:
    !
    !    B(0,0,X) = 1
    !
    !    B(1,0,X) =      1-X
    !    B(1,1,X) =               X
    !
    !    B(2,0,X) =     (1-X)^2
    !    B(2,1,X) = 2 * (1-X)   * X
    !    B(2,2,X) =               X^2
    !
    !    B(3,0,X) =     (1-X)^3
    !    B(3,1,X) = 3 * (1-X)^2 * X
    !    B(3,2,X) = 3 * (1-X)   * X^2
    !    B(3,3,X) =               X^3
    !
    !    B(4,0,X) =     (1-X)^4
    !    B(4,1,X) = 4 * (1-X)^3 * X
    !    B(4,2,X) = 6 * (1-X)^2 * X^2
    !    B(4,3,X) = 4 * (1-X)   * X^3
    !    B(4,4,X) =               X^4
    !
    !  Special values:
    !
    !    B(N,I,1/2) = C(N,K) / 2^N
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    12 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    David Kahaner, Cleve Moler, Steven Nash,
    !    Numerical Methods and Software,
    !    Prentice Hall, 1989,
    !    ISBN: 0-13-627258-4,
    !    LC: TA345.K34.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the degree of the Bernstein basis 
    !    polynomials.  N must be at least 0.
    !
    !    Input, real(dp) X, the evaluation point.
    !
    !    Output, real(dp) BERN(0:N), the values of the N+1 Bernstein basis
    !    polynomials at X.
    !
    implicit none

    integer(i4) n

    real(dp) bern(0:n)
    integer(i4) i
    integer(i4) j
    real(dp) x

    if ( n == 0 ) then

       bern(0) = 1.0_dp

    else if ( 0 < n ) then

       bern(0) = 1.0_dp - x
       bern(1) = x

       do i = 2, n
          bern(i) = x * bern(i-1)
          do j = i-1, 1, -1
             bern(j) = x * bern(j-1) + ( 1.0_dp - x ) * bern(j)
          end do
          bern(0) = ( 1.0_dp - x ) * bern(0)
       end do

    end if

    return
  end subroutine bp01
  subroutine bpab ( n, a, b, x, bern )

    !*****************************************************************************80
    !
    !! BPAB evaluates the Bernstein basis polynomials for [A,B] at a point.
    !
    !  Discussion:
    !
    !    BERN(N,I,X) = [N!/(I!*(N-I)!)] * (B-X)^(N-I) * (X-A)^I / (B-A)^N
    !
    !    B(0,0,X) =   1
    !
    !    B(1,0,X) = (      B-X              ) / (B-A)
    !    B(1,1,X) = (                X-A    ) / (B-A)
    !
    !    B(2,0,X) = (     (B-X)^2           ) / (B-A)^2
    !    B(2,1,X) = ( 2 * (B-X    * (X-A)   ) / (B-A)^2
    !    B(2,2,X) = (               (X-A)^2 ) / (B-A)^2
    !
    !    B(3,0,X) = (     (B-X)^3           ) / (B-A)^3
    !    B(3,1,X) = ( 3 * (B-X)^2 * (X-A)   ) / (B-A)^3
    !    B(3,2,X) = ( 3 * (B-X)   * (X-A)^2 ) / (B-A)^3
    !    B(3,3,X) = (               (X-A)^3 ) / (B-A)^3
    !
    !    B(4,0,X) = (     (B-X)^4           ) / (B-A)^4
    !    B(4,1,X) = ( 4 * (B-X)^3 * (X-A)   ) / (B-A)^4
    !    B(4,2,X) = ( 6 * (B-X)^2 * (X-A)^2 ) / (B-A)^4
    !    B(4,3,X) = ( 4 * (B-X)   * (X-A)^3 ) / (B-A)^4
    !    B(4,4,X) = (               (X-A)^4 ) / (B-A)^4
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    12 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    David Kahaner, Cleve Moler, Steven Nash,
    !    Numerical Methods and Software,
    !    Prentice Hall, 1989,
    !    ISBN: 0-13-627258-4,
    !    LC: TA345.K34.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the degree of the Bernstein basis 
    !    polynomials.  There is a set of N+1 Bernstein basis polynomials, each 
    !    of degree N, which form a basis for polynomials of degree N on [A,B].  
    !    N must be at least 0.
    !
    !    Input, real(dp) A, B, the endpoints of the interval on which the
    !    polynomials are to be based.  A and B should not be equal.
    !    Input, real(dp) X, the point at which the polynomials are to be
    !    evaluated.  X need not lie in the interval [A,B].
    !
    !    Output, real(dp) BERN(0:N), the values of the N+1 Bernstein basis
    !    polynomials at X.
    !
    implicit none

    integer(i4) n

    real(dp) a
    real(dp) b
    real(dp) bern(0:n)
    integer(i4) i
    integer(i4) j
    real(dp) x

    if ( eq(b,a) ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'BPAB - Fatal error!'
       write ( *, '(a,g14.6)' ) '  A = B = ', a
       stop
    end if

    if ( n == 0 ) then

       bern(0) = 1.0_dp

    else if ( 0 < n ) then

       bern(0) = ( b - x ) / ( b - a )
       bern(1) = ( x - a ) / ( b - a )

       do i = 2, n
          bern(i) = ( x - a ) * bern(i-1) / ( b - a )
          do j = i-1, 1, -1
             bern(j) = ( ( b - x ) * bern(j) + ( x - a ) * bern(j-1) ) / ( b - a )
          end do
          bern(0) = ( b - x ) * bern(0) / ( b - a )
       end do

    end if

    return
  end subroutine bpab
  subroutine bpab_approx ( n, a, b, ydata, xval, yval )

    !*****************************************************************************80
    !
    !! BPAB_APPROX evaluates the Bernstein polynomial approximant to F(X) on [A,B].
    !
    !  Formula:
    !
    !    BERN(F)(X) = sum ( 0 <= I <= N ) F(X(I)) * B_BASE(I,X)
    !
    !    where
    !
    !      X(I) = ( ( N - I ) * A + I * B ) / N
    !      B_BASE(I,X) is the value of the I-th Bernstein basis polynomial at X.
    !
    !  Discussion:
    !
    !    The Bernstein polynomial BERN(F) for F(X) is an approximant, not an
    !    interpolant; in other words, its value is not guaranteed to equal
    !    that of F at any particular point.  However, for a fixed interval
    !    [A,B], if we let N increase, the Bernstein polynomial converges
    !    uniformly to F everywhere in [A,B], provided only that F is continuous.
    !    Even if F is not continuous, but is bounded, the polynomial converges
    !    pointwise to F(X) at all points of continuity.  On the other hand,
    !    the convergence is quite slow compared to other interpolation
    !    and approximation schemes.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    12 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    David Kahaner, Cleve Moler, Steven Nash,
    !    Numerical Methods and Software,
    !    Prentice Hall, 1989,
    !    ISBN: 0-13-627258-4,
    !    LC: TA345.K34.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the degree of the Bernstein polynomial
    !    to be used.  N must be at least 0.
    !
    !    Input, real(dp) A, B, the endpoints of the interval on which the
    !    approximant is based.  A and B should not be equal.
    !
    !    Input, real(dp) YDATA(0:N), the data values at N+1 equally
    !    spaced points in [A,B].  If N = 0, then the evaluation point should
    !    be 0.5 * ( A + B).  Otherwise, evaluation point I should be
    !    ( (N-I)*A + I*B ) / N ).
    !
    !    Input, real(dp) XVAL, the point at which the Bernstein polynomial
    !    approximant is to be evaluated.  XVAL does not have to lie in the
    !    interval [A,B].
    !
    !    Output, real(dp) YVAL, the value of the Bernstein polynomial
    !    approximant for F, based in [A,B], evaluated at XVAL.
    !
    implicit none

    integer(i4) n

    real(dp) a
    real(dp) b
    real(dp) bvec(0:n)
    real(dp) xval
    real(dp) ydata(0:n)
    real(dp) yval
    !
    !  Evaluate the Bernstein basis polynomials at XVAL.
    !
    call bpab ( n, a, b, xval, bvec )
    !
    !  Now compute the sum of YDATA(I) * BVEC(I).
    !
    yval = dot_product ( ydata(0:n), bvec(0:n) )

    return
  end subroutine bpab_approx
  subroutine chfev ( x1, x2, f1, f2, d1, d2, ne, xe, fe, next, ierr )

    !*****************************************************************************80
    !
    !! CHFEV evaluates a cubic polynomial given in Hermite form.
    !
    !  Discussion:
    !
    !    This routine evaluates a cubic polynomial given in Hermite form at an
    !    array of points.  While designed for use by SPLINE_PCHIP_VAL, it may
    !    be useful directly as an evaluator for a piecewise cubic
    !    Hermite function in applications, such as graphing, where
    !    the interval is known in advance.
    !
    !    The cubic polynomial is determined by function values
    !    F1, F2 and derivatives D1, D2 on the interval [X1,X2].
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    25 June 2008
    !
    !  Author:
    !
    !    Original FORTRAN77 version by Fred Fritsch.
    !    FORTRAN90 version by John Burkardt.
    !
    !  Reference:
    !
    !    Fred Fritsch, Ralph Carlson,
    !    Monotone Piecewise Cubic Interpolation,
    !    SIAM Journal on Numerical Analysis,
    !    Volume 17, Number 2, April 1980, pages 238-246.
    !
    !    David Kahaner, Cleve Moler, Steven Nash,
    !    Numerical Methods and Software,
    !    Prentice Hall, 1989,
    !    ISBN: 0-13-627258-4,
    !    LC: TA345.K34.
    !
    !  Parameters:
    !
    !    Input, real(dp) X1, X2, the endpoints of the interval of
    !    definition of the cubic.  X1 and X2 must be distinct.
    !
    !    Input, real(dp) F1, F2, the values of the function at X1 and
    !    X2, respectively.
    !
    !    Input, real(dp) D1, D2, the derivative values at X1 and
    !    X2, respectively.
    !
    !    Input, integer(i4) NE, the number of evaluation points.
    !
    !    Input, real(dp) XE(NE), the points at which the function is to
    !    be evaluated.  If any of the XE are outside the interval
    !    [X1,X2], a warning error is returned in NEXT.
    !
    !    Output, real(dp) FE(NE), the value of the cubic function
    !    at the points XE.
    !
    !    Output, integer(i4) NEXT(2), indicates the number of
    !    extrapolation points:
    !    NEXT(1) = number of evaluation points to the left of interval.
    !    NEXT(2) = number of evaluation points to the right of interval.
    !
    !    Output, integer(i4) IERR, error flag.
    !    0, no errors.
    !    -1, NE < 1.
    !    -2, X1 == X2.
    !
    implicit none

    integer(i4) ne

    real(dp) c2
    real(dp) c3
    real(dp) d1
    real(dp) d2
    real(dp) del1
    real(dp) del2
    real(dp) delta
    real(dp) f1
    real(dp) f2
    real(dp) fe(ne)
    real(dp) h
    integer(i4) i
    integer(i4) ierr
    integer(i4) next(2)
    real(dp) x
    real(dp) x1
    real(dp) x2
    real(dp) xe(ne)
    real(dp) xma
    real(dp) xmi

    if ( ne < 1 ) then
       ierr = -1
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'CHFEV - Fatal error!'
       write ( *, '(a)' ) '  Number of evaluation points is less than 1.'
       write ( *, '(a,i8)' ) '  NE = ', ne
       stop
    end if

    h = x2 - x1

    if ( eq(h,0.0_dp) ) then
       ierr = -2
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'CHFEV - Fatal error!'
       write ( *, '(a)' ) '  The interval [X1,X2] is of zero length.'
       stop
    end if
    !
    !  Initialize.
    !
    ierr = 0
    next(1) = 0
    next(2) = 0
    xmi = min ( 0.0_dp, h )
    xma = max ( 0.0_dp, h )
    !
    !  Compute cubic coefficients expanded about X1.
    !
    delta = ( f2 - f1 ) / h
    del1 = ( d1 - delta ) / h
    del2 = ( d2 - delta ) / h
    c2 = -( del1 + del1 + del2 )
    c3 = ( del1 + del2 ) / h
    !
    !  Evaluation loop.
    !
    do i = 1, ne

       x = xe(i) - x1
       fe(i) = f1 + x * ( d1 + x * ( c2 + x * c3 ) )
       !
       !  Count the extrapolation points.
       !
       if ( x < xmi ) then
          next(1) = next(1) + 1
       end if

       if ( xma < x ) then
          next(2) = next(2) + 1
       end if

    end do

    return
  end subroutine chfev
  subroutine data_to_dif ( ntab, xtab, ytab, diftab )

    !*****************************************************************************80
    !
    !! DATA_TO_DIF sets up a divided difference table from raw data.
    !
    !  Discussion:
    !
    !    Space can be saved by using a single array for both the DIFTAB and
    !    YTAB dummy parameters.  In that case, the divided difference table will
    !    overwrite the Y data without interfering with the computation.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    04 September 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Carl deBoor,
    !    A Practical Guide to Splines,
    !    Springer, 2001,
    !    ISBN: 0387953663.
    !
    !  Parameters:
    !
    !    Input, integer(i4) NTAB, the number of pairs of points
    !    (XTAB(I),YTAB(I)) which are to be used as data.  The
    !    number of entries to be used in DIFTAB, XTAB and YTAB.
    !
    !    Input, real(dp) XTAB(NTAB), the X values at which data was taken.
    !    These values must be distinct.
    !
    !    Input, real(dp) YTAB(NTAB), the corresponding Y values.
    !
    !    Output, real(dp) DIFTAB(NTAB), the divided difference coefficients
    !    corresponding to the input (XTAB,YTAB).
    !
    implicit none

    integer(i4) ntab

    real(dp) diftab(ntab)
    integer(i4) i
    integer(i4) j
    !logical r8vec_distinct
    real(dp) xtab(ntab)
    real(dp) ytab(ntab)

    if ( .not. r8vec_distinct ( ntab, xtab ) ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'DATA_TO_DIF - Fatal error!'
       write ( *, '(a)' ) '  Two entries of XTAB are equal!'
       return
    end if
    !
    !  Copy the data values into DIFTAB.
    !
    diftab(1:ntab) = ytab(1:ntab)
    !
    !  Compute the divided differences.
    !
    do i = 2, ntab
       do j = ntab, i, -1

          diftab(j) = ( diftab(j) - diftab(j-1) ) / ( xtab(j) - xtab(j+1-i) )

       end do
    end do

    return
  end subroutine data_to_dif
  subroutine dif_val ( ntab, xtab, diftab, xval, yval )

    !*****************************************************************************80
    !
    !! DIF_VAL evaluates a divided difference polynomial at a point.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    05 September 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Carl deBoor,
    !    A Practical Guide to Splines,
    !    Springer, 2001,
    !    ISBN: 0387953663.
    !
    !  Parameters:
    !
    !    Input, integer(i4) NTAB, the number of divided difference
    !    coefficients in DIFTAB, and the number of points XTAB.
    !
    !    Input, real(dp) XTAB(NTAB), the X values upon which the
    !    divided difference polynomial is based.
    !
    !    Input, real(dp) DIFTAB(NTAB), the divided difference
    !    polynomial coefficients.
    !
    !    Input, real(dp) XVAL, the value where the polynomial
    !    is to be evaluated.
    !
    !    Output, real(dp) YVAL, the value of the polynomial at XVAL.
    !
    implicit none

    integer(i4) ntab

    real(dp) diftab(ntab)
    integer(i4) i
    real(dp) xtab(ntab)
    real(dp) xval
    real(dp) yval

    yval = diftab(ntab)
    do i = 1, ntab-1
       yval = diftab(ntab-i) + ( xval - xtab(ntab-i) ) * yval
    end do

    return
  end subroutine dif_val
  subroutine least_set_old ( ntab, xtab, ytab, ndeg, ptab, b, c, d, eps, ierror )

    !*****************************************************************************80
    !
    !! LEAST_SET_OLD constructs the least squares polynomial approximation to data.
    !
    !  Discussion:
    !
    !    The least squares polynomial is not returned directly as a simple
    !    polynomial.  Instead, it is represented in terms of a set of
    !    orthogonal polynomials appopriate for the given data.  This makes
    !    the computation more accurate, but means that the user can not
    !    easily evaluate the computed polynomial.  Instead, the routine
    !    LEAST_EVAL should be used to evaluate the least squares polynomial
    !    at any point.  (However, the value of the least squares polynomial
    !    at each of the data points is returned as part of this computation.)
    !
    !
    !    A discrete unweighted inner product is used, so that
    !
    !      ( F(X), G(X) ) = sum ( 1 <= I <= NTAB ) F(XTAB(I)) * G(XTAB(I)).
    !
    !    The least squares polynomial is determined using a set of
    !    orthogonal polynomials PHI.  These polynomials can be defined
    !    recursively by:
    !
    !      PHI(0)(X) = 1
    !      PHI(1)(X) = X - B(1)
    !      PHI(I)(X) = ( X - B(I) ) * PHI(I-1)(X) - D(I) * PHI(I-2)(X)
    !
    !    The array B(1:NDEG) contains the values
    !
    !      B(I) = ( X*PHI(I-1), PHI(I-1) ) / ( PHI(I-1), PHI(I-1) )
    !
    !    The array D(2:NDEG) contains the values
    !
    !      D(I) = ( PHI(I-1), PHI(I-1) ) / ( PHI(I-2), PHI(I-2) )
    !
    !    Using this basis, the least squares polynomial can be represented as
    !
    !      P(X)(I) = sum ( 0 <= I <= NDEG ) C(I) * PHI(I)(X)
    !
    !    The array C(0:NDEG) contains the values
    !
    !      C(I) = ( YTAB(I), PHI(I) ) / ( PHI(I), PHI(I) )
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    15 May 2004
    !
    !  Author:
    !
    !    FORTRAN90 version by John Burkardt
    !
    !  Reference:
    !
    !    Gisela Engeln-Muellges, Frank Uhlig,
    !    Numerical Algorithms with C,
    !    Springer, 1996,
    !    ISBN: 3-540-60530-4.
    !
    !  Parameters:
    !
    !    Input, integer(i4) NTAB, the number of data points.
    !
    !    Input, real(dp) XTAB(NTAB), the X data.  The values in XTAB
    !    should be distinct, and in increasing order.
    !
    !    Input, real(dp) YTAB(NTAB), the Y data values corresponding
    !    to the X data in XTAB.
    !
    !    Input, integer(i4) NDEG, the degree of the polynomial which the
    !    program is to use.  NDEG must be at least 0, and less than or
    !    equal to NTAB-1.
    !
    !    Output, real(dp) PTAB(NTAB), the value of the least
    !    squares polynomial at the points XTAB(1:NTAB).
    !
    !    Output, real(dp) B(1:NDEG), C(0:NDEG), D(2:NDEG), arrays
    !    needed to evaluate the polynomial.
    !
    !    Output, real(dp) EPS, the root-mean-square discrepancy of the
    !    polynomial fit.
    !
    !    Output, integer(i4) IERROR, error flag.
    !    zero, no error occurred;
    !    nonzero, an error occurred, and the polynomial could not be computed.
    !
    implicit none

    integer(i4) ndeg
    integer(i4) ntab

    real(dp) b(1:ndeg)
    real(dp) c(0:ndeg)
    real(dp) d(2:ndeg)
    real(dp) eps
    integer(i4) i
    integer(i4) i0l1
    integer(i4) i1l1
    integer(i4) ierror
    integer(i4) it
    integer(i4) k
    integer(i4) mdeg
    real(dp) ptab(ntab)
    real(dp) rn0
    real(dp) rn1
    real(dp) s
    real(dp) sum2
    real(dp) xtab(ntab)
    real(dp) y_sum
    real(dp) ytab(ntab)
    real(dp) ztab(2*ntab)

    ierror = 0
    !
    !  Check NDEG.
    !
    if ( ndeg < 0 ) then
       ierror = 1
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'LEAST_SET_OLD - Fatal error!'
       write ( *, '(a)' ) '  NDEG < 0.'
       stop
    end if

    if ( ntab <= ndeg ) then
       ierror = 1
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'LEAST_SET_OLD - Fatal error!'
       write ( *, '(a)' ) '  NTAB <= NDEG.'
       stop
    end if
    !
    !  Check that the abscissas are strictly increasing.
    !
    do i = 1, ntab-1
       if ( le(xtab(i+1),xtab(i)) ) then
          ierror = 1
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'LEAST_SET_OLD - Fatal error!'
          write ( *, '(a)' ) '  XTAB must be strictly increasing, but'
          write ( *, '(a,i8,a,g14.6)' ) '  XTAB(', i, ') = ', xtab(i)
          write ( *, '(a,i8,a,g14.6)' ) '  XTAB(', i+1, ') = ', xtab(i+1)
          stop
       end if
    end do

    i0l1 = 0
    i1l1 = ntab
    !
    !  The polynomial is of degree at least 0.
    !
    y_sum = sum ( ytab(1:ntab) )
    rn0 = ntab
    c(0) = y_sum / real ( ntab, kind=dp )

    ptab(1:ntab) = y_sum / real ( ntab, kind=dp )

    if ( ndeg == 0 ) then
       eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) )**2 )
       eps = sqrt ( eps / real ( ntab, kind=dp ) )
       return
    end if
    !
    !  The polynomial is of degree at least 1.
    !
    b(1) = sum ( xtab(1:ntab) ) / real ( ntab, kind=dp )

    s = 0.0_dp
    sum2 = 0.0_dp
    do i = 1, ntab
       ztab(i1l1+i) = xtab(i) - b(1)
       s = s + ztab(i1l1+i)**2
       sum2 = sum2 + ztab(i1l1+i) * ( ytab(i) - ptab(i) )
    end do

    rn1 = s
    c(1) = sum2 / s

    do i = 1, ntab
       ptab(i) = ptab(i) + c(1) * ztab(i1l1+i)
    end do

    if ( ndeg == 1 ) then
       eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) )**2 )
       eps = sqrt ( eps / real ( ntab, kind=dp ) )
       return
    end if

    ztab(1:ntab) = 1.0_dp

    mdeg = 2
    k = 2

    do

       d(k) = rn1 / rn0

       sum2 = 0.0_dp
       do i = 1, ntab
          sum2 = sum2 + xtab(i) * ztab(i1l1+i) * ztab(i1l1+i)
       end do

       b(k) = sum2 / rn1

       s = 0.0_dp
       sum2 = 0.0_dp
       do i = 1, ntab
          ztab(i0l1+i) = ( xtab(i) - b(k) ) * ztab(i1l1+i) &
               - d(k) * ztab(i0l1+i)
          s = s + ztab(i0l1+i) * ztab(i0l1+i)
          sum2 = sum2 + ztab(i0l1+i) * ( ytab(i) - ptab(i) )
       end do

       rn0 = rn1
       rn1 = s

       c(k) = sum2 / rn1

       it = i0l1
       i0l1 = i1l1
       i1l1 = it

       do i = 1, ntab
          ptab(i) = ptab(i) + c(k) * ztab(i1l1+i)
       end do

       if ( ndeg <= mdeg ) then
          exit
       end if

       mdeg = mdeg + 1
       k = k + 1

    end do
    !
    !  Compute the RMS error.
    !
    eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) )**2 )
    eps = sqrt ( eps / real ( ntab, kind=dp ) )

    return
  end subroutine least_set_old
  subroutine least_val_old ( x, ndeg, b, c, d, value )

    !*****************************************************************************80
    !
    !! LEAST_VAL_OLD evaluates a least squares polynomial defined by LEAST_SET_OLD.
    !
    !  Discussion:
    !
    !    This is an "old" version of the routine.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    15 May 2004
    !
    !  Author:
    !
    !    FORTRAN90 version by John Burkardt
    !
    !  Reference:
    !
    !    Gisela Engeln-Muellges, Frank Uhlig,
    !    Numerical Algorithms with C,
    !    Springer, 1996,
    !    ISBN: 3-540-60530-4.
    !
    !  Parameters:
    !
    !    Input, real(dp) X, the point at which the polynomial is
    !    to be evaluated.
    !
    !    Input, integer(i4) NDEG, the degree of the least squares
    !    polynomial.
    !
    !    Input, real(dp) B(1:NDEG), C(0:NDEG), D(2:NDEG), arrays
    !    defined by LEAST_SET_OLD, and needed to evaluate the polynomial.
    !
    !    Output, real(dp) VALUE, the value of the polynomial at X.
    !
    implicit none

    integer(i4) ndeg

    real(dp) b(1:ndeg)
    real(dp) c(0:ndeg)
    real(dp) d(2:ndeg)
    integer(i4) k
    real(dp) sk
    real(dp) skp1
    real(dp) skp2
    real(dp) value
    real(dp) x

    if ( ndeg <= 0 ) then

       value = c(0)

    else if ( ndeg == 1 ) then

       value = c(0) + c(1) * ( x - b(1) )

    else

       skp2 = c(ndeg)
       skp1 = c(ndeg-1) + ( x - b(ndeg) ) * skp2

       do k = ndeg-2, 0, -1
          sk = c(k) + ( x - b(k+1) ) * skp1 - d(k+2) * skp2
          skp2 = skp1
          skp1 = sk
       end do

       value = sk

    end if

    return
  end subroutine least_val_old
  subroutine least_set ( point_num, x, f, w, nterms, b, c, d )

    !*****************************************************************************80
    !
    !! LEAST_SET defines a least squares polynomial for given data.
    !
    !  Discussion:
    !
    !    This routine is based on ORTPOL by Conte and deBoor.
    !
    !    The polynomial may be evaluated at any point X by calling LEAST_VAL.
    !
    !    Thanks to Andrew Telford for pointing out a mistake in the form of
    !    the check that there are enough unique data points, 25 June 2008.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    25 June 2008
    !
    !  Author:
    !
    !    Original FORTRAN77 version by Samuel Conte, Carl deBoor.
    !    FORTRAN90 version by John Burkardt
    !
    !  Reference:
    !
    !    Samuel Conte, Carl deBoor,
    !    Elementary Numerical Analysis,
    !    Second Edition,
    !    McGraw Hill, 1972,
    !    ISBN: 07-012446-4,
    !    LC: QA297.C65.
    !
    !  Parameters:
    !
    !    Input, integer(i4) POINT_NUM, the number of data values.
    !
    !    Input, real(dp) X(POINT_NUM), the abscissas of the data points.
    !    At least NTERMS of the values in X must be distinct.
    !
    !    Input, real(dp) F(POINT_NUM), the data values at the points X(*).
    !
    !    Input, real(dp) W(POINT_NUM), the weights associated with
    !    the data points.  Each entry of W should be positive.
    !
    !    Input, integer(i4) NTERMS, the number of terms to use in the
    !    approximating polynomial.  NTERMS must be at least 1.
    !    The degree of the polynomial is NTERMS-1.
    !
    !    Output, real(dp) B(NTERMS), C(NTERMS), D(NTERMS), are quantities
    !    defining the least squares polynomial for the input data,
    !    which will be needed to evaluate the polynomial.
    !
    implicit none

    integer(i4) point_num
    integer(i4) nterms

    real(dp) b(nterms)
    real(dp) c(nterms)
    real(dp) d(nterms)
    real(dp) f(point_num)
    integer(i4) i
    integer(i4) j
    real(dp) p
    real(dp) pj(point_num)
    real(dp) pjm1(point_num)
    real(dp) s(nterms)
    real(dp), parameter :: tol = 0.0_dp
    integer(i4) unique_num
    real(dp) w(point_num)
    real(dp) x(point_num)
    !
    !  Make sure at least NTERMS X values are unique.
    !
    call r8vec_unique_count ( point_num, x, tol, unique_num )

    if ( unique_num < nterms ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'LEAST_SET - Fatal error!'
       write ( *, '(a)' ) '  The number of distinct X values must be'
       write ( *, '(a,i8)') '  at least NTERMS = ', nterms
       write ( *, '(a,i8)' ) '  but the input data has only ', unique_num
       write ( *, '(a)' ) '  distinct entries.'
       return
    end if
    !
    !  Make sure all W entries are positive.
    !
    do i = 1, point_num
       if ( le(w(i),0.0_dp) ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'LEAST_SET - Fatal error!'
          write ( *, '(a)' ) '  All weights W must be positive,'
          write ( *, '(a,i8)' ) '  but weight ', i
          write ( *, '(a,g14.6)' ) '  is ', w(i)
          return
       end if
    end do
    !
    !  Start inner product summations at zero.
    !
    b(1:nterms) = 0.0_dp
    c(1:nterms) = 0.0_dp
    d(1:nterms) = 0.0_dp
    s(1:nterms) = 0.0_dp
    !
    !  Set the values of P(-1,X) and P(0,X) at all data points.
    !
    pjm1(1:point_num) = 0.0_dp
    pj(1:point_num) = 1.0_dp
    !
    !  Now compute the value of P(J,X(I)) as
    !
    !    P(J,X(I)) = ( X(I) - B(J) ) * P(J-1,X(I)) - C(J) * P(J-2,X(I))
    !
    !  where
    !
    !    S(J) = < P(J,X), P(J,X) >
    !    B(J) = < x*P(J,X), P(J,X) > / < P(J,X), P(J,X) >
    !    C(J) = S(J) / S(J-1)
    !
    !  and the least squares coefficients are
    !
    !    D(J) = < F(X), P(J,X) > / < P(J,X), P(J,X) >
    !
    do j = 1, nterms

       d(j) = d(j) + sum ( w(1:point_num) * f(1:point_num) * pj(1:point_num) )
       b(j) = b(j) + sum ( w(1:point_num) * x(1:point_num) * pj(1:point_num)**2 )
       s(j) = s(j) + sum ( w(1:point_num) * pj(1:point_num)**2 )

       d(j) = d(j) / s(j)

       if ( j == nterms ) then
          c(j) = 0.0_dp
          return
       end if

       b(j) = b(j) / s(j)

       if ( j == 1 ) then
          c(j) = 0.0_dp
       else
          c(j) = s(j) / s(j-1)
       end if

       do i = 1, point_num
          p = pj(i)
          pj(i) = ( x(i) - b(j) ) * pj(i) - c(j) * pjm1(i)
          pjm1(i) = p
       end do

    end do

    return
  end subroutine least_set
  subroutine least_val ( nterms, b, c, d, x, px )

    !*****************************************************************************80
    !
    !! LEAST_VAL evaluates a least squares polynomial defined by LEAST_SET.
    !
    !  Discussion:
    !
    !    The least squares polynomial is assumed to be defined as a sum
    !
    !      P(X) = sum ( 1 <= I <= NTERMS ) D(I) * P(I-1,X)
    !
    !    where the orthogonal basis polynomials P(I,X) satisfy the following
    !    three term recurrence:
    !
    !      P(-1,X) = 0
    !      P(0,X) = 1
    !      P(I,X) = ( X - B(I-1) ) * P(I-1,X) - C(I-1) * P(I-2,X)
    !
    !    Therefore, the least squares polynomial can be evaluated as follows:
    !
    !    If NTERMS is 1, then the value of P(X) is D(1) * P(0,X) = D(1).
    !
    !    Otherwise, P(X) is defined as the sum of NTERMS > 1 terms.  We can
    !    reduce the number of terms by 1, because the polynomial P(NTERMS,X)
    !    can be rewritten as a sum of polynomials;  Therefore, P(NTERMS,X)
    !    can be eliminated from the sum, and its coefficient merged in with
    !    those of other polynomials.  Repeat this process for P(NTERMS-1,X)
    !    and so on until a single term remains.
    !    P(NTERMS,X) of P(NTERMS-1,X)
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    24 May 2005
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Samuel Conte, Carl deBoor,
    !    Elementary Numerical Analysis,
    !    Second Edition,
    !    McGraw Hill, 1972,
    !    ISBN: 07-012446-4,
    !    LC: QA297.C65.
    !
    !  Parameters:
    !
    !    Input, integer(i4) NTERMS, the number of terms in the least
    !    squares polynomial.  NTERMS must be at least 1.  The input value of NTERMS
    !    may be reduced from the value given to LEAST_SET.  This will
    !    evaluate the least squares polynomial of the lower degree specified.
    !
    !    Input, real(dp) B(NTERMS), C(NTERMS), D(NTERMS), the information
    !    computed by LEAST_SET.
    !
    !    Input, real(dp) X, the point at which the least squares polynomial
    !    is to be evaluated.
    !
    !    Output, real(dp) PX, the value of the least squares
    !    polynomial at X.
    !
    implicit none

    integer(i4) nterms

    real(dp) b(nterms)
    real(dp) c(nterms)
    real(dp) d(nterms)
    integer(i4) i
    real(dp) prev
    real(dp) prev2
    real(dp) px
    real(dp) x

    px = d(nterms)
    prev = 0.0_dp

    do i = nterms - 1, 1, -1

       prev2 = prev
       prev = px

       if ( i == nterms-1 ) then
          px = d(i) + ( x - b(i) ) * prev
       else
          px = d(i) + ( x - b(i) ) * prev - c(i+1) * prev2
       end if

    end do

    return
  end subroutine least_val
  subroutine least_val2 ( nterms, b, c, d, x, px, pxp )

    !*****************************************************************************80
    !
    !! LEAST_VAL2 evaluates a least squares polynomial defined by LEAST_SET.
    !
    !  Discussion:
    !
    !    This routine also computes the derivative of the polynomial.
    !
    !    The least squares polynomial is assumed to be defined as a sum
    !
    !      P(X) = sum ( 1 <= I <= NTERMS ) D(I) * P(I-1,X)
    !
    !    where the orthogonal basis polynomials P(I,X) satisfy the following
    !    three term recurrence:
    !
    !      P(-1,X) = 0
    !      P(0,X) = 1
    !      P(I,X) = ( X - B(I-1) ) * P(I-1,X) - C(I-1) * P(I-2,X)
    !
    !    Therefore, the least squares polynomial can be evaluated as follows:
    !
    !    If NTERMS is 1, then the value of P(X) is D(1) * P(0,X) = D(1).
    !
    !    Otherwise, P(X) is defined as the sum of NTERMS > 1 terms.  We can
    !    reduce the number of terms by 1, because the polynomial P(NTERMS,X)
    !    can be rewritten as a sum of polynomials;  Therefore, P(NTERMS,X)
    !    can be eliminated from the sum, and its coefficient merged in with
    !    those of other polynomials.  Repeat this process for P(NTERMS-1,X)
    !    and so on until a single term remains.
    !    P(NTERMS,X) of P(NTERMS-1,X)
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    24 May 2005
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) NTERMS, the number of terms in the least 
    !    squares polynomial.  NTERMS must be at least 1.  The value of NTERMS
    !    may be reduced from the value given to LEAST_SET.
    !    This will cause LEAST_VAL to evaluate the least squares polynomial
    !    of the lower degree specified.
    !
    !    Input, real(dp) B(NTERMS), C(NTERMS), D(NTERMS), the information
    !    computed by LEAST_SET.
    !
    !    Input, real(dp) X, the point at which the least squares polynomial
    !    is to be evaluated.
    !
    !    Output, real(dp) PX, PXP, the value and derivative of the least
    !    squares polynomial at X.
    !
    implicit none

    integer(i4) nterms

    real(dp) b(nterms)
    real(dp) c(nterms)
    real(dp) d(nterms)
    integer(i4) i
    real(dp) px
    real(dp) pxm1
    real(dp) pxm2
    real(dp) pxp
    real(dp) pxpm1
    real(dp) pxpm2
    real(dp) x

    px = d(nterms)
    pxp = 0.0_dp
    pxm1 = 0.0_dp
    pxpm1 = 0.0_dp

    do i = nterms - 1, 1, -1

       pxm2 = pxm1
       pxpm2 = pxpm1
       pxm1 = px
       pxpm1 = pxp

       if ( i == nterms - 1 ) then
          px = d(i) + ( x - b(i) ) * pxm1
          pxp = pxm1 + ( x - b(i) ) * pxpm1
       else
          px = d(i) + ( x - b(i) ) * pxm1 - c(i+1) * pxm2
          pxp = pxm1 + ( x - b(i) ) * pxpm1 - c(i+1) * pxpm2
       end if

    end do

    return
  end subroutine least_val2
  subroutine parabola_val2 ( dim_num, ndata, tdata, ydata, left, tval, yval )

    !*****************************************************************************80
    !
    !! PARABOLA_VAL2 evaluates a parabolic interpolant through tabular data.
    !
    !  Discussion:
    !
    !    This routine is a utility routine used by OVERHAUSER_SPLINE_VAL.
    !    It constructs the parabolic interpolant through the data in
    !    3 consecutive entries of a table and evaluates this interpolant
    !    at a given abscissa value.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    26 January 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) DIM_NUM, the dimension of a single data point.
    !    DIM_NUM must be at least 1.
    !
    !    Input, integer(i4) NDATA, the number of data points.
    !    NDATA must be at least 3.
    !
    !    Input, real(dp) TDATA(NDATA), the abscissas of the data
    !    points.  The values in TDATA must be in strictly ascending order.
    !
    !    Input, real(dp) YDATA(DIM_NUM,NDATA), the data points
    !    corresponding to the abscissas.
    !
    !    Input, integer(i4) LEFT, the location of the first of the three
    !    consecutive data points through which the parabolic interpolant
    !    must pass.  1 <= LEFT <= NDATA - 2.
    !
    !    Input, real(dp) TVAL, the value of T at which the parabolic
    !    interpolant is to be evaluated.  Normally, TDATA(1) <= TVAL <= T(NDATA),
    !    and the data will be interpolated.  For TVAL outside this range,
    !    extrapolation will be used.
    !
    !    Output, real(dp) YVAL(DIM_NUM), the value of the parabolic
    !    interpolant at TVAL.
    !
    implicit none

    integer(i4) ndata
    integer(i4) dim_num

    real(dp) dif1
    real(dp) dif2
    integer(i4) i
    integer(i4) left
    real(dp) t1
    real(dp) t2
    real(dp) t3
    real(dp) tval
    real(dp) tdata(ndata)
    real(dp) ydata(dim_num,ndata)
    real(dp) y1
    real(dp) y2
    real(dp) y3
    real(dp) yval(dim_num)
    !
    !  Check.
    !
    if ( left < 1 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!'
       write ( *, '(a)' ) '  LEFT < 1.'
       write ( *, '(a,i8)' ) '  LEFT = ', left
       stop
    end if

    if ( ndata-2 < left ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!'
       write ( *, '(a)' ) '  NDATA-2 < LEFT.'
       write ( *, '(a,i8)' ) '  NDATA = ', ndata
       write ( *, '(a,i8)' ) '  LEFT =  ', left
       stop
    end if

    if ( dim_num < 1 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!'
       write ( *, '(a)' ) '  DIM_NUM < 1.'
       stop
    end if
    !
    !  Copy out the three abscissas.
    !
    t1 = tdata(left)
    t2 = tdata(left+1)
    t3 = tdata(left+2)

    if ( le(t2,t1) .or. le(t3,t2) ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'PARABOLA_VAL2 - Fatal error!'
       write ( *, '(a)' ) '  T2 <= T1 or T3 <= T2.'
       stop
    end if
    !
    !  Construct and evaluate a parabolic interpolant for the data
    !  in each dimension.
    !
    do i = 1, dim_num

       y1 = ydata(i,left)
       y2 = ydata(i,left+1)
       y3 = ydata(i,left+2)

       dif1 = ( y2 - y1 ) / ( t2 - t1 )
       dif2 = ( ( y3 - y1 ) / ( t3 - t1 ) &
            - ( y2 - y1 ) / ( t2 - t1 ) ) / ( t3 - t2 )

       yval(i) = y1 + ( tval - t1 ) * ( dif1 + ( tval - t2 ) * dif2 )

    end do

    return
  end subroutine parabola_val2
  function pchst ( arg1, arg2 )

    !*****************************************************************************80
    !
    !! PCHST: PCHIP sign-testing routine.
    !
    !  Discussion:
    !
    !    This routine essentially computes the sign of ARG1 * ARG2.
    !
    !    The object is to do this without multiplying ARG1 * ARG2, to avoid
    !    possible over/underflow problems.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    25 June 2008
    !
    !  Author:
    !
    !    Original FORTRAN77 version by Fred Fritsch.
    !    FORTRAN90 version by John Burkardt.
    !
    !  Reference:
    !
    !    Fred Fritsch, Ralph Carlson,
    !    Monotone Piecewise Cubic Interpolation,
    !    SIAM Journal on Numerical Analysis,
    !    Volume 17, Number 2, April 1980, pages 238-246.
    !
    !  Parameters:
    !
    !    Input, real(dp) ARG1, ARG2, two values to check.
    !
    !    Output, real(dp) PCHST,
    !    -1.0, if ARG1 and ARG2 are of opposite sign.
    !     0.0, if either argument is zero.
    !    +1.0, if ARG1 and ARG2 are of the same sign.
    !
    implicit none

    real(dp) arg1
    real(dp) arg2
    real(dp) pchst

    pchst = sign ( 1.0_dp, arg1 ) * sign ( 1.0_dp, arg2 )

    if ( eq(arg1,0.0_dp) .or. eq(arg2,0.0_dp) ) then
       pchst = 0.0_dp
    end if

    return
  end function pchst
  subroutine r8_swap ( x, y )

    !*****************************************************************************80
    !
    !! R8_SWAP swaps two real values.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    01 May 2000
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input/output, real(dp) X, Y.  On output, the values of X and
    !    Y have been interchanged.
    !
    implicit none

    real(dp) x
    real(dp) y
    real(dp) z

    z = x
    x = y
    y = z

    return
  end subroutine r8_swap
  function r8_uniform_01 ( seed )

    !*****************************************************************************80
    !
    !! R8_UNIFORM_01 is a portable pseudorandom number generator.
    !
    !  Discussion:
    !
    !    This routine implements the recursion
    !
    !      seed = 16807 * seed mod ( 2**31 - 1 )
    !      unif = seed / ( 2**31 - 1 )
    !
    !    The integer arithmetic never requires more than 32 bits,
    !    including a sign bit.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    11 August 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Paul Bratley, Bennett Fox, Linus Schrage,
    !    A Guide to Simulation,
    !    Springer Verlag, pages 201-202, 1983.
    !
    !    Bennett Fox,
    !    Algorithm 647:
    !    Implementation and Relative Efficiency of Quasirandom
    !    Sequence Generators,
    !    ACM Transactions on Mathematical Software,
    !    Volume 12, Number 4, pages 362-376, 1986.
    !
    !  Parameters:
    !
    !    Input/output, integer(i4) SEED, the "seed" value, which should
    !    NOT be 0. (Otherwise, the output values of SEED and UNIFORM will be zero.)
    !    On output, SEED has been updated.
    !
    !    Output, real(dp) R8_UNIFORM_01, a new pseudorandom variate,
    !    strictly between 0 and 1.
    !
    implicit none

    integer(i4) k
    integer(i4) seed
    real(dp) r8_uniform_01

    k = seed / 127773

    seed = 16807 * ( seed - k * 127773 ) - k * 2836

    if ( seed < 0 ) then
       seed = seed + 2147483647
    end if
    !
    !  Although SEED can be represented exactly as a 32 bit integer,
    !  it generally cannot be represented exactly as a 32 bit real number!
    !
    r8_uniform_01 = real ( seed, kind=dp ) * 4.656612875E-10

    return
  end function r8_uniform_01
  subroutine r83_mxv ( n, a, x, b )

    !*****************************************************************************80
    !
    !! R83_MXV multiplies an R83 matrix times a vector.
    !
    !  Discussion:
    !
    !    The R83 storage format is used for a tridiagonal matrix.
    !    The superdiagonal is stored in entries (1,2:N), the diagonal in
    !    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
    !    original matrix is "collapsed" vertically into the array.
    !
    !  Example:
    !
    !    Here is how an R83 matrix of order 5 would be stored:
    !
    !       *  A12 A23 A34 A45
    !      A11 A22 A33 A44 A55
    !      A21 A32 A43 A54  *
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    02 November 2003
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the order of the linear system.
    !
    !    Input, real(dp) A(3,N), the R83 matrix.
    !
    !    Input, real(dp) X(N), the vector to be multiplied by A.
    !
    !    Output, real(dp) B(N), the product A * x.
    !
    implicit none

    integer(i4) n

    real(dp) a(3,n)
    real(dp) b(n)
    real(dp) x(n)

    b(1:n)   =            a(2,1:n)   * x(1:n)
    b(1:n-1) = b(1:n-1) + a(1,2:n)   * x(2:n)
    b(2:n)   = b(2:n)   + a(3,1:n-1) * x(1:n-1)

    return
  end subroutine r83_mxv
  subroutine r83_np_fs ( n, a, b, x )

    !*****************************************************************************80
    !
    !! R83_NP_FS factors and solves an R83 system.
    !
    !  Discussion:
    !
    !    The R83 storage format is used for a tridiagonal matrix.
    !    The superdiagonal is stored in entries (1,2:N), the diagonal in
    !    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
    !    original matrix is "collapsed" vertically into the array.
    !
    !    This algorithm requires that each diagonal entry be nonzero.
    !    It does not use pivoting, and so can fail on systems that
    !    are actually nonsingular.
    !
    !    Here is how an R83 matrix of order 5 would be stored:
    !
    !       *  A12 A23 A34 A45
    !      A11 A22 A33 A44 A55
    !      A21 A32 A43 A54  *
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    02 November 2003
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the order of the linear system.
    !
    !    Input/output, real(dp) A(3,N).
    !    On input, the tridiagonal matrix.
    !    On output, the data in these vectors has been overwritten
    !    by factorization information.
    !
    !    Input, real(dp) B(N), the right hand side of the linear system.
    !
    !    Output, real(dp) X(N), the solution of the linear system.
    !
    implicit none

    integer(i4) n

    real(dp) a(3,n)
    real(dp) b(n)
    integer(i4) i
    real(dp) x(n)
    real(dp) xmult
    !
    !  The diagonal entries can't be zero.
    !
    do i = 1, n
       if ( eq(a(2,i),0.0_dp) ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'R83_NP_FS - Fatal error!'
          write ( *, '(a,i8,a)' ) '  A(2,', i, ') = 0.'
          return
       end if
    end do

    x(1:n) = b(1:n)

    do i = 2, n
       xmult = a(3,i-1) / a(2,i-1)
       a(2,i) = a(2,i) - xmult * a(1,i)
       x(i)   = x(i)   - xmult * x(i-1)
    end do

    x(n) = x(n) / a(2,n)
    do i = n - 1, 1, -1
       x(i) = ( x(i) - a(1,i+1) * x(i+1) ) / a(2,i)
    end do

    return
  end subroutine r83_np_fs
  subroutine r83_uniform ( n, seed, a )

    !*****************************************************************************80
    !
    !! R83_UNIFORM randomizes an R83 matrix.
    !
    !  Discussion:
    !
    !    The R83 storage format is used for a tridiagonal matrix.
    !    The superdiagonal is stored in entries (1,2:N), the diagonal in
    !    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
    !    original matrix is "collapsed" vertically into the array.
    !
    !    Here is how an R83 matrix of order 5 would be stored:
    !
    !       *  A12 A23 A34 A45
    !      A11 A22 A33 A44 A55
    !      A21 A32 A43 A54  *
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    19 September 2003
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the order of the linear system.
    !
    !    Input/output, integer(i4) SEED, a seed for the random number
    !    generator.
    !
    !    Output, real(dp) A(3,N), the R83 matrix.
    !
    implicit none

    integer(i4) n

    real(dp) a(3,n)
    integer(i4) seed

    a(1,1) = 0.0_dp
    call r8vec_uniform_01 ( n-1, seed, a(1,2:n) )

    call r8vec_uniform_01 ( n,   seed, a(2,1:n) )

    call r8vec_uniform_01 ( n-1, seed, a(3,1:n-1) )
    a(3,n) = 0.0_dp

    return
  end subroutine r83_uniform
  subroutine r85_np_fs ( n, a, b, x )

    !*****************************************************************************80
    !
    !! R85_NP_FS factors and solves an R85 linear system.
    !
    !  Discussion:
    !
    !    The R85 storage format represents a pentadiagonal matrix as a 5 
    !    by N array, in which each row corresponds to a diagonal, and 
    !    column locations are preserved.  Thus, the original matrix is 
    !    "collapsed" vertically into the array.
    !
    !    The factorization algorithm requires that each diagonal entry be nonzero.
    !
    !    No pivoting is performed, and therefore the algorithm may fail
    !    in simple cases where the matrix is not singular.
    !
    !  Example:
    !
    !    Here is how an R85 matrix of order 6 would be stored:
    !
    !       *   *  A13 A24 A35 A46
    !       *  A12 A23 A34 A45 A56
    !      A11 A22 A33 A44 A55 A66
    !      A21 A32 A43 A54 A65  *
    !      A31 A42 A53 A64  *   *
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license. 
    !
    !  Modified:
    !
    !    17 September 2003
    !
    !  Author:
    !
    !    Original FORTRAN77 version by Ward Cheney, David Kincaid.
    !    FORTRAN90 version by John Burkardt.
    !
    !  Reference:
    !
    !    Ward Cheney, David Kincaid,
    !    Numerical Mathematics and Computing,
    !    Brooks-Cole Publishing, 2004,
    !    ISBN: 0534201121.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the order of the linear system.
    !
    !    Input/output, real(dp) A(5,N),
    !    On input, the pentadiagonal matrix.
    !    On output, the data has been overwritten by factorization information.
    !
    !    Input/output, real(dp) B(N).
    !    On input, B contains the right hand side of the linear system.
    !    On output, B has been overwritten by factorization information.
    !
    !    Output, real(dp) X(N), the solution of the linear system.
    !
    implicit none

    integer(i4) n

    real(dp) a(5,n)
    real(dp) b(n)
    integer(i4) i
    real(dp) x(n)
    real(dp) xmult

    do i = 1, n
       if ( eq(a(3,i),0.0_dp) ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'R85_NP_FS - Fatal error!'
          write ( *, '(a,i8,a)' ) '  A(3,', i, ') = 0.'
          stop
       end if
    end do

    do i = 2, n - 1

       xmult = a(2,i) / a(3,i-1)
       a(3,i) = a(3,i) - xmult * a(4,i-1)
       a(4,i) = a(4,i) - xmult * a(5,i-1)

       b(i) = b(i) - xmult * b(i-1)

       xmult = a(1,i+1) / a(3,i-1)
       a(2,i+1) = a(2,i+1) - xmult * a(4,i-1)
       a(3,i+1) = a(3,i+1) - xmult * a(5,i-1)

       b(i+1) = b(i+1) - xmult * b(i-1)

    end do

    call r8vec_print ( n, b, '  Bpart1: ')

    xmult = a(2,n) / a(3,n-1)
    a(3,n) = a(3,n) - xmult * a(4,n-1)

    x(n) = ( b(n) - xmult * b(n-1) ) / a(3,n)
    x(n-1) = ( b(n-1) - a(4,n-1) * x(n) ) / a(3,n-1)

    do i = n - 2, 1, -1
       x(i) = ( b(i) - a(4,i) * x(i+1) - a(5,i) * x(i+2) ) / a(3,i)
    end do

    return
  end subroutine r85_np_fs
  subroutine r85_print ( n, a, title )

    !*****************************************************************************80
    !
    !! R85_PRINT prints an R85 matrix.
    !
    !  Discussion:
    !
    !    The R85 storage format represents a pentadiagonal matrix as a 5 
    !    by N array, in which each row corresponds to a diagonal, and 
    !    column locations are preserved.  Thus, the original matrix is 
    !    "collapsed" vertically into the array.
    !
    !  Example:
    !
    !    Here is how an R85 matrix of order 6 would be stored:
    !
    !       *   *  A13 A24 A35 A46
    !       *  A12 A23 A34 A45 A56
    !      A11 A22 A33 A44 A55 A66
    !      A21 A32 A43 A54 A65  *
    !      A31 A42 A53 A64  *   *
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license. 
    !
    !  Modified:
    !
    !    17 September 2003
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the order of the matrix.
    !    N must be positive.
    !
    !    Input, real(dp) A(5,N), the matrix.
    !
    !    Input, character ( len = * ) TITLE, a title.
    !
    implicit none

    integer(i4) n

    real(dp) a(5,n)
    character ( len = * ) title

    call r85_print_some ( n, a, 1, 1, n, n, title )

    return
  end subroutine r85_print
  subroutine r85_print_some ( n, a, ilo, jlo, ihi, jhi, title )

    !*****************************************************************************80
    !
    !! R85_PRINT_SOME prints some of an R85 matrix.
    !
    !  Discussion:
    !
    !    The R85 storage format represents a pentadiagonal matrix as a 5 
    !    by N array, in which each row corresponds to a diagonal, and 
    !    column locations are preserved.  Thus, the original matrix is 
    !    "collapsed" vertically into the array.
    !
    !  Example:
    !
    !    Here is how an R85 matrix of order 6 would be stored:
    !
    !       *   *  A13 A24 A35 A46
    !       *  A12 A23 A34 A45 A56
    !      A11 A22 A33 A44 A55 A66
    !      A21 A32 A43 A54 A65  *
    !      A31 A42 A53 A64  *   *
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license. 
    !
    !  Modified:
    !
    !    05 January 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the order of the matrix.
    !    N must be positive.
    !
    !    Input, real(dp) A(5,N), the matrix.
    !
    !    Input, integer(i4) ILO, JLO, IHI, JHI, the first row and
    !    column, and the last row and column, to be printed.
    !
    !    Input, character ( len = * ) TITLE, a title.
    !
    implicit none

    integer(i4), parameter :: incx = 5
    integer(i4) n

    real(dp) a(5,n)
    character ( len = 14 ) ctemp(incx)
    integer(i4) i
    integer(i4) i2hi
    integer(i4) i2lo
    integer(i4) ihi
    integer(i4) ilo
    integer(i4) inc
    integer(i4) j
    integer(i4) j2
    integer(i4) j2hi
    integer(i4) j2lo
    integer(i4) jhi
    integer(i4) jlo
    character ( len = * ) title

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
    !
    !  Print the columns of the matrix, in strips of 5.
    !
    do j2lo = jlo, jhi, incx

       j2hi = j2lo + incx - 1
       j2hi = min ( j2hi, n )
       j2hi = min ( j2hi, jhi )

       inc = j2hi + 1 - j2lo

       write ( *, '(a)' ) ' '

       do j = j2lo, j2hi
          j2 = j + 1 - j2lo
          write ( ctemp(j2), '(i7,7x)' ) j
       end do

       write ( *, '(''  Col:  '',5a14)' ) ( ctemp(j2), j2 = 1, inc )
       write ( *, '(a)' ) '  Row'
       write ( *, '(a)' ) '  ---'
       !
       !  Determine the range of the rows in this strip.
       !
       i2lo = max ( ilo, 1 )
       i2lo = max ( i2lo, j2lo - 2 )

       i2hi = min ( ihi, n )
       i2hi = min ( i2hi, j2hi + 2 )

       do i = i2lo, i2hi
          !
          !  Print out (up to) 5 entries in row I, that lie in the current strip.
          !
          do j2 = 1, inc

             j = j2lo - 1 + j2

             if ( 2 < i-j .or. 2 < j-i ) then
                ctemp(j2) = '              '
             else if ( j == i+2 ) then
                write ( ctemp(j2), '(g14.6)' ) a(1,j)
             else if ( j == i+1 ) then
                write ( ctemp(j2), '(g14.6)' ) a(2,j)
             else if ( j == i ) then
                write ( ctemp(j2), '(g14.6)' ) a(3,j)
             else if ( j == i-1 ) then
                write ( ctemp(j2), '(g14.6)' ) a(4,j)
             else if ( j == i-2 ) then
                write ( ctemp(j2), '(g14.6)' ) a(5,j)
             end if

          end do

          write ( *, '(i5,1x,5a14)' ) i, ( ctemp(j2), j2 = 1, inc )

       end do

    end do

    return
  end subroutine r85_print_some
  subroutine r8ge_fs ( n, a, b, info )

    !*****************************************************************************80
    !
    !! R8GE_FS factors and solves an R8GE system.
    !
    !  Discussion:
    !
    !    The R8GE storage format is used for a general M by N matrix.  A storage 
    !    space is made for each entry.  The two dimensional logical
    !    array can be thought of as a vector of M*N entries, starting with
    !    the M entries in the column 1, then the M entries in column 2
    !    and so on.  Considered as a vector, the entry A(I,J) is then stored
    !    in vector location I+(J-1)*M.
    !
    !    R8GE storage is used by LINPACK and LAPACK.
    !
    !    R8GE_FS does not save the LU factors of the matrix, and hence cannot
    !    be used to efficiently solve multiple linear systems, or even to
    !    factor A at one time, and solve a single linear system at a later time.
    !
    !    R8GE_FS uses partial pivoting, but no pivot vector is required.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license. 
    !
    !  Modified:
    !
    !    04 March 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the order of the matrix.
    !    N must be positive.
    !
    !    Input/output, real(dp) A(N,N).
    !    On input, A is the coefficient matrix of the linear system.
    !    On output, A is in unit upper triangular form, and
    !    represents the U factor of an LU factorization of the
    !    original coefficient matrix.
    !
    !    Input/output, real(dp) B(N).
    !    On input, B is the right hand side of the linear system.
    !    On output, B is the solution of the linear system.
    !
    !    Output, integer(i4) INFO, singularity flag.
    !    0, no singularity detected.
    !    nonzero, the factorization failed on the INFO-th step.
    !
    implicit none

    integer(i4) n

    real(dp) a(n,n)
    real(dp) b(n)
    integer(i4) i
    integer(i4) info
    integer(i4) ipiv
!    integer(i4) j
    integer(i4) jcol
    real(dp) piv
    real(dp) row(n)
    real(dp) temp

    info = 0

    do jcol = 1, n
       !
       !  Find the maximum element in column I.
       !
       piv = abs ( a(jcol,jcol) )
       ipiv = jcol
       do i = jcol + 1, n
          if ( piv < abs ( a(i,jcol) ) ) then
             piv = abs ( a(i,jcol) )
             ipiv = i
          end if
       end do

       if ( eq(piv,0.0_dp) ) then
          info = jcol
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'R8GE_FS - Fatal error!'
          write ( *, '(a,i8)' ) '  Zero pivot on step ', info
          stop
       end if
       !
       !  Switch rows JCOL and IPIV, and B.
       !
       if ( jcol /= ipiv ) then

          row(1:n) = a(jcol,1:n)
          a(jcol,1:n) = a(ipiv,1:n)
          a(ipiv,1:n) = row(1:n)

          temp = b(jcol)
          b(jcol) = b(ipiv)
          b(ipiv) = temp

       end if
       !
       !  Scale the pivot row.
       !
       a(jcol,jcol+1:n) = a(jcol,jcol+1:n) / a(jcol,jcol)
       b(jcol) = b(jcol) / a(jcol,jcol)
       a(jcol,jcol) = 1.0_dp
       !
       !  Use the pivot row to eliminate lower entries in that column.
       !
       do i = jcol + 1, n
          if ( ne(a(i,jcol),0.0_dp) ) then
             temp = - a(i,jcol)
             a(i,jcol) = 0.0_dp
             a(i,jcol+1:n) = a(i,jcol+1:n) + temp * a(jcol,jcol+1:n)
             b(i) = b(i) + temp * b(jcol)
          end if
       end do

    end do
    !
    !  Back solve.
    !
    do jcol = n, 2, -1
       b(1:jcol-1) = b(1:jcol-1) - a(1:jcol-1,jcol) * b(jcol)
    end do

    return
  end subroutine r8ge_fs
  subroutine r8vec_bracket ( n, x, xval, left, right )

    !*****************************************************************************80
    !
    !! R8VEC_BRACKET searches a sorted R8VEC for successive brackets of a value.
    !
    !  Discussion:
    !
    !    An R8VEC is an array of double precision real values.
    !
    !    If the values in the vector are thought of as defining intervals
    !    on the real line, then this routine searches for the interval
    !    nearest to or containing the given value.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    06 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, length of input array.
    !
    !    Input, real(dp) X(N), an array sorted into ascending order.
    !
    !    Input, real(dp) XVAL, a value to be bracketed.
    !
    !    Output, integer(i4) LEFT, RIGHT, the results of the search.
    !    Either:
    !      XVAL < X(1), when LEFT = 1, RIGHT = 2;
    !      X(N) < XVAL, when LEFT = N-1, RIGHT = N;
    !    or
    !      X(LEFT) <= XVAL <= X(RIGHT).
    !
    implicit none

    integer(i4) n

    integer(i4) i
    integer(i4) left
    integer(i4) right
    real(dp) x(n)
    real(dp) xval

    do i = 2, n - 1

       if ( xval < x(i) ) then
          left = i - 1
          right = i
          return
       end if

    end do

    left = n - 1
    right = n

    return
  end subroutine r8vec_bracket
  subroutine r8vec_bracket3 ( n, t, tval, left )

    !*****************************************************************************80
    !
    !! R8VEC_BRACKET3 finds the interval containing or nearest a given value.
    !
    !  Discussion:
    !
    !    An R8VEC is an array of double precision real values.
    !
    !    The routine always returns the index LEFT of the sorted array
    !    T with the property that either
    !    *  T is contained in the interval [ T(LEFT), T(LEFT+1) ], or
    !    *  T < T(LEFT) = T(1), or
    !    *  T > T(LEFT+1) = T(N).
    !
    !    The routine is useful for interpolation problems, where
    !    the abscissa must be located within an interval of data
    !    abscissas for interpolation, or the "nearest" interval
    !    to the (extreme) abscissa must be found so that extrapolation
    !    can be carried out.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    05 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, length of the input array.
    !
    !    Input, real(dp) T(N), an array sorted into ascending order.
    !
    !    Input, real(dp) TVAL, a value to be bracketed by entries of T.
    !
    !    Input/output, integer(i4) LEFT.
    !
    !    On input, if 1 <= LEFT <= N-1, LEFT is taken as a suggestion for the
    !    interval [ T(LEFT), T(LEFT+1) ] in which TVAL lies.  This interval
    !    is searched first, followed by the appropriate interval to the left
    !    or right.  After that, a binary search is used.
    !
    !    On output, LEFT is set so that the interval [ T(LEFT), T(LEFT+1) ]
    !    is the closest to TVAL; it either contains TVAL, or else TVAL
    !    lies outside the interval [ T(1), T(N) ].
    !
    implicit none

    integer(i4) n

    integer(i4) high
    integer(i4) left
    integer(i4) low
    integer(i4) mid
    real(dp) t(n)
    real(dp) tval
    !
    !  Check the input data.
    !
    if ( n < 2 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'R8VEC_BRACKET3 - Fatal error!'
       write ( *, '(a)' ) '  N must be at least 2.'
       stop
    end if
    !
    !  If LEFT is not between 1 and N-1, set it to the middle value.
    !
    if ( left < 1 .or. n - 1 < left ) then
       left = ( n + 1 ) / 2
    end if
    !
    !  CASE 1: TVAL < T(LEFT):
    !  Search for TVAL in [T(I), T(I+1)] for intervals I = 1 to LEFT-1.
    !
    if ( tval < t(left) ) then

       if ( left == 1 ) then
          return
       else if ( left == 2 ) then
          left = 1
          return
       else if ( le(t(left-1),tval) ) then
          left = left - 1
          return
       else if ( le(tval,t(2)) ) then
          left = 1
          return
       end if
       !
       !  ...Binary search for TVAL in [T(I), T(I+1)] for intervals I = 2 to LEFT-2.
       !
       low = 2
       high = left - 2

       do

          if ( low == high ) then
             left = low
             return
          end if

          mid = ( low + high + 1 ) / 2

          if ( le(t(mid),tval) ) then
             low = mid
          else
             high = mid - 1
          end if

       end do
       !
       !  CASE2: T(LEFT+1) < TVAL:
       !  Search for TVAL in {T(I),T(I+1)] for intervals I = LEFT+1 to N-1.
       !
    else if ( t(left+1) < tval ) then

       if ( left == n - 1 ) then
          return
       else if ( left == n - 2 ) then
          left = left + 1
          return
       else if ( le(tval,t(left+2)) ) then
          left = left + 1
          return
       else if ( le(t(n-1),tval) ) then
          left = n - 1
          return
       end if
       !
       !  ...Binary search for TVAL in [T(I), T(I+1)] for intervals I = LEFT+2 to N-2.
       !
       low = left + 2
       high = n - 2

       do

          if ( low == high ) then
             left = low
             return
          end if

          mid = ( low + high + 1 ) / 2

          if ( le(t(mid),tval) ) then
             low = mid
          else
             high = mid - 1
          end if

       end do
       !
       !  CASE3: T(LEFT) <= TVAL <= T(LEFT+1):
       !  T is in [T(LEFT), T(LEFT+1)], as the user said it might be.
       !
    else

    end if

    return
  end subroutine r8vec_bracket3
  function r8vec_distinct ( n, x )

    !*****************************************************************************80
    !
    !! R8VEC_DISTINCT is true if the entries in an R8VEC are distinct.
    !
    !  Discussion:
    !
    !    An R8VEC is an array of double precision real values.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    16 September 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of entries in the vector.
    !
    !    Input, real(dp) X(N), the vector to be checked.
    !
    !    Output, logical R8VEC_DISTINCT is TRUE if all N elements of X
    !    are distinct.
    !
    implicit none

    integer(i4) n

    integer(i4) i
    integer(i4) j
    logical r8vec_distinct
    real(dp) x(n)

    r8vec_distinct = .false.

    do i = 2, n
       do j = 1, i - 1
          if ( eq(x(i),x(j)) ) then
             return
          end if
       end do
    end do

    r8vec_distinct = .true.

    return
  end function r8vec_distinct
  subroutine r8vec_even ( n, alo, ahi, a )

    !*****************************************************************************80
    !
    !! R8VEC_EVEN returns N real values, evenly spaced between ALO and AHI.
    !
    !  Discussion:
    !
    !    An R8VEC is an array of double precision real values.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    17 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of values.
    !
    !    Input, real(dp) ALO, AHI, the low and high values.
    !
    !    Output, real(dp) A(N), N evenly spaced values.
    !    Normally, A(1) = ALO and A(N) = AHI.
    !    However, if N = 1, then A(1) = 0.5*(ALO+AHI).
    !
    implicit none

    integer(i4) n

    real(dp) a(n)
    real(dp) ahi
    real(dp) alo
    integer(i4) i

    if ( n == 1 ) then

       a(1) = 0.5_dp * ( alo + ahi )

    else

       do i = 1, n
          a(i) = ( real ( n - i,     kind=dp ) * alo   &
               + real (     i - 1, kind=dp ) * ahi ) &
               / real ( n     - 1, kind=dp )
       end do

    end if

    return
  end subroutine r8vec_even
  subroutine r8vec_indicator ( n, a )

    !*****************************************************************************80
    !
    !! R8VEC_INDICATOR sets an R8VEC to the indicator vector.
    !
    !  Discussion:
    !
    !    An R8VEC is an array of double precision real values.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    09 February 2001
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of elements of A.
    !
    !    Output, real(dp) A(N), the array to be initialized.
    !
    implicit none

    integer(i4) n

    real(dp) a(n)
    integer(i4) i

    do i = 1, n
       a(i) = real ( i, kind=dp )
    end do

    return
  end subroutine r8vec_indicator
  subroutine r8vec_order_type ( n, a, order )

    !*****************************************************************************80
    !
    !! R8VEC_ORDER_TYPE determines the order type of an R8VEC.
    !
    !  Discussion:
    !
    !    An R8VEC is an array of double precision real values.
    !
    !    We assume the array is increasing or decreasing, and we want to
    !    verify that.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    20 July 2000
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of entries of the array.
    !
    !    Input, real(dp) A(N), the array to be checked.
    !
    !    Output, integer(i4) ORDER, order indicator:
    !    -1, no discernable order;
    !    0, all entries are equal;
    !    1, ascending order;
    !    2, strictly ascending order;
    !    3, descending order;
    !    4, strictly descending order.
    !
    implicit none

    integer(i4) n

    real(dp) a(n)
    integer(i4) i
    integer(i4) order
    !
    !  Search for the first value not equal to A(1).
    !
    i = 1

    do

       i = i + 1

       if ( n < i ) then
          order = 0
          return
       end if

       if ( a(1) < a(i) ) then

          if ( i == 2 ) then
             order = 2
          else
             order = 1
          end if

          exit

       else if ( a(i) < a(1) ) then

          if ( i == 2 ) then
             order = 4
          else
             order = 3
          end if

          exit

       end if

    end do
    !
    !  Now we have a "direction".  Examine subsequent entries.
    !
    do

       i = i + 1
       if ( n < i ) then
          exit
       end if

       if ( order == 1 ) then

          if ( a(i) < a(i-1) ) then
             order = -1
             exit
          end if

       else if ( order == 2 ) then

          if ( a(i) < a(i-1) ) then
             order = -1
             exit
          else if ( eq(a(i),a(i-1)) ) then
             order = 1
          end if

       else if ( order == 3 ) then

          if ( a(i-1) < a(i) ) then
             order = -1
             exit
          end if

       else if ( order == 4 ) then

          if ( a(i-1) < a(i) ) then
             order = -1
             exit
          else if ( eq(a(i),a(i-1)) ) then
             order = 3
          end if

       end if

    end do

    return
  end subroutine r8vec_order_type
  subroutine r8vec_print ( n, a, title )

    !*****************************************************************************80
    !
    !! R8VEC_PRINT prints an R8VEC.
    !
    !  Discussion:
    !
    !    An R8VEC is an array of double precision real values.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    16 December 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of components of the vector.
    !
    !    Input, real(dp) A(N), the vector to be printed.
    !
    !    Input, character ( len = * ) TITLE, a title.
    !
    implicit none

    integer(i4) n

    real(dp) a(n)
    integer(i4) i
    character ( len = * ) title

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
    write ( *, '(a)' ) ' '
    do i = 1, n
       write ( *, '(i8,g14.6)' ) i, a(i)
    end do

    return
  end subroutine r8vec_print
  subroutine r8vec_sort_bubble_a ( n, a )

    !*****************************************************************************80
    !
    !! R8VEC_SORT_BUBBLE_A ascending sorts an R8VEC using bubble sort.
    !
    !  Discussion:
    !
    !    An R8VEC is an array of double precision real values.
    !
    !    Bubble sort is simple to program, but inefficient.  It should not
    !    be used for large arrays.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    01 February 2001
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of entries in the array.
    !
    !    Input/output, real(dp) A(N).
    !    On input, an unsorted array.
    !    On output, the array has been sorted.
    !
    implicit none

    integer(i4) n

    real(dp) a(n)
    integer(i4) i
    integer(i4) j

    do i = 1, n-1
       do j = i+1, n
          if ( a(j) < a(i) ) then
             call r8_swap ( a(i), a(j) )
          end if
       end do
    end do

    return
  end subroutine r8vec_sort_bubble_a
  subroutine r8vec_uniform_01 ( n, seed, r )

    !*****************************************************************************80
    !
    !! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
    !
    !  Discussion:
    !
    !    An R8VEC is an array of double precision real values.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    19 August 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Paul Bratley, Bennett Fox, Linus Schrage,
    !    A Guide to Simulation,
    !    Springer Verlag, pages 201-202, 1983.
    !
    !    Bennett Fox,
    !    Algorithm 647:
    !    Implementation and Relative Efficiency of Quasirandom
    !    Sequence Generators,
    !    ACM Transactions on Mathematical Software,
    !    Volume 12, Number 4, pages 362-376, 1986.
    !
    !    Peter Lewis, Allen Goodman, James Miller,
    !    A Pseudo-Random Number Generator for the System/360,
    !    IBM Systems Journal,
    !    Volume 8, pages 136-143, 1969.
    !
    !  Parameters:
    !
    !    Input, integer(i4) M, the number of entries in the vector.
    !
    !    Input/output, integer(i4) SEED, the "seed" value, which should
    !    NOT be 0.  On output, SEED has been updated.
    !
    !    Output, real(dp) R(N), the vector of pseudorandom values.
    !
    implicit none

    integer(i4) n

    integer(i4) i
    integer(i4) k
    integer(i4) seed
    real(dp) r(n)

    do i = 1, n

       k = seed / 127773

       seed = 16807 * ( seed - k * 127773 ) - k * 2836

       if ( seed < 0 ) then
          seed = seed + 2147483647
       end if

       r(i) = real ( seed, kind=dp ) * 4.656612875D-10

    end do

    return
  end subroutine r8vec_uniform_01
  subroutine r8vec_unique_count ( n, a, tol, unique_num )

    !*****************************************************************************80
    !
    !! R8VEC_UNIQUE_COUNT counts the unique elements in an unsorted R8VEC.
    !
    !  Discussion:
    !
    !    An R8VEC is an array of double precision real values.
    !
    !    Because the array is unsorted, this algorithm is O(N^2).
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    08 December 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of elements of A.
    !
    !    Input, real(dp) A(N), the unsorted array to examine.
    !
    !    Input, real(dp) TOL, a nonnegative tolerance for equality.
    !    Set it to 0.0 for the strictest test.
    !
    !    Output, integer(i4) UNIQUE_NUM, the number of unique elements.
    !
    implicit none

    integer(i4) n

    real(dp) a(n)
    integer(i4) i
    integer(i4) j
    integer(i4) unique_num
    real(dp) tol

    unique_num = 0

    do i = 1, n

       unique_num = unique_num + 1

       do j = 1, i - 1

          if ( le(abs( a(i) - a(j) ),tol) ) then
             unique_num = unique_num - 1
             exit
          end if

       end do

    end do

    return
  end subroutine r8vec_unique_count
  subroutine spline_b_val ( ndata, tdata, ydata, tval, yval )

    !*****************************************************************************80
    !
    !! SPLINE_B_VAL evaluates a cubic B spline approximant.
    !
    !  Discussion:
    !
    !    The cubic B spline will approximate the data, but is not
    !    designed to interpolate it.
    !
    !    In effect, two "phantom" data values are appended to the data,
    !    so that the spline will interpolate the first and last data values.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    11 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Carl deBoor,
    !    A Practical Guide to Splines,
    !    Springer, 2001,
    !    ISBN: 0387953663.
    !
    !  Parameters:
    !
    !    Input, integer(i4) NDATA, the number of data values.
    !
    !    Input, real(dp) TDATA(NDATA), the abscissas of the data.
    !
    !    Input, real(dp) YDATA(NDATA), the data values.
    !
    !    Input, real(dp) TVAL, a point at which the spline is
    !    to be evaluated.
    !
    !    Output, real(dp) YVAL, the value of the function at TVAL.
    !
    implicit none

    integer(i4) ndata

    real(dp) bval
    integer(i4) left
    integer(i4) right
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) u
    real(dp) ydata(ndata)
    real(dp) yval
    !
    !  Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
    !
    call r8vec_bracket ( ndata, tdata, tval, left, right )
    !
    !  Evaluate the 5 nonzero B spline basis functions in the interval,
    !  weighted by their corresponding data values.
    !
    u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) )
    yval = 0.0_dp
    !
    !  B function associated with node LEFT - 1, (or "phantom node"),
    !  evaluated in its 4th interval.
    !
    bval = ( ( (     - 1.0_dp   &
         * u + 3.0_dp ) &
         * u - 3.0_dp ) &
         * u + 1.0_dp ) / 6.0_dp

    if ( 0 < left-1 ) then
       yval = yval + ydata(left-1) * bval
    else
       yval = yval + ( 2.0_dp * ydata(1) - ydata(2) ) * bval
    end if
    !
    !  B function associated with node LEFT,
    !  evaluated in its third interval.
    !
    bval = ( ( (       3.0_dp   &
         * u - 6.0_dp ) &
         * u + 0.0_dp ) &
         * u + 4.0_dp ) / 6.0_dp

    yval = yval + ydata(left) * bval
    !
    !  B function associated with node RIGHT,
    !  evaluated in its second interval.
    !
    bval = ( ( (     - 3.0_dp   &
         * u + 3.0_dp ) &
         * u + 3.0_dp ) &
         * u + 1.0_dp ) / 6.0_dp

    yval = yval + ydata(right) * bval
    !
    !  B function associated with node RIGHT+1, (or "phantom node"),
    !  evaluated in its first interval.
    !
    bval = u**3 / 6.0_dp

    if ( right+1 <= ndata ) then
       yval = yval + ydata(right+1) * bval
    else
       yval = yval + ( 2.0_dp * ydata(ndata) - ydata(ndata-1) ) * bval
    end if

    return
  end subroutine spline_b_val


  subroutine spline_beta_val ( beta1, beta2, ndata, tdata, ydata, tval, yval )

    !*****************************************************************************80
    !
    !! SPLINE_BETA_VAL evaluates a cubic beta spline approximant.
    !
    !  Discussion:
    !
    !    The cubic beta spline will approximate the data, but is not
    !    designed to interpolate it.
    !
    !    If BETA1 = 1 and BETA2 = 0, the cubic beta spline will be the
    !    same as the cubic B spline approximant.
    !
    !    With BETA1 = 1 and BETA2 large, the beta spline becomes more like
    !    a linear spline.
    !
    !    In effect, two "phantom" data values are appended to the data,
    !    so that the spline will interpolate the first and last data values.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    11 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, real(dp) BETA1, the skew or bias parameter.
    !    BETA1 = 1 for no skew or bias.
    !
    !    Input, real(dp) BETA2, the tension parameter.
    !    BETA2 = 0 for no tension.
    !
    !    Input, integer(i4) NDATA, the number of data values.
    !
    !    Input, real(dp) TDATA(NDATA), the abscissas of the data.
    !
    !    Input, real(dp) YDATA(NDATA), the data values.
    !
    !    Input, real(dp) TVAL, a point at which the spline is
    !    to be evaluated.
    !
    !    Output, real(dp) YVAL, the value of the function at TVAL.
    !
    implicit none

    integer(i4) ndata

    real(dp) a
    real(dp) b
    real(dp) beta1
    real(dp) beta2
    real(dp) bval
    real(dp) c
    real(dp) d
    real(dp) delta
    integer(i4) left
    integer(i4) right
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) u
    real(dp) ydata(ndata)
    real(dp) yval
    !
    !  Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
    !
    call r8vec_bracket ( ndata, tdata, tval, left, right )
    !
    !  Evaluate the 5 nonzero beta spline basis functions in the interval,
    !  weighted by their corresponding data values.
    !
    u = ( tval - tdata(left) ) / ( tdata(right) - tdata(left) )

    delta = ( ( 2.0_dp   &
         * beta1 + 4.0_dp ) &
         * beta1 + 4.0_dp ) &
         * beta1 + 2.0_dp + beta2

    yval = 0.0_dp
    !
    !  Beta function associated with node LEFT - 1, (or "phantom node"),
    !  evaluated in its 4th interval.
    !
    bval = 2.0_dp * ( beta1 * ( 1.0_dp - u ) )**3  / delta

    if ( 0 < left-1 ) then
       yval = yval + ydata(left-1) * bval
    else
       yval = yval + ( 2.0_dp * ydata(1) - ydata(2) ) * bval
    end if
    !
    !  Beta function associated with node LEFT,
    !  evaluated in its third interval.
    !
    a = beta2 + ( 4.0_dp + 4.0_dp * beta1 ) * beta1

    b = - 6.0_dp * beta1 * ( 1.0_dp - beta1 ) * ( 1.0_dp + beta1 )

    c = ( (     - 6.0_dp   &
         * beta1 - 6.0_dp ) &
         * beta1 + 0.0_dp ) &
         * beta1 - 3.0_dp * beta2

    d = ( (     + 2.0_dp   &
         * beta1 + 2.0_dp ) &
         * beta1 + 2.0_dp ) &
         * beta1 + 2.0_dp * beta2

    bval = ( a + u * ( b + u * ( c + u * d ) ) ) / delta

    yval = yval + ydata(left) * bval
    !
    !  Beta function associated with node RIGHT,
    !  evaluated in its second interval.
    !
    a = 2.0_dp

    b = 6.0_dp * beta1

    c = 3.0_dp * beta2 + 6.0_dp * beta1 * beta1

    d = - 2.0_dp * ( 1.0_dp + beta2 + beta1 + beta1 * beta1 )

    bval = ( a + u * ( b + u * ( c + u * d ) ) ) / delta

    yval = yval + ydata(right) * bval
    !
    !  Beta function associated with node RIGHT+1, (or "phantom node"),
    !  evaluated in its first interval.
    !
    bval = 2.0_dp * u**3 / delta

    if ( right + 1 <= ndata ) then
       yval = yval + ydata(right+1) * bval
    else
       yval = yval + ( 2.0_dp * ydata(ndata) - ydata(ndata-1) ) * bval
    end if

    return
  end subroutine spline_beta_val
  subroutine spline_bezier_val ( dim_num, interval_num, data_val, point_num, &
       point_t, point_val )

    !*****************************************************************************80
    !
    !! SPLINE_BEZIER_VAL evaluates a cubic Bezier spline.
    !
    !  Discussion:
    !
    !    The cubic Bezier spline of N parts is defined by choosing
    !    3*N+1 equally spaced T-abscissa values in the interval [0,N].
    !    This defines N subintervals, each of length 1, and each containing
    !    4 successives T abscissa values.
    !
    !    At each abscissa value, a DIM_NUM-dimensional Bezier control
    !    value is assigned.  Over each interval, a Bezier cubic function
    !    is used to define the value of the Bezier spline.  To the left of
    !    the first interval, or to the right of the last interval,
    !    extrapolation may be used to extend the spline definition to
    !    the entire real line.
    !
    !    Note that the Bezier spline will pass through the 1st, 4th,
    !    and in general 3*I+1 control values exactly.  The other control
    !    values are not interpolating points.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    17 June 2006
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) DIM_NUM, the spatial dimension.
    !
    !    Input, integer(i4) INTERVAL_NUM, the number of intervals.
    !
    !    Input, real(dp) DATA_VAL(DIM_NUM,3*INTERVAL_NUM+1), the control
    !    values.
    !
    !    Input, integer(i4) POINT_NUM, the number of sample points at
    !    which the Bezier cubic spline is to be evaluated.
    !
    !    Input, real(dp) POINT_T(POINT_NUM), the "T" values associated
    !    with the points.  A value of T between 0 and 1, for instance,
    !    is associated with the first interval, and a value of T between
    !    INTERVAL_NUM-1 and INTERVAL_NUM is in the last interval.
    !
    !    Output, real(dp) POINT_VAL(DIM_NUM,POINT_NUM), the value
    !    of the Bezier cubic spline at the sample points.
    !
    implicit none

    integer(i4), parameter :: cubic = 3
    integer(i4) interval_num
    integer(i4) dim_num
    integer(i4) point_num

    real(dp) bernstein_val(0:cubic)
    real(dp) data_val(dim_num,cubic*interval_num+1)
    integer(i4) dim
    integer(i4) interval
    integer(i4) offset
    integer(i4) point
    real(dp) point_t(point_num)
    real(dp) point_val(dim_num,point_num)
    real(dp) t
    real(dp) t_01

    do point = 1, point_num

       t = point_t(point)

       interval = int ( t + 1 )

       interval = max ( interval, 1 )
       interval = min ( interval, interval_num )

       offset = 1 + ( interval - 1 ) * cubic

       t_01 = t - real ( interval - 1, kind=dp )

       call bp01 ( cubic, t_01, bernstein_val )

       do dim = 1, dim_num
          point_val(dim,point) = dot_product ( &
               data_val(dim,offset:offset+cubic), bernstein_val(0:cubic) )
       end do

    end do

    return
  end subroutine spline_bezier_val
  subroutine spline_constant_val ( ndata, tdata, ydata, tval, yval )

    !*****************************************************************************80
    !
    !! SPLINE_CONSTANT_VAL evaluates a piecewise constant spline at a point.
    !
    !  Discussion:
    !
    !    NDATA-1 points TDATA define NDATA intervals, with the first
    !    and last being semi-infinite.
    !
    !    The value of the spline anywhere in interval I is YDATA(I).
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    16 November 2001
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) NDATA, the number of data points defining
    !    the spline.  NDATA must be at least 1.
    !
    !    Input, real(dp) TDATA(NDATA-1), the breakpoints.  The values
    !    of TDATA should be distinct and increasing.
    !
    !    Input, real(dp) YDATA(NDATA), the values of the spline in
    !    the intervals defined by the breakpoints.
    !
    !    Input, real(dp) TVAL, the point at which the spline is
    !    to be evaluated.
    !
    !    Output, real(dp) YVAL, the value of the spline at TVAL.
    !
    implicit none

    integer(i4) ndata

    integer(i4) i
    real(dp) tdata(ndata-1)
    real(dp) tval
    real(dp) ydata(ndata)
    real(dp) yval
    !
    !  Check NDATA.
    !
    if ( ndata < 1 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_CONSTANT_VAL - Fatal error!'
       write ( *, '(a)' ) '  NDATA < 1.'
       stop
    end if

    do i = 1, ndata-1
       if ( le(tval,tdata(i)) ) then
          yval = ydata(i)
          return
       end if
    end do

    yval = ydata(ndata)

    return
  end subroutine spline_constant_val
  subroutine spline_cubic_set_full ( n, t, y, ibcbeg, ybcbeg, ibcend, ybcend, ypp )

    !*****************************************************************************80
    !
    !! SPLINE_CUBIC_SET_FULL computes the second derivatives of a piecewise cubic spline.
    !
    !  Discussion:
    !
    !    DELETE THIS FUNCTION WHEN YOU ARE SATISFIED WITH THE CURRENT
    !    VERSION OF SPLINE_CUBIC_SET.
    !
    !    For data interpolation, the user must call SPLINE_CUBIC_SET to
    !    determine the second derivative data, passing in the data to be
    !    interpolated, and the desired boundary conditions.
    !
    !    The data to be interpolated, plus the SPLINE_CUBIC_SET output,
    !    defines the spline.  The user may then call SPLINE_CUBIC_VAL to
    !    evaluate the spline at any point.
    !
    !    The cubic spline is a piecewise cubic polynomial.  The intervals
    !    are determined by the "knots" or abscissas of the data to be
    !    interpolated.  The cubic spline has continous first and second
    !    derivatives over the entire interval of interpolation.
    !
    !    For any point T in the interval T(IVAL), T(IVAL+1), the form of
    !    the spline is
    !
    !      SPL(T) = A(IVAL)
    !             + B(IVAL) * ( T - T(IVAL) )
    !             + C(IVAL) * ( T - T(IVAL) )^2
    !             + D(IVAL) * ( T - T(IVAL) )^3
    !
    !    If we assume that we know the values Y(*) and YPP(*), which represent
    !    the values and second derivatives of the spline at each knot, then
    !    the coefficients can be computed as:
    !
    !      A(IVAL) = Y(IVAL)
    !      B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
    !        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
    !      C(IVAL) = YPP(IVAL) / 2
    !      D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
    !
    !    Since the first derivative of the spline is
    !
    !      SPL'(T) =     B(IVAL)
    !              + 2 * C(IVAL) * ( T - T(IVAL) )
    !              + 3 * D(IVAL) * ( T - T(IVAL) )^2,
    !
    !    the requirement that the first derivative be continuous at interior
    !    knot I results in a total of N-2 equations, of the form:
    !
    !      B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
    !      + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))^2 = B(IVAL)
    !
    !    or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
    !
    !      ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
    !      - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
    !      + YPP(IVAL-1) * H(IVAL-1)
    !      + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
    !      =
    !      ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
    !      - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
    !
    !    or
    !
    !      YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
    !      + YPP(IVAL) * H(IVAL)
    !      =
    !      6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
    !      - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
    !
    !    Boundary conditions must be applied at the first and last knots.
    !    The resulting tridiagonal system can be solved for the YPP values.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    06 June 2013
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Carl deBoor,
    !    A Practical Guide to Splines,
    !    Springer, 2001,
    !    ISBN: 0387953663.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of data points; N must be
    !    at least 2.
    !
    !    Input, real(dp) T(N), the points where data is specified.
    !    The values should be distinct, and increasing.
    !
    !    Input, real(dp) Y(N), the data values to be interpolated.
    !
    !    Input, integer(i4) IBCBEG, the left boundary condition flag:
    !    0: the spline should be a quadratic over the first interval;
    !    1: the first derivative at the left endpoint should be YBCBEG;
    !    2: the second derivative at the left endpoint should be YBCBEG;
    !    3: Not-a-knot: the third derivative is continuous at T(2).
    !
    !    Input, real(dp) YBCBEG, the left boundary value, if needed.
    !
    !    Input, integer(i4) IBCEND, the right boundary condition flag:
    !    0: the spline should be a quadratic over the last interval;
    !    1: the first derivative at the right endpoint should be YBCEND;
    !    2: the second derivative at the right endpoint should be YBCEND;
    !    3: Not-a-knot: the third derivative is continuous at T(N-1).
    !
    !    Input, real(dp) YBCEND, the right boundary value, if needed.
    !
    !    Output, real(dp) YPP(N), the second derivatives of
    !    the cubic spline.
    !
    implicit none

    integer(i4) n

    real(dp) a(n,n)
    real(dp) b(n)
    integer(i4) i
    integer(i4) ibcbeg
    integer(i4) ibcend
    integer(i4) info
    real(dp) t(n)
    real(dp) y(n)
    real(dp) ybcbeg
    real(dp) ybcend
    real(dp) ypp(n)
    !
    !  Check.
    !
    if ( n <= 1 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
       write ( *, '(a)' ) '  The number of knots must be at least 2.'
       write ( *, '(a,i8)' ) '  The input value of N = ', n
       stop
    end if

    do i = 1, n - 1
       if ( le(t(i+1),t(i)) ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
          write ( *, '(a)' ) '  The knots must be strictly increasing, but'
          write ( *, '(a,i8,a,g14.6)' ) '  T(',  i,') = ', t(i)
          write ( *, '(a,i8,a,g14.6)' ) '  T(',i+1,') = ', t(i+1)
          stop
       end if
    end do
    !
    !  Zero out the matrix.
    !
    a(1:n,1:n) = 0.0_dp
    !
    !  Set the first equation.
    !
    if ( ibcbeg == 0 ) then
       b(1) = 0.0_dp
       a(1,1) =  1.0_dp
       a(1,2) = -1.0_dp
    else if ( ibcbeg == 1 ) then
       b(1) = ( y(2) - y(1) ) / ( t(2) - t(1) ) - ybcbeg
       a(1,1) = ( t(2) - t(1) ) / 3.0_dp
       a(1,2) = ( t(2) - t(1) ) / 6.0_dp
    else if ( ibcbeg == 2 ) then
       b(1) = ybcbeg
       a(1,1) = 1.0_dp
       a(1,2) = 0.0_dp
    else if ( ibcbeg == 3 ) then
       b(1) = 0.0_dp
       a(1,1) = - ( t(3) - t(2) )
       a(1,2) =   ( t(3)        - t(1) )
       a(1,3) = - (        t(2) - t(1) )
    else
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
       write ( *, '(a)' ) '  The boundary flag IBCBEG must be 0, 1, 2 or 3.'
       write ( *, '(a,i8)' ) '  The input value is IBCBEG = ', ibcbeg
       stop
    end if
    !
    !  Set the intermediate equations.
    !
    do i = 2, n - 1
       b(i) = ( y(i+1) - y(i) ) / ( t(i+1) - t(i) ) &
            - ( y(i) - y(i-1) ) / ( t(i) - t(i-1) )
       a(i,i-1) = ( t(i)   - t(i-1) ) / 6.0_dp
       a(i,i)   = ( t(i+1) - t(i-1) ) / 3.0_dp
       a(i,i+1) = ( t(i+1) - t(i)   ) / 6.0_dp
    end do
    !
    !  Set the last equation.
    !
    if ( ibcend == 0 ) then
       b(n) = 0.0_dp
       a(n,n-1) = -1.0_dp
       a(n,n) = 1.0_dp
    else if ( ibcend == 1 ) then
       b(n) = ybcend - ( y(n) - y(n-1) ) / ( t(n) - t(n-1) )
       a(n,n-1) = ( t(n) - t(n-1) ) / 6.0_dp
       a(n,n) = ( t(n) - t(n-1) ) / 3.0_dp
    else if ( ibcend == 2 ) then
       b(n) = ybcend
       a(n,n-1) = 0.0_dp
       a(n,n) = 1.0_dp
    else if ( ibcend == 3 ) then
       b(n) = 0.0_dp
       a(n,n-2) = - ( t(n) - t(n-1) )
       a(n,n-1) =   ( t(n)          - t(n-2) )
       a(n,n) =   - (        t(n-1) - t(n-2) )
    else
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
       write ( *, '(a)' ) '  The boundary flag IBCEND must be 0, 1, 2 or 3.'
       write ( *, '(a,i8)' ) '  The input value is IBCEND = ', ibcend
       stop
    end if
    !
    !  Special case:
    !    N = 2, IBCBEG = IBCEND = 0.
    !
    if ( n == 2 .and. ibcbeg == 0 .and. ibcend == 0 ) then

       ypp(1) = 0.0_dp
       ypp(2) = 0.0_dp
       !
       !  Solve the linear system.
       !
    else

       call r8ge_fs ( n, a, b, info )
       ypp = b

    end if

    return
  end subroutine spline_cubic_set_full
  subroutine spline_cubic_set ( n, t, y, ibcbeg, ybcbeg, ibcend, ybcend, ypp )

    !*****************************************************************************80
    !
    !! SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic spline.
    !
    !  Discussion:
    !
    !    For data interpolation, the user must call SPLINE_CUBIC_SET to
    !    determine the second derivative data, passing in the data to be
    !    interpolated, and the desired boundary conditions.
    !
    !    The data to be interpolated, plus the SPLINE_CUBIC_SET output,
    !    defines the spline.  The user may then call SPLINE_CUBIC_VAL to
    !    evaluate the spline at any point.
    !
    !    The cubic spline is a piecewise cubic polynomial.  The intervals
    !    are determined by the "knots" or abscissas of the data to be
    !    interpolated.  The cubic spline has continous first and second
    !    derivatives over the entire interval of interpolation.
    !
    !    For any point T in the interval T(IVAL), T(IVAL+1), the form of
    !    the spline is
    !
    !      SPL(T) = A(IVAL)
    !             + B(IVAL) * ( T - T(IVAL) )
    !             + C(IVAL) * ( T - T(IVAL) )^2
    !             + D(IVAL) * ( T - T(IVAL) )^3
    !
    !    If we assume that we know the values Y(*) and YPP(*), which represent
    !    the values and second derivatives of the spline at each knot, then
    !    the coefficients can be computed as:
    !
    !      A(IVAL) = Y(IVAL)
    !      B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
    !        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
    !      C(IVAL) = YPP(IVAL) / 2
    !      D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
    !
    !    Since the first derivative of the spline is
    !
    !      SPL'(T) =     B(IVAL)
    !              + 2 * C(IVAL) * ( T - T(IVAL) )
    !              + 3 * D(IVAL) * ( T - T(IVAL) )^2,
    !
    !    the requirement that the first derivative be continuous at interior
    !    knot I results in a total of N-2 equations, of the form:
    !
    !      B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
    !      + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))^2 = B(IVAL)
    !
    !    or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
    !
    !      ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
    !      - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
    !      + YPP(IVAL-1) * H(IVAL-1)
    !      + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
    !      =
    !      ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
    !      - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
    !
    !    or
    !
    !      YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
    !      + YPP(IVAL) * H(IVAL)
    !      =
    !      6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
    !      - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
    !
    !    Boundary conditions must be applied at the first and last knots.
    !    The resulting tridiagonal system can be solved for the YPP values.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    07 June 2013
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Carl deBoor,
    !    A Practical Guide to Splines,
    !    Springer, 2001,
    !    ISBN: 0387953663.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of data points; N must be
    !    at least 2.
    !
    !    Input, real(dp) T(N), the points where data is specified.
    !    The values should be distinct, and increasing.
    !
    !    Input, real(dp) Y(N), the data values to be interpolated.
    !
    !    Input, integer(i4) IBCBEG, the left boundary condition flag:
    !    0: the spline should be a quadratic over the first interval;
    !    1: the first derivative at the left endpoint should be YBCBEG;
    !    2: the second derivative at the left endpoint should be YBCBEG;
    !    3: Not-a-knot: the third derivative is continuous at T(2).
    !
    !    Input, real(dp) YBCBEG, the left boundary value, if needed.
    !
    !    Input, integer(i4) IBCEND, the right boundary condition flag:
    !    0: the spline should be a quadratic over the last interval;
    !    1: the first derivative at the right endpoint should be YBCEND;
    !    2: the second derivative at the right endpoint should be YBCEND;
    !    3: Not-a-knot: the third derivative is continuous at T(N-1).
    !
    !    Input, real(dp) YBCEND, the right boundary value, if needed.
    !
    !    Output, real(dp) YPP(N), the second derivatives of
    !    the cubic spline.
    !
    implicit none

    integer(i4) n

    real(dp) a1(n)
    real(dp) a2(n)
    real(dp) a3(n)
    real(dp) a4(n)
    real(dp) a5(n)
    real(dp) b(n)
    integer(i4) i
    integer(i4) ibcbeg
    integer(i4) ibcend
!    integer(i4) info
    real(dp) t(n)
    real(dp) y(n)
    real(dp) ybcbeg
    real(dp) ybcend
    real(dp) ypp(n)
    !
    !  Check.
    !
    if ( n <= 1 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
       write ( *, '(a)' ) '  The number of knots must be at least 2.'
       write ( *, '(a,i8)' ) '  The input value of N = ', n
       stop
    end if

    do i = 1, n - 1
       if ( le(t(i+1),t(i)) ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
          write ( *, '(a)' ) '  The knots must be strictly increasing, but'
          write ( *, '(a,i8,a,g14.6)' ) '  T(',  i,') = ', t(i)
          write ( *, '(a,i8,a,g14.6)' ) '  T(',i+1,') = ', t(i+1)
          stop
       end if
    end do
    !
    !  Zero out the matrix.
    !
    a1(1:n) = 0.0_dp
    a2(1:n) = 0.0_dp
    a3(1:n) = 0.0_dp
    a4(1:n) = 0.0_dp
    a5(1:n) = 0.0_dp
    !
    !  Set the first equation.
    !
    if ( ibcbeg == 0 ) then
       b(1) = 0.0_dp
       a3(1) =  1.0_dp
       a4(1) = -1.0_dp
    else if ( ibcbeg == 1 ) then
       b(1) = ( y(2) - y(1) ) / ( t(2) - t(1) ) - ybcbeg
       a3(1) = ( t(2) - t(1) ) / 3.0_dp
       a4(1) = ( t(2) - t(1) ) / 6.0_dp
    else if ( ibcbeg == 2 ) then
       b(1) = ybcbeg
       a3(1) = 1.0_dp
       a4(1) = 0.0_dp
    else if ( ibcbeg == 3 ) then
       b(1) = 0.0_dp
       a3(1) = - ( t(3) - t(2) )
       a4(1) =   ( t(3)        - t(1) )
       a5(1) = - (        t(2) - t(1) )
    else
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
       write ( *, '(a)' ) '  The boundary flag IBCBEG must be 0, 1, 2 or 3.'
       write ( *, '(a,i8)' ) '  The input value is IBCBEG = ', ibcbeg
       stop
    end if
    !
    !  Set the intermediate equations.
    !
    do i = 2, n - 1
       b(i) = ( y(i+1) - y(i) ) / ( t(i+1) - t(i) ) &
            - ( y(i) - y(i-1) ) / ( t(i) - t(i-1) )
       a2(i) = ( t(i+1) - t(i)   ) / 6.0_dp
       a3(i) = ( t(i+1) - t(i-1) ) / 3.0_dp
       a4(i) = ( t(i)   - t(i-1) ) / 6.0_dp
    end do
    !
    !  Set the last equation.
    !
    if ( ibcend == 0 ) then
       b(n) = 0.0_dp
       a2(n) = -1.0_dp
       a3(n) = 1.0_dp
    else if ( ibcend == 1 ) then
       b(n) = ybcend - ( y(n) - y(n-1) ) / ( t(n) - t(n-1) )
       a2(n) = ( t(n) - t(n-1) ) / 6.0_dp
       a3(n) = ( t(n) - t(n-1) ) / 3.0_dp
    else if ( ibcend == 2 ) then
       b(n) = ybcend
       a2(n) = 0.0_dp
       a3(n) = 1.0_dp
    else if ( ibcend == 3 ) then
       b(n) = 0.0_dp
       a1(n) = - ( t(n) - t(n-1) )
       a2(n) =   ( t(n)          - t(n-2) )
       a3(n) = - (        t(n-1) - t(n-2) )
    else
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_CUBIC_SET - Fatal error!'
       write ( *, '(a)' ) '  The boundary flag IBCEND must be 0, 1, 2 or 3.'
       write ( *, '(a,i8)' ) '  The input value is IBCEND = ', ibcend
       stop
    end if
    !
    !  Special case:
    !    N = 2, IBCBEG = IBCEND = 0.
    !
    if ( n == 2 .and. ibcbeg == 0 .and. ibcend == 0 ) then

       ypp(1) = 0.0_dp
       ypp(2) = 0.0_dp
       !
       !  Solve the linear system.
       !
    else

       call penta ( n, a1, a2, a3, a4, a5, b, ypp )

    end if

    return
  end subroutine spline_cubic_set
  subroutine penta ( n, a1, a2, a3, a4, a5, b, x )

    !*****************************************************************************80
    !
    !! PENTA solves a pentadiagonal system of linear equations.
    !
    !  Discussion:
    !
    !    The matrix A is pentadiagonal.  It is entirely zero, except for
    !    the main diagaonal, and the two immediate sub- and super-diagonals.
    !
    !    The entries of Row I are stored as:
    !
    !      A(I,I-2) -> A1(I)
    !      A(I,I-1) -> A2(I)
    !      A(I,I)   -> A3(I)
    !      A(I,I+1) -> A4(I)
    !      A(I,I-2) -> A5(I)
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    07 June 2013
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Cheney, Kincaid,
    !    Numerical Mathematics and Computing,
    !    1985, pages 233-236.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the order of the matrix.
    !
    !    Input, real(dp) A1(N), A2(N), A3(N), A4(N), A5(N), the nonzero
    !    elements of the matrix.  Note that the data in A2, A3 and A4
    !    is overwritten by this routine during the solution process.
    !
    !    Input, real(dp) B(N), the right hand side of the linear system.
    !
    !    Output, real(dp) X(N), the solution of the linear system.
    !
    implicit none

    integer(i4) n

    real(dp) a1(n)
    real(dp) a2(n)
    real(dp) a3(n)
    real(dp) a4(n)
    real(dp) a5(n)
    real(dp) b(n)
    integer(i4) i
    real(dp) x(n)
    real(dp) xmult

    do i = 2, n - 1
       xmult = a2(i) / a3(i-1)
       a3(i) = a3(i) - xmult * a4(i-1)
       a4(i) = a4(i) - xmult * a5(i-1)
       b(i) = b(i) - xmult * b(i-1)
       xmult = a1(i+1) / a3(i-1)
       a2(i+1) = a2(i+1) - xmult * a4(i-1)
       a3(i+1) = a3(i+1) - xmult * a5(i-1)
       b(i+1) = b(i+1) - xmult * b(i-1)
    end do

    xmult = a2(n) / a3(n-1)
    a3(n) = a3(n) - xmult * a4(n-1)
    x(n) = ( b(n) - xmult * b(n-1) ) / a3(n)
    x(n-1) = ( b(n-1) - a4(n-1) * x(n) ) / a3(n-1)
    do i = n - 2, 1, -1
       x(i) = ( b(i) - a4(i) * x(i+1) - a5(i) * x(i+2) ) / a3(i)
    end do

    return
  end subroutine penta
  subroutine spline_cubic_val ( n, t, y, ypp, tval, yval, ypval, yppval )

    !*****************************************************************************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:
    !
    !    20 November 2000
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Carl deBoor,
    !    A Practical Guide to Splines,
    !    Springer, 2001,
    !    ISBN: 0387953663.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of data values.
    !
    !    Input, real(dp) T(N), the knot values.
    !
    !    Input, real(dp) Y(N), the data values at the knots.
    !
    !    Input, real(dp) YPP(N), the second derivatives of the
    !    spline at the knots.
    !
    !    Input, real(dp) 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(dp) YVAL, YPVAL, YPPVAL, the value of the spline, and
    !    its first two derivatives at TVAL.
    !
    implicit none

    integer(i4) n

    real(dp) dt
    real(dp) h
    integer(i4) left
    integer(i4) right
    real(dp) t(n)
    real(dp) tval
    real(dp) y(n)
    real(dp) ypp(n)
    real(dp) yppval
    real(dp) ypval
    real(dp) yval
    !
    !  Determine the interval [T(LEFT), T(RIGHT)] that contains TVAL.
    !  Values below T(1) or above T(N) use extrapolation.
    !
    call r8vec_bracket ( n, t, tval, left, right )
    !
    !  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_dp + ypp(left) / 3.0_dp ) * h &
         + dt * ( 0.5_dp * ypp(left) &
         + dt * ( ( ypp(right) - ypp(left) ) / ( 6.0_dp * h ) ) ) )

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

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

    return
  end subroutine spline_cubic_val
  subroutine spline_cubic_val2 ( n, t, y, ypp, left, tval, yval, ypval, yppval )

    !*****************************************************************************80
    !
    !! SPLINE_CUBIC_VAL2 evaluates a piecewise cubic spline at a point.
    !
    !  Discussion:
    !
    !    This routine is a modification of SPLINE_CUBIC_VAL; it allows the
    !    user to speed up the code by suggesting the appropriate T interval
    !    to search first.
    !
    !    SPLINE_CUBIC_SET must have already been called to define the
    !    values of YPP.
    !
    !    In the LEFT interval, let RIGHT = LEFT+1.  The form of the spline is
    !
    !      SPL(T) =
    !          A
    !        + B * ( T - T(LEFT) )
    !        + C * ( T - T(LEFT) )^2
    !        + D * ( T - T(LEFT) )^3
    !
    !    Here:
    !      A = Y(LEFT)
    !      B = ( Y(RIGHT) - Y(LEFT) ) / ( T(RIGHT) - T(LEFT) )
    !        - ( YPP(RIGHT) + 2 * YPP(LEFT) ) * ( T(RIGHT) - T(LEFT) ) / 6
    !      C = YPP(LEFT) / 2
    !      D = ( YPP(RIGHT) - YPP(LEFT) ) / ( 6 * ( T(RIGHT) - T(LEFT) ) )
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    24 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Carl deBoor,
    !    A Practical Guide to Splines,
    !    Springer, 2001,
    !    ISBN: 0387953663.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of knots.
    !
    !    Input, real(dp) T(N), the knot values.
    !
    !    Input, real(dp) Y(N), the data values at the knots.
    !
    !    Input, real(dp) YPP(N), the second derivatives of the spline at
    !    the knots.
    !
    !    Input/output, integer(i4) LEFT, the suggested T interval to 
    !    search.  LEFT should be between 1 and N-1.  If LEFT is not in this range,
    !    then its value will be ignored.  On output, LEFT is set to the
    !    actual interval in which TVAL lies.
    !
    !    Input, real(dp) 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(dp) YVAL, YPVAL, YPPVAL, the value of the spline, and
    !    its first two derivatives at TVAL.
    !
    implicit none

    integer(i4) n

    real(dp) dt
    real(dp) h
    integer(i4) left
    integer(i4) right
    real(dp) t(n)
    real(dp) tval
    real(dp) y(n)
    real(dp) ypp(n)
    real(dp) yppval
    real(dp) ypval
    real(dp) yval
    !
    !  Determine the interval [T(LEFT), T(RIGHT)] that contains TVAL.
    !
    !  What you want from R8VEC_BRACKET3 is that TVAL is to be computed
    !  by the data in interval [T(LEFT), T(RIGHT)].
    !
    call r8vec_bracket3 ( n, t, tval, left )
    right = left + 1
    !
    !  In the interval LEFT, the polynomial is in terms of a normalized
    !  coordinate  ( DT / H ) between 0 and 1.
    !
    dt = tval - t(left)
    h = t(right) - t(left)

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

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

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

    return
  end subroutine spline_cubic_val2
  subroutine spline_hermite_set ( ndata, tdata, ydata, ypdata, c )

    !*****************************************************************************80
    !
    !! SPLINE_HERMITE_SET sets up a piecewise cubic Hermite interpolant.
    !
    !  Discussion:
    !
    !    Once the array C is computed, then in the interval
    !    (TDATA(I), TDATA(I+1)), the interpolating Hermite polynomial
    !    is given by
    !
    !      SVAL(TVAL) =                 C(1,I)
    !         + ( TVAL - TDATA(I) ) * ( C(2,I)
    !         + ( TVAL - TDATA(I) ) * ( C(3,I)
    !         + ( TVAL - TDATA(I) ) *   C(4,I) ) )
    !
    !    This is algorithm CALCCF of Conte and deBoor.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    11 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Samuel Conte, Carl deBoor,
    !    Elementary Numerical Analysis,
    !    Second Edition,
    !    McGraw Hill, 1972,
    !    ISBN: 07-012446-4,
    !    LC: QA297.C65.
    !
    !  Parameters:
    !
    !    Input, integer(i4) NDATA, the number of data points.
    !    NDATA must be at least 2.
    !
    !    Input, real(dp) TDATA(NDATA), the abscissas of the data points.
    !    The entries of TDATA are assumed to be strictly increasing.
    !
    !    Input, real(dp) Y(NDATA), YP(NDATA), the value of the
    !    function and its derivative at TDATA(1:NDATA).
    !
    !    Output, real(dp) C(4,NDATA), the coefficients of the
    !    Hermite polynomial.
    !    C(1,1:NDATA) = Y(1:NDATA) and C(2,1:NDATA) = YP(1:NDATA).
    !    C(3,1:NDATA-1) and C(4,1:NDATA-1) are the quadratic and cubic
    !    coefficients.
    !
    implicit none

    integer(i4) ndata

    real(dp) c(4,ndata)
    real(dp) divdif1
    real(dp) divdif3
    real(dp) dt
    integer(i4) i
    real(dp) tdata(ndata)
    real(dp) ydata(ndata)
    real(dp) ypdata(ndata)
    !
    !  Check NDATA.
    !
    if ( ndata < 2 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_HERMITE_SET - Fatal error!'
       write ( *, '(a)' ) '  NDATA < 2.'
       stop
    end if

    c(1,1:ndata) = ydata(1:ndata)
    c(2,1:ndata) = ypdata(1:ndata)

    do i = 1, ndata - 1
       dt = tdata(i+1) - tdata(i)
       divdif1 = ( c(1,i+1) - c(1,i) ) / dt
       divdif3 = c(2,i) + c(2,i+1) - 2.0_dp * divdif1
       c(3,i) = ( divdif1 - c(2,i) - divdif3 ) / dt
       c(4,i) = divdif3 / ( dt * dt )
    end do

    c(3,ndata) = 0.0_dp
    c(4,ndata) = 0.0_dp

    return
  end subroutine spline_hermite_set
  subroutine spline_hermite_val ( ndata, tdata, c, tval, sval, spval )

    !*****************************************************************************80
    !
    !! 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.
    !
    !    In the interval (TDATA(I), TDATA(I+1)), the interpolating
    !    Hermite polynomial is given by
    !
    !      SVAL(TVAL) =                 C(1,I)
    !         + ( TVAL - TDATA(I) ) * ( C(2,I)
    !         + ( TVAL - TDATA(I) ) * ( C(3,I)
    !         + ( TVAL - TDATA(I) ) *   C(4,I) ) )
    !
    !    and
    !
    !      SVAL'(TVAL) =                    C(2,I)
    !         + ( TVAL - TDATA(I) ) * ( 2 * C(3,I)
    !         + ( TVAL - TDATA(I) ) *   3 * C(4,I) )
    !
    !    This is algorithm PCUBIC of Conte and deBoor.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    11 February 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Samuel Conte, Carl deBoor,
    !    Elementary Numerical Analysis,
    !    Second Edition,
    !    McGraw Hill, 1972,
    !    ISBN: 07-012446-4,
    !    LC: QA297.C65.
    !
    !  Parameters:
    !
    !    Input, integer(i4) NDATA, the number of data points.
    !    NDATA must be at least 2.
    !
    !    Input, real(dp) TDATA(NDATA), the abscissas of the data points.
    !    The entries of TDATA are assumed to be strictly increasing.
    !
    !    Input, real(dp) C(4,NDATA), the coefficient data computed by
    !    SPLINE_HERMITE_SET.
    !
    !    Input, real(dp) TVAL, the point where the interpolant is to
    !    be evaluated.
    !
    !    Output, real(dp) SVAL, SPVAL, the value of the interpolant
    !    and its derivative at TVAL.
    !
    implicit none

    integer(i4) ndata

    real(dp) c(4,ndata)
    real(dp) dt
    integer(i4) left
    integer(i4) right
    real(dp) spval
    real(dp) sval
    real(dp) tdata(ndata)
    real(dp) tval
    !
    !  Check NDATA.
    !
    if ( ndata < 2 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_HERMITE_VAL - Fatal error!'
       write ( *, '(a)' ) '  NDATA < 2.'
       stop
    end if
    !
    !  Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains
    !  or is nearest to TVAL.
    !
    call r8vec_bracket ( ndata, tdata, tval, left, right )
    !
    !  Evaluate the cubic polynomial.
    !
    dt = tval - tdata(left)

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

    spval = c(2,left) + dt * ( 2.0_dp * c(3,left) + dt * 3.0_dp * c(4,left) )

    return
  end subroutine spline_hermite_val
  subroutine spline_linear_int ( ndata, tdata, ydata, a, b, int_val )

    !*****************************************************************************80
    !
    !! SPLINE_LINEAR_INT evaluates the integral of a piecewise linear spline.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    01 November 2001
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) NDATA, the number of data points defining
    !    the spline.  NDATA must be at least 2.
    !
    !    Input, real(dp) TDATA(NDATA), YDATA(NDATA), the values of
    !    the independent and dependent variables at the data points.  The
    !    values of TDATA should be distinct and increasing.
    !
    !    Input, real(dp) A, B, the interval over which the
    !    integral is desired.
    !
    !    Output, real(dp) INT_VAL, the value of the integral.
    !
    implicit none

    integer(i4) ndata

    real(dp) a
    real(dp) a_copy
    integer(i4) a_left
    integer(i4) a_right
    real(dp) b
    real(dp) b_copy
    integer(i4) b_left
    integer(i4) b_right
    integer(i4) i_left
    real(dp) int_val
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) ydata(ndata)
    real(dp) yp
    real(dp) yval
    !
    !  Check NDATA.
    !
    if ( ndata < 2 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_LINEAR_INT - Fatal error!'
       write ( *, '(a)' ) '  NDATA < 2.'
       stop
    end if

    int_val = 0.0_dp

    if ( eq(a,b) ) then
       return
    end if

    a_copy = min ( a, b )
    b_copy = max ( a, b )
    !
    !  Find the interval [ TDATA(A_LEFT), TDATA(A_RIGHT) ] that contains, or is
    !  nearest to, A.
    !
    call r8vec_bracket ( ndata, tdata, a_copy, a_left, a_right )
    !
    !  Find the interval [ TDATA(B_LEFT), TDATA(B_RIGHT) ] that contains, or is
    !  nearest to, B.
    !
    call r8vec_bracket ( ndata, tdata, b_copy, b_left, b_right )
    !
    !  If A and B are in the same interval...
    !
    if ( a_left == b_left ) then

       tval = ( a_copy + b_copy ) / 2.0_dp

       yp = ( ydata(a_right) - ydata(a_left) ) / &
            ( tdata(a_right) - tdata(a_left) )

       yval = ydata(a_left) + ( tval - tdata(a_left) ) * yp

       int_val = yval * ( b_copy - a_copy )

       return
    end if
    !
    !  Otherwise, integrate from:
    !
    !  A               to TDATA(A_RIGHT),
    !  TDATA(A_RIGHT)  to TDATA(A_RIGHT+1),...
    !  TDATA(B_LEFT-1) to TDATA(B_LEFT),
    !  TDATA(B_LEFT)   to B.
    !
    !  Use the fact that the integral of a linear function is the
    !  value of the function at the midpoint times the width of the interval.
    !
    tval = ( a_copy + tdata(a_right) ) / 2.0_dp

    yp = ( ydata(a_right) - ydata(a_left) ) / &
         ( tdata(a_right) - tdata(a_left) )

    yval = ydata(a_left) + ( tval - tdata(a_left) ) * yp

    int_val = int_val + yval * ( tdata(a_right) - a_copy )

    do i_left = a_right, b_left - 1

       tval = ( tdata(i_left+1) + tdata(i_left) ) / 2.0_dp

       yp = ( ydata(i_left+1) - ydata(i_left) ) / &
            ( tdata(i_left+1) - tdata(i_left) )

       yval = ydata(i_left) + ( tval - tdata(i_left) ) * yp

       int_val = int_val + yval * ( tdata(i_left + 1) - tdata(i_left) )

    end do

    tval = ( tdata(b_left) + b_copy ) / 2.0_dp

    yp = ( ydata(b_right) - ydata(b_left) ) / &
         ( tdata(b_right) - tdata(b_left) )

    yval = ydata(b_left) + ( tval - tdata(b_left) ) * yp

    int_val = int_val + yval * ( b_copy - tdata(b_left) )

    if ( b < a ) then
       int_val = - int_val
    end if

    return
  end subroutine spline_linear_int
  subroutine spline_linear_intset ( n, int_x, int_v, data_x, data_y )

    !*****************************************************************************80
    !
    !! SPLINE_LINEAR_INTSET: piecewise linear spline with given integral properties.
    !
    !  Discussion:
    !
    !    The user has in mind an interval, divided by INT_N+1 points into
    !    INT_N intervals.  A linear spline is to be constructed,
    !    with breakpoints at the centers of each interval, and extending
    !    continuously to the left of the first and right of the last
    !    breakpoints.  The constraint on the linear spline is that it is
    !    required that it have a given integral value over each interval.
    !
    !    A tridiagonal linear system of equations is solved for the
    !    values of the spline at the breakpoints.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    27 January 2004
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of intervals.
    !
    !    Input, real(dp) INT_X(N+1), the points that define the intervals.
    !    Interval I lies between INT_X(I) and INT_X(I+1).
    !
    !    Input, real(dp) INT_V(N), the desired value of the integral of the
    !    linear spline over each interval.
    !
    !    Output, real(dp) DATA_X(N), DATA_Y(N), the values of the
    !    independent and dependent variables at the data points.  The values
    !    of DATA_X are the interval midpoints.  The values of DATA_Y are
    !    determined in such a way that the exact integral of the linear
    !    spline over interval I is equal to INT_V(I).
    !
    implicit none

    integer(i4) n

    real(dp) a(3,n)
    real(dp) data_x(n)
    real(dp) data_y(n)
    real(dp) int_v(n)
    real(dp) int_x(n+1)
    !
    !  Set up the easy stuff.
    !
    data_x(1:n) = 0.5_dp * ( int_x(1:n) + int_x(2:n+1) )
    !
    !  Set up C, D, E, the coefficients of the linear system.
    !
    a(3,1:n-2) = 1.0_dp &
         - ( 0.5_dp * ( data_x(2:n-1) + int_x(2:n-1) ) - data_x(1:n-2) ) &
         / ( data_x(2:n-1) - data_x(1:n-2) )
    a(3,n-1) = 0.0_dp
    a(3,n) = 0.0_dp

    a(2,1) = int_x(2) - int_x(1)

    a(2,2:n-1) = 1.0_dp &
         + ( 0.5_dp * ( data_x(2:n-1) + int_x(2:n-1) ) &
         - data_x(1:n-2) ) &
         / ( data_x(2:n-1) - data_x(1:n-2) ) &
         - ( 0.5_dp * ( data_x(2:n-1) + int_x(3:n) ) - data_x(2:n-1) ) &
         / ( data_x(3:n) - data_x(2:n-1) )

    a(2,n) = int_x(n+1) - int_x(n)

    a(1,1) = 0.0_dp
    a(1,2) = 0.0_dp

    a(1,3:n) = ( 0.5_dp * ( data_x(2:n-1) + int_x(3:n) ) &
         - data_x(2:n-1) ) / ( data_x(3:n) - data_x(2:n-1) )
    !
    !  Set up DATA_Y, which begins as the right hand side of the linear system.
    !
    data_y(1) = int_v(1)
    data_y(2:n-1) = 2.0_dp * int_v(2:n-1) / ( int_x(3:n) - int_x(2:n-1) )
    data_y(n) = int_v(n)
    !
    !  Solve the linear system.
    !
    call r83_np_fs ( n, a, data_y, data_y )

    return
  end subroutine spline_linear_intset
  subroutine spline_linear_val ( ndata, tdata, ydata, tval, yval, ypval )

    !*****************************************************************************80
    !
    !! SPLINE_LINEAR_VAL evaluates a piecewise linear spline at a point.
    !
    !  Discussion:
    !
    !    Because of the extremely simple form of the linear spline,
    !    the raw data points ( TDATA(I), YDATA(I)) can be used directly to
    !    evaluate the spline at any point.  No processing of the data
    !    is required.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    06 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) NDATA, the number of data points defining
    !    the spline.  NDATA must be at least 2.
    !
    !    Input, real(dp) TDATA(NDATA), YDATA(NDATA), the values of
    !    the independent and dependent variables at the data points.  The
    !    values of TDATA should be distinct and increasing.
    !
    !    Input, real(dp) TVAL, the point at which the spline is
    !    to be evaluated.
    !
    !    Output, real(dp) YVAL, YPVAL, the value of the spline and
    !    its first derivative dYdT at TVAL.  YPVAL is not reliable if TVAL
    !    is exactly equal to TDATA(I) for some I.
    !
    implicit none

    integer(i4) ndata

    integer(i4) left
    integer(i4) right
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) ydata(ndata)
    real(dp) ypval
    real(dp) yval
    !
    !  Check NDATA.
    !
    if ( ndata < 2 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_LINEAR_VAL - Fatal error!'
       write ( *, '(a)' ) '  NDATA < 2.'
       stop
    end if
    !
    !  Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains, or is
    !  nearest to, TVAL.
    !
    call r8vec_bracket ( ndata, tdata, tval, left, right )
    !
    !  Now evaluate the piecewise linear function.
    !
    ypval = ( ydata(right) - ydata(left) ) / ( tdata(right) - tdata(left) )

    yval = ydata(left) +  ( tval - tdata(left) ) * ypval

    return
  end subroutine spline_linear_val
  subroutine spline_overhauser_nonuni_val ( ndata, tdata, ydata, tval, yval )

    !*****************************************************************************80
    !
    !! SPLINE_OVERHAUSER_NONUNI_VAL evaluates the nonuniform Overhauser spline.
    !
    !  Discussion:
    !
    !    The nonuniformity refers to the fact that the abscissa values
    !    need not be uniformly spaced.
    !
    !    Thanks to Doug Fortune for pointing out that the point distances
    !    used to define ALPHA and BETA should be the Euclidean distances
    !    between the points.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    08 May 2007
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) NDATA, the number of data points.
    !    3 <= NDATA is required.
    !
    !    Input, real(dp) TDATA(NDATA), the abscissas of the data points.
    !    The values of TDATA are assumed to be distinct and increasing.
    !
    !    Input, real(dp) YDATA(NDATA), the data values.
    !
    !    Input, real(dp) TVAL, the value where the spline is to
    !    be evaluated.
    !
    !    Output, real(dp) YVAL, the value of the spline at TVAL.
    !
    implicit none

    integer(i4) ndata

    real(dp) alpha
    real(dp) beta
    real(dp) d21
    real(dp) d32
    real(dp) d43
    integer(i4) left
    real(dp) mbasis(4,4)
    real(dp) mbasis_l(3,3)
    real(dp) mbasis_r(3,3)
    integer(i4) right
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) ydata(ndata)
    real(dp) yval
    !
    !  Check NDATA.
    !
    if ( ndata < 3 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_OVERHAUSER_NONUNI_VAL - Fatal error!'
       write ( *, '(a)' ) '  NDATA < 3.'
       stop
    end if
    !
    !  Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
    !
    call r8vec_bracket ( ndata, tdata, tval, left, right )
    !
    !  Evaluate the spline in the given interval.
    !
    if ( left == 1 ) then

       d21 = sqrt ( ( tdata(2) - tdata(1) )**2 &
            + ( ydata(2) - ydata(1) )**2 )

       d32 = sqrt ( ( tdata(3) - tdata(2) )**2 &
            + ( ydata(3) - ydata(2) )**2 )

       alpha = d21 / ( d32 + d21 )

       call basis_matrix_overhauser_nul ( alpha, mbasis_l )

       call basis_matrix_tmp ( left, 3, mbasis_l, ndata, tdata, ydata, tval, yval )

    else if ( left < ndata - 1 ) then

       d21 = sqrt ( ( tdata(left) - tdata(left-1) )**2 &
            + ( ydata(left) - ydata(left-1) )**2 )

       d32 = sqrt ( ( tdata(left+1) - tdata(left) )**2 &
            + ( ydata(left+1) - ydata(left) )**2 )

       d43 = sqrt ( ( tdata(left+2) - tdata(left+1) )**2 &
            + ( ydata(left+2) - ydata(left+1) )**2 )

       alpha = d21 / ( d32 + d21 )
       beta  = d32 / ( d43 + d32 )

       call basis_matrix_overhauser_nonuni ( alpha, beta, mbasis )

       call basis_matrix_tmp ( left, 4, mbasis, ndata, tdata, ydata, tval, yval )

    else if ( left == ndata - 1 ) then

       d32 = sqrt ( ( tdata(ndata-1) - tdata(ndata-2) )**2 &
            + ( ydata(ndata-1) - ydata(ndata-2) )**2 )

       d43 = sqrt ( ( tdata(ndata) - tdata(ndata-1) )**2 &
            + ( ydata(ndata) - ydata(ndata-1) )**2 )

       beta  = d32 / ( d43 + d32 )

       call basis_matrix_overhauser_nur ( beta, mbasis_r )

       call basis_matrix_tmp ( left, 3, mbasis_r, ndata, tdata, ydata, tval, yval )

    else

       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_OVERHAUSER_NONUNI_VAL - Fatal error!'
       write ( *, '(a,i8)' ) '  Nonsensical value of LEFT = ', left
       write ( *, '(a,i8)' ) '  but 0 < LEFT < NDATA = ', ndata
       write ( *, '(a)' ) '  is required.'
       stop

    end if

    return
  end subroutine spline_overhauser_nonuni_val
  subroutine spline_overhauser_uni_val ( ndata, tdata, ydata, tval, yval )

    !*****************************************************************************80
    !
    !! SPLINE_OVERHAUSER_UNI_VAL evaluates the uniform Overhauser spline.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    06 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) NDATA, the number of data points.
    !    NDATA must be at least 3.
    !
    !    Input, real(dp) TDATA(NDATA), the abscissas of the data points.
    !    The values of TDATA are assumed to be distinct and increasing.
    !    This routine also assumes that the values of TDATA are uniformly
    !    spaced; for instance, TDATA(1) = 10, TDATA(2) = 11, TDATA(3) = 12...
    !
    !    Input, real(dp) YDATA(NDATA), the data values.
    !
    !    Input, real(dp) TVAL, the value where the spline is to
    !    be evaluated.
    !
    !    Output, real(dp) YVAL, the value of the spline at TVAL.
    !
    implicit none

    integer(i4) ndata

    integer(i4) left
    real(dp) mbasis(4,4)
    real(dp) mbasis_l(3,3)
    real(dp) mbasis_r(3,3)
    integer(i4) right
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) ydata(ndata)
    real(dp) yval
    !
    !  Check NDATA.
    !
    if ( ndata < 3 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_OVERHAUSER_UNI_VAL - Fatal error!'
       write ( *, '(a)' ) '  NDATA < 3.'
       stop
    end if
    !
    !  Find the nearest interval [ TDATA(LEFT), TDATA(RIGHT) ] to TVAL.
    !
    call r8vec_bracket ( ndata, tdata, tval, left, right )
    !
    !  Evaluate the spline in the given interval.
    !
    if ( left == 1 ) then

       call basis_matrix_overhauser_uni_l ( mbasis_l )

       call basis_matrix_tmp ( left, 3, mbasis_l, ndata, tdata, ydata, tval, yval )

    else if ( left < ndata - 1 ) then

       call basis_matrix_overhauser_uni ( mbasis )

       call basis_matrix_tmp ( left, 4, mbasis, ndata, tdata, ydata, tval, yval )

    else if ( left == ndata - 1 ) then

       call basis_matrix_overhauser_uni_r ( mbasis_r )

       call basis_matrix_tmp ( left, 3, mbasis_r, ndata, tdata, ydata, tval, yval )

    end if

    return
  end subroutine spline_overhauser_uni_val
  subroutine spline_overhauser_val ( dim_num, ndata, tdata, ydata, tval, yval )

    !*****************************************************************************80
    !
    !! SPLINE_OVERHAUSER_VAL evaluates an Overhauser spline.
    !
    !  Discussion:
    !
    !    Over the first and last intervals, the Overhauser spline is a
    !    quadratic.  In the intermediate intervals, it is a piecewise cubic.
    !    The Overhauser spline is also known as the Catmull-Rom spline.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !   08 April 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    JA Brewer, DC Anderson,
    !    Visual Interaction with Overhauser Curves and Surfaces,
    !    SIGGRAPH 77,
    !    in Proceedings of the 4th Annual Conference on Computer Graphics
    !    and Interactive Techniques,
    !    ASME, July 1977, pages 132-137.
    !
    !    Edwin Catmull, Raphael Rom,
    !    A Class of Local Interpolating Splines,
    !    in Computer Aided Geometric Design,
    !    edited by Robert Barnhill, Richard Reisenfeld,
    !    Academic Press, 1974, pages 317-326,
    !    ISBN: 0120790505.
    !
    !    David Rogers, Alan Adams,
    !    Mathematical Elements of Computer Graphics,
    !    Second Edition,
    !    McGraw Hill, 1989,
    !    ISBN: 0070535299.
    !
    !  Parameters:
    !
    !    Input, integer(i4) DIM_NUM, the dimension of a single data point.
    !    DIM_NUM must be at least 1.
    !
    !    Input, integer(i4) NDATA, the number of data points.
    !    NDATA must be at least 3.
    !
    !    Input, real(dp) TDATA(NDATA), the abscissas of the data
    !    points.  The values in TDATA must be in strictly ascending order.
    !
    !    Input, real(dp) YDATA(DIM_NUM,NDATA), the data points
    !    corresponding to the abscissas.
    !
    !    Input, real(dp) TVAL, the abscissa value at which the spline
    !    is to be evaluated.  Normally, TDATA(1) <= TVAL <= T(NDATA), and
    !    the data will be interpolated.  For TVAL outside this range,
    !    extrapolation will be used.
    !
    !    Output, real(dp) YVAL(DIM_NUM), the value of the spline at TVAL.
    !
    implicit none

    integer(i4) ndata
    integer(i4) dim_num

    integer(i4) left
    integer(i4) order
    integer(i4) right
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) ydata(dim_num,ndata)
    real(dp) yl(dim_num)
    real(dp) yr(dim_num)
    real(dp) yval(dim_num)
    !
    !  Check.
    !
    call r8vec_order_type ( ndata, tdata, order )

    if ( order /= 2 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_OVERHAUSER_VAL - Fatal error!'
       write ( *, '(a)' ) '  The data abscissas are not strictly ascending.'
       stop
    end if

    if ( ndata < 3 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_OVERHAUSER_VAL - Fatal error!'
       write ( *, '(a)' ) '  NDATA < 3.'
       stop
    end if

    if ( dim_num < 1 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_OVERHAUSER_VAL - Fatal error!'
       write ( *, '(a)' ) '  DIM_NUM < 1.'
       stop
    end if
    !
    !  Locate the abscissa interval T(LEFT), T(LEFT+1) nearest to or
    !  containing TVAL.
    !
    call r8vec_bracket ( ndata, tdata, tval, left, right )
    !
    !  Evaluate the "left hand" quadratic defined at T(LEFT-1), T(LEFT), T(RIGHT).
    !
    if ( 0 < left-1 ) then
       call parabola_val2 ( dim_num, ndata, tdata, ydata, left-1, tval, yl )
    end if
    !
    !  Evaluate the "right hand" quadratic defined at T(LEFT), T(RIGHT), T(RIGHT+1).
    !
    if ( right+1 <= ndata ) then
       call parabola_val2 ( dim_num, ndata, tdata, ydata, left, tval, yr )
    end if
    !
    !  Average the quadratics.
    !
    if ( left == 1 ) then

       yval(1:dim_num) = yr(1:dim_num)

    else if ( right < ndata ) then

       yval(1:dim_num) =  &
            ( ( tdata(right) - tval               ) * yl(1:dim_num)   &
            + (                tval - tdata(left) ) * yr(1:dim_num) ) &
            / ( tdata(right)        - tdata(left) )

    else

       yval(1:dim_num) = yl(1:dim_num)

    end if

    return
  end subroutine spline_overhauser_val
  subroutine spline_pchip_set ( n, x, f, d )

    !*****************************************************************************80
    !
    !! SPLINE_PCHIP_SET sets derivatives for a piecewise cubic Hermite interpolant.
    !
    !  Discussion:
    !
    !    This routine computes what would normally be called a Hermite
    !    interpolant.  However, the user is only required to supply function
    !    values, not derivative values as well.  This routine computes
    !    "suitable" derivative values, so that the resulting Hermite interpolant
    !    has desirable shape and monotonicity properties.
    !
    !    The interpolant will have an extremum at each point where
    !    monotonicity switches direction.
    !
    !    The resulting piecewise cubic Hermite function may be evaluated
    !    by SPLINE_PCHIP_VAL.
    !
    !    This routine was originally named "PCHIM".
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    15 December 2008
    !
    !  Author:
    !
    !    Original FORTRAN77 version by Fred Fritsch.
    !    FORTRAN90 version by John Burkardt.
    !
    !  Reference:
    !
    !    Fred Fritsch, Ralph Carlson,
    !    Monotone Piecewise Cubic Interpolation,
    !    SIAM Journal on Numerical Analysis,
    !    Volume 17, Number 2, April 1980, pages 238-246.
    !
    !    Fred Fritsch, Judy Butland,
    !    A Method for Constructing Local Monotone Piecewise Cubic Interpolants,
    !    SIAM Journal on Scientific and Statistical Computing,
    !    Volume 5, Number 2, 1984, pages 300-304.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of data points.  N must be
    !    at least 2.
    !
    !    Input, real(dp) X(N), the strictly increasing independent
    !    variable values.
    !
    !    Input, real(dp) F(N), dependent variable values to be
    !    interpolated.  F(I) is the value corresponding to X(I).
    !    This routine is designed for monotonic data, but it will work for any
    !    F array.  It will force extrema at points where monotonicity switches
    !    direction.
    !
    !    Output, real(dp) D(N), the derivative values at the
    !    data points.  If the data are monotonic, these values will determine
    !    a monotone cubic Hermite function.
    !
    implicit none

    integer(i4) n

    real(dp) d(n)
    real(dp) del1
    real(dp) del2
    real(dp) dmax
    real(dp) dmin
    real(dp) drat1
    real(dp) drat2
    real(dp) dsave
    real(dp) f(n)
    real(dp) h1
    real(dp) h2
    real(dp) hsum
    real(dp) hsumt3
    integer(i4) i
    integer(i4) ierr
    integer(i4) nless1
    !real(dp) pchst
    real(dp) temp
    real(dp) w1
    real(dp) w2
    real(dp) x(n)
    !
    !  Check the arguments.
    !
    if ( n < 2 ) then
       ierr = -1
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_PCHIP_SET - Fatal error!'
       write ( *, '(a)' ) '  Number of data points less than 2.'
       stop
    end if

    do i = 2, n
       if ( le(x(i),x(i-1)) ) then
          ierr = -3
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'SPLINE_PCHIP_SET - Fatal error!'
          write ( *, '(a)' ) '  X array not strictly increasing.'
          stop
       end if
    end do

    ierr = 0
    nless1 = n - 1
    h1 = x(2) - x(1)
    del1 = ( f(2) - f(1) ) / h1
    dsave = del1
    !
    !  Special case N=2, use linear interpolation.
    !
    if ( n == 2 ) then
       d(1) = del1
       d(n) = del1
       return
    end if
    !
    !  Normal case, 3 <= N.
    !
    h2 = x(3) - x(2)
    del2 = ( f(3) - f(2) ) / h2
    !
    !  Set D(1) via non-centered three point formula, adjusted to be
    !  shape preserving.
    !
    hsum = h1 + h2
    w1 = ( h1 + hsum ) / hsum
    w2 = -h1 / hsum
    d(1) = w1 * del1 + w2 * del2

    if ( le(pchst( d(1), del1 ), 0.0_dp) ) then

       d(1) = 0.0_dp
       !
       !  Need do this check only if monotonicity switches.
       !
    else if ( pchst ( del1, del2 ) < 0.0_dp ) then

       dmax = 3.0_dp * del1

       if ( abs ( dmax ) < abs ( d(1) ) ) then
          d(1) = dmax
       end if

    end if
    !
    !  Loop through interior points.
    !
    do i = 2, nless1

       if ( 2 < i ) then
          h1 = h2
          h2 = x(i+1) - x(i)
          hsum = h1 + h2
          del1 = del2
          del2 = ( f(i+1) - f(i) ) / h2
       end if
       !
       !  Set D(I)=0 unless data are strictly monotonic.
       !
       d(i) = 0.0_dp

       temp = pchst ( del1, del2 )

       if ( temp < 0.0_dp ) then

          ierr = ierr + 1
          dsave = del2
          !
          !  Count number of changes in direction of monotonicity.
          !
       else if ( eq(temp,0.0_dp) ) then

          if ( ne(del2,0.0_dp) ) then
             if ( pchst ( dsave, del2 ) < 0.0_dp ) then
                ierr = ierr + 1
             end if
             dsave = del2
          end if
          !
          !  Use Brodlie modification of Butland formula.
          !
       else

          hsumt3 = 3.0_dp * hsum
          w1 = ( hsum + h1 ) / hsumt3
          w2 = ( hsum + h2 ) / hsumt3
          dmax = max ( abs ( del1 ), abs ( del2 ) )
          dmin = min ( abs ( del1 ), abs ( del2 ) )
          drat1 = del1 / dmax
          drat2 = del2 / dmax
          d(i) = dmin / ( w1 * drat1 + w2 * drat2 )

       end if

    end do
    !
    !  Set D(N) via non-centered three point formula, adjusted to be
    !  shape preserving.
    !
    w1 = -h2 / hsum
    w2 = ( h2 + hsum ) / hsum
    d(n) = w1 * del1 + w2 * del2

    if ( le(pchst( d(n), del2 ), 0.0_dp) ) then
       d(n) = 0.0_dp
    else if ( pchst ( del1, del2 ) < 0.0_dp ) then
       !
       !  Need do this check only if monotonicity switches.
       !
       dmax = 3.0_dp * del2

       if ( abs ( dmax ) < abs ( d(n) ) ) then
          d(n) = dmax
       end if

    end if

    return
  end subroutine spline_pchip_set
  subroutine spline_pchip_val ( n, x, f, d, ne, xe, fe )

    !*****************************************************************************80
    !
    !! SPLINE_PCHIP_VAL evaluates a piecewise cubic Hermite function.
    !
    !  Description:
    !
    !    This routine may be used by itself for Hermite interpolation, or as an
    !    evaluator for SPLINE_PCHIP_SET.
    !
    !    This routine evaluates the cubic Hermite function at the points XE.
    !
    !    Most of the coding between the call to CHFEV and the end of
    !    the IR loop could be eliminated if it were permissible to
    !    assume that XE is ordered relative to X.
    !
    !    CHFEV does not assume that X1 is less than X2.  Thus, it would
    !    be possible to write a version of SPLINE_PCHIP_VAL that assumes a strictly
    !    decreasing X array by simply running the IR loop backwards
    !    and reversing the order of appropriate tests.
    !
    !    The present code has a minor bug, which I have decided is not
    !    worth the effort that would be required to fix it.
    !    If XE contains points in [X(N-1),X(N)], followed by points less than
    !    X(N-1), followed by points greater than X(N), the extrapolation points
    !    will be counted (at least) twice in the total returned in IERR.
    !
    !    The evaluation will be most efficient if the elements of XE are
    !    increasing relative to X; that is, for all J <= K,
    !      X(I) <= XE(J)
    !    implies
    !      X(I) <= XE(K).
    !
    !    If any of the XE are outside the interval [X(1),X(N)],
    !    values are extrapolated from the nearest extreme cubic,
    !    and a warning error is returned.
    !
    !    This routine was originally called "PCHFE".
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    14 August 2005
    !
    !  Author:
    !
    !    Original FORTRAN77 version by Fred Fritsch.
    !    FORTRAN90 version by John Burkardt.
    !
    !  Reference:
    !
    !    Fred Fritsch, Ralph Carlson,
    !    Monotone Piecewise Cubic Interpolation,
    !    SIAM Journal on Numerical Analysis,
    !    Volume 17, Number 2, April 1980, pages 238-246.
    !
    !  Parameters:
    !
    !    Input, integer(i4) N, the number of data points.  N must be
    !    at least 2.
    !
    !    Input, real(dp) X(N), the strictly increasing independent
    !    variable values.
    !
    !    Input, real(dp) F(N), the function values.
    !
    !    Input, real(dp) D(N), the derivative values.
    !
    !    Input, integer(i4) NE, the number of evaluation points.
    !
    !    Input, real(dp) XE(NE), points at which the function is to
    !    be evaluated.
    !
    !    Output, real(dp) FE(NE), the values of the cubic Hermite
    !    function at XE.
    !
    implicit none

    integer(i4) n
    integer(i4) ne

    real(dp) d(n)
    real(dp) f(n)
    real(dp) fe(ne)
    integer(i4) i
    integer(i4) ierc
    integer(i4) ierr
    integer(i4) ir
    integer(i4) j
    integer(i4) j_first
    integer(i4) j_new
    integer(i4) j_save
    integer(i4) next(2)
    integer(i4) nj
    real(dp) x(n)
    real(dp) xe(ne)
    !
    !  Check arguments.
    !
    if ( n < 2 ) then
       ierr = -1
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!'
       write ( *, '(a)' ) '  Number of data points less than 2.'
       stop
    end if

    do i = 2, n
       if ( le(x(i),x(i-1)) ) then
          ierr = -3
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!'
          write ( *, '(a)' ) '  X array not strictly increasing.'
          stop
       end if
    end do

    if ( ne < 1 ) then
       ierr = -4
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!'
       write ( *, '(a)' ) '  Number of evaluation points less than 1.'
       return
    end if

    ierr = 0
    !
    !  Loop over intervals.
    !  The interval index is IL = IR-1.
    !  The interval is X(IL) <= X < X(IR).
    !
    j_first = 1
    ir = 2

    do
       !
       !  Skip out of the loop if have processed all evaluation points.
       !
       if ( ne < j_first ) then
          exit
       end if
       !
       !  Locate all points in the interval.
       !
       j_save = ne + 1

       do j = j_first, ne
          if ( le(x(ir),xe(j)) ) then
             j_save = j
             if ( ir == n ) then
                j_save = ne + 1
             end if
             exit
          end if
       end do
       !
       !  Have located first point beyond interval.
       !
       j = j_save

       nj = j - j_first
       !
       !  Skip evaluation if no points in interval.
       !
       if ( nj /= 0 ) then
          !
          !  Evaluate cubic at XE(J_FIRST:J-1).
          !
          call chfev ( x(ir-1), x(ir), f(ir-1), f(ir), d(ir-1), d(ir), &
               nj, xe(j_first:j-1), fe(j_first:j-1), next, ierc )

          if ( ierc < 0 ) then
             ierr = -5
             write ( *, '(a)' ) ' '
             write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!'
             write ( *, '(a)' ) '  Error return from CHFEV.'
             stop
          end if
          !
          !  In the current set of XE points, there are NEXT(2) to the right of X(IR).
          !
          if ( next(2) /= 0 ) then

             if ( ir < n ) then
                ierr = -5
                write ( *, '(a)' ) ' '
                write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!'
                write ( *, '(a)' ) '  IR < N.'
                stop
             end if
             !
             !  These are actually extrapolation points.
             !
             ierr = ierr + next(2)

          end if
          !
          !  In the current set of XE points, there are NEXT(1) to the left of X(IR-1).
          !
          if ( next(1) /= 0 ) then
             !
             !  These are actually extrapolation points.
             !
             if ( ir <= 2 ) then
                ierr = ierr + next(1)
             else

                j_new = -1

                do i = j_first, j - 1
                   if ( xe(i) < x(ir-1) ) then
                      j_new = i
                      exit
                   end if
                end do

                if ( j_new == -1 ) then
                   ierr = -5
                   write ( *, '(a)' ) ' '
                   write ( *, '(a)' ) 'SPLINE_PCHIP_VAL - Fatal error!'
                   write ( *, '(a)' ) '  Could not bracket the data point.'
                   stop
                end if
                !
                !  Reset J.  This will be the new J_FIRST.
                !
                j = j_new
                !
                !  Now find out how far to back up in the X array.
                !
                do i = 1, ir-1
                   if ( xe(j) < x(i) ) then
                      exit
                   end if
                end do
                !
                !  At this point, either XE(J) < X(1) or X(i-1) <= XE(J) < X(I) .
                !
                !  Reset IR, recognizing that it will be incremented before cycling.
                !
                ir = max ( 1, i-1 )

             end if

          end if

          j_first = j

       end if

       ir = ir + 1

       if ( n < ir ) then
          exit
       end if

    end do

    return
  end subroutine spline_pchip_val
  subroutine spline_quadratic_val ( ndata, tdata, ydata, tval, yval, ypval )

    !*****************************************************************************80
    !
    !! SPLINE_QUADRATIC_VAL evaluates a piecewise quadratic spline at a point.
    !
    !  Discussion:
    !
    !    Because of the simple form of a piecewise quadratic spline,
    !    the raw data points ( TDATA(I), YDATA(I)) can be used directly to
    !    evaluate the spline at any point.  No processing of the data
    !    is required.
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    24 October 1999
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    Input, integer(i4) NDATA, the number of data points defining
    !    the spline.  NDATA should be odd and at least 3.
    !
    !    Input, real(dp) TDATA(NDATA), YDATA(NDATA), the values of
    !    the independent and dependent variables at the data points.  The
    !    values of TDATA should be distinct and increasing.
    !
    !    Input, real(dp) TVAL, the point at which the spline is to
    !    be evaluated.
    !
    !    Output, real(dp) YVAL, YPVAL, the value of the spline and
    !    its first derivative dYdT at TVAL.  YPVAL is not reliable if TVAL
    !    is exactly equal to TDATA(I) for some I.
    !
    implicit none

    integer(i4) ndata

    real(dp) dif1
    real(dp) dif2
    integer(i4) left
    integer(i4) right
    real(dp) t1
    real(dp) t2
    real(dp) t3
    real(dp) tdata(ndata)
    real(dp) tval
    real(dp) y1
    real(dp) y2
    real(dp) y3
    real(dp) ydata(ndata)
    real(dp) ypval
    real(dp) yval

    if ( ndata < 3 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_QUADRATIC_VAL - Fatal error!'
       write ( *, '(a)' ) '  NDATA < 3.'
       stop
    end if

    if ( mod ( ndata, 2 ) == 0 ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_QUADRATIC_VAL - Fatal error!'
       write ( *, '(a)' ) '  NDATA must be odd.'
       stop
    end if
    !
    !  Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains, or is
    !  nearest to, TVAL.
    !
    call r8vec_bracket ( ndata, tdata, tval, left, right )
    !
    !  Force LEFT to be odd.
    !
    if ( mod ( left, 2 ) == 0 ) then
       left = left - 1
    end if
    !
    !  Copy out the three abscissas.
    !
    t1 = tdata(left)
    t2 = tdata(left+1)
    t3 = tdata(left+2)

    if ( le(t2,t1) .or. le(t3,t2) ) then
       write ( *, '(a)' ) ' '
       write ( *, '(a)' ) 'SPLINE_QUADRATIC_VAL - Fatal error!'
       write ( *, '(a)' ) '  T2 <= T1 or T3 <= T2.'
       stop
    end if
    !
    !  Construct and evaluate a parabolic interpolant for the data
    !  in each dimension.
    !
    y1 = ydata(left)
    y2 = ydata(left+1)
    y3 = ydata(left+2)

    dif1 = ( y2 - y1 ) / ( t2 - t1 )

    dif2 = ( ( y3 - y1 ) / ( t3 - t1 ) &
         - ( y2 - y1 ) / ( t2 - t1 ) ) / ( t3 - t2 )

    yval = y1 + ( tval - t1 ) * ( dif1 + ( tval - t2 ) * dif2 )
    ypval = dif1 + dif2 * ( 2.0_dp * tval - t1 - t2 )

    return
  end subroutine spline_quadratic_val
  subroutine timestamp ( )

    !*****************************************************************************80
    !
    !! TIMESTAMP prints the current YMDHMS date as a time stamp.
    !
    !  Example:
    !
    !    31 May 2001   9:45:54.872 AM
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    06 August 2005
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    None
    !
    implicit none

    character ( len = 8 ) ampm
    integer(i4) d
    integer(i4) h
    integer(i4) m
    integer(i4) mm
    character ( len = 9 ), parameter, dimension(12) :: month = (/ &
         'January  ', 'February ', 'March    ', 'April    ', &
         'May      ', 'June     ', 'July     ', 'August   ', &
         'September', 'October  ', 'November ', 'December ' /)
    integer(i4) n
    integer(i4) s
    integer(i4) values(8)
    integer(i4) y

    call date_and_time ( values = values )

    y = values(1)
    m = values(2)
    d = values(3)
    h = values(5)
    n = values(6)
    s = values(7)
    mm = values(8)

    if ( h < 12 ) then
       ampm = 'AM'
    else if ( h == 12 ) then
       if ( n == 0 .and. s == 0 ) then
          ampm = 'Noon'
       else
          ampm = 'PM'
       end if
    else
       h = h - 12
       if ( h < 12 ) then
          ampm = 'PM'
       else if ( h == 12 ) then
          if ( n == 0 .and. s == 0 ) then
             ampm = 'Midnight'
          else
             ampm = 'AM'
          end if
       end if
    end if

    write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
         d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

    return
  end subroutine timestamp

end module mo_spline