sysrap/tests/stmm_vs_sboundary_test.cc

0001 /**
0002 stmm_vs_sboundary_test.cc
0003 ============================
0004
0005 ::
0006
0007    ~/o/sysrap/tests/stmm_vs_sboundary_test.sh
0008
0009
0010 Specialize the stmm.h stack calc to 2 non-thin layers and compare with sboundary.h
0011
0012 Had to adapt the sboundary.h polarization calc as what is in junoPMTOpticalModel
0013 looks clearly wrong.
0014 sboundary.h is a riff on qsim:propagate_at_boundary which is based on G4OpBoundaryProcess.
0015
0016 * But might be missing some subtlety ?
0017
0018 BUT: can the leap be made for the mom and pol appropriate to the real usage
0019 of a stack with 4 layers (2 thin in the middle) ?
0020
0021
0022 :google:`refraction stack of layers`
0023
0024 * https://www.reading.ac.uk/infrared/technical-library/thin-film-multilayers
0025
0026 * https://physics.stackexchange.com/questions/420069/resulting-refractive-index-of-multiple-layers-with-different-indexes
0027
0028 Snell's Law::
0029
0030     n1 s1 = n2 s2 = n3 s3 = n4 s4 = n5 s5
0031
0032 To find the angle at which the ray emerges you only need to consider the 2
0033 mediums in which the ray enters and emerges. It does not matter what layers it
0034 has passed through in between. If these 2 mediums are the same (eg air) then
0035 the beam emerges parallel to the direction from which it entered.
0036
0037 Neither can the theorem tell you whether the ray undergoes total internal
0038 reflection at one of the interfaces between layers of material. In that case
0039 the ray will emerge from the face which it entered, or a side face if the slab
0040 of multi-layer material is not wide enough. To find out if this happens you
0041 need to trace the ray through the layers, checking what happens at each
0042 boundary.
0043
0044
0045
0046 * https://repository.tudelft.nl/islandora/object/uuid:3a6b5efa-99d9-481a-9185-6ef14d13429c/datastream/OBJ/download
0047 * ~/opticks_refs/multilayer_Thesis_Version_XII.pdf
0048
0049
0050
0051 **/
0052
0053 #include "srngcpu.h"
0054 using RNG = srngcpu ;
0055
0056 #include "sphoton.h"
0057 #include "sstate.h"
0058 #include "srec.h"
0059 #include "sseq.h"
0060 #include "stag.h"
0061 #include "sevent.h"
0062 #include "sctx.h"
0063 //#include "scurand.h"
0064 #include "sboundary.h"
0065 #include "stmm.h"
0066
0067 const char* FOLD = getenv("FOLD") ;
0068
0069 struct scf
0070 {
0071     float frac_twopi ;
0072     float spare01 ;
0073     float spare02 ;
0074     float spare03 ;
0075
0076     float c_T ;
0077     float b_TransCoeff ;
0078     float spare12 ;
0079     float spare13 ;
0080
0081     float c_transverse_length ;
0082     float b_transverse_length ;
0083     float spare22 ;
0084     float spare23 ;
0085
0086     float c_E1_perp ;
0087     float b_E1_perp ;
0088     float spare32 ;
0089     float spare33 ;
0090 };
0091
0092 inline std::ostream& operator<<(std::ostream& os, const scf& cf)
0093 {
0094     os
0095         << std::setw(25) << " frac_twopi "          << std::setw(10) << std::fixed << std::setprecision(10) << cf.frac_twopi
0096         << std::endl
0097         << std::setw(25) << " c_T "                  << std::setw(10) << std::fixed << std::setprecision(10) << cf.c_T
0098         << std::setw(25) << " b_TransCoeff  "       << std::setw(10) << std::fixed << std::setprecision(10)  << cf.b_TransCoeff
0099         << std::setw(25) << " c-b "                  << std::setw(10) << std::fixed << std::setprecision(10)  << ( cf.c_T-cf.b_TransCoeff )
0100         << std::endl
0101         << std::setw(25) << " c_transverse_length " << std::setw(10) << std::fixed << std::setprecision(10) << cf.c_transverse_length
0102         << std::setw(25) << " b_transverse_length " << std::setw(10) << std::fixed << std::setprecision(10) << cf.b_transverse_length
0103         << std::setw(25) << " c-b "                 << std::setw(10) << std::fixed << std::setprecision(10) << ( cf.c_transverse_length - cf.b_transverse_length )
0104         << std::endl
0105         << std::setw(25) << " c_E1_perp "  << std::setw(10) << std::fixed << std::setprecision(10) << cf.c_E1_perp
0106         << std::setw(25) << " b_E1_perp "  << std::setw(10) << std::fixed << std::setprecision(10) << cf.b_E1_perp
0107         << std::setw(25) << " c-b "         << std::setw(10) << std::fixed << std::setprecision(10) << ( cf.c_E1_perp - cf.b_E1_perp )
0108         << std::endl
0109         ;
0110     return os ;
0111 }
0112
0113
0114 struct Custom
0115 {
0116     static constexpr const float one = 1.f ;
0117     const float3& nrm ;
0118     const float3& mom ;
0119     const float wl ;
0120     const float mct ;
0121     const float orient ;
0122     const float3 _transverse ;
0123     const float  transverse_length ;  // same as s1 (sine of incident angle)
0124     const bool   normal_incidence ;
0125     const float3 transverse ;
0126
0127     StackSpec<float,2> ss ;
0128     Stack<float,2> stack ;
0129
0130     const float c1 ;
0131     const float c2 ;
0132     const float n1 ;
0133     const float n2 ;
0134     const float eta ;
0135     const float c2c2 ;
0136     const bool tir ;
0137
0138     const float2 fT ;
0139     const float2 fR ;
0140
0141     const float2 _E2_t ;
0142     const float2 _E2_r ;
0143
0144     Custom(const float3& nrm, const float3& mom, float n1, float n2, float wl );
0145     void set_polarization(const float3& polarization_) ;
0146
0147     float3 polarization ;
0148
0149     float EdotN ;
0150     float E1_perp ;
0151     float E_s2 ;
0152
0153     float2 E1 ;
0154     float2 E1E1 ;
0155
0156     float T ;
0157     float R ;
0158     float A ;
0159
0160     float2 E2_t ;
0161     float2 E2_r ;
0162
0163     float2 TT ;
0164     float2 RR ;
0165
0166     void set_reflect( bool reflect_ );
0167
0168     bool reflect ;
0169     float3 new_mom ;
0170     float3 parallel ;
0171     float3 new_pol ;
0172 };
0173
0174 inline Custom::Custom(const float3& nrm_, const float3& mom_, float n1_, float n2_, float wl_)
0175     :
0176     nrm(nrm_),
0177     mom(mom_),
0178     wl(wl_),
0179     mct(dot(mom, nrm)),
0180     orient( mct < 0.f ? 1.f : -1.f),
0181     _transverse(orient*cross(mom,nrm)),
0182     transverse_length(length(_transverse)),
0183     normal_incidence(transverse_length < 1e-6f ),
0184     transverse(normal_incidence ? make_float3(0.f, 0.f, 0.f) : _transverse/transverse_length),
0185     ss(StackSpec<float,2>::Create2(n1_,n2_)),
0186     stack(wl, mct, ss ),
0187     c1(real(stack.ll[0].ct)),
0188     c2(real(stack.ll[1].ct)),
0189     n1(real(stack.ll[0].n)),
0190     n2(real(stack.ll[1].n)),
0191     c2c2(1.f - eta*eta*(1.f - c1 * c1 )),   // Snells law and trig identity
0192     tir(c2c2 < 0.f),   // HMM: maybe can just check imag(stack.ll[1].ct) ?
0193     eta(n1/n2),
0194     fT(make_float2(stack.art.T_s, stack.art.T_p)),
0195     fR(make_float2(stack.art.R_s, stack.art.R_p)),
0196     _E2_t(make_float2(real(stack.comp.ts), real(stack.comp.tp))),
0197     _E2_r(make_float2(real(stack.comp.rs), real(stack.comp.rp)))
0198 {
0199     assert( n1 == n1_ );
0200     assert( n2 == n2_ );
0201 }
0202
0203 inline void Custom::set_polarization(const float3& polarization_)
0204 {
0205     polarization = polarization_ ;
0206
0207     EdotN = orient*dot(polarization, nrm) ;
0208     E1_perp = dot(polarization,transverse);
0209     E_s2 = E1_perp*E1_perp ;
0210     E1   = make_float2( E1_perp, sqrtf(1.f - E_s2)) ;
0211     E1E1 = make_float2( E_s2, 1.f - E_s2);
0212     //E1E1 = E1*E1 ;
0213
0214     T = dot( fT, E1E1 );    // equiv to fT.x*E_s2 + fT.y*(one-E_s2)
0215     R = dot( fR, E1E1 );    // equiv to fR.x*E_s2 + fR.y*(one-E_s2)
0216     A = one - (T + R );
0217
0218     E2_t = E1*_E2_t ;
0219     E2_r = E1*_E2_r ;
0220
0221     TT  = normalize(E2_t) ;
0222     RR  = normalize(E2_r) ;
0223 }
0224
0225 inline void Custom::set_reflect( bool reflect_ )
0226 {
0227     reflect = reflect_ ;
0228
0229     new_mom = reflect
0230                     ?
0231                        mom + 2.0f*c1*orient*nrm
0232                     :
0233                        eta*mom + (eta*c1 - c2)*orient*nrm
0234                     ;
0235
0236     parallel = normalize(cross(new_mom, transverse));   // A_paral is P-pol direction using new p.mom
0237
0238     new_pol =  normal_incidence
0239                     ?
0240                        ( reflect ?  polarization*(n2>n1? -1.f:1.f) : polarization )
0241                     :
0242                        ( reflect ?
0243                                     ( tir ?  -polarization + 2.f*EdotN*orient*nrm : RR.x*transverse + RR.y*parallel )
0244                                  :
0245                                     TT.x*transverse + TT.y*parallel
0246                        )
0247                     ;
0248 }
0249
0250 /**
0251 T matches TransCoeff from sboundary.h despite very different impl, hopefully acts as Rosetta stone
0252 -----------------------------------------------------------------------------------------------------
0253
0254 The below demonstrates the match in the expressions, as done numerically already:
0255
0256 stmm.h::
0257
0258     fT.x = stack.art.T_s = (n2c2/n1c1)*ts*ts
0259     fT.y = stack.art.T_p = (n2c2/n1c1)*tp*tp
0260
0261     E_s2 = E1_perp*E1_perp
0262
0263     T = fT.x*E_s2 + fT.y*(one-E_s2);
0264       = (n2c2/n1c1)[ ts*ts*E_s2 + tp*tp*(1-E_s2) ]
0265
0266 sboundary.h::
0267
0268     E1 = E1_chk(make_float2( E1_perp, sqrtf(1.f - E1_perp*E1_perp) ))
0269     E1 =    ( E1_perp, sqrt( 1-E1_perp*E1_perp ))
0270     E1 =    ( E1_perp, sqrt( 1-E1_perp*E1_perp ))
0271
0272     _E2_t =  (ts,tp)
0273     E2_t  = (_E2_t*E1)
0274
0275     TransCoeff = n2c2*dot(E2_t,E2_t)/n1c1
0276                = (n2c2/n1c1)*dot(E2_t, E2_t)
0277                = (n2c2/n1c1)*( ts*ts*E_perp*E_perp + tp*tp*(1-E1_perp*E1_perp) )
0278                = (n2c2/n1c1)*( ts*ts*E_s2 + tp*tp*(1-E_s2) )
0279
0280 **/
0281
0282
0283 int main(int argc, char** argv)
0284 {
0285     RNG rng ;
0286     rng.set_fake(0.);     // 0.:force transmit 1.:force reflect
0287
0288     float3 mom = normalize(make_float3(1.f, 0.f, -1.f)) ;
0289     float3 nrm = make_float3(0.f, 0.f, 1.f ); // surface normal in +z direction
0290
0291     quad2 prd ;
0292     prd.q0.f.x = nrm.x ;
0293     prd.q0.f.y = nrm.y ;
0294     prd.q0.f.z = nrm.z ;
0295
0296     sctx ctx ;
0297     ctx.prd = &prd ;
0298
0299     float n1 = 1.f ;
0300     float n2 = 1.5f ;
0301     float wl = 420.f ;
0302
0303     Custom c(nrm, mom, n1, n2, wl);
0304
0305     sstate& s = ctx.s ;
0306     s.material1.x = n1 ;
0307     s.material2.x = n2 ;
0308
0309     sphoton& p = ctx.p ;
0310
0311     const int N = 16 ;
0312     std::vector<sphoton> pp(N*2) ;
0313     std::vector<scf>    cfs(N);
0314
0315     for(int i=0 ; i < N ; i++)
0316     {
0317         sphoton& p0 = pp[i*2+0] ;
0318         sphoton& p1 = pp[i*2+1] ;
0319         scf& cf = cfs[i] ;
0320
0321         float frac_twopi = float(i)/float(N)  ;
0322
0323         p.zero();
0324         p.mom = mom ;
0325         p.set_polarization(frac_twopi) ;
0326         p0 = p ;
0327
0328         c.set_polarization(p.pol);
0329
0330
0331         sboundary b(rng, ctx);
0332         p1 = b.p ;
0333
0334         c.set_reflect(b.reflect) ;
0335
0336
0337         cf.frac_twopi = frac_twopi ;
0338
0339         cf.c_transverse_length = c.transverse_length ;
0340         cf.b_transverse_length = b.transverse_length ;
0341
0342         cf.c_T = c.T ;
0343         cf.b_TransCoeff = b.TransCoeff ;
0344
0345         cf.c_E1_perp = c.E1_perp ;
0346         cf.b_E1_perp = b.E1_perp ;
0347
0348
0349         std::cout
0350            << " b._E2_t " << b._E2_t
0351            << " length(b._E2_t) " << length(b._E2_t)
0352            << std::endl
0353            << " c._E2_t " << c._E2_t
0354            << " length(c._E2_t) " << length(c._E2_t)
0355            << std::endl
0356            << std::endl
0357            << " b._E2_r " << b._E2_r
0358            << " length(b._E2_r) " << length(b._E2_r)
0359            << std::endl
0360            << " c._E2_r " << c._E2_r
0361            << " length(c._E2_r) " << length(c._E2_r)
0362            << std::endl
0363            << std::endl
0364            << " b.E2_r "  << b.E2_r  << std::endl
0365            << " c.E2_r "  << c.E2_r  << std::endl
0366            << std::endl
0367            << " c.E2_t "  << c.E2_t  << std::endl
0368            << " b.E2_t "  << b.E2_t  << std::endl
0369            << std::endl
0370            << " b.E1   " << b.E1
0371            << " length(b.E1) " << length(b.E1)
0372            << std::endl
0373            << " c.E1   " << c.E1
0374            << " length(c.E1) " << length(c.E1)
0375            << std::endl
0376            << " b.E1_chk   " << b.E1_chk
0377            << std::endl
0378            << std::endl
0379            << " b.RR " << b.RR
0380            << " length(b.RR) " << length(b.RR)
0381            << std::endl
0382            << " c.RR " << c.RR
0383            << " length(c.RR) " << length(c.RR)
0384            << std::endl
0385            << std::endl
0386            << " b.TT   " << b.TT
0387            << " length(b.TT)   " << length(b.TT)
0388            << std::endl
0389            << " c.TT   " << c.TT
0390            << " length(c.TT)   " << length(c.TT)
0391            << std::endl
0392            << std::endl
0393            << " c.fT   " << c.fT  << std::endl
0394            << " c.fR   " << c.fR  << std::endl
0395            << std::endl
0396            << " c.E1E1   " << c.E1E1
0397            << " length(c.E1E1) " << length(c.E1E1)
0398            << std::endl
0399            << " c.E_s2 " << c.E_s2 << std::endl
0400            << " 1.f-c.E_s2 " << 1.f-c.E_s2 << std::endl
0401            << " c.T   " << c.T  << std::endl
0402            << " c.R   " << c.R  << std::endl
0403            << " c.A   " << c.A  << std::endl
0404            << " b.flag " << OpticksPhoton::Flag(b.flag)
0405            << std::endl
0406            << " p0 " << p0.descDir()
0407            << std::endl
0408            << " p1 " << p1.descDir()
0409            << std::endl
0410            << cf
0411            << std::endl
0412            << " c.new_mom " << c.new_mom << std::endl
0413            << " c.new_pol " << c.new_pol << std::endl
0414            << std::endl
0415            ;
0416
0417      }
0418
0419      //std::cout << " stack " << stack << std::endl ;
0420
0421      NP* a = NP::Make<float>(N,2,4,4) ;
0422      a->read2<float>( (float*)pp.data() );
0423      a->save(FOLD, "pp.npy");
0424      std::cout << " save to " << FOLD << "/pp.npy" << std::endl;
0425
0426      NP* a_cfs = NP::Make<float>(N, 4, 4 );
0427      a_cfs->read2<float>( (float*)cfs.data() );
0428      a_cfs->save(FOLD, "cf.npy");
0429      std::cout << " save to " << FOLD << "/cf.npy" << std::endl;
0430
0431
0432      return 0 ;
0433 }