#!/usr/bin/env python3
"""Example: run a plain ngspice netlist with kisim -- no KiCad involved.

Part 1: a single transient run -- 5 V step into an RC low-pass.
Part 2: tolerance-corner sweep of R1/C1 (kisim corners + run_cases), reading
        the pole of H(s) = 1/(1 + s*R*C) off each step response: a first-order
        step reaches 1 - 1/e of its final value at t = tau = R*C, so
        s = -1/tau and -3 dB at fc = 1/(2*pi*tau).
Part 3: the same as a control model with scipy.signal -- analytic step
        response over the ngspice curves, poles, and a Bode plot.

    micromamba run -n kicad_work python examples/rc_lp.py
"""
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parents[1]))

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

from kisim import Case, KiSim, corners, run_cases, tol, worst

sim = KiSim("rc_lp.cir", tran="1u 2m")   # .cir input: used as-is, no kicad-cli

# --------------------------------------------------------------- single run
r = sim.run()                            # -> results/run.cir + run.raw (cached)

tau = 1e3 * 100e-9                       # R1 * C1
vout = r.v("vout")
print(f"5 V step into 1 kohm + 100 nF: tau = {tau * 1e6:.0f} us")
print(f"V(vout) at 1 tau: {np.interp(tau, r.t, vout):.3f} V   (1 - 1/e of 5 V)")
print(f"V(vout) settled:  {vout[-1]:.3f} V")

r.plot(["vin", "vout"], title="RC low-pass: 5 V step", file=sim.workdir / "rc_lp.png")


# ---------------------------------------------------------- tolerance corners
def tau_fc(r):
    """Read tau off the step response (the 1 - 1/e crossing), then the pole:
    s = -1/tau, i.e. -3 dB at fc = 1/(2*pi*tau)."""
    v = r.v("vout")
    tau = float(np.interp(0.632 * v[-1], v, r.t))
    return {"tau (us)": tau * 1e6, "fc -3dB (kHz)": 1e-3 / (2 * np.pi * tau)}


nominal = Case("nominal", analyses=[tau_fc],
               limits={"fc -3dB (kHz)": lambda fc: 1.3 <= fc <= 2.0})
wca = corners(nominal, {
    "R1": tol("1k", 0.05),        # +/- 5 %  -> tau 85.5 .. 115.5 us,
    "C1": tol("100n", 0.10),      # +/- 10 %     fc 1.38 .. 1.86 kHz
}, mode="extremes")               # every min/max combination: 2^2 runs + nominal
res = run_cases(sim, [nominal] + wca)

f_lo = worst(res, "fc -3dB (kHz)", "min", prefix="nominal__")
f_hi = worst(res, "fc -3dB (kHz)", "max", prefix="nominal__")
print(f"\npole spread over the corners: fc = {f_lo[1]:.2f} .. {f_hi[1]:.2f} kHz  (s = -1/tau)")


def corner_label(name):
    return "nominal" if name == "nominal" else name.split("__", 1)[1].replace("__", ", ")


def rc_of(name):
    """R1/C1 back out of the corner tag, e.g. nominal__R1=950__C1=9e-08."""
    if name == "nominal":
        return 1e3, 100e-9
    vals = dict(kv.split("=") for kv in name.split("__")[1:])
    return float(vals["R1"]), float(vals["C1"])


# ngspice step responses (solid) with the analytic model overlaid (dashed)
fig, ax = plt.subplots(figsize=(8.5, 4.5))
for name, entry in res.items():
    rr = entry["result"]
    R1c, C1c = rc_of(name)
    t_model = np.linspace(0, rr.t[-1], 2001)
    line = ax.plot(rr.t * 1e6, rr.v("vout"), label=corner_label(name))[0]
    _, y_model = signal.step(signal.TransferFunction([1.0], [R1c * C1c, 1.0]), T=t_model)
    ax.plot(t_model * 1e6, 5 * y_model, "--", lw=1, alpha=0.65,  # unit step -> 5 V drive
            color=line.get_color())
ax.plot([np.nan], [np.nan], "k--", lw=1, alpha=0.65, label="scipy model")
ax.axhline(0.632 * 5, color="grey", ls=":", lw=1)
ax.text(1950, 3.35, "63.2 % = 1 - 1/e", fontsize=8, color="grey", ha="right")
ax.set(xlabel="time (us)", ylabel="V(vout) (V)",
       title="RC low-pass at R1/C1 tolerance corners (dashed: scipy.signal model)")
ax.grid(alpha=0.3)
ax.legend(fontsize=8)
fig.tight_layout()
fig.savefig(sim.workdir / "rc_lp_corners.png", dpi=110)
print("plot:", sim.workdir / "rc_lp_corners.png")

# --------------------------------------------------- s-domain / control view
print("\n  corner                          pole (rad/s)    fc (kHz)   tau sim   tau = RC (us)")
for name in res:
    R1c, C1c = rc_of(name)
    sysc = signal.TransferFunction([1.0], [R1c * C1c, 1.0])
    p = sysc.poles[0]
    print(f"  {name:30s} {p:13.1f} {abs(p) / (2 * np.pi * 1e3):9.3f}"
          f" {res[name]['metrics']['tau (us)']:8.1f} {R1c * C1c * 1e6:11.1f}")

w = 2 * np.pi * np.logspace(1.5, 5, 300)          # 30 Hz .. 100 kHz
fig, (axm, axp) = plt.subplots(2, 1, figsize=(8.5, 5.5), sharex=True)
for name in res:
    R1c, C1c = rc_of(name)
    wm, mag, phase = signal.bode(signal.TransferFunction([1.0], [R1c * C1c, 1.0]), w)
    f_khz = wm / (2 * np.pi) / 1e3
    axm.semilogx(f_khz, mag, label=corner_label(name))
    axp.semilogx(f_khz, phase)
axm.axhline(-3.01, color="grey", ls=":", lw=1)
axm.text(12, -4.0, "-3 dB", fontsize=8, color="grey", ha="center", va="top")
axm.set_ylabel("|H(s)|  (dB)")
axp.set_ylabel("phase  (deg)")
axp.set_xlabel("frequency  (kHz)")
axm.set_title("s-domain view: H(s) = 1 / (1 + s R C) per corner")
axm.grid(alpha=0.3, which="both")
axp.grid(alpha=0.3, which="both")
axm.legend(fontsize=8)
fig.tight_layout()
fig.savefig(sim.workdir / "rc_lp_bode.png", dpi=110)
print("plot:", sim.workdir / "rc_lp_bode.png")