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

188 3 Simulation Results

CPU time consumption rapidly increases with increasing particle size and asphericity.

The program is written in a modular form, so that modifications, if required, should be fairly straightforward. In this context, the integration routines, the routines for computing spherical and associated Legendre functions or the routines for solving linear systems of equations can be replaced by more e cient routines.

The main “drawback” of the code is that the computer time requirements might by higher in comparison to other more optimized T -matrix programs. Our intention was to cover a large class of electromagnetic scattering problem and therefore we sacrifice the speed in favor of the flexibility and code modularization. The code shares several modules which are of general use and are not devoted to a specific application.

3.2 Electromagnetics Programs

The on-line directories created by Wriedt [266] and Flatau [67] provide links to several electromagnetics programs. In order to validate our software package we consider several computer programs relying on di erent methods. This section briefly reviews these methods without entering into deep details. More information can be found in the literature cited.

3.2.1 T -matrix Programs

Computer programs using the null-field method to calculate the electromagnetic scattering by homogeneous, axisymmetric particles have been given by Barber and Hill [8]. These programs compute:

–The angular scattering over a designated scattering plane or in all directions for a particle in a fixed orientation

–The angular scattering and the optical cross-sections for an ensemble of particles randomly oriented in a 2D plane or in 3D

–The scattering matrix elements for an ensemble of particles randomly oriented in 3D

–The normalized scattering cross-sections versus size parameter for a particle in fixed or random orientation

–The internal intensity distribution in the equatorial plane of a spheroidal particle

We would like to mention that the general program structure and some programming solutions given by Barber and Hill have been adopted in our computer code.

The T -matrix codes developed by Mishchenko and Travis [167] and Mishchenko et al. [168, 169] are e cient numerical tools for computing electromagnetic scattering by homogeneous, axisymmetric particles with size

significantly larger than a wavelength. The following computer programs are available on the World Wide Web at www.giss.nasa.gov/ crimim:

–A Fortran code for computing the amplitude and phase matrices for a homogeneous, axisymmetric particle in an arbitrary orientation

–A Fortran code for computing the far-field scattering and absorption characteristics of a polydisperse ensemble of randomly oriented, homogeneous, axisymmetric particles

The last code is suitable for practical applications requiring the knowledge of size-, shape-, and orientation-averaged quantities such as the optical cross-sections and phase and scattering matrix elements. The implementation of the analytical orientation-averaging procedure makes this code the fastest computer program for randomly oriented, axisymmetric particles. The following computer programs are also provided at the w site www.giss.nasa.gov/ crimim:

–A Lorenz–Mie code for computing the scattering characteristics of an ensemble of polydisperse, homogeneous, spherical particles and

–The Fortran code SCSMTM for computing the T -matrix of a sphere cluster and the orientation-averaged scattering matrix and optical crosssections [153]

3.2.2 MMP Program

MMP is a Fortran code developed by Hafner and Bomholt [94] which uses the multiple multipole method for solving arbitrary 3D electromagnetic problems. In the multiple multipole method, the electromagnetic fields are expressed as linear combinations of spherical wave fields corresponding to multipole sources. By locating these sources away from the boundary, the multipole expansions are smooth on the surface and the boundary singularities are avoided. The multipoles describing the internal field are located outside the particle, while the multipoles describing the scattered field are positioned inside the particle. Not only spherical multipoles can be used for field expansions; other discrete sources as for instance distributed vector spherical wave functions, vector Mie potentials or magnetic and electric dipoles can be employed. Therefore, other names for similar concepts have been given, e.g., method of auxiliary sources [275], discrete sources method [49, 62], fictitious sources method [140, 158], or Yasuura method [113]. An overview of the discrete sources method has been given by Wriedt [267] and a review of the latest literature in this field has been published by Fairweather et al. [63]. By convention, we categorize the multiple multipole method as a discrete sources method with multiple vector spherical wave functions. Below we summarize the basic concepts of the discrete sources method for the transmission boundary-value problem.

Let {Ψ 3µ, Φ3µ} be a system of radiating solutions to the Maxwell equations in Ds with the properties × Ψ 3µ = ksΦ3µ and × Φ3µ = ksΨ 3µ. Analogously,

system of vector functions (3.2), the expansion coe cients aµs,iN and bsµ,iN , µ = 1, 2, . . . , N , can be determined by solving the minimization problem,

This minimization problem leads to a system of normal equations for the expansion coe cients. Experience has shown that this technique often gives inaccurate and numerically unstable results with highly oscillating error distributions along the boundary surface. Accurate and stable results can be obtained by using the point matching method with an overdetermined system of equations, that is, the boundary conditions are fulfilled at a set of surface points, while the number of surface points exceeds the number of unknowns (more than twice). The resulting system of equations is solved in the least-squares sense by using, for example, the QR-factorization.

The main problem in the multiple multipole method is the choice of the number, position and the order of multipoles, and the distribution of matching points. Some useful criteria can be summarized as follows. Firstly, the matching points on the boundary should be set with enough density, and the intervals between adjacent pairs of matching points must be far less than the wavelength, for suppressing errors of expansion on the boundary. Secondly, the poles should be neither too near nor too far from the boundary, since the former requirement may avoid alternative errors corresponding to nonphysical rough solution on that interval; and a well-conditioned matrix will be formulated with the latter requirement. In view of the quasi-local behavior of a multipole expansion, only a restricted area around its origin is influenced by this pole, so that a reasonable scheme is to let the area of influence of each pole fit a portion of the boundary, but the whole boundary must be covered by the areas of influence of all the poles. Finally, the highest order of a pole is limited by the density of matching points covered by its area of influence. It should be noted that the MMP code contains routines for automatically optimizing the number, location and order of multipoles. More detailed overviews on methodology, code and simulation technique can be found in the works of Ludwig [147, 148], Hafner [91] and Bomholt [18]. The attractiveness and simplicity of the physical idea of the MMP and the public availability of the code by Hafner [92, 93] have resulted in widespread applications of this technique.

In addition to the MMP code, the DSM (discrete sources method) code developed by Eremin and Orlov [61] will be used for computer simulations. This DSM code is devoted to the analysis of homogeneous, axisymmetric particles using distributed vector spherical wave functions. For highly elongated particles, the sources are distributed along the axis of symmetry of the particle, while for highly flattened particles, the sources are distributed in the complex plane (see Appendix B).

192 3 Simulation Results

3.2.3 DDSCAT Program

DDSCAT is a freely available software package which applies the discrete dipole approximation (DDA) to calculate scattering and absorption of electromagnetic waves by particles with arbitrary geometries and complex refractive indices [55]. The discrete dipole approximation model the particle as an array of polarizable points and DDSCAT allows accurate calculations for particles with size parameter ksa < 15, provided the refractive index mr is not large compared to unity, |mr − 1| < 1. The discrete dipole approximation (sometimes referred to as the coupled dipole method (CDM)) was apparently first proposed by Purcell and Pennypacker [197]. The theory was reviewed and developed further by Draine [53], Draine and Goodman [56] and Draine and Flatau [54]. An improvement of this method was given by Piller and Martin [194] applying concepts from sampling theory, while Varadan et al. [237], Lakhtakia [130] and Piller [192] extended the discrete dipole approximation to anisotropic, bi-anisotropic and high-permittivity materials, respectively. Since the discrete dipole approximation can be derived from the volume integral equation we follow the analysis of Lakhtakia and Mulholland [131] and review the basic concepts of the volume integral equation method and the discrete dipole approximation.

The volume integral equation method relies on the integral representation for the electric field

The domain Di is discretized into simply connected elements (cells) Di,m,

m = 1, 2, . . . , Ncells, and Di = Nmcells=1 Di,m. Each element Di,m is modeled as being homogeneous such that mr(r) = mr,m for r Di,m, even though this

assumption can lead to an artificial material discontinuity across the interface of two adjacent elements. As a consequence of the discretization process, the element equation reads as

Due to the singularity at r = r it is necessary to approximate the integral

analytically in a small volume around the singularity. With Dr being a sphere of radius r around the singularity, we transform the integral as follows: