我一直在尝试编写一个程序来分析电介质球对项目的电磁波散射。
参考文献是 Sadiku 的《电磁学中的数值方法》第 3 版。
现在这本书在 Matlab 中有代码，但我必须使用 Python。
我用 Python 重写了代码并运行了它，但是值 Ey/Einc 和 Ez/Einc 爆炸了，当我增加间隔数时更是如此。即使在 n = 50 时，它们也与结果不匹配。
这是我的代码
'''
import numpy as np
from matplotlib import pyplot as py
import warnings
'''
The followings constants are defined globally, as they are used in the code.
All are in SI units.
'''
pi = np.pi
c = 3e8
e0 = 8.85e-12
mu0 = 4*pi*1e-7
'''
DOCSTRING for Taflove(er, mu, sigma, nu, num)
This funciton is used to calculate some parameters and Taflove's constants.
The inputs are relative permittivity and permeability, conductivity,
frequency of incoming wave and grid number.
Three tuples are returned.
First tuple contains lambd, delta and delta_t parameters.
Second tuple contains R constants.
Third tuple contains C constants.
'''
def Taflove(er, mur, sigma, nu, num):
lambd = c/nu
delta = lambd/num
delta_t = delta/(2*c)
R = delta_t/e0
R_a = (c*delta_t/delta)**2
R_b = delta_t/(mu0*delta)
C_a = (1 - R*(sigma/er))
C_b = R_a/er
return np.array((lambd, delta, delta_t)), np.array((R, R_a, R_b)), np.array((C_a, C_b))
'''
DOCSTRING for
'''
er = (1, 4) # Relative Permittivity
mur = (1, 1) # Relative Permeability (1 for non-magnetic)
sigma = (0.1, 0) # Conductivity (0 for perfect dielectric/lossless medium)
# Conductivity at boundary is taken to be non zero to prevent abrupt results
nu = 2.5e9 # Frequency of EM wave = 2.5 GHz (for Microwave)
num = 40 # Due to computing power
radius = 4.5e-2 # Radius of sphere is 4.5 cm
P_list, R_list, C_list1 = Taflove(er[0], mur[0], sigma[0], nu, num)
print(P_list)
print(R_list)
print(C_list1)
input('check1')
_, _, C_list2 = Taflove(er[1], mur[1], sigma[1], nu, num)
C = [C_list1, C_list2]
R, R_a, R_b = R_list
lambd, delta, delta_t = P_list
#1 = air
#2 = sphere
#E = np.array([0, 0, 1]) # Incoming EM wave is linearly polarized in z direction
#B = np.array([1, 0, 0])/c # B = k cap cross E/c (in x direction)
#k = 2*pi*np.array([0, 1, 0])/P_list[0] # k = 2pi/lambda in direction of wave (in y direction)
Size = (19, 39, 19) # Grid/Lattice size. More in y because wave is propagating in y driection
Center = (19.5, 20, 19) # Center of sphere
WMax = 2 # Steps in one program looop
TMax =20*5 # Number of Timesteps
Media = 2 # Number of Media
JPW = 3 # y coordinate of wave
I, J, K = Size
Ex = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
Hx = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
Ey = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
Hy = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
Ez = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
Hz = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
x_ = np.linspace(0, I+1, I+2)
y_ = np.linspace(0, J+1, J+2)
z_ = np.linspace(0, K+1, K+2)
print(x_)
print(y_)
print(z_)
radius_unit = radius/P_list[1]
print(radius_unit)
input('check')
pos_func = lambda R : np.sqrt(sum(np.array(R) - np.array(Center))**2)
XMed = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2))
YMed = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2))
ZMed = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2))
for i in range(Size[0]+2):
for j in range(Size[1] + 2):
for k in range(Size[2]+2):
if pos_func((i + 0.5, j, k))<radius_unit:
XMed[i, j, k] = 2
else:
XMed[i, j, k] = 1
if pos_func((i, j+ 0.5, k))<radius_unit:
YMed[i, j, k] = 2
else:
YMed[i, j, k] = 1
if pos_func((i, j, k+ 0.5))<radius_unit:
ZMed[i, j, k] = 2
else:
ZMed[i, j, k] = 1
print(XMed)
Ey1 = np.zeros(J+2)
Ez1 = np.zeros(J+2)
print(len(Ey1))
print("Setting up completed...")
input('a')
###########################
t3 = 2
t2 = 1
t1 = 0
for n in range(TMax):
if n%10 == 0:
print(str(100*n/TMax) + " % complete.")
#Applying soft truncation conditons
for k in range(K+1):
for j in range(J+1):
for i in range(I+1):
if i == 0:
if k != K and k!=0:
Hy[0, j, k, t3] = (Hy[1, j, k-1, t1] + Hy[1, j, k, t1] + Hy[1, j, k+1, t1])/3
Hz[0, j, k, t3] = (Hz[1, j, k-1, t1] + Hz[1, j, k, t1] + Hz[1, j, k+1, t1])/3
elif k == K:
Hy[0, j, k, t3] = (Hy[1, j, k-1, t1] + Hy[1, j, k, t1])/2
Hz[0, j, k, t3] = (Hz[1, j, k-1, t1] + Hz[1, j, k, t1])/2
else:
Hy[0, j, k, t3] =(Hy[1, j, k, t1] + Hy[1, j , k+1, t1])/2
Hz[0, j, k, t3] = (Hz[1, j, k, t1] + Hz[1, j, k+1, t1])/2
if j == 0:
Ex[i, 0, k, t3] = Ex[i, 1, k, t1]
Ez[i, 0, k, t3] = Ez[i, 1, k, t1]
elif j == J:
Ex[i, j, k, t3] = Ex[i, j-1, k, t1]
Ez[i, j, k, t3] = Ez[i, j-1, k, t1]
if k == 0:
if i != I and i!= 0:
Ex[i, j, 0, t3] = (Ex[i-1, j, 1, t1] + Ex[i, j, 1, t1] + Ex[i+1, j, 1, t1])/3
Ey[i, j, 0, t3] = (Ey[i-1, j, 1, t2] + Ey[i, j, 1, t2] + Ey[i+1, j, 1, t1])/3
elif i == 0:
Ex[0, j, 0, t3] = (Ex[i, j, 1, t1] + Ex[i+1, j, 1, t1])/2
Ey[0, j, 0, t3] = (Ey[i, j, 1, t1] + Ey[i+1, j, 1, t1])/2
else:
Ex[i, j, 0, t3] = (Ex[i, j, 1, t1] + Ex[i-1, j, 1, t1])/2
Ey[i, j, 0, t3] = (Ey[i, j, 1, t1] + Ey[i-1, j, 1, t1])/2
#Applying Yee's algorithm
Hx[i, j, k, t3] = Hx[i, j, k, t2] + R_b*(Ey[i, j, k+1, t2] - Ey[i, j, k, t2]\
+ Ez[i, j, k, t2] - Ez[i, j+1, k, t2])
Hy[i, j, k, t3] = Hy[i, j, k, t2] + R_b*(Ex[i, j, k, t2] - Ex[i, j, k+1, t2]\
+ Ez[i+1, j, k, t2] - Ez[i, j, k, t2])
Hz[i, j, k, t3] = Hz[i, j, k, t2] + R_b*(Ex[i, j+1, k, t2] - Ex[i, j, k, t2]\
+ Ey[i, j, k, t2] - Ey[i+1, j, k, t2])
if k == K:
Hx[i, j, k, t3] = Hx[i, j, k-1, t3]
Hy[i, j, k, t3] = Hy[i, j, k-1, t3]
if j !=0 and j !=J and k !=0:
M = int(XMed[i, j, k] - 1)
Ex[i, j, k, t3] = C[M][0]*Ex[i, j, k, t2] + (C[M][1]/R_list[2]) *(Hz[i, j, k, t3]\
-Hz[i, j-1, k, t3] + Hy[i, j, k-1, t3] - Hy[i, j, k, t3])
if k!=0:
M = int(YMed[i, j, k] - 1)
if i!=0:
Ey[i, j, k, t3] = C[M][0]*Ey[i, j, k, t2] + (C[M][1]/R_list[2]) *(Hz[i-1, j, k, t3]\
-Hz[i, j, k, t3] + Hx[i, j, k, t3] - Hx[i, j, k-1, t3])
if Ey[i, j, k, t3]!=0:
#print(Ey[i, j, k, t3])
#input(str(i) +','+ str(j)+',' + str(k))
pass
else:
Ey[i, j, k, t3] = C[M][0]*Ey[i, j, k, t2] + (C[M][1]/R_list[2]) *(0\
-Hz[i, j, k, t3] + Hx[i, j, k, t3] - Hx[i, j, k-1, t3])
if j !=0 and j!=J:
M = int(ZMed [i, j, k])-1
if M == 0:
Cam = 1
else:
Cam = C[M][0]
if i !=0:
Ez [i, j, k,t3] = Cam*Ez[i, j, k, t2]\
+ (C[M][1]/R_list[2]) *(Hy[i, j, k, t3]\
-Hy[i-1, j, k, t3] + Hx[i, j-1, k, t3] - Hx[i, j, k, t3])
else:
Ez [i, j, k,t3] = Cam*Ez[i, j, k, t2]\
+ (C[M][1]/R_list[2]) \
*(Hy[i, j, k, t3]\
-0 + Hx[i, j-1, k, t3]- Hx[i, j, k, t3])
if j == JPW:
Ez[i,j,k, t3] = Ez[i,j,k,t3] + np.sin(2*pi*nu*delta_t*(n+1))
if k==K:
Ex[i, j, k + 1, t3] = Ex[i, j, k , t3]
Ey[i, j, k + 1, t3] = Ey[i, j, k , t3]
if k == K and n>TMax-num-1:
temp = abs(Ey[I, j, K-1, t3])
if temp>Ey1[j]:
Ey1[j] = temp#/(Ey1[j]+0.0000001)
temp = abs(Ez[I, j, K, t3])
if temp > Ez1[j]:
Ez1[j] =temp#/(Ez1[j]+1e-10)
t1 = int(t2 + 0.)
t2 = int(t3 + 0.)
t3 = int(t3%3)
Integr = np.linspace(1, 100, 100)
print(Ey1)
for i in range(len(Ey1)):
if Ey1[i]<1:
py.plot(Integr[i], Ey1[i], 'r*')
py.xlabel('j')
py.ylabel('Computer Ey/Einc ')
py.grid()
py.show()
print(Ez1)
for i in range(len(Ez1)):
if Ez1[i]<1:
py.plot(Integr[i], Ez1[i], 'r*')
py.xlabel('j')
py.ylabel('Computer Ez/Einc ')
py.grid()
py.show()
有谁知道可能是什么错误？我只在上周检查过这个。谢谢。