Astrophysics and Space Science, 139 (1987),
1–4
[Erratum: vol. 146, (No. 1, July 1988), 201]
TRANSFORMATION OF GEOCENTRIC TO GEODETIC
COORDINATES WITHOUT APPROXIMATIONS
KAZIMIERZ M. BORKOWSKI
Toruń Radio Astronomy Observatory, Nicolaus Copernicus
University, Toruń, Poland
(Received 27 July, 1987)
| Abstract. An exact and relatively simple analytical transform
of the rectangular coordinates to the geodetic coordinates is
presented. It does not involve any approximation and the accuracy
of practical calculations depends exclusively on the round-off
errors. The algorithm is based on one solution to the quartic
equation in tg(45° – ψ/2), where
ψ is the parametric (or
eccentric) latitude. |
Present-day astrogeodetic work requires a consistent pair of
transforms between geocentric and geodetic coordinates whose high
accuracy holds for any altitude and latitude. The geocentric
coordinates are known to be directly calculable from geodetic
coordinates through the use of simple formulae (e.g. Equations (1)
below). It is not so with the reverse transformation. Until
recently there prevailed a widespread believe that the reverse
transformation has not any elementary solution (e.g. Woolard and
Clemence 1966, Heiskanen and Moritz 1967, Izotov et al. 1974,
Lang 1974, Hlibowicki 1981, Taff 1985, Pick 1985, Baranov et al.
1986, The Astronomical Almanac for the Year 1987 (p. K11)). Thus
a rich variety of different algorithms have been devised to find
an approximation to the geodetic latitude, φ, and the vertical
hight (altitude) above a reference ellipsoid, h, given the
geocentric coordinates: the equatorial, r, and the polar, z,
component. Paul (1973) seems to have been the first to
demonstrate the existance of exact analytical solutions to the
geodetic coordinates. His algorithm was later considerably
improved by Heikkinen (1982) while Hedgley (1976) presented a
solution based on entirely different approach. The common feature
of these three algorithms is solving of a full quartic
(biquadratic) equation in some variable (see also Vaníček and
Krakiwsky 1982). This is normally tedious work though a final
practical procedure may sometimes be reduced to acceptable number
of formulas (e.g. Heikkinen 1982).
Unaware of the existance of the above mentioned publications
concerning the exact analytical solutions I set up to work out
one. Now it seems to be considerably simpler than the three of
Paul (1973), Hedgley (1976) and Heikkinen (1982) both in
derivation and in practical implementation. Thus it carries the
potential of becoming a standard in future references. The
simplicity of my solution comes from the fact that of the five
coefficients of the quartic equation only two are different from
zero or one (in the absolute value). In the following we shall
preset first a complete derivation and then give rules for
practical use.
The standard formulae relating the geodetic and geocentric
coordinates are:
and
| z =
(Cb2/a2 + h) sin φ, | (1b) |
where
| C =
a2/√[(a cos φ)2 + (b sin φ)2] | (2) |
is known as the radius of curveture of the reference ellipse at
the direction defined by φ, and a and b are the semiaxes of
the ellipse.
Eliminating the h parameter from Equations (1) directly leads to
(e.g. Vaníček and Krakiwsky, 1982)
| r – z/tg φ = (a2 –
b2)/√[a2 + (b tg φ)2]. | (3) |
This is the equation that, under somewhat different representations,
served to derive many of approximations met in practice.
It is easy to chack that the substitution:
| tg φ = a(1 – t2)/(2bt) | (4) |
transforms Equation (3) into the following quartic equation
| t4 + 2Et3 + 2Ft – 1 = 0, | (5) |
where
| E = [bz – (a2 – b2)]/(ar) | (6) |
and
| F = [bz + (a2 – b2)]/(ar). | (7) |
As is well known, the fourth-degree polynomials can always be
solved using only elementary functions, though these general
algorithms often prove impractical. In our case the simple form
of Equation (5) suits well the procedure known as the Ferrari's
solution (e.g. Korn and Korn 1961). It uses one of the three
roots of the qubic resolvent:
where
and
to solve for all four roots of the quartic equation. A real root
of Equation (8) is
| v = –(Q + √D)1/3
– (Q – √D)1/3 , | (11a) |
where
For D < 0 the expression (11a) is equivalent to
| v = 2 √(–P) cos{(1/3)arccos[–Q/(–P)3/2]}, | (11b) |
which avoids the complex arithmetic. Since the case of D < 0
envelopes points not further than about 45 km from the center of
geodetic reference ellipsoids, it is rather unimportant for
practical applications but remains a valuable tool in discussion
of the general properties of the transform.
Now, having calculated v, the four roots of Equation (5) can be
found from
| t = ±√[G2
+ (F – vG)/(2G – E)] – G, | (13) |
where
| G = [±√(E2 + v) + E]/2. | (14) |
Of the solutions defined by Equations (13) and (14) all are real for
D ≤ 0, and two are real for D > 0. Any one of the real solutions
can serve to find a distinct geodetic coordinate set through Equations (4) and (1):
|
φ = arctg[a(1 – t2)/(2bt)], | (15a) |
and e.g.
|
h = (r – at) cosφ + (z – b) sinφ, | (15b) |
and all of them satisfy Equations (1). This completes our derivation,
but a few more comments are now in order.
Among the four solutions there is always one real which falls
outside the range [–90°,+90°], the normal range of latitudes,
therefore generally the value of φ of Equation (15a) should be
placed in the right quarter by taking into account relative signs
of the numerator and denominator of the argument of the
arctangent function.
The rules can be given to select the single conventional
solution, which corresponds to that value of φ which lies in
the same quarter as does the geocentric latitude (arctg z/r), and
it can be demonstrated that there is only one such solution among
the four. However, we find it more practical to restrict the
considerations to the almost obvious case of a > b, and to the
northern latitudes (positive or zero z value), which together
greatly simplify the rules: just skip the lower signs of the
double signs in Equations (13) and (14) to take the positive square
roots only. On the southern hemisphere the solutions are of
course the same, save the sign reversed, thus the second
restriction does not really limit practical usefulness of our
algorithm.
Worth noticing is the fact that the presented algorithm lacks
symmetry about the equatorial plane though, of course, it gives
symmetric analytical results. It means that the numerical results
will have different round-off errors on either side of the
equator. This may be indicative of a potential for improvement of
the procedure, however, at this stage it is not clear how it
could be utilized.
The auxiliary parameter t can be shown to be the tangent
function of one half of the complement to 90° of the
eccentric latitude (called also the parametric or reduced
latitude), ψ. An algorithm exactly analogous to the presented
above can be obtained by constructing a quartic equation in
tg(ψ/2). Such a version exhibit a singularity at the equator
(for z = 0), instead of at the poles as in our case.
The fact that there are three other solutions besides the
desired one makes practically all iterative procedures sensitive
to the starting value for φ. If the geocentric latitude is
chosen for the start, a reasonable choice, there will always be
regions on the (r,z) plane (inside the curve D = 0) such that the
iterations will converge to a wrong solution. For example, the
point (r[m],z[m]) = (16000,2000) when referred to the present IAU
ellipsoid may be represented by the following (φ[deg],h[m])
coordinates: (69.1546512,–6351904.5), (–4.3033845,–6362215.0),
(–66.8170389,–6355613.9) and (–178.0477051,–6394174.1). In this
example the second solution for geodetic latitude lies much
closer to the geocentric latitude (about +7 degrees) than the
first (which is conventionally the only right).
The described algorithm, as embodied in Equations (6), (7) and
(9) through (15), has been implemented in FORTRAN programming
language, and extensively tested against the forward transform
(Equations (1)) for rectangular coordinates spanning a range of from
single meters to more than 500 000 km in both coordinates. Except
for the mentioned singularity at the poles no other practical
limitation was found. The program, in the form of a subroutine
named GEOD, can be obtained from the author on request.
Acknowledgements
The author wishes to express his gratitude
for financial support from the Nicolaus Copernicus Astronomical
Center of the Polish Academy of Sciences under contract CPBP-01.11.
References
Baranov V.N., Bojko, E.G., Krasnorylov, I.I., Mashimov, M.M.,
Plokhov, Yu.V., Urmaev, M.S., and Yashin, S.N.: 1986, Space
Geodesy, Nedra, Moscow, p. 27 (in Russian).
Hedgley, D.R.: 1976, NASA Techn. Rep. R-458.
Heikkinen, M.: 1982, Zeitschrift f. Vermess. 107, 207 (in German).
Heiskanen, W.A. and Moritz, H.: 1967, Physical Geodesy,
W.H. Freeman and Co., San Francisco, p. 181.
Hlibowicki, R. (ed.): 1981, Higher Geodesy and Geodetic
Astronomy, PWN, Warszawa, p. 273 (in Polish).
Izotov, A.A. (ed.): 1974, Fundamentals of Satellite
Geodesy, Nedra, Moscow (in Russian), p. 34.
Korn, G.A. and Korn, T.M.: 1961, Mathematical Handbook for
Scientists and Engineers, McGraw-Hill, New York, p. 24.
Lang, K.R.: 1974, Astrophysical Formulae, Springer-Verlag,
Berlin, p. 497.
Paul, M.K.: 1973, Bull. Géod., Nouv. Ser., No. 108, 135.
Pick, M.: 1985, Studia Geoph. Geod., 29, 112.
Taff, L.G.: 1985, Celestial Mechanics. A Computational Guide for
the Practitioner, J. Wiley and Sons, New York, Ch. 3.
Urmaev, M.S.: 1981, Orbital Methods of Space Geodesy,
Nedra, Moscow, p. 14 (in Russian).
Vaníček, P. and Krakiwsky, E.J.: 1982, Geodesy: The Concepts,
North Holland Publ. Co., Amsterdam, p. 324.
Woolard, E.W. and Clemence, G.M.: 1966, Spherical Astronomy,
Academic Press, New York, Ch. 3.