多数の点電荷による電位，電界¶
In [1]:
from json.encoder import INFINITY
import numpy as np
import matplotlib.pyplot as plt
import scienceplots
plt.style.use(['science', 'notebook'])
#plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['font.family'] = 'serif' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定
クーロンの法則は，$N$個の点電荷$Q_i$ [C]が位置ベクトル$\boldsymbol{r}_i$ [m]（$i=1,2, \cdots, N$）にあるとき，観測点$\boldsymbol{r}$ [m]の静電界$\boldsymbol{E}_i$ [N/C]は， \begin{gather} \boldsymbol{E} = \sum_{i=1}^N \boldsymbol{E}_i, \ \ \ \ \ \boldsymbol{E}_i = \frac{1}{4 \pi \epsilon_0} \frac{Q_i(\boldsymbol{r}-\boldsymbol{r}_i)}{|\boldsymbol{r}-\boldsymbol{r}_i|^3} \end{gather} ここで， \begin{gather} \frac{1}{4 \pi \epsilon_0} = 9 \times 10^9 \end{gather} ただし，$\epsilon_0= 8.842 \times 10^{-12}$ [F/m]は真空中の誘電率を示す．また，電位$V$は， \begin{gather} V = \frac{1}{4 \pi \epsilon_0} \sum_{i=1}^N \frac{Q_i}{R_i}, \ \ \ \ \ R_i = |\boldsymbol{r}-\boldsymbol{r}_i| \end{gather}
In [2]:
def p1charge(r, rd, q):
qk = q/4.0/np.pi
rr = r-rd
norm = np.linalg.norm(rr, axis=2, keepdims=True)
with np.errstate(divide="ignore", invalid="ignore"):
v = qk / norm
e = qk*rr/norm**3
return e, v
In [3]:
def pncharge(r, rdn, q):
rd = rdn[0].reshape(1, 1, -1)
e, v = p1charge(r, rd, q[0])
if len(q)>=2:
for i in range(len(q)-1):
rd = rdn[i+1].reshape(1, 1, -1)
ei, vi = p1charge(r, rd, q[i+1])
e = e + ei
v = v + vi
return e,v
In [4]:
# Fig.1.3(b) に対応
q1 = np.array([1.0, -1.0]) # 点電荷の電荷量
rd1_p = np.array([[0.0, 1.0, 0.0]]) # 正の点電荷の位置ベクトル
rd1_m = np.array([[0.0, -1.0, 0.0]]) # 負の点電荷の位置ベクトル
rd1 = np.concatenate([rd1_p, rd1_m])
q1, rd1
Out[4]:
(array([ 1., -1.]),
array([[ 0., 1., 0.],
[ 0., -1., 0.]]))
In [5]:
# 直角座標系(x,y)で表した観測点
x = np.linspace(-3, 3, 13)
y = np.linspace(-3, 3, 13)
xx, yy = np.meshgrid(x,y, indexing='xy')# xy座標の並び (デフォルト)
zz = np.zeros(xx.shape)
r = np.stack([xx,yy,zz], axis=-1)
In [6]:
r00 = np.array([0.0, 0.0, 0.0])
e00, vv00 = pncharge(r00, rd1, q1)
r0 = np.array([0.0, 0.9, 0.0])
e0, vv0 = pncharge(r0, rd1, q1)
e1, vv1 = pncharge(r, rd1, q1)
e_x = e1[:,:,0]/e00[0,0,1]
e_y = e1[:,:,1]/e00[0,0,1]
vvv = vv1[:,:,0]/vv0[0:,0:,0]
In [7]:
def table_xy(ss,x1,x2,aaa,nc):
n1,n2 = np.shape(aaa)
n_n2 = n2//nc+1
for m in range(n_n2):
print(f'({m+1:>3d}/{n_n2:>3d})')
if n2 > nc*(m+1):
m2 = nc
else:
m2 = np.mod(n2,nc)
print(f'{ss:>8s}',end='')
[print(f'{x2[j+m*nc]:>8.3f}',end='') for j in range(m2)]
print('')
for i in range(n1):
print(f'{x1[i]:>8.3f}',end='')
[print(f'{aaa[i,j+m*nc]:8.3f}',end='') for j in range(m2)]
print('')
print('')
return
In [8]:
table_xy('y | x',y,x,e_x,nc=11-7)
( 1/ 4) y | x -3.000 -2.500 -2.000 -1.500 -3.000 -0.020 -0.026 -0.033 -0.038 -2.500 -0.024 -0.035 -0.049 -0.065 -2.000 -0.028 -0.043 -0.068 -0.108 -1.500 -0.028 -0.047 -0.084 -0.159 -1.000 -0.024 -0.042 -0.081 -0.174 -0.500 -0.014 -0.025 -0.050 -0.111 0.000 -0.000 -0.000 -0.000 -0.000 0.500 0.014 0.025 0.050 0.111 1.000 0.024 0.042 0.081 0.174 1.500 0.028 0.047 0.084 0.159 2.000 0.028 0.043 0.068 0.108 2.500 0.024 0.035 0.049 0.065 3.000 0.020 0.026 0.033 0.038 ( 2/ 4) y | x -1.000 -0.500 0.000 0.500 -3.000 -0.038 -0.025 -0.000 0.025 -2.500 -0.075 -0.058 -0.000 0.058 -2.000 -0.161 -0.170 -0.000 0.170 -1.500 -0.332 -0.692 -0.000 0.692 -1.000 -0.455 -1.971 nan 1.971 -0.500 -0.272 -0.644 -0.000 0.644 0.000 -0.000 -0.000 -0.000 -0.000 0.500 0.272 0.644 -0.000 -0.644 1.000 0.455 1.971 nan -1.971 1.500 0.332 0.692 -0.000 -0.692 2.000 0.161 0.170 -0.000 -0.170 2.500 0.075 0.058 -0.000 -0.058 3.000 0.038 0.025 -0.000 -0.025 ( 3/ 4) y | x 1.000 1.500 2.000 2.500 -3.000 0.038 0.038 0.033 0.026 -2.500 0.075 0.065 0.049 0.035 -2.000 0.161 0.108 0.068 0.043 -1.500 0.332 0.159 0.084 0.047 -1.000 0.455 0.174 0.081 0.042 -0.500 0.272 0.111 0.050 0.025 0.000 -0.000 -0.000 -0.000 -0.000 0.500 -0.272 -0.111 -0.050 -0.025 1.000 -0.455 -0.174 -0.081 -0.042 1.500 -0.332 -0.159 -0.084 -0.047 2.000 -0.161 -0.108 -0.068 -0.043 2.500 -0.075 -0.065 -0.049 -0.035 3.000 -0.038 -0.038 -0.033 -0.026 ( 4/ 4) y | x 3.000 -3.000 0.020 -2.500 0.024 -2.000 0.028 -1.500 0.028 -1.000 0.024 -0.500 0.014 0.000 -0.000 0.500 -0.014 1.000 -0.024 1.500 -0.028 2.000 -0.028 2.500 -0.024 3.000 -0.020
In [9]:
table_xy('y | x',y,x,e_y,nc=11-7)
( 1/ 4) y | x -3.000 -2.500 -2.000 -1.500 -3.000 -0.005 -0.011 -0.022 -0.038 -2.500 -0.002 -0.008 -0.021 -0.047 -2.000 0.004 -0.000 -0.013 -0.046 -1.500 0.012 0.013 0.010 -0.013 -1.000 0.021 0.030 0.044 0.064 -0.500 0.029 0.045 0.077 0.142 0.000 0.032 0.051 0.089 0.171 0.500 0.029 0.045 0.077 0.142 1.000 0.021 0.030 0.044 0.064 1.500 0.012 0.013 0.010 -0.013 2.000 0.004 -0.000 -0.013 -0.046 2.500 -0.002 -0.008 -0.021 -0.047 3.000 -0.005 -0.011 -0.022 -0.038 ( 2/ 4) y | x -1.000 -0.500 0.000 0.500 -3.000 -0.061 -0.084 -0.094 -0.084 -2.500 -0.092 -0.150 -0.181 -0.150 -2.000 -0.129 -0.304 -0.444 -0.304 -1.500 -0.115 -0.632 -1.920 -0.632 -1.000 0.089 0.114 nan 0.114 -0.500 0.307 0.897 2.222 0.897 0.000 0.354 0.716 1.000 0.716 0.500 0.307 0.897 2.222 0.897 1.000 0.089 0.114 nan 0.114 1.500 -0.115 -0.632 -1.920 -0.632 2.000 -0.129 -0.304 -0.444 -0.304 2.500 -0.092 -0.150 -0.181 -0.150 3.000 -0.061 -0.084 -0.094 -0.084 ( 3/ 4) y | x 1.000 1.500 2.000 2.500 -3.000 -0.061 -0.038 -0.022 -0.011 -2.500 -0.092 -0.047 -0.021 -0.008 -2.000 -0.129 -0.046 -0.013 -0.000 -1.500 -0.115 -0.013 0.010 0.013 -1.000 0.089 0.064 0.044 0.030 -0.500 0.307 0.142 0.077 0.045 0.000 0.354 0.171 0.089 0.051 0.500 0.307 0.142 0.077 0.045 1.000 0.089 0.064 0.044 0.030 1.500 -0.115 -0.013 0.010 0.013 2.000 -0.129 -0.046 -0.013 -0.000 2.500 -0.092 -0.047 -0.021 -0.008 3.000 -0.061 -0.038 -0.022 -0.011 ( 4/ 4) y | x 3.000 -3.000 -0.005 -2.500 -0.002 -2.000 0.004 -1.500 0.012 -1.000 0.021 -0.500 0.029 0.000 0.032 0.500 0.029 1.000 0.021 1.500 0.012 2.000 0.004 2.500 -0.002 3.000 -0.005
In [10]:
# 直角座標系(x,y)で表した観測点
nxy = 201
x = np.linspace(-3, 3, nxy)
y = np.linspace(-3, 3, nxy)
xx, yy = np.meshgrid(x,y, indexing='xy')# xy座標の並び (デフォルト)
zz = np.zeros(xx.shape)
r = np.stack([xx,yy,zz], axis=-1)
In [11]:
r0 = np.array([0.0, 0.9, 0.0])
e0, vv0 = pncharge(r0, rd1, q1)
e1, vv1 = pncharge(r, rd1, q1)
e_x = e1[:,:,0]
e_y = e1[:,:,1]
vvv = vv1[:,:,0]/vv0[0:,0:,0]
In [12]:
fig = plt.figure(figsize=(8, 8)) # 図のサイズ指定
v = np.linspace(-1, 1, 21) # 電位の等高線のレベル設定
# 電場ベクトルの力線（オレンジ）
plt.streamplot(xx, yy, e_x, e_y, color='orange', density=1.0, linewidth=1.5)
# 電位の等高線（黒線）
contour = plt.contour(x, y, vvv, levels=v, linewidths=0.8, colors='black')
contour.clabel(fmt='%1.1f', inline=True, fontsize=7)
# 電位のカラーマップ（青〜赤：bwr）
plt.contourf(x, y, vvv, levels=v, cmap="seismic") # "bwr", "coolwarm", "seismic", "twilight"
# 軸のスケーリングとカラーバー
plt.axis('scaled')
plt.colorbar(pad=0.1, shrink=0.77, orientation="vertical")
# 点電荷の描画（緑：負電荷、赤：正電荷）
plt.plot(rd1_m[:,0], rd1_m[:,1], 'o', color='green') # 負電荷
plt.plot(rd1_p[:,0], rd1_p[:,1], 'o', color='red') # 正電荷
# 図の保存と表示
fig.savefig('e_v_fig_1_3_b.pdf')
plt.show()
In [13]:
fig = plt.figure() # グラフ領域の作成
plt.xlabel("$x, y$") # x軸のラベル設定
plt.ylabel("V") # y軸のラベル設定
plt.ylim(-4, 4) # y軸範囲の設定
plt.plot(r[:,nxy//2,1], vvv[:,nxy//2], label='$x=0$')
plt.plot(r[nxy//2,:,0], vvv[nxy//2,:], label='$y=0$')
plt.legend(ncol=1, loc='lower right', fancybox=False, frameon=True)
fig.savefig('e_v_fig_1_3_b_xy.pdf')
plt.show()
