From e4360a48435cc762855d356b208b5544e6f40c4b Mon Sep 17 00:00:00 2001 From: Eugeniy Mikhailov Date: Thu, 28 Feb 2013 14:59:55 -0500 Subject: Added python version of mode-matching code by Bain Bronner --- beam_tracing/python/prismDiskCoupling.py | 186 +++++++++++++++++++++++++++++++ 1 file changed, 186 insertions(+) create mode 100644 beam_tracing/python/prismDiskCoupling.py (limited to 'beam_tracing/python/prismDiskCoupling.py') 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 -- cgit v1.2.3