%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