Hi,
Below I attach the code and results (it can be improved) when experimenting with computing the input impedance of a patch antenna - assuming a feed point/probe model. The resistance part should follow cos^2(distanceFromEdge). The example uses the computed E and H fields to estimate the Input impedance. I also tried with the matrix impedance (Ele.weak_form()._impl) but I could not obtain reasonable values:
Full code and examples: Ph.D. Marc Jofre Cruanyes
def BemppApiGridUnion(grids, domain_indices=None, swapped_normals=None):
"""
Return the union of a given list of grids.
Parameters
----------
grids: list
A list of grid objects.
domain_indices : list
Attach a list of domain indices to the new grid such that grid[j] received the domain index domain_indices[j]
swapped_normals : list of boolean
A list of the form [False, True, ...], that specifies for each grid if the normals should be swapped (True) or not (False). This is helpful if one grid is defined to be inside another grid.
This method returns a new grid object, which is the union of the input grid objects.
"""
vertex_offset = 0
element_offset = 0
vertex_count = sum([grid.number_of_vertices for grid in grids])
element_count = sum([grid.number_of_elements for grid in grids])
vertices = np.empty((3, vertex_count), dtype="float64")
elements = np.empty((3, element_count), dtype="uint32")
all_domain_indices = np.empty(element_count, dtype="uint32")
if domain_indices is None:
domain_indices = range(len(grids))
if swapped_normals is None:
swapped_normals = len(grids) * [False]
for index, grid in enumerate(grids):
nelements = grid.number_of_elements
nvertices = grid.number_of_vertices
vertices[:, vertex_offset : vertex_offset + nvertices] = grid.vertices
if swapped_normals[index]:
current_elements = grid.elements[[0, 2, 1], :]
else:
current_elements = grid.elements
elements[:, element_offset : element_offset + nelements] = (current_elements + vertex_offset)
all_domain_indices[element_offset : element_offset + nelements] = domain_indices[index]
vertex_offset += nvertices
element_offset += nelements
return bempp.api.grid.grid.Grid(vertices, elements, all_domain_indices)
def BemppIdentifyDOFcellsFromPoint(pointV1aux,gridAux,div_spaceAux,GridSizehAux,CelldofAux=False):
## Seems better for coefficient excitation to excite a single edge, and correspondingly its DOF. Variable CelldofAux=False identifies the closests edge DOF. Variable CelldofAux=True identifies the cells DOFs
#Get edges
EdgeInd=[] # Number identifying the edge number. A list because maybe more than one point is found
v1EdgeInd=[] # List of vertice 1
v2EdgeInd=[] # List of vertice 2
Diff1EdgeInd=[] # List of distances of vertice 1
Diff2EdgeInd=[] # List of distances of vertice 2
for i in np.arange(0,gridAux.number_of_edges,1): # First identify a few candidate edges
v1=gridAux.edges[0,i] # vertice 1 of the edge
V1axis=np.abs(gridAux.vertices[:, v1]-pointV1aux)
DiffPoint1=np.sqrt(V1axis[0]**2+V1axis[1]**2+V1axis[2]**2)
v2=gridAux.edges[1,i] # vertice 2 of the edge
V2axis=np.abs(gridAux.vertices[:, v2]-pointV1aux)
DiffPoint2=np.sqrt(V2axis[0]**2+V2axis[1]**2+V2axis[2]**2)
#if (DiffPoint1<1.5*GridSizehAux and DiffPoint2<1.5*GridSizehAux):
EdgeInd.append(i)
Diff1EdgeInd.append(DiffPoint1)
Diff2EdgeInd.append(DiffPoint2)
v1EdgeInd.append(gridAux.vertices[:, v1])
v2EdgeInd.append(gridAux.vertices[:, v2])
# Order from minimum distance (the sum of the two distances between the two vertices of the edge) to maximum distance the identified edges
EdgeInd=np.array(EdgeInd)# Convert to numpy array
v1EdgeInd=np.array(v1EdgeInd)# Convert to numpy array
v2EdgeInd=np.array(v2EdgeInd)# Convert to numpy array
Diff1EdgeInd=np.array(Diff1EdgeInd)# Convert to numpy array
Diff2EdgeInd=np.array(Diff2EdgeInd)# Convert to numpy array
TotalDiffEdgeInd=Diff1EdgeInd+Diff2EdgeInd # Total distance of the edge vertices to the point
OrderIndicesEdgeInd=np.argsort(TotalDiffEdgeInd)
OrderedMinEdgeIndPot=EdgeInd[OrderIndicesEdgeInd]
OrderedMinv1EdgeIndPot=v1EdgeInd[OrderIndicesEdgeInd]
OrderedMinv2EdgeIndPot=v2EdgeInd[OrderIndicesEdgeInd]
OrderedMinDiff1EdgeIndPot=Diff1EdgeInd[OrderIndicesEdgeInd]
OrderedMinDiff2EdgeIndPot=Diff2EdgeInd[OrderIndicesEdgeInd]
# Consider the first two edges as the rulers. Now, order the other edges that share vertices
# Find the common vertices of the first two edges, and then find the opposite vertices that the third edge should have
OrderedMinEdgeInd=np.zeros_like(OrderedMinEdgeIndPot)
OrderedMinv1EdgeInd=np.zeros_like(OrderedMinv1EdgeIndPot)
OrderedMinv2EdgeInd=np.zeros_like(OrderedMinv2EdgeIndPot)
OrderedMinDiff1EdgeInd=np.zeros_like(OrderedMinDiff1EdgeIndPot)
OrderedMinDiff2EdgeInd=np.zeros_like(OrderedMinDiff2EdgeIndPot)
if (CelldofAux):
OrderedMinEdgeInd[0:2]=OrderedMinEdgeIndPot[0:2]
OrderedMinv1EdgeInd[:,0:2]=OrderedMinv1EdgeIndPot[:,0:2]
OrderedMinv2EdgeInd[:,0:2]=OrderedMinv2EdgeIndPot[:,0:2]
OrderedMinDiff1EdgeInd[0:2]=OrderedMinDiff1EdgeIndPot[0:2]
OrderedMinDiff2EdgeInd[0:2]=OrderedMinDiff2EdgeIndPot[0:2]
v1a=gridAux.edges[0,OrderedMinEdgeInd[0]] # vertice a of edge 1
v1b=gridAux.edges[1,OrderedMinEdgeInd[0]] # vertice b of edge 1
v2a=gridAux.edges[0,OrderedMinEdgeInd[1]] # vertice a of edge 2
v2b=gridAux.edges[1,OrderedMinEdgeInd[1]] # vertice b of edge 2
if (v1a==v2a):
v3a=v1b
v3b=v2b
elif(v1a==v2b):
v3a=v1b
v3b=v2a
elif(v1b==v2a):
v3a=v1a
v3b=v2b
else: # v1b==v2b
v3a=v1a
v3b=v2a
# Find the third edge
Edge3Ind=0
for iIter in range(2,len(OrderedMinEdgeInd),1):
if (((v3a==gridAux.edges[0,OrderedMinEdgeIndPot[iIter]]).all() and (v3b==gridAux.edges[1,OrderedMinEdgeIndPot[iIter]]).all() ) or ((v3b==gridAux.edges[0,OrderedMinEdgeIndPot[iIter]]).all() and (v3a==gridAux.edges[1,OrderedMinEdgeIndPot[iIter]]).all() )):
Edge3Ind=iIter
OrderedMinEdgeInd[2]=OrderedMinEdgeIndPot[Edge3Ind]
OrderedMinv1EdgeInd[2,:]=OrderedMinv1EdgeIndPot[Edge3Ind,:]
OrderedMinv2EdgeInd[2,:]=OrderedMinv2EdgeIndPot[Edge3Ind,:]
OrderedMinDiff1EdgeInd[2]=OrderedMinDiff1EdgeIndPot[Edge3Ind]
OrderedMinDiff2EdgeInd[2]=OrderedMinDiff2EdgeIndPot[Edge3Ind]
# Find a fourth edge, following the min distance criterion
if (Edge3Ind==2):
OrderedMinEdgeInd[3]=OrderedMinEdgeIndPot[3]
OrderedMinv1EdgeInd[3,:]=OrderedMinv1EdgeIndPot[3,:]
OrderedMinv2EdgeInd[3,:]=OrderedMinv2EdgeIndPot[3,:]
OrderedMinDiff1EdgeInd[3]=OrderedMinDiff1EdgeIndPot[3]
OrderedMinDiff2EdgeInd[3]=OrderedMinDiff2EdgeIndPot[3]
else:
OrderedMinEdgeInd[3]=OrderedMinEdgeIndPot[2]
OrderedMinv1EdgeInd[3,:]=OrderedMinv1EdgeIndPot[2,:]
OrderedMinv2EdgeInd[3,:]=OrderedMinv2EdgeIndPot[2,:]
OrderedMinDiff1EdgeInd[3]=OrderedMinDiff1EdgeIndPot[2]
OrderedMinDiff2EdgeInd[3]=OrderedMinDiff2EdgeIndPot[2]
#print('OrderedMinEdgeInd: '+str(OrderedMinEdgeInd))
# Should take the first three edges, that normally define a cell. In can happen that the point of inteest falls in and edge or vertice (then this situation should be handled)
# But take the first six edges, because usually at the grid border less dofs are identified
NpotentialEdges=4
EdgeIndCell=OrderedMinEdgeInd[0:NpotentialEdges]
v1IndCell=OrderedMinv1EdgeInd[0:NpotentialEdges]
v2IndCell=OrderedMinv2EdgeInd[0:NpotentialEdges]
EdgeDiff1IndCell=OrderedMinDiff1EdgeInd[0:NpotentialEdges]
EdgeDiff2IndCell=OrderedMinDiff2EdgeInd[0:NpotentialEdges]
# Checks
if (np.sum(EdgeDiff1IndCell==0.0)>0 or np.sum(EdgeDiff2IndCell==0.0)>0): # Hit vertice
print("Vertice hit. It should be dealt a side!")
DotProdArray1=np.zeros((3),dtype=np.float32)
DotProdArray2=np.zeros((3),dtype=np.float32)
AbsProdArray1=np.zeros((3),dtype=np.float32)
AbsProdArray2=np.zeros((3),dtype=np.float32)
for iIter in range(0,3,1):
DotProdArray1[iIter]=np.abs(np.dot(pointV1aux-v1IndCell[iIter],v2IndCell[iIter]-v1IndCell[iIter]))
DotProdArray2[iIter]=np.abs(np.dot(pointV1aux-v2IndCell[iIter],v2IndCell[iIter]-v1IndCell[iIter]))
AbsProdArray1[iIter]=np.sqrt((pointV1aux[0]-v1IndCell[iIter,0])**2+(pointV1aux[1]-v1IndCell[iIter,1])**2+(pointV1aux[2]-v1IndCell[iIter,2])**2)*np.sqrt((v2IndCell[iIter,0]-v1IndCell[iIter,0])**2+(v2IndCell[iIter,1]-v1IndCell[iIter,1])**2+(v2IndCell[iIter,2]-v1IndCell[iIter,2])**2)
AbsProdArray2[iIter]=np.sqrt((pointV1aux[0]-v2IndCell[iIter,0])**2+(pointV1aux[1]-v2IndCell[iIter,1])**2+(pointV1aux[2]-v2IndCell[iIter,2])**2)*np.sqrt((v2IndCell[iIter,0]-v1IndCell[iIter,0])**2+(v2IndCell[iIter,1]-v1IndCell[iIter,1])**2+(v2IndCell[iIter,2]-v1IndCell[iIter,2])**2)
# Probably now hitting an edge is well dealt
#if (np.sum(DotProdArray1==AbsProdArray1)>0 or np.sum(DotProdArray2==AbsProdArray2)>0): # Hit edge
# print("Edge hit. It should be dealt a side!")
# Get DOF
else:
EdgeIndCell=OrderedMinEdgeIndPot[0:1]
v1IndCell=OrderedMinv1EdgeIndPot[0:1]
v2IndCell=OrderedMinv2EdgeIndPot[0:1]
EdgeDiff1IndCell=OrderedMinDiff1EdgeIndPot[0:1]
EdgeDiff2IndCell=OrderedMinDiff2EdgeIndPot[0:1]
Potentialdof=[]
"""
dofAux=pointInd
for iIter in range(0,len(EdgeIndCell),1):
for edges, dofs in zip(gridAux.edges, div_spaceAux.local2global):
for e, d in zip(edges, dofs):
if e == EdgeIndCell[iIter]:
dof[iIter] = d
"""
#Set-Up Coefficients
for iIter in range(0,len(EdgeIndCell),1):
for edges, dofs in zip(gridAux.element_edges.T, div_spaceAux.local2global):#, div_spaceAux.local_multipliers):
for i, e in enumerate(edges):
if (e == EdgeIndCell[iIter]):
Potentialdof.append(dofs[i])
if (CelldofAux): # DOFs forming a cell
dofAux = ([*set(Potentialdof)])[0:3] # Remove duplicate DOFs from the list, and consider the first three dofs (to form a cell)
else: # DOF of the closest edge
dofAux = [([*set(Potentialdof)])[0]] # Remove duplicate DOFs from the list, and consider only the first dofs (of an edge)
#print('dof: '+str(dofAux))
return dofAux
#############################################################################
# Theoretical dimensioning of patch antenna
OperationalFrequency=3e9 # Hz
OperationWavelengthFreeSpace=c_0/OperationalFrequency
wavenumberk=2.0*np.pi/OperationWavelengthFreeSpace
# Typical microstrip line
# Rectangular patch antenna: https://www.emcos.com/?application-examples=modeling-of-patch-antennas
#Patch antenna dielectric
RelativePermittivityPA=2.25 # Duroid at GHz
OperationWavelengthPA=c_0/OperationalFrequency/np.sqrt(RelativePermittivityPA)
wavenumberkPA=2.0*np.pi/OperationWavelengthPA
EffectiveRelativePermittivityApprox=(RelativePermittivityPA+1.0)/2.0
Height=0.2*c_0/OperationalFrequency/np.sqrt(EffectiveRelativePermittivityApprox)
Width=c_0/(2.0*OperationalFrequency*np.sqrt(EffectiveRelativePermittivityApprox))
# The effective relative permittivity might need to be used for the solver
EffectiveRelativePermittivity=((RelativePermittivityPA+1.0)/2.0)+((RelativePermittivityPA-1.0)/2.0)*((1.0+(10.0*Height/Width))**(-0.5))
OperationWavelengthEff=c_0/OperationalFrequency/np.sqrt(EffectiveRelativePermittivity)
wavenumberkeff=2.0*np.pi/OperationWavelengthEff
DeltaL=0.412*Height*((EffectiveRelativePermittivity+0.3)/(EffectiveRelativePermittivity-0.258))*((Width/Height+0.262)/(Width/Height+0.813))
Length=c_0/(2.0*OperationalFrequency*np.sqrt(EffectiveRelativePermittivity))-2.0*DeltaL # axis of the feed-point
WidthGround=6.0*Height+Width
LengthGround=6.0*Height+Length
# Dimensions with respect dielectric wavelength
HeightWavelengthEff=Height/OperationWavelengthEff
WidthWavelengthEff=Width/OperationWavelengthEff
LengthWavelengthEff=Length/OperationWavelengthEff
WidthGroundWavelengthEff=WidthGround/OperationWavelengthEff
LengthGroundWavelengthEff=LengthGround/OperationWavelengthEff
print('HeightWavelengthEff: '+str(HeightWavelengthEff))
print('WidthWavelengthEff: '+str(WidthWavelengthEff))
print('LengthWavelengthEff: '+str(LengthWavelengthEff))
print('WidthGroundWavelengthEff: '+str(WidthGroundWavelengthEff))
print('LengthGroundWavelengthEff: '+str(LengthGroundWavelengthEff))
#########################################################
# Study impedance - edge excitation
# https://speag.swiss/assets/downloads/publications/Samaras2004.pdf
# Attention, impedance depends on the mode launched
# Some theory:
# https://www.linkedin.com/pulse/how-model-antenna-feed-point-correctly-tony-golden
# https://www.linkedin.com/pulse/what-best-method-compute-antenna-input-impedance-tony-golden
GridSizeh=0.075*np.min([OperationWavelengthFreeSpace,OperationWavelengthPA,EffectiveRelativePermittivity]) # h : float; A floating point number specifying the grid size. Very critically computation value. Dimensions and positions of interest should be larger than this value.
Xlength=4.0*OperationWavelengthFreeSpace
Ylength=Xlength
Zrange=8.0*OperationWavelengthFreeSpace
# Number of points has to go in accordance with grid size
NxPoints=int(Xlength/GridSizeh)
NyPoints=int(Ylength/GridSizeh)
NzPoints=int(Zrange/GridSizeh)
[xMatrix, yMatrix, zMatrix] = np.mgrid[-Xlength/2.0:Xlength/2.0:NxPoints*1j,-Ylength/2.0:Ylength/2.0:NyPoints*1j,-2.0*Zrange/8.0:6.0*Zrange/8.0:NzPoints*1j]#np.mgrid[-extent:extent:nx * 1j,0:0:1j,-extent:extent:nz * 1j]
xCoord=np.linspace(-Xlength/2.0,Xlength/2.0,NxPoints)
yCoord=np.linspace(-Ylength/2.0,Ylength/2.0,NyPoints)
zCoord=np.linspace(-2.0*Zrange/8.0,6.0*Zrange/8.0,NzPoints)
dx=(xCoord[1]-xCoord[0])
dy=(yCoord[1]-yCoord[0])
dz=(zCoord[1]-zCoord[0])
PointsVstack = np.vstack((xMatrix.ravel(), yMatrix.ravel(), zMatrix.ravel()))
lenPointsVstack=PointsVstack.shape[1]
###################################################
# Compute interior and exterior indices
RectangleMaskX=np.logical_and((PointsVstack[0]-0.0)<0.5*WidthGroundWavelengthEff*OperationWavelengthEff,(PointsVstack[0]-0.0)>-0.5*WidthGroundWavelengthEff*OperationWavelengthEff)
RectangleMaskY=np.logical_and((PointsVstack[1]-0.0)<0.5*LengthGroundWavelengthEff*OperationWavelengthEff,(PointsVstack[1]-0.0)>-0.5*LengthGroundWavelengthEff*OperationWavelengthEff)
RectangleMaskZ=np.logical_and((PointsVstack[2]+0.5*HeightWavelengthEff*OperationWavelengthEff)<0.5*HeightWavelengthEff*OperationWavelengthEff,(PointsVstack[2]+0.5*HeightWavelengthEff*OperationWavelengthEff)>-0.5*HeightWavelengthEff*OperationWavelengthEff)
RectangleMask=RectangleMaskX*RectangleMaskY*RectangleMaskZ
interior_indices = np.where(RectangleMask)
exterior_indices = np.where(RectangleMask==False) #~interior_indices0 & ~interior_indices1
int_points = PointsVstack[:, interior_indices]
ext_points = PointsVstack[:, exterior_indices]
##########################
# Geometry
plate1 = np.array([[0.5*WidthGroundWavelengthEff*OperationWavelengthEff, 0.5*LengthGroundWavelengthEff*OperationWavelengthEff, -HeightWavelengthEff*OperationWavelengthEff],[-0.5*WidthGroundWavelengthEff*OperationWavelengthEff, 0.5*LengthGroundWavelengthEff*OperationWavelengthEff, -HeightWavelengthEff*OperationWavelengthEff],[-0.5*WidthGroundWavelengthEff*OperationWavelengthEff, -0.5*LengthGroundWavelengthEff*OperationWavelengthEff, -HeightWavelengthEff*OperationWavelengthEff],[0.5*WidthGroundWavelengthEff*OperationWavelengthEff, -0.5*LengthGroundWavelengthEff*OperationWavelengthEff, -HeightWavelengthEff*OperationWavelengthEff]])
plate2 = np.array([[0.5*WidthWavelengthEff*OperationWavelengthEff, 0.5*LengthWavelengthEff*OperationWavelengthEff, 0.0*OperationWavelengthEff],[-0.5*WidthWavelengthEff*OperationWavelengthEff, 0.5*LengthWavelengthEff*OperationWavelengthEff, 0.0*OperationWavelengthEff],[-0.5*WidthWavelengthEff*OperationWavelengthEff, -0.5*LengthWavelengthEff*OperationWavelengthEff, 0.0*OperationWavelengthEff],[0.5*WidthWavelengthEff*OperationWavelengthEff, -0.5*LengthWavelengthEff*OperationWavelengthEff, 0.0*OperationWavelengthEff]])
# Do not define a grid for th edielectric because it overlaps with the plates (functionBemppApiGridUnion) gridDielectricRect = bempp.api.shapes.cuboid
grid1 = bempp.api.shapes.screen(plate1,h=GridSizeh)
grid2 = bempp.api.shapes.screen(plate2,h=GridSizeh)
# BemppApiGridUnion - Not good that elements touch each other
grid = BemppApiGridUnion([grid1,grid2])
grid.plot()
# We define the spaces of order 0 RWG div-conforming functions and order 0 scaled Nédélec curl-conforming functions.
div_space = bempp.api.function_space(grid, "RWG", 0)
curl_space = bempp.api.function_space(grid, "SNC", 0)
### Source driven
VphasorIn=1.0+1j*0.0
NImpCalc=6
ZIn0=np.zeros((NImpCalc-1),dtype=np.complex64)
ZIn1=np.zeros((NImpCalc-1),dtype=np.complex64)
ZIn2=np.zeros((NImpCalc-1),dtype=np.complex64)
for iNImpCalc in range(1,NImpCalc,1): # Do not hit the center
print('iNImpCalc: '+str(iNImpCalc))
# DOF: Degree of freedom; are the points in the middle of the grid triangular-mesh-edges which the excitation coefficients are set.
# In particular, might be of interst to set coefficient for the three DOFs in the grid cells of a RWC/SNC space
# Point to identify the closest edge. The closest edge is the one that has the two vertices closest to the point of interest
pointV1 = [0.0*OperationWavelengthEff, iNImpCalc/float(NImpCalc)*LengthWavelengthEff*0.5*OperationWavelengthEff, 0.0*OperationWavelengthEff]
dof1=BemppIdentifyDOFcellsFromPoint(pointV1,grid,div_space,GridSizeh,False)
#Set-Up Coefficients
coefficients=np.zeros((div_space.global_dof_count),dtype=np.complex64)
coefficients[dof1] = VphasorIn # Source phasor. I believe it is a voltage source since the operators are electric field operators
################################################
A0_int = bempp.api.operators.boundary.maxwell.electric_field(div_space, div_space, curl_space, wavenumberkeff)
#identity = bempp.api.operators.boundary.sparse.identity(div_space, div_space, curl_space)
# Operators might be linearly combined - For instance one can dominate over the other
#A = bempp.api.GeneralizedBlockedOperator([[A0_int,identity],[A0_ext,identity]])
A = A0_int
trace_fun_int = bempp.api.GridFunction(div_space, coefficients=coefficients, dual_space=curl_space)
#trace_fun_ext = bempp.api.GridFunction(div_space, coefficients=0.0*coefficients, dual_space=curl_space)
rhs = trace_fun_int# [trace_fun_int,trace_fun_ext]#
################################
# Solver
lambda_data_sol = bempp.api.linalg.lu(A, rhs)
# Compute using electric and magnetic fields
#Volume space - cylinder//rectangle
[XCOORD,YCOORD,ZCOORD]=np.meshgrid(xCoord,yCoord,zCoord,indexing='ij')
ProximityFeedPoint=0.5*np.min([OperationWavelengthFreeSpace,OperationWavelengthPA,EffectiveRelativePermittivity])
zPlate2MaskX=np.logical_and((XCOORD-pointV1[0])<0.5*ProximityFeedPoint,(XCOORD-pointV1[0])>-0.5*ProximityFeedPoint)
zPlate2MaskY=np.logical_and((YCOORD-pointV1[1])<0.5*ProximityFeedPoint,(YCOORD-pointV1[1])>-0.5*ProximityFeedPoint)
zPlate2MaskZ=np.logical_and((ZCOORD-pointV1[2])<0.5*GridSizeh,(ZCOORD+HeightWavelengthEff*OperationWavelengthEff)>-0.5*GridSizeh)
zPlate2Mask=zPlate2MaskX*zPlate2MaskY*zPlate2MaskZ
pointFeedX=[]
pointFeedY=[]
pointFeedZ=[]
phiAngle=[]
for iIterX in range(0,NxPoints,1):
for iIterY in range(0,NyPoints,1):
for iIterZ in range(0,NzPoints,1):
if (zPlate2Mask[iIterX,iIterY,iIterZ]): # Point belongs to the volume
pointFeedX.append(xCoord[iIterX])
pointFeedY.append(yCoord[iIterY])
pointFeedZ.append(zCoord[iIterZ])
phiAngle.append(np.arctan2(yCoord[iIterY]-pointV1[1],xCoord[iIterX]-pointV1[0]))
SinglePointVstackFeed=np.vstack((pointFeedX, pointFeedY, pointFeedZ))
slp_pot_elecFeed = bempp.api.operators.potential.maxwell.electric_field(div_space, SinglePointVstackFeed, wavenumberkeff)
slp_pot_magFeed = bempp.api.operators.potential.maxwell.magnetic_field(div_space, SinglePointVstackFeed, wavenumberkeff)
EphasorFeed=-slp_pot_elecFeed*lambda_data_sol
HphasorFeed=slp_pot_magFeed*lambda_data_sol
# Assumin vertical feed. E field in z direction and H field in xy-plane. Integration along the whole plate separation
ZIn1[iNImpCalc-1]=np.nansum((EphasorFeed[2])/(-HphasorFeed[0]*np.sin(phiAngle)+HphasorFeed[1]*np.cos(phiAngle)))/float(len(phiAngle))
# Assumin vertical feed. E field in rho direction and H field in z-axis. Integration along the whole plate separation
#ZIn1[iNImpCalc-1]=np.nansum((EphasorFeed[0]*np.cos(phiAngle)+EphasorFeed[1]*np.sin(phiAngle))/(HphasorFeed[2]))/float(len(phiAngle))
## Seems like Admitance is computed. Invert with respect 1 and normalize with respect the integrated line and Relative Permittivity
ZIn1[iNImpCalc-1]=(1.0-ZIn1[iNImpCalc-1])/(HeightWavelengthEff*OperationWavelengthEff)*np.sqrt(RelativePermittivityPA)
print('Impedance1: '+str(ZIn1[iNImpCalc-1]))
"""
############################################
# Compute using the impedance matrix
#Grid space
[XCOORD,YCOORD,ZCOORD]=np.meshgrid(xCoord,yCoord,zCoord,indexing='ij')
ProximityFeedPoint=3.0*np.min([OperationWavelengthFreeSpace,OperationWavelengthPA,EffectiveRelativePermittivity])
zPlate2MaskX=np.logical_and((XCOORD-pointV1[0])<0.5*ProximityFeedPoint,(XCOORD-pointV1[0])>-0.5*ProximityFeedPoint)
zPlate2MaskY=np.logical_and((YCOORD-pointV1[1])<0.5*ProximityFeedPoint,(YCOORD-pointV1[1])>-0.5*ProximityFeedPoint)
zPlate2MaskZ=np.logical_and((ZCOORD-pointV1[2])<0.5*GridSizeh,(ZCOORD-pointV1[2])>-0.5*GridSizeh)
zPlate2Mask=zPlate2MaskX*zPlate2MaskY*zPlate2MaskZ
zPlate1MaskX=np.logical_and((XCOORD-pointV1[0])<0.5*ProximityFeedPoint,(XCOORD-pointV1[0])>-0.5*ProximityFeedPoint)
zPlate1MaskY=np.logical_and((YCOORD-pointV1[1])<0.5*ProximityFeedPoint,(YCOORD-pointV1[1])>-0.5*ProximityFeedPoint)
zPlate1MaskZ=np.logical_and((ZCOORD+HeightWavelengthEff*OperationWavelengthEff)<0.5*GridSizeh,(ZCOORD+HeightWavelengthEff*OperationWavelengthEff)>-0.5*GridSizeh)
zPlate1Mask=zPlate1MaskX*zPlate1MaskY*zPlate1MaskZ
elec_eff = bempp.api.operators.boundary.maxwell.electric_field(div_space, div_space, curl_space, wavenumberkeff)
zAux = elec_eff.weak_form()
zCalcPhasor1=0.0+1j*0.0
zCalcPhasor2=0.0+1j*0.0
zCalcPhasor1Count=0.0
zCalcPhasor2Count=0.0
for iIterX in range(0,NxPoints,1):
for iIterY in range(0,NyPoints,1):
for iIterZ in range(0,NzPoints,1):
if (zPlate1Mask[iIterX,iIterY,iIterZ]): # Point belongs to the grid
pointCalc = [xCoord[iIterX],yCoord[iIterY],zCoord[iIterZ]]
dofCalc=BemppIdentifyDOFcellsFromPoint(pointCalc,grid,div_space,GridSizeh,False)
zCalcPhasor1+=zAux._impl[dofCalc,dofCalc]
zCalcPhasor1Count+=1.0
if (zPlate2Mask[iIterX,iIterY,iIterZ]): # Point belongs to the grid
pointCalc = [xCoord[iIterX],yCoord[iIterY],zCoord[iIterZ]]
dofCalc=BemppIdentifyDOFcellsFromPoint(pointCalc,grid,div_space,GridSizeh,False)
zCalcPhasor2+=zAux._impl[dofCalc,dofCalc]
zCalcPhasor2Count+=1.0
if (zCalcPhasor1Count>0.0):
zCalcPhasor1=zCalcPhasor1/zCalcPhasor1Count
print('zCalcPhasor1: '+str(zCalcPhasor1))
if (zCalcPhasor2Count>0.0):
zCalcPhasor2=zCalcPhasor2/zCalcPhasor2Count
print('zCalcPhasor2: '+str(zCalcPhasor2))
ZIn1[iNImpCalc-1]=zCalcPhasor1+zCalcPhasor2
print('Impedance1: '+str(ZIn1[iNImpCalc-1]))
"""
print('Impedance1 array: '+str(ZIn1))
plt.figure()
ax = plt.gca()
line1,=plt.plot(np.arange(1,NImpCalc,1)/float(NImpCalc)*LengthWavelengthEff*0.5*OperationWavelengthEff,np.real(ZIn1),'royalblue',linestyle='solid',marker='x')
plt.xlabel('Position',fontsize=10)
plt.ylabel('Real Impedance',fontsize=10)
#plt.xscale('log')
#plt.yscale('log')
#plt.xticks(np.arange(lModeNumScan[0], lModeNumScan[-1]+1, 1))
#plt.setp(ax.get_xticklabels(),fontsize = 10)
#plt.setp(ax.get_yticklabels(),fontsize = 10)
#plt.xlim([0.0,1000])
#plt.ylim([-30,25])#np.min(GeneratedSignal)*1e1,np.max([np.max(IlluminationSignalFreqPump1MetallicState),np.max(IlluminationSignalFreqPump1SeCoState)])])
#plt.legend((line1,line2,line3), ('1','2','3'),loc='best', shadow = False, fancybox = False, frameon = False,fontsize=10)# 'best'
#plt.title('Channel width: '+str(DiameterChannel*1e6)+' [um], Channel height: '+str(DiameterChannel*1e6)+' [um]')
plt.show()
plt.figure()
ax = plt.gca()
line1,=plt.plot(np.arange(1,NImpCalc,1)/float(NImpCalc)*LengthWavelengthEff*0.5*OperationWavelengthEff,np.imag(ZIn1),'royalblue',linestyle='solid',marker='x')
plt.xlabel('Position',fontsize=10)
plt.ylabel('Imag Impedance',fontsize=10)
#plt.xscale('log')
#plt.yscale('log')
#plt.xticks(np.arange(lModeNumScan[0], lModeNumScan[-1]+1, 1))
#plt.setp(ax.get_xticklabels(),fontsize = 10)
#plt.setp(ax.get_yticklabels(),fontsize = 10)
#plt.xlim([0.0,1000])
#plt.ylim([-30,25])#np.min(GeneratedSignal)*1e1,np.max([np.max(IlluminationSignalFreqPump1MetallicState),np.max(IlluminationSignalFreqPump1SeCoState)])])
#plt.legend((line1,line2,line3), ('1','2','3','4','5'),loc='best', shadow = False, fancybox = False, frameon = False,fontsize=10)# 'best'
#plt.title('Channel width: '+str(DiameterChannel*1e6)+' [um], Channel height: '+str(DiameterChannel*1e6)+' [um]')
plt.show()
plt.figure()
ax = plt.gca()
line1,=plt.plot(np.arange(1,NImpCalc,1)/float(NImpCalc)*LengthWavelengthEff*0.5*OperationWavelengthEff,np.abs(ZIn1),'royalblue',linestyle='solid',marker='x')
plt.xlabel('Position',fontsize=10)
plt.ylabel('Mag. Impedance',fontsize=10)
#plt.xscale('log')
#plt.yscale('log')
#plt.xticks(np.arange(lModeNumScan[0], lModeNumScan[-1]+1, 1))
#plt.setp(ax.get_xticklabels(),fontsize = 10)
#plt.setp(ax.get_yticklabels(),fontsize = 10)
#plt.xlim([0.0,1000])
#plt.ylim([-30,25])#np.min(GeneratedSignal)*1e1,np.max([np.max(IlluminationSignalFreqPump1MetallicState),np.max(IlluminationSignalFreqPump1SeCoState)])])
#plt.legend((line1,line2,line3), ('1','2','3'),loc='best', shadow = False, fancybox = False, frameon = False,fontsize=10)# 'best'
#plt.title('Channel width: '+str(DiameterChannel*1e6)+' [um], Channel height: '+str(DiameterChannel*1e6)+' [um]')
plt.show()
plt.figure()
ax = plt.gca()
line1,=plt.plot(np.arange(1,NImpCalc,1)/float(NImpCalc)*LengthWavelengthEff*0.5*OperationWavelengthEff,np.angle(ZIn1),'royalblue',linestyle='solid',marker='x')
plt.xlabel('Position',fontsize=10)
plt.ylabel('Phase Impedance',fontsize=10)
#plt.xscale('log')
#plt.yscale('log')
#plt.xticks(np.arange(lModeNumScan[0], lModeNumScan[-1]+1, 1))
#plt.setp(ax.get_xticklabels(),fontsize = 10)
#plt.setp(ax.get_yticklabels(),fontsize = 10)
#plt.xlim([0.0,1000])
#plt.ylim([-30,25])#np.min(GeneratedSignal)*1e1,np.max([np.max(IlluminationSignalFreqPump1MetallicState),np.max(IlluminationSignalFreqPump1SeCoState)])])
#plt.legend((line1,line2,line3), ('1','2','3','4','5'),loc='best', shadow = False, fancybox = False, frameon = False,fontsize=10)# 'best'
#plt.title('Channel width: '+str(DiameterChannel*1e6)+' [um], Channel height: '+str(DiameterChannel*1e6)+' [um]')
plt.show()
