%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name                  : pn_poisson_calc.m
% Author                : Patrick Gsoels (11771814)
% Contact               : patrick.gsoels@student.tugraz.at
% Creation Date         : 16.07.2022
% Version               : 1.1
% Description           : This 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              : script
% To Do                 : -
% -------------------------------------------------------------------------
% 1:01: initial version [PG]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;
clc;

%% 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 = -4e-6; % in m

% upper grid bound
xmax =  4e-6; % in m

% number of grid points
N = 1000;

%% Doping Profile

profile_list = {'abrupt',...
                'linear',...
                'hyper_abrupt',...
                'abrupt_pnp',...
                'abrupt_npn',...
                'abrupt_pnpn',...
                'abrupt_npnp',...
                'arbitrary'};

doping_profile = profile_list{4};

%% Newton-Raphson Parameters

maxIter = 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;

%% Functions

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