# Computational most efficient way to convert Cartesian to Geodetic Coordinates

As far as I understand the literature there are several ways to convert a set of geoecentric cartesian coordinates to geodetic coordinates

What would be the most efficient way to convert a cartesian coordinates to geodetic coordinates - most efficient meaning fastest (and most direct) when being implemented on a computer?

-
Do you need to perform conversions for arbitrary geodetic heights or only at nominal ground surface (zero heights)? – whuber Jun 28 '12 at 14:00
Ground surface will do for now, but I'm also interested in arbitrary geodetic heights in the long run. – oschrenk Jun 28 '12 at 14:15

Depends on your accuracy requirements. One iteration of Heiskanen is efficient, but it takes 3 iterations to match the accuracy of a single iteration of Bowring's 1985 method. Even Lin and Wang's algorithm cannot match the efficiency of Bowring, if the trigonometric optimizations included below are used. Therefore, for all around performance I would recommend Bowring's 1985 algorithm.

-- Based on B. R. Bowring's 1985 algorithm (single iteration)
-----------------------------------------------------------------------------
function ECEF_to_Geo(ECEF  : Vector3;
Datum : Datum_Id_Type) return Geographic_Type is
Spheroid : Spheroid_Type := EARTH.DATUMS.Spheroid(Datum);
x    : Meters   renames ECEF(1);
y    : Meters   renames ECEF(2);
z    : Meters   renames ECEF(3);
a    : Meters   renames Models(Spheroid).Semi_Major_Axis;
b    : Meters   renames Models(Spheroid).Semi_Minor_Axis;
e2   : Unitless renames Models(Spheroid).Eccentricity_Squared;
eb2  : Unitless renames Models(Spheroid).Second_Eccen_Squared;
a2   : Scalar := a**2;
d    : Meters := eb2*b;
Ome2 : Unitless := 1.0 - e2;
p2, p, r, tu, tu2, su3, cu3, Phi, Lam, tp, cp, sp, h : Scalar;
Alt : Meters;
Geo : Geographic_Type;
begin -- ECEF_to_Geo
if ((x = 0.0) and (y = 0.0)) then
r   := abs(z);
Phi := Scalar'copy_sign(Half_Pi, z);
Lam := 0.0;
h   := r - b;
elsif (z = 0.0) then
Phi := 0.0;
Lam := arctan(y, x);
p   := sqrt(x**2 + y**2);
h   := p - a;
else
p2  := x**2 + y**2;
p   := sqrt(p2);
r   := sqrt(p2 + z**2);
tu  := b*z*(1.0 + d/r)/(a*p);
tu2 := tu**2;
cu3 := (1.0/sqrt(1.0 + tu2))**3;
su3 := cu3*tu2*tu;
tp  := (z + d*su3)/(p - e2*a*cu3);
Phi := arctan(tp);
Lam := arctan(y, x);
cp  := 1.0/sqrt(1.0 + tp**2);
sp  := cp*tp;
h   := p*cp + z*sp - a*sqrt(1.0 - e2*sp**2);
end if;
Alt := h;
Geo := Make_Geo(Lat, Lon, Alt, Datum);
return Geo;
end ECEF_to_Geo;
-