Problem description

We will use the traditional notations used in geodesy to describe the problem. We assume the case of longitude $\lambda$ has already been solved using $\lambda=\mathrm{atan2}(y,x)$ and that the remaining problem has been restricted to the 2D meridian plane using the coordinates $r=\sqrt{x^2+y^2}$ and $z$.

We consider a test point $A$ and an ellipse $\mathcal{E}$. After a suitable coordinates change, we can consider the coordinates of the test point in the canonical reference frame of the ellipse. This frame is centered on the ellipse center, has its absissae axis along the major axis and has its ordinates axis along the minor axis. In this reference frame, the coordinates of $A$ are $(r,z)$. Using suitable axes orientations, we can arrange to have $r>=0$. This assumption is used in the following equations.

Figure 1: coordinates definition

The ellipse is defined by its semi-major axis is $a_e$, and either its flattening $f$ ($0\le f<1$) or its semi-minor axis $a_p=a_e(1-f)$.

The signed distance between the point and the ellipse $h$ and the inclination $\varphi$ of the projection of the point on the ellipse are related to the cartesian coordinates:

$$\left\{\begin{aligned}r &= \left(\frac{a_e}{\sqrt{1-f(2-f)\sin^... ...{1-f(2-f)\sin^2\varphi}}+h\right)\sin\varphi \end{aligned}\right.$$

In the geopositioning domain, $a_e$ is the equatorial radius of the body, $a_p$ is the polar radius, $h$ is the altitude above the ellipsoid (negative when the point is below the surface of the ellipsoid) and $\varphi$ is the geodesic latitude.

The equation (1) is easy to apply when $h$ and $\varphi$ are known and $r$ and $z$ are desired, but it is impossible to reverse in the general case. It can be reversed in the specific case where $h$ is known to be null:

$$h = 0 \Rightarrow \varphi = \mathrm{atan2}\left(z,r(1-f)^2\right)$$ (2)