
clf;	
load shear_corr.out;
fprintf(1, '<sigma_{xy}(0)sigma_{xy}(0)>/Kb/T/V = %f W/M/K/ps\n', ...
shear_corr(1,2)/(shear_corr(2,1)-shear_corr(1,1)));

plot(shear_corr(:,1), shear_corr(:,2)/shear_corr(1,2));
% hold on; plot(shear_corr(:,1), shear_corr(:,2)/shear_corr(1,2), 'ro');
xlabel ('\tau [ps]');
ylabel ('<\sigma_{xy}(0)\sigma_{xy}(\tau)>/<\sigma_{xy}(0)\sigma_{xy}(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 (['<\sigma_{xy}(0)\sigma_{xy}(\tau)>/<\sigma_{xy}(0)\sigma_{xy}(0)> of Argon at T = ' T 'K, \rho = ' rho ' mol/m^3']); 

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

% Negative data ignored

clf;
% imax=1; 
% while shear_corr(imax,2)>0,  imax=imax+1; end; 
% imax=imax-1;
imax = size(shear_corr(:,1)); imax = imax(1,1);

semilogy (shear_corr(1:imax, 1), shear_corr(1:imax, 2)/shear_corr(1,2));
CHGY;
CHGX;

xlabel ('\tau [ps]');
ylabel ('<\sigma_{xy}(0)\sigma_{xy}(\tau)>/<\sigma_{xy}(0)\sigma_{xy}(0)>')
title (['<\sigma_{xy}(0)\sigma_{xy}(\tau)>/<\sigma_{xy}(0)\sigma_{xy}(0)> of Argon at T = ' T 'K, \rho = ' rho ' mol/m^3']); 

CPRINT;


