geodetic_to_cartesian Subroutine

public subroutine geodetic_to_cartesian(a, b, glat, lon, h, r)

Geodetic latitude, longitude, and height to Cartesian position vector.

References

  1. E. D. Kaplan, "Understanding GPS: Principles and Applications", Artech House, 1996.

Arguments

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

geoid semimajor axis [km]

real(kind=wp), intent(in) :: b

geoid semiminor axis [km]

real(kind=wp), intent(in) :: glat

geodetic latitude [rad]

real(kind=wp), intent(in) :: lon

longitude [rad]

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

geodetic altitude [km]

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

Cartesian position vector [x,y,z]


Source Code

    subroutine geodetic_to_cartesian(a,b,glat,lon,h,r)

    implicit none

    real(wp),intent(in) :: a                !! geoid semimajor axis [km]
    real(wp),intent(in) :: b                !! geoid semiminor axis [km]
    real(wp),intent(in) :: glat             !! geodetic latitude [rad]
    real(wp),intent(in) :: lon              !! longitude [rad]
    real(wp),intent(in) :: h                !! geodetic altitude [km]
    real(wp),dimension(3),intent(out) :: r  !! Cartesian position vector [x,y,z]

    real(wp) :: e2,slat,clat,slon,clon,tlat,ome2,d,q,aod

    slat    = sin(glat)
    clat    = cos(glat)
    tlat    = tan(glat)
    slon    = sin(lon)
    clon    = cos(lon)
    e2      = 1.0_wp - (b*b)/(a*a)
    ome2    = 1.0_wp - e2
    d       = sqrt( 1.0_wp + ome2*tlat*tlat )
    q       = sqrt( 1.0_wp - e2*slat*slat   )
    aod     = a/d

    r(1) = aod*clon + h*clon*clat
    r(2) = aod*slon + h*slon*clat
    r(3) = a*ome2*slat/q + h*slat

    end subroutine geodetic_to_cartesian