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