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

simplify dioload

This commit is contained in:
dwarning 2025-09-01 11:15:36 +02:00
parent 64a433ebaa
commit 3161e06584
3 changed files with 544 additions and 181 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -2,7 +2,7 @@
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
Modified by Paolo Nenzi 2003, Dietmar Warning 2012 and Ste Kulov 2021
**********/
#include "ngspice/ngspice.h"
@ -31,7 +31,6 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
double csat; /* area-scaled saturation current */
double czero;
double czof2;
double delvd; /* change in diode voltage temporary */
double evd;
double evrev;
@ -50,6 +49,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
int error;
int SenCond=0; /* sensitivity condition */
double diffcharge, deplcharge, diffcap, deplcap;
double tt;
double vp;
@ -81,6 +81,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
vtebrk = model->DIObrkdEmissionCoeff * vt;
tt = here->DIOtTransitTime;
vp = model->DIOsoftRevRecParam;
/*
* initialization
*/
@ -151,10 +152,10 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
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;
vd= *(ckt->CKTstate0 + here->DIOvoltage);
cd= *(ckt->CKTstate0 + here->DIOcurrent);
gd= *(ckt->CKTstate0 + here->DIOconduct);
goto load;
}
}
}
@ -178,52 +179,34 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
/*
* compute dc current and derivitives
*/
next1:
if (vd >= -3*vte) {
next1:
if (vd >= -3*vte) { /* forward */
evd = exp(vd/vte);
cd = csat*(evd-1);
gd = csat*evd/vte;
cd = csat*(evd-1) + ckt->CKTgmin*vd;
gd = csat*evd/vte + ckt->CKTgmin;
} else if((!(model->DIObreakdownVoltageGiven)) ||
vd >= -here->DIOtBrkdwnV) {
vd >= -here->DIOtBrkdwnV) { /* reverse */
arg = 3*vte/(vd*CONSTe);
arg = arg * arg * arg;
cd = -csat*(1+arg);
gd = csat*3*arg/vd;
cd = -csat*(1+arg) + ckt->CKTgmin*vd;
gd = csat*3*arg/vd + ckt->CKTgmin;
} else {
} else { /* breakdown */
evrev = exp(-(here->DIOtBrkdwnV+vd)/vtebrk);
cd = -csat*evrev;
gd = csat*evrev/vtebrk;
cd = -csat*evrev + ckt->CKTgmin*vd;
gd = csat*evrev/vtebrk + ckt->CKTgmin;
}
/*
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;
@ -237,105 +220,55 @@ next1:
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;
if (model->DIOsoftRevRecParamGiven) {
//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);
*/
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;
/*
//Newton's method
double diffcharge_old, diffcap_old;
//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));
//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);
//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);
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);
*/
//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));
//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));
//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);
//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);
} else {
//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));
diffcharge = tt*cd;
diffcap = tt*gd;
//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->DIOcapCharge) = deplcharge + diffcharge;
*(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;
capd = deplcap + diffcap;
here->DIOcap = capd;
@ -394,7 +327,7 @@ next1:
* check convergence
*/
if ( (!(ckt->CKTmode & MODEINITFIX)) || (!(here->DIOoff)) ) {
if (Check == 1) {
if (Check == 1) {
ckt->CKTnoncon++;
ckt->CKTtroubleElt = (GENinstance *) here;
}

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

@ -0,0 +1,372 @@
/**********
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);
}

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

@ -2,7 +2,7 @@
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
Modified by Paolo Nenzi 2003 and Dietmar Warning 2012
**********/
#include "ngspice/ngspice.h"
@ -23,6 +23,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
DIOmodel *model = (DIOmodel*)inModel;
DIOinstance *here;
double arg;
double capd;
double cd;
double cdeq;
double cdhat;
@ -30,6 +31,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
double csat; /* area-scaled saturation current */
double czero;
double czof2;
double delvd; /* change in diode voltage temporary */
double evd;
double evrev;
@ -48,15 +50,12 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
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)) {
vp = model->DIOsoftRevRecParam;
/* loop through all the instances of the model */
for (here = DIOinstances(model); here != NULL ;
here=DIOnextInstance(here)) {
@ -81,6 +80,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
vte = model->DIOemissionCoeff * vt;
vtebrk = model->DIObrkdEmissionCoeff * vt;
tt = here->DIOtTransitTime;
vp = model->DIOsoftRevRecParam;
/*
* initialization
*/
@ -151,10 +151,10 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
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;
vd= *(ckt->CKTstate0 + here->DIOvoltage);
cd= *(ckt->CKTstate0 + here->DIOcurrent);
gd= *(ckt->CKTstate0 + here->DIOconduct);
goto load;
}
}
}
@ -178,22 +178,23 @@ 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);
gd = csat*evd/vte;
} else if((!(model->DIObreakdownVoltageGiven)) ||
vd >= -here->DIOtBrkdwnV) { /* reverse */
vd >= -here->DIOtBrkdwnV) {
arg = 3*vte/(vd*CONSTe);
arg = arg * arg * arg;
cd = -csat*(1+arg);
gd = csat*3*arg/vd;
} else { /* breakdown */
} else {
evrev = exp(-(here->DIOtBrkdwnV+vd)/vtebrk);
cd = -csat*evrev;
@ -201,14 +202,28 @@ next1:
}
gd = gd + ckt->CKTgmin;
cd = cd + ckt->CKTgmin*vd;
/*
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;
@ -222,56 +237,107 @@ next1:
deplcap = czof2*(here->DIOtF3+here->DIOtGradingCoeff*vd/here->DIOtJctPot);
}
if (model->DIOsoftRevRecParamGiven) {
/* 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;
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 = 1*ckt->CKTdelta;
//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);
*/
//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));
/*
//Newton's method
double diffcharge_old, diffcap_old;
//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);
//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);
//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));
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])));
*(ckt->CKTstate0 + here->DIOdiffCap) = diffcap;
*(ckt->CKTstate0 + here->DIOoldCurr) = cd;
*(ckt->CKTstate0 + here->DIOoldCond) = gd;
//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));
//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);
//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);
} else {
//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));
diffcharge = tt*cd;
diffcap = tt*gd;
//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->DIOcapCharge) = deplcharge;
*(ckt->CKTstate0 + here->DIOdiffCharge) = diffcharge;
*(ckt->CKTstate0 + here->DIOdiffCap) = diffcap;
//*(ckt->CKTstate0 + here->DIOoldCurr) = cd;
//*(ckt->CKTstate0 + here->DIOoldCond) = gd;
here->DIOcap = deplcap + diffcap;
//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
@ -279,7 +345,7 @@ next1:
if( (!(ckt->CKTmode & MODETRANOP)) ||
(!(ckt->CKTmode & MODEUIC)) ) {
if (ckt->CKTmode & MODEINITSMSIG){
*(ckt->CKTstate0 + here->DIOcapCurrent) = deplcap + diffcap;
*(ckt->CKTstate0 + here->DIOcapCurrent) = capd;
if(SenCond){
*(ckt->CKTstate0 + here->DIOcurrent) = cd;
@ -287,7 +353,7 @@ next1:
#ifdef SENSDEBUG
printf("storing small signal parameters\n");
printf("cd = %.7e,vd = %.7e\n",cd,vd);
printf("capd = %.7e ,gd = %.7e \n",deplcap+diffcap,gd);
printf("capd = %.7e ,gd = %.7e \n",capd,gd);
#endif /* SENSDEBUG */
}
continue;
@ -302,7 +368,7 @@ next1:
printf("storing parameters for transient sensitivity\n"
);
printf("qd = %.7e, capd = %.7e,cd = %.7e\n",
*(ckt->CKTstate0 + here->DIOcapCharge),deplcap+diffcap,cd);
*(ckt->CKTstate0 + here->DIOcapCharge),capd,cd);
#endif /* SENSDEBUG */
continue;
}
@ -310,22 +376,14 @@ next1:
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);
error = NIintegrate(ckt,&geq,&ceq,capd,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);
}
}
}
@ -336,7 +394,7 @@ next1:
* check convergence
*/
if ( (!(ckt->CKTmode & MODEINITFIX)) || (!(here->DIOoff)) ) {
if (Check == 1) {
if (Check == 1) {
ckt->CKTnoncon++;
ckt->CKTtroubleElt = (GENinstance *) here;
}