diff --git a/examples/various/diode-soa-sh.cir b/examples/various/diode-soa-sh.cir
new file mode 100644
index 000000000..042f6d5ad
--- /dev/null
+++ b/examples/various/diode-soa-sh.cir
@@ -0,0 +1,117 @@
+SOA test for generic diode, including self-heating
+* forward direction
+* ngspice-35
+
+*.temp  200
+
+vtamb tamb 0 27
+
+v1 1 0 0.7
+R1 1 0 100 ; just a parallel resistor
+D1 1 0 tj dmod thermal
+Rtj tj tamb 100 ; the thermal resistance junction to ambient
+* diode model parameters include all SOA parameters
+.model dmod d (rs=200m bv=21 rth0=1e6 tnom=25 fv_max=1.5 bv_max=20 id_max=1.5 pd_max=1 te_max=175)
+
+.option warn=1 maxwarns=2
+
+.control
+save  @d1[id] all
+*dc v1 3 -22.5 -0.5
+dc v1 0.02 2 0.02
+*display
+set xbrushwidth=3
+
+* get data from diode model
+let pdmax =  @dmod[pd_max]
+let idmax = @dmod[id_max]
+let vmax = @dmod[fv_max]
+let tmax = @dmod[te_max]
+let tnom = @dmod[tnom]
+
+let iid = @d1[id]
+let ilen = length(iid)
+let soa = unitvec(ilen) * idmax
+* the current power dissipation in the diode
+let pd=@d1[id]*v(1) + 1p ; 1p for log scale, avoid 0
+
+* plot the static SOA diagram
+* no self heating
+let i = 0
+while i < ilen
+* power limit
+  let pp = soa[i] * v(1)[i]
+  if pp > pdmax
+    let soa[i] = soa[i] * pdmax / pp
+  end
+* voltage limit
+  if v(1)[i] > vmax
+    let soa[i] = 1p
+  end
+* temperature limit
+  let tcur = pp * @Rtj[r] + v(tamb)
+  if tcur[i] > tmax
+    let soa[i] = 1p
+  end  
+  let i = i + 1
+end
+
+settype current iid soa
+plot iid soa loglog ylimit 10m 10 xlimit 0.1 1 title 'Diode SOA (safe operating area, no self-heating)' ylabel 'Diode current' xlabel 'Diode forward voltage'
+
+
+
+unlet pdmax
+let pdmax = @dmod[pd_max] - (v(tj) - tnom) / @Rtj[r]
+let tdio = v(tj)
+
+echo
+*echo pdmax $&pdmax 
+*echo temp $&tdio 
+*echo tnom $&tnom
+echo
+
+let plen = length(pdmax)
+let i = 0
+while i < plen
+  if pdmax[i] < 0
+    let pdmax[i] = 1p
+  end
+  let i = i + 1
+end
+
+* plot the static SOA diagram
+* now with self heating
+let i = 0
+while i < ilen
+* power limit
+  let pp = soa[i] * v(1)[i]
+  if pp > pdmax[i]
+    let soa[i] = soa[i] * pdmax[i] / pp
+  end
+* voltage limit
+  if v(1)[i] > vmax
+    let soa[i] = 1p
+  end
+* temperature limit
+  let tcur = pp * @Rtj[r] + v(tamb)
+  if tcur[i] > tmax
+    let soa[i] = 1p
+  end  
+  let i = i + 1
+end
+
+settype current iid soa
+plot iid soa loglog ylimit 10m 10 xlimit 0.1 1  title 'Diode SOA (safe operating area, including self-heating)' ylabel 'Diode current' xlabel 'Diode forward voltage'
+
+*settype power pd pdmax
+*plot pd pdmax loglog ylimit 1m 10 xlimit 0.1 1
+
+*settype temperature tj
+*plot tj
+
+*plot pd vs tj pdmax vs tj
+
+.endc
+
+.end
diff --git a/src/spicelib/devices/dio/diosoachk.c b/src/spicelib/devices/dio/diosoachk.c
index c9e1f23e7..27d2dc0df 100644
--- a/src/spicelib/devices/dio/diosoachk.c
+++ b/src/spicelib/devices/dio/diosoachk.c
@@ -20,6 +20,7 @@ DIOsoaCheck(CKTcircuit *ckt, GENmodel *inModel)
     double vd;  /* current diode voltage */
     double id;  /* current diode current */
     double pd;  /* current diode power */
+    double pd_max;  /* maximum diode power */
     double te;  /* current diode temperature */
     int maxwarns;
     static int warns_fv = 0, warns_bv = 0, warns_id = 0, warns_pd = 0, warns_te = 0;
@@ -45,7 +46,7 @@ DIOsoaCheck(CKTcircuit *ckt, GENmodel *inModel)
             if (vd > model->DIOfv_max)
                 if (warns_fv < maxwarns) {
                     soa_printf(ckt, (GENinstance*) here,
-                               "Vd=%g V has exceeded Fv_max=%g V\n",
+                               "Vd=%.4g V has exceeded Fv_max=%.4g V\n",
                                vd, model->DIOfv_max);
                     warns_fv++;
                 }
@@ -53,7 +54,7 @@ DIOsoaCheck(CKTcircuit *ckt, GENmodel *inModel)
             if (-vd > model->DIObv_max)
                 if (warns_bv < maxwarns) {
                     soa_printf(ckt, (GENinstance*) here,
-                               "Vd=%g V has exceeded Bv_max=%g V\n",
+                               "Vd=%.4g V has exceeded Bv_max=%.4g V\n",
                                vd, model->DIObv_max);
                     warns_bv++;
                 }
@@ -62,7 +63,7 @@ DIOsoaCheck(CKTcircuit *ckt, GENmodel *inModel)
             if (id > fabs(model->DIOid_max))
                 if (warns_id < maxwarns) {
                     soa_printf(ckt, (GENinstance*) here,
-                               "Id=%g A at Vd=%g V has exceeded Id_max=%g A\n",
+                               "Id=%.4g A at Vd=%.4g V has exceeded Id_max=%.4g A\n",
                                id, vd, model->DIOid_max);
                     warns_id++;
                 }
@@ -72,25 +73,65 @@ DIOsoaCheck(CKTcircuit *ckt, GENmodel *inModel)
                       *(ckt->CKTstate0 + here->DIOvoltage) +
                       *(ckt->CKTstate0 + here->DIOcurrent) *
                       *(ckt->CKTstate0 + here->DIOcurrent) / here->DIOtConductance);
-            if (pd > fabs(model->DIOpd_max))
-                if (warns_pd < maxwarns) {
-                    soa_printf(ckt, (GENinstance*) here,
-                               "Pd=%g W at Vd=%g V has exceeded Pd_max=%g W\n",
-                               pd, vd, model->DIOpd_max);
-                    warns_pd++;
-                }
 
-            te = here->DIOtemp - CONSTCtoK;
-            if (te > model->DIOte_max)
-                if (warns_te < maxwarns) {
-                    soa_printf(ckt, (GENinstance*) here,
-                               "Te=%g C at Vd=%g V has exceeded te_max=%g C\n",
-                               te, vd, model->DIOte_max);
-                    warns_te++;
+            /* calculate max power including derating:
+               up to tnom the derating is zero,
+               at maximum temp allowed the derating is 100%.
+               Device temperature by self-heating or given externally. */
+            if (here->DIOthermal && model->DIOrth0Given && model->DIOpd_maxGiven
+                && model->DIOte_maxGiven && model->DIOnomTempGiven) {
+                te = ckt->CKTrhsOld[here->DIOtempNode];
+                if (te < model->DIOnomTemp)
+                    pd_max = model->DIOpd_max;
+                else {
+                    pd_max = model->DIOpd_max - (te - model->DIOnomTemp) / model->DIOrth0;
+                    pd_max = (pd_max > 0) ? pd_max : 0.;
                 }
+                if (pd > pd_max)
+                    if (warns_pd < maxwarns) {
+                        soa_printf(ckt, (GENinstance*)here,
+                            "Pd=%.4g W at Vd=%.4g V and Te=%.4g C has exceeded Pd_max=%.4g W\n",
+                            pd, vd, te, pd_max);
+                        warns_pd++;
+                    }
+                if (te > model->DIOte_max)
+                    if (warns_te < maxwarns) {
+                        soa_printf(ckt, (GENinstance*)here,
+                            "Te=%.4g C at Vd=%.4g V has exceeded te_max=%.4g C\n",
+                            te, vd, model->DIOte_max);
+                        warns_te++;
+                    }
 
+            }
+            /* derating without self-heating, external temp given */
+            else if (!here->DIOthermal && here->DIOtempGiven && model->DIOrth0Given && model->DIOpd_maxGiven 
+                                 && model->DIOte_maxGiven && model->DIOnomTempGiven) {
+                if (here->DIOtemp < model->DIOnomTemp)
+                    pd_max = model->DIOpd_max;
+                else {
+                    pd_max = model->DIOpd_max - (here->DIOtemp - model->DIOnomTemp) / model->DIOrth0;
+                    pd_max = (pd_max > 0) ? pd_max : 0.;
+                }
+                if (pd > pd_max)
+                    if (warns_pd < maxwarns) {
+                        soa_printf(ckt, (GENinstance*)here,
+                            "Pd=%.4g W at Vd=%.4g V and Te=%.4g C has exceeded Pd_max=%.4g W\n",
+                            pd, vd, here->DIOtemp - CONSTCtoK, pd_max);
+                        warns_pd++;
+                    }
+            }
+            /* no derating */
+            else {
+                pd_max = model->DIOpd_max;
+                if (pd > pd_max)
+                    if (warns_pd < maxwarns) {
+                        soa_printf(ckt, (GENinstance*)here,
+                            "Pd=%.4g W at Vd=%.4g V has exceeded Pd_max=%.4g W\n",
+                            pd, vd, pd_max);
+                        warns_pd++;
+                    }
+            }
         }
-
     }
 
     return OK;