summaryrefslogtreecommitdiff
path: root/Windows/spice/examples/Monte_Carlo
diff options
context:
space:
mode:
Diffstat (limited to 'Windows/spice/examples/Monte_Carlo')
-rw-r--r--Windows/spice/examples/Monte_Carlo/MC_2_circ.sp109
-rw-r--r--Windows/spice/examples/Monte_Carlo/MC_2_control.sp43
-rw-r--r--Windows/spice/examples/Monte_Carlo/MC_ring.sp251
-rw-r--r--Windows/spice/examples/Monte_Carlo/MonteCarlo.sp68
-rw-r--r--Windows/spice/examples/Monte_Carlo/OpWien.sp87
5 files changed, 558 insertions, 0 deletions
diff --git a/Windows/spice/examples/Monte_Carlo/MC_2_circ.sp b/Windows/spice/examples/Monte_Carlo/MC_2_circ.sp
new file mode 100644
index 00000000..bacd07b9
--- /dev/null
+++ b/Windows/spice/examples/Monte_Carlo/MC_2_circ.sp
@@ -0,0 +1,109 @@
+Circuit to perform Monte Carlo simulation in ngspice
+* 25 stage Ring-Osc. using inverters with BSIM3
+
+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
+
+.subckt inv1 dd ss sub well in out
+mn1 out in ss sub n1 w=2u l=0.35u as=3p ad=3p ps=4u pd=4u
+mp1 out in dd well p1 w=4u l=0.35u as=7p ad=7p ps=6u pd=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
+* output is buf
+cout buf ss 0.2pF
+*
+.options noacct
+
+* The following model parameters are varying statistically:
+* vth0, u0, tox
+* see the AGAUSS function used to define the parameter
+* the deviation is 10%, just for example, not measured
+
+********************************************************************************
+.model n1 nmos
++level=8
++version=3.3.0
++tnom=27.0
++nch=2.498e+17 tox=AGAUSS(9e-09, 9e-09, 10) xj=1.00000e-07
++lint=9.36e-8 wint=1.47e-7
++vth0=AGAUSS(.6322,.6322,10) k1=.756 k2=-3.83e-2 k3=-2.612
++dvt0=2.812 dvt1=0.462 dvt2=-9.17e-2
++nlx=3.52291e-08 w0=1.163e-6
++k3b=2.233
++vsat=86301.58 ua=6.47e-9 ub=4.23e-18 uc=-4.706281e-11
++rdsw=650 u0=AGAUSS(388.3203,388.3203,10) wr=1
++a0=.3496967 ags=.1 b0=0.546 b1=1
++dwg=-6.0e-09 dwb=-3.56e-09 prwb=-.213
++keta=-3.605872e-02 a1=2.778747e-02 a2=.9
++voff=-6.735529e-02 nfactor=1.139926 cit=1.622527e-04
++cdsc=-2.147181e-05
++cdscb=0 dvt0w=0 dvt1w=0 dvt2w=0
++cdscd=0 prwg=0
++eta0=1.0281729e-02 etab=-5.042203e-03
++dsub=.31871233
++pclm=1.114846 pdiblc1=2.45357e-03 pdiblc2=6.406289e-03
++drout=.31871233 pscbe1=5000000 pscbe2=5e-09 pdiblcb=-.234
++pvag=0 delta=0.01
++wl=0 ww=-1.420242e-09 wwl=0
++wln=0 wwn=.2613948 ll=1.300902e-10
++lw=0 lwl=0 lln=.316394 lwn=0
++kt1=-.3 kt2=-.051
++at=22400
++ute=-1.48
++ua1=3.31e-10 ub1=2.61e-19 uc1=-3.42e-10
++kt1l=0 prt=764.3
++noimod=2
++af=1.075e+00 kf=9.670e-28 ef=1.056e+00
++noia=1.130e+20 noib=7.530e+04 noic=-8.950e-13
+**** PMOS ***
+.model p1 pmos
++level=8
++version=3.3.0
++tnom=27.0
++nch=3.533024e+17 tox=AGAUSS(9e-09,9e-09,10) xj=1.00000e-07
++lint=6.23e-8 wint=1.22e-7
++vth0=AGAUSS(-.6732829,-.6732829,10) k1=.8362093 k2=-8.606622e-02 k3=1.82
++dvt0=1.903801 dvt1=.5333922 dvt2=-.1862677
++nlx=1.28e-8 w0=2.1e-6
++k3b=-0.24 prwg=-0.001 prwb=-0.323
++vsat=103503.2 ua=1.39995e-09 ub=1.e-19 uc=-2.73e-11
++rdsw=460 u0=AGAUSS(138.7609,138.7609,10)
++a0=.4716551 ags=0.12
++keta=-1.871516e-03 a1=.3417965 a2=0.83
++voff=-.074182 nfactor=1.54389 cit=-1.015667e-03
++cdsc=8.937517e-04
++cdscb=1.45e-4 cdscd=1.04e-4
++dvt0w=0.232 dvt1w=4.5e6 dvt2w=-0.0023
++eta0=6.024776e-02 etab=-4.64593e-03
++dsub=.23222404
++pclm=.989 pdiblc1=2.07418e-02 pdiblc2=1.33813e-3
++drout=.3222404 pscbe1=118000 pscbe2=1e-09
++pvag=0
++kt1=-0.25 kt2=-0.032 prt=64.5
++at=33000
++ute=-1.5
++ua1=4.312e-9 ub1=6.65e-19 uc1=0
++kt1l=0
++noimod=2
++af=9.970e-01 kf=2.080e-29 ef=1.015e+00
++noia=1.480e+18 noib=3.320e+03 noic=1.770e-13
+.end
+
+.end
diff --git a/Windows/spice/examples/Monte_Carlo/MC_2_control.sp b/Windows/spice/examples/Monte_Carlo/MC_2_control.sp
new file mode 100644
index 00000000..9b9f3606
--- /dev/null
+++ b/Windows/spice/examples/Monte_Carlo/MC_2_control.sp
@@ -0,0 +1,43 @@
+* Perform Monte Carlo simulation in ngspice
+* script for use with 25 stage Ring-Osc. BSIM3
+* circuit is in MC_2_circ.sp
+* edit 'set sourcepath' for your path to circuit file
+* start script by 'ngspice -o MC_2_control.log MC_2_control.sp'
+*
+.control
+ save buf $ we just need output vector buf, save memory by more than 10x
+ let mc_runs = 100 $ number of runs for monte carlo
+ let run = 1 $ number of the actual run
+
+* Where to find the circuit netlist file MC_2_circ.sp
+ set sourcepath = ( D:\Spice_general\ngspice\examples\Monte_Carlo )
+
+* create file for frequency information
+ echo Monte Carlo, frequency of R.O. > MC_frequ.log
+
+* run the simulation loop
+ dowhile run <= mc_runs
+ * without the reset switch there is some strange drift
+ * towards lower and lower frequencies
+ reset
+ set run ="$&run" $ create a variable from the vector
+ set rndseed = $run $ set the rnd seed value to the loop index
+ source MC_2_circ.sp $ load the circuit, including model data
+ tran 15p 200n 0
+ write mc_ring{$run}.out buf $ write each sim output to its own rawfile
+ linearize buf $ lienarize buf to allow fft
+ 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.1G to=0.7G
+ print fft_max >> MC_frequ.log $ print frequency to file
+ destroy all $ delete all output vectors
+ remcirc $ delete circuit
+ let run = run + 1 $ increase loop counter
+ end
+
+ quit
+
+.endc
+
+.end
diff --git a/Windows/spice/examples/Monte_Carlo/MC_ring.sp b/Windows/spice/examples/Monte_Carlo/MC_ring.sp
new file mode 100644
index 00000000..58e5c141
--- /dev/null
+++ b/Windows/spice/examples/Monte_Carlo/MC_ring.sp
@@ -0,0 +1,251 @@
+Perform Monte Carlo simulation in ngspice
+* 25 stage Ring-Osc. BSIM3
+
+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
+
+.subckt inv1 dd ss sub well in out
+mn1 out in ss sub n1 w=2u l=0.35u as=3p ad=3p ps=4u pd=4u
+mp1 out in dd well p1 w=4u l=0.35u as=7p ad=7p ps=6u pd=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
+*
+.options noacct
+.control
+ save buf $ we just need buf, save memory by more than 10x
+ 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
+*
+* define distributions for random numbers:
+* unif: uniform distribution, deviation relativ to nominal value
+* aunif: uniform distribution, deviation absolut
+* gauss: Gaussian distribution, deviation relativ to nominal value
+* agauss: Gaussian distribution, deviation absolut
+ define unif(nom, var) (nom + (nom*var) * sunif(0))
+ define aunif(nom, avar) (nom + avar * sunif(0))
+ define gauss(nom, var, sig) (nom + (nom*var)/sig * sgauss(0))
+ define agauss(nom, avar, sig) (nom + avar/sig * sgauss(0))
+*
+* We want to vary the model parameters vth0, u0, tox, lint, and wint
+* of the BSIM3 model for the NMOS and PMOS transistors.
+* We may obtain the nominal values (nom) by manually extracting them from
+* the parameter set. Here we get them automatically and store them into
+* variables. This has the advantage that you may change the parameter set
+* without having to look up the values again.
+ set n1vth0=@n1[vth0]
+ set n1u0=@n1[u0]
+ set n1tox=@n1[tox]
+ set n1lint=@n1[lint]
+ set n1wint=@n1[wint]
+ set p1vth0=@p1[vth0]
+ set p1u0=@p1[u0]
+ set p1tox=@p1[tox]
+ set p1lint=@p1[lint]
+ set p1wint=@p1[wint]
+*
+* run the simulation loop
+ dowhile run <= mc_runs
+ * without the reset switch there is some strange drift
+ * towards lower and lower frequencies
+ reset
+ * run=0 simulates with nominal parameters
+ if run > 0
+ altermod @n1[vth0]=gauss($n1vth0, 0.1, 3)
+ altermod @n1[u0]=gauss($n1u0, 0.05, 3)
+ altermod @n1[tox]=gauss($n1tox, 0.1, 3)
+ altermod @n1[lint]=gauss($n1lint, 0.1, 3)
+ altermod @n1[wint]=gauss($n1wint, 0.1, 3)
+ altermod @p1[vth0]=gauss($p1vth0, 0.1, 3)
+ altermod @p1[u0]=gauss($p1u0, 0.1, 3)
+ altermod @p1[tox]=gauss($p1tox, 0.1, 3)
+ altermod @p1[lint]=gauss($p1lint, 0.1, 3)
+ altermod @p1[wint]=gauss($p1wint, 0.1, 3)
+ end
+ tran 15p 50n 0
+* 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.
+*
+* We have to figure out what to do if a single simulation will not converge.
+* Is there a variable which may be set if there is no convergence?
+* Then we might skip this run and continue with a new run. It does not exist for now.
+* So we have to rely on the robustness of the following steps not leading
+* to a seg fault if the tran data are missing.
+*
+ set run ="$&run" $ create a variable from the vector
+ set mc_runs ="$&mc_runs" $ create a variable from the vector
+ echo simulation run no. $run of $mc_runs
+ * save the linearized data for having equal time scales for all runs
+ linearize buf $ linearize only buf, no other vectors needed
+ set dt = $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={$dt}.time
+ end
+ let vout{$run}={$dt}.buf $ store the output vector to plot 'plot_out'
+ setplot $dt $ 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.1G 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.1G to=0.7G
+ * store the fft vector
+ set dt = $curplot $ store the current plot to dt (spec i)
+ setplot $plot_fft $ make 'plot_fft' the active plot
+ if run=0
+ let frequency={$dt}.frequency
+ end
+ let fft{$run}={$dt}.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}]={$dt}.fft_max
+ let halfffts[{$run}]={$dt}.fft_40
+* setplot $plot_out
+* The following command does not work here. Why not? Probably not a real copy.
+* destroy $dt $ save memory, we don't need this plot (spec) any more
+ setplot $dt $ go back to the previous plot
+ let run = run + 1
+ end
+***** plotting **********************************************************
+* plot {$plot_out}.allv
+ plot {$plot_out}.vout0 $ just plot the tran output with nominal parameters
+* setplot $plot_fft
+* plot db(mag(ally)) xlimit .1G 1G ylimit -80 10
+ plot db(mag({$plot_fft}.ally)) xlimit .1G 1G ylimit -80 10
+*
+* create a histogram from vector maxffts
+ setplot $max_fft $ make 'max_fft' the active plot
+ set startfreq=400MEG
+ set bin_size=5MEG
+ set bin_count=20
+ compose xvec start=$startfreq step=$bin_size lin=$bin_count $ requires variables as parameters
+ settype frequency xvec
+ 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
+ plot yvec-1 vs xvec $ subtract 1 because with started with unitvec containing ones
+* calculate jitter
+ let diff40 = (vecmax(halfffts) - vecmin(halfffts))*1e-6
+ echo
+ echo Max. jitter is "$&diff40" MHz
+ rusage
+.endc
+********************************************************************************
+.model n1 nmos
++level=8
++version=3.3.0
++tnom=27.0
++nch=2.498e+17 tox=9e-09 xj=1.00000e-07
++lint=9.36e-8 wint=1.47e-7
++vth0=.6322 k1=.756 k2=-3.83e-2 k3=-2.612
++dvt0=2.812 dvt1=0.462 dvt2=-9.17e-2
++nlx=3.52291e-08 w0=1.163e-6
++k3b=2.233
++vsat=86301.58 ua=6.47e-9 ub=4.23e-18 uc=-4.706281e-11
++rdsw=650 u0=388.3203 wr=1
++a0=.3496967 ags=.1 b0=0.546 b1=1
++dwg=-6.0e-09 dwb=-3.56e-09 prwb=-.213
++keta=-3.605872e-02 a1=2.778747e-02 a2=.9
++voff=-6.735529e-02 nfactor=1.139926 cit=1.622527e-04
++cdsc=-2.147181e-05
++cdscb=0 dvt0w=0 dvt1w=0 dvt2w=0
++cdscd=0 prwg=0
++eta0=1.0281729e-02 etab=-5.042203e-03
++dsub=.31871233
++pclm=1.114846 pdiblc1=2.45357e-03 pdiblc2=6.406289e-03
++drout=.31871233 pscbe1=5000000 pscbe2=5e-09 pdiblcb=-.234
++pvag=0 delta=0.01
++wl=0 ww=-1.420242e-09 wwl=0
++wln=0 wwn=.2613948 ll=1.300902e-10
++lw=0 lwl=0 lln=.316394 lwn=0
++kt1=-.3 kt2=-.051
++at=22400
++ute=-1.48
++ua1=3.31e-10 ub1=2.61e-19 uc1=-3.42e-10
++kt1l=0 prt=764.3
++noimod=2
++af=1.075e+00 kf=9.670e-28 ef=1.056e+00
++noia=1.130e+20 noib=7.530e+04 noic=-8.950e-13
+**** PMOS ***
+.model p1 pmos
++level=8
++version=3.3.0
++tnom=27.0
++nch=3.533024e+17 tox=9e-09 xj=1.00000e-07
++lint=6.23e-8 wint=1.22e-7
++vth0=-.6732829 k1=.8362093 k2=-8.606622e-02 k3=1.82
++dvt0=1.903801 dvt1=.5333922 dvt2=-.1862677
++nlx=1.28e-8 w0=2.1e-6
++k3b=-0.24 prwg=-0.001 prwb=-0.323
++vsat=103503.2 ua=1.39995e-09 ub=1.e-19 uc=-2.73e-11
++rdsw=460 u0=138.7609
++a0=.4716551 ags=0.12
++keta=-1.871516e-03 a1=.3417965 a2=0.83
++voff=-.074182 nfactor=1.54389 cit=-1.015667e-03
++cdsc=8.937517e-04
++cdscb=1.45e-4 cdscd=1.04e-4
++dvt0w=0.232 dvt1w=4.5e6 dvt2w=-0.0023
++eta0=6.024776e-02 etab=-4.64593e-03
++dsub=.23222404
++pclm=.989 pdiblc1=2.07418e-02 pdiblc2=1.33813e-3
++drout=.3222404 pscbe1=118000 pscbe2=1e-09
++pvag=0
++kt1=-0.25 kt2=-0.032 prt=64.5
++at=33000
++ute=-1.5
++ua1=4.312e-9 ub1=6.65e-19 uc1=0
++kt1l=0
++noimod=2
++af=9.970e-01 kf=2.080e-29 ef=1.015e+00
++noia=1.480e+18 noib=3.320e+03 noic=1.770e-13
+.end
diff --git a/Windows/spice/examples/Monte_Carlo/MonteCarlo.sp b/Windows/spice/examples/Monte_Carlo/MonteCarlo.sp
new file mode 100644
index 00000000..608dc4cf
--- /dev/null
+++ b/Windows/spice/examples/Monte_Carlo/MonteCarlo.sp
@@ -0,0 +1,68 @@
+* Effecting a Monte Carlo calculation in ngspice
+V1 N001 0 AC 1 DC 0
+R1 N002 N001 141
+*
+C1 OUT 0 1e-09
+L1 OUT 0 10e-06
+C2 N002 0 1e-09
+L2 N002 0 10e-06
+L3 N003 N002 40e-06
+C3 OUT N003 250e-12
+*
+R2 0 OUT 141
+.control
+ let mc_runs = 5
+ let run = 0
+ set curplot=new $ create a new plot
+ set scratch=$curplot $ store its name to 'scratch'
+ setplot $scratch $ make 'scratch' the active plot
+ let bwh=unitvec(mc_runs) $ create a vector in plot 'scratch' to store bandwidth data
+
+* define distributions for random numbers:
+* unif: uniform distribution, deviation relativ to nominal value
+* aunif: uniform distribution, deviation absolut
+* gauss: Gaussian distribution, deviation relativ to nominal value
+* agauss: Gaussian distribution, deviation absolut
+* limit: if unif. distributed value >=0 then add +avar to nom, else -avar
+ define unif(nom, rvar) (nom + (nom*rvar) * sunif(0))
+ define aunif(nom, avar) (nom + avar * sunif(0))
+ define gauss(nom, rvar, sig) (nom + (nom*rvar)/sig * sgauss(0))
+ define agauss(nom, avar, sig) (nom + avar/sig * sgauss(0))
+* define limit(nom, avar) (nom + ((sgauss(0) ge 0) ? avar : -avar))
+ define limit(nom, avar) (nom + ((sgauss(0) >= 0) ? avar : -avar))
+*
+*
+ dowhile run < mc_runs $ loop starts here
+*
+* alter c1 = unif(1e-09, 0.1)
+* alter c1 = aunif(1e-09, 100e-12)
+* alter c1 = gauss(1e-09, 0.1, 3)
+* alter c1 = agauss(1e-09, 100e-12, 3)
+*
+ alter c1 = unif(1e-09, 0.1)
+ alter l1 = unif(10e-06, 0.1)
+ alter c2 = unif(1e-09, 0.1)
+ alter l2 = unif(10e-06, 0.1)
+ alter l3 = unif(40e-06, 0.1)
+ alter c3 = limit(250e-12, 25e-12)
+*
+ ac oct 100 250K 10Meg
+*
+* measure bandwidth at -10 dB
+ meas ac bw trig vdb(out) val=-10 rise=1 targ vdb(out) val=-10 fall=1
+*
+ set run ="$&run" $ create a variable from the vector
+ set dt = $curplot $ store the current plot to dt
+ setplot $scratch $ make 'scratch' the active plot
+ let vout{$run}={$dt}.v(out) $ store the output vector to plot 'scratch'
+ let bwh[run]={$dt}.bw $ store bw to vector bwh in plot 'scratch'
+ setplot $dt $ go back to the previous plot
+ let run = run + 1
+ end $ loop ends here
+*
+ plot db({$scratch}.allv)
+ echo
+ print {$scratch}.bwh
+.endc
+
+.end
diff --git a/Windows/spice/examples/Monte_Carlo/OpWien.sp b/Windows/spice/examples/Monte_Carlo/OpWien.sp
new file mode 100644
index 00000000..91b45d08
--- /dev/null
+++ b/Windows/spice/examples/Monte_Carlo/OpWien.sp
@@ -0,0 +1,87 @@
+OPWIEN.CIR - OPAMP WIEN-BRIDGE OSCILLATOR
+* http://www.ecircuitcenter.com/circuits/opwien/opwien.htm
+* single simulation run
+* 2 resistors and 2 capacitors of Wien bridge a varied statistically
+* number of variations: varia
+
+* Simulation time
+.param ttime=12000m
+.param varia=100
+.param ttime10 = 'ttime/varia'
+
+* nominal resistor and capacitor values
+.param res = 10k
+.param cn = 16NF
+
+* CURRENT PULSE TO START OSCILLATIONS
+IS 0 3 dc 0 PWL(0US 0MA 10US 0.1MA 40US 0.1MA 50US 0MA 10MS 0MA)
+*
+* RC TUNING
+VR2 r2 0 dc 0 trrandom (2 'ttime10' 0 1) $ Gauss controlling voltage
+*
+*VR2 r2 0 dc 0 trrandom (1 'ttime10' 0 3) $ Uniform within -3 3
+*
+* If Gauss, factor 0.033 is 10% equivalent to 3 sigma
+* if uniform, uniform between +/- 10%
+R2 4 6 R = 'res + 0.033 * res*V(r2)' $ behavioral resistor
+*R2 4 6 'res' $ constant R
+
+VC2 c2 0 dc 0 trrandom (2 'ttime10' 0 1)
+*C2 6 3'cn' $ constant C
+C2 6 3 C = 'cn + 0.033 * cn*V(c2)' $ behavioral capacitor
+
+VR1 r1 0 dc 0 trrandom (2 'ttime10' 0 1)
+*VR1 r1 0 dc 0 trrandom (1 'ttime10' 0 3)
+R1 3 0 R = 'res + 0.033 * res*V(r1)'
+*R1 3 0 'res'
+
+VC1 c1 0 dc 0 trrandom (2 'ttime10' 0 1)
+C1 3 0 C = 'cn + 0.033 * cn*V(c2)'
+*C1 3 0 'cn'
+
+* NON-INVERTING OPAMP
+R10 0 2 10K
+R11 2 5 18K
+XOP 3 2 4 OPAMP1
+* AMPLITUDE STABILIZATION
+R12 5 4 5K
+D1 5 4 D1N914
+D2 4 5 D1N914
+*
+.model D1N914 D(Is=0.1p Rs=16 CJO=2p Tt=12n Bv=100 Ibv=0.4n)
+*
+* OPAMP MACRO MODEL, SINGLE-POLE
+* connections: non-inverting input
+* | inverting input
+* | | output
+* | | |
+.SUBCKT OPAMP1 1 2 6
+* INPUT IMPEDANCE
+RIN 1 2 10MEG
+* DC GAIN (100K) AND POLE 1 (100HZ)
+EGAIN 3 0 1 2 100K
+RP1 3 4 1K
+CP1 4 0 1.5915UF
+* OUTPUT BUFFER AND RESISTANCE
+EBUFFER 5 0 4 0 1
+ROUT 5 6 10
+.ENDS
+*
+* ANALYSIS
+.TRAN 0.05MS 'ttime'
+*
+* VIEW RESULTS
+.control
+option noinit
+run
+plot V(4) 5*V(r1) 5*V(r2) 5*V(c1) 5*V(c2)
+linearize v(4)
+fft v(4)
+let v4mag = mag(v(4))
+plot v4mag
+plot v4mag xlimit 500 1500
+*wrdata histo v4mag
+rusage
+.endc
+
+.END