diff --git a/src/frontend/Makefile.am b/src/frontend/Makefile.am
index 0d9a622af..5eee39fce 100644
--- a/src/frontend/Makefile.am
+++ b/src/frontend/Makefile.am
@@ -55,6 +55,9 @@ libfte_la_SOURCES = \
     com_loadpull.h  \
     com_optimize.c  \
     com_optimize.h  \
+    snp2va.c  \
+    snp2va.h  \
+    com_presnp.c  \
     com_qpss.c  \
     com_qpss.h  \
     com_qpac.c  \
diff --git a/src/frontend/commands.c b/src/frontend/commands.c
index 1cfa7002c..b161e0973 100644
--- a/src/frontend/commands.c
+++ b/src/frontend/commands.c
@@ -93,6 +93,7 @@
 #include "spiceif.h" /* for com_snload() and com_snsave() */
 
 #include "com_dl.h"
+#include "snp2va.h"
 
 #ifdef XSPICE
 /* gtri - begin - wbk - add include files */
@@ -299,6 +300,10 @@ struct comm spcp_coms[] = {
       NULL,
       "library library ... : Loads a osdi library." } ,
 #endif
+    { "snp", com_pre_snp, FALSE, TRUE,          /* Enhancement-200 (use as `pre_snp`) */
+      { 040000, 040000, 040000, 040000 }, E_BEGINNING, 1, LOTS,
+      NULL,
+      "file.sNp [module] : compile a Touchstone S-parameter file to a Verilog-A n-port OSDI model. Use as `pre_snp file.sNp` (runs before circuit parsing); load the .osdi it writes with `pre_osdi file.osdi`." },
 #ifdef DEVLIB
     { "use", com_use, FALSE, TRUE,
       { 040000, 040000, 040000, 040000 }, E_BEGINNING, 1, LOTS,
diff --git a/src/frontend/inp.c b/src/frontend/inp.c
index f0de1289d..971696e95 100644
--- a/src/frontend/inp.c
+++ b/src/frontend/inp.c
@@ -842,13 +842,24 @@ inp_spsource(FILE *fp, bool comfile, char *filename, bool intfile)
         /* Now that the deck is loaded, do the pre commands, if there are any,
            before the circuit structure is set up */
         if (pre_controls) {
+            int pass;
             pre_controls = wl_reverse(pre_controls);
-            for (wl = pre_controls; wl; wl = wl->wl_next){
+            /* Enhancement-200: run every pre_snp (stripped to "snp") before any
+             * other pre_ command -- notably pre_osdi (stripped to "osdi") -- so the
+             * .osdi that pre_snp generates already exists when pre_osdi loads it.
+             * Pass 0 executes the snp commands (in deck order), pass 1 the rest. */
+            for (pass = 0; pass < 2; pass++) {
+                for (wl = pre_controls; wl; wl = wl->wl_next) {
+                    int is_snp = ciprefix("snp ", wl->wl_word) ||
+                                 strcasecmp(wl->wl_word, "snp") == 0;
+                    if (pass == 0 ? !is_snp : is_snp)
+                        continue;
 #ifdef OSDI
-                inputdir = dir_name;
+                    inputdir = dir_name;
 #endif
-                /* process each pre_xxx command */
-                cp_evloop(wl->wl_word);
+                    /* process each pre_xxx command */
+                    cp_evloop(wl->wl_word);
+                }
             }
 
 #ifdef OSDI
diff --git a/src/frontend/inpcom.c b/src/frontend/inpcom.c
index 0f4c905da..d07706896 100644
--- a/src/frontend/inpcom.c
+++ b/src/frontend/inpcom.c
@@ -1809,6 +1809,10 @@ static struct inp_read_t inp_read(FILE* fp, int call_depth, const char* dir_name
                     ciprefix("codemodel", buffer) ||
                     ciprefix("osdi", buffer) ||
                     ciprefix("pre_osdi", buffer) ||
+
+                    ciprefix("snp", buffer) ||
+
+                    ciprefix("pre_snp", buffer) ||
                     ciprefix("echo", buffer) || ciprefix("shell", buffer) ||
                     ciprefix("source", buffer) || ciprefix("cd", buffer) ||
                     ciprefix("load", buffer) || ciprefix("setcs", buffer) ||
diff --git a/src/frontend/snp2va.c b/src/frontend/snp2va.c
new file mode 100644
index 000000000..4404f769b
--- /dev/null
+++ b/src/frontend/snp2va.c
@@ -0,0 +1,566 @@
+/* Enhancement-200: Touchstone (.sNp) -> Verilog-A n-port converter, in C.
+ *
+ * A C port of the pure-Python snp2va.py (Enhancement-199): parse a Touchstone
+ * S-parameter file, convert S -> Y, fit every Y_ij(f) with a common-pole rational
+ * (Gustavsen vector fitting), and emit a Verilog-A n-port realized with laplace_nd,
+ * so through OpenVAF/OSDI it works in AC and transient. Used by the `pre_snp`
+ * front-end command, which then invokes openvaf-r to compile the emitted .va.
+ *
+ * The numerical core is self-contained (only stdio/stdlib/string/math/complex),
+ * so it does not pull in ngspice's own complex.h. Public entry point:
+ *     int snp2va_convert(const char *snp, const char *va, const char *module,
+ *                        char *msg, int msglen);
+ * returns 0 on success, non-zero on failure (msg gets a one-line status).
+ */
+#ifdef _MSC_VER
+#include <iostream>
+#include <complex>
+
+extern "C" {
+#include "snp2va.h"
+#include <string.h>
+};
+#else
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include <complex.h>
+#include "snp2va.h"
+#endif
+
+#ifdef _MSC_VER
+typedef std::complex<double> cplx;
+#define cabs abs
+#define cpow pow
+#define cimag imag
+#define creal real
+#define cexp exp
+#define I cplx(0.0, 1.0)
+#define strcasecmp _stricmp
+#else
+typedef double _Complex cplx;
+#endif
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+/* ============================ linear algebra ============================ */
+
+/* Least-squares min||A x - b|| for a REAL overdetermined system (m>=n) via
+ * Householder QR. A is row-major m*n, b length m, x length n. Returns 0 on ok. */
+static int lstsq_real(double *A, double *b, int m, int n, double *x)
+{
+    int i, j, k;
+    for (k = 0; k < n; k++) {
+        double norm = 0.0;
+        for (i = k; i < m; i++) norm += A[i*n+k]*A[i*n+k];
+        norm = sqrt(norm);
+        if (norm == 0.0) continue;
+        double alpha = (A[k*n+k] >= 0.0) ? -norm : norm;
+        double *v = (double*) calloc((size_t) m, sizeof(double));
+        v[k] = A[k*n+k] - alpha;
+        for (i = k+1; i < m; i++) v[i] = A[i*n+k];
+        double vn2 = 0.0;
+        for (i = k; i < m; i++) vn2 += v[i]*v[i];
+        if (vn2 == 0.0) { free(v); continue; }
+        for (j = k; j < n; j++) {
+            double s = 0.0;
+            for (i = k; i < m; i++) s += v[i]*A[i*n+j];
+            s = s*2.0/vn2;
+            for (i = k; i < m; i++) A[i*n+j] -= s*v[i];
+        }
+        double s = 0.0;
+        for (i = k; i < m; i++) s += v[i]*b[i];
+        s = s*2.0/vn2;
+        for (i = k; i < m; i++) b[i] -= s*v[i];
+        free(v);
+    }
+    for (i = n-1; i >= 0; i--) {
+        double acc = b[i];
+        for (j = i+1; j < n; j++) acc -= A[i*n+j]*x[j];
+        x[i] = (A[i*n+i] != 0.0) ? acc/A[i*n+i] : 0.0;
+    }
+    return 0;
+}
+
+/* In-place inverse of an n x n complex matrix (row-major), Gauss-Jordan w/ pivot.
+ * Returns 0 on ok, 1 if singular. */
+static int mat_inv_c(cplx *M, int n)
+{
+    int i, j, c, p;
+    cplx *A = (cplx*) malloc((size_t) n*2*n*sizeof(cplx));
+    for (i = 0; i < n; i++) {
+        for (j = 0; j < n; j++) A[i*2*n+j] = M[i*n+j];
+        for (j = 0; j < n; j++) A[i*2*n+n+j] = (i==j) ? 1.0 : 0.0;
+    }
+    for (c = 0; c < n; c++) {
+        p = c; double best = cabs(A[c*2*n+c]);
+        for (i = c+1; i < n; i++) { double v = cabs(A[i*2*n+c]); if (v > best) { best = v; p = i; } }
+        if (best == 0.0) { free(A); return 1; }
+        if (p != c) for (j = 0; j < 2*n; j++) { cplx t = A[c*2*n+j]; A[c*2*n+j]=A[p*2*n+j]; A[p*2*n+j]=t; }
+        cplx d = A[c*2*n+c];
+        for (j = 0; j < 2*n; j++) A[c*2*n+j] /= d;
+        for (i = 0; i < n; i++) if (i != c) {
+            cplx f = A[i*2*n+c];
+            if (f != 0.0) for (j = 0; j < 2*n; j++) A[i*2*n+j] -= f*A[c*2*n+j];
+        }
+    }
+    for (i = 0; i < n; i++) for (j = 0; j < n; j++) M[i*n+j] = A[i*2*n+n+j];
+    free(A);
+    return 0;
+}
+
+/* Monic polynomial (DESCENDING, len nr+1) from roots. coef must hold nr+1. */
+static void poly_from_roots(const cplx *roots, int nr, cplx *coef)
+{
+    int i, k;
+    coef[0] = 1.0;
+    for (i = 1; i <= nr; i++) coef[i] = 0.0;
+    for (k = 0; k < nr; k++) {
+        for (i = k+1; i >= 1; i--) coef[i] = coef[i] - roots[k]*coef[i-1];
+    }
+}
+
+static cplx poly_eval(const cplx *c, int deg, cplx x)
+{
+    cplx r = 0.0; int i;
+    for (i = 0; i <= deg; i++) r = r*x + c[i];
+    return r;
+}
+
+/* All roots of a DESCENDING-coeff polynomial (deg = len-1) via Durand-Kerner.
+ * roots[] must hold deg. Returns 0 on ok. */
+static int poly_roots(const cplx *cin, int deg, cplx *roots)
+{
+    int i, j, it;
+    if (deg <= 0) return 0;
+    cplx *c = (cplx*) malloc((size_t)(deg+1)*sizeof(cplx));
+    for (i = 0; i <= deg; i++) c[i] = cin[i] / cin[0];   /* monic */
+    cplx seed = 0.4 + 0.9*I;
+    for (i = 0; i < deg; i++) roots[i] = cpow(seed, (double) i);
+    for (it = 0; it < 500; it++) {
+        double maxstep = 0.0;
+        for (i = 0; i < deg; i++) {
+            cplx num = poly_eval(c, deg, roots[i]);
+            cplx den = 1.0;
+            for (j = 0; j < deg; j++) if (j != i) den *= (roots[i]-roots[j]);
+            cplx step = (cabs(den) > 1e-300) ? num/den : 0.0;
+            roots[i] -= step;
+            if (cabs(step) > maxstep) maxstep = cabs(step);
+        }
+        if (maxstep < 1e-14) break;
+    }
+    free(c);
+    return 0;
+}
+
+/* ============================ Touchstone I/O ============================ */
+
+typedef struct { double *freqs; cplx *S; int nf; int N; double z0; char ptype; } TS;
+
+static void ts_free(TS *t) { free(t->freqs); free(t->S); }
+
+/* returns 0 on ok */
+static int parse_touchstone(const char *fn, TS *out, char *msg, int msglen)
+{
+    FILE *f = fopen(fn, "r");
+    if (!f) { snprintf(msg, (size_t) msglen, "cannot open '%s'", fn); return 1; }
+    double fmul = 1e9, z0 = 50.0; char ptype = 'S'; char fmt[3] = "MA";
+    /* collect all numeric tokens after the '#' options line(s) */
+    double *nums = NULL; long ncap = 0, nn = 0;
+    char line[4096];
+    while (fgets(line, sizeof line, f)) {
+        char *h = strchr(line, '!'); if (h) *h = '\0';
+        char *p = line;
+        while (*p && isspace((unsigned char)*p)) p++;
+        if (*p == '\0') continue;
+        if (*p == '#') {
+            char *tok = strtok(p+1, " \t\r\n");
+            while (tok) {
+                if (!strcasecmp(tok,"HZ")) fmul=1.0;
+                else if (!strcasecmp(tok,"KHZ")) fmul=1e3;
+                else if (!strcasecmp(tok,"MHZ")) fmul=1e6;
+                else if (!strcasecmp(tok,"GHZ")) fmul=1e9;
+                else if (!strcasecmp(tok,"S")||!strcasecmp(tok,"Y")||!strcasecmp(tok,"Z")) ptype=(char)toupper((unsigned char)tok[0]);
+                else if (!strcasecmp(tok,"MA")||!strcasecmp(tok,"DB")||!strcasecmp(tok,"RI")) { fmt[0]=(char)toupper((unsigned char)tok[0]); fmt[1]=(char)toupper((unsigned char)tok[1]); }
+                else if (!strcasecmp(tok,"R")) { char *z=strtok(NULL," \t\r\n"); if (z) z0=atof(z); }
+                tok = strtok(NULL, " \t\r\n");
+            }
+            continue;
+        }
+        /* numeric data line */
+        char *tok = strtok(p, " \t\r\n");
+        while (tok) {
+            char *end; double v = strtod(tok, &end);
+            if (end != tok) {
+                if (nn >= ncap) { ncap = ncap ? ncap*2 : 1024; nums = (double*) realloc(nums, (size_t) ncap*sizeof(double)); }
+                nums[nn++] = v;
+            }
+            tok = strtok(NULL, " \t\r\n");
+        }
+    }
+    fclose(f);
+    /* infer port count N: try the file extension .sNp, else brute force */
+    int N = 0;
+    const char *dot = strrchr(fn, '.');
+    if (dot && (dot[1]=='s'||dot[1]=='S') && (fn[strlen(fn)-1]=='p'||fn[strlen(fn)-1]=='P')) {
+        N = atoi(dot+2);
+    }
+    if (N <= 0) {
+        int c;
+        for (c = 1; c <= 16; c++) if (nn % (1 + 2*c*c) == 0) { N = c; break; }
+    }
+    if (N <= 0) { free(nums); snprintf(msg,(size_t)msglen,"cannot determine port count"); return 1; }
+    int rec = 1 + 2*N*N;
+    int nf = (int)(nn / rec);
+    if (nf < 2) { free(nums); snprintf(msg,(size_t)msglen,"too few frequency points (%d)", nf); return 1; }
+    out->freqs = (double*) malloc((size_t) nf*sizeof(double));
+    out->S = (cplx*) malloc((size_t) nf*N*N*sizeof(cplx));
+    out->nf = nf; out->N = N; out->z0 = z0; out->ptype = ptype;
+    int r, kk;
+    for (r = 0; r < nf; r++) {
+        double *chunk = nums + (long) r*rec;
+        out->freqs[r] = chunk[0]*fmul;
+        double *vals = chunk+1;
+        cplx pv[16*16];
+        for (kk = 0; kk < N*N; kk++) {
+            double a = vals[2*kk], b = vals[2*kk+1];
+            if (!strcmp(fmt,"MA")) pv[kk] = a*cexp(I*b*M_PI/180.0);
+            else if (!strcmp(fmt,"DB")) pv[kk] = pow(10.0,a/20.0)*cexp(I*b*M_PI/180.0);
+            else pv[kk] = a + I*b;
+        }
+        /* Touchstone: N=2 order is S11 S21 S12 S22; general is row-major */
+        cplx *M = out->S + (long) r*N*N;
+        if (N == 2) { M[0]=pv[0]; M[2]=pv[1]; M[1]=pv[2]; M[3]=pv[3]; }
+        else for (kk = 0; kk < N*N; kk++) M[kk] = pv[kk];
+    }
+    free(nums);
+    return 0;
+}
+
+/* S/Y/Z -> Y (row-major per frequency), Yout must hold nf*N*N. */
+static int to_Y(const TS *t, cplx *Yout)
+{
+    int N = t->N, r, i, j;
+    cplx *tmp = (cplx*) malloc((size_t) N*N*sizeof(cplx));
+    for (r = 0; r < t->nf; r++) {
+        const cplx *M = t->S + (long) r*N*N;
+        cplx *Y = Yout + (long) r*N*N;
+        if (t->ptype == 'Y') { for (i=0;i<N*N;i++) Y[i]=M[i]; }
+        else if (t->ptype == 'Z') { for (i=0;i<N*N;i++) Y[i]=M[i]; if (mat_inv_c(Y,N)) { free(tmp); return 1; } }
+        else { /* S -> Y = (1/z0)(I-S)(I+S)^-1 */
+            cplx *IpS = tmp;
+            for (i=0;i<N;i++) for (j=0;j<N;j++) IpS[i*N+j] = ((i==j)?1.0:0.0) + M[i*N+j];
+            if (mat_inv_c(IpS,N)) { free(tmp); return 1; }
+            for (i=0;i<N;i++) for (j=0;j<N;j++) {
+                cplx acc = 0.0; int k;
+                for (k=0;k<N;k++) acc += (((i==k)?1.0:0.0) - M[i*N+k]) * IpS[k*N+j];
+                Y[i*N+j] = acc / t->z0;
+            }
+        }
+    }
+    free(tmp);
+    return 0;
+}
+
+/* ============================ vector fitting ============================ */
+/* layout: for each pole, 0='real', 1='cc-start' (its conjugate is the next). */
+static int build_layout(const cplx *poles, int Np, int *lay)
+{
+    int i = 0, m = 0;
+    while (i < Np) {
+        if (fabs(cimag(poles[i])) < 1e-9*fabs(creal(poles[i]))+1e-30) { lay[m++] = 0; i += 1; }
+        else { lay[m++] = 1; i += 2; }
+    }
+    return m;   /* number of blocks */
+}
+
+/* complex partial-fraction basis (real-valued cc combos), Ns x Np, row-major */
+static void build_basis(const cplx *s, int Ns, const cplx *poles, int Np, cplx *A)
+{
+    int r, i;
+    for (r = 0; r < Ns; r++) {
+        i = 0;
+        while (i < Np) {
+            if (fabs(cimag(poles[i])) < 1e-9*fabs(creal(poles[i]))+1e-30) {
+                A[r*Np+i] = 1.0/(s[r]-poles[i]); i += 1;
+            } else {
+                cplx p = poles[i];
+                A[r*Np+i]   = 1.0/(s[r]-p) + 1.0/(s[r]-conj(p));
+                A[r*Np+i+1] = I/(s[r]-p) - I/(s[r]-conj(p));
+                i += 2;
+            }
+        }
+    }
+}
+
+/* complex residues from real ctil coeffs, per real/cc layout */
+static void ctil_to_cres(const double *ctil, const cplx *poles, int Np, cplx *cres)
+{
+    int i = 0;
+    while (i < Np) {
+        if (fabs(cimag(poles[i])) < 1e-9*fabs(creal(poles[i]))+1e-30) { cres[i] = ctil[i]; i += 1; }
+        else { cres[i] = ctil[i] + I*ctil[i+1]; cres[i+1] = conj(cres[i]); i += 2; }
+    }
+}
+
+/* One vector-fit run (fixed pole count). s,F normalized. Returns fit in poles/res/d/e. */
+static void vector_fit(const cplx *s, int Ns, const cplx *F, int Nf, int Np,
+                       cplx *poles, cplx *res, double *d, double *e, int n_iter)
+{
+    int iter, i, j, k, r;
+    cplx *A = (cplx*) malloc((size_t) Ns*Np*sizeof(cplx));
+    for (iter = 0; iter < n_iter; iter++) {
+        build_basis(s, Ns, poles, Np, A);
+        int ncol = Nf*(Np+2) + Np;
+        int nrow = Ns*Nf;
+        /* real-stacked LS: (2*nrow) x ncol */
+        double *M = (double*) calloc((size_t)(2*nrow)*ncol, sizeof(double));
+        double *b = (double*) calloc((size_t)(2*nrow), sizeof(double));
+        for (k = 0; k < Nf; k++) {
+            for (r = 0; r < Ns; r++) {
+                int row = k*Ns + r;
+                cplx Fkr = F[(long)k*Ns+r];
+                for (j = 0; j < Np; j++) {
+                    cplx a = A[r*Np+j];
+                    M[(row)*ncol + k*(Np+2)+j]        = creal(a);
+                    M[(row+nrow)*ncol + k*(Np+2)+j]   = cimag(a);
+                    cplx neg = -Fkr*a;
+                    M[(row)*ncol + Nf*(Np+2)+j]       = creal(neg);
+                    M[(row+nrow)*ncol + Nf*(Np+2)+j]  = cimag(neg);
+                }
+                M[(row)*ncol + k*(Np+2)+Np]        = 1.0;   /* d (real) */
+                M[(row)*ncol + k*(Np+2)+Np+1]      = creal(s[r]);   /* e*s */
+                M[(row+nrow)*ncol + k*(Np+2)+Np+1] = cimag(s[r]);
+                b[row]      = creal(Fkr);
+                b[row+nrow] = cimag(Fkr);
+            }
+        }
+        double *x = (double*) calloc((size_t) ncol, sizeof(double));
+        lstsq_real(M, b, 2*nrow, ncol, x);
+        double *ctil = x + Nf*(Np+2);
+        cplx *cres = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        ctil_to_cres(ctil, poles, Np, cres);
+        /* relocate: roots of  D(s) + sum cres_i * D(s)/(s-a_i) */
+        cplx *D = (cplx*) malloc((size_t)(Np+1)*sizeof(cplx));
+        poly_from_roots(poles, Np, D);
+        cplx *numsig = (cplx*) malloc((size_t)(Np+1)*sizeof(cplx));
+        for (i = 0; i <= Np; i++) numsig[i] = D[i];
+        cplx *Di = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        cplx *sub = (cplx*) malloc((size_t) Np*sizeof(cplx));   /* poles minus i */
+        for (i = 0; i < Np; i++) {
+            int t = 0; for (j = 0; j < Np; j++) if (j != i) sub[t++] = poles[j];
+            poly_from_roots(sub, Np-1, Di);                     /* len Np, degree Np-1 */
+            for (j = 0; j < Np; j++) numsig[j+1] += cres[i]*Di[j]; /* align: Di is degree Np-1 (len Np), numsig degree Np (len Np+1) */
+        }
+        cplx *newp = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        poly_roots(numsig, Np, newp);
+        for (i = 0; i < Np; i++) if (creal(newp[i]) > 0) newp[i] = -creal(newp[i]) + I*cimag(newp[i]);
+        /* sort: real poles first, then by real, imag (keeps cc pairs adjacent-ish) */
+        for (i = 0; i < Np; i++) for (j = i+1; j < Np; j++) {
+            int swap = 0;
+            double ki = (fabs(cimag(newp[i]))>1e-6)?1:0, kj = (fabs(cimag(newp[j]))>1e-6)?1:0;
+            if (kj < ki) swap = 1;
+            else if (kj == ki) { if (creal(newp[j]) < creal(newp[i]) - 1e-30) swap = 1;
+                                 else if (fabs(creal(newp[j])-creal(newp[i]))<1e-30 && cimag(newp[j])<cimag(newp[i])) swap = 1; }
+            if (swap) { cplx t = newp[i]; newp[i]=newp[j]; newp[j]=t; }
+        }
+        for (i = 0; i < Np; i++) poles[i] = newp[i];
+        free(M); free(b); free(x); free(cres); free(D); free(numsig); free(Di); free(sub); free(newp);
+    }
+    /* final residues (fixed poles) */
+    build_basis(s, Ns, poles, Np, A);
+    for (k = 0; k < Nf; k++) {
+        int ncol = Np+2, nrow = Ns;
+        double *M = (double*) calloc((size_t)(2*nrow)*ncol, sizeof(double));
+        double *b = (double*) calloc((size_t)(2*nrow), sizeof(double));
+        for (r = 0; r < Ns; r++) {
+            cplx Fkr = F[(long)k*Ns+r];
+            for (j = 0; j < Np; j++) { M[r*ncol+j]=creal(A[r*Np+j]); M[(r+nrow)*ncol+j]=cimag(A[r*Np+j]); }
+            M[r*ncol+Np]=1.0;
+            M[r*ncol+Np+1]=creal(s[r]); M[(r+nrow)*ncol+Np+1]=cimag(s[r]);
+            b[r]=creal(Fkr); b[r+nrow]=cimag(Fkr);
+        }
+        double *x = (double*) calloc((size_t) ncol, sizeof(double));
+        lstsq_real(M, b, 2*nrow, ncol, x);
+        cplx *cr = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        ctil_to_cres(x, poles, Np, cr);
+        for (j = 0; j < Np; j++) res[(long)k*Np+j] = cr[j];
+        d[k] = x[Np]; e[k] = x[Np+1];
+        free(M); free(b); free(x); free(cr);
+    }
+    free(A);
+}
+
+static void seed_poles(double fmin, double fmax, int npair, cplx *p)
+{
+    int k;
+    for (k = 0; k < npair; k++) {
+        double beta = (npair==1) ? 2*M_PI*sqrt(fmin*fmax)
+                                 : 2*M_PI*fmin*pow(fmax/fmin, (double)k/(npair-1));
+        double a = beta/100.0;
+        p[2*k]   = -a + I*beta;
+        p[2*k+1] = -a - I*beta;
+    }
+}
+
+/* ============================ emit VA ============================ */
+/* strictly-proper numerator polynomial (ASCENDING, real) for Y_k = d + sum res/(s-p).
+ * num = d*D + sum_i res_i * D/(s-p_i). out[] holds Np+1 (degree Np). */
+static void num_proper(const cplx *poles, int Np, const cplx *res_k, double d_k, double *out_asc)
+{
+    int i, j, t;
+    cplx *D = (cplx*) malloc((size_t)(Np+1)*sizeof(cplx));
+    poly_from_roots(poles, Np, D);
+    cplx *num = (cplx*) calloc((size_t)(Np+1), sizeof(cplx));
+    for (i = 0; i <= Np; i++) num[i] = d_k*D[i];
+    cplx *Di = (cplx*) malloc((size_t) Np*sizeof(cplx));
+    cplx *sub = (cplx*) malloc((size_t) Np*sizeof(cplx));
+    for (i = 0; i < Np; i++) {
+        t = 0; for (j = 0; j < Np; j++) if (j != i) sub[t++] = poles[j];
+        poly_from_roots(sub, Np-1, Di);           /* degree Np-1, len Np */
+        for (j = 0; j < Np; j++) num[j+1] += res_k[i]*Di[j];
+    }
+    for (i = 0; i <= Np; i++) out_asc[Np-i] = creal(num[i]);   /* descending->ascending, real */
+    free(D); free(num); free(Di); free(sub);
+}
+
+static void emit_arr(FILE *f, const double *v, int n)
+{
+    int i; fprintf(f, "'{");
+    for (i = 0; i < n; i++) fprintf(f, "%s%.12g", i?", ":"", v[i]);
+    fprintf(f, "}");
+}
+
+/* ============================ public API ============================ */
+int snp2va_convert(const char *snpfile, const char *vafile, const char *module,
+                   char *msg, int msglen)
+{
+    TS ts; int i, j, k, r;
+    if (parse_touchstone(snpfile, &ts, msg, msglen)) return 1;
+    int N = ts.N, nf = ts.nf;
+    cplx *Y = (cplx*) malloc((size_t) nf*N*N*sizeof(cplx));
+    if (to_Y(&ts, Y)) { snprintf(msg,(size_t)msglen,"S->Y conversion failed (singular)"); ts_free(&ts); free(Y); return 1; }
+    /* Y entries in row-major i*N+j order, each a length-nf vector */
+    int Nf = N*N;
+    cplx *s = (cplx*) malloc((size_t) nf*sizeof(cplx));
+    for (r = 0; r < nf; r++) s[r] = 2*M_PI*ts.freqs[r]*I;
+    double wn = sqrt(cabs(s[0])*cabs(s[nf-1]));
+    cplx *sn = (cplx*) malloc((size_t) nf*sizeof(cplx));
+    for (r = 0; r < nf; r++) sn[r] = s[r]/wn;
+    cplx *F = (cplx*) malloc((size_t) Nf*nf*sizeof(cplx));
+    for (i = 0; i < N; i++) for (j = 0; j < N; j++) for (r = 0; r < nf; r++)
+        F[(long)(i*N+j)*nf+r] = Y[(long)r*N*N + i*N+j];
+
+    /* ---- order selection: climb, keep best STABLE fit, knee near a floor ---- */
+    double fmin = ts.freqs[0], fmax = ts.freqs[nf-1], tol = 1e-3;
+    cplx *bestP=NULL,*bestRes=NULL; double *bestD=NULL,*bestE=NULL; int bestNp=0; double bestErr=1e300;
+    cplx *prevP=NULL,*prevRes=NULL; double *prevD=NULL,*prevE=NULL; int prevNp=0; double prevErr=-1, firstErr=-1;
+    int chosenP=0; cplx *chP=NULL,*chRes=NULL; double *chD=NULL,*chE=NULL;
+    int npair;
+    for (npair = 1; npair <= 12; npair++) {
+        int Np = 2*npair;
+        cplx *P = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        cplx *res = (cplx*) malloc((size_t) Nf*Np*sizeof(cplx));
+        double *dd = (double*) malloc((size_t) Nf*sizeof(double));
+        double *ee = (double*) malloc((size_t) Nf*sizeof(double));
+        cplx *p0 = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        seed_poles(fmin, fmax, npair, p0);
+        for (i = 0; i < Np; i++) P[i] = p0[i]/wn;
+        vector_fit(sn, nf, F, Nf, Np, P, res, dd, ee, 12);
+        /* un-normalize */
+        for (i = 0; i < Np; i++) P[i] *= wn;
+        for (k = 0; k < Nf; k++) { for (i = 0; i < Np; i++) res[(long)k*Np+i] *= wn; ee[k] /= wn; }
+        /* rms rel error + stability */
+        double err = 0.0; int stable = 1;
+        for (i = 0; i < Np; i++) if (creal(P[i]) > 1e-6) stable = 0;
+        for (k = 0; k < Nf; k++) {
+            double numr=0, denr=0;
+            for (r = 0; r < nf; r++) {
+                cplx fit = dd[k] + s[r]*ee[k];
+                for (i = 0; i < Np; i++) fit += res[(long)k*Np+i]/(s[r]-P[i]);
+                cplx dif = fit - F[(long)k*nf+r];
+                numr += creal(dif)*creal(dif)+cimag(dif)*cimag(dif);
+                cplx fv = F[(long)k*nf+r];
+                denr += creal(fv)*creal(fv)+cimag(fv)*cimag(fv);
+            }
+            double e2 = sqrt(numr)/(sqrt(denr)+1e-300);
+            if (e2 > err) err = e2;
+        }
+        if (!(err==err)) stable = 0;   /* NaN */
+        free(p0);
+                if (firstErr < 0) firstErr = err;
+        int keep_best = stable && err < bestErr;
+        if (keep_best) {
+            free(bestP);free(bestRes);free(bestD);free(bestE);
+            bestP=P;bestRes=res;bestD=dd;bestE=ee;bestNp=Np;bestErr=err;
+        }
+        if (!stable) { if(!keep_best){free(P);free(res);free(dd);free(ee);} break; }
+        if (err < tol) { chosenP=Np; chP=P;chRes=res;chD=dd;chE=ee; if(keep_best){/*owned by best too*/} break; }
+        int near_floor = (err < 0.1*firstErr) || (err < 0.05);
+        if (prevErr >= 0 && err > 0.7*prevErr && near_floor) {  /* knee at floor -> use prev */
+            chosenP=prevNp; chP=prevP;chRes=prevRes;chD=prevD;chE=prevE;
+            if (!keep_best) { free(P);free(res);free(dd);free(ee); }
+            break;
+        }
+        /* shift prev <- current (free old prev unless it is the best) */
+        if (prevP && prevP!=bestP) { free(prevP);free(prevRes);free(prevD);free(prevE); }
+        prevP=P;prevRes=res;prevD=dd;prevE=ee;prevNp=Np;prevErr=err;
+    }
+    /* pick chosen, else best, else prev */
+    cplx *P; cplx *res; double *dd,*ee; int Np;
+    if (chP) { P=chP;res=chRes;dd=chD;ee=chE;Np=chosenP; }
+    else if (bestP) { P=bestP;res=bestRes;dd=bestD;ee=bestE;Np=bestNp; }
+    else { P=prevP;res=prevRes;dd=prevD;ee=prevE;Np=prevNp; }
+
+    /* ---- emit VA ---- */
+    FILE *fo = fopen(vafile, "w");
+    if (!fo) { snprintf(msg,(size_t)msglen,"cannot write '%s'", vafile); ts_free(&ts); return 1; }
+    double *den = (double*) malloc((size_t)(Np+1)*sizeof(double));
+    { cplx *D = (cplx*) malloc((size_t)(Np+1)*sizeof(cplx)); poly_from_roots(P, Np, D);
+      for (i = 0; i <= Np; i++) den[Np-i] = creal(D[i]); free(D); }
+    fprintf(fo, "`include \"disciplines.vams\"\n\n");
+    fprintf(fo, "// Generated by pre_snp from %s\n", snpfile);
+    fprintf(fo, "// %d-port, %d common poles; realized with laplace_nd (AC + transient).\n", N, Np);
+    fprintf(fo, "module %s(", module);
+    for (i = 0; i < N; i++) fprintf(fo, "%sp%d", i?", ":"", i+1);
+    fprintf(fo, ");\n    inout ");
+    for (i = 0; i < N; i++) fprintf(fo, "%sp%d", i?", ":"", i+1);
+    fprintf(fo, ";\n    electrical ");
+    for (i = 0; i < N; i++) fprintf(fo, "%sp%d", i?", ":"", i+1);
+    fprintf(fo, ";\n    analog begin\n");
+    double *nm = (double*) malloc((size_t)(Np+1)*sizeof(double));
+    for (i = 0; i < N; i++) {
+        fprintf(fo, "        I(p%d) <+ ", i+1);
+        for (j = 0; j < N; j++) {
+            int idx = i*N+j;
+            num_proper(P, Np, res + (long) idx*Np, dd[idx], nm);
+            if (j) fprintf(fo, "\n                 + ");
+            fprintf(fo, "laplace_nd(V(p%d), ", j+1); emit_arr(fo, nm, Np+1);
+            fprintf(fo, ", "); emit_arr(fo, den, Np+1); fprintf(fo, ")");
+            if (fabs(ee[idx]) > 1e-30) fprintf(fo, "\n                 + (%.12g)*ddt(V(p%d))", ee[idx], j+1);
+        }
+        fprintf(fo, ";\n");
+    }
+    fprintf(fo, "    end\nendmodule\n");
+    fclose(fo);
+    snprintf(msg,(size_t)msglen,"%d-port, %d poles, rms rel err %.2e", N, Np, bestErr<1e300?bestErr:0.0);
+    /* frees (leak-tolerant: one-shot tool) */
+    free(den); free(nm); free(Y); free(s); free(sn); free(F); ts_free(&ts);
+    return 0;
+}
+
+#ifdef SNP2VA_TEST
+int main(int argc, char **argv)
+{
+    setbuf(stderr, NULL);
+    if (argc < 3) { fprintf(stderr, "usage: %s file.sNp out.va [module]\n", argv[0]); return 2; }
+    char msg[256];
+    int rc = snp2va_convert(argv[1], argv[2], argc>3?argv[3]:"nport", msg, sizeof msg);
+    fprintf(stderr, "snp2va: %s\n", msg);
+    return rc;
+}
+#endif
diff --git a/src/frontend/snp2va.cpp b/src/frontend/snp2va.cpp
new file mode 100644
index 000000000..4404f769b
--- /dev/null
+++ b/src/frontend/snp2va.cpp
@@ -0,0 +1,566 @@
+/* Enhancement-200: Touchstone (.sNp) -> Verilog-A n-port converter, in C.
+ *
+ * A C port of the pure-Python snp2va.py (Enhancement-199): parse a Touchstone
+ * S-parameter file, convert S -> Y, fit every Y_ij(f) with a common-pole rational
+ * (Gustavsen vector fitting), and emit a Verilog-A n-port realized with laplace_nd,
+ * so through OpenVAF/OSDI it works in AC and transient. Used by the `pre_snp`
+ * front-end command, which then invokes openvaf-r to compile the emitted .va.
+ *
+ * The numerical core is self-contained (only stdio/stdlib/string/math/complex),
+ * so it does not pull in ngspice's own complex.h. Public entry point:
+ *     int snp2va_convert(const char *snp, const char *va, const char *module,
+ *                        char *msg, int msglen);
+ * returns 0 on success, non-zero on failure (msg gets a one-line status).
+ */
+#ifdef _MSC_VER
+#include <iostream>
+#include <complex>
+
+extern "C" {
+#include "snp2va.h"
+#include <string.h>
+};
+#else
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include <complex.h>
+#include "snp2va.h"
+#endif
+
+#ifdef _MSC_VER
+typedef std::complex<double> cplx;
+#define cabs abs
+#define cpow pow
+#define cimag imag
+#define creal real
+#define cexp exp
+#define I cplx(0.0, 1.0)
+#define strcasecmp _stricmp
+#else
+typedef double _Complex cplx;
+#endif
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+/* ============================ linear algebra ============================ */
+
+/* Least-squares min||A x - b|| for a REAL overdetermined system (m>=n) via
+ * Householder QR. A is row-major m*n, b length m, x length n. Returns 0 on ok. */
+static int lstsq_real(double *A, double *b, int m, int n, double *x)
+{
+    int i, j, k;
+    for (k = 0; k < n; k++) {
+        double norm = 0.0;
+        for (i = k; i < m; i++) norm += A[i*n+k]*A[i*n+k];
+        norm = sqrt(norm);
+        if (norm == 0.0) continue;
+        double alpha = (A[k*n+k] >= 0.0) ? -norm : norm;
+        double *v = (double*) calloc((size_t) m, sizeof(double));
+        v[k] = A[k*n+k] - alpha;
+        for (i = k+1; i < m; i++) v[i] = A[i*n+k];
+        double vn2 = 0.0;
+        for (i = k; i < m; i++) vn2 += v[i]*v[i];
+        if (vn2 == 0.0) { free(v); continue; }
+        for (j = k; j < n; j++) {
+            double s = 0.0;
+            for (i = k; i < m; i++) s += v[i]*A[i*n+j];
+            s = s*2.0/vn2;
+            for (i = k; i < m; i++) A[i*n+j] -= s*v[i];
+        }
+        double s = 0.0;
+        for (i = k; i < m; i++) s += v[i]*b[i];
+        s = s*2.0/vn2;
+        for (i = k; i < m; i++) b[i] -= s*v[i];
+        free(v);
+    }
+    for (i = n-1; i >= 0; i--) {
+        double acc = b[i];
+        for (j = i+1; j < n; j++) acc -= A[i*n+j]*x[j];
+        x[i] = (A[i*n+i] != 0.0) ? acc/A[i*n+i] : 0.0;
+    }
+    return 0;
+}
+
+/* In-place inverse of an n x n complex matrix (row-major), Gauss-Jordan w/ pivot.
+ * Returns 0 on ok, 1 if singular. */
+static int mat_inv_c(cplx *M, int n)
+{
+    int i, j, c, p;
+    cplx *A = (cplx*) malloc((size_t) n*2*n*sizeof(cplx));
+    for (i = 0; i < n; i++) {
+        for (j = 0; j < n; j++) A[i*2*n+j] = M[i*n+j];
+        for (j = 0; j < n; j++) A[i*2*n+n+j] = (i==j) ? 1.0 : 0.0;
+    }
+    for (c = 0; c < n; c++) {
+        p = c; double best = cabs(A[c*2*n+c]);
+        for (i = c+1; i < n; i++) { double v = cabs(A[i*2*n+c]); if (v > best) { best = v; p = i; } }
+        if (best == 0.0) { free(A); return 1; }
+        if (p != c) for (j = 0; j < 2*n; j++) { cplx t = A[c*2*n+j]; A[c*2*n+j]=A[p*2*n+j]; A[p*2*n+j]=t; }
+        cplx d = A[c*2*n+c];
+        for (j = 0; j < 2*n; j++) A[c*2*n+j] /= d;
+        for (i = 0; i < n; i++) if (i != c) {
+            cplx f = A[i*2*n+c];
+            if (f != 0.0) for (j = 0; j < 2*n; j++) A[i*2*n+j] -= f*A[c*2*n+j];
+        }
+    }
+    for (i = 0; i < n; i++) for (j = 0; j < n; j++) M[i*n+j] = A[i*2*n+n+j];
+    free(A);
+    return 0;
+}
+
+/* Monic polynomial (DESCENDING, len nr+1) from roots. coef must hold nr+1. */
+static void poly_from_roots(const cplx *roots, int nr, cplx *coef)
+{
+    int i, k;
+    coef[0] = 1.0;
+    for (i = 1; i <= nr; i++) coef[i] = 0.0;
+    for (k = 0; k < nr; k++) {
+        for (i = k+1; i >= 1; i--) coef[i] = coef[i] - roots[k]*coef[i-1];
+    }
+}
+
+static cplx poly_eval(const cplx *c, int deg, cplx x)
+{
+    cplx r = 0.0; int i;
+    for (i = 0; i <= deg; i++) r = r*x + c[i];
+    return r;
+}
+
+/* All roots of a DESCENDING-coeff polynomial (deg = len-1) via Durand-Kerner.
+ * roots[] must hold deg. Returns 0 on ok. */
+static int poly_roots(const cplx *cin, int deg, cplx *roots)
+{
+    int i, j, it;
+    if (deg <= 0) return 0;
+    cplx *c = (cplx*) malloc((size_t)(deg+1)*sizeof(cplx));
+    for (i = 0; i <= deg; i++) c[i] = cin[i] / cin[0];   /* monic */
+    cplx seed = 0.4 + 0.9*I;
+    for (i = 0; i < deg; i++) roots[i] = cpow(seed, (double) i);
+    for (it = 0; it < 500; it++) {
+        double maxstep = 0.0;
+        for (i = 0; i < deg; i++) {
+            cplx num = poly_eval(c, deg, roots[i]);
+            cplx den = 1.0;
+            for (j = 0; j < deg; j++) if (j != i) den *= (roots[i]-roots[j]);
+            cplx step = (cabs(den) > 1e-300) ? num/den : 0.0;
+            roots[i] -= step;
+            if (cabs(step) > maxstep) maxstep = cabs(step);
+        }
+        if (maxstep < 1e-14) break;
+    }
+    free(c);
+    return 0;
+}
+
+/* ============================ Touchstone I/O ============================ */
+
+typedef struct { double *freqs; cplx *S; int nf; int N; double z0; char ptype; } TS;
+
+static void ts_free(TS *t) { free(t->freqs); free(t->S); }
+
+/* returns 0 on ok */
+static int parse_touchstone(const char *fn, TS *out, char *msg, int msglen)
+{
+    FILE *f = fopen(fn, "r");
+    if (!f) { snprintf(msg, (size_t) msglen, "cannot open '%s'", fn); return 1; }
+    double fmul = 1e9, z0 = 50.0; char ptype = 'S'; char fmt[3] = "MA";
+    /* collect all numeric tokens after the '#' options line(s) */
+    double *nums = NULL; long ncap = 0, nn = 0;
+    char line[4096];
+    while (fgets(line, sizeof line, f)) {
+        char *h = strchr(line, '!'); if (h) *h = '\0';
+        char *p = line;
+        while (*p && isspace((unsigned char)*p)) p++;
+        if (*p == '\0') continue;
+        if (*p == '#') {
+            char *tok = strtok(p+1, " \t\r\n");
+            while (tok) {
+                if (!strcasecmp(tok,"HZ")) fmul=1.0;
+                else if (!strcasecmp(tok,"KHZ")) fmul=1e3;
+                else if (!strcasecmp(tok,"MHZ")) fmul=1e6;
+                else if (!strcasecmp(tok,"GHZ")) fmul=1e9;
+                else if (!strcasecmp(tok,"S")||!strcasecmp(tok,"Y")||!strcasecmp(tok,"Z")) ptype=(char)toupper((unsigned char)tok[0]);
+                else if (!strcasecmp(tok,"MA")||!strcasecmp(tok,"DB")||!strcasecmp(tok,"RI")) { fmt[0]=(char)toupper((unsigned char)tok[0]); fmt[1]=(char)toupper((unsigned char)tok[1]); }
+                else if (!strcasecmp(tok,"R")) { char *z=strtok(NULL," \t\r\n"); if (z) z0=atof(z); }
+                tok = strtok(NULL, " \t\r\n");
+            }
+            continue;
+        }
+        /* numeric data line */
+        char *tok = strtok(p, " \t\r\n");
+        while (tok) {
+            char *end; double v = strtod(tok, &end);
+            if (end != tok) {
+                if (nn >= ncap) { ncap = ncap ? ncap*2 : 1024; nums = (double*) realloc(nums, (size_t) ncap*sizeof(double)); }
+                nums[nn++] = v;
+            }
+            tok = strtok(NULL, " \t\r\n");
+        }
+    }
+    fclose(f);
+    /* infer port count N: try the file extension .sNp, else brute force */
+    int N = 0;
+    const char *dot = strrchr(fn, '.');
+    if (dot && (dot[1]=='s'||dot[1]=='S') && (fn[strlen(fn)-1]=='p'||fn[strlen(fn)-1]=='P')) {
+        N = atoi(dot+2);
+    }
+    if (N <= 0) {
+        int c;
+        for (c = 1; c <= 16; c++) if (nn % (1 + 2*c*c) == 0) { N = c; break; }
+    }
+    if (N <= 0) { free(nums); snprintf(msg,(size_t)msglen,"cannot determine port count"); return 1; }
+    int rec = 1 + 2*N*N;
+    int nf = (int)(nn / rec);
+    if (nf < 2) { free(nums); snprintf(msg,(size_t)msglen,"too few frequency points (%d)", nf); return 1; }
+    out->freqs = (double*) malloc((size_t) nf*sizeof(double));
+    out->S = (cplx*) malloc((size_t) nf*N*N*sizeof(cplx));
+    out->nf = nf; out->N = N; out->z0 = z0; out->ptype = ptype;
+    int r, kk;
+    for (r = 0; r < nf; r++) {
+        double *chunk = nums + (long) r*rec;
+        out->freqs[r] = chunk[0]*fmul;
+        double *vals = chunk+1;
+        cplx pv[16*16];
+        for (kk = 0; kk < N*N; kk++) {
+            double a = vals[2*kk], b = vals[2*kk+1];
+            if (!strcmp(fmt,"MA")) pv[kk] = a*cexp(I*b*M_PI/180.0);
+            else if (!strcmp(fmt,"DB")) pv[kk] = pow(10.0,a/20.0)*cexp(I*b*M_PI/180.0);
+            else pv[kk] = a + I*b;
+        }
+        /* Touchstone: N=2 order is S11 S21 S12 S22; general is row-major */
+        cplx *M = out->S + (long) r*N*N;
+        if (N == 2) { M[0]=pv[0]; M[2]=pv[1]; M[1]=pv[2]; M[3]=pv[3]; }
+        else for (kk = 0; kk < N*N; kk++) M[kk] = pv[kk];
+    }
+    free(nums);
+    return 0;
+}
+
+/* S/Y/Z -> Y (row-major per frequency), Yout must hold nf*N*N. */
+static int to_Y(const TS *t, cplx *Yout)
+{
+    int N = t->N, r, i, j;
+    cplx *tmp = (cplx*) malloc((size_t) N*N*sizeof(cplx));
+    for (r = 0; r < t->nf; r++) {
+        const cplx *M = t->S + (long) r*N*N;
+        cplx *Y = Yout + (long) r*N*N;
+        if (t->ptype == 'Y') { for (i=0;i<N*N;i++) Y[i]=M[i]; }
+        else if (t->ptype == 'Z') { for (i=0;i<N*N;i++) Y[i]=M[i]; if (mat_inv_c(Y,N)) { free(tmp); return 1; } }
+        else { /* S -> Y = (1/z0)(I-S)(I+S)^-1 */
+            cplx *IpS = tmp;
+            for (i=0;i<N;i++) for (j=0;j<N;j++) IpS[i*N+j] = ((i==j)?1.0:0.0) + M[i*N+j];
+            if (mat_inv_c(IpS,N)) { free(tmp); return 1; }
+            for (i=0;i<N;i++) for (j=0;j<N;j++) {
+                cplx acc = 0.0; int k;
+                for (k=0;k<N;k++) acc += (((i==k)?1.0:0.0) - M[i*N+k]) * IpS[k*N+j];
+                Y[i*N+j] = acc / t->z0;
+            }
+        }
+    }
+    free(tmp);
+    return 0;
+}
+
+/* ============================ vector fitting ============================ */
+/* layout: for each pole, 0='real', 1='cc-start' (its conjugate is the next). */
+static int build_layout(const cplx *poles, int Np, int *lay)
+{
+    int i = 0, m = 0;
+    while (i < Np) {
+        if (fabs(cimag(poles[i])) < 1e-9*fabs(creal(poles[i]))+1e-30) { lay[m++] = 0; i += 1; }
+        else { lay[m++] = 1; i += 2; }
+    }
+    return m;   /* number of blocks */
+}
+
+/* complex partial-fraction basis (real-valued cc combos), Ns x Np, row-major */
+static void build_basis(const cplx *s, int Ns, const cplx *poles, int Np, cplx *A)
+{
+    int r, i;
+    for (r = 0; r < Ns; r++) {
+        i = 0;
+        while (i < Np) {
+            if (fabs(cimag(poles[i])) < 1e-9*fabs(creal(poles[i]))+1e-30) {
+                A[r*Np+i] = 1.0/(s[r]-poles[i]); i += 1;
+            } else {
+                cplx p = poles[i];
+                A[r*Np+i]   = 1.0/(s[r]-p) + 1.0/(s[r]-conj(p));
+                A[r*Np+i+1] = I/(s[r]-p) - I/(s[r]-conj(p));
+                i += 2;
+            }
+        }
+    }
+}
+
+/* complex residues from real ctil coeffs, per real/cc layout */
+static void ctil_to_cres(const double *ctil, const cplx *poles, int Np, cplx *cres)
+{
+    int i = 0;
+    while (i < Np) {
+        if (fabs(cimag(poles[i])) < 1e-9*fabs(creal(poles[i]))+1e-30) { cres[i] = ctil[i]; i += 1; }
+        else { cres[i] = ctil[i] + I*ctil[i+1]; cres[i+1] = conj(cres[i]); i += 2; }
+    }
+}
+
+/* One vector-fit run (fixed pole count). s,F normalized. Returns fit in poles/res/d/e. */
+static void vector_fit(const cplx *s, int Ns, const cplx *F, int Nf, int Np,
+                       cplx *poles, cplx *res, double *d, double *e, int n_iter)
+{
+    int iter, i, j, k, r;
+    cplx *A = (cplx*) malloc((size_t) Ns*Np*sizeof(cplx));
+    for (iter = 0; iter < n_iter; iter++) {
+        build_basis(s, Ns, poles, Np, A);
+        int ncol = Nf*(Np+2) + Np;
+        int nrow = Ns*Nf;
+        /* real-stacked LS: (2*nrow) x ncol */
+        double *M = (double*) calloc((size_t)(2*nrow)*ncol, sizeof(double));
+        double *b = (double*) calloc((size_t)(2*nrow), sizeof(double));
+        for (k = 0; k < Nf; k++) {
+            for (r = 0; r < Ns; r++) {
+                int row = k*Ns + r;
+                cplx Fkr = F[(long)k*Ns+r];
+                for (j = 0; j < Np; j++) {
+                    cplx a = A[r*Np+j];
+                    M[(row)*ncol + k*(Np+2)+j]        = creal(a);
+                    M[(row+nrow)*ncol + k*(Np+2)+j]   = cimag(a);
+                    cplx neg = -Fkr*a;
+                    M[(row)*ncol + Nf*(Np+2)+j]       = creal(neg);
+                    M[(row+nrow)*ncol + Nf*(Np+2)+j]  = cimag(neg);
+                }
+                M[(row)*ncol + k*(Np+2)+Np]        = 1.0;   /* d (real) */
+                M[(row)*ncol + k*(Np+2)+Np+1]      = creal(s[r]);   /* e*s */
+                M[(row+nrow)*ncol + k*(Np+2)+Np+1] = cimag(s[r]);
+                b[row]      = creal(Fkr);
+                b[row+nrow] = cimag(Fkr);
+            }
+        }
+        double *x = (double*) calloc((size_t) ncol, sizeof(double));
+        lstsq_real(M, b, 2*nrow, ncol, x);
+        double *ctil = x + Nf*(Np+2);
+        cplx *cres = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        ctil_to_cres(ctil, poles, Np, cres);
+        /* relocate: roots of  D(s) + sum cres_i * D(s)/(s-a_i) */
+        cplx *D = (cplx*) malloc((size_t)(Np+1)*sizeof(cplx));
+        poly_from_roots(poles, Np, D);
+        cplx *numsig = (cplx*) malloc((size_t)(Np+1)*sizeof(cplx));
+        for (i = 0; i <= Np; i++) numsig[i] = D[i];
+        cplx *Di = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        cplx *sub = (cplx*) malloc((size_t) Np*sizeof(cplx));   /* poles minus i */
+        for (i = 0; i < Np; i++) {
+            int t = 0; for (j = 0; j < Np; j++) if (j != i) sub[t++] = poles[j];
+            poly_from_roots(sub, Np-1, Di);                     /* len Np, degree Np-1 */
+            for (j = 0; j < Np; j++) numsig[j+1] += cres[i]*Di[j]; /* align: Di is degree Np-1 (len Np), numsig degree Np (len Np+1) */
+        }
+        cplx *newp = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        poly_roots(numsig, Np, newp);
+        for (i = 0; i < Np; i++) if (creal(newp[i]) > 0) newp[i] = -creal(newp[i]) + I*cimag(newp[i]);
+        /* sort: real poles first, then by real, imag (keeps cc pairs adjacent-ish) */
+        for (i = 0; i < Np; i++) for (j = i+1; j < Np; j++) {
+            int swap = 0;
+            double ki = (fabs(cimag(newp[i]))>1e-6)?1:0, kj = (fabs(cimag(newp[j]))>1e-6)?1:0;
+            if (kj < ki) swap = 1;
+            else if (kj == ki) { if (creal(newp[j]) < creal(newp[i]) - 1e-30) swap = 1;
+                                 else if (fabs(creal(newp[j])-creal(newp[i]))<1e-30 && cimag(newp[j])<cimag(newp[i])) swap = 1; }
+            if (swap) { cplx t = newp[i]; newp[i]=newp[j]; newp[j]=t; }
+        }
+        for (i = 0; i < Np; i++) poles[i] = newp[i];
+        free(M); free(b); free(x); free(cres); free(D); free(numsig); free(Di); free(sub); free(newp);
+    }
+    /* final residues (fixed poles) */
+    build_basis(s, Ns, poles, Np, A);
+    for (k = 0; k < Nf; k++) {
+        int ncol = Np+2, nrow = Ns;
+        double *M = (double*) calloc((size_t)(2*nrow)*ncol, sizeof(double));
+        double *b = (double*) calloc((size_t)(2*nrow), sizeof(double));
+        for (r = 0; r < Ns; r++) {
+            cplx Fkr = F[(long)k*Ns+r];
+            for (j = 0; j < Np; j++) { M[r*ncol+j]=creal(A[r*Np+j]); M[(r+nrow)*ncol+j]=cimag(A[r*Np+j]); }
+            M[r*ncol+Np]=1.0;
+            M[r*ncol+Np+1]=creal(s[r]); M[(r+nrow)*ncol+Np+1]=cimag(s[r]);
+            b[r]=creal(Fkr); b[r+nrow]=cimag(Fkr);
+        }
+        double *x = (double*) calloc((size_t) ncol, sizeof(double));
+        lstsq_real(M, b, 2*nrow, ncol, x);
+        cplx *cr = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        ctil_to_cres(x, poles, Np, cr);
+        for (j = 0; j < Np; j++) res[(long)k*Np+j] = cr[j];
+        d[k] = x[Np]; e[k] = x[Np+1];
+        free(M); free(b); free(x); free(cr);
+    }
+    free(A);
+}
+
+static void seed_poles(double fmin, double fmax, int npair, cplx *p)
+{
+    int k;
+    for (k = 0; k < npair; k++) {
+        double beta = (npair==1) ? 2*M_PI*sqrt(fmin*fmax)
+                                 : 2*M_PI*fmin*pow(fmax/fmin, (double)k/(npair-1));
+        double a = beta/100.0;
+        p[2*k]   = -a + I*beta;
+        p[2*k+1] = -a - I*beta;
+    }
+}
+
+/* ============================ emit VA ============================ */
+/* strictly-proper numerator polynomial (ASCENDING, real) for Y_k = d + sum res/(s-p).
+ * num = d*D + sum_i res_i * D/(s-p_i). out[] holds Np+1 (degree Np). */
+static void num_proper(const cplx *poles, int Np, const cplx *res_k, double d_k, double *out_asc)
+{
+    int i, j, t;
+    cplx *D = (cplx*) malloc((size_t)(Np+1)*sizeof(cplx));
+    poly_from_roots(poles, Np, D);
+    cplx *num = (cplx*) calloc((size_t)(Np+1), sizeof(cplx));
+    for (i = 0; i <= Np; i++) num[i] = d_k*D[i];
+    cplx *Di = (cplx*) malloc((size_t) Np*sizeof(cplx));
+    cplx *sub = (cplx*) malloc((size_t) Np*sizeof(cplx));
+    for (i = 0; i < Np; i++) {
+        t = 0; for (j = 0; j < Np; j++) if (j != i) sub[t++] = poles[j];
+        poly_from_roots(sub, Np-1, Di);           /* degree Np-1, len Np */
+        for (j = 0; j < Np; j++) num[j+1] += res_k[i]*Di[j];
+    }
+    for (i = 0; i <= Np; i++) out_asc[Np-i] = creal(num[i]);   /* descending->ascending, real */
+    free(D); free(num); free(Di); free(sub);
+}
+
+static void emit_arr(FILE *f, const double *v, int n)
+{
+    int i; fprintf(f, "'{");
+    for (i = 0; i < n; i++) fprintf(f, "%s%.12g", i?", ":"", v[i]);
+    fprintf(f, "}");
+}
+
+/* ============================ public API ============================ */
+int snp2va_convert(const char *snpfile, const char *vafile, const char *module,
+                   char *msg, int msglen)
+{
+    TS ts; int i, j, k, r;
+    if (parse_touchstone(snpfile, &ts, msg, msglen)) return 1;
+    int N = ts.N, nf = ts.nf;
+    cplx *Y = (cplx*) malloc((size_t) nf*N*N*sizeof(cplx));
+    if (to_Y(&ts, Y)) { snprintf(msg,(size_t)msglen,"S->Y conversion failed (singular)"); ts_free(&ts); free(Y); return 1; }
+    /* Y entries in row-major i*N+j order, each a length-nf vector */
+    int Nf = N*N;
+    cplx *s = (cplx*) malloc((size_t) nf*sizeof(cplx));
+    for (r = 0; r < nf; r++) s[r] = 2*M_PI*ts.freqs[r]*I;
+    double wn = sqrt(cabs(s[0])*cabs(s[nf-1]));
+    cplx *sn = (cplx*) malloc((size_t) nf*sizeof(cplx));
+    for (r = 0; r < nf; r++) sn[r] = s[r]/wn;
+    cplx *F = (cplx*) malloc((size_t) Nf*nf*sizeof(cplx));
+    for (i = 0; i < N; i++) for (j = 0; j < N; j++) for (r = 0; r < nf; r++)
+        F[(long)(i*N+j)*nf+r] = Y[(long)r*N*N + i*N+j];
+
+    /* ---- order selection: climb, keep best STABLE fit, knee near a floor ---- */
+    double fmin = ts.freqs[0], fmax = ts.freqs[nf-1], tol = 1e-3;
+    cplx *bestP=NULL,*bestRes=NULL; double *bestD=NULL,*bestE=NULL; int bestNp=0; double bestErr=1e300;
+    cplx *prevP=NULL,*prevRes=NULL; double *prevD=NULL,*prevE=NULL; int prevNp=0; double prevErr=-1, firstErr=-1;
+    int chosenP=0; cplx *chP=NULL,*chRes=NULL; double *chD=NULL,*chE=NULL;
+    int npair;
+    for (npair = 1; npair <= 12; npair++) {
+        int Np = 2*npair;
+        cplx *P = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        cplx *res = (cplx*) malloc((size_t) Nf*Np*sizeof(cplx));
+        double *dd = (double*) malloc((size_t) Nf*sizeof(double));
+        double *ee = (double*) malloc((size_t) Nf*sizeof(double));
+        cplx *p0 = (cplx*) malloc((size_t) Np*sizeof(cplx));
+        seed_poles(fmin, fmax, npair, p0);
+        for (i = 0; i < Np; i++) P[i] = p0[i]/wn;
+        vector_fit(sn, nf, F, Nf, Np, P, res, dd, ee, 12);
+        /* un-normalize */
+        for (i = 0; i < Np; i++) P[i] *= wn;
+        for (k = 0; k < Nf; k++) { for (i = 0; i < Np; i++) res[(long)k*Np+i] *= wn; ee[k] /= wn; }
+        /* rms rel error + stability */
+        double err = 0.0; int stable = 1;
+        for (i = 0; i < Np; i++) if (creal(P[i]) > 1e-6) stable = 0;
+        for (k = 0; k < Nf; k++) {
+            double numr=0, denr=0;
+            for (r = 0; r < nf; r++) {
+                cplx fit = dd[k] + s[r]*ee[k];
+                for (i = 0; i < Np; i++) fit += res[(long)k*Np+i]/(s[r]-P[i]);
+                cplx dif = fit - F[(long)k*nf+r];
+                numr += creal(dif)*creal(dif)+cimag(dif)*cimag(dif);
+                cplx fv = F[(long)k*nf+r];
+                denr += creal(fv)*creal(fv)+cimag(fv)*cimag(fv);
+            }
+            double e2 = sqrt(numr)/(sqrt(denr)+1e-300);
+            if (e2 > err) err = e2;
+        }
+        if (!(err==err)) stable = 0;   /* NaN */
+        free(p0);
+                if (firstErr < 0) firstErr = err;
+        int keep_best = stable && err < bestErr;
+        if (keep_best) {
+            free(bestP);free(bestRes);free(bestD);free(bestE);
+            bestP=P;bestRes=res;bestD=dd;bestE=ee;bestNp=Np;bestErr=err;
+        }
+        if (!stable) { if(!keep_best){free(P);free(res);free(dd);free(ee);} break; }
+        if (err < tol) { chosenP=Np; chP=P;chRes=res;chD=dd;chE=ee; if(keep_best){/*owned by best too*/} break; }
+        int near_floor = (err < 0.1*firstErr) || (err < 0.05);
+        if (prevErr >= 0 && err > 0.7*prevErr && near_floor) {  /* knee at floor -> use prev */
+            chosenP=prevNp; chP=prevP;chRes=prevRes;chD=prevD;chE=prevE;
+            if (!keep_best) { free(P);free(res);free(dd);free(ee); }
+            break;
+        }
+        /* shift prev <- current (free old prev unless it is the best) */
+        if (prevP && prevP!=bestP) { free(prevP);free(prevRes);free(prevD);free(prevE); }
+        prevP=P;prevRes=res;prevD=dd;prevE=ee;prevNp=Np;prevErr=err;
+    }
+    /* pick chosen, else best, else prev */
+    cplx *P; cplx *res; double *dd,*ee; int Np;
+    if (chP) { P=chP;res=chRes;dd=chD;ee=chE;Np=chosenP; }
+    else if (bestP) { P=bestP;res=bestRes;dd=bestD;ee=bestE;Np=bestNp; }
+    else { P=prevP;res=prevRes;dd=prevD;ee=prevE;Np=prevNp; }
+
+    /* ---- emit VA ---- */
+    FILE *fo = fopen(vafile, "w");
+    if (!fo) { snprintf(msg,(size_t)msglen,"cannot write '%s'", vafile); ts_free(&ts); return 1; }
+    double *den = (double*) malloc((size_t)(Np+1)*sizeof(double));
+    { cplx *D = (cplx*) malloc((size_t)(Np+1)*sizeof(cplx)); poly_from_roots(P, Np, D);
+      for (i = 0; i <= Np; i++) den[Np-i] = creal(D[i]); free(D); }
+    fprintf(fo, "`include \"disciplines.vams\"\n\n");
+    fprintf(fo, "// Generated by pre_snp from %s\n", snpfile);
+    fprintf(fo, "// %d-port, %d common poles; realized with laplace_nd (AC + transient).\n", N, Np);
+    fprintf(fo, "module %s(", module);
+    for (i = 0; i < N; i++) fprintf(fo, "%sp%d", i?", ":"", i+1);
+    fprintf(fo, ");\n    inout ");
+    for (i = 0; i < N; i++) fprintf(fo, "%sp%d", i?", ":"", i+1);
+    fprintf(fo, ";\n    electrical ");
+    for (i = 0; i < N; i++) fprintf(fo, "%sp%d", i?", ":"", i+1);
+    fprintf(fo, ";\n    analog begin\n");
+    double *nm = (double*) malloc((size_t)(Np+1)*sizeof(double));
+    for (i = 0; i < N; i++) {
+        fprintf(fo, "        I(p%d) <+ ", i+1);
+        for (j = 0; j < N; j++) {
+            int idx = i*N+j;
+            num_proper(P, Np, res + (long) idx*Np, dd[idx], nm);
+            if (j) fprintf(fo, "\n                 + ");
+            fprintf(fo, "laplace_nd(V(p%d), ", j+1); emit_arr(fo, nm, Np+1);
+            fprintf(fo, ", "); emit_arr(fo, den, Np+1); fprintf(fo, ")");
+            if (fabs(ee[idx]) > 1e-30) fprintf(fo, "\n                 + (%.12g)*ddt(V(p%d))", ee[idx], j+1);
+        }
+        fprintf(fo, ";\n");
+    }
+    fprintf(fo, "    end\nendmodule\n");
+    fclose(fo);
+    snprintf(msg,(size_t)msglen,"%d-port, %d poles, rms rel err %.2e", N, Np, bestErr<1e300?bestErr:0.0);
+    /* frees (leak-tolerant: one-shot tool) */
+    free(den); free(nm); free(Y); free(s); free(sn); free(F); ts_free(&ts);
+    return 0;
+}
+
+#ifdef SNP2VA_TEST
+int main(int argc, char **argv)
+{
+    setbuf(stderr, NULL);
+    if (argc < 3) { fprintf(stderr, "usage: %s file.sNp out.va [module]\n", argv[0]); return 2; }
+    char msg[256];
+    int rc = snp2va_convert(argv[1], argv[2], argc>3?argv[3]:"nport", msg, sizeof msg);
+    fprintf(stderr, "snp2va: %s\n", msg);
+    return rc;
+}
+#endif
diff --git a/src/frontend/snp2va.h b/src/frontend/snp2va.h
new file mode 100644
index 000000000..010d58bba
--- /dev/null
+++ b/src/frontend/snp2va.h
@@ -0,0 +1,10 @@
+/* Enhancement-200: Touchstone (.sNp) -> Verilog-A converter + `pre_snp` command. */
+#ifndef SNP2VA_H
+#define SNP2VA_H
+#include "ngspice/wordlist.h"
+/* Convert a Touchstone file to a Verilog-A n-port model. Returns 0 on success;
+ * msg gets a one-line status/error. (No ngspice deps in the converter core.) */
+int snp2va_convert(const char *snpfile, const char *vafile, const char *module,
+                   char *msg, int msglen);
+void com_pre_snp(wordlist *wl);
+#endif
diff --git a/visualc/vngspice.vcxproj b/visualc/vngspice.vcxproj
index 7447299a7..ecb7bcf7f 100644
--- a/visualc/vngspice.vcxproj
+++ b/visualc/vngspice.vcxproj
@@ -972,6 +972,7 @@
     <ClInclude Include="..\src\frontend\shyu.h" />
     <ClInclude Include="..\src\frontend\signal_handler.h" />
     <ClInclude Include="..\src\frontend\sndprint.h" />
+    <ClInclude Include="..\src\frontend\snp2va.h" />
     <ClInclude Include="..\src\frontend\spec.h" />
     <ClInclude Include="..\src\frontend\spiceif.h" />
     <ClInclude Include="..\src\frontend\streams.h" />
@@ -1519,6 +1520,7 @@
     <ClCompile Include="..\src\frontend\com_optimize.c" />
     <ClCompile Include="..\src\frontend\com_option.c" />
     <ClCompile Include="..\src\frontend\com_plot.c" />
+    <ClCompile Include="..\src\frontend\com_presnp.c" />
     <ClCompile Include="..\src\frontend\com_pyplot.c" />
     <ClCompile Include="..\src\frontend\com_qpac.c" />
     <ClCompile Include="..\src\frontend\com_qpnoise.c" />
@@ -1614,6 +1616,7 @@
     <ClCompile Include="..\src\frontend\shyu.c" />
     <ClCompile Include="..\src\frontend\signal_handler.c" />
     <ClCompile Include="..\src\frontend\sndprint.c" />
+    <ClCompile Include="..\src\frontend\snp2va.cpp" />
     <ClCompile Include="..\src\frontend\spec.c" />
     <ClCompile Include="..\src\frontend\spiceif.c" />
     <ClCompile Include="..\src\frontend\streams.c" />