Addressing some comments from @paulromano

include/openmc/constants.h

@@ -319,8 +319,10 @@ enum class Interpolation {
   lin_log = 3,
   log_lin = 4,
   log_log = 5,
-  quadratic = 6,
-  cubic = 7
+  // skip 6 b/c ENDF-6 reserves this value for
+  // "special one-dimensional interpolation law"
+  quadratic = 7,
+  cubic = 8
 };
 
 enum class RunMode {

include/openmc/interpolate.h

@@ -33,10 +33,12 @@ inline double interpolate_log_log(
   return y0 * std::exp(f * std::log(y1 / y0));
 }
 
-inline double interpolate_lagrangian(const std::vector<double>& xs,
-  const std::vector<double>& ys, int idx, double x, int order)
+inline double interpolate_lagrangian(gsl::span<const double> xs,
+  gsl::span<const double> ys, int idx, double x, int order)
 {
-  std::vector<double> coeffs;
+  Expects(order <= 3);
+  std::array<double, 4> coeffs;
+  coeffs.fill(0.0);
 
   for (int i = 0; i < order + 1; i++) {
     double numerator {1.0};
@@ -47,14 +49,14 @@ inline double interpolate_lagrangian(const std::vector<double>& xs,
       numerator *= (x - xs[idx + j]);
       denominator *= (xs[idx + i] - xs[idx + j]);
     }
-    coeffs.push_back(numerator / denominator);
+    coeffs[i] = numerator / denominator;
   }
 
   return std::inner_product(
     coeffs.begin(), coeffs.end(), ys.begin() + idx, 0.0);
 }
 
-double interpolate(const std::vector<double>& xs, const std::vector<double>& ys,
+double interpolate(gsl::span<const double> xs, gsl::span<const double> ys,
   double x, Interpolation i = Interpolation::lin_lin)
 {
   int idx = lower_bound_index(xs.begin(), xs.end(), x);
@@ -63,6 +65,8 @@ double interpolate(const std::vector<double>& xs, const std::vector<double>& ys,
     idx--;
 
   switch (i) {
+  case Interpolation::histogram:
+    return ys[idx];
   case Interpolation::lin_lin:
     return interpolate_lin_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
   case Interpolation::log_log:

openmc/filter.py

@@ -1887,7 +1887,9 @@ class EnergyFunctionFilter(Filter):
     y : iterable of Real
         A grid of interpolant values in [eV]
     interpolation : str
-        The type of interpolation to be used.
+        The type of interpolation to be used. One of
+        ('histogram', 'linear-linear', 'linear-log', 'log-linear',
+        'log-log', 'quadratic', 'cubic')
     id : int
         Unique identifier for the filter
     num_bins : Integral
@@ -1896,9 +1898,12 @@ class EnergyFunctionFilter(Filter):
     """
 
     # keys selected to match those in function.py where possible
-    INTERPOLATION_SCHEMES = {2: 'linear-linear', 3: 'linear-log',
-                             4: 'log-linear', 5: 'log-log',
-                             6 : 'quadratic', 7 : 'cubic'}
+    # skip 6 b/c ENDF-6 reserves this value for
+    # "special one-dimensional interpolation law"
+    INTERPOLATION_SCHEMES = {1: 'histogram', 2: 'linear-linear',
+                             3: 'linear-log', 4: 'log-linear',
+                             5: 'log-log', 7: 'quadratic',
+                             8: 'cubic'}
 
     def __init__(self, energy, y, interpolation='linear-linear', filter_id=None):
         self.energy = energy
@@ -1991,22 +1996,11 @@ class EnergyFunctionFilter(Filter):
         if tab1d.n_regions > 1:
             raise ValueError('Only Tabulated1Ds with a single interpolation '
                              'region are supported')
-        if tab1d.interpolation[0] not in (2, 3, 4, 5):
-            raise ValueError('Only linear-linear, linear-log, log-linear, and '
+        interpolation = tab1d.interpolation[0]
+        if interpolation not in cls.INTERPOLATION_SCHEMES.keys():
+            raise ValueError('Only histogram, linear-linear, linear-log, log-linear, and '
                              'log-log Tabulated1Ds are supported')
-        out = cls(tab1d.x, tab1d.y)
-
-        # set interpolation type
-        if tab1d.interpolation[0] == 2:
-            out.interpolation = 'linear-linear'
-        elif tab1d.interpolation[0] == 3:
-            out.interpolation = 'linear-log'
-        elif tab1d.interpolation[0] == 4:
-            out.interpolation = 'log-linear'
-        elif tab1d.interpolation[0] == 5:
-            out.interpolation = 'log-log'
-
-        return out
+        return cls(tab1d.x, tab1d.y, interpolation)
 
     @property
     def energy(self):

src/tallies/filter_energyfunc.cpp

@@ -30,7 +30,9 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node)
   interpolation_ = Interpolation::lin_lin;
   if (check_for_node(node, "interpolation")) {
     std::string interpolation = get_node_value(node, "interpolation");
-    if (interpolation == "linear-linear") {
+    if (interpolation == "histogram") {
+      interpolation_ = Interpolation::histogram;
+    } else if (interpolation == "linear-linear") {
       interpolation_ = Interpolation::lin_lin;
     } else if (interpolation == "linear-log") {
       interpolation_ = Interpolation::lin_log;

tests/regression_tests/filter_energyfun/test.py

@@ -31,7 +31,7 @@ def model():
     assert filt1.interpolation == 'linear-linear'
 
     with pytest.raises(ValueError):
-        filt1.interpolation = '5th order polynomial'
+        filt1.interpolation = '🥏'
 
     # Also make a filter with the .from_tabulated1d constructor.  Make sure
     # the filters are identical.