Hello,
I am attempting to solve the following problem where I have a microstrip, and I excite the edge m using method of moments (MoM) . This is related to a previous post but it is a different question, which is why I am making a new thread.
I create a set of zero coefficients and place a voltage at the DOF point of my edge. Using these coefficients, I define a Gridfunction and solve the system. Where I am solving for the Electric Field Integral Equation (EFIE), represented by the maxwell single layer operator.
Once solved, I try to find the input impedance of the microstrip. Following the delta-gap excitation definition, we can relate the voltage to the edge’s basis function, the edge length, and current across the edge. Because the normal component of the RWG function is always one, then we obtain Z_{in} = V_{in} / l_m^2 I_m. Giving us the following result for the impedance.
Compared to the literature (from W.C.Gibson MoM in Electromagnetics book):

The magnitude is off, and the imaginary component is flipped. Something is wrong, but I do not know what. I need help identifying the problem.
Thank you!
Here is the code for the problem:
import bempp
import numpy as np
from bempp.api.linalg import lu
from scipy.constants.constants import speed_of_light
#Set-up problem
stepSize = 55; #Solution Step
fMax = 1000E6 #1Ghz
fMin = 100E6 #100Mhz
fStep = (fMax - fMin) / (stepSize-1)
frequency = np.arange(stepSize)*fStep + fMin
vacuum_permittivity = 8.854187817E-12
vacuum_permeability = 4 * np.pi * 1E-7
w = 2 * np.pi * frequency
k = w * np.sqrt(vacuum_permittivity * vacuum_permeability)
wavelength = speed_of_light / fMax
Vin = 1.0 #Excitation
#Rectangle Microstrip
W = 0.003 # 3mm
L = 1 #1m
strip = np.array([[-L/2, W/2, 0], [L/2, W/2, 0], [L/2, -W/2, 0], [-L/2, -W/2, 0]])
#Create empty Zin
Zin = np.zeros(stepSize, dtype = 'complex')
#Mesh
grid = bempp.api.shapes.shapes.screen(strip, h=(1/10)*wavelength)
edgeInd = 70 #Edge of interest
#Basis Functions to represent problem
div_space = bempp.api.function_space(grid, "RWG", 0)
curl_space = bempp.api.function_space(grid, "SNC", 0)
#Set-Up Coefficients
coefficients=np.zeros(div_space.global_dof_count, dtype = 'complex')
for edges, dofs, mults in zip(grid.element_edges.T, div_space.local2global, div_space.local_multipliers):
for i, e in enumerate(edges):
if e == edgeInd:
place = dofs[i]
coefficients[place] = Vin #Place excitation
#Excitation GridFunction
trace_fun = bempp.api.GridFunction(div_space, coefficients= coefficients, dual_space=curl_space)
for ite in range(stepSize):
#Set-up Appropriate Operators - EFIE
efie = bempp.api.operators.boundary.maxwell.electric_field(div_space, div_space, curl_space, k[ite])
#Solve System
solved_system = lu(efie, trace_fun)
#Impedance
I = solved_system.coefficients[place]*W*W
Zin[ite] = Vin / I
#Plot
import matplotlib
from matplotlib import pylab as plt
matplotlib.rcParams['figure.figsize'] = (10.0, 8.0)
plt.plot(frequency*1E-9, np.real(Zin), label = 'Resistance')
plt.plot(frequency*1E-9, np.imag(Zin), '--', label = 'Reactance')
plt.xlim((np.min(frequency)*1E-9,np.max(frequency)*1E-9))
plt.xlabel('Frequency')
plt.ylabel('Impedance')
plt.title('Impedance of Microstrip')
plt.grid(linewidth=1)
plt.legend()
plt.show()