%Fig2_2a.m
%Essential Electron Transport for Device Physics
%Calculation of 2D square lattice density of states (dos) from nearest 
%neigbor tight binding dispersion relation

clear all;
clf;
%plotting parameters + fontsizes
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);

natoms=100;				%number of lattice points (atoms) in each direction
npoints=100;			%number of points in plot
emax=6.5;				%max energy in units of t
emin=-6.5;				%min energy in units of t
if emax<=emin;error('Error: require emax > emin');end;
   estep=(emax-emin)/(npoints-1);

dos=zeros(1,npoints);	%set dos vector to zero
w=emin:estep:emax;		%energy 
t=-1;					%set energy scale - for s-orbital, + for p-orbital
gamma = 0.1;			%Energy broadening in units of t
gamma=gamma*abs(t);     %use units of t
gamma2=(gamma/2)^2;
const=gamma/pi;

kx0=pi/natoms;			%kx constant
ky0=pi/natoms;			%ky constant
kx=0:kx0:kx0*natoms;
ky=0:ky0:ky0*natoms;
	for ix=1:natoms+1
	for iy=1:natoms+1
% Nearest neigbor dispersion for the square lattice
		Ek(ix,iy)=2.0*t*(cos(kx(ix))+cos(ky(iy)));
      for j=1:npoints
			dos(j)=dos(j)+const/(((w(j)-Ek(ix,iy))^2.)+gamma2); %Lorentzian broadening
		end
	end
	end
figure(1);
plot(w,dos,'b','LineWidth',LW);
ttl=['\rmFig2.2a, nearest-neighbor tight-binding density of states, \Gamma=',...
    num2str(gamma),'\times|t|, n_{atoms}=',num2str(natoms)];
title(ttl);
xlabel('Energy, E (t)'),ylabel('Density of states, N(E)');

figure(2);
surf(Ek);
ttl=['\rmFig2.2a, nearest-neighbor tight-binding dispersion, \Gamma=',...
    num2str(gamma),'\times|t|, n_{atoms}=',num2str(natoms)];
shading interp;
title(ttl);
xlabel('Wave vector, k_x ( \Deltax \pi/L)');
ylabel('Wave vector, k_y (\Deltax \pi/L)');
zlabel('E_k');