eic-opticks/ana/tprism.py

0001 #!/usr/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022 tprism.py : Comparison of simulation with analytic expectation
0023 ================================================================
0024
0025 Analysis of "prism" and "newton" event categories
0026
0027 * "prism" uses all incident angles, see `tprism-`
0028 * "newton" uses one incident angle, see `tnewton-`
0029
0030
0031 TODO:
0032
0033 * handle multiple wavelengths
0034 * make comparison after histogramming, like reflection.py does
0035 * but here have incident angle as well as the deviation, so need 2d
0036 * http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.histogram2d.html
0037
0038
0039
0040 Prism Deviation Angle Calculation
0041 ------------------------------------
0042
0043 Hecht p163, two refractions thru a prism, CB and CD are normal to surface planes::
0044
0045    .
0046
0047                     A
0048                    / \
0049                   /   \
0050                  /     \
0051                 /       \
0052                /         \
0053               /           \
0054              /             \
0055             B. . . . . . . .D
0056        .   /                 \
0057           /          C        \
0058          /                     \
0059         /                       \
0060        --------------------------
0061
0062
0063 Polygon ABCD
0064
0065    BAD : alpha   (apex angle)
0066    CBA : 90 degrees
0067    CDA : 90 degrees
0068    BCD : 180 - alpha
0069
0070 Triangle BCD has angles:
0071
0072    CBD: t1
0073    CDB: i2
0074    BCD: 180 - alpha
0075
0076   ==>  alpha = t1 + i2
0077
0078
0079 Ray path thru prism, BD
0080
0081 Deviation at B:   (i1 - t1)
0082 Deviation at D:   (t2 - i2)
0083 Total
0084          delta:   (i1 - t1) + (t2 - i2)
0085
0086          delta = i1 + t2 - alpha
0087
0088 Aiming for expression providing delta as function of theta_i1,
0089 apex angle alpha and refractive index n
0090
0091 Snells law at 2nd interface, prism refractive index n in air
0092
0093    sin(t2) = n sin(i2)
0094
0095        t2  = arcsin( n sin(i2) )
0096            = arcsin( n sin(alpha - t1) )
0097            = arcsin( sin(alpha) n cos(t1) - cos(alpha) n sin(t1) )
0098
0099
0100       [ n cos(t1) ]^2  = (n^2 - n^2 sin(t1)^2 )
0101
0102                        =  (n^2 - sin(i1)^2 )           n sin(t1) = sin(i1)
0103
0104
0105        t2 = arcsin(  sin(alpha) sqrt(n^2 - sin^2(i1)) - sin(i1) cos(alpha) )
0106
0107
0108 2nd refraction has a critical angle where t2 = pi above which TIR will occur
0109
0110        n sin(i2) = 1*sin(t2) =  1
0111
0112              i2c = arcsin( 1./n )
0113
0114
0115 Propagate that back to 1st refraction
0116
0117          sin(i1) = n sin(t1) = n sin(alpha - i2)
0118
0119              i1  = arcsin( n sin(alpha - arcsin(1/n) ) )
0120
0121
0122 But there is another constraint in there with edge
0123
0124              n sin(alpha - arcsin(1/n)) = 1
0125                    alpha - arcsin(1/n) = arcsin(1/n)
0126
0127                            alpha/2 = arcsin(1/n) = i2c
0128
0129                            alpha = 2*i2c      (i2_c  43.304 n = 1.458)
0130
0131 This indicates that a 90 degree apex angle is not a good choice
0132 for dispersing prism, use 60 degree instead.
0133
0134
0135
0136 At minimum deviation delta, ray are parallel to base and have symmetric ray
0137
0138     i1 = t2
0139     t1 = i2
0140
0141     alpha = t1 + i2         ==>  t1 = i2 = alpha/2
0142
0143     delta = i1 + t2 - alpha
0144
0145
0146     sin(delta + alpha) = sin( i1 + t2 )
0147
0148     sin(i1) = n sin(t1)
0149
0150         i1 = arcsin( n sin(alpha/2) )
0151
0152
0153 Where to shoot from to get minimum deviation ?
0154
0155 * Use intersect frame coordinate with the transform explicitly specifified
0156
0157
0158 """
0159
0160 import os, sys, logging, numpy as np
0161 log = logging.getLogger(__name__)
0162
0163 import matplotlib.pyplot as plt
0164 from mpl_toolkits.mplot3d import Axes3D
0165
0166 from opticks.ana.base import opticks_main
0167 from opticks.ana.nbase import count_unique
0168 from opticks.ana.evt import Evt, costheta_
0169 from opticks.ana.ana import Rat, theta
0170 from opticks.ana.geometry import Shape, Plane, Boundary, Ray, Intersect, IntersectFrame, mat4_tostring, mat4_fromstring
0171
0172
0173 rad = np.pi/180.
0174 deg = 180./np.pi
0175
0176
0177 np.set_printoptions(suppress=True, precision=3)
0178 np.seterr(divide="ignore", invalid="ignore")
0179
0180 rat_ = lambda n,d:float(len(n))/float(len(d))
0181
0182 X,Y,Z,W = 0,1,2,3
0183
0184
0185 class Box(Shape):
0186     def __init__(self, parameters, boundary, wavelength=380 ):
0187         Shape.__init__(self, parameters, boundary)
0188
0189
0190 class Prism(Shape):
0191     def __init__(self, parameters, boundary):
0192         Shape.__init__(self, parameters, boundary)
0193
0194         alpha = self.parameters[0]
0195         height = self.parameters[1]
0196         depth = self.parameters[2]
0197
0198         a  = alpha*np.pi/180.
0199         pass
0200         self.a = a
0201         self.sa = np.sin(a)
0202         self.ca = np.cos(a)
0203
0204         self.alpha = alpha
0205
0206         ymax =  height/2.
0207         ymin = -height/2.
0208
0209         hwidth = height*np.tan(a/2.)
0210         apex = [0,ymax,0]
0211         base = [0,ymin,0]
0212         front = [0,ymin,depth/2.]
0213         back  = [0,ymin,-depth/2.]
0214
0215         lhs = Plane([-height, hwidth, 0], apex )
0216         rhs = Plane([ height, hwidth, 0], apex )
0217         bot = Plane([      0,     -1, 0], base )
0218         bck = Plane([      0,      0, -1], back )
0219         frt = Plane([      0,      0,  1], front )
0220
0221         self.lhs = lhs
0222         self.rhs = rhs
0223         self.bot = bot
0224         self.ymin = ymin
0225
0226         self.planes = [rhs, lhs, bot, frt, bck ]
0227         self.height = height
0228         self.depth = depth
0229         self.apex = np.array(apex)
0230
0231
0232     def __repr__(self):
0233         return "Prism(%s,%s) alpha %s " % (repr(self.parameters), repr(self.boundary), self.alpha)
0234
0235
0236     def intersectframe(self, ray):
0237         """
0238         Form coordinate system based centered on intersection point,
0239         with surface normal along Y.
0240         Another point in the plane of the intersected face,
0241         is used to pick the X direction.
0242         """
0243         isect = self.intersect(ray)
0244         assert len(isect)>0
0245         i0, i1 = isect
0246
0247         ifr = IntersectFrame( i0, a=prism.apex)
0248
0249         ## check transforming from world frame into intersect frame
0250         p_if = ifr.world_to_intersect(ifr.p)
0251         n_if = ifr.world_to_intersect(ifr.n, w=0)     # direction, not coordinate
0252         a_if = ifr.world_to_intersect(ifr.a)
0253
0254         pa_wf = ifr.a - ifr.p
0255         pa_if = ifr.world_to_intersect(pa_wf, w=0)
0256         pa_if /= np.linalg.norm(pa_if)
0257
0258         log.info("intersect position       ifr.p %s in p_if  %s " % (ifr.p, p_if ))
0259         log.info("intersect surface normal ifr.n %s in n_if  %s " % (ifr.n, n_if ))
0260         log.info("                         ifr.a %s in a_if  %s " % (ifr.a, a_if ))
0261         log.info("                         pa_wf %s in pa_if %s " % (pa_wf, pa_if ))
0262
0263         assert np.allclose( p_if, np.array([0,0,0,1]))
0264         assert np.allclose( n_if, np.array([0,1,0,0]))
0265         assert np.allclose( pa_if, np.array([1,0,0,0]))
0266
0267         ## check transforming from intersect frame into world frame
0268         o_if = np.array([0,0,0])
0269         o_wf = ifr.intersect_to_world(o_if)
0270         log.info("                        o_if %s o_wf %s " % (o_if, o_wf ))
0271         assert np.allclose( o_wf[:3], ifr.p )
0272
0273         return ifr
0274
0275
0276
0277     def intersect(self, ray):
0278          t0 = -np.inf
0279          t1 =  np.inf
0280          n0 = np.array([0,0,0])
0281          n1 = np.array([0,0,0])
0282
0283          for i, pl in enumerate(self.planes):
0284              n = pl.n
0285
0286              denom, t = pl.intersect(ray)
0287
0288              if denom < 0.:
0289                 if t > t0:
0290                     i0 = i
0291                     t0 = t
0292                     n0 = n
0293              else:
0294                  if t < t1:
0295                     i1 = i
0296                     t1 = t
0297                     n1 = n
0298
0299              log.info("i %2d denom %10.4f t %10.4f t0 %10.4f t1 %10.4f n0 %25s n1 %25s " % (i, denom, t, t0, t1, n0, n1 ))
0300
0301          if t0 > t1:
0302              log.into("no intersect")
0303              return []
0304          else:
0305              p0 = ray.position(t0)
0306              p1 = ray.position(t1)
0307              return Intersect(i0, t0, n0, p0),Intersect(i1, t1, n1, p1)
0308
0309     def intersect_unrolled(self, ray):
0310         """
0311         http://tavianator.com/fast-branchless-raybounding-box-intersections-part-2-nans/
0312
0313         # 0 * inf = nan   comparisons with nan always false
0314
0315         """
0316         pl = self.planes[0]
0317         t0 = pl
0318
0319
0320
0321 class PrismExpected(object):
0322     def __init__(self, a, n):
0323         """
0324         :param a: apex angle of prism
0325         :param n: refractive index of prism
0326         """
0327         self.a = a
0328         self.sa = np.sin(a)
0329         self.ca = np.cos(a)
0330
0331         self.single = type(n) is np.float64
0332
0333         self.n = n
0334         self.nn = n*n
0335         self.ni = 1./n
0336
0337         self.i2c = self.i2c_()
0338         self.i1c = self.i1c_()
0339
0340     def i2c_(self):
0341         """
0342         Critical angle
0343         """
0344         return np.arcsin(1./self.n)
0345
0346     def i1c_(self):
0347         """
0348         Other Critical angle
0349         """
0350         return np.arcsin( self.n*np.sin(self.a - np.arcsin(self.ni)))
0351
0352     def st2_(self, i1):
0353         return self.sa*np.sqrt(self.nn - np.sin(i1)*np.sin(i1)) - np.sin(i1)*self.ca
0354
0355     def t2_(self, i1):
0356         return np.arcsin(self.sa*np.sqrt(self.nn - np.sin(i1)*np.sin(i1)) - np.sin(i1)*self.ca)
0357
0358     def delta_(self, i1):
0359         return i1 + self.t2_(i1) - self.a
0360
0361     def i1mindev_(self):
0362         return np.arcsin( self.n*np.sin(self.a/2.) )  # incident angle with minimum deviation
0363
0364     def i1_domain(self):
0365         assert self.single
0366         _i1c = self.i1c/rad
0367         _i2c = self.i2c/rad
0368         log.debug( "_i1c %s " % _i1c)
0369         log.debug( "_i2c %s " % _i2c)
0370         log.info("plt.plot(dom,dl)")
0371         return np.linspace(_i1c+1e-9,90,200)
0372
0373     def expected(self):
0374         assert self.single
0375         dom = self.i1_domain()
0376         dl = self.delta_(dom*rad)/rad
0377         return dom, dl
0378
0379     def spawn_singles(self):
0380         """
0381         For plotting need a domains and corresponding values, BUT must
0382         split up by refractive index to avoid spagetti
0383         """
0384         un = np.unique(self.n)
0385         if len(un) > 10:
0386            log.warn("too many distinct indices");
0387            return []
0388
0389         return [PrismExpected(self.a, n) for n in un]
0390
0391
0392
0393 class PrismCheck(object):
0394     """
0395     Assumes canonical passage thru prism (ie not out the bottom)
0396     """
0397     def title(self):
0398         return "prism.py deviation vs incident angle"
0399
0400     def __init__(self, prism, xprism, sel, mask=True):
0401
0402         self.prism = prism
0403         self.xprism = xprism
0404         self.sel = sel
0405
0406         p0 = sel.rpost_(0)[:,:3]  # light source position
0407         p1 = sel.rpost_(1)[:,:3]  # 1st refraction point
0408         p2 = sel.rpost_(2)[:,:3]  # 2nd refraction point
0409         p3 = sel.rpost_(3)[:,:3]  # light absorption point
0410
0411         assert len(p0) == len(p1) == len(p2) == len(p3)
0412         N = len(p0)
0413
0414         w0 = sel.recwavelength(0)
0415         w1 = sel.recwavelength(1)
0416         w2 = sel.recwavelength(2)
0417         w3 = sel.recwavelength(3)
0418
0419         assert len(w0) == len(w1) == len(w2) == len(w3)
0420         assert len(w0) == len(p0)
0421         assert np.all(w0 == w1)
0422         assert np.all(w0 == w2)
0423         assert np.all(w0 == w3)
0424         self.wx = w0
0425         # prism does spatial sorting, no wavelength changes as no reemission
0426
0427         self.p0 = p0
0428         self.p1 = p1
0429         self.p2 = p2
0430         self.p3 = p3
0431
0432
0433         p01 = p1 - p0
0434         p12 = p2 - p1
0435         p23 = p3 - p2
0436
0437         self.p01 = p01
0438         self.p12 = p12
0439         self.p23 = p23
0440
0441         cdv = costheta_( p01, p23 )    # total deviation angle
0442         dv = np.arccos(cdv)
0443
0444         self.cdv = cdv
0445         self.dv = dv    # simulated deviation
0446
0447         lno = prism.lhs.ntile(N)  # prism lhs normal, repeated
0448         rno = prism.rhs.ntile(N)  # prism rhs normal, repeated
0449
0450         ci1 = costheta_(-p01, lno )    # incident 1
0451         ct1 = costheta_(-p12, lno )    # transmit 1
0452         ci2 = costheta_( p12, rno )    # incident 2
0453         ct2 = costheta_( p23, rno )    # transmit 2
0454
0455         i1 = np.arccos(ci1)
0456         t1 = np.arccos(ct1)
0457         i2 = np.arccos(ci2)
0458         t2 = np.arccos(ct2)
0459
0460         self.i1 = i1
0461         self.t1 = t1
0462         self.i2 = i2
0463         self.t2 = t2
0464
0465         # Snell check : refractive index from the angles at the 2 refractions
0466         n1 = np.sin(i1)/np.sin(t1)
0467         n2 = np.sin(t2)/np.sin(i2)
0468         dn = n1 - n2
0469
0470         #assert dn.max() < 1e-3
0471         #assert dn.min() > -1e-3
0472
0473
0474         self.n1 = n1
0475         self.n2 = n2
0476
0477         self.expected_deviation()
0478         self.compare_expected_with_simulated()
0479
0480     def expected_deviation(self):
0481         """
0482         Get 8 nan close to critical angle
0483
0484         In [108]: pc.i1[pc.inan]*deg
0485         Out[108]:
0486         array([ 24.752,  24.752,  24.752,  24.752,  24.752,  24.752,  24.752,
0487                 24.752])
0488
0489         simulating just those around critical angle,
0490         see that they are just sneaking out the prism grazing the face
0491         """
0492         i1 = self.i1
0493
0494         xdv = self.xprism.delta_(i1)
0495         ina = np.arange(len(xdv))[np.isnan(xdv)]
0496         msk = np.arange(len(xdv))[~np.isnan(xdv)]
0497
0498         if len(ina)>0:
0499            log.warning("expected_deviation needs nan masking ina:%s len(ina):%s " % (ina, len(ina)))
0500
0501         self.xdv = xdv
0502         self.ina = ina
0503         self.msk = msk
0504
0505         i1c = self.xprism.i1c_()
0506
0507         log.info("ina: %s " % self.ina)
0508         log.info("i1[ina]/rad: %s " % (i1[self.ina]/rad) )
0509         log.info("i1c/rad : %s " % (i1c/rad) )
0510
0511     def compare_expected_with_simulated(self):
0512         """
0513         """
0514         msk = self.msk
0515         mdf = self.dv[msk] - self.xdv[msk]
0516         log.info("  dv[msk]/rad %s " % (self.dv[msk]/rad) )
0517         log.info(" xdv[msk]/rad %s " % (self.xdv[msk]/rad) )
0518         log.info("mdf/rad:%s max:%s min:%s " % (mdf/rad, mdf.min()/rad, mdf.max()/rad ))
0519
0520         self.mdf = mdf
0521         self.mdv = self.dv[msk]
0522         self.mi1 = self.i1[msk]
0523
0524
0525
0526
0527 def test_intersectframe(prism):
0528     """
0529     Establish a frame at midface lhs intersection point with the prism
0530     """
0531     ray = Ray([-600,0,0], [1,0,0])
0532     ifr = prism.intersectframe(ray)
0533     i1m = prism.i1mindev()
0534     ti1m = np.tan(i1m)
0535     pm_if =  np.array([-1, 1./ti1m,0])*400
0536     pm_wf =  ifr.intersect_to_world(pm_if)
0537     log.info(" mindev position pm_if %s pm_wf %s " % (pm_if, pm_wf ))
0538     s_i2w = ifr.i2w_string()
0539     log.info("  s_i2w %s " % s_i2w );
0540
0541
0542
0543 def scatter_plot(xq, yq, sl):
0544     if sl is not None:
0545         x = xq[sl]
0546         y = yq[sl]
0547     else:
0548         x = xq
0549         y = yq
0550     pass
0551     plt.scatter(x*deg, y*deg)
0552
0553
0554 def vanity_plot(pc, sl=None):
0555     for xprism in pc.xprism.spawn_singles():
0556         dom, dl = xprism.expected()
0557         plt.plot(dom,dl)
0558
0559     scatter_plot(pc.i1, pc.dv, sl)
0560
0561 def deviation_plot(pc,sl=None):
0562     """
0563     masked plotting to avoid critical angle nan
0564
0565     large (+-0.5 degree) deviations between expected and simulated delta
0566     all occuring close to critical angle
0567
0568     away from critical angle, deviations less that 0.1 degree
0569     """
0570     scatter_plot(pc.mi1, pc.mdf, sl)
0571
0572
0573 def oneplot(pc, log_=False):
0574     fig = plt.figure()
0575     plt.title(pc.title())
0576     ax = fig.add_subplot(111)
0577
0578     vanity_plot(pc, sl=slice(0,10000))
0579     #deviation_plot(pc, sl=slice(0,10000))
0580     #deviation_plot(pc, sl=None)
0581
0582     fig.show()
0583
0584
0585
0586 def spatial(pc):
0587     """
0588     """
0589     w = pc.wx
0590     x = pc.p3[:,0]
0591     y = pc.p3[:,1]
0592     z = pc.p3[:,2]
0593
0594     #assert np.all(x == 1200.)
0595     #off = x != 1200.
0596     #print pc.p3[off]
0597
0598     from matplotlib.colors import LogNorm
0599
0600     ax.hist2d( w, y, bins=100, norm=LogNorm())
0601
0602
0603
0604
0605
0606
0607 if __name__ == '__main__':
0608
0609     args = opticks_main(tag="1", det="prism")
0610
0611
0612     plt.ion()
0613
0614     #wl = np.arange(10, dtype=np.float32)*70. + 100.
0615     # low range is off edge of the refractive index values
0616
0617     seqs = ["TO BT BT SA"]
0618
0619     try:
0620         sel = Evt(tag=args.tag, det=args.det, seqs=seqs, args=args)
0621     except IOError as err:
0622         log.fatal(err)
0623         sys.exit(args.mrc)
0624
0625     log.info("sel %s " % sel.brief)
0626
0627
0628     if not sel.valid:
0629         log.fatal("failed to load tag %s det %s " % (args.tag, args.det))
0630         sys.exit(1)
0631
0632     boundary = Boundary("Vacuum///GlassSchottF2")
0633
0634     prism = Prism("60.,300,300,0", boundary)
0635
0636     log.info("prism %s " % repr(prism))
0637
0638
0639     n = boundary.imat.refractive_index(sel.wl)
0640
0641     xprism = PrismExpected(prism.a, n)
0642
0643     pc = PrismCheck(prism, xprism, sel )
0644
0645
0646     oneplot(pc, log_=False)
0647     #spatial(pc)
0648
0649
0650     plt.show()
0651
0652
0653