Adom Fields

//! Partial-element equivalent circuit (PEEC) loop inductance from the board's copper (0.10.2, used by the pdn module).
//! The copper of two nets in a window is cut into cells on every copper layer; neighbouring cells of one net are
//! bars (x- or y-directed branches) and vias are z-directed branches between layers. Each branch has a partial self
//! inductance and a partial mutual inductance with every parallel branch (Ruehli: parallel filaments, the self at the
//! bar's geometric mean distance). At the frequencies where mounting inductance matters, jωL dwarfs the resistance,
//! so the inductive network is solved alone: [L Aᵀ; A 0][I; −φ] = [0; J], 1 A in at the port; the port's inductance
//! is its potential difference. A part's two pins can be shorted (the chip's rail pin to its ground pin, where the
//! decoupling loop closes inside the package). Validated against the 2-D field solve of a strip over a plane (tests).
use crate::board::Board;

const MU0_4PI: f64 = 1e-7; // H/m

/// mutual inductance of two parallel filaments along one axis: [a1, a2] and [b1, b2] (m), d apart (m)
pub fn m_par(a1: f64, a2: f64, b1: f64, b2: f64, d: f64) -> f64 {
    let f = |z: f64| z * (z / d).asinh() - (z * z + d * d).sqrt();
    MU0_4PI * (f(a2 - b1) - f(a2 - b2) - f(a1 - b1) + f(a1 - b2))
}
/// partial self inductance of a bar of length l, width w, thickness t (m): a filament at the cross-section's GMD
pub fn l_self(l: f64, w: f64, t: f64) -> f64 { m_par(0.0, l, 0.0, l, 0.2235 * (w + t)) }

#[derive(Clone, Copy, PartialEq)] pub enum Dir { X, Y, Z }
#[derive(Clone)] pub struct Branch { pub n1: usize, pub n2: usize, pub dir: Dir, pub a: (f64, f64, f64), pub len: f64, pub w: f64, pub t: f64 }

/// the loop inductance seen at a port (node p, node q) of a branch network: dense LU on the MNA system
pub fn solve(_nn: usize, br: &[Branch], shorts: &[(usize, usize)], p: usize, q: usize) -> Option<f64> {
    // only the nodes a branch or a short touches are unknowns (empty cells would leave the system singular)
    // keep only the piece of copper the port is on: floating islands leave their potentials undefined
    let mut adj: std::collections::BTreeMap<usize, Vec<usize>> = Default::default();
    for (a, b) in br.iter().map(|b| (b.n1, b.n2)).chain(shorts.iter().copied()) { adj.entry(a).or_default().push(b); adj.entry(b).or_default().push(a); }
    let mut keep: std::collections::BTreeSet<usize> = Default::default(); let mut stack = vec![p];
    while let Some(u) = stack.pop() { if keep.insert(u) { for v in adj.get(&u).into_iter().flatten() { if !keep.contains(v) { stack.push(*v); } } } }
    if !keep.contains(&q) { return None; }
    let br: Vec<Branch> = br.iter().filter(|b| keep.contains(&b.n1)).cloned().collect();
    let shorts: Vec<(usize, usize)> = shorts.iter().filter(|s| keep.contains(&s.0)).copied().collect();
    let mut used: std::collections::BTreeMap<usize, usize> = Default::default();
    for (a, b) in br.iter().map(|b| (b.n1, b.n2)).chain(shorts.iter().copied()) { let k = used.len(); used.entry(a).or_insert(k); let k = used.len(); used.entry(b).or_insert(k); }
    let (Some(&p), Some(&q)) = (used.get(&p), used.get(&q)) else { return None };
    let br: Vec<Branch> = br.iter().map(|b| Branch { n1: used[&b.n1], n2: used[&b.n2], ..b.clone() }).collect();
    let shorts: Vec<(usize, usize)> = shorts.iter().map(|(a, b)| (used[a], used[b])).collect();
    let (br, shorts, nn) = (&br[..], &shorts[..], used.len());
    let nb = br.len(); let ns = shorts.len();
    // unknowns: branch currents (nb), short currents (ns), node potentials except node q (ground) (nn - 1)
    let idx = |n: usize| -> Option<usize> { if n == q { None } else if n < q { Some(nb + ns + n) } else { Some(nb + ns + n - 1) } };
    let m = nb + ns + nn - 1;
    let mut a = vec![0.0f64; m * m]; let mut rhs = vec![0.0f64; m];
    for i in 0..nb { for j in i..nb {
        let (bi, bj) = (&br[i], &br[j]); if bi.dir != bj.dir { continue; }
        let v = if i == j { l_self(bi.len, bi.w, bi.t) } else {
            // the axis coordinate and the perpendicular distance of two parallel bars
            let (s1, s2, perp) = match bi.dir { Dir::X => (bi.a.0, bj.a.0, ((bi.a.1 - bj.a.1).powi(2) + (bi.a.2 - bj.a.2).powi(2)).sqrt()), Dir::Y => (bi.a.1, bj.a.1, ((bi.a.0 - bj.a.0).powi(2) + (bi.a.2 - bj.a.2).powi(2)).sqrt()), Dir::Z => (bi.a.2, bj.a.2, ((bi.a.0 - bj.a.0).powi(2) + (bi.a.1 - bj.a.1).powi(2)).sqrt()) };
            let d = perp.max(0.2235 * (bi.w + bi.t));
            m_par(s1, s1 + bi.len, s2, s2 + bj.len, d)
        };
        a[i * m + j] = v; a[j * m + i] = v;
    } }
    // KVL rows: L·I − (φ(n1) − φ(n2)) = 0 → L·I + φ(n2) − φ(n1); signs folded into the incidence
    let inc = |a: &mut Vec<f64>, row: usize, n1: usize, n2: usize| { if let Some(k) = idx(n1) { a[row * m + k] -= 1.0; } if let Some(k) = idx(n2) { a[row * m + k] += 1.0; } };
    for (i, b) in br.iter().enumerate() { inc(&mut a, i, b.n1, b.n2); }
    for (s, (n1, n2)) in shorts.iter().enumerate() { inc(&mut a, nb + s, *n1, *n2); }
    // KCL rows (one per node except q): currents leaving through branches and shorts = injected
    for (i, b) in br.iter().enumerate() { if let Some(r) = idx(b.n1) { a[r * m + i] += 1.0; } if let Some(r) = idx(b.n2) { a[r * m + i] -= 1.0; } }
    for (s, (n1, n2)) in shorts.iter().enumerate() { let c = nb + s; if let Some(r) = idx(*n1) { a[r * m + c] += 1.0; } if let Some(r) = idx(*n2) { a[r * m + c] -= 1.0; } }
    let rp = idx(p)?; rhs[rp] = 1.0; // 1 A into p, out at q
    let x = lu_solve(&mut a, &mut rhs, m)?;
    Some(x[rp])
}

fn lu_solve(a: &mut [f64], b: &mut [f64], m: usize) -> Option<Vec<f64>> {
    for k in 0..m {
        let piv = (k..m).max_by(|i, j| a[i * m + k].abs().partial_cmp(&a[j * m + k].abs()).unwrap())?;
        if a[piv * m + k].abs() < 1e-30 { return None; }
        if piv != k { for c in 0..m { a.swap(k * m + c, piv * m + c); } b.swap(k, piv); }
        let d = a[k * m + k];
        let (top, rest) = a.split_at_mut((k + 1) * m); let rowk = &top[k * m..(k + 1) * m];
        for i in 0..(m - k - 1) {
            let row = &mut rest[i * m..(i + 1) * m]; let f = row[k] / d; if f == 0.0 { continue; }
            row[k] = 0.0; for c in (k + 1)..m { row[c] -= f * rowk[c]; } b[k + 1 + i] -= f * b[k];
        }
    }
    let mut x = vec![0.0; m];
    for k in (0..m).rev() { let mut s = b[k]; for c in (k + 1)..m { s -= a[k * m + c] * x[c]; } x[k] = s / a[k * m + k]; }
    Some(x)
}

/// A grid of copper cells: for each layer and cell, the net there (or none). Builds the branch network.
pub struct Grid { pub x0: f64, pub y0: f64, pub d: f64, pub w: usize, pub h: usize, pub z: Vec<f64>, pub t: f64, pub cell: Vec<Vec<Option<u16>>> }
impl Grid {
    pub fn node(&self, k: usize, r: usize, c: usize) -> usize { (k * self.h + r) * self.w + c }
    /// branches between neighbouring cells of one net on one layer, and z branches at the given via points (layer pairs)
    pub fn network(&self, vias: &[(f64, f64, u16)]) -> Vec<Branch> {
        let mut br = Vec::new(); let (d, t) = (self.d * 1e-3, self.t);
        for k in 0..self.cell.len() { for r in 0..self.h { for c in 0..self.w {
            let Some(n) = self.cell[k][r * self.w + c] else { continue };
            let (x, y, z) = ((self.x0 + (c as f64 + 0.5) * self.d) * 1e-3, (self.y0 + (r as f64 + 0.5) * self.d) * 1e-3, self.z[k] * 1e-3);
            if c + 1 < self.w && self.cell[k][r * self.w + c + 1] == Some(n) { br.push(Branch { n1: self.node(k, r, c), n2: self.node(k, r, c + 1), dir: Dir::X, a: (x, y, z), len: d, w: d, t }); }
            if r + 1 < self.h && self.cell[k][(r + 1) * self.w + c] == Some(n) { br.push(Branch { n1: self.node(k, r, c), n2: self.node(k, r + 1, c), dir: Dir::Y, a: (x, y, z), len: d, w: d, t }); }
        } } }
        for (vx, vy, n) in vias {
            let (c, r) = (((vx - self.x0) / self.d) as isize, ((vy - self.y0) / self.d) as isize); if c < 0 || r < 0 || c as usize >= self.w || r as usize >= self.h { continue; }
            let (r, c) = (r as usize, c as usize);
            let ks: Vec<usize> = (0..self.cell.len()).filter(|k| self.cell[*k][r * self.w + c] == Some(*n)).collect();
            for p in ks.windows(2) { let (za, zb) = (self.z[p[0]], self.z[p[1]]); br.push(Branch { n1: self.node(p[0], r, c), n2: self.node(p[1], r, c), dir: Dir::Z, a: ((self.x0 + (c as f64 + 0.5) * self.d) * 1e-3, (self.y0 + (r as f64 + 0.5) * self.d) * 1e-3, za.min(zb) * 1e-3), len: (za - zb).abs() * 1e-3, w: 0.3e-3, t: 0.3e-3 }); }
        }
        br
    }
}

/// Rasterise the board's copper for the named nets inside a window (mm) at pitch d (mm).
pub fn grid_from_board(b: &Board, nets: &[String], win: [f64; 4], d: f64, z_layers: &[(String, f64)], t_mm: f64) -> Grid {
    let w = ((win[2] - win[0]) / d).ceil() as usize; let h = ((win[3] - win[1]) / d).ceil() as usize;
    let inside = |p: (f64, f64), poly: &[(f64, f64)]| { let mut c = false; let n = poly.len(); for i in 0..n { let (a, bb) = (poly[i], poly[(i + n - 1) % n]); if (a.1 > p.1) != (bb.1 > p.1) && p.0 < (bb.0 - a.0) * (p.1 - a.1) / (bb.1 - a.1) + a.0 { c = !c; } } c };
    let seg = |p: (f64, f64), a: (f64, f64), bb: (f64, f64)| { let (dx, dy) = (bb.0 - a.0, bb.1 - a.1); let l2 = dx * dx + dy * dy; let t = if l2 == 0.0 { 0.0 } else { (((p.0 - a.0) * dx + (p.1 - a.1) * dy) / l2).clamp(0.0, 1.0) }; (p.0 - a.0 - t * dx).hypot(p.1 - a.1 - t * dy) };
    let mut cell: Vec<Vec<Option<u16>>> = vec![vec![None; w * h]; z_layers.len()];
    for (k, (ln, _)) in z_layers.iter().enumerate() { for r in 0..h { for c in 0..w {
        let p = (win[0] + (c as f64 + 0.5) * d, win[1] + (r as f64 + 0.5) * d);
        for (ni, n) in nets.iter().enumerate() {
            let hit = b.polys.iter().chain(b.pad_polys.iter()).any(|q| q.net == *n && (q.layer == *ln || q.layer == "*.Cu") && q.pts.len() > 2 && inside(p, &q.pts))
                || b.tracks.iter().any(|t| t.net == *n && t.layer == *ln && seg(p, t.a, t.b) <= t.w / 2.0);
            if hit { cell[k][r * w + c] = Some(ni as u16); break; }
        }
    } } }
    // copper narrower than a cell must not vanish (a 0.25 mm trace at 0.6 mm cells would break the loop): every cell a
    // track passes through, and the cell under every pad and via, belongs to its net unless another net holds it
    let mut mark = |k: usize, x: f64, y: f64, ni: u16, cell: &mut Vec<Vec<Option<u16>>>| { let (c, r) = (((x - win[0]) / d).floor(), ((y - win[1]) / d).floor()); if c >= 0.0 && r >= 0.0 && (c as usize) < w && (r as usize) < h { let i = r as usize * w + c as usize; if cell[k][i].is_none() { cell[k][i] = Some(ni); } } };
    for (k, (ln, _)) in z_layers.iter().enumerate() { for (ni, n) in nets.iter().enumerate() {
        for t in b.tracks.iter().filter(|t| t.net == *n && t.layer == *ln) { let len = (t.b.0 - t.a.0).hypot(t.b.1 - t.a.1); let steps = (len / (d / 3.0)).ceil().max(1.0) as usize; for i in 0..=steps { let f = i as f64 / steps as f64; mark(k, t.a.0 + f * (t.b.0 - t.a.0), t.a.1 + f * (t.b.1 - t.a.1), ni as u16, &mut cell); } }
        for p in b.pads.iter().filter(|p| p.net == *n && (p.layers.iter().any(|l| l == ln || l == "*.Cu"))) { mark(k, p.x, p.y, ni as u16, &mut cell); }
        for v in b.vias.iter().chain(b.plated_holes.iter()).filter(|v| v.net == *n) { mark(k, v.x, v.y, ni as u16, &mut cell); }
    } }
    Grid { x0: win[0], y0: win[1], d, w, h, z: z_layers.iter().map(|x| x.1).collect(), t: t_mm * 1e-3, cell }
}

/// The loop inductance (nH) between two nets: 1 A in at port_a (net a), out at port_b (net b), with a short from
/// short_a (net a) to short_b (net b) closing the loop, over the copper in the window. None when the loop is not
/// closed by the copper in the window, or the window would take more than max_branches.
pub fn loop_nh(b: &Board, layers: &[(String, f64)], t_mm: f64, nets: [&str; 2], port: [(f64, f64); 2], short: [(f64, f64); 2], win: [f64; 4], pitch: f64, max_branches: usize, resolve_mm: f64) -> Result<(f64, usize, f64), String> {
    let nets_v = vec![nets[0].to_string(), nets[1].to_string()];
    let vias: Vec<(f64, f64, u16)> = b.vias.iter().chain(b.plated_holes.iter()).filter_map(|v| Some((v.x, v.y, nets_v.iter().position(|q| *q == v.net)? as u16))).collect();
    // the dense solve grows as the cube of the branches: coarsen the cells until the window fits the budget (the
    // rasterising is cheap; at 0.5 mm the buck's loops are within about 5 % of their converged values)
    let (mut pitch, mut tries) = (pitch, 0);
    let (g, br) = loop {
        let g = grid_from_board(b, &nets_v, win, pitch, layers, t_mm); let br = g.network(&vias);
        if br.len() <= max_branches || tries >= 4 { break (g, br); }
        pitch *= ((br.len() as f64 / max_branches as f64).sqrt() * 1.05).max(1.1); tries += 1;
    };
    if br.is_empty() { return Err("no copper of the two nets in the window".into()); }
    // a plane closer than the cells can resolve gives a wrong loop (the validation: −37 % at cells 2.5 times the dielectric, −17 % at
    // 1.25, a few % at 1 or under): a thin four-layer stackup needs a faster solver than this dense one
    if pitch > resolve_mm { return Err(format!("the nearest plane is {resolve_mm:.2} mm away, closer than the {pitch:.2} mm cells this solve can afford in the window can resolve (thin stackups need a faster solver, not yet built)")); }
    if br.len() > max_branches || pitch > 1.0 { return Err(format!("the window needs {} branches even at {pitch:.2} mm cells", br.len())); }
    // a point's node: the cell under it on the first layer (top down) where that net's copper is
    let at = |p: (f64, f64), ni: u16| -> Option<usize> { let (c, r) = (((p.0 - g.x0) / g.d) as isize, ((p.1 - g.y0) / g.d) as isize); if c < 0 || r < 0 || c as usize >= g.w || r as usize >= g.h { return None; } let (r, c) = (r as usize, c as usize); (0..g.cell.len()).find(|k| g.cell[*k][r * g.w + c] == Some(ni)).map(|k| g.node(k, r, c)) };
    let (Some(p), Some(q)) = (at(port[0], 0), at(port[1], 1)) else { return Err("a capacitor pad is not on its net's copper in the grid".into()) };
    let (Some(sa), Some(sz)) = (at(short[0], 0), at(short[1], 1)) else { return Err("a load or ground pin is not on its net's copper in the grid".into()) };
    let l = solve(2 * g.w * g.h * layers.len(), &br, &[(sa, sz)], p, q).ok_or("the loop does not close through the copper in the window")?;
    if !l.is_finite() || l <= 0.0 { return Err("the solve gave no inductance".into()); }
    Ok((l * 1e9, br.len(), pitch))
}

#[cfg(test)]
/// Validation: a strip (wc cells wide, lc cells long) over a full plane h below, shorted at the far end, port at the near
/// end. Returns (PEEC nH, 2-D field-solve nH from the strip's L per metre over the same length).
pub fn validate_strip(d_mm: f64, wc: usize, lc: usize, h_mm: f64, plane_cells: usize) -> (f64, f64) {
    let w = lc + 2; let hh = plane_cells.max(wc + 2);
    let r0 = (hh - wc) / 2;
    let mut cell = vec![vec![None; w * hh]; 2];
    for r in 0..hh { for c in 0..w { cell[1][r * w + c] = Some(0u16); } }
    for r in r0..r0 + wc { for c in 1..1 + lc { cell[0][r * w + c] = Some(0u16); } }
    let g = Grid { x0: 0.0, y0: 0.0, d: d_mm, w, h: hh, z: vec![h_mm, 0.0], t: 35e-6, cell };
    // the far-end short: one via per strip row at the last strip column
    let vias: Vec<(f64, f64, u16)> = (r0..r0 + wc).map(|r| ((lc as f64 + 0.5) * d_mm, (r as f64 + 0.5) * d_mm, 0u16)).collect();
    let br = g.network(&vias);
    // the port: the near-end strip column against the plane under it, all rows tied by shorts to one node each side
    let mut shorts = Vec::new(); let p = g.node(0, r0, 1); let q = g.node(1, r0, 1);
    for r in r0 + 1..r0 + wc { shorts.push((p, g.node(0, r, 1))); shorts.push((q, g.node(1, r, 1))); }
    let nn = 2 * w * hh;
    let l = solve(nn, &br, &shorts, p, q).unwrap_or(f64::NAN) * 1e9;
    let x = crate::xsec::Xs { w: wc as f64 * d_mm * 1e-3, t: 35e-6, gap: None, below: h_mm * 1e-3, above: None, er: 1.0, mask: None };
    let z = crate::xsec::solve(&x); let lpm = z.z0 / 299_792_458.0;
    (l, lpm * (lc as f64 - 1.0) * d_mm * 1e-3 * 1e9)
}

#[cfg(test)] mod tests {
    use super::*;
    #[test] fn wire_mutual_closed_form() {
        let (l, d) = (0.01, 0.001); let m = m_par(0.0, l, 0.0, l, d);
        let want = 2e-7 * l * ((l / d).asinh() - (1.0 + d * d / (l * l)).sqrt() + d / l);
        assert!((m - want).abs() / want < 1e-12);
    }
    #[test] fn strip_over_plane_matches_the_field_solve() {
        // a 1 mm strip 0.5 mm over a plane, shorted at the far end: within 3 % of the 2-D solve's L per metre × length
        let (pe, fs) = validate_strip(0.5, 2, 20, 0.5, 12); assert!((pe / fs - 1.0).abs() < 0.03, "PEEC {pe} nH, field solve {fs} nH");
    }
}