%% Simple demo of Boltzmann distribution

%Each particle can have energy 0,1E0, 2E0, ...9E0. 

E0=1;
N=6;
Etot=5*E0;

%intialize counters to zero
Num_E=zeros(10,1);
Omega=0;

for ss=0:10^N-1
    
    for nn=1:N  %put 'digits' representing energy into a vector
        E_State_Vector(N-nn+1)=mod(floor(ss/10^(nn-1)),10);
    end
    
    if sum(E_State_Vector)==Etot %State matches correct total energy
        E_State_Vector
        for ii=1:N %count up occurances of the individual particle energies
            Num_E(E_State_Vector(ii)+1)=Num_E(E_State_Vector(ii)+1)+1;
        end
        Omega=Omega+1;%keep track of total number of states with Etot
    end
end

S_divKb5=log(Omega)

figure
Eparticle=(0:9);
plot(Eparticle,Num_E/sum(Num_E),'linewidth',2)
xlabel('Energy of particle (E/E0)')
ylabel('Probabitily')
legend('Observed Energy Distribution')


%% Increase total energy, measure change in Entropy

E0=1;
N=6;
Etot=6*E0;

%intialize counters to zero
Num_E=zeros(10,1);
Omega=0;

for ss=0:10^N-1
    
    for nn=1:N  %put 'digits' representing energy into a vector
        E_State_Vector(N-nn+1)=mod(floor(ss/10^(nn-1)),10);
    end
    
    if sum(E_State_Vector)==Etot %State matches correct total energy
        E_State_Vector
        for ii=1:N %count up occurances of the individual particle energies
            Num_E(E_State_Vector(ii)+1)=Num_E(E_State_Vector(ii)+1)+1;
        end
        Omega=Omega+1;%keep track of total number of states with Etot
    end
end

S_divKb6=log(Omega)

figure
Eparticle=(0:9);
plot(Eparticle,Num_E/sum(Num_E),'linewidth',2)
xlabel('Energy of particle (E/E0)')
ylabel('Probabitily')
legend('Observed Energy Distribution')


%% Make and Estimate of Temperature and Compare to Boltzmann Distribution

%T=dS/dE

%make a crude calculation of dS/dE
KbT=1/(S_divKb6-S_divKb5)

figure
Eparticle=(0:9);
plot(Eparticle,Num_E/sum(Num_E),'linewidth',2)
hold on
plot(Eparticle,exp(-Eparticle/KbT)/sum(exp(-Eparticle/KbT)),'r','linewidth',2)
legend('Observed  Energy Distribution','Boltzmann Distribution')
xlabel('Energy of particle (E/E0)')
ylabel('Probabitily')






