*
Astrophysics and Space Science*, **139** (1987),
1–4

[Erratum: vol. **146**, (No. 1, July 1988), 201]

COORDINATES WITHOUT APPROXIMATIONS

*
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,

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

and

z =
(Cb^{2}/a^{2} + h) sin φ, | (1b) |

C =
a^{2}/√[(a cos φ)^{2} + (b sin φ)^{2}] | (2) |

r – z/tg φ = (a^{2} –
b^{2})/√[a^{2} + (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 – t^{2})/(2bt) | (4) |

transforms Equation (3) into the following quartic equation

t^{4} + 2Et^{3} + 2Ft – 1 = 0, | (5) |

E = [bz – (a^{2} – b^{2})]/(ar) | (6) |

F = [bz + (a^{2} – b^{2})]/(ar). | (7) |

v^{3} + 3Pv + 2Q = 0, | (8) |

P = (4/3)(EF + 1) | (9) |

Q = 2(E^{2} – F^{2}), | (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) |

D = P^{3} + Q^{2}. | (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 = ±√[G^{2}
+ (F – vG)/(2G – E)] – G, | (13) |

G = [±√(E^{2} + 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 – t^{2})/(2bt)], | (15a) |

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**

**References**

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.]