
% fitting Y = A+B*cos(C*X+D).*exp(E*X);
% to the trajectory of a damped oscillator

format compact
options=optimset('MaxFunEvals',100000,'MaxIter',1000);

A=0.7; B=0.15; C=5; D=1.3; E=-0.004;
A =0.7; B = 0.15; C = 5.27; D = 1.22; E = -0.0037;
param0=[A,B,C,D,E];


load -ascii pendulumCR.dat
X=pendulumCR(:,1);
Y1=pendulumCR(:,2);
Y2=pendulumCR(:,3);

param0=[-0.0003,15.5343,5.2766,1.2207,-0.0037];
Px=lsqcurvefit('damped_oscill',param0,X,Y1,[],[],options)
param0=[-0.0011,80.3684,5.2766,1.2207+pi/2,-0.0037];
Pv=lsqcurvefit('damped_oscill',param0,X,Y2,[],[],options);

omega=C
gamma_o_m=E

t0=pendulumCR(1,1);
t0=0;
t1=X(length(X));
T=t0:(t1-t0)/1000:t1;

figure(1)

subplot(2,1,1)
plot(X,Y1,'-',T,damped_oscill(Px,T))
axis tight

subplot(2,1,2)
plot(X,Y2,'-',T,damped_oscill(Pv,T))
axis tight


figure(2)

span=1:2:length(X);


tmax=30;

subplot(2,1,1)
plot(X(span),Y1(span),'o',T,damped_oscill(Px,T))
axis tight;v=axis;axis([0 tmax v(3) v(4)]);
xlabel('t (sec)'); ylabel('y (cm)');

subplot(2,1,2)
plot(X(span),Y2(span),'o',T,damped_oscill(Pv,T))
axis tight;v=axis;axis([0 tmax v(3) v(4)]);
xlabel('t (sec)'); ylabel('v_y (cm/sec)');

%print -djpeg mio
