clc; clear
couplingNumber = 0.5;
characteristicLength = 4e-3;
microInertia = 0.01;
pres = 4e6;
resFolder = sprintf('Result\\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 = diff(crackPos)./diff(time);
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)
