sysrap/tests/sboundary_test.cc

0001 /**
0002 sboundary_test.cc
0003 ====================
0004
0005 Build and run::
0006
0007     ~/o/sysrap/tests/sboundary_test.sh
0008
0009
0010
0011            +-st-+
0012             \   :
0013              \  ct
0014               \ :
0015                \:
0016          -------+--------
0017
0018 **/
0019
0020 #include "srngcpu.h"
0021 using RNG = srngcpu ;
0022
0023 #include "sphoton.h"
0024 #include "sstate.h"
0025 #include "srec.h"
0026 #include "sseq.h"
0027 #include "stag.h"
0028 #include "sevent.h"
0029 #include "sctx.h"
0030
0031 #include "sboundary.h"
0032
0033 const char* FOLD = getenv("FOLD") ;
0034
0035
0036
0037 int main(int argc, char** argv)
0038 {
0039     RNG rng ;
0040
0041     float3 nrm = make_float3(0.f, 0.f, 1.f ); // surface normal in +z direction
0042
0043     const int N = U::GetEnvInt("N",16) ;
0044     float n1 = U::GetE<float>("N1", 1.f) ;
0045     float n2 = U::GetE<float>("N2", 1.5f) ;
0046     const char* AOI = U::GetEnv("AOI", "45" );
0047     const char force = U::GetE<char>("FORCE", 'N') ;
0048     switch(force)
0049     {
0050         case 'R':rng.set_fake(1.f) ; break ;
0051         case 'T':rng.set_fake(0.f) ; break ;
0052     }
0053
0054     float aoi = 0.f ;
0055     if(strcmp(AOI, "BREWSTER") == 0)
0056     {
0057         aoi =  atanf(n2/n1) ;
0058     }
0059     else if(strcmp(AOI, "CRITICAL") == 0)
0060     {
0061         assert( n2/n1 <= 1.f );
0062         aoi = asinf(n2/n1);
0063     }
0064     else
0065     {
0066         aoi = std::atof(AOI)*M_PIf/180.f ;
0067     }
0068     float3 mom = normalize(make_float3(sinf(aoi), 0.f, -cosf(aoi))) ;
0069
0070
0071     std::cout
0072         << " N " << N
0073         << " N1 " << n1
0074         << " N2 " << n2
0075         << " AOI " << AOI
0076         << " FORCE " << force
0077         << std::endl
0078         << " aoi " << aoi
0079         << " aoi/M_PIf " << aoi/M_PIf
0080         << " aoi/M_PIf*180 " << aoi/M_PIf*180.f
0081         << " mom " << mom
0082         << std::endl
0083         ;
0084
0085
0086     quad2 prd ;
0087     prd.q0.f.x = nrm.x ;
0088     prd.q0.f.y = nrm.y ;
0089     prd.q0.f.z = nrm.z ;
0090
0091     sctx ctx ;
0092     ctx.prd = &prd ;
0093
0094     sstate& s = ctx.s ;
0095     s.material1.x = n1 ;
0096     s.material2.x = n2 ;
0097
0098     sphoton& p = ctx.p ;
0099
0100     const int pnum = 3 ;
0101     std::vector<sphoton> pp(N*pnum) ;
0102
0103     for(int i=0 ; i < N ; i++)  // loop over N photons with different polarization directions
0104     {
0105         sphoton& p0 = pp[i*3+0] ;
0106         sphoton& p1 = pp[i*3+1] ;
0107         sphoton& p2 = pp[i*3+2] ;
0108
0109         float frac_twopi = float(i)/float(N); // does not reach 1. to avoid including both 0 and 2pi
0110
0111         p.zero();
0112
0113         p.mom = mom ;                     // all N with same momentum direction
0114         p.set_polarization(frac_twopi) ;
0115         p.time = frac_twopi ;
0116
0117         p0 = p ;
0118         p0.wavelength = 0.1f ;  // used for scaling the pol vector
0119
0120         sboundary b(rng, ctx); // changes ctx.p, p, b.p  (all the same reference)
0121
0122         p1 = b.p ;
0123         p1.wavelength = b.Coeff ;
0124
0125         p2 = p1 ;
0126         p2.pol = b.alt_pol ;   // p2: same as p1 but with alt_pol
0127
0128         std::cout
0129             << " b.flag " << OpticksPhoton::Flag(b.flag)
0130             << std::endl
0131             << " p0 " << p0.descDir()
0132             << std::endl
0133             << " p1 " << p1.descDir()
0134             << std::endl
0135             << " p2 " << p2.descDir()
0136             << std::endl
0137             ;
0138
0139         std::cout << b ;
0140
0141      }
0142
0143      NP* a = NP::Make<float>(N,pnum,4,4) ;
0144      a->read2<float>( (float*)pp.data() );
0145      a->save(FOLD, "pp.npy");
0146      std::cout << " save to " << FOLD << "/pp.npy" << std::endl;
0147
0148      return 0 ;
0149 }