data = readdump_all('dump.lammpstrj');
% data.atom_data (time) = atom_id atom_type x y z 
t = data.timestep;
N = data.Natoms;
nt = length(t);
nleft = zeros(nt,1);
xbox = data.x_bound;
Lx=xbox(1,2)-xbox(1,1);
ybox = data.y_bound;
Ly=ybox(1,2)-ybox(1,1);
rho=N(1)/(Lx*Ly);
halfsize = 0.5*xbox(:,2);
for it = 1:nt
    xit = data.atom_data(:,3,it);
    jj = find(xit<halfsize(it));
    nleft(it) = length(jj);
end
figure(1), plot(t,nleft), xlabel('t'), ylabel('n')
std(nleft)
mean(nleft)

vx=squeeze(data.atom_data(:,6,:));
vy=squeeze(data.atom_data(:,7,:));
vxsq=mean(vx.^2);
vysq=mean(vy.^2);
T=(vxsq+vysq)/2;
figure(2), plot(t,vxsq,'-r',t,vysq,'-b',t,T,'-k')
xlabel('t'), ylabel('T')

Pkin=rho*T;
figure(3), plot(t,Pkin,'-k'), xlabel('t'), ylabel('Kinetic pressure')