optical/interactor/RayleighInteractor.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/interactor/RayleighInteractor.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <type_traits>
0010
0011 #include "corecel/Assert.hh"
0012 #include "corecel/math/ArrayOperators.hh"
0013 #include "corecel/math/ArraySoftUnit.hh"
0014 #include "corecel/math/ArrayUtils.hh"
0015 #include "corecel/random/distribution/BernoulliDistribution.hh"
0016 #include "corecel/random/distribution/IsotropicDistribution.hh"
0017 #include "corecel/random/distribution/RejectionSampler.hh"
0018 #include "corecel/random/engine/RanluxppRngEngine.hh"
0019 #include "celeritas/optical/Interaction.hh"
0020 #include "celeritas/optical/ParticleTrackView.hh"
0021
0022 namespace celeritas
0023 {
0024 namespace optical
0025 {
0026 //---------------------------------------------------------------------------//
0027 /*!
0028  * Sample optical Rayleigh scattering.
0029  *
0030  * Optical Rayleigh scattering is the elastic scattering of optical photons
0031  * in a material. The scattered polarization is guaranteed to be in the
0032  * same plane as the original polarization and new direction.
0033  */
0034 class RayleighInteractor
0035 {
0036   public:
0037     // Construct interactor from an optical track
0038     inline CELER_FUNCTION RayleighInteractor(ParticleTrackView const& particle,
0039                                              Real3 const& direction);
0040
0041     // Sample an interaction with the given RNG
0042     template<class Engine>
0043     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0044
0045   private:
0046     Real3 const& inc_dir_;  //!< Direction of incident photon
0047     Real3 const& inc_pol_;  //!< Polarization of incident photon
0048     IsotropicDistribution<real_type> sample_direction_{};
0049     RejectionSampler<real_type> reject_angle_{1};
0050 };
0051
0052 //---------------------------------------------------------------------------//
0053 // INLINE DEFINITIONS
0054 //---------------------------------------------------------------------------//
0055 /*!
0056  * Construct the interactor for the given optical track.
0057  */
0058 CELER_FUNCTION RayleighInteractor::RayleighInteractor(
0059     ParticleTrackView const& particle, Real3 const& direction)
0060     : inc_dir_(direction), inc_pol_(particle.polarization())
0061 {
0062     CELER_EXPECT(is_soft_unit_vector(inc_dir_));
0063     CELER_EXPECT(is_soft_unit_vector(inc_pol_));
0064     CELER_EXPECT(is_soft_orthogonal(inc_dir_, inc_pol_));
0065 }
0066
0067 //---------------------------------------------------------------------------//
0068 /*!
0069  * Sample a single optical Rayleigh interaction.
0070  */
0071 template<class Engine>
0072 CELER_FUNCTION Interaction RayleighInteractor::operator()(Engine& rng)
0073 {
0074     if constexpr (std::is_same_v<Engine, RanluxppRngEngine>)
0075     {
0076         // If using Ranlux++, force the RNG to create a new block of bits
0077         // in preparation for taking a bunch of samples
0078         rng.advance();
0079     }
0080
0081     Real3 new_dir, new_pol;
0082     do
0083     {
0084         do
0085         {
0086             new_dir = sample_direction_(rng);
0087
0088             // Project polarization onto plane perpendicular to new direction
0089             new_pol = make_unit_vector(make_orthogonal(inc_pol_, new_dir));
0090
0091             // Reject rare case of polarization and new direction being
0092             // coincident leading to loss of orthogonality
0093         } while (CELER_UNLIKELY(!is_soft_orthogonal(new_pol, new_dir)));
0094
0095         if (!BernoulliDistribution{0.5}(rng))
0096         {
0097             // Flip direction with 50% probability: there are two polarizations
0098             // perpendicular to the new direction and the original polarization
0099             new_pol = -new_pol;
0100         }
0101         // Accept with the probability of the scattered polarization overlap
0102         // squared
0103     } while (reject_angle_(ipow<2>(dot_product(new_pol, inc_pol_)), rng));
0104
0105     CELER_ENSURE(is_soft_unit_vector(new_dir));
0106     CELER_ENSURE(is_soft_unit_vector(new_pol));
0107     Interaction result;
0108     result.direction = new_dir;
0109     result.polarization = new_pol;
0110     return result;
0111 }
0112
0113 //---------------------------------------------------------------------------//
0114 }  // namespace optical
0115 }  // namespace celeritas