geodetic_to_cartesian_triaxial Subroutine

public subroutine geodetic_to_cartesian_triaxial(ax, ay, b, phi, lambda, h, r)

Function computes the Cartesian coordinates given the geodetic latitude (phi), longitude (lambda) and height (h) of a point related to an ellipsoid defined by its three semiaxes ax, ay and b

History

  • Jacob Williams, 10/29/2022 : Fortran verison of this algorithm, based on the Matlab (v1.0 01/03/2019) code.

Arguments

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

semiaxes (0 < b <= ay <= ax)
real(kind=wp), intent(in) :: ay

semiaxes (0 < b <= ay <= ax)
real(kind=wp), intent(in) :: b

semiaxes (0 < b <= ay <= ax)
real(kind=wp), intent(in) :: phi

geodetic latitude (radians)
real(kind=wp), intent(in) :: lambda

geodetic longitude (radians)
real(kind=wp), intent(in) :: h

geodetic height
real(kind=wp), intent(out), dimension(3) :: r

Cartesian position vector [x,y,z]

Source Code

subroutine geodetic_to_cartesian_triaxial(ax, ay, b, phi, lambda, h, r)

    real(wp),intent(in) :: ax !! semiaxes (0 < b <= ay <= ax)
    real(wp),intent(in) :: ay !! semiaxes (0 < b <= ay <= ax)
    real(wp),intent(in) :: b  !! semiaxes (0 < b <= ay <= ax)
    real(wp),intent(in) :: phi !! geodetic latitude (radians)
    real(wp),intent(in) :: lambda !! geodetic longitude (radians)
    real(wp),intent(in) :: h !! geodetic height
    real(wp),dimension(3),intent(out) :: r  !! Cartesian position vector [x,y,z]

    real(wp) :: ee2,ex2,N,cp,sp,cl,sl

    cp  = cos(phi)
    sp  = sin(phi)
    cl  = cos(lambda)
    sl  = sin(lambda)
    ee2 = (ax*ax-ay*ay)/(ax*ax)
    ex2 = (ax*ax-b*b)/(ax*ax)
    N   = ax/sqrt(one - ex2*sp*sp - ee2*cp*cp*sl*sl)

    r = [(N+h)*cp*cl, &
         (N*(one-ee2)+h)*cp*sl, &
         (N*(one-ex2)+h)*sp ]

end subroutine geodetic_to_cartesian_triaxial