2. Definitions of the Legendre functions

To examine the EGM2008 normalisation unambiguously, we add ‘EGM’ superscripts to the generic Equation (1.1) to indicate the NGA-normalised $P_{nm}$ and their fitting coefficients $C_{nm}$, $S_{nm}$:

(2.1)\begin{equation} \Phi = \frac{-GM}{r}\left[1 + \sum_{n\, =\, 2}^{\infty} \left(\frac{a}{r}\right)^{\!n} \sum_{m\, =\, 0}^{n} P^{\rm EGM}_{nm}(\cos\theta) \big(C^{\rm EGM}_{nm}\cos m\phi + S^{\rm EGM}_{nm}\sin m\phi\big)\right]. \end{equation}

Our task is to compute the gravitational field ${\boldsymbol {g}}$ from Equation (2.1) by calculating the associated Legendre functions $P^{\rm EGM}_{nm}$, and using the coefficients $C^{\rm EGM}_{nm}, S^{\rm EGM}_{nm}$ provided in the EGM2008 database.

We begin with a clean slate, by defining the following Legendre polynomials and associated functions. We avoid the more standard ‘$P$’ notation so that nothing will be confused with normalisations used by others.

(2.2)\begin{align} \text{Legendre polynomial:} \quad & {K}_n(x) \equiv \frac{1}{2^{n} n\textit{!}}\frac{\mathrm{d}^{n}}{\mathrm{d} x^{n}}(x^{2}-1)^{n} \quad (n = 0,1,2,\ldots), \end{align}

(2.3)\begin{align} \text{Associated Legendre function:} \quad & {K}_{nm}(x) \equiv (1-x^{2})^{m/2}\frac{\mathrm{d}^{m}}{\mathrm{d} x^{m}}{K}_n(x) \quad (m = 0,1,2,\ldots, n). \end{align}

(Some authors extend their Legendre functions to negative values of $m$, but we have no need to do so.) It follows trivially that

(2.4)\begin{equation} {K}_{n0}(x) = {K}_n(x). \end{equation}

In their standard geodesy reference, Hofmann–Wellenhof and Moritz (Reference Hofmann–Wellenhof and Moritz2006) define

(2.5)\begin{equation} \text{Hofmann--Wellenhof & Moritz:} \quad P_n(x) \equiv {K}_n(x),\quad P_{nm}(x) \equiv {K}_{nm}(x). \end{equation}

Gradshteyn and Ryzhik's table of functions and integrals (Gradshteyn and Ryzhik, Reference Gradshteyn and Ryzhik1994), and the software packages Matlab and Mathematica, define

(2.6)\begin{equation} \text{Grad.--Ryzh., Matlab, Mathematica:} \quad P_n(x) \equiv {K}_n(x),\quad P_n^{m}(x) \equiv ({-}1)^{m} {K}_{nm}(x). \end{equation}

Our aim is to express $P^{\rm EGM}_{nm}$ in terms of ${K}_{nm}$. This requires a knowledge of some identities of the various functions, as discussed next.

2.1 Useful Legendre identities

Legendre polynomials satisfy the following recurrence relation (Gradshteyn and Ryzhik, Reference Gradshteyn and Ryzhik1994):

(2.7)\begin{equation} (n+1){K}_{n+1}(x) = (2n+1)\,x\, {K}_n(x) - n {K}_{n-1}(x). \end{equation}

We will require the derivative of an associated Legendre function:

(2.8)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} x}{K}_{nm}(x) = \frac{(n-m+1){K}_{n+1,m}(x)-(n+1) x {K}_{nm}(x)}{x^{2}-1}. \end{equation}

This derivative can be verified quickly using Mathematica. Do this by first referring to Equation (2.6) to write

(2.9)\begin{equation} {K}_{nm}(x) = ({-}1)^{m} P_n^{m}(x), \end{equation}

where $P_n^{m}(x)$ is Mathematica's notation, returned by Mathematica's function LegendreP[n,m,x]. Then, mixing maths (in italic) with Mathematica code (in typewriter font), we obtainFootnote ¹

(2.10)\begin{align} \frac{\mathrm{d}}{\mathrm{d} x}\,{K}_{nm}(x) & = \frac{\mathrm{d}}{\mathrm{d} x}\,(-1)^{m} P_n^{m}(x) = \texttt{D[(-1)$^{{\texttt{m}}}$ LegendreP[n,m,x], x]}\nonumber\\[1ex] & = \frac{\texttt{(-1)$^{{\texttt{m}}}$ (-1-n) x LegendreP[n,m,x] + (-1)$^{{\texttt{m}}}$ (1-m+n) LegendreP[1+n,m,x]}}{\texttt{x$^{{\texttt{2}}}$-1}}\nonumber\\[1ex] & = \frac{-(n+1)\,x\, ({-}1)^{m} P_n^{m}(x) + (n-m+1) ({-}1)^{m} P_{n+1}^{m}(x)}{x^{2}-1}\nonumber\\[1ex] & = \frac{(n-m+1){K}_{n+1,m}(x)-(n+1)\,x\,{K}_{nm}(x)}{x^{2}-1}. \end{align}

This last expression is the right-hand side of Equation (2.8), as required to be shown.

2.2 Normalisation of associated Legendre functions

Our associated Legendre functions are orthogonal, satisfying the following identity (Gradshteyn and Ryzhik, Reference Gradshteyn and Ryzhik1994):

(2.11)\begin{equation} \int_{{-}1}^{1} {K}_{nm}(x){K}_{n'm}(x)\,\mathrm{d} x = \frac{1}{n+1/2}\frac{(n+m)\textit{!}}{(n-m)\textit{!}}\;{\rm \Delta}_{nn'}. \end{equation}

This suggests that we define a normalised associated Legendre function as follows, such as used in Davies’ book on quantum mechanics (Davies, Reference Davies1984):

(2.12)\begin{equation} \widehat{{K}}_{nm}(x) \equiv \sqrt{{(n+1/2)(n-m)\textit{!}\over (n+m)\textit{!}}}\;{K}_{nm}(x). \end{equation}

The $\widehat {{K}}_{nm}(x)$ are thus orthonormal:

(2.13)\begin{equation} \int_{{-}1}^{1} \widehat{{K}}_{nm}(x)\widehat{{K}}_{n'm}(x)\,\mathrm{d} x = {\rm \Delta}_{nn'}. \end{equation}

By a change of variables, it is straightforward to infer from Equation (2.13) that

(2.14)\begin{equation} \int_0^{\pi} \widehat{{K}}_{nm}(\cos\theta)\widehat{{K}}_{n'm}(\cos\theta)\sin\theta\,\mathrm{d}\theta = {\rm \Delta}_{nn'}. \end{equation}

This last integral is required when integrating over a solid angle $\mathrm {d}\Omega \equiv \sin \theta \,\mathrm {d}\theta \,\mathrm {d}\phi$, which we do next. Note that the values of $\widehat {{K}}_{nm}(x)$ for $m=0,\ldots,n$ are returned as an array by the Matlab command legendre(n, x, ’norm’).

2.3 Construction of spherical harmonics

Working with potentials in spherical polar coordinates is routinely accomplished by a series expansion in spherical harmonic functions. These functions are typically the associated Legendre functions multiplied by complex exponentials, namely ${K}_{nm}(\cos \theta )\, {e}^{{i} m\phi }$, with some normalisation included. The EGM2008 formula for potential, Equation (2.1), separates these spherical harmonics into their real and imaginary parts, ${K}_{nm}(\cos \theta )\cos m\phi$ and ${K}_{nm}(\cos \theta )\sin m\phi$. We must examine these parts to investigate the normalisation of the ${K}_{nm}$ used by the NGA in its tables of coefficients labelled $C^{\rm EGM}_{nm}, S^{\rm EGM}_{nm}$ in Equation (2.1).

To proceed, define two families of these real-valued spherical harmonics as follows, using real normalisation constants $\mathcal {N}_c, \mathcal {N}_s$:

(2.15)\begin{align} {F_{nm}}(\theta,\phi) & \equiv \mathcal{N}_c \widehat{{K}}_{nm}(\cos\theta)\cos m\phi \quad (m \text{ can equal } 0), \end{align}

(2.16)\begin{align} {G_{nm}}(\theta,\phi) & \equiv \mathcal{N}_s \widehat{{K}}_{nm}(\cos\theta)\sin m\phi \quad (m \not= 0). \end{align}

These normalisation constants involve an integration over a solid angle of $4\pi$. We use the following shorthand for an integral of any function $f(\theta,\phi )$ over that solid angle:

(2.17)\begin{equation} \iint f(\theta,\phi)\,\mathrm{d}\Omega \equiv \int_0^{2\pi}\int_0^{\pi} f(\theta,\phi)\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi. \end{equation}

We now calculate these constants $\mathcal {N}_c, \mathcal {N}_s$ by demanding that

(2.18)\begin{equation} \iint{F_{nm}}(\theta,\phi) {F_{n'm}}(\theta,\phi)\,\mathrm{d}\Omega = \iint{G_{nm}}(\theta,\phi) {G_{n'm}}(\theta,\phi)\,\mathrm{d}\Omega = {\rm \Delta}_{nn'}. \end{equation}

2.3.1 Calculation of $\mathcal {N}_c$

For the case of ${F_{nm}}$ in Equation (2.15), Equation (2.18) says that

(2.19)\begin{equation} \int_0^{2\pi} \int_0^{\pi} \mathcal{N}_c^{\, 2} \widehat{{K}}_{nm}(\cos\theta) \widehat{{K}}_{n'm}(\cos\theta)\cos^{2} m\phi \sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi = {\rm \Delta}_{nn'}. \end{equation}

This separates into a product of integrals over $\theta$ and $\phi$:

(2.20)

It follows that

(2.21)\begin{equation} \mathcal{N}_c = \begin{cases} 1/\sqrt{2\pi} & :m = 0,\\ 1/\sqrt{\pi} & :m \not= 0. \end{cases} \end{equation}

2.3.2 Calculation of $\mathcal {N}_s$

The calculation of $\mathcal {N}_s$, in Equation (2.16), proceeds in the same way as in Equations (2.19)–(2.21), but now only the $m\not = 0$ case exists. The analogous expression for Equation (2.20) replaces $\mathcal {N}_c$ with $\mathcal {N}_s$ and $\cos m\phi$ with $\sin m\phi$:

(2.22)

Hence,

(2.23)\begin{equation} \mathcal{N}_s = 1/\sqrt{\pi} = [\mathcal{N}_c \text{ for } m \not= 0]. \end{equation}

Finally, Equations (2.15) and (2.16) become

(2.24)\begin{align} {F_{nm}}(\theta,\phi) &= \genfrac{\{}{\}}{0pt}{}{1/\sqrt{2\pi} \quad (m=0)}{1/\sqrt{\pi} \quad\enspace (m\not= 0)}\times \widehat{{K}}_{nm}(\cos\theta)\cos m\phi, \end{align}

(2.25)\begin{align}G_{nm}(\theta ,\phi ) &= 1/\sqrt \pi \;\widehat{K}_{nm}(\cos \theta )\sin m\phi \quad (m\not{ = }0).\end{align}

Now introduce a factor of $\sqrt {4\pi }$ (we will see why shortly):

(2.26)

(2.27)

Treat ${G_{nm}}$ in the same way:

(2.28)

Hofmann–Wellenhof and Moritz (Reference Hofmann–Wellenhof and Moritz2006) defined the following two functions as their spherical harmonics:

(2.29)\begin{equation} \bar{\mathcal{R}}_{nm}(\theta,\phi)\equiv \sqrt{4\pi}{F_{nm}}(\theta,\phi),\quad \bar{\mathcal{S}}_{nm}(\theta,\phi)\equiv \sqrt{4\pi}{G_{nm}}(\theta,\phi). \end{equation}

They did this to enforce the following normalisation of their spherical harmonics:

(2.30)\begin{equation} \iint \bar{\mathcal{R}}_{nm}(\theta,\phi)\bar{\mathcal{R}}_{n'm}(\theta,\phi) \,\mathrm{d}\Omega = \iint \bar{\mathcal{S}}_{nm}(\theta,\phi)\bar{\mathcal{S}}_{n'm}(\theta,\phi) \,\mathrm{d}\Omega = 4\pi{\rm \Delta}_{n'n}, \end{equation}

since doing so ensures that the average values of $\bar {\mathcal {R}}^{2}_{nm}$ and $\bar {\mathcal {S}}^{2}_{nm}$ over the whole sphere both equal one (not that this carries any significance). Next, they defined a ‘fully normalised’ associated Legendre function $\bar {P}_{nm}(x)$ such that

(2.31)\begin{equation} \bar{\mathcal{R}}_{nm}(\theta,\phi) = \bar{P}_{nm}(\cos\theta)\cos m\phi,\quad \bar{\mathcal{S}}_{nm}(\theta,\phi) = \bar{P}_{nm}(\cos\theta)\sin m\phi. \end{equation}

It follows that their ‘fully normalised associated Legendre function’ is

(2.32)\begin{align} \bar{P}_{nm}(x) &= \sqrt{\frac{(n+1/2)(n-m)\textit{!}}{(n+m)\textit{!}}}{K}_{nm}(x) \times \genfrac{\{}{\}}{0pt}{}{\sqrt{2}\quad (m=0)}{\;\;2 \quad\; (m\not= 0)}\nonumber\\ &= \sqrt{\frac{k_m(2n+1)(n-m)\textit{!}}{(n+m)\textit{!}}}{K}_{nm}(x),\quad \text{where } k_m \equiv \begin{cases} 1 & :m = 0,\\ 2 & :m \not= 0. \end{cases} \end{align}

The square-root coefficient above appears in Hofmann–Wellenhof and Moritz (Reference Hofmann–Wellenhof and Moritz2006), which is a standard reference in the world of geodesy, and so its ‘fully normalised’ $\bar {P}_{nm}$ in Equation (2.32) are now a staple of that subject. EGM2008 uses these $\bar {P}_{nm}$ (see NGA-WGS84, 2014, p. 5-3), and so we must do likewise to apply that model. For transparency, I will denote Hofmann–Wellenhof & Moritz's $\bar {P}_{nm}$ of Equation (2.32) by $P^{\rm EGM}_{nm}$, as used in Equation (2.1). Geodesy's ‘fully normalised’ $P^{\rm EGM}_{nm}(x)$ can be related to the more standard ‘fully normalised’ associated Legendre function $\widehat {{K}}_{nm}(x)$ used outside geodesy by comparing Equation (2.32) with Equation (2.12), to obtain

(2.33)\begin{equation} P^{\rm EGM}_{nm}(x) = \widehat{{K}}_{nm}(x) \times\begin{cases} \sqrt{2} & :m = 0,\\ \;\;2 & :m \not= 0. \end{cases} \end{equation}

Note that Equation (2.32) is relevant to spherical harmonics, because the square-root coefficient contains a contribution from integrating $\cos ^{2} m\phi$ and $\sin ^{2} m\phi$. But these last two functions of $\phi$ have nothing to do with the associated Legendre function ${K}_{nm}(\cos \theta )$, which is a function of $\theta$ only. It follows that geodesy's ‘fully normalised’ associated Legendre function $\bar {P}_{nm}$ in Equation (2.32) is not, in fact, a true normalisation of an associated Legendre function ${K}_{nm}$. Indeed, it is clear from Equations (2.13) and (2.33) that

(2.34)\begin{equation} \int_{{-}1}^{1} P^{\rm EGM}_{nm}(x) P^{\rm EGM}_{n'm}(x)\,\mathrm{d} x = {\rm \Delta}_{nn'}\times\begin{cases} 2 & :m = 0,\\ 4 & :m \not= 0. \end{cases} \end{equation}

This is not a normalisation in the usual sense of the word; but because it was termed ‘full normalisation’ by Hofmann–Wellenhof & Moritz, most geodesists have followed suit, including the NGA. Of course, Equation (2.1) does use spherical harmonics. Nonetheless, associated Legendres are functions in their own right, and they already have their own meaningful normalisation – Equations (2.12) and (2.13) – which does not, and should not, contain any contribution from integrating the $\cos ^{2} m\phi$ and $\sin ^{2} m\phi$ functions that have nothing to do with associated Legendres.

Conventions such as normalisations should be meaningful and useful. If you require a normalised associated Legendre function outside the EGM2008 realm, use the mathematically meaningful normalisation in Equation (2.12) in place of Equation (2.32).

2.4 An example: the $J_2$ model of gravity

We can use the above expressions to generate a close approximation to the ‘$J_2$ model’ of Earth's gravity from Equation (2.1). The $J_2$ model suffices for most applications today, and is a fit of the ‘$n=2$, $m=0$’ truncation of Equation (2.1) to Earth's potential. That truncation itself is

(2.35)\begin{equation} \frac{\Phi}{-GM/r} \simeq 1 + \left(\frac{a}{r}\right)^{\!2} P^{\rm EGM}_{20}(\cos\theta) C^{\rm EGM}_{20}. \end{equation}

Use Equation (2.32) (with $P^{\rm EGM}_{nm}$ identical to $\bar {P}_{nm}$, as noted above) to write

(2.36)

Equation (2.35) then becomes

(2.37)\begin{equation} \frac{\Phi}{-GM/r} \simeq 1 + \left(\frac{a}{r}\right)^{\!2} C^{\rm EGM}_{20}\sqrt{5}(3\cos^{2}\theta-1)/2. \end{equation}

The $J_2$ model writes a constant ‘$-J_2$’ in place of the $C^{\rm EGM}_{20}\sqrt {5}$ in Equation (2.37):

(2.38)\begin{equation} J_2 \text{ model}: \frac{\Phi}{-GM/r} \simeq 1 - J_2 \left(\frac{a}{r}\right)^{\!2}(3\cos^{2}\theta-1)/2. \end{equation}

It evaluates $J_2$ by best-fitting Equation (2.38) to Earth's measured potential. In fact, the $(n,m)= (2,0)$ spherical harmonic gives the dominant contribution to $\Phi$, and so it is valid to approximate $-J_2$ by $C^{\rm EGM}_{20}\sqrt {5}$. The value of $C^{\rm EGM}_{20}$ is listed as approximately $-4{\cdot }842 \times {10^{-4}}$ in the file available on the NGA's web page (NGA, 2021). It follows that

(2.39)\begin{equation} J_2 \simeq{-}C^{\rm EGM}_{20}\sqrt{5} \simeq 4{\cdot}842\times {10^{{-}4}} \sqrt{5} \simeq 1{\cdot}083\times {10^{{-}3}}. \end{equation}

This value matches the one found in Table 3.5 in NGA-WGS84 (2014). Equation (2.38) matches equation (38) in Siouris (Reference Siouris2009), provided we:

• • realise that the $V$ in Siouris's equation (38) equals $-\Phi$;

• • truncate Siouris's equation (38) after its second term;

• • realise that $\phi _c$ in Siouris's equation (38) is geocentric latitude, equal to our $\pi /2-\theta$;

• • use the value for $C_{2,0}$ in Siouris's equation (38) that is listed after Siouris's equation (40).

3. Calculation of the gravity vector

The potential $\Phi$ is given by the EGM2008 formula, Equation (2.1), but calculating the gravity vector ${\boldsymbol {g}} = -\nabla \Phi$ requires a little more work. Begin with the following vector basis in spherical polar coordinates $r,\theta,\phi$:

• ${\boldsymbol {u}}_r =$ unit vector in the direction of increasing $r$ (that is, geocentric ‘up’);

• ${\boldsymbol {u}}_\theta =$ unit vector in the direction of increasing $\theta$ (geocentric ‘south’);

• ${\boldsymbol {u}}_\phi =$ unit vector in the direction of increasing $\phi$ (‘east’).

Note that these vectors use a geocentric view, not geodetic. In terms of these, the gravity vector is

(3.1)\begin{align} {\boldsymbol{g}} ={-}\nabla \Phi & = {-\atop}\!\frac{\partial\Phi}{\partial r}{\boldsymbol{u}}_r - \frac{1}{r}\frac{\partial\Phi}{\partial\theta}{\boldsymbol{u}}_\theta -\frac{1} {r\sin\theta}\frac{\partial\Phi}{\partial\phi}{\boldsymbol{u}}_\phi \end{align}

(3.2)\begin{align} & \equiv {g_r} {\boldsymbol{u}}_r + {g_\theta} {\boldsymbol{u}}_\theta + {g_\phi} {\boldsymbol{u}}_\phi. \end{align}

We require the partial derivatives of $\Phi$ with respect to $r,\theta,\phi$. Calculating these from Equation (2.1) with respect to $r$ and $\phi$ presents no problems. They are

(3.3)\begin{align} \frac{\partial\Phi}{\partial r} & = \frac{GM}{r^{2}}\left[1 + \sum_n(1+n)\left(\frac{a}{r}\right)^{n}\sum_m P^{\rm EGM}_{nm}(\cos\theta) (C^{\rm EGM}_{nm}\cos m\phi + S^{\rm EGM}_{nm}\sin m\phi)\right], \end{align}

(3.4)\begin{align} \frac{\partial\Phi}{\partial\phi} & = \frac{-GM}{r}\sum_n\left(\frac{a}{r}\right)^{n}\sum_m mP^{\rm EGM}_{nm}(\cos\theta) (\mathord{-}C^{\rm EGM}_{nm}\sin m\phi + S^{\rm EGM}_{nm}\cos m\phi). \end{align}

The derivative $\partial \Phi /\partial \theta$ must be tackled with care, because it requires $\mathrm {d} P^{\rm EGM}_{nm}(x)/\mathrm {d} x$, whose value depends on the normalisation of $P^{\rm EGM}_{nm}(x)$. To proceed, denote a generic choice of normalisation by $\widetilde {P}_{nm}$:

(3.5)\begin{equation} \widetilde{P}_{nm}(x) \equiv \alpha_{nm}{K}_{nm}(x), \end{equation}

for some normalisation $\alpha _{nm}$. We will calculate $\mathrm {d}\widetilde {P}_{nm}(\cos \theta )/\mathrm {d}\theta$, and then set $\alpha _{nm}$ to be the EGM normalisation in Equation (2.32). First, call on Equation (2.8) to write

(3.6)\begin{align} \frac{\mathrm{d}}{\mathrm{d} x}\widetilde{P}_{nm}(x) & = \alpha_{nm}\frac{\mathrm{d}}{\mathrm{d} x}{K}_{nm}(x)\nonumber\\[1ex] & = \alpha_{nm}\frac{(n-m+1){K}_{n+1,m}(x)-(n+1)\,x\,{K}_{nm}(x)}{x^{2}-1}\nonumber\\[1ex] & = \frac{\displaystyle\frac{\alpha_{nm}}{\alpha_{n+1,m}}(n-m+1)\alpha_{n+1,m}{K}_{n+1,m}(x)-(n+1)\,x\, \alpha_{nm}{K}_{nm}(x)}{x^{2}-1}\nonumber\\[1ex] & =\frac{\displaystyle\frac{\alpha_{nm}}{\alpha_{n+1,m}}(n-m+1)\widetilde{P}_{n+1,m}(x)-(n+1)\,x\, \widetilde{P}_{nm}(x)} {x^{2}-1}. \end{align}

It follows that

(3.7)\begin{align} \frac{\mathrm{d}}{\mathrm{d}\theta}\widetilde{P}_{nm}(\cos\theta) & = \left.\frac{\mathrm{d}\widetilde{P}_{nm}(x)}{\mathrm{d} x}\right|_{x\,=\,\cos\theta}\!\!\!\!\!\!\!\!\!\!\!\times{-}\!\sin\theta\nonumber\\[1ex] & = \frac{\displaystyle\frac{\alpha_{nm}}{\alpha_{n+1,m}}(n-m+1)\widetilde{P}_{n+1,m}(\cos\theta) - (n+1)\cos\theta\,\widetilde{P}_{nm}(\cos\theta)}{\sin\theta}. \end{align}

In particular, what is Equation (3.7) when $\widetilde {P}_{nm}$ is chosen to be $P^{\rm EGM}_{nm}$? The EGM2008 value of $\alpha _{nm}$ is, from Equation (2.32),

(3.8)\begin{equation} \alpha^{\rm EGM}_{nm} \equiv \sqrt{\frac{k_m(2n+1)(n-m)\textit{!}}{(n+m)\textit{!}}},\quad \text{where } k_m \equiv \begin{cases} 1 & :m = 0,\\ 2 & :m \not= 0. \end{cases} \end{equation}

Hence, the relevant term in Equation (3.7) becomes

(3.9)\begin{align} \frac{\alpha^{\rm EGM}_{nm}}{\alpha^{\rm EGM}_{n+1,m}}(n-m+1) & =\sqrt{\frac{k_m(2n+1)(n-m)\textit{!}}{(n+m)\textit{!}}} \sqrt{\frac{(n+1+m)\textit{!}}{k_m(2n+3)(n+1-m)\textit{!}}} \;\;(n-m+1)\nonumber\\[1ex] & = \sqrt{\frac{(n+1/2)(n+1+m)(n+1-m)}{n+3/2}} \equiv \beta_{nm}. \end{align}

We can now return to Equation (3.7) to write, with $\widetilde {P}_{nm}=P^{\rm EGM}_{nm}$,

(3.10)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}\theta}P^{\rm EGM}_{nm}(\cos\theta) = \frac{\beta_{nm}P^{\rm EGM}_{n+1,m}(\cos\theta) - (n+1)\cos\theta \,P^{\rm EGM}_{nm}(\cos\theta)}{\sin\theta}. \end{equation}

This is our sought-after derivative of EGM2008's Legendre functions. Finally, it follows from Equation (2.1) that

(3.11)\begin{align} \frac{\partial\Phi}{\partial\theta} & = \frac{-GM}{r\sin\theta}\sum_n \left(\frac{a} {r}\right)^{n} \sum_m [\beta_{nm} P^{\rm EGM}_{n+1,m}(\cos\theta) - (n+1)\cos\theta\, P^{\rm EGM}_{nm}(\cos\theta)]\nonumber\\ & \quad \times (C^{\rm EGM}_{nm}\cos m\phi + S^{\rm EGM}_{nm}\sin m\phi). \end{align}

We can now insert the partial derivatives in Equations (3.3), (3.4), (3.11) into Equation (3.1) to extract the sought-after components in Equation (3.2). The geocentric up component of ${\boldsymbol {g}}$ is

(3.12)\begin{equation} {g_r} = \frac{-GM}{r^{2}}\left[1 + \sum_n(1+n)\left(\frac{a}{r}\right)^{n}\sum_m P^{\rm EGM}_{nm}(\cos\theta) (C^{\rm EGM}_{nm}\cos m\phi + S^{\rm EGM}_{nm}\sin m\phi)\right]. \end{equation}

The geocentric south component of ${\boldsymbol {g}}$ is

(3.13)\begin{align} {g_\theta} & = \frac{GM}{r^{2}\sin\theta}\sum_n \left(\frac{a}{r}\right)^{n} \sum_m \left[\beta_{nm}P^{\rm EGM}_{n+1,m}(\cos\theta) - (n+1)\cos\theta\,P^{\rm EGM}_{nm}(\cos\theta)\right]\nonumber\\ & \quad \times (C^{\rm EGM}_{nm}\cos m\phi + S^{\rm EGM}_{nm}\sin m\phi). \end{align}

The east component of ${\boldsymbol {g}}$ is

(3.14)\begin{equation} {g_\phi} = \frac{GM}{r^{2}\sin\theta}\sum_n\left(\frac{a}{r}\right)^{n}\sum_m m P^{\rm EGM}_{nm}(\cos\theta) \left(\mathord{-}C^{\rm EGM}_{nm}\sin m\phi + S^{\rm EGM}_{nm}\cos m\phi\right). \end{equation}

How do these expressions compare with those produced by others? Examples in the literature are difficult to find. The analytical expressions above for our ${g_r}$ and ${g_\phi }$ equal their counterparts $G_r, G_\ell$ listed in equation (10) of Wu et al. (Reference Wu2016). Our ${g_\theta }$ superficially resembles its counterpart $G_\phi$ in their equation (10), but Wu et al. do not write out the derivatives of their associated Legendre functions with respect to co-latitude [our Equation (3.10) above]; unfortunately, it is precisely that quantity which needs the careful attention that we gave it in this section when we produced the non-trivial Equation (3.10). It is not clear if Wu et al. actually use their expressions for $G_r, G_\ell, G_\phi$. Wang et al. (Reference Wang2016) do something similar, by writing ‘$\mathrm {d} P_{nm}(\cos \theta )/\mathrm {d}\theta \,$’ in their equation (5) without giving an explicit expression for that derivative.

Chatfield's textbook (Chatfield, Reference Chatfield1997) gives analytical expressions for its equivalents of the components ${g_r}, {g_\theta }, {g_\phi }$. Chatfield defines his $P_n(x)$ recursively in his equation (7.18); comparing that with our Equations (2.2), (2.3), (2.7) and noting his equation (7.21), it is clear that he has

(3.15)\begin{equation} P_n(x) \equiv {K}_n(x),\quad P_{nm}(x) \equiv {K}_{nm}(x). \end{equation}

This contradicts his appendix D, which clearly uses $P_{nm}(x) \equiv (-1)^{m}{K}_{nm}(x)$ without explicitly saying so. Chatfield refers to p. 328 of a NASA report (NASA SP-365, 1977). Although that NASA report uses somewhat different language to define its Legendre functions, its functions $P_n^{m}(x)$ appear to equal $K_{nm}(x)$. That report gives a derivative of its $P_n^{m}(x)$ that Chatfield places into his equation (7.35). This is well and good, but Chatfield's equations do not use the EGM2008 normalisation and so cannot be used directly with the EGM2008 coefficients. These examples suggest that while spherical harmonics are treated as tools of the trade by geodesy authors, their correct use is far from clear.

On a side note, applying Equations (3.1) and (3.2) to the $J_2$ potential in Equation (2.38) yields the $J_2$ approximation of ${\boldsymbol {g}}$:

(3.16)\begin{align} {g_r}(J_2) & = \frac{-GM}{r^{2}}\left[1 + \frac{3J_2}{2}\left(\frac{a}{r}\right)^{\!2}(1-3\cos^{2}\theta)\right],\nonumber\\ {g_\theta}(J_2) & = \frac{GM}{r^{2}} 3J_2 \left(\frac{a}{r}\right)^{\!2} \sin\theta\cos\theta,\nonumber\\ {g_\phi}(J_2) & = 0. \end{align}

Of course, these components also result from evaluating Equations (3.12)–(3.14) for $n=2$, $m=0$.

In practice, a location's geodetic latitude $\lambda$ and height $h$ above the WGS-84 ellipsoid are typically provided in place of its geocentric co-latitude $\theta$ of Equation (2.1). The conversion from $\lambda$ to $\theta$ can be done in two steps. First, construct the radius vector ${\boldsymbol {r}}$ pointing from the ellipsoid's centre to the point at lat/long/height $(\lambda, \phi, h)$, expressed in ECEF Cartesian coordinates $(x,y,z)$ via the well-known transform:

(3.17)\begin{align} \genfrac{\{}{\}}{0pt}{}{r_x}{r_y} &= \left[\frac{a^{2}}{\sqrt{a^{2}\cos^{2}\lambda+b^{2}\sin^{2}\lambda}}+h\right]\cos\lambda \genfrac{\{}{\}}{0pt}{}{\cos\phi}{\sin\phi},\nonumber\\ r_z &= \left[\frac{b^{2}}{\sqrt{a^{2}\cos^{2}\lambda+b^{2}\sin^{2}\lambda}}+h\right]\sin\lambda, \end{align}

where $a$ and $b$ were listed just after Equation (1.1). Because of the WGS-84 ellipse's rotational symmetry, it suffices to set $\phi =0$ in Equation (3.17). The value of $r$ in Equation (2.1) is

(3.18)\begin{equation} r = \left\vert {{\boldsymbol{r}}} \right\vert = \sqrt{r_x^{2} + r_y^{2} + r_z^{2}}. \end{equation}

The geocentric co-latitude $\theta$ is the angle between ${\boldsymbol {r}}$ and the ECEF's $z$ axis:

(3.19)\begin{equation} \theta = \cos^{{-}1} \frac{r_z}{r}. \end{equation}

The geodetic north-east-down (NED) components of gravity, $g_N, g_E, g_D$, tend to be required rather than the geocentric components $g_r, g_\theta, g_\phi$ in Equations (3.12)–(3.14). Use the following, which can be derived by analysing a sequence of rotations that maps a copy of the normalised spherical-polar basis vectors onto the NED basis vectors:

(3.20)\begin{equation} \begin{bmatrix}\,{g_N}\,\\ \,{g_E}\,\\ \,{g_D}\,\end{bmatrix} =\begin{bmatrix} \phantom{+}\cos(\theta+\lambda) & -\!\sin(\theta+\lambda) & 0\\ 0 & 0 & 1\\ -\!\sin(\theta+\lambda) & -\!\cos(\theta+\lambda) & 0 \end{bmatrix} \begin{bmatrix}\,{g_r}\,\\ \,{g_\theta}\,\\ \,{g_\phi}\,\end{bmatrix}. \end{equation}

How various authors treat a point's geocentric (co-)latitude can be difficult to establish. The angle is useful in the context of spherical-polar coordinates and hence must be defined as above, independent of the ellipsoid. The EGM2008 model does this, as in figure 4.1 of NGA-WGS84 (2014) and figure 5.6 of Hofmann–Wellenhof and Moritz (Reference Hofmann–Wellenhof and Moritz2006). In contrast, Siouris (Reference Siouris2009) ignores the point's height when defining geocentric latitude – despite this then being inconsistent with his definition of the deviation from the normal in his equation (6), as the geometry in his figure 4 shows. The code in the geodesy package GeographicLib (Karney, 2020) certainly uses a correct equation analogous to Equation (3.17), despite its documentation giving an expression for what it calls geocentric latitude ‘$\theta$’ that contains no height above the ellipsoid and which certainly agrees with Equation (3.17) when $h=0$, as it should. A comment in its code calls that angle $\psi$ in at least one place, even though $\psi$ represents a different type of latitude in the documentation. The bottom line is that one must take care to ascertain what individual authors are doing.

4. Sensitivity to Legendre degree cut-off

Approximating the sum in Equation (2.1) requires it to be truncated at some maximum value of $n$, denoted $n_{\text {max}}$. It is instructive to plot the estimates of $\Phi$ and ${\boldsymbol {g}}$ returned by Equation (2.1) as a function of $n_{\text {max}}$, from $n_{\text {max}}=2,3,4,\ldots, 2190$, the largest value included in the EGM2008 database. In other words, we wish to examine the behaviour of the partial sums of Equation (2.1), as increasingly more of its terms are summed. I have used the ‘tide-free’ version of the EGM2008 database in this paper.

Begin with an arbitrary geodetic lat/long/height of $(\lambda, \phi, h) = (30^{\circ }, 30^{\circ }, 0\,\text {m})$ on the WGS-84 ellipsoid, and compute the corresponding value of $r$ and the geocentric co-latitude $\theta$ from Equations (3.17)–(3.19). The $P^{\rm EGM}_{nm}(\cos \theta )$ are provided by Equation (2.33), using the $\widehat {{K}}_{nm}(\cos \theta )$ returned by Matlab's legendre function, as mentioned shortly after Equation (2.14).

Figure 1(a) shows Equation (2.1)'s calculated values of $\Phi$ versus $n_{\text {max}}$. Note that each plot point was generated from the EGM2008 model, whereas the continuous jagged line joining the points is drawn only to guide the eye. As expected, the value of $\Phi$ fluctuates initially, but settles down as $n_{\text {max}}$ increases: it seems that $n_{\text {max}}$ should be chosen to be at least 100.

Figure 1. (a) Estimates of potential $\Phi$ found from summing Equation (2.1) to some largest value of $n$ called $n_{\text {max}}$, at geodetic lat/long/height 30$^{\circ }$, 30$^{\circ }$, 0 m. The estimates of $\Phi$ are a discrete set, indicated by the points, and the solid line is drawn only to guide the eye. (b) Deviation of each estimated component of ${\boldsymbol {g}}$ from its average over all 2189 values of $n_{\text {max}}$. Again, the solid lines are drawn only to guide the eye

Next, we require to plot ${\boldsymbol {g}}$, or rather its geodetic north-east-down components from Equations (3.20) and (3.12)–(3.14). For example, when $n_{\text {max}} = 100$, these estimates areFootnote ²

(4.1)\begin{equation} [\;g_N,\; g_E,\; g_D] \simeq [\;0{\cdot}014695,\;\; 4{\cdot}57208\times {10^{{-}5}},\;\; 9{\cdot}81866\;]\;\; \text{N/kg}. \end{equation}

To compare the behaviour of these estimated components as functions of $n_{\text {max}}$ on one graph, we do the following. Calculate 2189 estimates of $g_N$, one for each $n_{\text {max}}$, from $n_{\text {max}} = 2$ to $n_{\text {max}} = 2190$. Find the mean of these 2189 estimates and subtract that mean from each of the estimates to form a set of 2189 deviations $\Delta g_N$ from the mean. Now plot those deviations versus $\log _{10} n_{\text {max}}$, in red. Do similarly for $\Delta g_E, \Delta g_D$ in green and blue, respectively. The resulting deviations are plotted in Figure 1(b). It is apparent that the estimate of ${\boldsymbol {g}}$ does not settle down until approximately $n_{\text {max}} = 300$. This does not conflict with the fact that the estimate of $\Phi$ settles down for smaller $n_{\text {max}}$, since ${\boldsymbol {g}}$ is determined not only by $\Phi$ at the location of interest, but also by $\Phi$ at nearby locations, and all of these estimates of $\Phi$ are fluctuating as functions of $n_{\text {max}}$. Figures 2–4 repeat these calculations for a sample of other locations.

Figure 2. Similar to Figure 1, but for geodetic lat/long/height 20$^{\circ }$, 40$^{\circ }$, 3000 m

Figure 3. Similar to Figure 1, but for geodetic lat/long/height $-\textit{50}^{\circ }$, 100$^{\circ }$, 200 m

Figure 4. Similar to Figure 1, but for geodetic lat/long/height $\textit{11}{\cdot }\textit{35}^{\circ }$, $\textit{142}{\cdot }\textit{2}^{\circ }$, 0 m, near the Marianas Trench

The fluctuation of the partial sums of Equation (2.1) is evident in these plots. This reflects the fact that using higher values of $n$ (and $m$) in Equation (2.1) introduces spherical harmonics that model shorter-wavelength details in the potential. Because the partial sums are seen to fluctuate at the level of approximately $10^{-4}$ N/kg, we cannot expect that if values of $n$ greater than 2190 were present in the EGM2008 database, the fluctuations would suddenly vanish at those higher values – after all, the number 2190 is not special in this model. We expect similar-sized fluctuations to be present for at least some number of higher-$n$ partial sums. Hence, this limits the precision to which gravity can be determined from the model to approximately $10^{-4}$ N/kg.

This lack of settling down with higher values of $n_{\text {max}}$ is to be expected: the EGM2008 coefficients are presumably a global fit of a series of spherical harmonics to gravity measurements. This fit is determined by minimising some cost function, such as the sum of squares of the measurements’ deviations from the series. We cannot expect such a fit to improve simultaneously at all locations as $n_{\text {max}}$ increases.  It might well improve in some particular location at the expense of others, with the same sort of non-trivial behaviour that we see in the convergence of Fourier series. Hence, there is no fundamental reason to prefer the largest value of $n$ provided by EGM2008 (i.e., 2190) to calculate ${\boldsymbol {g}}$ for any given location. The subject of the numerics of the many computations involved in Equations (3.12)–(3.14) is an old one and a large one. Alternative approaches to ours of working with the associated Legendre functions exist, such as described by Holmes and Featherstone (Reference Holmes and Featherstone2002). It is not clear whether any of these alternative algorithms reduce the fluctuations appearing in Figures 1 to 4 by being numerically better behaved. In any given algorithm, these fluctuations set the limits to our ability to estimate ${\boldsymbol {g}}$, irrespective of whether the algorithm computes exactly to some large number of decimal places. So, we expect that the value of gravity generated from EGM2008 is trustworthy to no more than four decimal places.

5. Comparison of our analysis with other work

How closely do our estimates of ${\boldsymbol {g}}$ match those generated by others from the same EGM2008 database? Rather than estimate the physicist's gravity/gravitation ${\boldsymbol {g}}$, geodesists usually estimate what they call ‘gravity’ (and which they also denote ${\boldsymbol {g}}$), which they define to be minus the support force per unit mass on a body that is held at a fixed altitude. I will call this ‘apparent gravity’ ${\boldsymbol {g}}^{a}$, being the downward force we feel when supported at a fixed altitude. (One might also call it ‘specific weight’, where specific denotes ‘per unit mass’.) For example, if Earth spun once every 84${\cdot }$5 minutes in an inertial frame, then every object on its equator would have the same speed in that frame as a surface-skimming satellite. The apparent gravity on the equator would then be zero, since the object would hover over the ground, needing no support and hence feeling no weight.

To see how apparent gravity relates to the ${\boldsymbol {g}}$ that we calculated in Equation (3.2), consider a mass $m$ fixed to the geoid. In an inertial frame, two forces act on the mass: the actual gravity force $m{\boldsymbol {g}}$ pulling it down, and the support force ${\boldsymbol {S}}$ that holds it at a fixed altitude. These forces sum to give ‘mass times acceleration’ in the inertial frame in which Earth spins with angular velocity ${\boldsymbol {\Omega }}$. The mass's acceleration in that inertial frame is ${\boldsymbol {\Omega }}\times ({\boldsymbol {\Omega }}\times {\boldsymbol {r}})$, where ${\boldsymbol {r}}$ is the position vector of the mass relative to any point on Earth's spin axis; Earth's centre is usually chosen as that point. Hence ‘force equals mass times acceleration’ in the inertial frame becomes

(5.1)\begin{equation} m{\boldsymbol{g}} + {\boldsymbol{S}} = m {\boldsymbol{\Omega}}\times({\boldsymbol{\Omega}}\times{\boldsymbol{r}}). \end{equation}

The apparent gravity is minus the support force ${\boldsymbol {S}}$ per unit mass:

(5.2)\begin{equation} \text{apparent gravity } {\boldsymbol{g}}^{a} \equiv{-}{\boldsymbol{S}}/m = {\boldsymbol{g}} - {\boldsymbol{\Omega}}\times({\boldsymbol{\Omega}}\times{\boldsymbol{r}}). \end{equation}

Set Earth's spin rate to be the WGS-84 value of $\Omega = 7{\cdot }29211\;58553 \times {10^{-5}}$ rad/s as in table 3.4 of NGA-WGS84 (2014), and ${\boldsymbol {r}}$ to be the vector from Earth's centre to the point of interest, where Earth is modelled as the WGS-84 ellipsoid.

In Table 1, we compare our results for the north-east-down components of ${\boldsymbol {g}}^{a}$ at six locations with those of the geodesy package Tracker Component Library (US Naval Research Laboratory, 2021), quoting as many decimal places as needed to show where the package's estimates differ from our own. Both we and Tracker Component Library use the full sum of 2190 terms in the EGM2008 ‘tide-free’ database. Throughout this paper, I have used Matlab's legendre command to calculate the associated Legendre functions. Tracker Component Library uses the approach of Holmes and Featherstone (Reference Holmes and Featherstone2002) to calculate these functions. The results in Table 1 agree to typically one part in $10^{9}$.

Table 1. A comparison of our results with those of the Tracker Component Library software package. North-east-down coordinates are used, with units N/kg. Sufficient digits are quoted to make a clear comparison

For a partial comparison with the values predicted by the above theory, the International Gravity formula (‘IG’) (Moritz, Reference Moritz2000; Li and Götze, Reference Li and Götze2001) returns a best fit of the magnitude of apparent gravity on the GRS-80 ellipsoid (Featherstone and Dentith, Reference Featherstone and Dentith1997). It takes geodetic latitude $\lambda$ as its sole input:Footnote ³

(5.3)\begin{equation} \left\vert {{\boldsymbol{g}}^{a}} \right\vert_{\text{IG}} = 9{\cdot}780327 (1 + 5{\cdot}3024\times {10^{{-}3}}\sin^{2}\lambda - 5{\cdot}8\times {10^{{-}6}}\sin^{2} 2\lambda) \text{ N/kg}. \end{equation}

Table 2 compares our values of apparent gravity from Equation (5.2) for some arbitrary locations and using the full sum to $n=2190$, with Equation (5.3) (which, being for the GRS-80 ellipsoid, only permits a slightly imperfect comparison). We have also included the prediction of the $J_2$ model, Equation (3.16). For this, recall that it requires geocentric co-latitude $\theta$, not geodetic latitude $\lambda$, and that the WGS-84 value of $J_2$ is approximately $1{\cdot }08263 \times {10^{-3}}$.

Table 2. Estimates of apparent gravity (N/kg) at some arbitrary locations (latitude is geodetic), specified on the WGS-84 ellipsoid. The IG values are actually for the GRS-80 ellipsoid, but serve as a comparison

6. Concluding comments

Estimating Earth's gravity from the EGM2008 database requires intricate preliminary work. The following list of conclusions are evident.

• • Extracting values of gravity from EGM2008 is quick and straightforward if one is able to calculate the associated Legendre functions. But the exercise requires attention to the obscure details of the model, such as I have laid out in this paper. For example, due to a non-standard choice of associated-Legendre ‘normalisation’ that is now embedded in geodesy literature, care is needed when computing the associated Legendre functions used in EGM2008.

• • Because EGM2008 is a global fit of spherical harmonics to Earth's potential in Equation (1.1), we can never be sure that summing some given number of terms in Equation (1.1) will give a better estimate at any given point than summing a smaller number of terms. This is akin to the idea that when a given function $f(x)$ is approximated by, say, a 10-term Fourier series, then at any particular value $x=x_0$ at which we are interpolating the function, the sum of the first 5 terms might be closer to $f(x_0)$ than the sum of the first 10 terms. Because of that, summing the full EGM2008 set will not necessarily give a better value of ${\boldsymbol {g}}$ than summing a smaller number of terms.

• • Estimates of gravity found from summing even the complete EGM2008 database – with its almost five million coefficients, each specified to 15 decimal places – should not be accorded any precision finer than $10^{-4}$ N/kg.

On the second conclusion above, it might well be that in a small region of interest, an appropriate approach would measure a reliable value of gravity in that region's centre, and then determine the value of $n_{\text {max}}$ that returns a gravity value closest to this number. That value of $n_{\text {max}}$ might be much less than 2190. We might then use that $n_{\text {max}}$ for all calculations of ${\boldsymbol {g}}$ in that region, rather than always using $n_{\text {max}}=2190$.

The second and third conclusions above apply to any model that fits a single smooth function to a set of measurements. Also, with any such fit – whether a sum of spherical harmonics in three dimensions or a simple polynomial in one dimension – even a change to a single measurement will usually alter the value of the fitting function at every other point. In other words, if the gravity at Paris had been used as a measurement for constructing the EGM2008 database, and the gravity at Adelaide had not been used as a measurement, then an improvement in the gravity measurement at Paris would affect the gravity computed at Adelaide.

Perhaps that situation is less than optimal. When attention is confined to more localised regions, other approaches come into their own, such as the use of other functions discussed by Li (Reference Li2018) and his references.

Hypothetically, if measured values of gravity and/or potential were stored for a grid of locations, then might these values be interpolated in some small region local to the point of interest, independently of any changes made to the database in remote locations?

The problem with creating such a grid is that it must be three dimensional. Consider a grid in just two dimensions: say, for zero height above the ellipsoid. What resolution would such a grid have, if it were created from a set of gravity vectors that is as large as the EGM2008 database? The EGM2008 text file of coefficients has a size of about 240 megabytes, and is not claimed to be accurate for heights that depart appreciably from zero. Approximately half of the file is error terms that probably no one uses, and so the useful data occupies roughly 120 megabytes. Suppose this number of bytes were used to store individual values of the gravity vector for a spread of locations on Earth's surface. If each gravity vector occupied, say, 20 bytes, then the number of vectors that could be encoded with 120 megabytes would be 6 million. Earth has a surface area of $4\pi (6370\;\text {km})^{2}$. If this were divided into 6 million square tiles, each tile would have a side length of 9 km. That grid spacing could be acceptable in some applications. And, because each grid point holds the actual value of gravity measured at that point, this set of points could be interpolated to estimate gravity at intermediate points to a high accuracy. Or, rather than store gravity values, we might store potentials. If each potential occupies 10 bytes, the 120 megabytes can store potentials for 12 million locations. Each of the resulting 12 million square tiles will have a side length of 6${\cdot }$5 km.

Extending this grid idea to three dimensions greatly increases the required file space. Nonetheless, whatever is stored, an interpolation will need only a handful of calculations and a handful of data points. Contrast that with the probably hundreds of millions of multiplications, and the five million coefficients, that are currently involved in implementing EGM2008 at each point of interest. In practice, gravity values for a specific application that requires to estimate gravity in real time are probably best computed in advance, and any interpolations can be made from those pre-computed values.

One point worth remembering is that perhaps most systems that require gravity values have a limit on the use that they can make of those values. For example, the limiting factor for inertial navigation systems is not typically a knowledge of gravity, but of the vehicle's spatial orientation, which is required to subtract gravity's effect on the inertial measurement unit. These systems are typically unable to benefit from anything more accurate than a zeroth- or first-order gravity model.

EGM2008 is often held up as the premier model of Earth's gravity, to be called on to supply ultra-accurate estimates of gravity. But some knowledge of its accuracy is necessary; or, in the absence of that, we might be content with knowing its precision. It is clear that even in an idealised ‘tide-free’ environment, we can expect at most a $10^{-4}$ N/kg precision from EGM2008, at the cost of processing perhaps a hundred megabytes of data for each estimate made. In a given application, even this precision can be swamped by errors in equipment calibration, vehicle orientation and environmental noise. We might add to that list the difficulties of implementation caused by the model's non-standard normalisation of its associated Legendre functions.

I thank Stephen Thompson of the Defence Science and Technology Group for helping with the analysis of the examples, and Tim Sills of the Defence Science and Technology Group for discussions of the content of this paper. I also thank the referees for their constructive comments.