FDTDを実装する.

誰もが知ってるマクスウェル方程式

• $\bm{\nabla}\cdot \bm{D} =\rho$
• $\bm{\nabla}\cdot \bm{B} = 0$
• $\bm{\nabla}\times \bm{E} = -\frac{\partial \bm{B}}{\partial t}$
• $\bm{\nabla}\times \bm{H} = \bm{j} + \frac{\partial \bm{D}}{\partial t}$
  よく, $E,B$のみでマクスウェル方程式を書く人がいるがこちらを出発点にする.
  もちろん,
  $\rho, 0$は実数であり, $\bm{B} = \mu\bm{H}$, $\bm{D} = \epsilon \bm{E}$
  によって, $\mu, \epsilon$を定義する. ここで, $\mu, \epsilon$はテンソルである.

その他条件(補足)

• $\bm{j} = \sigma \bm{E}$
  $\sigma$はテンソルである.
  オームの法則もどきである.
  これはよく使われる式だが, $\sigma$が単純な実数だとするならば異様なことが起きる.
  ご存知のとおり, $\bm{E}$は横波である. であるならば$\sigma\bm{E}$は横波である.
  であるならば, $\bm{j}$は横波である.
  異様である.
  導線内の$\bm{j}$について考えると, 導線内では導線の接線と$\bm{j}$(の向き)が必ず一致してしまう. つまり, 普通に考えれば$\bm{j}$は縦波と考えられるのだ.

• $\bm{\nabla}\cdot (\bm{j}+\frac{\partial \bm{D}}{\partial t})=0$
  電荷保存を表すが, これはマクスウェル方程式の最後の式の発散を取ればいいので鼻くそほどの価値もない.

• $\bar{\bm{B}}=\bm{B}+\bm{B}_0$
  $\bm{B}$に定数$\bm{B}_0$を上乗せして$\bm{B}$の代わりに$\bar{\bm{B}}$を使ったって, マクスウェル方程式って変わらない.
  (ここで, ベクトルポテンシャルやゲージ変換の話はしない)

• $\bm{B}, \bm{H}$
  今更感あるけど$\bm{B}$, $\bm{H}$は別物です.
  $\bm{B}$は$\bm{f}=\bm{j}\times\bm{B}$によって定義され, $\bm{H}$はマクスウェル方程式の最後の式で定義されます.

マクスウェル方程式の差分化・離散化

$$\bm{\nabla}\times \bm{E} = -\frac{\partial \bm{B}}{\partial t}$$

を差分化する.
式番号の付け方がわからないのでこの式をファラデーの式と呼ぶことにする.
空間を格子状に分割し, 各点を$(i,j,k)$のように三つの整数と半整数の組みで表現する.
格子点間の距離は$x,y,z$軸方向それぞれ$\Delta x, \Delta y, \Delta z$とし,
計算ステップ幅を$\Delta t$とする.
各格子点に$\bm{E}$を割り当てると, 点$(i,j,k)$の$\bm{E}$は$\bm{E}(i,j,k)$として表す. また, $\bm{E}$の$x,y,z$成分を$E_x, E_y, E_z$とする.
$n$ステップ目, 点(i,j,k)の$E_x, E_y, E_z$を$E^n_x(i,j,k), E^n_y(i,j,k), E^n_z(i,j,k)$とする.

なんか式見にくくね? そんなことない?

$\bm{E}$の更新式は

$$E^n_x(i+1/2,j,k)=CEX(i+1/2,j,k)E^{n-1}_x(i+1/2,j,k)\\ +CEXLY(i+1/2,j,k)\{H^{n-1/2}_z(i+1/2,j+1/2,k)-H^{n-1/2}_z(i+1/2,j-1/2,k)\}\\ -CEZLZ(i+1/2,j,k)\{H^{n-1/2}_y(i+1/2,j,k+1/2)-H^{n-1/2}_y(i+1/2,j,k-1/2)\}$$

$$E^n_y(i,j+1/2,k)=CEY(i,j+1/2,k)E^{n-1}_y(i,j+1/2,k)\\ +CEYLZ(i,j+1/2,k)\{H^{n-1/2}_x(i,j+1/2,k+1/2)-H^{n-1/2}_x(i,j+1/2,k-1/2)\}\\ -CEYLX(i,j+1/2,k)\{H^{n-1/2}_z(i+1/2,j+1/2,k)-H^{n-1/2}_z(i-1/2,j+1/2,k)\}$$

$$E^n_z(i,j,k+1/2)=CEZ(i,j,k+1/2)E^{n-1}_z(i,j,k+1/2)\\ +CEZLX(i,j,k+1/2)\{H^{n-1/2}_y(i+1/2,j,k+1/2)-H^{n-1/2}_y(i-1/2,j,k+1/2)\}\\ -CEZLY(i,j,k+1/2)\{H^{n-1/2}_x(i,j+1/2,k+1/2)-H^{n-1/2}_x(i,j-1/2,k+1/2)\}$$

ただし,

$$CEX(i+1/2,j,k)=\frac{1-\frac{\sigma(i+1/2,j,k)\Delta t}{2\sigma(i+1/2,j,k)}}{1+\frac{\sigma(i+1/2,j,k)\Delta t}{2\sigma(i+1/2,j,k)}}$$

$$CEXLY(i+1/2,j,k)=CEXLZ(i+1/2,j,k)\Delta z$$

$$CEY(i,j+1/2,k)=\frac{1-\frac{\sigma(i,j+1/2,k)\Delta t}{2\sigma(i,j+1/2,k)}}{1+\frac{\sigma(i,j+1/2,k)\Delta t}{2\sigma(i,j+1/2,k)}}$$

$$CEYLZ(i,j+1/2,k)=CEXLZ(i,j+1/2,k)\Delta x$$

$$CEZ(i,j,k+1/2)=\frac{1-\frac{\sigma(i,j,k+1/2)\Delta t}{2\sigma(i,j,k+1/2)}}{1+\frac{\sigma(i,j,k+1/2)\Delta t}{2\sigma(i,j,k+1/2)}}$$

$$CEZLX(i,j,k+1/2)=CEXLZ(i,j,k+1/2)\Delta y$$

また$\bm{H}$の更新式は

$$H^{n+1/2}_x(i,j+1/2,k+1/2) = H^{n-1/2}_x(i,j+1/2,k+1/2)\\ - CHXLY(i,j+1/2,k+1/2)\{E^{n}_z(i,j+1,k+1/2) - E^{n}_z(i,j,k+1/2)\} \\+ CHXLZ(i,j+1/2,k+1/2)\{E^{n}_y(i,j+1/2,k+1) - E^{n}_y(i,j+1/2,k)\}$$

$$H^{n+1/2}_y(i+1/2,j,k+1/2) = H^{n-1/2}_y(i+1/2,j,k+1/2)\\ - CHYLZ(i+1/2,j,k+1/2)\{E^{n}_x(i+1/2,j,k+1) - E^{n}_x(i+1/2,j,k)\} \\+ CHYLX(i+1/2,j,k+1/2)\{E^{n}_z(i+1,j,k+1/2) - E^{n}_z(i,j,k+1/2)\}$$

$$H^{n+1/2}_z(i+1/2,j+1/2,k) = H^{n-1/2}_z(i+1/2,j+1/2,k)\\ - CHZLX(i+1/2,j+1/2,k)\{E^{n}_y(i+1,j+1/2,k) - E^{n}_y(i,j+1/2,k)\} \\+ CHZLY(i+1/2,j+1/2,k)\{E^{n}_x(i+1/2,j+1,k) - E^{n}_x(i+1/2,j,k)\}$$

ただし,

$$CHXLY(i,j+1/2,k+1/2)=\frac{\Delta t}{\mu(i,j+1/2,k+1/2)}\frac{1}{\Delta y}$$

$$CHXLZ(i,j+1/2,k+1/2)=\frac{\Delta t}{\mu(i,j+1/2,k+1/2)}\frac{1}{\Delta z}$$

$$CHYLZ(i+1/2,j,k+1/2)=\frac{\Delta t}{\mu(i+1/2,j,k+1/2)}\frac{1}{\Delta z}$$

$$CHYLX(i+1/2,j,k+1/2)=\frac{\Delta t}{\mu(i+1/2,j,k+1/2)}\frac{1}{\Delta x}$$

$$CHZLX(i+1/2,j+1/2,k)=\frac{\Delta t}{\mu(i+1/2,j+1/2,k)}\frac{1}{\Delta x}$$

$$CHZLY(i+1/2,j+1/2,k)=\frac{\Delta t}{\mu(i+1/2,j+1/2,k)}\frac{1}{\Delta y}$$

こんなん絶対コーディングしている最中に間違えるやん.

• $E_x(i+1/2,j,k)$ は Ex(i,j,k)

• $E_y(i,j+1/2,k)$ は Ey(i,j,k)

• $E_z(i,j,k+1/2)$ は Ez(i,j,k)

• $H_x(i,j+1/2,k+1/2)$ は Hx(i,j,k)

• $H_y(i+1/2,j,k+1/2)$ は Hy(i,j,k)

• $H_z(i+1/2,j+1/2,k)$ は Hz(i,j,k)

• $CEX(i+1/2,j,k)$ は CEX(i,j,k)

• $CEXLY(i+1/2,j,k)$ は CEXLY(i,j,k)

• $CEY(i,j+1/2,k)$ は CEY(i,j,k)

• $CEYLZ(i,j+1/2,k)$ は CEYLZ(i,j,k)

• $CEZ(i,j,k+1/2)$ は CEZ(i,j,k)

• $CEZLX(i,j,k+1/2)$ は CEZLX(i,j,k)

• $CHXLY(i,j+1/2,k+1/2)$ は CHXLY(i,j,k)

• $CHXLZ(i,j+1/2,k+1/2)$ は CHXLZ(i,j,k)

• $CHYLZ(i+1/2,j,k+1/2)$ は CHYLZ(i,j,k)

• $CHYLX(i+1/2,j,k+1/2)$ は CHYLX(i,j,k)

• $CHZLX(i+1/2,j+1/2,k)$ は CHZLX(i,j,k)

• $CHZLY(i+1/2,j+1/2,k)$ は CHZLY(i,j,k)

とする.

ただし, CEXなどは媒質に依存する.
媒質の種類は有限なので, 格子点の数だけCEXなどを用意するより媒質の種類を元にCEXなどを決定した方がメモリが省ける.
つまり, 各媒質にナンバリングし, idというN X N X Nの配列に媒質のナンバーを格納すれば, 個別にCEX, CEY, CEZなどを用意する必要がなくなる.
CEX(id)などとし, CEXという関数の中でif id==AIR　とすれば
return (AIRでのCEX) else if id==COPPER return (COPPERでのCEX)
のような関数を用意する.
ただし, idはE用のCEX, CEY, CEZとH用のCHLY, CHLZ, CHLX...に対してそれぞれ用意する必要がある.
つまり, ide(i,j,k)とidh(i,j,k)を用意する必要がある.
結果

• $CEX(i+1/2,j,k)$ は CEX(i,j,k) : idex(i,j,k)

• $CEXLY(i+1/2,j,k)$ は CEXLY(i,j,k) : idex(i,j,k)

• $CEY(i,j+1/2,k)$ は CEY(i,j,k) : idey(i,j,k)

• $CEYLZ(i,j+1/2,k)$ は CEYLZ(i,j,k) : idey(i,j,k)

• $CEZ(i,j,k+1/2)$ は CEZ(i,j,k) : idez(i,j,k)

• $CEZLX(i,j,k+1/2)$ は CEZLX(i,j,k) : idez(i,j,k)

• $CHXLY(i,j+1/2,k+1/2)$ は CHXLY(i,j,k) : idhx(i,j,k)

• $CHXLZ(i,j+1/2,k+1/2)$ は CHXLZ(i,j,k) : idhx(i,j,k)

• $CHYLZ(i+1/2,j,k+1/2)$ は CHYLZ(i,j,k) : idhy(i,j,k)

• $CHYLX(i+1/2,j,k+1/2)$ は CHYLX(i,j,k) : idhy(i,j,k)

• $CHZLX(i+1/2,j+1/2,k)$ は CHZLX(i,j,k) : idhz(i,j,k)

• $CHZLY(i+1/2,j+1/2,k)$ は CHZLY(i,j,k) : idhz(i,j,k)

のようにidを割り当てる.