"""v1: datasheet-anchored column sim.
AD7528 channel = datasheet Thevenin equivalent: V_thev = VREF*N/256, R_eq(code) in [0.8R, 2R],
R = 11k typ; C_OUT 50-120 pF code-dependent -> 64 outputs share one summing bus:
C_bus = 64 x C_OUT ~ 3.2-7.7 nF  (the v0 60 pF assumption was ~100x optimistic).
TIA = MCP6022-class linear model (GBW 10 MHz, A0 1e5, Vos up to 500 uV).
E1: settle vs feedback R_f and compensation C_f -> find the (R_f, C_f) that meets the
    <=3 us commit gate, or report the gate FAILS at the naive design point.
E2: two-phase offset rejection at the chosen operating point with Vos = 500 uV (datasheet max).
Sources: AD7528 datasheet (R typ 11k, Req 0.8R-2R, COUT 50-120 pF, settling 350/400 ns max);
MCP6022 (GBW 10 MHz, Vos max +-500 uV).
"""
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

N = 64
R_LAD = 11e3
R_EQ = 1.3 * R_LAD        # mid-code equivalent (0.8R..2R) -> 14.3k
C_OUT = 85e-12            # mid-code per-chip output capacitance
C_BUS = N * C_OUT         # 5.44 nF on the summing bus
A0 = 1e5
GBW = 10e6
F_P = GBW / A0            # 100 Hz dominant pole
VREF_STEP = 0.1           # per-branch Thevenin step (VREF*code/256 scale)

class OpAmp(SubCircuit):
    NODES = ('inp', 'inn', 'out')
    def __init__(self, name, vos=0.0):
        super().__init__(name, *self.NODES)
        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 * F_P * 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', 25)

def build(rf, cf, vos=0.0, nudge=False):
    c = Circuit('ep_col_v1')
    c.subcircuit(OpAmp('opamp', vos=vos))
    for i in range(N):
        c.PulseVoltageSource(f'in{i}', f'n{i}', c.gnd, initial_value=0, pulsed_value=VREF_STEP,
                             delay_time=1e-6, rise_time=50e-9, fall_time=50e-9,
                             pulse_width=1, period=2)
        c.R(f'w{i}', f'n{i}', 'sum', R_EQ)
    c.C('bus', 'sum', c.gnd, C_BUS)
    c.R('f', 'out', 'sum', rf)
    if cf > 0:
        c.C('f', 'out', 'sum', cf)
    c.X('amp', 'opamp', c.gnd, 'sum', 'out')
    if nudge:
        c.PulseVoltageSource('nud', 'nn', c.gnd, initial_value=0, pulsed_value=1.0,
                             delay_time=1e-6, rise_time=50e-9, fall_time=50e-9,
                             pulse_width=1, period=2)
        c.R('nudge', 'nn', 'sum', 500e3)   # 2 uA nudge
    return c

def run(c, end=40e-6, step=4e-9):
    sim = c.simulator(temperature=27, nominal_temperature=27)
    an = sim.transient(step_time=step, end_time=end)
    return np.array(an.time), np.array(an['out'])

def settle_time(t, v, tol=1e-3):
    vf = v[-1]
    band = tol * abs(vf) if abs(vf) > 1e-9 else tol
    outside = np.where(np.abs(v - vf) > band)[0]
    if len(outside) == 0: return 0.0
    idx = outside[-1] + 1
    return (t[idx] - 1e-6) if idx < len(t) else np.inf

# ---------- E1: (R_f, C_f) design sweep ----------
print('E1: settle vs (R_f, C_f)   [commit <=3us | kill >7us]  C_bus = %.2f nF' % (C_BUS * 1e9))
results = []
for rf in (10e3, 3e3, 1e3):
    # optimal TIA comp: C_f ~ sqrt(C_bus / (2*pi*GBW*R_f))
    cf_opt = np.sqrt(C_BUS / (2 * np.pi * GBW * rf))
    for cf in (0.0, cf_opt, 2 * cf_opt):
        t, v = run(build(rf, cf))
        ts = settle_time(t, v)
        ideal = -N * VREF_STEP * rf / R_EQ
        err = abs(v[-1] - ideal) / abs(ideal)
        results.append((rf, cf, ts, v[-1], err))
        print(f'  R_f={rf/1e3:4.0f}k  C_f={cf*1e12:7.1f}pF  t_settle={ts*1e6:7.2f}us  '
              f'Vout={v[-1]*1e3:8.2f}mV  static_err={err*100:.3f}%')

ok = [r for r in results if r[2] <= 3e-6 and r[4] < 0.01]
best = min(ok, key=lambda r: r[2]) if ok else min(results, key=lambda r: r[2])
verdict = 'COMMIT GATE MET' if ok else 'GATE FAILED at all tested points'
print(f'E1 verdict: {verdict} -> best (R_f={best[0]/1e3:.0f}k, C_f={best[1]*1e12:.0f}pF, t={best[2]*1e6:.2f}us)')

# ---------- E2: offset rejection at the chosen point ----------
rf_b, cf_b = best[0], best[1]
_, v_free = run(build(rf_b, cf_b, vos=500e-6))
t2, v_nud = run(build(rf_b, cf_b, vos=500e-6, nudge=True))
ideal = -N * VREF_STEP * rf_b / R_EQ
ideal_diff = -2e-6 * rf_b
single_err = abs(v_free[-1] - ideal)
diff_err = abs((v_nud[-1] - v_free[-1]) - ideal_diff)
print(f'E2: Vos=500uV | single-read err {single_err*1e3:.3f} mV | two-phase diff err '
      f'{diff_err*1e6:.2f} uV | rejection x{single_err/max(diff_err,1e-12):.0f} | '
      f'contrast floor 0.1-1 mV -> {"OK" if diff_err < 1e-4 else "MARGINAL/FAIL"}')

# ---------- figure ----------
fig, axes = plt.subplots(1, 2, figsize=(12, 4.4))
for rf in (10e3, 3e3, 1e3):
    cf_opt = np.sqrt(C_BUS / (2 * np.pi * GBW * rf))
    t, v = run(build(rf, cf_opt))
    axes[0].plot(t * 1e6, v * 1e3, label=f'R_f={rf/1e3:.0f}k, C_f={cf_opt*1e12:.0f}pF '
                 f'({settle_time(t,v)*1e6:.2f}µs)')
axes[0].axvspan(1, 4, color='#2e7d32', alpha=0.07)
axes[0].axvline(8, color='#b03a2e', ls=':', lw=1)
axes[0].set_xlabel('t (µs)'); axes[0].set_ylabel('V_out (mV)')
axes[0].set_title(f'E1: 64-chip column, C_bus={C_BUS*1e9:.1f}nF (datasheet C_OUT) — comp-C TIA')
axes[0].legend(fontsize=8)
axes[1].bar(['single read\n(500 µV Vos, ×noise-gain)', 'two-phase\ndifferential'],
            [single_err * 1e3, diff_err * 1e3], color=['#b03a2e', '#2e7d32'])
axes[1].set_yscale('log'); axes[1].set_ylabel('|error| (mV)')
axes[1].axhline(0.1, color='#888', ls='--', lw=1)
axes[1].text(0.5, 0.105, 'contrast signal floor 0.1 mV', fontsize=8, color='#666')
axes[1].set_title(f'E2 @ (R_f={rf_b/1e3:.0f}k): rejection ×{single_err/max(diff_err,1e-12):.0f}')
fig.suptitle('EP column v1 — AD7528 Thevenin + bus C + MCP6022-class TIA (datasheet-anchored)', fontsize=11)
fig.tight_layout()
fig.savefig('/home/yurenh2/ept/assets/figs/fig_spice_v1.png', dpi=150)
print('DONE_SPICE_V1')