gnuplot で磁力線を数値的に解いて描く

gnuplot だけで，磁力線を数値的に解き，磁力線に沿った磁場ベクトルを描く。磁場ベクトルが磁力線の接ベクトルであることから，微分方程式をルンゲ・クッタ法で解き，数値解をファイルに保存し，それを読み込んで磁場ベクトルの始点の座標としてベクトル場を描く。これを全て gnuplot だけでやる。

参考：

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

平行2直線電流がつくる磁場

$$\boldsymbol{B} = \frac{1}{2\pi\varepsilon_0 c^2} \frac{\boldsymbol{I}_1\times \left(\boldsymbol{\rho} – \boldsymbol{r}_1\right)}{|\boldsymbol{\rho} – \boldsymbol{r}_1|^2}
+ \frac{1}{2\pi\varepsilon_0 c^2} \frac{\boldsymbol{I}_2\times \left(\boldsymbol{\rho} – \boldsymbol{r}_2\right)}{|\boldsymbol{\rho} – \boldsymbol{r}_2|^2}$$\begin{eqnarray}
\boldsymbol{I}_1 &=& (0, 0, \pm I)\\
\boldsymbol{I}_2 &=& (0, 0, I) \\
\boldsymbol{\rho} &=& (x, y, 0) \\
\boldsymbol{r}_1 &=& (-1, 0, 0) \\
\boldsymbol{r}_2 &=& (1, 0, 0) \\
\end{eqnarray}

とし，適宜定数部分を規格化して
\begin{eqnarray}
B_x(x,y) &=& \pm \frac{-y}{\left(x-x_1\right)^2 + \left(y-y_1\right)^2 }
+\frac{-y}{\left(x-x_2\right)^2 + \left(y-y_2\right)^2} \\
B_y(x,y) &=& \pm \frac{x}{\left(x-x_1\right)^2 + \left(y-y_1\right)^2 }
+\frac{x}{\left(x-x_2\right)^2 + \left(y-y_2\right)^2} \\
\end{eqnarray}

規格化した $\hat{\boldsymbol{B}}$ は

\begin{eqnarray}
B(x,y) &=& \sqrt{\left\{B_x(x,y) \right\}^2 + \left\{B_y(x,y) \right\}^2}\\
\hat{B}_x(x,y) &=& \frac{B_x(x,y)}{B(x,y)} \\
\hat{B}_y(x,y) &=& \frac{B_y(x,y)}{B(x,y)} \\
\end{eqnarray}

平行2直線電流による磁力線

磁場ベクトル $\boldsymbol{B}$ が磁力線の接ベクトルであるから，微分方程式

$$\frac{d\boldsymbol{r}}{ds} = \hat{\boldsymbol{B}}$$

を（ルンゲ・クッタ法で）数値的に解いて，磁力線 $\boldsymbol{r}(s)$ を求める。

\begin{eqnarray}
\frac{dx}{ds} &=& \hat{B}_x(x, y) \equiv F_1(x,y) \\
\frac{dy}{ds} &=& \hat{B}_y(x, y) \equiv F_2(x,y)
\end{eqnarray}

同じ向きに電流 $I$ が流れている場合：plot 版

In [1]:
%gnuplot inline svg size 480,480 fixed enhanced font 'Noto Sans CJK JP,14'
In [2]:
# 磁場 B
Bx(x, y) = - y/(x**2 + y**2)
By(x, y) =   x/(x**2 + y**2)
# 平行2直線電流がつくる磁場
# 電線の位置
x1 = -1
y1 = 0
x2 = 1
y2 = 0

B2x(x, y) = Bx(x-x1, y-y1) + Bx(x-x2, y-y2)
B2y(x, y) = By(x-x1, y-y1) + By(x-x2, y-y2)

# 磁場 B の大きさ
B(x, y) = sqrt(B2x(x, y)**2 + B2y(x, y)**2)
# 規格化された hat B
F1(x, y) = B2x(x, y)/B(x, y)
F2(x, y) = B2y(x, y)/B(x, y)

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

# ルンゲ・クッタ法のパラメータ
h = 0.02
N = 300
array X[N]
array Y[N]

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

# 初期条件
# x 軸上からはじめる
do for [j=2:11]{
  x0 = 4./j + 1
  y0 = 0

  # 4次のルンゲ・クッタ法
  do for [i=1: N]{
    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 "para-line.txt" append
  # 第1象限
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  print "         "
  # 第2象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  print "         "
  # 第3象限
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  print "         "
  # 第4象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
  print "         "
  set print  
  
  if (j%2  ==0) {
  # ベクトル始点の座標データファイル
  # 適宜間引く
  set print "para.txt" append
  # 第1象限
  do for [i=30:30]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  # 第2象限
  do for [i=iend-20:iend-20]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  # 第3象限
  do for [i=30:30]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  # 第4象限
  do for [i=iend-20:iend-20]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
  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 "同じ向きに流れる平行2直線電流による磁場"
unset colorbox
set palette cubehelix start 2 cycles 0 saturation 3 negative
set key inside sample 2

# 矢印の長さ調整用
scaling = 0.1
# 電線の位置
set label 1 point pt 6 ps 2 lc "red" at x1, y1
set label 2 center at first x1,y1 "・" tc "red" font ",24"
set label 3 point pt 6 ps 2 lc "red" at x2, y2
set label 4 center at first x2,y2 "・" tc "red" font ",24"

plot "para-line.txt" \
  u 1:2:(B($1,$2)) w l lw 3 lc palette notitle, \
  "para.txt"   \
  u 1:2:(scaling*F1($1,$2)):(scaling*F2($1,$2)):(B($1,$2)) \
  w vec lc palette lw 6 filled head title "磁場"

In [3]:
set output "./vec-para1.svg"
replot
set output

同じ向きに電流 $I$ が流れている場合：splot 版

In [4]:
%gnuplot inline svg size 640,480 fixed enhanced font 'Noto Sans CJK JP,14'
In [5]:
# 磁場 B
Bx(x, y) = - y/(x**2 + y**2)
By(x, y) =   x/(x**2 + y**2)
# 平行2直線電流がつくる磁場
# 電線の位置
x1 = -1
y1 = 0
x2 = 1
y2 = 0

B2x(x, y) = Bx(x-x1, y-y1) + Bx(x-x2, y-y2)
B2y(x, y) = By(x-x1, y-y1) + By(x-x2, y-y2)

# 磁場 B の大きさ
B(x, y) = sqrt(B2x(x, y)**2 + B2y(x, y)**2)
# 規格化された hat B
F1(x, y) = B2x(x, y)/B(x, y)
F2(x, y) = B2y(x, y)/B(x, y)

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

# ルンゲ・クッタ法のパラメータ
h = 0.02
N = 300
array X[N]
array Y[N]

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

do for [k=1:3:2] {
  z = k
# 初期条件
# x 軸上からはじめる
do for [j=2:10:2]{
  x0 = 4./j + 1
  y0 = 0

  # 4次のルンゲ・クッタ法
  do for [i=1: N]{
    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 "para-line.txt" append
  # 第1象限
  do for [i=1:iend]{
    if (abs(X[i]-x2) < 0.07) {print "       "}
    else {
      print sprintf("%8.4f %8.4f %8.4f", X[i], Y[i], z)
      }
  }
  print "         "
  print "         "
  # 第2象限
  do for [i=iend:1:-1]{
    if (abs(X[i]-x2) < 0.07) {print "       "}
    else {
    print sprintf("%8.4f %8.4f %8.4f", -X[i], Y[i], z)}
  }
  print "         "
  print "         "
  # 第3象限
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f %8.4f", -X[i], -Y[i], z)
  }
  print "         "
  print "         "
  # 第4象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f %8.4f", X[i], -Y[i], z)
  }
  print "         "
  print "         "
  set print
  
  # ベクトル始点の座標データファイル
  # 適宜間引く
  if ((j==2)||(j==6)||(j==10)) {
  set print "para.txt" append
  # 第1象限
  do for [i=12:12]{
    print sprintf("%8.4f %8.4f %8.4f", X[i], Y[i], z)
  }
  # 第2象限
  do for [i=iend-20:iend-20]{
    print sprintf("%8.4f %8.4f %8.4f", -X[i], Y[i], z)
  }
  # 第3象限
  do for [i=30:30]{
    print sprintf("%8.4f %8.4f %8.4f", -X[i], -Y[i], z)
  }
  # 第4象限
  do for [i=iend-20:iend-20]{
    print sprintf("%8.4f %8.4f %8.4f", X[i], -Y[i], z)
  }
  set print
  }
}
}

reset
unset xtics
unset ytics
unset ztics
unset border

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

set title "同じ向きに流れる平行2直線電流による磁場"
unset colorbox
set palette cubehelix start 2 cycles 0 saturation 3 negative
set key inside sample 2
# 矢印の長さ調整用
scaling = 0.1
# 電線
set arrow 1 from x1,y1,0 to x1,y1,4.3 lw 4 lc "red" lt 3 filled head
set arrow 2 from x2,y2,0 to x2,y2,4.3 lw 4 lc "red" lt 3 filled head

set view 45,0,,

splot "para-line.txt" \
  u 1:2:3:(B($1,$2)) w l lw 3 lc palette notitle,\
  "para.txt" \
  u 1:2:3:(scaling*F1($1,$2)):(scaling*F2($1,$2)):(0):(B($1,$2)) \
  w vec lw 6 lc palette filled head title "磁場"

In [6]:
set output "./vec-para2.svg"
replot
set output

反対向きに電流 $I$ が流れている場合：plot 版

In [7]:
%gnuplot inline svg size 480,480 fixed enhanced font 'Noto Sans CJK JP,14'
In [8]:
# 磁場 B
Bx(x, y) = - y/(x**2 + y**2)
By(x, y) =   x/(x**2 + y**2)
# 平行2直線電流がつくる磁場
# 電線の位置
x1 = -1
y1 = 0
x2 = 1
y2 = 0

B2x(x, y) = -Bx(x-x1, y-y1) + Bx(x-x2, y-y2)
B2y(x, y) = -By(x-x1, y-y1) + By(x-x2, y-y2)

# 磁場 B の大きさ
B(x, y) = sqrt(B2x(x, y)**2 + B2y(x, y)**2)
# 規格化された hat B
F1(x, y) = B2x(x, y)/B(x, y)
F2(x, y) = B2y(x, y)/B(x, y)

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

# ルンゲ・クッタ法のパラメータ
h = 0.02
N = 300
array X[N]
array Y[N]

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

# 初期条件
# x 軸上からはじめる
do for [j=2:8]{
  x0 = 4./j + 1
  y0 = 0

  # 4次のルンゲ・クッタ法
  do for [i=1: N]{
    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 "para-line.txt" append
  # 第1象限
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  # 第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])
  }
  # 第3象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  print "        "
  set print
  
  if (j%2  ==1) {
  # ベクトル始点の座標データファイル
  # 適宜間引く
  set print "para.txt" append
  # 第4象限
  do for [i=25:25]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
  # 第3象限
  do for [i=25:25]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  set print
  }
}

# 内側の磁力線を追加で数値計算
h = 0.2
N = 30
array X[N]
array Y[N]
# 初期条件
# x 軸上からはじめる
do for [j=0:2]{
  x0 = 0.1*j
  y0 = 0

  # 4次のルンゲ・クッタ法
  do for [i=1: N]{
    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 (abs(Y1) < ymin) {break}
    if (abs(Y1) > ymax) {break}
    x0 = X1
    y0 = Y1
  }

  # 磁力線の座標データファイル
  set print "para-line.txt" append
  # 第4象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  # 第1象限
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f", X[i], -Y[i])
  }
 print "     "
  # 第3象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  # 第2象限
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f", -X[i], -Y[i])
  }
  print "     "
  set print 
  
  if (j%2  ==0) {
  # ベクトル始点の座標データファイル
  # 適宜間引く
  set print "para.txt" append
  # 1/4
  do for [i=3:3]{
    print sprintf("%8.4f %8.4f", X[i], Y[i])
  }
  # 3/4
  do for [i=3:3]{
    print sprintf("%8.4f %8.4f", -X[i], Y[i])
  }
  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 "反対向きに流れる平行2直線電流による磁場"
unset colorbox
set palette cubehelix start 2 cycles 0 saturation 3 negative

# 矢印の長さ調整用
scaling = 0.1
# 電線の位置
set label 1 point pt 6 ps 2 lc "red" at x1, y1
set label 2 center at first x1,y1 "×" tc "red" font ",24"
set label 3 point pt 6 ps 2 lc "red" at x2, y2
set label 4 center at first x2,y2 "・" tc "red" font ",24"

plot "para-line.txt" \
  u 1:2:(B($1,$2)) w l lw 3 lc palette notitle, \
  "para.txt"   \
  u 1:2:(scaling*F1($1,$2)):(scaling*F2($1,$2)):(B($1,$2)) \
  w vec lc palette lw 6 filled head title "磁場"

In [9]:
set output "./vec-para3.svg"
replot
set output

反対向きに電流 $I$ が流れている場合：splot 版

In [10]:
%gnuplot inline svg size 640,480 fixed enhanced font 'Noto Sans CJK JP,14'
In [11]:
# 磁場 B
Bx(x, y) = - y/(x**2 + y**2)
By(x, y) =   x/(x**2 + y**2)
# 平行2直線電流がつくる磁場
# 電線の位置
x1 = -1
y1 = 0
x2 = 1
y2 = 0

B2x(x, y) = -Bx(x-x1, y-y1) + Bx(x-x2, y-y2)
B2y(x, y) = -By(x-x1, y-y1) + By(x-x2, y-y2)

# 磁場 B の大きさ
B(x, y) = sqrt(B2x(x, y)**2 + B2y(x, y)**2)
# 規格化された hat B
F1(x, y) = B2x(x, y)/B(x, y)
F2(x, y) = B2y(x, y)/B(x, y)

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

# ルンゲ・クッタ法のパラメータ
h = 0.02
N = 300
array X[N]
array Y[N]

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

do for [k=1:3:2] {
  z = k
# 初期条件
# x 軸上からはじめる
do for [j=2:10:2]{
  x0 = 4./j + 1
  y0 = 0

  # 4次のルンゲ・クッタ法
  do for [i=1: N]{
    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 "para-line.txt" append
  # 第1象限
  do for [i=1:iend]{
    if (abs(X[i]-x2) < 0.07) {print "       "}
    else {
      print sprintf("%8.4f %8.4f %8.4f", X[i], Y[i], z)}
  }
  # 第4象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f %8.4f", X[i], -Y[i], z)
  }
  print "        "
  print "        "
  # 第2象限
  do for [i=1:iend]{
    if (abs(X[i]-x2) < 0.07) {print "       "}
    else {
    print sprintf("%8.4f %8.4f %8.4f", -X[i], Y[i], z)}
  }
  # 第3象限
  do for [i=iend:1:-1]{
    print sprintf("%8.4f %8.4f %8.4f", -X[i], -Y[i], z)
  }
  print "        "
  print "        "
  set print
  
  # ベクトル始点の座標データファイル
  # 適宜間引く
  if ((j==2)||(j==6)) {
  set print "para.txt" append
  # 第3象限
  do for [i=iend-55:iend-55]{
    print sprintf("%8.4f %8.4f %8.4f", X[i], -Y[i], z)
  }
  # 第4象限
  do for [i=iend-55:iend-55]{
    print sprintf("%8.4f %8.4f %8.4f", -X[i], -Y[i], z)
  }
  set print
  }
}
}

# 内側の磁力線を追加で数値計算
h = 0.02
N = 300
array X[N]
array Y[N]

do for [k=1:3:2] {
  z = k
# 初期条件
# x 軸上からはじめる
do for [j=0:2]{
  x0 = 0.1*j
  y0 = 0

  # 4次のルンゲ・クッタ法
  do for [i=1: N]{
    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 "para-line.txt" append
  # 1/4
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f %8.4f", X[i], Y[i], z)
  }
  print "        "
  print "        "
  # 2/4
  do for [i=iend:1:-1]{
    if (abs(X[i]-x2) < 0.1) {print "        "}
    else {
    print sprintf("%8.4f %8.4f %8.4f", X[i], -Y[i], z)
    }
  }
  print "        "
  print "        "
  # 3/4
  do for [i=1:iend]{
    print sprintf("%8.4f %8.4f %8.4f", -X[i], Y[i], z)
  }
  print "        "
  print "        "
  # 4/4
  do for [i=iend:1:-1]{
     if (abs(X[i]-x2) < 0.1) {print "        "}
    else {
   print sprintf("%8.4f %8.4f %8.4f", -X[i], -Y[i], z)
   }
  }
  print "        "
  print "        "
  set print  
  
  if ((j%2  ==0)) {
  # ベクトル始点の座標データファイル
  # 適宜間引く
  set print "para.txt" append
  # 第4象限
  do for [i=16:16]{
    print sprintf("%8.4f %8.4f %8.4f", X[i], Y[i], z)
  }
  # 第3象限
  do for [i=16:16]{
    print sprintf("%8.4f %8.4f %8.4f", -X[i], Y[i], z)
  }
  set print
  }
}
}

reset
unset xtics
unset ytics
unset ztics
unset border

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

set title "反対向きに流れる平行2直線電流による磁場"
unset colorbox
set palette cubehelix start 2 cycles 0 saturation 3 negative

# 矢印の長さ調整用
scaling = 0.2
# 電線
set arrow 1 from x1,y1,4 to x1,y1,0 lw 4 lc "red" lt 3 filled head
set arrow 2 from x2,y2,0 to x2,y2,4 lw 4 lc "red" lt 3 filled head
set key inside sample 2
set view 45,0,,

splot "para-line.txt" \
  u 1:2:3:(B($1,$2)) w l lw 3 palette notitle,\
  "para.txt" \
  u 1:2:3:(scaling*F1($1,$2)):(scaling*F2($1,$2)):(0):(B($1,$2)) \
  w vec lw 6 lc palette filled head title "磁場"

In [12]:
set output "./vec-para4.svg"
replot
set output