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

clear
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/1000;         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/20;     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
   xpmax=max(xp);

Emax=0.6;			%maximum particle energy (eV)
npoints=141;		%number of points in energy plot
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)
% 
%**************************************************************************
%Calculate electron wave function
%**************************************************************************
%loop energy 

E1=[Emax/npoints:Emax/npoints:Emax]; %set up energy array of npoints
psi2=zeros(npoints,N);%initialize psi2 to zero
C=zeros(1,N);D=zeros(1,N);
for j=1:npoints
    
 bigP=[1,0;0,1];					%default value of matrix bigP
      for i=1:N
        k(i)=sqrt(2*echarge*m*(E1(j)-V(i)))/hbar;%wave vector at energy E1
      end
     for n=1:(N-1)					%multiply out propagation matrix
         a=(1+k(n+1)/k(n));
         b=(1-k(n+1)/k(n));
         c=(exp(-eye*k(n)*dL(n)));
         d=(exp(eye*k(n)*dL(n)));
         p(1,1)=a*c;
         p(1,2)=b*c;
         p(2,1)=b*d;
         p(2,2)=a*d;
         bigP=0.5*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
         a=(1+k(n+1)/k(n));
         b=(1-k(n+1)/k(n));
         c=(exp(-eye*k(n)*dL(n)));
         d=(exp(eye*k(n)*dL(n)));
         p(1,1)=a*c;
         p(1,2)=b*c;
         p(2,1)=b*d;
         p(2,2)=a*d;
         bigP=0.5*bigP*p;
         Prop=inv(bigP)*[A(1);B(1)];
         C(n)=Prop(1);
         D(n)=Prop(2);
      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(j,:)=conj(psi).*psi;     %calculate modulus wave function squared
end
psi2max=max(max(psi2));       %find maximum value of psi2
%**************************************************************************
%Plot modulus of wave function squared as function of position and energy
%**************************************************************************
ttl1=['\rmFig3.8 (',num2str(nbarrier),'-barriers), m^*_e='...
    ,num2str(meff),'\timesm_0'];
figure(1)%surface plot with lighting
surfl(xp,E1,psi2);
shading interp;
colormap(copper);
axis square;
xlabel('Position, x (m)'),ylabel('Energy, E (eV)'),zlabel('|\psi(x)|^2');

title(ttl1);
axis([ 0 xpmax 0 Emax 0 psi2max]);

ttl2=['\rmFig3.8, -ln(|\psi|^2) contour (',num2str(nbarrier),...
    '-barriers), m^*_e= ',num2str(meff),'\timesm_0'];
figure(2)%contour plot of -ln(psi2)
contour(xp,E1,-log(psi2),40);
axis square;
xlabel('Position, x (m)'),ylabel('Energy, E (eV)'),zlabel('|\psi(x)|^2');
title(ttl2);
axis([ 0 xpmax 0 Emax]);

e=cputime-t;
e
