diff --git a/5_Simulations/Antenna/probe_fed_aeris10_v3.py b/5_Simulations/Antenna/probe_fed_aeris10_v3.py new file mode 100644 index 0000000..3483f17 --- /dev/null +++ b/5_Simulations/Antenna/probe_fed_aeris10_v3.py @@ -0,0 +1,316 @@ +#!/usr/bin/env python3 +# probe_fed_aeris10_v3.py +# +# Single-element probe-fed patch antenna sim for AERIS-10 — true 2-layer +# stackup (L1 patch / 0.508 mm RO4350B / L2 ground plane). Probe via through +# ground plane feeds the patch; LumpedPort with R=50 Ω across the substrate +# at the probe location models the coax launch. +# +# Why this topology: aperture-coupled v2 (4-layer Stack_Hybrid) capped at +# ~60 MHz BW because the 0.11 mm L4 backshort acted as a near-short reflector +# — wider BW is fundamentally coupling-limited there. Probe-fed patch on the +# same 0.508 mm patch substrate has no slot bottleneck; physics BW from +# 3.77·(εr-1)/εr²·(W/L)·(h/λ₀) is ~1.6% ≈ 170 MHz at 10.5 GHz. +# +# Patch geometry preserved from the existing 8x16 Gerber (Antenna_16_8.top): +# W = 7.854 mm (D10 first dimension; sets X-pitch 14.27 mm in the array) +# L = 6.56 mm (tuned at balanced profile to land f_res = 10.51 GHz; old +# Gerber's 7.356 mm at 0.102 mm sub gave f_res ~10.6 GHz, the +# thicker substrate adds ~1 mm of fringing-edge ΔL each side) +# +# Probe location: +# y_off = 2.14 mm from -y radiating edge → R_in = 41 Ω, VSWR = 1.18 at 10.5 +# GHz. R_edge fitted from sim ≈ 152 Ω; cos²(π·y_off/L) gives R_in. +# +# Verified design point (PROFILE=balanced, λ/25 mesh, 13 s/run): +# f_res = 10.510 GHz, S11 = -21.79 dB at 10.5 GHz, Zin = 42.7 + j2.0 Ω +# -10 dB BW = 180 MHz (10.40 – 10.58 GHz, 1.71%) +# Compare 4-layer Stack_Hybrid + cap: 60 MHz BW, -19 dB. 3× wider, no cap. +# +# Stackup: +# L1 Cu 0.035 mm ← patch +# -- RO4350B 0.508 mm εr=3.48 (patch substrate) +# L2 Cu 0.035 mm ← ground plane (with antipad clearance) +# air below; coax launches up through +# to the probe via from L2 ground. +# +# Run: +# cd /tmp && DYLD_LIBRARY_PATH=/Users/ganeshpanth/opt/openEMS/lib \ +# PROFILE=balanced PATCH_L_MM=7.54 FEED_OFFSET_MM=2.5 \ +# /Users/ganeshpanth/radar_venv/bin/python \ +# /Users/ganeshpanth/PLFM_RADAR/5_Simulations/Antenna/probe_fed_aeris10_v3.py +# +# Profiles: +# sanity — λ/18 mesh, fast, borderline convergence +# balanced — λ/25 mesh, slower, recommended for design verification +# +# Env overrides (all optional): +# PATCH_W_MM PATCH_L_MM FEED_OFFSET_MM (mm from -y radiating edge) +# FEED_X_MM (default 0; offset along W-axis, normally 0 for centred feed) + +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 (mm) — true 2-layer probe-fed +# ============================================================================ +T_CU = 0.035 +H_PATCH_SUB = 0.508 # RO4350B between L1 patch and L2 ground +EPS_RO4350B = 3.48 +TAN_RO4350B = 0.0037 + +# Z layers (L2 ground at z=0, patch on top) +Z_GND = 0.0 +Z_PATCH = Z_GND + T_CU + H_PATCH_SUB +Z_TOP = Z_PATCH + T_CU + +# ============================================================================ +# GEOMETRY (mm) — defaults preserve old Gerber's W; L recomputed for 0.508 mm +# ============================================================================ +PATCH_W = float(os.environ.get("PATCH_W_MM", "7.854")) +PATCH_L = float(os.environ.get("PATCH_L_MM", "6.56")) + +# Probe: feed offset along the L-axis (y), measured from -y radiating edge +# inward. R_in(y_off) = R_edge·cos²(π·y_off/L). y_off=2.14 mm with iter#3 L +# (R_edge≈152 Ω fitted) lands R_in=41 Ω, VSWR=1.18 at 10.5 GHz. +FEED_OFFSET_MM = float(os.environ.get("FEED_OFFSET_MM", "2.14")) +FEED_X_MM = float(os.environ.get("FEED_X_MM", "0.0")) + +# Substrate / ground extents (~λ/2 margin around patch) +GND_X_MARGIN = 14.3 +GND_Y_MARGIN = 14.3 +GND_X_HALF = PATCH_W/2 + GND_X_MARGIN +GND_Y_HALF = PATCH_L/2 + GND_Y_MARGIN + +# Air box (λ/2 above patch, λ/2 below ground) +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_probefed_v3" +os.makedirs(OUT_DIR, exist_ok=True) + + +# ============================================================================ +# Build + run a single FDTD case +# ============================================================================ +def run_case(patch_w, patch_l, feed_offset, feed_x, sim_path, profile_cfg, label=""): + 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 ---- + patch_sub.AddBox([-GND_X_HALF, -GND_Y_HALF, Z_GND + T_CU], + [+GND_X_HALF, +GND_Y_HALF, Z_PATCH], priority=1) + + # ---- L1: patch (centred on origin, L along y, W along x) ---- + copper.AddBox([-patch_w/2, -patch_l/2, Z_PATCH], + [+patch_w/2, +patch_l/2, Z_PATCH + T_CU], priority=10) + + # ---- L2: full ground plane ---- + # Single-element sim — antipad clearance around the probe is implicit + # in the LumpedPort box (FDTD treats the port column as the metal probe + # and the surrounding cells as substrate). For a multi-element array + # with real coax launches a physical clearance hole would be added. + copper.AddBox([-GND_X_HALF, -GND_Y_HALF, Z_GND], + [+GND_X_HALF, +GND_Y_HALF, Z_GND + T_CU], priority=10) + + # ---- mesh ---- + lambda_min_mm = (C0 / F_STOP) * 1000.0 + res = lambda_min_mm / profile_cfg["mesh_lambda_div"] + + # Probe location in patch frame + feed_y = -patch_l/2 + feed_offset # offset from -y radiating edge + feed_x_pos = feed_x + + xlines = [-AIR_X_HALF, -GND_X_HALF, -patch_w/2, feed_x_pos, +patch_w/2, + +GND_X_HALF, +AIR_X_HALF] + ylines = [-AIR_Y_HALF, -GND_Y_HALF, -patch_l/2, feed_y, +patch_l/2, + +GND_Y_HALF, +AIR_Y_HALF] + # Z mesh: substrate gets ≥6 interior cells for accurate field capture + 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) + + if os.environ.get("MESH_DEBUG"): + z_diff = np.diff(zlines) + print(f"[mesh] x={len(xlines)} y={len(ylines)} z={len(zlines)} cells={n_cells:,}") + print(f"[mesh] z min/max/avg cell: {z_diff.min()*1e3:.1f}/{z_diff.max()*1e3:.1f}/{z_diff.mean()*1e3:.1f} um") + + # ---- LumpedPort: vertical 50 Ω port across substrate at (feed_x, feed_y) ---- + # The lumped port replaces the coax+source: the 50 Ω resistor sits in the + # box, the metal column from L2 to L1 is implicit. Excitation is z-direction + # E-field across the substrate. + port_start = [feed_x_pos, feed_y, Z_GND + T_CU] + port_stop = [feed_x_pos, feed_y, Z_PATCH] + port = fdtd.AddLumpedPort(1, 50, port_start, port_stop, 'z', + excite=1.0, priority=5) + + # ---- run ---- + print(f"[case {label}] patch={patch_w:.2f}x{patch_l:.2f}mm " + f"feed=({feed_x_pos:.2f},{feed_y:.2f})mm cells={n_cells:,}") + 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) + port.CalcPort(sim_path, freq) + s11 = port.uf_ref / port.uf_inc + s11_dB = 20.0 * np.log10(np.abs(s11) + 1e-30) + zin = port.uf_tot / port.if_tot + vswr = (1 + np.abs(s11)) / (1 - np.abs(s11) + 1e-30) + return freq, s11_dB, zin, vswr, dt + + +def find_resonance(freq, s11_dB, zin=None): + """Resonance: where R peaks AND Im(Z)=0 nearby. Falls back to min(S11).""" + f_res, s11_min = None, None + if zin is not None: + # Find peak R in the band + mask = (freq >= 9.0e9) & (freq <= 11.5e9) + idx_band = np.where(mask)[0] + if len(idx_band) > 1: + r_band = np.real(zin[idx_band]) + i_pk = idx_band[int(np.argmax(r_band))] + # Use the X=0 crossing closest to the R peak + x = np.imag(zin) + sign = np.sign(x) + crossings = np.where(np.diff(sign) != 0)[0] + crossings_in_band = [c for c in crossings if mask[c]] + if crossings_in_band: + k = min(crossings_in_band, key=lambda c: abs(c - i_pk)) + t = -x[k] / (x[k+1] - x[k]) if x[k+1] != x[k] else 0 + f_res = freq[k] + t * (freq[k+1] - freq[k]) + s11_min = s11_dB[k] + t * (s11_dB[k+1] - s11_dB[k]) + + if f_res is None: + imin = int(np.argmin(s11_dB)) + f_res = freq[imin] + s11_min = float(s11_dB[imin]) + + # walk outward to find -10 dB crossings around f_res + below = s11_dB <= -10.0 + if not below.any(): + return f_res, s11_min, 0.0, 0.0, 0.0 + i_f = int(np.argmin(np.abs(freq - f_res))) + if not below[i_f]: + return f_res, s11_min, 0.0, 0.0, 0.0 + lo = i_f + while lo > 0 and below[lo-1]: + lo -= 1 + hi = i_f + while hi < len(below)-1 and below[hi+1]: + hi += 1 + f_lo, f_hi = freq[lo], freq[hi] + bw = f_hi - f_lo + bw_pct = bw / f_res * 100.0 + return f_res, s11_min, f_lo, f_hi, bw_pct + + +# ============================================================================ +# MAIN +# ============================================================================ +sim_path = os.path.join(OUT_DIR, "single") +freq, s11_dB, zin, vswr, dt = run_case( + PATCH_W, PATCH_L, FEED_OFFSET_MM, FEED_X_MM, sim_path, cfg) +f_res, s11_min, f_lo, f_hi, bw_pct = find_resonance(freq, s11_dB, zin) +i_res = int(np.argmin(np.abs(freq - f_res))) +i_op = int(np.argmin(np.abs(freq - 10.5e9))) + +print() +print("=" * 70) +print(f" Resonance (R peak + Im=0): {f_res/1e9:.3f} GHz (target 10.5 GHz)") +print(f" S11 at resonance : {s11_min:.2f} dB") +print(f" Zin at resonance : {zin[i_res].real:.1f} + j{zin[i_res].imag:.1f} Ω") +print(f" ── at 10.500 GHz exactly:") +print(f" S11 @ 10.5GHz : {s11_dB[i_op]:.2f} dB") +print(f" Zin @ 10.5GHz : {zin[i_op].real:.1f} + j{zin[i_op].imag:.1f} Ω") +print(f" VSWR @ 10.5GHz : {vswr[i_op]:.2f}") +print(f" -10 dB bandwidth : {(f_hi-f_lo)/1e6:.0f} MHz " + f"({f_lo/1e9:.3f} – {f_hi/1e9:.3f} GHz, {bw_pct:.2f}%)") +print(f" Sim time : {dt:.1f} s") +print("=" * 70) + +fig, ax = plt.subplots(figsize=(8.5, 4.5)) +ax.plot(freq/1e9, s11_dB, "b-", lw=1.6, label="S11") +ax.axhline(-10, color="r", ls="--", lw=0.8, label="-10 dB") +ax.axvline(f_res/1e9, color="g", ls=":", lw=0.8, + label=f"resonance {f_res/1e9:.3f} GHz") +if (f_hi-f_lo) > 0: + ax.axvspan(f_lo/1e9, f_hi/1e9, color="g", alpha=0.10, + label=f"BW {(f_hi-f_lo)/1e6:.0f} MHz ({bw_pct:.2f}%)") +ax.set_xlabel("Frequency (GHz)") +ax.set_ylabel("S11 (dB)") +ax.set_title(f"AERIS-10 Probe-Fed Patch v3 — 2-layer 0.508 mm RO4350B " + f"(W={PATCH_W} L={PATCH_L} y_off={FEED_OFFSET_MM}mm)") +ax.set_xlim(F_START/1e9, F_STOP/1e9) +ax.set_ylim(-40, 0) +ax.grid(True, alpha=0.3) +ax.legend(loc="lower right") +fig.tight_layout() +fig.savefig(os.path.join(OUT_DIR, "S11.png"), dpi=140) +plt.close(fig) + +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", "VSWR"]) + for k in range(len(freq)): + w.writerow([freq[k], s11_dB[k], zin[k].real, zin[k].imag, vswr[k]]) + +print(f"[out] {OUT_DIR}/S11.png") +print(f"[out] {OUT_DIR}/S11_data.csv")