Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-09 07:49:31 

0001 #pragma once
0002 /**
0003 sboundary.h
0004 =============
0005 
0006 Riffing on qsim.h:propagate_at_boundary to assist with comparison with stmm.h 
0007 
0008 comparing stmm.h with sboundary.h 
0009 -------------------------------------
0010 
0011 Q: pol comes in much later with stmm.h, how does it manage that ? 
0012 A: By calculating the s and p coeffs and then only applying them to the actual pol as the final step 
0013 
0014    * in the sboundary.h expressions below that is factoring off the E1.x E1.y 
0015    * PERHAPS CAN FACTORIZE sboundary.h ANALOGOUSLY ?
0016 
0017 
0018 From u4/CustomBoundary.h:doIt::
0019 
0020     267     const G4ThreeVector& surface_normal = theRecoveredNormal ;
0021     268     const G4ThreeVector& direction      = OldMomentum ;
0022     269     const G4ThreeVector& polarization   = OldPolarization ;
0023     ...
0024     271     G4double minus_cos_theta = direction*surface_normal ;
0025     272     G4double orientation = minus_cos_theta < 0. ? 1. : -1.  ;
0026     273     G4ThreeVector oriented_normal = orientation*surface_normal ;
0027     ...
0028     315     const double _si = stack.ll[0].st.real() ;
0029     ...
0030     327     // E_s2 : S-vs-P power fraction : signs make no difference as squared
0031     328     double E_s2 = _si > 0. ? (polarization*direction.cross(oriented_normal))/_si : 0. ;
0032   
0033     ///
0034     ///                                                                              ^^^  trans_length
0035     ///                                            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^        trans
0036     ///                                            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A_trans
0037     ///                               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ E1_perp
0038     ///
0039 
0040     329     E_s2 *= E_s2;
0041 
0042     ////    ^^^^^   E1_perp*E1_perp
0043 
0044     ...
0045     340     double fT_s = stack.art.T_s ;
0046     341     double fT_p = stack.art.T_p ;
0047     342     double fR_s = stack.art.R_s ;
0048     343     double fR_p = stack.art.R_p ;
0049     344     double one = 1.0 ;
0050     345     double T = fT_s*E_s2 + fT_p*(one-E_s2);  // THIS MATCHES TransCoeff from sboundary.h 
0051     346     double R = fR_s*E_s2 + fR_p*(one-E_s2);
0052     347     double A = one - (T+R);
0053 
0054 **/
0055 
0056 
0057 struct sboundary
0058 {
0059     sphoton& p ; 
0060 
0061     const sstate& s ; 
0062     const float& n1 ;
0063     const float& n2 ;    
0064     const float eta ;  
0065     const float3* normal ; // geometrical outwards normal 
0066     const float mct ; 
0067     const float orient ; 
0068     const float c1 ; 
0069     const float c2c2 ; 
0070     const bool tir ;
0071     const float c2 ; 
0072     const float n1c1 ; 
0073     const float n2c2 ;
0074     const float n2c1 ;
0075     const float n1c2 ;
0076     const float2 _E2_t ;
0077 
0078     const float3 transverse ; 
0079     const float transverse_length ; 
0080     const bool  normal_incidence ; 
0081     const float3 A_transverse ; 
0082     const float E1_perp ; 
0083 
0084     const float2 E1   ;
0085     const float2 E1_chk   ;
0086     const float2 E2_t ;
0087     const float2 _E2_r ;
0088     const float2 E2_r ;
0089     const float2 RR ;
0090     const float2 TT ; 
0091 
0092     const float TransCoeff ; 
0093     const float ReflectCoeff ; 
0094     const float u_reflect ;
0095     const bool reflect ; 
0096     const unsigned flag ; 
0097     const float Coeff ; 
0098     const float EdotN ; 
0099     float3 A_parallel ; 
0100     float3 alt_pol ;  // check an alternative polarization expression  
0101 
0102     sboundary(RNG& rng, sctx& ctx );  
0103 };
0104 
0105 inline sboundary::sboundary( RNG& rng, sctx& ctx ) 
0106     :
0107     p(ctx.p),
0108     s(ctx.s),
0109     n1(s.material1.x),
0110     n2(s.material2.x),
0111     eta(n1/n2),
0112     normal((float3*)&ctx.prd->q0.f.x),
0113     mct(dot(p.mom, *normal )),
0114     orient( mct < 0.f ? 1.f : -1.f ),
0115     c1(fabs(mct)),
0116     c2c2(1.f - eta*eta*(1.f - c1 * c1 )),   // Snells law and trig identity 
0117     tir(c2c2 < 0.f),
0118     c2( tir ? 0.f : sqrtf(c2c2)),
0119     n1c1(n1*c1),
0120     n2c2(n2*c2),
0121     n2c1(n2*c1),
0122     n1c2(n1*c2),
0123     _E2_t(make_float2( 2.f*n1c1/(n1c1+n2c2)   , 2.f*n1c1/(n2c1+n1c2) )),  // (ts,tp)  matches: real(c.stack.comp.ts,tp) 
0124     _E2_r(make_float2(  _E2_t.x - 1.f         , n2*_E2_t.y/n1 - 1.f  )),  // (rs,rp)  matches: real(c.stack.comp.rs,rp) 
0125     transverse(cross(p.mom, orient*(*normal))),
0126     transverse_length(length(transverse)),
0127     normal_incidence(transverse_length < 1e-6f ),
0128     A_transverse(normal_incidence ? p.pol : transverse/transverse_length),
0129     E1_perp(dot(p.pol, A_transverse)),
0130     E1(normal_incidence ? make_float2( 0.f, 1.f) : make_float2( E1_perp , length( p.pol - (E1_perp*A_transverse) ) )),   // ( S, P )
0131     E1_chk(make_float2( E1_perp, sqrtf(1.f - E1_perp*E1_perp) )),  // HMM: .x can be -ve, .y chosen +ve 
0132     E2_t(_E2_t*E1),        // deferred using pol to make closer to stmm.h 
0133     E2_r(_E2_r*E1),
0134     RR(normalize(E2_r)),   // elementwise multiplication is like non-uniform scaling, so cannot factor out E1 here 
0135     TT(normalize(E2_t)), 
0136     TransCoeff(tir || n1c1 == 0.f ? 0.f : n2c2*dot(E2_t,E2_t)/n1c1),
0137     ReflectCoeff(1.f - TransCoeff),
0138     u_reflect(curand_uniform(&rng)),
0139     reflect(u_reflect > TransCoeff), 
0140     flag(reflect ? BOUNDARY_REFLECT : BOUNDARY_TRANSMIT),
0141     Coeff(reflect ? ReflectCoeff : TransCoeff),
0142     EdotN(orient*dot(p.pol, *normal))
0143 {
0144     p.mom = reflect
0145                     ?
0146                        p.mom + 2.0f*c1*orient*(*normal)
0147                     :
0148                        eta*(p.mom) + (eta*c1 - c2)*orient*(*normal)
0149                     ;
0150 
0151     A_parallel = normalize(cross(p.mom, A_transverse));   // A_paral is P-pol direction using new p.mom
0152 
0153 
0154     alt_pol = reflect ?
0155                           p.pol - 2.f*orient*dot(p.pol,*normal)*(*normal)   // opposite sign to below expression
0156                       :
0157                           normalize((p.pol-(dot(p.pol,p.mom)*p.mom)));      // 
0158                       ;
0159 
0160 
0161     p.pol =  normal_incidence 
0162                     ?
0163                        ( reflect ?  p.pol*(n2>n1? -1.f:1.f) : p.pol )
0164                     :
0165                        ( reflect ?
0166                                     ( tir ?  -p.pol + 2.f*EdotN*orient*(*normal) : RR.x*A_transverse + RR.y*A_parallel )
0167                                  :
0168                                     TT.x*A_transverse + TT.y*A_parallel
0169                        )
0170                     ;
0171 
0172     
0173 
0174 
0175 
0176 }
0177 
0178 inline std::ostream& operator<<(std::ostream& os, const sboundary& b)  
0179 {   
0180     os 
0181        << std::setw(20) << " n1 "           << std::setw(10) << std::fixed << std::setprecision(4) << b.n1  << std::endl 
0182        << std::setw(20) << " n2 "           << std::setw(10) << std::fixed << std::setprecision(4) << b.n2  << std::endl 
0183        << std::setw(20) << " eta "          << std::setw(10) << std::fixed << std::setprecision(4) << b.eta  << std::endl 
0184        << std::setw(20) << " normal "       << std::setw(10) << std::fixed << std::setprecision(4) << *b.normal  << std::endl
0185        << std::setw(20) << " mct "           << std::setw(10) << std::fixed << std::setprecision(4) << b.mct  << std::endl 
0186        << std::setw(20) << " orient "           << std::setw(10) << std::fixed << std::setprecision(4) << b.orient  << std::endl 
0187        << std::setw(20) << " c1 "           << std::setw(10) << std::fixed << std::setprecision(4) << b.c1  << std::endl 
0188        << std::setw(20) << " c2c2 "           << std::setw(10) << std::fixed << std::setprecision(4) << b.c2c2  << std::endl 
0189        << std::setw(20) << " tir "           << std::setw(10) << b.tir  << std::endl 
0190        << std::setw(20) << " c2 "           << std::setw(10) << std::fixed << std::setprecision(4) << b.c2  
0191        << std::endl 
0192        << std::setw(20) << " _E2_t "        <<  std::setw(10) << b._E2_t  
0193        << std::setw(20) << " length(_E2_t) " << std::setw(10) << std::fixed << std::setprecision(4) << length(b._E2_t)
0194        << std::endl
0195        << std::setw(20) << " _E2_r "        <<  std::setw(10) << b._E2_r 
0196        << std::setw(20) << " length(_E2_r) " << std::setw(10) << std::fixed << std::setprecision(4) << length(b._E2_r)
0197        << std::endl
0198        << std::setw(20) << " transverse "        << std::setw(10) << b.transverse
0199        << std::setw(20) << " length(transverse) " << std::setw(10) << std::fixed << std::setprecision(4) << length(b.transverse)  
0200        << std::endl
0201        << std::setw(20) << " transverse_length " << std::setw(10) << std::fixed << std::setprecision(4) <<  b.transverse_length  << std::endl
0202        << std::setw(20) << " normal_incidence " << std::setw(10) << std::fixed << std::setprecision(4) <<  b.normal_incidence  
0203        << std::endl
0204        << std::setw(20) << " A_transverse "      << std::setw(10) <<  b.A_transverse  
0205        << std::setw(20) << " length(A_transverse) " << std::setw(10) << std::fixed << std::setprecision(4) <<  length(b.A_transverse)  
0206        << std::endl
0207        << std::setw(20) << " E1_perp "      << std::setw(10) << std::fixed << std::setprecision(4) <<  b.E1_perp  
0208        << std::endl
0209        << std::setw(20) << " E1 "           <<  std::setw(10) << b.E1  
0210        << std::setw(20) << " length(E1) "   <<  std::setw(10) << std::fixed << std::setprecision(4) << length(b.E1)
0211        << std::endl
0212        << std::setw(20) << " E1_chk "       << std::setw(10) << b.E1_chk  
0213        << std::setw(20) << " length(E1) "   <<  std::setw(10) << std::fixed << std::setprecision(4) << length(b.E1_chk)
0214        << std::endl
0215        << std::setw(20) << " E2_t "         << std::setw(10)<<  b.E2_t  
0216        << std::setw(20) << " length(E2_t) " << std::setw(10) << std::fixed << std::setprecision(4) <<  length(b.E2_t)  
0217        << std::endl
0218        << std::setw(20) << " E2_r "         << std::setw(10) << b.E2_r  
0219        << std::setw(20) << " length(E2_r) " << std::setw(10) << std::fixed << std::setprecision(4) <<  length(b.E2_r)  
0220        << std::endl
0221        << std::setw(20) << " RR "           << std::setw(10) <<  b.RR 
0222        << std::setw(20) << " length(RR) "   << std::setw(10) << std::fixed << std::setprecision(4) <<  length(b.RR)  
0223        << std::endl
0224        << std::setw(20) << " TT "           << std::setw(10) <<  b.TT  
0225        << std::setw(20) << " length(TT) "   << std::setw(10) << std::fixed << std::setprecision(4) <<  length(b.TT)  
0226        << std::endl
0227        << std::setw(20) << " TransCoeff "   << std::setw(10) << std::fixed << std::setprecision(4) << b.TransCoeff  << std::endl
0228        << std::setw(20) << " ReflectCoeff " << std::setw(10) << std::fixed << std::setprecision(4) << b.ReflectCoeff  << std::endl
0229        << std::setw(20) << " u_reflect "    << std::setw(10) << std::fixed << std::setprecision(4) << b.u_reflect  << std::endl
0230        << std::setw(20) << " reflect "      << std::setw(10) << std::fixed << std::setprecision(4) << b.reflect  << std::endl
0231        << std::setw(20) << " flag "         << std::setw(10) << std::fixed << std::setprecision(4) << OpticksPhoton::Flag(b.flag)  << std::endl 
0232        << std::setw(20) << " Coeff "        << std::setw(10) << std::fixed << std::setprecision(4) << b.Coeff  << std::endl
0233        << std::setw(20) << " EdotN "        << std::setw(10) << std::fixed << std::setprecision(4) << b.EdotN  
0234        << std::endl
0235        << std::setw(20) << " p.mom "        << std::setw(10) << b.p.mom  
0236        << std::setw(20) << " length(p.mom) " << std::setw(10) << std::fixed << std::setprecision(4) << length(b.p.mom)  
0237        << std::endl 
0238        << std::setw(20) << " A_parallel "   << std::setw(10) << b.A_parallel 
0239        << std::setw(20) << " length(A_parallel) " << std::setw(10) << std::fixed << std::setprecision(4) << length(b.A_parallel) 
0240        << std::endl
0241        << std::setw(20) << " p.pol "        << std::setw(10) << b.p.pol  
0242        << std::setw(20) << " length(p.pol) " << std::setw(10) << std::fixed << std::setprecision(4) << length(b.p.pol)  
0243        << std::endl 
0244        << std::setw(20) << " dot(p.mom,p.pol) " << std::setw(10) << std::fixed << std::setprecision(4) << dot(b.p.mom,b.p.pol)  
0245        << std::endl 
0246        << std::setw(20) << " alt_pol "        << std::setw(10) << b.alt_pol  
0247        << std::setw(20) << " length(alt_pol) " << std::setw(10) << std::fixed << std::setprecision(4) << length(b.alt_pol)  
0248        << std::endl 
0249        << std::setw(20) << " dot(p.mom,alt_pol) " << std::setw(10) << std::fixed << std::setprecision(4) << dot(b.p.mom,b.alt_pol)  
0250        << std::endl 
0251        ;
0252     return os; 
0253 }   
0254
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!