clf; clear;
name = {
'z', '\sigma'
'count', ''
'\rho', '\sigma^{-3}'
'v_x', '(\epsilon/m)^{1/2}'
'T', 'K'
'\tau_{xz}', '\epsilon/\sigma^{3}'
'P', '\epsilon/\sigma^{3}'
}; n = size (name,1);

fn = 'me2.out'; fid = fopen (fn, 'r');
kb=fscanf(fid,'%d',1); 
kstart=fscanf(fid,'%d',1); 
delta=fscanf(fid,'%f',1);
UENERGY_IN_J=fscanf(fid,'%f',1);
ULENGTH_IN_A=fscanf(fid,'%f',1);
ARGON_MASS_IN_KG=fscanf(fid,'%f',1);
nmb=fscanf(fid,'%d',1);
H=fscanf(fid,'%f',3);
nc=fscanf(fid,'%f',4);
df=fscanf(fid,'%f',4);
cf=fscanf(fid,'%f',2);
zwall = H(3) * nc(4)/nc(3)/2.;
zopc  = H(3) * df(3); ztfe  = H(3) * df(4);
me = fscanf(fid, '%f', [n inf]); me = me'; fclose(fid);

c = 1; d = 3; clf;

plot (me(:,c), me(:,d), 'o'); hold on;
plot (me(:,c), me(:,d));
xlabel([name{c,1} ' [' name{c,2} ']']);
ylabel([name{d,1} ' [' name{d,2} ']']);

title (['Density profile of ' num2str(nmb) ' Ar atoms ' ...
      'avg. over ' num2str(kb-kstart) ' steps (' ...
      num2str(delta*(kb-kstart)) ' ps)']);

v = axis; plot ([zopc zopc], [v(3) v(4)], 'w-.');
axis ([0 H(3) 0 0.6]);
text(16, 0.1, 'Boundary OPC');

print -deps density.eps;
