From 0d551cbb21f76894ea70f9e712f81d09f4b3c7f2 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Thu, 28 Jul 2022 11:14:30 -0500
Subject: [PATCH] Adding support for log-log interpolation on the C++ side

---
 include/openmc/tallies/filter_energyfunc.h |  2 ++
 openmc/filter.py                           |  5 +++--
 src/tallies/filter_energyfunc.cpp          | 26 ++++++++++++++++++----
 3 files changed, 27 insertions(+), 6 deletions(-)

diff --git a/include/openmc/tallies/filter_energyfunc.h b/include/openmc/tallies/filter_energyfunc.h
index c066dfa17..f650976ef 100644
--- a/include/openmc/tallies/filter_energyfunc.h
+++ b/include/openmc/tallies/filter_energyfunc.h
@@ -1,6 +1,7 @@
 #ifndef OPENMC_TALLIES_FILTER_ENERGYFUNC_H
 #define OPENMC_TALLIES_FILTER_ENERGYFUNC_H
 
+#include "openmc/constants.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
 
@@ -39,6 +40,7 @@ public:
 
   const vector<double>& energy() const { return energy_; }
   const vector<double>& y() const { return y_; }
+  Interpolation interpolation_;
   void set_data(gsl::span<const double> energy, gsl::span<const double> y);
 
 private:
diff --git a/openmc/filter.py b/openmc/filter.py
index 0d8458be5..f44ac99f8 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -1957,8 +1957,9 @@ class EnergyFunctionFilter(Filter):
         filter_id = int(group.name.split('/')[-1].lstrip('filter '))
 
         out = cls(energy, y, filter_id=filter_id)
-        if 'interpolation' in group:
-            out.interpolation = group['interpolation'][()]
+        if 'interpolation' in group.attrs:
+            out.interpolation =  \
+                openmc.data.INTERPOLATION_SCHEME[group.attrs['interpolation']]
 
         return out
 
diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp
index fd595e332..2fcca7110 100644
--- a/src/tallies/filter_energyfunc.cpp
+++ b/src/tallies/filter_energyfunc.cpp
@@ -22,9 +22,15 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node)
 
   if (!check_for_node(node, "y"))
     fatal_error("y values not specified for EnergyFunction filter.");
-
+  
   auto y = get_node_array<double>(node, "y");
 
+  // use linear-linear interpolation by default
+  interpolation_ = Interpolation::lin_lin;
+  if (check_for_node(node, "interpolation")) {
+    std::string interpolation = get_node_value(node, "interpolation");
+  }
+
   this->set_data(energy, y);
 }
 
@@ -58,12 +64,23 @@ void EnergyFunctionFilter::get_all_bins(
     // Search for the incoming energy bin.
     auto i = lower_bound_index(energy_.begin(), energy_.end(), p.E_last());
 
-    // Compute the interpolation factor between the nearest bins.
-    double f = (p.E_last() - energy_[i]) / (energy_[i + 1] - energy_[i]);
+    double f, w;
+    switch (interpolation_) {
+      case Interpolation::lin_lin:
+        f = (p.E_last() - energy_[i]) / (energy_[i + 1] - energy_[i]);
+        w = (1 - f) * y_[i] + f * y_[i + 1];
+        break;
+      case Interpolation::log_log:
+        f = log(p.E_last() / energy_[i]) / log(energy_[i + 1] / energy_[i]);
+        w = y_[i] * exp(f * log(y_[i + 1] / y_[i]));
+        break;
+      default:
+        fatal_error(fmt::format("Invalid interpolation scheme found on EnergyFunctionFilter {}", id()));
+    }
 
     // Interpolate on the lin-lin grid.
     match.bins_.push_back(0);
-    match.weights_.push_back((1 - f) * y_[i] + f * y_[i + 1]);
+    match.weights_.push_back(w);
   }
 }
 
@@ -72,6 +89,7 @@ void EnergyFunctionFilter::to_statepoint(hid_t filter_group) const
   Filter::to_statepoint(filter_group);
   write_dataset(filter_group, "energy", energy_);
   write_dataset(filter_group, "y", y_);
+  write_attr_int(filter_group, "interpolation", interpolation_);
 }
 
 std::string EnergyFunctionFilter::text_label(int bin) const