Mie 散乱

Gustov Mie は 1908 年に光が波長と同程度の球状物体に当たったときの散乱についての解析的な研究成果を発表しています。

これは散乱に関しては最も単純な問題といえます。

>	with(plots):with(plottools):

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

>	l := display(sphere([0,0,0],2), style=wireframe, color=red):

>	l:=l,seq(arrow(vector([0, -4, k/2+0.5]), vector([0, 10, 0]), .1, .2, 0.2, cylindrical_arrow,color=blue),k=1..5):

>	l:=l,seq(arrow(vector([0, -4, -k/2-0.5]), vector([0, 10, 0]), .1, .2, 0.2, cylindrical_arrow,color=blue),k=1..5):

>	l:=l,seq(arrow(vector([0, -4, k/4-0.8]), vector([0, 2.5, 0]), .1, .2, 0.2, cylindrical_arrow,color=blue),k=1..5):

>	l:=l,seq(arrow(vector([0, -1.5, k/4-0.8]), vector([0,7, 3*k-10.5]), .1, .2, 0.2, cylindrical_arrow,color=blue),k=1..6):

>	display(l,orientation=[-20,80],scaling=constrained);

Mie 理論

Mie 理論は光の粒子による散乱について説明しています。ここで”粒子”は周りの媒質 (medium) とは違う屈折率の球状の材質からなるものです。

そのような粒子の分子内の電子の振動からの双極子放射のパターンは強い散乱放射の根源となります。またすべての双極子からの放射パターンはすべて打ち消すことはなく、均一媒質に対してのときのように入射光の順方向以外ですべて打ち消しあうわけではありません。 むしろ放射ハ゜ターンに関して建設的にも破壊的にも干渉します。こうして粒子による光の散乱は方向によって効率をかえることになります。

Gustav Mie は 1908に任意のサイズの均一材質の球状粒子の散乱問題の解を発表しました。  Mieの論文では解は二つのパラメータ n r と x で説明しています。
粒子と媒質間の屈折率の相違の大きさは粒子の屈折率 n pと媒質の屈折率 n medの比であらわされています。

n r  = n p /n med

屈折率の違いの境界面は電磁波エネルギの再放射のアンテナとなりますがサイズパラメータ(x)であらわします。これは境界面の球の直径 (2a, ここで半径は a)と波長 (λ/ n med)の比です。

 Mie理論の計算は、粒子の真の幾何的断面面積A=π [ ]に対する散乱断面積(cross-sectional area of scattering） 　[ ]に比例する散乱の効率の計算を行います。

特定の幾何サイズの球として理想化された粒子の散乱を考えます。この球は入射した光子を新しい方向に曲げてしまうと考えてください。そして光子の軸に沿った透過を妨げます。従って影ができます。この説明では大雑把すぎます。しかしこれは散乱係数の本質を掴む簡単な概念を与えてくれます。これは吸収係数の考えと良く似ています。散乱の影のサイズは有効断面積（effective cross-section） ( σ s [ ])と呼ばれます。 これは散乱粒子の幾何サイズ (A [ ])よりは小さくなります。散乱効率(scattering efficiency） Qs [単位は無次元]と呼ばれる比例定数に比例しています。

>	sigma[s]=Q[s]*A[s];

[ ]   [-][ ]

散乱係数 (scattering coefficient)    [ ] は体積密度(volume density) s [ ]として記述された集中した個所に多くの散乱粒子を含む媒質を記述します。散乱係数は本質的に媒質の単位体積あたりの断面積です。

>	mu[s]=rho[s]*sigma[s];

[ ]   [ ][ ]

>

restart:

>	with(plots):with(plottools):

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

>	l := display(sphere([0,0,0],2),style=wireframe,color=red):

>	l:=l,seq(arrow(vector([0, -4, k/2+0.5]), vector([0, 10, 0]), .1, .2, 0.2, cylindrical_arrow,color=yellow),k=1..5):

>	l:=l,seq(arrow(vector([0, -4, -k/2-0.5]), vector([0, 10, 0]), .1, .2, 0.2, cylindrical_arrow,color=yellow),k=1..5):

>	l:=l,seq(arrow(vector([0, -4, k/4-0.8]), vector([0, 2.5, 0]), .1, .2, 0.2, cylindrical_arrow,color=yellow),k=1..5):

>	l:=l,seq(arrow(vector([0, -1.5, k/4-0.8]), vector([0,7, 3*k-10.5]), .1, .2, 0.2, cylindrical_arrow,color=yellow),k=1..6):

>	l:=l,arrow(vector([4, 0,-4.5]), vector([-3,0,3]), .1, .2, 0.2, cylindrical_arrow,color=black):

>	l:=l,textplot3d([4,0,-5,`A 幾何的散乱断面積`],color=black):

>	l:=l,arrow(vector([0,5,0]), vector([0,-1,0]), .1, .2, 0.2, cylindrical_arrow,color=black):

>	l:=l,textplot3d([0,8,0,`σｓ=Qｓ・A 有効断面積`],color=black):

>	l:=l,plot3d( [r*cos(phi),0,r*sin(phi)], r=0..2, phi=0..2*Pi,style=patchnogrid,color=green):

>	l:=l,plot3d( [r*cos(phi),3,r*sin(phi)], r=0..1, phi=0..2*Pi,style=patchnogrid,color=blue):

>	display(l,orientation=[-20,80],scaling=constrained);

実際上   に対する単位 [ ] は積 L が無次元となる長さの逆数となります。ここでL [cm]は光子の媒質内での伝搬の光路長です。 光路L後の散乱によって方向を変えられない光子の透過確率T はつぎのようになります。

>	T=exp(-mu[s]*L);

Mie散乱の数学

光源、球状散乱粒子そして観測点を考えます。この三点で散乱面と呼ばれる面が定義されます。入射光および散乱光は散乱面に平行および直交する要素に既約することができます。以下の図に示されるように、平行および直交成分は、散乱面に平行および直交する線形偏光フィルタによって実際的に選択することができます。

>

restart:

>	with(plots):with(plottools):

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

>	l := display(sphere([0,0,0],0.1),style=patchnogrid,color=red):

>	f := cuboid([-1,-1,-1],[1,1,1]): l:=l,display(f,style=wireframe,thickness=3,color=black):

>	l :=l, polygon([[0,1,1],[0,-1,1],[0,-1,-1],[0,1,-1]], color=green): l :=l,line([0,-3,0], [0,1.1,0],color=black,thickness=3):

>	l:=l,arrow(vector([0, -3, 0]), vector([0, 1, 0]), .1, .2, 0.2, cylindrical_arrow,color=red):

>	l:=l,textplot3d([0,-3,0.2,`光源より`],color=black):

>	l :=l,line([0,0,0], [0,-2,1.5],color=black,thickness=3):

>	l:=l,arrow(vector([0, -2, 1.5]), vector([0,-2/3,1.5/3]), .1, .2, 0.2, cylindrical_arrow,color=red):

>	l:=l,textplot3d([0,-3.2,2.2,`検出器へ`],color=black):

>	l:=l,textplot3d([0.2,0,0.2,`粒子`],color=black):

>	l:=l,textplot3d([0.2,0,1.5,`散乱角θ`],color=black):

>	l:=l,textplot3d([0.6,0.6,-0.2,`散乱面`],color=black):

>	cc:= polygon([[0,0.2,0.2],[0,0.2,-0.2],[0,-0.2,-0.2],[0,-0.2,0.2]], color=yellow):

>	cc:=rotate(cc,0,0,Pi/2):cc:=translate(cc,0,-1.9,0):

>

l:=l,cc:

>	cc:= polygon([[0,0.2,0.2],[0,0.2,-0.2],[0,-0.2,-0.2],[0,-0.2,0.2]], color=yellow):

>	cc:=rotate(cc,0,Pi/4,Pi/2):cc:=translate(cc,0,-1.8,1.4):

>

l:=l,cc:

>	l:=l,textplot3d([0,-1.5,-0.5,`線形偏光器`],color=black):

>	l:=l,textplot3d([-1.3,0,1.3,`線形偏光器`],color=black):

>	display(l,orientation=[-50,80],scaling=constrained);

散乱行列は、入射光とファーフィールド (far-field) として観測される散乱面に平行および直交成分の電界フィールドの、散乱面に平行および直交成分の入射光と散乱光の関係を記述します。 (Bohren and Huffman を参照):

>	restart:with(plots):with(plottools):

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

>	print(<<E[s,par]>,<E[s,per]>>,"=",exp(-ik*(r-z))/(-ik*r),<<S[2]|S[3]>,<S[4]|S[1]>>,<<E[i,par]>,<E[i,per]>>);

上の式は実用的に単純化されています。
指数関数項  は、散乱と観測点間の距離に依存する伝達 (transport) 因子 です。 散乱光を、

例えば偏光の角度あるいは方向の関数として散乱体から一定距離 r の地点で測定すると、 伝達因子は定数になります。
トータルフィールド (Etot) は入射フィールド (Ei)、散乱フィールド (Es) そしてこれらのフィールドの相互作用 (Eint) に依存します。 もし Ei を避ける位置から散乱を観測すると、 Ei と Eint の両方がゼロであり Es のみが観測されます。

kr >> 、k = , n c =  であるような直径 d の粒子から r の距離で Es のファーフィールドの観測に対して散乱要素 S3 および S4 はゼロに等しくなります。 (Bohren and Huffman の Eq. 4.75 を参照)。

実際の実験の測定強度 I = <E E*> =   です。ここで E = a exp(-iδ), a は電界フィールドの振幅で  は位相です。こうして実際的な散乱の測定では上の式は次のように単純化されます。

>	print(<<I[i,par]>,<I[i,per]>>,"=",Constant*<<abs(S[2])^2|0>,<0|abs(S[1])^2>>,<<I[s,par]>,<I[s,per]>>):

Maxwell 方程式により境界条件の線形性は、任意の粒子から散乱される平面波の入射光と散乱光の電界フィールド間の関係は行列形式で簡潔に表現されます。

 ここで行列は "振幅散乱行列" (van de Hulst 1957, Bohren and Huffman 1983,および Goody and Yung 1989) として知られています。この行列の要素のラベルは Hansen および Travis (1974) と同じであることに注意してください。
Hansen および Travis (1974) がそれぞれ垂直 (perpendicula) および平行 (paralle) の略の r および l を使ったことを除いて同じです。複素指数関数にマイナスの符号を持っているものがあります(これは前にあるマイナスの記号を消去します)。行列方程式は、原点から十分離れていれば任意の距離に対して有効です。

散乱粒子が球状であるとき特に単純です。この場合、球のアイソトロピ (等方性) は、散乱電界フィールドがその初期フィールドに直行する成分を一切持たないことを課します。散乱係数は複素数で角度依存の量です。従って振幅と位相の両方は、完全に特性化するためには振幅と位相の両方ともがすべての角度に対し決定されなければなりません。従って 球状粒子の場合そのような測定は可能です。ここで対称性は半分だけ測定すれば良いことを意味します。しかし非球状粒子の場合経験的に難しいです。しかしながら どの角度の一つのゼロでない測定でさえ、散乱粒子が非球状であることを示すのはあきらかです。

実験的に散乱係数を決定する問題は電界フィールドの変換ではなく簡単に測定できる Stokes パラメータを考えることによってより素直に形成されます。 この式は散乱行列を使って書くことができます。
 

偏光光の角度散乱パターン

Mie理論は散乱行列の二つの要素 S1() と S2() の角度依存性を計算します。これにより偏光光の散乱強度が計算されます。散乱パターンはまた粒子による散乱のアンアイソトロピー g を計算するのに使用されます。

Mie散乱の角度パターン

水の中の直径 0.304 ミクロンの非吸収ポリスチレンの球の HeNe レーザビームによる散乱パターンを考えます。

 = /  = 1.5721/1.3316 = 1.1806  相対屈折率
"x = 2a/(λ/ ) = (2)(3.1415)(0.152)/(0.6328/1.3316) = 2.0097"  サイズ・パラメータ

Mie 理論のアルゴリズムを実行:
Mie(nr, x) ---> S1(), S1() 複素数として計算します。
Mie(1.1806, 2.0097) ---> S1.re + jS1.im, S2.re + jS2.im

結果を見るためには変更された光源/検出器のペアの平行と直行の方向について散乱の強度を計算します。

Ipar = S2・S2* = Re{(S2.re + jS2.im)(S2.re - jS2.im)}

Iper = S1・S1* = Re{(S1.re + jS1.im)(S1.re - jS1.im)}

以下の図に示され説明されているように実際的に測定できます。

以下の図は散乱強度Ipar および Iperの角度依存性を説明する実際的な測定を記述しています。
Iper
水の中の球のdilute（弱い）解のテーブルに直交するレーザビーム偏光光を放射します。テーブルに平行な面内で散乱光を角度の関数として測定します。検出器の前にテーブルに垂直な線形偏光フィルターを置きます。
Ipar
水の中の球のdilute解のテーブルに平行なレーザービームの偏光光を放射します。テーブルに平行な面内で散乱光を角度の関数として測定します。 検出器の前にテーブルに平行な線形偏光フィルターを置きます。

>	l := display(sphere([0,0,0],1),style=wireframe,color=red):

>	l :=l,line([-5, 5,0], [5, 5, 0],color=black,thickness=1): l :=l,line([-5,-5,0], [5,-5, 0],color=black,thickness=1): l :=l,line([-5, 5,0], [-5,-5,0],color=black,thickness=1): l :=l,line([5,-5, 0], [5, 5, 0],color=black,thickness=1):

>	l :=l, polygon([[-4,1,1],[-4,-1,1],[-4,-1,-1],[-4,1,-1]], color=green): l :=l,line([-6,0,0], [4,0,0],color=black,thickness=2):

>	l:=l,arrow(vector([-6,0, 0]), vector([2, 0, 0]), .2, .4, 0.2, cylindrical_arrow,color=red):

>	l :=l,line([0,0,0], [3,-3,0],color=black,thickness=3): l:=l,arrow(vector([3, -3,0]), vector([1,-1,0]), .2, .4, 0.2, cylindrical_arrow,color=red):

>	cc:= polygon([[0,1,1],[0,1,-1],[0,-1,-1],[0,-1,1]], color=green):

>	cc:=rotate(cc,0,0,Pi/4):cc:=translate(cc,3,-3,0):

>

l:=l,cc:

>	l:=l,textplot3d([0,0,1.1,`粒子`],color=black):

>	l:=l,textplot3d([-6,0,0.5,`入射レーザ光`],color=black):

>	l:=l,textplot3d([-3,0,0.5,`線形偏光子`],color=black):

>	l:=l,textplot3d([5,-2,0,`線形偏光子`],color=black):

>	l:=l,textplot3d([5,-5,0,`検出器`],color=black):

>	l:=l,textplot3d([1.5,-1,0,`θ`],color=black):

>	display(l,orientation=[-80,70],scaling=constrained);

例：Ipar および Iperに対する散乱ハ゜ターンの極座標およびｘｙ座標プロット。

ランダムに変更された光源に対し全散乱光の強度は の項で与えられます。

>	S[11](theta)=1/2*abs(S[1](theta))^2+1/2*abs(S[2](theta))^2;

はMueller行列と呼ばれる行列の最初の要素です。行列は 4*4 の行列で Stokes パラメータ (Ii, Qi, Ui, Vi) の入力ベクトルに関連しており複素光源を記述します。 出力ベクトル (Is, Qs, Us, Vs) は伝達した光の性質を記述します。ランダムに変更された光に対し は伝達 (transport) のトータル強度を記述します。

>	I[s] = S[1,1]*I[i];

球状粒子の散乱

球状粒子による散乱は二つの散乱係数が対称性からゼロになるために最も単純な散乱です。数式的な厳密解は Gustav Mie (1908) によるものです。

Mie散乱の解は微視的な Maxwell 方程式に始まります。電界フィールド E と磁界フィールド H の複素数表現を使って表されます。

>

restart:

>	grad(E)=0;grad(H)=0;curl(E)=i*omega*mu*H;curl(H)=-i*omega*epsilon*E;

ここで

>	k^2=omega^2*epsilon*mu;

ベクトル波動方程式は以下のようになります。

>	Delta*E+k^2*E=0;

>	Delta*H+k^2*H=0;

つぎの段階はスカラー関数Ψを定義することです。二つのベクトル M と N から適切なベクトル微分を適用することで構成されます。関数と二つのベクトルがこのように定義されることから、もし  がスカラー波動方程式の解であるなら、M と Nは同じベクトル波動方程式の解です。さらに適切に定義されるならM と N はまた電磁フィールドの必要性を満足します。この関数は数理物理においてよく遭遇するスカラーの球状調和関数との類似性からベクトル球状調和関数と呼ばれます。

これはスカラーの波動方程式の解を得るための直接的な実験です。粒子は球状対称と仮定しているので、球面座標は表現すべき座標系としては明らかな選択です。ところで、他の形状の散乱問題も、この考えで、特定の粒子が特に単純な式を仮定できるように座標系を選択することにより、他の粒子形状の散乱問題が多分解けるであろうということは明らかです。

球状の問題に戻って球状波動方程式の直接的変数分離はすぐに一般解を導きます。

>	psi=sum(Sum(A[l]^m*cos(m*phi)*P[l]^m*(cos(theta))*z[n](kr)+B[l]^m*sin(m*phi)*P[l]^m*(cos(theta)),m=-l..l),l=0..infinity);

ここで  は Legendre 多項式で  は球状 Bessel 関数です。この無限級数はまだ問題の解にはなっていません。難しい点は、境界条件が球座標で指定されて (粒子の対称性を保持しつつ） しかし入射波は平面波で 直交 (Cartesian) 座標で与えられている点です。入射波と境界条件を同じ座標系で表現するためには平面波をベクトル球状調和関数で展開するか球面境界条件を直行座標で表現するかです。因習的ですが意味のない明示的なアプローチは平面波をベクトル球状調和関数で展開することです。計算はかなり複雑になります。従ってここでは解のみ示します。複素平面波の球状調和関数による目的の展開は次のように与えられます。

>	E[c]=sum(sum(B[m,n]^e*M[m,n]^e+B[m,n]^o*M[m,n]^o+A[m,n]^e*N[m,n]^e+A[m,n]^o*N[m,n]^o,'n'=0..infinity),m=-l..l);

もし平面波が x 軸方向に偏光されていて z 軸にそって進んでいる場合、m = 1 のみ寄与します。従って展開は次のようになります。

>	E[c]=sum(i^n*(2*n+1)/n/(n+1)*(M[1,n]^o(1)-i*N[1,n]^e(1)),'n'=1..infinity);

この展開に発生するベクトル球状調和関数の組み合わせは二つのスカラー成分に分かれます。それぞれは第一種の球状ベッセル関数を含みます。球の外部の初期フィールドはいまや既知です。

さらに二つのフィールドを計算する必要があります。球の外側の散乱波と球自体の内部に発生する波です。両方とも入射波の展開の式に似ています。しかし球の外側の散乱波は第一種の球状 Hankel 関数を含みます。一方内部では第一種の球状 Bessel 関数を含みます。この明らかな非対称性の理由はすでに説明されています。Hankel 関数は第一種 Bessel 関数の線形結合です。この関数は無限大で特異です。第二種の Bessel 関数 (Weber あるいは Neumann 関数として知られています) は原点で特異です。内部の波が第一種の Bessel 関数で構成されなければならないのは第二種の Bessel 関数の非物理学的な原点での特異性によるためです。

境界条件は、電磁フィールドの横成分のトータルは (入射光と散乱光のフィールドの和から内部フィールドを引いたものに等しい) 境界面においてゼロに等しいということが必要となります。上で決定されたフィールドは、いま境界条件に組み込まれます。4 つの式が得られます。内部および散乱フィールドの各成分の式です。しかしながら、もし粒子と媒質の透過性が等しければ、4 つの式は 2 つの式のみになります。 得えられた定数の級数は  と  と記されます。これらの式は、少し簡単な式で Goody と Yung により与えられました (1989)。補遺の関数を定義することにより、E と H フィールドは合理的に簡潔な式で表現されます。これらの式から振幅関数  と  が最終的に得られます。

>	S[1](theta)=sum((2*n+1)/n/(n+1)*(a[n]*Pi[n](cos(theta))+b[n]*tau[n](cos(theta))),'n'=1..infinity);

>	S[2](theta)=sum((2*n+1)/n/(n+1)*(b[n]*Pi[n](cos(theta))+a[n]*tau[n](cos(theta))),'n'=1..infinity);

吸収と散乱係数の式はもう少し単純です。 散乱振幅関数に使われる E と H に対する級数は適切な積分に組み込まれます。そして積分は項ごとに引き継がれます。この計算の結果は以下のようになります。

>	Q[ext]=2/x^2*sum((2*n+1)*Re(a[n]+b[n]),'n'=0..infinity);

>	Q[sca]=2/x^2*sum((2*n+1)*Re(abs(a[n])^2+abs(b[n])^2),'n'=0..infinity);

ここの方程式は Goody と Yung (1989) によって修正されたものです (erroneous equation (7.67a)）。これらの解を得るための努力は多大なるものでした。しかしながら散乱問題の解析的な解は既知のものとなったのです。不幸なことに計算すべき無限級数は実用的でないくらい収束が非常に遅いものです。実際、サイズ因子 (size factor) の球に対する正確な結果を得るためには近似的 x 項が計算されなければなりません。

>	x=2*Pi*a/lambda;

これらの級数の実際的な計算はコンピュータの助けを必要とする必要があります。r>>λにたいしては、現在改良されたMie理論を使った計算でさえもさらに現実的でなくなります。

参考文献：
Born, M. and Wolf, E. "Diffraction by a Conducting Sphere; Theory of Mie." §13.5 in Principles of Optics: Electromagnetic Theory of Propagation, Interference, and Diffraction of Light, 7th ed. Cambridge, England: Cambridge University Press, p. 633-644, 1999.
Goody, R. M. and Yung, Y. L. Atmospheric Radiation: Theoretical Basis, 2nd ed. New York: Oxford University Press, pp. 315-316, 1989.
Kerker, M. The Scattering of Light and Other Electromagnetic Radiation. New York: Academic Press, pp. 54-59, 1969.
van de Hulst, H. C. Light Scattering by Small Particles. New York: Dover, 1981.
Author: Eric W. Weisstein

対数微分（Logarithmic Derivative）

Mie 散乱の振幅 (amplitude) 関数 S1, S2 はつぎの式で計算されます。

>	S[1](theta)=sum((2*n+1)/n/(n+1)*(a[n]*Pi[n](cos(theta))+b[n]*tau[n](cos(theta))),'n'=1..infinity);

>	S[2](theta)=sum((2*n+1)/n/(n+1)*(b[n]*Pi[n](cos(theta))+a[n]*tau[n](cos(theta))),'n'=1..infinity);

散乱係数 a[n] と b[n] (一般的には複素数) はサイズパラメータ x と複素数の相対屈折率 m に依存しそして様々なマルチモードの寄与を特性化します（Dave １９６９を参照)

>	a[n]=(((D[n](z)/m)+n/x)*psi[n](x)-psi[n-1](x))/(((D[n](z)/m)+n/x)*xi[n](x)*xi[n](x)-xi[n-1](x));

>	b[n]=(((m*D[n](z))+n/x)*psi[n](x)-psi[n-1](x))/((m*(D[n](z))+n/x)*xi[n](x)*xi[n](x)-xi[n-1](x));

ここで量  と  は Ricatti - Bessel 関数で球状ベッセル関数のそれぞれ第一種の J[n] と 2 種の Y[n] です。（Bronstein＆Semendjajew１９９１参照）

>	psi[n](x)=x*J[n](x);

>	xi[n](x)=psi[n](x)-I*x*Y[n](x);

関数 Dn はつぎの様に定義される対数微分 (Logarithmic Derivative) です。

>	D[n](z)=Diff(ln(psi[n](z)),z);

>	D[n](z)=-n/z+psi[n-1](z)/psi[n](z);

Voshchinnikov（2002) によると実引数のベッセル関数の計算は第二種の Yn+1/2 の関数の下方向繰り返し計算に基づきます。

第一種ベッセル関数は次の関係式で決定されます。

>	J[n+1/2](x)=(2/Pi/x+Y[n+1/2](x)*J[n-1/2](x))/Y[n-1/2](x);

複素引数のベッセル関数の対数微分（Logarithmic Derivative）は逆方向（バックワード）で計算されます。

関数  及び  は角度依存の情報を持っています。

>	Pi[n](cos(theta))=1/sin(theta)*P1[n](cos(theta));

>	P1[n](z)=simplify(LegendreP(1,n,z),`LegendreP`);

>	tau[n](cos(theta))=Diff(P1[n](cos(theta)),theta);

ここで P1[n] は第一種ルジャンドル関数です。

その他の効率因子

消光 (extinction) 効率因子は次のように与えられます。

>	Q[ext]=2/x^2*sum((2*n+1)*Re(a[n]+b[n]),'n'=0..infinity);

散乱 (scattering) 効率因子はつぎのように与えられます。

>	Q[sca]=2/x^2*sum((2*n+1)*Re(abs(a[n])^2+abs(b[n])^2),'n'=0..infinity);

後方散乱有功因子は次のように与えられます。

>	Q[back]=1/x*abs(sum((2*n+1)*(-1)^n*(a[n]-b[n]),'n'=0..infinity))^2;

吸収効率は次のように導かれます。

>	Q[abs]=Q[ext]-Q[sca];

放射圧効率 (Radiation pressure efficiency) は次のように与えられます。

>	Q[pr]=Q[ext]-g*Q[sca];

アルベド (Albedo) はつぎのようになります。

>	Lambda=Q[sca]/Q[ext];

散乱非対称性 (scattering asymmetry) パラメータ g

>	g=<cos(theta)>;

>	g=Int(cos(theta)*F[11]/(4*Pi),Omega);

散乱断面積 (cross section）

>

C=G*Q;

ここで G は粒子の幾何学的断面積です。

>

G=Pi*a^2;

Mapleによる計算

ここでは Maple を使って Mie 散乱を実際に計算してみましょう。

>

restart;

>	with(plots):with(plottools):

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

球上粒子の直径を定義します。単位はミクロン。

Sphere diameter  microns

>	diameter:=1;

>	radius:=diameter/2;

真空中での波長。単位はミクロン。

Wavelength in Vacuum  microns

>	lambda:=0.6328;

媒質の屈折率。（Index of Refraction in Medium）

Real（実部） 屈折率（Index of Sphere）  
Imag（虚部） 消光係数

>	m:=1.5-0.000*I;

密度 (Density)

立方ミクロンあたりの粒子数 (spheres per cubic micron)

>	density:=2;

角度の数 (分割数：Number of angles)

>	nangles:=50;

サイズパラメータの計算

>	x:=2*Pi*radius/lambda;evalf(x);

 のベクトル宣言と初期設定

>	mu:=Vector(nangles):

>	for i to nangles do mu[i]:=evalf(cos(Pi/nangles*(i-1))) end do:

級数展開の停止項

>	nstop:=floor(x+4.05*x^(1/3)+2);

s1, s2, , ,  のベクトル宣言と初期設定

>	s1:=Vector(nangles,0):

>	s2:=Vector(nangles,0):

>	pi0:=Vector(nangles,0):

>	pi1:=Vector(nangles,1):

>	tau:=Vector(nangles,0):

ｚ の計算

>	z:=x*m:evalf(z);

Lenz の対数微分 (logarithmic derivative) 計算の関数定義

>	Lenz_Dn:=proc(z,n)local zinv,alpha,aj,alpha_j1,alpha_j2,ratio,runratio,result;

>	zinv:=2/z;

>	alpha:=(n+1/2)*zinv;

>	aj:=-(n+1+1/2)*zinv;

>	alpha_j1:=aj+1/alpha;

>	alpha_j2:=aj;

>	ratio:=alpha_j1/alpha_j2;

>	runratio:=alpha*ratio;

>	while abs(abs(ratio)-1) >1e-12 do

>	aj:=zinv-aj;

>	alpha_j1:=1/alpha_j1+aj;

>	alpha_j2:=1/alpha_j2+aj;

>	ratio:=alpha_j1/alpah_j2;

>	zinv:=-zinv;

>	runratio:=ratio*runratio;

>

end do;

>	result:=-n/z+runratio;

>	return result;

>

end:

対数微分 (Logarthsmic Derivative)  のベクトル宣言と初期設定

>	DD:=Vector(nstop+1,0):

>	print(abs(Im(m)*x)),((13.78*Re(m)-10.8)*Re(m)+3.9);

 の計算

>	if abs(Im(m)*x)<((13.78*Re(m)-10.8)*Re(m)+3.9) then

>	DD[1]:=evalf(1/tan(z));

>	zinv:=1/z;

>	for k from 2 to nstop+1 do

>	k_over_z:=(k-1)*zinv;

>	DD[k]:=evalf(1/(k_over_z-DD[k-1])-k_over_z);

>

end do:

>

else

>	DD[nstop+1]:=Lentz_Dn(z,nstop);

>	zinv:=1/z;

>	for k from nstop by -1 to 1 do

>	k_over_z:=k*zinv;

>	DD[k]:=evalf(k_over_z-1/(DD[k+1]+k_over_z));

>

end do:

>

end if:

 の実部と虚部をプロットしてみましょう。

>	RD:=array(1..nstop):ID:=array(1..nstop):

>	for k to nstop do

>	RD[k]:=evalf(Re(DD[k])):ID[k]:=evalf(Im(DD[k])):

>

end do:

>	PP:=array(1..2):PP[1]:=listplot(RD,color=red):PP[2]:=listplot(ID,color=blue):display(PP);

 と  の初期設定

>	psi0:=sin(x);

>	psi1:=psi0/x-cos(x);

 と  の初期設定

>	xi0:=psi0-cos(x)*I;

>	xi1:=psi1-(cos(x)/x+sin(x))*I;

各種初期設定

>	Q[sca]:=0;

>

g:=0;

>	Q[ext]:=0;

>

sign1:=1;

>	Q[bcalc]:=0;

>	a[nn-1]:=0;

>	b[nn-1]:=0;

繰り返し計算（アルゴリズムの中心部）

>	for n to nstop do if Re(m) = 0 then

>	a[nn]:=(n*psi1/x-psi0)/(n/x*xi1-xi0);

>	b[nn]:=psi1/xi1;

>	elif Im(m) = 0 then

>	z1:=Re(DD[n+1])/Re(m)+n/x;

>	a[nn]:=evalf((Re(z1)*psi1-psi0)/(Re(z1)*xi1-xi0));

>	z1:=Re(DD[n+1])*Re(m)+n/x;

>	b[nn]:=evalf((Re(z1)*psi1-psi0)/(Re(z1)*xi1-xi0));

>

else

>	z1:=DD[n+1]/m+n/x;

>	a[nn]:=evalf((Re(z1)*psi1-psi0+(Im(z1)+psi1)*I)/(z1*xi1-xi0));

>	z1:=DD[n+1]*m+n/x;

>	b[nn]:=evalf((Re(z1)*psi1-psi0+(Im(z1)+psi1)*I)/(z1*xi1-xi0));

>	end if; # 級数計算

>	for k from 1 to nangles do

>	factor0:=(2*n+1)/(n+1)/n;

>	tau[k]:=evalf(n*mu[k]*pi1[k]-(n+1)*pi0[k]);

>	alpha:=factor0*pi1[k];

>	beta:=factor0*tau[k]; # S1、S2の計算

>	s1[k]:=evalf(s1[k]+alpha*a[nn]+beta*b[nn]);

>	s2[k]:=evalf(s2[k]+alpha*b[nn]+beta*a[nn]);

>

od:

>

>	for k from 1 to nangles do

>	factor0:=pi1[k];

>	pi1[k]:=evalf(((2*n+1)*mu[k]*pi1[k]-(n+1)*pi0[k])/n);

>	pi0[k]:=factor0;

>

od:

>	factor0:=2*n+1;

>	g:=evalf(g+(n-1/n)*(Re(a[nn-1])*Re(a[nn])+Im(a[nn-1])*Im(a[nn])+Re(b[nn-1])*Re(b[nn])+Im(b[nn-1])*Im(b[nn])));

>	g:=evalf(g+factor0/n/(n+1)*(Re(a[nn])*Re(b[nn])+Im(a[nn])*Im(b[nn])));

>	Q[sca]:=evalf(Q[sca]+factor0*(abs(a[nn])^2+abs(b[nn])^2));

>	Q[ext]:=evalf(Q[ext]+factor0*(Re(a[nn])+Re(b[nn])));

>

>	sign1:=-sign1;

>	Q[bcalc]:=evalf(Q[bcalc]+sign1*factor0*(a[nn]-b[nn]));

>	factor0:=(2*n+1)/x;

>	xi:=factor0*xi1-xi0;

>

xi0:=xi1;

>

xi1:=xi;

>	psi:=factor0*psi1-psi0;

>	psi0:=psi1;

>	psi1:=Re(xi1);

>	a[nn-1]:=a[nn];

>	b[nn-1]:=b[nn];

>	print(n,DD[n+1],Q[ext],a[nn],b[nn],Q[sca]);

>

end do:

s1 の実部と虚部をプロット

>	Rs1:=array(1..nangles):Is1:=array(1..nangles):

>	for k to nangles do

>	Rs1[k]:=evalf(Re(s1[k])):Is1[k]:=evalf(Im(s1[k])):

>

end do:

>	p1:=listplot(Rs1,color=red):p2:=listplot(Is1,color=blue):display([p1,p2]);

 のプロット

>	listplot(tau);

 と  のプロット

>	p1:=listplot(pi0,color=red):p2:=listplot(pi1,color=blue):display([p1,p2]);

以下のページでJavaアプレットで計算できますので比較してみましょう。

Mie計算のプログラム

これをプロシージャのプログラムにします。

>

restart;

>	with(plots):with(plottools):

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

各パラメータの値を設定しておきます。

>	diameter:=1;

>	radius:=diameter/2;

>	lambda:=0.6328;

>	m:=1.5-0.000*I;

>	density:=2;

>	nangles:=50;

以下、プロシージャ関数 Mie の定義です。

>	Mie:=proc(radius,lambda,m,density,nangles) local x,mu,i,nstop,s1,s2,pi0,pi1,tau,z,Lenz_Dn,DD,zinv,k,k_over_z,RD,ID,psi0,psi1,xi0,xi1,Q,g,sign1,a,b,n,z1,factor0,alpha,beta,xi,psi,p1,p2,parallel,perpen,phasefn:

>	x:=2*Pi*radius/lambda;x:=evalf(x);

>	mu:=Vector(nangles):

>	for i to nangles do mu[i]:=evalf(cos(Pi/nangles*(i-1))) end do:

>	nstop:=floor(x+4.05*x^(1/3)+2):

>	s1:=Vector(nangles,0):

>	s2:=Vector(nangles,0):

>	pi0:=Vector(nangles,0):

>	pi1:=Vector(nangles,1):

>	tau:=Vector(nangles,0):

>	z:=x*m:evalf(z):

>	# print('z',z);

>	Lenz_Dn:=proc(z,n)local zinv,alpha,aj,alpha_j1,alpha_j2,ratio,runratio,result;

>	zinv:=2/z;

>	alpha:=(n+1/2)*zinv;

>	aj:=-(n+1+1/2)*zinv;

>	alpha_j1:=aj+1/alpha;

>	alpha_j2:=aj;

>	ratio:=alpha_j1/alpha_j2;

>	runratio:=alpha*ratio;

>	while abs(abs(ratio)-1) >1e-12 do

>	aj:=zinv-aj;

>	alpha_j1:=1/alpha_j1+aj;

>	alpha_j2:=1/alpha_j2+aj;

>	ratio:=alpha_j1/alpah_j2;

>	zinv:=-zinv;

>	runratio:=ratio*runratio;

>

end do;

>	result:=-n/z+runratio;

>	return result;

>	end proc: # print('nstop',nstop);

>	DD:=Vector(nstop+1,0):

>	# print(m,abs(Im(m)*x),(13.78*Re(m)-10.8)*Re(m)+3.9);

>	if abs(Im(m)*x)<((13.78*Re(m)-10.8)*Re(m)+3.9) then

>	# print('z',z); # print('z',z,tan(z)); # if tan(z) = 0 then DD[1]:=0 else DD[1]:=evalf(1/tan(z)) end if; DD[1]:=arctan(z);

>	zinv:=1/z; # print(DD[1],'zinv',zinv);

>	for k from 2 to nstop+1 do

>	k_over_z:=(k-1)*zinv; # print(n,k_over_z-DD[k-1]);

>	DD[k]:=evalf(1/(k_over_z-DD[k-1])-k_over_z);

>

end do:

>

else

>	DD[nstop+1]:=Lentz_Dn(z,nstop); # print('z',z);

>	zinv:=1/z;

>	for k from nstop by -1 to 1 do

>	k_over_z:=k*zinv; # print(DD[k+1]+k_over_z);

>	DD[k]:=evalf(k_over_z-1/(DD[k+1]+k_over_z));

>

end do:

>

end if:

>	# print(seq(DD[k],k=1..5));

>	RD:=array(1..nstop):ID:=array(1..nstop):

>	for k to nstop do

>	RD[k]:=evalf(Re(DD[k])):ID[k]:=evalf(Im(DD[k])):

>

end do:

>

>	seq(RD[k],k=1..5);seq(ID[k],k=1..5);

>

>	psi0:=sin(x); # print('x',x);

>	psi1:=psi0/x-cos(x);

>	xi0:=psi0-cos(x)*I;

>	xi1:=psi1-(cos(x)/x+sin(x))*I;

>	Q[sca]:=0;

>

g:=0;

>	Q[ext]:=0;

>

sign1:=1;

>	Q[bcalc]:=0;

>	a[nn-1]:=0;

>	b[nn-1]:=0;

>	for n to nstop do # print(xi1); if Re(m) = 0 then a[nn]:=(n*psi1/x-psi0)/(n/x*xi1-xi0); b[nn]:=psi1/xi1; elif Im(m) = 0 then z1:=Re(DD[n+1])/Re(m)+n/x; a[nn]:=evalf((Re(z1)*psi1-psi0)/(Re(z1)*xi1-xi0)); z1:=Re(DD[n+1])*Re(m)+n/x; b[nn]:=evalf((Re(z1)*psi1-psi0)/(Re(z1)*xi1-xi0)); else

>

>	z1:=DD[n+1]/m+n/x;

>	a[nn]:=evalf((Re(z1)*psi1-psi0+(Im(z1)+psi1)*I)/(z1*xi1-xi0));

>	z1:=DD[n+1]*m+n/x;

>	b[nn]:=evalf((Re(z1)*psi1-psi0+(Im(z1)+psi1)*I)/(z1*xi1-xi0));

>	end if; for k from 1 to nangles do

>	factor0:=(2*n+1)/(n+1)/n;

>	tau[k]:=evalf(n*mu[k]*pi1[k]-(n+1)*pi0[k]);

>	alpha:=factor0*pi1[k];

>	beta:=factor0*tau[k];

>	s1[k]:=evalf(s1[k]+alpha*a[nn]+beta*b[nn]);

>	s2[k]:=evalf(s2[k]+alpha*b[nn]+beta*a[nn]);

>

od:

>	for k from 1 to nangles do

>	factor0:=pi1[k];

>	pi1[k]:=evalf(((2*n+1)*mu[k]*pi1[k]-(n+1)*pi0[k])/n);

>	pi0[k]:=factor0;

>

od:

>	factor0:=2*n+1;

>	g:=evalf(g+(n-1/n)*(Re(a[nn-1])*Re(a[nn])+Im(a[nn-1])*Im(a[nn]) +

>	Re(b[nn-1])*Re(b[nn])+Im(b[nn-1])*Im(b[nn])));

>	g:=evalf(g+factor0/n/(n+1)*(Re(a[nn])*Re(b[nn])+Im(a[nn])*Im(b[nn])));

>	Q[sca]:=evalf(Q[sca]+factor0*(abs(a[nn])^2+abs(b[nn])^2));

>	Q[ext]:=evalf(Q[ext]+factor0*(Re(a[nn])+Re(b[nn])));

>	sign1:=-sign1;

>	Q[bcalc]:=evalf(Q[bcalc]+sign1*factor0*(a[nn]-b[nn]));

>	factor0:=(2*n+1)/x;

>	xi:=factor0*xi1-xi0;

>

xi0:=xi1;

>

xi1:=xi;

>	psi:=factor0*psi1-psi0;

>	psi0:=psi1;

>	psi1:=Re(xi1);

>	a[nn-1]:=a[nn];

>	b[nn-1]:=b[nn];

>	# print(n,DD[n+1],Q[ext],a[nn],b[nn],Q[sca]);

>

end do:

>

>

>	p1:=listplot(pi0,color=red):p2:=listplot(pi1,color=blue):display([p1,p2]);

>

evalf(x);

>

Q[sca];

>	Q[sca]:=evalf(Q[sca]*2/(x*x));

>	Q[ext]:=evalf(Q[ext]*2/(x*x));

>	g:=evalf(g*4/Q[sca]/(x*x));

>	density*Q[ext];

>	density*Q[sca];

>	Q[bcalc];s1[nangles];

>	Q[back]:=evalf(abs(Q[bcalc])^2/(x*x));

>

>	return s1,s2,Q[sca],Q[ext],Q[back],x:

>

end proc:

関数値 (return) は s1, s2, Q[sca], Q[ext], Q[back], x です。

波長を変えて関数を実際に呼んでみましょう。

>	lambda:=0.6;

>	ans:=Mie(radius,lambda,m,density,nangles);

6 番目がサイズパラメータです。偏光成分を計算します。

>	x:=ans[6];

>	parallel:=Vector(nangles):

>	perpen:=Vector(nangles):

>	phasefn:=Vector(nangles):

>	Q[sca]:=ans[3]:Q[ext]:=ans[4]:Q[back]:=ans[5]:

>	for i from 1 to nangles do

>	parallel[i]:=4*abs(ans[2][i])/(x*x*Q[sca]);

>	perpen[i]:=4*abs(ans[1][i])/(x*x*Q[sca]);

>	phasefn[i]:=(parallel[i]+perpen[i])/2;

>

end do:

偏光成分をプロットしてみましょう。

>	p1:=listplot(parallel,color=red):p2:=listplot(perpen,color=green):p3:=listplot(phasefn,color=blue):display([p1,p2,p3]);

波長依存の計算とアニメーション

波長を変えながらこの関数 Mie を呼び波長応答を計算してみましょう。

0.34 ミクロンから 0.78 ミクロンまで変えています。

波長毎の対数スケールの極座標系式のプロットをアニメーション表示させます。

>	pQ[sca]:=NULL:pQ[ext]:=NULL:pQ[back]:=NULL:

>

pp:=NULL:

>	for lambda from 0.34 by 0.02 to 0.78 do

>	ans:=Mie(radius,lambda,m,density,nangles);

>	print(lambda,ans[3],ans[4],ans[5]);

>	x:=ans[6];

>	parallel:=Vector(nangles):

>	perpen:=Vector(nangles):

>	phasefn:=Vector(nangles):

>	Q[sca]:=ans[3]:Q[ext]:=ans[4]:Q[back]:=ans[5]:

>	for i from 1 to nangles do

>	parallel[i]:=4*abs(ans[2][i])/(x*x*Q[sca]);

>	perpen[i]:=4*abs(ans[1][i])/(x*x*Q[sca]);

>	phasefn[i]:=(parallel[i]+perpen[i])/2;

>

end do:

>	pQ[sca]:=pQ[sca],ans[3]:

>	pQ[ext]:=pQ[ext],ans[4]:

>	pQ[back]:=pQ[back],ans[5]:

>

>	lnpar:=Vector(nangles,0):lnper:=Vector(nangles,0):

>

bias:=0:

>	for i from 1 to nangles do

>	lnpar[i]:=evalf(log10(4*abs(ans[2][i])/(x*x*Q[sca])));

>	lnper[i]:=evalf(log10(4*abs(ans[1][i])/(x*x*Q[sca])));

>	if lnpar[i] < bias then bias:=lnpar[i] end if:

>	if lnper[i] < bias then bias:=lnper[i] end if:

>

end do:

>

>

bias:=4:

>	for i from 1 to nangles do

>	rpar:=lnpar[i]+bias:rper:=lnper[i]+bias:

>	if rpar <0 then rpar:=0 end if:if rper <0 then rper:=0 end if:

>	xpar[i]:=rpar*cos((i-1)/nangles*Pi):

>	ypar[i]:=rpar*sin((i-1)/nangles*Pi):

>	xper[i]:=rper*cos((i-1)/nangles*Pi):

>	yper[i]:=rper*sin((i-1)/nangles*Pi):

>

end do:

>

>	p1:=pointplot({seq([xpar[k],ypar[k]],k=1..nangles),seq([xpar[k],-ypar[k]],k=1..nangles)},color=red,symbol=cross):

>	p2:=pointplot({seq([xper[k],yper[k]],k=1..nangles),seq([xper[k],-yper[k]],k=1..nangles)},color=blue,symbol=cross):

>

c1:=circle([0,0], 1, color=black): c2:=circle([0,0], 2, color=black): c3:=circle([0,0], 3, color=black): c4:=circle([0,0], 4, color=black): t1 := textplot([sqrt(2),sqrt(2),`0.01`],color=green): t2 := textplot([1.5*sqrt(2),1.5*sqrt(2),`0.1`],color=green): t3 := textplot([2*sqrt(2),2*sqrt(2),`1`],color=green): t4 := textplot([0.5*sqrt(2),0.5*sqrt(2),`0.001`],color=green): t5 := textplot([0,0,`0.0001`],color=green):

>	pp:=pp,display({p1,p2,c1,c2,c3,c4,t1,t2,t3,t4,t5},scaling=constrained); end do:

アニメーション表示します。

>	display(pp,insequence=true);

Q[sca] と Q[back] の波長依存のプロットをします。

>	p1:=listplot([pQ[sca]],color=red):p2:=listplot([pQ[back]],color=blue):display([p1,p2]);

粒子のサイズを 1.1 ミクロンに変更して同じ計算を行います。

>	diameter:=1.1:

>	radius:=diameter/2;

>	pQ[sca]:=NULL:pQ[ext]:=NULL:pQ[back]:=NULL:

>

pp:=NULL:

>	for lambda from 0.34 by 0.02 to 0.78 do

>	ans:=Mie(radius,lambda,m,density,nangles);

>	print(lambda,ans[3],ans[4],ans[5]);

>	x:=ans[6];

>	parallel:=Vector(nangles):

>	perpen:=Vector(nangles):

>	phasefn:=Vector(nangles):

>	Q[sca]:=ans[3]:Q[ext]:=ans[4]:Q[back]:=ans[5]:

>	for i from 1 to nangles do

>	parallel[i]:=4*abs(ans[2][i])/(x*x*Q[sca]);

>	perpen[i]:=4*abs(ans[1][i])/(x*x*Q[sca]);

>	phasefn[i]:=(parallel[i]+perpen[i])/2;

>

end do:

>	pQ[sca]:=pQ[sca],ans[3]:

>	pQ[ext]:=pQ[ext],ans[4]:

>	pQ[back]:=pQ[back],ans[5]:

>

>	lnpar:=Vector(nangles,0):lnper:=Vector(nangles,0):

>

bias:=0:

>	for i from 1 to nangles do

>	lnpar[i]:=evalf(log10(4*abs(ans[2][i])/(x*x*Q[sca])));

>	lnper[i]:=evalf(log10(4*abs(ans[1][i])/(x*x*Q[sca])));

>	if lnpar[i] < bias then bias:=lnpar[i] end if:

>	if lnper[i] < bias then bias:=lnper[i] end if:

>

end do:

>

>

bias:=4:

>	for i from 1 to nangles do

>	rpar:=lnpar[i]+bias:rper:=lnper[i]+bias:

>	if rpar <0 then rpar:=0 end if:if rper <0 then rper:=0 end if:

>	xpar[i]:=rpar*cos((i-1)/nangles*Pi):

>	ypar[i]:=rpar*sin((i-1)/nangles*Pi):

>	xper[i]:=rper*cos((i-1)/nangles*Pi):

>	yper[i]:=rper*sin((i-1)/nangles*Pi):

>

end do:

>

>	p1:=pointplot([seq([xpar[k],ypar[k]],k=1..nangles),seq([xpar[k],-ypar[k]],k=1..nangles)],color=red,symbol=cross):

>	p2:=pointplot([seq([xper[k],yper[k]],k=1..nangles),seq([xper[k],-yper[k]],k=1..nangles)],color=blue,symbol=cross):

>

c1:=circle([0,0], 1, color=black): c2:=circle([0,0], 2, color=black): c3:=circle([0,0], 3, color=black): c4:=circle([0,0], 4, color=black): tt:=textplot([3.2,3.2,cat(`λ＝`,convert(lambda,string))]): tt:=tt,textplot([3.2,3.6,cat(`直径＝`,convert(diameter,string))]): t1 := textplot([sqrt(2),sqrt(2),`0.01`],color=green): t2 := textplot([1.5*sqrt(2),1.5*sqrt(2),`0.1`],color=green): t3 := textplot([2*sqrt(2),2*sqrt(2),`1`],color=green): t4 := textplot([0.5*sqrt(2),0.5*sqrt(2),`0.001`],color=green): t5 := textplot([0,0,`0.0001`],color=green):

>	pp:=pp,display({p1,p2,c1,c2,c3,c4,tt,t1,t2,t3,t4,t5},scaling=constrained); end do:

アニメーション表示します。

>	display(pp,insequence=true);

サイズの変化

今度は波長を固定し粒子の半径を変えたときの計算をしてみましょう。

波長は 632 ナノメータです。

>	lambda:=0.632:

>	Q[sca]:=NULL:Q[ext]:=NULL:Q[back]:=NULL:

>	for radius from 0.8 by 0.05 to 1.6 do

>	ans:=Mie(radius,lambda,m,density,nangles);

>	print(radius,ans[3],ans[4],ans[5]);

>	Q[sca]:=Q[sca],ans[3]:

>	Q[ext]:=Q[ext],ans[4]:

>	Q[back]:=Q[back],ans[5]:

>

end do:

この結果（ 0.8 から 1.6 ミクロンまで 0.5 ミクロンづつ変化させた場合の Q[sca] と Q[back] ) のプロットです。

>	p1:=listplot([Q[sca]],color=red):p2:=listplot([Q[back]],color=blue):display([p1,p2]);

>