diff options
| -rw-r--r-- | beam_tracing/python/README.txt | 20 | ||||
| -rw-r--r-- | beam_tracing/python/TODO.txt | 7 | ||||
| -rw-r--r-- | beam_tracing/python/beam2FaceDistance.py | 53 | ||||
| -rw-r--r-- | beam_tracing/python/beamTrace.py | 49 | ||||
| -rw-r--r-- | beam_tracing/python/classes.py | 37 | ||||
| -rw-r--r-- | beam_tracing/python/colorGradient.py | 32 | ||||
| -rw-r--r-- | beam_tracing/python/demoPrismDiskCoupling.py | 21 | ||||
| -rw-r--r-- | beam_tracing/python/drawPrismDisk.py | 102 | ||||
| -rw-r--r-- | beam_tracing/python/exampleBeamTrace.py | 122 | ||||
| -rw-r--r-- | beam_tracing/python/faceBeamInteraction.py | 151 | ||||
| -rw-r--r-- | beam_tracing/python/fresnelReflection.py | 33 | ||||
| -rw-r--r-- | beam_tracing/python/gui_test.py | 11 | ||||
| -rw-r--r-- | beam_tracing/python/prismDiskCoupling.py | 186 |
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 |