StructureRad.hpp

#pragma once
#include "Coulomb/Coulomb.hpp"
#include "Coulomb/meTable.hpp"
#include "IO/FRW_fileReadWrite.hpp"
#include "Wavefunction/DiracSpinor.hpp"
#include <memory>
#include <vector>
class Wavefunction;
class DiracSpinor;
namespace DiracOperator {
class TensorOperator;
}
namespace ExternalField {
class CorePolarisation;
}
 
namespace MBPT {
 
class StructureRad {
 
public:
  StructureRad(const std::vector<DiracSpinor> &basis, double en_core,
               std::pair<int, int> nminmax = {0, 999},
               const std::string &Qk_fname = "",
               const std::vector<double> &fk = {},
               const std::vector<double> &etak = {});
 
private:
  bool m_use_Qk;
  std::vector<double> m_root_fk; // sqrt - add to both sides!
  std::vector<double> m_etak;
  Coulomb::YkTable mY{};
  std::optional<Coulomb::QkTable> mQ{std::nullopt}; //?
  // nb: it seems conter-intuative, but this copy makes it FASTER!
  std::vector<DiracSpinor> mCore{}, mExcited{};
 
public:
  const std::vector<DiracSpinor> &core() const { return mCore; }
  const std::vector<DiracSpinor> &excited() const { return mExcited; }
 
  const Coulomb::YkTable &Yk() const { return mY; }
  const std::optional<Coulomb::QkTable> &Qk() const { return mQ; }
 
  std::pair<double, double>
  srTB(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
       const DiracSpinor &v, double omega = 0.0,
       const ExternalField::CorePolarisation *const dV = nullptr) const;
 
  std::pair<double, double>
  srC(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
      const DiracSpinor &v,
      const ExternalField::CorePolarisation *const dV = nullptr) const;
 
  std::pair<double, double>
  norm(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
       const DiracSpinor &v,
       const ExternalField::CorePolarisation *const dV = nullptr) const;
 
  std::pair<double, double>
  srn(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
      const DiracSpinor &v, double omega = 0.0,
      const ExternalField::CorePolarisation *const dV = nullptr) const {
    const auto [tb, dvtb] = srTB(h, w, v, omega, dV);
    const auto [c, dvc] = srC(h, w, v, dV);
    const auto [n, dvn] = norm(h, w, v, dV);
    return {tb + c + n, dvtb + dvc + dvn};
  }
 
  double f_root_scr(int k) const {
    const auto sk = std::size_t(k);
    return sk < m_root_fk.size() ? m_root_fk[sk] : 1.0;
  }
 
  double eta_hp(int k) const {
    const auto sk = std::size_t(k);
    return sk < m_etak.size() ? m_etak[sk] : 1.0;
  }
 
  std::pair<double, double>
  BO(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
     const DiracSpinor &v,
     const ExternalField::CorePolarisation *const dV = nullptr, double fw = 1.0,
     double fv = 1.0) const;
 
  Coulomb::meTable<std::pair<double, double>>
  srn_table(const DiracOperator::TensorOperator *const h,
            const std::vector<DiracSpinor> &as,
            const std::vector<DiracSpinor> &tbs = {}, double omega = 0.0,
            const ExternalField::CorePolarisation *const dV = nullptr) const;
 
private:
  // "Top" diagrams
  double t1234(int k, const DiracSpinor &w, const DiracSpinor &r,
               const DiracSpinor &v, const DiracSpinor &c) const;
  double t1(int k, const DiracSpinor &w, const DiracSpinor &r,
            const DiracSpinor &v, const DiracSpinor &c) const;
  double t2(int k, const DiracSpinor &w, const DiracSpinor &r,
            const DiracSpinor &v, const DiracSpinor &c) const;
  double t3(int k, const DiracSpinor &w, const DiracSpinor &r,
            const DiracSpinor &v, const DiracSpinor &c) const;
  double t4(int k, const DiracSpinor &w, const DiracSpinor &r,
            const DiracSpinor &v, const DiracSpinor &c) const;
 
  // "Bottom" diagrams
  double b1234(int k, const DiracSpinor &w, const DiracSpinor &c,
               const DiracSpinor &v, const DiracSpinor &r) const;
 
  // "Centre" diagrams (most important)
  double c1(int k, const DiracSpinor &w, const DiracSpinor &a,
            const DiracSpinor &v, const DiracSpinor &c) const;
  double c2(int k, const DiracSpinor &w, const DiracSpinor &a,
            const DiracSpinor &v, const DiracSpinor &c) const;
  double d1(int k, const DiracSpinor &w, const DiracSpinor &r,
            const DiracSpinor &v, const DiracSpinor &m) const;
  double d2(int k, const DiracSpinor &w, const DiracSpinor &r,
            const DiracSpinor &v, const DiracSpinor &m) const;
 
  // "Normalisation" terms (most important)
  double n1(const DiracSpinor &v) const;
  double n2(const DiracSpinor &v) const;
  double dSigma_dE(const DiracSpinor &v, const DiracSpinor &i,
                   const DiracSpinor &j, const DiracSpinor &k) const;
 
  std::pair<double, double>
  z_bo(const DiracOperator::TensorOperator *const h, const DiracSpinor &w,
       const DiracSpinor &v, bool transpose = false,
       const ExternalField::CorePolarisation *const dV = nullptr) const;
 
  //
  double Q(int k, const DiracSpinor &a, const DiracSpinor &b,
           const DiracSpinor &c, const DiracSpinor &d) const {
    // inlcudes effective Coulomb screening
    const auto fk = f_root_scr(k);
    return m_use_Qk ? fk * mQ->Q(k, a, b, c, d) : fk * mY.Q(k, a, b, c, d);
  }
 
  double P(int k, const DiracSpinor &a, const DiracSpinor &b,
           const DiracSpinor &c, const DiracSpinor &d) const {
    // inlcudes effective Coulomb screening
    return m_use_Qk ? mQ->P2(k, a, b, c, d, mY.SixJ(), m_root_fk) :
                      mY.P2(k, a, b, c, d, m_root_fk);
  }
 
  double W(int k, const DiracSpinor &a, const DiracSpinor &b,
           const DiracSpinor &c, const DiracSpinor &d) const {
    // inlcudes effective Coulomb screening
    return Q(k, a, b, c, d) + P(k, a, b, c, d);
  }
};
 
} // namespace MBPT