PyMieScatt > Mie_ab() > 実装と式の関係, 虚数部の符号の実装が分かりにくい

```txt:動作環境
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](https://omlc.org/software/mie/maetzlermie/Maetzler2002.pdf) by Christian Mätzler


### Bohren and Huffman (1983), Eq. (4.88)

```math
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)}
```

```math
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)

```math
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:

```math
\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より

```math
sx = \sqrt{\frac{\pi}{2x}}x
```

```math
px = \sqrt{\frac{\pi}{2x}} xj_\nu(x)
```

```math
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
(追記: 以下の式は虚数部の符号に関して間違っていた）。

```math
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では[こちら](http://pymiescatt.readthedocs.io/en/latest/forward.html#functions-for-single-particles)のMieQ()で複素屈折率はm=n+ikとしている。

Mie_cd()を見ていると以下がある。こちらの方がHankel functionのように思える。

```
hx = jnx+(1.0j)*yx
```

gsxというのは$\xi$のことかもしれない。
Mätzler(2002)の(4.13)の下の式より

```math
...; \xi_\nu(x) = x h_\nu^{(1)}(x) ... (p.101, 183)
```

Hankel functionをMätzler(2002)のp3の式(4.13)に基づいてj()とy()で記載すると

```math
\xi_\nu(x) = x (j_\nu(x) + iy_\nu(x))
```

係数$\sqrt{\frac{\pi}{2x}}$を付けると、gsxの実装と合うようだ。

虚数部の符号の不一致はある。

再度、Mie_ab()の定義を確認する。

```py
...
px = sx*jv(nu,x)
...
chx = -sx*yv(nu,x)
...
gsx = px-(0+1j)*chx
```

chxの定義がマイナスであり、それにマイナスをかけてgsxの虚数部として定義しているので、虚数部の値はプラスということだ。


### gsx (正しい)

虚数部の符号を再確認して、以下のようだ。

```math
gsx = \sqrt{\frac{\pi}{2x}}xj_\nu(x) + i\sqrt{\frac{\pi}{2x}}xy_\nu(x)
```

以下としても表すことができる。

```math
gsx = \sqrt{\frac{\pi}{2x}}x (j_\nu(x) + iy_\nu(x))
```

元となる式と実装の間に余分な変換があり、読み取りにくい。
プログラムの例としては良くない。

以下のようにすると分かりやすい。

```py
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より

```math
..., \xi_n(z) = zh^{(1)}_n(z),...(p101., 183)
```

```math
\chi_n(z) = -zy_n(z)
```

```math
h_n^{(1)}(z) = j_n(z) + iy_n(z) ... (4.13)
```

上記より

```math
zh^{(1)}_n(z) = zj_n(z) + izy_n(z)
```

```math
= zj_n(z) - i\chi_n(z)
```

以下のようだ。

chx: $\chi(x)$