luke
/
ngspice
mirror of https://git.code.sf.net/p/ngspice/ngspice
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Set a framework to integrate Harmonic Balance

This commit is contained in:
Holger Vogt 2025-12-05 15:55:22 +01:00
parent 83e2b284f3
commit c5891f137c
8 changed files with 1249 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/include/ngspice/cktdefs.h
View File

@ -473,6 +473,12 @@ extern int CKTspCalcPowerWave(CKTcircuit* ckt);
extern int CKTspCalcSMatrix(CKTcircuit* ckt);
#endif
#ifdef WITH_HB
extern int HBan(CKTcircuit*, int);
extern int HBaskQuest(CKTcircuit*, JOB*, int, IFvalue*);
extern int HBsetParm(CKTcircuit*, JOB*, int, IFvalue*);
#endif
#ifdef __cplusplus
extern "C"
{

51
src/include/ngspice/hbardefs.h Normal file
View File

@ -0,0 +1,51 @@
/**********
Copyright 1990 Regents of the University of California. All rights reserved.
Author: 1985 Thomas L. Quarles
**********/
#ifndef ngspice_SPDEFS_H
#define ngspice_SPDEFS_H
#include "ngspice/jobdefs.h"
#ifdef WITH_HB
/* structure used to describe an AC analysis to be performed */
typedef struct {
int JOBtype;
JOB *JOBnextJob; /* pointer to next thing to do */
char *JOBname; /* name of this job */
double SPstartFreq;
double SPstopFreq;
double SPfreqDelta; /* multiplier for decade/octave stepping, */
/* step for linear steps. */
double SPsaveFreq; /* frequency at which we left off last time*/
int SPstepType; /* values described below */
int SPnumberSteps;
unsigned SPdoNoise : 1; /* Flag to indicate if SP noise must be calculated*/
int SPnoiseInput;
int SPnoiseOutput;
} HBAN;
/* available step types: XXX should be somewhere else */
#ifndef ngspice_ACDEFS_H
enum {
DECADE = 1,
OCTAVE,
LINEAR,
};
#endif
enum {
SP_DEC = 1,
SP_OCT,
SP_LIN,
SP_START,
SP_STOP,
SP_STEPS,
SP_DONOISE,
};
#endif
#endif

989
src/spicelib/analysis/hban.c Normal file
View File

@ -0,0 +1,989 @@
/*
****
* Holger Vogt 2025
* derived from span.c, Alessio Cacciatori 2021
****
*/
#include "ngspice/ngspice.h"
#include "ngspice/cktdefs.h"
#include "ngspice/hbardefs.h"
#include "ngspice/devdefs.h"
#include "ngspice/sperror.h"
#ifdef XSPICE
#include "ngspice/evt.h"
#include "ngspice/enh.h"
/* gtri - add - wbk - 12/19/90 - Add headers */
#include "ngspice/mif.h"
#include "ngspice/evtproto.h"
#include "ngspice/ipctiein.h"
/* gtri - end - wbk */
#endif
#define SQR(x) ((x) * (x))
#ifdef WITH_HB
#include "vsrc/vsrcdefs.h"
#include "../maths/dense/dense.h"
#include "../maths/dense/denseinlines.h"
#if (0)
int CKTspnoise(CKTcircuit* ckt, int mode, int operation, Ndata* data, NOISEAN* noisean);
int NInspIter(CKTcircuit* ckt, VSRCinstance* port);
NOISEAN* SPcreateNoiseAnalysys(CKTcircuit* ckt);
#endif
int initHBmatrix(CKTcircuit* ckt, int doNoise);
void deleteHBmatrix(CKTcircuit* ckt);
#define INIT_STATS() \
do { \
startTime = SPfrontEnd->IFseconds(); \
startdTime = ckt->CKTstat->STATdecompTime; \
startsTime = ckt->CKTstat->STATsolveTime; \
startlTime = ckt->CKTstat->STATloadTime; \
startkTime = ckt->CKTstat->STATsyncTime; \
} while(0)
#define UPDATE_STATS(analysis) \
do { \
ckt->CKTcurrentAnalysis = analysis; \
ckt->CKTstat->STATacTime += SPfrontEnd->IFseconds() - startTime; \
ckt->CKTstat->STATacDecompTime += ckt->CKTstat->STATdecompTime - startdTime; \
ckt->CKTstat->STATacSolveTime += ckt->CKTstat->STATsolveTime - startsTime; \
ckt->CKTstat->STATacLoadTime += ckt->CKTstat->STATloadTime - startlTime; \
ckt->CKTstat->STATacSyncTime += ckt->CKTstat->STATsyncTime - startkTime; \
} while(0)
#if(0)
/*----------------------------------
* Auxiliary data for S-Y-Z matrix
* conversion
*-----------------------------------
*/
CMat* eyem = NULL;
CMat* zref = NULL;
CMat* gn = NULL;
CMat* gninv = NULL;
CMat* vNoise = NULL;
CMat* iNoise = NULL;
// Aux data for Noise Calculation
double NF = 0;
double Rn = 0;
cplx Sopt;
double Fmin = 0;
double refPortY0;
int
CKTspnoise(CKTcircuit* ckt, int mode, int operation, Ndata* data, NOISEAN* noisean)
{
// Temporarily assign current job as a (dummy) NOISEAN analysis
// This is needed to avoid
HBAN* oldJob = (HBAN*)ckt->CKTcurJob;
ckt->CKTcurJob = (JOB*)noisean;
double outNdens;
int i;
int error;
outNdens = 0.0;
/* let each device decide how many and what type of noise sources it has */
for (i = 0; i < DEVmaxnum; i++) {
if (DEVices[i] && DEVices[i]->DEVnoise && ckt->CKThead[i]) {
int a = 0;
a++;
if (a == 0) a = 2;
error = DEVices[i]->DEVnoise(mode, operation, ckt->CKThead[i],
ckt, data, &outNdens);
if (error)
{
ckt->CKTcurJob = (JOB*)oldJob;
return (error);
}
}
}
switch (operation) {
case N_OPEN:
// Init all matrices
cinit(ckt->CKTNoiseCYmat, 0.0, 0.0);
cinit(ckt->CKTadjointRHS, 0.0, 0.0);
break;
case N_CALC:
{
// We have the Cy noise matrix,
// Equations from Stephen Maas 'Noise'
double knorm = 4.0 * CONSTboltz * (ckt->CKTtemp);
CMat* tempCy = cscalarmultiply(ckt->CKTNoiseCYmat, 1.0 / knorm); // cmultiply(, YConj);
#ifdef TRACE
printf("spnoise: CKTNoiseCYmat / (4*k*T)\n");
showcmat(tempCy);
#endif
if (ckt->CKTportCount == 2)
{
double Y21mod = cmodsqr(ckt->CKTYmat->d[1][0]);
Rn = (tempCy->d[1][1].re / Y21mod);
if (fabs(Rn) < 1e-30)
Rn = 1e-30;
if ((fabs(tempCy->d[1][1].re) < 1e-30) && (fabs(tempCy->d[1][1].im) < 1e-30))
{
tempCy->d[1][1].re = 1e-30;
}
cplx Ycor = csubco(ckt->CKTYmat->d[0][0],
cmultco(
cdivco(tempCy->d[0][1], tempCy->d[1][1]),
ckt->CKTYmat->d[1][0]
));
double Y11_Ycor = cmodsqr(csubco(ckt->CKTYmat->d[0][0], Ycor));
double Gu = tempCy->d[0][0].re - Rn * Y11_Ycor;
cplx Ysopt; Ysopt.re = sqrt(SQR(Ycor.re) + Gu / Rn); Ysopt.im = -Ycor.im;
cplx Y0; Y0.re = refPortY0; Y0.im = 0.0;
Sopt = cdivco(csubco(Y0, Ysopt),
caddco(Y0, Ysopt));
Fmin = 1.0 + 2.0 * Rn * (Ycor.re + Ysopt.re);
double Ysoptmod = cmodu(csubco(Y0, Ysopt));
NF = Fmin + (Rn / Y0.re) * SQR(Ysoptmod);
Fmin = 10.0 * log10(Fmin);
NF = 10.0 * log10(NF);
}
freecmat(tempCy);
}
break;
case N_CLOSE:
SPfrontEnd->OUTendPlot(data->NplotPtr);
FREE(data->namelist);
FREE(data->outpVector);
FREE(data->squared_value);
freecmat(ckt->CKTNoiseCYmat);
freecmat(ckt->CKTadjointRHS);
ckt->CKTNoiseCYmat = NULL;
ckt->CKTadjointRHS = NULL;
break;
default:
ckt->CKTcurJob = (JOB*)oldJob;
return (E_INTERN);
}
ckt->CKTcurJob = (JOB*)oldJob;
return (OK);
}
#endif
int
NInhbIter(CKTcircuit* ckt, VSRCinstance* port)
{
int i;
/* clear out the right hand side vector */
for (i = 0; i <= SMPmatSize(ckt->CKTmatrix); i++) {
ckt->CKTrhs[i] = 0.0;
ckt->CKTirhs[i] = 0.0;
}
ckt->CKTrhs[port->VSRCposNode] = 1.0; /* apply unit current excitation */
ckt->CKTrhs[port->VSRCnegNode] = -1.0;
SMPcaSolve(ckt->CKTmatrix, ckt->CKTrhs, ckt->CKTirhs, ckt->CKTrhsSpare,
ckt->CKTirhsSpare);
ckt->CKTrhs[0] = 0.0;
ckt->CKTirhs[0] = 0.0;
return (OK);
}
int initHBmatrix(CKTcircuit* ckt, int doNoise)
{
if (ckt->CKTAmat != NULL) freecmat(ckt->CKTAmat);
if (ckt->CKTBmat != NULL) freecmat(ckt->CKTBmat);
if (ckt->CKTSmat != NULL) freecmat(ckt->CKTSmat);
if (ckt->CKTYmat != NULL) freecmat(ckt->CKTYmat);
if (ckt->CKTZmat != NULL) freecmat(ckt->CKTZmat);
#if (0)
if (eyem != NULL) freecmat(eyem);
if (zref != NULL) freecmat(zref);
if (gn != NULL) freecmat(gn);
if (gninv != NULL) freecmat(gninv);
#endif
ckt->CKTAmat = newcmat(ckt->CKTportCount, ckt->CKTportCount, 0.0, 0.0);
if (ckt->CKTAmat == NULL)
return (E_NOMEM);
ckt->CKTBmat = newcmat(ckt->CKTportCount, ckt->CKTportCount, 0.0, 0.0);
if (ckt->CKTBmat == NULL)
return (3);
ckt->CKTSmat = newcmat(ckt->CKTportCount, ckt->CKTportCount, 0.0, 0.0);
if (ckt->CKTSmat == NULL)
return (E_NOMEM);
ckt->CKTYmat = newcmat(ckt->CKTportCount, ckt->CKTportCount, 0.0, 0.0);
if (ckt->CKTYmat == NULL)
return (E_NOMEM);
ckt->CKTZmat = newcmat(ckt->CKTportCount, ckt->CKTportCount, 0.0, 0.0);
if (ckt->CKTZmat == NULL)
return (E_NOMEM);
#if (0)
eyem = ceye(ckt->CKTportCount);
if (eyem == NULL)
return (E_NOMEM);
zref = newcmat(ckt->CKTportCount, ckt->CKTportCount, 0.0, 0.0);
if (zref == NULL)
return (E_NOMEM);
gn = newcmat(ckt->CKTportCount, ckt->CKTportCount, 0.0, 0.0);
if (gn == NULL)
return (E_NOMEM);
gninv = newcmat(ckt->CKTportCount, ckt->CKTportCount, 0.0, 0.0);
if (gninv == NULL)
return (E_NOMEM);
// Now that we have found the model, we may init the Zref and Gn ports
if (ckt->CKTVSRCid >= 0)
VSRCspinit(ckt->CKThead[ckt->CKTVSRCid], ckt, zref, gn, gninv);
if (doNoise)
{
// Allocate matrices and vector
if (ckt->CKTNoiseCYmat != NULL) freecmat(ckt->CKTNoiseCYmat);
ckt->CKTNoiseCYmat = newcmatnoinit(ckt->CKTportCount, ckt->CKTportCount);
if (ckt->CKTNoiseCYmat == NULL) return (E_NOMEM);
// Use CKTadjointRHS as a convenience storage for all solutions (each solution per each
// port excitation)
if (ckt->CKTadjointRHS != NULL) freecmat(ckt->CKTadjointRHS);
ckt->CKTadjointRHS = newcmatnoinit(ckt->CKTportCount, ckt->CKTmaxEqNum);
if (ckt->CKTadjointRHS == NULL) return (E_NOMEM);
if (vNoise != NULL) freecmat(vNoise);
if (iNoise != NULL) freecmat(iNoise);
vNoise = newcmatnoinit(1, ckt->CKTportCount);
iNoise = newcmatnoinit(1, ckt->CKTportCount);
VSRCinstance* refPort = (VSRCinstance*)(ckt->CKTrfPorts[0]);
refPortY0 = refPort->VSRCportY0;
}
#endif
return (OK);
}
void deleteHBmatrix(CKTcircuit* ckt)
{
if (ckt->CKTAmat != NULL) freecmat(ckt->CKTAmat);
if (ckt->CKTBmat != NULL) freecmat(ckt->CKTBmat);
if (ckt->CKTSmat != NULL) freecmat(ckt->CKTSmat);
if (ckt->CKTYmat != NULL) freecmat(ckt->CKTYmat);
if (ckt->CKTZmat != NULL) freecmat(ckt->CKTZmat);
#if (0)
if (eyem != NULL) freecmat(eyem);
if (zref != NULL) freecmat(zref);
if (gn != NULL) freecmat(gn);
if (gninv != NULL) freecmat(gninv);
eyem = NULL;
zref = NULL;
gn = NULL;
gninv = NULL;
ckt->CKTAmat = NULL;
ckt->CKTBmat = NULL;
ckt->CKTSmat = NULL;
ckt->CKTZmat = NULL;
ckt->CKTYmat = NULL;
if (ckt->CKTNoiseCYmat != NULL) freecmat(ckt->CKTNoiseCYmat);
if (ckt->CKTadjointRHS != NULL) freecmat(ckt->CKTadjointRHS);
if (vNoise != NULL) freecmat(vNoise);
if (iNoise != NULL) freecmat(iNoise);
vNoise = NULL;
iNoise = NULL;
ckt->CKTNoiseCYmat = NULL;
ckt->CKTadjointRHS = NULL;
#endif
}
#if (0)
NOISEAN* SPcreateNoiseAnalysys(CKTcircuit* ckt)
{
NOISEAN* internalNoiseAN = TMALLOC(NOISEAN, 1);
if (internalNoiseAN==NULL) return NULL;
SPAN* span = (SPAN*)ckt->CKTcurJob;
internalNoiseAN->NstartFreq = span->SPstartFreq;
internalNoiseAN->NstopFreq = span->SPstopFreq;
internalNoiseAN->NStpsSm = 1; // Force to output noise at every freq step
internalNoiseAN->JOBnextJob = NULL;
internalNoiseAN->JOBtype = span->JOBtype;
internalNoiseAN->JOBname = NULL;
internalNoiseAN->NfreqDelta = span->SPfreqDelta;
internalNoiseAN->NstpType = span->SPstepType;
internalNoiseAN->NnumSteps = span->SPnumberSteps;
return internalNoiseAN;
}
#endif
int
HBan(CKTcircuit* ckt, int restart)
{
HBAN* job = (HBAN*)ckt->CKTcurJob;
double freq;
double freqTol; /* tolerence parameter for finding final frequency */
double startdTime;
double startsTime;
double startlTime;
double startkTime;
double startTime;
int error;
int numNames;
int i;
IFuid* nameList; /* va: tmalloc'ed list of names */
IFuid freqUid;
static runDesc* spPlot = NULL;
runDesc* plot = NULL;
double* rhswoPorts = NULL;
double* irhswoPorts = NULL;
NOISEAN* internalNoiseAN = NULL;
// Noise analysis is performed at each freq of the SP Analysis
// A temporary dummy job is therefore created
/* variable must be static, for continuation of interrupted (Ctrl-C),
longer lasting noise anlysis */
static Ndata* data = NULL;
if (job->SPdoNoise)
{
data = TMALLOC(Ndata, 1);
}
if (ckt->CKTportCount == 0)
{
fprintf(stderr, "\nError: No RF Port is present, cannot run sp analysis\n");
controlled_exit(EXIT_BAD);
}
if (ckt->CKTportCount == 1)
{
fprintf(stderr, "\nError: Only one RF Port is found, we need at least two!\n");
controlled_exit(EXIT_BAD);
}
#ifdef XSPICE
/* Tell the code models what mode we're in */
g_mif_info.circuit.anal_type = MIF_DC;
g_mif_info.circuit.anal_init = MIF_TRUE;
#endif
/* start at beginning */
if (job->SPsaveFreq == 0 || restart) {
if (job->SPnumberSteps < 1)
job->SPnumberSteps = 1;
switch (job->SPstepType) {
case DECADE:
if (job->SPstartFreq <= 0) {
fprintf(stderr, "ERROR: AC startfreq <= 0\n");
return E_PARMVAL;
}
job->SPfreqDelta =
exp(log(10.0) / job->SPnumberSteps);
break;
case OCTAVE:
if (job->SPstartFreq <= 0) {
fprintf(stderr, "ERROR: AC startfreq <= 0\n");
return E_PARMVAL;
}
job->SPfreqDelta =
exp(log(2.0) / job->SPnumberSteps);
break;
case LINEAR:
if (job->SPnumberSteps - 1 > 1)
job->SPfreqDelta =
(job->SPstopFreq -
job->SPstartFreq) /
(job->SPnumberSteps - 1);
else
/* Patch from: Richard McRoberts
* This patch is for a rather pathological case:
* a linear step with only one point */
job->SPfreqDelta = 0;
break;
default:
return(E_BADPARM);
}
if (job->SPdoNoise)
{
data->lstFreq = job->SPstartFreq - 1;
data->delFreq = 0.0;
}
#ifdef XSPICE
/* gtri - begin - wbk - Call EVTop if event-driven instances exist */
if (ckt->evt->counts.num_insts != 0) {
error = EVTop(ckt,
(ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITJCT,
(ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITFLOAT,
ckt->CKTdcMaxIter,
MIF_TRUE);
EVTdump(ckt, IPC_ANAL_DCOP, 0.0);
EVTop_save(ckt, MIF_TRUE, 0.0);
}
else
#endif
/* If no event-driven instances, do what SPICE normally does */
if (!ckt->CKTnoopac) { /* skip OP if option NOOPAC is set and circuit is linear */
error = CKTop(ckt,
(ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITJCT,
(ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITFLOAT,
ckt->CKTdcMaxIter);
if (error) {
fprintf(stdout, "\nAC operating point failed -\n");
CKTncDump(ckt);
return(error);
}
}
else
fprintf(stdout, "\n Linear circuit, option noopac given: no OP analysis\n");
ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG;
error = CKTload(ckt);
if (error) return(error);
error = CKTnames(ckt, &numNames, &nameList);
if (error) return(error);
if (ckt->CKTkeepOpInfo) {
/* Dump operating point. */
error = SPfrontEnd->OUTpBeginPlot(ckt, ckt->CKTcurJob,
"AC Operating Point",
NULL, IF_REAL,
numNames, nameList, IF_REAL,
&plot);
if (error) return(error);
CKTdump(ckt, 0.0, plot);
SPfrontEnd->OUTendPlot(plot);
plot = NULL;
}
int extraSPdataLength = 3 * ckt->CKTportCount * ckt->CKTportCount;
if (job->SPdoNoise)
{
extraSPdataLength += ckt->CKTportCount * ckt->CKTportCount; // Add Cy
if (ckt->CKTportCount == 2)
extraSPdataLength += 4;
}
nameList = (IFuid*)TREALLOC(IFuid, nameList, numNames + extraSPdataLength);
// Create UIDs
for (int dest = 1; dest <= ckt->CKTportCount; dest++)
for (int j = 1; j <= ckt->CKTportCount; j++)
{
char tmpBuf[32];
sprintf(tmpBuf, "S_%d_%d", dest, j);
SPfrontEnd->IFnewUid(ckt, &(nameList[numNames++]), NULL, tmpBuf, UID_OTHER, NULL);
}
// Create UIDs
for (int dest = 1; dest <= ckt->CKTportCount; dest++)
for (int j = 1; j <= ckt->CKTportCount; j++)
{
char tmpBuf[32];
sprintf(tmpBuf, "Y_%d_%d", dest, j);
SPfrontEnd->IFnewUid(ckt, &(nameList[numNames++]), NULL, tmpBuf, UID_OTHER, NULL);
}
// Create UIDs
for (int dest = 1; dest <= ckt->CKTportCount; dest++)
for (int j = 1; j <= ckt->CKTportCount; j++)
{
char tmpBuf[32];
sprintf(tmpBuf, "Z_%d_%d", dest, j);
SPfrontEnd->IFnewUid(ckt, &(nameList[numNames++]), NULL, tmpBuf, UID_OTHER, NULL);
}
// Add noise related output, if needed
if (job->SPdoNoise)
{
// Create UIDs
for (int dest = 1; dest <= ckt->CKTportCount; dest++)
for (int j = 1; j <= ckt->CKTportCount; j++)
{
char tmpBuf[32];
sprintf(tmpBuf, "Cy_%d_%d", dest, j);
SPfrontEnd->IFnewUid(ckt, &(nameList[numNames++]), NULL, tmpBuf, UID_OTHER, NULL);
}
// Add NFMin, SOpt, Rn (related to port 1)
if (ckt->CKTportCount == 2)
{
SPfrontEnd->IFnewUid(ckt, &(nameList[numNames++]), NULL, "NF", UID_OTHER, NULL);
SPfrontEnd->IFnewUid(ckt, &(nameList[numNames++]), NULL, "SOpt", UID_OTHER, NULL);
SPfrontEnd->IFnewUid(ckt, &(nameList[numNames++]), NULL, "NFmin", UID_OTHER, NULL);
SPfrontEnd->IFnewUid(ckt, &(nameList[numNames++]), NULL, "Rn", UID_OTHER, NULL);
}
}
SPfrontEnd->IFnewUid(ckt, &freqUid, NULL, "frequency", UID_OTHER, NULL);
error = SPfrontEnd->OUTpBeginPlot(ckt, ckt->CKTcurJob,
ckt->CKTcurJob->JOBname,
freqUid, IF_REAL,
numNames, nameList, IF_COMPLEX,
&spPlot);
tfree(nameList);
if (error) return(error);
if (job->SPstepType != LINEAR) {
SPfrontEnd->OUTattributes(spPlot, NULL, OUT_SCALE_LOG, NULL);
}
freq = job->SPstartFreq;
}
else { /* continue previous analysis */
freq = job->SPsaveFreq;
job->SPsaveFreq = 0; /* clear the 'old' frequency */
/* fix resume? saj, indeed !*/
error = SPfrontEnd->OUTpBeginPlot(NULL, NULL,
NULL,
NULL, 0,
666, NULL, 666,
&spPlot);
/* saj*/
}
switch (job->SPstepType) {
case DECADE:
case OCTAVE:
freqTol = job->SPfreqDelta *
job->SPstopFreq * ckt->CKTreltol;
break;
case LINEAR:
freqTol = job->SPfreqDelta * ckt->CKTreltol;
break;
default:
return(E_BADPARM);
}
#ifdef XSPICE
/* gtri - add - wbk - 12/19/90 - Set anal_init and anal_type */
g_mif_info.circuit.anal_init = MIF_TRUE;
/* Tell the code models what mode we're in */
g_mif_info.circuit.anal_type = MIF_AC;
/* gtri - end - wbk */
#endif
INIT_STATS();
ckt->CKTcurrentAnalysis = DOING_AC | DOING_SP;
if (initSPmatrix(ckt, job->SPdoNoise))
return (E_NOMEM);
// Create Noise UID, if needed
if (job->SPdoNoise)
{
internalNoiseAN = SPcreateNoiseAnalysys(ckt);
if (internalNoiseAN == NULL)
return (E_NOMEM);
data->numPlots = 0; /* we don't have any plots yet */
data->freq = freq;
error = CKTspnoise(ckt, N_DENS, N_OPEN, data, internalNoiseAN);
if (error) {
tfree(internalNoiseAN);
return(error);
}
}
#ifdef KLU
if (ckt->CKTmatrix->CKTkluMODE)
{
/* Conversion from Real Matrix to Complex Matrix */
if (!ckt->CKTmatrix->SMPkluMatrix->KLUmatrixIsComplex)
{
for (i = 0; i < DEVmaxnum; i++)
if (DEVices[i] && DEVices[i]->DEVbindCSCComplex && ckt->CKThead[i])
DEVices[i]->DEVbindCSCComplex(ckt->CKThead[i], ckt);
ckt->CKTmatrix->SMPkluMatrix->KLUmatrixIsComplex = KLUMatrixComplex;
}
}
#endif
ckt->CKTactivePort = 0;
/* main loop through all scheduled frequencies */
while (freq <= job->SPstopFreq + freqTol) {
int activePort = 0;
//
if (SPfrontEnd->IFpauseTest()) {
/* user asked us to pause via an interrupt */
job->SPsaveFreq = freq;
return(E_PAUSE);
}
ckt->CKTomega = 2.0 * M_PI * freq;
/* Update opertating point, if variable 'hertz' is given */
if (ckt->CKTvarHertz) {
#ifdef KLU
if (ckt->CKTmatrix->CKTkluMODE)
{
/* Conversion from Complex Matrix to Real Matrix */
for (i = 0; i < DEVmaxnum; i++)
if (DEVices[i] && DEVices[i]->DEVbindCSCComplexToReal && ckt->CKThead[i])
DEVices[i]->DEVbindCSCComplexToReal(ckt->CKThead[i], ckt);
ckt->CKTmatrix->SMPkluMatrix->KLUmatrixIsComplex = KLUmatrixReal;
}
#endif
#ifdef XSPICE
/* Call EVTop if event-driven instances exist */
if (ckt->evt->counts.num_insts != 0) {
error = EVTop(ckt,
(ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITJCT,
(ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITFLOAT,
ckt->CKTdcMaxIter,
MIF_TRUE);
EVTdump(ckt, IPC_ANAL_DCOP, 0.0);
EVTop_save(ckt, MIF_TRUE, 0.0);
}
else
#endif
// If no event-driven instances, do what SPICE normally does
error = CKTop(ckt,
(ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITJCT,
(ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITFLOAT,
ckt->CKTdcMaxIter);
if (error) {
fprintf(stdout, "\nAC operating point failed -\n");
CKTncDump(ckt);
return(error);
}
ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG;
error = CKTload(ckt);
if (error) {
tfree(data); return(error);
}
}
#ifdef KLU
if (ckt->CKTmatrix->CKTkluMODE)
{
/* Conversion from Real Matrix to Complex Matrix */
for (i = 0; i < DEVmaxnum; i++)
if (DEVices[i] && DEVices[i]->DEVbindCSCComplex && ckt->CKThead[i])
DEVices[i]->DEVbindCSCComplex(ckt->CKThead[i], ckt);
ckt->CKTmatrix->SMPkluMatrix->KLUmatrixIsComplex = KLUMatrixComplex;
}
#endif
// Store previous rhs
if (rhswoPorts == NULL)
rhswoPorts = (double*)TMALLOC(double, ckt->CKTmaxEqNum);
else
rhswoPorts = (double*)TREALLOC(double, rhswoPorts, ckt->CKTmaxEqNum);
if (rhswoPorts == NULL) {
tfree(data); return (E_NOMEM);
}
if (irhswoPorts == NULL)
irhswoPorts = (double*)TMALLOC(double, ckt->CKTmaxEqNum);
else
irhswoPorts = (double*)TREALLOC(double, irhswoPorts, ckt->CKTmaxEqNum);
if (irhswoPorts == NULL) {
tfree(rhswoPorts);
tfree(data); return (E_NOMEM);
}
ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODESP;
// Let's sweep thru all available ports to build Y matrix
// Y_ij = I_i / V_j | V_k!=j = 0
// (we have only to modify rhs)
int vsrcLookupType = CKTtypelook("Vsource");
int vsrcRoot = -1;
// Get VSRCs root model
if (ckt->CKTVSRCid == -1)
{
for (i = 0; i < DEVmaxnum; i++) {
if (DEVices[i] && DEVices[i]->DEVacLoad && ckt->CKThead[i] && ckt->CKThead[i]->GENmodType == vsrcLookupType) {
vsrcRoot = i;
break;
}
}
if (vsrcRoot == -1)
return (E_NOMOD);
ckt->CKTVSRCid = vsrcRoot;
#if (0)
// Now that we have found the model, we may init the Zref and Gn ports
VSRCspinit(ckt->CKThead[vsrcRoot], ckt, zref, gn, gninv);
#endif
}
else
vsrcRoot = ckt->CKTVSRCid;
// Pre-load everything but RF Ports (these will be updated in the next cycle).
error = NIspPreload(ckt);
if (error) return (error);
// error = VSRCsaveNPData(ckt->CKThead[vsrcRoot]);
// if (error) return (error);
//Keep a backup copy
memcpy(rhswoPorts, ckt->CKTrhs, (size_t)ckt->CKTmaxEqNum * sizeof(double));
memcpy(rhswoPorts, ckt->CKTirhs, (size_t)ckt->CKTmaxEqNum * sizeof(double));
for (activePort = 1; activePort <= ckt->CKTportCount; activePort++)
{
// Copy the backup RHS into CKT's RHS
memcpy(ckt->CKTrhs, rhswoPorts, (size_t)ckt->CKTmaxEqNum * sizeof(double));
memcpy(ckt->CKTirhs, irhswoPorts, (size_t)ckt->CKTmaxEqNum * sizeof(double));
ckt->CKTactivePort = activePort;
// Update only VSRCs
error = VSRCspupdate(ckt->CKThead[vsrcRoot], ckt);
if (error)
{
tfree(rhswoPorts);
tfree(irhswoPorts);
tfree(data);
deleteSPmatrix(ckt);
return(error);
}
error = NIspSolve(ckt);
if (error) {
tfree(rhswoPorts);
tfree(irhswoPorts);
tfree(data);
deleteSPmatrix(ckt);
UPDATE_STATS(DOING_AC);
return(error);
}
#ifdef WANT_SENSE2
if (ckt->CKTsenInfo && (ckt->CKTsenInfo->SENmode & ACSEN)) {
long save;
int save1;
save = ckt->CKTmode;
ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG;
save1 = ckt->CKTsenInfo->SENmode;
ckt->CKTsenInfo->SENmode = ACSEN;
if (freq == job->SPstartFreq) {
ckt->CKTsenInfo->SENacpertflag = 1;
}
else {
ckt->CKTsenInfo->SENacpertflag = 0;
}
error = CKTsenAC(ckt);
if (error)
return (error);
ckt->CKTmode = save;
ckt->CKTsenInfo->SENmode = save1;
}
#endif
// We have done 1 activated port.
error = CKTspCalcPowerWave(ckt);
} //active ports cycle
#ifdef TRACE
printf("HBan: CKTAmat\n");
showcmat(ckt->CKTAmat);
printf("HBan: CKTBmat\n");
showcmat(ckt->CKTBmat);
#endif
// Now we can calculate the full S-Matrix
CKTspCalcSMatrix(ckt);
/*
* Now go with noise cycle, if required
*/
if (job->SPdoNoise)
{
data->freq = freq;
cinit(ckt->CKTNoiseCYmat, 0.0, 0.0);
for (activePort = 0; activePort < ckt->CKTportCount; activePort++)
{
/* the frequency will NOT be stored in array[0] as before; instead,
* it will be given in refVal.rValue (see later)
*/
ckt->CKTactivePort = activePort + 1;
NInspIter(ckt, (VSRCinstance*)(ckt->CKTrfPorts[activePort])); /* solve the adjoint system */
/* put the solution of the current adjoint system into the storage matrix*/
int j;
for (j = 0; j < ckt->CKTmaxEqNum; j++)
{
cplx temp;
temp.re = ckt->CKTrhs[j];
temp.im = ckt->CKTirhs[j];
ckt->CKTadjointRHS->d[activePort][j] = temp;
}
}
/*
now we have all the solutions of the adjoint system, we may look into actual
noise sourches
*/
error = CKTspnoise(ckt, N_DENS, N_CALC, data, internalNoiseAN);
if (error)
{
tfree(internalNoiseAN);
tfree(data);
tfree(rhswoPorts);
tfree(irhswoPorts);
deleteSPmatrix(ckt);
return(error);
}
data->lstFreq = freq;
}
error = CKTspDump(ckt, freq, spPlot, job->SPdoNoise);
if (error) {
UPDATE_STATS(DOING_AC);
tfree(internalNoiseAN);
tfree(rhswoPorts);
tfree(irhswoPorts);
tfree(data);
deleteSPmatrix(ckt);
return(error);
}
/* increment frequency */
switch (job->SPstepType) {
case DECADE:
case OCTAVE:
/* inserted again 14.12.2001 */
#ifdef HAS_PROGREP
{
double endfreq = job->SPstopFreq;
double startfreq = job->SPstartFreq;
endfreq = log(endfreq);
if (startfreq == 0.0)
startfreq = 1e-12;
startfreq = log(startfreq);
if (freq > 0.0)
SetAnalyse("sp", (int)((log(freq) - startfreq) * 1000.0 / (endfreq - startfreq)));
}
#endif
freq *= job->SPfreqDelta;
if (job->SPfreqDelta == 1) goto endsweep;
break;
case LINEAR:
#ifdef HAS_PROGREP
{
double endfreq = job->SPstopFreq;
double startfreq = job->SPstartFreq;
SetAnalyse("sp", (int)((freq - startfreq) * 1000.0 / (endfreq - startfreq)));
}
#endif
freq += job->SPfreqDelta;
if (job->SPfreqDelta == 0) goto endsweep;
break;
default:
tfree(internalNoiseAN);
tfree(rhswoPorts);
tfree(irhswoPorts);
tfree(data);
deleteSPmatrix(ckt);
return(E_INTERN);
}
}
endsweep:
SPfrontEnd->OUTendPlot(spPlot);
spPlot = NULL;
UPDATE_STATS(0);
tfree(internalNoiseAN);
tfree(rhswoPorts);
tfree(irhswoPorts);
deleteSPmatrix(ckt);
tfree(data);
#ifdef KLU
if (ckt->CKTmatrix->CKTkluMODE)
{
/* Conversion from Complex Matrix to Real Matrix */
for (i = 0; i < DEVmaxnum; i++)
if (DEVices[i] && DEVices[i]->DEVbindCSCComplexToReal && ckt->CKThead[i])
DEVices[i]->DEVbindCSCComplexToReal(ckt->CKThead[i], ckt);
ckt->CKTmatrix->SMPkluMatrix->KLUmatrixIsComplex = KLUmatrixReal;
}
#endif
return(0);
}
#endif

74
src/spicelib/analysis/hbaskq.c Normal file
View File

@ -0,0 +1,74 @@
/**********
Copyright 1990 Regents of the University of California. All rights reserved.
Author: 1985 Thomas L. Quarles
**********/
/*
*/
#include "ngspice/ngspice.h"
#include "ngspice/ifsim.h"
#include "ngspice/iferrmsg.h"
#include "ngspice/hbardefs.h"
#include "ngspice/cktdefs.h"
#ifdef WITH_HB
/* ARGSUSED */
int
HBaskQuest(CKTcircuit *ckt, JOB *anal, int which, IFvalue *value)
{
HBAN *job = (HBAN *) anal;
NG_IGNORE(ckt);
switch(which) {
case SP_START:
value->rValue = job->SPstartFreq;
break;
case SP_STOP:
value->rValue = job->SPstopFreq ;
break;
case SP_STEPS:
value->iValue = job->SPnumberSteps;
break;
case SP_DEC:
if (job->SPstepType == DECADE) {
value->iValue=1;
} else {
value->iValue=0;
}
break;
case SP_OCT:
if (job->SPstepType == OCTAVE) {
value->iValue=1;
} else {
value->iValue=0;
}
break;
case SP_LIN:
if (job->SPstepType == LINEAR) {
value->iValue=1;
} else {
value->iValue=0;
}
break;
case SP_DONOISE:
if (job->SPdoNoise)
value->iValue = 1;
else
value->iValue = 0;
break;
default:
return(E_BADPARM);
}
return(OK);
}
#endif

112
src/spicelib/analysis/hbsetp.c Normal file
View File

@ -0,0 +1,112 @@
/**********
Copyright 1990 Regents of the University of California. All rights reserved.
Author: 1985 Thomas L. Quarles
**********/
#include "ngspice/ngspice.h"
#include "ngspice/ifsim.h"
#include "ngspice/iferrmsg.h"
#include "ngspice/hbardefs.h"
#include "ngspice/cktdefs.h"
#include "analysis.h"
#ifdef RFSPICE
/* ARGSUSED */
int
HBsetParm(CKTcircuit *ckt, JOB *anal, int which, IFvalue *value)
{
HBAN *job = (HBAN *) anal;
NG_IGNORE(ckt);
switch(which) {
case SP_START:
if (value->rValue < 0.0) {
errMsg = copy("Frequency of < 0 is invalid for AC start");
job->SPstartFreq = 1.0;
return(E_PARMVAL);
}
job->SPstartFreq = value->rValue;
break;
case SP_STOP:
if (value->rValue < 0.0) {
errMsg = copy("Frequency of < 0 is invalid for AC stop");
job->SPstartFreq = 1.0;
return(E_PARMVAL);
}
job->SPstopFreq = value->rValue;
break;
case SP_STEPS:
job->SPnumberSteps = value->iValue;
break;
case SP_DEC:
if(value->iValue) {
job->SPstepType = DECADE;
} else {
if (job->SPstepType == DECADE) {
job->SPstepType = 0;
}
}
break;
case SP_OCT:
if(value->iValue) {
job->SPstepType = OCTAVE;
} else {
if (job->SPstepType == OCTAVE) {
job->SPstepType = 0;
}
}
break;
case SP_LIN:
if(value->iValue) {
job->SPstepType = LINEAR;
} else {
if (job->SPstepType == LINEAR) {
job->SPstepType = 0;
}
}
break;
case SP_DONOISE:
job->SPdoNoise = value->iValue == 1;
break;
default:
return(E_BADPARM);
}
return(OK);
}
static IFparm HBparms[] = {
{ "f1", SP_START, IF_SET|IF_ASK|IF_REAL, "fundamental frequency" },
{ "f2", SP_STOP, IF_SET|IF_ASK|IF_REAL, "second frequency" }
};
SPICEanalysis HBinfo = {
{
"HB",
"Harmonic Balance analysis",
NUMELEMS(HBparms),
HBparms
},
sizeof(HBAN),
FREQUENCYDOMAIN,
1,
HBsetParm,
HBaskQuest,
NULL,
HBan
};
#endif

46
src/spicelib/parser/inp2dot.c
View File

@ -742,61 +742,35 @@ dot_sp(char* line, void* ckt, INPtables* tab, struct card* current,
}
#ifdef WITH_HB
/*SP: Steady State Analyis */
/*HB: Harmonic Balance Analyis */
static int
dot_hb(char* line, void* ckt, INPtables* tab, struct card* current,
void* task, void* gnode, JOB* foo)
{
int error; /* error code temporary */
IFvalue ptemp; /* a value structure to package resistance into */
// IFvalue ptemp; /* a value structure to package resistance into */
IFvalue* parm; /* a pointer to a value struct for function returns */
char* nname; /* the oscNode name */
CKTnode* nnode; /* the oscNode node */
int which; /* which analysis we are performing */
char* word; /* something to stick a word of input into */
NG_IGNORE(gnode);
/* .pss Fguess StabTime OscNode <UIC>*/
which = ft_find_analysis("PSS");
/* .hb frequ */
which = ft_find_analysis("HB");
if (which == -1) {
LITERR("Periodic steady state analysis unsupported.\n");
LITERR("Harmonic Balance analysis unsupported.\n");
return (0);
}
IFC(newAnalysis, (ckt, which, "Harmonic Balance State Analysis", &foo, task));
IFC(newAnalysis, (ckt, which, "Harmonic Balance Analysis", &foo, task));
parm = INPgetValue(ckt, &line, IF_REALVEC, tab); /* Fguess */
GCA(INPapName, (ckt, which, foo, "freq", parm));
GCA(INPapName, (ckt, which, foo, "freq1", parm));
parm = INPgetValue(ckt, &line, IF_INTVEC, tab); /* StabTime */
GCA(INPapName, (ckt, which, foo, "harmonics", parm));
INPgetNetTok(&line, &nname, 0);
INPtermInsert(ckt, &nname, tab, &nnode);
ptemp.nValue = nnode;
GCA(INPapName, (ckt, which, foo, "oscnode", &ptemp)); /* OscNode given as string */
parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* PSS points */
GCA(INPapName, (ckt, which, foo, "points", parm));
parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* PSS harmonics */
GCA(INPapName, (ckt, which, foo, "harmonics", parm));
parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* SC iterations */
GCA(INPapName, (ckt, which, foo, "sc_iter", parm));
parm = INPgetValue(ckt, &line, IF_REAL, tab); /* Steady coefficient */
GCA(INPapName, (ckt, which, foo, "steady_coeff", parm));
parm = INPgetValue(ckt, &line, IF_REALVEC, tab); /* Fguess */
GCA(INPapName, (ckt, which, foo, "freq2", parm));
if (*line) {
INPgetTok(&line, &word, 1); /* uic? */
if (strcmp(word, "uic") == 0) {
ptemp.iValue = 1;
GCA(INPapName, (ckt, which, foo, "uic", &ptemp));
}
else {
fprintf(stderr, "Error: unknown parameter %s on .pss - ignored\n", word);
}
fprintf(stderr, "Error: unknown parameter %s on .hb - ignored\n", word);
}
return (0);
}

3
visualc/src/include/ngspice/config.h
View File

@ -11,6 +11,9 @@
/* Define if you want PSS analysis */
#define WITH_PSS /**/
/* Define if you want Harmonic Balance analysis */
#define WITH_HB /**/
/* Name of package */
#define PACKAGE "ngspice"

4
visualc/vngspice.vcxproj
View File

@ -975,6 +975,7 @@
<ClInclude Include="..\src\include\ngspice\cmconstants.h" />
<ClInclude Include="..\src\include\ngspice\cmproto.h" />
<ClInclude Include="..\src\include\ngspice\cmtypes.h" />
<ClInclude Include="..\src\include\ngspice\hbardefs.h" />
<ClInclude Include="..\src\include\ngspice\osdiitf.h" />
<ClInclude Include="..\src\include\ngspice\logicexp.h" />
<ClInclude Include="..\src\include\ngspice\spardefs.h" />
@ -1789,6 +1790,9 @@
<ClCompile Include="..\src\spicelib\analysis\dkerproc.c" />
<ClCompile Include="..\src\spicelib\analysis\dloadfns.c" />
<ClCompile Include="..\src\spicelib\analysis\dsetparm.c" />
<ClCompile Include="..\src\spicelib\analysis\hban.c" />
<ClCompile Include="..\src\spicelib\analysis\hbaskq.c" />
<ClCompile Include="..\src\spicelib\analysis\hbsetp.c" />
<ClCompile Include="..\src\spicelib\analysis\naskq.c" />
<ClCompile Include="..\src\spicelib\analysis\nevalsrc.c" />
<ClCompile Include="..\src\spicelib\analysis\ninteg.c" />