Source code

%% ScriptTutorial % Example script providing 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. %% %% 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> % 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: %% % * 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 % 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); %% 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(' = %.20g\n', stCoa.Cext); fprintf(' = %.20g\n', stCoa.Csca); fprintf(' = %.20g\n', stCoa.Cabs); %% 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); %% 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); %% Example of near-field calculation away from the surface tic disp('Near-field calculation away from the surface...'); % Calculates field at points along the $x$-axis % Note that the method does not work if one gets too close to the surface % (or it needs A LOT of integration points) thetaNF=pi/2; phiNF=0; nNbPtsNF=30; % number of points rNF=linspace(1.05*a,40*a,nNbPtsNF); stPtsNF.r=rNF; % [1 x RP] stPtsNF.theta=thetaNF+zeros(1,nNbPtsNF); % [1 x RP] stPtsNF.phi=phiNF+zeros(1,nNbPtsNF); % [1 x RP] nNbThetaNF = 180; % number of angles (phi and theta) for integral in near field % note that the number of integration points is the square of that. stRtfuncNF = sphMakeGeometry(nNbThetaNF, a, c); % new geometry with points % It is also necessary to extend the range of theta over [0;pi] instead of % [0;pi/2] stRtfuncNF=rvhGetThetasForAverage(stRtfuncNF); % get thetas over entire range [0,pi] % Calculates scattered near-field using integral expression stEscaNF = pstGetNearField(stResE, stRtfuncNF, stPtsNF); ENF=stEscaNF.E; % field amplitude % Get the scattered field from the VSH series expansion with % coefficients pnm and qnm and irregular VSH (with 'h1') stEtheta=vshEgenThetaAllPhi(stResE.lambda,stResE.epsilon1,stResE.pnm,stResE.qnm, ... stPtsNF.r,stPtsNF.theta,'h1'); stEPhi=vshEthetaForPhi(stEtheta,phiNF); % Get result in cartesian coordinates cost=cos(thetaNF); sint=sin(thetaNF); cosf=cos(phiNF); sinf=sin(phiNF); ExyzSeries = [stEPhi.Er .*sint*cosf + stEPhi.Et .*cost*cosf - stEPhi.Ef * sinf; ... stEPhi.Er .*sint*sinf + stEPhi.Et .*cost*sinf + stEPhi.Ef *cosf; ... stEPhi.Er .*cost - stEPhi.Et .*sint]; % [Ex; Ey; Ez] is [3 x RP] ESeries = realsqrt(sum(abs(ExyzSeries).^2,1)); % field amplitude figure; subplot(1,2,1) semilogy(rNF,[ENF;ESeries]); legend('Integral expression','Series (divergent close to the surface)'); title('Scattered field amplitude') xlabel('Distance from surface'); ylim([0.1,10]); subplot(1,2,2) plot(rNF,log10(abs((ENF-ESeries)./ESeries))); title('Relative errors between integral expression and scattered field amplitude') xlabel('Distance from surface'); ylim([-16,1]); disp(['... done in ', num2str(toc), ' seconds']);

Execution results

octave>ScriptTutorial Initialization and finding NQ and NB... sphCalculatePQ: Calculating P,Q for 1 m-values with N_Q = 80, N_B = 80, N_Theta = 120 ... done in 5 seconds. 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.4 seconds. Cross sections for orientation-averaged excitation: = 18.737732166309857718 = 16.601204532914803025 = 2.1365276333950546928 Computing fixed orientation... Cross sections for fixed excitation (KxEz): Cext = 25.506059713694494206 Csca = 23.011763940346138924 Cabs = 2.4942957733483552829 warning: legend: 'best' not yet implemented for location specifier Postprocessing for fixed orientation (Ntheta = 720) over [0;pi] ... done in 2 seconds. Near-field calculation away from the surface... ... done in 7.8365 seconds

Generated graphics