random/distribution/InverseSquareDistribution.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file celeritas/random/distribution/InverseSquareDistribution.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/Macros.hh"
0012 #include "corecel/Types.hh"
0013
0014 #include "UniformRealDistribution.hh"
0015
0016 namespace celeritas
0017 {
0018 //---------------------------------------------------------------------------//
0019 /*!
0020  * Sample \f$ 1/x^2 \f$ over a given domain.
0021  *
0022  * This distribution is defined on a positive range \f$ [a, b) \f$ and has the
0023  * normalized PDF:
0024  * \f[
0025    f(x; a, b) = \frac{a b}{x^2 (b - a)} \quad \mathrm{for} \  a \le x < b
0026    \f]
0027  * which integrated into a CDF and inverted gives a sample:
0028  * \f[
0029   x = \frac{a b}{a + \xi (b - a)}
0030    \f]
0031  */
0032 template<class RealType = ::celeritas::real_type>
0033 class InverseSquareDistribution
0034 {
0035   public:
0036     //!@{
0037     //! \name Type aliases
0038     using real_type = RealType;
0039     using result_type = real_type;
0040     //!@}
0041
0042   public:
0043     // Construct on an interval
0044     inline CELER_FUNCTION InverseSquareDistribution(real_type a, real_type b);
0045
0046     // Sample a random number according to the distribution
0047     template<class Generator>
0048     inline CELER_FUNCTION result_type operator()(Generator& rng) const;
0049
0050   private:
0051     RealType product_;
0052     UniformRealDistribution<RealType> sample_denom_;
0053 };
0054
0055 //---------------------------------------------------------------------------//
0056 // INLINE DEFINITIONS
0057 //---------------------------------------------------------------------------//
0058 /*!
0059  * Construct on the interval [a, b).
0060  *
0061  * As with UniformRealDistribution, it is allowable for the two bounds to be
0062  * out of order.
0063  */
0064 template<class RealType>
0065 CELER_FUNCTION
0066 InverseSquareDistribution<RealType>::InverseSquareDistribution(real_type a,
0067                                                                real_type b)
0068     : product_{a * b}, sample_denom_{a, b}
0069 {
0070     CELER_EXPECT(a > 0);
0071     CELER_EXPECT(b >= a);
0072 }
0073
0074 //---------------------------------------------------------------------------//
0075 /*!
0076  * Sample a random number according to the distribution.
0077  */
0078 template<class RealType>
0079 template<class Generator>
0080 CELER_FUNCTION auto
0081 InverseSquareDistribution<RealType>::operator()(Generator& rng) const -> result_type
0082 {
0083     return product_ / sample_denom_(rng);
0084 }
0085
0086 //---------------------------------------------------------------------------//
0087 }  // namespace celeritas