FDTD法の基礎

# 0. はじめに
## 0.1 まえがき
今までは偏微分方程式の数値解析を研究していたのだが、突然FDTD法を利用する必要が出てきた。
これまではどちらかというと静電界で現れる方程式の研究だったため、FDTD法のような(動)電磁界解析は触れていなかった。

これまでの研究のおかげで、数値計算の部分はすんなりと理解できるのだが、電磁気学の知識が若干乏しいのもあって、たまに詰まったりもする。
そこで、今回は自分のための知識整理、という意味で、ここにFDTD法の基礎を書き留めようと思う。

## 0.2 アウトライン
この記事では、1次元FDTD法の実装を目標とする。
2次元や3次元のほうがアニメーションで可視化したときにグラフィカルで面白いかもしれないが、
1次元のほうが、可視化したときの電界や磁界の動きがイメージしやすいため、1次元を選んだ。

といっても、実際に利用するときは3次元がほとんどであるため、3次元の定式化について簡単に説明したあと、プログラムは1次元で実装する、という流れにしようと思う。
なお、プログラムはpythonで作る。

また、**吸収境界条件**や**フーリエ変換**など、FDTD法を学習するにあたって欠かせないテクニックも簡単に解説していく。

電磁界解析のイメージが沸かない人は、3.4節の解析結果の図を軽く眺めて見てから読んだほうが良いかもしれない。

実装したプログラムは以下に置いている。
[atily17/FDTD1d(github)](https://github.com/atily17/FDTD1d)

# 1. FDTD法のアルゴリズム
## 1.1 FDTD法とは
**FDTD法**は、**Finite Difference Time Domain Method**の略で、日本語化すると、**有限差分時間領域法**である。
「有限差分」は簡単に言うと[差分法](https://qiita.com/atily17/items/b4e0976867e81a181709)であることを意味する。微分を差に置き換えるという意味だ。
「時間領域」は時間空間について解くことを意味する。対義語は「周波数領域」だ。

「有限差分」を「[有限要素](https://qiita.com/atily17/items/fa8abcc4d778c16fa11a)」に変えたFETD(Finite Element Time Domain)法や、
「時間領域」を「周波数領域」に変えた[FDFD(Finite Difference Frequency Domain)法](https://en.wikipedia.org/wiki/Finite-difference_frequency-domain_method)などもある。ただ、FDTD法と比べるとマイナー。

FDTD法は有限差分法を基礎においているため、定式化は割と簡単だ。
ただ、収束条件など別の注意事項があるため、必要な知識はできる限り説明していきたい。


## 1.2 マクスウェル方程式

まず、簡単に電磁気学の基礎方程式を解説しよう。

電磁気学では、マクスウェル方程式がすべてを表しているため、どんな現象もマクスウェル方程式で説明できる。
マクスウェル方程式は以下の四式だ。

```math
\begin{alignat}{2}
&\displaystyle\nabla \cdot \boldsymbol{D} & =& \rho \\\\
&\displaystyle\nabla \cdot \boldsymbol{B} & =& 0 \\\\
&\displaystyle\nabla \times \boldsymbol{E} + \frac{\partial \boldsymbol{B}}{\partial t} & = & 0 \\\\
&\displaystyle\nabla \times \boldsymbol{H} - \frac{\partial \boldsymbol{D}}{\partial t} & = & \boldsymbol{j}
\end{alignat}
```

ここで、$\boldsymbol{D}$は電束密度、$\boldsymbol{E}$は電界、$\boldsymbol{B}$は磁束密度、$\boldsymbol{H}$は磁界を表し、$\rho$は電荷密度、$\boldsymbol{j}$は電流密度である。[^1]

ただ、実はこれだけでは足りなく、$\boldsymbol{D}$と$\boldsymbol{E}$の関係と、$\boldsymbol{B}$と$\boldsymbol{H}$の関係の二式も必要になる。
材料に物性値である誘電率$\varepsilon$と、透磁率$\mu$の比例関係が成り立つと*仮定される*。これは、実験で与えられた近似値であることを覚えておこう。

```math
\begin{align}
\boldsymbol{D} & \approx \varepsilon \boldsymbol{E} \\\\
\boldsymbol{B} & \approx \mu \boldsymbol{H}
\end{align}
```

あえて、ここではイコールではなく、$\approx$で書くことにした。
材料によってはテンソルになったり、単純な比例関係じゃなかったりもするが、今回は、この比例関係が成り立つことを前提に進めていく。

ただし、真空の場合は真空の誘電率$\varepsilon_0$と真空の透磁率$\mu_0$から

```math
\begin{align}
\boldsymbol{D} & = \varepsilon_0 \boldsymbol{E} \\\\
\boldsymbol{B} & = \mu_0 \boldsymbol{H}
\end{align}
```

が厳密に成り立つことも覚えておこう。

さらに、FDTD法では、電界と電流を関係付ける式である、
$$
\boldsymbol{j} = \sigma \boldsymbol{E}
$$
という式も利用する。この$\sigma$は**導電率**、もしくは**電気伝導率**と呼ばれる。
これは、オームの法則を微分形で表現したものである。[^2]

これで、FDTD法の定式化で必要な式は揃った。

## 1.3 FDTD法の定式化
### 1.3.1 マクスウェル方程式の変形
FDTD法で使うマクスウェル方程式は以下の二つだ。

```math
\begin{align}
\displaystyle\nabla \times \boldsymbol{E} + \frac{\partial \boldsymbol{B}}{\partial t} = 0 \\\\
\displaystyle\nabla \times \boldsymbol{H} - \frac{\partial \boldsymbol{D}}{\partial t} = \boldsymbol{j}
\end{align}
```

ここに、誘電率、透磁率、導電率の関係式である、

```math
\begin{align}
\boldsymbol{D} &= \varepsilon \boldsymbol{E} \\\\
\boldsymbol{B} &= \mu \boldsymbol{H} \\\\
\boldsymbol{j} &= \sigma \boldsymbol{E}
\end{align}
```

を代入しよう。

#### 式(1) マクスウェル方程式の変形1
$$
\displaystyle\nabla \times \boldsymbol{E} = - \mu \frac{\partial \boldsymbol{H}}{\partial t} \tag{1}
$$
#### 式(2) マクスウェル方程式の変形2
$$
\displaystyle\nabla \times \boldsymbol{H} = \varepsilon \frac{\partial \boldsymbol{E}}{\partial t} + \sigma \boldsymbol{E} \tag{2}
$$

この、$\mu$と$\epsilon$と$\sigma$は材料によって決まる既知の値である。
よって、これでたくさんあった変数がたった二つになった。
あとはこれを差分式にするだけだ！　簡単！

### 1.3.2 Yee格子

さて差分近似！　と行きたいところだが、そのまえに**Yee格子**について説明する必要がある。

FDTD法は差分法のため、計算領域を格子分割する必要がある。
Yee格子は簡単に言うと、$\boldsymbol{E}$と$\boldsymbol{H}$で格子のとり方をずらすことで、計算しやすく、そして、高精度にできるというものだ。
そこまで難しくないし、簡単に説明しよう。

$E$の格子は、下のような直方体の格子とする。

![Yee2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/05c2596a-a08d-dfc3-b91d-6d0c62a8cc25.png)


ここで、直方体の大きさは、$\Delta x \times \Delta y \times \Delta z$とする。

直方体の頂点の座標は$(i\Delta x, j\Delta y, k \Delta z)$である。
毎度毎度、$\Delta x$などと書くのは面倒なので、$[i,j,k]$と表現することにする。
例えば、上の図の、濃い緑の直方体の手前左下の頂点は、$[3,2,0]$となる。

そして、この表現を使って、上の図の$E_x$, $E_y$, $E_z$は、
$$
\begin{matrix}
E_x:& [&3,     & \displaystyle3-\frac{1}{2}, &\displaystyle 1-\frac{1}{2} &] \\\\
E_y:& [&\displaystyle 3-\frac{1}{2}, &3,     &\displaystyle 1-\frac{1}{2} &] \\\\
E_z:& [&\displaystyle 3-\frac{1}{2}, &\displaystyle 3-\frac{1}{2}, &1     &]
\end{matrix}
$$
の位置にあることになる。
つまり、それぞれの電界は、
$$
\begin{matrix}
E_x:& [&i,     &\displaystyle j-\frac{1}{2}, &\displaystyle k-\frac{1}{2} &] \\\\
E_y:& [&\displaystyle i-\frac{1}{2}, &j,     &\displaystyle k-\frac{1}{2} &] \\\\
E_z:& [&\displaystyle i-\frac{1}{2}, &\displaystyle j-\frac{1}{2}, &k     &]
\end{matrix}
$$
の位置で定義されることになる。
この格子に関しては、特になんの工夫もない。

次に、磁界$\boldsymbol{H}$の格子について考えよう。
これも、特に違いはない。

![Yee3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/b3567b5a-eab1-e925-f28a-0b0686f76f3b.png)

直方体の大きさは、電界と同じ$\Delta x \times \Delta y \times \Delta z$に設定する。
そして、ここがYee格子の工夫点で、$\boldsymbol{H}$の格子の頂点が、$\boldsymbol{E}$の格子のそれぞれの直方体の中心に来るようにずらして置くのである。
以下に図を書いた[^3]。

![Yee4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/d60f99c7-0029-f604-51df-d473d2d45e3c.png)

このようにずらした格子のため、磁界の位置は、
$$
\begin{matrix}
H_x:& [&\displaystyle i-\frac{1}{2},     &j, & k &] \\\\
H_y:& [& i, & \displaystyle j -\frac{1}{2},     &k &] \\\\
H_z:& [& i, & j, &\displaystyle k -\frac{1}{2}&]
\end{matrix}
$$
となる。

説明の都合上、EとHは以上の座標のとり方になったが、原点を$(\Delta x/2, \Delta y/2, \Delta z/2)$ずらせば、以下のようになる。

#### 式(3): Eの位置
$$
\begin{matrix}
E_x:& [&\displaystyle i+\frac{1}{2},     &j, & k &] \\\\
E_y:& [& i, & \displaystyle j +\frac{1}{2},     &k &] \\\\
E_z:& [& i, & j, &\displaystyle k +\frac{1}{2}&]
\end{matrix} \tag{3}
$$

#### 式(4): Hの位置
$$
\begin{matrix}
H_x:& [&i,     &\displaystyle j+\frac{1}{2}, &\displaystyle k+\frac{1}{2} &] \\\\
H_y:& [&\displaystyle i+\frac{1}{2}, &j,     &\displaystyle k+\frac{1}{2} &] \\\\
H_z:& [&\displaystyle i+\frac{1}{2}, &\displaystyle j+\frac{1}{2}, &k     &]
\end{matrix}\tag{4}
$$

多くの解説書では、このように設定している場合が多いため、今後はこっちで説明することにする。

2次元の場合は、3次元格子から適切な$yz$平面を切り出すことで、

![Yee2d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/aa65f8ff-79b0-0c27-d709-cee6d18d1cdc.png)

となり、
1次元の場合は、2次元格子のある直線を切り出すことで、

![Yee1d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/b4280003-79d5-9fa0-0bf5-e0f62f834094.png)

となる。

また、マクスウェル方程式には時間微分も現れるため、時間についても離散化する必要がある。
時間刻み幅を$\Delta t$とすると、$n$ステップ目の時間は$t = n\Delta t$と表すことができる。
Yee格子では、空間と同様に時間も$\boldsymbol{E}$の計算と$\boldsymbol{H}$の計算において$\Delta t/2$だけずらす。
つまり、$\boldsymbol{E}$の値は、$t = n\Delta$の時間に、$\boldsymbol{H}$の値は、$t=(n+1/2)\Delta t$の時間のものを求める。

時間$t=n\Delta t$における$E_{x}[i,j-1/2,k-1/2]$を$E_{x}^{n}[i,j-1/2,k-1/2]$と、時間$t=(n+1/2)\Delta t$における$H_{x}[i-1/2,j,k]$を$H_{x}^{n+1/2}[i-1/2, j, k]$と表現することにしよう。

まあ、こんなふうに、Yee格子は$\boldsymbol{E}$の格子から空間や時間を半分ずらしたものを$\boldsymbol{H}$の格子と定義されているだけだ。

### 1.3.3 差分式
さて、ついに差分近似だ。

今更説明することでもないかもしれないが、差分近似とは、
(偏)微分の定義
$$
\frac{d f}{d x} = \lim_{\Delta x\rightarrow 0} \frac{f(x+\Delta x) - f(x)}{\Delta x}
$$
において、$\Delta x$が十分に小さいとすると、
$$
\frac{df}{dx} \approx \frac{f(x+\Delta x) - f(x)}{\Delta x}
$$
のように近似できる、というものだ。
$x \rightarrow x-\Delta x/2$と平行移動すると、上の式は、
$$
\frac{df}{dx} \approx \frac{\displaystyle f(x+\frac{\Delta x}{2}) - f(x-\frac{\Delta x}{2})}{\Delta x}
$$
となる。これは、中心差分と呼ばれ、このほうが一般に精度がよくなる。[^4]

さて、もう一度、(変形した)マクスウェル方程式を書き出しておこう。

$$
\displaystyle\nabla \times \boldsymbol{E} = - \mu \frac{\partial \boldsymbol{H}}{\partial t} \tag{1}
$$
$$
\displaystyle\nabla \times \boldsymbol{H} = \varepsilon \frac{\partial \boldsymbol{E}}{\partial t} + \sigma \boldsymbol{E} \tag{2}
$$
まず、式(1)の$x$成分

$$
\frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z} = - \mu \frac{\partial H_x}{\partial t}
$$
のみについて考えてみよう。

(3)と(4)を思い出すと、
$$
\begin{matrix}
E_y:& [& i, &\displaystyle j+\frac{1}{2}, &k &] \\\\
E_z:& [& i, & j, &\displaystyle k +\frac{1}{2} &]
\end{matrix} \tag{3}
$$
$$
\begin{matrix}
H_x:& [& i, &\displaystyle j+\frac{1}{2}, &\displaystyle k+\frac{1}{2} &]
\end{matrix} \tag{4}
$$
だった。

整数と半整数(整数-1/2)がどの座標かに注意して、
$$
\begin{array}{l}
\text{時間} : n\Delta t,\\\\
\text{空間} : [i,j+\frac{1}{2}, k+\frac{1}{2}]
\end{array}
$$
について、全て中心差分をとると、
$$
\begin{array}{ll}
\text{左辺第一項}:&\displaystyle \frac{\partial E_z}{\partial y} \approx \frac{E_z^n [i, j+1, k+\displaystyle\frac{1}{2}] - E_z^n [i, j, k+\displaystyle\frac{1}{2}]}{\Delta y} \\\\
\text{左辺第二項}:&\displaystyle \frac{\partial E_y}{\partial z} \approx \frac{E_y^n [i, j+\displaystyle\frac{1}{2}, k+1] - E_y^n [i, j+\displaystyle\frac{1}{2}, k]}{\Delta z} \\\\
\text{右辺第一項}:&\displaystyle \frac{\partial H_x}{\partial t} \approx \frac{H_x^{n+1/2}[i,j+\displaystyle\frac{1}{2}, k+\frac{1}{2}] - H_x^{n-1/2}[i,j+\displaystyle\frac{1}{2}, k+\frac{1}{2}]}{\Delta t}
\end{array}
$$
と近似できるため、これを代入することで、

```math
\begin{multline}
\frac{E_z^n [i, j+1, k+\displaystyle\frac{1}{2}] - E_z^n [i, j, k+\displaystyle\frac{1}{2}]}{\Delta y} - 
\frac{E_y^n [i, j+\displaystyle\frac{1}{2}, k+1] - E_y^n [i, j+\displaystyle\frac{1}{2}, k]}{\Delta z} \\\\
= - \mu \frac{H_x^{n+1/2}[i,j+\displaystyle\frac{1}{2}, k+\frac{1}{2}] - H_x^{n-1/2}[i,j+\displaystyle\frac{1}{2}, k+\frac{1}{2}]}{\Delta t}
\end{multline}
```

という式ができる。式が長く複雑に見えるが、全然難しいことはしていない。
$y$成分、$z$成分に関しても同様に導出できる。
ぜひ自分で計算してみよう。

次に、式(2)の$x$成分も差分化してみよう。
$$
\frac{\partial H_z}{\partial y}-\frac{\partial H_y}{\partial z} = \varepsilon\frac{\partial E_x}{\partial t} + \sigma E_x
$$

式(3)と式(4)から
$$
\begin{matrix}
E_x:& [& i+\displaystyle\frac{1}{2}, &j, &k &] \\\\
\end{matrix} \tag{3}
$$
$$
\begin{matrix}
H_y:& [&\displaystyle i+\frac{1}{2}, &j,     &\displaystyle k+\frac{1}{2} &] \\\\
H_z:& [&\displaystyle i+\frac{1}{2}, &\displaystyle j+\frac{1}{2}, &k     &]
\end{matrix} \tag{4}
$$
だった。

式(1)のときと半分だけズラした、
$$
\begin{array}{l}
\text{時間} : \displaystyle(n+\frac{1}{2})\Delta t,\\\\
\text{空間} : [i+\frac{1}{2},j, k]
\end{array}
$$
について中心差分を取ると、次のようになる。
$$
\begin{array}{ll}
\text{左辺第一項}:&\displaystyle \frac{\partial H_z}{\partial y} \approx \frac{H_z^{n+1/2} [\displaystyle i+\frac{1}{2}, j+\frac{1}{2}, k] - H_z^{n+1/2} [i + \frac{1}{2}, j - \frac{1}{2}, k]}{\Delta y} \\\\
\text{左辺第二項}:&\displaystyle \frac{\partial H_y}{\partial z} \approx \frac{H_y^{n+1/2} [i + \displaystyle\frac{1}{2}, j, k+\frac{1}{2}] - H_y^{n+1/2} [i+ \displaystyle\frac{1}{2}, j, k-\frac{1}{2}]}{\Delta z} \\\\
\text{右辺第一項}:&\displaystyle \frac{\partial E_x}{\partial t} \approx \frac{E_x^{n+1}[i+\displaystyle\frac{1}{2},j, k] - E_x^{n}[i+\displaystyle\frac{1}{2},j, k]}{\Delta t} \\\\
\text{右辺第二項}:&\displaystyle E_x = E_x^{n+1/2}[i+\frac{1}{2},j, k] \approx \frac{E_x^{n+1}[i+\displaystyle\frac{1}{2},j, k] + E_x^{n}[i+\displaystyle\frac{1}{2},j, k]}{2} 
\end{array}
$$
特筆すべきは右辺第二項だ。
これは、もともと微分じゃないので、差分近似する必要はない。
しかし、Yee格子では電界の値は時間が$n\Delta t$のときしか使わないと決めているので、$(n+\frac{1}{2})\Delta t$が出てきては困る。ということで、前後の時間($n\Delta t$と$(n+1)\Delta t$)の二つの値の平均をとってしまおう！　という発想だ。[^5]

長くなるので、わざわざ式は書き下さないが、これらを式(2)に代入すれば、離散式の完成。
こちらも$y$, $z$成分も同様に求めることができる。
これで、式(1)と合わせて計6個の式が得られることになる。


### 1.3.4 1次元FDTD法の定式化
式が煩雑になるので、ここからは実装すると決めている1次元で考えてみよう。

$\boldsymbol{E}$, $\boldsymbol{H}$の電磁波があった時、進行方向は、
$$
\boldsymbol{E} \times \boldsymbol{H}
$$
の方向になる。
つまり、電磁波が$z$方向に進行していたとして、電界は$E_x$成分のみとすると、磁界は$H_y$成分のみとなる。

よって、一次元FDTDでは、$E_x$, $H_y$以外は0と考えればよい。

例えば、式(1)の$y$成分は、
$$
\frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x} = - \mu \frac{\partial H_y}{\partial t}
$$
だが、$E_z$は0のため、
$$
\frac{\partial E_x}{\partial z} = - \mu \frac{\partial H_y}{\partial t}
$$

となる。[^6]
面倒なので、今後は、$E_x=E$, $H_y=H$と書くことにし、$x$座標、$y$座標はもう必要ないため、$E(z) =E(i\Delta z)$のとき、$E[i]$のように表現することにしよう。
1次元のときのYee格子は

- $E$: 時間 $n\Delta t$, 空間 $[i]$
- $H$: 時間 $(n+1/2)\Delta t$, 空間 $[i+1/2]$

の値が与られる。

上の式(1)の$y$成分は

#### 式(5): 1次元マクスウェル方程式(1)

$$
\frac{\partial {E}}{\partial z} = - \mu \frac{\partial H}{\partial t} \tag{5}
$$

であり、式(2)の$x$成分を考えると、

#### 式(6): 1次元マクスウェル方程式(2)

$$
-\frac{\partial {H}}{\partial z} = \varepsilon \frac{\partial E}{\partial t} + \sigma E \tag{6}
$$

である。
この式(5)と式(6)の二つを利用して、1次元FDTD法を定式化する。
導出は省略しよう。

$$
\frac{E^{n}[i+1] - E^{n}[i]}{\Delta z} = -\mu \frac{H^{n+1/2}[i+\frac{1}{2}] + H^{n-1/2}[i+\frac{1}{2}]}{\Delta t} \\\\
-\frac{H^{n+1/2}[i+1/2] - H^{n+1/2}[1-1/2]}{\Delta z} = \varepsilon \frac{E^{n+1}[i] E^n[i]}{\Delta t} + \sigma\frac{E^{n+1}[i] + E^n[i]}{2}
$$

アルゴリズムにしやすいように整理すると、次のようになる。

#### 式(7): 1次元FDTD(1)

$$
H^{n+1/2}[i+\frac{1}{2}] = H^{n-1/2}[i+\frac{1}{2}]-\frac{1}{\mu}\frac{\Delta t}{\Delta z}(E^{n}[i+1] - E^{n}[i]) \tag{7}
$$

#### 式(8): 1次元FDTD(2)

$$
E^{n+1}[i] = \frac{2\varepsilon - \sigma \Delta t}{2\varepsilon + \sigma \Delta t}E^{n}[i]-\frac{2}{2\varepsilon + \sigma \Delta t}\frac{\Delta t}{\Delta z}(H^{n+1/2}[i+\frac{1}{2}] - H^{n+1/2}[i-\frac{1}{2}]) \tag{8}
$$

これで、1次元FDTD法に必要な式は揃った。

### 1.3.5 アルゴリズム
FDTD法のアルゴリズムは以下のようになる。

1. $E^{0}[i]$と$H^{-1/2}[i+1/2]$を初期条件として与える。
2. 式(8)で、$E^{1}[i]$を求める。
3. 式(7)で、$H^{1+1/2}[i]$を求める。
4. 2,3を繰り返して、$E^{n}[i]$と$H^{n+1/2}[i]$を交互に求めていく。
5. 適当なところで打ち切って終了

基本的に以上の流れで計算できる。
ただし、FDTD法は、格子分割するため、領域は無限に大きくできないことは注意しよう。
例えば、$[0]$から$[I]$までの有限領域で計算しようと思っても、$E^{n+1}[0]$を求めようとした場合、領域外の値$H^{n+1/2}[-1/2]$が必要になってしまう。
そのため、**境界条件**というものを与える必要がある。

例えば、**完全導体**(電界が0になる)で囲まれた空間内で計算しようと思った場合は、$E^{n}[i]$を計算する段階で、$E[0]=0$, $E[I]=0$とすればいいだけだ。

ただ、無限に大きい空間を計算しようとした場合は、ちょっと面倒で、後述する**吸収境界条件**というものを使う必要がある。

3次元の場合も、上と同様の流れで計算していけばいい。

# 2. FDTD法で必要な技術・知識
上記まででFDTD法のアルゴリズムは実装できる。
しかし、FDTD法や電磁界解析には、その他いろんな技術や知識が必要になる。

ここからは、ちょっと雑多な情報になるかもしれないが、全て必須とも言える内容だ。
全部しっかりと理解しておこう。

## 2.1. フーリエ変換

3次元FDTDを解くと、時間$t$, 位置$x,y,z$, 電界$E_x, E_y, E_z$, 磁界$H_x, H_y, H_z$の10個のパラメータの関係が得られる。多い…。
1次元だとしても、時間$t$, 位置$z$, 電界$E(=E_x)$, 磁界$H(=H_y)$の4つのパラメータの関係が得られることになる。

これから紹介するフーリエ変換をすると、更にそこに**周波数**$f$とその**周波数成分の大きさ(スペクトル)**$\mathcal{F}$が出てくる。
しかも、$\mathcal{F}$は複素数のため、より複雑だ。
こんなに多いと、混乱しそうだ。整理しながらフーリエ変換について話していこう。

**フーリエ変換**はFDTD法で得られた結果を**時間領域**から**周波数領域**に変換する手法である。
これはどちらかというと、FDTD法というよりは、FDTD法で得られた結果から、その性質を見て検証する段階で使うものだ。

ところで、時間領域とか周波数領域とか何なんだ、と思う人もいると思う。
ここでは、あまり深いところには触れずに、簡潔に説明する。より深く知りたい人は、色んな所で詳しく解説されているので**フーリエ変換**、**離散フーリエ変換**、**高速フーリエ変換**あたりのキーワードでググろう。

まず、以下のグラフを見てほしい。

(gif動画のため、動かない場合はクリック)
![time(loss).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/bc65e681-7ae3-b586-5992-ff9c20e44f8d.gif)


上図の上のグラフは$z$方向に進む電界である。単純化のため、速度1で進んでいるとしよう。
ここで、観測点(青い点線)の電界の時間変化を見てみよう。これが、下のグラフである。
この下のグラフが時間領域のグラフである。

さて、これを周波数領域について見てみよう。
時間領域のグラフに注目して、(観測点上の)電界$E$は時間$t$の関数のため、$E(t)$と表す。

これを、以下の式でフーリエ変換する。[^7]
$$
\mathcal{F}(f) = \int_{-\infty}^\infty E(t)e^{-j t (2\pi f)}dt
$$
$j$は虚数単位である。ここで、$2\pi f=\omega$とすれば、
#### 式(9): フーリエ変換
$$
\mathcal{F}(\omega) = \int_{-\infty}^\infty E(t)e^{-j t \omega}dt
$$
となる。この$\omega$は**角周波数**と呼ばれる。
角周波数もまとめて周波数と表現している資料も多いのだが、ここでは厳密に言葉を使い分けることにする。つまり、周波数と言ったら$f$だし、各周波数と行ったら$\omega$である。

式(9)は、無限領域で積分しているが、例えば、上のグラフでは、$0.5\le t \le 0.6$以外では0になっている。
実際にフーリエ変換が必要になるときも、大抵は十分な時間$T$がたてば0になると考えることができて、
$$
\mathcal{F}(\omega) = \int_{0}^{T} E(t)e^{-j t \omega}dt
$$
とすることができる。
これを以下のように数値積分することで、コンピュータでも簡単に計算できるようになる。
#### 式(10): 離散フーリエ変換(DFT)
$$
\mathcal{F}(\omega) = \sum_{n=0}^{N} E(n\Delta t)e^{-j n\Delta t  \omega}\Delta t \tag{10}
$$
ここで、$t=n\Delta t, T=N\Delta t$で、$N$は分割数である。
この式を、**離散フーリエ変換**(Discrete Fourier Transform; **DFT**)という。[^8][^9]

それでは、実際にフーリエ変換してみよう。

![fourier.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/b61b22d5-70c0-f532-0e69-d3a66b9d7141.png)


下のグラフがフーリエ変換した結果だ。
え？　こんなの見ても全然意味がわからない？　その意味をこれから説明していく。

まず、電界(上のグラフ)が

$$
E(t) = A_0\cos{(0\times 2\pi)} + A_1\cos{(1\times 2\pi)} + A_2 \cos{(2\times 2\pi)} + \cdots + A_N(\cos{N \times 2\pi}) \\\\
+B_1\sin{(1\times 2\pi)} + B_2\sin{(2\times 2\pi)} + \cdots + B_N\sin{(N\times 2\pi)}
$$

って感じで、三角関数の和で表すことができると仮定しよう。
すると、この$A_n, B_n$(の絶対値)が大きければ大きいほど、「その周波数の波が強い」ことがわかる。
たとえば、上の図だったら、なんとなく周波数$f=8$くらいまでは結構強くて、それ以降は徐々に小さくなっている、という事が見て取れる。つまり、なんとなく低い周波数の寄与が大きいのかな～って事がわかるだろう。
また、よくみると、ある周波数の周期で($15, 25, 35,\cdots$)で強くなっていることがわかる。
これも特性の一つだろう。

前述の通り、スペクトルは複素数の値を取る。
このとき、実部が$\cos$成分の大きさ$A_n$を表し、虚部が$\sin$成分の大きさ$B_n$に対応しているのである。

たとえば、周波数$f=1$のとき、スペクトルはだいたい、$-0.19+6.02j$である。
つまり、$A_1 = -0.19$, $B_1=6.02$ということがわかるのである。
以下のアニメーションを見てほしい。もとの形状に戻っていくことがわかるだろうか？

(gif動画のため、動かない場合はクリック)
![ifourier(loss).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/66397637-1150-f943-d066-dbed16a3139e.gif)

その周波数帯の強さ(振幅)は、$\sqrt{(\text{Re}(\mathcal{F}))^2+(\text{im}(\mathcal{F}))^2}$で表される。$\text{Re}$と$\text{Im}$は実数部と虚数部を表す。
これは、三角関数の合成公式
$$
A_n\cos{(2\pi ft)}+B_n\sin{(2\pi ft)} =\sqrt{A_n^2 + B_n^2} \sin{(2\pi t + \alpha)}
$$
から、理解できる。また、$\alpha$に対応するところは、位相差と呼ばれたりする。
普通、実数部・虚数部よりも、こちらの強さと位相差のほうで特性を見ることが多い。

で、こんな事がわかって何が嬉しいか？　という話だが、電磁波というものは周波数帯で全然違う振る舞いをする。
例えば、電磁波が鉄などの導体にぶつかった時、周波数が高ければ高いほど、急激に減衰する**表皮効果**というものが発生する。
もし、とても薄い鉄の板に、低い周波数から高い周波数まで一様に分布している電磁波があったとすると、高い周波数成分は減衰して消えてしまい、低い周波数成分のみが残るのである。
つまり、外部に電磁波を漏らさないように、鉄で囲っていたとして、もし低い周波数成分が強い！　と言うことがわかると、外に漏れてしまうことが解る。
逆に、高い周波数成分が大きくて、低い周波数成分がそうでもないなら、鉄の板をもっと薄くして軽量化できそう、などいろいろなことが考えられるのだ。

それ以外にも、周波数が重要なのは**アンテナ**だ。
例えば、bluetoothなんかは、2.4GHz帯の電磁波で通信しており、この2.4GHz帯に対応する周波数帯はどこまで届いているだろうか？　ということを計測することができる。

ところで、FDTD法は電界$E^t[i]$や磁界$H^{t+\frac{1}{2}}[i+\frac{1}{2}]$を順次解いていく方法のため、式(10)のDFTをFDTDをときながら計算していくことが可能である。
ここに関しては、実装のところで詳しく話そう。

## 2.2 波源モデル(初期値)
FDTD法を解くには、初期値を与える必要がある。この初期値は波源ともいう。
初期値として強さ1の電界を与えるとなると、どうやって与えればいいだろうか？

最も簡単なのは、一つの適当な要素に1という値を与えて、それ以外を0にする、つまり、**デルタ関数**として与える手法だ。
実際には、格子分割して有限領域で求めるため、上のフーリエ変換で説明したような、矩形パルスになる。

しかし、このように設定してしまうと、分割数が大きくなるにつれて高周波成分が大きくなってしまうし、上のパルスを例に取ると、上述したように$f=15, 25, 35,\cdots$と、周期的に強い成分ができているし、逆に、$f=20, 30, 40,\cdots$あたりはほぼ0になっていることがわかる。
このようなスペクトルのパルスを波源とすると、観測したスペクトルが波源によるものなのか、その他の要素のものなのかがわからず、不便な時がある。

そこで、よく使われる**ガウスパルス**について説明しよう。

ガウスパルスは、ガウス関数
$$
p(\tau) = \exp{\left(-\frac{\tau^2}{2w^2}\right)}
$$
ここで、$w$はパルス幅である[^10]。
$\tau$は光の速度を$c$として、$\tau = t + z/c$で計算できる。$z/c$の項は、速度$c$で進んでいる電磁波が位置$z$のとき、$z=0$よりどれくらい進んでいるかを表す。

$\tau=0$の時、ガウスパルスは最大値となり$p(\tau)=1$になる。

ガウスパルスの便利なところは、フーリエ変換してもガウス関数のままである。

![gaussian.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/3451f299-0ad9-8f79-adb5-eaf5f784fbc8.png)


上の図は、ガウスパルスとそのスペクトル特性だ。
スペクトル特性は、ガウス関数の右半分になっている。
$w$を小さくすればするほど、ガウスパルスはより鋭く、スペクトル特性は、よりなだらかに、高い周波成に伸びていく。
逆に、$w$を小さくすると、ガウスパルスはなだらかになり、スペクトル特性は、低い周波帯のみとなる。

例えば、1GHzオーダーの解析をしたいのであれば、ガウスパルスの幅を狭く鋭くして、高周波成分にまで残せばいいし、逆に1MHzオーダーだったら、ガウスパルスの幅はなだらかにして、余計な高周波成分を含めないようにすればいい。

ガウスパルスの他にも、ガウス関数を微分した微分ガウスパルスや、sin波で変調したガウスパルスなどもよく用いられる。これらは、ガウスパルスと違い、平均が0になるため、直流成分をなくしたい場合によく用いられる。

なお、ガウス関数は、$\tau$がどれだけ大きくなっても、0にならない。
そこで、以下のようにある程度の時刻で打ち切るのが普通である。

$$
p(\tau)=\left\\{
\begin{array}{ll}
\exp{\left(\displaystyle-\frac{\tau^2}{2w^2}\right)} & \text{if }(-N_t w < \tau < N_tw)\\\\
0 & \text{else}
\end{array}
\right.
$$

ここで、$N_t=2$のときは、$p(\tau)\simeq 0.135=(13.5 \\%)$, $N_t=3$のときは、$p(\tau)\simeq 0.0111(=1.11 \\% )$となる。
この$N_t$を適当に決めて、打ち切ればいい。

## 2.3 吸収境界条件
次に、**吸収境界条件**について解説する。

簡単のため、以下の図の1次元Yee格子をもとに説明していく。

![Yee1d_I.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/67bd60f0-69b4-0393-71a2-7bddf2b8510b.png)


前述したとおり、FDTD法は、有限領域を解析するため、領域の境界には**境界条件**を設定する必要がある。

境界条件には、以下のように、境界をある値と固定するというものだ。
$$
E^{n}[0]=A, E^{n}[I] =B
$$
このような条件をディリクレ境界条件という。
$A=B=0$のとき、これは完全導体に囲まれていることを意味する。

しかし、電磁界解析では、無限に広いと仮定して解析したいことも多い。
このとき、使うのが吸収境界条件だ。

一言でいうと、吸収境界条件は「電磁界が反射しないような境界」である。
この吸収境界条件はいろいろな手法が提案されており、**Murの吸収境界**と**BerengerのPML吸収境界**が有名だ。
**Murの吸収境界**は実装が簡単だが精度に劣る。一方、**BerengerのPML吸収境界**は高精度だが、実装が難しく、また少しメモリも必要になる。
他にも、LiaoやHigdonの吸収境界などがあるが、上2つと比べるとあまり有名じゃないのかもしれない。詳しく知りたい人は、参考文献の「FDTD法による電磁界およびアンテナ解析」を読もう。

今回は、Murの吸収境界条件だけ説明する。
一次元FDTDではこれで十分だ。

### 2.3.1 Murの吸収境界条件
電磁波というのは、光速$c$の速度で進んでいる。
それならば、境界の外側へそのまま光速$c$で進んでいく式を立てる、というのが、Murの境界条件のアイデアだ。

ところで、FDTD法の定式化をしたマクスウェル方程式に速度$c$が出てこないのに、突然$c$が出てくるのに不安を覚える人もいるかも知れない。これについては、マクスウェル方程式をいじると、$1/\sqrt{\varepsilon_0 \mu_0}$の速度で電磁波が伝搬しているというのがわかる。
そして、この値が光速と同じだ！　ってことから、
$$
c = \frac{1}{\sqrt{\varepsilon_0 \mu_0}}
$$
ってことが解ったという歴史がある。
そのため、直接は見えないが、真空の誘電率と真空の透磁率の値で関連付けられており、ちゃんとまくすうぇう方程式に則ってるので安心しよう。

まず、上の一次元Yee格子の、$E^n[I]$側について考えよう。
もし、$E$が(形状を変えずに)速度$c$で$z$方向に進んでいるとすれば、
$$
E(z,t)=E(z-ct)
$$
と表すことができる。
なぜか？　まず、位置が$z$, 時間が$t_0$の電界は$E(z,t_0)$と表せる。
それから$t$たったとき、$E(z,t_0)$は$z$方向に$ct$進んでるのだから、
$$
E(z+ct,t_0+t) = E(z,t_0)
$$
と表せるだろう。さらに、$z+ct\rightarrow z$, $t_0+ t \rightarrow t$と置き直そう。すると、
$$
E(z,t) = E(z - ct, t_0)
$$
となるが、$t_0$は定数のため、結局、$E(z,t) = E(z-ct)$となる。

つまり、$E(z,t) = E(z-ct)$を境界条件としてしまえば、形状を変えず、反射もせず電磁波は境界を通り抜けていくことになる。

上の式の両辺を$z$で偏微分してみよう。
$$
\begin{array}{ll}
\displaystyle\frac{\partial E(z,t)}{\partial z} & =\displaystyle\frac{\partial E(z-ct)}{\partial z}\\\\
&\displaystyle= \frac{\partial (z-ct)}{\partial z} \frac{\partial E(z-ct)}{\partial (z-ct)} \\\\
&\displaystyle= \frac{\partial E(z-ct)}{\partial (z-ct)}
\end{array}
$$
次に、$t$でも偏微分する。
$$
\begin{array}{ll}
\displaystyle\frac{\partial E(z,t)}{\partial t} & =\displaystyle\frac{\partial E(z-ct)}{\partial t}\\\\
&\displaystyle= \frac{\partial (z-ct)}{\partial t} \frac{\partial E(z-ct)}{\partial (z-ct)} \\\\
&\displaystyle= -c\frac{\partial E(z-ct)}{\partial (z-ct)}
\end{array}
$$
二つの式の右辺が似ていることに注目して、代入すると、
$$
\displaystyle\frac{\partial E(z,t)}{\partial t} = -c\displaystyle\frac{\partial E(z,t)}{\partial z} 
\iff \displaystyle\frac{\partial E(z,t)}{\partial z} +\frac{1}{c}\displaystyle\frac{\partial E(z,t)}{\partial t} = 0
$$

という式が得られる。

あとは差分近似して、
$$
\frac{E^{n-\frac{1}{2}}[I] - E^{n-\frac{1}{2}}[I-1]}{\Delta z} + \frac{1}{c}\frac{ \displaystyle E^{n}[I-\frac{1}{2}] - \displaystyle E^{n-1}[I-\frac{1}{2}] }{\Delta t} = 0
$$
となり、FDTD法の導出のときと同じく、Yee格子に合わせるために、平均を取り、

```math
\begin{multline}
\frac{1}{\Delta z} \left( \frac{E^{n}[I] + E^{n-1}[I]}{2} - \frac{E^{n}[I-1] + E^{n-1}[I-1]}{2} \right) \\\\
+ \frac{1}{c\Delta t}\left( \frac{E^n[I]+E^n[I-1]}{2} - \frac{E^{n-1}[I]+E^{n-1}[I-1]}{2} \right) = 0
\end{multline}
```

となり、整理すると、
#### 式(11): Murの吸収境界(+z側)
$$
E^n[I] =  E^{n-1}[I-1] + \frac{c\Delta t - \Delta z}{c\Delta t + \Delta z}\left( E^n[I-1] - E^{n-1}[I] \right) \tag{11}
$$
となる。これがMurの吸収境界である。
$E^{n}[0]$の方は、電磁波の進行方向が逆のため、
$$
E(z,t)=E(z+ct)
$$
として計算すれば良い。
計算過程は省略するが(自分で計算してみよう)、結果は、式(11)の$I\rightarrow 0$, $I-1 \rightarrow 1$と書き換えた、

#### 式(12): Murの吸収境界(-z側)
$$
E^n[0] =  E^{n-1}[1] + \frac{c\Delta t - \Delta z}{c\Delta t + \Delta z}\left( E^n[1] - E^{n-1}[0] \right) \tag{12}
$$
となる。

Murの吸収境界は、一次元の問題の時、(差分近似の部分を除けば)反射係数は0になり、理想的な吸収境界となる。これが、上で一次元の場合はMurの境界条件で十分といった理由だ。
しかし、二次元以上の問題のときは、境界に対して斜めに入射してくる電磁波も無数に存在する。
このような問題では、$E=E(z-ct)$は成り立たず、精度が落ちてしまう。
その誤差が無視できない場合には、PML吸収境界など、他の吸収境界を使った方が良い。

## 2.4 Courantの安定条件(収束条件)
FDTD法は$\Delta z$を小さくすると、$\Delta t$も小さくしないと発散して正しくない解になる可能性がある。
結論から言おう。1次元FDTDのCourantの安定条件は、
#### 式(13): Courantの安定条件(1次元)
$$
c\Delta t \le \Delta z \tag{13}
$$
である。ここで$c$は速度を表す。

この式の意味を簡単に説明しよう。

電磁波が速度$c$で進んでいるとする。
すると、$t=n \Delta t$のとき、$z=0$の電界$E^n[0]$は、$t=(n+1)\Delta t$は$z=c\Delta t$まで進んでいることになる。

ここで、例えば、Courantの安定条件を満たさない、$c \Delta t = 2\Delta z$と設定したとする。
すると、速度$c$は
$$
c=(2\Delta z)/\Delta t
$$
となる。

以下にイメージしやすいように図を用意した。

![Courant.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/d98eba49-090c-eee0-6b7c-14010c8d2328.png)


Yee格子で表すと、$\Delta t$秒で、$z=[0]$のデータは$z=[2]$まで進むことになる(図の白抜きの赤丸の部分)。

しかし、1次元FDTDの式(7)と式(8)を思い出すと、$E^{n}[0]$を利用して求めることができる値は、$H^{n+1/2}[1/2]$(図の$t=(n+1/2)\Delta t$の行の赤丸)である。
さらに、$H^{n+1/2}[1/2]$で求めることができる値は、$E^{n}[1]$(図の$t=(n+1)\Delta t$の行の赤丸)である。

つまり、どう頑張っても、式(7)と式(8)では、$E^{n}[0]$にある値を$E^{n}[2]$まで送ることができないのだ。

そのため、Courant安定条件$c \le \Delta z/\Delta t$を満たさなければおかしな値が得られることになるのだ。
以上は感覚的な説明で、もう少し厳密な証明もあるが、ここでは説明しない。

なお、1次元FDTDではなく、$d$次元FDTDでは

#### 式(14): Courantの安定条件(d次元)
$$
c\Delta t \le \frac{\Delta}{\sqrt{d}} \tag{14}
$$
である。ただし、上式は立方(正方)格子の場合で、$\Delta=\Delta x = \Delta y = \Delta z$である。

# 3. 実装
説明がとても長くなってしまったが、実装に入ろう。

作るプログラムは何度も行っているように一次元の電磁界だ。
今回は、$z$方向に進む平面波($xy$平面に一様な波)がどのように伝搬していくか？
というシミュレーションをしてみよう。

ソースコードは以下のURLにおいた。

必要なパッケージはnumpy, scipy, matplotlibがあれば動くと思う。

- [atily17/FDTD1d(github)](https://github.com/atily17/FDTD1d)

以下では、このソースの中から抜粋して重要なところを解説していく。

## 3.1 問題設定
上述の通り、平面波のシミュレーションをする。
パルスは、わかりやすくガウスパルスとしよう。
以下のようなグラフの電界と磁界を右に移動させていくイメージだ。

![EH0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/deb91a40-8f89-9c76-8ee2-3298186a5bb8.png)

ガウスパルスは大きさ、パルス幅などを設定できるようにし、
境界条件も、ディリクレ境界条件と吸収境界条件の二つを与えられるようにしよう。

また、媒質の位置やパラメータも与えられるようにする。
これで結構面白い波の波形が得られる。

## 3.2 Yee格子
まず、Yee格子を作っていこう。

Yee格子の前に、媒質(Medium)クラスを作ることにした。

```python
class Medium:
    def __init__(self):
        self.setMediumInit()
        self.left = 0
        self.right = 0

    def setRegion(self, left, right):
        self.left=left
        self.right=right

    def setMediumInit(self):
        self.eps = constants.epsilon_0
        self.mu = constants.mu_0
        self.sgm = 0

    def setMedium(self, reps, rmu, sgm):
        self.setMediumInit()
        self.eps = reps*self.eps
        self.mu = rmu*self.mu
        self.sgm = sgm
```

これらから、媒質の各パラメータを設定できる。
簡単に説明すると、`self.eps`は真空の誘電率を初期値として入れている。透磁率も同様だ。
そして、`setMedium`関数に、`reps`(比誘電率$\varepsilon_r$), `rmu`(比透磁率$\mu$), `sgm`(導電率)を引数にとって、誘電率$\varepsilon=\varepsilon_r \varepsilon_0$, 透磁率$\mu = \mu_r \mu_0$と、導電率$\sigma$の値を設定できる。

それではYeeGridクラスを作ろう。
メンバ変数は以下の通りだ。

```python
class YeeGrid:
    def __init__(self, directory="output"):
        self.domain = 0 # 領域長さ
        self.n = 0      # 分割数
        self.dz = 0     # 空間分割幅
        self.dt = 0     # 時間分割幅
        self.mediums = []   # 媒質
        self.e_element = np.array([])   # E要素の座標(numpy)
        self.h_element = np.array([])   # H要素の座標(numpy)
```

まず、解析する幅だが、左の境界の幅は$z=0$で固定して、右の境界だけを設定するようにした。
これは、初期値なのでとりあえず全部に0を入れている。
媒質はリストで、Mediumオブジェクトを追加していく。

次に、グリッドを作る関数は以下の通りだ。

```python
class YeeGrid:
    def setGrid(self, domain, N, ratio=0.99):
        self.domain = domain
        self.n = N
        self.dz = self.domain/self.n
        self.e_element=np.linspace(0,self.domain, N+1)
        self.h_element=np.linspace(0+self.dz/2,self.domain-self.dz/2, self.n)
        c = constants.c
        self.dt = self.dz/c*ratio
```

引数の`domain`は領域の範囲で、`N`は分割数である。
これで、`self.e_element`には`N+1`個の要素が、`self.h_element`には`N`個の要素ができる。


`setGrid`関数で`self.e_element`と`self.h_element`に各格子の要素の座標を入れる。
`self.e_element=np.linspace(0,self.domain, N+1)`と
`self.h_element=np.linspace(0+dz/2,self.domain-dz/2, N)`の部分がそれだ。
磁界$H$の格子は$E$から半分ずらすため、格子間隔`dz`の半分をずらしている。
なお、`self.e_element`の要素数が`self.h_element`の要素数より1多いのも、Yee格子の図を見れば容易にわかる。

たとえば、`setGrid(1, 10)`なら、

```
self.e_element = [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
self.h_element = [0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95]
```

となる。
また、ここでCourantの安定条件から、$\Delta t$も決定した。
式は、以下の通りだ。
$$
\Delta t = r\frac{\Delta z}{c}
$$
理屈上、$\Delta t = \Delta z / c$でもいいのだが、丸め誤差によって発散する可能性が高いため、$r=0.99$をかけることで、1%だけ小さくした。

以下の`addMedium`関数は、空間に媒質を追加する関数だ。

```python
class YeeGrid:
    def addMedium(self, left, right, reps, rmu, sgm):
        m = Medium()
        m.setRegion(left, right)
        m.setMedium(reps, rmu, sgm)
        self.mediums.append(m)
```

特に説明することはないと思う。

## 3.3 FDTDソルバー
それでは、FDTDのメイン箇所を書いていこう。

### 3.3.1 初期条件
まず、初期条件を作っていこう。
Pulseクラスを作った。メンバ変数は以下の通りだ。

```python
class Pulse:
    def __init__(self):
        self.w = 100e-12    # パルス幅
        self.threshold = 4  # 打ち切る位置
        self.E = 1          # パルスの大きさ
```

`self.threshold`はパルスの式
$$
p(\tau)=
\left\\{
\begin{array}{ll}
\exp{\left(\displaystyle-\frac{\tau^2}{2w^2}\right)} & \text{if }(-N_t w < \tau < N_tw)\\\\
0 & \text{else}
\end{array}
\right.
$$
の$N_t$に値する。

この計算は、YeeGridオブジェクトを引数とした。

```python
class Pulse:
    def calc(self,grid):
        c = constants.c
        mu0 = constants.mu_0
        eps0 = constants.epsilon_0
        impedance = np.sqrt(mu0/eps0)
        z0 = c*self.w*self.threshold+2*grid.dz
        taues = (grid.e_element-z0)/c
        tauhs = grid.dt/2 + (grid.h_element-z0)/c
        return {"E": self.E*self._gaussian(taues, self.w), 
                "H":self.E
/impedance*self._gaussian(tauhs, self.w)}
```

ここで、注意すべきは、$p(\tau)$の$\tau$だ。$\tau = t + z/c$と書いたが、まず、$z$は左境界を0にしたため、そのままでは、ガウス関数の右半分しか計算できない。
そのため、$p(\tau)=0$となるように、$z$を$N_tw$だけずらす必要がある。
そこで、$z \rightarrow z-N_tw$と設定する。
上のプログラムでは、安全のためさらに$2\Delta z$分ずらしている。

`self._gaussian`関数は、上の式を再現しているだけで、ここでは載せない。

### 3.3.2 1次元FDTDの式(7), (8)を実装する。
次に、FDTDを実装しよう。

以下の流れになっている。

```python
class Solver:
# __init__などは省略
    def calc(self):
        self.calcInit() # 初期値
        self.save(0)    # そのときの電界、磁界をセーブ
        self.time = 0   # 今の時間
        self.calcFourier()  # フーリエ変換の計算
        for t in range(1, self.NTime+1):
            self.nextstep(t)    # 次のステップへ
            self.save(t)        # そのときの電界、磁界をセーブ
            if t % 10 == 0:
                print("iter="+str(t)+", ", end="")
                if self.judgeConvergence(bPrint=True):  # 収束判定
                    break
```

文字で書くと次のとおりだ。
1. まず、初期値を決める。
2. FDTD法の計算式で次のステップへ
3. 電界、磁界の値のセーブ
4. ループ10回につき収束判定をする。収束しているなら、ループから抜ける。

フーリエ変換とセーブの部分は後々説明する。

次に、2.のところを細かく見てみよう。

```python
class Solver:
    def nextstep(self, t):
        preE = self.E.copy()    # 電界Eのコピー
        self.calcH()            # Hを式(7)から計算
        self.time += self.grid.dt/2   # dt/2 時間足す
        self.htimes.append(self.time) # 計算に使った時間を保管
        self.calcFourier("H")   # Hのフーリエ変換
        self.calcE()            # Eを式(8)から計算
        self.calcBorder(preE)   # Eの境界条件を計算
        self.time += self.grid.dt/2   # dt/2 時間足す
        self.etimes.append(self.time) # 計算に使った時間を保管
        self.calcFourier("E")   # Eのフーリエ変換
```

`self.E`と`self.H`が電界、磁界の値が入ったnumpy配列である。
重要なのは、`calcH()`と`calcE()`, そして`calcBorder()`で、これが式(7)と式(8)と境界条件に値する。
最初に電界`E`をコピーしているのは、Murの吸収境界条件、式(11), 式(12)で、「一つ前の時間の電界」の値が必要だからだ。

`calcH()`は以下のようになる。

```python
class Solver:
    def calcH(self):
        mu = self.properties["mu"]
        dt = self.grid.dt
        dz = self.grid.dz
        c1 = 1/mu*dt/dz
        self.H = self.H - c1 * (self.E[1:] - self.E[:len(self.E)-1]) # 式(7)
```

ここで、`mu = self.properties["mu"]`の、`self.properties`には各要素の媒質の透磁率や誘電率、導電率などの値がnumpy配列で入っている。
高速化のため、forループではなく、numpyで一括計算している。式(7)と見比べてみよう。

`calcE`は以下の様になる。

```python
class Solver:
    def calcE(self):
        I=len(self.E)
        eps = self.properties["eps"][1:I-1]
        sgm = self.properties["sgm"][1:I-1]
        dt = self.grid.dt
        dz = self.grid.dz
        c1 = (2*eps-sgm*dt)/(2*eps+sgm*dt)
        c2 = 2*dt/((2*eps+sgm*dt)*dz)
        self.E[1:I-1] = c1*self.E[1:I-1] - c2*(self.H[1:]-self.H[:len(self.H)-1]) # 式(8)
```

これも、式(8)と見比べておこう。

そして、最後に境界条件を求めよう。

```python
class Solver:
    def calcBorder(self, preE):
        c = constants.c
        dz = self.grid.dz
        dt = self.grid.dt
        if self.borderL == BoundaryCond.ABSORB:
            c1 = (c*dt - dz)/(c*dt+dz)
            self.E[0] = preE[1] +  c1 * (self.E[1] - preE[0]) # 式(12)
        elif self.borderL == BoundaryCond.DIRICHLET:
            self.E[0] = 0

        I = len(self.E) - 1
        if self.borderR == BoundaryCond.ABSORB:
            c1 = (c*dt - dz)/(c*dt+dz)
            self.E[I] = preE[I-1] +  c1 * (self.E[I-1] - preE[I]) # 式(11)
        elif self.borderR == BoundaryCond.DIRICHLET:
            self.E[I] = 0
```

ここで、引数の`preE`は一つ前の電界だ。

このプログラムは、吸収境界条件と、ディリクレ境界条件を表している。
`self.borderL`は左境界、`self.borderR`は右境界で、Enumクラス`BoundaryCond`に、`ABSORB`, `DIRICHLET`から判断できる。

まず、`BoundaryCond.ABSORB`は吸収境界条件である。
これは、式(11), 式(12)を再現している。

そして、`BoundaryCond.DIRICHLET`はディリクレ境界条件で、右辺を0にしている。
これは、`E`の値を0にしているだけだ。

### 3.3.3 収束判定する
吸収境界条件でFDTDを回し続けていると、すべての要素で`E`,`H`の値が0になる。
このあとも、FDTDを回し続けるのは無意味なので、ある程度収束したら、打ち切るようにしよう。

このために、以下のような`judgeConvergence`関数を作った。

```python
class Solver:
    def judgeConvergence(self, cc = 1e-3, bPrint=False):
        emax = self.pulse.E
        eave = np.linalg.norm(self.E, 1)/len(self.E)
        conv = eave/emax
        bConvergence = conv < cc
        if bPrint:
            print("time="+str(self.time)+", ave/max="+str(conv))
        return bConvergence
```

簡単にいうと、最大値と平均の比がある値以下になったら打ち切るというものだ。
$$
\frac{\text{max}({E^0}[k])}{\text{average}(E^t[k])} < \epsilon
$$
最大値は、$t=0$のときの最大値、平均は今の時間の平均値だ。

収束判定条件に関しては自分の好きなように設計しても良いとおもう。

### 3.3.4. フーリエ離散変換
フーリエ離散変換の式を再掲しよう。
$$
\mathcal{F}(\omega) = \sum_{n=0}^{N} E(n\Delta t)e^{-j n\Delta t  \omega}\Delta t \tag{10}
$$

```python
class Solver:
    def calcFourier(self, eh=""):
        j = 1j
        f = np.exp(-j*self.time*2*np.pi*self.freqs)*self.grid.dt # 式(10)
        for observe in self.observes:
            if eh != "H":
                self.fourierE[observe["name"]] += self.E[observe["E"]]*f
            if eh != "E":
                self.fourierH[observe["name"]] += self.H[observe["H"]]*f
```

ここで、`self.freqs`は計算する周波数(numpy)、`self.observes`は観測点(リスト)だ。
今回は、全要素で計算せずに、観測点で設定した部分だけをフーリエ変換することにした。

各ステップごとに、どんどん足していくと勝手に4行目の式(10)を計算してることになっていることになる。
引数`eh`は、`E`に設定すれば電界のフーリエ変換を、`H`に設定すれば磁界のフーリエ変換をする。

他にも細かいことは色々あるが、ここでは省略しよう。
詳しくはソースコードを見てほしい

## 3.4 解析結果
計算結果が得られたら、次は可視化だ。
これについては、matplotlibが基本になるので解説は省略しよう。

ここでは、解析結果の例をあげる。

### 3.4.1 媒質なし、両側とも吸収境界条件

解析領域は、0~1、分割数は500で計算した。
両方の境界条件は吸収境界条件にし、媒質は置かない。

```python
N=500

grid = YeeGrid.YeeGrid()
grid.setGrid(1, N)

solver = Solver.Solver(grid)
solver.setInit(1,100.0e-12)
solver.setBorder("absorb", "absorb")
solver.calc()

graph = Graph.SpaceGraph(solver)
graph.animate(intervl=50)

graph = Graph.FourierGraph(solver)
graph.plotE(bSave=True)
```
下の`Graph`がグラフを出力するモジュールだ。

結果を見てみよう

- アニメーション

(gif動画のため、動かない場合はクリック)
![anime1(loss).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/a34b1605-54c8-d78e-a5e9-787ceff7b2b9.gif)


左側の軸が電界の値、右側の軸が磁界の値だ。

いい感じに流れているのが見て取れる。右境界で吸収されているのもわかるだろう。


- 周波数数特性(観測点: 0.5)

![fourierE1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/7d3c8323-97e2-db6c-148c-9128d0f30de0.png)


こちらもいい感じだ。
前に説明したように、ガウス関数のフーリエ変換はガウス関数だ。
ログスケールのためちょっと見にくいかもしれないが、滑らかになっているのがわかる。

### 3.4.2 媒質なし、右側境界をディリクレ境界条件

(gif動画のため、動かない場合はクリック)
![anime2(loss).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/81ecfe84-4c73-4320-ae19-5481ef1b3087.gif)


上手く反射しているのがわかる。
上のグラフより速度が早いように見えるが、Qiitaの制限である1MB以下にしようとしたらこうなってしまった。
同じ速度であることに注意しよう。

### 3.4.3 媒質あり、両側境界を吸収境界条件

比誘電率$\varepsilon_r=3$の媒質を中心に入れてみたものだ。
ガウスパルスが誘電体に入る動きが面白い

(gif動画のため、動かない場合はクリック)
![anime3(loss).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/485802/63dff883-81c7-c8f3-052d-e75d138e99a0.gif)


# 4. まとめ
今回はFDTD法の基礎について解説した。
FDTD法のアルゴリズムと、フーリエ変換、波源モデル、そしてMurの吸収境界条件を解説し、そしてこれらを用いてpythonで実装した。

まだまだ2次元以上のFDTD法の実装やPML吸収境界条件、遠方界の解析など、書ききれないことが多いが、ひとまずここまでにしておこう。予定より遥かに長くなってしまった。

筆者もまだ学んで数ヶ月のため、まだまだ勉強中だ。
なにか指摘があれば、コメントしてほしい。


# 参考文献
- [FDTD法による電磁界およびアンテナ解析](https://www.coronasha.co.jp/np/isbn/9784339006896/)
    - 式展開などをかなり細かく(ちょっと冗長に感じるほど)してくれるので、これ一冊読めば、FDTD法はだいたい理解したと思っていいと思う。
- [株式会社EEM](http://www.e-em.co.jp/)
    - [OpenFDTD](http://www.e-em.co.jp/OpenFDTD/index.html)の説明や、[EEMドキュメントの理論説明書, セミナー資料](http://www.e-em.co.jp/document.htm)はめちゃくちゃ参考にさせてもらった。
    - 株式会社と書いているが、社員がたった一人らしい。それで４つの製品と４つのOSS(フリーソフト)を出してるのはすごい。OpenFDTDはフリーのFDTDソフトとして有名だと思う。


[^1]: $\nabla \cdot$や$\nabla \times$の説明もしようと思ったが、冗長になると思ったのでやめた。
[^2]: 手元にある電磁気学の本にあまり詳しい説明がなかったのだが、おそらくこれも誘電率と同様に近似式だと思う。
[^3]: この図を作るのめっちゃ頑張った(どうでもいい)
[^4]: 中心差分は二次精度と呼ばれ、$\Delta x$が10倍小さくなると、誤差が$10^2$倍小さくなる。これは厳密にはテイラー展開で導出できる。
[^5]: これも二次精度であることをテイラー展開で導出できる。
[^6]: マクスウェル方程式で$x$,$y$方向に一様と考えてもいい。電磁気学の知識を使わず、式のみで考えたいときは、こっちで考えたほうがわかりやすい。
[^7]: なぜ式(9)がこんな形になるかは説明しないが、色んな所で解説されてるはずだ。
[^8]: このとき、$\Delta t$が省略されている流儀が多いが、今回は直感的なわかりやすさを重視して、$\Delta t$を利用しよう。全体のスケールが変わるである。
[^9]: さらに、DFTを高速にした高速フーリエ変換(FastFourierTransform;FFT)もある。簡単に概要を説明すると、計算の順番をうまいこと工夫して、DFTの計算量を減らしたものだ。
[^10]: いろいろ調べたが、パルス幅というのはいろんな定義があるようだ。[こちらのサイト](https://lr.ttu.ee/irm0120/2016kevad/width_table.pdf)を参考に、RMSwidthを採用した。