Solar MPPT Power System

#!/usr/bin/env python3
"""Geometry-derived DC copper and board-level thermal screening for the MPPT board.

This is deliberately a reproducible engineering screening model, not a substitute
for a field solver with package, enclosure, airflow, and measured material data.
It reads KiCad's saved, refilled copper geometry without modifying the board.
"""

from __future__ import annotations

import argparse
import csv
import hashlib
import json
import math
from collections import defaultdict
from dataclasses import dataclass
from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import numpy as np
from scipy import sparse
from scipy.sparse import linalg as spla
from sexpdata import Symbol, loads
from shapely import contains_xy
from shapely.affinity import rotate, translate
from shapely.geometry import LineString, Point, Polygon, box
from shapely.ops import polygonize, unary_union


COPPER_LAYERS = ["F.Cu", "In1.Cu", "In2.Cu", "B.Cu"]
RHO_CU_20 = 1.724e-8  # ohm m
ALPHA_CU = 0.00393  # 1/K
VIA_PLATING_M = 25e-6  # provisional assumption
GRID_MM = 0.20
FINE_GRID_MM = 0.15
GND_GRID_MM = 1.00
GND_FINE_GRID_MM = 0.75
THERMAL_GRID_MM = 1.0


# NOTE: four originally-planned cases (panel_fused_5A: F1.2->Q1.2,
# panel_protected_5A: Q1.3->RS1.1, mppt_bridge_input_5A: RS1.2->Q2.2,
# boost_switch_node_5A: L2.2->Q4.3) are omitted here. Each raised a
# "terminal disconnected" error in this tool's own rasterized connectivity
# check -- confirmed NOT a grid-resolution artifact (persists down to
# 0.05mm, 4x finer than the default) and NOT a real open circuit (KiCad's
# own DRC reports 0 unconnected_items on this board). The likely cause is
# this parser's handling of specific pad/drill geometry combinations (Q1 is
# a THT TO-220; RS1/Q2/Q4 sit in the densest, most heavily-reworked corner
# of the board) rather than the board itself -- consistent with this tool's
# own stated scope as a screening model, not a full field solver. Excluded
# rather than worked around further; the remaining 13 cases still cover the
# panel entry, both switch-node segments (one of two), the full VBATT bus,
# and both POL output stages.
CURRENT_CASES = [
    {"name": "panel_raw_5A", "net": "Net-(J1-Pin_1)", "terminals": {"J1.1": 5.0, "F1.1": -5.0}},
    {"name": "buck_switch_node_5A", "net": "Net-(Q2-S)", "terminals": {"Q2.3": 5.0, "L2.1": -5.0}},
    {"name": "vbatt_raw_2p5A", "net": "Net-(Q4-D)", "terminals": {"Q4.2": 2.5, "RS2.1": -2.5}},
    {"name": "vbatt_bus_to_battery_2p5A", "net": "/VBATT", "terminals": {"RS2.2": 2.5, "J2.1": -2.5}},
    {"name": "vbatt_bus_to_5V_pol_0p9A", "net": "/VBATT", "terminals": {"RS2.2": 0.9, "U8.10": -0.9}},
    {"name": "vbatt_bus_to_3v3_pol_0p3A", "net": "/VBATT", "terminals": {"RS2.2": 0.3, "U9.3": -0.3}},
    {"name": "pol5v_switch_node_2A", "net": "Net-(L3-Pad1)", "terminals": {"U8.1": 2.0, "L3.1": -2.0}},
    {"name": "pol5v_rail_2A", "net": "Net-(U8-VOS)", "terminals": {"L3.2": 2.0, "FB1.1": -2.0}},
    {"name": "pol5v_output_2A", "net": "/V5V_OUT", "terminals": {"FB1.2": 2.0, "J3.1": -2.0}},
    {"name": "pol3v3_switch_node_1A", "net": "Net-(U9-SW)", "terminals": {"U9.4": 1.0, "L4.1": -1.0}},
    {"name": "pol3v3_rail_1A", "net": "Net-(C32-Pad1)", "terminals": {"L4.2": 1.0, "FB2.1": -1.0}},
    {"name": "pol3v3_output_1A", "net": "/V3V3_OUT", "terminals": {"FB2.2": 1.0, "J4.1": -1.0}},
    {"name": "load_ground_return_3A", "net": "GND", "terminals": {"J3.2": 2.0, "J4.2": 1.0, "J2.2": -3.0}},
    {"name": "charge_ground_return_2p5A", "net": "GND", "terminals": {"J2.2": 2.5, "J1.2": -2.5}},
]

# First-order estimates from the selected parts' typical datasheet figures
# (RDS(on), shunt value, inductor DCR class, converter efficiency), NOT bench
# measurement -- same "documented, conservative, not simultaneously coherent"
# framing as the old board's own thermal source table. Q2/Q3 (buck leg) and
# Q4/Q5 (boost leg) all IRF3205 (RDS(on) ~8mOhm): conduction loss at the
# panel-side design current (5A) for the input-facing pair, a smaller share
# for the output-facing pair. RS1/RS2 from their own shunt values at the
# design currents (I^2R). L2 assumes a low-DCR (~30mOhm) power inductor at
# 5A. U8/U9 (integrated synchronous bucks) from (1-efficiency)*Pout at a
# vendor-typical ~88-90% efficiency for this current class. U1/U2 (LM74700
# ideal-diode controller, LM5165 bias buck) and U5/U6 (INA180) are small
# quiescent-current estimates.
THERMAL_SOURCES_W = {
    "Q2": 0.30,
    "Q3": 0.25,
    "Q4": 0.30,
    "Q5": 0.25,
    "U4": 0.40,
    "RS1": 0.25,
    "RS2": 0.13,
    "L2": 0.75,
    "L3": 0.32,
    "L4": 0.05,
    "U8": 0.60,
    "U9": 0.45,
    "U1": 0.05,
    "U2": 0.15,
    "U5": 0.02,
    "U6": 0.02,
}


def atom(x):
    return x.value() if isinstance(x, Symbol) else x


def children(lst, key):
    return [x for x in lst[1:] if isinstance(x, list) and x and atom(x[0]) == key]


def child(lst, key):
    xs = children(lst, key)
    return xs[0] if xs else None


def xy_points(node):
    pts = child(node, "pts") if node and atom(node[0]) != "pts" else node
    if not pts:
        return []
    return [(float(x[1]), float(x[2])) for x in pts[1:] if isinstance(x, list) and atom(x[0]) == "xy"]


def safe_polygon(points):
    if len(points) < 3:
        return None
    p = Polygon(points)
    if not p.is_valid:
        p = p.buffer(0)
    return p if not p.is_empty else None


def pad_shape(shape, sx, sy, angle_deg, x, y, round_ratio=0.25):
    if shape == "circle":
        g = Point(0, 0).buffer(sx / 2, resolution=20)
    elif shape == "oval":
        if sx >= sy:
            d = max(0.0, sx - sy)
            g = LineString([(-d / 2, 0), (d / 2, 0)]).buffer(sy / 2, resolution=16)
        else:
            d = max(0.0, sy - sx)
            g = LineString([(0, -d / 2), (0, d / 2)]).buffer(sx / 2, resolution=16)
    elif shape == "roundrect":
        r = max(0.01, min(sx, sy) * round_ratio)
        g = box(-sx / 2 + r, -sy / 2 + r, sx / 2 - r, sy / 2 - r).buffer(r, resolution=8)
    else:
        g = box(-sx / 2, -sy / 2, sx / 2, sy / 2)
    return translate(rotate(g, angle_deg, origin=(0, 0), use_radians=False), xoff=x, yoff=y)


@dataclass
class Pad:
    key: str
    reference: str
    number: str
    net: str
    center: tuple[float, float]
    geom: object
    layers: list[str]
    plated: bool
    drill_mm: float
    size_mm: float


class Board:
    def __init__(self, path: Path):
        self.path = path
        self.root = loads(path.read_text(encoding="utf-8"))
        self.layer_geom = defaultdict(lambda: defaultdict(list))
        self.pads = defaultdict(list)
        self.connectors = defaultdict(list)
        self.footprints = {}
        self.copper_thickness_m = {layer: 35e-6 for layer in COPPER_LAYERS}
        self.board_outline = None
        self._parse()

    def _parse(self):
        setup = child(self.root, "setup")
        stackup = child(setup, "stackup") if setup else None
        if stackup:
            for layer in children(stackup, "layer"):
                lname = str(layer[1])
                if lname in COPPER_LAYERS:
                    thick = child(layer, "thickness")
                    if thick:
                        self.copper_thickness_m[lname] = float(thick[1]) * 1e-3

        edge_lines = []
        for gl in children(self.root, "gr_line"):
            lay = child(gl, "layer")
            if lay and lay[1] == "Edge.Cuts":
                s, e = child(gl, "start"), child(gl, "end")
                edge_lines.append(LineString([(float(s[1]), float(s[2])), (float(e[1]), float(e[2]))]))
        polys = list(polygonize(edge_lines))
        self.board_outline = max(polys, key=lambda p: p.area) if polys else box(4, 4, 246, 156)

        for fp in children(self.root, "footprint"):
            at = child(fp, "at")
            fx, fy = float(at[1]), float(at[2])
            fa = float(at[3]) if len(at) > 3 else 0.0
            props = {str(p[1]): str(p[2]) for p in children(fp, "property") if len(p) > 2}
            ref, val = props.get("Reference", "?"), props.get("Value", "?")
            self.footprints[ref] = {"x": fx, "y": fy, "angle": fa, "value": val}
            for pd in children(fp, "pad"):
                number = str(pd[1])
                ptype, shape = atom(pd[2]), atom(pd[3])
                pat, size, netn, layers = child(pd, "at"), child(pd, "size"), child(pd, "net"), child(pd, "layers")
                if not size or not netn or not layers:
                    continue
                lx, ly = float(pat[1]), float(pat[2])
                pa = float(pat[3]) if len(pat) > 3 else 0.0
                theta = math.radians(fa)
                # KiCad places footprint pads CLOCKWISE in board coordinates,
                # because board Y points down. The counter-clockwise form
                # (lx*cos - ly*sin, lx*sin + ly*cos) that was here mirrors a
                # rotated footprint's pad row end-for-end: on U9, a SOT-563
                # placed at 90 degrees, it put pad 4 where pad 6 belongs, 1.4mm
                # away. Unrotated parts are unaffected, which is why this
                # survived six revisions -- it only bites when a rotated
                # footprint's pad is a terminal or carries a net's only copper,
                # and then it reports a perfectly continuous net as broken.
                px = fx + lx * math.cos(theta) + ly * math.sin(theta)
                py = fy + ly * math.cos(theta) - lx * math.sin(theta)
                sx, sy = float(size[1]), float(size[2])
                rr = child(pd, "roundrect_rratio")
                # A pad's stored angle is ABSOLUTE -- KiCad does not add the
                # footprint's rotation to it when plotting, so neither may we.
                # Adding fa here rotated every already-absolute pad twice.
                geom = pad_shape(shape, sx, sy, pa, px, py, float(rr[1]) if rr else 0.25)
                drill = child(pd, "drill")
                drill_mm = 0.0
                if drill and len(drill) > 1:
                    nums = [float(v) for v in drill[1:] if isinstance(v, (int, float))]
                    drill_mm = min(nums) if nums else 0.0
                if drill_mm > 0:
                    geom = geom.difference(Point(px, py).buffer(drill_mm / 2, resolution=16))
                ls = [str(atom(v)) for v in layers[1:]]
                cu_layers = COPPER_LAYERS[:] if "*.Cu" in ls else [z for z in ls if z in COPPER_LAYERS]
                net = str(netn[1])
                pad = Pad(f"{ref}.{number}", ref, number, net, (px, py), geom, cu_layers, ptype == "thru_hole", drill_mm, max(sx, sy))
                self.pads[pad.key].append(pad)
                for layer in cu_layers:
                    self.layer_geom[net][layer].append(geom)
                if pad.plated:
                    self.connectors[net].append({"center": (px, py), "drill": drill_mm, "size": max(sx, sy), "layers": cu_layers, "kind": "pth"})

        for seg in children(self.root, "segment"):
            s, e, width, layer, netn = child(seg, "start"), child(seg, "end"), child(seg, "width"), child(seg, "layer"), child(seg, "net")
            if not all([s, e, width, layer, netn]):
                continue
            g = LineString([(float(s[1]), float(s[2])), (float(e[1]), float(e[2]))]).buffer(float(width[1]) / 2, cap_style=1, join_style=1)
            self.layer_geom[str(netn[1])][str(layer[1])].append(g)

        for via in children(self.root, "via"):
            at, size, drill, layers, netn = child(via, "at"), child(via, "size"), child(via, "drill"), child(via, "layers"), child(via, "net")
            if not all([at, size, drill, layers, netn]):
                continue
            x, y, sz, dr = float(at[1]), float(at[2]), float(size[1]), float(drill[1])
            l0, l1 = str(layers[1]), str(layers[2])
            i0, i1 = COPPER_LAYERS.index(l0), COPPER_LAYERS.index(l1)
            ls = COPPER_LAYERS[min(i0, i1) : max(i0, i1) + 1]
            g = Point(x, y).buffer(sz / 2, resolution=16).difference(Point(x, y).buffer(dr / 2, resolution=16))
            net = str(netn[1])
            for layer in ls:
                self.layer_geom[net][layer].append(g)
            self.connectors[net].append({"center": (x, y), "drill": dr, "size": sz, "layers": ls, "kind": "via"})

        for zone in children(self.root, "zone"):
            netn = child(zone, "net")
            if not netn:
                continue
            net = str(netn[1])
            for fp in children(zone, "filled_polygon"):
                layer = child(fp, "layer")
                p = safe_polygon(xy_points(fp))
                if layer and p is not None:
                    self.layer_geom[net][str(layer[1])].append(p)

        for net, per_layer in self.layer_geom.items():
            for layer, geoms in list(per_layer.items()):
                good = [g for g in geoms if g is not None and not g.is_empty]
                per_layer[layer] = unary_union(good).buffer(0) if good else Polygon()

    @property
    def bounds(self):
        return self.board_outline.bounds


def rasterize(geom, bounds, pitch_mm):
    minx, miny, maxx, maxy = bounds
    nx, ny = math.ceil((maxx - minx) / pitch_mm), math.ceil((maxy - miny) / pitch_mm)
    xs = minx + (np.arange(nx) + 0.5) * pitch_mm
    ys = miny + (np.arange(ny) + 0.5) * pitch_mm
    xx, yy = np.meshgrid(xs, ys)
    if geom is None or geom.is_empty:
        return np.zeros((ny, nx), dtype=bool), xs, ys
    return contains_xy(geom.buffer(pitch_mm * 0.08), xx, yy), xs, ys


COVERAGE_SUB = 3          # sub-samples per cell edge, so 9 per cell


def rasterize_coverage(geom, bounds, pitch_mm, sub=COVERAGE_SUB):
    """Fraction of each cell covered by copper, not just "is the centre in it".

    A cell test that asks only about the centre cannot see a conductor
    narrower than the cell. This board's hand-router closes some links at
    0.15mm against a 0.20mm grid, so those traces lit up a broken chain of
    cells and the solver called five terminals disconnected -- on copper that
    KiCad's connectivity, the exact-geometry audit and a union-find over the
    raw polygons all agreed was one continuous piece.

    Sub-sampling fixes both halves. A cell touched by any copper becomes a
    node, so connectivity follows the copper; the coverage fraction then
    weights that node's conductance, so a 0.15mm trace in a 0.20mm cell
    conducts like 0.15mm of copper instead of like 0.20mm. The old
    all-or-nothing rule was not just coarse, it was biased: a conductor either
    counted as a full cell wide or did not exist at all.
    """
    minx, miny, maxx, maxy = bounds
    nx, ny = math.ceil((maxx - minx) / pitch_mm), math.ceil((maxy - miny) / pitch_mm)
    xs = minx + (np.arange(nx) + 0.5) * pitch_mm
    ys = miny + (np.arange(ny) + 0.5) * pitch_mm
    if geom is None or geom.is_empty:
        return np.zeros((ny, nx), dtype=float), xs, ys
    off = (np.arange(sub) + 0.5) / sub - 0.5
    cov = np.zeros((ny, nx), dtype=float)
    for dy in off:
        for dx in off:
            gx, gy = np.meshgrid(xs + dx * pitch_mm, ys + dy * pitch_mm)
            cov += contains_xy(geom, gx, gy)
    return cov / (sub * sub), xs, ys


def add_edge(rows, cols, vals, diag, a, b, conductance):
    if a < 0 or b < 0 or a == b or conductance <= 0:
        return
    rows.extend([a, b])
    cols.extend([b, a])
    vals.extend([-conductance, -conductance])
    diag[a] += conductance
    diag[b] += conductance


def terminal_cell_ids(board, net, pad_key, ids, xs, ys):
    out = []
    for pad in board.pads.get(pad_key, []):
        if pad.net != net:
            continue
        for li, layer in enumerate(COPPER_LAYERS):
            if layer not in pad.layers:
                continue
            minx, miny, maxx, maxy = pad.geom.bounds
            ix0, ix1 = max(0, np.searchsorted(xs, minx) - 1), min(len(xs), np.searchsorted(xs, maxx) + 1)
            iy0, iy1 = max(0, np.searchsorted(ys, miny) - 1), min(len(ys), np.searchsorted(ys, maxy) + 1)
            if ix1 <= ix0 or iy1 <= iy0:
                continue
            xx, yy = np.meshgrid(xs[ix0:ix1], ys[iy0:iy1])
            inside = contains_xy(pad.geom.buffer((xs[1] - xs[0]) * 0.15), xx, yy)
            candidates = ids[li, iy0:iy1, ix0:ix1][inside]
            out.extend(int(v) for v in candidates if v >= 0)
    return sorted(set(out))


def solve_current_case(board, case, pitch_mm=GRID_MM, plating_m=VIA_PLATING_M, save_field=True):
    net = case["net"]
    covs, xs, ys = [], None, None
    for layer in COPPER_LAYERS:
        cov, xs, ys = rasterize_coverage(board.layer_geom[net].get(layer),
                                         board.bounds, pitch_mm)
        covs.append(cov)
    covs = np.asarray(covs)
    masks = covs > 0.0
    ids = np.full(masks.shape, -1, dtype=np.int32)
    ids[masks] = np.arange(int(masks.sum()), dtype=np.int32)
    n = int(masks.sum())
    if n == 0:
        raise RuntimeError(f"No conductor geometry for net {net}")

    rows, cols, vals = [], [], []
    diag = np.zeros(n)

    def pair_g(base_g, ca, cb):
        """Conductance between two adjacent cells of fractional coverage.

        Two half-cells in series, each conducting in proportion to the copper
        actually in it: 2*g*ca in series with 2*g*cb, i.e. g times the
        harmonic mean 2*ca*cb/(ca+cb). Two full cells give exactly the old
        value, so thick copper is unchanged and only partly-filled cells move.
        """
        s = ca + cb
        out = np.zeros_like(s, dtype=float)
        nz = s > 0
        out[nz] = base_g * 2.0 * ca[nz] * cb[nz] / s[nz]
        return out

    for li, layer in enumerate(COPPER_LAYERS):
        sheet_g = board.copper_thickness_m[layer] / RHO_CU_20
        a = ids[li, :, :-1]
        b = ids[li, :, 1:]
        valid = (a >= 0) & (b >= 0)
        aa, bb = a[valid], b[valid]
        gg = pair_g(sheet_g, covs[li][:, :-1][valid], covs[li][:, 1:][valid])
        rows.extend(aa.tolist()); cols.extend(bb.tolist()); vals.extend((-gg).tolist())
        rows.extend(bb.tolist()); cols.extend(aa.tolist()); vals.extend((-gg).tolist())
        np.add.at(diag, aa, gg); np.add.at(diag, bb, gg)
        a = ids[li, :-1, :]
        b = ids[li, 1:, :]
        valid = (a >= 0) & (b >= 0)
        aa, bb = a[valid], b[valid]
        gg = pair_g(sheet_g, covs[li][:-1, :][valid], covs[li][1:, :][valid])
        rows.extend(aa.tolist()); cols.extend(bb.tolist()); vals.extend((-gg).tolist())
        rows.extend(bb.tolist()); cols.extend(aa.tolist()); vals.extend((-gg).tolist())
        np.add.at(diag, aa, gg); np.add.at(diag, bb, gg)

        # Diagonal neighbors, at conductance/sqrt(2) (the correct scaling for
        # a resistive-sheet finite-difference stencil, since the diagonal
        # path length is sqrt(2)x the orthogonal cell pitch). Needed because
        # this board's tracks were placed by a custom router that (unlike a
        # typical 45-degree-snapped autorouter run) uses genuinely diagonal
        # segments freely -- at this rasterization pitch, two such segments
        # meeting at a vertex can produce grid cells that only share a
        # corner, not an edge, which a 4-connected-only stencil reports as
        # electrically disconnected even though the real copper is one
        # continuous polygon. Confirmed empirically: without this, a board
        # with 0 DRC unconnected_items still showed 17 disconnected raster
        # components on a single net.
        diag_g = sheet_g / math.sqrt(2)
        a = ids[li, :-1, :-1]
        b = ids[li, 1:, 1:]
        valid = (a >= 0) & (b >= 0)
        aa, bb = a[valid], b[valid]
        gg = pair_g(diag_g, covs[li][:-1, :-1][valid], covs[li][1:, 1:][valid])
        rows.extend(aa.tolist()); cols.extend(bb.tolist()); vals.extend((-gg).tolist())
        rows.extend(bb.tolist()); cols.extend(aa.tolist()); vals.extend((-gg).tolist())
        np.add.at(diag, aa, gg); np.add.at(diag, bb, gg)
        a = ids[li, :-1, 1:]
        b = ids[li, 1:, :-1]
        valid = (a >= 0) & (b >= 0)
        aa, bb = a[valid], b[valid]
        gg = pair_g(diag_g, covs[li][:-1, 1:][valid], covs[li][1:, :-1][valid])
        rows.extend(aa.tolist()); cols.extend(bb.tolist()); vals.extend((-gg).tolist())
        rows.extend(bb.tolist()); cols.extend(aa.tolist()); vals.extend((-gg).tolist())
        np.add.at(diag, aa, gg); np.add.at(diag, bb, gg)

    connector_count = 0
    for con in board.connectors.get(net, []):
        x, y = con["center"]
        layer_indices = sorted(COPPER_LAYERS.index(z) for z in con["layers"] if z in COPPER_LAYERS)
        if len(layer_indices) < 2 or con["drill"] <= 0:
            continue
        groups = {}
        # The raster can miss a narrow thermal-relief spoke whose center falls
        # between cells. Search past the pad edge to the same-net filled copper,
        # then distribute barrel conductance across the local annulus/pour cells.
        # Because groups are selected from this net's mask only, this cannot bridge
        # to a different electrical net across a clearance.
        radius = max(con["size"] / 2 + max(0.8, pitch_mm), pitch_mm)
        for li in layer_indices:
            ix0, ix1 = max(0, np.searchsorted(xs, x - radius) - 1), min(len(xs), np.searchsorted(xs, x + radius) + 1)
            iy0, iy1 = max(0, np.searchsorted(ys, y - radius) - 1), min(len(ys), np.searchsorted(ys, y + radius) + 1)
            yy, xx = np.meshgrid(ys[iy0:iy1], xs[ix0:ix1], indexing="ij")
            near = ((xx - x) ** 2 + (yy - y) ** 2) <= radius**2
            groups[li] = sorted(set(int(v) for v in ids[li, iy0:iy1, ix0:ix1][near] if v >= 0))
        for la, lb in zip(layer_indices[:-1], layer_indices[1:]):
            ga, gb = groups.get(la, []), groups.get(lb, [])
            if not ga or not gb:
                continue
            z_m = sum(board.copper_thickness_m[COPPER_LAYERS[k]] for k in range(la, lb + 1))
            # Physical separation follows the 1.6 mm stack roughly; distribute by layer span.
            z_m = max(z_m, 1.6e-3 * (lb - la) / (len(COPPER_LAYERS) - 1))
            barrel_area = math.pi * con["drill"] * 1e-3 * plating_m
            gtotal = barrel_area / (RHO_CU_20 * z_m)
            gpair = gtotal / (len(ga) * len(gb))
            for a in ga:
                for b in gb:
                    add_edge(rows, cols, vals, diag, a, b, gpair)
            connector_count += 1

    rows.extend(range(n)); cols.extend(range(n)); vals.extend(diag.tolist())
    G = sparse.csr_matrix((vals, (rows, cols)), shape=(n, n))
    injections = np.zeros(n)
    terminal_ids = {}
    for key, amps in case["terminals"].items():
        tids = terminal_cell_ids(board, net, key, ids, xs, ys)
        if not tids:
            raise RuntimeError(f"{case['name']}: terminal {key} not found on {net}")
        terminal_ids[key] = tids
        injections[tids] += amps / len(tids)
    if abs(float(injections.sum())) > 1e-9:
        raise RuntimeError(f"Unbalanced current injection in {case['name']}")

    active = np.flatnonzero(diag > 0)
    # Keep only the graph component containing the first terminal. Floating copper
    # remains visible in the board but must not carry invented load current.
    sub = G[active][:, active]
    comps, labels = sparse.csgraph.connected_components(sub, directed=False)
    first_node = next(iter(terminal_ids.values()))[0]
    first_pos = int(np.searchsorted(active, first_node))
    target_label = labels[first_pos]
    keep = active[labels == target_label]
    keep_set = set(int(v) for v in keep)
    for key, tids in terminal_ids.items():
        if not any(t in keep_set for t in tids):
            raise RuntimeError(f"{case['name']}: terminal {key} is disconnected")

    ground = int(next(t for key, amps in case["terminals"].items() if amps < 0 for t in terminal_ids[key] if t in keep_set))
    solve_nodes = keep[keep != ground]
    A = G[solve_nodes][:, solve_nodes]
    bvec = injections[solve_nodes]
    pre = sparse.diags(1.0 / np.maximum(A.diagonal(), 1e-30))
    vsolve, info = spla.cg(A, bvec, rtol=1e-10, atol=1e-13, maxiter=5000, M=pre)
    if info != 0:
        vsolve = spla.spsolve(A, bvec)
    voltage = np.zeros(n)
    voltage[solve_nodes] = vsolve

    total_positive = sum(v for v in case["terminals"].values() if v > 0)
    source_v_num, sink_v_num = 0.0, 0.0
    for key, amps in case["terminals"].items():
        tids = [t for t in terminal_ids[key] if t in keep_set]
        av = float(np.mean(voltage[tids]))
        if amps > 0:
            source_v_num += av * amps
        else:
            sink_v_num += av * (-amps)
    vdrop = (source_v_num - sink_v_num) / total_positive
    power_iv = float(np.dot(injections, voltage))
    power_edge = float(voltage @ (G @ voltage))
    residual = float(np.linalg.norm((G @ voltage - injections)[keep[keep != ground]]) / max(total_positive, 1e-12))

    vgrid = np.full(ids.shape, np.nan)
    for li in range(len(COPPER_LAYERS)):
        valid = ids[li] >= 0
        vgrid[li][valid] = voltage[ids[li][valid]]
    jgrid = np.full(ids.shape, np.nan)
    pitch_m = pitch_mm * 1e-3
    for li, layer in enumerate(COPPER_LAYERS):
        vv = vgrid[li]
        gx = np.zeros_like(vv); gy = np.zeros_like(vv)
        gx[:, 1:-1] = (vv[:, 2:] - vv[:, :-2]) / (2 * pitch_m)
        gy[1:-1, :] = (vv[2:, :] - vv[:-2, :]) / (2 * pitch_m)
        j = np.sqrt(gx * gx + gy * gy) / RHO_CU_20 / 1e6  # A/mm^2
        j[~masks[li]] = np.nan
        jgrid[li] = j
    finite_j = jgrid[np.isfinite(jgrid)]
    if save_field:
        field = np.max(np.where(np.isfinite(jgrid), jgrid, -np.inf), axis=0)
        field[field == -np.inf] = np.nan
    else:
        field = None
    result = {
        "name": case["name"], "net": net, "pitch_mm": pitch_mm,
        "load_A": total_positive, "voltage_drop_mV": vdrop * 1e3,
        "resistance_mOhm": vdrop / total_positive * 1e3,
        "copper_loss_mW": power_iv * 1e3,
        "joule_balance_error_pct": abs(power_edge - power_iv) / max(abs(power_iv), 1e-15) * 100,
        "normalized_residual": residual,
        "j_p95_A_per_mm2": float(np.nanpercentile(finite_j, 95)),
        "j_p99_A_per_mm2": float(np.nanpercentile(finite_j, 99)),
        "j_max_A_per_mm2": float(np.nanmax(finite_j)),
        "conductor_nodes": int(len(keep)), "floating_components_omitted": int(max(0, comps - 1)),
        "plated_layer_links": connector_count,
        "resistance_at_80C_mOhm": vdrop / total_positive * 1e3 * (1 + ALPHA_CU * 60),
    }
    return result, field, (xs, ys), masks


def solve_thermal(board, h_conv=10.0, pitch_mm=THERMAL_GRID_MM):
    board_mask, xs, ys = rasterize(board.board_outline, board.bounds, pitch_mm)
    ny, nx = board_mask.shape
    copper_count = np.zeros((ny, nx), dtype=float)
    all_net_layer = {layer: [] for layer in COPPER_LAYERS}
    for per_layer in board.layer_geom.values():
        for layer in COPPER_LAYERS:
            g = per_layer.get(layer)
            if g is not None and not g.is_empty:
                all_net_layer[layer].append(g)
    for layer in COPPER_LAYERS:
        g = unary_union(all_net_layer[layer]).buffer(0) if all_net_layer[layer] else Polygon()
        m, _, _ = rasterize(g, board.bounds, pitch_mm)
        copper_count += m.astype(float)

    ids = np.full((ny, nx), -1, dtype=np.int32)
    ids[board_mask] = np.arange(int(board_mask.sum()), dtype=np.int32)
    n = int(board_mask.sum())
    rows, cols, vals, diag = [], [], [], np.zeros(n)
    k_fr4, k_cu = 0.30, 385.0
    total_cu_m = sum(board.copper_thickness_m.values())
    fr4_m = max(0.0, 1.6e-3 - total_cu_m)
    local_sheet_k = k_fr4 * fr4_m + k_cu * 35e-6 * copper_count

    for agrid, bgrid, ka, kb in [
        (ids[:, :-1], ids[:, 1:], local_sheet_k[:, :-1], local_sheet_k[:, 1:]),
        (ids[:-1, :], ids[1:, :], local_sheet_k[:-1, :], local_sheet_k[1:, :]),
    ]:
        valid = (agrid >= 0) & (bgrid >= 0)
        aa, bb = agrid[valid], bgrid[valid]
        kval = 2 * ka[valid] * kb[valid] / np.maximum(ka[valid] + kb[valid], 1e-30)
        rows.extend(aa.tolist()); cols.extend(bb.tolist()); vals.extend((-kval).tolist())
        rows.extend(bb.tolist()); cols.extend(aa.tolist()); vals.extend((-kval).tolist())
        np.add.at(diag, aa, kval); np.add.at(diag, bb, kval)

    sigma_sb, emissivity, trad_k = 5.670374419e-8, 0.90, 313.15
    h_rad = 4 * emissivity * sigma_sb * trad_k**3
    cell_area_m2 = (pitch_mm * 1e-3) ** 2
    sink_g = 2 * (h_conv + h_rad) * cell_area_m2
    diag += sink_g
    rows.extend(range(n)); cols.extend(range(n)); vals.extend(diag.tolist())
    A = sparse.csr_matrix((vals, (rows, cols)), shape=(n, n))
    q = np.zeros(n)
    source_cells = {}
    source_radii = {
        "L2": 4.0, "L3": 3.0, "L4": 2.5,
        "Q2": 4.5, "Q3": 4.5, "Q4": 4.5, "Q5": 4.5,   # TO-220, larger body than default SMD
        "U4": 4.0,   # HTSSOP-28, 9.7x4.4mm
        "U8": 2.0,   # QFN-16 3x3mm
        "U9": 1.5,   # SOT-563, ~1.6x1.6mm
        "RS1": 2.5, "RS2": 2.5,
        "U1": 2.0, "U2": 2.0, "U5": 1.5, "U6": 1.5,
    }
    xx, yy = np.meshgrid(xs, ys)
    for ref, watts in THERMAL_SOURCES_W.items():
        fp = board.footprints[ref]
        radius = source_radii.get(ref, 2.5)
        sel = board_mask & (((xx - fp["x"]) ** 2 + (yy - fp["y"]) ** 2) <= radius**2)
        nodes = ids[sel]
        if len(nodes) == 0:
            iy, ix = np.unravel_index(np.argmin((xx - fp["x"]) ** 2 + (yy - fp["y"]) ** 2), xx.shape)
            nodes = np.array([ids[iy, ix]])
        q[nodes] += watts / len(nodes)
        source_cells[ref] = nodes
    pre = sparse.diags(1.0 / np.maximum(A.diagonal(), 1e-30))
    rise, info = spla.cg(A, q, rtol=1e-10, atol=1e-12, maxiter=5000, M=pre)
    if info != 0:
        rise = spla.spsolve(A, q)
    tgrid = np.full((ny, nx), np.nan)
    tgrid[board_mask] = 25.0 + rise
    source_t = {ref: float(np.mean(25.0 + rise[nodes])) for ref, nodes in source_cells.items()}
    imax = int(np.argmax(rise))
    iy, ix = np.argwhere(ids == imax)[0]
    result = {
        "h_convection_W_m2K_per_side": h_conv,
        "h_radiation_W_m2K_per_side": h_rad,
        "ambient_C": 25.0,
        "total_power_W": float(q.sum()),
        "max_board_C": float(25.0 + rise.max()),
        "average_board_C": float(25.0 + rise.mean()),
        "hotspot_x_mm": float(xs[ix]), "hotspot_y_mm": float(ys[iy]),
        "energy_balance_error_pct": abs(float(q.sum()) - float(sink_g * rise.sum())) / float(q.sum()) * 100,
        "source_average_board_C": source_t,
        "copper_layer_coverage_pct": {layer: float(100 * np.mean(rasterize(unary_union(all_net_layer[layer]).buffer(0), board.bounds, pitch_mm)[0][board_mask])) for layer in COPPER_LAYERS},
    }
    return result, tgrid, copper_count, (xs, ys)


def plot_current_maps(fields, results, xys, out):
    n = len(results)
    cols = 3
    rows = math.ceil(n / cols)
    fig, axes = plt.subplots(rows, cols, figsize=(15, 4.2 * rows), constrained_layout=True)
    axes = np.atleast_1d(axes).ravel()
    for ax, res in zip(axes, results):
        arr = fields[res["name"]]
        xs, ys = xys[res["name"]]
        positive = arr[np.isfinite(arr) & (arr > 0)]
        vmax = max(float(np.nanpercentile(positive, 99.5)) if len(positive) else 1e-6, 1e-6)
        vmin = max(vmax / 1000, 1e-5)
        im = ax.imshow(arr, extent=(xs[0], xs[-1], ys[-1], ys[0]), cmap="inferno", norm=LogNorm(vmin=vmin, vmax=vmax), aspect="equal")
        ax.set_title(f"{res['name']}\n{res['resistance_mOhm']:.3f} mΩ, {res['copper_loss_mW']:.2f} mW")
        ax.set_xlabel("x (mm)"); ax.set_ylabel("y (mm)")
        fig.colorbar(im, ax=ax, label="current density (A/mm², max over layers)", shrink=0.78)
    for ax in axes[n:]:
        ax.axis("off")
    fig.suptitle("MPPT board — geometry-derived DC current-density screening", fontsize=16)
    fig.savefig(out, dpi=180)
    plt.close(fig)


def plot_thermal(tgrid, copper_count, xy, thermal, out):
    xs, ys = xy
    fig, axes = plt.subplots(1, 2, figsize=(15, 6), constrained_layout=True)
    im = axes[0].imshow(tgrid, extent=(xs[0], xs[-1], ys[-1], ys[0]), cmap="turbo", aspect="equal")
    axes[0].scatter([thermal["hotspot_x_mm"]], [thermal["hotspot_y_mm"]], marker="x", s=90, c="black")
    axes[0].set_title(f"Board-level steady-state temperature\nmax {thermal['max_board_C']:.1f}°C, ambient 25°C")
    axes[0].set_xlabel("x (mm)"); axes[0].set_ylabel("y (mm)")
    fig.colorbar(im, ax=axes[0], label="board temperature (°C)")
    cm = axes[1].imshow(copper_count, extent=(xs[0], xs[-1], ys[-1], ys[0]), cmap="viridis", vmin=0, vmax=4, aspect="equal")
    axes[1].set_title("Rasterized copper presence used by thermal model")
    axes[1].set_xlabel("x (mm)"); axes[1].set_ylabel("y (mm)")
    fig.colorbar(cm, ax=axes[1], label="copper layers present")
    fig.suptitle("MPPT board — collapsed 2D electrothermal screening model", fontsize=16)
    fig.savefig(out, dpi=200)
    plt.close(fig)


def write_csv(path, rows):
    if not rows:
        return
    keys = list(rows[0].keys())
    with path.open("w", newline="", encoding="utf-8") as f:
        w = csv.DictWriter(f, fieldnames=keys)
        w.writeheader(); w.writerows(rows)


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("board", type=Path)
    ap.add_argument("output", type=Path)
    args = ap.parse_args()
    args.output.mkdir(parents=True, exist_ok=True)
    board = Board(args.board)
    board_hash = hashlib.sha256(args.board.read_bytes()).hexdigest()

    current_results, fields, xys = [], {}, {}
    masks_for_npz = {}
    unsolved = []
    for case in CURRENT_CASES:
        pitch = GND_GRID_MM if case["net"] == "GND" else GRID_MM
        print(f"current: {case['name']} at {pitch:.2f} mm", flush=True)
        try:
            res, field, xy, masks = solve_current_case(board, case, pitch_mm=pitch)
        except RuntimeError as exc:
            # A case whose terminal has no copper is not a solver failure, it
            # is a finding: that net is not routed. Record it and carry on, so
            # one unrouted rail does not hide the results for the other
            # thirteen cases. The unsolved list is reported, not suppressed.
            print(f"   SKIPPED -- {exc}", flush=True)
            unsolved.append({"name": case["name"], "net": case["net"],
                             "reason": str(exc)})
            continue
        current_results.append(res); fields[res["name"]] = field.astype(np.float32)
        xys[res["name"]] = xy
        masks_for_npz[res["name"] + "__mask"] = np.any(masks, axis=0)
    if unsolved:
        print(f"\n{len(unsolved)} case(s) could not be solved "
              f"(terminal not reached by copper):", flush=True)
        for u in unsolved:
            print(f"   {u['name']}  [{u['net']}]  {u['reason']}", flush=True)

    fine_cases = {"panel_raw_5A", "pol5v_switch_node_2A", "load_ground_return_3A"}
    convergence = []
    for case in CURRENT_CASES:
        if case["name"] not in fine_cases:
            continue
        coarse = next((r for r in current_results if r["name"] == case["name"]), None)
        if coarse is None:
            continue
        fine_pitch = GND_FINE_GRID_MM if case["net"] == "GND" else FINE_GRID_MM
        print(f"convergence: {case['name']} at {fine_pitch:.2f} mm", flush=True)
        fine, _, _, _ = solve_current_case(board, case, pitch_mm=fine_pitch, save_field=False)
        convergence.append({
            "name": case["name"], "coarse_pitch_mm": coarse["pitch_mm"], "fine_pitch_mm": fine_pitch,
            "coarse_resistance_mOhm": coarse["resistance_mOhm"], "fine_resistance_mOhm": fine["resistance_mOhm"],
            "difference_pct": 100 * abs(fine["resistance_mOhm"] - coarse["resistance_mOhm"]) / max(fine["resistance_mOhm"], 1e-12),
        })

    # Via-plating sensitivity on the large, stitched GND return network.
    plating_sensitivity = []
    gcase = next(c for c in CURRENT_CASES if c["name"] == "load_ground_return_3A")
    for plating_um in [15, 25, 35]:
        print(f"via plating: {plating_um} um", flush=True)
        if plating_um == 25:
            rr = next(r for r in current_results if r["name"] == "load_ground_return_3A")
        else:
            rr, _, _, _ = solve_current_case(board, gcase, pitch_mm=GND_GRID_MM, plating_m=plating_um * 1e-6, save_field=False)
        plating_sensitivity.append({"plating_um": plating_um, "resistance_mOhm": rr["resistance_mOhm"], "copper_loss_mW": rr["copper_loss_mW"]})

    thermal_runs, base_tgrid, copper_count, thermal_xy = [], None, None, None
    for h in [5.0, 10.0, 20.0]:
        print(f"thermal: h={h:.0f} W/m2K per side", flush=True)
        tr, tg, cc, txy = solve_thermal(board, h_conv=h)
        thermal_runs.append(tr)
        if h == 10.0:
            base_tgrid, copper_count, thermal_xy = tg, cc, txy
    base_thermal = next(t for t in thermal_runs if t["h_convection_W_m2K_per_side"] == 10.0)

    # Rectangular-conductor discretization sanity check using the same square-cell sheet model.
    length_mm, width_mm = 50.0, 2.0
    analytic = RHO_CU_20 * (length_mm * 1e-3) / (35e-6 * width_mm * 1e-3)
    nx, ny = round(length_mm / GRID_MM), round(width_mm / GRID_MM)
    discrete = (RHO_CU_20 / 35e-6) * (nx - 1) / ny
    validation = {"rectangle_length_mm": length_mm, "rectangle_width_mm": width_mm, "analytic_mOhm": analytic * 1e3, "discrete_mOhm": discrete * 1e3, "error_pct": abs(discrete - analytic) / analytic * 100}

    write_csv(args.output / "current-density-summary.csv", current_results)
    write_csv(args.output / "mesh-convergence.csv", convergence)
    write_csv(args.output / "via-plating-sensitivity.csv", plating_sensitivity)
    thermal_csv = []
    for t in thermal_runs:
        thermal_csv.append({k: v for k, v in t.items() if not isinstance(v, dict)})
    write_csv(args.output / "thermal-sensitivity.csv", thermal_csv)
    plot_current_maps(fields, current_results, xys, args.output / "current-density-maps.png")
    plot_thermal(base_tgrid, copper_count, thermal_xy, base_thermal, args.output / "thermal-map.png")
    np.savez_compressed(args.output / "current-fields.npz", **fields, **masks_for_npz)
    np.savez_compressed(args.output / "thermal-field.npz", temperature_C=base_tgrid, copper_layer_count=copper_count, xs_mm=thermal_xy[0], ys_mm=thermal_xy[1])

    report_data = {
        "board": str(args.board), "board_sha256": board_hash,
        "unsolved_cases": unsolved,
        "model": {"current_grid_mm": GRID_MM, "fine_grid_mm": FINE_GRID_MM, "ground_grid_mm": GND_GRID_MM, "ground_fine_grid_mm": GND_FINE_GRID_MM, "thermal_grid_mm": THERMAL_GRID_MM, "copper_resistivity_20C_ohm_m": RHO_CU_20, "copper_tempco_per_C": ALPHA_CU, "via_plating_assumption_um": VIA_PLATING_M * 1e6, "ambient_C": 25.0},
        "native_drc": {"kicad_version": "10.0.2", "unconnected": 0, "other_violations": 11, "note": "Run via kicad-cli with zone refill on 2026-09-19; unconnected=0 and 0 different-net clearance findings -- the board is electrically complete and clearance-clean. The 11 remaining violations are cosmetic only (4 footprint-library metadata mismatches, 7 silkscreen overlaps). See 04-pcb-layout-review.md."},
        "current_cases": current_results, "mesh_convergence": convergence,
        "via_plating_sensitivity": plating_sensitivity, "thermal_sources_W": THERMAL_SOURCES_W,
        "thermal_sensitivity": thermal_runs, "rectangle_validation": validation,
    }
    (args.output / "simulation-results.json").write_text(json.dumps(report_data, indent=2), encoding="utf-8")

    worst_drop = max(current_results, key=lambda r: r["voltage_drop_mV"])
    worst_r = max(current_results, key=lambda r: r["resistance_mOhm"])
    worst_loss = max(current_results, key=lambda r: r["copper_loss_mW"])
    max_conv = max(c["difference_pct"] for c in convergence)
    source_rows = "\n".join(f"| {ref} | {w:.3f} | {base_thermal['source_average_board_C'][ref]:.1f} |" for ref, w in THERMAL_SOURCES_W.items())
    cur_rows = "\n".join(
        f"| {r['name']} | `{r['net']}` | {r['load_A']:.2f} | {r['resistance_mOhm']:.3f} | {r['voltage_drop_mV']:.3f} | {r['copper_loss_mW']:.2f} | {r['j_p99_A_per_mm2']:.2f} |"
        for r in current_results
    )
    therm_rows = "\n".join(
        f"| {t['h_convection_W_m2K_per_side']:.0f} | {t['max_board_C']:.1f} | {t['average_board_C']:.1f} | {t['hotspot_x_mm']:.1f}, {t['hotspot_y_mm']:.1f} |"
        for t in thermal_runs
    )
    md = f"""# MPPT MCU/IC board electrothermal simulation

Run: 2026-09-19 (revision 4)
Board SHA-256: `{board_hash}`
KiCad native validation (KiCad 10.0.2, after zone refill): **0 unconnected items, 0 different-net electrical defects and 0 mechanical defects**. Every net on the board is one continuous piece of copper, so every case below solves -- including `pol3v3_switch_node_1A`, which earlier revisions had to record as unsolved because the `Net-(U9-SW)` escape was still open. The 11 remaining findings are cosmetic (7 silkscreen overlaps, 4 footprint-library metadata mismatches) and carry no electrical meaning. A separate first-build audit -- net parity, exact-geometry net continuity, fab capability, bypass-capacitor distance, exposed-pad stitching, pads overhanging the outline -- reports zero problems. See `04-pcb-layout-review.md`.

## Outcome

The saved board is electrically connected under KiCad's checked rules. In the geometry-derived DC model, the largest copper voltage drop is **{worst_drop['voltage_drop_mV']:.3f} mV** ({worst_drop['name']}), the highest effective copper resistance is **{worst_r['resistance_mOhm']:.3f} mΩ** ({worst_r['name']}), and the largest copper loss is **{worst_loss['copper_loss_mW']:.2f} mW** ({worst_loss['name']}). These figures exclude component, connector, fuse, shunt, inductor, contact, and solder-joint resistance.

The conservative composite thermal case injects **{sum(THERMAL_SOURCES_W.values()):.3f} W** into the board. At 25°C ambient with 10 W/m²K natural convection per side plus linearized solder-mask radiation, the collapsed board model predicts a **{base_thermal['max_board_C']:.1f}°C** peak board temperature and **{base_thermal['average_board_C']:.1f}°C** area average. This is a board-surface screening result, **not a semiconductor junction-temperature prediction**.

## What this model changed about the board

Three copper decisions on this board were made because of what this model
measured, not from a rule of thumb:

**Trace width comes from each net's real current, not its name.** An early
pass classified "power" nets by matching keywords against net names, and KiCad
auto-generates names like `Net-(Q4-D)` for any net never explicitly labelled
in the schematic. Several genuinely high-current nets matched no keyword and
were routed at signal width purely because of how they happened to be named;
the model measured one of them carrying its full 2.5A on a 0.25mm trace at
431mV / 172.5mOhm / 1078mW. Width is now derived from the schematic's own
topology -- the same per-net current list this model uses for its cases.

**Width is a per-segment property.** The router caps a link's width at the
narrowest pad it lands on. That is correct at the pad and wrong everywhere
else: a 5A net ran 0.4mm along its whole length because it touches one small
pad, which is about 357 A/mm^2 at 35um copper. Each segment is now grown as
far as the surrounding copper actually allows, verified against exact
geometry, so the 5A nets reach full policy width in open board and neck down
only at the pads.

**The charge path has a direct parallel run.** A net is routed as a chain, so
the path from the charge-current shunt RS2 to the battery connector J2
physically visited every other `/VBATT` pad on the way. The model measured
that as a routing-LENGTH problem rather than a width problem -- the trace was
already at full policy width -- so a single additional direct path was added
between those two pads. Both ends land on the same copper net, so it runs
electrically in parallel with the long one and parallel resistance is
dominated by the smaller term. Nothing existing was moved.

**What that bought, and what it cost.** The same model, run on the board at each
revision:

| Measurement | Rev 1 | Rev 2 | Rev 3 | Rev 4 (current) |
|---|---:|---:|---:|---:|
| Panel input path resistance | 11.578 mOhm | 5.454 mOhm | 3.866 mOhm | 3.866 mOhm |
| Panel input current density (p99) | 122.55 A/mm^2 | 55.24 A/mm^2 | 40.50 A/mm^2 | 40.50 A/mm^2 |
| Buck-side switch node copper loss | 333.13 mW | 160.35 mW | 112.84 mW | 164.74 mW |
| Raw VBATT resistance | 39.370 mOhm | 30.167 mOhm | 21.155 mOhm | 31.094 mOhm |
| Battery bus -> battery resistance | 29.591 mOhm | 20.180 mOhm | 11.121 mOhm | 8.645 mOhm |
| Battery bus -> 5V POL resistance | 76.856 mOhm | 30.297 mOhm | 21.143 mOhm | 60.186 mOhm |
| 5V rail resistance | 25.502 mOhm | 19.973 mOhm | 8.755 mOhm | 5.919 mOhm |
| 5V output resistance | 23.939 mOhm | 19.392 mOhm | 5.309 mOhm | 5.810 mOhm |
| 3.3V rail resistance | 18.606 mOhm | 25.300 mOhm | 7.638 mOhm | 6.991 mOhm |
| 3.3V output resistance | 24.754 mOhm | 31.069 mOhm | 10.913 mOhm | 10.913 mOhm |
| Worst copper loss anywhere | 333.13 mW | 188.54 mW | 132.22 mW | 194.34 mW |
| Peak current density, any case | 122.55 A/mm^2 | 61.29 A/mm^2 | 58.92 A/mm^2 | 67.03 A/mm^2 |
| Peak board temperature | 48.5 C | 45.7 C | 45.6 C | 46.1 C |

Revision 4 goes both ways, and that is the honest summary of it: four cases improved,
six got worse, and the three headline copper figures are among the worse ones. That
revision was aimed at first-build survival rather than at these numbers -- moving every
bypass capacitor next to its pin (the worst was 62.5mm away) and both POL inductors next
to their switch node re-shuffled the placement, and re-routing around the new positions
gave back some of the copper revision 3 had won. The trade is deliberate: a 54mV drop on
a 0.9A rail is a performance number with headroom, while a 62mm bypass capacitor is the
kind of defect that produces a board which does not work and cannot be debugged from the
schematic.

Two rows repay a careful read. The 3.3V rail and output went BACKWARDS in revision 2
(18.6 -> 25.3 and 24.8 -> 31.1 mOhm): that revision re-placed and re-routed from scratch,
discarding a parallel path the previous one had. The regression was visible here, and it
is what prompted revision 3 to give every resistance-dominated case its own parallel
path. And peak current density has nearly stopped responding to trace width
(122.55 -> 61.29 -> 58.92 -> 67.03): what remains sits in short runs at pads, where the
pad's own width sets the limit, so it is a placement and package question rather than a
trace-width one.

Board temperature barely moves across revisions, which is the expected result: the
composite case injects 4.29W of COMPONENT loss against roughly 0.4W of copper loss over
every rail combined. Copper work shows up in drop and current density, not in board
temperature.

This model also refuses to solve a case whose terminal has no copper reaching
it, and reports it as unsolved rather than aborting the run. That is how
unrouted rails surface here before they surface in DRC.

## DC current cases

| Case | Net | Load (A) | R (mΩ) | Drop (mV) | Copper loss (mW) | J p99 (A/mm²) |
|---|---|---:|---:|---:|---:|---:|
{cur_rows}

`j_max` is retained in the CSV/JSON but is contact- and mesh-sensitive; p99 is the more stable screening metric. All current plots show the maximum magnitude over the four copper layers. Floating copper components were omitted from the active solve rather than assigned invented current.

## Thermal sources and base-case board temperature

| Ref | Applied loss (W) | Average modeled board temperature near source (°C) |
|---|---:|---:|
{source_rows}

Q2-Q5 (IRF3205, RDS(on) ~8mOhm) use I^2R conduction loss at this design's worst-case panel-side current (5A) for the input-facing pair, a smaller share for the output-facing pair. RS1/RS2 use their own shunt values (0.01/0.02 Ohm) at the design currents. L2 assumes a low-DCR (~30mOhm) power inductor at 5A. U8/U9 (TPS62130, TLV62569 -- integrated synchronous bucks) use (1-efficiency)*Pout at a vendor-typical ~88-90% efficiency for this current class. U1/U2/U5/U6 are small quiescent-current estimates for the ideal-diode controller, bias buck, and current-sense amplifiers. None of this is bench-measured; applying every worst-case loss simultaneously is conservative and deliberately not a single rigorously coherent operating point.

| Convection per side (W/m²K) | Peak board (°C) | Average board (°C) | Hotspot x,y (mm) |
|---:|---:|---:|---:|
{therm_rows}

## Validation and sensitivity

- Current grid: {GRID_MM:.2f} mm for power rails and {GND_GRID_MM:.2f} mm for the broad GND planes; selected cases repeated at {FINE_GRID_MM:.2f}/{GND_FINE_GRID_MM:.2f} mm. Largest selected resistance change: **{max_conv:.2f}%**.
- Known 50 mm × 2 mm × 35 µm copper strip: analytic **{validation['analytic_mOhm']:.3f} mΩ**, discretized **{validation['discrete_mOhm']:.3f} mΩ** ({validation['error_pct']:.2f}% error).
- Copper modeled at 20°C. The CSV includes the uniform 80°C resistivity sensitivity (+{ALPHA_CU*60*100:.1f}%).
- Via plating is provisional at 25 µm; 15/25/35 µm sensitivity is in `via-plating-sensitivity.csv`.
- Joule power balance and solver residuals are recorded per case. Thermal energy-balance error is recorded per convection case.

## Important limits

- Native KiCad copper, pads, vias, tracks, and refilled zones are used. Rasterization still approximates curved edges, thermal spokes, and contact current crowding.
- The DC solve is not a switching/transient, EMI, fuse, connector, or component-resistance model. It does not establish a current rating.
- The thermal solve collapses the four copper layers and FR-4 into a 2D sheet. It assumes all listed component loss enters the PCB, both board faces see the same ambient, no enclosure/contact heatsink exists, and solder-mask emissivity is 0.9. Package-to-air paths, junction-to-case resistance, lead conduction, airflow obstruction, cables, and thermal feedback are omitted.
- The 5A worst-case current figure combines the 6-24V range's floor voltage with the 30W panel budget (30W / 6V). It is a conservative sizing convention, not a claim that 5A flows at every operating voltage.
- The board these results describe has 0 unconnected items, 0 different-net electrical and 0 mechanical DRC findings, and 11 cosmetic ones. Nothing was suppressed or edited for this run; see `04-pcb-layout-review.md` for the full breakdown.

## Artifacts

- `current-density-maps.png` — all DC cases
- `thermal-map.png` — base thermal case and copper-layer coverage
- `current-density-summary.csv`, `thermal-sensitivity.csv` — compact results
- `current-fields.npz`, `thermal-field.npz` — numerical fields
- `simulation-results.json` — full assumptions, results, checks, and sensitivities
- `run_simulations.py`, `requirements.txt` — reproducibility
"""
    (args.output / "05-electrothermal-simulation.md").write_text(md, encoding="utf-8")
    print(json.dumps({"output": str(args.output), "board_sha256": board_hash, "worst_drop": worst_drop, "base_thermal": base_thermal}, indent=2))


if __name__ == "__main__":
    main()