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 #include "celeritas/phys/InteractionUtils.hh"
0020
0021 namespace celeritas
0022 {
0023 namespace optical
0024 {
0025 //---------------------------------------------------------------------------//
0026 /*!
0027  * Sample optical Rayleigh scattering.
0028  *
0029  * Optical Rayleigh scattering is the elastic scattering of optical photons
0030  * in a material. The scattered polarization is guarenteed to be in the
0031  * same plane as the original polarization and new direction if the latter
0032  * vectors are not parallel. Otherwise it will just be perpendicular to
0033  * the new direction.
0034  */
0035 class RayleighInteractor
0036 {
0037   public:
0038     // Construct interactor from an optical track
0039     inline CELER_FUNCTION RayleighInteractor(ParticleTrackView const& particle,
0040                                              Real3 const& direction);
0041
0042     // Sample an interaction with the given RNG
0043     template<class Engine>
0044     inline CELER_FUNCTION Interaction operator()(Engine& rng) const;
0045
0046   private:
0047     Real3 const& inc_dir_;  //!< Direction of incident photon
0048     Real3 const& inc_pol_;  //!< Polarization of incident photon
0049 };
0050
0051 //---------------------------------------------------------------------------//
0052 // INLINE DEFINITIONS
0053 //---------------------------------------------------------------------------//
0054 /*!
0055  * Construct the interactor for the given optical track.
0056  */
0057 CELER_FUNCTION
0058 RayleighInteractor::RayleighInteractor(ParticleTrackView const& particle,
0059                                        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(soft_zero(dot_product(inc_dir_, inc_pol_)));
0065 }
0066
0067 //---------------------------------------------------------------------------//
0068 /*!
0069  * Sample the scattered direction and polarization for a single optical
0070  * Rayliegh interaction.
0071  */
0072 template<class Engine>
0073 CELER_FUNCTION Interaction RayleighInteractor::operator()(Engine& rng) const
0074 {
0075     Interaction result;
0076     Real3& new_dir = result.direction;
0077     Real3& new_pol = result.polarization;
0078
0079     IsotropicDistribution sample_direction{};
0080     do
0081     {
0082         new_dir = sample_direction(rng);
0083
0084         auto projected_pol = dot_product(new_dir, inc_pol_);
0085
0086         if (soft_zero(projected_pol))
0087         {
0088             // If new direction is parallel to incident polarization, then
0089             // randomly sample azimuthal direction
0090             new_pol = ExitingDirectionSampler{0, new_dir}(rng);
0091         }
0092         else
0093         {
0094             // Project polarization onto plane perpendicular to new direction
0095             new_pol = make_unit_vector(inc_pol_ - projected_pol * new_dir);
0096             new_pol *= BernoulliDistribution{real_type{0.5}}(rng) ? 1 : -1;
0097         }
0098
0099         // Accept with the probability of the scattered polarization overlap
0100         // squared
0101     } while (RejectionSampler{ipow<2>(dot_product(new_pol, inc_pol_))}(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