% 
% Written by:
% -- 
% John L. Weatherwax, Ph.D.         2002-09-08
% 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.
% 
%-----

% 
% Code begins:
% --

close all;

N = 50;              % The number of ossilators.
q = 0:100;

nQ = length( q );

% Get omega:
omega = zeros( 1, nQ  );
for i=1:nQ,
  % omega = (exp(1)*q/N).^N; % The high energy limit ...
  omega(i) = nchoosek( q(i)+N-1, q(i) );
end
clear i;

sOverK = log( omega );

kTOverEpsilon = central_diff( q, sOverK );

figure; plot( q, sOverK, '-o' ); grid on;
title( 'S v.s. q' );

cOverNK = central_diff( q, kTOverEpsilon )/N;

figure; plot( kTOverEpsilon, cOverNK, '-o' ); grid on;
title( 'C v.s. T' );