p3_tlt

直線電荷による電気力線

直線電荷による電界，電位

　$N$個の$z$方向の線電荷 $Q_i$ [C]が位置ベクトル $\boldsymbol{\rho}_i$ [m]（$i=1,2, \cdots, N$）にあるとき，観測点 $\boldsymbol{\rho}$ [m]の静電界$\boldsymbol{E}_i$ [N/C]は，クーロンの法則より次式で与えられる． \begin{gather} \boldsymbol{E} = \sum_{i=1}^N \boldsymbol{E}_i, \ \ \ \ \ \boldsymbol{E}_i = \frac{1}{2 \pi \epsilon_0} \frac{Q_i(\boldsymbol{\rho}-\boldsymbol{\rho}_i)}{|\boldsymbol{\rho}-\boldsymbol{\rho}_i|^2} \end{gather} ここで，$\epsilon_0$ [F/m]は真空中の誘電率を示し， \begin{gather} \frac{1}{4 \pi \epsilon_0} = 9 \times 10^9, \ \ \ \ \ \epsilon_0= 8.842 \times 10^{-12} \end{gather} また，電位$V$は， \begin{gather} V = \frac{1}{2 \pi \epsilon_0} \sum_{i=1}^N Q_i \log_e \left( \frac{1}{R_i} \right), \ \ \ \ \ R_i = |\boldsymbol{\rho}-\boldsymbol{\rho}_i| \end{gather}

モジュールのインポート

import numpy as np  
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

import scienceplots
#Warning : As of version 2.0.0, you need to add import scienceplots before setting the style (plt.style.use('science')).
plt.style.use(['science', 'notebook'])

ユーザ関数

　任意のベクトルから単位ベクトル，ノルムを求めるため，

def uvector(aa): # 多次元配列で定義したベクトルの振幅，単位ベクトル
    norm = np.linalg.norm(aa)
    uv = aa / norm
    return norm, uv

1つの直線電荷による電界，電位，および観測点と電荷のある点との距離は，

def p1charge_line_2d(r, rd, q):# 1つの直線電荷による電界，電位
    qk = q/2.0/np.pi
    rr = r-rd
    norm = np.linalg.norm(rr)
    v = -qk * np.log(norm)
    with np.errstate(divide="ignore", invalid="ignore"):
        e = qk*rr/norm**2
    return e, v, norm

複数の直線電荷による電界，電位は，

def func_eline_q(s, r, q, rd):
    e = np.zeros(2)
    v = 0
    R = np.empty(0)
    for i in range(len(q)):
        ei, vi, Ri = p1charge_line_2d(r, rd[i], q[i])
        e = e + ei
        v = v + vi
        R = np.append(R, Ri)
    eamp, vs = uvector(e)
    R_min = np.min(R)
    if(R_min < 5.0e-3):
        vs = vs*0.0
    return vs

電気力線

　微分方程式の変数，初期の設定に必要な値を準備して，

s_ini, s_end, ns = 0.0, 5, 101
s = np.linspace(s_ini, s_end, ns)
nphi = 25
phi = np.linspace(0.0, 2.0*np.pi, nphi)
d_eps = 0.05
d_x, d_y = d_eps*np.cos(phi), d_eps*np.sin(phi)

微分方程式の計算を含めて，電気力線をプロットする関数は，

def e_field_lines_2(q, rd, s, d_x, d_y, nphi):
    nq = len(q)
    ns = len(s)
    s_ini = s[0]
    x_end = s[-1]
    ode = np.empty([ns,2,nphi*nq])
    ij = 0
    for j in range(nq):
        y0 = rd[j,:]
        for i in range(nphi):
            yeps = np.array([d_x[i], d_y[i]])
            y_ini = y0 + yeps
            if q[j]>0: # プラスの電荷から計算した電気力線の軌跡
                ode2 = solve_ivp(func_eline_q, [s_ini, s_end], y_ini, \
                                 dense_output=True, args=(q,rd,))
                y = ode2.sol(s)
            else: # マイナスの電荷から計算した電気力線の軌跡
                ode2 = solve_ivp(func_eline_q, [s_ini, -s_end], y_ini, \
                                 dense_output=True, args=(q,rd,))
                y = ode2.sol(-s)
            yt = np.transpose(y)
            plt.plot(yt[:,0], yt[:,1], linewidth=0.5, color="RebeccaPurple")
            ode[:,:,ij] = yt
            ij = ij+1
    return ode

テキストの図1.3(b)の点電荷を直線電荷を置き換えて，

# Fig.1.3(b) に対応
q = np.array([1.0, -1.0]) # 点電荷の電荷量
rd_p = np.array([[0.0, 1.0]]) # 正の点電荷の位置ベクトル
rd_m = np.array([[0.0, -1.0]]) # 負の点電荷の位置ベクトル
rd = np.concatenate([rd_p, rd_m])
rd

array([[ 0.,  1.],
       [ 0., -1.]])

プロットすると，

fig = plt.figure() # グラフ領域の作成
plt.axes().set_aspect('equal')
plt.xlim(-3.0, 3.0) # x軸範囲の設定
plt.ylim(-3.0, 3.0) # y軸範囲の設定
plt.xticks(np.linspace(-3.0, 3.0, 7))
plt.yticks(np.linspace(-3.0, 3.0, 7))
plt.xlabel("$x$") # 横軸のラベル設定
plt.ylabel("$y$") # 縦軸のラベル設定
#ode = e_field_lines_1(q, rd, s, d_x, d_y, nphi) # odeint
ode = e_field_lines_2(q, rd, s, d_x, d_y, nphi) # solve_ivp
plt.plot(rd_m[:,0], rd_m[:,1],'o')
plt.plot(rd_p[:,0], rd_p[:,1],'o')
fig.savefig('e_lines_1_3_b.pdf')
plt.show()

テキストの図1.4(a)の点電荷を直線電荷を置き換えて，

# Fig.1.4(a) に対応
q = np.array([1.0, 1.0, 1.0, -1.0, -1.0, -1.0])
# 正の点電荷の位置ベクトル
rd_p = np.array([[-1.0, 1.0],[0.0, 1.0],[1.0, 1.0]])
# 負の点電荷の位置ベクトル
rd_m = np.array([[-1.0, -1.0],[0.0, -1.0],[1.0, -1.0]])
rd = np.concatenate([rd_p, rd_m])
rd

array([[-1.,  1.],
       [ 0.,  1.],
       [ 1.,  1.],
       [-1., -1.],
       [ 0., -1.],
       [ 1., -1.]])

プロットして，

fig = plt.figure() # グラフ領域の作成
plt.axes().set_aspect('equal')
plt.xlim(-3.0, 3.0) # x軸範囲の設定
plt.ylim(-3.0, 3.0) # y軸範囲の設定
plt.xticks(np.linspace(-3.0, 3.0, 7))
plt.yticks(np.linspace(-3.0, 3.0, 7))
ode = e_field_lines_2(q, rd, s, d_x, d_y, nphi) # solve_ivp
plt.plot(rd_m[:,0], rd_m[:,1],'o')
plt.plot(rd_p[:,0], rd_p[:,1],'o')
fig.savefig('e_lines_1_4_a.pdf')

テキストの図1.4(b)の点電荷を直線電荷を置き換えて，

# Fig.1.4(b) に対応
q = np.array([1.0, 1.0, -1.0, -1.0])
rd_p = np.array([[2.0, 0.0],[1.0, 0.0]])
rd_m = np.array([[-1.0, 0.0],[-2.0, 0.0]])
rd = np.concatenate([rd_p, rd_m])
rd

array([[ 2.,  0.],
       [ 1.,  0.],
       [-1.,  0.],
       [-2.,  0.]])

fig = plt.figure() # グラフ領域の作成
plt.axes().set_aspect('equal')
plt.xlim(-3.0, 3.0) # x軸範囲の設定
plt.ylim(-3.0, 3.0) # y軸範囲の設定
plt.xticks(np.linspace(-3.0, 3.0, 7))
plt.yticks(np.linspace(-3.0, 3.0, 7))
ode = e_field_lines_2(q, rd, s, d_x, d_y, nphi) # solve_ivp
plt.plot(rd_m[:,0], rd_m[:,1],'o')
plt.plot(rd_p[:,0], rd_p[:,1],'o')
fig.savefig('e_lines_1_4_b.pdf')