Transform geocentric coordinates to geodetic for a reference ellipsoid of specified form.
Status: support routine.
This routine is closely based on the GCONV2H subroutine by Toshio Fukushima (see reference).
The equatorial radius, A, can be in any units, but meters is the conventional choice.
The flattening, F, is (for the Earth) a value around 0.00335, i.e. around 1/298.
The equatorial radius, A, and the geocentric vector, XYZ, must be given in the same units, and determine the units of the returned height, HEIGHT.
If an error occurs (J<0), ELONG, PHI and HEIGHT are unchanged.
The inverse transformation is performed in the routine GD2GCE.
The transformation for a standard ellipsoid (such as WGS84) can more conveniently be performed by calling GC2GD, which uses a numerical code (1 for WGS84) to identify the required A and F values.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in) | :: | a | equatorial radius (Notes 2,4) |
||
| real(kind=wp), | intent(in) | :: | f | flattening (Note 3) |
||
| real(kind=wp), | intent(in), | dimension(3) | :: | xyz | geocentric vector (Note 4) |
|
| real(kind=wp), | intent(out) | :: | elong | longitude (radians, east +ve) |
||
| real(kind=wp), | intent(out) | :: | phi | latitude (geodetic, radians) |
||
| real(kind=wp), | intent(out) | :: | height | height above ellipsoid (geodetic, Note 4) |
||
| integer, | intent(out) | :: | j |
|
subroutine GC2GDE ( a, f, xyz, elong, phi, height, j )
implicit none
real(wp),intent(in) :: a !! equatorial radius (Notes 2,4)
real(wp),intent(in) :: f !! flattening (Note 3)
real(wp),dimension(3),intent(in) :: xyz !! geocentric vector (Note 4)
real(wp),intent(out) :: elong !! longitude (radians, east +ve)
real(wp),intent(out) :: phi !! latitude (geodetic, radians)
real(wp),intent(out) :: height !! height above ellipsoid (geodetic, Note 4)
integer,intent(out) :: j !! status:
!! * 0 = OK
!! * -1 = illegal F
!! * -2 = illegal A
real(wp) :: aeps2, e2, e4t, ec2, ec, b, x, y, z, p2, absz, p, &
s0, pn, zc, c0, c02, c03, s02, s03, a02, a0, a03, &
d0, f0, b0, s1, cc, s12, cc2
! -------------
! Preliminaries
! -------------
! Validate ellipsoid parameters.
if ( f<0.0_wp .or. f>=1.0_wp ) then
j = -1
return
else if ( a <= 0.0_wp ) then
j = -2
return
end if
! Functions of ellipsoid parameters (with further validation of F).
aeps2 = a*a*1.0e-32_wp
e2 = (2.0_wp-f)*f
e4t = e2*e2*1.5_wp
ec2 = 1.0_wp-e2
if ( ec2 <= 0.0_wp ) then
j = -1
return
end if
ec = sqrt(ec2)
b = a*ec
! Cartesian components.
x = xyz(1)
y = xyz(2)
z = xyz(3)
! Distance from polar axis squared.
p2 = x*x + y*y
! Longitude.
if ( p2>0.0_wp ) then
elong = atan2(y,x)
else
elong = 0.0_wp
end if
! Unsigned z-coordinate.
absz = abs(z)
! Proceed unless polar case.
if ( p2>aeps2 ) then
! Distance from polar axis.
p = sqrt(p2)
! Normalization.
s0 = absz/a
pn = p/a
zc = ec*s0
! Prepare Newton correction factors.
c0 = ec*pn
c02 = c0*c0
c03 = c02*c0
s02 = s0*s0
s03 = s02*s0
a02 = c02+s02
a0 = sqrt(a02)
a03 = a02*a0
d0 = zc*a03 + e2*s03
f0 = pn*a03 - e2*c03
! Prepare Halley correction factor.
b0 = e4t*s02*c02*pn*(a0-ec)
s1 = d0*f0 - b0*s0
cc = ec*(f0*f0-b0*c0)
! Evaluate latitude and height.
phi = atan(s1/cc)
s12 = s1*s1
cc2 = cc*cc
height = (p*cc+absz*s1-a*sqrt(ec2*s12+cc2))/sqrt(s12+cc2)
else
! Exception: pole.
phi = dpi/2.0_wp
height = absz-b
end if
! Restore sign of latitude.
if ( z<0.0_wp ) phi = -phi
! OK status.
j = 0
end subroutine GC2GDE