PNC.hpp

#pragma once
#include "DiracOperator/TensorOperator.hpp"
#include "IO/InputBlock.hpp"
#include "Wavefunction/Wavefunction.hpp"
 
namespace DiracOperator {
 
//==============================================================================
 
class PNCnsi final : public ScalarOperator {
public:
  PNCnsi(double c, double t, const Grid &rgrid, double factor = 1.0,
         const std::string &in_units = "iQw*e-11")
      : ScalarOperator(Parity::odd, factor * PhysConst::GFe11 / std::sqrt(8.0),
                       Nuclear::fermiNuclearDensity_tcN(t, c, 1.0, rgrid),
                       {0, -1, +1, 0}, 0, Realness::imaginary),
        m_unit(in_units) {}
  std::string name() const override final { return "pnc-nsi"; }
  std::string units() const override final { return m_unit; }
 
private:
  const std::string m_unit{"iQw*e-11"};
};
 
class PNCnsi_const final : public ScalarOperator {
public:
  PNCnsi_const(double Rnuc_au, double factor = 1.0,
               const std::string &in_units = "iQw*e-11")
      : ScalarOperator(Parity::odd,
                       (3.0 / (4.0 * M_PI * Rnuc_au * Rnuc_au * Rnuc_au)) *
                           factor * PhysConst::GFe11 / std::sqrt(8.0),
                       {}, {0, -1, +1, 0}, 0, Realness::imaginary),
        m_unit(in_units) {}
  std::string name() const override final { return "pnc_const"; }
  std::string units() const override final { return m_unit; }
 
private:
  const std::string m_unit{"iQw*e-11"};
};
 
//==============================================================================
inline std::unique_ptr<DiracOperator::TensorOperator>
generate_pnc(const IO::InputBlock &input, const Wavefunction &wf) {
  using namespace DiracOperator;
  input.check(
      {{{"c",
         "Half-density radius for Fermi rho(r). [defaut: from wavefunction]"},
        {"t", "skin thickness [2.3]"},
        {"N", "Neutron number, for units [default: from wavefunction]"},
        {"print", "Write details to screen [true]"}}});
  if (input.has_option("help")) {
    return nullptr;
  }
  const auto r_rms = wf.get_rrms();
  const auto c = input.get("c", Nuclear::c_hdr_formula_rrms_t(r_rms));
  const auto t = input.get("t", Nuclear::default_t);
  const auto N = input.get("N", wf.Anuc() - wf.Znuc());
  std::string units =
      N == 1 ? "i(-Qw)e-11" : "i(-Qw/" + std::to_string(N) + ")e-11";
  if (input.get("print", true))
    std::cout << "pnc: with c=" << c << ", t=" << t << " [" << units << "]\n";
  return std::make_unique<PNCnsi>(c, t, wf.grid(), -1.0 * N, "units");
}
 
} // namespace DiracOperator