module pchsp_95_m

  implicit none

contains

    ! PURPOSE: Set derivatives needed to determine the Hermite

    ! representation of the cubic spline interpolant to given data,

    ! with specified boundary conditions.

    ! Part of the "pchip" package.

    ! CATEGORY: E1A

    ! KEYWORDS: cubic hermite interpolation, piecewise cubic

    ! interpolation, spline interpolation

    ! DESCRIPTION: "pchsp" stands for "Piecewise Cubic Hermite Spline"

    ! Computes the Hermite representation of the cubic spline

    ! interpolant to the data given in X and F satisfying the boundary

    ! conditions specified by Ibeg, iend, vc_beg and VC_end.

    ! The resulting piecewise cubic Hermite function may be evaluated

    ! by "pchfe" or "pchfd".

    ! NOTE: This is a modified version of C. de Boor's cubic spline

    ! routine "cubspl".

    ! REFERENCE: Carl de Boor, A Practical Guide to Splines, Springer,

    ! 2001, pages 43-47

    use assert_eq_m, only: assert_eq

    real, intent(in):: x(:)

    ! independent variable values

    ! The elements of X must be strictly increasing:

    !                X(I-1) < X(I),  I = 2...N.

    !           (Error return if not.)

    ! (error if size(x) < 2)

    real, intent(in):: f(:)

    !     dependent variable values to be interpolated

    !  F(I) is value corresponding to X(I).

    INTEGER, intent(in):: ibeg

    !     desired boundary condition at beginning of data

    !        IBEG = 0  to set pchsp_95(1) so that the third derivative is con-

    !              tinuous at X(2).  This is the "not a knot" condition

    !              provided by de Boor's cubic spline routine CUBSPL.

    !              This is the default boundary condition.

    !        IBEG = 1  if first derivative at X(1) is given in VC_BEG.

    !        IBEG = 2  if second derivative at X(1) is given in VC_BEG.

    !        IBEG = 3  to use the 3-point difference formula for pchsp_95(1).

    !              (Reverts to the default boundary condition if size(x) < 3 .)

    !        IBEG = 4  to use the 4-point difference formula for pchsp_95(1).

    !              (Reverts to the default boundary condition if size(x) < 4 .)

    !          NOTES:

    !           1. An error return is taken if IBEG is out of range.

    !           2. For the "natural" boundary condition, use IBEG=2 and

    !              VC_BEG=0.

    INTEGER, intent(in):: iend

    !           IC(2) = IEND, desired condition at end of data.

    !  IEND may take on the same values as IBEG, but applied to

    !  derivative at X(N). In case IEND = 1 or 2, The value is given in VC_END.

    !          NOTES:

    !           1. An error return is taken if IEND is out of range.

    !           2. For the "natural" boundary condition, use IEND=2 and

    !              VC_END=0.

    REAL, intent(in), optional:: vc_beg

    ! desired boundary value, as indicated above.

    !           VC_BEG need be set only if IBEG = 1 or 2 .

    REAL, intent(in), optional:: vc_end

    ! desired boundary value, as indicated above.

    !           VC_END need be set only if Iend = 1 or 2 .

    real pchsp_95(size(x))

    ! derivative values at the data points

    !           These values will determine the cubic spline interpolant

    !           with the requested boundary conditions.

    !           The value corresponding to X(I) is stored in

    !                PCHSP_95(I),  I=1...N.

    ! LOCAL VARIABLES:

    INTEGER n ! number of data points

    integer ierr, ic(2)

    real vc(2)

    !-------------------------------------------------------------------

    end if

    end if

  END function pchsp_95

end module pchsp_95_m