eic-opticks/ana/fresnel.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 fresnel.py : analytic reflection expectations
0023 ==================================================
0024
0025
0026 """
0027 import os, logging
0028 log = logging.getLogger(__name__)
0029 import numpy as np
0030 import matplotlib.pyplot as plt
0031 from mpl_toolkits.mplot3d import Axes3D
0032
0033 from opticks.ana.base import opticks_environment
0034 from opticks.ana.nbase import count_unique
0035
0036 np.set_printoptions(suppress=True, precision=3)
0037
0038
0039
0040 def fresnel(x, n1, n2, spol=True):
0041     """
0042     https://en.wikipedia.org/wiki/Fresnel_equations
0043     """
0044     cx = np.cos(x)
0045     sx = np.sin(x)
0046     disc = 1. - np.square(n1*sx/n2)
0047     qdisc = np.sqrt(disc)
0048     pass
0049     if spol:
0050         num = (n1*cx - n2*qdisc)
0051         den = (n1*cx + n2*qdisc)
0052     else:
0053         num = (n1*qdisc - n2*cx)
0054         den = (n1*qdisc + n2*cx)
0055     pass
0056     return np.square(num/den)
0057
0058
0059 def fresnel_factor(seqhis, i, n1, n2, spol=True):
0060     """
0061     :param seqhis: history sequence string eg "TO BT BR BT SA"
0062     :param n1: refractive index of initial material
0063     :param n2: refractive index of material that is transmitted into
0064
0065     Not aiming for generality, only works for simple geometries like raindrops, prisms, lens
0066     """
0067
0068     seqs = seqhis.split(" ")
0069
0070     rx = fresnel(i, n1, n2, spol=spol )
0071     tx = 1 - rx
0072
0073     ff = np.ones(len(i))
0074     for step in seqs:
0075         #print step
0076         if step in ("TO", "SA"):continue
0077         if step == "BT":
0078             ff *= tx
0079         elif step == "BR":
0080             ff *= rx
0081         else:
0082             assert 0, step
0083         pass
0084     pass
0085     return ff
0086
0087
0088
0089
0090 def fresnel_s( i, n, method=0):
0091     """
0092     sin(i-r)    si cr - ci sr
0093     -------- =  -------------
0094     sin(i+r)    si cr + ci sr
0095
0096     This form whilst pretty, gives nan at normal incidence, 0/0
0097     """
0098     si = np.sin(i)
0099     sr = si/n
0100
0101     if method == 0:
0102         ci = np.sqrt( 1 - si*si )
0103         cr = np.sqrt( 1 - sr*sr )
0104         num = si*cr - ci*sr
0105         den = si*cr + ci*sr
0106     else:
0107         i = np.arcsin(si)
0108         r = np.arcsin(sr)
0109         num = np.sin(i - r)
0110         den = np.sin(i + r)
0111         #log.info("i %s r %s num %s den %s " % (i,r,num,den))
0112     pass
0113     return np.square(num/den)
0114
0115
0116 def fresnel_p( i, n):
0117     """
0118     tan(i-r)
0119     --------
0120     tan(i+r)
0121     """
0122     si = np.sin(i)
0123     sr = si/n
0124     i = np.arcsin(si)
0125     r = np.arcsin(sr)
0126     num = np.tan(i - r)
0127     den = np.tan(i + r)
0128     return np.square(num/den)
0129
0130
0131
0132 class Fresnel(object):
0133     def __init__(self, n1, n2, dom=None ):
0134         if dom is None:
0135             dom = np.linspace(0,90,91)
0136
0137         n1 = np.asarray(n1)
0138         n2 = np.asarray(n2)
0139
0140         th = dom*np.pi/180.
0141         spol = fresnel(th, n1, n2, True)
0142         ppol = fresnel(th, n1, n2, False)
0143
0144         pass
0145         self.n1 = n1
0146         self.n2 = n2
0147         self.dom = dom
0148         self.th = th
0149
0150         self.spol_0 = spol
0151         self.ppol_0 = ppol
0152         #self.alternative_check()
0153
0154         self.cen  = (dom[:-1] + dom[1:])/2.
0155         # avg of bin edge values
0156         self.spol = (spol[:-1] + spol[1:])/2.
0157         self.ppol = (ppol[:-1] + ppol[1:])/2.
0158         self.upol = (self.spol+self.ppol)/2.     # unpol?
0159
0160         self.brewster = np.arctan(n2/n1)*180./np.pi
0161         self.critical = np.arcsin(n1/n2)*180./np.pi
0162
0163
0164     def alternative_check(self):
0165         """
0166         Alternative angle difference forms, misbehave at normal incidence
0167         Otherwise they match
0168         """
0169         th = self.th
0170         n1 = self.n1
0171         n2 = self.n2
0172
0173         spol_0 = self.spol_0
0174         ppol_0 = self.ppol_0
0175
0176         spol_2 = fresnel_s( th, n2/n1, method=1)
0177         spol_3 = fresnel_s( th, n2/n1, method=0)
0178         assert np.allclose( spol_0[1:], spol_2[1:] ), np.dstack([spol_0,spol_2, spol_3])
0179         assert np.allclose( spol_0[1:], spol_3[1:] ), np.dstack([spol_0,spol_2, spol_3])
0180
0181         ppol_2 = fresnel_p( th, n2/n1)
0182         assert np.allclose( ppol_0[1:], ppol_2[1:] ), np.dstack([ppol_0, ppol_2])
0183
0184
0185
0186     def __call__(self, xd, n):
0187         x = xd*np.pi/180.
0188         n1 = self.n1
0189         n2 = self.n2
0190         cx = np.cos(x)
0191         sx = np.sin(x)
0192         disc = 1. - np.square(n1*sx/n2)
0193         qdisc = np.sqrt(disc)
0194         pass
0195         spol = np.square((n1*cx - n2*qdisc)/(n1*cx + n2*qdisc))
0196         ppol = np.square((n1*qdisc - n2*cx)/(n1*qdisc + n2*cx))
0197         return n*(spol*f + (1.-f)*ppol)
0198
0199
0200     def pl(self):
0201         plt.plot(self.cen, self.spol, label="S (perp)", c="r")
0202         plt.plot(self.cen, self.ppol, label="P (para)", c="b")
0203
0204     def title(self):
0205         return "Fresnel %4.3f/%4.3f " % (self.n1, self.n2  )
0206
0207     def plot(self, fig, ny=1, nx=1, n=1, log_=False):
0208         plt.title(self.title())
0209         ax = fig.add_subplot(ny,nx,n)
0210         self.pl()
0211         self.angles()
0212         legend = ax.legend(loc='upper left', shadow=True)
0213         if log_:
0214             ax.set_yscale('log')
0215
0216
0217     def angles(self):
0218         a = self.brewster
0219         plt.plot([a, a], [1e-6, 1], 'k-', c="b", lw=2, label="Brewster")
0220         a = self.critical
0221         plt.plot([a, a], [1e-6, 1], 'k-', c="r", lw=2, label="Critical")
0222
0223
0224 if __name__ == '__main__':
0225     logging.basicConfig(level=logging.INFO)
0226     opticks_environment()
0227
0228     n1 = np.array([1.])
0229     n2 = np.array([1.458])
0230
0231     fr = Fresnel(n1,n2)
0232     fig = plt.figure()
0233     fr.plot(fig, log_=True)
0234     fig.show()
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244