%  Least squares estimation of parameters for the spring model

%  Load and pre-process data -- for example data1.mat
Cvec=[];
Kvec=[];

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);
d_data = disp;

 
close all;


%  starting guesses for parameter values
C0=1;
K0=1500;
q0=[C0;K0];

%  minimize the least squares function (cost_beam)
[q,cost]=fminsearch(@cost_beam,q0,[],time,d_data);

%  least squares solutions for C and K
C=q(1);
K=q(2);

Cvec=[Cvec; C];
Kvec=[Kvec; K];
%  calculate and plot the fitted model
x0=[d_data(1) 0]; %initial conditions
[t,x]=ode23(@ode_model,time,x0,[],C,K); %solve the spring model
d_model=x(:,1);
dydt=x(:, 2);
figure(3)
plot(time,d_data,'--',time, d_model,'r')
set(gca, 'FontSize', 18);
legend('experimental data','model displacement')
title('Two-parameter model estimation');
xlabel('time (s)');
ylabel('displacement');
%axis tight;
axis([0 0.5 -0.00006 0.00006])