%Fig4_3.m
%Essential Electron Transport for Device Physics
%Input:  Potential profile, doping, temperature
%Output: Current density as function of applied voltage bias

%This program uses 
% - CalculateChemicalPotential.m to calculate the chemical potential Mu 
%   for given temperature and doping concentration
% - PotentialProfileWithPoissonB.m to calculate the form of the potential
%   profile in the charge accumulation/depletion region by solving Poisson
%   equation in depletion approximation
% - option input parameters from file 'TunnelCurrentParameterInputB.txt'
% - option input potential profile from file
%   'TunnelCurrentPotentialProfileInputB.txt'

clear all;
clf;

FS      = 12;	%label fontsize 18
FSN     = 12;	%number fontsize 16
LW      = 1;	%linewidth

% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times');
set(0,'DefaultAxesFontSize', FSN);

% Change default text fonts.
set(0,'DefaultTextFontname', 'Times');
set(0,'DefaultTextFontSize', FSN);

%t = cputime;                    %start timer
eye = complex(0.,1.);			%square root of -1

hbar    = 1.05457159e-34;		%Planck's constant (Js)
echarge = 1.6021764e-19;		%electron charge (C)
kB      = 1.3806505e-23/echarge;%Boltzmann constant (eV/K)

m0      = 9.109382e-31;         %bare electron mass (kg)
eps0    = 8.854188e-12;         %permittivity of free space (F/m)

VbiasMin    = (pi*10^-9);       %min. applied voltage bias (eV) 
VbiasMax    = 0.40;             %max. applied voltage bias (eV)
VbiasPoints = 26;               %number of discrete voltage bias points, 26   
Vbias = linspace(VbiasMin,VbiasMax,VbiasPoints);%Voltage bias array 
                        
% input parameters from file 'ParameterInput.txt':
%[ParameterList]=textread('TunnelCurrentParameterInputB.txt','%*s %f %*s'); 
% m = ParameterList(1) * m0;         %effective electron mass (kg)
% eps1 = ParameterList(2) * eps0;    %permittivity of GaAs 13.1 (F/m)
% eps2 = ParameterList(3) * eps0;    %permittivity of AlAs 10.1 (F/m)
% monolayer = ParameterList(4);      %thickness of one monolayer GaAs (m) 
% Temperature = ParameterList(5);    %temperature (K)
% NDopingcm = ParameterList(6) 
% meff         0.07        (m/m0)  
% eps1         13.1        (eps(GaAs)/eps0) 
% eps2         13.1        (eps(AlAs)/eps0) 
% monolayer    0.282665e-9 (m)
% Temperature  300.0         (K)
% NDoping      1.0         (10^18/cm^(-3))
xi = 0.36;                 %Al fraction in GaAs
V0 = 0.84*xi;             %heterostructure conduction band offset (eV)
m  = 0.07 * m0;           %effective electron mass (kg) 0.07
eps1 = 13.1 * eps0;       %permittivity of GaAs 13.1 (F/m)
eps2 = 13.1 * eps0;       %permittivity of AlAs 10.1 (F/m)
monolayer = 0.282665e-9;  %thickness of one monolayer GaAs (m) 0.282665nm
Temperature = 300.0;      %temperature (K) 300
NDopingcm = 1.0;          %doping x10^18 cm^-3 

NDoping = NDopingcm*10^24;  %doping concentration in contacts (m^-3)
%end reading in parameters
if Temperature < 1; Temperature = 1; end    %never allow T=0 K

beta=1/(kB*Temperature);      %inverse thermal energy (eV^-1)

EFermi = ((hbar^2)*((3.0*(pi^2)*NDoping)^(2.0/3.0)))/(2.0*m*echarge); 
Emax = EFermi + (15.0/beta);  %cutoff for energy integral (eV)
nEnergy=2*VbiasPoints;        %number of energy steps in integration
Emin = (3*pi*10^-8);          %avoid divide by zero   
dEnergy = Emax/nEnergy;
const0  = (m*(echarge^3)*dEnergy*10^-4)/(beta*2*(pi^2)*(hbar^3));%10^-4 to convert from A/m^2 to A/cm^2
const2  = 2.0*echarge*m/(hbar^2);

%chemical potential for given temperature and doping 
Mu = CalculateChemicalPotential(Temperature,NDoping,m,EFermi);  
%set up energy array and precalculate bias independent terms
energyindex=[1:nEnergy];
Energy=Emin+dEnergy*double(energyindex);
expa=(beta*(Mu-Energy));
const1a = 1.0+(exp(expa));
nEnergyx2=2*nEnergy;
eindex=2*energyindex;

%setup potential profile:
%read data L(#monolayer), V(eV) from file 'PotentialProfileInput.txt'
%L(#monolayer) denotes the width of barrier region in the potential 
%profile with a certain value of potential V(eV)
%e.g.: for a single 4 monolayer thick 1 eV energy rectangular barrier, 
%we have 1 barrier region with L(1) = 4 monolayer and V(1) = 1 eV
%[L, V] = textread('TunnelCurrentPotentialProfileInputB.txt','%f %f','headerlines',2);    
%L = [0.01 0.01 0.01 0.01 5 5 5 4 4 5 5 5 0.01 0.01 0.01 0.01 0.01];
%V = [0 0 0 0 V0 V0 0 0 0 0 V0 V0 0 0 0 0 0];
L = [0.01 5 5 5 4 4 5 5 5 0.01];
V = [0 0 0 V0 V0 V0 V0 0 0 0];

L      = monolayer*L;%thickness of barrier in the potential profile (m)
NBarrier = length(L);%NBarrier = # of barriers in the potential profile

Lx(1) = 0.0;
for i = 2 : 1 : NBarrier+1
   Lx(i)      = Lx(i-1) + L(i-1);  
end

d2 = Lx(NBarrier+1);          %total thickness of potential profile (m) 

d1Max = -(eps1*d2)/(2.0*eps2)+sqrt(VbiasMax/echarge*eps1/NDoping+...
    (eps1*d2)/(2.0*eps2)*(eps1*d2)/(2.0*eps2)); 
%Width of charge accumulation/depletion region
%width of depletion region for maximum value of voltage bias,
%depending on NDoping

xmin = -(d1Max+2.0e-9);  %lower cutoff in x-direction, depending on NDoping                        
xmax = d2+d1Max+2.0e-9;  %upper cutoff in x-direction, depending on NDoping
xpoints = round((xmax-xmin)*1.0e10)*2; %# of points in x-direction for propagation matrix method
dx = (xmax-xmin)/xpoints;                 
E2=Energy.'*ones(1,xpoints);

for i = 1 : 1 : xpoints     %set up x-axis
   x(i)      = double(i-1) * dx + xmin;  
end

% fprintf(1,'\n');
% fprintf(1,'Voltage bias (V)  Current (A/cm^2)\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%outer loop calculates potential profile at each voltage bias Vbias(j)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Current=zeros(1,VbiasPoints);%initialize
t = cputime;                    %start timer
for j = 1:1:VbiasPoints
%calculate potential for each Vbias(j)
Vpot = PotentialProfileWithPoissonB(xpoints, x, dx, Lx, NBarrier, V, Vbias(j), NDoping, eps1, eps2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%loop to integrate over energy starting from Emin and stepping by dEnergy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 expb=(beta*(Mu-Energy-Vbias(j)));
 const1b =1.0+(exp(expb));
 const1ab=const1a./const1b;
 const1=log(const1ab);
  
V2=ones(nEnergy,1)*Vpot;
        
bigP2=speye(nEnergyx2); 
k2=sqrt(const2*(E2-V2));
for n=1:1:xpoints-1

    k2tmp1=exp(-eye*k2(:,n)*dx);
    k2tmp2=exp(eye*k2(:,n)*dx);
   p11=(1+k2(:,n+1)./k2(:,n)).*k2tmp1;
   p12=(1-k2(:,n+1)./k2(:,n)).*k2tmp1;
   p21=(1-k2(:,n+1)./k2(:,n)).*k2tmp2;
   p22=(1+k2(:,n+1)./k2(:,n)).*k2tmp2;
   p2=sparse(eindex-1,eindex-1,p11,nEnergyx2,nEnergyx2);
   p2=p2+sparse(eindex-1,eindex,p12,nEnergyx2,nEnergyx2);
   p2=p2+sparse(eindex,eindex-1,p21,nEnergyx2,nEnergyx2);
   p2=p2+sparse(eindex,eindex,p22,nEnergyx2,nEnergyx2);
   bigP2=0.5*bigP2*p2;
end
Trans=diag(bigP2);
Trans=Trans(eindex-1);
Trans=abs(1./Trans).^2;
Trans=Trans.*const1';

Current(j)=const0*sum(Trans);
  
%    fprintf(1,'    %1.3f              %5.9f \n', Vbias(j), Current(j));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%end of outer Vbias loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
run_time = cputime-t       %CPU run time
%plot results
CurrentMax=max(Current);

ttl=(['\rmFig4.3, \itn\rm=',num2str(NDopingcm,'%5.1f'),...
    '\times10^{18} cm^{-3}, \itT\rm=',...
    num2str(Temperature,'%5.0f'),' K, \mu=',num2str(Mu*1000,'%5.1f'),...
    ' meV, \itm\rm^*_e=',num2str(m/m0,'%7.2f'),'\timesm_0, \itV\rm_0=',...
    num2str(V0,'%7.2f'),' eV']); 

figure(1);
plot(Vbias,Current,'b','LineWidth',LW);
xlabel('Applied voltage, \itV\rm_{bias} (V)','FontSize',FS);
ylabel('Current density, \itJ\rm_{tot} (A/cm^{2})','FontSize',FS)
axis([0, VbiasMax, 0, 1.1*CurrentMax]);
grid on
title(ttl,'FontSize',FS);

figure(2);
plot(Vbias,log10(Current),'b','LineWidth',LW);
xlabel('Applied voltage, \itV\rm_{bias} (V)','FontSize',FS);
ylabel('log_{10} (\itJ\rm_{tot}) (A cm^{-2})','FontSize',FS);
axis([0, VbiasMax, log10(Current(2)), log10(1.1*CurrentMax)]);
grid on
title(ttl,'FontSize',FS);

figure(3);
plot(x*1e9,Vpot,'b','LineWidth',LW);
xlabel('Position, \itx\rm (nm)','FontSize',FS);
ylabel('Potential, \itV\rm (eV)','FontSize',FS);
axis([min(x*1e9),max(x*1e9), 1.1*min(Vpot),1.1*max(Vpot)]);
title(ttl,'FontSize',FS);

