1
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Potentially troublesome: 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 want sinpi(2*nu ...) instead of sin(2*pi*nu...). Check whether there is sinpi() in Python or a suitable Python library. Otherwise, I would suggest filing a feature request with the Python developers to add sinpi(). $\endgroup$
    – njuffa
    Mar 8 at 5:38
  • $\begingroup$ The argument cannot be that big, nu*delta t is of order of ~0.1, n goes to ~100, so the argument is reasonable I think. There is no sinpi() function :( $\endgroup$
    – Hridey
    Mar 8 at 5:45
  • $\begingroup$ Probably harmless then. I was thinking this could be the kind of case encountered in this question. I looked over the code and couldn't find anything obviously out of place like an incorrect variable name, or incorrect indexing. Is this a straight port of the MATLAB code in the book, or a heavily modified version? $\endgroup$
    – njuffa
    Mar 8 at 6:08
  • $\begingroup$ This is 90% the same version. Should I add the MATLAB code? (it is image not text :( ) $\endgroup$
    – Hridey
    Mar 8 at 6:12
  • $\begingroup$ I do not think that is a good idea as it would seem to involve the unauthorized propagation of copyrighted material. $\endgroup$
    – njuffa
    Mar 8 at 6:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.