Initialization

When the $A$ point is not at the center of the ellipse, we initialize the algorithm using a first line containing both the test point $A$ and the center of the ellipse. The cosine and sine of this line slope $\zeta_0$ are defined by equation (6). Note that we do not compute $\zeta_0$ itself, to avoid numerical problems when $r\ll z$ (i.e. when we are close to the minor axis), we also compute $t_0=\tan{\zeta_0/2}=\sin\zeta_0/(1-\cos\zeta_0)$ directly from $\cos\zeta_0$ and $\sin\zeta_0$ (using a stable formula when $\cos\zeta_0>=0$) as we will need it soon.

$$\left\{\begin{aligned}\cos\zeta_0 &= \frac{r}{\sqrt{r^2+z^2}}\\... ...^2+z^2}}\\ t_0 &= \frac{z}{r+\sqrt{r^2+z^2}} \end{aligned}\right.$$

Regardless of its slope, this line is known to have exactly two intersections with the ellipse. We choose the closest one to the $A$ point, which corresponds to the smallest root of binomial equation (5) in absolute value.

Since $1-f>0$, $\cos\zeta_0>0$ and $z \sin\zeta_0\ge0$, we can deduce $a>0$ and $b>0$. This implies that when the binomial has real roots, their sum is strictly positive, so the smallest root in absolute value is also the smallest root in algebric value.

The classical expression for the smallest root is:

$$k_0 = \frac{b-\sqrt{b^2-ac}}{a}$$

However, this expression is not numerically accurate when $b^2 \gg ac$ because in involves computing something similar to $b-(b-\varepsilon)$, leading to a numerical cancellation when $\varepsilon$ is too small. This case occurs when the test point $A$ is very close to the ellipse $\mathcal{E}$. The limit case is for a test points lying exactly on the ellipse, for which $k=0$ is a root so $c=0$.

The cancellation can be avoided even at the limit case by using the dual expression of the root:

$$k_0 = \frac{c}{b+\sqrt{b^2-ac}}$$ (7)

Since we know $b>0$, we know this expression is more stable than the other one in all cases.

The first point $P(\varphi_0)$ is computed from $\zeta_0$ and $k_0$ using equation (4) and the corresponding value $\varphi_0$ is computed from equation (2).

If the test point is very close to the ellipse, we have $\vert k_0\vert<\varepsilon\sqrt{r^2+z^2}$, for a given threshold ratio $\varepsilon$. In fact, in this case any slope would have given the same point since $P(\varphi_0)=A$. There is no need to perform any iteration and we can directly return the values $\varphi=\varphi_0$ and $h=k_0$.