1. はじめに

　Python三次元ポアソン方程式による静電場解析により，基本的課題として空間上の点電荷の電場解析が可能なことを先の投稿において報告した(1)．
　本稿では，応用課題として一様に帯電した球を対象に下記を実施する．
(1)三次元ポアソン方程式を差分化（再掲）
(2)繰り返し法によるプログラムを作成．
①メッシュ幅は1.0と0.5を選択できるようにする．
②電位を計算
③得られた電位から電場を計算
(3)計算結果を理論値と比較することにより解析の妥当性を検証

2. 例題の問題設定

各軸のメッシュ数は50とし，中心に一様に帯電した球を配置する．

　　　　　　図1　問題設定-一様に帯電した球の配置

3. ポアソン方程式の離散化(再掲)

　下記に手順を記す．
　ここで注意が必要な点はρは電荷密度であることである．

4. プログラム

import pprint
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.cm as cm

delta_L=0.5
#delta_L=1.0
LN=50
HLN=int(LN/2)
nx = LN
ny = LN
nz = LN
xmin = 0
xmax = LN*delta_L
ymin = 0
ymax = LN*delta_L
zmin = 0
zmax = LN*delta_L

p  = np.zeros((nz, ny, nx))
pd = np.zeros((nz, ny, nx))
charge = np.zeros((nz, ny, nx))
x  = np.linspace(xmin, xmax, nx)
y  = np.linspace(ymin, ymax, ny)
z  = np.linspace(zmin, zmax, nz)

eps0=1
y00=1
Q1=1/delta_L**3

for i in range(0, LN):
    for j in range(0,LN):
        for k in range(0,LN):
            rr=np.sqrt((i-HLN)**2+(j-HLN)**2+(k-HLN)**2)
            if rr <= 3.5 :
                charge[i,j,k] = Q1
            elif rr <= 4.5:
                charge[i,j,k] = Q1*(4.5-rr)
#
print("Riron")
for i in range(1,5):
    r=i*delta_L
    x2=-Q1*r**2/6+Q1*(4*delta_L)**2/2
    print('z=',r,' V=',x2)
for i in range(4,10):
    r=i*delta_L
    x2=Q1*(4*delta_L)**3/(3*r)
    print('z=',r,' V=',x2)    
for i in range(1,5):
    r=i*delta_L
    x2=Q1*r/3
    print('z=',r,'Ez=',x2)
for i in range(4,10):
    r=i*delta_L
    x2=Q1*(4*delta_L)**3/3/r**2
    print('z=',r,'Ez=',x2)

for I in range(500): #2000
    pd = p.copy()
    p[1:-1,1:-1,1:-1] = (pd[2:,1:-1,1:-1] + pd[:-2,1:-1,1:-1] +
                         pd[1:-1,2:,1:-1] + pd[1:-1,:-2,1:-1] + 
                         pd[1:-1,1:-1,2:] + pd[1:-1,1:-1,:-2] + (delta_L**2)*charge[1:-1, 1:-1, 1:-1]) /6 

    p[0][:][:]=0
    p[nx-1][:][:]=0
    p[:][0][:]=0
    p[:][ny-1][:]=0
    p[:][:][0]=0
    p[:][:][nz-1]=0

def plot2D(x, y, p):
    fig = plt.figure(figsize=(11, 7), dpi=100)
    ax = fig.gca(projection='3d')
    X, Y = np.meshgrid(x, y)
    surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False)
    ax.view_init(30, 225)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('p')

plot2D(x,y,p[:,:,HLN])

def plot_surface(pp):
    nx = LN
    ny = LN
    fig = plt.figure(figsize=(11,7), dpi=100)
    fig    
    xmin = 0
    xmax = LN-1
    ymin = 0
    ymax = LN-1
    x  = np.linspace(xmin, xmax, nx)
    y  = np.linspace(ymin, ymax, ny)
    X, Y = np.meshgrid(x, y)
    plt.contourf(X, Y, pp, alpha=0.1, cmap=cm.viridis)
    plt.colorbar()
    plt.xlabel('X')
    plt.ylabel('Y')

plot_surface(p[:,:,HLN])
plt.xlim([HLN-5,HLN+5])
plt.ylim([HLN-5,HLN+5])

plt.style.use('ggplot')
plt.rcParams["axes.facecolor"] = 'white'
fig = plt.figure()
ax = fig.gca(projection='3d')

L=9
cs=p[HLN-4:HLN+8,HLN-4:HLN+8,HLN-4:HLN+8]

X, Y, Z= np.meshgrid(np.arange(0, L,1), np.arange(0, L,1), np.arange(0, L,1))  

Ex = np.zeros([L,L,L])
Ey = np.zeros([L,L,L])
Ez = np.zeros([L,L,L])
EE = np.zeros([L,L,L])

for i in range(L):
  for j in range(L):
    for k in range(L):
      Ex[i,j,k]=-(cs[i,j+1,k]-cs[i,j-1,k])
      Ey[i,j,k]=-(cs[i+1,j,k]-cs[i-1,j,k])
      Ez[i,j,k]=-(cs[i,j,k+1]-cs[i,j,k-1]) 

for i in range(L):
  for j in range(L):
    for k in range(L):
      ii=i*2+HLN-8
      jj=j*2+HLN-8
      kk=k*2+HLN-8     

      Ex[i,j,k]=-(p[ii,jj+1,kk]-p[ii,jj-1,kk])/2/delta_L#*2#*0.004#*10
      Ey[i,j,k]=-(p[ii+1,jj,kk]-p[ii-1,jj,kk])/2/delta_L#*2#*0.004#*10
      Ez[i,j,k]=-(p[ii,jj,kk+1]-p[ii,jj,kk-1])/2/delta_L#*2#*0.004#*10

print("V=")
for i in range(1,12):
  print('z=',i*delta_L,'Vz=',p[HLN,HLN,HLN+i])

i=25
k=25
for j in range(26,35):
    print("x=",(j-25)*delta_L,"Ex=",-(p[i,j+1,k]-p[i,j-1,k])/2/delta_L)

ax.quiver(X,Y,Z, Ex, Ey, Ez, color='red',length=0.5*delta_L**2, normalize=False)
#ax.quiver(X,Y,Z, Ex, Ey, Ez, color='red',length=0.2*(delta_L**2), normalize=False)

for i in range(-2,2):
    for j in range(-2,2):
        for k in range(-2,2):
            X1,Y1, Z1=int(L/2)+i,int(L/2)+j,int(L/2)+k
            ax.scatter3D(X1,Y1,Z1,"o", color='blue')

for i in range(5):
    X1,Y1, Z1=int(L/2)-i,int(L/2),int(L/2)
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    X1,Y1, Z1=int(L/2)+i,int(L/2),int(L/2)
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    X1,Y1, Z1=int(L/2),int(L/2)-i,int(L/2)
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    X1,Y1, Z1=int(L/2),int(L/2)+i,int(L/2)
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    X1,Y1, Z1=int(L/2),int(L/2),int(L/2)-i
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    X1,Y1, Z1=int(L/2),int(L/2),int(L/2)+i
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')

# グラフ描画
plt.grid()
plt.draw()

plt.show()