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


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;

figure
plot(time,disp,'--',time, d_model,'r')
%set(gca, 'FontSize', 18);
legend('experimental data','model displacement')
title('Two-parameter model estimation with OLS');
xlabel('time (s)');
ylabel('displacement');
axis tight;
%axis([0 0.5 -0.00006 0.00006])


figure;
clf;
plot(time,d_res, 'g.');
%axis tight;
axis([0.0 0.5 -0.00006 0.00006]);
ylabel('Residual', 'FontSize', 18);
xlabel('Time', 'FontSize', 18);
grid on;


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

figure
qqplot(d_res)   
ylabel('residuals')