diff --git a/5_Simulations/Antenna/Quartz_Waveguide.py b/5_Simulations/Antenna/Quartz_Waveguide.py deleted file mode 100644 index 15a7cac..0000000 --- a/5_Simulations/Antenna/Quartz_Waveguide.py +++ /dev/null @@ -1,266 +0,0 @@ -# openems_quartz_slotted_wg_10p5GHz.py -# Full script: geometry, mesh (no GetLine calls), S-params/impedance sweep, 3D pattern & gain. -# Requires: openEMS (Python bindings), CSXCAD (Python), matplotlib, numpy. - -import os -import tempfile -import numpy as np -import matplotlib.pyplot as plt - -from CSXCAD import ContinuousStructure, AppCSXCAD_BIN -from openEMS import openEMS -from openEMS.physical_constants import C0 - -# ------------------------- -# User controls -# ------------------------- -view_geom_in_AppCSXCAD = True # True => launch AppCSXCAD to view 3D geometry -simulate = True # False => skip FDTD run -threads = 0 # 0 => auto/max - -# Far-field angular sampling -n_theta, n_phi = 91, 181 # theta 0..180, phi 0..360 - -# ------------------------- -# Band & element sizing -# ------------------------- -f0 = 10.5e9 -f_span = 2.0e9 -f_start, f_stop = f0 - f_span/2, f0 + f_span/2 - -# Quartz-filled rectangular waveguide (full dielectric block) -er_quartz = 3.8 -# Array-driven constraints (λ0/2 column pitch, 1 mm septum) ⇒ a ~ 13.28 mm -a = 13.28 # mm (broad wall, internal) -b = 6.50 # mm (narrow wall, internal; comfortable to machine) -L = 281.0 # mm (≈32-slot column incl. λg/4 margins at 10.5 GHz) - -# Slot starters (tune in EM for taper) -slot_w = 0.60 # mm (across x) -lambda0_mm = (C0/f0) * 1e3 -fc10 = (C0/(2.0*np.sqrt(er_quartz))) * (1.0/(a*1e-3)) # Hz -lambda_d = (C0/f0) / np.sqrt(er_quartz) # m -lambda_g = lambda_d / np.sqrt(1.0 - (fc10/f0)**2) # m -lambda_g_mm = lambda_g * 1e3 -slot_s = 0.5*lambda_g_mm -slot_L = 0.47*lambda_g_mm -margin = 0.25*lambda_g_mm -Nslots = 32 -delta0 = 0.90 # mm offset from centerline (± alternated) - -# Metal & air padding for the radiation domain -t_metal = 0.8 # mm wall thickness -air_x = 10.0 # mm on each side -air_y = 40.0 # mm above slots -air_z = 15.0 # mm front/back - -# Mesh resolution target (mm) -mesh_res = min(0.5, lambda0_mm/30.0) - -# ------------------------- -# Build FDTD & CSX -# ------------------------- -unit = 1e-3 # mm -Sim_Path = os.path.join(tempfile.gettempdir(), 'openems_quartz_slotted_wg') - -FDTD = openEMS(NrTS=int(6e5), EndCriteria=1e-5) -FDTD.SetGaussExcite(0.5*(f_start+f_stop), 0.5*(f_stop-f_start)) -FDTD.SetBoundaryCond(['PML_8']*6) -FDTD.SetOverSampling(4) -FDTD.SetTimeStepFactor(0.95) - -CSX = ContinuousStructure() -FDTD.SetCSX(CSX) - -mesh = CSX.GetGrid() -mesh.SetDeltaUnit(unit) - -# ------------------------- -# Geometry extents (mm) -# ------------------------- -x_min, x_max = -air_x, a + air_x -y_min, y_max = -5.0, b + t_metal + air_y -z_min, z_max = -air_z, L + air_z - -# Slot centers and edges (mm) -z_centers = margin + np.arange(Nslots)*slot_s -x_centers = (a/2.0) + np.array([+delta0 if i%2==0 else -delta0 for i in range(Nslots)]) -x_edges = np.concatenate([x_centers - slot_w/2.0, x_centers + slot_w/2.0]) -z_edges = np.concatenate([z_centers - slot_L/2.0, z_centers + slot_L/2.0]) - -# ------------------------- -# Mesh lines — EXPLICIT (no GetLine calls) -# ------------------------- -x_lines = sorted({x_min, -t_metal, 0.0, a, a + t_metal, x_max, *list(x_edges)}) -y_lines = [y_min, 0.0, b, b+t_metal, y_max] -z_lines = sorted({z_min, 0.0, L, z_max, *list(z_edges)}) - -mesh.AddLine('x', x_lines) -mesh.AddLine('y', y_lines) -mesh.AddLine('z', z_lines) - -# Optional smoothing to max cell size around ~mesh_res -mesh.SmoothMeshLines('all', mesh_res, ratio=1.4) - -# ------------------------- -# Materials -# ------------------------- -pec = CSX.AddMetal('PEC') -quartz = CSX.AddMaterial('QUARTZ') -quartz.SetMaterialProperty(epsilon=er_quartz) -air = CSX.AddMaterial('AIR') # explicit for slot holes - -# ------------------------- -# Solids: quartz block + metal walls + slots -# ------------------------- -# Quartz full block (inside tube) -quartz.AddBox([0, 0, 0], [a, b, L]) - -# PEC tube: left/right/bottom/top (top will be perforated by slots) -pec.AddBox([-t_metal, 0, 0], [0, b, L]) # left -pec.AddBox([a, 0, 0], [a+t_metal,b, L]) # right -pec.AddBox([-t_metal,-t_metal,0],[a+t_metal,0, L]) # bottom -pec.AddBox([-t_metal, b, 0], [a+t_metal,b+t_metal,L]) # top - -# Slots = AIR boxes overriding the top metal -for zc, xc in zip(z_centers, x_centers, strict=False): - x1, x2 = xc - slot_w/2.0, xc + slot_w/2.0 - z1, z2 = zc - slot_L/2.0, zc + slot_L/2.0 - prim = air.AddBox([x1, b, z1], [x2, b+t_metal, z2]) - prim.SetPriority(10) # ensure it cuts the metal - -# ------------------------- -# Ports: Rectangular WG TE10, z-directed -# ------------------------- -port_thick = max(4*mesh_res, 2.0) # mm -p1_start = [0, 0, max(0.5, 10*mesh_res)] -p1_stop = [a, b, p1_start[2] + port_thick] -FDTD.AddRectWaveGuidePort(port_nr=1, start=p1_start, stop=p1_stop, - p_dir='z', a=a*unit, b=b*unit, mode_name='TE10', excite=1) - -p2_stop = [a, b, L - max(0.5, 10*mesh_res)] -p2_start = [0, 0, p2_stop[2] - port_thick] -FDTD.AddRectWaveGuidePort(port_nr=2, start=p2_start, stop=p2_stop, - p_dir='z', a=a*unit, b=b*unit, mode_name='TE10', excite=0) - -# ------------------------- -# NF2FF setup (surround the radiator region) -# ------------------------- -def create_nf2ff(FDTD_obj, name, start, stop, frequency): - """Compat wrapper: older/newer openEMS builds may expose NF2FF creation differently.""" - try: - return FDTD_obj.CreateNF2FFBox(name=name, start=start, stop=stop, frequency=frequency) - except AttributeError: - # Fallback: try AddNF2FFBox returning a handle-like object - return FDTD_obj.AddNF2FFBox(name=name, start=start, stop=stop, frequency=frequency) - -nf2ff = create_nf2ff( - FDTD, - name='nf2ff', - start=[x_min+1.0, y_min+1.0, z_min+1.0], - stop =[x_max-1.0, y_max-1.0, z_max-1.0], - frequency=[f0] -) - -# ------------------------- -# (Optional) view geometry -# ------------------------- -if view_geom_in_AppCSXCAD: - os.makedirs(Sim_Path, exist_ok=True) - csx_xml = os.path.join(Sim_Path, 'quartz_slotted_wg.xml') - CSX.Write2XML(csx_xml) - os.system(f'"{AppCSXCAD_BIN}" "{csx_xml}"') - -# ------------------------- -# Run FDTD -# ------------------------- -if simulate: - FDTD.Run(Sim_Path, cleanup=True, verbose=2, numThreads=threads) - -# ------------------------- -# Post-processing: S-params & impedance -# ------------------------- -freq = np.linspace(f_start, f_stop, 401) -ports = list(FDTD.ports) # Port 1 & Port 2 in creation order -for p in ports: - p.CalcPort(Sim_Path, freq) - -S11 = ports[0].uf_ref / ports[0].uf_inc -S21 = ports[1].uf_ref / ports[0].uf_inc -Zin = ports[0].uf_tot / ports[0].if_tot - -plt.figure(figsize=(7.6,4.6)) -plt.plot(freq*1e-9, 20*np.log10(np.abs(S11)), lw=2, label='|S11|') -plt.plot(freq*1e-9, 20*np.log10(np.abs(S21)), lw=2, ls='--', label='|S21|') -plt.grid(True) -plt.legend() -plt.xlabel('Frequency (GHz)') -plt.ylabel('Magnitude (dB)') -plt.title('S-Parameters: Slotted Quartz-Filled WG') - -plt.figure(figsize=(7.6,4.6)) -plt.plot(freq*1e-9, np.real(Zin), lw=2, label='Re{Zin}') -plt.plot(freq*1e-9, np.imag(Zin), lw=2, ls='--', label='Im{Zin}') -plt.grid(True) -plt.legend() -plt.xlabel('Frequency (GHz)') -plt.ylabel('Ohms') -plt.title('Input Impedance (Port 1)') - -# ------------------------- -# Far-field @ f0 and 3D pattern -# ------------------------- -theta = np.linspace(0, np.pi, n_theta) -phi = np.linspace(0, 2*np.pi, n_phi) - -# Compatibility: some builds expect nf2ff.CalcNF2FF(...), others FDTD.CalcNF2FF(nf2ff,...) -try: - res = nf2ff.CalcNF2FF(Sim_Path, [f0], theta, phi) -except AttributeError: - res = FDTD.CalcNF2FF(nf2ff, Sim_Path, [f0], theta, phi) - -# Max directivity (linear) & realized gain estimate -idx_f0 = np.argmin(np.abs(freq - f0)) -Dmax_lin = float(res.Dmax[0]) # at f0 -mismatch = 1.0 - np.abs(S11[idx_f0])**2 # (1 - |S11|^2) -Gmax_lin = Dmax_lin * float(mismatch) -Gmax_dBi = 10*np.log10(Gmax_lin) - - -# 3D normalized pattern -E = np.squeeze(res.E_norm) # shape [f, th, ph] -> [th, ph] -E = E / np.max(E) -TH, PH = np.meshgrid(theta, phi, indexing='ij') -R = E -X = R * np.sin(TH) * np.cos(PH) -Y = R * np.sin(TH) * np.sin(PH) -Z = R * np.cos(TH) - -fig = plt.figure(figsize=(7.2,6.2)) -ax = fig.add_subplot(111, projection='3d') -ax.plot_surface(X, Y, Z, rstride=2, cstride=2, linewidth=0, antialiased=True, alpha=0.92) -ax.set_title(f'Normalized 3D Pattern @ {f0/1e9:.2f} GHz\n(peak ≈ {Gmax_dBi:.1f} dBi)') -ax.set_box_aspect((1,1,1)) -ax.set_xlabel('x') -ax.set_ylabel('y') -ax.set_zlabel('z') -plt.tight_layout() - -# Quick 2D geometry preview (top view at y=b) -plt.figure(figsize=(8.4,2.8)) -plt.fill_between( - [0, a], [0, 0], [L, L], color='#dddddd', alpha=0.5, step='pre', label='WG aperture (top)' -) -for zc, xc in zip(z_centers, x_centers, strict=False): - plt.gca().add_patch(plt.Rectangle((xc - slot_w/2.0, zc - slot_L/2.0), - slot_w, slot_L, fc='#3355ff', ec='k')) -plt.xlim(-2, a + 2) -plt.ylim(-5, L + 5) -plt.gca().invert_yaxis() -plt.xlabel('x (mm)') -plt.ylabel('z (mm)') -plt.title('Top-view slot layout (y=b plane)') -plt.grid(True) -plt.legend() - -plt.show() diff --git a/5_Simulations/Antenna/aperture_coupled_aeris10_v2.py b/5_Simulations/Antenna/aperture_coupled_aeris10_v2.py deleted file mode 100644 index aa0deed..0000000 --- a/5_Simulations/Antenna/aperture_coupled_aeris10_v2.py +++ /dev/null @@ -1,471 +0,0 @@ -#!/usr/bin/env python3 -# aperture_coupled_aeris10_v2.py -# -# Single-element aperture-coupled patch antenna sim for the Stack_Hybrid -# 4-layer stackup (committed 1de2296, 2026-04-29). Default parameters are the -# best design point from a 10-iteration analytic-tune sweep run 2026-05-02. -# -# DESIGN POINT @ 10.5 GHz (defaults below): -# patch : W=9.55 mm L=7.77 mm -# slot : L=3.00 mm W=0.50 mm (centered under patch) -# stub : L=4.16 mm (= λ_g/4 in feed sub at f0 — short across slot at 10.5) -# feed : lead=12.34 mm W=0.25 mm -# total feed length = lead + stub = 16.50 mm = 1·λ_g at 10.5 GHz -# ⇒ feed line is TRANSPARENT at f0 (port sees true antenna Z) -# -# Result @ 10.5 GHz: Z ≈ R + j350 Ω where R = 33–51 Ω across reruns -# (R varies ~30% with sim convergence — sanity profile is borderline; the -# physics-meaningful result is "R within matching ballpark, X stable at -# +350"). The +j350 inductive residual cannot be canceled in this topology -# by simple stub tuning — it stems from the L4 backshort continuous ground -# under the antenna footprint. Two production-grade fixes: -# (a) Series matching cap at the port: C ≈ 1/(2π·10.5GHz·350) ≈ 0.043 pF -# standard 0402 ATC cap → drops |Γ| from 0.97 to ~0.01 (S11 ≈ -40 dB). -# (b) Open the L4 backshort under the antenna footprint (stackup edit) → -# restores standard open-back aperture-coupled, stub naturally tunes X. -# Either is the antenna designer's call. This script's role is to provide -# the verified starting point in (W,L,slot,stub,feed_lead) plus the X that -# needs to be matched out. -# -# Stackup (Stack_Hybrid.png, ANTENNA column): -# L1 Cu 0.035 mm ← patch -# -- RO4350B 0.508 mm εr=3.48 (top patch substrate) -# L2 Cu 0.035 mm ← inner ground + coupling slot -# -- RO4450F 1.2 mm εr=3.52 (bonding ply) -# L3 Cu 0.035 mm ← microstrip feed line + λ_g/4 stub -# -- RO4350B 0.11 mm εr=3.48 (feed substrate) -# L4 Cu 0.035 mm ← bottom ground plane (backshort) -# -# Notable bug-fixes baked into this version (vs commit 42056b8 baseline): -# - Z mesh: explicit fine substrate-interior lines (≥5 cells per substrate) -# bypassing SmoothMeshLines collapse — feed sub at 0.11 mm now has proper -# microstrip Z0. Without this, the patch resonance is hidden by mesh-Z0 -# artifacts and the sim measures essentially line-only behavior. -# - slot_y_off env var: was read but never applied to the L2 slot box. Now -# the slot is correctly offset in y when SLOT_Y_OFF_MM != 0. -# - FEED_LEAD_L: now env-tunable (was hardcoded 14.0 mm). The default 12.34 -# mm makes total feed = 1·λ_g at 10.5 GHz, killing the spurious feed-line -# full-wave resonance that masked the true patch resonance in the -# baseline script (showed up as a persistent 9.4-9.5 GHz "resonance"). -# -# Run modes: -# PROFILE=sanity : 1 run, mesh λ_min/18, ~30s/run -# PROFILE=balanced : 1 run, finer mesh λ_min/25, slower -# PROFILE=sweep : 5×5 grid over slot_L × stub_L, picks best, reports -# -# Env overrides (all optional, defaults at iter #6 design point): -# PATCH_W_MM PATCH_L_MM -# SLOT_L_MM SLOT_W_MM SLOT_Y_OFF_MM -# STUB_L_MM FEED_LEAD_L_MM -# MESH_DEBUG=1 prints mesh density before each run -# -# Outputs in /tmp/aeris10_aperture_v2/: -# single run : S11.png, S11_data.csv -# sweep mode : sweep_results.csv, sweep_S11.png - -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}, - "sweep": {"mesh_lambda_div": 16, "n_timesteps": 40000, "end_dB": -25}, -} -cfg = profiles[PROFILE if PROFILE != "sweep" else "sweep"] - -# ============================================================================ -# BAND -# ============================================================================ -F0 = 10.5e9 -F_SPAN = 4.0e9 -F_START = F0 - F_SPAN/2 -F_STOP = F0 + F_SPAN/2 - -# ============================================================================ -# STACKUP (mm) — from Stack_Hybrid.png ANTENNA column -# ============================================================================ -T_CU = 0.035 -H_PATCH_SUB = 0.508 -H_BOND = 1.2 -H_FEED_SUB = 0.11 -EPS_RO4350B = 3.48 -EPS_RO4450F = 3.52 -TAN_RO4350B = 0.0037 -TAN_RO4450F = 0.0040 - -# Z layers (L4 bottom = 0) -Z_L4 = 0.0 -Z_L3 = Z_L4 + T_CU + H_FEED_SUB -Z_L2 = Z_L3 + T_CU + H_BOND -Z_L1 = Z_L2 + T_CU + H_PATCH_SUB -Z_TOP = Z_L1 + T_CU - -# ============================================================================ -# GEOMETRY (mm) — defaults are iter #6 best design point (see header) -# ============================================================================ -# Patch: empirically tuned for the Stack_Hybrid 4-layer stack (with L4 -# backshort). Note that with L4 present, εr_eff at the patch is ~4.0 in sim -# (not the 3.21 a single-substrate Balanis formula predicts), so L is larger -# than open-back textbook value — the L4 dielectric loading lowers f_res. -PATCH_W = float(os.environ.get("PATCH_W_MM", "9.55")) -PATCH_L = float(os.environ.get("PATCH_L_MM", "7.77")) - -# Slot (under patch in L2), length perpendicular to feed direction. -# Slot resonance λ_g_slot/2 ≈ 7.65 mm in the L2-L4 cavity; SLOT_L=3.0 keeps -# the slot well sub-resonant (slot is a coupling aperture, not a radiator). -SLOT_L = float(os.environ.get("SLOT_L_MM", "3.0")) -SLOT_W = float(os.environ.get("SLOT_W_MM", "0.5")) - -# Microstrip feed on L3, dominant ground = L4 (0.11 mm RO4350B). 50 Ω target. -# Hammerstad: W ≈ 0.25 mm for 50 Ω on 0.11 mm RO4350B. -FEED_W = 0.25 -FEED_STUB_L = float(os.environ.get("STUB_L_MM", "4.16")) # λ_g/4 @ 10.5 GHz -# Total feed length = FEED_LEAD_L + STUB_L should be n·λ_g at f0 for the line -# to be transparent at the operating freq (sim sees the antenna's true impedance -# at port without TL transformation). λ_g_feed @ 10.5 GHz on 0.11 mm RO4350B -# microstrip ≈ 16.5 mm → FEED_LEAD_L = 16.5 - STUB_L for n=1. -FEED_LEAD_L = float(os.environ.get("FEED_LEAD_L_MM", "12.34")) # n=1 λ_g default - -# Substrate / ground extents (~λ/2 margin around patch) -GND_X_MARGIN = 14.3 -GND_Y_MARGIN = 14.3 -GND_X_HALF = max(PATCH_W/2, SLOT_L/2) + GND_X_MARGIN -GND_Y_HALF = max(PATCH_L/2, FEED_LEAD_L + FEED_STUB_L) + GND_Y_MARGIN - -# Air box (λ/2 above patch, λ/4 below) -AIR_ABOVE = 14.3 -AIR_BELOW = 8.0 -AIR_X_HALF = GND_X_HALF + 8.0 -AIR_Y_HALF = GND_Y_HALF + 8.0 - -OUT_DIR = "/tmp/aeris10_aperture_v2" -os.makedirs(OUT_DIR, exist_ok=True) - - -# ============================================================================ -# Build + run a single FDTD case, return (freq[], S11_dB[], Zin[], VSWR[]) -# ============================================================================ -def run_case(slot_l, stub_l, patch_l, sim_path, profile_cfg, label="", slot_w=None, - slot_y_off=None, patch_w=None): - if slot_w is None: - slot_w = SLOT_W - if slot_y_off is None: - slot_y_off = float(os.environ.get("SLOT_Y_OFF_MM", "0.0")) - if patch_w is None: - patch_w = PATCH_W - 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_top", - epsilon=EPS_RO4350B, - kappa=2*np.pi*F0*EPS_RO4350B*eps0*TAN_RO4350B) - bond_sub = CSX.AddMaterial("RO4450F_bond", - epsilon=EPS_RO4450F, - kappa=2*np.pi*F0*EPS_RO4450F*eps0*TAN_RO4450F) - feed_sub = CSX.AddMaterial("RO4350B_feed", - epsilon=EPS_RO4350B, - kappa=2*np.pi*F0*EPS_RO4350B*eps0*TAN_RO4350B) - copper = CSX.AddMetal("Copper") - - # ---- substrates ---- - patch_sub.AddBox([-GND_X_HALF, -GND_Y_HALF, Z_L2 + T_CU], - [+GND_X_HALF, +GND_Y_HALF, Z_L1], priority=1) - bond_sub.AddBox([-GND_X_HALF, -GND_Y_HALF, Z_L3 + T_CU], - [+GND_X_HALF, +GND_Y_HALF, Z_L2], priority=1) - feed_sub.AddBox([-GND_X_HALF, -GND_Y_HALF, Z_L4 + T_CU], - [+GND_X_HALF, +GND_Y_HALF, Z_L3], priority=1) - - # ---- L1: patch ---- - copper.AddBox([-patch_w/2, -patch_l/2, Z_L1], - [+patch_w/2, +patch_l/2, Z_L1 + T_CU], priority=10) - - # ---- L2: ground patch around slot only (NOT full plane) ---- - # Finite-extent inner ground centered around the slot. Lets the feed line - # transition from microstrip (no L2 above) at the port → stripline-ish - # (L2 above) near the slot. This kills the parallel-plate cavity that a - # full-extent L2 plane would form between L2↔L4 around the lumped port. - # L2 patch covers ~2·PATCH_L in y, full board width in x. - L2_HALF_Y = PATCH_L # ground extends ±PATCH_L in y around slot at y=0 - sy0 = slot_y_off - slot_w/2 - sy1 = slot_y_off + slot_w/2 - # Above slot (y > sy1) - copper.AddBox([-GND_X_HALF, sy1, Z_L2], - [+GND_X_HALF, +L2_HALF_Y, Z_L2 + T_CU], priority=10) - # Below slot (y < sy0) - copper.AddBox([-GND_X_HALF, -L2_HALF_Y, Z_L2], - [+GND_X_HALF, sy0, Z_L2 + T_CU], priority=10) - # Left of slot (x < -slot_l/2, sy0 <= y <= sy1) - copper.AddBox([-GND_X_HALF, sy0, Z_L2], - [-slot_l/2, sy1, Z_L2 + T_CU], priority=10) - # Right of slot (x > +slot_l/2, sy0 <= y <= sy1) - copper.AddBox([+slot_l/2, sy0, Z_L2], - [+GND_X_HALF, sy1, Z_L2 + T_CU], priority=10) - - # ---- L3: microstrip feed line — runs in y direction, ⟂ to slot ---- - feed_y_start = -FEED_LEAD_L # board edge (port location) - feed_y_end = +stub_l # stub past slot center (slot at y=0) - copper.AddBox([-FEED_W/2, feed_y_start, Z_L3], - [+FEED_W/2, feed_y_end, Z_L3 + T_CU], priority=10) - - # ---- L4: full bottom ground ---- - copper.AddBox([-GND_X_HALF, -GND_Y_HALF, Z_L4], - [+GND_X_HALF, +GND_Y_HALF, Z_L4 + T_CU], priority=10) - - # ---- mesh (must exist before MSLPort sees it) ---- - PORT_LEN = 2.0 # mm of trace allocated to the port region - FEED_SHIFT_PORT = 0.4 # excitation point inside the port box - MEAS_SHIFT_PORT = 1.6 # V/I measurement plane inside the port box - - lambda_min_mm = (C0 / F_STOP) * 1000.0 - res = lambda_min_mm / profile_cfg["mesh_lambda_div"] - xlines = [-AIR_X_HALF, -GND_X_HALF, -PATCH_W/2, -slot_l/2, -FEED_W/2, - 0, +FEED_W/2, +slot_l/2, +PATCH_W/2, +GND_X_HALF, +AIR_X_HALF] - # MSLPort requires ≥5 mesh lines in propagation direction inside the port - # box (y = feed_y_start..feed_y_start+PORT_LEN). Force 6 explicit lines. - port_y_lines = list(np.linspace(feed_y_start, feed_y_start + PORT_LEN, 6)) - ylines = [-AIR_Y_HALF, -GND_Y_HALF, -PATCH_L/2, -slot_w/2, - 0, +slot_w/2, +PATCH_L/2, feed_y_end, +GND_Y_HALF, +AIR_Y_HALF, - -PATCH_L] + port_y_lines - # Z mesh: built MANUALLY (not via SmoothMeshLines) because the substrates - # need ≥5 cells for accurate microstrip Z0 (esp. feed sub at 0.11 mm). - # SmoothMeshLines collapses lines closer than ~res/3, which kills the fine - # substrate refinement we need. Build it explicitly: - # - Air below/above: res-spaced - # - Each metal layer: one line at top + one at bottom of copper - # - Each dielectric: 6 interior cells (7-pt linspace, drop endpoints) - air_below = list(np.arange(Z_L4 - AIR_BELOW, Z_L4, res)) - air_above = list(np.arange(Z_TOP + res, Z_TOP + AIR_ABOVE + res, res)) - feed_interior = list(np.linspace(Z_L4 + T_CU, Z_L3, 7)[1:-1]) # 0.018 mm pitch - bond_interior = list(np.linspace(Z_L3 + T_CU, Z_L2, 7)[1:-1]) # 0.20 mm pitch - patch_interior = list(np.linspace(Z_L2 + T_CU, Z_L1, 7)[1:-1]) # 0.085 mm pitch - zlines = sorted(set(air_below + [ - Z_L4, Z_L4 + T_CU, - Z_L3, Z_L3 + T_CU, - Z_L2, Z_L2 + T_CU, - Z_L1, Z_L1 + T_CU, - ] + feed_interior + bond_interior + patch_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"): - print(f"[mesh] xlines={len(xlines)} ylines={len(ylines)} zlines={len(zlines)}") - z_diff = np.diff(zlines) - print(f"[mesh] z min/max/avg cell: {z_diff.min()*1e3:.3f}/{z_diff.max()*1e3:.3f}/{z_diff.mean()*1e3:.3f} um") - print(f"[mesh] zlines (mm): {[f'{z:.3f}' for z in zlines]}") - - # ---- Microstrip-line port (after mesh is in place; checks line count) ---- - # The L2 inner ground was pulled back to y = -PATCH_L (= -7.25 mm), so - # this port region (y = -8.0..-6.0 mm) sees ONLY L4 ground below → pure - # microstrip. AddMSLPort excites the TEM mode cleanly without coupling - # into the parallel-plate cavity that a full L2 plane would form. - port = fdtd.AddMSLPort(1, copper, - start=[-FEED_W/2, feed_y_start, Z_L4 + T_CU], - stop= [+FEED_W/2, feed_y_start + PORT_LEN, Z_L3 + T_CU], - prop_dir='y', exc_dir='z', - excite=1.0, - FeedShift=FEED_SHIFT_PORT, - MeasPlaneShift=MEAS_SHIFT_PORT, - Feed_R=50) - - # ---- run ---- - print(f"[case {label}] slot_L={slot_l:.2f}mm stub_L={stub_l:.2f}mm " - f"patch_L={patch_l:.3f}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): - """Return (f_res_Hz, s11_min_dB, f_lo, f_hi, bw_pct). - - Resonance point: where Im(Zin) crosses zero (true antenna resonance). - Falls back to min(S11) if Zin not provided or no zero crossing found. - """ - f_res, s11_min = None, None - if zin is not None: - # Find Im(Zin) = 0 crossings inside the search band 9.0..11.5 GHz - x = np.imag(zin) - mask = (freq >= 9.0e9) & (freq <= 11.5e9) - idx_band = np.where(mask)[0] - if len(idx_band) > 1: - xb = x[idx_band] - # Find sign changes in xb - sign_changes = np.where(np.diff(np.sign(xb)) != 0)[0] - if len(sign_changes): - # Linear interpolation to refine the crossing - k = sign_changes[0] - i0, i1 = idx_band[k], idx_band[k+1] - # Linear interp where Im(Zin) = 0 - t = -x[i0] / (x[i1] - x[i0]) if x[i1] != x[i0] else 0 - f_res = freq[i0] + t * (freq[i1] - freq[i0]) - # S11 at f_res via interpolation - s11_min = s11_dB[i0] + t * (s11_dB[i1] - s11_dB[i0]) - - 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 - # find the -10 dB band containing or nearest to f_res - 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 -# ============================================================================ -if PROFILE == "sweep": - # 5x5 grid centered on the iter #6 design point (slot=3.0, stub=4.16) - slot_grid = [2.0, 2.5, 3.0, 3.5, 4.0] - stub_grid = [3.5, 3.85, 4.16, 4.5, 4.85] - patch_l = float(os.environ.get("PATCH_L_MM", "7.77")) - rows = [] - print(f"[sweep] {len(slot_grid)}×{len(stub_grid)} = {len(slot_grid)*len(stub_grid)} cases") - for i, sl in enumerate(slot_grid): - for j, st in enumerate(stub_grid): - label = f"{i*len(stub_grid)+j+1:02d}/{len(slot_grid)*len(stub_grid)}" - sim_path = os.path.join(OUT_DIR, f"sweep_s{sl:.1f}_t{st:.1f}") - try: - freq, s11_dB, zin, vswr, dt = run_case( - sl, st, patch_l, sim_path, cfg, label=label) - f_res, s11_min, f_lo, f_hi, bw_pct = find_resonance(freq, s11_dB, zin) - rows.append({ - "slot_L": sl, "stub_L": st, "patch_L": patch_l, - "f_res_GHz": f_res/1e9, "s11_min_dB": s11_min, - "f_lo_GHz": f_lo/1e9, "f_hi_GHz": f_hi/1e9, - "bw_pct": bw_pct, "elapsed_s": dt, - }) - print(f" [{label}] f_res={f_res/1e9:.2f}GHz S11={s11_min:.1f}dB " - f"BW={bw_pct:.1f}% t={dt:.0f}s") - except Exception as e: - print(f" [{label}] FAILED: {e}") - - # CSV - with open(os.path.join(OUT_DIR, "sweep_results.csv"), "w", newline="") as f: - if rows: - w = csv.DictWriter(f, fieldnames=rows[0].keys()) - w.writeheader() - for r in rows: - w.writerow(r) - - # Score = (closeness to 10.5 GHz) + (BW) + (S11 dip depth) - def score(r): - if r["bw_pct"] == 0: - return -1e9 - f_off = abs(r["f_res_GHz"] - 10.5) - return -10.0*f_off + r["bw_pct"] - 0.5*r["s11_min_dB"] - rows.sort(key=score, reverse=True) - print() - print("=" * 78) - print("Top 5 sweep results (best score first):") - print(f" {'slot':>6} {'stub':>6} {'f_res':>8} {'S11':>7} {'BW':>7} {'lo–hi':>15}") - for r in rows[:5]: - print(f" {r['slot_L']:>6.2f} {r['stub_L']:>6.2f} " - f"{r['f_res_GHz']:>6.2f}GHz {r['s11_min_dB']:>5.1f}dB " - f"{r['bw_pct']:>5.2f}% {r['f_lo_GHz']:.2f}–{r['f_hi_GHz']:.2f}GHz") - print("=" * 78) - sys.exit(0) - -# ---- single run ---- -sim_path = os.path.join(OUT_DIR, "single") -freq, s11_dB, zin, vswr, dt = run_case(SLOT_L, FEED_STUB_L, PATCH_L, 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))) - -# Also report at 10.5 GHz exactly so we see the impedance at the operating freq -i_op = int(np.argmin(np.abs(freq - 10.5e9))) - -print() -print("=" * 70) -print(f" Resonance (Im(Zin)=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(" ── 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 Aperture-Coupled Patch v2 — Stack_Hybrid 4-layer " - f"(slot={SLOT_L}mm stub={FEED_STUB_L}mm patch_L={PATCH_L}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") diff --git a/5_Simulations/Antenna/array_factor_adar1000_aeris10.py b/5_Simulations/Antenna/array_factor_adar1000_aeris10.py deleted file mode 100644 index e992acd..0000000 --- a/5_Simulations/Antenna/array_factor_adar1000_aeris10.py +++ /dev/null @@ -1,422 +0,0 @@ -#!/usr/bin/env python3 -# array_factor_adar1000_aeris10.py -# -# Phased-array beam-forming verification using the ADAR1000 firmware's actual -# phase-shifter codes. Combines: -# * The single-row 1x8 series-fed embedded element pattern (from -# edge_fed_row_nf2ff_aeris10_v3.py at 10.520 GHz, cached) for the y-axis -# row pattern. -# * A 16-element x-axis array factor at d = λ/2 = 14.286 mm pitch (matches -# the firmware's `element_spacing = wavelength/2` constant). -# * The firmware phase computation EXACTLY (ADAR1000_Manager.cpp:714-729): -# `calculatePhaseSettings()` only fills 4 phases (one per chip channel), -# and the broadcast loop applies the same 4-phase pattern to all 4 chips. -# * The 7-bit (128-state, 2.8125 deg/step) ADAR1000 phase quantization. -# -# Verifications a radar engineer would run at this stage: -# 1. Beam steering accuracy (commanded vs simulated peak angle). -# 2. Sidelobe and grating-lobe levels at multiple scan angles. -# 3. Scan loss (peak gain vs scan angle). -# 4. Null steering: place a deep null at a chosen angle. -# 5. Compare FIRMWARE behaviour vs CORRECT 16-element progressive phasing -# to expose the per-chip-broadcast bug. -# -# Inputs: -# /tmp/aeris10_edgefed_row_nf2ff_v3/farfield.csv (cached single-row pattern) -# -# Outputs: -# /tmp/aeris10_array_factor/scan_*.png (1D cuts at scan angles) -# /tmp/aeris10_array_factor/scan_loss.png (peak gain vs scan) -# /tmp/aeris10_array_factor/null_steering.png -# /tmp/aeris10_array_factor/firmware_vs_correct.png - -import os -import csv -import numpy as np -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt - -# ============================================================================ -# Constants (match firmware ADAR1000_Manager.cpp:714-729 exactly) -# ============================================================================ -F0 = 10.5e9 -C0 = 3.0e8 -LAMBDA = C0 / F0 # 28.5714 mm -D_X = LAMBDA / 2 # 14.2857 mm — element_spacing in firmware -N_TOTAL = 16 # 4 chips × 4 channels -N_PER_CHIP = 4 -N_CHIPS = 4 - -# ADAR1000 7-bit phase resolution -PHASE_STATES = 128 -PHASE_LSB_DEG = 360.0 / PHASE_STATES # 2.8125 deg/code - -OUT_DIR = "/tmp/aeris10_array_factor" -os.makedirs(OUT_DIR, exist_ok=True) - - -# ============================================================================ -# Embedded element pattern (single-row 1x8 at 10.520 GHz) -# ============================================================================ -def load_single_row_pattern(path="/tmp/aeris10_edgefed_row_nf2ff_v3/farfield.csv"): - """Returns theta_deg, h_pat_lin, e_pat_lin (linear |E|, peak normalised to 1).""" - th, h_dB, e_dB = [], [], [] - with open(path) as f: - r = csv.reader(f); next(r) - for row in r: - th.append(float(row[0])) - h_dB.append(float(row[1])) - e_dB.append(float(row[2])) - th = np.array(th); h_dB = np.array(h_dB); e_dB = np.array(e_dB) - # CSV stores normalised dB rel peak. Convert to linear amplitude. - h_lin = 10 ** (h_dB / 20.0) - e_lin = 10 ** (e_dB / 20.0) - return th, h_lin, e_lin - - -# ============================================================================ -# Phase code generators -# ============================================================================ -def firmware_phase_codes(angle_deg): - """Replicate ADAR1000Manager::calculatePhaseSettings() + the broadcast - loop in setBeamAngle(): 4 phases computed, same pattern applied to all - 4 chips. Returns 16 ADAR1000 phase codes (uint8 values 0..127).""" - angle_rad = np.deg2rad(angle_deg) - # firmware: phase_shift = 2π·d·sin(θ)/λ - phase_shift = (2 * np.pi * D_X * np.sin(angle_rad)) / LAMBDA - codes_4 = np.zeros(N_PER_CHIP, dtype=int) - for i in range(N_PER_CHIP): - ph = (i * phase_shift) % (2 * np.pi) - codes_4[i] = int(round(ph / (2 * np.pi) * PHASE_STATES)) % PHASE_STATES - # Broadcast: same 4-element pattern repeated to all 4 chips - return np.tile(codes_4, N_CHIPS) - - -def correct_phase_codes(angle_deg): - """Proper 16-element progressive phase shift (what the firmware should do).""" - angle_rad = np.deg2rad(angle_deg) - phase_shift = (2 * np.pi * D_X * np.sin(angle_rad)) / LAMBDA - codes = np.zeros(N_TOTAL, dtype=int) - for n in range(N_TOTAL): - ph = (n * phase_shift) % (2 * np.pi) - codes[n] = int(round(ph / (2 * np.pi) * PHASE_STATES)) % PHASE_STATES - return codes - - -def codes_to_radians(codes): - return np.asarray(codes) * (2 * np.pi / PHASE_STATES) - - -# ============================================================================ -# Array factor at the φ=0 (H-plane) cut, for an arbitrary 16-element phase set -# ============================================================================ -def array_factor_hplane(theta_deg_arr, phase_codes, amplitudes=None): - """Compute |AF(θ)| at φ=0 (H-plane). x_n = n*d. Element 0 at x=0, element - 15 at x=15·d. AF(θ) = Σ a_n · exp(j·k·x_n·sin(θ)) · exp(j·φ_n).""" - if amplitudes is None: - amplitudes = np.ones(len(phase_codes)) - k = 2 * np.pi / LAMBDA - th_rad = np.deg2rad(theta_deg_arr) - phases_rad = codes_to_radians(phase_codes) - af = np.zeros(len(th_rad), dtype=complex) - for n in range(len(phase_codes)): - xn = n * D_X - af += amplitudes[n] * np.exp(1j * (k * xn * np.sin(th_rad) + phases_rad[n])) - return np.abs(af) - - -def total_pattern_dB(theta_deg_arr, phase_codes, h_pat_lin): - """Single-row pattern × |AF| → normalised dB rel peak.""" - af = array_factor_hplane(theta_deg_arr, phase_codes) - pat_lin = af * h_pat_lin - pat_dB = 20 * np.log10(pat_lin / np.max(pat_lin) + 1e-30) - return pat_dB - - -def find_peak(theta_deg_arr, pat_dB): - i = int(np.argmax(pat_dB)) - return theta_deg_arr[i], pat_dB[i] - - -def find_main_lobe(theta_deg_arr, pat_dB, search_window=None): - """Find the deepest dip / peak in a window. Returns (peak_angle, peak_dB, - bw3dB, sll_dB, sll_angle).""" - if search_window is None: - mask = np.ones(len(theta_deg_arr), dtype=bool) - else: - lo, hi = search_window - mask = (theta_deg_arr >= lo) & (theta_deg_arr <= hi) - idx = np.where(mask)[0] - i_pk_local = idx[int(np.argmax(pat_dB[idx]))] - peak_angle = theta_deg_arr[i_pk_local] - peak_dB = pat_dB[i_pk_local] - # 3 dB beamwidth around peak - half = peak_dB - 3.0 - lo_i = i_pk_local - while lo_i > 0 and pat_dB[lo_i] > half: - lo_i -= 1 - hi_i = i_pk_local - while hi_i < len(pat_dB) - 1 and pat_dB[hi_i] > half: - hi_i += 1 - bw3 = theta_deg_arr[hi_i] - theta_deg_arr[lo_i] - # First null walk - null_lo, null_hi = lo_i, hi_i - while null_lo > 0 and pat_dB[null_lo - 1] < pat_dB[null_lo]: - null_lo -= 1 - while null_hi < len(pat_dB) - 1 and pat_dB[null_hi + 1] < pat_dB[null_hi]: - null_hi += 1 - # Sidelobes outside the null-bracketed main lobe - side_mask = np.ones(len(pat_dB), dtype=bool) - side_mask[null_lo:null_hi + 1] = False - if side_mask.any(): - i_sll = int(np.argmax(np.where(side_mask, pat_dB, -100))) - sll_dB = pat_dB[i_sll] - peak_dB - sll_angle = theta_deg_arr[i_sll] - else: - sll_dB, sll_angle = -np.inf, np.nan - return peak_angle, peak_dB, bw3, sll_dB, sll_angle - - -# ============================================================================ -# Verifications -# ============================================================================ -def main(): - theta_deg, h_pat_lin, e_pat_lin = load_single_row_pattern() - print(f"[load] embedded element pattern: {len(theta_deg)} samples, " - f"theta {theta_deg.min():.0f}..{theta_deg.max():.0f}°") - print(f"[const] λ={LAMBDA*1e3:.3f} mm, d=λ/2={D_X*1e3:.3f} mm, " - f"N={N_TOTAL} (4 chips × 4 ch)") - print(f"[const] phase LSB = {PHASE_LSB_DEG:.4f} deg/code (7-bit)") - print() - - # ------------------------------------------------------------------ - # 1) Steering accuracy at multiple commanded angles - # ------------------------------------------------------------------ - angles_to_test = [0, 5, 10, 15, 20, 30, 45] - print("=" * 90) - print(" TEST 1: Beam steering accuracy — firmware vs correct (φ=0 H-plane)") - print("=" * 90) - print(f"{'cmd':>5} | {'firmware':>40} | {'correct':>40}") - print(f"{'deg':>5} | {'peak deg':>10} {'BW3':>6} {'SLL dB':>7} {'SLL deg':>8} | " - f"{'peak deg':>10} {'BW3':>6} {'SLL dB':>7} {'SLL deg':>8}") - print("-" * 90) - - rows_for_csv = [] - for ang in angles_to_test: - codes_fw = firmware_phase_codes(ang) - codes_co = correct_phase_codes(ang) - pat_fw = total_pattern_dB(theta_deg, codes_fw, h_pat_lin) - pat_co = total_pattern_dB(theta_deg, codes_co, h_pat_lin) - pk_fw = find_main_lobe(theta_deg, pat_fw) - pk_co = find_main_lobe(theta_deg, pat_co) - print(f"{ang:>5} | {pk_fw[0]:>+10.1f} {pk_fw[2]:>5.1f}° {pk_fw[3]:>+6.1f} " - f"{pk_fw[4]:>+7.1f}° | {pk_co[0]:>+10.1f} {pk_co[2]:>5.1f}° " - f"{pk_co[3]:>+6.1f} {pk_co[4]:>+7.1f}°") - rows_for_csv.append((ang, pk_fw[0], pk_fw[2], pk_fw[3], pk_fw[4], - pk_co[0], pk_co[2], pk_co[3], pk_co[4])) - print() - - with open(os.path.join(OUT_DIR, "steering_table.csv"), "w", newline="") as f: - w = csv.writer(f) - w.writerow(["cmd_deg", "fw_peak_deg", "fw_bw3", "fw_sll_dB", "fw_sll_deg", - "co_peak_deg", "co_bw3", "co_sll_dB", "co_sll_deg"]) - for r in rows_for_csv: - w.writerow(r) - print(f"[out] {OUT_DIR}/steering_table.csv") - - # ------------------------------------------------------------------ - # 2) Side-by-side patterns at scan angles 0°, 15°, 30°, 45° - # ------------------------------------------------------------------ - show_angles = [0, 15, 30, 45] - fig, axes = plt.subplots(len(show_angles), 1, figsize=(11, 3.3*len(show_angles)), - sharex=True) - for ax, ang in zip(axes, show_angles): - codes_fw = firmware_phase_codes(ang) - codes_co = correct_phase_codes(ang) - pat_fw = total_pattern_dB(theta_deg, codes_fw, h_pat_lin) - pat_co = total_pattern_dB(theta_deg, codes_co, h_pat_lin) - ax.plot(theta_deg, pat_co, "g-", lw=1.4, label="correct 16-elem (gold)") - ax.plot(theta_deg, pat_fw, "r-", lw=1.4, label="firmware (4-elem broadcast)") - ax.axvline(ang, color="k", ls=":", lw=0.8, label=f"commanded θ={ang}°") - ax.axvline(-ang, color="grey", ls=":", lw=0.6, alpha=0.5, - label=f"-cmd θ={-ang}° (sign-flip)") - ax.set_xlim(-90, 90) - ax.set_ylim(-40, 2) - ax.set_ylabel("Pattern (dB)") - ax.set_title(f"setBeamAngle({ang}°) — H-plane (φ=0, x-scan) at 10.520 GHz") - ax.grid(True, alpha=0.3) - ax.legend(loc="lower right", fontsize=8) - axes[-1].set_xlabel("θ (deg)") - fig.tight_layout() - fig.savefig(os.path.join(OUT_DIR, "firmware_vs_correct.png"), dpi=140) - plt.close(fig) - print(f"[out] {OUT_DIR}/firmware_vs_correct.png") - - # ------------------------------------------------------------------ - # 3) Scan loss curve (peak gain vs commanded angle, both fw and correct) - # ------------------------------------------------------------------ - scan_angles = np.arange(-60, 61, 2) - peak_dB_fw = [] - peak_dB_co = [] - actual_peak_fw = [] - actual_peak_co = [] - # Reference broadside peak for absolute scan-loss — peak in dB rel peak is 0; - # we want amplitude relative to broadside, so compute |E|² without normalisation. - - def total_amp_lin(theta_deg_arr, phase_codes): - af = array_factor_hplane(theta_deg_arr, phase_codes) - return af * h_pat_lin - - ref_lin = total_amp_lin(theta_deg, np.zeros(N_TOTAL, dtype=int)) - ref_peak = float(np.max(ref_lin)) - for ang in scan_angles: - codes_fw = firmware_phase_codes(ang) - codes_co = correct_phase_codes(ang) - amp_fw = total_amp_lin(theta_deg, codes_fw) - amp_co = total_amp_lin(theta_deg, codes_co) - peak_dB_fw.append(20*np.log10(np.max(amp_fw)/ref_peak + 1e-30)) - peak_dB_co.append(20*np.log10(np.max(amp_co)/ref_peak + 1e-30)) - i_fw = int(np.argmax(amp_fw)) - i_co = int(np.argmax(amp_co)) - actual_peak_fw.append(theta_deg[i_fw]) - actual_peak_co.append(theta_deg[i_co]) - - fig, axes = plt.subplots(1, 2, figsize=(13, 4.5)) - ax = axes[0] - ax.plot(scan_angles, peak_dB_co, "g-", lw=1.6, label="correct 16-elem") - ax.plot(scan_angles, peak_dB_fw, "r-", lw=1.6, label="firmware (4-elem broadcast)") - # Theoretical scan loss = cos(θ) (single-element factor) → in dB: 20·log10(cos) - th_th = np.linspace(-60, 60, 121) - ax.plot(th_th, 20*np.log10(np.cos(np.deg2rad(th_th))), - "k--", lw=1.0, alpha=0.6, label="cos(θ) ideal scan loss") - ax.set_xlabel("Commanded scan angle (deg)") - ax.set_ylabel("Peak gain rel broadside (dB)") - ax.set_title("Scan loss vs commanded angle") - ax.set_xlim(-60, 60) - ax.set_ylim(-25, 2) - ax.grid(True, alpha=0.3) - ax.legend() - - ax = axes[1] - ax.plot(scan_angles, actual_peak_co, "g-", lw=1.6, label="correct 16-elem") - ax.plot(scan_angles, actual_peak_fw, "r-", lw=1.6, label="firmware") - ax.plot(scan_angles, scan_angles, "k--", lw=1.0, alpha=0.6, - label="ideal (peak = cmd)") - ax.plot(scan_angles, -scan_angles, "k:", lw=1.0, alpha=0.4, - label="sign-flipped (peak = -cmd)") - ax.set_xlabel("Commanded scan angle (deg)") - ax.set_ylabel("Actual peak angle (deg)") - ax.set_title("Beam pointing accuracy") - ax.set_xlim(-60, 60) - ax.set_ylim(-90, 90) - ax.grid(True, alpha=0.3) - ax.legend() - fig.tight_layout() - fig.savefig(os.path.join(OUT_DIR, "scan_loss.png"), dpi=140) - plt.close(fig) - print(f"[out] {OUT_DIR}/scan_loss.png") - - # ------------------------------------------------------------------ - # 4) Null-steering: place a null at a chosen angle (LCMV-style minimal) - # ------------------------------------------------------------------ - # Set main beam at θ=0, with an explicit null at θ_null=20° using simple - # phase-only synthesis: subtract a unit-amplitude vector pointed at θ_null. - th_null = 20.0 - k = 2*np.pi/LAMBDA - n = np.arange(N_TOTAL) - a_main = np.exp(1j * 0.0 * n) - a_null = np.exp(1j * k * n * D_X * np.sin(np.deg2rad(th_null))) - # Project: a' = a_main - / * a_null - proj = (np.vdot(a_null, a_main) / np.vdot(a_null, a_null)) * a_null - a_steered = a_main - proj - # Convert complex weights to phase codes (drop amplitude variation — - # ADAR1000 phase shifters are constant-amplitude; we keep amplitude=1 and - # use phase only for an honest sim). - phases_null = np.angle(a_steered) % (2*np.pi) - codes_null = np.round(phases_null / (2*np.pi) * PHASE_STATES).astype(int) % PHASE_STATES - pat_null = total_pattern_dB(theta_deg, codes_null, h_pat_lin) - # Reference: no null - pat_bs = total_pattern_dB(theta_deg, np.zeros(N_TOTAL, dtype=int), h_pat_lin) - - # Find depth of null in the steered pattern at θ=20° - i_null = int(np.argmin(np.abs(theta_deg - th_null))) - null_depth = pat_null[i_null] - 0.0 # rel peak - - fig, ax = plt.subplots(figsize=(11, 4.5)) - ax.plot(theta_deg, pat_bs, "k-", lw=1.2, alpha=0.5, label="broadside, no null") - ax.plot(theta_deg, pat_null, "b-", lw=1.6, - label=f"phase-only null @ θ={th_null}°") - ax.axvline(th_null, color="r", ls=":", lw=0.8, label=f"target null θ={th_null}°") - ax.axhline(-30, color="grey", ls="--", lw=0.6, alpha=0.5) - ax.set_xlim(-90, 90) - ax.set_ylim(-50, 2) - ax.set_xlabel("θ (deg)") - ax.set_ylabel("Pattern (dB)") - ax.set_title(f"Null-steering — broadside main beam with null at θ={th_null}° " - f"(achieved depth: {null_depth:.1f} dB rel peak)") - ax.grid(True, alpha=0.3) - ax.legend() - fig.tight_layout() - fig.savefig(os.path.join(OUT_DIR, "null_steering.png"), dpi=140) - plt.close(fig) - print(f"[out] {OUT_DIR}/null_steering.png") - - # ------------------------------------------------------------------ - # 5) Phase-quantization effect (compare unquantized continuous phase - # to 7-bit quantized phase at θ=15°) - # ------------------------------------------------------------------ - ang = 15 - angle_rad = np.deg2rad(ang) - phase_shift = (2 * np.pi * D_X * np.sin(angle_rad)) / LAMBDA - phases_continuous = np.array([(n*phase_shift) % (2*np.pi) for n in range(N_TOTAL)]) - codes_quant = correct_phase_codes(ang) - - def total_pattern_dB_continuous(theta_deg_arr, phases_rad, h_pat_lin): - k = 2*np.pi/LAMBDA - th_rad = np.deg2rad(theta_deg_arr) - af = np.zeros(len(th_rad), dtype=complex) - for n in range(len(phases_rad)): - af += np.exp(1j*(k*n*D_X*np.sin(th_rad) + phases_rad[n])) - amp = np.abs(af) * h_pat_lin - return 20*np.log10(amp / np.max(amp) + 1e-30) - - pat_cont = total_pattern_dB_continuous(theta_deg, phases_continuous, h_pat_lin) - pat_quant = total_pattern_dB(theta_deg, codes_quant, h_pat_lin) - pk_cont = find_main_lobe(theta_deg, pat_cont) - pk_quant = find_main_lobe(theta_deg, pat_quant) - print() - print("=" * 90) - print(" TEST 2: 7-bit phase quantization vs continuous (at cmd 15°)") - print(f" Continuous phase: peak θ={pk_cont[0]:+.1f}°, BW3={pk_cont[2]:.1f}°, " - f"SLL={pk_cont[3]:+.1f} dB") - print(f" Quantized 7-bit : peak θ={pk_quant[0]:+.1f}°, BW3={pk_quant[2]:.1f}°, " - f"SLL={pk_quant[3]:+.1f} dB") - print(f" → quantization adds {pk_cont[3] - pk_quant[3]:+.2f} dB to the SLL " - f"(positive = quantized has worse SLL)") - print("=" * 90) - - # ------------------------------------------------------------------ - # 6) Grating-lobe envelope check - # ------------------------------------------------------------------ - # Theoretical: grating lobes at sin(θ_g) = ±λ/d - sin(θ_0). At d=λ/2, NO - # grating lobes for any scan angle (since |λ/d - sin(θ_0)| ≥ 1 always). - # The firmware's 4-element broadcast effectively makes super-pitch d_super - # = 4d = 2λ → grating lobes at sin(θ_g) = ±λ/(4d) ± sin(θ_0) = ±0.5 ± sin(θ_0). - print() - print(" TEST 3: Grating-lobe geometry") - print(" Element pitch d = λ/2 → no real-space grating lobes at any scan ✓") - print(" Firmware's 4-elem broadcast → super-pitch d_super = 4d = 2λ") - print(" → grating lobes appear at sin(θ_g) = ±0.5 ± sin(θ_0)") - for ang in [0, 15, 30, 45]: - sin0 = np.sin(np.deg2rad(ang)) - gl = [] - for sign in [+1, -1]: - sin_g = sign*0.5 + sin0 # firmware steers to -ang due to sign convention - if abs(sin_g) <= 1: - gl.append(np.rad2deg(np.arcsin(sin_g))) - print(f" cmd {ang:+d}° → grating lobes at: {[f'{g:+.1f}°' for g in gl]}") - - -if __name__ == "__main__": - main() diff --git a/5_Simulations/Antenna/openems_quartz_slotted_wg_10p5GHz.py b/5_Simulations/Antenna/openems_quartz_slotted_wg_10p5GHz.py deleted file mode 100644 index 374902a..0000000 --- a/5_Simulations/Antenna/openems_quartz_slotted_wg_10p5GHz.py +++ /dev/null @@ -1,329 +0,0 @@ -# openems_quartz_slotted_wg_10p5GHz.py -# Slotted rectangular waveguide (quartz-filled, εr=3.8) tuned to 10.5 GHz. -# Builds geometry, meshes (no GetLine calls), sweeps S-params/impedance over 9.5-11.5 GHz, -# computes 3D far-field, and reports estimated max realized gain. - -import os -import tempfile -import numpy as np -import matplotlib.pyplot as plt -import time - -# --- openEMS / CSXCAD bindings --- -from openEMS import openEMS -from openEMS.physical_constants import C0 -try: - from CSXCAD import ContinuousStructure, AppCSXCAD_BIN - HAVE_APP = True -except ImportError: - from CSXCAD import ContinuousStructure - AppCSXCAD_BIN = None - HAVE_APP = False - -#Set PROFILE to "sanity" first; run and check [mesh] cells: stays reasonable. - -#If it's small, move to "balanced"; once happy, go "full". - -#Toggle VIEW_GEOM=True if you want the 3D viewer (requires AppCSXCAD_BIN available). - -# ========================= -# USER SETTINGS / PROFILES -# ========================= -PROFILE = "sanity" # choose: "sanity" | "balanced" | "full" -VIEW_GEOM = False # True => launch AppCSXCAD viewer (if available) -SIMULATE = True # False => skip FDTD (for quick post-proc dev) -THREADS = 0 # 0 => all cores - -# --- profiles tuned for i5-1135G7 + 16 GB --- -profiles = { - "sanity": { - "Nslots": 12, "mesh_res": 0.8, - "air_x": 6.0, "air_y": 20.0, "air_z": 10.0, - "n_theta": 61, "n_phi": 121, "freq_pts": 201, "pml": 6 - }, - "balanced": { - "Nslots": 24, "mesh_res": 0.6, - "air_x": 8.0, "air_y": 30.0, "air_z": 12.0, - "n_theta": 91, "n_phi": 181, "freq_pts": 301, "pml": 8 - }, - "full": { - "Nslots": 32, "mesh_res": 0.5, - "air_x": 10.0, "air_y": 40.0, "air_z": 15.0, - "n_theta": 91, "n_phi": 181, "freq_pts": 401, "pml": 8 - } -} -cfg = profiles[PROFILE] - -# ===================== -# BAND & WAVEGUIDE SPEC -# ===================== -f0 = 10.5e9 -f_span = 2.0e9 -f_start, f_stop = f0 - f_span/2, f0 + f_span/2 -er_quartz = 3.8 # fused silica/quartz - -# Array constraint (λ0/2 pitch, 1 mm septum) => internal a ~ 13.28 mm -lambda0_mm = (C0/f0) * 1e3 -a = 13.28 # mm (broad wall, internal) -b = 6.50 # mm (narrow wall, internal) <-- comfortable machining - -# Slot starters (tune later for taper) -slot_w = 0.60 # mm across x - -# --- guide wavelength at 10.5 GHz (TE10) --- -fc10 = (C0/(2.0*np.sqrt(er_quartz))) * (1.0/(a*1e-3)) # Hz -lambda_d = (C0/f0) / np.sqrt(er_quartz) # m -lambda_g = lambda_d / np.sqrt(1.0 - (fc10/f0)**2) # m -lambda_g_mm = lambda_g * 1e3 - -# --- slots geometry (from λg) --- -slot_s = 0.5*lambda_g_mm -slot_L = 0.47*lambda_g_mm -margin = 0.25*lambda_g_mm - -# =================================== -# FDTD / CSX / MESH (explicit lines) -# =================================== -unit_mm = 1e-3 -Sim_Path = os.path.join(tempfile.gettempdir(), f'openems_quartz_slotted_wg_{PROFILE}') - -FDTD = openEMS(NrTS=int(6e5), EndCriteria=1e-5) -FDTD.SetGaussExcite(0.5*(f_start+f_stop), 0.5*(f_stop-f_start)) -FDTD.SetBoundaryCond([f'PML_{cfg["pml"]}']*6) -FDTD.SetOverSampling(4) -FDTD.SetTimeStepFactor(0.95) - -CSX = ContinuousStructure() -FDTD.SetCSX(CSX) -mesh = CSX.GetGrid() -mesh.SetDeltaUnit(unit_mm) - -# Pads & extent -t_metal = 0.8 # mm metal wall thickness -air_x = cfg["air_x"] -air_y = cfg["air_y"] -air_z = cfg["air_z"] -mesh_res = cfg["mesh_res"] - -# Length from Nslots -Nslots = cfg["Nslots"] -guide_length_mm = margin + (Nslots-1)*slot_s + margin - -# Simulation extents (mm) -x_min, x_max = -air_x, a + air_x -y_min, y_max = -5.0, b + t_metal + air_y -z_min, z_max = -air_z, guide_length_mm + air_z - -# Slot centers and edges (mm) -z_centers = margin + np.arange(Nslots)*slot_s -delta0 = 0.90 # mm offset from centerline (± alternated) -x_centers = (a/2.0) + np.array([+delta0 if i%2==0 else -delta0 for i in range(Nslots)]) - -x_edges = np.concatenate([x_centers - slot_w/2.0, x_centers + slot_w/2.0]) -z_edges = np.concatenate([z_centers - slot_L/2.0, z_centers + slot_L/2.0]) - -# Mesh lines: explicit (NO GetLine calls) -x_lines = sorted({x_min, -t_metal, 0.0, a, a + t_metal, x_max, *list(x_edges)}) -y_lines = [y_min, 0.0, b, b+t_metal, y_max] -z_lines = sorted({z_min, 0.0, guide_length_mm, z_max, *list(z_edges)}) - -mesh.AddLine('x', x_lines) -mesh.AddLine('y', y_lines) -mesh.AddLine('z', z_lines) - -# Print complexity and rough memory (to help stay inside 16 GB) -Nx, Ny, Nz = len(x_lines)-1, len(y_lines)-1, len(z_lines)-1 -Ncells = Nx*Ny*Nz -mem_fields_bytes = Ncells * 6 * 8 # rough ~ (Ex,Ey,Ez,Hx,Hy,Hz) doubles -dx_min = min(np.diff(x_lines)) -dy_min = min(np.diff(y_lines)) -dz_min = min(np.diff(z_lines)) - -# Optional smoothing to limit max cell size -mesh.SmoothMeshLines('all', mesh_res, ratio=1.4) - -# ================= -# MATERIALS & SOLIDS -# ================= -pec = CSX.AddMetal('PEC') -quartzM = CSX.AddMaterial('QUARTZ') -quartzM.SetMaterialProperty(epsilon=er_quartz) -airM = CSX.AddMaterial('AIR') - -# Quartz full block -quartzM.AddBox([0, 0, 0], [a, b, guide_length_mm]) - -# PEC tube walls -pec.AddBox([-t_metal, 0, 0], [0, b, guide_length_mm]) # left -pec.AddBox([a, 0, 0], [a+t_metal,b, guide_length_mm]) # right -pec.AddBox([-t_metal,-t_metal,0],[a+t_metal,0, guide_length_mm]) # bottom -pec.AddBox( - [-t_metal, b, 0], [a + t_metal, b + t_metal, guide_length_mm] -) # top (slots will pierce) - -# Slots (AIR) overriding top metal -for zc, xc in zip(z_centers, x_centers, strict=False): - x1, x2 = xc - slot_w/2.0, xc + slot_w/2.0 - z1, z2 = zc - slot_L/2.0, zc + slot_L/2.0 - prim = airM.AddBox([x1, b, z1], [x2, b+t_metal, z2]) - prim.SetPriority(10) # ensure cut - -# ========= -# WG PORTS -# ========= -port_thick = max(4*mesh_res, 2.0) # mm -p1_start = [0, 0, max(0.5, 10*mesh_res)] -p1_stop = [a, b, p1_start[2] + port_thick] -FDTD.AddRectWaveGuidePort(port_nr=1, start=p1_start, stop=p1_stop, - p_dir='z', a=a*unit_mm, b=b*unit_mm, mode_name='TE10', excite=1) - -p2_stop = [a, b, guide_length_mm - max(0.5, 10*mesh_res)] -p2_start = [0, 0, p2_stop[2] - port_thick] -FDTD.AddRectWaveGuidePort(port_nr=2, start=p2_start, stop=p2_stop, - p_dir='z', a=a*unit_mm, b=b*unit_mm, mode_name='TE10', excite=0) - -# ========= -# NF2FF BOX -# ========= -def create_nf2ff(FDTD_obj, name, start, stop, frequency): - try: - return FDTD_obj.CreateNF2FFBox(name=name, start=start, stop=stop, frequency=frequency) - except AttributeError: - return FDTD_obj.AddNF2FFBox(name=name, start=start, stop=stop, frequency=frequency) - -nf2ff = create_nf2ff( - FDTD, - name='nf2ff', - start=[x_min+1.0, y_min+1.0, z_min+1.0], - stop =[x_max-1.0, y_max-1.0, z_max-1.0], - frequency=[f0] -) - -# ========== -# VIEW GEOM -# ========== -if VIEW_GEOM and HAVE_APP and AppCSXCAD_BIN: - os.makedirs(Sim_Path, exist_ok=True) - csx_xml = os.path.join(Sim_Path, f'quartz_slotted_wg_{PROFILE}.xml') - CSX.Write2XML(csx_xml) - os.system(f'"{AppCSXCAD_BIN}" "{csx_xml}"') - -# ... right before the FDTD run: -t0 = time.time() -FDTD.Run(Sim_Path, cleanup=True, verbose=2, numThreads=THREADS) -t1 = time.time() - -# ... right before NF2FF (far-field): -t2 = time.time() -try: - res = nf2ff.CalcNF2FF(Sim_Path, [f0], theta, phi) # noqa: F821 -except AttributeError: - res = FDTD.CalcNF2FF(nf2ff, Sim_Path, [f0], theta, phi) # noqa: F821 -t3 = time.time() - -# ... S-parameters postproc timing (optional): -t4 = time.time() -for p in ports: # noqa: F821 - p.CalcPort(Sim_Path, freq) # noqa: F821 -t5 = time.time() - - -# ======= -# RUN FDTD -# ======= -if SIMULATE: - FDTD.Run(Sim_Path, cleanup=True, verbose=2, numThreads=THREADS) - -freq = np.linspace(f_start, f_stop, profiles[PROFILE]["freq_pts"]) -ports = list(FDTD.ports) # Port 1 & 2 in creation order -for p in ports: - p.CalcPort(Sim_Path, freq) - -S11 = ports[0].uf_ref / ports[0].uf_inc -S21 = ports[1].uf_ref / ports[0].uf_inc -Zin = ports[0].uf_tot / ports[0].if_tot - -plt.figure(figsize=(7.6,4.6)) -plt.plot(freq*1e-9, 20*np.log10(np.abs(S11)), lw=2, label='|S11|') -plt.plot(freq*1e-9, 20*np.log10(np.abs(S21)), lw=2, ls='--', label='|S21|') -plt.grid(True) -plt.legend() -plt.xlabel('Frequency (GHz)') -plt.ylabel('Magnitude (dB)') -plt.title(f'S-Parameters (profile: {PROFILE})') - -plt.figure(figsize=(7.6,4.6)) -plt.plot(freq*1e-9, np.real(Zin), lw=2, label='Re{Zin}') -plt.plot(freq*1e-9, np.imag(Zin), lw=2, ls='--', label='Im{Zin}') -plt.grid(True) -plt.legend() -plt.xlabel('Frequency (GHz)') -plt.ylabel('Ohms') -plt.title('Input Impedance (Port 1)') - -# ========================== -# POST: 3D FAR-FIELD / GAIN -# ========================== -n_theta, n_phi = cfg["n_theta"], cfg["n_phi"] -theta = np.linspace(0, np.pi, n_theta) -phi = np.linspace(0, 2*np.pi, n_phi) - -try: - res = nf2ff.CalcNF2FF(Sim_Path, [f0], theta, phi) -except AttributeError: - res = FDTD.CalcNF2FF(nf2ff, Sim_Path, [f0], theta, phi) - -idx_f0 = np.argmin(np.abs(freq - f0)) -Dmax_lin = float(res.Dmax[0]) -mismatch = 1.0 - np.abs(S11[idx_f0])**2 -Gmax_lin = Dmax_lin * float(mismatch) -Gmax_dBi = 10*np.log10(Gmax_lin) - - -# Normalized 3D pattern -E = np.squeeze(res.E_norm) # [th, ph] -E = E / np.max(E) -TH, PH = np.meshgrid(theta, phi, indexing='ij') -R = E -X = R * np.sin(TH) * np.cos(PH) -Y = R * np.sin(TH) * np.sin(PH) -Z = R * np.cos(TH) - -fig = plt.figure(figsize=(7.2,6.2)) -ax = fig.add_subplot(111, projection='3d') -ax.plot_surface(X, Y, Z, rstride=2, cstride=2, linewidth=0, antialiased=True, alpha=0.92) -ax.set_title(f'Normalized 3D Pattern @ {f0/1e9:.2f} GHz\n(peak ≈ {Gmax_dBi:.1f} dBi)') -ax.set_box_aspect((1,1,1)) -ax.set_xlabel('x') -ax.set_ylabel('y') -ax.set_zlabel('z') -plt.tight_layout() - -# ========================== -# QUICK 2D GEOMETRY PREVIEW -# ========================== -plt.figure(figsize=(8.4,2.8)) -plt.fill_between( - [0, a], - [0, 0], - [guide_length_mm, guide_length_mm], - color='#dddddd', - alpha=0.5, - step='pre', - label='WG top aperture', -) -for zc, xc in zip(z_centers, x_centers, strict=False): - plt.gca().add_patch(plt.Rectangle((xc - slot_w/2.0, zc - slot_L/2.0), - slot_w, slot_L, fc='#3355ff', ec='k')) -plt.xlim(-2, a + 2) -plt.ylim(-5, guide_length_mm + 5) -plt.gca().invert_yaxis() -plt.xlabel('x (mm)') -plt.ylabel('z (mm)') -plt.title(f'Top-view slot layout (N={Nslots}, profile={PROFILE})') -plt.grid(True) -plt.legend() - - - -plt.show() diff --git a/5_Simulations/Antenna/probe_fed_array_aeris10_v3.py b/5_Simulations/Antenna/probe_fed_array_aeris10_v3.py deleted file mode 100644 index d095378..0000000 --- a/5_Simulations/Antenna/probe_fed_array_aeris10_v3.py +++ /dev/null @@ -1,384 +0,0 @@ -#!/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 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))) - -# DRIVEN_PORTS overrides DRIVEN_X/Y — comma/semicolon-separated list of "i,j" -# pairs all excited in-phase with equal amplitude. Models perfect 1:8 corporate -# splitter feeding an 8-patch sub-array. Example: "0,0;1,0;2,0;3,0;0,1;1,1;2,1;3,1" -# is a 4-cols × 2-rows sub-array anchored in the corner. -DRIVEN_PORTS_STR = os.environ.get("DRIVEN_PORTS", "") -if DRIVEN_PORTS_STR: - pairs = [] - for tok in DRIVEN_PORTS_STR.replace(';', ',').split(','): - tok = tok.strip() - if tok: - pairs.append(int(tok)) - if len(pairs) % 2 != 0: - raise ValueError("DRIVEN_PORTS must be even count of integers (i,j pairs)") - DRIVEN_SET = set((pairs[k], pairs[k+1]) for k in range(0, len(pairs), 2)) -else: - DRIVEN_SET = {(DRIVEN_X, DRIVEN_Y)} - -# 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, j) in DRIVEN_SET 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) - for (idx, p) in ports: - p.CalcPort(sim_path, freq) - - # For each excited port: S = uf_ref/uf_inc (active reflection with all - # driven ports excited). For each non-excited port: S = uf_ref / - # where is the incident wave from any one driven port (used as - # reference; all driven ports have equal amplitude so any works). - ref_inc = next(p for (idx, p) in ports if idx in DRIVEN_SET).uf_inc - S = {} - for (idx, p) in ports: - if idx in DRIVEN_SET: - S[idx] = p.uf_ref / p.uf_inc # active S11 - else: - S[idx] = p.uf_ref / ref_inc # coupling out - - 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))) - -N_DRIVEN = len(DRIVEN_SET) - -# Print coupling grid -print() -print("=" * 70) -if N_DRIVEN == 1: - dx, dy = list(DRIVEN_SET)[0] - print(f" {N_X}x{N_Y} probe-fed array — driven port at ({dx},{dy})") -else: - sub_str = ' '.join(f"({i},{j})" for (i, j) in sorted(DRIVEN_SET)) - print(f" {N_X}x{N_Y} probe-fed array — {N_DRIVEN}-patch sub-array driven in-phase") - print(f" Sub-array: {sub_str}") -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| at {F0/1e9:.2f} GHz, dB (driven ports show active S11):") -print() -# Layout grid as visual array -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) - marker = "*" if (i, j) in DRIVEN_SET else " " - row += f"{marker}{dB:>6.1f}" - print(row) -print(" (* = driven port)") -print() - -# For each driven port, report active S11 + Zin + per-port BW -print(f" Active S11 per driven port at {F0/1e9:.2f} GHz:") -S11_at_op = [] -zin_at_op = [] -for (i, j) in sorted(DRIVEN_SET): - s = S[(i, j)] - s_dB_op = 20*np.log10(abs(s[i_op]) + 1e-30) - zin_p = 50.0 * (1 + s) / (1 - s) - print(f" ({i},{j}) : S11 = {s_dB_op:>6.2f} dB Z = {zin_p[i_op].real:5.1f} + j{zin_p[i_op].imag:+5.1f} Ω") - S11_at_op.append(s_dB_op) - zin_at_op.append(zin_p[i_op]) - -if N_DRIVEN > 1: - print() - print(" Sub-array uniformity:") - print(f" S11 min/max/avg : {min(S11_at_op):>6.2f} / {max(S11_at_op):>6.2f} / " - f"{sum(S11_at_op)/N_DRIVEN:>6.2f} dB") - R_vals = [z.real for z in zin_at_op] - X_vals = [z.imag for z in zin_at_op] - print(f" R min/max/avg : {min(R_vals):>5.1f} / {max(R_vals):>5.1f} / " - f"{sum(R_vals)/N_DRIVEN:>5.1f} Ω") - print(f" X min/max/avg : {min(X_vals):+5.1f} / {max(X_vals):+5.1f} / " - f"{sum(X_vals)/N_DRIVEN:+5.1f} Ω") - # Average port (what the ADAR channel "sees" through ideal 1:8 splitter) - Z_avg = sum(zin_at_op) / N_DRIVEN - print(" Z avg (= what ADAR channel sees through ideal 1:8 splitter):") - print(f" Z = {Z_avg.real:.1f} + j{Z_avg.imag:+.1f} Ω, " - f"VSWR = {abs((Z_avg-50)/(Z_avg+50)) and (1+abs((Z_avg-50)/(Z_avg+50)))/(1-abs((Z_avg-50)/(Z_avg+50))):.2f}") - -# Coupling out (top non-driven ports) -nondriven_couplings = [(idx, abs(S[idx][i_op])) for idx in S.keys() - if idx not in DRIVEN_SET] -nondriven_couplings.sort(key=lambda x: -x[1]) -print() -print(f" Top-5 strongest couplings OUT of sub-array at {F0/1e9:.2f} GHz:") -for idx, val in nondriven_couplings[:5]: - dB = 20*np.log10(val + 1e-30) - print(f" ({idx[0]},{idx[1]}) |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) - -# Coupling heatmap at 10.5 GHz -fig, ax = plt.subplots(figsize=(7, 6.5)) -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) -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 {N_X}x{N_Y} probe-fed array — |S| at {F0/1e9:.2f} GHz') -for j in range(N_Y): - for i in range(N_X): - if (i, j) in DRIVEN_SET: - ax.text(i, j, f"DRIVEN\n{grid[j,i]:.1f} dB", ha='center', va='center', - color='red', fontsize=8, fontweight='bold') - else: - ax.text(i, j, f"{grid[j,i]:.1f}\ndB", ha='center', va='center', - color='white', fontsize=7) -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")