DDA_numercial_method_1d

import numpy as np
import math as m
import cmath as cm
import scipy.integrate as spint
import scipy
import matplotlib.pyplot as plt
import time

# dda problem for dielectric.
# finding polarization for a dieletric with incedent eliectric field
# considered cube NxNxN
N = 10 # number of segments !note, N<20!

l = 0.01 # length [m]
eps = 9 # epslion
# freq = 3e8/l / m.sqrt(eps)
# freq = 3e8/l
freq = 12.8e9
k = freq*2*m.pi/(3e8) # k-vector
a = l/N # length of each segment
betta = ( (3*a**3)*(eps - 1) / (4*m.pi *(eps + 2)) )**(-1) + 1j*k**3/(6*m.pi) #alpha^-1

# grin function for xx component
def grin(a,b,size):
    i1 = (a % N)
    j1 = (a//N) % N
    k1 = a // N**2
    i2 = (b % N)
    j2 = (b//N) % N
    k2 = b // N**2
    # print(i1,j1,k1)
    # print(i2,j2,k2)
    # print(k)
    R = size*m.sqrt( (i1-i2)**2 + (j1-j2)**2 + (k1-k2)**2)
    mul1 = ( (j2-j1)**2+(k2-k1)**2 + ( 2*(i2-i1)**2-(j2-j1)**2-(k2-k1)**2 )*( -1j/(k*R)+ (k*R)**(-2) ) )*size**2
    return k**2*cm.exp(1j*k*R)/(4*m.pi*R**3) * mul1



start_time = time.time()

A = np.zeros((N**3,N**3),dtype = complex)
for i in range(N**3):
    for j in range(N**3):
        if (i == j ):
            A[i][j] = betta
        else:
            A[i][j] = - grin(i,j,a)


E = np.ones((N,N,N),dtype=complex)
for i in range(N):
    E[:,:,i]= cm.exp(1j*k*a*i)
E_new = np.zeros(N**3,dtype  = complex)
for i in range(N**3):
    E_new[i] = E[i % N][i//N % N][i // N**2]

answer = np.linalg.solve(A,E_new)
p = np.zeros((N,N,N),dtype = complex )
for i in range(N**3):
    p[i % N][i//N % N][i // N**2] = answer[i]

p = p/a
# plt.plot(scipy.real((E[0,0,:])))
# plt.show()
end_time = time.time()
print(end_time-start_time,"sec")
plt.subplot(311)
# plt.imshow(scipy.real(p[:,:,m.ceil(N/2)-1]))
# plt.imshow(abs(p[:,:,m.ceil(N/2)-1]),aspect = 'auto')
plt.imshow(abs(p[:,:,5]),aspect = 'auto')
plt.title('Z const, Plane (XOY)')
plt.colorbar()
plt.subplot(312)
plt.imshow(abs(p[:,m.ceil(N/2),:]),aspect = 'auto')
plt.colorbar()
plt.title('Y const, Plane (XOZ)')
plt.subplot(313)
plt.imshow(abs(p[m.ceil(N/2),:,:]),aspect = 'auto')
plt.colorbar()
plt.title('X const, Plane (YOZ)')
plt.show()

plt.subplot(311)
plt.imshow(scipy.real(p[:,:,5]),aspect = 'auto')
plt.title('Z const, Plane (XOY)')
plt.colorbar()
plt.subplot(312)
plt.imshow(scipy.real(p[:,m.ceil(N/2),:]),aspect = 'auto')
plt.colorbar()
plt.title('Y const, Plane (XOZ)')
plt.subplot(313)
plt.imshow(scipy.real(p[m.ceil(N/2),:,:]),aspect = 'auto')
plt.colorbar()
plt.title('X const, Plane (YOZ)')
plt.show()

plt.subplot(311)
plt.imshow(scipy.imag(p[:,:,5]),aspect = 'auto')
plt.title('Z const, Plane (XOY)')
plt.colorbar()
plt.subplot(312)
plt.imshow(scipy.imag(p[:,m.ceil(N/2),:]),aspect = 'auto')
plt.colorbar()
plt.title('Y const, Plane (XOZ)')
plt.subplot(313)
plt.imshow(scipy.imag(p[m.ceil(N/2),:,:]),aspect = 'auto')
plt.colorbar()
plt.title('X const, Plane (YOZ)')
plt.show()
# print(m.ceil(N/2))
print(freq/1e9,"GHz")