Source code
%% ScriptSolveForNearFieldSpectrum
% Example script showing how to obtain the field expansion coefficients,
% far-field cross-sections and surface field properties for a spheroid in
% fixed orientation, as a function of wavelength.
% Plots the wavelength-dependent spectra for extinction, scattering, and
% absorption cross-sections (fixed orientation as well as orientation-averaged),
% along with lambda-dependent surface-averaged surface field properties.
%%
%% 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
% * N: number of multipoles for T-matrix
% * nNbTheta: number of thetas for quadratures
% * lambda: wavelength (in same unit as a and c)
% * k1: wavevector in embedding medium (of refractive index nM) (k1=2*pi*nM/lambda)
% * s: relative refractive index (s=n_Particle / nM)
% lambda, k1, and s must be vectors of identical size [L x 1]
clear
close all
%% Parameters of the scattering problem
% We define parameters for a gold nanorod in water, modeled as a prolate
% spheroid
%
% <<../fig/schematicp.png>>
%
a = 15; % in nm
c = 45; % in nm, i.e. 30 x 90nm full-axes
lambda = (400:5:900).'; % 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 = 20;
% Number of points for Gaussian quadratures to compute integrals in P and Q matrices
nNbTheta = 50;
% Number of points for post-processing (computing the surface field 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.lambda=lambda;
stParams.sIncType = sIncType;
% For surface fields, the following parameters are also needed:
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;
stOptions.bOutput = false; % false to suppress messages in lambda-loop
%% Solving for the T-matrix (all wavelengths)
tic;
[stC, stAbcdnm, stEsurf] = slvForNearFieldSpectrum(stParams,stOptions);
fprintf('\nT/R-matrices and near fields (N = %d) ... done in %.g seconds.\n', N, toc);
% To test for convergence and accuracy, we choose the wavelength with the largest
% k1|s| and repeat the calculation with N=N+5 and nNbTheta=nNbTheta+5
[~,indWorst]=max(abs(stParams.k1 .* stParams.s));
stParams2 = pstGetParamsStructOneLambda(stParams,lambda(indWorst));
stParams2.N=stParams2.N+5;
stParams2.nNbTheta=stParams2.nNbTheta+5;
% Also add more theta to post-processing to test accuracy of surface averages
stParams2.nNbThetaPst=stParams2.nNbThetaPst+5;
fprintf('Convergence testing for lambda = %.f\n', lambda(indWorst));
tic;
[stC2, stAbcdnm2, stEsurf2] = slvForNearField(stParams2,stOptions);
relerrExt = (abs(stC.Cext(indWorst)./stC2.Cext-1));
relerrSca = (abs(stC.Csca(indWorst)./stC2.Csca-1));
relerrAbs = (abs(stC.Cabs(indWorst)./stC2.Cabs-1));
relerrExtoa = (abs(stC.Cextoa(indWorst)./stC2.Cextoa-1));
relerrScaoa = (abs(stC.Cscaoa(indWorst)./stC2.Cscaoa-1));
relerrAbsoa = (abs(stC.Cabsoa(indWorst)./stC2.Cabsoa-1));
relerrM = (abs(stEsurf.MLocAve(indWorst)./stEsurf2.MLocAve-1));
relerrMperp = (abs(stEsurf.MLocPerpAve(indWorst)./stEsurf2.MLocPerpAve-1));
relerrF = (abs(stEsurf.F0E4Ave(indWorst)./stEsurf2.F0E4Ave-1));
fprintf('\nT-matrix (N = %d) ... done in %.g seconds.\n', N, toc);
%% Plotting the results
fh = figure('Name','ScriptSolveForNearFieldSpectrum');
set(fh, 'Position', [100, 100, 1000, 500]);
subplot(1,2,1)
plot(lambda,[stC.Cext,stC.Csca,stC.Cabs,stC.Cextoa,stC.Cscaoa,stC.Cabsoa]);
legend({['Cext (err. ', num2str(relerrExt,3),')'], ...
['Csca (err. ', num2str(relerrSca,3),')'], ...
['Cabs (err. ', num2str(relerrAbs,3),')'], ...
[' (err. ', num2str(relerrExtoa,3),')'], ...
[' (err. ', num2str(relerrScaoa,3),')'], ...
[' (err. ', num2str(relerrAbsoa,3),')']}, ...
'Location','Best');
title(['a=', num2str(a), ', c=',num2str(c),', N=', int2str(N), ', Nt=', int2str(nNbTheta)]);
xlabel('Wavelength [nm]')
ylabel('Cross-section [nm^2]')
subplot(1,2,2)
semilogy(lambda,[stEsurf.MLocAve,stEsurf.MLocPerpAve,stEsurf.F0E4Ave]);
legend({['<|E|^2> (err. ', num2str(relerrM,3),')'], ...
['<|E_{perp}|^2> (err. ', num2str(relerrMperp,3),')'], ...
['<|E|^4> (err. ', num2str(relerrF,3),')']}, ...
'Location','Best');
title('Surface-averaged surface field properties');
xlabel('Wavelength [nm]')
ylabel('Field Enhancement Factor')
Execution results
octave>ScriptSolveForNearFieldSpectrum
Loop over 101 lambda values...
Calculating surface-fields for all wavelengths...
... done in 14.686 seconds.
T/R-matrices and near fields (N = 20) ... done in 1e+001 seconds.
Convergence testing for lambda = 900
T-matrix (N = 20) ... done in 1 seconds.
warning: legend: 'best' not yet implemented for location specifier
warning: legend: 'best' not yet implemented for location specifier
Generated graphics
