pa-137

src/frontend/Makefile.am

@@ -55,6 +55,8 @@ libfte_la_SOURCES = \
     com_optimize.h  \
     com_qpss.c  \
     com_qpss.h  \
+    com_qpac.c  \
+    com_qpac.h  \
     com_hb.c    \
     com_hb.h    \
     com_checkpoint.c    \

src/frontend/com_commands.h

@@ -8,6 +8,7 @@ void com_altermod(wordlist *wl);
 void com_alterparam(wordlist *wl);
 void com_optimize(wordlist *wl);   /* Enhancement-130 */
 void com_qpss(wordlist *wl);       /* Enhancement-133 */
+void com_qpac(wordlist *wl);       /* Enhancement-137 */
 void com_hb(wordlist *wl);         /* Enhancement-134 */
 void com_savestate(wordlist *wl);  /* Enhancement-131 */
 void com_loadstate(wordlist *wl);  /* Enhancement-131 */

src/frontend/com_qpac.c

@@ -0,0 +1,58 @@
+/**********
+Enhancement-137: two-tone small-signal QPAC (quasi-periodic AC) -- `qpac <f_in>`.
+
+Injects a small signal at frequency f_in around the QPSS operating point retained by a
+prior `qpss <expr> <f1> <f2> hb`, and reports the response at every sideband
+f_in + k1*f1 + k2*f2 -- the two-tone analogue of PAC. The quasi-periodic operating point
+mixes the small signal to the sidebands through the same 2-D conversion matrix the QPSS
+Newton used as its Jacobian. The heavy lifting is QPACanalyze() (spicelib/analysis/
+dcpss.c); this command parses f_in and runs it.
+**********/
+
+#include "ngspice/ngspice.h"
+#include "ngspice/cpdefs.h"
+#include "ngspice/cktdefs.h"
+#include "ngspice/ftedefs.h"
+#include "ngspice/fteext.h"
+#include "ngspice/wordlist.h"
+#include "ngspice/cpextern.h"
+
+#include "com_qpac.h"
+
+static double qpacnum(const char *w)
+{
+    char *s = (char *) w;
+    double v = 0.0;
+    if (ft_numparse(&s, FALSE, &v) < 0)
+        v = atof(w);
+    return v;
+}
+
+void
+com_qpac(wordlist *wl)
+{
+    CKTcircuit *ckt;
+    double f_in;
+    int    verbose, err;
+
+    if (!ft_curckt || !ft_curckt->ci_ckt) {
+        fprintf(cp_err, "Error: qpac: there is no circuit loaded.\n");
+        return;
+    }
+    ckt = ft_curckt->ci_ckt;
+
+    if (!wl || !wl->wl_word) {
+        fprintf(cp_err, "Usage: qpac <f_in>   (run `qpss <expr> <f1> <f2> hb` first)\n");
+        return;
+    }
+    f_in = qpacnum(wl->wl_word);
+    if (f_in <= 0.0) {
+        fprintf(cp_err, "Error: qpac: need f_in > 0.\n");
+        return;
+    }
+
+    verbose = cp_getvar("qpac_verbose", CP_BOOL, NULL, 0);
+    err = QPACanalyze(ckt, f_in, verbose ? 1 : 0);
+    if (err != OK)
+        fprintf(cp_err, "qpac: quasi-periodic AC did not complete (error %d).\n", err);
+}

src/frontend/com_qpac.h

@@ -0,0 +1,7 @@
+#ifndef ngspice_COM_QPAC_H
+#define ngspice_COM_QPAC_H
+
+/* Enhancement-137: two-tone small-signal QPAC (quasi-periodic AC). */
+void com_qpac(wordlist *wl);
+
+#endif

src/frontend/commands.c

@@ -436,6 +436,10 @@ struct comm spcp_coms[] = {
       { 040, 040, 040, 040 }, E_DEFHMASK, 3, LOTS,
       NULL,
       "expr f1 f2 [periods] [maxorder] | expr f1 f2 hb [K1] [K2] : two-tone quasi-periodic steady state (transient / frequency-domain HB, incommensurate-capable)." },
+    { "qpac", com_qpac, TRUE, FALSE,             /* Enhancement-137 */
+      { 040, 040, 040, 040 }, E_DEFHMASK, 1, LOTS,
+      NULL,
+      "f_in : two-tone small-signal QPAC -- response at sidebands f_in+k1f1+k2f2 around the `qpss ... hb` operating point." },
     { "hb", com_hb, TRUE, FALSE,                 /* Enhancement-134 */
       { 040, 040, 040, 040 }, E_DEFHMASK, 2, LOTS,
       NULL,

src/include/ngspice/cktdefs.h

@@ -477,6 +477,7 @@ struct hbspectrum {
 };
 extern int HBanalyze(CKTcircuit *, double f0, int K, int P, int maxiter, double tol, int verbose, struct hbspectrum *out); /* E-134; E-209 out */
 extern int QPSShb(CKTcircuit *, double f1, double f2, int K1, int K2, int P1, int P2, int maxiter, double tol, int verbose); /* E-136 */
+extern int QPACanalyze(CKTcircuit *, double f_in, int verbose); /* E-137 */
 
 #endif
 

src/spicelib/analysis/dcpss.c

@@ -1667,24 +1667,26 @@ struct qp_harm {
     double *Gmr, *Gmi, *Cmr, *Cmi;   /* [nnz*Dsz] 2-D difference spectra */
     double f1, f2;
     double *Vr, *Vi;           /* [Ntot] retained operating point (for qpac); else NULL */
+    double *B0r, *B0i;         /* [N] AC-source stimulus for qpac (E-137); else NULL */
+    int     has_src;           /* 1 if a netlist AC source stamped B0 */
 };
 
-struct qp_harm *qpss_hb_saved = NULL;   /* retained QPSS op-point for qpac (E-137) */
+static struct qp_harm *qpss_hb_saved = NULL;   /* retained QPSS op-point for qpac (E-137) */
 
 static int qp_didx(const struct qp_harm *h, int d1, int d2)
 { return (d1 + 2*h->K1) * h->D2c + (d2 + 2*h->K2); }
 
 /* free the arrays a qp_harm owns (not the struct itself) */
-void qp_free(struct qp_harm *hd)
+static void qp_free(struct qp_harm *hd)
 {
     if (!hd) return;
     FREE(hd->rr); FREE(hd->cc); FREE(hd->h1); FREE(hd->h2);
     FREE(hd->Gmr); FREE(hd->Gmi); FREE(hd->Cmr); FREE(hd->Cmi);
-    FREE(hd->Vr); FREE(hd->Vi);
+    FREE(hd->Vr); FREE(hd->Vi); FREE(hd->B0r); FREE(hd->B0i);
 }
 
 /* assemble the dense Ntot x Ntot 2-D conversion matrix at input freq f_in */
-void qp_build_matrix(struct qp_harm *hd, double f_in, double *Ar, double *Ai)
+static void qp_build_matrix(struct qp_harm *hd, double f_in, double *Ar, double *Ai)
 {
     int Nh = hd->Nh, N = hd->N, Ntot = hd->Ntot, nnz = hd->nnz, ni, mi, e;
     memset(Ar, 0, (size_t)Ntot * (size_t)Ntot * sizeof(double));
@@ -2106,6 +2108,20 @@ QPSShb(CKTcircuit *ckt, double f1, double f2, int K1, int K2, int P1, int P2,
                 sv->Vr = TMALLOC(double, Ntot); sv->Vi = TMALLOC(double, Ntot);
                 memcpy(sv->Vr, Vr, (size_t)Ntot*sizeof(double));
                 memcpy(sv->Vi, Vi, (size_t)Ntot*sizeof(double));
+                /* capture the AC-source stimulus B0 (a netlist `AC`-flagged source's
+                 * RHS -- bias-independent) at the op-point, for qpac (E-137). */
+                sv->B0r = TMALLOC(double, N); sv->B0i = TMALLOC(double, N); sv->has_src = 0;
+                for (i = 0; i <= N; i++) { ckt->CKTrhsOld[i] = 0.0; ckt->CKTrhs[i] = 0.0; }
+                ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG;
+                CKTload(ckt);
+                ckt->CKTomega = 1.0;
+                ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEAC;
+                CKTacLoad(ckt);
+                for (i = 0; i < N; i++) {
+                    sv->B0r[i] = ckt->CKTrhs[i+1];
+                    sv->B0i[i] = ckt->CKTirhs[i+1];
+                    if (sv->B0r[i] != 0.0 || sv->B0i[i] != 0.0) sv->has_src = 1;
+                }
                 if (qpss_hb_saved) { qp_free(qpss_hb_saved); FREE(qpss_hb_saved); }
                 qpss_hb_saved = sv;
             } else {
@@ -2120,6 +2136,80 @@ QPSShb(CKTcircuit *ckt, double f1, double f2, int K1, int K2, int P1, int P2,
 }
 
 
+/* ======================================================================
+ * Enhancement-137: two-tone small-signal QPAC -- the quasi-periodic
+ * analogue of PAC. Around the QPSS operating point retained by
+ * `qpss ... hb` (qpss_hb_saved), inject a small signal at f_in; the
+ * quasi-periodic operating point converts it to the sidebands
+ * f_in + k1*f1 + k2*f2 through the SAME 2-D conversion matrix
+ * (qp_build_matrix at f_in), solved by pss_csolve. Same construction
+ * as pac_solve_at (E-121/122), on the two-tone harmonic set.
+ * ====================================================================== */
+int
+QPACanalyze(CKTcircuit *ckt, double f_in, int verbose)
+{
+    struct qp_harm *hd = qpss_hb_saved;
+    int    N, Nh, Ntot, i, hi, i00, numNames, error, k1, k2, ord;
+    IFuid *nameList = NULL;
+    double *Ar, *Ai, *Xr, *Xi;
+
+    if (!hd) {
+        fprintf(stderr, "qpac: no QPSS operating point -- run `qpss <expr> <f1> <f2> hb` first.\n");
+        return E_NOTFOUND;
+    }
+    N = hd->N; Nh = hd->Nh; Ntot = hd->Ntot;
+    Ar = TMALLOC(double, (size_t)Ntot*(size_t)Ntot);
+    Ai = TMALLOC(double, (size_t)Ntot*(size_t)Ntot);
+    Xr = TMALLOC(double, Ntot); Xi = TMALLOC(double, Ntot);
+    (void) verbose;
+
+    qp_build_matrix(hd, f_in, Ar, Ai);
+
+    /* stimulus in the (0,0) sideband: netlist AC source RHS, or unit current at node 1 */
+    memset(Xr, 0, (size_t)Ntot*sizeof(double));
+    memset(Xi, 0, (size_t)Ntot*sizeof(double));
+    i00 = hd->K1 * (2*hd->K2 + 1) + hd->K2;          /* harmonic index of (0,0) */
+    if (hd->has_src) {
+        for (i = 0; i < N; i++) {
+            Xr[(size_t)i00*(size_t)N + (size_t)i] = hd->B0r[i];
+            Xi[(size_t)i00*(size_t)N + (size_t)i] = hd->B0i[i];
+        }
+    } else {
+        Xr[(size_t)i00*(size_t)N + 0] = 1.0;
+    }
+    if (pss_csolve(Ntot, Ar, Ai, Xr, Xi)) {
+        fprintf(stderr, "qpac: singular conversion matrix.\n");
+        FREE(Ar); FREE(Ai); FREE(Xr); FREE(Xi);
+        return E_SINGULAR;
+    }
+    FREE(Ar); FREE(Ai);
+
+    error = CKTnames(ckt, &numNames, &nameList);
+    fprintf(stdout,
+            "\nQPAC: two-tone small-signal response (f_in = %g Hz, f1 = %g Hz, f2 = %g Hz)\n"
+            "  node      (k1,k2)   sideband f_in+k1f1+k2f2 [Hz]   |response|      phase [deg]\n",
+            f_in, hd->f1, hd->f2);
+    for (i = 0; i < N; i++) {
+        const char *nm = (!error && i < numNames) ? (const char *) nameList[i] : "?";
+        for (ord = 0; ord <= hd->K1 + hd->K2; ord++)
+            for (k1 = -hd->K1; k1 <= hd->K1; k1++)
+                for (k2 = -hd->K2; k2 <= hd->K2; k2++) {
+                    double fsb, xr, xi;
+                    if (abs(k1) + abs(k2) != ord) continue;
+                    hi = (k1 + hd->K1) * (2*hd->K2 + 1) + (k2 + hd->K2);
+                    fsb = f_in + k1*hd->f1 + k2*hd->f2;
+                    xr = Xr[(size_t)hi*(size_t)N + (size_t)i];
+                    xi = Xi[(size_t)hi*(size_t)N + (size_t)i];
+                    fprintf(stdout, "  %-8s  (%2d,%2d)   %18.6e   %14.6e   %10.3f\n",
+                            nm, k1, k2, fsb, hypot(xr, xi), atan2(xi, xr) * 180.0/M_PI);
+                }
+    }
+    if (nameList) tfree(nameList);
+    FREE(Xr); FREE(Xi);
+    return OK;
+}
+
+
 int
 DCpss(CKTcircuit *ckt,
        int restart)   /* forced restart flag */

visualc/vngspice.vcxproj

@@ -895,6 +895,7 @@
     <ClInclude Include="..\src\frontend\com_option.h" />
     <ClInclude Include="..\src\frontend\com_plot.h" />
     <ClInclude Include="..\src\frontend\com_pyplot.h" />
+    <ClInclude Include="..\src\frontend\com_qpac.h" />
     <ClInclude Include="..\src\frontend\com_qpss.h" />
     <ClInclude Include="..\src\frontend\com_rehash.h" />
     <ClInclude Include="..\src\frontend\com_set.h" />
@@ -1513,6 +1514,7 @@
     <ClCompile Include="..\src\frontend\com_option.c" />
     <ClCompile Include="..\src\frontend\com_plot.c" />
     <ClCompile Include="..\src\frontend\com_pyplot.c" />
+    <ClCompile Include="..\src\frontend\com_qpac.c" />
     <ClCompile Include="..\src\frontend\com_qpss.c" />
     <ClCompile Include="..\src\frontend\com_rehash.c" />
     <ClCompile Include="..\src\frontend\com_set.c" />