I have been trying to write a program that analyses EM wave scattering by a dielectric sphere for a project.
The reference is Sadiku's book Numerical Methods in electromagnetics Edition 3.
Now the book has code in Matlab but I have to use Python.
I rewrote the code in Python and ran it, but the values Ey/Einc and Ez/Einc blow up, much more so when I increase number of intervals. Even at n = 50 they do not match the results.
Here is my code
'''
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()
Does anyone know what could be the bug? I have been checking this only for the last week. Thanks.
np.sin(2*pi*nu*delta_t*(n+1)). How big in magnitude can the input to the sine function get there? Generally speaking, for maximum accuracy, you would wantsinpi(2*nu ...)instead ofsin(2*pi*nu...). Check whether there issinpi()in Python or a suitable Python library. Otherwise, I would suggest filing a feature request with the Python developers to addsinpi(). $\endgroup$sinpi()function :( $\endgroup$