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