SpECTRE
v2024.05.11

Defines the C++ interface to SPHEREPACK. More...
#include <Spherepack.hpp>
Classes  
struct  InterpolationInfo 
Struct to hold cached information at a set of target interpolation points. More...  
Public Types  
using  FirstDeriv = tnsr::i< DataVector, 2, Frame::ElementLogical > 
Type returned by gradient function.  
using  SecondDeriv = tnsr::ij< DataVector, 2, Frame::ElementLogical > 
Type returned by second derivative function.  
Public Member Functions  
Spherepack (size_t l_max, size_t m_max)  
Here l_max and m_max are the largest fullyrepresented l and m in the Ylm expansion.  
std::array< size_t, 2 >  physical_extents () const 
void  gradient (const std::array< double *, 2 > &df, gsl::not_null< const double * > collocation_values, size_t physical_stride=1, size_t physical_offset=0) const 
Computes Pfaffian derivative (df/dtheta, csc(theta) df/dphi) at the collocation values. To act on a slice of the input and output arrays, specify stride and offset (assumed to be the same for input and output).  
void  gradient_from_coefs (const std::array< double *, 2 > &df, gsl::not_null< const double * > spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0, size_t physical_stride=1, size_t physical_offset=0) const 
Same as gradient , but takes the spectral coefficients (rather than collocation values) of the function. This is more efficient if one happens to already have the spectral coefficients. To act on a slice of the input and output arrays, specify strides and offsets.  
void  scalar_laplacian (gsl::not_null< double * > scalar_laplacian, gsl::not_null< const double * > collocation_values, size_t physical_stride=1, size_t physical_offset=0) const 
Computes Laplacian in physical space. To act on a slice of the input and output arrays, specify stride and offset (assumed to be the same for input and output).  
void  scalar_laplacian_from_coefs (gsl::not_null< double * > scalar_laplacian, gsl::not_null< const double * > spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0, size_t physical_stride=1, size_t physical_offset=0) const 
Same as scalar_laplacian above, but the input is the spectral coefficients (rather than collocation values) of the function. This is more efficient if one happens to already have the spectral coefficients. To act on a slice of the input and output arrays, specify strides and offsets.  
void  second_derivative (const std::array< double *, 2 > &df, gsl::not_null< SecondDeriv * > ddf, gsl::not_null< const double * > collocation_values, size_t physical_stride=1, size_t physical_offset=0) const 
Computes Pfaffian first and second derivative in physical space. The first derivative is \(df(i) = d_i f\), and the second derivative is \(ddf(i,j) = d_i (d_j f)\), where \(d_0 = d/d\theta\) and \(d_1 = csc(\theta) d/d\phi\). ddf is not symmetric. To act on a slice of the input and output arrays, specify stride and offset (assumed to be the same for input and output).  
std::pair< FirstDeriv, SecondDeriv >  first_and_second_derivative (const DataVector &collocation_values) const 
Simpler, less general interface to second_derivative.  
double  definite_integral (gsl::not_null< const double * > collocation_values, size_t physical_stride=1, size_t physical_offset=0) const 
Computes the integral over the sphere.  
const std::vector< double > &  integration_weights () const 
Returns weights \(w_i\) such that \(sum_i (c_i w_i)\) is the definite integral, where \(c_i\) are collocation values at point i.  
template<typename T >  
InterpolationInfo< T >  set_up_interpolation_info (const std::array< T, 2 > &target_points) const 
Sets up the InterpolationInfo structure for interpolating onto a set of target \((\theta,\phi)\) points. Does not depend on the function being interpolated.  
template<typename T >  
void  interpolate (gsl::not_null< T * > result, gsl::not_null< const double * > collocation_values, const InterpolationInfo< T > &interpolation_info, size_t physical_stride=1, size_t physical_offset=0) const 
Interpolates from collocation_values onto the points that have been passed into the set_up_interpolation_info function. To interpolate a different function on the same spectral grid, there is no need to recompute interpolation_info . If you specify stride and offset, acts on a slice of the input values. The output has unit stride.  
template<typename T , typename R >  
void  interpolate_from_coefs (gsl::not_null< T * > result, const R &spectral_coefs, const InterpolationInfo< T > &interpolation_info, size_t spectral_stride=1, size_t spectral_offset=0) const 
Same as interpolate , but assumes you have spectral coefficients. This is more efficient if you already have the spectral coefficients available. If you specify stride and offset, acts on a slice of the input coefs. The output has unit stride.  
template<typename T >  
T  interpolate (const DataVector &collocation_values, const std::array< T, 2 > &target_points) const 
Simpler interface to interpolate . If you need to call this repeatedly on different spectral_coefs or collocation_values for the same target points, this is inefficient; instead use set_up_interpolation_info and the functions that use InterpolationInfo .  
template<typename T >  
T  interpolate_from_coefs (const DataVector &spectral_coefs, const std::array< T, 2 > &target_points) const 
DataVector  prolong_or_restrict (const DataVector &spectral_coefs, const Spherepack &target) const 
Takes spectral coefficients compatible with *this , and either prolongs them or restricts them to be compatible with target . This is done by truncation (restriction) or padding with zeros (prolongation).  
size_t  l_max () const 
Sizes in physical and spectral space for this instance.  
size_t  m_max () const 
Sizes in physical and spectral space for this instance.  
size_t  physical_size () const 
Sizes in physical and spectral space for this instance.  
size_t  spectral_size () const 
Sizes in physical and spectral space for this instance.  
const std::vector< double > &  theta_points () const 
Collocation points theta and phi. More...  
const std::vector< double > &  phi_points () const 
Collocation points theta and phi. More...  
std::array< DataVector, 2 >  theta_phi_points () const 
Collocation points theta and phi. More...  
void  phys_to_spec (gsl::not_null< double * > spectral_coefs, gsl::not_null< const double * > collocation_values, size_t physical_stride=1, size_t physical_offset=0, size_t spectral_stride=1, size_t spectral_offset=0) const 
Spectral transformations. To act on a slice of the input and output arrays, specify strides and offsets.  
void  spec_to_phys (gsl::not_null< double * > collocation_values, gsl::not_null< const double * > spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0, size_t physical_stride=1, size_t physical_offset=0) const 
Spectral transformations. To act on a slice of the input and output arrays, specify strides and offsets.  
void  phys_to_spec_all_offsets (gsl::not_null< double * > spectral_coefs, gsl::not_null< const double * > collocation_values, size_t stride) const 
Spectral transformations where collocation_values and spectral_coefs are assumed to point to 3dimensional arrays (I1 x S2 topology), and the transformations are done for all 'radial' points at once by internally looping over all values of the offset from zero to stride 1 (the physical and spectral strides are equal and are called stride ).  
void  spec_to_phys_all_offsets (gsl::not_null< double * > collocation_values, gsl::not_null< const double * > spectral_coefs, size_t stride) const 
Spectral transformations where collocation_values and spectral_coefs are assumed to point to 3dimensional arrays (I1 x S2 topology), and the transformations are done for all 'radial' points at once by internally looping over all values of the offset from zero to stride 1 (the physical and spectral strides are equal and are called stride ).  
DataVector  phys_to_spec (const DataVector &collocation_values, size_t physical_stride=1, size_t physical_offset=0) const 
Simpler, less general interfaces to phys_to_spec and spec_to_phys . Acts on a slice of the input and returns a unitstride result.  
DataVector  spec_to_phys (const DataVector &spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0) const 
Simpler, less general interfaces to phys_to_spec and spec_to_phys . Acts on a slice of the input and returns a unitstride result.  
DataVector  phys_to_spec_all_offsets (const DataVector &collocation_values, size_t stride) const 
Simpler, less general interfaces to phys_to_spec_all_offsets and spec_to_phys_all_offsets . Result has the same stride as the input.  
DataVector  spec_to_phys_all_offsets (const DataVector &spectral_coefs, size_t stride) const 
Simpler, less general interfaces to phys_to_spec_all_offsets and spec_to_phys_all_offsets . Result has the same stride as the input.  
void  gradient_all_offsets (const std::array< double *, 2 > &df, gsl::not_null< const double * > collocation_values, size_t stride=1) const 
Same as gradient but pointers are assumed to point to 3dimensional arrays (I1 x S2 topology), and the gradient is done for all 'radial' points at once by internally looping over all values of the offset from zero to stride 1.  
void  gradient_from_coefs_all_offsets (const std::array< double *, 2 > &df, gsl::not_null< const double * > spectral_coefs, size_t stride=1) const 
Same as gradient but pointers are assumed to point to 3dimensional arrays (I1 x S2 topology), and the gradient is done for all 'radial' points at once by internally looping over all values of the offset from zero to stride 1.  
FirstDeriv  gradient (const DataVector &collocation_values, size_t physical_stride=1, size_t physical_offset=0) const 
Simpler, less general interfaces to gradient . Acts on a slice of the input and returns a unitstride result.  
FirstDeriv  gradient_from_coefs (const DataVector &spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0) const 
Simpler, less general interfaces to gradient . Acts on a slice of the input and returns a unitstride result.  
FirstDeriv  gradient_all_offsets (const DataVector &collocation_values, size_t stride=1) const 
Simpler, less general interfaces to gradient_all_offsets . Result has the same stride as the input.  
FirstDeriv  gradient_from_coefs_all_offsets (const DataVector &spectral_coefs, size_t stride=1) const 
Simpler, less general interfaces to gradient_all_offsets . Result has the same stride as the input.  
DataVector  scalar_laplacian (const DataVector &collocation_values, size_t physical_stride=1, size_t physical_offset=0) const 
Simpler, less general interfaces to scalar_laplacian . Acts on a slice of the input and returns a unitstride result.  
DataVector  scalar_laplacian_from_coefs (const DataVector &spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0) const 
Simpler, less general interfaces to scalar_laplacian . Acts on a slice of the input and returns a unitstride result.  
Static Public Member Functions  
static void  add_constant (const gsl::not_null< DataVector * > spectral_coefs, const double c) 
Adds a constant (i.e. \(f(\theta,\phi)\) += \(c\)) to the function given by the spectral coefficients, by modifying the coefficients.  
static double  average (const DataVector &spectral_coefs) 
Returns the average of \(f(\theta,\phi)\) over \((\theta,\phi)\).  
static constexpr size_t  physical_size (const size_t l_max, const size_t m_max) 
Static functions to return the correct sizes of vectors of collocation points and spectral coefficients for a given l_max and m_max. Useful for allocating space without having to create a Spherepack.  
static constexpr size_t  spectral_size (const size_t l_max, const size_t m_max) 
Defines the C++ interface to SPHEREPACK.
The class Spherepack
defines the C++ interface to the fortran library SPHEREPACK used for computations on the surface of a sphere.
Given a realvalued, scalar function \(g(\theta, \phi)\), SPHEREPACK expands it as:
\begin{align} g(\theta, \phi) &=\frac{1}{2}\sum_{l=0}^{l_{\max}}\bar P_l^0(\cos\theta) a_{l0} +\sum_{l=1}^{l_{\max}}\sum_{m=1}^{\min(l, m_{\max})}\bar P_l^m(\cos\theta)\{ a_{lm}\cos m\phi b_{lm}\sin m\phi\}\label{eq:spherepack_expansion} \end{align}
where \(a_{lm}\) and \(b_{lm}\) are realvalued spectral coefficient arrays used by SPHEREPACK, \(P_l^m(x)\) are defined as
\begin{align} \bar P_l^m(x)&=\sqrt{\frac{(2l+1)(lm)!}{2(l+m)!}}\;P_{lm}(x) \end{align}
and \(P_{nm}(x)\) are the associated Legendre polynomials as defined, for example, in Jackson's "Classical Electrodynamics".
The standard expansion of \(g(\theta, \phi)\) in terms of scalar spherical harmonics is
\begin{align} g(\theta, \phi) &= \sum_{l=0}^{l_{\max}}\sum_{m=\min(l, m_{\max})}^{\min(l, m_{\max})} A_{lm} Y_{lm}(\theta,\phi), \end{align}
where \(Y_{lm}(\theta,\phi)\) are the usual complexvalued scalar spherical harmonics (as defined, for example, in Jackson's "Classical Electrodynamics") and \(A_{lm}\) are complex coefficients.
The relationship between the complex coefficients \(A_{lm}\) and SPHEREPACK's realvalued \(a_{lm}\) and \(b_{lm}\) is
\begin{align} a_{l0} & = \sqrt{\frac{2}{\pi}}A_{l0}&\qquad l\geq 0,\\ a_{lm} & = (1)^m\sqrt{\frac{2}{\pi}} \mathrm{Re}(A_{lm}) &\qquad l\geq 1, m\geq 1, \\ b_{lm} & = (1)^m\sqrt{\frac{2}{\pi}} \mathrm{Im}(A_{lm}) &\qquad l\geq 1, m\geq 1. \end{align}
Sometimes it is useful to expand a realvalued function in the form
\begin{align} g(\theta, \phi) &= \sum_{l=0}^\infty\sum_{m=0}^l \left[ c_{lm}\mathrm{Re}(Y_{lm}(\theta, \phi))+ d_{nm}\mathrm{Im}(Y_{lm}(\theta, \phi)) \right]. \end{align}
The coefficients here are therefore
\begin{align} c_{l0} &= A_{l0},\\ c_{lm} &= 2\mathrm{Re}(A_{lm}) \qquad m\geq 1,\\ d_{lm} &=2\mathrm{Im}(A_{lm}). \end{align}
Internally, SPHEREPACK can represent its expansion in two ways which we will refer to as modal and nodal representations:
spectral_coefs
in the methods below. For this C++ interface, they are saved in a single DataVector
. To help you index the coefficients as expected by this interface, use the class SpherepackIterator
.collocation_values
in the methods below. This is an array of the expanded function \(g(\theta,\phi)\) evaluated at collocation values \((\theta_i,\phi_j)\), where \(\theta_i\) are GaussLegendre quadrature nodes in the interval \((0, \pi)\) with \(i = 0, ..., l_{\max}\), and \(\phi_j\) is distributed uniformly in \((0, 2\pi)\) with \(i = 0, ..., 2m_{\max}\). The angles of the collocation points can be computed with the method theta_phi_points
.To convert between the two representations the methods spec_to_phys
and phys_to_spec
can be used. For internal calculations SPHEREPACK will usually convert to spectral coefficients first, so it is in general more efficient to use these directly.
Most methods of SPHEREPACK will compute the requested values of e.g. gradient
or scalar_laplacian
at the collocation points, effectively returning an expansion in nodal form as defined above. To evaluate the function at arbitrary angles \(\theta\), \(\phi\), these values have to be "interpolated" (i.e. the new expansion evaluated) using interpolate
.
Spherepack stores two types of quantities:

inline 
Collocation points theta and phi.
The phi points are uniform in phi, with the first point at phi=0.
The theta points are GaussLegendre in \(\cos(\theta)\), so there are no points at the poles.

inlinestaticconstexpr 
spectral_size
is the size of the buffer that holds the coefficients; it is not the number of coefficients (which is \(m_{\rm max}^2+(l_{\rm max}m_{\rm max})(2m_{\rm max}+1)\)). To simplify its internal indexing, SPHEREPACK uses a buffer with more space than necessary. See SpherepackIterator for how to index the coefficients in the buffer. std::array< DataVector, 2 > ylm::Spherepack::theta_phi_points  (  )  const 
Collocation points theta and phi.
The phi points are uniform in phi, with the first point at phi=0.
The theta points are GaussLegendre in \(\cos(\theta)\), so there are no points at the poles.

inline 
Collocation points theta and phi.
The phi points are uniform in phi, with the first point at phi=0.
The theta points are GaussLegendre in \(\cos(\theta)\), so there are no points at the poles.