/
accuracy.html
80 lines (61 loc) · 7.27 KB
/
accuracy.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN"
"http://www.w3.org/TR/html4/strict.dtd">
<html>
<head>
<title>SHTOOLS - Accuracy</title>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<link rel="stylesheet" type="text/css" href="css/sh.css">
<link rel="shortcut icon" type="image/vnd.microsoft.icon" href="images/favicon.ico">
<link rel="icon" type="image/vnd.microsoft.icon" href="images/favicon.ico">
<meta name="description" content="The spherical harmonic transforms in SHTOOLS are proven to be accurate to approximately degree 2800, corresponding to a spatial resolution of better than 4 arc minutes.">
</head>
<body>
<div class="main">
<p class="centeredimage"><img src="images/logo.jpg" width=894 height=135 alt="SHTOOLS --- Tools for working with spherical harmonics"></p>
<table class="menu">
<tbody>
<tr>
<td><a href="http://shtools.ipgp.fr/">HOME</a></td>
<td><a href="https://github.com/SHTOOLS/SHTOOLS/releases">DOWNLOAD</a></td>
<td class="selected"><a href="documentation.html">DOCUMENTATION</a></td>
<td><a href="faq.html">FAQ</a> </td>
</tr>
</tbody>
</table>
<p class="dir">
> <a href="../index.html" class="dir">Home</a> > <a href="documentation.html" class="dir">Documentation</a> </p>
<h1>Accuracy</h1>
<p>Three different types of algorithms are provided for calculating the spherical harmonic coefficients of a function. <tt>SHExpandLSQ</tt> performs a least squares inversion given data at arbitrary points on the sphere. This algorithm is slow for large spherical harmonic degrees, but may be desirable in the case where the data are noisy. <tt>SHExpandDH</tt> and <tt>SHExpandDHC</tt> perform an expansion using the sampling theorem of <i>Driscoll and Healy</i> (1994, Adv. Applied. Math., 15, 202-250) for real and complex data, respectively. For a bandlimited function that is sampled on <tt>N</tt> equally spaced points in both longitude and latitude (<tt>N</tt> being even), this transform is exact if the spherical harmonic bandwidth is less than or equal to <tt>N/2-1</tt>. <tt>SHExpandGLQ</tt> and <tt>SHExpandGLQC</tt> integrate numerically either a real or complex function over each <tt>P<sub>lm</sub></tt> using Gauss-Legendre quadrature, and it is this algorithm that will be discussed in detail below.</p>
<p class="nomarginbot">
The accuracy of the <tt>SHExpandGLQ</tt> spherical harmonic transform depends upon a number of factors, including:</p>
<ul>
<li> Accuracy of the Legendre functions</li>
<li> Accuracy of the Gauss-Legendre quadrature</li>
<li> Nature of the data set</li>
<li> Accumulation of roundoff errors</li>
</ul>
<p>
When using double precision floating point numbers, the standard Legendre functions overflow near spherical harmonic degree 15. Hence, all high-order expansions should make use of the geodesy, orthonormal, or Schmidt semi-normalized functions. It is well known that the standard three-term recursion formula that is utilized in calculating these leads to underflows near the poles for large values of <tt>m</tt>. To circumvent this problem, the normalized Legendre functions are here calculated using the method presented in <i>Holmes and Featherston</i> (2002, J. Geodesy, 76, 279-299), where the functions are first scaled by 10<sup>280</sup> <tt>sin</tt> θ prior to performing the recursions, and then unscaled at the end of the recursion. This ensures that the initial values of the recursion do not underflow, which would cause the <tt>P<sub>lm</sub></tt>s to be equal to 0 for all values of <tt>l</tt> for a given θ and <tt>m</tt>. The recursion used is to first calculate <tt>P<sub>mm</sub></tt>, and then to obtain <tt>P<sub>lm</sub></tt> for all values of <tt>l</tt>. The scaled portion of this routine does not lose precision and overflow until about degree 2800.</p>
<p>
Integration by Gauss-Legendre quadrature is exact when the integrand is a terminating polynomial. The zonal Legendre functions, <tt>P<sub>l0</sub></tt>, are polynomials of degree <tt>l</tt>, but the Legendre functions for <tt>m>0</tt> are not. Nevertheless, when <tt>m+m'</tt> is even, the product <tt>P<sub>lm</sub> P<sub>l'm'</sub></tt> is a polynomial of order <tt>l+l'</tt>. The present algorithm makes the assumption that this product is also a polynomial when <tt>m+m'</tt> is odd. Tests using 20% more Gauss-Legendre integration points than required by this assumption do not give rise to any improvement (utilizing the criterion presented below), indicating that this assumption is valid for all practical purposes.</p>
<p>
A set of synthetic spherical harmonic coefficients was created to test the accuracy of this routine. For a given spherical harmonic bandwidth, these coefficients were expanded in the space domain using <tt>MakeGridGLQ</tt>, and then re-expanded into spherical harmonics using <tt>SHExpandGLQ</tt>. The maximum and rms relative errors between the initial and final set of coefficients are plotted below in Figure 1 as a function of the bandwidth of the initial function. Two cases are shown. In the first, 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 l<sup>2</sup>. In the second, these coefficients were scaled such that the power spectrum was proportional to l<sup>-2</sup>. This second power spectrum is typical of geophysical observables, such as gravity and topography.</p>
<img class="imgcenter" src="images/rel_error.jpg" alt="rel_error.jpg" width=600 height=392 >
<p class="smaller">
<b>Figure 1.</b> Maximum and rms relative errors of the spherical harmonic coefficients as a function of spherical harmonic bandwidth, following a call to <tt>MakeGridGLQ</tt> and <tt>SHExpandGLQ</tt>.</p>
<p>
As is seen in the above plot, the errors associated with the transform and inverse pair are negligible to about degree 2600, and then grow somewhat between degrees 2700 and 2800. The errors increase in an quasi-exponential 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. RMS errors for the coefficients as a function of degree l are typically three orders of magnitude smaller than the maximum relative errors. While the errors are slightly larger for the set of coefficients that possessed a "red" power spectrum, the difference between the "red" and "blue" spectra is modest, implying that the form of the data do not strongly affect the accuracy of these routines. Relative errors for the routine <tt>SHExpandDH</tt> are nearly identical as those from <tt>SHExpandGLQ</tt>. The errors associated with the complex routines are lower by a few orders of magnitude.</p>
<p class="dir">
> <a href="../index.html" class="dir">Home</a> > <a href="documentation.html" class="dir">Documentation</a> </p>
<table class="footer2" summary = "SHTOOLS; Fortran and Python spherical harmonic transform software package">
<tbody>
<tr>
<td class="c1"><a href="http://www.ipgp.fr/">Institut de Physique du Globe de Paris</a></td>
<td class="c2"><a href="http://www.sorbonne-paris-cite.fr/index.php/en">University of Sorbonne Paris Cité</a></td>
<td class="c3">© 2016 <a href="https://github.com/SHTOOLS">SHTOOLS</a></td>
</tr>
</tbody>
</table>
</div>
</body>
</html>