平行2本導線による電界，電位

モジュールのインポート

import numpy as np
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'])

平行2本導線に対する2本の直線電荷

　半径$a$の円筒導体が，円筒の中心軸の間隔を$d$として平行におかれている．このような平行２本導線の静電界および電位は，一方の円筒中心軸から$b_1$，$b_2$偏位した位置に線電荷をおいて求めることができ，次のような関係を導くことがきる． \begin{gather} b_1 = \frac{a}{k}, \ \ \ \ \ b_2 = ka, \ \ \ \ \ d=b_1+b_2 \end{gather} ただし，$k$ は導体表面の電位の式に表れる定数である．これより，$b_1$を求めると， \begin{align} &b_1 b_2 = \frac{a}{k} \cdot ka = a^2 = b_1(d-b_1) \\ &b_1^2 - d b_1 + a^2 = 0 \end{align} 上式に$b_1 = d-b_2$を代入すると，次式が得られる． \begin{gather} b_2^2 - d b_2 + a^2 = 0 \end{gather} よって，$b_1$，$b_2$は次のように２次方程式の根によって得られる． \begin{gather} b_2, b_1 = \frac{d \pm \sqrt{d^2+4a^2}}{2} \end{gather} 規格化した半径が次のような場合を計算しよう． 
# Fig. 1.8(a)
a = 0.2 # 円筒導体の半径0.2, 0.5
d = 2.0 # 2つの円筒導体の中心間距離
2*a/d
0.2
このとき，$b_1$，$b_2$は， 
dd = np.sqrt(d**2-4*a**2)
b1 = (d - dd)/2.0
b2 = (d + dd)/2.0
b1, b2
(0.020204102886728803, 1.9797958971132712)

関数

def uvector3(aa):
    norm = np.linalg.norm(aa, axis=2, keepdims=True)
    uv = aa / norm
    return norm, uv
def p1charge2(r, rd, q):
    qk = q/2.0/np.pi
    rr = r-rd
    norm = np.linalg.norm(rr, axis=2, keepdims=True)
    v = -qk * np.log(norm)
    with np.errstate(divide="ignore", invalid="ignore"):
        e = qk*rr/norm**2
    return e, v

計算

2次元データを描く準備として， 
# 直角座標系(x,y)で表した観測点
x_min, x_max, n_x = -2.0, 2.0, 200
y_min, y_max, n_y = -2.0, 2.0, 200
x = np.linspace(x_min, x_max, n_x)
y = np.linspace(y_min, y_max, n_y)
xx, yy = np.meshgrid(x,y, indexing='xy') # xy座標の並び (デフォルト)
r = np.stack([xx,yy], axis=-1)
2つの平行な正・負の直線電荷による電位，電界は， 
q = np.array([1.0, -1.0])# 電荷量
nq = len(q)# 点電荷の数
rdi = np.array([-d/2+b1, 0.0])# 1つ目の点電荷の位置ベクトル(x',y')
rd = rdi.reshape(1, 1, -1)
rr = r-rd
e, vv = p1charge2(r, rd, q[0])
rdi = np.array([-d/2+b2, 0.0])# 2つ目の点電荷の位置ベクトル
rd = rdi.reshape(1, 1, -1)
ei, vi = p1charge2(r, rd, q[1])
e = e + ei
vv = vv + vi
vvv = vv.reshape(xx.shape)
eamp, vs = uvector3(e)
emax = e.max()
e_x = e[:,:,0]/emax
e_y = e[:,:,1]/emax
定数$k$は， 
k = a/b1
v1 = np.log(k)/2/np.pi
vvvn = vvv/v1
円柱導体表面の軌跡は， 
phi_min, phi_max, n_phi = 0.0, 2.0*np.pi, 181
phi = np.linspace(phi_min, phi_max, n_phi)
sp = a*np.sin(phi)
cp = a*np.cos(phi)
x_1, y_1 = -d/2+cp, sp
x_2, y_2 = d/2+cp, sp
プロットして， 
fig = plt.figure(figsize=(8, 8)) # グラフ領域の作成
v = np.linspace(-1, 1, 11)
plt.streamplot(xx, yy, e_x, e_y,  color='orange', density=0.58, linewidth=1.5)
ax = plt.contour(x, y, vvvn, levels=v, linewidths=0.8)# 等高線表示
ax.clabel(fmt='%1.1f', inline=1, fontsize=7)# 等高線の値を表示
plt.contourf(x, y, vvvn, levels=v, cmap="Blues")# coolwarm
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.axis('scaled')
plt.colorbar(pad=0.1, shrink=0.77, orientation="vertical")
plt.plot(-d/2+b1, 0.0, 'o', color = 'green')
plt.plot(-d/2+b2, 0.0, 'o', color = 'red')
plt.plot(x_1, y_1, color = 'gray')
plt.plot(x_2, y_2, color = 'gray')
plt.plot(-d/2, 0.0, '.', color = 'gray')
plt.plot(d/2, 0.0, '.', color = 'gray')
fig.savefig('Fig_1_8_a_streamplot.pdf')
plt.show()

同様にして， 
# Fig. 1.8(b)
a = 0.5 # 円筒導体の半径0.2, 0.5
d = 2.0 # 2つの円筒導体の中心間距離
2*a/d
0.5
これより，最初から計算しなおしてプロットすると，