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/README.txt | 20 +++ beam_tracing/python/TODO.txt | 7 + beam_tracing/python/beam2FaceDistance.py | 53 ++++++++ beam_tracing/python/beamTrace.py | 49 +++++++ beam_tracing/python/classes.py | 37 ++++++ beam_tracing/python/colorGradient.py | 32 +++++ beam_tracing/python/demoPrismDiskCoupling.py | 21 +++ beam_tracing/python/drawPrismDisk.py | 102 +++++++++++++++ beam_tracing/python/exampleBeamTrace.py | 122 ++++++++++++++++++ beam_tracing/python/faceBeamInteraction.py | 151 ++++++++++++++++++++++ beam_tracing/python/fresnelReflection.py | 33 +++++ beam_tracing/python/gui_test.py | 11 ++ beam_tracing/python/prismDiskCoupling.py | 186 +++++++++++++++++++++++++++ 13 files changed, 824 insertions(+) create mode 100644 beam_tracing/python/README.txt create mode 100644 beam_tracing/python/TODO.txt create mode 100644 beam_tracing/python/beam2FaceDistance.py create mode 100644 beam_tracing/python/beamTrace.py create mode 100644 beam_tracing/python/classes.py create mode 100644 beam_tracing/python/colorGradient.py create mode 100644 beam_tracing/python/demoPrismDiskCoupling.py create mode 100644 beam_tracing/python/drawPrismDisk.py create mode 100644 beam_tracing/python/exampleBeamTrace.py create mode 100644 beam_tracing/python/faceBeamInteraction.py create mode 100644 beam_tracing/python/fresnelReflection.py create mode 100644 beam_tracing/python/gui_test.py create mode 100644 beam_tracing/python/prismDiskCoupling.py 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 -- cgit v1.2.3