Simplify get_energy_index and use it in place of duplicated lookup logic

src/distribution_angle.cpp

@@ -88,13 +88,8 @@ double AngleDistribution::evaluate(double E, double mu) const
   int i;
   double r;
   get_energy_index(energy_, E, i, r);
-
-  double pdf = 0.0;
-  if (r > 0.0)
-    pdf += r * distribution_[i + 1]->evaluate(mu);
-  if (r < 1.0)
-    pdf += (1.0 - r) * distribution_[i]->evaluate(mu);
-  return pdf;
+  return r * distribution_[i + 1]->evaluate(mu) +
+         (1.0 - r) * distribution_[i]->evaluate(mu);
 }
 
 } // namespace openmc

src/distribution_energy.cpp

@@ -11,7 +11,6 @@
 #include "openmc/math_functions.h"
 #include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
-#include "openmc/search.h"
 
 namespace openmc {
 
@@ -159,19 +158,9 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const
 
   // Find energy bin and calculate interpolation factor -- if the energy is
   // outside the range of the tabulated energies, choose the first or last bins
-  auto n_energy_in = energy_.size();
   int i;
   double r;
-  if (E < energy_[0]) {
-    i = 0;
-    r = 0.0;
-  } else if (E > energy_[n_energy_in - 1]) {
-    i = n_energy_in - 2;
-    r = 1.0;
-  } else {
-    i = lower_bound_index(energy_.begin(), energy_.end(), E);
-    r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);
-  }
+  get_energy_index(energy_, E, i, r);
 
   // Sample between the ith and [i+1]th bin
   int l;

src/math_functions.cpp

@@ -991,13 +991,22 @@ double cyl_bessel_j(int n, double x)
 void get_energy_index(
   const vector<double>& energies, double E, int& i, double& f)
 {
-  // Get index and interpolation factor for linear-linear energy grid
+  // Get index and interpolation factor for linear-linear energy grid. The index
+  // is kept within the topmost interval so that both energies[i] and
+  // energies[i + 1] are valid for callers.
+  const int n = energies.size();
   i = 0;
   f = 0.0;
-  if (E >= energies.front()) {
-    i = lower_bound_index(energies.begin(), energies.end(), E);
-    if (i + 1 < energies.size())
-      f = (E - energies[i]) / (energies[i + 1] - energies[i]);
+  if (n < 2 || E < energies.front())
+    return;
+
+  i = lower_bound_index(energies.begin(), energies.end(), E);
+  if (i < n - 1) {
+    f = (E - energies[i]) / (energies[i + 1] - energies[i]);
+  } else {
+    // E lies above the top of the grid; use the topmost interval
+    i = n - 2;
+    f = 1.0;
   }
 }