Fast ultrahigh-degree global spherical harmonic synthesis on nonequispaced grid points at irregular surfaces

Abstract

In this paper, we presented a fast unified method to compute the gravity field functionals and their directional derivatives up to arbitrary orders on nonequispaced grid points at irregular surfaces using ultrahigh-degree models. The direct spherical harmonic synthesis (SHS) for computing the gravity field functionals at arbitrary locations through the Legendre functions is a time-consuming task for high-order and -degree models. Besides, to compute the derivatives of SHS in terms of latitude, the derivatives of the Legendre functions are needed. Therefore, we used Fourier coefficients of Wigner d-functions to compute the directional derivatives of SHS up to arbitrary orders. We also showed that these functions and their derivatives up to order 2 are stable up to ultrahigh-degree $2^{14} = 16{,}384$ using extended double precision (i.e., 80 bits variables). Although 2D-FFT can accelerate the computation of global SHS (GSHS), it restricts the results on equispaced grid points. Hence, we used the nonequispaced FFT (NFFT) for computing GSHS on irregular grid points on the sphere that it is the fast nonequispaced GSHS (NGSHS). For maximum degree N and computing points of ${\mathcal {O}}(N^2)$ with arbitrary locations, the direct computation methods have the complexity of ${\mathcal {O}}(N^4)$. But the presented algorithm with and without precomputed Fourier coefficients of Wigner d-functions has the complexity of ${\mathcal {O}}(N^2 \log ^2 N + N^2 s^2)$ and ${\mathcal {O}}(N^3 + N^2 s^2)$, respectively, where s is cutoff parameter of convolution in NFFT. Using a convolution technique in frequency domain, the NGSHS on the ellipsoid was computed. For computation the gravity field functionals by the NGSHS at irregular surfaces, we defined the Taylor expansion and the Padé approximation both on the sphere and on the ellipsoid. The results showed that the constructed Padé approximation on the ellipsoid provides better accuracy. Finally, we showed that the introduced unified algorithm achieves the required accuracy and that it is faster than direct computations.

Appendix

To prove Eq. (31), let the signal f be given on the ellipsoid by

$$\begin{aligned} f(r(\vartheta ), \vartheta , \lambda )= & {} \sum _{\ell =0}^{N} \left( \frac{R_{\mathrm {E}}}{r(\vartheta )}\right) ^{\ell +1} \nonumber \\&\times \sum _{m=-N}^{N} {\mathrm {e}}^{\jmath m \lambda } \sum _{m'=-N}^{N} {\mathrm {e}}^{\jmath m' \vartheta } \; {\bar{a}}_{m' m}^\ell \; \bar{f}_{\ell m}. \end{aligned}$$

(74)

Considering the Fourier expansion of $(R_{\mathrm {E}}/r(\vartheta ))^{\ell +1}$ as

$$\begin{aligned} \left( \frac{R_{\mathrm {E}}}{r(\vartheta )}\right) ^{\ell +1} = \sum _{m'=-N}^{N} {\mathrm {e}}^{\jmath m' \vartheta } \; R_{\ell m'}, \end{aligned}$$

(75)

the signal on the ellipsoid can be presented by

$$\begin{aligned} f(r(\vartheta ), \vartheta , \lambda ) = \sum _{\ell =0}^{N} r_{\ell }(\vartheta , \lambda ) \; f_{\ell }(\vartheta , \lambda ), \end{aligned}$$

(76)

where

$$\begin{aligned} f_{\ell }(\vartheta , \lambda ) = \sum _{m=-N}^{N} {\mathrm {e}}^{\jmath m \lambda } \sum _{m'=-N}^{N} {\mathrm {e}}^{\jmath m' \vartheta } \; {\bar{a}}_{m' m}^\ell \; \bar{f}_{\ell m} \end{aligned}$$

(77)

is the signal on the sphere $r=R_{\mathrm {E}}$ per degree and

$$\begin{aligned} r_{\ell }(\vartheta , \lambda )= & {} \left( \frac{R_{\mathrm {E}}}{r(\vartheta )}\right) ^{\ell +1}\nonumber \\= & {} \sum _{m=-N}^{N} {\mathrm {e}}^{\jmath m \lambda } \sum _{m'=-N}^{N} {\mathrm {e}}^{\jmath m' \vartheta } \; R_{\ell m m'}, \end{aligned}$$

(78)

in which

$$\begin{aligned} R_{\ell m m'} = {\left\{ \begin{array}{ll} R_{\ell m'} &{}\text {if}\ m = 0\\ 0 &{}\text {otherwise}. \end{array}\right. } \end{aligned}$$

(79)

Since $r_{\ell }(\vartheta , \lambda )$ is multiplied by $f_{\ell }(\vartheta , \lambda )$ in space domain in Eq. (76), considering Eqs. (77) and (78), one can convolve $R_{\ell m m'}$ with ${\bar{a}}_{m' m}^\ell \; \bar{f}_{\ell m}$ in frequency domain to achieve the same result by

$$\begin{aligned} f(r(\vartheta ), \vartheta , \lambda )= & {} \sum _{\ell =0}^{N} \sum _{m=-N}^{N} {\mathrm {e}}^{\jmath m \lambda } \sum _{m'=-N}^{N} {\mathrm {e}}^{\jmath m' \vartheta } \nonumber \\&\times R_{\ell m m'} \circledast {\bar{a}}_{m' m}^\ell \; \bar{f}_{\ell m}. \end{aligned}$$

(80)

Finally, rearranging the summations, the signal is given by

$$\begin{aligned} f(r(\vartheta ), \vartheta , \lambda )= & {} \sum _{m=-N}^{N} {\mathrm {e}}^{\jmath m \lambda } \sum _{m'=-N}^{N} {\mathrm {e}}^{\jmath m' \vartheta } \nonumber \\&\times \sum _{\ell =0}^{N} R_{\ell m m'} \circledast {\bar{a}}_{m' m}^\ell \; \bar{f}_{\ell m}. \end{aligned}$$

(81)