diff --git a/src/frontend/Makefile.am b/src/frontend/Makefile.am index 8dd94a014..3aac34495 100644 --- a/src/frontend/Makefile.am +++ b/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 \ diff --git a/src/frontend/com_commands.h b/src/frontend/com_commands.h index 0fcbfb23b..aa29a374a 100644 --- a/src/frontend/com_commands.h +++ b/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 */ diff --git a/src/frontend/com_qpac.c b/src/frontend/com_qpac.c new file mode 100644 index 000000000..e6d811c25 --- /dev/null +++ b/src/frontend/com_qpac.c @@ -0,0 +1,58 @@ +/********** +Enhancement-137: two-tone small-signal QPAC (quasi-periodic AC) -- `qpac `. + +Injects a small signal at frequency f_in around the QPSS operating point retained by a +prior `qpss 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 (run `qpss 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); +} diff --git a/src/frontend/com_qpac.h b/src/frontend/com_qpac.h new file mode 100644 index 000000000..bc91afdfa --- /dev/null +++ b/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 diff --git a/src/frontend/commands.c b/src/frontend/commands.c index 35a7f16fc..5c61229b7 100644 --- a/src/frontend/commands.c +++ b/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, diff --git a/src/include/ngspice/cktdefs.h b/src/include/ngspice/cktdefs.h index c967f0296..f6c695ab8 100644 --- a/src/include/ngspice/cktdefs.h +++ b/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 diff --git a/src/spicelib/analysis/dcpss.c b/src/spicelib/analysis/dcpss.c index 63f6ff4e2..eedfff318 100644 --- a/src/spicelib/analysis/dcpss.c +++ b/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 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 */ diff --git a/visualc/vngspice.vcxproj b/visualc/vngspice.vcxproj index b45c11426..4a49b3b71 100644 --- a/visualc/vngspice.vcxproj +++ b/visualc/vngspice.vcxproj @@ -895,6 +895,7 @@ + @@ -1513,6 +1514,7 @@ +