corecel/grid/SplineInterpolator.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file corecel/grid/SplineInterpolator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Assert.hh"
0012 #include "corecel/Macros.hh"
0013 #include "corecel/Types.hh"
0014 #include "corecel/cont/Array.hh"
0015 #include "corecel/math/Algorithms.hh"
0016 #include "corecel/math/PolyEvaluator.hh"
0017
0018 #include "detail/InterpolatorTraits.hh"
0019
0020 namespace celeritas
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Interpolate using a cubic spline.
0025  *
0026  * Given a set of \f$ n \f$ data points \f$ (x_i, y_i) \f$ such that \f$ x_0 <
0027  * x_1 < \dots < x_{n - 1} \f$, a cubic spline \f$ S(x) \f$ interpolating on
0028  * the points is a piecewise polynomial function consisting of \f$ n - 1 \f$
0029  * cubic polynomials \f$ S_i \f$ defined on \f$ [x_i, x_{i + 1}] \f$. The \f$
0030  * S_i \f$ are joined at \f$ x_i \f$ such that both the first and second
0031  * derivatives, \f$ S'_i \f$ and \f$ S''_i \f$, are continuous.
0032  *
0033  * The \f$ i^{\text{th}} \f$ piecewise polynomial \f$ S_i \f$ is given by:
0034  * \f[
0035    S_i(x) = a_0 + a_1(x - x_i) + a_2(x - x_i)^2 + a_3(x - x_i)^3,
0036  * \f]
0037  * where \f$ a_i \f$ are the polynomial coefficients, expressed in terms of the
0038  * second derivatives as:
0039  * \f{align}{
0040    a_0 &= y_i \\
0041    a_1 &= \frac{\Delta y_i}{\Delta x_i} - \frac{\Delta x_i}{6} \left[ S''_{i +
0042           1} + 2 S''_{i} \right] \\
0043    a_2 &= \frac{S''_i}{2} \\
0044    a_3 &= \frac{1}{6 \Delta x_i} \left[ S''_{i + 1} - S''_i \right]
0045  * \f}
0046  */
0047 template<typename T = ::celeritas::real_type>
0048 class SplineInterpolator
0049 {
0050   public:
0051     //!@{
0052     //! \name Type aliases
0053     using real_type = T;
0054     //!@}
0055
0056     //! (x, y) point and second derivative
0057     struct Spline
0058     {
0059         real_type x;
0060         real_type y;
0061         real_type ddy;
0062     };
0063
0064   public:
0065     // Construct with left and right values for x, y, and the second derivative
0066     inline CELER_FUNCTION SplineInterpolator(Spline left, Spline right);
0067
0068     // Interpolate
0069     inline CELER_FUNCTION real_type operator()(real_type x) const;
0070
0071   private:
0072     real_type x_lower_;  //!< Lower grid point \f$ x_i \f$
0073     Array<T, 4> a_;  //!< Spline coefficients
0074 };
0075
0076 //---------------------------------------------------------------------------//
0077 // INLINE DEFINITIONS
0078 //---------------------------------------------------------------------------//
0079 /*!
0080  * Construct with left and right values for x, y, and the second derivative.
0081  */
0082 template<class T>
0083 CELER_FUNCTION
0084 SplineInterpolator<T>::SplineInterpolator(Spline left, Spline right)
0085 {
0086     CELER_EXPECT(left.x < right.x);
0087
0088     x_lower_ = left.x;
0089     real_type h = right.x - left.x;
0090     a_[0] = left.y;
0091     a_[1] = (right.y - left.y) / h - h / 6 * (right.ddy + 2 * left.ddy);
0092     a_[2] = T(0.5) * left.ddy;
0093     a_[3] = 1 / (6 * h) * (right.ddy - left.ddy);
0094 }
0095
0096 //---------------------------------------------------------------------------//
0097 /*!
0098  * Interpolate using a cubic spline.
0099  */
0100 template<class T>
0101 CELER_FUNCTION auto SplineInterpolator<T>::operator()(real_type x) const
0102     -> real_type
0103 {
0104     return PolyEvaluator{a_}(x - x_lower_);
0105 }
0106
0107 //---------------------------------------------------------------------------//
0108 }  // namespace celeritas