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              : };