diff --git a/examples/Monte_Carlo/mc_ring_circ.net b/examples/Monte_Carlo/mc_ring_circ.net
new file mode 100644
index 000000000..cd058a5eb
--- /dev/null
+++ b/examples/Monte_Carlo/mc_ring_circ.net
@@ -0,0 +1,56 @@
+Perform Monte Carlo simulation in ngspice
+* 25 stage Ring-Osc. BSIM3 or 4 with statistical variation of model parameters
+* Model parameters are varied according to the PDK selection.
+* Tested with 3 different commercial HSPICE libraries from 2 vendors.
+* To be started with script MC_ring_ts.sp
+
+.options noacct seedinfo
+
+vin in out dc 0.5 pulse 0.5 0 0.1n 5n 1 1 1
+vdd dd 0 dc 3.3
+vss ss 0 dc 0
+ve  sub  0 dc 0
+vpe well 0 dc 3.3
+
+* transistors to be selected according to the library (here: p33ll and n33ll or pch_5_mac and nch_5_mac
+* or pe3 and ne3 or p1 and n1 (these models see below))
+.subckt inv1 dd ss sub well in out
+*XMP1 out in  dd  well p33ll w=5u l=800n m=3 nf=1 ad=1.35p as=1.35p pd=9.6u ps=9.6u mosmis_mod=1
+*XMN1 out in  ss  sub n33ll w=5u l=800n m=1 nf=3 ad=0.9p as=0.9p pd=6.6u ps=6.6u mosmis_mod=1
+*XMP1 out in  dd  well pch_5_mac w=5u l=800n m=3 nf=1 ad=1.35p as=1.35p pd=9.6u ps=9.6u mosmis_mod=1
+*XMN1 out in  ss  sub nch_5_mac w=5u l=800n m=1 nf=3 ad=0.9p as=0.9p pd=6.6u ps=6.6u mosmis_mod=1
+*XMP1 out in  dd  well pe3 w=5u l=800n m=3 nf=1 ad=1.35p as=1.35p pd=9.6u ps=9.6u mosmis_mod=1
+*XMN1 out in  ss  sub ne3 w=5u l=800n m=1 nf=3 ad=0.9p as=0.9p pd=6.6u ps=6.6u mosmis_mod=1
+MP1 out in  dd  well p1 w=5u l=800n m=3 ad=1.35p as=1.35p pd=9.6u ps=9.6u
+MN1 out in  ss  sub n1 w=5u l=800n m=1 ad=0.9p as=0.9p pd=6.6u ps=6.6u
+.ends inv1
+
+.subckt inv5 dd ss sub well in out
+xinv1 dd ss sub well in 1 inv1
+xinv2 dd ss sub well 1  2 inv1
+xinv3 dd ss sub well 2  3 inv1
+xinv4 dd ss sub well 3  4 inv1
+xinv5 dd ss sub well 4 out inv1
+.ends inv5
+
+xinv1 dd ss sub well in out5 inv5
+xinv2 dd ss sub well out5 out10 inv5
+xinv3 dd ss sub well out10 out15 inv5
+xinv4 dd ss sub well out15 out20 inv5
+xinv5 dd ss sub well out20 out inv5
+xinv11 dd 0 sub well out buf inv1
+cout  buf ss 0.2pF
+
+ *** Model library files.
+* Add your library here
+* Chose the transistors for XMP1 and XMN1 accordingly
+*.lib "jc_usage.l" MC_LIB
+*.lib "my_ts_usage.l" MC_LIB
+*.lib "x_usage.l" MC_LIB
+
+* or use the BSIM3 model with internal parameters except Vth0
+* that varies the threshold voltage +-3 sigma around a mean of +-0.6V
+.model p1 PMOS version=3.3.0 Level=8 Vth0=agauss(-0.6, 0.1, 3)
+.model n1 NMOS version=3.3.0 Level=8 Vth0=agauss(0.6, 0.1, 3)
+
+.end
diff --git a/examples/Monte_Carlo/mc_ring_lib_complete_actual.cir b/examples/Monte_Carlo/mc_ring_lib_complete_actual.cir
new file mode 100644
index 000000000..9fd239e47
--- /dev/null
+++ b/examples/Monte_Carlo/mc_ring_lib_complete_actual.cir
@@ -0,0 +1,212 @@
+Perform Monte Carlo simulation in ngspice
+* 25 stage Ring-Osc. BSIM3 or 4 with statistical variation of model parameters
+* Model parameters are varied according to the PDK selection.
+* Tested with 3 different commercial HSPICE libraries from 2 vendors.
+* Add your library to mc_ring_circ.net and choose transistors accordingly.
+* Add the library path to the .LIB statement.
+* A simple BSIM3 inverter R.O. serves as an MC example.
+
+.options noacct
+
+vin in out dc 0.5 pulse 0.5 0 0.1n 5n 1 1 1
+vdd dd 0 dc 3.3
+vss ss 0 dc 0
+ve  sub  0 dc 0
+vpe well 0 dc 3.3
+
+* transistors to be selected according to the library (here: p33ll and n33ll or pch_5_mac and nch_5_mac
+* or pe3 and ne3 or p1 and n1 (these models see below))
+.subckt inv1 dd ss sub well in out
+*XMP1 out in  dd  well p33ll w=5u l=800n m=3 nf=1 ad=1.35p as=1.35p pd=9.6u ps=9.6u mosmis_mod=1
+*XMN1 out in  ss  sub n33ll w=5u l=800n m=1 nf=3 ad=0.9p as=0.9p pd=6.6u ps=6.6u mosmis_mod=1
+XMP1 out in  dd  well pch_5_mac w=5u l=800n m=3 nf=1 ad=1.35p as=1.35p pd=9.6u ps=9.6u mosmis_mod=1
+XMN1 out in  ss  sub nch_5_mac w=5u l=800n m=1 nf=3 ad=0.9p as=0.9p pd=6.6u ps=6.6u mosmis_mod=1
+*XMP1 out in  dd  well pe3 w=5u l=800n m=3 nf=1 ad=1.35p as=1.35p pd=9.6u ps=9.6u mosmis_mod=1
+*XMN1 out in  ss  sub ne3 w=5u l=800n m=1 nf=3 ad=0.9p as=0.9p pd=6.6u ps=6.6u mosmis_mod=1
+*MP1 out in  dd  well p1 w=5u l=800n m=3 ad=1.35p as=1.35p pd=9.6u ps=9.6u
+*MN1 out in  ss  sub n1 w=5u l=800n m=1 ad=0.9p as=0.9p pd=6.6u ps=6.6u
+.ends inv1
+
+.subckt inv5 dd ss sub well in out
+xinv1 dd ss sub well in 1 inv1
+xinv2 dd ss sub well 1  2 inv1
+xinv3 dd ss sub well 2  3 inv1
+xinv4 dd ss sub well 3  4 inv1
+xinv5 dd ss sub well 4 out inv1
+.ends inv5
+
+xinv1 dd ss sub well in out5 inv5
+xinv2 dd ss sub well out5 out10 inv5
+xinv3 dd ss sub well out10 out15 inv5
+xinv4 dd ss sub well out15 out20 inv5
+xinv5 dd ss sub well out20 out inv5
+xinv11 dd 0 sub well out buf inv1
+cout  buf ss 0.2pF
+
+ *** Model library files.
+* Add your library here (full path required, or path relative to path
+* of ngspice executable (interactive mode), or relative to path of
+* input file (batch mode))
+* Chose the transistors for XMP1 and XMN1 according to the library
+*.lib "jc_usage.l" MC_LIB
+*.lib "../../../various/lib-test/my_usage.l" MC_LIB
+.lib "D:\Spice_general\tests\lib-test\ts14\my_ts_usage.l" MC_LIB
+*.lib "x_usage.l" MC_LIB
+
+* or use the BSIM3 model with internal parameters except Vth0
+* that varies the threshold voltage +-3 sigma around a mean of +-0.6V
+*.model p1 PMOS version=3.3.0 Level=8 Vth0=agauss(-0.6, 0.1, 3)
+*.model n1 NMOS version=3.3.0 Level=8 Vth0=agauss(0.6, 0.1, 3)
+
+.control
+  let mc_runs = 10               $ number of runs for monte carlo
+  let run = 0                     $ number of actual run
+  set curplot = new               $ create a new plot
+  set curplottitle = "Transient outputs"
+  set plot_out = $curplot         $ store its name to 'plot_out'
+  set curplot = new               $ create a new plot
+  set curplottitle = "FFT outputs"
+  set plot_fft = $curplot         $ store its name to 'plot_fft'
+  set curplot = new               $ create a new plot
+  set curplottitle = "Oscillation frequency"
+  set max_fft = $curplot          $ store its name to 'max_fft'
+  let mc_runsp = mc_runs + 1
+  let maxffts = unitvec(mc_runsp) $ vector for storing max measure results
+  let halfffts = unitvec(mc_runsp)$ vector for storing measure results at -40dB rising
+  unlet mc_runsp
+
+  set mc_runs = $&mc_runs        $ create a variable from the vector
+  let seeds = mc_runs + 2
+  setseed $&seeds
+  unlet seeds
+
+  save buf                        $ we just need buf, save memory by more than 10x
+
+* run the simulation loop
+
+* We have to figure out what to do if a single simulation will not converge.
+* There is now the variable sim_status, that is 0 if simulation ended regularly,
+* and 1 if the simulation has been aborted with error message '...simulation(s) aborted'.
+* Then we skip the rest of the run and continue with a new run.
+
+  dowhile run <= mc_runs
+
+    set run = $&run              $ create a variable from the vector
+
+    * run=0 simulates with nominal parameters
+    if run > 0
+      echo
+      echo * * * * * *
+      echo Source the circuit again internally for run no. $run
+      echo * * * * * *
+      setseed $run
+      mc_source  $ re-source the input file
+    else
+      echo run no. $run
+    end
+    echo simulation run no. $run of $mc_runs
+    tran 100p 1000n 0
+    echo Simulation status $sim_status
+    let simstat = $sim_status
+    if simstat = 1
+      if run = mc_runs
+        echo go to end
+      else
+        echo go to next run
+      end
+      destroy $curplot
+      goto next
+    end
+
+* select stop and step so that number of data points after linearization is not too
+* close to 8192, which would yield varying number of line length and thus scale for fft.
+*
+    set dt0 = $curplot
+    * save the linearized data for having equal time scales for all runs
+    linearize buf                 $ linearize only buf, no other vectors needed
+    set dt1 = $curplot             $ store the current plot to dt (tran i+1)
+    setplot $plot_out             $ make 'plt_out' the active plot
+    * firstly save the time scale once to become the default scale
+    if run=0
+       let time={$dt1}.time
+    end
+    let vout{$run}={$dt1}.buf     $ store the output vector to plot 'plot_out'
+    setplot $dt1                  $ go back to the previous plot (tran i+1)
+    fft buf $ run fft on vector buf
+    let buf2=db(mag(buf))
+    * find the frequency where buf has its maximum of the fft signal
+    meas sp fft_max MAX_AT buf2 from=0.05G to=0.7G
+    * find the frequency where buf is -40dB at rising fft signal
+    meas sp fft_40 WHEN buf2=-40 RISE=1 from=0.05G to=0.7G
+    * store the fft vector
+    set dt2 = $curplot            $ store the current plot to dt (spec i)
+    setplot $plot_fft             $ make 'plot_fft' the active plot
+    if run=0
+       let frequency={$dt2}.frequency
+    end
+    let fft{$run}={$dt2}.buf      $ store the output vector to plot 'plot_fft'
+    * store the measured value
+    setplot $max_fft              $ make 'max_fft' the active plot
+    let maxffts[{$run}]={$dt2}.fft_max
+    let halfffts[{$run}]={$dt2}.fft_40
+    destroy $dt0  $dt1  $dt2      $ save memory, we don't need this plot (spec) any more
+
+    label next
+    remcirc
+    let run = run + 1
+  end
+***** plotting **********************************************************
+if $?batchmode
+  echo
+  echo Plotting not available in batch mode
+  echo Write linearized vout0 to vout{$mc_runs} to rawfile $rawfile
+  echo
+  write $rawfile {$plot_out}.allv
+  rusage
+  quit
+else
+  plot {$plot_out}.vout0          $ just plot the tran output with run 0 parameters
+  setplot $plot_fft
+  plot db(mag(ally)) xlimit 0 1G ylimit -80 10
+*
+* create a histogram from vector maxffts
+  setplot $max_fft                $ make 'max_fft' the active plot
+  set startfreq=50MEG
+  set bin_size=1MEG
+  set bin_count=100
+  compose osc_frequ start=$startfreq step=$bin_size lin=$bin_count $ requires variables as parameters
+  settype frequency osc_frequ
+  let bin_count=$bin_count        $ create a vector from the variable
+  let yvec=unitvec(bin_count)     $ requires vector as parameter
+  let startfreq=$startfreq
+  let bin_size=$bin_size
+  * put data into the correct bins
+  let run = 0
+  dowhile run < mc_runs
+    set run = $&run             $ create a variable from the vector
+    let val = maxffts[{$run}]
+    let part = 0
+    * Check if val fits into a bin. If yes, raise bin by 1
+    dowhile part < bin_count
+      if ((val < (startfreq + (part+1)*bin_size)) & (val > (startfreq + part*bin_size)))
+        let yvec[part] = yvec[part] + 1
+                break
+      end
+      let part = part + 1
+    end
+    let run = run + 1
+  end
+  * plot the histogram
+  set plotstyle=combplot
+  let counts = yvec - 1             $ subtract 1 because we started with unitvec containing ones
+  plot counts vs osc_frequ
+* calculate jitter
+  let diff40 = (vecmax(halfffts) - vecmin(halfffts))*1e-6
+  echo
+  echo Max. jitter is "$&diff40" MHz
+end
+  rusage
+*  quit
+.endc
+
+.end