ECEF coordinates to geodetic coordinates: Where is error in calculation?

Question

I know, there are several questions to this topic, but I am just searching for the error in my own equations. I have no deeper geodetic background, but I am trying to understand these formulas instead just copy them.

So, I have a XYZ coordinate and want to get a geodetic coordinate. This is what I tried:

tan(longitude) = Y/X

For the latitude I tried this approach. On a sphere you could calculate it that way:

diameter_xy = sqrt(y^2 + x^2); // => d
tan(latitude) = z / d

Now, I have no sphere but an ellipsoid. So the latitude I already calculated is the same as the geocentric latitude of an ellipsoid, right? (let's call it p)

I tried to transform this into the geodetic latitude (q). The book "Geometric Geodesy I" from R. Rapp (p. 25, formula 3.62) creates this relation between geocentric and geodetic latitudes:

tan(p) = (1 - e^2) * tan(q)

So:

tan(q) = tan(p) / (1 - e^2) 
       = (z / d) / (1 - e^2)                // taking the spherical formula above
       = (z / sqrt(y^2 + x^2)) / (1 - e^2)
       = z / ((1 - e^2) * sqrt(y^2 + x^2))

Well, now I don't come to an expected result. When I have a look a several solutions, the formula differs:

For example: https://gis.stackexchange.com/a/292635 The initial equation for the latitude seems to be the same as mine (beside the point, that atan2 is used instead of atan, which is also not that clear to me). But later the height h becomes part of the equation. It absolutely not clear to my, why.

Well, I don't understand were is error and where do the real solutions (e.g. Ferrari or so) get that complicated.

Edit:

In my case the values result in a greater height. The XYZ point is not on the surface. But since I am not interested in the height above ground, but only in the angles, I thought, this wouldn't make any difference. Maybe I am wrong? Does my calculation only work for h == 0m above ground? If so, how do I recognize the height and how do I eliminate this?

My use case is:

X = 7000000m
Y = 2000000m
Z = 10000000m
e = 0.00669438 // WGS84

So, I am getting:

lat: 53.95°
lon: 15.95°

The lat real value is something like:

lat: 54.04° (h: 6,005 km)

Answer 1

About the ATAN2 function:

Your definitions for lambda (longitude) and psi (geocentric latitude) have a domain restriction:

lambda = ATAN(y/x) , with x > 0. (i.e., just for -PI/2 < lambda < PI/2.)

Because: if x = 0 lambda must be PI/2 or -PI/2 (depends on the sign of y), and for x = -1, y = -1 lambda must be -3*PI/4, not PI/4. It is similar to x = -1, y = 1. With latitudes, the only restriction is at the poles (x = 0 , y = 0).

That domain restriction is solved with ATAN2 function.

---

About the conversion from, following the book nomenclature, psi (geocentric latitude) to phi (geodetic latitude):

Imagine a point P in the y=0 plane, geocentric and geodetic latitudes over the WGS84 ellipsoid may be represented as follows:

Note the following problem: The geocentric latitude must be calculated for P', the orthogonal projection of P over the ellipsoid surface. But the angle obtained by ATAN(z/x) is for the dashed line instead.

The formula provided by R. Rapp in that part of the book is for points that belong to the ellipsoid surface, and can't be used with the coordinates of points far away from it.

---

How can you know if the point belongs to the ellipsoid surface? Verify if your coordinates satisface the oblated ellipsoid equation:

x^2/a^2 + y^2/a^2 + z^2/b^2 = 1

How can you eliminate this? Well, perform the "complicated" calculation.