ScriptSolveForFixedSpectrum

An example script showing how to obtain the field expansion coefficients and far-field cross-sections for a spheroid in a fixed orientation, as a function of wavelength. Plots the wavelength-dependent spectra for extinction, scattering, and absorption cross-sections for fixed orientation, as well as orientation-averaged.

Contents

Instructions

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:

clear all
close all

Parameters of the scattering problem

We define parameters for a gold nanorod in water, modeled as a prolate spheroid

a = 15; % in nm
c = 45; % in nm, i.e. 30 x 90nm full-axes
lambda = (400:5:1100).'; % in nm
epsilon2 = epsAu(lambda);
epsilon1 = 1.33^2; % for water
% Define incident field 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;

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.epsilon2=epsilon2;
stParams.epsilon1=epsilon1;
stParams.sIncType = sIncType;

% Optional parameters may also be defined as follows:
stOptions.bGetR = false;
stOptions.Delta = 0;
stOptions.NB = 0; % NB will be estimated automatically
stOptions.bGetSymmetricT = false;
stOptions.bOutput = false; % suppress messages in lambda-loop

T-matrix calculation

Solve for T (all wavelengths)

tic;
stC = slvForFixedSpectrum(stParams,stOptions);
fprintf('\nT-matrix (N = %d) ... done in %.2f seconds.\n', N, toc);
% use this call instead to keep the expansion coefficients for further
% processing
%[stC, stAbcdnm] = slvForT(stParams,stOptions);
Loop over 141 lambda values...

T-matrix (N = 20) ... done in 22.92 seconds.

Convergence test

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;
fprintf('\nConvergence testing for lambda = %.2f.\n', lambda(indWorst));
tic;
stC2 = slvForFixed(stParams2,stOptions);
fprintf('\nT-matrix (N = %d) ... done in %.2f seconds.\n\n', N, toc);

% relative errors
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));
Convergence testing for lambda = 1100.00.

T-matrix (N = 20) ... done in 0.37 seconds.

Plotting the results

figure('Name','ScriptSolveForFixedSpectrum');
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),')'], ...
    ['<Cext> (err. ', num2str(relerrExtoa,3),')'], ...
    ['<Csca> (err. ', num2str(relerrScaoa,3),')'], ...
    ['<Cabs> (err. ', num2str(relerrAbsoa,3),')']}, ...
    'Location','Best');
title(['a=', num2str(a), ', c=',num2str(c),', N=', int2str(N), ', Nt=', int2str(nNbTheta)]);
xlabel('Wavelength [nm]') % x-axis label
ylabel('Cross-section [nm^2]') % y-axis label