j_也W（ojq、；血血卩；（cpW伽cp
h_Hl6,（（X）卩：
伽ex}•壬：
（瓦9?
血0（）
式中：
n（）和n（）分别是贝塞尔函数和第一类汉克尔函数；n（）和n（）是n（）和n（）的导数；为无因次直径，D,D为颗粒的
实际直径；是入射光的波长；m是散射颗粒相对于周围介质的折射率，它是一
个复数，虚部是颗粒对光的吸收的量化。
由以上公式可见，Mie散射计算的关键
是振幅函数§（）和S2（）,它们是一个无穷求和的过程，理论上无法计算。
求解振幅函数的关键是计算an和bn,所以Mie散射的计算难点是求解an和bn。
Mie散射理论的数值计算
通过以上分析可知，Mie散射计算的核心是求解an和bn,我们编制程序也是围绕它进行编写。
在an和bn的表达式中n（）,n（）,n（）和n（）满足下列递推关系：
W如咏i（oO-WHod
LA
0；（«）=-十卬（06+«［（«）
&（0£）=咕1（（X）°?
（0）
c（oi=亠oc&+G-1Cod
这些函数的初始值为
l（Of）=coscc
尿®二sni0（
Ci（0（）=cosft-ismft&（tt）=呂inWlcostX
与散射角有关的n（）和n（）满足下列递推公式
Ti=TtCOsO-TTnSttf0
\：
_1八J7
iT-=-JT-icoit%-f?
F-in-1n-1
7Tw=L—-1.*|f~ui-lH-7Trr-5
m=0
7Ti=0
7T：
=7T‘】=0
有了这些递推公式可以很方便地通过计算机程序求解。
但是对于n的大小,
因为计算机不可能计算无穷个数据，所以n在计算之前就要被确定。
散射理论基础与Matlab实现
若散射体为均匀球体，如图1所示，照射光为线偏振平面波，振幅为E,光强10,沿z轴传播，其电场矢量沿x轴振动。
散射体位于坐标原点0,P为观测点。
散射光方向（0P方向）与照射光方向（z轴）所组成的平面称为散射面，照射光方向至散射光方向之间的夹角B称为散射角，而x轴至0P在xy平面上投影线（0P'）之间的夹角©称为极化角。
观测点与散射体相距r。
根据经典的Mie散射理论,
散射粒子的尺度参数为a=2na/入，其中a为球形粒子的半径，散射粒子相对周
围介质的折射率为m=m1+i*m2。
则散射光垂直于散射面和平行于散射面的两
个分量的振幅函数为
护;f也纲-択叭何巴（“皿丿务=5（久丿貯；r也列-用£；f戊丿y（m盘）
監（QV；（m（xj-'卩；（隨）0^（ma）b"~m^n（a）貯；（mA）•E：
（a）爭,m^）
P：
‘fCO£0丿
Tn二sine
以上式中：
J=需胃仃8詢
0；心丿=^n-1（Z）-卫收（Z）
■W
£；（刃=1（Z）-卫£，刃
Jn+1/2（z）和Yn+1/2（z）分别为半整数阶的第一类，第二类贝塞尔函数
P⑴n（cos9）为一阶n次第一类缔合勒让德函数；Pn（cos9）为第一类勒让德函数。
在数值模拟过程中选取初始下：
巧=cosO
=1
微粒子对光的散射和吸收是电磁波与微粒子相互作用的重要特征，而微粒对电磁辐射的吸收与散射与粒子的线度有密切关系，对于不同线度的粒子必须应用不同的散射理论。
Mie散射理论主要用于从亚微米至微米的尺寸段；在微米以下至纳米的光散射则近似为形式更明晰简单的瑞利散射定律，散射光强烈依赖于光波长入（I〜入-4）；而对大于微米至毫米的大粒子则近似为意义明确的夫朗和费衍射规律了。
Mie散射理论给出了球型粒子在远场条件下的散射场振幅an、bn以及粒子
内部电磁场振幅cn、dn的计算表达式,通常称为Mie散射系数
*肿片（m戈）[xh?
']（x）[-h（n'（x）[mxjn（fnx）]
jrf1
卜_"7X丿『XWJ1-7廿（兀丿『爪力n（丿片龙丿1
n^\jn（rnx）[xh^]（x）J-hnf（x）[mxjn（fnx）]
=就贰（[英讥'（英）1-汨（工）fxj戳（刈]
5一一\r）（1?
zi'7
（1）j-\r■/\J
mx）[xh„（x）]-h„（x）[mxjn（mx）J
TX）]•g讥¥"（X）『艾gX）i
"狀y汎x）{xh计（x）]-jUj（x）[mxy„（mx）]
式中m表示微粒子外部介质的相对折射率,x=ka,a为球的半径,k=2n/入称为波数,卩为相对磁导率，即球的磁导率与介质磁导率的比值,jn（x）和h⑴n（x）分别为第一类虚宗量球Bessel函数和Hankell函数。
散射系数，消光系数及偏振状态下散射相位函数：
g=飞刃2n+Vf|an\2+||2）
an^l
散射截面Cscs（散射率Qsca）＞吸收截面Cabs（吸收率Qabs）、消光截面Cext（消光率Qext）、后向散射截面Cb（后向散射率Qb）以及辐射压力Cpr（辐射压力效率Qpr）。
其表达式如下：
其中i为sca、abs、ext、pr分别表示散射、吸收、消光、辐射压力。
按照能量守恒定律有：
OeV—Q宾$亠Otjfcr—脊TO+迁7裁
Qpr（辐射压力效率的计算公式）
Qpr=-Q=eaV心0>
这些都是无穷级数求和，在实际计算过程中必须取有限项，Bohren和Huffman给出了级数项最大值取舍的标准：
“昨=x+4a1J+2
对于单位振幅入射波经微粒散射后，其散射场振幅的大小与散射角有关，在球坐标系下，远场散射振幅的大小为：
ikr
£>=讥cos^*Sy（cosO）
-ikr
ikr
^7sin
ikr
其中S1和S2为散射辐射电场在垂直及平行于散射面的两个偏振分量。
微球内部场振幅计算公式
颗粒内部电场强度为：
DO
\~*―—1—_1rrbJXrfb
7]/丄i1n~□“A«>1m/
其中M
（1）o1n和N
（1）e1n为矢量波球谐函数，在球坐标系中定义如下：
0
M；i；=
COS1#*・T„j訂FH1X）
sin"・rn（cosO）j„（rmx）
t“小『rmxj*（Fmx）1
cos*p•rwfco$p>
rmx
吸收截面Qabs
具有损耗介质颗粒的吸收截面为:
其中&〃是粒子相对介电常数的虚部，经整理可得:
式中mn、nn为:
„2n（ln+f”十D
实际上由Mie散射理论可知，上式中的积分项为电场强度的平方对角度B、©全空间积分的平均值，即：
于是吸收效率为:
El耳x^dx
式中X=rk=z/m。
当xn1时即瑞利散射情况，颗粒的内部平均场强为常数，其值
9
为：
-_
ImprovedMiescatteringalgorithmsWJ.Wiscombe
Mie计算存在的问题就是如何最有效地构造Mie计算，同时保证准确性和避
免数值的不稳定性和病态。
Mie计算以耗时著称，首先无穷项级数N的求和，例如：
100m的水滴在0.5m的可见光散射情况下，大约需1260项求和。
其次，典型
t.=-：
L+1IRetafl+bn）.
的计算都希望能对一系列半径（如对尺寸分布求积分）、一系列波长（如对太阳光谱求积分）及一系列折射率求和（如通过散射参量反推折射率）。
Qeit.2_
Qa刍£（2»+l）（|an[3+|iJ^x2m
■只£（＜1用£1门+|十L）
+1
4理[n（n4-2）
:
V
r»-.-
2n+t
+
n（n+1）
N+L
必3■E*7—9訶訂从）+如％
n\n+I）
N2n+1・
=L=Mn丁狙川+%声“g几
n-m（n+1）
Ai（jnx）=*；（耐训M耐人
当折射率虚部mIm很大时，用向后循环法求An很不稳定。
而向前递推总是稳定的（但向后递推安全时，总是优先选择，因为其计算速度很快）。
得岀允许向后递推的经验标准：
mlraJC£/（marK
用正确的向前地推与相对应的向后地推做比较，当发现对Q©山Qaen•和g的相对误差超过10-6
时，认为计算失败。
对于一对确定的（x,mRe）,我们采用向后递推寻找第一个循环失败的"1'研究表明：
对于确定的皿只刖，，肌臥工的值随着x的增加很快趋向于一个确定值
xmfm>min（xmira）=/（mRj>*ux
伽二—8+筑.22耐社一0.4474^1.|（^
+0W04mg#-0.000175rp»L»
如果在任意角度下S,、S2的实部和虚部的相对误差超过105时，认为对S,和
S2的向后递推失败。
（而此时，QscaQxt并不受影响，因为当S1，S2的相对误
差达到105时，QScaQext的相对误差总维持在1010以下。
）
=13,78^1*-10,8mr*+3.9.
1200
■■00»000
900
600
700
€00
500
400
300
200
■00
打伽“=16.35mL+对散射强度|£屮+|$屮|和偏正度
（|S2p"|Si|2）/（|S2P+S1p）
连分式算法总结：
Mie散射计算的核心是计算an和bn
収/釦仇（md）-刖X（疣）认（搐a丿
°n和3叽（曲）-（a）叽（ma）
用心丿虹fh但）・虹（罠）认（肋彰
bt\—y1/'zI尸J、f/.
mkn（a）%（ma）-行f、a）tpn（z）
其中®n（a）=aJn（a）,En（a）=aJn（a）+iaYn（a）,Jn和Yn分别是第一和「类贝塞耳函数,a
称为当量直径,a=2nr/入,r是球形颗粒的真实半径，入是入射光的波长，m为折射率
♦wW=
tnU）=（m/2）叫几也（z）+卜1）啦”说@）]
r
—-tL<仇3
DM八訣仏八叽（p）
式中p为函数任一自变量。
贝塞耳函数递推关系式：
tf/Lf.HUMY饶丿f卍Hf口丿
Ca）=叽.、（a）・;（a）=x（a）-
a（X
『Dm『加抚丿/删十h/o7fo丿-也」亠ifa丿
5[Dn（/m+ji/Ct）-n_if厌）
r『血D划（me）+k/a]W刃-C—if加
斤一[mDn（m^）+n/a]^n（a）-Zn-\（a）
Mie散射计算中Jn、Yn、Dn的计算是关键和难点。
对于Dn,我们采用的是Lentz的连分式的算法：
H丁—15加
Dn（ma）=-+T
J（!
fn闯
Lentz证明有如下关系：
J蛊.1fW）_『g]『awdf业…Vi1
J,ma）[az][g…s]
其中,%=（・】）10m+k•S5丿。
我们注意到当时，
'。
所以可以利用上式累积相乘直到满足精度要求。
（可根据精度要求例如10-7来确定所要达到的k值）对于J八Yn的生成本文也采用连分式的算法。
具体方案如下：
令Cn=Jn-1（a）/Jn（a）,根据贝塞耳差积公式：
J訂勿Yz・JzM勿丿=
由以上二式整理得：
n養.、（広）/d）j（Q•a"
Yjg=
JM-1（^7
上式中Cn的计算是采用类似于Dn的连分式的形式，计算中可调用同一函数计算。
若已知初值：
Jo=sinif/口；Jifa丿=sin^/a1•cos^/场；
Yo=-cos^/a；Yi（a）=-cos^/a2-sina/coa；这样就可计算出各级Jn和Yn。
WilliamJ丄entz关于连分式的文章：
*WU）=佃z/2）山几
SnU）=5“2）叫+H）啦®@）]
“世上）左丿”诜（/其中n+%=勺N=mx0
以九仗）为基础，采用贝塞尔函数比值的连分式表示法：
几（刃人儿“⑵，利用此法可产生所有的乩】3），尽管耗时，但能减少存储需求。
同时可通过计算高阶：
：
；VUe值，使用下面的递推公式，从后往前算出其他值。
=纽T-M九I&）e九i0）
不像一般的函数，贝塞尔函数的比值一旦超过可控制的边界，就不再增长，初始的高阶值决定了所有低阶值的准确性，因此，采用新方法计算准确的初始比值是必要的。
/=+1
a24-t
旳+1
5十I
1111
一角+运+石+石+运+•…
Jy-f）__齐“［I11
几⑷=2M+—2®+1）"+沙+加」
1+-2®+3比“十……
处于分母位置的+号表示分母上加上一个特殊的连分式。
类似于上式f中的表示
形式。
定义一种新的符号：
”-［映2,…1=叫+2+£+2+…
/rrtv）=［碍厲.+£
Lentz给出了n阶部分收敛值为：
/a=
1］•t^rl-1J1•"陌］
■-'・
心卅…•卫］
-旳］…Sn-1*…，旳］
.码.♦….旳
例如：
实变量H=1・°，虚数"=乩5:
计算过程:
a„=（-l）n*l2（v+n—n=112,3>.・・
19［a2］=—21［阴卫J=旳+»=—20*94丁36842
（19）（-20J4736842）（22.95226131）（^24.95643131）（26.95993017）
（-21）（22.95238095）（-24.95643154）（26.95993017）
=16.95228198,
米散射学习目前所遇到的困难：
到底怎样的计算结果才算正确，如何能找到一个米散射计算结果准确又有效的数据库，来验证自己算法及程序的正确性。
倒退式算法的总结：
由于Mi.级数崗也般速应隔特狂侑x的增大而减慢讥仮此不同的*佰即使在楝同的计專轴屋娈求下疔需级数顼號叫也不一悴•①川坪决冈学者W.J.在兀變计莽的基砒上，塞右
前人的工作总结出一个\的经脸计算公式.利用该公式给出的项数可以使帮个计算误差小于107.Wiscomlrf公式尼
丄
Af+4^3+】
0*02w工壬盘
叫=4
j+4.05x3+1
81
■X丰4x5十2
4200三工壬20000
当J小于042时*可以利用lluyleigh公式计算，而当X比20000更大时•几何光学即可适用•
Dn的计算采用Dave的倒推式:
[n/切口+Dn（血僅丿]
由于Dn函数有很强的收敛性对于Dn的倒推计算的初值的选取有很强的随意性。
因为当n—x时Dn（ma）—0所以可以取0作为初值。
倒推起点选取大一些，可以保证Dn函数的收敛完全，但是同时却增加了计算时间。
所以必须选取一个最佳的选择标准。
通过试算，作者认为最佳的上限为
Arstop=1*5ft；i^+10
这里ml是复折射率的实部.
同样，对于贝塞耳函数Jn的计算也可以用倒推的方法计算产生：
*丿?
JF■1*,*..
Jn_2心丿=J「1（^）"Jn（久丿
flL
上式是一个普通的Jn的递推式,知道了Jn和Jn-1，可以顺利地计算出所有的Jn序列值。
为了避免计算Jn的繁琐而又能发挥递推式的快速的优点，采用下面的
**
办法:
假设N—x时，取某一个递推初始值为：
JN（）0,JN1（）,
其中&是一个很小的数，如可取10-6。
将初值代入上式，就可以算出所有的J*。
观察同一自变量的J*和J序列,发现它们对应项之间有固定的倍数关系。
如定义这
个倍数为B，那么
■
J畀血=0旳川（a）
由于J1（a）的计算是非常便利的（J1=Sina/a2-COSa/a），所以B=J1/J1*,计算出Jn*（a）可以算出Jn（a）。
和Dn的计算一样,Jn的倒推起始点的公式为：
N^top=h+10
关于贝塞尔函数的倒退过程在另一文献中的描述
Bessel函数的计算厉来被视为数值计算的难题，用其两项递推公式从低阶往往得不出所碍正确的高阶值.经脸袤明「当阶散血大于*时•从零阶与一阶函数向上递推出得人完全不可靠•此时找们说递推时产生了所谓的循环谋劉川*因此对虚豪■的球Bom!
砸数只能用逆推法或级数展开怯来计算计算溟差随起始高阶函数的阶数的増大而减少•文献［7］给岀了计算起始阶数的经驗公式
I・O3jc+65
10<45
45sr<2M
250s»
此时起始高阶函数近似值选为
I=oA-io-*2
用两项递推公式几」=纽亍七-為和得出人之后，再由下列2试得汗£臥的I函数
X-亠
:
sinx：
•-九
smx护U
前口耳二0
■■
Jo人G、=’
**
Jo-1
苴中nIf-1,Af-2,M-3,.J.0
利用初始值
Jnsin
Jo
JnCOS
J1
sin
0
sin
0
升级会员 