Fix ACE law 4/44/61 sampling for discrete lines and for missing continuum

src/distribution_energy.cpp

@@ -213,8 +213,13 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const
   double E_l_k = distribution_[l].e_out[k];
 
   if (k < n_discrete) {
-    // Discrete case
-    return E_l_k;
+    // Discrete case. Discrete lines correspond by index between adjacent
+    // incident energies, so interpolate the line energy rather than returning
+    // the value tabulated at grid l alone. Histogram interpolation on the
+    // incident energy grid pins l to i, so no interpolation is applied there.
+    const auto& e_lo = distribution_[i].e_out;
+    const auto& e_hi = distribution_[i + 1].e_out;
+    return e_lo[k] + (histogram_interp ? 0.0 : r) * (e_hi[k] - e_lo[k]);
   } else {
     // Continuous case
     double p_l_k = distribution_[l].p[k];
@@ -251,30 +256,36 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const
         "distribution."};
     }
 
-    // Now interpolate between incident energy bins i and i + 1
-    if (!histogram_interp && n_energy_out > 1) {
-      // Interpolation for energy E1 and EK
-      n_energy_out = distribution_[i].e_out.size();
-      n_discrete = distribution_[i].n_discrete;
-      const double E_i_1 = distribution_[i].e_out[n_discrete];
-      const double E_i_K = distribution_[i].e_out[n_energy_out - 1];
-
-      n_energy_out = distribution_[i + 1].e_out.size();
-      n_discrete = distribution_[i + 1].n_discrete;
-      const double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
-      const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
-
-      const double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
-      const double E_K = E_i_K + r * (E_i1_K - E_i_K);
-
-      if (l == i) {
-        return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
-      } else {
-        return E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
-      }
-    } else {
+    // Histogram interpolation on the incident energy grid means the table at i
+    // applies as tabulated, with no interpolation between grids
+    if (histogram_interp)
       return E_out;
-    }
+
+    // Unit-base interpolation between the continua at i and i + 1. The ACE
+    // format requires the number of discrete lines to match at every incident
+    // energy, but does not require a continuum to be present at each one, so
+    // check for each side and fall back to whichever one has a continuum.
+    const int n_i = distribution_[i].e_out.size();
+    const int nd_i = distribution_[i].n_discrete;
+    const int n_i1 = distribution_[i + 1].e_out.size();
+    const int nd_i1 = distribution_[i + 1].n_discrete;
+    const bool cont_i = nd_i < n_i;
+    const bool cont_i1 = nd_i1 < n_i1;
+
+    const double r_c = cont_i ? (cont_i1 ? r : 0.0) : 1.0;
+    const double E_i_1 = cont_i ? distribution_[i].e_out[nd_i] : 0.0;
+    const double E_i_K = cont_i ? distribution_[i].e_out[n_i - 1] : 0.0;
+    const double E_i1_1 = cont_i1 ? distribution_[i + 1].e_out[nd_i1] : 0.0;
+    const double E_i1_K = cont_i1 ? distribution_[i + 1].e_out[n_i1 - 1] : 0.0;
+
+    const double E_1 = E_i_1 + r_c * (E_i1_1 - E_i_1);
+    const double E_K = E_i_K + r_c * (E_i1_K - E_i_K);
+
+    // Now interpolate between incident energy bins i and i + 1
+    const double E_lo = (l == i) ? E_i_1 : E_i1_1;
+    const double E_hi = (l == i) ? E_i_K : E_i1_K;
+    return (E_hi != E_lo) ? E_1 + (E_out - E_lo) * (E_K - E_1) / (E_hi - E_lo)
+                          : E_out;
   }
 }
 

src/secondary_correlated.cpp

@@ -166,23 +166,9 @@ Distribution& CorrelatedAngleEnergy::sample_dist(
   // Sample between the ith and [i+1]th bin
   int l = r > prn(seed) ? i + 1 : i;
 
-  // Interpolation for energy E1 and EK
-  int n_energy_out = distribution_[i].e_out.size();
-  int n_discrete = distribution_[i].n_discrete;
-  double E_i_1 = distribution_[i].e_out[n_discrete];
-  double E_i_K = distribution_[i].e_out[n_energy_out - 1];
-
-  n_energy_out = distribution_[i + 1].e_out.size();
-  n_discrete = distribution_[i + 1].n_discrete;
-  double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
-  double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
-
-  double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
-  double E_K = E_i_K + r * (E_i1_K - E_i_K);
-
   // Determine outgoing energy bin
-  n_energy_out = distribution_[l].e_out.size();
-  n_discrete = distribution_[l].n_discrete;
+  int n_energy_out = distribution_[l].e_out.size();
+  int n_discrete = distribution_[l].n_discrete;
   double r1 = prn(seed);
   double c_k = distribution_[l].c[0];
   int k = 0;
@@ -211,8 +197,12 @@ Distribution& CorrelatedAngleEnergy::sample_dist(
 
   double E_l_k = distribution_[l].e_out[k];
   if (k < n_discrete) {
-    // Discrete case
-    E_out = E_l_k;
+    // Discrete case. Discrete lines correspond by index between adjacent
+    // incident energies, so interpolate the line energy rather than using the
+    // value tabulated at grid l alone.
+    const auto& e_lo = distribution_[i].e_out;
+    const auto& e_hi = distribution_[i + 1].e_out;
+    E_out = e_lo[k] + r * (e_hi[k] - e_lo[k]);
   } else {
     // Continuous case
     double p_l_k = distribution_[l].p[k];
@@ -240,11 +230,31 @@ Distribution& CorrelatedAngleEnergy::sample_dist(
       }
     }
 
+    // Unit-base interpolation between the continua at i and i + 1. The ACE
+    // format requires the number of discrete lines to match at every incident
+    // energy, but does not require a continuum to be present at each one, so
+    // check for each side and fall back to whichever one has a continuum.
+    const int n_i = distribution_[i].e_out.size();
+    const int nd_i = distribution_[i].n_discrete;
+    const int n_i1 = distribution_[i + 1].e_out.size();
+    const int nd_i1 = distribution_[i + 1].n_discrete;
+    const bool cont_i = nd_i < n_i;
+    const bool cont_i1 = nd_i1 < n_i1;
+
+    const double r_c = cont_i ? (cont_i1 ? r : 0.0) : 1.0;
+    const double E_i_1 = cont_i ? distribution_[i].e_out[nd_i] : 0.0;
+    const double E_i_K = cont_i ? distribution_[i].e_out[n_i - 1] : 0.0;
+    const double E_i1_1 = cont_i1 ? distribution_[i + 1].e_out[nd_i1] : 0.0;
+    const double E_i1_K = cont_i1 ? distribution_[i + 1].e_out[n_i1 - 1] : 0.0;
+
+    const double E_1 = E_i_1 + r_c * (E_i1_1 - E_i_1);
+    const double E_K = E_i_K + r_c * (E_i1_K - E_i_K);
+
     // Now interpolate between incident energy bins i and i + 1
-    if (l == i) {
-      E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
-    } else {
-      E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
+    const double E_lo = (l == i) ? E_i_1 : E_i1_1;
+    const double E_hi = (l == i) ? E_i_K : E_i1_K;
+    if (E_hi != E_lo) {
+      E_out = E_1 + (E_out - E_lo) * (E_K - E_1) / (E_hi - E_lo);
     }
   }
 

src/secondary_kalbach.cpp

@@ -125,23 +125,9 @@ void KalbachMann::sample_params(
   // Sample between the ith and [i+1]th bin
   int l = r > prn(seed) ? i + 1 : i;
 
-  // Interpolation for energy E1 and EK
-  int n_energy_out = distribution_[i].e_out.size();
-  int n_discrete = distribution_[i].n_discrete;
-  double E_i_1 = distribution_[i].e_out[n_discrete];
-  double E_i_K = distribution_[i].e_out[n_energy_out - 1];
-
-  n_energy_out = distribution_[i + 1].e_out.size();
-  n_discrete = distribution_[i + 1].n_discrete;
-  double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
-  double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
-
-  double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
-  double E_K = E_i_K + r * (E_i1_K - E_i_K);
-
   // Determine outgoing energy bin
-  n_energy_out = distribution_[l].e_out.size();
-  n_discrete = distribution_[l].n_discrete;
+  int n_energy_out = distribution_[l].e_out.size();
+  int n_discrete = distribution_[l].n_discrete;
   double r1 = prn(seed);
   double c_k = distribution_[l].c[0];
   int k = 0;
@@ -170,8 +156,12 @@ void KalbachMann::sample_params(
 
   double E_l_k = distribution_[l].e_out[k];
   if (k < n_discrete) {
-    // Discrete case
-    E_out = E_l_k;
+    // Discrete case. Discrete lines correspond by index between adjacent
+    // incident energies, so interpolate the line energy rather than using the
+    // value tabulated at grid l alone.
+    const auto& e_lo = distribution_[i].e_out;
+    const auto& e_hi = distribution_[i + 1].e_out;
+    E_out = e_lo[k] + r * (e_hi[k] - e_lo[k]);
     km_r = distribution_[l].r[k];
     km_a = distribution_[l].a[k];
 
@@ -215,11 +205,31 @@ void KalbachMann::sample_params(
                (distribution_[l].a[k + 1] - distribution_[l].a[k]);
     }
 
+    // Unit-base interpolation between the continua at i and i + 1. The ACE
+    // format requires the number of discrete lines to match at every incident
+    // energy, but does not require a continuum to be present at each one, so
+    // check for each side and fall back to whichever one has a continuum.
+    const int n_i = distribution_[i].e_out.size();
+    const int nd_i = distribution_[i].n_discrete;
+    const int n_i1 = distribution_[i + 1].e_out.size();
+    const int nd_i1 = distribution_[i + 1].n_discrete;
+    const bool cont_i = nd_i < n_i;
+    const bool cont_i1 = nd_i1 < n_i1;
+
+    const double r_c = cont_i ? (cont_i1 ? r : 0.0) : 1.0;
+    const double E_i_1 = cont_i ? distribution_[i].e_out[nd_i] : 0.0;
+    const double E_i_K = cont_i ? distribution_[i].e_out[n_i - 1] : 0.0;
+    const double E_i1_1 = cont_i1 ? distribution_[i + 1].e_out[nd_i1] : 0.0;
+    const double E_i1_K = cont_i1 ? distribution_[i + 1].e_out[n_i1 - 1] : 0.0;
+
+    const double E_1 = E_i_1 + r_c * (E_i1_1 - E_i_1);
+    const double E_K = E_i_K + r_c * (E_i1_K - E_i_K);
+
     // Now interpolate between incident energy bins i and i + 1
-    if (l == i) {
-      E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
-    } else {
-      E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
+    const double E_lo = (l == i) ? E_i_1 : E_i1_1;
+    const double E_hi = (l == i) ? E_i_K : E_i1_K;
+    if (E_hi != E_lo) {
+      E_out = E_1 + (E_out - E_lo) * (E_K - E_1) / (E_hi - E_lo);
     }
   }
 }

tests/regression_tests/weightwindows/shared/results_true.dat

@@ -1 +1 @@
-19b676a5e0f2784747d95fe884a44af3b31b4f5193e5ac514378d6d83cfa8ae9f4722da408a98b9b9cd3e369517d0b29bc523fd6f5492856f457e14d15007fa0
+3916cc59d580f8caf1dac2ae30d7a5926f650ab1e314af4ecedd3521c87aa7ab60f3ddfee940ed221dbfd2a40deb0f726eb6252ca3333445c7aa8155e7003094