Source driven problem using Method of Moments

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.
Microstrip
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):
StripImp
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()

Hi. Without further information I won’t be able to help. We are happy to to give advise on the software. But I am not familiar enough with your application to give competent advise on what you are doing there.

Okay, could you help me with computing the surface current at a certain point(s)/ or edge of my grid once I have solved the system?

For example, in this case, I am interested in obtaining the current, I_m, across the edge m. Because with this current I can solve for the impedance, using the equation I posted earlier, Z_{in} = V_{in} / l_m^2 I_m.

Hello,

What operator does BEMPP have that can evaluate the current of my grid?

For example, in the problem where I have \boldsymbol{F} = \int f \cdot u, where f is the mathematical description of my geometry, u is the complex scalar values of my solution data, and \boldsymbol{F} is the quantity I am solving for, in my case current.