optical/interactor/MieInteractor.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/MieInteractor.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Assert.hh"
0012 #include "corecel/Constants.hh"
0013 #include "corecel/Macros.hh"
0014 #include "corecel/Types.hh"
0015 #include "corecel/math/Algorithms.hh"
0016 #include "corecel/math/ArrayOperators.hh"
0017 #include "corecel/math/ArraySoftUnit.hh"
0018 #include "corecel/math/ArrayUtils.hh"
0019 #include "corecel/math/SoftEqual.hh"
0020 #include "corecel/random/distribution/BernoulliDistribution.hh"
0021 #include "corecel/random/distribution/RejectionSampler.hh"
0022 #include "corecel/random/distribution/UniformRealDistribution.hh"
0023 #include "celeritas/Types.hh"
0024 #include "celeritas/optical/Interaction.hh"
0025 #include "celeritas/optical/MieData.hh"
0026 #include "celeritas/optical/ParticleTrackView.hh"
0027
0028 namespace celeritas
0029 {
0030 namespace optical
0031 {
0032 //---------------------------------------------------------------------------//
0033 /*!
0034  * Sample optical Mie scattering using the Henyey-Greenstein distribution.
0035  *\f[
0036  *   P(\cos\theta) \propto \frac{1 - g^2}{(1 + g^2 - 2g\cos\theta)^{3/2}}
0037  * \f].
0038  * Parameters:
0039  * - forward_ratio: probability of using forward vs backward lobe
0040  * - forward_g, backward_g: HG asymmetry parameters for each lobe
0041  */
0042 class MieInteractor
0043 {
0044   public:
0045     // Construct with shared and state data
0046     inline CELER_FUNCTION MieInteractor(NativeCRef<MieData> const& shared,
0047                                         ParticleTrackView const& particle,
0048                                         Real3 const& direction,
0049                                         OptMatId const& mat_id);
0050
0051     // Sample an interaction with the given RNG
0052     template<class Engine>
0053     inline CELER_FUNCTION Interaction operator()(Engine& rng) const;
0054
0055   private:
0056     Real3 const& inc_dir_;  //!< Incident photon direction
0057     Real3 const& inc_pol_;  //!< Incident polarization
0058     MieMaterialData const& mie_params_;  //!< Mie scattering params
0059
0060     // Choose forward/backward scattering using random generator
0061     BernoulliDistribution sample_forward_;
0062 };
0063
0064 //---------------------------------------------------------------------------//
0065 // INLINE DEFINITIONS
0066 //---------------------------------------------------------------------------//
0067 /*!
0068  * Construct with shared and state data.
0069  */
0070 CELER_FUNCTION MieInteractor::MieInteractor(NativeCRef<MieData> const& shared,
0071                                             ParticleTrackView const& particle,
0072                                             Real3 const& direction,
0073                                             OptMatId const& mat_id)
0074     : inc_dir_(direction)
0075     , inc_pol_(particle.polarization())
0076     , mie_params_(shared.mie_record[mat_id])
0077     , sample_forward_(mie_params_.forward_ratio)
0078 {
0079     CELER_EXPECT(shared);
0080     CELER_EXPECT(mat_id < shared.mie_record.size());
0081     CELER_EXPECT(is_soft_unit_vector(inc_dir_));
0082     CELER_EXPECT(is_soft_unit_vector(inc_pol_));
0083     CELER_EXPECT(is_soft_orthogonal(inc_dir_, inc_pol_));
0084 }
0085
0086 //---------------------------------------------------------------------------//
0087 /*!
0088  * Sample a single optical Mie scattering event.
0089  */
0090 template<class Engine>
0091 CELER_FUNCTION Interaction MieInteractor::operator()(Engine& rng) const
0092 {
0093     using namespace celeritas::literals;
0094
0095     Interaction result;
0096     Real3& new_dir = result.direction;
0097     Real3& new_pol = result.polarization;
0098
0099     using UniformRealDist = UniformRealDistribution<real_type>;
0100     UniformRealDist sample_r(0, 1);
0101     do
0102     {
0103         real_type r = sample_r(rng);
0104         bool is_forward = sample_forward_(rng);
0105
0106         // Select forward or backward lobe
0107         real_type g
0108             = (is_forward ? mie_params_.forward_g : mie_params_.backward_g);
0109
0110         // Compute cos(theta) distribution for optical photons that undergo mie
0111         // scattering
0112         // NOTE: floating point cancellation means the expression can result in
0113         // a cosine slightly greater than unity
0114         real_type costheta = 2 * r * ipow<2>((1 + g) / (1 - g + 2 * g * r))
0115                                  * (1 - g + g * r)
0116                              - 1;
0117         costheta = celeritas::min(costheta, 1_r);
0118
0119         CELER_ASSERT(costheta >= -1 && costheta <= 1);
0120
0121         // Reverse cos theta if backward_g is chosen
0122         if (!is_forward)
0123             costheta = -costheta;
0124
0125         // Sample azimuthal angle
0126         UniformRealDist sample_phi(0, real_type(2 * constants::pi));
0127         real_type phi = sample_phi(rng);
0128
0129         // Creating new direction
0130         new_dir = from_spherical(costheta, phi);
0131
0132         // Project polarization onto plane perpendicular to new direction
0133         new_pol = make_unit_vector(make_orthogonal(inc_pol_, new_dir));
0134     } while (CELER_UNLIKELY(!is_soft_orthogonal(new_pol, new_dir)));
0135
0136     if (!BernoulliDistribution{0.5}(rng))
0137     {
0138         // Flip direction with 50% probability: there are two polarizations
0139         // perpendicular to the new direction and the original polarization
0140         new_pol = -new_pol;
0141     }
0142
0143     CELER_ENSURE(is_soft_unit_vector(new_dir));
0144     CELER_ENSURE(is_soft_unit_vector(new_pol));
0145     CELER_ENSURE(is_soft_orthogonal(new_pol, new_dir));
0146
0147     return result;
0148 }
0149
0150 //---------------------------------------------------------------------------//
0151 }  // namespace optical
0152 }  // namespace celeritas