clf;	
load corr.out;
fprintf(1, '<J(0)J(0)>/Kb/T^2/V = %f W/M/K/ps\n', ...
corr(1,2)/(corr(2,1)-corr(1,1)));

% plot(corr(:,1), corr(:,2)/corr(1,2), 'ro');
% hold on; 
plot(corr(:,1), corr(:,2)/corr(1,2));
xlabel ('\tau [ps]');
ylabel ('<J(0)J(\tau)>/<J(0)J(0)>')
v = axis;
axis ([v(1) v(2) 0 1]);
CHGY;
CHGX;

if exist('T'), 
 a = input(['T = ' T ' K? '], 's');
 if ~isempty(a), T = a; end;
else
 T = input('T = ', 's');
end;

if exist('rho'), 
 a = input(['density = ' rho ' mol/m^3? '], 's');
 if ~isempty(a), rho = a; end;
else
  rho = input('density(mol/m^3) = ', 's');
end;

title (['<J(0)J(\tau)>/<J(0)J(0)> of Argon at T = ' T 'K, \rho = ' rho ' mol/m^3']); 

CPRINT;
a = input('press a key..', 's');

clf;
imax=1; 
while corr(imax,2)>0,  imax=imax+1; end; 
imax=imax-1;
semilogy (corr(1:imax, 1), corr(1:imax, 2)/corr(1,2), 'o');
xlabel ('\tau [ps]');
ylabel ('<J(0)J(\tau)>/<J(0)J(0)>')
title (['<J(0)J(\tau)>/<J(0)J(0)> of Argon at T = ' T 'K, \rho = ' rho ' mol/m^3']); 

CPRINT;
