U = load('data1');
t = U.trace_x';
disp = U.trace_y(4,:);
N = length(t);
time = t;
[maxd,index] = max(disp);
time = time(index:N)-time(index);
disp = disp(index:N);
disp = disp';
disp = disp - mean(disp);
plot(time, disp);
xlabel('t');
ylabel('y(t)');

C0 = [0:.05:2];
K0 = [1000:2:2000];
LCK = zeros(length(C0), length(K0));
for i = 1:length(C0)
    for j = 1:length(K0)
        CK = [C0(i);K0(j)];
        [c,d] = cost_beam(CK,time,disp);
        LCK(i,j) = c;
    end
end

surf(C0,K0, LCK')
xlabel('C')
ylabel('K')

C0 = 1;
K0 = 1500;
CK0 = [C0;K0];
[CK, cost] = fminsearch(@cost_beam, CK0, [], time, disp);
C = CK(1)
K = CK(2)

[c, d_model] =cost_beam(CK, time,disp);
d_res = disp - d_model;
m = length(disp);
sigma2 = sum(d_res.^2)/m;


plot(time,d_res)
xlabel('time')
ylabel('residuals')

plot(d_model,d_res,'k.')
xlabel('fitted')
ylabel('residuals')

qqplot(d_resid,'g*')
ylabel('residuals')