Source code
%% ScriptTutorialSpectrum
% Example script providing a step-by-step tutorial to solve the scattering
% problem as a function of wavelength, with explicit calls to the low-level
% functions used in the calculations of intermediate quantities.
%
% As such, this script provides a more in-depth understanding of the code
% than provided by the other example scripts. For a more application-oriented
% perspective, the "ScriptSolve" family of scripts use pre-defined functions
% performing the same steps internally, but invisible to the end-user.
%%
%% Description
%
% The example considers a silver prolate spheroid of semi-axes 20 x 80 nm
% (aspect ratio 4), in water.
% The script calculate the T-matrix up to multipole order N and for all 0<=m>
%
a=10;
c=40;
h = c/a; % aspect ratio, h=c/a for prolate spheroids
lambda = (400:2:900)'; % [L x 1] in nm
epsilon2 = epsAg(lambda); % [L x 1] Dielectric function of particle
epsilon1 = 1.33^2; % scalar or % [L x 1] Dielectric function of medium (water here)
% Incident field properties (if required)
sIncType='KxEz'; % incident along x and polarized along z
% type doc vshMakeIncidentParameters for more options (plane wave excitation only)
% note that only k1 and s are needed for the T-matrix calculation
% but other parameters are needed for E-field calculations
k1 = 2*pi*sqrt(epsilon1)./lambda; % [L x 1]
s = sqrt(epsilon2)./sqrt(epsilon1); % [L x 1] relative refractive index of particle
xmax = max(a,c)*k1; % maximumum size parameter xmax= k1 * max(a,c)
% For convenience, k1 and s are stored in a struct
stParamsAll.k1=k1;
stParamsAll.s=s;
%% Parameters governing convergence
% The following parameters will be needed:
%%
% * N: Number of multipoles required for T-matrix
% * abmvec: Vector containing the values of |m| for which T is computed
% Only m>=0 is needed and m<=N.
% For all m, use absmvec=0:N (= [0,1,2 ..., N] )
% * NQ: Number of multipoles for the P and Q matrix, NQ>=N
% * NB: Number of multipoles to compute the Bessel functions
% in the improved algorithm, NB>=NQ.
% * nNbTheta: Number of theta's for the Gaussian quadrature
% In this example script, these are chosen constant (independent of lambda)
% Maximum multipole order for T-matrix and series expansions of fields
N = 20;
% m-numbers used in the calculations
% For most incident excitations and for orientation-averaged properties,
% all |m|<=N need to be considered:
absmvec = (0:1:N)';
% Advanced users can define the stIncPar first and use the following instead
% absmvec = stIncPar.absmvec.';
% or specify a single m-value for testing for example, i.e.
% absmvec = 1; % m=1 only
% Number of points for Gaussian quadratures to compute integrals in P and Q matrices
% By symmetry, points are only computed from theta=0 to pi/2
nNbTheta = 50;
% Make structure describing spheroidal geometry and quadrature points for
% numerical integrations
stGeometry = sphMakeGeometry(nNbTheta, a, c);
% Make structure with incident field parameters
stIncPar = vshMakeIncidentParameters(sIncType, N);
%% Defining number of multipoles for each step
%
% The T-matrix and corresponding field expansion coefficients will be
% calculated up to n=N for all m in absmvec (note that m>=0).
% Here, we use simply use NQ = N and check convergence of the
% results by running a second calculation for a larger N
% [see JQSRT 160, 29 (2015)].
NQ = N; % Maximum multipole order for computing P and Q matrices
% P and Q are calculated using the stable and accurate algorithm in
% [JQSRT 123, 153 (2013)].
% For this algorithm we need to specify how many extra order are needed
% to compute Bessel function products to the highest accuracy
% This can be estimated with the following function
NB=sphEstimateNB(NQ, stGeometry, stParamsAll);
% or can be specified by user (advanced), for example
% NB=NQ;
fprintf('\n... done in %.g seconds\n', toc);
%% Calculation of T-matrix and optical properties
fprintf('\nLoop over lambda...\n');
tic;
% Coefficients of incident wave do not change with wavelength
stIncEabnm=vshGetIncidentCoefficients(N,stIncPar);
% Initialize variables to store results for each wavelength
stAbcdnm.anm = stIncEabnm.anm; % [1 x P] where P=N(N+2)
stAbcdnm.bnm = stIncEabnm.bnm; % [1 x P] where P=N(N+2)
P=N*(N+2); % number of elements in p-index
L=length(lambda);
stAbcdnm.pnm = zeros(L,P);
stAbcdnm.qnm = zeros(L,P);
stAbcdnm.cnm = zeros(L,P);
stAbcdnm.dnm = zeros(L,P);
stCoa = struct();
stCoa.Csca = zeros(L,1);
stCoa.Cext = zeros(L,1);
stCoa.Cabs = zeros(L,1);
%% loop over lambda to calculate the wavelength-dependent T-matrix
% The T/R matrices are not stored to avoid memory issues, but all the
% physical properties and field expansion coefficients can be calculated in
% the loop and stored
bGetR = true; % To calculate R and internal field
stParams1 = struct(); % temp structure for simulation
stParams1.bOutput = false; % less verbose output
fprintf('lambda = ');
for lInd=1:length(lambda)
% uncomment next line for real-time progress
% fprintf('\b\b\b%.3g', lambda(lInd));
stParams1.k1 = stParamsAll.k1(lInd);
stParams1.s = stParamsAll.s(lInd);
% This calculates P and Q using the algorithm of [JQSRT 123, 153 (2013)]
CstPQa = sphCalculatePQ(NQ, absmvec, stGeometry, stParams1, NB);
% Get T=-PQ^{-1} and R=Q^{-1} for all m using the inversion procedures
% described in [JQSRT 123, 153 (2013)].
CstTRa = rvhGetTRfromPQ(CstPQa,bGetR);
% If only the T-matrix is required, use instead bGetR=false
% If needed, discard higher order multipoles
% (which are affected by the finite size of P and Q)
if NQ>N
CstTRa = rvhTruncateMatrices(CstTRa, N);
end
% T and R matrices now include N multipoles
% If required, one may symmetrize the T-matrix (this assumes that the upper
% triangular parts of the matrices are correct, see JQSRT 160, 29 (2015))
% CstTRa = rvhGetSymmetricMat(CstTRa, {'st4MT'});
% Calculate the (Ext, Abs, Sca) orientation-averaged cross-sections
stQoaOneLambda = rvhGetAverageCrossSections(stParams1.k1, CstTRa);
stCoa.Csca(lInd) = stQoaOneLambda.Csca;
stCoa.Cext(lInd) = stQoaOneLambda.Cext;
stCoa.Cabs(lInd) = stQoaOneLambda.Cabs;
% Get the field expansion coefficients from T and R for a given incident
% excitation (defined earlier in stIncPar)
stAbcdnmOneLambda = rvhGetFieldCoefficients(N, CstTRa, stIncPar,stIncEabnm);
stAbcdnm.pnm(lInd,:) = stAbcdnmOneLambda.pnm;
stAbcdnm.qnm(lInd,:) = stAbcdnmOneLambda.qnm;
if bGetR % internal fields only if R has been calculated
stAbcdnm.cnm(lInd,:) = stAbcdnmOneLambda.cnm;
stAbcdnm.dnm(lInd,:) = stAbcdnmOneLambda.dnm;
end
end
fprintf('\n\nT-matrices (N = %d) ... done in %.g seconds.\n', N, toc);
%% Further post-processing for well-defined orientation
%
fprintf('\nPost-processing and plotting ...\n', N, toc);
% Calculate the (Ext, Abs, Sca) cross-sections for fixed orientation
% This can be done for all wavelengths in one go
stQ = pstGetCrossSections(k1, stAbcdnm);
% For surface-properties, we need to re-define the geometry with theta over [0;pi]
% and with more integration points, to ensure finer evaluation of the
% surface fields and of their averages
nNbThetaPst = 360; % number of theta for evaluating fields
stRtfuncPst = sphMakeGeometry(nNbThetaPst, a, c); % new geometry with more points
% It is also necessary to extend the range of theta over [0;pi] instead of
% [0;pi/2]
stRtfuncPst=rvhGetThetasForAverage(stRtfuncPst); % get thetas over entire range [0,pi]
stGeometryPts=sphMakeGeometry(0, a, c, [0; pi/2]); % This is to evaluate the field at theta=0 and pi/2
% For convenience this will prepare a result structure for postprocessing
stResE=pstMakeStructForField(stAbcdnm, N, lambda, epsilon2, epsilon1, stIncPar,a,c);
% Calculates the surface electric field E partial series expansion (for each m)
% on the surface as well as average values M=|E|^2, F=|E|^4
% This is done for all wavelengths
stEsurf=pstSurfaceField(stResE,stRtfuncPst);
%% Example of plots at fixed wavelength
lambda0 = 754; % Position of resonance in extinction
% The following plots the surface field at one or more given
% phi (as a function of theta) for a fixed lambda0
phi0=[0,pi/4,pi/2];
M = pstGetThetaDepFieldIntensity(stEsurf,phi0,lambda0); % [3 x T]
figure('Name',['Theta-dependence of surface-field intensity M=|E|^2 for fixed phi at lambda=',num2str(lambda0)]);
semilogy(stEsurf.theta, M);
legend({['phi=', num2str(phi0(1))], ...
['phi=', num2str(phi0(2))], ...
['phi=', num2str(phi0(3))]}, ...
'Location', 'Best');
xlabel('Theta [radians]');
ylabel('Field Intensity Enhancement Factor');
% The following makes a 3D surface plot of the surface
% field everywhere on the surface (note that this requires to recompute
% the surface fields). Use 90x90 pts here:
pstPlotAllSurfaceField(90,stResE,lambda0);
%% Example of plots as a function of wavelength
% Calculates field intensity at theta=0, pi/2, and phi=0
stEsurfPts=pstSurfaceField(stResE,stGeometryPts);
stEphi=vshEthetaForPhi(stEsurfPts,0);
Mpts=abs(stEphi.Er).^2+abs(stEphi.Et).^2+abs(stEphi.Ef).^2;
fh = figure('Name','Wavelength-dependent properties');
set(fh, 'Position', [100, 100, 1000, 800]);
subplot(2,2,1)
plot(lambda, [stCoa.Cext, stCoa.Cabs, stCoa.Csca]);
legend({'Extinction-OA', 'Absorption-OA', 'Scattering-OA'},'Location','Best');
xlabel('Wavelength [nm]');
ylabel('Orientation-averaged Cross-section [nm^2]');
title('Far-field orientation averaged cross-sections');
subplot(2,2,2)
plot(lambda, [stQ.Cext, stQ.Cabs, stQ.Csca]);
legend({'Extinction', 'Absorption', 'Scattering'},'Location','Best');
xlabel('Wavelength [nm]');
ylabel('Cross-section [nm^2]');
title(['Far-field cross-sections for fixed orientation: ', stIncPar.type]);
subplot(2,2,3)
semilogy(lambda, [stEsurf.MLocAve,stEsurf.MLocPerpAve,Mpts]);
legend({'', '', 'M at \theta=0', 'M at \theta=\pi/2,\phi=0'},'Location','Best');
xlabel('Wavelength [nm]');
ylabel('Field-Intensity Enhancement Factors');
Execution results
octave>ScriptTutorialSpectrum
Initialization and finding NQ and NB...
... done in 0.2 seconds
Loop over lambda...
lambda =
T-matrices (N = 20) ... done in 8e+001 seconds.
Post-processing and plotting ...
warning: legend: 'best' not yet implemented for location specifier
warning: legend: 'best' not yet implemented for location specifier
warning: legend: 'best' not yet implemented for location specifier
warning: legend: 'best' not yet implemented for location specifier
Generated graphics


