pixalynx/nrf9151-gps-tracker

scripts/supply-resistance.py

#!/usr/bin/env python3
"""DC resistance between two lands through the routed copper (nodal analysis; every trace is a resistor between its
end lands, 1 oz outer / 0.5 oz inner copper, 1 mOhm per via). Parallel paths count, so it answers "how much does the
supply drop" rather than "how long is the shortest trace".
Usage: python3 scripts/supply-resistance.py U3.VSYS:FB1.A J6.BAT_P:U3.VBAT ...
"""
import json, math, sys
d=json.load(open('dist/index/circuit.json'))
src={c['source_component_id']:c['name'] for c in d if c['type']=='source_component'}
sp={p['source_port_id']:p for p in d if p['type']=='source_port'}; pp={p['pcb_port_id']:p for p in d if p['type']=='pcb_port'}
def pname(pid): s=sp[pp[pid]['source_port_id']]; return f"{src[s['source_component_id']]}.{s['name']}"
edges=[]
for t in d:
    if t['type']!='pcb_trace': continue
    r=t['route']; ends=[pname(p[k]) for p in r for k in ('start_pcb_port_id','end_pcb_port_id') if p.get(k)]
    if len(ends)!=2 or ends[0]==ends[1]: continue
    R=0
    for a,b in zip(r,r[1:]):
        layer = a.get('layer') if a['route_type']=='wire' else b.get('layer') if b['route_type']=='wire' else None
        if not layer: continue
        L=math.hypot(a['x']-b['x'],a['y']-b['y']); w=(a.get('width') if a['route_type']=='wire' else None) or b.get('width')
        R+=1.72e-8*L*1e-3/(w*1e-3*(17.5e-6 if layer.startswith('inner') else 35e-6))
    vias=len({(round(p['x'],3),round(p['y'],3)) for p in r if p['route_type']=='via'})
    R+=0.001*vias
    edges.append((ends[0],ends[1],R))
def eff(s,g):
    # nodal analysis: inject 1 A at s, ground g, solve G v = i
    nodes=sorted({e[0] for e in edges}|{e[1] for e in edges}); comp={}
    # keep only the connected component of s
    adj={}
    for a,b,R in edges: adj.setdefault(a,[]).append(b); adj.setdefault(b,[]).append(a)
    seen={s}; st=[s]
    while st:
        u=st.pop()
        for v in adj.get(u,[]):
            if v not in seen: seen.add(v); st.append(v)
    if g not in seen: return None
    idx={n:i for i,n in enumerate(sorted(seen - {g}))}; n=len(idx)
    G=[[0.0]*n for _ in range(n)]; I=[0.0]*n; I[idx[s]]=1.0
    for a,b,R in edges:
        if a not in seen: continue
        c=1/max(R,1e-6)
        for x,y in ((a,b),(b,a)):
            if x in idx:
                G[idx[x]][idx[x]]+=c
                if y in idx: G[idx[x]][idx[y]]-=c
    for i in range(n):  # Gauss-Jordan
        p=max(range(i,n),key=lambda r: abs(G[r][i])); G[i],G[p]=G[p],G[i]; I[i],I[p]=I[p],I[i]
        for r in range(n):
            if r!=i and G[r][i]:
                f=G[r][i]/G[i][i]
                for k in range(i,n): G[r][k]-=f*G[i][k]
                I[r]-=f*I[i]
    return I[idx[s]]/G[idx[s]][idx[s]]
for s,g in [a.split(':') for a in sys.argv[1:]]:
    r=eff(s,g); print(f"{s} -> {g}: {'n/a' if r is None else f'{r*1000:.0f} mOhm'}")