# Let's define the temperature as well as k*T expressed in electron volts
math/def T 300
math/def kT kB*T/e

# Effecive electron and hole masses for silicon
math/def m_e 1.08*me
math/def m_h 0.56*me

# Band gap EG is for silicon, in electron volts
# EA is the energy level of acceptor states, relative above the valence
# band. ED is the energy level of the donor states, expressed as an amount
# below the conduction band
math/def EG 1.11
math/def EA 0.045
math/def ED 0.045

# From the effective masses, we can calculate the effecive number of states
# in the conduction and valence bands
math/def NC 2*(2*PI*m_e*kT*e/(h*h))^(3/2)
math/def NV 2*(2*PI*m_h*kT*e/(h*h))^(3/2)

# Let's specify the density of acceptors and donors at the n-side of the
# device.
math/def NA_n 0
math/def ND_n 1e22

# For a block of n-type material, we can calculate the chemical potential
# by solving the equation for charge neutrality:
#
#  n + N_A(-) = p + N_D(+)
#
# where
#
#  n = N_C*exp((mu-E_G)/(kT))
#  p = N_V*exp(-mu/(kT))
#  N_A(-) = N_A*f(E_A)
#  N_D(+) = N_D*(1-f(E_G-E_D))
#
# and
#
#  f(E) = 1/(exp((E-mu)/(kT)) + 1)
#
# is the Fermi distribution
#
# The syntax of this command is:
#
#   math/calc/mu E_min E_max E_A E_D E_gap kT N_C N_V N_A N_D variable
#
# where the units used don't really matter, as long as they are consistent.
# This commands solves the neutrality equation in terms of the chemical
# potential. It is assumed to lie in [E_min,E_max] (here [0,EG]), and
# the rest of the parameters have already been encountered. When the
# chemical potential could be found, it's value is stored in the specified
# variable, here 'mu_n'
#
math/calc/mu 0 EG EA ED EG kT NC NV NA_n ND_n mu_n

# We also calculate the chemical potential for the p-type material
math/def NA_p 0.5e22
math/def ND_p 0
math/calc/mu 0 EG EA ED EG kT NC NV NA_p ND_p mu_p

# When the n- and p-type materials are brought into contact, charge
# redistribution will occur in the contact region, and their chemical
# potentials will align. The built-in voltage that arises is precisely
# the difference between the chemical potentials.
math/def Vint mu_n-mu_p

# Let's define a variable for the applied voltage
math/def Vapp 0

# Defining some variables based on the p-side settings. The first two
# will just be used to initialize densities, the last two are used
# to set the background charge in the Poisson equation. 
math/def p_p NV*exp(-mu_p/kT)
math/def n_p NC*exp((mu_p-EG)/kT)
math/def NAmin_p NA_p/(exp((EA-mu_p)/kT)+1)
math/def NDplus_p ND_p/(1+exp((mu_p-EG+ED)/kT))

# Same for the n-side settings
math/def p_n NV*exp(-mu_n/kT)
math/def n_n NC*exp((mu_n-EG)/kT)
math/def NAmin_n NA_n/(exp((EA-mu_n)/kT)+1)
math/def NDplus_n ND_n/(1+exp((mu_n-EG+ED)/kT))

# Defining variables for mobilities, diffusion constants (using Einstein
# relation) and of relative permittivity
math/def cm 0.01
math/def eMob 1200.0*cm*cm
math/def hMob 400.0*cm*cm
math/def De eMob*kT
math/def Dh hMob*kT
math/def epsilonRel 11.7

# We'll be using generation/recombination of type G-R = G-r*(n*p-n_i^2)
# In case no current flows (dark situation), the final drift-diffusion
# equations will become 
#
#   0 = G - r*(n*p - n_i^2) - 0
#
# so in the dark (G = 0), we'd find n_i = sqrt(n*p). This will differ
# on each side of the device, and well define variables containing these
# values for later use.
math/def ni_p sqrt(p_p*n_p)
math/def ni_n sqrt(p_n*n_n)
# We'll use this for the generation of electron-hole pairs
math/def GEN 1e27

# Now that we've defined a whole lot of useful variables, we're going
# to create an actual 1D simulation. Here, our device will we 2 micrometers
# long, and the equations will be discretized in 256 points.
sim1/new 256 2e-6

# We'll define two region names: pside, for the left part of the device,
# and nside for the right part of the device.
reg1/new pside
reg1/append/line pside 1 NX1/2

reg1/new nside
reg1/append/line nside NX1/2+1 NX1

# Initialize p, n and the background charge at the p- and n-sides of the 
# material. Note that the p and n values at the contacts (at grid point 0 
# and NX1) will also be set this way (and these won't change during the
# simulation)
sim1/reg/set p pside p_p
sim1/reg/set n pside n_p
sim1/reg/set bg pside NDplus_p-NAmin_p

sim1/reg/set p nside p_n
sim1/reg/set n nside n_n
sim1/reg/set bg nside NDplus_n-NAmin_n

# Set generation, intrinsic electron density and recombination prefactors.
# For this last one, we'll use the Lagevin prefactor for bimolecular
# recombination.
sim1/reg/set g ALL GEN
sim1/reg/set ni pside ni_p
sim1/reg/set ni nside ni_n
sim1/reg/set rf ALL e*(eMob+hMob)/(eps0*epsilonRel)

# Set mobilities, diffusion constants and relative permittivity on all
# grid points
sim1/reg/set nmob ALL eMob
sim1/reg/set pmob ALL hMob
sim1/reg/set dn ALL De
sim1/reg/set dp ALL Dh
sim1/reg/set eps ALL epsilonRel

# Set the voltage drop over the device as being the difference between
# the internal voltage and the applied voltage, as usual. We'll initialize
# the electrostatic potential in the interior of the device by linearly
# interpolating the boundaries. We'll also store the temperature (not used
# by the equations themselves, but can be used to determine some scale
# factors)
sim1/phi/set Vint-Vapp yes
sim1/temp/set T

# We're going to 'run' the simulation for 1 step, with dt=0, this is just
# so that the currents have been calculated throughout the device, which
# we can then plot
sim1/run 1 0

# Save this starting situation to a file
sim1/save pn1d_start.dat

sim1/grid/wplot n
sim1/grid/wplot p
sim1/grid/wplot V
sim1/grid/wplot jx