##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
##PNJunctionAbrupt.py
##
##PURPOSE:
##  To serve as a starting point for calculating the field and potential
##  in an abrupt PN junction, i.e one in which the donor and acceptor 
##  concentrations are described by step functions
##
##KEYWORDS:
##  Mode: int
##    0 - use a Newta Newton Iteration methodd
##    1 - use an implicit method
##  TvFlag: int
##    0 - Don't plot
##    1 - Plot according to TvPPT
##  TvPPT: int
##    0 - Plot every step
##    1 - Plot every 200th step
##  NSteps: int
##    The number of steps to execute
##  SaveFig: int
##    0 - Don't save any figures
##    1 - Save the figures at the end
##  V_bias: float
##    The bias voltage applied to the n side of the junction
##
##EXAMPLES
##  run PNJunctionAbrupt.py
##
###################################################################
import matplotlib.pylab as mplot
from numpy import *
from pylab import *
from LU_decomposition import *
from FDGrid_1d import *
import pdb

execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
PlotDir = '/nfs/slac/g/ki/ki03/lances/PNSimulations/'

##KEYWORDS
keywords    = ['Mode=', 'TvFlag=', 'NSteps=', 'SaveFig=','SaveMultFigs=',\
               'TvPPT=', 'V_bias=', 'Verbosity=']
Mode = 0
TvFlag  = 0
NSteps  = 85
SaveMultFigs=0
SaveFig = 0
TvPPT   = 0
V_bias  = 0
Verbosity = 0

opts, extraparams = getopt.getopt(sys.argv[1:],'',keywords)
for o,p in opts:
  if o in ['--Mode']:
    Mode = int(p)
  elif o in ['--TvFlag']:
    TvFlag = int(p)
  elif o in ['--NSteps']:
    NSteps = int(p)
  elif o in ['--SaveFig']:
    SaveFig= int(p)
  elif o in ['--SaveMultFigs']:
    SaveMultFigs=int(p)
  elif o in ['--TvPPT']:
    TvPPT = int(p)
  elif o in ['--V_bias']:
    V_bias = float(p)
  elif o in ['--Verbosity']:
    Verbosity = int(p)

##MODE STRING
if Mode == 0:
   ModeStr = 'Newton'
else:
   ModeStr = 'Implicit'

##############################################################################
##CONSTANTS:
e_0   = 8.85419e-12          # F * m^-1 - 8.85419 e-14 F*cm^-1
e_si  = 11.8                 # Relative Permittivity for Si	
q     = double(1.6022e-19)   # Coulombs per electron
k     = double(1.38066e-23)  # Joules per Kelvin
T     = double(300)          # Kelvin
Theta = double(q/(k*T))      # Factor relating energies

##Effective Masses
m_e  = 9.1e-31             # Kg
m_ee = 0.26                # effective mass of electrons/mass of electron
m_eh = 0.386               # effective mass of holes/mass of electron

##Mobilities
mu_n = double(1200.e-4)   # m^2/V s 
mu_p = double(400.e-4)    # m^2/V s

##Diffusion constants
D_n  = k*T*mu_n/q
D_p  = k*T*mu_p/q

##############################################################################
##Parameters for the PIN didoe - 100 um thick
n_i     = 1.e16
n_bulk  = 1.e18
N_d_max = 1.e24
N_a_max = 1.e24

##The changes in dx between i and i +/- 1/2
TP   = 300
GR   = FDGrid_1d(TP, N_a_max, n_bulk, N_d_max, 1, Simple = 1)
X    = GR.X;    dx = GR.dx;     dx_p = GR.dx_p 
N    = 4*TP-1;  NN = N-1;       MN   = N/2

##Difference is doping density
DC = GR.DC
#########################################################################

##Assume Phi at x=-inf is 0, calculate Phi_max according to 
##Ashcroft Mermin p. 575
E_g = 1.12*q               # Joules = eV*q 
N_c  = 2.5*(m_ee)**1.5*(T/300)**1.5*10.**18
P_v  = 2.5*(m_eh)**1.5*(T/300)**1.5*10.**18

##Built-in-potential calculated in Ashcroft-Mermin
phi_Max = (E_g+k*T*log(N_d_max*N_a_max/(N_c*P_v)))/q

##Built-in potential calculated in Neudeck 
V_bi   = k*T*log(N_d_max*N_a_max/n_i**2)/q

##Calculation of depletion region width
x_p     = ((2*e_si*e_0*V_bi/q)*(N_d_max/(N_a_max*(N_a_max+N_d_max))))**.5
x_n     = ((2*e_si*e_0*V_bi/q)*(N_a_max/(N_d_max*(N_a_max+N_d_max))))**.5
E_Max   = (q*N_d_max/(e_si*e_0))*x_n

##Now setup and integrate
phi   = zeros(N, dtype=double)
E     = zeros(N, dtype=double)
E_bar = zeros(N, dtype = double)
E_del = zeros(N, dtype = double)
n     = GR.n
p     = GR.p

##Boundaries potentails calculated from Kurata p. 36
phi_Min = -(1/Theta)*log(p[0]/n_i)
phi_Max =  (1/Theta)*log(n[NN]/n_i)+V_bias
phi[0]  = phi_Min
phi[NN] = phi_Max
phi[:]  = frange(NN)*(phi_Max-phi_Min)/NN + phi_Min
y       = zeros([N,3],dtype=double)
y[:,0]  = p; y[:,1] = n; y[:,2] = phi

##########################################################
dt      = 1.e-22        #Time step

#############################################################################
for Sn in arange(NSteps):

  ##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ##Calculate the density, currents, etc.
  ##y[:,0] = p, y[:,1] = n, y[:,2] = phi
  rho  = q*(y[:,0] - y[:,1] + DC)
  Be = Theta*(y[0:NN,2]-y[1:N,2])
  phi_del = y[1:N,2]-y[0:NN,2]
  lambda_p1 = mu_p*(-phi_del)/(1-exp(-Be))
  lambda_p2 = mu_p*(-phi_del)/(1-exp(Be))
  lambda_n1 = mu_n*(-phi_del)/(1-exp(Be))
  lambda_n2 = mu_n*(-phi_del)/(1-exp(-Be))
  Jp = q*(lambda_p1*y[0:NN,0]+lambda_p2*y[1:N,0])/dx
  Jn = q*(lambda_n1*y[0:NN,1]+lambda_n2*y[1:N,1])/dx

  ##++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ##Twelve Matrix elements
  ##P Terms
  MA1  = q*lambda_p1/dx                                    # dJp_i+1/2/dp_i
  MA2  = q*lambda_p2/dx                                    # dJp_i+1/2/dp_i+1
  MA3  = mu_p*(1-exp(-Be)-Be*exp(-Be))/((1-exp(-Be))**2)   # dLam_p_1/d_phi_i
  MA4  = mu_p*(1-exp(Be) +Be*exp(Be) )/((1-exp(Be))**2)    # dLam_p_2/d_phi_i
  MA5  = q*(y[0:NN,0]*MA3+y[1:N,0]*MA4)/dx                 # dJp_i+1/2/d_phi_i
  MA6  = -MA5                                              # dJp_i+1/2/d_phi_i+1
  ##N Terms
  MA7  = q*lambda_n1/dx                                    # dJn_i+1/2/dn_i
  MA8  = q*lambda_n2/dx                                    # dJn_i+1/2/dn_i+1
  MA9  = (mu_n/mu_p)*MA4                                   # dLam_n_1/d_phi_i
  MA10 = (mu_n/mu_p)*MA3                                   # dLam_n_2/d_phi_i
  MA11 = q*(y[0:NN,1]*MA9+y[1:N,1]*MA10)/dx                # dJn_i+1/2/d_phi_i
  MA12 = -MA11                                             # dJn_i+1/2/d_phi_i+1

  ##Form the matrix elements
  ##A*del_y_i-1 + B*del_y + C*del_y_i+1 = F_i
  ##A Matrix####################
  A    = zeros([NN,3,3],dtype=double)
  A[0:NN-1,0,0]  = -MA1[0:NN-1]/(q*dx_p)   
  A[0:NN-1,0,1]  = 0   
  A[0:NN-1,0,2]  = -MA5[0:NN-1]/(q*dx_p)

  A[0:NN-1,1,0] = 0                 
  A[0:NN-1,1,1] = -MA7[0:NN-1]/(q*dx_p) 
  A[0:NN-1,1,2] = -MA11[0:NN-1]/(q*dx_p)

  A[0:NN-1,2,0] = 0             
  A[0:NN-1,2,1] = 0         
  A[0:NN-1,2,2] = 1/(dx[0:NN-1]*dx_p)

  ##B Matrix#####################
  if Mode == 1:
    ImpBool = 1
  else:
    ImpBool = 0
  del_t = 1.e-12
  B    = zeros([N,3,3],dtype=double)
  B[1:NN,0,0]  = (MA1[1:NN]-MA2[0:NN-1])/(q*dx_p) +1./del_t*ImpBool 
  B[1:NN,0,1]  = 0  
  B[1:NN,0,2]  = (MA5[1:NN]-MA6[0:NN-1])/(q*dx_p)

  B[1:NN,1,0]  = 0                      
  B[1:NN,1,1]  = (MA7[1:NN]-MA8[0:NN-1])/(q*dx_p) -1./del_t*ImpBool
  B[1:NN,1,2]  = (MA11[1:NN]-MA12[0:NN-1])/(q*dx_p)
  
  B[1:NN,2,0]  = q/(e_si*e_0)
  B[1:NN,2,1]  = -q/(e_si*e_0)
  B[1:NN,2,2]  = (-1/dx[0:NN-1] - 1/dx[1:NN])*(1/dx_p)

  #Boundary Conditions are in B[0,:,:] and B[N,:,:]
  B[0,:,:]  = eye(3,dtype=double)   # del_y = 0
  B[NN,:,:] = eye(3,dtype=double)   # del_y = 0

  ##C Matrix######################
  C    = zeros([NN,3,3],dtype=double)
  C[1:NN,0,0] = MA2[1:NN]/(q*dx_p)        
  C[1:NN,0,1] = 0           
  C[1:NN,0,2] = MA6[1:NN]/(q*dx_p)

  C[1:NN,1,0] = 0             
  C[1:NN,1,1] = MA8[1:NN]/(q*dx_p)      
  C[1:NN,1,2] = MA12[1:NN]/(q*dx_p)

  C[1:NN,2,0] = 0             
  C[1:NN,2,1] = 0           
  C[1:NN,2,2] = 1/(dx[1:NN]*dx_p)

  ##F Matrix######################
  F = zeros([N,3],dtype=double)
  F[1:NN,0]  = -(Jp[1:NN]-Jp[0:NN-1])/(dx_p*q)
  F[1:NN,1]  = -(Jn[1:NN]-Jn[0:NN-1])/(dx_p*q)
  F[1:NN,2]  = -rho[1:NN]/(e_si*e_0)\
               -1/(dx[0:NN-1]*dx_p)*y[0:NN-1,2]\
               +1/dx_p*(1/dx[0:NN-1]+1/dx[1:NN])*y[1:NN,2]\
               -1/(dx[1:NN]*dx_p)*y[2:N,2]
  d_y   = LU_solve_3x3_1d(0, NN, A, B, C, F) 
  y = y + d_y

  ##Calculate the electric field
  E = (y[1:N,2]-y[0:NN,2])/dx

  ##Find the residual
  p_res = d_y[:,0].max() ;p_res_loc = (where(d_y[:,0] == p_res)[0])[0]
  n_res = d_y[:,1].max() ;n_res_loc = (where(d_y[:,1] == n_res)[0])[0]
  phi_res = d_y[:,2].max() ; phi_res_loc = (where(d_y[:,2] == phi_res)[0])[0]
  if Verbosity == 1:
    print 'Step Number: ' +str(Sn)
    print 'p residual : ' +str(p_res) + ' at ' + str(p_res_loc)
    print 'n residual : ' +str(n_res) + ' at ' + str(n_res_loc)
    print 'phi residual : ' + str(phi_res) + ' at ' + str(phi_res_loc)
 
  ##Plot Every Step
  if TvFlag == 1 and TvPPT == 0:
    mplot.figure(0)
    mplot.clf()
    mplot.semilogy(X, y[:,0])
    mplot.semilogy(X, y[:,1])
    mplot.title('Number of Holes/Electrons')
    mplot.figure(1)
    mplot.clf()
    mplot.plot(X, rho)
    mplot.title('Charge Density')
    mplot.figure(2)
    mplot.clf()
    mplot.plot(X, y[:,2])
    mplot.figure(3)
    mplot.clf()
    mplot.plot(X[1:N], E)

PlotLast = 1
if PlotLast == 1:

  #Get the full range
  XFullRange = X.max()-X.min()
  XJunction  = XFullRange/2

  #Plot the hole/electron densities
  mplot.figure(0)
  mplot.clf()
  mplot.semilogy(X, y[:,0]*1e-6, label = '$p_{v}$', lw=2.5, color='r')
  mplot.semilogy(X, y[:,1]*1e-6, label = '$n_{c}$', lw=2.5, color='b')
  mplot.axvline(X[N/2]-x_p, color='k', linestyle='--')
  mplot.axvline(X[N/2]+x_n, color='k', linestyle='--')
  mplot.xlabel('x (meters)', fontsize=18)
  mplot.ylabel('Number Density (cm$^{-3}$)', fontsize=18)
  YRange = mplot.ylim()
  YMin   = YRange[0] ; YMax = YRange[1]
  mplot.text(XJunction-x_p-.05*XFullRange, sqrt(YRange[1]), r'$x_{p}$', \
             fontsize=16,color='k')
  mplot.text(XJunction+x_n+.02*XFullRange, sqrt(YRange[1]), r'$x_{n}$', \
             fontsize=16,color='k')
  mplot.legend()
  #Fill in the P and N regions
  mplot.fill([X[0],X[N/2],X[N/2],X[0]], [YMin,YMin,YMax,YMax],\
              facecolor="#ffdede")
  mplot.fill([X[N/2],X[NN],X[NN],X[N/2]], [YMin,YMin,YMax,YMax],\
              facecolor="#dedeff")
  Ax = gca()
  Ax.set_ylim([YMin, YMax])
  mplot.draw()

  if SaveFig == 1:
    mplot.savefig(PlotDir+'PNJunction_'+ModeStr+'_Potential_'+\
                           str(NSteps)+'Steps.png')
    mplot.savefig(PlotDir+'PNJunction_'+ModeStr+'_Potential_'+\
                           str(NSteps)+'Steps.eps') 

  #Plot the charge density
  mplot.figure(1)
  mplot.clf()
  mplot.plot(X, rho*1.e-6, lw=2.5, color='k')
  mplot.axvline(X[N/2]-x_p, color='k', linestyle='--')
  mplot.axvline(X[N/2]+x_n, color='k', linestyle='--')
  mplot.xlabel('x (meters)', fontsize=18)
  mplot.ylabel(r'Space Charge Density (C\ cm$^{-3}$)',
               fontsize=18)
  #Here is the proper text for the charge density 
  #r'$\left[\rho = q(p_{v}-n_{c}+N_{d}^{+}-N_{a}^{-})\right]$'
  YRange = mplot.ylim()
  YMin   = YRange[0] ; YMax = YRange[1]
  mplot.text(XJunction-x_p-.05*XFullRange,  (YMin+YMax)/2, r'$x_{p}$', \
             fontsize=16,color='k')
  mplot.text(XJunction+x_n+.02*XFullRange,  (YMin+YMax)/2, r'$x_{n}$', \
             fontsize=16,color='k')
  #Fill in the P and N regions
  YRange = mplot.ylim()
  YMin   = YRange[0] ; YMax = YRange[1]
  mplot.fill([X[0],X[N/2],X[N/2],X[0]], [YMin,YMin,YMax,YMax],\
              facecolor="#ffdede")
  mplot.fill([X[N/2],X[NN],X[NN],X[N/2]], [YMin,YMin,YMax,YMax],\
              facecolor="#dedeff")
  Ax = gca()
  Ax.set_ylim([YMin, YMax])
  mplot.draw()

  if SaveFig == 1:
    mplot.savefig(PlotDir+'PNJunction_'+ModeStr+'_pv_nc_'+\
                  str(NSteps)+'Steps.png')
    mplot.savefig(PlotDir+'PNJunction_'+ModeStr+'_pv_nc_'+\
                  str(NSteps)+'Steps.eps') 

  #Plot the potential
  mplot.figure(2)
  mplot.clf()
  mplot.plot(X, y[:,2], lw=2.5, color='k')
  mplot.xlabel('x (meters)', fontsize=18)
  mplot.axvline(X[N/2]-x_p, color='k', linestyle='--')
  mplot.axvline(X[N/2]+x_n, color='k', linestyle='--')
  mplot.ylabel(r'$\phi$ (Volts)', fontsize=18)
  YRange = mplot.ylim()
  YMin   = YRange[0] ; YMax = YRange[1]
  mplot.text(XJunction-x_p-.05*XFullRange,  (YMin+YMax)/2, r'$x_{p}$', \
             fontsize=16,color='k')
  mplot.text(XJunction+x_n+.02*XFullRange,  (YMin+YMax)/2, r'$x_{n}$', \
             fontsize=16,color='k')
  #Fill in the P and N regions
  YRange = mplot.ylim()
  YMin   = YRange[0] ; YMax = YRange[1]
  mplot.fill([X[0],X[N/2],X[N/2],X[0]], [YMin,YMin,YMax,YMax],\
              facecolor="#ffdede")
  mplot.fill([X[N/2],X[NN],X[NN],X[N/2]], [YMin,YMin,YMax,YMax],\
              facecolor="#dedeff")
  Ax = gca()
  Ax.set_ylim([YMin, YMax])
  mplot.draw()

  if SaveFig == 1:
    mplot.savefig(PlotDir+'PNJunction_'+ModeStr+'_Rho_'+str(NSteps)+'Steps.png')
    mplot.savefig(PlotDir+'PNJunction_'+ModeStr+'_Rho_'+str(NSteps)+'Steps.eps')