wigner_functions.h

Go to the documentation of this file.

 
#ifndef wigner_functions_h
#define wigner_functions_h
 
#include "rational.h"
 
#include <wigner/wigxjpf/inc/wigxjpf.h>
 
#include <algorithm>
#include <array>
 
#ifdef FAST_WIGNER_PATH_3J
#define DO_FAST_WIGNER 1
#include <wigner/fastwigxj/inc/fastwigxj.h>
#else
#define DO_FAST_WIGNER 0
#endif
 
 
std::pair<Rational, Rational> wigner_limits(std::pair<Rational, Rational> a, std::pair<Rational, Rational> b);
 
template<Index pos> constexpr
std::pair<Rational, Rational> wigner3j_limits([[maybe_unused]] const Rational a=0,
                                              [[maybe_unused]] const Rational b=0,
                                              [[maybe_unused]] const Rational c=0,
                                              [[maybe_unused]] const Rational d=0,
                                              [[maybe_unused]] const Rational e=0) {
  static_assert(pos < 7 and pos > 0, "Only valid for pos := 1, 2, 3, 4, 5, and 6");
  
  if constexpr (pos == 1 or pos == 2) {
    const Rational maxX = b - a;
    std::pair<Rational, Rational> out {-maxX, a + b};
    if (out.first > out.second) swap(out.first, out.second);
    if (out.second > maxX) out.second = maxX;
    if (out.first > maxX) out = {RATIONAL_UNDEFINED, RATIONAL_UNDEFINED};
    return out;
  } else if constexpr (pos == 3) {
    const Rational maxX = a + b;
    const Rational minX = abs(a - b);
    if (maxX >= minX) return {minX, maxX};
    return {RATIONAL_UNDEFINED, RATIONAL_UNDEFINED};
  } else {
    const Rational lim = pos == 4 ? abs(a) :
                         pos == 5 ? abs(b) :
                       /*pos == 6*/ abs(c);
    const Rational val = - e - d;
    if (-lim <= val and val <= lim) return {val, val};
    return {RATIONAL_UNDEFINED, RATIONAL_UNDEFINED};
  }
}
 
Numeric wigner3j(const Rational j1,
                 const Rational j2,
                 const Rational j3,
                 const Rational m1,
                 const Rational m2,
                 const Rational m3);
 
Numeric wigner6j(const Rational j1,
                 const Rational j2,
                 const Rational j3,
                 const Rational l1,
                 const Rational l2,
                 const Rational l3);
 
Index make_wigner_ready(int largest, int fastest, int size);
 
bool is_wigner_ready(int j);
 
bool is_wigner3_ready(const Rational& J);
 
bool is_wigner6_ready(const Rational& J);
 
template <class ... Integer> constexpr
int temp_init_size(Integer... vals) noexcept {
  constexpr auto N = sizeof...(Integer);
  static_assert(N > 0);
  const std::array<int, N> v{int(vals)...};
  return 1 + 2 * (*std::max_element(v.begin(), v.end()));
}
 
Numeric dwigner3j(Index M, Index J1, Index J2, Index J);
 
Numeric dwigner6j(Index A, Index B, Index C, Index D, Index F);
 
#endif  // wigner_functions_h