1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
#!/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")
|
