mie 散射系数的新算法
王少清 任中京　张希明　何芳 江海鹰
(山东建筑材料工业学院，济南　250022)
提要：介绍了一种计算mie散射系数的新方法，给出了计算实例。
关键词：米氏散射，mie系数，mie计算
new algorithm of mie scattering parameters
wang shaoqing ren zhongjing　zhang ximing 　hefang jiang haiying
(shandong institute o f building materials，jinan　250022)
abstract：a new algorithm of calculating mie parameters was introduced. some calculation results done with this algorithm were given.
key words: mie scattering ，mie parameters，mie calculation
1　引言
mie 理论是球形颗粒对单色光的散射场分布的严格解析解[1]，目前在环保、动力、气象、天文、两相流及粉体颗粒尺寸分布测量等工程技术领域中有广泛的应用。利用单一颗粒或颗粒群光散射场的测量数据。可以反推得散射颗粒或颗粒群的很多物理性质，如颗粒的尺寸、颗粒的折射率等[2]。但反推必须事先计算出各种尺寸的颗粒在各种复折射率下的散射场分布数据。1968 年dave[3]zui先发表了完整的mie散射计算方法，以后lentz[4]和wiscomb[5]又针对部分计算提出了新算法。国内也有人发表了他们自己的算法[6、7、8]。但总的看来，这些算法均有各自的局限性。尤其是当颗粒尺寸或折射率的虚部值较大时，往往计算速度过慢或产生溢出和不收敛等现象。本文介绍作者发展的mie 散射新算法。该算法的特点是不受颗粒尺寸及折射率的限制，不会产生溢出和不收敛的现象，且具有较快的计算速度。
2　mie散射系数的计算公式
mie 散射计算的中心问题是计算mie 散射系数an 和bn，其表达式为[9]
其中为颗粒的尺寸参数，定义为α=πd/λ，d为颗粒直径，λ为入射光在颗粒周围介质中的波长，而m为颗粒在周围介质中的相对复折射率，即m= m1 + im2 ( m2< 0) , 式中i 为虚数单位。而ψn( z) 和ξn( z) ( z 既表示α，又表示mα)的表达式为
ψn( z)=（πz/2）jn+1(z)
ξn( z)=ψn( z)+iχn(z)
χn(z)= ? (πz/2)( ? 1)n-1j-(n-1)(z)
= ? (πz/2)nn+1(z)
3 计算中产生溢出的原因
计算mie 散射系数须先计算ψn 和χn。一般采用递推的方法。递推又分为向前递推( 即从n= 0 开始) 与向后递推( 即从n= n 开始至n=0, n 为预先设定值) 。实验表明，向前递推总是快于向后递推。ψn 和χn的初值为
分析以上两式可知，当m2≠0 时, 若颗粒尺寸d很大, 或复折射率的虚部值m2很大，将使得乘积m2d很大，可使两式中的项exp( - m2 α) =exp ( - πm2d/ λ) 的值超过计算机的数据限，从而产生溢出，这是产生溢出的重要原因。另外，在递推过程中，不恰当的算法也可能造成溢出。
4　an和bn的新算法
为解决上述问题，作者提出了新的算法。将an和bn的公式变形如下: 令
其中lnr、lnj分别表示ln(m)的实部与虚部。将（2）式代入（1）式，并用anr、anj和bnr、bnj分别表示的实部与虚部。如此可推得
在上述四个公式中采用比值的形式是非常重要的，这样可避免递推过程中当ai与bi 较大时乘法运算可能产生的溢出。这是本算法的一个重要特点。在以上四式中
由于均为实变量函数，计算不产生溢出。关键是和的算法如何处理，才能保证计算中不产生溢出。在lentz的算法中是采用连分式计算ln的值，其精度的保证是由在大量计算基础上得出一个截断项数n与参数a与m的经验公式而实现的。这样的经验公式，一则有实用上的局限性，再则也会带来截断误差。文献（6）对此经验公式做了改进，但仍陷于a=1~100，m1=1~2，m2=0~1的范围。下面介绍本文作者发展的关于ln的新算法。该算法的特点是不受a及m值的限制，不会产生溢出或不收敛等病态现象，且具有较快的计算速度。令
上面导出的(3) - (20)诸式，构成了mie系数an和bn的完整算法。由于an和bn是从n= 1开始计算，利用初值公式(16) - (20)即可算得任意级数的an和bn的值，故没有舍入误差的问题。从(16)式可见，因为y= m2ɑ≤0, 故无论m2和ɑ取何值均不会产生溢出，再加上(3)的各式中采用了比的形式，又避免了计算过程中的溢出，这就从根本上解决了溢出的问题。
5　计算实例
利用以上算法编制了计算颗粒散射场强度和消光系数的计算机程序。当波长为λ的单位振幅的平面自然光入射颗粒时，颗粒的散射光强为[9]
ξ为计算机在双精度下的zui小数据限。
图1示出散射光强的一组计算实例。其中取m1=1.33，m2=-0.4，λ=0.6328三图分别对应于颗粒的直径为d=0.001,1.0和30μm。而d图为颗粒直径d=100μm时散射花样的局部放大。可见，随着颗粒尺寸的增大，前向散射迅速加强，并且出现了复杂的旁瓣。
实部(a)与虚部(b)的变化情况。可见随着m1和m2的增大，虽然颗粒的尺寸保持不变，但散射也加强，且后向散射随着m1和m2的增大而加强。
图3给出有关消光系数的计算结果，其中a)和b)分别表示消光系数随折射率实部与虚部的变化情况。可见随着颗粒直径的增大，消光系数趋近于2; 折射率的增大，尤其是折射率虚部的增大，使这一趋近变得更快和更明显。另外当折射率的虚部m2= 0时，消光系数随颗粒直径的增大而振荡；但当m2≠0 时，振荡迅速消失。
光强zui大值处所对应的fm(z)的值即由上式确定。由上式也可见：此极限情况下的焦移大小主要由s0/f及na所决定。
参考文献：
[1]m. bor n and e. wolf, principles of optics ( sixt hedition ) , pergaman press ( oxford, new york, paris) , pp611(1980)
[2]peter chylek, v. ramasw amy , a. ashkin and j.m . dziedzic,“simultaneous determination o f refractiv endex and size o f spherical dielectric particles from light scattering data”,appl. opt . ,vo l. 22, no . 15, p2302- 2307(1983)
[3]j. v. dave, repo rt 320 - 3237 ( ibm scientific center , 1968)
[4]w. j. lentz, appl. opt. , 15, 668( 1976)w. j. wiscome, appl. opt, vol. 19, no. 9, pp1505( 1980)
[6]顾冠亮等, 有关光散射物理量的数值计算, 上海机械学院学报, 1984 年4 期pp21
[7]余其铮等, mie 散射算法的改进, 哈尔滨工业大学学报, 1987 年3 期, pp21
[8]郑刚等, m ie 散射的数值计算, 应用激光, vol.12, no . 5, pp220( 1992)
[9]h. c. v an de hulst , lig ht scatter ing by small particles,do ver publica tion lnc. , new yor k, chp. 13( 1981)