1 files changed, 151 insertions, 0 deletions
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