geodetic_to_cartesian_triaxial_2 Subroutine

public pure subroutine geodetic_to_cartesian_triaxial_2(a, b, c, lat, long, h, r)

Geodetic to Cartesian for Triaxial Ellipsoid.

References

  • S. Bektas, "Geodetic Computations on Triaxial Ellipsoid", International Journal of Mining Science (IJMS), Volume 1, Issue 1, June 2015, p 25-34

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) :: lat

latitude (rad)
real(kind=wp), intent(in) :: long

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

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

Cartesian coordinates (x,y,z)

Source Code

    pure subroutine geodetic_to_cartesian_triaxial_2(a,b,c,lat,long,h,r)

    implicit none

    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) :: lat  !! latitude (rad)
    real(wp),intent(in) :: long !! longitude (rad)
    real(wp),intent(in) :: h    !! altitude
    real(wp),dimension(3),intent(out) :: r  !! Cartesian coordinates (x,y,z)

    real(wp) :: ex2,ee2,v,a2,clat,slat,clon,slon,omee2,omex2

    a2    = a * a
    ex2   = (a2-c**2)/a2
    ee2   = (a2-b**2)/a2
    clat  = cos(lat)
    slat  = sin(lat)
    clon  = cos(long)
    slon  = sin(long)
    omee2 = 1.0_wp-ee2
    omex2 = 1.0_wp-ex2
    v     = a/sqrt(1.0_wp-ex2*slat**2-ee2*clat**2*slon**2)

    r = [(v+h)*clon*clat, &
         (v*omee2+h)*slon*clat, &
         (v*omex2+h)*slat ]

    end subroutine geodetic_to_cartesian_triaxial_2