The problem of the transformation is reduced to solving of the equation 2 sin(ψ – Ω) = c sin 2ψ, where Ω = arctg [bz/(ar)], c = (a2 – b2)/√{(ar)2 + (bz)2}, a and b are the semi-axes of the reference ellipsoid, and z and r are the polar and equatorial, respectively, components of the position vector in the Cartesian system of coordinates. Then, the geodetic latitude is found as φ = arctg [(a/b) tg ψ], and the height above the ellipsoid as h = (r – a cosψ) cosφ + (z – b sinψ) sinφ. Two accurate closed solutions are proposed of which one is approximate in nature and the other is exact. They are shown to be superior to others, found in literature and in practice, in both or either accuracy and/or simplicity.

The transformation of coordinates between the geodetic and Cartesian system is a matter of everyday practice in geodesy. While conversion of the geodetic coordinates, the latitude φ, the longitude λ and the height above a reference ellipsoid h, to the geocentric Cartesian x, y and z coordinates is done straightforwardly the reverse conversion poses certain difficulties. Usually it is being performed by successive approximations in an iterative process (e.g. Heiskanen and Moritz, 1967, Bartelme and Meissl, 1975, Urmaev, 1981, Hlibowicki, 1981), or using closed approximative formulae cast as series expansions (e.g. Morrison and Pines, 1961, Izotov et al., 1974, Long, 1976, Taff, 1985) or other expressions (e.g. Pick, 1985, Baranov et al., 1986). Even though there exist closed exact solutions for transformation of the Cartesian coordinates onto the geodetic coordinates (Paul, 1973, Hedgley, 1976, Heikkinen, 1982, Vaniček and Krakiwsky, 1982) many authors seem to believe to the contrary (e.g. Woolard and Clemence, 1966, Heiskanen and Moritz, 1967, Izotov et al., 1974, Taff, 1985, Baranov et al., 1986, Astronomical Almanac, 1987). The exact solutions are in principle error free, however, in practice errors arise due to rounding-off in relatively complicated algorithms. Nevertheless, in certain applications, they might be preferred to simpler approximative solutions both due to their higher accuracy and completeness (range of latitudes and heights).

In this report we offer two solutions: a very simple approximation algorithm, which gives accuracies well above any practical needs, and in fact approaching the precision of exact algorithms, and one truly exact, independent of the so far published. The derivation of the latter is so simple that it may serve as a standard reference. The simplicity comes from the fact that in our quartic (biquadratic) equation only two of the five terms have coefficients different from zero or ±1 (other authors solve a general form of the quartic equation). We shall also present a comparison of performance of a few of more popular algorithms, along with ours.

To start with, observe that the Cartesian coordinates can be obtained from the geodetic ones through the following formulae

r = a cosψ + h cosφ     and

(1a)

z = b sinψ + h sinφ,

(1b)

where a (b) is the semimajor (semiminor) axis of the reference ellipsoid and r, the equatorial component, is split into x and y by multiplying Eq.(1a) by the cosine and sine, respectively, function of the longitude. The ψ parameter, sometimes called parametric or reduced latitude, is calculated from

ψ = arctg([b/a] tgφ).

(2)

The Eqs. (1) are equivalent to the more common ones

r = (N + h) cosφ

(3a)

z = [N(1 – e2) + h] sinφ,

(3b)

where N = a/(1 – e² sin²φ)^1/2, and e² = (a² – b²)/a². Now, given r and z, the problem of finding φ and h can be shown to reduce to solving of the following equation:

2 sin(ψ – Ω) – c sin 2ψ = 0,

(4)

where Ω = arctg [bz/(ar)] and c = (a² – b²)/[(ar)² + (bz)²]^1/2. This equation obtains from geometrical considerations. Equivalently, it can be derived directly from Eqs. (1) by eliminating the h parameter. A ψ that solves Eq. (4), by virtue of Eq. (2), is used to calculate the geodetic latitude:

φ = arctg[(a/b) tgψ].

(5)

This, in turn, allows the use of both or any one of the two equations (1) to find the height, e.g.

h = (r – a cosψ) cosφ + (z – b sinψ) sinφ.

(6)

We prefer this form for it has similar accuracy over all range of latitudes unlike the formulas met in practice and involving division by the sine or cosine of the latitude. Solving of the Eq. (4) should not present any difficulty in view of the fact that nowadays even handheld calculators have a built-in capability of finding the zeroes of transcendental functions of the type of the left side of (4). A convenient starting value for ψ can be chosen to equal arctg [az/(br)], which, as evidenced by Eqs. (1), is exact on the surface of the ellipsoid (with h = 0). Applying the Newton's method results in a rapidly converging series of approximations. This method uses the derivative of a function to be solved with respect to the variable (the ψ here), which in our case is

w = 2[cos(ψ′ – Ω) – c cos 2ψ′],

(7)

to find next approximation:

ψ = ψ′ – [2sin(ψ′ – Ω) – c sin 2ψ′]/w,

(8)

where the prime indicates the starting value or previous approximation. Computer tests have shown that just two iterations provide accuracies by far exceeding the needs of even very demanding practical user. The accuracies are of the order of 1E–9 rad in the geodetic latitude for all points further than a thousand of kilometers away from the origin of the Cartesian system. Another way of saying roughly the same is to claim that it is some six orders of magnitude more accurate than other algorithms met in literature, save the few truly exact solutions.

Table 1 presents errors of position calculated with the use of 10 different algorithms. The errors are defined to be the distances between the true position and the one calculated. To make the comparison meaningful the first three iterative algorithms (based on Heiskanen and Moritz (1967), Astronomical Almanac (1987), and Bartelme and Meissl (1975), respectively) were restricted to only two iterations (their performance for the third iteration is shown in brackets). The algorithms numbered 4 to 7 were taken from Izotov at al. (1974), Long (1976; Eq. (8) and (15) therein), Baranov et al. (1986) and Pick (1985), respectively. The Baranov et al. formula has been corrected for a misprint in the cited source (the a parameter in their expression for tgφ has been inadvertadly squared). It is clear from this table that our relatively simple algorithm, whose performance is listed in the column headed 8, compares favourably even with the exact solutions (the two rightmost columns in Table 1, of which the first corresponds to the algorithm given by Heikkinen (1982) and the second is ours, still to be presented, see also Borkowski 1987). Although our approximative solution is iterative in principle, due to its exceptionally fast convergence, all the process of iterations is reduced to two formulas which are essentially the first (Eq. (8) with ψ′ = Ω) and the second approximation (Eq. (8) itself). With such an understanding this algorithm can be safely treated also as a closed solution. It was used for calculation of entries to Table 1.

Eq.(4), when expressed in tg(π/4 – ψ/2) ≡ t, leads to the fourth-degree polynomial (biquadratic or quartic equation):

t4 + 2Et3 + 2Ft – 1 = 0,

(9)

where

E = [bz – (a2 – b2)]/(ar)     and

(10)

F = [bz + (a2 – b2)]/(ar).

(11)

Applying the Ferrari's solution (e.g. Korn and Korn, 1961) to Eq. (9) one obtains:

t = ±√{G2 + (F – vG)/(2G – E)} – G,

(12)

where

G = [±(E2 + v)1/2 + E]/2,

(13)

v = (D1/2 – Q)1/3 – (D1/2 + Q)1/3,

(14a)

D = P3 + Q2,

(15)

P = (4/3)(EF + 1)    and

(16)

Q = 2(E2 – F2).

(17)

For D < 0, to avoid the complex arithmetic, it is recommended to use the following equivalent of Eq.(14a)

v = 2 √{–P} cos{(1/3)arccos[Q/(–P)3/2]}.

(14b)

The case of D < 0 will be of marginal use in practice since D is a positive quantity anywhere further than about 45 km from the Earth center. Now, any real of the four solutions (corresponding to the four combinations of +/– signs in Eqs. (12) and (13)) represents a distinctive set of geodetic coordinates obtainable through:

φ = arctg [a(1 – t2)/(2bt)]    and

(18)

h = (r – at) cosφ + (z – b) sinφ .

(19)

As a practical matter, one does not need all the four solutions even though each of them satisfies Eqs. (1) (they are all real for D < 0 and two are real for D > 0). To obtain readily programmable algorithm yielding the single desired solution just skip out the double signs in Eqs. (12) and (13) to leave only positive square roots there. Then it will work provided a > b and φ > 0. Noteworthy, this algorithm will work properly also for southern latitudes (z < 0) if one additionally preassigns the z-coordinate sign to the minor semiaxis, b.

Table 1

Position displacement errors (in mm) for ten algorithms. The numbers in brackets are for 3 iterations with procedures of Heiskanen and Moritz (1967), Astronomical Almanac (1987) and Bartelme and Meissl (1975), algorithms No. 1, 2 and 3, respectively. Other algorithms are based on approximate closed solutions: Izotov et al., (1974; No. 4), Long (1975; 5), Baranov et al., (1986; 6), Pick (1985; 7) and this paper (8), and exact solutions: Heikkinen (1982; 9) and this paper (10). The bottom line lists respective timing measurements in units of 0.1 ms.

φ [deg]	h [km]	1	2	3	4	5	6	7	8	9	10
89.000	100000	109(0)	0(0)	1(0)	721	190	1	925	.000000	.000015	.000000
89.000	10000	460(1)	0(0)	22(0)	721	158	0	126	.000004	.000002	.000000
89.000	1000	228(1)	0(0)	25(0)	721	857	0	0	.000000	.000000	.000000
89.000	0	0(0)	0(0)	–(–)	721	1219	0	0	.000000	.000000	.000000
89.000	–1000	431(3)	0(0)	64(0)	721	1836	0	1	.000000	.000001	.000001
70.000	100000	89(0)	0(0)	0(0)	550	112	248	13722	.000021	.000017	.000000
70.000	10000	377(1)	1(0)	10(0)	635	105	105	2094	.000000	.000000	.000000
70.000	1000	187(1)	6(0)	11(0)	660	370	5	6	.000000	.000001	.000000
70.000	0	0(0)	9(0)	–(–)	646	506	0	0	.000000	.000001	.000000
70.000	–1000	353(2)	12(0)	29(0)	611	734	10	34	.000000	.000000	.000001
45.000	100000	30(0)	1(0)	0(0)	167	24	452	10075	.000000	.000000	.000000
45.000	10000	128(0)	37(0)	0(0)	365	5	191	2683	.000000	.000000	.000000
45.000	1000	63(0)	181(1)	0(0)	587	215	9	14	.000000	.000001	.000001
45.000	0	0(0)	242(1)	–(–)	636	327	0	0	.000000	.000000	.000000
45.000	–1000	120(0)	340(1)	0(0)	695	519	18	81	.000000	.000000	.000000
20.000	100000	2(0)	2(0)	0(0)	13	22	91	13641	.000007	.000000	.000015
20.000	10000	7(0)	73(0)	10(0)	118	25	38	2090	.000000	.000000	.000002
20.000	1000	3(0)	358(2)	11(0)	424	72	2	7	.000000	.000001	.000001
20.000	0	0(0)	478(3)	–(–)	546	121	0	0	.000000	.000001	.000001
20.000	–1000	6(0)	672(5)	28(0)	736	205	4	35	.000000	.000001	.000000
1.000	100000	0(0)	0(0)	1(0)	0	1	0	922	.000015	.000015	.000009
1.000	10000	0(0)	5(0)	22(0)	5	1	0	126	.000000	.000000	.000000
1.000	1000	0(0)	25(0)	24(0)	25	10	0	0	.000000	.000001	.000000
1.000	0	0(0)	33(0)	–(–)	33	17	0	0	.000000	.000000	.000000
1.000	–1000	0(0)	47(0)	63(0)	47	27	0	1	.000000	.000001	.000000
Exec. time:		56(78)	51(70)	50(63)	44	34	46	146	74	56	63

We have tested our exact algorithm, as embodied in Eqs. (10) through (19), over wide range of (r,z) plane (save the z-axis, where the quantities E and F are obviously indeterminate) without encountering any problems. As to the round-off errors, they are certainly not serious for practical applications (we have always used the DOUBLE PRECISION arithmetic in FORTRAN programming language), though we found them considerably larger in a close vicinity of the r and z axes than elsewhere. The main source of these errors belongs to evaluation of Eq. (14a), which involves the difference of cubic roots (computed as the exponent of one third of the natural logarithm). The simple way to circumvent this problem is to correct a result obtained from Eq. (14a), say v′, by using the cubic resolvent (of which v′ is one of the three roots):

v = –(v′3 + 2Q)/(3P).

(20)

It can be easily proved that v from the above equation is more accurate than v′ if v² < P, which is the case everywhere outside the radial distance from the Earth's centre of about 70 km. The solution of Heikkinen (1982), which also involves the cubic root operation, is free from the considered weakness, however, it is not designed to handle the case of D < 0. For D > 0 and with the correction (20) in our algorithm the two solutions can be considered equally good.

To compare the effectiveness of the discussed algorithms we have measured the amount of CPU time required on an IBM PC type machine using the Lahey Computer Systems FORTRAN compiler. The results, given in units of 0.1 ms in the last row of Table 1, show that all times are comparable despite considerable differences in complexity of the algorithms. This is so, because the main contribution to the computation time is due to the trigonometric functions and they are a few in number in each algorithm. For example, the evaluation of the cosine function takes about 0.6 ms, or 6 units of Table 1. Since the computing is really cheap, no special attempt was undertaken to optimize the execution times. Also, sizes of the algorithms in our implementation are comparable and range from 7 to 15 FORTRAN cards, with the single exception of Pick's solution, which anyway seems to be the worst in every respect.

To conclude, we find that the iterative algorithms, such like these published by Heiskanen and Moritz (1967) or in Astronomical Almanac (1987), are the simplest to program and give acceptable results with 3 to 4 iterations. With one or two additional iterations they can also be used when higher accuracy is required. However, due to round-off errors attempting more than 6 iterations may not improve the result obtained so far, which may be slightly worse than the one obtainable with algorithms based on exact solutions. In this case the two algorithms presented in this paper, or the one of Heikkinen (1982), seem to be more appropriate. They guarantee accuracies a few orders below 1 mm over all practically useful area on the (r,z)-plane, from as close as a few tens of kilometers from the very centre of the Earth to many Earth's radii out into the space. For that matter the correction (20) to our mathematically exact solution may be safely omitted. For more general solutions of the transformation between the Cartesian and geodetic coordinates the only readily available algorithm is our exact solution (presented originally in full detail, unfortunately with an annoying misprint in the equation corresponding to Eq. (12), in Borkowski (1987)).

FORTRAN listings of all algorithms analysed in this report may be requested from the author.

Acknowledgements

I would like to thank dr B . Benciolini for helpful comments on an earlier version of this report. This work was financially supported through the Polish Government research problem RPBR No. RR.I.11/2.

REFERENCES

Astronomical Almanac for the Year 1987, USNO and RGO, Washington and London, p. K11.

V.N. BARANOV et al., 1986, Space Geodesy, Nedra, Moskva (in Russian), p. 27.

N . BARTELME, P . MEISSL, 1975, Ein einfaches, rasches und numerisch stabiles Verfahren zur Bestimmung des kürzesten Abstandes eines Punktes von einen sphäroidischen Rotationsellipsoid, AVN, Heft 12, p. 436–439 (Wichmann, Karlsruhe) (in German).

W. BENNING, 1974, Der Kürzeste Abstand eines in rechtwinkeligen Koordinaten gegebenen Aussenpunktes vom Ellipsoid, AVN, Heft 11, (Wichmann, Karlsruhe) (in German).

K.M. BORKOWSKI, 1987, Transformation of Geocentric to Geodetic Coordinates without Approximations, Astrophys. Space Sci., 139, 1–4.

D.R. HEDGLEY, 1976, An Exact Transformation from Geocentric to Geodetic Coordinates for Nonzero Altitudes, NASA Techn. Rep. R–458.

M. HEIKKINEN, 1982, Geschlossene Formeln zur Berechnung raumlicher geodätischer Koordinaten aus rechtwinkligen Koordinaten, Zeitschrift Vermess. (in German), 107, 207–211 .

W.A. HEISKANEN and H. MORITZ, 1967, Physical Geodesy, W.H. Freeman and Co., San Francisco, p. 181.

R.HLIBOWICKI (ed.), 1981, Higher Geodesy and Geodetic Astronomy, PWN, Warszawa (in Polish), p. 273.

A.A. IZOTOV et al., 1974, Fundamentals of Satellite Geodesy, Nedra, Moskva (in Russian), p. 34.

G.A. KORN and T.M. KORN, 1961, Mathematical Handbook for Scientists and Engineers, McGraw–Hill, New York, p. 24.

S.A.T. LONG, 1975, General Altitude Transformation between Geocentric and Geodetic Coordinates, Celestial Mech., 12, 225–230.

J. MORRISON and S. PINES, 1961, The Reduction from Geocentric to Geodetic Coordinates, Astron. J., 66, 15–16.

M.K. PAUL, 1973, A Note on Computation of Geodetic Coordinates from Geocentric (Cartesian) Coordinates, Bull. Géod., Nouv. Ser., No 108, 135–139.

M . PICK, 1985, Closed Formulae for Transformation of the Cartesian Coordinate System into a System of Geodetic Coordinates, Studia geoph. et geod., 29, 112–119.

L.G. TAFF, 1985, Celestial Mechanics. A Computational Guide for the Practitioner, J . Wiley and Sons, New York, ch. 3.

M.S. URMAEV, 1981, Orbital Methods of Space Geodesy, Nedra, Moskva (in Russian), p . 14.

P . VANIČEK and E.J. KRAKIWSKI, 1982, Geodesy: The Concepts, North Holland Publ. Co., Amsterdam, p . 324.

E.W. WOOLARD and G.M. CLEMENCE, 1966, Spherical Astronomy, Academic Press, New York, ch. 3.

   Received : 01.06.1988.
   Accepted : 23.09.1988.