##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^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 *
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  = 15
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
TP     = 3000                           # Total number of points on grid
Num_P  = double(TP/3)                   # Total number of points in P+ region
Num_I  = double(TP/3)                   # Total number of points in I  region
Num_N  = double(TP/3)                   # Total number of points in N+ region
Num_Pf = double(TP/3)
Num_PIf = double(TP/24)
Num_Ic  = double(TP/4)
Num_INf = double(TP/24)
Num_Nf  = double(TP/3)

#Get the indices of the boundaries
PI_B   = Num_Pf                         # P/I Interface index
IN_B   = Num_Pf+Num_PIf+Num_Ic+Num_INf  # I/N Interface index
I_fc_B = Num_Pf+Num_PIf
I_cf_B = Num_Pf+Num_PIf+Num_Ic

##Spacing and actual physical locations of coarse/fine transitions
X_l    = 100e-6                                      #Total length of diode
dx_pf  = 1e-6/Num_P                                  #spacing of fine region
dx_nf  = 1e-6/Num_N                                  #spacing of fine region
X_c    = X_l - (2e-6+dx_pf*Num_PIf+dx_nf*Num_INf)    #Total length of coarse region
dx_c   = X_c/Num_Ic                                  #spacing of coarse region

##Locations of transition regions
x_o      = 0
x_p_b    = x_o + Num_Pf*dx_pf
x_i_fc_b = x_p_b+Num_PIf*dx_pf
x_i_cf_b = x_i_fc_b+X_c
x_n_b    = x_i_cf_b+Num_INf*dx_nf
x_f      = X_l

##Full array of x coordinates is X
X = zeros(TP, dtype=double)                                 #In meters
x_pf   = frange(x_o,      x_p_b-dx_pf,    dx_pf)         #P Region
X[0:PI_B] = x_pf                                             #P Region
x_pi_f = frange(x_p_b,    x_i_fc_b-dx_pf, dx_pf)    #I Region
X[PI_B:I_fc_B] = x_pi_f                                  #I Region
x_i_c  = frange(x_i_fc_b, x_i_cf_b-dx_c,  dx_c)
X[I_fc_B:I_cf_B] = x_i_c
x_in_f = frange(x_i_cf_b, x_n_b-dx_nf,    dx_nf)
X[I_cf_B:IN_B] = x_in_f 
x_nf   = frange(x_n_b,    x_f-dx_nf,      dx_nf)                #N Region
X[IN_B: TP] = x_nf                                   #N Region

##The changes in dx between i and i +/- 1/2
dx   = X[1:TP]-X[0:TP-1]
dx_p = (dx[1:TP-1]+dx[0:TP-2])/2
N   = len(X)
NN  = len(X)-1
MN  = N/2

##Intrinsic Concentration
n_i = 1.e16
n_bulk = 1.e17

##Donor Concentration
BOff = 200
N_d_max = 1.e24
N_d = zeros(TP, dtype=double)
N_d[0:PI_B-BOff] = 0
N_d[PI_B-BOff:IN_B+BOff] = double(n_bulk)
N_d[IN_B+BOff:TP] = double(N_d_max)

##Acceptor Concentration
N_a_max = 1.e24
N_a = zeros(TP, dtype=double)
N_a[0:PI_B-BOff] = double(N_a_max)
N_a[PI_B-BOff:IN_B+BOff] = 0
N_a[IN_B+BOff:TP] = 0

##Difference is doping density
DC = N_d-N_a
#########################################################################

##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   = zeros(N, dtype=double)
p   = zeros(N, dtype=double)

##Set up concentrations
n[0:PI_B-BOff] = n_i**2/N_a_max
n[PI_B-BOff:IN_B+BOff] = double(n_bulk) 
n[IN_B+BOff:TP]        = N_d_max

p[0:PI_B-BOff] = N_a_max
p[PI_B-BOff:IN_B+BOff] = n_i**2/n_bulk
p[IN_B+BOff:TP] = n_i**2/N_d_max

##Boundary conditions according to Kurata p 37
n[NN]    = DC[NN]/2*(1+(1+(2*n_i/DC[NN])**2)**.5)
p[0]     = -DC[0]/2*(1+(1+(2*n_i/DC[0])**2)**.5)
n[0]     = n_i**2/p[0]
p[NN]    = n_i**2/n[NN]

##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
  pdb.set_trace()
  ##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])

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')