surface/model/FresnelCalculator.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/surface/model/FresnelCalculator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Macros.hh"
0010 #include "celeritas/optical/MaterialView.hh"
0011 #include "celeritas/optical/ParticleTrackView.hh"
0012 #include "celeritas/optical/Types.hh"
0013 #include "celeritas/optical/surface/SurfacePhysicsUtils.hh"
0014
0015 #include "SurfaceInteraction.hh"
0016
0017 namespace celeritas
0018 {
0019 namespace optical
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Helper class for calculating quantities for Fresnel equations.
0024  *
0025  * Incident waves are decomposed into transverse-electric (TE) and
0026  * transverse-magnetic (TM) polarizations. The interaction plane is defined as
0027  * the span of the photon direction and the surface normal. The TE (TM)
0028  * polarization has the electric (magnetic) field polarization transverse to
0029  * the interaction plane. If the direction is parallel to the surface normal,
0030  * then interaction plane is degenerate and the incident photon is defined to
0031  * be entirely in the TE polarization.
0032  *
0033  * This calculator helps handle this degenerate case, and also handles
0034  * calculating reflectivity and transmission in the total internal reflection
0035  * case.
0036  *
0037  * Convention follows \cite{fowles-optics-1975} section 2.6, except that
0038  * photon phase is not tracked in Celeritas.
0039  */
0040 class FresnelCalculator
0041 {
0042   public:
0043     // Construct from initial state
0044     explicit inline CELER_FUNCTION FresnelCalculator(
0045         Real3 const& direction,
0046         Real3 const& polarization,
0047         Real3 const& normal,
0048         real_type relative_r_index);
0049
0050     // Construct from track views
0051     explicit inline CELER_FUNCTION FresnelCalculator(
0052         Real3 const& inc_direction,
0053         ParticleTrackView const& photon,
0054         Real3 const& normal,
0055         MaterialView const& pre_material,
0056         MaterialView const& post_material);
0057
0058     // Whether the interaction will be total internal reflection
0059     inline CELER_FUNCTION bool is_total_internal_reflection() const;
0060
0061     // Calculate total reflectivity
0062     inline CELER_FUNCTION real_type calc_reflectivity() const;
0063
0064     // Calculate interaction for refracted wave
0065     inline CELER_FUNCTION SurfaceInteraction refracted_interaction() const;
0066
0067   private:
0068     // Get refracted photon direction
0069     inline CELER_FUNCTION Real3 refracted_direction() const;
0070
0071     // Calculate transmission coefficients
0072     inline CELER_FUNCTION real_type calc_transmission_te() const;
0073     inline CELER_FUNCTION real_type calc_transmission_tm() const;
0074
0075     // Calculate reflectivity coefficients
0076     inline CELER_FUNCTION real_type calc_reflectivity_te() const;
0077     inline CELER_FUNCTION real_type calc_reflectivity_tm() const;
0078
0079     // Polarization axes
0080     inline CELER_FUNCTION Real3 const& te_axis() const;
0081     inline CELER_FUNCTION Real3 tm_axis(Real3 const& direction) const;
0082
0083     // Incident photon polarization components
0084     inline CELER_FUNCTION real_type inc_te_component() const;
0085     inline CELER_FUNCTION real_type inc_tm_component() const;
0086
0087     Real3 const& direction_;
0088     Real3 const& polarization_;
0089     Real3 const& normal_;
0090     real_type relative_r_index_;
0091
0092     real_type cos_theta_;
0093     real_type cosine_ratio_;
0094     Real3 p_axis_;
0095
0096     // Helper function for calculating reflectivity
0097     inline CELER_FUNCTION real_type reflectivity_ratio(real_type x) const;
0098 };
0099
0100 //---------------------------------------------------------------------------//
0101 // FREE FUNCTIONS
0102 //---------------------------------------------------------------------------//
0103
0104 inline CELER_FUNCTION real_type calc_relative_r_index(
0105     units::MevEnergy energy,
0106     MaterialView const& pre_material,
0107     MaterialView const& post_material)
0108 {
0109     auto calc_r = [energy](MaterialView const& mat) {
0110         return mat.make_refractive_index_calculator()(
0111             value_as<units::MevEnergy>(energy));
0112     };
0113
0114     return calc_r(post_material) / calc_r(pre_material);
0115 }
0116
0117 //---------------------------------------------------------------------------//
0118 // INLINE DEFINITIONS
0119 //---------------------------------------------------------------------------//
0120 /*!
0121  * Construct calculator from initial photon and surface physics data.
0122  */
0123 CELER_FUNCTION FresnelCalculator::FresnelCalculator(Real3 const& direction,
0124                                                     Real3 const& polarization,
0125                                                     Real3 const& normal,
0126                                                     real_type relative_r_index)
0127     : direction_(direction)
0128     , polarization_(polarization)
0129     , normal_(normal)
0130     , relative_r_index_(relative_r_index)
0131 {
0132     CELER_EXPECT(is_soft_unit_vector(direction_));
0133     CELER_EXPECT(is_soft_unit_vector(polarization_));
0134     CELER_EXPECT(is_soft_unit_vector(normal_));
0135     CELER_EXPECT(relative_r_index > 0);
0136     CELER_EXPECT(is_entering_surface(direction_, normal));
0137
0138     // Sometimes dot product of normalized parallel vectors is the next
0139     // representation after 1. Round down to exactly 1 to avoid errors.
0140     using namespace celeritas::literals;
0141     CELER_EXPECT(-dot_product(direction_, normal_) <= std::nextafter(1_r, 2_r));
0142     cos_theta_ = min(-dot_product(direction_, normal_), 1_r);
0143
0144     // Snell's law
0145     real_type sin_phi = sqrt(1 - ipow<2>(cos_theta_)) / relative_r_index_;
0146
0147     // If undergoing total internal reflection, set the cosine ratio to exactly
0148     // zero. This gives the correct reflectivity and transmission coefficients.
0149     cosine_ratio_ = sin_phi >= 1 ? 0 : sqrt(1 - ipow<2>(sin_phi)) / cos_theta_;
0150
0151     Real3 s_axis = make_orthogonal(direction_, normal);
0152     s_axis = soft_zero(norm(s_axis)) ? cross_product(polarization_, normal)
0153                                      : make_unit_vector(s_axis);
0154     p_axis_ = cross_product(normal_, s_axis);
0155 }
0156
0157 //---------------------------------------------------------------------------//
0158 /*!
0159  * Construct from track views.
0160  */
0161 CELER_FUNCTION FresnelCalculator::FresnelCalculator(
0162     Real3 const& inc_direction,
0163     ParticleTrackView const& photon,
0164     Real3 const& normal,
0165     MaterialView const& pre_material,
0166     MaterialView const& post_material)
0167     : FresnelCalculator(
0168           inc_direction,
0169           photon.polarization(),
0170           normal,
0171           calc_relative_r_index(photon.energy(), pre_material, post_material))
0172 {
0173 }
0174
0175 //---------------------------------------------------------------------------//
0176 /*!
0177  * Whether the photon is subject to total internal reflection.
0178  */
0179 CELER_FUNCTION bool FresnelCalculator::is_total_internal_reflection() const
0180 {
0181     // In the constructor, the cosine ratio is set to exactly zero for total
0182     // internal reflection.
0183     return cosine_ratio_ == 0;
0184 }
0185
0186 //---------------------------------------------------------------------------//
0187 /*!
0188  * Calculate total reflectivity for the incident photon.
0189  */
0190 CELER_FUNCTION real_type FresnelCalculator::calc_reflectivity() const
0191 {
0192     real_type te_comp_sq = ipow<2>(this->inc_te_component());
0193     real_type tm_comp_sq = ipow<2>(this->inc_tm_component());
0194     real_type total_reflectivity
0195         = (te_comp_sq * ipow<2>(this->calc_reflectivity_te())
0196            + tm_comp_sq * ipow<2>(this->calc_reflectivity_tm()))
0197           / (te_comp_sq + tm_comp_sq);
0198
0199     CELER_ENSURE(0 <= total_reflectivity && total_reflectivity <= 1);
0200
0201     return total_reflectivity;
0202 }
0203
0204 //---------------------------------------------------------------------------//
0205 /*!
0206  * Calculate interaction for refracted wave.
0207  */
0208 CELER_FUNCTION SurfaceInteraction FresnelCalculator::refracted_interaction() const
0209 {
0210     CELER_ASSERT(!this->is_total_internal_reflection());
0211
0212     SurfaceInteraction result;
0213     result.action = SurfaceInteraction::Action::refracted;
0214     result.direction = this->refracted_direction();
0215
0216     result.polarization = {0, 0, 0};
0217     axpy(this->calc_transmission_te() * this->inc_te_component(),
0218          this->te_axis(),
0219          &result.polarization);
0220     axpy(this->calc_transmission_tm() * this->inc_tm_component(),
0221          this->tm_axis(result.direction),
0222          &result.polarization);
0223     result.polarization = make_unit_vector(result.polarization);
0224
0225     CELER_ENSURE(result.is_valid());
0226
0227     return result;
0228 }
0229
0230 //---------------------------------------------------------------------------//
0231 /*!
0232  * Calculate direction of refracted photon.
0233  */
0234 CELER_FUNCTION Real3 FresnelCalculator::refracted_direction() const
0235 {
0236     CELER_EXPECT(!this->is_total_internal_reflection());
0237     Real3 dir = direction_;
0238     axpy(cos_theta_ * (1 - relative_r_index_ * cosine_ratio_), normal_, &dir);
0239     return make_unit_vector(dir);
0240 }
0241
0242 //---------------------------------------------------------------------------//
0243 /*!
0244  * Calculate transmission coefficient of the TE component.
0245  *
0246  * Derived from equation 2.52 of \cite{fowles-optics-1975}.
0247  */
0248 CELER_FUNCTION real_type FresnelCalculator::calc_transmission_te() const
0249 {
0250     return this->calc_reflectivity_te() + 1;
0251 }
0252
0253 //---------------------------------------------------------------------------//
0254 /*!
0255  * Calculate transmission coefficient of the TM component.
0256  *
0257  * Derived from equation 2.53 of \cite{fowles-optics-1975}.
0258  */
0259 CELER_FUNCTION real_type FresnelCalculator::calc_transmission_tm() const
0260 {
0261     return relative_r_index_ * (this->calc_reflectivity_tm() + 1);
0262 }
0263
0264 //---------------------------------------------------------------------------//
0265 /*!
0266  * Calculate reflectivity coefficient of the TE component.
0267  *
0268  * Equivalent to equation 2.54 of \cite{fowles-optics-1975}.
0269  */
0270 CELER_FUNCTION real_type FresnelCalculator::calc_reflectivity_te() const
0271 {
0272     return -this->reflectivity_ratio(cosine_ratio_ * relative_r_index_);
0273 }
0274
0275 //---------------------------------------------------------------------------//
0276 /*!
0277  * Calculate reflectivity coefficient of the TM component.
0278  *
0279  * Equivalent to equation 2.55 of \cite{fowles-optics-1975}.
0280  */
0281 CELER_FUNCTION real_type FresnelCalculator::calc_reflectivity_tm() const
0282 {
0283     return this->reflectivity_ratio(cosine_ratio_ / relative_r_index_);
0284 }
0285
0286 //---------------------------------------------------------------------------//
0287 /*!
0288  * Get the polarization vector for the TE component.
0289  */
0290 CELER_FUNCTION Real3 const& FresnelCalculator::te_axis() const
0291 {
0292     return p_axis_;
0293 }
0294
0295 //---------------------------------------------------------------------------//
0296 /*!
0297  * Get the polarization vector for the TM component based on the given
0298  * direction.
0299  */
0300 CELER_FUNCTION Real3 FresnelCalculator::tm_axis(Real3 const& direction) const
0301 {
0302     return cross_product(this->te_axis(), direction);
0303 }
0304
0305 //---------------------------------------------------------------------------//
0306 /*!
0307  * Get the incident photon TE polarization component.
0308  */
0309 CELER_FUNCTION real_type FresnelCalculator::inc_te_component() const
0310 {
0311     return dot_product(polarization_, this->te_axis());
0312 }
0313
0314 //---------------------------------------------------------------------------//
0315 /*!
0316  * Get the incident photon TM polarization component.
0317  */
0318 CELER_FUNCTION real_type FresnelCalculator::inc_tm_component() const
0319 {
0320     return dot_product(polarization_, this->tm_axis(direction_));
0321 }
0322
0323 //---------------------------------------------------------------------------//
0324 /*!
0325  * Helper function for calculating reflectivity coefficients.
0326  */
0327 CELER_FUNCTION real_type FresnelCalculator::reflectivity_ratio(
0328     real_type x) const
0329 {
0330     return (x - 1) / (x + 1);
0331 }
0332
0333 //---------------------------------------------------------------------------//
0334 }  // namespace optical
0335 }  // namespace celeritas