%% Finite Difference based 1D Nonlinerar Possion Equation Solver for pn-Junctions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name                  : pn_poisson_calc.mlx
% Author                : Patrick Gsoels (11771814)
% Contact               : patrick.gsoels@student.tugraz.at
% Creation Date         : 16.07.2022
% Version               : 1.1
% Description           : This live script computes the charge carrier
%                         concentrations, charge-density, electric field,
%                         elctric potential and band-structure of arbitrary
%                         pn-junction doping profiles.
% Language              : MATLAB R2020b
% Category              : live script
% To Do                 : -
% -------------------------------------------------------------------------
% 1:01: initial version [PG]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

% Table of Contents
% Poisson Equation for Electrostatics
% Newton-Raphson Method﻿
% Finite Difference Formulation
% Discretization of the 1D Poisson Equation
% Declaration of Constants
% Grid Specifications
% Doping Profile
% Newton-Raphson Parameters
% Implementation of a Newton-Raphson Assisted Finite Difference 1D Poisson Solver
% Plotting the Results
%     Doping Profile
%     Charge Carrier Concentrations
%     Charge Density
%     Internal Electric Field
%     Build-in Potential
%     Band Diagram
% References
%% Poisson Equation for Electrostatics
% In electrodynamics Gauss's law states the relation between the electric displacement field \textbf{D} and the free charge density \rho within a volume. The differential form of Gauss's law can be written as
% \nabla \cdot \textbf{D} = \rho.    (1)
% For a homogenious medium the following relation between electric displacement field \textbf{D} and electric field \textbf{E} holds:
% \textbf{D} = \varepsilon \cdot \textbf{E}    (2)
% In electrostatics we do not have any time-dependent fields. Hence, Faraday's law of induction will simplify to the following expression:
% \nabla \times \textbf{E} = \textbf{0}     (3)
% The simplified version of Faraday's law of induction states that the electric field \textbf{E} is free of curl. Thus, the electric field \textbf{E} is a conservative field for which one can define a potential.
% \textbf{E} = -\nabla V    (4)
% Inserting Equation (2) and (4) into Equation (1) will give us the Poisson equation for electrostatics.
% \nabla \cdot \left( \varepsilon \cdot \nabla V \right) = -\rho    (5)
% For a constant dielectric media Equation (5) simplifies to
%  \nabla^2 V\left(x,y,z \right) = -\frac{\rho\left(x,y,z \right)}{\varepsilon}.    (6)
% pn-junctions do typically have doping profiles which intentionally vary along the latteral axis which is perpendicular to the junction interface. Thus, the Poisson equation from (6) reduces to one spacial coordinate.
%  \nabla^2 V\left(x \right) = \frac{\text{d}^2V}{\text{d}x^2} = -\frac{\rho\left(x \right)}{\varepsilon}.    (7)
% Due to the interconnection of differently doped materials, charge carriers diffuse accros the pn-junction. At the same time an internal electric field establishes due to the movement of charge carriers. These facts lead to the necessity of introducing a potential-dependent charge density \rho, which reformulates the Poisson equation to
%  \nabla^2  V\left(x \right) = \frac{\text{d}^2V}{\text{d}x^2} = -\frac{\rho\left( V\left(x \right) \right)}{\varepsilon}.    (8)
% ﻿According to [1], assuming the dopants to be totally ionised, the charge density inside a semiconductor is given by
% \rho\left( x \right) = q\cdot \left[ p\left( V\left(x \right) \right) - n\left( V\left(x \right) \right) + N_D(x) - N_A(x) \right].    (9)
% According to [2], the electron n(x) and hole density p(x) can be formulated as
% n(x) = n_i \cdot \text{exp} \left( \frac{E_{Fn} - E_i}{k_B \cdot T} \right) = n_i \cdot \text{exp} \left( \frac{q\cdot V(x) }{k_B \cdot T} \right)    (10)
% and
% p(x) = n_i \cdot \text{exp} \left( \frac{E_i - E_{Fp}}{k_B \cdot T} \right) = n_i \cdot \text{exp} \left( \frac{-q\cdot V(x) }{k_B \cdot T} \right).    (11)
% Inserting Equation (10) and (11) into (9) will provide us with the full formulation of the potential-dependent charge density \rho.
% \rho\left( x \right) = q\cdot \left[ n_i \cdot \text{exp} \left( \frac{-q\cdot V(x) }{k_B \cdot T} \right) -  n_i \cdot \text{exp} \left( \frac{q\cdot V(x) }{k_B \cdot T} \right) + N_D(x) - N_A(x) \right]    (12)
% with n_i being the intrinsic carrier concentration of the semiconductor, which is defined as
% n_i^2 = n\cdot p = N_c(T_{300K}) \cdot \left(\frac{T}{T_{300K}} \right)^{\frac{3}{2}} \cdot N_v(T_{300K}) \cdot \left(\frac{T}{T_{300K}} \right)^{\frac{3}{2}}  \cdot \text{exp}\left( \frac{-E_g}{k_B \cdot T} \right).    (13)
% Bear in mind that the formulars for n,p and n_i are derived by using the Boltzmann approximation instead of the real fermi-dirac distribution. Hence, the condition \left| E - E_F \right| > 3 \cdot k_B \cdot T has to hold.
% With the help of some trigonometric identity (e^x - e^{-x} = 2\cdot \sinh(x)), Equation (12) can be reformulated as
% \rho\left( x  \right) = q\cdot \left[ N_D(x) - N_A(x)  - 2\cdot n_i \cdot \sinh \left( \frac{q\cdot V(x) }{k_B \cdot T} \right) \right]
%
% .    (14)
% Since our one-dimensional Poisson equation is non-linear, one will need at some point the partial derivative of the charge density \rho with respect to the electric potential V.
% \frac{\partial \rho}{\partial V} =  - \frac{2\cdot q^2 \cdot n_i}{k_B \cdot T} \cdot \cosh \left( \frac{q\cdot V(x) }{k_B \cdot T} \right)
%
%     (15)
%% Newton-Raphson Method﻿
% The Newton-Raphson method is mathematically the most preferred of the several known methods for the solution of systems of nonlinear equations. It is a successive approximation technique based on an initial estimate of the solution and the Taylor series expansion [1]. ﻿Consider the set of equations
% \textbf{f} \left( \textbf{V} \right)  = \left[ \begin{array}{c}
% f_1\left( \textbf{V} \right)\\
% \vdots\\
% f_N\left( \textbf{V} \right)\\
% \end{array} \right] =\textbf{0}    (16)
% A Taylor approximation of (16), which is terminated after the linear term, leads to
% \textbf{f}\left( \textbf{V}^{(k+1)} \right) \approx \textbf{f}\left( \textbf{V}^{(k)} \right) + \mathcal{J}\left( \textbf{V}^{(k)}  \right)  \cdot \left(\textbf{V}^{(k+1)} - \textbf{V}^{(k)}  \right) = \textbf{f}\left( \textbf{V}^{(k)} \right) + \mathcal{J}\left( \textbf{V}^{(k)}  \right)  \cdot \Delta\textbf{V}^{(k)} \overset{!}{=} \textbf{0}.    (17)
% \textbf{f}\left( \textbf{V}^{(k+1)} \right) \approx  \textbf{f}\left( \textbf{V}^{(k)} \right) + \mathcal{J}\left( \textbf{V}^{(k)}  \right)  \cdot \Delta\textbf{V}^{(k)} =
% \left[ \begin{array}{c}
% f_1\left( \textbf{V}^{(k)} \right) + \sum\limits_{j=1}^{N} \left.\frac{\partial f_1}{\partial V_j^{(k)}}\right|_{\textbf{V}^{(k)}} \cdot \Delta V_j^{(k)}  \\
% \vdots\\
% f_N\left( \textbf{V}^{(k)} \right) + \sum\limits_{j=1}^{N} \left.\frac{\partial f_N}{\partial V_j^{(k)}}\right|_{\textbf{V}^{(k)}} \cdot \Delta V_j^{(k)} \\
% \end{array} \right]
%
%
% \overset{!}{=} \textbf{0}.    (18)
% with
% \mathcal{J}\left( \textbf{V}^{(k)}  \right)  = \left[ \left.\frac{\partial \textbf{f}\left(\textbf{V}\right) }{\partial V_1}\right|_{\textbf{V} = \textbf{V}^{(k)}} \ \ldots\ \left.\frac{\partial \textbf{f}\left(\textbf{V}\right) }{\partial V_N}\right|_{\textbf{V} = \textbf{V}^{(k)}} \right] =
% \left[ \begin{array}{c}
% \nabla^T \left. f_1\left( \textbf{V} \right) \right|_{\textbf{V}=\textbf{V}^{(k)}}   \\
% \vdots\\
% \nabla^T \left. f_M\left( \textbf{V} \right) \right|_{\textbf{V}=\textbf{V}^{(k)}} \\
% \end{array} \right] =
%
% \left[ \begin{array}{ccc}
% \left.\frac{\partial f_1\left(\textbf{V}\right) }{\partial V_1}\right|_{\textbf{V} = \textbf{V}^{(k)}} & \cdots & \left.\frac{\partial f_1\left(\textbf{V}\right) }{\partial V_N}\right|_{\textbf{V} = \textbf{V}^{(k)}}  \\
% \vdots &\ddots & \vdots \\
% \left.\frac{\partial f_M\left(\textbf{V}\right) }{\partial V_1}\right|_{\textbf{V} = \textbf{V}^{(k)}} & \cdots & \left.\frac{\partial f_M\left(\textbf{V}\right) }{\partial V_N}\right|_{\textbf{V} = \textbf{V}^{(k)}} \\
% \end{array} \right]    (19)
% being the Jaccobian.
% After rearranging (17), one gets the expression for the Newton-Raphson step
% \Delta\textbf{V}^{(k)} = - \mathcal{J}\left( \textbf{V}^{(k)}  \right) \cdot \textbf{f}\left( \textbf{V}^{(k)} \right)    (20)
% and the (k+1)-th iterate
% \textbf{V}^{(k+1)} = \textbf{V}^{(k)}  + \Delta\textbf{V}^{(k)} .    (21)
%% Finite Difference Formulation
% The finite difference formulation can be again derived from the Taylor approximation of V(x) as depicted below.
% V\left( x_j + \Delta x_j \right) = V\left( x_j \right) + \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \cdot \Delta x_j +
% \frac{1}{2!} \cdot \left. \frac{\text{d}^2 V\left( x \right)}{\text{d}x^2} \right|_{x=x_j} \cdot \Delta x_j^2 +
% \frac{1}{3!} \cdot \left. \frac{\text{d}^3 V\left( x \right)}{\text{d}x^3} \right|_{x=x_j} \cdot \Delta x_j^3 + \ldots = V\left( x_j \right) + \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \cdot \Delta x_j + \mathcal{O} \left(\Delta x_j^2 \right) \approx V\left( x_j \right) + \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \cdot \Delta x_j      (22)
% \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \approx {V_j^{(FD)}}' = \frac{V\left( x_j + \Delta x_j \right) - V\left( x_j \right)}{\Delta x_j} = \frac{V_{j+1} - V_{j}}{\Delta x_j} [forward difference]    (23)
% V\left( x_j - \Delta x_j \right) = V\left( x_j \right) - \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \cdot \Delta x_j +
% \frac{1}{2!} \cdot \left. \frac{\text{d}^2 V\left( x \right)}{\text{d}x^2} \right|_{x=x_j} \cdot \Delta x_j^2 -
% \frac{1}{3!} \cdot \left. \frac{\text{d}^3 V\left( x \right)}{\text{d}x^3} \right|_{x=x_j} \cdot \Delta x_j^3 + \ldots = V\left( x_j \right) - \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \cdot \Delta x_j + \mathcal{O} \left(\Delta x_j^2 \right) \approx V\left( x_j \right) + \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \cdot \Delta x_j       (24)
% \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \approx {V_j^{(BD)}}' = \frac{V\left( x_j \right) - V\left( x_j - \Delta x_j\right)}{\Delta x_j} = \frac{V_{j} - V_{j-1}}{\Delta x_j} [backward difference]    (25)
% One gets the central difference formulation for the first order derivative if one subtracts Equation (24) from (22) like shown below.
% V\left( x_j + \Delta x_j \right)  - V\left( x_j - \Delta x_j \right) = 2 \cdot \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \cdot \Delta x_j + 2\cdot \mathcal{O} \left(\Delta x_j^3 \right) \approx 2 \cdot \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \cdot \Delta x_j    (26)
% \left. \frac{\text{d}V\left( x \right)}{\text{d}x} \right|_{x=x_j} \approx {V_j^{(CD)}}' = \frac{V\left( x_j + \Delta x_j \right) - V\left( x_j - \Delta x_j \right)}{2 \Delta x_j} = \frac{V_{j+1} - V_{j-1}}{2\Delta x_j} [central difference]    (27)
% The second order derivative can be approximated by inserting Equation (23) and (25) into the expression depicted below.
% \left. \frac{\text{d}^2 V\left( x \right)}{\text{d}x^2} \right|_{x=x_k} \approx {V_j}'' = \frac{{V_j^{(FD)}}' - {V_j^{(BD)}}' }{\Delta x_j} = \frac{V_{j+1} - 2\cdot V_j + V_{j-1}}{\Delta x_j^2}    (28)
%% Discretization of the 1D Poisson Equation
% With the aid of the finite difference formulation one can approximate Equation (8), which results in
% \frac{V_{j+1} - 2\cdot V_j + V_{j-1}}{\Delta x_j^2} = -\frac{\rho\left(V_j \right)}{\varepsilon}, j=1,\ldots,N.    (29)
% After rearranging Equation (29) one gets the non-linear algebraic difference formulation of the Poisson Equation.
% f_i\left( \textbf{V} \right) = \frac{V_{j+1} - 2\cdot V_j + V_{j-1}}{\Delta x_j^2}
% + \frac{\rho\left(V_j \right)}{\varepsilon}  = 0, j=1,\ldots,N.    (30)
% By using the Newton-Raphson technique for solving non-linear systems of equations, the non-linear algebraic difference equation (30) formulates to
% f_i\left( \textbf{V}^{(k+1)} \right) \approx \left[ f_i\left( \textbf{V}^{(k)} \right) \right] + \left\lbrace  \sum\limits_{j=1}^{N} \left.\frac{\partial f_i}{\partial V_j}\right|_{\textbf{V} = \textbf{V}^{(k)}} \cdot \Delta V_j^{(k)} \right\rbrace  = \left[ \frac{V_{j+1}^{(k)} - 2\cdot V_j^{(k)} + V_{j-1}^{(k)}}{\Delta x_j^2}
% + \frac{\rho\left(V_j^{(k)} \right)}{\varepsilon} \right] + \left\lbrace \frac{\Delta V_{j+1}^{(k)} - 2\cdot \Delta V_j^{(k)} + \Delta V_{j-1}^{(k)}}{\Delta x_j^2} + \left. \frac{\partial \rho}{\partial V}\right|_{V=V_j^{(k)}} \cdot \Delta V_j^{(k)} \right\rbrace = 0.    (31)
% The known parts of Equation (31) will define our right hand-side and the unknown parts will define the left-hand side of the system of equations which needs to be solved for \Delta V_j^{(k)} .
% - \left\lbrace \frac{\Delta V_{j+1}^{(k)} - 2\cdot \Delta V_j^{(k)} + \Delta V_{j-1}^{(k)}}{\Delta x_j^2} + \left. \frac{\partial \rho}{\partial V}\right|_{V=V_j^{(k)}} \cdot \Delta V_j^{(k)} \right\rbrace = \left[ \frac{V_{j+1}^{(k)} - 2\cdot V_j^{(k)} + V_{j-1}^{(k)}}{\Delta x_j^2}
% + \frac{\rho\left(V_j^{(k)} \right)}{\varepsilon} \right]    (32)
% ﻿By defining the right-hand side of Equation (32) as
% R_j^{(k)} = \frac{V_{j+1}^{(k)} - 2\cdot V_j^{(k)} + V_{j-1}^{(k)}}{\Delta x_j^2}
% + \frac{\rho\left(V_j^{(k)} \right)}{\varepsilon}    (33)
% and the left-hand side as
% M_{ji}^{(k)} \cdot \Delta V_j^{(k)} = \left\{ \begin{array}{cl}
% \left( \frac{2}{\Delta x_j^2} - \left. \frac{\partial \rho}{\partial V}\right|_{V=V_j^{(k)}}\right) \cdot \Delta V_j^{(k)} & \textrm{, } i=j\\
% -\frac{1}{\Delta x_j^2} \cdot \Delta V_j^{(k)} & \textrm{, } i = \left| j \pm 1 \right| \\
% 0 & \textrm{, else} \\
% \end{array}\right..    (34)
% With Equation (33) and (34) one can easily write Equation (32) in matrix form
% \left[ \begin{array}{ccccccc}
% M_{1,1}^{(k)} & M_{1,2}^{(k)} & 0 & & & &  \\
% M_{2,1}^{(k)} & M_{2,2}^{(k)} & M_{2,3}^{(k)} & 0 & & &  \\
% 0 & M_{3,2}^{(k)} & M_{3,3}^{(k)} & M_{3,4}^{(k)} & 0 & &  \\
%  & \ddots & \ddots & \ddots & \ddots & \ddots &  \\
%  & & 0 & M_{(N-2),(N-3)}^{(k)} & M_{(N-2),(N-2)}^{(k)} & M_{(N-2),(N-1)}^{(k)} & 0  \\
%  & & & 0 & M_{(N-1),(N-2)}^{(k)} & M_{(N-1),(N-1)}^{(k)} & M_{(N-1),N}^{(k)}  \\
%  & & & & 0 & M_{N,(N-1)}^{(k)} & M_{N,N}^{(k)}  \\
% \end{array} \right] \cdot
% \left[ \begin{array}{c}
% \Delta V_1^{(k)}  \\
% \Delta V_2^{(k)}  \\
% \Delta V_3^{(k)}  \\
% \vdots \\
% \Delta V_{N-2}^{(k)}  \\
% \Delta V_{N-1}^{(k)}  \\
% \Delta V_N^{(k)}  \\
% \end{array} \right] =
% \left[ \begin{array}{c}
% R_1^{(k)}  \\
% R_2^{(k)}  \\
% R_3^{(k)}  \\
% \vdots \\
% R_{N-2}^{(k)}  \\
% R_{N-1}^{(k)}  \\
% R_N^{(k)}  \\
% \end{array} \right]    (35)
%  with \textbf{V}^{(0)} = \textbf{V}_{init}, V_0^{(k)} = V_{init}\left( x_0 \right), V_{N+1}^{(k)}  = V_{init}\left( x_{N+1} \right), \Delta V_0^{(k)} = \Delta V_{N+1}^{(k)} = 0 and
% V_{init}\left( x_j \right) = \left\{ \begin{array}{cl}
% \textrm{sign}\left( N_D^{(j)} - N_A^{(j)} \right) \cdot \frac{k_B \cdot T}{q} \cdot \ln \left( \frac{\left| N_D^{(j)} - N_A^{(j)} \right|}{n_i} \right) & \textrm{, } \left| N_D^{(j)} - N_A^{(j)} \right| \geq n_i \\
% 0 & \textrm{, else} \\
% \end{array}\right..    (36)
% Finally, one has to solve Equation (35) for the Newton-Raphson step \Delta \textbf{V}^{(k)}  in order to generate the next iterate of the potential \textbf{V}^{(k+1)} as defined by Equation (21).
% \textbf{V}^{(k+1)} = \textbf{V}^{(k)}  + \Delta\textbf{V}^{(k)}     (37)
%% Declaration of Constants

% temperature
% (temperature-dependent band-gap and intrinsic carrier concentration is not implemented)
T = 300;                 % in degree K

% permittivity of vacuum
eps_0 = 8.854187813e-12; % in As/Vm

% relative permittivity
eps_r = 12;

% absolute permittivity
epsilon = eps_0 * eps_r; % in As/Vm

% Boltzmann constant
kB = 1.380658e-23;       % in J/K

% elementary charge
q = 1.60217733e-19;      % in As

% intrinsic carrier concentration of SI
ni = 1.45e16;            % in 1/m^3

% band gap energy of SI
Eg = 1.12;               % in eV

%% Grid Specifications

% lower grid bound
xmin = str2double("-4e-6"); % in m

% upper grid bound
xmax = str2double("4e-6"); % in m

% number of grid points
N = round(str2double("1000"));

%% Doping Profile

doping_profile = "abrupt_pnp";

%% Newton-Raphson Parameters

maxIter = round(str2double("1000"));
alpha = 1e-6;

%% Implementation of a Newton-Raphson Assisted Finite Difference 1D Poisson Solver

% grid points (sampling)
x = linspace(xmin,xmax,N).';

% step-size
dx = x(2)-x(1);

% grid extension (needed for central differences)
xg = [ min(x)-dx;x;max(x)+dx ];

% load doping profile
[NA,ND] = doping_profiles(doping_profile,x);

% initializing electron concentration
n = ND;

% initializing hole concentration
p = NA;

% initializing electric potential
V = zeros(size(xg));
Vinit = band_initialization(NA,ND,ni,q,kB,T);
V(2:N+1) = Vinit;
V(1) = V(2);       % boundary
V(end) = V(end-1); % boundary

V_old = V;

% figure;
% plot(x,V(2:N+1))

% loop for Newton-Raphson Iterations
for k=1:maxIter 
  
  % updating charge carrier concentrations
  n = ni*exp(q*V(2:N+1)/(kB*T));
  p = ni*exp(-q*V(2:N+1)/(kB*T));
  
  % charge density
  % (assuming that the dopants are totally ionized)
  rho = q*(p - n + ND - NA);
  
  % charge density change due to the internal potential
  drho_dV = -2*q^2*ni/(kB*T) * cosh(q*V(2:N+1)/(kB*T));
  
  % finite difference formulation for the potential
  d2V_dx2 = ( V(1:N,1) - 2*V(2:N+1,1) + V(3:N+2,1) )/(dx^2);
  
  % non-linear equation system
  Mjj_k = 2/(dx^2) - 1/epsilon*drho_dV;
  M_off = 1/(dx^2)*ones(N,1);
  M_k = spdiags([-M_off Mjj_k -M_off],-1:1,N,N);
  R_k = 1/epsilon*rho + d2V_dx2;
  
  % solving for potential changes
  delta_V = M_k\R_k;

  % Newton-Raphson Step (updating the potential)
  %V(2:N+1) = V(2:N+1)+(0.98^k)*delta_V;
  V(2:N+1) = V(2:N+1)+delta_V;
  
  % stopping criteria
  if( norm( (V - V_old), "inf" ) < alpha )
    break;
  end  
  
  V_old = V;
end

% calculating the electric field
E = -diff(V(1:N+1))/dx;

% build in voltage
Vbi = abs(V(end)-V(1));

disp(['Number of Newton-Raphson Iterations: ',num2str(k)])
disp(['build in voltage: Vbi = ',num2str(Vbi),'V'])

%% Plotting the Results

set(groot, ...
'DefaultTextInterpreter', 'LaTeX',...
'DefaultAxesTickLabelInterpreter', 'LaTeX',...
'DefaultLegendInterpreter', 'LaTeX',...
'DefaultAxesFontName', 'LaTeX',...
'DefaultFigureColor', 'w', ...
'DefaultAxesLineWidth', 0.5, ...
'DefaultAxesXColor', 'k', ...
'DefaultAxesYColor', 'k', ...
'DefaultAxesFontUnits', 'points', ...
'DefaultAxesFontSize', 15, ...
'DefaultAxesFontName', 'Helvetica', ...
'DefaultLineLineWidth', 1.5, ...
'DefaultTextFontUnits', 'Points', ...
'DefaultTextFontSize', 20, ...
'DefaultTextFontName', 'Helvetica', ...
'DefaultAxesBox', 'on', ...
'DefaultAxesTickLength', [0.01 0.01],...
'defaultAxesXGrid','on',...
'defaultAxesYGrid','on',...
'defaultfigureposition',[50 50 1080 400],...
'defaultAxesTitleFontSize',1.3);

%% Doping Profile

% plot chosen doping profile
figure, hold on, grid minor;
plot(x,NA,'b','DisplayName','$N_A$','LineWidth',1.5)
plot(x,ND,'r','DisplayName','$N_D$','LineWidth',1.5)
plot(x,(ND-NA),'g--','DisplayName','$N_D - N_A$','LineWidth',1.5)
title('doping concentrations')
xlabel('$x$ in $m$')
ylabel('$N(x)$ in $1$/$m^3$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

%% Charge Carrier Concentrations

figure, hold on, grid minor;
plot(x,p,'b','DisplayName','$p$','LineWidth',1.5)
plot(x,n,'r','DisplayName','$n$','LineWidth',1.5)
title('carrier concentrations')
xlabel('$x$ in $m$')
ylabel('$n(x),p(x)$ in $1$/$m^3$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

%% Charge Density

figure, hold on, grid minor;
plot(x,rho,'b','DisplayName','$\rho$','LineWidth',1.5)
title('charge density')
xlabel('$x$ in $m$')
ylabel('$\rho(x)$ in $As$/$m^3$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

%% Internal Electric Field

figure, hold on, grid minor;
plot(x,E,'b','DisplayName','$E$','LineWidth',1.5)
title('electric field')
xlabel('$x$ in $m$')
ylabel('$E(x)$ in $V$/$m$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

%% Build-in Potential

figure, hold on, grid minor;
plot(x,V(2:N+1),'b','DisplayName','$V$','LineWidth',1.5)
title('potential')
xlabel('$x$ in $m$')
ylabel('$V(x)$ in $V$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

%% Band Diagram

figure, hold on, grid minor;
plot(x,-V(2:N+1)+Eg/2,'r','DisplayName','conduction band','LineWidth',1.5)
plot(x,-V(2:N+1)-Eg/2,'b','DisplayName','valence band','LineWidth',1.5)
title('band diagram')
xlabel('$x$ in $m$')
ylabel('energy in $eV$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

%

function V = band_initialization(NA,ND,ni,q,kB,T)

  % if needed, add bias voltage to this function

  dN = ND - NA;
  V = zeros(size(dN));

  for i=1:length(V)
    if( abs(dN) >= ni )
      V(i) = sign(dN(i))*kB*T/q*log(abs(dN(i))/ni);
    end
  end
end

function [NA,ND] = doping_profiles(str,x)
  switch str
    case 'abrupt' 
      NA = 1E14*ones(size(x));                                  % 1/cm^3
      NA = NA*1e6;                                              % 1/m^3
      ND = [ zeros(size(x(x<0))); ones(size(x(x>=0))) ] * 2E14; % 1/cm^3
      ND = ND*1e6;                                              % 1/m^3
    case 'linear'
      alpha  = 1e19*100e6;  % 1/m^4
      NA = zeros(size(x));
      ND = zeros(size(x));

      for i=1:length(x)
        if( x(i) <= 0 )
          NA(i) = -alpha*x(i);
        elseif( x(i) > 0  )  
          ND(i) = alpha*x(i);
        end
      end  
    case 'hyper_abrupt'
      NA=1E15*exp(0.5*(x/1e-6)).*heaviside(-x); % 1/cm^3
      NA = NA*1e6;                        % 1/m^3
      ND=1E15*exp(0.5*(-x/1e-6)).*heaviside(x); % 1/cm^3
      ND = ND*1e6;                        % 1/m^3
    case 'abrupt_pnp'
      NA = ( heaviside(-x-1e-6)*5E15 + heaviside(x+1e-6)*1E15 )*1e6; % 1/m^3
      ND = ( heaviside(x+1e-6)*3E15 - heaviside(x-1e-6)*3E15 )*1e6;  % 1/m^3
    case 'abrupt_npn'
      ND = ( heaviside(-x-1e-6)*5E15 + heaviside(x+1e-6)*1E15 )*1e6; % 1/m^3
      NA = ( heaviside(x+1e-6)*3E15 - heaviside(x-1e-6)*3E15 )*1e6;  % 1/m^3
    case 'abrupt_pnpn'
      NA = ( heaviside(-x-2e-6)*1E18 + ( heaviside(x-2e-6) - heaviside(x-3e-6) )*0.8E18 )*1e6; % 1/m^3
      ND = ( ( heaviside(x+2e-6) - heaviside(x-2e-6) )*5E17 + heaviside(x-3e-6)*1.2E18 )*1e6;  % 1/m^3
    case 'abrupt_npnp'
      ND = ( heaviside(-x-2e-6)*1E18 + ( heaviside(x-2e-6) - heaviside(x-3e-6) )*0.8E18 )*1e6; % 1/m^3
      NA = ( ( heaviside(x+2e-6) - heaviside(x-2e-6) )*5E17 + heaviside(x-3e-6)*1.2E18 )*1e6;  % 1/m^3
    case 'arbitrary'  
      [NA,ND] = myProfile(x);
    otherwise
      error(['[ERROR]: ',str,' is not a valid doping profile!'])
  end
end

function [NA,ND] = myProfile(x)
  NA = 1E14*ones(size(x));       % 1/cm^3
  NA = NA*1e6;                   % 1/m^3

  ND = exp( -(x/1e-6).^2 )*1e15; % 1/cm^3
  ND = ND*1e6;                   % 1/m^3
end

%% References
% [1] Jabr, R.A. & Hamad, Mohammed & Mohanna, Yasser. (2007). Newton-Raphson Solution of Poisson's Equation in a Pn Diode. International Journal of Electrical Engineering Education. 44. 23-33. 10.7227/IJEEE.44.1.3. https://www.researchgate.net/publication/233646108_Newton-Raphson_Solution_of_Poisson%27s_Equation_in_a_Pn_Diode
% [2] SZE, S. M. (1985). Semiconductor devices, physics and technology. New York, Wiley.
% [3] F. R. Shapiro, "The numerical solution of Poisson's equation in a pn diode using a spreadsheet," in IEEE Transactions on Education, vol. 38, no. 4, pp. 380-384, Nov. 1995, doi: 10.1109/13.473161.
% [4] Chapra, Steven C, and Raymond P. Canale. Numerical Methods for Engineers. Boston: McGraw-Hill Higher Education, 2006. Print.
% [5] Lecture 9 (Computational Electromagnetics) - Finite Difference Method, Dr. Raymond Rumpf, https://www.youtube.com/watch?v=v-exTNOSG3g