An improved algorithm for calculation of light scattering by a coated sphere has been recently presented by Liu Lei et al 2007. This algorithm improves the speed of the algorithm of Toon OB and Ackerman 1981 when the relative particle size exceeds ~100.
In the improved algorithm, the an and bn coefficients of a series expansion of the scattering properties (see Equations 3 in Mie theory: Overview), are calculated as follows:
|
an = Mn (x1) Cn (x1, x2, m1, m2)
|
(1)
|
|
bn = Mn (x1) Dn (x1, x2, m1, m2)
|
(2)
|
where, n is the index, running from 1 to nmax, x1, x2, m1, m2 are respectively the relative sizes (x) and complex refractive indices (m) of the shell and core of the coated sphere, and
|
Cn =
|
|
[Ln(y1) - m1 Ln(x1)] + Qn(y1, y3) [m1 Ln(x1) - Nn(y1)]
|
|
[Ln(y1) - m1 Nn(x1)] + Qn(y1, y3) [m1 Nn(x1) - Nn(y1)]
|
|
|
|
|
|
[m1 Ln(y2) - m2 Ln(y3)] / [m1 Ln(y2) - m2 Nn(y3)]
|
|
[m1 Ln(y2) - m2 Ln(y3)] / [m1 Ln(y2) - m2 Nn(y3)]
|
|
(3)
|
|
Dn =
|
|
[m1 Ln(y1) - Ln(x1)] + Qn(y1, y3) [Ln(x1) - m1 Nn(y1)]
|
|
[m1 Ln(y1) - Nn(x1)] + Qn(y1, y3) [Nn(x1) - m1 Nn(y1)]
|
|
|
|
|
|
[m2 Ln(y2) - m1 Ln(y3)] / [m2 Ln(y2) - m1 Nn(y3)]
|
|
[m2 Ln(y2) - m1 Ln(y3)] / [m2 Ln(y2) - m1 Nn(y3)]
|
|
(4)
|
where y1 = m1x1, y2 = m2x2, y3 = m1x2.
Function Ln(z) is a logarithmic derivative of the Ricatti-Bessel function, for example, http://mathworld.wolfram.com/Riccati-BesselFunctions.html, which can be calculated with the downward recursion [see Mie theory: Algorithms for functions of the refractive index and/or size of the sphere, Ln(z) is defined as An(z) there]. The final value Lnmax for the downward recursion can be found by using a continued-fraction algorithm of Lentz WJ 1976.
Function Nn(z) is also a logarithmic derivative of the Ricatti-Bessel function, which can be calculated with the upward recursion:
The initial value of Nn(z) for the upward recursion is N0(z) = -i (when the imaginary part of the variable z is negative). The results of calculation of Nn(z) indicate that the upward recursion for this function is stable and accurate when the upward recurrence is executed with the argument z being real or complex (Mackowski DW et al 1990).
Function Mn(z) is the ratio of the Ricatti-Bessel functions, which is calculated by using the upward recursion:
|
Mn(z) = Mn-1(z)
|
|
Ln-1(z) - n / z
|
|
Nn-1(z) - n / z
|
|
(6)
|
beginning with M0(z) = sinz / (sinz + i cosz). When the argument z is real, the case of Equations 1 and 2, the function Mn(z) is bounded. When the numerical index n is large enough, Mn(z) sharply decreases with increasing n. This property of the function Mn(z) determines the number of terms in Mie scattering calculation. Similar to the formula given by Wiscombe WJ 1980 cited in Mie theory: Number of terms in Mie series, the number of terms, nmax, can be determined depending on the required precision of the result:
|
nmax = {
|
| x1 + 7.5x10.34 + 2 for Mnmax ≤ 10-18
|
| x1 + 6x11/3 + 2 for Mnmax ≤ 10-12
|
| x1 + 4x11/3 + 2 for Mnmax ≤ 10-6
|
| x1 + 4.88x10.31 for Mnmax ≤ 10-5
|
|
(7)
|
The function Qn(z1, z2) can be calculated with the upward recursion:
|
Qn(z1, z2) = Qn-1(z1, z2) +
|
|
[Ln(z1) + n / z1] [Nn(z2) + n / z2]
|
|
[Ln(z2) + n / z2] [Nn(z1) + n / z1]
|
|
(8)
|
by using Q0(z1, z2) = (1 + i cot z1) / (1 + i cot z2). While | z1 | > | z2 | and z1 and z2 are real or complex, the recurrence relation is stable and Qn(z1, z2) converges with the increasing numerical index n (Yang Wen 2003). Numerical results indicate that the method is stable and accurate even when the particle size parameter x1 is as large as 100000 (Liu Lei et al 2007). See also the program code.
