bisection_special_2 Subroutine

private subroutine bisection_special_2(cz, Xo, Yo, tol, n, m, Gm)

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: cz
real(kind=wp), intent(in) :: Xo
real(kind=wp), intent(in) :: Yo
real(kind=wp), intent(in) :: tol
integer, intent(out) :: n
real(kind=wp), intent(out) :: m
real(kind=wp), intent(out) :: Gm

Called by

proc~~bisection_special_2~~CalledByGraph proc~bisection_special_2 geodesy_module::bisection_special_2 proc~cartesian_to_geodetic_triaxial geodesy_module::cartesian_to_geodetic_triaxial proc~cartesian_to_geodetic_triaxial->proc~bisection_special_2

Source Code

subroutine bisection_special_2(cz, Xo, Yo, tol, n, m, Gm)

    real(wp),intent(in) :: cz, Xo, Yo, tol
    integer,intent(out) :: n
    real(wp),intent(out) :: m, Gm

    real(wp) :: d1,Gd1,d2,d,MM

    d1 = -one+Yo
    Gd1 = (cz*Xo*cz*Xo)/((cz+d1)*(cz+d1))
    d2 = -one+sqrt(cz*Xo*cz*Xo+Yo*Yo)
    d = (d2 - d1)/two
    n = 0
    m = -two

    do while (d > tol)
        n = n + 1
        MM = m
        m = d1 + d
        Gm = ((cz*Xo*cz*Xo)/((cz+m)*(cz+m)))+((Yo*Yo)/((one+m)**2))-one
        if (MM == m + tol .or. Gm == zero) return
        if (sign(one,Gm) == sign(one,Gd1)) then
            d1 = m
            Gd1 = Gm
        else
            d2 = m
        end if
        d = (d2 - d1)/two
    end do

    n = n + 1
    m = d1 + d
    Gm = ((cz*Xo*cz*Xo)/((cz+m)*(cz+m)))+((Yo*Yo)/((one+m)**2))-one

end subroutine bisection_special_2