Line data Source code
1 1 : // Distributed under the MIT License.
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 <optional>
13 : #include <ostream>
14 : #include <vector>
15 :
16 : #include "Utilities/ErrorHandling/Assert.hpp"
17 :
18 : namespace ylm {
19 :
20 : /*!
21 : * \ingroup SpectralGroup
22 : * \brief
23 : * Iterates over spectral coefficients stored in SPHEREPACK format.
24 : *
25 : * \details
26 : * The internal SPHEREPACK ordering is not intuitive, so
27 : * SpherepackIterator exists for the purpose of iterating over an
28 : * array containing SPHEREPACK coefficients and determining the
29 : * (l,m) of each entry in the array.
30 : *
31 : * SPHEREPACK expands \f$f(\theta,\phi)\f$ as
32 : * \f[
33 : * f(\theta,\phi) =
34 : * \frac{1}{2} \sum_{l=0}^{l_{max}} \bar{P}_l^0 a(0,l)
35 : * + \sum_{m=1}^{m_{max}} \sum_{l=m}^{l_{max}} \bar{P}_l^m
36 : * \left( a(m,l) \cos(m \phi) - b(m,l) \sin(m \phi)\right)
37 : * \f]
38 : * where \f$a(m,l)\f$ and \f$b(m,l)\f$ are the SPHEREPACK
39 : * spectral coefficients, and \f$\bar{P}_l^m\f$ are
40 : * unit-orthonormal associated Legendre polynomials:
41 : * \f[
42 : * \bar{P}_l^m =
43 : * (-1)^m \sqrt{\frac{(2l+1)(l-m)!}{2(l+m)!}} P_l^m,
44 : * \f]
45 : * where \f$P_l^m\f$ are the associated Legendre polynomials as defined
46 : * for example in Jackson "Classical Electrodynamics".
47 : * \example
48 : * \snippet Test_SpherepackIterator.cpp spherepack_iterator_example
49 : */
50 1 : class SpherepackIterator {
51 : public:
52 : /// SPHEREPACK has two coefficient variables, 'a' and 'b', that hold
53 : /// the cos(m*phi) and sin(m*phi) parts of the spectral
54 : /// coefficients.
55 0 : enum class CoefficientArray { a, b };
56 :
57 : /// Construct a SpherepackIterator.
58 : ///
59 : /// For all uses of SpherepackIterator by the SPHEREPACK routines,
60 : /// zero_m_is_real should be true, indicating that the m=0 values of
61 : /// the coefficients are real (or in other words, the 'b' array has
62 : /// no m=0 values).
63 : ///
64 : /// However, we sometimes use the same storage scheme (and hence
65 : /// SpherepackIterator) to hold spin-weighted-spherical-harmonic
66 : /// coefficients of complex quantities, in which case the m=0 values
67 : /// of the coefficients can be complex; in that case, zero_m_is_real
68 : /// should be false.
69 1 : SpherepackIterator(size_t l_max_input, size_t m_max_input, size_t stride = 1,
70 : bool zero_m_is_real = true);
71 :
72 0 : size_t l_max() const { return l_max_; }
73 0 : size_t m_max() const { return m_max_; }
74 0 : size_t n_th() const { return n_th_; }
75 0 : size_t n_ph() const { return n_ph_; }
76 0 : size_t stride() const { return stride_; }
77 0 : bool zero_m_is_real() const { return zero_m_is_real_; }
78 :
79 : /// Size of a SPHEREPACK coefficient array (a and b combined), not
80 : /// counting stride. For non-unit stride, the size of the array
81 : /// should be spherepack_array_size()*stride.
82 1 : size_t spherepack_array_size() const {
83 : return (l_max_ + 1) * (m_max_ + 1) * 2;
84 : }
85 :
86 0 : SpherepackIterator& operator++() {
87 : ASSERT(current_compact_index_ != offset_into_spherepack_array.size(),
88 : "Incrementing an invalid iterator: "
89 : << current_compact_index_ << " "
90 : << offset_into_spherepack_array.size());
91 : ++current_compact_index_;
92 : return *this;
93 : }
94 :
95 0 : explicit operator bool() const {
96 : return current_compact_index_ < offset_into_spherepack_array.size();
97 : }
98 :
99 : /// Current index into a SPHEREPACK coefficient array.
100 1 : size_t operator()() const {
101 : return offset_into_spherepack_array[current_compact_index_];
102 : }
103 :
104 : /// Given an offset into a SPHEREPACK coefficient array, return the compact
105 : /// index corresponding to that offset.
106 : ///
107 : /// Essentially the inverse of operator(). If the offset points to an element
108 : /// that SPHEREPACK doesn't actually use (i.e. no compact index can reach the
109 : /// given offset), then a std::nullopt is returned.
110 1 : std::optional<size_t> compact_index(const size_t offset) const {
111 : return offset_to_compact_index_[offset];
112 : }
113 :
114 : /// Returns the current compact index that SpherepackIterator uses internally.
115 : /// This does not index a SPHEREPACK coefficient array. \see operator() for an
116 : /// index into a SPHEREPACK coefficient array
117 1 : size_t current_compact_index() const { return current_compact_index_; }
118 :
119 : /// Current values of l and m.
120 1 : size_t l() const { return compact_l_[current_compact_index_]; }
121 0 : size_t m() const { return compact_m_[current_compact_index_]; }
122 : /// Whether the iterator points to an element of 'a' or 'b',
123 : /// i.e. points to the cos(m*phi) or sin(m*phi) part of a spectral
124 : /// coefficient.
125 1 : CoefficientArray coefficient_array() const {
126 : return (current_compact_index_ < number_of_valid_entries_in_a_
127 : ? CoefficientArray::a
128 : : CoefficientArray::b);
129 : }
130 :
131 : /// Reset iterator back to beginning value. Returns *this.
132 1 : SpherepackIterator& reset();
133 :
134 : /// Set iterator to specific value of l, m, array. Returns *this.
135 1 : SpherepackIterator& set(size_t l_input, size_t m_input,
136 : CoefficientArray coefficient_array_input);
137 :
138 : /// Same as 'set' above, but assumes CoefficientArray is 'a' for m>=0 and
139 : /// 'b' for m<0. This is useful when converting between true
140 : /// spherical harmonics (which allow negative values of m) and
141 : /// SPHEREPACK coefficients (which have only positive values of m,
142 : /// but two arrays for sin(m*phi) and cos(m*phi) parts).
143 1 : SpherepackIterator& set(size_t l_input, int m_input) {
144 : return set(l_input, size_t(abs(m_input)),
145 : m_input >= 0 ? CoefficientArray::a : CoefficientArray::b);
146 : }
147 :
148 : /// Set iterator to a specific compact index. Returns *this.
149 1 : SpherepackIterator& set(size_t compact_index);
150 :
151 : private:
152 0 : size_t l_max_, m_max_, n_th_, n_ph_, stride_;
153 0 : bool zero_m_is_real_;
154 0 : size_t number_of_valid_entries_in_a_;
155 0 : size_t current_compact_index_;
156 0 : std::vector<size_t> offset_into_spherepack_array;
157 0 : std::vector<size_t> compact_l_, compact_m_;
158 0 : std::vector<std::optional<size_t>> offset_to_compact_index_;
159 : };
160 :
161 0 : inline bool operator==(const SpherepackIterator& lhs,
162 : const SpherepackIterator& rhs) {
163 : return lhs.l_max() == rhs.l_max() and lhs.m_max() == rhs.m_max() and
164 : lhs.stride() == rhs.stride() and lhs() == rhs() and
165 : lhs.zero_m_is_real() == rhs.zero_m_is_real();
166 : }
167 :
168 0 : inline bool operator!=(const SpherepackIterator& lhs,
169 : const SpherepackIterator& rhs) {
170 : return not(lhs == rhs);
171 : }
172 :
173 0 : inline std::ostream& operator<<(
174 : std::ostream& os,
175 : const SpherepackIterator::CoefficientArray& coefficient_array) {
176 : return os << (coefficient_array == SpherepackIterator::CoefficientArray::a
177 : ? 'a'
178 : : 'b');
179 : }
180 : } // namespace ylm
|
