MoM_thinWire

import numpy as np
import math as m
import cmath as cm
import scipy.integrate as spint
import scipy
import matplotlib.pyplot as plt
#Parameter
lamda=300
L=150
a=1
eps=1
mu=1
#number of segments
N = 15
#Computer constants
k0=2*m.pi/lamda
k=k0*m.sqrt(mu*eps)
step_z=float(L/(N))
#functions
#integration
def fun(value):
    R = m.sqrt( value**2 + a**2)
    return -cm.exp(-1j*k*R)/(4*R*m.pi)

def complex_quadrature(func, a, b, **kwargs):
    def real_func(x):
        return scipy.real(func(x))
    def imag_func(x):
        return scipy.imag(func(x))
    real_integral = spint.quad(real_func, a, b, **kwargs)
    imag_integral = spint.quad(imag_func, a, b, **kwargs)
    return (real_integral[0] + 1j*imag_integral[0])

#start field
V = np.ones((N,1))
#Build Impedance Matrix
Z_tmp = np.zeros((N,N),dtype=complex)
Z = np.zeros((N,N),dtype=complex)

for i in range(N):
    Z_tmp[i][i] = 1/(4*m.pi)*m.log( (m.sqrt(1+(2*a/step_z)**2)+1)/ (m.sqrt(1+(2*a/step_z)**2)-1))-1j*k*step_z/(4*m.pi)
for i in range(N):
    for j in range(N):
        r1= cm.sqrt(((i-j)*step_z+step_z/2)**2+a**2)
        r2= cm.sqrt(((i-j)*step_z-step_z/2)**2+a**2)
        t1= ((i-j)*step_z+step_z/2)*(1+1j*k*r1)/r1**3 * cm.exp(-1j*k*r1)
        t2= ((i-j)*step_z-step_z/2)*(1+1j*k*r2)/r2**3 * cm.exp(-1j*k*r2)
        if i != j :
            Z_tmp[i][j] = complex_quadrature(fun,(i-j)*step_z+step_z/2,(i-j)*step_z-step_z/2)
        Z[i][j] = k**2*Z_tmp[i][j] + t2 - t1
Z = 1j*step_z/k * Z
I = np.dot(Z,V)
np.savetxt('out.txt',Z,fmt='%.2e')
print(np.size(I))
plt.plot(np.linspace(0,L,N),abs(scipy.real(I)))
plt.title('Distribution of current on thin wire. Method of Moments')
plt.ylabel('Amplitude I, A')
plt.xlabel('L, mm')
plt.show()