% 
% Written by:
% -- 
% John L. Weatherwax, Ph.D.         07-23-2002
% MIT Lincoln Laboratory
% 244 Wood Street
% Lexington, MA 02420-9176 USA
% 
% email: weatherwax@ll.mit.edu
% Voice: (781) 981-5370       Fax: (781) 981-0721
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

RUN_NUMBER = 4;

switch RUN_NUMBER
 case 1
  qTotal = 6; % Total energy of system
  NA = 3;     % Number of oscillators in solid "A"
  NB = 3;     % Number of oscillators in solid "B"
 case 2
  qTotal = 6; % Total energy of system
  NA = 6;     % Number of oscillators in solid "A"
  NB = 4;     % Number of oscillators in solid "B"
 case 3
  qTotal = 100; % Total energy of system
  NA = 200;     % Number of oscillators in solid "A"
  NB = 100;     % Number of oscillators in solid "B"
 case 4
  qTotal = 80; % Total energy of system
  NA = 100;     % Number of oscillators in solid "A"
  NB = 100;     % Number of oscillators in solid "B"
 otherwise
  error( 'Unknown value of RUN_NUMBER' );
end

% 
% Code begins:
% ---

close all;

qA = 0:qTotal;     % The full          range of q_a.
qB = qTotal - qA;  % The corresponding range of q_b.

nMacrostates = qTotal+1;

omegaA = zeros( 1, nMacrostates );
for i=1:nMacrostates,
  omegaA(i) = nchoosek( qA(i) + NA - 1, qA(i) );
end

omegaB = zeros( 1, nMacrostates );
for i=1:nMacrostates,
  omegaB(i) = nchoosek( qB(i) + NB - 1, qB(i) );
end

omegaTotal = omegaA .* omegaB;

% Plot a bar chart of number of states:
figure; bar( qA, omegaTotal );
xlabel( 'q_A' ); ylabel( '\Omega_{Total}' );

% Plot a bar chart of probability:
figure; bar( qA, omegaTotal/sum( omegaTotal ) );
xlabel( 'q_A' ); ylabel( 'Probability' );



clear i