diff --git a/src/frontend/com_hb.c b/src/frontend/com_hb.c index f69b8d00f..b87996516 100644 --- a/src/frontend/com_hb.c +++ b/src/frontend/com_hb.c @@ -17,6 +17,8 @@ makes sure the circuit is built, and runs it. #include "ngspice/fteext.h" #include "ngspice/wordlist.h" #include "ngspice/cpextern.h" +#include "ngspice/dvec.h" /* Enhancement-209: dvec_alloc + VF_/SV_ for result vectors */ +#include "ngspice/sim.h" #include "circuits.h" #include "com_hb.h" @@ -30,6 +32,63 @@ static double hbnum(const char *w) return v; } +/* Enhancement-209: publish the harmonic-balance spectrum as nutmeg vectors so the + user can plot / print / wrdata it directly. A fresh "hb" plot holds a real scale + `hbfrequency` (0, f0, 2f0, ..., K f0) and one COMPLEX vector per node, carrying + the single-sided amplitude (|.| and phase match the printed table). The plot is + left current, so `plot mag(out)` / `print out` / `wrdata sp v(out)` work at once. */ +static void +hb_publish(CKTcircuit *ckt, const struct hbspectrum *sp) +{ + int numNames = 0, error, j, k, nv = 0; + int N = sp->N, K = sp->K; + IFuid *nameList = NULL; + struct plot *pl; + struct dvec *fv; + + if (!sp->Vr || !sp->Vi || N <= 0 || K < 1) + return; + error = CKTnames(ckt, &numNames, &nameList); + if (error || numNames <= 0) { + if (nameList) + tfree(nameList); + return; + } + + pl = plot_alloc("hb"); + pl->pl_name = copy("Harmonic Balance"); + pl->pl_title = copy((ft_curckt && ft_curckt->ci_name) ? ft_curckt->ci_name : "hb"); + plot_new(pl); + plot_setcur(pl->pl_typename); + + /* the harmonic-frequency scale, created first so it becomes the plot scale */ + fv = dvec_alloc(copy("hbfrequency"), SV_FREQUENCY, + (short) (VF_REAL | VF_PERMANENT), K + 1, NULL); + for (k = 0; k <= K; k++) + fv->v_realdata[k] = k * sp->f0; + vec_new(fv); + + for (j = 0; j < numNames && j < N; j++) { + const char *nm = (const char *) nameList[j]; + int isI = (nm && strstr(nm, "#branch") != NULL); + struct dvec *v = dvec_alloc(copy(nm), isI ? SV_CURRENT : SV_VOLTAGE, + (short) (VF_COMPLEX | VF_PERMANENT), K + 1, NULL); + for (k = 0; k <= K; k++) { + double sc = (k == 0) ? 1.0 : 2.0; /* single-sided amplitude */ + size_t idx = (size_t) (k + K) * (size_t) N + (size_t) j; + v->v_compdata[k].cx_real = sc * sp->Vr[idx]; + v->v_compdata[k].cx_imag = sc * sp->Vi[idx]; + } + vec_new(v); + nv++; + } + tfree(nameList); + + fprintf(cp_out, "hb: spectrum stored in the current 'hb' plot -- 'hbfrequency' + " + "%d node vector%s (try plot mag() or wrdata out ).\n", + nv, nv == 1 ? "" : "s"); +} + void com_hb(wordlist *wl) { @@ -81,7 +140,15 @@ com_hb(wordlist *wl) ft_curckt->ci_curTask = ft_curckt->ci_defTask; ckt->CKTcurJob = ft_curckt->ci_defTask ? ft_curckt->ci_defTask->jobs : NULL; - err = HBanalyze(ckt, f0, K, P, maxiter, tol, verbose ? 1 : 0); - if (err != OK) - fprintf(cp_err, "hb: harmonic balance did not complete (error %d).\n", err); + { + struct hbspectrum sp; + memset(&sp, 0, sizeof sp); + err = HBanalyze(ckt, f0, K, P, maxiter, tol, verbose ? 1 : 0, &sp); + if (err != OK) + fprintf(cp_err, "hb: harmonic balance did not complete (error %d).\n", err); + else + hb_publish(ckt, &sp); + FREE(sp.Vr); + FREE(sp.Vi); + } } diff --git a/src/include/ngspice/cktdefs.h b/src/include/ngspice/cktdefs.h index 433a94275..d9d22327b 100644 --- a/src/include/ngspice/cktdefs.h +++ b/src/include/ngspice/cktdefs.h @@ -460,7 +460,18 @@ extern int PSSaskQuest(CKTcircuit *, JOB *, int , IFvalue *); extern int PSSsetParm(CKTcircuit *, JOB *, int , IFvalue *); extern int PSSinit(CKTcircuit *, JOB *); extern int DCpss(CKTcircuit *, int); -extern int HBanalyze(CKTcircuit *, double f0, int K, int P, int maxiter, double tol, int verbose); /* E-134 */ +/* Enhancement-209: HBanalyze hands its converged spectrum back to the frontend + (com_hb) so it can publish nutmeg vectors. When `out` is non-NULL and the run + converges, ownership of the Vr/Vi arrays passes to the caller (which frees them + after publishing); pass NULL to keep the old behaviour (printed table only). */ +struct hbspectrum { + int N; /* number of circuit unknowns (solution row stride) */ + int K; /* highest harmonic */ + double f0; /* fundamental frequency (Hz) */ + double *Vr, *Vi; /* [(2K+1)*N] two-sided Fourier coefficients, index (k+K)*N+i */ +}; +extern int HBanalyze(CKTcircuit *, double f0, int K, int P, int maxiter, double tol, int verbose, struct hbspectrum *out); /* E-134; E-209 out */ + #endif #ifdef RFSPICE diff --git a/src/spicelib/analysis/dcpss.c b/src/spicelib/analysis/dcpss.c index 73a14e5c0..61a9a7d94 100644 --- a/src/spicelib/analysis/dcpss.c +++ b/src/spicelib/analysis/dcpss.c @@ -1421,7 +1421,7 @@ hb_extract(CKTcircuit *ckt, const double *vsamp, int N, int P, int K, } int -HBanalyze(CKTcircuit *ckt, double f0, int K, int Pin, int maxiter, double tol, int verbose) +HBanalyze(CKTcircuit *ckt, double f0, int K, int Pin, int maxiter, double tol, int verbose, struct hbspectrum *out) { int N = SMPmatSize(ckt->CKTmatrix); int P = Pin > 0 ? Pin : ((8 * K < 32) ? 32 : 8 * K); @@ -1565,6 +1565,16 @@ HBanalyze(CKTcircuit *ckt, double f0, int K, int Pin, int maxiter, double tol, i } (void) w0; + /* Enhancement-209: hand the converged two-sided spectrum to the frontend so + `com_hb` can publish it as nutmeg vectors (hbfrequency + one complex vector + per node). Ownership of Vr/Vi passes to the caller; NULL them here so the + cleanup below does not free them. */ + if (out && rc == OK) { + out->N = N; out->K = K; out->f0 = f0; + out->Vr = Vr; out->Vi = Vi; + Vr = NULL; Vi = NULL; + } + FREE(Vr); FREE(Vi); FREE(IRr); FREE(IRi); FREE(Isr); FREE(Isi); FREE(Fr); FREE(Fi); FREE(Jr); FREE(Ji); FREE(Kr); FREE(Ki); FREE(vsamp); return rc;