This commit is contained in:
Meisam Bahadori 2026-07-26 16:59:56 +02:00 committed by Holger Vogt
parent 3ab718aabd
commit 648dd4bb8f
3 changed files with 93 additions and 5 deletions

View File

@ -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(<node>) or wrdata out <node>).\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);
}
}

View File

@ -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

View File

@ -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;