Journal ArticleDOI

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

01 May 2002-Journal of Geodesy (Springer-Verlag)-Vol. 76, Iss: 5, pp 279-299
TL;DR: In this article, Horner's method is used to compute a scaled version of the Legendre function, complete to degree and order 2700 for all latitudes (except at the poles for first derivatives).
Abstract: Spherical harmonic expansions form partial sums of fully normalised associated Legendre functions (ALFs). However, when evaluated increasingly close to the poles, the ultra-high degree and order (e.g. 2700) ALFs range over thousands of orders of magnitude. This causes existing recursion techniques for computing values of individual ALFs and their derivatives to fail. A common solution in geodesy is to evaluate these expansions using Clenshaw's method, which does not compute individual ALFs or their derivatives. Straightforward numerical principles govern the stability of this technique. Elementary algebra is employed to illustrate how these principles are implemented in Clenshaw's method. It is also demonstrated how existing recursion algorithms for computing ALFs and their first derivatives are easily modified to incorporate these same numerical principles. These modified recursions yield scaled ALFs and first derivatives, which can then be combined using Horner's scheme to compute partial sums, complete to degree and order 2700, for all latitudes (except at the poles for first derivatives). This exceeds any previously published result. Numerical tests suggest that this new approach is at least as precise and efficient as Clenshaw's method. However, the principal strength of the new techniques lies in their simplicity of formulation and implementation, since this quality should simplify the task of extending the approach to other uses, such as spherical harmonic analysis.

### 1 Introduction

• Current geodetic practice is witnessing an increase in the construction and use of ultrahigh degree spherical harmonic expansions of the geopotential or topography.
• These have taken the form of simpleeffects models (eg. Featherstone, 1999; Nov´ak et al., 2001) for which synthetic geopotential coefficients up to degree and order2700 and2160, respectively, were produced without reference to a mass distribution.
• Hybrids of source and effects models also exist.
• Haagmans (2000) combines empirically determined coefficients with synthetic ones derived from numerical integration over isostatically compensated source masses to degree and order2160.
• The numerical means for including the necessary ALFs constitutes the principal chal- lenge to evaluating ultra-high degree spherical harmonic expansions.

### 1.1 Spherical Harmonic Expansions

• The general notationS (d), (d)m , X (d) m andP (d) nm( ) is used whenever a textual or mathematical reference applies to both the undifferentiated quantities and the first derivatives simultaneously.
• Underflows in the computation of anyP (d) nm( ) excludes the corresponding (matching degree and order) coefficients from contributing toS(d), whereas an overflow in the computation of one or moreP (d) nm( ) prevents any result forS (d) from being achieved at all.
• Thus, IEEE double precision only permits a maximum range of 620 orders of magnitude within which to compute and store the requiredP (d) nm( ) values.
• A numerically more stable alternative to these standard recursion relations is the Clenshaw (1955) method (cf. Tscherning and Poder, 1982; Gleason, 1985; Deakin, 1998).
• These principles are quite simple, both in concept and in application.

### 2 Forward Recursions for the Calculation of ALFs

• The most direct approach for evaluatingS (d) (Eq. 1) employs a recursive algorithm to computeP nm( ).
• These values ofP (d) nm( ) are multiplied by the correspondingEnm terms to yield the required series values ofX (d) m (Eq. 3), which subsequently yield (d)m (Eq. 2) and henceS(d) (Eq. 1).
• UnlikeP nm( ) cosm sinm , P (d6=0) nm ( ) cosm sinm is not normalised in 9 the correct usage of the word since its average squared value integrated over the unit sphere is not unity.
• This paper focuses solely on the computation of fully normalised ALFs and their first derivatives.
• No numerical tests were conducted for the above quasinormalisations.

### 2.1 Standard Forward Column Methods

• The most popular recursive algorithm used for computingP nm( ) in geodesy can be obtained by fully normalising, for example, Magnus et al. (1966, Eq. 4.3.3(2)).
• Note that in Fig. 2, the degree increases in rows down, the order increases in columns to the right, and the diagonal elements of the matrix are the sectoral values.
• This nomenclature will be employed throughout the paper.

### 2.2 Standard Forward Row Methods

• The next approach is termed thestandard forward row recursion (Fig. 3), and appears to be rarely used in geodesy.
• As with the standard forward column recursion (Section 2.1), the sectoralPmm( ) serve as seed values for the forward row recursion, and can be computed using Eq. (13).
• The standard forward row recursion computes nonsectoralP nm( ) of constantn (a ‘row’ in Fig. 3) and sequentially decreasingm (to the left (ie., ‘forward’) from the diagonal in Fig. 3).
• Using the same argument to that introduced for the forward column recursion, the nonexistant value ofP n;n+1( ) required to computeP n;n 1( ) in Eq. (18) may be disregarded because the corresponding recursion coefficient,hn;n 1, is always zero.
• Note that, to computeP nm( ) using the forward row recursion, Eq. (18) uses the cor- responding sectoral values of the samen, rather than the samem, as seed values.

### 2.3 Numerical Problems with the Standard Forward Methods

• Even when applied in IEEE double precision, both the standard forward column (Eq. 11) and standard forward row (Eq. 18) recursions will underflow forM > 1900 in the colatitude range 20Æ 160Æ.
• The numerical instability of both these forward recursions is noted in the geodetic literature (eg. Gleason, 1985) and elsewhere (eg. Libbrecht, 1985).
• 0 (ie., towards the poles) and asm increases.
• Accordingly, the high degree and order values ofPmm( ) will exceed the range of magnitudes capable being stored in IEEE double precision, thereby resulting in an underflow.

### 2.4 Other Normalisations and the Edmonds Recursion

• Belikov (1991) and Belikov and Taybatorov (1991) present a suite of recursive algorithms for computing the quantitieŝP (d)nm, whereP̂ (d) nm are related to un-normalisedP (d) nm according to the modified normalisation P̂nm( ) = 2 m n! (n+m)! Pnm( ) (25) However, this approach is also subject to numerical limitations.
• The description of the test 14 results which support this claim indicate that these computations were only performed at the equator, although the point is not clear.
• Nevertheless, Risbo (1996) includes a Fortran 77 subroutine for this recursion.
• Therefore, although it may be possible to incorporate the new approaches presented here into the Risbo (1996) approach, it is clear that the Edmonds recursion alone does not solve the numerical problems encountered towards the poles as described in Section 1.2.

### 2.5 The Modified Forward Row Method

• A simple, yet effective, method by which this problem of underflowingP mm( ) may be avoided is to eliminate theum term from the recursion process in Eq. (13).
• This will be demonstrated numerically in Section 4. 16 2.6 Modified Forward Column Method Values ofP (d) nm( ) um may also be computed using what will be termed thefirst modified forward column recursion.
• Thus, for the remainder of this paper, no partial sumsS(1) will be computed at the poles.

### 3.2 Reverse Column Method

• Results from timing tests presented in Section 4.3 show the reverse column methods, described below, to be highly inefficient in comparison with the other approaches presented in this paper for evaluating the required partial sums.
• Thus, the reverse column methods are used here to highlight the basic similarities and differences between these two approaches.
• To compute any value ofPnm( ) Pmm( ) , the second modified forward column recursion (eg. 34) aggregates the necessary(almt) andblm recursive terms in the sequence of increasing degreel (sequentially down each column in Fig. 2).
• The effect of using the recursion in Eq. (46) in this way is to sequentially aggregate the (almt) andblm recursive terms, in the sequence of decreasingl, until the recursion terminates at the computation ofsmm = Pnm( ) Pmm( ) .

### 3.3 Standard Clenshaw Methods

• The standard Clenshaw methods, summarised below, closely resemble the reverse column recursions (Section 3.2).
• The Clenshaw (1955) approach, which was formulated originally to evaluate partial sums of Chebyshev polynomials, was adapted for use in geodesy by Gulick (1970) to compute partial sums ofP (d) nm( ).
• Section 3.3.1 introduces a sim- ple implementation of the Clenshaw (1955) approach, whilst Section 3.3.2 presents the implementation that is used more commonly in geodesy (cf. Gleason, 1985).

### 3.3.1 The first Clenshaw method

• The simplest implementation of the Clenshaw (1955) technique, herein termed thefirst Clenshaw method, uses the recursions in Eqs. (46) and (47) to compute, directly, the intermediate sumsX (d) m Pmm( ) , without evaluating individual values ofP (d) nm( ) Pmm( ) .
• As in the reverse column recursion (Section 3.2), the recursive process terminates at the computation ofsmm (on the diagonal in Fig. 7), except that, in this case, the sectoralsmm = Xm Pmm( ) .
• This study is confined to the computation ofS (1).
• These seed values allow the recursive computation of all_slm , of constantm and sequentially decreasingl, from _snm to _smm .

### 3.3.2 The second Clenshaw method

• This allows the recursive computation of all slm , of constantm (a column in Fig. 7), and sequentially increasingl (upwards and towards the diagonal in Fig. 7), fromsnm to smm , wheresmm = Xm Pmm( )qm .
• The standard approach (cf. Gleason, 1985; Deakin, 1998) is to combine these values using the implementation of Horner’s scheme given in Eq. (50).

### 4.1 Viable Methods

• These algorithms are summarised in Table 1.
• As mentioned in Section 2.6, the first (MFC-1) and second (MFC-2) modified forward column recursions are, essentially, a single method.
• This leaves, to this point, five separate methods for computing spherical harmonic expansions.
• The tests of numerical accuracy will use analytic solutions for the sum of the square ofP (d) nm( ), to compare the modified forward row and modified forward column algorithms only.

### 4.2 Relative Numerical Precision

• The first step in comparing different methods for computing the partial sumsS (d) in Eq. (1) is to choose some appropriate values forEnm and .
• The relative precision (RP) for each was calculated using RP = s(d) s(d)s(d) (55) wheres(d) is the value of Eq. (54) for the summation computed in precision by each respective method ands(d) is the result for the same sum computed using the second Clenshaw method in precision.
• That is, all combinations of the five algorithms with the three implementations of Horner’s scheme showed that the choice of implementation of Horner’s scheme was irrelevant to the observed precision of the algorithm.
• One interesting feature is that the relative precision signatures obtained from all of the column methods for the northern latitudes are almost identical, whilst the signatures for the southern latitudes are not (cf. Figs. 8 through 10).
• This phenomenon is not observed in the accuracy tests for the two modified forward algorithms (Section 4.4).

### 4.3 Numerical Efficiency

• The five methods that successfully computedS (d), for M = 2700 and for integer values of to the poles, were tested for their relative numerical efficiency.
• The square roots and inverted square roots required to construct the recursion coefficients were computed once by each algorithm and then stored for multiple use in the synthesis 32 subroutines.
• The CPU times for the reverse column algorithm are excessively large and so were extrapolated from the computation times for a single parallel.
• All computations were performed, once again, on aSun Ultra 10 (333MHz) workstation that uses a virtual or ‘swapped’ RAM configuration, which is slower than actual RAM.
• It should be noted that these CPU times, in addition to showing the relative efficiency of each approach, are also functions of the computer architecture, compiler and programming language employed, as well as the programmer’s implementation of these algorithms.

### 4.4 Accuracy

• As mentioned, the numerical evaluations presented in Section 4.2 are tests of precision only, since both the tested algorithms and the ‘control’ algorithm computed in IEEE extended double precision may contain shared systematic errors.
• These could be due to any one of compiler, computer architecture or programming errors, for example.
• This necessitates that the squaring ofP (d) nm( ) um 10 280, as well as their combination using Horner’s scheme, be performed in IEEE extended double precision.
• The modified forward column algorithm becomes increasingly less accurate than the modified forward row algorithm as the computation approaches the poles.
• These results support those obtained for these two algorithms in the precision tests (Section 4.2).

### 5 Summary, Conclusion and Recommendation

• This paper has shown that standard Clenshaw methods for evaluating high degree spherical harmonic expansions derive their stability from simple numerical principles.
• An efficient technique for the computation of the gravimetric quantities from geoptential earth models.
• 547–555. Libbrecht KG (1985) Practical considerations for the generation of large order spherical harmonics, also known as Geophys J Int 139.
• The modified forward column and modified forward row recursions are immediately more versatile than the standard Clenshaw methods, since they do not automatically combine quantities of the same orderm.

Did you find this useful? Give us your feedback

Content maybe subject to copyright    Report

A uniﬁed approach to the Clenshaw summation
and the recursive computation of very high
degree and order fully normalised associated
Legendre functions
S. A. Holmes
1
, W. E. Featherstone
2
1
Department of Spatial Sciences, Curtin University of Technology, GPO Box U1987, Perth, WA 6845, Australia.
e-mail: holmes@vesta.curtin.edu.au, Tel: +61 8 9487 3838 Fax: +61 8 9266 2703
2
Department of Spatial Sciences, Curtin University of Technology, GPO Box U1987, Perth, WA 6845, Australia.
e-mail: W.Featherstone@cc.curtin.edu.au, Tel: +61 8 9266 2734 or 0401 103 734 Fax: +61 8 9266 2703
Send offprint requests to: S. A. Holmes

2
Abstract. Spherical harmonic expansions form partial sums of fully normalised asso-
ciated Legendre functions (ALFs). However, when evaluated increasingly close to the
poles, the ultra-high degree and order (eg. 2700) ALFs range over thousands of orders
of magnitude. This causes existing recursion techniques for computing values of indi-
vidual ALFs and their derivatives to fail. A common solution in geodesy is to evaluate
these expansions using Clenshaw’s (1955) method, which does not compute individual
ALFs or their derivatives. Straightforward numerical principles govern the stability of
this technique. This paper employs elementary algebra to illustrate how these principles
are implemented in Clenshaw’s method. It also demonstrates how existing recursion al-
gorithms for computing ALFs and their ﬁrst derivatives are easily modiﬁed to incorporate
these same numerical principles. These modiﬁed recursions yield scaled ALFs and ﬁrst
derivatives, which can then be combined using Horner’s scheme to compute partial sums,
complete to degree and order 2700, for all latitudes (except at the poles for ﬁrst deriva-
tives). This exceeds any previously published result. Numerical tests suggest that this new
approach is at least as precise and efﬁcient as Clenshaw’s method. However, the principal
strength of the new techniques lies in their simplicity of formulation and implementation,
since this quality should simplify the task of extending the approach to other uses, such
as spherical harmonic analysis.
Key words. Spherical harmonic expansions, Fully normalised associated Legendre Func-
tions, Clenshaw summation, Recursion, Horner’s scheme

3
1 Introduction
Current geodetic practice is witnessing an increase in the construction and use of ultra-
high degree spherical harmonic expansions of the geopotential or topography. For ex-
ample, Wenzel (1998) released coefﬁcients up to degree
1800
, which were empirically
derived to describe the gravitational potential of the Earth. Wenzel (1998) states that the
maximum degree of
1800
for the spherical harmonic model was set by the numerical
stability of the recursion algorithm adopted to compute the required fully normalised as-
sociated Legendre functions (ALFs).
The recent interest in synthetic Earth gravity models, used for comparing and validat-
ing gravity ﬁeld determination techniques, has already seen the use of ultra-high degree
spherical harmonic expansions. These have taken the form of simple effects models (eg.
Featherstone, 1999; Nov´ak et al., 2001) for which synthetic geopotential coefﬁcients up to
degree and order
2700
and
2160
, respectively, were produced without reference to a mass
distribution. There is also interest in source models in which synthetic geopotential coef-
ﬁcients are generated by analytical or numerical Newtonian integration over a synthetic
global density distribution and topography (eg. Pail, 1999). Hybrids of source and effects
models also exist. For example, Haagmans (2000) combines empirically determined co-
efﬁcients with synthetic ones derived from numerical integration over isostatically com-
pensated source masses to degree and order
2160
. Lastly, other scientiﬁc disciplines, such
as meteorology, quantum physics and electronic engineering, are also also showing in-
creased interest in high degree spherical harmonic modelling and analysis.
The numerical means for including the necessary ALFs constitutes the principal chal-
lenge to evaluating ultra-high degree spherical harmonic expansions.Therefore, it is timely

4
to critically examine the accuracy and numerical efﬁciency of algorithms that compute in-
dividual ALFs and their partial sums.
1.1 Spherical Harmonic Expansions
Truncated spherical harmonic expansions of a function, or its derivatives, reduce to sums
S
(
d
)
of ALFs or their
d
-th derivatives, respectively. These are
S
(
d
)
=
c
M
X
m
=0
(
d
)
m
(1)
where
(
d
)
m
=
2
X
=1
8
>
>
<
>
>
:
X
(
d
)
m
cos
m
for
=1
X
(
d
)
m
sin
m
for
=2
(2)
and
X
(
d
)
m
=
M
X
n
=
E
nm
P
(
d
)
nm
(
)
(3)
For arguments of spherical polar coordinates (
r
,
,
) and for integer degree
n
0
and
order
0
m
n
:
M
is the maximum ﬁnite degree of the spherical harmonic expansion;
is an integer that may vary with
m
;
c
is a real numbered constant;
E
nm
is a real num-
ber incorporating the fully normalised spherical harmonic coefﬁcients,
C
nm
1
and
C
nm
2
;
P
nm
(
)
are the fully normalised ALFs; the superscript
(
d
)
indicates the
d
-th derivative
with respect to
, or deﬁnite integration (
d
=
1
) between two parallels. This paper deals
only with undifferentiated functions (
d
=0
) or ﬁrst derivatives of these functions (
d
=1
).
For
d
= 0
, the superscript
(
d
)
is omitted. Thus
S
(0)
,
(0)
m
,
X
(0)
m
and
P
(0)
nm
(
)
are written
S
,
m
,
X
m
and
P
nm
, respectively. First derivatives of these quantities are written
S
(1)
,
(1)
m
,
X
(1)
m
and
P
(1)
nm
(
)
, respectively. The general notation
S
(
d
)
,
(
d
)
m
,
X
(
d
)
m
and
P
(
d
)
nm
(
)
is used whenever a textual or mathematical reference applies to both the undifferentiated
quantities and the ﬁrst derivatives simultaneously.

5
The example of a truncated spherical harmonic expansion of the gravitational potential
V
(
r;;
)
is instructive here. Often, it is written as
V
(
r;;
)=
GM
r
+
GM
r
M
X
n
=2
a
r
n
n
X
m
=0
(
C
nm
1
cos
m
+
C
nm
2
sin
m
)
P
nm
(
)
(4)
where
GM
is the product of the Universal gravitational constant and the mass of the
Earth. Alternatively, Eq. (4) may be written as
V
(
r;;
) =
GM
r
+
GM
r
M
X
m
=0
"
cos
m
M
X
n
=
a
r
n
C
nm
1
P
nm
(
)
+sin
m
M
X
n
=
a
r
n
C
nm
2
P
nm
(
)
#
(5)
where
is either
2
or
m
; whichever is the greater. Relating Eq. (5) to the form of Eqs. (1)
to (3) yields
E
nm
=
8
>
>
<
>
>
:
a
r
n
C
nm
1
;
for
=1
a
r
n
C
nm
2
;
for
=2
(6)
and
X
m
=
M
X
n
=
a
r
n
C
nm
P
nm
(
)
(7)
When evaluating gravimetric quantities (eg., disturbing potential, geoid heights, grav-
ity anomalies, etc.) in a sequence of points for which
r
and
are constant (ie., along a
geodetic parallel), the form of Eq. (5) is numerically more efﬁcient than that of Eq. (4) (cf.
Tscherning et al., 1983). This is because each
X
m
in Eq. (3) is independent of
, and thus
need only be evaluated once for each parallel. If all such computation points are equally
spaced in longitude, further numerical efﬁciencies can be achieved through application
of the recursion algorithm developed by Rizos (1979). Abd-Elmotaal (1997) contains a
re-derivation of this algorithm which demonstrates that, contrary to the approach of Ri-
zos (1979), the algorithm can be applied in full without prior rotation of the geopotential
coefﬁcients.

##### Citations
More filters
Book ChapterDOI
01 Jan 2015
TL;DR: A general review of the mathematical formalism that is used in describing gravity and topography of the terrestrial planets is given in this article, where the basic properties of Earth, Venus, Mars, Mercury, and the Moon are characterized.
Abstract: This chapter reviews our current knowledge of the gravity and topography of the terrestrial planets and describes the methods that are used to analyze these data. A general review of the mathematical formalism that is used in describing gravity and topography is first given. Next, the basic properties of Earth, Venus, Mars, Mercury, and the Moon are characterized. Following this, the relationship between gravity and topography is quantified, and techniques by which geophysical parameters can be constrained are detailed. Analysis methods include crustal thickness modeling, geoid/topography ratios, spectral admittance and correlation functions, and localized spectral analysis and wavelet techniques. Finally, the major results that have been obtained by modeling the gravity and topography of Earth, Venus, Mars, Mercury, and the Moon are summarized.

301 citations

### Cites methods from "A unified approach to the Clenshaw ..."

• ...If starting values used in the recursion are appropriately scaled, as is summarized by Holmes and Featherstone (2002), these can be computed to high accuracy up to a maximum spherical harmonic degree of about 2700....

[...]

Journal ArticleDOI
TL;DR: In this paper, the Earth's topography harmonic expansion was used to obtain a spherical harmonic model of the topography of the Earth, which was then used for the estimation of induced gravity perturbations.
Abstract: The availability of high-resolution global digital elevation data sets has raised a growing interest in the feasibility of obtaining their spherical harmonic representation at matching resolution, and from there in the modelling of induced gravity perturbations. We have therefore estimated spherical Bouguer and Airy isostatic anomalies whose spherical harmonic models are derived from the Earth’s topography harmonic expansion. These spherical anomalies differ from the classical planar ones and may be used in the context of new applications. We succeeded in meeting a number of challenges to build spherical harmonic models with no theoretical limitation on the resolution. A specific algorithm was developed to enable the computation of associated Legendre functions to any degree and order. It was successfully tested up to degree 32,400. All analyses and syntheses were performed, in 64 bits arithmetic and with semi-empirical control of the significant terms to prevent from calculus underflows and overflows, according to IEEE limitations, also in preserving the speed of a specific regular grid processing scheme. Finally, the continuation from the reference ellipsoid’s surface to the Earth’s surface was performed by high-order Taylor expansion with all grids of required partial derivatives being computed in parallel. The main application was the production of a 1′ × 1′ equiangular global Bouguer anomaly grid which was computed by spherical harmonic analysis of the Earth’s topography–bathymetry ETOPO1 data set up to degree and order 10,800, taking into account the precise boundaries and densities of major lakes and inner seas, with their own altitude, polar caps with bedrock information, and land areas below sea level. The harmonic coefficients for each entity were derived by analyzing the corresponding ETOPO1 part, and free surface data when required, at one arc minute resolution. The following approximations were made: the land, ocean and ice cap gravity spherical harmonic coefficients were computed up to the third degree of the altitude, and the harmonics of the other, smaller parts up to the second degree. Their sum constitutes what we call ETOPG1, the Earth’s TOPography derived Gravity model at 1′ resolution (half-wavelength). The EGM2008 gravity field model and ETOPG1 were then used to rigorously compute 1′ × 1′ point values of surface gravity anomalies and disturbances, respectively, worldwide, at the real Earth’s surface, i.e. at the lower limit of the atmosphere. The disturbance grid is the most interesting product of this study and can be used in various contexts. The surface gravity anomaly grid is an accurate product associated with EGM2008 and ETOPO1, but its gravity information contents are those of EGM2008. Our method was validated by comparison with a direct numerical integration approach applied to a test area in Morocco–South of Spain (Kuhn, private communication 2011) and the agreement was satisfactory. Finally isostatic corrections according to the Airy model, but in spherical geometry, with harmonic coefficients derived from the sets of the ETOPO1 different parts, were computed with a uniform depth of compensation of 30 km. The new world Bouguer and isostatic gravity maps and grids here produced will be made available through the Commission for the Geological Map of the World. Since gravity values are those of the EGM2008 model, geophysical interpretation from these products should not be done for spatial scales below 5 arc minutes (half-wavelength).

284 citations

### Cites background or methods from "A unified approach to the Clenshaw ..."

• ...Wenzel single normalization (1998), or Horner’s schemes (in cos φ) on partial sums such as in Holmes and Featherstone (2002),...

[...]

• ...Figure 3a shows that Plm(x) becomes nonsignificant above some maximum order, a phenomenon used by Jekeli et al. (2007)—see below. On Fig. 3b it is clear that the Hlm’s increase regularly with l before becoming significant, then oscillate around a stable value; and the Plm’s behave similarly (starting at ∼cos mφ), which we will use in our algorithm. When using recursive relations on the Hlm(x), as it has become customary in geodesy, usual tricks, e.g. Wenzel single normalization (1998), or Horner’s schemes (in cos φ) on partial sums such as in Holmes and Featherstone (2002),...

[...]

• ...Figure 3a shows that Plm(x) becomes nonsignificant above some maximum order, a phenomenon used by Jekeli et al. (2007)—see below....

[...]

• ...When using recursive relations on the Hlm(x), as it has become customary in geodesy, usual tricks, e.g. Wenzel single normalization (1998), or Horner’s schemes (in cos ϕ) on partial sums such as in Holmes and Featherstone (2002), do not work at such level....

[...]

Journal ArticleDOI
TL;DR: In this article, the authors developed a compact and robust representation of the transformation from geodetic to Quasi-Dipole (QD), Apex, and Modified Apex coordinates, by fitting the QD coordinates to spherical harmonics.
Abstract: [1] Many structural and dynamical features of the ionized and neutral upper atmosphere are strongly organized by the geomagnetic field, and several magnetic coordinate systems have been developed to exploit this organization. Quasi-Dipole coordinates are appropriate for calculations involving horizontally stratified phenomena like height-integrated currents, electron densities, and thermospheric winds; Modified Apex coordinates are appropriate for calculations involving electric fields and magnetic field-aligned currents. The calculation of these coordinates requires computationally expensive tracing of magnetic field lines to their apexes. Interpolation on a precomputed grid provides faster coordinate conversions, but requires the overhead of a sufficiently fine grid, as well as finite differencing to obtain coordinate base vectors. In this paper, we develop a compact and robust representation of the transformation from geodetic to Quasi-Dipole (QD), Apex, and Modified Apex coordinates, by fitting the QD coordinates to spherical harmonics in geodetic longitude and latitude. With this representation, base vectors may be calculated directly from the expansion coefficients. For an expansion truncated at order 6, the fitted coordinates deviate from the actual coordinates by a maximum of 0.4°, and typically by 0.1°. The largest errors occur in the equatorial Atlantic region. Compared to interpolation on a pre-computed grid, the spherical harmonic representation is much more compact and produces smooth base vectors. An algorithm for efficiently and concurrently computing scalar and vector spherical harmonic functions is provided in the appendix. Computer code for producing the expansion coefficients and evaluating the fitted coordinates and base vectors is included in the auxiliary material.

167 citations

### Cites background or methods from "A unified approach to the Clenshaw ..."

• ...Holmes and Featherstone [2002] noted that the standard forward column method, applied at double precision, will produce numerical underflow for m > 1900; we therefore expect that our algorithm will produce reliable output for orders considerably higher than m = 72....

[...]

• ...The recursion relations given in equations (A6), (A7), (A8), and (A12) are adapted fromHolmes and Featherstone [2002, equations (13), (11), (15), and (21)]. The recursion coefficients are largely the same as those given in Holmes and Featherstone [2002]; we have relabeled their f as d, and named a new coefficient c....

[...]

Journal ArticleDOI
TL;DR: SHTools as mentioned in this paper is an open-source archive of both Fortran 95 and Python routines for performing spherical harmonic analyses, including spherical harmonic transforms and spectral analysis of global gravity and magnetic fields.
Abstract: Geophysical analyses are often performed in spherical geometry and require the use of spherical harmonic functions to express observables or physical quantities. When expanded to high degree, the accuracy and speed of the spherical harmonic transforms and reconstructions are of paramount importance. SHTools is a time and user-tested open-source archive of both Fortran 95 and Python routines for performing spherical harmonic analyses. The routines support all spherical-harmonic normalization conventions used in the geosciences, including 4p-normalized, Schmidt seminormalized, orthonormalized, and unnormalized harmonics, along with the option of employing the Condon-Shortley phase factor of ð21Þ m. Data on the sphere can be sampled on a variety of grid formats, including equally spaced cylindrical grids and grids appropriate for integration by Gauss-Legendre quadrature. The spherical-harmonic transforms are proven to be fast and accurate for spherical harmonic degrees up to 2800. Several tools are provided for the geoscientist, including routines for performing localized spectral analyses and basic operations related to global gravity and magnetic fields. In the Python environment, operations are very simple to perform as a result of three class structures that encompass all operations on grids, spherical harmonic coefficients, and spatiospectral localization windows. SHTools is released under the unrestrictive BSD 3-clause license.

158 citations

Journal ArticleDOI

141 citations

##### References
More filters
01 Jan 1960
TL;DR: In this paper, the angular momentum, one of the most fundamental quantities in all of quantum mechanics, is introduced and a concise introduction to its application in atomic, molecular, and nuclear physics is provided.
Abstract: This book offers a concise introduction to the angular momentum, one of the most fundamental quantities in all of quantum mechanics. Beginning with the quantization of angular momentum, spin angular momentum, and the orbital angular momentum, the author goes on to discuss the Clebsch-Gordan coefficients for a two-component system. After developing the necessary mathematics, specifically spherical tensors and tensor operators, the author then investigates the 3-j, 6-j, and 9-j symbols. Throughout, the author provides practical applications to atomic, molecular, and nuclear physics. These include partial-wave expansions, the emission and absorption of particles, the proton and electron quadrupole moment, matrix element calculation in practice, and the properties of the symmetrical top molecule.

5,050 citations

Book
21 Sep 1957
TL;DR: In this article, the angular momentum, one of the most fundamental quantities in all of quantum mechanics, is introduced and a concise introduction to its application in atomic, molecular, and nuclear physics is provided.
Abstract: This book offers a concise introduction to the angular momentum, one of the most fundamental quantities in all of quantum mechanics. Beginning with the quantization of angular momentum, spin angular momentum, and the orbital angular momentum, the author goes on to discuss the Clebsch-Gordan coefficients for a two-component system. After developing the necessary mathematics, specifically spherical tensors and tensor operators, the author then investigates the 3-j, 6-j, and 9-j symbols. Throughout, the author provides practical applications to atomic, molecular, and nuclear physics. These include partial-wave expansions, the emission and absorption of particles, the proton and electron quadrupole moment, matrix element calculation in practice, and the properties of the symmetrical top molecule.

4,377 citations

### "A unified approach to the Clenshaw ..." refers background in this paper

• ...Risbo (1996) claims that the Edmonds (1957) recursion for D-matricies can be used to compute fully normalised ALFs up to degree 200,000....

[...]

Book
01 Jan 1966

3,150 citations

01 Jul 1998
TL;DR: The Earth Gravitational Model 1996 (EGM96) as discussed by the authors was developed by the NASA Goddard Space Flight Center (GSFC), the National Imagery and Mapping Agency (NIMA), and The Ohio State University (OSU) to develop an improved spherical harmonic model of the Earth's gravitational potential.
Abstract: The NASA Goddard Space Flight Center (GSFC), the National Imagery and Mapping Agency (NIMA), and The Ohio State University (OSU) have collaborated to develop an improved spherical harmonic model of the Earth's gravitational potential to degree 360. The new model, Earth Gravitational Model 1996 (EGM96), incorporates improved surface gravity data, altimeter-derived gravity anomalies from ERS-1 and from the GEOSAT Geodetic Mission (GM), extensive satellite tracking data-including new data from Satellite Laser Ranging (SLR), the Global Postioning System (GPS), NASA's Tracking and Data Relay Satellite System (TDRSS), the French DORIS system, and the US Navy TRANET Doppler tracking system-as well as direct altimeter ranges from TOPEX/POSEIDON (T/P), ERS-1, and GEOSAT. The final solution blends a low-degree combination model to degree 70, a block-diagonal solution from degree 71 to 359, and a quadrature solution at degree 360. The model was used to compute geoid undulations accurate to better than one meter (with the exception of areas void of dense and accurate surface gravity data) and realize WGS84 as a true three-dimensional reference system. Additional results from the EGM96 solution include models of the dynamic ocean topography to degree 20 from T/P and ERS-1 together, and GEOSAT separately, and improved orbit determination for Earth-orbiting satellites.

832 citations

### "A unified approach to the Clenshaw ..." refers methods in this paper

• ...For example, all methods will report underflows (8 ) when EGM96 coefficients are employed for the lower degrees....

[...]

• ...Lemoine FG, Kenyon SC, Factor JK, Trimmer RG, Pavlis NK, Chinn DS, Cox CM, Klosko SM, Luthcke SB, Torrence MH, Wang YM, Williamson RG, Pavlis EC, Rapp RH, Olson TR (1998) The development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96, NASA/TP-1998206861, National Aeronautics and Space Administration, Washington....

[...]

• ...For Enm , one might use empirically generated coefficients such as EGM96 (Lemoine et al., 1998) and/or GPM98B (Wenzel, 1998) to compute the lower degreeEnm ....

[...]

• ...For Enm , one might use empirically generated coefficients such as EGM96 (Lemoine et al., 1998) and/or GPM98B (Wenzel, 1998) to compute the lower degree Enm ....

[...]

Journal ArticleDOI

820 citations