From 3f06a42abb89a9b09cb3f2e86680a31bae90cc08 Mon Sep 17 00:00:00 2001
From: GuySten <62616591+GuySten@users.noreply.github.com>
Date: Tue, 23 Dec 2025 22:18:59 +0200
Subject: [PATCH] Refactor get_energy_index to prevent repetition (#3686)

---
 include/openmc/math_functions.h | 11 +++++++++++
 src/distribution_angle.cpp      | 18 +++---------------
 src/math_functions.cpp          | 15 +++++++++++++++
 src/secondary_correlated.cpp    | 16 +++-------------
 src/secondary_kalbach.cpp       | 16 +++-------------
 src/secondary_thermal.cpp       | 15 +--------------
 6 files changed, 36 insertions(+), 55 deletions(-)

diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index c30ef7558..0d960c33d 100644
--- a/include/openmc/math_functions.h
+++ b/include/openmc/math_functions.h
@@ -9,6 +9,7 @@
 #include <cstdlib>
 
 #include "openmc/position.h"
+#include "openmc/search.h"
 
 namespace openmc {
 
@@ -200,5 +201,15 @@ std::complex<double> faddeeva(std::complex<double> z);
 //! \return Derivative of Faddeeva function evaluated at z
 std::complex<double> w_derivative(std::complex<double> z, int order);
 
+//! Helper function to get index and interpolation function on an incident
+//! energy grid
+//!
+//! \param energies energy grid
+//! \param E incident energy
+//! \param i grid index
+//! \param f interpolation factor
+void get_energy_index(
+  const vector<double>& energies, double E, int& i, double& f);
+
 } // namespace openmc
 #endif // OPENMC_MATH_FUNCTIONS_H
diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp
index 5e92b794a..a571b2a5f 100644
--- a/src/distribution_angle.cpp
+++ b/src/distribution_angle.cpp
@@ -7,6 +7,7 @@
 
 #include "openmc/endf.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/math_functions.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 #include "openmc/vector.h" // for vector
@@ -64,23 +65,10 @@ AngleDistribution::AngleDistribution(hid_t group)
 
 double AngleDistribution::sample(double E, uint64_t* seed) const
 {
-  // Determine number of incoming energies
-  auto n = energy_.size();
-
-  // Find energy bin and calculate interpolation factor -- if the energy is
-  // outside the range of the tabulated energies, choose the first or last bins
+  // Find energy bin and calculate interpolation factor
   int i;
   double r;
-  if (E < energy_[0]) {
-    i = 0;
-    r = 0.0;
-  } else if (E > energy_[n - 1]) {
-    i = n - 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
   if (r > prn(seed))
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 5469b56c8..9473f928b 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -919,4 +919,19 @@ std::complex<double> w_derivative(std::complex<double> z, int order)
   }
 }
 
+// Helper function to get index and interpolation function on an incident energy
+// grid
+void get_energy_index(
+  const vector<double>& energies, double E, int& i, double& f)
+{
+  // Get index and interpolation factor for linear-linear energy grid
+  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]);
+  }
+}
+
 } // namespace openmc
diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp
index 0e4891dd3..3f26fdb01 100644
--- a/src/secondary_correlated.cpp
+++ b/src/secondary_correlated.cpp
@@ -10,6 +10,7 @@
 
 #include "openmc/endf.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/math_functions.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 
@@ -156,21 +157,10 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group)
 void CorrelatedAngleEnergy::sample(
   double E_in, double& E_out, double& mu, 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();
+  // Find energy bin and calculate interpolation factor
   int i;
   double r;
-  if (E_in < energy_[0]) {
-    i = 0;
-    r = 0.0;
-  } else if (E_in > energy_[n_energy_in - 1]) {
-    i = n_energy_in - 2;
-    r = 1.0;
-  } else {
-    i = lower_bound_index(energy_.begin(), energy_.end(), E_in);
-    r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]);
-  }
+  get_energy_index(energy_, E_in, i, r);
 
   // Sample between the ith and [i+1]th bin
   int l = r > prn(seed) ? i + 1 : i;
diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp
index 4a7878343..6ac91e665 100644
--- a/src/secondary_kalbach.cpp
+++ b/src/secondary_kalbach.cpp
@@ -9,6 +9,7 @@
 #include "xtensor/xview.hpp"
 
 #include "openmc/hdf5_interface.h"
+#include "openmc/math_functions.h"
 #include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
@@ -117,21 +118,10 @@ KalbachMann::KalbachMann(hid_t group)
 void KalbachMann::sample(
   double E_in, double& E_out, double& mu, 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();
+  // Find energy bin and calculate interpolation factor
   int i;
   double r;
-  if (E_in < energy_[0]) {
-    i = 0;
-    r = 0.0;
-  } else if (E_in > energy_[n_energy_in - 1]) {
-    i = n_energy_in - 2;
-    r = 1.0;
-  } else {
-    i = lower_bound_index(energy_.begin(), energy_.end(), E_in);
-    r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]);
-  }
+  get_energy_index(energy_, E_in, i, r);
 
   // Sample between the ith and [i+1]th bin
   int l = r > prn(seed) ? i + 1 : i;
diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp
index 8b9e8737c..030d398aa 100644
--- a/src/secondary_thermal.cpp
+++ b/src/secondary_thermal.cpp
@@ -1,6 +1,7 @@
 #include "openmc/secondary_thermal.h"
 
 #include "openmc/hdf5_interface.h"
+#include "openmc/math_functions.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 
@@ -11,20 +12,6 @@
 
 namespace openmc {
 
-// Helper function to get index on incident energy grid
-void get_energy_index(
-  const vector<double>& energies, double E, int& i, double& f)
-{
-  // Get index and interpolation factor for elastic grid
-  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]);
-  }
-}
-
 //==============================================================================
 // CoherentElasticAE implementation
 //==============================================================================