#Script to illustrate using Sage to analyze data concerning the decomposition of butadiene.

#Times at which data was taken (seconds)
times = [0, 1000, 1800, 2800, 3600, 4400, 5200, 6200]; #Times in seconds

#Butadiene concentration, moles per liter at each time above
data = [0.01, 0.00625, 0.00476, 0.0037, 0.00313, 0.0027, 0.00241, 0.00208];

#Plot the raw data versus time. Call the plot "plot1".
pdata = list(zip(times,data))
plt1 = scatter_plot(pdata)
show(plt1)

#Does not look 0th order. Is it first order? Try a logarithmic transformation of the data (as was done for H2O2).
log_of_data = [log(data[i]) for i in range(len(data))]
pdatalog = list(zip(times,log_of_data))
plt2 = scatter_plot(pdatalog)
show(plt2)

#Not first order either (not a straight line). But fit a line y = -k*t+b to this data anyway.
var('k, b, t') #Declare k,b,t as symbolic variables
model(t) = k*t+b #Specify the model to be fit, with "t" as the independent variable
sol = find_fit(pdatalog,model,parameters=[k,b]) #Fit the model by adjusting k and b
f(t) = model(k=sol[0].rhs(),b=sol[1].rhs()) #Define the best fit function of the form in "model"

#Now plot best-fit model overlayed on data
plt3 = plot(f(t),t,[0,6200],color='red')
pp = plt2+plt3
pp.axes_labels(['time (seconds)','log(concentration)'])
show(pp)

#Hmm, not too good. Doesn't appear to be first order...