diff --git a/.github/workflows/bun-formatcheck.yml b/.github/workflows/bun-formatcheck.yml
new file mode 100644
index 0000000..be449cf
--- /dev/null
+++ b/.github/workflows/bun-formatcheck.yml
@@ -0,0 +1,26 @@
+# Created using @tscircuit/plop (npm install -g @tscircuit/plop)
+name: Format Check
+
+on:
+  push:
+    branches: [main]
+  pull_request:
+    branches: [main]
+
+jobs:
+  format-check:
+    runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Setup bun
+        uses: oven-sh/setup-bun@v2
+        with:
+          bun-version: latest
+
+      - name: Install dependencies
+        run: bun install
+
+      - name: Run format check
+        run: bun run format:check
diff --git a/.github/workflows/bun-pver-release.yml b/.github/workflows/bun-pver-release.yml
new file mode 100644
index 0000000..a13f49e
--- /dev/null
+++ b/.github/workflows/bun-pver-release.yml
@@ -0,0 +1,70 @@
+# Created using @tscircuit/plop (npm install -g @tscircuit/plop)
+name: Publish to npm
+on:
+  push:
+    branches:
+      - main
+  workflow_dispatch:
+
+env:
+  UPSTREAM_REPOS: "" # comma-separated list, e.g. "eval,tscircuit,docs"
+  UPSTREAM_PACKAGES_TO_UPDATE: "" # comma-separated list, e.g. "@tscircuit/core,@tscircuit/protos"
+
+jobs:
+  publish:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          token: ${{ secrets.TSCIRCUIT_BOT_GITHUB_TOKEN }}
+      - name: Setup bun
+        uses: oven-sh/setup-bun@v2
+        with:
+          bun-version: latest
+      - uses: actions/setup-node@v3
+        with:
+          node-version: 20
+          registry-url: https://registry.npmjs.org/
+      - run: npm install -g pver
+      - run: bun install --frozen-lockfile
+      - run: bun run build
+      - run: pver release
+        env:
+          NODE_AUTH_TOKEN: ${{ secrets.NPM_TOKEN }}
+          GITHUB_TOKEN: ${{ secrets.TSCIRCUIT_BOT_GITHUB_TOKEN }}
+
+      # - name: Create Pull Request
+      #   id: create-pr
+      #   uses: peter-evans/create-pull-request@v5
+      #   with:
+      #     commit-message: "chore: bump version"
+      #     title: "chore: bump version"
+      #     body: "Automated package update"
+      #     branch: bump-version-${{ github.run_number }}
+      #     base: main
+      #     token: ${{ secrets.TSCIRCUIT_BOT_GITHUB_TOKEN }}
+      #     committer: tscircuitbot <githubbot@tscircuit.com>
+      #     author: tscircuitbot <githubbot@tscircuit.com>
+
+      # - name: Enable auto-merge
+      #   if: steps.create-pr.outputs.pull-request-number != ''
+      #   run: |
+      #     gh pr merge ${{ steps.create-pr.outputs.pull-request-number }} --rebase --delete-branch
+      #   env:
+      #     GH_TOKEN: ${{ secrets.TSCIRCUIT_BOT_GITHUB_TOKEN }}
+
+      # - name: Trigger upstream repo updates
+      #   if: env.UPSTREAM_REPOS && env.UPSTREAM_PACKAGES_TO_UPDATE
+      #   run: |
+      #     IFS=',' read -ra REPOS <<< "${{ env.UPSTREAM_REPOS }}"
+      #     for repo in "${REPOS[@]}"; do
+      #       if [[ -n "$repo" ]]; then
+      #         echo "Triggering update for repo: $repo"
+      #         curl -X POST \
+      #           -H "Accept: application/vnd.github.v3+json" \
+      #           -H "Authorization: token ${{ secrets.TSCIRCUIT_BOT_GITHUB_TOKEN }}" \
+      #           -H "Content-Type: application/json" \
+      #           "https://api.github.com/repos/tscircuit/$repo/actions/workflows/update-package.yml/dispatches" \
+      #           -d "{\"ref\":\"main\",\"inputs\":{\"package_names\":\"${{ env.UPSTREAM_PACKAGES_TO_UPDATE }}\"}}"
+      #       fi
+      #     done
\ No newline at end of file
diff --git a/.github/workflows/bun-test.yml b/.github/workflows/bun-test.yml
new file mode 100644
index 0000000..e659651
--- /dev/null
+++ b/.github/workflows/bun-test.yml
@@ -0,0 +1,31 @@
+# Created using @tscircuit/plop (npm install -g @tscircuit/plop)
+name: Bun Test
+
+on:
+  pull_request:
+  push:
+    branches:
+      - main
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    timeout-minutes: 5
+
+    # Skip test for PRs that not chore: bump version
+    if: "${{ github.event_name != 'pull_request' || github.event.pull_request.title != 'chore: bump version' }}"
+
+    steps:
+      - name: Checkout code
+        uses: actions/checkout@v4
+
+      - name: Setup bun
+        uses: oven-sh/setup-bun@v2
+        with:
+          bun-version: latest
+
+      - name: Install dependencies
+        run: bun install
+
+      - name: Run tests
+        run: bun test
diff --git a/.github/workflows/bun-typecheck.yml b/.github/workflows/bun-typecheck.yml
new file mode 100644
index 0000000..68acb7d
--- /dev/null
+++ b/.github/workflows/bun-typecheck.yml
@@ -0,0 +1,26 @@
+# Created using @tscircuit/plop (npm install -g @tscircuit/plop)
+name: Type Check
+
+on:
+  push:
+    branches: [main]
+  pull_request:
+    branches: [main]
+
+jobs:
+  type-check:
+    runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Setup bun
+        uses: oven-sh/setup-bun@v2
+        with:
+          bun-version: latest
+
+      - name: Install dependencies
+        run: bun i
+
+      - name: Run type check
+        run: bunx tsc --noEmit
diff --git a/README.md b/README.md
index 195df0c..f1e6dfc 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,9 @@
 # spicey
 
+Run [SPICE](https://en.wikipedia.org/wiki/SPICE) simulations in native javascript. An alternative to [ngspice](https://ngspice.sourceforge.io/)
+
+[![npm version](https://img.shields.io/npm/v/spicey.svg)](https://www.npmjs.com/package/spicey)
+
 ```tsx
 import { simulate, formatAcResult } from "spicey"
 
@@ -31,3 +35,20 @@ formatAcResult(result1.ac)
   ..."
 `)
 ```
+
+## Proposed directory structure
+
+To make it easy to extend the simulator (for example to add transistor models later), the library is now organized into focused modules:
+
+```
+lib/
+  analysis/            # High-level simulation entry points (simulate, simulateAC, simulateTRAN)
+  constants/           # Shared numeric constants
+  formatting/          # Result formatting helpers
+  math/                # Numeric utilities such as Complex arithmetic and matrix solvers
+  parsing/             # Netlist parsing and circuit data structures
+  stamping/            # Matrix/RHS stamping helpers for modified nodal analysis
+  utils/               # Generic helpers (e.g., logarithmic sweeps)
+```
+
+Each exported function lives in its own file with a matching name, so new capabilities can be added without creating monolithic modules.
diff --git a/lib/analysis/simulate.ts b/lib/analysis/simulate.ts
new file mode 100644
index 0000000..09a624c
--- /dev/null
+++ b/lib/analysis/simulate.ts
@@ -0,0 +1,12 @@
+import { parseNetlist } from "../parsing/parseNetlist"
+import { simulateAC } from "./simulateAC"
+import { simulateTRAN } from "./simulateTRAN"
+
+function simulate(netlistText: string) {
+  const circuit = parseNetlist(netlistText)
+  const ac = simulateAC(circuit)
+  const tran = simulateTRAN(circuit)
+  return { circuit, ac, tran }
+}
+
+export { simulate }
diff --git a/lib/analysis/simulateAC.ts b/lib/analysis/simulateAC.ts
new file mode 100644
index 0000000..b930652
--- /dev/null
+++ b/lib/analysis/simulateAC.ts
@@ -0,0 +1,114 @@
+import { EPS } from "../constants/EPS"
+import { Complex } from "../math/Complex"
+import { solveComplex } from "../math/solveComplex"
+import type { ParsedCircuit } from "../parsing/parseNetlist"
+import { logspace } from "../utils/logspace"
+import { stampAdmittanceComplex } from "../stamping/stampAdmittanceComplex"
+import { stampVoltageSourceComplex } from "../stamping/stampVoltageSourceComplex"
+
+function simulateAC(ckt: ParsedCircuit) {
+  if (!ckt.analyses.ac) return null
+
+  const { mode, N, f1, f2 } = ckt.analyses.ac
+  const nNodeVars = ckt.nodes.count() - 1
+  const nVsrc = ckt.V.length
+  const Nvar = nNodeVars + nVsrc
+
+  const freqs =
+    mode === "dec"
+      ? logspace(f1, f2, N)
+      : (() => {
+          const arr: number[] = []
+          const npts = Math.max(2, N)
+          const step = (f2 - f1) / (npts - 1)
+          for (let i = 0; i < npts; i++) arr.push(f1 + i * step)
+          return arr
+        })()
+
+  const nodeVoltages: Record<string, Complex[]> = {}
+  ckt.nodes.rev.forEach((name, id) => {
+    if (id !== 0) nodeVoltages[name] = []
+  })
+  const elementCurrents: Record<string, Complex[]> = {}
+
+  const twoPi = 2 * Math.PI
+
+  for (const f of freqs) {
+    const A = Array.from({ length: Nvar }, () =>
+      Array.from({ length: Nvar }, () => Complex.from(0, 0)),
+    )
+    const b = Array.from({ length: Nvar }, () => Complex.from(0, 0))
+
+    for (const r of ckt.R) {
+      if (r.R <= 0) throw new Error(`R ${r.name} must be > 0`)
+      const Y = Complex.from(1 / r.R, 0)
+      stampAdmittanceComplex(A, ckt.nodes, r.n1, r.n2, Y)
+    }
+
+    for (const c of ckt.C) {
+      const Y = Complex.from(0, twoPi * f * c.C)
+      stampAdmittanceComplex(A, ckt.nodes, c.n1, c.n2, Y)
+    }
+
+    for (const l of ckt.L) {
+      const denom = Complex.from(0, twoPi * f * l.L)
+      const Y =
+        denom.abs() < EPS ? Complex.from(0, 0) : Complex.from(1, 0).div(denom)
+      stampAdmittanceComplex(A, ckt.nodes, l.n1, l.n2, Y)
+    }
+
+    for (const vs of ckt.V) {
+      const Vph = Complex.fromPolar(vs.acMag || 0, vs.acPhaseDeg || 0)
+      stampVoltageSourceComplex(A, b, ckt.nodes, vs, Vph)
+    }
+
+    const x = solveComplex(A, b)
+
+    for (let id = 1; id < ckt.nodes.count(); id++) {
+      const idx = id - 1
+      const nodeName = ckt.nodes.rev[id]
+      if (!nodeName) continue
+      const series = nodeVoltages[nodeName]
+      if (!series) continue
+      series.push(x[idx] ?? Complex.from(0, 0))
+    }
+
+    for (const r of ckt.R) {
+      const v1 =
+        r.n1 === 0 ? Complex.from(0, 0) : (x[r.n1 - 1] ?? Complex.from(0, 0))
+      const v2 =
+        r.n2 === 0 ? Complex.from(0, 0) : (x[r.n2 - 1] ?? Complex.from(0, 0))
+      const Y = Complex.from(1 / r.R, 0)
+      const i = Y.mul(v1.sub(v2))
+      ;(elementCurrents[r.name] ||= []).push(i)
+    }
+    for (const c of ckt.C) {
+      const v1 =
+        c.n1 === 0 ? Complex.from(0, 0) : (x[c.n1 - 1] ?? Complex.from(0, 0))
+      const v2 =
+        c.n2 === 0 ? Complex.from(0, 0) : (x[c.n2 - 1] ?? Complex.from(0, 0))
+      const Y = Complex.from(0, twoPi * f * c.C)
+      const i = Y.mul(v1.sub(v2))
+      ;(elementCurrents[c.name] ||= []).push(i)
+    }
+    for (const l of ckt.L) {
+      const v1 =
+        l.n1 === 0 ? Complex.from(0, 0) : (x[l.n1 - 1] ?? Complex.from(0, 0))
+      const v2 =
+        l.n2 === 0 ? Complex.from(0, 0) : (x[l.n2 - 1] ?? Complex.from(0, 0))
+      const denom = Complex.from(0, twoPi * f * l.L)
+      const Y =
+        denom.abs() < EPS ? Complex.from(0, 0) : Complex.from(1, 0).div(denom)
+      const i = Y.mul(v1.sub(v2))
+      ;(elementCurrents[l.name] ||= []).push(i)
+    }
+    for (const vs of ckt.V) {
+      const i = x[vs.index] ?? Complex.from(0, 0)
+      ;(elementCurrents[vs.name] ||= []).push(i)
+    }
+  }
+
+  return { freqs, nodeVoltages, elementCurrents }
+}
+
+export { simulateAC }
diff --git a/lib/analysis/simulateTRAN.ts b/lib/analysis/simulateTRAN.ts
new file mode 100644
index 0000000..60e3986
--- /dev/null
+++ b/lib/analysis/simulateTRAN.ts
@@ -0,0 +1,108 @@
+import { EPS } from "../constants/EPS"
+import { solveReal } from "../math/solveReal"
+import type { ParsedCircuit } from "../parsing/parseNetlist"
+import { stampAdmittanceReal } from "../stamping/stampAdmittanceReal"
+import { stampCurrentReal } from "../stamping/stampCurrentReal"
+import { stampVoltageSourceReal } from "../stamping/stampVoltageSourceReal"
+
+function simulateTRAN(ckt: ParsedCircuit) {
+  if (!ckt.analyses.tran) return null
+  const { dt, tstop } = ckt.analyses.tran
+  const steps = Math.max(1, Math.ceil(tstop / Math.max(dt, EPS)))
+
+  const nNodeVars = ckt.nodes.count() - 1
+  const nVsrc = ckt.V.length
+  const Nvar = nNodeVars + nVsrc
+
+  const times: number[] = []
+  const nodeVoltages: Record<string, number[]> = {}
+  ckt.nodes.rev.forEach((name, id) => {
+    if (id !== 0) nodeVoltages[name] = []
+  })
+  const elementCurrents: Record<string, number[]> = {}
+
+  const blankA = () =>
+    Array.from({ length: Nvar }, () => new Array(Nvar).fill(0))
+  const blankB = () => new Array(Nvar).fill(0)
+
+  let t = 0
+  for (let step = 0; step <= steps; step++, t = step * dt) {
+    times.push(t)
+    const A = blankA()
+    const b = blankB()
+
+    for (const r of ckt.R) {
+      const G = 1 / r.R
+      stampAdmittanceReal(A, ckt.nodes, r.n1, r.n2, G)
+    }
+
+    for (const c of ckt.C) {
+      const Gc = c.C / Math.max(dt, EPS)
+      stampAdmittanceReal(A, ckt.nodes, c.n1, c.n2, Gc)
+      const Ieq = -Gc * c.vPrev
+      stampCurrentReal(b, ckt.nodes, c.n1, c.n2, Ieq)
+    }
+
+    for (const l of ckt.L) {
+      const Gl = Math.max(dt, EPS) / l.L
+      stampAdmittanceReal(A, ckt.nodes, l.n1, l.n2, Gl)
+      stampCurrentReal(b, ckt.nodes, l.n1, l.n2, l.iPrev)
+    }
+
+    for (const vs of ckt.V) {
+      const Vt = vs.waveform ? vs.waveform(t) : vs.dc || 0
+      stampVoltageSourceReal(A, b, ckt.nodes, vs, Vt)
+    }
+
+    const x = solveReal(A, b)
+
+    for (let id = 1; id < ckt.nodes.count(); id++) {
+      const idx = id - 1
+      const nodeName = ckt.nodes.rev[id]
+      if (!nodeName) continue
+      const series = nodeVoltages[nodeName]
+      if (!series) continue
+      series.push(x[idx] ?? 0)
+    }
+
+    for (const r of ckt.R) {
+      const v1 = r.n1 === 0 ? 0 : (x[r.n1 - 1] ?? 0)
+      const v2 = r.n2 === 0 ? 0 : (x[r.n2 - 1] ?? 0)
+      const i = (v1 - v2) / r.R
+      ;(elementCurrents[r.name] ||= []).push(i)
+    }
+    for (const c of ckt.C) {
+      const v1 = c.n1 === 0 ? 0 : (x[c.n1 - 1] ?? 0)
+      const v2 = c.n2 === 0 ? 0 : (x[c.n2 - 1] ?? 0)
+      const i = (c.C * (v1 - v2 - c.vPrev)) / Math.max(dt, EPS)
+      ;(elementCurrents[c.name] ||= []).push(i)
+    }
+    for (const l of ckt.L) {
+      const v1 = l.n1 === 0 ? 0 : (x[l.n1 - 1] ?? 0)
+      const v2 = l.n2 === 0 ? 0 : (x[l.n2 - 1] ?? 0)
+      const Gl = Math.max(dt, EPS) / l.L
+      const i = Gl * (v1 - v2) + l.iPrev
+      ;(elementCurrents[l.name] ||= []).push(i)
+    }
+    for (const vs of ckt.V) {
+      const i = x[vs.index] ?? 0
+      ;(elementCurrents[vs.name] ||= []).push(i)
+    }
+
+    for (const c of ckt.C) {
+      const v1 = c.n1 === 0 ? 0 : (x[c.n1 - 1] ?? 0)
+      const v2 = c.n2 === 0 ? 0 : (x[c.n2 - 1] ?? 0)
+      c.vPrev = v1 - v2
+    }
+    for (const l of ckt.L) {
+      const v1 = l.n1 === 0 ? 0 : (x[l.n1 - 1] ?? 0)
+      const v2 = l.n2 === 0 ? 0 : (x[l.n2 - 1] ?? 0)
+      const Gl = Math.max(dt, EPS) / l.L
+      l.iPrev = Gl * (v1 - v2) + l.iPrev
+    }
+  }
+
+  return { times, nodeVoltages, elementCurrents }
+}
+
+export { simulateTRAN }
diff --git a/lib/constants/EPS.ts b/lib/constants/EPS.ts
new file mode 100644
index 0000000..a9a2c70
--- /dev/null
+++ b/lib/constants/EPS.ts
@@ -0,0 +1 @@
+export const EPS = 1e-15
diff --git a/lib/formatting/formatAcResult.ts b/lib/formatting/formatAcResult.ts
new file mode 100644
index 0000000..5565ba2
--- /dev/null
+++ b/lib/formatting/formatAcResult.ts
@@ -0,0 +1,27 @@
+import { Complex } from "../math/Complex"
+
+function formatAcResult(
+  ac: {
+    freqs: number[]
+    nodeVoltages: Record<string, Complex[]>
+  } | null,
+) {
+  if (!ac) return "No AC analysis.\n"
+  const nodes = Object.keys(ac.nodeVoltages)
+  const lines: string[] = []
+  lines.push(`f(Hz), ` + nodes.map((n) => `${n}:|V|,∠V(deg)`).join(", "))
+  for (let k = 0; k < ac.freqs.length; k++) {
+    const freq = ac.freqs[k]
+    if (freq == null) continue
+    const parts = [freq.toPrecision(6)]
+    for (const n of nodes) {
+      const z = ac.nodeVoltages[n]?.[k]
+      if (!z) continue
+      parts.push(`${z.abs().toPrecision(6)},${z.phaseDeg().toPrecision(6)}`)
+    }
+    lines.push(parts.join(", "))
+  }
+  return lines.join("\n")
+}
+
+export { formatAcResult }
diff --git a/lib/formatting/formatTranResult.ts b/lib/formatting/formatTranResult.ts
new file mode 100644
index 0000000..957f67b
--- /dev/null
+++ b/lib/formatting/formatTranResult.ts
@@ -0,0 +1,25 @@
+function formatTranResult(
+  tran: {
+    times: number[]
+    nodeVoltages: Record<string, number[]>
+  } | null,
+) {
+  if (!tran) return "No TRAN analysis.\n"
+  const nodes = Object.keys(tran.nodeVoltages)
+  const header = ["t(s)", ...nodes.map((n) => `${n}:V`)]
+  const lines = [header.join(", ")]
+  for (let k = 0; k < tran.times.length; k++) {
+    const time = tran.times[k]
+    if (time == null) continue
+    const row = [time.toPrecision(6)]
+    for (const n of nodes) {
+      const value = tran.nodeVoltages[n]?.[k]
+      if (value == null) continue
+      row.push((+value).toPrecision(6))
+    }
+    lines.push(row.join(", "))
+  }
+  return lines.join("\n")
+}
+
+export { formatTranResult }
diff --git a/lib/index.ts b/lib/index.ts
index 45f812f..15350ef 100644
--- a/lib/index.ts
+++ b/lib/index.ts
@@ -1,814 +1,7 @@
-/*------------------------------------------------------------------------------
-  Simple SPICE in native JavaScript
-  - AC sweep (.ac dec/lin) for R, L, C, independent V sources
-  - Transient (.tran) with backward-Euler for R, L, C, independent V sources
-  - PULSE() and DC for time-domain independent V sources; AC <mag> [phase] for AC
-  - Minimal, educational, not a full SPICE (no nonlinear devices, no .op, etc.)
-------------------------------------------------------------------------------*/
-
-// ------------------------ Utilities ------------------------
-
-const EPS = 1e-15
-
-const unitMul = {
-  // SPICE-style suffixes (case-insensitive)
-  t: 1e12,
-  g: 1e9,
-  meg: 1e6,
-  k: 1e3,
-  m: 1e-3,
-  u: 1e-6,
-  n: 1e-9,
-  p: 1e-12,
-  f: 1e-15,
-}
-
-function parseNumberWithUnits(raw) {
-  if (raw == null) return NaN
-  let s = String(raw).trim()
-  if (s === "") return NaN
-
-  // If it's a plain number (possibly scientific notation), parse directly.
-  if (/^[+-]?\d*\.?\d+(?:[eE][+-]?\d+)?$/.test(s)) return parseFloat(s)
-
-  // E.g., "100u", "10k", "1.8V", "0.09u", "1meg", "100uF"
-  const m = s.match(/^([+-]?\d*\.?\d+(?:[eE][+-]?\d+)?)([a-zA-Z]+)$/)
-  if (!m) return parseFloat(s) // last-ditch attempt; may produce NaN
-
-  let val = parseFloat(m[1])
-  let suf = m[2].toLowerCase()
-
-  // Strip dimension letters like 'f','h','v','a','s','ohm' if appended to suffix
-  // so "100uF" -> "u", "1kOhm" -> "k".
-  suf = suf.replace(/(ohm|v|a|s|h|f)$/g, "")
-
-  if (suf === "meg") return val * unitMul.meg
-  if (suf.length === 1 && unitMul.hasOwnProperty(suf)) return val * unitMul[suf]
-
-  // Unknown suffix -> just parse number part.
-  return val
-}
-
-// Tokenizer that keeps parenthesis groups intact, e.g. PULSE(...) as one token
-function smartTokens(line) {
-  const re = /"[^"]*"|\([^()]*\)|\S+/g
-  const out = []
-  let m
-  while ((m = re.exec(line)) !== null) out.push(m[0])
-  return out
-}
-
-// Log-spaced frequencies for .ac dec
-function logspace(f1, f2, pointsPerDec) {
-  if (f1 <= 0 || f2 <= 0) throw new Error(".ac frequencies must be > 0")
-  if (f2 < f1) [f1, f2] = [f2, f1]
-  const decades = Math.log10(f2 / f1)
-  const n = Math.max(1, Math.ceil(decades * pointsPerDec))
-  const arr = []
-  for (let i = 0; i <= n; i++) {
-    arr.push(f1 * Math.pow(10, i / pointsPerDec))
-  }
-  // Ensure last point ≤ f2 (and include f2 if needed)
-  if (arr[arr.length - 1] < f2 * (1 - 1e-12)) arr.push(f2)
-  return arr
-}
-
-// ------------------------ Complex numbers ------------------------
-
-class Complex {
-  constructor(re = 0, im = 0) {
-    this.re = re
-    this.im = im
-  }
-  static from(re, im = 0) {
-    return new Complex(re, im)
-  }
-  static fromPolar(mag, deg = 0) {
-    const ph = (deg * Math.PI) / 180
-    return new Complex(mag * Math.cos(ph), mag * Math.sin(ph))
-  }
-  clone() {
-    return new Complex(this.re, this.im)
-  }
-  add(b) {
-    return new Complex(this.re + b.re, this.im + b.im)
-  }
-  sub(b) {
-    return new Complex(this.re - b.re, this.im - b.im)
-  }
-  mul(b) {
-    return new Complex(
-      this.re * b.re - this.im * b.im,
-      this.re * b.im + this.im * b.re,
-    )
-  }
-  div(b) {
-    const d = b.re * b.re + b.im * b.im
-    if (d < EPS) throw new Error("Complex divide by ~0")
-    return new Complex(
-      (this.re * b.re + this.im * b.im) / d,
-      (this.im * b.re - this.re * b.im) / d,
-    )
-  }
-  inv() {
-    const d = this.re * this.re + this.im * this.im
-    if (d < EPS) throw new Error("Complex invert by ~0")
-    return new Complex(this.re / d, -this.im / d)
-  }
-  abs() {
-    return Math.hypot(this.re, this.im)
-  }
-  phaseDeg() {
-    return (Math.atan2(this.im, this.re) * 180) / Math.PI
-  }
-}
-
-// ------------------------ Linear solvers ------------------------
-
-function solveReal(A, b) {
-  // Gaussian elimination with partial pivoting
-  const n = A.length
-  // Build augmented matrix
-  for (let i = 0; i < n; i++) {
-    A[i] = A[i].slice() // clone row
-    A[i].push(b[i])
-  }
-
-  for (let k = 0; k < n; k++) {
-    // Pivot
-    let imax = k
-    let vmax = Math.abs(A[k][k])
-    for (let i = k + 1; i < n; i++) {
-      const v = Math.abs(A[i][k])
-      if (v > vmax) {
-        vmax = v
-        imax = i
-      }
-    }
-    if (vmax < EPS) throw new Error("Singular matrix (real)")
-
-    if (imax !== k) {
-      const tmp = A[k]
-      A[k] = A[imax]
-      A[imax] = tmp
-    }
-
-    // Eliminate
-    const pivot = A[k][k]
-    for (let i = k + 1; i < n; i++) {
-      const f = A[i][k] / pivot
-      if (Math.abs(f) < EPS) continue
-      for (let j = k; j <= n; j++) A[i][j] -= f * A[k][j]
-    }
-  }
-
-  // Back substitution
-  const x = new Array(n).fill(0)
-  for (let i = n - 1; i >= 0; i--) {
-    let s = A[i][n]
-    for (let j = i + 1; j < n; j++) s -= A[i][j] * x[j]
-    x[i] = s / A[i][i]
-  }
-  return x
-}
-
-function solveComplex(A, b) {
-  const n = A.length
-  // Augment
-  for (let i = 0; i < n; i++) {
-    A[i] = A[i].map((z) => (z.clone ? z.clone() : Complex.from(z, 0)))
-    A[i].push(b[i].clone ? b[i].clone() : Complex.from(b[i], 0))
-  }
-  // GE with partial pivoting on |pivot|
-  for (let k = 0; k < n; k++) {
-    let imax = k,
-      vmax = A[k][k].abs()
-    for (let i = k + 1; i < n; i++) {
-      const v = A[i][k].abs()
-      if (v > vmax) {
-        vmax = v
-        imax = i
-      }
-    }
-    if (vmax < EPS) throw new Error("Singular matrix (complex)")
-    if (imax !== k) {
-      const tmp = A[k]
-      A[k] = A[imax]
-      A[imax] = tmp
-    }
-
-    const pivot = A[k][k]
-    for (let i = k + 1; i < n; i++) {
-      const f = A[i][k].div(pivot)
-      if (f.abs() < EPS) continue
-      for (let j = k; j <= n; j++) A[i][j] = A[i][j].sub(f.mul(A[k][j]))
-    }
-  }
-  // Back substitute
-  const x = new Array(n)
-  for (let i = n - 1; i >= 0; i--) {
-    let s = A[i][n]
-    for (let j = i + 1; j < n; j++) s = s.sub(A[i][j].mul(x[j]))
-    x[i] = s.div(A[i][i])
-  }
-  return x
-}
-
-// ------------------------ Circuit & parsing ------------------------
-
-class NodeIndex {
-  constructor() {
-    this.map = new Map([["0", 0]]) // ground
-    this.rev = ["0"]
-  }
-  getOrCreate(name) {
-    const key = String(name)
-    if (this.map.has(key)) return this.map.get(key)
-    const idx = this.rev.length
-    this.map.set(key, idx)
-    this.rev.push(key)
-    return idx
-  }
-  get(name) {
-    return this.map.get(String(name))
-  }
-  count() {
-    return this.rev.length
-  }
-  // For MNA unknowns, ground (0) is removed; nodes 1..N-1 map to indices 0..N-2
-  matrixIndexOfNode(nodeId) {
-    if (nodeId === 0) return -1
-    return nodeId - 1
-  }
-}
-
-// Waveforms (time-domain)
-function parsePulseArgs(s) {
-  // PULSE (V1 V2 TD TR TF TON PERIOD [Ncycles])
-  const clean = s.trim().replace(/^pulse\s*\(/i, "(")
-  const inside = clean.replace(/^\(/, "").replace(/\)$/, "").trim()
-  const parts = inside.split(/[\s,]+/).filter((x) => x.length)
-  if (parts.length < 7) throw new Error("PULSE(...) requires 7 or 8 args")
-  const vals = parts.map(parseNumberWithUnits)
-  return {
-    v1: vals[0],
-    v2: vals[1],
-    td: vals[2],
-    tr: vals[3],
-    tf: vals[4],
-    ton: vals[5],
-    period: vals[6],
-    ncycles: parts[7] != null ? vals[7] : Infinity,
-  }
-}
-
-function pulseValue(p, t) {
-  if (t < p.td) return p.v1
-  const tt = t - p.td
-  const cyclesDone = Math.floor(tt / p.period)
-  if (cyclesDone >= p.ncycles) return p.v1
-  const tc = tt - cyclesDone * p.period
-
-  if (tc < p.tr) {
-    // rise
-    const a = tc / Math.max(p.tr, EPS)
-    return p.v1 + (p.v2 - p.v1) * a
-  } else if (tc < p.tr + p.ton) {
-    // on
-    return p.v2
-  } else if (tc < p.tr + p.ton + p.tf) {
-    // fall
-    const a = (tc - (p.tr + p.ton)) / Math.max(p.tf, EPS)
-    return p.v2 + (p.v1 - p.v2) * a
-  } else {
-    // off
-    return p.v1
-  }
-}
-
-class Circuit {
-  constructor() {
-    this.nodes = new NodeIndex()
-    this.R = [] // {name,n1,n2,R}
-    this.C = [] // {name,n1,n2,C,vPrev}
-    this.L = [] // {name,n1,n2,L,iPrev}
-    this.V = [] // {name,n1,n2,dc,acMag,acPhaseDeg,waveform, index}
-    this.analyses = { ac: null, tran: null }
-    this.skipped = [] // lines/devices not supported
-  }
-  nodeId(n) {
-    return this.nodes.getOrCreate(n)
-  }
-  nodeNameById(id) {
-    return this.nodes.rev[id]
-  }
-
-  addRes(name, n1, n2, val) {
-    this.R.push({ name, n1, n2, R: val })
-  }
-  addCap(name, n1, n2, val) {
-    this.C.push({ name, n1, n2, C: val, vPrev: 0 })
-  }
-  addInd(name, n1, n2, val) {
-    this.L.push({ name, n1, n2, L: val, iPrev: 0 })
-  }
-  addVsrc(name, n1, n2, spec) {
-    const vs = Object.assign(
-      {
-        name,
-        n1,
-        n2,
-        dc: 0,
-        acMag: 0,
-        acPhaseDeg: 0,
-        waveform: null,
-        index: -1,
-      },
-      spec || {},
-    )
-    this.V.push(vs)
-  }
-}
-
-function parseNetlist(text) {
-  const ckt = new Circuit()
-  const lines = text.split(/\r?\n/)
-  let seenTitle = false
-
-  for (let raw of lines) {
-    let line = raw.trim()
-    if (!line) continue
-
-    // Kill full-line '*' comments
-    if (/^\*/.test(line)) continue
-
-    // Stop at .end (still allow trailing whitespace)
-    if (/^\s*\.end\b/i.test(line)) break
-
-    // Strip inline comments after '//' or ';' (not SPICE standard, but convenient)
-    line = line.replace(/\/\/.*$/, "")
-    line = line.replace(/;.*$/, "")
-
-    const tokens = smartTokens(line)
-    if (tokens.length === 0) continue
-
-    const first = tokens[0]
-
-    // Title line: allowed to be any free text if first non-empty line and not element/directive
-    if (!seenTitle && !/^[rclvgmiqd]\w*$/i.test(first) && !/^\./.test(first)) {
-      seenTitle = true
-      continue
-    }
-
-    // Directives
-    if (/^\./.test(first)) {
-      const dir = first.toLowerCase()
-      if (dir === ".ac") {
-        // .ac {dec|lin} N fstart fstop
-        const mode = tokens[1]?.toLowerCase()
-        if (mode !== "dec" && mode !== "lin")
-          throw new Error(".ac supports 'dec' or 'lin'")
-        const N = parseInt(tokens[2], 10)
-        const f1 = parseNumberWithUnits(tokens[3])
-        const f2 = parseNumberWithUnits(tokens[4])
-        ckt.analyses.ac = { mode, N, f1, f2 }
-      } else if (dir === ".tran") {
-        // .tran tstep tstop [tstart ...] -> we use first two
-        const dt = parseNumberWithUnits(tokens[1])
-        const tstop = parseNumberWithUnits(tokens[2])
-        ckt.analyses.tran = { dt, tstop }
-      } else if (
-        dir === ".include" ||
-        dir === ".lib" ||
-        dir === ".model" ||
-        dir === ".options" ||
-        dir === ".op"
-      ) {
-        ckt.skipped.push(line)
-      } else {
-        ckt.skipped.push(line)
-      }
-      continue
-    }
-
-    // Elements: R, C, L, V (independent V)
-    const typeChar = first[0].toLowerCase()
-    const name = first
-
-    try {
-      if (typeChar === "r") {
-        // Rname n1 n2 value
-        const n1 = ckt.nodeId(tokens[1])
-        const n2 = ckt.nodeId(tokens[2])
-        const val = parseNumberWithUnits(tokens[3])
-        ckt.addRes(name, n1, n2, val)
-      } else if (typeChar === "c") {
-        const n1 = ckt.nodeId(tokens[1])
-        const n2 = ckt.nodeId(tokens[2])
-        const val = parseNumberWithUnits(tokens[3])
-        ckt.addCap(name, n1, n2, val)
-      } else if (typeChar === "l") {
-        const n1 = ckt.nodeId(tokens[1])
-        const n2 = ckt.nodeId(tokens[2])
-        const val = parseNumberWithUnits(tokens[3])
-        ckt.addInd(name, n1, n2, val)
-      } else if (typeChar === "v") {
-        // Vname n+ n- [<DC value>] [dc <val>] [ac <mag> [phase]] [pulse (...)]
-        const n1 = ckt.nodeId(tokens[1])
-        const n2 = ckt.nodeId(tokens[2])
-        let spec = { dc: 0, acMag: 0, acPhaseDeg: 0, waveform: null }
-
-        // Flexible scan of the remainder
-        let i = 3
-        if (i < tokens.length && !/^[a-zA-Z]/.test(tokens[i])) {
-          // Bare value = DC
-          spec.dc = parseNumberWithUnits(tokens[i])
-          i++
-        }
-        while (i < tokens.length) {
-          const key = tokens[i].toLowerCase()
-          if (key === "dc") {
-            spec.dc = parseNumberWithUnits(tokens[i + 1])
-            i += 2
-          } else if (key === "ac") {
-            spec.acMag = parseNumberWithUnits(tokens[i + 1])
-            if (i + 2 < tokens.length && /^[+-]?\d/.test(tokens[i + 2])) {
-              spec.acPhaseDeg = parseNumberWithUnits(tokens[i + 2])
-              i += 3
-            } else {
-              i += 2
-            }
-          } else if (key.startsWith("pulse")) {
-            // Could be "pulse(...)" in one token or "pulse" + "(...)"
-            let argToken = key.includes("(") ? key : tokens[i + 1]
-            if (!argToken || !/\(.*\)/.test(argToken)) {
-              throw new Error("Malformed PULSE() specification")
-            }
-            const p = parsePulseArgs(argToken)
-            spec.waveform = (t) => pulseValue(p, t)
-            i += key.includes("(") ? 1 : 2
-          } else if (/^\(.*\)$/.test(key)) {
-            // Lone "(...)" immediately after pulse
-            i++
-          } else {
-            // Unrecognized -> stop scanning (to be permissive)
-            i++
-          }
-        }
-        ckt.addVsrc(name, n1, n2, spec)
-      } else {
-        // Unknown device (M, D, I, etc.)
-        ckt.skipped.push(line)
-      }
-    } catch (e) {
-      throw new Error(`Parse error on line: "${line}"\n${e.message}`)
-    }
-  }
-
-  // Assign MNA indices to voltage sources (branch currents)
-  const nNodes = ckt.nodes.count() - 1 // excluding ground
-  for (let i = 0; i < ckt.V.length; i++) {
-    ckt.V[i].index = nNodes + i
-  }
-  return ckt
-}
-
-// ------------------------ Stamping helpers ------------------------
-
-// Add admittance Y (real) between n1-n2 (node ids) into real matrix
-function stampY_real(A, nidx, n1, n2, Y) {
-  const i1 = nidx.matrixIndexOfNode(n1)
-  const i2 = nidx.matrixIndexOfNode(n2)
-  if (i1 >= 0) A[i1][i1] += Y
-  if (i2 >= 0) A[i2][i2] += Y
-  if (i1 >= 0 && i2 >= 0) {
-    A[i1][i2] -= Y
-    A[i2][i1] -= Y
-  }
-}
-
-// Add admittance Y (complex) between n1-n2 into complex matrix
-function stampY_cplx(A, nidx, n1, n2, Y) {
-  const i1 = nidx.matrixIndexOfNode(n1)
-  const i2 = nidx.matrixIndexOfNode(n2)
-  if (i1 >= 0) A[i1][i1] = A[i1][i1].add(Y)
-  if (i2 >= 0) A[i2][i2] = A[i2][i2].add(Y)
-  if (i1 >= 0 && i2 >= 0) {
-    A[i1][i2] = A[i1][i2].sub(Y)
-    A[i2][i1] = A[i2][i1].sub(Y)
-  }
-}
-
-// Stamp current source I from nPlus to nMinus into real RHS
-function stampI_real(b, nidx, nPlus, nMinus, I) {
-  const iPlus = nidx.matrixIndexOfNode(nPlus)
-  const iMinus = nidx.matrixIndexOfNode(nMinus)
-  if (iPlus >= 0) b[iPlus] -= I
-  if (iMinus >= 0) b[iMinus] += I
-}
-
-// Stamp current source (complex) into complex RHS
-function stampI_cplx(b, nidx, nPlus, nMinus, I) {
-  const iPlus = nidx.matrixIndexOfNode(nPlus)
-  const iMinus = nidx.matrixIndexOfNode(nMinus)
-  if (iPlus >= 0) b[iPlus] = b[iPlus].sub(I)
-  if (iMinus >= 0) b[iMinus] = b[iMinus].add(I)
-}
-
-// Stamp independent voltage source (MNA) into real A,b, with value V (real)
-function stampVsrc_real(A, b, nidx, vs, V, totalUnknowns) {
-  const i1 = nidx.matrixIndexOfNode(vs.n1)
-  const i2 = nidx.matrixIndexOfNode(vs.n2)
-  const j = vs.index // global unknown index (includes node vars)
-  // Ensure matrix has room
-  // KCL rows:
-  if (i1 >= 0) A[i1][j] += 1
-  if (i2 >= 0) A[i2][j] -= 1
-  // Source equation row
-  A[j][j] += 0 // keep for clarity
-  if (i1 >= 0) A[j][i1] += 1
-  if (i2 >= 0) A[j][i2] -= 1
-  b[j] += V
-}
-
-// Stamp independent voltage source into complex A,b, with value V (Complex)
-function stampVsrc_cplx(A, b, nidx, vs, V) {
-  const i1 = nidx.matrixIndexOfNode(vs.n1)
-  const i2 = nidx.matrixIndexOfNode(vs.n2)
-  const j = vs.index
-  if (i1 >= 0) A[i1][j] = A[i1][j].add(Complex.from(1, 0))
-  if (i2 >= 0) A[i2][j] = A[i2][j].sub(Complex.from(1, 0))
-  if (i1 >= 0) A[j][i1] = A[j][i1].add(Complex.from(1, 0))
-  if (i2 >= 0) A[j][i2] = A[j][i2].sub(Complex.from(1, 0))
-  b[j] = b[j].add(V)
-}
-
-// ------------------------ Analyses ------------------------
-
-function simulateAC(ckt) {
-  if (!ckt.analyses.ac) return null
-
-  const { mode, N, f1, f2 } = ckt.analyses.ac
-  const nNodeVars = ckt.nodes.count() - 1 // exclude ground
-  const nVsrc = ckt.V.length
-  const Nvar = nNodeVars + nVsrc
-
-  const freqs =
-    mode === "dec"
-      ? logspace(f1, f2, N)
-      : (function () {
-          const arr = []
-          const npts = Math.max(2, N)
-          const step = (f2 - f1) / (npts - 1)
-          for (let i = 0; i < npts; i++) arr.push(f1 + i * step)
-          return arr
-        })()
-
-  // Results: nodeName -> array of Complex; element currents map likewise
-  const Vout = {}
-  ckt.nodes.rev.forEach((name, id) => {
-    if (id !== 0) Vout[name] = []
-  })
-  const Ielts = {} // element currents (Res, Cap, Ind, Vsrc) by name
-
-  const twoPi = 2 * Math.PI
-
-  for (const f of freqs) {
-    // Build complex matrix A (Nvar x Nvar) and RHS b (Nvar)
-    const A = Array.from({ length: Nvar }, () =>
-      Array.from({ length: Nvar }, () => Complex.from(0, 0)),
-    )
-    const b = Array.from({ length: Nvar }, () => Complex.from(0, 0))
-
-    // R stamp: Y = 1/R
-    for (const r of ckt.R) {
-      if (r.R <= 0) throw new Error(`R ${r.name} must be > 0`)
-      const Y = Complex.from(1 / r.R, 0)
-      stampY_cplx(A, ckt.nodes, r.n1, r.n2, Y)
-    }
-
-    // C stamp: Y = j*omega*C
-    for (const c of ckt.C) {
-      const Y = Complex.from(0, twoPi * f * c.C) // jωC
-      stampY_cplx(A, ckt.nodes, c.n1, c.n2, Y)
-    }
-
-    // L stamp: Y = 1/(j*omega*L)
-    for (const l of ckt.L) {
-      const denom = Complex.from(0, twoPi * f * l.L) // jωL
-      const Y =
-        denom.abs() < EPS ? Complex.from(0, 0) : Complex.from(1, 0).div(denom) // 1/(jωL)
-      stampY_cplx(A, ckt.nodes, l.n1, l.n2, Y)
-    }
-
-    // V sources: include ALL voltage sources in AC. DC-only become zero amplitude (=short in small-signal).
-    for (const vs of ckt.V) {
-      const Vph = Complex.fromPolar(vs.acMag || 0, vs.acPhaseDeg || 0)
-      stampVsrc_cplx(A, b, ckt.nodes, vs, Vph)
-    }
-
-    // Solve
-    const x = solveComplex(A, b)
-
-    // Extract node voltages
-    for (let id = 1; id < ckt.nodes.count(); id++) {
-      const idx = id - 1
-      Vout[ckt.nodeNameById(id)].push(x[idx])
-    }
-
-    // Element currents (phasors)
-    // R, C, L via I = Y * (v1 - v2); Vsrc via branch current x[j]
-    for (const r of ckt.R) {
-      const v1 = r.n1 === 0 ? Complex.from(0, 0) : x[r.n1 - 1]
-      const v2 = r.n2 === 0 ? Complex.from(0, 0) : x[r.n2 - 1]
-      const Y = Complex.from(1 / r.R, 0)
-      const i = Y.mul(v1.sub(v2)) // current from n1 -> n2
-      ;(Ielts[r.name] ||= []).push(i)
-    }
-    for (const c of ckt.C) {
-      const v1 = c.n1 === 0 ? Complex.from(0, 0) : x[c.n1 - 1]
-      const v2 = c.n2 === 0 ? Complex.from(0, 0) : x[c.n2 - 1]
-      const Y = Complex.from(0, twoPi * f * c.C)
-      const i = Y.mul(v1.sub(v2))
-      ;(Ielts[c.name] ||= []).push(i)
-    }
-    for (const l of ckt.L) {
-      const v1 = l.n1 === 0 ? Complex.from(0, 0) : x[l.n1 - 1]
-      const v2 = l.n2 === 0 ? Complex.from(0, 0) : x[l.n2 - 1]
-      const Y = (function () {
-        const d = Complex.from(0, twoPi * f * l.L)
-        return d.abs() < EPS ? Complex.from(0, 0) : Complex.from(1, 0).div(d)
-      })()
-      const i = Y.mul(v1.sub(v2))
-      ;(Ielts[l.name] ||= []).push(i)
-    }
-    for (const vs of ckt.V) {
-      const i = x[vs.index] // branch current from n1->n2
-      ;(Ielts[vs.name] ||= []).push(i)
-    }
-  }
-
-  return { freqs, nodeVoltages: Vout, elementCurrents: Ielts }
-}
-
-function simulateTRAN(ckt) {
-  if (!ckt.analyses.tran) return null
-  const { dt, tstop } = ckt.analyses.tran
-  const steps = Math.max(1, Math.ceil(tstop / Math.max(dt, EPS)))
-
-  const nNodeVars = ckt.nodes.count() - 1
-  const nVsrc = ckt.V.length
-  const Nvar = nNodeVars + nVsrc
-
-  const times = []
-  const Vout = {}
-  ckt.nodes.rev.forEach((name, id) => {
-    if (id !== 0) Vout[name] = []
-  })
-  const Ielts = {} // element currents by name vs time
-
-  // Pre-allocate matrices (real)
-  function blankA() {
-    const A = Array.from({ length: Nvar }, () => new Array(Nvar).fill(0))
-    return A
-  }
-  function blankB() {
-    return new Array(Nvar).fill(0)
-  }
-
-  // Time stepping (Backward Euler companion models)
-  let t = 0
-  for (let step = 0; step <= steps; step++, t = step * dt) {
-    times.push(t)
-    const A = blankA()
-    const b = blankB()
-
-    // Resistors
-    for (const r of ckt.R) {
-      const G = 1 / r.R
-      stampY_real(A, ckt.nodes, r.n1, r.n2, G)
-    }
-
-    // Capacitors: Norton Gc = C/dt, Ieq = -Gc*vPrev (from n1->n2)
-    for (const c of ckt.C) {
-      const Gc = c.C / Math.max(dt, EPS)
-      stampY_real(A, ckt.nodes, c.n1, c.n2, Gc)
-      const Ieq = -Gc * c.vPrev
-      stampI_real(b, ckt.nodes, c.n1, c.n2, Ieq)
-    }
-
-    // Inductors: Norton Gl = dt/L, Ieq = iPrev (from n1->n2)
-    for (const l of ckt.L) {
-      const Gl = Math.max(dt, EPS) / l.L // conductance
-      stampY_real(A, ckt.nodes, l.n1, l.n2, Gl)
-      stampI_real(b, ckt.nodes, l.n1, l.n2, l.iPrev)
-    }
-
-    // Voltage sources (independent): set branch equation and KCL couplings
-    for (const vs of ckt.V) {
-      // Instantaneous value
-      const Vt = vs.waveform ? vs.waveform(t) : vs.dc || 0
-      stampVsrc_real(A, b, ckt.nodes, vs, Vt, Nvar)
-    }
-
-    // Solve
-    const x = solveReal(A, b)
-
-    // Store node voltages
-    for (let id = 1; id < ckt.nodes.count(); id++) {
-      const idx = id - 1
-      Vout[ckt.nodeNameById(id)].push(x[idx])
-    }
-
-    // Element currents (from n1 -> n2)
-    for (const r of ckt.R) {
-      const v1 = r.n1 === 0 ? 0 : x[r.n1 - 1]
-      const v2 = r.n2 === 0 ? 0 : x[r.n2 - 1]
-      const i = (v1 - v2) / r.R
-      ;(Ielts[r.name] ||= []).push(i)
-    }
-    for (const c of ckt.C) {
-      const v1 = c.n1 === 0 ? 0 : x[c.n1 - 1]
-      const v2 = c.n2 === 0 ? 0 : x[c.n2 - 1]
-      // Backward-Euler i = C*(v - vPrev)/dt
-      const i = (c.C * (v1 - v2 - c.vPrev)) / Math.max(dt, EPS)
-      ;(Ielts[c.name] ||= []).push(i)
-    }
-    for (const l of ckt.L) {
-      const v1 = l.n1 === 0 ? 0 : x[l.n1 - 1]
-      const v2 = l.n2 === 0 ? 0 : x[l.n2 - 1]
-      const Gl = Math.max(dt, EPS) / l.L
-      const i = Gl * (v1 - v2) + l.iPrev // Norton: i = Gl*v + iPrev
-      ;(Ielts[l.name] ||= []).push(i)
-    }
-    for (const vs of ckt.V) {
-      const i = x[vs.index] // branch current, orientation n1->n2
-      ;(Ielts[vs.name] ||= []).push(i)
-    }
-
-    // Update dynamic states for next step
-    for (const c of ckt.C) {
-      const v1 = c.n1 === 0 ? 0 : x[c.n1 - 1]
-      const v2 = c.n2 === 0 ? 0 : x[c.n2 - 1]
-      c.vPrev = v1 - v2
-    }
-    for (const l of ckt.L) {
-      const v1 = l.n1 === 0 ? 0 : x[l.n1 - 1]
-      const v2 = l.n2 === 0 ? 0 : x[l.n2 - 1]
-      const Gl = Math.max(dt, EPS) / l.L
-      l.iPrev = Gl * (v1 - v2) + l.iPrev // update using computed i
-    }
-  }
-
-  return { times, nodeVoltages: Vout, elementCurrents: Ielts }
-}
-
-// ------------------------ Facade ------------------------
-
-function simulate(netlistText) {
-  const ckt = parseNetlist(netlistText)
-  const ac = simulateAC(ckt)
-  const tran = simulateTRAN(ckt)
-  return { circuit: ckt, ac, tran }
-}
-
-// Pretty format helpers for quick inspection
-function formatAcResult(ac) {
-  if (!ac) return "No AC analysis.\n"
-  const nodes = Object.keys(ac.nodeVoltages)
-  const lines = []
-  lines.push(`f(Hz), ` + nodes.map((n) => `${n}:|V|,∠V(deg)`).join(", "))
-  for (let k = 0; k < ac.freqs.length; k++) {
-    const parts = [ac.freqs[k].toPrecision(6)]
-    for (const n of nodes) {
-      const z = ac.nodeVoltages[n][k]
-      parts.push(`${z.abs().toPrecision(6)},${z.phaseDeg().toPrecision(6)}`)
-    }
-    lines.push(parts.join(", "))
-  }
-  return lines.join("\n")
-}
-
-function formatTranResult(tran) {
-  if (!tran) return "No TRAN analysis.\n"
-  const nodes = Object.keys(tran.nodeVoltages)
-  const header = ["t(s)", ...nodes.map((n) => `${n}:V`)]
-  const lines = [header.join(", ")]
-  for (let k = 0; k < tran.times.length; k++) {
-    const row = [tran.times[k].toPrecision(6)]
-    for (const n of nodes) row.push((+tran.nodeVoltages[n][k]).toPrecision(6))
-    lines.push(row.join(", "))
-  }
-  return lines.join("\n")
-}
-
-export {
-  parseNetlist,
-  simulate,
-  simulateAC,
-  simulateTRAN,
-  formatAcResult,
-  formatTranResult,
-  // Expose classes if you want to extend
-  Complex,
-}
+export { parseNetlist } from "./parsing/parseNetlist"
+export { simulate } from "./analysis/simulate"
+export { simulateAC } from "./analysis/simulateAC"
+export { simulateTRAN } from "./analysis/simulateTRAN"
+export { formatAcResult } from "./formatting/formatAcResult"
+export { formatTranResult } from "./formatting/formatTranResult"
+export { Complex } from "./math/Complex"
diff --git a/lib/math/Complex.ts b/lib/math/Complex.ts
new file mode 100644
index 0000000..508b79d
--- /dev/null
+++ b/lib/math/Complex.ts
@@ -0,0 +1,64 @@
+import { EPS } from "../constants/EPS"
+
+class Complex {
+  re: number
+  im: number
+
+  constructor(re = 0, im = 0) {
+    this.re = re
+    this.im = im
+  }
+
+  static from(re: number, im = 0) {
+    return new Complex(re, im)
+  }
+
+  static fromPolar(mag: number, deg = 0) {
+    const ph = (deg * Math.PI) / 180
+    return new Complex(mag * Math.cos(ph), mag * Math.sin(ph))
+  }
+
+  clone() {
+    return new Complex(this.re, this.im)
+  }
+
+  add(b: Complex) {
+    return new Complex(this.re + b.re, this.im + b.im)
+  }
+
+  sub(b: Complex) {
+    return new Complex(this.re - b.re, this.im - b.im)
+  }
+
+  mul(b: Complex) {
+    return new Complex(
+      this.re * b.re - this.im * b.im,
+      this.re * b.im + this.im * b.re,
+    )
+  }
+
+  div(b: Complex) {
+    const d = b.re * b.re + b.im * b.im
+    if (d < EPS) throw new Error("Complex divide by ~0")
+    return new Complex(
+      (this.re * b.re + this.im * b.im) / d,
+      (this.im * b.re - this.re * b.im) / d,
+    )
+  }
+
+  inv() {
+    const d = this.re * this.re + this.im * this.im
+    if (d < EPS) throw new Error("Complex invert by ~0")
+    return new Complex(this.re / d, -this.im / d)
+  }
+
+  abs() {
+    return Math.hypot(this.re, this.im)
+  }
+
+  phaseDeg() {
+    return (Math.atan2(this.im, this.re) * 180) / Math.PI
+  }
+}
+
+export { Complex }
diff --git a/lib/math/solveComplex.ts b/lib/math/solveComplex.ts
new file mode 100644
index 0000000..0a12cd9
--- /dev/null
+++ b/lib/math/solveComplex.ts
@@ -0,0 +1,75 @@
+import { EPS } from "../constants/EPS"
+import { Complex } from "./Complex"
+
+function solveComplex(A: Complex[][], b: Complex[]) {
+  const n = A.length
+  for (let i = 0; i < n; i++) {
+    const row = A[i]
+    const bi = b[i]
+    if (!row || !bi) throw new Error("Matrix dimensions mismatch")
+    const copy = row.map((z) => z.clone())
+    copy.push(bi.clone())
+    A[i] = copy
+  }
+
+  for (let k = 0; k < n; k++) {
+    let imax = k
+    const pivotRow = A[k]
+    if (!pivotRow) throw new Error("Matrix row missing")
+    let vmax = pivotRow[k]?.abs() ?? 0
+    for (let i = k + 1; i < n; i++) {
+      const row = A[i]
+      if (!row) throw new Error("Matrix row missing")
+      const v = row[k]?.abs() ?? 0
+      if (v > vmax) {
+        vmax = v
+        imax = i
+      }
+    }
+    if (vmax < EPS) throw new Error("Singular matrix (complex)")
+    if (imax !== k) {
+      const tmp = A[k]
+      A[k] = A[imax]!
+      A[imax] = tmp!
+    }
+
+    const pivotRowUpdated = A[k]
+    if (!pivotRowUpdated) throw new Error("Pivot row missing")
+    const pivot = pivotRowUpdated[k]
+    if (!pivot) throw new Error("Zero pivot encountered")
+    for (let i = k + 1; i < n; i++) {
+      const row = A[i]
+      if (!row) throw new Error("Matrix row missing")
+      const entry = row[k]
+      if (!entry) continue
+      const f = entry.div(pivot)
+      if (f.abs() < EPS) continue
+      for (let j = k; j <= n; j++) {
+        const target = row[j]
+        const source = pivotRowUpdated[j]
+        if (!target || !source) continue
+        row[j] = target.sub(f.mul(source))
+      }
+    }
+  }
+
+  const x = new Array(n)
+  for (let i = n - 1; i >= 0; i--) {
+    const row = A[i]
+    if (!row) throw new Error("Matrix row missing")
+    let s = row[n]
+    if (!s) throw new Error("Augmented column missing")
+    for (let j = i + 1; j < n; j++) {
+      const coeff = row[j]
+      const sol = x[j]
+      if (!coeff || !sol) continue
+      s = s.sub(coeff.mul(sol))
+    }
+    const pivot = row[i]
+    if (!pivot) throw new Error("Zero pivot on back-substitution")
+    x[i] = s.div(pivot)
+  }
+  return x
+}
+
+export { solveComplex }
diff --git a/lib/math/solveReal.ts b/lib/math/solveReal.ts
new file mode 100644
index 0000000..01479c6
--- /dev/null
+++ b/lib/math/solveReal.ts
@@ -0,0 +1,75 @@
+import { EPS } from "../constants/EPS"
+
+function solveReal(A: number[][], b: number[]) {
+  const n = A.length
+  for (let i = 0; i < n; i++) {
+    const row = A[i]
+    const bi = b[i]
+    if (!row || bi == null) throw new Error("Matrix dimensions mismatch")
+    const copy = row.slice()
+    copy.push(bi)
+    A[i] = copy
+  }
+
+  for (let k = 0; k < n; k++) {
+    let imax = k
+    const pivotRow = A[k]
+    if (!pivotRow) throw new Error("Matrix row missing")
+    let vmax = Math.abs(pivotRow[k] ?? 0)
+    for (let i = k + 1; i < n; i++) {
+      const row = A[i]
+      if (!row) throw new Error("Matrix row missing")
+      const v = Math.abs(row[k] ?? 0)
+      if (v > vmax) {
+        vmax = v
+        imax = i
+      }
+    }
+    if (vmax < EPS) throw new Error("Singular matrix (real)")
+
+    if (imax !== k) {
+      const tmp = A[k]
+      A[k] = A[imax]!
+      A[imax] = tmp!
+    }
+
+    const pivotRowUpdated = A[k]
+    if (!pivotRowUpdated) throw new Error("Pivot row missing")
+    const pivot = pivotRowUpdated[k]
+    if (pivot == null) throw new Error("Zero pivot encountered")
+    for (let i = k + 1; i < n; i++) {
+      const row = A[i]
+      if (!row) throw new Error("Matrix row missing")
+      const entry = row[k]
+      if (entry == null) continue
+      const f = entry / pivot
+      if (Math.abs(f) < EPS) continue
+      for (let j = k; j <= n; j++) {
+        const target = row[j]
+        const source = pivotRowUpdated[j]
+        if (target == null || source == null) continue
+        row[j] = target - f * source
+      }
+    }
+  }
+
+  const x = new Array(n).fill(0)
+  for (let i = n - 1; i >= 0; i--) {
+    const row = A[i]
+    if (!row) throw new Error("Matrix row missing")
+    let s = row[n]
+    if (s == null) throw new Error("Augmented column missing")
+    for (let j = i + 1; j < n; j++) {
+      const coeff = row[j]
+      const sol = x[j]
+      if (coeff == null || sol == null) continue
+      s -= coeff * sol
+    }
+    const pivot = row[i]
+    if (pivot == null) throw new Error("Zero pivot on back-substitution")
+    x[i] = s / pivot
+  }
+  return x
+}
+
+export { solveReal }
diff --git a/lib/parsing/parseNetlist.ts b/lib/parsing/parseNetlist.ts
new file mode 100644
index 0000000..abff1ad
--- /dev/null
+++ b/lib/parsing/parseNetlist.ts
@@ -0,0 +1,380 @@
+import { EPS } from "../constants/EPS"
+
+type Waveform = ((t: number) => number) | null
+
+type PulseSpec = {
+  v1: number
+  v2: number
+  td: number
+  tr: number
+  tf: number
+  ton: number
+  period: number
+  ncycles: number
+}
+
+type ParsedResistor = { name: string; n1: number; n2: number; R: number }
+type ParsedCapacitor = {
+  name: string
+  n1: number
+  n2: number
+  C: number
+  vPrev: number
+}
+type ParsedInductor = {
+  name: string
+  n1: number
+  n2: number
+  L: number
+  iPrev: number
+}
+type ParsedVoltageSource = {
+  name: string
+  n1: number
+  n2: number
+  dc: number
+  acMag: number
+  acPhaseDeg: number
+  waveform: Waveform
+  index: number
+}
+
+type ParsedACAnalysis = {
+  mode: "dec" | "lin"
+  N: number
+  f1: number
+  f2: number
+} | null
+
+type ParsedTranAnalysis = {
+  dt: number
+  tstop: number
+} | null
+
+type ParsedCircuit = {
+  nodes: CircuitNodeIndex
+  R: ParsedResistor[]
+  C: ParsedCapacitor[]
+  L: ParsedInductor[]
+  V: ParsedVoltageSource[]
+  analyses: {
+    ac: ParsedACAnalysis
+    tran: ParsedTranAnalysis
+  }
+  skipped: string[]
+}
+
+class NodeIndex {
+  private map: Map<string, number>
+  rev: string[]
+
+  constructor() {
+    this.map = new Map([["0", 0]])
+    this.rev = ["0"]
+  }
+
+  getOrCreate(name: string) {
+    const key = String(name)
+    if (this.map.has(key)) return this.map.get(key) as number
+    const idx = this.rev.length
+    this.map.set(key, idx)
+    this.rev.push(key)
+    return idx
+  }
+
+  get(name: string) {
+    return this.map.get(String(name))
+  }
+
+  count() {
+    return this.rev.length
+  }
+
+  matrixIndexOfNode(nodeId: number) {
+    if (nodeId === 0) return -1
+    return nodeId - 1
+  }
+}
+
+type CircuitNodeIndex = NodeIndex
+
+function parseNumberWithUnits(raw: unknown) {
+  if (raw == null) return NaN
+  let s = String(raw).trim()
+  if (s === "") return NaN
+  if (/^[+-]?\d*\.?\d+(?:[eE][+-]?\d+)?$/.test(s)) return parseFloat(s)
+  const unitMul = {
+    t: 1e12,
+    g: 1e9,
+    meg: 1e6,
+    k: 1e3,
+    m: 1e-3,
+    u: 1e-6,
+    n: 1e-9,
+    p: 1e-12,
+    f: 1e-15,
+  } as const
+  const m = s.match(/^([+-]?\d*\.?\d+(?:[eE][+-]?\d+)?)([a-zA-Z]+)$/)
+  if (!m) return parseFloat(s)
+  const [, numberPart, suffixPart] = m
+  if (numberPart == null) return parseFloat(s)
+  let val = parseFloat(numberPart)
+  let suf = (suffixPart ?? "").toLowerCase()
+  suf = suf.replace(/(ohm|v|a|s|h|f)$/g, "")
+  if (suf === "meg") return val * unitMul.meg
+  if (suf.length === 1 && suf in unitMul) {
+    const key = suf as keyof typeof unitMul
+    return val * unitMul[key]
+  }
+  return val
+}
+
+function smartTokens(line: string) {
+  const re = /"[^"]*"|\([^()]*\)|\S+/g
+  const out: string[] = []
+  let m: RegExpExecArray | null
+  while ((m = re.exec(line)) !== null) out.push(m[0])
+  return out
+}
+
+function requireToken(tokens: string[], index: number, context: string) {
+  const token = tokens[index]
+  if (token == null) throw new Error(context)
+  return token
+}
+
+function parsePulseArgs(token: string): PulseSpec {
+  const clean = token.trim().replace(/^pulse\s*\(/i, "(")
+  const inside = clean.replace(/^\(/, "").replace(/\)$/, "").trim()
+  const parts = inside.split(/[\s,]+/).filter((x) => x.length)
+  if (parts.length < 7) throw new Error("PULSE(...) requires 7 or 8 args")
+  const vals = parts.map((value) => parseNumberWithUnits(value))
+  if (vals.some((v) => Number.isNaN(v)))
+    throw new Error("Invalid PULSE() numeric value")
+  return {
+    v1: vals[0]!,
+    v2: vals[1]!,
+    td: vals[2]!,
+    tr: vals[3]!,
+    tf: vals[4]!,
+    ton: vals[5]!,
+    period: vals[6]!,
+    ncycles: parts[7] != null ? vals[7]! : Infinity,
+  }
+}
+
+function pulseValue(p: PulseSpec, t: number) {
+  if (t < p.td) return p.v1
+  const tt = t - p.td
+  const cyclesDone = Math.floor(tt / p.period)
+  if (cyclesDone >= p.ncycles) return p.v1
+  const tc = tt - cyclesDone * p.period
+  if (tc < p.tr) {
+    const a = tc / Math.max(p.tr, EPS)
+    return p.v1 + (p.v2 - p.v1) * a
+  }
+  if (tc < p.tr + p.ton) {
+    return p.v2
+  }
+  if (tc < p.tr + p.ton + p.tf) {
+    const a = (tc - (p.tr + p.ton)) / Math.max(p.tf, EPS)
+    return p.v2 + (p.v1 - p.v2) * a
+  }
+  return p.v1
+}
+
+function parseNetlist(text: string): ParsedCircuit {
+  const ckt: ParsedCircuit = {
+    nodes: new NodeIndex(),
+    R: [],
+    C: [],
+    L: [],
+    V: [],
+    analyses: { ac: null, tran: null },
+    skipped: [],
+  }
+
+  const lines = text.split(/\r?\n/)
+  let seenTitle = false
+
+  for (const raw of lines) {
+    let line = raw.trim()
+    if (!line) continue
+    if (/^\*/.test(line)) continue
+    if (/^\s*\.end\b/i.test(line)) break
+    line = line.replace(/\/\/.*$/, "")
+    line = line.replace(/;.*$/, "")
+
+    const tokens = smartTokens(line)
+    if (tokens.length === 0) continue
+
+    const first = tokens[0]!
+    if (first.length === 0) continue
+
+    if (!seenTitle && !/^[rclvgmiqd]\w*$/i.test(first) && !/^\./.test(first)) {
+      seenTitle = true
+      continue
+    }
+
+    if (/^\./.test(first)) {
+      const dir = first.toLowerCase()
+      if (dir === ".ac") {
+        const mode = requireToken(tokens, 1, ".ac missing mode").toLowerCase()
+        if (mode !== "dec" && mode !== "lin")
+          throw new Error(".ac supports 'dec' or 'lin'")
+        const N = parseInt(
+          requireToken(tokens, 2, ".ac missing point count"),
+          10,
+        )
+        const f1 = parseNumberWithUnits(
+          requireToken(tokens, 3, ".ac missing start frequency"),
+        )
+        const f2 = parseNumberWithUnits(
+          requireToken(tokens, 4, ".ac missing stop frequency"),
+        )
+        ckt.analyses.ac = { mode, N, f1, f2 }
+      } else if (dir === ".tran") {
+        const dt = parseNumberWithUnits(
+          requireToken(tokens, 1, ".tran missing timestep"),
+        )
+        const tstop = parseNumberWithUnits(
+          requireToken(tokens, 2, ".tran missing stop time"),
+        )
+        ckt.analyses.tran = { dt, tstop }
+      } else {
+        ckt.skipped.push(line)
+      }
+      continue
+    }
+
+    const typeChar = first.charAt(0).toLowerCase()
+    const name = first
+
+    try {
+      if (typeChar === "r") {
+        const n1 = ckt.nodes.getOrCreate(
+          requireToken(tokens, 1, "Resistor missing node"),
+        )
+        const n2 = ckt.nodes.getOrCreate(
+          requireToken(tokens, 2, "Resistor missing node"),
+        )
+        const val = parseNumberWithUnits(
+          requireToken(tokens, 3, "Resistor missing value"),
+        )
+        ckt.R.push({ name, n1, n2, R: val })
+      } else if (typeChar === "c") {
+        const n1 = ckt.nodes.getOrCreate(
+          requireToken(tokens, 1, "Capacitor missing node"),
+        )
+        const n2 = ckt.nodes.getOrCreate(
+          requireToken(tokens, 2, "Capacitor missing node"),
+        )
+        const val = parseNumberWithUnits(
+          requireToken(tokens, 3, "Capacitor missing value"),
+        )
+        ckt.C.push({ name, n1, n2, C: val, vPrev: 0 })
+      } else if (typeChar === "l") {
+        const n1 = ckt.nodes.getOrCreate(
+          requireToken(tokens, 1, "Inductor missing node"),
+        )
+        const n2 = ckt.nodes.getOrCreate(
+          requireToken(tokens, 2, "Inductor missing node"),
+        )
+        const val = parseNumberWithUnits(
+          requireToken(tokens, 3, "Inductor missing value"),
+        )
+        ckt.L.push({ name, n1, n2, L: val, iPrev: 0 })
+      } else if (typeChar === "v") {
+        const n1 = ckt.nodes.getOrCreate(
+          requireToken(tokens, 1, "Voltage source missing node"),
+        )
+        const n2 = ckt.nodes.getOrCreate(
+          requireToken(tokens, 2, "Voltage source missing node"),
+        )
+        const spec: Omit<
+          ParsedVoltageSource,
+          "name" | "n1" | "n2" | "index"
+        > & { index?: number } = {
+          dc: 0,
+          acMag: 0,
+          acPhaseDeg: 0,
+          waveform: null,
+          index: -1,
+        }
+        let i = 3
+        if (i < tokens.length && !/^[a-zA-Z]/.test(tokens[i]!)) {
+          spec.dc = parseNumberWithUnits(tokens[i]!)
+          i++
+        }
+        while (i < tokens.length) {
+          const key = tokens[i]!.toLowerCase()
+          if (key === "dc") {
+            const valueToken = requireToken(tokens, i + 1, "DC value missing")
+            spec.dc = parseNumberWithUnits(valueToken)
+            i += 2
+          } else if (key === "ac") {
+            const magToken = requireToken(tokens, i + 1, "AC magnitude missing")
+            spec.acMag = parseNumberWithUnits(magToken)
+            const phaseToken = tokens[i + 2]
+            if (phaseToken != null && /^[+-]?\d/.test(phaseToken)) {
+              spec.acPhaseDeg = parseNumberWithUnits(phaseToken)
+              i += 3
+            } else {
+              i += 2
+            }
+          } else if (key.startsWith("pulse")) {
+            const argToken = key.includes("(")
+              ? key
+              : requireToken(tokens, i + 1, "PULSE() missing arguments")
+            if (!argToken || !/\(.*\)/.test(argToken))
+              throw new Error("Malformed PULSE() specification")
+            const p = parsePulseArgs(argToken)
+            spec.waveform = (t: number) => pulseValue(p, t)
+            i += key.includes("(") ? 1 : 2
+          } else if (/^\(.*\)$/.test(key)) {
+            i++
+          } else {
+            i++
+          }
+        }
+        ckt.V.push({
+          name,
+          n1,
+          n2,
+          dc: spec.dc,
+          acMag: spec.acMag,
+          acPhaseDeg: spec.acPhaseDeg,
+          waveform: spec.waveform,
+          index: spec.index ?? -1,
+        })
+      } else {
+        ckt.skipped.push(line)
+      }
+    } catch (err) {
+      if (err instanceof Error) {
+        throw new Error(`Parse error on line: "${line}"\n${err.message}`)
+      }
+      throw err
+    }
+  }
+
+  const nNodes = ckt.nodes.count() - 1
+  for (let i = 0; i < ckt.V.length; i++) {
+    const vs = ckt.V[i]
+    if (!vs) continue
+    vs.index = nNodes + i
+  }
+  return ckt
+}
+
+export type {
+  ParsedCircuit,
+  ParsedACAnalysis,
+  ParsedTranAnalysis,
+  ParsedResistor,
+  ParsedCapacitor,
+  ParsedInductor,
+  ParsedVoltageSource,
+  CircuitNodeIndex,
+}
+export { parseNetlist }
diff --git a/lib/stamping/stampAdmittanceComplex.ts b/lib/stamping/stampAdmittanceComplex.ts
new file mode 100644
index 0000000..f98949e
--- /dev/null
+++ b/lib/stamping/stampAdmittanceComplex.ts
@@ -0,0 +1,32 @@
+import type { CircuitNodeIndex } from "../parsing/parseNetlist"
+import { Complex } from "../math/Complex"
+
+function stampAdmittanceComplex(
+  A: Complex[][],
+  nidx: CircuitNodeIndex,
+  n1: number,
+  n2: number,
+  Y: Complex,
+) {
+  const i1 = nidx.matrixIndexOfNode(n1)
+  const i2 = nidx.matrixIndexOfNode(n2)
+  if (i1 >= 0) {
+    const row1 = A[i1]
+    if (!row1) throw new Error("Matrix row missing while stamping")
+    row1[i1] = row1[i1]?.add(Y) ?? Y
+  }
+  if (i2 >= 0) {
+    const row2 = A[i2]
+    if (!row2) throw new Error("Matrix row missing while stamping")
+    row2[i2] = row2[i2]?.add(Y) ?? Y
+  }
+  if (i1 >= 0 && i2 >= 0) {
+    const row1 = A[i1]
+    const row2 = A[i2]
+    if (!row1 || !row2) throw new Error("Matrix row missing while stamping")
+    row1[i2] = row1[i2]?.sub(Y) ?? Complex.from(0, 0).sub(Y)
+    row2[i1] = row2[i1]?.sub(Y) ?? Complex.from(0, 0).sub(Y)
+  }
+}
+
+export { stampAdmittanceComplex }
diff --git a/lib/stamping/stampAdmittanceReal.ts b/lib/stamping/stampAdmittanceReal.ts
new file mode 100644
index 0000000..ad04dc3
--- /dev/null
+++ b/lib/stamping/stampAdmittanceReal.ts
@@ -0,0 +1,31 @@
+import type { CircuitNodeIndex } from "../parsing/parseNetlist"
+
+function stampAdmittanceReal(
+  A: number[][],
+  nidx: CircuitNodeIndex,
+  n1: number,
+  n2: number,
+  Y: number,
+) {
+  const i1 = nidx.matrixIndexOfNode(n1)
+  const i2 = nidx.matrixIndexOfNode(n2)
+  if (i1 >= 0) {
+    const row1 = A[i1]
+    if (!row1) throw new Error("Matrix row missing while stamping")
+    row1[i1] = (row1[i1] ?? 0) + Y
+  }
+  if (i2 >= 0) {
+    const row2 = A[i2]
+    if (!row2) throw new Error("Matrix row missing while stamping")
+    row2[i2] = (row2[i2] ?? 0) + Y
+  }
+  if (i1 >= 0 && i2 >= 0) {
+    const row1 = A[i1]
+    const row2 = A[i2]
+    if (!row1 || !row2) throw new Error("Matrix row missing while stamping")
+    row1[i2] = (row1[i2] ?? 0) - Y
+    row2[i1] = (row2[i1] ?? 0) - Y
+  }
+}
+
+export { stampAdmittanceReal }
diff --git a/lib/stamping/stampCurrentComplex.ts b/lib/stamping/stampCurrentComplex.ts
new file mode 100644
index 0000000..a9edfa9
--- /dev/null
+++ b/lib/stamping/stampCurrentComplex.ts
@@ -0,0 +1,17 @@
+import type { CircuitNodeIndex } from "../parsing/parseNetlist"
+import { Complex } from "../math/Complex"
+
+function stampCurrentComplex(
+  b: Complex[],
+  nidx: CircuitNodeIndex,
+  nPlus: number,
+  nMinus: number,
+  current: Complex,
+) {
+  const iPlus = nidx.matrixIndexOfNode(nPlus)
+  const iMinus = nidx.matrixIndexOfNode(nMinus)
+  if (iPlus >= 0) b[iPlus] = (b[iPlus] ?? Complex.from(0, 0)).sub(current)
+  if (iMinus >= 0) b[iMinus] = (b[iMinus] ?? Complex.from(0, 0)).add(current)
+}
+
+export { stampCurrentComplex }
diff --git a/lib/stamping/stampCurrentReal.ts b/lib/stamping/stampCurrentReal.ts
new file mode 100644
index 0000000..1f5693f
--- /dev/null
+++ b/lib/stamping/stampCurrentReal.ts
@@ -0,0 +1,16 @@
+import type { CircuitNodeIndex } from "../parsing/parseNetlist"
+
+function stampCurrentReal(
+  b: number[],
+  nidx: CircuitNodeIndex,
+  nPlus: number,
+  nMinus: number,
+  current: number,
+) {
+  const iPlus = nidx.matrixIndexOfNode(nPlus)
+  const iMinus = nidx.matrixIndexOfNode(nMinus)
+  if (iPlus >= 0) b[iPlus] = (b[iPlus] ?? 0) - current
+  if (iMinus >= 0) b[iMinus] = (b[iMinus] ?? 0) + current
+}
+
+export { stampCurrentReal }
diff --git a/lib/stamping/stampVoltageSourceComplex.ts b/lib/stamping/stampVoltageSourceComplex.ts
new file mode 100644
index 0000000..641a722
--- /dev/null
+++ b/lib/stamping/stampVoltageSourceComplex.ts
@@ -0,0 +1,37 @@
+import type { ParsedVoltageSource } from "../parsing/parseNetlist"
+import type { CircuitNodeIndex } from "../parsing/parseNetlist"
+import { Complex } from "../math/Complex"
+
+function stampVoltageSourceComplex(
+  A: Complex[][],
+  b: Complex[],
+  nidx: CircuitNodeIndex,
+  source: ParsedVoltageSource,
+  voltage: Complex,
+) {
+  const i1 = nidx.matrixIndexOfNode(source.n1)
+  const i2 = nidx.matrixIndexOfNode(source.n2)
+  const j = source.index
+  const one = Complex.from(1, 0)
+  if (i1 >= 0) {
+    const row1 = A[i1]
+    if (!row1)
+      throw new Error("Matrix row missing while stamping voltage source")
+    row1[j] = row1[j]?.add(one) ?? one
+  }
+  if (i2 >= 0) {
+    const row2 = A[i2]
+    if (!row2)
+      throw new Error("Matrix row missing while stamping voltage source")
+    row2[j] = row2[j]?.sub(one) ?? Complex.from(0, 0).sub(one)
+  }
+  const branchRow = A[j]
+  if (!branchRow)
+    throw new Error("Branch row missing while stamping voltage source")
+  if (i1 >= 0) branchRow[i1] = branchRow[i1]?.add(one) ?? one
+  if (i2 >= 0)
+    branchRow[i2] = branchRow[i2]?.sub(one) ?? Complex.from(0, 0).sub(one)
+  b[j] = (b[j] ?? Complex.from(0, 0)).add(voltage)
+}
+
+export { stampVoltageSourceComplex }
diff --git a/lib/stamping/stampVoltageSourceReal.ts b/lib/stamping/stampVoltageSourceReal.ts
new file mode 100644
index 0000000..cf6613a
--- /dev/null
+++ b/lib/stamping/stampVoltageSourceReal.ts
@@ -0,0 +1,34 @@
+import type { ParsedVoltageSource } from "../parsing/parseNetlist"
+import type { CircuitNodeIndex } from "../parsing/parseNetlist"
+
+function stampVoltageSourceReal(
+  A: number[][],
+  b: number[],
+  nidx: CircuitNodeIndex,
+  source: ParsedVoltageSource,
+  voltage: number,
+) {
+  const i1 = nidx.matrixIndexOfNode(source.n1)
+  const i2 = nidx.matrixIndexOfNode(source.n2)
+  const j = source.index
+  if (i1 >= 0) {
+    const row1 = A[i1]
+    if (!row1)
+      throw new Error("Matrix row missing while stamping voltage source")
+    row1[j] = (row1[j] ?? 0) + 1
+  }
+  if (i2 >= 0) {
+    const row2 = A[i2]
+    if (!row2)
+      throw new Error("Matrix row missing while stamping voltage source")
+    row2[j] = (row2[j] ?? 0) - 1
+  }
+  const branchRow = A[j]
+  if (!branchRow)
+    throw new Error("Branch row missing while stamping voltage source")
+  if (i1 >= 0) branchRow[i1] = (branchRow[i1] ?? 0) + 1
+  if (i2 >= 0) branchRow[i2] = (branchRow[i2] ?? 0) - 1
+  b[j] = (b[j] ?? 0) + voltage
+}
+
+export { stampVoltageSourceReal }
diff --git a/lib/utils/logspace.ts b/lib/utils/logspace.ts
new file mode 100644
index 0000000..f42f316
--- /dev/null
+++ b/lib/utils/logspace.ts
@@ -0,0 +1,17 @@
+import { EPS } from "../constants/EPS"
+
+function logspace(f1: number, f2: number, pointsPerDecade: number) {
+  if (f1 <= 0 || f2 <= 0) throw new Error(".ac frequencies must be > 0")
+  if (f2 < f1) [f1, f2] = [f2, f1]
+  const decades = Math.log10(f2 / f1)
+  const n = Math.max(1, Math.ceil(decades * pointsPerDecade))
+  const arr: number[] = []
+  for (let i = 0; i <= n; i++) {
+    arr.push(f1 * Math.pow(10, i / pointsPerDecade))
+  }
+  const last = arr[arr.length - 1]
+  if (last == null || last < f2 * (1 - EPS)) arr.push(f2)
+  return arr
+}
+
+export { logspace }
diff --git a/package.json b/package.json
index 5e3f25d..63c5112 100644
--- a/package.json
+++ b/package.json
@@ -3,6 +3,9 @@
   "main": "dist/index.js",
   "type": "module",
   "version": "0.0.1",
+  "files": [
+    "dist"
+  ],
   "scripts": {
     "build": "tsup-node ./lib/index.ts",
     "format": "biome format --write .",