# Minimax polynomial ECEF to geodetic coordinate transformation approximations

John-Olof Nilsson  The author is with Avioniq Sweden AB, e-mail<EMAIL_ADDRESS>
###### Abstract.

Minimax polynomial ECEF to geodetic coordinate transformation approximations
are presented, including often preferable n-vector versions. The
approximations provide tunable computational-cost-to-accuracy trade-off and an
unprecedented low latency down to an accuracy of $\sim 10^{-5}$m, which is
demonstrated in an extensive benchmark. This sets a new standard for fast ECEF
to geodetic coordinate transformations and opens up a new realm of further
improvement opportunities and extensions to other geodetic quantities.

##  1\. Introduction

Cartesian earth-centered earth-fixed (ECEF) coordinate $\\{x,y,z\\}$ to
geodetic coordinates $\\{\phi,\lambda,h\\}$ (latitude, longitude and altitude,
_lla_) or $\\{\bar{n},h\\}$ (n-vector and altitude, _nva_) transformations are
integral to geodesy and related modeling and simulations and, for many
applications, constitute a significant computational cost. The significance of
the problem is demonstrated by well over 200 publications spanning more than 6
decades [67]. However, with its apparent simplicity and many publications, it
may also appear to be a closed chapter, but, somewhat surprisingly, methods
based on minimax polynomials, the cornerstone of _approximation theory_ , are
completely lacking. This report fills in this blank and demonstrates how the
_complete_ transformations can be implemented with minimax polynomials,
including often preferable n-vector versions. (See [32] for the benefits of
the n-vector representation.) This results in high accuracy-to-computational-
cost tunability, straight-forward arithmetic-only separable implementations
and unprecedented low latency down to an accuracy of at least $\sim 10^{-5}$m,
at which the Bowring’s method is equally good. The performance is demonstrated
in an extensive benchmark of more than 40 methods. This sets a new standard
for fast ECEF to geodetic coordinate transformations (and evaluations) and
opens up a new realm of further improvement and extension opportunities.

The improvement opportunities primarily comes from the resemblance to
elementary function approximations, for which a vast literature and tool sets
exist. For a starting point, see [64]. See also the _potential improvement_ in
the discussion Section 9. The extension opportunities primarily comes from the
presented approximation method being trivially extendable to other (in
addition to latitude and altitude) ellipsoidal $z$-axis rotation symmetric and
$xy$-plan reflection symmetric or anti-symmetric quantities. The one caveat to
it all is that the approximations are only valid for a _preselected_ altitude
range. The range can be large enough to cover the vast majority of
applications with performance figures provided for up to the altitude range
$[-5000,500\;000]$m. However, significantly beyond that, one may have to
resort to range splitting, which, on the other hand, is normal for minimax
approximations. Further, the presented benchmark is valuable in itself. It is
the most extensive benchmark to date and the first benchmark including _nva_
transformations, which arguably provides a superior view of the performance of
different methods. The results indicate that, in contrast to frequent claims
of the opposite; in terms of computational efficiency and for altitude ranges
and conversion accuracies reasonable for most applications; there have only
been a few marginal improvements beyond [9] \+ [10] method, with [68], [30]
and [79] being noteworthy. Similar conclusions have also been reached in a
smaller benchmark [17].

If you just want a fast polynomial approximation, skip to Table 1 for pseudo-
code implementations, then skip to the benchmark results in Figure 1 to select
polynomial orders and, finally, get the polynomials from Appendix G. Or take
the long route and start in Section 2 for some history and continue with
Section 3 where the enabling initial problem transformation is presented.
Section 4 and 5 follows with polynomial approximations of the transformation
components and the complete transformations. Section 6 deals with the
implementation, Section 7 with the benchmark and Section 8 discusses the
benchmark results. Finally, general discussions, caveats, potential
improvements and final recommendations are provided in Section 9. Many details
are referred to the Appendices A-G.

##  2\. Background and Novelty

As for most minimax approximations, there are four fundamental components of
the presented method; _range reduction_ , _series expansion_ , _polynomial
approximation_ and _range reconstruction_ ; for which some background and
general comments are provided here. First a general comment, there is a
significant Soviet Bloc literature on the coordinate transformation subject.
This literature is hard to access and has largely been ignored but references
are listed in [67].

The _range reduction_ aim at reducing the essentially unbounded ECEF
coordinate values to some bounded interval. For polynomial methods, this
essentially entail transforming the ECEF coordinates to the variables of the
series expansion and the polynomial approximation. Similar steps are present
in most iterative and close-form methods. Here, the range reduction is split
into a commonly appearing, but here implicit, geocentric latitude part for the
derivation and an explicit but significantly less expensive coordinate
normalization for the actual approximation.

_Series expansions_ are not strictly necessary for polynomial approximation
but often helpful. Polynomial approximation, and particularly series
expansions, of ECEF (geocentric latitude) to geodetic coordinate
transformations is nothing new. Rather it appears to have been the normal way
to transform coordinates in the pre-digital-computers era. Expansions can be
found on p40-43 in the standard reference [43]. Helmert, in turn, cite even
older literature. Further, an extensive treatment of series expansions can be
found in [1]. With the launch of Sputnik 1 on the 4th of October 1957 and the
heightened interest and requirements of the space race, numerous publications
of altitude dependent series expansions can be found from 1958 and up to 1976
[7, 45, 62, 33, 46, 74, 61, 73, 34, 59, 60, 19, 85]. However, after the
seminal [9] iterative method, they appear to wane in popularity. Some notable
exceptions are [22, 68, 93, 38, 23]. There has also recently been some work on
non-altitude dependent expansions [69, 56, 51]. However, the series expansions
of the literature typically have three fundamental short-comings: 1) They
partially rely on inefficient (local) derivative based series expansions. 2)
The series variables are often not directly available. 3) The series are only
for latitude. Here novel general series expansions are presented, directly in
ECEF coordinate norm and normalized coordinates for latitude, altitude and the
n-vector. The series expansions contains no derivatives.

The _polynomial approximation_ of the literature methods are, without
exceptions, done by truncating the series expansions. This is convenient but
it is well known in the area of _approximation theory_ that this is typically
suboptimal. Rather minimax approximations should normally be used. The Remez
algorithm for computing such polynomials has been around since [76] and the
field can be described as well established since the 60s [28]. However, there
is a simple reason minimax approximations has not been used so far, the
coordinate transformation is not a 1-dimensional problem, making the
techniques not directly applicable. Nonetheless, here it is shown how minimax
approximations can be used. Only series which can be truncated with close to
minimax performance are truncated. In all other places, optimal minimax
approximations are used.

Finally, the _range reconstruction_ aim at transforming some intermediate
approximated quantity to the final quantities of interest. For most
approximation methods this means transforming some quantity, proportional to
trigonometric function of the latitude, to the final latitude and altitude,
i.e. the range reconstruction come in the form of inverse trigonometric
functions. These inverse trigonometric functions are just as much a part of
the coordinate transformation as any other part and should be approximated
with suitable accuracy. Most likely, this is occasionally done in practice but
with some rare exceptions [88, 89, 103] this is not discussed, despite the
range reconstruction making up half the computational cost (as we will see)
and suitable approximations being readily available since at least the 50s,
see [39]. Further, in many cases, n-vector representation is preferable [32],
and they can be though of as just a different range reconstruction. Most
iterative and closed form methods are straight-forward to convert. However,
most polynomial methods are not, in a sensible way. Here matched inverse
trigonometric function approximations are provided and dedicated (minimax)
approximations for the n-vector representation.

Altogether, the mentioned (novel) components make polynomial approximations
competitive again. Further, the computer development have probably contributed
here too. Modern super-scalar vector CPUs are very good at evaluating
polynomials, giving polynomial methods and edge over iterative methods.

Note that iterative methods appear to have been around for a long time as
well. The classical method by Hirvonen was first published in [44]111I have
not yet managed to access the frequently cited 1958 publication by K. Rinner,
or found any description of its results. Hence, it is not cited just yet..
[62] states _“Various procedures for obtaining a solution of these equations
by iteration techniques exist.”_ without any further references. In turn,
(published) exact solutions has been around since at least [21, 84] with more
following soon after [87, 72, 6, 40]. (Also [21] states that the
transformation is normally carried out by iterative methods.) However, [7]
refers to solutions by _“Lagrange multiplier”_ and _“high-degree algebraic
equation”_. Further, [97] claims the earliest closed form solution can be
found in [20] but, as pointed out by [4], the latitude equation can be found
e.g. on p49 in the book _Analytical Conics_ [83]. Finally, methods for solving
quartic equations are known since Ferrari’s 1545 solution, so exact solutions
have probably been applied before Ecker’s and Sugai’s work.

##  3\. Problem transformation

The multivariate and essentially unbounded nature of the coordinate
transformation makes polynomial approximation techniques not directly
applicable. The unboundedness is handled by the _range reduction_ comprised of
the ECEF to geocentric coordinates transformations

(1) $\begin{matrix}\phi_{c}\\\ \lambda_{c}\\\
h_{c}\end{matrix}\;\begin{matrix}=\\\ =\\\
=\end{matrix}\;\begin{matrix}[l]\sin^{-1}(z/p)\\\ \operatorname{atan2}(y,x)\\\
p-h_{0}\end{matrix}\quad\quad\text{and}\quad\quad\bar{n}_{c}=\begin{bmatrix}x/p\\\
y/p\\\ z/p\end{bmatrix}=\begin{bmatrix}r\\\ s\\\ t\end{bmatrix}$

where $p=\sqrt{x^{2}+y^{2}+z^{2}}$ and where $\phi_{c}$, $\lambda_{c}$,
$h_{c}$, $h_{0}$ and $\bar{n}_{c}$ are the geocentric latitude, longitude,
altitude, reference radius and n-vector, respectively. Note that

(2) $\bar{n}=\begin{bmatrix}\cos(\phi)\cos(\lambda)\\\
\cos(\phi)\sin(\lambda)\\\
\sin(\phi)\end{bmatrix}\quad\text{and}\quad\bar{n}_{c}=\begin{bmatrix}\cos(\phi_{c})\cos(\lambda_{c})\\\
\cos(\phi_{c})\sin(\lambda_{c})\\\ \sin(\phi_{c})\end{bmatrix}$

and that $\lambda_{c}=\lambda$. Also note, later on, $h_{0}$ will be taken to
be 0 but $h_{c}$ and $p$ are still kept separate since a non-zero $h_{0}$
could be used, or even some other expression for $h_{c}$. Clearly, $\phi_{c}$,
$\lambda_{c}$, $r$, $s$ and $t$ are bounded and $h_{c}$ is bounded in the
sense that an interval of interest can typically be set. In turn, the
multivariate nature is handled by the following key _series expansions_
derived in Appendix A

(3) $\displaystyle\phi$
$\displaystyle=\phi_{c}+t\sqrt{1-t^{2}}\omega(h_{c},t^{2})\quad\quad$
$\displaystyle\sin(\phi)$
$\displaystyle=\sin(\phi_{c})\eta(t^{2},\omega(h_{c},t^{2}))$ $\displaystyle
h$ $\displaystyle=\mu(h_{c},t^{2})$ $\displaystyle\cos(\phi)$
$\displaystyle=\cos(\phi_{c})\rho(t^{2},\omega(h_{c},t^{2}))$

where, with the customary (abuse of) notation $\sum^{\infty}_{n=0}(\cdot)$
meaning $\lim_{N\to\infty}\sum^{N}_{n=0}(\cdot)$ and $\sum^{\prime}$ meaning
that the zeroth term is to be halved,

$\begin{split}\omega(u,v)&=\sum^{\infty}_{n=1}\sum_{k=0}^{n-1}(-1)^{k}\binom{2n}{2k+1}v^{k}(1-v)^{n-k-1}b_{n}(u)\\\
\mu(u,v)&=u+\mathop{\smash{\sideset{}{{}^{\prime}}{\sum}}\vphantom{\sum}}\limits^{\infty}_{n=0}\sum_{k=0}^{n}(-1)^{k}\binom{2n}{2k}v^{k}(1-v)^{n-k}c_{n}(u)\\\
\eta(v,w)&=\sigma(\delta(v,w))+(1-v)\tau(w,\delta(v,w))\\\
\rho(v,w)&=\sigma(\delta(v,w))-v\tau(w,\delta(v,w))\end{split}$

where $b_{n}(\cdot)$ and $c_{n}(\cdot)$ are altitude dependent Fourier
coefficients, encapsulating the _reference ellipsoid_ properties, and

$\delta(v,w)=v(1-v)w^{2}$
$\sigma(\delta)=\sum_{l\in\mathbb{N}_{e}}(-1)^{l/2}\frac{1}{l!}\delta^{l/2}\quad\text{and}\quad\tau(w,\delta)=w\sum_{l\in\mathbb{N}_{o}}(-1)^{(l-1)/2}\frac{1}{l!}\delta^{(l-1)/2}$

where $\mathbb{N}_{e}$ and $\mathbb{N}_{o}$ are the sets of all positive even
and odd natural numbers. Note,
$\delta(t^{2},\omega(h_{c},t^{2}))=(\phi-\phi_{c})^{2}$, i.e. $\sigma(\delta)$
and $\tau(w,\delta)$ are polynomials in $(\phi-\phi_{c})^{2}$. Further, note,
the only assumption on $\phi$ and $h$ in the derivation is that they are
rotation symmetric with respect to the z-axis and that they are reflection
anti-symmetric and symmetric, respectively, with respect to the $xy$-plane,
i.e. the expressions are valid for any such quantities. The quantity specific
attributes are carried by $b_{n}(\cdot)$ and $c_{n}(\cdot)$.

For $h$ no _range reconstruction_ is required while, from (2)-(3), it follows
for $n$ as

(4) $\bar{n}\approx\begin{bmatrix}r\,\rho(t^{2},\omega_{,}(h_{c},t^{2}))\\\
s\,\rho(t^{2},\omega_{,}(h_{c},t^{2}))\\\
t\,\eta(t^{2},\omega_{,}(h_{c},t^{2}))\end{bmatrix}$

while for $\phi$ thre are six reconstruction alternatives

(5)
$\begin{split}\phi&=\sin^{-1}\left(t\,\eta(t^{2},\omega(h_{c},t^{2}))\right)\\\
\phi&=\sin^{-1}(t)+t\sqrt{1-t^{2}}\omega(h_{c},t^{2})\\\
\phi&=\text{sgn}(z)\cos^{-1}\\!\left(q/p\,\rho(t^{2},\omega(h_{c},t^{2}))\right)\end{split}$

$\begin{split}\phi&=\text{sgn}(z)\cos^{-1}\\!\left(q/p\right)\\!+\\!t\sqrt{1\\!-\\!t^{2}}\omega(h_{c},t^{2})\\\
\phi&=\tan^{-1}(z/q)\\!+\\!t\sqrt{1-t^{2}}\omega(h_{c},t^{2})\\\
\phi&=\tan^{-1}\left(\frac{z\,\eta(t^{2},\omega(h_{c},t^{2}))}{q\rho(t^{2},\omega(h_{c},t^{2}))}\right)\end{split}$

where $q=\sqrt{x^{2}+y^{2}}$. $\sin^{-1}$ is the least expensive to
approximate but is ill-conditioned around the poles (the derivative goes to
infinity) and using it anywhere near the poles requires a square root argument
reduction. $\cos^{-1}$ requires a square root and a divide and is similarly
ill-conditioned around the equator. $\tan^{-1}$ is well conditioned for all
$\phi_{c}$ but requires a square root, a divide and some more arithmetic
operations compared to $\sin^{-1}$ and $\cos^{-1}$. Similarly, the additive
correction is ill-conditioned around the poles (the derivative of
$t\sqrt{1-t^{2}}$ goes to infinity) whereas the multiplicative correction is
well conditioned for all $\phi_{c}$. See Appendix B for some figures of what
this ill-conditioning means.

As will be shown, for the least accurate approximations,
$\sin^{-1}$-expressions are suitably used. Combining the $\sin^{-1}$ and
$\cos^{-1}$ expressions for values around the poles gives potentially good
performance but is a hard combination to tune, among other things to avoid
discontinuous $\phi$ errors, and gives undesirable conditional code, and is
not further investigated here. See Section 9 for some more discussion on it.
For more accurate approximations, the $\tan^{-1}$ expressions are suitably
used.

##  4\. Partial polynomial approximations

The problem transformation of the previous sections reduces the approximation
problem to that of $\omega(u,v)$, $\mu(u,v)$, $\sigma(\delta)$,
$\tau(w,\delta)$, $\sin^{-1}$, $\tan^{-1}$ and $\operatorname{atan2}$.
Naturally, minimax approximation are sought. The challenge is the two
dimensional nature of $\omega(u,v)$ and $\mu(u,v)$. Fortunately, $\mu(u,v)$,
being a Fourier cosine series, is a Chebyshev series with respect to $v$ and
truncated Chebyshev series give close to minimax approximations [11].
Similarly, truncation of $\omega(u,v)$ gives a close to minimax optimal
approximation with respect to the weighting $v\sqrt{1-v^{2}}$. For
$b_{n}(h_{c})$, $c_{n}(h_{c})$ and $\sigma(\delta)$ and $\tau(w,\delta)$,
explicit minimax approximations are used. Let

$\begin{split}\\{b^{(M)}_{n,m}:m\in[0,M]\\}&=\operatorname*{arg\,min}_{\\{b_{m}:m\in[0,M]\\}}\max_{h_{c}\in[{h_{c}}_{\min},{h_{c}}_{\max}]}\left|b_{n}(h_{c})-\sum_{m=0}^{M}b_{m}h_{c}^{m}\right|\\\
\\{c^{(M)}_{n,m}:m\in[0,M]\\}&=\operatorname*{arg\,min}_{\\{c_{m}:m\in[0,M]\\}}\max_{h_{c}\in[{h_{c}}_{\min},{h_{c}}_{\max}]}\left|c_{n}(h_{c})-\sum_{m=0}^{M}c_{m}h_{c}^{m}\right|\\\
\\{\varsigma^{(L)}_{l}:l\in[0,\lfloor{L/2}\rfloor]\\}&=\operatorname*{arg\,min}_{\\{\varsigma_{l}:l\in[0,\lfloor{L/2}\rfloor]\\}}\max_{\delta\in[0,\delta_{\max}]}\left|\sigma(\delta)-\sum_{l=0}^{\lfloor{L/2}\rfloor}\varsigma_{l}\delta^{l}\right|\\\
\\{\vartheta^{(L)}_{l}:l\in[0,\lfloor{(L-1)/2}\rfloor]\\}&=\operatorname*{arg\,min}_{\\{\vartheta_{l}:l\in[0,\lfloor{(L-1)/2}\rfloor]\\}}\max_{\delta\in[0,\delta_{\max}]}\left|\tau(1,\delta)-\sum_{l=0}^{\lfloor{(L-1)/2}\rfloor}\vartheta_{l}\delta^{l}\right|\end{split}$

where $b^{(M)}_{n,m}$ and $c^{(M)}_{n,m}$ are the $m^{th}$ coefficients of the
minimax approximation of degree $M$ with respect to $h_{c}$ of the $n^{th}$
Fourier coefficients. Similarly, $\varsigma^{(L)}_{l}$ and
$\vartheta^{(L)}_{l}$ are the $l^{th}$ coefficients of the minimax
approximation of degree $\lfloor{L/2}\rfloor$ and $\lfloor{(L-1)/2}\rfloor$,
respectively, i.e. a total degree of the multiplicative corrections of $L$.
The ranges $[{h_{c}}_{\min},{h_{c}}_{\max}]$ and $[0,\delta_{\max}]$ are the
ranges over which the polynomial approximations are minimax optimal. See
Section 6 for a discussions on $[{h_{c}}_{\min},{h_{c}}_{\max}]$. In turn,
$\delta_{\max}$ is defined by $[{h_{c}}_{\min},{h_{c}}_{\max}]$ via

$\begin{split}\delta_{\max}&=\max_{\begin{subarray}{c}h_{c}\in[{h_{c}}_{\min},{h_{c}}_{\max}]\\\
\phi_{c}\in[-\pi/2,\pi/2]\end{subarray}}\delta(\sin^{2}(\phi_{c}),\omega(h_{c},\sin^{2}(\phi_{c})))\\\
&=\max_{\begin{subarray}{c}h_{c}\in[{h_{c}}_{\min},{h_{c}}_{\max}]\\\
\phi_{c}\in[0,\pi/2]\end{subarray}}(\phi-\phi_{c})^{2}\end{split}$

Altogether, this gives the (essentially) minimax power series approximations

(6) $\displaystyle\omega(u,v)$ $\displaystyle\approx\omega_{N,M}(u,v)$
$\displaystyle\sigma(\delta)$ $\displaystyle\approx\sigma_{L}(\delta)$
$\displaystyle\mu(u,v)$ $\displaystyle\approx\mu_{N,M}(u,v)$
$\displaystyle\quad\quad\quad\tau(w,\delta)$
$\displaystyle\approx\tau_{L}(w,\delta)$

where

$\begin{split}\omega_{N,M}(u,v)&=\sum^{N}_{n=1}\sum_{k=0}^{n-1}(-1)^{k}\binom{2n}{2k+1}v^{k}(1-v)^{n-k-1}\sum_{m=0}^{M}b^{(N)}_{n,m}u^{m}\\\
\mu_{N,M}(u,v)&=u+\mathop{\smash{\sideset{}{{}^{\prime}}{\sum}}\vphantom{\sum}}\limits^{N}_{n=0}\sum_{k=0}^{n}(-1)^{k}\binom{2n}{2k}v^{k}(1-v)^{n-k}\sum_{m=0}^{M}c^{(M)}_{n,m}u^{m}\\\
\sigma_{L}(\delta)&=\sum_{l=0}^{\lfloor{L/2}\rfloor}\varsigma^{(L)}_{l}\delta^{l}\quad\text{and}\quad\tau_{L}(w,\delta)=w\sum_{l=0}^{\lfloor{(L-1)/2}\rfloor}\vartheta^{(L)}_{l}\delta^{l}\end{split}$

Minimax approximations are also naturally used for $\sin^{-1}$, $\tan^{-1}$
and $\operatorname{atan2}$; math library functions are typically minimax
approximations. Many different approximations are conceivable. _The focus here
is not these approximations_ for which a vast literature and tool sets exist,
see e.g. [64, 65, 18, 12]. However, to exemplify the performance, basic
approximations are used based on the sub-range ($[0,1]$ and $[-1,1]$)
approximations

$\begin{split}\\{\alpha^{(I)}_{i}:i\in[0,I]\\}&=\operatorname*{arg\,min}_{\\{\alpha_{i}:i\in[0,I]\\}}\max_{x\in[0,1]}\left|\sin^{-1}(x)-\tfrac{\pi}{2}+\sqrt{1-x}\sum_{i=0}^{I}\alpha_{i}x^{i}\right|\\\
\\{\beta^{(J)}_{j}:j\in[0,J]\\}&=\operatorname*{arg\,min}_{\\{\beta_{j}:j\in[0,J]\\}}\max_{x\in[-1,1]}\left|\tan^{-1}(x)-\sum_{j=0}^{J}\beta_{j}x^{2j+1}\right|\end{split}$

and the resulting polynomials

$\chi_{I}^{\prime}(x)=\sum^{I}_{i=1}\alpha^{(I)}_{i}x^{i}\quad\text{and}\quad\xi_{J}^{\prime}(x)=\sum^{J}_{j=0}\beta^{(J)}_{j}x^{2j+1}$

These approximations are mapped to the full ranges by the following
conditional-free bit-twiddling transformations

(7)
$\begin{split}\sin^{-1}(x)&\approx\chi_{I}(x)=\texttt{copysign}\left(\tfrac{\pi}{2}-\sqrt{1-|x|}\chi_{I}^{\prime}(|x|),x\right)\\\
\tan^{-1}(y,x)&\approx\psi_{I}(y,x)=\texttt{copysign}\left(\tfrac{\pi}{4}+\xi_{J}^{\prime}\left(\tfrac{|y|-x}{|y|+x}\right),y\right)\quad\forall\quad
x\geq 0\\\ \operatorname{atan2}(y,x)&\approx\xi_{J}(y,x)=\pi
d(y,x)+\texttt{copysign}\left(\tfrac{\pi}{4}+\xi_{J}^{\prime}\left(\tfrac{|y|-|x|}{|y|+|x|}\right),x\oplus
y\right)\end{split}$

where $d(y,x)=(((x\land y)\ggg(\texttt{sizeof}(x)\cdot
8-2))\land-2)\lor(x\gg(\texttt{sizeof}(x)\cdot 8-1))$, and -2 within it, are
Two’s complement integers of the size of $x$ and $y$ implicitly converted to
real in the multiplication with $\pi$ and where
$\texttt{copysign}(\cdot,\cdot)$ and $\texttt{sizeof}(\cdot)$ have the obvious
(C/C++) meanings. Further, all logical operators ($\oplus$ is xor) are bit-
wise and $\ggg$ and $\gg$ are bit shift with and without sign extension,
respectively. Note, $\approx\psi_{I}(0,0)$ contains a $0/0$. Most
$\operatorname{atan2}$ implementations would return $0$ but $\lambda$ is not
defined for the poles so it is acceptable. Note the two argument version of
arctangent $\tan^{-1}(y,x)$. In contrast to $\operatorname{atan2}(y,x)$, the
second argument is limited and the output is in the range $[-\pi/2,\pi/2]$.
The fact that the $\tan^{-1}$ argument can be split in a dividend and a
_positive_ divisor is instead used to make the mapping from $[-1,1]$
conditional free at the cost of two extra additions. Further, note that the
bit manipulations of the real values, which are typically implementation
defined for most programming languages, only rely on the real value
representation having an initial leading sign bit, so it is most likely
possible to implement on almost all platforms but may require some extra care.

##  5\. Complete polynomial approximations

Combining (2), (3), (5), (6), (4) and (7) gives the four latitude
approximations

(8) $\displaystyle\phi$
$\displaystyle\approx\chi_{I}(t)+t\sqrt{1-t^{2}}\,\omega_{N_{\phi},M_{\phi}}(h_{c},t^{2})$
(9) $\displaystyle\phi$
$\displaystyle\approx\psi_{I}(z,q)+t\sqrt{1-t^{2}}\,\omega_{N_{\phi},M_{\phi}}(h_{c},t^{2})$
(10) $\displaystyle\phi$
$\displaystyle\approx\chi_{I}(t\,\eta_{L_{\phi}}(t^{2},\omega_{N_{\phi},M_{\phi}}(h_{c},t^{2})))$
(11) $\displaystyle\phi$
$\displaystyle\approx\psi_{I}\left(z\,\eta_{L_{\phi}}(t^{2},\omega_{N_{\phi},M_{\phi}}(h_{c},t^{2})),q\rho_{L_{\phi}}(t^{2},\omega_{N_{\phi},M_{\phi}}(h_{c},t^{2}))\right)$

and the longitude, n-vector and altitude approximations

(12) $\displaystyle\lambda$ $\displaystyle\approx\xi_{J}(y,x)$ (13)
$\displaystyle\bar{n}$
$\displaystyle\approx\begin{bmatrix}r\,\rho_{L_{n}}(t^{2},\omega_{N_{n},M_{n}}(h_{c},t^{2}))\\\
s\,\rho_{L_{n}}(t^{2},\omega_{N_{n},M_{n}}(h_{c},t^{2}))\\\
t\,\eta_{L_{n}}(t^{2},\omega_{N_{n},M_{n}}(h_{c},t^{2}))\end{bmatrix}$ (14)
$\displaystyle h$ $\displaystyle\approx\mu_{N_{h},M_{h}}(h_{c},t^{2})$

From these follow the four _lla_ ($\\{\phi,\lambda,h\\}$) approximations

(15)
$\begin{split}f_{\sin\\!+,\\{N_{\phi},M_{\phi},N_{h},M_{h},I,J\\}}(\\{x,y,z\\})&=\left\\{\eqref{eq:lat_approx1},\eqref{eq:lon_approx},\eqref{eq:alt_approx}\right\\}\\\
f_{\tan\\!+,\\{N_{\phi},M_{\phi},N_{h},M_{h},I,J\\}}(\\{x,y,z\\})&=\left\\{\eqref{eq:lat_approx2},\eqref{eq:lon_approx},\eqref{eq:alt_approx}\right\\}\\\
f_{\sin\\!\times,\\{L_{\phi},N_{\phi},M_{\phi},N_{h},M_{h},I,J\\}}(\\{x,y,z\\})&=\left\\{\eqref{eq:lat_approx3},\eqref{eq:lon_approx},\eqref{eq:alt_approx}\right\\}\\\
f_{\tan\\!\times,\\{L_{\phi},N_{\phi},M_{\phi},N_{h},M_{h},I,J\\}}(\\{x,y,z\\})&=\left\\{\eqref{eq:lat_approx4},\eqref{eq:lon_approx},\eqref{eq:alt_approx}\right\\}\end{split}$

and the _nva_ ($\\{\bar{n},h\\}$) approximation

(16)
$f_{\bar{n},\\{L_{n},N_{n},M_{n},N_{h},M_{h}\\}}(\\{x,y,z\\})=\left\\{\eqref{eq:n_approx},\eqref{eq:alt_approx}\right\\}$

for finite index limits $I$, $J$, $L_{\phi}$, $L_{n}$, $N_{\phi}$, $M_{\phi}$,
$N_{n}$, $M_{n}$, $N_{h}$ and $M_{h}$.

##  6\. Implementation

Implementing (15) and (16) entail computation of the coefficients, suitable
polynomial transformations and arithmetic sequences for the polynomials and
auxiliary quantities. For an application, a specific approximation with
specific index limits, giving the best performance for a given platform for
the desirable accuracy or computational cost design constraints, also have to
be selected, i.e. a benchmark as exemplified in the next section should be
carried out. (Even though many application will probably get away with just
selecting an approximation from this report.)

Further, the reference radius $h_{0}$ and the range
$[{h_{c}}_{\min},{h_{c}}_{\max}]$ are parameters of the approximations. Given
minimax approximations, $h_{0}$ is of less importance and can be chosen as
$h_{0}=0$, which eliminates an addition for each transformation. If other
polynomial approximations, i.e. Taylor, Padé, etc., are used, $h_{0}$ will be
of importance. $[{h_{c}}_{\min},{h_{c}}_{\max}]$ should cover all $h_{c}$
values of an application. For most applications, it is easier to set an
altitude range $[h_{\min},h_{\max}]$ which gives
$[{h_{c}}_{\min},{h_{c}}_{\max}]=[h_{\min}+b-h_{0},h_{\max}+a-h_{0}]$, where
$a$ and $b$ are the reference ellipsoid semi-major and semi-minor axes. Note,
even a zero altitude range still requires a geocentric altitude range of
$a-b\approx 21385$m. For the benchmark, the limits $h_{\min}=-5000$m and
$h_{\max}=100\,000$m are primarily used but results for larger and smaller
ranges are also provided in Appendix D. Given minimax approximations, $h_{0}$
is of less importance and can be chosen as $h_{0}=0$, which eliminates an
addition for each transformation.

The computations of the coefficients $\varsigma^{(L)}_{l}$,
$\vartheta^{(L)}_{l}$, $\alpha^{(I)}_{i}$ and $\beta^{(J)}_{i}$ only depend on
arithmetic and elementary functions. Hence, they were trivially computed with
the Remez exchange algorithm implementation of Sollya 8.0 [14]. In contrast,
the computations of the coefficients $b^{(M)}_{n,m}$ and $c^{(M)}_{n,m}$ were
surprisingly hard. They appear numerically ill-conditioned and include
numerical solutions of integrals making corresponding proper interval
arithmetic challenging to implement. Instead, they were computed with a
straight-forward basic Remez exchange algorithm implementation. However, to
ensure that the linear system of the algorithm iterations were non-singular
for the higher polynomial degrees and higher order Fourier coefficients, it
had to be implemented with arbitrary precision arithmetic. Consequently, the
integrals of $b_{n}(h_{c})$ and $c_{n}(h_{c})$ had to be solved numerically
with arbitrary precision. This in turn requires reference transformations
implemented in arbitrary precision. For the arbitrary precision arithmetic,
the GMP [36] (linked against), the MPFR [27] (basic arbitrary precision
arithmetic) and the MPLAPACK [66] (overloaded operators and arbitrary
precision linear algebra) libraries were used. For the reference
transformation, the closed-form transformation by [95] was used. For the
integrals, the Composite Simpson’s 1/3 rule was used with 128 sub-intervals.
(In principle something more suitable such as Clenshaw-Curtis integration
should be used but, in practice, the integrals are benign enough to use
Simpson’s.) All computations were done with 200 bits of precision. The Remez
algorithm termination criteria was set to ratio of the maximum and minimum
error being less than $1+1e^{-5}$. This gave coefficient values with $\sim$20
decimal digits of precision, giving some more precision than the $\leq$17
digits of double precision to handle final polynomial coefficients
computations of $\omega_{N,M}(u,v)$ and $\mu_{N,M}(u,v)$, which were done with
[35].

Beyond the coefficients, the polynomial evaluation has many degrees of freedom
and even if arithmetic trees of minimal size could be found, in practice, they
will typically be suboptimal due to instruction parallelism of modern
superscalar processors and, in general, a large number of permutations has to
be benchmarked for each polynomial order to find the best evaluation, a
daunting task. See [24] for a glimpse of what throughout optimization entails.
And even an optimal evaluation, put in the overall arithmetic sequence, may
not be optimal anymore. That said, Estrin’s schema has been used as a middle-
of-the-road approach. Only marginal differences has been observed in
comparison with other general evaluation schema such as direct evaluation or
Horner-2. The power factors of the schema were constructed by starting from
the smallest power and solving the minimal _change-making_ problem with
respect to the already computed powers.

With the polynomials at hand, a reasonable overall arithmetic tree and an
overall arithmetic sequence (traversal of the overall arithmetic tree/graph)
of the polynomials and the auxiliary quantities is straight-forward to find.
Pseudo-code of the approximations, defining such sequences are found in Table
1. The approximations of Table 1, as well as the literature methods, have been
implemented in C++, using IEEE-754 double precision. (Note, an ulp of IEEE-754
single precision at sea level is roughly 0.5m making it insufficient for most
applications.) Polynomial coefficients were pre-computed and other derived
constants computed at compile time. The implementations were compiled with gcc
11.4 with -O3 -march=native. -O3 enables auto-vectorization and -march=native
enables AVX instructions for the benchmark platform, which is beneficial for
the polynomial evaluations. Results without AVX instructions are also provided
in Appendix D.

$f_{\sin\\!+,\\{N_{\phi},M_{\phi},N_{h},M_{h},I,J\\}}(\\{x,y,z\\})$ | $f_{\tan\\!+,\\{N_{\phi},M_{\phi},N_{h},M_{h},I,J\\}}(\\{x,y,z\\})$
---|---
$\begin{aligned} p&=\sqrt{x^{2}+y^{2}+z^{2}}\\\ t&=z/p\\\ t^{2}&=tt\\\\[-2.84526pt] \phi&=\chi_{I}(t)+t\sqrt{1-t^{2}}\,\omega_{N_{\phi},M_{\phi}}(p,t^{2})\\\ \lambda&=\xi_{J}(y,x)\\\ h&=\mu_{N_{h},M_{h}}(p,t^{2})\\\\[2.84526pt] &\text{return}\,\\{\phi,\lambda,h\\}\end{aligned}$ |  $\begin{aligned} p&=\sqrt{x^{2}+y^{2}+z^{2}}\\\\[-2.84526pt] q&=\sqrt{x^{2}+y^{2}}\\\ t&=z/p\\\\[-1.42262pt] t^{2}&=tt\\\\[-2.84526pt] \phi&=\psi_{I}(z,q)+t\sqrt{1-t^{2}}\,\omega_{N_{\phi},M_{\phi}}(p,t^{2})\\\ \lambda&=\xi_{J}(y,x)\\\ h&=\mu_{N_{h},M_{h}}(p,t^{2})\\\\[2.84526pt] &\text{return}\,\\{\phi,\lambda,h\\}\end{aligned}$
$f_{\sin\\!\times,\\{L_{\phi},N_{\phi},M_{\phi},N_{h},M_{h},I,J\\}}(\\{x,y,z\\})$ | $f_{\tan\\!\times,\\{L_{\phi},N_{\phi},M_{\phi},N_{h},M_{h},I,J\\}}(\\{x,y,z\\})$
$\begin{aligned} p&=\sqrt{x^{2}+y^{2}+z^{2}}\\\ t&=z/p\\\ t^{2}&=tt\\\\[-1.42262pt] \omega&=\omega_{N_{\phi},M_{\phi}}(p,t^{2})\\\ \delta&=t^{2}(1-t^{2})w^{2}\\\ \sigma&=\sigma_{L_{\phi}}(\delta)\\\ \tau&=\tau_{L_{\phi}}(\omega,\delta)\\\ \phi&=\chi_{I}(t(\sigma+(1-t^{2})\tau))\\\ \lambda&=\xi_{J}(y,x)\\\ h&=\mu_{N_{h},M_{h}}(p,t^{2})\\\\[2.84526pt] &\text{return}\,\\{\phi,\lambda,h\\}\end{aligned}$ |  $\begin{aligned} p&=\sqrt{x^{2}+y^{2}+z^{2}}\\\\[-2.84526pt] q&=\sqrt{x^{2}+y^{2}}\\\ t^{2}&=(z/p)^{2}\\\ \omega&=\omega_{N_{\phi},M_{\phi}}(p,t^{2})\\\ \delta&=t^{2}(1-t^{2})w^{2}\\\ \sigma&=\sigma_{L_{\phi}}(\delta)\\\ \tau&=\tau_{L_{\phi}}(\omega,\delta)\\\ \phi&=\psi_{I}\left(z(\sigma+(1-t^{2})\tau),q(\sigma-t^{2}\tau)\right)\\\ \lambda&=\xi_{J}(y,x)\\\ h&=\mu_{N_{h},M_{h}}(p,t^{2})\\\\[2.84526pt] &\text{return}\,\\{\phi,\lambda,h\\}\end{aligned}$
$f_{\bar{n},\\{L_{n},N_{n},M_{n},N_{h},M_{h}\\}}(\\{x,y,z\\})$ |
$\begin{aligned} p&=\sqrt{x^{2}+y^{2}+z^{2}}\\\ t&=z/p\\\ t^{2}&=tt\\\\[-1.42262pt] \omega&=\omega_{N_{n},M_{n}}(p,t^{2})\\\ \delta&=t^{2}(1-t^{2})w^{2}\\\ \sigma&=\sigma_{L_{n}}(\delta)\\\ \tau&=\tau_{L_{n}}(\omega,\delta)\\\ \rho^{\prime}&=(\sigma-t^{2}\tau)/p\\\ \bar{n}&=\left[x\rho^{\prime},y\rho^{\prime},t(\sigma+(1-t^{2})\tau)\right]\\\ h&=\mu_{N_{h},M_{h}}(p,t^{2})\\\\[2.84526pt] &\text{return}\,\\{\bar{n},h\\}\end{aligned}$ |
Table 1. Pseudo-code implementation for the five different ECEF
($\\{x,y,z\\}$) to geodetic ($\\{\phi,\lambda,h\\}$ or $\\{\bar{n},h\\}$)
coordinate transformations. Equality sign imply assignment. Powers are assumed
implemented inline unless available as a left hand side variable. Note that
$h_{0}=0$ and has been eliminated. A selection of polynomials can be found in
Appendix G.

##  7\. Benchmark

The performance of the coordinate transformation methods are quantified with
pairs of computational cost and error measurements. The desirable way to
measure error may vary significantly for different applications. Hence, the
result of multiple error measures are presented. In contrast, for most
applications, the most interesting measure of computational cost, and the only
one presented, is latency.

Normally, measuring anything but mean latency for algorithms is difficult,
especially for the current algorithms since the latencies are in many cases
data dependent. Hence, the mean has to be with respect to data. To separate
the algorithm from the system influence, the minimum mean latency of multiple
tests is used. The latency results were obtained with Nanobench [55] using the
minimum of the mean latency per sample of sets of 1000 samples repeated
$10\,000$ times. This gave stable results for all methods. For further details
of the latency measurements, see Appendix F.

Error measurements were separately computed with $5\cdot 10^{8}$ uniformly
distributed ECEF samples (ignoring elliptical effects)
$\\{x^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k},y^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k},z^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k}\\}$
and geodetic references $\\{\phi_{k},\lambda_{k},h_{k}\\}$ and $\bar{n}_{k}$
over the volume of the respective altitude ranges, see [77] for details. The
surface area of the earth is $\sim 5\cdot 10^{8}\text{km}^{2}$ so it means one
sample per $\text{km}^{2}$. Extra samples were added around the poles, where
methods frequently have sharp error peaks. All error computations were done in
the Intel 80-bit extended precision format (long double on Linux), ensuring
that the limits of the IEEE-754 double precision can be observed. Let
$f_{lla/nva}(x,y,z)$ denote the transformation to be evaluated and the
corresponding _lla_ and _nva_ transformation samples

$\\{\tilde{\phi}_{k},\tilde{\lambda}_{k},\tilde{h}_{k}\\}=f_{lla}(x^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k},y^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k},z^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k})\quad\textrm{and}\quad\\{\tilde{n}_{k},\tilde{h}_{k}\\}=f_{nva}(x^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k},y^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k},z^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k})$

respectively. Further, for the _nva_ , latitude and longitude values are
computed with

$\begin{split}\tilde{\phi}_{k}&=\operatorname{atan2}({n_{k}}_{z},\sqrt{{n_{k}}_{x}^{2}+{n_{k}}_{y}^{2}})\\\
\tilde{\lambda}_{k}&=\operatorname{atan2}({n_{k}}_{y},{n_{k}}_{x})\end{split}$

with extended precision, where
$\tilde{n}_{k}=[{n_{k}}_{x};{n_{k}}_{y};{n_{k}}_{z}]$. Let
$e_{\alpha}(\phi_{k},\tilde{\phi}_{k},\lambda_{k},\tilde{\lambda}_{k},h_{k},\tilde{h}_{k})$
and
$e_{\alpha}(\phi_{k},\tilde{\phi}_{k},\lambda_{k},\tilde{\lambda}_{k},h_{k},\tilde{h}_{k},\bar{n}_{k},\tilde{n}_{k})$
denote $\alpha$ error measures for _lla_ and _nva_ transformations,
respectively. See Appendix C for the definitions of the errors measures. The
accuracy of $f_{lla/nva}(x,y,z)$ is then quantified with the maximum error

$\begin{split}&\max_{k}(|e_{\alpha}(\phi_{k},\tilde{\phi}_{k},\lambda_{k},\tilde{\lambda}_{k},h_{k},\tilde{h}_{k})|)\\\
&\max_{k}(|e_{\alpha}(\phi_{k},\tilde{\phi}_{k},\lambda_{k},\tilde{\lambda}_{k},h_{k},\tilde{h}_{k},\bar{n}_{k},\tilde{n}_{k})|)\end{split}$

respectively. (However, with the uniform samples, any statistics of the errors
could be computed.)

For the polynomial approximations, the number of possible index limit
combinations are too great to be exhaustively benchmarked. Initial narrowing
of index limits is discussed in Appendix B. This still gives a test set of
$\sim 3000$ index limit combination. Therefore, only results for selected
approximations are labeled. The _nva_ approximations have distinct error
levels corresponding to different $N$. Consequently, the selection process is
rather straight forward. In contrast, the _lla_ approximations demonstrates an
almost continuous range of errors. Consequently, a simple selection of _lla_
approximations matching the selected _nva_ approximations is used. The
selected, in the charts labeled and in Appendix G provided (corresponding
polynomials) approximations are:

$f_{\bar{n},\\{0,0,0,0,0\\}}$ |  The simplest possible _nva_ approximation corresponding to a spherical earth. The maximum Euclidean, horizontal and altitude errors are $2.1\cdot 10^{4}$m, $2.1\cdot 10^{4}$m and $1.1\cdot 10^{4}$m, respectively.
---|---
$f_{\bar{n},\\{1,1,0,1,0\\}}$ |  $N=1$ approximation with the maximum Euclidean, horizontal and altitude errors are $85$m, $85$m and $13$m. Note that there is a significant n-vector magnitude error of $2.8\cdot 10^{-6}$. This gives an equivalent gravity error which is order of magnitude the same as typical gravity anomalies.
$f_{\bar{n},\\{2,2,1,2,0\\}}$ |  $N=2$ approximation with the maximum Euclidean, horizontal and altitude errors are $0.44$m, $0.44$m and $0.10$m, respectively. The magnitude error is $9.9\cdot 10^{-12}$. Note that there are faster approximations with roughly the same Euclidean error but with a significant magnitude error. The selected approximation gives the lowest relative gravity error of $6.8\cdot 10^{-8}$ for the range, well below typical gravity anomalies.
$f_{\bar{n},\\{3,3,2,3,2\\}}$ |  $N=3$ approximation with the maximum Euclidean, horizontal and altitude errors are $1.2\cdot 10^{-3}$m, $1.2\cdot 10^{-3}$ and $4.8\cdot 10^{-4}$m, respectively. The magnitude error is $6.7\cdot 10^{-13}$.
$f_{\bar{n},\\{3,4,3,4,2\\}}$ |  $N=4$ approximation with the maximum Euclidean, horizontal and altitude errors are $4.7\cdot 10^{-6}$m, $4.7\cdot 10^{-6}$m and $1.2\cdot 10^{-6}$m, respectively. The magnitude error is again $6.7\cdot 10^{-13}$.
$f_{\bar{n},\\{4,5,4,5,3\\}}$ |  $N=5$ approximation with the maximum Euclidean, horizontal and altitude errors are $2.1\cdot 10^{-8}$m, $2.1\cdot 10^{-8}$m and $4.1\cdot 10^{-9}$m, respectively. The magnitude error is at the numerical limit at $2.9\cdot 10^{-16}$. Note, $L=3$ gives a marginally faster approximation with essentially the same Euclidean error but the magnitude error become significant giving worse $g$ approximation. (However, as noted, such $g$ errors are below typical gravity anomalies.)
$f_{\bar{n},\\{4,6,4,5,3\\}}$ |  $N=6$ approximation with the maximum Euclidean, horizontal and altitude errors are just above the numerical limits at $5.5\cdot 10^{-9}$m, $5.3\cdot 10^{-8}$m and $4.1\cdot 10^{-9}$m, respectively. The magnitude error is at the numerical limit at $2.9\cdot 10^{-16}$. Note that the altitude approximation only has $N=5$.
$f_{\sin\\!\times,\\{0,0,0,0,0,2,3\\}}$ |  The simplest possible _lla_ approximation. The Euclidean, horizontal and altitude errors are $2.0\cdot 10^{4}$m, $2.0\cdot 10^{4}$m and $1.1\cdot 10^{4}$m, respectively.
---|---
$f_{\sin\\!+,\\{1,0,1,0,4,5\\}}$ |  $N=1$ approximation with Euclidean, horizontal and altitude errors of $114$m, $113$m and $13$m, respectively.
$f_{\sin\\!+,\\{2,1,2,0,7,8\\}}$ |  $N=2$ approximation with Euclidean, horizontal and altitude errors of $0.41$m, $0.40$m and $0.10$m, respectively.
$f_{\tan\\!+,\\{3,2,3,1,12,12\\}}$ |  $N=3$ approximation with Euclidean, horizontal and altitude errors of $1.2\cdot 10^{-3}$m, $1.2\cdot 10^{-3}$m and $4.8\cdot 10^{-4}$m, respectively.
$f_{\tan\\!\times,\\{3,4,3,4,2,14,14\\}}$ |  $N=4$ approximation with Euclidean, horizontal and altitude errors of $7.9\cdot 10^{-6}$m, $7.9\cdot 10^{-6}$m and $1.2\cdot 10^{-6}$m, respectively.
$f_{\tan\\!\times,\\{4,5,4,5,3,17,17\\}}$ |  $N=5$ approximation with Euclidean, horizontal and altitude errors of $2.5\cdot 10^{-8}$, $2.5\cdot 10^{-8}$ and $4.1\cdot 10^{-9}$, respectively.
$f_{\tan\\!\times,\\{4,6,4,5,3,18,18\\}}$ |  $N=6$ approximation with Euclidean, horizontal and altitude error of $8.2\cdot 10^{-9}$, $7.7\cdot 10^{-9}$ and $4.1\cdot 10^{-9}$, respectively. Note that the altitude approximation only has $N=5$.

The results of the literature methods, against which the approximations are
benchmarked, are labeled with numbers. These numbers and the corresponding
methods are listed in chronological order below, together with short comments
about the implementation and the results. Note that _all iterative methods
have been implemented without trigonometric function evaluations in the
iteration loop!_ Methods of which trigonometric functions within the iteration
loop has not been avoidable, e.g. [49], simply have much longer latencies and
have not been included in the benchmark.

0 |  Spherical earth approximation. Not including any elementary function approximations, i.e. direct implementation of (1). This is in many ways the fastest (and least accurate) conceivable approximation and a special case of the presented method.
---|---
1-6 |  Iterative method by [44] with increasing number of iterations. The results start with zero iterations, using only the suggested initialization. Note that the method is often referred to as Hirvonen and Mortiz’s since it appear in publications bearing both names [46]. Further, it is also often referred to as Heiksanen and Mortiz’s since it appears in their 1967 book [42]. Similarly, the same methods is also referred to as Torge’s, e.g. [98], since it appears in [90]. Implemented without trigonometric function evaluations in the iteration loop.
7 |  Series expansion method by [62]. All but the necessary arctan:s have been replaced with multi-angle formulas and small angle approximations. The altitude formula has been replaced by that of Bowring.
8 |  Method by [33].
9-11 |  Method by [3]. It performs poorly around the poles and, therefore, the maximum errors levels out. Closer than 6cm from the poles, it returns NaN and this area has been replaced with a spherical approximation.
12 |  Method by [72]. It performs poorly around the equator, hence, the poor performance.
13-14 |  Original method by [9] with one and two iterations.
15 |  Closed form method by [41].
16-18 |  [10] iterative method with refined altitude calculations and increasing number of iterations.
19 |  Method by [71]. Poor conditioning makes the method perform poorly around the poles, hence the poor results. Except for the poles it achieves the numerical precision.
20-26 |  Iterative method by [101] starting from zero iterations.
27-33 |  Second iterative method by [101] starting from zero iterations.
34-37 |  Method by Goad (1987) retrieved from [98] with increasing number of iterations starting from zero.
38 |  Closed formula method by [8] excluding inside evolute handling.
39 |  Closed formula method by [104].
40-42 |  Iterative method by [58] with increasing number of iterations. The results start with zero iterations, using only the suggested initialization.
43 |  Approximation method by [68]. No n-vector method is implementable since the algorithm works directly on the latitude.
44 |  Single iteration Bowrings’s method optimized by [88]. The method is not suitable for n-vector implementation.
45-50 |  Iterative method by [80] with increasing number of iterations starting from zero. The suggested polar region iterations are not implemented but rather only the adapted altitude formula [99].
51-56 |  Iterative method by [80] with altitude computations according to the Bowring’s 1985 update [10].
57-62 |  Iterative method by [29] with increasing number of iterations.
63 |  Sofair’s revised 2000 method [82]. It does perform particularly poor at the poles but not particularly good anywhere.
64-67 |  Iterative altitude first method by [75] with increasing number of iterations. The results start with two iterations. Less iterations gives error outside the plot.
68-71 |  Iterative latitude first method by [75] with increasing number of iterations. The results start with zero iterations, using only the suggested initialization.
72-74 |  Method by [102] with 0 to 2 iterations. One square root is avoided by reformulating the inverse cotangent to an inverse tangent and for the n-vector version, the altitude and n-vector normalization is done jointly, avoiding an extra square root.
75 |  Closed formula method by [95] with n-vector version by [31].
76 |  Iterative method by [30]. Note, with more than one iteration it is possible to find points for which it diverges. It also been noted for high altitudes by Ward [100] but it may also happen for low altitudes.
77 |  Close form method by [81].
78-83 |  Iterative method _a.2_ by [25] with increasing number of iterations starting from zero iterations.
84-89 |  Iterative method _b.1_ by [25] with increasing number of iterations starting from zero iterations.
90-91 |  Method by [79].
92 |  Closed formula method by [96] excluding inside evolute handling.
93-95 |  Method I by [57] with one to three iterations.
96-97 |  Method II by [57] with one and two iterations.
98 |  Method by [50] provided in GeographicLib but with special handling for non-WGS84 spheroids removed. Note, this is an adaptation of [95].
99 |  Closed formula method by[70].
100-102 |  Perturbation method by [47] with increasing ($3^{rd}$, $4^{th}$ and $5^{th}$) order of the approximation.
103,104 |  Method by [78] and modified version by Claessens. Implementations from [17].
105,106 |  Method by [94] and modified version by Claessens. Implementations from [17].
107 |  Bowring’s method modified by [17].
108 |  Spherical by [17].
109-111 |  Method by Dave Knopp found in the Peridetic code base [52] with 1-3 iterations.

In Figure 1 to 3, the latency versus the different error measures are shown
for a nominal altitude range of $[-5000,100\,000]$m. Additional results for
smaller and larger altitude ranges as well as result without AVX instructions
are provided in Appendix D.

Figure 1. Total (Euclidean distance), horizontal (latitude error arch length),
and vertical (altitude) maximum absolute error over the volume of the altitude
interval $[-5000,100\,000]$. See Appendix C for detailed definitions. Grey
points are the complete set of tested polynomial approximations. Circled and
labeled points are approximation of particular interest. Black points are the
indicated literature methods. Clearly, the suggested approximations improves
on the achievable computational cost over the full accuracy range.

Figure 2. Latitude and longitude (arch length) and n-vector (magnitude and
direction) maximum absolute error over the volume of the altitude interval
$[-5000,100\,000]$. See Appendix C for detailed definitions. Grey points are
the complete set polynomial approximations. Black points are the indicated
literature methods. Clearly, even though the polynomial approximations
introduces a magnitude error of the n-vector, they still improve on the
achievable performance on the gravity computations.

Figure 3. Fractional gravity (Euclidean distance) maximum absolute error over
the volume of the altitude interval $[-5000,100\,000]$. See Appendix C for
detailed definitions. Grey points are the complete set polynomial
approximations. Black points are the indicated literature methods. Clearly,
even though the polynomial approximations introduces a magnitude error of the
n-vector, they still improve on the achievable performance on the gravity
computations.

##  8\. Benchmark discussion

The most important aspect of the benchmark is naturally correctness. As
discussed in Appendix F, there are many factors influencing the results which
together warrant a level of doubt in any benchmark results of this kind.
Therefore, to minimize this doubt, in Appendix E, the benchmark machinery is
cross-validated by reproducing results published elsewhere. It also contains
some discussions about other benchmarks found in the literature, specifically
note [100][98].

About the results, the first thing to note are the outer bounds. The apparent
lower numerical limits of the Euclidean, horizontal, altitude, latitude,
longitude errors are comparable to the numerical precision (ulp) of IEEE-754
double precision of $\operatorname{ulp}_{64}(a)\approx 9.3\cdot 10^{-10}$. A
further detailed numerical analysis is beyond the scope of this report. The
upper limits are comparable to the difference between the major and minor
axes, $a-b\approx 2.14\cdot 10^{4}$. The lowest latency measurements $\sim 10$
clock cycles can be expected from a square root and a handful of arithmetic
operations of the spherical transformation approximation. The upper values of
$\sim 200$ clock cycles can be expected from numerous square and cubic roots
and a large number of arithmetic operations. Naturally, all these limits only
changes marginally with varying altitude range.

The second thing to note is that the _nva_ latencies are roughly half that of
the _lla_ latencies, i.e. inverse trigonometric functions make up roughly half
the cost of an _lla_ method. The vast majority of _lla_ methods contains the
same inverse trigonometric functions, see [68] for an exception, and, hence,
the most relevant latencies are those for the _nva_ methods. The differences
between methods are simply easier to observe when not diluted by common
trigonometric components. However, the methods by Feltens [25], Gade [31] and
Knopp [52] are actually described in _nva_ form but converting most of the
methods is rather straight-forward but, as noted, not possible for all of
them, so _lla_ latencies are still relevant when comparing such methods, and
for determining the computational cost for applications requiring _lla_
representation. Note that latencies for trigonometric functions may vary
significantly from implementation to implementation, e.g. see CORE-MATH
benchmark [86], requiring care when making comparison of absolute _lla_
latency measurements from different sources.

From the _lla_ results, it may appear that the presented approximations
provide superior performance over the full accuracy range. However, the
matched trigonometric approximations may be introduced to any method and,
considering the _nva_ result, a more balanced conclusion is rather that
matched trigonometric approximations should be employed for any transformation
method. Instead, the _nva_ results, suggests that the presented approximations
provide significantly improved performance down to $\sim 10^{-5}$m accuracy.
Of course, the exact number depend on the considered altitude range and may
vary from platform to platform. For at hint of the variability of the results,
see [17].

Note, for the current approximations, the pattern that 1) the most preferable
combination have identical $N_{\phi}$ and $N_{h}$ values, 2) the Euclidean
errors are normally almost identical to the horizontal errors and 3) the
altitude errors are smaller than the horizontal errors, can be expected. From
Fig. 4, the altitude is less expensive to compute and, hence, for a given
minimum maximum Euclidean error for a computational latency, the altitude
error can be expected to be smaller. Further, for identical $N_{\phi}$ and
$N_{h}$, the latitude and altitude errors are out-of-phase, which is
preferable from an Euclidean error point of view. Combined with the fact that
the altitude error can be expected to be smaller, the maximum Euclidean error
can be expected to be dominated by the horizontal error.

Apart from the presented methods, below $10^{-5}$, the following set of
methods are within expected platform variability in terms of computational
cost: [9] with modifications of [10], [68] (only _lla_), [30] and [79]. Beyond
that, the benchmark results speaks for themselves. The horizontal, vertical,
latitude and longitude errors can be viewed as basic error components and the
Euclidean and g-errors can be viewed as examples of common application
specific combined error measures. The n-vector errors are special errors
primarily relevant for the presented method (containing explicit n-vector
approximations).

##  9\. Discussions

A set of fundamentally new ECEF to geodetic coordinate transformation
approximations based on minimax polynomials have been presented. The benchmark
results show that these approximations improves the achievable computational
cost down to an accuracy of $\sim 10^{-5}$m. Corresponding improvements are
also seen in other presented error measures and are large enough and shown in
a convincingly large benchmark to be labeled _significant_ despite natural
uncertainties associated with the latency measurements. The improvements comes
with the caveats that

  * •

The approximations are valid over a preselected altitude range
$[h_{\min},h_{\max}]$ which has to be selected _before_ the transformation are
implemented. One cannot just iterate until convergence. Note that, a zero-
range still imply an approximation range of $a-b\approx 21385$m.

  * •

For the n-vector approximations, a magnitude error is also introduced.

  * •

The performance of the methods is relatively improved by vector (AVX)
instructions, especially for high order polynomials. Non-AVX results are
provided in Appendix D.

On the other hand, for the vast majority of applications, the exemplified
ranges will likely suffice and _potential improvements_ below contains several
remedies for the problem. Further, even considering the magnitude errors, the
approximation improves the achievable performance for e.g. gravity
computations for which the n-vector magnitude enters directly. In addition to
the improved performance, the approximations are attractive in that

  * •

They resemble traditional polynomial approximations of elementary functions
and related tools and techniques may be used to improve implementations and
even to implement fixed-point versions, see e.g. [63].

  * •

They are easy to specify and implement with straight-forward unbranching code
only requiring regular arithmetic operations and the simplest elementary
function, i.e. square-root, giving minimal dependencies on standard libraries.
No special handling around the poles, the equator, etc. is required.

  * •

The latitude, longitude, altitude and n-vector approximations are naturally
separable making it possible to obtain only individual values at a lower
computational cost.

  * •

They are highly tunable in that a specific approximation may be selected to
fulfill a computational cost or accuracy requirement. Especially, the
horizontal and vertical accuracy can be set independently.

  * •

Although only exemplified with WGS84 latitude and altitude, the method is
extendable to a wider set of rotation symmetric and equatorial plan reflection
symmetric and antisymmetric properties and the underlying series expansions
may have a wider applicability within other geodesy problems enabling similar
results to be obtained for these.

  * •

The polynomial, separable and extendable nature of the components of the
transformation approximation makes them amenable to vectorization (SIMD
instructions) making them suitable for modern computational platforms.

In contrast, the greatest drawbacks of the approximations are that

  * •

A large number of coefficients need to be precomputed to arrive at an actual
approximation implementation, making such an approximation somewhat hard to
use out-of-the-box. Changing the preselected altitude range requires all
coefficients to be recomputed.

  * •

The convergence rate with respect to a greater number of coefficients
indicates that the approximations becomes inefficient for accuracies beyond
double precision, i.e. other methods have better error convergence rates.

However, coefficients covering most regular use-cases are provided in this
report and applications requiring more than double precision, i.e. subatomic
accuracy, are rare.

Clearly, the presented approximations sets a new standard for _fast_ ECEF to
geodetic coordinate transformations. Further, their polynomial nature opens up
a new realm of further improvements. Specifically, note the following overall
areas of potential improvement

  * •

_Lookup tables_ It is well known that many elementary functions are often
implemented with a combination of minimax polynomials and coefficient lookup
tables rather than just a single polynomial. Diverging from clean arithmetic
code, for many platforms, such strategies could probably be employed to
further improve the approximation performance for the geodetic power series
approximations. For the inverse trigonometric function approximations, this is
definitely the case.

  * •

_Evaluation optimization_ The expression trees of the individual (multivariat)
polynomials have numerous degrees of freedom and the naive employed Estrin’s
schema are probably suboptimal, both in a number of operations sense and in a
latency sense for a given platform. Hence, it may be improved upon both in a
general sense, see e.g. [53], and in a platform specific sense, see e.g. [24].
The fact that the $\omega_{N,M}(\cdot,\cdot)$ and $\mu_{N,M}(\cdot,\cdot)$
polynomials are Bernstein polynomials with respect to $t^{2}\in[0,1]$ could
potentially be exploited, see e.g. [15]. However, changing the individual
evaluations schema may come in conflict with the following improvement option.

  * •

_Cross polynomial evaluation optimization_ The approximations are comprised of
multiple polynomials with common factors and terms. This is only exploited
through the polynomial arguments. However, more common components may be
separated out and computed jointly in order to reduce the overall number of
required arithmetic operations. For example, for specific $N$ and $M$, joint
evaluation schema for $\mu_{N,M}(u,v)$ and $\omega_{N,M}(u,v)$ with less
arithmetic operations than the Horner scheme can be constructed. Further, as
the results with and without AVX show, vectorization is important for
performance. For platforms with vector operations, any evaluation schema has
to take this into account and manual vectorization could potentially improve
the latency. Furthermore, since the presented method can be extended to other
geodetic quantities, vectorized method computing numerous such quantities at
no extra cost may be constructed.

  * •

_Alternative polynomial approximations_ In general, regular minimax
approximations provide good accuracy-to-computational-cost performance.
However, weighted (relative) minimax, rational minimax, Padé or similar
approximations may provide better performance. See [92][93] for examples.
Other approximations may also give a different error growth relative some
nominal $h_{0}$ which may be preferable for some applications.

  * •

_Full coefficient optimization_ The polynomial approximations are split over
inner and outer polynomials, with one assuming no errors in the other.
Further, latitude, longitude and altitude errors are only connected via the
index limits selection. This is clearly suboptimal and coefficients could
potentially be adjusted to provide better overall approximations, i.e. higher-
dimensional minimax approximations. Potentially, a low hanging fruit would be
to replace the minimax optimization over $h_{c}$ with the same over a first
order altitude approximation.

  * •

_Combined sin-cos approximations_ Polynomial _lla_ transformation
approximations with accuracies below $\sim 10^{-4}$m combining $\sin^{-1}$ and
$\cos^{-1}$ approximation for different $t$ ranges with better performance
than the presented approximations can easily be constructed. Note, in this
case, neither the $\sin^{-1}$ nor the $\cos^{-1}$ approximations require a
square root argument reduction. However, for a proper implementation, errors
at the range extrema have to be matched to avoid error discontinuities.
Further, approximation orders and ranges have to be optimized. This is,
overall, a challenging but doable task.

  * •

_Series expansions for coefficient calculations_ The coefficients $b_{n}(u)$
and $c_{n}(u)$ are rather cumbersome to compute and, therefore, computing the
corresponding minimax approximations is also cumbersome. The problem is that
the integrals cannot be solved analytically and one has to resort to numerical
integration. This could be circumvented by using series expansions for the
differences. See for example [59]. The series sum could then be moved out of
the integral and each term solved analytically resulting in a new series which
would be much easier to implement and approximate with plain truncation. In
turn, this would enable use of standard tools for computing the minimax
approximations.

However, these are all very different technical areas from the work presented
here and are left for future improvements.

Based on all the results of this report, my final recommendations are: Before
selecting a ECEF to geodetic coordinate transformation method, review the
actual transformation requirements of the application. Favor _nva_
representation. In comparison with _lla_ , the computational cost is roughly
half. If _lla_ representation is required, make sure to use matched
trigonometric function approximations. When it comes to the method selection,
it is easier to list when possibly not to prefer the presented methods

  * •

If the altitude range cannot be bounded at all, instead consider methods by
[9] \+ [10] or [79] and make a look-up table in $p$ for number of iterations
until convergence. Avoid [30] since it has numerical problems at very high
altitudes. Alternatively, especially if handling of points inside the earth’s
evolute is required consider closed-form methods by [96] and [50]
(GeographicLib).

  * •

If the required accuracy is below $~{}\sim 10^{-5}m$, for _nva_
representation, also consider methods by [9] \+ [10], [30] or [79] and, for
_lla_ representation, also consider the mentioned methods and [68].

If computational resources are scarce and there is plenty of time and human
resources, consider the areas of potential improvement mentioned above. High
altitude ranges can be handled and there is likely more performance to be
found.

## References

References are sorted in alphabetical order. References which have not been
accessed are labeled accordingly together with a citation to one ore more
references citing the inaccessible source. To the extent possible, DOI:s or
URL:s are provided.

## References

  * [1] Oscar S. Adams  “Latitude developments connected with geodesy and cartography”, Special Publication 67, 1921  URL: https://geodesy.noaa.gov/library/pdfs/Special_Publication_No_67.pdf
  * [2] A. Amiri-Simkooei  “Comparison of different algorithms to transform geocentric to geodetic coordinates”  In _Survey Review_ 36, 2002, pp. 627–633  DOI: https://doi.org/10.1179/003962602791482966
  * [3] R.. Baird  “Cartesian Coordinate to Geodetic Coordinate Conversion”, White Sands Missile Range Data Reduction Directorate, 1964
  * [4] Maciej Bajorek, Marek Kulczycki and Marcin Ligas  “A comparison of iterative methods of the cubic rate convergence in the problem of transformation between Cartesian and geodetic coordinates”  In _Geomatics and Environmental Engineering_ 8, 2014, pp. 15  DOI: https://doi.org/10.7494/geom.2014.8.2.15
  * [5] Maurizio Barbarella and Marco Gatti  “On the transformation from geocentric to ellipsoidal system”  In _Bollettino di Geodesia e Scienze Affini_ 52, 1993, pp. 109–132  URL: https://www.researchgate.net/publication/322884617_On_the_transformation_from_geocentric_to_ellipsoidal_system
  * [6] W. Benning  “Der kürzeste Abstand eines in rechtwinkligen Koordinaten gegebenen Außenpunktes vom Ellipsoid”  In _Allgemeine Vermessungs-Nachrichten_ 81.11, 1974, pp. 429–433
  * [7] W.J. Berger and J.R. Ricupito  “Geodetic Latitude and Altitude of a Satelite”  In _ARS Journal_ 30.9, 1960, pp. 901–902  DOI: https://doi.org/10.2514/8.5262
  * [8] K. Borkowski  “Transformation of Geocentric to Geodetic Coordinates Without Approximations”  In _Astrophysics and Space Science_ 139, 1987, pp. 1–4  DOI: https://doi.org/10.1007/BF00643807
  * [9] B.. Bowring  “Transformation from spatial to geographical coordinates”  In _Survey Review_ 23, 1976, pp. 323–327  DOI: https://doi.org/10.1179/sre.1976.23.181.323
  * [10] B.. Bowring  In _Survey Review_ 28.218, 1985, pp. 202–206  DOI: https://doi.org/10.1179/sre.1985.28.218.202
  * [11] John P. Boyd  “Chebyshev and Fourier Spectral Methods”  Mineola, NY: Dover Publications, 2001
  * [12] Nicolas Brunie, Florent De Dinechin, Olga Kupriianova and Christoph Lauter  “Code Generators for Mathematical Functions”  In _2015 IEEE 22nd Symposium on Computer Arithmetic_ , 2015, pp. 66–73  DOI: https://doi.org/10.1109/ARITH.2015.22
  * [13] Robert C. Burtch  “A comparison of methods used in rectangular to geodetic coordinate transformations”  In _ACSM Annual Conference and Technology Exhibition_ , 2006
  * [14] Sylvain Chevillard, Mioara Maria Joldeş and Christoph Lauter  “Sollya: An Environment for the Development of Numerical Codes”  In _Mathematical Software - ICMS 2010_ , 2010, pp. 28–31  DOI: https://doi.org/10.1007/978-3-642-15582-6˙5
  * [15] Filip Chudy and Pawel Woźny  “Fast and accurate evaluation of dual Bernstein polynomials”  In _Numerical Algorithms_ 87, 2021, pp. 1001–1015  DOI: https://doi.org/10.1007/s11075-020-00996-5
  * [16] James Nohl Churchyard  “Comment on “A new algorithm for the computation of the geodetic coordinates as a function of earth-centered earth-fixed coordinates””  In _Journal of Guidance, Control and Dynamics_ 9.4, 1986, pp. 511  DOI: https://doi.org/10.2514/3.20141
  * [17] S.J. Claessens  “Efficient transformation from Cartesian to geodetic coordinates”  In _Computers & Geosciences_ 133, 2019, pp. 104307  DOI: https://doi.org/10.1016/j.cageo.2019.104307
  * [18] Eva Darulova and Anastasia Volkova  “Sound Approximation of Programs with Elementary Functions”  In _Computer Aided Verification_ Springer International Publishing, 2019, pp. 174–183  DOI: https://doi.org/10.1007/978-3-030-25543-5˙11
  * [19] Andre Deprit and Andree Deprit-Bartholome  “Conversion from geocentric to geodetic coordinates”  In _Celestial mechanics_ 12, 1975, pp. 489–493  DOI: https://doi.org/10.1007/BF01595392
  * [20] Heinrich Dörrie  “Kubische und biquadratische Gleichungen”  Berlin, Boston: Oldenbourg Wissenschaftsverlag, 1948  DOI: doi:10.1515/9783486775990
  * [21] Erhart Ecker  “Die Normalenfällung auf das Ellipsoid”  In _Österreichische Zeitschrift für Vermessungswesen_ 55.3, 1967, pp. 82–85  URL: https://www.ovg.at/de/bibliothek/vgi-die-zeitschrift/files/pdf/4053
  * [22] Bernd Eissfeller  “A Taylor Series Expansion for the Transformation Problem of Cartesian Baseline Components into Ellipsoidal Coordinate Differences”  In _Schriftenreihe, Universitärer Studiengang Vermessungswesen, Universitaet der Bundeswehr München_ , 1985, pp. 47–64  URL: https://www.unibw.de/geodaesie/bau-9-1-ingenieurgeodaesie/downloads/sonstige/heft-19.pdf
  * [23] Alberto Escapa and Toshio Fukushima  “Application of Lagrange Expansion Formula in Obtaining an Analytical Approximate Solution of Cartesian to Geodetic Transformation”, Presentation at: IUGG G05 Multi-Signal Positioning, Remote Sensing and Applications, 2019  DOI: https://doi.org/10.13140/RG.2.2.17838.25921
  * [24] Timothée Ewart, Francesco Cremonesi, Felix Schürmann and Fabien Delalondre  “Polynomial Evaluation on Superscalar Architecture, Applied to the Elementary Function $e^{x}$”  In _ACM Transactions on Mathematical Software_ 46.3, 2020  DOI: https://doi.org/10.1145/3408893
  * [25] J. Feltens  “Vector methods to compute azimuth, elevation, ellipsoidal normal, and the Cartesian ($X$, $Y$, $Z$) to geodetic ($\psi$, $\lambda$, $h$) transformation”  In _Journal of Geodesy_ 82, 2008, pp. 493–504  DOI: https://doi.org/10.1007/s00190-008-0246-5
  * [26] Hok Sum Fok and H. Iz  “A Comparative Analysis of the Performance of Iterative and Non-iterative Solutions to the Cartesian to Geodetic Coordinate Transformation”  In _Journal of Geospatial Engineering_ 5.2, 2003, pp. 61–74  URL: https://www.researchgate.net/publication/237622491_A_Comparative_Analysis_of_the_Performance_of_Iterative_and_Non-iterative_Solutions_to_the_Cartesian_to_Geodetic_Coordinate_Transformation
  * [27] Laurent Fousse et al.  “MPFR: A Multiple-Precision Binary Floating-Point Library with Correct Rounding”  New York, NY, USA: Association for Computing Machinery, 2007, pp. 13–28  DOI: https://doi.org/10.1145/1236463.1236468
  * [28] W. Fraser  “A Survey of Methods of Computing Minimax and Near-Minimax Polynomial Approximations for Functions of a Single Independent Variable”  In _Journal of ACM_ 12.3, 1965, pp. 295–314  DOI: https://doi.org/10.1145/321281.321282
  * [29] Toshio Fukushima  “Fast transform from geocentric to geodetic coordinates”  In _Journal of Geodesy_ 73, 1999, pp. 603–610  DOI: https://doi.org/10.1007/s001900050271
  * [30] Toshio Fukushima  “Transformation from Cartesian to Geodetic Coordinates Accelerated by Halley’s Method”  In _Journal of Geodesy_ 79, 2006, pp. 689–693  DOI: https://doi.org/10.1007/s00190-006-0023-2
  * [31] Gade G.H. and Gade K.  “n-vector – formulas with derivations”, 2007
  * [32] Kenneth Gade  “A Non-singular Horizontal Position Representation”  In _The Journal of Navigation_ 63.3 Cambridge University Press, 2010, pp. 395–417  DOI: https://doi.org/10.1017/S0373463309990415
  * [33] Robert H. Gersten  “Geodetic Sub-latitude and Altitude of a Space Vehicle”  In _The Journal of the Astronautical Sciences_ 3.1, 1961, pp. 28–29
  * [34] Bassford C. Getchell  “Geodetic Latitude and Altitude from Geocentric Coordinates”  In _Celestial Mechanics_ 5.3, 1972, pp. 300–302  DOI: https://doi.org/10.1007/BF01228431
  * [35] “GiNaC 1.8.7”, https://www.ginac.de, 2023
  * [36] Torbjörn Granlund and the GMP development team  “GNU MP: The GNU Multiple Precision Arithmetic Library” http://gmplib.org/, 2012
  * [37] George P. Gredan  “Transforming Cartesian coordinates X,Y,Z to Geographical coordinates $\phi$, $\lambda$, $h$”  In _The Australian Surveyor_ 44.1, 1999, pp. 55–63  DOI: https://doi.org/10.1080/00050351.1999.10558773
  * [38] J.. Guo, X.. Zhao and Y.. Wu  “Transformation from Cartesian to geodetic coordinates using Lagrange inversion theorem”  In _Acta Geodaetica Et Cartographica Sinica_ 43.10, 2014, pp. 998–1004  URL: http://xb.chinasmp.com/EN/Y2014/V43/I10/998
  * [39] Cecil Hastings, Jeanne T. Wayward and James P. Wong  “Approximations for Digital Computers”  Princeton: Princeton University Press, 1955  DOI: doi:10.1515/9781400875597
  * [40] Jr. Hedgley  “An exact transformation from geocentric to geodetic coordinates for nonzero altitudes”, 1976  URL: https://ntrs.nasa.gov/api/citations/19760012748/downloads/19760012748.pdf
  * [41] M. Heikkinen  “Geschlossene Formeln zur Berechnung räumlicher geodätischer Koordinaten aus rechtwinkligen Koordinaten”  In _Zeitschrift für Vermessungswesen_ 107.5, 1982, pp. 207–211
  * [42] W. Heiskanen and H. Moritz  “Physical Geodesy”  W.H. FreemanCompany, 1967  URL: https://archive.org/details/HeiskanenMoritz1967PhysicalGeodesy
  * [43] Friedrich Robert Helmert  “Die mathematischen und physikalischen theorieen der höheren geodäsie”  B. G. Teubner, 1880
  * [44] Reino Antero Hirvonen  “New theory of the gravimetric geodesy”, 1959
  * [45] Reino Antero Hirvonen  “New theory of the gravimetric geodesy”, 1960
  * [46] Reino Antero Hirvonen and Helmut Moritz  “Practical Computation of Gravity at High Altitudes”, 1963  URL: https://archive.org/details/DTIC_AD0420541
  * [47] Hatem Hmam  “Approximation solutions to the Cartesian to geodetic coordinate transformation problem”  In _Proceeding of the IGNSS Symposium_ , 2018
  * [48] ISO  “ISO/IEC 9899:2011 Information technology — Programming languages — C”  Geneva, Switzerland: International Organization for Standardization, 2011, pp. 683 (est.)  URL: http://www.iso.org/iso/iso_catalogue/catalogue_tc/catalogue_detail.htm?csnumber=57853
  * [49] G.. Jones  “New solutions for the geodetic coordinate transformation”, 2002, pp. 437–446  DOI: https://doi.org/10.1007/s00190-002-0267-4
  * [50] Charles F.. Karney  “Geodesics on an ellipsoid of revolution”, 2011  DOI: https://doi.org/10.48550/arXiv.1102.1215
  * [51] Charles F.. Karney  “On auxiliary latitudes”  In _Survey Review_ 56.395 Informa UK Limited, 2023, pp. 165–180  DOI: http://dx.doi.org/10.1080/00396265.2023.2217604
  * [52] Dave Knopp  “Peridetic Transformations – Math Summary”, 2021  URL: https://github.com/Stellacore/peridetic
  * [53] J. Kuipers, A. Plaat, J.A.M. Vermaseren and H.J. van den Herik  “Improving multivariate Horner schemes with Monte Carlo tree search”  In _Computer Physics Communications_ 184.11, 2013, pp. 2391–2395  DOI: https://doi.org/10.1016/j.cpc.2013.05.008
  * [54] Miljenko Lapaine  “A New Direct Solution of the Transformation Problem of Cartesian into Ellipsoidal Coordinates”  In _Determination of the Geoid_ New York, NY: Springer New York, 1991, pp. 395–404  DOI: https://doi.org/10.1007/978-1-4612-3104-2˙46
  * [55] Martin Leitner-Ankerl  “Nanobench”, https://github.com/martinus/nanobench, 2022
  * [56] Xiaoyong Li et al.  “Simplified Expansions of Common Latitudes with Geodetic Latitude and Geocentric Latitude as Variables”  In _Applied Sciences_ 12.15, 2022  DOI: https://doi.org/10.3390/app12157818
  * [57] Marcin Ligas and Piotr Banasik  “Conversion between Cartesian and geodetic coordinates on a rotational ellipsoid by solving a system of nonlinear equations”  In _Geodesy and Cartography_ 60, 2011, pp. 145–159  DOI: https://doi.org/10.2478/v10277-012-0013-x
  * [58] Kuo-Chi Lin and Jie Wang  “Transformation from geocentric to geodetic coordinates using Newton’s iteration”  In _Belletin géodésique_ 69, 1995, pp. 300–303  DOI: https://doi.org/10.1007/BF00806742
  * [59] Sheila Ann T. Long  “Derivation of transformation formulas between geodentric oand goedetic coordinates for nonzero altitudes”, 1974  URL: https://ntrs.nasa.gov/api/citations/19740021140/downloads/19740021140.pdf
  * [60] Sheila Ann T. Long  “General-altitude Transformations between Geocentric and Geodetic Coordinates”  In _Celestial Mechanics_ 12.2, 1975, pp. 225–230  DOI: https://doi.org/10.1007/BF01230214
  * [61] A.. Mikhailov  “The Reduction of Geocentric Latitude to Sea Level”  In _Soviet Astronomy_ 10.6, 1967  URL: https://articles.adsabs.harvard.edu/pdf/1967SvA....10.1052M
  * [62] John Morrison and Samuel Pines  “The reduction from geocentric coordinates”  In _The Astronomical Journal_ 66.1, 1961, pp. 15–16  DOI: https://doi.org/10.1086/108351
  * [63] Christophe Mouilleron, Amine Najahi and Guillaume Revy  “Automated Synthesis of Target-Dependent Programs for Polynomial Evaluation in Fixed-Point Arithmetic”  In _2014 16th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing_ , 2014, pp. 141–148  DOI: https://doi.org/10.1109/SYNASC.2014.27
  * [64] Jean-Michel Muller  “Elementary functions - algorithms and implementation”  Birkhäuser Boston, MA, 2006  DOI: https://doi.org/10.1007/b137928
  * [65] Jean-Michel Muller  “Elementary Functions and Approximate Computing”  In _Proceedings of the IEEE_ 108.12, 2020, pp. 2136–2149  DOI: https://doi.org/10.1109/JPROC.2020.2991885
  * [66] Maho Nakata  “MPLAPACK version 2.0.1 user manual”, 2022  DOI: https://doi.org/10.48550/arXiv.2109.13406
  * [67] John-Olof Nilsson  “A comprehensive Cartesian to geodetic coordinate transformation reference listing”, 2024  arXiv:2405.05352 [physics.geo-ph]
  * [68] Donald .K. Olson  “Converting Earth-centered, Earth-fixed coordinates to geodetic coordinates”  In _IEEE Transactions on Aerospace and Electronic Systems_ 32.1, 1996, pp. 473–476  DOI: http://doi.org/10.1109/7.481290
  * [69] Sebastián Orihuela  “Funciones de Latitud”, 2013  URL: https://www.academia.edu/7580468/Funciones_de_Latitud
  * [70] Karl Osen  “Accurate Conversion of Earth-Fixed Earth-Centered Coordinates to Geodetic Coordinates”, 2017  URL: https://hal.science/hal-01704943v2
  * [71] M.. Ozone  “Non-Iterative Solution of the $\phi$ Equation”  In _Surveying and Mapping_ 45.2, 1985, pp. 169–171
  * [72] M.. Paul  “A note on computation of Geodetic coordinates from geocentric (Cartesian) coordinates”  In _Bulletin Géodésique_ 108, 1973, pp. 135–139  DOI: https://doi.org/10.1007/BF02522075
  * [73] Kiril Pavlov  “Transformation de coordonnées rectangulaires spaciales $X$, $Y$, $Z$ en coordonnées géographiques $\phi$, $\lambda$, $H$ pour les hauteurs inférieures a 10 km”  In _Schweizerische Zeitschrift für Vermessung, Photogrammetrie und Kulturtechnik_ 66.12, 1968, pp. 381–388  DOI: https://doi.org/10.5169/seals-222321
  * [74] M. Pick  “Transformation of the spatial rectangular coordinates into the geodetic coordinates”  In _Bulletin Géodésique_ 83, 1967, pp. 21–26  DOI: https://doi.org/10.1007/BF02526103
  * [75] J. Pollard  “Iterative vector methods for computing geodetic latitude and height from rectangular coordinates”  In _Journal of Geodesy_ 76, 2002, pp. 36–40  DOI: https://doi.org/10.1007/s001900100220
  * [76] Evgeny Yakovlevich Remez  “Sur la détermination des polynômes d’approximation de degré donnée”  In _Comm. Soc. Math. Kharkov_ 10, 1934
  * [77] D. Rosca  “New uniform grids on the sphere”  In _Astronomy & Astrophysics_ 520.A63, 2010  DOI: https://doi.org/10.1051/0004-6361/201015278
  * [78] Paul Sampson  “Fitting conic sections to “very scattered” data: An iterative refinement of the Bookstein algorithm”  In _Computer Graphics and Image Processing_ 18, 1982, pp. 97–108  DOI: https://doi.org/10.1016/0146-664X(82)90101-0
  * [79] Chanfang Shu and Fei Li  “An iterative algorithm to compute geodetic coordinates”  In _Computers & Geosciences_ 36.9, 2010, pp. 1145–1149  DOI: https://doi.org/10.1016/j.cageo.2010.02.004
  * [80] Lars E. Sjöberg  “An efficient iterative solution to transform rectangular geocentric coordinates to geodetic coordinates”  In _Zeitschrift für Vermessungswesen_ 124, 1999, pp. 295–297
  * [81] Lars E. Sjöberg  “A strict transformation from Cartesian to geodetic coordinates”  In _Survey Review_ 40.308, 2008, pp. 156–163  DOI: https://doi.org/10.1179/003962608X290942
  * [82] Isaac Sofair  “Improved Method for Calculating Exact Geodetic Latitude and Altitude Revisited”  In _Journal of Guidance, Control and Dynamics_ 23.2, 2000, pp. 369  DOI: https://doi.org/10.2514/2.4534
  * [83] Duncan M’Laren Young Sommerville  “Analytical Conics”, Bell’s mathematical Series  London: G. BellSons, 1924  URL: https://archive.org/details/dli.ernet.247333
  * [84] I. Sugai  “Note on Exact Geodetic Sub-Latitude and Geodetic Altitude of a Space Vehicle”  In _The Journal of the Astronautical Sciences_ 14.3, 1967, pp. 134–137
  * [85] Hans Sünkel  “Ein nicht-iteratives Verfahren zur Transformation geodätischer Koordinaten”  In _Österreichische Zeitschrift für Vermessungswesen und Photogrammetrie_ 64.1, 1976, pp. 29–33  URL: https://www.ovg.at/de/bibliothek/vgi-die-zeitschrift/files/pdf/4216
  * [86] “The CORE-MATH project” Accessed: 2024-01-11, ”https://core-math.gitlabpages.inria.fr/”, 2022
  * [87] V. Tomelleri  “L’asteriode tetracuspide retto ellittico bi-alettato e Ie coordiante alti-geografiche ellissoidi che di un punto assegnato”  In _Bollettino di geodesia e scienze affini_ 29.1 and 2, 1970, pp. 25-74 (1) and 137–166 (2)
  * [88] Ralph Toms  “An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion”  In _Proceedings on Standards for the Interoperability of Distributed Simulations_ , 1996  URL: https://www.osti.gov/servlets/purl/231228
  * [89] Ralph Marvin Toms  “New efficient procedures for geodetic coordinate transformations”  In _Proceedings on Standards for the Interoperability of Distributed Simulations_ , 1998, pp. 1024–1033  URL: https://api.semanticscholar.org/CorpusID:115239162
  * [90] Wolfgang Torge  “Geodäsie”  De Gruyter, 1975
  * [91] Wolfgang Torge  “Geodesy”  De Gruyter, 2001  DOI: https://doi.org/10.1515/9783110879957.fm
  * [92] James D. Turner  “A non-iterative and non-singular perturbation solution for transforming Cartesian to geodetic coordinates”  In _Journal of Geodesy_ 83.2, 2009, pp. 139–145  DOI: https://doi.org/10.1007/s00190-008-0247-4
  * [93] James D. Turner and Tarek Elgohary  “A Simple Perturbation Algorithm for Inverting the Cartesian to Geodetic Transformation”  In _Mathematical Problems in Engineering_ , 2013  DOI: https://doi.org/10.1155/2013/712729
  * [94] Alexei Yu. Uteshev and Marina V. Goncharova  “Point-to-ellipse and point-to-ellipsoid distance equation analysis”  In _Journal of Computational and Applied Mathematics_ 328, 2018, pp. 232–251  DOI: https://doi.org/10.1016/j.cam.2017.07.021
  * [95] Hugues Vermeille  “Computing geodetic coordinates from geocentric coordinates”  In _Journal of Geodesy_ 78, 2004, pp. 94–95  DOI: https://doi.org/10.1007/s00190-004-0375-4
  * [96] Hugues Vermeille  “An analytical method to transform geocentric into geodetic coordinates”  In _Journal of Geodesy_ 85, 2011, pp. 105–117  DOI: https://doi.org/10.1007/s00190-010-0419-x
  * [97] T Vincenty  “Comment and Discussion, Re: Non-Iterative Solution of the $\phi$ Equation, by Mohammed Id Ozone”  In _Surveying and Mapping_ 45.3, 1985, pp. 265–267
  * [98] Robert W. Voll  “Comparing The Accuracy and Efficiency of Algorithms For Converting Cartesian to Geodetic Coordinates”, 1990  URL: http://hdl.handle.net/10919/41139
  * [99] Thomas Wahlberg  “Undersökning av strikta och iterativa metoder för omvandling från kartesiska till geodetiska koordinater”, 2009  URL: http://kth.diva-portal.org/smash/get/diva2:1064960/FULLTEXT01.pdf
  * [100] Steven Ward  “Comparison of Earth-Centered, Earth-Fixed to Geodetic Coordinate Conversion Algorithms”  In _Simulation Innovation Workshop (SIW)_ , 2020  URL: https://github.com/planet36/ecef-geodetic/blob/main/docs/2020_SIW_16_modified.pdf
  * [101] Ziqing Wei  “Positioning with NAVSTAR, The Global Positioning System”, 1986
  * [102] Yuanxin Wu, Ping Wang and Xiaoping Hu  “Algorithm of Earth-centered Earth-fixed coordinates to geodetic coordinates”  In _IEEE Transactions on Aerospace and Electronic Systems_ 39.4, 2003, pp. 1457–1461  DOI: https://doi.org/10.1109/TAES.2003.1261144
  * [103] D. Zanevičius and F. Keršys  “Technologies for calculating geodetic coordinates applying h-geometry functions”  In _Geodesy and Cartography_ 36.4, 2010, pp. 160–163  DOI: https://doi.org/10.3846/gc.2010.26
  * [104] J. Zhu  “Conversion of Earth-centered Earth-fixed coordinates to geodetic coordinates”  In _IEEE Transactions on Aerospace and Electronic Systems_ 30.3, 1994, pp. 957–961  DOI: https://doi.org/10.1109/7.303772

##  Appendix A Power series expansions

To handle the multivariate nature of the coordinate transformation, the
approach is to express the difference between (1) and $\phi$, $h$ and
$\bar{n}$ in Fourier (Chebyshev) series, with respect to $\phi_{c}$ and with
$h_{c}$ dependent coefficients. This separates the dependencies and enable
truncation of the series with close to minimax properties.To start, relations
between $\phi$ and $\phi_{c}$ and $h$ and $h_{c}$ are required. The actual
analytical relations are of less importance and it is sufficient that they are
numerically evaluable. Hence, let

$\begin{split}\phi&=f(x,y,z)\\\ h&=g(x,y,z)\end{split}$

be some reference ECEF-to-geodetic coordinate transformation of choice, e.g.
see [95]. (Note, the opaque nature of the reference means that the following
derivation may be carried out for any $x=y=0$ rotation symmetric and $z$-plane
reflection symmetric and anti-symmetric properties.) The geocentric to
geodetic coordinate transformations follows as

$\begin{split}\phi&=f_{c}(\phi_{c},h_{c})\\\
&=f((h_{c}+h_{0})\cos(\phi_{c}),0,(h_{c}+h_{0})\sin(\phi_{c}))\\\
h&=g_{c}(\phi_{c},h_{c})\\\
&=g((h_{c}+h_{0})\cos(\phi_{c}),0,(h_{c}+h_{0})\sin(\phi_{c}))\end{split}$

$f_{c}(\cdot)$ and $g_{c}(\cdot)$ are naturally defined for
$\phi_{c}\in[-\pi/2,\pi/2]$ and, hence, the differences between the geocentric
and the geodetic coordinates, $f_{c}(\phi_{c},h_{c})-\phi_{c}=\phi-\phi_{c}$
and $g_{c}(\phi_{c},h_{c})-h_{c}=h-h_{c}$, are continuous and may be viewed as
periodic with period $\pi$. Obviously, $\phi-\phi_{c}$ is antisymmetric and
$h-h_{c}$ is symmetric around $\phi_{c}=0$. Hence, they may be expressed as
the Fourier $\sin$ and $\cos$ series

$\phi-\phi_{c}=\sum^{\infty}_{n=1}b_{n}(h_{c})\sin\left(2n\phi_{c}\right)$

where

$b_{n}(h_{c})=\frac{2}{\pi}\int^{\pi/2}_{-\pi/2}(f_{c}(\phi_{c},h_{c})-\phi_{c})\sin\left(2n\phi_{c}\right)\textrm{d}\phi_{c}$

and

$h-h_{c}=\mathop{\smash{\sideset{}{{}^{\prime}}{\sum}}\vphantom{\sum}}\limits^{\infty}_{n=0}c_{n}(h_{c})\cos\left(2n\phi_{c}\right)$

where

$c_{n}(h_{c})=\frac{2}{\pi}\int^{\pi/2}_{-\pi/2}(g_{c}(\phi_{c},h_{c})-h_{c})\cos\left(2n\phi_{c}\right)\textrm{d}\phi_{c}$

Further, from the multi-angle formulas

$\begin{split}\sin(2n\phi_{c})&=\sum_{k=0}^{n-1}(-1)^{k}\binom{2n}{2k+1}\sin^{2k+1}(\phi_{c})\cos^{2n-2k-1}(\phi_{c})\\\
&=\sin(\phi_{c})\sqrt{1-\sin^{2}(\phi_{c})}\sum_{k=0}^{n-1}(-1)^{k}\binom{2n}{2k+1}\sin^{2k}(\phi_{c})(1-\sin^{2}(\phi_{c}))^{n-k}\\\
\cos(2n\phi_{c})&=\sum_{k=0}^{n}(-1)^{k}\binom{2n}{2k}\sin^{2k}(\phi_{c})\cos^{2n-2k}(\phi_{c})\\\
&=\sum_{k=0}^{n}(-1)^{k}\binom{2n}{2k}\sin^{2k}(\phi_{c})(1-\sin^{2}(\phi_{c}))^{n-k}\end{split}$

It may not be obvious why a Fourier series expansion would be suitable but the
fact that the differences are smooth and periodic gives a hint that a few
series terms could be sufficient to make a good approximation of the
differences. Further, from the multi-angle formulas, it can be seen that the
Fourier series can be expressed as a power series in $\sin(\phi_{c})=z/p$,
which can easily be computed. Combining it all gives the expressions for the
additive corrections to $\phi$ and $h$

(17)
$\begin{split}\phi&=\phi_{c}+\sin(\phi_{c})\sqrt{1-\sin^{2}(\phi_{c})}\omega(h_{c},\sin^{2}(\phi_{c}))\\\
&=\phi_{c}+t\sqrt{1-t^{2}}\omega(h_{c},t^{2})\\\
h&=\mu(h_{c},\sin^{2}(\phi_{c}))\\\ &=\mu(h_{c},t^{2})\end{split}$

where $\omega(\cdot,\cdot)$ and $\mu(\cdot,\cdot)$ are the power series

$\begin{split}\omega(u,v)&=\sum^{\infty}_{n=1}\sum_{k=0}^{n-1}(-1)^{k}\binom{2n}{2k+1}v^{k}(1-v)^{n-k-1}b_{n}(u)\\\
\mu(u,v)&=u+\mathop{\smash{\sideset{}{{}^{\prime}}{\sum}}\vphantom{\sum}}\limits^{\infty}_{n=0}\sum_{k=0}^{n}(-1)^{k}\binom{2n}{2k}v^{k}(1-v)^{n-k}c_{n}(u)\end{split}$

Further, from (17) and the Taylor series of $\sin(\phi)$ and $\cos(\phi)$
around $\phi_{c}$

$\begin{split}\frac{\sin(\phi)}{\sin(\phi_{c})}\\!\\!\\!&\,\,\,=\frac{\sin\left(\phi_{c}+\sin(\phi_{c})\sqrt{1-\sin^{2}(\phi_{c})}\omega(h_{c},\sin^{2}(\phi_{c}))\right)}{\sin(\phi_{c})}\\\
=&\sum^{\infty}_{l=0}\left.\frac{\partial^{l}}{\partial\varphi^{l}}\frac{\sin(\varphi)}{l!}\right|_{\phi_{c}}\frac{\left(\sin(\phi_{c})\sqrt{1-\sin^{2}(\phi_{c})}\omega(h_{c},\sin^{2}(\phi_{c}))\right)^{l}}{\sin(\phi_{c})}\\\
=&\sum_{l\in\mathbb{N}_{e}}(-1)^{l/2}\frac{1}{l!}\left(\sin(\phi_{c})\sqrt{1-\sin^{2}(\phi_{c})}\omega(h_{c},\sin^{2}(\phi_{c}))\right)^{l}\\\
&+\sum_{l\in\mathbb{N}_{o}}(-1)^{(l-1)/2}\frac{1}{l!}\frac{\cos(\phi_{c})}{\sin(\phi_{c})}\left(\sin(\phi_{c})\sqrt{1-\sin^{2}(\phi_{c})}\omega(h_{c},\sin^{2}(\phi_{c}))\right)^{l}\\\
=&\sum_{l\in\mathbb{N}_{e}}(-1)^{l/2}\frac{1}{l!}(\sin^{2}(\phi_{c}))^{l/2}(1-\sin^{2}(\phi_{c}))^{l/2}\omega(h_{c},\sin^{2}(\phi_{c}))^{l}\\\
&+\sum_{l\in\mathbb{N}_{o}}(-1)^{(l-1)/2}\frac{1}{l!}(\sin^{2}(\phi_{c}))^{(l-1)/2}(1-\sin^{2}(\phi_{c}))^{(l+1)/2}\omega(h_{c},\sin^{2}(\phi_{c}))^{l}\end{split}$
$\begin{split}\frac{\cos(\phi)}{\cos(\phi_{c})}\\!\\!\\!&\,\,\,=\frac{\cos\left(\phi_{c}+\sin(\phi_{c})\sqrt{1-\sin^{2}(\phi_{c})}\omega(h_{c},\sin^{2}(\phi_{c}))\right)}{\cos(\phi_{c})}\\\
=&\sum^{\infty}_{l=0}\left.\frac{\partial^{l}}{\partial\varphi^{l}}\frac{\cos(\varphi)}{l!}\right|_{\phi_{c}}\frac{\left(\sin(\phi_{c})\sqrt{1-\sin^{2}(\phi_{c})}\omega(h_{c},\sin^{2}(\phi_{c}))\right)^{l}}{\cos(\phi_{c})}\\\
=&\sum_{l\in\mathbb{N}_{e}}(-1)^{l/2}\frac{1}{l!}\left(\sin(\phi_{c})\sqrt{1-\sin^{2}(\phi_{c})}\omega(h_{c},\sin^{2}(\phi_{c}))\right)^{l}\\\
&+\sum_{l\in\mathbb{N}_{o}}(-1)^{(l+1)/2}\frac{1}{l!}\frac{\sin(\phi_{c})}{\cos(\phi_{c})}\left(\sin(\phi_{c})\sqrt{1-\sin^{2}(\phi_{c})}\omega(h_{c},\sin^{2}(\phi_{c}))\right)^{l}\\\
=&\sum_{l\in\mathbb{N}_{e}}(-1)^{l/2}\frac{1}{l!}(\sin^{2}(\phi_{c}))^{l/2}(1-\sin^{2}(\phi_{c}))^{l/2}\omega(h_{c},\sin^{2}(\phi_{c}))^{l}\\\
&+\sum_{l\in\mathbb{N}_{o}}(-1)^{(l+1)/2}\frac{1}{l!}(\sin^{2}(\phi_{c}))^{(l+1)/2}(1-\sin^{2}(\phi_{c}))^{(l-1)/2}\omega(h_{c},\sin^{2}(\phi_{c}))^{l}\end{split}$

where $\mathbb{N}_{e}$ and $\mathbb{N}_{o}$ are the sets of all positive even
and odd natural numbers. Let

$\delta(v,w)=v(1-v)w^{2}$
$\sigma(\delta)=\sum_{l\in\mathbb{N}_{e}}(-1)^{l/2}\frac{1}{l!}\delta^{l/2}\quad\text{and}\quad\tau(w,\delta)=w\sum_{l\in\mathbb{N}_{o}}(-1)^{(l-1)/2}\frac{1}{l!}\delta^{(l-1)/2}$

and

$\eta(v,w)=\sigma(\delta(v,w))+(1-v)\tau(w,\delta(v,w))$
$\rho(v,w)=\sigma(\delta(v,w))-v\tau(w,\delta(v,w))$

Then, multiplicative sin and cos corrections follow as

$\begin{split}\sin(\phi)&=\sin(\phi_{c})\eta(t^{2},\omega(h_{c},t^{2}))\\\
\cos(\phi)&=\cos(\phi_{c})\rho(t^{2},\omega(h_{c},t^{2}))\end{split}$

##  Appendix B Reasonable index limit combinations

The polynomial orders determines the absolute errors and computational costs.
Hence, they will, to some extent, be application and platform specific and
will ultimately have to be set by benchmarking different combinations, which
is done in the next section. Unfortunately, the combinations
$\\{I,J,L_{\phi},N_{\phi},M_{\phi},N_{h},M_{h}\\}\in\mathbb{N}^{7}$ or
$\\{I,J,N_{\phi},M_{\phi},N_{h},M_{h}\\}\in\mathbb{N}^{6}$, for the _lla_
approximations, and $\\{L_{n},N_{n},M_{n},N_{h},M_{h}\\}\in\mathbb{N}^{5}$,
for the _nva_ approximations, are too many to be tested exhaustively.
Therefore, the combinations have to be narrowed down based on their relative
errors which is here made based on the three assumptions

  * •

Applications typically have a horizontal accuracy requirement rather than
separate latitude and longitude accuracy requirements. Hence,
$\\{I,J,L_{\phi/n},$ $N_{\phi/n},M_{\phi/n}\\}$ are narrowed down jointly.

  * •

Applications may have vastly different requirements on horizontal and vertical
accuracy. Hence, $\\{I,J,L_{\phi/n},N_{\phi/n},M_{\phi/n}\\}$ and
$\\{N_{h},M_{h}\\}$ are narrowed down independently.

  * •

A part from position error, $L_{n}$ also affect n-vector magnitude errors,
which are hard to compare and significant for $L_{\phi}\leq 3$. Hence, for
$L_{\phi}\leq 3$ the limitations set by $L_{\phi}$ are ignored.

Starting from sufficiently high limits, to attain horizontal (latitude and
longitude) and altitude errors at the numerical limits of the IEEE-754 double
precision implementation, and lowering one limit at at a time gives the
maximum error shown in Fig. 4. Let $E^{f_{i}}_{K}(\cdot)$ denote the values
shown in Fig. 4 for the respective approximation $f_{i}$ and index limit $K$.
Start with complete sets of index limit value combinations $\\{L,N,...\\}$.
Split according to the assumptions. Let
$E^{f_{i}}_{\max}=\max(\\{E^{f_{i}}_{K}(X):X\in\\{L,N,...\\}\\}$ be the
maximum error for such a set. Then reject the combinations if either

  * •

$K>K_{\min}$ and $E^{f_{i}}_{K}(K-1)<E^{f_{i}}_{\max}$ and
$E^{f_{i}}_{K}(K)<E^{f_{i}}_{\max}0.25$

  * •

$K>K_{\min}$ and $E^{f_{i}}_{K}(K)<0.95E^{f_{i}}_{K}(K-1)$

The first conditions handles the potential constructively or destructively
inference between errors, meaning the combinations cannot be reduced to the
index limits providing limiting errors closest to each other. The second
condition handles increasing index limits not reducing the errors more. 0.25
and 0.95 are trade-offs between ensuring that all reasonable combinations are
captured and limiting the total number of tested combination. 0.25 and 0.95
give $\sim 3000$ combinations to test for the benchmark.

Figure 4. Maximum error (latitude arc error, longitude arc error and altitude
error) of the labeled approximation and index limit set to the indicated value
and the remaining index limits set to $I=20$, $J=20$, $L_{\phi/n}=6$,
$N_{\phi/h/n}=8$, $M_{\phi/h/n}=6$. This provides an approximate achievable
lower maximum error limit given individual index limits. The leveling out of
some of the curves above the numerical limit of $10^{-1}$ is caused by poor
numerical conditioning and corresponding errors of $\sin^{-1}$, i.e.
$\chi_{I}$, and $t\sqrt{1-t^{2}}$ around the poles.

##  Appendix C Error definitions

The errors used in the comparison of different approximations and literature
methods are defined as follows. Throughout, the _Earth’s prime-vertical radius
of curvature_ is
$\tilde{N}_{k}=a^{2}(a^{2}\cos^{2}(\tilde{\phi}_{k})+b^{2}\sin^{2}(\tilde{\phi}_{k}))^{-1/2}$
and _Earth’s meridional radius of curvature_
$\tilde{M}_{k}=(1-e^{2})/a^{2}\tilde{N}^{3}_{k}$. Note, for
$\\{x,y,z\\}=\\{0,0,0\\}$ only altitude error is well defined but this is not
repeated for every measure.

  * •

Euclidean error in meter. This will be

$e_{\leftrightarrow}=\sqrt{(\tilde{x}_{k}-x^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k})^{2}+(\tilde{y}_{k}-y^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k})^{2}+(\tilde{z}_{k}-z^{\mathchoice{\scalebox{1.5}[1.0]{$\displaystyle\circ$}}{\scalebox{1.5}[1.0]{$\textstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptstyle\circ$}}{\scalebox{1.5}[1.0]{$\scriptscriptstyle\circ$}}}_{k})^{2}}$

where

$\begin{split}\tilde{x}_{k}&=\left(\tilde{N}_{k}+\tilde{h}_{k}\right)\cos(\tilde{\phi}_{k})\cos(\tilde{\lambda}_{k})\\\
\tilde{y}_{k}&=\left(\tilde{N}_{k}+\tilde{h}_{k}\right)\cos(\tilde{\phi}_{k})\sin(\tilde{\lambda}_{k})\\\
\tilde{z}_{k}&=\left(\tfrac{b^{2}}{a^{2}}\tilde{N}_{k}+\tilde{h}_{k}\right)\sin(\tilde{\phi}_{k})\end{split}$

  * •

Geodetic horizontal error in meter. Computing the error is equivalent to the
_inverse geodetic problem_. However, since the errors (distances) can be
assumed small, sufficient accuracy can be achieved with a local spherical
approximation which can be further simplified

$\begin{split}e_{\frown}&=(\sqrt{1+e^{2}}/(1+e^{2}\cos^{2}(\phi_{k}))+h_{k})\operatorname{atan2}(|\tilde{n}_{k}\times\bar{n}_{k}|,\tilde{n}_{k}\cdot\bar{n}_{k})\\\
&\approx(\sqrt{1+e^{2}}/(1+e^{2}\cos^{2}(\phi_{k}))+h_{k})|\tilde{n}_{k}\times\bar{n}_{k}|/(|\tilde{n}_{k}||\bar{n}_{k}|)\end{split}$

where $\sqrt{1+e^{2}}/(1+e^{2}\cos^{2}(\phi_{k}))+h_{k}$ is the local
spherical curvature radius and $e^{2}=(1-b^{2}/a^{2})$ is the Earth’s
eccentricity.

  * •

Latitude error arch length in meter. This in general requires a numerical
solution of the _meridian arc_ but for small error
$\tilde{\phi}_{k}-\phi_{k}$, relative the change in curvature of the earth, it
can be approximated with

$e_{\phi}=(\tilde{\phi}_{k}-\phi_{k})(a+h_{k})/a\tilde{M}_{k}$

  * •

Longitude error arch length in meter

$e_{\lambda}=\begin{cases}\phi_{k}=\pi/2&0\\\
\text{otherwise}&(\tilde{\lambda}_{k}-\lambda_{k})\cos(\phi_{k})(a+h_{k})/a\tilde{N}_{k}\end{cases}$

Longitude is not defined for the polar axis and the errors there is always
$0$.

  * •

Altitude error in meter

$e_{h}=\tilde{h}_{k}-h_{k}$

  * •

Fractional gravity vector error based on Somigliana’s (WGS84) gravity formula
with second order free-air anomaly correction [91]

$\quad\quad\textsl{g}(\phi,h)=\gamma_{a}\frac{1+\kappa_{\textsl{g}}\sin(\phi)}{(1-e^{2}\sin(\phi))^{-1/2}}(1-(\kappa_{1}-\kappa_{2}\sin(\phi))h+\kappa_{3}h^{2})$

where

    * –

$\kappa_{\textsl{g}}=\frac{b\gamma_{b}}{a\gamma_{a}}-1$

    * –

$\gamma_{a}=9.780\,326\,771\,5\;\textrm{m/s}^{2}$ is the gravity at the
equator

    * –

$\gamma_{b}=9.832\,186\,368\,5\;\textrm{m/s}^{2}$ is the gravity the poles

    * –

$\kappa_{1}=2(1+f+m)/a=3.157\,04\cdot 10^{-7}\textrm{m}^{-1}$

    * –

$\kappa_{2}=4f/a=2.102\,69\cdot 10^{-9}\textrm{m}^{-1}$

    * –

$\kappa_{3}/a^{2}=7.374\,52\cdot^{-14}\textrm{m}^{-1}$

giving the error

$e_{g}=||\textsl{g}(\tilde{\phi}_{k},\tilde{h}_{k})\tilde{n}_{k}-\textsl{g}(\phi_{k},h_{k})\bar{n}_{k}||/|\textsl{g}(\phi_{k},h_{k})|$

  * •

N-vector magnitude error

$e_{|n|}=1-||\tilde{n}_{k}||$

  * •

N-vector orientation error in radians

$e_{\arg(n)}=\operatorname{atan2}(|\tilde{n}_{k}\times\bar{n}_{k}|,\tilde{n}_{k}\cdot\bar{n}_{k})$

This may be implemented as

$\begin{split}e_{\arg(n)}&=\sin^{-1}(|\tilde{n}_{k}\times\bar{n}_{k}|/|\tilde{n}_{k}|\bar{n}_{k}|)\\\
&\approx|\tilde{n}_{k}\times\bar{n}_{k}|/(|\tilde{n}_{k}||\bar{n}_{k}|)\end{split}$

for angles below $\pi/2$ and small angles, respectively.

##  Appendix D More benchmark results

Additional results for smaller and larger altitude ranges as well as for
methods compiled without AVX support are provided in Fig. 5 and Fig. 6,
respectively.

Figure 5. Euclidean error for the altitude ranges $[-5000,32\,000]$m and
$[-5000,500\,000]$m. Comparing with the Euclidean errors for the altitude
range $[-5000,100\,000]$m shown in Fig. 1, it can be observed that the errors
are stretched relative the maximum error of $\sim 10^{4}$ but that the overall
picture is the same. Note that some methods have numerical problems at very
large altitudes [100].

Figure 6. Euclidean error for all methods compiled without AVX support
(compiler flag -mno-avx). The error results are essentially identical but the
latency results are worse, especially for the polynomial methods, i.e. the
polynomial methods can to a large extent be vectorized and the lack of SIMD
operations favor other methods. However, SIMD operations has been prevalent on
CPUs over more than a decade and CPUs without SIMD might have different
arithmetic performance profiles. Therefore, too certain conclusions shouldn’t
be drawn from these results.

##  Appendix E Other benchmarks and benchmark cross-validation

The results can only partially be compared with other results found in the
literature. Foremost, there are no other benchmarks considering _nva_
performance. Further, the only extensive modern benchmark found in the
literature is [100]. There is also an early interesting benchmark [98] but it
is somewhat too old and has too poor code standard to be of general interest.
Worth mentioning is also the solid but narrow benchmark [4]. Frequently cited
is the benchmark [37] but it is too small to be of general interest. In
addition, a few more benchmarks can be found [2] and [26]. However, they are
of such poor quality that they can largely be ignored. A small benchmark can
also be found in [99] (Swedish) and a review of methods in [13]. Further, most
articles presenting transformation methods comes with some benchmark but they
are generally too small to draw far reaching conclusions.

Absolute timing results are scarce in the literature. A comparison of latency
measurements of transformation methods overlapping with those of [100] is
provided in the first plot of Fig. 7, together with a least squares line fit.
Performance variation figures presented in [17] indicate that differences of
5-10% can be expected even for identical implementations. Hence, the deviation
from the linear fit can be expected. Further, the latencies was measured on an
Intel Core i7 1185G7 (current) and on an Intel Core i9 13950HX (Ward). From
CPU single core benchmark results, the difference can be expected to be +45%
which is in line with 0.65 (+54%). Finally, the small constant part indicate
that the measurement overhead is low. Altogether, this provides an additional
confidence in the latency measurements. Some other timing results can also be
found in the literature. In [30], latencies are provided for 8 methods but
they are in the range $[500,5000]$ and clearly too old to be comparable.

Maximum error results are more common and results from Ward [100], Claessens
[17], Ligas and Banasik [57] and Lin and Wang [58] are compared against error
results produced by the current benchmark machinery in the second plot of Fig.
7. The results are a mix of Euclidean distance, latitude errors and altitude
errors, all converted to distances in meter. Ideally all results should fall
on the $y=x$ line included in the plot. Most results are spot on while some
results deviate considerably. Specifically

  1. -

All altitude error results, except those of [42], presented in [57] appear to
be off by a factor 1000. From comments in the code from [100], others appear
to have hade the same experience. The results are included in the plot
corrected for this.

  2. -

Altitude error results of [42] presented in [57] are roughly a factor of 8
off. There is no obvious explanation for this and the results have been
excluded not to bloat the plot.

  3. -

Results of [29] presented in [57] appear to have one more iteration than
stated. The results are plotted against the apparent correct number of
iterations in the plots.

  4. -

Multiple error results for Borkowski [8] are considerably higher. This is
likely because the presented benchmarks ignore or miss the poor performance of
the method very close to the poles.

  5. -

Methods by [78] adapted in [17] have a considerably lower performance than
expected. Code is provided so likely something is wrong in the original
description.

Other small variations can be explained by slight implementation differences
and the sampling not capturing the largest errors. Altogether, this provides
an additional confidence in the error measurements.

Figure 7. Latency measurement comparison of the current benchmark and latency
measurements for [100]. The latency measurements are not those of the article
but rather the latest measurements retrieved from
https://github.com/planet36/ecef-geodetic at 2023-12-11. Further, a latency of
$\operatorname{atan2}$ of 75 clock cycles has been subtracted from the latency
measurements of the current benchmark to compensate for longitude not being
computed. A least squares fit line is also shown.

##  Appendix F Some more comments about the benchmark

The intention of this report is _not_ to be a review and benchmark of methods
for ECEF to geodetic coordinate transformations. Hence, it is rather brief
when it comes to benchmark details. Nonetheless, it contains the most
extensive published benchmark of said methods which makes the benchmark
valuable in itself, and due to the scarcity of benchmarks in the literature,
some further comments about the benchmark are appropriate:

  * •

The benchmark methods have been selected from a combination of commonly cited
methods and methods with expected good performance. From [67], it can clearly
be seen not to be an exhaustive set. However, common methods are included.
Poorly described methods without a clear implementation have been excluded.
For example, the promising method [92] is lacking some coefficients. I hope I
am doing better in this report.

  * •

Micro-benchmarking on modern superscalar CPUs with modern compilers is
challenging to say the least. To mention a few challenges:

    * –

Different cores typically provide different performance. Scheduling is
correlated over at least seconds. Hence, just looping over benchmark methods
from first to last will not give stable results. Rather the benchmarks methods
will have to be looped over multiple times and results aggregated. We loop 100
times over the benchmark methods.

    * –

Many modern CPUs have dynamic frequency scaling. This makes latency
measurements sensitive in terms of time dependent on the thermal and power
supply state of the CPU and make conversion to clock cycles more difficult.
Further, some instructions like AVX may be more impacted than others.

    * –

To avoid branch prediction, the provided data have to be unpredictably drawn
from a pre-generated set. Unfortunately, to make it unpredictable, the base
set may not fit in the L1 cache which will introduce infrequent L2 fetch
latencies of 10 clock cycles. So the data set will have to be large enough but
not too large to introduce significant L2/L3 cache misses and the sample
drawing has to be randomly generated at runtime to avoid the compiler seeing
through it. We have used $2^{14}=16384$ base samples with randomly shuffled
indices to retrieve them.

    * –

Regular CPU clocks have a typical accuracies of tenth to hundreds of ns,
roughly the time of a transformation. This means sets of transformations has
to be measured. On the other hand, these measurement sets cannot be too large
since then system interruption may corrupt a large portion of the
measurements. We have used measurement sets of $1000$ samples repeated
$10\,000$ times (with different samples) for every benchmark methods loop.

    * –

Modern compilers are very good at removing unused code and it may interleave
multiple transformations. To avoid this, some benchmark loop logic need to be
made dependent on the result. We have used the micro benchmarking framework
Nanobench [55] to handle this. Coding your own benchmark harness is a
challenging task. Don’t do it!

  * •

Many iterative methods only require one or a handful of iterations to reach
numerical accuracy for reasonable altitude ranges. Iterating such methods
until convergence, rather than a fixed number of iterations, possibly
dependent on $p$, will incur a significantly higher computational cost.
Iterating until convergence and using an altitude range tuned for the
convergence of a particular method is a good way to make that method appear
better than all other methods. Don’t iterate until convergence when
benchmarking!

  * •

Clearly, different transformation methods do demonstrate different accuracies,
both in terms of absolute Euclidean error and in terms of different latitude
and altitude errors, and applications do have varying accuracy requirements.
Hence, the accuracy has to be taken into account when comparing the
computational cost of methods. Further, since the relation between
computational cost and latency is roughly logarithmic, it should be plotted in
logarithmic scale. Interpreting just commonly provided tabulated values is
hard!

  * •

Comparing _lla_ and _nva_ coordinate transformations, as noted, it can be seen
that _lla_ methods are dominated by the evaluation of inverse trigonometric
functions. This means that the core computational cost of methods are mixed by
less relevant common components. As a result, _nva_ versions should preferably
be compared.

  * •

Standard library trigonometric functions are typically implemented with some
lookup tables. This means that they are cache sensitive. Consequently,
different sample orders can give significantly different performance values.
The same goes for branch-prediction, as discussed above, which to a lesser
extent can affect values. As a result, one has to be very clear on what and in
what order samples are provided to the methods, and simply be aware of that
the results will depend on how the benchmarking is done, which is another
reason to favor _nva_ methods.

  * •

Longitude has a singularity at the poles, which is another reason for favoring
comparison of the n-vector versions. However, in practice, this is of less
concern when using double precision. First, the likelihood of randomly
sampling a problematic point is minimal (the pole would naturally be sampled
with quasi-random sequences but $\operatorname{atan2}(0,0)$ is well defined
for most implementations). Second, the problem is inherent to longitude and
can trivially be solved by simply ignoring samples too close to the poles. On
the other hand, some methods have poor latitude accuracy at the poles, see
e.g. [29], and some methods demonstrate narrow maximum error peaks at other
points, see e.g. [70]. Hence, when benchmarking with respect to maximum error,
pure random or quasi-random sampling cannot be used but some points will have
to be manually identified and added.

  * •

It may appear surprising that the accuracy for longitude vary between methods.
This is probably because the C standard (since C99) allows floating point
contraction, see Section 6.5 of [48]. GCC enables contraction by default and
different contractions will result in different numerical results which is
manifested in different maximum errors. These differences are of little
practical importance but it demonstrates to what length the compiler goes to
optimize the compiled methods. Consequently, benchmarking should be done in
some low level language such as C/C++ and one has to be clear on the compiler
and compiler options.

Few, if any, benchmarks found in the literature comes close to considering the
above points. See Appendix E for some more comments.

##  Appendix G Labelled approximations

In the following Table 5 to Table 7, the polynomials of the key
approximations, labeled in the result plots, are provided.

$\sigma_{0}(d)$ |  = $9.9999716626081590\\!\cdot\\!10^{-1}$
---|---
$\sigma_{1}(d)$ |  = $9.9999716626081590\\!\cdot\\!10^{-1}$
$\sigma_{2}(d)$ |  = $9.9999999999933083\\!\cdot\\!10^{-1}-4.9999952770986831\\!\cdot\\!10^{-1}d$
$\sigma_{3}(d)$ |  = $9.9999999999933083\\!\cdot\\!10^{-1}-4.9999952770986831\\!\cdot\\!10^{-1}d$
$\sigma_{4}(d)$ |  = $1.0000000000000000-4.9999999999989962\\!\cdot\\!10^{-1}d+4.1666643052156936\\!\cdot\\!10^{-2}d^{2}$
$\tau_{0}(w,d)$ |  = $0.0$
$\tau_{1}(w,d)$ |  = $w9.9999905541991507\\!\cdot\\!10^{-1}$
$\tau_{2}(w,d)$ |  = $w9.9999905541991507\\!\cdot\\!10^{-1}$
$\tau_{3}(w,d)$ |  = $w(9.9999999999986617\\!\cdot\\!10^{-1}-1.6666657220863013\\!\cdot\\!10^{-1}d)$
$\tau_{4}(w,d)$ |  = $w(9.9999999999986617\\!\cdot\\!10^{-1}-1.6666657220863013\\!\cdot\\!10^{-1}d)$
Table 5. Approximation polynomials of labelled ECEF to geodetic coordinate approximations as implemented with Horner-2. The coefficients are computed for $h\in[-5000,100\,000]$ and $h_{0}=0$. The inverse trigonometric functions approximations are provided in Table 7.  $\mu_{0,0}(u,v)$ |  = $-6367431.3222291581+u$
---|---
$\mu_{1,0}(u,v)$ |  = $-6378123.6318397466+u+21384.619221178389v$
$\mu_{2,0}(u,v)$ |  = $-6378136.9666263694+u+21491.297514157781v-106.67829297939183v^{2}$
$\mu_{3,1}(u,v)$ |  = $-6378136.9994595172+0.99999999993245003u+(21635.897371183597$
|
$-2.2448844415581047\\!\cdot\\!10^{-5}u)v+(-255.12110929990945+2.2892370728809709\\!\cdot\\!10^{-5}u$
|  $+(3.9094878910461248-4.4352564908397368\\!\cdot\\!10^{-7}u)v)v^{2}$
$\mu_{4,2}(u,v)$ |  = $-6378136.9999909624+0.99999999999763889u+1.5802065572606207\\!\cdot\\!10^{-19}u^{2}$
|
$+(21780.894277501458-6.7654346288251749\\!\cdot\\!10^{-5}u+3.523327586518059\\!\cdot\\!10^{-12}u^{2})v$
|
$+(-405.06255735505755+6.9626677884681548\\!\cdot\\!10^{-5}u-3.6418186823541659\\!\cdot\\!10^{-12}u^{2}$
|
$+(9.0523504657800462-2.0213005835652036\\!\cdot\\!10^{-6}u+1.2163741069533121\\!\cdot\\!10^{-13}u^{2})v)v^{2}$
|
$+(-0.19833376155550902+4.897368096169858\\!\cdot\\!10^{-8}u-3.1466289533477168\\!\cdot\\!10^{-15}u^{2})v^{4}$
$\mu_{5,3}(u,v)$ |  = $-6378136.999999837+0.99999999999993427u+(9.0086894408412461\\!\cdot\\!10^{-21}$
|
$-4.1721771358431153\\!\cdot\\!10^{-28}u)u^{2}+(21926.846232683445-0.0001359137163728073u$
|
$+(1.4164341421699455\\!\cdot\\!10^{-11}-5.529324196314687\\!\cdot\\!10^{-19}u)u^{2})v+(-557.90712954128696$
|
$+0.00014110897616511272u+(-1.4785186246301743\\!\cdot\\!10^{-11}+5.7903314015836112\\!\cdot\\!10^{-19}u)u^{2}$
|
$+(16.200215643865732-5.3621764869014012\\!\cdot\\!10^{-6}u+(6.4221743608384831\\!\cdot\\!10^{-13}$
|
$-2.7041524201192779\\!\cdot\\!10^{-20}u)u^{2})v)v^{2}+(-0.46569261099818271+1.7163926304246266\\!\cdot\\!10^{-7}u$
|
$+(-2.2001911719684914\\!\cdot\\!10^{-14}+9.6932580539601789\\!\cdot\\!10^{-22}u)u^{2}+(0.01212864332278238$
|
$-4.7225675478058751\\!\cdot\\!10^{-9}u+(6.2930011258784072\\!\cdot\\!10^{-16}-2.8522125189628692\\!\cdot\\!10^{-23}u)u^{2})v)v^{4}$
$\mu_{5,4}(u,v)$ |  = $-6378136.9999995893+0.99999999999977973u+(4.5142727491369904\\!\cdot\\!10^{-20}$
|
$-4.1725151386930229\\!\cdot\\!10^{-27}u)u^{2}+1.4635120955670929\\!\cdot\\!10^{-34}u^{4}+(22073.775734502189$
|
$-0.00022753572483565161u+(3.5589046523527045\\!\cdot\\!10^{-11}-2.7795229385686731\\!\cdot\\!10^{-18}u)u^{2}$
|
$+8.6774096768703165\\!\cdot\\!10^{-26}u^{4})v+(-713.77257463984938+0.00023830325013161956u$
|
$+(-3.7512900428756296\\!\cdot\\!10^{-11}+2.9410407690609742\\!\cdot\\!10^{-18}u)u^{2}-9.2051543802841124\\!\cdot\\!10^{-26}u^{4}$
|
$+(25.535665117187818-1.1183565979926647\\!\cdot\\!10^{-5}u+(2.0034813315984979\\!\cdot\\!10^{-12}$
|
$-1.6851285921853015\\!\cdot\\!10^{-19}u)u^{2}+5.5133913606023886\\!\cdot\\!10^{-27}u^{4})v)v^{2}+(-0.87962551332934436$
|
$+4.2975946940615239\\!\cdot\\!10^{-7}u+(-8.2360392773041159\\!\cdot\\!10^{-14}+7.2421780503870395\\!\cdot\\!10^{-21}u)u^{2}$
|
$-2.4446464398246533\\!\cdot\\!10^{-28}u^{4}+(0.026555348121184361-1.3718782035050124\\!\cdot\\!10^{-8}u$
|
$+(2.7329656904472876\\!\cdot\\!10^{-15}-2.4714925716231548\\!\cdot\\!10^{-22}u)u^{2}+8.5203151372608724\\!\cdot\\!10^{-30}u^{4})v)v^{4}$
$\omega_{0,0}(u,v)$ |  = $0.0$
$\omega_{1,0}(u,v)$ |  = $0.0066677813753770136$
$\omega_{2,1}(u,v)$ |  = $0.013446184736230014-1.0515236264437181\\!\cdot\\!10^{-9}u+(-0.00022196483792195034$
|  $+2.4231357903357331\\!\cdot\\!10^{-11}u)v$
$\omega_{3,2}(u,v)$ |  = $0.02023800884030541-3.1689947923483929\\!\cdot\\!10^{-9}u+1.6503606600562022\\!\cdot\\!10^{-16}u^{2}$
|
$+(-0.00047856006021247579+1.0339367946779417\\!\cdot\\!10^{-10}u-6.1276159848868637\\!\cdot\\!10^{-18}u^{2})v$
|
$+(1.2255186842036269\\!\cdot\\!10^{-5}-2.9789880393313049\\!\cdot\\!10^{-12}u+1.8956055600526669\\!\cdot\\!10^{-19}u^{2})v^{2}$
$\omega_{4,3}(u,v)$ |  = $0.027074537946419158-6.3663290353839142\\!\cdot\\!10^{-9}u+(6.6347135933514977\\!\cdot\\!10^{-16}$
|
$-2.5899885739079588\\!\cdot\\!10^{-23}u)u^{2}+(-0.00082614540139308665+2.6591180000141378\\!\cdot\\!10^{-10}u$
|
$+(-3.1458081288862306\\!\cdot\\!10^{-17}+1.3160531764922644\\!\cdot\\!10^{-24}u)u^{2})v+(2.7326402303242933\\!\cdot\\!10^{-5}$
|
$-9.9255846855924077\\!\cdot\\!10^{-12}u+(1.2611748377621174\\!\cdot\\!10^{-18}-5.524714330136761\\!\cdot\\!10^{-26}u)u^{2}$
|
$+(-8.3406214314548869\\!\cdot\\!10^{-7}+3.2207409051958092\\!\cdot\\!10^{-13}u+(-4.2673520748401045\\!\cdot\\!10^{-20}$
|  $+1.9262444612240382\\!\cdot\\!10^{-27}u)u^{2})v)v^{2}$
$\omega_{5,4}(u,v)$ |  = $0.033956858269743512-1.0657992633478044\\!\cdot\\!10^{-8}u+(1.6670252730549704\\!\cdot\\!10^{-15}$
|
$-1.3019553810283226\\!\cdot\\!10^{-22}u)u^{2}+4.0645825278533089\\!\cdot\\!10^{-30}u^{4}+(-0.0012704734728483523$
|
$+5.4298359021842561\\!\cdot\\!10^{-10}u+(-9.6247756886895796\\!\cdot\\!10^{-17}+8.0493882503509578\\!\cdot\\!10^{-24}u)u^{2}$
|
$-2.6240930039793388\\!\cdot\\!10^{-31}u^{4})v+(5.0407946329043319\\!\cdot\\!10^{-5}-2.4311236288225156\\!\cdot\\!10^{-11}u$
|
$+(4.6237683172662731\\!\cdot\\!10^{-18}-4.046026381146591\\!\cdot\\!10^{-25}u)u^{2}+1.3611682445895033\\!\cdot\\!10^{-32}u^{4}$
|
$+(-1.9142270212250502\\!\cdot\\!10^{-6}+9.8416784098310913\\!\cdot\\!10^{-13}u+(-1.9547420630659842\\!\cdot\\!10^{-19}$
|
$+1.7642687942146882\\!\cdot\\!10^{-26}u)u^{2}-6.074031842567329\\!\cdot\\!10^{-34}u^{4})v)v^{2}+(6.3806316097197935\\!\cdot\\!10^{-8}$
|
$-3.4050334493087929\\!\cdot\\!10^{-14}u+(6.9509879641250328\\!\cdot\\!10^{-21}-6.4057707012819846\\!\cdot\\!10^{-28}u)u^{2}$
|  $+2.2416300605419168\\!\cdot\\!10^{-35}u^{4})v^{4}$
$\omega_{6,5}(u,v)$ |  = $0.040885295839059019-1.6058512798514048\\!\cdot\\!10^{-8}u+(3.3508256373764818\\!\cdot\\!10^{-15}$
|
$-3.9268407514331456\\!\cdot\\!10^{-22}u)u^{2}+(2.452407337947089\\!\cdot\\!10^{-29}-6.3787271347333624\\!\cdot\\!10^{-37}u)u^{4}$
|
$+(-0.0018135837610434001+9.6632292183719524\\!\cdot\\!10^{-10}u+(-2.2823864407320369\\!\cdot\\!10^{-16}$
|
$+2.8625533791616239\\!\cdot\\!10^{-23}u)u^{2}+(-1.8662042533586911\\!\cdot\\!10^{-30}+5.0002116900156677\\!\cdot\\!10^{-38}u)u^{4})v$
|
$+(8.3959005641911052\\!\cdot\\!10^{-5}-5.0462939089357817\\!\cdot\\!10^{-11}u+(1.2777379582642763\\!\cdot\\!10^{-17}$
|
$-1.6756621695205854\\!\cdot\\!10^{-24}u)u^{2}+(1.1268286476356885\\!\cdot\\!10^{-31}-3.0887598999489422\\!\cdot\\!10^{-39}u)u^{4}$
|
$+(-3.6754339736361957\\!\cdot\\!10^{-6}+2.3558477370497502\\!\cdot\\!10^{-12}u+(-6.2286650409724249\\!\cdot\\!10^{-19}$
|
$+8.4234421419378366\\!\cdot\\!10^{-26}u)u^{2}+(-5.7956335293000243\\!\cdot\\!10^{-33}+1.6169761725178128\\!\cdot\\!10^{-40}u)u^{4})v)v^{2}$
|
$+(1.5042738265816568\\!\cdot\\!10^{-7}-1.0030850101437109\\!\cdot\\!10^{-13}u+(2.7300583918580043\\!\cdot\\!10^{-20}$
|
$-3.7741601218646458\\!\cdot\\!10^{-27}u)u^{2}+(2.6418278049580857\\!\cdot\\!10^{-34}-7.473104708913223\\!\cdot\\!10^{-42}u)u^{4}$
|
$+(-5.2615868937848073\\!\cdot\\!10^{-9}+3.5969760005163239\\!\cdot\\!10^{-15}u+(-9.9797774889619664\\!\cdot\\!10^{-22}$
|
$+1.4008066097662379\\!\cdot\\!10^{-28}u)u^{2}+(-9.926763165672736\\!\cdot\\!10^{-36}+2.8366680915535596\\!\cdot\\!10^{-43}u)u^{4})v)v^{4}$
Table 6. Approximation polynomials of labeled ECEF to geodetic coordinate approximations as implemented with Horner-2. The coefficients are computed for $h\in[-5000,100\,000]$ and $h_{0}=0$. The inverse trigonometric functions approximations are provided in Table 7.  $\xi_{3}(x)$ |  = $x(9.9535795470534027\\!\cdot\\!10^{1}-2.8869023791774234\\!\cdot\\!10^{1}x^{2}+7.9339041370619872\\!\cdot\\!10^{2}x^{4})$
---|---
$\xi_{5}(x)$ |  = $x(9.9986632904710004\\!\cdot\\!10^{1}-3.3030478126799267\\!\cdot\\!10^{1}x^{2}+(1.8015928178377459\\!\cdot\\!10^{1}$
|
$-8.5156335459610334\\!\cdot\\!10^{2}x^{2})x^{4}+2.0845107810858079\\!\cdot\\!10^{2}x^{8})$
$\xi_{8}(x)$ |  = $x(9.9999933557895379\\!\cdot\\!10^{1}-3.3329860785866898\\!\cdot\\!10^{1}x^{2}+(1.9946565664948636\\!\cdot\\!10^{1}$
|
$-1.3908629610769453\\!\cdot\\!10^{1}x^{2})x^{4}+(9.6421974777848821\\!\cdot\\!10^{2}-5.5912328807716874\\!\cdot\\!10^{2}x^{2}$
|
$+(2.1862959296299865\\!\cdot\\!10^{2}-4.0545676075818409\\!\cdot\\!10^{3}x^{2})x^{4})x^{8})$
$\xi_{12}(x)$ |  = $x(9.9999999943022189\\!\cdot\\!10^{1}-3.3333327040718965\\!\cdot\\!10^{1}x^{2}+(1.9999793539915971\\!\cdot\\!10^{1}$
|
$-1.4282551389892477\\!\cdot\\!10^{1}x^{2})x^{4}+(1.1083671099533618\\!\cdot\\!10^{1}-8.9411184337301530\\!\cdot\\!10^{2}x^{2}$
|
$+(7.1430746289771899\\!\cdot\\!10^{2}-5.2514556999956749\\!\cdot\\!10^{2}x^{2})x^{4})x^{8}+(3.2232566609514774\\!\cdot\\!10^{2}$
|
$-1.4721216068835910\\!\cdot\\!10^{2}x^{2}+(4.2936462844677191\\!\cdot\\!10^{3}-5.8769992093623553\\!\cdot\\!10^{4}x^{2})x^{4})x^{16})$
$\xi_{14}(x)$ |  = $x(9.9999999998328136\\!\cdot\\!10^{1}-3.3333333086677668\\!\cdot\\!10^{1}x^{2}+(1.9999989164783083\\!\cdot\\!10^{1}$
|
$-1.4285491033770566\\!\cdot\\!10^{1}x^{2})x^{4}+(1.1108488389521914\\!\cdot\\!10^{1}-9.0713495063830683\\!\cdot\\!10^{2}x^{2}$
|
$+(7.5933101852628750\\!\cdot\\!10^{2}-6.3109521418777435\\!\cdot\\!10^{2}x^{2})x^{4})x^{8}+(4.9445357766039380\\!\cdot\\!10^{2}$
|
$-3.3981220148545598\\!\cdot\\!10^{2}x^{2}+(1.8820607710988513\\!\cdot\\!10^{2}-7.6115686823698953\\!\cdot\\!10^{3}x^{2})x^{4}$
|
$+(1.9542989485087399\\!\cdot\\!10^{3}-2.3593188960412083\\!\cdot\\!10^{4}x^{2})x^{8})x^{16})$
$\xi_{17}(x)$ |  = $x(9.9999999999991586\\!\cdot\\!10^{1}-3.3333333331539649\\!\cdot\\!10^{1}x^{2}+(1.9999999885901307\\!\cdot\\!10^{1}$
|
$-1.4285710866597597\\!\cdot\\!10^{1}x^{2})x^{4}+(1.1111052308663949\\!\cdot\\!10^{1}-9.0902614649300906\\!\cdot\\!10^{2}x^{2}$
|
$+(7.6874146803011537\\!\cdot\\!10^{2}-6.6400996804419749\\!\cdot\\!10^{2}x^{2})x^{4})x^{8}+(5.7751981894433748\\!\cdot\\!10^{2}$
|
$-4.9340121701547388\\!\cdot\\!10^{2}x^{2}+(3.9762413378979828\\!\cdot\\!10^{2}-2.8629853703447561\\!\cdot\\!10^{2}x^{2})x^{4}$
|
$+(1.7310365496208254\\!\cdot\\!10^{2}-8.2096187905995513\\!\cdot\\!10^{3}x^{2}+(2.8093457170694520\\!\cdot\\!10^{3}$
|
$-6.0940031581488599\\!\cdot\\!10^{4}x^{2})x^{4})x^{8})x^{16}+6.2436108681934464\\!\cdot\\!10^{5}x^{32})$
$\xi_{18}(x)$ |  = $x(9.9999999999998558\\!\cdot\\!10^{1}-3.3333333332990472\\!\cdot\\!10^{1}x^{2}+(1.9999999975665970\\!\cdot\\!10^{1}$
|
$-1.4285713471275124\\!\cdot\\!10^{1}x^{2})x^{4}+(1.1111095442381502\\!\cdot\\!10^{1}-9.0907156394328137\\!\cdot\\!10^{2}x^{2}$
|
$+(7.6906649156526998\\!\cdot\\!10^{2}-6.6566100086412207\\!\cdot\\!10^{2}x^{2})x^{4})x^{8}+(5.8364603332581473\\!\cdot\\!10^{2}$
|
$-5.1030977240766378\\!\cdot\\!10^{2}x^{2}+(4.3267805093829844\\!\cdot\\!10^{2}-3.4100063401437874\\!\cdot\\!10^{2}x^{2})x^{4}$
|
$+(2.3697852281270090\\!\cdot\\!10^{2}-1.3702546920513891\\!\cdot\\!10^{2}x^{2}+(6.1837955998122361\\!\cdot\\!10^{3}$
|
$-2.0101395373502504\\!\cdot\\!10^{3}x^{2})x^{4})x^{8})x^{16}+(4.1436296328711983\\!\cdot\\!10^{4}-4.0407586855442829\\!\cdot\\!10^{5}x^{2})x^{32})$
$\chi_{2}(x)$ |  = $1.5702116976458603-2.0212058400388336\\!\cdot\\!10^{1}x+4.6707077880155363\\!\cdot\\!10^{2}x^{2}$
$\chi_{4}(x)$ |  = $1.5707878616372565-2.1412466400217967\\!\cdot\\!10^{1}x+(8.4666692360295800\\!\cdot\\!10^{2}$
|
$-3.5756539685967685\\!\cdot\\!10^{2}x)x^{2}+8.6486772213302613\\!\cdot\\!10^{3}x^{4}$
$\chi_{7}(x)$ |  = $1.5707963049952714-2.1459880383418343\\!\cdot\\!10^{1}x+(8.8979049893930584\\!\cdot\\!10^{2}$
|
$-5.0174715212327148\\!\cdot\\!10^{2}x)x^{2}+(3.0893053197591073\\!\cdot\\!10^{2}-1.7089810809721586\\!\cdot\\!10^{2}x$
|
$+(6.6712932596193130\\!\cdot\\!10^{3}-1.2628309167103800\\!\cdot\\!10^{3}x)x^{2})x^{4}$
$\chi_{12}(x)$ |  = $1.5707963267933006-2.1460183603240527\\!\cdot\\!10^{1}x+(8.9048588758505052\\!\cdot\\!10^{2}$
|
$-5.0792024609074237\\!\cdot\\!10^{2}x)x^{2}+(3.3671622816568150\\!\cdot\\!10^{2}-2.4303153862137010\\!\cdot\\!10^{2}x$
|
$+(1.8329796462243492\\!\cdot\\!10^{2}-1.3751038185856593\\!\cdot\\!10^{2}x)x^{2})x^{4}+(9.5272136428286426\\!\cdot\\!10^{3}$
|
$-5.5320203601637119\\!\cdot\\!10^{3}x+(2.4008869075915833\\!\cdot\\!10^{3}-6.6857373783667137\\!\cdot\\!10^{4}x)x^{2}$
|  $+8.7773781127018793\\!\cdot\\!10^{5}x^{4})x^{8}$
$\chi_{14}(x)$ |  = $1.5707963267948586-2.1460183658456538\\!\cdot\\!10^{1}x+(8.9048621136695723\\!\cdot\\!10^{2}$
|
$-5.0792770716881518\\!\cdot\\!10^{2}x)x^{2}+(3.3680571282738509\\!\cdot\\!10^{2}-2.4367246446961202\\!\cdot\\!10^{2}x$
|
$+(1.8626098028204800\\!\cdot\\!10^{2}-1.4676637467188599\\!\cdot\\!10^{2}x)x^{2})x^{4}+(1.1531109639172660\\!\cdot\\!10^{2}$
|
$-8.5675083245346239\\!\cdot\\!10^{3}x+(5.6003920887711281\\!\cdot\\!10^{3}-2.9601736597984054\\!\cdot\\!10^{3}x)x^{2}$
|
$+(1.1466969584531041\\!\cdot\\!10^{3}-2.8304614323065713\\!\cdot\\!10^{4}x+3.2965787399023498\\!\cdot\\!10^{5}x^{2})x^{4})x^{8}$
$\chi_{17}(x)$ |  = $1.5707963267948965-2.1460183660245246\\!\cdot\\!10^{1}x+(8.9048622536943916\\!\cdot\\!10^{2}$
|
$-5.0792814046672612\\!\cdot\\!10^{2}x)x^{2}+(3.3681275570828086\\!\cdot\\!10^{2}-2.4374172925348297\\!\cdot\\!10^{2}x$
|
$+(1.8670862830280291\\!\cdot\\!10^{2}-1.4876918953706647\\!\cdot\\!10^{2}x)x^{2})x^{4}+(1.2172747383436669\\!\cdot\\!10^{2}$
|
$-1.0071291531222582\\!\cdot\\!10^{2}x+(8.2100842886296591\\!\cdot\\!10^{3}-6.3278993837482987\\!\cdot\\!10^{3}x)x^{2}$
|
$+(4.3671279910859448\\!\cdot\\!10^{3}-2.5369588366614161\\!\cdot\\!10^{3}x+(1.1590307122361844\\!\cdot\\!10^{3}$
|
$-3.8309521446820528\\!\cdot\\!10^{4}x)x^{2})x^{4})x^{8}+(8.0476117495016973\\!\cdot\\!10^{5}-8.0043584568231787\\!\cdot\\!10^{6}x)x^{16}$
$\chi_{18}(x)$ |  = $1.5707963267948966-2.1460183660253421\\!\cdot\\!10^{1}x+(8.9048622545903351\\!\cdot\\!10^{2}$
|
$-5.0792814434974340\\!\cdot\\!10^{2}x)x^{2}+(3.3681284430401487\\!\cdot\\!10^{2}-2.4374295730838710\\!\cdot\\!10^{2}x$
|
$+(1.8671988159799984\\!\cdot\\!10^{2}-1.4884115235275848\\!\cdot\\!10^{2}x)x^{2})x^{4}+(1.2206039463975348\\!\cdot\\!10^{2}$
|
$-1.0185410804411290\\!\cdot\\!10^{2}x+(8.5042859077692046\\!\cdot\\!10^{3}-6.9026076875376219\\!\cdot\\!10^{3}x)x^{2}$
|
$+(5.2181956820067259\\!\cdot\\!10^{3}-3.4855347796590377\\!\cdot\\!10^{3}x+(1.9414549204238104\\!\cdot\\!10^{3}$
|
$-8.4603434241738747\\!\cdot\\!10^{4}x)x^{2})x^{4})x^{8}+(2.6620763505161977\\!\cdot\\!10^{4}-5.3242183459453063\\!\cdot\\!10^{5}x$
|  $+5.0486339748452772\\!\cdot\\!10^{6}x^{2})x^{16}$
Table 7. Inverse trigonometric function approximation polynomials of labeled
ECEF to geodetic coordinate (lla) approximations as implemented with Horner-2.