■概要

物理現象のシミュレーションは、対象の物理現象の視覚的な理解や拡張にとって重要なものです。
様々な物理現象の中で、電場／磁場中で運動する電荷をシミュレーションしてみます。
力を受けて運動をする電荷は下の図のように電場、磁場によって影響を受けます。
影響を受けた電荷は速度が変化し、その結果経路が変化します。

「※－左が電場の例、右が直線電流に対する磁場の例」

例えば、青い色の静止した電荷があり、力を受け運動する赤い色の電荷があると仮定すると、
赤い色の電荷は青い色の電荷の電場に影響を受けるため、移動方向や速さが変化します。
（逆に言えば、静止した電荷の位置や大きさなどを変えることで赤い色が動く軌道も変わります。）

従って、その特徴を利用し青い色の電荷の属性を調整しながら、
動き始めた赤色の電荷がマーク（x表示）された位置を通るようにシミュレーションしてみました。

■理論の背景

電磁場上の電荷の動きを計算するためには電場、磁場、速度を含む式を解かなくてはなりません。
それに最も適した公式であるローレンツ力（Lorentz force）公式を利用すると
電磁場上に存在している電荷に適用される力を電場と磁場に分けて表現することができます。

$$
F = q\ (E\ +\ v\ \times B)
$$

■シミュレーション実装

Google Colabで実装しffmpegを使って動画を保存しました。

・必要なパラメーター

変数名	単位	意味
charge_site	Meter(m)	静止した電荷の位置
e	Coulomb(C)	電気素量
m_P	mass of Proton	陽子の質量
B	Tesla(T)	磁束密度
x0	Meter(m)	電子の初期位置（x,y）
v0	Micro Per Sec(μm/s)	電子の初期速度
hs	Second(s)	時間間隔
xRange	Meter(m)	軌道を表すxの範囲
yRange	Meter(m)	軌道を表すyの範囲

#設定値のソースコード
#e：Quantity of Electric Charge
#m_p：mass of proton
#epsilon_0：permeability
from scipy.constants import e, m_p, epsilon_0
import numpy as np

cmr = e/m_p # cmr = charge to mass 比率 (q/m)
k = e/(4*np.pi*epsilon_0) # k = 1/(4pi*epsilon_0)
scale = 0.000001

charge_site = [[-1*scale,-2*scale],[2*scale,-0.5*scale],[0,1*scale]]
charge = [1*e,-1*e,1*e]
Bz = 2.5e-5
x0 = [-4.*scale,-4.*scale]
v0 = [0.,np.sqrt(cmr)*scale]

hs = 0.05/np.sqrt(cmr)
xRange = (-5*scale,15*scale)
yRange = (-8*scale,8*scale)

・静止した電荷から受ける力を計算するためにローレンツ力（LorentzForce）を実装しました。
ここで電場（E）を計算しますが、電場は各々の静止した電荷との距離と電荷の大きさを用いて計算出来ます。
磁束密度（B）の場合は、静止した電荷が置いてある平面に対し、垂直になるように定義しました。（磁場が計算に入ると少々複雑になるため）
変数cmrはq/mで

F=q(E+v \times B) → ma=q(E+v \times B) → a=\frac{q}{m}(E+v \times B)

と言う式から、加速度（a）を求めることができます。

#ソースコード
#動いている静止した電荷の速度と位置を基にして、
#静止した電荷の加速度を返り値として渡す
def lorentzForce(v,x): # v = 速度, x = 位置
  E = np.array([0.,0.,0.])

  #静止した電荷たちとの相対的な距離を計算し
  #入力した静止した電荷が受ける電力を求める
  for q, cs in zip(charge, charge_sites):
    E += q*(cs-x)/np.linalg.norm((cs-x))**3
  return cmr*(E+np.cross(v, B)) # cmr = charge to mass 比率 (q/m)

・加速度から位置を計算するための微分方程式の解
簡単な方法としてはオイラー法（Euler Method）が存在しますが、
誤差がどんどん大きくなる短所があるため、ルンゲクッタ法（Runge-Kutta 4th Method）を使いました。

‐Euler Methodの場合

$a_{t_i}\ =\ \frac{q}{m}\ (E_{t_i}\ +\ v_{v_i}\ \times\ B)$
$v_{t+1}\ =\ v_{t}\ +\ a_{t_i}\ dt$
$x_{t+1}\ =\ x_{t_i}\ +\ v_{t_{i+1}}\ dt$

‐Runge-Kutta 4th Methodの場合

$Suppose\ x_t\ =\ x_0,\ v_t\ =\ v_0,\ x_{t+1}\ =\ g(v_t,t)\ Then\ RK\ 4'th\ method\ is$
$f_0\ =\ f(x_0,t_0)$
$g_0\ =\ g(v_0,t_0)$
$f_1\ =\ f(x_0\ + \frac{1}{2}h_x,\ t_0\ + \frac{1}{2}h_xf_0)$
$g_1\ =\ g(v_0\ + \frac{1}{2}h_v,\ t_0\ + \frac{1}{2}h_vg_0)$
$f_2\ =\ f(x_0\ + \frac{1}{2}h_x,\ t_0\ + \frac{1}{2}h_xf_1)$
$g_2\ =\ g(v_0\ + \frac{1}{2}h_v,\ t_0\ + \frac{1}{2}h_vg_1)$
$f_3\ =\ f(x_0\ +\ h_x,\ t_0\ +\ h_xf_2)$
$g_3\ =\ g(v_0\ +\ h_v,\ t_0\ +\ h_vg_2)$
$x_1\ =\ x_0\ +\ \frac{1}{6}(f_0\ +\ 2f_1\ +\ 2f_2\ +\ f_3)h_x$
$v_1\ =\ v_0\ +\ \frac{1}{6}(g_0\ +\ 2g_1\ +\ 2g_2\ +\ g_3)h_v$

Runge-Kutta 4th Methodのイメージは↓になります。

#Runge-Kutta 4th Method
def rk2(a, b, fa, fb, hs):
  a1 = fa(a,b)*hs
  b1 = fb(a,b)*hs
  ak = a+a1*0.5
  bk = b+b1*0.5
  a2 = fa(ak, bk)*hs
  b2 = fb(ak, bk)*hs
  ak = a+a2*0.5
  bk = b+b2*0.5
  a3 = fa(ak, bk)*hs
  b3 = fb(ak, bk)*hs
  ak = a+a3
  bk = b+b3
  a4 = fa(ak, bk)*hs
  b4 = fb(ak, bk)*hs
  a = a + (a1 + 2*a2 + 2*a3 + a4)/6
  b = b + (b1 + 2*b2 + 2*b3 + b4)/6

  return (a, b)

・加速度から速度と位置を求めるために微分方程式を解くルンゲクッタ法（Runge-Kutta 4th Method）を使いました。
1回のRunge-Kutta 4th Methodの実行で1ステップに対する解が求められるため、
ステップの経過に伴う電荷の速度や位置を計算するためには複数回反復しなければなりません。
この場合は、4x400=1600回のRunge-Kutta 4th Methodを実行し、
4回ごとに画像を保存して、動画を作りました。
その他、グラフを確認するためmatplotlibを使いました。
ちなみに、位置を求めるためには速度と加速度が含まれたLorentz Force式と位置を含む式が必要ですが、
関数NextXで位置を定義しました（そのままvを戻してますが、指数関数などの関数だと思ってください）

#電荷軌道計算ソースコード
import matplotlib.animation as manimation
import gc

verticalTarget = (0,0) #ターゲット1の位置
horizontalTarget = (0.000002,0.000002) #ターゲット2の位置

FFMpegWriter = manimation.writers['ffmpeg']
metadata = dict(title='電荷の軌道', artist='Matplotlib', comment='Obenef')
writer = FFMpegWriter(fps=15, metadata=metadata)

def find_trajectory(charge_site, charge, Bz):
  import matplotlib.pyplot as plt
  global charge_sites, B
  x = x0[:]
  v = v0[:]
  x.append(0,)
  v.append(0,)
  x = np.array(x)
  v = np.array(v)

  z = np.zeros((len(charge_site),1), dtype=float)
  charge_sites = np.append(charge_site, z, axis=1)
  B = np.array([0.,0.,Bz])
  (datX, datY) = (charge_sites[:,0],charge_sites[:,1])
  traj = np.array([x])
  fig = plt.figure()
  l = plt.plot([],[],'ro',[],[],'k-',[],[],'r-',[],[],'bo',[],[],'mx',[],[],'mx')
  lst = -30

  plt.title('trajectory of charge')
  plt.xlabel('x(m)')
  plt.ylabel('y(m)')
  plt.grid('on')
  plt.xlim(xRange)
  plt.ylim(yRange)

  #動画でセーブ
  with writer.saving(fig, "電荷の軌道.mp4", 100):
    for i in range(400):
      for j in range(4):
        v, x = rk2(v, x, lorentzForce, nextX, hs) #RK2のステップ
        traj = np.append(traj, [x], axis=0)
      l[0].set_data(x[0],x[1])
      l[1].set_data(traj[:,0],traj[:,1])
      l[2].set_data(traj[lst:,0],traj[lst:,1])
      l[3].set_data(datX,datY)
      l[4].set_data(verticalTarget)
      l[5].set_data(horizontalTarget)
      writer.grab_frame()
  gc.collect()

def nextX(v,x):
  return v

#実行コード
charge_site = [[-1*scale,-2*scale],[2*scale,-0.5*scale],[0,1*scale]] #青い色の電荷の位置
charge = [1*e,-1*e,1*e] # calibrate electric charge
find_trajectory(charge_site, charge, Bz)

■シミュレーション結果

シミュレーションの結果は以下になります。

■感想

物理の理論はあらゆる場所で用いられています。
その物理現象をシミュレーションを通して理解を深め、
応用にも繋げることができれば良いなと思いました。
次回また機会があれば、さらに面白いシミュレーションを発信してみたいです。
以上です、ご覧いただきありがとうございました。

電磁場上で動く電荷のシミュレーションについて

■概要

■理論の背景

■シミュレーション実装

■シミュレーション結果

■感想