SpherepackIterator.hpp
Go to the documentation of this file.
2 // See LICENSE.txt for details.
3
4 ///
5 /// \file
6 /// Defines class SpherepackIterator.
7
8 #pragma once
9
10 #include <cmath>
11 #include <cstddef>
12 #include <ostream>
13 #include <vector>
14
15 #include "ErrorHandling/Assert.hpp"
16
17 /*!
18  * \ingroup SpectralGroup
19  * \brief
20  * Iterates over spectral coefficients stored in SPHEREPACK format.
21  *
22  * \details
23  * The internal SPHEREPACK ordering is not intuitive, so
24  * SpherepackIterator exists for the purpose of iterating over an
25  * array containing SPHEREPACK coefficients and determining the
26  * (l,m) of each entry in the array.
27  *
28  * SPHEREPACK expands \f$f(\theta,\phi)\f$ as
29  * \f[
30  * f(\theta,\phi) =
31  * \frac{1}{2} \sum_{l=0}^{l_{max}} \bar{P}_l^0 a(0,l)
32  * + \sum_{m=1}^{m_{max}} \sum_{l=m}^{l_{max}} \bar{P}_l^m
33  * \left( a(m,l) \cos(m \phi) - b(m,l) \sin(m \phi)\right)
34  * \f]
35  * where \f$a(m,l)\f$ and \f$b(m,l)\f$ are the SPHEREPACK
36  * spectral coefficients, and \f$\bar{P}_l^m\f$ are
37  * unit-orthonormal associated Legendre polynomials:
38  * \f[
39  * \bar{P}_l^m =
40  * (-1)^m \sqrt{\frac{(2l+1)(l-m)!}{2(l+m)!}} P_l^m,
41  * \f]
42  * where \f$P_l^m\f$ are the associated Legendre polynomials as defined
43  * for example in Jackson "Classical Electrodynamics".
44  * \example
45  * \snippet Test_SpherepackIterator.cpp spherepack_iterator_example
46  */
48  public:
49  /// SPHEREPACK has two coefficient variables, 'a' and 'b', that hold
50  /// the cos(m*phi) and sin(m*phi) parts of the spectral
51  /// coefficients.
52  enum class CoefficientArray { a, b };
53
54  SpherepackIterator(size_t l_max_input, size_t m_max_input, size_t stride = 1);
55
56  size_t l_max() const noexcept { return l_max_; }
57  size_t m_max() const noexcept { return m_max_; }
58  size_t n_th() const noexcept { return n_th_; }
59  size_t n_ph() const noexcept { return n_ph_; }
60  size_t stride() const noexcept { return stride_; }
61
62  /// Size of a SPHEREPACK coefficient array (a and b combined), not
63  /// counting stride. For non-unit stride, the size of the array
64  /// should be spherepack_array_size()*stride.
65  size_t spherepack_array_size() const noexcept {
66  return (l_max_ + 1) * (m_max_ + 1) * 2;
67  }
68
69  SpherepackIterator& operator++() noexcept {
70  ASSERT(current_compact_index_ != offset_into_spherepack_array.size(),
71  "Incrementing an invalid iterator: "
72  << current_compact_index_ << " "
73  << offset_into_spherepack_array.size());
74  ++current_compact_index_;
75  return *this;
76  }
77
78  explicit operator bool() const noexcept {
79  return current_compact_index_ < offset_into_spherepack_array.size();
80  }
81
82  /// Current index into a SPHEREPACK coefficient array.
83  size_t operator()() const noexcept {
84  return offset_into_spherepack_array[current_compact_index_];
85  }
86
87  /// Current values of l and m.
88  size_t l() const noexcept { return compact_l_[current_compact_index_]; }
89  size_t m() const noexcept { return compact_m_[current_compact_index_]; }
90  /// Whether the iterator points to an element of 'a' or 'b',
91  /// i.e. points to the cos(m*phi) or sin(m*phi) part of a spectral
92  /// coefficient.
94  return (current_compact_index_ < number_of_valid_entries_in_a_
95  ? CoefficientArray::a
96  : CoefficientArray::b);
97  }
98
99  /// Reset iterator back to beginning value. Returns *this.
100  SpherepackIterator& reset() noexcept;
101
102  /// Set iterator to specific value of l, m, array. Returns *this.
103  SpherepackIterator& set(size_t l_input, size_t m_input,
104  CoefficientArray coefficient_array_input) noexcept;
105
106  /// Same as 'set' above, but assumes CoefficientArray is 'a' for m>=0 and
107  /// 'b' for m<0. This is useful when converting between true
108  /// spherical harmonics (which allow negative values of m) and
109  /// SPHEREPACK coefficients (which have only positive values of m,
110  /// but two arrays for sin(m*phi) and cos(m*phi) parts).
111  SpherepackIterator& set(size_t l_input, int m_input) noexcept {
112  return set(l_input, size_t(abs(m_input)),
113  m_input >= 0 ? CoefficientArray::a : CoefficientArray::b);
114  }
115
116  private:
117  size_t l_max_, m_max_, n_th_, n_ph_, stride_;
118  size_t number_of_valid_entries_in_a_;
119  size_t current_compact_index_;
120  std::vector<size_t> offset_into_spherepack_array;
121  std::vector<size_t> compact_l_, compact_m_;
122 };
123
124 inline bool operator==(const SpherepackIterator& lhs,
125  const SpherepackIterator& rhs) {
126  return lhs.l_max() == rhs.l_max() and lhs.m_max() == rhs.m_max() and
127  lhs.stride() == rhs.stride() and lhs() == rhs();
128 }
129
130 inline bool operator!=(const SpherepackIterator& lhs,
131  const SpherepackIterator& rhs) {
132  return not(lhs == rhs);
133 }
134
135 inline std::ostream& operator<<(
136  std::ostream& os,
138  return os << (coefficient_array == SpherepackIterator::CoefficientArray::a
139  ? 'a'
140  : 'b');
141 }
SpherepackIterator & reset() noexcept
Reset iterator back to beginning value. Returns *this.
Definition: SpherepackIterator.cpp:129
CoefficientArray coefficient_array() const noexcept
Whether the iterator points to an element of &#39;a&#39; or &#39;b&#39;, i.e. points to the cos(m*phi) or sin(m*phi) ...
Definition: SpherepackIterator.hpp:93
size_t operator()() const noexcept
Current index into a SPHEREPACK coefficient array.
Definition: SpherepackIterator.hpp:83
size_t spherepack_array_size() const noexcept
Size of a SPHEREPACK coefficient array (a and b combined), not counting stride. For non-unit stride...
Definition: SpherepackIterator.hpp:65
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
size_t l() const noexcept
Current values of l and m.
Definition: SpherepackIterator.hpp:88
CoefficientArray
SPHEREPACK has two coefficient variables, &#39;a&#39; and &#39;b&#39;, that hold the cos(m*phi) and sin(m*phi) parts ...
Definition: SpherepackIterator.hpp:52
Defines macro ASSERT.
Iterates over spectral coefficients stored in SPHEREPACK format.
Definition: SpherepackIterator.hpp:47