function [N,t,Lx,Ly,Lz,vx,vy,vz,T,Ep,Et,P]=readLammps(dumpfilename,logfilename,vcolumns)
% Only works in 3D
% vcolumns = [ivx ivy ivz] column numbers of velocity data
data = readdump_all(dumpfilename);
t = data.timestep;
nt=length(t);
N = data.Natoms;
N = N(1);
% box size
xbox = data.x_bound;
Lx=xbox(1,2)-xbox(1,1);
ybox = data.y_bound;
Ly=ybox(1,2)-ybox(1,1);
zbox = data.z_bound;
Lz=zbox(1,2)-zbox(1,1);
% velocities
vx=squeeze(data.atom_data(:,vcolumns(1),:));
vy=squeeze(data.atom_data(:,vcolumns(2),:));
vz=squeeze(data.atom_data(:,vcolumns(3),:));
logdata = readlog(logfilename);
dum=size(logdata.Chead);
nlast=dum(1);
lastarray=str2num(logdata.data{nlast});
%Step Temp E_pair TotEng Press
T=lastarray(1:nt,2);
Ep=lastarray(1:nt,3);
Et=lastarray(1:nt,4);
P=lastarray(1:nt,5);

