"""v1.1: (A) ENOB-after-per-cell-cal Monte Carlo (resistor lot spread + code-dependent R_eq),
(B) op-amp noise integrated over the closed-loop bandwidth at the v1 operating point
    (R_f=1k, C_f=294pF) vs the 0.1-1 mV contrast floor.
DC transfer is algebraic (settled state); SPICE noise via ngspice .noise on the TIA.
"""
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

rng = np.random.default_rng(7)
N = 64
R0 = 11e3
RF, CF = 1e3, 294e-12
EN_V = 8.7e-9          # MCP6022 input voltage noise, V/rtHz (flatband)
FS_CODE = 255

def req_of_code(code, rbase):
    # datasheet: equivalent output R varies 0.8R..2R with code (excluding all-0s).
    # smooth monotone proxy between the endpoints:
    frac = max(code, 1) / 255.0
    return rbase * (2.0 - 1.2 * frac)

def settled_out(codes, rbases, vref=0.1, rf=RF):
    # inverting summing node: Vout = -rf * sum(vref*code_i/256 / Req_i(code_i))
    return -rf * sum(vref * (c / 256.0) / req_of_code(c, rb) for c, rb in zip(codes, rbases))

# ---------- A: per-BIT ladder mismatch -> INL not removable by single-scale cal ----------
# each DAC's 8 bit-legs carry independent relative errors eps_k; per-cell cal removes only the
# overall gain; the residual code-shape error IS the INL. sigma_bit set so typical INL ~ datasheet
# relative accuracy (+-1 LSB class).
NMC = 400
SIG_BIT = 0.004
enob = []
for _ in range(NMC):
    eps = rng.normal(0, SIG_BIT, (N, 8))
    w = (2 ** np.arange(8))[None, :] * (1 + eps)          # per-cell bit weights
    codes = np.arange(1, 256)
    bits = ((codes[:, None] >> np.arange(8)[None, :]) & 1)  # 255 x 8
    val = bits @ w.T                                      # 255 x N raw values
    ideal = codes.astype(float)[:, None]
    cal = ideal[-1] / val[-1]                             # per-cell single-scale cal at code 255
    resid = val * cal[None, :] - ideal                    # LSB units
    rms = np.sqrt(np.mean(resid ** 2))
    enob.append(8 - np.log2(max(2 * np.sqrt(3) * rms / 1.0, 1e-9)) - 0)
enob = np.array(enob)
inl_p95 = None
print(f'A: per-bit mismatch sigma {SIG_BIT*100:.1f}% -> ENOB after per-cell cal: '
      f'mean {enob.mean():.2f} b, p5 {np.percentile(enob,5):.2f} b  [screen bar >=7.0]')

# ---------- B: ngspice noise analysis at the operating point ----------
class OpAmp(SubCircuit):
    NODES = ('inp', 'inn', 'out')
    def __init__(self, name):
        super().__init__(name, *self.NODES)
        self.B('gain', 'x', self.gnd, v=f'1e5*(v(inp)-v(inn))')
        rp = 1e6
        cp = 1.0 / (2 * np.pi * 100.0 * 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)

c = Circuit('noise_col')
c.subcircuit(OpAmp('opamp'))
c.SinusoidalVoltageSource('vn', 'np', c.gnd, amplitude=1.0)   # noise injection point (input-referred)
for i in range(N):
    c.R(f'w{i}', c.gnd, 'sum', req_of_code(128, R0))
c.C('bus', 'sum', c.gnd, 64 * 85e-12)
c.R('f', 'out', 'sum', RF)
c.C('f', 'out', 'sum', CF)
c.X('amp', 'opamp', 'np', 'sum', 'out')   # inject at the + input = input-referred source
sim = c.simulator(temperature=27, nominal_temperature=27)
an = sim.ac(start_frequency=10, stop_frequency=100e6, number_of_points=30, variation='dec')
f = np.array(an.frequency)
gain = np.abs(np.array(an['out']))        # noise gain vs frequency
# integrate (EN_V * gain)^2 df
integrand = (EN_V * gain) ** 2
vn_rms = np.sqrt(np.trapz(integrand, f))
nbw_gain = np.sqrt(np.trapz(gain ** 2, f) / max(f[-1] - f[0], 1))
print(f'B: output-referred op-amp noise RMS {vn_rms*1e6:.2f} uV over closed-loop band '
      f'(peak noise gain {gain.max():.1f}) | contrast floor 0.1 mV -> '
      f'{"OK" if vn_rms < 1e-4 else "MARGINAL"} (ratio {1e-4/vn_rms:.1f}x headroom)')

# ---------- figure ----------
fig, axes = plt.subplots(1, 2, figsize=(11.5, 4.2))
axes[0].hist(enob, bins=30, color='#2c6fbb', alpha=0.85)
axes[0].axvline(7.0, color='#b03a2e', ls='--')
axes[0].text(7.0, axes[0].get_ylim()[1]*0.9, ' screen bar 7.0', color='#b03a2e', fontsize=9)
axes[0].set_xlabel('ENOB proxy after per-cell cal (bits)'); axes[0].set_ylabel('MC count')
axes[0].set_title(f'A: 400-trial MC, per-bit mismatch 0.4% — mean {enob.mean():.2f} b')
axes[1].loglog(f, EN_V * gain * 1e9, color='#7b5aa6')
axes[1].set_xlabel('Hz'); axes[1].set_ylabel('output noise density (nV/√Hz)')
axes[1].set_title(f'B: op-amp noise → output; integrated RMS {vn_rms*1e6:.1f} µV (floor 100 µV)')
fig.suptitle('EP column v1.1 — ENOB-after-cal Monte Carlo + noise budget (datasheet-anchored)', fontsize=11)
fig.tight_layout()
fig.savefig('/home/yurenh2/ept/assets/figs/fig_spice_v11.png', dpi=150)
print('DONE_SPICE_V11')