Practice_Final_Solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practice Final Exam, ChE 263\n",
    "\n",
    "This exam covers the following topics:\n",
    "\n",
    "* Problem 1: Variable types and formats\n",
    "* Problem 2: Debugging\n",
    "* Problem 3: Solve Linear Equations\n",
    "* Problem 4: Differentiate and Integrate\n",
    "* Problem 5: Interpolation\n",
    "* Problem 6: General Nonlinear Regression and Plotting\n",
    "* Problem 7: Equation Solution and Plotting\n",
    "* Problem 8: Parameter Estimation\n",
    "* Problem 9: Implicit Equation Solve\n",
    "\n",
    "A template for each problem is provided with instructions on what is needed. In some cases, the problem is to debug a section of code to produce a specific result. In other cases, the problem is to import a module and produce a result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1: Variable types and formats\n",
    "\n",
    "#### Part A\n",
    "Add **x=2.0** (float) and **y='3'** (string) together to produce a number 5 (integer). Show the type of variable with **type()** to confirm that the result is an integer.\n",
    "\n",
    "Hint: A string is converted to a number with **int()** or **float()**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5\n",
      "<class 'int'>\n"
     ]
    }
   ],
   "source": [
    "x = 2.0\n",
    "y = '3'\n",
    "\n",
    "# add two integers together (5)\n",
    "z = int(x) + int(y)\n",
    "# show the result of adding\n",
    "print(z)\n",
    "# show that z is an integer\n",
    "print(type(z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part B\n",
    "\n",
    "Show $\\pi$ to 30 decimal places."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pi = 3.141592653589793115997963468544\n"
     ]
    }
   ],
   "source": [
    "print('pi = {:.30f}'.format(np.pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2: Debugging\n",
    "\n",
    "You are working on the launch sequence of a rocket. The following code is intended to count backwards from 5 to 0 with steps of -1 in 1 second intervals. Print `blast off` after reaching zero.\n",
    "\n",
    "```python\n",
    "import time\n",
    "for i in range(5,1)\n",
    "    time.wait(1.0)\n",
    "    print(i)\n",
    "print('Blast off\")\n",
    "```\n",
    "\n",
    "The code has a few bugs (errors) that prevent it from running or producing the correct result. Find the errors in the code to produce:\n",
    "\n",
    "```\n",
    "5\n",
    "4\n",
    "3\n",
    "2\n",
    "1\n",
    "0\n",
    "Blast off\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5\n",
      "4\n",
      "3\n",
      "2\n",
      "Blast off\n"
     ]
    }
   ],
   "source": [
    "# corrected code\n",
    "import time\n",
    "for i in range(5,1,-1):\n",
    "    time.sleep(1.0)\n",
    "    print(i)\n",
    "print('Blast off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3: Solve Linear Equations\n",
    "\n",
    "Find the values of $x_0$, $x_1$ that satisfy the following equations:\n",
    "\n",
    "$5 x_0 + x_1 = 2$\n",
    "\n",
    "$3 x_0 + 12 x_1 = 1$\n",
    "\n",
    "Put equations into matrix form with $A \\;x = b$ and solve $x = A^{-1}\\; b$. The $A$ matrix is given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.40350877 -0.01754386]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[5,1],[3,12]])\n",
    "b = np.array([2,1])\n",
    "\n",
    "x = np.linalg.solve(A,b)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4: Differentiate and Integrate\n",
    "\n",
    "#### Part A\n",
    "\n",
    "Compute the derivative $\\frac{d f(x)}{dx}$ for the following function:\n",
    "\n",
    "$f(x) = \\frac{1}{1+e^{-x}}$\n",
    "\n",
    "Find a solution with scipy (numerically) and sympy (symbolically). For the numeric solution (scipy), use value of $x=0.2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.229249632313315\n"
     ]
    }
   ],
   "source": [
    "# numeric method\n",
    "from scipy.misc import derivative\n",
    "def f(x):\n",
    "    return 1.0/(1.0+np.exp(-x))\n",
    "print(derivative(f,0.2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{e^{- x}}{\\left(1 + e^{- x}\\right)^{2}}$"
      ],
      "text/plain": [
       "exp(-x)/(1 + exp(-x))**2"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# analytic solution with sympy\n",
    "import sympy as sym\n",
    "x = sym.Symbol('x')\n",
    "dx = sym.diff(1/(1+sym.exp(-x)),x)\n",
    "dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part B\n",
    "\n",
    "Compute the integral for the function:\n",
    "\n",
    "$f(x) = \\frac{1}{1+e^{-x}}$\n",
    "\n",
    "Find a solution with scipy (numerically) and sympy (symbolically). For the numeric solution, use limits of integration for $x$ between $-1$ and $1$:\n",
    "\n",
    "$\\int_{-1}^{1} \\left(\\frac{1}{1+e^{-x}}\\right) dx$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x + \\log{\\left(1 + e^{- x} \\right)}$"
      ],
      "text/plain": [
       "x + log(1 + exp(-x))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution 1: Numeric\n",
    "from scipy.integrate import quad\n",
    "def f(x):\n",
    "    return (1.0/(1.0+np.exp(-x)))\n",
    "result = quad(f,-1,1)[0]\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x + \\log{\\left(1 + e^{- x} \\right)}$"
      ],
      "text/plain": [
       "x + log(1 + exp(-x))"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution 2: Symbolic\n",
    "import sympy as sym\n",
    "x = sym.Symbol('x')\n",
    "fs = 1/(1+sym.exp(-x))\n",
    "result = sym.integrate(fs)\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 5: Interpolation\n",
    "\n",
    "Data for `x` and `y` is shown below:\n",
    "\n",
    "```python\n",
    "x = [1,3,7,10]\n",
    "y = [2,9,20,35]\n",
    "```\n",
    "\n",
    "Create a cubic spline interpolation to approximate the value of `y` at `x=5`. Show a cubic spline interpolation on a plot together with the data points. Label the plot with appropriate $x$ and $y$ axis labels, and a legend. Save the plot as `plot.png`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14.20634920634921\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAna0lEQVR4nO3deXxU1f3/8deHEMQQBAWNKEpYXBB+BgyLisrihqgoKCJSBFvF1rq0WhSKX7VfBLXwVXFpLYqFltQAFoFaS2shoBVRiAZks0ghCFKRnRgDAc7vjzMsgSAhyZ07k3k/H495ZO65M3M/OeInZz733HPNOYeIiCSOamEHICIi0aXELyKSYJT4RUQSjBK/iEiCUeIXEUkw1cMOoCzq16/v0tPTww6jQr799ltq1aoVdhgxQ/1xgPqiJPVHSRXpj9zc3I3OuZMPbY+LxJ+ens6CBQvCDqNCZs+eTadOncIOI2aoPw5QX5Sk/iipIv1hZvmltavUIyKSYJT4RUQSjBK/iEiCiYsaf2mKi4tZu3YtRUVFYYdSJnXq1GHZsmVhh0HNmjVp2LAhycnJYYciIiGJ28S/du1aateuTXp6OmYWdjhHtWPHDmrXrh1qDM45Nm3axNq1a2ncuHGosYhIeAIr9ZhZTTP72MwWmtkSM/tVpH2cma0ys7zIo1V5Pr+oqIh69erFRdKPFWZGvXr14uZbkkgiy8qC9HTo0qUj6el+u7IEOeLfCXRxzhWYWTLwLzP7W2TfIOfcmxU9gJL+sVOficS+rCwYOBAKCwGM/Hy/DdC3b8U/P7ARv/MKIpvJkYfWgBYROYqhQ2Fv4Xcks2t/W2Ghb68MFuR6/GaWBOQCzYCXnXOPmNk44CL8N4KZwGDn3M5S3jsQGAiQlpaWmZ2dXWJ/nTp1aNasWWCxH6sRI0aQmprK/fffX+r+6dOnc/bZZ3PuuedGObLDffHFF2zbti3UGAoKCkhNTQ01hlihvihJ/eHLO5e7f/IIz9CVGeyJFGfMHLNmzSnz53Tu3DnXOdfmsB3OucAfQF0gB2gJNAAMOA4YDzx2tPdnZma6Qy1duvSwtu8zYYJzjRo5Z+Z/TphwTG8/qscff9yNHDnyiPtvu+02N3ny5Mo9aDkda98FIScnJ+wQYob6oiT1h89R4NzpfOnA7X80anRsnwMscKXk1KjM43fObY0k/q7OufWRmHYCvwfaBX38ffWy/HzfffvqZRU9WTJ8+HDOPvtsLrnkEj7//HMAXn31Vdq2bUtGRgY33XQThYWFzJ07l3feeYdBgwbRqlUrVq5cWerrRESYMYPxPaaSkgLraLi/OSUFhg+vnEMEOavnZDOrG3l+PHAlsNzMGkTaDLgRWBxUDPsMHbrvJMkBFa2X5ebmkp2dTV5eHu+88w7z588HoGfPnsyfP5+FCxfSvHlzxo4dy8UXX0y3bt0YOXIkeXl5NG3atNTXiUiC+/xzuPVWOr43jFdf2UOjRr6806gRjBlTOSd2IdhZPQ2A8ZE6fzVgknPubTObZWYn48s9ecCPA4wBgDVrjq29LN5//3169OhBSkoKAN27dwdg8eLFPProo2zdupWCggKuvvrqUt9f1teJSILYvBmuuw6OOw6mTOG2Rknc1g9mz55T6YvWBZb4nXOLgNaltHcJ6phHcuaZvrxTWntlGzBgAFOnTiUjI4Nx48Yxe/bsCr1ORBJAcTH06uVHozk50KhRoIdLiLV6hg/39bGDVbRedtlllzF16lS+++47duzYwV/+8hfAX6HboEEDiouLyTroJEJqaio7duzYv32k14lIApoyBWbNgldfhYsvDvxwCZH4+/b19TFfL6NS6mUXXHABvXv3JiMjg2uuuYa2bdsCMGzYMNq3b0+HDh1KTN28+eabGTlyJK1bt2blypVHfJ2IJKDevWHuXLj99qgcLm7X6jlWfftW3omRfYYOHcrQUs4Q/+QnPzms7cILL2Tp0qUlXlPa60Qkgbz7Lpx6Kvy//wcXXRS1wybEiF9EJOYsWgQ9e8J99/l55lGkxC8iEm3//a+fwXPCCf6CoiivoZUwpR4RkZjw7bc+6W/aBO+/D6efHvUQlPhFRKJp1Cj49FOYOhUuuCCUEJT4RUSiacgQ6NABrrgitBBU4xcRiYbsbPjmG6hRI9SkD0r8UTF79mx69epV6r5u3bqxdevW6AYkItE1cSL06QMjRoQdCaBST+jeeeedsEMQkSDNnu0vzLrkEnjqqbCjATTir5A//OEPnH/++WRkZNCvXz8GDBjAm28euKPkwTeT2LFjB9deey3nnHMOP/7xj9m7dy8A6enpbNy4sdTPE5E4t3gx3HgjNG0K06ZBzZphRwRUpRF/aavX3XIL3HOPX4O5W7fD9w8Y4B8bN8LNN5fcd5RF05YsWcKTTz7J3LlzqV+/Pps3b+bBBx884utzc3NZunQpjRo1omvXrkyZMoWbDzpmaZ8nInHuvvv8wmB/+xucdFLY0eynEX85zZo1i169elG/fn0ATjrKf9TMzEyaNGlCUlISffr04V//+leFPk9E4sDEiX5ZhoBX2zxWVWfE/30j9JSU799fv/5RR/hlUb169f0lnL1797Jr14EbJdshV+Ydui0iVcS338Lzz8PDD8Mpp/hHjNGIv5y6dOnC5MmT2bRpEwCbN28mPT2d3NxcwN9cvbi4eP/rc3NzWbVqFXv37mXixIlccsklR/08EYkzu3b5dfUfewzmzQs7miOqOiP+KGvRogVDhw6lY8eOJCUl0bp1a5555hluuOEGMjIy6Nq1K7Vq1dr/+gsuuIB7772XL774gs6dO9OjR4+jft64ceOi/FuJSLnt2QP9+/t6/pgxcOmlYUd0REr8FdC/f3/69+9fom3eQX/ln3nmGQA6derEjBkzqF279mGfsXr16u/9PBGJA875E7nZ2fD003DXXWFH9L1U6hERqahVq2DCBBg0CB55JOxojkojfhGRimrSBPLyoHHjsCMpE434RUTK67e/heee88+bNIn6uvrlFVjiN7OaZvaxmS00syVm9qtIe2Mz+8jMvjCziWZWI6gYREQC88c/+gtEZ82CyDTueBHkiH8n0MU5lwG0Arqa2YXAM8BzzrlmwBbgRwHGICJS+SZP9lf9d+nin1eLr+JJYNE6ryCymRx5OKALsG9Bm/HAjUHFICJS6aZPh9tug4sv9s9jZP2dYxHoyV0zSwJygWbAy8BKYKtzbnfkJWuBUu87ZmYDgYEAaWlpzD7kyto6deqwY8eOYAIPwJ49e0qNt2XLlsyZM4d69eod82e+/fbbNGvWjHPPPfeY3ldUVHRYf0ZbQUFB6DHECvVFSbHeH6fNmkXa2WezaPBg9syfH/jxguiPQBO/c24P0MrM6gJvAWXOUM65McAYgDZt2rhOhyzCtmzZslLnxR/RqixYOBQK10DKmZAxHBr3Lfv7K2jHjh2lxmtmpKamHtvvEvH3v/+d5ORk2rZte0zvq1mzJq1btz7m41Wm2bNnc+h/00SlvigpZvujoABSU/2CkLt3c2n16EyKDKI/olKYcs5tBXKAi4C6ZravxxoC6wIPYFUWfDwQCvMB539+PNC3V8CECRNo164drVq14u677+bll19m0KBB+/ePGzeOe++9F4A+ffqQmZlJixYtGDNmzGGftXr1alq2bLl/e9SoUTzxxBMAvPrqq7Rt25aMjAxuuukmCgsLmTt3LtOnT2fQoEG0atWKlStXsnLlSrp27UpmZiaXXnopy5cvr9DvJyIRM2f6qZpz5/rtKCX9oAQ5q+fkyEgfMzseuBJYhv8DsG894v7AtKBi2G/hUNhTWLJtT6FvL6dly5YxceJEPvjgA/Ly8khKSiI1NZW33npr/2smTpzIrbfeCsDLL79Mbm4uCxYs4IUXXti/Jk9Z9OzZk/nz57Nw4UKaN2/O2LFjufjii+nevTsjR44kLy+Ppk2bMnDgQF588UVyc3MZNWoU99xzT7l/PxGJyMmB66+HU0+Fs84KO5pKEeSfrQbA+EidvxowyTn3tpktBbLN7EngU2BsgDF4hWuOrb0MZs6cSW5u7v4yy3fffccpp5xCkyZNmDdvHmeddRbLly+nQ4cOALzyyiv777b15ZdfsmLFijLX9RcvXsyjjz7K1q1bKSgo4Oqrrz7sNQUFBcydO7fELR537txZ7t9PRPBJ/7rr/Bz9mTPh5JPDjqhSBJb4nXOLgMMKyc65/wDtgjpuqVLOjJR5SmkvJ+cc/fv356lDbqX2+uuvM2nSJM4991x69OiBmTF79mxmz57Nhx9+SEpKCp06daKoqKjE+w5e0hkosX/AgAFMnTqVjIwMxo0bV+qJnr1791K3bl3y8vLK/TuJyEEWLYJrr/VJf9asmFxeubzia/JpeWUMh6SUkm1JKb69nC6//HLefPNNNmzYAPhllPPz8+nRowfTpk3jjTfe2F/m2bZtG3Xr1iUlJYXly5eXWMhtn7S0NDZs2MCmTZvYuXMnb7/99v59O3bsoEGDBhQXF5OVdeC8RO3atffPFDrhhBNo3LgxkydPBvwfpoULF5b79xNJeC1awEMP+VF/FUr6kCiJv3FfaDcGUhoB5n+2G1OhWT3nnXceTz75JFdddRXnn38+V155JevXr+fEE0+kefPm5Ofn066d/2LTtWtXdu/eTfPmzRk8eDAXXnjhYZ+XnJzMY489Rrt27bjyyitLTNEcNmwY7du3p0OHDiXab731VkaOHEnr1q1ZuXIlWVlZjB07loyMDFq0aMG0acGfPhGpct59F9atg6QkGDasypR3SnDOxfwjMzPTHWrp0qWHtcWy7du3hx3CfrHQdzk5OWGHEDPUFyWF2h9vveVccrJzvXuHF8MhKtIfwAJXSk5NjBG/iMjRTJwIN98MmZnwyithRxMoJX4Rkd//3i/D0KED/OMfULdu2BEFKq4Tv/8mI8dCfSZyiOJieOEFuOIKf9vEclxFH2/i9vKzmjVrsmnTJurVq4fFyRrYYXPOsWnTJmrG4aJSIpXOOdi9G5KT/Qnd2rXhuOPCjioq4jbxN2zYkLVr1/LNN9+EHUqZFBUVxUTCrVmzJg0bNgw7DJFwOQe/+AX8+98wZQrUrx92RFEVt4k/OTmZxnFymzPwCy2FvTCaiOBH+XfdBePGwb33+mmbCSaua/wiIsekqAh69fJJ/4knfG0/zm6iUhnidsQvInLM+vWDqVN9wr/vvrCjCY0Sv4gkjocfhptugshyKokq8b7jiEhi+fxzeO45/7xt24RP+qARv4hUZR98AN27+xO4P/hB1Vx3pxw04heRqumtt/xFWfXqwYcfKukfRIlfRKqe3/zG1/JbtfK3S2zaNOyIYooSv4hUPSecADfc4O+alWAXZ5WFEr+IVA2FhfDee/75D37gr8hNSfn+9yQoJX4RiX9ffQWXXQbXXANff+3btIbXEWlWj4jEt7w8uP562LLFr6mflhZ2RDFPI34RiV/Tp8Mll/jnH3wA110XbjxxIrDEb2ZnmFmOmS01syVm9kCk/QkzW2dmeZFHt6BiEJEq7oMPoHlz+PhjyMgIO5q4EWSpZzfwkHPuEzOrDeSa2buRfc8550YFeGwRqaqKiiA/H845B0aMgF274Pjjw44qrgQ24nfOrXfOfRJ5vgNYBpwe1PFEJAF89RV07AhdusC33/orcpX0j5lF41Z8ZpYOvAe0BB4EBgDbgQX4bwVbSnnPQGAgQFpaWmZ2dnbgcQapoKCA1NTUsMOIGeqPA9QXJR2pP05YsoQWjz1G0nffsXzIEDZeemkI0UVfRf59dO7cOdc51+awHc65QB9AKpAL9IxspwFJ+G8bw4HXj/YZmZmZLt7l5OSEHUJMUX8coL4oqdT+eO0152rUcK5pU+cWL456TGGqyL8PYIErJacGOqvHzJKBPwNZzrkpkT80Xzvn9jjn9gKvAu2CjEFE4pxz/mKsjh39SdwWLcKOKO4FdnLX/B3QxwLLnHPPHtTewDm3PrLZA1gcVAwiEse++gr27oWGDf38/Jo1obouPaoMQfZiB6Af8JmZ5UXafgn0MbNWgANWA3cHGIOIxKM5c6B3bz+6nzkTdA6kUgWW+J1z/wJKu2b6naCOKSJxzjnOyM6G116DZs38LRKl0unKXRGJDdu3w0030fR3v4MePVTPD5ASv4jEBjNYsYIvfvITmDTJL60sgVDiF5HwOAdvvOGXVK5dGz75hLW33KKVNQOmxC8i4di+Hfr0gdtug9/9zrclJ4cbU4LQ3CgRib7cXD9rZ/VqeOopeOCBsCNKKEr8IhJd2dlw++1+3fw5c6BDh7AjSjgq9YhIdLVtC716+RuoKOmHQolfRII3Ywb85Cf+ZG7TppCVBfXqhR1VwlLiF5HgFBXBz37m74X7/vuweXPYEQlK/CISlEWLoF07GD0a7r8f5s/XKD9G6OSuiFS+4mJ/A/SdO+Gvf4VuusNqLFHiF5HK8+WX0KCBn48/aRI0aQInnxx2VHIIlXpEpOKcg7Fj/do6//d/vq19eyX9GKXELyIV89VXvqxz553Qpg3cemvYEclRKPGLSPlNn35gzfzRo+Gf/4RGjcKOSo5CNX4RKb+0NDj/fL9+/llnhR2NlJFG/CJSds7B66/D4MF+u317mD1bST/OKPGLSNmsXg1XXw0/+hHMmwe7dvl2LaEcd5T4ReT77dkDzz/va/kffggvvwyzZkGNGmFHJuWkGr+IfL9162DoUOjUCX77WzjzzLAjkgrSiF9EDldY6E/YOucT/aefwttvK+lXEYElfjM7w8xyzGypmS0xswci7SeZ2btmtiLy88SgYhCRo8vKgvR0qFbN/5z18Axo2RLuusvf8Bzg7LNVy69Cghzx7wYecs6dB1wI/NTMzgMGAzOdc2cBMyPbIhKCrCwYOBDy86GBW8ev82+hy8hr2FZUw8/Wad8+7BAlAIElfufceufcJ5HnO4BlwOnADcD4yMvGAzcGFYOIfL+hQ31Vx9jLP7mC6/kLQ3mStskLoWPHsMOTgJhzLviDmKUD7wEtgTXOubqRdgO27Ns+5D0DgYEAaWlpmdnZ2YHHGaSCggJSU1PDDiNmqD8OCLMvhnQ+jlwyKaYGlzGHtTTkPzTFzDFr1pxQYtK/jZIq0h+dO3fOdc61OWyHcy7QB5AK5AI9I9tbD9m/5WifkZmZ6eJdTk5O2CHEFPXHAaH0xfr1zt1+u3Pg7mO082dxDzwaNYp+SPvo30ZJFekPYIErJacGOqvHzJKBPwNZzrkpkeavzaxBZH8DYEOQMYjIQYqL4bnn/MnaN95g8fVDeOP4H5V4SUoKDB8eUnwSFUHO6jFgLLDMOffsQbumA/0jz/sD04KKQUQOMWAAPPggXHIJLF5My+kjeP7VWjRq5CftNGoEY8ZA375hBypBCvICrg5AP+AzM8uLtP0SeBqYZGY/AvKBWwKMQURWroS6df1tD3/2M+jd2y+jHJme2bevEn2iCSzxO+f+BRxp4u/lQR1XRCK2bfM1m9Gj4e674YUXoG3bsKOSGKAlG0Sqmt27/Qqajz4KGzf68s6QIWFHJTFESzaIVDWDB/sR/rnnwvz5/o9AgwZhRyUx5KgjfjO7D5jgnNsShXhEpDw++wxq1vTr4t9zD1x0EfTsqWUWpFRlGfGnAfPNbJKZdY3M1hGRWLBunV8fv1Ur+J//8W1NmsBNNynpyxEdNfE75x4FzsJPzRwArDCzEWbWNODYRORItm3z6y2cdRZMmAA//zn85jdhRyVxokw1/sgVYP+NPHYDJwJvmtmvA4xNRI7k2WdhxAi48UZYvhxGjYKTTgo7KokTZanxPwDcDmwEXgMGOeeKzawasAJ4ONgQRYQ9e/zI/owzoEsXP8K/8UZo3TrsyCQOlWXEfxJ+nZ2rnXOTnXPFAM65vcB1gUYnkuicg6lTISPDT8scH1nYtm5dJX0pt7LU+B93zuUfYd+yyg9JRAB47z2/Hn6PHn5u/uTJMG5c2FFJFaALuERijXN+Rs6yZfD1134efr9+UF3/u0rl0AVcIrFiwQLo2tXf0Bzghz+Ef/8b7rhDSV8qlRK/SMhSV6yA7t39OjoLFkByst+RnAzHHRducFIlaRghEqZHH6XN8OFw4ol+QbX77oPatcOOSqo4JX6RaMvL89My69WDjh1ZvXYt6aNHQ506YUcmCUKlHpFoyc2FG27w0zCfe863XXklqwcMUNKXqFLiFwnaRx/BdddBmzZ+iuavfgW/+EXYUUkCU6lHJGi//jXMmwfDhvkavkb3EjKN+EUqk3MwYwZ07OjX0AF48UVYvdrfGEVJX2KAEr9IZdi7F6ZM8VMyr7kG/vMfWLvW7zvtNEhNDTc+kYOo1CNSUXv2+KUVcnOhaVMYMwb694caNcKOTKRUGvGLlEdhIWRn+9JOUhL06QNvvOHLO3fdpaQvMS2wxG9mr5vZBjNbfFDbE2a2zszyIo9uQR1fJBBbtsCTT0J6uk/2ubm+/aGH4NZbtbSCxIUgR/zjgK6ltD/nnGsVebwT4PFFKs/WrT65n3GGv8Vh27YwZ46foikSZwJL/M6594DNQX2+SFQUFPifNWr4Us6NN8LChfDXv8Jll4Uamkh5mb+rYkAfbpYOvO2caxnZfgJ/397twALgIefcliO8dyAwECAtLS0zOzs7sDijoaCggFTN7NgvpvvDOeosXMiZEydy/Jdf8vH48ZCURLWiIvbWrFnph4vpvgiB+qOkivRH586dc51zh38tdc4F9gDSgcUHbacBSfhvGsOB18vyOZmZmS7e5eTkhB1CTInJ/ti927nJk51r1845cK5+fed+9SvnCgsDPWxM9kWI1B8lVaQ/gAWulJwa1Vk9zrmvnXN7nL9t46tAu2geXwSAVVkwNR3+VM3/XJXl2995B3r1go0b4Te/gfx8eOwxOP74MKMVqXRRnYJgZg2cc+sjmz2Axd/3epFKtyoLPh4Iewr9dmE+fHCHf97tVpg2Da691k/RFKmiAkv8ZvYG0Amob2ZrgceBTmbWCnDAauDuoI4vUqqFQw8k/X2qFfv2xn39DVFEqrjAEr9zrk8pzWODOp5ImXybD1ZKe+GaqIciEhZduStVW3Ex/OlPsHKl304+tfTXpZwZvZhEQqbEL1XT1q0wciQ0aQJ9+8L48b697ShISin52qQUyBge9RBFwqLry6XqGTIEXnrJX3zVpQu88opfMRN8HR98Tb9wjR/pZww/0C6SAJT4pWpYuBAyMvzz7dv9FbYPPuhvc3ioxn2V6CWhKfFL/Nq9G956C5591t/hat48vzzySy+BlXYGV0RANX6JR4WF/mblzZrBLbfAN9/4ZN+ihd+vpC/yvTTil/ixcyccd5xfA3/YMGjZEp5/Hq6/XhdciRwDJX6JfR995Ms5S5bAokVQqxYsXQqnHmFqpoh8L5V6JDbt2ePvYXvJJXDhhfD3v0O3blBU5Pcr6YuUm0b8EpumTYObboLGjWH0aLjjDqhdO+yoRKoEJX6JDevWwYsvwmmnwf33+zVzpk6F665T/V6kkqnUI+FauBBuv93fw3bkSPj8c99evTrccIOSvkgAlPglPP/zP9Cqla/l//Sn8MUX8PLLYUclUuWp1CPRs3Onv29tp05+++qrITUVBg6EE08MNTSRRKIRvwRvyxZ4+ml/ovaOO/xqmeBn7DzyiJK+SJQp8UuwBg+GM8/0C6e1bOmnZQ4ZEnZUIglNiV8q3xdfHHi+aZNfMO3TT+Ef/4CrrtKSCiIhU41fKodzkJMDzzzjE/yCBZCZCWPGKNGLxBiN+KVi9uyByZOhbVu4/HI/PXPECGja1O9X0heJORrxS8UUFMCdd0Jamh/d9+sHNWuGHZWIfA8lfjk227f7O1r985/+RG2dOjB3Lpx7ri62EokTgZV6zOx1M9tgZosPajvJzN41sxWRn5rHFy82bIChQ/0MnUce8TX9LVv8vhYtlPRF4kiQNf5xQNdD2gYDM51zZwEzI9sS6z7+2C+p8NRTcMUVMH8+vPsunHRS2JGJSDkElvidc+8Bmw9pvgEYH3k+HrgxqONLBa1Y4ZM7+PvW3n23XwP/zTehTZtwYxORCjHnXHAfbpYOvO2caxnZ3uqcqxt5bsCWfdulvHcgMBAgLS0tMzs7O7A4o6GgoIDU1NSwwziqlFWraDRhAqfMnk3Rqafy0R//CNUqf3wQL/0RDeqLktQfJVWkPzp37pzrnDt8pOacC+wBpAOLD9reesj+LWX5nMzMTBfvcnJywg7h+y1Z4lzPns6Bc7VqOTdokHPr1wd2uJjvjyhSX5Sk/iipIv0BLHCl5NRoz+r52swaOOfWm1kDYEOUjy+H2rvXj+jXrvUzdR59FH72M6hXL+zIRCQg0b6AazrQP/K8PzAtysdPWFlZ/vxstWr+54z//RiuvdavpQNw5ZXw5Zf+JuZK+iJVWpDTOd8APgTOMbO1ZvYj4GngSjNbAVwR2ZaAZWX5lY/z86GN+5iX86+l6+Pt2fnePGjQwL/IDE44IdxARSQqAiv1OOf6HGHX5UEdU0o3dCgUFsLDPMMzDGYTJzGEEUw/8V6W/Fz3sRVJNLpyt6pbuJDq+alAU97mOpLYw4vcRwG1sbVhByciYdAibVXVsmVwyy3QqhVP1xoGwFJa8BS/pAA/yj/zzDADFJGwKPFXNatWQf/+/qYnf/sbPPoo7tnnSEkp+bKUFBg+PJwQRSRcKvVUNS++CJMm+SmZgwfDySfTC9hVy9f616zxI/3hw6Fv37CDFZEwKPHHu61bYeRIPx2zUyef3R96CE4/vcTL+vZVohcRT4k/Xn33Hbz0kl84bcsWvwZ+p06agy8iR6XEH48mTfKj+rVroWtXn/xbtQo7KhGJE0r88WLfYnpmsH69v/DqD3+Azp3DjUtE4o5m9cSDTz/1NfzxkRWtf/pT+OgjJX0RKRcl/li2bh0MGACZmZCXd+DG5dWr6ybmIlJuKvXEqldegQcfhD17YNAgGDIE6tYNOyoRqQKU+GPJ3r2wezfUqAGnnQbdu/sTt40bhx2ZiFQhKvXEig8/hAsvhCef9Nvdu0N2tpK+iFQ6Jf6wffUV9OsHF1/sp2c2bx52RCJSxanUE6bsbLjzTiguhl/+0tfxda9REQmYEn8Ydu3ydfzmzeGqq2DUKGjSJOyoRCRBqNQTTfn50KOHH+UDZGTAlClK+iISVUr80bBzJ2dOmOBH+P/4B7RoceBKXBGRKFOpJ2gLF0Lv3jT5/HPo2ROee053QBGRUGnEH7RTToGUFBY9/TT8+c9K+iISOiX+yuYcjB3ra/l79/rF1HJz2dy+fdiRiYgAISV+M1ttZp+ZWZ6ZLQgjhkAsX+7XxL/zTti8GbZt8+1aV0dEYkiYI/7OzrlWzrk2IcZQOXbt8lfcZmTAZ5/Ba69BTg6ceGLYkYmIHEYndytDcfGB8s7o0ZCWFnZEIiJHZC6EaYVmtgrYAjjgd865MaW8ZiAwECAtLS0zOzs7ukEeRbWiIhr++c+s7dWLvTVqUH3bNnbXqXPE1xcUFJCqq3L3U38coL4oSf1RUkX6o3PnzrmlVlWcc1F/AKdHfp4CLAQu+77XZ2ZmupgyZ45zzZo5B85NmVKmt+Tk5AQbU5xRfxygvihJ/VFSRfoDWOBKyamh1Pidc+siPzcAbwHtwojjmH37LTzwAHTs6GfszJrlyzsiInEk6onfzGqZWe19z4GrgMXRjqNcfvhDeOEFuO8+WLRItz4UkbgUxsndNOAt81McqwN/cs7NCCGOsikq8rN2TjgBHn8c7rnHj/hFROJU1BO/c+4/QEa0j1suCxbA7bdDq1bwpz/BeeeFHZGISIXpyt3S7N4Nw4bBRRfB9u3+huciIlWE5vEfavVq6NMH5s2D226Dl17ShVgiUqUo8R+qRg3YuNHfHat377CjERGpdCr1gF9X58kn/RTN006DZcuU9EWkylLinzMHzj8f/vd//clcgOr6IiQiVVfiJv7du+Gxx/xc/Fq14MMPoV18XEcmIlIRiTu07dfP1/HvuMNflKW1QUQkQSRe4nfOr49/773QvbufwSMikkASJ/Hv2gVDhkBSEvz619Chg3+IiCSYKlvjz8qC9HSoVg0ubriGjeddBs8+65dgCGEpahGRWFElR/xZWTBwIBQWwjW8wx/X9SOZYt67fzKXjb457PBEREJVJUf8Q4f6pH8yG3iTm/mSM7iAT7h9mpK+iEiVHPGvWeN/fsMpdGUG82lLEcdja8KNS0QkFlTJEf+ZZx54/j6XUcTxh7WLiCSqKpn4hw+HlJSSbSkpvl1EJNFVycTfty+MGQONGvkp+40a+e2+fcOOTEQkfFWyxg8+ySvRi4gcrkqO+EVE5MiU+EVEEowSv4hIglHiFxFJMEr8IiIJxlwcLFhmZt8A+WHHUUH1gY1hBxFD1B8HqC9KUn+UVJH+aOScO/nQxrhI/FWBmS1wzrUJO45Yof44QH1RkvqjpCD6Q6UeEZEEo8QvIpJglPijZ0zYAcQY9ccB6ouS1B8lVXp/qMYvIpJgNOIXEUkwSvwiIglGiT9gZnaGmeWY2VIzW2JmD4QdU9jMLMnMPjWzt8OOJWxmVtfM3jSz5Wa2zMwuCjumsJjZzyP/jyw2szfMrGbYMUWTmb1uZhvMbPFBbSeZ2btmtiLy88TKOJYSf/B2Aw85584DLgR+ambnhRxT2B4AloUdRIwYDcxwzp0LZJCg/WJmpwP3A22ccy2BJODWcKOKunFA10PaBgMznXNnATMj2xWmxB8w59x659wnkec78P9jnx5uVOExs4bAtcBrYccSNjOrA1wGjAVwzu1yzm0NNahwVQeON7PqQArwVcjxRJVz7j1g8yHNNwDjI8/HAzdWxrGU+KPIzNKB1sBHIYcSpueBh4G9IccRCxoD3wC/j5S+XjOzWmEHFQbn3DpgFLAGWA9sc879I9yoYkKac2595Pl/gbTK+FAl/igxs1Tgz8DPnHPbw44nDGZ2HbDBOZcbdiwxojpwAfBb51xr4Fsq6at8vInUrm/A/zE8DahlZj8IN6rY4vzc+0qZf6/EHwVmloxP+lnOuSlhxxOiDkB3M1sNZANdzGxCuCGFai2w1jm37xvgm/g/BInoCmCVc+4b51wxMAW4OOSYYsHXZtYAIPJzQ2V8qBJ/wMzM8DXcZc65Z8OOJ0zOuSHOuYbOuXT8ibtZzrmEHdU55/4LfGlm50SaLgeWhhhSmNYAF5pZSuT/mctJ0BPdh5gO9I887w9Mq4wPVeIPXgegH350mxd5dAs7KIkZ9wFZZrYIaAWMCDeccES+9bwJfAJ8hs9NCbV0g5m9AXwInGNma83sR8DTwJVmtgL/rejpSjmWlmwQEUksGvGLiCQYJX4RkQSjxC8ikmCU+EVEEowSv4hIglHiFxFJMEr8IiIJRolfpBzMrK2ZLTKzmmZWK7KOfMuw4xIpC13AJVJOZvYkUBM4Hr/mzlMhhyRSJkr8IuVkZjWA+UARcLFzbk/IIYmUiUo9IuVXD0gFauNH/iJxQSN+kXIys+n45aUbAw2cc/eGHJJImVQPOwCReGRmtwPFzrk/mVkSMNfMujjnZoUdm8jRaMQvIpJgVOMXEUkwSvwiIglGiV9EJMEo8YuIJBglfhGRBKPELyKSYJT4RUQSzP8Hk6brbufQyysAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.interpolate import interp1d\n",
    "x = np.array([1,3,7,10])\n",
    "y = np.array([2,9,20,35])\n",
    "plt.plot(x,y,'bo',label='data')\n",
    "\n",
    "xp = np.linspace(1,10,100)\n",
    "y1 = interp1d(x,y,kind='cubic')\n",
    "print(y1(5))\n",
    "plt.plot(xp,y1(xp),'r--',label='cubic')\n",
    "plt.plot(5,y1(5),marker='o',linestyle='',color='orange',label='evaluate')\n",
    "plt.legend(); plt.grid(); plt.xlabel('x'); plt.ylabel('y')\n",
    "plt.savefig('plot.png',dpi=600)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 6: Nonlinear Regression and Plotting\n",
    "\n",
    "Pressure in a tank is recorded as $P = [1.5,2.6,3.5,10.2,20.3,30.2]$ at respective times of $[0,1,2,3,4,5]$ min. Create a nonlinear approximation of the pressure trend as:\n",
    "\n",
    "$P = A \\, e^{n\\,t}$\n",
    "\n",
    "where $t$ is time, $A$ is an unknown pre-exponential factor, and $n$ is also an unknown parameter. Show the optimized parameter values as well as a plot with the data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = 1.8338873466478285\n",
      "n = 0.5669820526775803\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2a2f98f8708>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcoUlEQVR4nO3deXSU5fnG8e8NghipS0uwqEBccRc1uKIiVnGrorjjjlJbF1DrBvWntoLUXWtdgihRgrigBXEvglRRMCCyI4qAKEoEFyjKEu7fH8+kgCZkkszknXfm+pzznsm8M2TuOR4uHp/V3B0REYmfBlEXICIitaMAFxGJKQW4iEhMKcBFRGJKAS4iElMb1eeHNWvWzAsKCurzI0VEYm/ChAnfuHv+z+/Xa4AXFBRQWlpanx8pIhJ7ZjavsvvqQhERiSkFuIhITCnARURiSgEuIhJTCnARkZhSgIuIxJQCXEQkphTgIiLp9OOP0KMHzJ+f8l+tABcRSad//hMeeADmzk35r1aAi4ikyw8/wO23Q6dOcNhhKf/1CnARkXS55x5YsgT69EnLr1eAi4ikwzffwN13w6mnwn77peUjFOAiIunQrx8sX85L7f5KQQE0aAAFBVBSkrqPqDbAzayJmY03s4/MbJqZ3Zq4v52ZjTOzT8zsGTNrnLqyRERibMECePBBPj3kPM68dVfmzQN3mDcPundPXYgn0wJfAXR0972BtsAxZnYg8HfgXnffEfgW6JaakkREYu6WW8Cd8+fczPLl67+0fDn07p2aj6k2wD1YlnjaKHE50BF4PnG/GOicmpJERGJs+nR44gm47DLGfllQ6VtSNSU8qT5wM2toZpOARcCbwKfAd+6+OvGWBcA2VfzZ7mZWamalZWVlKShZRCSD9eoFTZtCr160alX5W6q6X1NJBbi7l7t7W2BbYH9gl2Q/wN2L3L3Q3Qvz839xIpCISPZ4910YNgyuvx6aNaNPH8jLW/8teXmpm1VYo1ko7v4dMAo4CNjCzCqOZNsW+CI1JYmIxJB7CO4WLaBnTwC6doWiImjdGszCY1FRuJ8K1Z6JaWb5wCp3/87MNgGOIgxgjgJOBYYA5wPDUlOSiEgMDR8eWuBFRes1u7t2TV1g/1wyhxq3AIrNrCGhxf6su48ws+nAEDO7DfgQGJCeEkVEMtyqVXDddbDrrnDhhfX2sdUGuLtPBvap5P4cQn+4iEhuKyqCjz+GESNgo2TaxamhlZgiInXx/fdh3nfHjnDccfX60QpwEZG6uP12WLwY7rorjFTWIwW4iEhtzZsH990H554L+/yipzntFOAiIrV1442h1X3bbZF8vAJcRKQ2xo6Fp5+Ga6+Fli0jKUEBLiJSU2vWhMU6W28dFu9EpP7mu4iIZIuSEvjgA3jySdh008jKUAtcRKQmli2DG26Adu3St8QySWqBi4jUxN//Dl9+Cc89F47ZiZBa4CIiyZozB+68E84+Gw4+OOpqFOAiIkm7+uqwVP6OO6KuBFAXiohIcl5/Pez1ffvtsE2l59fUO7XARUSqs3Il9OgBO+0EV10VdTX/oxa4iEh1HngAZs2Cl1+GjTeOupr/UQtcRGRDFiwIuw2ecEK97zZYHQW4iMiGXH01lJeHVniGUYCLiFTljTfCfO9evWC77aKu5hcU4CIilVmxAq64AnbYIWxYlYE0iCkiUpm77w7HpL36KjRpEnU1lVILXETk5z79FP72N+jSBY45JupqqqQAFxFZlztcdhk0agT33x91NRukLhQRkXU980xYdfnAAxmz4rIqaoGLiFT47rtwUENhIfzpT1FXUy21wEVEKtx4I5SVwSuvQMOGUVdTLbXARUQA3nkHHnkk7Hmy775RV5OUagPczFqa2Sgzm25m08ysR+L+LWb2hZlNSlyZtcZURCRZP/0El1wCBQVh9klMJNOFshq4xt0nmtmvgAlm9mbitXvd/a70lSciUg/69oWZM8PgZYRnXNZUtQHu7guBhYmfl5rZDCCzh2ZFRJI1ZUrY4/u88+Doo6OupkZq1AduZgXAPsC4xK3LzWyymT1uZltW8We6m1mpmZWWlZXVrVoRkVRavRouvhi23BLuuSfqamos6QA3s6bAUKCnu/8APAzsALQltNDvruzPuXuRuxe6e2F+fn7dKxYRSZV774Xx48Oc79/8JupqaiypADezRoTwLnH3FwDc/Wt3L3f3NUB/YP/0lSkikmKzZsFNN0HnznDGGVFXUyvJzEIxYAAww93vWed+i3XedjIwNfXliYikQXk5XHQR5OXBQw+BWdQV1Uoys1AOAc4FppjZpMS9XsBZZtYWcGAu8Ic01CciknoPPghjx0JxMbRoUf37M1Qys1DeASr75+mV1JcjIpJms2eHFZfHHQfnnht1NXWilZgikjvKy+H888P+3v37x7brpIL2QhGR3HH33fDeezBoEGy9ddTV1Jla4CKSG6ZNC7NOTjkFzj476mpSQgEuItlv1aqw0nKzzeDhh2PfdVJBXSgikv1uvRUmToShQ6F586irSRm1wEUku40dG/Y6ueCC0H2SRRTgIpK9li4NUwVbtcr48y1rQ10oIpK9rroK5s6Ft98O/d9ZRi1wEclOQ4fCgAFw/fXQvn3U1aSFAlxEss+CBeGEnXbtwgBmllKAi0h2KS8P/d4rV8LgwdCoUdQVpY36wEUku9xxB4weDQMHwo47Rl1NWqkFLiLZ4/334f/+D848MyzcyXIKcBHJDt9+G4K7ZcusWm25IepCEZH4c4du3eCLL+Ddd2GLLaKuqF4owEUk/h56CF58Mew2uH/unO6oLhQRibeJE+Hqq+H448PCnRyiABeR+PruOzj11LBB1cCBOdHvvS51oYhIPLnDhRfC55/DmDHQrFnUFdU7BbiIxNO998K//gX33AMHHRR1NZFQF4qIxM+774Y9Tk4+GXr2jLqayCjARSRevvoKTjsNCgrg8cdzrt97XQpwEak3JSUhdxs0CI8lJTX8BatWwemnh8HLF17ImfneVVEfuIjUi5IS6N4dli8Pz+fNC88BunZN8pdcdx385z/hl+25Z1rqjJNqW+Bm1tLMRpnZdDObZmY9Evd/bWZvmtnsxOOW6S9XROKqd++14V1h+fJwPymDB8N998GVV2bNqfJ1lUwXymrgGnffDTgQuMzMdgNuAEa6+07AyMRzEZFKzZ9fs/vrmTgxLJU/9FC4666U1hVn1Qa4uy9094mJn5cCM4BtgJOA4sTbioHOaapRRLJAq1Y1u/8/ixaF2SbNmsHzz2f1/t41VaNBTDMrAPYBxgFbufvCxEtfAVtV8We6m1mpmZWWlZXVpVYRibE+fSAvb/17eXnhfpUqBi0XLQpzvps3T2eJsZN0gJtZU2Ao0NPdf1j3NXd3wCv7c+5e5O6F7l6Yn59fp2JFJL66doWiImjdOsz8a906PN/gAGbPnuFA4v79Yb/96qvU2EhqFoqZNSKEd4m7v5C4/bWZtXD3hWbWAliUriJFJDt07VqDGScPPRSua6+Fc85Ja11xlcwsFAMGADPc/Z51XhoOnJ/4+XxgWOrLE5Gc9NZbYbbJCSfA7bdHXU3GSqYFfghwLjDFzCYl7vUC+gHPmlk3YB5weloqFJHc8sknYYfBXXYJ870bNoy6ooxVbYC7+ztAVWtVj0xtOSKS05YsCft6N2gAw4fDZptFXVFG00pMEckMK1dCly4wdy6MHAnbbx91RRlPAS4i0XMP6+pHj4ZBg6B9+6grigVtZiUi0evTB4qL4eabazBNRRTgIhKtp56Cm24KwX3zzVFXEysKcBGJzsiRcNFFcMQROb+3d20owEUkGlOmwCmnQJs2YW/vxo2jrih2FOAiUv/mz4djj4WmTeHVV3P+YIba0iwUEalfixdDp06wbFk4Tb5ly6grii0FuIjUn+XLw/L4zz6DN96AvfaKuqJYU4CLSP2o2Bp2/Hh47jk47LCoK4o9BbiIpN+aNXDBBfDyy/Doo2HwUupMg5gikl7ucMUV4UzLfv3WnmQsdaYAF5H0uummsK/3ddfB9ddHXU1WUYCLSPr06xeWyV9ySfhZUkoBLiLp8cADcOONcPbZ8PDDWmWZBgpwEUm9xx6DHj3CafLFxTqUIU0U4CKSWsXFYaDymGPg6adhI012SxcFuIikzqBBcOGFcOSRYX+TjTeOuqKspgAXkdQYMgTOPx86dIBhw2CTTaKuKOspwEWk7gYPDvt5t28PL70EeXlRV5QTFOAiUjeDBsG558Khh4aVlptuGnVFOUMBLiK1V1wM550Hhx8ewrtp06gryikKcBGpnaKitQOWI0ao5R0BBbiI1Nx998Ef/hAOZRg+XH3eEak2wM3scTNbZGZT17l3i5l9YWaTEtdx6S1TRDJG375w1VXQpQu8+KJmm0QomRb4QOCYSu7f6+5tE9crqS1LRDKOe9iMqndvOOecMG1Q51hGqtoAd/cxwJJ6qEVEMlV5eegyueMO+NOfwuClVlhGri594Jeb2eREF8uWVb3JzLqbWamZlZaVldXh40QkEitWhA2p+veHXr3gwQehgYbPMkFt/ys8DOwAtAUWAndX9UZ3L3L3QncvzM/Pr+XHiUgkli4NZ1g++yzceWfYGla7CmaMWv0/kLt/XfGzmfUHRqSsIhHJDIsWhVkmH30EAweGZfKSUWoV4GbWwt0XJp6eDEzd0PtFJGY+/TTsJvjFF2Ffk+OPj7oiqUS1AW5mTwMdgGZmtgC4GehgZm0BB+YCf0hfiSJSr8aNg9//PhxEPHIkHHRQ1BVJFaoNcHc/q5LbA9JQi4hEbdgwOOssaNECXn0Vdt456opkAzSULCLBgw/CKafAHnvAe+8pvGNAAS6S68rL4cor4YorwoyTUaOgefOoq5IkKMBFctnSpXDSSfCPf8DVV4dTdLQpVWxoKZVIrpo7F048EaZPh4cegj/+MeqKpIYU4CK56J13Qn/3ypVhsPKoo6KuSGpBXSgiuebxx6FjR9hiizBlUOEdWwpwkVyxahVcfjl06xZO0Bk3Dtq0iboqqQMFuEguKCsLLe1//hOuuSZ0m2xZ5R50EhPqAxfJduPHh8MXvvkmHEDctWvUFUmKqAUuks369w+nxW+0Ebz7rsI7yyjARbLR8uVw0UXQvTsccQSUlsK++0ZdlaSYAlwk28yeHTageuIJ+Mtf4OWX4Te/iboqSQP1gYtkk+eeC7NMGjcOA5XHVHacrWQLtcBFssGPP4aVlKefDrvvDh9+qPDOAQpwkbibORMOPBAeeQSuuw7GjIGWLaOuSuqBulBE4sodBgyAHj0gLw9eeSUcgSY5Qy1wkTj69tvQXXLJJaH1PWmSwjsHKcBF4mbUKNh7b/jXv6BfP3jzTdhmm6irkggowEXiYsUK+POf4cgjYZNNYOxYuP56aKC/xrlKfeAicTBpEpx3HkyZEmab3HmnDl4QtcBFMtrq1XDbbdCuXdiQasSIcPiCwltQC1wkc02bBhdeCB98EE6K/8c/tKJS1qMWuEimWbUqtLr32Qc++wyeeQYGD1Z4yy8owEUiUlICBQVhDLKgIDznww/hgAPgppvCkWfTp4fpgiKVqDbAzexxM1tkZlPXufdrM3vTzGYnHrUzvEgNlJSEjQLnzQvrcb6e9yNfXXADawrbwcKFMHQoDBkC+flRlyoZLJkW+EDg55sq3ACMdPedgJGJ5yKSpN69w46vAB0ZyWT24prVf+fZvAtCq/uUUyKtT+Kh2gB39zHAkp/dPgkoTvxcDHRObVki2W3+fMhnEU9yLiP5HRCC/Oz/PqajziRpte0D38rdFyZ+/grYqqo3mll3Mys1s9KysrJafpxIFikv54YtH2Umu3AGz/BXbmJPpjCKjrRqFXVxEid1HsR0dwd8A68XuXuhuxfmqz9Pct2ECXDQQfRdcilTG+zF3nzEzfyVFTQhLw/69Im6QImT2gb412bWAiDxuCh1JYlkocWL4dJLw4Kc+fNh0CA+Lx7Fj613xQxat4aiIh1ZKTVT24U8w4HzgX6Jx2Epq0gkm6xeHZL5L3+BH36AK6+EW2+FzTenK9D1nKgLlDhLZhrh08B7QBszW2Bm3QjBfZSZzQZ+l3guIut6661wkPBll0HbtmE/k/vug803j7gwyRbVtsDd/awqXjoyxbWIZIdPPoFrrw3bvRYUhHMqu3QBs6grkyyjlZgiqbJ4MfTsCbvuCv/+N/TtCzNmwKmnKrwlLbSZlUhd/fRT2Giqb9/Qz33xxaGf+7e/jboyyXJqgYvUVnk5DBwIO+8cDhM++GCYPBkefVThLfVCAS5SU+7w4ouw115hu9ff/jYcc/byy7D77lFXJzlEAS6SLPfQt33ggWGvkvJyePZZGDcOOnSIujrJQQpwkWS8/XYI6aOOCrsFDhgAU6fCaadpgFIiowAX2ZAxY8Ihwh06wOzZYbBy9my46CLYSHMAJFoKcJGfc4fRo+GII+Dww8PRZnffDZ9+CpdfDhtvHHWFIoCmEYqs5Q6vvRaOMxs7NgxO3ntvOHkhLy/q6kR+QQEuUl4eTsDp1y8cadaqFTz4YOgm2WSTqKsTqZICXHLXTz/Bk0/CnXeG5e877xwGJ885Bxo3jro6kWopwCX3LF4MDz0UWtmLFsF++8Hzz0PnztCwYdTViSRNAS65Y9assBtgcTH8+CMcdxz8+c9hhommAkoMKcAlu1Usvrn//rBScuONw6kJV10Fe+wRdXUidaIAl+y0dCkMGhTmbc+YAc2bwy23wB//GH4WyQIKcMkuM2aE/u3i4hDihYXw1FNhxaTmb0uWUYBL/K1cGTaXeuSRsACncWM4/fRwEs4BB6h/W7KWAlzia/ZseOyxsKXrokXh9Ju+faFbN3WTSE5QgEu8/PgjvPBCCO7Ro8O0vxNOCCe+H300NNDuEJI7FOCS+dzhgw/giSfg6afh++9h++2hTx+44ALYeuuoKxSJhAJcMteCBVBSErpIZs4My9pPPTUscT/sMLW2JecpwCWz/PBD6CJ56qlwyo07tG8P/fuHmSSbbx51hSIZQwEu0VuxAl59FQYPhpdeCnuU7LAD3HxzWHSz445RVyiSkRTgEo1Vq+Ctt2DIkDAF8PvvIT8/nOh+9tnh2DJN/xPZoDoFuJnNBZYC5cBqdy9MRVGSpVatCt0izz0XQnvxYthsMzj5ZDjjDPjd76BRo6irFImNVLTAj3D3b1LweyQb/fRT2Itk6FAYPhyWLIGmTeHEE8Nim06doEmTqKsUiSV1oUjqff89vPJKaGW/+iosWxYGH088Ebp0UWiLpEhdA9yBN8zMgUfdvSgFNUkczZsXBiCHDQsLbFavhq22Cv3ZJ58MHTvqkASRFKtrgLd39y/MrDnwppnNdPcx677BzLoD3QFatWpVx4+TjLF6NYwbByNGhGvq1HB/l13CVq2dO4eBSM3VFkkbc/fU/CKzW4Bl7n5XVe8pLCz00tLSlHyeRODrr+GNN0L3yOuvw7ffwkYbwaGHwvHHw+9/H44lE5GUMrMJlU0SqXUL3Mw2BRq4+9LEz0cDf61DjZJpVqyA994Lof3aa+HAXwhdIyedBMceG/Yf2WKLSMsUyVV16ULZCnjRwlzdjYDB7v5aSqqSaLiHrpB//ztco0fD8uVhw6hDDgl7j3TqBPvso64RkQxQ6wB39znA3imsReqbO8yZExbUvPVWmKP99dfhtZ13hgsvDC3sDh3CfG0RySiaRphL3OHTT+Htt0PrevTosGEUQIsWcOSRYTHNkUeCBpxFMp4CPJutWRO6RN55B/7zHxgzBr78MrzWvHloWR9+eJji16aNlq6LxIwCPJv8979h3+x33w3X2LFhUQ3ANtuEsD7ssHDtumtGBXZJCfTuDfPnh8Z/nz5hHysRqZoCPK4qukPefz/MFHnvPZg8GcrLw+u77RaWqh96aNiOtaAgowJ7XSUl0L17GC+FsCaoe/fws0JcpGopmweeDM0Dr4OystC6Hj9+7bV4cXitaVPYf3846KC1169/HW29NVBQEEL751q3hrlz67sakcyT8nngkkZLlsCECeEqLQ1XRcKZhdb1SSeFlY4HHAC77x6m+sXU/Pk1uy8igQI8Su6wcGFYIFNxTZy4frNzhx1CUF9+ObRrB/vuC7/6VWQlp0OrVpW3wDURRmTDFOD1ZeVKmDEj9FN/9NHaq6xs7Xt22il0hVx6Key3XwjrGHWF1FafPuv3gQPk5YX7IlI1BXiqrVkTWtBTp669Jk+GWbPCBlAAG28Me+wRtldt2xb23jtcicUyJSXQ++LcmZFR8d00C0WkZjSIWVvl5SGop08PLetp08I1Y8b6TclWrWCvvWDPPcO1995hleNGlf/b+fMZGRBao0VFCjSRXFXVIKYCvDpLl8LHH4cW9KxZMHNmCOmPPw6bPVVo0SIMJlZce+4ZBhtruARdMzJE5Oc0C2VDfvop7Akye/baqyK0Fy5c+74GDWC77cIimE6dwt7Xu+0WnqdoRz7NyBCRZOVOgH/7bQjpOXPgk0/CIpiKa8GCMCOkQrNmoZujU6fw2KZNuHbcMfRfp5FmZIhIsrInwJctC8k3dy589tnaxzlzwmPFkvIKzZuHKXodOoRg3mGHMAtkp51gyy0j+AKBZmSISLLiEeBr1sBXX8Hnn4e+hIrHefPC49y5YfHLupo0CR3K220HBx8cAnr77cPz7bfP2LnUmpEhIsnK+AAvKYGGl3bnzGUD1n+hadMwsteqVZg73bp1uLbbLgT3Vltl7N4f1enaVYEtItXL6ACvmFJ34PKzeJv9+JyWlDVpxZ/vb8lpl2wR24AWEUmFjJ5GqCl1IiJVTyPM6IMNNaVORKRqGR3gVU2d05Q6EZEMD/A+fcIUunVpSp2ISJDRAd61a9gDpHXrMF7ZurX2BBERqZDRs1BAU+pERKqS0S1wERGpmgJcRCSmFOAiIjGlABcRiSkFuIhITNXrUnozKwMqWRyflGbANyksJw70nXODvnNuqMt3bu3u+T+/Wa8BXhdmVlrZXgDZTN85N+g754Z0fGd1oYiIxJQCXEQkpuIU4EVRFxABfefcoO+cG1L+nWPTBy4iIuuLUwtcRETWoQAXEYmpWAS4mR1jZrPM7BMzuyHqetLNzB43s0VmNjXqWuqDmbU0s1FmNt3MpplZj6hrSjcza2Jm483so8R3vjXqmuqLmTU0sw/NbETUtdQHM5trZlPMbJKZJX+mZDK/O9P7wM2sIfAxcBSwAPgAOMvdp0daWBqZ2WHAMuBJd98j6nrSzcxaAC3cfaKZ/QqYAHTO8v/GBmzq7svMrBHwDtDD3d+PuLS0M7OrgUJgM3c/Iep60s3M5gKF7p7yhUtxaIHvD3zi7nPcfSUwBDgp4prSyt3HAEuirqO+uPtCd5+Y+HkpMAPYJtqq0suDZYmnjRJXZremUsDMtgWOBx6LupZsEIcA3wb4fJ3nC8jyv9y5zMwKgH2AcRGXknaJroRJwCLgTXfP+u8M3AdcB6yJuI765MAbZjbBzLqn8hfHIcAlR5hZU2Ao0NPdf4i6nnRz93J3bwtsC+xvZlndXWZmJwCL3H1C1LXUs/buvi9wLHBZoos0JeIQ4F8ALdd5vm3inmSRRD/wUKDE3V+Iup765O7fAaOAYyIuJd0OAU5M9AkPATqa2aBoS0o/d/8i8bgIeJHQLZwScQjwD4CdzGw7M2sMnAkMj7gmSaHEgN4AYIa73xN1PfXBzPLNbIvEz5sQBulnRlpUmrn7je6+rbsXEP4ev+Xu50RcVlqZ2aaJgXnMbFPgaCBls8syPsDdfTVwOfA6YXDrWXefFm1V6WVmTwPvAW3MbIGZdYu6pjQ7BDiX0CKblLiOi7qoNGsBjDKzyYRGypvunhPT6nLMVsA7ZvYRMB542d1fS9Uvz/hphCIiUrmMb4GLiEjlFOAiIjGlABcRiSkFuIhITCnARURiSgEuIhJTCnARkZj6f7z+duVKyYcwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "x = np.array([0,1,2,3,4,5])\n",
    "y = np.array([1.5,2.6,3.5,10.2,20.3,30.2])\n",
    "\n",
    "# Step 1: define nonlinear function\n",
    "# Step 2: use curve_fit to optimize the parameters\n",
    "# Step 3: plot the function with the optimized parameters and data\n",
    "\n",
    "\n",
    "## Solution\n",
    "def f(x,A,n):\n",
    "    return A*np.exp(n*x)\n",
    "\n",
    "popt, pcov = curve_fit(f,x,y)\n",
    "A = popt[0]\n",
    "n = popt[1]\n",
    "print('A = '+str(A))\n",
    "print('n = '+str(n))\n",
    "\n",
    "plt.plot(x,y,'bo')\n",
    "xp = np.linspace(0,5,100)\n",
    "plt.plot(xp,f(xp,A,n),'r-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 7: Equation Solution and Plotting\n",
    "\n",
    "#### Part A\n",
    "\n",
    "The ideal gas law is shown below:\n",
    "\n",
    "$P V_m=R_g T$\n",
    "\n",
    "where $P$ is the pressure, $V_m$ is the molar volume, $T$ is the temperature, and $R_g$ is the universal gas constant. Use $R_g=8.314 \\frac{J}{mol\\,K}$.\n",
    "\n",
    "Create the function $P(V_m,T)$ such that the pressure ($P$) is a function of molar volume $\\left(V_m\\right)$ and temperature ($T$). Convert the pressure value from $Pa$ to $bar$.\n",
    "\n",
    "Use the function and range variables to create a $P$ vs $T$ plot where $T$ ranges from $100K$ to $1200K$ and $V_m$ = 2.24 $\\frac{L}{mol}$ = 0.00224 $\\frac{m^3}{mol}$ . Include $x$ and $y$ labels on the plot and a legend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvO0lEQVR4nO3dd5gUVdbH8e8hKC5ZGBFFxSwMmUFFxMCa0xpXXROLgIoKGIFXXTDgLiuKYkIUI7ogmBADiwKC2QHJGNCVFUQZUdIShznvH7cGRxiGHpiemp7+fZ6nH/pWhzpFwenbt2+da+6OiIikjwpxByAiIqVLiV9EJM0o8YuIpBklfhGRNKPELyKSZirFHUAi6tat6w0bNow7DBGRlDJ16tSf3T1j8+0pkfgbNmxIdnZ23GGIiKQUM1tQ2HYN9YiIpBklfhGRNKPELyKSZpI+xm9mFYFsYJG7n2ZmTwNHA8ujp3R09+nFfd8NGzawcOFC1q5dW2KxSmqoUqUKDRo0oHLlynGHIpKSSuPH3R7APKBGgW03ufvoHXnThQsXUr16dRo2bIiZ7VCAkjrcnaVLl7Jw4UL23XffuMMRSUlJHeoxswbAqcATJf3ea9eupU6dOkr6acbMqFOnjr7pieyAZI/x3w/cDORttr2/mc00s0FmtnNhLzSzrmaWbWbZOTk5hb65kn560nmXtLF+fVLeNmmJ38xOA5a4+9TNHuoDHAK0AXYFehX2encf6u5Z7p6VkbHF9QciIuXbW2/BgQfCe++V+Fsns8ffDjjDzL4DRgAdzGy4uy/2YB3wFHBoEmNIqmrVqhW6vWPHjowevX0/YfTr14+BAwcW+tjw4cNp1qwZmZmZNG/enM6dO7Ns2bLt2k++SZMmUbNmTVq0aEGjRo24/fbbd+j9RGQHLV0Kl14Kp5wC//0vPPpoie8iaYnf3fu4ewN3bwhcAExw94vNrD6Ahe/rZwKzkxVDefL2228zaNAg3nrrLebMmcO0adM44ogj+Omnn3b4vdu3b8/06dPJzs5m+PDhTJs27XeP5+bm7vA+ElWa+xIpU9zhxRehUSN47jmoUgXuuQeGDy/xXcUxj/95M5sFzALqAnfFEEOJcneuueYaDj74YI477jiWLFmy6bGpU6dy9NFH07p1a0488UQWL14MwOOPP06bNm1o3rw555xzDqtXry5yH/3792fgwIHsueeeAFSsWJFOnTpx8MEHA3DHHXfQpk0bmjRpQteuXclfWW3w4ME0btyYZs2accEFFxS5j6pVq9K6dWvmz59Pv379uOSSS2jXrh2XXHIJOTk5nHPOObRp04Y2bdrwwQcfAPDee+/RokULWrRoQcuWLVm5ciWLFy/mqKOOokWLFjRp0oQpU6YAv/+GNHr0aDp27AiEb0hXXnklhx12GDfffDPffPMNJ510Eq1bt6Z9+/Z88cUXiZ4KkdT0ww9w9tlw/vmQkwNHHw2zZsGNN0KlJEy+dPcyf2vdurVvbu7cub81wmdlyd+2oWrVqu7u/tJLL/lxxx3nubm5vmjRIq9Zs6aPGjXK169f723btvUlS5a4u/uIESP8r3/9q7u7//zzz5ve55ZbbvHBgwe7u3vfvn39nnvu2WJftWvX9mXLlm01lqVLl266f/HFF/uYMWPc3b1+/fq+du1ad3f/9ddft3jdxIkT/dRTT90U0z777OOzZ8/2vn37eqtWrXz16tXu7n7hhRf6lClT3N19wYIFfsghh7i7+2mnnebvv/++u7uvXLnSN2zY4AMHDvS77rrL3d1zc3N9xYoVv/v7cncfNWqUX3bZZe7uftlll/mpp57qubm57u7eoUMH/+qrr9zd/eOPP/Zjjz12i7h/d/5FUlVenvuwYe41a4acU726+5Ah7hs3lsjbA9leSE5NiSJtZd3kyZO58MILqVixInvssQcdOnQA4Msvv2T27Nkcf/zxAGzcuJH69esDMHv2bG699VaWLVvGqlWrOPHEExPe36xZs7jkkktYuXIld999N+effz4TJ07kn//8J6tXr+aXX34hMzOT008/nWbNmnHRRRdx5plncuaZZxb6flOmTKFly5ZUqFCB3r17k5mZyahRozjjjDPYZZddAHjnnXeYO3fuptesWLGCVatW0a5dO66//nouuugizj77bBo0aECbNm3o1KkTGzZs4Mwzz6RFixbbPKbzzjuPihUrsmrVKj788EPOO++8TY+tW7cu4b8bkZTx7bfQtSu8+25on3oqDBkCDRokfdflI/GX0QXj3Z3MzEw++uijLR7r2LEjr776Ks2bN+fpp59m0qRJRb5XZmYm06ZN49hjj6Vp06ZMnz6da665hjVr1rB27Vq6detGdnY2e+21F/369ds0z/2NN95g8uTJvP766/Tv359Zs2ZRabOvju3bt2fs2LFb7LNq1aqb7ufl5fHxxx9TpUqV3z2nd+/enHrqqbz55pu0a9eOcePGcdRRRzF58mTeeOMNOnbsyPXXX8+ll176u2mYm8/Dz99XXl4etWrVYvr06UX+fYikrI0b4cEH4ZZbYPVqqFMHHngA/vIXKKWpyqrVUwKOOuooRo4cycaNG1m8eDETJ04E4OCDDyYnJ2dT4t+wYQNz5swBYOXKldSvX58NGzbw/PPPb3Mfffr04cYbb2ThwoWbtq1Zswb4LYnWrVuXVatWbZpRlJeXx/fff8+xxx7LgAEDWL58OatWrdquYzzhhBN48MEHN7XzE/M333xD06ZN6dWrF23atOGLL75gwYIF1KtXjy5dutC5c+dNPxbXq1ePefPmkZeXxyuvvFLofmrUqMG+++7LqFGjgPDhOWPGjO2KWaTMmTMH2rWD664LSf+CC2DePLjoolJL+lBeevwxO+uss5gwYQKNGzdm7733pm3btgDstNNOjB49mu7du7N8+XJyc3Pp2bMnmZmZ3HnnnRx22GFkZGRw2GGHsXLlyiL3ccopp5CTk8PJJ5/Mxo0bqVWrFk2aNOHEE0+kVq1adOnShSZNmrD77rvTpk0bIAwtXXzxxSxfvhx3p3v37tSqVWu7jnHw4MFcffXVNGvWjNzcXI466iiGDBnC/fffz8SJE6lQoQKZmZmcfPLJjBgxgnvuuYfKlStTrVo1nn32WQD+8Y9/cNppp5GRkUFWVtZWP4Sef/55rrrqKu666y42bNjABRdcQPPmzbcrbpEyYf16GDAA7rwTNmyAPfYI0zTPOCOWcMzL6DBJQVlZWb75Qizz5s2jUaNGMUUkcdP5l5Tx2Wdw+eVhlg5Aly5hmmbNmknftZlNdfeszbdrqEdEJBlWr4abb4bDDw9Jf7/9wg+5Q4eWStIvioZ6RERK2nvvQefOMH8+VKgAN9wAd9wBf/hD3JEBKZ743V0Fu9JQKgxPSppasQJ69QrTMgGaNIFhw+DQslWZJmWHeqpUqcLSpUuVBNKMR/X4N59WKhK7N96Axo1D0q9cGW6/HaZOLXNJH1K4x9+gQQMWLlzI1ko2S/mVvwKXSJmQkwM9e8ILL4T2oYeGXn6TJrGGVZSUTfyVK1fWCkwiEh93GDkSrr0Wfv4ZdtkF+veH7t2hYsW4oytSyiZ+EZHYLFoEV10Fr78e2sceC48/DvvvH29cCUrZMX4RkVLnHhJ848Yh6deoEdrvvpsySR/U4xcRScw334SLr6KSLJx+erj6NiqVnkrU4xcRKcrGjXDffdC0aUj6GRkwYgS89lpKJn0ohcRvZhXN7HMzGxu19zWzT8xsvpmNNLOdkh2DiMh2mT0bjjgiXIC1Zk0opjZ3blgwJYWvISqNHn8PYF6B9gBgkLsfAPwKXF4KMYiIJG79eujXD1q1gk8/DTXyx44NyyDWrRt3dDssqYnfzBoApwJPRG0DOgD5K5E/Q1h3V0SkbPj005Dwb789VNK84opQTvnUU+OOrMQku8d/P3AzkBe16wDL3D1/Re2FQKGDZGbW1cyyzSxbF2mJSNKtXh2GdNq2DYn+gANg0qRwJW6NGnFHV6KSlvjN7DRgibtP3Z7Xu/tQd89y96yMjIwSjk5EpICJE8OPt/fdF9o33QQzZoRFz8uhZE7nbAecYWanAFWAGsADQC0zqxT1+hsAi5IYg4jI1i1bFkonP/54aDdtCk8+CVlblLAvV5LW43f3Pu7ewN0bAhcAE9z9ImAicG70tMuA15IVg4jIVo0ZA5mZIelXrhzKJmdnl/ukD/HM4+8FXG9m8wlj/sNiiEFE0tWSJWGt2z/9CX74ISyUMn063HYb7JQes8tL5cpdd58ETIrufwuUvTqlIlK+uYcKmj16wNKlYVGUu++Ga64p80XVSppKNohI+ff993DllfDmm6F93HFhCcQ0rfCrkg0iUn7l5YV6OpmZIenXrBlq5f/732mb9EE9fhEpr77+Oqx7O3lyaJ95Jjz8MOyxR6xhlQXq8YtI+ZKbC/fcA82ahaS/224wahS8/LKSfkQ9fhEpP2bMgMsvD2vdAlx6abgoq06deOMqY9TjF5HUt25dmI6ZlRWS/l57wVtvwTPPKOkXQj1+EUltH30UevnzoiLAV18Nf/87VK8eb1xlmHr8IpKaVq2Cnj2hXbuQ9A86KIzpP/SQkv42KPGLSOoZPz7U1XngAahQAfr0CeP77dvHHVlK0FCPiKSOX3+FG28MhdQAWrQI8/JbtYo1rFSjHr+IpIZXXoHGjUPS33nnUG4hf9EUKRb1+EWkbPvpp1BPZ3S0cF+7dvDEE3DIIfHGlcLU4xeRsskdnn0WGjUKSb9qVXjwwfADrpL+DlGPX0TKngULwlq348aF9gknhKJq++wTb1zlhHr8IlJ25OWFejpNmoSkX7s2PP00vP22kn4JSlqP38yqAJOBnaP9jHb3vmb2NHA0sDx6akd3n56sOEQkRXz5ZSiq9v77oX3uuWFoZ/fd442rHErmUM86oIO7rzKzysD7ZvZW9NhN7j46ifsWkVSxYQPcey/06xdKL9SrB488AmefHXdk5VbSEr+7O7AqalaObp6s/YlICvr881Bu4fPPQ7tjx1BUrXbtWMMq75I6xm9mFc1sOrAEGO/un0QP9TezmWY2yMx23spru5pZtpll5+TkJDNMESlta9fCLbdAmzYh6e+zTxjTf+opJf1SkNTE7+4b3b0F0AA41MyaAH2AQ4A2wK6ExdcLe+1Qd89y96yMjIxkhikipenDD6Fly3ABVl4edO8Os2eHmTtSKkplVo+7LwMmAie5+2IP1gFPoYXXRdLDqlUhyR95JHzxRZiL//77od5OtWpxR5dWkpb4zSzDzGpF93cBjge+MLP60TYDzgRmJysGESkjxo0LUzQffBAqVoT/+78wxHPEEXFHlpaSOaunPvCMmVUkfMC86O5jzWyCmWUABkwHrkxiDCISp19+geuvDwuiQKirM2xYKK4msUnmrJ6ZQMtCtndI1j5FpAx56aWwKMpPP4WiarffDjfcAJVUMCBuOgMiUrIWLw5F1V5+ObTbtw9F1Q46KN64ZBOVbBCRkuEepmM2bhySfrVqofzCpElK+mWMevwisuO++w66dg0rYwGcdBI89hjsvXesYUnh1OMXke23cSMMHhxm7IwfD7vuCs89B2++qaRfhqnHLyLbZ968UFTtww9D+89/DtM1d9st3rhkm9TjF5Hi2bAB+vcPUzI//BDq1w/LIo4cqaSfItTjF5HETZsGnTrBjBmh3bkz3HMP1KoVa1hSPOrxi8i2rVkDvXvDoYeGpL/vvvDOO/D440r6KUg9fhEp2pQpoWf/1VdgBtddB3feGdbAlZSkxC8ihVu5MvTyH3kktBs3DuUWDj883rhkh2moR0S29NZbkJkZkn6lStC3bxjfV9IvF9TjF5HfLF0ahnKeey60s7LgySehadN445ISpR6/iIRyCy++CI0ahaRfpQoMHAgffaSkXw6pxy+S7n74Abp1g9deC+2jjw5F1Q44IN64JGkSSvxmVgFoDuwBrAFmu/uSZAYmIknmHoZxbrgBli+H6tVDL79zZ6igwYDyrMjEb2b7E9bEPQ74GsgBqgAHmdlq4DHgGXfPK+S1VYDJwM7Rfka7e18z2xcYAdQBpgKXuPv6kjskEdmmb7+FLl1gwoTQPvVUGDIEGjSINy4pFdv6WL8LeA7Y391PdPeL3f1cd28GnAHUBC7ZymvXAR3cvTnQAjjJzA4HBgCD3P0A4Ffg8hI4DhFJxMaNMGhQKKo2YQLUrQsvvACvv66kn0aKTPzufiHwAdC2kMeWuPv97v7MVl7r7r4qalaObg50AEZH258hrLsrIsk2Zw60axeWQlyzBi64AObOhQsvDBdmSdrY5kBeNIzz8Pa8uZlVNLPpwBJgPPANsMzdc6OnLAT23J73FpEErV8Pd9wBLVvCJ5/AnnvCmDHwr39BRkbc0UkMEv0F510zO8eseN0Cd9/o7i2ABsChwCGJvtbMuppZtpll5+TkFGe3IpLvs8/CXPy+fUNVzSuuCD3/00+POzKJUaKJ/wpgFLDOzFaY2UozW5HoTtx9GTCRMGRUy8zyf1RuACzaymuGunuWu2dlqFciUjyrV8NNN4UrbWfNgv33D2P6Q4ZAzZpxRycxSyjxu3t1d6/g7ju5e42oXaOo15hZhpnViu7vAhwPzCN8AJwbPe0y4LXtjl5EtjRpEjRvHqZmQpiuOXMmHHtsrGFJ2ZHwBVxmVhs4kDCdEwB3n1zES+oDz5hZRcIHzIvuPtbM5gIjzOwu4HNg2HZFLiK/t3w59OoV1rqFcMXtsGHQpk28cUmZk+gFXJ2BHoShmenA4cBHhBk6hXL3mUDLQrZ/SxjvF5GS8sYbYfx+0SKoXBluvTVU1txpp7gjkzIo0TH+HkAbYIG7H0tI6MuSFZSIJCgnBy66CE47LST9ww6Dzz+Hv/1NSV+2KtHEv9bd1wKY2c7u/gVwcPLCEpEiuYfpmI0bhwuwdtkF7rsPPvgglFMWKUKiY/wLox9qXwXGm9mvwIJkBSUiRVi4EK66CsaODe0OHcISiPvtF29ckjISSvzuflZ0t5+ZTSSUang7aVGJyJby8kLVzJtughUroEYNuPdeuPxyXXkrxVKcWT2tgCMJZRc+UGE1kVI0f34oqjZpUmifcUZYHWtPXfguxZfQGL+Z/Y1QV6cOUBd4ysxuTWZgIgLk5ob5+E2bhqSfkQEjRsCrryrpy3ZLtMd/EdC8wA+8/yBM67wrSXGJyKxZYRjns89C++KLQ2XNunXjjUtSXqKzen6gwIVbhBr7hZZaEJEdtG5dqK3TqlVI+g0ahB9yn3tOSV9KxLYWYnmQMKa/HJhjZuOj9vHAp8kPTyTNfPJJ6OXPmRPaV10F//hH+CFXpIRsa6gnO/pzKvBKge2TkhKNSLr63//gttvg/vvDHP0DDggzeI4+Ou7IpBwqMvFvbZEVESlBEyaEGTvffhvWur3pJujXL1yUJZIERY7xm9nrZna6mVUu5LH9zOwOM+uUvPBEyrFly0LC/+MfQ9Jv1iwM9QwYoKQvSbWtoZ4uwPXA/Wb2C78ttr4vMB94yN1VVlmkuMaMCeP3P/wQaur87W9w882hwJpIkm1rqOdH4GbgZjNrSCi1vAb4yt1XJz88kXJmyRLo3h1Gjgzttm1D6eRGjeKNS9JKwlfuuvt3wHdJi0SkPHMPxdR69IClS6FqVfj736FbN6hYMe7oJM0knPhFZDt9/z1ceSW8+WZoH388DB0KDRvGGpakr0Qv4Co2M9vLzCaa2Vwzm2NmPaLt/cxskZlNj26nJCsGkVjl5YU1bjMzQ9KvVQueegrGjVPSl1gVp0jbLsDe7v5lgi/JBW5w92lmVh2YGl0ABjDI3QcWM1aR1PHVV2HGzuRoddKzzoKHH4b69eONS4TEi7SdTqjN83bUbmFmY4p6jbsvdvdp0f2VhIXWVVVKyrfcXPjnP8Ni55MnQ716MHo0vPyykr6UGYkO9fQjrJO7DMDdpxOmdCYkmhHUEvgk2nSNmc00syejRdwLe01XM8s2s+ycnJxEdyUSnxkzwtKHvXrB2rVw2WUwdy6cc07ckYn8TqKJf4O7L99smyfyQjOrBrwE9HT3FcCjwP5AC2AxcG9hr3P3oe6e5e5ZGRkZCYYpEoN168Li5llZMG0a7L03vP02PP007Lpr3NGJbCHRMf45ZvYXoKKZHQh0Bz7c1ouiK35fAp5395cB3P2nAo8/DowtdtQiZcWHH4aial98EdpXXx2maVavHm9cIkVItMd/LZAJrANeIFTr7FnUC8zMgGHAPHe/r8D2ggOdZwGzixGvSNmwalWYk3/kkSHpH3wwTJkCDz2kpC9l3jZ7/GZWEXjD3Y8FbinGe7cDLgFmmdn0aNv/AReaWQvCUNF3wBXFeE+R+I0fD127wnffhYuvbr45lFyoUmWbLxUpC7aZ+N19o5nlmVnNQsb5i3rd+0BhK0C/WZwARcqMX3+FG24Ic/EBWrQI5RZatYo1LJHiSnSMfxWh5z4e+F/+RnfvnpSoRMqaV14J5RV+/DEUVevXD268UUXVJCUlmvhfjm4i6eXHH+Haa8NcfIB27cICKYccEm9cIjsgocSvBVkk7biHNW579gxDPFWrhiUQu3ULi6WIpLCEEr+Z/YdC5u27+34lHpFI3BYsCEXV3n47tE88ER57DPbZJ964REpIokM9WQXuVwHOA3RlipQveXnw6KPQu3eYrlm7dlgD95JLwAqbpyCSmhId6lm62ab7zWwq8LeSD0kkBl9+CZ07w/vvh/a554Y5+fXqxRuXSBIkOtRTcL5aBcI3ANXyl9S3YQPce2+YpbNuHey+e6iiefbZcUcmkjSJJu+C9XRyCRde/bnEoxEpTdOnQ6dO8Pnnof3Xv4YPgdqF1g0UKTcSHeo5NtmBiJSatWvhzjthwADYuDEsijJ0aFgZSyQNJFqPv4eZ1bDgCTObZmYnJDs4kRL3wQfhitu77w4/5nbvDrNmKelLWkl0QnKnqKTyCUAdQg2efyQtKpGStnJluBCrffvwQ+4hh4Qfch94AKpVizs6kVKVaOLPn8t2CvCsu8+h8Do8ImXPuHHQpEmYpVOxItxySxjXP+KIuCMTiUWiP+5ONbN/E1bd6hOtoZuXvLBESsAvv8B118Gzz4Z2q1ahqFqLFrGGJRK3RBP/5YQVs75199Vmtivw16RFJbKjRo8Oi6IsWQI77wy33x4qa1bSLGSRRP8XtAWmu/v/zOxioBXwQPLCEtlOixfDNdeExc0hjOk/8QQcdFC8cYmUIYmO8T8KrDaz5sANwDfAs0W9wMz2MrOJZjbXzOaYWY9o+65mNt7Mvo7+1KRp2XHuoU5+48Yh6VerBo88ApMmKemLbCbRxJ/r7g78CXjI3R8GtrW+XC5wg7s3Bg4HrjazxkBv4F13PxB4N2qLbL///CcUUuvUCZYtg5NPhjlz4KqrVElTpBCJ/q9YaWZ9CNM43zCzCkCRK1C4+2J3nxbdXwnMA/YkfHjkl3l+BjhzO+IWCRdfDR4cZuyMHw+77hp+yH3jDdh777ijEymzEk385xMWWu/k7j8CDYB7Et2JmTUEWgKfAPXcfXH00I9AoVWwzKyrmWWbWXZOTk6iu5J0MW9eGL/v0QNWr4bzzw/bVElTZJsSSvxRsn8J2Dna9DPwSiKvNbNq0Wt7RheBFXxfp5A6/9FjQ909y92zMjIyEtmVpIMNG6B//zAl86OPoH59ePVVGDECdtst7uhEUkKiJRu6AKOBx6JNewKvJvC6yoSk/7y75y/d+JOZ1Y8erw8sKWbMkq6mToWsLLj1Vli/Hrp0gblz4U9/ijsykZSS6FDP1UA7YAWAu38NFNm9MjMDhgHz3P2+Ag+NAS6L7l8GvFacgCUNrVkTFkc57DCYORP22w/eeScUVqtVK+7oRFJOovP417n7eovGTs2sElsZoimgHeHH4FlmNj3a9n+EGj8vmtnlwAJU3lmKMnlyWCDl66/DDJ3rr4c77ghr4IrIdkk08b9nZv8H7GJmxwPdgNeLeoG7v8/W6/n8MfEQJS2tWAF9+oS5+BDm5z/5ZOj1i8gOSXSopxeQA8wCrgDeBG5NVlCS5t58M0zRfOSRUGKhb1+YNk1JX6SEbLPHb2YVgTnufgjwePJDkrT188+hqNrw4aGdlRV6+U2bxhuXSDmzzR6/u28EvjQzXREjyeEOI0eG4Zzhw6FKFRg4MEzXVNIXKXGJjvHXBuaY2afA//I3uvsZSYlK0scPP0C3bvBaNLnr6KNDUbUDDog3LpFyLNHEf1tSo5D04x5q4994IyxfDtWrh15+586qryOSZEUmfjOrAlwJHED4YXeYu+eWRmBSjn37bbj4asKE0D7tNHj0UWjQIN64RNLEtrpWzwBZhKR/MnBv0iOS8mvjRhg0KMzYmTAB6taFF16AMWOU9EVK0baGehq7e1MAMxsGfJr8kKRcmjMHLr8cPvkktP/yF7j/flAdJpFSt60e/4b8Oxrike2yfn240rZly5D099wTXn8dnn9eSV8kJtvq8Tc3s/yKmka4cndFdN/dvUZSo5PU9tlnoZc/a1ZoX3EFDBgANWvGG5dImisy8bt7xdIKRMqR1avD1bb33Qd5ebD//mGK5jHHxB2ZiJB4yQaRxEyaBM2bh6mZEKZrzpyppC9ShiQ6j1+kaMuXQ69e8Fi0ZEPTpmGefps28cYlIltQj1923BtvQGZmSPqVK0O/fpCdraQvUkapxy/bLycHevYMc/EhVM8cNix8CIhImZW0Hr+ZPWlmS8xsdoFt/cxskZlNj26nJGv/kkTu8K9/haJqL7wAu+wSfsj94AMlfZEUkMyhnqeBkwrZPsjdW0S3N5O4f0mGhQvhjDPCBVg//wwdOsDs2aGcckVNAhNJBUlL/O4+GfglWe8vpSwvL6xxm5kJY8dCjRrw+ONh7dv99os7OhEphjh+3L3GzGZGQ0G1t/YkM+tqZtlmlp2Tk1Oa8cnm5s+HP/4xXIC1YkXo8c+dGypp2tZW1xSRsqq0E/+jwP5AC2AxRRR9c/eh7p7l7lkZurQ/Hrm5YT5+06Zhfn5GRlgw5dVXQ+kFEUlJpTqrx91/yr9vZo8DY0tz/1IMM2eGcgvZ2aF98cWhqFqdOrGGJSI7rlR7/GZWv0DzLGD21p4rMVm3LpRbaN06JP299grz9J97TklfpJxIWo/fzP4FHAPUNbOFQF/gGDNrATjwHXBFsvYv2+Hjj0Mvf+7c0O7WDf7+9/BDroiUG0lL/O5+YSGbhyVrf7ID/vc/uO22MJTjDgceGIqqHXVU3JGJSBKoZEO6e/fd8OPtoEFhrdtevWDGDCV9kXJMJRvS1bJlcNNNoWcPoaLmsGFhbF9EyjX1+NPRa6+FcgtPPAE77QT9+4dFU5T0RdKCevzpZMkS6N49zMUHaNs29PIbNYo3LhEpVerxpwN3GD48JPiRI6FqVRg8GKZMUdIXSUPq8Zd3//0vXHklvPVWaB9/fKibv+++8cYlIrFRj7+8ysuDRx8NRdXeegtq1YKnnoJx45T0RdKcevzl0ddfhwJqkyeH9llnwcMPQ/36Rb9ORNKCevzlSW4u/POf0KxZSPr16sHo0fDyy0r6IrKJevzlxYwZ0KkTTJsW2pdeGi7K2nXXeOMSkTJHPf5Ut3Yt3HorZGWFpL/33vD22/DMM0r6IlIo9fhT2YcfhqJqX3wR2tdcA3ffDdWrxxuXiJRpSvypaNUquOUWePDBMEf/4IPDVbhHHhl3ZCKSAjTUk2rGjw9F1QYPDkXVeveG6dOV9EUkYerxp4pff4Ubbghz8QFatIAnn4SWLWMNS0RST9J6/NFi6kvMbHaBbbua2Xgz+zr6c6uLrUsBr7wSiqo99RTsvHNYHOXTT5X0RWS7JHOo52ngpM229QbedfcDgXejtmzNjz/CeefB2WeH+0ceGaZt9u4NlSvHHZ2IpKikJX53nwz8stnmPwHPRPefAc5M1v5TmnuYjtm4cbgAq1q18EPue++FH3JFRHZAaY/x13P3xdH9H4F6W3uimXUFugLsvffepRBaGbFgAVxxRaipA3DiiaGo2j77xBuXiJQbsc3qcXcnLLq+tceHunuWu2dlZGSUYmQxycuDhx4KRdXGjYPatUOv/623lPRFpESVdo//JzOr7+6Lzaw+sKSU9182fflluBDrgw9C+9xzw4dAva1+IRIR2W6l3eMfA1wW3b8MeK2U91+2bNgQZug0bx6S/u67w0svwahRSvoikjRJ6/Gb2b+AY4C6ZrYQ6Av8A3jRzC4HFgB/Ttb+y7xp00Ivf/r00O7UCQYODEM8IiJJlLTE7+4XbuWhPyZrnylh7Vq4/Xa45x7YuBEaNoTHH4fjjos7MhFJE7pytzS9/37o5X/1FZhBjx5w111huqaISClR4i8NK1dCnz5hFSwIC5wPGwZt28Ybl4ikJRVpS7a334YmTULSr1QJbrsNPv9cSV9EYqMef7IsXQrXXw/PPhvarVuHXn7z5vHGJSJpTz3+kuYeyiw0bhySfpUqYR3cjz9W0heRMkE9/pK0eDFcfXWopgnQvn1YIOWgg+KNS0SkAPX4S4J7qI3fqFFI+tWrw6OPwqRJSvoiUuaox7+j/vMf6NoV3nkntE85BYYMgb32ijcuEZGtUI9/e23cCA88EGbsvPMO1KkDw4fD2LFK+iJSpqnHvz3mzoXOneGjj0L7/PPDGri77RZvXCIiCVCPvzg2bAhX2rZsGZL+HnvAa6/BiBFK+iKSMtTjT9TUqaGQ2syZod2lS5imWatWrGGJiBSXevzbsmYN9OoFhx4akv5++8G778LQoUr6IpKS1OMvynvvhZ79119DhQpw3XVw551QtWrckYmIbDcl/sKsWBF6+UOGhHZmZii3cNhh8cYlIlICNNSzuTffDIl+yBCoXBn69g2Lpijpi0g5EUuP38y+A1YCG4Fcd8+KI47f+fln6NkTnn8+tNu0Cb38pk1jDUtEpKTFOdRzrLv/HOP+A3d48UW49lrIyYFddgnj+D17QsWKcUcnIlLi0nuMf9Ei6NYNxowJ7WOOCcsgHnBArGGJiCRTXGP8DvzbzKaaWdfCnmBmXc0s28yyc3JySnjvHhJ848Yh6deoAY89FqZpKumLSDkXV4//SHdfZGa7AePN7At3n1zwCe4+FBgKkJWV5SW252++CUXVJkwI7dNOC5U0GzQosV2IiJRlsfT43X1R9OcS4BXg0KTvdONGuO++8GPthAlQty688ELo8Svpi0gaKfXEb2ZVzax6/n3gBGB2Unc6ezYccQTccEO4EvcvfwmF1i68EMySumsRkbImjqGeesArFhJuJeAFd387KXtavx7+/nfo3z8UWNtzzzCsc/rpSdmdiEgqKPXE7+7fAslffHbFCmjXLvT2Aa64AgYMgJo1k75rEZGyrPxO56xRI4znr1kT1r095pi4IxIRKRPKb+IHeOQR2Gkn+MMf4o5ERKTMKN+JX2WTRUS2oCJtIiJpRolfRCTNKPGLiKQZJX4RkTSjxC8ikmaU+EVE0owSv4hImjH3kqt4nCxmlgMsiDuOBNUF4l9ZLDnK87FB+T4+HVvq2pHj28fdMzbfmBKJP5WYWXaZWEM4CcrzsUH5Pj4dW+pKxvFpqEdEJM0o8YuIpBkl/pI3NO4Akqg8HxuU7+PTsaWuEj8+jfGLiKQZ9fhFRNKMEr+ISJpR4i8GM9vLzCaa2Vwzm2NmPaLtu5rZeDP7OvqzdrTdzGywmc03s5lm1ireI9g2M6toZp+b2diova+ZfRIdw0gz2ynavnPUnh893jDWwBNgZrXMbLSZfWFm88ysbXk5d2Z2XfRvcraZ/cvMqqTyuTOzJ81siZnNLrCt2OfKzC6Lnv+1mV0Wx7FsbivHdk/073Kmmb1iZrUKPNYnOrYvzezEAttPirbNN7PexQrC3XVL8AbUB1pF96sDXwGNgX8CvaPtvYEB0f1TgLcAAw4HPon7GBI4xuuBF4CxUftF4ILo/hDgquh+N2BIdP8CYGTcsSdwbM8AnaP7OwG1ysO5A/YE/gPsUuCcdUzlcwccBbQCZhfYVqxzBewKfBv9WTu6X7uMHtsJQKXo/oACx9YYmAHsDOwLfANUjG7fAPtF/5ZnAI0TjiHuv4RUvgGvAccDXwL1o231gS+j+48BFxZ4/qbnlcUb0AB4F+gAjI3+I/1c4B9kW2BcdH8c0Da6Xyl6nsV9DEUcW80oOdpm21P+3EWJ//sowVWKzt2JqX7ugIabJcdinSvgQuCxAtt/97yydGybPXYW8Hx0vw/Qp8Bj46Jzuel8Fva8bd001LOdoq/HLYFPgHruvjh66EegXnQ//z9kvoXRtrLqfuBmIC9q1wGWuXtu1C4Y/6Zjix5fHj2/rNoXyAGeioaynjCzqpSDc+fui4CBwH+BxYRzMZXyc+7yFfdcpcw53EwnwjcYSNKxKfFvBzOrBrwE9HT3FQUf8/Dxm3JzZM3sNGCJu0+NO5YkqUT4ev2ou7cE/kcYLtgkhc9dbeBPhA+3PYCqwEmxBpVkqXqutsXMbgFygeeTuR8l/mIys8qEpP+8u78cbf7JzOpHj9cHlkTbFwF7FXh5g2hbWdQOOMPMvgNGEIZ7HgBqmVml6DkF4990bNHjNYGlpRlwMS0EFrr7J1F7NOGDoDycu+OA/7h7jrtvAF4mnM/ycu7yFfdcpdI5xMw6AqcBF0UfbJCkY1PiLwYzM2AYMM/d7yvw0Bggf8bAZYSx//ztl0azDg4Hlhf4qlqmuHsfd2/g7g0JP/hNcPeLgInAudHTNj+2/GM+N3p+me2BufuPwPdmdnC06Y/AXMrBuSMM8RxuZn+I/o3mH1u5OHcFFPdcjQNOMLPa0beiE6JtZY6ZnUQYZj3D3VcXeGgMcEE0E2tf4EDgU+Az4MBo5tZOhP+zYxLeYdw/cqTSDTiS8PVyJjA9up1CGB99F/gaeAfYNXq+AQ8Tfn2fBWTFfQwJHucx/DarZ7/oH9p8YBSwc7S9StSeHz2+X9xxJ3BcLYDs6Py9SpjpUS7OHXA78AUwG3iOMAskZc8d8C/C7xUbCN/WLt+ec0UYL58f3f4a93EVcWzzCWP2+XllSIHn3xId25fAyQW2n0KYWfgNcEtxYlDJBhGRNKOhHhGRNKPELyKSZpT4RUTSjBK/iEiaUeIXEUkzSvxSpplZHTObHt1+NLNFBdo7xR1fQWZ2jJkdkcT338XM3rNQQbXhZtUdu5jZ1GjO+kAz65CsOCT1Vdr2U0Ti4+5LCfPvMbN+wCp3HxhXPGZWyX+rf7O5Y4BVwIcl9H6b6wS87O4bw3Vam97jEuBaoIO7/2pmDwKPAxMSjUPSi3r8knLMrHXU851qZuMKXMY/ycwGmVm2hXr7bczs5agW+13RcxpGdc+fj54z2sz+kMD73m9m2UAPMzvdQh37z83sHTOrFxXtuxK4Lvo20t7MnjazcwvEvSr68xgzm2JmY4C5UQ/+HjP7zEI99iu2cugX8dvVqvnv+WdCzaET3P1nAHdfANQxs91L6u9cyhclfkk1BjwInOvurYEngf4FHl/v7lmE+vOvAVcDTYCOZpZfgfJg4BF3bwSsALpFNZiKet+d3D3L3e8F3gcO91DsbQRws7t/F+1zkLu3cPcp2ziOVkAPdz+IcOXmcndvA7QBukSX5/920GFYa79oP/n2AR4iJP0fN3v/aYR6PSJb0FCPpJqdCYl8fDTcUZFw+Xu+/Hols4A5HtXXMbNvCUWtlgHfu/sH0fOGA92Bt7fxviML3G8AjIy+EexEqPNfXJ+6e/7rTgCaFfh2UJNQk6Xg+9aNYi8oB/gF+DMwaLPHlhAqdYpsQYlfUo0REnrbrTy+Lvozr8D9/Hb+v/fN65R4Au/7vwL3HwTuc/cxZnYM0G8rr8kl+lZtZhUIHxKFvZ8B17p7UQXE1hBq7BS0mlCvZYqZLXH3gqV8q0SvEdmChnok1awDMsysLYQy2WaWWcz32Dv/9cBfCEM3XxbjfWvyWwncguu4riQsyZnvO6B1dP8MoPJW3m8ccFU03ISZHWRhkZhN3P1XoKKZVdls+xJC7f27rcB6rMBBhIJtIltQ4pdUk0coJTzAzGYQKhkWdwrll8DVZjaPUKHzUXdfX4z37QeMMrOphGUL870OnJX/4y5hZs3R0fu15fe9/IKeIJRRnhZN0XyMwr+N/5tQIfZ3oiGjM4AnzezQ6APkAEIlUpEtqDqnpJVo9s1Yd28SdyzFZWatgOvc/ZJtPO8soJW731Y6kUmqUY9fJEW4+zRgoplV3MZTKwH3lkJIkqLU4xcRSTPq8YuIpBklfhGRNKPELyKSZpT4RUTSjBK/iEia+X9sehvt2xFRrwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# part A\n",
    "Rg = 8.314  # J/mol-K\n",
    "Vm = 0.00224 # m^3/mol\n",
    "\n",
    "def P(Vm,T):\n",
    "    # return pressure in bar\n",
    "    return 1e-5 * Rg * T / Vm\n",
    "T = np.linspace(100,1200)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure()\n",
    "plt.plot(T,P(Vm,T),'r-',linewidth=2,label='Ideal Gas Pressure')\n",
    "plt.xlabel('Temperature (K)')\n",
    "plt.ylabel('Pressure (bar)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part B\n",
    "\n",
    "Repeat part A but use the non-ideal gas Van der Waals Equation of State:\n",
    "\n",
    "$P+\\frac{a}{V_m^2}=\\frac{R_g T}{V_m-b}$\n",
    "\n",
    "with constants $a$ and $b$ for ethanol with critical properties $T_c=514 K$ and $P_c=6.137\\mathrm{x}10^6 Pa$:\n",
    "\n",
    "$a = 25 \\frac{\\left(R_g T_c\\right)^2}{64 P_c} = 25 \\frac{\\left(8.314 \\frac{J}{mol\\,K} 514 K\\right)^2}{64 \\,\\mathrm{x}\\,6.137\\mathrm{x}10^6 Pa}$ = $1.1623 \\frac{Pa\\,mol^2}{m^6}$\n",
    "\n",
    "$b = \\frac{R_g T_c}{8 P_c} = \\frac{8.314 \\frac{J}{mol\\,K} 514 K}{8 \\,\\mathrm{x}\\,6.137\\mathrm{x}10^6 Pa}$ = $8.70416\\mathrm{x}10^{-5} \\frac{m^3}{mol}$\n",
    "\n",
    "Compare the ideal gas and non-ideal gas results on the same plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9J0lEQVR4nO3deZyN5fvA8c+FYYrsQilLi20Y21gjISotSipkq6isRQstXyq+KUqR0KZEX0qlxC+VLIUWo5E9FCLZyr7Ncv/+uJ5ZMMsZ5syZM+d6v17n5TzPeZ5z7jNnXHOf+7nu6xbnHMYYY0JHnkA3wBhjTPaywG+MMSHGAr8xxoQYC/zGGBNiLPAbY0yIyRfoBviiZMmSrkKFCoFuhjHGBJXo6Oi9zrlSp+8PisBfoUIFli9fHuhmGGNMUBGRrantt6EeY4wJMRb4jTEmxFjgN8aYEBMUY/ypiY2NZfv27Rw/fjzQTTG5THh4OOXKlSMsLCzQTTHGL4I28G/fvp0LLriAChUqICKBbo7JJZxz7Nu3j+3bt1OxYsVAN8cYvwjaoZ7jx49TokQJC/omS4kIJUqUsG+SJlcL2sAPWNA3fmG/VybHOHnSL08b1IHfGGNyJefg3XehUiXYsCHLn94C/zkoVKhQqvu7d+/OzJkzz+o5hw0bxujRo1N9bOrUqdSsWZPq1asTGRnJfffdx/79+8/qdRItXLiQIkWKUKtWLapWrcozzzxzTs9njDlHa9dC8+bQowfs2AGTJ2f5S1jgDxJffvklY8aM4f/+7/9Ys2YNK1asoHHjxuzateucn7tp06bExMSwfPlypk6dyooVK055PC4u7pxfw1fZ+VrG5ChHj8KQIRAZCYsXc6DkZbj3p8Lzz2f5S1ngzwLOOfr27UvlypVp1aoVu3fvTnosOjqaq6++mrp169KmTRt27twJwJtvvklUVBSRkZG0b9+eo0ePpvsaI0aMYPTo0Vx88cUA5M2bl3vuuYfKlSsD8OyzzxIVFUVERAS9evUicWW1sWPHUq1aNWrWrMldd92V7msULFiQunXrsmnTJoYNG0aXLl1o0qQJXbp0Yc+ePbRv356oqCiioqJYsmQJAIsWLaJWrVrUqlWL2rVrc+jQIXbu3EmzZs2oVasWERERfPfdd8Cp35BmzpxJ9+7dAf2G9MADD9CgQQMee+wxNm/ezHXXXUfdunVp2rQp69ev9/WjMCY4zZkD1avDyJEQH8+x+/pRPew3ZhboDP645uScy/G3unXrutOtXbs2eUNHxLL+loGCBQs655z7+OOPXatWrVxcXJzbsWOHK1KkiPvoo4/cyZMnXaNGjdzu3budc85Nnz7d9ejRwznn3N69e5Oe58knn3Rjx451zjk3dOhQN2rUqDNeq1ixYm7//v1ptmXfvn1J9++++273+eefO+ecK1u2rDt+/Lhzzrl///33jPMWLFjg2rZtm9Sm8uXLu9WrV7uhQ4e6OnXquKNHjzrnnOvYsaP77rvvnHPObd261VWpUsU559yNN97ovv/+e+ecc4cOHXKxsbFu9OjRbvjw4c455+Li4tzBgwdP+Xk559xHH33kunXr5pxzrlu3bq5t27YuLi7OOedcixYt3G+//eacc+6HH35w11xzTZrv219O+f0yxl/+/NO5225LjjmRkc4tW+acc27YMOc6djy3pweWu1RiatDm8eckixcvpmPHjuTNm5eLLrqIFi1aALBhwwZWr17NtddeC0B8fDxly5YFYPXq1Tz11FPs37+fw4cP06ZNG59fb9WqVXTp0oVDhw7x3//+lzvvvJMFCxbw4osvcvToUf755x+qV6/OTTfdRM2aNencuTPt2rWjXbt2qT7fd999R+3atcmTJw+DBw+mevXqfPTRR9x8882cd955AHzzzTesXbs26ZyDBw9y+PBhmjRpwsCBA+ncuTO33XYb5cqVIyoqinvuuYfY2FjatWtHrVq1MnxPHTp0IG/evBw+fJilS5fSoUOHpMdOnDjh88/GmKAQFwfjxsF//gOHDxN//gW83uYzavdvylUNNSw/+STkzeufl88dgT+HLhjvnKN69eosW7bsjMe6d+/OrFmziIyM5N1332XhwoXpPlf16tVZsWIF11xzDTVq1CAmJoa+ffty7Ngxjh8/Tu/evVm+fDmXXHIJw4YNS8pDnzNnDosXL2b27NmMGDGCVatWkS/fqR9706ZN+eKLL854zYIFCybdT0hI4IcffiA8PPyUYwYPHkzbtm2ZO3cuTZo0Yd68eTRr1ozFixczZ84cunfvzsCBA+natespaZKn58knvlZCQgJFixYlJiYm3Z+HMUHrhx/ggQdg5Urdvu023qj7Jv2fLE7ltbBqFYSFQT4/Rmcb488CzZo1Y8aMGcTHx7Nz504WLFgAQOXKldmzZ09S4I+NjWXNmjUAHDp0iLJlyxIbG8u0adMyfI0hQ4bwyCOPsH379qR9x44dA5KDaMmSJTl8+HBSRlFCQgJ//vkn11xzDS+88AIHDhzg8OHDZ/UeW7duzbhx45K2EwPz5s2bqVGjBo8//jhRUVGsX7+erVu3Urp0aXr27Ml9992XdLG4dOnSrFu3joSEBD799NNUX6dw4cJUrFiRjz76CNA/nisT/4MYE8z+/VcDfuPGGvTLl4fZs+Hjj+kxsDgtWsALL2jQ97fc0eMPsFtvvZVvv/2WatWqcemll9KoUSMA8ufPz8yZM+nfvz8HDhwgLi6Ohx56iOrVq/Pcc8/RoEEDSpUqRYMGDTh06FC6r3HDDTewZ88err/+euLj4ylatCgRERG0adOGokWL0rNnTyIiIihTpgxRUVGADi3dfffdHDhwAOcc/fv3p2jRomf1HseOHUufPn2oWbMmcXFxNGvWjIkTJ/LKK6+wYMEC8uTJQ/Xq1bn++uuZPn06o0aNIiwsjEKFCjFlyhQARo4cyY033kipUqWoV69emn+Epk2bxoMPPsjw4cOJjY3lrrvuIjIy8qzabUzAOQfTpsGgQbB7N+TLx9K7xvLSoZ58cG0+CgDh4TB/fvY1SVwOHSZJqV69eu70hVjWrVtH1apVA9Qik9vZ75fJEhs2wIMPgjcKwFVXET9+IhF3Vmf9ehgzBh56yH8vLyLRzrl6p++3Hr8xxmS1Y8c0//6FF7TsQokSuFGjke7dyCvCpEkwb56O/ASCBX5jjMlKX30FvXvD5s0A7O70EAOOjaT6jgI85eU3NGumt0Cxi7vGGJMV/voL7roL2rTRoB8RAd9/z5r7xjD90wK8/DJkcCkv21jgN8aYcxEfD6+9BlWrwowZcP75HHnuZVixApo04Zpr4NVXIToaLrgg0I1VFviNMeZsLV8ODRpAv35w8CDuxpsYM/BPLh3zML/9kZyX2b8/5KR1fSzwG2NMZh04oMG+fn3tyl9yCcyahcz+nJhtxfnnH/j440A3Mm0W+M+BiDBo0KCk7dGjRzNs2LAsee60yjNv2bKFiIiIs37etEpJ79q1i06dOlGpUiXq1q1Lo0aN0pxklRnNmzencuXKREZG0qRJEzb4oba4MdnGOZg+HapU0eGdPHk4/tBgdi1cB7fcAsDLL2vNtSFDAtzWdFjgPwcFChTgk08+Ye/evYFuyjlxztGuXTuaNWvG77//TnR0NNOnTz9llvC5mDZtGitXrqRbt248+uijZzweHx+fJa+TEeccCQkJ2fJaJhfatEkv3HbsCH//DY0asXrGGiLnPk+nngWTKseUKAE33BDYpmbEAv85yJcvH7169WLMmDFnPLZlyxZatGhBzZo1admyJdu2bQO0Rk///v1p3LgxlSpV8mnBlujoaCIjI4mMjGT8+PFJ++Pj43n00UeJioqiZs2aTJo0CYDDhw/TsmVL6tSpQ40aNfjss8/Sff5vv/2W/Pnz80CKpOLy5cvTr1+/pPfStGlT6tSpQ506dVi6dClAmuWX09KsWTM2bdoE6DePQYMGERkZybJly5g6dSr169enVq1a3H///cTHxxMfH0/37t2JiIigRo0aST/n1EpNn/4NKSIigi1btrBlyxYqV65M165diYiI4M8//2TUqFFJP7OhQ4dm+PM3Ie7ECXj2Wc3S+fprKFYM3ngDvv+e0s0qs28f7Nypk3KDhd8Dv4jkFZFfROQLb7uiiPwoIptEZIaI5M+a1zmzbPVNN+m+2bOT973xhu7r1St5319/6b6LLsr86/bp04dp06Zx4MCBU/b369ePbt268euvv9K5c2f69++f9NjOnTv5/vvv+eKLLxg8eHCGr9GjRw/GjRt3Rs2at99+myJFivDzzz/z888/8+abb/LHH38QHh7Op59+yooVK1iwYAGDBg0ivRnaa9asoU6dOmk+fuGFF/L111+zYsUKZsyYkfRePvjgA9q0aUNMTAwrV67MsArn7NmzqVGjBgBHjhyhQYMGrFy5khIlSjBjxgyWLFlCTEwMefPmZdq0acTExLBjxw5Wr17NqlWr6NGjB6ClH3755Rd+/fVXJk6cmOHPb+PGjfTu3Zs1a9awYcMGNm7cyE8//URMTAzR0dEsXrw4w+cwIWr+fKhZE4YO1T8AXbvy3bubcff1hDx5KFVKD/nlFyhdOtCN9V129PgHAOtSbL8AjHHOXQ78C9ybDW3wm8KFC9O1a1fGjh17yv5ly5bRqVMnALp06cL333+f9Fi7du3IkycP1apVy3AFrf3797N//36aebM9unTpkvTYV199xZQpU6hVqxYNGjRg3759bNy4EeccTzzxBDVr1qRVq1bs2LEjUyt19enTh8jIyKSaP7GxsfTs2ZMaNWrQoUOHpPLMUVFRTJ48mWHDhrFq1SouSCNXrXPnztSqVYslS5Yk9crz5s1L+/btAZg/fz7R0dFERUVRq1Yt5s+fz++//06lSpX4/fff6devH19++SWFCxcGSCo1PXXq1DMqjaamfPnyNGzYMOln9tVXX1G7dm3q1KnD+vXr2bhxo88/GxMidu2Cu++GVq3gt990TP/bb7k//D2a3VKM999PPjQyEgoUCFxTz4ZfZ+6KSDmgLTACGChal7cF0Mk75D1gGDDhXF8rtQ5typ5+ol69Tu3tg/b0z6Vk0UMPPUSdOnWSeqQZKZDitySxJ/7kk08yZ84cAJ9LEjvnGDdu3Bm1/N9991327NlDdHQ0YWFhVKhQ4YwyyClVr16dj1OkIIwfP569e/dSr56W+BgzZgylS5dm5cqVJCQkJJVmTqv88ummTZuW9FyJwsPDyesVG3fO0a1bN55PZYm5lStXMm/ePCZOnMiHH37IO++8k2ap6ZTj9ynfb8ry0s45hgwZwv3335/mz8OEsIQEHRYYPFgzd8LD4emn4ZFHIH9+Gm/TemuxsYFu6Lnxd4//FeAxIPF/ZAlgv3MucWHV7cDFqZ0oIr1EZLmILN+zZ4+fm3luihcvzh133MHbb7+dtK9x48ZMnz4d0MDXtGnTdJ9jxIgRxMTEnBH0ixYtStGiRZO+MaQs4dymTRsmTJhArPdb+Ntvv3HkyBEOHDjAhRdeSFhYGAsWLGDr1q3pvnaLFi04fvw4EyYk//1NuRTkgQMHKFu2LHny5OH9999PuhibVvnlzGrZsiUzZ85MWrLyn3/+YevWrezdu5eEhATat2/P8OHDWbFiRZqlpitUqJD0+itWrOCPP/5I9bXatGnDO++8k1QZdMeOHacslWlCWEyMlkx+8EEN+tddx6Y5G/im/hOQX0eku3aFjRvh3qAep/Bjj19EbgR2O+eiRaR5Zs93zr0BvAFanTNrW5f1Bg0axGuvvZa0PW7cOHr06MGoUaMoVaoUkydPPuvnnjx5Mvfccw8iQuvWrZP233fffWzZsoU6dergnKNUqVLMmjWLzp07c9NNN1GjRg3q1atHlSpV0n1+EWHWrFk8/PDDvPjii5QqVYqCBQvywgsvANC7d2/at2/PlClTuO6665J60AsXLky1/HJmVatWjeHDh9O6dWsSEhIICwtj/PjxnHfeefTo0SOpJ//888+nWWo6sX3Vq1enQYMGXHnllam+VuvWrVm3bl1S6exChQoxdepULrzwwrNqu8kFDh3SlbDGjtUe/0UXwauvsurK9tRvIBQqBOvWQcmSei3QW0QvqPmtLLOIPA90AeKAcKAw8CnQBijjnIsTkUbAMOdcuusOWllmk93s9ysEOAeffAIDBsCOHZAnj07KevZZKFwY53SI3/s7QPHigW5w5qVVltlvQz3OuSHOuXLOuQrAXcC3zrnOwALgdu+wbkD6uYbGGJPV/vgDbrwRbr9dg35UFEe/i2ZokVfYF6tJBCI6Eev994Mz6KcnEHn8j6MXejehY/5vZ3C8McZkjZMntU5+9eowdy4UKQKvvw7LlvHAxFo8+6xex0102hLTuUa21ON3zi0EFnr3fwfqZ9HznrKAtzFZIRhWpTNnYfFiXflknZdd3qkTvPQSlCkDwFNPwfr1gVscJTsF7czd8PBw9u3bZ/9JTZZyzrFv376klFWTC+zZAz16wNVXa9C/4grcV18z89ZpDHm1TNJhV14JP/6oxTZzu6BdgatcuXJs376dnJ7qaYJPeHg45cqVC3QzzLlKSIDJk+Gxx+CffzQl84kn4PHH2b4nnM6X68jPrbdqkU04c/Z/bhW0gT8sLIyKOanAtTEm51i1SvPxlyzR7VatcONfR668AtAqys8/D+efD/XOyHnJ/YJ2qMcYY85w5Ag8/jjUqaNBv3Rp+OADNoz7imt6XUHKskwDB+p4fp4QjIIh+JaNMbnS559DtWrw4ou6HGKfPnq1tmNHps8QFi3SC7gmiId6jDEGgG3bdG3DxPLjtWvDpEmcqBmVVDxt8GAdz0+xblJIsx6/MSY4xcbC6NG6yPlnn+lK5q++yrFFPzFoehR162olZdDqmSNG5L6JWGfLAr8xJvgsWaLj+I8+CkePQocOmqrZvz8Slo85c3Rz0aJANzRnsqEeY0zw2LdPx23eeku3K1aE8ePZV/96LrgA8qOzbadMgbx5oW7dgLY2x7IevzEm53MO3n1XF0R56y0IC9MrtWvWMCfheqpWhZEjkw+vX9+Cfnqsx2+MydnWrtWc/MRczObNYcIE/SOA5uLv2aOjPwkJoZmemVn2IzLG5ExHj+pM28hIDfqlSsGUKSR88y3RR5LXmLjmGh3L/7//s6DvK/sxGWNynrlztYLm889DXJyul7p+Pcc7dKFpM6FJE9iwIfnwZs0s6GeG/aiMMTnH9u1aI79tW9iyBWrWhKVLYdIkKF6c8HCoXFnTMnfuDHRjg5cFfmNM4MXFwSuvaE7+xx9DwYJaMjk6mh+kEZs3Jx/68ss67N+8eaAaG/ws8BtjAuvHHyEqCh5+GA4f1nKZ69bBwIFMm5GPxo2hZ09N7AEoWlRv5uxZ4DfGBMb+/dC7NzRqBDExUL681tv55BMtnwm0aaPrpNSvr18KTNawwG+MyV7OwbRpOlg/YYLOtHr8cVizhj0Nb2LkyOTefcmSsHGj5uiHhQW22bmJ5fEbY7LPb79pL3/+fN2+6ioN/hERJCRA07qarVOmDHTvrocULBiw1uZa1uM3xvjf8eMwdCjUqKFBv3hxePttTcCPiAA0HfPJJ6FVK03PNP5jPX5jjH999ZXWxt+0Sbd79IAXXyS+WEleG6e9+zvv1IfuvltvobIEYqBY4DfG+MfOnbrM1fTpul29ug7rNG0KwOefwkMPQYkScP31ULiwBfzsYkM9xpisFR8P48drLZ3p0+G88/Tq7IoVSUEfoF076NoV3nlHg77JPhb4jTFZJzoaGjaEvn3h4EG48UadbfX443z/U36uvhr27tVDReC99+DmmwPb5FBkgd8Yc+4OHNDlD+vXh+XLoVw5+PRTzcuvUAHn9Nru4sUwalSgG2tsjN8Yc/acg48+0sH6nTs1J3/QIBg2DAoVIi4O8uXT3v2kSfD++zBkSKAbbSzwG2POzubNmq0zb55uN2wIEydCZCS7dkH/e6FQIc3aBLj8cnjmmcA11ySzoR5jTOacOAHPPadZOvPmaeGcSZN0JZTISECrMcyaBTNmwF9/BbKxJjXW4zfG+O7bb3XmbWIx/C5dYPRouPBC9u3T1EzQagzvv6+11y66KHDNNamzHr8xJmO7d2uQb9lSg37lyvpHYMoUuPBCXn0VLr0UFi5MPuWOO3QtdJPzWOA3xqQtIUGHcSpXhqlTITwchg+HlSt1zUPP/v26UmLicL/J2WyoxxiTupUr4YEH4IcfdLtNG52YddllHDsGOzbpBVvQTJ2rrtIvBCbnsx6/MeZUhw5pSmbduhr0y5aFDz/U1cwvu4xNm/Qa7o03au01gPz5LegHEwv8xhjlnE66qlZN1zd0DgYMgPXroUOHpEI6l1yi6fp589q6t8HKhnqMMfDHH9CvH8yZo9tRUZqTX6cOAHPnao++QAG9zZ0LF1+sPX0TfKzHb0woO3lSC6hVr65Bv3BheO01WLYsKegPGgRt2+phiSpWtKAfzCzwGxOqFi+G2rX1yuyxY9Cpk6Zq9umj4zieW27RvwdlywawrSZL+S3wi0i4iPwkIitFZI2IPOPtrygiP4rIJhGZISLWbzAmO+3dC/fcA1dfrZUzr7gCvv5a18EtU4aNG+GDD5IPb9YMtm2DXr0C12STtfzZ4z8BtHDORQK1gOtEpCHwAjDGOXc58C9wrx/bYIxJlJCghXMqV4bJk3WsZtgw+PVXXe8Q2LEDatbURbLWr08+tUiRwDTZ+IdPF3dFJA8QCVwEHANWO+d2p3eOc84Bh73NMO/mgBZAJ2//e8AwYEJmG26MyYTVqzUnf8kS3W7VSnPyr7zylMMuvhg6d4a4OChVKgDtNNki3cAvIpcBjwOtgI3AHiAcuFJEjgKTgPeccwlpnJ8XiAYuB8YDm4H9zrk475DtwMVpnNsL6AVw6aWXZu5dGWPUkSPw7LOanhkXB6VLw5gxcNddIMLRo1ox8957k/8GTJp0yhC/yYUy6vEPB14H7vd68ElE5EK0594F7bmfwTkXD9QSkaLAp0AVXxvmnHsDeAOgXr16LoPDjTGnmz1bUzS3btUc/N69YcQIrabpGTpUa6wtXw7z5+s+C/q5X7qB3znX0RvmaQQsPe2x3cArvryIc26/iCzwnqeoiOTzev3lgB1n03BjTBq2bdOJV7Nm6XatWtqNr1//jEMHD4aYGPjvf7OzgSbQMry46w3jjM/sE4tIKa+nj4icB1wLrAMWALd7h3UDPsvscxtjUhEbCy+9pDNvZ83SVVBeeQV+/hnq109aLKtTJ52UC1pG+euvdb6WCR2+ZvXMF5H2It6cbd+UBRaIyK/Az8DXzrkv0GsGA0VkE1ACeDtTLTbGnGnZMqhXDx55RMf1b79d03IGDNC1D9G1z3v3hv/9T5fCNaHL15IN9wMDgTgROQ4ImrhTOK0TnHO/ArVT2f87cOZ3TmNM5v3zj07AeuMN3a5YUWfe3nADoBmcInorUgRef11PuemmALbZBJxPPX7n3AXOuTzOufzOucLedppB3xjjZ87pIihVqmjQDwuDJ57QtE0v6K9fr5Ov3n8/+bQOHeD++yGPzdkPaT4XaRORYsAVaDonAM65xf5olDEmHevWwYMPwqJFun311TBhAlStesphP/2kafv798Pdd1uwN8l8ncB1HzAAzcKJARoCy9DJWMaY7HD0qKZjjhqlF3JLltSLuV26JJVMPnhQ6+qA7t67V2fhWtA3Kfn66zAAiAK2OueuQcfu9/urUcaY0/zf/0FEhOZdxsZCz55aUK1rVxDh5El4+GGdhLV3r54iAgMHQrFigW26yXl8DfzHnXPHAUSkgHNuPVDZf80yxgBaPKdDBx23/+MPqFFDx2/eeAOKF086LCxMS+7s3atroBuTHl/H+Ld7OfmzgK9F5F9gq78aZUzIi4vTWjpPPQWHD0PBglpboX9/jfJokM+bV3v0IvDmmzqe75XRNyZNPgV+59yt3t1h3gzcIsCXfmuVMaHsp5+0oNovv+h2u3bw6quQombVN99Ax4667u3kybqvUqXsb6oJTpnJ6qkDXIVW2FzinDvpt1YZE4r274cnn9QMHec00I8bBzfffMahl16qa6Jv2wYnTuhyiMb4yqcxfhH5D1qIrQRQEpgsIk/5s2HGhAzndDptlSo6wypvXnjsMV0kxQv68fHay0905ZX6xeCbbyzom8zztcffGYhMcYF3JJrWOdxP7TImNGzcqHUUEqN6kya6yHlERNIhCQlaPn/hQq2g2cJLoq5ZM/uba3IHX7N6/iLFxC2gAFZV05izd/y4XqyNiNCgX7y4ro61ePEpQR80B79lS7joIr3ma8y5ymghlnHomP4BYI2IfO1tXwv85P/mGZMLff219vI3bdLtHj3gxRd1QpZn6VIN+A0b6vZjj2lpfVsC0WSFjIZ6lnv/RqMLqSRa6JfWGJOb/f23zqj63/90u1o1vZDbrNkph82dq9k6V1wBK1dCeLguj5s/fwDabHKljBZiSXVlLWNMJsTH60IoTzwBBw7AeefBf/6jfwRSieYtW+o8Laugafwlo6Ge2ejyh18652JPe6wS0B3Y4px7x28tNCaYrVihOfk//6zbbdtqimbFikmH7N6t5XeGD9cMnQIFdClEb56WMVkuo6Genmgd/ldE5B+SF1uvCGwCXnPO2Qpaxpzu4EF4+mmtjZ+QAOXKwdixOhnrtPWMbrtNqzAULAjDhuk+C/rGnzIa6vkbeAx4TEQqoKtqHQN+c84d9X/zjAkyzsHMmbry1c6dmpM/cKBG9AsuSPWU55/XopvdumVvU03o8nnmrnNuC7DFby0xJtht3gx9+8KXXjWTBg00J79WraRD4uO1+kJcnGbqADRtmnyKMdnB58BvjEnDiRM6SD9ihObnFy0KI0dq6eTTCuEvXw6DBuk13U6ddATImOxmgd+Yc7Fgga6GtWGDbnfpon8ESpdOOsS55GH9Bg204GaDBhb0TeD4vC6PiJwnIlaD3xjQVJyuXbV+woYNULmy1lOYMuWUoL94sY70JP5dAHjuOc3TNyZQfC3SdhNam+dLb7uWiHzux3YZkzMlJOgiKJUr6yrmBQpoJF+5MrmITgrvvqsLpIwcmf1NNSYtvg71DAPq483Ydc7FiEjF9E4wJtdZuVKHdZYt0+3WrXWxlMsvP+Wwo0fh/PP1/ksv6QzcgQOzua3GpMPXoZ5Y59yB0/a5rG6MMTnS4cPwyCNQt64G/TJlYMYMTcVJEfT37NFVEtu21XF90NWxhgyx0skmZ/G1x79GRDoBeUXkCqA/sNR/zTImB3AOPvtMq6Nt364ZOv366dBOKtXS8uSBRYu0x79unZbiMSYn8rXH3w+oDpwAPkCrdT7kpzYZE3hbt8Itt8Ctt2rQr1sXfvxRZ9+mCPrbtiX37kuUgA8/hDVrLOibnC3DwC8ieYE5zrknnXNR3u2pxEVZjMlVYmO1RHK1ajB7NhQurGUXfvwR6tU75dAJE/Qa73spShk2bw7ly2dvk43JrAwDv3MuHkgQEasEbnK3776D2rXh8cd1vObOO2H9eujTR0svnKZQIZ2vtWJFANpqzDnwdYz/MLDKW4jlSOJO51x/v7TKmOy0d68G+3e8IrOXXaZr37ZufcphiWP3devq9t1369q3DRpkc3uNOUe+Bv5PvJsxuUdCgo7TPPoo7NundRSGDIHBg3X1kxR27tSaOvv3a/AvVUpn41rQN8HIp8BvC7KYXGfNGs3J/+473W7ZUnv5V16Z6uFlykCFCrBrl35BKFUq+5pqTFbzKfCLyB+kkrfvnKuU5S0yxp+OHNF0zJde0hKZF14IL7+sFdNS1Ml3Dj79VFdFLFlSH/rf/zShx5ZANMHO16GelOkM4UAHoHjWN8cYP/riCy2bvHWrRvIHH9SKmsWKnXHoc8/B0KFajicxa8d6+Sa38CmP3zm3L8Vth3PuFaCtf5tmTBb5809d5uqmmzTo16qlM3Bffz3VoA/6BaB0aWjYMDlP35jcwtehnjopNvOg3wCspLPJ2eLidMLVf/6jQzyFCmlXvm9fyHfqr++GDTBnTnJNncsvhy1bzrjGa0yu4GvwfinF/Th0Ja47srw1xmSVH37QRc5XrtTt9u3hlVdSLYJ/8KBm5xw4oF8GEotsWtA3uZWvWT3X+LshxmSJf//VlMw33tAxmooVdebtDTekeUrhwroM4qZNp6ySaEyu5Ws9/gEiUljUWyKyQkRaZ3DOJSKyQETWisgaERng7S8uIl+LyEbv39QHWY3JDOdg6lStoTBpkg7lPPEErF59RtA/ckSXP1ywIHnfkCE6f6u4pSyYEOBrkbZ7nHMHgdZACaALkNHSEnHAIOdcNaAh0EdEqgGDgfnOuSuA+d62MWdv/XrNw+/SRWsjX321DvGMGJFcGD+FSZM0g7N3b138HE7J5DQm1/M18Cf+t7gBmOKcW5NiX6qcczudcyu8+4eAdcDFwC1A4oSw94B2mWyzMerYMXj6aahZU7vvJUvqklcLFkDVqqccmjIzp29fuOsu/YKQSgkeY3I9XwN/tIh8hQb+eSJyAZDg64uISAWgNvAjUNo5t9N76G+gdBrn9BKR5SKyfM+ePb6+lAkVX34JEREwfLhW1OzZU1NzunU7YyLW9On6JeC4V082f36djJVYc8eYUONr4L8XHZKJcs4dBcKAHr6cKCKFgI+Bh7zhoiTOOUcaK3k5595wztVzztUrZTNnTKIdO+COO+D66+H336FGDViyRC/mpjJAHxsLzz6rlRmmTAlAe43JgXxN52wExDjnjojI3UAd4NWMThKRMDToT3POJRZ52yUiZZ1zO0WkLLD7bBpuQkxcnE64euopOHRIx+6feQYGDICwsFMOTUjQsfuwMO3dv/WWlua5994Atd2YHMbXHv8E4KiIRAKDgM1Auv0nERHgbWCdc+7lFA99DnTz7ncDPstUi03o+flnTbQfMECD/i23aInMRx45I+ivX69VNEeMSN7XuLGOBOXx9bfdmFzO1/8Kcd6wzC3Aa8658cAFGZzTBM3+aSEiMd7tBjQb6FoR2Qi0IuPsIBOq9u/XRVAaNNDVTi65BGbN0tull6Z6yu7dsHSpXuM9bmvEGZMqX4d6DonIEDSQNxWRPOg4f5qcc9+TduZPS9+baEJO4hXZgQPh77819WbgQC29UKjQGYfv2qV1dUCraU6dCm3b2sxbY9Lia4//TnSh9Xucc38D5YBRfmuVCV0bN+rKV506adBv3Bh++UXXwT0t6DsHDz2kk3PXr0/e37kzFC2ara02Jqj4Wp3zb/QibQFv117gU381yoSg48f1Ym2NGvDNN5qh89Zbmo5To0aqp4jA4cNw8qQm9hhjfONryYaewExgkrfrYmCWn9pkQs033+gkrGHD4MQJ6N5du/D33nvGFdk9e7SycqJRo2D5csvYMSYzfB3q6YNerD0I4JzbCFzor0aZEPH33zouc+21OsRTtSosXAiTJ6e66smSJXpIt26asglaTt8KqxmTOb4G/hPOuZOJGyKSjzQmXhmTofh4mDABqlSBDz7Qq7AjRkBMjE6xTUOVKnqdNyxMSykbY86Or1k9i0TkCeA8EbkW6A3M9l+zTK61YoXWyf/5Z92+/notm1zpzOWb4+Nh5kydqCsCJUrAjz9C+fJWVM2Yc+Frj/9xYA+wCrgfmAs85a9GmVzo4EFNwYmK0qB/8cUa1efMSTXoA7Rrp8XUJk9O3lehggV9Y85Vhj1+EckLrHHOVQHe9H+TTK7iHHz8sc66/esvvVj78MOawXNB+nMAO3bUTM4yZbKprcaEiAwDv3MuXkQ2iMilzrlt2dEok0v8/rvOvP3yS92uXx8mToTatVM9/PvvYd8+rcgAGvhvvjnVOVvGmHPg6xh/MWCNiPwEHEnc6Zy72S+tMsHtxAkYPVpLJh8/DkWKwMiRWjAnjQL4P/6oNXaKF9cyPBdeqEM6FvSNyXq+Bv6n/doKk3ssWqQXbxOn0nbuDC+9lFxTIQ3168ONN+qXgcKFs6GdxoSwdAO/iIQDDwCXoxd233bOxWVHw0yQ2b0bHn00uej9lVdqGeWWqZdl2rVLKyz/97+asi8Cn31mFTSNyQ4Z9fjfA2KB74DrgWrAAH83ygSRhAR4+214/HH4918oUEAXOX/8cb2fht694ZNPNGXznXd0nwV9Y7JHRoG/mnOuBoCIvA385P8mmaDx6686rLNsmW5fe6328i+/PMNTX3xRg/5//uPnNhpjzpBRHys28Y4N8Zgkhw/rsE6dOhr0y5TRMsrz5qUa9OPi9Fpv797J+y67TMvqV6iQba02xngy6vFHikji5HhBZ+4e9O4755xdhgs1n30G/frBn3/qwHzfvpq9U6RImqds26bj+SdO6BeEmjWzsb3GmDOkG/idc6nn3pnQs3Ur9O8Pn3+u23Xrak5+vXqpHh4XB/m8365KlWDsWChXzoK+MTmBXU4z6YuN1drH1app0L/gAhg3ThPv0wj6ixZB9erw7bfJ+3r1ghtuyKY2G2PSZYHfpG3JEh3Hf+wxOHoU7rxT8/P79k1zIhbo2im//QavvJJ9TTXG+M7XCVwmlOzbp+mYb7+t25ddptk6rVuneco//+isW9C/E0WKwP33Z0NbjTGZZj1+k8w5LYVZubIG/fz5Nd9y1ao0g/6//0L79tCwoVZnAD2tXz/91xiT81jgN2rNGl0E5Z57tMffooXm6T/zDJx3XpqnnX++jv7s3KmVNI0xOZ8N9YS6o0fhuec00T4uTqujvfwydOqUZuH7TZvgkkt0Ym6BAprCX7So7jPG5HzW4w9lc+Zo+s3IkTqN9v77tfveuXOaQf+ddyAiQldKTFSjhgV9Y4KJ9fhD0fbtujDKJ5/odmSk5uQ3bJjhqVdcASdP6tCOc7YaljHByHr8oSQuDsaMgapVNegXLKjDOsuXpxn0jx6Fr79O3m7aFNauhTfftKBvTLCyHn+o+OEHrZewcqVu33abJtqnM0Zz+LDWx9+6FWJidA4XQJUqfm+tMcaPrMef2/37rwb8xo016JcvD7Nn6zq4GQzMFyoE11yjXxDirESfMbmG9fhzK+dg2jQYNEgXScmXDx55BJ5+WnMw0zjlo490yL9yZd338suauRMWlo1tN8b4lQX+3GjDBnjwQViwQLebNoUJEzSDJx3jx+vEq6ZNYeFCXRjF1rw1JvexoZ7c5NgxnWlbs6YG/RIldCZuYtW0DHTqpCsmdu6cDW01xgSM9fhzi3nzoE8f2LxZt++9F154QYN/Gtav1+ycUaO0d1+8uGbspFN/zRiTC1jgD3Z//QUPPwwffqjbERGak9+kSbqnxcZq+Z0//9RTevTQ/Rb0jcn9LPAHq/h4rZj55JNw6JBesB02DB56yKcrsWFhuu7tN99Au3b+bqwxJiexwB+Mli/XFM3oaN2+6SZdHKV8+TRPOXxYlz+sXRu6ddN9d92lN2NMaLGLu8HkwAFNu6lfX4P+JZfAp5/qyljpBH2AuXPh1Vd1jfRjx7KpvcaYHMl6/MHAOZgxQ8fy//5bB+IffhiGDk033zIhQS/aAnTooJN377473SrLxpgQ4Lcev4i8IyK7RWR1in3FReRrEdno/VvMX6+fa2zaBNddBx07atBv1AhWrNBUnDSCvnPwwQeawbl7t+4T0clYdepkY9uNMTmSP4d63gWuO23fYGC+c+4KYL63bVJz4gQ8+6ym3Hz1FRQrBm+8Ad9/r3n6GZg8WdM133orG9pqjAkqfhvqcc4tFpEKp+2+BWju3X8PWAg87q82BK3586F3b12xHKBrV+3hX3hhmqfEx+vShwULau9+0iSdw3XPPdnUZmNM0MjuMf7Szrmd3v2/gdJpHSgivYBeAJdeemk2NC0H2LVLa+tMm6bbVapoqYXmzdM97bffNFPniitgyhTdV6mS3owx5nQBy+pxzjnApfP4G865es65eqVKlcrGlgVAQoJOuqpSRYN+eLgucbVyZYZBH/QCbkyMflHYt8/vrTXGBLnsDvy7RKQsgPfv7mx+/ZwnJkZLJj/4IOzfrxdy16yBJ56A/PnTPC2xMgPA5ZfDZ59puYV0KjQYYwyQ/YH/c8CbPkQ34LNsfv2c49AhGDgQ6taFH3+Eiy6CmTM14T6DMZpHHtFhnW++Sd7XujUUKeLnNhtjcgV/pnP+D1gGVBaR7SJyLzASuFZENgKtvO3Q4pwuglK1qi6DCFpmYd06aN/ep/UMixfXVP61a/3bVGNM7iQ61J6z1atXzy1fvjzQzTh3f/wBfftqrx50Bu7EiVpHIR27d+vi5pGRuh0bqxd0fai0bIwJYSIS7Zyrd/p+K9mQHU6ehOef10Vr587VMZnXX4elSzMM+itX6peD229PLrUQFmZB3xhz9qxkg78tWqQXbtet0+3OnWH0aChTxqfTq1bV4f+yZeHgQSu3YIw5d9bj95c9e6B7d03HXLcu+Wrs1KnpBv34eF0c5fhx3c6fH779VtdZKZ3mrAdjjPGdBf6slpCgdRIqV4b33tOVyp95Bn79FVq2zPD07t2hVy947rnkfaVK+XTN1xhjfGJDPVlp1Sqtk790qW63aqVj+Vdc4fNTPPAALF6sqf3GGOMPFvizwpEj2qt/+WUdqylTRlM177wzw6764sWalvnAA7rdpAls3Jju3C1jjDknFvjP1eef6+Io27ZpkO/TB4YPh6JFMzx182a45hrNyW/WTJN+wIK+Mca/LPCfrW3boH9/rZUAWuh+4kSIivL5KS67TNP6ixXT+8YYkx3s4m5mxcZqieSqVTXoX3ABvPKKll3IIOj//beup5KY2Qm6HOKwYXoN2BhjsoP1+DNjyRIdjF/tLSp2xx06ln/RRT6dPmIETJ8O//4LX37px3YaY0w6LPD7Yt8+GDw4eTmrSpVg/HitpJkB55Kv7z77rK6XnjJV0xhjspsN9aTHOc3Fr1JFg35YGDz1lPb4Mwj6cXHwwgtw/fWa2g86lj9lCpQvnw1tN8aYNFiPPy1r1+ryh4sW6Xbz5roaVpUqPp1+8KBmd+7erUsg+jB3yxhjsoUF/tMdParpmKNGabe9VCmN4J07Z5iTf+yYLp4loqWT334b8uWzoG+MyVlsqCeluXO17OXzz2vQv/9+2LAB7r47w6C/aBHUqAHvvJO878YbfboMYIwx2coCP8D27Vr3uG1b2LIFatbUsgsTJ+rAvA927NAJWe+9p5cGjDEmpwrtoZ64OHjtNXj6aTh8GAoW1NILAwboGE06nNO/F5dcotsdO+q+Dh2soJoxJmcL3R5/4oSrhx/WoH/bbTqzatCgDIP+kSNw6626ItauXbpPRC8DWLkFY0xOF3qBf/9+zdZp1AhiYjS3cvZsXQc3sfuegfPP1wu5cXFabdkYY4JJ6Az1OAcffAADB2qOZb582rt/+mkd4snAb79ppk7Jktq7f+st/bdcuWxouzHGZKHQ6PFv2KC18e++W4P+VVdpb3/kSJ+C/owZer33oYeS911yiQV9Y0xwyt2B//hxGDpUo/a330KJEppvuWhRplYrr1dPSyeHhenwjjHGBLPcO9SzbRu0aKE5lgD33KM1FEqWzPDUI0e0zH7Hjrp92WU61HPxxX5srzHGZJPcG/gvvlgH5cPDtdRC06Y+nRYXB/Xra8WGEiWgdevkpzPGmNwg9wb+vHnh00+15EImcizz5YMuXbR8cqlSfmyfMcYEiLggmGZar149t3z5cr88t3N68bZsWbj6at0XF6f7w8L88pLGGJMtRCTaOVfv9P25t8fvo+nToVMnLbG/ejWcd16G87eMMSaohXyIa99eh/+7dLHlD40xoSF3p3OmYu1aTec/dky38+fX7M6ePSFPyP00jDGhKKRCnXPas582DUaPTt5vRdWMMaEkJIZ6Ete9FdFinJMnQ79+gW6VMcYERq4O/IcO6RK5xYrBsGG6r1EjvRljTKjK1YF/1SoYO1Yzdfr29WnSrjHG5Hq5eoy/cWMdy1+61IK+McYkytU9ftDKy8YYY5Ll6h6/McaYM1ngN8aYEBOQwC8i14nIBhHZJCKDA9EGY4wJVdke+EUkLzAeuB6oBnQUkWrZ3Q5jjAlVgejx1wc2Oed+d86dBKYDtwSgHcYYE5ICEfgvBv5Msb3d23cKEeklIstFZPmePXuyrXHGGJPb5diLu865N5xz9Zxz9UrZiijGGJNlAhH4dwCXpNgu5+0zxhiTDbJ9BS4RyQf8BrREA/7PQCfn3Jp0ztkDbM2eFp6zksDeQDfCT+y9Ba/c/P7svaWtvHPujCGTbJ+565yLE5G+wDwgL/BOekHfOydoxnpEZHlqS53lBvbegldufn/23jIvICUbnHNzgbmBeG1jjAl1OfbirjHGGP+wwJ/13gh0A/zI3lvwys3vz95bJmX7xV1jjDGBZT1+Y4wJMRb4jTEmxFjgzwQRuUREFojIWhFZIyIDvP3FReRrEdno/VvM2y8iMtarQvqriNQJ7DvImIjkFZFfROQLb7uiiPzovYcZIpLf21/A297kPV4hoA33gYgUFZGZIrJeRNaJSKPc8tmJyMPe7+RqEfmfiIQH82cnIu+IyG4RWZ1iX6Y/KxHp5h2/UUS6BeK9nC6N9zbK+738VUQ+FZGiKR4b4r23DSLSJsX+s69y7Jyzm483oCxQx7t/AToRrRrwIjDY2z8YeMG7fwPwf4AADYEfA/0efHiPA4EPgC+87Q+Bu7z7E4EHvfu9gYne/buAGYFuuw/v7T3gPu9+fqBobvjs0FpXfwDnpfjMugfzZwc0A+oAq1Psy9RnBRQHfvf+LebdL5ZD31trIJ93/4UU760asBIoAFQENqPzn/J69yt5v8srgWo+tyHQP4RgvgGfAdcCG4Cy3r6ywAbv/iSgY4rjk47LiTe0fMZ8oAXwhfcfaW+KX8hGwDzv/jygkXc/n3ecBPo9pPPeinjBUU7bH/SfHcmFD4t7n8UXQJtg/+yACqcFx0x9VkBHYFKK/accl5Pe22mP3QpM8+4PAYakeGye91kmfZ6pHZfRzYZ6zpL39bg28CNQ2jm303vob6C0d9+nSqQ5yCvAY0CCt10C2O+ci/O2U7Y/6b15jx/wjs+pKgJ7gMneUNZbIlKQXPDZOed2AKOBbcBO9LOIJvd8doky+1kFzWd4mnvQbzDgp/dmgf8siEgh4GPgIefcwZSPOf3zG3Q5siJyI7DbORcd6Lb4ST706/UE51xt4Ag6XJAkiD+7YuiaFhWBi4CCwHUBbZSfBetnlREReRKIA6b583Us8GeSiIShQX+ac+4Tb/cuESnrPV4W2O3tD6ZKpE2Am0VkC7o4TgvgVaCoV1gPTm1/0nvzHi8C7MvOBmfSdmC7c+5Hb3sm+ocgN3x2rYA/nHN7nHOxwCfo55lbPrtEmf2sgukzRES6AzcCnb0/bOCn92aBPxNERIC3gXXOuZdTPPQ5kJgx0A0d+0/c39XLOmgIHEjxVTVHcc4Ncc6Vc85VQC/4feuc6wwsAG73Djv9vSW+59u943NsD8w59zfwp4hU9na1BNaSCz47dIinoYic7/2OJr63XPHZpZDZz2oe0FpEinnfilp7+3IcEbkOHWa92Tl3NMVDnwN3eZlYFYErgJ/QqsZXeJlb+dH/s5/7/IKBvsgRTDfgKvTr5a9AjHe7AR0fnQ9sBL4BinvHC7q+8GZgFVAv0O/Bx/fZnOSsnkreL9om4COggLc/3Nve5D1eKdDt9uF91QKWe5/fLDTTI1d8dsAzwHpgNfA+mgUStJ8d8D/0ekUs+m3t3rP5rNDx8k3erUeg31c6720TOmafGFcmpjj+Se+9bQCuT7H/BjSzcDPwZGbaYCUbjDEmxNhQjzHGhBgL/MYYE2Is8BtjTIixwG+MMSHGAr8xxoQYC/wmRxOREiIS493+FpEdKbbzB7p9KYlIcxFp7MfnP09EFolWUK1wWnXHniIS7eWsjxaRFv5qhwl+AVls3RhfOef2ofn3iMgw4LBzbnSg2iMi+Vxy/ZvTNQcOA0uz6PlOdw/wiXMuXudpJT1HF6Af0MI596+IjAPeBL71tR0mtFiP3wQdEanr9XyjRWReimn8C0VkjIgsF623HyUin3i12Id7x1Tw6p5P846ZKSLn+/C8r4jIcmCAiNwkWsf+FxH5RkRKe0X7HgAe9r6NNBWRd0Xk9hTtPuz921xEvhORz4G1Xg9+lIj8LFqP/f403npnkmerJj7nHWjNodbOub0AzrmtQAkRKZNVP3OTu1jgN8FGgHHA7c65usA7wIgUj590ztVD689/BvQBIoDuIpJYgbIy8LpzripwEOjt1WBK73nzO+fqOedeAr4HGjot9jYdeMw5t8V7zTHOuVrOue8yeB91gAHOuSvRmZsHnHNRQBTQ05uen/ymdVirkvc6icoDr6FB/+/Tnn8FWq/HmDPYUI8JNgXQQP61N9yRF53+niixXskqYI3z6uuIyO9oUav9wJ/OuSXecVOB/sCXGTzvjBT3ywEzvG8E+dE6/5n1k3Mu8bzWQM0U3w6KoDVZUj5vSa/tKe0B/gHuAMac9thutFKnMWewwG+CjaABvVEaj5/w/k1IcT9xO/H3/fQ6Jc6H5z2S4v444GXn3Oci0hwYlsY5cXjfqkUkD/pHIrXnE6Cfcy69AmLH0Bo7KR1F67V8JyK7nXMpS/mGe+cYcwYb6jHB5gRQSkQagZbJFpHqmXyOSxPPBzqhQzcbMvG8RUgugZtyHddD6JKcibYAdb37NwNhaTzfPOBBb7gJEblSdJGYJM65f4G8IhJ+2v7daO39/0qK9ViBK9GCbcacwQK/CTYJaCnhF0RkJVrJMLMplBuAPiKyDq3QOcE5dzITzzsM+EhEotFlCxPNBm5NvLiLZtZc7T1fI07t5af0FlpGeYWXojmJ1L+Nf4VWiD2FN2R0M/COiNT3/oBcjlYiNeYMVp3ThBQv++YL51xEoNuSWSJSB3jYOdclg+NuBeo4557OnpaZYGM9fmOChHNuBbBARPJmcGg+4KVsaJIJUtbjN8aYEGM9fmOMCTEW+I0xJsRY4DfGmBBjgd8YY0KMBX5jjAkx/w9dqFwKpTcwvgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# part B\n",
    "Rg = 8.314    # J/mol-K\n",
    "Vm = 0.00224  # m^3/mol\n",
    "Tc = 514      # K\n",
    "Pc = 6.137e6  # Pa\n",
    "a = 25*(Rg*Tc)**2/(64*Pc)\n",
    "b = Rg*Tc/(8*Pc)\n",
    "\n",
    "def P2(Vm,T):\n",
    "    # return pressure in bar\n",
    "    return 1e-5 * (Rg * T / (Vm-b) - a / Vm**2)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(T,P(Vm,T),'r-',linewidth=2,label='Ideal Gas Pressure')\n",
    "plt.plot(T,P2(Vm,T),'b:',linewidth=2,label='Non-Ideal Gas Pressure')\n",
    "plt.xlabel('Temperature (K)')\n",
    "plt.ylabel('Pressure (bar)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 8: Parameter Estimation and Statistics\n",
    "\n",
    "The yield (concentration) of green fluorescent protein produced in a reaction over time has been determined and is shown below.\n",
    "\n",
    "```\n",
    "Time (minutes)  Concentration (mg/mL)\n",
    "15              107.32\n",
    "30              203.05\n",
    "45              341.26\n",
    "60              401.24\n",
    "180             844.01\n",
    "480             1135.12\n",
    "720             1374.70\n",
    "1200            1651.26\n",
    "```\n",
    "\n",
    "Fit the data to the equation \n",
    "\n",
    "$G = A t^2 + B\\ln{t}+D$\n",
    "\n",
    "Where $G$ is the concentration of GFP in mg/mL at a given time ($t$) in minutes and $A$, $B$, and $D$ are adjustable parameters to minimize the difference between predicted and measured $G$.\n",
    "\n",
    "Solve for $A$, $B$, and $D$ to fit the data. Create a plot of the measured and predicted values and determine the $R^2$ value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = 0.00018907725436787096\n",
      "B = 311.39815072936636\n",
      "D = -815.2247831298712\n",
      "R^2 = 0.9932723267947231\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqW0lEQVR4nO3deZgU1dXH8e8BFcENVFQWYVhN0CjqCG5xF9dgFDQoRowLcYtroiLxRZOQKBrX+BoxElFHXAAFcWFTouZFFNCERdFBBSEgA0ZQh53z/nFroBln6Znpnurl93mefqb7Vk3XKRv7TNW991xzd0RERKrSIO4AREQk8ylZiIhItZQsRESkWkoWIiJSLSULERGp1jZxB5Auu+++uxcUFMQdhohI1pgxY8Zyd29e0bacTRYFBQVMnz497jBERLKGmS2obJtuQ4mISLWULEREpFpKFiIiUi0lCxERqZaShYiIVEvJQkREqqVkISIi1VKyEBHJAc889h1XNH+em+1OCgqgqCi175+zk/JERHJeaSm8+ioL7nqOn0wbRx9KWcje/HnB9fTvvy0Affum5lC6shARySalpTBqFPTpA82bQ+/e7PDeGwynH8fwBu34jA1sS2kpDByYusPqykJEJNNFVxA8/zyMGwfffRcSxQUXwNln0+L4o9hQwdf5woWpC0HJQkQkE333HbzySkgQL78cEkbz5nD++XDOOXDUUbBN+Apv1RYWVFDVqU2b1IWjZCEikim++y4khuefD4mitBT22GPzFURigkg0eDD07x92L9OkSWhPFSULEZE4ffvt1gli9eqQIPr1Cwnixz+uMEEkKuvEHjgw3Hpq0yYkilR1boOShYhI/Vu1KvQ9jBwZ+iLWrIG99oKLLoLevUOCaNiwRm/Zt29qk0N5ShYiIvXh66/hpZfCFcT48bBuHbRsCZdeGhLEEUfUOEHUJyULEZF0+eorGDs2JIiJE2H9emjdGq64Anr1gsMPhwbZMYNByUJEJJVKSmDMmHCLafJk2LAB2raFq68OVxDdumVNgkikZCEiUldLl8ILL4QEMWUKbNoEHTrADTeEK4jCQjCLO8o6UbIQEamNRYtg9OiQIN5+G9xhn31gwICQILp2zfoEkUjJQkQkWZ99FkptjBoF77wT2vbbDwYNCgli331zKkEkUrIQEanKxx+H5DByJMycGdoOOihMZOjVK1xN5AElCxGRRO4wZ86WBDF7dmjv3h2GDAkJon37eGOMgZKFiIg7zJgR+iBGjQpXE2Zw5JFw//1w5pmw995xRxkrJQsRyU+bNsHUqSE5jB4dKvE1bAjHHgvXXQc//WmYVS1AGtezMLNhZrbMzGYntN1mZovN7IPocWrCtgFmVmxm88zspIT2k6O2YjO7OV3xikge2LABXn8drrwyTI478kh46KHQST1sGHz5ZZg8d9llShTlpPPK4nHgL8AT5drvdfe7ExvMrAvQB9gXaAlMMrPO0eaHgBOBRcB7ZjbW3eemMW4RySXr1oXJcaNGhclyy5dD48Zw8smh/+H002GXXeKOMuOlLVm4+5tmVpDk7mcAz7j7WuAzMysGukXbit39UwAzeybaV8lCRCpXWhrqL40aFeoxrVoFO+0UEkOvXiFR7LBD3FFmlTj6LK4yswuA6cAN7v5foBXwTsI+i6I2gC/KtXev7I3NrD/QH6BNKlf9EJHMt2pVKPE9atSWtSB23TUkh1694IQToFGjuKPMWvWdLB4Gfg949PPPwEWpenN3HwoMBSgsLPRUva+IZKgVK8KVw6hRMGFCuOW0115hLYizzoKjj4Ztt407ypxQr8nC3b8se25mjwLjopeLgcRxaa2jNqpoF5F8tHQpvPhiSBBvvAEbN4bVfq68MiSIww7L6FLf2apek4WZtXD3JdHLM4GykVJjgafN7B5CB3cn4F3AgE5m1o6QJPoA59VnzCKSARYsCMNbR4+Gf/4zzIvo3BluvDEkiIMPztkyG5kibcnCzEYAxwC7m9kiYBBwjJl1JdyG+hz4JYC7zzGz5wgd1xuAK919Y/Q+VwHjgYbAMHefk66YRSSDlJXZGDUqTJgDOOAAuO22kCByuA5TJjL33Ly1X1hY6NOnT487DBFJljvMmrUlQcyJ/i7s3j0kh7POgo4d440xx5nZDHcvrGibZnCLSHzc4b33tsyiLi4OCwP9+Mcqs5FhlCxEpH5t3Bj6HcoSxKJFsM02cNxx8JvfhDIbe+wRd5RSjpKFiKTf+vVh5NKoUWEk07JlYc7DSSeFUt+nnx7mREjGUrIQkfRYsybMfRg9GsaOhf/+N8yaPu20MEnulFPCrGrJCkoWIpI6334Lr74ariBefjm8btoUevYMHdQ9eoS6TJJ1lCxEpG5Wrtwyi/q118IVRfPmcO654Qri2GNhu+3ijlLqKG0lykUkfkVFUFAQBhgVFITXKbF8OTz2GJx6akgMP/85vPsuXHIJTJkCS5bA0KGhT0KJIifoykIkRxUVQf/+oZ4ehEnQ/fuH53371uINly6FF14IVxBTpoRRTQUFcPXV4Qqie/eQlSQnaVKeSI4qKAgJory2beHzz5N8k0WLtiw1+tZbW8ps9O4dEsSBB2oWdQ7RpDyRPLRwYc3aN/v885AcRo6Ed6KVA/bbDwYNCglCZTbykpKFSI5q06biK4sKl3opLt6SIMquyA88MMyB6NUL9tknrbFK5lOyEMlRgwdv3WcB0KRJaAdg3ryQHEaOhA8+CG3dusGQISFBtG9f3yFLBlOyEMlRZZ3YAweGW09t2sCDV8zlJ/NHwo+eh9nRCgGHHw733BPmQbRtG1/AktGULERyWN/znL4HzIHnnw+Pmz4M/Q0//jE88EBIEK1aVf9GkveULERyjXu4aihLEB99FBLEUUdtWU2uRYu4o5Qso2Qhkgvcw/oPzz0XHvPmhTkPRx0Fv/pVSBB77RV3lJLFlCxEslVigii7gmjQAI4+Gq65JiSIPfeMO0rJEUknCzPbAVhTttypiMTkww/h2WdDkvjwwy0J4uqrlSAkbSpNFmbWAOgD9AUOAdYCjcxsOfAy8Ii7F9dLlCL57uOPQ3J49tnQH1HWB3HVVWGYqxKEpFlVVxZvAJOAAcBsd98EYGa7AscCd5rZC+7+VPrDFMlDn30WksOzz26ZB3HkkWEUU69e0LJlrOFJfqkqWZzg7uvLN7r7V8AoYJSZbZu2yETy0aJFof/hmWdCFVcIBfruuQfOPhtat443PslblSaLihJFGTNb6O5tqtlnGHA6sMzd94va7gJ+AqwD5gO/cPevzawA+BCYF/36O+5+WfQ7BwOPA42BV4BrPFerH0p+WrYszKJ+5plQrA/goIPgzjvhnHNCRUCRmNW2nnAyVcQeB04u1zYR2M/d9wc+JtziKjPf3btGj8sS2h8GLgU6RY/y7ymSfb7+Gv7+97DeQ8uWYf7DihXwu9+F/okZM+DGG5UoJGPUduhstX/Zu/ub0RVDYtuEhJfvAL2reg8zawHs7O7vRK+fAH4KvFrDeEXit3o1jBsHI0aEJUfXrYN27UJSOPfcUNlV1VwlQ1U1Gur6yjYBO6bg2BcBzya8bmdm7wOrgN+6+1tAK2BRwj6LoraKAzPrD/QHaFNhaU2RerZhA0yeDE8/HRYO+uabMDnu8stDgujWTQlCskJVVxY7VbHt/roc1MwGAhuAskUelwBt3H1F1EfxopntW9P3dfehwFAIix/VJUaRWnOHadPCUnXPPgslJdC0aeh/OO+8MCeiYcO4oxSpkao6uG9PxwHN7EJCx/fxZR3V7r6WMI8Dd59hZvOBzsBiIHH4R+uoTSTzzJsXEkRREXz6KWy/PZx+eij/esop0KhR3BGK1Fq1fRZm1g74FVCQuL+796zpwczsZOBG4Gh3L01obw585e4bzaw9oSP7U3f/ysxWmdmhwDTgAuDBmh5XJG2+/DKMYnrqqbBoUIMGcNxxcOutYTb1zjvHHaFISiTTwf0i8BjwErAp2Tc2sxHAMcDuZrYIGEQY/dQImGjhPm3ZENmjgN+Z2froGJdF8zkArmDL0NlXUee2xK20FF58EZ58EiZOhI0bw1DXe+6BPn1U0VVyklU3ZcHMprl793qKJ2UKCwt9etnykCJ1tWkTvPkmPPFEmDT37bdhNaHzzw+3mbp0iTtCkTozsxnuXljRtmSuLO43s0HABKJ+BQB3n5mi+EQy18cfhwTx5JNhubmddgod1RdcEBYQalDbqUoi2SWZZPEj4OfAcWy5DeXRa5Hc8/XXYRTT8OEwdWpICD16hBnVPXuGhaxF8kwyyeJsoL27r0t3MCKx2bgxzIf4+9/DfIi1a2HffWHIkHCbSUX7JM8lkyxmA02BZekNRSQGxcXw+OPhKmLRImjWDC69FPr1g4MP1oQ5kUgyyaIp8JGZvcfWfRY1HjorkhFKS0PhvmHD4B//CLeZTjopjGbq2VPzIUQqkEyyGJT2KETSrOgpZ8RvZnL60r9xnj3Nzr4KOnaEwYPDVUSrSqvIiAhV14YaD7wGvOruH9VfSCIptHIl715TxH5PPMo4/4DVbM/zfjZFjS7mgkFH0fd83WYSSUZV4/76Af8FbjOzmWb2sJmdEa3FLZK53OGdd+Cii6BFC7oNv5JNDlfwEC1YQj+eYMLaoxn4WyUKkWRVVRtqKWHm9OPRetzdgVOAG81sNTDB3YfUS5QiyVi1KpTd+OtfYdYs2GEHOP98uj16Ke9RSPllWBYujCdMkWyU1HoW0frbU6PH/5jZ7sBJ6QxMJGnvvw8PPxzKgH/3HRx4YEgY550HO+3EsgnAgu//mqrYiyQvmUKCD/L9xY5WAqqlIfFZsyaU3fjf/w23nBo3DutDXHYZFBZuNeR18GDo3z8MgirTpEloF5HkJFOroBHQFfgkeuxPKBV+sZndl7bIRCqycCEMGAB77x1Kbnz1Fdx3HyxeDI89Bocc8r25EX37wtCh0LZt2NS2bXjdt288pyCSjZK5DbU/cIS7bwQws4eBt4AjgVlpjE0kcIcpU+DBB2HMmNDWs2dYt/q445Kqz9S3r5KDSF0kkyyaEZZRXRm93gHYNVp7Ym3lvyZSR6tXh4WEHnggdFjvtltYr/ryy9XhIFLPkkkWQ4APzGwKYTjJUcAfoyG0k9IYm+SrJUvgoYdCJ/WKFbD//uEW07nnhr4JEal31SYLd3/MzF4BukVNt7j7f6Lnv0lbZJJ//vWvUHJjxAjYsCHcarr22rBmtWo0icQqqaGzQHPCiKhtgEPNDHcfnb6wJG+4w4QJcPfdMGlSmBtx2WVwzTXQoUPc0YlIJJmhs8MIndxz2Ho9CyULqb3168OaEUOGhP6Ili3hjjvCGNdmzeKOTkTKSebK4lB315qRkhrffQd/+1u43bRwYVgz4vHHQ3/EdtvFHZ2IVCKZZDHVzLq4+9y0RyO56+uvQ6f1fffB8uVw5JHh9amnamlSkSyQTLJ4gpAwlhLWszDA3X3/tEYmuWH5crj3XvjLX0LtptNOC5Pqjjgi7shEpAaS+ZPuMcIa3CcDPwFOj35Wy8yGmdkyM5ud0LarmU00s0+in82idjOzB8ys2Mz+bWYHJfxOv2j/T8ysX01OUGLy5Zfwm9+E6dJ/+lNYXOj992HcOCUKkSyUTLIocfex7v6Zuy8oeyT5/o8Tkkyim4HJ7t4JmBy9hlDRtlP06A88DCG5EBZg6k4YvjuoLMFIBvryS/j1r6Fdu9AvcdZZMGcOPPccdO0ad3QiUkvJ3IZ638yeBl5i62VVqx0N5e5vmllBueYzgGOi58OBKcBNUfsT7u7AO2bW1MxaRPtOdPevAMxsIiEBjUgidqkvK1aEkU0PPghr18L558NvfwudOsUdmYikQDLJojEhSfRIaKvL0Nk93X1J9HwpsGf0vBXwRcJ+i6K2ytolE6xcGa4g7r0Xvv02FGC69Vbo3DnuyEQkhZKZwf2LdB3c3d3Mypc/rzUz60+4hUUb1Q5Kr9WrQ3nwP/4xVH7t3Rtuvx26aJS1SC6qtM/CzH4b9RdUtv04Mzu9Fsf8Mrq9RPRzWdS+GNg7Yb/WUVtl7d/j7kPdvdDdC5s3b16L0KRaGzeGeRGdO4e+iUMOgRkzwtoSShQiOauqK4tZwEtmtgaYCZQA2xM6oLsSigj+sRbHHEtY3/uO6OeYhParzOwZQmf2SndfYmbjCYULyzq1ewADanFcqQt3GD8+VH2dNQu6dYMnn4Rjjok7MhGpB1WtwT0GGGNmnYAjgBbAKuApoL+7r67uzc1sBKGDenczW0QY1XQH8JyZXUxY7PKcaPdXgFOBYqAU+EUUx1dm9nvgvWi/35V1dks9mT0brr8eJk4M9Zqeey7cdlJxP5G8YWHwUe4pLCz06dO18mudlJSEzupHH4VddoFBg8JaEirLIZKTzGyGuxdWtC3ZqrOST9avD53XgwaFWk5XXRWe71ppF5aI5DglC9nalCkhOcyZAyeeCPffDz/8YdxRiUjMVMFNgqVLw0S6Y48NVxMvvhg6tJUoRITk1rNoDlwKFCTu7+4XpS8sqTebNsEjj4TifqtXhz6KAQO0fKmIbCWZ21BjgLcIQ2U3pjccqVezZ4fFhqZOheOPD/0UmnktIhVIJlk0cfeb0h6J1J916+APfwjVYHfZBZ54ItyC0lBYEalEMn0W48zs1LRHIvVjxgwoLITf/x5+9jP48EP4+c+VKESkSskki2sICWONmX0TPValOzBJsXXrQn9E9+6hQuxLL8FTT4HKoohIEpIpJLhTfQQiaTR3brh6mDkT+vULFWKbaUkQEUleUvMszKwncFT0coq7j0tfSJIy7mE50xtvhB13hNGj4cwz445KRLJQtbehzOwOwq2oudHjGjP7U7oDk8oVFUFBATRoEH4WFVWwU0kJ9OwJV18dRjrNmqVEISK1lsyVxalAV3ffBGBmw4H3UeXXWBQVhdGupaXh9YIF4TWEdYcAeP31MLppxQp44IEwI1sd2CJSB8nO4G6a8HyXNMQhSRo4cEuiKFNaGtrZtCkMiT3hhDAk9t134Ve/UqIQkTpL5sriT4R1uN8AjNB3cXNao5JKLVxYcfu3C1bA6T+HV1+F884Ls7J33LF+gxORnJXMaKgRZjYFOCRqusndl6Y1KqlUmzbh1lOiA/iAcQ1/CpOXhFnYl12mqwkRSamqllX9QfTzIMLCR4uiR8uoTWIweDA0abLl9dk8xz85gl133gBvvRXWm1CiEJEUq+rK4nqgP/DnCrY5cFxaIpIqlXVi//aWTVyy8H8YyGBKOh1O8zdHwV57xRuciOSsqpZVjcbYcIq7r0ncZmbbpzUqqVLfs1bTd0w/WPg8XHIJzf/yF2jUKO6wRCSHJTMa6v+SbJP6sGwZHHccjBwJd98NQ4cqUYhI2lV6ZWFmewGtgMZmdiBhJBTAzkCTyn5P0mj+fOjRA5YsCcnirLPijkhE8kRVfRYnARcCrYF7Etq/AW5JY0xSkQ8+gJNPhg0b4I03QkFAEZF6UlWfxXBguJn1cvdR9RiTlPfmm/CTn8DOO4dEoaVORaSeJTPPYpSZnQbsC2yf0P672hzQzPYBnk1oag/8D2GW+KVASdR+i7u/Ev3OAOBiwkp9V7v7+NocOytNmABnnBGKQE2YAHvvHXdEIpKHklmD+6+EPopjgb8BvYF3a3tAd58HdI3euyGwGHgB+AVwr7vfXe74XYA+hGTVEphkZp3dPfeXeH3tNfjpT+EHP4CJE7X2hIjEJpnRUIe7+wXAf939duAwIFULNR8PzHf3BVXscwbwjLuvdffPgGKgW4qOn7lefjlcUXTpApMnK1GISKySSRZlcyxKzawlsJ4wozsV+gAjEl5fZWb/NrNhZla2Ok8r4IuEfRZFbd9jZv3NbLqZTS8pKalol+wwfnwoJ/6jH8GkSbDbbnFHJCJ5Lplk8ZKZNQXuAmYCnwNP1/XAZrYd0BN4Pmp6GOhAuEW1hIpnjlfJ3Ye6e6G7FzbP1r/E3347JIouXcKtp113jTsiEZGq+yzMrAEw2d2/BkaZ2Thge3dfmYJjnwLMdPcvAcp+Rsd9FChbjW8xkNir2zpqyz3vvw+nnRY6sSdM0NKnIpIxqryyiBY8eijh9doUJQqAc0m4BWVmibe2zgRmR8/HAn3MrJGZtQM6UYcO9oxVXAwnnRTWoZg4EfbYI+6IREQ2S2Y9i8lm1gsY7e6eioOa2Q7AicAvE5qHmFlXQpHCz8u2ufscM3uOsKTrBuDKnBsJtWIFnHpqWLxo0qRQh1xEJINYdd//ZvYNsAPhi3oNoeyHu/vO6Q+v9goLC3369Olxh1G9tWvhxBNh2rSwHOoRR8QdkYjkKTOb4e6FFW1LZlLeTqkPSQBwh4svDutQjBihRCEiGava0VBmNjmZNqmFO++EoqKwolGfPnFHIyJSqaqqzm5PmLm9ezTnIbHqbIXzHKQGJk2CgQPhZz+DAQPijkZEpEpV3Yb6JXAtocTGDLYki1XAX9IbVo5bsCBcSfzwh/C3v2kZVBHJeFVVnb0fuN/MfuXuD9ZjTLltzRro1QvWr4fRo2HHHeOOSESkWsl0cD9oZocDBYn7u/sTaYwrd910E8yYAS++CJ1TVWJLRCS9kqk6+yShDMcHhBLhEOZCKFnU1Pjx8MADcPXVoUigiEiWSGZSXiHQJVUT8vJWSQlceCHsuy/ccUfc0YiI1EgyyWI2sBehuJ/Uhjv07w9ffRXWqGjcOO6IRERqJJlksTsw18zeBdaWNbp7z7RFlWuGDw99FHffDQccEHc0IiI1lkyyuC3dQeS0khK44QY48ki47rq4oxERqZVkRkP9w8zaAp3cfZKZNQEapj+0HHHDDfDNN/DII9AgmeVDREQyTzLlPi4FRgKPRE2tgBfTGFPumDwZnnwyDJft0iXuaEREai2ZP3WvBI4gzNzG3T8BtNhCddasgcsug44d4ZZb4o5GRKROkumzWOvu6ywqSWFm2xDmWUhVhgwJCxpNmqTRTyKS9ZK5sviHmd0CNDazEwlrZr+U3rCy3NKlIVn07g3HHx93NCIidZZMsrgZKAFmEYoLvgL8Np1BZb3bbw+LGv3pT3FHIiKSEsnchmoMDHP3RwHMrGHUVprOwLLWRx/Bo4/C5ZeH/goRkRyQzJXFZEJyKNMYmJSecHLAgAHQpAncemvckYiIpEwyyWJ7d/+27EX0vEn6Qspib78dZmrfdBPsoQFjIpI7kkkW35nZQWUvzOxgYHX6QspigwbBXnvBtdfGHYmISEol02dxLfC8mf2HsFreXsDP6npgM/sc+IZQ9nyDuxea2a7As4S1Mz4HznH3/1oYt3s/cCqhr+RCd59Z1xhS6r334PXX4a67YIcd4o5GRCSlkin38Z6Z/QDYJ2qa5+7rU3T8Y919ecLrm4HJ7n6Hmd0cvb4JOAXoFD26Aw9HPzPHnXfCLruE6rIiIjkm2WJFhwD7AwcB55rZBWmK5wxgePR8OPDThPYnPHgHaGpmLdIUQ819/HFYIvXKK2HnneOORkQk5eJcKc+BCWbmwCPuPhTY093L1s1YCuwZPW8FfJHwu4uitq3W2DCz/kB/gDZt2tQxvBq46y7YbruwAp6ISA6Kc6W8I919sZntAUw0s48SN7q7R4kkaVHCGQpQWFhYPyVJ/vMfeOIJuPhi2HPP6vcXEclCydyGKlspL6XcfXH0cxnwAtAN+LLs9lL0c1m0+2Jg74Rfbx21xe/++2HDBvj1r+OOREQkbZJJFmUr5Y03s7Flj7oc1Mx2MLOdyp4DPQhJaSzQL9qtHzAmej4WuMCCQ4GVCber4rNuHQwbBmecAe3bxx2NiEjaxLVS3p7AC1El222Ap939NTN7D3jOzC4GFgDnRPu/Qhg2W0wYOvuLNMRUc+PGwfLlcMklcUciIpJWya6UtydhRBTAu9Gto1pz90+B7y1G7e4rgO+VaY36S66syzHT4rHHoGVL6NEj7khERNIqmZXyzgHeBc4m/KU/zcx6pzuwjLd4Mbz2Glx4IWyTzAWaiEj2SuZbbiBwSNnVhJk1JxQSHJnOwDLe8OGwaRP8IjPuiImIpFMyHdwNyt12WpHk7+Uu99CxffTRKkMuInkhmS/916KRUBea2YXAy8Cr6Q0rw735Jsyfz3WzL6ZBAygogKKiuIMSEUmfZDq4f2NmZwFHRk1D3f2F9IaV2ebf+nd2Z2ceWdELBxYs2FISqm/fWEMTEUmLSq8szKyjmR0B4O6j3f16d78eKDGzDvUWYaZZv57d/zmG0ZzF6oRlPUpLYeDAGOMSEUmjqm5D3QesqqB9ZbQtP739Nrts+pqx9PzepoULY4hHRKQeVJUs9nT3WeUbo7aCtEWU6caOZS2NmMiJ39tUn7ULRUTqU1XJomkV2xpXsS13ucOYMSw/4Hi8yY5bbWrSBAYPjikuEZE0qypZTDezS8s3mtklwIz0hZTB5s6Fzz6j1eU9GToU2rYFs/Bz6FB1botI7qpqNNS1hPpNfdmSHAqB7YAz0xxXxikqgsVXjeVGoPvvT+fqO+Hzz+OOSkSkflSaLNz9S+BwMzsW2C9qftndX6+XyDJIUVEYGjupdCzvUci7i1tpqKyI5JVqJ+W5+xvu/mD0yLtEAWFI7E6lS+nOtM2joDRUVkTySX6X7UjSwoVwGi/TAN9qyKyGyopIvlCySEKbNtCTsSygDf9m/63aRUTygZJFEgb/fhNH8w8m0AMwQENlRSS/KFkkoW/hPJqyko93O1xDZUUkL2nVnmRMnQrAXW8fxl0/iDkWEZEY6MoiGVOnQrNm0Llz3JGIiMRCySIZU6fCoYdCA/3nEpH8pG+/6qxcGcp8HHZY3JGIiMSm3pOFme1tZm+Y2Vwzm2Nm10Ttt5nZYjP7IHqcmvA7A8ys2MzmmdlJ9RrwtGmhgOChh9brYUVEMkkcHdwbgBvcfaaZ7QTMMLOJ0bZ73f3uxJ3NrAvQB9gXaAlMMrPO7r6xXqKdOjVUC+zevV4OJyKSier9ysLdl7j7zOj5N8CHQKsqfuUM4Bl3X+vunwHFQLf0RxqZOhX23Rd23rneDikikmli7bMwswLgQGBa1HSVmf3bzIaZWbOorRXwRcKvLaKS5GJm/c1suplNLykpqXuAmzbBO++ov0JE8l5sycLMdgRGAde6+yrgYaAD0BVYAvy5pu/p7kPdvdDdC5s3b173IOfNCx3cShYikudiSRZmti0hURS5+2gIJdHdfaO7bwIeZcutpsXA3gm/3jpqS79oMp6ShYjkuzhGQxnwGPChu9+T0N4iYbczgdnR87FAHzNrZGbtgE7Au/URa/GTU/m6QTMa/rAzBQVhXQsRkXwUx2ioI4CfA7PM7IOo7RbgXDPrCjjwOfBLAHefY2bPAXMJI6murI+RUEVF0PUf7/B/fiibaMCCBWjBIxHJW+bucceQFoWFhT59+vRa/377thv5aGFj7uF6BnDH5va2bbWcqojkJjOb4e6FFW3TDO5K+MIv2I71FNNxq3YteCQi+UjJohKH7TEfgPl02KpdCx6JSD5SsqjENacVA2x1ZaEFj0QkXylZVKL7bsVs3LYR27RppQWPRCTvafGjysyfT8OO7flsrvKpiIi+CStTXAwdO1a/n4hIHlCyqIg7zJ8PHTpUv6+ISB5QsqjI0qVQWqorCxGRiJJFgqIiKCiAo1qGkVCvL1SyEBEBJYvNiopCOY8FC6ADIVlc80AH1YMSEUHJYrOBA8OdJ4AOzGcDDfloTVsGDow3LhGRTKBkEUks49GRYj6ngA1sq/IeIiIoWWyWWMajI8Wby3yovIeIiJLFZoMHh3Ie4HSkmGI6qryHiEhEM7gjZWU87rr5K5ouWslXzToy9EGV9xARASWLrfTtC307FsOhcOvjHaBn3BGJiGQG3YYqb34oTa4JeSIiWyhZlFcc5ljQrl28cYiIZBAli/KKi6F1a2jcOO5IREQyhpJFefPn6xaUiEg5ShblFRer2qyISDlZkyzM7GQzm2dmxWZ2c1oOsnEjnHwyHHNMWt5eRCRbZcXQWTNrCDwEnAgsAt4zs7HuPjelB2rYEIYPT+lbiojkgmy5sugGFLv7p+6+DngGOCPmmERE8ka2JItWwBcJrxdFbVsxs/5mNt3MppeUlNRbcCIiuS5bkkVS3H2ouxe6e2Hz5s3jDkdEJGdkS7JYDOyd8Lp11CYiIvUgW5LFe0AnM2tnZtsBfYCxMcckIpI3smI0lLtvMLOrgPFAQ2CYu8+JOSwRkbyRFckCwN1fAV6JOw4RkXyULbehREQkRubucceQFmZWAixIcvfdgeVpDKc+5dK5QG6dj84lc+XS+dTlXNq6e4VDSXM2WdSEmU1398K440iFXDoXyK3z0blkrlw6n3Sdi25DiYhItZQsRESkWkoWwdC4A0ihXDoXyK3z0blkrlw6n7Sci/osRESkWrqyEBGRailZiIhItfI6WdTL6nspZmZ7m9kbZjbXzOaY2TVR+65mNtHMPol+NovazcweiM7x32Z2ULxn8H1m1tDM3jezcdHrdmY2LYr52ageGGbWKHpdHG0viDXwcsysqZmNNLOPzOxDMzssyz+X66J/Y7PNbISZbZ8tn42ZDTOzZWY2O6Gtxp+FmfWL9v/EzPrFcS5RHBWdz13Rv7V/m9kLZtY0YduA6HzmmdlJCe21/85z97x8EGpMzQfaA9sB/wK6xB1XEnG3AA6Knu8EfAx0AYYAN0ftNwN3Rs9PBV4FDDgUmBb3OVRwTtcDTwPjotfPAX2i538FLo+eXwH8NXreB3g27tjLncdw4JLo+XZA02z9XAjrxXwGNE74TC7Mls8GOAo4CJid0FajzwLYFfg0+tkset4sg86nB7BN9PzOhPPpEn2fNQLaRd9zDev6nRf7P8oY/zEdBoxPeD0AGBB3XLU4jzGE5WbnAS2ithbAvOj5I8C5Cftv3i8THoRy85OB44Bx0f+wyxP+J9j8OREKSR4WPd8m2s/iPoconl2iL1cr156tn0vZgmO7Rv+txwEnZdNnAxSU+3Kt0WcBnAs8ktC+1X5xn0+5bWcCRdHzrb7Lyj6bun7n5fNtqKRW38tk0aX+gcA0YE93XxJtWgrsGT3P9PO8D7gR2BS93g342t03RK8T4918LtH2ldH+maAdUAL8Pbql9jcz24Es/VzcfTFwN7AQWEL4bz2D7PxsytT0s8joz6iciwhXR5Cm88nnZJHVzGxHYBRwrbuvStzm4c+GjB8TbWanA8vcfUbcsaTANoTbBA+7+4HAd4RbHZtly+cCEN3PP4OQBFsCOwAnxxpUCmXTZ1EdMxsIbACK0nmcfE4WWbv6npltS0gURe4+Omr+0sxaRNtbAMui9kw+zyOAnmb2OfAM4VbU/UBTMysrn58Y7+ZzibbvAqyoz4CrsAhY5O7TotcjCckjGz8XgBOAz9y9xN3XA6MJn1c2fjZlavpZZPpnhJldCJwO9I0SIKTpfPI5WWTl6ntmZsBjwIfufk/CprFA2WiNfoS+jLL2C6IRH4cCKxMuxWPl7gPcvbW7FxD++7/u7n2BN4De0W7lz6XsHHtH+2fEX4fuvhT4wsz2iZqOB+aShZ9LZCFwqJk1if7NlZ1P1n02CWr6WYwHephZs+hKq0fUlhHM7GTCLdye7l6asGks0CcaodYO6AS8S12/8+LsgIr7QRgF8TFhhMDAuONJMuYjCZfP/wY+iB6nEu4PTwY+ASYBu0b7G/BQdI6zgMK4z6GS8zqGLaOh2kf/uIuB54FGUfv20eviaHv7uOMudw5dgenRZ/MiYQRN1n4uwO3AR8Bs4EnC6Jqs+GyAEYS+lvWEq76La/NZEPoCiqPHLzLsfIoJfRBl3wN/Tdh/YHQ+84BTEtpr/Z2nch8iIlKtfL4NJSIiSVKyEBGRailZiIhItZQsRESkWkoWIiJSLSULEcDMdjOzD6LHUjNbHD3/1sz+N03HvNbMLqjh7/xfEvs8Y2adah+ZyPdp6KxIOWZ2G/Ctu9+dxmNsA8wkVBDeUN3+NXzvo4Hz3f3SVL6v5DddWYhUwcyOsS3rbNxmZsPN7C0zW2BmZ5nZEDObZWavRWVYMLODzewfZjbDzMaXlZgo5zhgZlmiMLMpZnavmU23sBbGIWY2OlpH4Q8J8XybENcU27J+RlE00xrgLeCEhLIcInWmZCFSMx0IX/Q9gaeAN9z9R8Bq4LQoYTwI9Hb3g4FhwOAK3ucIQhXXROvcvZCwTsQY4EpgP+BCM6uoguuBwLWE9QvaR++Ju28izO49oPanKbI1/eUhUjOvuvt6M5tFWEzmtah9FmG9gX0IX/AToz/0GxLKNJTXAviwXFtZnZ5ZwByPakWZ2aeEAnDlC/O96+6Lon0+iI7/drRtGaFabC5U9JUMoGQhUjNrIfz1bmbrfUun3ybC/09G+KI/rJr3WU2op/S9947ea21Ce9l7VxhLZGO5fbaPjiGSEroNJZJa84DmZnYYhHLyZrZvBft9CHRMYxydCQUARVJCyUIkhdx9HaFE951m9i9CNdDDK9j1VcK6yilnZnsCqz2UTRdJCQ2dFYmJmb0A3Ojun6T4fa8DVrn7Y6l8X8lvurIQic/NhI7uVPsaGJ6G95U8pisLERGplq4sRESkWkoWIiJSLSULERGplpKFiIhUS8lCRESq9f8heSM9ucw/ygAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "x = np.array([15,30,45,60,180,480,720,1200])\n",
    "y = np.array([107.32,203.5,341.26,401.24,844.01,1135.12,1374.70,1651.26])\n",
    "\n",
    "def f(t,A,B,D):\n",
    "    return A*t**2 + B*np.log(t) + D\n",
    "\n",
    "popt, pcov = curve_fit(f,x,y)\n",
    "A = popt[0]\n",
    "B = popt[1]\n",
    "D = popt[2]\n",
    "print('A = '+str(A))\n",
    "print('B = '+str(B))\n",
    "print('D = '+str(D))\n",
    "\n",
    "plt.plot(x,y,'bo')\n",
    "xp = np.linspace(15,1200,100)\n",
    "plt.plot(xp,f(xp,A,B,D),'r-')\n",
    "plt.xlabel('Time (min)')\n",
    "plt.ylabel('Concentration (mg/mL)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R^2 = 0.9932723267947231\n"
     ]
    }
   ],
   "source": [
    "# calculate R^2\n",
    "from sklearn.metrics import r2_score\n",
    "rsq2 = r2_score(y, f(x,A,B,D))\n",
    "print('R^2 = '+str(rsq2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. Implicit Equation Solve\n",
    "\n",
    "Determine the two roots of the polynomial using `fsolve` from `scipy.optimize`.\n",
    "\n",
    "$x^2+3x+2=0$\n",
    "\n",
    "Verify the solutions from the quadratic formula:\n",
    "\n",
    "$\\frac{-b\\pm \\sqrt{b^2-4ac}}{2a}$\n",
    "\n",
    "with solutions:\n",
    "\n",
    "$x_0=\\frac{-3+\\sqrt{9-8}}{2} = \\frac{-2}{2} = -1$\n",
    "\n",
    "$x_1=\\frac{-3-\\sqrt{9-8}}{2} = \\frac{-4}{2} = -2$\n",
    "\n",
    "A better way to find polynomial roots is with `numpy` such as\n",
    "\n",
    "```python\n",
    "import numpy as np\n",
    "np.polynomial.polynomial.polyroots([2,3,1])\n",
    "```\n",
    "\n",
    "but use this exercise to practice `fsolve` and how to use a different initial guess to find alternate solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "def fcn(x):\n",
    "    return x**2+3*x+2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.]\n"
     ]
    }
   ],
   "source": [
    "z = fsolve(fcn,0)  # guess 0\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-2.]\n"
     ]
    }
   ],
   "source": [
    "z = fsolve(fcn,-3) # guess -3\n",
    "print(z)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}