Ek.hpp

#pragma once
#include "DiracOperator/TensorOperator.hpp"
#include "IO/InputBlock.hpp"
#include "Wavefunction/Wavefunction.hpp"
 
namespace DiracOperator {
 
//==============================================================================
 
class Ek : public TensorOperator {
public:
  Ek(const Grid &gr, const int k)
      : TensorOperator(k, Angular::evenQ(k) ? Parity::even : Parity::odd, -1.0,
                       gr.rpow(k), 0),
        m_k(k) {}
  double angularF(const int ka, const int kb) const override final {
    return Angular::Ck_kk(m_k, ka, kb);
  }
  std::string name() const override {
    return std::string("E") + std::to_string(m_k);
  }
  std::string units() const override {
    return m_k == 1 ? "|e|aB" : std::string("|e|aB^") + std::to_string(m_k);
  }
 
private:
  int m_k;
};
 
//==============================================================================
 
class E1 final : public Ek {
public:
  E1(const Grid &gr) : Ek(gr, 1) {}
};
 
//==============================================================================
class sigma_r final : public ScalarOperator {
public:
  sigma_r(const Grid &rgrid) : ScalarOperator(Parity::odd, -1.0, rgrid.r()) {}
  std::string name() const override final { return "s.r"; }
  std::string units() const override final { return "aB"; }
};
 
//==============================================================================
 
class E1v final : public TensorOperator
// d_v = (ie/w alpha) v{alpha}   [v{a} = g0v{g}]\f$
// <a||dv||b> = -2e/(w alpha) Int[ fagb <ka||s||-kb> - gafb <-ka||s||kb>]
{
public:
  E1v(const double alpha, const double omega = 0.0)
      : TensorOperator(1, Parity::odd, -0.0, {}, 0, Realness::real, true),
        m_alpha(alpha) {
    updateFrequency(omega);
  }
  std::string name() const override final { return "E1v"; }
  std::string units() const override final { return "|e|aB"; }
 
  double angularF(const int, const int) const override final { return 1.0; }
 
  double angularCff(int, int) const override final { return 0; }
  double angularCgg(int, int) const override final { return 0; }
  double angularCfg(int ka, int kb) const override final {
    return Angular::S_kk(ka, -kb);
  }
  double angularCgf(int ka, int kb) const override final {
    return -Angular::S_kk(-ka, kb);
  }
 
  void updateFrequency(const double omega) override final {
    m_constant = std::abs(omega) > 1.0e-10 ? -2.0 / (m_alpha * omega) : 1.0;
  }
 
private:
  double m_alpha; // (including var-alpha)
};
 
//==============================================================================
 
inline std::unique_ptr<DiracOperator::TensorOperator>
generate_sigma_r(const IO::InputBlock &input, const Wavefunction &wf) {
  using namespace DiracOperator;
  input.check({{"no options", ""}});
  if (input.has_option("help")) {
    return nullptr;
  }
  return std::make_unique<sigma_r>(wf.grid());
}
 
inline std::unique_ptr<DiracOperator::TensorOperator>
generate_E1(const IO::InputBlock &input, const Wavefunction &wf) {
  using namespace DiracOperator;
  input.check({{"no options", ""}});
  if (input.has_option("help")) {
    return nullptr;
  }
  return std::make_unique<E1>(wf.grid());
}
 
inline std::unique_ptr<DiracOperator::TensorOperator>
generate_E1v(const IO::InputBlock &input, const Wavefunction &wf) {
  using namespace DiracOperator;
  input.check({{"no options", ""}});
  if (input.has_option("help")) {
    return nullptr;
  }
  return std::make_unique<E1v>(wf.alpha(), 0.0);
}
 
inline std::unique_ptr<DiracOperator::TensorOperator>
generate_E2(const IO::InputBlock &input, const Wavefunction &wf) {
  using namespace DiracOperator;
  input.check({{"no options", ""}});
  if (input.has_option("help")) {
    return nullptr;
  }
  return std::make_unique<Ek>(wf.grid(), 2);
}
 
//------------------------------------------------------------------------------
inline std::unique_ptr<DiracOperator::TensorOperator>
generate_Ek(const IO::InputBlock &input, const Wavefunction &wf) {
  using namespace DiracOperator;
  input.check({{"k", "Rank: k=1 for E1, =2 for E2 etc. [1]"}});
  if (input.has_option("help")) {
    return nullptr;
  }
  const auto k = input.get("k", 1);
  return std::make_unique<Ek>(wf.grid(), k);
}
 
} // namespace DiracOperator