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:

r = (C + h) cos φ(1a)

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:

v3 + 3Pv + 2Q = 0,(8)
where
P = (4/3)(EF + 1)(9)
and
Q = 2(E2 – F2),(10)

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
D = P3 + Q2.(12)

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.


[See also this paper.]