diff --git a/examples/soa/diode-soa-sh.cir b/examples/soa/diode-soa-sh.cir
new file mode 100644
index 000000000..042f6d5ad
--- /dev/null
+++ b/examples/soa/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/examples/soa/diode-soa-simple.cir b/examples/soa/diode-soa-simple.cir
new file mode 100644
index 000000000..f61d21639
--- /dev/null
+++ b/examples/soa/diode-soa-simple.cir
@@ -0,0 +1,37 @@
+SOA test for generic diode, including self-heating
+* forward direction
+* ngspice-35
+
+.temp  100
+
+vtamb tamb 0 27
+
+v1 1 0 0.7
+R1 1 0 100 ; just a parallel resistor
+Dth1 1 0 tj dmod1 thermal
+Rtj tj tamb 100 ; the thermal resistance junction to ambient
+* diode model parameters include all SOA parameters
+.model dmod1 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)
+
+D2 1 0 dmod2 temp=125
+* diode model parameters include all SOA parameters
+.model dmod2 d (rs=200m bv=21 rth0=100 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  @dth1[id] @d2[id] all
+*dc v1 3 -22.5 -0.5
+dc v1 0.02 2 0.02
+*display
+
+set xbrushwidth=3
+plot @dth1[id] @d2[id] loglog ylimit 10m 10 xlimit 0.1 1
+
+
+settype temperature tj
+plot tj xlimit 0.5 1 ylimit 0 300
+
+.endc
+
+.end
diff --git a/examples/soa/nic_soa.cir b/examples/soa/nic_soa.cir
new file mode 100644
index 000000000..5752d99d7
--- /dev/null
+++ b/examples/soa/nic_soa.cir
@@ -0,0 +1,39 @@
+CMOS NIC
+*
+.subckt osc_cmos ib_osz lc ra vdd vss
+m16 ib_osz ib_osz vss vss n1 w=20u l=1u m=8
+m15 ra ib_osz vss vss n1 w=20u l=1u m=2
+m8 net99 net95 ra ra n1 w=20u l=1u m=2
+m1 net95 net95 net93 net93 n1 w=20u l=1u m=2
+m25 net99 net99 vdd vdd p1 w=3.3u l=0.5u m=1
+m5 net99 net99 vdd vdd p1 w=20u l=1u m=5
+m4 net95 net99 vdd vdd p1 w=20u l=1u m=5
+r23 net99 vss r=38K
+r18 net93 lc r=10
+.ends osc_cmos
+*
+.subckt psens LC
+R1 LC P001 40K
+L1 LC P002 14.9u
+R2 P002 0 0.55
+L2 P001 0 1.4m
+.ends psens
+*
+xi36 bias lc ra vdd 0 osc_cmos
+v39 vdd 0 dc=3.5 pulse ( 0 3.5 10u 10n 10n 1 2 )
+r4 ra 0 3.972K
+c23 lc 0 1.8n
+i37 vdd bias dc=1u
+*
+xi18 lc psens
+*
+.option warn=1
+.control
+tran 1u 1m 0 50n
+plot v(LC)
+.endc
+*
+.include modelcard.nmos
+.include modelcard.pmos
+*
+.end