special_cases Subroutine

private subroutine special_cases(a, b, c, x, y, z, phi, lambda, h, done)

Special cases for lat/lon/altitude

Arguments

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

ellipsoid radii a >= b >= c
real(kind=wp), intent(in) :: b

ellipsoid radii a >= b >= c
real(kind=wp), intent(in) :: c

ellipsoid radii a >= b >= c
real(kind=wp), intent(in) :: x

Cartesian x coordinate
real(kind=wp), intent(in) :: y

Cartesian y coordinate
real(kind=wp), intent(in) :: z

Cartesian z coordinate
real(kind=wp), intent(out) :: phi

latitude (rad)
real(kind=wp), intent(out) :: lambda

longitude (rad)
real(kind=wp), intent(out) :: h

altitude
logical, intent(out) :: done

true if one of the special cases was computed

Called by

proc~~special_cases~~CalledByGraph proc~special_cases geodesy_module::special_cases proc~cartesian_to_geodetic_triaxial_2 geodesy_module::cartesian_to_geodetic_triaxial_2 proc~cartesian_to_geodetic_triaxial_2->proc~special_cases proc~cartesianintogeodetici geodesy_module::CartesianIntoGeodeticI proc~cartesianintogeodetici->proc~special_cases proc~cartesianintogeodeticii geodesy_module::CartesianIntoGeodeticII proc~cartesianintogeodeticii->proc~special_cases

Source Code

    subroutine special_cases(a,b,c,x,y,z,phi,lambda,h,done)

    real(wp),intent(in)  :: a      !! ellipsoid radii `a >= b >= c`
    real(wp),intent(in)  :: b      !! ellipsoid radii `a >= b >= c`
    real(wp),intent(in)  :: c      !! ellipsoid radii `a >= b >= c`
    real(wp),intent(in)  :: x      !! Cartesian x coordinate
    real(wp),intent(in)  :: y      !! Cartesian y coordinate
    real(wp),intent(in)  :: z      !! Cartesian z coordinate
    real(wp),intent(out) :: phi    !! latitude (rad)
    real(wp),intent(out) :: lambda !! longitude (rad)
    real(wp),intent(out) :: h      !! altitude
    logical,intent(out)  :: done   !! true if one of the special cases was computed

    logical :: x0, y0, z0

    real(wp),parameter :: zero_tol  = 10.0_wp * epsilon(1.0_wp)  !! zero tolerance for singularities

    x0 = abs(x) <= zero_tol
    y0 = abs(y) <= zero_tol
    z0 = abs(z) <= zero_tol

    if (x0 .and. y0 .and. z0) then ! center of the body
        phi    = zero
        lambda = zero
        h      = -c    ! just pick this value
        done = .true.
        return
    else if (x0 .and. y0) then ! (on the z-axis)
        if (z>=zero) then
            phi    = halfpi
            lambda = zero
            h      = z-c
        else
            phi    = -halfpi
            lambda = zero
            h      = -(z+c)
        end if
        done = .true.
        return
    else if (x0 .and. z0) then  ! on the y-axis
        if (y>=zero) then
            phi    = zero
            lambda = halfpi
            h      = y-b
        else
            phi    = zero
            lambda = -halfpi
            h      = -(y+b)
        end if
        done = .true.
        return
    else if (y0 .and. z0) then  ! on the x-axis
        if (x>=zero) then
            phi    = zero
            lambda = zero
            h      = x-a
        else
            phi    = zero
            lambda = pi
            h      = -(x+a)
        end if
        done = .true.
        return
    end if

    phi    = zero
    lambda = zero
    h      = zero
    done = .false.

    end subroutine special_cases