Adom Fields

//! Noise couplings: plane cavity resonance and crosstalk, each finding an issue with an action.
//! Cavities: a power or switching pour facing the ground plane is a parallel-plate cavity; its modes solve
//! the Helmholtz equation on the real overlap with open (magnetic-wall) edges: the Neumann Laplacian's
//! lowest eigenpairs (Lanczos, fully reorthogonalised, then Jacobi on the small tridiagonal), k² = λ/h²,
//! f = c k / (2π √εr). A mode near a strong harmonic of the switching (or a named clock) rings.
//! Crosstalk: each quiet net next to a switching node or clock: closest approach, run within 1 mm, the side
//! coupling and the overlap on the next layer as capacitance, the coupled step ΔV Cm / (Cm + Cg).
use crate::board::{Board, LAYERS};
use crate::fields::{Fields, Issue, Spec};
use crate::flow::{is_ground, is_switch_name, rail_from_name};
use serde_json::{json, Value};

const C_MM: f64 = 2.99792458e11;   // mm/s
const E0: f64 = 8.854e-15;         // F/mm
fn r2(v: f64) -> f64 { (v * 100.0).round() / 100.0 }

/// The lowest `k` eigenpairs of a symmetric operator given as a matvec on n unknowns (Lanczos + Jacobi).
pub fn lowest_modes(n: usize, k: usize, apply: &dyn Fn(&[f64], &mut [f64])) -> Vec<(f64, Vec<f64>)> {
    if n == 0 { return Vec::new(); }
    let m = n.min(260);
    let mut vs: Vec<Vec<f64>> = Vec::with_capacity(m);
    let (mut alpha, mut beta) = (vec![0.0; m], vec![0.0; m]);
    let mut v: Vec<f64> = (0..n).map(|i| ((i as f64 * 12.9898).sin() * 43758.5453).fract() - 0.5).collect();
    let nv = v.iter().map(|x| x * x).sum::<f64>().sqrt(); for x in v.iter_mut() { *x /= nv; }
    let mut w = vec![0.0; n]; let mut steps = 0;
    for j in 0..m {
        vs.push(v.clone()); steps = j + 1;
        apply(&v, &mut w);
        alpha[j] = w.iter().zip(&v).map(|(a, b)| a * b).sum();
        // full reorthogonalisation keeps the small modes clean
        for _ in 0..2 { for q in &vs { let c: f64 = w.iter().zip(q).map(|(a, b)| a * b).sum(); for (wi, qi) in w.iter_mut().zip(q) { *wi -= c * qi; } } }
        let b = w.iter().map(|x| x * x).sum::<f64>().sqrt();
        if j + 1 < m { beta[j] = b; if b < 1e-10 { break; } v = w.iter().map(|x| x / b).collect(); }
    }
    let mm = steps;
    // the tridiagonal, then Jacobi
    let mut a = vec![vec![0.0; mm]; mm]; let mut z = vec![vec![0.0; mm]; mm];
    for i in 0..mm { a[i][i] = alpha[i]; z[i][i] = 1.0; if i + 1 < mm { a[i][i + 1] = beta[i]; a[i + 1][i] = beta[i]; } }
    for _sweep in 0..60 {
        let mut off = 0.0; for p in 0..mm { for q in p + 1..mm { off += a[p][q] * a[p][q]; } } if off < 1e-22 { break; }
        for p in 0..mm { for q in p + 1..mm {
            if a[p][q].abs() < 1e-14 { continue; }
            let th = (a[q][q] - a[p][p]) / (2.0 * a[p][q]); let t = th.signum() / (th.abs() + (th * th + 1.0).sqrt()); let t = if th == 0.0 { 1.0 } else { t };
            let c = 1.0 / (t * t + 1.0).sqrt(); let s = t * c;
            for r in 0..mm { let (arp, arq) = (a[r][p], a[r][q]); a[r][p] = c * arp - s * arq; a[r][q] = s * arp + c * arq; }
            for r in 0..mm { let (apr, aqr) = (a[p][r], a[q][r]); a[p][r] = c * apr - s * aqr; a[q][r] = s * apr + c * aqr; }
            for r in 0..mm { let (zrp, zrq) = (z[r][p], z[r][q]); z[r][p] = c * zrp - s * zrq; z[r][q] = s * zrp + c * zrq; }
        } }
    }
    let mut order: Vec<usize> = (0..mm).collect(); order.sort_by(|x, y| a[*x][*x].partial_cmp(&a[*y][*y]).unwrap());
    order.into_iter().take(k).map(|c| { let mut vec = vec![0.0; n]; for j in 0..mm { let y = z[j][c]; for i in 0..n { vec[i] += y * vs[j][i]; } } (a[c][c], vec) }).collect()
}

pub fn review(f: &Fields, b: &Board, spec: &Spec, issues: &mut Vec<Issue>) -> Value {
    let (w, h, n) = (f.w, f.h, f.w * f.h);
    let rail = |nm: &str| spec.voltages.get(nm).copied().or_else(|| rail_from_name(nm));
    let real: Vec<usize> = (0..LAYERS.len()).filter(|k| f.net_ids[*k].iter().filter(|x| **x != u16::MAX).count() > n * 15 / 100).collect();
    let gid = crate::fields::board_ground(spec, b, &f.net_names).and_then(|g| f.net_names.iter().position(|x| *x == g.0)).map(|i| i as u16);
    // the strong tones, each a frequency range (a tuned LO sweeps): switching harmonics (trapezoid envelope within 40 dB
    // of the fundamental) and the named clocks: a crystal's pins, or an LC oscillator no chip drives, swing nearly as a
    // sine, so only the first five harmonics count as strong; a logic clock's run on (or the spec's harmonics)
    let sine_net = |net: &str| b.pads.iter().any(|p| p.net == net && (p.reference.starts_with('Y') || p.reference.starts_with('X'))) || !b.pads.iter().any(|p| p.net == net && p.reference.starts_with('U'));
    let mut tones: Vec<(f64, f64, String)> = Vec::new();
    for (net, sw) in &spec.switching { let fs = sw.freqHz.unwrap_or(500e3); let tr = sw.riseNs.unwrap_or(8.0) * 1e-9; let knee = 1.0 / (std::f64::consts::PI * tr); let top = (knee * 10.0).min(spec.emcFmaxHz.unwrap_or(1e9)); let mut k = 1.0; while k * fs <= top { tones.push((k * fs, k * fs, format!("{net} harmonic {}", k as u64))); k += if k * fs < 100e6 { 1.0 } else { (k * 0.01).max(1.0).floor() }; } }
    for (net, c) in &spec.clocks { let sine = sine_net(net); let nh = c.harmonics().unwrap_or(if sine { 5 } else { ((spec.emcFmaxHz.unwrap_or(1e9) * 2.0) / c.top_hz().max(1.0)) as usize }).min(4000); for k in 1..=nh.max(1) { let kf = k as f64; tones.push((kf * c.lo_hz(), kf * c.top_hz(), format!("{net} {} harmonic {k}", if sine { "oscillator" } else { "clock" }))); } }
    let hz_txt = |a: f64, z: f64| -> String { let one = |x: f64| if x >= 1e9 { format!("{:.2} GHz", x / 1e9) } else if x >= 1e6 { format!("{:.1} MHz", x / 1e6) } else { format!("{:.0} kHz", x / 1e3) }; if (z - a).abs() < a * 1e-3 { one(a) } else { format!("{}–{}", one(a), one(z)) } };
    let mut out = json!({});
    // ---- cavities --------------------------------------------------------------------------------------------
    let mut cavities = Vec::new();
    // every pour facing ground on the next copper layer, in either direction (issue 9: not just the first two layers)
    let layer_pairs: Vec<(usize, usize)> = real.windows(2).flat_map(|w| [(w[0], w[1]), (w[1], w[0])]).collect();
    if let Some(g) = gid { for (k0, k1) in layer_pairs {
        let (hmm, er) = crate::fields::gap_between(spec, b, k0, k1);
        let mut pairs: std::collections::BTreeMap<u16, usize> = Default::default();
        for i in 0..n { let a = f.net_ids[k0][i]; if a != u16::MAX && a != g && f.net_ids[k1][i] == g { *pairs.entry(a).or_default() += 1; } }
        for (id, cnt) in pairs {
            let name = f.net_names[id as usize].clone();
            if !(rail(&name).is_some() || spec.switching.contains_key(&name)) || (cnt as f64 * f.d * f.d) < 20.0 { continue; }
            // coarsen to at most ~5000 cells
            let s = ((cnt as f64 / 5000.0).sqrt().ceil() as usize).max(1); let (cw, ch) = (w.div_ceil(s), h.div_ceil(s)); let hc = f.d * s as f64;
            let mut mask = vec![false; cw * ch];
            for i in 0..n { if f.net_ids[k0][i] == id && f.net_ids[k1][i] == g { mask[(i / w / s) * cw + (i % w) / s] = true; } }
            // the largest connected piece
            let mut comp = vec![usize::MAX; cw * ch]; let mut best = (0usize, 0usize);
            for st in 0..cw * ch { if !mask[st] || comp[st] != usize::MAX { continue; } let mut stack = vec![st]; comp[st] = st; let mut sz = 0; while let Some(v) = stack.pop() { sz += 1; let (r, c) = (v / cw, v % cw); for (dr, dc) in [(0i64, 1i64), (0, -1), (1, 0), (-1, 0)] { let (rr, cc) = (r as i64 + dr, c as i64 + dc); if rr < 0 || cc < 0 || rr >= ch as i64 || cc >= cw as i64 { continue; } let u = rr as usize * cw + cc as usize; if mask[u] && comp[u] == usize::MAX { comp[u] = st; stack.push(u); } } } if sz > best.1 { best = (st, sz); } }
            let cells: Vec<usize> = (0..cw * ch).filter(|i| comp[*i] == best.0).collect(); let nn = cells.len(); if nn < 12 { continue; }
            let mut idx = vec![usize::MAX; cw * ch]; for (j, c) in cells.iter().enumerate() { idx[*c] = j; }
            let nb: Vec<Vec<usize>> = cells.iter().map(|&c| { let (r, cc) = (c / cw, c % cw); let mut v = Vec::new(); if cc + 1 < cw && idx[c + 1] != usize::MAX { v.push(idx[c + 1]); } if cc > 0 && idx[c - 1] != usize::MAX { v.push(idx[c - 1]); } if r + 1 < ch && idx[c + cw] != usize::MAX { v.push(idx[c + cw]); } if r > 0 && idx[c - cw] != usize::MAX { v.push(idx[c - cw]); } v }).collect();
            let apply = |x: &[f64], y: &mut [f64]| { for i in 0..nn { let mut s_ = nb[i].len() as f64 * x[i]; for &j in &nb[i] { s_ -= x[j]; } y[i] = s_; } };
            let modes = lowest_modes(nn, 5, &apply);
            let mut ms = Vec::new();
            for (lam, vec) in modes.iter().skip(1) {   // the first is the constant: the plate's own static charge
                if *lam <= 1e-9 { continue; }
                let k = lam.sqrt() / hc; let fr = C_MM * k / (2.0 * std::f64::consts::PI * er.sqrt());
                let (mut amax, mut at) = (0.0, 0usize); for (j, v) in vec.iter().enumerate() { if v.abs() > amax { amax = v.abs(); at = j; } }
                let atxy = [f.x0 + ((cells[at] % cw) as f64 + 0.5) * hc, f.y0 + ((cells[at] / cw) as f64 + 0.5) * hc];
                let near = tones.iter().filter(|t| fr >= t.0 * 0.95 && fr <= t.1 * 1.05).min_by(|a, z| ((a.0 + a.1) / 2.0 - fr).abs().partial_cmp(&((z.0 + z.1) / 2.0 - fr).abs()).unwrap());
                ms.push(json!({"fHz": fr.round(), "antinode": [r2(atxy[0]), r2(atxy[1])], "near": near.map(|t| t.2.clone())}));
                if let Some(t) = near { issues.push(Issue { severity: "medium".into(), kind: "cavity-resonance".into(), what: format!("{} on {} over {} on {}: a plane mode at {:.0} MHz sits within 5 % of the {} ({})", name, LAYERS[k0], f.net_names[g as usize], LAYERS[k1], fr / 1e6, t.2, hz_txt(t.0, t.1)), at: Some([r2(atxy[0]), r2(atxy[1])]), layer: Some(LAYERS[k0].into()), net: Some(name.clone()), reference: None, fix: format!("Damp it where it rings loudest, ({:.1}, {:.1}): a 1 to 10 nF capacitor from {} to ground there, or ground vias through the pour.", atxy[0], atxy[1], name), action: Some(json!({"type": "add-part", "what": "capacitor", "value": "4.7nF", "between": [name, f.net_names[g as usize]], "at": [r2(atxy[0]), r2(atxy[1])], "why": "damp a plane resonance"})) }); }
            }
            // the shapes over the cavity's own bounding box, three decimals: the app animates them
            let (r0, r1) = (cells.iter().map(|c| c / cw).min().unwrap(), cells.iter().map(|c| c / cw).max().unwrap()); let (c0, c1) = (cells.iter().map(|c| c % cw).min().unwrap(), cells.iter().map(|c| c % cw).max().unwrap()); let (bw, bh) = (c1 - c0 + 1, r1 - r0 + 1);
            let shapes: Vec<Vec<f64>> = modes.iter().skip(1).take(3).map(|(_, v)| { let mx = v.iter().fold(0.0f64, |a, b| a.max(b.abs())).max(1e-12); let mut sub = vec![0.0f64; bw * bh]; for (j, c) in cells.iter().enumerate() { sub[(c / cw - r0) * bw + (c % cw - c0)] = (v[j] / mx * 1000.0).round() / 1000.0; } sub }).collect();
            let mask_sub: Vec<u8> = (0..bw * bh).map(|i| { let c = (r0 + i / bw) * cw + c0 + i % bw; if comp[c] == best.0 { 1 } else { 0 } }).collect();
            cavities.push(json!({"net": name, "layer": LAYERS[k0], "overLayer": LAYERS[k1], "er": r2(er), "over": f.net_names[g as usize], "areaMm2": r2(cnt as f64 * f.d * f.d), "gapMm": hmm, "grid": {"cols": bw, "rows": bh, "cellMm": hc, "x0": f.x0 + c0 as f64 * hc, "y0": f.y0 + r0 as f64 * hc, "mask": mask_sub}, "modes": ms, "shapes": shapes}));
        }
    }
    }
    out["cavities"] = json!(cavities);
    // ---- crosstalk ----------------------------------------------------------------------------------------------
    // aggressors: a switching node is a step of its swing; a clock or oscillator is a sine at its amplitude (the spec's,
    // else 1 V for a crystal or LC oscillator, the logic rail for a chip-driven clock) and its mid frequency
    struct Agg { net: String, amp: f64, hz: Option<f64> }
    let logic_v = rail("+3V3").or_else(|| rail("VCC")).unwrap_or(3.3);
    let mut aggressors: Vec<Agg> = spec.switching.iter().map(|(k, s)| Agg { net: k.clone(), amp: s.highV - s.lowV, hz: None }).collect();
    for nm in f.net_names.iter().filter(|nm| is_switch_name(nm) && !spec.switching.contains_key(*nm)) { aggressors.push(Agg { net: nm.clone(), amp: f.net_names.iter().filter_map(|x| rail(x)).fold(0.0, f64::max), hz: None }); }
    for (k, c) in &spec.clocks { aggressors.push(Agg { net: k.clone(), amp: c.amplitude().unwrap_or(if sine_net(k) { 1.0 } else { logic_v }), hz: Some(c.mid_hz()) }); }
    let mut xt = Vec::new();
    let budget = spec.crosstalkMv.unwrap_or(50.0);
    let thick = |k: usize| spec.copperUmByLayer.get(LAYERS[k]).copied().unwrap_or(spec.copperUm) / 1000.0;
    // the nearest real copper layer to k: the plane a victim on k sits over
    let plane_gap = |k: usize| -> (f64, f64) { let i = real.iter().position(|x| *x == k).unwrap_or(0); let nb: Vec<usize> = [i.checked_sub(1), Some(i + 1)].into_iter().flatten().filter_map(|j| real.get(j).copied()).collect(); nb.iter().map(|&o| crate::fields::gap_between(spec, b, k, o)).min_by(|a, z| a.0.partial_cmp(&z.0).unwrap()).unwrap_or((1.465, 4.3)) };
    // a net's shunt to AC ground: the capacitors and resistors from it to ground or any rail (rails are AC ground)
    let shunt = |net: &str| -> (f64, f64) { let mut cs = 0.0; let mut g = 0.0; let mut refs: Vec<&str> = b.pads.iter().filter(|p| p.net == net).map(|p| p.reference.as_str()).collect(); refs.sort(); refs.dedup();
        for r in refs { let ns: Vec<&str> = b.pads.iter().filter(|p| p.reference == r).map(|p| p.net.as_str()).collect(); if ns.len() != 2 { continue; } let other = if ns[0] == net { ns[1] } else { ns[0] }; if !(is_ground(other) || rail(other).is_some()) { continue; }
            let val = b.values.get(r).map(String::as_str).unwrap_or("");
            if r.starts_with('C') { if let Some(c) = crate::emc::parse_value(val, 'F') { cs += c; } }
            else if r.starts_with('R') { let u = val.to_ascii_uppercase().replace("OHM", ""); let t = u.trim(); let ohm = if let Some(i) = t.find(|c: char| c == 'K' || c == 'M' || c == 'R') { let m = match &t[i..i + 1] { "K" => 1e3, "M" => 1e6, _ => 1.0 }; format!("{}.{}", &t[..i], &t[i + 1..]).trim_end_matches('.').parse::<f64>().unwrap_or(0.0) * m } else { t.parse::<f64>().unwrap_or(0.0) }; if ohm > 0.0 { g += 1.0 / ohm; } } }
        (cs, g) };
    for ag in &aggressors {
        let agg = &ag.net; let dv = ag.amp;
        let Some(aid) = f.net_names.iter().position(|x| x == agg) else { continue }; let aid = aid as u16;
        let is_clock = ag.hz.is_some();
        // one discrete part away is the aggressor's own circuit (a bootstrap capacitor, an oscillator's crystal or
        // transistor): coupled through that part far more than through stray copper; its circuit paths are reported below
        let disc = |r: &str| r.starts_with('C') || (is_clock && !r.starts_with('U') && !r.starts_with('J') && !r.starts_with("TP"));
        let one: std::collections::BTreeSet<String> = { let mut out = std::collections::BTreeSet::new(); for p in b.pads.iter().filter(|p| p.net == *agg && disc(&p.reference)) { for q in b.pads.iter().filter(|q| q.reference == p.reference && q.net != *agg) { out.insert(q.net.clone()); } } out };
        for (vid, vname) in f.net_names.iter().enumerate() {
            let vid = vid as u16; if vid == aid || Some(vid) == gid || f.flow.contains_key(vname) || rail(vname).is_some() || aggressors.iter().any(|a| a.net == *vname) || one.contains(vname.as_str()) { continue; }
            // side coupling within 3 mm on each layer; a plane h below shields it, roughly as 1 / (1 + (s / 2h)²)
            let (mut cm, mut closest, mut at, mut run) = (0.0f64, f64::MAX, [0.0, 0.0], 0usize); let mut cg_cu = 0.0;
            let rad = (3.0 / f.d).ceil() as i64;
            for &k in &real { let (hk, erk) = plane_gap(k); for i in 0..n {
                if f.net_ids[k][i] != vid { continue; } cg_cu += E0 * erk * f.d * f.d / hk.max(0.05);
                let (r, c) = ((i / w) as i64, (i % w) as i64); let mut dmin = f64::MAX;
                for dr in -rad..=rad { for dc in -rad..=rad { let (rr, cc) = (r + dr, c + dc); if rr < 0 || cc < 0 || rr >= h as i64 || cc >= w as i64 { continue; } if f.net_ids[k][rr as usize * w + cc as usize] == aid { dmin = dmin.min(((dr * dr + dc * dc) as f64).sqrt() * f.d); } } }
                if dmin <= 3.0 { let s_ = dmin.max(0.05); let eff = (erk + 1.0) / 2.0; let shield = 1.0 / (1.0 + (s_ / (2.0 * hk)).powi(2)); cm += E0 * eff * (thick(k) / s_ + 0.6 * (1.0 + f.d / s_).ln()) * f.d * shield; if dmin <= 1.0 { run += 1; } if dmin < closest { closest = dmin; at = [f.x0 + (c as f64 + 0.5) * f.d, f.y0 + (r as f64 + 0.5) * f.d]; } }
            } }
            // overlap with the aggressor on each adjacent copper layer, at that pair's gap and εr
            let mut over = 0.0; for wpair in real.windows(2) { let (ga, gb) = (wpair[0], wpair[1]); let (hg, eg) = crate::fields::gap_between(spec, b, ga, gb); let mut o = 0.0; for i in 0..n { if (f.net_ids[ga][i] == vid && f.net_ids[gb][i] == aid) || (f.net_ids[gb][i] == vid && f.net_ids[ga][i] == aid) { o += f.d * f.d; } } over += o; cm += E0 * eg * o / hg; }
            if cm <= 0.0 { continue; }
            // the victim to ground: its copper over its nearest plane, about 3 pF of pins, and the parts on it (issue 9)
            let (cpart, gpart) = shunt(vname);
            let cg = cg_cu + 3e-12 + cpart;
            let v = match ag.hz { None => dv * cm / (cm + cg), Some(hz) => { let wv = 2.0 * std::f64::consts::PI * hz; dv * wv * cm / ((wv * (cm + cg)).powi(2) + gpart * gpart).sqrt() } };
            let mv = v * 1e3;
            let vb = spec.victimBudgetsMv.get(vname).copied().unwrap_or(budget);
            xt.push(json!({"aggressor": agg, "victim": vname, "kind": if is_clock { "sine" } else { "step" }, "atHz": ag.hz, "amplitudeV": dv, "cmPf": (cm * 1e15).round() / 1e3, "victimPf": r2(cg * 1e12), "victimShuntOhm": if gpart > 0.0 { json!(r2(1.0 / gpart)) } else { Value::Null }, "closestMm": if closest < f64::MAX { r2(closest) } else { -1.0 }, "runMm": r2(run as f64 * f.d), "overlapMm2": r2(over), "coupledMv": (mv * 1e4).round() / 1e4, "coupledUv": (mv * 1e3 * 10.0).round() / 10.0, "budgetMv": vb, "at": [r2(at[0]), r2(at[1])]}));
            if mv > vb {
                let amt = if mv >= 1.0 { format!("{:.0} mV", mv) } else { format!("{:.0} µV", mv * 1e3) };
                let what = if is_clock { format!("{} picks up about {} of {} ({:.0} MHz) through {:.0} fF (closest {:.2} mm; {:.0} pF to ground on the victim)", vname, amt, agg, ag.hz.unwrap_or(0.0) / 1e6, cm * 1e15, closest, cg * 1e12) } else { format!("{} picks up about {} from each {} edge ({:.2} pF, closest {:.2} mm, {:.1} mm run)", vname, amt, agg, cm * 1e12, closest, run as f64 * f.d) };
                issues.push(Issue { severity: if mv > 4.0 * vb { "medium".into() } else { "low".into() }, kind: "crosstalk".into(), what, at: Some([r2(at[0]), r2(at[1])]), layer: None, net: Some(vname.clone()), reference: None, fix: format!("Route {} at least {:.1} mm from {} (three times the spacing), or run a ground guard between them, and keep {} off the copper under {}.", vname, (closest * 3.0).max(1.0), agg, vname, agg), action: Some(json!({"type": "reroute", "net": vname, "awayFrom": agg, "minSpacingMm": r2((closest * 3.0).max(1.0)), "at": [r2(at[0]), r2(at[1])], "or": "ground guard between"})) });
            }
        }
    }
    xt.sort_by(|a, z| z["coupledMv"].as_f64().partial_cmp(&a["coupledMv"].as_f64()).unwrap());
    out["crosstalk"] = json!(xt);
    out["budgetMv"] = json!(budget);
    out["victimBudgetsMv"] = json!(spec.victimBudgetsMv);
    // ---- the frequency plan: a tone or harmonic inside a receive band or its image (issue 9) --------------------------
    let mut plan = Vec::new();
    for (bn, bd) in &spec.bands {
        let image = bd.ifHz.map(|fi| if bd.lowSideLo { (bd.minHz - 2.0 * fi, bd.maxHz - 2.0 * fi) } else { (bd.minHz + 2.0 * fi, bd.maxHz + 2.0 * fi) });
        for (net, c) in &spec.clocks {
            // every harmonic up to the band, weak ones included: a receiver hears microvolts
            let nh = ((bd.maxHz.max(image.map_or(0.0, |x| x.1)) / c.lo_hz().max(1.0)).floor() as usize).min(c.harmonics().unwrap_or(400).max(1) * 1000).min(400);
            for k in 1..=nh.max(1) {
                if k == 1 && c.swept() { continue; }   // a tuned LO sits beside its own band by design: its leakage is a circuit path
                let (a, z) = (k as f64 * c.lo_hz(), k as f64 * c.top_hz());
                for (what, lo, hi) in [("band", bd.minHz, bd.maxHz)].into_iter().chain(image.map(|x| ("image band", x.0, x.1))) {
                    if z >= lo && a <= hi {
                        let sev = if k <= 5 { "medium" } else { "low" };
                        plan.push(json!({"tone": net, "harmonic": k, "fromHz": a, "toHz": z, "band": bn, "where": what, "bandHz": [lo, hi]}));
                        issues.push(Issue { severity: sev.into(), kind: "tone-in-band".into(), what: format!("harmonic {} of {} ({}) lands in the {} {} ({})", k, net, hz_txt(a, z), bn, what, hz_txt(lo, hi)), at: None, layer: None, net: Some(net.clone()), reference: None, fix: format!("Keep {} and its traces away from the {} front end and shield or filter it, or move its frequency so no harmonic falls in {} ({}).", net, bn, hz_txt(lo, hi), what), action: Some(json!({"type": "review", "net": net, "why": "a tone in a receive band", "harmonic": k, "band": bn, "where": what})) });
                    }
                }
            }
        }
    }
    // aggressors with no oscillator net of their own (issue 9): a full-speed USB pair's 12 MHz edges and harmonics, and an FM
    // IF limiter's odd harmonics of the IF (a hard limiter squares the IF, so 3x, 5x, 7x... are strong)
    let mut extra: Vec<(String, f64, bool)> = Vec::new();
    let usb = f.net_names.iter().any(|n| { let u = n.to_ascii_uppercase(); u.contains("USB_D") || u.ends_with("D+") || u.ends_with("D-") || u.ends_with("USB_DP") || u.ends_with("USB_DM") });
    if usb { extra.push(("USB full-speed edges (12 MHz)".into(), 12e6, false)); if spec.rf.usbHighSpeed == Some(true) { extra.push(("USB high-speed (480 MHz)".into(), 480e6, false)); } }
    for (bn, bd) in &spec.bands { if let Some(fi) = bd.ifHz { if spec.rf.limiter.unwrap_or_else(|| bn.to_ascii_uppercase().contains("FM")) { extra.push((format!("{bn} IF limiter ({:.1} MHz, odd harmonics)", fi / 1e6), fi, true)); } } }
    for (name, base, odd) in &extra {
        for (bn, bd) in &spec.bands {
            let image = bd.ifHz.map(|fi| if bd.lowSideLo { (bd.minHz - 2.0 * fi, bd.maxHz - 2.0 * fi) } else { (bd.minHz + 2.0 * fi, bd.maxHz + 2.0 * fi) });
            let top = bd.maxHz.max(image.map_or(0.0, |x| x.1));
            let mut k = if *odd { 3 } else { 1 };
            while (k as f64) * base <= top && k <= 400 {
                let fk = k as f64 * base;
                for (what, lo, hi) in [("band", bd.minHz, bd.maxHz)].into_iter().chain(image.map(|x| ("image band", x.0, x.1))) {
                    if fk >= lo && fk <= hi {
                        plan.push(json!({"tone": name, "harmonic": k, "fromHz": fk, "toHz": fk, "band": bn, "where": what, "bandHz": [lo, hi]}));
                        issues.push(Issue { severity: if k <= 5 { "medium".into() } else { "low".into() }, kind: "tone-in-band".into(), what: format!("harmonic {} of {} ({:.2} MHz) lands in the {} {} ({})", k, name, fk / 1e6, bn, what, hz_txt(lo, hi)), at: None, layer: None, net: None, reference: None, fix: if name.starts_with("USB") { "Keep the USB pair and its connector away from the front end and over unbroken ground; a common-mode choke on the pair cuts the edge harmonics.".into() } else { "Filter and shield the limiter's output, keep its traces away from the front end, and consider an IF that keeps its odd harmonics out of the band.".into() }, action: None });
                    }
                }
                k += if *odd { 2 } else { 1 };
            }
        }
    }
    out["frequencyPlan"] = json!(plan);
    // ---- circuit paths: from each oscillator through the parts to each sensitive net (issue 9: report, don't hide) -----
    let mut paths = Vec::new();
    for (net, _) in &spec.clocks { for target in &spec.sensitiveNets {
        // breadth-first through two-terminal and transistor parts, never through a rail or ground, at most 8 parts
        let mut prev: std::collections::BTreeMap<String, (String, String)> = Default::default(); let mut q = std::collections::VecDeque::new(); q.push_back((net.clone(), 0usize)); prev.insert(net.clone(), (String::new(), String::new()));
        let mut found = false;
        while let Some((cur, d)) = q.pop_front() { if cur == *target { found = true; break; } if d >= 8 { continue; }
            let mut refs: Vec<&str> = b.pads.iter().filter(|p| p.net == cur && !p.reference.starts_with('U') && !p.reference.starts_with('J') && !p.reference.starts_with("TP")).map(|p| p.reference.as_str()).collect(); refs.sort(); refs.dedup();
            for r in refs { for p in b.pads.iter().filter(|p| p.reference == r && p.net != cur) { if is_ground(&p.net) || rail(&p.net).is_some() || p.net.is_empty() || prev.contains_key(&p.net) { continue; } prev.insert(p.net.clone(), (cur.clone(), r.to_string())); q.push_back((p.net.clone(), d + 1)); } } }
        if found { let mut hops = Vec::new(); let mut cur = target.clone(); while let Some((p, r)) = prev.get(&cur).cloned() { if p.is_empty() { break; } hops.push(json!({"part": r, "value": b.values.get(&r), "from": p, "to": cur})); cur = p; } hops.reverse();
            let txt = hops.iter().map(|h| format!("{} {}", h["part"].as_str().unwrap_or(""), h["value"].as_str().unwrap_or(""))).collect::<Vec<_>>().join(" → ");
            issues.push(Issue { severity: "low".into(), kind: "circuit-path".into(), what: format!("{} reaches {} through {} part(s): {}", net, target, hops.len(), txt), at: None, layer: None, net: Some(target.clone()), reference: None, fix: format!("Simulate this path (ngspice AC with the stages' small-signal models) for the level of {} at {}; add isolation (a buffer, a lower-value coupling capacitor, a shield) if it is over the target.", net, target), action: Some(json!({"type": "review", "why": "circuit path from an oscillator to a sensitive net", "from": net, "to": target, "hops": hops})) });
            paths.push(json!({"from": net, "to": target, "hops": hops})); }
    } }
    out["circuitPaths"] = json!(paths);
    out["tones"] = json!(tones.iter().take(40).map(|t| json!({"fHz": t.0, "toHz": t.1, "what": t.2})).collect::<Vec<_>>());
    let _ = b; out
}

#[cfg(test)] mod tests {
    use super::*;
    /// a 20 x 10 mm rectangle: the first mode is (1,0), f = c / (2 a √εr)
    #[test] fn rectangle_modes() {
        let (cw, ch, hc) = (40usize, 20usize, 0.5); let n = cw * ch;
        let apply = |x: &[f64], y: &mut [f64]| { for i in 0..n { let (r, c) = (i / cw, i % cw); let mut s = 0.0; let mut k = 0.0; if c + 1 < cw { s += x[i + 1]; k += 1.0; } if c > 0 { s += x[i - 1]; k += 1.0; } if r + 1 < ch { s += x[i + cw]; k += 1.0; } if r > 0 { s += x[i - cw]; k += 1.0; } y[i] = k * x[i] - s; } };
        let m = lowest_modes(n, 3, &apply);
        assert!(m[0].0.abs() < 1e-6, "constant mode {}", m[0].0);
        let f1 = C_MM * (m[1].0.sqrt() / hc) / (2.0 * std::f64::consts::PI * 4.3f64.sqrt());
        let exact = C_MM / (2.0 * 20.0 * 4.3f64.sqrt());
        assert!((f1 - exact).abs() / exact < 0.01, "{f1} vs {exact}");
    }
}