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