From d18e9e524b1b5293c5d3ec4aa39982ff818cd714 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 1 Sep 2026 22:16:50 +0300
Subject: [PATCH 1/2] simplify get-energy-index

---
 src/distribution_angle.cpp  | 10 ++++------
 src/distribution_energy.cpp | 13 +------------
 src/math_functions.cpp      | 20 +++++++++++++++-----
 3 files changed, 20 insertions(+), 23 deletions(-)

diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp
index ecb5961f632..bf2f2b54f28 100644
--- a/src/distribution_angle.cpp
+++ b/src/distribution_angle.cpp
@@ -89,12 +89,10 @@ double AngleDistribution::evaluate(double E, double mu) const
   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;
+  // Both distributions are always valid since get_energy_index keeps i within
+  // the topmost interval
+  return r * distribution_[i + 1]->evaluate(mu) +
+         (1.0 - r) * distribution_[i]->evaluate(mu);
 }
 
 } // namespace openmc
diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index 7712c0d763f..4e124e42f39 100644
--- a/src/distribution_energy.cpp
+++ b/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;
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index ddacc2bd9b6..2a7e434d7b7 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -991,13 +991,23 @@ 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, matching the handling in
+  // ContinuousTabular::sample.
+  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;
   }
 }
 

From f6f95afb7d5b83c4ec38d9842e6e1529e22d22e1 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 1 Sep 2026 22:23:10 +0300
Subject: [PATCH 2/2] simplify

---
 src/distribution_angle.cpp | 3 ---
 src/math_functions.cpp     | 3 +--
 2 files changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp
index bf2f2b54f28..f3efa27f507 100644
--- a/src/distribution_angle.cpp
+++ b/src/distribution_angle.cpp
@@ -88,9 +88,6 @@ double AngleDistribution::evaluate(double E, double mu) const
   int i;
   double r;
   get_energy_index(energy_, E, i, r);
-
-  // Both distributions are always valid since get_energy_index keeps i within
-  // the topmost interval
   return r * distribution_[i + 1]->evaluate(mu) +
          (1.0 - r) * distribution_[i]->evaluate(mu);
 }
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 2a7e434d7b7..6be79fc0fd8 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -993,8 +993,7 @@ void get_energy_index(
 {
   // 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, matching the handling in
-  // ContinuousTabular::sample.
+  // energies[i + 1] are valid for callers.
   const int n = energies.size();
   i = 0;
   f = 0.0;