Implementation

An implementation in the Java language of the algorithms explained in this paper is available in source form as a stand-alone file at http://www.spaceroots.org/documents/distance/Ellipsoid.java. It is included in a class modelling an ellipsoid.

The language-independant form of this implementation follows. In order to improve readability, some code optimizations have been removed (use of precomputed values like $(1-f)$, $(1-f)^{2}(r^2-a^2)+z^2$ and others). The control parameters of the algorithm are the maximal number of iterations and the two thresholds $\varepsilon$ and $\varepsilon_c$.

if 
$(\sqrt{r^2+z^2}<\varepsilon a_e)$ return 
$\varphi=\pi/2,\quad h=-(1-f)a_e$




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)$



// distance to the ellipse along the current line as the root of a 2nd degree

// 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))$



if 
$(\vert k\vert<\varepsilon\sqrt{r^2+z^2})$ return 
$\varphi,\quad h=k$



begin loop for at most N iterations



// 
$4^$th degree normalized polynom describing circle/ellipse intersections

// 
$\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$



// reduce the polynom to degree 3 by removing the already known real root

// 
$\tau^3 + b \tau^2 + c \tau + d = 0$


$b \leftarrow b + t$


$c \leftarrow c + t$


$d \leftarrow d + t$



// 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$

if $(D >= 0)$ then

// beware that some programming languages do not provide a cubic root function,

// and that the power function does not accept negative arguments

// (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)$

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)$

if 
$(\tilde{\varphi} \times \varphi < 0)$ 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)$

if 
$(\tilde{\varphi} \times \varphi < 0)$ 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)$

end if

end if

end if



// $\varphi$-wise midpoint on the ellipse


$\Delta\varphi \leftarrow \vert\tilde{\varphi} - \varphi\vert/2$


$\varphi \leftarrow (\varphi + \tilde{\varphi})/2$



if 
$(\Delta\varphi < \varepsilon_c)$ return 
$\varphi,\quad h=r\cos\varphi+z\sin\varphi-a_e\sqrt{1-f(2-f)\sin^2\varphi}$




$\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}$

if 
$(\mathrm{inside})$ 
$k \leftarrow -k$


$t \leftarrow \Delta z/(\Delta r+k)$



end loop



generate a convergence error