Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers: II first-, second-, and third-order derivatives (original) (raw)
Related papers
By extending the exponent of floating point numbers with an additional integer as the power index of a large radix, we compute fully normalized associated Legendre functions (ALF) by recursion without underflow problem. The new method enables us to evaluate ALFs of extremely high degree as 2 32 = 4,294,967,296, which corresponds to around 1 cm resolution on the Earth's surface. By limiting the application of exponent extension to a few working variables in the recursion, choosing a suitable large power of 2 as the radix, and embedding the contents of the basic arithmetic procedure of floating point numbers with the exponent extension directly in the program computing the recurrence formulas, we achieve the evaluation of ALFs in the doubleprecision environment at the cost of around 10% increase in computational time per single ALF. This formulation realizes meaningful execution of the spherical harmonic synthesis and/or analysis of arbitrary degree and order.
On computation and use of Fourier coefficients for associated Legendre functions
Journal of Geodesy, 2016
The computation of spherical harmonic series in very high resolution is known to be delicate in terms of performance and numerical stability. A major problem is to keep results inside a numerical range of the used data type during calculations as under-/overflow arises. Extended data types are currently not desirable since the arithmetic complexity will grow exponentially with higher resolution levels. If the associated Legendre functions are computed in spectral domain then regular grid transformations can be applied highly efficiently and convenient for derived quantities as well. In this article we compare three recursive computations of the associated Legendre functions as trigonometric series, thereby ensuring a defined numerical range for each constituent wave-number, separately. The results to high degree and order show the numerical strength of the proposed method. First, the evaluation of Fourier coefficients of the associated Legendre functions has been done with respect to the floating-point precision requirements. Secondly, the numerical accuracy in the cases of standard Double and long Double precision arithmetic is demonstrated. Following Bessel's inequality the obtained accuracy estimates of the Fourier coefficients are directly transferable to the associated Legendre functions themselves and to derived functionals as well. Therefore, they can provide an essential insight to modern geodetic applications that depend on efficient spherical harmonic analysis and synthesis beyond [5 × 5] arcmin resolution.
Ultra-high degree spherical harmonic analysis and synthesis using extended-range arithmetic
Journal of Geodesy, 2008
We present software for spherical harmonic analysis (SHA) and spherical harmonic synthesis (SHS), which can be used for essentially arbitrary degrees and all co-latitudes in the interval (0 • , 180 •). The routines use extended-range floating-point arithmetic, in particular for the computation of the associated Legendre functions. The price to be paid is an increased computation time; for degree 3, 000, the extended-range arithmetic SHS program takes 49 times longer than its standard arithmetic counterpart. The extended-range SHS and SHA routines allow us to test existing routines for SHA and SHS. A comparison with the publicly available SHS routine GEOGFG18 by Wenzel and HARMONIC SYNTH by Holmes and Pavlis confirms what is known about the stability of these programs. GEOGFG18 gives errors <1 mm for latitudes [−89 • 57.5 , 89 • 57.5 ] and maximum degree 1, 800. Higher degrees significantly limit the range of acceptable latitudes for a given accuracy. HARMONIC SYNTH gives good results up to degree 2, 700
A recursive method is developed to compute the ratios of the oblate spheroidal harmonics of the second kind and their first-, second-, and third-order derivatives. The recurrence formulas consist of three kinds: (1) fixeddegree increasing-order, (2) mixed-degree increasing-order, and (3) fixed-order decreasing-degree. The three seed values are evaluated by rapidly convergent series. The derivatives of the ratios are recursively obtained from the values and lower-order derivatives of the same harmonic order and of the same or higher degrees. The new method precisely and quickly computes the ratios and their low-order derivatives. It provides 13 correct digits of the ratios of degree as high as 262,000 and runs 20-100 times faster than the existing methods.
Journal of Geodesy, 2018
n order to accelerate the spherical harmonic synthesis and/or analysis of arbitrary function on the unit sphere, we developed a pair of procedures to transform between a truncated spherical harmonic expansion and the corresponding two-dimensional Fourier series. First, we obtained an analytic expression of the sine/cosine series coefficient of the 4π fully normalized associated Legendre function in terms of the rectangle values of the Wigner d function. Then, we elaborated the existing method to transform the coefficients of the surface spherical harmonic expansion to those of the double Fourier series so as to be capable with arbitrary high degree and order. Next, we created a new method to transform inversely a given double Fourier series to the corresponding surface spherical harmonic expansion. The key of the new method is a couple of new recurrence formulas to compute the inverse transformation coefficients: a decreasing-order, fixed-degree, and fixed-wavenumber three-term formula for general terms, and an increasing-degree-and-order and fixed-wavenumber two-term formula for diagonal terms. Meanwhile, the two seed values are analytically prepared. Both of the forward and inverse transformation procedures are confirmed to be sufficiently accurate and applicable to an extremely high degree/order/wavenumber as 230 ≈109. The developed procedures will be useful not only in the synthesis and analysis of the spherical harmonic expansion of arbitrary high degree and order, but also in the evaluation of the derivatives and integrals of the spherical harmonic expansion.
Discrete Spherical Harmonic Transforms: Numerical Preconditioning and Optimization
Lecture Notes in Computer Science, 2008
Spherical Harmonic Transforms (SHTs) which are essentially Fourier transforms on the sphere are critical in global geopotential and related applications. Among the best known strategies for discrete SHTs are Chebychev quadratures and least squares. The numerical evaluation of the Legendre functions are especially challenging for very high degrees and orders which are required for advanced geocomputations. The computational aspects of SHTs and their inverses using both quadrature and least-squares estimation methods are discussed with special emphasis on numerical preconditioning that guarantees reliable results for degrees and orders up to 3800 in REAL*8 or double precision arithmetic. These numerical results of spherical harmonic synthesis and analysis using simulated spectral coefficients are new and especially important for a number of geodetic, geophysical and related applications with ground resolutions approaching 5 km.
Journal of Geodetic Science, 2018
In order to accelerate the spherical/spheroidal harmonic synthesis of any function, we developed a new recursive method to compute the sine/cosine series co-e cient of the π fully-and Schmidt quasi-normalized associated Legendre functions. The key of the method is a set of increasing-degree/order mixed-wavenumber two-to four-term recurrence formulas to compute the diagonal terms. They are used in preparing the seed values of the decreasing-order xed-degree, and xed-wavenumber two-and three-term recurrence formulas, which are obtained by modifying the classic relations. The new method is accurate and capable to deal with an arbitrary high de-gree/order/wavenumber. Also, it runs signi cantly faster than the previous method of ours utilizing the Wigner d function, say around 20 times more when the maximum degree exceeds 1,000.
Vector spherical harmonics: Concepts and applications to the single centre expansion method
Computer Physics Communications, 2000
In this paper the basic formulation of the Vector Spherical Harmonics (VSH) is recalled from the literature and the details on the various notations clarified. By means of compact definitions of differential operators in terms of VSH, a new formulation for the gradient and Laplacian of a generic Single Centre Expanded (SCE) function is developed. The new formulae are so derived to be directly computable and implementable in numerical codes to face problems in atomic and molecular physics. A set of numerical computations have been carried out to assess the proposed formulae where the behaviour of the gradient and Laplacian of the electron density and static potential of the water molecule is exploited in detail within the SCE framework of reference. : S 0 0 1 0 -4 6 5 5 ( 0 0 ) 0 0 1 3 7 -5
Fortran-and C-codes for higher order and degree geopotential and derivatives computation
Geopotential computation in general did not present numerical problems in older models where order and degree were as low as 30. However with the appearance of much higher order and degree the older conventional algorithms could not be used anymore. The kernel of such algorithms is the recursive computation of the Legendre polynomials. They must be computed in their normalized recursions and variants to avoid numerical flaws at a cost of more expensive processing times mainly due to the need of square roots computations. This paper outlines the main implementation aspects of the geopotential computation and derivatives to higher order and degree. The codes were tested up to 2159 order and 2190 degree without any noticeable flaw. It is expected some numerical degradation near the poles, however to the 0.000001° of proximity, i.e. ±89.999999° latitude, no problems are reported. Computer codes were developed in double precision and are freely available in Fortran-77 and ANSI-C languages.