corecel/grid/FindInterp.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/FindInterp.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Assert.hh"
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012
0013 namespace celeritas
0014 {
0015 //---------------------------------------------------------------------------//
0016 /*!
0017  * Result of finding a point on a grid for interpolating.
0018  *
0019  * The resulting index will be in [0, grid.size() - 1)
0020  * and the fraction will be in [0, 1).
0021  */
0022 template<class T>
0023 struct FindInterp
0024 {
0025     size_type index{};  //!< Lower index into the grid
0026     T fraction{};  //!< Fraction of the value between its neighbors
0027 };
0028
0029 template<class T>
0030 CELER_FUNCTION FindInterp(size_type, T) -> FindInterp<T>;
0031
0032 //---------------------------------------------------------------------------//
0033 /*!
0034  * Find the index of the value and its fraction between neighboring points.
0035  *
0036  * The grid class should have a floating point value and must have methods \c
0037  * find, \c front, \c back, and \c operator[] .
0038  *
0039  * The value must be bounded by the grid and less than the final value. The
0040  * result will always have an index such that its neighbor to the right is a
0041  * valid point on the grid, and the fraction between neghbors may be zero (in
0042  * the case where the value is exactly on a grid point) but is always less than
0043  * one. If the requested point is exactly on a coincident grid point, the lower
0044  * point and a fraction of zero will result.
0045  */
0046 template<class Grid>
0047 inline CELER_FUNCTION auto
0048 find_interp(Grid const& grid, typename Grid::value_type value)
0049 {
0050     CELER_EXPECT(value >= grid.front() && value < grid.back());
0051
0052     auto index = grid.find(value);
0053     CELER_ASSERT(index + 1 < grid.size());
0054     auto const lower_val = grid[index];
0055     auto const upper_val = grid[index + 1];
0056     return FindInterp{index, (value - lower_val) / (upper_val - lower_val)};
0057 }
0058
0059 //---------------------------------------------------------------------------//
0060 }  // namespace celeritas