diff --git a/5_Simulations/Antenna/array_factor_adar1000_aeris10.py b/5_Simulations/Antenna/array_factor_adar1000_aeris10.py
new file mode 100644
index 0000000..3074a11
--- /dev/null
+++ b/5_Simulations/Antenna/array_factor_adar1000_aeris10.py
@@ -0,0 +1,426 @@
+#!/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
+    ref_pat = total_pattern_dB(theta_deg, np.zeros(N_TOTAL, dtype=int), h_pat_lin)
+    # Find peak amplitude (ref) on linear scale by recomputing without normalisation
+    # Use a simpler proxy: peak in dB rel peak is 0; we want amplitude relative to
+    # broadside. Compute |E|² broadside vs scanned 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, a_main>/<a_null, a_null> * 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)
+    phases_quant = codes_to_radians(codes_quant)
+
+    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(f"  Element pitch d = λ/2 → no real-space grating lobes at any scan ✓")
+    print(f"  Firmware's 4-elem broadcast → super-pitch d_super = 4d = 2λ")
+    print(f"  → 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/edge_fed_row_nf2ff_aeris10_v3.py b/5_Simulations/Antenna/edge_fed_row_nf2ff_aeris10_v3.py
new file mode 100644
index 0000000..f042c03
--- /dev/null
+++ b/5_Simulations/Antenna/edge_fed_row_nf2ff_aeris10_v3.py
@@ -0,0 +1,331 @@
+#!/usr/bin/env python3
+# edge_fed_row_nf2ff_aeris10_v3.py
+#
+# Far-field analysis of the 1x8 series-fed row at f=10.520 GHz to confirm
+# whether the operating mode at the TX-centered design point is broadside
+# (peak in +z, θ=0°) or scanned (peak at θ ≠ 0).
+#
+# Same geometry as edge_fed_row_aeris10_v3.py at the verified design point
+# (CONN_LEN=8.15, no inset on patch 0). Adds an nf2ff probe box around the
+# simulation domain so the radiation pattern can be computed post-sim.
+#
+# Coordinate convention (openEMS/standard physics):
+#   z is normal to the board (+z = upward, away from ground plane)
+#   y is the array axis (patches arranged along +y)
+#   x is the patch-W direction
+#   theta = angle from +z (broadside is θ=0°)
+#   phi=0  → xz plane (H-plane cut, perpendicular to array)
+#   phi=90 → yz plane (E-plane cut, ALONG array axis — this is where the
+#            array factor lives; broadside check happens HERE)
+#
+# Run:
+#   cd /tmp && DYLD_LIBRARY_PATH=/Users/ganeshpanth/opt/openEMS/lib \
+#     PROFILE=balanced \
+#     /Users/ganeshpanth/radar_venv/bin/python \
+#     /Users/ganeshpanth/PLFM_RADAR/5_Simulations/Antenna/edge_fed_row_nf2ff_aeris10_v3.py
+
+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
+
+PROFILE = os.environ.get("PROFILE", "balanced")
+# end_dB is the ENERGY ratio threshold reported in the engine's "(-X dB)" line
+# (10·log10 convention). nf2ff needs at least -40 dB (energy ratio 1e-4) for a
+# valid radiated-power integral.
+profiles = {
+    "sanity":   {"mesh_lambda_div": 18, "n_timesteps": 150000, "end_dB": -35},
+    "balanced": {"mesh_lambda_div": 25, "n_timesteps": 300000, "end_dB": -40},
+}
+cfg = profiles[PROFILE]
+
+F0      = 10.5e9
+F_SPAN  = 4.0e9
+F_TX    = float(os.environ.get("F_TX_GHZ", "10.520")) * 1e9
+F_START = F0 - F_SPAN/2
+F_STOP  = F0 + F_SPAN/2
+
+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
+
+# Verified TX-centered design point
+N_PATCHES    = int(os.environ.get("N_PATCHES", "8"))
+PATCH_W      = float(os.environ.get("PATCH_W_MM", "7.854"))
+PATCH_L      = float(os.environ.get("PATCH_L_MM", "6.95"))
+FEED_W       = float(os.environ.get("FEED_W_MM", "1.16"))
+INSET_DEPTH  = float(os.environ.get("INSET_DEPTH_MM", "0.0"))
+INSET_GAP    = float(os.environ.get("INSET_GAP_MM", "0.30"))
+FEED_LEAD_L  = float(os.environ.get("FEED_LEAD_MM", "15.5"))
+CONN_LEN     = float(os.environ.get("CONN_LEN_MM", "8.15"))
+
+PITCH = PATCH_L + CONN_LEN
+
+Y_FEED_BOARD_EDGE = -PATCH_L/2 - FEED_LEAD_L
+Y_LAST_PATCH_TOP  = (N_PATCHES - 1) * PITCH + PATCH_L/2
+
+GND_X_MARGIN = 14.3
+GND_Y_MARGIN = 14.3
+GND_X_HALF = max(PATCH_W/2, FEED_W/2 + INSET_GAP) + GND_X_MARGIN
+GND_Y_NEG  = Y_FEED_BOARD_EDGE - GND_Y_MARGIN
+GND_Y_POS  = Y_LAST_PATCH_TOP + GND_Y_MARGIN
+
+AIR_ABOVE  = 14.3
+AIR_BELOW  = 14.3
+AIR_X_HALF = GND_X_HALF + 8.0
+AIR_Y_NEG  = GND_Y_NEG - 8.0
+AIR_Y_POS  = GND_Y_POS + 8.0
+
+OUT_DIR = "/tmp/aeris10_edgefed_row_nf2ff_v3"
+os.makedirs(OUT_DIR, exist_ok=True)
+
+
+# ============================================================================
+# Build sim
+# ============================================================================
+fdtd = openEMS(NrTS=cfg["n_timesteps"],
+               EndCriteria=10**(cfg["end_dB"]/10.0))
+fdtd.SetGaussExcite(F0, F_SPAN/2.0)
+fdtd.SetBoundaryCond(["MUR"]*6)
+
+CSX = ContinuousStructure()
+fdtd.SetCSX(CSX)
+mesh = CSX.GetGrid()
+mesh.SetDeltaUnit(1e-3)
+
+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")
+
+patch_sub.AddBox([-GND_X_HALF, GND_Y_NEG, Z_GND + T_CU],
+                 [+GND_X_HALF, GND_Y_POS, Z_PATCH], priority=1)
+copper.AddBox([-GND_X_HALF, GND_Y_NEG, Z_GND],
+              [+GND_X_HALF, GND_Y_POS, Z_GND + T_CU], priority=10)
+
+notch_half_w = FEED_W/2 + INSET_GAP
+
+for i in range(N_PATCHES):
+    py0 = i * PITCH - PATCH_L/2
+    py1 = i * PITCH + PATCH_L/2
+    if i == 0 and INSET_DEPTH > 0.001:
+        copper.AddBox([-PATCH_W/2, py0 + INSET_DEPTH, Z_PATCH],
+                      [+PATCH_W/2, py1, Z_PATCH + T_CU], priority=10)
+        copper.AddBox([-PATCH_W/2, py0, Z_PATCH],
+                      [-notch_half_w, py0 + INSET_DEPTH, Z_PATCH + T_CU],
+                      priority=10)
+        copper.AddBox([+notch_half_w, py0, Z_PATCH],
+                      [+PATCH_W/2, py0 + INSET_DEPTH, Z_PATCH + T_CU],
+                      priority=10)
+    else:
+        copper.AddBox([-PATCH_W/2, py0, Z_PATCH],
+                      [+PATCH_W/2, py1, Z_PATCH + T_CU], priority=10)
+
+for i in range(N_PATCHES - 1):
+    cy0 = i * PITCH + PATCH_L/2
+    cy1 = (i + 1) * PITCH - PATCH_L/2
+    copper.AddBox([-FEED_W/2, cy0, Z_PATCH],
+                  [+FEED_W/2, cy1, Z_PATCH + T_CU], priority=10)
+
+feed_y_start = Y_FEED_BOARD_EDGE
+feed_y_end   = (-PATCH_L/2 + INSET_DEPTH) if INSET_DEPTH > 0.001 else -PATCH_L/2
+copper.AddBox([-FEED_W/2, feed_y_start, Z_PATCH],
+              [+FEED_W/2, feed_y_end, Z_PATCH + T_CU], priority=10)
+
+# Mesh
+lambda_min_mm = (C0 / F_STOP) * 1000.0
+res = lambda_min_mm / cfg["mesh_lambda_div"]
+PORT_LEN = 2.0
+
+xlines = [-AIR_X_HALF, -GND_X_HALF, -PATCH_W/2, -notch_half_w, -FEED_W/2,
+          0, +FEED_W/2, +notch_half_w, +PATCH_W/2, +GND_X_HALF, +AIR_X_HALF]
+
+ylines = [AIR_Y_NEG, GND_Y_NEG, feed_y_start]
+for i in range(N_PATCHES):
+    ylines.append(i * PITCH - PATCH_L/2)
+    ylines.append(i * PITCH)
+    ylines.append(i * PITCH + PATCH_L/2)
+ylines += [GND_Y_POS, AIR_Y_POS]
+port_y_lines = list(np.linspace(feed_y_start, feed_y_start + PORT_LEN, 6))
+ylines += port_y_lines
+
+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(sorted(set(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)
+
+port = fdtd.AddMSLPort(1, copper,
+    start=[-FEED_W/2, feed_y_start, Z_GND + T_CU],
+    stop= [+FEED_W/2, feed_y_start + PORT_LEN, Z_PATCH + T_CU],
+    prop_dir='y', exc_dir='z',
+    excite=1.0,
+    FeedShift=0.4, MeasPlaneShift=1.6,
+    Feed_R=50)
+
+# NF2FF box (created AFTER mesh is set up, BEFORE running)
+nf2ff = fdtd.CreateNF2FFBox()
+
+sim_path = os.path.join(OUT_DIR, "sim")
+print(f"[run] N={N_PATCHES} patch={PATCH_W}x{PATCH_L} CONN={CONN_LEN} "
+      f"PROFILE={PROFILE} cells={n_cells:,}")
+t0 = time.time()
+fdtd.Run(sim_path, verbose=0, cleanup=True)
+dt_sim = time.time() - t0
+print(f"[run] sim time {dt_sim:.1f} s")
+
+# Verify S11 dip is where we expect
+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)
+i_tx = int(np.argmin(np.abs(freq - F_TX)))
+print(f"[s11] S11 @ {F_TX/1e9:.3f} GHz = {s11_dB[i_tx]:.2f} dB")
+
+# ============================================================================
+# Far-field
+# ============================================================================
+print(f"[ff] Computing far-field at {F_TX/1e9:.3f} GHz")
+# CalcNF2FF expects theta/phi in DEGREES (it converts internally).
+# NOTE: Prad and Dmax are buggy in this openEMS version (Prad ≈ 0,
+# sign-dependent). E_norm pattern shape is valid — that's enough to identify
+# broadside vs scanned mode. Plot NORMALIZED pattern (peak at 0 dB).
+theta_deg_arr = np.arange(-180, 181, 2.0)
+phi_deg_arr   = np.array([0.0, 90.0])
+y_array_center = ((N_PATCHES - 1) * PITCH) / 2
+center = [0.0, y_array_center, Z_PATCH/2]
+
+t0 = time.time()
+nf2ff_res = nf2ff.CalcNF2FF(sim_path, F_TX, theta_deg_arr, phi_deg_arr,
+                            center=center, verbose=0)
+dt_ff = time.time() - t0
+print(f"[ff] far-field calc time {dt_ff:.1f} s")
+
+# Result shape: E_norm[freq] is [n_theta, n_phi]
+E_norm = np.array(nf2ff_res.E_norm[0])
+
+# Normalised pattern (peak = 0 dB). Pattern shape is what matters for
+# broadside vs scanned identification.
+E_max_h = float(np.max(np.abs(E_norm[:, 0])))
+E_max_e = float(np.max(np.abs(E_norm[:, 1])))
+E_max_global = max(E_max_h, E_max_e)
+
+pat_dB_h_norm = 20*np.log10(np.abs(E_norm[:, 0]) / E_max_global + 1e-30)
+pat_dB_e_norm = 20*np.log10(np.abs(E_norm[:, 1]) / E_max_global + 1e-30)
+
+theta_deg = theta_deg_arr
+i_bs   = int(np.argmin(np.abs(theta_deg)))
+i_pk_h = int(np.argmax(np.abs(E_norm[:, 0])))
+i_pk_e = int(np.argmax(np.abs(E_norm[:, 1])))
+
+# -3 dB beamwidth in E-plane (array axis — main lobe)
+def beamwidth_3dB(theta_deg, pat_dB, i_peak):
+    half_db = pat_dB[i_peak] - 3.0
+    # Walk left
+    lo = i_peak
+    while lo > 0 and pat_dB[lo] > half_db:
+        lo -= 1
+    # Walk right
+    hi = i_peak
+    while hi < len(pat_dB) - 1 and pat_dB[hi] > half_db:
+        hi += 1
+    return theta_deg[hi] - theta_deg[lo]
+
+bw_e = beamwidth_3dB(theta_deg, pat_dB_e_norm, i_pk_e)
+bw_h = beamwidth_3dB(theta_deg, pat_dB_h_norm, i_pk_h)
+
+print()
+print("=" * 78)
+print(f"  Far-field NORMALISED pattern at f = {F_TX/1e9:.3f} GHz")
+print(f"  ── E-plane (φ=90°, along array axis y — array factor lives here) ──")
+print(f"  Broadside (θ=0°) level : {pat_dB_e_norm[i_bs]:.2f} dB (rel peak)")
+print(f"  Peak direction         : θ = {theta_deg[i_pk_e]:+.1f}° (peak = {pat_dB_e_norm[i_pk_e]:.2f} dB)")
+print(f"  -3 dB beamwidth        : {bw_e:.1f}°")
+print(f"  ── H-plane (φ=0°, perpendicular to array) ──")
+print(f"  Broadside (θ=0°) level : {pat_dB_h_norm[i_bs]:.2f} dB (rel peak)")
+print(f"  Peak direction         : θ = {theta_deg[i_pk_h]:+.1f}° (peak = {pat_dB_h_norm[i_pk_h]:.2f} dB)")
+print(f"  -3 dB beamwidth        : {bw_h:.1f}°")
+print()
+if abs(theta_deg[i_pk_e]) < 5.0:
+    print(f"  → BROADSIDE CONFIRMED: E-plane main lobe at θ={theta_deg[i_pk_e]:+.1f}° (within ±5° of normal)")
+else:
+    print(f"  → SCANNED MODE: E-plane main lobe at θ={theta_deg[i_pk_e]:+.1f}° (NOT broadside)")
+print("=" * 78)
+
+# Plot (normalized to peak = 0 dB)
+fig, axes = plt.subplots(1, 2, figsize=(14, 5))
+for ax, pat_dB, title, peak_deg in [
+    (axes[0], pat_dB_h_norm, f"H-plane (φ=0°, xz cut, perp. to array)",
+        theta_deg[i_pk_h]),
+    (axes[1], pat_dB_e_norm, f"E-plane (φ=90°, yz cut, ALONG array — array factor)",
+        theta_deg[i_pk_e]),
+]:
+    ax.plot(theta_deg, pat_dB, "b-", lw=1.6)
+    ax.axvline(0, color="r", ls=":", lw=0.8, label=f"broadside (θ=0°)")
+    if abs(peak_deg) > 1.0:
+        ax.axvline(peak_deg, color="g", ls=":", lw=0.8,
+                   label=f"peak at θ={peak_deg:+.1f}°")
+    ax.set_xlabel("θ (deg)")
+    ax.set_ylabel("Pattern (dB rel peak)")
+    ax.set_title(f"{title}\nat {F_TX/1e9:.3f} GHz")
+    ax.set_xlim(-180, 180)
+    ax.set_ylim(-40, 2)
+    ax.grid(True, alpha=0.3)
+    ax.legend(loc="lower center")
+fig.tight_layout()
+fig.savefig(os.path.join(OUT_DIR, "farfield_cuts.png"), dpi=140)
+plt.close(fig)
+print(f"[out] {OUT_DIR}/farfield_cuts.png")
+
+# Polar plot (top half-sphere only, θ from -90 to +90)
+fig = plt.figure(figsize=(12, 5))
+ax_h = fig.add_subplot(121, projection="polar")
+ax_e = fig.add_subplot(122, projection="polar")
+mask = (theta_deg >= -90) & (theta_deg <= 90)
+for ax, pat_dB, title in [
+    (ax_h, pat_dB_h_norm, "H-plane (φ=0°)"),
+    (ax_e, pat_dB_e_norm, "E-plane (φ=90°, array axis)"),
+]:
+    th_pol = np.deg2rad(theta_deg[mask])
+    g_norm = np.clip(pat_dB[mask] - np.max(pat_dB[mask]), -40, 0)
+    ax.plot(th_pol, g_norm, "b-", lw=1.6)
+    ax.set_theta_zero_location("N")     # θ=0 at top (broadside)
+    ax.set_theta_direction(-1)          # clockwise
+    ax.set_rlim(-40, 0)
+    ax.set_rticks([-30, -20, -10, 0])
+    ax.set_title(f"{title}\n(normalised to peak)")
+fig.tight_layout()
+fig.savefig(os.path.join(OUT_DIR, "farfield_polar.png"), dpi=140)
+plt.close(fig)
+print(f"[out] {OUT_DIR}/farfield_polar.png")
+
+with open(os.path.join(OUT_DIR, "farfield.csv"), "w", newline="") as f:
+    w = csv.writer(f)
+    w.writerow(["theta_deg", "pattern_hplane_dB_norm", "pattern_eplane_dB_norm"])
+    for k in range(len(theta_deg)):
+        w.writerow([theta_deg[k], pat_dB_h_norm[k], pat_dB_e_norm[k]])
+print(f"[out] {OUT_DIR}/farfield.csv")