Added python version of mode-matching code by Bain Bronner

1 files changed, 186 insertions, 0 deletions

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