Source code
%% ScriptSolveForNearField
% Example script showing how to obtain the field expansion coefficients,
% far-field cross-sections and surface field properties for a spheroid in a
% fixed orientation and at a single wavelength.
% Outputs the cross-sections and surface-averaged properties with accuracy
% estimates.
% Plots the field enhancement factors as a function of theta for three values of phi.
% Also produces a 3D plot of the surface field intensity on the particle.
%%
%% Initialization
%
% Note that you need to run InitPath in the root folder first to add
% required folders to the Matlab path so that functions can be called
% Alternatively, uncomment the following line
%
% run('..\InitPath');
%
% The following parameters should be defined:
%%
% * a: semi-axis along x,y
% * c: semi-axis along z
% * k1: wavevector in embedding medium (of refractive index nM) (k1=2*pi*nM/lambda)
% * s: relative refractive index (s=n_Particle / nM)
% * N: number of multipoles for T-matrix
% * nNbTheta: number of thetas for quadratures
% * sIncType: string such a 'KxEz' defining the incident field)
% * or stIncPar: struct defining the incident field (alternative
% to sIncType)
% * nNbThetaPst: number of thetas for surface field calculations
clear
close all
%% Parameters of the scattering problem
% We define parameters for a gold nanodisc in water, modeled as an oblate
% spheroid
%
% <<../fig/schematico.png>>
%
c = 15; % in nm
a = 45; % in nm, i.e. 20 x 100nm full-axes
lambda = 650; % in nm
epsilon2 = epsAu(lambda);
epsilon1 = 1.33^2; % for water
% Define incident excitation along main axis
sIncType = 'KxEz';
%% Convergence parameters
% Maximum multipole order for T-matrix and series expansions of fields
N = 40;
% Number of points for Gaussian quadratures to compute integrals in P and Q matrices
nNbTheta = 120;
% Number of points for post-processing (computing the surface fields
% averages)
nNbThetaPst = 360;
%% Collect simulation parameters in a structure
k1 = 2*pi/lambda * sqrt(epsilon1);
s = sqrt(epsilon2)/sqrt(epsilon1);
stParams.a=a; stParams.c=c;
stParams.k1=k1; stParams.s=s;
stParams.N=N; stParams.nNbTheta=nNbTheta;
stParams.sIncType = sIncType;
% For surface fields, the following parameters are also needed:
stParams.lambda = lambda;
stParams.epsilon2= epsilon2;
stParams.epsilon1= epsilon1;
stParams.nNbThetaPst = nNbThetaPst;
% Optional parameters may also be defined as follows:
stOptions.bGetR = true; % This is needed for near fields and will be overridden in any case
stOptions.Delta = 0; % Use Delta=-1 to estimate Delta automatically
stOptions.NB = 0; % NB will be estimated automatically
stOptions.bGetSymmetricT = false;
%% T-matrix calculation
tic;
[stC, stAbcdnm, stEsurf] = slvForNearField(stParams,stOptions);
fprintf('\nT/R-matrices and near fields (N = %d) ... done in %.f seconds.\n', N, toc);
% To test for convergence and accuracy for a given set of parameters, one
% can for example repeat the calculation with N=N+5 and nNbTheta=nNbTheta+5
% as illustrated below
fprintf('Convergence testing...');
tic;
stParams2=stParams;
stParams2.N=stParams2.N+5;
stParams2.nNbTheta=stParams2.nNbTheta+5;
[stC2, stAbcdnm2, stEsurf2] = slvForNearField(stParams2,stOptions);
fprintf('\nT/R-matrices and near fields (N = %d) ... done in %.f seconds.\n', N, toc);
fprintf('Results for a=%g, c=%g, k1=%g, s=%g+%gi, N=%d, Nt=%d\n',...
a, c, k1, real(s),imag(s), N, nNbTheta);
fprintf('\nCross sections for fixed excitation (and estimated accuracy):\n');
fprintf('Cext = %.20g, relative error: %.2g\n', stC.Cext, abs(stC.Cext./stC2.Cext-1))
fprintf('Csca = %.20g, relative error: %.2g\n', stC.Csca, abs(stC.Csca./stC2.Csca-1))
fprintf('Cabs = %.20g, relative error: %.2g\n', stC.Cabs, abs(stC.Cabs./stC2.Cabs-1))
fprintf('\nCross sections for orientation-averaged excitation (and estimated accuracy):\n');
fprintf(' = %.20g, relative error: %.2g\n', stC.Cextoa, abs(stC.Cextoa./stC2.Cextoa-1))
fprintf(' = %.20g, relative error: %.2g\n', stC.Cscaoa, abs(stC.Cscaoa./stC2.Cscaoa-1))
fprintf(' = %.20g, relative error: %.2g\n', stC.Cabsoa, abs(stC.Cabsoa./stC2.Cabsoa-1))
fprintf('\nSurface-averaged surface-field properties (accuracy not tested):\n');
fprintf('<|E|^2> = %.20g\n', stEsurf.MLocAve);
fprintf('<|E_{perp}|^2> = %.20g\n', stEsurf.MLocPerpAve);
fprintf('<|E|^4> = %.20g\n', stEsurf.F0E4Ave);
%% Examples of postprocessing for surface fields
%
% The following gets the surface field intensity enhancement factor (M) at
% one or more given phi (as a function of theta)
phi0=[0,pi/4,pi/2];
M = pstGetThetaDepFieldIntensity(stEsurf,phi0); % [3 x T]
M2 = pstGetThetaDepFieldIntensity(stEsurf2,phi0); % [3 x T]
relerrM = max(abs(M-M2),[],2)./max(abs(M2),[],2); % [3 x 1]
% Make plots to show results
figure('Name','Theta-dependence of surface-field intensity M=|E|^2 for fixed phi');
plot(stEsurf.theta, M);
legend({['phi=', num2str(phi0(1)), ' (err=', num2str(relerrM(1)),')'], ...
['phi=', num2str(phi0(2)), ' (err=', num2str(relerrM(2)),')'], ...
['phi=', num2str(phi0(3)), ' (err=', num2str(relerrM(3)),')']}, ...
'Location', 'Best');
xlabel('Theta [radians]')
ylabel('M=|E/E0|^2')
% The following makes a 3D surface plot of the surface
% field everywhere on the surface
% NOTE that this requires to recompute all
% the surface fields from scratch. Using 90x90 pts here:
tic;
stResE=pstMakeStructForField(stAbcdnm, stParams);
pstPlotAllSurfaceField(90, stResE);
fprintf('\n3D plot ... done in %.f seconds.\n', toc);
Execution results
octave>ScriptSolveForNearField
sphCalculatePQ: Calculating P,Q for 41 m-values with N_Q = 40, N_B = 41, N_Theta = 120
T/R-matrices and near fields (N = 40) ... done in 3 seconds.
Convergence testing...sphCalculatePQ: Calculating P,Q for 46 m-values with N_Q = 45, N_B = 45, N_Theta = 125
T/R-matrices and near fields (N = 40) ... done in 3 seconds.
Results for a=45, c=15, k1=0.0128564, s=0.119065+2.67664i, N=40, Nt=120
Cross sections for fixed excitation (and estimated accuracy):
Cext = 162.29599945903115099, relative error: 2e-015
Csca = 88.598051846348710114, relative error: 8.9e-016
Cabs = 73.697947612682440877, relative error: 3.6e-015
Cross sections for orientation-averaged excitation (and estimated accuracy):
= 23273.007567191281851, relative error: 5e-015
= 13699.998155139272058, relative error: 4.8e-015
= 9573.0094120520097931, relative error: 5.3e-015
Surface-averaged surface-field properties (accuracy not tested):
<|E|^2> = 2.4940806742084613568
<|E_{perp}|^2> = 2.4537899387743351554
<|E|^4> = 7.1825788909828434115
warning: legend: 'best' not yet implemented for location specifier
3D plot ... done in 1 seconds.
Generated graphics

