CoulombIntegrals.hpp

#pragma once
#include "Angular/Angular.hpp"
#include "Wavefunction/DiracSpinor.hpp"
#include <optional>
#include <vector>
 
namespace Coulomb {
 
std::vector<double> yk_ab(const int k, const DiracSpinor &Fa,
                          const DiracSpinor &Fb, const std::size_t maxi = 0);
 
void yk_ab(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
           std::vector<double> &ykab, const std::size_t maxi = 0);
 
void bk_ab(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
           std::vector<double> &b0, std::vector<double> &binf,
           const std::size_t maxi = 0);
 
void gk_ab(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
           std::vector<double> &g0, std::vector<double> &ginf,
           const std::size_t maxi = 0);
 
//==============================================================================
 
double Rk_abcd(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
               const DiracSpinor &Fc, const DiracSpinor &Fd);
 
double Rk_abcd(const DiracSpinor &Fa, const DiracSpinor &Fc,
               const std::vector<double> &ykbd);
 
DiracSpinor Rkv_bcd(const int k, const int kappa_v, const DiracSpinor &Fb,
                    const DiracSpinor &Fc, const DiracSpinor &Fd);
 
DiracSpinor Rkv_bcd(const int kappa_v, const DiracSpinor &Fc,
                    const std::vector<double> &ykbd);
 
void Rkv_bcd(DiracSpinor *const Rkv, const DiracSpinor &Fc,
             const std::vector<double> &ykbd);
 
//==============================================================================
 
double Qk_abcd(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
               const DiracSpinor &Fc, const DiracSpinor &Fd);
 
DiracSpinor Qkv_bcd(const int k, int kappa_v, const DiracSpinor &Fb,
                    const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//==============================================================================
double g_abcd(const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c,
              const DiracSpinor &d, int tma, int tmb, int tmc, int tmd);
 
//==============================================================================
bool Qk_abcd_SR(int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
                const DiracSpinor &Fc, const DiracSpinor &Fd);
bool Qk_abcd_SR(int k, int ka, int kb, int kc, int kd);
 
bool Pk_abcd_SR(int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
                const DiracSpinor &Fc, const DiracSpinor &Fd);
 
bool Pk_abcd_SR(int k, int ka, int kb, int kc, int kd);
 
//==============================================================================
 
double Pk_abcd(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
               const DiracSpinor &Fc, const DiracSpinor &Fd);
 
DiracSpinor Pkv_bcd(const int k, int kappa_v, const DiracSpinor &Fb,
                    const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//==============================================================================
 
 
double Wk_abcd(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb,
               const DiracSpinor &Fc, const DiracSpinor &Fd);
 
DiracSpinor Wkv_bcd(const int k, int kappa_v, const DiracSpinor &Fb,
                    const DiracSpinor &Fc, const DiracSpinor &Fd);
 
//==============================================================================
 
std::pair<int, int> k_minmax_Ck(const DiracSpinor &a, const DiracSpinor &b);
std::pair<int, int> k_minmax_Ck(int kappa_a, int kappa_b);
 
std::pair<int, int> k_minmax_tj(int tja, int tjb);
 
std::pair<int, int> k_minmax_Q(const DiracSpinor &a, const DiracSpinor &b,
                               const DiracSpinor &c, const DiracSpinor &d);
 
std::pair<int, int> k_minmax_Q(int kappa_a, int kappa_b, int kappa_c,
                               int kappa_d);
 
std::pair<int, int> k_minmax_P(const DiracSpinor &a, const DiracSpinor &b,
                               const DiracSpinor &c, const DiracSpinor &d);
 
std::pair<int, int> k_minmax_W(const DiracSpinor &a, const DiracSpinor &b,
                               const DiracSpinor &c, const DiracSpinor &d);
 
//==============================================================================
 
template <class A>
static int twojk(const A &a) {
  if constexpr (std::is_same_v<A, DiracSpinor>) {
    return a.twoj();
  } else {
    static_assert(std::is_same_v<A, int>);
    return 2 * a;
  }
}
 
template <class A>
static std::optional<int> twojknull(const A &a) {
  if constexpr (std::is_same_v<A, DiracSpinor>) {
    return a.twoj();
  } else if constexpr (std::is_same_v<A, int>) {
    static_assert(std::is_same_v<A, int>);
    return 2 * a;
  } else {
    return std::nullopt;
  }
}
 
template <class A, class B, class C, class D, class E, class F>
static double sixj(const A &a, const B &b, const C &c, const D &d, const E &e,
                   const F &f) {
  return Angular::sixj_2(twojk(a), twojk(b), twojk(c), twojk(d), twojk(e),
                         twojk(f));
}
 
template <class A = std::optional<int>, class B = std::optional<int>,
          class C = std::optional<int>, class D = std::optional<int>,
          class E = std::optional<int>, class F = std::optional<int>>
static bool sixjTriads(const A &a, const B &b, const C &c, const D &d,
                       const E &e, const F &f) {
  return Angular::sixjTriads(twojknull(a), twojknull(b), twojknull(c),
                             twojknull(d), twojknull(e), twojknull(f));
}
 
template <class A, class B, class C>
static bool triangle(const A &a, const B &b, const C &c) {
  return Angular::triangle(twojk(a), twojk(b), twojk(c));
}
 
} // namespace Coulomb