\documentclass[a4paper,11pt,leqno]{article}
\usepackage{color}
\usepackage{graphics}
\usepackage{amsmath}
\usepackage{times}
\usepackage{mathptmx}
\usepackage{url}
\begin{document}
\title{Quick computation of the distance between a point and an ellipse}
\author{L. \textsc{Maisonobe}\thanks{luc@spaceroots.org}}
\date{September 2003, revised May 2005, February 2006, June 2006}
\maketitle
\tableofcontents
\clearpage

\section{Introduction}
We consider the following 2D problem: given a test point $A$ on a
plane and an ellipse $\mathcal{E}$, find the point of the ellipse
which is the closest to the test point.

This problem occurs in several contexts. The first one is to
geopositioning. Planet bodies are often roughly modeled as biaxial
ellipsoids, i.e. ellipsoids having a rotational symetry axis (the
polar axis). For such shapes, the meridian planes are ellipses. Given
a point position in 3D cartesian coordinates, we want to compute the
longitude, geodetic latitude and altitude. The last two data are
obtained by solving this kind of problem. Another context was
encountered while computing approximations of graphical
representations of 3D circles on a 2D display. In order to build an
error model of a B\'ezier-based approximation of the ellipse, the
distance of thousands of approximated curves had to be computed. This
later problem is described in another technical note which is
available on the same place as this one.

The latest version of this document is always available in the
SpaceRoots site downloads page at
\url{http://www.spaceroots.org/downloads.html}. It can be browsed
on-line or retrieved as a PDF, compressed PostScript or LaTeX source
file.

\section{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
\emph{is} used in the following equations.

\begin{figure}[htbp]
\begin{center}\begin{minipage}{100mm}
\caption{\label{fig:definition}coordinates definition}
\setlength{\unitlength}{0.1mm}\begin{picture}(1000,600)
  \textcolor[gray]{0.5}{
    \put(0,50){\vector(1,0){1000}}
    \put(800,20){\mbox{$a_e$}}
    \put(950,60){\mbox{$r$}}
    \put(160,0){\vector(0,1){600}}
    \put(0,400){\mbox{$(1-f)a_e$}}
    \put(170,550){\mbox{$z$}}
  }%

  \textcolor[rgb]{0.15,0.13,0.73}{
    \put(800,150){\mbox{$\mathcal{E}$}}
    \qbezier[0](800,50)(800,230)(572,344)
    \qbezier[0](572,344)(460,400)(160,400)
  }%

  \textcolor[rgb]{0.21,0.39,0.54}{
    \put(572,344){\line(-1,-2){170}}
    \put(606,412){\vector(1,2){34}}
    \put(606,412){\vector(-1,-2){34}}
    \put(575,415){\mbox{$h$}}
    \put(472,394){\line(2,-1){200}}
    \put(590,380){\line(2,-1){36}}
    \put(608,326){\line(1,2){18}}
    \put(640,480){\circle{15}}
    \put(660,480){\mbox{$A$}}
    \put(480,85){\mbox{$\varphi$}}
    \qbezier(475,50)(475,81)(447,94)
  }%
\end{picture}\end{minipage}\end{center}
\end{figure}

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:
\begin{equation}\label{eqn:r-z-expression}
\left\{\begin{aligned}
r &= \left(\frac{a_e}{\sqrt{1-f(2-f)\sin^2\varphi}}+h\right)\cos\varphi\\
z &= \left(\frac{(1-f)^2a_e}{\sqrt{1-f(2-f)\sin^2\varphi}}+h\right)\sin\varphi
\end{aligned}\right.
\end{equation}

In the geopositioning field, $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 geodetic latitude.

The equation~(\ref{eqn:r-z-expression}) 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 when $h$ is known
to be null:
\begin{equation}\label{eqn:phi-for-null-h}
h = 0 \Rightarrow \varphi = \mathrm{atan2}\left(z,r(1-f)^2\right)
\end{equation}

\section{Preliminary results}
\subsection{Special case handling: center or the ellipse}
A special case we will discard in the algorithm is the ellipse
center. This point is easily detected in a preliminary check using a
test like $\sqrt{r^2+z^2}<a_e \varepsilon_\text{close}$, where
$\varepsilon_\text{close}$ is an arbitrary close approach threshold
ratio. The ellipse points closest to the center are both endpoints of
the minor axis. We arbitrarily select the point with positive
inclination $\varphi$: $\varphi=+\frac{\pi}{2}$, $h=-(1-f)a_e$.

\subsection{Inside/outside check}
Given a test point $A$, it is very simple to check if it lies inside
the ellipse or outside of it. This check is done by computing the sign
of expression:
\begin{equation}\label{eqn:inside-outside}
\left(\frac{r}{a_e}\right)^2+\left(\frac{z}{(1-f)a_e}\right)^2 - 1
\end{equation}
the point is outside if the sign is positive, inside otherwise.

\subsection{Intersection between a line and an ellipse}
Lets consider a line having slope $\zeta$ and containing the test
point $A$.

\begin{figure}[htbp]
\begin{center}\begin{minipage}{100mm}
\caption{\label{fig:intersection}line/ellipse intersection}
\setlength{\unitlength}{0.1mm}\begin{picture}(1000,600)
  \textcolor[gray]{0.5}{
    \put(0,50){\vector(1,0){1000}}
    \put(950,60){\mbox{$r$}}
    \put(160,0){\vector(0,1){600}}
    \put(170,550){\mbox{$z$}}
  }%

  \textcolor[rgb]{0.15,0.13,0.73}{
    \qbezier[0](800,50)(800,230)(572,344)
    \qbezier[0](572,344)(460,400)(160,400)
  }%

  \textcolor[rgb]{0.00,0.39,0.00}{
    \put(672,394){\vector(2,1){100}}
    \put(672,394){\vector(-2,-1){99}}
    \put(662,404){\mbox{$k$}}
    \put(572,344){\line(-2,-1){200}}
    \put(372,259){\line(1,0){140}}
    \put(480,270){\mbox{$\zeta$}}
    \qbezier(452,259)(452,271)(447,281)
    \put(772,444){\circle{15}}
    \put(740,460){\mbox{$A$}}
  }%

  \textcolor[rgb]{0.21,0.39,0.54}{
    \put(400,0){\line(1,2){185}}
    \put(492,384){\line(2,-1){160}}
    \put(572,344){\circle{15}}
    \put(530,380){\mbox{$P(\varphi)$}}
    \put(480,85){\mbox{$\varphi$}}
    \qbezier(475,50)(475,81)(447,94)
  }%
\end{picture}\end{minipage}\end{center}
\end{figure}

For well chosen $\zeta$ values, this line intersects the ellipse. Let
$(r',z')$ be the cartesian coordinates of the intersection point
$P(\varphi)$:
\begin{equation}\label{eqn:intersection-point}
\left\{\begin{aligned}
r' &= r-k\cos\zeta\\
z' &= z-k\sin\zeta
\end{aligned}\right.
\end{equation}

where $k$ is the signed distance between the test point $A$ and the
intersection point $P(\varphi)$. $k$ is positive if the test point is
outside of the ellipse and negative otherwise. Note that the sign of
$k$ implies the choice of the line orientation, so depending on the
test point location inside or outside of the ellipse, we have to
choose either $\zeta$ or $\pi-\zeta$; the following equations take
care of this choice. Since the intersection point $P(\varphi)$ belongs
to the ellipse, the following equation holds:
\begin{equation*}
\left(\frac{r'}{a_e}\right)^2+\left(\frac{z'}{(1-f)a_e}\right)^2 = 1
\end{equation*}
introducing the test point coordinates and the signed distance $k$
\begin{align}
&\left(\frac{r-k\cos\zeta}{a_e}\right)^2+\left(\frac{z-k\sin\zeta}{(1-f)a_e}\right)^2
= 1\nonumber\\
\Leftrightarrow\quad
&(1-f)^2[(r-k\cos\zeta)^2-a_e^2] + (z-k\sin\zeta)^2 = 0\nonumber\\
\Leftrightarrow\quad
&a k^2 - 2b k + c = 0\label{eqn:k-zeta-relation}
\quad\text{where}\quad
\left\{\begin{aligned}
a &= (1-f)^2\cos^2\zeta+\sin^2\zeta\\
b &= (1-f)^2r\cos\zeta+z\sin\zeta\\
c &= (1-f)^2(r^2-a_e^2)+z^2
\end{aligned}\right.
\end{align}

\section{Iterative method}
We intend to compute $h$ and $\varphi$ from $r$ and $z$ using a very
quick iterative method based on a suite of lines containing the test
point $A$ and intersecting the ellipse on varying points
$P(\varphi_n)$. The lines will converge to the projection line.

\subsection{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 (\ref{eqn:initialization}). 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\ge0$) as we will need it soon.

\begin{equation}\label{eqn:initialization}
\left\{\begin{aligned}
\cos\zeta_0 &= \frac{r}{\sqrt{r^2+z^2}}\\
\sin\zeta_0 &= \frac{z}{\sqrt{r^2+z^2}}\\
t_0 \quad   &= \frac{z}{r+\sqrt{r^2+z^2}}
\end{aligned}\right.
\end{equation}

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~(\ref{eqn:k-zeta-relation}) in absolute value.

Since $1-f>0$, and both $r \cos\zeta_0$ and $z \sin\zeta_0$ are
non-negative and cannot be simultaneously null, we can deduce $a>0$
and $b>0$. Hence since the binomial has real roots (the line does
intersect the ellipse), their sum is strictly positive, therefore the
smallest root in absolute value is also the smallest root in algebraic
value.

The classical expression for the smallest root in algebraic value is:
\begin{equation*}
k_0 = \frac{b-\sqrt{b^2-ac}}{a}
\end{equation*}

However, this expression is not numerically accurate when $b^2 \gg ac$
because it 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 is avoided even at the limit case by using the dual
expression of the root which is more stable than the classical one
since $b>0$:
\begin{equation}\label{eqn:k-root}
k_0 = \frac{c}{b+\sqrt{b^2-ac}}
\end{equation}

The first point $P(\varphi_0)$ is computed from $\zeta_0$ and $k_0$
using equation~(\ref{eqn:intersection-point}) and the corresponding
value $\varphi_0$ is computed from equation~(\ref{eqn:phi-for-null-h}).

If the test point is very close to the ellipse, we have
$|k_0|<\sqrt{r^2+z^2}\varepsilon_\text{close}$, using the same
arbitrary close approach threshold ratio $\varepsilon_\text{close}$
already used for handling the special case of ellipse center. 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$.

\subsection{Iterations}
At the start of each iteration, we have an ellipse point
$P(\varphi_n)$ and the corresponding inclination $\varphi_n$, the
signed distance $k_n\ne0$ between $P(\varphi_n)$ and the test point
$A$ and the value of $t_n=\tan{\zeta_n/2}$ where $\zeta_n$ is the
slope of the line joining $P(\varphi_n)$ to the test point.

If the estimate $P(\varphi_n)$ were perfect, we would have
$\varphi_n=\varphi$ and $k_n=h$. If $h\ne0$ the circle centered on the
test point $A$ with radius $|k_n|$ would be tangent to the ellipse. In
our iterative process, the estimate is not perfect and the circle
intersects the ellipse at least at one other point
$P(\tilde{\varphi}_n)$,
$\tilde{\varphi}_n\ne\varphi_n$. Figures~\ref{fig:interior-point}
and~\ref{fig:exterior-point} show these circle/ellipse intersections
in two different situations. As shown by these figures, there can be
up to three intersection points in addition to the already known point
$P(\varphi_n)$.

\begin{figure}[htbp]
\begin{center}\begin{minipage}{60mm}
\caption{\label{fig:interior-point}interior point}
\setlength{\unitlength}{0.1mm}\begin{picture}(350,600)
  \textcolor[gray]{0.5}{
    \put(30,250){\vector(1,0){540}}
    \put(550,260){\mbox{$r$}}
    \put(60,100){\vector(0,1){280}}
    \put(70,390){\mbox{$z$}}
  }%

  \textcolor[rgb]{0.15,0.13,0.73}{
    \qbezier[0](460,250)(460,265.9)(429.6,280.6)
    \qbezier[0](429.6,280.6)(399.1,295.3)(342.8,306.6)
    \qbezier[0](342.8,306.6)(286.6,317.8)(213.1,323.9)
    \qbezier[0](213.1,323.9)(139.6,330)(60,330)
    \qbezier[0](60,170)(139.6,170)(213.1,176.1)
    \qbezier[0](213.1,176.1)(286.6,182.2)(342.8,193.4)
    \qbezier[0](342.8,193.4)(399.1,204.7)(429.6,219.4)
    \qbezier[0](429.6,219.4)(460,234.1)(460,250)
  }%

  \textcolor[rgb]{0.00,0.39,0.00}{
    \put(300,265){\line(6,1){160}}
    \put(300,265){\circle{15}}
    \put(290,265){\mbox{$A$}}
    \put(419.8,284.9){\circle{15}}
    \put(420,310){\mbox{$P(\varphi_n)$}}
    \put(194.6,325.3){\circle{15}}
    \put(110,350){\mbox{$P(\tilde{\varphi}_n)$}}
    \put(216.9,176.4){\circle{15}}
    \put(408.7,210.8){\circle{15}}
    \qbezier[0](421.4,265)(421.4,289.2)(412.2,311.5)
    \qbezier[0](412.2,311.5)(402.9,333.8)(385.9,350.9)
    \qbezier[0](385.9,350.9)(368.8,367.9)(346.5,377.2)
    \qbezier[0](346.5,377.2)(324.2,386.4)(300,386.4)
    \qbezier[0](300,386.4)(275.8,386.4)(253.5,377.2)
    \qbezier[0](253.5,377.2)(231.2,367.9)(214.1,350.9)
    \qbezier[0](214.1,350.9)(197.1,333.8)(187.8,311.5)
    \qbezier[0](187.8,311.5)(178.6,289.2)(178.6,265)
    \qbezier[0](178.6,265)(178.6,240.8)(187.8,218.5)
    \qbezier[0](187.8,218.5)(197.1,196.2)(214.1,179.1)
    \qbezier[0](214.1,179.1)(231.2,162.1)(253.5,152.8)
    \qbezier[0](253.5,152.8)(275.8,143.6)(300,143.6)
    \qbezier[0](300,143.6)(324.2,143.6)(346.5,152.8)
    \qbezier[0](346.5,152.8)(368.8,162.1)(385.9,179.1)
    \qbezier[0](385.9,179.1)(402.9,196.2)(412.2,218.5)
    \qbezier[0](412.2,218.5)(421.4,240.8)(421.4,265)
  }%
\end{picture}\end{minipage}\hspace{3mm}\begin{minipage}{60mm}
\caption{\label{fig:exterior-point}exterior point}
\setlength{\unitlength}{0.1mm}\begin{picture}(350,600)
  \textcolor[gray]{0.5}{
    \put(30,250){\vector(1,0){540}}
    \put(550,260){\mbox{$r$}}
    \put(60,100){\vector(0,1){280}}
    \put(70,390){\mbox{$z$}}
  }%

  \textcolor[rgb]{0.15,0.13,0.73}{
    \qbezier[0](60,170)(139.6,170)(213.1,176.1)
    \qbezier[0](213.1,176.1)(286.6,182.2)(342.8,193.4)
    \qbezier[0](342.8,193.4)(399.1,204.7)(429.6,219.4)
    \qbezier[0](429.6,219.4)(460,234.1)(460,250)
    \qbezier[0](460,250)(460,265.9)(429.6,280.6)
    \qbezier[0](429.6,280.6)(399.1,295.3)(342.8,306.6)
    \qbezier[0](342.8,306.6)(286.6,317.8)(213.1,323.9)
    \qbezier[0](213.1,323.9)(139.6,330)(60,330)
  }%

  \textcolor[rgb]{0.00,0.39,0.00}{
    \put(560,50){\line(-5,2){350}}
    \put(560,50){\circle{15}}
    \put(570,60){\mbox{$A$}}
    \put(238.9,178.4){\circle{15}}
    \put(200,120){\mbox{$P(\varphi_n)$}}
    \put(332.1,310.1){\circle{15}}
    \put(260,340){\mbox{$P(\tilde{\varphi}_n)$}}
    \qbezier[0](361.6,333.3)(314.2,300.1)(280.2,253.3)
    \qbezier[0](280.2,253.3)(246.2,206.5)(229.3,151.1)
  }%
\end{picture}\end{minipage}\end{center}
\end{figure}

The additional intersection points can be computed by introducing
$\tau=\tan{\zeta/2}$ in equation (\ref{eqn:k-zeta-relation}) an
multiplying by $(1+\tau^2)^2$. We get:
\begin{align*}
&\begin{aligned}[t]
k^2 &[(1+\tau^2)^2-f(2-f)(1-\tau^2)^2]\\
-2k &[(1-f)^2r(1+\tau^2)(1-\tau^2)+2z\tau(1+\tau^2)]\\
+   &[(1-f)^2(r^2-a_e^2)+z^2](1+\tau^2)^2 = 0
\end{aligned}\\
\Leftrightarrow&
\phantom{+}\alpha\tau^4+\beta\tau^3+\gamma\tau^2+\delta\tau+\varepsilon
= 0\\
\text{where}\quad&
\left\{\begin{aligned}
\alpha &= (1-f)^2((r+k)^2-a_e^2)+z^2\\
\beta  &= -4kz\\
\gamma &= 2[(1-f)^2(r^2-k^2-a_e^2)+2k^2+z^2]\\
\delta &= -4kz\\
\varepsilon &=(1-f)^2((r-k)^2-a_e^2)+z^2
\end{aligned}\right.
\end{align*}

This is a quartic equation for which we already know one root:
$t_n=\tan{\zeta_n/2}$.  We can therefore reduce this equation to a
normalized cubic polynomial by a simple coefficients identification
(the value of the constant term $\varepsilon$ of the quartic is of
course not needed for this degree reduction):
\begin{equation}\label{eqn:cubic}
 \tau^3+a_1\tau^2+a_2\tau+a_3=0
\quad\text{where}\quad
\left\{\begin{aligned}
a_1 &= \frac{\beta}{\alpha} + t_n\\
a_2 &= \frac{\gamma}{\alpha} + a_1 t_n\\
a_3 &= \frac{\delta}{\alpha} + a_2 t_n
\end{aligned}\right.
\end{equation}

The other circle/ellipse intersections are given by the real roots of
equation~(\ref{eqn:cubic}). In order to compute them, we compute first
the discriminant $D$:
\begin{equation*}
Q=\frac{3a_2-a_1^2}{9}
\qquad
R=\frac{9a_1a_2-27a_3-2a_1^3}{54}
\qquad
D=Q^3+R^2
\end{equation*}

If $D$ is positive, the following expressions compute the two real
numbers $S$ et $T$ and allow to deduce the unique real root
$\tilde{t}_{n,1}$. We do not need the two complex roots
$\tilde{t}_{n,2}$ and $\tilde{t}_{n,3}$. This corresponds to the case
depicted in figure~\ref{fig:exterior-point}. It also occurs for
interior points when we are close to the solution.

\begin{gather}
S=\sqrt[3]{R+\sqrt{Q^3+R^2}}
\qquad
T=\sqrt[3]{R-\sqrt{Q^3+R^2}}\nonumber\\
\left\{\begin{aligned}
\tilde{t}_{n,1} &= S+T-\frac{a_1}{3}\\
\tilde{t}_{n,2} &= -\frac{S+T}{2}-\frac{a_1}{3}+i\frac{\sqrt{3}(S-T)}{2}\\
\tilde{t}_{n,3} &= -\frac{S+T}{2}-\frac{a_1}{3}-i\frac{\sqrt{3}(S-T)}{2}\\
\end{aligned}\right.
\end{gather}

If $D$ is negative, the three roots are all real ones. This
corresponds to the case depicted in
figure~\ref{fig:interior-point}. In this case, the trigonometric
expressions of the root are simpler to use:
\begin{gather}
\left\{\begin{aligned}
\tilde{t}_{n,1} &= 2\sqrt{-Q}\cos{\frac{\theta}{3}}-\frac{a_1}{3}\\
\tilde{t}_{n,2} &= 2\sqrt{-Q}\cos{\frac{\theta+2\pi}{3}}-\frac{a_1}{3}\\
\tilde{t}_{n,3} &= 2\sqrt{-Q}\cos{\frac{\theta+4\pi}{3}}-\frac{a_1}{3}\\
\end{aligned}\right.
\quad\text{with}\quad
\cos\theta = \frac{R}{\sqrt{-Q^3}}
\end{gather}

For each $\tilde{t}_{n,i}$ representing a circle/ellipse intersection,
we compute $\tilde{\zeta}_{n,i}$ and $P(\tilde{\varphi}_{n,i})$ using
equation~(\ref{eqn:intersection-point}). Since this point belongs to
the ellipse we deduce $\tilde{\varphi}_{n,i}$ using
equation~(\ref{eqn:phi-for-null-h}).

When there are three roots, only one of them\footnote{a proper
demonstration would be welcome here ...} leads to a point on the same
half plane ($z<0$ or $z>0$) as the test point $A$ and the already
known intersection point $P(\varphi_n)$. We select this root and
disgard the two other ones. This selection is done by comparing the
signs of $\varphi_n$ and $\tilde{\varphi}_{n,i}$.

\subsection{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 $P(\varphi_{n+1})$ ensures $|k_{n+1}|<|k_n|$.

As can be noted from figure~\ref{fig:exterior-point}, when we are far
from the solution, the elliptical arc is not restricted to the angular
sector bounded by the two lines $\left(A,P(\varphi_n)\right)$ with
slope $\zeta_n$ and $\left(A,P(\tilde{\varphi}_n)\right)$ with slope
$\tilde{\zeta}_n$. Hence we cannot use a simple interpolation between
$\zeta_n$ and $\tilde{\zeta}_n$. We can however interpolate between
$\varphi_n$ and $\tilde{\varphi}_n$.

Interpolating with all the available information ($\zeta_n$,
$\varphi_n$, $\tilde{\zeta}_n$ and $\tilde{\varphi}_n$) might be
tempting at first sight. However, some tests show that in corner cases
with extremely flat ellipses (say $f = 0.99$), computing
$\varphi_{n+1}$ by weighting $\varphi_n$ and $\tilde{\varphi}_n$ with
$|\tilde{\varphi}_n-\tilde{\zeta}_n|$ and $|\varphi_n-\zeta_n|$
respectively fails to converge. Even when convergence can be reached,
the accuracy gain with respect to the following simpler method is
almost impossible to spot.

A safer and simpler solution is based on the following observation:
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. 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:
\begin{equation}
\varphi_{n+1} = \frac{\varphi_{n}+\tilde{\varphi}_n}{2}
\end{equation}

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:
\begin{equation}\left\{\begin{aligned}
\Delta r = r-\frac{a_e\cos\varphi_{n+1}}{\sqrt{1-f(2-f)\sin^2\varphi_{n+1}}}\\
\Delta z = z-\frac{(1-f)^2a_e\sin\varphi_{n+1}}{\sqrt{1-f(2-f)\sin^2\varphi_{n+1}}}
\end{aligned}\right.\end{equation}
and to deduce the signed distance and the tangent of the half slope:
\begin{equation}\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.\end{equation}
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~(\ref{eqn:inside-outside}). 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.

\subsection{Convergence}
Since at the end of each iteration we have an interval known to
contain the solution, and since we choose the middle of this interval
($\varphi$-wise), we can use its half-size as a convergence
criterion. This bound is very pessimistic. In fact it corresponds
almost to the error of the previous iteration as it is based on the
previous $\varphi$.

We stop the iterations when:
\begin{equation}
\frac{|\varphi_{n}-\tilde{\varphi}_n|}{2}<\varepsilon_\text{angle}
\end{equation}
where $\varepsilon_\text{angle}$ is the angular convergence threshold.

When this threshold is met, we use $\varphi=\varphi_{n+1}$ and we
compute $h$ using the following expression which is true when
$\varphi$ has been found:
\begin{equation}\label{eqn:h}
h = r\cos\varphi+z\sin\varphi-a_e\sqrt{1 - f(2-f)\sin^2\varphi}
\end{equation}

\section{Order of the method}
The numerical example in appendix shows that the error decreases very
quickly. Some experiments have been performed with a large number of
significant digits (100 digits). These experiments seem to show the
method would be of order 2. This order is not very
impressive. However, the main strength of the method is that it starts
very quickly, even for very flat ellipses.  The numerical example in
appendix shows that in a quite difficult case we already have
$0.013^\circ$ after only one iteration (but we don't know it, as the
upper bound of the error we can compute during this first iteration is
much larger).

This order should be analyzed more precisely.

\section{Robustness}
The method also seems quite robust. The cancellation problem when
computing $k_0$ is avoided by using the proper expression. Another
cancellation problem occurs while computing $R-\sqrt{Q^3+ R^2}$ when $R
\gg Q$. However, this has no consequence as it affects only the $T$
term which in this case is much smaller than $S$ and only $S+T$ is
used to compute the solution. As the algorithm converges, all the
coefficients also converge ($k$, $\alpha$, $\beta$, $\gamma$,
$\delta$, $a_1$, $a_2$, $a_3$, $Q$, $R$, $S$, $T$). There is no
numerical explosion anywhere.

\appendix
\clearpage\section{Approximate formulas}
The technical manual for \emph{Geocentric Datum of Australia} provides
direct transformation
formulas\footnote{\url{http://www.anzlic.org.au/icsm/gdatm/xyzcd.htm},
but this URL seems not valid anymore} that we reproduce here, after
having renamed some parameters to be consistent with our notation and
to correct an ambiguity\footnote{$f$ was used in the document for both
the flatness and the latitude}, and after having replaced an
equality sign by an approximation sign:
\begin{equation}
\tan{\varphi} \approx \frac{z(1-f) + f(2-f) a_e sin^3u}
                           {(1-f)(r - f(2-f) a_e \cos^3u)}
\end{equation}
where
\begin{equation}
\tan u = \frac{z}{r}\left[(1-f) + f(2-f) \frac{a_e}{\sqrt{r^2+z^2}}\right] 
\end{equation}

Some numerical checks show that these expressions are very good
approximation but that they can be applied only for nearly spherical
bodies. For the Earth, the maximal error are about $10^{-11}$ degrees
in latitude and $2\times 10^{-9}$ metres in altitude!

However, when we consider very flat ellipses (for example $f=0.9$),
the errors become tremendous (more than $100$ degrees on $\varphi$).

\clearpage\section{Numerical example}
The following example shows the first iterations of the algorithm for a
very flat ellipse:
\begin{align*}\left\{\begin{aligned}
a_e &= 100\\
f   &= 0.9
\end{aligned}\right.
\end{align*}

We consider a test point $A$ having the following coordinates:
\begin{equation*}A\left\{\begin{aligned}
\varphi &= 75^\circ\\
h       &= 0.1
\end{aligned}\right.
\Rightarrow
\left\{\begin{aligned}
r &= 93.713969911344535171\\
z &= \phantom{0}3.593079627683806165
\end{aligned}\right.\end{equation*}

We assume we know only $r$ and $z$ and want to compute the unknown
$\varphi$ and $h$.

Equations (\ref{eqn:initialization}) give us the tangent of the half
slope:
\begin{equation*}
t_0=0.0191634
\end{equation*}

Applying equation (\ref{eqn:k-root}), we get the signed distance
between $A$ and the intersection of the first line and the ellipse:
\begin{equation*}
k_0= 0.3419763
\Rightarrow
\varphi_0=75.3818944^\circ
\end{equation*}

From this value of $k$, we deduce the first coefficients $\alpha$,
$\beta$, $\gamma$ and $\delta$ of the quartic (we don't compute the
uneeded constant term $\varepsilon$) and reduce it to a cubic using
the already known root $t_0$:
\begin{equation*}
 \tau^3 + a_1 \tau^2 + a_2 \tau + a_3=0
\quad\text{where}\quad
\left\{\begin{aligned}
a_1 &= -3.5542558\\
a_2 &= \phantom{-}1.3365805\\
a_3 &= -3.5478057\\
\end{aligned}\right.
\end{equation*}
This polynomial has only one real root:
\begin{equation*}
\tilde{t}_0=3.4640705
\Rightarrow
\tilde{\varphi}_0=74.5916410^\circ
\end{equation*}

The middle point between $\varphi_0$ and $\tilde{\varphi}_0$ is the
best approximation we have at the end of this first iteration:
\begin{equation*}
\varphi_1=74.9867677^\circ
\end{equation*}

We know an upper bound for the error of this approximation is the
half-width of the interval, which is $0.3951267^\circ$ (in fact the
error is about 30 times smaller). As we consider this is too much, we
compute the starting values for a second iteration:
\begin{equation*}
k_1=0.1005980
\quad
t_1=0.8577741
\end{equation*}

The steps of the second iteration are similar to the steps of the
first one. They lead to the following result:
\begin{equation*}
\tilde{\varphi}_1=75.0132026^\circ
\end{equation*}

So at the end of the second iteration, we have an estimate as the
middle point between $\varphi_1$ and $\tilde{\varphi}_1$:
$\varphi_2=74.9999851^\circ$ and an upper bound of its error of
$0.0132175^\circ$ (the real error being almost 900 times smaller).

A third iteration would lead to $\varphi_3=74.99999999998141^\circ$
and an upper bound of its error of $1.48\times 10^{-5}{}^\circ$ (the
real error being almost 800000 times smaller).

\clearpage\section{Implementation}

This appendix presents a pseudo-code implementation of the algorithm
described in this paper. In order to improve readability, there are no
code optimizations. The control parameters of the algorithm are the
maximal number of iterations and the two thresholds
$\varepsilon_\text{close}$ and $\varepsilon_\text{angle}$.

An implementation of this pseudo-code in the Java language is
available in source form as a stand-alone file on the SpaceRoots
site\footnote{\url{http://www.spaceroots.org/documents/distance/Ellipsoid.java}}.
It is designed as a class modelling an ellipsoid. This implementation
includes some code optimizations (use of precomputed values like
$(1-f)$, $(1-f)^{2}(r^2-a^2)+z^2$ and others). The control parameters
have been arbitrarily set to $\varepsilon_\text{close} = 10^{-10}$,
$\varepsilon_\text{angle}=10^{-14}$ radians and $N = 100$.
\begin{tabbing}
xx\=xx\=xx\=xx\=xx\=xx\=\kill
\textbf{if} $(\sqrt{r^2+z^2}<a_e \varepsilon_\text{close})$ \textbf{then}\+\\
\textbf{return} $\varphi=\pi/2,\quad h=-(1-f)a_e$\-\\
\textbf{end if}\\
\\
$\text{inside} \leftarrow (1-f)^{2}(r^2-a^2)+z^2\le 0$\\
$\cos\zeta\leftarrow r/\sqrt{r^2+z^2}$\\
$\sin\zeta\leftarrow z/\sqrt{r^2+z^2}$\\
$t\leftarrow z/\left(r+\sqrt{r^2+z^2}\right)$\\
\\
\textsl{// distance to the ellipse along the current line as the root of a 2nd degree}\\
\textsl{// polynom closest to zero: $a k^2 - 2 b k + c = 0$}\\
$a \leftarrow (1-f)^2\cos^2\zeta+\sin^2\zeta$\\
$b \leftarrow (1-f)^2r\cos\zeta+z\sin\zeta$\\
$c \leftarrow (1-f)^2(r^2-a_e^2)+z^2$\\
$k \leftarrow c / (b + \sqrt{b^2 - a c})$\\
$\varphi \leftarrow \mathrm{atan2}(z - k\sin{\zeta}, (1-f)^2 (r - k\cos\zeta))$\\
\\
\textbf{if} $(|k|<\sqrt{r^2+z^2}\varepsilon_\text{close})$ \textbf{then}\+\\
  \textbf{return} $\varphi,\quad h=k$\-\\
\textbf{end if}\\
\\
\textbf{begin loop for at most N iterations}\+\\
\\
\textsl{// $4^\text{th}$ degree normalized polynom describing circle/ellipse intersections}\\
\textsl{// $\tau^4 + b \tau^3 + c \tau^2 + d \tau + e = 0$ (there is no need to compute e here)}\\
$a \leftarrow (1-f)^2((r+k)^2-a_e^2)+z^2$\\
$b \leftarrow -4kz / a$\\
$c \leftarrow 2\left[(1-f)^2(r^2-k^2-a_e^2)+2k^2+z^2\right] / a$\\
$d \leftarrow b$\\
\\
\textsl{// reduce the polynom to degree 3 by removing the already known real root}\\
\textsl{// $\tau^3 + b \tau^2 + c \tau + d = 0$}\\
$b \leftarrow b + t$\\
$c \leftarrow c + t b$\\
$d \leftarrow d + t c$\\
\\
\textsl{// find the other real root}\\
$Q \leftarrow (3 c - b^2) / 9$\\
$R \leftarrow (b (9 c - 2 b^2) - 27 d) / 54$\\
$D \leftarrow Q^3+R^2$\\
\textbf{if} $(D >= 0)$ \textbf{then}\+\\
  \textsl{// beware that some programming languages do not provide a cubic root function,}\\
  \textsl{// and that the power function does not accept negative arguments}\\
  \textsl{// (both $R+\sqrt{D}$ and $R-\sqrt{D}$ can be negative)}\\
  $\tilde{t} \leftarrow \sqrt[3]{R+\sqrt{D}} + \sqrt[3]{R-\sqrt{D}} - b/3$\\
  $\tilde{\varphi} \leftarrow \mathrm{atan2}\left(z(1+\tilde{t}^2)-2k\tilde{t},
                                                 (1-f)^2 [r(1+\tilde{t}^2) - k(1-\tilde{t}^2)]\right)$\-\\
\textbf{else}\+\\
  $Q \leftarrow -Q$\\
  $\theta \leftarrow \arccos{R/(Q\sqrt{Q})}$\\
  $\tilde{t} \leftarrow 2\sqrt{Q} \cos(\theta/3) - b/3$\\
  $\tilde{\varphi} \leftarrow \mathrm{atan2}\left(z(1+\tilde{t}^2)-2k\tilde{t}, (1-f)^2 [r(1+\tilde{t}^2) - k(1-\tilde{t}^2)]\right)$\\
  \textbf{if} $(\tilde{\varphi} \times \varphi < 0)$ \textbf{then}\+\\
    $\tilde{t} \leftarrow 2\sqrt{Q} \cos((\theta+2\pi)/3) - b/3$\\
    $\tilde{\varphi} \leftarrow \mathrm{atan2}\left(z(1+\tilde{t}^2)-2k\tilde{t}, (1-f)^2 [r(1+\tilde{t}^2) - k(1-\tilde{t}^2)]\right)$\\
    \textbf{if} $(\tilde{\varphi} \times \varphi < 0)$ \textbf{then}\+\\
      $\tilde{t} \leftarrow 2\sqrt{Q} \cos((\theta+4\pi)/3) - b/3$\\
      $\tilde{\varphi} \leftarrow \mathrm{atan2}\left(z(1+\tilde{t}^2)-2k\tilde{t}, (1-f)^2 [r(1+\tilde{t}^2) - k(1-\tilde{t}^2)]\right)$\-\\
  \textbf{end if}\-\\
\textbf{end if}\-\\
\textbf{end if}\\
\\
\textsl{// $\varphi$-wise midpoint on the ellipse}\\
$\Delta\varphi \leftarrow |\tilde{\varphi} - \varphi|/2$\\
$\varphi \leftarrow (\varphi + \tilde{\varphi})/2$\\
\\
\textbf{if} $(\Delta\varphi < \varepsilon_\text{angle})$ \textbf{then}\+\\
\textbf{return} $\varphi,\quad h=r\cos\varphi+z\sin\varphi-a_e\sqrt{1-f(2-f)\sin^2\varphi}$\-\\
\textbf{end if}\\
\\
$\Delta_r \leftarrow r - \frac{a_e\cos\varphi}{\sqrt{1-f(2-f)\sin^2\varphi}}$\\
$\Delta_z \leftarrow z - \frac{a_e(1-f)^2\sin\varphi}{\sqrt{1-f(2-f)\sin^2\varphi}}$\\
$k \leftarrow \sqrt{\Delta_r^2+\Delta_z^2}$\\
\textbf{if} $(\mathrm{inside})$ \textbf{then}\+\\
  $k \leftarrow -k$\-\\
\textbf{end if}\\
$t \leftarrow \Delta z/(\Delta r+k)$\-\\
\\
\textbf{end loop}\\
\\
\textbf{generate a convergence error}\\
\end{tabbing}
\end{document}