一种改进的mie散射级数算法

一种改进的mie散射级数算法

mie散射理论的计算方法mie色散理论在光散射法中的测量颗粒大小方面发挥着极其重要的作用。这是激光测量的一项重要技术。Mie散射理论中Mie散射级数的计算很重要,它牵扯复杂的贝塞尔无穷级数,计算非常麻烦。几十年来,国内外学者陆续发表许多有关Mie散射理论的数值计算方法,特别是近几年随着计算机技术的迅速发展,Mie散射理论大大提高了计算速度和准确度。然而,Mie散射理论的计算方法还不够完善,各个学者的计算结果也有很大差异,所以Mie散射理论的计算仍然在继续探讨中。笔者在前人的基础上通过Matlab进行了相关计算。1mie的振幅函数Mie散射理论是麦克斯韦方程对处在均匀介质中的均匀颗粒在平面单色波照射下的严格数学解。由Mie散射知道,距离散射体r处p点的散射光强为Isca=I0gλ28π2r2gI(θ‚φ)(1)I(θ‚φ)=|S1(θ)|2sin2φ+|S2(θ)|2cos2φ(2)Ιsca=Ι0gλ28π2r2gΙ(θ‚φ)(1)Ι(θ‚φ)=|S1(θ)|2sin2φ+|S2(θ)|2cos2φ(2)式中:λ为光波波长;I0为入射光强;Isca为散射光强;θ为散射角;φ为偏振光的偏振角。S1(θ)=∑n=1∞2n+1n(n+1)[anπn+bnτn](3)S2(θ)=∑n=1∞2n+1n(n+1)[anτn+bnπn](4)S1(θ)=∑n=1∞2n+1n(n+1)[anπn+bnτn](3)S2(θ)=∑n=1∞2n+1n(n+1)[anτn+bnπn](4)式中:S1(θ)和S2(θ)是振幅函数;an和bn是与贝塞尔函数和汉克尔函数有关的函数;πn和τn是连带勒让得函数的函数,仅与散射角θ有关。其中an=φn(α)φ′n(mα)−mφ′n(α)φn(mα)εn(α)φ′n(mα)−mε′n(α)φn(mα)(5)bn=mφn(α)φ′n(mα)−φ′n(α)φn(mα)mεn(α)φ′n(mα)−ε′n(α)φn(mα)(6)an=φn(α)φ′n(mα)-mφ′n(α)φn(mα)εn(α)φ′n(mα)-mε′n(α)φn(mα)(5)bn=mφn(α)φ′n(mα)-φ′n(α)φn(mα)mεn(α)φ′n(mα)-ε′n(α)φn(mα)(6)式中:φn(α)和εn(α)分别是贝塞尔函数和第一类汉克尔函数;φ′n(mα)和ε′n(mα)是φn(α)和εn(α)的导数;α为无因次直径,α=πDλ‚Dα=πDλ‚D为颗粒的实际直径;λ是入射光的波长;m是散射颗粒相对于周围介质的折射率,它是一个复数,虚部是颗粒对光的吸收的量化。由以上公式可见,Mie散射计算的关键是振幅函数S1(θ)和S2(θ),它们是一个无穷求和的过程,理论上无法计算。求解振幅函数的关键是计算an和bn,所以Mie散射的计算难点是求解an和bn。2n函数中散射角的材料.通过以上分析可知,Mie散射计算的核心是求解an和bn,我们编制程序也是围绕它进行编写。在an和bn的表达式中φn(α),φ′n(α),εn(α)和ε′n(α)满足下列递推关系:φn(α)=2n−1αφn−1(α)−φφ−2(α)(7)φ′n(α)=−nαφn(α)+φn−1(α)(8)εn(α)=2n−1αεn−1(α)−εn−2(α)(9)ε′n(α)=−nαεn(α)+εn−1(α)(10)φn(α)=2n-1αφn-1(α)-φφ-2(α)(7)φ′n(α)=-nαφn(α)+φn-1(α)(8)εn(α)=2n-1αεn-1(α)-εn-2(α)(9)ε′n(α)=-nαεn(α)+εn-1(α)(10)这些函数的初始值为φ-1(α)=cosα(11)φ0(α)=sinα(12)ε-1(α)=cosα-isinα(13)ε0(α)=sinα+icosα(14)与散射角有关的πn和τn满足下列递推公式:τn=πncosθ-π′nsin2θ(15)πn=2n−1n−1πn−1cosθ−nn−1πn−2(16)π′n=(2n−1)πn−1+π′n−2(17)π0=0(18)π0=0(19)π′0=π′1=0(20)πn=2n-1n-1πn-1cosθ-nn-1πn-2(16)π′n=(2n-1)πn-1+π′n-2(17)π0=0(18)π0=0(19)π′0=π′1=0(20)有了这些递推公式可以很方便地通过计算机程序求解。但是对于n的大小,因为计算机不可能计算无穷个数据,所以n在计算之前就要被确定。3计算公式与计算程序结果对比关于n的估计可按Wiscombe给出的经验公式得出。本文将an适当处理为an=φn(α)Dn−mφ′n(α)εn(α)Dn−mε′n(α)(21)Dn=φ′n(mα)φn(mα)(22)an=φn(α)Dn-mφ′n(α)εn(α)Dn-mε′n(α)(21)Dn=φ′n(mα)φn(mα)(22)这样,对an的求解变为对φn(α),φ′n(α),εn(α),ε′n(α)和Dn的求解。计算机计算程序流程如图1所示。下面图2～图7中纵坐标是an的模(即Mie散射级数中的1个,是无量纲的),横坐标是nstop。图2～图7分别是在α=4,m=0.75针对非吸收性颗粒和α=100,m=1.2-1.2i针对吸收性颗粒的条件下,按照前向递推,后向递推和连分式法计算得到的级数模。表1是本程序计算的消光系数kext与Wiscombe(MIEVO)和Du(MIECPP)计算结果的比较。表1表明计算结果基本一致,说明本文算法正确。从以上图中可以看出|an|的震荡是很剧烈的,特别是在α较大的情况下(α=100),n的选择是比较正确的,每幅图在接近nstop时|an|衰减得很快。然而,3种方法的计算结果不尽相同,在α=4,m=0.75时后向递推和连分式结果很接近,前向递推结果整体偏小,但变化趋势与其他方法大致相同;在α=100,m=1.2时前向递推和连分式结果相似,但后向递推结果和它们差别较大,变化趋势却大致相同。4后向递推的算法利用Matlab编写的程序计算不同颗粒直径及其相对折射率下的Mie散射级数可以看到,3种方法的计算结果不尽相同,但是Mie散射级数震荡很剧烈。通过计算多种α和m发现,在它们较小的情况下前向递推还是比较稳定的,随着α和m的增大,特别是α的增大前向递推计算结果极易溢出。单纯的后向递推法由于选择的初始值很难保证