Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <vector>
       8             : 
       9             : /*!
      10             :  * \brief Computes Wigner 3J symbols.
      11             :  *
      12             :  * Given $l_2$, $m_2$, $l_3$, and $m_3$,
      13             :  * computes and caches the Wigner 3J symbols
      14             :  *
      15             :  * \begin{align}
      16             :  * \left(
      17             :  * \begin{array}{ccc}
      18             :  *  l_1 & l_2 & l_3 \\
      19             :  * -m_3-m_2 & m_2 & m_3
      20             :  * \end{array}\right)
      21             :  * \end{align}
      22             :  *
      23             :  * for all $l_1$.
      24             :  *
      25             :  * Internally, uses recursion relations to efficiently compute many 3J symbols
      26             :  * at once.  This should be more efficient than having a single function
      27             :  * that takes $l_2$, $m_2$, $l_3$, $m_3$, $l_1$
      28             :  * and computes a single 3J symbol, assuming
      29             :  * that the user creates WignerThreeJ objects infrequently, and that the user
      30             :  * calls each WignerThreeJ object multiple times.
      31             :  *
      32             :  * Here we limit the quantum numbers to integer values.  It is trivial to
      33             :  * extend to half-integer values because the underlying routine (the
      34             :  * drc3jj.f function from the slatec library) supports it.
      35             :  *
      36             :  * Uses lazy evaluation: the constructor computes nothing, but the first time
      37             :  * operator() is called for a nonzero 3J symbol, it computes all the
      38             :  * nonzero 3J symbols.
      39             :  *
      40             :  * Internally, allocates a std::vector.  A future optimization may be to
      41             :  * pass in a buffer of the appropriate size, so that WignerThreeJ would not
      42             :  * need to do memory allocations.
      43             :  */
      44           1 : class WignerThreeJ {
      45             :  public:
      46           0 :   WignerThreeJ(size_t l2, int m2, size_t l3, int m3);
      47             :   /*!
      48             :    * Returns the Wigner 3J symbol
      49             :    * \begin{align}
      50             :    * \left(
      51             :    * \begin{array}{ccc}
      52             :    *  l_1 & l_2 & l_3 \\
      53             :    * -m_3-m_2 & m_2 & m_3
      54             :    * \end{array}\right).
      55             :    * \end{align}
      56             :    *
      57             :    * Internally, computes and caches all nonzero 3J symbols for all $l_1$.
      58             :    * Returns zero for $l_1 < l_1^\mathrm{min}$ and
      59             :    * $l_1 > l_1^\mathrm{max}$.
      60             :    */
      61           1 :   double operator()(size_t l1);
      62             :   /// Returns the smallest value of $l_1$ that gives a nonzero 3J symbol.
      63           1 :   size_t l1_min() const { return l1_min_; }
      64             :   /// Returns the largest value of $l_1$ that gives a nonzero 3J symbol.
      65           1 :   size_t l1_max() const { return l1_max_; }
      66             : 
      67             :  private:
      68           0 :   void recompute();
      69             :   // The underlying routine takes doubles (to deal with half-integer spins),
      70             :   // so we need to store doubles here.
      71           0 :   double l2_, m2_, l3_, m3_;
      72           0 :   size_t l1_min_, l1_max_;
      73           0 :   bool up_to_date_;
      74           0 :   std::vector<double> coefs_;
      75             : };