Numerical computations
Spherical harmonic reconstructions
If the spherical harmonic coefficients of a function were known up to a maximum degree , the function could be evaluated at arbitrary colatitude and longitude using the equation
\begin{equation} f\left(\theta,\phi\right) = \sum_{l=0}^{L} \sum_{m=l}^l f_{lm} \, Y_{lm}\left(\theta,\phi \right). \end{equation}
Using the separate variables and for the cosine and sine spherical harmonic coefficients, respectively, and after interchanging the order of summations over and , the function can be expressed equivalently as
Defining the two component vector
The function can be written more simply as
For a given latitude band, the function can thus be evaluated on a series of grid nodes all at the same time using an inverse fast Fourier transform. For this operation, SHTOOLS makes use of the highly optimized software package FFTW [Frigo and Johnson 2005] that supports grids of arbitrary length and both real and complex data. To adequately sample the function in longitude, a minimum of data points should be employed.
The coefficients and depend upon the associated Legendre functions, which are calculated efficiently using standard threeterm recursion relations over adjacent spherical harmonic degrees. For a given colatitude, the sectoral term is first calculated using an analytic equation, and then is calculated for all values of . To circumvent numerical underflows near the poles for large values of , the associated Legendre functions are calculated using the approach of Holmes and Featherstone [2002], where the sectoral terms are multiplied by prior to performing the recursions, and then appropriately unscaled at the end of the recursion. This ensures that the Legendre functions are accurate to about degree 2800.
Spherical harmonic transforms
The spherical harmonic coefficients of a function can be calculated from the relation
Defining the two intermediary variables and
it is seen that for a given colatitude and degree , all of the angular orders can be calculated at once by making use of a fast Fourier transform of the function . Replacing the integral over latitude with a numerical quadrature rule, the spherical harmonic coefficients can be calculated as
Here, is the latitudinal weight, and is the number of latitudinal points over which the integration is performed.
Supported grid formats
It is possible to choose the weights and the locations of the latitudinal sampling points such that the quadrature in equation \eqref{eq:quadrature} is exact. SHTOOLS makes use of two grid formats that accommodate exact quadrature using GaussLegendre quadrature [e.g., Press et al. 1992], and exact quadrature using regular grids sampled according to the Driscoll and Healy (1994) theorem. In both techniques, the quadrature is exact only when the function being integrated is a terminating polynomial. The functions and can be approximated as polynomials of maximum degree , and when multiplied by the associated Legendre function, the integrand is approximately a polynomial of maximum degree . The following table summarizes the properties of the different types of grids supported by SHTOOLS.
Acronym  Name  Shape ()  

DH  Driscoll and Healy  
DH2  Driscoll and Healy  
GLQ  GaussLegendre Quadrature 
GaussLegendre Quadrature
For the case of GaussLegendre quadrature (GLQ
), the quadrature is exact when the function is sampled in latitude at the zeros of the Legendre Polynomial of degree . Since the function also needs to be sampled on equally space grid nodes for the Fourier transforms in longitude, the function is sampled on a grid of size .
Driscoll and Healy [1994]
The second type of grid is for data that are sampled on regular grids. As shown by Driscoll and Healy [1994], an exact quadrature exists when the function is sampled at equally spaced nodes in latitude and equally spaced nodes in longitude. For this sampling (DH
), the grids make use of the longitude band at 90 N, but not 90 S, and the number of samples is , which is always even. Given that the sampling in latitude was imposed a priori, these grids contain almost twice as many samples in latitude as the grids used with GaussLegendre quadrature.
For geographic data, it is common to work with grids that are equally spaced in degrees latitude and longitude. SHTOOLS provides the option of using grids of size , and when performing the Fourier transforms for this case (DH2
), the coefficients and with are discarded.
Accuracy
To test the accuracy of the spherical harmonic transforms and reconstructions, two sets of synthetic spherical harmonic coefficients were created. Each coefficient was chosen to be a random Gaussian distributed number with unit variance, and the coefficients were then scaled such that the power spectrum was proportional to either or . For a given spherical harmonic bandwidth , the function was reconstructed in the space domain using the GaussLegendre quadrature implementation and then reexpanded into spherical harmonics. The maximum and rootmean square (rms) relative errors between the initial and final set of coefficients were then computed as a function of spherical harmonic degree, as plotted in the figure below.
The errors associated with the transform and inverse pair increase in an quasiexponential manner, with the maximum relative error being approximately 1 part in a billion for degrees close to 400, and about 1 part in a million for degrees close to 2600. Though the errors are negligible to about degree 2600, they then grow somewhat between degrees 2700 and 2800. The rms errors for the coefficients as a function of the bandwidth are typically three orders of magnitude smaller than the maximum relative errors.
The relative errors are nearly the same for the two tested power spectra, implying that the form of the data do not strongly affect the accuracy of the routines. Relative errors using the Driscoll and Healy [1994] sampled grids are nearly identical to those using GaussLegendre quadrature. The errors associated with the routines for complex data are lower by a few orders of magnitude.
Timing tests
The speed of the spherical harmonic reconstructions and transforms were tested for both real and complex data using the GaussLegendre and Driscoll and Healy [1994] quadrature implementations. The amount of time in seconds required to perform these operations is plotted in the figure below as a function of the spherical harmonic bandwidth of the function. These calculations were performed on a modern Mac Pro 2.7 GHz 12Core Intel Xeon E5 using 64 bit executables and level 3 optimizations.
For the real GaussLegendre quadrature routines, the transform time is on the order of one second for degrees close to 800 and about 30 seconds for degree 2600. For the real Driscoll and Healy [1994] sampled grids, the transform time is close to a second for degree 600 and about 1 minute for degrees close to 2600. The complex routines are slower by a factor of about 1.4.
References

Driscoll, J. R. and D. M. Healy, Computing Fourier transforms and convolutions on the 2sphere, Adv. Appl. Math., 15, 202250, doi:10.1006/aama.1994.1008, 1994.

Frigo, M., and S. G. Johnson, The design and implementation of FFTW3, Proc. IEEE, 93, 216–231, doi:10.1109/JPROC.2004.840301, 2005.

Holmes, S. A., and W. E. Featherstone, A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions, J. Geodesy, 76, 279 299, doi:10.1007/s0019000202162, 2002.

Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, “Numerical Recipes in FORTRAN: The Art of Scientific Computing,” 2nd ed., Cambridge Univ. Press, Cambridge, UK, 1992.