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": 1,
   "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": 2,
   "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": 3,
   "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": 4,
   "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": "iVBORw0KGgoAAAANSUhEUgAAAOwAAADfCAYAAAD1PT8XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJztnXl4DNn6x7+dfRUhshCCCHHFFkSWIZEwDBKRETH8LmYY\nMa6xC0OQa72xb3ExY5g7jEHsLmYQjCaWmBnbRBAROrJ09s7SSS/1+yO3OzrpNV2d7krO53n6edJV\ndeq8XalvnVPnvOd9WRRFUSAQCIzASN8GEAgE9SGCJRAYBBEsgcAgiGAJBAZBBEsgMAgiWAKBQTRL\nwVZXV2PXrl2orq4mdZO6mVU31Qzh8XhU165dKR6PR+omdTOq7mbZwhIITMWgBFtWVobQ0FC8f/8e\nAPDy5UtERUVh5MiRmDt3Lvh8vp4tJBD0i8EI9s8//8Rnn32GN2/eSLctXrwYixcvxsWLF+Hu7o49\ne/boz0ACwQAwGMEeP34ccXFxcHR0BADk5OSgrKwM/fv3BwBERkbi4sWL+jSRQNA7Jvo2QML69esB\nANT/1iLk5ubCyclJut/R0RG5ubl6sY1AMBQMRrB1EYvF9bYZGanfIaiurlY4lD548GAAQFhYGFgs\nVsMMbCCSBxKpu3nUnZeXB6BmfEYRZmZmMDMzU+t8BitYZ2dn6Y8FAC6XC2dnZ7XL79u3D7t371a4\nn6IoZGRkSL9bWlrCxcWlYcZqSIsWLRqlHlK3/usWiURgsVgIDAxUeMzs2bPx9ddfq3U+gxWsi4sL\nLC0tkZKSgv79+yMxMVHpj65LdHQ0Pv/8c7n7wsLCkJGRUW+Aa+PGjdqaTSDIEBISAoqicO7cOYXH\nqNu6AgYo2A+7LFu3bkVsbCzKysrg6uqKzZs3q30eZd2Mut0iIyMjIlaCzmCxWLCxsaHlXAYn2GvX\nrkn/7tKlC37++Wed1sdisZCYmKjTOggEujCYaZ3GxtLSEosXL4ZYLMbYsWP1bQ6BoBbNVrAuLi6k\nG0xgHM1WsAQCEyGCJRAYBBEsgcAgiGAJBAZhcNM6BGaSmlGI9/llaOtgg+6dWtF23tjYWDx9+hQC\ngQCZmZnw8PAAULMYZOLEiXLL3L59G1u2bIFYLIa9vT3Wr18PFxcXZGZmIjY2FiUlJbC3t8c///lP\ndOzYEQCwYcMGsNlsGBkZITo6GqNHj0ZWVhb+/ve/IykpCUCNu2xMTAx4PB527dqlkcMDXRDBErQi\nr6gC6w/dRzqnRLrN3dUOy6b6wNHeSuvzr127FgCQlZWFyZMn4/Tp00qPFwgEWLJkCX766Sd06NAB\nJ06cwNq1a5GQkIBvvvkG48aNQ0REBB49eoR58+bhzJkzSE5OxqNHj/Df//4XBQUF+OSTTzBs2DAA\ntU42ErFWVVUhISEBJib6kQ4RLEEr6ooVANI5JVh/6D62zw/SWb07d+7E9evXZbYFBwdj2rRpWL58\nOTp06AAA6N69O44ePQoASE1NxciRIwEAvXv3Rl5eHjgcDsRisXSxCJ/Pr9dySsQqFouxY8cOjRah\n0A0RLEEpY2POY9anvTBsoJt027nf0vHzlRdY8cXAemKVkM4pQWpGIQp5fGw+/BBHVo+AlYUpbXbN\nmTMHc+bMkbvvk08+AVAjtN27d2Po0KEAgB49euD8+fOIjIxEcnIySkpKkJ+fj4CAABw7dgyDBg0C\nn8/HokWLYG5uDqDGeX/JkiW4dOkSzp07p1exAkSwBC14n694yZhkv4W5bm6xnTt3St8tgZqua3Bw\nsHTVS1VVFRYtWgSKohAdHQ2g5j11zZo1OHLkCAYPHoxu3brB1NQUx48fh6mpKZKTk1FYWIjJkyej\nd+/eaN26NXJycmBubo4lS5Zg4cKFOHHiBExN6XvwaAoRLEEppzeG1tsWNtgdYYPdkZpRqLSsZAAq\nYGNb2u1S1sKWlpZixowZcHNzw/bt22FsbAwAEAqFSEhIgKmpKSiKwrFjx+Dq6oqdO3di4sSJMDIy\ngoODA4KCgpCSkoLhw4fDyclJ+h598+ZNxMfHIzY2lvbfoy5kWofQYLp3agV3Vzu5+7q42tE6WgzU\nLkRXxT/+8Q/06dMH8fHxUrECwLZt2/Drr78CqAlJ5OXlBTs7O3h6euKXX34BAFRUVODu3bvo0aMH\nAMiU37BhAy5cuCDTsjc2RLAakppRiGsP3qpsXZoLy6b61BOtu6sdvpnqQ3td6kSLuHPnDlJSUpCc\nnIzw8HCEh4dj+vTpAGrWPB8+fBihoaH49ddf8a9//QsAMHPmTIjFYnzyySeIiopCWFgYBg4cWK9O\nR0dHxMXFYdmyZXoLV8Si1H1sNSFCQkIAyC7lU4Wupy+Yjq7mYZlOQ+41ZZAWVk2UTV8QarrHIQM6\nELHqGCJYNUjNKFQ5fUEgNAZEsGqgzvQFgdAYEMGqQVsH5fF4VO0nEOiCCFYNGnv6gkBQBBGsmiia\nvljy9/54kp6vJ6sIzQ3i6aQmjvZW2D4/SGb6opubPVYfuIvHr/Lx7bKhaG1nqW8zCU0cRgj27Nmz\n+Pbbb8FisTBo0CDExMTozZbunVrJdIF9ejhjbFCXZi/WtPx0ZPPy4GLriG4O7rSdtyHrYS9evIh/\n//vfAAAvLy+sWbMGJiYmKCsrQ0xMDDIzM2FlZYXNmzfDza1mUQNT1sMafAb2yspKysfHhyosLKRE\nIhEVGRlJ3blzR6tzBgcHU8HBwTRZ2LzhlhVQMb+soyJ/nin9xPyyjuKWFdBaD4fDUet/VlJSQgUE\nBFAFBTX1z58/nzp+/DhFURS1bt06avfu3RRFUVRycjIVFRVFURRF3blzR/p3fn4+NWDAAIrP58vU\nKRKJqIULF1KzZ8+mBAKB2nbTfa8ZfAsrEolAURT4fD6srKwgFAphYWGhb7MUwq8WwtTYCMbGzWN4\nYNPtvcgoeiezLaPoHTbd3ov4j5fprF5F62G//vpr3Lx5E8bGxqioqEBBQQHs7GrGHq5fv47//Oc/\nAABfX18UFBQgJyeHrIelE2tra8yZMwcjRoyApaUlfHx80LdvX32bJRduUSXWHryHXl0cMC3MS9/m\n0MLEE19jer/PENzZX7rt4osknHx2ETGDvqonVgkZRe+Qlp+OosoS7Lx7EN+Fb4SVKX2vDcpW6xgb\nG+P69etYunQpnJ2d4e9fY3vdFKZt2rRBTk4Oo9bDGnwzkJaWhlOnTuHGjRu4desWWCwWvv/+e5Xl\nqqurUVZWJvdDUZTaKz80oayyGgKhGEP6taf93IZINi9Pq/3asHPnTqlzf3h4OMaOHYtdu3ZJ9w8Z\nMgT37t3DoEGDEBcXB0BxCtMP18Neu3YNR48exePHjwGg3npYgUCgsa0URSm8F8vKyhSmRZWHwbew\nbDYbvr6+sLe3BwBERETg6NGj+OKLL5SWU5VuUhfpBzu1tcPuRUNgZNS4OUh1yU+Ru+ptG9k1GCO7\nBiMtP11pWckAlG97b9rtUtTCFhYWIi0tDX5+fgCA8PBwzJ8/v8YeFxdwuVxpKyv5OyEhQafrYXk8\nHvr166dwvybpJg2+hfX09ASbzUZFRQUoikJSUpJ0raIyoqOj8fDhQ7mfdu3awdbWVif2NiWxqqKb\ngzs62cvvTXS270DraDGg3npYoVCIRYsWgcvlAgD++9//YsCAAQCAwMBAaeKze/fuwdraGk5OTjpf\nD2tra6vwXnz48KE0IoY6GLxgAwICEBYWhoiICIwZMwYCgQAzZsxQWc7MzAw2NjZyPywWq1EycWfm\nlGLNgXuo4GvejWIKiwNm1hNtJ/v2WBSg/k2oLur8zxwdHbFixQp88cUXCA8PB4fDwcKFCwEAc+fO\nRVpaGkJDQ7Fx40bEx8cD0P16WEm6SUUfTaaHyHpYHRK79zZKyqqxarovHFo27XlaXc3DMh267zWD\nf4dlMgsn9YO5qTGt0QINlW4O7kSojQARrA6xt60/X0wiMxC0gQi2kcgrqsDKfXeQxS2XbiMhZgia\nYvCDTk2FVfuTZcQKkBAzBM0hgtWQtPx03MhIVjkH+SGpGYXg5MmPSpHOKcGjV1y6zCM0cUiXWE3y\nywvr+c12sm+PxQEz4WCt/F1UVQiZq/feoneXNiptIO+/BCJYNdHGyV1VCJmBXs4y3y/dyYBAKEbY\n4JpRVyaEWC1NfQ5+djYsXFzQorsnbeelc3kdh8NBbGwsioqKYGxsjJiYGPj6+kqX3b179w6mpqZY\nvHgx/Pz8DHJ5HRGsGqTlp6t0clc2pSEJMSMv8mIXVzt81LudzLYLtzPQxbWl9Lu+MsSpQxWXi9QN\nG1Ge/lq6zdq9M7p/EwPzNqp7DarQNN1kaWkp1q9fj3PnzqFVq1ZYsGABTp8+jcjISKxZswYREREI\nCwtDeno6Jk+eDDabjT179sDDwwN79uzB69evMWXKFNy6dQsASTfJSNRxclc1B7lsqo/cVlJehPxZ\nn/aGlUXNv0adEKv67B7XFSsAlKe/RuqGjeizdZPO6tVkeZ3Eb3zMmDFSR4ZOnTpBIBCgvLxcunQO\nAN69e4eWLVvKnJcsr2MYLraOWu0H5IeYUSS0Hp1bS/9WJ8SqLgV759MouH81A05DQ2rrPH8B744l\novvypfXEKqE8/TVKU5+juqgIL7Zsh8+PB2FiRV/3XZPldQEBAQAgzQ0LAN999x26d+8OG5ua1xUj\nIyNMnToVKSkpWL16tfQ4Q1teRwSrBhInd3ndYk2d3OuGmFGFqvdfFwdrtc9FN/zsbJX7jXQUbEBV\nuknJ8rrNmzdj1apV2LJli/TY/fv3IzExET/++KPMOQ8dOoSsrCxMmDABffr0gbm5ucGlmyTTOmqi\nyMl9nv90sDPv62R9LaA8xKqNlSluPOTopF4J/iePybSuANA2dDQGHj4ECxcXpWUtXFzg4O8H/5PH\naG1dgZoW9syZM9LP6dOn8fXXX6OwsBDJycnS48LDw/HixQvp97i4OFy6dAlHjx6VLrO7e/cuCgtr\nsje0a9cO3t7eePXqFQBIl9dNnjwZrVu3li4Y0BdEsGriYN0K8R8vw5qQRZjlMxlrQhYh/uNl+Onx\nGSTc+wHvSt7rrG55IVbdnG1BiSmZ7nNj06K7J6zdO8vdZ+3uTutoMaD98rqdO3fi1atXOHLkCFq3\nrr1uV65cwYEDBwDURKV4+vQpevbsCcDw0k2S1Tpa8jgnFSwWCz2d6L055VH3/ZdXUQ1bKz1E7vsA\nXY8SS5CMEqvzP7t8+TISEhJgbGwMDw8PxMXFgcViwdfXF05OTrCxsQFFUWCxWNi7dy+srKywbNky\nZGZmwsTEBHPmzEFQUJDcOi9fvoy4uDicPXtWJtyMIuherUME24TgVVTjcvIbRAR1afQgcLqah2U6\nBrO8rrCwENeuXUN6ejpYLBbc3d0xdOjQekPizQ2+gI+y6gqV3k+6YPeJP/HkVT6G9Gvf6OtvW3T3\nJEJtBDR+DFdXV2Pjxo0ICwsDm82WjpjduHEDo0ePRnx8PPh8Pu2GMoESfilWJW3Fljv7dTYIpYyo\nod2wZPKAJr9YvjmjcQv71VdfYeTIkbh27Zo0FKSE6upqnD17FrNmzVIrsmFTw9LUEi62jhj7txGN\nEoKmLp3b1R9N5hZVwqGlBVgsFvFFbgJo/A7L4/FUBjBT5xh90lTfYetSUFKJeVtvInhAezx6yTVo\nX+SmCt33msZdYmVCZLPZKo9pblQJ1Y85Szfc4kq0amGB35/nKfRFJjALWocSFbmKNVc4JdlYcHk1\n2Jn6EYanWytER/TEm+xSufslvsjqkJpRiGsP3qp9PEE30Oqa2AxniJTytiQLViYW8GzTRW82ZOeX\nK90v8UX+9uwTuLdrieD+st5cTFja15ygVbC6GmhJSkrC7t27wefzERAQgOXLl+ukHrrx79AfPq59\nYWJU6y3T2OFAVfkiS/ZfT+GgRR0njIfPc7H+4H1UC2VTXBjK0r7miME7/7979w5xcXFITExE69at\nMWXKFNy8eROBgYH6Nk0tJGLVJmKFNqhai9u9UytQFIVvpg6AQ50ct/ef5dYTqwRNl/aREWp60Fiw\nPXr0kNuSUhQlN9mQtly9ehWjRo2Co2PNErZt27bpJ5GulugrLSOgei0ui8VCT3eHeuUULTqQ8OHS\nvsevuHjwVy6mjvqbjJcV6VLTi8aCPXToEJydnRttXWBmZibMzMwwffp0cLlcDBkyBPPmzWuUuulC\n24gV2qLJWtwPae+ofLT/w+72vWc5uPc0p16azaUJbHCLKmW2adqlJq1zLRoL9syZM3j06BH69OmD\nYcOGwc/PT6ctnkgkApvNxpEjR2BtbY2vvvoKZ86cQXh4uNJykgS98mjswTE6IlbQgaZrcdXpTksY\n5uOGvl1lF/Lff5ZdT6wS1OlSN5XWWZJuUhFmZmZqa0hjwa5btw5isRgpKSm4cuUKNm3aBA8PDwwb\nNgyBgYGwtqZ3QbWDgwP8/Pyk6SZDQkLw+PFjlYLVR7pJRaiKSFFYUdxIlmiOuqFtOrq0QEcX2Wua\nXaDeCDUACIQirD/0AFFDu8KzY802Q45lpQl0ppts0KCTkZERfHx84ONT8097+vQprly5gv379+PM\nmTMNOaVChgwZIo1UZ2VlBTabjeDgYJXloqOj8fnnn8vdFxYWRquNqlAWsQIAWlvZN6o9mtDQ7jQA\ndG2v/LgPu9S5hRVISc3F2KCanoahx7JSF4FQDFNzKzx8eF3hMU0ue92pU6dw8OBBCIVC+Pv7IzY2\nVqspJH24JioaJR7j+TH82vfTi+9xYzBv2w2FXeptH7SS/Coh0t4WwaN9S1hZmOLag7fY/vMfCs/r\n4mCNhRO90c1NtWj18Q4s6c5f/c83AICnfySrKKEeGrewH3/8sdKbS5IYl04iIiIQERFB+3kbE0nE\nClXzsGXV5djE3ovP+45HRwXJkpmEul1qC3MT9PaoXfCuav44O7+8XlZA9qMsmJkYw6dHTZxnfb4D\ny+vO04HGgpXEiSU0DFVpGX999RveFHPQwqL+CC0Tc7A2tEutbMDLvZ0dZkb0gqujrKhPJr1Eeydb\nqWDpeAduSOv8V0aBTsQKNECwkvdWoOZl+s8//4RIJIK3t3ejDuQ0VcI9h6N/215oZVkbCEBfThd0\noukINaC4dVbUQg737SiNIqnqHfinX57Du5ujdICrLuq0zmWVAthY1rbyYjGF6euvwMtdd3G2Guzp\n9OTJE8ycORMODg4Qi8XIycnBvn374O3tTad9zQ4jIyN0aCmbCeCf17chtzxfZltDnC6Y1kJr2jqP\n8Oso/VtVPOcL7NcQCMUygn3xtginbrzC7Mg+Klvns7+l48jl5zi2bqT0FdHIiAULMxOIRbobFmqw\nYOPj47Fx40ZpkGY2m434+HgcO3aMNuMINSKrK1YJGUXvkMp9he4qFhcwvYVuSOus6h3Yv6cLvLvJ\nTrelZ5Ug+Uk2Rgd0UjlC7eJgDWsLE/AqBGhhXTvKu3XuYFiYm4DDLdNJt7jB7kplZWVSsQLARx99\nhIqKClqMItSiyumCUyobzPvWm/t4VfBGZpsyt8imirJ4zl1c7TB7fF/07CLrjjmwhzNWTfdFbqHy\n+/h9fhkGdHfCwZXDZcQK1AyeAfJD09JBgwUrFouRn1/75M/Ly2uyUxMfUpr6HHlJ11Ga+rxR6lPl\ndNHBrq3M98OPT+Hh+yfS7+q4RTZV5IlGUT4jAGjVwgLe3RzVWuGk6l6XdOdbtbCoJ2ptaHCXePLk\nyRg/fjyGDx8OoCZe66xZs2gzzNBorPi7ddEkTQhFURja+SN4OXWTblPVQv/25q5a77NMe/8FdDNC\nXdclUxWmJkYwNaHP777Bgh03bhzat2+PW7duQSwWY8OGDfD19aXNMENDX1nagJo0IfLeQRcFRMsc\nx2KxEOk1WmabqhbaxEj2Fvgr7yXuvvsdk/t8ChNjE9ref/UpeG1GqDOK34BlUQGKb4VOLTsqbJ0V\nIRALIRKLNCqjjAYLtrq6GjweDx07dgQAcDgcJCYmYty4cXTZZjCUpj5XmaVNlzF51XW6kIeyFrqt\nrRPCPD+W2fY4NxV3Ob/ji35RALRfFkiH4LUVe0PKG5nxYd7jDsw/sNvc/j2MzHoBUO10IfndRZX0\n+ok3WLCzZ88Gh8NB27a171AsFqtJCladLG2NEURbldOFIhS10IsDZtbzYx7k5oOurTsBUP3+e/jR\nabzIT8fqkEUy+7JKc9DasiUsTC20Ery2YtemvKZ2iykxQEG67FReeTposGAzMzNx8eJFmWRBTRV1\nsrQZMpq00O1aOKNdixpPIVXvv2+K3oJXJbsih6IoLPl1PT7rOQZdWndUOeAlpsRoY9VaRkBiSgwW\nWGqJhqIoiCgxjFlGMgNB+eWF2HAroV6Ssoyid4i/tQebRsRKt7EzH6CsuhwjPIIAqH5QLbi0Gls/\nWSmz/f8S52Jyn08xwiNIaXltafDbcKdOnZpNhH9VWdpsPbvh9bcHwHv5qpEt04xuDu4I6uSndiut\n6v13TPcRWB4kuyyML6yCW0tXtGvhotY64M3sfWC/fSCznZ35AFHHZ6k1uv1H9jNMPDEbRXzZAaLN\nt/cpzCiYWZIlMzr+ODcV7MxaG1TZXbcuALA1t4bR/x4YqsprQ4Nb2M8++wxhYWHw9vaGiUntaTZs\n2ECLYYZG929iFIwSL0bRw9+RfeEiWvbto0cL6UfVCHXPD0ajJViaWmDd0BgAUDll5GLriJ7O3aUt\nugQxJYapkSkEYoHCspJF/+YmZrA2s6oXnsjGTPm67A+DBgR19EO5oHbuVdWDaoH/l/W27Qv7l9rl\ntaHBgt24cSN8fHzQoUMHOu0xWMzbtEGfrZvkZmkzbdkS3RYvgH0/WbdMSiwGq5FC6egKdUeo5aHO\nlJS81j6okx9cbB2x4tpmheeWiKKHY1ccHLul3v5Ir1F4nJuqsjwA/M3RQyO7VaUW/bC8sQgwEtPn\nqqhV1MSm2poqQ16WNiNTUzh8FCCzraqgAE+Xr0TXBfNg27X2htAmLaM+UjpqM0INNFzwmsw/qyrv\nwhWgJU+IYlsTZLcxVau8Ng8qAJjXLRK/r1mN9ZX0TekAWgi2a9euePz4MXr16kWnPU0G7s1bEFVU\nwLJdzSi6No4X+nLa+JCGjlBLBP/07jUUZL5Ga7fO8PINUausRDT8F6+lgrPo2llj0bTk1nZ3i9tY\nwXuF6pkMbR9UuTv2y9RLFw0WbFZWFiZMmAAnJyeZEBe6WMDORNqNHYM2gYNg8r8YV9o4XtDltKGP\nFvrDh40JgBIAf7pfVuthY1shwmeXi1CeXiTdZv26CLZ9RYAaocPkiaYltwK5O/bDRc3r5sIVwD6b\nDwsXAVA/Emw9xEIh3h79WeG8vbY0WLCLFi1SfVAzhsViwbx1zbpIVY4Xf8xdgPZRkXDw95NuL/r9\nD+RdvwGXkZ+odNqw8egCUBSMTE3lHkdHC91QsSt62Py17l/os3WT0nd8VQ+qumMElEiEnMu/okWP\nv0FUWanyuplYW8G0ZUuYylnHrc41K0t/Dd7z53AZNVJ6DMvYGO/PXlByRbRD4xGRq1evAoA0CFvd\nDwBcuXKFXisZjirHi4rMtxCW8mS2VWa9R+Hd+2o5bRQ9/APJ4yaguljWq+b9hYvIv31H6Y2viiou\nF38uWIwnS5fj5Y7deLJ0Of5csBhVXK7McaKqKvBevoJYUDuyq+xBVZHxBtkXLsps4ySewvP4zSrL\nlqe/xosdu/DH7LmyO4yMkHHwB5Q8fabWdftrzXq8P3teZnvZq3Sk7/0Wqev+pfKa8dJeIOPAIZnf\nzGKx4BiiOkhgQ9FYsM+fP8e0adNw9epVmbi/FRUVuHr1KqZOnYq//vqLViOZjirHCo95c9A6wE9m\nW6uBA9B9+VK1nDYoQTVMbG1gWifNZ97Va+D+xlbZ0gBA1umzeH9OtmUQV1fj2ao1cm/cpyviZLdl\nvMHjRUtQycmSblMlGgFP9iHFz82VPghUlRVXVUNUJRt3msViof34cbBx76zWdbN0dYWVm+wsR9nr\n18i5/AvKM97ILffhNbPr6YUOEyeAEgpljuny1QyF8/baonGXePbs2Xj27Bl27dqFhQsXwt7eHiKR\nCDweD76+vli4cCF69uypC1sZi8TxQp5wrN3d4Rg0uN52C0dHWPwvPYmyspLuad1RagBwjRyH8owM\nFN69p9A2iVtlYcpDmLWSdVPMv52MyqwsBeVyZHyoLdu2RYu/dYf4g5tXlWjsvfvKfHf7+/9B/L9G\nQFVZl9CR8IxZWG97+/G1A0qqrluPVbH19tn17Amn4cOQe/lXhXVLrplVe1dYtXeVe4xk3h7paUp/\nh6Y06B22R48e2Lt3L8rKyvDmzZuasCYdOsDGRvk6wuaMMscLXZV1CPCDWSt7cE6cVHiMRBhuf58E\nIxPZ26EqnyuviJQPfahNW9ii5wbZAH2qHlR134VNW9T2EFSVteveXaltQMOum6WLMxyDApUKVh1X\nVMm8vdnge6CEBrBaBwBsbGzg5eWl+kAaiI+PR3FxMWPnfpU5XuiyrLqiaeFZ32vJTsX/Vp0bVx8P\nKgkNvW6aPmiUYWRqCigYDGwIjAgknpycjAULFiAoKIgWweojkLg+0WaU+M8FixXeuH22qh60ksA0\nhxG65r7pvtcMXrDFxcWIjo7GyJEj8fz5cyJYLWjIjW8IThv6RNuHBd33msEndF61ahUWLFiA9+/l\nr7wgqI88t0pVaNMdbwo05JrpEp0Ils/nw8LCQuvznDhxAm3btsXAgQNx+vRpjcoaUrrJpoCh3bhM\nQq/pJuVRVVUFLpcr/dy6dQtr1qzR+ryXLl0Cl8tFcnIySkpKUFFRgXXr1mH58uUqyxpSuklC80bv\n6SY/ZN68eSgtLUW7du1QUVEBb29v2lqw77//Xvr36dOncf/+fbXEChhWuklC88bW1hbXr9OTblJr\nwW7fvh3p6elo1aoVHj58iKFDh6Kqqkrb02qNsm5Gc4ifTDAcWCwWbT4KtHSJ3d3dIRAI0KZNG1RW\nVkIsFsPc3JyOU0sZO3Ysxo4dS+s5CQSmobYvMUVRePbsmcLurqmpKXr37o3i4mLExtZ3+SIQCNqj\ntmBZLBYgdMZdAAAL70lEQVRevnwp9Rd++PCh3ONcXFwwadIk2gwkEAi1aNQlFggE2LdvH/r0qQ02\n9vbtW9y7dw8hISFo1aomVGX//v3ptZJAIADQcHldcXGxjFgBoEOHDoiMjMSVK1dQWFhIq3EEAkEW\njQRbWVmpcF9UVBQuX76stUEEAkExGgm2tLRUqWgFAsVxZAkEgvZoJNhx48Zh6dKlENZZYS+huJje\nxD8EAkEWjQTr6ekJf39/fPrpp0hKSpKJtp6WliaT4JlAINCPxo4TUVFRcHZ2xtq1a7F48WJ07twZ\nYrEYFEVh//79urCRQCD8jwZ5OgUGBmLw4MF48uQJOBwOHB0d0a9fP+LyRyDomAa7JrJYLPTq1YtE\n/icQGhFmZ2oiEJoZRLAEAoMggiUQGAQRLIHAIIhgCQQGQQRLIDAIIlgCgUEQwRIIDIIIlkBgEESw\nBAKDIIIlEBiEwefWAYCDBw/i1KlTAICePXti9erVMDFhhOkEAq0YfAv7+PFjnD59GomJiTh//jyE\nQiGOHDmib7MIBL1g8IK1s7PDypUrpYHJPT09kZ2drWerCAT9YPCCdXNzk4ZNLSgowOHDhzF06FA9\nW0Ug6AfGvAhyOBzMnDkTEyZMUCvuMUk3STAUDC7dpK5JTU1FdHQ0oqOj1c4qQNJNEgwFOtNNsigD\nb24KCwsRFhaGuLg4jbrCylpYSbrJpKQkWmwkEBQREhICiqJw7tw5hcc0qRb2hx9+QHl5ORISErB7\n926wWCwEBgZi3rx5SsuRdJMEQ8Hg0k3qkvnz52P+/Pn6NoNAMAgMfpSYQCDUQgRLIDAIIlgCgUEQ\nwRIIDIIIlkBgEESwBAKDIIIlEBgEESyBwCCIYAkEBkEESyAwCCJYAoFBEMESCAyCCJZAYBBEsAQC\ngyCCJRAYBBEsgcAgiGAJBAZBBEsgMAgiWAKBQRDBEggMggiWQGAQRLAEAoNghGAvXbqE0aNHY/jw\n4UhISNC3OQSC3jB4webn52PTpk04fPgwLl68iJSUFNy+fVvfZhEIesHgBXv79m34+vqiZcuWMDY2\nxpgxY3Dx4kV9m0Ug6AWDF2xubi6cnJyk352cnJCTk6NHiwgE/WHwqTrk5eoyMlL9nFGWDCsvLw8i\nkQghISFa20cgKCM7OxvGxsbNJ92kk5MTHjx4IP2el5cHZ2dnleVUpZtksVigKKrRE2NRFAUejwdb\nW1tSdzOo29jYGCKRiLZ0k6AMnNzcXCo4OJgqKCigqqurqS+++IK6cuWKynJVVVUUj8eT+8nOzqa6\ndu1K8Xi8RvgFsvB4PFJ3M6w7Oztb4f1YVVWl9vkMvoV1dHRETEwMpkyZAoFAgKFDh6qVJ1aTbgaB\noGtsbGxoSTlp8IIFgOHDh2P48OH6NoNA0DsGP0pMIBBqIYIlEBgEESyBwCCIYAkEBmEcFxcXp28j\n9IGxsTEGDhwIY2NjUjepmzF1syhKjisRgUAwSEiXmEBgEESwBAKDIIIlEBgEESyBwCCIYAkEBkEE\nSyAwiGYn2MYO6Hbw4EGEhoYiNDQUy5Ytg1AoxIsXLxAVFYWRI0di7ty54PP5OrUhPj4e33zzDQA0\nWt1JSUmIiIjAyJEjsW7dukat++zZsxg9ejRCQ0OxceNGndddVlaG0NBQvH//HgDw8uVLuXWVlZVh\n1qxZGDVqFCIjI/H27VvNK9PZQkADhMvlUkOGDKGKiooooVBITZ06lWKz2Tqr79GjR1RoaCjF5/Mp\niqKomJgY6uDBg9SYMWOoBw8eUBRFUTt27KC2bNmiMxvu3LlD+fr6UkuXLqUoimqUut++fUsNGjSI\nys3NpYRCITVp0iTqxo0bjVJ3ZWUl5ePjQxUWFlIikYiKjIyk7ty5o7O6//jjD2r06NGUl5cXlZWV\nRVGU4mu8bt06avfu3RRFUVRycjIVFRWlcX3NqoVt7IBudnZ2WLlyJczNzQEA3bp1Q1paGsrKytC/\nf38AQGRkpM5sKC4uxvbt2zFz5kwAQE5OTqPUffXqVYwaNQqOjo4wNjbGtm3b4OHh0Sh1i0QiUBQF\nPp8PgUAAkUgEExMTndV9/PhxxMXFwdHREYD8a3zp0iUAwPXr1xEREQEA8PX1RUFBgcbxyRixHpYu\nGjugm5ubG9zc3AAABQUFOHz4MD777DOZrpCjoyNyc3N1Uv+qVauwYMECaVet7u/XVd2ZmZkwMzPD\n9OnTweVyMWTIEAQFBTVK3dbW1pgzZw5GjBgBS0tL+Pj4wNTUVGd1r1+/HkBt7DF511hyj9Xd16ZN\nG+Tk5KgV8khCs2phqQYGdNMWDoeDKVOmYMKECdInr65tOHHiBNq2bYuBAwdKt4nF4kapWyQSgc1m\nY9OmTThx4gSePHmCe/fuNUrdaWlpOHXqFG7cuIFbt27ByMgIbDa7UeoGlF9jOq5/s2phGxrQTRtS\nU1MRHR2N6OhoTJo0CdnZ2cjLy5Pu53K5OrHh0qVL4HK5SE5ORklJCSoqKmBkZAQul6vzuh0cHODn\n5wd7e3sAQEhICJ4+fdoodbPZbPj6+krrHjt2LA4cONAodQOAs7Ozwv+vs7MzuFyutJX98G91aVYt\nrL+/P+7evYvCwkIIBAKcO3cOgYGBOquvsLAQX375JVauXIlJkyYBAFxcXGBpaYmUlBQAQGJiok5s\n+P7773H+/HmcOXMGc+bMQXBwMNatWwcLCwud1z1kyBDcvn0bPB5P2tp6e3s3St2enp5gs9moqKgA\nRVFISkrCgAEDGqVuQPn/NygoCImJiQCAe/fuwdraWmPBNqsWtqEB3RrKDz/8gPLyciQkJGD37t1g\nsVgIDAzE1q1bERsbi7KyMri6umLz5s06s6EuW7ZswYoVK3Rad69evTBjxgxMnDgRQqEQ/v7++PTT\nT9G7d2+d1x0QEICwsDBERETAzMwMXl5emDFjBkaMGKHTuj8Mn6ro/zt37lwsX74coaGhMDMzQ3x8\nvOb1UPJe7AgEgkHSrLrEBALTIYIlEBgEESyBwCCIYAkEBkEESyAwCCJYAoFBEMESCAyCCJZAYBBE\nsE2c6dOnS31bY2Nj8dtvv+msLj6fj2nTpqG8vFzhMRwOB7Nnz9aZDU0d4unUxPH09MTNmzc19llt\nCOvXr4eXlxfCwsKUHrdz5060bt1a6l9NUB8i2CbMokWLcOHCBXh4eGDnzp1YuXIlxo8fj9DQUPTo\n0QMzZ87ElStXwOPxMGfOHPz+++94+vQpxGIxduzYgY4dO4LL5WLNmjXIysqCUCjERx99hIULF9Zb\nFpaXl4fx48fj6tWrMDExwV9//YXVq1dDIBCAoiiMGzcOEydOBFCzLjQyMhJJSUkwMWlW7uxaQ7rE\nTRiJ0/mBAwfQqVMnmX0ikQhWVlY4d+6c1Cl90qRJOH36NHr16oWffvoJABATE4PRo0fj5MmTOHny\nJDgcDg4dOlSvrqtXr2LAgAFSAe7duxfh4eE4efIk9u/fjwcPHkjXgzo5OcHJyQn379/X4a9vmpDH\nWzNAUSfq448/BgC0b98ebdq0gaenp/T7q1evUFlZibt376K4uBh79uwBAFRVVcHCwqLeudLT0+Hq\n6ir9Pnz4cKxevRrJycnw9fXFihUrZFrl9u3bIz09Hf7+/rT9zuYAEWwz4MOlXx9iZmYm/Vte11Qk\nEgEAfvzxR9jY2AAASktL5UZJYLFYMhEVRo0aBT8/PyQnJ+Pu3bvYtWsXfvrpJ3Ts2BEAIBQKGyXa\nR1ODXLEmjomJCYRCYYPK2tjYwNvbG9999x2AmlHgL7/8EkePHq13bOfOnWViVc2YMQM3b97EqFGj\nsGrVKtja2oLD4Uj3v3v3Dp07d26QXc0ZItgmzrBhwzBlyhQ8efJEpqVV1OrWZfPmzXj58iVCQ0Mx\nduxY9O3bF9OmTZNbT0pKivThMH/+fBw9ehTh4eEYP348goOD8dFHHwGoGaDKz8+XiTdFUA8ySkyg\njbVr18LLywvh4eFKj9u2bRscHR3JtE4DIC0sgTbmz5+P8+fPq3ScePnypXSKh6AZpIUlEBgEaWEJ\nBAZBBEsgMAgiWAKBQRDBEggMggiWQGAQRLAEAoMggiUQGAQRLIHAIIhgCQQG8f8V1ur6W2C/WwAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c6d4c88>"
      ]
     },
     "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": 5,
   "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 = 45e3 , # Ea in J/mol\n",
    "    dH = -10e3 , # ∆H in J/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": [
    "Given $C_a, T, t$\n",
    "\n",
    "$\\frac{dC_a}{dt}=Ra$\n",
    "\n",
    "$R_a=k_r \\times C_b-k_f \\times C_a$\n",
    "\n",
    "$k_f=A*\\exp\\left(\\frac{-E_A}{R T}\\right)$\n",
    "\n",
    "$Kc=\\exp\\left(\\frac{\\Delta G}{R T}\\right)=\\exp\\left(\\frac{\\Delta H - T \\times \\Delta S}{R T}\\right)=k_f/k_r$\n",
    "\n",
    "$C_b=C_0-C_a$\n",
    "\n",
    "$k_r = \\frac{k_f}{K_c}$\n",
    "\n",
    "$\\frac{dC_a}{dt}=A*\\exp\\left(\\frac{-E_A}{R T}\\right) *(C_a-\\frac{1}{\\exp\\left(\\frac{\\Delta H - T \\times \\Delta S}{R T}\\right)} *(C_0-C_a))$\n",
    "\n",
    "From Wolfram the solution to this separable differential equation is \n",
    "\n",
    "$C_a = \\frac{C_0}{\\exp\\left(\\frac{\\Delta H - T \\times \\Delta S}{R T}\\right) +1}+C_1 * \\exp\\left(t * (A *\\exp\\left(\\frac{-E_A}{R T}\\right) * (-1) (1+\\frac{1}{\\exp\\left(\\frac{\\Delta H - T \\times \\Delta S}{R T}\\right)} )\\right)$\n",
    "\n",
    "$C_1 = C_0 - \\frac{C_0 * \\frac{1}{\\exp\\left(\\frac{\\Delta H - T \\times \\Delta S}{R T}\\right)}}{\\frac{1}{\\exp\\left(\\frac{\\Delta H - T \\times \\Delta S}{R T}\\right) }+ 1}=C_0 - \\frac{C_0}{\\exp\\left(\\frac{\\Delta H - T \\times \\Delta S}{R T}\\right) +1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could approach this trying to use curve fit..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Ca(t,A,E,dH,dS):\n",
    "    R=8.314 #J/mol/K\n",
    "    C0=10.\n",
    "    T=298.15\n",
    "    C1 = C0-C0/(np.exp((dH-T*dS)/(R*T))+1)       \n",
    "    b = C0/(np.exp((dH-T*dS)/(R*T))+1)     \n",
    "    exp = np.exp(t * 10**(A) * np.exp(-E/(R*T)) * (-1)*(1 + 1/(np.exp((dH-T*dS)/(R*T)))))\n",
    "    return b+C1*exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.odr import Model, Data, RealData, ODR\n",
    "import scipy.odr as odr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter 0: 1.0 +/- inf (1 st. dev.)\n",
      "Parameter 1: 1.0 +/- inf (1 st. dev.)\n",
      "Parameter 2: -38471.4639949151 +/- inf (1 st. dev.)\n",
      "Parameter 3: -124.10792789431792 +/- inf (1 st. dev.)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/envs/py3/lib/python3.5/site-packages/scipy/optimize/minpack.py:715: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  category=OptimizeWarning)\n"
     ]
    }
   ],
   "source": [
    "optimal_parameters, covariance = scipy.optimize.curve_fit(Ca,times,cA)\n",
    "def report(optimal_parameters, covariance):\n",
    "    \"Make this a function so we can reuse it in cells below\"\n",
    "    parameter_errors = np.sqrt(np.diag(covariance))\n",
    "    for i in range(len(optimal_parameters)):\n",
    "        print(\"Parameter {}: {} +/- {} (1 st. dev.)\".format(i,\n",
    "                                                            optimal_parameters[i],\n",
    "                                                            parameter_errors[i]))\n",
    "\n",
    "report(optimal_parameters, covariance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "hmmmm these parameters look funky...."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANgAAADOCAYAAABcr8xrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAHNpJREFUeJzt3X1UU/f9B/B3woNPoBHEB+bmWltWXFuOqyCoVRCVTAgE\nmEp1FKiKYp3tFJlYKx4qMJUHu9q1tj6ws3JEkaf0cNx2FHGNMtCtHUydrg61AgVKhBIRYpLv7w9/\nZCIPyQ33QgKf1zmcA7k3Hz835pN7c+/3+7kixhgDIUQQ4qFOgJDhjAqMEAFRgREiICowQgREBUaI\ngKjACBGQVRSYRqPBBx98AI1GQzEsLBdLiWFpuRgwK9DW1sbc3NxYW1sbxbCwXCwlhqXl0sUq9mCE\nWCvOBfbJJ59AKpUiJCQEhw8f7rG8oaEBr7/+OgIDAxEdHQ2VSsVLooRYI04FVl5ejs8//xz5+fko\nKirCV199hbNnz3ZbJzk5GeHh4SgpKYFMJkNKSgqvCRNiTTgV2LVr17Bw4UKMGzcOIpEIr776arcC\n02q1qKioQFBQEABALpejrKwMOp2O36wJsRKcCmzWrFlQKpVobW1FZ2cnSktL0dTUZFje0tICBwcH\n2NjYAABsbGzg6OhIh4lkxLLlsrKPjw9CQ0MRGRkJiUSCefPm4auvvjIs1+v1PZ7DGINYbLyONRpN\nn6dGFy5cCAAIDg6GSCTiknK3PIZLDEvKxVJi8BWnsbERAKBWq/tcx97eHvb29ibF41RgDx48wNKl\nSxEdHQ0AyM7Oxo9+9CPDcicnJ6jVauj1eojFYuh0OrS3t0MikRiNffjwYRw6dIhLOpyNHz9+2MTg\nK85wisFHHJ1OB5FIhEWLFvW5zubNm/GrX/3KpHicCqy2thbbtm1DYWEhHj58iLy8POzdu/d/wWxt\n4enpCYVCAblcDoVCAS8vL8MhY382bNiAmJiYXpcFBwcDAEpLS7mkSwhn/v7+YIxBoVD0uY6pey+A\nY4G5ublBJpMhJCQEer0eMTExmD17Nnbt2gV/f3/4+fkhKSkJiYmJOHLkCCQSCdLT002K3d9udyCH\nDYRwJRKJ4ODgwE8sxix/RrO/vz8A4Ny5c0OcCRnu+H6v0UgOQgREBUaIgKjACBEQFRghAqICI0RA\nVlNgDzu1uF5DQ66IdbGaAvv+gQYJh77A21llaLzfPtTpjHjXa1Q4d/muoB96iYmJKCoq6nd5fX29\nYP8+HzhdaLYEt+61IjW7Egd/7TvUqYxIjffbkZpdiVv3Wg2PzZw+ATujvTB54thBzaWiogKWfhnX\navZgT7p1r5UOF4fI08UF/O9Djw9paWkICAhAZGQkvvnmGwBAVlYWVq5ciZ///OeIiIhAU1MTPv74\nYzQ2NiI2NhYqlQpnzpzBqlWrIJfLERAQgMpKfvIZKKssMACo+67v0c5EGNdrVD2KqwsfH3p//vOf\nce3aNZw5cwbvv/8+7ty5A61Wi1u3buHkyZM4c+YMfvzjH6OkpAQbN27E5MmT8emnn2LixIk4efIk\nPvroIxQVFWHDhg04fvz4gHLhC+dDxOLiYnz66aeGCZcJCQndliuVSsTHx2PatGkAAHd3d6SmpvKT\n7RNcJ/EzVoyYztiHWt13arg/42R2/MrKSgQEBEAsFsPJyQkLFy6Era0tEhMTkZeXh5qaGnz55Zf4\n4Q9/aHgOYwwikQiHDh3C+fPnUVNTg8rKSpMGmA8GTgXW0dGB1NRU/OlPf8KECRMQERGB8vJy+Pj4\nGNapqqpCXFwcoqKieE+2y3PTJwzoP5KYx9iH2kA/9EQiUbfvVLa2trh//z5iYmLwxhtvQCqVQiwW\n9/je1d7ejvDwcMjlcnh6euInP/kJcnJyBpQLXzgdIup0OjDG0NHRgUePHkGr1WL06NHd1qmursaF\nCxcQGhqKN998Ew0NDbwmPHP6BCRGe/Eak5jG/RknzJw+oddlfHzo+fj4oKSkBBqNBt9//z2++OIL\niEQieHt7IyIiAs8++ywuXrxomNhra2sLnU6H27dvw8bGBnFxcfD29sZf//rXXif/DgVOe7Bx48Zh\ny5YtkEqlGDNmDLy8vDB79uxu60gkEgQFBWH+/Pk4ceIEtm3bhs8++2zAiY4fZ4/9m1+lPdcQ2xnt\n1etZRD4+9Pz9/VFdXQ2ZTAYXFxc899xz6OzsxI0bNxAcHAw7Ozu4u7sbTn74+/tj/fr1+Pjjj+Hu\n7o6AgACMHTsWnp6eqK2tHXA+fOA0XeXGjRtITEzE0aNH4eDggPj4eHh4eOCNN97o8zmenp44f/68\n0fk1/bUMoAmXlud6jQp136nhOslhWH3omTrhUpCWAUqlEt7e3pg4cSIAICwsDCdOnDAUmEajQXZ2\nNmJjY7s9z87OzmhsYy0D+JpSTvjh/ozTsCqsJ7W1teGVV17pc7lgLQNeeOEF7Nu3D5s3b8aYMWNQ\nWlqKn/70p4bl9vb2KCgogJubG3x9fZGfnw8PDw+MGjXKaGxTWgYQMhgcHR1x/vz5PpcL1jJg/vz5\nCA4ORlhYGOzt7fHSSy8hNja2W8uArKwsJCUlISMjA87Ozti3b59JsallALEU1DKAEIFQywBCrAgV\nGCECogIjREBUYMSi7Ny5EwEBAZg1axYA4Pz588jOzgbweBieqX02uxibUyY0q5sPRoa3wsJCVFVV\nGa6dXr161bDs66+/RnNz81ClZhYqMGIx4uLiAADz5s3Dw4cPUVhYiNzcXADA6NGjkZOTg/b2dvz+\n97/Hxo0bkZ6ejr/97W/Q6XSQSqWG56elpaGsrAyTJ0+GXq/H3Llzh2ybqMCIQcbFT/D3umpBYr/i\n+hK2zY/td52PPvoI7u7uKC4uRmRkJJ5//nlEREQAANatWwdnZ2dUVlZi06ZNOHXqFLRaLQoKCqDV\narFp0yZcuHABHR0dhjllLS0tQz5IgQqMWBxTLs0qlUr8+9//NsxcfvjwIW7cuIGGhoZuc8r6u0vK\nYKACIwbG9jCDxZSRO3q9HgkJCViyZAkAoLW1FaNGjUJ6enqPOWVDic4iEovCGDP8AI/vktp1C+In\nf/f29sbJkyeh1Wrx8OFDREVF4dKlS73OKRtKtAcjFkUkEhl+AGDu3LlISEiARCKBn58fPvzwQxw4\ncABbt27FnTt3EBoaCq1Wi8DAQCxevBgAeswpG0qcxyIa68nR0NCA7du3o7m5GS4uLsjMzIST08Cm\nNdBYRDJYhnQsYldPjj/+8Y8oLi7GlStXUF5e3m2d5ORkhIeHo6SkBDKZDCkpKbwkSog14rUnh1ar\nRUVFBYKCggAAcrkcZWVlhuNmQkYaXntytLS0wMHBwdAyy8bGBo6OjlCpVHBxceE3c0KsAKc92I0b\nN1BQUICysjJDx59jx44ZlvfWyYcxBrHY+D+j0WigVqt7/XnyrBIhQmOM9fleVKvVffaO6Q2vPTmc\nnJygVquh1+shFouh0+nQ3t4OiURiNDb15CCWwmJ7ctja2sLT0xMKhQJyuRwKhQJeXl4mdVmlnhzE\nUlh0T46kpCQkJibiyJEjkEgkJk8voJ4cxFJQTw5CBEI9OQixIlRghAiICowQAVGBESIgKjBCBEQF\nRoiAqMAIERAVGCECogIjREBUYIQIiNNYxNzcXOTm5hruBl9XVwc/P79u9wBTKpWIj4/HtGnTAADu\n7u5ITU3lN2tCrASnAouIiDA0grx9+zbWrVuH+Pj4butUVVUhLi4OUVFR/GVJiJUyu6tUcnIytmzZ\n0mOmcnV1NTo7O1FUVARXV1fs3r0bU6ZMGXCihFgjs76DXblyBc3Nzb3O05JIJFi7di0KCwuxYMEC\nbNu2bcBJ8uF6jQrnLt/F9RrVUKdCRhCz9mAnTpzoc3JkWlqa4ffXXnsNmZmZUKvVRufXaDSaPqdi\nD2RGTeP9dqRmV+LWvVbDYzOnT8DOaC9MnjjW7Lhk+OpqGdCX/uYuPo1zgWm1WpSXl+O9997rsUyj\n0SA7Oxuxsd1bMHfdiqY/QrUMeLq4AODWvVakZlfi4K99zYpJhrchaxkAADdv3sSMGTMwdmzPT397\ne3sUFBTAzc0Nvr6+yM/Ph4eHB0aNGmU0rhAtA67XqHoUV5db91pxvUYF92cG1hSVDD9D1jIAAO7e\nvQtXV9dujz3ZMiArKwtJSUnIyMiAs7Nzt1P4/RGiZUDdd33v5ruWU4GRp/HZMoBzgUmlUkil0m6P\n7d271/C7u7s7Tp06NfDMeOA6qf8XydhyQgZqWI/kcH/GCTOnT+h12XPTJ9DeiwhuWBcYAOyM9upR\nZDOnT0BitNcQZURGkmF/+6LJE8fi4K99cb1Ghbrv1HCd5EB7LjJohn2BdXF/xokKiwy6YX+ISMhQ\nogIjREAj5hCRD/Q9jnBFBWYCGs9IzEWHiCbobzwjIf2hAjPClPGMhPSFCswIU8YzEtIX3ntyNDQ0\nYPv27WhuboaLiwsyMzPh5GS9JwRoPCMZCE57sIiICBQVFaGwsBAHDx7E+PHje/TkSE5ORnh4OEpK\nSiCTyZCSksJrwoONxjOSgTD7ELG3nhxarRYVFRUICgoCAMjlcpSVlUGn0w080yFE4xmJucw6Td9X\nT46WlhY4ODgY7slsY2MDR0dHqFSqHs1xniZUywA+0HjGkWVIWwYAfffk0Ov1PR5jjEEsNr6jFKpl\nAJ9oPOPIMKQtA/rryeHk5AS1Wg29Xg+xWAydTof29nZIJBKjcYVoGUCIOYa0ZUB/PTlsbW3h6ekJ\nhUIBuVwOhUIBLy8vwyFjf4RoGUCIOfhsGcD5JEdfPTm6Kj4pKQnFxcUICgpCfn4+3n33XV4SJcQa\nidhQn0Ewgb+/PwDg3LlzQ5wJGe74fq/RSA5CBEQFRoiAqMAIERAVGCECogIjREBUYIQIiAqMEAFR\ngREiICowQgREBUaIgDgN9i0tLcWhQ4fQ0dGB+fPn45133um2XKlUIj4+HtOmTQPw+FZGqamp/GVL\niJUxucC++eYb7NmzB6dPn4azszOioqJw4cIFLFq0yLBOVVUV4uLiEBUVJUiyhFgbkwvs7NmzCAwM\nxOTJkwEAWVlZPaaXVFdXo7OzE0VFRXB1dcXu3bsxZcoUfjMmxIqY/B3szp070Ol0WLduHUJCQpCT\nk4MJE7r3qZBIJFi7di0KCwuxYMECbNu2jfeECbEmJu/BdDodlEolcnJyMG7cOMTFxaGoqAhyudyw\nTlpamuH31157DZmZmVCr1SZNXrPknhxkZBmSnhyTJk2Cj48PJk6cCODxvJmqqipDgWk0GmRnZyM2\nNrbb8+zs7EyKbw09OcjIMCQ9Ofz8/JCQkIC2tjaMHTsWSqUSixcvNiy3t7dHQUEB3Nzc4Ovri/z8\nfHh4eGDUqFEmxTfWk6OpXYXVeaZtlCVjjIExQCSiVgiWqPFBM1wcnQa/J8fLL7+M2NhYrF69Glqt\nFvPmzUN4eDh27doFf39/+Pn5ISsrC0lJScjIyICzs3O3jr/GDPueHIxBq2fdDndFIhFsxaLH1UYs\nBp89OahlwCB5O6us15tIzJw+AQd/7Tv4CZFeUcsAK0R3aBm5qMAGAd2hZeSiAhsEdIeWkYsKbBAI\ncYeW6zUqnLt8lw4vLRzdo3mQ7Iz26vU+z1zv0EL3i7YuVGCDhK87tPR3v2g6G2l5qMAG2UDu0GLK\n2Ui6+4tloe9gVoTORlofKjArQmcjrQ8VmBWh+0VbH04FVlpairCwMCxfvrzXm5s3NDTg9ddfR2Bg\nIKKjo6FS0SlkvtH9oq0Lry0DkpOTER4ejpCQEOTn5yMlJQUZGRmCJD5S0f2irQtvLQO0Wi0qKirw\nu9/9DgAgl8uRmpoKnU5n0h0uCTd83C+ailR4JhfYnTt3YG9vj3Xr1qGpqQl+fn54++23DctbWlrg\n4OBgKCYbGxs4OjpCpVLBxcWF/8yJ2ehi9eAx+TtYV8uAAwcOIC8vD9XV1SgqKjIs1+v1PZ7DGINY\nbNo/odFooFare/15PEnR4mfVWI3+LlaT/7UM6Ounr9YWveGtZYCTkxPUajX0ej3EYjF0Oh3a29sh\nkUhMik8tAwYHXaw2ziJbBtja2sLT0xMKhQJyuRwKhQJeXl4mf/8y1jKA8MOUi9UjvcAcHR0ts2VA\nUlISEhMTceTIEUgkEqSnp5ucyLBvGWAh+L5YPRxPlFDLADIgfbUveG76BGSZOGB4uJ4ooZYBZMD4\nuFhNJ0pMQ6PpR6CBXqymEyWmowIbwcy9WE0nSkxHh4iEMxrVbzoqMMIZjeo3HRUYMQuN6jcNfQcj\nZhmuo/ofafXQ6noO+zMXFRgZkOEyqr/rup7q+w5e41KBkSFjSRere7uuxwf6DkaGjKVcrO7vut5A\nUYGRIWFJN8QQshsX50PErVu34tq1axgzZgwA4M0338SSJUsMy5VKJeLj4zFt2jQAgLu7O1JTU3lK\nlwwXlnSxWsjrdpwL7OrVq8jLy+tzflZVVRXi4uIQFRU14OTI8CXExWpzT5Z0XdcT4jCRU4G1tLRA\npVIhISEBdXV1WLZsGTZv3txtnerqanR2dqKoqAiurq7YvXs3pkyZwmvSxPr196bmerGaj5MlXfcO\nqDH5XzUNp+9gzc3NmD9/Pvbv349Tp07h8uXLKCgo6LaORCLB2rVrUVhYiAULFmDbtm28JkyGD74u\nVvNxsqTrup7T+NEYP870CZXGDGg+2NmzZ6FQKAydpHrj6emJ8+fPG53AptFo+ux10DWjubS01NxU\niQUbyHWw6zUqJBz6os/l+ze/yimmv78/GGNQKBR9rtPf5OCncTpE/Ne//mXoKAU8bnTzZEsAjUaD\n7OxsxMbGdnuenZ2d0djUk2PkGsjFaiFOlgxJTw4AePToEdLS0uDl5QU7Ozvk5ubiF7/4hWG5vb09\nCgoK4ObmBl9fX+Tn58PDwwOjRo0yGpt6chBzCHGyZEh6cgDA7NmzsXr1aqxYsQI6nQ5SqRTLly/v\n1pcjKysLSUlJyMjIgLOzM/bt22dSbOrJQczB58mSLtSTg5An8Dnkiu/3Go1FJFbPkkf2U4GRYYOP\nkf18o7GIhAiICowQAVGBESIgKjBCBEQFRoiAqMAIERAVGCECogIjREC8twxoaGjA9u3b0dzcDBcX\nF2RmZsLJybIu/hEyWHhvGZCcnIzw8HCEhIQgPz8fKSkpyMjIGHCihFgjToeIT7YMCA4O7jF/S6vV\noqKiAkFBQQAAuVyOsrIy6HQ6/jImxIrw2jKgpaUFDg4OhkmYNjY2cHR0hEo1eC24CLEknA4RZ86c\niYMHDxr+joyMhEKhQFhYGIDHM5yfxhiDWGy8jvtrGdDY2AidTmeYSkCIUOrr62FjYwO1uu+Z0kPW\nMsDJyQlqtRp6vR5isRg6nQ7t7e2QSCRGYxtrGSASicAYM3vyJWMMbW1tcHR0tPoYlpSLpcTgK46N\njQ10Oh1vLQPAOPjHP/7Bli5dytRqNevs7GQxMTGspKSk2zobN25khYWFjDHGCgoK2KZNm0yK3dnZ\nydra2nr9qa+vZ25ubqytrY1Lut20tbUNmxiWlIulxOA7l/r6+j7fj52dnSbH471lQFJSEhITE3Hk\nyBFIJBKkp6ebFJvLbpcQoTk4OPDSNoDzafro6GhER0d3e2zv3r2G36dOnYrjx48PODFChgMayUGI\ngKjACBEQFRghAqICI0RANnv27Nkz1EmYwsbGBnPnzu123W0kx7CkXCwlhqXlAlhJ41FCrBUdIhIi\nICowQgREBUaIgKjACBEQFRghAqICI0RAFltgarUaMpkMdXV1AID//Oc/WLVqFZYvX4633noLHR0d\n/T7/+PHjkMlkkMlk2LlzJ7RaLW7evMkpRnp6OgIDAyGTyZCdnQ0AnGN02bdvHxITE82OsXXrVkil\nUoSGhiI0NBRnz57l/JqUlpYiLCwMy5cvR0pKCudccnNzIZfLERoaCrlcDi8vL/zmN7/hnAcAFBcX\nIygoCDKZDPv37zfrdfnkk08glUoREhKCw4cPc4ph6vtLrVZj06ZNCAwMxIoVK3D37l2j29aN2RNn\nBPTll1+yoKAg9uKLL7La2lrGGGMhISHs8uXLjDHG3n//fZaRkdHn8//5z38ymUzGOjo6GGOMJSQk\nsOPHj3OKUVZWxn75y18yvV7POjo62OLFi9l///tfTjG6XLp0iXl7e7MdO3Zw3pYuy5YtY62trd0e\n4xLn7t277NVXX2UNDQ1Mq9WyNWvWsLKyMrNyYYyxmpoa5u/vzxoaGjjHePjwIfPy8mIqlYrpdDq2\nYsUKdunSJU5xLl26xIKCgpharWZ6vZ7FxcWxv/zlLybF4PL+SklJYYcOHWKMMVZeXs5WrVpl0uvT\nxSILLDExkV25coUtXryY1dbWsvr6eubv729YXldX1+3vp92+fdvwYjHG2NGjR9mOHTs4xWCMMa1W\nyxhj7N69e8zX15dzHowxdv/+fbZy5UqWnZ3NduzYYXaMOXPmsA0bNjCZTMY++OADznGOHTvGfvvb\n3xr+bmxsZLW1tZxz6RITE8OKi4vN2h61Ws08PT1ZXV0d6+joYGFhYayyspJTnCNHjrD9+/cb/s7J\nyWGRkZEmxTDl/bVkyRLGGGNLlixhdXV1hmVLlixh9fX1/W7fkyzyEDE1NRWvvPIK2P8PMmloaMCU\nKVMMyydPnoyGhoY+nz9jxgzMmTMHwONGPZ999hmeffZZTjGAx0NmDh48iKCgIHh7e+Pbb7/lHCMp\nKQlbt241tLnjui1d2/B0s6Hi4mJOce7cuQOdTod169YhJCQEOTk5aGxs5JwLAFy5cgXNzc0IDg42\na3vGjRuHLVu2QCqVYtGiRfjBD34AOzs7TnFmzZoFpVKJ1tZWdHZ2orS01OQYpry/vv32216Xubi4\nGJaZwiIL7Gm9NdMxpZHOvXv3EBUVhYiICEPBcY3x9ttvo7y8HHV1dbh48SKnGHl5eXB1dcXcuXMN\nj5mzLV3NhsaPH4/Ro0cjMjISly5d4hRHp9NBqVTiwIEDyMvLQ3V1NSoqKjjnAgAnTpxATEwMAPO2\n58aNGygoKEBZWRm++OILiMViKJVKTnF8fHwQGhqKyMhIrF+/HnPmzIFWq+Wci7FtMPe918UqbiE7\ndepUNDY2Gv5uamrC1KlT+33O9evXsWHDBmzYsAFr1qxBfX09pxhff/019Ho93NzcMHr0aCxbtgxX\nr15FU1OTyTHOnDmDpqYmlJeXo7W1Fe3t7RCLxZxiAL03G+p6rqlxJk2aBB8fH0ycOBHA45t9d8Xl\nkotWq0V5eTnee+89AOb93yiVSnh7extyCQ0NxdGjRznl8uDBAyxdutQwuz47Oxuurq74+9//zikX\nY9swdepUNDU1GfZiT/5uCqvYg02bNg1jxozBlStXAACnT5/GokWL+lxfpVJh/fr12L17N9asWWNW\njFu3bmHPnj3QarXQaDQ4e/YsFi5ciNGjR5sc49ixY/j8889RVFSELVu2YPHixUhJSeEUAwAePXqE\ntLQ0PHjwABqNBrm5uVi1ahWnOH5+frh48SLa2toMe7Of/exnnHO5efMmZsyYgbFjxwLg/roCwAsv\nvAClUon29nYwxlBaWgpPT09OudTW1mLjxo3QarVoa2vD6dOnsXLlSs7bY2wbfH19cfr0aQBARUUF\nxo0bx6nALHoP9mTrrczMTOzatQtqtRrTp0/vt5nOH/7wBzx48AAffvghDh06BJFIhEWLFnGKERAQ\ngGvXriEkJAQ2NjZYvnw5pFIpZs6ciXfffdekGH3JyMjgFKOvZkPPP/+8yXFefvllxMbGYvXq1dBq\ntZg3bx7Cw8Ph4eHBKZe7d+/C1dW122NcXlcAmD9/PoKDgxEWFgZ7e3u8+OKLiI2NhVQqNTkXNzc3\nyGQyhISEQK/XIyYmBrNnz+b02pry/nrrrbfwzjvvQCaTwd7eHvv27et323r8G4zRdBVChGIVh4iE\nWCsqMEIERAVGiICowAgREBUYIQKiAiNEQFRghAiICowQAVGBESIgKjBCBEQFRoiA/g/S+LfVv+1J\nIQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11183f400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the data\n",
    "plt.plot(times, cA, 'o', label='data')\n",
    "# Make a new x array with 50 points for smoother lines\n",
    "x_many_points = np.linspace(times.min(),times.max(),50)\n",
    "# Plot the fitted curve\n",
    "plt.plot(x_many_points, Ca(x_many_points, *optimal_parameters), label='fitted')\n",
    "# Add the legend, in the \"best\" location to avoid hiding the data\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot looks odd... Also How do we introduce the additional data sets.  Curve_fit will only take 1 data set...\n",
    "\n",
    "A search online found the following resource  http://stackoverflow.com/questions/26242979/fit-multiple-parametric-curves-with-scipy\n",
    "\n",
    "Will attempt to implement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def function(t,C0,T,logA,E,dH,dS):\n",
    "    R=8.314 #J/mol/K\n",
    "    C1 = C0-C0/(np.exp((dH-T*dS)/(R*T))+1)       \n",
    "    b = C0/(np.exp((dH-T*dS)/(R*T))+1)     \n",
    "    exp = np.exp(t * 10**(logA) * np.exp(-E/(R*T)) * (-1)*(1 + 1/(np.exp((dH-T*dS)/(R*T)))))\n",
    "    return b+C1*exp\n",
    "\n",
    "t=array([10,  20,  30,  40,  50,  60,  70,  80,  90, 100,10,  20,  30,  40,  50,  60,  70,  80,  90, 100,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,7.230,  6.073,  5.452,  5.317,  5.121, 4.998,  4.951,  4.978,  5.015,  5.036,5.137,  4.568,  4.548,  4.461,  4.382, 4.525,  4.483,  4.565,  4.459,  4.635])\n",
    "n1=10\n",
    "n2=20\n",
    "\n",
    "def leastsq_function(params, *args):\n",
    "    logA = params[0]\n",
    "    E = params[1]\n",
    "    dH = params[2]\n",
    "    dS = params[3]\n",
    "    t = args[0]\n",
    "    cA = args[1]\n",
    "    C0 = args[2]\n",
    "    T = args[3]\n",
    "\n",
    "    Ca_fit = np.empty(t.shape)\n",
    "    Ca_fit[:n1] = function(t[:n1], C0[0], T[0], logA, E, dH, dS)\n",
    "    Ca_fit[n1:n2] = function(t[n1:n2], C0[1], T[1], logA, E, dH, dS)\n",
    "    Ca_fit[n2:] = function(t[n2:], C0[2], T[2], logA, E, dH, dS)\n",
    "\n",
    "    return cA - Ca_fit\n",
    "\n",
    "\n",
    "\n",
    "params0 = [6, 45e3, -10e3,-50]\n",
    "args = (t, cA, (10, 10, 10), (298.15, 308.15, 323.15), (n1, n1+n2))\n",
    "pfit, pcov, infodict, errmsg, success = scipy.optimize.leastsq(leastsq_function, params0, args=args,full_output=1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are reasonable parameters given the initial guesses!  Lets work to plot them against the 3 data sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "logA=pfit[0]\n",
    "E=pfit[1]\n",
    "dH=pfit[2]\n",
    "dS=pfit[3]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM4AAADOCAYAAAB2Hz3EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJztnXl8U2W6x78nS9M23feW0tIWSssOsqNsRekAZXMQrsoI\njrI4DnqZkRGdudc7c9URcdSRGVBHZRZGGJXFZVwuskihsgmydINSSvd9SdqmaZJz/0hboXRJ0rRJ\n2/P9fPJpk5Pz5Elyfnnf855nEURRFJGQkLAKmaMdkJDojUjCkZCwAUk4EhI2IAlHQsIGJOFISNiA\nJBwJCRtwqHD0ej1vvPEGer3e4XYkG87ri7PYuAXRgWg0GjE2NlbUaDQOtyPZcF5fnMXGzUhTNQkJ\nG7BaOG+99RaJiYksWrSIN998szt8kpBweqwSTkpKCp988gkfffQR+/fv5/z58xw8eLC7fJOQcFqs\nEk5qairTp09HrVYjCAJ33XWXJByJfolVwhk2bBjJyclUV1fT0NDAoUOHKC0t7S7fJCScFkEUrYuO\n3rlzJ3v37sXHx4epU6dy/vx5duzY0e7z9Xp9u0uA06dPp7a2lgEDBiAIgnWe34QoiuTn53fJjmTD\neX2xh42SkhL0ej1Hjx7Fw8Ojzee4uLjg4uJikT2FNS9eW1vL3XffzapVqwCziCIiIjrc580332Tb\ntm3WvIxNeHl5STbsbMNedpzBhtFoRBAEZsyY0e5zHn/8cX7+859bZM+qESczM5Nf/OIX7Nu3j/r6\nelasWMH//u//Mnbs2Hb36WjEWbhwIQCHDh2y1AUJCZtISEhAFEU+/vjjdp/TbSNObGwsSUlJLFq0\nCJPJxOrVqzsUTWfOdGUaISFhLYIgtDtNsxarhAOwZs0a1qxZY5cXl5DorUiRAxISNiAJR0LCBiTh\nSEjYgMOFU99gIC27wtFuSEhYhcOFU1OrZ9O2Yzz56hFKKusc7Y6EhEU4XDjNZOVV88LOU452Q0LC\nIpxGOGAWjzRtcyxp2RV8ffpGt38PmzdvZv/+/R1uLyws7FYfuoLV13G6wwEFYGi6X1CmJT7Kz4Ee\n9U9KKut4YecpsvKqWx6LCffmmVUTCfJ173F/Tp48iZVhlD2Kw0ccORDBDxEEYQH2ubIrYR2tRQP2\nnz6/+OKLzJ07l5UrV5KbmwvAq6++yn333cePfvQjVqxYQWlpKTt27KCkpIQ1a9ZQUVHB559/zvLl\ny1m8eDFz587l1CnHT+kdLhwR8EdADQwO95ZGGweQll1xm2iasdf0+csvvyQ1NZXPP/+c119/nZyc\nHAwGA1lZWezZs4fPP/+cQYMG8dlnn7Fu3TqCgoJ4++238fX1Zc+ePWzfvp39+/ezdu1a3nvvvS77\n01UcPlUzNv2NdVHws4cmONSX/kpBmbbT7V39QTt16hRz585FJpPh5+fH9OnTUSgUbN68mQ8++IDs\n7GzOnTvHwIEDW/YRRRFBENi2bRuHDx8mOzubU6dOIZfLu+SLPXD4iKNWu6D2ckWhN1FyvdLR7vRL\nOpse22P6LAjCLecsCoWCyspKVq9ejclkIjExkTlz5tx2XlNXV8e9995LXl4eEyZMYOXKlU5x7uNw\n4bipFMxfOgKArz9No1Fv6GQPCXsTH+VHTLh3m9vsNX2eMmUKn332GXq9npqaGo4dO4YgCEyePJkV\nK1YQHR3N8ePHMZlMgFlYRqOR69evI5fLWb9+PZMnT+abb75peY4jcbhwAIaOCCEyxp+aah0nDmc5\n2p1+yTOrJt4mnphwbzavmmgX+wkJCUyePJmkpCQee+wxBg8eTENDAxkZGSxcuJBVq1YRHx/fsmiQ\nkJDAo48+iru7O/Hx8cydO5elS5eiVqvJz8+3i09dwerUaXuSkJAAwNdff01RfjVvvfoNCoWMx5+e\njZePm6Pc6tekZVdQUKYlLMCjTy3U3Hys2QOnGHEAQgZ4M3ZiBIZGE1//O83R7vRb4qP8SJgQ0adE\n0x04jXAAZv0oDheVgotn88nLkRYKJJwXpxKOh6eKu+YMAeDLA5edYvVEQqItnEo4AJOmR+Hr705+\nTiWXzxU42h0JiTZxOuEoFHLmLBgGwMHPUqXlaQmnxOmEAxA3sml5ukrHiSPXHO2OhMRtOKVwBEHg\nnoXDQYATh69SU1XvaJckJG7BauEcOHCABQsWkJSUxJYtW7rDJwBCw70ZOyGCRr2RQ/9O77bXkbiV\njLIsjmSnkFHWvRei+1U+jk6n44UXXuCLL77A29ubFStWkJKSwpQpU7rFuVnz4rj8fQEXzuYx4c5B\nDIjw7ZbXkYCy2gpePr6D7MrclseifAfy1LR1BKh7/ppOn8rHMRqNiKKITqejsbERg8GAq6trd/l2\n6/L0fml5ujtpLRqA7MpcXj7efkF9a+m3+ThqtZoNGzaQmJjIjBkzCA8P77QEbleZdFcUPn7u5OVU\ncvm8tDzdHWSUZd0mmmayK3PtMm3ra/k4VgknIyODvXv3cuTIkZbo1nfffbfDffR6PVqtts2bKIqd\njiIKpZy7k+IBOPhpKo2Nxg6fL2E9hZqSLm23hM7ycV566SXOnTtHbW1tyz435+McP36cP/7xj+zd\nu5e6OtuqIYmi2O6xqNVqrepIbdU5TnJyMpMnT8bX13yusXTpUt5//30efvjhdvfprM2HJe0b4kaG\nEhHtx41rFaQcyWL63bHWuC3RCaGeQV3abgkd5eM8/PDDJCYmIpPJ2s3HWbx4MRMmTGDo0KHs2rXL\nJh80Gg133HFHu9utafNhlXDi4uJ46aWXePzxx3Fzc+PQoUMMHz68w33Wrl3L6tWr29zW3OajMwRB\nYO6i4bz92jGOH7rKmIkD8fKWoqftxdCAGKJ8B7Y5XYv2jWBoQEyXX2PKlCm8/fbbLF++HJ1Ox7Fj\nx4iKimrJx9FoNBw/fpzZs2cDbefjiKLIr3/9a5vzcTw9PTl8+HC72y1t8QFWCmfatGksXLiQpUuX\n4uLiwogRIzrtXGCvNh+h4T6MmTCQ86dyOfTvdBb/R/eeW/U3npq2rs1VtV9OW2sX+wkJCVy8eJGk\npCQCAwNvy8dRKpVt5uPs2LGjJR/H3d2dCRMm2JyPY882H06Tj2MJ2hod235/CH2DkZ8+cRcDIny6\n071+SUZZFoWaEkI9g+wy0jgLfTYfxxI8vFy5M6E5evqStDzdDQwNiGFm1JQ+JZruoFcJB2Dy9Gh8\n/NzIu15JqrQ8LeEgep1wFMqbo6fTpOVpCYfQ64QDED/KvDxdXVlPyhGpuIdEz9MrhXNz9PTxQ1fR\nVOsc7ZJEP6NXCgcgbKAPY8YPbIqelop7SPQsvVY4YI6edlHJ+f5MHvk3pOIe9qAmLZ2SQ4epSeve\nVI5+lVbgbHg2LU8f+nc6e//xHT994i7c1ZZf/ZX4gYbSUtJe3EJt1g8Zt+qYaOI3b0IVGNjj/vSp\ntAJnZMrMGCJj/Kksr+PDv53FaHR8edTeSGvRANRmXSPtRfslK/bbtAJnRC6XMXpmDCq1C9evlvF/\nH6c62qVeR01a+m2iaaY265pdpm19La2gV0/Vbu4i5gbEI3AqORtXLxUzmyIMJDpH18m5hK6wEK/4\nuC69htTmw4m4uYtYPXAN85z46L/TyblW7kDPeheuoaFd2m4JUpsPJ6GtLmJVQB4mBGD3u6epqpDa\nv1uCV3wc6pjoNrepY2K6PNqA1ObDaWivi1ghUIFIQ30je947jb5BKmhoCfGbN90mHvOq2lN2sS+1\n+bAjXQn1TsuuYNO2Y21ukwEJgZ5UldYSPyqUH6+8A0Fmee5Pf6YmLR1dYSGuoaF2GWmchX6dVnAz\nHXURiw735qG1k1F7uJB2oZBvDl7pYe96L17xcQTNntWnRNMdOFw4DQ01pBddsmnfjrqIefu6s2zV\nBGRygaNfZpB2wXmvQkv0Phy/HG0ycOjU27znGmx18bsgX3de+8+Z7XYRi4jyY/69o/jkX9+z//1z\n+AWoCQ7rvDiIhERnOHzEARijUpLfheJ3HXURGzspgkl3RdGoN7L73VPUahu66q6EhOOFYxBFVILA\nKJXSbsXvWnN30jCihgRQXVnPB389g9Fw+3JmWnYFX5++QVp2hd1fX6Lv4fCpmq5pTW+8SsnZhkYK\nNSV2z3eXyWX8+Cd38M7rydy4VsEX+y8x/8ejgFujD5qJCffmmVUTCfJ1t6sfEn0Hh484RkQKDUa8\n5TKGKhV2KX7XFm7uLix/eAIuKgVnU3I4c+I6wG2iAcjKq+aFnY4PJJRwXqwSzu7du1m8eDFLlixh\n8eLFTJw4kV/96lddduK0rhGAuzzUxPq3fQXbHgQGe7L0wXEgwBf7LnH02LXbRNNMVl61NG2zI888\n8wxz585l2DBzvYjDhw+zc+dOAC5cuMDWrVutstdZPk93Y9VUbcWKFaxYsQKA69ev88gjj/DLX/6y\ny06kNxpIQIYfRrRV2Xj6dp94YocFkzAvnq8/S+P4v9NRAe0tFxSUaaW25XZi3759XLhwAaVSCcDl\ny5dbtl29epXy8t4VW2jzOc5vf/tbNmzYQGAXk5y8VJ78NuGX+GjzyMv8lOLrR7tVOABTZ8VQXFDD\npXP5DEYgDZG2op/CAuxT9bG/s379egCmTp1KfX09+/btY/fu3QC4urqya9cu6urq+POf/8y6devY\nunUr3377LUajkcTExJb9X3zxRY4cOUJQUBAmk4lJkyY57D3ZJJwzZ85QXl5uce3njnBVqBgaEIPR\nZwAF1w5SXZqGrrYUV3X3ZR0KgkDS8tFUlGkpyK0mGrjKrZFHg8O9+8xo88Ffz5B5ubhbbMcOD2bZ\nQ+M7fM727duJj4/nwIEDrFy5kiFDhrTMXB555BH8/f05deoUjz32GP/6178wGAzs3bsXg8HAY489\nxtGjR9HpdC35PFVVVXY59rqCTcJ5//332y2k3hq9Xt9u+4Sbw+TkClcCB0ykOOcbinO+IXLYvba4\nZjFKpZz7Vk/grT98A1o9YUBBk3iaow8k7IslYZHJycmkp6e3ZHnW19eTkZFBcXHxLfk8M2bMsOn1\ntdq2g4Oh4zrnrbFaOAaDgZSUFH73u99Z9Hxr2nwERdxJ8Y1kygvOMGBwIgoXtbXuWYWXtxsrHp7I\nX/98ggEGE7MmRzJ6/MA+M9I009mI0FNYUmTfZDKxadMm5syZA0B1dTUqlYqtW7fels9jLQ5r8wGQ\nmZlJZGQk7u6WXeOwps2Hi5svvsGjqCw6T2nuCUJj7rbWPasJj/RlwbJRHHj/PJknbxAZ5IE4yNeq\nTgoSndPcRKz54JfL5S0zEblcjtForsg6efJk9uzZw8yZM2lsbOShhx5iw4YNbbYJsbYboMPafADc\nuHGDsLAwi59vbZuP4EEzqCw6T0nuCYIHzUQmV1rrotWMbqrP9vm+S/zfx6mUl2j50dKRyOUOv8zV\nZxAEoeUGMGnSJDZt2oSPjw+zZs3iT3/6Ey+//DIbN24kJyeHJUuWYDAYmD9/fkvPnNZtQmzxoU+3\n+cg4vR1t5TUihy0jILznzjWuppfw0d/P0qAzEDUkgGUPjcfVrfuFK9H99It8nOBI84lfcc43PZpf\nPjguiNU/vxMfPzeyr5Tx7hvJVJbXdr6jRL/DKYXjHRiHyj0QXW0xNWUZPfraQSGe/HTDXQyI9KWs\nWGuOb5MiCCRa4ZTCEQQZwZF3AVCcc7THX1/tqeIn66cwfEwYdbV6/r49hYtn83rcDwnnxSmFA+Af\nNh6FUo2m4ip1NT1fnEGplLP0wXFMvzsWo9HEvn+e48iXGU5RmkjC8TitcGRyJYEDpwDmcx1HIAgC\nMxOHsvj+scjlMr75KpN9u85hkJpZ9XucVjgAgQOnIcgUVBSdR69rO4q5Jxh1Rzgr103GXe3CpXP5\n/G17CrUaKZO0P+PUwlGqPPAPHQeiiZIbxx3qS0S0Pw9vuJOAIA/ycip554/HKC3SONSn3kRfSytw\nauEABEVOB6As71uMBsf+yvsFqHl4w51EDQmgqqKed99IJiujxKE+9Rb27dvHp59+SmqquSj+5cuX\nW+LG+lVaQU/h5hGMV0AcNWXplOefJijyTof64+qm5P5HJ/H53ot89+0N/vmXU/xoyQjGTx3kUL+c\nGSmtwEEER86gpiyd4hvHCIyYiiA4dqCUy2XM//Eo/IM8+L9PUvn3RxcpL9Vyd9JwZE5YMTTr/N+o\nLu2e9ifegcOIGfOTDp8jpRU4CE+/GNw8w6jXFFBVfAnfkFGOdglBEJgyIwY/fzV7d33HyW+yqSir\n494Hx+Gi6hUfa4/j6LQCe9IrvmFBEAiOnM71S7spzjnqFMJpZuiIEFb9bCq73znNldRi3nn9GD+6\ndySDYgIc7VoLnY0IPYWj0wrsidMvDjTjFzIGpcqb2uobaKuuO9qdWwgN9+GnT96Jb5AHpcVa/vbn\nFD76+1lqquod7ZrT0FZaQXMqQVtpBQaDgfr6eh566CFOnDjRZpsQR9IrRhwAQSYnKOJO8q98RvH1\no3iMGeRol1poqc1WUkMQMACBy+cLyLhczPS7hzB5RjQKheO7iDkSZ0grsCdOmVbQHobGei5+8zwm\no54Rd25C5e4c06EnXz1yS5kpBRCOQCDmg8QvQM3cxcMZEh/sIA8l+kVaQXsolG4EDJgIiBTnJDva\nHaDtznAG4DoiqZjwDfKgoqyW9/9yit3vnKKiTEpT6Av0KuEATddxBMrzT2FodHyrwvY6wwHUAvGz\nokm6bzTuHi5kphaz/eUjHP48nUZ9253ipBrWvYNec47TjMrNz1yXoPh7SnO/JTR6tkP96az22oBA\nT+Kj/IgfFcqRLzM4ffw6xw5e4fszudyzcDjxo0IRBEGqYd3L6HUjDkDwIHMYTmnucUwmx/b47Kgz\n3M212VzdlCQuHsGajdOJjPGnpkrHh387y993fEtJkUaqYd3L6JXCUXtH4OETRWNDDZWF5x3tToed\n4VoTHOrFT9ZP4d4Hx+Hl7cr1q2W8ufUI+rwa2lp3k2pYOye9bqrWTPCg6WjPZ1OccxS/sDscWs6p\ns85wrREEgeFjBzBkWDDJX1/h+OEsQkTwB3IRaR3uKNWwdj56rXC8A4ehcg+gXltEavZRSgUVoZ5B\ndu+tYw3xUX5WHeAuKgWz58XjGerF+/84iw8C0QgEIVKESBUgItWwdkasmqodOnSIpUuXMm/ePJ5/\n/vnu8skiBEGGOmQcAKlpB/jzqb/xm6+38quvXqCstndNbSaMHYAp3ItMTOgQ8UBgMDJGIzDC05Uw\nH1dHu9hl+m0+Tm5uLs899xw7duzgk08+IS0tjaNHe76Qxs28fe0U9SaRaKWCAJn5rWR3oZeoI3lm\n1UQCwr25hEgWJjSIKBFw0+j54/Nfs/vdU1xNL0E0dX692hmXtPttPs7BgweZP38+QUHmjmmvvvqq\nVSVD7U1GWRZXK/M45+rCVDcXFqhV7NbWoxNp6SXqyGmbtbR1nuSnUnA25ToXzuaRebmYzMvF+Pq7\nc8eUSMZMGIi7h+oWG866pN2v83FycnJwcXHhkUceobS0lFmzZvHkk092p28dUqgxZ16e1OmJVsoJ\nUci5z8ONPdp6GkS6pZdoT9D6PGnevaNImD+Mi9/lceb4dUqKNBz8NI3DX2QwfHQYd0yNJDzSXOu6\nvSXt//rqIoJ39/zIjQryZv24jvsZ9et8HKPRSHJyMrt27UKtVrN+/Xr279/P4sWLO9zP0jYf1tLc\nK1QP7NHW8x8eboQq5CzzcONfmvpu6yXqCFSuCsZPHcQdUyLJza7gzIkcUi8UcOFsHhfO5hES5sWA\n2ECy22nLWFvfiNpL6fBC8o7Ox3FIm4+AgACmTJmCr68vYA6au3DhQqfCsabNhzUMDYghyncg2ZW5\n6ETYo9XxH56uDFDIedDXh8E+4TbZdWYEQSAi2p+IaH/maoZz7tQNzqbkUFRQQ1FBDaMRKAdKENHd\ntF/VxXJWDY8gYUKEVa9n6fK6Nf53Rp9r8zFr1iw2bdqERqPB3d2d5OTklnDvjrCmzYe1PDVtHS8f\n30F2ZS51oshujY6f+HgSiIGr595jyLiHkckddx7Wnag9VdyZMISpswZzNb2EIwevUJhTQTACwQjU\nIFKJSDXmHqfWLGl3x7lSv23zMWrUKNasWcP999+PwWBg6tSp3Htv513TrG3zYQ0Baj9euucZMsqy\nKNSUEOoZRJRHAJmnt6OtzOLquZ0MHru6R1qFOAqZTCB2WDCxw4LZ+PIh6ou0BABeCHg1pTUYFQI5\n5wtQNhiIjPFHqew4N6ij8J/X/nOmTX46Qz5On2/z0VX09ZVknN6OXleJV0AcMWMeQibrtdd6LaZ5\npLiWV4034I2Av1yG3HjzFEfGoMEBDI4PYnBcEH4Bt3a9S8uuYNO29rMrtzx+V6+MYrD3sdYnjyYX\nN19ix68j48x2asrSufb934kevbLPi6etJe24Qb6UFmu5mlbC1fRiblyr4Gp6CVfTzauSfgFqhsQH\nERMXxKAY/w7TJEAK/2mmzx5JKnc/YsevJfP0dqpLU8m+sIvoUQ8iyPp+CnPrJe2gEE+CQjyZOiuG\nBl0j2VfKzOJJK6GirJaTx7I5eSwbhVJGYJg3QUAV5hXL1kjhP2b6rHAAXN0DiB2/lozT26kquUT2\npd1EjVjRL8TTHipXJXEjQ4kbGYooipQUaZpGoxJysysozKkkEhmRQD0iWqAWkTogdICXNNo00aeF\nA+CqDiJ2/DoyT2+nsug8giBj0IjlDi9q6AwIgkBwqBfBoV5Mmz0YXb15NLr4fQGpFwpxM4EbtNRO\nkBXW8tYrRwkN9yF0oDeh4d4Eh3qh6GSxwRloNJgwGE12s9fnhQPmMrpDxq8h88ybVBR+hyCTEzns\nx5J4WuHqpiR+VCjxo8yj0cmzeVzLKkesb0RX00BxfnXLNaNzTfl1MplAYIgnoeHeZkGFexMS9oOY\n7H0tyFqaF0wqanSdP9kK+oVwANw9w4i941Eyz7xFef5pBEFORPxSh19Nd1YEQWDy+IFMHj+w5TGT\n0URpiZbC3GoK86oozKumqKCa4oIaigtqOH8q17yvTMAvQE2BVkdhnZ56zNeSIgZ48czqST0aN9fW\n0ro96DfCAXD3CmfIHY+SefYtyvK+RRDkDIxbJInHQmRyWcvUbsxEs6DaE1N5iRYVMOimAHxjvoZt\nvz/M8LggfP3d8QtQ4+uvxi9AjbePKzK5fWcAbVUgshf9SjgAau+BDBn3U66c/QuluccRZDLCY5Mk\n8dhIW2K6nFXO//w5GXdAjYAroGq6CQYTGZeKbrcjE/DxaxZT098ANX7+7vj4u9tU0LGzpfWu0O+E\nA+DhM4jB4x7m6tm/UJJzDEFQMGDIjyTx2ImiilrqgXqgnB8uvgqACyIrZsQQ6aemoryWirI6Kstr\nqSyro6Kstt26cy5uStSeKgID1Hh4qfDwdMXDS4Wnl+sP9z1VyBU/jFrduXTeL4UD4OkbzeBxD3Pl\nu3covn6YCl0NFZ6RDk+/7gu0d8CKmM91ho8Mu22hQDSJ1FTrqCivpbLMLKiiwmoys8oRGo3o6xvR\n1zdSWdLxKOLmrrxFTCM8XCnV6sgB7Lem1o+FA+DpN5jAuHspTN0DRWe5mvMtf6vXE+gTzlPT1hGg\nlq5Z2EJzyay2zi9uLpl1M4JMwNvXDW9fN6IGm0sbP/nqEbIazeW/5IgoASUQ7uvO4mlRaDUNaGsa\n0NTo0Nbo0GoaqK9rpL6ukZKmNpNuQAQyzmEWrb3o18IB2Jb6FYKmniUeroxVKRnlouBSQwnbj/+J\n39zzG0e712t5ZtXENiOs2yqZ1RatT+yNTTcdkFZZx+pBfkxtQ4CNegNaTQOamgZqNTo0NQ1oa3R8\ndkyOysV+15v6tXAyyrLIrjQvob5TU8ckVxdGuCgYrVIySqzm/KntDIlLQu3V93J7uhtrS2a1xtaY\nOaWLAl9/Bb7+twaveryiuu25XaFfC6c5/RqgyiTyZV0Dx+v1jHdVMkalRKi6Rvq3r+PpP4TQqNl4\n+MZ0uIBwc3qDdJ5kxtqSWc10dmLv6Ji5fi2cttKrtaLIkXo9KTo9m+JnYSpLRVN+BU35FdTeEYRE\nzcI7cNgtUQdltRUtCXXNRPkOlM6TuoAt50k9Sb+OOWlOv26LAT4RjBx5HyOnP8PAuEW4uPpQW32D\nrPN/JfXEK5Tnn0E0mbMWW4sGem+ZKmfCmtLCPU2/HnHg1vTrZqJ8B/LLaWsBkMldCIq4k8DwKVQU\nnaMo+zC62hKuX95DQdaXEDCc3FaiaaY3lqlyJrp6ntSd9HvhtJV+3daBLsjk+IeNxy90HNWlqRRl\nH6a2+gbkHWe9t5ozDXq+a2ikoVU+bW8tU+VM2Hqe1J30e+E0MzQgxqIDXBBk+ASNwDtwONrKa2Rl\nfIq7Jo/pbiomubqQpjeQZzBSaDBSYRL7VJkqiR9w+DmOsV5HTVq6o92wGkEQ8PSLYcyUJ/gKT9L1\njbgAY1RKFqhdedRbzX/6eiLLOUTB1S+pLk3HoO+8jWFGWRZHslPIKMvq/jchYTMOH3Eaa6q5+PSz\nqGOiid+8CVVgoKNdsppHpz3Oy8d3cLQ6j0iFnAEKOREqFd6YqCnPpKY8s+W5KvcA1N4RLTc3z1Bk\nMoW0MtfLcHiVG11xMVtihgKgjolmzB9edpQ7Xab1eZKhsZ66mlxqq3Korc5FW52DsVXfUkGmwN1r\nAN9VFZGurabAaKTmpsLqUb4DeemeZ3r6rfQ5HF7lZuPGjaSmpuLm5gbAz372s5aqi12lNusaNWnp\neMXH2cVeT9P6PEmhdMPLPxYv/1jAXJRPX1+OtvoGtdU3qK26Qb2mgNqqHIYCQz3M7TzqTCJakwmt\nKKJtKOHChfcJ8YlAqfJCqfLCReWFQuXZ56v2NGOPC8uNJgPGpssH9sDqT/7y5ct88MEHNpeubU2j\nUoVBoUBhMAfz6QoLe61wOkMQBFTuAajcA/APNff2MRkbSc74km8zvyJUISNMLsdHLsNdJqd5WaGx\n6Dtyi767zZ5CqUap8mwRVPOttKGO8sYGAjyDifGPQaZwQS5XIcgUvSp1wh7T12YblfVVdvXNKuFU\nVVVRUVGVniJQAAARK0lEQVTBpk2bKCgo4J577uHxxx/vkgMNru58Ne9+Er7Yg0rfgGtoaJfs9TZk\nciXBQfGcufhZS/iuElDLBDxkAh6CjGVxs/FVKGls0NDYUNNyMzTWYmispV57e2KYGnM+zKVbHhWQ\nyV2QK1TI5C7m/+U3/a9wQdZ0Xy53oUynoUZfh7ebNyGeIchkCgRBhiCTIwjyzv8KMq5V5VKkLSfE\nM4gh/k0hS4KAgACCrEMhd3Rh2dLpa1s27IFVwikvL2fatGn89re/xcXFhbVr17J3716WLl1qswMy\nk5GS0IF8kbSSxakn+uxo0xE3F5AHaMQcO1dlEon2DWfksCVt7mcy6ptEZBbURxc+okFXjYcgoJYJ\nuAgCSsBdocTHxR2ToQGT0XyzFFfMes7pwvvzBGqBdtscC7ImIf0gKBMiSaYGRG81pqZkuOYzP9FU\nzneH/welXAFNFXh+2LdZiAJ6YyMJxgpEL3f+RyagsaApl6VYJZyYmBhee+21lvsrV67k448/7lA4\nnbX5cKvT4ldaREVgCJ/evYzYugYC3O0bydob6CyCoS1kcpeWqV9GWRZfVZa288x6fpewhqEBMYii\nCZOxEZOxAaNRj8mox2Ro/r8Bk1HPR5c+pqauChdBQC6AHPN1Cx9XDyaGjUYUjYgm421/TS33TRRp\nCjEYDeZ9BfPhLQAKmRy10g1EERERRJO5ELtoQmz+v0kibrIfRNAasVGLvrHzz9W/qY6BDAFwQJsP\ngEuXLrU0lQJzSwa5vOMch87afPi4q9k4Poq/aeRkVmh56dtM/nPCYMI83axxrddjaQRDe9wc6d3e\n9qEBMQiCDLlChVyhoq1S9BllWRxsT4D1FYwed0enfmWUZbHj6/Z7ev4uYW2HNkRRJKPsKs8d+kOL\n4JpvYB5Unp3+cwb7DwKxeTQSm/RmarGRVXGdrclvAlBjMj/e420+ABobG3nxxReZOHEiSqWS3bt3\n8+Mf/7jDfSxp8xE8PJ4njCbePHeNCyU1bPk2kycmDCbKR93mfn0ZSyMYWtNZhIKlEQyWCrA7bQiC\nQFzgECJvmr7eTLRPBHEhwzp8DYD4sNEE+ISTXZnbMs07e/Zsu8+3ps2HVZEDY8eO5f7772fZsmUk\nJSUxcuRI5s2b1+E+Li4ueHh4tHm7ue2Di1zG+nExTArzpbbRyCsnr5BerrHGvX5NR5He0b4RFovR\nHgK0l4ifmrbutvfU2fS1MxvtHYseHh5WCcfhF0Dh1otSJlFkd2oeh3NKUcgE1o6NYkywj6Nc7FXY\nK/rgV1+90PYvvW8Ev79nc4/ZaMYe13Gmz5qB0WTk+NFkm/ZvjdMJB8zz04+vFPLp1SJkAjw0MpKp\n4f6OcLFX0tUDzZ7XT5wlhMjekQNOKZxmDmaXsCctD4Dl8eHMiZIijXsSe/zSO0s6ucNDbnqSOVFB\nuCvl/PViDnvS8qhtNLBwSGivuvrdm7F1ocLeNpwRh6cVdMbUcH/WjY1GIRP49GoR76fmYXLcICnR\nSzE1NmKst1/HAqcXDsDYEB82jB+MSi7jcE4p735/HYMdrwJL9F0aSks5v/Ep9BUVNNbYrwB7rxAO\nQHyAJ7+YNAQPpZyTBZVs/y4LfVOjoJq0dEoOHe6VCXF9HUd/N2kvbqE265rd7Tr1OU5ronzUPDU5\nlldPX+VCSQ2vnkhj5se7MGRmtDzH1oS4mrR0dIWFuIaG9st4OXvTUFp620Hb08mKNWnp3SIa6EUj\nTjNhnm48PSWWIHcVVzUN7BsxDZ3rD42KarOukfbiFovtNQ/lF59+liuvb+Pi089yfuNTNJS2F/fl\n/Dj6Vx7a/qW39rvpKrrCwm6z3atGnGb83VQ8FiDj9e+KKQ8M5d8LVzLs0hnc6rS41WlwLy6jIjUN\nv2Hxndrq6Avubdmo9v6Vt3UU7uiXvieTFbszRaVXCgdAWVJE4id/5+vE5ZSEDuTbuxJv2f7h9XrU\n+d/jrVKab67mvz4qJd4qBd6uSuR5uVTl5LYZ7GjLF+zo6Z69fgS6KsDOfumtTVa09XP1io9DHRPd\nLdM1p74A2hE1aelcfPpZDAoFV4eMQuPtQ727B/XuHtS5e6D3D6TewoYoCn0D7nVaFIZGBFFsuplw\nCwlB5eNjDosXhKbwePNfmSAgCOZwdVHfgCYtHZNGA01RugpPD7xHDEfu6tryOh1dfWq9zdprVfqq\nKirP3J4l2ozv+HG4+FgWulR+6jSGmtvjBBVenvhPnNBjvhh1OqouXLzFF4WXJz6jRt7yuVpi4587\ntuCqq+P4pUud72QBvXbEufnXJC7t1i9JHRPDmD9sQW80Ud3QaL7pGn/4v8FAdUMjFdUaKqo06NzU\n1Pi0EdIjApUWtsMLiYSQVo+V1wF1bT27e4gb3f42rRG05ZbZCYuGsHa25Vloozt9KavFnBpnIWHR\nGJQqlC5tzS1so9cKByB+86Z2phRPAeaI60B3FYEdJMad3/gUmmvZ6NzUGBTmjEKTTIZrRARDn/oF\nJkRE0Rx8ahLNcXQmmv6KItobuVx7dyeiIANBQGw1dkT85AHcw39oE3L78G5+pLG6hry9+2koLm55\njmtIMAOWLEZpQX2H+vx8cv7xz5seaeXHg/+B+4ABndqpvnSJws++aHd76PxEvEeM6NROY00NefsO\n0FBU3PKYKiSY8CWLLHo/dfn53PjH++1ut/T9NJOmUiC3Y3PeXi0cVWAgY/7wcpfOLZrFJ2stvg1r\nUXl1nkxXcqEEfW77c+ghNaUEhXR+oJ3f8gIhrefi2Rm4luRZdn4S6ovHPxvbnM+rY2IYM75zHwBq\nqgK5mN3+atzI8JV4hfpa5A9Pb7D5uylJO4+sAz+G1JQSZOF7AlDZuaN1rxZOM17xcTafjHdVfJ2t\n3FiysmOvVajORmBL6OiEWh0TY/XnbOt3Y4/PtTvpE8KxB7Z+wfY40Oy1CmWPERjsI8CuYm8B2xtJ\nOHagqweavX9duzICg/0E2FWcQcDtIQnHDnT1QHPWX9euCrCrOIuA20ISjh3pyoHmzL+ujsbRAm4L\nSThOgjP/ukrcjiQcJ8MZf10lbqfXRUdLSDgDNgnnpZdeYvNm60r8SEj0JawWTkpKCvv37+8OXyQk\neg1WCaeqqorXXnuNdevWdZc/EhK9AquE89///d9s3LjRbk2lJCR6Kxavqn3wwQeEhYUxadIk9u3b\nZ/ELdNTmo6SkBKPR2JKXIyHRXRQWFiKXy+3W5sPiRLaHH36Y0tJS5HI51dXV1NXVsXDhQp599tkO\n93vjjTc6bPMhCAJhYWFdKjIoiiIajQZPT0+b7Ug2nNcXe9ho/pE2GtvvA2pNmw9EG9i7d6/49NNP\nW/TchoYGUaPRtHkrLCwUY2NjRY1GY4sbLWg0mi7bkWw4ry/2tFFYWNju8djQ0GCxvW6/AGrN8Cch\n0d00t/ToKjYJZ8mSJSxZ0nZfSgmJ/oAUOSAhYQOScCQkbEASjoSEDUjCkZCwAflzzz33nEMdkMuZ\nNGlSp23fe8KOZMN5fXEWG804tJKnhERvRZqqSUjYgCQcCQkbkIQjIWEDknAkJGxAEo6EhA1IwpGQ\nsIEeF45WqyUpKYmCggIArly5wvLly5k3bx5PPPEEOl3nvejfe+89kpKSSEpK4plnnsFgMJCZmWmV\nna1btzJ//nySkpLYuXMngNU2mrm5eIktNjZu3EhiYmJL8OzBgwet/lwOHTrE0qVLmTdvHs8//7zV\nvuzevZvFixezZMkSFi9ezMSJE/nVr35ltR8HDhxgwYIFJCUlsWXLFps+k7feeovExEQWLVrEm2++\naZUNS48vrVbLY489xvz581m2bBk3btzo0KfbsDnBwQbOnTsnLliwQBwxYoSYn58viqIoLlq0SDx9\n+rQoiqL4+uuvi6+88kqHNr7//nsxKSlJ1Ol0oiiK4qZNm8T33nvPKjtHjhwRH3zwQdFkMok6nU6c\nPXu2eO3aNat9EUVRPHHihDh58uSW/CRbbNxzzz1idXX1LY9ZY+fGjRviXXfdJRYXF4sGg0F84IEH\nxCNHjtjkiyiKYnZ2tpiQkCAWFxdbZaO+vl6cOHGiWFFRIRqNRnHZsmXiiRMnrLJx4sQJccGCBaJW\nqxVNJpO4fv168auvvrLIhjXH1/PPPy9u27ZNFEVRTElJEZcvX27RZ9NMjwpn8+bN4pkzZ8TZs2eL\n+fn5YmFhoZiQkNCyvaCg4Jb7bXH9+vWWD0IURfGdd94Rn376aavtGAwGURRFMS8vT5w5c6ZNvlRW\nVor33XefuHPnTvHpp5+22cb48ePFtWvXiklJSeIbb7xhtZ13331X/P3vf99yv6SkRMzPz7fal2ZW\nr14tHjhwwGo/tFqtOGHCBLGgoEDU6XTi0qVLxVOnTlll4y9/+Yu4ZcuWlvu7du0SV65caZENS46v\nOXPmiKIoinPmzBELCgpats2ZM0csLCxs16/W9OhU7YUXXuCOO+5AbApWKC4uJjg4uGV7UFAQxcXF\n7e0OQGRkJOPHjwegvLycf/zjH0RHR1ttRy6X89prr7FgwQImT55MUVGR1TZaFy+x5f2Ul5czbdo0\ntmzZwr/+9S9Onz7NgQMHrLKTk5OD0WjkkUceYdGiRezatYuSkhKrfQE4c+YM5eXlLFy40Or3o1ar\n2bBhA4mJicyYMYMBAwagVCqtsjFs2DCSk5Oprq6moaGBQ4cOWWzDkuOrqKiozW2BgYEt2yzBoYsD\nJtPt3W1lMstcysvL46GHHmLFihUtQrLWzpNPPklKSgoFBQUcP37cKhs3Fy9pxpb3ExMTw2uvvYaX\nlxeurq6sXLmSEydOWGXHaDSSnJzMyy+/zAcffMDFixc5efKk1b4AvP/++6xevRqw/v1kZGSwd+9e\njhw5wrFjx5DJZCQnJ1tlY8qUKSxZsoSVK1fy6KOPMn78eAwGg03vpSP/u3LsgYNrR4eEhFBSUtJy\nv7S0lJCQ1h1obyctLY21a9eydu1aHnjgAQoLC62yc/XqVUwmE7Gxsbi6unLPPfdw+fJlSktLLbbx\n+eefU1paSkpKSkvxEplMZpUNgEuXLlFaWsqsWbOAH75Qa+wEBAQwZcoUfH3NLQYTEhJa7Frji8Fg\nICUlhd/97neA9d9PcnIykydPbvFjyZIlvPPOO1b5UVtby913382qVasA2LlzJ2FhYZw9e9aq99KZ\n/yEhIZSWlraMOjf/bwkOHXFCQ0Nxc3PjzJkzAHz44YfMmDGjw30qKip49NFH+a//+i8eeOABm+xk\nZWXx3HPPYTAY0Ov1HDx4kOnTp+Pq6mqxjXfffZdPPvmE/fv3s2HDBmbPns3zzz9vlQ2AxsZGXnzx\nRWpra9Hr9ezevZvly5dbZWfWrFkcP34cjUbTMvqMGzfOal8yMzOJjIzE3d0dsP5zjYuLIzk5mbq6\nOkRR5NChQ0yYMMEqP/Lz81m3bh0GgwGNRsOHH37IfffdZ/V76cz/mTNn8uGHHwJw8uRJ1Gq1VcJx\nyIhzc4mfP/zhD/z6179Gq9USHh7O1q1bO9z3r3/9K7W1tfzpT39i27ZtCILAjBkzrLIzd+5cUlNT\nWbRoEXK5nHnz5pGYmEhMTAy/+c1vLPalLV555RWrbIwdO5b777+fZcuWYTQaSUxMZN68eQwZMsRi\nO6NGjWLNmjXcf//9GAwGpk6dyr333svo0aOt8uXGjRuEhd3aG92az3XatGksXLiQpUuX4uLiwogR\nI1izZg2JiYkW+xEbG0tSUhKLFi3CZDKxevVqxo4da9Xnasnx9cQTT/Dss8+SlJSEi4sLL730Uoef\nzW2vIYpSWoGEhLVIkQMSEjYgCUdCwgYk4UhI2IAkHAkJG5CEIyFhA5JwJCRsQBKOhIQNSMKRkLAB\nSTgSEjYgCUdCwgYk4UhI2MD/A32yxN79U5HMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11188e6a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Rewrite the data to check\n",
    "t1=np.array([ 10,  20,  30,  40,  50,  60,  70,  80,  90, 100])\n",
    "cA1=np.array([ 8.649,  7.441,  7.141,  6.366,  6.215, 5.990, 5.852,  5.615,  5.481 ,  5.644])\n",
    "t2=np.array([ 10,  20,  30,  40,  50,  60,  70,  80,  90, 100])\n",
    "cA2=np.array([ 7.230,  6.073,  5.452,  5.317,  5.121, 4.998,  4.951,  4.978,  5.015,  5.036])\n",
    "t3=np.array([ 10,  20,  30,  40,  50,  60,  70,  80,  90, 100])\n",
    "cA3=np.array([ 5.137,  4.568,  4.548,  4.461,  4.382, 4.525,  4.483,  4.565,  4.459,  4.635])\n",
    "\n",
    "# Plot the data\n",
    "plt.plot(t1, cA1, 'o', label='data')\n",
    "plt.plot(t2, cA2, 'o', label='data')\n",
    "plt.plot(t3, cA3, 'o', label='data')\n",
    "# Plot the fitted curve\n",
    "plt.plot(t1, function(t1,10.,298.15,logA,E,dH,dS ), label='fitted')\n",
    "plt.plot(t2, function(t2,10.,308.15,logA,E,dH,dS ), label='fitted')\n",
    "plt.plot(t3, function(t3,10.,323.15,logA,E,dH,dS ), label='fitted')\n",
    "# Add the legend, in the \"best\" location to avoid hiding the data\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fitted curves matched all the data sets!  Time to calculate the errors..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Used help from http://stackoverflow.com/questions/14581358/getting-standard-errors-on-fitted-parameters-using-the-optimize-leastsq-method-i\n",
    "\n",
    "In order to solve for the errors we first need to convert the covariance matrix by multiplying it by the residual variance.  We then need to get the square root of the diagonals of the matrix for the erros associated with the fitted parameters.  \n",
    "\n",
    "This code below is written do this with the previously defined leastsq_function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.367232258793\n",
      "2137.71097888\n",
      "790.492861472\n",
      "2.51289798174\n",
      "\n",
      "Fit paramters and parameter errors from leastsq method :\n",
      "pfit =  [  9.48853174e+00   6.43077887e+04  -1.15315215e+04  -3.73495283e+01]\n",
      "perr =  [0.36723225879332022, 2137.7109788831367, 790.49286147190332, 2.5128979817355352]\n"
     ]
    }
   ],
   "source": [
    "args = (t, cA, (10, 10, 10), (298.15, 308.15, 323.15), (n1, n1+n2))\n",
    "\n",
    "s_sq = (leastsq_function(pfit,*args)**2).sum()/(len(cA)-len(params0))\n",
    "pcov = pcov * s_sq\n",
    "error = [] \n",
    "for i in range(len(pfit)):\n",
    "    error.append(np.absolute(pcov[i][i])**0.5)\n",
    "    print(np.absolute(pcov[i][i])**0.5)\n",
    "pfit_leastsq = pfit\n",
    "perr_leastsq = error\n",
    "\n",
    "print(\"\\nFit paramters and parameter errors from leastsq method :\")\n",
    "print(\"pfit = \", pfit_leastsq)\n",
    "print(\"perr = \", perr_leastsq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "starting_guess=params0\n",
    "optimized_parameters=pfit_leastsq\n",
    "standard_errors=perr_leastsq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6, 45000.0, -10000.0, -50]\n",
      "[  9.48853174e+00   6.43077887e+04  -1.15315215e+04  -3.73495283e+01]\n",
      "[0.36723225879332022, 2137.7109788831367, 790.49286147190332, 2.5128979817355352]\n"
     ]
    }
   ],
   "source": [
    "# Finish your notebook with this cell\n",
    "print(params0)\n",
    "print(pfit_leastsq)\n",
    "print(perr_leastsq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [py3]",
   "language": "python",
   "name": "Python [py3]"
  },
  "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
}