sysrap/tests/sboundary_test.py

0001 #!/usr/bin/env python
0002 """
0003 sboundary_test.py
0004 ===================
0005
0006 Start from qudarap/tests/qsim_test.py
0007
0008 Explain those two, at Brewster angle::
0009
0010     In [5]: np.c_[pp[:,0,2],pp[:,1,2]]
0011     Out[5]:
0012     array([[-0.555,  0.   , -0.832,  0.   , -0.555, -0.   ,  0.832,  0.   ],
0013            [-0.552,  0.098, -0.828,  0.   , -0.   , -1.   ,  0.   ,  0.   ],
0014            [-0.544,  0.195, -0.816,  0.   , -0.   , -1.   ,  0.   ,  0.   ],
0015            ..
0016            [ 0.   ,  1.   ,  0.   ,  0.   , -0.   , -1.   ,  0.   ,  0.   ],
0017            [ 0.054,  0.995,  0.082,  0.   , -0.   , -1.   ,  0.   ,  0.   ],
0018            [ 0.552,  0.098,  0.828,  0.   , -0.   , -1.   ,  0.   ,  0.   ],
0019            [ 0.555, -0.   ,  0.832,  0.   , -0.483,  0.491,  0.725,  0.   ],
0020            [ 0.552, -0.098,  0.828,  0.   , -0.   ,  1.   ,  0.   ,  0.   ],
0021            [ 0.512, -0.383,  0.769,  0.   , -0.   ,  1.   ,  0.   ,  0.   ],
0022
0023
0024 Those are from initial pure P polarization that are forced to reflect
0025 which means that even when the coeff for reflection is zero for a particulr
0026 polarization are forced to calculate an outgoing polarization for it.
0027
0028 What happens in that situation with the calculation is that a
0029 near zero length E2_r (S,P) vector gets dignified by normalization
0030 into a vector with highly imprecise values that gets used
0031 to construct the new pol::
0032
0033               E2_r          (-0.000,-0.000)        length(E2_r)     0.0000
0034                  RR          (-0.491,-0.871)          length(RR)     1.0000
0035
0036
0037 This is not a problem with the calc, just with forcing the reflect/transmit.
0038
0039
0040 To avoid this messing up the illustration switched to scaling the the length
0041 of the polarization vector with the relevant coefficient,
0042 so these polarizations that would be quashed dissappear into the origin
0043
0044
0045 """
0046
0047 import os, numpy as np
0048 from opticks.ana.fold import Fold
0049 from opticks.ana.pvplt import *
0050
0051 try:
0052     import pyvista as pv
0053 except ImportError:
0054     pv = None
0055 pass
0056
0057 COLORS = "cyan red green blue cyan magenta yellow pink orange purple lightgreen".split()
0058
0059 import matplotlib as mp
0060 CMAP = os.environ.get("CMAP", "hsv")  # Reds
0061 cmap = mp.cm.get_cmap(CMAP)
0062
0063
0064 if __name__ == '__main__':
0065     t = Fold.Load(symbol="t")
0066     print(repr(t))
0067     pp = t.pp
0068
0069     polscale = float(os.environ.get("POLSCALE","1"))
0070     label = os.environ.get("TOPLINE", "sboundary_test.py")
0071     B = int(os.environ.get("B","1"))
0072
0073
0074     lim = slice(None)
0075
0076     mom0 = pp[:,0,1,:3]
0077     pol0 = pp[:,0,2,:3]
0078     mom00 = mom0[0]     # all mom0 expected to be the same, just varying pol
0079
0080     mom1 = pp[:,B,1,:3]
0081     pol1 = pp[:,B,2,:3]
0082     mom10 = mom1[0]     # all mom1 expected to be the same, just varying pol
0083
0084     pp[:,0,0,:3] = -mom0    # illustrative choice incident position on unit hemisphere
0085     pp[:,B,0,:3] =  mom1    # illustrative choice reflected/transmitted position on unit hemisphere
0086
0087     print(" mom1.shape %s " % (str(mom1.shape)))  #  (16, 3)
0088
0089     ppv = np.c_[pp[:,0,2],pp[:,B,2]]
0090     print("ppv = np.c_[pp[:,0,2],pp[:,B,2]] \n", repr(ppv))
0091
0092
0093
0094     pos = np.array( [[0,0,0]] , dtype=np.float32 )
0095     rvec = np.array( [[mom00[0],mom00[1],-mom00[2] ]], dtype=np.float32 )
0096     vec1 = np.array( [[mom10]], dtype=np.float32 )
0097
0098     lhs  = np.array( [[-1,0,0]], dtype=np.float32 )
0099     rhs  = np.array( [[ 1,0,0]], dtype=np.float32 )
0100
0101
0102
0103     if not pv is None:
0104         pl = pvplt_plotter(label=label)
0105         pvplt_viewpoint( pl )
0106
0107         N = len(pp)
0108         ii = list(range(N))
0109
0110         wscale = False
0111         wcut = True
0112
0113         for i in ii:
0114             frac = float(i)/float(N)
0115             polcol = cmap(frac)
0116             #polcol = COLORS[ i % len(COLORS)]
0117             pvplt_photon( pl, pp[i:i+1,0], polcol=polcol, polscale=polscale, wscale=wscale, wcut=wcut )
0118             pvplt_photon( pl, pp[i:i+1,B], polcol=polcol, polscale=polscale, wscale=wscale, wcut=wcut )
0119         pass
0120
0121
0122         pvplt_lines( pl, pos, rvec )
0123         pvplt_lines( pl, pos, vec1 )
0124
0125         pvplt_add_line_a2b( pl, lhs, rhs )
0126
0127         cp = pl.show()
0128
0129         print(cp)
0130
0131