{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Script to illustrate using Sage to analyze data concerning the decomposition of butadiene.\n",
    "\n",
    "#Times at which data was taken (seconds)\n",
    "times = [0, 1000, 1800, 2800, 3600, 4400, 5200, 6200]; #Times in seconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Butadiene concentration, moles per liter at each time above\n",
    "data = [0.01, 0.00625, 0.00476, 0.0037, 0.00313, 0.0027, 0.00241, 0.00208];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the raw data versus time. Call the plot \"plot1\".\n",
    "pdata = list(zip(times,data))\n",
    "plt1 = scatter_plot(pdata)\n",
    "show(plt1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Does not look 0th order. Is it first order? Try a logarithmic transformation of the data (as was done for H2O2).\n",
    "log_of_data = [log(data[i]) for i in range(len(data))]\n",
    "pdatalog = list(zip(times,log_of_data))\n",
    "plt2 = scatter_plot(pdatalog)\n",
    "show(plt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Doesn't look first order (not a straight line). But fit a line y = -k*t+b to this data anyway.\n",
    "var('k, b, t') #Declare k,b,t as symbolic variables\n",
    "model(t) = k*t+b #Specify the model to be fit, with \"t\" as the independent variable\n",
    "sol = find_fit(pdatalog,model,parameters=[k,b]) #Fit the model by adjusting k and b\n",
    "f(t) = model(k=sol[0].rhs(),b=sol[1].rhs()) #Define the best fit function of the form in \"model\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now plot best-fit model overlayed on data\n",
    "plt3 = plot(f(t),t,[0,6200],color='red')\n",
    "pp = plt2+plt3\n",
    "pp.axes_labels(['time (seconds)','log(concentration)'])\n",
    "show(pp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Hmm, not too good. Doesn't appear to be first order..."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}