%Fig3_7.m
%Essential Electron Transport for Device Physics
%resonant tunneling through 3 barriers

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);
hbar_txt=['\fontname{MT Extra}h\fontname{Arial}'];

t=cputime;
nbarrier=3;					%number of potential barriers
bx=0.4e-9;					%barrier width (m)
wx=0.6e-9;					%well width (m)
V0=0.5;						%barrier energy (eV)
   
N=140+1;                    %number of samples of potential
for j=1:1:20                %set up position and potential array
   dL(j)=wx/5;         V(j)=0;
   dL(j+20)=bx/20;      V(j+20)=V0;   %barrier 1
   dL(j+40)=wx/20;      V(j+40)=0;
   dL(j+60)=bx/20;      V(j+60)=V0;   %barrier 2
   dL(j+80)=wx/20;      V(j+80)=0;
   dL(j+100)=bx/20;     V(j+100)=V0;  %barrier 3
   dL(j+120)=wx/5;     V(j+120)=0;   
end
dL(N)=wx/20;   V(N)=0;
xp(1)=dL(1);
for j=2:N;xp(j)=xp(j-1)+dL(j);end;   %set up x-plot axis
Emin=pi*1e-5;				%add (pi*1.0e-5) to energy avoid divide by zero
Emax=0.8;					%maximum particle energy (eV)
npoints=1500;				%number of points in energy plot
dE=Emax/npoints;			%energy increment (eV)
hbar=1.0545715e-34;			%Planck's constant (J s)
eye=complex(0.,1.);			%square root of -1
m0=9.109382e-31;			%bare electron mass (kg)
meff=1.0;					%effective electron mass / m0
m=meff*m0;					%effective electron mass (kg)
echarge=1.6021764e-19;		%electron charge (C)

% vectorized form of outer for-loop below, computing all energylevels
% simultaneously
energyindex=[1:npoints];
E=Emin+dE*energyindex;
E2=E.'*ones(1,N);
V2=ones(npoints,1)*V;
npointsx2=2*npoints;
bigP2=speye(npointsx2); 
k2=sqrt(2*echarge*m*(E2-V2))/hbar;
for n=1:(N-1)
   p11=0.5*(1+k2(:,n+1)./k2(:,n)).*exp(-eye*k2(:,n)*dL(n));
   p12=0.5*(1-k2(:,n+1)./k2(:,n)).*exp(-eye*k2(:,n)*dL(n));
   p21=0.5*(1-k2(:,n+1)./k2(:,n)).*exp(eye*k2(:,n)*dL(n));
   p22=0.5*(1+k2(:,n+1)./k2(:,n)).*exp(eye*k2(:,n)*dL(n));
   p2=sparse(2*energyindex-1,2*energyindex-1,p11,npointsx2,npointsx2);
   p2=p2+sparse(2*energyindex-1,2*energyindex,p12,npointsx2,npointsx2);
   p2=p2+sparse(2*energyindex,2*energyindex-1,p21,npointsx2,npointsx2);
   p2=p2+sparse(2*energyindex,2*energyindex,p22,npointsx2,npointsx2);
   bigP2=bigP2*p2;
end
Trans2=diag(bigP2);
Trans2=Trans2(2*energyindex-1);
Trans2=abs(1./Trans2).^2;
%**************************************************************************
%Plot potential and transmission coefficient
%**************************************************************************
figure(1);                %plot potential 
subplot(1,2,1)
plot(xp*1e9,V,'LineWidth',LW)
axis([0,xp(end)*1e9,0,Emax])
xlabel('Position, x (nm)')
ylabel('Potential energy, V (eV)')
ttl=['\rmFig3.7, m^*_e=',...
    num2str(meff),'\timesm_0, L_b=',num2str(bx*1e9),...
    ' nm, L_w=',num2str(wx*1e9),' nm'];
title(ttl);
subplot(1,2,2),plot(Trans2,E,'LineWidth',LW);
axis([0,1,0,Emax]);
xlabel('Transmission coefficient');
ylabel('Energy, E (eV)');
  
figure(2);                %plot transmission coefficient
subplot(1,2,1),plot(xp*1e9,V,'LineWidth',LW);
axis([0,xp(end)*1e9,0,Emax]);
xlabel('Position, x (nm)');
ylabel('Potential energy, V (eV)');
title(ttl);
subplot(1,2,2);plot(-log(Trans2),E,'LineWidth',LW);axis([0,20,0,Emax]);
xlabel('-ln trans. coeff.');ylabel('Energy, E (eV)');

%**************************************************************************
%Calculate electron wave function
%**************************************************************************
E1=0.3218;              %fixed electron energy (eV), 0.3218 0.4016

bigP=[1,0;0,1];			%default value of matrix bigP
for i=1:N
 k(i)=sqrt(2*echarge*m*(E1-V(i)))/hbar;%wave vector at energy E1
end
% vectorized construction of k
k2=sqrt(2*echarge*m*(E1-V))/hbar;
for n=1:(N-1)						%multiply out propagation matrix
 p(1,1)=0.5*(1+k(n+1)/k(n))*exp(-eye*k(n)*dL(n));
 p(1,2)=0.5*(1-k(n+1)/k(n))*exp(-eye*k(n)*dL(n));
 p(2,1)=0.5*(1-k(n+1)/k(n))*exp(eye*k(n)*dL(n));
 p(2,2)=0.5*(1+k(n+1)/k(n))*exp(eye*k(n)*dL(n));
 bigP=bigP*p;
end

C(1)=1/bigP(1,1);                    %transmission amplitude
B(1)=bigP(2,1)*C(1);                 %reflection amplitude
A(1)=1;                              %incident amplitude

bigP=[1,0;0,1];                      %default value of matrix bigP
for n=1:(N-1)                        %multiply out propagation matrix
 p(1,1)=0.5*(1+k(n+1)/k(n))*exp(-eye*k(n)*dL(n));
 p(1,2)=0.5*(1-k(n+1)/k(n))*exp(-eye*k(n)*dL(n));
 p(2,1)=0.5*(1-k(n+1)/k(n))*exp(eye*k(n)*dL(n));
 p(2,2)=0.5*(1+k(n+1)/k(n))*exp(eye*k(n)*dL(n));
 bigP=bigP*p;
 Prop=inv(bigP)*[A(1);B(1)];
 C(n)=Prop(1);
 D(n)=Prop(2);
 A(n+1)=C(n);
 B(n+1)=D(n);
end
for n=1:(N-1)
  psi(n)=(C(n)+D(n));          %calculate wave function
end
psi(N)=psi(N-1);               %dummy value for N-th value of wave function

psi2=conj(psi).*psi;           %calculate wave function squared
%**************************************************************************
%Plot modulus of wave function squared
%**************************************************************************
figure(3);
plot(xp*1e9,psi2,'LineWidth',LW);
xlabel('Position, x (nm)'),ylabel('|\psi(x)|^2');
ttl2=['\rm',num2str(nbarrier),'-barrier, m^*_e=',...
    num2str(meff),' \times m_0, L_b=',num2str(bx*1e9),...
    ' nm, L_w=',num2str(wx*1e9),' nm, E=',num2str(E1),' eV'];
title(ttl2);
%**************************************************************************
%Plot real and imaginary part of wave function
%**************************************************************************
figure(4);
plot(xp*1e9,real(psi),'b','LineWidth',LW); 
hold on
plot(xp*1e9,imag(psi),'r','LineWidth',LW); 
xlabel('Position, x (nm)'),ylabel('Real (blue), Imag (red), (\psi(x))');
title(ttl2);      
hold off
e=cputime-t;
e