01_Getting-started-project.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computational Mechanics Project #01 - Heat Transfer in Forensic Science\n",
    "\n",
    "We can use our current skillset for a macabre application. We can predict the time of death based upon the current temperature and change in temperature of a corpse. \n",
    "\n",
    "Forensic scientists use Newton's law of cooling to determine the time elapsed since the loss of life, \n",
    "\n",
    "$\\frac{dT}{dt} = -K(T-T_a)$,\n",
    "\n",
    "where $T$ is the current temperature, $T_a$ is the ambient temperature, $t$ is the elapsed time in hours, and $K$ is an empirical constant. \n",
    "\n",
    "Suppose the temperature of the corpse is 85$^o$F at 11:00 am. Then, 2 hours later the temperature is 74$^{o}$F. \n",
    "\n",
    "Assume ambient temperature is a constant 65$^{o}$F.\n",
    "\n",
    "1. Use Python to calculate $K$ using a finite difference approximation, $\\frac{dT}{dt} \\approx \\frac{T(t+\\Delta t)-T(t)}{\\Delta t}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Empirical constant k = 0.61\n"
     ]
    }
   ],
   "source": [
    "Ti = 85\n",
    "T = 74\n",
    "deltat = 2\n",
    "Ta = 65\n",
    "\n",
    "dTdt = ((T - Ti)/deltat)\n",
    "k = -dTdt/(T - Ta)\n",
    "print('Empirical constant k =', round(k,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Change your work from problem 1 to create a function that accepts the temperature at two times, ambient temperature, and the time elapsed to return $K$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "metadata": {},
   "outputs": [],
   "source": [
    "def measure_K(Temp_t1,Temp_t2,Temp_ambient,delta_t):\n",
    "    ''' Determine the value of K based upon temperature of corpse \n",
    "    when discovered, Temp_t1\n",
    "    after time, delta_t, Temp_t2\n",
    "    with ambient temperature, Temp_ambient\n",
    "    Arguments\n",
    "    ---------\n",
    "    Temp_t1,Temp_t2,Temp_ambient,delta_t\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    measure_K \n",
    "    \n",
    "    '''\n",
    "     \n",
    "    measure_K = -(((Temp_t2 - Temp_t1)/delta_t)/(Temp_t2 - Temp_ambient))\n",
    "    return(measure_K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 309,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Empirical constant K =  0.61\n"
     ]
    }
   ],
   "source": [
    "Temp_t1 = 85\n",
    "Temp_t2 = 74\n",
    "Temp_ambient = 65\n",
    "delta_t = 2\n",
    "\n",
    "print('Empirical constant K = ', round(measure_K(Temp_t1,Temp_t2,Temp_ambient,delta_t), 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. A first-order thermal system has the following analytical solution, \n",
    "\n",
    "    $T(t) =T_a+(T(0)-T_a)e^{-Kt}$\n",
    "\n",
    "    where $T(0)$ is the temperature of the corpse at t=0 hours i.e. at the time of discovery and $T_a$ is a constant ambient temperature. \n",
    "\n",
    "    a. Show that an Euler integration converges to the analytical solution as the time step is decreased. Use the constant $K$ derived above and the initial temperature, T(0) = 85$^o$F. \n",
    "\n",
    "    b. What is the final temperature as t$\\rightarrow\\infty$?\n",
    "    \n",
    "    c. At what time was the corpse 98.6$^{o}$F? i.e. what was the time of death?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[85.         72.52816327 67.83366211 66.06661355 65.40148205 65.15112112\n",
      " 65.05688322 65.02141131 65.00805939 65.00303362 65.00114188 65.00042981\n",
      " 65.00016179 65.0000609  65.00002292 65.00000863 65.00000325 65.00000122\n",
      " 65.00000046 65.00000017 65.00000007 65.00000002 65.00000001 65.\n",
      " 65.         65.         65.         65.         65.         65.\n",
      " 65.         65.         65.         65.         65.         65.\n",
      " 65.         65.         65.         65.         65.         65.\n",
      " 65.         65.         65.         65.         65.         65.\n",
      " 65.         65.        ]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "t = np.linspace(0,50,50)\n",
    "T_numerical = np.zeros(len(t))\n",
    "T_a = 65\n",
    "T_o = 85\n",
    "K = 0.61112\n",
    "dt = t[1] - t[0]\n",
    "T_numerical[0] = T_o\n",
    "\n",
    "for i in range(1,len(t)):\n",
    "    T_numerical[i] = T_numerical[i-1] + (-K*(T_numerical[i-1] - T_a)*dt)\n",
    "print(T_numerical)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[85.         75.7203139  70.7462565  68.08008367 66.65097319 65.88494754\n",
      " 65.47434577 65.25425678 65.13628562 65.07305123 65.03915661 65.02098856\n",
      " 65.0112502  65.00603028 65.00323233 65.00173258 65.00092869 65.00049779\n",
      " 65.00026682 65.00014302 65.00007666 65.00004109 65.00002203 65.00001181\n",
      " 65.00000633 65.00000339 65.00000182 65.00000097 65.00000052 65.00000028\n",
      " 65.00000015 65.00000008 65.00000004 65.00000002 65.00000001 65.00000001\n",
      " 65.         65.         65.         65.         65.         65.\n",
      " 65.         65.         65.         65.         65.         65.\n",
      " 65.         65.        ]\n"
     ]
    }
   ],
   "source": [
    "T_analytical = np.zeros(len(t))\n",
    "\n",
    "T_analytical = T_a + ((T_o - T_a)*np.exp(-K*t))\n",
    "T_analytical[0] = T_o\n",
    "\n",
    "print(T_analytical)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU5dn/8c83+woECAoETKxWRRGUqBSXSlGL1oUuWKm2alvpplZ9/KnVWm37ax/rY+2jj1YfxLW1WrXy01a0blVrLSooiigKIkIE2ffsyfX745yEIUySE8hkwsz1fr3GmbnnLNcJL+eaezn3LTPDOeecaysj2QE455zrnTxBOOeci8sThHPOubg8QTjnnIvLE4Rzzrm4PEE455yLyxOESyuSjpb0fjccZ4mk47ojpmTYlfgl3S7p6l08/7GSqnblGC7xPEG4bhV+8ayUVBhT9l1JLyQxrFZm9k8z2y/ZcUQlqUJSs6TfJ+n850h6ObbMzL5vZr9MRjyuZ3mCcImQBfw42UG0JSkr2THshG8B64EzJOUmOxiXXjxBuET4L+BSSf3afiCpXJLFfllLekHSd8PX50j6l6TfSdogabGkcWH5MkmrJJ0ds2+upBskLQ1rLrdLyg8/O1ZSlaTLJX0K3N22aUPSMEmPSlotaa2kW8Lyz0h6PixbI+n+eNcT5/rGSvpUUmZM2ZclvR2+PlzSbEmbwnhv7OSQ3wJ+CjQAp7Q5l0n6vqSFktZLulWSuhK/pD0lVUsaEFM2Jvx7jARuBz4naYukDeHn90j6vzHbnyZpbnhNH0qaGJafK+k9SZvDf8fvdfb3c72LJwiXCLOBF4BLd3L/I4C3gQHAn4AHgcOAfYCzgFskFYXb/gb4LDA6/Hwo8LOYY+0J9Af2AqbGniT8Ev8b8DFQHu77YMvHwH8CQ4ADgGHAtZ0FbmazgK3AF2KKvxFeB8BNwE1m1gf4DPBQe8eSdDRQFsb0EEGyaOtkgr/NKOB04Itdid/MPiX4tzo9pvgs4EEzmwd8H/i3mRWZWbwEczhwH/B/gH7AMcCS8ONVYXx9gHOB30k6tL3rdb2PJwiXKD8DLpBUuhP7fmRmd5tZE/Bngi+3X5hZnZk9DdQD+4S/ls8DLjazdWa2Gfg1cEbMsZqBa8J9a9qc53CCL9D/Y2ZbzazWzF4GMLNFZvZMuN9q4Ebg8xHjfwCYAiCpGDgpLIOgJrCPpIFmtiVMKO05G3jSzNYTJJgTJQ1qs811ZrbBzJYC/yBIlF2N/16CpNCSNKcAf4h4rd8B7grP1Wxmn5jZgjCGJ8zsQwu8CDwNHB3xuK4X8AThEsLM3iH4dX7FTuy+MuZ1TXi8tmVFQClQAMwJm6M2AE+F5S1Wm1ltO+cZBnxsZo1tP5A0SNKDkj6RtAn4IzAwYvx/Ar4S9hl8BXjDzD4OP/sOQY1ngaTXJZ0c7wBhM9lk4H4AM/s3sJSgNhLr05jX1QR/l67G/xgwQtLewPHARjN7LeK1DgM+bOcaTpQ0S9K68N/mpA5icL2QJwiXSNcQ/MIfGlO2NXwuiCnbcyePv4YgWRxoZv3CR18zK4rZpqPpipcBw9vpvP7PcN+Dw+agswiabTplZu8SNFudyPbNS5jZQjObAgwiaB57JHbEV4wvEzTN/D7s0/iU4O8Yr5kpnsjxhwn0IeBM4JtsX3vobLrnZQRNZdsJk+NfgBuAPcLmqZntxeB6J08QLmHMbBFBE9GFMWWrgU+AsyRlSvo2cb5gIh6/GbiDoG17EICkoZK+2PGerV4DVgDXSSqUlCfpyPCzYmALsEHSUII29q74E8F1HwM83FIo6SxJpWHsG8Lipjj7nw3cBYwkaDYaDRwJjA47jzvT1fjvA84BTiWobbRYCZRJymlnvzuBcyVNkJQR/v33B3KAXGA10CjpROCECHG7XsQThEu0XwBtfyGfR/CFtRY4EHhlF45/ObAImBU2pTwLRLrPIezjOIWgc3spUAV8Pfz458ChwEbgCeDRLsb1AHAs8LyZrYkpnwjMl7SFoMP6jLZNYOEX+gTgv83s05jHHIImtLPpXJfiN7N/EfTXvGFmS2I+eh6YD3wqaU2c/V4j7IAOz/UisFfYH3QhQc1kPUFN6vEIcbteRL5gkHMOQNLzwJ/MbHqyY3G9gycI5xySDgOeAYaFv/6d8yYm59KdpHsJmuYu8uTgYnkNwjnnXFxeg3DOORfX7jh5WbsGDhxo5eXlyQ7DOed2G3PmzFljZnFnPEipBFFeXs7s2bOTHYZzzu02JH3c3mfexOSccy4uTxDOOefi8gThnHMurpTqg3DOpYaGhgaqqqqorW1vIl7XVXl5eZSVlZGdnR15H08Qzrlep6qqiuLiYsrLywkXyXO7wMxYu3YtVVVVVFRURN4voU1Mki6WNF/SO5IeCGfLvDaco35u+DipnX0nSnpf0iJJO7OmgHNuN1VbW8uAAQM8OXQTSQwYMKDLNbKEJYhwRsoLgUozOwjIZNtKX78zs9HhY2acfTOBWwnm0x8BTJE0IlGxOud6H08O3Wtn/p6J7qTOAvLDBVkKgOUR9zscWGRmi82snmBN3tMSEaCZ8e+7L+ftF/6SiMM759xuK2EJwsw+IVhNainBoiwbw/WEAc6X9LakuySVxNl9KMFKVS2q2H5VslaSpkqaLWn26tWruxynJEYuuZet7/69y/s651xn7rnnHs4///xOt1m+fNvv5+9+97u8++67XT7XCy+8wMknx13FdqcksomphOBXfwXBwvCFks4CbiNYQWw0QeL4bbzd45TFnVXQzKaZWaWZVZaWxr1bvFO1ykP1W3ZqX+ec21VtE8T06dMZMSL5reqJbGI6DvjIzFabWQPBilbjzGylmTXFLBd5eJx9qwgWQ29RRvTmqS6rzSggo35r5xs659LKpEmTGDNmDAceeCDTpk0DoKioiKuuuopRo0YxduxYVq5cCcBf//pXjjjiCA455BCOO+641vIWmzdvpqKigoaGBgA2bdpEeXk5Dz/8MLNnz+bMM89k9OjR1NTUcOyxx7ZOG/TUU09x6KGHMmrUKCZMmADAa6+9xrhx4zjkkEMYN24c77//fkKuP5HDXJcCYyUVECwsPwGYLWmwma0It/ky8E6cfV8H9pVUQbB+8RkESxYmRH1mAZmNniCc641+/tf5vLt8U7cec8SQPlxzyoGdbnfXXXfRv39/ampqOOyww/jqV7/K1q1bGTt2LL/61a+47LLLuOOOO/jpT3/KUUcdxaxZs5DE9OnTuf766/ntb7c1kBQXF3PsscfyxBNPMGnSJB588EG++tWvMnnyZG699VZuuOEGKisrtzv/6tWrOe+883jppZeoqKhg3bp1AOy///689NJLZGVl8eyzz3LllVfyl790fz9qwhKEmb0q6RHgDaAReBOYBkyXNJqgyWgJ8D0ASUOA6WZ2kpk1Sjof+DvB6Ke7zGx+omJtyCok22sQzrk2br75ZmbMmAHAsmXLWLhwITk5Oa3t/GPGjOGZZ54Bgns3vv71r7NixQrq6+vj3m/w3e9+l+uvv55JkyZx9913c8cdd3R4/lmzZnHMMce0Hqt///4AbNy4kbPPPpuFCxciqbVW0t0SeqOcmV0DXNOm+JvtbLscOCnm/UxghyGwidCUVUhuzcaeOJVzroui/NJPhBdeeIFnn32Wf//73xQUFHDsscdSW1tLdnZ265DRzMxMGhsbAbjgggu45JJLOPXUU3nhhRe49tprdzjmkUceyZIlS3jxxRdpamrioIMO6jAGM4s7PPXqq69m/PjxzJgxgyVLlnDsscfu8vXG43MxAZZTSJ5VJzsM51wvsnHjRkpKSigoKGDBggXMmjWr0+2HDg0GW957773tbvetb32LKVOmcO6557aWFRcXs3nzjqu9fu5zn+PFF1/ko48+AmhtYoo91z333NOl6+oKTxAAOcUUUENtQ1OyI3HO9RITJ06ksbGRgw8+mKuvvpqxY8d2uP21117L5MmTOfrooxk4cGC725155pmsX7+eKVOmtJadc845fP/732/tpG5RWlrKtGnT+MpXvsKoUaP4+te/DsBll13GT37yE4488kiamhL3vZVSa1JXVlbaziwY9O69F1Kx+AE2XrKMPfvmJSAy51xXvPfeexxwwAHJDiMhHnnkER577DH+8Ic/9Pi54/1dJc0xs8p42/tkfUBmXh/yVc/SrTWeIJxzCXPBBRfw5JNPMnNmj3Sv7jJPEEB2fh8AtmzeCMS7sds553bd//zP/yQ7hC7xPgggp6AYgC2bNiQ5Euec6z08QQB5RX0BqNniCcI551p4ggDywwRRt7V779Z0zrndmScIIL8wSBD11X6znHPOtfAEAWTkBX0QDTU73qjinHPJ8Pjjj3Pdddft1L7l5eWsWbNml2PwUUwAOUUANNV6gnDOJV9jYyOnnnoqp556alLj6LQGocBISV+UdIykAT0RWI/KDWoQzZ4gnHOhJUuWcMABB3Deeedx4IEHcsIJJ+wwFfeaNWsoLy8HgikvJk2axCmnnEJFRQW33HILN954I4cccghjx45tnSbjww8/ZOLEiYwZM4ajjz6aBQsWAMHd1Jdccgnjx4/n8ssv326hoZUrV/LlL3+ZUaNGMWrUKF555RUg/nTk3andGoSkcuAyYCLwEbAayCOYhnsDcDvwR0uFW7HDGgS+aJBzvc+TV8Cn87r3mHuOhBM7b75ZuHAhDzzwAHfccQenn356p1Nqv/POO7z55pvU1tayzz778Jvf/IY333yTiy++mPvuu4+LLrqIqVOncvvtt7Pvvvvy6quv8sMf/pDnn38egA8++IBnn32WzMzM7eZYuvDCC/n85z/PjBkzaGpqYsuW4Lsq3nTkAwZ032/4jpqYridY/e38cHGfVuHU3N8Azgbu2XHX3Ux2Ps1kkOkJwjkXo6KigtGjRwPB1N5LlizpcPvx48dTXFxMcXExffv25ZRTTgFg5MiRvP3222zZsoVXXnmFyZMnt+5TV1fX+nry5MlkZmbucNznn3+e++67DwhmkO3bNxhYE2868h5JEGZ2OoCkLKC5zcerzOyGbosi2STqMgrI8EWDnOt9IvzST5Tc3NzW15mZmdTU1JCVlUVzc/CVWFtb2+72GRkZre8zMjJobGykubmZfv36MXfu3LjnKywsjBxbe9ORd6coo5hei1i2W2vMyiersZpUaDFzziVOeXk5c+bMAYKJ97qiT58+VFRU8PDDDwPBeg9vvfVWp/tNmDCB2267DYCmpiY2bdrU5enId0a7CULSIEmjgPywk/rg8HEUUBDl4JIuljRf0juSHpCUJ+m/JC2Q9LakGZL6tbPvEknzJM2V1PUpWruoMauQAmrYWu9Tfjvn2nfppZdy2223MW7cuJ0aSnr//fdz5513MmrUKA488EAee+yxTve56aab+Mc//sHIkSMZM2YM8+fP7/J05Duj3em+JZ0LfBsYDcTWhzYDd5vZwx0eWBoKvAyMMLMaSQ8RrBC3HHg+XFb0NwBmdnmc/ZcAlWYW+V9gZ6f7Blj7uyOZty6TfS55irKSSPnPOZcgqTzddzJ123TfZnY3cLek083soZ2MJ4ugBtJAUOtYbmZPx3w+C/jaTh67e+UUUah1bKxpoMwndHXOuQ6HuU4xsweAwZIubPu5md3c0YHN7BNJNwBLgRrg6TbJAYIayp/bOwTwtCQD/tfM4g7ylTQVmAowfPjwjkLqkHKLKWI562sSs/i3c87tbjrqpG75HT0QKI3z6JCkEuA0oAIYAhRKOivm86uARuD+dg5xpJkdCpwI/EjSMfE2MrNpZlZpZpWlpZ2G1a6MvGIKqWFjtScI53oDHzDSvXbm79lRE9Pvw+erdzKe44CPzGw1gKRHgXHAHyWdDZwMTGjvRjszWx4+r5I0AzgceGknY+lUdkExhaplo9cgnEu6vLw81q5dy4ABA5CU7HB2e2bG2rVrycvr2oqZnc7FJGkf4FZgTzMbJelg4Etm9p+d7LoUGCupgKCJaQIwW9JE4HLg82ZW3c45C4EMM9scvj4B+EXkq9oJ2QV9yMIThHO9QVlZGVVVVaxevTrZoaSMvLw8ysrKurRPlMn6pgNXEiQJgHnAA0CHCcLMXpX0CPAGQVPSm8A0YD6QCzwT/jKYZWbfD+/Onm5mJwF7ADPCz7OAP5nZU126si7Kzu+D1MCmrTWJPI1zLoLs7GwqKiqSHUbai5IgCs3slZZqnplZOCqpU2Z2DXBNm+J92tl2OXBS+HoxMCrKObqLwvmYarf6mhDOOQfR7qReK6mCYFQRkiYBnyY0qmTIDRJEXbWvKueccxCtBnE+cCewv6SPgRXAlIRGlQxhDaLBE4RzzgEREoSZLQK+IKkvwZ3XGxIfVhKEa0I0+apyzjkHRBvFlANMAsqBrJi+iF8nNLKeFtYgmus8QTjnHERrYpoB1AJzgNSdyS7sgzBPEM45B0RLEHuZ2UEJjyTZwhpERsNWmpuNjAy/Occ5l96ijGKaJWlEwiNJtrAPopAaNtc2JjkY55xLvo4m63uTYGhrNjBX0kKgDhDB7RCH9kyIPSSsQRRRy4aaevoWZCc5IOecS66Omph6xzTcPSUrl2ZlUagan27DOefoeLK+D1teKxi6VNrR9rs9iebsAgobfD4m55yDaMNcf0gwUd5aoDksNiDl+iUsp4iimlo2+JTfzjkXqUZwCXBAy7TdqUy5xRRQyzqvQTjnXKRRTFXAukQH0htk5BZRhPdBOOccRKtBLAKel/Q3glFMQOdLju6OMvKKKcr4xBOEc84RLUGsCB99EhxL8uUU0Sejzpcddc45ok3WdzWApFwzq+ts+91abjFFvqqcc84BEfogJB0uaR6wMHw/StL/RDm4pIslzZf0jqQHJOVJ6i/pGUkLw+eSdvadKOl9SYskXdGlq9pZOUUUUsOGmvoeOZ1zzvVmUTqpbwZOJhjmipm9BYzvbCdJQ4ELgcpwLqdM4AzgCuA5M9sXeC5833bfTIIlTk8kGE47pUem+8gtIt9q2FjjU20451yUBJFhZh+3KYs6q2sWkC8pCygAlgOnAfeGn99LMJV4W4cDi8xssZnVAw+G+yVWThFZNFJTXZ3wUznnXG8XJUEsk3Q4YJIyJV0EfNDZTmb2CXADsJSgk3ujmT0N7GFmK8JtVgCD4uw+FFgW874qLNuBpKmSZkuavXr1Lt6q0bKqXI2vS+2cc1ESxA8IbpYbDqwExoZlHQr7Fk4DKoAhQKGksyLGFW+ubYu3oZlNM7NKM6ssLS2NePh2hGtCqGErDU3NnWzsnHOprcNRTGFfwBlmdsZOHPs44KOWO7AlPQqMA1ZKGmxmKyQNBlbF2bcKGBbzvoygeSqxwhpEIbVsqmlgQFFuwk/pnHO9VYc1CDNrAr66k8deCoyVVBBO9jcBeA94HDg73OZs4LE4+74O7CupIlzy9Ixwv8TK3ZYgNvhQV+dcmotyo9w/Jd1E0FG8taXQzN7uaCcze1XSI8AbQCPwJjANKAIekvQdgiQyGUDSEGC6mZ1kZo2Szgf+TjD66S4zm9/lq+uqnGDRoCKf8ts55yIliM+Hz7ELBBlwTGc7mtk1wDVtiusIahNtt10OnBTzfiYwM0J83SemBuEJwjmX7qLcSX10TwTSK7SsKqcan27DOZf2Ii0AJOmLwIFAXkuZmf06UUElTeu61F6DcM65KAsG/R7oR9CkdDdBp/WsBMeVHDnexOSccy2i3AdxlJl9A1gbTtx3BMGw09STlQOZOZRk1fmqcs65tBclQdSEz7WS9gRqgfKERZRsOYWUZNV5DcI5l/ai9EE8KakfwbQZcwnmYbq34112YznF9G3yBOGcc1FGMV0bvnw4XFUu38xSdwnS3CL61NWx0af8ds6luXabmCT9Mub1FwDMrCalkwNAThHF8hqEc8511AfxpZjXNyQ6kF4jN1g0yBOEcy7dRemkTi85ReRT46OYnHNpr6M+iEGSLiSYervldSszuzmhkSVLbjF5zdXUNTZT29BEXnZmsiNyzrmk6ChB3A2Uxnmd2nKKyGkKVpTbVNPgCcI5l7baTRDhTXHpJ7eI7KZqwNhQ08CgPnmd7uKcc6nI+yDayikiw5rIpcE7qp1zac0TRFstM7riM7o659JbpNlcd4ak/YA/xxTtDfwM+BywX1jWD9hgZqPj7L8E2Exw53ajmVUmKtbttKwJIV9VzjmX3qLM5no+cJ+ZbZL0v8AhwE/M7LmO9jOz94HR4TEygU+AGWb23zHH/i2wsYPDjDezNZ1fRjeKrUF4gnDOpbEoTUxTw+RwAjAU+AFwfRfPMwH40Mw+bikI16k+HXigi8dKrJgahCcI51w6i5IgLHw+EbjbzOZE3C/WGeyYCI4GVprZwg7O+7SkOZKmtndgSVMlzZY0e/Xq1V0MK45wXerSnAY2Vvt8TM659BXli/4tSTOBUwhmdi1iW9LolKQc4FTg4TYfTaHj2sORZnYoQWL6kaS4a2Cb2TQzqzSzytLSbrhVI6xBlOb4KCbnXHqL0kl9LjAGWGRm1ZIGAt/pwjlOBN4ws5UtBZKygK+Ex43LzJaHz6skzQAOB17qwnl3TtgHMSC7nqWeIJxzaazTGoSZNRGMQPpBWJQfZb8Y8WoKxwELzKwq3g6SCiUVt7wGTgDe6cI5d15YgyjJqvdRTM65tNbpF72kW4DxwFlh0Vbg9igHl1QAHA882uajHfokJA0Jm7IA9gBelvQW8BrwhJk9FeWcuyzsg+jnq8o559JclCamcWZ2qKQ3AcxsXdiv0CkzqwYGxCk/J07ZcuCk8PViYFSUc3S7zCzIyqNvRh2bPEE459JYlKaiBkkZhB3TkgYAzQmNKtlyCilWLRuqGzCL3B/vnHMpJUqCuBX4C1Aq6efAy8BvEhpVsuUUUaQaGpuN6vqmZEfjnHNJEWVN6vskzSHoWBYw2cx6psM4WXKLybcaADbWNFCYm7AZSZxzrteKOhqpGFgfTpOxQtLwBMaUfDlF2yUI55xLR1FGMf0UuAb4aViUB/wpkUElXW4ROc1BgvClR51z6SpKDeJrBKOLtgKY2SdAn0QGlXQxq8p5DcI5l66iJIg6C4bytIxiKkhsSL1AbhGZjVsBfKircy5tRUkQj0q6Fegr6VzgaeCuxIaVZDnFZDZsAWBDjU/Y55xLT1FGMf1G0olAPcHNa78ysycTHlky5RZB/VYyM7yJyTmXvjpMEOFCPzPN7ItAaieFWDmFyJrZI6/ZE4RzLm112MQUTtRXLym1O6XbCmd0HZzX6KOYnHNpK8odYFsI1oR4mnAkE4CZXZKwqJItN5iwb1CurwnhnEtfURLEs+EjfYQ1iKEFjby0qTbJwTjnXHJE6aS+sycC6VXCNSGGFTZTtawGMyNYQts559JHpwkinOa77ZSmG4HZwH+a2bpEBJZU4ZoQgwuaqK5vYn11A/0LI81w7pxzKSNKE9MzBJP0tUyvcQbBdN+bgXsI1ptOLWENYs+8RgCWrav2BOGcSztRFww6Kub9m5JeNrOjJM1rbydJ+wF/jinaG/gZ0A84D1gdll9pZjPb7I6kicBNQCYw3cyuixBr9wj7IEpzgg7qqvU1jBrWr8dO75xzvUGUO6mLJY1peSPpULbNxdTY3k5m9r6ZjTaz0cAYoBqYEX78u5bP2kkOmQTrUJwIjACmSBoR6Yq6Q+u61HUAVK2v7rFTO+dcbxGlBvE94A+SsgmamuqB70gqBK6PeJ4JwIdm9nHEzt7DgUXh0qNIehA4DXg34vl2TViDyGuuoV9BNss8QTjn0lCUUUyzgBHhUqMyszUxHz8Q8TxntNn2fEnfIujo/g8zW99m+6HAspj3VcAR8Q4saSowFWD48G5apiIjE7LyoW4zZSX5VK2v6Z7jOufcbiTKehClkv4XuNfM1kgaIemcqCeQlEPQkf1wWHQb8BlgNLAC+G283eKUxV0c2symmVmlmVWWlpZGDatzuUVQv4WyfgWeIJxzaSlKH8Q9wIvAsPD9QuA/unCOE4E3zGwlgJmtNLMmM2sG7iBoTmqrKuZ8AGXA8i6cc9flFEHdFob1z6dqfTXBjOfOOZc+oiSIQWb2J4KhrZhZA9DUhXNMIaZ5SdLgmM++DMRb3/p1YF9JFWEN5Azg8S6cc9e11CBKCqhtaGbNFp/22zmXXqIkiK2S+rNtwaDDCO6B6FS4uNDxwKMxxddLmifpbWA8cHG47RBJMwHMrBE4H/g78B7wkJnNj3ZJ3SSnGOq2UFaSD+Ad1c65tBNlFNOlwF+BvSW9SNCB/LUoBzezamBAm7JvtrPtcoKlTVvezwR2GALbY3KLYMtKhvUPFtCrWl/DocNLkhaOc871tCijmGZLGg8cQNB5/K6ZpX57S04R1C9maL+gBuH3Qjjn0k27CUJSe1NoDJeEmfVsn0BPyw06qQtzs+hfmMOydT6SyTmXXjqqQUwOnwcC44B/ENQgPk8wqim1E0ROMdQH61IPK8n3GoRzLu20myBa+gokPQ6MMLNPwvdDgZt7JrwkCkcx0dxMWUkB763YlOyInHOuR0UZxbR3S3IILQf2S1A8vUdOYfDcsLX1burmZr8XwjmXPqKMYnpJ0hME9zIYwT0JLyU0qt4gnI+Jui2U9S+gvqmZ1Vvq2KNPXnLjcs65HhIlQfyIYFjrMeH7+4BHEhZRbxGuSx3cLNcXCEYyeYJwzqWLKMNcjWAepYc72zaltNYgNjOsZE8Alq2rYcxeSYzJOed6ULt9EJL+IekHkoa0Kc+SdIykOyWdm/gQkyRcE6Jlug3weyGcc+mloxrEl4DvAjPCkUvrgHwgF3gOuNXMZic+xCSJ6YPIy85kYFGuz+rqnEsrHQ1zrSYYznqzpFxgEFDTZj2I1NXaB7EVgLKSfJ+PyTmXVqIMc8XM6sxsWdokB9hWg6gP5iUc1t/XhXDOpZdICSIt5W5rYoKgBrF8Qw1Nfi+Ecy5NeIJoT3Z4o1zrdBsFNDQZKzfVJjEo55zrOZEShKSycEZXJOVKKkxsWL1ARkaQJGJqEIA3Mznn0kaUNam/TTAx3/SwaC/gsUQG1WvkFrX2QbQuHLTOO6qdc+khyp3UFxKsG/0qgJl9IPumU0AAABhzSURBVGlQZztJ2g/4c0zR3sDPCBYcOgWoBz4EzjWzDXH2X0Kwcl0T0GhmlRFi7V7hutQAQ70G4ZxLM1GamGpjFwiSlEkw7XeHzOx9MxttZqOBMUA1MAN4BjjIzA4GPgB+0sFhxofH6PnkANtmdAVyszLZo0+u3yznnEsbURLEvyRdBuSF/RB/Bv7WxfNMAD40s4/N7OlwzWmAWUBZF4/Vc3L7QO3G1rdlJQV+L4RzLm1ESRCXETT1LAB+THAX9VVdPM8ZBLPBtvVt4Ml29jHgaUlzJE1t78CSpkqaLWn26tWruxhWJ/oNhw1LW98OC6f9ds65dNBhH0TYnHSXmZ0N3LYzJ5CUA5xKm6YkSVcBjcD97ex6pJktD/s7npG0wMx2mGbczKYB0wAqKyu79yaFknLYvAIaaiA7n7KSAv769goam5rJyvQRws651Nbht5yZNQGDJWXvwjlOBN4ws5UtBZLOBk4Gzgxni4137uXh8yqCvovDdyGGnVNSHjyHtYiyknyamo0VG/1eCOdc6osyimkx8E9JjwFbWwrNLOqyo1OIaV6SNBG4HPh8ON/TDsL7LDLMbHP4+gTgFxHP131KKoLndR9B6X4M698yq2tN62vnnEtVUdpJVhOMPCoASmMenZJUABwPPBpTfAtQTNBsNFfS7eG2QyTNDLfZA3hZ0lvAa8ATZvZUlHN2q5YaxPolQMy9EN5R7ZxLA1EWDLp6Zw8e1hAGtCnbp51tlwMnha8XA6N29rzdpnBgcDd1mCAG981H8nshnHPpodMEIekZghFF2zGzExISUW8iQf+K1gSRk5XB4D55fi+Ecy4tROmD+GnM6zzgq0BdYsLphUrKYe2i1rdlJQVUrfMahHMu9UVpYnq1TdGLkl5MUDy9T0k5LHoWzECirCSfWYvXJjsq55xLuCiT9fWJefSTNAEY3AOx9Q4l5dBYC1uCUbpl/Qv4dFMt9Y3NyY3LOecSLEoT03yCPggR3Nj2EXBeIoPqVWKHuhbvSVlJPs0GKzbWsNeA1J/13DmXvqIkiL3NrCG2QFKU/VJD7FDXvT633boQniCcc6ksyn0QbfsgILg3IT30GwaodSTTsJKWm+V8JJNzLrW1WxMI50AaDORLGsm2Kb77ENw0lx6ycqHP0Jh7IfLIzBDLfCSTcy7FddRU9CWC2VbLgN/HlG8Gdvrmud1S/wpY/xEAWZkZDO7r90I451JfuwnCzO4G7pZ0upk91IMx9T4le8HCZ1rflpXks8zvpnbOpbgo90E8JOmLwIEEN8q1lP86kYH1KiXlwTDX+mrIKaCspICXPujmtSecc66XiXIfxO+Bs4FLgHzgLCDufEopq2Wo64aPAfhMaRGrNtexdkv63FDunEs/UUYxHWVm3wDWhhP3HUFvXiY0EWLvhQAqy0sAmP3x+mRF5JxzCRclQbSsjlMrac/wfXnCIuqN2kz7fXBZX3KyMpi9ZF3SQnLOuUSLcsPbTEn9gBuAuUATcG9Co+ptCvpDTnFrgsjNymRUWV9eX+I1COdc6uqwBiEpA3jSzDaY2cNABTDSzK7s7MCS9gsXBGp5bJJ0kaT+kp6RtDB8Lmln/4mS3pe0SNIVO3V13UWC/uWtCQKgsrw/73yyker6xqSF5ZxzidTZmtTNwE0x72vMLFK7ipm9b2ajzWw0MAaoJlhb+grgOTPbF3gufL8dSZnArQTrWY8ApkgaEe2SEqSkvPVeCIDDy/vT2GzMXbYheTE551wCRemDeEbSabt4ngnAh2b2MXAa25qo7gUmxdn+cGCRmS02s3rgwXC/5Ckph/UfQ3Mwi+uhw0uQYLY3MznnUlSUPojzgb6S6oAagik3zMz6d+E8ZwAPhK/3MLMVBAdZEU7p0dZQYFnM+yqC0VM7kDQVmAowfPjwLoTURSXl0FQHWz6FPkPoW5DNfnsU87p3VDvnUlSUGsRAIBsoAkrD96VRTyApBzgVeLgLcSlO2Q7LngKY2TQzqzSzytLSyGF1XctQ15h+iMPK+/PGx+tpbPK1IZxzqafTBGFmTcBk4PLw9WBgdBfOcSLwhpmtDN+vlDQYIHxeFWefKmBYzPsyYHkXztn9Woa6rtvWD1FZXsLW+iYWfLo5OTE551wCRbmT+hZgPPDNsKgauL0L55jCtuYlgMcJ7swmfH4szj6vA/tKqghrIGeE+yVP32GgjO1qEIdXBK1sr33kzUzOudQTpYlpnJl9j/CGuXAUU06Ug0sqAI4HHo0pvg44XtLC8LPrwm2HSJoZnqORoO/j78B7wENmNj/SFSVKVg70KdsuQQzum8/QfvnM/tgThHMu9UTppG4I74cwAEkDgEiN7mZWDQxoU7aWYFRT222XAyfFvJ8JzIxynh7Tv3y7oa4Ah5WX8K8P12JmSPG6TpxzbvcUpQZxK/AXoFTSz4GXgd8kNKreqqR8uxoEwGEV/Vm9uY6P1/r6EM651BJluu/7JM0BjguLJpvZO4kNq5cqKYetq6FuC+QWAcFIJoDXl6yjfKCvUe2cSx1RahAAmUADUN+FfVJPy0imcNpvgH1Ki+ibn+03zDnnUk6UUUxXEYxCGkIw3PRPkn6S6MB6pTbTfgNkZIjDykv8hjnnXMqJUhs4CzjMzH5qZlcRTIPxrcSG1Uu1mfa7RWV5fxav2coaX0DIOZdCoiSIj9m+ryILWJyYcHq5/BLI7btjR3XYD+HNTM65VBIlQVQD8yVNl3QHMA/YIOlGSTcmNrxeJs603wAHDe1Dri8g5JxLMVHug3gifLSYlaBYdg8l5bBy+3v2crMyGTWsn/dDOOdSSpRhrnf2RCC7jZJyeP9JaG6CjMzW4sPL+3Pbix9SXd9IQU6UvOucc71blFFMEyW9LmmVpHWS1ktK35/KJeXQVA+bV2xXXFleQlOzMXepLyDknEsNUfogbgG+R7BGQ5en+045cab9Bjh0r2ABode8mck5lyKiJIgqYK6ZNZhZU8sj0YH1WnGm/Qbok5fNAXv28ZFMzrmUEaWx/DLgr5JeAFoH+pvZzYkKqlfrWwbK3KEGAcHEfQ/PqaKxqZmszPS94dw5lxqifIv9HGgC+hE0LbU80lNmdpAk4iSII/YeQHV9E7MWezOTc273F6UGMcjMxiQ8kt1J/4q4CeIL+w9iQGEO97zyEUftO7Dn43LOuW4UpQbxnKQvJDyS3UlJOaz7EGz7ZbLzsjM5c+xePLdgFR+t2Zqc2JxzrptESRDnAc9K2tLVYa6S+kl6RNICSe9J+pykP0uaGz6WSJrbzr5LJM0Lt5vdlYtKuGFHQM16qHp9h4/OGjuc7IwM7v7XR3F2dM653UeUBDEQyAb60vVhrjcBT5nZ/sAo4D0z+7qZjTaz0QQLET3awf7jw20rI56vZ+x/MmTlwdsP7fDRoOI8Thk1hIdnV7GxuiEJwTnnXPfoNEGEQ1onA5eHrwcDozvbT1If4BjgzvA49Wa2IeZzAacTTCW+e8nrA5+dCPMfhaYdk8C3jyqnpqGJB19fmoTgnHOue0S5k/oWYDzwzbCoGrg9wrH3BlYDd0t6M5zsL3bJtaOBlWa2sJ39DXha0hxJUzuIb6qk2ZJmr169OkJY3eTg06F6LXz4jx0+OnBIX8bu3Z97X1lCY1Ok5budc67XidLENM7MvgfUApjZOiAnwn5ZwKHAbWZ2CLAVuCLm8yl0XHs40swOBU4EfiTpmHgbmdk0M6s0s8rS0h4cfbvP8ZDXD+bt2MwE8J2j9mb5xlqemv9pz8XknHPdKEqCaJCUQfCLHkkDgCg/i6uAKjN7NXz/CEHCQFIW8BXgz+3tbGbLw+dVwAyChYp6j6wcOHASLHgiWKO6jQn7D2KvAQXc+bJ3Vjvndk/tJojwSxzgVoLO5FJJPwdeBn7T2YHN7FNgmaT9wqIJwLvh6+OABWZW1c65CyUVt7wGTgDe6fxyetjBX4eGanh/5g4fZWSIc8eV8+bSDbyx1KffcM7tfjqqQbwGYGb3AT8FbgDWA5PN7MGIx78AuF/S2wQd278Oy8+gTfOSpCGSWr5p9wBelvRWGMcTZvZUxHP2nGFjoe+wuKOZACZXDqM4L4u7vBbhnNsNdXQntVpemNl8YH4H28ZlZnOBHYaomtk5ccqWAyeFrxcTDIvt3TIyYOTX4F83w5bVULR9H0hhbhZnHDaMu/61hOUbahjSLz9JgTrnXNd1VIMolXRJe48ei7C3G3k6WFMw5DWOs8eVY2bc++8lPRqWc87tqo4SRCZQBBS383AAe4yAPQ5qt5mprKSAiQftyQOvLmVrXWMPB+ecczuvoyamFWb2ix6LZHc2cjI8ew2s/RAGfGaHj79zVAUz533Kzc8t5CcnHZCEAJ1zrus6qkGog89crJFfAwTzHon78Zi9+vONI4bzvy8t5rG5n/RsbM45t5M6ShATeiyK3V3fMtjryOCmuTYzvLa49pQDObyiP5c98jZvLfN1q51zvV+7CSK8Y9pFdfBkWLsIlr8Z9+OcrAxuO/NQBhblMvUPs1m1qbaHA3TOua7xdTG7y4jTIDMH5j3c7iYDinKZfnYlm2sbmfqHOdQ2pO/S3s653s8TRHfJL4F9T4B3/gJN7Y9WOmBwH248fRRzl23gyhnzsHaapJxzLtk8QXSn0WfClpXw1BXt9kUATDxoMBcf91kefeMTn6vJOddreYLoTvudCOMugNfvgOc6HiF8wRf24aSRe/Lrme/x17eW91CAzjkXXUf3QbiukuD4X0LdZnj5RsgtgqP/I+6mGRnihsmjWL6hlgseeJNn3l3JNaeMYEBRbg8H7Zxz8XkNortJ8KUbg5vnnvsFvHZHu5sW5GTx0Pc+xyXHf5Yn31nB8b97icffWu79Es65XsETRCJkZMKk22C/k2DmpTC3/XWRcrIyuHDCvvztgqMZ1r+ACx94k/Pum8NKHwbrnEsyTxCJkpkNX7sbKo6Bx34I7z7e4eb77VnMoz8Yx1UnHcA/F67muBtf5M6XP2LtlroeCtg557anVGrOqKystNmzZyc7jO3VbYE/TILlc+HIC+HQb0FJeYe7fLRmKz959G1mLV5HZoY4ap+BnDZ6CCccuCdFud5t5JzrPpLmmNkOyzKAJ4ieUbMeHjs/WHnODD7zBRhzTjDqKTM77i5mxoJPN/P4W8t5fO5yPtlQQ25WBscdsAcTD9qTUWX9GNY/H8mnzHLO7bykJQhJ/YDpwEEEa1p/G/gicB6wOtzsSjPbYc1OSROBmwimHZ9uZtd1dr5emyBabKyCN/8Ib9wHmz6BwkFwyFmw7/EwYB8oLA06udswM95Yup7H5y7nb2+vYO3WegCKc7M4YHAfRgwJHgfs2Ych/fIoKcghI8MTh3Ouc8lMEPcC/zSz6ZJygALgImCLmd3QwX6ZwAfA8UAV8DowxczebW8f2A0SRIvmJlj0LMy5Bz54Cqw5KM/tE0wXPmDfIGH0HRqU5fWB3GLI7UtjdiHvrTXmr6ph/opq3v10M++t2ER1/bZpO7IyxKDiXAb1yQufc+mXn0NRXhZFuVkU52VRmJNFUV4W+dmZ5GZnkJOZQW52JjmZGeRkZZCblUFmhsiUPNk4l8I6ShAJa9CW1Ac4BjgHwMzqgfqITSKHA4vCpUeR9CBwGtBhgthtZGTCZ78YPDavhJXzYM2iYLK/tQth6b+DmWHjyAJGho/gWFlYbg7NeVk0kUkTGTSZaGwSjeugcY1oNGgyMAv+9gYYwmJmdDegFqjZrmzb65Z/NrX+p+UpZpsdXrRbsN0xnXO7pjqzLyOu+le3HzeRPZ57EzQj3S1pFDAH+HH42fmSvgXMBv7DzNa32XcosCzmfRVwRLyTSJoKTAUYPnx490XfU4r3CB77HLd9eX01bF0V3HRXuyl4rtsMdRuDju/mhmDOp6Z61NxAZlMjmU31QW2k9WHhcxNmRlOz0djURGNTc/hoosmMpmZoNqO5uZlma3kdpJBmC5q4zIIk0txS47SWRGOtr1ufadnE4k45Yq3/ia939or1zqicA2jM7pOQ4yYyQWQBhwIXmNmrkm4CrgBuAX5J8H/cL4HfEvRNxIr32zLu/6FmNg2YBkETU/eE3gvkFEBOebcdTgT/ID4GyjkXVSLvg6gCqszs1fD9I8ChZrbSzJrMrBm4g6A5Kd6+w2LelwE+YZFzzvWghCUIM/sUWCZpv7BoAvCupMExm30ZeCfO7q8D+0qqCDu3zwA6vtPMOedct0p0i8MFwP3hl/xi4FzgZkmjCZqMlgDfA5A0hGA460lm1ijpfODvBMNc7zKz+QmO1TnnXAy/Uc4559JYR8NcfS4m55xzcXmCcM45F5cnCOecc3F5gnDOORdXSnVSS1oNfLyTuw8E1nRjOLsDv+bUl27XC37NXbWXmZXG+yClEsSukDS7vZ78VOXXnPrS7XrBr7k7eROTc865uDxBOOeci8sTxDbTkh1AEvg1p750u17wa+423gfhnHMuLq9BOOeci8sThHPOubjSPkFImijpfUmLJF2R7HgSQdJdklZJeiemrL+kZyQtDJ9Lkhljd5M0TNI/JL0nab6kH4flKXvdkvIkvSbprfCafx6Wp+w1Q7CGvaQ3Jf0tfJ/S1wsgaYmkeZLmSpodlnX7dad1gpCUCdwKnAiMAKZIGpHcqBLiHmBim7IrgOfMbF/gufB9KmkkWM72AGAs8KPw3zaVr7sO+IKZjQJGAxMljSW1rxmCpYzfi3mf6tfbYryZjY65/6HbrzutEwTBanaLzGyxmdUDDwKnJTmmbmdmLwHr2hSfBtwbvr4XmNSjQSWYma0wszfC15sJvkCGksLXbYEt4dvs8GGk8DVLKgO+BEyPKU7Z6+1Et193uieIocCymPdVYVk62MPMVkDwZQoMSnI8CSOpHDgEeJUUv+6wuWUusAp4JlzyN5Wv+b+By4DmmLJUvt4WBjwtaY6kqWFZt193uq9hrzhlPu43hUgqAv4CXGRmm6R4/+Spw8yagNGS+gEzJB2U7JgSRdLJwCozmyPp2GTH08OONLPlkgYBz0hakIiTpHsNogoYFvO+DFiepFh62sqW9cHD51VJjqfbScomSA73m9mjYXHKXzeAmW0AXiDoe0rVaz4SOFXSEoLm4S9I+iOpe72tzGx5+LwKmEHQXN7t153uCeJ1YF9JFeG62WcAjyc5pp7yOHB2+Pps4LEkxtLtFFQV7gTeM7MbYz5K2euWVBrWHJCUDxwHLCBFr9nMfmJmZWZWTvD/7vNmdhYper0tJBVKKm55DZwAvEMCrjvt76SWdBJBO2YmcJeZ/SrJIXU7SQ8AxxJMCbwSuAb4f8BDwHBgKTDZzNp2ZO+2JB0F/BOYx7b26SsJ+iFS8rolHUzQOZlJ8OPvITP7haQBpOg1twibmC41s5NT/Xol7U1Qa4Cgm+BPZvarRFx32icI55xz8aV7E5Nzzrl2eIJwzjkXlycI55xzcXmCcM45F5cnCOecc3F5gnBpQdKAcObLuZI+lfRJzPtXEnTOQyRND19fK+nSRJynnXOPlHRPT53PpaZ0n2rDpQkzW0swwymSrgW2mNkNCT7tlcD/TeQJJGWZWWPbcjObJ6lM0nAzW5rIGFzq8hqES3uStoTPx0p6UdJDkj6QdJ2kM8M1FuZJ+ky4Xamkv0h6PXwcGeeYxcDBZvZWTPEISS9IWizpwphtL5H0Tvi4KCwr1/brd1waJjbCY/xa0ovAjyVNDvd9S9JLMef7K8Edxs7tFK9BOLe9UcABBNOjLwamm9nhChYcugC4CLgJ+J2ZvSxpOPD3cJ9YlQTTH8TaHxgPFAPvS7oNOBg4FziCYPLIV8Mv/vWdxNnPzD4PIGke8EUz+6Rlqo3QbII1Aa6PfPXOxfAE4dz2Xm+ZMlnSh8DTYfk8gi93COY4GhEzM2wfScXhuhMtBgOr2xz7CTOrA+okrQL2AI4CZpjZ1vCcjwJH0/mcYH+Oef0v4B5JDwGPxpSvAoZ0chzn2uUJwrnt1cW8bo5538y2/18ygM+ZWU0Hx6kB8jo4dlN4vPbmH29k+ybgtsfa2vLCzL4v6QiChXPmShod9rnkhXE4t1O8D8K5rnsaOL/ljaTRcbZ5D9gnwrFeAiZJKghn5vwywSSDK4FB4eirXODk9g4g6TNm9qqZ/QxYw7Yp7D/Ljs1czkXmNQjnuu5C4FZJbxP8P/QS8P3YDcxsgaS+cZqeaLPdG+Fw1NfCoulm9iaApF8QzD77EcG03e35L0n7EtRGngNaOsbHA0909eKca+GzuTqXIJIuBjab2fRON+7+c+cCLwJHxRsG61wU3sTkXOLcxvb9Dj1pOHCFJwe3K7wG4ZxzLi6vQTjnnIvLE4Rzzrm4PEE455yLyxOEc865uDxBOOeci+v/Awg828DTmEBSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(t, T_analytical, '-', label = 'analytical')\n",
    "plt.plot(t, T_numerical, '-', label = 'numerical')\n",
    "plt.legend()\n",
    "plt.title('Numerical vs Analytical')\n",
    "plt.xlabel('Time (hours)')\n",
    "plt.ylabel('Temperature (degrees Farenheit)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) The temperature as time goes to infinity is 65 degrees Farenheit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 302,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The corpse was 98.6 degrees Farenheit 51.0 minutes ago from 11:00am.\n",
      "The time of death was 10:09am\n"
     ]
    }
   ],
   "source": [
    "Temp_initial = 98.6\n",
    "T_t = 85\n",
    "\n",
    "t_hours = ((np.log((T_t - T_a)/(Temp_initial - T_a)))/(-K))\n",
    "t_mins = t_hours*60\n",
    "\n",
    "print('The corpse was 98.6 degrees Farenheit', round(t_mins,0), 'minutes ago from 11:00am.')\n",
    "print('The time of death was 10:09am')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Now that we have a working numerical model, we can look at the results if the\n",
    "ambient temperature is not constant i.e. T_a=f(t). We can use the weather to improve our estimate for time of death. Consider the following Temperature for the day in question. \n",
    "\n",
    "    |time| Temp ($^o$F)|\n",
    "    |---|---|\n",
    "    |8am|55|\n",
    "    |9am|58|\n",
    "    |10am|60|\n",
    "    |11am|65|\n",
    "    |noon|66|\n",
    "    |1pm|67|\n",
    "\n",
    "    a. Create a function that returns the current temperature based upon the time (0 hours=11am, 65$^{o}$F) \n",
    "    *Plot the function $T_a$ vs time. Does it look correct? Is there a better way to get $T_a(t)$?\n",
    "\n",
    "    b. Modify the Euler approximation solution to account for changes in temperature at each hour. \n",
    "    Compare the new nonlinear Euler approximation to the linear analytical model. \n",
    "    At what time was the corpse 98.6$^{o}$F? i.e. what was the time of death? \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [],
   "source": [
    "Temp = np.array([55,58,60,65,66,67])\n",
    "Time = np.array([-3,-2,-1,0,1,2])\n",
    "\n",
    "def current_temperature(time_a):\n",
    "    Temp_a = Temp[time_a]\n",
    "    \n",
    "    return Temp_a\n",
    "    print(Temp_a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[55. 58. 60. 65. 66. 67.]\n",
      "[-3. -2. -1.  0.  1.  2.]\n"
     ]
    }
   ],
   "source": [
    "current_temp = np.zeros(6)\n",
    "time = np.zeros(6)\n",
    "\n",
    "for i in range(0,6):\n",
    "    current_temp[i] = current_temperature(i)\n",
    "    time[0] = -3\n",
    "    if 1 <= i <= 3:\n",
    "        time[i-1] = time[0] + [i-1]\n",
    "    elif 4 <= i <= 6:\n",
    "        time[i] = time[0] + [i]\n",
    "    \n",
    "\n",
    "print(current_temp)\n",
    "print(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5xdVX338c83CZH7LQTlIolRwVIETAY7IqJcCnhDW4WqaFGkVCuoVYpavIA+WkVKBbkoBQLYqCCIoiIFLEV5+gyYgSAgoDGQEgEJMcg9yWS+zx97jzmMc9mTmX3OzJzv+/Wa15y9zt5r/9YEfmeftddeS7aJiIj2MaXVAURERHMl8UdEtJkk/oiINpPEHxHRZpL4IyLaTBJ/RESbSeKPiGgzSfwxJiS9Q9JCSU9IelDSjyXt0+q4+ki6T9KBg7x3RBn3E5KeltTbsP1Es2MdS5IOkbS41XHE+JLEH6Mm6SPAV4AvAM8FdgLOBt60HnVNq1I2lmwvsL2p7U2B1wIP9G2XZeOSpCmSav1/uO6/fbRGEn+MiqQtgM8CH7D9XdtP2l5j+we2/6nc50JJ/6fhmNdIWtawfZ+kj0n6BfCkpGmDlG0v6XJJyyXdK+mDDXWcJOlSSRdLelzSnZI6yve+QfFh9IPyKv6E9Wjn8yV9X9IjkpZIel/De1+UtEDSJWX9iyS9QNJnyv3vk7Rfw/5dkj4nqVvSH8o2bdHw/qsk3STpUUm3SHplv2M/K+km4Clge0l/L+nust2LJR1V7jsDuAKY0/ANZoakb0v6ZEOdz/pWIOkhScdLuhN4bLj2x8STxB+j9QpgQ4oEMxpvB14PbGm7p38Z0Av8ALgN2AE4APiwpIMb6jgU+Ha5/5XAmQC23wX8L/DG8ir+lJEEJmkqcBXwP8D2wCHAP0t6dcNufwV8rTz3PcB/AU8CzwP+leIbUKO/BY4o2zK93AdJs4HvAScCWwOfBL4naauGY99ZHr8Z8BDwIMU3lc2B9wFnSfpz2yvKuJY0fINZUbHZfwP8JTCjYvtjAknij9GaATzSkKzX1xm277f99CBlewEzbX/W9mrbS4B/B97WsP+Ntq+yvRb4BrDHKGPqsw+woe0vlef+FTC/37l/Yvv68u9wGUUS/tdy+9vASyRt1LD/fNt3234C+AzFhxzAkcB3bV9nu9f2VcAvgYMajj3P9j3lN6se21favteF64AbyphH499sP1D+7au0PyaQ9N/FaK0AtpE0bZTJ//5hymZRdGs82lA2FfhZw/ZDDa+fAjYcg7j6zj17gHNf17D9u4bXTwPLvW4GxL4Ps00aXje2bSmwcdndMwt4u6TDGt7fgOJKmwGORdKhFN8MXkRxMbcxz/67rI/+f/vh2h8TSBJ/jNb/A54B3kxxpTuQJymSUZ/nDbDPQNPENpbdD9xr+8XrE+Qg9Vd1P3C37ZeOoo7+nt/weifgKdt/kHQ/xRX9cUMc+8e2SNoE+A7wVuDHtnskXQ2o/74NRvrvUUf7o4XS1ROjYvsPwKcp+pXfLGljSRtIeq2kvr70RcDrJG0t6XnAh9fjVDcDj5U3fDeSNFXSbpL2qnj874A563FegBsBJH1Y0obljebdJc1dz/oA3i1pZ0mbAicBl5TlFwGHSTqgbONG5euBkjPARhTfCB4Gesur/9c0vP87YNvyPH0WAW+QtKWkHYChPmSgnvZHCyXxx6jZPg34CEV3w3KKK8RjKW5SQtHffhtwH3AN65LcSM6xFngjsCdwL/AIcB6wxVDHNfgX4JPlSJnjR3juNcDrgL0pumWWA+cAoxnq+Q3gW8BvKW5cf7Q81xLgLcDJFG1cCnyIQf5ftf0IcDzFje8VFN+8rmrY5TaKG91Ly7ZvDVwALKa44f3DMo5B1dT+aCFlIZaI5pLUBZxp+z9aHUu0p1zxR0S0mST+iIg2k66eiIg2kyv+iIg2MyHG8W+zzTaePXt2q8OIiJhQuru7H7E9s3/5hEj8s2fPZuHCha0OIyJiQpG0dKDydPVERLSZJP6IiDaTxB8R0WaS+CMi2kwSf0REm0nij4hoM0n8ERHjVPfSlZx1/WK6l64c03onxDj+iIh20710JUec18Xqnl6mT5vCgqM7mTdrq+EPrCBX/BER41DXkhWs7uml17Cmp5euJSvGrO4k/oiIcahzzgymT5vCVMEG06bQOWfGmNWdrp6IiHFo3qytWHB0J11LVtA5Z8aYdfNAEn9ExLg1b9ZWY5rw+6SrJyKizSTxR0S0mST+iIg2k8QfEdFmkvgjItpMEn9ERJtJ4o+IaDNJ/BERbSaJPyKizSTxR0S0mVoTv6QtJV0m6W5Jd0l6RVl+nKR7JN0p6ZQ6Y4iIiGere66e04Grbb9V0nRgY0n7AW8Cdre9StK2NccQEZNA99KVtUxY1o5qS/ySNgf2Bd4NYHs1sFrS+4Ev2l5Vlj9cVwwRMTnUuShJOxq2q0eFl0o6WNK+kqpOCj0HWA7Ml3SrpPMkbQLsDLxK0k2SbpC01yDnPUbSQkkLly9fXrlBETH51LkoSTsa9Ipf0mzgBOAQ4F6KJL4h8GJJjwJfA/7Dtoeoey5wnO2bJJ0OfLws3wroBPYCLpU0p389ts8FzgXo6OgY7BwR0Qb6FiVZ09M75ouStKOhunpOAc4BjrXd2/iGpO2BdwBHAhcOcvwyYJntm8rtyygS/zLgu2Wiv1lSL7ANxQdLRMSfqHNRknY0aOK3fTiApGlAb7+3H7Z96lAV235I0v2SdrF9D3AA8EvgN8D+wH9L2hmYDjwyijZERBuoa1GSdlTl5u7NFF02w5UN5DhgQTmiZwnwHuBJ4AJJdwCrgSOH6C6KiIgxNlQf/7bAdsBGkl4KqHxrc2DjKpXbXgR0DPDWO0cYZ0REjJGhrvhfDxwF7Aic3VD+OPCpOoOKiIj6DNXHP59iKObhti9tYkwREVGjobp63m77W8B2kj7Y/33bZ9QaWURE1GKorp6+2+fbNCOQiIhojqG6es4uf6c/PyJiEqkyZcOLJP2npNvK7d0lfaL+0CIiog5VpmU+DziZdQ9x3U6GY0ZETFhVEv8mtv+nb6N82GpNfSFFRESdqiT+FZJeABhA0puBh2qNKiKG1L10JWddv5jupStbHUpMQFWmbDgWOB94iaSlwIPA22uNKiIGlbnpY7SGTfy2FwP7S9oCkO1H6w8rIgYz0Nz0SfwxEsMm/nKCtTcDs4FpUjFlj+0v1BpZRAwoc9PHaFXp6rkCeAboBtbWG05EDCdz08doVUn8s2zvVnskEVFZ5qaP0agyqqdL0q61RxIREU0x1CRtt1IM4dwAWCTp18Aqinn5bbvKQiwRETHODNXV89amRREREU0z1CRtv+l7rWIoz8yh9o+IiImhynDOfwA+C6xg3Xw9BtLvHxExAVW5gv8I8Ge2l9cdTERE1K/KqJ5lwO/rDiQiIpqjyhX/YuC/JP2QYlQPkKUXIyImqiqJ/8HyZ/OaY4mIiCaoMknbpwAkPcf2quH2j4iI8a3K0osvl3Q78Otyew9JX61SuaQtJV0m6W5Jd0l6RcN7x0uypCzmHhHRRFVu7p4BvIFiOCe2bwP2q1j/6cDVtl8C7AHcBSDp+cBfAv870oAjImJ0qiT+KbaX9isbdpZOSZsD+1Is4oLt1Q1z+f8bcALlql4REdE8VRL//ZJeDljSVEkfBn5V4bg5wHJgvqRbJZ0naRNJhwK/Lb85DErSMZIWSlq4fHkeIYiIGCtVEv/7KR7i2gn4HdBZlg1nGjAXOMf2y4AngZOAE4FPD3ew7XNtd9jumDlzZoXTRUREFUOO6pE0FXib7betR93LgGW2byq3L6NI/C8AbitX8toRuEXSy21nAfeIiCYY8orf9lrgLetTcZnI75e0S1l0AHCL7W1tz7Y9m+LDYW6SfkRE81R5gOtnkk4Hvk3RXQOA7V9UOPY4YEG5bu8S4D3rFWVERIyZKon/1eXvxoVXTDFiZ0i2FwEdQ7w/u8L5IyJiDFV5cvdVzQgkIiKao9LCKpIOBv4c2LCvzPYX6goqIiLqU2UhlrOBLSm6duZT3OztqjmuiIioSZVx/PvYfgewopyw7S8ohmFGRMQEVCXxP13+fkbS84BngNm1RRQREbWq0sf/Y0lbAqcCiyjm6bmo1qgiIqI2VUb1nFS+/E65CtdGtrMUY0TEBDVoV4+kzzW83h/A9tNJ+hERE9tQffyvb3h9at2BREREc1S5uRsREZPIUH3820r6IKCG139k+4xaI4uIiFoMlfjnAzMHeB0RERPYoIm/fFgrIiImmfTxR0S0mST+iIg2k8QfEdFmhk38ko6VtHn5+uuSbpZ0QP2hRUREHapc8R9j+zFJBwE7AO8HTqk3rIiIqEuVxO/y92uB+ba7Kx4XERHjUJUEfpukq4A3UszUuSnrPgwiImKCqTIt83uAecBi209J2gZ4b71hRUREXYa94re9FphD0bcPsFGV4yIiYnyqMqrnTGA/4J1l0ZPA1+oMKiIi6lPlyn1v239PseQi5Xz806tULmlLSZdJulvSXZJeIenL5fYvJF1Rru4VERFNUiXxr5E0hfKGrqQZQG/F+k8Hrrb9EmAP4C7gWmA327sDvwI+MeKoIyJivVVJ/GcBlwMzJZ0M3Ah8abiDyoe+9gXOB7C92vajtq+x3VPu1gXsuF6RR0TEeqmy5u7FkrqBAynm5j/M9h0V6p4DLAfmS9oD6AY+ZPvJhn2OAi4ZedgREbG+qo7O2QxYafsrwIOSdqpwzDRgLnCO7ZdR3BT+eN+bkk4EeoAFAx0s6RhJCyUtXL58ecUwIyJiOFVG9XwS+AzwybJoQ+CbFepeBiyzfVO5fRnFBwGSjgTeABxhe8CHwWyfa7vDdsfMmVkDJiJirFS54n8r8DqKK3Zs/xbYfLiDbD8E3C9pl7LoAOCXkg4BPgYcavup9Yo6IiLWW5Und1fZtqS+UT0bj6D+44AFkqYDSyieAv458BzgWkkAXbbfN7KwIyJifVVJ/N+VdBawhaT3UEzXcEGVym0vAjr6Fb9oZCFGRMRYqjKq50uSXgusphiL/3nbP649soiIqMWQiV/SVOAq2wcDSfYREZPAkDd3ywnaVvetwBURERNflT7+Jyjm5L+GcmQPgO2P1BZVRETUpkriv678iYiISaDKzd3zmxFIREQ0x7CJX9Kt/OlSi38AFgL/Uk7THBERE0SVrp5rKSZn65um4W0U0zI/DlwIHFpLZBERUYsqiX9v2/s0bN8q6Ubb+0i6va7AIiKiHlXm6tlM0ry+DUlzWTdXT8/Ah0RExHhV5Yr/74FvSNqAostnNfBeSZsAp9QZXEREjL0qo3q6gF3LJRdl+5GGt79VW2QREVGLKvPxz5T0deAi249I2lXSu+sPLSIi6lClj/9C4Abg+eX2r4GP1hVQRETUq0ri39b2NymGcGJ7DbC21qgiIqI2VRL/k5K2pnyIS9JeFGP4IyJiAqoyqud44AfAHEk3ADtQLMcYERETUJVRPQsl7Qf8GcVwzl/aXl17ZBERUYtBE7+kwaZi2EkStq+sKaaIiKjRUFf8h5W/twH2Bq6nuOJ/NcUonyT+iIgJaNDEb/tdAJKuBHa1/dtyewfgjOaEFxERY63KqJ45fUm/9ACwS03xREREzaqM6vmppB9RTM9gimmZf1prVBERUZsqif8DFMM39y23LwYuqy2iiIioVZXhnAa+U/6MiKQtgfOA3Si+LRwF3ANcAswG7gMOt71ypHVHRMT6GbSPX9L1kt4vaft+5dMk7SvpfEnvGab+04Grbb8E2AO4C/g48BPbLwZ+Um5HRESTDHVz9/XABsAVkpZJ+oWkXwNLgPcA59ieP9jBkjan6B46H8D2atuPAm8CLip3uwh48+ibEe2se+lKzrp+Md1L88UxooqhhnM+RTFs8wxJzwG2BZ7uNx//UOYAy4H5kvYAuoEPAc+1/WB5jgclbTuaBkR76166kiPO62J1Ty/Tp01hwdGdzJu1VavDihjXqgznxPYq2/ePIOlD8aEyl+KbwcuAJxlBt46kYyQtlLRw+fLlIzhttJOuJStY3dNLr2FNTy9dS1a0OqSIca9S4l9Py4Bltm8qty+j+CD4naTtAMrfDw90sO1zbXfY7pg5c2aNYcZE1jlnBtOnTWGqYINpU+icM6PVIUWMe1WGc64X2w9Jul/SLrbvAQ4Afln+HAl8sfz9/bpiiMlv3qytWHB0J11LVtA5Z0a6eSIqqJT4Je0IvNj29WV//zTbT1Y49DhggaTprLspPAW4VNJ7gf9l3ZxAEetl3qytkvAjRmDYxC/pKOBYYAvghcAs4GzgwOGOtb0I6BjgrQNGFmZERIyVKn38HwQ6gccAbP+KYoRPRERMQFUS/zONC69ImkoxPXNERExAVRL//5V0ArBhuRLXJcAP6w0rIiLqUiXxn0CxuPrdFA9g/QQ4sc6gIiKiPkPe3C27dS6wfSRwTnNCioiIOg15xW97LbCdpA2aFE9ERNSsyjj+JcDPJH2fYtoFAGxn+cWIiAmoSuJfDlwLbFz+RETEBFZlIZZPNSOQiIhojipP7l5LsXrWs9g+qJaIIiKiVlW6ej7Z8HpD4C3AqnrCiYiIulXp6rmpX9ENkm6oKZ6IiKhZla6ezRs2pwDzgO1qiygiImpVpavnToo+fgE9wL3A39UZVERE1KdK4p9je01jgaTaFnCJiIh6VZmrp38fP8DNYx1IREQ0x6BX7pK2pejL30jSS1k3FfPm5EGuiIgJa6gum9cDRwE7Uqy41edxIA91RURMUIMmftvzgfmSDrd9aRNjioiIGlUZx3+ppIOBP6d4gKuv/At1BhYREfWoMo7/bGBLYF9gPsWTu101xxURETWpMqpnH9vvAFaUE7b9BUW/f0RETECVFlvv+y3peeX27NoiioiIWlV5EOsqSVsCpwKLgLXARbVGFRERtRluzd0pwI9tPwp8R9IPgY1s/75K5ZLuoxj+uRbosd0haU/gaxQ3inuAf7CdB8IiIppkuDV3e4HTG7afrpr0G+xne0/bHeX2KcDJtvcEPl1uxxjpXrqSs65fTPfSla0OJSLGqSpdPddKepPt74/ROU3x9C/AFsADY1Rv2+teupIjzutidU8v06dNYcHRncybtVWrw4qIcaZK4j8W2ELSKuBpiqkbbHvrCscauEaSga/bPhf4MPCfkk6l+Max90AHSjoGOAZgp512qnCq6FqygtU9vfQa1vT00rVkRRJ/RPyJKol/m1HU/0rbD5Tz/lwr6W7grcA/2r5c0uHA+cCB/Q8sPyTOBejo6PiTpR/jT3XOmcH0aVNY09PLBtOm0DlnRqtDiohxSPbwOVXS2yimZ/6CpB2B59ruHtGJpJOAJyjm+dnStiUJ+IPtzYc6tqOjwwsXLhzJ6dpW99KVdC1ZQeecGbnaj2hzkrob7q/+0bDj+CWdCewHvKsseopiVM5wx20iabO+18BBwB0UffqvLnfbH/h1lQZENfNmbcUH9ntRkn5EDKpKV8/etudKuhXA9u8lTa9w3HOBK4qLeqYB37R9taQngNPLxVyeoezHj4iI5qiS+NeU4/kNIGkG0DvcQbaXAHsMUH4jxbq9ERHRAlWmbDgLuByYKelk4EbgS7VGFRERtakyLfPFkrpZN/LmMNt31BtWRETUpeqi6VOBNRTdPVW+JURExDhVZVTPicC3gO0ppmP+pqRP1B1YRETUo8oV/zuBebafApD0eaAb+Jc6A4uIiHpU6bZZyrM/IKYBS+oJJyIi6lbliv8p4E5J/0nRx38QcKOk0wBsf6TG+CIiYoxVSfw/Kn/6ZL3diIgJrMpwzvObEUhERDRHlVE9h0j6uaSHJf1e0kpJI12MJSIixokqXT1nAocDt1NhqoaIiBjfqiT+ZcCichnGiIiY4Kok/hOAH0j6b2BVX6HtM+oKKiIi6lMl8Z9MMV3DlqSrJyJiwquS+Le1nWmUIyImiSpP7v5E0v61RxIREU1RJfH/HXCdpCcynDMiYuKr0tWzTe1RRERE0wx7xW97LXAY8LHy9XbAnnUHFhER9ajy5O6ZwH7Au8qip4Cv1RlURETUp0pXz96250q6FcD27yVNrzmuiIioSZWbu2skTaGYkhlJM8h4/oiICWvQxC+p79vAWcDlwExJJwM3Al9qQmwREVGDobp6bgbm2r5YUjdwICDgMNt3NCW6iIgYc0MlfvW9sH0ncOdIK5d0H/A4sBbosd1Rlh8HHAv0AD+yfcJI646IiPUzVOKfKWnQZRVtn1bxHPvZfqRvQ9J+wJuA3W2vkrRtxXoiImIMDJX4pwKb0nDlP0beD3zR9ioA2w+Pcf0RETGEoRL/g7Y/O8r6DVwjycDXbZ8L7Ay8StLngWeA423/vP+Bko4BjgHYaaedRhlGRET0qdTHPwqvtP1A2Z1zraS7y3NuBXQCewGXSppj240Hlh8S5wJ0dHSYiIgYE0ON4z9gtJXbfqD8/TBwBfByihW9vuvCzRTPBGQ+oIiIJhk08dse1QyckjaRtFnfa+Ag4A7ge8D+ZfnOwHTgkcHqiYiIsVVlyob19VzgCkl95/mm7avL6R4ukHQHsBo4sn83T0RE1Ke2xG97CbDHAOWrgXfWdd6IiBhalbl6IiJiEknij4hoM0n8ERFtJok/IqLNJPFHRLSZJP6IiDaTxB8R0WaS+CMi2kwSf0REm5nUib976UrOun4x3UtXtjqUiIhxo865elqqe+lKjjivi9U9vUyfNoUFR3cyb9ZWrQ4rIqLlJu0Vf9eSFazu6aXXsKanl64lK1odUkTEuDBpE3/nnBlMnzaFqYINpk2hc86MVocUETEuTNqunnmztmLB0Z10LVlB55wZ6eaJiChN2sQPRfJPwo+IeLZJ29UTEREDS+KPiGgzSfwREW0miT8ios0k8UdEtJkk/oiINiPbrY5hWJKWA0vX8/BtgEfGMJyJIG1uD2lzexhNm2fZntm/cEIk/tGQtNB2R6vjaKa0uT2kze2hjjanqycios0k8UdEtJl2SPzntjqAFkib20Pa3B7GvM2Tvo8/IiKerR2u+CMiokESf0REm2mLxC/pc5J+IWmRpGskbd/qmOom6cuS7i7bfYWkLVsdU90kHSbpTkm9kibtkD9Jh0i6R9JiSR9vdTzNIOkCSQ9LuqPVsTSDpOdLul7SXeV/0x8ay/rbIvEDX7a9u+09gR8Cn251QE1wLbCb7d2BXwGfaHE8zXAH8NfAT1sdSF0kTQXOAl4L7Aq8XdKurY2qKS4EDml1EE3UA3zU9p8BncAHxvLfuS0Sv+3HGjY3ASb9HW3b19juKTe7gB1bGU8z2L7L9j2tjqNmLwcW215iezXwbeBNLY6pdrZ/Cvy+1XE0i+0Hbd9Svn4cuAvYYazqn9QrcDWS9Hngb4E/APu1OJxmOwq4pNVBxJjYAbi/YXsZ8BctiiWaQNJs4GXATWNV56RJ/JKuA543wFsn2v6+7ROBEyV9AjgW+ExTA6zBcG0u9zmR4mvjgmbGVpcqbZ7kNEDZpP8G264kbQpcDny4X8/FqEyaxG/7wIq7fhP4EZMg8Q/XZklHAm8ADvAkeWBjBP/Ok9Uy4PkN2zsCD7QolqiRpA0okv4C298dy7rboo9f0osbNg8F7m5VLM0i6RDgY8Chtp9qdTwxZn4OvFjSCyRNB94GXNnimGKMSRJwPnCX7dPGvP5JciE4JEmXA7sAvRTTO7/P9m9bG1W9JC0GngOsKIu6bL+vhSHVTtJfAV8FZgKPAotsH9zaqMaepNcBXwGmAhfY/nyLQ6qdpG8Br6GYovh3wGdsn9/SoGokaR/gZ8DtFHkL4J9tXzUm9bdD4o+IiHXaoqsnIiLWSeKPiGgzSfwREW0miT8ios0k8UdEtJkk/hgzkmaUM6AukvSQpN82bP9Pk2P5cjmr4ZclvXm0E1yVbbte0hOSzuz33ucl3S/pidFFPWwMX5G0b/n6Pknb1Hm+fuc+VdL+zTpf1CvDOaMWkk4CnrB9aovO/xgw0/YqSRcCP7R92QiOn9YwyR2SNqGYL2U3illPj214r5Pi+ZBf2950rNrQL56tgatsd5bb9wEdth8Zw3OIIif0DvDeLODfbR80VueL1skVfzRF39WwpNdIukHSpZJ+JemLko6QdLOk2yW9sNxvpqTLJf28/HnlAHXOlvQzSbeUP3uX5VdSzMJ6k6TPUDyt/eXym8cLy5+rJXWXx7+kPO5CSadJuh74UuO5bD9p+0bgmf5x2O6y/eAA8b1R0k2SbpV0naTnluUnSbqoXBviPkl/LemUsv1Xl4/q9/dW4Op+ZceV7b69oQ1bS/qeinUYuiTt3nDO4xtiu6P8+80u53w/G7gFeH75d7ijrPcfyzYuBWZIGmiepJhgkvijFfYAPgS8FHgXsLPtlwPnAceV+5wO/JvtvYC3lO/19zDwl7bnAn8DnAFg+1Dgadt72j6ZYkqDfyq3f0OxePVxtucBxwNnN9S5M3Cg7Y+OQTtvBDptv4xi+uQTGt57IfB6iimV/wO43vZLgafL8v5eCXT3K3ukbPs5ZTsATgZuLddh+Gfg4gpx7gJcXMa5DbCD7d3KeOY37HdLGUdMcJNmkraYUH7ed4Us6TfANWX57aybMvtAYNei9wGAzSVtVs5N3mcD4ExJewJrKZL2kMrZDvcGvtNQ93MadvmO7bUjb9KAdgQukbQdMB24t+G9H9teI+l2iqkX+q7mbwdmD1DXdsDyfmV9E3d1UyxAA7APxQcltv+rvDexxTBxLrXdVb5eAsyR9FWKyQyvadjvYWDSr17XDpL4oxVWNbzubdjuZd1/k1OAV9h+eoh6/pFi3pY9yv3/pBtmAFOAR8vV2AbyZIU6qvoqcJrtKyW9Bjip4b1VALZ7Ja1pmD218W/Q6Glgw35lfX+3tQ3HDDZtcw/P/obfWNcf22x7paQ9gIOBDwCHU6zn0HfMUP8eMUGkqyfGq2so1k0AoLyq728L4MHyZuS7KK6cB/I4sBn8cTW2eyUdVtarMtHVYQugbzLAI0dZ113Aiyrs91PgCCjup1B0Bz0G3AfMLcvnAi8Y6OBypNAU25cDn+o7prQzxfKWMcEl8cd49UGgo7xJ+d1Gn+sAAADLSURBVEtgoJlFzwaOlNRFkZQGu1r/NvBP5U3WF1IkxvdKug24k4pLF5YjaU4D3i1pWd8Q0fLG7DJg47L8pPKQkyi6lH4GjHb0zY8oZqcczkmUfzfgi6z7wLkc2FrSIuD9FOswD2QH4L/L/S6kXKu5vOH8ImDh+oUf40mGc0ZMEJJuBN5g+9EWnPuvgLm2P9Xsc8fYyxV/xMTxUWCnFp17GvCvLTp3jLFc8UdEtJlc8UdEtJkk/oiINpPEHxHRZpL4IyLaTBJ/RESb+f8PHHYLuOIUnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(time, current_temp, '.')\n",
    "plt.title('Current Temperature')\n",
    "plt.xlabel('Time after 11am (hours)')\n",
    "plt.ylabel('Temperature (degrees Farenheit)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 295,
   "metadata": {},
   "outputs": [],
   "source": [
    "ambient_temp = np.array([55,58,60,65,66,67])\n",
    "hours = np.array([-3,-2,-1,0,1,2])\n",
    "hours_before = -3\n",
    "hours_after = 3\n",
    "N = 50\n",
    "\n",
    "time1 = np.linspace(0,hours_after,N)\n",
    "dt1 = time1[1] - time1[0]\n",
    "temp1 = np.zeros(len(time1))\n",
    "\n",
    "for i in range(0,len(time1)):\n",
    "    if time1[i] <= -3:\n",
    "        temp1[i] = ambient_temp[0]\n",
    "    elif -3 < time1[i] <= -2:\n",
    "        temp1[i] = ambient_temp[1]\n",
    "    elif -2 < time1[i] <= -1:\n",
    "        temp1[i] = ambient_temp[2]\n",
    "    elif -1 < time1[i] <= 0:\n",
    "        temp1[i] = ambient_temp[3]\n",
    "    elif 0 < time1[i] <= 1:\n",
    "        temp1[i] = ambient_temp[4]\n",
    "    elif 1 < time1[i] <= 2:\n",
    "        temp1[i] = ambient_temp[5]\n",
    "    elif 2 < time1[i]:\n",
    "        temp1[i] = ambient_temp[5]\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 291,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dfbwdVX3v8c83T/IkEkLAgEoMAl5UCOHADSIURAViK3oLFrRcRL1YrqBURdFqG+ytIj4VyoPFQIQ2ihhBqUUUKVC5vYmekPAQAoKRSCDAIYA8Sgj53j9mHdgc9jlndnL2ecj5vl+v/doza2bNrHUG9i9rzcxask1ERERdY4a6ABERMbIkcEREREsSOCIioiUJHBER0ZIEjoiIaEkCR0REtCSBIyIiWpLAEcOCpPdJ6pT0hKRVkn4q6S1DXa5uku6W9LZetr2/lPsJSU9LWtew/sRgl3UgSTpU0l1DXY4YXhI4YshJ+gTwj8CXgO2A1wDnAoevx7HG1UkbSLbn2d7C9hbAYcB93eslbViSNEZSW38D2v23j6GRwBFDStIrgC8CH7V9me0nbT9r+99sn1L2+Y6k/9OQ50BJKxvW75b0GUk3A09KGtdL2vaSfiipS9LvJH2s4RizJV0q6WJJj0taKqmjbPsXqmD2b6UV8en1qOerJf1Y0kOSlkv6q4Ztp0uaJ+n75fhLJL1W0t+V/e+WdFDD/gsk/b2kRZL+UOr0iobt+0taKOlRSTdK2q9H3i9KWgg8BWwv6SOSbi/1vkvSB8u+k4DLgWkNLahJki6R9PmGY76oVSLpfkmfkrQUeKy/+sfIk8ARQ21fYBOqH6gNcTTwTmAr22t7pgHrgH8DbgJ2AA4GTpZ0SMMx3gVcUva/AjgbwPYxwO+BPyutiDNaKZikscCVwH8B2wOHAp+T9CcNu70H+FY59x3AfwBPAq8Evk7VAmv0P4H3l7pMKPsgaSrwI+BvgK2BzwM/kjSxIe9flvwvB+4HVlG1lLYE/go4R9IbbK8u5Vre0IJaXbPafwG8HZhUs/4xgiRwxFCbBDzU8GO/vs6yfY/tp3tJ2xuYbPuLttfYXg58GziqYf8bbF9p+zngX4A9NrBM3d4CbGL7K+XcvwHm9jj3NbavLX+H+VQ/4l8v65cAr5e0acP+c23fbvsJ4O+ogiTAscBltn9he53tK4HbgHc05J1j+47Ssltr+wrbv3PlF8D1pcwb4pu27yt/+zr1jxEk/Y8x1FYD20gat4HB455+0nak6pZ5tCFtLPDLhvX7G5afAjYZgHJ1n3tqk3P/omH9gYblp4EuvzACaXcw3LxhubFuK4DNSnfVjsDRko5s2D6e6l/6NMmLpHdRtUxeR/WPyc148d9lffT82/dX/xhBEjhiqP0/4I/Au6n+pd3Mk1Q/Zt1e2WSfZsM8N6bdA/zO9s7rU8hejl/XPcDttt+0Acfo6dUNy68BnrL9B0n3ULUoTuoj7/N1kbQ58APgCOCnttdKugpQz30btHo92lH/GELpqoohZfsPwN9S9au/W9JmksZLOkxS972EJcAsSVtLeiVw8nqc6lfAY+WG+aaSxkp6o6S9a+Z/AJi2HucFuAFA0smSNik36neXNGM9jwfwAUm7SNoCmA18v6RfBBwp6eBSx03LcrMfd4BNqVokDwLrSuvjwIbtDwDblvN0WwL8qaStJO0A9BWkoD31jyGUwBFDzvY3gE9QdZd0Uf0L9USqm7xQ3W+4Cbgb+Dkv/Ei2co7ngD8DpgO/Ax4C5gCv6Ctfgy8Dny9PKn2qxXM/C8wC3kzVrdQFnAdsyKO6/wJ8D7iX6sb/J8u5lgN/DpxGVccVwMfp5f912w8Bn6J6cGA1VcvvyoZdbqJ6UGBFqfvWwIXAXVQPDPyklKNXbap/DCFlIqeIkUXSAuBs2/861GWJ0SktjoiIaEkCR0REtCRdVRER0ZK0OCIioiWj4j2ObbbZxlOnTh3qYkREjCiLFi16yPbknumjInBMnTqVzs7OoS5GRMSIImlFs/R0VUVEREsSOCIioiUJHBER0ZIEjoiIaEkCR0REtKTd8w1vJWl+mZZymaR9y/SYS8rnbklLesl7qKQ7ylSWpzakby3pakl3lu+JzfJHRER7tLvFcSZwle3XU82mtsz2X9iebns68EPgsp6ZylST51BNZ7kb1cQ0u5XNp1LNlrYzcE1Zj1Fs0YpHOOfau1i04pG2pA/GOYby3Bt7/UbrufvbtiHa9h6HpC2BA4APANheA6xp2C7gvcBbm2TfB7irDBGNpEuAw6mmwDycF+YLuAi4DvhMG6oQI8CiFY/w/jkLWLN2HRPGjWHeh2ey144TByx9MM4xlOfe2Os3Ws/d1zkGQjtbHNOoxt2fK2mxpDlltrFu+wMP2L6zSd4dePHUkytLGsB2tlcBlO9tm51c0vGSOiV1dnV1bWhdYphasHw1a9auY53h2bXrWLB89YCmD8Y5hvLcG3v9Ruu5+9u2odoZOMYBM4DzbO9JNd1kY7fS0fQ+AYyapLU0GqPt82132O6YPPklb8zHRmLmtElMGDeGsYLx48Ywc9qkAU0fjHMM5bk39vqN1nP3t21DtW103DJV5QLbU8v6/sCptt8paRzVzGV72V7ZJO++wGzbh5T1zwLY/rKkO4ADba+SNAW4zvaufZWlo6PDGXJk47VoxSMsWL6amdMmvagpPlDpg3GOoTz3xl6/0Xru/rbVIWmR7Y6XpLdzWHVJvwQ+bPsOSbOBzW2fIulQ4LO2/6SXfOOA3wAHUwWYXwPvs71U0leB1bZPL09bbW37032VI4EjIqJ1vQWOdj9VdRIwT9LNVHM9f6mkH0WPbipJ20u6EsD2Wqo5p38GLAMutb207Ho68HZJdwJvL+sRETFIRsVETmlxRES0bqhaHBERsZFJ4IiIiJYkcEREREsSOCIioiUJHBER0ZIEjoiIaEkCR0REtCSBIyIiWpLAERERLUngiIiIliRwRERES/qdAbDM1PdGYHvgaWCp7YGbESQiIkaUXgOHpKnAp4FDgd9Rzea3CbCzpEeBbwH/6tEwSmJERDyvrxbHGcB5wIm21zVukLQ98D7gWOA7bStdREQMO70GDtvvhecnVVrXY/ODtr/WzoJFRMTwVOfm+K9qpkVExCjQa+CQtK2kPYBNJb1J0u7l8xZgszoHl7SVpPmSbpe0rMwljqSTJN0haamkM5rk21XSkobPY5JOLttmS7q3Ydus9at6RESsj77ucbwT+CDwKuDchvTHgS/UPP6ZwFW2j5A0AdhM0kHA4cDutp+RtG3PTLbvoJpqFkljqeYdv7xhl2+mqywiYmj0dY9jLjBX0nttX9rqgSVtCRwAfKAcbw2wRtIJwOm2nynpD/ZzqIOB39pe0WoZIiJi4PXVVXV0WZwi6WM9PzWOPY3qEd65khZLmiNpc2AXYH9JCyVdL2nvfo5zFPC9HmknSrpZ0oWSJvZS/uMldUrq7OrqqlHciIioo6+b490/yNsAk5t8+jMOmAGcZ3tP4Eng1JI+EZgJnAJcWl4yfInSvfUu4AcNyecBO1F1Za0Cvt4sr+3zbXfY7pg8uU5xIyKijr66qs4t33XvZ/S0Elhpe2FZn08VOFYCl5UXB38laR1VcGrWLDgMuNH2Aw3len5Z0reBn6xn+SIiYj30+ziupNdJ+pmkm8r67pI+218+2/cD90jatSQdDNwG/Ah4aznWLsAE4KFeDnM0PbqpJE1pWH0PcGt/ZYmIiIFT5z2OOcBpvPAS4C3AX9Y8/knAPEk3U3UtfQm4EJgm6VbgEuBY25a0vaQruzNK2gx4O3BZj2OeIemWcsyDgL+uWZaIiBgA/Q5yCGxu+7+6b0OUH/ln6xzc9hKgo8mmlwQe2/cBsxrWnwImNdnvmDrnjoiI9qjT4lgt6bWAASS9G7i/raWKiIhhq06L40TgAuD1klZQPcl0dN9ZIiJiY9Vv4LB9F/BWSa8AZPvR9hcrIiKGqzoTOU0A3g1MBcY13Ov4UltLFhERw1KdrqrLgT8Ci4Dn2luciIgY7uoEjh1tv7HtJYmIiBGhzlNVCyTt1vaSRETEiNDXnOOLqR7BHQ8skXQn8Awgqtc5ZgxOESMiYjjpq6vqiEErRUREjBh9DXL42+7lMnrt5L72j4iI0aHO47j/G/gisJoXxqsykPseERGjUJ0WxCeA/2Y7syFFREStp6pWAg+3uyARETEy1Glx3AX8h6SfUD1VBYDts9pWqoiIGLbqBI5V5bNlm8sSEREjQJ1BDr8AIOlltp/pb/+IiNi41Zk6dh9JtwB3lvU9JP1TnYNL2krSfEm3S1omad+SfpKkOyQtlXRGL3nvLjP9LZHU2ZC+taSrJd1ZvifWqmlERAyIOjfHzwL+lOpxXGzfRDVlax1nAlfZfj2wB7BM0kHA4cDutt8AfK2P/AfZnm67cRbBU4FrbO8MXFPWIyJikNQJHGNsr+iR1u8ouZK2BA6gmgQK22vKXB4nAKd3d3vZfrC1InM4cFFZvohqyPeIiBgkdQLHPZL2ASxprKSTgd/UyDcN6ALmSlosaY6kzYFdgP0lLZR0vaS9e8lv4OeSFkk6viF9O9urAMr3ts0ySzpeUqekzq6uvIISETFQ6gSOE6heAnwN8AAws6T1ZxwwAzjP9p7Ak1TdSuOAieU4pwCXqnt2qBfbrwykeBjwUUkH1Djn82yfb7vDdsfkyZNbyRoREX3oM3BIGgscZfso29uUz1G2H6px7JXAStsLy/p8qkCyErjMlV9RDWOyTc/Mtu8r3w9STSa1T9n0gKQppXxTgFa7uiIiYgP0GThsPwf8+foc2Pb9VN1cu5akg4HbgB8BbwWQtAswAXhRIJK0uaSXdy8D7wBuLZuvAI4ty8cCP16f8kVExPqp8wLgLyWdCVxC1d0EgO2ba+Q9CZhX5i1fDhxXjnGhpFuBNcCxti1pe2CO7VnAdsDlpQdrHPBd21eVY55O1b31IeD3wJE1yhEREQNEtvveQfplk2Tbbumew1Dq6OhwZ2dn/ztGRMTzJC3q8ToEUO/N8f3bU6SIiBiJak3MJOkQ4A3AJt1ptr/UrkJFRMTwVWcip3OBrahe5ptLdbN8QZvLFRERw1Sd9zjeYvt9wOoy4OF/B17V3mJFRMRwVSdwPF2+/yjplcAfgaltK1FERAxrde5x/FTSVlSDES6hGqfqor6zRETExqrOU1Wzy+IPyiyAm9rOVLIREaNUr11Vkv6+YfmtALafTtCIiBjd+rrH8c6G5b7mzIiIiFGkzs3xiIiI5/V1j2NbSR8D1LD8PNtntbVkERExLPUVOOYCk5ssR0TEKNZr4Cgv+0VERLxI7nFERERLEjgiIqIlbQ0ckraSNF/S7ZKWSdq3pJ8k6Q5JSyWd0STfqyVdW/IslfTxhm2zJd0raUn5zGpnHSIi4sXqjI57InCx7cck/TOwJ/BZ29fUOP6ZwFW2jyizAG4m6SDgcGB3289I2rZJvrXAJ23fWKaQXSTpatu3le3ftJ13SyIihkCdFsfxJWi8A9gBOAF4SSuhJ0lbUg3FfgGA7TW2Hy35T7f9TEl/sGde26ts31iWHweWlXNHRMQQqxM4uueWPQyYa3tRzXzTgC5grqTFkuZI2hzYBdhf0kJJ10vau6+DSJpK1cpZ2JB8oqSbJV0oaWIv+Y6X1Cmps6urq0ZxIyKijjoB4CZJVwJ/RjVS7ha8EEz6Mg6YAZxne0/gSeDUkj4RmAmcAlwqSc0OUM71Q+Bk24+V5POAnYDpwCrg683y2j7fdoftjsmT8wpKRMRAqRM4jgNmA/vYfopq+tgP1ci3Elhpu7ulMJ8qkKwELnPlV8A6YJuemSWNpwoa82xf1p1u+wHbz9leB3wb2KdGWSIiYoD0GzhsP0fV7XRCSdq0Zr77gXsk7VqSDgZuA34EvBVA0i7ABOChxrylBXIBsMz2N3psm9Kw+h7g1v7KEhERA6fOU1VnA+OpbnT/A1WX07eAPu9NFCcB88oTVcupWi9PAhdKuhVYAxxr25K2B+bYngXsBxwD3CJpSTnW52xfCZwhaTpVd9ndwEfqVjYiIjZcnRkA32x7hqTFALYfLoGgX7aXAB1NNv1lk33vA2aV5RuoBldsdsxj6pw7IiLao849jmcljaHcEJc0ieq+REREjEJ1Asc5VDepJ0s6DbgB+EpbSxUREcNWnTnHL5a0CHgbVffRkbZzQzoiYpSqO1bVy4FHbP8jsErSa9pYpoiIGMbqPFX1eaqnnHYCLqZ6j+O7wFvaW7SIiBiO6rQ4jqB62ulJANv3Alu2s1ARETF81Qkcz9g2LzxVtVl7ixQREcNZncBxmaRzgFdIOg74OXBhe4sVERHDVZ2nqr4i6TCqt7z3AP7B9k/bXrKIiBiW+gwcksYCV9o+BEiwiIiIvruqygCHa8qkTBEREbXGqnqCak6On1OerAKw/Ym2lSoiIoatOoHjF+UTERFR6+b4BYNRkIiIGBnqvDm+mJdOFfsHoBP4su2H21GwiIgYnup0VV1NNbjhd8v6UVTDqj8OfAd4V1tKFhERw1KdFwDfbPsU24vL5zPA/ra/BLy2r4yStpI0X9LtkpZJ2reknyTpDklLJZ3RS95Dyz53STq1IX1rSVdLurN8T2yhvhERsYHqBI6XS9qre0XSDF4Yq2ptP3nPBK6y/XqqlweXSToIOBzY3fYbgK/1zFTeHzkHOAzYDTha0m5l86nANbZ3Bq4p64Nq0YpHOOfau1i04pHa29qdvrGfu79tETF46nRVfQT4F0njqbqs1gAfkrQ50LS1AFDe/TgA+ACA7TVU74ScAJxu+5mS/mCT7PsAd9leXo51CVWwua18H1j2uwi4DvhMjXoMiEUrHuH9cxawZu06Jowbw7wPz2SvHSf2ua3d6Rv7ufv7u0fE4Oq3xWF7ge3dgJnATNu72f5/tp+0/b0+sk4DuoC5khZLmlOCzS7A/pIWSrpe0t5N8u4A3NOwvrKkAWxne1Up2ypg22Ynl3S8pE5JnV1dXf1Vs7YFy1ezZu061hmeXbuOBctX97ut3ekb+7n72xYRg6vfwCFpsqR/Bi6y/ZCk3SR9oMaxxwEzgPNs70n18uCpJX0iVSA6BbhUknqetsnxej7Z1Sfb59vusN0xefLkVrL2aea0SUwYN4axgvHjxjBz2qR+t7U7fWM/d3/bImJwqRoxvY8dpH8H5gGfsb1H6bK60fab+sn3SmCB7allfX+qwDGWqqvqupL+W6qWTFdD3n2B2WWMLCR9FsD2lyXdARxoe5WkKcB1tnftqywdHR3u7Ozss56tWLTiERYsX83MaZNe0l3S27Z2p2/s5+5vW0QMPEmLbHe8JL1G4Pi17b0lLS4tByQtsT29xkl/CXzY9h2SZgObA78Ftrf9t5J2obrB/Ro3FETSOOA3wMHAvcCvgffZXirpq8Bq26eXp622tv3pvsox0IEjImI06C1w1Lk5/qSkrXlhIqe9qd7hqOMkYJ6kCcBy4DiqLqsLJd1KdaP9WNuWtD0wx/Ys22slnQj8jKqFcqHtpeWYp1N1b30I+D1wZM2yRETEAKjT4uigeqz2DcBNVDepj7C9pP3FGxhpcUREtG69Wxy2O8u7F/+N6qb1beXR2oiIGIV6DRySehtK5DWSsH1Fm8oUERHDWF8tju57B9sAbwaupWpx/AlwPZDAERExCvUaOGwfAyDpCmA32/eW9R2AswaneBERMdzUGatqWnfQKO4D+nxvIiIiNl51Hsf9z/IS4PeoHsk9CvjPtpYqIiKGrTqB46PAEVQDFgJcDMxvW4kiImJYq/M4roEflE9ERIxyvd7jkHStpBPKG92N6eMkHSDpAknHtb+IERExnPTV4ngn8GHg8vIk1cPApsDLqMaXOsd2XseOiBhl+noc9ymqx27PkvQyqnkvnrb90GAVLiIihp86N8cps/Xd0++OERGx0avzHkdERMTzEjgiIqIltQKHpFeVEXKR9LIyd3hERIxCdeYc/yDVgIZzStKOwI/bWaiIiBi+6rQ4PgbMBB4DsP0bqies+iVpK0nzJd0uaZmkfSXNlnSvpCXlM6tJvl0bti+R9Jikk8u2fvNHRET71Hmq6o+210gCQNJYquHV6zgTuMr2EWX62M2AQ4Bv2v5ab5ls3wFMbzjfvcDlDbv0mT8iItqnTovj/0r6NLBJuc/xfeAn/WWStCXV+FYXANheY/vR9SjjwcBvba9Yj7wRETHA6gSOTwOPA7cDH6d6a/xvauSbBnQBcyUtljSn4ab6iZJulnShpIn9HOcoqpF5G/WbX9LxkjoldXZ1ddUobkRE1KFqDMNeNlbdRBfaPrblA0sdwAJgP9sLJZ1JdZ/kbOAhqiHa/x6YYvuDvRxjAtX8H2+w/UBJ265u/m4dHR3u7MzoKBERrZC0yHZHz/Q+Wxy2nwOmSBq/HudcCay0vbCszwdm2H7A9nO21wHfBvbp4xiHATd2B41SplbyR0TEAKtzc3w58EtJPwae7E603ef0sbbvl3SPpF3Lze6DgdskTbG9quz2HuDWPg5zND26qVrMHxERA6xO4OgCrqZ6ImqzFo9/EjCvdDktB46jGjRxOlVX093ARwDK8O1zbM8q65sBb+/e3uCMZvkjImJw9HmPY2ORexwREa3r7R5Hvy0OSVdT/ev+RWy/Y4DKFhERI0idrqrPNyxvAvw58Ex7ihMREcNdnTnHF/ZIul7S9W0qT0REDHN1uqq2bFgdA+wFTGlbiSIiYlir01W1lOoeh4C1wO+A/9XOQkVExPBVJ3BMs/1sY4KkWlPORkTExqfOWFU973EA/GqgCxIRESNDry0HSdtS3cvYVNKbeGEo9S1p/UXAiIjYSPTV5fRO4IPAq4BzG9IfB77QzkJFRMTw1WvgsD2Xakj099q+dBDLFBERw1id9zgulXQI8AaqFwC707/UzoJFRMTwVOc9jnOBrahm85tL9eb4gjaXKyIihqk6T1W9xfb7gNW2vwD8d6r7HhERMQrVCRx/7P6W9MqyPrVtJYqIiGGtzot8V0raCvgasAR4DrioraWKiIhhq8/AIWkM8FPbjwI/kPQTYFPbD9c5eAk4c4A3Ug1b8kHgEKohS7rKbp+zfWWTvHdTPfr7HLC2e0x4SVsD36dq9dwNvNf2I3XKExERG66/OcfXAWc2rD9dN2gUZwJX2X49sAewrKR/0/b08nlJ0GhwUNmncSKRU4FrbO8MXFPWIyJikNS5x3G1pMNbPXAZVfcA4AIA22tKy2VDHc4LXWUXAe8egGNGRERNdQLHicDlkp6W9LCkRyTVaXVMo+qOmitpsaQ5kjbvPqakmyVdKGliL/kN/FzSIknHN6RvZ3sVQPnetllmScdL6pTU2dXV1WyXiIhYD3UCxzbAeGALYHJZn1wj3zhgBnCe7T2BJ6m6lc4DdgKmA6uAr/eSfz/bM4DDgI9KOqDGOZ9n+3zbHbY7Jk+uU9yIiKij38Bh+zngSOAzZXkK1Y9+f1YCKxtmEJwPzLD9gO3nyv2TbwP79HLe+8r3g8DlDfs9IGkKQPl+sEZZIiJigPQbOCSdDRwEHFOSngK+1V8+2/cD90jatSQdDNzW/aNfvAe4tck5N5f08u5l4B0N+10BHFuWjwV+3F9ZIiJi4NR5j+PNtmdIWgxg+2FJE2oe/yRgXtl/OXAccJak6VT3MO4GPgIgaXtgju1ZwHZU91W6y/hd21eVY54OXCrpQ8DvqVpDERExSOoEjmfL+xwGkDQJWFfn4LaXAB09ko/pZd/7gFlleTnV47vN9ltN1XqJiIghUOfm+DnAD4HJkk4DbgC+0tZSRUTEsFVnWPWLJS0C3laSjrT9kvsSERExOtTpqgIYCzxL1V1Vp5USEREbqTpPVf0N8D1ge6rh1L8r6bPtLlhERAxPdVocfwnsZfspAEn/ACwCvtzOgkVExPBUp9tpBS8OMOOoHq2NiIhRqE6L4ylgqaSfUd3jeAdwg6RvANj+RBvLFxERw0ydwPHv5dMt841HRIxidR7HvWAwChIRESNDnaeqDpX0a0kPtjisekREbITqdFWdDbwXuIWaQ41ERMTGq07gWAksKcOgR0TEKFcncHwa+DdJ1wHPdCfaPqtdhYqIiOGrTuA4jWq4ka1IV1VExKhXJ3Bsa3uvtpckIiJGhDpvjl8j6a1tL0lERIwIdQLH/wJ+IemJVh/HlbSVpPmSbpe0TNK+kmZLulfSkvKZ1STfqyVdW/IslfTxhm395o+IiPap01W1zQYc/0zgKttHlOljNwMOAb5p+2t95FsLfNL2jWXu8UWSrrZ9W9neX/6IiGiTflsctp+jmtf7M2V5CjC9v3yStgQOAC4ox1lj+9E6hbK9yvaNZflxYBmwQ528ERHRXnXeHD8bOIgX5gp/CvhWjWNPA7qAuZIWS5ojafOy7URJN0u6UNLEfs4/FdgTWNiQ3G9+ScdL6pTU2dXVVaO4ERFRR517HG+2/RHgjwC2HwYm1Mg3DpgBnGd7T+BJ4FTgPGAnqlbLKuDrvR1A0hZU852fbPuxklwrv+3zbXfY7pg8eXKN4kZERB11AsezksZQDamOpEnUe59jJbDSdndLYT4ww/YDtp8rb6J/G9inWWZJ46mCxjzbl3Wn180fERHt0WvgkNR94/wcqh/wyZJOA24AvtLfgW3fD9wjadeSdDBwm6QpDbu9B7i1yblFdW9kme1v9NjWb/6IiGifvp6q+hVVC+FiSYuAtwECjrRd98f6JGBeeaJqOXAccJak6VQtmLuBjwBI2h6YY3sWsB/VPZVbJC0px/qc7SuBM5rlj4iIwSHbzTdIi8u9iRGvo6PDnZ2dQ12MiIgRRdIi2x090/tqcUyW1Ou0sD27kCIiYnToK3CMBbag6p6KiIgA+g4cq2x/cdBKEhERI0Jfj+OmpRERES/RV+A4eNBKERERI0avgaO8IR4REfEidd4cj4iIeF4CR0REtCSBIyIiWpLAERERLUngiIiIliRwRERESxI4IiKiJQkcERHRkgSOiIhoSQJHRES0pK2BQ9JWkuZLul3SMkn7Spot6V5JS8pnVi95D5V0h6S7JJ3akL61pKsl3Vm+J7azDhER8WLtbnGcCVxl+/XAHiYsIykAAAnQSURBVMCykv5N29PL58qemSSNpZrr/DBgN+BoSbuVzacC19jeGbimrLfFohWPcM61d7FoxSPtOkVExIjT13wcG0TSlsABwAcAbK8B1ki1RmvfB7jL9vJyrEuAw4HbyveBZb+LgOuAzwxcySuLVjzC++csYM3adUwYN4Z5H57JXjumcRMR0c4WxzSgC5grabGkOZI2L9tOlHSzpAt76WraAbinYX1lSQPYzvYqgPK9bbOTSzpeUqekzq6urpYLv2D5atasXcc6w7Nr17Fg+eqWjxERsTFqZ+AYB8wAzrO9J/AkVbfSecBOwHRgFfD1JnmbNUvcysltn2+7w3bH5MmTWyo4wMxpk5gwbgxjBePHjWHmtEktHyMiYmPUtq4qqlbCStsLy/p84FTbD3TvIOnbwE96yfvqhvVXAfeV5QckTbG9StIU4MGBLzrsteNE5n14JguWr2bmtEnppoqIKNrW4rB9P3CPpF1L0sHAbeXHvtt7gFubZP81sLOk10qaABwFXFG2XQEcW5aPBX484IUv9tpxIh896HUJGhERDdrZ4gA4CZhXfvyXA8cBZ0maTtX1dDfwEQBJ2wNzbM+yvVbSicDPgLHAhbaXlmOeDlwq6UPA74Ej21yHiIhoILulWwcjUkdHhzs7O4e6GBERI4qkRbY7eqbnzfGIiGhJAkdERLQkgSMiIlqSwBERES0ZFTfHJXUBK9Yz+zbAQwNYnKGUugw/G0s9IHUZrjakLjvafskb1KMicGwISZ3NnioYiVKX4WdjqQekLsNVO+qSrqqIiGhJAkdERLQkgaN/5w91AQZQ6jL8bCz1gNRluBrwuuQeR0REtCQtjoiIaEkCR0REtCSBo5B0qKQ7JN0l6SXzmKtyVtl+s6QZQ1HOOmrU5UBJf5C0pHz+dijK2Z8yQ+SDkpoNvT9irkmNeoyI6wEg6dWSrpW0TNJSSR9vss9IuS516jLsr42kTST9StJNpR6nNdlnYK+J7VH/oRq6/bdU091OAG4Cduuxzyzgp1SzE84EFg51uTegLgcCPxnqstaoywFUs0je2sv2kXJN+qvHiLgepaxTgBll+eXAb0bw/yt16jLsr035O29RlscDC4GZ7bwmaXFU9gHusr3c9hrgEuDwHvscDlzsygJgqx6TUg0XdeoyItj+T+DhPnYZEdekRj1GDNurbN9Ylh8HlgE79NhtpFyXOnUZ9srf+YmyOr58ej71NKDXJIGjsgNwT8P6Sl76H1CdfYaDuuXctzRtfyrpDYNTtAE3Uq5JHSPuekiaCuxJ9S/cRiPuuvRRFxgB10bSWElLqKbSvtovTNndbUCvSbtnABwp1CStZ8Sus89wUKecN1KNQfOEpFnAj4Cd216ygTdSrkl/Rtz1kLQF8EPgZNuP9dzcJMuwvS791GVEXBvbzwHTJW0FXC7pjbYb76kN6DVJi6OyEnh1w/qrgPvWY5/hoN9y2n6su2lr+0pgvKRtBq+IA2akXJM+jbTrIWk81Q/tPNuXNdllxFyX/uoy0q6N7UeB64BDe2wa0GuSwFH5NbCzpNeW+dGPAq7osc8VwP8sTyfMBP5ge9VgF7SGfusi6ZWSVJb3ofrvYPWgl3TDjZRr0qeRdD1KOS8Altn+Ri+7jYjrUqcuI+HaSJpcWhpI2hR4G3B7j90G9JqkqwqwvVbSicDPqJ5KutD2Ukl/VbZ/C7iS6smEu4CngOOGqrx9qVmXI4ATJK0FngaOcnn0YjiR9D2qp1q2kbQS+DuqG38j6prUqMeIuB7FfsAxwC2lTx3gc8BrYGRdF+rVZSRcmynARZLGUgW2S23/pJ2/XxlyJCIiWpKuqoiIaEkCR0REtCSBIyIiWpLAERERLUngiIiIliRwxLAhaVLDKKT3S7q3Yf2/BrksXy0jjX5V0rsl7baBx5tURmJ9QtLZPbb9g6R7JD3RW/6BIOkfJR1Qlu8ezBfZJH1N0lsH63zRXnkcN4YlSbOBJ2x/bYjO/xgw2fYzkr5DNULq/Bbyj7O9tmF9c6qxkN4IvNH2iQ3bZgIrgDttbzFQdehRnq2BK23PLOt3Ax22HxrAc4jqN2Vdk207At+2/Y6BOl8MnbQ4YkTo/te4qvkRrpd0qaTfSDpd0vtVzUdwi6Sdyn6TJf1Q0q/LZ78mx5wq6ZeSbiyfN5f0K4DNgYWS/g54F/DV0vLZqXyukrSo5H99yfcdSd+QdC3wlcZz2X7S9g3AH3uWw/aCZm/xSvozSQslLZb0C0nblfTZki6S9PPScvgfks4o9b9K1TAaPR0BXNUj7aRS71sa6rC1pB+pmrNhgaTdG875qYay3Vr+flNVzWdxLtW4Tq8uf4dby3H/utRxBTBJ0iublC1GmASOGIn2AD4OvInqzd9dbO8DzAFOKvucCXzT9t7An5dtPT0IvN32DOAvgLMAbL8LeNr2dNunUQ3XcEpZ/y1wPnCS7b2ATwHnNhxzF+Bttj85APW8gWpehT2phsf/dMO2nYB3Ug2X/a/AtbbfRPV28zubHGs/YFGPtIdK3c8r9QA4DVhse3eqt6gvrlHOXamG7N4T2AbYwfYbS3nmNux3YylHjHAZciRGol93/wtd0m+Bn5f0W4CDyvLbgN2q3hMAtpT08jLvQrfxwNmSpgPPUf3o90nVSKpvBn7QcOyXNezygzJS6UB4FfB9VfMmTAB+17Dtp7aflXQL1dAy3a2JW4CpTY41BejqkdY9qN8i4H+U5bdQBVps/0e5N/OKfsq5oszxALAcmCbpn4B/54VrA1Wg3r6fY8UIkMARI9EzDcvrGtbX8cJ/02OAfW0/3cdx/hp4gKoFM4Ym3UhNjAEetT29l+1P1jhGXf8EfMP2FZIOBGY3bHsGwPY6Sc82jJ/U+Ddo9DSwSY+07r/bcw15eht+ey0v7qFoPNbzdbb9iKQ9gEOAjwLvBT7YkKev6xEjRLqqYmP1c6DxBnSzH/pXAKvKzdxjqP7l3szjVFOLUuZr+J2kI8txVX4o2+EVwL1l+dgNPNYy4HU19vtP4P1Q3U+i6s56DLibavpbVM1X/dpmmcuTWmNs/xD4QneeYheg6bzrMbIkcMTG6mNAR7nJexvwV032ORc4VtICqh+13loLlwCnlJvUO1H9sH5I0k3AUmpOzVueZPoG8AFJK7sf8S03tlcCm5X02SXLbKousV8CG/r0079TjdDbn9mUvxtwOi8ErB8CW6saRfYEqvm5m9kBuK7s9x3gs/D8vBevAzrXr/gxnORx3IhRQtINwJ+WyX4G+9zvAWbY/sJgnzsGXlocEaPHJylzTQyBccDXh+jcMcDS4oiIiJakxRERES1J4IiIiJYkcEREREsSOCIioiUJHBER0ZL/D49OsUinQvbxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(time1, temp1, '.')\n",
    "plt.title('Current Temperature')\n",
    "plt.xlabel('Time after 11am (hours)')\n",
    "plt.ylabel('Temperature (degrees Farenheit)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 303,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_numerical = np.zeros(len(time1))\n",
    "temp_initial = 85\n",
    "temp_numerical[0] = temp_initial\n",
    "\n",
    "for i in range(1,len(time1)):\n",
    "    k1 = -dTdt/(74 - temp1[i])\n",
    "    temp_numerical[i] = temp_numerical[i-1] + (-k1*(temp_numerical[i-1] - temp1[i-1])*dt1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 304,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_numerical1 = np.zeros(len(time1))\n",
    "temp_numerical1[0] = temp_initial\n",
    "\n",
    "for i in range(1,len(time1)):\n",
    "    temp_numerical1[i] = temp_numerical1[i-1] + (-k*(temp_numerical1[i-1] - Temp_ambient)*dt1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 305,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3gUVffA8e9JIySG0JuUACJCSAiQhN6RIoj0rhQRUEDE/voqYHktiAVFQKQrXUGKqIBSBKT3XgNShIBSQk05vz9mkt8SUjYkYQPcz/Psw065M2dCsmfnzsy5oqoYhmEYRmJurg7AMAzDyJpMgjAMwzCSZBKEYRiGkSSTIAzDMIwkmQRhGIZhJMkkCMMwDCNJJkEY6SYiESLSMJlltURk352OyTCM9DMJ4j4mIstF5F8RyZZZ+1DVP1S1TEZsK5VE1EVEouzXVRGJc5iOyoj9u4qINBGRg3d4nzNE5Lr98/tHRH4RkdJOtPMWERWRIg7z7nj8RsYwCeI+JSIBQC1AgRYuDSYDqOpUVX1AVR8AmgIn46fteVmSiLiJSKb+HYqIx202fdf+2RUDLgHfZFxUzktH/EY6mQRx/3oKWAtMAro5LhCRSSIySkR+tr9BrhaRgiLyuX3GsVdEKibaXpiI7LaXTxQRb3tbdUXkuMO2C4vIDyISKSJHROR5h2VDRWSWiEwRkUsisktEQu1l32J9UC2wY3o1rQcsIkVFZJ6InBWRwyLS12HZhyIyVURm2tvfKiIlRGSIvX6EiNRzWH+tiLwrIptE5IJ9TP4Oy2uJyDoROS8im0WkRqK274jIOuAKUFhE+tg/10siclBEetrr5gHmAiUdzojy2N/w33TY5k3f0kXkbxF5WUR2ARdTO/6UqGoUMAMob2+nhn1sF0TkpIh85vAhvtL+d58da5tk4ncXkbfsOM7aP/uc9vYfEZEYEXlGRP4CFjnM6yEix+3fn1ccjreGiGwRkYv2sX/gzLEZqVBV87oPX8BB4DmgMhANFHBYNgk4ay/zBn4HjmAlFXfgPWCZw/oRwE6gKJAbWA28Zy+rCxy337sBm4DBgBdQEjgMNLaXDwWuAY/Z+/kAWJtoPw2dOLaEfTrMcwd2AK/Z+34YOAbUsZd/iPVhXQ/wAGbax/yyPT0A2OOwvbXAUeAR4AFgATDOXhYAnAMa2sf8GBAJ5HJoexgoA3ja228BlADEbncVCLTXbwIcTHQ8M4A3HaZvWgf4G9gAFAayp3b8SfwME7YP5AC+B5bY0+FAmL3NUli/S33tZd5YZ6VFkovNnvc68IcdnzfW79xEe9kj9jbGAT52/PHzvrLXDwNuACXtNluAdvZ7P6CKq//G7oWXywMwLxf8p0NNrKSQ157eCwxyWD4J+MZhOvGHYxBw3mE6Iv4Dwp5+DDhkv6/L/yeIKsCxRLH8x+GDYSiw1GFZOeBqov3cboKoAxxINO9tYLT9/kNggcOydlgf8mJP57M/oLLb02uBoQ7rVwIu2++HOP787HkrgA4Obd9I5Rh+AfrY7283QXR29viT2P8MrCR1HjiFdRZQPJl1Xwem2++dTRBHgBoO0yWwErTw/8mgsMPy+Hl5HeZtB1ra79cD/wXyuPrv6156mb69+1M3YLGqnrWnp9nzPnNY57TD+6tJTCfu1//L4f1RrG+GiRXH6k457zDPHeubZLy/Hd5fAbxFxENVY5I5FmcVBwKS2PdSh+nExxip9qePPQ3g6/A+8TH72N1MxYFOItLOYbknN/9MHNsiIi2AN4GHsM46fLj553I7HPfhzPEn9j9VfS/xTBEpB3yClRSzY50BrXY2KBERrLPNRSLiWC3UDchjv49T1ZOJmsY6/M6C9fsR/3vYDesLxn67q22wqv7qbExG0kyCuM+ISHagPeAuIvEfxtmAnCJSQVW33eamizq8LwYk/uMG6wPriKqmejdMMtJTevgvYK+qBqVjG4klPuYrqnrB7jcfp6oDUmibcCwi4gvMBtoCP6tqjIj8gvVt+qZ1HVzGSiLxCqa0DzL2+L8BlmN16USJyOtY3WLJxXrTPFVVETkBtFbVTYlXFpG8yWwnWaq6B+ggIu5AR2COiORS1Rtp2Y5xM3OR+v7TEojF6r4JsV9lsb6tPpWO7fYTkSIikht4A6sPP7H1wEUReU1EstsXKsuLSJiT+ziNdd3idqwCEJEXxLoV00NEgkWk0m1uD6C7iDwsIg9gfXuNP+bJQDsRaWAfY3b7fVIf4mB9C/cEzgBx9tlEXYflp4H89n7ibQWai0hOEXkQqxswJRl5/H7ABTs5BALPxC9Q1evABW7+f0oq/jHAhyJS1I4rv4g8fhuxYLd/SkTyqGqsvX8F4m53e4bFJIj7TzesPv9jqvp3/AsYCXSR27+lcBqwGOvi62GsC9k3sf94H8dKSkewLoSPA/wTr5uMD4A37TuDXk5LcKoajXVtpDpWd1AkMJpbu8rS4ltgOnAC68PoJXtfh4E2WH38Z+39DSSZvze72+RlrAvd57CS+CKHVbYB84Gj9rHnBiZgXRw+Biy040hWBh//IKCXWM+XfMWtXwYGA7PtWFskE/8wrO6t30XkErAGq8vqdjXHunPqEtbvSfsM6Ja878VfgDMMIw1EZC0wUlW/c3UshpFZzBmEYRiGkSSTIAzDMIwkmS4mwzAMI0nmDMIwDMNI0j31HETevHk1ICDA1WEYhmHcNTZt2nRWVfMlteyeShABAQFs3LjR1WEYhmHcNUTkaHLLTBeTYRiGkSSTIAzDMIwkmQRhGIZhJOmeugZhGI6io6M5fvw4165dc3UohuFy3t7eFClSBE9PT6fbmARh3LOOHz+On58fAQEBWBWmDeP+pKqcO3eO48ePU6JECafbZfZYuIPEGjZyp4hMt6tIDhWRE2IN6bhVRB5Lpm0TEdkn1vCLr2dmnMa96dq1a+TJk8ckB+O+JyLkyZMnzWfTmZYg7BLEzwOhqloea3CSjvbiz1Q1xH4tSqKtO1aVyKZYZak72YOUGEaamORgGJbb+VvI7IvUHkB2u4S0D0kPIpOUcKwhCg/bA37MAJ7IjABjo2+wZvJbHNz8e2Zs3jAM466VaQlCVU8Aw7Hq1Z/CGmBksb24v4hsF5EJIpIrieYPcvNwicftebcQkd4islFENkZGRqY5zqjLUZQ6MhW3hS9w5erV1BsYRhr8/fffdOzYkVKlSlGuXDkee+wx9u/fn2Hb//HHH9m9e/dtt4+IiGDatGkprvPZZ5/h7e3NhQsX0rz9unXrJvnw6vz58/nwww/TvD2A8+fPM2rUqDQvy2o+//xzrly54uowUpSZXUy5sL71l8Aai9dXRLpiDVJSCmvQmFNYY9ve0jyJeUlWFVTVsaoaqqqh+fIl+bR4ivxz5ubfuh9QMu4oKycPTnN7w0iOqtKqVSvq1q3LoUOH2L17N++//z6nT59OvbGT7kSCmD59OmFhYcydO/e295NYixYteP3127u0eLckCFUlLi75Qe1uJ0HExNzhMZBUNVNeQDtgvMP0U8CoROsEADuTaFsN+NVh+j/Af1LbZ+XKlfV27R3xhF4dnEdX/Ln2trdhZC27d+926f5/++03rVWrVpLL4uLi9OWXX9bAwEAtX768zpgxQ1VVly1bpnXq1NE2bdpomTJltHPnzhoXF6eqqq+99pqWLVtWg4KC9KWXXtLVq1drrly5NCAgQCtUqKAHDx7UsWPHamhoqAYHB2vr1q318uXLqqrarVs3HTBggFarVk1LlCihs2fPVlXVKlWqaI4cObRChQr66aef3hLnwYMHtWzZsrp8+XJt1KhRwvyJEyfqE088oc2bN9eAgAD98ssv9ZNPPtGQkBCtUqWKnjt3TlVV69SpowMHDtRq1appYGCgrlu3LqF9v379VFX1zJkz2rp1aw0NDdXQ0FBdtWqVqqoOGTJEe/TooXXq1NESJUroiBEjVFW1Q4cO6u3trRUqVNCXX375pniTWjZs2DANDQ3VoKAgHTx4sKqqHjlyRMuUKaNPP/20BgYGaufOnXXJkiVavXp1feihhxLiHDJkiHbt2lXr1aunDz30kI4dOzZhX8lt95FHHtFnn31WQ0JCNCIiQvv27auVK1fWcuXKJaw3YsQI9fT01PLly2vdunVVVdXX1zdh27Nnz9Zu3bol/N8NGjRI69atqy+++KJGRUVpjx49NDQ0VENCQvTHH39M8ncsKUn9TQAbNZnP1My8zfUYUFVEfICrQANgo4gUUtVT9jqtgJ1JtN0AlBaREljDOXYEOmdirJR48iuivwjF65eXOF12CQX8s2fm7ow77O0Fu9h98mKGbrNc4RwMeTww2eU7d+6kcuXKSS6bM2cOW7duZdu2bZw9e5awsDBq164NwJYtW9i1axeFCxemRo0arF69mnLlyjF37lz27t2LiHD+/Hly5sxJixYtaN68OW3btgUgZ86cPPOMNUT0m2++yfjx4xkwwBqu+tSpU6xatYq9e/fSokUL2rZty4cffsjw4cNZuHBhknFOnz6dTp06UatWLfbt28eZM2fInz9/wvFt2bKFa9eu8dBDD/HRRx+xZcsWBg0axJQpU3jhhRcAuHz5MmvWrGHlypX07NmTnTtv/pMfOHAggwYNombNmhw7dozGjRuzZ88eAPbu3cuyZcu4dOkSZcqU4dlnn+XDDz9k586dbN269ZZ4Ey9bvHgxBw4cYP369agqLVq0YOXKlRQrVoyDBw8ye/Zsxo4dS1hYGNOmTWPVqlXMnz+f999/nx9//BGA7du3s3btWi5fvkzFihVp1qwZO3fuTHa7+/btY+LEiQlnMv/73//InTs3sbGxNGjQgO3bt/P888/z6aefsmzZMvLmzZvs71C8/fv3s3TpUtzd3XnjjTeoX78+EyZM4Pz584SHh9OwYUN8fX1T3U5aZeY1iHXA98BmYIe9r7HAMBHZISLbgXpY49siIoVFZJHdNgboD/wK7AFmqequzIoVwCvXg1yp/SZV2cGcSZ8SF2fGyTAyz6pVq+jUqRPu7u4UKFCAOnXqsGHDBgDCw8MpUqQIbm5uhISEEBERQY4cOfD29qZXr17MmTMHHx+fJLe7c+dOatWqRVBQEFOnTmXXrv//s2nZsiVubm6UK1fO6W6uGTNm0LFjR9zc3GjdujWzZ89OWFavXj38/PzIly8f/v7+PP744wAEBQURERGRsF6nTp0AqF27NhcvXuT8+fM37WPp0qX079+fkJAQWrRowcWLF7l06RIAzZo1I1u2bOTNm5f8+fOnuXtu8eLFLF68mIoVK1KpUiX27t3LgQMHAChRogRBQUG4ubkRGBhIgwYNEJFb4n/iiSfInj07efPmpV69eqxfvz7F7RYvXpyqVasmtJ81axaVKlWiYsWK7Nq167a6BNu1a4e7u3vCMX344YeEhIRQt25drl27xrFjx9K8TWdk6oNyqjoEGJJo9pPJrHsSa1D1+OlF3Dxwe6bLV+dZIrfPosM/o5i6rAlPNkj6259x90npm35mCQwM5Pvvv09ymaYwUFe2bNkS3ru7uxMTE4OHhwfr16/nt99+Y8aMGYwcOZLff7/1zrvu3bvz448/UqFCBSZNmsTy5cuT3G5K+4+3fft2Dhw4wKOPPgrAjRs3KFmyJP369btle25ubgnTbm5uN/WVJ769MvF0XFwcf/75J9mz33rWntTPIi1Ulf/85z/06dPnpvkRERHpij+l7Tp+kz9y5AjDhw9nw4YN5MqVi+7duyf7LILjfhKv47hNVeWHH36gTJkyKR57RjC1mBy5uZG34yj85SoPrBia4V0Sxv2lfv36XL9+nW+++SZh3oYNG1ixYgW1a9dm5syZxMbGEhkZycqVKwkPD092W1FRUVy4cIHHHnuMzz//PKELxc/PL+HbNsClS5coVKgQ0dHRTJ06NdUYE7d3NH36dIYOHUpERAQRERGcPHmSEydOcPRostWhkzRz5kzAOmvy9/fH39//puWNGjVi5MiRCdNJdR05G3PiZY0bN2bChAlERUUBcOLECc6cOZOm+OfNm8e1a9c4d+4cy5cvJywszOntXrx4EV9fX/z9/Tl9+jQ///xzsrEWKFCAPXv2EBcXl+INAY0bN+bLL79MSPJbtmxJ0/GkhUkQiUiBQG5UfZ5WbiuZ9N0krkXHujok4y4lIsydO5clS5ZQqlQpAgMDGTp0KIULF6ZVq1YEBwdToUIF6tevz7BhwyhYsGCy27p06RLNmzcnODiYOnXq8NlnnwHQsWNHPv74YypWrMihQ4d49913qVKlCo8++iiPPPJIqjEGBwfj4eFBhQoVErYZb8aMGbRq1eqmea1atWLGjBlp+jnkypWL6tWr07dvX8aPH3/L8i+++IKNGzcSHBxMuXLlGDNmTIrby5MnDzVq1KB8+fK88sorKS5r1KgRnTt3plq1agQFBdG2bdtkk0tywsPDadasGVWrVuWtt96icOHCTm+3QoUKVKxYkcDAQHr27EmNGjUSlvXu3ZumTZtSr149wLp+0rx5c+rXr0+hQoWSjeett94iOjqa4OBgypcvz1tvvZWm40mLe2pM6tDQUM2QAYOir3HliyqcuXCVySHTGdLadDXdjfbs2UPZsmVdHYZxFxs6dCgPPPAAL7/8sqtDyRBJ/U2IyCZVDU1qfXMGkRRPb3xafUGA22nybf6cX3b+7eqIDMMw7jhTzTU5JesQW6ErfbZNo9v34ZR/8GmK5Er6zhHDMO5NQ4cOdXUILmXOIFLg3vQD1K8w7+uXvDp1DdGxyT8VaRiGca8xCSIl3jnwaDOWonKG5n9/xadLMq6GjmEYRlZnEkRqAmogNZ6ns8fv7Fs5m5X7014Q0DAM425kEoQz6v2XuPyBfJLtG96ZuYIzl8wQloZh3PtMgnCGRzbc2nyDv1zltZgxDJqxhVhTisNwQnLlvpcvX07z5s1dHR4AgwcPZunSpWluN3DgQB588MEUK5YmJyAggLNnz94yf8yYMUyZMiXN24OUK9M6U7U2q3j//fddHUICkyCcVSAQt4aDeVQ2UDhiDqOXH3R1REYWp3eg3HdGeOedd2jYsGGa2sQ/7Vu0aFFWrlyZYbH07duXp5566rba3i0JIjY25YdvbydBZFYZ8FQThFiCRKSxiNQWkTyZEsndoGo/NKAm73p9y6wlq/jz0DlXR2RkYcuWLcPT05O+ffsmzAsJCaFWrVqAVT6jbdu2PPLII3Tp0iWhdMI777xDWFgY5cuXp3fv3gnz69aty2uvvUZ4eDgPP/wwf/zxBwBXrlyhffv2BAcH06FDB6pUqZIwSM/ixYupVq0alSpVol27dgmlIRx17949oWZUQEAAQ4YMoVKlSgQFBbF3795kj618+fI8++yzTJ8+PWH+0KFD6datG40aNSIgIIA5c+bw6quvEhQURJMmTYiOjk5Y9+OPPyY8PJzw8HAOHjyY0H748OEAHDp0iCZNmlC5cmVq1aqVEEv37t15/vnnqV69OiVLlkyI/fXXX+ePP/4gJCTklqfCEy+LjY3llVdeISwsjODgYL7++msAli9fTp06dWjfvj0PP/wwr7/+OlOnTiU8PJygoCAOHTqUEEPfvn2pVasWDz/8cEI13JS2W69ePTp37kxQUBBgFU+sXLkygYGBjB07NiHOq1evEhISQpcuXYiIiKB8+fIJxzF8+PCEW2/r1q3LG2+8QZ06dRgxYgSRkZG0adOGsLAwwsLCWL16dZL/d2mR7HMQIhIAvAo0AY4AkYA3Vhnu88AY4Du9lx7FTo2bG9JyNNlGV2ek21h6TSvMwoF1yJ/D29WRGan5+XX4e0fGbrNgEDRNflS0lMp9Q9JlvWvWrEn//v0ZPNgavOrJJ59k4cKFCZVSY2JiWL9+PYsWLeLtt99m6dKljBo1ily5crF9+3Z27txJSEgIAGfPnuW9995j6dKl+Pr68tFHH/Hpp58mbDs5efPmZfPmzYwaNYrhw4czbty4W9aJLwP+xBNP8MYbbxAdHY2npydgfbAvW7aM3bt3U61aNX744QeGDRtGq1at+Omnn2jZsiUAOXLkYP369QmlwROXHO/duzdjxoyhdOnSrFu3jueeey6hQGFaS5cnXjZ27Fj8/f3ZsGED169fp0aNGjRq1AiAbdu2sWfPHnLnzk3JkiXp1asX69evZ8SIEXz55Zd8/vnngHVWsmLFCg4dOkS9evU4ePAgU6ZMSXa769evZ+fOnZQoUQKACRMmkDt3bq5evUpYWBht2rThww8/ZOTIkQn1qByryibl/PnzrFixAoDOnTsnWzb9dqX0oNwwrNHf+qvqTZ2MIlIYa3yGbsCkdEVwt8lZDGn2KcFznuHJuNn0n56Dab2q4OFueuuMtIkv6w0klPWuWbMmy5YtY9iwYVy5coV//vmHwMDAhATRunVrACpXrpzw4bFq1SoGDhwIQPny5QkODgZg7dq17N69O6H+z40bN6hWrVqqcTnuY86cObcsv3HjBosWLeKzzz7Dz8+PKlWqsHjxYpo1awZA06ZN8fT0JCgoiNjYWJo0aQIkXwa8U6dODBo06KZ9REVFsWbNGtq1a5cw7/r16wnvb6d0uaPFixezffv2hLOPCxcucODAAby8vAgLC0uohVSqVKmED/igoCCWLVuWsI327dvj5uZG6dKlKVmyJHv37k1xu+Hh4QnJAawaVPFF+f766y8OHDhAnjxp66Dp0KFDwvulS5feVEo8vmy6n59fmrbpKNkEoartAUTEA0h8FeqMqg6/7b3e7YLbw8Hf6L99FisiyjJ8cS5eb5p6YTTDhVL4pp9ZUir3DUmXsr527RrPPfccGzdupGjRogwdOvSm0s/xbRxLXyd3Eq+qPProozd1ATkjqX04+uWXX7hw4UJCV8mVK1fw8fFJSBCOZbM9PT0TylinVEY7qRLgOXPmTLaya1pLlyemqnz55Zc0btz4pvnLly9Pdxnw5LbrWLJ7+fLlLF26lD///BMfH5+EcR0S8/DwuOkmgJTKgKdUNv12OfO1d72T8+4vzYYjuYoz7oExTFuxjSW7s9aFR8P1Uir3nZz4D4C8efMSFRWVYoKJV7NmTWbNmgXA7t272bHD6kqrWrUqq1evTujfv3LlCvv3p/9hz+nTpzNu3LiEMuBHjhxh8eLFaR5fOb4M+MyZM285s8mRIwclSpRIGKBIVdm2bVuK20trGfDRo0cnXBPZv38/ly9fTlP8s2fPJi4ujkOHDnH48GHKlCnj9HYvXLhArly58PHxYe/evaxduzZhmaenZ0L7AgUKcObMGc6dO8f169eTHfkP0l423RnJJggRyS8iFYDs9kXqYPtVEzBFibL5Qdvx+Mf+w+gck3hp1haOnUvbH4hxb0up3Hdy4ocMDQoKomXLloSFhaW6n+eee47IyEiCg4P56KOPCA4Oxt/fn3z58jFp0iQ6depEcHAwVatWTfais7OuXLnCr7/+mnC2ANa32Jo1a7JgwYI0bev69etUqVKFESNG3HJRGWDq1KmMHz+eChUqEBgYyLx581LcXkqlyxMv69WrF+XKlaNSpUqUL1+ePn36pPlOoDJlylCnTh2aNm3KmDFjEkb8c2a7TZo0ISYmhuDgYN56662bRqDr3bs3wcHBdOnSBU9PTwYPHkyVKlVo3rx5iiXc01o23RnJlvsWkR5ATyAEcExFl4CJqjo7yYYulGHlvtNi9QhYMpi36c2GPC34vm91vD3d72wMRpLul3LfsbGxREdH4+3tzaFDh2jQoAH79+/Hy8vL1aHds7p3737TWOB3i7SW+07pGsREYKKItFfVWRkb5j2k2gA4tIw3j06h8cmHeHtBTj5oHeTqqIz7yJUrV6hXrx7R0dGoKqNHjzbJwcgQKd3m2klVpwOFROT5xMtV9YtMjexu4eYGrb7GfXR1pmUbS631+alYLCftQ4u6OjLjPuHn58cdP3O+z02aNMnVIdwRKV2kzmX/mxfIl8QrVSIySER2ichOEZkuIt4i8rGI7BWR7SIyV0RyJtM2QkR2iMhWEcnav/1+BaDlaPJfOcgXeebw5o872fbXeVdHZXB7d7gYxr3odv4WMm3IURF5EFgFlFPVqyIyC1gEnAR+V9UYEfkIQFVfS6J9BBCqqrcWbEmGS65BOPrlDVj7FW94vsoyt6osGFCTvA9kS72dkSmOHDmCn58fefLkueWWRMO4n6gq586d49KlSzc9iwG3eQ3CofFDwFdAQVWtICLBQDNV/cCJuDyw7oKKxrrz6aSqLnZYvha4u67ypKThEDj2J++eHUPjy4XpN3Uz3/Wqgqd5iM4lihQpwvHjx4mMNCXaDcPb2zvhwUxnpXoGISLLgTeAr1S1olhfxXaqamCqGxcZCPwPuAosVtUuiZYvAGaq6ndJtD0C/Aso8LWqjk1mH72B3gDFihWrfPTo0dTCylz/HoWva/FvtgepevpVOtd4mCGPp/qjMgzDcImUziCc+Wrrq6pr4ifs2kvRKawfv9NcwBNACaAw4CsiXR2W/xeIAaYms4kaqloJaAr0E5HaSa2kqmNVNVRVQ/Plc+rSSObKVRxajiHXhd18V3Q+E1dHMHfLcVdHZRiGkWbOJIhzIlIC65s8ItIS+NuJdg2BI6oaqarRwBygur2NbkBzoEtyxf5U9aT97xlgLhDuxD6zhkceg+rPExY5h0EFt/H6DzvYeeKCq6MyDMNIE2cSRH9gPPCIiBwFXgeedaLdMaCqiPjY3VINgD0i0gR4DWihqkk+eiwiviLiF/8eaATsdGKfWUeDwVC0Ks9f/pKK2c/Q59tN/HP5hqujMgzDcFqqCUJVD6pqfaAQUEFVq6rqESfarQO+BzYDO+x9jQVGAn7AEvsW1jFgVYgVkUV28wLAKhHZhlX36SdV/SXth+dC7p7QdgLimZ2JD3zFpaiLPPvdJm7EpH30LcMwDFdw5iK1F9ASCMDhridVzTrj4tlcfptrUg7+Bt+14WixltTZ345O4cV4v1V5c9ulYRhZQnovUs8FOmAlh1iHl+GMhxpA7Vcofmwuo8vtYvr6Y0z508V3WhmGYTgh1ecggOKqWj711Yxk1X0djq+nydHh9Cr5Ce8s3E3JfL7UKp0F7royDMNIhjNnEGtFpFymRzYHitIAACAASURBVHIvc3OHNhOQBwrwRtQHVM4TQ7+pmzkceev4wIZhGFlFSuNBbBGRzUAVYKtdU2mzw3wjLXzzQIdvcbtyjin+Y8jmpvSavJELV1N9pMQwDMMlUupiundKYGQVhUOg+Wd4//gs88uVpdbWBgyYvoUJ3ULNmNaGYWQ5yX4qqeqh+BdwGGugoKsOL+N2hHSGsGcotHsck8OPs3J/JO/9tMfVURmGYdwi1a+tIvIcEAn8Afxmv5Zmclz3tsbvQ9Gq1Ng1hNcrxTJpTQSTVqf6aIlhGMYd5Uy/xotAWVUto6pl7Ze5aJ0eHl7QfjJk86PPqcE8/rAv7yzczW97Trs6MsMwjATOJIjjwD+ZHch9x68gtJ+CXPiLzz2+oHwhXwZM32JqNhmGkWU4kyAOAr+LyCsi8nz8K7MDuy8UqwqPfYz74d+YXuJncmb35OnJGzh1wVziMQzD9ZxJEKeAlUAO0jjkqOGE0J4Q3hvfTWP4vtoRLl+PpeekjURdj3F1ZIZh3OecHnJURLKp6vVMjiddsmQtJmfExsB3reHYn2yp/y1tf4qjdum8fPOUuf3VMIzMla5aTCISLiI7gAP2dAUR+TKDY7y/uXtAu0ngX4SKa/ozvFFulu2L5O0Fu29roHHDMIyM4MzX0y+wBvc5B6Cq24B6mRnUfcknN3SaCTHXabXnZfrXKMi3a48yesUhV0dmGMZ9ypkE4aaqicuPmmqumSHfw9B2ApzZxUuXP6NlhYIM+2UfP2wyQ5YahnHnOZMg/hKRcEBFxF1EXgD2Z3Jc96/Sj8Kj7yJ7F/BJvp+p8VAeXvthOyv2R7o6MsMw7jPOJIhnsR6WKwacBqri3JCjxu2q1g8qdsV91ceMCzlE6QJ+PPvdJnYcN89IGIZx56SYIETEHeioqh1VNa/96qiqZ+9QfPcnEWj2GZSoTfZFA5na8Aa5fLzoMWk9R89ddnV0hmHcJ1JMEKoaC7S5Q7EYjjy8oP23kLskuRf0YFqr3MTEKd0mrOdcVJa+29gwjHuEM11Mf4jICBGpJiLB8a9Mj8yA7Dmhyyxw96L4z08xuUMJTl24Rs9JG7hsHqQzDCOTOZMg6gCVgGHAV/ZrZGYGZTjIFQCdZkDUaSr88Syj2pdjx4kL9Pl2E9djzM1khmFknlQThKrWSuJV25mNi8ggeyS6nSIyXUS8RSS3iCwRkQP2v7mSadtERPaJyEEReT2tB3ZPKRIKrb+B4xtpsPctPmpdnlUHz/LCjK3ExpkH6QzDyBxO1XEQkcYi8qKIvBH/cqLNg8DzQKiqlgfcgY7A68Bvqloaa2yJWz787YvjXwFNgXJAp/t+XOxyLaDRu7B7Hu3Oj+fNZmX5eeffvDFnh3na2jCMTOFMqY1RQDesW12zA12Bh5zcvgeQXUQ8AB/gJPAEMNlePhlomUS7cOCgqh5W1RvADLvd/a1afwh9GlaPoJfXUgbUf4iZG//iw5/3ujoywzDuQc6cQdRU1c7AOVV9C6gCFEmtkaqeAIYDx7Aqwl5Q1cVAAVU9Za9zCsifRPMHgb8cpo/b824hIr1FZKOIbIyMvMcfJhOBpsOgTDP4+VVeLLybJ6sW5+uVhxm93JTkMAwjYzmTIOIHJ7gmIgWBa0BAao3sawtPACWAwoCviHR1Mi5JYl6S/SiqOlZVQ1U1NF+++6AKubsHtB0PRcORub15O/gfWlQozEe/7GXaumOujs4wjHuIMwniZxHJiXU2sBWIAL53ol1D4IiqRqpqNDAHqA6cFpFCAPa/Z5Joexwo6jBdBKt7ygDwzG7d2ZS7JG4zu/BJbTfqlcnHf3/cwfxt5sdkGEbGcOYupqGqel5VZ2OdDQSpaqoXqbG6lqqKiI+ICNAA2APMx7qmgf3vvCTabgBKi0gJEfHCurg934l93j98ckPXH8DrATynt2d0s3yEBeRm0Myt/LLzb1dHZxjGPSDZBCEi7zq8rw+gqldV1anxqVV1HdaZxmZgh72vscCHwKMicgB41J5GRAqLyCK7bQzQH/gVK6nMUtVdaT66e51/EXhyDsRcxXtmOya0L0lwEX8GTN/Msn1JnZgZhmE4L9kR5URks6pWSvw+K7trR5RLr6N/wrctoUB5LrT/gc6Td3DwTBQTu4dR/aG8ro7OMIwsLF0jyhl3geLVoM14OLkZ/3nd+bZbCAF5fHl68kY2Rjh1wmcYhnGLlBJEfhF5XkQGOrxPeN2pAA0nlW0OLb6Ew8vI/XNfvu1ZiUL+3nSfuIFtf513dXSGYdyFUkoQE4F8QF6H944vI6up2BWafAh7F5L/91eY2iuMXL6ePDVhPbtOmrEkDMNIG4/kFtgPxRl3m6rPwrWLsPx9CmXzY9rTQ+gwdi1dxq1jWq+qlCucw9URGoZxlzDXIO5FdV61ynKs/5qi2z5neu+q+Hi603ncWnafvOjq6AzDuEuYBHEvEoFG70HFJ2HlMIrvm3BTkjDdTYZhOMMkiHuVCDw+Asq1hMVvUvzITGb0roaPpztdxq0zScIwjFQ5U821v4jksN9/LSLrRaRB5odmpJubuzWOROnGsHAQxY7+YJKEYRhOc+YMoreqXhSRRlgVVZ/FGl3OuBt4eEH7KVCqAcwfQLG/5t2UJHaeMEnCMIykOZMg4h+1bgpMVNVNTrYzsgpPb+g4FUrWgXnPUezET8zoXQ1fLw86f7OWreY5CcMwkuDMB/02u0bS41iVXR8gmdLbRhbmmR06TofiNWBub4qd+oUZvauS08eLruPWsf6IeeLaMIybOZMgegBDgXBVvQJ4A09nZlBGJvHyscqEF60CP/Si6N9LmdWnGvlzZKPbhPWsOnDW1REahpGFOFPuOxYoiXXtAaxhR00X090q2wPQZTY8WBm+70HBU78xs3c1iufxoefkDfy+97SrIzQMI4tw5i6mkUA9rLGoAS4DYzIzKCOTZfODrt9DoQow6yny/WV1Nz1S0I8+327i5x2nXB2hYRhZgDNnAtVVtQ/WUKPY40F4ZWpURubz9ocn50LhSjC7BzkPL+C7XlUILpKTftM2M3fLcVdHaBiGizmTIKJFxA37wrSI5AHiMjUq487w9rcGHLKvSeTYN4cpPcOpWjIPL87axpQ/I1wdoWEYLuRMgvgK+AHIJyJvA6uAjzI1KuPOie9uKl4D5vbBd/cMJnQPo2HZAgyet4sRSw+Q3KBShmHc25Kt5hpPVaeIyCagISBAO1XdmemRGXeOly90ngUzu8C8fnjHRjO6Szden7ODz5bu598rNxjcvBxubuLqSA3DuINSTRA2P+BfO1nkEZFiqnosMwMz7jAvH+s5iVlPwcIX8Ii9wbA2vfHP7sn4VUe4cDWaYW2D8XQ3N7AZxv0i1QQhIm8CNYBSwBSs5yCmATUzNzTjjvP0hg7fwfc94OdXcbt+iTcfe5Hcvl58/Os+Ll2LZmTnSnh7urs6UsMw7gBnziDaAhWBzQCqeiK+eF9KRKQMMNNhVklgMFANKGPPywmcV9WQJNpHAJeAWCAmuUG1jQzm4QXtJsG8fvD7u8i18/R79F38s3vy1rydPDV+Pd90C8U/u6erIzUMI5M5kyCuq6qKSPxdTD7ObFhV9wEhdht34AQwV1U/j19HRD4BUqoWV09VzeO9d5q7J7QcA9lywJov4doFujb/HP/snrw4ayvtx/zJpJ5hFPLP7upIDcPIRM50KM8Rka8AfxHpASwGJqRxPw2AQ6p6NH6GiAjQHpiexm0Zd4KbGzz2MdR6GTZPgR+e5vHAvEzuEc6J81dpPWoN+/6+5OooDcPIRM6U2vgIWAjMByoA/3M8C3BSR25NBLWA06p6ILldA4tFZJOI9E5uwyLSW0Q2isjGyMjINIZlpEgEGrxljU63ay7M6ET1Yj7M6lON2Dil3Zg1rDt8ztVRGoaRSSSle9ztrqFFqtr4tncg4gWcBAJV9bTD/NHAQVX9JJl2hVX1pIjkB5YAA1R1ZUr7Cg0N1Y0bN95uqEZKNk+BBQOth+o6Tef4tWx0n7iBY+eu8FmHEJoFF3J1hIZh3AYR2ZTcNd4UzyDsQn03nLkonYKmwOZEycEDaM3NF7ET7/uk/e8ZYC4Qno4YjPSq9BS0nQAnNsGEphRx+5fv+1YjuIg//advZuLqI66O0DCMDObMNYgorDEhvhaRT+NfadhHJ27tXmoI7FXVJAv+iIiviPjFvwcaAebhPFcLbAVdvocLx2H8o+SMOsR3varQuFxB3l6wm3cX7iY2zjx1bRj3CmcSxFLgPWA9sMvhlSr7jqdHgTmJFt1yTUJECtsDEwEUAFaJyDZ7vz+p6i/O7NPIZCXrQM+fIS4WJjTG++Q6vupSiR41Ahi/6gh9vt3ElRsxro7SMIwMkOI1iLuNuQZxB/17FL5rA+ePQZtxUK4Fk9dE8PaCXZQrnIPx3cIokMPb1VEahpGK274GYTfeIiKbE72WicjHIpI748M17gq5ikPPX6FQsFWeY/03dKsewLhuoRyOvEzLr1az59RFV0dpGEY6ONPFtAT4DWuY0aft6bXAv8CkTIvMyPp888BT8+HhxrDoZVgymPoP52N232rEqdJ29BqW7Tvj6igNw7hNzg4Y9IqqbrFfrwG1VPV9oEQmx2dkdV4+0GEqhPaE1SNgdjcC83kxr19Niufx5elJG5i8JsKUDDeMu5AzCcJPRCrHT4hIJSD+tldzNdIAdw9o9qn1QN2eBTCpOQXdLzK7bzXqP5KfIfN38cbcndyIMeNMGcbdxJkE0Qf4VkQOiMhB4Dugj3376bBMjc64e4hA9QHQ4Vs4vQvGNcD3wkG+fjKUZ+uWYvr6Yzw5fh3/XL7h6kgNw3CS03cx2UONSlYunmfuYsoiTmyCaR0h5jp0mAIl6zJ3y3Fe+2EHBXJkY9xTYZQp6OfqKA3DIP13MeUTka+Byap6VkTKiUj3jA7SuIc8WBme+Q38H7Ruhd00iVYVizCzd1WuRcfRetRqlu4+nfp2DMNwKWe6mCYBK4Ci9vQB4KXMCsi4R+QsBj1/gRJ1rBpOi16l4oN+zO9fg5L5HuCZbzfy1bKD5uK1YWRhziSI/Ko6DYgDUNVorEF8DCNl3v7WWNfV+sP6r2FqGwp5XmVWn2o0Dy7Mx7/uo9+0zURdN/c6GEZW5EyCuGw/EBc/YFAY1khvhpE6dw9o/D944is4ugbGNSD7hYN80TGE/z5Wll92/k2rr1Zz5OxlV0dqGEYiziSIl4EFQEkRWYFVQ2lApkZl3HsqdoVuC+F6FHzTADmwhGdql+Tbp6twNuo6LUau4rc95rqEYWQlzgwYtBGoB9QBBgLlVHVrZgdm3IOKVYHeyyB3CZjWHlZ9Ro1SeVgwoCbFcvvw9OSNjFh6gDhTEdYwsoRkb3MVkRYpNVTV+ZkSUTqY21zvEjeuwLznrFHqyraAlqO45ubDG3N3MGfzCRqWzc8n7ULw9/F0daSGcc9L6TbXlBLEt/bbvEB1YBkgWGcSK1T1iUyINV1MgriLqMKfI2HJYMhTGjp8h+YtzeQ1Ebz30x4K5fRmdJfKlH/Q39WRGsY97baeg1DVJ1X1SSAaq1uppZ0UAjElNoz0in/y+skf4cpZ+KY+snch3WuUYGafasTEKq1Hr2HaumPmVljDcBFnLlKXVNUTDtMngTKZFI9xvylZB/qshLylYWZXWPo2lYvm4Kfna1GlRG7emLuDF2dtM4MQGYYLOJMgVorITyLSVUS6APOBlZkcl3E/8S8CPX6GSt1g1acwtR253a8yqUc4gxo+zI9bT9Dyq9UcPBPl6kgN477iTILoh/U0dRWgKjDFnmcYGcfTG1p8AY+PgCMrYEIT3C8eZ2DD0kzpGc7ZqBu0GLmKuVuSHMbcMIxMYIYcNbKew8th5lNW0ug8EwpX5NSFqwycvpX1Ef/QplIR3nkiEN9sHq6O1DDuerd1kdoeVvRZESmcaL6HiNQWkfEi0iOjgzUMStaFpxeDezaY+Bjs/YlC/tmZ9kwVnq//EHO2HOfxkavYfdIMaWoYmSmlLqZmgCcwV0SOi8h2ETkAHAZ6AKNVdWJyjUWkjIhsdXhdFJEXRGSoiJxwmP9YMu2biMg+ETkoIq+n5yCNu1D+R6yKsPkegRldYO1oPNzdeLFRGab2qkLUtRhajlrNt3+a0eoMI7M41cUkItmA/MDV2xkPQkTcgRNY1zF6AFGqOjyV9fcDjwLHgQ1AJ1XdndJ+TBfTPejGFZjzDOxdCOG9ofH74O7J2ajrvDx7G8v3RdI4sAAftQkmp4+Xq6M1jLtOusaDAFDV66r6VzoGC2oAHFLVo06uHw4cVNXDqnoDmAFkuQfzjDvAywfaT7Erwo6FiU3h3wjyPpCNCd3C+O9jZfltzxmafP4Haw5l2bGsDOOu5FSCyAAdsYr8xetvd1lNEJFcSaz/IPCXw/Rxe94tRKS3iGwUkY2RkZEZF7GRdbi5WxVh206EyH0wphbsnIObm/BM7ZLMea46Pl7udBm3jg8W7TFjXxtGBsn0BCEiXkALYLY9azRQCggBTgGfJNUsiXlJ9oWp6lhVDVXV0Hz58mVAxEaWVb419P0D8j4M3/eA+QPgxhWCi+Rk4fM16RhWjK9XHqbVqNUcPGMq0htGejmVIESkiIjUs99nExHfNOyjKbBZVU8DqOppVY1V1TjgG6zupMSO8/8j2AEUwXqC27jf5QqwRqqrOQg2fwtj68LpXfh4efBB6yDGPlmZk+ev0vzLVeYCtmGkkzNjUvfEenp6nD2rODAvDfvohEP3kogUcljWCtiZRJsNQGkRKWGfgXS0YzAMcPeEhkPhyblw7TyMrQdrx0BcHI0CC/LrC7UJL5GHt+btouekDZy+eM3VERvGXcmZM4jnsZ6gvgigqvux7mhKlYj4YN2JNMdh9jAR2SEi27HGmRhkr1tYRBbZ+4gB+gO/AnuAWaq6y6kjMu4fpepB39VWPadfXoPvWsGFE+TP4c2k7mEMfbwcaw6do9FnK5m/7aQ5mzCMNEr1NlcRWauqVUVki6pWtG9B3aqqQXcmROeZ21zvU6qwaRL8+oZ1dtH8MyjfBoBDkVG8NGsbW/86T7OgQrzbsjy5fc3tsIYRL723ua4WkVcBb/s6xExgYUYGaBjpIgKhPaDvKvsCdk/4/mm4+i+l8j3A932r8UrjMize/TeNPlvJ0t1maFPDcIYzCeJV4BKwF2vI0d+A/2ZmUIZxW/KUgh6/QL03YfePMKo6HPodD3c3+tV7iHn9apL3AS96TdnIy7O3ceFKtKsjNowsLcUuJrs7aYKqdrtzId0+08VkJDi5Beb0hrP7odJT0Og98PbnRkwcX/x2gNErDpHH14v/tQri0XIFXB2tYbjMbXcxqWosUEhEzODAxt2lcEXo8wfUeAG2fAejqsGBJXh5uPFy4zLM61eD3L5ePDNlI89P38I/l2+4OmLDyHKcuUg9BuuhtnnA5fj5qvpF5oaWduYMwkjS8U0w7zmI3AshXaynsrPn4kZMHKOXH2LksgPk8PbknSfK81hQQUSSek7TMO5N6b1IHQksAXyAfA4vw7g7FKlsDWta62XYNgO+qgp7F+Hl4cbAhqVZMKAmD+bKTr9pm+n73Sbz3IRh2MyAQcb95eQW+LEfnNkFZVtA02GQoxAxsXF888cRPl+6Hy93N15r+gidw4vh5mbOJox7W0pnEM50MS0hiTpIqtooY8LLOCZBGE6JuQF/fgnLPwKPbNBwCFTuCW5uHDl7mf/O3cGaQ+eoXDwXH7QO4uECfq6O2DAyTXoTRBWHSW+gDXBdVV/JuBAzhkkQRpqcOwQLB1ljYBcJt8bDLlAOVWXO5hO899Nuoq7H0LdOKfrVewhvT3dXR2wYGS5dCSKZDa5Q1TrpjiyDmQRhpJkqbJ8Jv/wHrl+E6s9D7VfAy4dzUdf53097mLPlBCXy+vLuE+WpWTqvqyM2jAyVrovUIpLD4ZVTRBoAhVJrZxh3BRGo0BH6b4Sg9rDqU/gqHPYsII+vF592COHbp8OJU6Xr+HX0n7bZXMQ27hvOdDH9hXUNQoAY4AjwtqquyPzw0sacQRjpFrEaFr1iXcQu1QAe+xjylOJadCxjVhxi1PJDeLm78ULD0nSvHoCH+50ac8swMkd6r0F4qmp0onkedsXVLMUkCCNDxMbAhm9g2fsQcw2qD4BaL4GXL0fPXWbwvF2s2B/JIwX9eK9leUIDcrs6YsO4bel9DmJdEvPWpy8kw8jC3D2g6rNWt1Nga/jjE/iqCuz8geK5fZjUI4wxXStx4Wo0bcf8yUuztnHmkul2Mu49ySYIEckvIhWA7CISJCLB9qsm1kNzhnFv8ysArb+GHj+Dd06rSuyEJsjJLTQpX4ilL9ahb51SzN92gvrDVzB25SEzHrZxT0m2i0lEegA9scpsbHVYdAmYqKqzk2zoQqaLycg0cbFWTaff34XLkVChMzQYDDkKceTsZd5duJvf956hZD5fBjcvR90yTo2pZRgul95rEO1VdVamRJbBTIIwMt21i1aX09pR4OZpjY1drR94+bBs7xneWbibI2cv07Bsft5sVo6AvGkZvt0w7rx0PwchIo2BQKwH5QBQ1fczLMIMYhKEccf8cwSWDIY988GvMNR7A0I6cyNOmLj6CF/8doAbsXF0qxbAgPql8fcxBZGNrCm9ZxCjgJxAbWAi1pPUa1W1Z0YHml4mQRh33NE1sPgtOLER8pWFhkPh4cacuXSdTxbvZ9amv8iZ3ZMXGj5M5yrF8DS3xRpZTHoTxHZVDRaRbapaQUT8gB9MLSbDsKlaZxJL34Z/DkHxGvDoO1AklN0nL/LeT7tZc+gcJfP58t/HylL/kfympLiRZaT3Ntf4+/euiUhBezrAiZ2WEZGtDq+LIvKCiHwsIntFZLuIzBWRnMm0jxCRHXZb86lvZF0iUO4J6LcOmn1ijWI3rgHM7Eo5jxNM7VWFcU+FgsLTkzfSdfw6dhy/4OqoDSNVzpxBDAU+Bx4FvgRigcmq+obTO7GGLj0BVAHKAL+raoyIfASgqq8l0SYCCFXVs87ux5xBGFnC9Sj4cySsGQk3oiCoHdR9neicJfhu7VG++O0A/16JpkWFwrzSuAxFc5u7xg3Xue0uJhFxA8JUdZ09nR3Irqr/pDGARsAQVa2RaH4roK2qdkmiTQQmQRh3syv/wOoRsO5riL0BFbtA7Ve56F2Qr1ccYvyqI8TGKV2rFmdA/dLk9vVydcTGfSi91yDWqmrVdAYwAdisqiMTzV8AzFTV75JocwT4F6sO1NeqOjaZbfcGegMUK1as8tGjR9MTqmFkvEunrSKAGydY05W6Qc1B/E0ePl+6n1kb/8LXy4O+dUvRo0YAPl4ero3XuK+kN0G8C2xU1Xm3uXMv4CQQqKqnHeb/FwgFWmsSQYhIYVU9KSL5sYY8HaCqK1PalzmDMLK0C8dhxTDYOhXEDSp2hZqDOHA9Fx/9so+le06T94Fs9K9Xik5VipHNw4w/YWS+9CaIfwF/4DpwFauqq6qqUxXKROQJoJ/jXU8i0g3oCzRQ1StObGMoEKWqw1NazyQI465w/his+gw2fwsohHSGmi+y6VJOPv51L2sP/8ODObMzsEFpWld60FSMNTJVehNEkl9jVDXWyZ3PAH5V1Yn2dBPgU6COqkYm08YXcFPVS/b7JcA7qvpLSvsyCcK4q1w4bl2j2DQZ4mIguANaYyCrL+Tl41/3su34BUrm9WXQow/TLKiQGR/byBQZ8SR1R6Ckqr4vIkWAAqq6yYl2PsBfdtsL9ryDQDbgnL3aWlXtKyKFgXGq+piIlATm2ss9gGmq+r/U9mcShHFXungK1nwBGyda5cUfaYbWeIHFF4vyyeJ97D8dRZkCfgxsWJomgQVNojAyVHrPIEYCnkBtVS0rIrmxzgjCMj7U9DEJwrirXT4L68dadz1dOw8BtYitMYiFUWUY8ftBDkde5pGCfrzQsDSNyplEYWSM9CaIzapaSUS2qGpFe942Va2QCbGmi0kQxj3hehRsnmw9R3HpJBQMJq7aABbEhDNiWQSHz16mbKEcDGxQmkblCphEYaRLep+kjrafh1B7Y3kAU/TeMDJLtgesCrEDt0GLkRBzDbe5z/DE8qYsqbKFL1sGcPVGDH2/28RjX/zBwu0niY1LvavYMNLKmTOIp4BWWLekTgDaY41JPSPzw0sbcwZh3JPi4uDgUvjzSziyEjx9iQvpwmK/1gzbcJ3DkZcpmc+XfnUfokVIYVMQ0EiTjLhIHQg0tCd/U9WdGRhfhjEJwrjnndpujUWx43uIi0HLNGVdvnYM3ZGHvaejKJo7O8/WeYg2lR80z1EYTsmIBBEM1MTqZlqtqtszNsSMYRKEcd+4eAo2fAObJsGVc2jeMuwr1pEhx4JZd/w6BXJko2eNEnSuUgw/bzMWhZG89F6k/i/QGeu2UwGeAKaq6gcZHWh6mQRh3Heir8GuOdadT6e2otlycDKgNZ/8W4s5x7Lj5+3Bk1WL06NGCfL5ZXN1tEYWlN4EsQeoHP/Es/1swyZVLZvhkaaTSRDGfUv/r707j6r6PBM4/n3YNG4oqysobsRdUapRo8a41EzNNE3aZLpOO81pJ02X6TLtnNPGTKdzMk2b6ZJJ57Rp06ZLYtMsTZNo0o4xVRM3QMV9BUVAQFBRFIH7zB/vj+SWXODCvRe48HzOuYf729+XV+/D7/e+93kVSna7YbIHngdfA7UjFrBeV/Dd4gkQm8CdOaO5d3GWTYNq/kaoAWIj8EFVveQtDwGeUtXbwl7SEFmAMAa4XAEFv3aPny6cpvGGFLYOWsW3S3M56Utl5ZR0Pr04i5zMYTZxkQk5QDwHzANexfVBrAS2AucAVPVfwlraEFiAMMaPzwcnNkHeE3BkA2gTRYm5PHbxJl64NpspY1L59OIsVk1Nt3xPfVioAeJTbW1X1Z+HULawsgBhTCsulUL+k1DwG7h4hvr4ZJm4CgAAE85JREFURF7SRTx+ZRG1Q7P55MJx3DV3tHVo90Ehj2KKFhYgjGmHrwlOvQH5v0YPv4Q0XedE3ESeuLqITXGLWZmTzcdvGss466foM0K9g1gNfBvIxCXO61C6765kAcKYDqirhsJnXNrxc4U0Escm32z+0LgI34QVfGzRJBZPTLF+il4u1ABxHPft6UL8UmwEm+67K1mAMKaTyvbBvvU07V1PbF0lFxjEi40L2DFkBbkLV3JHjj1+6q1CDRCbgVtUtcfnX7IAYUyImhrh5Os0FfwODr9MrK+eIl86G2Qh9dnvZ/WypWQPH9LdpTRhFGqAyAUeADbjZpUDQFV/FMYyhoUFCGPC6NpFOPgitbufZmDpm8Tg45Avg/why0m76cMsyc0hIc5GP0W7UAPEBqCBdz9i+mY4CxkOFiCMiZDac9QVPEPt7qdJv1QIQCETqchYzaSlH2VM1uRuLqDprFADRJ6q5kSkZGFmAcKYyPOdP0XRG78m7vAfybh+HIDj8ZNpmLyW8Us/QkLK2O4toOmQUAPEd4GNqropEoULJwsQxnStquJDHH39NyQVv0K2ngSgbOCNJExbS/LcD0Cq3Vn0dKEGiBogEagDrmPDXI0xLTT5lJ35uyl982myqjYzO8bdWVwaOI7+M9aSMPV2GDkbYqzPoqcJNUAETCpvw1yNMYFU1F7j1bcKqM57jpy6bcyPOUSc+Lg+IJ347Pcik98LWUsg/obuLqohPPNB3A1kqep/ishoIF1V89o5ZjKw3m9VFvAt4Elv/VigCJcIsCbA8auBHwKxwOOq+lB75bQAYUzPoarkn67hxbcOUH/wFW7W3SyN3ccAruGL609M1jKYvBomrYbBw7u7uH1WqHcQjwLxwM2qeqOIJAGvquq8DhQgFjgLvAe4D6hW1YdE5OvAMFX91wD7HwVWACXALuAeVT3Y1nUsQBjTM12ub2RDYRl/zDuJFL/J8ph8buu3h9SmCrfD8BkwcQVMWAGj50FsXPcWuA8JNUDkq+ocESlQ1dneur2qOrMDBVgJPKCqC0XkCLBUVctEZASwWVUnt9h/AbBOVVd5y98AaG+SIgsQxvR8Z6rreL7gLM/lnaFfzRFWxe9h7cCDZF07QIw2Qf9EGH+LCxbjb4EhI7q7yL1aWwEimDDdICIxuFTfiEgyft+HCNLdwFPe+3RVLQPwgkRagP1HAWf8lktwdx/GmCg3JmkAn18+kftvmUD+6Vk8m7+ID+wrQ69eYPWAw9w94AjTT24j/sDz7oC0KS5QZC2DzJsgYUD3VqAPaTVAiEicqjYC/wM8C6SKyIO4vEwPBnsBEUkA1gLf6EC5AmUHC3irIyL3AvcCZGRkdOASxpjuJCLkZCaRk5nEuvdN5Y2jlbxQMIF7Ds2jvvEelg89x8fSTzG3aQ8Dd/4M3noUYhMgYwFkLYVxS2DkLIgJOI7GhEGrj5iaHy1576cCt+I+uP+iqvuDvoDI7cB9qrrSW7ZHTMaYVtVea2Dj/nL+uKeUN09U4VOYkRbPP2WUsyRuP4mlW6DC647slwhjF7pgMe5mSLsRLPtsh3SqD8K/zyHEiz+N69R+wlt+GDjv10mdpKpfa3FMHK6Tejmuc3sX8A+qeqCta1mAMKZ3qaytZ8P+Mv60t5RdRW6w4/RRidyVncCawcdJqdgOp/4KNafcAQNSXMDIXOR+pt5o371oR2cDRAnwSGsnVdVWt/mdYwCuLyFLVS9665KB3wMZwGngLlWtFpGRuOGsa7z91gA/wA1z/YWqfqe961mAMKb3Kr1wlZf3lfGnfaXsK7kIuGCxZvoI3pfZyOgLu6Boq3td9Lowb0hy/RaZCyFzAaRPtxFSLXQ2QJQBPyFwfwCqGnQ/RFexAGFM33Cmuo4N+8t4ubCcvWcuADB15BDWTB/BqqnpTIivhuJtULQNirbAhWJ3YPxAGDPP9WNkzHdDahP69ux5nQ0Qb/dBRAsLEMb0PSU1dWzcX87LhWUUnHbBYnzqQFZNHc6qqcOZMToRuVQKZ7bD6e1w+i0o3w8oSCwMnwajc2FMrgsYw8b2qX6Mbu2D6EoWIIzp28ovXuO1g+W8eqCc7SerafIpIxP7s3LqcFZMSSd3XBLxsTFuroszu1ywKNkJZ/Ph+mV3koFpLliMyoHRc10OqX6Du7diEdTZAJGkqtURLVmYWYAwxjSruXKd/ztcwcb95Ww5Vkl9o4/B/eNYNjmNW6eks3RyKkOap1H1NbmRUWd2uMBRshOqT3pnEkjN9gJGDoyc476bEZfQbXULp5BzMUULCxDGmEDqrjey9VgVfz54jk2HKzh/5TpxMcL8rGRuyU5j+Y1pZCa36Iuoq3Z3Fmfz4OxuKNkNV72/mWMTIH2au7tofqVmR2UHuAUIY4zxNPmUgtM1/PnQOf5y8BwnKq8AkJU6kOXZaSzLTmPeWO9RlD9VN5y2dA+UFrhX2V6ov+S2x/V3dxYjZsKIGe5n2lSI79/FNewYCxDGGNOK4vNX2HS4gk2HK9hxsprrTT4G94tj0cQUlk5OZcmkNIYntvIh7/O5R1Gl+S5YlO2Fsn1Q74bhIrFu0qT0aTB8uusQT58Og1K7roLtsABhjDFBuFLfyNbjVbx+uILNRyopv3QNgOzhg1kyOZUlk1KZm5lEQlwbX75TdcNqy/a5gFFeCOf2w6Wz7+wzKN0FjfQp7i4j7UYXSLphjgwLEMYY00GqypFztbxxpJLNRyrZXVxNQ5MyICGWBVnJLJ6YwuJJqWSlDESCGRZbV+0CRfl+72chVB6Bpnq3XWIgabwLFmlTXMBIu9Gti2CHuAUIY4wJ0eX6RrYdr2LrsSq2HKuk6HwdAKOG3sDiiSksmpjCTeNTSBrYgQ/zpkb3iKriAFQcgnMH3Giq6lO8nZ80Js4LHNmuIzxlknslTwhLZlsLEMYYE2anz9ex5XglW45Wse1EFbXXGgGYMmIICycks3BCCrnjkhiQ0ImRTQ1XoeoYVB52gaPyCFQegpoi0ObZFgSGjnHBIjUbVv5Hp77gZwHCGGMiqLHJx76zF3nzeBXbjp8nr7iG600+4mOF2WOGMX98MguykpmdMZT+8SGkJ2+4BtUnoOooVB51P6uOuM7yf36zU6e0AGGMMV3o6vUmdhdXs+34ed46UUXh2Yv4FBLiYpiTMZQFWSksGJ/MzDGJ9IsLw3wWqp1OD2IBwhhjutGlaw3sOlXNWyfO89bJ8xwsu4Qq9IuLYXbGUHLHJfOecUnMzhjauUdSIbAAYYwxPciFuuvsPFXNzlPV7DhVzYFSd4cRFyPMGJ3IvHFJzMtMIidzGMM60undCRYgjDGmB6u91kBecQ07vKCxr+QCDU3us3li2iDmjk1i3thhzM1MYkzSDcENqw2SBQhjjIki1xqa2FdykV1F1ewqqiavqIbaejdKKmVQP3Iyh5KTOYyczGFMHZkYUsd3WwEi+jJLGWNML9c/PpbccUnkjksCXP6oI+W15J+uIb+4hrzTNbx64BwACbExzByTyPp7FxATE955LCxAGGNMDxcbI0wZOYQpI4fwkfmZgJuvu+C0CxYX6xrCHhzAAoQxxkSl1MH9WDl1OCunDo/YNdrIOGWMMaYvswBhjDEmoIg+YhKRocDjwDRc5qlPAl8EJnu7DAUuqOqsAMcWAbVAE9DYWi+7McaYyIh0H8QPgY2qeqeIJAADVPVDzRtF5PvAxTaOX6aqVREuozHGmAAiFiBEZAhwM/AJAFW9Dlz32y7AB4FbIlUGY4wxnRfJPogsoBJ4QkQKRORxEfGfFXwxcE5Vj7VyvAKviUieiNzb2kVE5F4R2S0iuysrK8NXemOM6eMiGSDigDnAT1R1NnAF+Lrf9nuAp9o4fqGqzgHeC9wnIjcH2klVf6qqc1V1bmpqz5nn1Rhjol0kA0QJUKKqO7zlP+ACBiISB9wBrG/tYFUt9X5WAM8DuREsqzHGmBYi1gehquUickZEJqvqEWA5cNDbfCtwWFVLAh3rPYqKUdVa7/1K4N/bu2ZeXl6ViBR3ssgpQG/pEO8tdekt9QCrS0/UW+oBodUls7UNkR7FdD/wW28E00ngH731d9Pi8ZKIjAQeV9U1QDrwvJexMA74napubO9iqtrpZ0wisru3DKXtLXXpLfUAq0tP1FvqAZGrS0QDhKruAd5VaFX9RIB1pcAa7/1JYGYky2aMMaZt9k1qY4wxAVmAeMdPu7sAYdRb6tJb6gFWl56ot9QDIlSXXjVhkDHGmPCxOwhjjDEBWYAwxhgTUJ8KECKyWkSOiMhxEfl6gO0iIj/ytu8TkTndUc5gBFGXpSJyUUT2eK9vdUc52yMivxCRChHZ38r2aGqT9uoSLW0yRkReF5FDInJARL4QYJ+oaJcg6xIt7dJfRHaKyF6vLg8G2Ce87aKqfeIFxAIncDmiEoC9wJQW+6wBNgACzAd2dHe5Q6jLUuCl7i5rEHW5GfcN+/2tbI+KNgmyLtHSJiOAOd77wcDRKP6/EkxdoqVdBBjkvY8HdgDzI9kufekOIhc4rqon1WWWfRq4vcU+twNPqrMdGCoiI7q6oEEIpi5RQVX/ClS3sUu0tEkwdYkKqlqmqvne+1rgEDCqxW5R0S5B1iUqeL/ry95ivPdqOcoorO3SlwLEKOCM33IJ7/6HEsw+PUGw5Vzg3Y5uEJGpXVO0sIuWNglWVLWJiIwFZuP+WvUXde3SRl0gStpFRGJFZA9QAfxZ38l11yys7RLpVBs9iQRY1zL6BrNPTxBMOfOBTFW9LCJrgBeAiREvWfhFS5sEI6raREQGAc8CX1TVSy03Bzikx7ZLO3WJmnZR1SZglrjZOp8XkWmq6t/nFdZ26Ut3ECXAGL/l0UBpJ/bpCdotp6pear4dVdVXgHgRSem6IoZNtLRJu6KpTUQkHveB+ltVfS7ALlHTLu3VJZrapZmqXgA2A6tbbApru/SlALELmCgi47zkgXcDL7bY50XgY95IgPnARVUt6+qCBqHduojIcPGyHYpILq6tz3d5SUMXLW3SrmhpE6+MPwcOqeojrewWFe0STF2iqF1SvTsHROQGvKzYLXYLa7v0mUdMqtooIp8DXsWNAvqFqh4Qkc942/8XeAU3CuA4UMc72Wd7lCDrcifwWRFpBK4Cd6s3zKEnEZGncKNIUkSkBHgA1/kWVW0CQdUlKtoEWAh8FCj0nncD/BuQAVHXLsHUJVraZQTwKxGJxQWx36vqS5H8DLNUG8YYYwLqS4+YjDHGdIAFCGOMMQFZgDDGGBOQBQhjjDEBWYAwxhgTkAUI0+VEJNkvc2a5iJz1W36zi8vysJcZ82ER+XsRmRLi+ZK97KGXReTRFtu+IyJnRORya8eHg4j8QERu9t4XdeWXvkTkeyJyS1ddz0SWDXM13UpE1gGXVfV73XT9S0CqqtaLyC9xWT3/0IHj41S10W95IC7fzzRgmqp+zm/bfKAYOKaqg8JVhxblSQJeUdX53nIRMFdVq8J4DcF9dvgCbMsEfqaqK8N1PdN97A7C9CjNf12Ly9H/hoj8XkSOishDIvJhcfnwC0VkvLdfqog8KyK7vNfCAOccKyJbRCTfe93krX8RGAjsEJEHgLXAw96dzHjvtVFE8rzjs73jfikij4jI68B/+V9LVa+o6lbgWstyqOr2QN9qFZH3icgOESkQkb+ISLq3fp2I/EpEXvPuBO4Qke969d8oLoVES3cCG1usu9+rd6FfHZJE5AVxcwZsF5EZftf8il/Z9nu/v7Hi5lR4DJe7aIz3e9jvnfdLXh2LgWQRGR6gbCbKWIAwPdlM4AvAdNy3YSepai7wOHC/t88Pgf9W1XnAB7xtLVUAK1R1DvAh4EcAqroWuKqqs1T1QVyagq96yydwE8Hfr6o5wFeAx/zOOQm4VVW/HIZ6bsXl9Z+NS93+Nb9t44HbcGmcfwO8rqrTcd/4vS3AuRYCeS3WVXl1/4lXD4AHgQJVnYH7ZvGTQZRzMi6V9GwgBRilqtO88jzht1++Vw4T5fpMqg0TlXY1/8UtIieA17z1hcAy7/2twBT31AOAISIy2Mv93yweeFREZgFNuA/3NonL/nkT8Izfufv57fKMl1kzHEYD68Xl7U8ATvlt26CqDSJSiEur0nx3UAiMDXCuEUBli3XNCerygDu894twARVV3eT1nSS2U85ib44BgJNAloj8GHiZd9oGXEAe2c65TBSwAGF6snq/9z6/ZR/v/NuNARao6tU2zvMl4BzujiSGAI9/AogBLqjqrFa2XwniHMH6MfCIqr4oIkuBdX7b6gFU1SciDX45gvx/B/6uAv1brGv+vTX5HdNaWuhG/vbJgv+53q6zqtaIyExgFXAf8EHgk37HtNUeJkrYIyYT7V4D/DuCA32gJwJlXqfqR3F/iQdSi5uWEm/OgFMicpd3XvE+ECMhETjrvf94iOc6BEwIYr+/Ah8G19+Dewx1CSjCTZuKuPmMxwU62BsZFaOqzwLfbD7GMwkIOC+3iS4WIEy0+zww1+tsPQh8JsA+jwEfF5HtuA+v1v76fxr4qtdZPB73AfopEdkLHCDIaV29kUOPAJ8QkZLmobNeB3MJMMBbv847ZB3uUdYWINTRRi/jMsq2Zx3e7w14iHcC07NAkrjMp5/FzeEcyChgs7ffL4FvwNtzL0wAdneu+KYnsWGuxvQyIrIV+DtvUpmuvvb7gTmq+s2uvrYJP7uDMKb3+TLefAfdIA74fjdd24SZ3UEYY4wJyO4gjDHGBGQBwhhjTEAWIIwxxgRkAcIYY0xAFiCMMcYE9P9sGEgbrgo8NAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(time1, temp_numerical1, '-', label = 'Constant Ambient temperature')\n",
    "plt.plot(time1, temp_numerical, '-', label = 'Change in Ambient temperature')\n",
    "plt.title('Ambient Temperature Patterns')\n",
    "plt.xlabel('Time after 11am (hours)')\n",
    "plt.ylabel('Temperature (degrees Farenheit)')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 307,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The corpse was 98.6 degrees Farenheit 43.0 minutes ago from 11:00am.\n",
      "The time of death was 10:17am\n"
     ]
    }
   ],
   "source": [
    "temp_living = 98.6\n",
    "temp_found = 85\n",
    "\n",
    "time_found_hours = ((np.log((temp_found - 60)/(temp_living - 60)))/(-k))\n",
    "time_found_mins = time_found_hours*60\n",
    "\n",
    "print('The corpse was 98.6 degrees Farenheit', round(time_found_mins,0), 'minutes ago from 11:00am.')\n",
    "print('The time of death was 10:17am')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}