PCHIP Subroutine

public subroutine PCHIP(x, y, N, PolyCoef, csWork)

Arguments

Type IntentOptional Attributes Name
real(kind=8), intent(in) :: x(:)
complex(kind=8), intent(in) :: y(:)
integer, intent(in) :: N
complex(kind=8), intent(inout) :: PolyCoef(4,N)
complex(kind=8), intent(inout) :: csWork(4,N)

Calls

proc~~pchip~~CallsGraph proc~pchip PCHIP proc~cspline CSPLINE proc~pchip->proc~cspline proc~fprime_interior_cmplx fprime_interior_Cmplx proc~pchip->proc~fprime_interior_cmplx proc~fprime_left_end_cmplx fprime_left_end_Cmplx proc~pchip->proc~fprime_left_end_cmplx proc~fprime_right_end_cmplx fprime_right_end_Cmplx proc~pchip->proc~fprime_right_end_cmplx proc~h_del h_del proc~pchip->proc~h_del proc~fprime_interior fprime_interior proc~fprime_interior_cmplx->proc~fprime_interior proc~fprime_left_end fprime_left_end proc~fprime_left_end_cmplx->proc~fprime_left_end proc~fprime_right_end fprime_right_end proc~fprime_right_end_cmplx->proc~fprime_right_end

Called by

proc~~pchip~~CalledByGraph proc~pchip PCHIP proc~cpchip cPCHIP proc~cpchip->proc~pchip proc~evaluatessp EvaluateSSP proc~evaluatessp->proc~cpchip proc~bellhopcore~2 BellhopCore proc~bellhopcore~2->proc~evaluatessp proc~influencecervenycart InfluenceCervenyCart proc~bellhopcore~2->proc~influencecervenycart proc~traceray2d~2 TraceRay2D proc~bellhopcore~2->proc~traceray2d~2 proc~influencecervenycart->proc~evaluatessp proc~readenvironment ReadEnvironment proc~readenvironment->proc~evaluatessp proc~reflect2d~2 Reflect2D proc~reflect2d~2->proc~evaluatessp proc~step2d Step2D proc~step2d->proc~evaluatessp proc~traceray2d~2->proc~evaluatessp proc~traceray2d~2->proc~reflect2d~2 proc~traceray2d~2->proc~step2d proc~bellhopcore BellhopCore proc~bellhopcore->proc~influencecervenycart program~bellhop BELLHOP program~bellhop->proc~bellhopcore~2 program~bellhop->proc~readenvironment program~bellhop3d BELLHOP3D program~bellhop3d->proc~readenvironment program~bellhop3d->proc~bellhopcore

Source Code

  SUBROUTINE PCHIP( x, y, N, PolyCoef, csWork )

    ! This implements the monotone piecewise cubic Hermite interpolating
    ! polynomial (PCHIP) algorithm. This is a new variant of monotone PCHIP
    ! (paper submitted to JASA). Also see;
    !
    ! F. N. Fritsch and J. Butland, "A Method for Constructing Local Monotone
    ! Piecewise Cubic Interpolants", SIAM Journal on Scientific and Statistical
    ! Computing, 5(2):300-304, (1984) https://doi.org/10.1137/0905021
    !
    ! F. N. Fritsch and R. E. Carlson. "Monotone Piecewise Cubic Interpolation",
    ! SIAM Journal on Numerical Analysis, 17(2):238-246, (1980) https://doi.org/10.1137/0717021
    !
    ! N is the number of nodes
    ! x is a vector of the abscissa values
    ! y is a vector of the associated ordinate values
    ! PolyCoef are the coefficients of the standard polynomial
    ! csWork is a temporary work space for the cubic spline

    INTEGER,          INTENT( IN  )   :: N
    REAL    (KIND=8), INTENT( IN  )   :: x(:)
    COMPLEX (KIND=8), INTENT( IN  )   :: y(:)
    COMPLEX (KIND=8), INTENT( INOUT ) :: PolyCoef( 4, N ), csWork( 4, N )

    INTEGER           :: ix, iBCBeg, iBCEnd
    REAL     (KIND=8) :: h1, h2
    COMPLEX  (KIND=8) :: del1, del2, f1, f2, f1prime, f2prime, fprimeT

    !  Precompute estimates of the derivatives at the nodes
    !  The vector PolyCoef(1,*) holds the ordinate values at the nodes
    !  The vector PolyCoef(2,*) holds the ordinate derivatives at the nodes

    IF ( N == 2 ) THEN

    ! handle special case of two data points seperately (linear interpolation)

    PolyCoef( 1, 1 ) = y( 1 )
    PolyCoef( 2, 1 ) = ( y( 2 ) - y( 1 ) ) / ( x( 2 ) - x( 1 ) )
    PolyCoef( 3, 1 ) = 0.0D0
    PolyCoef( 4, 1 ) = 0.0D0

    ELSE

    ! general case of more than two data points

    PolyCoef( 1, 1 : N ) = y( 1 : N )

    ! left endpoint (non-centered 3-point difference formula)

    CALL h_del( x, y, 2, h1, h2, del1, del2 )
    fprimeT = ( ( 2.0D0 * h1 + h2 ) * del1 - h1 * del2 ) / ( h1 + h2 )
    PolyCoef( 2, 1 ) = fprime_left_end_Cmplx( del1, del2, fprimeT )

    ! right endpoint (non-centered 3-point difference formula)

    CALL h_del( x, y, N - 1, h1, h2, del1, del2 )
    fprimeT = ( -h2 * del1 + ( h1 + 2.0D0 * h2 ) * del2 ) / ( h1 + h2 )
    PolyCoef( 2, N ) = fprime_right_end_Cmplx( del1, del2, fprimeT )

    ! compute coefficients of the cubic spline interpolating polynomial

    iBCBeg = 1   ! specified derivatives at the end points
    iBCEnd = 1
    csWork( 1, 1 : N ) = PolyCoef( 1, 1 : N )
    csWork( 2,     1 ) = PolyCoef( 2, 1 )
    csWork( 2,     N ) = PolyCoef( 2, N )
    CALL CSpline( x, csWork( 1, 1 ), N, iBCBeg, iBCEnd, N )

    ! interior nodes (use derivatives from the cubic spline as initial estimate)

    DO ix = 2, N - 1
       CALL h_del( x, y, ix, h1, h2, del1, del2 )
       ! check if the derivative from the cubic spline satisfies monotonicity
       PolyCoef( 2, ix ) = fprime_interior_Cmplx( del1, del2, csWork( 2, ix ) )
    END DO

    !                                                               2      3
    ! compute coefficients of std cubic polynomial: c0 + c1*x + c2*x + c3*x
    !

    DO ix = 1, N - 1
       h  = x( ix + 1 ) - x( ix )

       f1 = PolyCoef( 1, ix )
       f2 = PolyCoef( 1, ix + 1 )

       f1prime = PolyCoef( 2, ix )
       f2prime = PolyCoef( 2, ix + 1 )

       PolyCoef( 3, ix ) = ( 3.0D0 * ( f2 - f1 )  - h * ( 2.0D0 * f1prime + f2prime ) ) / h**2
       PolyCoef( 4, ix ) = ( h * ( f1prime + f2prime ) - 2.0D0 * ( f2 - f1 ) ) / h**3
    END DO

    END IF

    RETURN
  END SUBROUTINE PCHIP