%% Spring-Mass Parameter Estimation
%% Kurt Bryan and SIMIODE
%
% A script to estimate spring and damping constants from experimental data.
%
% The Data: First load in plotting commands, and a command to read the data from an Excel spreadsheet, which should be in the current directory.

udat0 = readmatrix('spring_mass_data_clean.xls');

% The array "udat0" is an N x 2 matrix containing (time, position) pairs, N = 1460, with time starting at t = 0.82 seconds, position in meters

[N,~] = size(udat0)

% A plot of the position data over time.

plot(udat0(:,1),udat0(:,2))

% Recenter the data values so that the mean position is 0, also so that time starts at t = 0. First compute mean value of data over the total time interval (might be best to do this over an integer number of cycles):

uave = sum(udat0(:,2))/N;
udat = [udat0(:,1)-0.82 udat0(:,2)-uave];
Tmax = udat(N,1); %Maximum time

% Plot the recentered data

plot(udat(:,1),udat(:,2))

% Estimating the spring and damping constants: We will fit a function y(t) of the form

syms y(t);
syms t;
syms d1;
syms alpha;
syms omega;
y(t,d1,alpha,omega) = d1*exp(-alpha*t).*cos(omega*t)

% to this data.
% Based on the plot above we can guess that the initial amplitude d1 is about 0.05. We can estimate alpha by the rate of decay of the amplitude of the oscillations. For example, at time t = 25 the amplitude is down to about 0.035, so 0.05*exp(-alpha*25) = 0.035, which leads to alpha equal to about 0.014 (solve 0.05*exp(-alpha*25) = 0.035 for alpha).
% We can estimate omega by estimating the period of the motion, e.g., count how many complete oscillations the mass undergoes during the approximate 29 second data set. We can then plot y(t) with these estimated values, as

plot(udat(:,1),udat(:,2))
hold on
fplot(y(t,0.05,0.014,7.0),[0 Tmax],'-r')
hold off

% It may help to plot on a smaller time range, at first:

plot(udat(:,1),udat(:,2))
hold on
fplot(y(t,0.05,0.014,7.0),[0 Tmax],'-r')
hold off
axis([0 5 -0.05 0.05])

% Obviously some adjustment in omega and perhaps the other parameters is in order.
% Adjust d1, alpha, and omega to obtain the best (visual) fit possible, then use the formulas in Modeling Exercise 6.3.5 to estimate the spring and damping constant.
%