doxy/html/pchsp__95__m_8_f90_source.html

module pchsp_95_m


  implicit none


contains


  function pchsp_95(x, f, ibeg, iend, vc_beg, vc_end)


    ! 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:

    real wk(2, size(x)) ! real array of working storage

    INTEGER n ! number of data points

    integer ierr, ic(2)

    real vc(2)


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


    n = assert_eq(size(x), size(f), "pchsp_95 n")

    if ((ibeg == 1 .or. ibeg == 2) .and. .not. present(vc_beg)) then

       print *, "vc_beg required for IBEG = 1 or 2"

       stop 1

    end if

    if ((iend == 1 .or. iend == 2) .and. .not. present(vc_end)) then

       print *, "vc_end required for IEND = 1 or 2"

       stop 1

    end if

    ic = (/ibeg, iend/)

    if (present(vc_beg)) vc(1) = vc_beg

    if (present(vc_end)) vc(2) = vc_end

    call pchsp(ic, vc, n, x, f, pchsp_95, 1, wk, size(wk), ierr)

    if (ierr /= 0) stop 1


  END function pchsp_95


end module pchsp_95_m