From c2b9aa353d8f2f77ac0cd00663a8d6f633fcb7bd Mon Sep 17 00:00:00 2001 From: dmayerich Date: Tue, 23 Apr 2013 14:41:49 -0500 Subject: [PATCH] Added scattering amplitude calculation. --- msinput.inp | 2 +- mstm-gui.py | 36 +++++++++++++----------------------- mstm-main-v2.2.f90 | 60 +++++++++++++++++++++++++++++++++++++++++++++--------------- mstm_simparser.py | 50 +++++++++++++++++++++++++++++++++++++++++++++----- 4 files changed, 104 insertions(+), 44 deletions(-) diff --git a/msinput.inp b/msinput.inp index 5012c29..6abca35 100755 --- a/msinput.inp +++ b/msinput.inp @@ -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 diff --git a/mstm-gui.py b/mstm-gui.py index 36a2a25..04f787e 100755 --- a/mstm-gui.py +++ b/mstm-gui.py @@ -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 diff --git a/mstm-main-v2.2.f90 b/mstm-main-v2.2.f90 index 3f82e92..0f68fbb 100755 --- a/mstm-main-v2.2.f90 +++ b/mstm-main-v2.2.f90 @@ -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'')') diff --git a/mstm_simparser.py b/mstm_simparser.py index 9babb7d..5fd282e 100755 --- a/mstm_simparser.py +++ b/mstm_simparser.py @@ -1,3 +1,4 @@ +#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 -- libgit2 0.21.4