#A worksheet to estimate spring and damping constants from experimental data.

#The Data: Load in the data, 1460 data points, the position (meters) of the
#mass every 1/50 of a second, starting at time t = 0.82 seconds.
load('spring_mass_data_clean.sage');
N = len(udat0);

#The data is in the array "udat0". A plot:
scatter_plot(udat0,facecolor='red',markersize=10).show();

#Center the data vertically and rescale so time starts at t = 0. First
#compute average of data (might be better to use integral number
#of cycles):
uave = add(udat0[i][1] for i in range(N))/N;

#Now subtract uave, make data start at t = 0:
udat = [[udat0[i][0]-0.82,udat0[i][1]-uave] for i in range(N)];
Tmax = udat[N-1][0]; #Time of last sample.

#A plot of the recentered data:
plt1 = scatter_plot(udat,facecolor='red',markersize=10);
show(plt1);

#Estimating the spring and damping constants: We will fit a
#function y(t) of the form
var('t d1 alpha omega');
y(t,d1,alpha,omega) = d1*exp(-alpha*t)*cos(omega*t);
#to the 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
plt2 = plot(y(t,0.05,0.014,7.0),(t,0,Tmax),color='red');
show(plt1+plt2);

#It may help to plot on a smaller time range, at first:
show(plt1+plt2,xmin=0,xmax=5,ymin=-0.05,ymax=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.