%Fig2_1.m
%Essential Electron Transport for Device Physics
%Calculation of 1D density of states (dos) from nearest and next nearest
%neigbor tight binding dispersion relation

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);

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


w=emin:estep:emax;      %energy 
t=-1;					%set energy scale - for s-orbital, + for p-orbital
tprime=0.2*t			%tprime
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
kx=0:kx0:kx0*natoms;
dosa=zeros(1,npoints); %initialize dos for nearest neighbors
dosb=zeros(1,npoints); %initialize dos for nearest and next nearest neighbors
	for ix=1:natoms+1
% Nearest neigbor and next nearest neighbor dispersion for the linear lattice
		Ekprime(ix)=(2.0*t*cos(kx(ix)))+(2.0*tprime*cos(2*kx(ix)));
		Ek(ix)=2.0*t*cos(kx(ix));
      for j=1:npoints
dosa(j)=dosa(j)+const/((w(j)-Ek(ix))^2.+gamma2); %Lorentzian broadening
dosb(j)=dosb(j)+const/((w(j)-Ekprime(ix))^2.+gamma2); %Lorentzian broadening
		end
	end
ttl=['Nearest neighbor tight binding dispersion, \gamma =',...
    num2str(gamma),'\timest, n_{atoms}=',num2str(natoms)];   
figure(1);
plot(Ek,'b','LineWidth',LW);
hold on;
plot(Ekprime,'r','LineWidth',LW);
grid on;
title(ttl);
axis([0 natoms 1.1*min(Ekprime) 1.1*max(Ek)]);
xlabel('Wave vector, k_x ( \Deltax \pi/L)'),ylabel('Energy, E_k(t)');
hold off;

figure(2);
plot(w,dosa,'b','LineWidth',LW);
hold on
plot(w,dosb,'r','LineWidth',LW);
title(ttl);
grid on
xlabel('Energy, E (t)'),ylabel('Density of states, N(E)');


	
