gnuplot で電場ベクトルの方向場を描く

正負の点電荷がつくる電場の向きを表す方向場を gnuplot で描く。以下のページの抜粋。

正負の点電荷がつくる電場

$\boldsymbol{r}_1$ にある正電荷 $q\ (>0)$ と,$\boldsymbol{r}_2$ にある負電荷 $-q$ がつくる電場ベクトル $\boldsymbol{E}$ は

\begin{eqnarray}
\boldsymbol{E}
&=& \frac{q}{4\pi \varepsilon_0} \frac{\boldsymbol{r} – \boldsymbol{r}_1}{|\boldsymbol{r} – \boldsymbol{r}_1|^3}
– \frac{q}{4\pi \varepsilon_0} \frac{\boldsymbol{r} – \boldsymbol{r}_2}{|\boldsymbol{r} – \boldsymbol{r}_2|^3}
\end{eqnarray}

$z=0$ の $xy$ 平面で2次元ベクトル場を描く。適宜無次元化を行い,
$\displaystyle \frac{q}{4\pi \varepsilon_0}\rightarrow 1$ とする。また,

$$\boldsymbol{r} = (x, y, 0), \quad \boldsymbol{r}_1 = (0, 1, 0),
\quad \boldsymbol{r}_2 = (0, -1, 0)$$

とすると,

$$E_x = \frac{x}{\left(x^2 + (y-1)^2\right)^{3/2}}-\frac{x}{\left(x^2 + (y+1)^2\right)^{3/2}}$$$$E_y = \frac{y-1}{\left(x^2 + (y-1)^2\right)^{3/2}}-
\frac{y+1}{\left(x^2 + (y+1)^2\right)^{3/2}}$$

$\boldsymbol{E}$ の向きを表す方向場(単位ベクトル場)

一般にベクトルは大きさと向きを表す。電場 $\boldsymbol{E}$ の大きさは電荷からの距離の2乗に反比例するので,]ベクトルの長さが大きく変化する。見やすくするために,以下のように $\boldsymbol{E}$ の向きを表す方向場,すなわち大きさが $1$ の単位ベクトル場を定義して,ベクトルの向きのみを正しく描画する。

$$\hat{\boldsymbol{E}} \equiv \frac{\boldsymbol{E}}{\sqrt{\boldsymbol{E}\cdot\boldsymbol{E}}}$$$$\hat{E}_x = \frac{E_x}{\sqrt{E_x^2 + E_y^2}}$$$$\hat{E}_y = \frac{E_y}{\sqrt{E_x^2 + E_y^2}}$$

正負の点電荷がつくる電場の方向場を描く

グラフのサイズ等の設定

グラフのサイズを設定。弘大 JupyterHub の gnuplot_kernel ではインライン SVG。sizefont の設定例。

In [1]:
%gnuplot inline svg size 640,640 fixed enhanced font 'Noto Sans CJK JP,16'

電場ベクトルの定義

電場ベクトルを x, y の2変数関数として定義。

In [2]:
# 電場 E
Ex(x, y) = x/sqrt(x**2 + y**2)**3
Ey(x, y) = y/sqrt(x**2 + y**2)**3

# 正負の点電荷がつくる電場
E2x(x, y) = Ex(x, y-1) - Ex(x, y+1)
E2y(x, y) = Ey(x, y-1) - Ey(x, y+1)

# 電場の大きさ
E(x, y) = sqrt(E2x(x, y)**2 + E2y(x, y)**2)

# 電場の方向場。規格化された hat E
hE2x(x, y) = scaling * E2x(x, y)/E(x, y)
hE2y(x, y) = scaling * E2y(x, y)/E(x, y)

ベクトル場を描く

gnuplot でベクトルを描く基本は ここ を参照。

点電荷が $y$ 軸上にあるので,ベクトルの始点が $y$ 軸に来ないように urange, vrange, samples, isosamples を設定しておく。(点電荷上では電場の分母がゼロになってしまうから。)
以下の例では,set urange [-5.5:5.5], set samples 12 なので,ベクトルの始点の $x$ 座標は

$$-5.5, -4.5, \cdots, -0.5, 0.5, \cdots, 5.5$$

となり,$y$ 軸上 ($x = 0$) に始点がこない。

special filename “++” を使ってベクトルをプロットする際の書式は…

plot "++" using 1:2:(Vx($1,$2)):(Vy($1,$2)) with vectors
In [3]:
reset
# special filename "++" のときは
# urange, vrange で描画範囲の設定
set urange [-5.5:5.5]
set vrange [-5.5:5.5]
set size ratio 1

# x 軸方向のサンプル数
set samples 12
# y 軸方向のサンプル数
set isosamples 12

# 表示範囲の設定
set xrange [-6.5:6.5]
set yrange [-6.5:6.5]
set xtics 1
set ytics 1
set grid

# ベクトルの長さ調整用
scaling = 0.7

plot "++" using 1:2:(hE2x($1,$2)):(hE2y($1,$2)) \
     with vectors lw 2 filled head notitle

点電荷の位置にラベル

正電荷,負電荷の位置にラベルをつけてみる。

In [4]:
# 正電荷の位置に青丸と +
# 
# set label 1 point pt 7 ps 2 lc "white" at 0, 1 front
set label 2 point pt 6 ps 2 lc "blue" at 0, 1 front
set label 3 center at first 0, 1 "+" tc "blue" front

# 負電荷の位置に赤丸と -
# 一旦白く塗りつぶす
set label 11 point pt 7 ps 2 lc "white" at 0, -1 front
# あらためて赤丸を描く
set label 12 point pt 6 ps 2 lc "red" at 0, -1 front
set label 13 center at first 0, -1 "ー" tc "red" front

replot

linecolor palette

lc palette (linecolor palette) オプションを指定すると,using の5つ目の値(以下では電場ベクトル $\boldsymbol{E}$ の大きさ (E($1,$2)))に対応する line color をつけてくれる。

In [5]:
plot "++" u 1:2:(hE2x($1,$2)):(hE2y($1,$2)):(E($1,$2)) \
     w vec lc palette lw 4 filled head notitle

palette オプションの設定事項

set palette cubehelix start 0 cycles 0 saturation 3 negative として,電場が弱い($E$ すなわち (E($1,$2)) が小さい)ところは薄い色になるように palette を設定してみます。

  • start 0:青,1:赤,2:緑
  • cycles 0:濃さだけが変化
  • saturation 彩度。小さい値ほど淡い。特に 0 はモノクロ。
  • negative 小さい値ほど薄く。
In [6]:
set palette cubehelix start 0 cycles 0 saturation 3 negative
plot "++" u 1:2:(hE2x($1,$2)):(hE2y($1,$2)):(E($1,$2)) \
     w vec notitle lc palette lw 4 filled head

電荷の位置から離れると,電場の大きさは急激に小さくなるので,このままでは薄くてほとんど見えないベクトルばかりになる。これだと少し寂しいので saturation を思い切って大きな値にしてみる。

In [7]:
set palette cubehelix start 0 cycles 0 saturation 300 negative
plot "++" u 1:2:(hE2x($1,$2)):(hE2y($1,$2)):(E($1,$2)) \
     w vec notitle lc palette lw 4 filled head

一応の完成作

以上の設定をまとめてみたのが以下の例。サンプル数を増やし,少し短めのベクトルをたくさん描いてみた。

In [8]:
reset
# special filename "++" のときは
# urange, vrange で描画範囲の設定
set urange [-5.25:5.25]
set vrange [-5.25:5.25]
set size ratio 1

# x 軸方向のサンプル数
set samples 22
# y 軸方向のサンプル数
set isosamples 22

# 表示範囲の設定
set xrange [-6:6]
set yrange [-6:6]

# 目盛と外枠の非表示
unset xtics
unset ytics
unset border

# 正電荷の位置に青丸と +
# set label 1 point pt 7 ps 2 lc "white" at 0, 1 front
set label 2 point pt 6 ps 2 lc "blue" at 0, 1 front
set label 3 center at first 0, 1 "+" tc "blue" front

# 負電荷の位置に赤丸と -
# 一旦白く塗りつぶす
set label 11 point pt 7 ps 2 lc "white" at 0, -1 front
# あらためて赤い輪を描く
set label 12 point pt 6 ps 2 lc "red" at 0, -1 front
set label 13 center at first 0, -1 "ー" tc "red" front

set title "正負の点電荷がつくる電場の向き"

# グラフ右のカラーボックス非表示
unset colorbox

# ベクトルの長さ調整用
scaling = 0.3
# colorbox の範囲 show cbrange で確認可能
set cbrange [0:11]

set palette cubehelix start 0 cycles 0 saturation 500 negative
plot "++" u 1:2:(hE2x($1,$2)):(hE2y($1,$2)):(E($1,$2)) \
     w vec notitle lc palette lw 4 filled head