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