FittingExercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting of experimental chemical kinetics data\n",
    "You perform some experiments in a batch reactor to determine the rate expression and thermochemistry for the reversible chemical reaction\n",
    "$\\require{mhchem}$\n",
    "$$\\ce{A <=> B}$$\n",
    "\n",
    "Recall from thermodynamics that \n",
    "$\\Delta G = \\Delta H - T \\Delta S$\n",
    "and $\\Delta G = R T \\ln K_a$\n",
    "where $K_a$ is the activity-based equilibrium constant of the chemical reaction, $R$ is the molar gas constant (8.314 J/mol/K) and $T$ is the temperature in Kelvin.\n",
    "\n",
    "If we assume ideal solution (unit fugacity coefficients) then $K_a = K_c$ giving us the concentration-based equilibrium constant $K_c$. \n",
    "From kinetics you recall\n",
    "$K_c = \\frac{k_f}{k_r}$ \n",
    "where \n",
    "$k_f$ is the forward rate coefficient and $k_r$ is the reverse rate coefficient.\n",
    "i.e. the rate of the reaction $\\ce{A->B}$ is $k_f \\times C_A$\n",
    "and the reverse reaction $\\ce{B->A}$ is $k_r \\times C_B$\n",
    "where $C_A$ and $C_B$ are the concentrations of species A and B respectively. \n",
    "In a batch reactor $\\frac{dN_A}{dt} = r_{A(net)} V$, so (dividing through by the reactor volume $V$) $\\frac{dC_A}{dt} = r_{A(net)}$ where $r_{A(net)}$ is the net rate of formation of species A, i.e. $r_{A(net)} = k_r C_B - k_f C_A$.\n",
    "Assume the forward rate coefficient $k_f$ follows Arrhenius form, $k_f = A \\exp\\left(\\frac{-E_A}{R T}\\right)$ where $A$ is the \"pre-exponential factor\" and $E_A$ is the activation energy.\n",
    "\n",
    "Fortunately, in this case you have good reason to believe that species A and B have very similar temperature-dependent heat capacities, so that $\\Delta H_{rxn}$ and $\\Delta S_{rxn}$ are independent of temperature.\n",
    "\n",
    "You start the experiment with no B ($C_B=0$), and at time zero have some way to initiate the reaction, starting with a set concentration of $C_A$. \n",
    "\n",
    "You wish to determine the four paramaters:\n",
    "$log_{10} A$,\n",
    "$E_A$,\n",
    "$\\Delta H_{rxn}$,\n",
    "$\\Delta S_{rxn}$.\n",
    "\n",
    "Based on a literature search, quantum chemistry calculations, and prior experience, your current estimates are as follows:\n",
    "```\n",
    "logA = 6.  # base-ten logarithm of A in s^-1\n",
    "Ea = 45.   # Ea in kJ/mol\n",
    "dH = -10.  # ∆H in kJ/mol\n",
    "dS = -50.  # ∆S in J/mol/K\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy.integrate\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import random\n",
    "import SALib as sa\n",
    "import SALib.sample\n",
    "# from SALib.sample import morris as ms\n",
    "# from SALib.analyze import morris as ma\n",
    "# from SALib.plotting import morris as mp\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This cell just tries to make graphs look nicer\n",
    "try:\n",
    "    import seaborn as sns\n",
    "except ImportError:\n",
    "    # This block will be run if there's an ImportError, i.e you don't have seaborn installed.\n",
    "    sns = False\n",
    "    print (\"If you want to try different figure formatting, \"\n",
    "           \"type 'conda install seaborn' at an anaconda command prompt or terminal. \"\n",
    "           \"See https://stanford.edu/~mwaskom/software/seaborn/ for details\")\n",
    "    # If not using seaborn, we can still control the size of the figures this way\n",
    "    from pylab import rcParams\n",
    "    rcParams['figure.figsize'] = 3, 3\n",
    "else:\n",
    "    # This block will be run if there is no ImportError\n",
    "    sns.set_style(\"ticks\")\n",
    "    sns.set_context(\"paper\",rc={\"figure.figsize\": (2, 2)})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a \"named tuple\" data type to store the exprimental data in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "ExperimentData = namedtuple('ExperimentData', ['T', 'cA_start', 'times', 'cA'])\n",
    "\n",
    "def plot_experiment(e):\n",
    "    \"\"\"\n",
    "    Plots the experimental data provided in 'e' \n",
    "    which should be of the type ExperimentData.\n",
    "    \"\"\"\n",
    "    plt.plot(0, e.cA_start, 'ko')\n",
    "    plt.plot(e.times, e.cA,':o', label=\"T={:.0f}K\".format(e.T))\n",
    "    plt.ylim(0,)\n",
    "    plt.ylabel('$C_A$ (mol/L)')\n",
    "    plt.xlabel('time (s)')\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now here are the data from your three experiments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment 0 was at T=298.15K and ran for 100 seconds\n",
      "Experiment 1 was at T=308.15K and ran for 100 seconds\n",
      "Experiment 2 was at T=323.15K and ran for 100 seconds\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKsAAACgCAYAAABzL5zaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHfJJREFUeJztnXtcVOW6x7/DDLdhYAYBQVBBMa+omFnm1nCreUuy1NKT\nWzTNWxfPTtuV5SXbptt27crKraIbyywvpRaH8piZHnWLJWgXcKeJgFxDbgMkAjNz/hhnCTIDM8wM\n48L1/Xz81Fq886x3zfrNu971PM/7LJnBYDAgISEC3FzdAQkJa5HEKiEaJLFKiAbRinXXrl2SXSfb\nvtXsSmIVuV1n2r7V7IpWrBK3H60q1qysLB5++GEAtm7dyquvvsrzzz9PaWlpa3ZDQqS0mlivXLnC\np59+ilKppKamhu+//54VK1YwefJkdu/e3VrdkBAxitY6UGBgIEuWLGHu3LmUl5cTGBgIQEhICEVF\nRTbZqq6upqqqiqysLORyuUP7WVNTQ05OjkNtOtOuM207w65Op6Oqqorq6mq8vLxs+myribU+AQEB\nlJWVAVBQUED79u0ttt21a1ejCXlVVRWZmZmMHj3aKf0bOXKkqOw607az7E6cOBEfH58G+6ZOncrU\nqVMtfkbW2uHWuXPnEh8fz/bt28nIyECr1bJq1SpUKpXVNrKyshg9ejTZ2dnU1dXx8ccfc++99zqx\n1xKOoqCggOnTp3Pw4EHCw8Nt+myrj6zx8fEAzJgxo8U2TLf+uro6+vTpwyOPPOKQvkm0Hi2Zvona\ndfXxxx9z9uxZV3dDopUQtVilW//thajFKnF7IYlVQjRIYhUJpdpqUv5TSKm22tVdcRku8bNK2MYX\nxy6SkJhGnc6AQi7j8dg+PDgs0iYbGzZs4NKlSyQnJ3PvvfcSGRnJ/PnzG7RZvXo1CoWCwsJClixZ\nQnFxMQkJCQQFBTFgwADGjx/Ptm3b+O233ygsLGT+/PmkpaXh6enJ+PHj2b17N+fOnWPZsmUOD9aA\nJNZbiuLyqwAEqL2F7bKKa4JQAep0BhIS0xjWPwz9dRe5qX1TPPnkk4DRz/3666+zc+dOli1bJvz9\nmWeeYdiwYcTExHDw4EFOnDjBxYsXWbhwIT169GDx4sWMGDGCr7/+muDgYAoKCoQopMFg4O2330aj\n0bBy5UrHfSE3IYn1FuJ/k7MAeGxMT2E770qVIFQTdToDGXnl/JJV2qC9Leh0OnQ6nbCtVCqJiYkh\nOzubr776itdee43S0lLef/991Go1v//+OxUVFSxcuJChQ4dy4sQJPvnkE0JDQ9m8eTP+/v4sWrSo\npaduFZJYbyHGDA5vtF1WeY0TP+Q2EKxCLqNrqJqIDn42H8MUsJw+fXqjvx06dIjk5GT+9re/4enp\nSXp6Ok8++SQdO3ZkwYIFyOVyPv74Y4YOHUpgYCBVVVXIZDJmzpzJ2LFjmT9/PmvWrKFTp04298sa\n2vwDlpgeTALU3g1u6QFqbyLDNDwe2weFXAYYhTo7Ngp/P69G7a1BJpOZ3Z+dnc3y5cvRarUsX76c\nw4cPExwczLp163j55ZcZOXIk7dq1IzQ0lBUrVrB161bi4uIwGAx4eXmhVCpZvXo1zz//PJWVlS3/\nEprquxiXYufk5DBy5Ei++eYbOnbsaLGdIx5MbhVKtdVk5JXTNVSNv59t2Uq3EtZeO3O02ZG1RFtt\n9sFEDCOsOfz9vBjYM1jUQrWXNivWS3nlFh9MJMRJmxVr11C1MM8zYXowkRAnbVas/n5ejR5MxgyO\nIL+4ysU9axllV8s5m59G2dXb987Qpl1XDw6LZFj/MDLyyunSwY/N+39GpxPd8yRfnj/M9h/2otPr\nkLvJmdF/EuO7j3B1t1qdNi1WuP5gcv2h5MWZg1zcm6YpuWpc6tPOWyNsl1drBaEC6PQ6tv+wlyGd\nBqLH0KB9U1gTbn3nnXeoqqqisrKSl156iaqqKtatW4dGoyEyMpLp06ffnuHWwsJC3nvvPWE5ywsv\nvNBqx9brDRxJzWH4nR1xczPvd3QFhy4eB+DRqAnCdkHFb4JQTej0OjLLcjhffKlB+6ZoLtz6wgsv\n0KtXL0aPHs3WrVtJS0vj1KlTxMXFER0dzfz585k6dertGW7NyMggOTmZqKgoeva0PVxoD79klXLo\nu2wGR4Wg9HJv1WM3xajIoY22y6srOJmT2kCwcjc5EZqOdNaEtfhYN4dbAUaPHs3OnTtJSkriwQcf\nJCkpiQ4dOgDg6+tLZWXl7RluDQkJISEhgY4dOzJnzhxmzZqFp6dno3bmVrfW1NTYdexeXdrx2sIh\nFqM5ruLm23k7bw3tvDXM6D+pwZw1rv9kNN4t82pYCrfW1dVx+vRppk2bxoABA4iPjyc0NJT8/HyC\ng4PRarXo9XqHhVsXLFiAh4dHg33NrW51mVh37NjB5MmTAVCpVNTV1ZkVq7kTMEVB7MEk1BJtNWUV\n1+gaduu6tMZ3H8GQTgPJLMshQtOxxUIFy+FWhULBV199xVdffUVVVRUzZ84kKCiItWvXsm/fPu6/\n//4G4dbq6moWL17MyZMnG4Vb4+Pjm12tvHHjRpsjWC4Lt6anp7N+/XrCwsIICwtj9uzZVn/WnpDd\nzazZ9h19IwOJHdbVLjsS1mHPtXPZyNq7d282btzoqsMLPD/jLhRyo7u5rcTf2ypt3nXVHCahfnHs\nIv/6Ig2dXvxJL22VNhvBsoUSbbUgVBB/0ktbpc2L1Zow5aW8ckGoJup0BvYe+dXZ3bOamtJSSlPP\nUHMblwdt09MAa8OUpqSX+llabjIZnYJ9he0TP+YR7K+kWyeje6k157d5iUlkbvsQQ10dMoWCiFlx\nhMY+YJONlkSwLly4wM6dO1GpVLRr146nnnqqQZulS5dy6NChth/BcjalV8sthilvdv2Ykl7qJ2rP\njo1i9D03lpkk/5wvbH9x7CJbv0hD7+D57bXiEgA8A9oJ2zVlZYJQAQx1dWRu+5DAoUMwXL8bmNo3\nRUsiWNXV1axcuRKlUsmcOXMAGrRJT0839qmtR7CcTVZZjsUwZbQZP2X9pBdzo+WSxwYCN5K69frG\nq03tHWELD34NQOf/mipsX83PF4RqwlBXR9WlTCp+Od+gvS00F8GaOHGiEE7duHEjsbGxZtvk5eW1\n/QiWs4nQdETuJjcbprRE/aQXS1hK6j6XWcKQfqF29Tl49P2NtmvLyyk+cbKBYGUKBT5dIlDaWDIS\nbItgLVq0iDVr1hAbG8vgwYOpq6vj7NmzQpvNmzfTu3dvacGgvWi81czoPwm5m3HuJHeTc09YNP++\nnGKXXXNJ3TIZpF8qscsuGG/n9W/pngHtUHXtQsSsOGQK47giUyiIeHwmHv7+jdpbQ3MRrL/+9a8k\nJCTw0EMP8dprr5GVlcXevXtZunRpozam90NICwabwJYoSNnVciFM+e/LKdzZIYoQX8uVtq3h5oWI\nk/7YjamjeuDh7viHChM1paVUXcrEp0sEHv7+TjuOsxFlBKu10HirhTmqoxKWm5rfXswpI/nnAqaP\ndWwmmYe/v6hF6gisFmtKSgo//PAD5eXl+Pv7c+edd9KvXz9n9s1pVNdWszV1F49ETaC9T0CLbFia\n317K09I5xNfMJyTspVmxHjx4kMTERCIjI+natSs9evSgoqKCb775hi1btvDAAw8wZsyY1uirw6jV\n1xnT77wcn2k16u7Owv8bDAayCyvwU3pIOQcOoFmxqlQq3n333Ub7x44dC8CxY8cc3ysn4+up4r/6\nTXT6cf79Yz5bPv+JssprbaLQhqtp1hswZMgQs/v37NkDwLBhwxzbo1YmoySbbWf2OMV2984aSiuu\nNSi08S8rcg7EVPKoNWmx68qZL85tTQqriugZaBzpHL3cObuwolHOga5eoY3Xt5/mStlV4W+5RZXs\nP/ors1cf5JX4ZGavPsgXxy46pC9tgTbvDWiOezsZI1POWO5sLuegfqENpZcCldK4BsxgMPDsW0eo\nqdMLy8WtiY7dTjm4zYp1x44dZveXlNjnBM/NzWXDhg2oVCo0Gg0LFy60y5492JJHYAuWcg5Monr6\nkWihrcEAj43uydbEtAY2TCWP+vt4sP3Lc8ya0Ftw7LelwnPW0KxYlUql2f3PPPOMXQdOSEigc+fO\nZGZmOvU1kdZgax6BLTSXc2DCzU1GzJ0d+eDLdLMj8ZWyq2TklQtCzcgtZ8vnP2MK6dwOo3CzYlWr\n1dx99902va7SGrKyspgyZQpz5sxh9uzZfPjhh2bbOWN1682YzSOQydFec0zY0JqcA1O7pkbiv86/\n8bB7LrOYm2OPplF4oJ8Xv5X8zi9ZpQwbYFyufauNwk5Z3arX61m/fj0VFRW0b9+ewYMHM3DgwEYH\nspWgoCB8fHxQKBRN/hCctbq1PqY8gvpz1qj2PSh1QV0pa0fiIX1D2fL5zw1zcN1uzIcz87Ukp+Uz\nbEAYJdpqIaURbnglBkd1oL2/+TunsyivvAa0bHVrs2IdNWoUo0aNAqCoqIjk5GTWrl1LdXU1a9eu\nbUF3jcyZM4c333wTX19fxo8f32I7juLm5c5qLz/qp02czU+jb3BPITHGmVgzEjc3Ct/dJ4RBvYMB\nY6aY3oxX4v09P7BqnvEtjeWV1/Bwl+PtaZSEM6YMXxy7yKZdLffLNyvWo0ePNtj28/Nj+PDh2Jv/\nEhkZydtvv22XDUdTP48AbmQo5VUUEp/yCW+OWSaItX6CjD0PYfbQ3Chs6r85r4TcTcbkEd2E7c++\n/RWll4Jp9/ewefGkNcIuKK4iITHNrsJ4zYr1p59+AozTga+//prs7GzCwsIYN24cw4cPb/GBxUSo\nbzBvjFmGl7vxQuz++X/Yl34AncH1Vf3sGYX7dQsS2oy9NxwPhVxILr958WTl77UM6N6eXl2MKYkl\n2mpU3u4cSM40Oxf+vzM53NUrWCjPtPjto43ygG2lWbE+/fTTACxdupSYmBj69u3Ljz/+SFFRkV0H\nFhve14VaerWcz9KSMH3tOr2O7WebdnOJYRQODTQ+N6T8p9BscvnlwgoG9ryRWvnOzjPE3Blm8R1d\n+49eJCTAh+6djZlicydG8c7us9TacQ5WBwVycnKEOeqYMWN47LHH7DiseMkqy+Hm8UFn0JFedIEh\nne8CYP+5/2XcHX/EU+FxS9VWtWYUthTImPdQ3wYCnzGuFwUllt/R9fozw4SaDAB/vKszFVdr7Zqz\nWh1ulcvlHD9+nPLyco4dO4aXlzj8dI5ewmxyc9XHTeZGr0Dj/K9Or+Pri8dwd1MYgw1nGwcbTl0+\ni96gt3gMV1a5NlcxvP6Dm4lunTT06RJgsRR+faGaeHBYJK8/3fJcEqtH1r///e9s2rSJDz74gPDw\ncLs8Aa2FI5Yw34w5N1dc/8n4K41LtBVuct6fsBq4HmwwNA42/PP7D7i74z8AuFZXw9cX/48JPYwe\nF1tHYmdMMax1nwlz4a9S0XuV41atZva4gU16D2SKlvvIrRZrUFBQg6W7tzo1JaUWlzDbm3FvbVU/\nS8GGZTH/LTypF1YWcTr3Ryb0GHVjJDbUG4nP7qWrf2d6BnVrZN+Zwpa5X0OuvoLM3QuwLD5FcBZe\n0UeND5syOfLgAMB8kbsvzx9myxHz4XtrsFqs27dvJykpCTAmXchkMnbu3NniAzubqsxMi0uYHbE8\n5GY3l6U25kbhbgERQpvOmjBeGbEYsDASG3RsTdnF38e+DBhTGpNzUhl3xx/N5jMMCuuPr4eP4Lkw\nYYuwk375ho9+3GdsK5Mzpc8DTO4zDjCWjtdWVxLh3/H6j+szdNenNDqDju1nP0PhJmd0txgALhRf\nIqssh4Gh/dj+w170N4W1bcFqse7bt4+dO3faHblqLXy6RCBTKBotYa4pLqY09Qz+dw5olX7YUlvV\n0kj8xF3ThO06fR1Kd2+L+QyHLh4nV1vAc0ON1VZS837meNb3nMxJMZuoc6EkkzP5acy7y/jAfPDX\no3x49jMM1x8jdQYdu39OZGTXIWi81ZwrusB3OT/w7JAnrv+4Gs69dQY9Z/LSBLFera3myu8lZvtr\nK1aLdejQoaSkpBBeb616aKh96+SdiYe/PxGz4hrOWR+fScWFi7R30rp2S1gzCpvamRuJewTecMh3\nD+xK98CulF0tN1sX4e6O0Q32BSg1BCr9LSbqhKiCiGrfXdjv7uYuCNWEAYOQ1DMwtB99g3sBlmsz\nzB90oyZBv5Be9AvpJfTXq9ryg2VzWC3W3NxctmzZgo+PD2CMjrzzzjstPnBrEBr7AIFDh1hcwqyv\nqSHzg+1EPD4TN4XCpuXOzloabe1IbEnYke0aFr4I13RE7elL4vlDZgt+aLzVdFLfGHQGdOiD3E2O\nZ1UNQaV1FPkruObjIRQH8VJ44qXwbLIP5vqs8VYzpyIS7deHOdPC78ZqsZaVlbF169YWHsZ1NLWE\n+fecHPS1tbgpFDZ5Dmz1MtgqbGtH4vHdRzDItxtZaamE97mToA6dzbYzieqz5D20K75GSYAnUwY3\nLSqPxOPI9aBzg9rYYRZ/NNb2oaakFGVSMlUtH1itF6uvry+vvvoqkZE3bknm3lkvJlRdu9LtyQUW\nPQcKXxXqPr3xDDKGJcvT0pErvW3yMrTEfWatuOvbvqD4hNombEf/chXN/mIMOh0yuZwI/9/hxt0f\ng8HA1dxc5N5KlEnJmKaicj0okk5S8/BM3DWaRhVdLPVBX1vLtaIivK9PFYtPfdfogddWrA4KmEKt\nSqVS+NdWsOQ5yPs8kd8v5wj7cnZ/Sun3KRa9DOffWk95WrqwP/fzL8hM+KChsBM+aDJAkZeYxOkn\nFpC+ajWnn1hAXmKS8Lfyn9PQX7dVU1LKpS3/amD70pZ/UZWVLbT/5Y1/UFdZeePHeL0Qm0GnI3Pb\nh6QuWozumjFlTyaT8cPi56k8f97s+VVmXOLUY3FCe4CLG+PJ3PZBox9uTWkptWVl/PTScqGtul8U\nuNlXrarZT3/77bcAPPzww43+ARw+fNiuDtwKmDwH9ZEpFPRe8TKaATeWnvR44TmChseYbevTJQLf\nnj3wDLxRNOPKv5MFgZgw6HRUXcoE4D/r3qDiwo2Cxdm79jQp7l/f3UB1QSFg/IGZo7qg4Ea/3OTI\n5HKLP0afTmEN+tfn1ZX4dOtm/vzCwwmPm478+ht19DU1lKakYqjTNbJbdSkTj8BAei19QcjOU4aF\n0WX2LGSKlqdYNitWd3d3Fi1axLp16/j00085dOgQ+/fvZ+3atTzzzDNOKRrb2pg8B+aKn9W/7SmU\nSrzaB1ls22HcGLyCg4X2vV78S+MLL5fj0yUCAL/evfBod+M2X5qS2qS4o1avwivEaN/SD8y3+x3C\ndvfF/43c29ti2y5PzEZR7w7p17MHXkGBZs/PMzCADuPGCm3dPDzo9/oaiz9cmUyGb4/uDb6/0NgH\n6Pu312gpVhdmO3fuHGfPnkWr1aJWq4mOjm71NwOacOSrherjDG9Aoznr4zMJnWA+2bymtJTTTyxo\n5Bu+a8tG6+bDTdi2pa2zzg/svHYGF7NkyRJDUlKSTZ+5fPmyoXv37obLly87qVeO5VpJiaEkJdVw\nraSk2ba5X/yP4cSkRw3HH5xkODHpUUNuYtPfjS22bWlrC7bYtefaubRuwLZt2wS/bVvGlgqAzfmG\n7bHtrEqErVXh0C6xFhUVERQU1HxDMxw+fBhfX1+io6Obb3ybIZW3NI9NYv3nP/9JWVkZBoOBiooK\nysvL2bBhQ4sOnJiYiFqtJiMjA3d3d/7whz+gVjd2PLfGUmyJ1qclS7FtqnxdWVlJamoqSqUSd3d3\nunXrZvdtfP/+/Xh4eNi0wtVZD1gSzqfVKl+rVCruu+8+ysrK+Pbbb8nNzbV7GfVDDz1k1+clbh+s\nCikkJCSQlJQkvNhAo9Hw8MMPo9FonNo5CYn6WCVWf39/unTpgpubG1qtlpMnTwKWa7dKSDgDq8Sq\nVCrp3bs3SqUSPz8/+vbtS0qKfa/okZCwFavEWnpT4oVKpXLau44kJCxhlVh9fHw4cOBAg335+flO\n6ZCEhCWs8gZMmDCBt956iz179jBkyBDJzynhEqx2XT377LOkpaVx8uRJwsPDb4nKfxK3Fzb5Wfv0\n6UOfPn2c1RcJiSZpsy8almh7SGKVEA2SWCVEgyRWCdEgiVVCNEhilRANklglRIMkVgnR4LIFg2fO\nnGHnzp2oVCratWvHU0895aquSIgEl4lVq9WycuVKlEolc+bMcVU3JESEy8QaE2MsNrtx40ZiY2Nd\n1Q0JEeEysVZVVbFmzRpiY2MZPHiwxXbS6ta2idNXtzqSl156iezsbEJDQ5HL5Ta9/UVa3SpeWm11\nqyNZs2aNqw4tIVIk15WEaJDEKiEaJLFKiAZJrBKiQRKrhGiQxCohGiSxSogGSawSokESq4RokMQq\nIRoksUqIBkmsEqJBEquEaJDEKiEaJLFKiAZJrBKiQRKrhGhw2UqBwsJC1q1bh0ajITIykunTp7uq\nKxIiwWVi3bVrF3FxcURHRzNv3jymTZuGXC636rM6nQ6AgoICZ3ZRwgmYrpnpGtqCy8R65coVOnTo\nAICfnx8VFRVmXwJnbnVrVVUVgDQai5h58+Y1epVqc6tbXSbW0NBQCgoKCA4ORqvV4ufnZ7aduROo\nrq5m4sSJbN682erR2FoWLFjAxo0bHWrTmXadadsZdnU6HfPmzePzzz/Hy8vLps+6TKxTpkxh7dq1\nqFQq7r//ftzcrH/W8/LywsfHh/DwcIf3y8PDwynLu51l15m2nWXXx8fHZqGCC8UaGBjIm2++6arD\nS4gQyXUlIRoksUqIBvkrr7zyiqs70VKioqIku062fSvZdVmtKwkJW5GmARKiQRKrhGiQxCohGiSx\nSogGlwUFWoqjs7VufhGHl5cXubm5VFZWsnTpUvz9/e2y/9xzzzFixAjy8/MdZjc3N5cNGzagUqlQ\nq9V4eno6xHZhYSHvvfceKpUKMAZu7LWblZXFn//8Z/bt28fWrVsb2KupqbHpWopuZDVla61YsYKj\nR4+2KHunPqYXcSxfvpzU1FROnz7NihUrmDx5Mrt377bL9rZt24RkDUfaTUhIoHPnzmi1WqKiohxm\nOyMjg+TkZAoKClCr1XbbvXLlCp9++ilKpZKamhq+//57VqxYwZQpU9i9e7fN11J0YjWXrWUPMTEx\nKJVK4UUcAQEBAISEhFBUVNRiu4cPH8bX15fo6Gj0er3D7IJxtIqJiWH16tVs2bLFYbZDQkJISEjg\nrbfe4rvvvhNG0pbaDQwMZMmSJSiVSsrLywkMDAQgODiY3377zeZrKbppgLXZWtZS/0UcgwYN4tCh\nQ4Ax77J9+/YttpuYmIharSYjIwNAGGHttQsQFBSEj48PcrkcLy8vysrKHGJ7x44dTJ48GQCVSkVe\nXp7D+hwQENCgn8HBwej1evLz862+lqILCly5ckXI1oqKiuKRRx6xy97NL+Lo3bs3GRkZaLVaVq1a\nJczfWsr+/fvx8PCguLjYYXYvXrzIu+++i6+vL/fccw+lpaUOsZ2ens769esJCwsjLCwMd3d3h9id\nO3cu8fHxbN++vYG96upqm66l6MQqcfsiujmrxO2LJFYJ0SCJVUI0SGKVEA2SWCVEgyRWJ1JWVsaR\nI0cAWLlypd32Lly4wCeffNJof21tLatXr7bb/q2OJFYncv78eU6dOgXAqlWr7LYXHx/PpEmTGu13\nd3cnLCyM1NRUu49xKyOJ1Yl89NFHHDhwgJSUFCEyFBcXx1/+8hcmT57Mjh07mD9/PpMmTUKr1VJc\nXMy8efOYMWNGoxcxl5aWUltbi6enJ6dOnWLatGk8+uijJCUlATBixAg+++yzVj/H1kQSqxP505/+\nxLhx4xg4cCAymQwwhhpffPFFli1bxpdffsmmTZsYPnw4qampbNq0iVmzZrF9+3YMBkODkfKnn34i\nIiICgCNHjrBgwQI++ugjod5CeHg46enprX6OrYkk1lbG19eXgIAAVCoVnTt3Boxx+GvXrnHp0iXe\nf/994uLiOHPmDPn5+cLntFqtUF5p7ty5HDlyhLlz53L16lWhjaOr09xqiC6RRUzIZDJM0Wxrotrh\n4eE8+OCD9OvXjwMHDnDHHXcIf1Or1WRlZQGQlJTE/Pnz6dChAw899JAwj1Uo2vblbNtn52I6derE\n8ePHue+++4RpgOm/NyOTyZg3bx4vv/wylZWVBAUF8cYbbwh/79evH3v27AGgZ8+eLFy4EF9fXyZM\nmAAY0wZ79uzp5DNyLVIii4h47rnnWL16tdk6UVu2bGHQoEH079/fBT1rHaQ5q4iYO3cue/fubbS/\ntraW3NzcNi1UkEZWCREhjawSokESq4RokMQqIRoksUqIBkmsEqJBEquEaPh/wmqRXnUgalQAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11639cb00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from numpy import array\n",
    "experiments = [ExperimentData(T=298.15, \n",
    "                              cA_start=10.0, \n",
    "                              times=array([ 10,  20,  30,  40,  50,  60,  70,  80,  90, 100]), \n",
    "                              cA=array([ 8.649,  7.441,  7.141,  6.366,  6.215, 5.990, 5.852,  5.615,  5.481 ,  5.644])), \n",
    "               ExperimentData(T=308.15, \n",
    "                              cA_start=10.0, \n",
    "                              times=array([ 10,  20,  30,  40,  50,  60,  70,  80,  90, 100]), \n",
    "                              cA=array([ 7.230,  6.073,  5.452,  5.317,  5.121, 4.998,  4.951,  4.978,  5.015,  5.036])), \n",
    "               ExperimentData(T=323.15, \n",
    "                              cA_start=10.0, \n",
    "                              times=array([ 10,  20,  30,  40,  50,  60,  70,  80,  90, 100]), \n",
    "                              cA=array([ 5.137,  4.568,  4.548,  4.461,  4.382, 4.525,  4.483,  4.565,  4.459,  4.635])),\n",
    "              ]\n",
    "\n",
    "for i,e in enumerate(experiments):\n",
    "    print(\"Experiment {} was at T={}K and ran for {} seconds\".format(i, e.T, e.times[-1]))\n",
    "    plot_experiment(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ParameterSet = namedtuple('ParameterSet', ['logA', 'Ea', 'dH', 'dS'])\n",
    "\n",
    "# This is a sensible starting guess for your fitting\n",
    "starting_guess = ParameterSet(\n",
    "    logA = 6. , # base-ten logarithm of A in s^-1\n",
    "    Ea = 45. , # Ea in kJ/mol\n",
    "    dH = -10. , # ∆H in kJ/mol\n",
    "    dS = -50.  # ∆S in J/mol/K\n",
    "    )\n",
    "\n",
    "# This should end up with your optimized parameters\n",
    "optimized_parameters = ParameterSet(0,0,0,0)\n",
    "\n",
    "\n",
    "# This should end up with your standard errors (one sigma)\n",
    "# for the uncertainties in the fitted parameters.\n",
    "# i.e. there should be a 68% chance the true value is \n",
    "# at least this close to your optimized parameter.\n",
    "standard_errors = ParameterSet(0,0,0,0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, now insert some cells to determine the optimized_parameters and their standard_errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def func(x, logA, Ea, dH, dS):\n",
    "    return a * np.exp(-b * x) + c #need to come up with what form this should be...\n",
    "\n",
    "coefficients, cov = curve_fit(func, xdata, ydata)\n",
    "\n",
    "# fit for each given temp \n",
    "\n",
    "# for T=298.15K\n",
    "x = [0, 10,  20,  30,  40,  50,  60,  70,  80,  90, 100] #time\n",
    "y= [10, 8.649,  7.441,  7.141,  6.366,  6.215, 5.990, 5.852,  5.615,  5.481 ,  5.644] #Ca\n",
    "\n",
    "\n",
    "def func(x, a, b, c):\n",
    "  return a*x**b + c\n",
    "\n",
    "# determining the model form with using the designated parameters\n",
    "\n",
    "# dH-T*dS= R*T*ln(Kc)\n",
    "# Kc=Kf/Kr\n",
    "# Kf=A*exp(-Ea/(R*T))  --> log(Kf/A)=-Ea/(R*T) ... log(Kf)=logA*(-E/(R*T))\n",
    "\n",
    "# dNa/dt= ra_net*Vol --> integrate --> Na = ra_net*Vol*t + Na_0\n",
    "# ra_net=Kr*Cb-Kf*Ca\n",
    "\n",
    "# Na = (Kr*Cb - Kf*Ca)*Vol*t + Na_0\n",
    "# Na= (Kr*Cb - A*exp(-Ea/(R*T)))*vol*t +Na_0\n",
    "\n",
    "# dCa/dt = ra_net = Kr*Cb - Kf*Ca\n",
    "# integrate --> Ca = Kr*Cb*t - Kf*Ca*t + Ca_0 --> need this in terms of logA, Ea, dH, dS (and time, t)\n",
    "# dH-T*dS = R*T*ln(Kf/Kr)\n",
    "# exp((dH-T*dS)/(R*T)) = Kf/Kr --> Kr = Kf*exp(dH-T*dS)/(R*T)\n",
    "# Ca = Kf*exp(dH-T*dS)/(R*T)*Cb*t - Kf*Ca*t + Ca_0 --> take log\n",
    "# log(Ca) = log(Kf)*(dH-T*dS)/(R*T)*Cb*t - log(Kf)*log(Ca)*t + Ca_0 --> factor out log(Ca)\n",
    "# log(Ca)*(1+log(Kf)*t) = log(Kf)*(dH-T*dS)/(R*T)*Cb*t + Ca_0\n",
    "# log(Ca) = (log(Kf)*(dH-T*dS)/(R*T)*Cb*t + Ca_0) / (1+log(Kf)*t)\n",
    "# log(Ca) = (logA*(-E/(R*T))*(dH-T*dS)/(R*T)*Cb*t + Ca_0) / (1+logA*(-E/(R*T))*t)\n",
    "\n",
    "# still have Cb ... shouldn't be there\n",
    "# Cb = Ca_0-Ca ??? I think this makes sense...\n",
    "# log(Ca) = log(Kf)*(dH-T*dS)/(R*T)*(Ca_0-Ca)*t - log(Kf)*log(Ca)*t + Ca_0\n",
    "\n",
    "# new problem : have both log(Ca) and Ca... can't factor that ...\n",
    "\n",
    "# also the graph is just Ca, not log(Ca) ... getting rid of that results in A not logA which is a parameter...\n",
    "\n",
    "\n",
    "\n",
    "#yaj = func(x, coefficients[0], coefficients[1], coefficients[2])\n",
    "\n",
    "print(coeffs)\n",
    "print(cov)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ParameterSet(logA=6.0, Ea=45.0, dH=-10.0, dS=-50.0)\n",
      "ParameterSet(logA=0, Ea=0, dH=0, dS=0)\n",
      "ParameterSet(logA=0, Ea=0, dH=0, dS=0)\n"
     ]
    }
   ],
   "source": [
    "# Finish your notebook with this cell\n",
    "print(starting_guess)\n",
    "print(optimized_parameters)\n",
    "print(standard_errors)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}