ScriptTutorial

This script provides a step-by-step tutorial to solve the scattering problem at a single 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.

Contents

Description

The example considers an absorbing dielectric (s=1.5+0.02i) in the shape of prolate spheroid of aspect ratio 10, and size parameter 10. The script calculate the T-matrix up to multipole order N and for all 0<=m<N From it various optical properties are also derived:

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');
clear all
close all

fprintf('Initialization and finding NQ and NB...\n');
tic; % for timer
Initialization and finding NQ and NB...

Scatterer properties

The following parameters should be defined:

In addition, the following parameter will be needed if field expansion coefficients and/or the cross-sections for a given fixed orientation are sought:

For actual field calculations (such as surface fields), we will also need

Parameters of the scattering problem

We define parameters for a prolate spheroid

h = 10; % aspect ratio, h=c/a for prolate spheroids
k1 = 1; % incident wavenumber k1=2pi/lambda * nM
xmax = 10; % maximumum size parameter xmax= k1 * max(a,c)
% From which we deduce
c = xmax / k1;
a = c / h;

epsilon2 = (1.5+0.02i)^2; % dielectric constant of scatterer
epsilon1 = 1;
lambda = 2*pi*sqrt(epsilon1)/k1;
s = sqrt(epsilon2)/sqrt(epsilon1);
% note that only s is needed for T-matrix but other parameters are needed for E-field calculations

% For convenience, k1 and s are stored in a struct
stk1s.k1=k1;
stk1s.s=s;
% Incident field properties (if required)
sIncType='KxEz'; % incident along x  and polarized along z
% type  help vshMakeIncidentParameters for more options (plane wave excitation only)

Parameters governing convergence

The following parameters will be needed:

% Maximum multipole order for T-matrix and series expansions of fields
N = 50;
% 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 that=0 to pi/2
nNbTheta = 120;

% 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).

% For this, we may need to calculate P and Q with more multipole,
% NQ=N+Delta [see JQSRT 160, 29 (2015)].
% This can ensure the accuracy of the entire T-matrix [see JQSRT 160, 29 (2015)].
% Delta can be estimated from the convergence of T_{11}^{22,m=1}
% with the following call (advanced users)
[Delta, T2211err]= sphEstimateDelta(stGeometry, stk1s);

% Alternatively one may simply use NQ = N and check convergence of the
% results by running a second calculation for a larger N
% [see JQSRT 160, 29 (2015)].
% Delta=0;

NQ = N+Delta;% 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, stk1s);
% or can be specified by user (advanced), for example
% NB=NQ;

fprintf('\n ... done in %.g seconds.\n', toc);
sphCalculatePQ: Calculating P,Q for 1 m-values with N_Q = 80, N_B = 80, N_Theta = 120

 ... done in 3 seconds.

Calculation of the T-matrix

Calculating P and Q

fprintf('\nMain business...\n');
tic;
% This calculates P and Q using the algorithm of [JQSRT 123, 153 (2013)]
CstPQa = sphCalculatePQ(NQ, absmvec, stGeometry, stk1s, NB);
% CstPQa is a cell {1 x M} of structures, one for each m in absmvec. Each stPQa structure
% contains structures describing P (st4MPeo, st4MPoe) and Q (st4MQeo,
% st4MQoe). These make use of the reflection symmetry and do not contain
% the zero elements.

% To obtain the full matrices in standard form, use the following call
% jj=2; % jj=1 corresponds to m=0, jj=2 to m=1, etc... (see absmvec definition)
jj = 1;
[Pm, nvecP] = rvhGetFullMatrix(CstPQa{jj},'st4MP');
[Qm, nvecQ] = rvhGetFullMatrix(CstPQa{jj},'st4MQ');
fprintf('Pm and Qm are full matrices for m=%d\n', absmvec(jj));
% nvec contains the multipole orders for each block of Pm and Qm
% Note that P and Q are square matrices of dimension NQ+1-m (or NQ for m=0)
% since only elements P_{nk} with n,k>=m (or 1 for m=0) are needed.
% Therefore for a given m, P(i,j) correspond to P_{n=i+m-1,k=j+m-1}
% or P_{n=i+m,k=j+m} for m=0

fprintf('P and Q matrices... done in %.g seconds.\n', toc);

% Calculating T and R
%
tic;
% Get T=-PQ^{-1} and R=Q^{-1} for all m using the inversion procedures
% described in [JQSRT 123, 153 (2013)].
CstTRa = rvhGetTRfromPQ(CstPQa,true);
% CstTRa has the same format as CstPQa and contains T and R matrices for
% all m in absmvec
% These T and R matrices go up to NQ multipoles

% If only the T-matrix is required, use instead
% CstTRa = rvhGetTRfromPQ(CstPQa,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'});

fprintf('\nT-matrix (N = %d) ... done in %.g seconds.\n', N, toc);

% Calculate the (Ext, Abs, Sca) orientation-averaged cross-sections
stCoa = rvhGetAverageCrossSections(k1, CstTRa);

fprintf('\nCross sections for orientation-averaged excitation:\n');
fprintf('<Cext> = %.20g\n', stCoa.Cext);
fprintf('<Csca> = %.20g\n', stCoa.Csca);
fprintf('<Cabs> = %.20g\n', stCoa.Cabs);
Main business...
sphCalculatePQ: Calculating P,Q for 51 m-values with N_Q = 83, N_B = 83, N_Theta = 120
Pm and Qm are full matrices for m=0
P and Q matrices... done in 6 seconds.

T-matrix (N = 50) ... done in 0.3 seconds.

Cross sections for orientation-averaged excitation:
<Cext> = 18.737732166309676529
<Csca> = 16.601204532914721312
<Cabs> = 2.1365276333949552168

Further post-processing for well-defined orientation

fprintf('\nComputing fixed orientation...\n\n');
tic;
% For this section, it is assumed that the incident field properties are
% defined in stIncPar

% Get the field expansion coefficients from T and R for a given incident
% excitation (defined earlier in stIncPar)
stAbcdnm = rvhGetFieldCoefficients(N, CstTRa, stIncPar);

% Calculate the (Ext, Abs, Sca) cross-sections for this excitation
stC = pstGetCrossSections(k1, stAbcdnm);

format long
fprintf('\nCross sections for fixed excitation (%s):\n', sIncType);
fprintf('Cext = %.20g\n', stC.Cext);
fprintf('Csca = %.20g\n', stC.Csca);
fprintf('Cabs = %.20g\n', stC.Cabs);

% For convenience this will prepare a result structure for postprocessing
stResE=pstMakeStructForField(stAbcdnm, N, lambda, epsilon2, epsilon1, stIncPar, a, c);

% For surface-properties, we may want to use more points to ensure finer
% evaluation of the surface fields and of their averages
nNbThetaPst = 360; % number of theta for evaluating fields

% The surface fields can be caculated as
% stEsurf=pstSurfaceField(stResE,nNbThetaPst);
% To avoid repeated calculation of the geometry, one may also do it
% separately like so:
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]
% Calculate 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
 stEsurf=pstSurfaceField(stResE,stRtfuncPst);
Computing fixed orientation...


Cross sections for fixed excitation (KxEz):
Cext = 25.506059713694384072
Csca = 23.011763940346192214
Cabs = 2.4942957733481918581

Plots of the surface field

as a function of theta, for one or more values of phi

phi0=[0,pi/4,pi/2];
M = pstGetThetaDepFieldIntensity(stEsurf,phi0); % [3 x T]
figure('Name','Theta-dependence of surface-field intensity M=|E|^2 for fixed phi');
plot(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')

% 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:
pstPlotAllSurfaceField(90,stResE);

fprintf('\nPostprocessing for fixed orientation (Ntheta = %d) over [0;pi] ... done in %.g seconds.\n', ...
        2*nNbThetaPst, toc);
Postprocessing for fixed orientation (Ntheta = 720) over [0;pi] ... done in 3 seconds.