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

S to Z matrix conversion by formula instead of Y inversion 

low limiting Rn and Cy to prevent division by 0, fix provided by Alessio Cacciatori

there are still problems in Z matrix conversion in specific networks
This commit is contained in:
dwarning 2025-01-03 18:21:46 +01:00
parent caa0a7f4a8
commit ea33459ba9
2 changed files with 54 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
src/spicelib/analysis/cktspdum.c
View File

@ -29,6 +29,9 @@ extern double Rn;
extern cplx Sopt;
extern double Fmin;
#ifdef TRACE
cplx cdet(CMat* M);
#endif
int CKTmatrixIndex(CKTcircuit* ckt, int source, int dest)
{
@ -39,26 +42,42 @@ int CKTspCalcSMatrix(CKTcircuit* ckt)
{
CMat* Ainv = cinverse(ckt->CKTAmat);
if (Ainv == NULL) return (E_NOMEM);
cmultiplydest(ckt->CKTBmat, Ainv, ckt->CKTSmat);
cmultiplydest(ckt->CKTBmat, Ainv, ckt->CKTSmat); // S = B * A-1
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);
// Calculate Z matrix 1. Formula
CMat* SxZ0 = cmultiply(ckt->CKTSmat, zref); // S * Z0
CMat* SxZ0pZ0 = csum(SxZ0, zref); // S * Z0 + Z0
CMat* SxZ0pZ0xGn = cmultiply(SxZ0pZ0, gn); // (S * Z0 + Z0) * Gn
CMat* EmS = cminus(eyem, ckt->CKTSmat); // E - S
CMat* EmSinv = cinverse(EmS); // (E - S)-1
CMat* EmSinvxSxZ0pZ0xGn = cmultiply(EmSinv, SxZ0pZ0xGn); // (E - S)-1 * (S * Z0 + Z0) * Gn
cmultiplydest(gninv, EmSinvxSxZ0pZ0xGn, ckt->CKTZmat); // Z = Gn-1 * (E - S)-1 * (S * Z0 + Z0) * Gn
#ifdef TRACE
cplx de = cdet(ckt->CKTZmat);
printf("spCalcSMatrix: CKTZmat det: %g %g\n", de.re, de.im);
showcmat(ckt->CKTZmat);
#endif
// Calculate Y matrix 1. Formula
CMat* EmSxGn = cmultiply(EmS, gn); // (E - S) * Gn
CMat* SxZ0pZ0inv = cinverse(SxZ0pZ0); // (S * Z0 + Z0)-1
CMat* SxZ0pZ0invxEmSxGn = cmultiply(SxZ0pZ0inv, EmSxGn); // (S * Z0 + Z0)-1 * (E - S) * Gn
cmultiplydest(gninv, SxZ0pZ0invxEmSxGn, ckt->CKTYmat); // Y = Gn-1 * (S * Z0 + Z0)-1 * (E - S) * Gn
#ifdef TRACE
de = cdet(ckt->CKTYmat);
printf("spCalcSMatrix: CKTYmat det: %g %g\n", de.re, de.im);
showcmat(ckt->CKTYmat);
#endif
freecmat(temp);
freecmat(temp2);
freecmat(temp3);
freecmat(temp4);
freecmat(temp5);
freecmat(SxZ0);
freecmat(SxZ0pZ0);
freecmat(SxZ0pZ0xGn);
freecmat(EmS);
freecmat(EmSinv);
freecmat(EmSinvxSxZ0pZ0xGn);
freecmat(EmSxGn);
freecmat(SxZ0pZ0inv);
freecmat(SxZ0pZ0invxEmSxGn);
return (OK);
}
@ -96,7 +115,6 @@ int CKTspCalcPowerWave(CKTcircuit* ckt)
return (OK);
}
int
CKTspDump(CKTcircuit *ckt, double freq, runDesc *plot, int doNoise)
{

21
src/spicelib/analysis/span.c
View File

@ -116,13 +116,23 @@ CKTspnoise(CKTcircuit* ckt, int mode, int operation, Ndata* data, NOISEAN* noise
// 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]),
@ -869,7 +879,12 @@ SPan(CKTcircuit* ckt, int restart)
} //active ports cycle
#ifdef TRACE
printf("SPan: CKTAmat\n");
showcmat(ckt->CKTAmat);
printf("SPan: CKTBmat\n");
showcmat(ckt->CKTBmat);
#endif
// Now we can calculate the full S-Matrix
CKTspCalcSMatrix(ckt);