Continuation of Constant Solution Spherical Harmonic

Abstract

A new analytical method for the computation of a truncated series of solid spherical harmonic coefficients (HCs) from data on a spheroid (i.e. an oblate ellipsoid of revolution) is derived, using a transformation between surface and solid spherical HCs. A two-step procedure is derived to extend this transformation beyond degree and order (d/o) 520. The method is compared to the Hotine–Jekeli transformation in a numerical study based on the EGM2008 global gravity model. Both methods are shown to achieve submicrometre precision in terms of height anomalies for a model to d/o 2239. However, both methods result in spherical harmonic models that are different by up to 7.6 mm in height anomalies and 2.5 mGal in gravity disturbances due to the different coordinate system used. While the Hotine–Jekeli transformation requires the use of an ellipsoidal coordinate system, the new method uses only spherical polar coordinates. The Hotine–Jekeli transformation is numerically more efficient, but the new method can more easily be extended to cases where (a linear combination of) normal derivatives of the function under consideration are given on the surface of the spheroid. It therefore provides a solution to many types of ellipsoidal boundary-value problems in the spectral domain.

1 INTRODUCTION

Spherical harmonic expansions are used in many branches of Earth and planetary science to construct global models of various quantities of interest. They are especially useful for the representation of harmonic functions. For example, global models of the Earth's gravity field are almost exclusively represented by a series of spherical harmonic coefficients (HCs). The International Centre for Global Earth Models currently lists more than 150 spherical harmonic models of the Earth's gravitational potential, as well as spherical harmonic models of the gravitational fields of the Moon, Mars and Venus. Spherical harmonic expansions are also used to describe the geomagnetic field (e.g. Cain et al. 1989; Maus et al. 2005; Lesur et al. 2008; Chulliat et al. 2015). Other applications include heat flow models (e.g. Chapman & Pollack 1980; Pollack et al. 1993; Hofmeister & Criss 2005; Hamza et al. 2008), isostatic-topographic models (e.g. Rummel et al. 1988; Panasyuk & Hager 2000; Kaban et al. 2004; Balmino et al. 2012; Claessens & Hirt 2013), seismic wave speed models (e.g. Su et al. 1994) and post-seismic Earth surface deformation models (e.g. Pollitz 1996; Riva & Vermeersen 2002).

Spherical HCs can easily be computed if data are available on the surface of a sphere that is completely within the region where the function is harmonic. However, since the Earth (as well as many other bodies in the solar system) is to a much higher level of accuracy approximated by a spheroid (i.e. an oblate ellipsoid of revolution), the need arises for a methodology to compute spherical HCs from data on a spheroid.

One possible methodology is a numerical one, using a statistical technique such as least-squares estimation or collocation (e.g. Sansò & Tscherning 2003; Abd-Elmotaal et al. 2014), but here an analytical solution is sought. Analytical methodologies have the advantage that they allow for computation of coefficients above the maximum degree defined by the resolution of the data grid, and that they can avoid problems associated with the inversion of an ill-conditioned matrix that burden the least-squares estimation procedure.

Several analytical methodologies have been proposed, such as the upward-continuation method (Cruz 1986) and a method based on Green's second integral theorem (Sjöberg 1988). Both these methods do not only require the function values on the spheroid, but also their radial or normal derivatives, and they are now rarely used in practice if at all.

The most popular method is here called the spheroidal harmonics method. In this method, the function on the surface of the spheroid is first expanded into a series of spheroidal harmonics, which can then be transformed into a spherical harmonic expansion using a transformation between spherical and spheroidal HCs (e.g. Buchdahl et al. 1977; Jekeli 1988; Dechambre & Scheeres 2002). The Hotine–Jekeli transformation between spherical and spheroidal HCs (Jekeli 1988) was used, for example, in the creation of the global gravity model EGM2008 (Pavlis et al. 2012) and the global model of the Earth's lithospheric magnetic field NGDC-720 (Maus 2010).

An alternative method is here called the surface harmonics method. It follows from a transformation between solid spherical HCs and surface spherical HCs (Claessens 2006; Claessens & Featherstone 2008). Like the spheroidal harmonics method, it does not require any other information than the function values on the spheroidal surface. Claessens & Featherstone (2008) show that an efficient transformation from surface to solid spherical HCs can only be achieved up to degree and order (d/o) ∼520, which severely hampers the surface harmonics method. However, in Section 4, it will be shown that an efficient transformation can be achieved up to ultrahigh d/o using a new two-step procedure.

Apart from this improvement to the surface harmonics method, this paper also provides a numerical comparison between the spheroidal harmonics method and the surface harmonics method. The major advantage of the surface harmonics method is that it can easily be extended to cases where the function itself is not given on the surface of the spheroid, but (a linear combination of) its surface normal derivatives. This provides an advantage, for example, when computing spherical HCs of the gravitational potential from gravity disturbances or gravity anomalies (see Section 6). An additional but only very minor advantage of the surface harmonics method over the spheroidal harmonics method is that the data are completely treated in spherical (polar) coordinates. Spheroidal coordinates (Jacobi ellipsoidal coordinates) that are used in spheroidal harmonics are not required. Note that many global data sets are provided in terms of geodetic coordinates (Gauss ellipsoidal coordinates), but these coordinates cannot directly be used for harmonic analysis or synthesis, so in these cases a re-parametrization of the data is always required. See Gruber & Abrykosov (2014) for a further discussion.

The determination of spherical HCs from function values on a spheroid is a solution to an ellipsoidal boundary-value problem. Several types of boundary-value problems can be distinguished, based on the form of the data provided. When the function itself is known on the spheroid, a so-called Dirichlet boundary-value problem arises, and the solution of this case is the main focus of this paper. However, in Section 6, it is shown that other types of boundary-value problems, such as the Neumann, Robin and second-order (gradiometric) boundary-value problem (e.g. Mackie 1965) can also be solved using the surface harmonics method.

2 SPHERICAL AND SPHEROIDAL HARMONICS

The series of spherical harmonics provides a solution to Laplace's differential equation, which all harmonic functions obey (e.g. Hobson 1931; Sigl 1985). Since the Earth's gravitational and magnetic fields are harmonic vector fields exterior to the masses of the Earth, they can be represented anywhere in the exterior space by a solid spherical harmonic expansion

\begin{equation} f(r,\theta ,\lambda ) = \sum _{n=0}^{\infty } \bigg ( \frac{R}{r} \bigg )^{n+1} \sum _{m=-n}^n \overline{f}_{nm}^R \overline{Y}_{nm} (\theta , \lambda ) \end{equation}

(1)

where f is a function that is harmonic external to some closed domain, (r, θ, λ) is the set of spherical polar coordinates, R is the radius of a reference sphere, |$\overline{f}_{nm}^R$| are the fully normalized solid spherical HCs of degree n and order m and |$\overline{Y}_{nm}$| are the fully normalized spherical harmonic functions. The spherical harmonic functions form a set of orthogonal base functions on the sphere, and are defined by

\begin{equation} \overline{Y}_{nm} (\theta , \lambda ) = \overline{P}_{n|m|} (\cos \theta ) \bigg \lbrace \begin{array}{lll}\cos m\lambda &\quad \mbox{for} &\quad m \le 0 \\ \sin m\lambda &\quad \mbox{for} &\quad m > 0 \end{array} \end{equation}

(2)

where |$\overline{P}_{n|m|} (\cos \theta )$| are fully normalized (4π-normalized) associated Legendre functions. These form the solution to the associated Legendre differential equation (e.g. Abramowitz & Stegun 1972) and numerically stable algorithms to compute them up to high d/o are provided in, for example, Masters & Richards-Dinger (1998), Holmes & Featherstone (2002), Jekeli et al. (2007) and Fukushima (2012).

A more general alternative to the spherical harmonic expansion is an ellipsoidal harmonic expansion. Since the geodetic reference ellipsoid is a spheroid, this paper is only concerned with the spheroidal harmonic expansion, which is a special case of the ellipsoidal harmonic expansion for a triaxial ellipsoid (the spherical harmonic expansion in eq. (1) is in turn a special case of the spheroidal harmonic expansion). The spheroidal expansion provides a representation of the function f in spheroidal coordinates (u, β, λ)

\begin{equation} f(u,\beta ,\lambda ) = \sum _{n=0}^{\infty } \sum _{m=-n}^n \frac{Q_{nm} \big ( i \frac{u}{E} \big )}{Q_{nm} \big ( i \frac{b}{E} \big ) } \overline{f}_{nm}^{\mathrm{u}} \overline{Y}_{nm} (\beta , \lambda ) \end{equation}

(3)

where b is the semi-minor axis of the spheroid, E is the linear eccentricity of the spheroid and Q_nm are associated Legendre functions of the second kind. The computation of the associated Legendre functions of the second kind Q_nm is hampered by numerical instabilities (Sona 1995), but Fukushima (2013) provides a stable algorithm for the computation of the ratio of two of these functions.

When the spherical harmonic expansion is synthesized on the reference sphere (r = R), eq. (1) reduces to

\begin{equation} f(R,\theta ,\lambda ) = \sum _{n=0}^{\infty } \sum _{m=-n}^n \overline{f}_{nm}^R \overline{Y}_{nm} (\theta , \lambda ) \end{equation}

(4)

and when the spheroidal harmonic expansion is synthesized on the reference spheroid (u = b), eq. (3) reduced to

\begin{equation} f(u,\beta ,\lambda ) = \sum _{n=0}^{\infty } \sum _{m=-n}^n \overline{f}_{nm}^{\mathrm{u}} \overline{Y}_{nm} (\beta , \lambda ). \end{equation}

(5)

Eqs (4) and (5) represent the inverse Legendre transform (Jekeli 1988). A Legendre transform can be applied to any continuous function, harmonic or not, defined on any surface for which every point is uniquely defined by its latitude and longitude, in a suitable coordinate system. Thus, a function defined on the surface of a spheroid can also be represented by a series of coefficients using spherical polar latitude and longitude

\begin{equation} f(r_{\mathrm{e}},\theta ,\lambda ) = \sum _{n=0}^{\infty } \sum _{m=-n}^n \overline{f}_{nm}^{\mathrm{e}} \overline{Y}_{nm} (\theta , \lambda ) \end{equation}

(6)

where r _e is the spheroidal radius, which for a spheroidal surface can be expressed as a function of the spherical polar latitude θ, the semi-major axis of the spheroid a and the first numerical eccentricity of the spheroid e ²

\begin{equation} r_{\mathrm{e}}=a\sqrt{\frac{1-e^2}{1-e^2 \sin ^2 \theta }} \end{equation}

(7)

The coefficients |$\overline{f}_{nm}^{\mathrm{e}}$| are called surface spherical HCs. They can be computed from an integration over the unit sphere σ, the (direct) Legendre transform

\begin{equation} \overline{f}_{nm}^{\mathrm{e}} = \frac{1}{4 \pi } \int \limits _{\sigma } f(r_{\mathrm{e}}, \theta , \lambda ) \overline{Y}_{nm}(\theta , \lambda ) d\sigma . \end{equation}

(8)

The coefficients |$\overline{f}_{nm}^R$| (eq. 4) and |$\overline{f}_{nm}^{\mathrm{u}}$| (eq. 5) can also be computed via the Legendre transform, but they are called solid spherical HCs and solid spheroidal HCs, respectively, to indicate that they can be used more generally in the solid harmonic expansions in eqs (1) and (3).

While the coefficients |$\overline{f}_{nm}^{\mathrm{e}}$| and |$\overline{f}_{nm}^{\mathrm{u}}$| may represent the same function on the same surface, they will attain different numerical values due to the use of a different coordinate system.

3 SPHEROIDAL HARMONICS METHOD

For the sake of completeness, the well-known formulae of the Hotine–Jekeli transformation between spherical and spheroidal HCs are shown here in the current notation. The forward transformation, that is, the transformation from solid spherical HCs |$\overline{f}_{nm}^R$| to solid spheroidal HCs |$\overline{f}_{nm}^{u}$|⁠, is given by Jekeli (1988)

\begin{equation} \overline{f}_{nm}^u = \overline{S}_{n|m|}\left(\frac{b}{E}\right) \sum _{i=0}^w \delta _{nmi} \overline{f}_{n-2i,m}^R \end{equation}

(9)

where |$\overline{S}_{n|m|}(\frac{b}{E})$| is Jekeli's renormalized Legendre function of the second kind, δ _nmi are weights that are computed recursively (Gleason 1988, eq. 1.18), and

\begin{equation} w = \mathrm{int}\left(\frac{n-m}{2}\right) \end{equation}

(10)

The reverse transformation is given by Jekeli (1988)

\begin{equation} \overline{f}_{nm}^R = \sum _{i=0}^w \frac{\Delta _{nmi}}{\overline{S}_{n|m|}\left(\frac{b}{E}\right)} \overline{f}_{n-2i,m}^u \end{equation}

(11)

where the weights Δ _nmi can also be computed recursively (Gleason 1988, eq. 1.24). Sebera et al. (2012, eq. 20) provide optimised recurrence relations for Jekeli's renormalized Legendre functions of the second kind. The reverse transformation can be used to generate a series of solid spherical HCs from data on a spheroid, after the data on the spheroid has been expanded into solid spheroidal HCs. This is the transformation of most interest here.

4 SURFACE HARMONICS METHOD

4.1 Forward transformation

Claessens & Featherstone (2008) provide a transformation from solid spherical HCs |$\overline{f}_{nm}^R$| to surface spherical HCs |$\overline{f}_{nm}^{\mathrm{e}}$| (the forward transformation). Just as in the spheroidal harmonic method, the reverse transformation is of most interest, because this provides a method to generate a set of solid spherical HCs from data on a spheroid. The reverse transformation is treated in Section 4.2, but the derivation of the forward transformation is also briefly repeated here, primarily to facilitate the subsequent explanation of the reverse transformation.

Assuming that the function f is harmonic everywhere outside and on the surface of the spheroid, the solid and surface spherical harmonic expansions in eqs (1) and (6) can be equated

\begin{equation} \sum _{n=0}^{\infty } \sum _{m=-n}^n \bigg ( \frac{R}{r_{\mathrm{e}}} \bigg )^{n+1} \overline{f}_{nm}^R \overline{Y}_{nm} (\theta , \lambda ) = \sum _{n=0}^{\infty } \sum _{m=-n}^n \overline{f}_{nm}^{\mathrm{e}} \overline{Y}_{nm} (\theta , \lambda ) \end{equation}

(12)

The forward transformation relies on a binomial series expansion of the latitude-dependent term |$(R/ r_{\mathrm{e}})^{n+1}$|

\begin{eqnarray} \bigg ( \frac{R}{r_{\mathrm{e}}} \bigg )^{n+1} = \bigg ( \frac{c}{\sqrt{1-e^2}} \bigg )^{n+1} (1-e^2 \sin ^2 \theta )^{\frac{n+1}{2}} = \sum _{j=0}^{\infty } \alpha _{nj} \sin ^{2j} \theta \!\!\!\nonumber\\ \end{eqnarray}

(13)

where |$c=R/ a$|⁠, that is, the ratio between the reference sphere radius and the semi-major axis of the spheroid, and

\begin{equation} \alpha _{nj} = \bigg ( \frac{c}{\sqrt{1-e^2}} \bigg )^{n+1} (-1)^j {\frac{n+1}{2} \atopwithdelims ()j} e^{2j} \end{equation}

(14)

Note that the summation over j will have only a finite number of non-zero terms for odd n, whereas for even n, the number of non-zero terms is infinite.

The transformation then utilizes a relation among spherical harmonic functions of the same order m (Claessens 2005)

\begin{equation} \sin ^{2j} \theta \overline{Y}_{nm} (\theta , \lambda ) = \sum _{i=-j}^j \overline{K}_{nm}^{2i,2j} \overline{Y}_{n+2i,m} (\theta , \lambda ) \end{equation}

(15)

where the weights |$\overline{K}_{nm}^{2i,2j}$| depend only on the spherical harmonic degree n, order m and the summation indices i and j. The formulae for their computation are shown in Claessens and Featherstone (2008, eqs 18 and 19). Using eqs (13) and (15), the transformation formula from solid to surface spherical HCs takes the form of an infinite weighted summation over solid spherical HCs of equal order m (Claessens & Featherstone 2008)

\begin{equation} \overline{f}_{nm}^{\mathrm{e}} = \sum _{i=-\infty }^{\infty } \lambda _{nmi} \overline{f}_{n-2i,m}^{R} \end{equation}

(16)

where the weights also follow from an infinite summation

\begin{equation} \lambda _{nmi} = \sum _{j=|i|}^{\infty } \alpha _{n-2i,j} \overline{K}_{n-2i,m}^{2i,2j} \end{equation}

(17)

Upon comparison of eqs (9) and (16), it can be seen that both the spheroidal harmonics method and the surface harmonics method include a simple weighted summation over spherical HCs. While eq. (9) contains a finite summation and eq. (16) an infinite summation, both summations are in practice truncated after a number of terms, which is permissible as both are convergent.

4.2 Reverse transformation

The reverse transformation from surface spherical HCs to solid spherical HCs is more complicated, but a solution can be found when it is acknowledged that, in practice, the spherical harmonic series is only evaluated up to a certain maximum d/o M. This means that eq. (16) can be written in matrix form

\begin{equation} \mathcal {F}^{{\mathrm{e}}} = \Lambda \mathcal {F}^{R} \end{equation}

(18)

where |$\mathcal {F}^{{\mathrm{e}}}$| and |$\mathcal {F}^{R}$| are vectors containing all surface and solid HCs up to the maximum d/o, respectively and Λ is a matrix containing the weights λ _nmi . Since Λ is square, eq. (18) can easily be inverted, which results in the reverse transformation formula

\begin{equation} \mathcal {F}^{R} = \Lambda ^{-1} \mathcal {F}^{{\mathrm{e}}} \end{equation}

(19)

The inversion of the matrix Λ faces two potential problems: (1) the matrix may not be well conditioned for a high M and (2) the inversion is computationally intensive for high M. Since a spherical harmonic expansion up to d/o M contains (M + 1)² coefficients, the dimension of Λ is (M + 1)² by (M + 1)². It can be seen from eq. (16) that, if the coefficients in the vectors are first sorted by order m, and subsequently sorted by even and odd degrees n, the matrix Λ becomes block diagonal, which facilitates its inversion. Still, it requires the inversion of 4M matrices of dimension up to |$(M+2)/ 2$| by |$(M+2)/ 2$| for M is even and |$(M+1)/ 2$| by |$(M+1)/ 2$| for M is odd, which is a considerable computational effort for large M. A more efficient method is presented hereafter.

Claessens & Featherstone (2008) have shown that Λ is well conditioned at least up to a certain d/o, thanks to the diagonal dominance of the matrix for low M. For low degrees and orders, the weights λ _nmi are close to 1 for i = 0 and close to zero for i ≠ 0. Since the weights of index i = 0 are all situated on the diagonal of the matrix Λ, for spherical harmonic expansions up to a certain d/o the matrix will be strictly diagonally dominant, that is, for all rows the absolute value of the diagonal element will be larger than the sum of the absolute values of all other elements in the row. This comes down to the condition

\begin{equation} |\lambda _{nm0}| > \sum _{i=-\infty , i \ne 0}^{\infty } |\lambda _{nmi}| \mbox{for} 0 \le m \le n \le M. \end{equation}

(20)

In the case of a spheroid with the eccentricity of the Earth (e ² ≈ 0.00669), the condition in eq. (20) is met for M ≤ 520 (Claessens & Featherstone 2008). Since strictly diagonally dominant matrices are always well conditioned, the inverse of the matrix Λ can be found when this condition is met.

More interestingly, a set of linear equations, such as the one in eq. (18), can be solved using an iterative approach which will always converge if the condition of diagonal dominance is met (e.g. Golub & Van Loan 1996, p. 120). An iterative approach, such as the Jacobi method, Gauss–Seidel method, or successive overrelaxation (SOR) method (e.g. Strang 1986; Golub & Van Loan 1996), is computationally much more efficient than a complete matrix inversion. However, once diagonal dominance is lost, that is, beyond d/o 520, these iterative algorithms quickly diverge.

Nowadays, a maximum spherical harmonic d/o of M = 520 is no longer adequate for many practical applications. Spherical harmonic models of the Earth's gravity field up to higher d/o are already in existence, such as EGM2008 (Pavlis et al. 2012, M = 2190). Also, recent harmonic models of the Earth's magnetic field, such as EMM2015 (Chulliat et al. 2015, M = 720), exceed d/o 520. Therefore, a two-step transformation procedure that overcomes the problems with the inversion of Λ is proposed, as follows.

In the two-step procedure for the reverse transformation, an approximate solution for the solid spherical HCs |$\overline{f}_{nm}^R$| is computed first. The rationale behind this is the following. If the approximate coefficients are a close approximation of the true coefficients, the transformation from approximate to true coefficients potentially gives rise to a diagonally dominant transformation matrix, since the diagonal elements will be close to 1 and the off-diagonal elements will be small.

A good approximation for the solid spherical HCs can be found by approximating each term in the summation on the left-hand side of eq. (12) by each term in the summation on the right-hand side in eq. (12), even though the summation does not strictly hold term-by-term

\begin{equation} \overline{f}_{nm}^R \overline{Y}_{nm} \approx \bigg ( \frac{r_e}{R} \bigg )^{n+1} \overline{f}_{nm}^{\mathrm{e}} \overline{Y}_{nm} \end{equation}

(21)

where the term |$(R/ r_{\mathrm{e}})^{n+1}$| has been moved to the right-hand side of the equation. The term |$(r_{\mathrm{e}}/ R)^{n+1}$| is then expanded into a binomial series similar to eq. (13), and eq. (15) is inserted

\begin{equation} \hat{f}_{nm}^R \overline{Y}_{nm} = \sum _{j=0}^{\infty } \omega _{nj} \sum _{i=-j}^j K_{nm}^{2i,2j} \overline{f}_{nm}^{\mathrm{e}} \overline{Y}_{n+2i,m} \end{equation}

(22)

where the hat in |$\skew4\hat{f}_{nm}^R$| indicates it is an approximation of |$\overline{f}_{nm}^R$|⁠, and

\begin{equation} \omega _{nj} = \bigg ( \frac{\sqrt{1-e^2}}{c} \bigg )^{n+1} (-1)^j {-\frac{n+1}{2} \atopwithdelims ()j} e^{2j}. \end{equation}

(23)

Reorganization of the summation order results in a practical formula for the computation of |$\skew4\hat{f}_{nm}^R$|

\begin{equation} \skew4\hat{f}_{nm}^R = \sum _{i=-\infty }^{\infty } \mu _{nmi} \overline{f}_{n-2i,m}^{\mathrm{e}} \end{equation}

(24)

where the weight factors μ _nmi are computed in a very similar fashion as the weight factors λ _nmi in eq. (17)

\begin{equation} \mu _{nmi} = \sum _{j=|i|}^{\infty } \omega _{n-2i,j} K_{n-2i,m}^{2i,2j}. \end{equation}

(25)

Thus, the first step in the two-step procedure is the computation of the approximate coefficients |$\skew4\hat{f}_{nm}^R$| through eq. (24). In the second step, the 'true' coefficients |$\overline{f}_{nm}^R$| are computed from the approximate coefficients. This second step is a reverse transformation that follows from inserting eq. (16) into eq. (24)

\begin{eqnarray} \skew4\hat{f}_{nm}^R = \sum _{i=-\infty }^{\infty } \mu _{nmi} \sum _{j=-\infty }^{\infty } \lambda _{n-2i,m,j} f_{n-2(i+j),m}^R = \sum _{i=-\infty }^{\infty } \nu _{nmi} f_{n-2i,m}^R \nonumber\\ \end{eqnarray}

(26)

where the weight factors ν _nmi follow from a rearrangement of summation orders

\begin{equation} \nu _{nmi} = \sum _{j=-\infty }^{\infty } \mu _{nmj} \lambda _{n-2j,m,i-j}. \end{equation}

(27)

The inversion of eq. (26) is computationally stable and efficient, since a transformation matrix containing the weights ν _nmi is strongly diagonally dominant up to very high d/o M, that is,

\begin{equation} |\nu _{nm0}| > \sum _{i=-\infty , i \ne 0}^{\infty } |\nu _{nmi}| \mbox{for} 0 \le m \le n \le M \end{equation}

(28)

for very large M. This is shown numerically in Section 5.1.

Due to the diagonal dominance shown by the weights ν _nmi , eq. (26) can be solved iteratively. Various methods can be utilized. Using the Jacobi iteration method (e.g. Strang 1986), the solution after k + 1 iterations is given by

\begin{equation} \overline{f}_{nm(k+1)}^R = \frac{1}{\nu _{nm0}} \bigg ( \hat{f}_{nm}^R - \sum _{i={{\rm int}}(\frac{n-M}{2}),i\ne 0}^{{{\rm int}}(\frac{n}{2})} \nu _{nmi} \overline{f}_{n-2i,m(k)}^R \bigg ) \end{equation}

(29)

where |$\overline{f}_{nm(k)}^R$| are the solid spherical HCs after k iterations. The approximate coefficients can be used as initial values for the iteration: |$\overline{f}_{nm(0)}^R=\skew4\hat{f}_{nm}^R$|⁠.

In the Gauss–Seidel method, the coefficients |$\overline{f}_{nm}^R$| are computed consecutively for increasing degree n, and this method utilizes the fact that when computing |$\overline{f}_{nm}^R$|⁠, the coefficients of lower degree n have already been computed. Eq. (29) is then replaced by

\begin{eqnarray} \overline{f}_{nm(k+1)}^R &=& \frac{1}{\nu _{nm0}} \Bigg ( \hat{f}_{nm}^R - \sum _{i={{\rm int}}(\frac{n-M}{2})}^{-1} \nu _{nmi} \overline{f}_{n-2i,m(k)}^R \nonumber\\ &&- \sum _{i=1}^{{{\rm int}}(\frac{n}{2})} \nu _{nmi} \overline{f}_{n-2i,m(k+1)}^R \Bigg ). \end{eqnarray}

(30)

Faster convergence can be obtained with the SOR method, where an overrelaxation factor ω is introduced. The solution in this case reads (cf. Strang 1986)

\begin{eqnarray} \overline{f}_{nm(k+1)}^R &=& \overline{f}_{nm(k)}^R + \frac{\omega }{\nu _{nm0}} \Bigg ( \hat{f}_{nm}^R - \sum _{i={{\rm int}}(\frac{n-M}{2})}^{0} \nu _{nmi} \overline{f}_{n-2i,m(k)}^R \nonumber\\ &&- \sum _{i=1}^{{{\rm int}}(\frac{n}{2})} \nu _{nmi} \overline{f}_{n-2i,m(k+1)}^R \Bigg ) \end{eqnarray}

(31)

where the factor ω should be set to a value between 0 and 2. The SOR method is equal to the Gauss–Seidel method for ω = 1, but may converge faster for other choices of the overrelaxation factor.

5 NUMERICAL SIMULATION

The surface harmonics method and the spheroidal harmonics method are here tested and compared numerically. The harmonic function used in the tests is the disturbing potential T computed from the EGM2008 global gravity field model (Pavlis et al. 2012) to M = 2190, where the Geodetic Reference System parameters adopted in the creation of EGM2008 (Pavlis et al. 2012, appendix A) were used to define the reference potential. The tests were performed using a modified version of the harmonic_synth software (Holmes & Pavlis 2006) for spherical harmonic synthesis, the SHTools software version 3.1 (Wieczorek 2015) for spherical harmonic analysis and in-house software for the coefficient transformations.

Because the surface harmonics method requires a grid of function values in terms of geocentric latitude θ and the spheroidal harmonics method requires a grid in terms of reduced latitude β, two regular global grids of the disturbing potential were generated through spherical harmonic synthesis on the EGM2008 reference ellipsoid with a resolution of 2.5 arcmin (see Fig. 1). These grids contain the data from which a solid spherical harmonic expansion is to be computed by the spheroidal and surface harmonics methods.

Figure 1.

Overview of numerical computations and comparisons (SHA: spherical/spheroidal harmonic analysis; SHS: spherical/spheroidal harmonic synthesis and CT: coefficient transformation).

The first step in this process is a spherical harmonic analysis on both grids. Through the Driscoll & Healy (1994) algorithm, as implemented in SHTools, harmonic series up to d/o 2159 were obtained. Note that the Driscoll & Healy (1994) algorithm requires a grid resolution of 2.5 arcmin to obtain a series to d/o 2159, but algorithms that require a lower resolution also exist (e.g. McEwen & Wiaux 2011). Analysis of the grid in terms of geocentric latitude results in a surface spherical harmonic series, while the grid in terms of reduced latitude results in a solid spheroidal harmonic series. Fig. 2 shows the spectra of these two series, indicating that they have similar power in all wavelengths. However, the coefficients are significantly different for high d/o, as is shown by the degree variances of the coefficient differences.

Figure 2.

Degree variances of (1) surface spherical HCs of the disturbing potential |$\overline{T}_{nm}^{\mathrm{e}}$| (black), (2) solid spheroidal HCs of the disturbing potential |$\overline{T}_{nm}^{u}$| (blue) and (3) the differences between those two (red), indicating that while both series have very similar power (the latter obscures the former), the differences between coefficients are significant.

5.1 Coefficient transformation

To test the spheroidal harmonics method, the Hotine–Jekeli transformation (eq. 11) was performed on the solid spheroidal harmonic series using the recurrence relations for Jekeli's renormalized Legendre functions of the second kind by Sebera et al. (2012, eq. 20). The summation was truncated at i = 40, resulting in a set of solid spherical HCs to d/o 2239. Due to the convergence of the series in eq. (11), the truncation at i = 40 is more than sufficient. For comparison, EGM2008 which is also computed using the Hotine–Jekeli transformation is truncated at d/o 2190.

To test the surface harmonics method, the reverse transformation detailed in Section 4.2 was performed on the surface spherical harmonic series. As in the Hotine–Jekeli transformation, all infinite series were truncated after 40 terms, that is, the summations in eqs (17) and (25) were performed up to j = 40, the summation in eq. (24) from i = −40 to 40 and the summation in eq. (27) was performed from j = −40 to 40. This is also more than sufficient (see Section 5.2).

The success of the reverse transformation hinges on the question whether eq. (28) holds up to M = 2239. Fig. 3 shows that this is clearly the case, as the dashed red line remains well below the solid red line for all degrees n. Only the case m = 0 is shown in this figure, as this provides a 'worst case'.

Figure 3.

Values of |λ _n00| (solid blue) and |${\sum _{i=-40, i\ne 0}^{40} |\lambda _{n0i}|}$| (dashed blue) from (17), as well as |μ _n00| (solid black) and |$ {\sum _{i=-40, i\ne 0}^{40} |\mu _{n0i}|}$| (dashed black) from (25) and |ν _n00| (solid red) and |${\sum _{i=-40, i\ne 0}^{40} |\nu _{n0i}|}$| (dashed red) from (27) based on GRS80 parameters.

The Gauss–Seidel method (eq. 30) was used, because it is difficult to predict an optimum overrelaxation factor for the SOR method. The convergence is rapid for all degrees, as can be seen in Fig. 4. A total of 10 iterations were applied to compute the final coefficients. The computational effort required for each iteration is insignificant compared to the effort required to compute the weights ν _nmi .

Figure 4.

Degree variances of approximate coefficients of the disturbing potential (24) (black) and of additions to those coefficients from iterations 1 to 6 (red through to blue) using Gauss–Seidel iteration (30), indicating rapid convergence for all degrees.

5.2 Comparison of solid spherical harmonic series

The spheroidal harmonics method and the surface harmonics method each provide a solid spherical harmonic expansion of the disturbing potential. These could easily be transformed into an expansion of the gravitational potential, like EGM2008, using the parameters of the reference ellipsoid. From a theoretical point of view, it can be expected that these two expansions are not identical, and neither will be identical to EGM2008. The reason for this is that each of the models is based on a different truncated series, with as a result a different omission error.

A comparison between the two solid spherical harmonic series resulting from the surface harmonics method and the spheroidal harmonics method is performed in the spectral domain and the space domain. First, Fig. 5 shows the differences between the power spectra of the two series. It can be seen that the differences between the coefficients are very small except for the very highest d/o close to M. This is as expected, because the coefficient transformation formulae (for both the surface harmonics method and the spheroidal harmonics method) show that the omission error in the series truncated at d/o 2159 will only affect solid spherical HCs close to that degree (2079 < n < 2239 based on i = 40, but any effect outside 2130 < n < 2190 is insignificant).

Figure 5.

Degree variances of (1) surface spherical HCs of the disturbing potential |$\overline{T}_{nm}^{\mathrm{e}}$| (black), (2) solid spheroidal HCs of the disturbing potential |$\overline{T}_{nm}^{u}$| (blue) and (3) the differences between those two (red), indicating that both series are nearly identical (the latter obscures the former) except for the coefficients close to the maximum degree.

Fig. 6 shows a close up of the range from degree 2120 to 2200. This shows that the coefficients resulting from the surface harmonics method decrease in value more rapidly than those resulting from the spheroidal harmonics method. Therefore, while the latter can safely be truncated after d/o 2190, like EGM2008, the former can already be truncated after d/o 2180. This also shows that in practice it is sufficient to evaluate significantly less than 40 terms. The large number of terms were only used in this study to highlight the high precision of the methods (see Section 5.3). Fig. 6 also helps to explain the differences between the highest degree coefficients of the topographic potential model dV_ELL_RET2012 and EGM2008, which was noted in Claessens & Hirt (2013) but not explained.

Figure 6.

Zoom-in on the 'tail' end of Fig. 5: degree variances of (1) surface spherical HCs of the disturbing potential |$\overline{T}_{nm}^{\mathrm{e}}$| (black), (2) solid spheroidal HCs of the disturbing potential |$\overline{T}_{nm}^{u}$| (blue) and (3) the differences between those two (red), highlighting the differences close to the maximum degree.

The two solid spherical harmonic series were also compared in the space domain. Fig. 7 shows the differences between regular global 2.5 arcmin grids of height anomalies in terms of geodetic latitude. The maximum difference is 7.6 mm and the rms of differences is 0.47 mm. Fig. 8 shows the same, but for gravity disturbances. The maximum difference is 2.5 mGal and the rms of differences is 0.15 mGal.

Figure 7.

Differences between height anomalies computed from solid spherical HCs that were generated with (1) the surface harmonics method and (2) the spheroidal harmonics method (units in metres; Robinson projection).

Figure 8.

Differences between gravity disturbances computed from solid spherical HCs that were generated with (1) the surface harmonics method and (2) the spheroidal harmonics method (units in milligal; Robinson projection).

A very similar pattern can be observed in Figs 7 and 8. In both figures, the differences are of a short-wavelength nature (in agreement with Figs 5 and 6), and are largest in mountainous regions near the equator. An explanation for these short-wavelength differences is that the two series to M = 2159, one of surface spherical HCs and the other of solid spheroidal HCs, do not contain the exact same signal. This can, for example, be seen from the fact that a transformation of surface spherical HCs up to d/o M to solid spheroidal HCs using the transformation formulae described here would result in (incomplete) non-zero solid spheroidal HCs beyond the maximum degree M (and the same argument holds vice versa). Hence, the two series of solid spherical HCs to M = 2239 do not contain the exact same signal either. Figs 7 and 8 give an indication of the differences caused by this.

5.3 Validation of coefficient transformations

The numerical precision of the surface and spheroidal harmonics methods is analysed by a comparison of height anomalies in the space domain. For validation of the spheroidal harmonics method, a regular global 2.5 arcmin grid of the disturbing potential in terms of reduced latitude was synthesized from the solid spheroidal harmonic series and subsequently transformed into height anomalies using Bruns's equation. This is compared to height anomalies on the same grid synthesized from the solid spherical harmonic series (see Fig. 1). Any discrepancies between these two grids must be the result of rounding errors and of the truncation of the transformation series. However, note that these grids are not identical to grids of height anomalies that would be synthesized directly from EGM2008 to d/o 2190.

The same process is applied to validate the surface harmonics method, using the surface harmonic series instead and synthesising on a grid in terms of geocentric latitude (see Fig. 1). Table 1 shows that both the spheroidal harmonics method and the surface harmonic methods achieve a very high precision with differences in height anomalies no larger than a few micrometres. The rms of differences is smallest for the surface harmonics method, which may be explained by the more rapid decrease of coefficients shown in Fig. 6.

Table 1.

Global statistics of transformation validation: differences between height anomalies synthesized from (1) surface spherical HCs and approximate solid spherical HCs (24), (2) surface spherical HCs and solid spherical HCs from the surface harmonics method (30) and (3) solid spheroidal HCs and solid spherical HCs from the spheroidal harmonics method (11) ((±1.0 × 10⁻¹²) indicates a value in between −10⁻¹² and +10⁻¹², that is, a value smaller than can accurately be determined given the number of digits used in the computations) (units in metres).

	Approximate	Surface CT	Spheroidal CT
Minimum	− 1.7 × 10−02	− 1.9 × 10−07	(±1.0 × 10−12)
Maximum	1.1 × 10−02	1.6 × 10−06	6.2 × 10−06
Mean	1.8 × 10−05	(±1.0 × 10−12)	8.0 × 10−11
rms	0.7 × 10−03	2.7 × 10−10	1.3 × 10−08

	Approximate	Surface CT	Spheroidal CT
Minimum	− 1.7 × 10−02	− 1.9 × 10−07	(±1.0 × 10−12)
Maximum	1.1 × 10−02	1.6 × 10−06	6.2 × 10−06
Mean	1.8 × 10−05	(±1.0 × 10−12)	8.0 × 10−11
rms	0.7 × 10−03	2.7 × 10−10	1.3 × 10−08

Table 1.

Global statistics of transformation validation: differences between height anomalies synthesized from (1) surface spherical HCs and approximate solid spherical HCs (24), (2) surface spherical HCs and solid spherical HCs from the surface harmonics method (30) and (3) solid spheroidal HCs and solid spherical HCs from the spheroidal harmonics method (11) ((±1.0 × 10⁻¹²) indicates a value in between −10⁻¹² and +10⁻¹², that is, a value smaller than can accurately be determined given the number of digits used in the computations) (units in metres).

	Approximate	Surface CT	Spheroidal CT
Minimum	− 1.7 × 10−02	− 1.9 × 10−07	(±1.0 × 10−12)
Maximum	1.1 × 10−02	1.6 × 10−06	6.2 × 10−06
Mean	1.8 × 10−05	(±1.0 × 10−12)	8.0 × 10−11
rms	0.7 × 10−03	2.7 × 10−10	1.3 × 10−08

	Approximate	Surface CT	Spheroidal CT
Minimum	− 1.7 × 10−02	− 1.9 × 10−07	(±1.0 × 10−12)
Maximum	1.1 × 10−02	1.6 × 10−06	6.2 × 10−06
Mean	1.8 × 10−05	(±1.0 × 10−12)	8.0 × 10−11
rms	0.7 × 10−03	2.7 × 10−10	1.3 × 10−08

Table 1 also shows the precision of the approximate solid spherical HCs that are computed as an intermediate step in the surface harmonics method (eq. 24). These approximate coefficients result in errors up to a few centimetres in height anomalies, and more significant errors for quantities that have more power in the higher degrees. As such, direct use of the approximate coefficients is not recommended.

In terms of numerical efficiency, the spheroidal harmonics method is more efficient than the surface harmonics method, mostly because of the two-step procedure required in the surface harmonics method. However, neither method is very demanding. Both coefficient transformations are less time-consuming than the generation of a global full-resolution grid of function values through spherical harmonic synthesis. If required, higher efficiency could be achieved in both methods—in particular in the surface harmonics method—by a lower truncation of convergent series.

6 SOLUTIONS FOR OTHER BOUNDARY CONDITIONS

Up to here it was assumed that values of the function of interest are given on the surface of the spheroid, which results in a Dirichlet-type boundary-value problem (e.g. Mackie 1965; Sigl 1985). However, a solid spherical harmonic expansion may also be sought if (a linear combination of) first- and second-order derivatives of the function are given on the spheroid.

If the derivatives are in the radial direction, both the spheroidal harmonics method and the surface harmonics method can easily be adapted to this situation, since there is a one-to-one relation between spherical HCs of a function and the spherical HCs of its radial derivatives. However, when derivatives in the normal direction to the spheroid are given, the spheroidal harmonics method cannot be used. The surface harmonics method, on the other hand, can be used in this situation. This is of significance, because common quantities in geodesy, such as gravity disturbances and gravity anomalies can be accurately related to normal derivatives of the disturbing potential.

In the generation of EGM2008 and earlier global gravity models, approximate ellipsoidal corrections were applied to gravity anomalies to obtain an expression that uses radial derivatives, following Gleason (1988) and Rapp & Pavlis (1990). These ellipsoidal corrections are unnecessary when using the surface harmonic method, as follows.

Let |${\overline{f}^{\prime }}_{nm}^{\mathrm{e}}$| be the surface spherical HCs of the first-order normal derivative f ′ with respect to the spheroidal surface

\begin{equation} f^{\prime }(r_{\mathrm{e}},\theta ,\lambda ) = \sum _{n=0}^{\infty } \sum _{m=-n}^n {\overline{f}^{\prime }}_{nm}^{\mathrm{e}} \overline{Y}_{nm} (\theta , \lambda ). \end{equation}

(32)

Likewise, let |${\overline{f}^{\prime \prime }}_{nm}^{\mathrm{e}}$| be the surface spherical HCs of the second-order normal derivative f ″ with respect to the spheroidal surface

\begin{equation} f^{\prime \prime }(r_{\mathrm{e}},\theta ,\lambda ) = \sum _{n=0}^{\infty } \sum _{m=-n}^n {\overline{f}^{\prime \prime }}_{nm}^{\mathrm{e}} \overline{Y}_{nm} (\theta , \lambda ). \end{equation}

(33)

These series can be used to solve Neumann and second-order boundary-value problems, respectively. Claessens & Featherstone (2008) give formulae for the computation of |${\overline{f}^{\prime }}_{nm}^{\mathrm{e}}$| and |${\overline{f}^{\prime \prime }}_{nm}^{\mathrm{e}}$| from a set of solid spherical HCs |$\overline{f}_{nm}^R$|⁠, assuming the function f is harmonic anywhere on and outside the spheroidal surface. These take the same form as (eq. 16), a weighted summation over coefficients of equal order m

\begin{equation} {\overline{f}^{\prime }}_{nm}^{\mathrm{e}} = \sum _{i=-\infty }^{\infty } \lambda ^{\prime }_{nmi} \overline{f}_{n-2i,m}^{R} \end{equation}

(34)

\begin{equation} {\overline{f}^{\prime \prime }}_{nm}^{\mathrm{e}} = \sum _{i=-\infty }^{\infty } \lambda ^{\prime \prime }_{nmi} \overline{f}_{n-2i,m}^{R}. \end{equation}

(35)

The equations for the weights |$\lambda ^{\prime }_{nmi}$| and |$\lambda ^{\prime \prime }_{nmi}$| are given in Claessens and Featherstone (2008, eqs 41 and 42).

The transformation matrices Λ′ and Λ ″, containing the weights |$\lambda ^{\prime }_{nmi}$| and |$\lambda ^{\prime \prime }_{nmi}$| respectively, are both diagonally dominant up to d/o 520, just as Λ (Claessens & Featherstone 2008). The reverse transformation can, however, be performed efficiently using a two-step procedure similar to the one described in Section 4.2.

Approximating normal derivatives by radial derivatives gives the approximate relations (cf. eq. 21)

\begin{equation} \overline{f}_{nm}^R \overline{Y}_{nm} \approx - \bigg ( \frac{r_e}{R} \bigg )^{n+2} \frac{R}{n+1} {\overline{f}^{\prime }}_{nm}^{\mathrm{e}} \overline{Y}_{nm} \end{equation}

(36)

\begin{equation} \overline{f}_{nm}^R \overline{Y}_{nm} \approx \bigg ( \frac{r_e}{R} \bigg )^{n+3} \frac{R^2}{(n+1)(n+2)} {\overline{f}^{\prime \prime }}_{nm}^{\mathrm{e}} \overline{Y}_{nm}. \end{equation}

(37)

This leads to expressions for the approximate solid spherical HCs (cf. eq. 24)

\begin{equation} \hat{f}_{nm}^R = \sum _{i=-\infty }^{\infty } \mu ^{\prime }_{nmi} {\overline{f}^{\prime }}_{n-2i,m}^{\mathrm{e}} \end{equation}

(38)

\begin{equation} \hat{f}_{nm}^R = \sum _{i=-\infty }^{\infty } \mu ^{\prime \prime }_{nmi} {\overline{f}^{\prime \prime }}_{n-2i,m}^{\mathrm{e}} \end{equation}

(39)

where the weight factors |$\mu ^{\prime }_{nmi}$| and |$\mu ^{\prime \prime }_{nmi}$| are computed in a very similar fashion as the weight factors μ _nmi in eq. (25)

\begin{equation} \mu ^{\prime }_{nmi} = - \frac{R}{n+1-2i} \sum _{j=|i|}^{\infty } \omega _{n+1-2i,j} K_{n-2i,m}^{2i,2j} \end{equation}

(40)

\begin{equation} \mu ^{\prime \prime }_{nmi} = \frac{R^2}{(n+1-2i)(n+2-2i)} \sum _{j=|i|}^{\infty } \omega _{n+2-2i,j} K_{n-2i,m}^{2i,2j}. \end{equation}

(41)

Eqs (40) and (41) contain singularities for n − 2i = −1 and −2, but it can be seen from eqs (38) and (39) that |$\mu ^{\prime }_{nmi}$| and |$\mu ^{\prime \prime }_{nmi}$| are only required for n − 2i ≥ 0, since spherical HCs of negative degree are discarded in the practical evaluation of eqs (38) and (39).

Once approximate solid spherical HCs are available, the second step of the reverse transformation is identical to that described in Section 4.2. The relation between the approximate coefficients |$\hat{f}_{nm}^R$| and the 'true' coefficients |$\overline{f}_{nm}^R$| is given by eqs (26) and (27), but with primes added to the weights λ _nmi , μ _nmi and ν _nmi . The reverse transformation can then again be solved iteratively using the Jacobi, Gauss–Seidel or SOR method (eqs 29–31).

As in the Dirichlet problem, the weights |$\nu ^{\prime }_{nmi}$| and |$\nu ^{\prime \prime }_{nmi}$| for the Neumann and second-order problems are strongly diagonally dominant up to high d/o (at least up to n = 2190). This diagonally dominant behaviour of the weights in all cases is similar, because the magnitude of the weight functions λ _nmi , |$\lambda ^{\prime }_{nmi}$| and |$\lambda ^{\prime \prime }_{nmi}$| is dominated by the contribution from the terms |$(R / r_{\mathrm{e}})^{n+1}$|⁠, |$(R / r_{\mathrm{e}})^{n+2}$| and |$(R / r_{\mathrm{e}})^{n+3}$|⁠, respectively, for large values of n. Since the inverse of these power terms are incorporated in μ _nmi , |$\mu ^{\prime }_{nmi}$| and |$\mu ^{\prime \prime }_{nmi}$|⁠, the cause of the non-diagonal dominance in the high degrees is 'counter-balanced'.

Boundary-value problems where the boundary condition is a linear combination of zero-, first- and second-order derivatives of a harmonic function, such as the Robin boundary-value problem, can also be solved using the surface harmonic approach presented here, by simply applying a linear combination of the coefficient transformations. A special example is the computation of a global gravity model from gravity anomalies on the surface of a spheroid. The forward transformation for that problem is given in Claessens & Hirt (2015), and the two-step procedure outlined here can be used for the reverse transformation.

7 CONCLUDING REMARKS

In this paper, it is shown how various boundary-value problems where the boundary is a spheroid can be solved using a transformation between solid and surface spherical HCs. It is shown that surface HCs of the harmonic function or its first- or second-order derivative, defined on the surface of any spheroid, can be expressed as a weighted summation over solid spherical HCs of that function. This provides a generalization of the well-known 'spherical' case, where surface HCs are defined on the surface of a sphere, as an alternative to existing methods. This is of particular interest to Earth sciences, since the surface of the Earth can accurately be described by a spheroid.

The method was compared numerically to the Hotine–Jekeli transformation between solid spheroidal and spherical HCs. These test show that both methods can achieve submicrometre precision in terms of height anomalies for a model to d/o 2239. However, both methods result in spherical harmonic models that are different by up to 7.6 mm in height anomalies and 2.5 mGal in gravity disturbances due to the different coordinate system used. The main advantage of the Hotine–Jekeli transformation is that it is numerically more efficient. The new method for transformation from surface to solid spherical HCs has as its main advantage that it can more easily be customized to different boundary quantities, such as gravity disturbances or gravity anomalies for global gravity modelling.

The software packages harmonic_synth (Holmes & Pavlis 2006) and SHTools (Wieczorek 2015) were used in the numerical study and Generic Mapping Tools (Wessel et al. 2013) was used for the creation of Figs 2–8. The authors are hereby acknowledged for making their software freely available.

REFERENCES

Comparison among three harmonic analysis techniques on the sphere and the ellipsoid

J. Appl. Geod.

2014

8

1

19

Handbook of Mathematical Functions

1972

Dover Publications

Spherical harmonic modelling to ultra-high degree of Bouguer and isostatic anomalies

J. Geod.

2012

86

499

520

On a relation between spherical and spheroidal harmonics

J. Phys. A: Math. Gen.

1977

10

1833

1836

Derivation of geomagnetic model to n = 63

Geophys. J. Int.

1989

97

421

441

Global heat flow: spherical harmonic representation

EOS, Trans. Am. geophys. Un.

1980

61

383

The Enhanced Magnetic Model 2015–2020

2015

National Centers for Environmental Information

NOAA

doi:10.7289/V56971HV, last accessed 10 December 2015

New relations among associated Legendre functions and spherical harmonics

J. Geod.

2005

79

398

406

Solutions to ellipsoidal boundary value problems for gravity field modelling

PhD thesis

2006

Perth

Curtin University of Technology

The Meissl scheme for the geodetic ellipsoid

J. Geod.

2008

82

513

522

Ellipsoidal topographic potential: new solutions for spectral forward gravity modeling of topography with respect to a reference ellipsoid

J. geophys. Res.

2013

118

5991

6002

A surface spherical harmonic expansion of gravity anomalies on the ellipsoid

J. Geod.

2015

89

10

1035

1048

Ellipsoidal corrections to potential coefficients obtained from gravity anomaly data on the ellipsoid, Rep. 371

Department of Geodetic Science and Surveying

1986

Columbus

Ohio State University

Transformation of spherical harmonic coefficients to ellipsoidal harmonic coefficients

Astronom. Astrophys.

2002

387

1114

1122

Computing Fourier transforms and convolutions on the 2-Sphere

Adv. Appl. Math.

1994

15

202

250

Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers

J. Geod.

2012

86

271

285

Recursive computation of oblate spheroidal harmonics of the second kind and their first-, second-, and third-order derivatives

J. Geod.

2013

87

303

309

Comparing ellipsoidal corrections to the transformation between the geopotential's spherical and ellipsoidal spectrums

Manuscr. Geod.

1988

13

114

129

Matrix Computations

1996

John Hopkins University Press

High resolution spherical and ellipsoidal harmonic expansions by fast Fourier transform

Stud. Geophys. Geod.

2014

58

595

608

Spherical harmonic analysis of earth's conductive heat flow

Int. J. Earth Sci.

2008

97

205

226

The Theory of Spherical and Ellipsoidal Harmonics

1931

Cambridge Univ. Press

Earth's heat flux revised and linked to chemistry

Tectonophysics

2005

395

159

177

A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalized associated Legendre functions

J. Geod.

2002

76

279

299

A Fortran program for very-high-degree harmonic synthesis (version 05/01/2006)

2006

The exact transformation between ellipsoidal and spherical harmonic expansions

Manuscr. Geod.

1988

13

106

113

On the computation and approximation of ultra-high-degree spherical harmonic series

J. Geod.

2007

81

9

603

615

A new isostatic model of the lithosphere and gravity field

J. Geod.

2004

78

368

385

GRIMM: the GFZ Reference Internal Magnetic Model based on vector satellite and observatory data

Geophys. J. Int.

2008

173

382

394

Boundary Value Problems

1965

Oliver and Boyd Ltd

On the efficient calculation of ordinary and generalized spherical harmonics

Geophys. J. Int.

1998

135

307

309

An ellipsoidal harmonic representation of Earth's lithospheric magnetic field to degree and order 720

Geochem. Geophys. Geosys.

2010

11

Q06015

doi:10.1029/2010GC003026

et al.

The 10th generation international geomagnetic reference field

Phys. Earth planet Inter.

2005

151

320

322

A novel sampling theorem on the sphere

IEEE Trans. Signal Process.

2011

59

5876

5887

Models of isostatic and dynamic topography, geoid anomalies, and their uncertainties

J. geophys. Res.

2000

105

28 199

28 211

The development and evaluation of the Earth Gravitational Model 2008 (EGM2008)

J. geophys. Res.

2012

117

B04406

doi:10.1029/2011JB008916

Heat flow from the Earth's interior: analysis of the global data set

Rev. Geophys.

1993

31

267

280

Coseismic deformation from earthquake faulting on a layered spherical earth

Geophys. J. Int.

1996

125

1

14

The development and analysis of geopotential coefficient models to spherical harmonic degree 360

J. Geophys. Res.

1990

95

B13

21 885

21 911

Approximation method for high-degree harmonics in normal mode modelling

Geophys. J. Int.

2002

151

309

313

Comparisons of global topographic-isostatic models to the Earth's observed gravity field, Rep. 388

Department of Geodetic Science and Surveying

1988

Columbus

Ohio State University

Fast spherical collocation: theory and examples

J. Geod.

2003

77

101

112

On computing ellipsoidal harmonics using Jekeli's renormalization

J. Geod.

2012

86

713

726

Introduction to Potential Theory

1985

Abacus Press

A new integral approach to geopotential coefficient determinations from terrestrial gravity or satellite altimetry data

Bull. Geod.

1988

62

93

101

Numerical problems in the computation of ellipsoidal harmonics

J. Geod.

1995

70

117

126

Linear Algebra and its Applications

1986

3rd edn

Harcourt Brace and Company

Degree 12 model of shear velocity heterogeneity in the mantle

J. geophys. Res.

1994

99

6945

6980

Generic mapping tools: improved version released

EOS, Trans. Am. geophys. Un.

2013

94

409

410

SHTOOLS, tools for working with spherical harmonics, version 3.1

2015

© The Author 2016. Published by Oxford University Press on behalf of The Royal Astronomical Society.

© The Author 2016. Published by Oxford University Press on behalf of The Royal Astronomical Society.

floresharse1939.blogspot.com

Source: https://academic.oup.com/gji/article/206/1/142/2606501

Share this post