TensorOperator.hpp

#pragma once
#include "Angular/Wigner369j.hpp"
#include "Maths/Grid.hpp"
#include "Maths/NumCalc_quadIntegrate.hpp"
#include "Maths/SphericalBessel.hpp"
#include "Physics/PhysConst_constants.hpp"
#include "Potentials/NuclearPotentials.hpp"
#include "Wavefunction/DiracSpinor.hpp"
#include "qip/Vector.hpp"
#include <cmath>
#include <string>
#include <vector>
 
//! Dirac Operators: General + derived
namespace DiracOperator {
 
//! Parity of operator
enum class Parity { even, odd, Error };
 
//! Realness of matrix element; impacts symmetry only
enum class Realness { real, imaginary, Error };
 
//! Type of matrix element returned
/*! @details
 Specifies which form of matrix element is used or returned.
 
  @value Reduced    : Reduced matrix element
  @value Stretched  : Matrix element for stretched states (m = j)
  @value HFConstant : Hyperfine (A,B,C...) constant
 
  For off-diagonal, j:=min(ja,jb)
*/
enum class MatrixElementType { Reduced, Stretched, HFConstant, Error };
 
//! Convert string to Parity
Parity parse_Parity(const std::string &s);
 
//! Convert Parity to string
std::string parse_Parity(Parity p);
 
//! Convert string to Realness
Realness parse_Realness(const std::string &s);
 
//! Convert Realness to string
std::string parse_Realness(Realness r);
 
//! Convert string to MatrixElementType
MatrixElementType parse_MatrixElementType(const std::string &s);
 
//! Convert MatrixElementType to string
std::string parse_MatrixElementType(MatrixElementType t);
 
//==============================================================================
//! @brief General operator (virtual base class); operators derive from this.
/*! @details
  - k is rank, c is multiplicative constant, d_order is derivative order,
  - pi is parity, may be Parity::even or ::odd.
  - RorI may be Realness::real or Realness::imaginary.
  - Note: You may not construct a TensorOperator. Instead, you must construct
    one of the derived 'operators' (there are some general ones); see
    operators.hpp for list of operators. Operators work by overrideing the
    angularCxx() functions and angularF().
  - c, v, and Cxx are included in radial integral.
*/
class TensorOperator {
protected:
  //! @brief Constructs a tensor operator description.
  /*! @details
    Initialises the basic properties of a tensor operator.
 
    @param rank_k     Tensor rank k of the operator.
    @param pi         Parity of the operator (Parity::even or Parity::odd).
    @param constant   Multiplicative constant c, included in the radial integral.
    @param vec        Overall v=f(r) radial function, as std::vector array
    @param RorI       Specifies whether the matrix element is real or imaginary
                      (Realness::real or Realness::imaginary).
    @param freq_dep   Indicates whether the operator depends on frequency
                      (e.g., dynamic external fields).
 
    @note TensorOperator is a virtual base class and may not be constructed
    directly. It is intended to be initialised by derived operator classes.
  */
  TensorOperator(int rank_k, Parity pi, double constant = 1.0,
                 const std::vector<double> &vec = {}, int diff_order = 0,
                 Realness RorI = Realness::real, bool freq_dep = false)
    : m_rank(rank_k),
      m_parity(pi),
      m_diff_order(diff_order),
      opC(RorI),
      m_freqDependantQ(freq_dep),
      m_constant(constant),
      m_vec(vec){};
 
public:
  virtual ~TensorOperator() = default;
  TensorOperator(const TensorOperator &) = default;
  TensorOperator &operator=(const TensorOperator &) = default;
  TensorOperator(TensorOperator &&) = default;
  TensorOperator &operator=(TensorOperator &&) = default;
 
protected:
  int m_rank;
  Parity m_parity;
  int m_diff_order;
  Realness opC;
  bool m_freqDependantQ{false};
 
protected:
  // these may be updated for frequency-dependant operators
  double m_constant; // included in radial integral
  std::vector<double> m_vec;
 
public:
  bool freqDependantQ() const { return m_freqDependantQ; }
 
public:
  //! If matrix element <a|h|b> is zero, returns true
  bool isZero(const int ka, int kb) const;
  bool isZero(const DiracSpinor &Fa, const DiracSpinor &Fb) const;
 
  bool selectrion_rule(int twoJA, int piA, int twoJB, int piB) const {
    if (twoJA == twoJB && twoJA == 0.0)
      return false;
 
    if (Angular::triangle(twoJA, twoJB, 2 * m_rank) == 0)
      return false;
 
    // if (2 * m_rank < std::abs(twoJA - twoJB))
    //   return false;
    return (m_parity == Parity::even) == (piA == piB);
  }
 
  //! Update frequency for frequency-dependant operators.
  virtual void updateFrequency(const double){};
 
  //! Permanently re-scales the operator by constant, lambda
  void scale(double lambda);
 
  //! Returns a const ref to vector v
  const std::vector<double> &getv() const { return m_vec; }
 
  //! Returns a const ref to constant c
  double getc() const { return m_constant; }
 
  int get_d_order() const { return m_diff_order; }
 
  //! returns true if operator is imaginary (has imag MEs)
  bool imaginaryQ() const { return (opC == Realness::imaginary); }
 
  //! Rank k of operator
  int rank() const { return m_rank; }
 
  //! returns parity, as integer (+1 or -1)
  int parity() const { return (m_parity == Parity::even) ? 1 : -1; }
 
  //! returns relative sign between <a||x||b> and <b||x||a>
  int symm_sign(const DiracSpinor &Fa, const DiracSpinor &Fb) const {
    const auto sra_i = imaginaryQ() ? -1 : 1;
    const auto sra = Angular::neg1pow_2(Fa.twoj() - Fb.twoj());
    return sra_i * sra;
  }
 
  //! Returns "name" of operator (e.g., 'E1')
  virtual std::string name() const { return "Operator"; };
  //! Returns units of operator (usually au, may be MHz, etc.)
  virtual std::string units() const { return "au"; };
 
public:
  // These are needed for radial integrals
  // Usually just constants, but can also be functions of kappa
 
  //! Angular factor for f_a*f_b part of radial integral
  /*! @details 
    - Often constant, though sometimes depends on \f$\kappa_a,\kappa_b\f$
    \f[
      \int_0^\infty v(r)\left(
        C_{ff}f_a(r)f_b(r) + 
        C_{fg}f_a(r)g_b(r) + 
        C_{gf}g_a(r)f_b(r) + 
        C_{gg}g_a(r)g_b(r)
      \right)\,{\rm d}r
    \f]
  */
  virtual double angularCff(int /*k_a*/, int /*k_b*/) const { return 1.0; }
  //! Angular factor for g_a*g_b part of radial integral
  virtual double angularCgg(int, int) const { return 1.0; }
  //! Angular factor for f_a*g_b part of radial integral
  virtual double angularCfg(int, int) const { return 0.0; }
  //! Angular factor for g_a*f_b part of radial integral
  virtual double angularCgf(int, int) const { return 0.0; }
 
public:
  //! @brief angularF: links radiation integral to RME.
  //! RME = <a||h||b> = angularF(a,b) * radial_int(a,b)
  virtual double angularF(const int, const int) const = 0;
 
  //! radial_int = Fa * radial_rhs(a, Fb) (a needed for angular factor)
  /*! @details
  
  By default, is defined as:
 
    \f[
      \delta F_{\rm b} = 
      C_{\rm overall}\,v(r)
      \begin{pmatrix}
        C_{ff}f_b(r) + C_{fg}g_b(r) \\
        C_{gf}f_b(r) + C_{gg}g_b(r)
      \end{pmatrix}
    \f]
 
    \f$ C_{f/g} \f$ are angular coeficients - see \ref TensorOperator::angularCff
 
    See also: \ref TensorOperator::radial_rhs
  */
  virtual DiracSpinor radial_rhs(const int kappa_a,
                                 const DiracSpinor &Fb) const;
 
  //! Radial part of integral R_ab = (Fa|t|Fb).
  /*! @details
 
    <a||t||b> = angularF(a,b) * radialIntegral(a,b)
 
    The 'radial'  integral \f$R_{ab}\f$ for operator \f$t\f$ is defined as
    \f[
      \langle a ||t|| b \rangle \equiv R_{ab} * A_{ab}
    \f]
 
    where \f$A_{ab}\f$ is TensorOperator::angularF()
 
    By default, is implemented as:
    \f[
      R = \int_0^\infty F_a\cdot\delta F_{\rm b} \,{\rm d}r
    \f]
 
    where \f$ \delta F_{\rm b} \f$ is defined in \ref TensorOperator::radial_rhs
 
    So, if \ref TensorOperator::radial_rhs is defined normally, we have:
    \f[
      C_{\rm overall}\int_0^\infty v(r)\left(
        C_{ff}f_a(r)f_b(r) + 
        C_{fg}f_a(r)g_b(r) + 
        C_{gf}g_a(r)f_b(r) + 
        C_{gg}g_a(r)g_b(r)
      \right)\,{\rm d}r
    \f]
 
    \f$ C_{f/g} \f$ are angular coeficients - see \ref TensorOperator::angularCff
 
    @warning
    - This is defined for the "standard" case, which doesn't always suit
    - \ref TensorOperator::radial_rhs may be overridden for special cases
    - If radial_rhs is overridden, then radialIntegral must also be_
  */
  virtual double radialIntegral(const DiracSpinor &Fa,
                                const DiracSpinor &Fb) const;
 
  //! ME = rme3js * RME
  double rme3js(const int twoja, const int twojb, int two_mb = 1,
                int two_q = 0) const;
 
  //! ME = rme3js * RME
  double rme3js(const DiracSpinor &Fa, const DiracSpinor &Fb, int two_mb = 1,
                int two_q = 0) const {
    return rme3js(Fa.twoj(), Fb.twoj(), two_mb, two_q);
  }
 
  //! <a||h||b> = Fa * reduced_rhs(a, Fb) (a needed for angular factor)
  DiracSpinor reduced_rhs(const int ka, const DiracSpinor &Fb) const;
 
  //! <b||h||a>  = Fa * reduced_lhs(a, Fb) (a needed for angular factor)
  DiracSpinor reduced_lhs(const int ka, const DiracSpinor &Fb) const;
 
  //! The reduced matrix element
  double reducedME(const DiracSpinor &Fa, const DiracSpinor &Fb) const;
 
  //! Returns "full" matrix element, for optional (ma, mb, q) [taken as int 2*].
  //! If not specified, returns z-component (q=0), with ma=mb=min(ja,jb)
  double fullME(const DiracSpinor &Fa, const DiracSpinor &Fb,
                std::optional<int> two_ma = std::nullopt,
                std::optional<int> two_mb = std::nullopt,
                std::optional<int> two_q = std::nullopt) const;
 
  //! Converts reduced matrix element to different "type" (MatrixElementType)
  double matel_factor(MatrixElementType type, int twoJa, int twoJb) const;
 
  double matel_factor(MatrixElementType type, const DiracSpinor &Fa,
                      const DiracSpinor &Fb) const {
    return matel_factor(type, Fa.twoj(), Fb.twoj());
  }
};
 
//============================================================================
//============================================================================
//! Speacial case for scalar operator
class ScalarOperator : public TensorOperator {
public:
  ScalarOperator(Parity pi, double in_coef,
                 const std::vector<double> &in_v = {},
                 const std::array<int, 4> &in_g = {1, 0, 0, 1}, int in_diff = 0,
                 Realness rori = Realness::real)
    : TensorOperator(0, pi, in_coef, in_v, in_diff, rori),
      c_ff(in_g[0]),
      c_fg(in_g[1]),
      c_gf(in_g[2]),
      c_gg(in_g[3]) {}
 
  ScalarOperator(const std::vector<double> &in_v = {}, double in_coef = 1.0)
    : TensorOperator(0, Parity::even, in_coef, in_v, 0),
      c_ff(1.0),
      c_fg(0.0),
      c_gf(0.0),
      c_gg(1.0) {}
 
public:
  virtual double angularF(const int ka, const int kb) const override {
    // For scalar operators, <a||h||b> = RadInt / 3js
    // 3js:= 1/(Sqrt[2j+1]) ... depends on m???
    // |k| = j+1/2
    return (std::abs(ka) == std::abs(kb)) ? std::sqrt(2.0 * std::abs(ka)) : 0.0;
  }
 
private:
  const double c_ff, c_fg, c_gf, c_gg;
 
protected:
  double virtual angularCff(int, int) const override { return c_ff; }
  double virtual angularCgg(int, int) const override { return c_gg; }
  double virtual angularCfg(int, int) const override { return c_fg; }
  double virtual angularCgf(int, int) const override { return c_gf; }
};
 
//------------------------------------------------------------------------------
//! Speacial operator: 0
class NullOperator final : public ScalarOperator {
public:
  NullOperator() : ScalarOperator(Parity::even, 0, {}) {}
 
protected:
  double virtual angularCff(int, int) const override final { return 0.0; }
  double virtual angularCgg(int, int) const override final { return 0.0; }
  double virtual angularCfg(int, int) const override final { return 0.0; }
  double virtual angularCgf(int, int) const override final { return 0.0; }
};
 
//******************************************************************************
// Helper functions: Useful for several operators
//******************************************************************************
 
//! Pab function: Int[ (fa*gb + pm*ga*fb) * t(r) , dr]. pm = +/-1 (usually)
/*! @details
Note: does not know selection rules, so only evaluate if non-zero
*/
double Pab(double pm, const std::vector<double> &t, const DiracSpinor &Fa,
           const DiracSpinor &Fb);
 
//! Rab function: Int[ (fa*fb + pm*ga*gb) * t(r) , dr]. pm = +/-1 (usually)
/*! @details
Note: does not know selection rules, so only evaluate if non-zero
*/
double Rab(double pm, const std::vector<double> &t, const DiracSpinor &Fa,
           const DiracSpinor &Fb);
 
//! Pab_rhs function: dF_ab += A * t(r) * (g, pm*f) , pm=+/-1 (usually).
//! NOTE: uses +=, so can combine. Ensure empty to begin.
/*! @details
Note: does not know selection rules, so only evaluate if non-zero.
Should have Fa*Pab_rhs = A * Pab.
*/
void Pab_rhs(double pm, const std::vector<double> &t, DiracSpinor *dF,
             const DiracSpinor &Fb, double A = 1.0);
 
//! Rab_rhs function: dF_ab += A * t(r) * (f, pm*g) , pm=+/-1 (usually).
//! NOTE: uses +=, so can combine. Ensure empty to begin.
/*! @details
Note: does not know selection rules, so only evaluate if non-zero.
Should have Fa*Rab_rhs = A * Rab.
*/
void Rab_rhs(double pm, const std::vector<double> &t, DiracSpinor *dF,
             const DiracSpinor &Fb, double A = 1.0);
 
//! Vab function: Int[ (fa*gb ) * t(r) , dr].
double Vab(const std::vector<double> &t, const DiracSpinor &Fa,
           const DiracSpinor &Fb);
//! Wab function: Int[ (ga*fb ) * t(r) , dr].
double Wab(const std::vector<double> &t, const DiracSpinor &Fa,
           const DiracSpinor &Fb);
 
//! Gab function: Int[ (ga*gb ) * t(r) , dr].
double Gab(const std::vector<double> &t, const DiracSpinor &Fa,
           const DiracSpinor &Fb);
 
void Gab_rhs(const std::vector<double> &t, DiracSpinor *dF,
             const DiracSpinor &Fb, double a);
 
// Same - for constant t(r)=c
 
//! Pab[1] function: Int[ (fa*gb + pm*ga*fb) , dr]. pm = +/-1 (usually)
double Pab(double pm, const DiracSpinor &Fa, const DiracSpinor &Fb);
 
//! Rab[1] function: Int[ (fa*fb + pm*ga*gb) , dr] = Int[ (pm-1)*ga*gb) , dr].
//! NOTE: assumes NOT diagonal, using orthogonality condition.
double Rab(double pm, const DiracSpinor &Fa, const DiracSpinor &Fb);
 
//! Pab_rhs[1] function: dF_ab += A  * (g, pm*f) , pm=+/-1 (usually).
void Pab_rhs(double pm, DiracSpinor *dF, const DiracSpinor &Fb, double A = 1.0);
 
//! Rab_rhs[1] function: dF_ab += A * (f, pm*g)  = dF_ab += A * (0, (pm-1)*g).
//! NOTE: assumes NOT diagonal, using orthogonality condition.
void Rab_rhs(double pm, DiracSpinor *dF, const DiracSpinor &Fb, double A = 1.0);
 
//! Gab = Int[ ga*gb , dr] - (just relativistic correction part of integral)
double Gab(const DiracSpinor &Fa, const DiracSpinor &Fb);
 
//! Gab_rhs(r) += a*g_b(r). Note: uses += so may be sumulative
void Gab_rhs(DiracSpinor *dF, const DiracSpinor &Fb, double a);
 
} // namespace DiracOperator