clc; clear
couplingNumber = 0.0;
characteristicLength = 4e-3;
microInertia = 0.01;
pres = 2e6;
resFolder = sprintf('Result_pr_1_3\\CN_%g_LB_%g_J0_%g_S0_%g\\', couplingNumber, ...
    characteristicLength*1e3, microInertia, pres/1e6);
load([resFolder, 'init'])
totStep = init.nrecord;
timeSteps = 1:1:totStep;
x = init.x;
y = init.y;

time = zeros(size(timeSteps));
crackPos = time;
for i = 1:numel(timeSteps)
    step = timeSteps(i);
    stepFile = [resFolder,sprintf('step%d.mat',step)];
    if(~exist(stepFile, 'file'))
        break;
    end
     load(stepFile);
%     h = figure(200); clf
%     set(h, 'Position', [100, 100, 600, 320*3/4])
    time(i) = step.currentTime;
    dmg = step.dmg(:);
    crackPoints = dmg>0.4;
    crackPos(i) = max(x(crackPoints));
end
% speed = conv(crackPos, [-1,-2,-3,-4,4,3,3,1], 'endpoints', 'no')
% speed = diff(crackPos)./diff(time);
plot(time, crackPos ,'r--')
crackPos = smooth(crackPos, 80);

dt = 0.5e-7;
ngp = 10;
%time0 = zeros(numel(time)-2*ngp+2,1);
for i = ngp:numel(crackPos)-ngp+1
   speed(i-ngp+1) = 0;
   time0(i-ngp+1) = time(i-ngp+1);
   for j = 1:ngp-1
     speed(i-ngp+1) =  speed(i-ngp+1) + (crackPos(i+j)-crackPos(i-j)) / (2*j*dt);
    
%     time0(i-ngp+1) = time0(i-ngp+1) + (time(i+j)+time(i-j));
   end
    speed(i-ngp+1) = speed(i-ngp+1) / (ngp-1);
%    time0(i-ngp+1) = time0(i-ngp+1)/(2*(ngp-1)+1);
end

plot(time0([1:1:end,end])*1e6, speed([1:1:end, end])/3102 )
axis([0, 40, 0, 0.72])


% init.dens = 2440;
% microInertia = init.microInertia;
% dens = 2440;
% lambda = init.nu * init.emod / (1+init.nu) / (1-2*init.nu);
% mu_r = init.mu_r;
% mu = init.mu;
% gamma = init.gamma;
% c1 = sqrt((lambda+2*mu+2*mu_r)/dens)
% c2 = sqrt((mu+mu_r)/dens)
% c3 = sqrt(gamma/microInertia/dens)
% c4 = sqrt(4*mu_r*microInertia/dens/microInertia)
% c0 = sqrt(max(lambda+2*init.mu+2*init.mu_r, init.mu+init.mu_r)/init.dens);
% figure(300); clf
% plot((time(2:end)+time(1:end-1))/2*1e6, speed/3102)
% time0 = (time(1:end-1)+time(2:end) )/ 2;
timeSpeed = [time0'*1e6, speed'/3102];
