optical/surface/GaussianRoughnessSampler.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/optical/surface/GaussianRoughnessSampler.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/math/Algorithms.hh"
0012 #include "corecel/random/distribution/NormalDistribution.hh"
0013 #include "corecel/random/distribution/RejectionSampler.hh"
0014 #include "corecel/random/distribution/TruncatedDistribution.hh"
0015 #include "celeritas/phys/InteractionUtils.hh"
0016
0017 namespace celeritas
0018 {
0019 namespace optical
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Sample a facet normal from a Gaussian roughness model.
0024  *
0025  * The Gaussian roughness model was introduced in
0026  * \citet{levin-morephysical-1996, https://doi.org/10.1109/NSSMIC.1996.591410}
0027  * . The "facet slope", an angle \f$ \alpha \f$ along a linear slice of a
0028  * crystal surface, is approximated as a normal distribution standard deviation
0029  * \f$ \sigma_\alpha \f$ .  (The paper justifies this distribution based on
0030  * surface roughness measurements with a bismuth germanate [BGO] crystal.)
0031  * Assuming an azimuthally isotropic surface, the polar distribution must be
0032  * expressed in terms of the tilt \f$ \theta \f$. The Jacobian factor for
0033  * spherical coordinates contributes a \f$ \sin \theta \f$ term, leading to the
0034  * spherical PDF \f[
0035    p(\sigma_\alpha; \theta,\phi) = \frac{1}{2\pi}\,\frac{1}{\mathcal N}\,
0036    \exp\!\left(-\frac{\theta^{2}}{2\sigma_\alpha^{2}}\right)\sin\theta \,.
0037  * \f]
0038  *
0039  * The polar angle \f$ \theta \f$ is sampled using rejection:
0040  * - Draw \f$ \alpha \f$ from the positive half of a normal distribution
0041  * - Reject angles greater than 90 degrees (for physicality) or 4 sigma (for
0042  *   sampling efficiency)
0043  * - Use an acceptance function \f$ \sin \theta \f$ bounded by the maximum
0044  *   theta
0045  *
0046  * The extra limitation of the angle being less than 4 sigma reduces the
0047  * rejection fraction by a factor of ~25 for smooth crystals (sigma = 0.01).
0048  */
0049 class GaussianRoughnessSampler
0050 {
0051   public:
0052     // Construct from sigma_alpha, global normal, and incident direction
0053     inline CELER_FUNCTION GaussianRoughnessSampler(Real3 const& normal,
0054                                                    real_type sigma_alpha);
0055
0056     // Sample facet normal
0057     template<class Engine>
0058     inline CELER_FUNCTION Real3 operator()(Engine& rng);
0059
0060   private:
0061     Real3 normal_;
0062     TruncatedDistribution<NormalDistribution<real_type>> sample_alpha_;
0063     real_type f_max_;
0064     RejectionSampler<real_type> reject_;
0065 };
0066
0067 //---------------------------------------------------------------------------//
0068 // INLINE DEFINITIONS
0069 //---------------------------------------------------------------------------//
0070 /*!
0071  * Construct from sigma_alpha and global normal.
0072  */
0073 CELER_FUNCTION GaussianRoughnessSampler::GaussianRoughnessSampler(
0074     Real3 const& normal, real_type sigma_alpha)
0075     : normal_(normal)
0076     , sample_alpha_(-constants::pi / 2, constants::pi / 2, 0, sigma_alpha)
0077     , f_max_(fmin(real_type{1}, 4 * sigma_alpha))
0078     , reject_(f_max_)
0079 {
0080     CELER_EXPECT(sigma_alpha > 0);
0081     CELER_EXPECT(is_soft_unit_vector(normal_));
0082 }
0083
0084 //---------------------------------------------------------------------------//
0085 /*!
0086  * Sample a facet normal via the Gaussian roughness model.
0087  */
0088 template<class Engine>
0089 CELER_FUNCTION Real3 GaussianRoughnessSampler::operator()(Engine& rng)
0090 {
0091     using std::fabs;
0092
0093     real_type cos_alpha{};
0094     real_type sin_alpha{};
0095     do
0096     {
0097         // Sample positive angle according to gaussian (chances of having a
0098         // nonpositive slope are generally vanishingly small)
0099         real_type alpha = std::fabs(sample_alpha_(rng));
0100         sincos(alpha, &sin_alpha, &cos_alpha);
0101
0102         // Transform to polar angle using rejection
0103     } while (sin_alpha < f_max_ && reject_(sin_alpha, rng));
0104
0105     // Rotate normal by alpha and then sample azimuthal rotation uniformly
0106     return ExitingDirectionSampler{cos_alpha, normal_}(rng);
0107 }
0108
0109 //---------------------------------------------------------------------------//
0110 }  // namespace optical
0111 }  // namespace celeritas