With respect to what "normal" is geodetic latitude defined for the WGS84 elliptical model of the Earth?

There are some features of geodetic latitude specifically with respect the WGS84 that I do not quite understand. In short:

With respect to what "normal" is geodetic latitude defined for the WGS84 elliptical model of the Earth?

I'll ask the question again in a little more detail at the bottom after giving some context and a few examples.

In earth centered coordinates, the WGS84 ellipse is defined by the equation

``````(x/r)^2 + (y/r)^2 + (z/br)^2 = 1
``````

where `r ~= 6378137` meters is the radius of earth and `b = (1-f)r` where `f ~= 3.35e-3` is the flattening factor. According to the internet geodetic latitude is:

The angle between the normal and the equatorial plane. This is the definition assumed when the word latitude is used without qualification. The definition must be accompanied with a specification of the ellipsoid.

This seems sensible enough. The normal, that is, the vector perpendicular to the tangent plane of the WGS84 ellipse at the earth centered, earth fixed (ECEF) point `(x,y,z)` is

``````N(x,y,z) = ( x/r, y/r, z/br )
``````

It follows that should I want a point on the ellipse at a given geodetic latitude `l` and longitude `g`, I first find a normal vector pointing in this direction, e.g,

``````n = ( cos(l) cos(g), cos(l) * sin(g), sin(l) )
``````

and then map it to a corresponding point on the ellipse

``````p =  ( r n[1], r n[2], r b n[3] )
``````

where `n[1]`, `n[2]`, `n[3]` are the coordinates of the normal vector `n`, e.g. `n[3] = sin(l)`. One can check that the normal `N(p)` is in the direction `n` and so does point in the desired direction described by the latitude `l` and longitude `g`. Specially the angle between the normal `n` and the equatorial plane is `l` as desired.

To confirm this consider the following example. Suppose we want to know the ECEF point corresponding with latitude 10 degrees and longitude 0 degrees. We construct a normal of the form

``````n = ( cos(10*pi/180), 0, sin(10*pi/180) ) = (0.984807753012208, 0, 0.1736481776669303)
``````

The corresponding ECEF point with this normal is

``````( r cos(10*pi/180), 0, rb sin(10*pi/180) ) = (6281238.767374026,0,1103838.45524868)
``````

We can confirm that this point is on the WGS84 ellipse with

``````(6281238.767374026/r)^2 + (1103838.45524868/rb)^2 = 0.9999999999999999
``````

which seems close enough to 1. One can confirm that the normal of this point is in the direction of n as required and that the angle between the normal and the equatorial plane is

``````atan(1103838.45524868/6281238.767374026) = 0.1745329251994329 radians = 9.999999999999997 degrees
``````

which is close enough to the desired 10 degrees.

This all seems sensible and simple. However, it is not what is computed by standard, typically iterative, algorithms (e.g. here) for determining location on the WGS84 ellipse in ECEF coordinates given geodetic latitude and longitude. Online calculators exist too (e.g. here) that all appear use equivalent algorithms, iterative or otherwise, and produce consistent results. For example, when asked for the ECEF location corresponding with latitude 10 degrees and longitude 0 degrees these algorithms return

``````( 6281872.83, 0.0, 1100248.55 )
``````

You can reproduce this calculation at here. One can confirm that this point is indeed on the WGS84 ellipse

``````(6281872.83/r)^2 + (1100248.55/br)^2 = 1.000000000245795
``````

which seems close enough to 1. The normal at this point is

``````( 6281872.83/r, 0.0, 1100248.55/br ) = (0.9849071648978377, 0, 0.1719247524673192)
``````

The angle between this normal and the equatorial plane is

``````atan(0.1719247524673192/0.9849071648978377) = 0.1739595032375784 radian = 9.967145341705625 degree
``````

This is not 10 degrees and so it appears that this point is not at a geodetic latitude of 10 degrees.

This is rather confusing and leads to the more detailed version of my question.

Question:

If the typical definition of the normal `N(x,y,z)` as the vector perpendicular to the tangent plane is not what is used to define geodetic latitude for WGS84, then what "normal" is used? Is there an intuitive description of this normal? For example, is it related to the gravitational normal, or some other entity with a physical or geometric meaning?

Note:

It's possible that I've made an error, either conceptual, arithmetical, or otherwise, and that the latitude is defined using the normal perpendicular to the tangent plane once that error is corrected. In that case identification of the error would be the accepted answer to this question.

``````atan(1103838.45524868/6281238.767374026)
``````atan(1103838.45524868*r^2/(6281238.767374026*(br)^2))