%%%% This part of script comes from the last lecure. If you haven't cleaned it from
%%%% your matlab, no need to run step 1, go to step 2 directly.

%%step 1
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 = 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;

%% step 2: weighted least square
C0 = 1;
K0 = 1500;
CK0 = [C0;K0];
[CK2, cost2] = fminsearch(@cost2_beam, CK0, [], time, disp);
C2 = CK2(1)
K2 = CK2(2)

[c, d_model2] =cost2_beam(CK2, time,disp);
d_res2 = disp - d_model2;
m = length(disp);
sigma22 = sum(d_res2.^2)/m;

%% Comparison
subplot(2,2,1)
plot(d_model)
xlabel('time')
ylabel('fitted value')
title('fitted value vs. time')

subplot(2,2,2)
plot(d_model2)
xlabel('time')
ylabel('fitted value')
title('fitted value vs. time in weighted LS')

subplot(2,2,3)
plot(time,d_res)
xlabel('time')
ylabel('residuals')
title('residual vs. time')

subplot(2,2,4)
plot(time,d_res2)
xlabel('time')
ylabel('residuals')
title('residual vs. time in weighted LS')