Added scattering amplitude calculation.

dmayerich
1 parent c0a2a7bb
Showing 4 changed files with 104 additions and 44 deletions Show diff stats
msinput.inp
mstm-gui.py
mstm-main-v2.2.f90
mstm_simparser.py
@@ -39,7 +39,7 @@ max_scattering_angle_deg
 number_scattering_angles
 181
 normalize_scattering_matrix
-1
+0
 gaussian_beam_constant
 0.0d0
 gaussian_beam_focal_point
@@ -98,6 +98,7 @@ class GuiWindow(QtGui.QMainWindow):
 		#plot results of interest
 		wl = self.results['lambda']
+		#for a fixed orientation
 		if int(self.params['fixed_or_random_orientation']) == 0:
 			unpol = self.results['extinction_unpolarized']
 			para = self.results['extinction_parallel']
@@ -132,33 +133,19 @@ class GuiWindow(QtGui.QMainWindow):
 		self.results = RunSimulation(False)
 		if self.params['calculate_near_field']:
-			#verts = self.params['near_field_plane_vertices']
-			#dx = (verts[2] - verts[0])/(self.params.nSteps)
-			#x = arange(verts[0], verts[2], dx)
-			#print(len(x))
-			#y = arange(verts[1], verts[3], dx)
-			#X, Y = meshgrid(x, y)
+
 			E = array(self.results.gridNearField)
-			#pcolor(X, Y, E, cmap=cm.RdBu)
-			#colorbar()
-			#axis([verts[0], verts[2], verts[1], verts[3]])
+
 			pcolor(E, cmap=cm.RdBu)
 			colorbar()
+			
+			#compute the maximum field magnitude in the plane
 			print("Maximum enhancement: " + str(abs(E).max()))
+			
+			#get the scattering amplitude matrix
+			print(self.results.scatAmpMatrix[0])
-		# make these smaller to increase the resolution
-		#dx, dy = 0.05, 0.05
-
-		#x = arange(-3.0, 3.0001, dx)
-		#y = arange(-3.0, 3.0001, dy)
-		#X,Y = meshgrid(x, y)
-
-		#Z = self.func3(X, Y)
-		#pcolor(X, Y, Z, cmap=cm.RdBu, vmax=abs(Z).max(), vmin=-abs(Z).max())
-		#colorbar()
-		#axis([-3,3,-3,3])
-
 		show()
 	def saveresults(self):
@@ -184,12 +171,12 @@ class GuiWindow(QtGui.QMainWindow):
 		a = self.ui.spinRadius.value()
 		self.ui.tblSpheres.setItem(0, 0, QtGui.QTableWidgetItem(str(a)))
-		self.ui.tblSpheres.setItem(0, 1, QtGui.QTableWidgetItem(str(-(d + 2*a)/2)))
+		self.ui.tblSpheres.setItem(0, 1, QtGui.QTableWidgetItem(str(-(d/2.0 + a))))
 		self.ui.tblSpheres.setItem(0, 2, QtGui.QTableWidgetItem(str(0.0)))
 		self.ui.tblSpheres.setItem(0, 3, QtGui.QTableWidgetItem(str(0.0)))
 		self.ui.tblSpheres.setItem(1, 0, QtGui.QTableWidgetItem(str(a)))
-		self.ui.tblSpheres.setItem(1, 1, QtGui.QTableWidgetItem(str((d + 2*a)/2)))
+		self.ui.tblSpheres.setItem(1, 1, QtGui.QTableWidgetItem(str((d/2.0 + a))))
 		self.ui.tblSpheres.setItem(1, 2, QtGui.QTableWidgetItem(str(0.0)))
 		self.ui.tblSpheres.setItem(1, 3, QtGui.QTableWidgetItem(str(0.0)))
@@ -315,6 +302,9 @@ def RunSimulation(spectralSim = True):
 		if parameters['calculate_near_field']:
 			results.parseNearField('nf-temp.dat')
+			
+		#get the scattering amplitude matrix
+		results.calcScatteringAmp()
 		#update the progress bar
@@ -348,31 +348,61 @@
             write(1,'('' perpendicular total ext, abs, scat efficiencies'')')
             write(1,'(6e13.5)') qexttotper,qabstotper,qscatotper
          endif
-
+         
          write(1,'('' scattering matrix elements'')')
          if(normalizesm.eq.0) then
-            write(1,'('' theta    s11         s22         s33'',&
-                   &''         s44'',&
-                   &''         s21         s32         s43         s31'',&
-                   &''         s42         s41'')')
+            write(1,'('' theta    s11         s12         s13'',&
+                   &''         s14'',&
+                   &''         s21         s22         s23         s24'',&
+                   &''         s31         s32         s33         s34'',&
+                   &''         s41         s42         s43         s44'')')
             do i=1,numtheta
                thetad=theta1d+(theta2d-theta1d)*(i-1)/max(1.d0,dble(numtheta-1))
-               write(1,'(f8.2,10e12.4)') thetad,smt(1,1,i),smt(2,2,i),smt(3,3,i), &
-                     smt(4,4,i),smt(1,2,i),smt(2,3,i),smt(3,4,i),smt(1,3,i), &
-                     smt(2,4,i),smt(1,4,i)
+               write(1,'(f8.2,16e12.4)') thetad,smt(1,1,i),smt(2,1,i),smt(3,1,i), &
+                     smt(4,1,i),smt(1,2,i),smt(2,2,i),smt(3,2,i),smt(4,2,i), &
+                     smt(1,3,i),smt(2,3,i),smt(3,3,i),smt(4,3,i),smt(1,4,i), &
+                     smt(2,4,i),smt(3,4,i),smt(4,4,i)
+
+!         write(1,'('' scattering matrix elements'')')
+!         if(normalizesm.eq.0) then
+!            write(1,'('' theta    s11         s22         s33'',&
+!                   &''         s44'',&
+!                   &''         s21         s32         s43         s31'',&
+!                   &''         s42         s41'')')
+!            do i=1,numtheta
+!               thetad=theta1d+(theta2d-theta1d)*(i-1)/max(1.d0,dble(numtheta-1))
+!               write(1,'(f8.2,10e12.4)') thetad,smt(1,1,i),smt(2,2,i),smt(3,3,i), &
+!                     smt(4,4,i),smt(1,2,i),smt(2,3,i),smt(3,4,i),smt(1,3,i), &
+!                     smt(2,4,i),smt(1,4,i)
             enddo
          else
-            write(1,'('' theta    s11         s22/s11     s33'',&
-                   &''/s11     s44'',&
-                   &''/s11     s21/s11     s32/s11     s43/s11     s31'',&
-                   &''/s11     s42/s11     s41/s11'')')
+!This code was changed by David Mayericha and Thomas van Dijk
+            write(1,'('' theta    s11         s12/s11     s13'',&
+                   &''/s11     s14'',&
+                   &''/s11     s21/s11     s22/s11     s23/s11     s24'',&
+                   &''/s11     s31/s11     s32/s11     s33/s11     s34'',&
+                   &''/s11     s41/s11     s42/s11     s43/s11     s44/s11'')')
             do i=1,numtheta
                thetad=theta1d+(theta2d-theta1d)*(i-1)/max(1.d0,dble(numtheta-1))
                s11=smt(1,1,i)
-               write(1,'(f8.2,10e12.4)') thetad,smt(1,1,i),smt(2,2,i)/s11,smt(3,3,i)/s11, &
-                     smt(4,4,i)/s11,smt(1,2,i)/s11,smt(2,3,i)/s11,smt(3,4,i)/s11,smt(1,3,i)/s11, &
-                     smt(2,4,i)/s11,smt(1,4,i)/s11
+               write(1,'(f8.2,16e12.4)') thetad,smt(1,1,i),smt(2,1,i)/s11,smt(3,1,i)/s11, &
+                     smt(4,1,i)/s11,smt(1,2,i)/s11,smt(2,2,i)/s11,smt(3,2,i)/s11,smt(4,2,i)/s11, &
+                     smt(1,3,i)/s11,smt(2,3,i)/s11,smt(3,3,i)/s11,smt(4,3,i)/s11,smt(1,4,i)/s11, &
+                     smt(2,4,i)/s11,smt(3,4,i)/s11,smt(4,4,i)/s11
             enddo
+
+!original code here
+!            write(1,'('' theta    s11         s22/s11     s33'',&
+!                   &''/s11     s44'',&
+!                   &''/s11     s21/s11     s32/s11     s43/s11     s31'',&
+!                   &''/s11     s42/s11     s41/s11'')')
+!            do i=1,numtheta
+!               thetad=theta1d+(theta2d-theta1d)*(i-1)/max(1.d0,dble(numtheta-1))
+!               s11=smt(1,1,i)
+!               write(1,'(f8.2,10e12.4)') thetad,smt(1,1,i),smt(2,2,i)/s11,smt(3,3,i)/s11, &
+!                     smt(4,4,i)/s11,smt(1,2,i)/s11,smt(2,3,i)/s11,smt(3,4,i)/s11,smt(1,3,i)/s11, &
+!                     smt(2,4,i)/s11,smt(1,4,i)/s11
+!            enddo
          endif
          if(fixedorrandom.eq.1) then
             write(1,'('' scattering matrix expansion coefficients'')')
+#this code parses the results of a simulation and stores them in a SimParserClass structure
 from pylab import *
 class SimParserClass:
@@ -12,6 +13,10 @@ class SimParserClass:
     gridNearField = []
     maxNearField = []
+    #the stokes matrix is read from the output
+    stokesMatrix = []
+    scatAmpMatrix = []
+    
     def __init__(self, parameters):
@@ -23,6 +28,7 @@ class SimParserClass:
         self.simResults['extinction_parallel'] = list()
         self.simResults['extinction_perpendicular'] = list()
         self.simResults['extinction_total'] = list()
+        self.simResults['detector_field'] = list()
         self.gridNearField = []
         self.maxNearField = []
@@ -36,10 +42,17 @@ class SimParserClass:
         while True:
             line = inFile.readline().strip()
-
+			
+			#if the simulation is for a single plane wave
             if int(self.params['fixed_or_random_orientation']) == 0:
                 if line == 'scattering matrix elements':
-                    break
+					#empty the stokes matrix
+					self.stokesMatrix = []
+					inFile.readline()
+					for s in range(0, 181):
+						values = map(float, inFile.readline().strip().split())
+						self.stokesMatrix.append(values)
+					break;
                 elif line == 'unpolarized total ext, abs, scat efficiencies, w.r.t. xv, and asym. parm':
                     values = inFile.readline().strip().split(' ')
                     self.simResults['extinction_unpolarized'].append(values[0])
@@ -49,7 +62,8 @@ class SimParserClass:
                 elif line == 'perpendicular total ext, abs, scat efficiencies':
                     values = inFile.readline().strip().split(' ')
                     self.simResults['extinction_perpendicular'].append(values[0])
-            
+                
+            #if the simulation is for random orientations
             else:
                 if line == 'scattering matrix elements':
                     break
@@ -83,6 +97,31 @@ class SimParserClass:
         E = array(self.gridNearField)
         self.maxNearField.append(abs(E).max())
+        
+    #calculate and return the scattering amplitude matrix
+    def calcScatteringAmp(self):
+		#compute the number of entries in the stokes matrix
+		nEntries = len(self.stokesMatrix)
+		
+		#initialize the scattering amplitude matrix to empty
+		self.scatAmpMatrix = []
+		
+		for s in range(0, nEntries):
+			Z = self.stokesMatrix[s]
+			
+			scatEntry = []
+			s11 = complex(sqrt(0.5 * (Z[1] - Z[2] - Z[5] + Z[6])), 0.0)
+			scatEntry.append(s11)
+			scatEntry.append(complex(-0.5 * (Z[3] + Z[7]) / s11, 0.5 * (Z[4] + Z[8]) / s11))
+			scatEntry.append(complex(-0.5 * (Z[9] + Z[10]) / s11, -0.5 * (Z[13] + Z[14]) / s11))
+			scatEntry.append(complex(0.5 * (Z[11] + Z[12]) / s11, -0.5 * (Z[12] - Z[15]) / s11))
+			
+			self.scatAmpMatrix.append(scatEntry)
+			
+		S = self.scatAmpMatrix[0]
+		E = [S[0], S[2]]
+		self.simResults['detector_field'].append(E)
+		print(E)
     def saveFile(self, fileName):
@@ -101,10 +140,11 @@ class SimParserClass:
             result += '\t' + str(self.simResults['extinction_unpolarized'][i])
             result += '\t' + str(self.simResults['extinction_parallel'][i])
             result += '\t' + str(self.simResults['extinction_perpendicular'][i])
+            result += '\t' + str(self.simResults['detector_field'][i][0]) + '\t' + str(self.simResults['detector_field'][i][1])
             #parse the near field if it is included in the simulation
-            if int(parameters['calculate_near_field']) == 1:
-                result += '\t' + str(maxNearField)
+            #if int(parameters['calculate_near_field']) == 1:
+            #    result += '\t' + str(maxNearField)
             result += '\n'
         return result