From eb9be337b181d9da997cd9c50f07e5ade7d8d7f3 Mon Sep 17 00:00:00 2001 From: Jason <83615043+JJassonn69@users.noreply.github.com> Date: Sun, 3 May 2026 16:18:42 +0545 Subject: [PATCH] =?UTF-8?q?sim(antenna):=20add=204x4=20probe-fed=20array?= =?UTF-8?q?=20model=20=E2=80=94=20mutual=20coupling=20characterisation?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Single-port-driven, 15-port-terminated FDTD sim built on the same v3 patch design point (W=7.854, L=6.56, FEED_OFFSET=2.14) at the Gerber pitch (14.27 mm X / 15.01 mm Y). Reads out S_jd for all elements when port (DRIVEN_X, DRIVEN_Y) is excited. Configurable array size + driven port via env vars; default = 4x4, drive inner element (1,1). Result at PROFILE=balanced, drive (1,1): Active S11 @ 10.5 GHz : -16.7 dB, Z = 45.7 + j13.5 Ω Active -10 dB BW : 190 MHz (10.44-10.63 GHz; ≈ single-element 180) Nearest-neighbour coupling: -x edge neighbour (0,1): -18.3 dB (strongest — edge element loads less) -y/+y in-array : -22 dB +x interior neighbour : -24.7 dB Diagonal : -32 dB Far corners : -45 to -53 dB Outputs: coupling_grid.png heatmap + S_matrix.csv (all S_jd full sweep) + S11_data.csv (driven-port only). Known limitation: port box must align with mesh lines; SmoothMeshLines may sub-cell-shift seed lines, causing some DRIVEN_X/DRIVEN_Y choices to give "Unused primitive" warning + zero excitation. Default (1,1) is verified; other choices are best-effort. Documented inline. --- .../Antenna/probe_fed_array_aeris10_v3.py | 358 ++++++++++++++++++ 1 file changed, 358 insertions(+) create mode 100644 5_Simulations/Antenna/probe_fed_array_aeris10_v3.py diff --git a/5_Simulations/Antenna/probe_fed_array_aeris10_v3.py b/5_Simulations/Antenna/probe_fed_array_aeris10_v3.py new file mode 100644 index 0000000..4478721 --- /dev/null +++ b/5_Simulations/Antenna/probe_fed_array_aeris10_v3.py @@ -0,0 +1,358 @@ +#!/usr/bin/env python3 +# probe_fed_array_aeris10_v3.py +# +# 4x4 (default; configurable) probe-fed patch array sim for AERIS-10. Built +# on the same single-element design point as probe_fed_aeris10_v3.py but +# placed on the 8x16 Gerber pitch (14.27 mm X / 15.01 mm Y). +# +# Purpose: characterise mutual coupling between elements. Each patch has its +# own probe-via port; only one port is excited per sim run, the other 15 are +# terminated in 50 Ω. From this we read: +# - S_dd (active S11 of the driven element with array loaded) +# - S_jd for all other ports j (coupling driven → j) +# Pattern of |S_jd| dB values across the array tells us nearest-neighbour vs +# diagonal vs skip-one coupling, edge vs interior asymmetry, etc. +# +# Per-element design (matches probe_fed_aeris10_v3.py iter#3): +# PATCH_W = 7.854 mm PATCH_L = 6.56 mm FEED_OFFSET = 2.14 mm +# Substrate: 0.508 mm RO4350B (εr=3.48, tanδ=0.0037) +# Pitch: 14.27 mm × 15.01 mm (X-pitch ~λ₀/2 at 10.5 GHz, Y-pitch ~1.05·λ₀/2) +# +# Run: +# cd /tmp && DYLD_LIBRARY_PATH=/Users/ganeshpanth/opt/openEMS/lib \ +# PROFILE=sanity \ +# /Users/ganeshpanth/radar_venv/bin/python \ +# /Users/ganeshpanth/PLFM_RADAR/5_Simulations/Antenna/probe_fed_array_aeris10_v3.py +# +# Env overrides: +# ARRAY_NX ARRAY_NY (default 4, 4) +# PITCH_X_MM PITCH_Y_MM (default 14.27, 15.01 from Gerber) +# DRIVEN_X DRIVEN_Y (0-indexed; default = inner element 1,1) +# PATCH_W_MM PATCH_L_MM FEED_OFFSET_MM (default v3 design point) +# +# Output (in /tmp/aeris10_array_v3/): +# S_matrix.csv — driven-column S parameters (mag dB + phase deg) at 10.5 GHz +# S11_data.csv — driven port full sweep +# coupling_grid.png — heatmap of |S_jd| dB at 10.5 GHz across array + +import os +import sys +import time +import csv +import numpy as np +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +from openEMS import openEMS +from openEMS.physical_constants import C0 +from CSXCAD import ContinuousStructure +from CSXCAD.SmoothMeshLines import SmoothMeshLines + +# ============================================================================ +# PROFILES +# ============================================================================ +PROFILE = os.environ.get("PROFILE", "sanity") +profiles = { + "sanity": {"mesh_lambda_div": 18, "n_timesteps": 50000, "end_dB": -30}, + "balanced": {"mesh_lambda_div": 25, "n_timesteps": 80000, "end_dB": -40}, +} +cfg = profiles[PROFILE] + +# ============================================================================ +# BAND +# ============================================================================ +F0 = 10.5e9 +F_SPAN = 4.0e9 +F_START = F0 - F_SPAN/2 +F_STOP = F0 + F_SPAN/2 + +# ============================================================================ +# STACKUP +# ============================================================================ +T_CU = 0.035 +H_PATCH_SUB = 0.508 +EPS_RO4350B = 3.48 +TAN_RO4350B = 0.0037 + +Z_GND = 0.0 +Z_PATCH = Z_GND + T_CU + H_PATCH_SUB +Z_TOP = Z_PATCH + T_CU + +# ============================================================================ +# PATCH (per-element, from v3 iter#3) +# ============================================================================ +PATCH_W = float(os.environ.get("PATCH_W_MM", "7.854")) +PATCH_L = float(os.environ.get("PATCH_L_MM", "6.56")) +FEED_OFFSET_MM = float(os.environ.get("FEED_OFFSET_MM", "2.14")) + +# ============================================================================ +# ARRAY +# ============================================================================ +N_X = int(os.environ.get("ARRAY_NX", "4")) +N_Y = int(os.environ.get("ARRAY_NY", "4")) +PITCH_X = float(os.environ.get("PITCH_X_MM", "14.27")) +PITCH_Y = float(os.environ.get("PITCH_Y_MM", "15.01")) +DRIVEN_X = int(os.environ.get("DRIVEN_X", str(N_X // 2 - (N_X+1) % 2))) # inner element +DRIVEN_Y = int(os.environ.get("DRIVEN_Y", str(N_Y // 2 - (N_Y+1) % 2))) + +# Array footprint extent (centre patch on origin) +ARRAY_X_HALF = (N_X-1)/2 * PITCH_X + PATCH_W/2 +ARRAY_Y_HALF = (N_Y-1)/2 * PITCH_Y + PATCH_L/2 + +# Substrate / ground extents (~λ/2 margin around array) +GND_MARGIN = 14.3 +GND_X_HALF = ARRAY_X_HALF + GND_MARGIN +GND_Y_HALF = ARRAY_Y_HALF + GND_MARGIN + +# Air box +AIR_ABOVE = 14.3 +AIR_BELOW = 14.3 +AIR_X_HALF = GND_X_HALF + 8.0 +AIR_Y_HALF = GND_Y_HALF + 8.0 + +OUT_DIR = "/tmp/aeris10_array_v3" +os.makedirs(OUT_DIR, exist_ok=True) + + +def patch_center(i, j): + """Centre coordinate of patch at array index (i,j), origin at array centre.""" + x = -(N_X-1)/2 * PITCH_X + i * PITCH_X + y = -(N_Y-1)/2 * PITCH_Y + j * PITCH_Y + return x, y + + +# ============================================================================ +# Build + run +# ============================================================================ +def run_case(sim_path, profile_cfg): + fdtd = openEMS(NrTS=profile_cfg["n_timesteps"], + EndCriteria=10**(profile_cfg["end_dB"]/20.0)) + fdtd.SetGaussExcite(F0, F_SPAN/2.0) + fdtd.SetBoundaryCond(["MUR"]*6) + + CSX = ContinuousStructure() + fdtd.SetCSX(CSX) + mesh = CSX.GetGrid() + mesh.SetDeltaUnit(1e-3) + + # ---- materials ---- + eps0 = 8.854e-12 + patch_sub = CSX.AddMaterial("RO4350B", + epsilon=EPS_RO4350B, + kappa=2*np.pi*F0*EPS_RO4350B*eps0*TAN_RO4350B) + copper = CSX.AddMetal("Copper") + + # ---- substrate (full board extent) ---- + patch_sub.AddBox([-GND_X_HALF, -GND_Y_HALF, Z_GND + T_CU], + [+GND_X_HALF, +GND_Y_HALF, Z_PATCH], priority=1) + + # ---- L2: full ground plane ---- + copper.AddBox([-GND_X_HALF, -GND_Y_HALF, Z_GND], + [+GND_X_HALF, +GND_Y_HALF, Z_GND + T_CU], priority=10) + + # ---- L1: 4x4 patch array + ports ---- + ports = [] + feed_locs = [] + for i in range(N_X): + for j in range(N_Y): + cx, cy = patch_center(i, j) + copper.AddBox([cx - PATCH_W/2, cy - PATCH_L/2, Z_PATCH], + [cx + PATCH_W/2, cy + PATCH_L/2, Z_PATCH + T_CU], + priority=10) + feed_x = cx + feed_y = cy - PATCH_L/2 + FEED_OFFSET_MM + feed_locs.append((feed_x, feed_y, i, j)) + + # ---- mesh ---- + lambda_min_mm = (C0 / F_STOP) * 1000.0 + res = lambda_min_mm / profile_cfg["mesh_lambda_div"] + + # X mesh: array extent + air, plus patch edges + feed locations + xlines = [-AIR_X_HALF, -GND_X_HALF, +GND_X_HALF, +AIR_X_HALF] + for i in range(N_X): + cx, _ = patch_center(i, 0) + xlines += [cx - PATCH_W/2, cx, cx + PATCH_W/2] + # Y mesh + ylines = [-AIR_Y_HALF, -GND_Y_HALF, +GND_Y_HALF, +AIR_Y_HALF] + for j in range(N_Y): + _, cy = patch_center(0, j) + ylines += [cy - PATCH_L/2, cy, cy + PATCH_L/2, + cy - PATCH_L/2 + FEED_OFFSET_MM] # feed y location + + # Z mesh: 6 cells in substrate + air_below = list(np.arange(Z_GND - T_CU - AIR_BELOW, Z_GND - T_CU, res)) + air_above = list(np.arange(Z_TOP + res, Z_TOP + AIR_ABOVE + res, res)) + sub_interior = list(np.linspace(Z_GND + T_CU, Z_PATCH, 7)[1:-1]) + zlines = sorted(set(air_below + [ + Z_GND - T_CU, Z_GND, Z_GND + T_CU, + Z_PATCH, Z_PATCH + T_CU, + ] + sub_interior + air_above)) + + xlines = SmoothMeshLines(np.array(xlines), res) + ylines = SmoothMeshLines(np.array(ylines), res) + zlines = np.array(zlines) + mesh.AddLine("x", xlines) + mesh.AddLine("y", ylines) + mesh.AddLine("z", zlines) + n_cells = len(xlines) * len(ylines) * len(zlines) + + # ---- ports (one excited, 15 terminated 50Ω) ---- + # NOTE: each port box must land exactly on existing mesh lines. The seed + # mesh includes feed_x (= patch centre cx) and feed_y for each (i,j) — so + # this normally works — but SmoothMeshLines can sub-cell-shift seed lines + # in some configurations and a port box ends up between two mesh lines, + # leaving openEMS without an excitation cell ("Unused primitive" warning, + # zero energy). The DRIVEN=(1,1) inner-element case has been verified to + # land on the mesh; other driven-port choices are best-effort. If you see + # NaN/zero results for a different DRIVEN_X/DRIVEN_Y, that's the cause. + for (feed_x, feed_y, i, j) in feed_locs: + port_num = i * N_Y + j + 1 + excite_amp = 1.0 if (i == DRIVEN_X and j == DRIVEN_Y) else 0.0 + port = fdtd.AddLumpedPort(port_num, 50, + [feed_x, feed_y, Z_GND + T_CU], + [feed_x, feed_y, Z_PATCH], + 'z', excite=excite_amp, priority=5) + ports.append(((i, j), port)) + + # ---- run ---- + print(f"[case] {N_X}x{N_Y} array, driven=({DRIVEN_X},{DRIVEN_Y}), " + f"cells={n_cells:,}, sub={H_PATCH_SUB}mm, pitch={PITCH_X}x{PITCH_Y}mm") + t0 = time.time() + fdtd.Run(sim_path, verbose=0, cleanup=True) + dt = time.time() - t0 + + # ---- post-process ---- + freq = np.linspace(F_START, F_STOP, 401) + driven_port = next(p for (idx, p) in ports if idx == (DRIVEN_X, DRIVEN_Y)) + for (idx, p) in ports: + p.CalcPort(sim_path, freq) + + # S_jd (j is each port, d is driven). For driven port, S_dd = uf_ref/uf_inc. + # For other ports, S_jd = uf_ref_j / uf_inc_d (no incident wave at j from + # its own source since port j is unexcited). + S = {} + for (idx, p) in ports: + if idx == (DRIVEN_X, DRIVEN_Y): + S[idx] = p.uf_ref / p.uf_inc + else: + S[idx] = p.uf_ref / driven_port.uf_inc + + return freq, S, dt, ports + + +# ============================================================================ +# MAIN +# ============================================================================ +sim_path = os.path.join(OUT_DIR, "single") +freq, S, dt, ports = run_case(sim_path, cfg) + +# At 10.5 GHz +i_op = int(np.argmin(np.abs(freq - F0))) + +# Print coupling grid +print() +print("=" * 70) +print(f" 4x4 probe-fed array — driven port at ({DRIVEN_X},{DRIVEN_Y})") +print(f" Substrate: {H_PATCH_SUB} mm RO4350B, pitch {PITCH_X}x{PITCH_Y} mm") +print(f" Sim time: {dt:.1f} s") +print("=" * 70) +print() +print(f" S parameters at {F0/1e9:.2f} GHz (|S_j,driven| in dB):") +print() +# Layout grid as visual array (i is x-direction, j is y-direction) +# Print y high to low so it matches usual visual orientation +header = " " + "".join(f" i={i:1d} " for i in range(N_X)) +print(header) +for j in reversed(range(N_Y)): + row = f" j={j:1d}: " + for i in range(N_X): + val = abs(S[(i, j)][i_op]) + dB = 20*np.log10(val + 1e-30) + row += f"{dB:>7.1f}" + print(row) +print() +# Driven port S11 vs frequency +S_dd = S[(DRIVEN_X, DRIVEN_Y)] +S_dd_dB = 20*np.log10(np.abs(S_dd) + 1e-30) +zin_d = 50.0 * (1 + S_dd) / (1 - S_dd) # Z = Z0·(1+S)/(1-S) +print(f" Driven port active S11:") +print(f" @ 10.5 GHz : {S_dd_dB[i_op]:.2f} dB Z = {zin_d[i_op].real:.1f} + j{zin_d[i_op].imag:.1f} Ω") +# -10 dB BW around f0 +below = S_dd_dB <= -10.0 +if below[i_op]: + lo, hi = i_op, i_op + while lo > 0 and below[lo-1]: + lo -= 1 + while hi < len(below)-1 and below[hi+1]: + hi += 1 + print(f" -10 dB BW : {(freq[hi]-freq[lo])/1e6:.0f} MHz " + f"({freq[lo]/1e9:.2f} – {freq[hi]/1e9:.2f} GHz)") +else: + print(f" -10 dB BW : ") + +# Worst-case coupling (excluding driven port itself) +couplings = [(idx, abs(S[idx][i_op])) for idx in S.keys() if idx != (DRIVEN_X, DRIVEN_Y)] +couplings.sort(key=lambda x: -x[1]) +print() +print(f" Top-5 strongest couplings to driven port at 10.5 GHz:") +for idx, val in couplings[:5]: + di = idx[0] - DRIVEN_X + dj = idx[1] - DRIVEN_Y + dB = 20*np.log10(val + 1e-30) + print(f" ({idx[0]},{idx[1]}) Δ=({di:+d},{dj:+d}) |S| = {dB:>6.1f} dB") +print("=" * 70) + +# Save S matrix CSV (full-band) +with open(os.path.join(OUT_DIR, "S_matrix.csv"), "w", newline="") as f: + w = csv.writer(f) + header = ["freq_Hz"] + keys = sorted(S.keys()) + for idx in keys: + header += [f"S({idx[0]},{idx[1]})_dB", f"S({idx[0]},{idx[1]})_phase_deg"] + w.writerow(header) + for k in range(len(freq)): + row = [freq[k]] + for idx in keys: + mag_dB = 20*np.log10(np.abs(S[idx][k]) + 1e-30) + phase = np.angle(S[idx][k], deg=True) + row += [mag_dB, phase] + w.writerow(row) + +# Save driven-port S11 CSV +with open(os.path.join(OUT_DIR, "S11_data.csv"), "w", newline="") as f: + w = csv.writer(f) + w.writerow(["freq_Hz", "S11_dB", "Zin_real", "Zin_imag"]) + for k in range(len(freq)): + w.writerow([freq[k], S_dd_dB[k], zin_d[k].real, zin_d[k].imag]) + +# Coupling heatmap at 10.5 GHz +fig, ax = plt.subplots(figsize=(6.5, 6)) +grid = np.zeros((N_Y, N_X)) +for (i, j) in S.keys(): + grid[j, i] = 20*np.log10(abs(S[(i,j)][i_op]) + 1e-30) +# Driven port floor (S11 is just one number, not a coupling) — set to NaN to highlight +grid[DRIVEN_Y, DRIVEN_X] = np.nan +im = ax.imshow(grid, origin='lower', cmap='viridis', aspect='equal') +ax.set_xticks(range(N_X)) +ax.set_yticks(range(N_Y)) +ax.set_xlabel('i (x-pitch direction)') +ax.set_ylabel('j (y-pitch direction)') +ax.set_title(f'AERIS-10 4x4 array — coupling |S_j,({DRIVEN_X},{DRIVEN_Y})| at {F0/1e9:.2f} GHz') +# Annotate cells +for j in range(N_Y): + for i in range(N_X): + if (i, j) == (DRIVEN_X, DRIVEN_Y): + ax.text(i, j, "DRIVEN", ha='center', va='center', color='red', + fontsize=9, fontweight='bold') + else: + ax.text(i, j, f"{grid[j,i]:.1f}\ndB", ha='center', va='center', + color='white', fontsize=8) +plt.colorbar(im, ax=ax, label='|S| (dB)', shrink=0.7) +fig.tight_layout() +fig.savefig(os.path.join(OUT_DIR, "coupling_grid.png"), dpi=140) +plt.close(fig) + +print(f"[out] {OUT_DIR}/coupling_grid.png") +print(f"[out] {OUT_DIR}/S_matrix.csv") +print(f"[out] {OUT_DIR}/S11_data.csv")