gnuplot で円電流による磁場を数値的に解いて描く

参考：

gnuplot だけでプログラミングしたり，数値計算したりする際の参考として

円電流がつくる xz 平面内の磁場

\begin{eqnarray}
B_x(x,0,z) &=& \frac{I}{4\pi \varepsilon_0 c^2}
\int_0^{2\pi} \frac{a z \cos \phi’ }{\left\{ x^2 + z^2 + a^2 – 2 a x \cos\phi’ \right\}^{\frac{3}{2}}}\,d\phi’ \\ \ \\
B_y(x,0,z) &=& 0 \\ \ \\
B_z(x,0,z) &=& \frac{I }{4\pi \varepsilon_0 c^2}
\int_0^{2\pi} \frac{a^2 – a x \cos\phi’}{\left\{x^2 + z^2 + a^2 – 2 a x \cos\phi’ \right\}^{\frac{3}{2}}}\, d\phi’
\end{eqnarray}

適宜，定数部分を規格化して

\begin{eqnarray}
B_x(x,z) &\Rightarrow&
\int_0^{2\pi} \frac{a z \cos \phi’ }{\left\{ x^2 + z^2 + a^2 – 2 a x \cos\phi’ \right\}^{\frac{3}{2}}}\,d\phi’ \\ \ \\
B_y(x,z) &=& 0 \\ \ \\
B_z(x,z) &\Rightarrow&
\int_0^{2\pi} \frac{a^2 – a x \cos\phi’}{\left\{x^2 + z^2 + a^2 – 2 a x \cos\phi’ \right\}^{\frac{3}{2}}}\, d\phi’
\end{eqnarray}

磁場ベクトルが磁力線の接ベクトルであることから，定積分を含む微分方程式

\begin{eqnarray}
B_x(x,z) &=& \int_0^{2\pi} f_1(x,z,\phi)\,d\phi \\
B_z(x,z) &=&\int_0^{2\pi} f_2(x,z,\phi)\,d\phi \\
f_1(x,z,\phi) &\equiv& \frac{a z \cos \phi}{\left\{ x^2 + z^2 + a^2 – 2 a x \cos\phi\right\}^{\frac{3}{2}}} \\
f_2(x,z,\phi) &\equiv& \frac{a^2 – a x \cos\phi}{\left\{x^2 + z^2 + a^2 – 2 a x \cos\phi \right\}^{\frac{3}{2}}}\ \\
B(x, z) &\equiv& \sqrt{B_x(x,z)^2 + B_z(x,z)^2} \\
F_1(x, z) &\equiv& \frac{B_x(x,z)}{B(x,z)} \\
F_2(x, z) &\equiv& \frac{B_z(x,z)}{B(x,z)} \\
\frac{dx}{ds} &=& F_1(x, z)\\
\frac{dz}{ds} &=& F_2(x, z)
\end{eqnarray}

を数値的に解いて磁力線を求める。

円電流による磁場：plot 版

In [1]:
%gnuplot inline svg size 480,480 fixed enhanced font 'Noto Sans CJK JP,14'
In [2]:
# 被積分関数
f1(x,z,ph) = a*z*cos(ph)/sqrt(x**2+z**2+a**2-2*a*x*cos(ph))**3
f2(x,z,ph) = (a**2-a*x*cos(ph))/sqrt(x**2+z**2+a**2-2*a*x*cos(ph))**3

# 2*pi 積分区間。
# N  分割数 N は偶数とすること（シンプソン法の決め事）
# 三項演算子による再帰的定義関数で和をとる。

simp_f1(x, z, N, i) = (i == 2 ? \
  (2*pi)/(N*3.0)* \
    (f1(x,z,(i-2)*2*pi/N) + 4.0*f1(x,z,(i-1)*2*pi/N) + f1(x,z,i*2*pi/N)) : \
  simp_f1(x, z, N, i-2) + \
  (2*pi)/(N*3.0)* \
    (f1(x,z,(i-2)*2*pi/N) + 4.0*f1(x,z,(i-1)*2*pi/N) + f1(x,z,i*2*pi/N)))

simp_f2(x, z, N, i) = (i == 2 ? \
  (2*pi)/(N*3.0)* \
    (f2(x,z,(i-2)*2*pi/N) + 4.0*f2(x,z,(i-1)*2*pi/N) + f2(x,z,i*2*pi/N)) : \
  simp_f2(x, z, N, i-2) + \
  (2*pi)/(N*3.0)* \
    (f2(x,z,(i-2)*2*pi/N) + 4.0*f2(x,z,(i-1)*2*pi/N) + f2(x,z,i*2*pi/N)))

# 円電流の半径
a = 1.

# 実際の計算の際は，simp_f(a, b, N, N) と同じ N を2回書く
# N = 40 くらいでよさそう

Bx(x,z) = simp_f1(x, z, 40, 40)
Bz(x,z) = simp_f2(x, z, 40, 40)

B(x,z) = sqrt(Bx(x,z)**2 + Bz(x,z)**2)

F1(x,z) = Bx(x,z)/B(x,z)
F2(x,z) = Bz(x,z)/B(x,z)

# 計算範囲
xmin = 0
xmax = 3
ymin = 0
ymax = 3

# ルンゲ・クッタ法のパラメータ
h = 0.05
Nend = 100
array X[Nend]
array Y[Nend]

# データファイルの削除
# あとで append で開くので
! rm -f en.txt
! rm -f en-line.txt

# 初期条件
# x 軸上からはじめる
do for [j=1:4]{
  x0 = 0.2*j
  y0 = 0

  # 4次のルンゲ・クッタ法
  do for [i=1: Nend]{
    X[i] = x0
    Y[i] = y0
    k1 = h*F1(x0, y0)
    m1 = h*F2(x0, y0)
    k2 = h*F1(x0+0.5*k1, y0+0.5*m1)
    m2 = h*F2(x0+0.5*k1, y0+0.5*m1)
    k3 = h*F1(x0+0.5*k2, y0+0.5*m2)
    m3 = h*F2(x0+0.5*k2, y0+0.5*m2)
    k4 = h*F1(x0+k3, y0+m3)
    m4 = h*F2(x0+k3, y0+m3)
    X1 = x0 + (k1 + 2.*k2 + 2.*k3 + k4)/6
    Y1 = y0 + (m1 + 2.*m2 + 2.*m3 + m4)/6
    iend = i
    # 計算範囲を逸脱したら do loop から break
    if (X1 < xmin) {break}
    if (X1 > xmax) {break}
    if (Y1 < ymin) {break}
    if (Y1 > ymax) {break}
    x0 = X1
    y0 = Y1
  }
  # 磁力線の座標データファイル
  set print "en-line.txt" append
  # 第1象限
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  if (Y[iend] > 0.5) {print "         "}
  # 第4象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
  print "         "
  # 第2象限
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  if (Y[iend] > 0.5) {print "         "}
  # 第3象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  print "         "
  set print

  # ベクトル始点の座標データファイル
  # 適宜間引く
  set print "en.txt" append
  # 第1象限
    print sprintf("%8.4f %8.4f", X[15], Y[15])
  # 第2象限
    print sprintf("%8.4f %8.4f", -X[15], Y[15])
  # 第3象限
    print sprintf("%8.4f %8.4f", -X[20], -Y[20])
  # 第4象限
    print sprintf("%8.4f %8.4f", X[20], -Y[20])
  set print
}

# 磁力線の座標データファイル
# z 軸上のデータ追加
set print "en-line.txt" append
  do for [i=0:100] {
    z = -3 + 6./100*i
    print sprintf("%8.4f %8.4f", 0, z)
  }
set print

set print "en.txt" append
  do for [i=0:10:5] {
    z = -2 + 4./10*i
    print sprintf("%8.4f %8.4f", 0, z)
  }
set print

reset
unset xtics
unset ytics
unset border

set xrange [-3:3]
set yrange [-3:3]
set size ratio 1
set zeroaxis

set title "円電流による磁場（xz面）"
unset colorbox
# green.pal 改
set palette defined (\
  0 '#E5F5E0',\
  1 '#C7E9C0',\
  2 '#A1D99B',\
  3 '#74C476',\
  4 '#41AB5D',\
  5 '#238B45',\
  6 '#005A32',\
  7 '#002000' \
)
set key inside sample 1
# 電線の位置
set label 1 point pt 6 ps 1 lc "red" at a, 0
set label 2 center at first a, 0 "×" tc "red" font ",14"
set label 3 point pt 6 ps 1 lc "red" at -a, 0
set label 4 center at first -a, 0 "・" tc "red" font ",14"
scaling=0.05
plot "en-line.txt" \
  u 1:2:(log(B($1,$2))) w l lw 3 lc palette notitle, \
  "en.txt"   \
  u 1:2:(scaling*F1($1,$2)):(scaling*F2($1,$2)):(log(B($1,$2))) \
  w vec lc palette lw 8 filled head title "磁場"

In [3]:
set output "./en-jiba2.svg"
replot
set output

円電流回路による磁場：splot 版

In [4]:
%gnuplot inline svg size 640,640 fixed enhanced font 'Noto Sans CJK JP,14'
In [5]:
# 半円の座標データファイル
# theta の分割数
Nth = 100
# 半径
a = 1

set print "han-en.txt"
  do for [i=0:Nth]{
    th = pi/Nth * i
    print sprintf("%8.4f %8.4f", a*cos(th), a*sin(th))
  }
  print "      "
set print

# 半円柱の電流ベクトルデータファイル
set print "han-en-vec2.txt"
#  print 0, 1, 0, -1, 0, 0
  print 0, -1, 0, 1,  0, 0
set print

reset
unset xtics
unset ytics
unset ztics
unset border

set xrange [-3:3]
set yrange [-3:3]
set zrange [-3:3]
set view 60,15,,

set key inside sample 1
set title "円電流による磁場"
unset colorbox
# green.pal 改
set palette defined (\
  0 '#E5F5E0',\
  1 '#C7E9C0',\
  2 '#A1D99B',\
  3 '#74C476',\
  4 '#41AB5D',\
  5 '#238B45',\
  6 '#005A32',\
  7 '#002000' \
)
scaling=0.05
set parametric
splot [u=-2:2] u,0,0 dt(5,10) lw 0.5 lc 'black' notitle, \
               0,u,0 dt(5,10) lw 0.5 lc 'black' notitle, \
  "han-en.txt" \
    u 1:2:(0):(-$2) w l lw 4 lc 'red' notitle, \
  "en-line.txt" \
    u 1:(0):2:(log(B($1,$2))) w l lw 3 lc palette notitle, \
  "en.txt" \
    u 1:(0):2:(scaling*F1($1,$2)):(0):(scaling*F2($1,$2)):(log(B($1,$2))) \
    w vec filled head lc palette lw 6 title "磁場", \
  "han-en-vec2.txt" \
    u 1:2:3:(scaling*$4):5:6 w vec filled head lw 8  lc 'red' title "電流"  , \
  "han-en.txt" \
    u 1:(-$2):(0):(-$2) w l lw 4 lc 'red' notitle

Out[5]:
	dummy variable is t for curves, u/v for surfaces
In [6]:
set output "./en-jiba3.svg"
replot
set output