From abd5b5ea0424ab9505ba9da5b81904fa14f02b59 Mon Sep 17 00:00:00 2001 From: Alessio Cacciatori Date: Tue, 29 Mar 2022 09:06:10 +0200 Subject: [PATCH] S-parameters in ngspice With this commit the patch provided by Alessio Cacchiatori the S-parameter is completed: Noise simulation added with C matrix output Y and Z matrix output enabled To allow compiling with gcc, the dense.h inline functions have been put into denseinlines.h --- src/frontend/commands.c | 5 +- src/frontend/runcoms.c | 2 +- src/include/ngspice/cktdefs.h | 7 +- src/include/ngspice/noisedef.h | 20 +- src/include/ngspice/{spdefs.h => spardefs.h} | 0 src/include/ngspice/tskdefs.h | 3 + src/maths/dense/dense.c | 119 +--- src/maths/dense/dense.h | 291 ++++++-- src/maths/dense/denseinlines.h | 144 ++++ src/maths/ni/niaciter.c | 1 - src/spicelib/analysis/Makefile.am | 2 + src/spicelib/analysis/cktdest.c | 2 + src/spicelib/analysis/cktspdum.c | 117 +++- src/spicelib/analysis/cktspnoise.c | 154 +++++ src/spicelib/analysis/nevalsrc.c | 180 ++++- src/spicelib/analysis/span.c | 681 ++++++++++++------- src/spicelib/analysis/spaskq.c | 2 +- src/spicelib/analysis/spsetp.c | 2 +- src/spicelib/devices/cktinit.c | 3 + src/spicelib/devices/vsrc/vsrcacld.c | 40 +- src/spicelib/devices/vsrc/vsrcext.h | 5 +- src/spicelib/devices/vsrc/vsrcload.c | 4 +- src/spicelib/devices/vsrc/vsrctemp.c | 21 +- src/spicelib/parser/inp2dot.c | 5 +- visualc/vngspice.vcxproj | 5 +- 25 files changed, 1367 insertions(+), 448 deletions(-) rename src/include/ngspice/{spdefs.h => spardefs.h} (100%) create mode 100644 src/maths/dense/denseinlines.h create mode 100644 src/spicelib/analysis/cktspnoise.c diff --git a/src/frontend/commands.c b/src/frontend/commands.c index c18ef7fc8..c6754ef08 100644 --- a/src/frontend/commands.c +++ b/src/frontend/commands.c @@ -307,12 +307,11 @@ struct comm spcp_coms[] = { /* SP */ #endif #ifdef RFSPICE -/* SP: S parameter Analysis */ +/* S parameter Analysis */ { "sp", com_sp, TRUE, TRUE, { 0, 0, 0, 0 }, E_DEFHMASK, 0, LOTS, NULL, - "[.sp line args] : Do a S parameter analysis." }, - /* SP */ + "[.sp line args] : Do an S-parameter analysis." }, #endif { "ac", com_ac, TRUE, TRUE, { 0, 0, 0, 0 }, E_DEFHMASK, 0, LOTS, diff --git a/src/frontend/runcoms.c b/src/frontend/runcoms.c index c50cbb0b8..8ea318bf9 100644 --- a/src/frontend/runcoms.c +++ b/src/frontend/runcoms.c @@ -186,7 +186,7 @@ com_pss(wordlist *wl) #endif #ifdef RFSPICE -/* S parameter Analysis*/ +/* S-parameter Analysis*/ void com_sp(wordlist* wl) { diff --git a/src/include/ngspice/cktdefs.h b/src/include/ngspice/cktdefs.h index 95d2ec6a8..03fe07fe4 100644 --- a/src/include/ngspice/cktdefs.h +++ b/src/include/ngspice/cktdefs.h @@ -292,6 +292,11 @@ struct CKTcircuit { CMat* CKTSmat; CMat* CKTYmat; CMat* CKTZmat; + // Data for RF Noise Calculations + double* CKTportY; + CMat* CKTNoiseCYmat; + unsigned int CKTnoiseSourceCount; + CMat* CKTadjointRHS; // Matrix where Znj are stored. Znj = impedance from j-th noise source to n-th port #endif #ifdef WITH_PSS /* SP: Periodic Steady State Analysis - 100609 */ @@ -455,7 +460,7 @@ extern int DCpss(CKTcircuit *, int); extern int SPan(CKTcircuit*, int); extern int SPaskQuest(CKTcircuit*, JOB*, int, IFvalue*); extern int SPsetParm(CKTcircuit*, JOB*, int, IFvalue*); -extern int CKTspDump(CKTcircuit*, double, runDesc*); +extern int CKTspDump(CKTcircuit*, double, runDesc*, unsigned int); extern int CKTspLoad(CKTcircuit*); extern unsigned int CKTmatrixIndex(CKTcircuit*, unsigned int, unsigned int); extern int CKTspCalcPowerWave(CKTcircuit* ckt); diff --git a/src/include/ngspice/noisedef.h b/src/include/ngspice/noisedef.h index 49821d024..6140a504a 100644 --- a/src/include/ngspice/noisedef.h +++ b/src/include/ngspice/noisedef.h @@ -118,7 +118,25 @@ typedef struct { /* misc constants */ +#ifdef RFSPICE +#define NOISE_ADD_OUTVAR(ckt, data, fmt, aname, bname) \ + if (ckt->CKTcurrentAnalysis & DOING_SP) { \ + ckt->CKTnoiseSourceCount++; \ + } \ + else\ + do { \ + data->namelist = TREALLOC(IFuid, data->namelist, data->numPlots + 1); \ + if (!data->namelist) \ + return E_NOMEM; \ + char *name = tprintf(fmt, aname, bname); \ + if (!name) \ + return E_NOMEM; \ + SPfrontEnd->IFnewUid(ckt, &(data->namelist[data->numPlots++]), \ + NULL, name, UID_OTHER, NULL); \ + tfree(name); \ + } while(0) +#else #define NOISE_ADD_OUTVAR(ckt, data, fmt, aname, bname) \ do { \ data->namelist = TREALLOC(IFuid, data->namelist, data->numPlots + 1); \ @@ -131,7 +149,7 @@ typedef struct { NULL, name, UID_OTHER, NULL); \ tfree(name); \ } while(0) - +#endif void NevalSrc (double *noise, double *lnNoise, CKTcircuit *ckt, int type, int node1, int node2, double param); void NevalSrc2 (double *, double *, CKTcircuit *, int, int, int, double, int, int, double, double); diff --git a/src/include/ngspice/spdefs.h b/src/include/ngspice/spardefs.h similarity index 100% rename from src/include/ngspice/spdefs.h rename to src/include/ngspice/spardefs.h diff --git a/src/include/ngspice/tskdefs.h b/src/include/ngspice/tskdefs.h index 43bbd784f..e5167dce0 100644 --- a/src/include/ngspice/tskdefs.h +++ b/src/include/ngspice/tskdefs.h @@ -29,6 +29,9 @@ struct TSKtask { #define DOING_TRCV 2 #define DOING_AC 4 #define DOING_TRAN 8 +#ifdef RFSPICE +#define DOING_SP 16 +#endif int TSKbypass; int TSKdcMaxIter; /* iteration limit for dc op. (itl1) */ diff --git a/src/maths/dense/dense.c b/src/maths/dense/dense.c index 097da80ad..b9a33e147 100644 --- a/src/maths/dense/dense.c +++ b/src/maths/dense/dense.c @@ -1,5 +1,6 @@ #include "dense.h" +#include "denseinlines.h" #include "ngspice/ngspice.h" #include #include @@ -7,124 +8,6 @@ #include "ngspice/bool.h" #include "ngspice/iferrmsg.h" -inline double cmod(cplx a); - -inline void setcplx(cplx* d, double r, double i) -{ - d->re = r; d->im = i; -} - -inline void cmultc(cplx* res, cplx a, cplx b) -{ - res->re = a.re * b.re - a.im * b.im; - res->im = a.im * b.re + a.re * b.im; -} - -inline cplx cmultco(cplx a, cplx b) -{ - cplx res; - res.re = a.re * b.re - a.im * b.im; - res.im = a.im * b.re + a.re * b.im; - return res; -} - - - -inline cplx cdivco(cplx a, cplx b) -{ - cplx res; - double dmod = cmod(b); - dmod *= dmod; - res.re = (a.re * b.re + a.im * b.im) / dmod; - res.im = (a.im * b.re - a.re * b.im) / dmod; - return res; -} - - -inline cplx cmultdo(cplx a, double d) -{ - cplx res; - res.re = a.re * d; - res.im = a.im * d; - return res; -} -inline void cmultd(cplx* res, cplx a, double d) -{ - res->re = a.re * d; - res->im = a.im * d; -} - -inline void caddc(cplx* res, cplx a, cplx b) -{ - res->re = a.re + b.re; - res->im = a.im + b.im; -} - -inline cplx caddco( cplx a, cplx b) -{ - cplx res; - res.re = a.re + b.re; - res.im = a.im + b.im; - return res; -} - - -inline void caddd(cplx* res, cplx a, double d) -{ - res->re = a.re + d; -} - -inline void csubc(cplx* res, cplx a, cplx b) -{ - res->re = a.re - b.re; - res->im = a.im - b.im; -} - -inline cplx csubco(cplx a, cplx b) -{ - cplx res; - res.re = a.re - b.re; - res.im = a.im - b.im; - return res; -} - -inline void csubd(cplx* res, cplx a, double d) -{ - res->re = a.re - d; -} - -inline double cmodinv(cplx a) -{ - return 1.0 / cmod(a); -} - -inline double cmod(cplx a) -{ - return (a.re * a.re + a.im * a.im); -} - -inline int ciszero(cplx a) -{ - return (a.re == 0) && (a.im == 0); -} -inline cplx cinv(cplx a) -{ - cplx res; - double cpmod = cmodinv(a); - res.re = a.re * cpmod; - res.im = -a.im * cpmod; - return res; -} - -inline cplx conju(cplx a) -{ - cplx res; - res.re = a.re; - res.im = -a.im; - return res; -} - - void showmat(Mat* A) { if (A->row > 0 && A->col > 0) { diff --git a/src/maths/dense/dense.h b/src/maths/dense/dense.h index d0dd5dad2..49b0c3fb3 100644 --- a/src/maths/dense/dense.h +++ b/src/maths/dense/dense.h @@ -1,6 +1,7 @@ #ifndef ngspice_DENSE_MATRIX_H #define ngspice_DENSE_MATRIX_H +#include typedef struct cplx { @@ -26,78 +27,222 @@ typedef struct MatList{ struct MatList* next; }MatList; -extern void showmat(Mat* A); -extern void showcmat(CMat* A); -extern CMat* newcmat(int r, int c, double dr, double di); -extern CMat* newcmatnoinit(int r, int c); -extern Mat* newmat(int r, int c, double d); -extern Mat* newmatnoinit(int r, int c); -extern void freecmat(CMat* A); -extern void freemat(Mat* A); -extern CMat* ceye(int n); -extern Mat* eye(int n); -extern Mat* zeros(int r, int c); -extern CMat* czeros(int r, int c); -extern Mat* ones(int r, int c); -extern CMat* cones(int r, int c); -extern Mat* randm(int r, int c, double l, double u); -extern CMat* randcm(int r, int c, double l, double u); -extern double get(Mat* M, int r, int c); -extern cplx getcplx(CMat* M, int r, int c); -extern void set(Mat* M, int r, int c, double d); -extern void setc(CMat* M, int r, int c, cplx d); -extern Mat* scalarmultiply(Mat* M, double c); -extern CMat* cscalarmultiply(CMat* M, double c); -extern CMat* complexmultiply(CMat* M, cplx c); -extern Mat* sum(Mat* A, Mat* B); -extern CMat* csum(CMat* A, CMat* B); -extern Mat* minus(Mat* A, Mat* B); -extern CMat* cminus(CMat* A, CMat* B); -extern Mat* submat(Mat* A, int r1, int r2, int c1, int c2); -extern void submat2(Mat* A, Mat* B, int r1, int r2, int c1, int c2); -extern CMat* subcmat(CMat* A, int r1, int r2, int c1, int c2); -extern void subcmat2(CMat* A, CMat* B, int r1, int r2, int c1, int c2); -extern Mat* multiply(Mat* A, Mat* B); -extern CMat* cmultiply(CMat* A, CMat* B); -extern Mat* removerow(Mat* A, int r); -extern Mat* removecol(Mat* A, int c); -extern void removerow2(Mat* A, Mat* B, int r); -extern void removecol2(Mat* A, Mat* B, int c); -extern CMat* cremoverow(CMat* A, int r); -extern CMat* cremovecol(CMat* A, int c); -extern void cremoverow2(CMat* A, CMat* B, int r); -extern void cremovecol2(CMat* A, CMat* B, int c); -extern Mat* transpose(Mat* A); -extern CMat* ctranspose(CMat* A); -extern double trace(Mat* A); -extern cplx ctrace(CMat* A); -extern Mat* adjoint(Mat* A); -extern CMat* cadjoint(CMat* A); -extern CMat* ctransposeconj(CMat* source); -extern Mat* inverse(Mat* A); -extern CMat* cinverse(CMat* A); -extern Mat* copyvalue(Mat* A); -extern CMat* copycvalue(CMat* A); -extern Mat* triinverse(Mat* A); -extern CMat* ctriinverse(CMat* A); -extern Mat* rowechelon(Mat* A); -extern CMat* crowechelon(CMat* A); -extern Mat* hconcat(Mat* A, Mat* B); -extern CMat* chconcat(CMat* A, CMat* B); -extern Mat* vconcat(Mat* A, Mat* B); -extern CMat* cvconcat(CMat* A, CMat* B); -extern double norm(Mat* A); -extern double cnorm(CMat* A); -extern Mat* nullmat(Mat* A); -extern MatList* lu(Mat* A); -extern double innermultiply(Mat* a, Mat* b); -extern MatList* qr(Mat* A); -extern int complexmultiplydest(CMat* M, cplx c, CMat* dest); -extern void cinversedest(CMat* A, CMat* dest); -extern int copycvaluedest(CMat* A, CMat* dest); -extern int cmultiplydest(CMat* A, CMat* B, CMat* dest); -extern void init(Mat* A, double d); -extern void cinit(CMat* A, double dr, double di); -extern void resizemat(Mat* A, int r, int c); + void showmat(Mat* A); + void showcmat(CMat* A); + CMat* newcmat(int r, int c, double dr, double di); + CMat* newcmatnoinit(int r, int c); + Mat* newmat(int r, int c, double d); + Mat* newmatnoinit(int r, int c); + void freecmat(CMat* A); + void freemat(Mat* A); + CMat* ceye(int n); + Mat* eye(int n); + Mat* zeros(int r, int c); + CMat* czeros(int r, int c); + Mat* ones(int r, int c); + CMat* cones(int r, int c); + Mat* randm(int r, int c, double l, double u); + CMat* randcm(int r, int c, double l, double u); + double get(Mat* M, int r, int c); + cplx getcplx(CMat* M, int r, int c); + void set(Mat* M, int r, int c, double d); + void setc(CMat* M, int r, int c, cplx d); + Mat* scalarmultiply(Mat* M, double c); + CMat* cscalarmultiply(CMat* M, double c); + CMat* complexmultiply(CMat* M, cplx c); + Mat* sum(Mat* A, Mat* B); + CMat* csum(CMat* A, CMat* B); + Mat* minus(Mat* A, Mat* B); + CMat* cminus(CMat* A, CMat* B); + Mat* submat(Mat* A, int r1, int r2, int c1, int c2); + void submat2(Mat* A, Mat* B, int r1, int r2, int c1, int c2); + CMat* subcmat(CMat* A, int r1, int r2, int c1, int c2); + void subcmat2(CMat* A, CMat* B, int r1, int r2, int c1, int c2); + Mat* multiply(Mat* A, Mat* B); + CMat* cmultiply(CMat* A, CMat* B); + Mat* removerow(Mat* A, int r); + Mat* removecol(Mat* A, int c); + void removerow2(Mat* A, Mat* B, int r); + void removecol2(Mat* A, Mat* B, int c); + CMat* cremoverow(CMat* A, int r); + CMat* cremovecol(CMat* A, int c); + void cremoverow2(CMat* A, CMat* B, int r); + void cremovecol2(CMat* A, CMat* B, int c); + Mat* transpose(Mat* A); + CMat* ctranspose(CMat* A); + double trace(Mat* A); + cplx ctrace(CMat* A); + Mat* adjoint(Mat* A); + CMat* cadjoint(CMat* A); + CMat* ctransposeconj(CMat* source); + Mat* inverse(Mat* A); + CMat* cinverse(CMat* A); + Mat* copyvalue(Mat* A); + CMat* copycvalue(CMat* A); + Mat* copyvalue(Mat* A); + CMat* copycvalue(CMat* A); + Mat* triinverse(Mat* A); + CMat* ctriinverse(CMat* A); + Mat* rowechelon(Mat* A); + CMat* crowechelon(CMat* A); + Mat* hconcat(Mat* A, Mat* B); + CMat* chconcat(CMat* A, CMat* B); + Mat* vconcat(Mat* A, Mat* B); + CMat* cvconcat(CMat* A, CMat* B); + double norm(Mat* A); + double cnorm(CMat* A); + Mat* nullmat(Mat* A); + MatList* lu(Mat* A); + double innermultiply(Mat* a, Mat* b); + MatList* qr(Mat* A); + int complexmultiplydest(CMat* M, cplx c, CMat* dest); + void cinversedest(CMat* A, CMat* dest); + int copycvaluedest(CMat* A, CMat* dest); + int cmultiplydest(CMat* A, CMat* B, CMat* dest); + void init(Mat* A, double d); + void cinit(CMat* A, double dr, double di); + void resizemat(Mat* A, int r, int c); + /* + inline void setcplx(cplx* d, double r, double i); + inline void cmultc(cplx* res, cplx a, cplx b); + inline cplx cmultco(cplx a, cplx b); + inline cplx cmultdo(cplx a, double d); + inline void cmultd(cplx* res, cplx a, double d); + inline void caddc(cplx* res, cplx a, cplx b); + inline cplx caddco(cplx a, cplx b); + inline void caddd(cplx* res, cplx a, double d); + inline void csubc(cplx* res, cplx a, cplx b); + inline cplx csubco(cplx a, cplx b); + inline void csubd(cplx* res, cplx a, double d); + inline double cmodinv(cplx a); + inline double cmodsqr(cplx a); + + inline int ciszero(cplx a); + inline cplx cinv(cplx a); + inline cplx conju(cplx a); + + extern inline double cmodu(cplx a); + extern inline cplx cdivco(cplx a, cplx b); + + inline void setcplx(cplx* d, double r, double i) + { + d->re = r; d->im = i; + } + + inline void cmultc(cplx* res, cplx a, cplx b) + { + res->re = a.re * b.re - a.im * b.im; + res->im = a.im * b.re + a.re * b.im; + } + + inline cplx cmultco(cplx a, cplx b) + { + cplx res; + res.re = a.re * b.re - a.im * b.im; + res.im = a.im * b.re + a.re * b.im; + return res; + } + + + + inline cplx cdivco(cplx a, cplx b) + { + cplx res; + double dmod = cmodinv(b); + + res.re = (a.re * b.re + a.im * b.im) * dmod; + res.im = (a.im * b.re - a.re * b.im) * dmod; + return res; + } + + + inline cplx cmultdo(cplx a, double d) + { + cplx res; + res.re = a.re * d; + res.im = a.im * d; + return res; + } + inline void cmultd(cplx* res, cplx a, double d) + { + res->re = a.re * d; + res->im = a.im * d; + } + + inline void caddc(cplx* res, cplx a, cplx b) + { + res->re = a.re + b.re; + res->im = a.im + b.im; + } + + inline cplx caddco(cplx a, cplx b) + { + cplx res; + res.re = a.re + b.re; + res.im = a.im + b.im; + return res; + } + + + inline void caddd(cplx* res, cplx a, double d) + { + res->re = a.re + d; + } + + inline void csubc(cplx* res, cplx a, cplx b) + { + res->re = a.re - b.re; + res->im = a.im - b.im; + } + + inline cplx csubco(cplx a, cplx b) + { + cplx res; + res.re = a.re - b.re; + res.im = a.im - b.im; + return res; + } + + inline void csubd(cplx* res, cplx a, double d) + { + res->re = a.re - d; + } + + inline double cmodinv(cplx a) + { + return 1.0 / cmodsqr(a); + } + + inline double cmodsqr(cplx a) + { + return (a.re * a.re + a.im * a.im); + } + + inline double cmodu(cplx a) + { + return sqrt(cmodsqr(a)); + } + inline int ciszero(cplx a) + { + return (a.re == 0) && (a.im == 0); + } + inline cplx cinv(cplx a) + { + cplx res; + double cpmod = cmodinv(a); + res.re = a.re * cpmod; + res.im = -a.im * cpmod; + return res; + } + + inline cplx conju(cplx a) + { + cplx res; + res.re = a.re; + res.im = -a.im; + return res; + } + */ + #endif diff --git a/src/maths/dense/denseinlines.h b/src/maths/dense/denseinlines.h new file mode 100644 index 000000000..23d05f199 --- /dev/null +++ b/src/maths/dense/denseinlines.h @@ -0,0 +1,144 @@ + +#include + + inline void setcplx(cplx* d, double r, double i); + inline void cmultc(cplx* res, cplx a, cplx b); + inline cplx cmultco(cplx a, cplx b); + inline cplx cmultdo(cplx a, double d); + inline void cmultd(cplx* res, cplx a, double d); + inline void caddc(cplx* res, cplx a, cplx b); + inline cplx caddco(cplx a, cplx b); + inline void caddd(cplx* res, cplx a, double d); + inline void csubc(cplx* res, cplx a, cplx b); + inline cplx csubco(cplx a, cplx b); + inline void csubd(cplx* res, cplx a, double d); + inline double cmodinv(cplx a); + inline double cmodsqr(cplx a); + + inline int ciszero(cplx a); + inline cplx cinv(cplx a); + inline cplx conju(cplx a); + + inline double cmodu(cplx a); + inline cplx cdivco(cplx a, cplx b); + + inline void setcplx(cplx* d, double r, double i) + { + d->re = r; d->im = i; + } + + inline void cmultc(cplx* res, cplx a, cplx b) + { + res->re = a.re * b.re - a.im * b.im; + res->im = a.im * b.re + a.re * b.im; + } + + inline cplx cmultco(cplx a, cplx b) + { + cplx res; + res.re = a.re * b.re - a.im * b.im; + res.im = a.im * b.re + a.re * b.im; + return res; + } + + + + inline cplx cdivco(cplx a, cplx b) + { + cplx res; + double dmod = cmodinv(b); + + res.re = (a.re * b.re + a.im * b.im) * dmod; + res.im = (a.im * b.re - a.re * b.im) * dmod; + return res; + } + + + inline cplx cmultdo(cplx a, double d) + { + cplx res; + res.re = a.re * d; + res.im = a.im * d; + return res; + } + inline void cmultd(cplx* res, cplx a, double d) + { + res->re = a.re * d; + res->im = a.im * d; + } + + inline void caddc(cplx* res, cplx a, cplx b) + { + res->re = a.re + b.re; + res->im = a.im + b.im; + } + + inline cplx caddco(cplx a, cplx b) + { + cplx res; + res.re = a.re + b.re; + res.im = a.im + b.im; + return res; + } + + + inline void caddd(cplx* res, cplx a, double d) + { + res->re = a.re + d; + } + + inline void csubc(cplx* res, cplx a, cplx b) + { + res->re = a.re - b.re; + res->im = a.im - b.im; + } + + inline cplx csubco(cplx a, cplx b) + { + cplx res; + res.re = a.re - b.re; + res.im = a.im - b.im; + return res; + } + + inline void csubd(cplx* res, cplx a, double d) + { + res->re = a.re - d; + } + + inline double cmodinv(cplx a) + { + return 1.0 / cmodsqr(a); + } + + inline double cmodsqr(cplx a) + { + return (a.re * a.re + a.im * a.im); + } + + inline double cmodu(cplx a) + { + return sqrt(cmodsqr(a)); + } + inline int ciszero(cplx a) + { + return (a.re == 0) && (a.im == 0); + } + inline cplx cinv(cplx a) + { + cplx res; + double cpmod = cmodinv(a); + res.re = a.re * cpmod; + res.im = -a.im * cpmod; + return res; + } + + inline cplx conju(cplx a) + { + cplx res; + res.re = a.re; + res.im = -a.im; + return res; + } + + diff --git a/src/maths/ni/niaciter.c b/src/maths/ni/niaciter.c index 5905aac56..eccac93a0 100644 --- a/src/maths/ni/niaciter.c +++ b/src/maths/ni/niaciter.c @@ -95,7 +95,6 @@ int NIspSolve(CKTcircuit* ckt) #endif - int NIacIter(CKTcircuit *ckt) { diff --git a/src/spicelib/analysis/Makefile.am b/src/spicelib/analysis/Makefile.am index 62ab58d1a..c6d98c35f 100644 --- a/src/spicelib/analysis/Makefile.am +++ b/src/spicelib/analysis/Makefile.am @@ -62,6 +62,7 @@ libckt_la_SOURCES = \ cktsetup.c \ cktsgen.c \ cktsopt.c \ + cktspnoise.c \ ckttemp.c \ cktterr.c \ ckttroub.c \ @@ -83,6 +84,7 @@ libckt_la_SOURCES = \ nevalsrc.c \ ninteg.c \ noisean.c \ + noisesp.c \ nsetparm.c \ optran.c \ com_optran.h \ diff --git a/src/spicelib/analysis/cktdest.c b/src/spicelib/analysis/cktdest.c index 6c4fb82b4..ca4c67786 100644 --- a/src/spicelib/analysis/cktdest.c +++ b/src/spicelib/analysis/cktdest.c @@ -112,6 +112,8 @@ CKTdestroy(CKTcircuit *ckt) freecmat(ckt->CKTSmat); ckt->CKTSmat = NULL; freecmat(ckt->CKTYmat); ckt->CKTYmat = NULL; freecmat(ckt->CKTZmat); ckt->CKTZmat = NULL; + freecmat(ckt->CKTNoiseCYmat); ckt->CKTNoiseCYmat = NULL; + freecmat(ckt->CKTadjointRHS); ckt->CKTadjointRHS = NULL; #endif FREE(ckt); diff --git a/src/spicelib/analysis/cktspdum.c b/src/spicelib/analysis/cktspdum.c index b8189bbd0..90dddfa78 100644 --- a/src/spicelib/analysis/cktspdum.c +++ b/src/spicelib/analysis/cktspdum.c @@ -15,7 +15,16 @@ Author: 1985 Thomas L. Quarles #include "ngspice/ifsim.h" #include "vsrc/vsrcdefs.h" -#ifdef RFSPICE +extern CMat* eyem; +extern CMat* zref; +extern CMat* gn; +extern CMat* gninv; + +extern double NF; +extern double Rn; +extern cplx Sopt; +extern double Fmin; + unsigned int CKTmatrixIndex(CKTcircuit* ckt, unsigned int source, unsigned int dest) { @@ -28,6 +37,25 @@ int CKTspCalcSMatrix(CKTcircuit* ckt) if (Ainv == NULL) return (E_NOMEM); cmultiplydest(ckt->CKTBmat, Ainv, ckt->CKTSmat); freecmat(Ainv); + // Calculate Y matrix + CMat* temp = cmultiply(ckt->CKTSmat, zref); + CMat* temp2 = csum(temp, zref); + CMat* temp3 = cmultiply(temp2, gn); + + CMat* temp4 = cminus(eyem, ckt->CKTSmat); + CMat* temp5 = cinverse(temp4); + + cmultiplydest(temp5, temp3, temp); + cmultiplydest(gninv, temp, ckt->CKTZmat); + + cinversedest(ckt->CKTZmat, ckt->CKTYmat); + + freecmat(temp); + freecmat(temp2); + freecmat(temp3); + freecmat(temp4); + freecmat(temp5); + return (OK); } @@ -64,8 +92,9 @@ int CKTspCalcPowerWave(CKTcircuit* ckt) return (OK); } + int -CKTspDump(CKTcircuit *ckt, double freq, runDesc *plot) +CKTspDump(CKTcircuit *ckt, double freq, runDesc *plot, unsigned int doNoise) { double *rhsold; double *irhsold; @@ -77,10 +106,20 @@ CKTspDump(CKTcircuit *ckt, double freq, runDesc *plot) rhsold = ckt->CKTrhsOld; irhsold = ckt->CKTirhsOld; freqData.rValue = freq; - unsigned int extraSPdataCount = ckt->CKTportCount * ckt->CKTportCount; + unsigned int extraSPdataCount = 3* ckt->CKTportCount * ckt->CKTportCount; valueData.v.numValue = ckt->CKTmaxEqNum - 1 + extraSPdataCount; - data = TMALLOC(IFcomplex, ckt->CKTmaxEqNum - 1 + extraSPdataCount); + unsigned int datasize = ckt->CKTmaxEqNum - 1 + extraSPdataCount; + + // Add Cy matrix, NF, Rn, SOpt, NFmin + if (doNoise) + { + datasize += ckt->CKTportCount * ckt->CKTportCount; + if (ckt->CKTportCount == 2) // account for NF, Sopt, NFmin, Rn + datasize += 4; + } + + data = TMALLOC(IFcomplex, datasize); valueData.v.vec.cVec = data; for (i=0;iCKTmaxEqNum-1;i++) { data[i].real = rhsold[i+1]; @@ -89,18 +128,85 @@ CKTspDump(CKTcircuit *ckt, double freq, runDesc *plot) if (ckt->CKTrfPorts ) { + unsigned int nPlot = ckt->CKTmaxEqNum - 1 ; // Cycle thru all ports for (unsigned int pdest = 0; pdest < ckt->CKTportCount; pdest++) { for (unsigned int psource = 0; psource < ckt->CKTportCount; psource++) { - unsigned int nPlot = ckt->CKTmaxEqNum - 1 + CKTmatrixIndex(ckt, pdest, psource); cplx sij = ckt->CKTSmat->d[pdest][psource]; data[nPlot].real = sij.re; data[nPlot].imag = sij.im; + nPlot++; } } + + // Put Y data + for (unsigned int pdest = 0; pdest < ckt->CKTportCount; pdest++) + { + for (unsigned int psource = 0; psource < ckt->CKTportCount; psource++) + { + //unsigned int nPlot = ckt->CKTmaxEqNum - 1 + CKTmatrixIndex(ckt, pdest, psource); + cplx yij = ckt->CKTYmat->d[pdest][psource]; + data[nPlot].real = yij.re; + data[nPlot].imag = yij.im; + nPlot++; + } + } + + // Put Z data + for (unsigned int pdest = 0; pdest < ckt->CKTportCount; pdest++) + { + for (unsigned int psource = 0; psource < ckt->CKTportCount; psource++) + { + //unsigned int nPlot = ckt->CKTmaxEqNum - 1 + CKTmatrixIndex(ckt, pdest, psource); + cplx zij = ckt->CKTZmat->d[pdest][psource]; + data[nPlot].real = zij.re; + data[nPlot].imag = zij.im; + nPlot++; + } + } + + if (doNoise) + { + // Put Cy data + for (unsigned int pdest = 0; pdest < ckt->CKTportCount; pdest++) + { + for (unsigned int psource = 0; psource < ckt->CKTportCount; psource++) + { + //unsigned int nPlot = ckt->CKTmaxEqNum - 1 + CKTmatrixIndex(ckt, pdest, psource); + cplx CYij = ckt->CKTNoiseCYmat->d[pdest][psource]; + data[nPlot].real = CYij.re; + data[nPlot].imag = CYij.im; + nPlot++; + } + } + + if (ckt->CKTportCount == 2) + { + // If we have two ports, put also NF, Sopt, NFmin, Rn + data[nPlot].real = NF; + data[nPlot].imag = 0.0; + nPlot++; + + data[nPlot].real = Sopt.re; + data[nPlot].imag = Sopt.im; + nPlot++; + + data[nPlot].real = Fmin; + data[nPlot].imag = 0.0; + nPlot++; + + data[nPlot].real = Rn; + data[nPlot].imag = 0.0; + nPlot++; + } + + + } } + + SPfrontEnd->OUTpData(plot, &freqData, &valueData); @@ -108,4 +214,3 @@ CKTspDump(CKTcircuit *ckt, double freq, runDesc *plot) FREE(data); return(OK); } -#endif diff --git a/src/spicelib/analysis/cktspnoise.c b/src/spicelib/analysis/cktspnoise.c new file mode 100644 index 000000000..546267873 --- /dev/null +++ b/src/spicelib/analysis/cktspnoise.c @@ -0,0 +1,154 @@ +/********** +Copyright 1990 Regents of the University of California. All rights reserved. +Author: 1987 Gary W. Ng +**********/ + +/* + * CKTnoise (ckt, mode, operation, data) + * + * This routine is responsible for naming and evaluating all of the + * noise sources in the circuit. It uses a series of subroutines to + * name and evaluate the sources associated with each model, and then + * it evaluates the noise for the entire circuit. + */ + +#include "ngspice/ngspice.h" +#include "ngspice/cktdefs.h" +#include "ngspice/devdefs.h" +#include "ngspice/iferrmsg.h" +#include "ngspice/noisedef.h" +#include "ngspice/sperror.h" + +#ifdef RFSPICE + +// Derived from CKTnoise + + +int +CKTSPnoise (CKTcircuit *ckt, int mode, int operation, Ndata *data) +{ + NOISEAN *job = (NOISEAN *) ckt->CKTcurJob; + + double outNdens; + int i; + IFvalue outData; /* output variable (points to list of outputs)*/ + IFvalue refVal; /* reference variable (always 0)*/ + 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] ) { + error = DEVices[i]->DEVnoise (mode, operation, ckt->CKThead[i], + ckt,data, &outNdens); + if (error) return (error); + } + } + + switch (operation) { + + case N_OPEN: + + /* take care of the noise for the circuit as a whole */ + + switch (mode) { + + case N_DENS: + + data->namelist = TREALLOC(IFuid, data->namelist, data->numPlots + 1); + + SPfrontEnd->IFnewUid (ckt, &(data->namelist[data->numPlots++]), + NULL, "onoise_spectrum", UID_OTHER, NULL); + + data->namelist = TREALLOC(IFuid, data->namelist, data->numPlots + 1); + + SPfrontEnd->IFnewUid (ckt, &(data->namelist[data->numPlots++]), + NULL, "inoise_spectrum", UID_OTHER, NULL); + + /* we've added two more plots */ + + data->outpVector = + TMALLOC(double, data->numPlots); + data->squared_value = + data->squared ? NULL : TMALLOC(char, data->numPlots); + break; + + case INT_NOIZ: + + data->namelist = TREALLOC(IFuid, data->namelist, data->numPlots + 1); + SPfrontEnd->IFnewUid (ckt, &(data->namelist[data->numPlots++]), + NULL, "onoise_total", UID_OTHER, NULL); + + data->namelist = TREALLOC(IFuid, data->namelist, data->numPlots + 1); + SPfrontEnd->IFnewUid (ckt, &(data->namelist[data->numPlots++]), + NULL, "inoise_total", UID_OTHER, NULL); + /* we've added two more plots */ + + data->outpVector = + TMALLOC(double, data->numPlots); + data->squared_value = + data->squared ? NULL : TMALLOC(char, data->numPlots); + break; + + default: + return (E_INTERN); + } + + break; + + case N_CALC: + + switch (mode) { + + case N_DENS: + if ((job->NStpsSm == 0) + || data->prtSummary) + { + data->outpVector[data->outNumber++] = outNdens; + data->outpVector[data->outNumber++] = + (outNdens * data->GainSqInv); + + refVal.rValue = data->freq; /* the reference is the freq */ + if (!data->squared) + for (i = 0; i < data->outNumber; i++) + if (data->squared_value[i]) + data->outpVector[i] = sqrt(data->outpVector[i]); + outData.v.numValue = data->outNumber; /* vector number */ + outData.v.vec.rVec = data->outpVector; /* vector of outputs */ + SPfrontEnd->OUTpData (data->NplotPtr, &refVal, &outData); + } + break; + + case INT_NOIZ: + data->outpVector[data->outNumber++] = data->outNoiz; + data->outpVector[data->outNumber++] = data->inNoise; + if (!data->squared) + for (i = 0; i < data->outNumber; i++) + if (data->squared_value[i]) + data->outpVector[i] = sqrt(data->outpVector[i]); + outData.v.vec.rVec = data->outpVector; /* vector of outputs */ + outData.v.numValue = data->outNumber; /* vector number */ + SPfrontEnd->OUTpData (data->NplotPtr, &refVal, &outData); + break; + + default: + return (E_INTERN); + } + break; + + case N_CLOSE: + SPfrontEnd->OUTendPlot (data->NplotPtr); + FREE(data->namelist); + FREE(data->outpVector); + FREE(data->squared_value); + break; + + default: + return (E_INTERN); + } + return (OK); +} + +#endif \ No newline at end of file diff --git a/src/spicelib/analysis/nevalsrc.c b/src/spicelib/analysis/nevalsrc.c index d7edc8895..b2acf900f 100644 --- a/src/spicelib/analysis/nevalsrc.c +++ b/src/spicelib/analysis/nevalsrc.c @@ -15,15 +15,81 @@ Author: 1987 Gary W. Ng */ +/* +Modified by Alessio Cacciatori: S-Params noise is calculated as noise current +correlation matrix. Per each noise source, the noise voltage at RF ports is converted +as an input noise current source by using (already availalble) Y matrix. +Outside RFSPICE declaration, code is legacy NGSPICE code. +*/ + #include "ngspice/ngspice.h" #include "ngspice/cktdefs.h" #include "ngspice/const.h" #include "ngspice/noisedef.h" +#ifdef RFSPICE +#include "../maths/dense/denseinlines.h" + +extern CMat* eyem; +extern CMat* zref; +extern CMat* gn; +extern CMat* gninv; +extern CMat* vNoise; +extern CMat* iNoise; + +#endif void NevalSrc (double *noise, double *lnNoise, CKTcircuit *ckt, int type, int node1, int node2, double param) { +#ifdef RFSPICE + if (ckt->CKTcurrentAnalysis & DOING_SP) + { + double inoise = 0.0; + + switch (type) { + + case SHOTNOISE: + inoise = 2 * CHARGE * fabs(param); /* param is the dc current in a semiconductor */ + break; + + case THERMNOISE: + inoise = 4 * CONSTboltz * ckt->CKTtemp * param; /* param is the conductance of a resistor */ + break; + + case N_GAIN: + inoise = 0.0; + break; + + } + inoise = sqrt(inoise); + // Calculate input equivalent noise current source (we have port impedance attached) + for (unsigned int s = 0; s < ckt->CKTportCount; s++) + vNoise->d[0][s] = cmultdo(csubco(ckt->CKTadjointRHS->d[s][node1], ckt->CKTadjointRHS->d[s][node2]), inoise); + + for (unsigned int d = 0; d < ckt->CKTportCount; d++) + { + cplx in; + double yport = 1.0 / zref->d[d][d].re; + + in.re = vNoise->d[0][d].re * yport; + in.im = vNoise->d[0][d].im * yport; + + for (unsigned int s = 0; s < ckt->CKTportCount; s++) + caddc(&in, in, cmultco(ckt->CKTYmat->d[d][s], vNoise->d[0][s])); + + iNoise->d[0][d] = in; + } + + + for (unsigned int d = 0; d < ckt->CKTportCount; d++) + for (unsigned int s = 0; s < ckt->CKTportCount; s++) + ckt->CKTNoiseCYmat->d[d][s] = caddco(ckt->CKTNoiseCYmat->d[d][s], cmultco(iNoise->d[0][d], conju(iNoise->d[0][s]))); + + return; + } + +#endif double realVal; double imagVal; double gain; @@ -87,6 +153,69 @@ double phi21) /* Phase of signal 2 relative to signal 1 */ double realOut, imagOut, param_gain; double T0, T1, T2, T3; +#ifdef RFSPICE + if (ckt->CKTcurrentAnalysis & DOING_SP) + { + + double knoise = 0.0; + + T0 = sqrt(param1); + T1 = sqrt(param2); + cplx cfact; + cfact.re = cos(phi21); + cfact.im = sin(phi21); + + + switch (type) { + + case SHOTNOISE: + knoise = 2 * CHARGE; /* param is the dc current in a semiconductor */ + break; + + case THERMNOISE: + knoise = 4 * CONSTboltz * ckt->CKTtemp; /* param is the conductance of a resistor */ + break; + + case N_GAIN: + knoise = 0.0; + break; + + } + + knoise = sqrt(knoise); + // Calculate input equivalent noise current source (we have port impedance attached) + for (unsigned int s = 0; s < ckt->CKTportCount; s++) + { + cplx vNoiseA = cmultdo(csubco(ckt->CKTadjointRHS->d[s][node1], ckt->CKTadjointRHS->d[s][node2]), knoise * sqrt(param1) ); + cplx vNoiseB = cmultco(cmultdo(csubco(ckt->CKTadjointRHS->d[s][node3], ckt->CKTadjointRHS->d[s][node4]), knoise * sqrt(param1)), cfact); + + vNoise->d[0][s] = caddco(vNoiseA, vNoiseB); + } + + for (unsigned int d = 0; d < ckt->CKTportCount; d++) + { + double yport = 1.0 / zref->d[d][d].re; + + cplx in; + in.re = vNoise->d[0][d].re * yport; + in.im = vNoise->d[0][d].im * yport; + + for (unsigned int s = 0; s < ckt->CKTportCount; s++) + caddc(&in, in, cmultco(ckt->CKTYmat->d[d][s], vNoise->d[0][s])); + + iNoise->d[0][d] = in; + } + + + for (unsigned int d = 0; d < ckt->CKTportCount; d++) + for (unsigned int s = 0; s < ckt->CKTportCount; s++) + ckt->CKTNoiseCYmat->d[d][s] = caddco(ckt->CKTNoiseCYmat->d[d][s], cmultco(iNoise->d[0][d], conju(iNoise->d[0][s]))); + + return; + } +#endif + + realVal1 = ckt->CKTrhs [node1] - ckt->CKTrhs [node2]; imagVal1 = ckt->CKTirhs [node1] - ckt->CKTirhs [node2]; realVal2 = ckt->CKTrhs [node3] - ckt->CKTrhs [node4]; @@ -130,10 +259,59 @@ void NevalSrcInstanceTemp (double *noise, double *lnNoise, CKTcircuit *ckt, int type, int node1, int node2, double param, double param2) { + + +#ifdef RFSPICE + if (ckt->CKTcurrentAnalysis & DOING_SP) + { + double inoise = 0.0; + + switch (type) { + + case SHOTNOISE: + inoise = 2 * CHARGE * fabs(param); /* param is the dc current in a semiconductor */ + break; + + case THERMNOISE: + inoise = 4.0 *CONSTboltz* (ckt->CKTtemp + param2)* param; /* param is the conductance of a resistor */ + break; + + case N_GAIN: + inoise = 0.0; + return; + + } + + inoise = sqrt(inoise); + // Calculate input equivalent noise current source (we have port impedance attached) + for (unsigned int s = 0; s < ckt->CKTportCount; s++) + vNoise->d[0][s] = cmultdo(csubco(ckt->CKTadjointRHS->d[s][node1], ckt->CKTadjointRHS->d[s][node2]),inoise); + + for (unsigned int d = 0; d < ckt->CKTportCount; d++) + { + cplx in; + double yport = 1.0 / zref->d[d][d].re; + + in.re = vNoise->d[0][d].re * yport; + in.im = vNoise->d[0][d].im * yport; + + for (unsigned int s = 0; s < ckt->CKTportCount; s++) + caddc(&in, in, cmultco(ckt->CKTYmat->d[d][s], vNoise->d[0][s])); + + iNoise->d[0][d] = in; + } + + + for (unsigned int d = 0; d < ckt->CKTportCount; d++) + for (unsigned int s = 0; s < ckt->CKTportCount; s++) + ckt->CKTNoiseCYmat->d[d][s] = caddco(ckt->CKTNoiseCYmat->d[d][s], cmultco(iNoise->d[0][d], conju(iNoise->d[0][s]))); + + return; + } +#endif double realVal; double imagVal; double gain; - realVal = ckt->CKTrhs [node1] - ckt->CKTrhs [node2]; imagVal = ckt->CKTirhs [node1] - ckt->CKTirhs [node2]; gain = (realVal*realVal) + (imagVal*imagVal); diff --git a/src/spicelib/analysis/span.c b/src/spicelib/analysis/span.c index ac9672333..a53c871dd 100644 --- a/src/spicelib/analysis/span.c +++ b/src/spicelib/analysis/span.c @@ -5,7 +5,7 @@ */ #include "ngspice/ngspice.h" #include "ngspice/cktdefs.h" -#include "ngspice/spdefs.h" +#include "ngspice/spardefs.h" #include "ngspice/devdefs.h" #include "ngspice/sperror.h" @@ -19,12 +19,16 @@ /* gtri - end - wbk */ #endif - +#define SQR(x) ((x) * (x)) #ifdef RFSPICE +#include "vsrc/vsrcdefs.h" #include "vsrc/vsrcext.h" #include "../maths/dense/dense.h" +#include "../maths/dense/denseinlines.h" + + #define INIT_STATS() \ do { \ @@ -45,7 +49,23 @@ do { \ ckt->CKTstat->STATacSyncTime += ckt->CKTstat->STATsyncTime - startkTime; \ } while(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) @@ -54,16 +74,21 @@ CKTspnoise(CKTcircuit * ckt, int mode, int operation, Ndata * data) double outNdens; int i; +#ifdef LEGACY IFvalue outData; /* output variable (points to list of outputs)*/ IFvalue refVal; /* reference variable (always 0)*/ +#endif int error; - outNdens = 0.0; + job->NStpsSm = 1; /* 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) return (error); @@ -73,56 +98,51 @@ CKTspnoise(CKTcircuit * ckt, int mode, int operation, Ndata * data) switch (operation) { case N_OPEN: - - /* take care of the noise for the circuit as a whole */ - - switch (mode) { - - case N_DENS: - - data->namelist = TREALLOC(IFuid, data->namelist, data->numPlots + 1); - - SPfrontEnd->IFnewUid(ckt, &(data->namelist[data->numPlots++]), - NULL, "onoise_spectrum", UID_OTHER, NULL); - - data->namelist = TREALLOC(IFuid, data->namelist, data->numPlots + 1); - - SPfrontEnd->IFnewUid(ckt, &(data->namelist[data->numPlots++]), - NULL, "inoise_spectrum", UID_OTHER, NULL); - - /* we've added two more plots */ - - data->outpVector = - TMALLOC(double, data->numPlots); - data->squared_value = - data->squared ? NULL : TMALLOC(char, data->numPlots); - break; - - case INT_NOIZ: - - data->namelist = TREALLOC(IFuid, data->namelist, data->numPlots + 1); - SPfrontEnd->IFnewUid(ckt, &(data->namelist[data->numPlots++]), - NULL, "onoise_total", UID_OTHER, NULL); - - data->namelist = TREALLOC(IFuid, data->namelist, data->numPlots + 1); - SPfrontEnd->IFnewUid(ckt, &(data->namelist[data->numPlots++]), - NULL, "inoise_total", UID_OTHER, NULL); - /* we've added two more plots */ - - data->outpVector = - TMALLOC(double, data->numPlots); - data->squared_value = - data->squared ? NULL : TMALLOC(char, data->numPlots); - break; - - default: - return (E_INTERN); - } - + // 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); + + + if (ckt->CKTportCount == 2) + { + + double Y21mod = cmodsqr(ckt->CKTYmat->d[1][0]); + Rn = (tempCy->d[1][1].re / Y21mod) ; + cplx Ycor = csubco(ckt->CKTYmat->d[0][0], + cmultco( + cdivco(tempCy->d[0][1], tempCy->d[1][1]), + tempCy->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 / Ysopt.re) * SQR(Ysoptmod); + Fmin = 10.0 * log10(Fmin); + NF = 10.0 * log10(NF); + } + + freecmat(tempCy); + } + + +#ifdef LEGACY switch (mode) { case N_DENS: @@ -159,6 +179,7 @@ CKTspnoise(CKTcircuit * ckt, int mode, int operation, Ndata * data) default: return (E_INTERN); } +#endif break; case N_CLOSE: @@ -166,6 +187,10 @@ CKTspnoise(CKTcircuit * ckt, int mode, int operation, Ndata * data) FREE(data->namelist); FREE(data->outpVector); FREE(data->squared_value); + freecmat(ckt->CKTNoiseCYmat); + freecmat(ckt->CKTadjointRHS); + ckt->CKTNoiseCYmat = NULL; + ckt->CKTadjointRHS = NULL; break; default: @@ -175,8 +200,8 @@ CKTspnoise(CKTcircuit * ckt, int mode, int operation, Ndata * data) } -void -NInspIter(CKTcircuit * ckt, int posDrive, int negDrive) +int +NInspIter(CKTcircuit * ckt, VSRCinstance* port) { int i; @@ -187,15 +212,130 @@ NInspIter(CKTcircuit * ckt, int posDrive, int negDrive) ckt->CKTirhs[i] = 0.0; } - ckt->CKTrhs[posDrive] = 1.0; /* apply unit current excitation */ - ckt->CKTrhs[negDrive] = -1.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 initSPmatrix(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 (eyem != NULL) freecmat(eyem); + if (zref != NULL) freecmat(zref); + if (gn != NULL) freecmat(gn); + if (gninv != NULL) freecmat(gninv); + + + 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); + + 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; + + } + return (OK); +} + +void deleteSPmatrix(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 (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; +} int SPan(CKTcircuit *ckt, int restart) @@ -221,8 +361,14 @@ SPan(CKTcircuit *ckt, int restart) double* rhswoPorts = NULL; double* irhswoPorts = NULL; - int* portPosNodes = NULL; - int* portNegNodes = NULL; + + /* 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) @@ -233,21 +379,6 @@ SPan(CKTcircuit *ckt, int restart) - if (ckt->CKTAmat != NULL) freecmat(ckt->CKTAmat); - if (ckt->CKTBmat != NULL) freecmat(ckt->CKTBmat); - if (ckt->CKTSmat != NULL) freecmat(ckt->CKTSmat); - - 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); - #ifdef XSPICE /* gtri - add - wbk - 12/19/90 - Add IPC stuff and anal_init and anal_type */ @@ -276,7 +407,7 @@ SPan(CKTcircuit *ckt, int restart) return E_PARMVAL; } job->SPfreqDelta = - exp(log(10.0)/job->SPnumberSteps); + exp(log(10.0) / job->SPnumberSteps); break; case OCTAVE: if (job->SPstartFreq <= 0) { @@ -284,108 +415,122 @@ SPan(CKTcircuit *ckt, int restart) return E_PARMVAL; } job->SPfreqDelta = - exp(log(2.0)/job->SPnumberSteps); + exp(log(2.0) / job->SPnumberSteps); break; case LINEAR: - if (job->SPnumberSteps-1 > 1) + if (job->SPnumberSteps - 1 > 1) job->SPfreqDelta = - (job->SPstopFreq - - job->SPstartFreq) / - (job->SPnumberSteps - 1); + (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 */ + /* 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); - } -#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"); - + + if (job->SPdoNoise) + { + data->lstFreq = job->SPstartFreq - 1; + data->delFreq = 1.0; + } + #ifdef XSPICE -/* gtri - add - wbk - 12/19/90 - Add IPC stuff */ + /* gtri - begin - wbk - Call EVTop if event-driven instances exist */ - /* Send the operating point results for Mspice compatibility */ - if(g_ipc.enabled) - { - /* Call CKTnames to get names of nodes/branches used by - BeginPlot */ - /* Probably should free nameList after this block since - called again... */ - error = CKTnames(ckt,&numNames,&nameList); - if(error) return(error); + 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); - /* We have to do a beginPlot here since the data to return is - * different for the DCOP than it is for the AC analysis. - * Moreover the begin plot has not even been done yet at this - * point... - */ - SPfrontEnd->OUTpBeginPlot (ckt, ckt->CKTcurJob, - ckt->CKTcurJob->JOBname, - NULL, IF_REAL, - numNames, nameList, IF_REAL, - &spPlot); - txfree(nameList); + 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"); - ipc_send_dcop_prefix(); - CKTdump(ckt, 0.0, spPlot); - ipc_send_dcop_suffix(); +#ifdef XSPICE + /* gtri - add - wbk - 12/19/90 - Add IPC stuff */ - SPfrontEnd->OUTendPlot (spPlot); - } -/* gtri - end - wbk */ + /* Send the operating point results for Mspice compatibility */ + if (g_ipc.enabled) + { + /* Call CKTnames to get names of nodes/branches used by + BeginPlot */ + /* Probably should free nameList after this block since + called again... */ + error = CKTnames(ckt, &numNames, &nameList); + if (error) return(error); + + /* We have to do a beginPlot here since the data to return is + * different for the DCOP than it is for the AC analysis. + * Moreover the begin plot has not even been done yet at this + * point... + */ + SPfrontEnd->OUTpBeginPlot(ckt, ckt->CKTcurJob, + ckt->CKTcurJob->JOBname, + NULL, IF_REAL, + numNames, nameList, IF_REAL, + &spPlot); + txfree(nameList); + + ipc_send_dcop_prefix(); + CKTdump(ckt, 0.0, spPlot); + ipc_send_dcop_suffix(); + + SPfrontEnd->OUTendPlot(spPlot); + } + /* gtri - end - wbk */ #endif ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG; error = CKTload(ckt); - if(error) return(error); + if (error) return(error); - error = CKTnames(ckt,&numNames,&nameList); - 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; + } + + unsigned int extraSPdataLength = 3 * ckt->CKTportCount * ckt->CKTportCount; + if (job->SPdoNoise) + { + extraSPdataLength += ckt->CKTportCount * ckt->CKTportCount; // Add Cy + if (ckt->CKTportCount == 2) + extraSPdataLength += 4; + } - 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; - } - - unsigned int extraSPdataLength = ckt->CKTportCount * ckt->CKTportCount; nameList = (IFuid*)TREALLOC(IFuid, nameList, numNames + extraSPdataLength); @@ -399,12 +544,57 @@ SPan(CKTcircuit *ckt, int restart) SPfrontEnd->IFnewUid(ckt, &(nameList[numNames++]), NULL, tmpBuf, 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); + // Create UIDs + for (unsigned int dest = 1; dest <= ckt->CKTportCount; dest++) + for (unsigned 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 (unsigned int dest = 1; dest <= ckt->CKTportCount; dest++) + for (unsigned 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 (unsigned int dest = 1; dest <= ckt->CKTportCount; dest++) + for (unsigned 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); @@ -455,7 +645,20 @@ SPan(CKTcircuit *ckt, int restart) INIT_STATS(); - ckt->CKTcurrentAnalysis = DOING_AC; + ckt->CKTcurrentAnalysis = DOING_AC | DOING_SP; + + if (initSPmatrix(ckt, job->SPdoNoise)) + return (E_NOMEM); + + // Create Noise UID, if needed + + if (job->SPdoNoise) + { + + data->numPlots = 0; /* we don't have any plots yet */ + error = CKTspnoise(ckt, N_DENS, N_OPEN, data); + if (error) return(error); + } ckt->CKTactivePort = 0; /* main loop through all scheduled frequencies */ @@ -499,8 +702,33 @@ SPan(CKTcircuit *ckt, int restart) } ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG; error = CKTload(ckt); - if (error) return(error); + if (error) { + tfree(data); return(error); + } } + + // 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) @@ -522,16 +750,14 @@ SPan(CKTcircuit *ckt, int restart) return (E_NOMOD); ckt->CKTVSRCid = vsrcRoot; + + // Now that we have found the model, we may init the Zref and Gn ports + VSRCspinit(ckt->CKThead[vsrcRoot], ckt, zref, gn, gninv); } else vsrcRoot = ckt->CKTVSRCid; - if (rhswoPorts == NULL) - rhswoPorts = (double*)TREALLOC(double, rhswoPorts, ckt->CKTmaxEqNum); - if (irhswoPorts == NULL) - irhswoPorts = (double*)TREALLOC(double, irhswoPorts, ckt->CKTmaxEqNum); - ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODESP; // Pre-load everything but RF Ports (these will be updated in the next cycle). error = NIspPreload(ckt); if (error) return (error); @@ -556,6 +782,8 @@ SPan(CKTcircuit *ckt, int restart) { tfree(rhswoPorts); tfree(irhswoPorts); + tfree(data); + deleteSPmatrix(ckt); return(error); } @@ -563,6 +791,8 @@ SPan(CKTcircuit *ckt, int restart) if (error) { tfree(rhswoPorts); tfree(irhswoPorts); + tfree(data); + deleteSPmatrix(ckt); UPDATE_STATS(DOING_AC); return(error); } @@ -600,83 +830,78 @@ SPan(CKTcircuit *ckt, int restart) // Now we can calculate the full S-Matrix CKTspCalcSMatrix(ckt); -#ifdef XSPICE - /* gtri - modify - wbk - 12/19/90 - Send IPC stuff */ - - if (g_ipc.enabled) - ipc_send_data_prefix(freq); - - error = CKTspDump(ckt, freq, spPlot); - - if (g_ipc.enabled) - ipc_send_data_suffix(); - - /* gtri - modify - wbk - 12/19/90 - Send IPC stuff */ -#else - error = CKTspDump(ckt, freq, spPlot); -#endif - if (error) { - UPDATE_STATS(DOING_AC); - tfree(rhswoPorts); - tfree(irhswoPorts); - return(error); - } - - /* * Now go with noise cycle, if required */ - -#ifdef NOISE_AVAILABLE - - // To be completed if (job->SPdoNoise) { - if (portPosNodes == NULL) - { - portPosNodes = TMALLOC(int, ckt->CKTportCount); - portNegNodes = TMALLOC(int, ckt->CKTportCount); - VSRCgetActivePortNodes(ckt->CKThead[vsrcRoot], ckt, portPosNodes, portNegNodes); - } - static Ndata* data; + data->delFreq = freq - data->lstFreq; + data->freq = freq; - double realVal; - double imagVal; - int error; - int posOutNode; - int negOutNode; - //Keep a backup copy - memcpy(rhswoPorts, ckt->CKTrhs, ckt->CKTmaxEqNum * sizeof(double)); - memcpy(rhswoPorts, ckt->CKTirhs, ckt->CKTmaxEqNum * sizeof(double)); + + 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) */ - // Copy the backup RHS into CKT's RHS - memcpy(ckt->CKTrhs, rhswoPorts, ckt->CKTmaxEqNum * sizeof(double)); - memcpy(ckt->CKTirhs, irhswoPorts, ckt->CKTmaxEqNum * sizeof(double)); ckt->CKTactivePort = activePort+1; - posOutNode = portPosNodes[activePort]; - negOutNode = portNegNodes[activePort]; - NInspIter(ckt, posOutNode, negOutNode); /* solve the adjoint system */ - - /* now we use the adjoint system to calculate the noise - * contributions of each generator in the circuit - */ - - error = CKTspnoise(ckt, N_DENS, N_CALC, data); - if (error) + 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++) { - tfree(portPosNodes); tfree(portNegNodes); - return(error); + 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); + if (error) + { + tfree(data); + tfree(rhswoPorts); + tfree(irhswoPorts); + deleteSPmatrix(ckt); + return(error); + } + data->lstFreq = freq; + } + + +#ifdef XSPICE + /* gtri - modify - wbk - 12/19/90 - Send IPC stuff */ + + if (g_ipc.enabled) + ipc_send_data_prefix(freq); + + error = CKTspDump(ckt, freq, spPlot, job->SPdoNoise); + + if (g_ipc.enabled) + ipc_send_data_suffix(); + + /* gtri - modify - wbk - 12/19/90 - Send IPC stuff */ +#else + error = CKTspDump(ckt, freq, acPlot, job->SPdoNoise)); +#endif + if (error) { + UPDATE_STATS(DOING_AC); + tfree(rhswoPorts); + tfree(irhswoPorts); + tfree(data); + deleteSPmatrix(ckt); + return(error); } -#endif /* increment frequency */ switch (job->SPstepType) { @@ -717,7 +942,8 @@ SPan(CKTcircuit *ckt, int restart) default: tfree(rhswoPorts); tfree(irhswoPorts); - tfree(portPosNodes); tfree(portNegNodes); + tfree(data); + deleteSPmatrix(ckt); return(E_INTERN); } @@ -728,7 +954,8 @@ endsweep: UPDATE_STATS(0); tfree(rhswoPorts); tfree(irhswoPorts); - tfree(portPosNodes); tfree(portNegNodes); + deleteSPmatrix(ckt); + tfree(data); return(0); } diff --git a/src/spicelib/analysis/spaskq.c b/src/spicelib/analysis/spaskq.c index b299b2338..b858b955a 100644 --- a/src/spicelib/analysis/spaskq.c +++ b/src/spicelib/analysis/spaskq.c @@ -8,7 +8,7 @@ Author: 1985 Thomas L. Quarles #include "ngspice/ngspice.h" #include "ngspice/ifsim.h" #include "ngspice/iferrmsg.h" -#include "ngspice/spdefs.h" +#include "ngspice/spardefs.h" #include "ngspice/cktdefs.h" #ifdef RFSPICE diff --git a/src/spicelib/analysis/spsetp.c b/src/spicelib/analysis/spsetp.c index 7a456717c..47a9cc0a4 100644 --- a/src/spicelib/analysis/spsetp.c +++ b/src/spicelib/analysis/spsetp.c @@ -6,7 +6,7 @@ Author: 1985 Thomas L. Quarles #include "ngspice/ngspice.h" #include "ngspice/ifsim.h" #include "ngspice/iferrmsg.h" -#include "ngspice/spdefs.h" +#include "ngspice/spardefs.h" #include "ngspice/cktdefs.h" #include "analysis.h" diff --git a/src/spicelib/devices/cktinit.c b/src/spicelib/devices/cktinit.c index cc57200ad..7fde41af6 100644 --- a/src/spicelib/devices/cktinit.c +++ b/src/spicelib/devices/cktinit.c @@ -141,6 +141,9 @@ CKTinit(CKTcircuit **ckt) /* new circuit to create */ sckt->CKTVSRCid = -1; sckt->CKTrfPorts = NULL; sckt->CKTSmat = sckt->CKTAmat = sckt->CKTBmat = sckt->CKTYmat = sckt->CKTZmat = NULL; + sckt->CKTNoiseCYmat = NULL; + sckt->CKTadjointRHS = NULL; + sckt->CKTnoiseSourceCount = 0; #endif return OK; diff --git a/src/spicelib/devices/vsrc/vsrcacld.c b/src/spicelib/devices/vsrc/vsrcacld.c index 03f4b785c..f39def77d 100644 --- a/src/spicelib/devices/vsrc/vsrcacld.c +++ b/src/spicelib/devices/vsrc/vsrcacld.c @@ -11,6 +11,32 @@ Author: 1985 Thomas L. Quarles #ifdef RFSPICE + +int VSRCspinit(GENmodel* inModel, CKTcircuit* ckt, CMat* zref, CMat* gn, CMat* gninv) +{ + if (!(ckt->CKTmode & MODESP) && !(ckt->CKTcurrentAnalysis & DOING_SP)) + return (OK); + VSRCmodel* model = (VSRCmodel*)inModel; + VSRCinstance* here; + + for (; model != NULL; model = VSRCnextModel(model)) { + + /* loop through all the instances of the model */ + for (here = VSRCinstances(model); here != NULL; + here = VSRCnextInstance(here)) { + + if (here->VSRCisPort) + { + int i = here->VSRCportNum - 1; + zref->d[i][i].re = here->VSRCportZ0; + gn->d[i][i].re = 2.0 * here->VSRCki; + gninv->d[i][i].re = 1.0 / gn->d[i][i].re; + } + } + } + return (OK); +} + int VSRCspupdate(GENmodel* inModel, CKTcircuit* ckt) { if (!(ckt->CKTmode & MODESP)) @@ -36,12 +62,14 @@ int VSRCspupdate(GENmodel* inModel, CKTcircuit* ckt) } -int VSRCgetActivePortNodes(GENmodel* inModel, CKTcircuit* ckt, int* posNodes, int* negNodes) + +int VSRCgetActivePorts(GENmodel* inModel, CKTcircuit* ckt, VSRCinstance** ports) { if (!(ckt->CKTmode & MODESP)) return (OK); for (unsigned int n = 0; n < ckt->CKTportCount; n++) - posNodes[n] = negNodes[n] = 0; + ports[n] = NULL; + VSRCmodel* model = (VSRCmodel*)inModel; VSRCinstance* here; @@ -54,9 +82,7 @@ int VSRCgetActivePortNodes(GENmodel* inModel, CKTcircuit* ckt, int* posNodes, in if (here->VSRCisPort) { int id = here->VSRCportNum - 1; - posNodes[id] = here->VSRCposNode; - negNodes[id] = here->VSRCnegNode; - + ports[id] = here; } } } @@ -78,9 +104,9 @@ VSRCacLoad(GENmodel* inModel, CKTcircuit* ckt) double acReal, acImag; + #ifdef RFSPICE - double g0; - g0 = 0; + double g0 = 0; acReal = 0.0; acImag = 0.0; /* diff --git a/src/spicelib/devices/vsrc/vsrcext.h b/src/spicelib/devices/vsrc/vsrcext.h index 7761b867a..0a0ace100 100644 --- a/src/spicelib/devices/vsrc/vsrcext.h +++ b/src/spicelib/devices/vsrc/vsrcext.h @@ -18,5 +18,8 @@ extern int VSRCpzSetup(SMPmatrix*,GENmodel*,CKTcircuit*,int*); extern int VSRCtemp(GENmodel*,CKTcircuit*); #ifdef RFSPICE extern int VSRCspupdate(GENmodel*, CKTcircuit*); -extern int VSRCgetActivePortNodes(GENmodel* inModel, CKTcircuit* ckt, int* posNodes, int* negNodes); +#include "vsrcdefs.h" + +extern int VSRCgetActivePorts(GENmodel* inModel, CKTcircuit* ckt, VSRCinstance** ports); +extern int VSRCspinit(GENmodel* inModel, CKTcircuit* ckt, CMat* zref, CMat* gn, CMat* gninv); #endif diff --git a/src/spicelib/devices/vsrc/vsrcload.c b/src/spicelib/devices/vsrc/vsrcload.c index 9c2bb52d3..96422b5f2 100644 --- a/src/spicelib/devices/vsrc/vsrcload.c +++ b/src/spicelib/devices/vsrc/vsrcload.c @@ -41,7 +41,7 @@ VSRCload(GENmodel *inModel, CKTcircuit *ckt) here=VSRCnextInstance(here)) { #ifndef RFSPICE - *(here->VSRCposIbrPtr) += 1.0; + * (here->VSRCposIbrPtr) += 1.0; *(here->VSRCnegIbrPtr) -= 1.0; *(here->VSRCibrPosPtr) += 1.0; *(here->VSRCibrNegPtr) -= 1.0; @@ -60,7 +60,7 @@ VSRCload(GENmodel *inModel, CKTcircuit *ckt) *(here->VSRCnegNegPtr) += g0; *(here->VSRCposNegPtr) -= g0; *(here->VSRCnegPosPtr) -= g0; - } + } else { *(here->VSRCposIbrPtr) += 1.0; diff --git a/src/spicelib/devices/vsrc/vsrctemp.c b/src/spicelib/devices/vsrc/vsrctemp.c index 2fb776bc4..730397f33 100644 --- a/src/spicelib/devices/vsrc/vsrctemp.c +++ b/src/spicelib/devices/vsrc/vsrctemp.c @@ -88,7 +88,26 @@ VSRCtemp(GENmodel *inModel, CKTcircuit *ckt) ckt->CKTrfPorts = (GENinstance**)TREALLOC(GENinstance*, ckt->CKTrfPorts, ckt->CKTportCount); ckt->CKTrfPorts[ckt->CKTportCount - 1] = (GENinstance*)here; - + // Reorder ports according to their PortNum + unsigned int done = 0; + while (!done) + { + int nMax = ckt->CKTportCount - 1; + done = 1; + for (int n = 0; n < nMax; n++) + { + VSRCinstance* a = (VSRCinstance*)ckt->CKTrfPorts[n]; + VSRCinstance* b = (VSRCinstance*)ckt->CKTrfPorts[n + 1]; + if (a->VSRCportNum > b->VSRCportNum) + { + // Swap a and b. Restart + done = 0; + ckt->CKTrfPorts[n] = (GENinstance*)b; + ckt->CKTrfPorts[n + 1] = (GENinstance*)a; + break; + } + } + } } #endif } diff --git a/src/spicelib/parser/inp2dot.c b/src/spicelib/parser/inp2dot.c index 625372d11..f1ba03623 100644 --- a/src/spicelib/parser/inp2dot.c +++ b/src/spicelib/parser/inp2dot.c @@ -616,7 +616,7 @@ dot_pss(char *line, void *ckt, INPtables *tab, struct card *current, #ifdef RFSPICE -/* S Parameter Analyis */ +/* S-Parameter Analyis */ static int dot_sp(char* line, void* ckt, INPtables* tab, struct card* current, void* task, void* gnode, JOB* foo) @@ -652,7 +652,7 @@ dot_sp(char* line, void* ckt, INPtables* tab, struct card* current, } #ifdef WITH_HB -/* HB */ +/*SP: Steady State Analyis */ static int dot_hb(char* line, void* ckt, INPtables* tab, struct card* current, void* task, void* gnode, JOB* foo) @@ -714,6 +714,7 @@ dot_hb(char* line, void* ckt, INPtables* tab, struct card* current, #endif + static int dot_options(char *line, CKTcircuit *ckt, INPtables *tab, struct card *current, TSKtask *task, CKTnode *gnode, JOB *foo) diff --git a/visualc/vngspice.vcxproj b/visualc/vngspice.vcxproj index af6add978..2721aad2c 100644 --- a/visualc/vngspice.vcxproj +++ b/visualc/vngspice.vcxproj @@ -952,7 +952,7 @@ - + @@ -1067,6 +1067,7 @@ + @@ -1686,6 +1687,7 @@ + @@ -1709,6 +1711,7 @@ +