Openems patch antenna array

[mudather2003]

Hi , I am working on project to design Smart mimo microstrip patch antenna array on python interface With little experience in programming, I now face many problems without finding any solutions to them
This is the code please help me
import sys
def keep_window_open_on_error(exc_type, exc_value, tb):
import traceback
traceback.print_exception(exc_type, exc_value, tb)
input("\n[حدث خطأ برمجي!] اضغط مفتاح Enter لإغلاق النافذة...")
sys.exit(1)

sys.excepthook = keep_window_open_on_error

import os, tempfile, csv
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

=========================================================================

إصلاح مسارات وقواعد بيانات الـ DLL (هام جداً لمستخدمي ويندوز)

=========================================================================

OPENEMS_PATH = os.environ.get('OPENEMS_PATH', r'C:\openEMS_v0.0.36\openEMS')
os.environ['OPENEMS_INSTALL_PATH'] = OPENEMS_PATH

if os.name == 'nt' and hasattr(os, 'add_dll_directory'):
if os.path.exists(OPENEMS_PATH):
os.add_dll_directory(OPENEMS_PATH)
else:
print(f"تحذير: لم يتم العثور على مسار openEMS المحدد: {OPENEMS_PATH}")
print("تأكد أن openEMS مثبت في المسار الصحيح، على سبيل المثال: D:\Refrince\openEMS")
sys.exit(1)

qt5_path = os.path.join(OPENEMS_PATH, 'qt5')
if os.path.exists(qt5_path):
    os.add_dll_directory(qt5_path)

from CSXCAD import ContinuousStructure
from openEMS import openEMS
from openEMS.physical_constants import *

=========================================================================

المعاملات والأبعاد الهندسية (PARAMETERS)

=========================================================================

Sim_Path = os.path.join(tempfile.gettempdir(), 'MIMO_Patch_2x2')

TARGET_FREQ_GHZ = 2.4
f0 = TARGET_FREQ_GHZ * 1e9
fc = 1e9

patch_width = 32
patch_length = 40

spacing_x = 60
spacing_y = 60

substrate_epsR = 3.38
substrate_thickness = 1.524
substrate_width = 160
substrate_length = 160

feed_pos = -8
feed_R = 50
SimBox = np.array([320, 320, 180])

element_positions = {
1: (-spacing_x/2, -spacing_y/2),
2: ( spacing_x/2, -spacing_y/2),
3: (-spacing_x/2, spacing_y/2),
4: ( spacing_x/2, spacing_y/2),
}
NPORTS = 4

sub_start = [-substrate_width/2, -substrate_length/2, 0]
sub_stop = [ substrate_width/2, substrate_length/2, substrate_thickness]
gnd_start = sub_start.copy()
gnd_stop = sub_stop.copy()
gnd_start[2] = 0
gnd_stop[2] = 0

إنشاء مجلدات الإخراج لحفظ النتائج والمخططات

plot_dir = os.path.join(Sim_Path, 'Plots')
data_dir = os.path.join(Sim_Path, 'Data')
pat_dir = os.path.join(Sim_Path, 'Patterns')
os.makedirs(plot_dir, exist_ok=True)
os.makedirs(data_dir, exist_ok=True)
os.makedirs(pat_dir, exist_ok=True)

=========================================================================

تشغيل المحاكاة الفردية لكل منفذ (أربعة محاكيات منفصلة)

=========================================================================

ports_run = {}
for active_port in range(1, NPORTS + 1):
print(f"\n--- إعداد وتشغيل المحاكاة للمنفذ {active_port} من {NPORTS} ---")

FDTD = openEMS(NrTS=30000, EndCriteria=1e-1)
FDTD.SetGaussExcite(f0, fc)
FDTD.SetBoundaryCond(['PML_8'] * 6)

CSX = ContinuousStructure()
FDTD.SetCSX(CSX)

mesh = CSX.GetGrid()
mesh.SetDeltaUnit(1e-3)
mesh_res = C0/(f0+fc)/1e-3/40

mesh.AddLine('x', [-SimBox[0]/2, SimBox[0]/2])
mesh.AddLine('y', [-SimBox[1]/2, SimBox[1]/2])
mesh.AddLine('z', [-SimBox[2]/3, SimBox[2]*2/3])

patch = CSX.AddMetal('patch')
for port_nr, (xc, yc) in element_positions.items():
    start = [xc - patch_width/2, yc - patch_length/2, substrate_thickness]
    stop  = [xc + patch_width/2, yc + patch_length/2, substrate_thickness]
    patch.AddBox(priority=10, start=start, stop=stop)
    mesh.AddLine('x', [start[0], stop[0], xc + feed_pos])
    mesh.AddLine('y', [start[1], stop[1], yc])

FDTD.AddEdges2Grid(dirs='xy', properties=patch, metal_edge_res=mesh_res/2)

substrate = CSX.AddMaterial('substrate', epsilon=substrate_epsR)
substrate.AddBox(priority=0, start=sub_start, stop=sub_stop)
mesh.AddLine('x', [sub_start[0], sub_stop[0]])
mesh.AddLine('y', [sub_start[1], sub_stop[1]])

gnd = CSX.AddMetal('gnd')
gnd.AddBox(gnd_start, gnd_stop, priority=10)
FDTD.AddEdges2Grid(dirs='xy', properties=gnd)

mesh.AddLine('z', np.linspace(0, substrate_thickness, 5))

ports = {}
for port_nr, (xc, yc) in element_positions.items():
    excitation = 1.0 if active_port == port_nr else 0.0
    port_start = [xc + feed_pos, yc, 0]
    port_stop  = [xc + feed_pos, yc, substrate_thickness]
    ports[port_nr] = FDTD.AddLumpedPort(port_nr, feed_R, port_start, port_stop, 'z', excitation, priority=5, edges2grid='xy')

mesh.SmoothMeshLines('all', mesh_res, 1.2)
nf2ff = FDTD.CreateNF2FFBox()

Sim_Path_case = os.path.join(Sim_Path, f'Port{active_port}')
os.makedirs(Sim_Path_case, exist_ok=True)

CSX.Write2XML(os.path.join(Sim_Path_case, 'mimo_2x2_antenna.xml'))
FDTD.Run(Sim_Path_case)
ports_run[active_port] = {'ports': ports, 'nf2ff': nf2ff, 'path': Sim_Path_case}

=========================================================================

استخراج ومعالجة معاملات S بالكامل (S-Parameters)

=========================================================================

print("\n" + "="*50 + "\nحساب معاملات S-Parameters...\n" + "="*50)
f = np.linspace(max(1e9, f0-fc), f0+fc, 401)
S = np.zeros((NPORTS, NPORTS, len(f)), dtype=complex)

for ap in range(1, NPORTS+1):
data = ports_run[ap]
ports_dict = data['ports']
sim_dir = data['path']
for pn in range(1, NPORTS+1):
ports_dict[pn].CalcPort(sim_dir, f)
inc = ports_dict[ap].uf_inc
for pn in range(1, NPORTS+1):
S[pn-1, ap-1, :] = ports_dict[pn].uf_ref / inc

S_dB = 20*np.log10(np.maximum(np.abs(S), 1e-10))
idx_f0 = np.argmin(np.abs(f - f0))

حفظ جدول قيم S-Parameters عند التردد المستهدف f0

with open(os.path.join(data_dir, 'S_parameters_at_f0.txt'), 'w') as fs:
fs.write(f'Frequency: {f[idx_f0]/1e9:.3f} GHz\n')
fs.write('Sij (dB)\n')
for i in range(NPORTS):
for j in range(NPORTS):
fs.write(f'S{i+1}{j+1} = {S_dB[i,j,idx_f0]:.2f} dB\n')

تصدير معاملات S بالكامل إلى ملف CSV

with open(os.path.join(data_dir, 'S_parameters_full.csv'), 'w', newline='') as fcsv:
writer = csv.writer(fcsv)
header = ['Freq(GHz)'] + [f'S{i+1}{j+1}(dB)' for i in range(NPORTS) for j in range(NPORTS)]
writer.writerow(header)
for k in range(len(f)):
row = [f[k]/1e9]
for i in range(NPORTS):
for j in range(NPORTS):
row.append(S_dB[i,j,k])
writer.writerow(row)

رسم بياني لمعاملات S

plt.figure(figsize=(12,6))
for i in range(NPORTS):
plt.plot(f/1e9, S_dB[i,i,:], label=f'S{i+1}{i+1}')
plt.plot(f/1e9, S_dB[1,0,:], '--', label='S21')
plt.plot(f/1e9, S_dB[0,1,:], '--', label='S12')
plt.xlabel('Frequency (GHz)')
plt.ylabel('S-Parameters (dB)')
plt.title('MIMO 2x2 S-Parameters')
plt.grid(True)
plt.legend()
plt.savefig(os.path.join(plot_dir, 'S_Parameters.png'), dpi=150)
plt.close()

=========================================================================

حساب معاملات ECC و Diversity Gain لمصفوفة الـ MIMO

=========================================================================

print("حساب ECC و Diversity Gain...")
ECC = np.zeros((NPORTS, NPORTS, len(f)))
DG = np.zeros((NPORTS, NPORTS, len(f)))

for i in range(NPORTS):
for j in range(i+1, NPORTS):
num = np.abs( np.sum( np.conj(S[:,i,:]) * S[:,j,:], axis=0) )**2
den = (1 - np.sum(np.abs(S[:,i,:])**2, axis=0)) * (1 - np.sum(np.abs(S[:,j,:])2, axis=0))
ecc = num / np.maximum(den, 1e-12)
ECC[i,j,:] = ecc
DG_lin = 10.0 * np.sqrt(np.maximum(1 - ecc2, 0))
DG[i,j,:] = 10.0 * np.log10(DG_lin + 1e-30)

رسم الـ ECC لجميع المنافذ المشتركة

plt.figure(figsize=(10,6))
for i in range(NPORTS):
for j in range(i+1, NPORTS):
plt.plot(f/1e9, ECC[i,j,:], label=f'ECC{i+1}{j+1}')
plt.xlabel('Frequency (GHz)')
plt.ylabel('ECC')
plt.title('Envelope Correlation Coefficient')
plt.grid(True)
plt.legend()
plt.savefig(os.path.join(plot_dir, 'ECC_all.png'), dpi=150)
plt.close()

حفظ قيم ECC و DG في ملف نصي مستقل

with open(os.path.join(data_dir, 'ECC_and_DG_at_f0.txt'), 'w') as fe:
fe.write(f'ECC and Diversity Gain at {f[idx_f0]/1e9:.3f} GHz\n')
for i in range(NPORTS):
for j in range(i+1, NPORTS):
fe.write(f'ECC{i+1}{j+1} = {ECC[i,j,idx_f0]:.5f}, DG{i+1}{j+1} = {DG[i,j,idx_f0]:.2f} dB\n')

=========================================================================

دوال استخراج أنماط الإشعاع (2D و 3D) والتعامل مع حسابات الحقل البعيد

=========================================================================

def calc_2d(nf2ff_obj, sim_dir, freq, theta_deg, phi_deg, center):
res = nf2ff_obj.CalcNF2FF(sim_dir, freq, theta_deg, [phi_deg], center=center)
E_norm = np.abs(res.E_norm[0]) / np.max(np.abs(res.E_norm[0]))
Dmax_lin = np.max(res.Dmax)
if Dmax_lin <= 0:
Dmax_lin = 1.0
return np.squeeze(20np.log10(E_norm) + 10np.log10(Dmax_lin)), Dmax_lin

def calc_3d_full(nf2ff_obj, sim_dir, freq, center):
th = np.arange(0, 181, 2)
ph_deg = np.linspace(0, 360, 37)
res = nf2ff_obj.CalcNF2FF(sim_dir, freq, th, ph_deg, center=center)
U = np.abs(res.E_theta[0])**2 + np.abs(res.E_phi[0])**2
U_norm = U / np.max(U)
P_rad = res.Prad[0]
Dmax_lin = res.Dmax[0]
return th, np.deg2rad(ph_deg), U_norm, P_rad, Dmax_lin

def plot_polar_2d(theta_deg, dBi_xz, dBi_yz, title, fname):
plt.figure(figsize=(8,8))
ax = plt.subplot(111, polar=True)
ax.plot(np.deg2rad(theta_deg), dBi_xz, 'b-', lw=2, label='xz-plane')
ax.plot(np.deg2rad(theta_deg), dBi_yz, 'r--', lw=2, label='yz-plane')
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)
ax.set_title(title)
plt.legend(loc='upper right')
plt.grid(True)
plt.savefig(fname, dpi=150)
plt.close()

def plot_3d_surf(th, ph, U_norm, title, fname):
TH, PH = np.meshgrid(np.deg2rad(th), ph, indexing='ij')
X = U_norm * np.sin(TH) * np.cos(PH)
Y = U_norm * np.sin(TH) * np.sin(PH)
Z = U_norm * np.cos(TH)
fig = plt.figure(figsize=(10,9))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='jet', linewidth=0, antialiased=True, alpha=0.9)
ax.set_title(title)
plt.colorbar(ax.collections[0], shrink=0.5, aspect=10, label='U / U_max')
plt.savefig(fname, dpi=150)
plt.close()

def save_pattern_csv(th, ph, U_norm, fname):
TH, PH = np.meshgrid(th, ph, indexing='ij')
with open(fname, 'w', newline='') as f:
writer = csv.writer(f)
writer.writerow(['theta_deg', 'phi_deg', 'U_norm'])
for i in range(len(th)):
for j in range(len(ph)):
writer.writerow([th[i], np.rad2deg(ph[j]), U_norm[i,j]])

=========================================================================

استخراج المخططات الإشعاعية الفردية + حساب الكفاءة والكسب الحقيقي الدقيق

=========================================================================

print("استخراج الأنماط الفردية، الاتجاهية، الكفاءة، الكسب المحقق...")
theta_2d = np.arange(0, 361, 2)
center = [0, 0, substrate_thickness/2]
Dmax_lin = {}
eff_rad = {}
eff_total = {}
G_real_dBi = {}

for ap in range(1, NPORTS+1):
nf2ff = ports_run[ap]['nf2ff']
pdir = ports_run[ap]['path']

dBi_xz, _ = calc_2d(nf2ff, pdir, f0, theta_2d, 0, center)
dBi_yz, _ = calc_2d(nf2ff, pdir, f0, theta_2d, 90, center)
plot_polar_2d(theta_2d, dBi_xz, dBi_yz,
              f'Port {ap} Individual Pattern ({f0/1e9:.2f} GHz)',
              os.path.join(plot_dir, f'Radiation_2D_Port{ap}.png'))

th3, ph3, Un, P_rad, dmax_val = calc_3d_full(nf2ff, pdir, f0, center)
Dmax_lin[ap] = dmax_val

plot_3d_surf(th3, ph3, Un,
             f'Port {ap} 3D Pattern ({f0/1e9:.2f} GHz)',
             os.path.join(plot_dir, f'Radiation_3D_Port{ap}.png'))

save_pattern_csv(th3, ph3, Un, os.path.join(pat_dir, f'Pattern_Port{ap}.csv'))

ports_dict = ports_run[ap]['ports']

# 1. حساب القدرة المغذية الداخلة الصافية من المنفذ النشط فقط
uf_inc_active = ports_dict[ap].uf_inc[idx_f0]
P_inc = 0.5 * np.abs(uf_inc_active)**2 / feed_R

# 2. حساب مجموع القدرات المرتدة والمتسربة المستهلكة في جميع المنافذ (النشطة والساكنة)
P_out_total = 0.0
for pn in range(1, NPORTS+1):
    uf_ref_pn = ports_dict[pn].uf_ref[idx_f0]
    P_out_total += 0.5 * np.abs(uf_ref_pn)**2 / feed_R

# 3. حساب القدرة الصافية المقبولة داخل الهيكل بالكامل
P_acc = P_inc - P_out_total
if P_acc <= 0:
    P_acc = 1e-12 

# 4. حساب الكفاءات المعتمدة على موازنة الطاقة الفعلية المقاسة بالوات (مصححة رقمياً لـ FDTD)
eff_rad[ap] = P_rad / P_acc
eff_total[ap] = P_rad / P_inc

if eff_rad[ap] > 1.0: eff_rad[ap] = 1.0
if eff_total[ap] > 1.0: eff_total[ap] = 1.0

# الكسب المحقق (Realized Gain) = الاتجاهية × الكفاءة الكلية
G_real_lin = Dmax_lin[ap] * eff_total[ap]
G_real_dBi[ap] = 10.0 * np.log10(G_real_lin) if G_real_lin > 0 else -100.0

print(f"  Port{ap}: P_rad={P_rad:.6g}, P_acc={P_acc:.6g}, Directivity={10*np.log10(Dmax_lin[ap]):.2f} dBi, G_realized={G_real_dBi[ap]:.2f} dBi")

=========================================================================

حساب عرض النطاق الترددي الفردي والمشترك (Bandwidth Sii < -10 dB)

=========================================================================

print("حساب عرض النطاق...")
bandwidth_ports = {}
bandwidth_common = (None, None, 0.0)

for ap in range(1, NPORTS+1):
s_dB = S_dB[ap-1, ap-1, :]
mask = s_dB < -10
if np.any(mask):
idx = np.where(mask)[0]
f_low = f[idx[0]]
f_high = f[idx[-1]]
bw = (f_high - f_low)/1e9
bandwidth_ports[ap] = (f_low/1e9, f_high/1e9, bw)
else:
bandwidth_ports[ap] = (None, None, 0.0)

mask_all = np.all(S_dB.diagonal(axis1=0, axis2=1) < -10, axis=0)
if np.any(mask_all):
idx_all = np.where(mask_all)[0]
f_low_all = f[idx_all[0]]/1e9
f_high_all = f[idx_all[-1]]/1e9
bw_common = f_high_all - f_low_all
bandwidth_common = (f_low_all, f_high_all, bw_common)

=========================================================================

نمط المصفوفة المجمع والكامل (Combined Array Pattern)

=========================================================================

print("حساب نمط المصفوفة المجمع...")
def array_2d(theta_vals, phi_val):
Et = np.zeros(len(theta_vals), dtype=complex)
Ep = np.zeros(len(theta_vals), dtype=complex)
for ap in range(1, NPORTS+1):
pdir = ports_run[ap]['path']
res = ports_run[ap]['nf2ff'].CalcNF2FF(pdir, f0, theta_vals, [phi_val], center=center)
Et += res.E_theta[0].flatten()
Ep += res.E_phi[0].flatten()
U = np.abs(Et)**2 + np.abs(Ep)**2
U_norm = U / np.max(U)
return 10*np.log10(U_norm + 1e-12)

array_xz = array_2d(theta_2d, 0)
array_yz = array_2d(theta_2d, 90)
plot_polar_2d(theta_2d, array_xz, array_yz,
f'Full Array Combined Pattern ({f0/1e9:.2f} GHz)',
os.path.join(plot_dir, 'Radiation_2D_Array_Combined.png'))

th3 = np.arange(0, 181, 2)
ph3 = np.linspace(0, 2*np.pi, 37)
U_arr = np.zeros((len(th3), len(ph3)))
for i, phi_rad in enumerate(ph3):
Etot = np.zeros(len(th3), dtype=complex)
Epot = np.zeros(len(th3), dtype=complex)
for ap in range(1, NPORTS+1):
pdir = ports_run[ap]['path']
cut = ports_run[ap]['nf2ff'].CalcNF2FF(pdir, f0, th3, [np.rad2deg(phi_rad)], center=center)
Etot += cut.E_theta[0].flatten()
Epot += cut.E_phi[0].flatten()
U_arr[:,i] = np.abs(Etot)**2 + np.abs(Epot)**2
U_arr_norm = U_arr / np.max(U_arr)
plot_3d_surf(th3, ph3, U_arr_norm,
f'Full Array Combined 3D Pattern ({f0/1e9:.2f} GHz)',
os.path.join(plot_dir, 'Radiation_3D_Array_Combined.png'))
save_pattern_csv(th3, ph3, U_arr_norm, os.path.join(pat_dir, 'Pattern_Array_Combined.csv'))

=========================================================================

حساب تشكيل الحزمة وتوجيهها (Beamforming) بين المنافذ

=========================================================================

print("حساب تشكيل الحزمة...")
def beamform(phase_shift_deg):
Et = np.zeros(len(theta_2d), dtype=complex)
Ep = np.zeros(len(theta_2d), dtype=complex)
for ap in [1, 2]:
pdir = ports_run[ap]['path']
res = ports_run[ap]['nf2ff'].CalcNF2FF(pdir, f0, theta_2d, [0], center=center)
w = 1.0 if ap == 1 else np.exp(1j * np.deg2rad(phase_shift_deg))
Et += res.E_theta[0].flatten() * w
Ep += res.E_phi[0].flatten() * w
U = np.abs(Et)**2 + np.abs(Ep)**2
U_norm = U / np.max(U)
return 10*np.log10(U_norm + 1e-12)

patterns_bf = {0: beamform(0), 45: beamform(45), 90: beamform(90)}
plt.figure(figsize=(10,6))
for ph, pat in patterns_bf.items():
plt.plot(theta_2d, pat, label=f'Phase shift = {ph}°')
plt.xlabel('Theta (deg)')
plt.ylabel('Normalized Directivity (dB)')
plt.title('Beamforming between Port1 and Port2 (xz-plane)')
plt.grid(True)
plt.legend()
plt.xlim([0, 360])
plt.ylim([-30, 5])
plt.savefig(os.path.join(plot_dir, 'Beamforming.png'), dpi=150)
plt.close()

=========================================================================

صياغة وحفظ التقرير النهائي (FINAL REPORT)

=========================================================================

print("حفظ التقرير النهائي...")
with open(os.path.join(data_dir, 'final_report.txt'), 'w', encoding='utf-8') as rep:
rep.write("="60 + "\n")
rep.write(f"نتائج هوائي MIMO 2x2 عند {f[idx_f0]/1e9:.3f} GHz\n")
rep.write("="60 + "\n\n")
rep.write(" معاملات S (dB) : **\n")
for i in range(NPORTS):
for j in range(NPORTS):
rep.write(f" S{i+1}{j+1} = {S_dB[i,j,idx_f0]:.2f} dB\n")

rep.write("\n** ECC و Diversity Gain : **\n")
for i in range(NPORTS):
    for j in range(i+1, NPORTS):
        rep.write(f"  ECC{i+1}{j+1} = {ECC[i,j,idx_f0]:.5f}, DG{i+1}{j+1} = {DG[i,j,idx_f0]:.2f} dB\n")
        
rep.write("\n** اتجاهية، كفاءة، وكسب محقق : **\n")
for ap in range(1, NPORTS+1):
    rep.write(f"  Port{ap}: Directivity = {10*np.log10(Dmax_lin[ap]):.2f} dBi, " +
              f"Radiation Eff. = {eff_rad[ap]*100:.1f}%, Total Eff. = {eff_total[ap]*100:.1f}%, " +
              f"Realized Gain = {G_real_dBi[ap]:.2f} dBi\n")
              
rep.write("\n** عرض النطاق (Sii < -10 dB) : **\n")
for ap in range(1, NPORTS+1):
    fl, fh, bw = bandwidth_ports[ap]
    if fl:
        rep.write(f"  Port{ap}: {fl:.3f} - {fh:.3f} GHz (BW = {bw:.3f} GHz)\n")
    else:
        rep.write(f"  Port{ap}: لم يتحقق الشرط\n")
        
if bandwidth_common[0]:
    rep.write(f"\n** عرض النطاق المشترك: {bandwidth_common[0]:.3f} - {bandwidth_common[1]:.3f} GHz (BW = {bandwidth_common[2]:.3f} GHz) **\n")
else:
    rep.write("\n** لا يوجد عرض نطاق مشترك يلبي الشرط **\n")

print(f"\n✅ تم الانتهاء بنجاح. جميع النتائج محفوظة في:\n{Sim_Path}")
print(f" - الرسوم البيانية: {plot_dir}")
print(f" - البيانات الرقمية والتقرير: {data_dir}")
print(f" - الأنماط الإشعاعية بصيغة (CSV): {pat_dir}")
plt.show()