%Fig1_10.m
%Essential Electron Transport for Device Physics
%Electron in periodic potential with rectangular barriers
%real (blue) and imaginary (red) band structure
clear all
clf;		
%plotting parameters + fontsizes
FS      = 12;	%label fontsize
FSN     = 12;	%number fontsize
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);

hbarChar=['\fontname{MT Extra}h\fontname{Times}'];

npoints=1750;			%number of points in energy plot 1750
hbar=1.0545715e-34;		%Planck's constant (Js)
echarge=1.6021764e-19;	%electron charge (C)
m0=9.109382e-31;		%bare electron mass (kg)

V0eV=3.0;               %barrier energy [eV] 3.0
V0=V0eV*echarge;        %barrier energy [J]
Lb=0.4e-9;              %barrier thickness [m] 0.4
Lw=0.6e-9;              %well thickness [m] 0.6
L=Lb+Lw;                %cell period [m]

k1=linspace(0,3.5*pi/L,npoints);%range of k-values considered determines E
E=(k1.^2)*(hbar^2)/(2*m0);      %electron energy, E [J]
EeV=E/echarge;max(EeV)
k2=sqrt(2*m0*(E-V0))/hbar; %k2 real for E<V0

k12=k1.^2;
k22=k2.^2;
Aconst=(k12+k22)./(2*k1.*k2);

theta=(cos(k2*Lb).*cos(k1*Lw))-(Aconst.*sin(k2*Lb).*sin(k1*Lw));

kLre=zeros(npoints);    %initialize real band structure
kLim=zeros(npoints);    %initialize imaginary band structure

if any(abs(theta) < 1.0)    
    kLre=real((acos(theta))/pi);
end

if any(abs(theta) >= 1.0)   
    kLim=real((acosh(theta))/pi);
end
%**************************************************************************
%plot figure
%**************************************************************************

ttl=['\rmFig1.10, V_0 = ',num2str(V0eV,'%3.1f'),' eV, L_b = ',...
 num2str(Lb*1e9,'%3.1f'),' nm, L_w = ',...
 num2str(Lw*1e9,'%3.1f'),' nm, L = ',num2str((Lw+Lb)*1e9,'%3.1f'),' nm'];
figure(1)
plot(kLim, EeV, 'r','LineWidth',LW);
hold on;
plot(kLre, EeV, 'b','LineWidth',LW);
axis([0,1.0,0,max(EeV)]);%/3 /5
grid on
xlabel('Wave vector, kL/\pi');
%yttl=['Energy, E (',hbarChar,'^2/2m_0)'];
yttl=['E (eV)'];
ylabel(yttl);
hold off;

title (ttl);