%Fig3_9.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}'];

nbarrier=3;						%number of potential barriers
bx=0.4e-9;						%barrier thickness (m)
wx=0.6e-9;						%well width (m)
V0=0.5;							%barrier energy (eV)

N=(2*nbarrier)+1;				%number of samples of potential
for j=1:2:(2*nbarrier)          %set up position and potential array
   dL(j)=wx;   V(j)=0;
   dL(j+1)=bx; V(j+1)=V0;
end
   dL(N)=wx;   V(N)=0;
   
Emin=0.1+(pi*1e-5);				%add (pi*1.0e-5) to energy to avoid divide by zero
Emax=0.5;						%maximum particle energy (eV)
npoints=500;					%number of points in energy plot
dE=(Emax-Emin)/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)
  
    bxmin=0.1e-9;               %minimum central barrier thickness
    bxmax=1.5e-9;               %maximum central barrier thickness
    mpoints=500;                %number of samples
    dbx=(bxmax-bxmin)/mpoints;

    for mcount=1:mpoints         %loop for barrier thickness
    dL(4)=bxmin+(mcount*dbx);    %central barrier thickness
    xp(mcount)=dL(4);            %central barrier thickness vector used in plot
    for j=1:npoints              %loop for energy values                
      E(j)=dE*j+Emin;
      bigP=[1,0;0,1];			%default value of matrix bigP
      	for i=1:N
         k(i)=sqrt(2*echarge*m*(E(j)-V(i)))/hbar;%wave vector at energy E
      	end
      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
      Trans(j,mcount)=(abs(1/bigP(1,1)))^2;     %transmission coefficient
      
    end
    end
ttl=(['\rmFig3.9 (',num2str(nbarrier),'-barriers), \itm\rm^*_e='...
    ,num2str(meff),'\times\itm\rm_0']);
figure(1);
surfl(xp,E,Trans);
axis square;
axis([bxmin bxmax Emin Emax 0 1]);
colormap(copper);
shading interp;
title(ttl);
xlabel('Central barrier thickness (m)');
ylabel('Energy, \itE\rm (eV)')
zlabel('Transmission')
