QFP58

#!/usr/bin/env python3
"""
QFP-58 openEMS models (FDTD).  Geometry comes from qfp_geom / qfp_layers, so parameter edits propagate.

  python3 qfp_sim.py patch   [--res fine|coarse]   patch + pins + disks + ground, 4 ports @ r_port (100 ohm)
                                                   one excitation; 4-fold symmetry -> full 4x4 S + embedded pattern
  python3 qfp_sim.py network [--hand RHCP|LHCP]    feed board only, SMA port (50) + 4 ports @ r_port (100)
  python3 qfp_sim.py cascade [--hand ...]          combines patch + network results (no FDTD run): S11, CP gain, AR
  python3 qfp_sim.py full    [--hand ...]          complete antenna, single SMA port: ground truth (slow)

Split plane for patch/network: r_port = rp + 2.0 mm on the radial 100-ohm stubs (identical lumped ports on both
sides, so port parasitics largely cancel in the cascade). The network sits behind a solid ground plane, so
coupling between the two sub-models is neglected; the 'full' model checks that assumption.
Frame: z=-1 network copper (feed B.Cu), z=0 ground (feed F.Cu), z=4 patch copper, z=5 coupling disks.
"""
import os, sys, math, argparse, json
import numpy as np
from shapely.geometry import Point
from qfp_geom import P, build, u
from qfp_layers import feed_layers, patch_layers, fracture, outline_feed

_trapz = getattr(np, "trapezoid", None) or getattr(np, "trapz")
C0, EPS0, ETA0 = 299792458.0, 8.8541878128e-12, 376.730313
F0 = P['f0']
Z_H, Z_GND, Z_PATCH, Z_DISK = -P['sub_h'], 0.0, P['gap_air'], P['gap_air'] + P['sub_h']
R_PORT = P['rp'] + P['match_len']      # interface: end of the 100-ohm matching section
PORT_ANG = [k*math.pi/2 for k in range(4)]          # port k at k*90 deg (CCW, top view)

RES = dict(fine=dict(fine=0.10, patch=0.30, pin=0.10, sub=0.25, gap=0.35),
           medium=dict(fine=0.13, patch=0.35, pin=0.12, sub=0.25, gap=0.4),
           coarse=dict(fine=0.20, patch=0.45, pin=0.15, sub=0.34, gap=0.5))


# ----------------------------------------------------------------------------- CSX helpers
def setup(sim_path, nrts, end=1e-4):
    from CSXCAD import ContinuousStructure
    from openEMS import openEMS
    fd = openEMS(NrTS=nrts, EndCriteria=end)
    fd.SetGaussExcite(F0, 1.2e9)
    fd.SetBoundaryCond(['MUR']*6)
    csx = ContinuousStructure(); fd.SetCSX(csx)
    csx.GetGrid().SetDeltaUnit(1e-3)
    os.makedirs(sim_path, exist_ok=True)
    # sim_er/sim_tand: the laminate actually simulated. sub_er/sub_tand stay the design values, because the line
    # widths and arc lengths are synthesized from them (i.e. "boards drawn for 4.4, built on something else").
    er, tand = P.get('sim_er', P['sub_er']), P.get('sim_tand', P['sub_tand'])
    kappa = 2*math.pi*F0*EPS0*er*tand                               # tan(delta) at f0
    M = dict(fr4=csx.AddMaterial('FR4', epsilon=er, kappa=kappa),
             # PEC sheets: openEMS' conducting-sheet model went unstable on this mesh (energy growth after
             # ~3 ns). Conductor loss is ~10 % of the FR4 dielectric loss here; add it analytically if needed.
             cu=csx.AddMetal('cu'),
             pec=csx.AddMetal('pins_vias'))
    return fd, csx, M

def add_geom(prop, geom, z, prio=10):
    for poly in fracture(geom):
        xy = np.array(poly.exterior.coords)[:-1].T
        prop.AddPolygon(points=xy, norm_dir='z', elevation=z, priority=prio)

def add_feed_substrate(M, f):
    xy = np.array(outline_feed(f).exterior.coords)[:-1].T
    M['fr4'].AddLinPoly(points=xy, norm_dir='z', elevation=Z_H, length=P['sub_h'], priority=1)

def add_patch_board(M):
    PL = patch_layers()
    M['fr4'].AddCylinder(start=[0, 0, Z_PATCH], stop=[0, 0, Z_DISK], radius=P['patch_d']/2 + P['patch_margin'], priority=1)
    add_geom(M['cu'], PL['B_Cu'], Z_PATCH)
    add_geom(M['cu'], PL['F_Cu'], Z_DISK)
    for a in PORT_ANG:                                   # 1.0 mm Cu feed pins, B.Cu (z=-1) to disk (z=5)
        p = P['rp']*u(a)
        M['pec'].AddCylinder(start=[p[0], p[1], Z_H], stop=[p[0], p[1], Z_DISK], radius=0.5, priority=20)

def add_vias(M, f):
    for v in f['vias']:
        M['pec'].AddCylinder(start=[v[0], v[1], Z_H], stop=[v[0], v[1], Z_GND], radius=P['via_drill']/2, priority=20)

def add_resistors(csx, f):
    """Wilkinson isolation resistors as openEMS lumped elements. Lumped elements must be axis-aligned, so each
    rotated 0603 is modelled as a 0.25 mm wide resistive strip along the better-aligned axis, spanning exactly
    the copper gap between its two pads (found with shapely). Returns extra mesh lines."""
    from shapely.geometry import LineString
    from qfp_layers import rect
    d, lines = P['res_pitch']/2, []
    for ref, c, a, val in f['res']:
        pads = [rect(c + k*d*u(a), P['res_pad'][0], P['res_pad'][1], a) for k in (-1, 1)]
        best = None
        for ax in (0, 1):
            e = np.eye(2)[ax]
            ln = LineString([tuple(c - 3*e), tuple(c + 3*e)])
            hit = [ln.intersection(pd) for pd in pads]
            if any(h.is_empty for h in hit):
                continue
            ends = sorted([max(h.coords, key=lambda q: -abs(np.dot(np.array(q) - c, e))) for h in hit],
                          key=lambda q: np.dot(np.array(q) - c, e))
            gap = abs(np.dot(np.array(ends[1]) - np.array(ends[0]), e))
            if best is None or gap < best[1]:
                best = (ax, gap, ends)
        assert best, f'{ref}: no axis-aligned path between pads'
        ax, gap, (p0, p1) = best
        w = 0.25
        lo, hi = min(p0[ax], p1[ax]), max(p0[ax], p1[ax])
        start, stop = [0, 0, Z_H], [0, 0, Z_H]
        start[ax], stop[ax] = lo, hi
        start[1-ax], stop[1-ax] = c[1-ax] - w/2, c[1-ax] + w/2
        le = csx.AddLumpedElement(ref, ny='xy'[ax], caps=True, R=float(val))
        le.AddBox(start, stop, priority=15)
        lines += [start, stop]
        print(f'{ref}: {val} ohm along {"xy"[ax]}, gap {gap:.2f} mm')
    return lines

def stub_width(f, r0=None):
    r0 = P['rp'] if r0 is None else r0
    return [p['w'] for p in f['prims'] if p['kind'] == 'seg' and abs(np.linalg.norm(p['p0']) - r0) < 1e-6][0]

def lumped_port(fd, nr, pos, w, R, excite, along):
    """vertical lumped port (trace z=-1 -> ground z=0) across a line of width w; `along` = line direction (rad)"""
    if abs(math.cos(along)) > 0.9:     # line along x -> sheet in yz
        s, e = [pos[0], pos[1] - w/2, Z_H], [pos[0], pos[1] + w/2, Z_GND]
    elif abs(math.sin(along)) > 0.9:
        s, e = [pos[0] - w/2, pos[1], Z_H], [pos[0] + w/2, pos[1], Z_GND]
    else:                              # oblique line: small square box on the line end
        h = w/3
        s, e = [pos[0] - h, pos[1] - h, Z_H], [pos[0] + h, pos[1] + h, Z_GND]
    return fd.AddLumpedPort(nr, R, s, e, 'z', excite, priority=5)


# ----------------------------------------------------------------------------- mesh
def seeds(centres, half, res):
    return np.concatenate([np.arange(c - half, c + half + 1e-9, res) for c in centres])

def merge(important, seed, dmin):
    """keep all `important` lines, add seed lines only if >= dmin from every kept line (avoids tiny cells)"""
    keep = []
    for x in np.unique(np.round(important, 4)):
        if not keep or x - keep[-1] >= dmin/2:
            keep.append(x)
    keep = np.array(keep)
    for x in np.sort(np.round(seed, 4)):
        if np.min(np.abs(keep - x)) >= dmin:
            keep = np.sort(np.append(keep, x))
    return keep

def make_mesh(csx, xy_seed, z_seed, xy_bound, z_bound, extra_xy=(), dmin=0.08):
    lam_air = C0/(F0 + 1.2e9)*1e3
    g = csx.GetGrid()
    imp = np.concatenate([[-xy_bound, xy_bound], extra_xy, -np.array(extra_xy, float)])
    xy = merge(imp, np.concatenate([xy_seed, -np.asarray(xy_seed)]), dmin)    # mirror-symmetric -> 90-deg symmetric
    for ax in 'xy':
        g.AddLine(ax, xy)
        g.SmoothMeshLines(ax, lam_air/20, 1.4)
    g.AddLine('z', merge(np.concatenate([z_bound, [Z_H, Z_GND, Z_PATCH, Z_DISK]]), z_seed, 0.1))
    g.SmoothMeshLines('z', lam_air/20, 1.4)
    n = [len(g.GetLines(a)) for a in 'xyz']
    print(f'mesh {n[0]}x{n[1]}x{n[2]} = {np.prod(n)/1e6:.2f} Mcells, min cell '
          f'{min(np.diff(g.GetLines(a)).min() for a in "xyz"):.3f} mm')
    return g

def nf2ff_box(fd, g, inset=4):
    lx, ly, lz = (g.GetLines(a) for a in 'xyz')
    return fd.CreateNF2FFBox(start=[lx[inset], ly[inset], lz[inset]], stop=[lx[-inset-1], ly[-inset-1], lz[-inset-1]])


# ----------------------------------------------------------------------------- models
def model_patch(sim_path, res, threads, nrts, run=True, ff=True):
    r = RES[res]
    fd, csx, M = setup(sim_path, nrts)
    f, _ = build('RHCP')
    FL = feed_layers(f)
    add_feed_substrate(M, f)
    add_geom(M['cu'], FL['F_Cu'], Z_GND)
    add_patch_board(M)
    wout = stub_width(f)
    cu = M['cu']
    ports = []
    for k, a in enumerate(PORT_ANG):
        pp = P['rp']*u(a)
        cu.AddPolygon(points=np.array(fracture(Point(*pp).buffer(P['pin_pad']/2, 64))[0].exterior.coords)[:-1].T,
                      norm_dir='z', elevation=Z_H, priority=10)
        e = R_PORT*u(a); t = u(a + math.pi/2)
        quad = np.array([pp - wout/2*t, e - wout/2*t, e + wout/2*t, pp + wout/2*t]).T
        cu.AddPolygon(points=quad, norm_dir='z', elevation=Z_H, priority=10)
        ports.append(lumped_port(fd, k + 1, e, wout, P['z_match'], 1.0 if k == 0 else 0.0, a))
    xy_seed = np.concatenate([np.arange(-P['patch_d']/2 - 1.5, P['patch_d']/2 + 1.5 + 1e-9, r['patch']),
                              seeds([-P['rp'], 0.0, P['rp']], 1.6, r['pin'])])
    z_seed = np.concatenate([np.arange(Z_H, Z_GND + 1e-9, r['sub']), np.arange(Z_GND, Z_PATCH + 1e-9, r['gap']),
                             np.arange(Z_PATCH, Z_DISK + 1e-9, r['sub'])])
    g = make_mesh(csx, xy_seed, z_seed, P['board_d']/2 + 17, [Z_H - 17, Z_DISK + 17],
                  extra_xy=[R_PORT, wout/2, P['rp']])
    nf = nf2ff_box(fd, g) if ff else None
    if run:
        fd.Run(sim_path, cleanup=True, numThreads=threads)
    return fd, ports, nf


def model_network(sim_path, hand, res, threads, nrts, run=True, full=False):
    r = RES[res]
    fd, csx, M = setup(sim_path, nrts)
    f, _ = build(hand)
    FL = feed_layers(f)
    add_feed_substrate(M, f)
    add_vias(M, f)
    rlines = add_resistors(csx, f)
    ports = []
    s = f['sma']
    sma_end = s['sig_c'] + (P['sma_sig_len']/2 - 0.5)*s['u']
    w_in = f['prims'][-1]['w']
    if full:
        add_geom(M['cu'], FL['B_Cu'], Z_H)
        add_geom(M['cu'], FL['F_Cu'], Z_GND)
        add_patch_board(M)
        ports.append(lumped_port(fd, 1, sma_end, w_in, P['z_in'], 1.0, math.atan2(*s['u'][::-1])))
    else:
        net = FL['B_Cu'].difference(Point(0, 0).buffer(R_PORT, 256))
        add_geom(M['cu'], net, Z_H)
        solid = outline_feed(f).buffer(-P['edge_pullback'])          # antipads belong to the patch sub-model
        add_geom(M['cu'], solid, Z_GND)
        ports.append(lumped_port(fd, 1, sma_end, w_in, P['z_in'], 1.0, math.atan2(*s['u'][::-1])))
        for k, a in enumerate(PORT_ANG):
            ports.append(lumped_port(fd, k + 2, R_PORT*u(a), stub_width(f, R_PORT), P['z_out'], 0.0, a))
    reach = 24.0
    xy_seed = np.arange(-reach, reach + 1e-9, r['fine'])
    if full:
        xy_seed = np.concatenate([xy_seed, seeds([-P['rp'], 0.0, P['rp']], 1.6, r['pin'])])
    zt = Z_DISK + 17 if full else Z_GND + 8
    z_seed = np.arange(Z_H, Z_GND + 1e-9, r['sub'])
    if full:
        z_seed = np.concatenate([z_seed, np.arange(Z_GND, Z_PATCH + 1e-9, r['gap']), np.arange(Z_PATCH, Z_DISK + 1e-9, r['sub'])])
    g = make_mesh(csx, xy_seed, z_seed, P['board_d']/2 + (17 if full else 10), [Z_H - (17 if full else 10), zt],
                  extra_xy=[R_PORT, stub_width(f, R_PORT)/2, stub_width(f)/2]
                  + [float(q[0]) for q in rlines] + [float(q[1]) for q in rlines])
    nf = nf2ff_box(fd, g) if full else None
    if run:
        fd.Run(sim_path, cleanup=True, numThreads=threads)
    return fd, ports, nf


# ----------------------------------------------------------------------------- post-processing
FREQ = np.linspace(5.0e9, 6.6e9, 161)
TH = np.arange(0, 181, 2.0)
PH = np.arange(0, 360, 5.0)          # 90 deg must be a multiple of the step (pattern rotation by index shift)

def cp(Et, Ep):
    ER, EL = (Et + 1j*Ep)/np.sqrt(2), (Et - 1j*Ep)/np.sqrt(2)     # IEEE, exp(+jwt), outgoing wave
    ar = (np.abs(ER) + np.abs(EL))/np.maximum(np.abs(np.abs(ER) - np.abs(EL)), 1e-30)
    return ER, EL, 20*np.log10(ar)

def gain_dbi(E, p_inc, r=1.0):
    return 10*np.log10(4*np.pi*np.abs(E)**2*r**2/(2*ETA0)/p_inc + 1e-30)

def prad(Et, Ep, r=1.0):
    U = (np.abs(Et)**2 + np.abs(Ep)**2)*r**2/(2*ETA0)          # phi grid is periodic [0, 360)
    return _trapz(U.sum(axis=1)*np.radians(PH[1] - PH[0])*np.sin(np.radians(TH)), np.radians(TH))

def rotate(E, k):
    return np.roll(E, int(round(90/(PH[1] - PH[0])))*k, axis=1)

def plot_3d(path, g, title, lim=(-10, 10)):
    """3D pattern: radius and colour = gain in dBi, clipped to `lim` (same convention as the TFP page)."""
    import matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
    th, ph = np.radians(TH), np.radians(np.append(PH, 360))
    T, F = np.meshgrid(th, ph, indexing='ij')
    gc = np.clip(np.concatenate([g, g[:, :1]], axis=1), *lim)
    r = gc - lim[0]
    X, Y, Z = r*np.sin(T)*np.cos(F), r*np.sin(T)*np.sin(F), r*np.cos(T)
    norm = matplotlib.colors.Normalize(*lim)
    fig = plt.figure(figsize=(7.2, 7.6))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(X, Y, Z, facecolors=plt.cm.plasma(norm(gc)), rstride=1, cstride=1,
                    linewidth=0.15, edgecolor=(0, 0, 0, 0.35), antialiased=True, shade=False)
    xy = max(np.abs(X).max(), np.abs(Y).max())                       # equal aspect, fitted to the pattern
    ax.set(xlim=(-xy, xy), ylim=(-xy, xy), zlim=(Z.min(), Z.max()))
    ax.set_box_aspect((2*xy, 2*xy, Z.max() - Z.min()), zoom=1.25)
    ax.view_init(elev=18, azim=-60); ax.set_axis_off()
    fig.colorbar(plt.cm.ScalarMappable(norm=norm, cmap='plasma'), ax=ax, orientation='horizontal',
                 fraction=0.04, pad=0.02, label='realized gain (dBi)')
    ax.set_title(title)
    plt.savefig(path, dpi=110, bbox_inches='tight'); plt.close(fig)

def plot_cuts(path, gco, gx, hand, title, step=15, rlim=(-20, 10)):
    """Vertical cross-sections: elevation +90 (boresight) .. -90 (back), one curve per azimuth, all overlaid."""
    import matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
    el = np.radians(90 - TH)
    fig = plt.figure(figsize=(5.2, 7.6))
    ax = fig.add_subplot(111, projection='polar')
    xh = 'LHCP' if hand == 'RHCP' else 'RHCP'
    for j in range(0, len(PH), int(round(step/(PH[1] - PH[0])))):
        ax.plot(el, np.clip(gco[:, j], *rlim), color='C1', lw=0.8, alpha=0.7, label=f'{hand} (co-pol)' if j == 0 else None)
        ax.plot(el, np.clip(gx[:, j], *rlim), color='C0', lw=0.8, alpha=0.7, label=f'{xh} (x-pol)' if j == 0 else None)
    ax.set_thetamin(-90); ax.set_thetamax(90); ax.set_rlim(*rlim); ax.set_rticks(np.arange(rlim[0], rlim[1] + 1, 5))
    ax.set_xlabel('realized gain (dBi)'); ax.grid(ls=':', alpha=0.6)
    ax.set_title(f'{title}\nelevation (deg), all azimuths every {step}°', fontsize=10)
    ax.legend(loc='upper right', fontsize=8, bbox_to_anchor=(1.12, 1.02))
    plt.savefig(path, dpi=110, bbox_inches='tight'); plt.close(fig)

def report_pattern(tag, Et, Ep, p_inc, outdir, f):
    ER, EL, AR = cp(Et, Ep)
    i90 = int(round(90/(PH[1]-PH[0])))
    gR, gL = gain_dbi(ER, p_inc), gain_dbi(EL, p_inc)
    hand = 'RHCP' if gR[0, 0] > gL[0, 0] else 'LHCP'
    gco = gR if hand == 'RHCP' else gL
    gx = gL if hand == 'RHCP' else gR
    hp = TH[np.argmax(gco[:, 0] < gco[0, 0] - 3)] if np.any(gco[:, 0] < gco[0, 0] - 3) else np.nan
    sel = TH <= (hp if np.isfinite(hp) else 30)
    # -3 dB half-angle per azimuth (linear interpolation), for the beamwidth spread over all cuts
    hps = [np.interp(-3, (gco[:, j] - gco[0, j])[::-1], TH[::-1]) for j in range(len(PH))
           if np.any(gco[:, j] < gco[0, j] - 3)]
    gtot = 10*np.log10(10**(gR/10) + 10**(gL/10))
    res = dict(freq_GHz=f/1e9, co_pol=hand, realized_gain_boresight_dBi=float(gco[0, 0]),
               xpol_gain_boresight_dBi=float(gx[0, 0]),
               AR_boresight_dB=float(AR[0, 0]), HPBW_half_angle_phi0_deg=float(hp),
               HPBW_deg_min_max=[float(2*min(hps)), float(2*max(hps))] if hps else None,
               AR_mean_within_HPBW_dB=float(AR[sel, :].mean()), AR_max_theta30_dB=float(AR[TH <= 30, :].max()),
               front_to_back_dB=float(gco[0, 0] - gtot[-1, :].max()),
               radiated_power_over_incident=float(prad(Et, Ep)/p_inc))
    np.savez_compressed(f'{outdir}/{tag}_farfield_{f/1e9:.2f}GHz.npz', theta_deg=TH, phi_deg=PH,
                        gain_RHCP_dBi=gR, gain_LHCP_dBi=gL, AR_dB=AR)
    try:
        plot_3d(f'{outdir}/{tag}_pattern3d_{f/1e9:.2f}GHz.png', gco,
                f'{tag}: {hand} realized gain @ {f/1e9:.2f} GHz (dBi)')
        plot_cuts(f'{outdir}/{tag}_cuts_{f/1e9:.2f}GHz.png', gco, gx, hand, f'{tag} @ {f/1e9:.2f} GHz')
    except Exception as ex:
        print('3D/cut plot skipped:', ex)
    try:
        import matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
        fig, ax = plt.subplots(1, 2, figsize=(11, 4))
        for j, lab in ((0, 'φ=0°'), (i90, 'φ=90°')):
            ax[0].plot(TH, gco[:, j], label=f'{hand} {lab}')
            ax[0].plot(TH, (gL if hand == 'RHCP' else gR)[:, j], '--', label=f'x-pol {lab}')
            ax[1].plot(TH, AR[:, j], label=lab)
        ax[0].set(xlabel='θ (deg)', ylabel='realized gain (dBi)', ylim=(-30, 12), title=f'{tag} @ {f/1e9:.2f} GHz')
        ax[1].set(xlabel='θ (deg)', ylabel='axial ratio (dB)', ylim=(0, 10), xlim=(0, 90))
        for a in ax: a.grid(alpha=.3); a.legend(fontsize=8)
        plt.tight_layout(); plt.savefig(f'{outdir}/{tag}_pattern_{f/1e9:.2f}GHz.png', dpi=110); plt.close(fig)
    except Exception as ex:
        print('plot skipped:', ex)
    return res

def nf2ff_fields(nf, sim_path, freqs):
    r = nf.CalcNF2FF(sim_path, freqs, TH, PH, center=[0, 0, P['gap_air']/2], read_cached=True)
    return [np.array(r.E_theta[i]) for i in range(len(freqs))], [np.array(r.E_phi[i]) for i in range(len(freqs))]

def post_patch(sim_path, fd, ports, nf):
    ZR = P['z_match']
    for p in ports: p.CalcPort(sim_path, FREQ, ref_impedance=ZR)
    a1 = ports[0].uf_inc
    s = np.array([p.uf_ref/a1 for p in ports])            # s[d] = S_{1+d,1}; ports share Z_ref -> V ratio == wave ratio
    S = np.zeros((len(FREQ), 4, 4), complex)
    for m in range(4):
        for n in range(4):
            S[:, m, n] = s[(m - n) % 4]
    g_act = s[0] - s[2]                                    # sequential excitation (0,-90,-180,-270): adjacent terms cancel
    z_act = ZR*(1 + g_act)/(1 - g_act)
    out = dict(freq=FREQ, s=s, S=S, g_act=g_act, z_act=z_act, z_ref=ZR)
    fl = [F0 - 0.15e9, F0, F0 + 0.15e9]
    summary = dict(f_ReZmax_GHz=float(FREQ[np.argmax(z_act.real)]/1e9), ReZmax=float(z_act.real.max()),
                   f_min_active_S11_GHz=float(FREQ[np.argmin(abs(g_act))]/1e9),
                   active_S11_at_f0_dB=float(20*np.log10(abs(g_act[np.argmin(abs(FREQ - F0))]))),
                   active_S11_vs_z_out_worst_5p65_5p95_dB=float((20*np.log10(abs((z_act - P['z_out'])/(z_act + P['z_out']))))[(FREQ >= 5.65e9) & (FREQ <= 5.95e9)].max()),
                   Z_active_at_f0=str(np.round(z_act[np.argmin(abs(FREQ - F0))], 1)),
                   S21_adjacent_f0_dB=float(20*np.log10(abs(s[1][np.argmin(abs(FREQ - F0))]))),
                   S31_opposite_f0_dB=float(20*np.log10(abs(s[2][np.argmin(abs(FREQ - F0))]))),
                   active_S11_below_m10dB_GHz=[float(FREQ[i]/1e9) for i in np.where(20*np.log10(abs(g_act)) < -10)[0][[0, -1]]]
                   if np.any(20*np.log10(abs(g_act)) < -10) else None)
    np.savez(f'{sim_path}/patch_sparams.npz', **out)
    plot_s(sim_path, 'patch', {'active S11 (seq. feed)': g_act, 'S11 (single port)': s[0],
                               'S21 adjacent': s[1], 'S31 opposite': s[2]})
    if nf is None:
        return summary
    Et, Ep = nf2ff_fields(nf, sim_path, fl)
    k1 = [np.argmin(abs(FREQ - x)) for x in fl]
    a1n = [ports[0].uf_inc[i]/np.sqrt(ZR) for i in k1]          # power-wave amplitude of the excitation
    np.savez(f'{sim_path}/patch_results.npz', **out, pat_f=np.array(fl), Et=np.array(Et)/np.array(a1n)[:, None, None],
             Ep=np.array(Ep)/np.array(a1n)[:, None, None], P_inc=ports[0].P_inc)
    # ideal sequential feed (perfect network): a_k = exp(-j k 90deg)/2  -> total incident power 4*|a|^2/2 = 1/2 ... use unit per port
    pats = []
    for i, f in enumerate(fl):
        E1t, E1p = Et[i]/a1n[i], Ep[i]/a1n[i]
        EtS = sum(np.exp(-1j*k*np.pi/2)*rotate(E1t, k) for k in range(4))
        EpS = sum(np.exp(-1j*k*np.pi/2)*rotate(E1p, k) for k in range(4))
        pats.append(report_pattern('patch_ideal_feed', EtS, EpS, 4*0.5, sim_path, f))
    summary['ideal_feed_patterns'] = pats
    return summary

def post_network(sim_path, ports):
    ports[0].CalcPort(sim_path, FREQ, ref_impedance=P['z_in'])
    for p in ports[1:]: p.CalcPort(sim_path, FREQ, ref_impedance=P['z_out'])
    a0 = ports[0].uf_inc/np.sqrt(P['z_in'])
    col = np.array([ports[0].uf_ref/np.sqrt(P['z_in'])/a0] + [p.uf_ref/np.sqrt(P['z_out'])/a0 for p in ports[1:]])
    np.savez(f'{sim_path}/network_results.npz', freq=FREQ, col=col)
    i0 = np.argmin(abs(FREQ - F0))
    ph = np.degrees(np.angle(col[1:, i0]/col[1, i0]))
    summary = dict(S11_f0_dB=float(20*np.log10(abs(col[0, i0]))),
                   S_k1_f0_dB=[float(20*np.log10(abs(c))) for c in col[1:, i0]],
                   rel_phase_f0_deg=[float(x) for x in ph],
                   network_loss_f0_dB=float(-10*np.log10(np.sum(abs(col[:, i0])**2))))
    plot_s(sim_path, 'network', {'S11 (SMA)': col[0], **{f'S{k+2}1 (pin {k*90}°)': col[k+1] for k in range(4)}})
    return summary

def cascade(patch_dir, net_dir, outdir):
    """network (ports at the interface, ref z_out) driving the patch (ref z_match). First order: the Wilkinson
    outputs are assumed matched (their own reflection S_kk of the network is neglected)."""
    pr, nr = np.load(f'{patch_dir}/patch_results.npz'), np.load(f'{net_dir}/network_results.npz')
    Zp, Zn = float(pr['z_ref']) if 'z_ref' in pr else 100.0, P['z_out']
    S, col = pr['S'], nr['col']                       # S: (f,4,4) ref Zp ; col: (5,f) ref Zn (0 = SMA)
    I4 = np.eye(4)
    Zm = Zp*np.linalg.solve(np.transpose(I4 - S, (0, 2, 1)), np.transpose(I4 + S, (0, 2, 1))).transpose(0, 2, 1)
    Sn = np.linalg.solve((Zm + Zn*I4).transpose(0, 2, 1), (Zm - Zn*I4).transpose(0, 2, 1)).transpose(0, 2, 1)
    a_n = col[1:]                                       # waves launched into the patch (ref Zn) per unit SMA wave
    b_n = np.einsum('fkl,lf->kf', Sn, a_n)              # reflected by the patch (ref Zn)
    gin = col[0] + np.sum(col[1:]*b_n, axis=0)          # back through the network (reciprocity)
    V, Icur = np.sqrt(Zn)*(a_n + b_n), (a_n - b_n)/np.sqrt(Zn)
    a_p = (V + Zp*Icur)/(2*np.sqrt(Zp))                 # incident waves in the patch reference (ref Zp)
    summary = dict(S11_SMA_f0_dB=float(20*np.log10(abs(gin[np.argmin(abs(FREQ - F0))]))),
                   S11_SMA_worst_5p65_5p95_dB=float((20*np.log10(abs(gin)))[(FREQ >= 5.65e9) & (FREQ <= 5.95e9)].max()),
                   S11_SMA_below_m10dB_GHz=[float(FREQ[i]/1e9) for i in np.where(20*np.log10(abs(gin)) < -10)[0][[0, -1]]]
                   if np.any(20*np.log10(abs(gin)) < -10) else None,
                   power_absorbed_in_patch_ports_f0=float(np.sum(abs(a_p[:, np.argmin(abs(FREQ - F0))])**2
                                                               - abs(np.einsum('kl,l->k', S[np.argmin(abs(FREQ - F0))], a_p[:, np.argmin(abs(FREQ - F0))]))**2)))
    pats = []
    for i, f in enumerate(pr['pat_f']):
        fi = np.argmin(abs(FREQ - f))
        EtS = sum(a_p[k, fi]*rotate(pr['Et'][i], k) for k in range(4))
        EpS = sum(a_p[k, fi]*rotate(pr['Ep'][i], k) for k in range(4))
        pats.append(report_pattern('cascade', EtS, EpS, 0.5, outdir, f))
    summary['patterns'] = pats
    plot_s(outdir, 'cascade', {'S11 at SMA': gin, **{f'wave into pin {k*90}°': a_p[k] for k in range(4)}})
    return summary

def post_full(sim_path, ports, nf):
    ports[0].CalcPort(sim_path, FREQ, ref_impedance=P['z_in'])
    g = ports[0].uf_ref/ports[0].uf_inc
    fl = [F0 - 0.15e9, F0, F0 + 0.15e9]
    Et, Ep = nf2ff_fields(nf, sim_path, fl)
    pats = []
    for i, f in enumerate(fl):
        fi = np.argmin(abs(FREQ - f))
        a0 = ports[0].uf_inc[fi]/np.sqrt(P['z_in'])
        pats.append(report_pattern('full', Et[i]/a0, Ep[i]/a0, 0.5, sim_path, f))
    plot_s(sim_path, 'full', {'S11 at SMA': g})
    return dict(S11_f0_dB=float(20*np.log10(abs(g[np.argmin(abs(FREQ - F0))]))), patterns=pats)

def plot_s(outdir, tag, curves):
    try:
        import matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
        fig, ax = plt.subplots(figsize=(7, 4))
        for k, v in curves.items(): ax.plot(FREQ/1e9, 20*np.log10(abs(v)), label=k)
        ax.axvspan(5.65, 5.95, color='g', alpha=.08)
        ax.set(xlabel='f (GHz)', ylabel='dB', ylim=(-40, 1), title=f'QFP-58 {tag}'); ax.grid(alpha=.3); ax.legend(fontsize=8)
        plt.tight_layout(); plt.savefig(f'{outdir}/{tag}_sparams.png', dpi=110); plt.close(fig)
        np.savetxt(f'{outdir}/{tag}_sparams.csv', np.column_stack([FREQ/1e9] + [f(v) for v in curves.values()
                   for f in (lambda x: 20*np.log10(abs(x)), lambda x: np.degrees(np.angle(x)))]),
                   delimiter=',', header='f_GHz,' + ','.join(f'{k} dB,{k} deg' for k in curves), fmt='%.5g')
    except Exception as ex:
        print('plot skipped:', ex)


# ----------------------------------------------------------------------------- main
if __name__ == '__main__':
    ap = argparse.ArgumentParser()
    ap.add_argument('mode', choices=['patch', 'network', 'cascade', 'full'])
    ap.add_argument('--hand', default='RHCP', choices=['RHCP', 'LHCP'])
    ap.add_argument('--res', default='fine', choices=list(RES))
    ap.add_argument('--threads', type=int, default=0)
    ap.add_argument('--nrts', type=int, default=60000)
    ap.add_argument('--out', default='sim_out')
    ap.add_argument('--post-only', action='store_true', help='skip FDTD, re-run post-processing')
    ap.add_argument('--no-ff', action='store_true', help='patch: skip far field (fast tuning runs)')
    ap.add_argument('--set', nargs='*', default=[], help='override parameters, e.g. --set patch_d=25.6 rp=3.2')
    ap.add_argument('--tag', default='')
    ap.add_argument('--patch-dir', default='', help='cascade: patch result dir (default <out>/patch_<res>)')
    a = ap.parse_args()
    for kv in a.set:
        k, v = kv.split('='); P[k] = float(v)
    R_PORT, Z_PATCH, Z_DISK = P['rp'] + P['match_len'], P['gap_air'], P['gap_air'] + P['sub_h']
    base = os.path.abspath(a.out)
    run = not a.post_only
    if a.mode == 'patch':
        d = f'{base}/patch_{a.res}{a.tag}'
        fd, ports, nf = model_patch(d, a.res, a.threads, a.nrts, run, ff=not a.no_ff)
        s = post_patch(d, fd, ports, nf)
    elif a.mode == 'network':
        d = f'{base}/network_{a.hand}_{a.res}{a.tag}'
        fd, ports, _ = model_network(d, a.hand, a.res, a.threads, a.nrts, run)
        s = post_network(d, ports)
    elif a.mode == 'full':
        d = f'{base}/full_{a.hand}_{a.res}{a.tag}'
        fd, ports, nf = model_network(d, a.hand, a.res, a.threads, a.nrts, run, full=True)
        s = post_full(d, ports, nf)
    else:
        d = f'{base}/cascade_{a.hand}_{a.res}{a.tag}'; os.makedirs(d, exist_ok=True)
        s = cascade(a.patch_dir or f'{base}/patch_{a.res}', f'{base}/network_{a.hand}_{a.res}{a.tag}', d)
    json.dump(s, open(f'{d}/summary.json', 'w'), indent=1, default=float)
    print(json.dumps(s, indent=1, default=float))