% Fig_1_14.m
% Essential Semiconductor Laser Device Physics
% Kronig-Penney band structure
% electron in periodic potential with rectangular barriers
% real and imaginary band structure
clear all
clf;		
%plotting parameters + fontsizes
clear all;
clf;

FS      = 28;	%label fontsize 18
FSN     = 28;	%number fontsize 16
LW      = 2;	%linewidth 2

% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times New Roman')
set(0,'DefaultAxesFontSize', FSN)

% Change default text fonts.
set(0,'DefaultTextFontname', 'Times New Roman')
set(0,'DefaultTextFontSize', FSN)

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

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=1.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=['\itV\rm_0 = ',num2str(V0eV,'%3.1f'),' eV, \itL\rm_b = ',...
 num2str(Lb*1e9,'%3.1f'),' nm, \itL\rm_w = ',...
 num2str(Lw*1e9,'%3.1f'),' nm, \itL\rm = ',...
 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, \itkL\rm/\pi');
%yttl=['Energy, E (',hbarChar,'^2/2m_0)'];
yttl=['Energy, \itE_k\rm (eV)'];
ylabel(yttl);
hold off;

title (ttl);