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

trial with qd node

This commit is contained in:
dwarning 2025-09-17 16:55:01 +02:00
parent 3161e06584
commit b409e7ab27
6 changed files with 82 additions and 856 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
src/spicelib/devices/dio/diodefs.h
View File

@ -40,6 +40,7 @@ typedef struct sDIOinstance {
const int DIOnegNode; /* number of negative node of diode */
const int DIOtempNode; /* number of the temperature node of the diode */
int DIOposPrimeNode; /* number of positive prime node of diode */
int DIOqdNode; /* number of qdif node */
double *DIOposPosPrimePtr; /* pointer to sparse matrix at
* (positive,positive prime) */
@ -65,6 +66,11 @@ typedef struct sDIOinstance {
double *DIOposPrimeTempPtr;
double *DIOnegTempPtr;
/* rev-rec */
double *DIOqdQdPtr;
double *DIOqdPosPrimePtr;
double *DIOqdNegPtr;
double DIOcap; /* stores the diode capacitance */
double *DIOsens; /* stores the perturbed values of geq and ceq in ac
@ -218,21 +224,18 @@ typedef struct sDIOinstance {
#define DIOconduct DIOstate+2
#define DIOcapCharge DIOstate+3
#define DIOcapCurrent DIOstate+4
#define DIOdiffCharge DIOstate+5
#define DIOdiffCurrent DIOstate+6
#define DIOdiffCap DIOstate+7
#define DIOoldCurr DIOstate+8
#define DIOoldCond DIOstate+9
#define DIOdifCharge DIOstate+5
#define DIOdifCurrent DIOstate+6
#define DIOqth DIOstate+10 /* thermal capacitor charge */
#define DIOcqth DIOstate+11 /* thermal capacitor current */
#define DIOqth DIOstate+7 /* thermal capacitor charge */
#define DIOcqth DIOstate+8 /* thermal capacitor current */
#define DIOdeltemp DIOstate+12 /* thermal voltage over rth0 */
#define DIOdIdio_dT DIOstate+13
#define DIOdeltemp DIOstate+9 /* thermal voltage over rth0 */
#define DIOdIdio_dT DIOstate+10
#define DIOnumStates 14
#define DIOnumStates 11
#define DIOsensxp DIOstate+14 /* charge sensitivities and their derivatives.
#define DIOsensxp DIOstate+11 /* charge sensitivities and their derivatives.
* +10 for the derivatives - pointer to the
* beginning of the array */

85
src/spicelib/devices/dio/dioload.c
View File

@ -2,7 +2,6 @@
Copyright 1990 Regents of the University of California. All rights reserved.
Author: 1985 Thomas L. Quarles
Modified: 2000 AlansFixes
Modified by Paolo Nenzi 2003, Dietmar Warning 2012 and Ste Kulov 2021
**********/
#include "ngspice/ngspice.h"
@ -48,8 +47,8 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
int Check=0;
int error;
int SenCond=0; /* sensitivity condition */
double diffcharge, deplcharge, diffcap, deplcap;
double deplcharge, deplcap;
double difcharge, difcap, cdif=0.0, gdif=0.0;
double tt;
double vp;
@ -108,17 +107,23 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
Check=1;
if(ckt->CKTmode & MODEINITSMSIG) {
vd= *(ckt->CKTstate0 + here->DIOvoltage);
difcharge= *(ckt->CKTstate0 + here->DIOqdNode);
} else if (ckt->CKTmode & MODEINITTRAN) {
vd= *(ckt->CKTstate1 + here->DIOvoltage);
difcharge= *(ckt->CKTstate1 + here->DIOqdNode);
} else if ( (ckt->CKTmode & MODEINITJCT) &&
(ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC) ) {
vd=here->DIOinitCond;
difcharge=0;
} else if ( (ckt->CKTmode & MODEINITJCT) && here->DIOoff) {
vd=0;
difcharge=0;
} else if ( ckt->CKTmode & MODEINITJCT) {
vd=here->DIOtVcrit;
difcharge=0;
} else if ( ckt->CKTmode & MODEINITFIX && here->DIOoff) {
vd=0;
difcharge=0;
} else {
#ifndef PREDICTOR
if (ckt->CKTmode & MODEINITPRED) {
@ -129,10 +134,12 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
*(ckt->CKTstate1 + here->DIOcurrent);
*(ckt->CKTstate0 + here->DIOconduct) =
*(ckt->CKTstate1 + here->DIOconduct);
difcharge = DEVpred(ckt,here->DIOdifCharge);
} else {
#endif /* PREDICTOR */
vd = *(ckt->CKTrhsOld+here->DIOposPrimeNode)-
*(ckt->CKTrhsOld + here->DIOnegNode);
difcharge = *(ckt->CKTrhsOld + here->DIOqdNode);
#ifndef PREDICTOR
}
#endif /* PREDICTOR */
@ -179,22 +186,21 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
/*
* compute dc current and derivitives
*/
next1:
if (vd >= -3*vte) { /* forward */
next1: if (vd >= -3*vte) {
evd = exp(vd/vte);
cd = csat*(evd-1) + ckt->CKTgmin*vd;
gd = csat*evd/vte + ckt->CKTgmin;
} else if((!(model->DIObreakdownVoltageGiven)) ||
vd >= -here->DIOtBrkdwnV) { /* reverse */
vd >= -here->DIOtBrkdwnV) {
arg = 3*vte/(vd*CONSTe);
arg = arg * arg * arg;
cd = -csat*(1+arg) + ckt->CKTgmin*vd;
gd = csat*3*arg/vd + ckt->CKTgmin;
} else { /* breakdown */
} else {
evrev = exp(-(here->DIOtBrkdwnV+vd)/vtebrk);
cd = -csat*evrev + ckt->CKTgmin*vd;
@ -220,55 +226,32 @@ next1:
deplcap = czof2*(here->DIOtF3+here->DIOtGradingCoeff*vd/here->DIOtJctPot);
}
*(ckt->CKTstate0 + here->DIOcapCharge) = deplcharge;
if (model->DIOsoftRevRecParamGiven) {
if (ckt->CKTmode & MODEINITTRAN) {
diffcharge = tt * cd;
diffcap = tt * gd;
*(ckt->CKTstate2 + here->DIOdiffCharge) = diffcharge;
*(ckt->CKTstate2 + here->DIOdiffCap) = diffcap;
*(ckt->CKTstate1 + here->DIOdiffCharge) = diffcharge;
*(ckt->CKTstate1 + here->DIOdiffCap) = diffcap;
*(ckt->CKTstate1 + here->DIOoldCurr) = cd;
*(ckt->CKTstate1 + here->DIOoldCond) = gd;
difcharge = tt * cd;
difcap = tt * gd;
}
else {
double dt = ckt->CKTdelta;
//Backward Euler
// Qk+1 = tt (VP Qk + h Ik+1) / (tt VP + h)
//diffcharge = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + dt*cd) / (tt*vp + dt));
//diffcap = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + dt*gd) / (tt*vp + dt));
//Trap
// Qk+1 = (Qk (2 tt VP - h) + h tt (Ik + Ik+1)) / (2 tt VP + h)
diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCurr)) + cd)) / (2.0*tt*vp + dt);
diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCond)) + gd)) / (2.0*tt*vp + dt);
//Gear
// Qk+1 = tt ((VP (4 Qk - Qk-1) + 2 h Ik+1) / (3 tt VP + 2 h))
//diffcharge = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) + 2.0*dt*cd) / (3.0*tt*vp + 2.0*dt));
//diffcap = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) + 2.0*dt*gd) / (3.0*tt*vp + 2.0*dt));
difcharge = *(ckt->CKTstate0 + here->DIOqdNode);
difcap = tt * gd;
}
//printf("time: %e vd: %e\n", ckt->CKTtime, vd);
//printf("%e %e %e %e %e %e %e\n", ckt->CKTtime, ckt->CKTdelta, diffcharge, diffcap, vd, cd, gd);
} else {
diffcharge = tt*cd;
diffcap = tt*gd;
difcharge = tt*cd;
difcap = tt*gd;
}
*(ckt->CKTstate0 + here->DIOdiffCharge) = diffcharge;
*(ckt->CKTstate0 + here->DIOdiffCap) = diffcap;
*(ckt->CKTstate0 + here->DIOcapCharge) = deplcharge + diffcharge;
*(ckt->CKTstate0 + here->DIOoldCurr) = cd;
*(ckt->CKTstate0 + here->DIOoldCond) = gd;
*(ckt->CKTstate0 + here->DIOdifCharge) = difcharge;
//printf("difcharge = %.7e, difcap = %.7e, cd = %.7e\n", difcharge,difcap,cd);
capd = deplcap + diffcap;
capd = deplcap + difcap;
here->DIOcap = capd;
@ -309,14 +292,24 @@ next1:
if (ckt->CKTmode & MODEINITTRAN) {
*(ckt->CKTstate1 + here->DIOcapCharge) =
*(ckt->CKTstate0 + here->DIOcapCharge);
*(ckt->CKTstate1 + here->DIOdifCharge) =
*(ckt->CKTstate0 + here->DIOdifCharge);
}
error = NIintegrate(ckt,&geq,&ceq,capd,here->DIOcapCharge);
error = NIintegrate(ckt,&geq,&ceq,deplcap,here->DIOcapCharge);
if(error) return(error);
gd=gd+geq;
cd=cd+*(ckt->CKTstate0 + here->DIOcapCurrent);
error = NIintegrate(ckt,&geq,&ceq,difcap,here->DIOdifCharge);
if(error) return(error);
gd=gd+geq;
cd=cd+*(ckt->CKTstate0 + here->DIOdifCurrent);
gdif=geq;
cdif=*(ckt->CKTstate0 + here->DIOdifCurrent);
if (ckt->CKTmode & MODEINITTRAN) {
*(ckt->CKTstate1 + here->DIOcapCurrent) =
*(ckt->CKTstate0 + here->DIOcapCurrent);
*(ckt->CKTstate1 + here->DIOdifCurrent) =
*(ckt->CKTstate0 + here->DIOdifCurrent);
}
}
}
@ -357,6 +350,14 @@ next2: *(ckt->CKTstate0 + here->DIOvoltage) = vd;
*(here->DIOnegPosPrimePtr) -= gd;
*(here->DIOposPrimePosPtr) -= gspr;
*(here->DIOposPrimeNegPtr) -= gd;
if (model->DIOsoftRevRecParamGiven) {
//printf("cd = %.7e, cdif = %.7e, gdif = %.7e\n", cd, cdif, gdif);
*(ckt->CKTrhs + here->DIOqdNode) += tt * (cd - vp * cdif) - tt*(gd-gdif)*vd;
*(here->DIOqdQdPtr) += tt*(gd-gdif);
*(here->DIOqdPosPrimePtr) -= tt*(gd-gdif);
*(here->DIOqdNegPtr) -= tt*(gd-gdif);
}
}
}
return(OK);

372
src/spicelib/devices/dio/dioload_separate.c
View File

@ -1,372 +0,0 @@
/**********
Copyright 1990 Regents of the University of California. All rights reserved.
Author: 1985 Thomas L. Quarles
Modified: 2000 AlansFixes
Modified by Paolo Nenzi 2003, Dietmar Warning 2012 and Ste Kulov 2021
**********/
#include "ngspice/ngspice.h"
#include "ngspice/devdefs.h"
#include "ngspice/cktdefs.h"
#include "diodefs.h"
#include "ngspice/const.h"
#include "ngspice/trandefs.h"
#include "ngspice/sperror.h"
#include "ngspice/suffix.h"
int
DIOload(GENmodel *inModel, CKTcircuit *ckt)
/* actually load the current resistance value into the
* sparse matrix previously provided
*/
{
DIOmodel *model = (DIOmodel*)inModel;
DIOinstance *here;
double arg;
double capd;
double cd;
double cdeq;
double cdhat;
double ceq;
double csat; /* area-scaled saturation current */
double czero;
double czof2;
double delvd; /* change in diode voltage temporary */
double evd;
double evrev;
double gd;
double geq;
double gspr; /* area-scaled conductance */
double sarg;
#ifndef NOBYPASS
double tol; /* temporary for tolerence calculations */
#endif
double vd; /* current diode voltage */
double vdtemp;
double vt; /* K t / Q */
double vte, vtebrk;
int Check=0;
int error;
int SenCond=0; /* sensitivity condition */
double diffcharge, deplcharge, diffcap, deplcap;
double tt;
double vp;
/* loop through all the diode models */
for( ; model != NULL; model = DIOnextModel(model)) {
/* loop through all the instances of the model */
for (here = DIOinstances(model); here != NULL ;
here=DIOnextInstance(here)) {
/*
* this routine loads diodes for dc and transient analyses.
*/
if(ckt->CKTsenInfo){
if((ckt->CKTsenInfo->SENstatus == PERTURBATION)
&& (here->DIOsenPertFlag == OFF))continue;
SenCond = here->DIOsenPertFlag;
#ifdef SENSDEBUG
printf("DIOload \n");
#endif /* SENSDEBUG */
}
csat = here->DIOtSatCur;
gspr = here->DIOtConductance;
vt = CONSTKoverQ * here->DIOtemp;
vte = model->DIOemissionCoeff * vt;
vtebrk = model->DIObrkdEmissionCoeff * vt;
tt = here->DIOtTransitTime;
vp = model->DIOsoftRevRecParam;
/*
* initialization
*/
if(SenCond){
#ifdef SENSDEBUG
printf("DIOsenPertFlag = ON \n");
#endif /* SENSDEBUG */
if((ckt->CKTsenInfo->SENmode == TRANSEN)&&
(ckt->CKTmode & MODEINITTRAN)) {
vd = *(ckt->CKTstate1 + here->DIOvoltage);
} else{
vd = *(ckt->CKTstate0 + here->DIOvoltage);
}
#ifdef SENSDEBUG
printf("vd = %.7e \n",vd);
#endif /* SENSDEBUG */
goto next1;
}
Check=1;
if(ckt->CKTmode & MODEINITSMSIG) {
vd= *(ckt->CKTstate0 + here->DIOvoltage);
} else if (ckt->CKTmode & MODEINITTRAN) {
vd= *(ckt->CKTstate1 + here->DIOvoltage);
} else if ( (ckt->CKTmode & MODEINITJCT) &&
(ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC) ) {
vd=here->DIOinitCond;
} else if ( (ckt->CKTmode & MODEINITJCT) && here->DIOoff) {
vd=0;
} else if ( ckt->CKTmode & MODEINITJCT) {
vd=here->DIOtVcrit;
} else if ( ckt->CKTmode & MODEINITFIX && here->DIOoff) {
vd=0;
} else {
#ifndef PREDICTOR
if (ckt->CKTmode & MODEINITPRED) {
*(ckt->CKTstate0 + here->DIOvoltage) =
*(ckt->CKTstate1 + here->DIOvoltage);
vd = DEVpred(ckt,here->DIOvoltage);
*(ckt->CKTstate0 + here->DIOcurrent) =
*(ckt->CKTstate1 + here->DIOcurrent);
*(ckt->CKTstate0 + here->DIOconduct) =
*(ckt->CKTstate1 + here->DIOconduct);
} else {
#endif /* PREDICTOR */
vd = *(ckt->CKTrhsOld+here->DIOposPrimeNode)-
*(ckt->CKTrhsOld + here->DIOnegNode);
#ifndef PREDICTOR
}
#endif /* PREDICTOR */
delvd=vd- *(ckt->CKTstate0 + here->DIOvoltage);
cdhat= *(ckt->CKTstate0 + here->DIOcurrent) +
*(ckt->CKTstate0 + here->DIOconduct) * delvd;
/*
* bypass if solution has not changed
*/
#ifndef NOBYPASS
if ((!(ckt->CKTmode & MODEINITPRED)) && (ckt->CKTbypass)) {
tol=ckt->CKTvoltTol + ckt->CKTreltol*
MAX(fabs(vd),fabs(*(ckt->CKTstate0 +here->DIOvoltage)));
if (fabs(delvd) < tol){
tol=ckt->CKTreltol* MAX(fabs(cdhat),
fabs(*(ckt->CKTstate0 + here->DIOcurrent)))+
ckt->CKTabstol;
if (fabs(cdhat- *(ckt->CKTstate0 + here->DIOcurrent))
< tol) {
vd= *(ckt->CKTstate0 + here->DIOvoltage);
cd= *(ckt->CKTstate0 + here->DIOcurrent);
gd= *(ckt->CKTstate0 + here->DIOconduct);
goto load;
}
}
}
#endif /* NOBYPASS */
/*
* limit new junction voltage
*/
if ( (model->DIObreakdownVoltageGiven) &&
(vd < MIN(0,-here->DIOtBrkdwnV+10*vtebrk))) {
vdtemp = -(vd+here->DIOtBrkdwnV);
vdtemp = DEVpnjlim(vdtemp,
-(*(ckt->CKTstate0 + here->DIOvoltage) +
here->DIOtBrkdwnV),vtebrk,
here->DIOtVcrit,&Check);
vd = -(vdtemp+here->DIOtBrkdwnV);
} else {
vd = DEVpnjlim(vd,*(ckt->CKTstate0 + here->DIOvoltage),
vte,here->DIOtVcrit,&Check);
}
}
/*
* compute dc current and derivitives
*/
next1:
if (vd >= -3*vte) { /* forward */
evd = exp(vd/vte);
cd = csat*(evd-1) + ckt->CKTgmin*vd;
gd = csat*evd/vte + ckt->CKTgmin;
} else if((!(model->DIObreakdownVoltageGiven)) ||
vd >= -here->DIOtBrkdwnV) { /* reverse */
arg = 3*vte/(vd*CONSTe);
arg = arg * arg * arg;
cd = -csat*(1+arg) + ckt->CKTgmin*vd;
gd = csat*3*arg/vd + ckt->CKTgmin;
} else { /* breakdown */
evrev = exp(-(here->DIOtBrkdwnV+vd)/vtebrk);
cd = -csat*evrev + ckt->CKTgmin*vd;
gd = csat*evrev/vtebrk + ckt->CKTgmin;
}
if ((ckt->CKTmode & (MODEDCTRANCURVE | MODETRAN | MODEAC | MODEINITSMSIG)) ||
((ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC))) {
/*
* charge storage elements
*/
czero=here->DIOtJctCap;
if (vd < here->DIOtDepCap){
arg=1-vd/here->DIOtJctPot;
sarg=exp(-here->DIOtGradingCoeff*log(arg));
deplcharge = here->DIOtJctPot*czero*(1-arg*sarg)/(1-here->DIOtGradingCoeff);
deplcap = czero*sarg;
} else {
czof2=czero/here->DIOtF2;
deplcharge = czero*here->DIOtF1+czof2*(here->DIOtF3*(vd-here->DIOtDepCap)+
(here->DIOtGradingCoeff/(here->DIOtJctPot+here->DIOtJctPot))*(vd*vd-here->DIOtDepCap*here->DIOtDepCap));
deplcap = czof2*(here->DIOtF3+here->DIOtGradingCoeff*vd/here->DIOtJctPot);
}
*(ckt->CKTstate0 + here->DIOcapCharge) = deplcharge;
if (model->DIOsoftRevRecParamGiven) {
if (ckt->CKTmode & MODEINITTRAN) {
diffcharge = tt * cd;
diffcap = tt * gd;
*(ckt->CKTstate2 + here->DIOdiffCharge) = diffcharge;
*(ckt->CKTstate2 + here->DIOdiffCap) = diffcap;
*(ckt->CKTstate1 + here->DIOdiffCharge) = diffcharge;
*(ckt->CKTstate1 + here->DIOdiffCap) = diffcap;
*(ckt->CKTstate1 + here->DIOoldCurr) = cd;
*(ckt->CKTstate1 + here->DIOoldCond) = gd;
}
else {
double dt = ckt->CKTdelta;
//Backward Euler
// Qk+1 = tt (VP Qk + h Ik+1) / (tt VP + h)
//diffcharge = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + dt*cd) / (tt*vp + dt));
//diffcap = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + dt*gd) / (tt*vp + dt));
//Trap
// Qk+1 = (Qk (2 tt VP - h) + h tt (Ik + Ik+1)) / (2 tt VP + h)
diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCurr)) + cd)) / (2.0*tt*vp + dt);
diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCond)) + gd)) / (2.0*tt*vp + dt);
//Gear
// Qk+1 = tt ((VP (4 Qk - Qk-1) + 2 h Ik+1) / (3 tt VP + 2 h))
//diffcharge = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) + 2.0*dt*cd) / (3.0*tt*vp + 2.0*dt));
//diffcap = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) + 2.0*dt*gd) / (3.0*tt*vp + 2.0*dt));
}
//printf("time: %e vd: %e\n", ckt->CKTtime, vd);
//printf("%e %e %e %e %e %e\n", ckt->CKTtime, ckt->CKTdelta, diffcharge, diffcap, cd, gd);
} else {
diffcharge = tt*cd;
diffcap = tt*gd;
}
*(ckt->CKTstate0 + here->DIOdiffCharge) = diffcharge;
*(ckt->CKTstate0 + here->DIOdiffCap) = diffcap;
//*(ckt->CKTstate0 + here->DIOoldCurr) = cd;
//*(ckt->CKTstate0 + here->DIOoldCond) = gd;
capd = deplcap + diffcap;
here->DIOcap = capd;
/*
* store small-signal parameters
*/
if( (!(ckt->CKTmode & MODETRANOP)) ||
(!(ckt->CKTmode & MODEUIC)) ) {
if (ckt->CKTmode & MODEINITSMSIG){
*(ckt->CKTstate0 + here->DIOcapCurrent) = capd;
if(SenCond){
*(ckt->CKTstate0 + here->DIOcurrent) = cd;
*(ckt->CKTstate0 + here->DIOconduct) = gd;
#ifdef SENSDEBUG
printf("storing small signal parameters\n");
printf("cd = %.7e,vd = %.7e\n",cd,vd);
printf("capd = %.7e ,gd = %.7e \n",capd,gd);
#endif /* SENSDEBUG */
}
continue;
}
/*
* transient analysis
*/
if(SenCond && (ckt->CKTsenInfo->SENmode == TRANSEN)){
*(ckt->CKTstate0 + here->DIOcurrent) = cd;
#ifdef SENSDEBUG
printf("storing parameters for transient sensitivity\n"
);
printf("qd = %.7e, capd = %.7e,cd = %.7e\n",
*(ckt->CKTstate0 + here->DIOcapCharge),capd,cd);
#endif /* SENSDEBUG */
continue;
}
if (ckt->CKTmode & MODEINITTRAN) {
*(ckt->CKTstate1 + here->DIOcapCharge) =
*(ckt->CKTstate0 + here->DIOcapCharge);
*(ckt->CKTstate1 + here->DIOdiffCharge) =
*(ckt->CKTstate0 + here->DIOdiffCharge);
}
error = NIintegrate(ckt,&geq,&ceq,deplcap,here->DIOcapCharge);
if(error) return(error);
gd=gd+geq;
cd=cd+*(ckt->CKTstate0 + here->DIOcapCurrent);
error = NIintegrate(ckt,&geq,&ceq,diffcap,here->DIOdiffCharge);
if(error) return(error);
gd=gd+geq;
cd=cd+*(ckt->CKTstate0 + here->DIOdiffCurrent);
if (ckt->CKTmode & MODEINITTRAN) {
*(ckt->CKTstate1 + here->DIOcapCurrent) =
*(ckt->CKTstate0 + here->DIOcapCurrent);
*(ckt->CKTstate1 + here->DIOdiffCurrent) =
*(ckt->CKTstate0 + here->DIOdiffCurrent);
}
}
}
if(SenCond) goto next2;
/*
* check convergence
*/
if ( (!(ckt->CKTmode & MODEINITFIX)) || (!(here->DIOoff)) ) {
if (Check == 1) {
ckt->CKTnoncon++;
ckt->CKTtroubleElt = (GENinstance *) here;
}
}
next2: *(ckt->CKTstate0 + here->DIOvoltage) = vd;
*(ckt->CKTstate0 + here->DIOcurrent) = cd;
*(ckt->CKTstate0 + here->DIOconduct) = gd;
if(SenCond) continue;
#ifndef NOBYPASS
load:
#endif
/*
* load current vector
*/
cdeq=cd-gd*vd;
*(ckt->CKTrhs + here->DIOnegNode) += cdeq;
*(ckt->CKTrhs + here->DIOposPrimeNode) -= cdeq;
/*
* load matrix
*/
*(here->DIOposPosPtr) += gspr;
*(here->DIOnegNegPtr) += gd;
*(here->DIOposPrimePosPrimePtr) += (gd + gspr);
*(here->DIOposPosPrimePtr) -= gspr;
*(here->DIOnegPosPrimePtr) -= gd;
*(here->DIOposPrimePosPtr) -= gspr;
*(here->DIOposPrimeNegPtr) -= gd;
}
}
return(OK);
}

430
src/spicelib/devices/dio/dioload_sk.c
View File

@ -1,430 +0,0 @@
/**********
Copyright 1990 Regents of the University of California. All rights reserved.
Author: 1985 Thomas L. Quarles
Modified: 2000 AlansFixes
Modified by Paolo Nenzi 2003 and Dietmar Warning 2012
**********/
#include "ngspice/ngspice.h"
#include "ngspice/devdefs.h"
#include "ngspice/cktdefs.h"
#include "diodefs.h"
#include "ngspice/const.h"
#include "ngspice/trandefs.h"
#include "ngspice/sperror.h"
#include "ngspice/suffix.h"
int
DIOload(GENmodel *inModel, CKTcircuit *ckt)
/* actually load the current resistance value into the
* sparse matrix previously provided
*/
{
DIOmodel *model = (DIOmodel*)inModel;
DIOinstance *here;
double arg;
double capd;
double cd;
double cdeq;
double cdhat;
double ceq;
double csat; /* area-scaled saturation current */
double czero;
double czof2;
double delvd; /* change in diode voltage temporary */
double evd;
double evrev;
double gd;
double geq;
double gspr; /* area-scaled conductance */
double sarg;
#ifndef NOBYPASS
double tol; /* temporary for tolerence calculations */
#endif
double vd; /* current diode voltage */
double vdtemp;
double vt; /* K t / Q */
double vte, vtebrk;
int Check=0;
int error;
int SenCond=0; /* sensitivity condition */
double diffcharge, deplcharge, diffcap, deplcap;
double tt;
double vp;
/* loop through all the diode models */
for( ; model != NULL; model = DIOnextModel(model)) {
/* loop through all the instances of the model */
for (here = DIOinstances(model); here != NULL ;
here=DIOnextInstance(here)) {
/*
* this routine loads diodes for dc and transient analyses.
*/
if(ckt->CKTsenInfo){
if((ckt->CKTsenInfo->SENstatus == PERTURBATION)
&& (here->DIOsenPertFlag == OFF))continue;
SenCond = here->DIOsenPertFlag;
#ifdef SENSDEBUG
printf("DIOload \n");
#endif /* SENSDEBUG */
}
csat = here->DIOtSatCur;
gspr = here->DIOtConductance;
vt = CONSTKoverQ * here->DIOtemp;
vte = model->DIOemissionCoeff * vt;
vtebrk = model->DIObrkdEmissionCoeff * vt;
tt = here->DIOtTransitTime;
vp = model->DIOsoftRevRecParam;
/*
* initialization
*/
if(SenCond){
#ifdef SENSDEBUG
printf("DIOsenPertFlag = ON \n");
#endif /* SENSDEBUG */
if((ckt->CKTsenInfo->SENmode == TRANSEN)&&
(ckt->CKTmode & MODEINITTRAN)) {
vd = *(ckt->CKTstate1 + here->DIOvoltage);
} else{
vd = *(ckt->CKTstate0 + here->DIOvoltage);
}
#ifdef SENSDEBUG
printf("vd = %.7e \n",vd);
#endif /* SENSDEBUG */
goto next1;
}
Check=1;
if(ckt->CKTmode & MODEINITSMSIG) {
vd= *(ckt->CKTstate0 + here->DIOvoltage);
} else if (ckt->CKTmode & MODEINITTRAN) {
vd= *(ckt->CKTstate1 + here->DIOvoltage);
} else if ( (ckt->CKTmode & MODEINITJCT) &&
(ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC) ) {
vd=here->DIOinitCond;
} else if ( (ckt->CKTmode & MODEINITJCT) && here->DIOoff) {
vd=0;
} else if ( ckt->CKTmode & MODEINITJCT) {
vd=here->DIOtVcrit;
} else if ( ckt->CKTmode & MODEINITFIX && here->DIOoff) {
vd=0;
} else {
#ifndef PREDICTOR
if (ckt->CKTmode & MODEINITPRED) {
*(ckt->CKTstate0 + here->DIOvoltage) =
*(ckt->CKTstate1 + here->DIOvoltage);
vd = DEVpred(ckt,here->DIOvoltage);
*(ckt->CKTstate0 + here->DIOcurrent) =
*(ckt->CKTstate1 + here->DIOcurrent);
*(ckt->CKTstate0 + here->DIOconduct) =
*(ckt->CKTstate1 + here->DIOconduct);
} else {
#endif /* PREDICTOR */
vd = *(ckt->CKTrhsOld+here->DIOposPrimeNode)-
*(ckt->CKTrhsOld + here->DIOnegNode);
#ifndef PREDICTOR
}
#endif /* PREDICTOR */
delvd=vd- *(ckt->CKTstate0 + here->DIOvoltage);
cdhat= *(ckt->CKTstate0 + here->DIOcurrent) +
*(ckt->CKTstate0 + here->DIOconduct) * delvd;
/*
* bypass if solution has not changed
*/
#ifndef NOBYPASS
if ((!(ckt->CKTmode & MODEINITPRED)) && (ckt->CKTbypass)) {
tol=ckt->CKTvoltTol + ckt->CKTreltol*
MAX(fabs(vd),fabs(*(ckt->CKTstate0 +here->DIOvoltage)));
if (fabs(delvd) < tol){
tol=ckt->CKTreltol* MAX(fabs(cdhat),
fabs(*(ckt->CKTstate0 + here->DIOcurrent)))+
ckt->CKTabstol;
if (fabs(cdhat- *(ckt->CKTstate0 + here->DIOcurrent))
< tol) {
vd= *(ckt->CKTstate0 + here->DIOvoltage);
cd= *(ckt->CKTstate0 + here->DIOcurrent);
gd= *(ckt->CKTstate0 + here->DIOconduct);
goto load;
}
}
}
#endif /* NOBYPASS */
/*
* limit new junction voltage
*/
if ( (model->DIObreakdownVoltageGiven) &&
(vd < MIN(0,-here->DIOtBrkdwnV+10*vtebrk))) {
vdtemp = -(vd+here->DIOtBrkdwnV);
vdtemp = DEVpnjlim(vdtemp,
-(*(ckt->CKTstate0 + here->DIOvoltage) +
here->DIOtBrkdwnV),vtebrk,
here->DIOtVcrit,&Check);
vd = -(vdtemp+here->DIOtBrkdwnV);
} else {
vd = DEVpnjlim(vd,*(ckt->CKTstate0 + here->DIOvoltage),
vte,here->DIOtVcrit,&Check);
}
}
/*
* compute dc current and derivitives
*/
next1:
if (vd >= -3*vte) {
evd = exp(vd/vte);
cd = csat*(evd-1);
gd = csat*evd/vte;
} else if((!(model->DIObreakdownVoltageGiven)) ||
vd >= -here->DIOtBrkdwnV) {
arg = 3*vte/(vd*CONSTe);
arg = arg * arg * arg;
cd = -csat*(1+arg);
gd = csat*3*arg/vd;
} else {
evrev = exp(-(here->DIOtBrkdwnV+vd)/vtebrk);
cd = -csat*evrev;
gd = csat*evrev/vtebrk;
}
/*
if (ckt->CKTmode & MODEINITTRAN) {
*(ckt->CKTstate1 + here->DIOoldCurr) = cd;
*(ckt->CKTstate1 + here->DIOoldCond) = gd;
*(ckt->CKTstate2 + here->DIOoldCurr) = cd;
*(ckt->CKTstate2 + here->DIOoldCond) = gd;
}
else {
//cd = (*(ckt->CKTstate1 + here->DIOoldCurr) - (ckt->CKTdelta)*cd) / (1 + (ckt->CKTdelta)*vp);
//gd = (*(ckt->CKTstate1 + here->DIOoldCond) - (ckt->CKTdelta)*gd) / (1 + (ckt->CKTdelta)*vp);
//cd = (1 - vp)*cd + vp*(*(ckt->CKTstate1 + here->DIOoldCurr))*(ckt->CKTdelta);
//gd = (1 - vp)*gd + vp*(*(ckt->CKTstate1 + here->DIOoldCond))*(ckt->CKTdelta);
}
*/
if ((ckt->CKTmode & (MODEDCTRANCURVE | MODETRAN | MODEAC | MODEINITSMSIG)) ||
((ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC))) {
/*
* charge storage elements
*/
/* junction(depletion) charge */
czero=here->DIOtJctCap;
if (vd < here->DIOtDepCap){
arg=1-vd/here->DIOtJctPot;
sarg=exp(-here->DIOtGradingCoeff*log(arg));
deplcharge = here->DIOtJctPot*czero*(1-arg*sarg)/(1-here->DIOtGradingCoeff);
deplcap = czero*sarg;
} else {
czof2=czero/here->DIOtF2;
deplcharge = czero*here->DIOtF1+czof2*(here->DIOtF3*(vd-here->DIOtDepCap)+
(here->DIOtGradingCoeff/(here->DIOtJctPot+here->DIOtJctPot))*(vd*vd-here->DIOtDepCap*here->DIOtDepCap));
deplcap = czof2*(here->DIOtF3+here->DIOtGradingCoeff*vd/here->DIOtJctPot);
}
/* diffusion charge */
if (ckt->CKTmode & MODEINITTRAN) {
diffcharge = tt * cd;
diffcap = tt * gd;
*(ckt->CKTstate2 + here->DIOdiffCharge) = diffcharge;
*(ckt->CKTstate2 + here->DIOdiffCap) = diffcap;
*(ckt->CKTstate1 + here->DIOdiffCharge) = diffcharge;
*(ckt->CKTstate1 + here->DIOdiffCap) = diffcap;
*(ckt->CKTstate1 + here->DIOoldCurr) = cd;
*(ckt->CKTstate1 + here->DIOoldCond) = gd;
}
else {
double dt = ckt->CKTdelta;
/*
//Fixed-point iteration
double diffcharge_old, diffcap_old;
//Forward Euler predictor
diffcharge = (*(ckt->CKTstate1 + here->DIOdiffCharge))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCurr)) / vp);
diffcap = (*(ckt->CKTstate1 + here->DIOdiffCap))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCond)) / vp);
do {
diffcharge_old = diffcharge;
diffcharge = (*(ckt->CKTstate1 + here->DIOdiffCharge)) + ((dt / vp)*(cd - (diffcharge_old / tt)));
} while (fabs(diffcharge - diffcharge_old) > 1e-12);
do {
diffcap_old = diffcap;
diffcap = (*(ckt->CKTstate1 + here->DIOdiffCap)) + ((dt / vp)*(gd - (diffcap_old / tt)));
} while (fabs(diffcap - diffcap_old) > 1e-12);
*/
/*
//Newton's method
double diffcharge_old, diffcap_old;
//Forward Euler predictor
diffcharge = (*(ckt->CKTstate1 + here->DIOdiffCharge))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCurr)) / vp);
diffcap = (*(ckt->CKTstate1 + here->DIOdiffCap))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCond)) / vp);
do {
diffcharge_old = diffcharge;
diffcharge = diffcharge_old + ((tt*vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + tt*dt*cd - diffcharge_old*(dt + tt*vp)) / (dt + tt*vp));
} while (fabs(diffcharge - diffcharge_old) > 1e-12);
do {
diffcap_old = diffcap;
diffcap = diffcap_old + ((tt*vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + tt * dt*gd - diffcap_old * (dt + tt*vp)) / (dt + tt*vp));
} while (fabs(diffcap - diffcap_old) > 1e-12);
*/
//Backward Difference Operator
//diffcharge = tt * (cd - vp * (((*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) / (ckt->CKTdeltaOld[1])));
//diffcap = tt * (gd - vp * (((*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) / (ckt->CKTdeltaOld[1])));
//Standard Diffy-Q w/ trap integral
//diffcharge = (dt / (2.0*vt))*(exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOoldCurr)) + cd) + exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOdiffCharge));
//diffcap = (dt / (2.0*vt))*(exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOoldCond)) + gd) + exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOdiffCap));
//Forward Euler
//diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(tt*vp - dt) + dt*tt*cd) / (tt*vp);
//diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(tt*vp - dt) + dt*tt*gd) / (tt*vp);
//Backward Euler
diffcharge = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + dt*cd) / (tt*vp + dt));
diffcap = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + dt*gd) / (tt*vp + dt));
//Trap
//diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCurr)) + cd)) / (2.0*tt*vp + dt);
//diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCond)) + gd)) / (2.0*tt*vp + dt);
//Gear
//diffcharge = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) + 2.0*dt*cd) / (3.0*tt*vp + 2.0*dt));
//diffcap = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) + 2.0*dt*gd) / (3.0*tt*vp + 2.0*dt));
//Nothing
//diffcharge = tt * cd;
//diffcap = tt * gd;
}
*(ckt->CKTstate0 + here->DIOdiffCharge) = diffcharge;
*(ckt->CKTstate0 + here->DIOdiffCap) = diffcap;
//*(ckt->CKTstate0 + here->DIOoldCurr) = cd;
//*(ckt->CKTstate0 + here->DIOoldCond) = gd;
//printf("Time: %e\n", ckt->CKTtime);
//printf("Charge: %e\n", diffcharge);
//printf("Cap: %e\n", diffcap);
//printf("Curr: %e\n", cd);
//printf("Cond: %e\n", gd);
//printf("DelCharge: %e\n", diffcharge-(*(ckt->CKTstate1 + here->DIOdiffCharge)));
//printf("DelCap: %e\n", diffcap-(*(ckt->CKTstate1 + here->DIOdiffCap)));
//printf("DelCurr: %e\n", cd-(*(ckt->CKTstate1 + here->DIOoldCurr)));
//printf("DelCond: %e\n\n", gd-(*(ckt->CKTstate1 + here->DIOoldCond)));
//diffcharge = here->DIOtTransitTime*cd;
*(ckt->CKTstate0 + here->DIOcapCharge) =
diffcharge + deplcharge;
//diffcap = here->DIOtTransitTime*gd;
capd = diffcap + deplcap + here->DIOcmetal + here->DIOcpoly;
here->DIOcap = capd;
/*
* store small-signal parameters
*/
if( (!(ckt->CKTmode & MODETRANOP)) ||
(!(ckt->CKTmode & MODEUIC)) ) {
if (ckt->CKTmode & MODEINITSMSIG){
*(ckt->CKTstate0 + here->DIOcapCurrent) = capd;
if(SenCond){
*(ckt->CKTstate0 + here->DIOcurrent) = cd;
*(ckt->CKTstate0 + here->DIOconduct) = gd;
#ifdef SENSDEBUG
printf("storing small signal parameters\n");
printf("cd = %.7e,vd = %.7e\n",cd,vd);
printf("capd = %.7e ,gd = %.7e \n",capd,gd);
#endif /* SENSDEBUG */
}
continue;
}
/*
* transient analysis
*/
if(SenCond && (ckt->CKTsenInfo->SENmode == TRANSEN)){
*(ckt->CKTstate0 + here->DIOcurrent) = cd;
#ifdef SENSDEBUG
printf("storing parameters for transient sensitivity\n"
);
printf("qd = %.7e, capd = %.7e,cd = %.7e\n",
*(ckt->CKTstate0 + here->DIOcapCharge),capd,cd);
#endif /* SENSDEBUG */
continue;
}
if (ckt->CKTmode & MODEINITTRAN) {
*(ckt->CKTstate1 + here->DIOcapCharge) =
*(ckt->CKTstate0 + here->DIOcapCharge);
}
error = NIintegrate(ckt,&geq,&ceq,capd,here->DIOcapCharge);
if(error) return(error);
gd=gd+geq;
cd=cd+*(ckt->CKTstate0 + here->DIOcapCurrent);
if (ckt->CKTmode & MODEINITTRAN) {
*(ckt->CKTstate1 + here->DIOcapCurrent) =
*(ckt->CKTstate0 + here->DIOcapCurrent);
}
}
}
if(SenCond) goto next2;
/*
* check convergence
*/
if ( (!(ckt->CKTmode & MODEINITFIX)) || (!(here->DIOoff)) ) {
if (Check == 1) {
ckt->CKTnoncon++;
ckt->CKTtroubleElt = (GENinstance *) here;
}
}
next2: *(ckt->CKTstate0 + here->DIOvoltage) = vd;
*(ckt->CKTstate0 + here->DIOcurrent) = cd;
*(ckt->CKTstate0 + here->DIOconduct) = gd;
if(SenCond) continue;
#ifndef NOBYPASS
load:
#endif
/*
* load current vector
*/
cdeq=cd-gd*vd;
*(ckt->CKTrhs + here->DIOnegNode) += cdeq;
*(ckt->CKTrhs + here->DIOposPrimeNode) -= cdeq;
/*
* load matrix
*/
*(here->DIOposPosPtr) += gspr;
*(here->DIOnegNegPtr) += gd;
*(here->DIOposPrimePosPrimePtr) += (gd + gspr);
*(here->DIOposPosPrimePtr) -= gspr;
*(here->DIOnegPosPrimePtr) -= gd;
*(here->DIOposPrimePosPtr) -= gspr;
*(here->DIOposPrimeNegPtr) -= gd;
}
}
return(OK);
}

24
src/spicelib/devices/dio/diosetup.c
View File

@ -319,6 +319,17 @@ DIOsetup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *states)
}
}
/* rev-rec */
if (model->DIOsoftRevRecParamGiven) {
if(here->DIOqdNode == 0) {
error = CKTmkVolt(ckt, &tmp, here->DIOname, "Qdiff");
if(error) return(error);
here->DIOqdNode = tmp->number;
}
} else {
here->DIOqdNode = 0;
}
int selfheat = ((here->DIOtempNode > 0) && (here->DIOthermal) && (model->DIOrth0Given));
/* macro to make elements with built in test for out of memory */
@ -345,6 +356,13 @@ do { if((here->ptr = SMPmakeElt(matrix, here->first, here->second)) == NULL){\
TSTALLOC(DIOnegTempPtr, DIOnegNode, DIOtempNode);
}
/* rev-rec */
if (model->DIOsoftRevRecParamGiven) {
TSTALLOC(DIOqdQdPtr , DIOqdNode, DIOqdNode);
TSTALLOC(DIOqdPosPrimePtr, DIOqdNode, DIOposPrimeNode);
TSTALLOC(DIOqdNegPtr , DIOqdNode, DIOnegNode);
}
}
}
return(OK);
@ -369,6 +387,12 @@ DIOunsetup(
&& here->DIOposPrimeNode != here->DIOposNode)
CKTdltNNum(ckt, here->DIOposPrimeNode);
here->DIOposPrimeNode = 0;
/* rev-rec */
if (model->DIOsoftRevRecParamGiven) {
CKTdltNNum(ckt, here->DIOqdNode);
here->DIOqdNode = 0;
}
}
}
return OK;

2
src/spicelib/devices/dio/diotrunc.c
View File

@ -21,7 +21,7 @@ DIOtrunc(GENmodel *inModel, CKTcircuit *ckt, double *timeStep)
for( ; model != NULL; model = DIOnextModel(model)) {
for(here=DIOinstances(model);here!=NULL;here = DIOnextInstance(here)){
CKTterr(here->DIOcapCharge,ckt,timeStep);
CKTterr(here->DIOdiffCharge,ckt,timeStep);
CKTterr(here->DIOdifCharge,ckt,timeStep);
}
}
return(OK);