% find dislocation velocity from cfg files

if (~exist('filepath'))
  filepath = '../../runs/ta-edge-10';
  %filepath = '../../runs/ta-edge-100';
  %filepath = '../../runs/ta-71-md-100';
end
nstart = 1;
nend   = 400;
nstep  = 1;
%x0 = 98;      %initial dislocation position
x0 = 74.34;      %initial dislocation position

frameno = nstart:nstep:nend;  nframe = length(frameno);
%t = frameno * 100 * 1e-3; % time in ps
t = frameno * 1000 * 1e-3; % time in ps
xavg = zeros(size(frameno));

for i=1:nframe,
    % load atom configuration
    ifile = frameno(i);
    filename = sprintf('%s/auto%04d.cfg',filepath,ifile);
    [np,data,h]=loadcfg(filename);
    
    % obtain scaled coordinates and cental symmetry paramter
    s = data(:,1:3); csd = data(:,4); 

    % view histogram of central symmetry parameter
    %figure(3); hist(csd);

    % select dislocation core atoms
    ind = find(csd>4.5 & s(:,2) > 0.5 & s(:,2) < 0.516);
    %ind = find(csd>4.5 & csd<10);
    s = data(ind,1:3); r = s*h'; 

    % make sure all atoms are nearest images
    sx = s(:,1); sx = sx - round(sx-sx(1)); s(:,1) = sx; r = s*h';
    
    [sz ind] = sort(s(:,3)); s = s(ind,:); r = r(ind,:);
    
    % compute average dislocation position
    xavg(i) = mean(r(:,1)) - x0;
    if(i>1) xavg(i) = xavg(i) - round((xavg(i)-xavg(i-1))/h(1,1))*h(1,1); end
        
    % plot dislocation position
    figure(1); plot(r(:,3),r(:,1),'.-'); axis equal
    xlim([0 h(3,3)]); ylim([0 h(1,1)]);
    drawnow; pause(0.2);
end

% estimate dislocation velocity 
%  (from second half of data)
nmid = round(nframe*0.7);
vavg = (xavg(nframe) - xavg(nmid)) / ...
       (t(nframe) - t(nmid));  % in A/ps
vavg_in_m_s = vavg * 100;  % in m/s

fs = 12;
figure(2);
plot(t, xavg, '.-', ...
     t(nmid:end), (t(nmid:end)-t(nmid))*vavg+xavg(nmid),'r--');
set(gca,'FontSize',fs);
xlabel('t (ps)');
ylabel('x (angstrom)');
title(sprintf('vavg = %.1f (m/s)',vavg_in_m_s));