%Fig2_10.m
%Essential Electron Transport for Device Physics
%Plotting parameters + fontsizes
clear all;
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=['\fontname{MT Extra}h\fontname{Arial}'];
%Energy levels equally spaced harmonic potential
%calculate combinations for three particles with Pauli exclusion
%Etot = total actual energy + num_particles/2 hbar omega (for convenience)
%example: lowest energy for three particles is 0.5+1.5+2.5=4.5hw
%this corresponds to Etot=4.5hw+3/2hw=6
%Etot=15 correpsonds to 13.5hw
%Etot=50 corresponds to 48.5hw
nplots=3;%3 number of plots
Etotmin=6;%6
Etotstep=4;%4, 14, 24, 34
Etotmax=Etotmin+nplots*Etotstep;
Etot=Etotmin-Etotstep;

for ij=1:nplots
    Etot=Etot+Etotstep;%total actual energy + (3 particles x 1/2 hbar omega) (for convenience)
N=Etot;
 Etotv(ij)=Etot
%Energy=zeros(1,N*(N-1));
Counter=zeros(1,N*(N-1));
Emax=N+(N-1)+(N-2);%for three particles
A=zeros(Emax,N);
for i=1:N
    for j=1:N-1
        if i ~= j
            for k=1:N-2
                if (k~=j) && (k~=i)
                Energy=i+j+k;
                Counter(1,Energy)=Counter(1,Energy)+1;%number of configs for E
                A(Energy,i)=A(Energy,i)+1;
                A(Energy,j)=A(Energy,j)+1;
                A(Energy,k)=A(Energy,k)+1;
                end
            end
        end
    end
end
%**************************************************************************
s=char('y','k','r','g','b','m','c');	%plot curves in different colors
 j=1+mod(ij,7); 
x=linspace(1,N,N);
pttl=(['-o',s(j)]);
ttl=(['\rmFig2.10, 3 particles, E_{tot}^{min}=',num2str(Etotmin-1.5),hbar,...
    '\omega, E_{tot}^{max}=',num2str(Etotmax-Etotstep-1.5),hbar,'\omega']);
figure(1)
plot(x-0.5,(A(Etot,:)'/Counter(Etot)),pttl,'MarkerFaceColor',s(j),'LineWidth',LW)
hold on
grid on
xlabel('Energy loss, $\hbar\omega$ (meV)', 'Interpreter', 'latex','FontSize',FS);
%xlabel(['Energy, E (',hbar,'\omega)'],'FontSize',FS);
ylabel('Occupation probability','FontSize',FS);
title (ttl);
ymax=max(A(Etot,:)'/Counter(Etot));
axis([0 Etot 0 1.1]);

end
hold off
