diff --git a/src/spicelib/devices/dio/diodefs.h b/src/spicelib/devices/dio/diodefs.h index 82894cc34..8a5137c9d 100644 --- a/src/spicelib/devices/dio/diodefs.h +++ b/src/spicelib/devices/dio/diodefs.h @@ -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 */ diff --git a/src/spicelib/devices/dio/dioload.c b/src/spicelib/devices/dio/dioload.c index dcc1a2554..2e50ddfd2 100644 --- a/src/spicelib/devices/dio/dioload.c +++ b/src/spicelib/devices/dio/dioload.c @@ -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); diff --git a/src/spicelib/devices/dio/dioload_separate.c b/src/spicelib/devices/dio/dioload_separate.c deleted file mode 100644 index 723afb089..000000000 --- a/src/spicelib/devices/dio/dioload_separate.c +++ /dev/null @@ -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); -} diff --git a/src/spicelib/devices/dio/dioload_sk.c b/src/spicelib/devices/dio/dioload_sk.c deleted file mode 100644 index 99fdc5add..000000000 --- a/src/spicelib/devices/dio/dioload_sk.c +++ /dev/null @@ -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); -} diff --git a/src/spicelib/devices/dio/diosetup.c b/src/spicelib/devices/dio/diosetup.c index 49723d8d7..31adbe132 100644 --- a/src/spicelib/devices/dio/diosetup.c +++ b/src/spicelib/devices/dio/diosetup.c @@ -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; diff --git a/src/spicelib/devices/dio/diotrunc.c b/src/spicelib/devices/dio/diotrunc.c index 548acfc75..3caac55be 100644 --- a/src/spicelib/devices/dio/diotrunc.c +++ b/src/spicelib/devices/dio/diotrunc.c @@ -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);