Diode model with selfheating option

src/frontend/inpcom.c

@@ -4277,7 +4277,17 @@ int get_number_terminals(char *c)
         case 'v':
         case 'i':
         case 'd':
-            return 2;
+            i = 0;
+            /* find the first token with "off" or "=" in the line*/
+            while ((i < 10) && (*c != '\0')) {
+                char *inst = gettok_instance(&c);
+                strncpy(nam_buf, inst, sizeof(nam_buf) - 1);
+                txfree(inst);
+                if (strstr(nam_buf, "off") || strstr(nam_buf, "thermal") || strchr(nam_buf, '='))
+                    break;
+                i++;
+            }
+            return i - 2;
             break;
         case 'u':
         case 'j':

src/frontend/subckt.c

@@ -1648,7 +1648,7 @@ numnodes(const char *line, struct subs *subs, wordlist const *modnames)
 
     /* If model names equal node names, this code will fail! */
 
-    if ((c == 'm') || (c == 'p') || (c == 'q')) { /* IF this is a mos, cpl or bjt*/
+    if ((c == 'm') || (c == 'p') || (c == 'q') || (c == 'd')) { /* IF this is a mos, cpl, bjt or diode */
         char *s = nexttok(line);       /* Skip the instance name */
         int gotit = 0;
         int i = 0;
@@ -1928,7 +1928,27 @@ devmodtranslate(struct card *s, char *subname, wordlist * const orig_modnames)
             name = gettok_node(&t);  /* get second attached netname */
             bxx_printf(&buffer, "%s ", name);
             tfree(name);
-            name = gettok(&t);
+            name = gettok_node(&t);  /* this can be either a model name or a node name. */
+            if (name == NULL) {
+                name = copy(""); /* allow 'tfree' */
+            } else {
+                found = 0;
+                while (!found) {
+                    wlsub = wl_find(name, orig_modnames);
+                    if (wlsub) {
+                        found = 1;
+                        break;
+                    } else {
+                        bxx_printf(&buffer, "%s ", name);
+                        tfree(name);
+                        name = gettok(&t);
+                        if (name == NULL) {  /* No token anymore - leave */
+                            name = copy(""); /* allow 'tfree' */
+                            break;
+                        }
+                    }
+                }  /* while  */
+            }
 
             translate_mod_name(&buffer, name, subname, orig_modnames);
 
@@ -2162,7 +2182,7 @@ inp_numnodes(char c)
     case 'c':
         return (2);
     case 'd':
-        return (2);
+        return (3);
     case 'e':
         return (2); /* changed from 4 to 2 by SDB on 4.22.2003 to enable POLY */
     case 'f':

src/spicelib/devices/dio/dio.c

@@ -23,6 +23,7 @@ IFparm DIOpTable[] = { /* parameters */
  IOPU("lp",    DIO_LP,    IF_REAL, "Length of polysilicon capacitor (level=3)"),
  IOPU("wm",    DIO_WM,    IF_REAL, "Width of metal capacitor (level=3)"),
  IOPU("wp",    DIO_WP,    IF_REAL, "Width of polysilicon capacitor (level=3)"),
+ IOP("thermal",DIO_THERMAL, IF_FLAG, "Self heating mode selector"),
 
  IP("sens_area",DIO_AREA_SENS,IF_FLAG,"flag to request sensitivity WRT area"),
  OP("vd",      DIO_VOLTAGE,IF_REAL, "Diode voltage"),
@@ -112,7 +113,10 @@ IFparm DIOmPTable[] = { /* model parameters */
  IOP( "bv_max",   DIO_MOD_BV_MAX,  IF_REAL, "maximum voltage in reverse direction"),
  IOP( "isr",  DIO_MOD_ISR,  IF_REAL, "Recombination saturation current"),
  IOP( "nr",   DIO_MOD_NR,   IF_REAL, "Recombination current emission coefficient"),
-
+/* self heating */
+ IOP("rth0",  DIO_MOD_RTH0,  IF_REAL, "Self-heating thermal resistance"),
+ IOP("cth0",  DIO_MOD_CTH0,  IF_REAL, "Self-heating thermal capacitance"),
+/* scaled parasitic capacitances level 3 model */
  IOP( "lm",  DIO_MOD_LM,  IF_REAL, "Length of metal capacitor (level=3)"),
  IOP( "lp",  DIO_MOD_LP,  IF_REAL, "Length of polysilicon capacitor (level=3)"),
  IOP( "wm",  DIO_MOD_WM,  IF_REAL, "Width of metal capacitor (level=3)"),
@@ -127,7 +131,8 @@ IFparm DIOmPTable[] = { /* model parameters */
 
 char *DIOnames[] = {
     "D+",
-    "D-"
+    "D-",
+    "Tj"
 };
 
 int DIOnSize = NUMELEMS(DIOnames);

src/spicelib/devices/dio/dioacld.c

@@ -28,7 +28,7 @@ DIOacLoad(GENmodel *inModel, CKTcircuit *ckt)
         /* loop through all the instances of the model */
         for (here = DIOinstances(model); here != NULL ;
                 here=DIOnextInstance(here)) {
-            gspr=here->DIOtConductance*here->DIOarea;
+            gspr=here->DIOtConductance;
             geq= *(ckt->CKTstate0 + here->DIOconduct);
             xceq= *(ckt->CKTstate0 + here->DIOcapCurrent) * ckt->CKTomega;
             *(here->DIOposPosPtr ) += gspr;
@@ -42,6 +42,25 @@ DIOacLoad(GENmodel *inModel, CKTcircuit *ckt)
             *(here->DIOposPrimePosPtr ) -= gspr;
             *(here->DIOposPrimeNegPtr ) -= geq;
             *(here->DIOposPrimeNegPtr +1 ) -= xceq;
+            int selfheat = ((here->DIOtempNode > 0) && (here->DIOthermal) && (model->DIOrth0Given));
+            if (selfheat) {
+                double dIth_dVrs = here->DIOdIth_dVrs;
+                double dIth_dVdio = here->DIOdIth_dVdio;
+                double dIth_dT = here->DIOdIth_dT;
+                double gcTt = here->DIOgcTt;
+                double dIrs_dT = here->DIOdIrs_dT;
+                double dIdio_dT = *(ckt->CKTstate0 + here->DIOdIdio_dT);
+                (*(here->DIOtempPosPtr)      += -dIth_dVrs);
+                (*(here->DIOtempPosPrimePtr) += -dIth_dVdio + dIth_dVrs);
+                (*(here->DIOtempNegPtr)      +=  dIth_dVdio);
+                (*(here->DIOtempTempPtr)     += -dIth_dT + 1/model->DIOrth0 + gcTt);
+                (*(here->DIOposTempPtr)      +=  dIrs_dT);
+                (*(here->DIOposPrimeTempPtr) +=  dIdio_dT - dIrs_dT);
+                (*(here->DIOnegTempPtr)      += -dIdio_dT);
+
+                double xgcTt= *(ckt->CKTstate0 + here->DIOcqth) * ckt->CKTomega;
+                (*(here->DIOtempTempPtr + 1) +=  xgcTt);
+            }
         }
     }
     return(OK);

src/spicelib/devices/dio/dioask.c

@@ -48,6 +48,9 @@ DIOask (CKTcircuit *ckt, GENinstance *inst, int which, IFvalue *value,
         case DIO_M:
             value->rValue = here->DIOm;
             return(OK);
+        case DIO_THERMAL:
+            value->iValue = here->DIOthermal;
+            return(OK);
 
         case DIO_TEMP:
             value->rValue = here->DIOtemp-CONSTCtoK;
@@ -81,7 +84,9 @@ DIOask (CKTcircuit *ckt, GENinstance *inst, int which, IFvalue *value,
                 return(E_ASKPOWER);
             } else {
                 value->rValue = *(ckt->CKTstate0 + here->DIOcurrent) *
-                        *(ckt->CKTstate0 + here->DIOvoltage);
+                        *(ckt->CKTstate0 + here->DIOvoltage) +
+                        *(ckt->CKTstate0 + here->DIOcurrent) *
+                        *(ckt->CKTstate0 + here->DIOcurrent) / here->DIOtConductance;
             }
             return(OK);
         case DIO_QUEST_SENS_DC:

src/spicelib/devices/dio/dioconv.c

@@ -21,6 +21,7 @@ DIOconvTest(GENmodel *inModel, CKTcircuit *ckt)
     DIOinstance *here;
     double delvd,vd,cdhat,cd;
     double tol;
+    double delTemp, deldelTemp;
     /*  loop through all the diode models */
     for( ; model != NULL; model = DIOnextModel(model)) {
 
@@ -36,8 +37,17 @@ DIOconvTest(GENmodel *inModel, CKTcircuit *ckt)
                     *(ckt->CKTrhsOld + here->DIOnegNode);
 
             delvd=vd- *(ckt->CKTstate0 + here->DIOvoltage);
+
+            int selfheat = ((here->DIOtempNode > 0) && (here->DIOthermal) && (model->DIOrth0Given));
+            if (selfheat)
+                delTemp = *(ckt->CKTrhsOld + here->DIOtempNode);
+            else
+                delTemp = 0.0;
+            deldelTemp = delTemp - *(ckt->CKTstate0 + here->DIOdeltemp);
+
             cdhat= *(ckt->CKTstate0 + here->DIOcurrent) + 
-                    *(ckt->CKTstate0 + here->DIOconduct) * delvd;
+                    *(ckt->CKTstate0 + here->DIOconduct) * delvd +
+                    *(ckt->CKTstate0 + here->DIOdIdio_dT) * deldelTemp;
 
             cd= *(ckt->CKTstate0 + here->DIOcurrent);
 
@@ -48,7 +58,7 @@ DIOconvTest(GENmodel *inModel, CKTcircuit *ckt)
                     MAX(fabs(cdhat),fabs(cd))+ckt->CKTabstol;
             if (fabs(cdhat-cd) > tol) {
                 ckt->CKTnoncon++;
-		ckt->CKTtroubleElt = (GENinstance *) here;
+                ckt->CKTtroubleElt = (GENinstance *) here;
                 return(OK); /* don't need to check any more device */
             }
         }

src/spicelib/devices/dio/diodefs.h

@@ -38,7 +38,8 @@ typedef struct sDIOinstance {
 
     const int DIOposNode;     /* number of positive node of diode */
     const int DIOnegNode;     /* number of negative node of diode */
-    int DIOposPrimeNode;    /* number of positive prime node of diode */
+    const int DIOtempNode;    /* number of the temperature node of the diode */
+    int DIOposPrimeNode;      /* number of positive prime node of diode */
 
     double *DIOposPosPrimePtr;      /* pointer to sparse matrix at 
                                      * (positive,positive prime) */
@@ -55,6 +56,15 @@ typedef struct sDIOinstance {
     double *DIOposPrimePosPrimePtr; /* pointer to sparse matrix at 
                                      * (positive prime,positive prime) */
 
+    /* self heating */
+    double *DIOtempPosPtr;
+    double *DIOtempPosPrimePtr;
+    double *DIOtempNegPtr;
+    double *DIOtempTempPtr;
+    double *DIOposTempPtr;
+    double *DIOposPrimeTempPtr;
+    double *DIOnegTempPtr;
+
     double DIOcap;   /* stores the diode capacitance */
 
     double *DIOsens; /* stores the perturbed values of geq and ceq in ac
@@ -86,6 +96,7 @@ typedef struct sDIOinstance {
     double DIOw;        /* width for the diode */
     double DIOl;        /* length for the diode */
     double DIOm;        /* multiplier for the diode */
+    int    DIOthermal;  /* flag indicate self heating on */
 
     double DIOlengthMetal;   /* Length of metal capacitor (level=3) */
     double DIOlengthPoly;    /* Length of polysilicon capacitor (level=3) */
@@ -101,16 +112,21 @@ typedef struct sDIOinstance {
     double DIOtJctSWCap;     /* temperature adjusted sidewall junction capacitance */
     double DIOtTransitTime;  /* temperature adjusted transit time */
     double DIOtGradingCoeff; /* temperature adjusted grading coefficient (MJ) */
-    double DIOtConductance;  /* temperature adjusted series conductance */
-
-    double DIOtDepCap;  /* temperature adjusted transition point in */
-                        /* the curve matching (Fc * Vj ) */
-    double DIOtDepSWCap;  /* temperature adjusted transition point in */
-                          /* the curve matching (Fcs * Vjs ) */
-    double DIOtSatCur;  /* temperature adjusted saturation current */
-    double DIOtSatSWCur;  /* temperature adjusted side wall saturation current */
-    double DIOtTunSatCur;        /* tunneling saturation current */
-    double DIOtTunSatSWCur;      /* sidewall tunneling saturation current */
+    double DIOtConductance;    /* temperature adjusted series conductance */
+    double DIOtConductance_dT; /* temperature adjusted series conductance temperature derivative */
+
+    double DIOtDepCap;       /* temperature adjusted transition point in */
+                             /* the curve matching (Fc * Vj ) */
+    double DIOtDepSWCap;     /* temperature adjusted transition point in */
+                             /* the curve matching (Fcs * Vjs ) */
+    double DIOtSatCur;         /* temperature adjusted saturation current */
+    double DIOtSatCur_dT;      /* temperature adjusted saturation current temperature derivative */
+    double DIOtSatSWCur;       /* temperature adjusted side wall saturation current */
+    double DIOtSatSWCur_dT;    /* temperature adjusted side wall saturation current temperature derivative */
+    double DIOtTunSatCur;      /* tunneling saturation current */
+    double DIOtTunSatCur_dT;   /* tunneling saturation current temperature derivative */
+    double DIOtTunSatSWCur;    /* sidewall tunneling saturation current */
+    double DIOtTunSatSWCur_dT; /* sidewall tunneling saturation current temperature derivative */
 
     double DIOtVcrit;   /* temperature adjusted V crit */
     double DIOtF1;      /* temperature adjusted f1 */
@@ -123,9 +139,17 @@ typedef struct sDIOinstance {
 
     double DIOforwardKneeCurrent; /* Forward Knee current */
     double DIOreverseKneeCurrent; /* Reverse Knee current */
-    double DIOjunctionCap;     /* geometry adjusted junction capacitance */
-    double DIOjunctionSWCap;     /* geometry adjusted junction sidewall capacitance */
-    double DIOtRecSatCur; /* temperature adjusted recombination saturation current */
+    double DIOjunctionCap;        /* geometry adjusted junction capacitance */
+    double DIOjunctionSWCap;      /* geometry adjusted junction sidewall capacitance */
+    double DIOtRecSatCur;         /* temperature adjusted recombination saturation current */
+    double DIOtRecSatCur_dT;      /* temperature adjusted recombination saturation current */
+
+    double DIOdIth_dVrs;
+    double DIOdIth_dVdio;
+    double DIOdIth_dT;
+    double DIOgcTt;
+    double DIOdIrs_dT;
+    double DIOdIdio_dT;
 
     double DIOcmetal; /* parasitic metal overlap capacitance */
     double DIOcpoly;  /* parasitic polysilicon overlap capacitance */
@@ -175,10 +199,16 @@ typedef struct sDIOinstance {
 #define DIOcapCharge DIOstate+3
 #define DIOcapCurrent DIOstate+4
 
-#define DIOnumStates 5
+#define DIOqth DIOstate+5     /* thermal capacitor charge */
+#define DIOcqth DIOstate+6    /* thermal capacitor current */
+
+#define DIOdeltemp DIOstate+7 /* thermal voltage over rth0 */
+#define DIOdIdio_dT DIOstate+8
+
+#define DIOnumStates 9
 
-#define DIOsensxp DIOstate+5    /* charge sensitivities and their derivatives.
-                                 * +6 for the derivatives - pointer to the
+#define DIOsensxp DIOstate+9    /* charge sensitivities and their derivatives.
+                                 * +10 for the derivatives - pointer to the
                                  * beginning of the array */
 
 #define DIOnumSenStates 2
@@ -248,6 +278,9 @@ typedef struct sDIOmodel {       /* model structure for a diode */
     unsigned DIOrecSatCurGiven : 1;
     unsigned DIOrecEmissionCoeffGiven : 1;
 
+    unsigned DIOrth0Given :1;
+    unsigned DIOcth0Given :1;
+
     unsigned DIOlengthMetalGiven : 1;     /* Length of metal capacitor (level=3) */
     unsigned DIOlengthPolyGiven : 1;      /* Length of polysilicon capacitor (level=3) */
     unsigned DIOwidthMetalGiven : 1;      /* Width of metal capacitor (level=3) */
@@ -312,6 +345,9 @@ typedef struct sDIOmodel {       /* model structure for a diode */
     double DIOrecSatCur; /* Recombination saturation current */
     double DIOrecEmissionCoeff; /* Recombination emission coefficient */
 
+    double DIOrth0;
+    double DIOcth0;
+
     double DIOlengthMetal;     /* Length of metal capacitor (level=3) */
     double DIOlengthPoly;      /* Length of polysilicon capacitor (level=3) */
     double DIOwidthMetal;      /* Width of metal capacitor (level=3) */
@@ -348,6 +384,7 @@ enum {
     DIO_L,
     DIO_M,
     DIO_DTEMP,
+    DIO_THERMAL,
     DIO_LM,
     DIO_LP,
     DIO_WM,
@@ -407,6 +444,8 @@ enum {
     DIO_MOD_BV_MAX,
     DIO_MOD_ISR,
     DIO_MOD_NR,
+    DIO_MOD_RTH0,
+    DIO_MOD_CTH0,
 
     DIO_MOD_LM,
     DIO_MOD_LP,
@@ -418,5 +457,7 @@ enum {
     DIO_MOD_XP,
 };
 
+void DIOtempUpdate(DIOmodel *inModel, DIOinstance *here, double Temp, CKTcircuit *ckt);
+
 #include "dioext.h"
 #endif /*DIO*/

src/spicelib/devices/dio/dioload.c

@@ -25,7 +25,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
     double arg;
     double argsw;
     double capd;
-    double cd, cdb, cdsw;
+    double cd, cdb, cdsw, cdb_dT, cdsw_dT;
     double cdeq;
     double cdhat;
     double ceq;
@@ -46,7 +46,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
     double evd;
     double evrev;
     double gd, gdb, gdsw, gen_fac, gen_fac_vd;
-    double t1, evd_rec, cdb_rec, gdb_rec;
+    double t1, evd_rec, cdb_rec, gdb_rec, cdb_rec_dT;
     double geq;
     double gspr;    /* area-scaled conductance */
     double sarg;
@@ -56,13 +56,17 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
     double vd;      /* current diode voltage */
     double vdtemp;
     double vt;      /* K t / Q */
-    double vte, vtesw, vtetun;
-    double vtebrk;
-    int Check = 0;
+    double vte, vtesw, vtetun, vtebrk;
+    int Check_dio=0, Check_th;
     int error;
     int SenCond=0;    /* sensitivity condition */
     double diffcharge, diffchargeSW, deplcharge, deplchargeSW, diffcap, diffcapSW, deplcap, deplcapSW;
 
+    double deldelTemp, delTemp, Temp;
+    double ceqqth=0.0, Ith=0.0, gcTt=0.0, vrs=0.0;
+    double dIdio_dT, dIth_dVdio=0.0, dIrs_dT=0.0, dIth_dVrs=0.0, dIth_dT=0.0;
+    double argsw_dT, csat_dT, csatsw_dT;
+
     /*  loop through all the diode models */
     for( ; model != NULL; model = DIOnextModel(model)) {
 
@@ -70,10 +74,16 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
         for (here = DIOinstances(model); here != NULL ;
                 here=DIOnextInstance(here)) {
 
+            int selfheat = ((here->DIOtempNode > 0) && (here->DIOthermal) && (model->DIOrth0Given));
+
             /*
              *     this routine loads diodes for dc and transient analyses.
              */
 
+            if (selfheat)
+                Check_th = 1;
+            else
+                Check_th = 0;
 
             if(ckt->CKTsenInfo){
                 if((ckt->CKTsenInfo->SENstatus == PERTURBATION)
@@ -88,15 +98,15 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
             cd = 0.0;
             cdb = 0.0;
             cdsw = 0.0;
+            cdsw_dT = 0.0;
             gd = 0.0;
             gdb = 0.0;
             gdsw = 0.0;
-            csat = here->DIOtSatCur;
-            csatsw = here->DIOtSatSWCur;
-            gspr = here->DIOtConductance * here->DIOarea;
+            delTemp = 0.0;
             vt = CONSTKoverQ * here->DIOtemp;
             vte = model->DIOemissionCoeff * vt;
             vtebrk = model->DIObrkdEmissionCoeff * vt;
+            gspr = here->DIOtConductance;
             /*
              *   initialization
              */
@@ -110,8 +120,10 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
                 if((ckt->CKTsenInfo->SENmode == TRANSEN)&&
                         (ckt->CKTmode & MODEINITTRAN)) {
                     vd = *(ckt->CKTstate1 + here->DIOvoltage);
+                    delTemp = *(ckt->CKTstate1 + here->DIOdeltemp);
                 } else{
                     vd = *(ckt->CKTstate0 + here->DIOvoltage);
+                    delTemp = *(ckt->CKTstate0 + here->DIOdeltemp);
                 }
 
 #ifdef SENSDEBUG
@@ -120,20 +132,25 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
                 goto next1;
             }
 
-            Check=1;
+            Check_dio=1;
             if(ckt->CKTmode & MODEINITSMSIG) {
                 vd= *(ckt->CKTstate0 + here->DIOvoltage);
+                delTemp = *(ckt->CKTstate0 + here->DIOdeltemp);
             } else if (ckt->CKTmode & MODEINITTRAN) {
                 vd= *(ckt->CKTstate1 + here->DIOvoltage);
+                delTemp = *(ckt->CKTstate1 + here->DIOdeltemp);
             } else if ( (ckt->CKTmode & MODEINITJCT) &&
                     (ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC) ) {
                 vd=here->DIOinitCond;
             } else if ( (ckt->CKTmode & MODEINITJCT) && here->DIOoff) {
                 vd=0;
+                delTemp = 0.0;
             } else if ( ckt->CKTmode & MODEINITJCT) {
                 vd=here->DIOtVcrit;
+                delTemp = 0.0;
             } else if ( ckt->CKTmode & MODEINITFIX && here->DIOoff) {
                 vd=0;
+                delTemp = 0.0;
             } else {
 #ifndef PREDICTOR
                 if (ckt->CKTmode & MODEINITPRED) {
@@ -144,16 +161,34 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
                             *(ckt->CKTstate1 + here->DIOcurrent);
                     *(ckt->CKTstate0 + here->DIOconduct) =
                             *(ckt->CKTstate1 + here->DIOconduct);
+                    *(ckt->CKTstate0 + here->DIOdeltemp) =
+                            *(ckt->CKTstate1 + here->DIOdeltemp);
+                    delTemp = DEVpred(ckt,here->DIOdeltemp);
+                    *(ckt->CKTstate0 + here->DIOdIdio_dT) =
+                            *(ckt->CKTstate1 + here->DIOdIdio_dT);
+                    *(ckt->CKTstate0+here->DIOqth) =
+                            *(ckt->CKTstate1+here->DIOqth);
                 } else {
 #endif /* PREDICTOR */
                     vd = *(ckt->CKTrhsOld+here->DIOposPrimeNode)-
                             *(ckt->CKTrhsOld + here->DIOnegNode);
+                    if (selfheat)
+                        delTemp = *(ckt->CKTrhsOld + here->DIOtempNode);
+                    else
+                        delTemp = 0.0;
+                    *(ckt->CKTstate0+here->DIOqth) = model->DIOcth0 * delTemp;
+                    if((ckt->CKTmode & MODEINITTRAN)) {
+                        *(ckt->CKTstate1+here->DIOqth) =
+                            *(ckt->CKTstate0+here->DIOqth);
+                    }
 #ifndef PREDICTOR
                 }
 #endif /* PREDICTOR */
                 delvd=vd- *(ckt->CKTstate0 + here->DIOvoltage);
+                deldelTemp = delTemp - *(ckt->CKTstate0 + here->DIOdeltemp);
                 cdhat= *(ckt->CKTstate0 + here->DIOcurrent) +
-                        *(ckt->CKTstate0 + here->DIOconduct) * delvd;
+                        *(ckt->CKTstate0 + here->DIOconduct) * delvd +
+                        *(ckt->CKTstate0 + here->DIOdIdio_dT) * deldelTemp;
                 /*
                  *   bypass if solution has not changed
                  */
@@ -167,10 +202,17 @@ 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;
+                            if ((here->DIOtempNode == 0) ||
+                                (fabs(deldelTemp) < (ckt->CKTreltol * MAX(fabs(delTemp),
+                                      fabs(*(ckt->CKTstate0+here->DIOdeltemp)))+
+                                      ckt->CKTvoltTol*1e4))) {
+                                vd= *(ckt->CKTstate0 + here->DIOvoltage);
+                                cd= *(ckt->CKTstate0 + here->DIOcurrent);
+                                gd= *(ckt->CKTstate0 + here->DIOconduct);
+                                delTemp = *(ckt->CKTstate0 + here->DIOdeltemp);
+                                dIdio_dT= *(ckt->CKTstate0 + here->DIOdIdio_dT);
+                                goto load;
+                            }
                         }
                     }
                 }
@@ -184,17 +226,39 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
                     vdtemp = DEVpnjlim(vdtemp,
                             -(*(ckt->CKTstate0 + here->DIOvoltage) +
                             here->DIOtBrkdwnV),vtebrk,
-                            here->DIOtVcrit,&Check);
+                            here->DIOtVcrit,&Check_dio);
                     vd = -(vdtemp+here->DIOtBrkdwnV);
                 } else {
                     vd = DEVpnjlim(vd,*(ckt->CKTstate0 + here->DIOvoltage),
-                            vte,here->DIOtVcrit,&Check);
+                            vte,here->DIOtVcrit,&Check_dio);
                 }
+                if (selfheat)
+                    delTemp = DEVlimitlog(delTemp,
+                        *(ckt->CKTstate0 + here->DIOdeltemp), 100, &Check_th);
+                else
+                    delTemp = 0.0;
             }
             /*
              *   compute dc current and derivitives
              */
-next1:      if (model->DIOsatSWCurGiven) {              /* sidewall current */
+next1:      
+            if (selfheat) {
+                Temp = here->DIOtemp + delTemp;
+                DIOtempUpdate(model, here, Temp, ckt);
+            } else {
+                Temp = here->DIOtemp;
+            }
+
+            csat = here->DIOtSatCur;
+            csat_dT = here->DIOtSatCur_dT;
+            csatsw = here->DIOtSatSWCur;
+            csatsw_dT = here->DIOtSatSWCur_dT;
+            gspr = here->DIOtConductance;
+            vt = CONSTKoverQ * Temp;
+            vte = model->DIOemissionCoeff * vt;
+            vtebrk = model->DIObrkdEmissionCoeff * vt;
+
+            if (model->DIOsatSWCurGiven) {              /* sidewall current */
 
                 if (model->DIOswEmissionCoeffGiven) {   /* current with own characteristic */
 
@@ -205,61 +269,87 @@ next1:      if (model->DIOsatSWCurGiven) {              /* sidewall current */
                         evd = exp(vd/vtesw);
                         cdsw = csatsw*(evd-1);
                         gdsw = csatsw*evd/vtesw;
+                        cdsw_dT = csatsw_dT * (evd - 1) - csatsw * vd * evd / (vtesw * Temp);
 
                     } else if((!(model->DIObreakdownVoltageGiven)) ||
                             vd >= -here->DIOtBrkdwnV) { /* reverse */
 
                         argsw = 3*vtesw/(vd*CONSTe);
                         argsw = argsw * argsw * argsw;
+                        argsw_dT = 3 * argsw / Temp;
                         cdsw = -csatsw*(1+argsw);
                         gdsw = csatsw*3*argsw/vd;
+                        cdsw_dT = -csatsw_dT - (csatsw_dT*argsw + csatsw*argsw_dT);
 
                     } else {                            /* breakdown */
+                        double evrev_dT;
 
                         evrev = exp(-(here->DIOtBrkdwnV+vd)/vtebrk);
+                        evrev_dT = (here->DIOtBrkdwnV+vd)*evrev/(vtebrk*Temp);
                         cdsw = -csatsw*evrev;
                         gdsw = csatsw*evrev/vtebrk;
+                        cdsw_dT = -(csatsw_dT*evrev + csatsw*evrev_dT);
 
                     }
 
                 } else { /* merge saturation currents and use same characteristic as bottom diode */
+
                     csat = csat + csatsw;
+                    csat_dT = csat_dT + csatsw_dT;
+                    cdsw_dT = 0.0;
 
                 }
 
             }
 
+            /*
+             *   temperature dependent diode saturation current and derivative
+             */
+
             if (vd >= -3*vte) {                 /* bottom current forward */
 
                 evd = exp(vd/vte);
                 cdb = csat*(evd-1);
                 gdb = csat*evd/vte;
+                cdb_dT = csat_dT * (evd - 1) - csat * vd * evd / (vte * Temp);
                 if (model->DIOrecSatCurGiven) { /* recombination current */
-                    evd_rec = exp(vd/(model->DIOrecEmissionCoeff*vt));
+                    double vterec = model->DIOrecEmissionCoeff*vt;
+                    evd_rec = exp(vd/(vterec));
                     cdb_rec = here->DIOtRecSatCur*(evd_rec-1);
-                    gdb_rec = here->DIOtRecSatCur*evd_rec/vt;
+                    gdb_rec = here->DIOtRecSatCur*evd_rec/vterec;
+                    cdb_rec_dT = here->DIOtRecSatCur_dT * (evd_rec - 1) 
+                                -here->DIOtRecSatCur * vd * evd_rec / (vterec*Temp);
                     t1 = pow((1-vd/here->DIOtJctPot), 2) + 0.005;
                     gen_fac = pow(t1, here->DIOtGradingCoeff/2);
-                    gen_fac_vd = here->DIOtGradingCoeff * (1-vd/here->DIOtJctPot) * pow(t1, (here->DIOtGradingCoeff/2-1));
+                    gen_fac_vd = -here->DIOtGradingCoeff * (1-vd/here->DIOtJctPot) 
+                                                         * pow(t1, (here->DIOtGradingCoeff/2-1));
                     cdb_rec = cdb_rec * gen_fac;
                     gdb_rec = gdb_rec * gen_fac + cdb_rec * gen_fac_vd;
                     cdb = cdb + cdb_rec;
                     gdb = gdb + gdb_rec;
+                    cdb_dT = cdb_dT + cdb_rec_dT*gen_fac;
                 }
 
             } else if((!(model->DIObreakdownVoltageGiven)) ||
                     vd >= -here->DIOtBrkdwnV) { /* reverse */
 
+                double darg_dT;
+
                 arg = 3*vte/(vd*CONSTe);
                 arg = arg * arg * arg;
+                darg_dT = 3 * arg / Temp; 
                 cdb = -csat*(1+arg);
                 gdb = csat*3*arg/vd;
+                cdb_dT = -csat_dT - (csat_dT*arg + csat*darg_dT);
 
             } else {                            /* breakdown */
+                double evrev_dT;
 
                 evrev = exp(-(here->DIOtBrkdwnV+vd)/vtebrk);
+                evrev_dT = (here->DIOtBrkdwnV+vd)*evrev/(vtebrk*Temp);
                 cdb = -csat*evrev;
                 gdb = csat*evrev/vtebrk;
+                cdb_dT = -(csat_dT*evrev + csat*evrev_dT);
 
             }
 
@@ -270,6 +360,8 @@ next1:      if (model->DIOsatSWCurGiven) {              /* sidewall current */
 
                 cdsw = cdsw - here->DIOtTunSatSWCur * (evd - 1);
                 gdsw = gdsw + here->DIOtTunSatSWCur * evd / vtetun;
+                cdsw_dT = cdsw_dT - here->DIOtTunSatSWCur_dT * (evd - 1) 
+                                  - here->DIOtTunSatSWCur * vd * evd / (vtetun * Temp);
 
             }
 
@@ -280,11 +372,14 @@ next1:      if (model->DIOsatSWCurGiven) {              /* sidewall current */
 
                 cdb = cdb - here->DIOtTunSatCur * (evd - 1);
                 gdb = gdb + here->DIOtTunSatCur * evd / vtetun;
+                cdb_dT = cdb_dT - here->DIOtTunSatCur_dT * (evd - 1)
+                                - here->DIOtTunSatCur * vd * evd / (vtetun * Temp);
 
             }
 
             cd = cdb + cdsw;
             gd = gdb + gdsw;
+            dIdio_dT = cdb_dT + cdsw_dT;
 
             if (vd >= -3*vte) { /* limit forward */
 
@@ -365,6 +460,7 @@ next1:      if (model->DIOsatSWCurGiven) {              /* sidewall current */
                         if(SenCond){
                             *(ckt->CKTstate0 + here->DIOcurrent) = cd;
                             *(ckt->CKTstate0 + here->DIOconduct) = gd;
+                            *(ckt->CKTstate0 + here->DIOdIdio_dT) = dIdio_dT;
 #ifdef SENSDEBUG
                             printf("storing small signal parameters\n");
                             printf("cd = %.7e,vd = %.7e\n",cd,vd);
@@ -400,6 +496,15 @@ next1:      if (model->DIOsatSWCurGiven) {              /* sidewall current */
                         *(ckt->CKTstate1 + here->DIOcapCurrent) =
                                 *(ckt->CKTstate0 + here->DIOcapCurrent);
                     }
+                    if (selfheat)
+                    {
+                        error = NIintegrate(ckt, &gcTt, &ceqqth, model->DIOcth0, here->DIOqth);
+                        if (error) return(error);
+                        if (ckt->CKTmode & MODEINITTRAN) {
+                            *(ckt->CKTstate1 + here->DIOcqth) =
+                                    *(ckt->CKTstate0 + here->DIOcqth);
+                        }
+                    }
                 }
             }
 
@@ -409,7 +514,7 @@ next1:      if (model->DIOsatSWCurGiven) {              /* sidewall current */
              *   check convergence
              */
             if ( (!(ckt->CKTmode & MODEINITFIX)) || (!(here->DIOoff))  ) {
-                if (Check == 1)  {
+                if ((Check_th == 1) || (Check_dio == 1)) {
                     ckt->CKTnoncon++;
                     ckt->CKTtroubleElt = (GENinstance *) here;
                 }
@@ -417,18 +522,44 @@ next1:      if (model->DIOsatSWCurGiven) {              /* sidewall current */
 next2:      *(ckt->CKTstate0 + here->DIOvoltage) = vd;
             *(ckt->CKTstate0 + here->DIOcurrent) = cd;
             *(ckt->CKTstate0 + here->DIOconduct) = gd;
+            *(ckt->CKTstate0 + here->DIOdeltemp) = delTemp;
+            *(ckt->CKTstate0 + here->DIOdIdio_dT) = dIdio_dT;
 
             if(SenCond)  continue;
 
 #ifndef NOBYPASS
             load:
 #endif
+            if (selfheat) {
+                double dIrs_dVrs, dIrs_dgspr, dIth_dIrs;
+                vrs = *(ckt->CKTrhsOld + here->DIOposNode) - *(ckt->CKTrhsOld + here->DIOposPrimeNode);
+                Ith = vd*cd + vrs*vrs*gspr; /* Diode dissipated power */
+                dIrs_dVrs = gspr;
+                dIrs_dgspr = vrs;
+                dIrs_dT = dIrs_dgspr * here->DIOtConductance_dT;
+                dIth_dVrs = vrs*gspr;
+                dIth_dIrs = vrs;
+                dIth_dVrs = dIth_dVrs + dIth_dIrs*dIrs_dVrs;
+                dIth_dT = dIth_dIrs*dIrs_dT + dIdio_dT*vd;
+                dIth_dVdio = cd + vd*gd;
+                here->DIOdIth_dVrs = dIth_dVrs;
+                here->DIOdIth_dVdio = dIth_dVdio;
+                here->DIOdIth_dT = dIth_dT;
+                here->DIOgcTt = gcTt;
+                here->DIOdIrs_dT = dIrs_dT;
+            }
             /*
              *   load current vector
              */
             cdeq=cd-gd*vd;
             *(ckt->CKTrhs + here->DIOnegNode) += cdeq;
             *(ckt->CKTrhs + here->DIOposPrimeNode) -= cdeq;
+            if (selfheat) {
+                *(ckt->CKTrhs + here->DIOposNode)      +=  dIrs_dT*delTemp;
+                *(ckt->CKTrhs + here->DIOposPrimeNode) +=  dIdio_dT*delTemp - dIrs_dT*delTemp;
+                *(ckt->CKTrhs + here->DIOnegNode)      += -dIdio_dT*delTemp;
+                *(ckt->CKTrhs + here->DIOtempNode)     +=  Ith - dIth_dVdio*vd - dIth_dVrs*vrs - dIth_dT*delTemp - ceqqth;
+            }
             /*
              *   load matrix
              */
@@ -439,6 +570,15 @@ next2:      *(ckt->CKTstate0 + here->DIOvoltage) = vd;
             *(here->DIOnegPosPrimePtr) -= gd;
             *(here->DIOposPrimePosPtr) -= gspr;
             *(here->DIOposPrimeNegPtr) -= gd;
+            if (selfheat) {
+                (*(here->DIOtempPosPtr)      += -dIth_dVrs);
+                (*(here->DIOtempPosPrimePtr) += -dIth_dVdio + dIth_dVrs);
+                (*(here->DIOtempNegPtr)      +=  dIth_dVdio);
+                (*(here->DIOtempTempPtr)     += -dIth_dT + 1/model->DIOrth0 + gcTt);
+                (*(here->DIOposTempPtr)      +=  dIrs_dT);
+                (*(here->DIOposPrimeTempPtr) +=  dIdio_dT - dIrs_dT);
+                (*(here->DIOnegTempPtr)      += -dIdio_dT);
+            }
         }
     }
     return(OK);

src/spicelib/devices/dio/diomask.c

@@ -178,6 +178,12 @@ DIOmAsk (CKTcircuit *ckt, GENmodel *inModel, int which, IFvalue *value)
         case DIO_MOD_NR:
             value->rValue = model->DIOrecEmissionCoeff;
             return(OK);
+        case DIO_MOD_RTH0:
+            value->rValue = model->DIOrth0; 
+            return(OK);
+        case DIO_MOD_CTH0:
+            value->rValue = model->DIOcth0; 
+            return(OK);
 
         case DIO_MOD_LM:
             value->rValue = model->DIOlengthMetal;

src/spicelib/devices/dio/diompar.c

@@ -217,6 +217,14 @@ DIOmParam(int param, IFvalue *value, GENmodel *inModel)
             model->DIOrecEmissionCoeff = value->rValue;
             model->DIOrecEmissionCoeffGiven = TRUE;
             break;
+        case  DIO_MOD_RTH0:
+            model->DIOrth0 = value->rValue;
+            model->DIOrth0Given = TRUE;
+            break;
+        case  DIO_MOD_CTH0:
+            model->DIOcth0 = value->rValue;
+            model->DIOcth0Given = TRUE;
+            break;
 
         case DIO_MOD_LM:
             model->DIOlengthMetal = value->rValue;

src/spicelib/devices/dio/dioparam.c

@@ -44,6 +44,9 @@ DIOparam(int param, IFvalue *value, GENinstance *inst, IFvalue *select)
             here->DIOm = value->rValue;
             here->DIOmGiven = TRUE;
             break;
+        case DIO_THERMAL:
+            here->DIOthermal = (value->iValue != 0);
+            break;
 
         case DIO_TEMP:
             here->DIOtemp = value->rValue+CONSTCtoK;

src/spicelib/devices/dio/diopzld.c

@@ -29,7 +29,7 @@ DIOpzLoad(GENmodel *inModel, CKTcircuit *ckt, SPcomplex *s)
         /* loop through all the instances of the model */
         for (here = DIOinstances(model); here != NULL ;
                 here=DIOnextInstance(here)) {
-            gspr=here->DIOtConductance*here->DIOarea;
+            gspr=here->DIOtConductance;
             geq= *(ckt->CKTstate0 + here->DIOconduct);
             xceq= *(ckt->CKTstate0 + here->DIOcapCurrent);
             *(here->DIOposPosPtr ) += gspr;

src/spicelib/devices/dio/diosacl.c

@@ -110,7 +110,7 @@ DIOsAcLoad(GENmodel *inModel, CKTcircuit *ckt)
             geq0 = *(here->DIOsenGeq);
             xceq0 = *(here->DIOsenCeq) * ckt->CKTomega;
             A0 = here->DIOarea;
-            gspr0=here->DIOtConductance*A0;
+            gspr0=here->DIOtConductance;
             cpos0 = gspr0 * vspr;
             icpos0 = gspr0 * ivspr;
             cposprm0 = geq0 * vd - xceq0 * ivd - cpos0;

src/spicelib/devices/dio/diosetup.c

@@ -171,6 +171,27 @@ DIOsetup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *states)
             model->DIOrecSatCur = 1e-14;
         }
 
+        /* set lower limit of saturation current */
+        if (model->DIOsatCur < ckt->CKTepsmin)
+            model->DIOsatCur = ckt->CKTepsmin;
+
+        if(!model->DIOnomTempGiven) {
+            model->DIOnomTemp = ckt->CKTnomTemp;
+        }
+
+        if((!model->DIOresistGiven) || (model->DIOresist==0)) {
+            model->DIOconductance = 0.0;
+        } else {
+            model->DIOconductance = 1/model->DIOresist;
+        }
+
+        if (!model->DIOrth0Given) {
+            model->DIOrth0 = 0;
+        }
+        if (!model->DIOcth0Given) {
+            model->DIOcth0 = 1e-5;
+        }
+
         if(!model->DIOlengthMetalGiven) {
             model->DIOlengthMetal = 0.0;
         }
@@ -286,6 +307,8 @@ DIOsetup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *states)
                 }
             }
 
+            int selfheat = ((here->DIOtempNode > 0) && (here->DIOthermal) && (model->DIOrth0Given));
+
 /* macro to make elements with built in test for out of memory */
 #define TSTALLOC(ptr,first,second) \
 do { if((here->ptr = SMPmakeElt(matrix, here->first, here->second)) == NULL){\
@@ -299,6 +322,17 @@ do { if((here->ptr = SMPmakeElt(matrix, here->first, here->second)) == NULL){\
             TSTALLOC(DIOposPosPtr,DIOposNode,DIOposNode);
             TSTALLOC(DIOnegNegPtr,DIOnegNode,DIOnegNode);
             TSTALLOC(DIOposPrimePosPrimePtr,DIOposPrimeNode,DIOposPrimeNode);
+
+            if (selfheat) {
+                TSTALLOC(DIOtempPosPtr,      DIOtempNode,     DIOposNode);
+                TSTALLOC(DIOtempPosPrimePtr, DIOtempNode,     DIOposPrimeNode);
+                TSTALLOC(DIOtempNegPtr,      DIOtempNode,     DIOnegNode);
+                TSTALLOC(DIOtempTempPtr,     DIOtempNode,     DIOtempNode);
+                TSTALLOC(DIOposTempPtr,      DIOposNode,      DIOtempNode);
+                TSTALLOC(DIOposPrimeTempPtr, DIOposPrimeNode, DIOtempNode);
+                TSTALLOC(DIOnegTempPtr,      DIOnegNode,      DIOtempNode);
+            }
+
         }
     }
     return(OK);

src/spicelib/devices/dio/diotemp.c

@@ -5,7 +5,7 @@ Modified: 2000 AlansFixes
 Modified by Paolo Nenzi 2003 and Dietmar Warning 2012
 **********/
 
-    /* perform the temperature update to the diode */
+/* perform the temperature update to the diode */
 
 #include "ngspice/ngspice.h"
 #include "ngspice/cktdefs.h"
@@ -15,269 +15,256 @@ Modified by Paolo Nenzi 2003 and Dietmar Warning 2012
 #include "ngspice/suffix.h"
 #include "ngspice/cpdefs.h"
 
-int
-DIOtemp(GENmodel *inModel, CKTcircuit *ckt)
-{
+void DIOtempUpdate(DIOmodel *inModel, DIOinstance *here, double Temp, CKTcircuit *ckt) {
+
     DIOmodel *model = (DIOmodel*)inModel;
     double xfc, xfcs;
-    double vte;
+    double vt, vte, vts, vtt, vtr;
     double cbv;
     double xbv;
     double xcbv;
     double tol;
-    double vt;
     double vtnom;
-    DIOinstance *here;
     int iter;
     double dt;
     double factor;
     double tBreakdownVoltage;
 
+    double egfet1,arg1,fact1,pbfact1,pbo,gmaold,pboSW,gmaSWold;
+    double fact2,pbfact,arg,egfet,gmanew,gmaSWnew;
+    double arg1_dT, arg2, arg2_dT;
     double gclimit;
 
     if (!cp_getvar("DIOgradingCoeffMax", CP_REAL, &gclimit, 0))
         gclimit = 0.9;
 
-    /*  loop through all the diode models */
-    for( ; model != NULL; model = DIOnextModel(model)) {
-        if(!model->DIOnomTempGiven) {
-            model->DIOnomTemp = ckt->CKTnomTemp;
-        }
-        vtnom = CONSTKoverQ * model->DIOnomTemp;
-        /* limit grading coeff to max of .9, set new limit with variable DIOgradingCoeffMax */
-        if(model->DIOgradingCoeff>gclimit) {
-            SPfrontEnd->IFerrorf (ERR_WARNING,
-                    "%s: grading coefficient too large, limited to %g",
-                    model->DIOmodName, gclimit);
-            model->DIOgradingCoeff=gclimit;
-        }
-        /* limit activation energy to min of .1 */
-        if(model->DIOactivationEnergy<.1) {
-            SPfrontEnd->IFerrorf (ERR_WARNING,
-                    "%s: activation energy too small, limited to 0.1",
-                    model->DIOmodName);
-            model->DIOactivationEnergy=.1;
+    vt = CONSTKoverQ * Temp;
+    vte = model->DIOemissionCoeff * vt;
+    vts = model->DIOswEmissionCoeff * vt;
+    vtt = model->DIOtunEmissionCoeff * vt;
+    vtr = model->DIOrecEmissionCoeff * vt;
+    vtnom = CONSTKoverQ * model->DIOnomTemp;
+    dt = Temp - model->DIOnomTemp;
+
+    /* Junction grading temperature adjust */
+    factor = 1.0 + (model->DIOgradCoeffTemp1 * dt)
+                 + (model->DIOgradCoeffTemp2 * dt * dt);
+    here->DIOtGradingCoeff = model->DIOgradingCoeff * factor;
+
+    /* limit temperature adjusted grading coeff
+     * to max of .9, or set new limit with variable DIOgradingCoeffMax
+     */
+    if(here->DIOtGradingCoeff>gclimit) {
+      SPfrontEnd->IFerrorf (ERR_WARNING,
+            "%s: temperature adjusted grading coefficient too large, limited to %g",
+            here->DIOname, gclimit);
+      here->DIOtGradingCoeff=gclimit;
+    }
+
+    /* this part gets really ugly - I won't even try to
+     * explain these equations */
+    fact2 = Temp/REFTEMP;
+    egfet = 1.16-(7.02e-4*Temp*Temp)/
+            (Temp+1108);
+    arg = -egfet/(2*CONSTboltz*Temp) +
+            1.1150877/(CONSTboltz*(REFTEMP+REFTEMP));
+    pbfact = -2*vt*(1.5*log(fact2)+CHARGE*arg);
+    egfet1 = 1.16 - (7.02e-4*model->DIOnomTemp*model->DIOnomTemp)/
+            (model->DIOnomTemp+1108);
+    arg1 = -egfet1/(CONSTboltz*2*model->DIOnomTemp) +
+            1.1150877/(2*CONSTboltz*REFTEMP);
+    fact1 = model->DIOnomTemp/REFTEMP;
+    pbfact1 = -2 * vtnom*(1.5*log(fact1)+CHARGE*arg1);
+
+    if (model->DIOtlevc == 0) {
+            pbo = (model->DIOjunctionPot-pbfact1)/fact1;
+            gmaold = (model->DIOjunctionPot-pbo)/pbo;
+            here->DIOtJctCap = here->DIOjunctionCap /
+                    (1+here->DIOtGradingCoeff*
+                    (400e-6*(model->DIOnomTemp-REFTEMP)-gmaold) );
+            here->DIOtJctPot = pbfact+fact2*pbo;
+            gmanew = (here->DIOtJctPot-pbo)/pbo;
+            here->DIOtJctCap *= 1+here->DIOtGradingCoeff*
+                    (400e-6*(Temp-REFTEMP)-gmanew);
+    } else if (model->DIOtlevc == 1) {
+            here->DIOtJctPot = model->DIOjunctionPot - model->DIOtpb*(Temp-REFTEMP);
+            here->DIOtJctCap = here->DIOjunctionCap *
+                    (1+model->DIOcta*(Temp-REFTEMP));
+    }
+
+    if (model->DIOtlevc == 0) {
+            pboSW = (model->DIOjunctionSWPot-pbfact1)/fact1;
+            gmaSWold = (model->DIOjunctionSWPot-pboSW)/pboSW;
+            here->DIOtJctSWCap = here->DIOjunctionSWCap /
+                    (1+model->DIOgradingSWCoeff*
+                    (400e-6*(model->DIOnomTemp-REFTEMP)-gmaSWold) );
+            here->DIOtJctSWPot = pbfact+fact2*pboSW;
+            gmaSWnew = (here->DIOtJctSWPot-pboSW)/pboSW;
+            here->DIOtJctSWCap *= 1+model->DIOgradingSWCoeff*
+                    (400e-6*(Temp-REFTEMP)-gmaSWnew);
+    } else if (model->DIOtlevc == 1) {
+            here->DIOtJctSWPot = model->DIOjunctionSWPot - model->DIOtphp*(Temp-REFTEMP);
+            here->DIOtJctSWCap = here->DIOjunctionSWCap *
+                    (1+model->DIOctp*(Temp-REFTEMP));
+    }
+
+    arg1 = ((Temp / model->DIOnomTemp) - 1) * model->DIOactivationEnergy / vte;
+    arg1_dT = model->DIOactivationEnergy / (vte*model->DIOnomTemp)
+              - model->DIOactivationEnergy*(Temp/model->DIOnomTemp -1)/(vte*Temp);
+    arg2 = model->DIOsaturationCurrentExp / model->DIOemissionCoeff * log(Temp / model->DIOnomTemp);
+    arg2_dT = model->DIOsaturationCurrentExp / model->DIOemissionCoeff / Temp;
+    here->DIOtSatCur = here->DIOm * model->DIOsatCur * here->DIOarea * exp(arg1 + arg2);
+    here->DIOtSatCur_dT = here->DIOm * model->DIOsatCur * here->DIOarea * exp(arg1 + arg2) * (arg1_dT + arg2_dT);
+
+    arg1 = ((Temp / model->DIOnomTemp) - 1) * model->DIOactivationEnergy / vts;
+    arg1_dT = model->DIOactivationEnergy / (vts*model->DIOnomTemp)
+              - model->DIOactivationEnergy*(Temp/model->DIOnomTemp -1)/(vts*Temp);
+    arg2 = model->DIOsaturationCurrentExp / model->DIOswEmissionCoeff * log(Temp / model->DIOnomTemp);
+    arg2_dT = model->DIOsaturationCurrentExp / model->DIOswEmissionCoeff / Temp;
+    here->DIOtSatSWCur = here->DIOm * model->DIOsatSWCur * here->DIOpj * exp(arg1 + arg2);
+    here->DIOtSatSWCur_dT = here->DIOm * model->DIOsatSWCur * here->DIOpj * exp(arg1 + arg2) * (arg1_dT + arg2_dT);
+
+    arg1 = ((Temp / model->DIOnomTemp) - 1) * model->DIOtunEGcorrectionFactor*model->DIOactivationEnergy / vtt;
+    arg1_dT = model->DIOtunEGcorrectionFactor*model->DIOactivationEnergy / (vtt*model->DIOnomTemp)
+              - model->DIOactivationEnergy*(Temp/model->DIOnomTemp -1)/(vtt*Temp);
+    arg2 = model->DIOtunSaturationCurrentExp / model->DIOtunEmissionCoeff * log(Temp / model->DIOnomTemp);
+    arg2_dT = model->DIOtunSaturationCurrentExp / model->DIOtunEmissionCoeff / Temp;
+    here->DIOtTunSatCur = here->DIOm * model->DIOtunSatCur * here->DIOarea * exp(arg1 + arg2);
+    here->DIOtTunSatCur_dT = here->DIOm * model->DIOtunSatCur * here->DIOarea * exp(arg1 + arg2) * (arg1_dT + arg2_dT);
+
+    arg1 = ((Temp / model->DIOnomTemp) - 1) * model->DIOtunEGcorrectionFactor*model->DIOactivationEnergy / vtt;
+    arg1_dT = model->DIOtunEGcorrectionFactor*model->DIOactivationEnergy / (vtt*model->DIOnomTemp)
+              - model->DIOactivationEnergy*(Temp/model->DIOnomTemp -1)/(vtt*Temp);
+    arg2 = model->DIOtunSaturationCurrentExp / model->DIOtunEmissionCoeff * log(Temp / model->DIOnomTemp);
+    arg2_dT = model->DIOtunSaturationCurrentExp / model->DIOtunEmissionCoeff / Temp;
+    here->DIOtTunSatSWCur = here->DIOm * model->DIOtunSatSWCur * here->DIOpj * exp(arg1 + arg2);
+    here->DIOtTunSatSWCur_dT = here->DIOm * model->DIOtunSatSWCur * here->DIOpj * exp(arg1 + arg2) * (arg1_dT + arg2_dT);
+
+    arg1 = ((Temp / model->DIOnomTemp) - 1) * model->DIOactivationEnergy / vtr;
+    arg1_dT = model->DIOactivationEnergy / (vtr*model->DIOnomTemp)
+              - model->DIOactivationEnergy*(Temp/model->DIOnomTemp -1)/(vtr*Temp);
+    arg2 = model->DIOsaturationCurrentExp / model->DIOrecEmissionCoeff * log(Temp / model->DIOnomTemp);
+    arg2_dT = model->DIOsaturationCurrentExp / model->DIOrecEmissionCoeff / Temp;
+    here->DIOtRecSatCur = here->DIOm * model->DIOrecSatCur * here->DIOarea * exp(arg1 + arg2);
+    here->DIOtRecSatCur_dT = here->DIOm * model->DIOrecSatCur * here->DIOarea * exp(arg1 + arg2) * (arg1_dT + arg2_dT);
+
+    xfc=log(1-model->DIOdepletionCapCoeff);
+    xfcs=log(1-model->DIOdepletionSWcapCoeff);
+
+    /* the defintion of f1, just recompute after temperature adjusting
+     * all the variables used in it */
+    here->DIOtF1=here->DIOtJctPot*
+            (1-exp((1-here->DIOtGradingCoeff)*xfc))/
+            (1-here->DIOtGradingCoeff);
+    /* same for Depletion Capacitance */
+    here->DIOtDepCap=model->DIOdepletionCapCoeff*
+            here->DIOtJctPot;
+    here->DIOtDepSWCap=model->DIOdepletionSWcapCoeff*
+            here->DIOtJctSWPot;
+    /* and Vcrit */
+    vte=model->DIOemissionCoeff*vt;
+
+    here->DIOtVcrit = vte * log(vte/(CONSTroot2*here->DIOtSatCur));
+
+    /* limit junction potential to max of 1/FC */
+    if(here->DIOtDepCap > 1.0) {
+        here->DIOtJctPot=1.0/model->DIOdepletionCapCoeff;
+        here->DIOtDepCap=model->DIOdepletionCapCoeff*here->DIOtJctPot;
+        SPfrontEnd->IFerrorf (ERR_WARNING,
+                "%s: junction potential VJ too large, limited to %f",
+                model->DIOmodName, here->DIOtJctPot);
+    }
+    /* limit sidewall junction potential to max of 1/FCS */
+    if(here->DIOtDepSWCap > 1.0) {
+        here->DIOtJctSWPot=1.0/model->DIOdepletionSWcapCoeff;
+        here->DIOtDepSWCap=model->DIOdepletionSWcapCoeff*here->DIOtJctSWPot;
+        SPfrontEnd->IFerrorf (ERR_WARNING,
+                "%s: junction potential VJS too large, limited to %f",
+                model->DIOmodName, here->DIOtJctSWPot);
+    }
+
+    /* and now to compute the breakdown voltage, again, using
+     * temperature adjusted basic parameters */
+    if (model->DIObreakdownVoltageGiven){
+        if (model->DIOtlev == 0) {
+            tBreakdownVoltage = model->DIObreakdownVoltage - model->DIOtcv * dt;
+        } else {
+            tBreakdownVoltage = model->DIObreakdownVoltage * (1 - model->DIOtcv * dt);
         }
-        /* limit depletion cap coeff to max of .95 */
-        if(model->DIOdepletionCapCoeff>.95) {
-            SPfrontEnd->IFerrorf (ERR_WARNING,
-                    "%s: coefficient Fc too large, limited to 0.95",
-                    model->DIOmodName);
-            model->DIOdepletionCapCoeff=.95;
+        if (model->DIOlevel == 1) {
+            cbv = here->DIOm * model->DIObreakdownCurrent;
+        } else { /* level=3 */
+            cbv = here->DIOm * model->DIObreakdownCurrent * here->DIOarea;
         }
-        /* limit sidewall depletion cap coeff to max of .95 */
-        if(model->DIOdepletionSWcapCoeff>.95) {
-            SPfrontEnd->IFerrorf (ERR_WARNING,
-                    "%s: coefficient Fcs too large, limited to 0.95",
-                    model->DIOmodName);
-            model->DIOdepletionSWcapCoeff=.95;
+        if (cbv < here->DIOtSatCur * tBreakdownVoltage/vt) {
+            cbv=here->DIOtSatCur * tBreakdownVoltage/vt;
+#ifdef TRACE
+        SPfrontEnd->IFerrorf (ERR_WARNING, "%s: breakdown current increased to %g to resolve", here->DIOname, cbv);
+        SPfrontEnd->IFerrorf (ERR_WARNING,
+        "incompatibility with specified saturation current");
+#endif
+        xbv=tBreakdownVoltage;
+    } else {
+        tol=ckt->CKTreltol*cbv;
+        xbv=tBreakdownVoltage-model->DIObrkdEmissionCoeff*vt*log(1+cbv/
+                (here->DIOtSatCur));
+        iter=0;
+        for(iter=0 ; iter < 25 ; iter++) {
+            xbv=tBreakdownVoltage-model->DIObrkdEmissionCoeff*vt*log(cbv/
+                    (here->DIOtSatCur)+1-xbv/vt);
+            xcbv=here->DIOtSatCur *
+                 (exp((tBreakdownVoltage-xbv)/(model->DIObrkdEmissionCoeff*vt))-1+xbv/vt);
+            if (fabs(xcbv-cbv) <= tol) goto matched;
         }
-        /* set lower limit of saturation current */
-        if (model->DIOsatCur < ckt->CKTepsmin)
-            model->DIOsatCur = ckt->CKTepsmin;
-        if((!model->DIOresistGiven) || (model->DIOresist==0)) {
-            model->DIOconductance = 0.0;
-        } else {
-            model->DIOconductance = 1/model->DIOresist;
+#ifdef TRACE
+                    SPfrontEnd->IFerrorf (ERR_WARNING, "%s: unable to match forward and reverse diode regions: bv = %g, ibv = %g", here->DIOname, xbv, xcbv);
+#endif
         }
-        xfc=log(1-model->DIOdepletionCapCoeff);
-        xfcs=log(1-model->DIOdepletionSWcapCoeff);
+        matched:
+            here->DIOtBrkdwnV = xbv;
+    }
+
+    /* transit time temperature adjust */
+    factor = 1.0 + (model->DIOtranTimeTemp1 * dt)
+                 + (model->DIOtranTimeTemp2 * dt * dt);
+    here->DIOtTransitTime = model->DIOtransitTime * factor;
+
+    /* Series resistance temperature adjust */
+    here->DIOtConductance = here->DIOm * model->DIOconductance * here->DIOarea;
+    if(model->DIOresistGiven && model->DIOresist!=0.0) {
+        factor = 1.0 + (model->DIOresistTemp1) * dt
+                 + (model->DIOresistTemp2 * dt * dt);
+        here->DIOtConductance = here->DIOm * model->DIOconductance * here->DIOarea / factor;
+        here->DIOtConductance_dT = here->DIOm * -model->DIOconductance * here->DIOarea *
+            (model->DIOresistTemp1 + model->DIOresistTemp2 * dt) / (factor*factor);
+    }
+
+    here->DIOtF2=exp((1+here->DIOtGradingCoeff)*xfc);
+    here->DIOtF3=1-model->DIOdepletionCapCoeff*
+            (1+here->DIOtGradingCoeff);
+    here->DIOtF2SW=exp((1+model->DIOgradingSWCoeff)*xfcs);
+    here->DIOtF3SW=1-model->DIOdepletionSWcapCoeff*
+            (1+model->DIOgradingSWCoeff);
+}
+
+int
+DIOtemp(GENmodel *inModel, CKTcircuit *ckt)
+{
+    DIOmodel *model = (DIOmodel*)inModel;
+    DIOinstance *here;
+
+    /*  loop through all the diode models */
+    for( ; model != NULL; model = DIOnextModel(model)) {
 
+        /* loop through all the instances */
         for(here=DIOinstances(model);here;here=DIOnextInstance(here)) {
-            double egfet1,arg1,fact1,pbfact1,pbo,gmaold,pboSW,gmaSWold;
-            double fact2,pbfact,arg,egfet,gmanew,gmaSWnew;
-            /* loop through all the instances */
 
             if(!here->DIOdtempGiven) here->DIOdtemp = 0.0;
 
             if(!here->DIOtempGiven)
                 here->DIOtemp = ckt->CKTtemp + here->DIOdtemp;
 
-            dt = here->DIOtemp - model->DIOnomTemp;
-
-            /* Junction grading temperature adjust */
-            factor = 1.0 + (model->DIOgradCoeffTemp1 * dt)
-                         + (model->DIOgradCoeffTemp2 * dt * dt);
-            here->DIOtGradingCoeff = model->DIOgradingCoeff * factor;
-
-            /* limit temperature adjusted grading coeff
-             * to max of .9, or set new limit with variable DIOgradingCoeffMax
-             */
-            if(here->DIOtGradingCoeff>gclimit) {
-              SPfrontEnd->IFerrorf (ERR_WARNING,
-                    "%s: temperature adjusted grading coefficient too large, limited to %g",
-                    here->DIOname, gclimit);
-              here->DIOtGradingCoeff=gclimit;
-            }
-
-            vt = CONSTKoverQ * here->DIOtemp;
-            /* this part gets really ugly - I won't even try to
-             * explain these equations */
-            fact2 = here->DIOtemp/REFTEMP;
-            egfet = 1.16-(7.02e-4*here->DIOtemp*here->DIOtemp)/
-                    (here->DIOtemp+1108);
-            arg = -egfet/(2*CONSTboltz*here->DIOtemp) +
-                    1.1150877/(CONSTboltz*(REFTEMP+REFTEMP));
-            pbfact = -2*vt*(1.5*log(fact2)+CHARGE*arg);
-            egfet1 = 1.16 - (7.02e-4*model->DIOnomTemp*model->DIOnomTemp)/
-                    (model->DIOnomTemp+1108);
-            arg1 = -egfet1/(CONSTboltz*2*model->DIOnomTemp) +
-                    1.1150877/(2*CONSTboltz*REFTEMP);
-            fact1 = model->DIOnomTemp/REFTEMP;
-            pbfact1 = -2 * vtnom*(1.5*log(fact1)+CHARGE*arg1);
-
-            if (model->DIOtlevc == 0) {
-                    pbo = (model->DIOjunctionPot-pbfact1)/fact1;
-                    gmaold = (model->DIOjunctionPot-pbo)/pbo;
-                    here->DIOtJctCap = here->DIOjunctionCap /
-                            (1+here->DIOtGradingCoeff*
-                            (400e-6*(model->DIOnomTemp-REFTEMP)-gmaold) );
-                    here->DIOtJctPot = pbfact+fact2*pbo;
-                    gmanew = (here->DIOtJctPot-pbo)/pbo;
-                    here->DIOtJctCap *= 1+here->DIOtGradingCoeff*
-                            (400e-6*(here->DIOtemp-REFTEMP)-gmanew);
-            } else if (model->DIOtlevc == 1) {
-                    here->DIOtJctPot = model->DIOjunctionPot - model->DIOtpb*(here->DIOtemp-REFTEMP);
-                    here->DIOtJctCap = here->DIOjunctionCap *
-                            (1+model->DIOcta*(here->DIOtemp-REFTEMP));
-            }
-
-            if (model->DIOtlevc == 0) {
-                    pboSW = (model->DIOjunctionSWPot-pbfact1)/fact1;
-                    gmaSWold = (model->DIOjunctionSWPot-pboSW)/pboSW;
-                    here->DIOtJctSWCap = here->DIOjunctionSWCap /
-                            (1+model->DIOgradingSWCoeff*
-                            (400e-6*(model->DIOnomTemp-REFTEMP)-gmaSWold) );
-                    here->DIOtJctSWPot = pbfact+fact2*pboSW;
-                    gmaSWnew = (here->DIOtJctSWPot-pboSW)/pboSW;
-                    here->DIOtJctSWCap *= 1+model->DIOgradingSWCoeff*
-                            (400e-6*(here->DIOtemp-REFTEMP)-gmaSWnew);
-            } else if (model->DIOtlevc == 1) {
-                    here->DIOtJctSWPot = model->DIOjunctionSWPot - model->DIOtphp*(here->DIOtemp-REFTEMP);
-                    here->DIOtJctSWCap = here->DIOjunctionSWCap *
-                            (1+model->DIOctp*(here->DIOtemp-REFTEMP));
-            }
-
-            here->DIOtSatCur = model->DIOsatCur * here->DIOarea * exp(
-                    ((here->DIOtemp/model->DIOnomTemp)-1) *
-                    model->DIOactivationEnergy/(model->DIOemissionCoeff*vt) +
-                    model->DIOsaturationCurrentExp/model->DIOemissionCoeff *
-                    log(here->DIOtemp/model->DIOnomTemp) );
-            here->DIOtSatSWCur = model->DIOsatSWCur * here->DIOpj * exp(
-                    ((here->DIOtemp/model->DIOnomTemp)-1) *
-                    model->DIOactivationEnergy/(model->DIOswEmissionCoeff*vt) +
-                    model->DIOsaturationCurrentExp/model->DIOswEmissionCoeff *
-                    log(here->DIOtemp/model->DIOnomTemp) );
-
-            here->DIOtTunSatCur = model->DIOtunSatCur * here->DIOarea * exp(
-                    ((here->DIOtemp/model->DIOnomTemp)-1) *
-                    model->DIOtunEGcorrectionFactor*model->DIOactivationEnergy/vt +
-                    model->DIOtunSaturationCurrentExp *
-                    log(here->DIOtemp/model->DIOnomTemp) );
-            here->DIOtTunSatSWCur = model->DIOtunSatSWCur * here->DIOpj * exp(
-                    ((here->DIOtemp/model->DIOnomTemp)-1) *
-                    model->DIOtunEGcorrectionFactor*model->DIOactivationEnergy/vt +
-                    model->DIOtunSaturationCurrentExp *
-                    log(here->DIOtemp/model->DIOnomTemp) );
-
-            here->DIOtRecSatCur = model->DIOrecSatCur * here->DIOarea * exp(
-                    ((here->DIOtemp/model->DIOnomTemp)-1) *
-                    model->DIOactivationEnergy/(model->DIOrecEmissionCoeff*vt) +
-                    model->DIOsaturationCurrentExp/model->DIOrecEmissionCoeff *
-                    log(here->DIOtemp/model->DIOnomTemp) );
-
-            /* the defintion of f1, just recompute after temperature adjusting
-             * all the variables used in it */
-            here->DIOtF1=here->DIOtJctPot*
-                    (1-exp((1-here->DIOtGradingCoeff)*xfc))/
-                    (1-here->DIOtGradingCoeff);
-            /* same for Depletion Capacitance */
-            here->DIOtDepCap=model->DIOdepletionCapCoeff*
-                    here->DIOtJctPot;
-            here->DIOtDepSWCap=model->DIOdepletionSWcapCoeff*
-                    here->DIOtJctSWPot;
-            /* and Vcrit */
-            vte=model->DIOemissionCoeff*vt;
-
-            here->DIOtVcrit = vte * log(vte/(CONSTroot2*here->DIOtSatCur));
-
-            /* limit junction potential to max of 1/FC */
-            if(here->DIOtDepCap > 1.0) {
-                here->DIOtJctPot=1.0/model->DIOdepletionCapCoeff;
-                here->DIOtDepCap=model->DIOdepletionCapCoeff*here->DIOtJctPot;
-                SPfrontEnd->IFerrorf (ERR_WARNING,
-                        "%s: junction potential VJ too large, limited to %f",
-                        model->DIOmodName, here->DIOtJctPot);
-            }
-            /* limit sidewall junction potential to max of 1/FCS */
-            if(here->DIOtDepSWCap > 1.0) {
-                here->DIOtJctSWPot=1.0/model->DIOdepletionSWcapCoeff;
-                here->DIOtDepSWCap=model->DIOdepletionSWcapCoeff*here->DIOtJctSWPot;
-                SPfrontEnd->IFerrorf (ERR_WARNING,
-                        "%s: junction potential VJS too large, limited to %f",
-                        model->DIOmodName, here->DIOtJctSWPot);
-            }
-
-            /* and now to compute the breakdown voltage, again, using
-             * temperature adjusted basic parameters */
-            if (model->DIObreakdownVoltageGiven){
-                if (model->DIOtlev == 0) {
-                    tBreakdownVoltage = model->DIObreakdownVoltage - model->DIOtcv * dt;
-                } else {
-                    tBreakdownVoltage = model->DIObreakdownVoltage * (1 - model->DIOtcv * dt);
-                }
-                if (model->DIOlevel == 1) {
-                    cbv = model->DIObreakdownCurrent;
-                } else { /* level=3 */
-                    cbv = model->DIObreakdownCurrent * here->DIOarea;
-                }
-                if (cbv < here->DIOtSatCur * tBreakdownVoltage/vt) {
-                    cbv=here->DIOtSatCur * tBreakdownVoltage/vt;
-#ifdef TRACE
-                    SPfrontEnd->IFerrorf (ERR_WARNING, "%s: breakdown current increased to %g to resolve", here->DIOname, cbv);
-                    SPfrontEnd->IFerrorf (ERR_WARNING,
-                    "incompatibility with specified saturation current");
-#endif
-                    xbv=tBreakdownVoltage;
-                } else {
-                    tol=ckt->CKTreltol*cbv;
-                    xbv=tBreakdownVoltage-model->DIObrkdEmissionCoeff*vt*log(1+cbv/
-                            (here->DIOtSatCur));
-                    iter=0;
-                    for(iter=0 ; iter < 25 ; iter++) {
-                        xbv=tBreakdownVoltage-model->DIObrkdEmissionCoeff*vt*log(cbv/
-                                (here->DIOtSatCur)+1-xbv/vt);
-                        xcbv=here->DIOtSatCur *
-                             (exp((tBreakdownVoltage-xbv)/(model->DIObrkdEmissionCoeff*vt))-1+xbv/vt);
-                        if (fabs(xcbv-cbv) <= tol) goto matched;
-                    }
-#ifdef TRACE
-                    SPfrontEnd->IFerrorf (ERR_WARNING, "%s: unable to match forward and reverse diode regions: bv = %g, ibv = %g", here->DIOname, xbv, xcbv);
-#endif
-                }
-                matched:
-                    here->DIOtBrkdwnV = xbv;
-            }
-
-            /* transit time temperature adjust */
-            factor = 1.0 + (model->DIOtranTimeTemp1 * dt)
-                         + (model->DIOtranTimeTemp2 * dt * dt);
-            here->DIOtTransitTime = model->DIOtransitTime * factor;
-
-            /* Series resistance temperature adjust */
-            here->DIOtConductance = model->DIOconductance;
-            if(model->DIOresistGiven && model->DIOresist!=0.0) {
-                factor = 1.0 + (model->DIOresistTemp1) * dt
-                         + (model->DIOresistTemp2 * dt * dt);
-                here->DIOtConductance = model->DIOconductance / factor;
-            }
-
-            here->DIOtF2=exp((1+here->DIOtGradingCoeff)*xfc);
-            here->DIOtF3=1-model->DIOdepletionCapCoeff*
-                    (1+here->DIOtGradingCoeff);
-            here->DIOtF2SW=exp((1+model->DIOgradingSWCoeff)*xfcs);
-            here->DIOtF3SW=1-model->DIOdepletionSWcapCoeff*
-                    (1+model->DIOgradingSWCoeff);
+            DIOtempUpdate(model, here, here->DIOtemp, ckt);
 
         } /* instance */
 

src/spicelib/parser/inp2d.c

@@ -15,78 +15,100 @@ Modified: 2001 Paolo Nenzi (Cider Integration)
 void INP2D(CKTcircuit *ckt, INPtables * tab, struct card *current)
 {
 
-/* Dname <node> <node> <model> [<val>] [OFF] [IC=<val>] */
+/* Dname <node> <node> [<temp>] <model> [<val>] [OFF] [IC=<val>] */
 
-    int mytype;			/* the type we looked up */
-    int type;			/* the type the model says it is */
-    char *line;			/* the part of the current line left to parse */
-    char *name;			/* the resistor's name */
-    char *nname1;		/* the first node's name */
-    char *nname2;		/* the second node's name */
-    CKTnode *node1;		/* the first node's node pointer */
-    CKTnode *node2;		/* the second node's node pointer */
-    int error;			/* error code temporary */
-    GENinstance *fast;		/* pointer to the actual instance */
-    IFvalue ptemp;		/* a value structure to package resistance into */
-    int waslead;		/* flag to indicate that funny unlabeled number was found */
-    double leadval;		/* actual value of unlabeled number */
-    char *model;		/* the name of the model */
-    INPmodel *thismodel;	/* pointer to model description for user's model */
-    GENmodel *mdfast;		/* pointer to the actual model */
-    IFuid uid;			/* uid of default model */
+    int mytype;         /* the type we looked up */
+    int type;           /* the type the model says it is */
+    char *line;         /* the part of the current line left to parse */
+    char *name;         /* the resistor's name */
+    const int max_i = 3;/* the maximum node count */
+    CKTnode *node[3];   
+    int error;          /* error code temporary */
+    int numnodes;       /* flag indicating 4 or 5 nodes */
+    GENinstance *fast;  /* pointer to the actual instance */
+    IFvalue ptemp;      /* a value structure to package resistance into */
+    int waslead;        /* flag to indicate that funny unlabeled number was found */
+    double leadval;     /* actual value of unlabeled number */
+    INPmodel *thismodel;/* pointer to model description for user's model */
+    GENmodel *mdfast;   /* pointer to the actual model */
+    IFuid uid;          /* uid of default model */
+    int i;
 
     mytype = INPtypelook("Diode");
     if (mytype < 0) {
-	LITERR("Device type Diode not supported by this binary\n");
-	return;
+    LITERR("Device type Diode not supported by this binary\n");
+    return;
     }
     line = current->line;
     INPgetNetTok(&line, &name, 1);
     INPinsert(&name, tab);
-    INPgetNetTok(&line, &nname1, 1);
-    INPtermInsert(ckt, &nname1, tab, &node1);
-    INPgetNetTok(&line, &nname2, 1);
-    INPtermInsert(ckt, &nname2, tab, &node2);
-    INPgetNetTok(&line, &model, 1);
-    INPinsert(&model, tab);
-    current->error = INPgetMod(ckt, model, &thismodel, tab);
-    if (thismodel != NULL) {
-	if ((mytype != thismodel->INPmodType)
 
+    for (i = 0; ; i++) {
+        char *token;
+        INPgetNetTok(&line, &token, 1);
+        if (i >= 2 && INPlookMod(token)) {
+            INPinsert(&token, tab);
+            txfree(INPgetMod(ckt, token, &thismodel, tab));
+            if (!thismodel) {
+                LITERR ("Unable to find definition of given model");
+                return;
+            }
+            break;
+        }
+        if (i >= max_i) {
+            LITERR ("could not find a valid modelname");
+            return;
+        }
+        INPtermInsert(ckt, &token, tab, &node[i]);
+    }
+
+    if (i > max_i) {
+        LITERR("Too many nodes for this model type");
+        return;
+    }
+
+    numnodes = i;
+
+    if (thismodel != NULL) {
+        if ((mytype != thismodel->INPmodType)
 #ifdef CIDER
-                      && (thismodel->INPmodType != INPtypelook("NUMD"))
-                      && (thismodel->INPmodType != INPtypelook("NUMD2"))
+            && (thismodel->INPmodType != INPtypelook("NUMD"))
+            && (thismodel->INPmodType != INPtypelook("NUMD2"))
 #endif
-       ){
-	    LITERR("incorrect model type");
-	    return;
-	}
-	type = thismodel->INPmodType; /*HT 050903*/
-	mdfast = (thismodel->INPmodfast);
+        ) {
+            LITERR("incorrect model type");
+            return;
+        }
+        type = thismodel->INPmodType; /*HT 050903*/
+        mdfast = (thismodel->INPmodfast);
     } else {
-	type = mytype;
-	if (!tab->defDmod) {
-	    /* create default D model */
-	    IFnewUid(ckt, &uid, NULL, "D", UID_MODEL, NULL);
-	    IFC(newModel, (ckt, type, &(tab->defDmod), uid));
-	}
-	mdfast = tab->defDmod;
+        type = mytype;
+        if (!tab->defDmod) {
+            /* create default D model */
+            IFnewUid(ckt, &uid, NULL, "D", UID_MODEL, NULL);
+            IFC(newModel, (ckt, type, &(tab->defDmod), uid));
+        }
+        mdfast = tab->defDmod;
     }
+
     IFC(newInstance, (ckt, mdfast, &fast, name));
-    IFC(bindNode, (ckt, fast, 1, node1));
-    IFC(bindNode, (ckt, fast, 2, node2));
+    for (i = 0; i < max_i; i++)
+        if (i < numnodes)
+            IFC (bindNode, (ckt, fast, i + 1, node[i]));
+        else
+            GENnode(fast)[i] = -1;
+
     PARSECALL((&line, ckt, type, fast, &leadval, &waslead, tab));
     if (waslead) {
-
 #ifdef CIDER
-    	if( type == INPtypelook("NUMD2") ) {
+        if( type == INPtypelook("NUMD2") ) {
             LITERR(" error:  no unlabelled parameter permitted on NUMD2\n");
-	} else {
+    } else {
 #endif
-	ptemp.rValue = leadval;
-	GCA(INPpName, ("area", &ptemp, ckt, type, fast));
+        ptemp.rValue = leadval;
+        GCA(INPpName, ("area", &ptemp, ckt, type, fast));
     }
 #ifdef CIDER    
-  }
+    }
 #endif
 }