dchfie Function

private pure function dchfie(x1, x2, f1, f2, d1, d2, a, b) result(integral)

Cubic Hermite Function Integral Evaluator

Called by dpchia to evaluate the integral of a single cubic (in Hermite form) over an arbitrary interval (A,B).

Note

There is no error return from this routine because zero is indeed the mathematically correct answer when X1==X2 .

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: x1

endpoints if interval of definition of cubic
real(kind=wp), intent(in) :: x2

endpoints if interval of definition of cubic
real(kind=wp), intent(in) :: f1

function values at the ends of the interval
real(kind=wp), intent(in) :: f2

function values at the ends of the interval
real(kind=wp), intent(in) :: d1

derivative values at the ends of the interval
real(kind=wp), intent(in) :: d2

derivative values at the ends of the interval
real(kind=wp), intent(in) :: a

endpoints of interval of integration
real(kind=wp), intent(in) :: b

endpoints of interval of integration

Return Value real(kind=wp)


Called by

proc~~dchfie~~CalledByGraph proc~dchfie pchip_module::dchfie proc~dpchia pchip_module::dpchia proc~dpchia->proc~dchfie

Source Code

    pure function dchfie (x1, x2, f1, f2, d1, d2, a, b) result(integral)

    real(wp),intent(in) :: x1 !! endpoints if interval of definition of cubic
    real(wp),intent(in) :: x2 !! endpoints if interval of definition of cubic
    real(wp),intent(in) :: f1 !! function values at the ends of the interval
    real(wp),intent(in) :: f2 !! function values at the ends of the interval
    real(wp),intent(in) :: d1 !! derivative values at the ends of the interval
    real(wp),intent(in) :: d2 !! derivative values at the ends of the interval
    real(wp),intent(in) :: a  !! endpoints of interval of integration
    real(wp),intent(in) :: b  !! endpoints of interval of integration
    real(wp) :: integral

    real(wp) :: dterm, fterm, h, phia1, phia2, &
                phib1, phib2, psia1, psia2, psib1, psib2, &
                ta1, ta2, tb1, tb2, ua1, ua2, ub1, ub2

    ! validity check input.
    if (x1 == x2) then
        integral = zero
    else
        h = x2 - x1
        ta1 = (a - x1) / h
        ta2 = (x2 - a) / h
        tb1 = (b - x1) / h
        tb2 = (x2 - b) / h

        ua1 = ta1**3
        phia1 = ua1 * (two - ta1)
        psia1 = ua1 * (three*ta1 - four)
        ua2 = ta2**3
        phia2 =  ua2 * (two - ta2)
        psia2 = -ua2 * (three*ta2 - four)

        ub1 = tb1**3
        phib1 = ub1 * (two - tb1)
        psib1 = ub1 * (three*tb1 - four)
        ub2 = tb2**3
        phib2 =  ub2 * (two - tb2)
        psib2 = -ub2 * (three*tb2 - four)

        fterm =   f1*(phia2 - phib2) + f2*(phib1 - phia1)
        dterm = ( d1*(psia2 - psib2) + d2*(psib1 - psia1) )*(h/six)

        integral = (half*h) * (fterm + dterm)
    endif

    end function dchfie