Source code
%% ScriptSolveForT
% Example script showing how to obtain the T-matrix and
% the scattering matrix for random orientation for a spheroid
% at a single wavelength. Prints the cross-sections with accuracy estimates,
% saves the T-matrix elements to an external text file, and plots the
% theta-dependent scattering matrix elements.
%%
%% 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
clear
close all
%% Parameters of the scattering problem
% We define aspect ratio, wavenumber, and size parameter for a
% prolate spheroid
%
% <<../fig/schematicp.png>>
%
h = 10; % aspect ratio, h=c/a for prolate spheroids
s = 1.5+0.02i; % relative refractive index
k1 = 1; % incident wavenumber k1=2pi/lambda * nM
xmax = 10; % maximum size parameter xmax= k1 * max(a,c)
% ... from which we deduce
c = xmax / k1;
a = c / h;
%% Convergence parameters
% Maximum multipole order for T-matrix and series expansions of fields
N = 30;
% Number of points for Gaussian quadratures to compute integrals in P and Q matrices
nNbTheta = 120;
%% Collect simulation parameters in a structure
stParams.a=a; stParams.c=c;
stParams.k1=k1; stParams.s=s;
stParams.N=N; stParams.nNbTheta=nNbTheta;
% 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;
%% Solving for the T-matrix
tic;
% For the Scattering matrix, we need to keep the entire T-matrix
[stCoa, CstTRa] = slvForT(stParams, stOptions);
fprintf('\nT-matrix (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...\n');
tic;
stParams2=stParams;
stParams2.N=stParams2.N+5;
stParams2.nNbTheta=stParams2.nNbTheta+5;
[stCoa2, CstTRa2] = slvForT(stParams2,stOptions);
fprintf('\nT-matrix (N = %d) ... done in %.f seconds.\n', N, toc);
%% Reshape the T-matrix to long format, and export to text file
T = exportTmatrix(CstTRa, 'Tmatrix.txt');
%% Display orientation-averaged results
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 orientation-averaged excitation (and estimated accuracy):\n');
fprintf(' = %.20g, relative error: %.2g\n', stCoa.Cext, ...
abs(stCoa.Cext./stCoa2.Cext-1));
fprintf(' = %.20g, relative error: %.2g\n', stCoa.Csca, ...
abs(stCoa.Csca./stCoa2.Csca-1));
fprintf(' = %.20g, relative error: %.2g\n', stCoa.Cabs, ...
abs(stCoa.Cabs./stCoa2.Cabs-1));
%% Calculate the scattering matrix, test its accuracy, and plot the results
lambda=2*pi/k1; % We need lambda here, so assume embedding medium is air
nNbThetaSM=360;
tic
stSM = pstScatteringMatrixOA(CstTRa,lambda,stCoa.Csca,nNbThetaSM);
fprintf('Scattering Matrix (N = %d) ... done in %.f seconds.\n', N, toc);
tic
stSM2 = pstScatteringMatrixOA(CstTRa2,lambda,stCoa2.Csca,nNbThetaSM);
fprintf('Scattering Matrix (N = %d) ... done in %.f seconds.\n', N+5, toc);
% errors in expansion coefficients alpha and beta
errRelAB = abs(stSM.AllAB./stSM2.AllAB(1:(2*N+1),:)-1);
errAbsAB = abs(stSM.AllAB-stSM2.AllAB(1:(2*N+1),:));
% errors in angle-dependent scattering matrix elements
errRelF = abs(stSM.AllF./stSM2.AllF-1);
errAbsF = abs(stSM.AllF-stSM2.AllF);
figure('Name','Scattering matrix')
plot(stSM.thetadeg,stSM.AllF(:,2:7));
xlabel('Theta [deg]')
ylabel('Scattering Matrix Element')
legend({'F_{11}','F{22}','F{33}','F{44}','F{12}','F{34}'})
title(['Maximum absolute error in F is ', num2str(max(max(errAbsF)))]);
Execution results
octave>ScriptSolveForT
sphCalculatePQ: Calculating P,Q for 31 m-values with N_Q = 30, N_B = 30, N_Theta = 120
T-matrix (N = 30) ... done in 1 seconds.
Convergence testing...
sphCalculatePQ: Calculating P,Q for 36 m-values with N_Q = 35, N_B = 35, N_Theta = 125
T-matrix (N = 30) ... done in 1 seconds.
Results for a=1, c=10, k1=1, s=1.5+0.02i, N=30, Nt=120
Cross sections for orientation-averaged excitation (and estimated accuracy):
= 18.737732166323013416, relative error: 7e-013
= 16.601204532910415423, relative error: 2.8e-013
= 2.1365276334125979929, relative error: 8.3e-012
Scattering Matrix (N = 30) ... done in 317 seconds.
Scattering Matrix (N = 35) ... done in 552 seconds.
Generated graphics
