%%%%%%%%%%%%%%%%%%%%%%%%%
%    Function Pac       %
%%%%%%%%%%%%%%%%%%%%%%%%%
function pac(i,j)
global dat T D;

clf;
if exist('dat')==0
load dat.out;
end

if exist('i')<5 T=i; end
if exist('j')<5 D=j; end
ss=input('T = ','s');
if isstr(ss) T = ss; end;
ss=input('Density = ','s');
if isstr(ss) D=ss; end;

minx = min(dat(:,2));
maxx = max(dat(:,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT THE MSD EVOLUTION %
%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,1);
plot(dat(:,2),dat(:,7));
title(['MSD at T = ' T 'K  Density = ' D]);
xlabel('ps');
CHG_AXIS(minx,maxx);

%%%%%%%%%%%%%%%%%%%%%%%%%
%   PLOT THE PRESSURE   %
%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,2);
plot(dat(:,2),dat(:,5));
title('Pressure');
ylabel('Bar');
xlabel('ps');
CHG_AXIS(minx,maxx);

%%%%%%%%%%%%%%%%%%%%%%%%%%
%  THE POTENTIAL ENERGY  %
%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,3);
plot(dat(:,2),dat(:,4));
title('Potential energy');
xlabel('ps');
ylabel('kJ/mol');
CHG_AXIS(minx,maxx);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TOTAL ENERGY CONSERVATION %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,4);
plot(dat(:,2),dat(:,6));
title('Total energy conservation');
xlabel('ps');
ylabel('kJ/mol');
CHG_AXIS(minx,maxx);

%%%%%%%%%%%%%
% PRINT OUT %
%%%%%%%%%%%%%
PRINT;

end;
