charm_shs

Module to perform solid spherical harmonic synthesis of point and mean values.

Note

This documentation is written for double precision version of CHarm.

Synthesis of point values

Functions to perform spherical harmonic synthesis of point values and of their first- and second-order gradients in LNOF.

void charm_shs_point(const charm_point *pnt, const charm_shc *shcs, unsigned long nmax, double *f, charm_err *err)

Performs the exact solid spherical harmonic synthesis of point values:

\[f(r, \varphi,\lambda) = \dfrac{\mu}{R} \, \sum_{n = 0}^{n_{\max}} \left( \dfrac{R}{r} \right)^{n + 1} \, \sum_{m = 0}^{n} \left( \bar{C}_{nm}\, \cos(m \, \lambda) + \bar{S}_{nm} \, \sin(m \, \lambda) \right) \, \bar{P}_{nm}(\sin\varphi){.}\]

The synthesis is done either at

• a grid if pnt->type is CHARM_CRD_POINT_GRID, CHARM_CRD_POINT_GRID_GL, CHARM_CRD_POINT_GRID_DH1 or CHARM_CRD_POINT_GRID_DH2 (see charm_point), or

• scattered points if pnt->type is CHARM_CRD_POINT_SCATTERED.

For grids, efficient 1D FFT-based algorithm is applied along the latitude parallels whenever possible. Otherwise, the Chebyshev recurrences (e.g., Balmino et al 2012) are employed along the parallels. The latter option is slower. The following conditions must be satisfied in order for CHarm to apply FFT:

• pnt->nlon > 1,

• (pnt->nlon - 1) / 2 >= nmax, where the division is rounded down,

• the longitudinal step dlon is constant,

• pnt->lon[0] == 0.0 within the numerical threshold (charm_glob_threshold), and

• pnt->lon[pnt->nlon - 1] + dlon == 2.0 * PI within the numerical threshold (charm_glob_threshold).

If any of these conditions is not satisfied, CHarm uses the Chebyshev recurrences. The algorithm is selected automatically by CHarm with the preference being FFT.

For scattered points, tedious point-by-point synthesis is applied.

The synthesis is parallelized using OpenMP.

Notes on evaluation points organized in grids (not relevant for scattered points)

• To improve the computational speed, exploited is the symmetry of Legendre functions with respect to the equator. This is done automatically, provided that the grid stored in pnt is recognized as symmetric with respect to the equator. Symmetric and non-symmetric grids are explained here.

• The longitudes of the grid must be sampled increasingly with equal spacing, that is, the difference between any two consecutive longitudes must be positive and constant. The function performs a check on this. If the longitudes do not pass the test, an error message is written to err and the program returns to the caller.

References:

• Sneeuw, N. (1994) Global spherical harmonic analysis by least-squares and numerical quadrature methods in historical perspective. Geophysical Journal International 118:707-716

• Heiskannen, W. A., Moritz, H. (1967) Physical Geodesy. W. H. Freeman and Company, San Francisco, 364 pp

• Fukushima, T. (2012) Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers. Journal of Geodesy 86:271-285.

• Balmino, Vales, Bonvalot, Briais, 2012. Spherical harmonic modelling to ultra-high degree of Bouguer and isostatic anomalies. Journal of Geodesy 86:499-520, doi: 10.1007/s00190-011-0533-4

Note

The synthesis always starts at degree nmin = 0. To use some other value of nmin (but satisfying 1 <= nmin <= nmax), you can set the respective coefficients in shcs to zero and use the newly modified structure shcs as the input to this function. The synthesis will still formally start at degree 0, but since the coefficients up to degree nmin - 1 are zero, they do not contribute to the output signal as required.

Note

The function does not use pnt->w, so no constrains are applied.

Parameters:

• pnt – [in] Evaluation points, organized either as a grid or as scattered points (see charm_point).

• shcs – [in] Spherical harmonic coefficients (see charm_shc).

• nmax – [in] Maximum spherical harmonic degree of the synthesis.

• f – [out] Pointer to an output array with the synthesized signal. If pnt holds a grid, the pointer f must have an access to pnt->nlat * pnt->nlon array elements. The value of the signal synthesized at pnt->lat[i] and pnt->lon[j] can be found as: f[i * pnt->nlon + j] with i = 0, 1, ..., pnt->nlat - 1 and j = 0, 1, ..., pnt->nlon 1. In case pnt represents scattered points, f must have an access to pnt->nlat == pnt->nlon elements and f[i] stands for the value synthesized at pnt->lat[i] and pnt->lon[i] with i = 0, 1, ..., pnt->nlat - 1 .

• err – [out] Error reported by the function (if any).

void charm_shs_point_grad1(const charm_point *pnt, const charm_shc *shcs, unsigned long nmax, double **f, charm_err *err)

Performs the synthesis of point values of the first-order gradient in LNOF,

\[\begin{split}\nabla f = % \begin{bmatrix} f_{\mathrm{x}}\\ f_{\mathrm{y}}\\ f_{\mathrm{z}}\\ \end{bmatrix} % = % \begin{bmatrix} \dfrac{1}{r} \, \dfrac{\partial f}{\partial \varphi}\\ -\dfrac{1}{r \, \cos\varphi} \, \dfrac{\partial f}{\partial \lambda}\\ \dfrac{\partial f}{\partial r}\\ \end{bmatrix} {.}\end{split}\]

If \(f\) is, for instance, the gravitational potential at \((r, \varphi, \lambda)\), then \(\nabla f\) is the vector of the gravitational acceleration at \((r, \varphi, \lambda)\).

Tip

It is more efficient to get the three elements of the first-order gradient in a single run by this function then by three separate calls of charm_shs_point_guru() yielding the three first-order derivatives with respect to the spherical coordinates. On the other hand, charm_shs_point_guru() is a better choice if you need only a single element of the gradient vector.

Warning

In polar areas, the output values of \(f_{\mathrm{x}}\) and \(f_{\mathrm{y}}\) are likely inaccurate or wrong, because (i) CHarm evaluates the latitudinal derivatives of Legendre functions by recurrences that are singular at the poles and (ii) the \(\frac{1}{\cos\varphi}\) term is singular at the poles.

Parameters:

• pnt – [in] Same as in charm_shs_point().

• shcs – [in] Same as in charm_shs_point().

• nmax – [in] Same as in charm_shs_point().

• f – [out] Double pointer to double with an access to at least three pointers to double. The x, y and z elements are stored in arrays pointed to by f[0], f[1] and f[2], respectively.

• err – [out] Same as in charm_shs_point().

void charm_shs_point_grad2(const charm_point *pnt, const charm_shc *shcs, unsigned long nmax, double **f, charm_err *err)

Performs the synthesis of point values of the second-order gradient in LNOF,

\[\begin{split}\nabla \otimes \nabla f = % \begin{bmatrix} f_{\mathrm{xx}}, f_{\mathrm{xy}}, f_{\mathrm{xz}}\\ f_{\mathrm{yx}}, f_{\mathrm{yy}}, f_{\mathrm{yz}}\\ f_{\mathrm{zx}}, f_{\mathrm{zy}}, f_{\mathrm{zz}}\\ \end{bmatrix} {.}\end{split}\]

If \(f\) is, for instance, the gravitational potential at \((r, \varphi, \lambda)\), then \(\nabla \otimes \nabla f\) is the gravitational tensor at \((r, \varphi, \lambda)\).

Due to the properties of harmonic fields, only six elements of the tensor are unique. In CHarm, these are computed after Reed (1973) by the following equations (note that the \(y\)-axis in Reed 1973 has the opposite direction as here which explains the sign swaps)

\[\begin{split}f_{\mathrm{xx}} &= \frac{1}{r^2} \, \frac{\partial^2 f}{\partial \varphi^2} + \frac{1}{r} \, \frac{\partial f}{\partial r},\\ % f_{\mathrm{xy}} &= -\frac{1}{r^2 \, \cos\varphi} \, \frac{\partial^2 f}{\partial \varphi \, \partial\lambda} - \frac{\tan\varphi}{r^2 \, \cos\varphi} \, \frac{\partial f}{\partial \lambda},\\ % f_{\mathrm{xz}} &= \frac{1}{r} \, \frac{\partial^2 f}{\partial r \, \partial \varphi} - \frac{1}{r^2} \, \frac{\partial f}{\partial \varphi},\\ % f_{\mathrm{yy}} &= \frac{1}{r^2 \, \cos^2\varphi} \, \frac{\partial^2 f}{\partial \lambda^2} + \frac{1}{r} \, \frac{\partial f}{\partial r} - \frac{\tan\varphi}{r^2} \, \frac{\partial f}{\partial \varphi},\\ % f_{\mathrm{yz}} &= -\frac{1}{r \, \cos\varphi} \, \frac{\partial^2 f}{\partial r \, \partial \lambda} + \frac{1}{r^2 \, \cos\varphi} \, \frac{\partial f}{\partial \lambda},\\ % f_{\mathrm{zz}} &= \frac{\partial^2 f}{\partial r^2}.\end{split}\]

Tip

It is significantly more efficient to get the six elements of the tensor in a single run by this function then by nine separate calls of charm_shs_point_guru() yielding the nine derivatives with respect to the spherical coordinates that are needed to get the derivatives in LNOF. On the other hand, charm_shs_point_guru() may be a better choice if you need only a single tensor element.

References:

• Reed, G. B. (1973) Application of kinematical geodesy for determining the short wavelength components of the gravity field by satellite gradiometry. Report No. 201, Ohio State University, Department of Geodetic Sciences, Columbus, USA

Warning

In polar areas, all output values except for \(f_{\mathrm{zz}}\) are likely inaccurate or wrong, because (i) CHarm evaluates the latitudinal derivatives of Legendre functions by recurrences that are singular at the poles and (ii) the \(\frac{1}{\cos\varphi}\) and \(\frac{1}{\cos^2\varphi}\) terms are singular at the poles.

Parameters:

• pnt – [in] Same as in charm_shs_point().

• shcs – [in] Same as in charm_shs_point().

• nmax – [in] Same as in charm_shs_point().

• f – [out] Double pointer to double with an access to at least six pointers to double. The xx, xy, xz, yy, yz and zz elements are stored in arrays pointed to by f[0], f[1], ..., f[5], respectively.

• err – [out] Same as in charm_shs_point().

void charm_shs_point_guru(const charm_point *pnt, const charm_shc *shcs, unsigned long nmax, unsigned dr, unsigned dlat, unsigned dlon, double *f, charm_err *err)

Performs the synthesis of

\[\frac{1}{r^{j + k} \, \cos^k\varphi} \, \frac{\partial^{i + j + k} f}{\partial r^i \, \partial \varphi^j \, \lambda^k}\]

for \(i = 0, 1, 2\), \(j = 0, 1, 2\) and \(k = 0, 1, 2\) satisfying \(i + j + k \leq 2\).

Computed can be, for instance, quantities like \(\dfrac{1}{r \, \cos\varphi} \, \dfrac{\partial f}{\partial \lambda}\) ( \(i = 0\), \(j = 0\), \(k = 1\)) or \(\dfrac{1}{r^2 \, \cos\varphi} \, \dfrac{\partial^2 f}{\partial \varphi \, \partial \lambda}\) ( \(i = 0\), \(j = 1\), \(k = 1\)).

The rational behind this function is twofold.

• It lowers memory requirements when compared with charm_shs_point_grad1() and charm_shs_point_grad2(). For instance, the evaluation of the full tensor by charm_shs_point_grad2() requires a huge amount of RAM if the point grid is really dense and/or maximum harmonic degree is large. With this function and some additional effort, the outputs from charm_shs_point_grad1() and charm_shs_point_grad2() can be reconstructed with much lower memory requirements. This will be significantly slower, but in some cases it may be the only practical solution.

• These quantities are sometimes needed to perform some specialized tasks, for instance, as in cap-modified spectral gravity forward modelling (Bucha et al, 2019).

References

Bucha, B., Hirt, C., Kuhn, M. (2019). Cap integration in spectral gravity forward modelling up to the full gravity tensor. Journal of Geodesy 93, 1707-1737, doi: 10.1007/s00190-019-01277-3

Note

The sum dr + dlat + dlon cannot be larger than 2.

Warning

In polar areas, the outputs are likely inaccurate or wrong if dlat > 0 and/or dlon > 0, because (i) CHarm evaluates the latitudinal derivatives of Legendre functions by recurrences that are singular at the poles and (ii) the \(\frac{1}{\cos\varphi}\) and \(\frac{1}{\cos^2\varphi}\) terms are singular at the poles.

Parameters:

• pnt – [in] Same as in charm_shs_point().

• shcs – [in] Same as in charm_shs_point().

• nmax – [in] Same as in charm_shs_point().

• dr – [in] Order of the radial derivative (variable \(i\) in the equation above).

• dlat – [in] Order of the latitudinal derivative (variable \(j\) in the equation above).

• dlon – [in] Order of the longitudinal derivative (variable \(k\) in the equation above).

• f – [out] Same as in charm_shs_point(), but the output quantity now represents one of the derivatives shows above.

• err – [out] Same as in charm_shs_point().

Synthesis of mean values

Functions to perform spherical harmonic synthesis of mean values.

void charm_shs_cell(const charm_cell *cell, const charm_shc *shcs, unsigned long nmax, double *f, charm_err *err)

Performs the exact solid spherical harmonic synthesis of mean values over computational cells:

\[\begin{split}\tilde{f}(r, \varphi_{\mathrm{min}}, \varphi_{\mathrm{max}}, \lambda_{\mathrm{min}}, \lambda_{\mathrm{max}}) =& \dfrac{1}{\Delta \sigma} \, \dfrac{\mu}{R} \, \sum_{n = 0}^{n_{\max}} \left( \dfrac{R}{r} \right)^{n + 1} \\ &\times\sum_{m = 0}^{n} \left( \bar{C}_{nm}\, \int\limits_{\lambda_{\mathrm{min}}}^{\lambda_{\mathrm{max}}} \cos(m \, \lambda) \, \mathrm{d}\lambda + \bar{S}_{nm} \, \int\limits_{\lambda_{\mathrm{min}}}^{\lambda_{\mathrm{max}}} \sin(m \, \lambda) \, \mathrm{d}\lambda \right) \\ &\times\int\limits_{\varphi_{\mathrm{min}}}^{\varphi_{\mathrm{max}}} \bar{P}_{nm}(\sin\varphi) \, \cos\varphi \, \mathrm{d}\varphi{,}\end{split}\]

The synthesis is done either at

• a grid if cell->type is CHARM_CRD_CELL_GRID (see charm_cell), or

• scattered cells if cell->type is CHARM_CRD_CELL_SCATTERED.

For grids, efficient 1D FFT-based algorithm is applied along the latitude parallels whenever possible. Otherwise, the Chebyshev recurrences are employed along the parallels. The latter option is slower. The following conditions must be satisfied in order for CHarm to apply FFT:

• cell->nlon > 1,

• (cell->nlon - 1) / 2 >= nmax, where the division is rounded down,

• the longitudinal step dlon is constant,

• cell->lonmin[0] == 0.0 within the numerical threshold (charm_glob_threshold),

• cell->lonmax[2 * cell->nlon - 1] == 2.0 * PI within the numerical threshold (charm_glob_threshold), and

• the maximum longitude of the jth cell must be equal to the minimum longitude of the j + 1th cell.

If any of these conditions is not satisfied, CHarm uses the Chebyshev recurrences. The algorithm is selected automatically by CHarm with the preference being FFT.

For scattered cells, tedious cell-by-cell synthesis is applied.

The synthesis is parallelized using OpenMP.

Notes on evaluation cells organized in grids (not relevant for scattered cells)

• To improve the computational speed, exploited is the symmetry of Legendre functions with respect to the equator. This is done automatically, provided that the grid stored in cell is recognized as symmetric with respect to the equator. Symmetric and non-symmetric grids are explained here.

• The longitudes of the grid must be sampled increasingly with equal spacing, that is, the difference between any two consecutive minimum cell longitudes must be positive and constant and the same holds for the maximum cell longitudes as well; moreover, the step in minimum cell longitudes must be equal to the step in maximum cell longitudes. Two examples of an acceptable longitudinal sampling (shown are minimum and maximum longitudes of grid cells):

  lonmin: 0.0, 1.0, 2.0, 3.0, 4.0, 5.0
  lonmax: 1.0, 2.0, 3.0, 4.0, 5.0, 6.0
  

  lonmin: 0.0, 0.5, 1.0, 1.5, 2.0, 2.5
  lonmax: 1.0, 1.5, 2.0, 2.5, 3.0, 3.5
  

  • The function performs all necessary checks. If the longitudes do not pass the test, an error message is written to err and program returns to the caller.

References:

• Colombo, O. L. (1981) Numerical methods for harmonic analysis on the sphere. Reports of the Department of Geodetic Science, Report No. 310, 140 pp

• Heiskannen, W. A., Moritz, H. (1967) Physical Geodesy. W. H. Freeman and Company, San Francisco, 364 pp

• Jekeli, C., Lee, J. K., Kwon, J. H. (2007) On the computation and approximation of ultra-high-degree spherical harmonic series. Journal of Geodesy 81:603-615

• Fukushima, T. (2012) Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers. Journal of Geodesy 86:271-285.

• Balmino, Vales, Bonvalot, Briais, 2012. Spherical harmonic modelling to ultra-high degree of Bouguer and isostatic anomalies. Journal of Geodesy 86:499–520, doi: 10.1007/s00190-011-0533-4

• Fukushima, T. (2014) Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers: III integral. Computers & Geosciences 63:17-21

Note

The synthesis automatically starts at degree nmin = 0. To use some other value of nmin, see the tip given in charm_shs_point().

Parameters:

• cell – [in] Evaluation cells, organized either as a grid or as scattered cells (see charm_cell).

• shcs – [in] Spherical harmonic coefficients (see charm_shc).

• nmax – [in] Maximum spherical harmonic degree of the synthesis.

• f – [out] Pointer to an output array with the synthesized mean values. If cell holds a grid, the pointer f must have an access to cell->nlat * cell->nlon array elements. The mean value of the signal synthesized at the ith cell in the latitudinal direction and jth cell in the longitudinal direction can be found as: f[i * cell->nlon + j] with i = 0, 1, ..., cell->nlat - 1 and j = 0, 1, ..., cell->nlon - 1. In case cell represents scattered cells, f must have an access to cell->nlat == cell->nlon elements and f[i] stands for the mean value synthesized at the ith cell with i = 0, 1, ..., cell->nlat - 1.

• err – [out] Error reported by the function (if any).

void charm_shs_cell_isurf(const charm_cell *cell, const charm_shc *shcs1, unsigned long nmax1, const charm_shc *shcs2, unsigned long nmax2, unsigned long nmax3, unsigned long nmax4, double *f, charm_err *err)

Similarly as charm_shs_cell(), this function computes mean values of a function given by shcs1. The difference is that here, the mean values are integrated not on the unit sphere, but instead on a band-limited irregular surface \(r(\varphi,\lambda)\) defined by spherical harmonic coefficients in shcs2:

\[\begin{split}\tilde{f}(r, \varphi_{\mathrm{min}}, \varphi_{\mathrm{max}},\lambda_{\mathrm{min}},\lambda_{\mathrm{max}}) =& \frac{1}{\Delta \sigma} \, \frac{\mu^{(1)}}{R^{(1)}} \, \int\limits_{\varphi = \varphi_{\mathrm{min}}}^{\varphi_{\mathrm{max}}} \int\limits_{\lambda = \lambda_{\mathrm{min}}}^{\lambda_{\mathrm{max}}}\sum_{n = 0}^{n_{\max_1}} \left( \frac{R^{(1)}}{r(\varphi, \lambda)}\right)^{(n + 1)} \\ &\times \sum_{m = 0}^{n} \left( \bar{C}_{nm}^{(1)}\, \cos(m \, \lambda) + \bar{S}_{nm}^{(1)} \, \sin(m \, \lambda) \right) \, \bar{P}_{nm}(\sin\varphi)\\ &\times \mathrm{d}\lambda \, \cos\varphi \, \mathrm{d}\varphi\end{split}\]

with

\[r(\varphi, \lambda) = \sum_{n = 0}^{n_{\max_2}} \sum_{m = 0}^{n} \left( \bar{C}_{nm}^{(2)}\, \cos(m \, \lambda) + \bar{S}_{nm}^{(2)} \, \sin(m \, \lambda) \right) \, \bar{P}_{nm}(\sin\varphi){.}\]

This function can only be used with evaluation cells organized in a grid, that is, cell->type is set to CHARM_CRD_CELL_GRID. Currently, only the Chebyshev recurrences are supported for computation along the latitude parallels.

The Fourier coefficients are evaluated after Fukushima (2018).

The function is parallelized using OpenMP.

Notes on evaluation cells organized in grids

• The symmetry of Legendre function is not employed in this function. Numerical experiments have shown that, in our implementation, the speed improvement is completely negligible if any.

• The structure of the longitudes must be the same as discussed in charm_shs_cell().

References:

• Fukushima, T. (2012) Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers. Journal of Geodesy 86:271-285.

• Fukushima, T. (2018) Fast computation of sine/cosine series coefficients of associated Legendre functions of arbitrary high degree and order.

Warning

Contrary to charm_shs_point() and charm_shs_cell(), the computation is not exact, but the output accuracy can be well controled by the input variables.

Note

The function does not use spherical radii cell->r, since these are synthesized from shcs2.

Note

The synthesis automatically starts at degree nmin = 0. To use some other value of nmin, see the tip given in charm_shs_point().

Warning

Both scaling constants of the shcs2 structure that defines the irregular surface (shcs2->mu and shcs2->r) have to be equal to 1.0.

Warning

This function is very demanding on RAM. More than 4 * (nmax1 + 1)^2 * (nmax3 + 1)^2 * 8 / 1024^3 GBs of memory are required in double precision. For instance, for nmax1 = 150 and nmax3 = 150, about 15 GBs of RAM are occupied during the computation.

Parameters:

• cell – [in] Evaluation cells organized as a grid (see charm_cell).

• shcs1 – [in] Spherical harmonic coefficients of the function (e.g., the gravitational potential), the block mean values of which are synthesized (see charm_shc).

• nmax1 – [in] Maximum harmonic degree of the function given by shcs1, the block mean values of which are synthesized.

• shcs2 – [in] Spherical harmonic coefficients of the surface, on which the block mean values are synthesized (see charm_shc). The surface must be defined by the spherical radii r, that is, the synthesis of shcs2 yields spherical radii of the surface, on which the mean values are sought.

• nmax2 – [in] Maximum harmonic degree of the surface defined by shcs2, on which the block mean values are synthesized.

• nmax3 – [in] Maximum harmonic degree of the powers (shcs1->r / r)^(n + 1) used in the synthesis of the mean values. The term (shcs1->r / r)^(n + 1) is a 2D sphere-like surface and nmax3 defines the maximum harmonic degree of this surface to be used in the synthesis. In theory, nmax3 should be infinite, since shcs1->r / r is non-band-limited for r band-limited to nmax2.

• nmax4 – [in] Maximum harmonic degree, up to which the powers (shcs1->r / r)^(n + 1) are harmonically analyzed. The higher is the value of nmax4, the smaller is the aliasing effect. In theory, nmax4 should be infinite, since shcs1->r / r is non-band-limited for band-limited r. Note that nmax4 cannot be smaller than nmax3. In other words, the term (shcs1->r / r)^(n + 1) is harmonically analyzed up to degree nmax4, but to synthesize it later, the nmax3 maximum degree is used.

• f – [out] Pointer to an output array with the synthesized mean values. The pointer f must have an access to cell->nlat * cell->nlon array elements. The mean value of the signal synthesized at the ith cell in the latitudinal direction and jth cell in the longitudinal direction can be found as: f[i * cell->nlon + j] with i = 0, 1, ..., cell->nlat - 1 and j = 0, 1, ..., cell->nlon - 1.

• err – [out] Error reported by the function (if any).