diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h index b88ab2f2e7..0683140339 100644 --- a/include/openmc/interpolate.h +++ b/include/openmc/interpolate.h @@ -60,6 +60,31 @@ inline double interpolate(const std::vector& xs, return interpolate(xs, ys, idx, x, i); } +inline double interpolate_lagrangian( + const std::vector& xs, const std::vector& ys, double x) +{ + int idx = lower_bound_index(xs.begin(), xs.end(), x); + + std::vector coeffs; + + int order = 3; + + for (int i = 0; i < order + 1; i++) { + double numerator {1.0}; + double denominator {1.0}; + for (int j = 0; j < order; j++) { + if (i == j) + continue; + numerator *= (x - xs[idx + j]); + denominator *= (xs[idx + i] - xs[idx + j]); + } + coeffs.push_back(numerator / denominator); + } + + return std::inner_product( + coeffs.begin(), coeffs.end(), ys.begin() + idx, 0.0); +} + } // namespace openmc #endif \ No newline at end of file