The following paper has been published in:
Bull. Géod., 63 (1989), pp. 50–56.
Torun Radio Astronomy Observatory,
Nicolaus Copernicus
University,
ul. Chopina 12/18, PL–87–100 Torun, Poland
ACCURATE ALGORITHMS TO TRANSFORM
GEOCENTRIC TO GEODETIC COORDINATES
Abstract 
The problem
of the transformation is reduced to solving of the equation

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 roundingoff
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 – e^{2}) + h] sinφ,  (3b) 
where N = a/(1 – e^{2} sin^{2}φ)^{1/2}, and e^{2} = (a^{2} – b^{2})/a^{2}. 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^{2} – b^{2})/[(ar)^{2} + (bz)^{2}]^{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 builtin 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 fourthdegree polynomial (biquadratic or quartic equation):
t^{4} + 2Et^{3} + 2Ft – 1 = 0,  (9) 
where
E = [bz – (a^{2} – b^{2})]/(ar) and  (10) 
F = [bz + (a^{2} – b^{2})]/(ar).  (11) 
Applying the Ferrari's solution (e.g. Korn and Korn, 1961) to Eq. (9) one obtains:
t = ±√{G^{2} + (F – vG)/(2G – E)} – G,  (12) 
where
G = [±(E^{2} + v)^{1/2} + E]/2,  (13) 
v = (D^{1/2} – Q)^{1/3} – (D^{1/2} + Q)^{1/3},  (14a) 
D = P^{3} + Q^{2},  (15) 
P = (4/3)(EF + 1) and  (16) 
Q = 2(E^{2} – F^{2}).  (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 – t^{2})/(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 zcoordinate 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 zaxis,
where the quantities E and F are obviously indeterminate) without
encountering any problems. As to the roundoff 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^{2} < 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 roundoff 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.
ADDENDUM to the paper The following FORTRAN code is included here for convenience of those interested in practical implementation of the exact algorithm presented. Be aware however that it does not belong to the original article — it comes from a Polish version of the above paper, where you can find also a figure (Rys. 1.) illustrating the geometry involved. 
subroutine GEOD(r,z,fi,h) c Program to transform Cartesian to geodetic coordinates c Input: r, z = equatorial [m] and polar [m] components c Output: fi, h = geodetic coord's (latitude [rad], height [m]) implicit real*8(ah,oz) c IAU ellipsoid: semimajor axis and inverse flattening data a,frec /6378140.d0,298.257d0/ b=dsign(aa/frec,z) if(r.eq.0d0) return E=((z+b)*b/aa)/r F=((zb)*b/a+a)/r P=(E*F+1.)*4d0/3. Q=(E*EF*F)*2. D=P*P*P+Q*Q if(D.ge.0d0) then s=dsqrt(D)+Q s=dsign(dexp(dlog(dabs(s))/3d0),s) v=P/ss v=(Q+Q+v*v*v)/(3*P) else v=2.*dsqrt(P)*dcos(dacos(Q/P/dsqrt(P))/3.) endif G=.5*(E+dsqrt(E*E+v)) t=dsqrt(G*G+(Fv*G)/(G+GE))G fi=datan((1.t*t)*a/(2*b*t)) h=(ra*t)*dcos(fi)+(zb)*dsin(fi) end Test data: r z fi h 4000000. 6000000. 0.985526645027216 847786.688189974 4000.000 6000.000 1.48883906081174 6350591.52477262 