% Two Einstein solids NA & NB oscillators
% in contact to exchange energy
NA = 20;   % Number of oscillators in crystal A
NB = 10;   % Number of oscillators in crystal B
q = 30;    % qA+qB total energy of crystals A and B
% Find multiplicity of all macrostates
nstate = q+1;   % number of states 
qA = 0:q;       % all possible values of qA
qB = q - qA;    % all possible values of qB
OmegaA = zeros(nstate,1);
OmegaB = OmegaA;
FontSize=20;
% Loop through all macrostates and find multiplicity
for i = 1:nstate
    OmegaA(i) = nchoosek(qA(i)+NA-1,qA(i)); %binomial 
    OmegaB(i) = nchoosek(qB(i)+NB-1,qB(i)); %binomial
end
OmegaTOT= OmegaA.*OmegaB;   % Total multiplicity
subplot(1,2,1), plot(qA,OmegaTOT/max(OmegaTOT),'-o');
subplot(1,2,2), semilogy(qA,OmegaTOT/max(OmegaTOT),'-o');
% pretty plotting...
maxO=sprintf('%0.3g',max(OmegaTOT));
title(['max(\Omega_{Tot}) = ',maxO],'FontSize',FontSize)
for i=1:2
    subplot(1,2,i)
    xlabel('q_A','FontSize',FontSize)
    ylabel('\Omega_{Tot}(N_A, N_B, q_A) / max(\Omega_{Tot})',...
        'FontSize',FontSize)
    h=gca;
    h.FontSize=FontSize;
end

%

P=OmegaTOT/sum(OmegaTOT);
subplot(1,2,1), plot(qA,P,'-o');
subplot(1,2,2), semilogy(qA,P,'-o');
% pretty plotting...
FontSize=20;
for i=1:2
    subplot(1,2,i)
    xlabel('q_A','FontSize',FontSize)
    ylabel('P(N_A, N_B, q_A)','FontSize',FontSize)
    h=gca;
    h.FontSize=FontSize;
end
[Pmax,qmax]=max(P);
subplot(1,2,1), text(60, Pmax+0.01,['q_A(P_{max})=',num2str(qmax)],'FontSize',FontSize-2)