M1.hpp

#pragma once
#include "DiracOperator/TensorOperator.hpp"
#include "IO/InputBlock.hpp"
#include "Wavefunction/Wavefunction.hpp"
#include "jls.hpp"
 
namespace DiracOperator {
 
 
class M1 final : public TensorOperator {
public:
  M1(const Grid &gr, const double alpha, const double omega = 0.0)
      : TensorOperator(1, Parity::even, +0.0, {}, 0, Realness::real, true),
        m_r(gr.r()),
        m_alpha(alpha) {
    updateFrequency(omega);
  }
  M1 &operator=(const M1 &) = delete;
  M1(const M1 &) = default;
  ~M1() = default;
  std::string name() const override final { return std::string("M1"); }
  std::string units() const override final { return std::string("|mu_B|"); }
 
  double angularF(const int ka, const int kb) const override final {
    return (ka + kb) * Angular::Ck_kk(1, -ka, kb);
  }
  double angularCff(int, int) const override final { return 0.0; }
  double angularCgg(int, int) const override final { return 0.0; }
  double angularCfg(int, int) const override final { return 1.0; }
  double angularCgf(int, int) const override final { return 1.0; }
 
  void updateFrequency(const double omega) override final {
    if (std::abs(omega) > 1.0e-10) {
      m_constant = -3.0 / std::abs((m_alpha * m_alpha * omega));
      m_vec =
          SphericalBessel::fillBesselVec_kr(1, std::abs(omega) * m_alpha, m_r);
    } else {
      // j1(kr) -> (r*k)/3, for k<<1
      m_constant = -1.0 / (m_alpha);
      m_vec = m_r;
    }
  }
 
private:
  const std::vector<double> m_r; // store radial vector (copy)
  const double m_alpha;
};
 
//==============================================================================
class M1nr final : public TensorOperator {
public:
  M1nr() : TensorOperator(1, Parity::even, 1.0) {}
  std::string name() const override final { return "M1nr"; }
  std::string units() const override final { return "|mu_B|"; }
 
  double angularF(const int, const int) const override final {
    // angular factor included in L and S
    return 1.0;
  }
 
  DiracSpinor radial_rhs(const int kappa_a,
                         const DiracSpinor &Fb) const override final {
    return m_l.reduced_rhs(kappa_a, Fb) + 2.0 * m_s.reduced_rhs(kappa_a, Fb);
  }
 
  double radialIntegral(const DiracSpinor &Fa,
                        const DiracSpinor &Fb) const override final {
    // nb: faster not to do this, but nicer this way
    return Fa * radial_rhs(Fa.kappa(), Fb);
  }
 
private:
  DiracOperator::s m_s{};
  DiracOperator::l m_l{};
};
 
//------------------------------------------------------------------------------
inline std::unique_ptr<DiracOperator::TensorOperator>
generate_M1(const IO::InputBlock &input, const Wavefunction &wf) {
  using namespace DiracOperator;
  input.check({{"", "No input options"}});
  if (input.has_option("help")) {
    return nullptr;
  }
  return std::make_unique<M1>(wf.grid(), wf.alpha(), 0.0);
}
 
inline std::unique_ptr<DiracOperator::TensorOperator>
generate_M1nr(const IO::InputBlock &input, const Wavefunction &) {
  using namespace DiracOperator;
  input.check({{"", "No input options"}});
  if (input.has_option("help")) {
    return nullptr;
  }
  return std::make_unique<M1nr>();
}
 
} // namespace DiracOperator