"""v0 SPICE demo of the clockless EP cell (single summing column, behavioral op-amp).
Two experiments against the Stage-1 kill criteria:
  E1: system settle with 64 weight branches -> t_settle(0.1%) vs the <=3us commit / >7us kill gate
  E2: two-phase (free vs nudged) differential read cancels a deliberate 5 mV op-amp offset
      -> the zero-reference story, quantified in SPICE.
Behavioral models (v0): op-amp = single-pole VCVS (A0=1e5, GBW ~ 8 MHz, MCP6022-class),
MDAC branch = 100k weight resistor, summing-node parasitic 60 pF lumped. v1 swaps in vendor
macromodels (TI/ADI .subckt) unchanged elsewhere.
"""
import os
os.environ.setdefault('NGSPICE_LIBRARY_PATH', '/home/yurenh2/miniconda3/lib/libngspice.so')
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from PySpice.Spice.Netlist import Circuit, SubCircuit
from PySpice.Unit import *

N_BRANCH = 64
R_W = 100e3        # per-branch weight resistor (MDAC ladder scale)
R_F = 10e3         # TIA feedback
C_NODE = 60e-12    # summing-node lumped parasitic (DAC outputs + bus)
A0 = 1e5           # op-amp DC gain
F_P = 80.0         # dominant pole (Hz) -> GBW = A0 * F_P = 8 MHz
V_OS = 5e-3        # deliberate op-amp input offset (E2)
I_NUDGE = 2e-6     # phase-2 nudge current into the node (beta-scaled)

class OpAmp(SubCircuit):
    NODES = ('inp', 'inn', 'out')
    def __init__(self, name, A0=A0, fp=F_P, vos=0.0):
        super().__init__(name, *self.NODES)
        # offset in series with inp; single-pole gain; 50-ohm output
        self.V('os', 'inp', 'inpo', vos)
        self.B('gain', 'x', self.gnd, v=f'{A0}*(v(inpo)-v(inn))')
        rp = 1e6
        cp = 1.0 / (2 * np.pi * fp * rp)
        self.R('p', 'x', 'p1', rp)
        self.C('p', 'p1', self.gnd, cp)
        self.B('buf', 'o', self.gnd, v='v(p1)')
        self.R('out', 'o', 'out', 50)

def build(vos=0.0, nudge=False, vin=0.05):
    c = Circuit('ep_cell_v0')
    c.subcircuit(OpAmp('opamp', vos=vos))
    # 64 branches step from 0 to vin at t=1us (worst case, all together)
    for i in range(N_BRANCH):
        c.PulseVoltageSource(f'in{i}', f'n{i}', c.gnd,
                             initial_value=0, pulsed_value=vin,
                             delay_time=1e-6, rise_time=5e-9, fall_time=5e-9,
                             pulse_width=1, period=2)
        c.R(f'w{i}', f'n{i}', 'sum', R_W)
    c.C('node', 'sum', c.gnd, C_NODE)
    c.R('f', 'out', 'sum', R_F)
    c.X('amp', 'opamp', c.gnd, 'sum', 'out')   # inverting TIA: inp=gnd, inn=sum
    if nudge:
        c.PulseVoltageSource('nud', 'nn', c.gnd, initial_value=0, pulsed_value=1.0,
                             delay_time=1e-6, rise_time=5e-9, fall_time=5e-9,
                             pulse_width=1, period=2)
        c.R('nudge', 'nn', 'sum', 1.0 / I_NUDGE)   # ~beta-scaled current into the node
    return c

def settled_value_and_time(t, v, tol_frac=1e-3):
    vf = v[-1]
    err = np.abs(v - vf)
    band = tol_frac * abs(vf) if abs(vf) > 1e-9 else tol_frac
    outside = np.where(err > band)[0]
    t_set = t[outside[-1] + 1] if len(outside) and outside[-1] + 1 < len(t) else t[0]
    return vf, t_set

def run(c):
    sim = c.simulator(temperature=27, nominal_temperature=27)
    an = sim.transient(step_time=2e-9, end_time=12e-6)
    t = np.array(an.time)
    vout = np.array(an['out'])
    return t, vout

# ---------- E1: settle ----------
t1, v1 = run(build(vos=0.0, nudge=False))
vf, t_settle = settled_value_and_time(t1, v1)
t_settle_us = (t_settle - 1e-6) * 1e6
ideal = -N_BRANCH * 0.05 * R_F / R_W
print(f'E1  settled Vout {vf*1000:.3f} mV (ideal {ideal*1000:.3f} mV) | '
      f't_settle(0.1%) = {t_settle_us:.3f} us  [commit <=3us, kill >7us]')

# ---------- E2: offset cancellation via two-phase differential ----------
_, v_free = run(build(vos=V_OS, nudge=False))
t2, v_nud = run(build(vos=V_OS, nudge=True))
vf_free, _ = settled_value_and_time(t1, v_free)
vf_nud, _ = settled_value_and_time(t2, v_nud)
ideal_diff = -I_NUDGE * R_F
single_err = abs(vf_free - ideal)
diff_err = abs((vf_nud - vf_free) - ideal_diff)
print(f'E2  offset {V_OS*1000:.1f} mV | single-read error {single_err*1e3:.4f} mV | '
      f'two-phase differential error {diff_err*1e6:.3f} uV | cancellation x{single_err/max(diff_err,1e-12):.0f}')

# ---------- figure ----------
fig, axes = plt.subplots(1, 2, figsize=(11, 4.2))
axes[0].plot(t1 * 1e6, v1 * 1e3, color='#2c6fbb')
axes[0].axvline(1 + t_settle_us, color='#d95f02', ls='--', lw=1)
axes[0].axvspan(1, 4, color='#2e7d32', alpha=0.07)
axes[0].text(1 + t_settle_us, v1.min() * 1e3, f'  settle {t_settle_us:.2f} µs', color='#d95f02', fontsize=9)
axes[0].set_xlabel('t (µs)'); axes[0].set_ylabel('V_out (mV)')
axes[0].set_title(f'E1: 64-branch column settle (gate: ≤3 µs commit)')
axes[1].bar(['single read\n(5 mV offset)', 'two-phase\ndifferential'],
            [single_err * 1e3, diff_err * 1e3], color=['#b03a2e', '#2e7d32'])
axes[1].set_ylabel('|error| (mV)'); axes[1].set_yscale('log')
axes[1].set_title(f'E2: offset rejection ×{single_err/max(diff_err,1e-12):.0f} — the free phase is the zero-reference')
fig.suptitle('Clockless EP cell v0 — PySpice/ngspice (behavioral op-amp, 100k/10k/60pF)', fontsize=11)
fig.tight_layout()
fig.savefig('/home/yurenh2/ept/assets/figs/fig_spice_v0.png', dpi=150)
print('DONE_SPICE_V0')