parabolic Function

private pure function parabolic(xa, xb, xc, fa, fb, fc) result(xp)

Parabolic interpolation, with the option to fall back to bisection if the parabolic interpolation fails. Used by rbp.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: xa
real(kind=wp), intent(in) :: xb
real(kind=wp), intent(in) :: xc
real(kind=wp), intent(in) :: fa
real(kind=wp), intent(in) :: fb
real(kind=wp), intent(in) :: fc

Return Value real(kind=wp)


Calls

proc~~parabolic~~CallsGraph proc~parabolic parabolic proc~bisect bisect proc~parabolic->proc~bisect

Called by

proc~~parabolic~~CalledByGraph proc~parabolic parabolic proc~rbp rbp_solver%rbp proc~rbp->proc~parabolic

Source Code

    pure function parabolic(xa,xb,xc,fa,fb,fc) result(xp)

    real(wp),intent(in) :: xa,xb,xc
    real(wp),intent(in) :: fa,fb,fc
    real(wp) :: xp

    real(wp) :: a,b,c,denom
    logical :: error  !! used to check for errors in parabolic interpolation

    a = (fa-fc)/(xa-xc)/(xa-xb) + (fc-fb)/(xb-xc)/(xa-xb)
    b = (fc-fa)*(xb-xc)/(xa-xc)/(xa-xb) - (fc-fb)*(xa-xc)/(xb-xc)/(xa-xb)
    c = fc

    ! root of the parabola:
    denom = b + sign(1.0_wp, b)*sqrt(b**2 - 4*a*c)
    error = denom == 0.0_wp
    if (.not. error) then
        xp = xc - 2.0_wp*c / denom
        error = xp==xa .or. xp==xb .or. xp==xc
    end if
    if (error) then
        ! failure in parabolic interpolation - fall back to bisection
        if (fa*fc<0.0_wp) then
            xp = bisect(xa,xc)
        else
            xp = bisect(xc,xb)
        end if
    end if

    end function parabolic