% Script fits a damped sine function to experimental oscillation data. 
% y= amplitude * exp^(nt)* sin(omega*t + phi*t)
clear all
close all

%% Load data
testnr=1;
%fname=['./data/filename' num2str(testnr) '.txt']; % loads a file called filename_1.txt i data subfolder
%M=load(fname); % M=matrix with testdata if no header in file
%t=M(200:end,1); Make sure to cut data at index with nice data, plot(t,v) first to check where
%v=M(200:end,2); 
% Dummy experimental data to run script first time t, v - typically loaded from file, see above
t=0:0.005:10; % [s] time
v= 5 * exp(-0.5*t) .* sin(3.23*2*pi*t + (-pi)) + 0.1.*sin(80*2*pi*t + (-pi/8)) + 5; % Exp y data
figure(1)
plot(t,v)

%% Guess for frequency and damping for damped sine function
x(1) = 3 *2*pi;  % x(1)= omega=f*2*pi, f is the frequency guess
x(2) = 0.0;     % x(2)= n, damping guess
x(3) = -pi;     % x(3)= phi, phase guess
x(4) = max(abs(v));  % x(4)=amplitude guess
x(5) = 0; % Vertical displacement, if not oscilating about zero
x = lsqcurvefit(@expsin, x, t, v); % Curvefit

%% Eigenvalue frequency and damping and fitted y-data
c = x(2) + 1i*x(1);
fprintf(1, '-------------------\nRESULT: \nFrequency: %.3f Hz \nDamping: %.3f\n-------------------\n', x(1)/(2*pi), x(2));

y = x(4) * exp(x(2)*t) .* sin(x(1)*t + x(3));

figure(2)
plot(t,v,'b',t,y,'r--')
xlabel('Time (s)')
ylabel('Amplitude')
legend('Experiment','Curve fit')



%% Function call to self-defined function: Sometimes works to have here at end of same script, 
% uncomment and test, or have this expsin function as separate file in the same folder
% function y = expsin(xd, td)
% 
% y = xd(4) * exp(xd(2)*td) .* sin(xd(1)*td + xd(3)) + xd(5);
% 
% end