function TunnelCurrentSimple
% Calculation of current density through simple planar tunnel
% barrier
clear all;
clf;
%plotting parameters + fontsizes
FS      = 18;	%label fontsize 18
FSN     = 16;	%number fontsize 16
LW      = 2;	%linewidth

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

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

%*******************
% Input Parameters *
%*******************
% Voltage Bias Parameters
VbiasMin = pi*1e-7;          % Minimum voltage bias (V)
VbiasMax = 0.3;              % Maximum voltage bias (V) 0.3, 2.3
NumVbias = 50;               % Number of samples in voltage range
Vbias = linspace(VbiasMin, VbiasMax, NumVbias);

% Potential
BarrEnergy = 0.3;            % Energy of potential barrier 0.3 (eV)
BarrLength = 5e-9;           % Thickness of potential barrier 5e-9 (m)

% Material Parameters
meff = 0.07;                 % Effective electron mass coefficient
nDopingc = 1e18;             % n-type doping concentration of contacts (cm^-3)

% Sampling Parameters
Numx = 200;                  % Number of samples in x-axis 100
NumE = 500;                  % Number of samples in energy range 500

% Temperature
T = 300;                     % Absolute temperature 300, 4.2 (K)

% Universal Constants
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)

%********************************
% Setup for current calculation *
%********************************
% Change units
m = meff*m0;                 % Effective electron mass (kg)
nDoping = nDopingc*1e6;      % Doping concentration to m^-3
beta = (kB*T)^-1;            % Inverse thermal energy (eV^-1)

% Create x-axis
xMin = -1e-9;
xMax = BarrLength + 1e-9;
x = linspace(xMin, xMax, Numx);
dx = x(2) - x(1);

% Create zero bias potential
Index1 = x >= 0 & x <= BarrLength;
Index2 = x > BarrLength;
V0 = zeros(1, Numx);
V0(Index1) = BarrEnergy;

% Energy range
EFermi = hbar^2*(3*pi^2*nDoping)^(2/3)/(2*m*echarge);
EMin = pi*1e-7;
EMax = EFermi + 15/beta;
Energy = linspace(EMin, EMax, NumE)';
dE = Energy(2) - Energy(1);
E = Energy*ones(1, Numx);

% Chemical potential
mu = CalculateChemicalPotential(T, nDoping, m, EFermi);

% Calculated constants
const0 = echarge^3*m*dE*1e-4/(2*pi^2*hbar^3*beta);
expNum = 1 + exp((mu - Energy)*beta);
constk = sqrt(2*m*echarge)/hbar;

%********************
% Calculate current *
%********************
Current = zeros(NumVbias, 1);
for BiasIndex = 1:NumVbias
    % Create potential
    V = zeros(1, Numx);
    V(Index2) = -Vbias(BiasIndex);
    V(Index1) = V0(Index1) - x(Index1)*Vbias(BiasIndex)/BarrLength;
    V2 = ones(NumE, 1)*V;
    
    % Integration constant
    expDen = 1 + exp((mu - Energy - Vbias(BiasIndex))*beta);
    IntConst = const0*log(expNum./expDen);
 
    k2 = constk*sqrt(E - V2);
    bigP11b = ones(NumE, Numx);
    bigP21b = zeros(NumE, Numx);
    for xIndex = Numx - 1:-1:1 
        ktmp1 = exp(-1i*k2(:, xIndex)*dx);
        ktmp2 = exp( 1i*k2(:, xIndex)*dx);
        kquotient = k2(:, xIndex + 1)./k2(:, xIndex);
        %Setting up transition matrix
        p11b = (1 + kquotient).*ktmp1/2;
        p12b = (1 - kquotient).*ktmp1/2;
        p21b = (1 - kquotient).*ktmp2/2;
        p22b = (1 + kquotient).*ktmp2/2;
        % matrix product
        bigP11b(:, xIndex) = (p11b.*bigP11b(:, xIndex + 1) + p12b.*bigP21b(:, xIndex + 1));
        bigP21b(:, xIndex) = (p21b.*bigP11b(:, xIndex + 1) + p22b.*bigP21b(:, xIndex + 1));
    end
    Trans = abs(1./bigP11b(:, 1)).^2;
    Current(BiasIndex) = sum(IntConst.*Trans);
end

%***************
% Plot Results *
%***************
% Potential
figure(1)
clf
plot(x*1e9, V0, x*1e9, V,'LineWidth',LW)
legend('Vbias = 0 V', ['Vbias = ',num2str(VbiasMax),' V'], 'Location', 'SouthWest')
ttl=['Potential profile, T = ',num2str(T),' K, n = ',...
    num2str(nDoping/1e6,'%7.1e'),' cm^{-3}, \mu = ',...
    num2str(mu,'%6.4f'),' eV'];
title(ttl);
xlabel('Position, x (nm)')
ylabel('Potential, V(x) (eV)')
xlim([xMin*1e9 xMax*1e9])
ylim([1.2*min(V) 1.2*BarrEnergy])

% Current
figure(2)
subplot(1,2,1);
plot(Vbias, Current,'LineWidth',LW);
title('Linear scale current density');
xlabel('Voltage bias, V_{bias} (V)');
ylabel('Current density, J (A/cm^2)');
xlim([0 VbiasMax]);
ylim([0 1.2*max(Current)]);

subplot(1,2,2);
plot(Vbias, log10(Current),'LineWidth',LW);
title('Log_{10} scale current density');
xlabel('Voltage bias, V_{bias} (V)');
ylabel('Current density, log_{10}(J) (A/cm^2)');
xlim([0 VbiasMax]);

ttl3=['n = ',num2str(nDopingc/1e18),'\times10^{18} cm^{-3}, T = ',...
    num2str(T),' K, \mu = ',num2str(mu*1000,'%3.1f'),' meV, L = ',...
    num2str(BarrLength*1e9),' nm, m_{eff} = ',...
    num2str(meff),'\timesm_0'];
% Current
figure(3)
plot(Vbias, Current*1e-5,'LineWidth',LW);
title(ttl3);
xlabel('Voltage bias, V_{bias} (V)');
ylabel('Current density, J (mA/\mum^2)');
xlim([0 VbiasMax]);
ylim([0 1.2*max(Current*1e-5)]);


end

function [Mu] = CalculateChemicalPotential(Temp, NDoping, m, EFermi)
hbar = 1.05457159e-34;                       %Planck's constant (Js)
echarge = 1.6021764e-19;                     %Electron charge (C)
kBoltz = 1.3806505e-23;                      %Boltzmann constant (eV/K)
const = (2.0*m/(hbar*hbar)).^1.5/(2.0*pi*pi); 

if abs(Temp) > 1.0e-10
    beta = 1.0/(kBoltz*Temp);
else
    beta = 1.0/(kBoltz*1.0e-10);
end

MuMin = kBoltz/echarge*Temp*log(NDoping/2.0*((2.0*pi*hbar*hbar)/m*beta).^1.5);
MuMax = EFermi;

EMin = 0.0;
EMax = EFermi + 15.0*kBoltz*Temp/echarge;   %Energy cutoff for given temperature T (K)

nSteps = 10000; %Number of energy points in integral
EStep = (EMin+EMax)/nSteps;                 
rerr = 1.0;

if Temp < 1.0e-10
    Mu = EFermi;
else    
    while rerr > 0.0001
        Mu = MuMin + (MuMax-MuMin)/2.0;
        nprime = 0.0;
        E = EMin;
        for j = 0 : nSteps
            nprime = nprime + sqrt(E*echarge)*1.0/(1.0+exp(-beta*echarge*(Mu-E)));        
            E = E + EStep;
        end
        nprime = const*nprime*EStep*echarge;
        if nprime > NDoping
            MuMax = Mu;
        else
            MuMin = Mu;
            rerr = abs(NDoping-nprime)/NDoping;
        end
    end
end
return
end