xyz2fl Subroutine

private subroutine xyz2fl(ax, ay, b, x, y, z, latitude, longitude)

Computes the transformation of Cartesian to geodetic coordinates on the surface of the ellipsoid assuming x,y,z are all non-negative Angular coordinates in radians

This is based on the C++ version

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: ax
real(kind=wp), intent(in) :: ay
real(kind=wp), intent(in) :: b
real(kind=wp), intent(in) :: x
real(kind=wp), intent(in) :: y
real(kind=wp), intent(in) :: z
real(kind=wp), intent(out) :: latitude
real(kind=wp), intent(out) :: longitude

Called by

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

Source Code

subroutine xyz2fl(ax, ay, b, x, y, z, latitude, longitude)

    real(wp),intent(in) :: ax, ay, b, x, y, z
    real(wp),intent(out) :: latitude, longitude

    real(wp) :: nom,den,dex,xme,rot
    real(wp) :: ax2,ay2,b2,Ex2,Ee2,lex2,lee2,mex,mee

    ! note: these could be precomputed:
    ax2  = ax*ax
    ay2  = ay*ay
    b2   = b*b
    Ex2  = ax2-b2
    Ee2  = ax2-ay2
    lex2 = Ex2/ax2
    lee2 = Ee2/ax2
    mex  = one-lex2
    mee  = one-lee2

    nom = mee*z
    xme = mee*x
    dex = xme*xme+y*y
    den = mex*sqrt(dex)
    rot = sqrt(dex)

    if (den==zero)  then
        latitude = halfpi
        longitude = zero
    else
        if (nom<=den) then
            latitude = atan(nom/den)
        else
            latitude = halfpi-atan(den/nom)
        end if
        if (y<=xme) then
            den = xme+rot
            longitude = two*atan(y/den)
        else
            den = y+rot
            longitude = halfpi - two*atan(xme/den)
        end if
    end if

end subroutine xyz2fl