From 6bafdf70d5965f31c17cf7a829d53e177221ce46 Mon Sep 17 00:00:00 2001
From: dwarning <dwarning>
Date: Wed, 5 May 2021 22:53:39 +0200
Subject: [PATCH] sub-coll junction calculation only if ISS given

---
 src/spicelib/devices/bjt/bjtload.c | 34 ++++++++++++++++++------------
 src/spicelib/devices/bjt/bjttemp.c | 23 ++++++++++----------
 2 files changed, 32 insertions(+), 25 deletions(-)

diff --git a/src/spicelib/devices/bjt/bjtload.c b/src/spicelib/devices/bjt/bjtload.c
index 9f9317f12..ca14ff620 100644
--- a/src/spicelib/devices/bjt/bjtload.c
+++ b/src/spicelib/devices/bjt/bjtload.c
@@ -434,8 +434,13 @@ BJTload(GENmodel *inModel, CKTcircuit *ckt)
                 vbc = DEVpnjlim(vbc,*(ckt->CKTstate0 + here->BJTvbc),vt,
                         here->BJTtVcrit,&ichk1);
                 if (ichk1 == 1) icheck=1;
-                vsub = DEVpnjlim(vsub,*(ckt->CKTstate0 + here->BJTvsub),vt,
-                        here->BJTtSubVcrit,&ichk1);
+                if (model->BJTsubSatCurGiven) {
+                    vsub = DEVpnjlim(vsub,*(ckt->CKTstate0 + here->BJTvsub),vt,
+                            here->BJTtSubVcrit,&ichk1);
+                } else {
+                    vsub = DEVpnjlim(vsub,*(ckt->CKTstate0 + here->BJTvsub),vt,
+                            50,&ichk1);
+                }
                 if (ichk1 == 1) icheck=1;
                 vrci = vbc - vbcx; /* in case vbc was limited */
             }
@@ -502,19 +507,22 @@ next1:      vtn=vt*here->BJTtemissionCoeffF;
             gbcn+=ckt->CKTgmin;
             cbcn+=ckt->CKTgmin*vbc;
 
-            vts=vt*here->BJTtemissionCoeffS;
-
-            if(vsub <= -3*vts) {
-                arg=3*vts/(vsub*CONSTe);
-                arg = arg * arg * arg;
-                gdsub = csubsat*3*arg/vsub+ckt->CKTgmin;
-                cdsub = -csubsat*(1+arg)+ckt->CKTgmin*vsub;
+            if (model->BJTsubSatCurGiven) {
+                vts=vt*here->BJTtemissionCoeffS;
+                if(vsub <= -3*vts) {
+                    arg=3*vts/(vsub*CONSTe);
+                    arg = arg * arg * arg;
+                    gdsub = csubsat*3*arg/vsub+ckt->CKTgmin;
+                    cdsub = -csubsat*(1+arg)+ckt->CKTgmin*vsub;
+                } else {
+                    evsub = exp(MIN(MAX_EXP_ARG,vsub/vts));
+                    gdsub = csubsat*evsub/vts + ckt->CKTgmin;
+                    cdsub = csubsat*(evsub-1) + ckt->CKTgmin*vsub;
+                }
             } else {
-                evsub = exp(MIN(MAX_EXP_ARG,vsub/vts));
-                gdsub = csubsat*evsub/vts + ckt->CKTgmin;
-                cdsub = csubsat*(evsub-1) + ckt->CKTgmin*vsub;
+                gdsub = ckt->CKTgmin;
+                cdsub = ckt->CKTgmin*vsub;
             }
-
             /*
              *   Kull's Quasi-Saturation model
              */
diff --git a/src/spicelib/devices/bjt/bjttemp.c b/src/spicelib/devices/bjt/bjttemp.c
index be4dfa7aa..d74217f73 100644
--- a/src/spicelib/devices/bjt/bjttemp.c
+++ b/src/spicelib/devices/bjt/bjttemp.c
@@ -170,12 +170,14 @@ BJTtemp(GENmodel *inModel, CKTcircuit *ckt)
             factlog = ratio1 * model->BJTenergyGap/vt +
                     model->BJTtempExpIS*ratlog;
             if ((model->BJTtlev == 0) || (model->BJTtlev == 1)) {
-              factor = exp(factlog);
-              here->BJTtSatCur = model->BJTsatCur * factor;
-              here->BJTtSubSatCur = model->BJTsubSatCur * factor;
+                factor = exp(factlog);
+                here->BJTtSatCur = model->BJTsatCur * factor;
+                if (model->BJTsubSatCurGiven)
+                    here->BJTtSubSatCur = model->BJTsubSatCur * factor;
             } else if (model->BJTtlev == 3) {
-              here->BJTtSatCur = pow(model->BJTsatCur,(1+model->BJTtis1*dt+model->BJTtis2*dt*dt));
-              here->BJTtSubSatCur = pow(model->BJTsubSatCur,(1+model->BJTtiss1*dt+model->BJTtiss2*dt*dt));
+                here->BJTtSatCur = pow(model->BJTsatCur,(1+model->BJTtis1*dt+model->BJTtis2*dt*dt));
+                if (model->BJTsubSatCurGiven)
+                    here->BJTtSubSatCur = pow(model->BJTsubSatCur,(1+model->BJTtiss1*dt+model->BJTtiss2*dt*dt));
             }
 
             if (model->BJTintCollResistGiven) {
@@ -230,7 +232,7 @@ BJTtemp(GENmodel *inModel, CKTcircuit *ckt)
                     (4e-4*(here->BJTtemp-REFTEMP)-gmanew);
             } else if (model->BJTtlevc == 1) {
                 here->BJTtBEcap = model->BJTdepletionCapBE*
-                    (1+model->BJTcte*dt);			
+                    (1+model->BJTcte*dt);           
                 here->BJTtBEpot = model->BJTpotentialBE - model->BJTtvje*dt;
             }
             if (model->BJTtlevc == 0) {
@@ -245,7 +247,7 @@ BJTtemp(GENmodel *inModel, CKTcircuit *ckt)
                     (4e-4*(here->BJTtemp-REFTEMP)-gmanew);
             } else if (model->BJTtlevc == 1) {
                 here->BJTtBCcap = model->BJTdepletionCapBC*
-                    (1+model->BJTctc*dt);			
+                    (1+model->BJTctc*dt);           
                 here->BJTtBCpot = model->BJTpotentialBC - model->BJTtvjc*dt;
             }
             if (model->BJTtlevc == 0) {
@@ -260,7 +262,7 @@ BJTtemp(GENmodel *inModel, CKTcircuit *ckt)
                     (4e-4*(here->BJTtemp-REFTEMP)-gmanew);
             } else if (model->BJTtlevc == 1) {
                 here->BJTtSubcap = model->BJTcapSub*
-                    (1+model->BJTcts*dt);			
+                    (1+model->BJTcts*dt);           
                 here->BJTtSubpot = model->BJTpotentialSubstrate - model->BJTtvjs*dt;
             }
 
@@ -274,12 +276,9 @@ BJTtemp(GENmodel *inModel, CKTcircuit *ckt)
                     (1 - here->BJTtjunctionExpBC);
             here->BJTtVcrit = vt *
                      log(vt / (CONSTroot2*here->BJTtSatCur*here->BJTarea));
-            if(model->BJTsubSatCurGiven) {
+            if (model->BJTsubSatCurGiven)
                 here->BJTtSubVcrit = vt *
                          log(vt / (CONSTroot2*here->BJTtSubSatCur*here->BJTarea));
-            } else {
-                here->BJTtSubVcrit = here->BJTtVcrit;
-            }
             here->BJTtf2 = exp((1 + here->BJTtjunctionExpBE) * xfc);
             here->BJTtf3 = 1 - model->BJTdepletionCapCoeff *
                     (1 + here->BJTtjunctionExpBE);