This module consists of a mix of high-level, core, low-level and
supplementary routines for solving a multiple scattering problem using
the T-matrix formalism. We list below the subroutines of the
multiscat
module with a brief explanation. A list of common
arguments and their brief description is at the end of this section. The
other arguments are explained after each subroutine.
mapNF(ncut, wavelen, inc,ehost, geometry, scheme, tfiles_, escat_, nselect_, verb_, noRTR_, dump_oaE2, dump_oaB2, field, Bfield, N_OC, orAvextEB_int, oa_ldoc, p_label)
Calculates the electric and magnetic near fields, and normalised optical
chirality ($\overline{\mathscr{C}}$)
for a multiple scattering problem, for different incidence directions
and wavelengths, as well as the orientation-averaged value of external
⟨|E|2⟩, ⟨|B|2⟩
and $\langle\overline{\mathscr{C}}\rangle$.
escat_, tfiles_, nselect_, verb_, noRTR_
are optional
inputs. dump_oaE2, dump_oaB2
are logical flags selecting
whether the orientation-averaged values ⟨|E|2⟩ and ⟨|B|2⟩ will be
calculated, respectively.
spectrumFF(ncut, wavelen, ehost, geometry, scheme, escat_, tfiles_, nselect_, noRTR_, verb_, sig_oa_, sig_, sig_abs_, jsig_abs_oa)
Calculates cross-section spectra for (multiple) fixed orientations,
partial absorptions, and orientation-averaged cross-sections for a
particle cluster. T-matrices for individual scatterers are
either constructed using Mie theory or read from an optional argument
tfiles_
.
escat_, tfiles_, nselect_, verb_, noRTR_
are optional
inputs. jsig_abs_oa
contains the orientation-averaged
absorption cross-section of each particle (valid for homogeneous spheres
only, at present).
solve(wavelen, ehost, geometry, nselect_, scheme_, verb_, noRTR_, TIJ, cJ_, cJint_, csAbs_, ierr_ )
This routine is the crux of terms and
solves a given multiple scattering problem by operating in a specified
scheme. TIJ
is an in/output argument,
cJ_, cJint_, csAbs_
are optional in/output arguments,
nselect_, scheme_, verb_, noRTR_
are optional inputs, and
ierr_
is an optional output. TIJ
(lmax × nscat,
lmax × nscat)
as the input argument stores the T-matrix of nonspherical
particles as the diagonal blocks of the matrix, or dielectric values of
different shells for spherical particles as the diagonal elements of the
matrix. nscat
is the number of scatterers. This subroutine
updates and returns TIJ
for the whole system as the output.
cJ_(nscat x l_{max}, 2, nfi)
as the input argument contains
details of the incident field and as a output argument contains incident
plane wave coefficients in the first column and scattering coefficients
in the second column. nfi
is the number of incident angles.
cJint_(nscat x l_{max}, 4, 2)
: contains the regular and
irregular field coefficients for each concentric region inside spherical
scatterers. csAbs_(nscat,4)
: contains absorption cross
section inside each shell of each spherical scatterer.
stageAmat(scatXYZ, scatMiet, rtr, right_, balance_, verb_, A, Tmats_)
Stages a pre-staged matrix A.
A (l_{max} x nscat, l_{max} x nscat)
: an in/output matrix,
which must contain 1-body T-matrices in the diagonal blocks on
input and is a pre-staged matrix on the output;
right_, balance_, verb_
are optional inputs;
Tmats_(l_{max}, l_{max}, nscat)
: an optional output matrix
which contains the 1-body T-matrix of each particle.
balance_
: a logical input argument which determines whether
balancing is applied or not.
calcTIJStout(scatXYZ, scatMiet, rtr, TIJ)
Calculates the scatterer-centred T-matrix using the recursive
scheme presented in Stout. TIJ
is an in/output argument.
TIJ(l_{max} x nscat, l_{max} x nscat)
as the input argument
stores the T-matrix of nonspherical particles as the diagonal
blocks, or dielectric values of different shells for spherical particles
as diagonal blocks.
calcTIMackowski(scatXYZ, scatMiet, rtr, TIJ)
Calculates the cluster’s T-matrix using Mackowski &
Mishchenko’s formulation. TIJ
is an in/output argument.
TIJ(l_{max} x nscat, l_{max} x nscat)
as the input argument
stores the T-matrix of nonspherical particles as the diagonal
blocks, or dielectric values of different shells for spherical particles
as diagonal blocks. The output TIJ
is the scatterer-centred
T-matrices calculated using Mackowski & Mishchenko’s
scheme.
balanceMatJI(j, jregt, iregt, i, rev_, mnq_, Mat)
Performs balancing on a matrix (Mat
) using two weights
(indexed by j and i). Mat
is here taken
as relating two vectors of VSWF coefficients, cj (centred at
j) and ci (centred at
i), such that cj = Mat ci.
Logical inputs jregt
and iregt
specify whether
cj and
ci are
regular or not. Mat
is an in/output argument.
balanceVecJ(j, jregt, rev_, Vec)
Performs balancing on a single vector (V
) with a weight
indexed by j. V
corresponds here to the VSWF coefficients of particle j. Vec
is an
unbalanced/ balanced vector as the in/output argument. j
specifies the scatterer.
calcCsStout(scatXYZR, aJ, fJ, nmax2_, tol_, verb_, sig)
Calculates the extinction, scattering and absorption cross-sections from
the incident and scattered coefficients using the Stout formulae.
nmax2_, tol_, verb_
are optional inputs and
sig
is an in/output matrix.
calcCs(scatXYZR, inc, fJ, nmax2_, tol_, verb_, sig)
Calculates the extinction, scattering and absorption cross-sections from
the incident and scattered coefficients which are collapsed to the
common origin. Depending on the dimension of the sig
, each
cross-section is either just a total sum, or resolved into contributions
from the multipole orders. inc
: a vector of incidence
angles.
calcOAprops(Tmat, rtol_, sigOA, verb_)
Calculates orientation-averaged cross-sections and circular dichroism
(CD) by transforming the T-matrix (Tmat
) from
"parity" (M–N) basis to "helicity" (L–R) basis, following Ref. .
rtol_
is an optional input, verb_
is an
optional output, and sigOA
is an in/output matrix
containing orientation-averaged cross-sections and CD in each column for
n = 1, …, nmax.
contractTmat(Tin, scatXYZR, rtr, mack_, Tout, verb_)
Combines the scatterer-centred T-matrices into a common origin;
the output will be the collective T-matrix (Tout
).
verb_
is an optional in/output, mack_
is an
optional logical input to calculate the collective T-matrix
based on Mackowski & Mishchenko’s scheme.
alphaTensor(T, Alpha)
Conversion of l ≤ 3 spherical
multipoles of the T-matrix T into cartesian multipoles (tensor
Alpha
), following formulas for ‘Higher-Order Polarizability
Tensors’ in Mun. This function is called when the keyword
DumpCollectiveTmatrix
is present, and outputs a file
alpha_col.txt
in the working directory.
diagnoseTmat(mode_, verb_, Tmat)
Determines the value of n ≤ nmax
when Tr (ℜ(Tcol)) converges to rtol_G := 10−ncut(3). If mode_
> 0, also tests for the general
symmetry, which applies to all T-matrices. (See equation 5.34
on p. 121 of Mishchenko).
calcOaStout(TIJ, scatXYZ, verb_, sigOA, cdOA_, jAbsOA)
Calculates the orientation-averaged extinction and scattering
cross-sections defined in equations 44 and 47 of Stout. The absorption
cross-section is then deduced as the difference. TIJ
is the
collective T-matrix, sigOA(3)
contains
orientation-averaged extinction and scattering cross-sections, and
cdOA_
is an optional output containing the corresponding
values of
jAbsOA
: contains the orientation-averaged absorption
cross-section for each particle.
applyRotTranzRotOnMat(vtacs, bigdOP, rightOP, mat)
Performs the factorised translation of T-matrices when changing
origin. Instead of a single multiplication of a T-matrix by a
dense matrix containing the general translation-addition coefficients,
this routine executes three multiplications by sparse matrices
representing 1) a rotation, 2) a translation along the z-axis; and 3) an
inverse rotation. This is meant to be more efficient when high multipole
orders are included.
vtacs(2x pmax,2 x pmax)
: axial VTACs with (m, n, q)
indexing, bigdOP(pmax, pmax)
: optional input for rotation,
rightOP
: an optional logical input argument for applying
the product from the right. mat
: a non/translated matrix as
the in/output.
calcField(r, geometry, ipwVec, ipwE0, scaCJ, intCJreg_, intCJirr_, scatK_, verb_, reE, imE, reB, imB, reE_sca,imE_sca, reB_sca,imB_sca, p_label)
Calculates the electric and magnetic near-field values at the determined
grid points.
r
: a matrix containing the coordinates of the grid points;
ipwVec(3)
,ipwE0(3)
: contain the wavevector and
amplitude of the incident field, respectively; scaCJ
: a
vector containing scattering coefficients,
intCJreg_, intCJirr_
: contain the regular and irregular
parts of the incident field coefficients transformed to the origin of
each particle, respectively, scatK_
is the wavenumber in
the host medium, and
reE, imE, reB, imB, reE_sca,imE_sca, reB_sca,imB_sca
:
contain real(re) and imaginary(im) parts of the total electric (E) and
magnetic (B) fields and the scattered field values at the grid
points.
dumpTmat(tmat, filename, lambda, eps_med, tol_, verb_)
Routine for dumping the collective T-matrix (tmat
)
to a file in the format:
s, s’, n, n’, m, m’, T_re, T_im
filename
is an argument of type character corresponding
to the name of the output file; lambda
: the value of
wavelength; eps_med
: the dielectric value of the host
medium.
dumpMatrix(mat, ofile, tolOP, verb_)
Outputs matrix mat
to a desired optional tolerance
(tolOP
). ofile
: the name of the output
file.
offsetTmat(off, miet, rtr, right, bigD_, useD_, balJI_, Tmat)
Offsets the supplied T-matrix Tmat
by
off
, which can be either a square matrix of VTACs or a
(note: complex!) displacement vector kr(3) from which
VTACs will be generated. Regular or irregular VTACs will be generated
depending on whether kr(3) is purely
real or purely imaginary. If the logical input miet
is
true, Tmat will be treated as diagonal. If the logical input
rtr
is true, then offsetting will be based on factorised
translation. If the logical input right
is true, then
offsetting will be done by post-multiplying Tmat from the right.
balJI_
triggers balancing of the VTACs and the
T-matrices individually, before offsetting, but currently works
only without factorised translation.
readTmatFile(filename, unit, wavelen, verb_, Tmat)
Reads a T-matrix from the input file (filename
)
and import it into the matrix Tmat
. unit
: an
integer indexing the name of the T-matrix file.
wavelen
is the value of the wavelength.
parseInc(inc, verb_, inc_dirn, inc_ampl)
Calculates the amplitude and direction vector of the incident plane wave
based on the input Euler angles (α, β, γ).
inc_dirn
and inc_ampl
are vectors containing
the wavevector and amplitude of the incident electric field in cartesian
coordinates, respectively. inc
is a vector consisting of
polarisation type and Euler angles of the incidence direction.
calcStokesScaVec(sca_angles, inc2, ncut, wavelen, ehost, geometry, scheme, tfiles_, escat_, nselect_, noRTR_, verb_, StokesPhaseMat, StokesScaVec, diff_sca)
Calculates the Stokes phase matrix (StokesPhaseMat
), Stokes
scattering vector (StokesScaVec
), and differential
scattering cross-sections (diff_sca
).
sca_angles
is a matrix of desired scattering angles; if it
is not specified in the input file, they are taken equal to the
incidence angles. inc2
is a matrix containing incidence
angles.
calcLDOC(Ef, Bf, verb_, N_OpC)
Calculates the normalised optical chirality ($\overline{\mathscr{C}}$) relative to the
optical chirality of circularly polarised light. Ef
,
Bf
, and N_OpC
are matrices containing the
electric and magnetic field, and $\overline{\mathscr{C}}$ values at the grid
points, respectively.
calcOaNFUnpol(r, geometry, TIJ, lambda, ehost, escat, p_label, verb_, orEB2)
Calculates the orientation average of the total external electric and
magnetic field intensities. r
is a matrix containing the
cartesian coordinates of the grid points. TIJ
is the
scatterer-centred T-matrix of the cluster. orEB2
is a vector containing the value of orientation-averaged external
electric and magnetic field intensities at the grid points.
calcOaNF(pol_type, r, geometry, TIJ, Oa_OC, ehost, p_label, scatK_, verb_, Oa_EB2)
Calculates the orientation average of normalised optical chirality $\langle\overline{\mathscr{C}}\rangle$, and
near field intensities $\langle
|\mathbf{E}_\mathrm{tot}(\it{k}\mathbf{r})|^2\rangle$, $\langle
|\mathbf{B}_\mathrm{tot}(\it{k}\mathbf{r})|^2\rangle$ for
circularly polarised incident light. pol_type
is the
polarisation type, r
is a matrix containing the cartesian
coordinate of the grid points. TIJ
is the scatterer-centred
T-matrix of the structure.
calcTrace(TRANSA, TRANSB, A, B, tr)
Calculates the trace of a product of two matrices,
op(A)*op(B)
. The input characters TRANSA
and
TRANSB
determine the operation op
, following
the convention of blas’
gemm
. Specifically, op = ’N’
corresponds to
o
p
(
A
)
=
A
(no operation), whereas op = ’C’
corresponds to o
p
(
A
)
=
A
†
.
acs_int_
: a matrix containing partial internal
absorption inside each scatterer and for each shell.
aJ(\text{nscat}\times l_{max}), fJ(\text{nscat}\times l_{max})
:
contains incident and scattering coefficients.
Bfield
: contains the real and imaginary parts of the
magnetic near field at the specified grid points, wavelengths, and
incidence.
ehost
: a vector of dielectric permittivity of the
host medium at specified wavelengths.
escat_(nscat, 4, size(wavelen))
: depending on the
number of wavelengths, it is a 2D or 3D array of dielectric values for
each scatterer, for each shell and wavelength.
field
: contains the real and imaginary parts of the
electric near field at the specified grid points, wavelengths, and
incidence.
geometry
: a matrix containing physical information
of different scatterers such as centre, dimensions and
direction.
ierr_
: an integer value (0 or 1 or 2); 0 indicates
solving was successful, 1 means there is an error in processing
arguments, and 2 means an error in prestaging, staging, or
solving/inverting Ax = b.
iregt
: logical input, specifies whether vectors are
regular or not.
jregt
: logical input, specifies whether vectors are
regular or not.
mnq_
: an optional logical argument which is false by
default, but if true will change the indexing convention from (q,n,m) to
(m,n,q), which is used to make the z-axial VTACs block-diagonal. Note
that index q corresponds to
s in this user guide.
ncut
: a vector in the form [n1, n2, tol], which contains
the values corresponding to the keyword "MultipoleCutoff
".
Default values: [8, 8, −8].
nmax2_
: an integer value equals to ncut(2).
noRTR_
: an optional input with logical value
.true.
or .false.
for the keyword
DisableRTR
. Default: .false.
.
nselect_
: an optional input matrix which includes
information about multipole selection for different scatterers.
oa_ldoc
(npts × 4 × nwavelen): contains the
orientation averaged value of $\overline{\mathscr{C}}$ at different grid
points and wavelengths.
orAvextE_int
(npts × nwavelen): contains the orientation
averaged electric field intensity values at different grid points and
wavelengths.
p_label
: a matrix determining the position of each
grid point, whether it is inside the surrounding medium or particles,
and in which layer.
rev_
: an optional logical input which is false by
default; triggers the reverse of balancing – "unbalancing".
right_
: a logical input. There are two ways for
obtaining the TIJ
matrix. This argument determines whether
the product is taken from the left or from the right.
rtr
: a logical input that is the reverse of
noRTR_
.
scatMiet
(nscat): a logical vector with
.true.
and .false.
values, determining whether
a scatterer is spherical or not.
scatXYZ
(3,nscat): a matrix containing the cartesian
coordinates (in lab frame) of the particle’s centre.
scatXYZR
(4,nscat): a matrix containing the cartesian
coordinates (in lab frame) and the radius of the smallest circumscribed
sphere of each particle.
scheme, scheme_
: an integer value specifying the
selected scheme.
sig_
: a matrix containing cross-sections
(Extinction, Scattering, Absorption) for different polarisation(s),
wavelength(s), and incidence(s).
sig_abs_
(4 × nscat × 4 × nwavelen × nfi): a 5D array
containing absorption cross-sections inside each shell for each
scatterer for 4 Jones vectors, different wavelengths and different
incidence directions.
sig_oa_
(6 × n × nwavelen): a matrix consisting of
orientation-averaged cross-sections and CD at different wavelengtha. The
first column gives the values for nmax
and other columns contain values for different value of n = 1, …, ncut(2).
tfiles_
: a matrix of character type, includes the
T-matrix filename and filepath for non-spherical
scatterers.
tol_, rtol_
: a real value
rtol_G = 10^{\textbf{ncut(3)}}
.
N_OC
: contains $\overline{\mathscr{C}}$ at the specified
grid points, wavelengths, and incidence.
verb_
: an integer variable containing the verbosity
value ( ∈ [0, 1, 2, 3]) (the default
value is 1).
wavelen
: a vector of specified
wavelength(s).