celeritas/grid/InverseRangeCalculator.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 celeritas/grid/InverseRangeCalculator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Assert.hh"
0012 #include "corecel/data/Collection.hh"
0013 #include "corecel/grid/Interpolator.hh"
0014 #include "corecel/grid/NonuniformGrid.hh"
0015 #include "corecel/grid/SplineInterpolator.hh"
0016 #include "corecel/grid/UniformGrid.hh"
0017 #include "corecel/grid/UniformGridData.hh"
0018 #include "corecel/math/Algorithms.hh"
0019 #include "corecel/math/Quantity.hh"
0020 #include "celeritas/Quantities.hh"
0021
0022 namespace celeritas
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Calculate the energy that would limit a particle to a particular range.
0027  *
0028  * This should provide the inverse of the result of \c RangeCalculator. The
0029  * given \c range is not allowed to be greater than the maximum range in the
0030  * physics data.
0031  *
0032  * The range must be monotonically increasing in energy, since it's defined as
0033  * the integral of the inverse of the stopping power (which is always
0034  * positive). For ranges shorter than the minimum energy in the table, the
0035  * resulting energy is scaled:
0036  * \f[
0037     E = E_\mathrm{min} \left( \frac{r}{r_\mathrm{min}} \right)^2
0038  * \f]
0039  * This scaling is the inverse of the off-the-end energy scaling in the
0040  * RangeCalculator.
0041  *
0042  * \todo Construct with \c UniformGridRecord
0043  */
0044 class InverseRangeCalculator
0045 {
0046   public:
0047     //!@{
0048     //! \name Type aliases
0049     using Energy = units::MevEnergy;
0050     using Values
0051         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0052     //!@}
0053
0054   public:
0055     // Construct from state-independent data
0056     inline CELER_FUNCTION InverseRangeCalculator(UniformGridRecord const& grid,
0057                                                  Values const& values);
0058
0059     // Find and interpolate from the energy
0060     inline CELER_FUNCTION Energy operator()(real_type range) const;
0061
0062   private:
0063     UniformGrid log_energy_;
0064     NonuniformGrid<real_type> range_;
0065     Span<real_type const> deriv_;
0066 };
0067
0068 //---------------------------------------------------------------------------//
0069 // INLINE DEFINITIONS
0070 //---------------------------------------------------------------------------//
0071 /*!
0072  * Construct from range data.
0073  *
0074  * The range is expected to be monotonically increasing with energy.
0075  * Lower-energy particles have shorter ranges.
0076  */
0077 CELER_FUNCTION
0078 InverseRangeCalculator::InverseRangeCalculator(UniformGridRecord const& grid,
0079                                                Values const& values)
0080     : log_energy_(grid.grid)
0081     , range_(grid.value, values)
0082     , deriv_(values[grid.derivative])
0083 {
0084     CELER_EXPECT(range_.size() == log_energy_.size());
0085 }
0086
0087 //---------------------------------------------------------------------------//
0088 /*!
0089  * Calculate the energy of a particle that has the given range.
0090  */
0091 CELER_FUNCTION auto InverseRangeCalculator::operator()(real_type range) const
0092     -> Energy
0093 {
0094     CELER_EXPECT(range >= 0 && range <= range_.back());
0095
0096     if (range < range_.front())
0097     {
0098         // Very short range:  this corresponds to "energy < emin" for range
0099         // calculation: range = r[0] * sqrt(E / E[0])
0100         return Energy{std::exp(log_energy_.front())
0101                       * ipow<2>(range / range_.front())};
0102     }
0103     // Range should *never* exceed the longest range (highest energy) since
0104     // that should have limited the step
0105     if (CELER_UNLIKELY(range >= range_.back()))
0106     {
0107         CELER_ASSERT(range == range_.back());
0108         return Energy{std::exp(log_energy_.back())};
0109     }
0110
0111     // Search for lower bin index
0112     auto idx = range_.find(range);
0113     CELER_ASSERT(idx + 1 < log_energy_.size());
0114
0115     real_type result;
0116     if (deriv_.empty())
0117     {
0118         // Interpolate: 'x' = range, y = log energy
0119         result = LinearInterpolator<real_type>(
0120             {range_[idx], std::exp(log_energy_[idx])},
0121             {range_[idx + 1], std::exp(log_energy_[idx + 1])})(range);
0122     }
0123     else
0124     {
0125         // Use cubic spline interpolation
0126         result = SplineInterpolator<real_type>(
0127             {range_[idx], std::exp(log_energy_[idx]), deriv_[idx]},
0128             {range_[idx + 1], std::exp(log_energy_[idx + 1]), deriv_[idx + 1]})(
0129             range);
0130     }
0131     return Energy{result};
0132 }
0133
0134 //---------------------------------------------------------------------------//
0135 }  // namespace celeritas