optical/interactor/RayleighInteractor.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/optical/interactor/RayleighInteractor.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/math/ArrayOperators.hh"
0012 #include "corecel/math/ArraySoftUnit.hh"
0013 #include "corecel/math/ArrayUtils.hh"
0014 #include "corecel/math/SoftEqual.hh"
0015 #include "celeritas/optical/Interaction.hh"
0016 #include "celeritas/optical/ParticleTrackView.hh"
0017 #include "celeritas/phys/InteractionUtils.hh"
0018 #include "celeritas/random/distribution/BernoulliDistribution.hh"
0019 #include "celeritas/random/distribution/IsotropicDistribution.hh"
0020 #include "celeritas/random/distribution/RejectionSampler.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 guarenteed to be in the
0032  * same plane as the original polarization and new direction if the latter
0033  * vectors are not parallel. Otherwise it will just be perpendicular to
0034  * the new direction.
0035  */
0036 class RayleighInteractor
0037 {
0038   public:
0039     // Construct interactor from an optical track
0040     inline CELER_FUNCTION RayleighInteractor(ParticleTrackView const& particle,
0041                                              Real3 const& direction);
0042
0043     // Sample an interaction with the given RNG
0044     template<class Engine>
0045     inline CELER_FUNCTION Interaction operator()(Engine& rng) const;
0046
0047   private:
0048     Real3 const& inc_dir_;  //!< Direction of incident photon
0049     Real3 const& inc_pol_;  //!< Polarization of incident photon
0050 };
0051
0052 //---------------------------------------------------------------------------//
0053 // INLINE DEFINITIONS
0054 //---------------------------------------------------------------------------//
0055 /*!
0056  * Construct the interactor for the given optical track.
0057  */
0058 CELER_FUNCTION
0059 RayleighInteractor::RayleighInteractor(ParticleTrackView const& particle,
0060                                        Real3 const& direction)
0061     : inc_dir_(direction), inc_pol_(particle.polarization())
0062 {
0063     CELER_EXPECT(is_soft_unit_vector(inc_dir_));
0064     CELER_EXPECT(is_soft_unit_vector(inc_pol_));
0065     CELER_EXPECT(soft_zero(dot_product(inc_dir_, inc_pol_)));
0066 }
0067
0068 //---------------------------------------------------------------------------//
0069 /*!
0070  * Sample the scattered direction and polarization for a single optical
0071  * Rayliegh interaction.
0072  */
0073 template<class Engine>
0074 CELER_FUNCTION Interaction RayleighInteractor::operator()(Engine& rng) const
0075 {
0076     Interaction result;
0077     Real3& new_dir = result.direction;
0078     Real3& new_pol = result.polarization;
0079
0080     IsotropicDistribution sample_direction{};
0081     do
0082     {
0083         new_dir = sample_direction(rng);
0084
0085         auto projected_pol = dot_product(new_dir, inc_pol_);
0086
0087         if (soft_zero(projected_pol))
0088         {
0089             // If new direction is parallel to incident polarization, then
0090             // randomly sample azimuthal direction
0091             new_pol = ExitingDirectionSampler{0, new_dir}(rng);
0092         }
0093         else
0094         {
0095             // Project polarization onto plane perpendicular to new direction
0096             new_pol = make_unit_vector(inc_pol_ - projected_pol * new_dir);
0097             new_pol *= BernoulliDistribution{real_type{0.5}}(rng) ? 1 : -1;
0098         }
0099
0100         // Accept with the probability of the scattered polarization overlap
0101         // squared
0102     } while (RejectionSampler{ipow<2>(dot_product(new_pol, inc_pol_))}(rng));
0103
0104     CELER_ENSURE(is_soft_unit_vector(new_dir));
0105     CELER_ENSURE(is_soft_unit_vector(new_pol));
0106     CELER_ENSURE(soft_zero(dot_product(new_pol, new_dir)));
0107
0108     return result;
0109 }
0110
0111 //---------------------------------------------------------------------------//
0112 }  // namespace optical
0113 }  // namespace celeritas