RITHESH LED Display

"""Voltage drop on the 5 V supply and ground copper, solved on the real filled copper.

The filled +5V copper (front) and GND copper (back plane + front fill + tracks) are rasterised on a
0.4 mm grid. Each cell is a square of copper: sheet resistance of 1 oz (35 um) copper is 0.49 mOhm per
square at 25 C (0.57 at 65 C). Every LED draws the same current from its VDD pad and returns it through
its GND pad; the whole current comes in at the eFuse output (+5V) and returns to the USB-C ground.
The resulting sparse network is solved exactly (Kirchhoff), so narrow necks and gaps are included.

Usage: python3 power_analysis.py [total LED current in A ...]      (default: 1.88 and 2.9)
"""
import math
import os
import sys

import numpy as np
import pcbnew
from scipy.sparse import coo_matrix

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import make_board as mb  # noqa: E402

TO = pcbnew.ToMM
G = 0.4                         # grid, mm
RS = 0.49e-3                    # ohm per square, 1 oz copper, 25 C
RS_HOT = 0.57e-3                # at 65 C


def xy(v):
    return (TO(v.x) - mb.X0, TO(v.y) - mb.Y0)


def raster(board, netname, layers):
    """Boolean grid of cells covered by copper of `netname` on the given layers (zones + tracks + pads)."""
    nx, ny = int(mb.W / G) + 1, int(mb.H / G) + 1
    grid = np.zeros((len(layers), ny, nx), bool)
    xs = (np.arange(nx) + 0.5) * G
    ys = (np.arange(ny) + 0.5) * G
    for li, layer in enumerate(layers):
        polys = []
        for z in board.Zones():
            if z.GetIsRuleArea() or z.GetNetname() != netname or not z.IsOnLayer(layer):
                continue
            fp = z.GetFilledPolysList(layer)
            for i in range(fp.OutlineCount()):
                o = fp.Outline(i)
                outer = [xy(o.CPoint(k)) for k in range(o.PointCount())]
                holes = []
                for h in range(fp.HoleCount(i)):
                    hh = fp.Hole(i, h)
                    holes.append([xy(hh.CPoint(k)) for k in range(hh.PointCount())])
                polys.append((outer, holes))
        import shapely
        from shapely.geometry import Polygon
        X, Y = np.meshgrid(xs, ys)
        m = np.zeros(X.shape, bool)
        for outer, holes in polys:
            m |= shapely.contains_xy(Polygon(outer, holes), X, Y)
        grid[li] = m
        for t in board.GetTracks():
            if t.GetNetname() != netname or isinstance(t, pcbnew.PCB_VIA) or t.GetLayer() != layer:
                continue
            (x0, y0), (x1, y1) = xy(t.GetStart()), xy(t.GetEnd())
            w = TO(t.GetWidth())
            L = max(math.hypot(x1 - x0, y1 - y0), 1e-6)
            for s in np.linspace(0, 1, int(L / (G / 2)) + 2):
                cx, cy = x0 + (x1 - x0) * s, y0 + (y1 - y0) * s
                r = w / 2
                i0, i1 = int((cy - r) / G), int((cy + r) / G)
                j0, j1 = int((cx - r) / G), int((cx + r) / G)
                grid[li, max(i0, 0):i1 + 1, max(j0, 0):j1 + 1] = True
    return grid


def vias_of(board, netname):
    return [xy(t.GetPosition()) for t in board.GetTracks() if isinstance(t, pcbnew.PCB_VIA) and t.GetNetname() == netname]


def solve(grid, vias, sources, sinks, rs):
    """grid: (L, ny, nx) bool. vias join layers (0.8 mOhm each). sources: [(layer, x, y)] held at 0 V;
    sinks: [(layer, x, y, amps)]. Returns potential map (V, negative = drop) and per-sink voltage."""
    L, ny, nx = grid.shape
    idx = -np.ones(grid.shape, int)
    idx[grid] = np.arange(grid.sum())
    n = int(grid.sum())
    rows, cols, vals = [], [], []
    g = 1.0 / rs
    for li in range(L):
        a = idx[li]
        for (da, db) in ((0, 1), (1, 0)):
            p = a[:ny - da, :nx - db]
            q = a[da:, db:]
            ok = (p >= 0) & (q >= 0)
            pp, qq = p[ok], q[ok]
            rows += [pp, qq, pp, qq]
            cols += [pp, qq, qq, pp]
            vals += [np.full(pp.size, g)] * 2 + [np.full(pp.size, -g)] * 2

    def cell(li, x, y):
        i, j = int(y / G), int(x / G)
        best = None
        for r in range(0, 8):
            for di in range(-r, r + 1):
                for dj in range(-r, r + 1):
                    ii, jj = i + di, j + dj
                    if 0 <= ii < ny and 0 <= jj < nx and idx[li, ii, jj] >= 0:
                        return idx[li, ii, jj]
        return best
    if L > 1:
        gv = 1.0 / 0.8e-3
        for (x, y) in vias:
            a, b = cell(0, x, y), cell(1, x, y)
            if a is None or b is None:
                continue
            rows += [np.array([a, b, a, b])]
            cols += [np.array([a, b, b, a])]
            vals += [np.array([gv, gv, -gv, -gv])]
    A = coo_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))), shape=(n, n)).tocsr()
    rhs = np.zeros(n)
    sink_cells = []
    for (li, x, y, amps) in sinks:
        c = cell(li, x, y)
        sink_cells.append(c)
        if c is not None:
            rhs[c] -= amps
    fixed = [cell(li, x, y) for (li, x, y) in sources]
    fixed = [f for f in fixed if f is not None]
    # the source cells are tied to 0 V through a very small resistance (1 micro-ohm)
    from scipy.sparse import diags
    from scipy.sparse.linalg import cg
    d = np.full(n, 1e-6)        # tiny leak so isolated specks of copper cannot make the system singular
    for f in fixed:
        d[f] += 1e6
    A = (A + diags(d)).tocsr()
    M = diags(1.0 / A.diagonal())
    v, info = cg(A, rhs, tol=1e-10, maxiter=40000, M=M)
    if info != 0:
        print("warning: solver did not fully converge", info)
    return v, [v[c] if c is not None else float("nan") for c in sink_cells], idx


# Power-entry copper, section by section (widths / lengths as routed in route_board.py).
# (name, width mm, length mm, share of the input current: (low, high))
RIGHT, LEFT = 1.55 / 0.5, 0.75 / 0.3 + 9.1 / 0.6           # squares: right-pad stub vs left-pad path
SECTIONS = [
    ("VBUS right pad stub (A9/B4)", 0.5, 1.55, None),
    ("VBUS left pad stub (A4/B9)", 0.3, 0.75, "left"),
    ("VBUS left pad, back-layer link", 0.6, 9.1, "left"),
    ("VBUS main, 0.8 mm part", 0.8, 1.1, 1.0),
    ("VBUS main, 1.2 mm to TVS", 1.2, 9.4, 1.0),
    ("TVS -> input capacitor", 0.8, 1.6, 1.0),
    ("input capacitor -> eFuse IN", 0.6, 1.45, 1.0),
    ("eFuse OUT pin neck", 0.25, 1.2, 1.0),
    ("5 V out -> +5V fill", 1.2, 3.05, 1.0),
]


def section_table(currents=(2.2, 2.9), contact_mohm=(0.0, 20.0)):
    """Per-section drop, dissipation and IPC-2221 temperature rise (outer 1 oz, long-track rule: short
    sections run cooler because heat spreads into pads and wide copper)."""
    def rise(i, w):
        a = w / 0.0254 * 1.378                          # cross-section in mil^2
        return (i / (0.048 * a ** 0.725)) ** (1 / 0.44)
    # how the input splits between the two VBUS pads: without and with ~20 mOhm contact resistance per pad
    splits = []
    for rc in contact_mohm:
        rr, rl = RIGHT * 0.49 + rc, LEFT * 0.49 + 1.6 + rc      # mOhm (two vias on the left path, 0.8 each)
        splits.append(rl / (rr + rl))
    lo, hi = min(splits), max(splits)
    print(f"VBUS current share: right pad {lo * 100:.0f}-{hi * 100:.0f} %, left pad {100 - hi * 100:.0f}-{100 - lo * 100:.0f} %")
    for I in currents:
        print(f"-- input current {I:.1f} A --")
        for name, w, L, share in SECTIONS:
            if share is None:
                i = I * hi
            elif share == "left":
                i = I * (1 - lo)
            else:
                i = I * share
            r = RS * L / w
            p = i * i * r
            # short neck held at both ends by pads / wide copper: peak rise = P * L / (8 k A), copper k = 390
            neck = p * (L * 1e-3) / (8 * 390 * (w * 1e-3) * 35e-6)
            est = f"short-neck est. {neck:4.1f} C" if L < 5 else "long track: use IPC"
            print(f"  {name:32s} {w:.2f} x {L:4.2f} mm  I {i:.2f} A  drop {i * r * 1000:5.1f} mV  "
                  f"heat {p * 1000:5.1f} mW  IPC rise {rise(i, w):5.1f} C  {est}")
        print(f"  (eFuse switch 34 mOhm typ: drop {I * 34:.0f} mV, heat {I * I * 34:.0f} mW; "
              f"USB-C VBUS contacts rated 5 A for the connector = 2.5 A per pad pair)")


def main():
    if len(sys.argv) > 1 and sys.argv[1] == "sections":
        return section_table()
    currents = [float(a) for a in sys.argv[1:]] or [1.88, 2.9]
    board = pcbnew.LoadBoard(mb.PCB)
    fps = {f.GetReference(): f for f in board.GetFootprints()}
    leds = [f"LED{n}" for n in range(1, 108)]

    def pad(ref, num):
        return xy([p for p in fps[ref].Pads() if p.GetNumber() == num][0].GetPosition())

    v5 = raster(board, "+5V", [pcbnew.F_Cu])
    gnd = raster(board, "GND", [pcbnew.F_Cu, pcbnew.B_Cu])
    gvias = vias_of(board, "GND")
    src5 = [(0,) + pad("U2", "5")]
    srcg = [(0,) + pad("J1", "A1"), (0,) + pad("J1", "A12"), (1,) + pad("J1", "A1"), (1,) + pad("J1", "A12")]
    print(f"+5V copper: {v5.sum() * G * G:.0f} mm^2 on the front; GND copper: front {gnd[0].sum() * G * G:.0f} "
          f"mm^2, back {gnd[1].sum() * G * G:.0f} mm^2, {len(gvias)} ground vias")
    results = {}
    for I in currents:
        per = I / len(leds)
        _, vv, _ = solve(v5, [], src5, [(0,) + pad(r, "1") + (per,) for r in leds], RS)
        _, vg, _ = solve(gnd, gvias, srcg, [(0,) + pad(r, "3") + (per,) for r in leds], RS)
        worst = max(range(len(leds)), key=lambda k: (-vv[k]) + (-vg[k]))
        tot = [(-vv[k]) + (-vg[k]) for k in range(len(leds))]
        results[I] = (max(-x for x in vv), max(-x for x in vg), max(tot), leds[worst])
        print(f"total LED current {I:.2f} A ({per * 1000:.1f} mA per LED): worst +5V drop {max(-x for x in vv) * 1000:.0f} mV, "
              f"worst GND rise {max(-x for x in vg) * 1000:.0f} mV, worst total {max(tot) * 1000:.0f} mV at {leds[worst]} "
              f"(at 65 C copper: {max(tot) * 1000 * RS_HOT / RS:.0f} mV); average {1000 * sum(tot) / len(tot):.0f} mV")
    return results


if __name__ == "__main__":
    main()