Selecting a new point for next iteration

We now have two points $P(\varphi_n)$ and $P(\tilde{\varphi}_n)$ known to be both at signed distance $k_n$ from the test point $A$. The ellipse points between $P(\varphi_n)$ and $P(\tilde{\varphi}_n)$ are all closer to $A$ than $P(\varphi_n)$ and $P(\tilde{\varphi}_n)$ (they are inside the circle). So selecting one point in this elliptical arc as the next point ensures $\vert k_{n+1}\vert<\vert k_n\vert$.

As can be noted from figure 4, when we are far from the solution, the elliptical arc is not restricted to the angular sector bounded by the two lines $\left(T,P(\varphi_n)\right)$ and $\left(T,P(\tilde{\varphi}_n)\right)$, so we cannot use a simple interpolation between $\zeta_n$ and $\tilde{\zeta}_n$. We can however interpolate between $\varphi_n$ and $\tilde{\varphi}_n$.

Once we are close to the solution, the circle is almost tangent to the ellipse, so the elliptical arc between the two intersection point is small and has an almost constant radius of curvature: it is very close to a circle. In this case, $P(\tilde{\varphi}_n)$ is close to the symetrical of $P(\varphi_n)$ with respect to the true projection. We rely on this asymptotical behavior to set up the selection scheme for the next iteration inclination:

$$\varphi_{n+1} = \frac{\varphi_{n}+\tilde{\varphi}_n}{2}$$ (11)

In order to start the new iteration, we also need to compute $k_{n+1}$ and $t_{n+1}$ from $\varphi_{n+1}$. The simplest way to do this is to compute first the vector joining the new intersection point $P(\varphi_{n+1})$ to the test point:

$$\left\{\begin{aligned}\Delta r = r-\frac{a_e\cos\varphi_{n+1}}{... ...i_{n+1}}{\sqrt{1-f(2-f)\sin^2\varphi_{n+1}}} \end{aligned}\right.$$

and to deduce the signed distance and the tangent of the half slope:

$$\left\{\begin{aligned}k_{n+1} &= \pm\sqrt{\Delta r^2+\Delta z^2... ...t_{n+1} &= \frac{\Delta z}{\Delta r+k_{n+1}} \end{aligned}\right.$$

The sign of $k_{n+1}$ is chosen to be negative when the test point is inside the ellipse and positive otherwise, which is checked using equation (3). Here again, we do not compute $\zeta_{n+1}$ but directly $t_{n+1}$ from the vector coordinates, again in order to avoid numerical problems for points close to the minor axis ( $\zeta\approx\pm\pi/2$) where the arc tangent function accuracy is very bad in typical computers.