# 1.1 define system of equations
system	name=pairdiff
+	sysvars=interstitial,pipair,phosphorus,intcluster,inttot
+	nterm=7
+	term0=box_div.lapflux(interstitial|interstitial)@{silicon}
+	term1=nodal.cluster(intcluster,interstitial|intcluster,interstitial)@{silicon}
+	term2=interface.radiation(interstitial|interstitial)@{silicon/exposed}
+	term3=nodal.integrate(interstitial|inttot)@{silicon}
+	term4=nodal.two2one(phosphorus,interstitial,pipair|phosphorus,interstitial,pipair)@{silicon}
+	term5=box_div.lapflux(pipair|pipair)@{silicon}
+	term6=transient.ddt(interstitial,pipair,phosphorus,intcluster,inttot|interstitial,pipair,phosphorus,intcluster,inttot)@{silicon}

# 1.2 change default timestep
dbase modify name=library/math/systems/default_numerical_parameters/time.min rval=1e-12
dbase modify name=library/math/systems/default_numerical_parameters/time.init rval=1e-11

# 2.1 define intcluster parameters
dbase prefix=library/physics/silicon/intcluster
dbase create name=scale rval=1e10
dbase create name=background rval=1e10
dbase create name=kf.cluster.interstitial sval="4*PI*a*di"
dbase create name=kr.cluster.interstitial sval="di/(a*a)*exp(-(3.54-1.77)/kT)"
dbase prefix=""

# 2.2 define interstitial parameters
dbase prefix=library/physics/silicon/interstitial
dbase create name=scale rval=1e5
dbase create name=background rval=1e5
dbase create name=Dix sval="di=51*exp(-1.77/kT); 1e8*di"
dbase create name=Cstar sval="cs=2.9e24*exp(-3.18/kT)"
dbase create name=Krad sval="di/1e-7 *1e4"
dbase prefix=""

# 2.3 define inttot parameters
dbase create name = library/physics/silicon/inttot/scale rval=1

# 2.4 define pipair parameters
dbase prefix =library/physics/silicon/pipair
dbase create name=scale rval=1e5
dbase create name=background rval=1e-3
dbase create name=Dix
+	sval="dp=2.31e2/60*exp(-3.66/kT); dip=dp*(5e22/cs)*exp(-0.1/kT)/(4*PI); dip*1e8"
dbase create name=kf.two2one.phosphorus.interstitial
+	sval="4*PI*di*a"
dbase create name=kr.two2one.phosphorus.interstitial
+	sval="di/a^2*exp(-0.1/kT)"
dbase prefix=""

# 2.5 define phosphorus parameters
dbase prefix = library/physics/silicon/phosphorus
dbase create name=scale rval=1e10
dbase create name=background rval=1e10
dbase prefix=""


# 2.6 install constants
dbase create name=dummy sval="a=5e22^-(1/3)"
dbase print name=dummy temper=800

# 3.1 set up the structure and solve
grid xloc=0,1,10 xdel=0.01,0.01,1

implant elem=interstitial conc=1e8
implant elem=intcluster conc=1e12
implant elem=inttot conc=.001
implant elem=pipair conc=1e5
implant elem=phosphorus dose=5e13 range=1 sigma=0.02 plusfactor=1

dbase create name=options/timestep ival=1
dbase create name=options/dump.system ival=1
dbase create name=options/movie.reset ival=1
dbase create name=options/movie sval=interstitial,intcluster,phosphorus,pipair

dbase dump=mydbase
solve min=30 temper=1000 system=pairdiff

dump solution outf=runclust.out