termsProgram
termsProgram
is the main module, with subroutines listed
below; it reads the keywords and their corresponding values in the input
file and then calls different subroutines of multiscat
module for the calculation of requested outputs.
readInputFile(inputfile)
Reads the input file containing specific keywords and the corresponding
parameter values.
errorParsingArguments(keyword)
If there is an error with the parameter values assigned to a keyword,
this subroutine informs the user and stops the program.
calcEpsilon()
Updates the escat
array storing the relative dielectric
function(s) for each scatterer evaluated at the specified
wavelengths.escat
is an array for which the number of rows, columns and
the 3\(^{rd}\) dimension correspond to
the number of shells, scatterers, and wavelengths,
respectively.
calcGridPoints(points)
For a near-field calculation, this subroutine calculates the grid points
based on the number of steps, lower and upper values along the desired
axes.points(3, nGridPoints)
is an in/output matrix storing the
cartesian coordinates \((x,y,z)\) of
the grid points.
sentence2words(sentence, words, nwords_)
Reads each line of the input file as a sentence and splits it into
space-separated words.sentence
, words
are in/output character
arrays, and nwords_
is an optional output integer
containing the number of words in the sentence.
dumpNFs2TXTFile(filename, incidences, Epower, wavelen, work, Ef, p_label)
Exports electric and magnetic near fields into a plain text file.
filename
is the name of the text file.
incidences, Epower, wavelen, work
are the incidence angles,
selected power for mapping fields, wavelength, and near-field
quantities, respectively. Ef
is a logical flag which
selects either electric field or magnetic field. p_label
is
an integer array indexing the position of each grid point, whether it is
inside the surrounding medium or a particle, and in which
layer.
dumpNFs2HDF5File(fname, groupname, filename, incidences, Epower, wavelen, work, p_label)
Exports electric and magnetic near fields into a HDF5 file.
filename
, fname
, groupname
are
the names of the HDF5 file, group, and subgroup name, respectively.
incidences, Epower, wavelen, work
are the incidence angles,
selected power for mapping fields, wavelength, and near-field
quantities, respectively. Ef
is a logical flag which
selects either electric field or magnetic field. p_label
is
an integer array indexing the position of each grid point, whether it is
inside the surrounding medium or a particle, and in which
layer.
countLines(filename) result(nlines)
Counts lines in a text file.
multiscat
module
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
\(\langle|\mathbf{E}|^2\rangle,
\langle|\mathbf{B}|^2\rangle\) 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 \(\langle|\mathbf{E}|^2\rangle\) and \(\langle|\mathbf{B}|^2\rangle\) 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
(\(l_{max}\times \text{nscat}\), \(l_{max}\times\text{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, \(c_j\) (centred at \(j\)) and \(c_i\) (centred at \(i\)), such that \(c_j = \text{Mat}\, c_i\). Logical inputs
jregt
and iregt
specify whether \(c_j\) and \(c_i\) 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, \dots, n_{max}\).
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\leq 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 \leq
n_{max}\) when \(\operatorname{Tr}(\Re{(\text{Tcol})})\)
converges to \(\text{rtol\_G}:=10^{-\text{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 \(k\mathbf{r}\)(3) from which VTACs will be
generated. Regular or irregular VTACs will be generated depending on
whether \(k\mathbf{r}\)(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 (\(\alpha,
\beta, \gamma\)). 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
\(\mathtt{op(A) = A}\) (no operation),
whereas op = ’C’
corresponds to \(\mathtt{op(A) = A^\dagger}\).
RotMatX(ang) result(rotMat)
Calculates a rotation matrix along the x axis using input argument
angle(ang
).
RotMatY(ang) result(rotMat)
Calculates a rotation matrix along the y axis using input argument
angle(ang
).
RotMatZ(ang) result(rotMat)
Calculates a rotation matrix along the z axis using input argument
angle(ang
).
rotZYZmat(angles) result(mat)
Calculates rotation matrix mat
for ZY’Z’, using the Euler
angles
=(\(\alpha,\beta,\gamma\))
List of common arguments
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 [\(n_1\), \(n_2\), 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
(\(\text{npts}\times 4\times
\text{nwavelen}\)): contains the orientation averaged value of
\(\overline{\mathscr{C}}\) at different
grid points and wavelengths.
orAvextE_int
(\(\text{npts}\times \text{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\times
\text{nscat}\times 4\times \text{nwavelen}\times \text{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\times
\text{n}\times \text{nwavelen}\)): a matrix consisting of
orientation-averaged cross-sections and CD at different wavelengtha. The
first column gives the values for \(n_{max}\) and other columns contain values
for different value of \(n=1, \dots,
\text{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 (\(\in [0,1,2,3]\)) (the default
value is 1).
wavelen
: a vector of specified
wavelength(s).
miet
module
This module contains routines for calculating one-body T-matrices (currently limited to spherical scatterers, using Mie theory).
calcMieTMat(x, s, zeropad_, tmat)
Calculates the diagonal T-matrix of a spherical scatterer for a
given size parameter x=kR
, relative refractive indices
(\(s = k_{in}/k_{out}\));
zeropad_=nmax
maximum value of the multipole index inferred
from tmat
’s dimensions.
calcMieCoeffs(x, s, gammas, deltas)
Calculates the Mie coefficients for a spherical scatterer as defined by
equations H.46 and H.47 of Ref. . The coefficients are interpreted as
magnetic and electric susceptibilities (\(\Gamma_n\) and \(\Delta_n\), respectively) of the scattered
field. Note the relation to standard Mie coefficients: \(a_n = -\Delta_n\) and \(b_n = -\Gamma_n\).
calcCoatMieCoeffs(x, s, gammas, deltas)
Calculates the Mie coefficients for a coated sphere based on the
equations H.110 and H.113 of Etchegoin and Le Ru.
calcStoutCoeffs(x, rri, nmax, Cn, Dn)
Calculates the Cn,Dn
coefficients as defined by equation
(50) in Stout. These coefficients are used to calculate absorption
cross-sections. rri
is the relative refractive index,
nmax
is the maximum value of the multipole index.
calcMieIntCoeffs(a, k, scaCoeffs, intCoeffsReg, intCoeffsIrr, csAbs)
Calculates the regular and irregular VSWF coefficients for the field
inside each concentric region of a (layered) Mie scatterer. The formulae
are based on Eqs. H.117–H.123 of Etchegoin and Le Ru. a
,
k
, and scaCoeffs
are vectors of the radius of
the concentric interfaces, relative refractive index, and scattered
field coefficients for the host medium, respectively.
intCoeffsReg
and intCoeffsIrr
are matrices of
regular and irregular field coefficients for each concentric region
inside the scatterers and csAbs
contains the partial
absorptions calculated using equation (29) in Mackowski.
swav
module
This module contains routines for calculating and transforming scalar
(SSWs) and vector spherical waves (VSWFs). It depends on
Amos (toms644.f)
to calculate spherical Bessel and Hankel
functions using recurrence. In order to limit redundancy, parameter
definitions are renewed only where they are changed.
calcVTACs(r0, k, regt, vtacs)
Calculates the irregular (if regt=.false.
) or the regular
(if regt=.true.
) vector translation-addition coefficients
for a given kr0
.r0
is a relative position vector, k
is the
wavenumber, regt
is a logical argument which determines the
type: regular or irregular, and vtacs(1:2*pmax,1:2*pmax)
is
the input/output array.
calcSTACs(r0, k, pmax, regt, scoeff)
Calculates the scalar translation-addition coefficients.(\(\alpha_{nu,mu;n,m}\) or \(\beta_{nu,mu;n,m}\)). The output
corresponds to the scalar translation-addition coefficients
\alpha
(irregular, for regt=.false.
) or
\beta
(regular, for regt=.true.
).pmax
is a maximal composite index and
scoeff(0:pmax,0:pmax)
is the coefficients matrix.
calcVTACsAxial(r0, k, pmax, regt, flip, mqn_, vtacs)
Calculates the irregular (if regt=.false.
) or the regular
(if regt=.true.
) vector translation-addition coefficients
for a given kr0
, along the z-axis.r0
is the z-axial displacement distance, flip is a logical
argument, mqn_
is a logical argument for changing from
qnm
to mqn
indexing, and
vtacs(1:2*pmax,1:2*pmax)
is the matrix of
coefficients.
calcSTACsAxial(r0, k, pmax, regt, flip, stacs)
Calculates the normalised scalar translation-addition coefficients along
the z-axis for a given kr0
.r0
is a displacement distance and
stacs(0:pmax,0:pmax)
is the coefficients matrix
corresponding to \(\alpha\) (irregular,
for regt=.false.
) or \(\beta\).
calcVSWs(r, k, pmax, regt, cart, waves, wavesB)
Calculates (at r
) the normalised vector spherical waves,
\(M_{nm}\) and \(N_{nm}\) for evaluation of electric and
magnetic fieldsr(3)
is the cartesian coordinate of a point in 3D;
cart
is a logical argument which triggers conversion to
cartesian coordinates; waves(2*pmax,3)
contains elements
(\(M_{nm}\) and \(N_{nm}\)) of the abstract column vector
defined in Eq. B1 of Ref. and wavesB(2*pmax,3)
is similar
to waves
, only swapping the position of \(M_{nm}\) and \(N_{nm}\) and multiplying by \(-ik\) for calculation of the magnetic
field.
calcSSWs(xyz, k, pmax, regt, psi)
Calculates (at xyz
) the scalar spherical waves
\psi_{nm}
.xyz
is the cartesian coordinates of a point in 3D;
psi(0:pmax)
contains elements of the spherical waves
\psi_{nm}
as defined by equation 13a in Chew.
calcJCoeffsPW(ipwE0, kVec, xyz, ipwCoeffsJ)
Translates the supplied ipwCoeffs
coefficients to different
centres for an incident plane wave.ipwE0(3)
is the incident plane wave’s amplitude vector,
kVec(3)
is the incident wave vector,
xyz(3,nscat)
is a matrix containing the centre of different
scatterers, and ipwCoeffsJ
(\(\text{nscat} \times \text{lmax}\)) is a
vector containing the translated incident plane wave coefficients to the
centre of different scatterers.
calcCoeffsPW(ipwE0, ipwDirn, ipwCoeffs)
Calculates the coefficients for expressing an incident plane wave in
terms of vector spherical waves \(M_{nm}\) and \(N_{nm}\).ipwDirn
(3) is the normalised direction vector of the
incident plane wave and ipwCoeffs
(2*pmax) contains
coefficients for expressing an incident plane wave in terms of vector
spherical waves \(M_{nm}\) and \(N_{nm}\), up to a maximum \(n_{max}\).
offsetCoeffsPW(a, kVec, xyzr, aJ)
Translates the VSWF coefficients (a
) of an incident plane
wave (centred at the origin) to another origin.a
(\(l_{max}\)) contains
coefficients for a regular VSWF expansion centred at the origin for an
incident plane wave, xyzr
includes centres of different
scatterers, and aJ
contains scatterer centred
coefficients.
calcWignerBigD(angles, pmax, bigD)
Calculates the Wigner D-functions (\(D^s_{m,n}(\alpha,\beta,\gamma)\)).angles
(3) includes (\(\alpha,\beta,\gamma\)) in radians and
bigD
(pmax,pmax) contains Wigner D-coefficients.
calcWignerLittled(theta, pmax, d)
Calculates the Wigner d-functions (\(d^s_{m,n}(\theta)\)).theta
is angle in radians and d
(0:pmax,0:pmax)
are values for \(d^s_{m,n}\) in block
diagonal matrix form.
calcWignerd0andMore(x, pmax, d, pi, tau)
Calculates the Wigner d-functions for \(n=0\) and also computes the derivative
functions for optional outputs pi
and
tau
.x
is \(\cos(\theta)\),
d
(0:pmax), pi
(0:pmax),
tau
(0:pmax) contain values for \(d^s_{m,0}\), \(\pi_{m,s}\), and \(\tau_{m,s}\) respectively.
calcRiccatiBessels(z, nmax, regt, f, df)
Calculates the Riccati-Bessel functions \(\psi_n\) (if regt=.true.
) or
\(\xi_n\) (regt=.false.
),
and their derivatives, for \(n=1,\dots,n_\text{max}\).z
is a scalar complex argument, f
(1:nmax) is a
matrix containing Riccati-Bessel functions \(\psi_n(z)=z*j_n(z)\) or \(\xi_n(z)=z*h_n(z)\) for \(n=1,\dots,n_\text{max}\), and
df
(1:nmax) are the corresponding derivatives of
f
.
calcSphBessels(z, nmax, regt, bes)
A wrapper routine for computing spherical Bessel/Hankel functions of the
first kind for a complex argument z
.bes(0:nmax)
contain Bessel (\(J_{n+1/2}\)) or Hankel (\(H_{n+1/2}\)) function (of \(1^{st}\) kind) values for \(n=0,\dots,n_\text{max}\) for a complex
argument z
.
xyz2rtp(xyz, rtp, cth)
Transforms the cartesian coordinates (x,y,z)
of a point in
3D space to spherical polar coordinates
(r,\theta,\phi)
xyz(3)
is a vector of cartesian coordinates,
rtp(3)
is a vector of spherical polar coordinates, and
cth
is cos\((\theta)\).
rtp2xyz(rtp, xyz)
The inverse of xyz2rtp
. Transforms the spherical polar
coordinates (r,\theta,\phi)
of a point in 3D space to
cartesian coordinates \((x,y,z)\).
calcVTrtp2xyz(rtp, transform)
Calculates the matrix of transformation from a vector in spherical
coordinates to a vector in cartesian coordinates at point
(r,\theta,\phi)
(in spherical polar coordinates).
calcVTxyz2rtp(rtp, transform)
The inverse of calcVTrtp2xyz
. Calculates the matrix of
transformation from a vector in cartesian coordinates to a vector in
spherical coordinates at point (r,\theta,\phi)
(in
spherical polar coordinates).
calcAbsMat(Xi, ro, mat)
Calculates the absorption matrix \(\Gamma_j =
\mathtt{mat(l_{max},l_{max})}\) for the input arguments
Xi
and ro
(Eq. (49) of Stout). \(\Gamma_j\) is used in the evaluation of the
orientation-averaged absorption cross-section inside each
particle.
calcLamMat(Xi, ro, mat)
Calculates the "Lambda" matrix \(\Lambda_j =
\mathtt{mat(l_{max},l_{max})}\) for the input arguments
Xi
and ro
(Eq. (53) of Stout). \(\Lambda_j\) is used in the evaluation of
the orientation-averaged internal electric field inside homogeneous
spheres.
nm2p(n, m, l)
Calculates a generalised index l=n(n+1)+m
, for a unique
(n,m)
, (Vector spherical harmonics are spanned by two
indices: \(n\) and \(m\), such that \(0 \leq n \leq n_{max}\) and \(-n \leq m \leq n\)).n,m,l
are integers.
p2nm(p, n, m)
Calculates unique (n,m)
from a given composite index
p
.p
is a real value.
nm2pv2(n, m, p)
Some recurrences are defined only for \(m \geq
0\), in which case we shall use a second version of the composite
index \(p_{v2}=n(n+1)/2+m\).
testPmax(name, pmax, nmax)
Tests pmax
for commensurability, i.e. is \(p_{max}==n_{max}(n_{max}+2)\) and \(n_{max}=m_{max}\)? If not, then the program
will be stopped.
sphmsv
module
This module contains routines for calculating Stokes incident vector, Stokes phase matrix and scattering matrix for an input T-matrix. The formulae are based on Mishchenko.
calcStokesIncVec(ehost_, ipwDirn_, ipwAmpl_, verb_, Stokes_Vec)
Calculates the Stokes incident vector Stokes_Vec
.
calcStokesPhaseMat(SMat, verb_, Z)
Calculates the Stokes phase matrix for the specified incident and
scattered angles. SMat(2,2)
and Z(4,4)
are the
scattering and Stokes phase matrices which follow Eqs. (5.11-14) and
(2.106-121) of Mishchenko.
calcScatMat(tmat, host_K, spwDirn_, ipwDirn_, verb_, SMat)
Calculates the scattering matrix using the T-matrix, for the
specified incident and scattering angles.
linalg
module
This module contains wrappers to drive LAPACK’s square-matrix inversion routines and linear solvers.
invSqrMat(trans_, verb_ A)
Calculates inverse of a complex-valued square matrix
A(n,n)
, using the ZGETRF and ZGETRI routines in LAPACK.
A
is overwritten by inv(A) on the output.
trans_
is an optional logical input, in case
.true.
the routine considers transpose of A and finally
returns the transpose of the inverted matrix as the output.
verb_
: an optional input of the verbosity value.
solLinSys(isol_, verb_, A, X)
Solves a complex-valued linear system of equations \(Ax=b\), where A
(n,n) is a
square matrix, b
(n) is a known vector, and
x
(n) is the vector to be determined. Depending on the value
isol_
, calls solLinSysV
or
solLinSysVX
. Both A
and X
are
overwritten on output.
solLinSysV(verb_, A, X)
For solving a linear system, uses LAPACK’s "simple" driver
ZGESV.
solLinSysVX(verb_, A, X)
For solving a linear system, uses LAPACK’s "simple" driver
ZGESVX.
eps
module
This module contains wavelength-dependent dielectric functions
epsXX(lambda)
for various materials including Au, Ag, Al,
Cr, Pd, Pt, Si, and Water).
interp1( x1, y1, x2, y2 )
y2
using the input values
x1,y1
at the points x2
.epsAu(wavelength) result(eps)
Returns the wavelength-dependent relative dielectric function of gold.
This function uses the analytical expression given in Eq. (E.2) of
Etchegoin and Le Ru.
epsAg(wavelength) result(eps)
Returns the wavelength-dependent relative dielectric function of silver.
This function uses the analytical expression given in Eq. (E.1) of
Etchegoin and Le Ru.
epsPt(wavelength) result(eps)
Returns the wavelength-dependent relative dielectric function of a
Lorentz-Drude metal, with the parameters for Pt from Rakic.
epsPd(wavelength) result(eps)
Returns the wavelength-dependent relative dielectric function of a
Lorentz-Drude metal, with the parameters for Pd from Rakic.
epsSi(wavelength) result(eps)
Returns the wavelength-dependent relative dielectric function of Silicon
in the range 206.6 nm to 1200.0 nm interpolated from Aspnes.
epsAl(wavelength) result(eps)
Returns the wavelength-dependent relative dielectric function of
Aluminum in the range 103.32 nm to 2755.2 nm from Rakic.
epsCr(wavelength) result(eps)
Returns the wavelength-dependent relative dielectric function of
Aluminum in the range 100.8 nm to 31 \(\mu
m\), from tabulated data.
epsWater(wavelength) result(eps)
Returns the wavelength-dependent relative dielectric function of Water
at temperature \(20^o\text{C}\) in the
range 200 nm to 3000 nm from Daimon.
epsDrude(wavelength, eps_infty, lambda_p, mu_p) result(eps)
Returns the wavelength-dependent relative dielectric function of a Drude
metal.
HDFfive
module
This module contains subroutines for reading and writing data in HDF5 format.
h5_crtgrp(filename_, main_grpname, subgrpsname)
This subroutine creates subgroups in an existing group.
h5_wrtvec2file(filename_, groupname, dsetname, dset_data)
This subroutine writes vector data in a dataset in an existing
group.
h5_wrt2file(filename_, groupname, dsetname, dset_data)
This subroutine writes data in a dataset in an existing group.
h5_wrt_attr(attribute, dataset_id)
This subroutine adds an attribute to an existing dataset, typically a
brief explanation about the contents of the dataset.