Файл: Doicu A., Wriedt T., Eremin Y.A. Light scattering by systems of particles (OS 124, Springer, 2006.pdf

182 2 Null-Field Method

with

The T -operator solution is a formal solution since we assumed the invertibility of the integral equation (2.219). The usefulness of the formalism in numerical applications depends upon the possibility of discretizing (2.219) in a suitable way. This can be done for a class of simple infinite surfaces as for example plane surfaces or surfaces with periodic roughness and small amplitude roughness. It should be noted that for plane surfaces we obtain exactly the results derived in our earlier analysis.

3

Simulation Results

In this chapter we present computer simulation results for the scattering problems considered in Chap. 2. For the numerical analysis, we use our own software package and some existing electromagnetic scattering programs. After a concise description of the T -matrix code, we present the theoretical bases of the electromagnetics programs which we used to verify the accuracy of the new implementation. We then present simulation results for homogeneous, axisymmetric and nonaxisymmetric particles, inhomogeneous, layered and composite particles, and clusters of particles. The last sections deal with scattering by a particle on or near a plane surface and the computation of the e ective wave number of a half-space with randomly distributed particles.

3.1 T -matrix Program

A Fortran computer program has been written to solve various scattering problems in the framework of the null-field method. This section gives a short description of the code, while more details concerning the significance of the input and output parameters are given in the documentation on CD-ROM.

The main program TMATRIX.f90 calls a T -matrix routine for solving a specific scattering problem. These routines computes the T -matrix of:

–Homogeneous, dielectric (isotropic, chiral) and perfectly conducting, axisymmetric particles (TAXSYM.f90),

–Homogeneous, dielectric (isotropic, uniaxial anisotropic, chiral) and perfectly conducting, nonaxisymmetric particles (TNONAXSYM.f90)

–Axisymmetric, composite particles (TCOMP.f90)

–Axisymmetric, layered particles (TLAY.f90)

–An inhomogeneous, dielectric, axisymmetric particle with an arbitrarily shaped inclusion (TINHOM.f90)

–An inhomogeneous, dielectric sphere with a spherical inclusion (TINHOM2SPH.f90)

184 3 Simulation Results

–An inhomogeneous, dielectric sphere with an arbitrarily shaped inclusion (TINHOMSPH.f90)

–An inhomogeneous, dielectric sphere with multiple spherical inclusions (TINHOMSPHREC.f90)

–Clusters of arbitrarily shaped particles (TMULT.f90)

–Two homogeneous, dielectric spheres (TMULT2SPH.f90)

–Clusters of homogeneous, dielectric spheres (TMULTSPH.f90 and TMULTSPHREC.f90)

–Concentrically layered sphere (TSPHERE.f90)

–A homogeneous, dielectric or perfectly conducting, axisymmetric particle on or near a plane surface (TPARTSUB.f90)

Three other routine are invoked by the main program

–SCT.f90 computes the scattering characteristics of a particle using the previously calculated T matrix

–SCTAVRGSPH.f90 computes the scattering characteristics of an ensemble of polydisperse, homogeneous spherical particles

–EFMED.f90 computes the e ective wave number of a half-space with randomly distributed particles

Essentially, the code performs a convergence test and computes the T -matrix and the scattering characteristics of particles with uniform orientation distribution functions.

An important part of the T -matrix calculation is the convergence procedure over the maximum expansion order Nrank, maximum azimuthal order Mrank and the number of integration points Nint. In fact, Nrank, Mrank and Nint are input parameters and their optimal values must be found by the user. This is accomplished by repeated convergence tests based on the analysis of the di erential scattering cross-section as discussed in Sect. 2.1.

The scattering characteristics depend on the type of the orientation distribution function. By convention, the uniform distribution function is called complete if the Euler angles αp, βp and γp are uniformly distributed in the intervals (0, 360◦), (0, 180◦) and (0, 360◦), respectively. The normalization constant is 4π for axisymmetric particles and 8π2 for nonaxisymmetric particles. The uniform distribution function is called incomplete if the Euler angles αp, βp and γp are uniformly distributed in the intervals (αp min, αp max), (βp min, βp max), and (γp min, γp max), respectively. For axisymmetric particles, the orientational average is performed over αp and βp, and the normalization constant is

(αp max − αp min) (cos βp min − cos βp max) ,

while for nonaxisymmetric particles, the orientational average is performed over αp, βp and γp, and the normalization constant is

(αp max − αp min) (cos βp min − cos βp max) (γp max − γp min) .

The scattering characteristics computed by the code are summarized later.

3.1.1 Complete Uniform Distribution Function

For the complete uniform distribution function, the external excitation is a vector plane wave propagating along the Z-axis of the global coordinate system and the scattering plane is the XZ-plane. The code computes the following orientation-averaged quantities:

– The scattering matrix F at a set of Nθ,RND scattering angles

–The extinction matrix K

–The scattering and extinction cross-sections Cscat and Cext and

–The asymmetry parameter cos Θ .

The scattering angles, at which the scattering matrix is evaluated, are uniformly spaced in the interval (θmin,RND, θmax,RND). The elements of the scattering matrix are expressed in terms of the ten average quantities

for macroscopically isotropic and mirror-symmetric media. Spq Sp1q1 are computed at Nθ,GS scattering angles, which are uniformly spaced in the interval (0, 180◦). The scattering matrix is calculated at the same sample angles and polynomial interpolation is used to evaluate the scattering matrix at any polar angle θ in the range (θmin,RND, θmax RND). The average quantities Spq Sp1q1 can be computed by using a numerical procedure or the analytical orientation-averaging approach described in Sect. 1.5. The numerical orientation-averaging procedure chooses the angles αp, βp and γp to sample the intervals (0◦, 360◦), (0◦, 180◦) and (0◦, 360◦), respectively. The prescription for choosing the angles is to

– Uniformly sample in αp

– Uniformly sample in cos βp or nonuniformly sample in βp

– Uniformly sample in γp

The integration over αp and γp are performed with Simpson’s rule, and the number of integration points Nα and Nγ must be odd numbers. The integration over βp can also be performed with Simpson’s rule, and in this specific case, the algorithm samples uniformly in cos βp, and the number of integration points Nβ is an odd number. Alternatively, Gauss–Legendre quadrature method can be used for averaging over βp, and Nβ can be any integer number. For the analytical orientation-averaging approach, the maximum expansion and azimuthal orders Nrank and Mrank (specifying the dimensions of the T

186 3 Simulation Results

matrix) can be reduced. In this case, a convergence test over the extinction and scattering cross-sections gives the e ective values Nranke and Mranke .

The orientation-averaged extinction matrix K is computed by using (1.125) and (1.126), and note that for macroscopically isotropic and mirrorsymmetric media the o -diagonal elements are zero and the diagonal elements are equal to the orientation-averaged extinction cross-section per particle.

For macroscopically isotropic and mirror-symmetric media, the orientationaveraged scattering and extinction cross-sections Cscat = Cscat I andCext = Cext I are calculated by using (1.124) and (1.122), respectively, while for macroscopically isotropic media, the code additionally computesCscat V accordingly to (1.133), and Cext V as Cext V = K14 .

The asymmetry parameter cos Θ is determined by angular integration over the scattering angle θ. The number of integration points is Nθ,GS and Simpson’s rule is used for calculation. For macroscopically isotropic and mirror-symmetric media, the asymmetry parameter cos Θ = cos Θ I is calculated by using (1.134), while for macroscopically isotropic media, the code supplementarily computes cos Θ V accordingly to (1.135).

The physical correctness of the computed results is tested by using the inequalities (1.138) given by Hovenier and van der Mee [103]. The message that the test is not satisfied means that the computed results may be wrong.

The code also computes the average di erential scattering cross-sectionsσdp and σds for a specific incident polarization state and at a set of Nθ,GS scattering angles (cf. (1.136) and (1.137)). These scattering angles are uniformly spaced in the interval (0◦, 180◦) and coincide with the sample angles at which the average quantities Spq Sp1q1 are computed. The polarization state of the incident vector plane wave is specified by the complex amplitudes Ee0,β and Ee0,α, and the di erential scattering cross-sections are calculated

|Ee0,β |2 + |Ee0,α|2

Note that Nθ,GS is the number of scattering angles at which Spq Sp1q1 , σdp and σds are computed, and also gives the number of integration points for calculating Cscat V, cos Θ I and cos Θ V. Because the integrals are computed with Simpson’s rule, Nθ,GS must be an odd number.

3.1.2 Incomplete Uniform Distribution Function

For the incomplete uniform distribution function, we use a global coordinate system to specify both the direction of propagation and the states of polarization of the incident and scattered waves, and the particle orientation. A special orientation with a constant orientation angle can be specified by setting Nδ = 1 and δmin = δmax, where δ stands for αp, βp and γp. For example,

uniform particle orientation distributions around the Z-axis can be specified by setting Nβ = 1 and βp min = βp max.

The scattering characteristics are averaged over the particle orientation by using a numerical procedure. The prescriptions for choosing the sample angles and the significance of the parameters are as in Sect. 3.1.1. For each particle orientation we compute the following quantities:

–The phase matrix Z at Nϕ scattering planes

–The extinction matrix K for a plane wave incidence

–The scattering and extinction cross-sections Cscat and Cext for incident parallel and perpendicular linear polarizations

–The mean direction of propagation of the scattered field g for incident parallel and perpendicular linear polarizations

The azimuthal angles describing the positions of the scattering planes at which the phase matrix is computed are ϕ(1), ϕ(2), . . . , ϕ(Nϕ). In each scattering plane i, i = 1, 2, . . . , Nϕ, the number of zenith angles is Nθ (i), while the zenith angle varies between θmin(i) and θmax(i). For each particle orientation, the phase and extinction matrices Z and K are computed by using (1.77) and (1.79), respectively, while the scattering and extinction cross-sections Cscat and Cext are calculated accordingly to (1.101) and (1.100), respectively.

The mean direction of propagation of the scattered field g is evaluated by angular integration over the scattering angles θ and ϕ (cf. (1.87)), and the numbers of integration points are Nθ,asym and Nϕ,asym. The integration over ϕ is performed with Simpson’s rule, while the integration over θ can be performed with Simpson’s rule or Gauss–Legendre quadrature method. We note that the optical cross-sections and the mean direction of propagation of the scattered field are computed for linearly polarized incident waves (vector plane waves and Gaussian beams) by choosing αpol = 0◦ and αpol = 90◦ (cf. (1.18)).

The code also calculates the average di erential scattering cross-sections σdp and σds in the azimuthal plane ϕGS and at a set of Nθ,GS scattering angles. The average di erential scattering cross-sections can be computed for scattering angles ranging from 0◦ to 180◦ in the azimuthal plane ϕGS and from 180◦ to 0◦ in the azimuthal plane ϕGS + 180◦, or for scattering angles ranging from 0◦ to 180◦ in the azimuthal plane ϕGS. The calculations are performed for elliptically polarized vector plane waves (characterized by the complex polarization unit vector epol as in (3.1)) and linearly polarized Gaussian beams (characterized by the polarization angle αpol).

Mrank = Nrank.

The codes perform calculations with doubleor extended-precision floating point variables. The extended-precision code is slower than the double-