Adom Fields

//! Trace impedance by a 2-D field solve of the cross-section (0.9.0): the trace (or pair) between its planes, Laplace's
//! equation on a grid with the stackup's dielectric, solved once with the dielectric and once in air:
//!   C  = Σ ε·ΔV² over the grid's edges (the field energy at 1 V),  C0 = the same in air,
//!   Z0 = 1 / (c·√(C·C0)),  εeff = C / C0,  delay = √εeff / c,  L = 1 / (c²·C0)
//! For a pair the even and odd modes give Zeven, Zodd and Zdiff = 2·Zodd. Any layer: a microstrip (one plane under, air or
//! solder mask over) or a stripline (planes on both sides). Validated against the closed forms in the tests.
use crate::board::Board;
use crate::fields::{Issue, Spec};
use crate::flow::{is_ground, rail_from_name};
use serde_json::{json, Value};

const E0: f64 = 8.854187817e-12; const C0: f64 = 299_792_458.0;

/// a cross-section: widths and positions in metres; `below`/`above` = distance to a plane under/over the trace bottom
#[derive(Clone, Debug)]
pub struct Xs { pub w: f64, pub t: f64, pub gap: Option<f64>, pub below: f64, pub above: Option<f64>, pub er: f64, pub mask: Option<(f64, f64)> }

#[derive(Clone, Debug)] pub struct Zr { pub z0: f64, pub eeff: f64, pub zodd: Option<f64>, pub zeven: Option<f64>, pub zdiff: Option<f64>, pub cpf_m: f64, pub lnh_m: f64 }

/// a tensor-product grid, fine at the copper edges and growing outward (finite volume: each edge's weight is
/// ε · (the dual cell's width across it) / (the node spacing along it))
struct Grid { xs: Vec<f64>, ys: Vec<f64>, eps: Vec<f64>, fixed: Vec<Option<f64>> }
impl Grid {
    fn nx(&self) -> usize { self.xs.len() }
    fn ny(&self) -> usize { self.ys.len() }
    /// the edges: (a, b, weight without ε0)
    fn edges(&self) -> Vec<(usize, usize, f64)> {
        let (nx, ny) = (self.nx(), self.ny()); let mut e = Vec::with_capacity(2 * nx * ny);
        for j in 0..ny { for i in 0..nx { let k = j * nx + i;
            // ε lives on the cells between the nodes: an edge's weight sums the two cells it borders, each by its half
            let ec = |ci: usize, cj: usize| self.eps[cj * (nx - 1) + ci];
            if i + 1 < nx { let m = k + 1; let mut w = 0.0; if j > 0 { w += ec(i, j - 1) * (self.ys[j] - self.ys[j - 1]) / 2.0; } if j + 1 < ny { w += ec(i, j) * (self.ys[j + 1] - self.ys[j]) / 2.0; } e.push((k, m, w / (self.xs[i + 1] - self.xs[i]))); }
            if j + 1 < ny { let m = k + nx; let mut w = 0.0; if i > 0 { w += ec(i - 1, j) * (self.xs[i] - self.xs[i - 1]) / 2.0; } if i + 1 < nx { w += ec(i, j) * (self.xs[i + 1] - self.xs[i]) / 2.0; } e.push((k, m, w / (self.ys[j + 1] - self.ys[j]))); }
        } }
        e
    }
    /// conjugate gradient (Jacobi preconditioned) on the free nodes
    fn solve(&self) -> (Vec<f64>, Vec<(usize, usize, f64)>) {
        let n = self.nx() * self.ny(); let edges = self.edges();
        let mut adj: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
        for &(a, b, w) in &edges { adj[a].push((b, w)); adj[b].push((a, w)); }
        let free = |k: usize| self.fixed[k].is_none();
        let mut rhs = vec![0.0; n]; let mut diag = vec![1.0; n];
        for k in 0..n { if !free(k) { continue; } diag[k] = adj[k].iter().map(|x| x.1).sum::<f64>().max(1e-300); for &(m, w) in &adj[k] { if let Some(v) = self.fixed[m] { rhs[k] += w * v; } } }
        let apply = |x: &[f64], y: &mut [f64]| { for k in 0..n { if !free(k) { y[k] = 0.0; continue; } let mut s = diag[k] * x[k]; for &(m, w) in &adj[k] { if free(m) { s -= w * x[m]; } } y[k] = s; } };
        let mut x = vec![0.0; n]; let mut r = rhs.clone(); let b2: f64 = rhs.iter().map(|v| v * v).sum::<f64>().max(1e-300);
        let mut z: Vec<f64> = (0..n).map(|k| r[k] / diag[k]).collect(); let mut p = z.clone(); let mut rz: f64 = (0..n).map(|k| r[k] * z[k]).sum();
        let mut ap = vec![0.0; n];
        for _ in 0..30000 {
            apply(&p, &mut ap); let pap: f64 = (0..n).map(|k| p[k] * ap[k]).sum(); if pap.abs() < 1e-300 { break; }
            let a = rz / pap; for k in 0..n { x[k] += a * p[k]; r[k] -= a * ap[k]; }
            if r.iter().map(|v| v * v).sum::<f64>() / b2 < 1e-22 { break; }
            for k in 0..n { z[k] = r[k] / diag[k]; }
            let rz2: f64 = (0..n).map(|k| r[k] * z[k]).sum(); let bt = rz2 / rz; rz = rz2; for k in 0..n { p[k] = z[k] + bt * p[k]; }
        }
        ((0..n).map(|k| self.fixed[k].unwrap_or(x[k])).collect(), edges)
    }
    /// C·V² per unit length = Σ w·ε0·ΔV²
    fn energy2(v: &[f64], edges: &[(usize, usize, f64)]) -> f64 { edges.iter().map(|&(a, b, w)| w * (v[a] - v[b]).powi(2)).sum::<f64>() * E0 }
}

/// nodes from lo to hi: every key point, and between them spacing that starts at `fine` at each key point and grows by
/// 25 % per cell up to `coarse`
fn axis(lo: f64, hi: f64, mut keys: Vec<f64>, fine: f64, coarse: f64) -> Vec<f64> {
    keys.push(lo); keys.push(hi); keys.retain(|k| *k >= lo - 1e-15 && *k <= hi + 1e-15);
    keys.sort_by(|a, b| a.partial_cmp(b).unwrap()); keys.dedup_by(|a, b| (*a - *b).abs() < fine * 0.25);
    let mut out = vec![keys[0]];
    for w in keys.windows(2) {
        let (a, b) = (w[0], w[1]); let len = b - a;
        // steps growing from both ends
        let mut left = Vec::new(); let mut step = fine; let mut pos = 0.0;
        while pos + step < len / 2.0 { pos += step; left.push(pos); step = (step * 1.25).min(coarse); }
        let mut pts: Vec<f64> = left.iter().map(|p| a + p).collect();
        let right: Vec<f64> = left.iter().rev().map(|p| b - p).filter(|q| *q > a + pos + fine * 0.5).collect();
        pts.extend(right); if pts.is_empty() && len > 2.0 * fine { pts.push(a + len / 2.0); }
        out.extend(pts); out.push(b);
    }
    out.sort_by(|a, b| a.partial_cmp(b).unwrap()); out.dedup_by(|a, b| (*a - *b).abs() < 1e-15); out
}

/// the grid for a cross-section; `air` replaces every dielectric by 1; `v1`, `v2` are the conductors' potentials
fn build(x: &Xs, air: bool, v1: f64, v2: f64) -> Grid {
    let s = x.gap.unwrap_or(0.0);
    let spans: Vec<(f64, f64, f64)> = if x.gap.is_some() { vec![(-(s / 2.0 + x.w), -s / 2.0, v1), (s / 2.0, s / 2.0 + x.w, v2)] } else { vec![(-x.w / 2.0, x.w / 2.0, v1)] };
    let reach = x.below.max(x.above.unwrap_or(0.0));
    let side = (8.0 * reach).max(4.0 * x.w); let half = spans.last().unwrap().1 + side;
    let (yb, yt) = (x.below, x.below + x.t);
    let top = match x.above { Some(a) => yt + a, None => yt + 10.0 * x.below.max(x.w) };
    let fine_x = (x.w / 24.0).min(if x.gap.is_some() { s / 8.0 } else { f64::MAX }).min(x.t);
    let fine_y = (x.t / 3.0).min(x.below / 16.0);
    let mut kx = vec![0.0]; for sp in &spans { kx.push(sp.0); kx.push(sp.1); }
    let mut ky = vec![yb, yt]; if let Some((mt, _)) = x.mask { ky.push(yt + mt); }
    let xs = axis(-half, half, kx, fine_x, reach.max(x.w) / 3.0);
    let ys = axis(0.0, top, ky, fine_y, reach.max(x.w) / 3.0);
    let (nx, ny) = (xs.len(), ys.len());
    let mut g = Grid { xs: xs.clone(), ys: ys.clone(), eps: vec![1.0; (nx - 1) * (ny - 1)], fixed: vec![None; nx * ny] };
    for cj in 0..ny - 1 { for ci in 0..nx - 1 { let yc = (ys[cj] + ys[cj + 1]) / 2.0;
        let in_diel = yc < yb || x.above.is_some();
        let in_mask = x.above.is_none() && x.mask.is_some_and(|(mt, _)| yc > yb && yc < yt + mt);
        g.eps[cj * (nx - 1) + ci] = if air { 1.0 } else if in_diel { x.er } else if in_mask { x.mask.unwrap().1 } else { 1.0 };
    } }
    for j in 0..ny { for i in 0..nx { let k = j * nx + i; let (xx, y) = (xs[i], ys[j]);
        if j == 0 || j == ny - 1 || i == 0 || i == nx - 1 { g.fixed[k] = Some(0.0); continue; }
        if y >= yb - 1e-15 && y <= yt + 1e-15 { for sp in &spans { if xx >= sp.0 - 1e-15 && xx <= sp.1 + 1e-15 { g.fixed[k] = Some(sp.2); } } }
    } }
    g
}

pub fn solve(x: &Xs) -> Zr {
    let c_of = |air: bool, v1: f64, v2: f64| { let g = build(x, air, v1, v2); let (v, e) = g.solve(); Grid::energy2(&v, &e) };
    if x.gap.is_none() {
        let (c, c0) = (c_of(false, 1.0, 0.0), c_of(true, 1.0, 0.0));
        let z0 = 1.0 / (C0 * (c * c0).sqrt());
        return Zr { z0, eeff: c / c0, zodd: None, zeven: None, zdiff: None, cpf_m: c * 1e12, lnh_m: 1.0 / (C0 * C0 * c0) * 1e9 };
    }
    // a pair: energy at (1,0), (1,1) and (1,-1); per line Codd = W(1,-1)/1 /2 ... with W = ½VᵀCV doubled: e = VᵀCV
    let e10 = c_of(false, 1.0, 0.0); let e11 = c_of(false, 1.0, 1.0); let e1m = c_of(false, 1.0, -1.0);
    let a10 = c_of(true, 1.0, 0.0); let a11 = c_of(true, 1.0, 1.0); let a1m = c_of(true, 1.0, -1.0);
    let (ce, co) = (e11 / 2.0, e1m / 2.0); let (ae, ao) = (a11 / 2.0, a1m / 2.0);
    let ze = 1.0 / (C0 * (ce * ae).sqrt()); let zo = 1.0 / (C0 * (co * ao).sqrt());
    let z0 = 1.0 / (C0 * (e10 * a10).sqrt());
    Zr { z0, eeff: co / ao, zodd: Some(zo), zeven: Some(ze), zdiff: Some(2.0 * zo), cpf_m: co * 1e12, lnh_m: 1.0 / (C0 * C0 * ao) * 1e9 }
}

fn r1(v: f64) -> f64 { (v * 10.0).round() / 10.0 }

/// protocol differential targets by net name (Ω): USB 90, PCIe 85, Ethernet/LVDS/HDMI/MIPI 100
fn diff_target(a: &str) -> Option<(f64, &'static str)> {
    let u = a.to_ascii_uppercase();
    if u.contains("USB") { Some((90.0, "USB")) } else if u.contains("PCIE") || u.contains("PCI_E") { Some((85.0, "PCIe")) } else if u.contains("ETH") || u.contains("MDI") || u.contains("LVDS") || u.contains("HDMI") || u.contains("TMDS") || u.contains("MIPI") || u.contains("CSI") || u.contains("DSI") { Some((100.0, "100 Ω differential")) } else { None }
}
/// the partner of a pair member by name: _P/_N, +/-, DP/DM, _H/_L
fn partner(n: &str) -> Option<String> {
    for (a, b) in [("_P", "_N"), ("+", "-"), ("DP", "DM"), ("_H", "_L"), ("P", "N")] { if let Some(s) = n.strip_suffix(a) { return Some(format!("{s}{b}")); } }
    None
}

pub fn review(b: &Board, spec: &Spec, stack: &Value, issues: &mut Vec<Issue>) -> Value {
    let real: Vec<String> = if spec.real_copper.is_empty() { b.copper.clone() } else { spec.real_copper.clone() };
    let gaps = if spec.real_gaps.is_empty() { crate::fields::slot_gaps(spec) } else { spec.real_gaps.clone() };
    if real.len() < 2 { return json!({"ran": false}); }
    let z: Vec<f64> = { let mut acc = 0.0; let mut v = vec![0.0]; for g in gaps.iter().take(real.len() - 1) { acc += g; v.push(acc); } v };
    let _ = stack; // the whole-layer plane classification is no longer the reference (issue 23): each segment's is found locally
    let er = spec.epsilonR.unwrap_or(4.3);
    let mask = if spec.rf.noSolderMask == Some(true) { None } else { Some((20e-6, 3.6)) };
    // the cross-section with its references given: `below` the nearest reference layer on one side, `above` the one on
    // the other side (a stripline), each a real-layer index
    let xs_for = |li: usize, w: f64, gap: Option<f64>, below_i: usize, above_i: Option<usize>| -> Xs {
        let t = spec.copperUmByLayer.get(&real[li]).copied().unwrap_or(spec.copperUm) * 1e-6;
        let d1 = (z[below_i] - z[li]).abs() * 1e-3;
        let outer = li == 0 || li == real.len() - 1;
        Xs { w: w * 1e-3, t, gap: gap.map(|g| g * 1e-3), below: d1.max(t), above: above_i.map(|j| ((z[j] - z[li]).abs() * 1e-3 - t).max(t)), er, mask: if outer && above_i.is_none() { mask } else { None } }
    };
    // ---- the reference copper under each trace, found locally (issue 23) -----------------------------------------------
    // ground and supply copper on every real layer, rasterised at 0.1 mm; a trace's reference on a layer is the share of
    // its length that layer's reference copper covers under the trace's centre and both edges
    let refnet = |n: &str| is_ground(n) || rail_from_name(n).is_some() || spec.voltages.contains_key(n);
    let mut ras = crate::raster::Raster::new(b, 0.1); ras.draw_board(b);
    let (rw, rh, rres, rx0, ry0) = (ras.w, ras.h, ras.res, ras.x0, ras.y0);
    let refmask: Vec<Vec<u16>> = real.iter().map(|l| { let mut m = vec![0u16; rw * rh]; let mut names: Vec<&String> = ras.masks.keys().filter(|(n, ll)| ll == l && refnet(n)).map(|(n, _)| n).collect(); names.sort();
        for (k, n) in names.iter().enumerate() { if let Some(mm) = ras.masks.get(&((*n).clone(), l.clone())) { for i in 0..rw * rh { if mm[i] != 0 && m[i] == 0 { m[i] = k as u16 + 1; } } } } m }).collect();
    let refnames: Vec<Vec<String>> = real.iter().map(|l| { let mut names: Vec<String> = ras.masks.keys().filter(|(n, ll)| ll == l && refnet(n)).map(|(n, _)| n.clone()).collect(); names.sort(); names }).collect();
    let at = |k: usize, x: f64, y: f64| -> u16 { let (c, r) = (((x - rx0) / rres).floor(), ((y - ry0) / rres).floor()); if c < 0.0 || r < 0.0 || c as usize >= rw || r as usize >= rh { 0 } else { refmask[k][r as usize * rw + c as usize] } };
    // transitions are not reference: a point within 0.3 mm of a via's or a through-hole's ring (any net: its antipad cuts
    // the planes) or of the trace's own pads is left out of the coverage, and a segment that is all transition gets no
    // verdict (its length is counted in `transitionMm`)
    let holes: Vec<(f64, f64, f64)> = b.vias.iter().chain(b.plated_holes.iter()).map(|v| (v.x, v.y, v.size / 2.0 + 0.3)).collect();
    let transition = |t: &crate::board::Track, x: f64, y: f64| -> bool { holes.iter().any(|(hx, hy, r)| (x - hx).hypot(y - hy) <= *r) || b.pads.iter().any(|p| p.net == t.net && (x - p.x).hypot(y - p.y) <= (p.area / std::f64::consts::PI).sqrt() + 0.3) };
    let samples = |t: &crate::board::Track| -> usize { let len = (t.b.0 - t.a.0).hypot(t.b.1 - t.a.1); let steps = (len / 0.25).ceil().max(1.0) as usize; (0..=steps).filter(|i| { let f = *i as f64 / steps as f64; !transition(t, t.a.0 + f * (t.b.0 - t.a.0), t.a.1 + f * (t.b.1 - t.a.1)) }).count() };
    // each other layer's coverage of a set of segments, and the net covering most of it
    let coverage = |segs: &[&crate::board::Track], k: usize| -> (f64, Option<String>) {
        let (mut n, mut hit) = (0usize, 0usize); let mut by: std::collections::BTreeMap<u16, usize> = Default::default();
        for t in segs { let (dx, dy) = (t.b.0 - t.a.0, t.b.1 - t.a.1); let len = dx.hypot(dy); let steps = (len / 0.25).ceil().max(1.0) as usize; let (nx, ny) = if len > 0.0 { (-dy / len, dx / len) } else { (0.0, 1.0) };
            for i in 0..=steps { let f = i as f64 / steps as f64; let (x, y) = (t.a.0 + f * dx, t.a.1 + f * dy); if transition(t, x, y) { continue; } n += 1;
                let ids = [at(k, x, y), at(k, x + nx * t.w / 2.0, y + ny * t.w / 2.0), at(k, x - nx * t.w / 2.0, y - ny * t.w / 2.0)];
                if ids.iter().all(|i| *i != 0) { hit += 1; *by.entry(ids[0]).or_default() += 1; } } }
        let net = by.iter().max_by_key(|x| x.1).map(|(id, _)| refnames[k][*id as usize - 1].clone());
        (if n == 0 { 0.0 } else { hit as f64 / n as f64 }, net)
    };
    // the references of a set of segments on layer li: on each side the nearest layer covering 90 % of them; a nearer layer
    // covering part of them makes the section mixed, nothing at all leaves it with no reference
    let reference = |segs: &[&crate::board::Track], li: usize| -> Result<(usize, Option<usize>, Value), (String, Value)> {
        let mut sides: Vec<Option<(usize, f64, Option<String>)>> = Vec::new(); let mut partial: Vec<String> = Vec::new(); let mut seen = Vec::new();
        for side in [-1.0f64, 1.0] {
            let mut cand: Vec<usize> = (0..real.len()).filter(|j| *j != li && (z[*j] - z[li]) * side > 0.0).collect(); cand.sort_by(|a, b2| (z[*a] - z[li]).abs().partial_cmp(&(z[*b2] - z[li]).abs()).unwrap());
            let mut chosen = None;
            for j in cand { let (c, net) = coverage(segs, j); seen.push(json!({"layer": real[j], "coverage": (c * 100.0).round() / 100.0, "net": net}));
                if c >= 0.9 { chosen = Some((j, c, net)); break; } if c >= 0.1 { partial.push(format!("{} covers {:.0} %", real[j], c * 100.0)); break; } }
            sides.push(chosen);
        }
        let mut refs: Vec<(usize, f64, Option<String>)> = sides.into_iter().flatten().collect(); refs.sort_by(|a, b2| (z[a.0] - z[li]).abs().partial_cmp(&(z[b2.0] - z[li]).abs()).unwrap());
        let rec = |r: &(usize, f64, Option<String>)| json!({"layer": real[r.0], "net": r.2, "distanceMm": ((z[r.0] - z[li]).abs() * 1000.0).round() / 1000.0, "coverage": (r.1 * 100.0).round() / 100.0});
        if !partial.is_empty() { return Err((format!("mixed reference: {}", partial.join(", ")), json!(seen))); }
        match refs.first() { None => Err(("no reference: no ground or supply copper under it on any layer".into(), json!(seen))), Some(r0) => Ok((r0.0, refs.get(1).map(|r| r.0), json!({"below": rec(r0), "above": refs.get(1).map(rec)}))) }
    };
    let signal = |n: &str| !n.is_empty() && !is_ground(n) && rail_from_name(n).is_none() && !spec.voltages.contains_key(n) && !spec.loads.contains_key(n) && !n.starts_with("unconnected-");
    let mut cache: std::collections::BTreeMap<String, Zr> = Default::default();
    let mut nets_out = Vec::new();
    let mut by_net: std::collections::BTreeMap<String, Vec<&crate::board::Track>> = Default::default();
    for t in b.tracks.iter().filter(|t| signal(&t.net)) { by_net.entry(t.net.clone()).or_default().push(t); }
    let targets = &spec.impedance;
    for (net, ts) in &by_net {
        // group the segments by layer, width and their own reference (a long net can cross several)
        let mut groups: std::collections::BTreeMap<(String, i64, String), Vec<&crate::board::Track>> = Default::default();
        let mut refs_of: std::collections::BTreeMap<String, Result<(usize, Option<usize>, Value), (String, Value)>> = Default::default();
        let mut transition_mm = 0.0;
        for t in ts { let Some(li) = real.iter().position(|l| *l == t.layer) else { continue }; if samples(t) == 0 { transition_mm += (t.a.0 - t.b.0).hypot(t.a.1 - t.b.1); continue; } let r = reference(&[*t], li); let key = match &r { Ok((a, b2, _)) => format!("{a}:{b2:?}"), Err((w, _)) => w.clone() }; refs_of.entry(key.clone()).or_insert(r); groups.entry((t.layer.clone(), (t.w * 1000.0).round() as i64, key)).or_default().push(t); }
        let mut rows = Vec::new();
        for ((layer, wum, key), segs) in &groups {
            let Some(li) = real.iter().position(|l| l == layer) else { continue };
            let len: f64 = segs.iter().map(|t| (t.a.0 - t.b.0).hypot(t.a.1 - t.b.1)).sum();
            match &refs_of[key] {
                Ok((bi, ai, rec)) => {
                    let ck = format!("{layer}:{wum}:{bi}:{ai:?}");
                    if !cache.contains_key(&ck) { cache.insert(ck.clone(), solve(&xs_for(li, *wum as f64 / 1000.0, None, *bi, *ai))); }
                    let zr = &cache[&ck];
                    rows.push(json!({"layer": layer, "widthMm": *wum as f64 / 1000.0, "lengthMm": r1(len), "status": "ok", "reference": rec, "z0": r1(zr.z0), "eeff": (zr.eeff * 100.0).round() / 100.0, "psPerMm": (zr.eeff.sqrt() / C0 * 1e9 * 100.0).round() / 100.0, "pfPerM": r1(zr.cpf_m), "nhPerM": r1(zr.lnh_m)}));
                }
                Err((why, seen)) => rows.push(json!({"layer": layer, "widthMm": *wum as f64 / 1000.0, "lengthMm": r1(len), "status": why, "layersUnder": seen, "z0": Value::Null})),
            }
        }
        if rows.is_empty() && transition_mm == 0.0 { continue; }
        let target = targets.get(net).copied();
        if let Some(tz) = target { let unresolved: f64 = rows.iter().filter(|r| r["z0"].is_null()).filter_map(|r| r["lengthMm"].as_f64()).sum();
            if unresolved > 0.0 { issues.push(Issue { severity: "low".into(), kind: "impedance".into(), what: format!("{net}: {unresolved:.1} mm of it has a mixed or missing reference, so its impedance is not evaluated there ({})", rows.iter().filter(|r| r["z0"].is_null()).filter_map(|r| r["status"].as_str()).collect::<std::collections::BTreeSet<_>>().into_iter().collect::<Vec<_>>().join("; ")), at: None, layer: None, net: Some(net.clone()), reference: None, fix: "Route it over continuous ground (or a plane), or check that stretch with a field solver that models the actual reference; review.impedance lists each stretch and the copper under it.".into(), action: Some(json!({"type": "review", "net": net, "why": "impedance reference"})) }); }
            for r in rows.iter().filter(|r| r["z0"].is_number()) { let z0 = r["z0"].as_f64().unwrap_or(0.0); if (z0 - tz).abs() > 0.1 * tz { issues.push(Issue { severity: "medium".into(), kind: "impedance".into(), what: format!("{}: {:.0} Ω on {} at {:.3} mm wide over {} {} {:.3} mm away, against a {:.0} Ω target", net, z0, r["layer"].as_str().unwrap_or(""), r["widthMm"].as_f64().unwrap_or(0.0), r["reference"]["below"]["layer"].as_str().unwrap_or("?"), r["reference"]["below"]["net"].as_str().unwrap_or(""), r["reference"]["below"]["distanceMm"].as_f64().unwrap_or(0.0), tz), at: None, layer: r["layer"].as_str().map(String::from), net: Some(net.clone()), reference: None, fix: "Change the trace width (wider lowers Z0) or the stackup's dielectric under it; the cross-section solve gives the width for the target.".into(), action: Some(json!({"type": "widen", "net": net, "why": "impedance target", "z0": z0, "target": tz})) }); } } }
        nets_out.push(json!({"net": net, "segments": rows, "targetOhm": target, "transitionMm": r1(transition_mm)}));
    }
    // pairs: by name, with the median spacing of their parallel segments on a shared layer
    let mut pairs = Vec::new(); let mut seen: std::collections::BTreeSet<String> = Default::default();
    for (net, ts) in &by_net {
        let Some(p) = partner(net) else { continue }; let Some(ps) = by_net.get(&p) else { continue };
        if seen.contains(net) || seen.contains(&p) { continue; } seen.insert(net.clone()); seen.insert(p.clone());
        let mut gaps_mm = Vec::new(); let mut layer = String::new(); let mut w = 0.0;
        for a in ts { for c in ps.iter().filter(|c| c.layer == a.layer && (c.w - a.w).abs() < 1e-6) {
            let (dx, dy) = (a.b.0 - a.a.0, a.b.1 - a.a.1); let (ex, ey) = (c.b.0 - c.a.0, c.b.1 - c.a.1);
            let la = dx.hypot(dy); let lc = ex.hypot(ey); if la < 0.5 || lc < 0.5 { continue; }
            if ((dx * ey - dy * ex) / (la * lc)).abs() > 0.05 { continue; }
            let mx = (c.a.0 + c.b.0) / 2.0 - a.a.0; let my = (c.a.1 + c.b.1) / 2.0 - a.a.1;
            let dist = (mx * dy - my * dx).abs() / la; let gap = dist - a.w;
            if gap > 0.0 && gap < 5.0 * a.w { gaps_mm.push(gap); layer = a.layer.clone(); w = a.w; }
        } }
        if gaps_mm.is_empty() { continue; }
        gaps_mm.sort_by(|x, y| x.partial_cmp(y).unwrap()); let gap = gaps_mm[gaps_mm.len() / 2];
        let Some(li) = real.iter().position(|l| *l == layer) else { continue };
        let segs: Vec<&crate::board::Track> = ts.iter().chain(ps.iter()).filter(|t| t.layer == layer && (t.w - w).abs() < 1e-6).copied().collect();
        let target0 = spec.impedance.get(&format!("{net}/{p}")).copied().map(|t| (t, "spec")).or_else(|| diff_target(net));
        let (bi, ai, refrec) = match reference(&segs, li) { Ok(x) => x, Err((why, seen)) => {
            pairs.push(json!({"pair": [net, p], "layer": layer, "widthMm": w, "gapMm": (gap * 1000.0).round() / 1000.0, "status": why, "layersUnder": seen, "zdiff": Value::Null, "targetOhm": target0.map(|t| t.0), "standard": target0.map(|t| t.1)}));
            if target0.is_some() { issues.push(Issue { severity: "low".into(), kind: "diff-impedance".into(), what: format!("{net}/{p}: impedance not evaluated on {layer}: {why}"), at: None, layer: Some(layer.clone()), net: Some(net.clone()), reference: None, fix: "Route the pair over continuous ground, or check it with a field solver that models the actual reference; review.impedance.pairs lists the copper under it.".into(), action: Some(json!({"type": "review", "pair": [net, p], "why": "impedance reference"})) }); }
            continue } };
        let zr = solve(&xs_for(li, w, Some(gap), bi, ai));
        let target = spec.impedance.get(&format!("{net}/{p}")).copied().map(|t| (t, "spec")).or_else(|| diff_target(net));
        let zd = zr.zdiff.unwrap_or(0.0);
        pairs.push(json!({"pair": [net, p], "layer": layer, "widthMm": w, "gapMm": (gap * 1000.0).round() / 1000.0, "status": "ok", "reference": refrec, "zdiff": r1(zd), "zodd": zr.zodd.map(r1), "zeven": zr.zeven.map(r1), "targetOhm": target.map(|t| t.0), "standard": target.map(|t| t.1)}));
        if let Some((tz, why)) = target { if (zd - tz).abs() > 0.1 * tz { issues.push(Issue { severity: "medium".into(), kind: "diff-impedance".into(), what: format!("{}/{}: {:.0} Ω differential on {} ({:.3} mm wide, {:.3} mm gap) over {} {} {:.3} mm away, against {:.0} Ω ({})", net, p, zd, layer, w, gap, refrec["below"]["layer"].as_str().unwrap_or("?"), refrec["below"]["net"].as_str().unwrap_or(""), refrec["below"]["distanceMm"].as_f64().unwrap_or(0.0), tz, why), at: None, layer: Some(layer.clone()), net: Some(net.clone()), reference: None, fix: format!("Adjust the pair's width and gap for {:.0} Ω (a narrower gap lowers it); keep the pair over unbroken ground.", tz), action: Some(json!({"type": "review", "pair": [net, p], "why": "differential impedance", "zdiff": r1(zd), "target": tz})) }); } }
    }
    json!({"ran": true, "method": "2-D Laplace field solve of each cross-section (grid, conjugate gradient), dielectric and air; solder mask 20 µm εr 3.6 on outer layers unless rf.noSolderMask. The reference is found under each segment: on each side the nearest layer whose ground or supply copper covers 90 % of it (centre and both edges, sampled every 0.25 mm); a nearer layer covering part of it marks the section mixed and a section with none has no reference, and neither gets an impedance", "nets": nets_out, "pairs": pairs})
}

#[cfg(test)] mod tests {
    use super::*;
    /// Hammerstad-Jensen microstrip with Wheeler's thickness correction (w, h, t in mm)
    fn hammerstad(w: f64, h: f64, t: f64, er: f64) -> f64 {
        let dw = t / std::f64::consts::PI * (1.0 + (2.0 * h / t).ln()); let u = (w + dw) / h;
        let z_air = if u >= 1.0 { 376.73 / (u + 1.393 + 0.667 * (u + 1.444).ln()) } else { 60.0 * (8.0 / u + u / 4.0).ln() };
        let eeff = (er + 1.0) / 2.0 + (er - 1.0) / 2.0 / (1.0 + 12.0 / u).sqrt();
        z_air / eeff.sqrt()
    }
    #[test] fn stripline_matches_closed_form() {
        // a centred stripline: w 0.2 mm, t 0.035 mm, b = 0.5 mm plane to plane, εr 4.3.
        // Wadell's thickness-corrected stripline: 60/√εr · ln(4b / (0.67π w (0.8 + t/w))) ≈ 45.8 Ω
        let x = Xs { w: 0.2e-3, t: 0.035e-3, gap: None, below: 0.2325e-3, above: Some(0.2325e-3), er: 4.3, mask: None };
        let z = solve(&x); eprintln!("stripline z0 {:.2} eeff {:.3}", z.z0, z.eeff); assert!(z.z0 > 44.5 && z.z0 < 47.5, "{}", z.z0); assert!((z.eeff - 4.3).abs() < 0.05, "{}", z.eeff);
    }
    #[test] fn microstrip_matches_closed_form() {
        // w 0.3 mm over 0.2 mm FR4 (εr 4.2), t 35 µm, no mask: Hammerstad-Jensen with Wheeler's thickness correction gives
        // 55.1 Ω and εeff 3.16 (the IPC-2141 formula in signals.rs says 54.0 but is known to be off by 5-10 %)
        let x = Xs { w: 0.3e-3, t: 0.035e-3, gap: None, below: 0.2e-3, above: None, er: 4.2, mask: None };
        let z = solve(&x); let hj = hammerstad(0.3, 0.2, 0.035, 4.2); eprintln!("microstrip z0 {:.2} (Hammerstad-Jensen {:.2}) eeff {:.3}", z.z0, hj, z.eeff);
        assert!((z.z0 - hj).abs() / hj < 0.03, "{} vs {}", z.z0, hj); assert!(z.eeff > 2.9 && z.eeff < 3.3, "{}", z.eeff);
    }
    #[test] fn pair_odd_below_single() {
        let x = Xs { w: 0.15e-3, t: 0.035e-3, gap: Some(0.15e-3), below: 0.1e-3, above: None, er: 4.2, mask: None };
        let z = solve(&x); eprintln!("pair z0 {:.1} zodd {:?} zeven {:?} zdiff {:?}", z.z0, z.zodd, z.zeven, z.zdiff); let zo = z.zodd.unwrap(); assert!(zo < z.z0 && z.zeven.unwrap() > z.z0, "{zo} {} {:?}", z.z0, z.zeven);
        assert!(z.zdiff.unwrap() > 70.0 && z.zdiff.unwrap() < 120.0, "{:?}", z.zdiff);
    }
}