eic-opticks/ana/droplet.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 droplet.py : Analytic All Order Rainbow angle calculations
0023 =============================================================
0024
0025 Geometrical calculation of deviation, incident and refracted angles
0026 at minimum deviation for k orders of rainbows.
0027
0028 """
0029 import os, logging, numpy as np
0030 log = logging.getLogger(__name__)
0031
0032 import matplotlib.pyplot as plt
0033 from matplotlib.patches import Rectangle
0034
0035 from opticks.ana.base import opticks_environment
0036 from opticks.ana.boundary import Boundary
0037 deg = np.pi/180.
0038
0039
0040 class Droplet(object):
0041     def __init__(self, boundary):
0042         self.boundary = boundary
0043
0044     @classmethod
0045     def seqhis(cls, arg, src=None):
0046         pp = [arg] if type(arg) is int else arg
0047         return map(lambda _:cls.seqhis_(_,src=src), pp)
0048
0049     @classmethod
0050     def seqhis_(cls, p, src=None):
0051         seq = "" if src is None else src + " "
0052         if p == 0:
0053             seq += "BR "
0054         elif p == 1:
0055             seq += "BT BT "
0056         elif p > 1:
0057             seq += "BT " + "BR " * (p-1) + "BT "
0058         else:
0059             assert 0
0060         pass
0061         seq += "SA"
0062         return seq
0063
0064
0065     def deviation_angle(self, w, k=1):
0066         d = self.deviation_angle_(w, k=k)
0067         return d['dv']
0068
0069     def deviation_angle_(self, w, k=1):
0070         """
0071         tot deviation, incident, refracted angles at the minimum deviation
0072         """
0073         if w is None:
0074             w = np.array([780., 380.])  # red, blue
0075
0076         n = self.boundary.imat.refractive_index(w)
0077
0078         i = np.arccos( np.sqrt((n*n - 1.)/(k*(k+2.)) ))
0079         r = np.arcsin( np.sin(i)/n )
0080         dv = ( k*np.pi + 2*i - 2*r*(k+1) ) % (2*np.pi)
0081
0082         return dict(n=n,i=i,r=r,dv=dv,k=k,w=w)
0083
0084     def rainbow_table(self):
0085         redblue = np.array([780., 380.])
0086         lfmt = "%3s " + " %10s " * 7
0087         print lfmt % ( "k", "th(red)", "th(blue)", "dth", "i(red)", "i(blue)", "r(red)", "r(blue)" )
0088
0089         for k in range(1,21):
0090             d = self.deviation_angle_(redblue, k=k)
0091             dv = d['dv']
0092             i = d['i']
0093             r = d['r']
0094             fmt = "%3d " + " %10.2f " * 7
0095             print fmt % (k, dv[0]/deg, dv[1]/deg, (dv[1]-dv[0])/deg, i[0]/deg, i[1]/deg, r[0]/deg, r[1]/deg )
0096
0097     def bow_angle_rectangles(self, ks=range(1,11+1), yleg=2):
0098         ax = plt.gca()
0099         ymin, ymax = ax.get_ylim()
0100         for k in ks:
0101             dvr = self.deviation_angle(None, k)/deg
0102             rect = Rectangle( (dvr[0], ymin), dvr[1]-dvr[0], ymax-ymin, alpha=0.1 )
0103             ax.add_patch(rect)
0104             ax.annotate( "%s" % k, xy=((dvr[0]+dvr[1])/2, yleg), color='red')
0105         pass
0106
0107
0108
0109
0110 if __name__ == '__main__':
0111     logging.basicConfig(level=logging.INFO)
0112     opticks_environment()
0113
0114     boundary = Boundary("Vacuum///MainH2OHale")
0115     droplet = Droplet(boundary)
0116
0117     droplet.rainbow_table()
0118
0119     plt.ion()
0120     fig = plt.figure()
0121     ax = plt.gca()
0122     ax.set_ylim([0,4])
0123     ax.set_xlim([0,360])
0124
0125     droplet.bow_angle_rectangles()
0126
0127
0128