Back to home page

EIC code displayed by LXR
 
 

    

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

0001 #pragma once
0002 /**
0003 stmm.h : Thin/Thick Multi-layer Stack TMM "Transfer Matrix Method" A,R,T calculation 
0004 ======================================================================================
0005 
0006 CAUTION : probably divergence from ~/j/Layr/JPMT.h 
0007 
0008 Developed in ~/j/Layr/Layr.h, moved into syrap to assist with comparison
0009 with qsim.h propagate_at_boundary 
0010 
0011 
0012 1. nothing JUNO specific here
0013 2. header-only implementation, installed with j/PMTSim 
0014 
0015 See Also
0016 -----------
0017 
0018 Layr.rst 
0019     notes/refs on TMM theory and CPU+GPU implementation 
0020 
0021 LayrTest.{h,cc,cu,py,sh} 
0022     build, run cpu+gpu scans, plotting, comparisons float/double cpu/gpu std/thrust 
0023 
0024 JPMT.h
0025     JUNO specific collection of PMT rindex and thickness into arrays 
0026 
0027 
0028 Contents of Layr.h : (persisted array shapes)
0029 ----------------------------------------------
0030 
0031 namespace Const 
0032     templated constexpr functions: zero, one, two, pi, twopi
0033 
0034 template<typename T> struct Matx : (4,2)
0035     2x2 complex matrix  
0036 
0037 template<typename T> struct Layr : (4,4,2)
0038     d:thickness, complex refractive index, angular and Fresnel coeff,  S+P Matx
0039 
0040     * d = zero : indicates thick (incoherent) layer 
0041 
0042 template<typename F> struct ART_ : (3,4) 
0043     results  
0044 
0045 template<typename T, int N> StackSpec :  (4,3) 
0046     4 sets of complex refractive index and thickness
0047 
0048     * HMM maybe pad to (4,4) if decide to keep ?
0049 
0050 template <typename T, int N> struct Stack : (constituents persisted separately) 
0051     N Layr stack : all calculations in ctor  
0052 
0053     Layr<T> ll[N] ;   
0054     Layr<T> comp ;  // composite for the N layers 
0055     ART_<T>  art ; 
0056 
0057     LAYR_METHOD Stack(T wl, T minus_cos_theta, const StackSpec4<T>& ss);
0058 
0059 **/
0060 
0061 #ifdef WITH_THRUST
0062 #include <thrust/complex.h>
0063 #else
0064 #include <complex>
0065 #include <cmath>
0066 #endif
0067 
0068 #if defined(__CUDACC__) || defined(__CUDABE__)
0069 #else
0070 #include <iostream>
0071 #include <iomanip>
0072 #include <sstream>
0073 #include <string>
0074 #include <cassert>
0075 #include <cstdlib>
0076 #include <vector>
0077 #endif
0078 
0079 #if defined(__CUDACC__) || defined(__CUDABE__)
0080 #    define LAYR_METHOD __host__ __device__ __forceinline__
0081 #else
0082 #    define LAYR_METHOD inline 
0083 #endif
0084 
0085 
0086 /**
0087 Const
0088 -------
0089 
0090 The point of these Const functions is to plant the desired type of constant
0091 into assembly code with no runtime transients of another type. Important for
0092 avoiding expensive unintentional doubles in GPU code. The constexpr means that
0093 the conversions and calculations happen at compile time, NOT runtime. 
0094 
0095 **/
0096 
0097 namespace Const
0098 {
0099     template<typename T>  
0100     LAYR_METHOD constexpr T zero(){ return T(0.0) ; } 
0101  
0102     template<typename T>
0103     LAYR_METHOD constexpr T one() { return T(1.0) ; } 
0104 
0105     template<typename T>
0106     LAYR_METHOD constexpr T two() { return T(2.0) ; } 
0107     
0108     template<typename T>
0109     LAYR_METHOD constexpr T pi() { return T(M_PI) ; } 
0110 
0111     template<typename T>
0112     LAYR_METHOD constexpr T twopi() { return T(2.0*M_PI) ; } 
0113 }
0114 
0115 
0116 #if defined(__CUDACC__) || defined(__CUDABE__)
0117 #else
0118     #ifdef WITH_THRUST
0119     template<typename T>
0120     inline std::ostream& operator<<(std::ostream& os, const thrust::complex<T>& z)  
0121     {
0122         os << "(" << std::setw(10) << std::fixed << std::setprecision(4) << z.real() 
0123            << " " << std::setw(10) << std::fixed << std::setprecision(4) << z.imag() << ")t" ; 
0124         return os; 
0125     }
0126     #else
0127     template<typename T>
0128     inline std::ostream& operator<<(std::ostream& os, const std::complex<T>& z)  
0129     {
0130         os << "(" << std::setw(10) << std::fixed << std::setprecision(4) << z.real() 
0131            << " " << std::setw(10) << std::fixed << std::setprecision(4) << z.imag() << ")s" ; 
0132         return os; 
0133     }
0134     #endif  // clarity is more important than brevity 
0135 #endif
0136 
0137 template<typename T>
0138 struct Matx
0139 {
0140 #ifdef WITH_THRUST
0141     thrust::complex<T> M00, M01, M10, M11 ;   
0142 #else
0143     std::complex<T>    M00, M01, M10, M11 ;       
0144 #endif
0145     LAYR_METHOD void reset();             
0146     LAYR_METHOD void dot(const Matx<T>& other); 
0147 };
0148 
0149 template<typename T>
0150 LAYR_METHOD void Matx<T>::reset()
0151 {
0152     M00.real(1) ; M00.imag(0) ; // conversion from int 
0153     M01.real(0) ; M01.imag(0) ; 
0154     M10.real(0) ; M10.imag(0) ; 
0155     M11.real(1) ; M11.imag(0) ; 
0156 }
0157 /**
0158 
0159       | T00  T01  |  |  M00   M01 | 
0160       |           |  |            | 
0161       | T10  T11  |  |  M10   M11 | 
0162 
0163 **/
0164 template<typename T>
0165 LAYR_METHOD void Matx<T>::dot(const Matx<T>& other)
0166 {
0167 #ifdef WITH_THRUST
0168     using thrust::complex ; 
0169 #else
0170     using std::complex ; 
0171 #endif
0172     complex<T> T00(M00) ; 
0173     complex<T> T01(M01) ; 
0174     complex<T> T10(M10) ; 
0175     complex<T> T11(M11) ; 
0176 
0177     M00 = T00*other.M00 + T01*other.M10;
0178     M01 = T00*other.M01 + T01*other.M11;
0179     M10 = T10*other.M00 + T11*other.M10;
0180     M11 = T10*other.M01 + T11*other.M11;
0181 }
0182 
0183 #if defined(__CUDACC__) || defined(__CUDABE__)
0184 #else
0185 template<typename T>
0186 inline std::ostream& operator<<(std::ostream& os, const Matx<T>& m)  
0187 {
0188     os 
0189        << "| " << m.M00 << " " << m.M01 << " |" << std::endl 
0190        << "| " << m.M10 << " " << m.M11 << " |" << std::endl  
0191        ;
0192     return os; 
0193 }
0194 #endif
0195 
0196 
0197 /**
0198 Layr : (4,4,2) 
0199 -----------------
0200 
0201 The comp layers do not have the 0th (4,2) filled::
0202 
0203    assert np.all( e.f.comps[:,0] == 0 ) 
0204 
0205 **/
0206 
0207 template<typename T>
0208 struct Layr
0209 {
0210     // ---------------------------------------- 0th (4,2)
0211     T  d ;
0212     T  pad=0 ;
0213 #ifdef WITH_THRUST 
0214     thrust::complex<T>  n, st, ct ; 
0215 #else
0216     std::complex<T>     n, st, ct ;
0217 #endif
0218     // ---------------------------------------- 1st (4,2)
0219 
0220 #ifdef WITH_THRUST 
0221     thrust::complex<T>  rs, rp, ts, tp ;    
0222 #else
0223     std::complex<T>     rs, rp, ts, tp ;    
0224 #endif
0225     // ---------------------------------------- 2nd (4,2)
0226     Matx<T> S ;                             
0227     // ---------------------------------------- 3rd (4,2)
0228     Matx<T> P ;                               
0229     // ---------------------------------------- 
0230 
0231     LAYR_METHOD void reset(); 
0232     LAYR_METHOD void load4( const T* vals ); 
0233 };
0234 
0235 template<typename T>
0236 LAYR_METHOD void Layr<T>::reset()
0237 {
0238     S.reset(); 
0239     P.reset(); 
0240 }
0241 
0242 template<typename T>
0243 LAYR_METHOD void Layr<T>::load4(const T* vals)
0244 {
0245     d   = vals[0] ; 
0246     pad = vals[1] ; 
0247     n.real(vals[2]) ; 
0248     n.imag(vals[3]) ; 
0249 }
0250 
0251 
0252 
0253 #if defined(__CUDACC__) || defined(__CUDABE__)
0254 #else
0255 template<typename T>
0256 inline std::ostream& operator<<(std::ostream& os, const Layr<T>& l)  
0257 {
0258     os 
0259        << "Layr"
0260        << std::endl 
0261        << std::setw(4) << "n:" << l.n  
0262        << std::setw(4) << "d:" << std::fixed << std::setw(10) << std::setprecision(4) << l.d  
0263        << std::endl 
0264        << std::setw(4) << "st:" << l.st 
0265        << std::setw(4) << "ct:" << l.ct
0266        << std::endl 
0267        << std::setw(4) << "rs:" << l.rs 
0268        << std::setw(4) << "rp:" << l.rp
0269        << std::endl 
0270        << std::setw(4) << "ts:" << l.ts 
0271        << std::setw(4) << "tp:" << l.tp
0272        << std::endl 
0273        << "S" 
0274        << std::endl 
0275        << l.S 
0276        << std::endl 
0277        << "P"
0278        << std::endl 
0279        << l.P
0280        << std::endl 
0281        ;
0282     return os; 
0283 }
0284 #endif
0285 
0286 template<typename F>
0287 struct ART_
0288 {   
0289     F R_s;     // R_s = a.arts[:,0,0]
0290     F R_p;     // R_p = a.arts[:,0,1]
0291     F T_s;     // T_s = a.arts[:,0,2]
0292     F T_p;     // T_p = a.arts[:,0,3]
0293 
0294     F A_s;     // A_s = a.arts[:,1,0]
0295     F A_p;     // A_p = a.arts[:,1,1]
0296     F R;       // R   = a.arts[:,1,2]
0297     F T;       // T   = a.arts[:,1,3]
0298 
0299     F A;       // A   = a.arts[:,2,0]
0300     F A_R_T ;  // A_R_T = a.arts[:,2,1] 
0301     F wl ;     // wl  = a.arts[:,2,2]
0302     F mct ;    // mct  = a.arts[:,2,3]   
0303 
0304     // persisted into shape (3,4) 
0305 };
0306 
0307 #if defined(__CUDACC__) || defined(__CUDABE__)
0308 #else
0309 template<typename T>
0310 inline std::ostream& operator<<(std::ostream& os, const ART_<T>& art )  
0311 {
0312     os 
0313         << "ART_" << std::endl 
0314         << " R_s " << std::fixed << std::setw(10) << std::setprecision(4) << art.R_s 
0315         << " R_p " << std::fixed << std::setw(10) << std::setprecision(4) << art.R_p 
0316         << std::endl 
0317         << " T_s " << std::fixed << std::setw(10) << std::setprecision(4) << art.T_s 
0318         << " T_p " << std::fixed << std::setw(10) << std::setprecision(4) << art.T_p 
0319         << std::endl 
0320         << " A_s " << std::fixed << std::setw(10) << std::setprecision(4) << art.A_s 
0321         << " A_p " << std::fixed << std::setw(10) << std::setprecision(4) << art.A_p 
0322         << std::endl 
0323         << " R   " << std::fixed << std::setw(10) << std::setprecision(4) << art.R   
0324         << " T   " << std::fixed << std::setw(10) << std::setprecision(4) << art.T   
0325         << " A   " << std::fixed << std::setw(10) << std::setprecision(4) << art.A  
0326         << " A_R_T " << std::fixed << std::setw(10) << std::setprecision(4) << art.A_R_T 
0327         << std::endl 
0328         << " wl  " << std::fixed << std::setw(10) << std::setprecision(4) << art.wl  << std::endl 
0329         << " mct " << std::fixed << std::setw(10) << std::setprecision(4) << art.mct << std::endl 
0330         ;
0331     return os; 
0332 }
0333 #endif
0334 
0335 
0336 #if defined(__CUDACC__) || defined(__CUDABE__)
0337 #else
0338 namespace sys
0339 {
0340     template<typename T>
0341     inline std::vector<T>* getenvvec(const char* ekey, const char* fallback = nullptr, char delim=',')
0342     {
0343         char* _line = getenv(ekey);
0344         const char* line = _line ? _line : fallback ; 
0345         if(line == nullptr) return nullptr ; 
0346 
0347         std::stringstream ss; 
0348         ss.str(line);
0349         std::string s;
0350 
0351         std::vector<T>* vec = new std::vector<T>() ; 
0352 
0353         while (std::getline(ss, s, delim)) 
0354         {   
0355             std::istringstream iss(s);
0356             T t ; 
0357             iss >> t ; 
0358             vec->push_back(t) ; 
0359         }   
0360         return vec ; 
0361     }
0362 }
0363 #endif
0364 
0365 
0366 template<typename T>
0367 struct LayrSpec
0368 {
0369     T nr, ni, d ; 
0370 #if defined(__CUDACC__) || defined(__CUDABE__)
0371 #else
0372     static int EGet(LayrSpec<T>& ls, int idx); 
0373 #endif
0374 };
0375 
0376 #if defined(__CUDACC__) || defined(__CUDABE__)
0377 #else
0378 template<typename T>
0379 LAYR_METHOD int LayrSpec<T>::EGet(LayrSpec<T>& ls, int idx)
0380 {
0381     std::stringstream ss ; 
0382     ss << "L" << idx ; 
0383     std::string ekey = ss.str(); 
0384     std::vector<T>* vls = sys::getenvvec<T>(ekey.c_str()) ; 
0385     if(vls == nullptr) return 0 ; 
0386     const T zero(0) ; 
0387     ls.nr = vls->size() > 0u ? (*vls)[0] : zero ; 
0388     ls.ni = vls->size() > 1u ? (*vls)[1] : zero ; 
0389     ls.d  = vls->size() > 2u ? (*vls)[2] : zero ; 
0390     return 1 ; 
0391 }
0392 
0393 template<typename T>
0394 LAYR_METHOD std::ostream& operator<<(std::ostream& os, const LayrSpec<T>& ls )  
0395 {
0396     os 
0397         << "LayrSpec<" << ( sizeof(T) == 8 ? "double" : "float" ) << "> "  
0398         << std::fixed << std::setw(10) << std::setprecision(4) << ls.nr << " "
0399         << std::fixed << std::setw(10) << std::setprecision(4) << ls.ni << " ; "
0400         << std::fixed << std::setw(10) << std::setprecision(4) << ls.d  << ")"
0401         << std::endl 
0402         ;
0403     return os ; 
0404 }
0405 #endif
0406 
0407 
0408 template<typename T, int N>  
0409 struct StackSpec
0410 {
0411     LayrSpec<T> ls[N] ; 
0412 
0413 #if defined(__CUDACC__) || defined(__CUDABE__)
0414 #else
0415     LAYR_METHOD static StackSpec<T,N> Create2(float n1, float n2) ; 
0416     LAYR_METHOD void eget() ; 
0417 #endif
0418 
0419 }; 
0420 
0421 #if defined(__CUDACC__) || defined(__CUDABE__)
0422 #else
0423 
0424 template<typename T, int N>
0425 LAYR_METHOD StackSpec<T,N> StackSpec<T,N>::Create2(float n1, float n2)
0426 {
0427     assert( N == 2 ); 
0428     StackSpec<T,N> ss ; 
0429 
0430     ss.ls[0].nr = n1 ; 
0431     ss.ls[0].ni = 0. ;
0432     ss.ls[0].d  = 0. ;
0433 
0434     ss.ls[1].nr = n2 ; 
0435     ss.ls[1].ni = 0. ;
0436     ss.ls[1].d  = 0. ;
0437 
0438     return ss ; 
0439 }
0440 
0441 
0442 template<typename T, int N>
0443 LAYR_METHOD void StackSpec<T,N>::eget()  
0444 {
0445     int count = 0 ; 
0446     for(int i=0 ; i < N ; i++) count += LayrSpec<T>::EGet(ls[i], i); 
0447     assert( count == N ) ; 
0448 }
0449 
0450 template<typename T, int N>
0451 LAYR_METHOD std::ostream& operator<<(std::ostream& os, const StackSpec<T,N>& ss )  
0452 {
0453     os 
0454         << "StackSpec<" 
0455         << ( sizeof(T) == 8 ? "double" : "float" ) 
0456         << "," << N
0457         << ">"  
0458         << std::endl ;
0459 
0460     for(int i=0 ; i < N ; i++) os << ss.ls[i] ; 
0461     return os ; 
0462 }
0463 
0464 #endif
0465 
0466 
0467 /**
0468 Stack
0469 -------
0470 
0471 **/
0472 
0473 template <typename T, int N>
0474 struct Stack
0475 {
0476     Layr<T> ll[N] ;  
0477     Layr<T> comp ;  // composite for the N layers 
0478     ART_<T>  art ; 
0479 
0480     LAYR_METHOD void zero();
0481     LAYR_METHOD Stack();
0482     LAYR_METHOD Stack(T wl, T minus_cos_theta, const StackSpec<T,N>& ss);
0483 };
0484 
0485 template<typename T, int N>
0486 LAYR_METHOD void Stack<T,N>::zero()
0487 {
0488     art = {} ; 
0489     comp = {} ; 
0490     for(int i=0 ; i < N ; i++) ll[i] = {} ; 
0491 }
0492 
0493 
0494 /**
0495 Stack::Stack
0496 ---------------
0497 
0498 Caution that StackSpec contains refractive indices that depend on wavelength, 
0499 so the wavelength dependency enters twice. 
0500 
0501 SO instead pass in a reference to the object "QPMT.hh/spmt.h" 
0502 that handles the PMT properties, and is responsible for:
0503 
0504 1. holding PMT properties (or textures)
0505 2. doing property lookup as function of wavelenth/energy
0506 3. populating Layr<T>* with the results 
0507 
0508 HMM: 
0509 
0510 * how to test the counterpair pair in a way that can work on both CPU and GPU ?
0511 * dont need to use NP::combined_interpolate_5 GPU compatible lookup code can be used on CPU also 
0512 
0513 
0514 HMM: more physical to use dot(photon_momentum,outward_surface_normal) 
0515 as "angle" parameter, the dot product is -cos(aoi)
0516 
0517 1. -1 at normal incidence against surface_normal, inwards going 
0518 2. +1 at normal incidence with the surface_normal, outwards going  
0519 3.  0 at glancing incidence (90 deg AOI) : potential for math blowouts here 
0520 4. sign of dot product indicates when must flip the stack of parameters
0521 5. angle scan plots can then use aoi 0->180 deg, which is -cos(aoi) -1->1   
0522    (will there be continuity across the turnaround ?)
0523 
0524 **/
0525 
0526 template<typename T, int N>
0527 LAYR_METHOD Stack<T,N>::Stack()
0528 {
0529     zero(); 
0530 }
0531 
0532 template<typename T, int N>
0533 LAYR_METHOD Stack<T,N>::Stack(T wl, T minus_cos_theta, const StackSpec<T,N>& ss ) 
0534 {
0535     // minus_cos_theta, aka dot(mom,normal)
0536 #ifdef WITH_THRUST
0537     using thrust::complex ; 
0538     using thrust::norm ; 
0539     using thrust::conj ; 
0540     using thrust::exp ; 
0541     using thrust::sqrt ; 
0542     using thrust::sin ; 
0543     using thrust::cos ; 
0544 #else
0545     using std::complex ; 
0546     using std::norm ; 
0547     using std::conj ; 
0548     using std::exp ; 
0549     using std::sqrt ; 
0550     using std::sin ; 
0551     using std::cos ; 
0552 #endif
0553 
0554 #if defined(__CUDACC__) || defined(__CUDABE__)
0555 #else
0556     assert( N >= 2); 
0557 #endif
0558 
0559     const T zero(Const::zero<T>()) ; 
0560     const T one(Const::one<T>()) ; 
0561     const T two(Const::two<T>()) ; 
0562     const T twopi(Const::twopi<T>()) ; 
0563 
0564     for(int i=0 ; i < N ; i++)
0565     {
0566         int j = minus_cos_theta < zero ? i : N - 1 - i ;  
0567         //  minus_cos_theta < zero  : against normal : ordinary stack  : j = i 
0568         //  minus_cos_theta >= zero : with normal    : backwards stack : j from end 
0569 
0570         ll[j].n.real(ss.ls[i].nr) ; 
0571         ll[j].n.imag(ss.ls[i].ni) ; 
0572         ll[j].d = ss.ls[i].d ; 
0573     }
0574 
0575     // ll[0]   is "top"     : start layer : incident
0576     // ll[N-1] is "bottom"  : end   layer : transmitted
0577 
0578 
0579     art.wl = wl ; 
0580     art.mct = minus_cos_theta ; 
0581 
0582     const complex<T> zOne(one,zero); 
0583     const complex<T> zI(zero,one); 
0584     const complex<T> mct(minus_cos_theta);  // simpler for everything to be complex
0585 
0586     // Snell : set st,ct of all layers (depending on indices(hence wl) and incident angle) 
0587     Layr<T>& l0 = ll[0] ; 
0588     l0.ct = minus_cos_theta < zero  ? -mct : mct ; 
0589     //
0590     //  flip picks +ve ct that constrains the angle to first quadrant 
0591     //  this works as : cos(pi-theta) = -cos(theta)
0592     //  without flip, the ART values are non-physical : always outside 0->1 for mct > 0 angle > 90  
0593     //
0594 
0595     l0.st = sqrt( zOne - mct*mct ) ;  
0596 
0597     for(int idx=1 ; idx < N ; idx++)  // for N=2 idx only 1, sets only ll[1] 
0598     {
0599         Layr<T>& l = ll[idx] ; 
0600         l.st = l0.n * l0.st / l.n  ; 
0601         l.ct = sqrt( zOne - l.st*l.st );
0602     }     
0603 
0604     // Fresnel : set rs/rp/ts/tp for N-1 interfaces between N layers
0605     // (depending on indices(hence wl) and layer ct) 
0606     // HMM: last layer unset, perhaps zero it ?
0607     for(int idx=0 ; idx < N-1 ; idx++)
0608     {
0609         // cf OpticalSystem::Calculate_rt  
0610         // https://en.wikipedia.org/wiki/Fresnel_equations
0611         Layr<T>& i = ll[idx] ; 
0612         const Layr<T>& j = ll[idx+1] ;  
0613 
0614         i.rs = (i.n*i.ct - j.n*j.ct)/(i.n*i.ct+j.n*j.ct) ;  // r_s eoe[12] see g4op-eoe
0615         i.rp = (j.n*i.ct - i.n*j.ct)/(j.n*i.ct+i.n*j.ct) ;  // r_p eoe[7]
0616         i.ts = (two*i.n*i.ct)/(i.n*i.ct + j.n*j.ct) ;       // t_s eoe[9]
0617         i.tp = (two*i.n*i.ct)/(j.n*i.ct + i.n*j.ct) ;       // t_p eoe[8]
0618     }
0619 
0620     // populate transfer matrices for both thick and thin layers  
0621     // for N=2 only one interface
0622 
0623     ll[0].reset();    // ll[0].S ll[0].P matrices set to identity 
0624 
0625     for(int idx=1 ; idx < N ; idx++)
0626     {
0627         const Layr<T>& i = ll[idx-1] ;            
0628         Layr<T>& j       = ll[idx] ;          
0629         // looking at (i,j) pairs of layers 
0630 
0631         complex<T> tmp_s = one/i.ts ; 
0632         complex<T> tmp_p = one/i.tp ;   
0633         // at glancing incidence ts, tp approach zero : blowing up tmp_s tmp_p
0634         // which causes the S and P matrices to blow up yielding infinities at mct zero
0635         //
0636         // thick layers indicated with d = 0. 
0637         // thin layers have thickness presumably comparable to art.wl (WITH SAME LENGTH UNIT: nm)
0638         complex<T> delta         = j.d == zero ? zero : twopi*j.n*j.d*j.ct/art.wl ; 
0639         complex<T> exp_neg_delta = j.d == zero ? one  : exp(-zI*delta) ; 
0640         complex<T> exp_pos_delta = j.d == zero ? one  : exp( zI*delta) ; 
0641 
0642         j.S.M00 = tmp_s*exp_neg_delta      ; j.S.M01 = tmp_s*i.rs*exp_pos_delta ; 
0643         j.S.M10 = tmp_s*i.rs*exp_neg_delta ; j.S.M11 =      tmp_s*exp_pos_delta ; 
0644 
0645         j.P.M00 = tmp_p*exp_neg_delta      ; j.P.M01 = tmp_p*i.rp*exp_pos_delta ; 
0646         j.P.M10 = tmp_p*i.rp*exp_neg_delta ; j.P.M11 =      tmp_p*exp_pos_delta ; 
0647 
0648         // NB: for thin layers the transfer matrices combine interface and propagation between them 
0649     }
0650 
0651     /*
0652     For d=0 (thick layer) 
0653 
0654 
0655     S     |  1/ts    rs/ts |    
0656           |                |
0657           |  rs/ts   1/ts  |     
0658 
0659 
0660     P     |  1/tp    rp/tp |    
0661           |                |
0662           |  rp/tp   1/tp  |     
0663 
0664     */
0665 
0666 
0667 
0668     // product of the transfer matrices
0669     comp.d = zero ; 
0670     comp.st = zero ; 
0671     comp.ct = zero ; 
0672     comp.S.reset(); 
0673     comp.P.reset(); 
0674 
0675     for(int idx=0 ; idx < N ; idx++) // TODO: start from idx=1 as ll[0].S ll[0].P always identity
0676     {
0677         const Layr<T>& l = ll[idx] ; 
0678         comp.S.dot(l.S) ; 
0679         comp.P.dot(l.P) ; 
0680     }
0681     // at glancing incidence the matrix from 
0682     // one of the layers has infinities, which 
0683     // yields nan in the matrix product 
0684     // and yields nan for all the below Fresnel coeffs 
0685     //
0686     // extract amplitude factors from the composite matrix
0687     comp.rs = comp.S.M10/comp.S.M00 ; 
0688     comp.rp = comp.P.M10/comp.P.M00 ; 
0689     comp.ts = one/comp.S.M00 ; 
0690     comp.tp = one/comp.P.M00 ; 
0691 
0692     Layr<T>& t = ll[0] ; 
0693     Layr<T>& b = ll[N-1] ; 
0694 
0695     // getting from amplitude to power relations for relectance R (same material and angle) and tranmittance T
0696     // https://www.brown.edu/research/labs/mittleman/sites/brown.edu.research.labs.mittleman/files/uploads/lecture13_0.pdf
0697 
0698     complex<T> _T_s = (b.n*b.ct)/(t.n*t.ct)*norm(comp.ts) ;  // aka n2c2/n1c1*ts*ts 
0699     complex<T> _T_p = (b.n*b.ct)/(t.n*t.ct)*norm(comp.tp) ;  // aka n2c2/n1c1*tp*tp
0700 
0701     art.T_s = _T_s.real() ; 
0702     art.T_p = _T_p.real() ; 
0703     art.R_s = norm(comp.rs) ;    // norm is squared magnitude, ie sum of squares of real and imag parts, so for real its just the square
0704     art.R_p = norm(comp.rp) ; 
0705 
0706     art.A_s = one-art.R_s-art.T_s;   // absorption factor by subtracting reflection and transmission
0707     art.A_p = one-art.R_p-art.T_p;
0708 
0709     //
0710     // complex<T> _T_p = (conj(b.n)*b.ct)/(conj(t.n)*t.ct)*norm(comp.tp) ; 
0711     // top and bot layers assumed to always have real index only, so remove the conj 
0712     //
0713     //
0714     // average of S and P : actually not very useful 
0715     art.R   = (art.R_s+art.R_p)/two ;
0716     art.T   = (art.T_s+art.T_p)/two ;
0717     art.A   = (art.A_s+art.A_p)/two ;
0718     art.A_R_T = art.A + art.R + art.T ;  
0719 }
0720 
0721 #if defined(__CUDACC__) || defined(__CUDABE__)
0722 #else
0723 template <typename T, int N>
0724 inline std::ostream& operator<<(std::ostream& os, const Stack<T,N>& stk )  
0725 {
0726     os << "Stack"
0727        << "<" 
0728        << ( sizeof(T) == 8 ? "double" : "float" )
0729        << ","
0730        << N 
0731        << ">" 
0732        << std::endl
0733        ; 
0734     for(int idx=0 ; idx < N ; idx++) os << "idx " << idx << std::endl << stk.ll[idx] ; 
0735     os << "comp" 
0736        << std::endl 
0737        << stk.comp 
0738        << std::endl 
0739        << stk.art 
0740        << std::endl 
0741        ; 
0742     return os; 
0743 }
0744 #endif
0745
[ 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!