RITHESH LED Display

"""Estimate the USB D+/D- differential impedance with a small 2D field solver (finite differences).

Cross-section of the JLCPCB standard 2-layer 1.6 mm board:
  ground plane (other layer) | FR-4 1.53 mm, er 4.5 | 35 um copper traces | ~20 um solder mask er 3.8 | air.
The pair: 0.30 mm traces, 0.20 mm gap (as routed). Two cases:
  A: pair over the ground plane only (no ground copper beside it),
  B: pair with ground fill on the same layer 0.25 mm from each trace (as along most of this board's USB route).
For each case the capacitance per metre is found from the field energy with and without the dielectrics,
giving odd/even-mode impedances; Zdiff = 2 * Zodd, Zcommon = Zeven / 2.

Usage: python3 usb_impedance.py
"""
import numpy as np
from scipy.sparse import coo_matrix, diags
from scipy.sparse.linalg import cg

E0 = 8.854e-12
C0 = 299792458.0
D = 0.02                     # grid step, mm (0.01 mm gives ~5 ohm lower; see notes)
H, T, MASK = 1.53, 0.035, 0.02
ER_FR4, ER_MASK = 4.5, 3.8
W, S, G = 0.30, 0.20, 0.25
XHALF, YTOP = 10.0, 8.0     # far boundaries: results converged at this size


def solve(case, with_dielectric, mode):
    nx, ny = int(2 * XHALF / D) + 1, int(YTOP / D) + 1
    x = -XHALF + np.arange(nx) * D
    y = np.arange(ny) * D                       # y = 0 is the ground plane surface
    X, Y = np.meshgrid(x, y)
    er = np.ones((ny, nx))
    if with_dielectric:
        er[Y < H] = ER_FR4
        top = H + T + MASK
        er[(Y >= H) & (Y < H + MASK)] = ER_MASK                # mask on the laminate surface
    fixed = np.zeros((ny, nx), bool)
    val = np.zeros((ny, nx))
    fixed[0, :] = True                                          # ground plane
    fixed[-1, :] = fixed[:, 0] = fixed[:, -1] = True            # far boundaries at 0 V
    in_cu = (Y >= H) & (Y <= H + T)
    a1 = (X >= -S / 2 - W) & (X <= -S / 2) & in_cu
    a2 = (X >= S / 2) & (X <= S / 2 + W) & in_cu
    if with_dielectric:                                         # mask over the trace tops
        over = (Y > H + T) & (Y <= H + T + MASK) & (((X >= -S / 2 - W) & (X <= -S / 2)) | ((X >= S / 2) & (X <= S / 2 + W)))
        er[over] = ER_MASK
    fixed[a1] = True
    fixed[a2] = True
    val[a1] = 1.0
    val[a2] = -1.0 if mode == "odd" else 1.0
    if case == "B":                                             # coplanar ground fill 0.25 mm away
        g = in_cu & ((X <= -S / 2 - W - G) | (X >= S / 2 + W + G))
        fixed[g] = True
        val[g] = 0.0
    idx = -np.ones((ny, nx), int)
    free = ~fixed
    idx[free] = np.arange(free.sum())
    n = int(free.sum())
    rows, cols, vals = [], [], []
    rhs = np.zeros(n)
    for di, dj in ((0, 1), (0, -1), (1, 0), (-1, 0)):
        i0, i1 = max(0, -di), ny - max(0, di)
        j0, j1 = max(0, -dj), nx - max(0, dj)
        a = (slice(i0, i1), slice(j0, j1))
        b = (slice(i0 + di, i1 + di), slice(j0 + dj, j1 + dj))
        eps = 2 * er[a] * er[b] / (er[a] + er[b])                # harmonic mean at the face
        m = free[a]
        ia = idx[a][m]
        e = eps[m]
        rows.append(ia)
        cols.append(ia)
        vals.append(e)
        nb_free = free[b][m]
        rows.append(ia[nb_free])
        cols.append(idx[b][m][nb_free])
        vals.append(-e[nb_free])
        np.add.at(rhs, ia[~nb_free], e[~nb_free] * val[b][m][~nb_free])
    A = coo_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))), shape=(n, n)).tocsr()
    phi_f, info = cg(A, rhs, tol=1e-9, maxiter=20000, M=diags(1.0 / A.diagonal()))
    phi = val.copy()
    phi[free] = phi_f
    # field energy per metre: W = 1/2 * sum(eps * |grad phi|^2) * cell area (grid in mm cancels: 2D)
    ex = np.diff(phi, axis=1) / D
    ey = np.diff(phi, axis=0) / D
    erx = 0.5 * (er[:, 1:] + er[:, :-1])
    ery = 0.5 * (er[1:, :] + er[:-1, :])
    Wm = 0.5 * E0 * ((erx * ex ** 2).sum() + (ery * ey ** 2).sum()) * D * D
    return Wm                                                  # joules per metre with 1 V on the traces


def impedance(case):
    out = {}
    for mode in ("odd", "even"):
        c = solve(case, True, mode)      # C per trace = W for +/-1 V (odd) or +1/+1 V (even), see docstring
        c0 = solve(case, False, mode)
        z = 1.0 / (C0 * np.sqrt(c * c0))
        out[mode] = (z, c / c0)
    return out


def main():
    print(f"cross-section: FR-4 {H} mm (er {ER_FR4}), 35 um copper, {MASK * 1000:.0f} um mask (er {ER_MASK}); "
          f"traces {W} mm, gap {S} mm")
    for case, label in (("A", "over the ground plane only"), ("B", f"with same-layer ground fill {G} mm away")):
        r = impedance(case)
        zo, ze = r["odd"][0], r["even"][0]
        print(f"case {case} ({label}): Zodd {zo:.1f} ohm, Zeven {ze:.1f} ohm -> "
              f"Zdiff {2 * zo:.0f} ohm, Zcommon {ze / 2:.0f} ohm (effective er {r['odd'][1]:.2f})")


if __name__ == "__main__":
    main()