Analytical algorithm to determine the spherical particle size distribution from spectral absorption measurements

Optics Communications 355 (2015) 391–397

Contents lists available at ScienceDirect

Optics Communications journal homepage: www.elsevier.com/locate/optcom

Analytical algorithm to determine the spherical particle size distribution from spectral absorption measurements Jian-Qi Zhao n, Jiangnan Li State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China Canadian Centre for Climate Modelling and Analysis, University of Victoria, Victoria, British Columbia, Canada

art ic l e i nf o

a b s t r a c t

Article history: Received 13 May 2015 Received in revised form 28 June 2015 Accepted 29 June 2015

A modified anomalous diffraction theory (MADT) by including the effects of reflection and refraction is proposed. With respect to MADT, an analytical technique for retrieval of spherical particle size distribution (PSD), based on absorption was developed. Also, an MADT transform pair between the size distribution and the absorption spectrum was constructed. This provides a new tool to study the related particle optical properties. By Gaver–Stehfest inversion method, the derived complex-inversion formula is finally converted into the new real-inversion formula so that the measured absorption data can be applied directly. The inversion experiments show that the use of extended precision instead of double precision arithmetic can produce more reliable results at the expense of CPU time. The effects of complex refractive index on retrieval of PSD were examined. Also it was found that an appropriate reduction of the truncation number with the smoothing technique improved the anti-noise ability for the algorithm. & 2015 Elsevier B.V. All rights reserved.

Keywords: Inverse problem Absorption spectrum Particle size distribution

1. Introduction One important characteristic of many-particle systems is the particle size distribution (PSD). The retrieval of PSD from light scattering measurements is a fundamental work, which has many applications in various fields of science and engineering, such as astrophysics, colloidal chemistry, biophysics, atmospheric optics, oceanophysics, and remote sensing, etc. Over the past several decades, there were many inversion schemes proposed for PSD. These schemes can be roughly classified into two categories: the numerical method e.g., [1–9], and the analytical method [10–16]. Recently, the integral transform technique (analytical-type method) has been developed to retrieve PSD of non-spherical particles based on the extinction and absorption spectrum e.g., [17–20]. In this method, the integral kernel function is expressed in a suitable closed-form by the anomalous diffraction theory (hereafter referred to as ADT) [21]. Thus the PSD for various types of particle shape can be retrieved analytically. The main advantage of the integral transform technique lies in its simplicity and practicality. In particular, this method can provide a global and quick view of the inversion results. However, the application of ADT requires n

Corresponding author. E-mail address: [email protected] (J.-Q. Zhao).

http://dx.doi.org/10.1016/j.optcom.2015.06.071 0030-4018/& 2015 Elsevier B.V. All rights reserved.

that the particle refractive index is close to unity. This constraint implies that the reflection and refraction are negligible as a ray passes through a particle, and so can bring the convenience in mathematical treatment, but will result in a loss of accuracy. The existing methods for calculating absorption based on modification to ADT [22–26] generally can not be applied to the analytical inversion. This paper is devoted to research on improving the accuracy of ADT through the inclusion of the effects of reflection and refraction, and applying the result to the analytical retrieval of spherical particle PSD from spectral absorption measurements. In the following section, a new modified ADT (MADT) is used to compute the absorption efficiency of spherical particles. In Section 3, the MADT integral transform pair between the PSD and the absorption spectrum is constructed, which allows the retrieval of PSD analytically. In addition the real-variable inverse formula is provided. In Section 4, the numerical simulations are implemented to test the performance of the suggested method. In Section 5, we conclude with a brief discussion.

2. Inversion problem and absorption approximation The spectral absorption coefficient βabs (k ) at wavenumber k for light scattering by a dilute suspension of polydisperse spherical

392

J.-Q. Zhao, J. Li / Optics Communications 355 (2015) 391–397

particles of similar optical properties is given by the relationship

βabs(k, m) =

∫0

∞

πa2Q abs(k, a, m)f (a)da

(1)

where a denotes the particle radius, m = mr + imi is the complex refractive index of the particle relative to the medium with mi ≥ 0, f (a) is the PSD of particle system with a radius at a , which needs to be determined, Q abs (k , a, m) is the absorption efficiency of a single particle, representing the fraction of the incident bean absorbed per unit cross-sectional area of a particle. In Eq. (1), the integral equation consists in retrieving unknown f (a) from the measured βabs (k , m) with a series of k . An efficient treatment of the absorption efficiency required in the integral equation is particularly important for the retrieval of PSD. In principle, given k ,a , and m, Q abs (k , a, m) can be calculated exactly on the basis of the Mie scattering algorithm [27–29]. However, the mathematical expression of Q abs (k , a, m) involves the summation of an infinite series of complicated terms including Bessel functions, which hinders a direct analytical inversion of the integral equation. The ADT is an alternate approximation to Mie theory for light scattering by spherical particles in a homogeneous medium, which can be used to build a closed-form expression of Q abs for large optically soft particles. In ADT, the absorption efficiency of a sphere (with radius a and complex refractive index m) is an integral over all local angles of incidence θ at the sphere’s front surface [21],

Q abs = 2

∫0

π /2

⎡⎣1 − exp(−4m ka cosθ )⎤⎦ sinθ cosθ dθ i

8mi ka , 3

Q abs = 1,

(3)

as ka → ∞

(7)

where Im indicates the imaginary part of a complex quantity; r0 is independent of sizes and stands part of the incident energy reflected by a sphere, which can be computed from

(|R1|2 + |R2|2 )cosθ sinθ dθ

R2 =

⎝ mr + 1 ⎠

(m

r

4

−1

) (m

r

2

)

+1

)

ln(mr )

3mr 7 − 7mr 6 − 13mr 5 − 9mr 4 + 7mr 3 − 3mr 2 − mr − 1

+

(

)(

)

3 m r 4 − 1 m r 2 + 1 (m r + 1 )

(10)

Typically, 1 − r0 changes very little with wavenumber, so it is regarded as a constant in this paper. In addition, 1 − r0 → 1, as m → 1; and 1 − r0 → 0, as |m| → ∞. Therefore, for a sufficiently small (compared to the wavelength) particle, ADT absorption formula is consistent with the Rayleigh limit, but for a sufficiently large (compared to the wavelength) particle, ADT absorption formula does not approaches the geometrical optics limit. The physical significance of Eq. (7) is that all the light which penetrates into a sufficiently large absorbing sphere will be absorbed and none of the unreflected light will be transmitted. Referring to the geometrical optics approximation [30,32], the ADT absorption scheme of Eq. (2) may be revised as

∫0 (1 − r0){1 − exp⎡⎣−(4mika cos θ )γθ(1 − r0)−1⎤⎦} π /2

(11)

sinθ cosθ dθ

where γθ represents the geometrical path length ratio of the refracted ray to the non-refracted ray. For simplicity we consider that mi ≪ mr . This condition holds for most substances in the optical range of the spectrum. Under this assumption, the angle of refraction is approximately real, and then γθ can be approximately expressed by

1 − mr −2 sin2 θ

γθ =

(12)

cos θ

From Eq. (11) the revised formula of Q abs is found to be

× exp ( − t 1 − mr −2 )]}

lim Q abs (k, a, m) = 1 − r0

ka →∞

cosθ − m cosθ′ , cosθ + m cosθ′

)

2

(5)

(6)

R1 =

+1

(

8m r 4 m r 4 + 1

{1 + 2mr 2t −2 [(1 + t ) exp ( − t ) − (1 + t 1 − mr −2 )

⎛ m2 − 1 ⎞ ⎟ lim Q abs (k, a, m) = 4kaIm ⎜ ka → 0 ⎝ m2 + 2 ⎠

π /2

2

r

(4)

However, the correctly asymptotic limits, which have universal acceptance based on the Rayleigh theory and the geometrical optics, are given by [21,22,30–32]

∫0

r

3

2

) ln⎛⎜ m −1 ⎞⎟ +

Q abs = (1 − r0 )

as ka → 0

r0 =

(m

(2)

where w = 4mi ka . Eq. (2) illustrates that the absorption by conventional ADT is dependent on the imaginary part of the refractive index; the real part of the refractive index provides no contribution. However, this is usually not accurate. Let us examine the asymptotic behavior of the ADT absorption formula. From Eq. (3), it can be found that

Q abs =

r0 =

(

m r 22 m r 2 − 1

Q abs = 2

which can be evaluated analytically and its result is given by

exp ( − w ) 1 − exp ( − w ) Q abs = 1 + 2 −2 w w2

Eqs. (6) and (7) are comparable with the results from Mie theory. In general, m is not real but complex, therefore r0 can not be evaluated analytically. However, in the case of mi ≪ mr , r0 can be approximated by ignoring mi as shown in references [31,33].

m cosθ − cosθ′ m cosθ + cosθ′

(8)

(13)

where

t=

4mi ka 1 − r0

(14)

However, Eq. (13) is not convenient for inverse-problem solution. The improved formula of Q abs has to be as simple as possible. It is found that the θ dependence of γθ is a second-order effect, and we can simplify the problem by making γθ as an equivalent constant γ in Eq. (11). This simplification will produce a new modified ADT denoted as MADT. The equivalent constant γ can be determined in the small particle limit as

∫0

π /2

(γ cosθ )sinθ cosθ dθ =

∫0

π /2

(γθ cosθ )sinθ cosθ dθ

(15)

As a result

(9)

where θ and θ′ satisfy Snell’s law of refraction: sin θ = m sin θ′; R1 and R2 are the Fresnel coefficients for light polarized parallel and perpendicular to the plane of incidence. Numerical evaluations of

3/2 γ = mr 2 ⎡⎣1 − (1 − mr −2) ⎤⎦

(16)

Replacing γθ with γ and performing integral operation in Eq. (11), the new formula of Q abs can be represented in the form

J.-Q. Zhao, J. Li / Optics Communications 355 (2015) 391–397 0.008 0.007

in evaluation of Q abs by using MADT instead of ADT, especially for very small and very large size parameters. Hence in the next section we will use the MADT absorption formula as a starting point for the analytical inversion.

Mie theory MADT ADT

0.006

393

Qabs

0.005

3. MADT transform and real inversion formula 0.004

For mathematical convenience, the modified effective wavenumber κ in Eq. (17) will be regarded as an independent variable. Substituting Eq. (17) into Eq. (1), the absorption coefficient of spherical particles becomes

0.003 0.002

m=1.33+10-6 i

0.001 0.000

0

500

1000

1500

βabs(κ ) = π (1 − r0) ×

2000

ka Fig. 1. Comparison of the modified anomalous diffraction theory (MADT) with the Mie theory and the anomalous diffraction theory (ADT) in the evaluations of Q abs with refractive index m = 1.33 + 10−6i .

⎡ exp(−κa) 1−exp(−κa) ⎤ Q abs(κa) = (1−r0)⎢1 + 2 ⎥ −2 κa ⎣ ⎦ (κa)2

κ 2⎡⎣βad (κ ) − βad (∞)⎤⎦ = 2π (1 − r0)

∫0

−2π (1 − r0)

∞

(κa + 1)exp( − κa)f (a)da

∫0

∞

f (a)da

lim Q abs =

(23)

∞

0

a2f (a) da is the high

frequency limit of βabs (κ ). Taking a derivative with respect to κ on the both sides of Eq. (23) leads to

(1 − r0 )

ka → 0

(22)

where the dependence of βabs on k and m is implicit. After simple algebra, we can get

κ→∞

(18)

8mi ka 2 ⎡ mr ⎣1 − (1 − mr −2)3/2⎤⎦ 3

d −1 κ 2⎡⎣βabs(κ ) − βabs(∞)⎤⎦ = 2π (1 − r0)κ dκ

{

(20)

ka →∞

lim Q abs = 0

(21)

|m |→∞

These asymptotic limits coincide with the well-known formulas [21,32]. Moreover, as m → 1, Eq. (17) reduces to traditional ADT scheme Eq. (3). It is therefore expected that Eq. (17) can give a reasonable solution for any intermediate values of ka and m. To confirm this we conducted numerical calculations of Q abs . Figs. 1 and 2 give some examples on comparisons of MADT, ADT, and Mie theory, which show that there is substantial improvement

Mie theory MADT ADT

1.2

}

∫0

∞

a2f (a)exp(−κa)da

(24)

Thus

(19)

lim Q abs = 1 − r0

1.4

d κ 2⎡⎣βabs(κ ) − βabs(∞)⎤⎦ = − 2π (1 − r0)κ dκ

{

From Eq. (17) we can obtain the following results:

1.0

Qabs

⎡ exp(−κa) 1 − exp(−κa) ⎤ −2 a2⎢1 + 2 ⎥f (a)da κa ⎣ ⎦ (κa)2

where βabs (∞) = lim βabs (κ ) = π (1 − r0 ) ∫

4mi mr 2 ⎡⎣1 − (1 − mr −2)3/2⎤⎦ k

0.8

} ∫

0

∞

a2f (a)exp(−κa)da

(25)

By the inverse Laplace transform, we can obtain finally an expression of PSD

f (a) =

−1 4π 2i(1 − r0)a2

∫γ κddκ {κ 2⎡⎣βabs(κ ) − βabs(∞)⎤⎦}exp(κa)dκ

(26)

where the integral is usually taken along the so-called Bromwich contour γ running vertically from σ − i∞ to σ + i∞ for a large real σ , and lying to the right of all the singularities of the integrand. The analytical inverse formula can yield the accurate PSD if the correct limit behavior of the absorption spectrum and refractive index are used. A constant refractive index is permissive in spectral inversion. However, in practice the refractive index could be predetermined inadequately for various reasons. Then to what extent will the difference in the retrieval result be if an inaccurate or unmatched m is employed? From Eq. (26), it can be found that the retrieved PSD is associated with the true PSD by the scaling relationship

f (a) = s f ftrue (sa a)

0.6

(27)

with

0.4

sf = sa2

0.2 0.0

∞

(17)

where

κ=

∫0

m=1.33+10-2i 0

500

1000

1500

2000

ka Fig. 2. Comparison of the modified anomalous diffraction theory (MADT) with the Mie theory and the anomalous diffraction theory (ADT) in the evaluations of Q abs with refractive index m = 1.33 + 10−2i .

(1 − r0)true , (1 − r0)

sa =

κ κtrue

(28)

where the subscript “true” means the accurate result or accurate parameter. Analytical The analytical expression of Eq. (27) is consistent with our previous findings [19,20], and it provides a quick view of the sensitivity of the retrieval method to the inversion parameter m, which can be used to cross-reference with other numerical results [34]. Also Eq. (27) implies that the retrieved PSD and true PSD belong to one family of functions.

394

J.-Q. Zhao, J. Li / Optics Communications 355 (2015) 391–397

Hence an unmatched input of m usually does not influence the type of mode of PSD in a retrieval process. Moreover, from Eq. (17), it is easy to obtain

−

1 d 2 {κ [Q abs (κa) − Q abs (∞)]} = exp ( − κa) 2 (1 − r0 ) κ dκ

(29)

inversion algorithm was programmed in FORTRAN and run on a PC with double precision (REAL*8) and extended precision (REAL*16) floating point variables, respectively. The input signals of spectral absorption is simulated in terms of MADT with complex refractive index m = 1.33 + 10−3i and the size distribution of ftrue (a). 2

where Q abs (∞) = lim Q abs (κa) = 1 − r0 . κa →∞

Defining

a

functional

ftrue (a) = f0

n= 1

F [ψ (ξ ) , ξ ] as

(30)

where ψ (∞) = lim ψ (ξ ), ψ (ξ ) is an arbitrary differentiable function ξ→∞

of variable ξ , then Eqs. (22) and (26) can be rearranged as a MADT transform pair

βabs(κ ) =

∫0

πa2f (a) =

∞

(31)

∫γ F[Q abs(−κa), − κa]F ⎡⎣βabs(κ ), κ⎤⎦dκ

(32)

This transform pair reveals the connection between the size distribution and the absorption spectrum, which is a new representation of the solution to the inverse problem. In practice, the absorption coefficient βabs (κ )can only be measured on the positive real axis of a complex plane; the formal solution to the inverse problem shown in Eq. (26) can not be directly applied since it invokes complex contour integrals. Therefore, it is desirable to obtain a real inversion formula which works with only the values of βabs (κ )on real axis (0, ∞). Fortunately, in computational mathematics, the real-inverse problem of the Laplace transform has been resolved, and several schemes have been suggested [35–42]. Among them, the Gaver– Stehfest inversion technique is remarkable, as it can produce stable and accurate results with fast convergence for most cases [43]. By the Gaver–Stehfest inversion

f (a) =

ln2 a

N

⎛ ln2 ⎞ ⎟ a ⎠

∑ vsF ⎜⎝s s= 1

min (s, N /2)

vs = ( − 1) s + N /2

∑ j =[(s + 1) /2]

(33) j N /2 (2j )! (N /2 − j )!j!(j − 1)!(s − j )!(2j − s )!

(34)

where F represents the Laplace transform of f ; the number of expansion terms N is referred to as the “Stehfest number” which must be even and it can be optimized by trial and error; [(s + 1) /2] is the greatest integer less than or equal to (s + 1) /2. The weighting coefficients of vs alternate in sign and grow rapidly. As shown in Eq. (34), the Gaver–Stehfest algorithm is linear to values of F and involves only real algebra. By the Gaver–Stehfest inversion technique, the derived complex-inversion formula of Eq. (26) can be rewritten as N

⎤ ⎡1 d ln2 κ 2⎡⎣βabs(κ ) − βabs(∞)⎤⎦ ⎥ f (a ) = − ∑ vs⎢⎣ 3 ⎦ κ κ d 2π (1 − r0 )a s = 1

{

}

(un a)cn exp (−un a)

(36)

p1 = 0.6,p2 = 1 − p1, u1 = 8μm−1,c1 = 8,u2 = 12 μm−1, and c2 = 30. This typical bimodal gamma distribution is usually appropriate for atmospheric particles. The accuracy of the algorithm is measured by the maximum relative error (MRE) defined as

MRE =

Q abs(κa)⎡⎣πa2f (a)⎤⎦da

1 2π i

pn un Γ (cn + 1)

where the adjustable parameters were chosen as f0 = 1cm−3,

d 2 1 {ξ [ψ (ξ ) − ψ (∞)]} 2 (1 − r0 ) ξ dξ

κ = s ln2 a

(35)

This is the general inverse formula for absorbing particles. It permits the computation of the particle size distribution f (a) from the absorption function βabs (κ ) − βabs (∞) for particles of a series of κ.

max

To assess the performance of the new inversion algorithm, we implemented simulations of the retrieval calculations. The

ftrue (a)

(37)

10 10

double precision extended precision

10 10 10 10 10 10 10 10 10 10 4

4. Retrieval simulant test

f (a) − ftrue (a)

a ∈(0,4μm)

where (0, 4 μm) is a significant size interval in which the properties of the function curve of ftrue (a) are shown adequately. Fig. 3 shows the results of MRE for retrieved PSD. It is shown that the Stehfest number N has a strong effect on the accuracy of PSD retrieval. For the extended precision arithmetic, the bigger Stehfest number generally corresponds to the higher accuracy. For the double precision arithmetic, the error of inversion will increase with increase of Stehfest number, when N exceeds a critical point. It is also found that the use of extended precision instead of double precision variables usually provides a more reliable result and requires only a negligibly small amount of extra memory, but at the expense of a significant increase in CPU time. In the following calculations, only the double precision arithmetic will be employed, and the threshold value of N will be regarded as an optimal value. Fig. 4 shows the effects of the complex refractive index on the retrieval of PSD. It is found that if the m is set properly the retrieved PSD agrees very well with the true PSD; on the contrary, if the m is set improperly the retrieved PSD differs significantly from the true PSD. Generally the retrieved PSDs are similar in shape and peak position for different m due to the scaling relation in Eq. (27). However, we also found that the inversion failed when sa exceeds around 10. With the intention of examining the noise tolerance level of an inversion algorithm, a random noise is added to the input signals

MRE

F [ψ (ξ ), ξ ] = −

∑

8

12

16

20

N

24

28

32

36

40

Fig. 3. The maximum relative error (MRE) of retrieved PSD vs. Stehfest number N in size interval (0, 4 μm) . The inversions are calculated with double precision (REALn8) and extended precision (REALn16) floating point variables, respectively.

J.-Q. Zhao, J. Li / Optics Communications 355 (2015) 391–397

395

2

Without smoothing With Smoothing

25

True Retrieved: m=1.33+1.0 ×10-3 i m=1.40+1.0×10-3 i m=1.20+1.0×10-3 i m=1.33+1.5×10-3 i m=1.33+0.6×10-3 i

20

15

N

PSD (µm-1cm-3)

3

10

1

5

0

0

0

1

2

a (µ m)

3

10

4

Fig. 4. Effects of complex refractive index m on PSD retrieval. The inversion algorithm is conducted with double precision arithmetic and N = 24 .

10

10

10

10

Noise level

10

10

10

10

10

Fig. 5. The optimized N for a given noise level. The inversion algorithm is conducted with double precision arithmetic. The “smoothing” in the legend means that the Savitzky–Golay smoothing technique is applied to both the input of spectral absorption and the output of retrieved result.

of the spectral absorption as

10

(38)

where R (0, 1) represents a normally distributed random number with mean of 0 and a standard deviation of 1, ε is a constant, indicating the noise level, β0 (κ ) is the undisturbed absorption data. On the other hand, we employed the Savitzky–Golay data smoothing method [44,45] to suppress fluctuations and increase noise tolerance for retrieval. Figs. 5 and 6 demonstrate the effects of noise level on the optimized N and the MRE of retrieved PSD. It is observed that: (i) with an increase in noise level, the optimized N decreases but the MRE increases; (ii) the optimized N with the application of the smoothing technique is greater than or equal to that without the application of the smoothing technique, and the MRE with the application of the smoothing technique is undoubtedly smaller than the MRE without the application of the smoothing technique. Therefore, an appropriate reduction of N with the smoothing technique indeed helps to improve the anti-noise ability of the algorithm. Figs. 7 and 8 provide the examples showing the reproducibility of PSD for different noise levels. It confirms that the best effect of the smoothing procedure is for ε = 10−8 and N = 16 (Fig. 8) where the retrieval errors are reduced by almost an order of magnitude.

Without smoothing With Smoothing 10

MRE

β (κ ) = β0 (κ )[1 + ε⋅R (0, 1)]

10

10

10

10 10

10

10

10

10

10

10

10

10

10

Noise level Fig. 6. The maximum relative error (MRE) of retrieved PSD for a given noise level. The inversion algorithm is conducted with double precision arithmetic and the optimized N .

0.8

True Retrieved:

0.7

ε =10-8, N=16, without smoothing ε=10-8, N=16, with somoothing

5. Summary and conclusions

PSD (µ m-1cm-3)

0.6

In this work we have proposed a modified ADT by including the effects of reflection and refraction. Compared to conventional ADT the large improvement in evaluating Q abs is shown especially for small and large size parameters. In the context of MADT, the analytical technique for retrieval of spherical PSD based on absorption data is therefore developed. Also an MADT transform pair between the size distribution and the absorption spectrum is constructed. This provides a new tool to study the related particle optical properties. By the Gaver–Stehfest inversion technique, the derived complex-inversion formula is converted into the new real-inversion formula so that the measured absorption data can be applied directly. The retrieval experiments show that the use of the extended precision instead of double precision arithmetic yields more reliable results at the expense of CPU time. Generally, the Stehfest number N (indicating truncation number) plays a role

10

0.5 0.4 0.3 0.2 0.1 0.0

0

1

2

a (µ m)

3

4

Fig. 7. Reproducibility of PSD for double precision arithmetic and noise level ε = 10−8.

396

J.-Q. Zhao, J. Li / Optics Communications 355 (2015) 391–397

2.0

True Retrieved:

ε =10-7, N=16, without smoothing ε =10-7, N=16, with somoothing

PSD (µ m-1cm-3)

1.5

1.0

0.5

0.0

0

1

2

a (µm)

3

4

Fig. 8. Reproducibility of PSD for double precision arithmetic and noise level ε = 10−7 .

similar to a regularization parameter. For the extended precision arithmetic, the bigger the Stehfest number almost means the higher the accuracy. For the double precision arithmetic, however, it shows that: (i) the Stehfest number have an optimal value, at which the highest accuracy is achieved; (ii) with the increase of the noise level, the optimized N decreases; (iii) the optimized N with the application of the smoothing technique is not less than that without the application of the smoothing technique. Also it is confirmed that the MRE increases with the increase of noise level, and the MRE with application of the smoothing technique is undoubtedly smaller than the MRE without application of the smoothing technique. Therefore, an appropriate reduction ofN with the smoothing technique indeed helps to improve the antinoise ability for the algorithm. The effects of complex refractive index m on retrieval of PSD are examined. It is found that if the complex refractive index m is set properly the retrieved PSD agrees very well with the true PSD; on the contrary, if the m is set improperly the retrieved PSD differs significantly from the true PSD. Generally the retrieved PSDs are similar in curve shape and peak position for different values of m. This similarity relation also presents a quick view of the sensitivity of a retrieval method to complex refractive index.

Acknowledgments This research was financially supported by the NSFC (Grant no. 41075015).

References [1] D.L. Phillips, A technique for the numerical solution of certain integral equations of the first kind, J. Assoc. Comput. Mach. 9 (1962) 84–97. [2] S. Twomey, On the numerical solution of Fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature, J. Assoc. Comput. Mach. 10 (1963) 97–101. [3] G. Yamamoto, M. Tanaka, Determination of aerosol size distribution function from spectral attenuation measurements, Appl. Opt. 8 (1969) 447–453. [4] H. Grassl, Determination of aerosol size distributions from spectral attenuation measurements, Appl. Opt. 10 (1971) 2534–2538. [5] G.E. Shaw, Inversion of optical scattering and spectral extinction measurements to recover aerosol size spectra, Appl. Opt. 18 (1979) 988–993. [6] G. Ramachandran, D. Leith, L. Todd, Extraction of spatial aerosol distribution from multispectral light extinction measurements with computed tomography, J. Opt. Soc. Am. A11 (1994) 144–154.

[7] M. Ye, Lu. Y. WangS, T. Hu, Z. Zhu, Y. Xu, Inversion of particle-size distribution from angular light-scattering data with genetic algorithms, Appl. Opt. 38 (1999) 2677–2685. [8] O. Dubovik, M.D. King, A flexible inversion algorithm for retrieval of aerosol optical properties from Sun and sky radiance measurements, J. Geophys. Res.Atmos. 105 (D16) (2000) 20673–20696 2000. [9] Y. Wang, S. Fan, X. Feng, G. Yan, Y. Guan, Regularized inversion method for retrieval of aerosol particle size distribution function in W1,2 space, Appl. Opt. 45 (2006) 7456–7467. [10] K.S. Shifrin, A.Y. Perelman, The determination of the spectrum of particles in a dispersed system from data on its transparency. I. The fundamental equation for the determination of the spectrum of the particles, Opt. Spectrosc. 15 (1963) 285–289. [11] M.A. Box, B.H. McKellar, Analytic inversion of multispectral extinction data in the anomalous diffraction approximation, Opt. Lett. 3 (1978) 91–93. [12] K.S. Shifrin, A.Y. Perelman, V.M. Volgin, Calculations of particle-radius distribution density from the integral characteristics of the spectral attenuation coefficient, Opt. Spectrosc. 51 (1981) 534–538. [13] A.L. Fymat, Analytical inversions in remote sensing of particle size distributions. 1. Multispectral extinctions in the anomalous diffraction approximation, Appl. Opt. 17 (1978) 1675–1676. [14] C.B. Smith, Inversion of the anomalous diffraction approximation for variable complex index of refraction near unity, Appl. Opt. 21 (1982) 3363–3366. [15] J.D. Klett, Anomalous diffraction model for inversion of multispectral extinction data including absorption effects, Appl. Opt. 23 (1984) 4499–4508. [16] J. Wang, F.R. Hallett, Spherical particle size determination by analytical inversion of the UV–visible–NIR extinction spectrum, Appl. Opt. 35 (1996) 193–197. [17] J.Q. Zhao, A new method for the determination of the spheroidal nonabsorbing particle size distribution, J. Grad. School Chin. Acad. Sci. 27 (3) (2010) 323–330. [18] J.Q. Zhao, Simple technique to evaluate effects of nonsphericity and orientation on particle size distribution retrieval, Int. J. Phys. Sci. 6 (31) (2011) 7100–7105. [19] J.Q. Zhao, J. Li, Analytical transform techniques to retrieve non-spherical particle size distribution, J. Quant. Spectrosc. Radiat. Transf. 129 (2013) 287–297. [20] J.Q. Zhao, F. Zhang, L. Li, Analytical inversion of the absorption spectrum to determine non-spherical particle size distribution, J. Quant. Spectrosc. Radiat. Transf. 149 (2014) 128–137. [21] H.C. van de Hulst, Light Scattering by Small Particles, Wiley, New York, 1957. [22] P.H. Levine, Absorption efficiencies for large spherical particles: a new approximation, Appl. Opt. 17 (24) (1978) 3861. [23] S.A. Ackerman, G.L. Stephens, The absorption of solar radiation by cloud droplets: an application of anomalous diffraction theory, J. Atmos. Sci. 44 (1987) 1574–1588. [24] D.L. Mitchell, Parameterization of the Mie extinction and absorption coefficients for water clouds, J. Atmos. Sci. 57 (2000) 1311–1326. [25] M. Xu, M. Lax, R.R. Alfano, Light anomalous diffraction using geometrical path statistics of rays and Gaussian ray approximation, Opt Lett. 28 (2003) 179–181. [26] P. Yang, Z. Zhang, B.A. Baum, H.L. Huang, Y. Hu, A new look at anomalous diffraction theory (ADT): algorithm in cumulative projected-area distribution domain and modified ADT, J. Quant. Spectrosc. Radiat. Transf. 89 (2004) 421–442. [27] G. Mie, Beitrage zur optik trüber medien speziell kolloidaler metallösungen (Contributions to the optics of turbid media, particulary of colloidal metal solutions), Ann. Phys. 25 (1908) 377–445. [28] W.J. Wiscombe, Improved Mie scattering algorithms, Appl. Opt 19 (9) (1980) 1505–1509. [29] M.I. Mishchenko, L.D. Travis, A.A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles, Cambridge University Press, Cambridge, 2002. [30] H.M. Nussenzveig, W.J. Wiscombe, Efficiency factors in Mie scattering, Phys. Rev. Lett. 45 (1980) 1490–1494. [31] C. Acquista, A. Cohen, J.A. Cooney, J. Wimp, Asymptotic behavior of the efficiencies in Mie scattering, J. Opt. Soc. Am. 70 (8) (1980) 1023–1025. [32] C.F. Bohren, B.M. Herman, Absorption and Scattering of Light by Small Particles, Wiley, New York, 1983. [33] A.A. Kokhanovsky, E.P. Zege, Local optical parameters of spherical polydispersions: simple approximations, App. Opt. 34 (24) (1995) 5513–5519. [34] Y. Liu, W.P. Arnott, J. Hallet, Particle size distribution retrieval from multispectral optical depth: influences of particle nonsphericity and refractive index, J. Geophys. Res. 104 (D24) (1999) 31753–31762. [35] E.L. Post, Generalized differentiation, Trans. Amer. Math. Soc. 32 (1930) 723–781. [36] D.V. Widder, The inversion of the Laplace integral and the related moment problem, Trans. Amer. Math. Soc. 36 (1934) 07–200. [37] G. Doetsch, Handbuch der Laplace, Transformation, Vol. 1, Birkhäuser, Verlag, Basel, 1950. [38] D.P. Gaver Jr., Observing stochastic processes ans approximate transform inversion, Oper. Res. 14 (3) (1966) 444–459. [39] H. Stehfest, Algorithm 368, Numerical Inversion of Laplace Transforms, CACM 13 (1) (1970) 47–49. [40] J. Peng, S.K. Chung, Laplace transforms and generators of semigroups of operators, Proc. Amer. Math. Soc. 126 (8) (1998) 2407–2416. [41] V.K. Tuan, D.T. Duc, A new real inversion formula for the Laplace transform and its convergence rate, Fract. Calc. Appl. Anal. 5 (4) (2002) 387–394. [42] V.K. Tuan, T. Tuan, A real-variable inverse formula for the Laplace transform, Integral Transfor. Spec. Funct. 23 (8) (2012) 551–555.

J.-Q. Zhao, J. Li / Optics Communications 355 (2015) 391–397

[43] P.P. Valko, J. Abate, Comparison of sequence accelerators forthe Gaver method of numerical Laplace transform inversion, Comp. Math. Appl. 48 (3) (2004) 629–636. [44] A. Savitzky, M.J.E. Golay, Smoothing and differentiation of data by simplified least squares procedures, Anal. Chem. 36 (1964) 1627–1639.

397

[45] H.H. Madden, Comments on the Savitzky-Golay convolution method for leastsquares-fit smoothing and differentiation of digital data, Anal. Chem. 50 (9) (1978) 1383–1386.