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.
-
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(\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
isCHARM_CRD_POINT_GRID
,CHARM_CRD_POINT_GRID_GL
,CHARM_CRD_POINT_GRID_DH1
orCHARM_CRD_POINT_GRID_DH2
(seecharm_point
), orscattered points if
pnt->type
isCHARM_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
), andpnt->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. Below are shown some examples of symmetric and non-symmetric grids. (While these examples use the unit degrees for a more intuitive understanding, note that the function inputs assume radians.)Examples of some symmetric grids (shown are only latitudes
pnt->lat
, the longitudespnt->lon
do not affect the grid symmetry):Equator included:
or-90.0, -60.0, -30.0, 0.0, 30.0, 60.0, 90.0
-80.0, -60.0, -40.0, -20.0, 0.0, 20.0, 40.0, 60.0, 80.0
Equator excluded:
or-35.0, -25.0, -15.0, -5.0, 5.0, 15.0, 25.0, 35.0
-90.0, -85.0, -80.0, 80.0, 85.0, 90.0
Varying spacing, equator excluded:
-90.0, -80.0, -75.0, -70.0, -69.0, 69.0, 70.0, 75.0, 80.0, 90.0
Examples of some non-symmetric grids (shown are only latitudes again):
The negative latitude of
-90.0 deg
does not have its positive counterpart-90.0, -60.0, -30.0, 0.0, 30.0, 60.0
The relative difference between
fabs(-4.99)
andfabs(5.0)
is larger thancharm_glob_threshold2
(see charm_glob):-35.0, -25.0, -15.0, -4.99, 5.0, 15.0, 25.0, 35.0
In other words, the grid is symmetric in the latitudinal direction, provided that all positive latitudes have their negative counterparts at the right places. The check is done within a given threshold to suppress numerical inaccuracies, see the
charm_glob_threshold2
variable in the charm_glob module. The zero latitude, i.e., the equator, may be present in the grid (again, within the specified numerical thresholdcharm_glob_threshold2
).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 ofnmin
(but satisfying1 <= nmin <= nmax
), you can set the respective coefficients inshcs
to zero and use the newly modified structureshcs
as the input to this function. The synthesis will still formally start at degree0
, but since the coefficients up to degreenmin - 1
are zero, they do not contribute to the output signal as required.Note
The function does not use
cell->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 pointerf
must have an access topnt->nlat * pnt->nlon
array elements. The value of the signal synthesized atpnt->lat[i]
andpnt->lon[j]
can be found as:f[i * pnt->nlon + j]
withi = 0, 1, ..., pnt->nlat - 1
andj = 0, 1, ..., pnt->nlon 1
. In casepnt
represents scattered points,f
must have an access topnt->nlat == pnt->nlon
elements andf[i]
stands for the value synthesized atpnt->lat[i]
andpnt->lon[i]
withi = 0, 1, ..., pnt->nlat - 1
.err – [out] Error reported by the function (if any).
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}(\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} \\ &\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) \\ &\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
isCHARM_CRD_CELL_GRID
(seecharm_cell
), orscattered cells if
cell->type
isCHARM_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
), andthe maximum longitude of the
j
th cell must be equal to the minimum longitude of thej + 1
th 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. Below are shown some examples of symmetric and non-symmetric grids. (While these examples use the unit degrees for a more intuitive understanding, note that the function inputs assume radians.)Examples of some symmetric grids (shown are only the minimum and the maximum latitudes
latmin
andlatmax
, respectively; the longitudes do not affect the grid symmetry):latmin: -90.0, -60.0, -30.0, 0.0, 30.0, 60.0 latmax: -60.0, -30.0, 0.0, 30.0, 60.0, 90.0
latmin: -35.0, -25.0, -15.0, -5.0, 5.0, 15.0, 25.0 latmax: -25.0, -15.0, -5.0, 5.0, 15.0, 25.0, 35.0
latmin: -90.0, -85.0, 80.0, 85.0 latmax: -85.0, -80.0, 85.0, 90.0
latmin: -90.0, -80.0, -75.0, -70.0, 69.0, 70.0, 75.0, 80.0 latmax: -80.0, -75.0, -70.0, -69.0, 70.0, 75.0, 80.0, 90.0
Examples of some grids that are not considered as symmetric (shown are again only latitudes):
latmin: -90.0, -60.0, -30.0, 0.0, 30.0 latmax: -60.0, -30.0, 0.0, 30.0, 60.0
latmin: -35.0, -25.0, -15.0, 5.0, 15.0, 25.0 latmax: -25.0, -15.0, -4.99, 15.0, 25.0, 35.0
In other words, the grid is symmetric in the latitudinal direction, provided that all positive latitudes have their negative counterparts at the right places. The check is done within a given threshold to suppress numerical inaccuracies, see the
charm_glob_threshold2
variable in the charm_glob module. The zero latitude, i.e., the equator, may be present in the grid (again, within the specified numerical thresholdcharm_glob_threshold2
).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 ofnmin
, see the tip given incharm_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 pointerf
must have an access tocell->nlat * cell->nlon
array elements. The mean value of the signal synthesized at thei
th cell in the latitudinal direction andj
th cell in the longitudinal direction can be found as:f[i * cell->nlon + j]
withi = 0, 1, ..., cell->nlat - 1
andj = 0, 1, ..., cell->nlon - 1
. In casecell
represents scattered cells,f
must have an access tocell->nlat == cell->nlon
elements andf[i]
stands for the mean value synthesized at thei
th cell withi = 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 byshcs1
. 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 inshcs2
:\[\begin{split}\tilde{f}(\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 toCHARM_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()
andcharm_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 fromshcs2
.Note
The synthesis automatically starts at degree
nmin = 0
. To use some other value ofnmin
, see the tip given incharm_shs_point()
.Warning
Both scaling constants of the
shcs2
structure that defines the irregular surface (shcs2->mu
andshcs2->r
) have to be equal to1.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, fornmax1 = 150
andnmax3 = 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 radiir
, that is, the synthesis ofshcs2
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 andnmax3
defines the maximum harmonic degree of this surface to be used in the synthesis. In theory,nmax3
should be infinite, sinceshcs1->r / r
is non-band-limited forr
band-limited tonmax2
.nmax4 – [in] Maximum harmonic degree, up to which the powers
(shcs1->r / r)^(n + 1)
are harmonically analyzed. The higher is the value ofnmax4
, the smaller is the aliasing effect. In theory,nmax4
should be infinite, sinceshcs1->r / r
is non-band-limited for band-limitedr
. Note thatnmax4
cannot be smaller thannmax3
. In other words, the term(shcs1->r / r)^(n + 1)
is harmonically analyzed up to degreenmax4
, but to synthesize it later, thenmax3
maximum degree is used.f – [out] Pointer to an output array with the synthesized mean values. The pointer
f
must have an access tocell->nlat * cell->nlon
array elements. The mean value of the signal synthesized at thei
th cell in the latitudinal direction andj
th cell in the longitudinal direction can be found as:f[i * cell->nlon + j]
withi = 0, 1, ..., cell->nlat - 1
andj = 0, 1, ..., cell->nlon - 1
.err – [out] Error reported by the function (if any).