GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04 LTS desktop amd64
TensorFlow v1.2.1
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)
scipy v0.19.1
geopandas v0.3.0
MATLAB R2017b (Home Edition)
ADDA v.1.3b6
gnustep-gui-runtime v0.24.0-3.1
Numpy実装による球粒子の光散乱コードPyMieScatt関連。
Mie_ab()
https://github.com/bsumlin/PyMieScatt/blob/master/PyMieScatt/Mie.py
http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_ab
Reference
- Bohren C.F. and D. R. Huffman, "Absorption and Scattering of Light by Small Particles", John Wiley, New York, NY (1983)
- PDF: MATLAB Functions for Mie Scattering and Absorption, Research Report No. 2002-08, June 2002 by Christian Mätzler
Bohren and Huffman (1983), Eq. (4.88)
a_n = \frac{[D_n(mx)/m + n/x]\psi_n(x) - \psi_{n-1}(x)}{[D_n(mx)/m + n/x]\xi(x) - \xi_{n-1}(x)}
a_n = \frac{[mD_n(mx) + n/x]\psi_n(x) - \psi_{n-1}(x)}{[mD_n(mx) + n/x]\xi(x) - \xi_{n-1}(x)}
Bohren and Huffman (1983), Eq. (4.89)
D_{n-1} = \frac{n}{\rho} - \frac{1}{D_n + n/\rho}
実装と式の関係
-
px: $\psi_n(x)$
-
p1x: $\psi_{n-1}(x)$
-
gsx: $\xi_n(x)$
-
gs1x: $\xi_{n-1}(x)$
-
da: $[D_n(mx)/m + n/x]$
-
db: $[mD_n(mx) + n/x]$
-
ch: ???
-
ch1x: ???
Bohren and Huffman (1983), p101
The scattering coefficients (4.53) can be simplified somewhat by introducing the Ricatti-Bessel functions:
\psi_n(\rho) = \rho j_n(\rho), \xi_n(\rho) = \rho h_n^{(1)}(\rho)
Bohren and Huffman (1983), p183
The Ricatti-Bessel functions $\chi_n(x)$ is $-zy_n(z)$.
実装からの式推測
実装およびMätzler(2002)のp3より
sx = \sqrt{\frac{\pi}{2x}}x
px = \sqrt{\frac{\pi}{2x}} xj_\nu(x)
chx = -\sqrt{\frac{\pi}{2x}} xy_\nu(x)
( https://en.wikipedia.org/wiki/Bessel_function#Definitions )
j: 第1種球ベッセル関数
y: 第2種球ベッセル関数
p1x: sin(x)を先頭の要素としてpxの前に加えた
ch1x: cos(x)を先頭の要素としてchxの前に加えた
係数が式変更されて実装されているので、文献との齟齬が生じる。
gsx
(追記: 以下の式は虚数部の符号に関して間違っていた）。
gsx = \sqrt{\frac{\pi}{2x}}xj_\nu(x) - i\sqrt{\frac{\pi}{2x}}xy_\nu(x)
gsxはHankel functionのようでもあるが、Mätzler(2002)の式(4.13)と比べた時、虚数部$iy_n(z)$の符号がプラスマイナスで異なる。
別途、PyMieScattのgsxはBohren and Huffman(1983) p481 line132のXIに対応するようだ。
XI = CMPLX(APSI, -CHI)
ここでも虚数の符号はマイナスであり、PyMieScattの実装と整合する。
数値計算ソフトで複素屈折率をm=n+ikと定義する場合とm=n-ikと定義する場合があるが、それと関係するのだろうか。
PyMieScattではこちらのMieQ()で複素屈折率はm=n+ikとしている。
Mie_cd()を見ていると以下がある。こちらの方がHankel functionのように思える。
hx = jnx+(1.0j)*yx
gsxというのは$\xi$のことかもしれない。
Mätzler(2002)の(4.13)の下の式より
...; \xi_\nu(x) = x h_\nu^{(1)}(x) ... (p.101, 183)
Hankel functionをMätzler(2002)のp3の式(4.13)に基づいてj()とy()で記載すると
\xi_\nu(x) = x (j_\nu(x) + iy_\nu(x))
係数$\sqrt{\frac{\pi}{2x}}$を付けると、gsxの実装と合うようだ。
虚数部の符号の不一致はある。
再度、Mie_ab()の定義を確認する。
...
px = sx*jv(nu,x)
...
chx = -sx*yv(nu,x)
...
gsx = px-(0+1j)*chx
chxの定義がマイナスであり、それにマイナスをかけてgsxの虚数部として定義しているので、虚数部の値はプラスということだ。
gsx (正しい)
虚数部の符号を再確認して、以下のようだ。
gsx = \sqrt{\frac{\pi}{2x}}xj_\nu(x) + i\sqrt{\frac{\pi}{2x}}xy_\nu(x)
以下としても表すことができる。
gsx = \sqrt{\frac{\pi}{2x}}x (j_\nu(x) + iy_\nu(x))
元となる式と実装の間に余分な変換があり、読み取りにくい。
プログラムの例としては良くない。
以下のようにすると分かりやすい。
coef = np.sqrt(np.pi/2.0/x)
gsx = coef*x*jv(nu,x) + (0+1j)*coef*x*yv(nu,x)
chx
Mätzler(2002)のp3より
..., \xi_n(z) = zh^{(1)}_n(z),...(p101., 183)
\chi_n(z) = -zy_n(z)
h_n^{(1)}(z) = j_n(z) + iy_n(z) ... (4.13)
上記より
zh^{(1)}_n(z) = zj_n(z) + izy_n(z)
= zj_n(z) - i\chi_n(z)
以下のようだ。
chx: $\chi(x)$