Updated - BUT WRONG???

src/DiracOperator/Operators/hfs.hpp

@@ -291,8 +291,8 @@ generate_hfs(const IO::InputBlock &input, const Wavefunction &wf) {
 
   const auto nuc = wf.nucleus();
   const auto isotope = Nuclear::findIsotopeData(nuc.z(), nuc.a());
-  const auto mu = input.get("mu", isotope.mu);
-  const auto I_nuc = input.get("I", isotope.I_N);
+  const auto mu = input.get("mu", isotope.mu ? *isotope.mu : 0.0);
+  const auto I_nuc = input.get("I", isotope.I_N ? *isotope.I_N : 0.0);
   const auto print = input.get("print", true);
   const auto k = input.get("k", 1);
 
@@ -372,7 +372,7 @@ generate_hfs(const IO::InputBlock &input, const Wavefunction &wf) {
   } else if (distro_type == DistroType::shell) {
     Fr = Hyperfine::sphericalShell_F();
   } else if (distro_type == DistroType::SingleParticle) {
-    const auto pi = input.get("parity", isotope.parity);
+    const auto pi = input.get("parity", isotope.parity ? *isotope.parity : 0);
     const auto l_tmp = int(I_nuc + 0.5 + 0.0001);
     auto l = ((l_tmp % 2 == 0) == (pi == 1)) ? l_tmp : l_tmp - 1;
     l = input.get("l", l); // can override derived 'l' (not recommended)

src/Modules/HFAnomaly.cpp

@@ -215,23 +215,6 @@ void HFAnomaly(const IO::InputBlock &input, const Wavefunction &wf) {
         "List: For 1D2_target only: minimum and maximum magnetic "
         "radii, as a fraction of charge radius, and number of "
         "steps in the fit. Default=[0.9,1.5,10]"}});
-
-  std::cout << "\n";
-  double res = 0.0;
-  for (const auto &element : Nuclear::NuclearDataTable) {
-
-    auto r0 = element.r_rms;
-    if (r0 <= 0.0)
-      continue;
-    // if (element.A % 2 == 0)
-    //   continue;
-    auto rt = Nuclear::approximate_r_rms(element.A, element.Z);
-    auto eps = std::pow((rt - r0) / r0, 2);
-    std::cout << element.Z << " " << element.A << " " << rt << "\n";
-    res += eps;
-  }
-  std::cout << res << "\n";
-
   if (input.has_option("help")) {
     return;
   }
@@ -333,17 +316,17 @@ void HFAnomaly(const IO::InputBlock &input, const Wavefunction &wf) {
       nucleus.t() = *t;
     }
     if (rrms) {
-      nucleus.r_rms() = *rrms;
+      nucleus.set_rrms(*rrms);
     }
     if (c_hdr) {
       // this will over-ride given rms
-      nucleus.r_rms() = Nuclear::rrms_formula_c_t(*c_hdr, nucleus.t());
+      nucleus.set_rrms(Nuclear::rrms_formula_c_t(*c_hdr, nucleus.t()));
     }
     // If A or given rrms are zero, explicitely set to pointlike nucleus
     // This isn't required, but makes output more explicit
     if (nucleus.a() == 0.0 || nucleus.r_rms() == 0.0) {
       nucleus.t() = 0.0;
-      nucleus.r_rms() = 0.0;
+      nucleus.set_rrms(0.0);
       nucleus.type() = Nuclear::ChargeDistro::point;
     }
 

src/Modules/isotopeShift.cpp

@@ -54,7 +54,7 @@ void fieldShift(const IO::InputBlock &input, const Wavefunction &wf) {
       const auto dr2 = rB * rB - r0 * r0;
 
       auto nuc_b = wf.nucleus();
-      nuc_b.r_rms() = rB;
+      nuc_b.set_rrms(rB);
 
       wfB.update_Vnuc(Nuclear::formPotential(nuc_b, wf.grid().r()));
 

src/Physics/NuclearData.cpp

@@ -4,12 +4,12 @@ namespace Nuclear {
 //==============================================================================
 Isotope findIsotopeData(int z, int a) {
   if (a == 0)
-    return Isotope{z, a, 0, 0, 0, 0};
+    return Isotope{z, 0, {}, {}, {}, {}, {}};
   for (const auto &nucleus : NuclearDataTable) {
     if (nucleus.Z == z && nucleus.A == a)
       return nucleus;
   }
-  return Isotope{z, a, -1, 0, 0, -1};
+  return Isotope{z, a, {}, {}, {}, {}, {}};
 }
 
 std::vector<Isotope> findIsotopeList(int z) {
@@ -22,31 +22,23 @@ std::vector<Isotope> findIsotopeList(int z) {
 }
 
 double find_rrms(int z, int a) {
-  auto nuc = findIsotopeData(z, a);
-  if (!nuc.r_ok())
-    return 0;
-  return nuc.r_rms;
+  const auto nuc = findIsotopeData(z, a);
+  return nuc.r_rms ? *nuc.r_rms : 0.0;
 }
 
 double find_mu(int z, int a) {
   auto nuc = findIsotopeData(z, a);
-  if (!nuc.mu_ok())
-    return 0;
-  return nuc.mu;
+  return nuc.mu ? *nuc.mu : 0.0;
 }
 
 int find_parity(int z, int a) {
-  auto nuc = findIsotopeData(z, a);
-  if (!nuc.parity_ok())
-    return 0;
-  return nuc.parity;
+  const auto nuc = findIsotopeData(z, a);
+  return nuc.parity ? *nuc.parity : 0;
 }
 
 double find_spin(int z, int a) {
-  auto nuc = findIsotopeData(z, a);
-  if (!nuc.I_ok())
-    return -1.0;
-  return nuc.I_N;
+  const auto nuc = findIsotopeData(z, a);
+  return nuc.I_N ? *nuc.I_N : 0.0;
 }
 
 //==============================================================================