Resistor (2D) - Quasi Fermi Method

   ###################################################################
   #                          Constants                              #
   ###################################################################
   
   set T 300.0
   set k 1.38066e-23
   set q 1.619e-19
   set Vt [expr {$k*$T/$q}]
   set ni 1.1e10
   set esi [expr 11.8 * 8.85418e-14]
   set eps [expr $esi / $q]
   set Emob 350.0
   set Hmob 150.0
   set small 1.0e-10
   
   ###################################################################
   #                     Define the Grid                             #
   ###################################################################
   line x loc=0.0 spac=0.05 tag=Top
   line x loc=0.5 spac=0.1
   line x loc=1.0 spac=0.1 tag=Bottom
   
   line y loc=-2.0 spac=0.05 tag=Lside
   line y loc=2.0 spac=0.05  tag=Rside
   
   mater add name=Silicon
   region Silicon xlo=Top xhi=Bottom ylo=Lside yhi=Rside
   init
   
   ###################################################################
   #                  Contacts and Doping                            #
   ###################################################################
   contact name=GND Silicon xlo=-0.05 xhi=0.05 ylo=-1.25 yhi=-1.00 add
   contact name=VSS Silicon xlo=-0.05 xhi=0.05 ylo=1.00 yhi=1.25 add
   window col=2 row=1
   plot2d grid graph= Grid
   plot2d contact=GND graph= Grid
   plot2d contact=VSS graph= Grid
   
   sel z=1.0e20 name=Nd
   sel z=1.0e17 name=Na
   sel z=(Nd-Na) name=Doping
   
   ###################################################################
   #                       Equations                                 #
   ###################################################################
   # Bulk Equations
   DevicePackage
   solution name=Potential nosolve
   solution add name=DevPsi solve negative damp
   solution add name=Qfn solve negative 
   solution add name=Qfp solve negative
   solution add name=Elec const solve val= "1.0e10*exp((DevPsi-Qfn)/0.025)"
   solution add name=Hole const solve val= "1.0e10*exp((Qfp-DevPsi)/0.025)"
   
   set eqnP "$eps * grad(DevPsi) + Doping - Elec + Hole"
   set eqnQfn "ddt(Elec) - $Emob * Elec * grad(Qfn)"
   set eqnQfp "ddt(Hole) - $Hmob * Hole * grad(Qfp)"
   
   pdbSetString Silicon DevPsi Equation $eqnP
   pdbSetString Silicon Qfn    Equation $eqnQfn
   pdbSetString Silicon Qfp    Equation $eqnQfp
   
   pdbSetDouble Silicon DevPsi DampValue $Vt
   pdbSetDouble Silicon DevPsi Abs.Error 1.0e-9
   pdbSetDouble Silicon Elec   Abs.Error 1.0e-5
   pdbSetDouble Silicon Hole   Abs.Error 1.0e-5
   
   # Contact Equations
   proc ohmic.contact {Contact} {
       set vt 0.02558357
       set ni 1.1e10
       pdbSetBoolean $Contact Qfn Flux 1
       pdbSetBoolean $Contact Qfp Flux 1
       pdbSetBoolean $Contact DevPsi Flux 1
       pdbSetBoolean $Contact Qfn Fixed 1
       pdbSetBoolean $Contact Qfp Fixed 1
       pdbSetBoolean $Contact DevPsi Fixed 1
       pdbSetDouble $Contact Qfn Flux.Scale 1.619e-19
       pdbSetDouble $Contact Qfp Flux.Scale 1.619e-19
       pdbSetString $Contact Qfn Equation "Qfn-$Contact"
       pdbSetString $Contact Qfp Equation "Qfp-$Contact"
       pdbSetString $Contact DevPsi Equation "Doping - $ni*exp((DevPsi-Qfn)/$vt) + $ni*exp((Qfp-DevPsi)/$vt)"
       pdbSetString $Contact Equation "1.619e-19 * (Flux_Hole - Flux_Elec)"
   }
   ohmic.contact VSS
   ohmic.contact GND
   
   #Bias contacts
   contact name=VSS voltage supply=0.0
   contact name=GND voltage supply=0.0
   
   ###################################################################
   #                         Solve                                   #
   ###################################################################
   #Initial Guess
   sel z= {(Doping>0.0)
       ?  (0.025*log( (Doping + 1e-10) / 1.1e10))
       :  (-0.025*log(-(Doping + 1e-10) / 1.1e10))} name = DevPsi
   sel z=$ni*exp(DevPsi/$Vt)  name=Elec
   sel z=$ni*exp(-DevPsi/$Vt) name=Hole
   
   #First DC solve at Equilibrium (0V)
   device
   
   puts "Electron Flux [contact name=VSS sol=Qfn flux]"
   puts "Hole Flux [contact name=VSS sol=Qfp flux]"
   
   #initialize an array, bias (Vb) and counter (i), and create and new graph window
   array set curr {}
   set i 1
   #Ramp from 0.0-1.2V and Plot
   for {set bias 0.0} {$bias < 1.2} {set bias [expr $bias+0.05]} {
       contact name=VSS supply=$bias
       device init
       set curr($i) [expr ([contact name=VSS sol=Qfn flux] - [contact name=VSS sol=Qfp flux])]
       chart graph=IV curve=IV xval=$bias yval=$curr($i)
       incr i
   }