Post-processing raw data using Matlab¶
This notebook introduces how to process the raw data from SIE calculations to find optical forces and additional physical information such as near-fields, surface charges, and scattering cross-sections. The raw data we process has mainly two types:
- Surface currents on the edges of a simulation mesh
- Electric field amplitudes evaluated from surface currents
For more information about the raw data itself, you can check the other notebook about raw data visualization.
0. Simulation and post-processing flow¶
Once we get simulation results from the SIE calculations, we process the raw data using Matlab. The image below gives an idea of how the raw data are processed to find the physical quantities that we are interested in. The flow in the upper branch gives results about the antenna itself (e.g., near-fields, scattering cross-sections) and is used to calculate analytical forces. The lower branch includes the particle geometry so that we can get full-wave solutions to find the optical force numerically as well as surface charges and near-fields when the particle geometry is taken into account.
0-1. Matlab file tree¶
You can check the folder tree in the PostProcessing folder with the cell below. Inside the folder PostProcessing, you find main scripts and a subfolder named functions for functions used in the main scripts.
from python.folderTree import list_files, list_folders
list_files('PostProcessing')
1. Optical force calculation using Matlab¶
1-1. The dipole approaximation (DA): force_DA.m¶
The analytical solution for the optical force assumes the particle as a dipole so that we don't need to put a particle geometry into the field calculation. Using the electric fields evaluated from the sole antenna, we perform the force computation using the force_DA.m file. We take the raw data in the folder path, RawData/Antenna/FieldEvaluation/..., define parameters and variables, and find the analytical solution for the optical force with the line below.
F = analyticalForceCalc(fdPath, wavelength, diameter, position, delta, lightIntensity);
1-2. Maxwell's stress tensor (MST): force_MST.m¶
The numerical solution for the optical force includes the particle geometry in the SIE calculation. The surface currents are first obtained with SIE, and we directly use the surface currents to evaluate the Maxwell's stress tensor and the optical force using the force_MST.m file.
Since the calculation includes the particle geometry, we have different meshes when changing the particle size or position. Depending on the parameter we choose to sweep, the folder structure differs as you can check in the following cell.
list_folders('RawData/Antenna+Particle/Simulation')
When we scan the position of the particle, we first defined the position vector, set the data path RawData/Antenna+Particle/Simulation/movingParticle/..., and do the force calculation with:
[Fx,Fy,Fz] = forceP(fdPath, domain_num, position, lightIntensity);
When we choose the size of the particle and the incident wavelength to be swept, we first decide the illumination direction (either the trasmission mode or the reflection mode), set the datapath, RawData/Antenna+Particle/Simulation/fixedParticle/sizeWavelengthSweep/..., and do the force calculation with the line:
[Fx,Fy,Fz,D,L]= forceDL(fdPath, domain_num, wavelength, lightIntensity);
The input domain_num indicates the domain index of the particle where we want to find the force by integrating the Maxwell's stress tensor on it. The function outputs force arrays in the x, y, and z directions, which can be either vector with a single variable or two-dimensional arrays with two variables.
The functions forceP and forceDL deal with different data structures, but they basically perform the same calculation with force.m function, which computes the time-averaged optical force on the particle using the stress tensor as introduced in the MST method in the paper.
2. Near-field plot: nearfield.m¶
The Matlab script nearfield.m plots electric field amplitudes that we get from the SIE calculations. In the paper, we have chosen four different parameter sets as representatives of each coupling regime. In the folder RawData/Antenna+Particle/FieldEvaluation/, we have evaluated fields for these conditions and put them into subfolders like below:
list_folders('RawData/Antenna+Particle/FieldEvaluation')
You can draw the near-field plots for all these conditions using the line:
plotNearfield(fdPath, letters, L, D, plt, save)
The function plots four figures as follows:
 |
 |
 |
 |
3. Surface polarization charges: surfacePolarizationCharge.m¶
Similar to the near-field plot. the script surfacePolarizationCharge.m plots surface charges for the chosen parameter sets. The syntax is also similar:
plotSurfaceCharge(fdPath, Pname, L, D, plt, save)
The image below shows an example:
4. Dipole moments: dipoleMoment.m¶
The electrical dipole moments of the studied structure can be evaluated with the dipoleMoment.m file. The core of the calculation is in the momentEall.m function. The syntax of this function is:
[ED, EQ, EO] = momentEall(mesh, solfilename, lambda, epsr, cm)
It takes mesh and solution file names as inputs as well as the wavelength, the relative permittivities of all domains, and the center of mass. Then it outputs total electrical moments, i.e., the electrical moments of the whole structure based on surface charges. We can use this function to calculate electrical moments for each domain; this can be done by setting the epsilon values (the relative permittivities) of domains as that of the background except for the domain of interest. The center of mass should also be set for the domain of interest. By doing so, we can separately compute the electrical moments only for the nanoparticle and the antenna.
5. Absorption cross-sections: absCrossSection.m¶
The absorption cross-sections of the studied structure can be evaluated with the absCrossSection.m file. Inside the main file, the function calcCabs.m computes the absorption cross-sections and the wavelength vector. The syntax of the function is as following:
[Cabs, wavelength] = calcCabs(fdPath, num_sol, num_domain)
The input fdPath indicates the folder path where we find the job, mesh, and solution files. The function cabsall.m finds the absorped energy per second (i.e., the absorbed power) of all domains. The output Cabs is a matrix whose first column gives the total absorption cross-secions and the rest of the columns gives the absorption cross-sections of each domain in the order of the domain indices. The unit is in $nm^2$.
6. Scattering cross-sections: scatCrossSection.m¶
To evaluate scattering cross-sections, we first need to find far-field amplitudes. In the folder RawData/Antenna/FieldEvaluation/farfield/, you can find far-fields for the antenna used in the paper. The fields are evaluated on a 50 micron-radius sphere. Then we find the total energy scattered per second from this sphere by integrating the Poynting vectors on the surface. In the subfolder scatteringCrossSection, the file scat.sh does this job. The output scat.dat has two columns: the first column shows the wavelength, and the second one shows the total energy scattered per second. By dividing it with the energy incident per unit area per second gives the cross-section for scattering. The Matlab script scatCrossSection.m computes the incident power and calculates the scattering cross-section from the original scat.dat file. The unit is in $nm^2$.