Added python version of mode-matching code by Bain Bronner

13 files changed, 824 insertions, 0 deletions

diff --git a/beam_tracing/python/README.txt b/beam_tracing/python/README.txt
new file mode 100644
index 0000000..75faf8a
--- /dev/null
+++ b/beam_tracing/python/README.txt
@@ -0,0 +1,20 @@
+Original author: Eugeniy Mikhailov

+Conversion from Octave to Python: Bain Bronner

+

+Language: Python 2.7xx

+Requires: numpy 1.6.x and scipy 0.10.x

+

+prismDiskCoupling.py calculates and draws the critical angle for a given prism geometry and indexes of refraction for the prism and disk

+demoPrismDiskCoupling.py is an example call of the above function

+

+beamTrace.py traces the beam according to Fresnel equations with lighter colors indicating less intensity

+exampleBeamTrace.py is an example call of the above function

+

+the rest of the functions serve as helpers to the above

+

+changes:

+	1) in converting to Python, the beams initial travel angle and start position for the exampleBeamTrace were determined by a slightly modified prismDiskCoupling.py

+	2) slightly different annotation styles, including the use of arrows, for the angles in prismDiskCoupling

+	3) variable and files converted to python naming standard, i.e. this_is_a_file.m --> thisIsAFile.py

+	4) more or less same logic and mathematics as previously
+	5) disk placement was added, so now the disk is not always assumed to be in the center

diff --git a/beam_tracing/python/TODO.txt b/beam_tracing/python/TODO.txt
new file mode 100644
index 0000000..e3ef167
--- /dev/null
+++ b/beam_tracing/python/TODO.txt
@@ -0,0 +1,7 @@
+TODO:

+

+1/23/13

+Add in the ability to easily change the disk location and subsequent beam tracking, i.e. not have the disk always at the center of the prism

+**2/3/13 Above done.
+

+GUI for easy input of different indexes of refraction, prism angle, etc.?

diff --git a/beam_tracing/python/beam2FaceDistance.py b/beam_tracing/python/beam2FaceDistance.py
new file mode 100644
index 0000000..8aeedec
--- /dev/null
+++ b/beam_tracing/python/beam2FaceDistance.py
@@ -0,0 +1,53 @@
+import numpy as np

+

+def beam2FaceDistance(beam, face):

+# return  distance to face if beam hits it or infinity otherwise

+# for beams and faces structure description see face_beam_interaction.m

+

+   # default return values for the case when beam misses the face

+   isFaceHit = False

+   hitPosition = [float('NaN'), float('NaN')]

+   hitDistance = float('Inf')

+

+   k = beam.k

+   # face direction or kf

+   kf = face.vertex2 - face.vertex1 # not a unit vector

+

+   # simple check for intersection of two vectors 

+   # if beam are parallel no intersection is possible

+   # we do this via cross product calculation

+   

+   angleTolerance = 1e-14;

+   if abs(np.cross(k, kf)) <= angleTolerance:

+      # beams misses face

+      return [hitDistance, hitPosition, isFaceHit]

+   

+   # let's find the intersection of lines passing through face and light beam

+   # we introduce parameters tb and tf 

+   # which define distance (in arb. units) along k vectors

+   # we are solving 

+   # xb_o+ k_x*tb = v1_x+ kf_x*tf

+   # yb_o+ k_y*tb = v1_y+ kf_y*tf

+   # A*t = B

+   A = np.array([[k[0], -kf[0]],\

+                [k[1], -kf[1]]])

+   B = np.array([[face.vertex1[0] - beam.origin[0]],\

+                [face.vertex1[1] - beam.origin[1]]])

+   A = np.linalg.inv(A)

+   t = np.dot(A, B)

+   tb = t[0] 

+   tf = t[1]

+

+   # beam intersects face only if 0<=tf<=1 and tb>=0

+   if 0 <= tf and tf <= 1 and tb >= angleTolerance:

+      # intersection,  beam hits the face

+      isFaceHit = True

+   else:

+      # beam misses the face

+      return [hitDistance, hitPosition, isFaceHit]

+

+   # calculating hit position

+   hitPosition = face.vertex1 + kf * tf

+   hitDistance = tb # distance along the light beam

+

+   return [hitDistance, hitPosition, isFaceHit]
\ No newline at end of file
diff --git a/beam_tracing/python/beamTrace.py b/beam_tracing/python/beamTrace.py
new file mode 100644
index 0000000..6f2cd96
--- /dev/null
+++ b/beam_tracing/python/beamTrace.py
@@ -0,0 +1,49 @@
+import pylab

+from faceBeamInteraction import *

+import numpy as np

+

+def beamTrace(beams, faces, borders, borderLimits):

+   # trace beam and all it reflections and refractions on faces

+   # for beams and faces structure description see faceBeamInteraction.m

+   # border similar to faces cell array which enclose the beam

+   # i.e. it does not leave the polygon consisting of border faces

+

+   intensityThreshold = 1e-2 # intensity of the weakest beam we still trace

+   processedBeams = []

+

+   nBeams = len(beams)

+   if nBeams == 0:

+      # no more beam to trace

+      return processedBeams

+

+   # trace each beam in beams cell array

+   intenseEnoughBeams = []

+   intenseBeamsCounter = 0

+   for beam in beams:

+      [isFaceHit, hitPosition, hitDistance, newBeams] =\

+       faceBeamInteraction(beam, faces)

+      if (isFaceHit):

+         beam.hitPosition = hitPosition

+         # remove beams which are to weak to trace

+	 nNewBeams = len(newBeams)

+	 for newBeam in newBeams:

+	    if newBeam.intensity >= intensityThreshold:

+	       intenseBeamsCounter += 1

+	       intenseEnoughBeams.append(newBeam)

+	 #try:

+	 newProcessedBeams = \

+	     beamTrace(intenseEnoughBeams, faces, borders, borderLimits)

+	 #except ValueError:

+	  #  break

+	 processedBeams.append(newProcessedBeams)

+      else:

+	 # beam does not hit face but it should stop at borders

+	 beam.face = None

+	 [isFaceHit, hitPosition, hitDistance, newBeams] = \

+	  faceBeamInteraction(beam, borders)

+	 if not isFaceHit:

+	    print 'borders are badly defined, the beam misses them'

+	    return False

+	 beam.hitPosition = hitPosition

+   processedBeams.append(newBeams)

+   return processedBeams

diff --git a/beam_tracing/python/classes.py b/beam_tracing/python/classes.py
new file mode 100644
index 0000000..d4ee631
--- /dev/null
+++ b/beam_tracing/python/classes.py
@@ -0,0 +1,37 @@
+import numpy as np

+

+class beam:

+   def __init__(self):

+      '''light beam class'''

+      # source of light beam; default (0,0)

+      beam.origin = np.array([])

+      # k-vector or direction of light beam; default (1,1)

+      beam.k = np.array([])

+      # intensity of light beam; default max

+      beam.intensity = 1.0

+      # if beam starts at a face, the face index; default none

+      beam.face = None

+      # 0 for s-polarization 1 for p-polarization; default s

+      beam.polarization = 0

+      # status: incoming, reflected, or refracted; default incoming

+      beam.status = 'incoming'

+      

+class face:

+   def __init__(self):

+      '''prism face class'''

+      # first vertex point

+      face.vertex1 = np.array([])

+      # second vertex point

+      face.vertex2 = np.array([])

+      # left (right) side is left (right) of line from vertex 1 to 2

+      # index of refraction for the left side

+      face.nLeft = np.array([])

+      # index of refraction for the right side

+      face.nRight = np.array([])

+      

+class border:

+   def __init__(self):

+      border.vertex1 = np.array([])

+      border.vertex2 = np.array([])

+      border.nLeft = np.array([])

+      border.nRight = np.array([])
\ No newline at end of file
diff --git a/beam_tracing/python/colorGradient.py b/beam_tracing/python/colorGradient.py
new file mode 100644
index 0000000..dd9c856
--- /dev/null
+++ b/beam_tracing/python/colorGradient.py
@@ -0,0 +1,32 @@
+import colorsys

+import numpy as np

+

+def colorGradient(pos, baseColorRGB, nColors):

+# returns rgb=[R,G,B] components from the brightness gradient of the base color

+# base_color_rgb - [R,G,B] of the base color

+# Ncolors - number of positions in gradient

+

+   # color map  as gradient of particular color

+   # starting values

+   hsv_s = colorsys.rgb_to_hsv(baseColorRGB[0], baseColorRGB[1], baseColorRGB[2])

+   hue_f = hsv_s[0]

+   sat_f = hsv_s[1]

+   val_f = hsv_s[2]

+   # stop values

+   hue_s = hue_f

+   sat_s = 0

+   val_s = 1

+   hue = np.linspace(hue_s, hue_f, nColors)

+   sat = np.linspace(sat_s, sat_f, nColors)

+   val = np.linspace(val_s, val_f, nColors)

+ 

+   # check that pos is with in the reach

+   pos = np.floor(pos)

+   if pos < 1:

+      pos = 1

+   if pos > nColors:

+      pos = nColors

+   

+   pos = int(pos)

+   RGB = colorsys.hsv_to_rgb(hue[pos-1], sat[pos-1], val[pos-1])

+   return RGB
\ No newline at end of file
diff --git a/beam_tracing/python/demoPrismDiskCoupling.py b/beam_tracing/python/demoPrismDiskCoupling.py
new file mode 100644
index 0000000..60ef351
--- /dev/null
+++ b/beam_tracing/python/demoPrismDiskCoupling.py
@@ -0,0 +1,21 @@
+from prismDiskCoupling import *

+

+# prismDiskCoupling(prismAngle, nDisk, nPrism, diskX, coupDesc)

+   # prismAngle - angle of bottom right prism corner

+   # nDisk - disk index of refraction

+   # nPrism - prism index of refraction

+   # diskX - x coordinate of disk contact, can be between (-1, 1)

+   # coupDesc - coupling description of prism/disk info

+

+prismAngle = 60

+#prismAngle = 45

+

+#nDisk = 1.430 # CaF2

+nDisk = 2.256 # LiNbO3

+#nDisk = 2.176

+

+#nPrism = 2.52 # Rutile (TiO2)

+nPrism = 2.400 # Diamond (C)

+

+prismDiskCoupling(prismAngle, nDisk, nPrism, 0.0, 'Blank')

+pylab.show() 
\ No newline at end of file
diff --git a/beam_tracing/python/drawPrismDisk.py b/beam_tracing/python/drawPrismDisk.py
new file mode 100644
index 0000000..a099096
--- /dev/null
+++ b/beam_tracing/python/drawPrismDisk.py
@@ -0,0 +1,102 @@
+import numpy as np

+import pylab

+import sys

+

+

+def drawPrismDisk(prismAngleInDegrees, nDisk, nPrism, diskX, coupDesc):

+   # Calculates incident angle for proper coupling into the disc via prism

+   # prismAngleInDegrees - angle of the prism faces in degrees

+   # diskX - x coordinate of the disk center

+   # coupDesc - short annotation of the situation

+

+   prismAngle = prismAngleInDegrees * np.pi / 180 # convert to radians

+

+   # critical angle for beam from prism to disk. we want thetaDisk to be 90

+   # degrees, i.e. total internal reflection -> sin(thetaDisk) = 1.0

+

+   thetaPrism1 = np.arcsin(nDisk / nPrism)

+   thetaPrismInDegrees1 = thetaPrism1 * 180 / np.pi

+

+   # find inside angle on the incident side of prism

+   thetaPrism2 = prismAngle - thetaPrism1

+   thetaPrismInDegrees2 = thetaPrism2 * 180 / np.pi #convert to degrees

+

+   # calculate refracted angle out of the prism with respect to the normal

+   # positive means above the normal; negative below

+   tempArcsin = nPrism * np.sin(thetaPrism2)

+

+   try:

+      thetaAir = np.arcsin(tempArcsin)

+   except ValueError:

+      print 'Error'

+      sys.exit('Total internal reflection at right prism face')

+   

+   thetaAir = np.arcsin(tempArcsin)

+   thetaAirInDegrees = thetaAir * 180 / np.pi #convert to degrees

+

+   # angle in the air relative to horizon

+   thetaAirHorizon = thetaAir + (np.pi/2 - prismAngle)

+   thetaAirHorizonInDegrees = thetaAirHorizon * 180 / np.pi #convert to degrees

+

+#   -----GRAPHICAL OUTPUT-----

+# PRISM

+   # bottom face will go from (-1,0) to (1,0) regardless of prismAngle

+   xFace1 = np.linspace(-1,1)

+   yFace1 = 0 * xFace1

+   

+   # second face to the right of origin 

+   # at prismAngle with respect to negative x-direction

+   xFace2 = np.linspace(1,0)

+   yFace2 = (1 - xFace2) * np.tan(prismAngle)

+   

+   # third face to the left of origin

+   # at prismAngle with respect to positive x-direction

+   xFace3 = np.linspace(-1,0)

+   yFace3 = (xFace3 + 1) * np.tan(prismAngle)

+

+   # draw prism

+   pylab.plot(xFace1, yFace1, 'black',\

+              xFace2, yFace2, 'black',\

+              xFace3, yFace3, 'black')

+   

+

+# BEAMS

+   # Inside Prism

+   # drawing line from disk contact point to face2 at angle thetaPrism

+   xCross = (np.tan(prismAngle) + diskX / np.tan(thetaPrism1)) /\

+             (1/np.tan(thetaPrism1) + np.tan(prismAngle))

+   yCross = -np.tan(prismAngle) * (xCross - 1)

+   xBeamPrism = np.linspace(diskX, xCross)

+   yBeamPrism = np.linspace(0, yCross)

+   pylab.plot(xBeamPrism, yBeamPrism, color = 'red', linewidth = 2)

+   

+   # Outside Prism

+   xOutStart = xCross

+   yOutStart = yCross

+   xOutStop = 1.5

+   xBeamOut = np.linspace(xOutStart, xOutStop)

+   yBeamOut = yOutStart + np.tan(thetaAirHorizon) * (xBeamOut - xOutStart)

+   # this keeps the beam from going farther than necessary

+   i = 0

+   for value in yBeamOut:

+      if value > 2.0:

+         break

+      else:

+         i +=1

+   # ensure same dimension

+   yBeamOut = yBeamOut[:i]

+   xBeamOut = xBeamOut[:i]

+   pylab.plot(xBeamOut, yBeamOut, color = 'red', linewidth = 2)

+   

+   

+# General annotation and axis set-up

+   pylab.text(-1.4, 1.7, coupDesc)

+   pylab.text(-1.4, 1.55, 'nDisk = %.3f' % nDisk)

+   pylab.text(-1.4, 1.4, 'nPrism = %.3f' % nPrism)   

+   pylab.xlim(-1.5, 1.5)

+   pylab.ylim(-0.5, 2.0)

+   pylab.title(coupDesc)

+

+

+   return [xBeamOut[-1], yBeamOut[-1], yFace2[-1], yFace3[-1],\

+           thetaPrism1, thetaPrism2, thetaAirHorizon]
\ No newline at end of file
diff --git a/beam_tracing/python/exampleBeamTrace.py b/beam_tracing/python/exampleBeamTrace.py
new file mode 100644
index 0000000..1e35118
--- /dev/null
+++ b/beam_tracing/python/exampleBeamTrace.py
@@ -0,0 +1,122 @@
+from beamTrace import *

+from classes import *

+from colorGradient import *

+from drawPrismDisk import *

+from faceBeamInteraction import *

+from fresnelReflection import *

+import pylab

+import numpy as np

+

+#prismAngle = 45 # rutile prism

+prismAngle = 60 # diamond prism

+

+nDisk = np.array([2.256, 2.176]) # LiNbO3

+#nDisk = np.array([1.375, 1.387]) #MgF2

+#nDisk = np.array([1.430, 1.430]) #CaF2

+

+nPrism = np.array([2.400, 2.400]) # diamond

+#nPrism = np.array([2.521, 2.793]) # rutile

+

+nAir = np.array([1.0, 1.0])

+

+coupDesc = 'LiNbO3 disk\nDiamond (C) Prism'

+

+# keep in mind the bottom of prism always goes from -1 to 1, regardless of

+# prismAngle

+diskX = -0.85

+

+

+# the following calls prismDiskCoupling from drawPrismDisk

+# it is the same as prismDiskCoupling function, with unnecessary parts removed

+# calling this function allows for different ideal beam start locations to be

+# easily put into beamTrace and its helpers

+

+[xBeamOrigin, yBeamOrigin, yFace2, yFace3,\

+ thetaPrism, thetaPrism2, thetaAirRth] =\

+ drawPrismDisk(prismAngle, nDisk[0], nPrism[0], diskX, coupDesc)

+

+

+

+# ------------------------------------------------------------

+# Should not have to alter below to acquire different outputs

+# ------------------------------------------------------------

+

+

+# prism faces

+face1 = face()

+face1.vertex1 = np.array([-1, 0])

+face1.vertex2 = np.array([1, 0])

+# left (right) side is left (right) of line from vertex 1 to 2

+face1.nLeft = nPrism

+face1.nRight = nDisk

+

+face2 = face()

+face2.vertex1 = np.array([1, 0])

+face2.vertex2 = np.array([0, yFace2])

+face2.nLeft = nPrism

+face2.nRight = nAir

+

+face3 = face()

+face3.vertex1 = np.array([-1, 0])

+face3.vertex2 = np.array([0, yFace3])

+face3.nLeft = nAir

+face3.nRight = nPrism

+

+faces = np.array([face1, face2, face3])

+

+#beams

+beamInitialTravelAngle = thetaAirRth

+

+beam1 = beam()

+beam1.k = np.array([-np.cos(beamInitialTravelAngle), \

+                    -np.sin(beamInitialTravelAngle)])

+beam1.origin = np.array([xBeamOrigin, yBeamOrigin])

+beam1.face = float('NaN')

+beam1.intensity = 1

+beam1.polarization = 1 # 0 for s  and 1 for p

+beam1.status = 'incoming' # could be reflected, refracted, incoming

+

+beam2 = beam()

+beam2.k = np.array([-np.cos(beamInitialTravelAngle), \

+                    -np.sin(beamInitialTravelAngle)])

+beam2.origin = np.array([xBeamOrigin, yBeamOrigin])

+beam2.face = float('NaN')

+beam2.intensity = 1

+beam2.polarization = 0 # 0 for s  and 1 for p

+beam2.status = 'incoming'

+

+beams = np.array([beam1, beam2])

+

+

+# image borders, we don't want a beam traced to infinity

+border1 = border()

+border1.vertex1 = np.array([-1.5, -0.5])

+border1.vertex2 = np.array([1.5, -0.5])

+border1.nRight = nAir

+border1.nLeft = nAir

+

+border2 = border()

+border2.vertex1 = np.array([1.5, -0.5])

+border2.vertex2 = np.array([1.5, 2.0])

+border2.nRight = nAir

+border2.nLeft = nAir

+

+border3 = border()

+border3.vertex1 = np.array([1.5, 2.0])

+border3.vertex2 = np.array([-1.5, 2.0])

+border3.nRight = nAir

+border3.nLeft = nAir

+

+border4 = border()

+border4.vertex1 = np.array([-1.5, 2.0])

+border4.vertex2 = np.array([-1.5, -0.5])

+border4.nRight = nAir

+border4.nLeft = nAir

+

+borders = np.array([border1, border2, border3, border4])

+borderLimits = np.array([-1.5, -0.5, 1.5, 2.0])

+# end parameters

+

+

+processedBeams = beamTrace(beams, faces, borders, borderLimits)

+pylab.show()

diff --git a/beam_tracing/python/faceBeamInteraction.py b/beam_tracing/python/faceBeamInteraction.py
new file mode 100644
index 0000000..a3b1d0b
--- /dev/null
+++ b/beam_tracing/python/faceBeamInteraction.py
@@ -0,0 +1,151 @@
+from classes import *

+from beam2FaceDistance import *

+from colorGradient import *

+from fresnelReflection import *

+import numpy as np

+import pylab

+

+def faceBeamInteraction(beamIn, faces):

+# calculates refracted and reflected beam after interaction with a face

+# beam - structure defining the light beam

+#    beam.origin - array [x,y] origin/soutce of the light beam

+#    beam.k - k vector i.e. direction [kx,ky]

+#    beam.intensity - intensity of the beam

+#    beam.face - if beam starts from face then its index is here

+#    beam.polarization - can be either 1 for 's' or 2 for 'p'

+#    beam.status -  could be 'reflected', 'refracted', 'incoming'

+# faces cell array of face structures

+# face - structure definiong the beam

+#    face.vertex2 - [x,y] of the 2nd point/vertex of the face

+#    face.n_left  - indexes of refraction on the left hand side 

+#                   with respect to 1st -> 2nd vertex direction

+#                   [ns,np] - for s and p polarizations

+#    face.n_right - indexes of refraction on the right hand side 

+#                   [ns,np] - for s and p polarizations

+

+   k = beamIn.k 

+   polarization = beamIn.polarization

+

+   # we go over all faces to find the closest which beam hits

+   nFaces = len(faces)

+   hitDistance = float('Inf')

+   isFaceHit = False

+   hitPosition = [float('NaN'), float('NaN')]

+   closestFaceIndex = float('NaN')

+   i = 0

+   for face in faces:

+      [hitDistanceTmp, hitPositionTmp, isFaceHitTmp] = \

+       beam2FaceDistance(beamIn, face)

+      if hitDistanceTmp < hitDistance:

+         # this is the closer face

+         isFaceHit = isFaceHitTmp

+	 hitPosition = hitPositionTmp

+	 hitDistance = hitDistanceTmp

+	 closestFaceIndex = i

+      i += 1

+   

+   if not isFaceHit:

+      newBeams = []

+      return [isFaceHit, hitPosition, hitDistance, newBeams]

+

+   # closest face

+   face = faces[closestFaceIndex]

+   kf = face.vertex2 - face.vertex1 # not a unit vector

+

+# hold on

+# draw face

+#line([face.vertex1(1), face.vertex2(1)], \

+#[face.vertex1(2), face.vertex2(2)] , 'linestyle', '-', 'color', [0,0,0])

+   # draw beam

+   nColors = 256

+   lineWidth = 2

+   if polarization == 0:

+      # s-polarization

+      lineBaseColor = [1,0,0] # RGB - red

+   else:

+      # p-polarization

+      lineBaseColor = [0,0,1] # RGB - blue

+   #plot(x,y,line_str)

+   lineColor = colorGradient(beamIn.intensity * nColors, lineBaseColor, nColors)

+   pylab.plot([beamIn.origin[0], hitPosition[0]],\

+              [beamIn.origin[1], hitPosition[1]],\

+              color = lineColor, linewidth = 2, linestyle = '-')

+

+# find is beam arriving from left or right using vectors cross product property

+# if z component of the product 'k x kf' is positive then beam arrives from 

+# the left

+   if (k[0]*kf[1] - k[1]*kf[0]) > 0:

+      # beam coming from the left

+      n1=face.nLeft[polarization]

+      n2=face.nRight[polarization]

+      # this means that notmal and beam k vector point to the same half plane

+      # relative to the face

+      nfAndkInTheSamePlane = True

+   else:

+      # beam coming from the right

+      n1 = face.nRight[polarization]

+      n2 = face.nLeft[polarization]

+      # this means that notmal and beam k vector point to the different half plane

+      # relative to the face

+      nfAndkInTheSamePlane = False

+

+   # normal vector to the face, looks to the left of it 

+   nf = [kf[1], -kf[0]] / np.linalg.norm(kf)

+   # incidence angle calculation

+   cosTheta_i = np.dot(k, nf) / (np.linalg.norm(k)*np.linalg.norm(nf))

+   sinTheta_i = -(k[0]*nf[1]-k[1]*nf[0]) /\

+              (np.linalg.norm(k) * np.linalg.norm(nf))

+   # positive angle to the right from normal before incidence to the face

+   theta_i = np.arctan2(sinTheta_i, cosTheta_i)

+

+   # we need to make sure that angle of incidence belong to [-pi/2,pi/2] interval

+   if theta_i > np.pi/2:

+      theta_i = np.pi-theta_i

+   if theta_i < -np.pi/2:

+      theta_i = -np.pi-theta_i

+	

+   # angle of the normal with respect to horizon

+   thetaNormal = np.arctan2(nf[1], nf[0])

+

+   # reflected beam direction

+   if nfAndkInTheSamePlane:

+      thetaReflected = thetaNormal - theta_i + np.pi

+   else:

+      thetaReflected = thetaNormal + theta_i

+

+   # coefficients of reflection and transmission for given polarization

+   [R, thetaRefractedRel2Normal] = fresnelReflection(n1, n2, theta_i)

+   reflectivity = R[polarization]

+   transmission = 1 - reflectivity

+	

+   beamReflected = beam()

+   beamReflected.origin = hitPosition

+   beamReflected.k = [np.cos(thetaReflected), np.sin(thetaReflected)]

+   beamReflected.face = closestFaceIndex

+   beamReflected.polarization = polarization

+   beamReflected.intensity = beamIn.intensity * reflectivity

+   beamReflected.status = 'reflected'

+   newBeams = [beamReflected]

+

+

+   # refracted beam direction

+   # refracted angle with respect to normal

+   if transmission < 0.001:

+      pass

+   else:

+      # beam refracts

+      if nfAndkInTheSamePlane:

+	 thetaRefracted = thetaNormal + thetaRefractedRel2Normal

+      else:

+	 thetaRefracted = thetaNormal - thetaRefractedRel2Normal + np.pi

+

+      beamRefracted = beam()

+      beamRefracted.origin = hitPosition

+      beamRefracted.k = [np.cos(thetaRefracted), np.sin(thetaRefracted)]

+      beamRefracted.face=closestFaceIndex

+      beamRefracted.polarization = polarization

+      beamRefracted.intensity = beamIn.intensity * transmission

+      beamRefracted.status = 'refracted'

+      newBeams.append(beamRefracted)

+      

+   return [isFaceHit, hitPosition, hitDistance, newBeams]
\ No newline at end of file
diff --git a/beam_tracing/python/fresnelReflection.py b/beam_tracing/python/fresnelReflection.py
new file mode 100644
index 0000000..2e4d4bc
--- /dev/null
+++ b/beam_tracing/python/fresnelReflection.py
@@ -0,0 +1,33 @@
+import numpy as np

+

+def fresnelReflection(n1, n2, thetaInc):

+# calculates intensity  reflection coefficient for s and p polarizations

+# for light travelling from material with index of refraction n1 to 

+# material with n2 

+# theta_i - incident angle in medium 1 with respect to normal, could be a vector

+# R       - coefficients of reflection array [Rs, Rp]

+# theta_t - transmitted/refracted angle in medium 2 with respect to normal

+

+   #if np.size(thetaInc, axis=0) != 1: 

+    #  error('thetaInc must be a vector or scalar')

+      

+   # see  http://en.wikipedia.org/wiki/Fresnel_equations

+   # refraction angle or angle of transmitted beam with respect to normal

+   sinThetaTrans = n1 / n2 * np.sin(thetaInc)

+   

+   # special cases: total internal reflection

+   if abs(sinThetaTrans) >= 1:

+      sinThetaTrans = np.copysign(1, sinThetaTrans)

+   thetaTrans = np.arcsin(sinThetaTrans) # angle of refraction

+

+

+   cosThetaTrans = np.cos(thetaTrans)

+   cosThetaInc = np.cos(thetaInc)

+

+   Rs = ((n1 * cosThetaInc - n2 * cosThetaTrans) /\

+         (n1 * cosThetaInc + n2 * cosThetaTrans)) ** 2

+   Rp = ((n1 * cosThetaTrans - n2 * cosThetaInc) /\

+         (n1 * cosThetaTrans + n2 * cosThetaInc)) ** 2

+   R = np.array([[Rs], [Rp]])

+

+   return [R, thetaTrans]
\ No newline at end of file
diff --git a/beam_tracing/python/gui_test.py b/beam_tracing/python/gui_test.py
new file mode 100644
index 0000000..f8edff4
--- /dev/null
+++ b/beam_tracing/python/gui_test.py
@@ -0,0 +1,11 @@
+import wx

+class MyFrame(wx.Frame):

+   """ We simply derive a new class of Frame. """

+   def __init__(self, parent, title):

+      wx.Frame.__init__(self, parent, title=title, size=(200,100))

+      self.control = wx.TextCtrl(self, style=wx.TE_MULTILINE)

+      self.Show(True)

+    

+app = wx.App(False)

+frame = MyFrame(None, 'Small editor')

+app.MainLoop()
\ No newline at end of file
diff --git a/beam_tracing/python/prismDiskCoupling.py b/beam_tracing/python/prismDiskCoupling.py
new file mode 100644
index 0000000..7cc9f70
--- /dev/null
+++ b/beam_tracing/python/prismDiskCoupling.py
@@ -0,0 +1,186 @@
+import numpy as np

+import pylab

+import sys

+

+

+def prismDiskCoupling(prismAngleInDegrees, nDisk, nPrism, diskX, coupDesc):

+   # Calculates incident angle for proper coupling into the disc via prism

+   # prismAngleInDegrees - angle of the prism faces in degrees

+   # diskX - x coordinate of the disk center

+   # coupDesc - short annotation of the situation

+

+   prismAngle = prismAngleInDegrees * np.pi / 180 # convert to radians

+

+   # critical angle for beam from prism to disk

+   # recall nDisk*sin(thetaDisk)=nPrism*sin(thetaPrism) where angles are

+   # counted from normal to the face and we want thetaDisk to be 90 degrees

+   # for total internal reflection, thus sin(thetaDisk)=1

+

+   thetaPrism1 = np.arcsin(nDisk / nPrism)

+   thetaPrismInDegrees1 = thetaPrism1 * 180 / np.pi

+

+   # find inside angle on the incident side of prism

+   thetaPrism2 = prismAngle - thetaPrism1

+   thetaPrismInDegrees2 = thetaPrism2 * 180 / np.pi #convert to degrees

+

+   # calculate refracted angle out of the prism with respect to the normal

+   # positive means above the normal; negative below

+   tempArcsin = nPrism * np.sin(thetaPrism2)

+

+   try:

+      thetaAir = np.arcsin(tempArcsin)

+   except ValueError:

+      print 'Error'

+      sys.exit('Total internal reflection at right prism face')

+   

+   thetaAir = np.arcsin(tempArcsin)

+   thetaAirInDegrees = thetaAir * 180 / np.pi #convert to degrees

+

+   # angle in the air relative to horizon

+   thetaAirHorizon = thetaAir + (np.pi/2 - prismAngle)

+   thetaAirHorizonInDegrees = thetaAirHorizon * 180 / np.pi #convert to degrees

+

+#   -----GRAPHICAL OUTPUT-----

+# PRISM

+   # bottom face will go from (-1,0) to (1,0) regardless of prismAngle

+   xFace1 = np.linspace(-1,1)

+   yFace1 = 0 * xFace1

+   

+   # second face to the right of origin 

+   # at prismAngle with respect to negative x-direction

+   xFace2 = np.linspace(1,0)

+   yFace2 = (1 - xFace2) * np.tan(prismAngle)

+   

+   # third face to the left of origin

+   # at prismAngle with respect to positive x-direction

+   xFace3 = np.linspace(-1,0)

+   yFace3 = (xFace3 + 1) * np.tan(prismAngle)

+

+   # draw prism

+   pylab.plot(xFace1, yFace1, 'black',\

+              xFace2, yFace2, 'black',\

+              xFace3, yFace3, 'black')

+

+# DISK 

+   # center will be located at (diskX, -radius)

+   radius = 0.25

+   diskCenterX = diskX

+   diskCenterY = -radius

+   xDisk = diskCenterX + (radius * np.cos(np.linspace(0, 2 * np.pi)))

+   yDisk = diskCenterY + (radius * np.sin(np.linspace(0, 2 * np.pi)))

+   pylab.plot(xDisk, yDisk, 'black')

+   

+

+# BEAMS

+   # Inside Prism

+   # drawing line from disk contact point to face2 at angle thetaPrism

+   xCross = (np.tan(prismAngle) + diskX / np.tan(thetaPrism1)) /\

+             (1/np.tan(thetaPrism1) + np.tan(prismAngle))

+   yCross = -np.tan(prismAngle) * (xCross - 1)

+   xBeamPrism = np.linspace(diskX, xCross)

+   yBeamPrism = np.linspace(0, yCross)

+   pylab.plot(xBeamPrism, yBeamPrism, color = 'red', linewidth = 2)

+   

+   # Outside Prism

+   xOutStart = xCross

+   yOutStart = yCross

+   xOutStop = 1.5

+   xBeamOut = np.linspace(xOutStart, xOutStop)

+   yBeamOut = yOutStart + np.tan(thetaAirHorizon) * (xBeamOut - xOutStart)

+   # this keeps the beam from going farther than necessary

+   i = 0

+   for value in yBeamOut:

+      if value > 2.0:

+         break

+      else:

+         i +=1

+   # ensure same dimension

+   yBeamOut = yBeamOut[:i]

+   xBeamOut = xBeamOut[:i]

+   pylab.plot(xBeamOut, yBeamOut, color = 'red', linewidth = 2)

+

+  

+   

+# NORMALS

+   # incident face normal

+   

+   thetaNorm = np.pi/2 - prismAngle

+   # coordinates of normal outside the prism

+   xNormOut1 = np.linspace(xCross, xCross + 0.25)

+   yNormOut1 = yCross + (xNormOut1 - xCross) * np.tan(thetaNorm)

+   # coordinates of normal inside the prism

+   xNormIn1 = np.linspace(xCross, xCross - 0.25)

+   yNormIn1 = yCross + (xNormIn1 - xCross) * np.tan(thetaNorm)

+   pylab.plot(xNormOut1, yNormOut1, 'b:', \

+              xNormIn1, yNormIn1, 'b:', linewidth = 1.5)

+   

+   # normal at the disk contact

+   yNorm2 = np.linspace(-0.25, 0.25)

+   xNorm2 = diskX + yNorm2 * 0.0

+   pylab.plot(xNorm2, yNorm2, 'b:', linewidth = 1.5)

+   

+   

+# ARCS AND LABELS

+   # radius for all drawn arcs

+   arcRadius = 0.1

+   

+   # arc to show thetaAir

+   # phiArcAir is the angle through which we are drawing

+   phiArcAir = np.linspace(thetaAirHorizon, thetaNorm)

+   xArcAir = xCross + arcRadius * np.cos(phiArcAir)

+   yArcAir = yCross + arcRadius * np.sin(phiArcAir)

+   pylab.plot(xArcAir, yArcAir, color = 'blue', linewidth = 1.0)

+   # annotatation of thetaAir

+   pylab.annotate('%.2f' % thetaAirInDegrees +  u'\N{DEGREE SIGN}',\

+                  xy = (xArcAir[29], yArcAir[29]),\

+                  xytext=(xCross, yCross + 0.3),\

+                  arrowprops = dict(facecolor = 'black', arrowstyle = '->'))

+

+   

+   # arc to show prismAngle

+   #phiArcPrism

+   phiArcPrism = np.linspace(np.pi, np.pi - prismAngle)

+   xArcPrism = 1 + arcRadius * np.cos(phiArcPrism)

+   yArcPrism = arcRadius * np.sin(phiArcPrism)

+   pylab.plot(xArcPrism, yArcPrism, color = 'blue', linewidth = 1.0)

+   # annotate prismAngle

+   pylab.annotate('%.2f' % prismAngleInDegrees + u'\N{DEGREE SIGN}',\

+                  xy = (xArcPrism[29], yArcPrism[29]),\

+                  xytext=(0.8, -0.2),\

+                  arrowprops = dict(facecolor = 'black', arrowstyle = '->'))   

+  

+

+   # arc to show thetaPrism1

+   # phiArcPrism1 is the angle through which we are drawing

+   phiArcPrism1 = np.linspace(np.pi/2 - thetaPrism1, np.pi/2)

+   xArcPrism1 = diskX + arcRadius * np.cos(phiArcPrism1)

+   yArcPrism1 = arcRadius * np.sin(phiArcPrism1)

+   pylab.plot(xArcPrism1, yArcPrism1, color = 'blue', linewidth = 1.0)

+   # annotatation of thetaPrism1

+   pylab.annotate('%.2f' % thetaPrismInDegrees1 + u'\N{DEGREE SIGN}',\

+                  xy = (xArcPrism1[29], yArcPrism1[29]),\

+                  xytext=(diskX - 0.2, 0.3),\

+                  arrowprops = dict(facecolor = 'black', arrowstyle = '->'))

+

+   

+   # arc to show thetaPrism2

+   # phiArcPrism2 is the angle through which we are drawing

+   phiArcPrism2 = np.linspace(thetaNorm + thetaPrism2 + np.pi,\

+                             thetaNorm + np.pi)

+   xArcPrism2 = xCross + arcRadius * np.cos(phiArcPrism2)

+   yArcPrism2 = yCross + arcRadius * np.sin(phiArcPrism2)

+   pylab.plot(xArcPrism2, yArcPrism2, color = 'blue', linewidth = 1.0)

+   # annotation of thetaPrism2

+   pylab.annotate('%.2f' % thetaPrismInDegrees2 + u'\N{DEGREE SIGN}',\

+                  xy = (xArcPrism2[29], yArcPrism2[29]),\

+                  xytext=(xCross - 0.5, yCross + 0.1),\

+                  arrowprops = dict(facecolor = 'black', arrowstyle = '->'))

+   

+   

+# General annotation and axis set-up

+   pylab.text(-1.4, 1.7, coupDesc)

+   pylab.text(-1.4, 1.55, 'nDisk = %.3f' % nDisk)

+   pylab.text(-1.4, 1.4, 'nPrism = %.3f' % nPrism)   

+   pylab.xlim(-1.5, 1.5)

+   pylab.ylim(-0.5, 2.0)

+   pylab.title(coupDesc)
\ No newline at end of file