summaryrefslogtreecommitdiff
path: root/beam_tracing/python/prismDiskCoupling.py
blob: aff59ff987adb1497fffcf6726bdf546739a216f (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
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
   yOutStop = 2.0
   yBeamOut = np.linspace(yOutStart, yOutStop)
   xBeamOut = xOutStart + (1/np.tan(thetaAirHorizon)) * (yBeamOut - yOutStart)
   # 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 + 0.12, 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.2, yCross + 0.5),\
                  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)