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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K = 0.6111111111111112\n"
     ]
    }
   ],
   "source": [
    "#assign variables to the given data values\n",
    "temp1=85 #temperature at time 1\n",
    "temp2=74 #temperature at time 2\n",
    "tempAmbient=65 #ambient temperature\n",
    "deltaT=2 #change in time\n",
    "#Solving Newton's Law of Cooling for K and plugging in the arguments (dT = temp2-temp1, dt = deltaT, T-Ta = temp2-TempAmbient)\n",
    "K=(-(temp2-temp1)/deltaT)/(temp2-tempAmbient)\n",
    "print(\"K =\",K)"
   ]
  },
  {
   "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": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K = 0.6111111111111112\n"
     ]
    }
   ],
   "source": [
    "#create function to solve for K\n",
    "def measure_K(Temp_t1,Temp_t2,Temp_ambient,delta_t): #function accepts four arguments, one for each value of Newton's Law of Cooling\n",
    "    K=-((Temp_t2-Temp_t1)/delta_t)/(Temp_t2-Temp_ambient) #solving Newton's Law of Cooling for K and plugging in the arguments (dT = temp_t2-temp_t1, dt = delta_t, T-Ta = Temp_2-Temp_ambient)\n",
    "    return K\n",
    "\n",
    "#assign variables to the given data values\n",
    "temp1=85 #temperature at time 1\n",
    "temp2=74 #temperature at time 2\n",
    "tempAmbient=65 #ambient temperature\n",
    "deltaT=2 #change in time\n",
    "#plug variables in the function \"measure_k\" to determine the corresponding K value\n",
    "print(\"K =\",measure_K(temp1,temp2,tempAmbient,deltaT))"
   ]
  },
  {
   "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": 516,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f496ea04438>"
      ]
     },
     "execution_count": 516,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEWCAYAAACnlKo3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3hUZfrw8e+TMqmQTicJRTooRXoVQUQEKTZUihRF133tiqsL/ta1oK5lFRWpi2JFQDqCgApI7x0k9BYIkNDS7vePMxlSJmVSmCTcn+uaKzPnPOecewaSe85TjYiglFJK5YeHuwNQSilVcmkSUUoplW+aRJRSSuWbJhGllFL5pklEKaVUvnm5O4DrKTw8XKKjo90dhlJKlSjr16+PFZEIZ/tuqCQSHR3NunXr3B2GUkqVKMaYg9nt0+ospZRS+aZJRCmlVL5pElFKKZVvbk0ixphAY8wrxpitxph4Y0ysMWalMWaQMcakKzfaGCPZPJ5353tQSqkbmdsa1o0xHsB8oDUwBfgv4A88CEwC6gIvZTrsGSA207b1RRupUkqp7Lizd1YLoC3woYg8k7bRGDMW2AU8RtYkMlNEYq5bhEoppXLkzuqssvafx9JvFJFErLuNi84OMsaUNcZc1+R3/nIS09cfYemuU9fzskopVey5805kDXAOeNEYEwOsBvyAQUBT4HEnx2wBygApxpg1wL9EZH5RBrlk50lGfLWBxJRUWlQLpVOdckV5uWLr/PnzxMbGkpiY6O5QlFIF5OnpSZkyZQgNDcXHx6dA53JbEhGROGNMT2A88H26XfFAXxGZmW7bOWAcsBKIA2oDTwNzjTGPisjk7K5jjBkODAeIjIx0Oc6GVYJITk0FYE3MWU5duEK5sr4un6cku3LlCidPnqRKlSr4+fmRrs+DUqqEERGSkpK4cOEChw4dIjIyskCJxN1dfBOAbcB7QB9gKLAPmGaM6ZJWSEQ+FJHHRGSKiPwsIu8CjYCTwAfGmMDsLiAi40SkmYg0i4hwOmo/R+XK+NKiWpj9XDB/2wmXz1HSnT59moiICPz9/TWBKFXCGWOw2WyEh4cTEhLC2bNnC3Q+tyURY0xDrDuLX0TkBRGZISITsBrbTwBfGmM8szteRM4AnwPBWD28isxdjSo6ns/dcrwoL1UsXblyhcDAbPO0UqqEKlu2LPHx8QU6hzvvRJ4BfIEf0m8UkUvAXCAKiM7lHDH2n+GFHFsG3RpUwMP+BXztwbOcOH+lKC9X7CQnJ+PldUNNs6bUDcHb25uUlJQCncOdSaSy/aezuw2vTD+zc5P958lCiSgb4YE+tKqRvkrrxrsb0WospUqfwvi9dmcS2WH/OSj9RmNMMNALqwF9vzHGyxgTlPlgY0xVYARwBqtarEjd1bCS4/mNWKWllFLOuLOO4kNgAPC2vX1kBRAKDAMqAk+KSLI9qRwwxswEdnKtd9ZQIBB4UEQuF3Wwd9Qvz2uztpGSKqw7GMexc5epFOxX1JdVSqlizW13IiJyEGgOTAU6YU178jJwGKuL71h70cvAdKAZMBIYCzwELAZai8gPXAdhgT60tldpAczbqncjqmgMGjSoyKoPo6Oj6dixY5GcO72OHTvirgXgRo8ejTGGmJiYIjm/MYZBgwYVyblLIrd28RWR/SIyUESqiIi3iJQVkfYi8lO6MldFZKiINBSREHu5iiLST0TWXM9472qYrpeWJhFVTI0ePZqZM2fmXlA5de7cOUaPHs2yZcvcHUqJ4O5xIiXKHfUr4GnvprXx0DmOxF1yc0RKZfX6669nm0R2797NokWLrnNEJcu5c+d4/fXXs00ily9f5ssvv7y+QRVjmkTyIiUZ9i0mJHY9bWpe6008f+uNN/BQlWw+Pj7YbDZ3h1Gi+fr64u3t7e4wig1NIrnZtwTerw1f9YXf36NHuiqtOVqlVarEx8fz6quv0qJFC8LDw/Hx8aFmzZq8/PLLXLqU8a5z2bJlGGOYPHkykyZNon79+vj4+BAVFcWYMWOynHvRokXcf//9VK9eHT8/P4KDg+natSvLly/PNa6///3vGGPYu3dvln3Hjx/Hy8uLIUOGEBMT42hLmTJlCsYYxyNNdm0iGzdu5N5776V8+fL4+PhQtWpVHnzwQfbv3+8o891339GzZ0/HNBnh4eHcc889bNmyJdf3kJO5c+fSoUMHwsPD8fPzIzIykj59+rBnz54M5bZs2ULv3r0JCwvD19eXevXqMWbMmDyNc8ipnSl9G8eyZcuoVq0aYN3RpX1+6dt3smsTGT9+PE2aNMHPz4+goCC6du3KH3/8ke31Vq1aRYcOHQgICCA8PJyhQ4eSkJCQ63spbnQEWW7CasIl+xIm+5fSrZsX//A0JKUImw+f4/DZS1QN9XdvjKpQHD16lPHjx9O3b1/69++Pl5cXy5cvZ8yYMWzcuJGFCxdmOebzzz/n5MmTDBkyhODgYL766iteeuklqlSpQv/+/R3lJk+ezNmzZxkwYABVqlRxXKtz584sXbqUdu3aZRvXY489xn//+18mTpzIW2+9lWHflClTSElJYciQIURERDB16lQeeeQR2rVrx/Dhw/P0vufMmUPfvn0JCAhg6NCh1KxZkxMnTrBw4UK2bdtGjRo1APjkk08IDQ1l+PDhVKhQgf379zNu3DjatGnDhg0buOmmm3K5UlbLly+nZ8+eNGzYkJEjRxIcHMyxY8dYvHgx+/bto1atWgCsW7eODh064O3tzZNPPkmFChWYPXs2L730Eps3b+brr792+drO1K1blw8++IBnnnmG3r1706dPH4BcZ2x46aWXGDNmDM2bN+fNN98kPj6ecePG0alTJ2bNmkX37t0zlN+0aRM9evRg8ODB9O/fn2XLljFhwgQ8PDwYN25cobyX68WIiLtjuG6aNWsm69atc/3ACXfA4T+t593fY/D2m1m6+zQAI++sw2MdahRilMXPzp07qVu3bpbt0S/PdUM0eRPz9l0uH5OYmIgxJktVxWuvvcYbb7zB6tWrad68OWB9Y+3UqRMVK1Zkx44dBAcHA3Dp0iWioqKoWbMmq1atcpzj4sWLBAQEZDjvyZMnqV+/Ps2bN2fevHmO7YMGDWLKlCmk/91s3bo1MTExHDp0KMPsAbVq1cLLy4sdO3Y4thljGDhwIJMnT87yHqOjo4mOjnbU96fFa4xh48aNVK5cOUP51NRUPDw8sn0PO3fu5JZbbmHIkCGMHTvWsb1jx47ExMTk2kPq2Wef5YMPPuDkyZOUK5f9DNlt2rRh9erVbNiwgUaNGgHWRIL3338/P/zwA4sXL6Zz586A1bHg9ddf58CBA447CGefaXafV0xMDNWqVWPUqFGMHj061/K7d++mbt26tG7dml9//dVRXXjs2DHq1atHcHAw+/fvx9PT03G8MYaVK1fSsmVLx3nvuusuFi1aRFxc3HWdZii73+/0jDHrRaSZs31anZUXDftde75tOnc1SjfwUKu0Sg2bzeZIIMnJycTFxREbG8vtt98OwOrVq7McM3jwYEcCAfD396dly5ZZqp7S//FNSEjgzJkzeHp60qJFC6fnzWz48OEcP36c+fOvrXzw22+/sXfvXoYMGeLaG01n4cKFxMbG8txzz2VJIIAjgaR/DyLChQsXiI2NJSIigtq1a+fpPTgTFGSNI54+fTrJyclOy5w6dYqVK1fSs2dPRwIB64/xK6+8AsCMGTPydf3CMGvWLESEF198MUN7U6VKlRg0aBAHDx5k48aNGY5p1apVhgQCcNttt5GcnFxkXZOLiiaRvKh3D6TNBXloFV0rJ+LtadWvbjlynkNntJdWaTF27FgaNWqEj48PoaGhREREONoQ4uLispSvXr16lm1hYWGcOXMmw7b9+/fzwAMPEBISQpkyZQgPDyciIoJ58+Y5PW9m999/P0FBQUyYMMGxbcKECdhsNgYMGODiu7wmLdk1btw417IbN26kR48elClThqCgICIiIoiIiGDr1q15eg/O/O1vf6Nx48Y88cQThIaG0r17dz7++GNOnz7tKHPgwAEA6tevn+X4evXq4eHhwV9//ZWv6xeGnOJr0KABQJb4svt/A2T5v1PcaZtIXgRGQPWOsH8JAGX3zaL9TW1ZYl/pcO7W44zoWLqrtJzJT5VRcfaf//yH5557jq5du/L3v/+dSpUqYbPZOHr0KIMGDSLVvq5MemlVFDlJSEigffv2XLx4kaeffpqGDRtSpkwZPDw8eOutt/j1119zPYefnx8PP/wwX3zxBSdOnMDf358ff/yRnj17kp8lDtKkVe/kNrjx0KFDtG/fnrJly/Laa69Ru3ZtAgICMMbw9NNP57tBOCwsjLVr1/L777/zyy+/8Ntvv/HMM88watQo5s2bR6tWrZxWQbkqu/eX3d2PK/ITX07/b0paE4MmkbxqeK8jibB1One1uDddEjl2QyaR0mbq1KlER0czf/78DNU4CxYsKNB5lyxZwrFjx5g4cSKDBw/OsO/VV1/N83mGDx/Op59+yv/+9z+CgoK4dOlSgaqyAGrXrg1YdxldunTJttyMGTNISEjg559/plOnThn2nTlzpkCLGnl6etKxY0fHHd+WLVto2rQpb7zxBnPnznV8a9++fXuWY3ft2kVqaqrTb/bphYaGAnD27FnHc8h6hwCuT0qY1vFg+/btjudp0tqqcouvJNPqrLyqcxd42Vc0PLmVruXOYfO0Pr5tRy8QE+t0SXhVgnh6emKMyfBNMDk5mbfffrvA54Ws3zAXLVrkUltCo0aNaN68ORMnTmTChAlERkbStWvXLOUCAwPzvNBQ165dCQ8P5/333+f48azte2kxZ/cevvzyS06cyP94qdjY2Czb6tSpg5+fn+M9lCtXjtatWzN79my2bduWIba03mq9e/fO8TppvbwWL16cYfv777+fpWxao3ZeP8OePXtijOHdd98lKSnJsf348eNMmjSJqKioPFUXllR6J5JXvmWh1h2wYxYAgXtm0r5WNxbvtGahn7v1OE92qunOCFUB9evXj5EjR3LnnXfSp08fLly4wLRp0wo8sKxt27ZUqFCB5557jpiYGKpUqcKmTZuYOnUqDRs2ZOvWrXk+1/Dhwxk6dCgAo0aNynDHlKZly5YsXryYd955h8jISIwxPPDAA07P5+/vz4QJE+jXrx8NGjRwdPE9ffo0Cxcu5Nlnn6VXr17ceeed+Pv788gjj/C3v/2NkJAQVqxYwbx586hRo0a+q4WGDRvGkSNH6Nq1K1FRUVy+fJnvvvuO+Pj4DG09H330ER06dKBdu3aOLr5z5sxh4cKF9O/f39EzKzsPPvggr7zyCsOHD2fXrl2EhYUxf/58p0ksLCyMmjVr8u2331KjRg3Kly9PQEAAd999t9Nz165dmxdeeIExY8bQvn177r//fkcX34SEBL7++us8VXuWWCJywzyaNm0qBbLjZ5FRZa3Hh41kxvrDEvXSHIl6aY7c+eFvBTt3MbZjxw53h3BdJCcny5tvvik1atQQm80mkZGR8sILL8iOHTsEkFGjRjnKLl26VACZNGlSlvMMHDhQrF+tazZv3ix33HGHBAcHS2BgoHTo0EF+++03p2WdbUuTkJAgZcuWFQ8PD4mJiXFaZs+ePdKlSxcpU6aMABnOFRUVJR06dMhyzOrVq6VXr14SFhYmNptNqlatKg8++KDs37/fUWb58uXSpk0bCQwMlKCgIOnevbts3bpVOnToIFFRURnO52ybM9OnT5e7775bKleuLDabTcLDw6V9+/by448/Zim7adMm6dWrl4SEhIjNZpM6derIO++8I8nJyRnKjRo1SgA5cOBAhu1//vmntG7dWnx8fCQsLEyGDRsmcXFxAsjAgQOzfB6tW7cWf39/ATK8F2flRUTGjRsnt9xyi/j4+EiZMmXk9ttvl99+y/p3IbvjJ02aJIAsXbo0u4+rSOTl9xtYJ9n8XdVxIq5IugLv3QRXLwBwccBCGk+MIzHZanD99bkOVI8ofcvI5qUfubo+rl69SsWKFbn11ludDn5UylU6TuR68vaFuj0dLwN2z6RjrWs9Y3SxKlXUvv76a+Li4njsscfcHYpSgCYR12UeeNjw2ihbHXioisrs2bP55JNPePHFF6lXrx69evVyd0hKAdqw7rpq7SGgHFw8BRdP0dVvDz5eHlxNTmXXiXj2nYqnZrky7o5SlTJPPfUUx44do2nTpowfP750N9SqEkXvRFzl4QkN+jhe+u2eQafa6e5Gtuj08KrwxcTEkJiYyKpVq5yOjFbKXTSJ5EeDdFVaO2Zzd/1ry+bO3XrMDQEppZR7aBLJjyrNIDjKen71PLd7b8bX2/oo95xMYM/JeDcGp5RS148mkfwwxpoGxc5n5090rlPe8Vp7aSmlbhSaRPIrfS+tPQvoVfdaY/rcrcdL3CRqSimVH5pE8qtcXShvTfNM8hU6ylr8vK0eM/tOJbDnZMlb5lIppVylSaQg0t2N2HZOp3Pd9L20tIFdKVX6aRIpiAZ9rz3fv5Teta9Nhz1Hq7SUUjcATSIFERwJVe1LXEoK7RNX4G+zqrT+On2RXSe0l5ZSqnRzKYkYY/yNMV2NMc8YY94yxrxpf97FGONfVEEWa+mqtLx3TOf2utpLS2U1aNAglxc7UgUTHR3tWOiqOIqJicEYw+jRo90dSoHkKYkYY+40xswC4oD5wPvAS8DL9ucLgDhjzExjTLeiCrZYSr/++uE/6Vvj2hKq2kurZFm2bBnGmGwfXl46S5ByTUxMDKNHj2bTpk3uDqXI5PhbYYxph5UkmgExwERgFbAfOAMYIBSoCbQCugHzjDHrgWdE5I8ii7y4CIyAGp1gn7ViWuvLywmw1eNiYgoHYi+y/dgFGlQOcnOQyhUPPvgg3bt3z7Ld2QJQqvjavXu32+/+YmJieP3114mOjuaWW27JsC9tEa6S/uUkt+iXATOB50Tk9xzKrQCmABhjOgBP248t2Z9OXjXo50gi3jt+oku9TszcZPXOmrv1uCaREqZJkyY8/PDD7g4jV0lJSaSkpODr6+vuUPIlPj6eMmWKbrLSgqz7fj0YY0rsv116uX21aiIifXNJIBmIyHIR6Q00LVhoJUim9dfvi7q23vrcLVqlVRqlVX1Nnjw5yz5X2j+OHz/OiBEjiIyMxGazUalSJYYPH86pU6cylBs9ejTGGLZv386zzz5LlSpV8PX15c8//8z23Kmpqfz73/+mffv2VKhQAZvNRmRkJCNGjODMmTMZyqavn//mm29o1KgRvr6+REZGMnr06CzL36a9x9OnTzNgwADCwsIICAigc+fObNy4Mdtzf/fddzRt2hQ/Pz+eeuopR5ktW7bQu3dvwsLC8PX1pV69eowZM4aUlBRHmcOHDxMWFkbDhg25fPlyhms89NBDeHh4ZFhD3VmbSNq2zZs3c/vttxMYGEi5cuV4/vnnSU5O5sqVKzz//PNUrlwZX19f2rdvz86dOzOcIz4+nldffZUWLVoQHh6Oj48PNWvW5OWXX+bSpUuOcpMnT6ZTp04ADB482FEtmhZTdm0iycnJvPPOO9SrVw9fX1/CwsLo3bt3lmWU0x8/Z84cbr31Vnx9falYsSIvvPBCvpcsdlWOdwoisjm/Jy7IsSVOpvXXmyf8ShmfFsRfTebQ2UtsO3qBhlX0bqSkuHTpktO1t202G2XLli206xw6dIhWrVqRmJjIkCFDqFGjBvv27eOzzz5j6dKlrFu3jqCgjP9vHnroIfz8/HjuuecwxlCxYsVsz5+YmMi7775L37596dWrFwEBAaxdu5YJEybwxx9/sH79emw2W4ZjZs+ezYcffuhYx/znn3/m9ddf5+DBg0yaNCnLNbp160ZoaCijR4/mxIkTfPLJJ7Rv355Vq1bRoEGDDGVnzpzJxx9/zIgRI3j88ccdn+W6devo0KED3t7ejuvOnj2bl156ic2bN/P1118DULVqVSZNmkSvXr14+umn+eKLLwCYOHEi06ZN4+WXX+b222/P9XM/cuQIXbp04f7776dfv34sWrSI999/H09PT7Zv387ly5d5+eWXiY2N5b333uOee+5h586djurMo0ePMn78ePr27Uv//v3x8vJi+fLljBkzho0bNzpWnGzfvj2vvPIKb775JsOHD6ddu3YAlC9fPtvYwPo3/v777+nSpQsjRozgxIkTfPrpp7Rq1Yrff/+dxo0bZyg/b948xo4dy+OPP86jjz7KrFmzeO+99wgJCeGVV17J9fMosOzWzU17AOUAW27lSsKjwGus5yTT+uvPfLPBsf76m/NK9hrl2a7BnPZ+i+MjH9LWTc/ucdddd2Upm9c11p1t69mzp0RERMjhw4czbF+7dq14enpmWNM9bd3wDh06SFJSUp7eT2pqqly6dCnL9vHjxwsg3333nWPbgQMHBBAPDw9Zv359hnPcc889AsiqVauyvJ/evXtLamqqY/u6devEGCN33HFHlnN7eXk5/b/UunVr8fT0lM2bN2e47r333iuALF68OEP5p556SgD54YcfZOfOnRIQECAtW7bM8rk4W08+KipKAPn+++8zbG/SpIkYY6Rnz54Z3s9HH30kgCxYsMCx7erVq5KYmJjlfbz66qsCyOrVqx3bcvp/kva5pP93XrRokQBy3333ZYhj8+bN4unpKW3bts1yvL+/f4b15FNTU6V+/fpSoUKFLNd0pqBrrOelpfA44OjHaozxs3ftrV6g7FXa1OwCPvZvjXEx9K962rFLq7RKluHDh/PLL79kefz73/8utGucP3+eOXPm0LNnT3x9fYmNjXU8oqOjqVmzJosWLcpy3NNPP53nhlhjDH5+fgCkpKRw7tw5YmNjue222wBYvXp1lmO6dOlCkyZNMpzjxRdfBGDGjBlZyr/44osZqu6aNm1Kly5dWLx4MQkJGaf+ueuuu7Ks5X3q1ClWrlxJz549adSoUYbrpn2Lznzdd999l8aNGzNs2DD69euHl5cX06ZNy/PnUrlyZe69994M29q2bYuI8NRTT2V4P2l3D3v37nVss9lseHt7A1bVU1xcHLGxsY67IGefa16lvdd//OMfGeJo1KgRPXr04I8//uD06dMZjrnnnnuIjo52vDbG0KlTJ06cOJHl36Ao5CWJZK7c9cfq3htd6NGUZN6+UPdux8vG5xZTxtf6T30k7jJbjpx3V2TKRTfddBO33357lsfNN99caNfYvXs3qampTJgwgYiIiCyP3bt3c/LkySzH1apVy6XrfP/997Ro0QI/Pz9CQkKIiIigenXr+19cXFyW8pn/yAPUq1cPgL/++ivP5VNSUjh48GCusR84cADA6UJb9erVw8PDI8t1fXx8+Oabb4iPj2f79u2MHTuWatWqZTk+O87KhoSEON2Xtj1zG9LYsWNp1KgRPj4+hIaGEhER4WjrcPa55tWBAwfw8PBw+rmmVQ+mfWZp0v490wsLC3Mad1HIb+8pHTXlTMN+sOkrADx3zKBbvQf4YYM14HDu1uPcXDXYndEVvtE3bmLMqeE8Lw2aaXemDz/8MAMHDnRaJu0uIj1//7yP6f3pp5+4//77ad68OR999BFVq1bF19eXlJQUunXrRmpqapZjCqNLbHZ33c5iz+8d+ty5cx2N7hs3bqR///55PjanpYWz25c+zv/85z8899xzdO3alb///e9UqlQJm83G0aNHGTRokNPPNa/y83nk9H6uRw3IjdEF93rJtP76Q+UP8gNWw+XcLccZeWcdt/dbV4UjNDQUgLNnz2bZ5+wbe2Y1a9bEGENiYmKeGoPzY+rUqfj6+rJ06dIMf8B37dqV7TE7duzIdpuzb7w7d+6kZcuWWbZ5enoSFRWVa4xp59y+fXuWfbt27SI1NTXLddevX8/IkSPp0qUL4eHhvP/++3Tp0oWuXbvmer3CMHXqVKKjo5k/f36GsUMLFizIUtbV3/caNWqwcOFCdu7cmaF6D679O7hy13U9uHX0lDEm0BjzijFmqzEm3hgTa4xZaYwZZDJ9+saY2vYR8XHGmIvGmN+NMbe5K3anMq2/3vDsYoL8rLrTo+cus+nwOXdFpgpZtWrV8PLyytClFGDlypU5drtNExYWRvfu3fnpp5+clheRLHXfrvL09MQYk+GbsYjwxhtvZHvML7/8woYNGzKUHzNmDGDVvWc2ZsyYDN92N2zYwOLFi+ncuTOBgYG5xliuXDlat27N7Nmz2bZtW4brvvXWWwD07t3bsT0hIYEHHniA4OBgpk6dyhdffEG1atUYMGBAlm7RRSXtc03/vpOTk3n77bezlE37DJx92XAm7TN+6623Mpx/27Zt/Pzzz7Rt25aIiIiChF/o8non0swYc8X+PG10UFtjjNP6GRH5KbcTGmM8sKZQaY01UPG/WO0tDwKTgLpYbS8YY2oAK4FkYAxwHhgGLDTG3Ckii7NcwF0a3gurPwfAc9ds7qo7kGkbrD8Gc7ccp3FkiDujU3mwYcMGvvrqK6f77rnnHgIDAwkMDGTQoEGMHz+eBx98kI4dO7J3714mTZpEo0aN2Lw59x7un332GW3btqV9+/YMGDCAxo0bk5qayl9//cWsWbMYMGBAgeZV6tevH9OnT+e2225jwIABJCUlMXPmzAxjGTK7+eabue2223jyySepWLEis2bNYvHixTzyyCO0atUqS/mDBw9yxx130LNnT44fP84nn3yCn58f7777bp7j/Oijj+jQoQPt2rVzdPGdM2cOCxcupH///nTu3NlRdsSIEezfv58FCxY4usp+++23tGnThoEDBzJv3rwiv9vv168fI0eO5M4776RPnz5cuHCBadOmORrb06tXrx5lypRh7Nix+Pv7ExwcTLly5RydGzLr0qUL9913H99++y1xcXH06NHD0cXX19eXjz/+uEjfW75k120r7QGkAimZHqk5bE/J7Zz287bC6jr5QabtNuAv4Fy6bd/bz39Lum2BwEFgN2Dycs0i7eKbJjVV5MNGjq6mWxd/5ejq2+rNxZKSkpr7OYqZvHQBLA1y6+ILyN69ex3l4+PjZejQoRIaGip+fn7Spk0bWbFiRZ67+IqInD59Wp5//nm56aabxMfHR4KCgqRBgwby97//XbZv3+4ol9bFN31Xzhu9OLAAACAASURBVLwYN26c1K1bV3x8fKRChQoybNgwOXPmjAAycOBAR7n03U2nTZsmDRs2FJvNJlWqVJHXXnstS5fWtPdz6tQpefjhhx2fQadOnWTdunUZyjrryprZpk2bpFevXhISEiI2m03q1Kkj77zzjiQnJzvKTJkyRQB58cUXsxw/ZswYAeS9995zbMuui2/mbSLZf77OYk9OTpY333xTatSoITabTSIjI+WFF16QHTt2OH2fc+fOlcaNG4uPj4+jm3ZOn0tSUpK8/fbbUqdOHbHZbBISEiK9evWSLVu25Bpbbu/HmYJ28TWSS8OLMcZ5q18ORGRKbmWMMXdgTdz4ooi8m2nfGqCyiFQ2xgRgzdO1QkQ6Zyr3GvB/QAsRWZPbNZs1aybr1q1z4Z3k05J/we/vAZBa9x6a7H6Ic5eSAJg+ojVNo0rW3cjOnTud9hZRpUdMTAzVqlVj1KhRebr7GTRoEFOmTNGu66VAXn6/jTHrRaSZs325VmflJSHk0xrgHPCiMSYGWA34AYOwpkx53F6uEeCDNfFjZmmVybfaz1c8NLzXkUQ89i6gV53HmbLBqhOdu+V4iUsiSimVHbc1rItIHNATOItVXXUQ2AU8CfQVkS/tRSvZfx51cpq0bZWzu44xZrgxZp0xZl1BGyrzrFydDOuv9w+6NufNvK3HSU3Vb29KqdIhxyRijHFtZFPGY2vnoVgCsA14D+gDDAX2AdOMMV3sZdL6Jl51cvyVTGWyEJFxItJMRJpd114N6RaruunUAkL8rUa3ExeusP5Q/gcjKaVUcZLbnch2Y8xEY0yDXMo5GGMaG2OmYiWHnMo1xOpx9YuIvCAiM0RkAtAWOAF8aYzxBNK6kjib1zltHuXsu5u4S7r11z3+WkbfOtemfNYVD1VxEx0djYjkuTfY5MmTtT1EAbknkbuBFsBmY8xGY8wbxphexpgGxphKxpjKxpiGxpg+xpgxxpitwDrgFqBHLud+BisJ/JB+o4hcAuYCUVhTqxyz73JWZZW2zVlVl3tlWn/9wYBrfe/nbT1OilZpKaVKgRyTiIgswGrYfhir6ukV4CdgM3AYOARsAn4EngdOAw8AN4vIwlyunZYAnI3Z90r3cytWVVbWTuqQNlT2OnS5yod0VVrVT8wjLMAavX4q/irrYvI2+EgppYqzXBvWRSRFRL4RkXZARazBgP8CPgc+w+piey9QTkRuE5EfRCQvk8ekza8wKP1G+wDGXljrue8XkQRgNtDRGHNzunKBWG0oeylOPbPSq9/bsf66ObyaB9K1MM3dWrKqtLTqQqnSpzB+r12aO0tETmL1pCoMHwIDgLft7SMrsNZrH4aVrJ4UkbSZ7EYCnYFFxpgPgAv2cpWBu6S4/oULCM+w/vp9vmv4FGud5XlbTzDq7vp4ehT/ubS8vLxITk52OiJXKVVyJSUl5TiBY164s4vvQaA5MBXohDXtyctY1WR9RWRsurL7gDZY40JexurNdRHolodqM/dqcK1KK/LYXMIDrf4BsQlXWXOgZFRp+fr6Xpd1CZRS19eFCxcKvM69y0nEGFPGGPNPY8wfxpi9xphW9u3h9u118nouEdkvIgNFpIqIeItIWRFpL07m3hKRnSLSS0SCRcRfRNpKcZozKzvp1l83J7czsOa1jmRztx7L7qhiJSIigtOnT3Pp0iWt1lKqhBMREhMTiY2NJS4uzjEjdX65VJ1ljIkA/gCqY43nqI41yhwRibVPkRIMPFugqEoT37JQqxvsmAlAH+8/eZ82ACzYdoLRd9fHy9OtkynnytfXl/Lly3PixAmuXnU2XEcpVZJ4enpSpkwZIiMj8fFxNnoi71xdT+QNoAJWt99DQOa5l2dhtV2o9Br2cySRSkfmUi6wE6cSEolNSGTNgbO0rhnu5gBzFxQURFBQkLvDUEoVM65+Be4BjBWRDVizmmb2F1C1wFGVNunWXzdxMQyrcW3E+pwS1ktLKaXSczWJhGNVY2UnlWujyFWaTOuv9/RY4Xi+YNsJklPyv5ymUkq5k6tJ5ARQI4f9jbGquVRm6QYeljs0j4plrJrEsxcT+fOvktFLSymlMnM1icwDhhhjKmbeYYxpgTXuY1ZhBFbqpK2/DpiLpxgRfa0aq6T00lJKqcxcTSKvYy1RuxF4C6tdZKAx5hvgN6x5rt4p1AhLi0zrr3fnD8fzBdtOkKRVWkqpEsilJCIiJ7Dmq1oNPAoY4BHgPmAR0E5EtG4mOw3vdTwNO7SQqLLWxx93KYmV+8+4KyqllMo3lwcoiMhhEemFNUVJC6ykEiEid4vIkcIOsFSp3BRCogEwVy/wZNUDjl1zt2iVllKq5MlzEjHGBNrXFrkXQEQuiMhaEVmjdx95ZEyGaVC6pFyr0lq4/SSJyVqlpZQqWfKcROyz6T4AlC26cG4A6aq0go8s4aYga7jN+ctJrNgf666olFIqX1ytztqBtVCUyq9066+b5Cs8VXm3Y5eueKiUKmlcTSJjgBEFWXtdkWHMSMfE5Y7nC7ef0CotpVSJ4urcWXWwpmrfaoyZg7UgVOb1zUVE/lUYwZVaDfrC4tEAlDn6Bw2CB7PtnI34K8n8se80t9Up7974lFIqj1xNIqPTPe+dTRnBWvlQZSdt/fXDf2IkhacqbOOxc00AmLPluCYRpVSJ4WoSqVYkUdyIGvaDw38C0PbyMsBKIr9sP8nV5BR8vAq22phSSl0Prg42PJiXR1EFW6qkW3894OQ6bg22Vg6Mv5rM73u0l5ZSqmQo3qshlWZp66/bPRmx2fF8rk4Pr5QqIVxd2XBiHoqJiAzJZzw3lob3wj5rhd8WF38F+4qHv+w4yZWkFHy9tUpLKVW8udomMigPZQTQJJIXaeuvJ1/B7+xOOobEsiwunISryfy25zRd61dwd4RKKZUjV9tEPDI/AG+gNvAl8CcQUgRxlk4+Zaz11+0eC93oeK5VWkqpkqDAbSIikiIie0XkMeAMOhW8a9INPGx6YQlpqw4vtldpKaVUcVbYDevzgb6FfM7SLd3667b4Q3QPOQrAxcQUlu0+5c7IlFIqV4WdRMKAwEI+Z+nm7Qv1rq2/PjR4veP5HJ1LSylVzBVKEjHGBBtj+gHPAOtzK68ySTc9fKNzv+KBNX/Wkp2nuJyoVVpKqeLLpSRijEk1xqRkfmC1hXwPpALPFkWgpVq69de9Lp+mb8h+AC4npbBUq7SUUsWYq118/0day+81ApwF9gDfiEh8YQR2Q/HwtCZlXP0ZAAPLrOOHuJsAa3r47g0rujM6pZTKlktJREQGFVEcqmE/RxKpe245PvTlKjaW7DrJpcRk/G2u5nullCp6rlZn/dMY0yCH/fWNMf8seFg3oHTrr3smXqB/qLVY1ZWkVH7dpVVaSqniydWG9dFAoxz2NwBG5TuaG1mm9df7+61xPNcVD5VSxVVhd/H1BZIL+Zw3jnTrr9c49weB9vW+ft11iotX9WNVShU/uSYRY0xZY0ykMSbSviks7XWmxy3AQ1grH6r8KFcHyjcEwCPlKgNDtwNwNTmVr1frDPtKqeInL3cizwAH7A8BPkz3Ov1jPXA78HmRRHqjaHhtwP+D6aq0PvhlL0fPXXZHREopla28dPlZZv9pgH8CM4AtmcoIkAD8KSIrCy26G1G69dcrn/2T5hGPsea0F5eTUhj983a+HNDMvfEppVQ6uSYREVkOLAcwxkQBn4vI6qIO7IaVaf31d+sdoMNya8zILztOsmj7CZ0iXilVbLg6FfzgwkwgxpjRxhjJ4ZGUx7LPF1ZMxUK6mX2jjs3ngVurOl6P/nm7NrIrpYqNfI1gM8Z4AnWw1g7JkohE5Lc8nuonYJ+T7Y2AF4DZTvY9A2RehLx0zddVvzfMfwkkBQ7/ycjuASzaYePsxUSOnb/Ch4v38I+76rk7SqWUcj2JGGNeAl4GyuZQLE/ruorIFrK2r2CM+cL+dIKTw2aKSExezl9ipa2/bl86N2j/z/yj+70894O1DvvEFTH0blyFepVy+idQSqmi5+qI9aHAW8Am4FWsxvYPgXex5s9aBzxakICMMf7AA8BRYEE2ZcoaY0r3PCDpxoyw9Uf6NKlMy+qhAKSkCv+YuZXU1MzTmCml1PXl6mDDx7F6YHUCxtm3zRWRl7GqoKLJ411IDu7DusuZJCLO5kHfApwHrhhjVhpj7izg9YqntPXXAU5uwxzfxBv3NMTb0wCw8dA5pq055MYAlVLK9SRSF/jB/jzta7AXgIgcx0os/6+AMQ2xn3tipu3n7Od/CugFjASigLnGmEHZncwYM9wYs84Ys+706dMFDO068ikDtdPlx1l/o2aIFyM61HBsemfBLk7HX3VDcEopZXE1iaQAF+3P036GptsfA9yU32CMMbWBtsCvInIg/T4R+VBEHhORKSLys4i8i3X3cxL4wBjjdEVFERknIs1EpFlERER+Q3OPjiMz3I2w5P94olNNosL8AYi/kswbc3e4MUCl1I3O1SRyCKgGICJXsaY4aZdu/61YbSP5NcT+c3xeCovIGawR8sFA6wJct3iKqA1d37j2+s9P8T24jDfuuTaR8qxNx/hjb+bOakopdX24mkR+A+5K9/oH4DFjzERjzGRgKDAvP4HYG8oHYCWhGS4cGmP/GZ6f6xZ7tw6Fm+649nrmCNpV8qDnzZUcm16btY0rSbqMrlLq+nM1iXwEfGqM8bO/HoWVNAYCjwC/YHX/zY+7gfLAVPtdTl6lVZ+dzOd1izdjoNenEGCviks4CT8/xat31aGMr9VB7UDsRT5btt+NQSqlblSujljfLSJfiMhl++uLItITq10kSETuFJH8VmelVWVlGRtijPEyxgQ52V4VGIG1xnvpnbMrMAJ6jb32evdcyu39jhe71XFs+mzZfvafTnBDcEqpG1mek4gxJtBebXVv5n0icl5E8v0XzBhTCegGrBGRrU6KBAIxxphJxpgXjTHDjDHvAVuBIGBEWmIrtWp1hebDr71eMJL+1a9yc9VgABJTUnlt5jZEdOyIUur6yXMSsSeJB8h5pHp+DcIaX5Jdg/plYDrQDKtr71istUsWA61F5IdsjitduvwfRNjvPpIu4TljGG/1rIWHNXSElfvPMHPTUffFp5S64bjaJrIDa0BhoRKRN0XEiMiX2ey/KiJDRaShiISIiLeIVBSRfiKyxtkxpZK3H/QdD5426/XxTdTb/SmD21RzFHljzk7OXUp0U4BKqRuNq0lkDDDCGFOrKIJReVChIXROt4z9Hx/yfK1TVAyyxpOcuZjIOwt2uyk4pdSNxtX5p+pgjQ3ZaoyZA+wF+0Lg14iI/KswglPZaPkE7PsF/loGCH5znuRfd/zE0O+tCZG/WXOIfk0r0zQqNMfTKKVUQRlXGmKNMal5KCYiUtD5s4pEs2bNZN26de4Oo3BcOA6ftYLLcQBI/d4Mu/gEi3dZU7vUqVCG2U+1xdvT1ZtNpZTKyBizXkScLqvq6l+Yanl4VM9/qCrPylaEnv91vDTbZzCm1g78vK38vetEPBP/OJDd0UopVShcHSdyMC+PogpWZVL3bmgy0PEydNk/+GcbX8frDxfv5Uhc5tpGpZQqPPmu6zDG1DTGtHE2CFBdR93eglD7zL6JCdx/5A3qlbcmaLyclMLon7fr2BGlVJFxOYkYY3oYY/YDu7Hm0mpq317OGLPPGNMvxxOowmULgL5fgofVR8LjyFrGV1vm2L145ykWbi+dM8IopdzP1ZUNO2JNjngWeB1rZUMAROQUsB9rQKK6nio3hU6vOF5W2vwxL9U/73j9+uztJFxNdkdkSqlSztU7kX8Cm4EWwKdO9q8CmhQ0KJUPbZ6GqDbWc0lleOzbRPpbieP4+St88MseNwanlCqtXE0izYCvRSS7rr5HgAoFC0nli4cn9P4CfKwmKs/zB/mq8o+O3ZNWHGDb0fPZHa2UUvniahLxBHKapj0c0Dk33CW4Ktz9geNl5OGfea6iNZ9lqsA/Zm4jJVUb2ZVShcfVJLKTjCsZZtYDq7pLuUuDvtDoWrPUExc/JcrTmp1/8+FzTFtzyF2RKaVKIVeTyASgnzFmSLpjxRjjb4z5GGgFjCvMAFU+dH8XgqMA8Ey8wFdhE/HAqoEcs2AXp+KvuDM6pVQp4upgw8+A74AvsebNEuAb4DzwN2CyiHxd2EEqF/mWhT5fgrH+eate2MDLZRcCEH8lmTfm7HRndEqpUsTlcSIi8jDQF1gC7MLq7jsPuFdEhuR0rLqOIltA+xcdL4cmfUND8xcAP28+xm97TrsrMqVUKZKvEesiMkNE+opIfRGpJyK9RGR6YQenCqj9C1DlVgA8JJkJgZ/jh1WV9dqsbVxJSnFndEqpUqBAU7waY/yMMX6FFYwqZJ5eVrWWLRCAcklH+JevVdt48Mwlxi7d587olFKlQH6mPSlnjBlrjDkGJAAJxpjj9m3lCz9EVSCh1ayGdrt+LOEOj7UAfLZ8P/tOJbgrMqVUKeDqtCfVgI3A41iN6bOAn4Fz9m0bjDE6FXxxc/ODUL+34+W7PhMoRxxJKcKrM7fqBI1KqXxz9U7kfSAM6CMidUWkj4j0FpG6WI3tYcB7hR2kKiBjoMcHULYyAGXlAu/bPseQyp9/neWnDUfdHKBSqqRyNYl0Bj4VkZmZd4jIDOAzexlV3PiFWNOi2OfMbOexlUc9FwDw73k7OXdJJxpQSrnO1SQiWONDsrPHXkYVR9XaQdunHS9f8v6WuuYgZy8m8vb8XW4MTClVUrmaRJYDnXLY3xFYlt9g1HXQ8RWoeDMANpL5yPsTfEjk27WHWRtz1s3BKaVKGleTyNNAC2PM+8aYcmkb7T22/oM1RfzT2R6t3M/LBn0ngJfVM7uWx1FGek0D4B8ztpKUkt0EzUoplZWrSWQJ4IeVKI4bY84YY2KB48D/s+/71RjzV7rH/sINWRVY+E3Wsrp2g7wW0dFjI3tOJjD+9wNuDEwpVdJ4uVj+ENrmUTo0HQR7f4HdcwF41/sLul19h4+W7KFHo4pUDfV3b3xKqRLB3EhjBJo1aybr1q1zdxjFx8Uz8FkrSLDWYF+S0pghSc9zW53yTBjYDGNMLidQSt0IjDHrRaSZs30FmvZElXABYXDPZ46XnT038rDnYn7ddYqF20+4MTClVEmR7yRiX0OkqjEmMvOjMANURaxmZ2j5hOPlq15fUdMcYfTPO0i4muzGwJRSJYGr0554GmNeMcYcBeKBGOCAk4cqSTqPgnL1AfA1SXzk/SlnL8Tzn0V73ByYUqq4c7Vh/T/AU8AG4AcgrtAjUtefty/0HQ/jOkLKVep7HOQ5r+95Z+VD9GlSmQaVg9wdoVKqmHI1iTwE/CQi/YoiGOVG5etB13/BfGshq8e85vJbaiMGT/Zh0qBbNZEopZxytU3EG1hUFIGoYqD5cKh5u+Pl+96fkxQfy31frGLp7lNuDEwpVVy5mkRWAvWKIhBVDBgDvcaCfxgAFUwcE2zv4ZMYx9Ap6/hmzSE3B6iUKm5cTSIvAv2NMb2KIhhVDJQpD70+dbxs6rGXGbZRRMoxRv60lfcW7tb1R5RSDi4lERHZCgwDphtjDhljlhtjfs30WFI0oarrpvad0O1t0qaNj/Y4yU+2UdxqdvHJ0n08+/1mEpN1ji2llOtdfLsD39uPKwtEAtUyPfK8sqExZrQxRnJ4JGUqX9sYM9MYE2eMuWiM+d0Yc5sr70HlUcsRcP9XjokaQ0wCX9nepJfHH8zYeJSBE9dw/nJSLidRSpV2Lk17YozZAvgDve13JQW7uDGNgEZOdjUCXgBmiEgfe9kawBogGfgQa3neYUAD4E4RWZzb9XTak3w4uh6mPQAXrzWsv5/Uj/+m9KZW+TJMGtycysF+bgxQKVXUcpr2xNUkchl4SUQ+LqzgsrnOF8BwoIeIzLVv+x5rCd6mIrLJvi0Q2A5cAepILm9Gk0g+xR2EaffB6WsLV/2Y0p6RSUMJKRPApMG3Ur+SdgFWqrQqzLmzDgK+BQ8pe8YYf+AB4CiwwL4tAOgJLEtLIAAikgCMB2oBtxZlXDe0kCh4dCFU7+jY1M/zN6Z4v82V+DPc9/kqlu857bbwlFLu42oS+RgYar8DKCr3YbW3TBKRFPu2RoAPsMpJ+T/tPzWJFCW/YHjoR2j8sGNTa88d/GQbTUjScR6dvJbv1moXYKVuNK6OWE8AzgE7jTGTsObJSslcSET+V4CYhmCtWTIx3bZK9p9HnZRP21bZ2cmMMcOxqsaIjNS5IQvE0xt6fgIh1eDXfwFQ0+MYM23/ZGji87w0XTgad5lnutTSaeSVukG4mkQmp3v+ajZlBMhXEjHG1AbaAktEJP1EjmkrJF11ctiVTGUyBiMyDhgHVptIfuJS6RgD7Z+HkGiYOQJSEgk3F/jW9i+eTnqSj3+FI+cu83afRti8dKUBpUo7V5NIpyKJ4poh9p/jM22/ZP/p4+QY30xl1PXQsB+UrQzf9ofLZ/E1SYz1/og3k/szfkN3Tl24ytiHm1DW19vdkSqlipBLSURElhdVIMYYL2AAcBaYkWn3MftPZ1VWaducVXWpohTVCoYuhq/vhbP78TDCq95fE21OMGrfIO77fBWTBt9KxSDtAqxUaVWQRal8jDGVjTG2QorlbqA8MFVEMldbbcWqymrl5LiW9p/ad9cdwmpYiSTy2j/Nw15LGO/9HodPnKL3pyvZceyCGwNUShUll5OIMaaJMeZXrEWpDmG1YWCMKWeMWWKMuT3HE2QvrSprQuYd9q68s4GOxpib08USCAwF9mINRFTu4B8Kj8yEBtdWCOjkuZkfba/DhaPc98Uqft+rXYCVKo1cnfbkFuB3oAaZGs9F5BTgBwx0NQhjTCWgG7Amh5HwI7FGqS8yxrxsjHnCHktl4KncBhqqIpa2sFX7Fxyb6nocYqbPP4lK3MfgSWv5Yd1hNwaolCoKrt6J/B9W+0R94GXSZui7ZgnQPB9xDAI8ydqg7iAi+4A2WONCXgbeAy4C3URkYT6uqQqbMXDbq9YswB5Wc1sFE8f3ttdpxwZe+HELHy7eo7MAK1WKuJpE2gFf2quXnP0lOMS1MR15JiJviogRkS9zKbdTRHqJSLCI+ItI27zMmaWus8YPw8PTwceaCiXAXGW893s84rmIDxfv5cUft5CUorMAK1UauJpEfLGqlLJTtgCxqNKkekcYsgiCrAGenkb4l/dkXvWayvT1h3h08lrir+gswEqVdK4mkf1A0xz23wbsyH84qlQpVweGLYFKTRybhnrN5zPvD1m39wj3ffEnJ85fyeEESqniztUkMg14JFMPLAEwxjyH1Tg+tZBiU6VBYDkYNBfq9HBsusNzHd/a3iD2+CF6j13BrhPaBVipksrVJPIeVsP2QuA3rATygTHmKDAG+AUYW6gRqpLP5g/3/Q9a/c2x6WaPv5jh808CL+zj3s9WsWJfrBsDVErll6vL4yYCXYDngctY81bVAmKx1l/vISLaYqqy8vCEO/4Nd70PxvpvV8XEMt02moZJmxg4cQ3T1x9xc5BKKVe5Ok5kItaiUB+ISDMRCbD3krpZRN4HmtjLKOXcrUPhwe/AZq0mUNZcYor3O/Q2S3nuh838d8le7QKsVAnianXWIKyBhtmpRj4GG6obTK2uMHg+lLF6g3ubFN71HsfzXt/xn192MfKnrdoFWKkSorDn6g4AtN+myl3FRtacW+UbOjb9zWsWH3l/yoy1+3noy9WsjTnrxgCVUnmR6yy+xphIIDrdpjrGmPZOioYCI4B9hROaKvWCKsOj8+GHwbDvFwB6eq6iojnD8JhnuffzszSPDuWJTjXoUCtCF7pSqhgyudU/G2NGAaNwPkI9Q1EgFRgsIsWym2+zZs1k3Tqd7LfYSUmG+S/Cumtzb56UYMYnd+eblNtIwJ/6lcryRMeadGtQAU8PTSZKXU/GmPUi0szpvjwkkZuBW7CSxESsVQIzr3UuWEvnrhWRYjvLniaRYkwEVn0Ki14l/feVC+LHtJTbmZjcjVOEUD08gMc71uCeWyrryolKXScFSiKZTjQKmC4i2woruOtJk0gJsHMOzH0WEk5m2JwonsxKacO4lB7slSpUCvJlWPvqPHBrJH42TzcFq9SNodCSSEmnSaSESL4KW76DFR/Dmb1Zdi9Oacy45B6skTqEBfjwaNtqPNwyiiA/XYpXqaKgScROk0gJk5oKexbAio/g8J9Zdm9MrckXyT1YlNqMAB8bj7SK4tG21QgP9HFDsEqVXppE7DSJlGCHVsPKj2HXXDL38TiQWp7xKXfxY0p78PLlgVurMqx9daqE+LsnVqVKGU0idppESoHYvbDyv7D5W0i5mnGXlGVKclempnQhwaMs9zSuzOMdalCzXKCbglWqdNAkYqdJpBSJPwlrvoC14+FKxiVuLokP36V0ZEJKd44SQbf6FXiiY00aVglyU7BKlWyaROw0iZRCV+Nhw1Sre/CFjBM4JosH81Jb8EVyD7ZLNdrXiuDJjjVoXi1UBy4q5QJNInaaREqxlCTYPsNqhD+ZtQf6Hyn1GZfSg99SG9E0KpQnO9WgU+1ymkyUygNNInaaRG4AIrD/VyuZHFieZfeO1Ci+SL6LuaktualiKE90rEH3hhV1FLxSOdAkYqdJ5AZzbJPVo2v7DMi0zM0RCWdi8p18m9KJcmGhPN6hBr2bVMbHSwcuKpWZJhE7TSI3qLgYWDUWNk6FpEsZdp0Xf6amdGFK8h14lq3AkLbV6FKvPFFh/lrVpZSdJhE7TSI3uEtnrd5cq7+ASxmX470qXkxPaceElO7sl8pUCvKlVY1wWtcIo3XNMCoG+bkpaKXcT5OInSYRBUDSZdg0DVZ9Amf/yrJ7f2pFVqXWY0VqA1al1uMcZagWHmAllBrhtKweSpiOilc3EE0idppEVAapKbBrjtUIf3S98yJi2CFRrEitz8rUBqxJrc1lfKlToQyt7XcqzauHPkdv/AAAEo1JREFUUtZX5+1SpZcmETtNIsopETi40hprsn8JJF/JtmiieLJRbmJlSn1WpNZns9QkxXjRsEowbex3Kk2jQnRmYVWqaBKx0ySicpV0BY6sgb+Ww4HfrDsUScm2+EXxYW1qHcedyg6JxNvTi8aRwdadSs0wbq4SrGufqBJNk4idJhHlsisXrLuUA8utxHJqe47Fz0ogq1LrsTK1AStS6xMjFfC3eXFrdKijTaVepbI6LkWVKJpE7DSJqAJLOGXdoaQllXMHcyx+VMJYmVKflan1WZHagFOEUNbXi5bVw2hdI4w2NcOpWS5QuxOrYk2TiJ0mEVXo4mKspPLXciuxXDydY/F9qZUcVV+rUutygUDCA31oVSOMehXLEh3mT1RYAFFh/gT4eF2f96BULjSJ2GkSUUVKBE7tvHaXcnAFXL2QbfEUMWyTao67lN2pVYglCMFqP4ko43MtqYT6ExUe4Hitqziq60mTiJ0mEXVdpSTDsY1wYJmVVA6vybIGSmZXxZujEsZRCeeIRHBUwh2PIxLBSUJIwZMQf2+iwgIy3LmkvQ4NsGn1mCpUmkTsNIkot0q6DIdXX6v6OrYxy5xeuUkWD04Q6kgqRzIlmeMSho+PH1Hh15JKVKiVZKLDAyhXxkcTjHKZJhE7TSKqWLl8zqry+mu5tYZ83EG4cq5Ap0wVw2mCstzJHLG/PutVnnJhoVZSCQsgMsyfysF+hAbYCPG3ERJgI8DmqYlGZaBJxE6TiCr2rsbDucNw7hCcz/Tz3GG4eKrAlzgrgRkSzEkJ4TwBnJcALhDAJROA+IXg5R+MzT+IkEBfgv1thAZ4W4nG30ZIhuc2yvp6aeIpxXJKIm7v/mGMCQVeAe4BqgDxwDbgnyLyu73MaGBUNqd4QUTeuw6hKlX0fMpA+XrWw5mky3D+qNW1+PzhTAnnMBJ/7P+3d+9BkpXlHce/v3O6e2Znlr0BkYvAKreAYcUgBOMNjLECJGqRkCoSDEiQcJMSE6FUogiIVZssUIGosAJrVqCMBWyglKu6EbkEgSBoFBaXXQQWdhf2PtMz032e/PG+zfa2PbMzPdN9dmeeT9Wp0/2e93Q/Zy799Dnve94XbecS2RxtZo42M48Xhq9UATZCdYPYRA8brWebRPOa9bKMnjefb1YvldJM6J5F2jOLQu9sunrnMGOXXmb3FJkdz3Tm9BZDQuopMWNa0e+XmQRyTSKS9gOWAtOBG4DngJnAPGDvJrtcAKxtKGs+6JFzk1FxGux2QFiaUHUINr4ckss2SSacydiGl1A2NOq3S2XMYguztAUYufsyGdAXl/hf2mddbKD3zST0hvWygh42xuRTSbqxtItq2g2Fbih2Y4VukuI0VJpGUuohLXVT6OqhUOqh2N1DqbuHrq4uekopPaWU7mJKT6nAtGJ4Pq2UMq0YlsSTVNvlfSbynRjDPDNbNYr6S8xsRXtDcm4nlhZh9tywNKEsg82v1iWZlbBlbWifKW+A8nqy/vVY/3oorycd2jKucHo0QA8D7Kk3hq9khDOfCjD8sGXbGLKUMiXKFBmgxIAVKVNiNSXKVnpz25C6qCRdZGkX1bSLLJ2Gpd1YoYSSAqQFSIqQFlFSQGkxLuFxUiiRpAWS2uNikTQtkhSKpMUSaVokLZQolEqkhQKFQheFYolisUAxTeIiSmlCIU0opKKQiEQiTUQq7fSJLrckIukDwPuA881slaQiUDSzvu3sNwPoM7NKJ+J0blJJEpixV1j4o+ZV6p9UK28mF8rrY7KJCafucbVvHdUt66j2r0flDSQD6ykMbiJhbL3PRquoKkX62YX+ULC9z+FqXDokM1EhoUKBCumbj8skZCRULSFDZPF5RoIpPDcSMiUYIiON5WFtJJgSTGnYrhRiOXXlpjQ+31qOEvb60/M46LCjJvRY8zwTOT6uX5R0F3AckEpaBlxqZt9pss/TwC5AVdJjwGVmdndnwnVuCkoL0LtrWEaqFpdtmIWOAk2SDv3hjKcy2E91oJ/qYD/ZYB/ZUD821B8GwqyUUaWMqmXS6gBptUyaDVDIBkjYsTsEJTJKVCkNl7nGcvJhDetx+PnavwAmTxI5OK4XAsuAU4Eu4LPAYklFM7sp1lkPXA88DKyL+34G+L6k081s0XBvIulM4EyAfffdtw2H4ZxrSoLuGWFpIgFKrbyuGVSHoFJLNluTDpVy6HxQKZMN9lMZ6GNwoI/KQB+VgX6qg31UB/qxyiBWHcKyIaxaCa+XVbBsCFUrkA2hrArZEIlVUBYXq5BkFRKrkFiVxCqk8XFKeJxS3WGTXJJM/BQFuXXxlfQA8CfAcuAQMxuM5bNjWRnY26x5VxNJuxJ6cXUD+5jZ5u29p3fxdc51RJZBNhST01CYAK322LKwZGGdVStUsypZtUqWValmFawayrMsI8tCucXtWbWK2dbtllWwLAvlWRWzKpZl4XFtbWH/PY46kV333n/Mh7OjdvGNFzO5tZZAAMxsnaQ7gb8jnHH8qtnOZva6pG8ClwB/DNzX3nCdc26UkgSSLihsfxrlhIZ2qJ1MnrG/FNevNtlW66k1ezuvsSKud5uIgJxzzo1Nnknksbh+a5NttbLt3Z57YFy/NiEROeecG5M8k8gSwt3pp0iaXiuUtCfh7vVlZva8pIKkmY07S9oHOBt4ndDg7pxzrsNyaxOJbR//BFwHPCrpRkJnjbPj+rxYdTrwgqQlhPaRWu+sM+K2k82sv/H1nXPOtV+ud6yb2fWS1gIXApcRBk54BPgbM3soVusHbiPcGfVxQuJYCzwAzDezx37nhZ1zznVE3sOeYGa3A7ePsH2AcNbhnHNuB7Mz9yxzzjmXsyk1n4ikNcDKFnffjd8dQXiy82OeGvyYp4bxHPN+ZrZ7sw1TKomMh6THh7tjc7LyY54a/JinhnYds1/Ocs451zJPIs4551rmSWT0rs87gBz4MU8NfsxTQ1uO2dtEnHPOtczPRJxzzrXMk4hzzrmWeRJxzjnXMk8iI5CUSLpA0q8llSX9VtICSb15x9YOkg6SdKmkRyWtkbRJ0lOSvjhZj7mRpB5JL0gySdfmHU87SZoj6V8lPR//vtdI+rGk9+cd20STNF3SFyQ9E/+u10p6WNJpksYy4/kOR9LnJX1P0vL4d7tiO/UPlrRE0jpJWyQ9KOlDrb5/7mNn7eCuAs4H7gAWAIfE5++S9OHhpu7diZ0OnAvcCdwMDAHHApcDfy3p6CkwYvKlTIFJziTtBywlDGh6A/AcMBOYB+ydX2QTT1IC3E2YAfXbwDVAD3AycBPh//qi3AIcvyuAN4AngVkjVZS0P2HqjAowH9gAfAq4V9JxZvbAmN/dzHxpsgDvIIwqfFtD+acBI4w0nHucE3zM7wZmNim/PB7zeXnH2Obj/8P4z/XZeLzX5h1TG4/1QeC3wJ55x9KBY31P/H1e1VBeApYD6/OOcZzH9/a6x78AVoxQ9z+BKnB4Xdl0wnBQzxJ77I5l8ctZwzsZEHB1Q/lCoA84peMRtZmZPW5mG5ps+m5c/0En4+kkSSnhd3sPI4wqPRlI+gDwPsJUCqskFSX15B1XG82I61fqC81skDCW1JaORzSBzGz5aOrFS9IfBZaa2VN1+28GvgUcBBw51vf3JDK8IwlnItvMV2JmZeApWvhh78Rq0xVP5mmILwB+n62ToU1mx8f1i5LuIszZs0XSc5Im3Zcjwv/weuBCSSdJ2je2C3wNOAK4JNfoOmce0EWYs6nRo3HtSWQC7QWstTCfSaOXgd0klTocU8fFb+hfIlzmuSXncNpC0tuArwCXmtmKnMPphIPjeiEwBzgV+HtgEFgs6ZN5BdYOZraO8A38DcLlnJXArwntf39pZgtzDK+T9orrl5tsq5WNuT3MG9aH1wM0SyAA5bo6g50JJzdXA0cDXzCzZ/MOpk2+AbwAXJl3IB2yS1xvAo6Nl3WQdAehjeAKSd+2ydVxZDOhveBOQsPyHEISuUXSx8zs/jyD65DaJctmn2vlhjqj5mciw+sjnPo1011XZ9KSdBnh8s71Zva1vONph3j55iPAWWY2lHc8HVLrYXdrLYHAm9/Y7wT2YOvZyk5P0mGExHG/mX3OzO4wsxsI7UKvAgvjGfdkV/u8ava51vJnmieR4b1CuGTV7Ae+N+FS16Q9C5F0CXAxoQvkWflG0x7xd3sl8APgVUkHSDoA2C9WmRnLRuw2uRN6Ka5fbbJtVVzP7lAsnXAB4UPye/WFZtYHfJ/w+57b+bA6rtaxoNklq1pZs0tdI/IkMryfEX4+R9UXSuoGDgcezyOoTpD0ZeDLwH8AZ1jsBzgJTQN2B04AltUtS+P2U+LzM/IIro1qnUXe2mRbrWx1h2LphNoHZLOzjULDejJ7hnAp6z1Nth0d12P+XPMkMrzvEvqWf6ah/FOE64Y3dzyiDpD0JUJvlcXAJyfZdfFGW4CTmiznxO33xOd35hJd+ywhtIecIml6rVDSnsDHgWVm9nxewbXB/8X1afWF8QzzY8A64DcdjqnjYlfeu4BjJL2zVh7/Bs4gfGF6bJjdh+VDwY9A0jWENoE7CJc8anesPwR8aLJ9wEo6F7gWeBH4Z0IX53qvTYUGSElzCQ3t/25mk7LLr6QzgeuAXwI3Em68OxvYE/hzM7svx/AmVLw7/0nCJbqbCf+/cwhfCOcC55rZ13MLcJwkfYKtl2A/TfhdLojPV5rZ4rq6BxASxRBhRI6NhJ/DYcAJZnbvmAPI+27LHXkhnP7+I+FOzgHC9cIrgel5x9am411EOPsablmad4wd+jnMZZLfsR6P80TC/QFbCGcm9wHvzTuuNh3r/oQhT16KH6AbgZ8AJ+Yd2wQc29Kx/M8Svgz/F+HemT7gp8CHW31/PxNxzjnXMm8Tcc451zJPIs4551rmScQ551zLPIk455xrmScR55xzLfMk4pxzrmWeRJxzzrXMk4hzoyDpGEkm6bS8YwGQdLakjZJ2HeN+h0vKJH2wXbG5qcWTiHN14ofsJXHokx2SpJmESbSuMrPXx7KvhWlRlwALJKkd8bmpxZOIc9s6nDCC8dyG8p8QRv1d3LhDDs4BZhHGOWvF1YRpYY/fXkXntseTiHOjYGaZmZXNrJpnHJIS4EzgbjNb0+LLPAisYJLOE+M6y5OIc1GciOum+PTHsQ3EJC1q1iZSXybpHEnPSipLekbSCbHOYZLuie0Xr0v6N0nFJu99oKTFklZJGpS0QtK/SOptqHoU4SzpB01eYx9JN0paKWlA0mpJD0s6tb6ehQHz7gX+rH4oeOdaMRUmYnFutG4nDIV+JnAF8KtY/huGnyoZwlzds4FvEeaqPh9YIukkYCFwK6Ed4iOEobpXA5fXdpZ0BPAjwqiq1xFGi35nfJ33SvqgbZ26t9Ygvs28D5IKwP2ECZi+DjwHzATmAe8njGBb7xHgHwhTxN4zwrE5NyJPIs5FZva0pEcISeR+M1ta2ybpmBF23Qs41Mw2xLo/An5OSEp/ZWa3x3rflPQEIelcXrf/jYRpaY80s0117/nD+Bp/SximH+DQuG6cROlQwrzoF5nZ/FEcbm3/d+BJxI2DX85ybvwW1RIIhGREmK/ilboEUvNTYI/aZSRJhxHOFm4BuiTtVlti3S2EM5ia3YGKmW1seN3a+x8r6fdGEXOtV9do6jo3LE8izo3f8iZl6wizIzYrB6jd33FIXH8FWNOwrAZ6gbfU7W+AGrvnmtlK4KuEhLNK0hOS5ks6cpiYa/v7hEJuXPxylnPjN1yPrZF6cqlhvYDhLyutq3u8hjDj5gy2nn0AYGYXS7oROIHQDnIG8DlJ883soobXnFP3es61zJOIc9vq9DfzZXFdNbMHRlH/F3F9IPB440YzWw5cA1wjqZvQC+tCSQvMbHVd1QMaXs+5lvjlLOe2tTmu54xYa+L8L+GD/CxJb2/cKKkgqT6WpXF9dEO9mY1dh82szNYeZrMbXvpooAI81HrozvmZiHONfgZkwBclzSY0bDdr25gQZmaSPkHo4vt0vBz1S6CHcLZwIvB5tvbOeoLQBnM8296xfixwvaTbgGcJyfAIwiWt/zGzZ2sVY3vKccA9ZrYZ58bBk4hzdczsRUmnAxcB3wCKhHssFrXxPZ+S9C5Csvgo4U7yTYS7yhcBP6yra5KuA66Q9BYzey1uqnUpPobQJTgFXiTc77Kg4S0/AOxLGD7FuXFRuHnVObezkDSD0Jay0MwubmH/OwhJ5N3mHwBunDyJOLcTknQWMB9421hG8pV0OPAkcKyZ/Xe74nNThycR55xzLfPeWc4551rmScQ551zLPIk455xrmScR55xzLfMk4pxzrmWeRJxzzrXMk4hzzrmW/T9VFpehX3/eJAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "#set up the analytical function to find T(t)\n",
    "def funcAnalytical(tempInit,tempAmbient,K,time): #function accepts 4 arguments: initial temperature (tempInit), ambient temperaturte(tempAbmbient), the constant K value (K), and the elapsed time (time)\n",
    "    currentTemp = tempAmbient + (tempInit-tempAmbient)*(math.exp((-K)*time)) #use the arguments to solve for Temperature as a function of time (currentTemp)\n",
    "    return currentTemp\n",
    "\n",
    "#function for euler formula\n",
    "def funcEuler(ambientTemp,initialTemp,initialTime,finalTime,K,stepSize): #function accepts 6 arguments: ambient temperature (ambientTemp), initial temperaturte(tempInit), the constant K value (K), and the elapsed time (time)\n",
    "    currentTime = initialTime\n",
    "    currentTemp = initialTemp\n",
    "    while currentTime <= finalTime: #increment the current time by the step size until it reaches the desired time\n",
    "        currentTime += stepSize #increment the current time by the step size until it reaches the desired time\n",
    "        currentTemp += stepSize * (-K*(currentTemp-ambientTemp)) #Add the previous value of the formula to the current value to get the total temperature\n",
    "    return (currentTemp)\n",
    "\n",
    "initialTemp=85\n",
    "ambientTemp=65\n",
    "K=0.611\n",
    "startTime = 0\n",
    "endTime = 10\n",
    "timeStep = .1\n",
    "timeArray = np.zeros(endTime+1)\n",
    "tempAnalytical = np.zeros(endTime+1)\n",
    "tempEuler = np.zeros(endTime+1,)\n",
    "\n",
    "for t in range (startTime,endTime+1,1):\n",
    "    timeArray[t] = t\n",
    "    tempAnalytical[t] = funcAnalytical(initialTemp,ambientTemp,K,t)\n",
    "    tempEuler[t] = funcEuler(ambientTemp,initialTemp,startTime,t,K,timeStep)\n",
    "    \n",
    "plt.rcParams.update({'font.size': 18})\n",
    "plt.rcParams['lines.linewidth'] = 3\n",
    "plt.xlabel('time(s)')\n",
    "plt.ylabel('temperature(F)')\n",
    "plt.plot(timeArray,tempAnalytical,label='analytical solution')\n",
    "plt.plot(timeArray,tempEuler,label='Euler approximation')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 432,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "83.778\n",
      "74.99637927230133\n",
      "70.66782949423398\n",
      "68.0172421597134\n",
      "66.60621455878542\n",
      "65.80281649611554\n",
      "65.42737555550276\n",
      "65.22751135075704\n",
      "65.12111458892966\n",
      "65.0644747772047\n",
      "65.03432284196589\n",
      "65.01827160219376\n",
      "65.00972680080102\n",
      "65.00517801629104\n",
      "65.00275649242322\n",
      "65.00146740567278\n",
      "65.00078116645285\n",
      "65.00041585025761\n",
      "65.00022137591307\n",
      "65.00011784841777\n",
      "65.00006681866704\n",
      "65.00003557060062\n",
      "65.0000189358406\n",
      "65.00001008040498\n",
      "65.00000536625579\n",
      "65.00000285670083\n",
      "65.00000152075114\n",
      "65.00000080956465\n",
      "65.00000043096789\n",
      "65.00000022942372\n",
      "65.00000012213265\n",
      "As t approaches infinity, temperature approaches 65 F\n"
     ]
    }
   ],
   "source": [
    "def funcEuler(ambientTemp,initialTemp,initialTime,finalTime,K,stepSize): #function accepts 6 arguments: ambient temperature (ambientTemp), initial temperaturte(tempInit), the constant K value (K), and the elapsed time (time)\n",
    "    currentTime = initialTime\n",
    "    currentTemp = initialTemp\n",
    "    while currentTime <= finalTime: #increment the current time by the step size until it reaches the desired time\n",
    "        currentTime += stepSize #increment the current time by the step size until it reaches the desired time\n",
    "        currentTemp += stepSize * (-K*(currentTemp-ambientTemp)) #Add the previous value of the formula to the current value to get the total temperature\n",
    "    print (currentTemp)\n",
    "    return (currentTemp)\n",
    "\n",
    "initialTemp=85\n",
    "ambientTemp=65\n",
    "K=0.611\n",
    "startTime = 0\n",
    "endTime = 30\n",
    "timeStep = .1\n",
    "\n",
    "tempEuler = np.zeros(endTime+1,)\n",
    "\n",
    "for t in range (startTime,endTime+1,1):\n",
    "    funcEuler(ambientTemp,initialTemp,startTime,t,K,timeStep)\n",
    "print ('As t approaches infinity, temperature approaches 65 F')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 438,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.8490896782572298\n",
      "The time of death was 0.85 hours before the body was found. ( About 10:09 am)\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "def timeOfDeath(ambientTemp,initialTemp,K):\n",
    "    time = (math.log(98.6-ambientTemp)-math.log(initialTemp-ambientTemp))/-K\n",
    "    return time\n",
    "\n",
    "tempA = 65\n",
    "tempI = 85\n",
    "K = 0.611\n",
    "print (timeOfDeath(tempA,tempI,K))\n",
    "print ('The time of death was about 0.85 hours before the body was found. ( Around 10:09 am)')"
   ]
  },
  {
   "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": 518,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The function looks correct, and although I am sure there is a better way to find the solution, I have no idea what it is.\n",
      "The Euler approximation does not work with times lower than 0 and I can not figure out why.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEWCAYAAACnlKo3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3hUxfrA8e9skk0nnU4IBOkgTToEEVAR6YqAUqQI13JVrt0rcK+Koth+ylU6IqhYEUFBlKIiSOhVaaG3QICEBJLNzu+Ps9lk0zfZZBN4P8+zDzlnZmffsyT77pkzZ0ZprRFCCCGKwuTuAIQQQpRfkkSEEEIUmSQRIYQQRSZJRAghRJFJEhFCCFFknu4OoDSFh4frqKgod4chhBDlyubNm+O11hG5ld1QSSQqKorY2Fh3hyGEEOWKUupIXmXSnSWEEKLIJIkIIYQoMkkiQgghikySiBBCiCKTJCKEEKLIJIkIIYQoshtqiG9RpFqsxMZdcHcYopREBHpzU6VAd4chRLkhSaQAiVfTGDJro7vDcLvmlc0MbBRERf/r+1fmylm4eMqLAO/r+zjFjc3Dw4PAwEBCQ0Px9vYuVlvylyIKVCXAg/Gtw6lUpSp4mlFKuTukEmVSitqVAvD29HB3KEK4nNaatLQ0Ll++zNGjR4mMjCxWIpEkUgBPDxPtaoe5Owy3uqu2FzWqVsbb39/doZSoaxYraelWrFpzIiGFWuH+133CFDcepRRms5nw8HAALly4QJUqVYrcniSRAgT5evHp2LbuDsOt9u/fT1S1cLy8vNwdSolKTrVw8GwSGki6ZuFichoh/mZ3hyVEialQoQJxcXHFSiIyOksUyGKx4Ol5/X/f8DN7EhaQeVp/6lIKlnSrGyMSomR5eXmRnp5erDYkiYhCuVG6dSpV8MHsYfxZWKyaU5euujkiIUqOK/6uJYkIkYWHSVE12Ne+nZCcSuLVNDdGJETZJklEiGwq+HoR7Jt5/efExRSsVu3GiIQouySJCJGLKsG+eJiMU/1Ui5WziaXXrTVixIgS6z6MioqiS5cuJdJ2Vl26dMFdC8BNmjQJpRRxcXEl0r5SihEjRpRI2+WRJBEhcuHlYaJKkI99+1xiKimpxbsAWVomTZrEt99+6+4wyq2LFy8yadIk1qxZ4+5QyoXrf8iNEEUU4mcmITmNK9csaDQnLqYQHVH27x2ZPHkyw4cPp2/fvjnK/vrrrzIfv7tdvHiRyZMnA+R61paSkoKHh9yImkHORITIg1KKasG+9g/d5FQL56+kujmq4vH29sZslntfisPHx+e6v2fKGZJEhLBJTEzkxRdfpE2bNoSHh+Pt7U3jBvWY8eZ/SUlJBuD0paukWqysWbMGpRTz5s1j7ty5NGrUCG9vb2rWrMnUqVNztL1y5UoGDRpE7dq18fX1JTg4mB49erB27doC43rsscdQSrF///4cZadOncLT05NRo0YRFxdnT3jz589HKWV/ZMjrmsjWrVu55557qFSpEt7e3tSoUYPBgwdz8OBBe53PP/+c3r1726fJCA8Pp2/fvuzYsaPAY8jPsmXLiImJITw8HF9fXyIjI+nfvz9///23Q70dO3bQr18/wsLC8PHxoWHDhkydOrVQ9znkd50p6zWONWvWUKtWLcA4o8t4/7Je38nrmsisWbNo0aIFvr6+BAUF0aNHD3777bc8X++PP/4gJiYGf39/wsPDGT16NElJSQUeS1kj3VlC2Jw4cYJZs2YxYMAAhgwZgqenJ2vXruWDd99i69atTP/kK6xac/JiClobo7U+/PBDzpw5w6hRowgODuaTTz7hmWeeoXr16gwZMsTe9rx587hw4QLDhg2jevXq9te67bbbWL16NZ06dcozroceeoj/+7//Y86cOUyZMsWhbP78+aSnpzNq1CgiIiJYsGABDzzwAJ06dWLs2LGFOu7vv/+eAQMG4O/vz+jRo6lTpw6nT59mxYoV7Nq1i+joaADef/99QkNDGTt2LJUrV+bgwYPMmDGDDh06sGXLFm666SZn33LWrl1L7969adKkCc899xzBwcGcPHmSVatWceDAAerWrQtAbGwsMTExeHl58fDDD1O5cmWWLl3KM888w/bt21m4cKHTr52bBg0a8Pbbb/PEE0/Qr18/+vfvD0BAQEC+z3vmmWeYOnUqrVu35tVXXyUxMZEZM2Zw6623smTJEnr27OlQf9u2bfTq1YuRI0cyZMgQ1qxZw+zZszGZTMyYMcMlx1JaVMYfw42gVatWOjY21t1hlDt79+6lQYMGOfZHPbvMDdEUTtxrdzn9nNTUVJRSOboq/v3vf/Pyyy/zyXeraNK8JQBHd2/i7jt6UKVKFfbs2UNwcDAAycnJ1KxZkzp16vDHH3/Y27hy5Qr+2eYeO3PmDI0aNaJ169YsX77cvn/EiBHMnz+frH+b7du3Jy4ujqNHjzrMHlC3bl08PT3Zs2ePfZ9SiuHDhzNv3rwcxxgVFUVUVJT9onFGvEoptm7dSrVq1RzqW61WTCZTnsewd+9emjVrxqhRo5g+fbp9f5cuXYiLiytwhNSTTz7J22+/zZkzZ6hYsWKe9Tp06MDGjRvZsmULTZs2BYyJBAcNGsQXX3zBqlWruO222wBjYMHkyZM5fPiw/Qwit/c0r/crLi6OWrVqMXHiRCZNmlRg/b/++osGDRrQvn17fvnlF3t34cmTJ2nYsCHBwcEcPHjQfh0l4+xm/fr1tG2bOaXSXXfdxcqVK0lISCgwablSXn/fWSmlNmutW+VWJt1ZQtiYzWZ7ArFYLCQkJBAfH0+3bt0AOLRnm73u+STj2sjIkSPtCQTAz8+Ptm3b5uh6yvrhm5SUxPnz5/Hw8KBNmzZs3FjwUgNjx47l1KlT/PDDD/Z969atY//+/YwaNaoIR2tYsWIF8fHxTJgwIUcCAewJJOsxaK25fPky8fHxREREUK9evUIdQ26CgoIA+Oqrr7BYLLnWOXv2LOvXr6d37972BALGh/Hzzz8PwDfffFOk13eFJUuWoLXm6aefdrjeVLVqVUaMGMGRI0fYunWrw3PatWvnkEAAunbtisViKbGhySVFkogQWUyfPp2mTZvi7e1NaGgoERER9msI6VeT8MwyJQpA7dq1c7QRFhbG+fPnHfYdPHiQ++67j5CQEAIDAwkPDyciIoLly5eTkJBQYFyDBg0iKCiI2bNn2/fNnj0bs9nMsGHDinq49mTXvHnzAutu3bqVXr16ERgYSFBQEBEREURERLBz585CHUNuHnnkEZo3b84//vEPQkND6dmzJ++99x7nzp2z1zl8+DAAjRo1yvH8hg0bYjKZOHToUJFe3xXyi69x48YAOeLL6/cGyPG7U9bJNRFRZEXpMirL3nrrLSZMmECPHj147LHHqFq1KmazmRMnThgXUrWmapAPRy8k259jKcT8jElJSXTu3JkrV67w+OOP06RJEwIDAzGZTEyZMoVffvmlwDZ8fX25//77+eijjzh9+jR+fn58+eWX9O7dm4iIiCIfc0b3TkHDfo8ePUrnzp2pUKEC//73v6lXrx7+/sZw58cff7zIF4TDwsLYtGkTv/76Kz/99BPr1q3jiSeeYOLEiSxfvpx27drl2gXlrLyOL6+zH2cUJb78hgiXt0sMkkSEsFmwYAFRUVH88MMPDt04P/74o/3nIF8vKvhkXjNJSE7DqjWmfD6Ef/75Z06ePMmcOXMYOXKkQ9mLL75Y6PjGjh3LBx98wMcff0xQUBDJycnF6soCqFevHmCcZXTv3j3Pet988w1JSUl899133HrrrQ5l58+fL9aiRh4eHnTp0sV+xrdjxw5atmzJyy+/zLJly+zf2nfv3p3jufv27cNqteb6zT6r0NBQwFg7I+NnyHmGAM5PSpgx8GD37t32nzNkXKsqKL7yTLqzhLDx8PBAKeXwTdBisfDaa6/Zt5UyJmjM+MNJS7dyLvFage1Czm+YK1eudOpaQtOmTWndujVz5sxh9uzZREZG0qNHjxz1AgICuHDhQqHa7NGjB+Hh4UybNo1Tp07lKM+IOa9jmDlzJqdPny70MWQXHx+fY1/9+vXx9fW1H0PFihVp3749S5cuZdeuXQ6xZYxW69evX76vkzHKa9WqVQ77p02blqNuxkXtwr6HvXv3RinFG2+8QVpa5mSdp06dYu7cudSsWbNQ3YXllZyJCGEzcOBAnnvuOe6880769+/P5cuXWbRoUY7RWmZPE6FZFqs6m3iNYF8vvL1y76Lo2LEjlStXZsKECcTFxVG9enW2bdvGggULaNKkCTt37ix0jGPHjmX06NEATJw40eGMKUPbtm1ZtWoVr7/+OpGRkSiluO+++3Jtz8/Pj9mzZzNw4EAaN25sH+J77tw5VqxYwZNPPkmfPn2488478fPz44EHHuCRRx4hJCSE33//neXLlxMdHV3kbqExY8Zw/PhxevToQc2aNUlJSeHzzz8nMTHR4VrPu+++S0xMDJ06dbIP8f3+++9ZsWIFQ4YMsY/MysvgwYN5/vnnGTt2LPv27SMsLIwffvgh1yQWFhZGnTp1+Oyzz4iOjqZSpUr4+/tz991359p2vXr1eOqpp5g6dSqdO3dm0KBB9iG+SUlJLFy48Pq+w11rfcM8WrZsqYXz9uzZ4+4QSoXFYtGvvvqqjo6O1mazWUdGRuqnnnpK79mzRwN64sSJ9rq//PKLBvR/pn2gtx9L0AfOJmqr1aq11nr48OHa+NPKtH37dn377bfr4OBgHRAQoGNiYvS6detyrZvbvgxJSUm6QoUK2mQy6bi4uFzr/P3337p79+46MDBQAw5t1axZU8fExOR4zsaNG3WfPn10WFiYNpvNukaNGnrw4MH64MGD9jpr167VHTp00AEBATooKEj37NlT79y5U8fExOiaNWs6tJfbvtx89dVX+u6779bVqlXTZrNZh4eH686dO+svv/wyR91t27bpPn366JCQEG02m3X9+vX166+/ri0Wi0O9iRMnakAfPnzYYf+GDRt0+/bttbe3tw4LC9NjxozRCQkJGtDDhw/P8X60b99e+/n5acDhWHKrr7XWM2bM0M2aNdPe3t46MDBQd+vWTa9bty5HvbyeP3fuXA3o1atX5/V2lYjC/H0DsTqPz1W5T0QUqDDjyG9EKakWDpy9gsb4G6oe4udwhlISrl27RpUqVbjllltYsWJFib6WuDHIfSJCuImv2ZPwwMykcepSCmklvJzuwoULSUhI4KGHHirR1xGisOSaiBDFUDHQh0spaaRarKTbltONDPVz+essXbqUI0eOMGnSJBo2bEifPn1c/hpCFIUkESGKwcNkzPR7OP4KABeTUwnx8yLQx7WzvD766KOcPHmSli1bMmvWrOv7Qq0oVySJCFFMgT5eBPuZuZhsTIVyIiGFmyp52ldGdIXyNhWGuHHINREhXKBKkE/mcrrppbucrhDuJElECBcwltP1tW/HJ14jJbX4U2oIUdZJEhHCRUL8vAjwNnqINXA8IXPdESGuV5JEhHCR7MvppqSlE59UvpfTFaIgkkSEcCFvLw8qBWZORnjm8lVSLQUv3ypEeSVJRAgXCw/0xsc2j5axnO5V6dYS1y1JIkK4mMnWrZXh8tU0LqWk5fMMIcovSSJClAB/b0/C/DO7tU5evIqlhKdEEcIdnEoiSik/pVQPpdQTSqkpSqlXbT93V0q5fq4HIcqxykHeeNmW033un+Pw8pS7zEtTVFSUfaGrsiguLg6lFJMmTXJ3KMVSqDvWlVJ3AuOAO2zPyX4rrgYsSqkfgA+11j8iRDmzZs2aHKv2ZeXh4eHUuhkeJhNVg305cv6KfV/SNYt9GLC4/sXFxTFv3jz69u1Ls2bN3B1Oicj3t1kp1QmYBrQC4oA5wB/AQeA8RjIJBeoA7TCSzHKl1GbgCa31byUWuRAlZPDgwfTs2TPH/twWgCpI9uV0jSlRAvJdTle4xl9//eX0UreuFhcXx+TJk4mKisqRRDIW4fL0LN9fKgqKfg3wLTBBa/1rPvV+B+YDKKVigMdtzy3f7464IbVo0YL777/fZe1VDfYl47PsmiWdc4nXqFTBp9jtpqWlkZ6ejo9P8dtyh8TERAIDA0us/eKs+14alFLl9v8uq4K+WrXQWg8oIIE40Fqv1Vr3A1oWLzQhyq41a9aglGLevHk5ykaMGOHwDdjsabIP+QVjOd2raca9I6dOnWL8+PFERkZiNpupWrUqY8eO5ezZsw5tTpo0CaUUu3fv5sknn6R69er4+PiwYcOGPGO0Wq288sordO7cmcqVK2M2m4mMjGT8+PGcP3/eoW7W/vlPP/2Upk2b4uPjQ2RkJJMmTcrRjZdxjOfOnWPYsGGEhYXh7+/PbbfdxtatW/Ns+/PPP6dly5b4+vry6KOP2uvs2LGDfv36ERYWho+PDw0bNmTq1Kmkp2feY3Ps2DHCwsJo0qQJKSkpDq8xdOhQTCaTwxrquV0Tydi3fft2unXrRkBAABUrVuRf//oXFouFq1ev8q9//Ytq1arh4+ND586d2bt3r0MbiYmJvPjii7Rp04bw8HC8vb2pU6cOzz77LMnJyfZ68+bNs3ePjhw5EqUUSil7THldE7FYLLz++us0bNgQHx8fwsLC6NevX45llLM+//vvv+eWW27Bx8eHKlWq8NRTTxV5yWJn5XumoLXeXtSGi/NcIdwpOTk517W3zWYzFSpUKFKbWZOI1poTCSl4ppynffv2pKamMmrUKKKjozlw4AD/+9//WL16NbGxsQQFBTm0M3ToUHx9fZkwYQJKKapUqZLna6ampvLGG28wYMAA+vTpg7+/P5s2bWL27Nn89ttvbN68GbPZcSXGpUuX8s4779jXMf/uu++YPHkyR44cYe7cuTle44477iA0NJRJkyZx+vRp3n//fTp37swff/xB48aNHep+++23vPfee4wfP55x48bZ38vY2FhiYmLw8vKyv+7SpUt55pln2L59OwsXLgSgRo0azJ07lz59+vD444/z0UcfATBnzhwWLVrEs88+S7du3Qr8vzh+/Djdu3dn0KBBDBw4kJUrVzJt2jQ8PDzYvXs3KSkpPPvss8THx/Pmm2/St29f9u7da+/OPHHiBLNmzWLAgAEMGTIET09P1q5dy9SpU9m6dat9xcnOnTvz/PPP8+qrrzJ27Fg6deoEQKVKlfKNb+jQoSxevJju3bszfvx4Tp8+zQcffEC7du349ddfad68uUP95cuXM336dMaNG8eDDz7IkiVLePPNNwkJCeH5558v8P0otrzWzc14ABUBc0H1ysND1lgvmjzXYJ5Yoew+imD16tX2Nclze9x111056s6dOzdHO/mtm77j2EW9/ViC3n4sQd/Rs5eOiIjQx44dc6i7adMm7eHh4bCme8a64TExMTotLa1Qx2O1WnVycnKO/bNmzdKA/vzzz+37Dh8+rAFtMpn05s2bHdro27evBvQff/yR43j69etnX1tea61jY2O1UkrffvvtOdr29PTM9Xepffv22sPDQ2/fvt3hde+55x4N6FWrVjnUf/TRRzWgv/jiC713717t7++v27Ztm+N9yW09+Zo1a2pAL1682GF/ixYttFJK9+7d2+F43n33XQ3oH3/80b7v2rVrOjU1NcdxvPjiixrQGzdutO/L7/ck433J+v+8cuVKDeh7773XIY7t27drDw8P3bFjxxzP9/Pzc1hP3mq16kaNGunKlSvneM3cFHeN9cJcKTwFDMzYUEr52ob21i5W9hKijBo7diw//fRTjscrr7xS7LYjbMvpJl6+xMofl9Pr7rvx8fEhPj7e/oiKiqJOnTqsXLkyx/Mff/zxQl+IVUrh62vc9Jiens7FixeJj4+na9euAGzcuDHHc7p3706LFi0c2nj66acB+Oabb3LUf/rppx267lq2bEn37t1ZtWoVSUlJDnXvuuuuHGt5nz17lvXr19O7d2+aNm3q8LoZ36Kzv+4bb7xB8+bNGTNmDAMHDsTT05NFixYV+n2pVq0a99xzj8O+jh07orXm0UcfdTiejLOH/fv32/eZzWa8vIzBEhaLhYSEBOLj4+1nQbm9r4WVcawvvPCCQxxNmzalV69e/Pbbb5w7d87hOX379iUqKsq+rZTi1ltv5fTp0zn+D0pCYd717MMb/IBngFXAIZdHJISb3XTTTYXqFimKjOV0dx48gNVqZe6cOcydMyfXurVr5/yeVrduXadeb/HixUybNo2tW7eSluZ413xCQkKO+tk/5AEaNmwIwKFDOf/c86q/cuVKjhw5QqNGjfKN/fDhwwAO9bK2YzKZcryut7c3n376KY0aNWL37t0sXLiQWrVq5Xh+XnKrGxISkmtZxv7s15CmT5/Ohx9+yO7du7FaHW8ize19LazDhw9jMplyfV8bN27MkiVLOHz4MBEREfb9uf2ehIWF2eMOCAgocjyFUdTRUzI+UcCkS+6OwG3yGzqa3wVNk2053U0Yc2nd1f9exjw4Ev9c7h3JOIvIys+v8Pf0fv311wwaNIjWrVvz7rvvUqNGDXx8fEhPT+eOO+7I8eEH+R9XYek85gnLLfa86hZk2bJl9ovuW7duZciQIYV+bn5LC+dVljXOt956iwkTJtCjRw8ee+wxqlatitls5sSJE4wYMSLX97WwivJ+5Hc8RX1/nSFDcIUogtDQUAAuXLiQoyy3b+xZBfh40bRBPZRSpKWmUbdFe+pWCnTpcroACxYswMfHh9WrVzt8gO/bty/P5+zZsyfPfbl94927dy9t27bNsc/Dw4OaNWsWGGNGm7t3785Rtm/fPqxWa47X3bx5M8899xzdu3cnPDycadOm0b17d3r06FHg67nCggULiIqK4ocffnC4d+jHH3PeY+1sUo6OjmbFihXs3bvXoXsPMv8fnDnrKg0yd5YQRVCrVi08PT0dhpQCrF+/Pt9htxka1qpG5649+OXHpWzetJEzlx2X09Va5+j7dpaHhwdKKYdvxlprXn755Tyf89NPP7FlyxaH+lOnTgWMvvfspk6d6vBtd8uWLaxatYrbbrutUN0oFStWpH379ixdupRdu3Y5vO6UKVMA6Nevn31/UlIS9913H8HBwSxYsICPPvqIWrVqMWzYsBzDoktKxvua9bgtFguvvfZajroZ70FuXzZyk/EeT5kyxaH9Xbt28d1339GxY0eHrqyyoLBnIq2UUhm/5Rl3B3VUSgXnVllr/XWxIxPCTbZs2cInn3ySa1nfvn0JCAggICCAESNGMGvWLAYPHkyXLl3Yv38/c+fOpWnTpmzfnv8Id08PE/83/QNu79qFkQPv4u4B99GpbSs8TYpDhw6xZMkShg0bVqx5lQYOHMhXX31F165dGTZsGGlpaXz77bcO9zJkd/PNN9O1a1cefvhhqlSpwpIlS1i1ahUPPPAA7dq1y1H/yJEj3H777fTu3ZtTp07x/vvv4+vryxtvvFHoON99911iYmLo1KmTfYjv999/z4oVKxgyZAi33Xabve748eM5ePAgP/74o32o7GeffUaHDh0YPnw4y5cvL/G71AcOHMhzzz3HnXfeSf/+/bl8+TKLFi2yX2zPqmHDhgQGBjJ9+nT8/PwIDg6mYsWK9sEN2XXv3p17772Xzz77jISEBHr16mUf4uvj48N7771XosdWJHkN28p4AFYgPdvDms/+9ILadNdDhvgWTWGGAF4PChriC+j9+/fb6ycmJurRo0fr0NBQ7evrqzt06KB///33fIf4ZmW1WvWmvYf18Ice0ZG1orXZ21sHBQXpxo0b68cee0zv3r3bXjdjiG/WoZyFMWPGDN2gQQPt7e2tK1eurMeMGaPPnz+vAT18+HB7vazDTRctWqSbNGmizWazrl69uv73v/+dY0hrxvGcPXtW33///fb34NZbb9WxsbEOdXMbyprdtm3bdJ8+fXRISIg2m826fv36+vXXX9cWi8VeZ/78+RrQTz/9dI7nT506VQP6zTfftO/La4hv9n1a5/3+5ha7xWLRr776qo6OjtZms1lHRkbqp556Su/ZsyfX41y2bJlu3ry59vb2tg/Tzu99SUtL06+99pquX7++NpvNOiQkRPfp00fv2LGjwNgKOp7cFHeIr9IFXHhRSg0vQmKa7+xzSkOrVq10bGysu8Mod/bu3ZvraBFRfNcs6ew/k4TV9ndYJciHiMDSnwojLi6OWrVqMXHixEKd/YwYMYL58+eXyoVbUbIK8/etlNqstW6VW1mB3VllNSEIcT3w9vSgYgVvTl8yeovPXL5GBV8vvGXaeFFOyIV1IdwsPMBxOd0TCSnyDV+UG/kmEaWUc3c2OT63XlGfK8SNxKQU1UMy7wlJumbhoiynK8qJgs5Ediul5iilGhdQz04p1VwptQDYVWBlIQQAfmZPwgMypy4/VcrL6UZFRaG1LvRosHnz5snZkgAKviZyN8aiVNuVUjuAZcAmjEWpLpC5KNVNQFvgTqAhsAfoVUIxC3FdqlTBmBIlLd2KxWrl1KWr1AiVVadF2VbQVPA/KqV+Au4F/gE8D+T29SNjYPYa4D7gK6116X2NEuI64GGbEiXOtpxuQnIqIX5eBPjkvP9AiLKiMKOz0oFPgU+VUpWAGIyzjQiMhHIOo+tqrdY65yIM4rqgtXb7UqM3ggq+XgT5enHJdk3kxMUUbqroicnFU6IIAa6ZW8upubO01meAxcV+VRtbUpoM3AVUAk4D3wATtdYXs9WtB7yOkcTMwBZbvV9cFY/InaenJxaLJdc7coXrVQ32JemahXSr5prFytnEa1QOKv/LqIqyJy0tLd8JHAvDbUN8lVIVgY3AgxjruD8KLAHGA6uVUn5Z6kYD64F2wFTgKSAAWKGUKpk5u4Wdj49PqaxLIAxeHiYqZ1mD/VyW5XSFcKXLly8Xe517p5OIUipQKfWSUuo3pdR+pVQ72/5w2/76hWzqeaAmMFxr/ajW+iOt9aPAcKAZ8GSWulOAYOB2rfUUrfV0oBNwEvhAST9LiYqIiODcuXMkJyfLiJxSEupvxt9sdBRoNMfl3hHhIlprUlNTiY+PJyEhwT4jdVE51Z2llIoAfgNqAwds//raAou3TZESjGMCyMutQArwWbb9nwNzgJHAy0opf6A3sEZrvS2jkkoUNxEAACAASURBVNY6SSk1C/gPcAvwpzPHIgrPx8eHSpUqcfr0aa5du+bucG4YaelGV1ZG7kg46UVALuuOCOEsDw8PAgMDiYyMxNvbu+An5MPZ38iXgcpAG+AokH3u5SXAbdmflAdv4KrO9vVKa21VSqUAtZVS4RjDh72BP3JpI2PObUkiJSwoKIigoCB3h3HDeeunv3nvZ2Np1gBvT1Y9GSPXR0SZ4mx3Vi9gutZ6C7kP9T0E1ChkW7uBEKVUs6w7bdshts1IoKrt5xO5tJGxr1peL6KUGquUilVKxRZ3fQYhSts/ukRTO8IfMO5kn/id3MMryhZnk0g4RjdWXqxAYb8mvWOrv1gp1VMpFamUuhOjOytjzgc/2wMgt36Uq1nq5UprPUNr3Upr3aqsLeYiREF8vDx4tV8T+/aK3Wf4cddpN0YkhCNnk8hpIDqf8uYY3VwF0lr/inFjYiDGnfBHgKXAauB7W7XLQMYKOrl13GUkrLxX2RGinGtbO4z7bsk8wZ/43S4Sr8rcWqJscDaJLAdGKaWqZC9QSrUBhmFcFykUrfUXQHWM5NMZqKq1HmfbZ8E46zlpq55bl1XGvty6uoS4bjx3ZwP73FpnLl/jjRV/uTkiIQzOJpHJGB/uWzGG3WpguFLqU2Adxgf+6840qLVO11pv01r/qrU+q5SqjJFU1mqtk4GdGF1ZOdfmNObrApCVpsR1LcjPi4l3N7RvL9hwhM1HEtwYkRAGp5KI1vo0xgd3xk2CCngAY26tlUAnrXXhVqTPhVLKBLwHeACv2F4zCaObq4tS6uYsdQOA0cB+ZGSWuAH0alqFW+sZ1/W0hue/3kmqRaaoE+7l9M2GWutjWus+GLP3tsFIKhFa67u11scL245SKkAptUcp9YpSarRSagJGMrgHeFFrvTpL9eeAS8BKpdSzSql/AL9idGc9mn2YsBDXI6UU/+3bGF/bAlZ/nUlk5q+H3ByVuNEVOonYPvTnKKXuAdBaX9Zab9Ja/1nEs49UYAcwBHgfeAE4D9yhtX41a0Wt9QGgA8Z9Ic8CbwJXbHVXFOG1hSiXqof4MaFH5lpx7/68n8PxV9wYkbjRFTqJ2LqV7gMquOKFtdapWuv7tNa1tNY+WutQrfXteSUFrfVerXUfrXWw1tpPa91Ra73KFbEIUZ6MaB9F42rGn2GqxcoL3+yUKVGE2zjbnbUHiCqBOIQQheTpYeK1/k3JmB1+/cHzfLVFBigK93A2iUwFxhdn7XUhRPE1rhbEqI617NsvL9vD+SSZ10yUPmfnzqoPHAN2KqW+xxgZlf1GP621/q8rghNC5O2J7nVZvvM0Jy6mcDE5jZeX7eXtQc0KfqIQLqSc6UtVShVmPKHWWhdvlZMS0qpVKx0bK7eUiOvH6r/OMnLuJvv2xw+2pnNdmd5HuJZSarPWulVuZc52Z9UqxKN20UMVQjjj1noV6X1zVfv2C9/uJCVVFrASpcfZmw2PFOZRUsEKIXL6d6+GBPkaSxcfu5DCOz//7eaIxI3EbcvjCiFcIyLQmxd6NrBvz/r1MLtPXnJjROJG4uzKhnMKUU1rrUcVMR4hRBHc06o6X289zoZDF0i3ap77eiff/KMDHiZZOVqULGdHZ40oRB0NSBIRohQppXi1XxPuePdXUi1Wdhy/xPz1cTyYZRiwECXB2WsipuwPwAuoB8zEmJYkJN9GhBAlonZEAI/eWse+/ebKvzhxMcWNEYkbQbGvidimct+vtX4IY+4rp6aCF0K4zkMx0dxUMQCA5NR0Xvp2l0yJIkqUqy+s/wAMcHGbQohCMnuamNI/czndn/ed5QdZTleUIFcnkTAgwMVtCiGc0CoqlKFtIu3bE7/bzaUUWU5XlAyXJBGlVLBSaiDwBLDZFW0KIYru6TvqUzHQWE73XOI1Xv9xn5sjEtcrp5KIUsqqlErP/sC4FrIYsAJPlkSgQojCC/L1YnLvRvbtRRuPsimuyIuOCpEnZ4f4fowxhDcrDVwA/gY+1VonuiIwIUTx3NG4Mt0aVGLV3jMAPPf1TpY91hFvzzI5tZ0op5xKIlrrESUUhxDCxZRS/KdPI/44GM+V1HQOnE3iwzWH+Ge3m9wdmriOONud9ZJSqnE+5Y2UUi8VPywhhCtUDfblqdvr2bc/WH2AA2eT3BiRuN44e2F9EtA0n/LGwMQiRyOEcLkH2kVxc41gAFLTrTz/9U6sVrl3RLiGq4f4+gAWF7cphCgGD5Pitf5N8LTNo/Vn3AUWxx5zc1TielFgElFKVVBKRSqlMgaeh2VsZ3s0A4ZirHwohChDGlSpwJjOmUv9vLp8L2cTr7oxInG9KMyZyBPAYdtDA+9k2c762Ax0Az4skUiFEMXyz9tuIjLUD4DLVy38Z+keN0ckrgeFGZ21xvavAl4CvgF2ZKujgSRgg9Z6vcuiE0K4jI+XB6/0a8wDs/8E4PsdpxjQ4iy31q/o5shEeVZgEtFarwXWAiilagIfaq03lnRgQgjX63RTBP2bV+PrrScAePHbXax8ojP+3s7eMiaEwdn7REaWVCBlltagre6OQpQWZQJ1fS/k9MJdDVj911kSktM4cTGFt3/6mxd7NXR3WKKcKtLXD6WUB1AfY+2QHNdVtNbrihlX2ZF8Ht6IdncUorR4mCEkCkKjIbQ2hNaCMNvPQTXAVP7v9g4L8ObFuxoy4YvtAMz5/TA9GlWmda1QN0cmyiOnk4hS6hngWaBCPtXK/1+auDGlp0L838YjO5OXLcHUzkwsobWMhBNUAzzKT5dQ/xbV+HrrcX4/cB6rhgfnbeLjUa1pESlrygnnOLvG+mhgCsY1kpXAK8DbQBrGkriHgOkujtH9lKtvpxFlVn5dl9Y0OL/feOzPVmbygpCatsRS2/FMJrhmmUswSile69+Uvh/8zvkrqSRdszB89p8sGN2GZrYbE4UoDOXMqmdKqVggVWvdXikVBpwDummtf1FKVQG2Ac9preeUTLjF06pVKx0bG+vuMERZdvUyJByG8wfhwiG4cBgu2H5OOlO0Nk2eEBzpmFwyzmSCI8HDy7XH4IS/TicyeOYGLlxJBSDQx5OFo9vQtLokEpFJKbVZa90q1zInk8gV4EWt9dtKqVAgHrhDa73SVv5foLfW+mYXxO1ykkREsVxLsiWWQ5mJ5YIt4SQVcfVA5ZElwWTtJos29nuaXXsMudh3+jKDZ2wgIdlYuKqCjyeLxrSlcbWgEn9tUT7kl0ScPcdOB67Yfs74N+vVuDhApggV1yfvAKjS1Hhkl3rF8azl/EHb9iFIPJl3mzrdOPNJOAwHf3YsUybjWos9sWQ5kwmpCZ7eLjms+pUrsHB0W4bM2sDF5DQuX7UwdNZGFo5uI4lEFMjZM5FdwFKt9XO27SPA91rrh23bHwK9tNbVSyLY4pIzEeEWqVcgIS5LcjmU+bh8omhtKhMEVc92/cV2JhNcE7x8nG5y98lLDJm50b6UbrCfF4tGt6Vh1fzG0IgbgSu7s6YDHbXWTW3bbwKPYyxWZQLuB+ZorccWO+oSIElElDlpKZlnLFm7yc4fgsvHi9ioMs5gQmtl6yKrDSG18k0wu05cYsjMDVy+asyjGuLnxaIxbWlQRRLJjcyVSaQe0AX4WGudopTyBz4F7rJVWQkM1VqXyXU4JYmIciUtBRKOZLn+kqWb7NIxci4yWhgKKlRzvP+lZkeo3tJeY+fxSwyZtYFEWyIJ9Tfz6Zi21Ksc6JrjEuWOy5JIPi8QBKRrrcv0ajeSRMR1I+0qXDySSxfZQbh03PlZFjr9C7q+aL9bf/uxi9w/ayOJ14xEEuZv5tOxbalbSRLJjcglSUQpFQC8B/ygtf7ChfGVGkki4oZguWY7g8k+kuwQXDyad4Lp8jx0eca+ufVoAsNm/2lPJOEBxhnJTZJIbjiu7M5KBh7VWs92VXClSZKIuOFZUo1EkpFY9i6FI79nlt/2EnSaYN/cfCSB4XP+JMmeSLz5bGxb6lQMKO3IhRvll0ScvRV7DxBV7IiEEO7haYbwOlD3dmg7Hu7/GqK7Zpb//B/4/T37ZsuaIcx/8Bb8zcZMRvFJ1xg8cwMHz5XpnmtRipxNIlOB8UqpuiURjBCilHn5wH2LoFbnzH0//Rv+yJy9qGXNUOY92Bo/WyI5l3iNwTM2cEgSicD57qyXgH5AQ+B7jBmEkrNV01rr/7osQheS7iwh8pB6BRbe49i11fNNaD3Gvrnx0HlGzN1ESlo6AJUqePP52HZEhfuXdrSilLnymkhhhnxorXWZnMVXkogQ+biWCJ8MgGNZ1pzr9Q60ylxGaMOh84zMkkiqBPnw2di21AyTRHI9c+U1kVqFeNQueqhCCLfxDoShX0K1LJ8V3z8OWz+xb7atHcbsEa3w8TI+Ok5dusrgGRs4ej57h4S4UTiVRLTWRwrzKKlghRAlzKcC3P8VVG2euW/JI7D9M/tm++hwZg+/BW9P4+Pj5KWrDJ65gWMXJJHciIq8UIZSqo5SqoPtRkMhxPXCN9gYtVW5iW2Hhm/Hw84v7VU61Aln1vBWmG2J5MTFFO6bsYHjCZJIbjROJxGlVC+l1EHgL2Ad0NK2v6JS6oBSaqCLYxRClDa/UHhgCVRsZGxrK3w9FnZ/a6/S6aYIZg5zTCSDZ27gxMUUd0Qs3MSpJKKU6gJ8A1wAJgMqo0xrfRY4CNznwviEEO7iHwbDlkBEfWNbp8NXo2DfMnuVmLoRzHigJWYP46Pk2IUUBs/YwElJJDcMZ89EXgK2A22AD3Ip/wNoUdyghBBlREAEDPsOwmzLBFktsHg4/PWjvUqXehX5KEsiOXohmcEzN3D60lV3RCxKmbNJpBWwUOs8Z3c7DlQuXkhCiDIlsBIMX2rM+AvGWvOLH4ADq+xVbq1fkf/d3wIvD6Nz4sh5I5GcuSyJ5HrnbBLxAK7lUx4OpBY9HCFEmVShipFIgmsa2+mp8NlQOLTGXuW2BpWYPrSlPZEcjr/C4BkbOCuJ5LrmbBLZC3TKp7wXRneXEOJ6E1TdSCRBNYxty1VYdB8c/tVepXvDSrw/pAWeJiORHIq/wuCZGzibKInkeuVsEpkNDFRKjcryXK2U8lNKvQe0A2a4MkAhRBkSUtNIJBWqGduWFFg0CI78Ya9ye6PKvD+kOR62RHLw3BWGzNzIucT8OjFEeeXszYb/Az4HZmLMm6UxVja8BDwCzNNaL3R1kEKIMiS0lpFIAmyXP9OuwMKBcOxPe5U7Glfh/wZnJpIDZ5MYOmsD8UmSSK43Tt8norW+HxgA/Azswxjuuxy4R2s9yrXhCSHKpLBoI5H4VzS2U5OMebdObLZX6dmkCu/e18yeSP4+k8TQmRs5L4nkulKkO9a11t9orQdorRtprRtqrftorb9ydXBCiDIsoi4M/w78wo3ta5dhQT84uc1epVfTqrw9qBm2PMJfZxIZOmsjF67I+JvrRZGnPQFQSvkqpXxdFYwQopyp2MC4IdE3xNi+egkW9IXTu+xVet/smEj2nTYSSYIkkutCUaY9qaiUmq6UOgkkAUlKqVO2fZVcH6IQokyr3NhIJD62afRSEuDj3nB2r71Kn2bVmHbvzShbItl76jL3z97IxWRJJOWds9Oe1AK2AuMwLqYvAb4DLtr2bVFKyVTwQtxoqtwMD3wL3hWM7eTzML83nPvbXqVf8+q8MTAzkew+aSSSS8lpbghYuIqzZyLTgDCgv9a6gda6v9a6n9a6AcbF9jDgzcI2ppQKUEo9r5TaqZRKVErFK6XWK6VGKKVUtrr1lFLfKqUSlFJXlFK/KqW65tW2EKKUVWthzP5rDjS2r5yF+XfD+YP2KgNbVuf1AU3tiWTXics8MGcjl1IkkZRXziaR24APtNbfZi/QWn8D/M9Wp0BKKRPwA/BfYBMwAXgZ4674ucBrWepGA+sx7kOZCjwFBAArlFLdnDwGIURJqXEL3P8leNlWOkw6DfN6wYVD9ir3tqrBa/2b2Ld3HL/EsNkbuXxVEkl55GwS0Rj3h+Tlb1udwmgDdATe01o/qLWeobV+B+OO+MPAQ1nqTgGCgdu11lO01tNt9U4CH2Q/axFCuFFkWxi6GDxtY24STxpdWwmZ69UNuiWSKVkSyfbjlxg2+08SJZGUO84mkbXArfmUdwHWFLItW+cpJ7Pu1FqnAvHAFQCllD/QG1ijtd6WpV4SMAuoC9xSyNcUQpSGqI4w5DPw9DG2Lx0zurYuHbdXGdw6klf6NbZvbzt2keFz/iTpmqW0oxXF4GwSeRxoo5SappSqmLHTNmLrLYyzi8cL2dafGBfkn1ZK3aOUirRd95iCsdDVJFu9poA3xjTz2W2w/StJRIiypnYXuG8heJiN7YtHjK6ty5nfG4e2qcl/+zSyb285epERkkjKFWeTyM+AL0aiOKWUOq+UigdOAf+0lf2ilDqU5XEwt4a01gkYZxgXgMXAEYw74B8GBmitZ9qqVrX9eyKXZjL2VcsrYKXUWKVUrFIq9ty5c84cqxCiuOp0g0ELweRlbCccNs5IEs/YqzzQLorJvTMTSeyRBEbO/ZMrkkjKBWeTyFFgD8ayuOuAHcDOLNt7MJJB1sfRfNpLAnZhjOjqD4wGDgCLlFLdbXX8bP/mNlfC1Wx1crBda2mltW4VERFR0PEJIVytbg+492MweRrb5w8YiSQp80vd8PZRTLy7oX17U1wCI+dtIjlVEklZ5+lMZa11F1e9sFKqCcaIqye01h9m2f8pRmKZaRuVlWwr8s6lGVuHq72OEKIsqt8TBs6BL0Yay+zG/wUf97HNvxUGwMgOtUi3al5eZtyk+OfhCzw4bxNzR7TG1+zhzuhFPoo17UkxPYGRBL7IulNrnQwsA2oCUWReeM+tyypjX25dXUKIsqRhHxgwE5TtY+fsbljQB5Iv2KuM7lSbF3o2sG9vOHSBUfM3kZKaXtrRikIqchKxrSFSw3ZB3OFRyCYyEkBuXzE8s/y7E6Mrq10u9dra/o0tdOBCCPdpPAD6fgjYRuWf3mlM2phy0V5lTOfaPHdnffv2+oPnGfNxLFfTJJGURc5Oe+Jhu8P8BJAIxGHc05H9URh7bP+OyPYawUAfIAE4aBvKuxToopS6OUu9AIxrKPsxRnoJIcqDmwdBnw+wJ5JT24xp5K9etld5KCaap++oZ9/+7UC8JJIySmld2HsDQSn1LvAosAX4DeODPget9eRCtFXT1k4IsBD4HQgFxmB0Yz1su6kQpVQdjESRBrwNXLbVawLcpbVeUZj4W7VqpWNj5aRFiDJh8zxY+s/M7Rpt4P6vwDvQvuuD1Qd4Y8Vf9u3OdSOY8UBLfLzkGklpUkpt1lq3yrXMySQSj3HT30AXBRYNvIQxVUolIAXYBryjtf46W90GGFOhxABmjAQ0SWu9qrCvJ0lEiDLmz5mw/F+Z25HtjWlTzP72Xe/9vJ+3fsqcyLFLvQg+eqAl3p6SSEqLK5PIJeAprXW5XEddkogQZdCGD+HHZzK3ozrBkMVgzhy5/86qv3lnVeaMS13rV+R/97eQRFJK8ksizl5YXw80LLCWEEIUVttx0OPlzO24X+GzIZB21b7r8W51eaxrHfv2L/vO8vDCLaRarKUZqciFs0nkaWCIUqpPSQQjhLhBtX8UbpuYuX1oNXx+P1gy7zF+ontdHrk1M5Gs2nuWhxdJInE3p5KI1nonxgXtr5RSR5VSa5VSv2R7/FwyoQohrmudnoQuz2duH/gJFg8Hi7H6oVKKCT3qMr5LtL3KT3vO8OinW0hLl0TiLs4O8e2JMc+VCWMW3kigVraHrGwohCiaLs9A56cyt//+Ab4cCenGFPFKKZ6+vR4PxWR+zKzYfYbHPt0qicRNnO3Oeg04BtystQ7WWtfK7VECcQohbhS3vgAdskwGvu97+HoMpBvzaCmlePaO+ozplPlR88Ou0zz+2TYskkhKnbNJ5CaMRaR2lkQwQgiBUtBtErR7JHPf7m/g23FgTbdVUTzfswGjOmYmkmU7T/HE4u2SSEqZs0nkCJmTHgohRMlQyhix1TrLAqc7v4Alj4DVaquiePGuBozsEGWvsnT7SSZ8sZ10a+FvXRDF42wSeQ8YbZtyRAghSo5ScOfr0OrBzH3bF8H3/3RIJC/1asjwdjXtVZZsO8m/JJGUGqemgsdY/+MisFcpNRdjnqwck9lorT92QWxCiBudUtBzGlgtsMX2sbLlY2ORq7umgVIopZjUuxFWDQs2GOu4f7P1BErBGwNvxsOk3HgA1z9nk8i8LD+/mEcdDUgSEUK4hskEvd41LqxvX2Tsi51tLHJ15+v2RDK5dyOsWrNwo7EO3tdbTmBSiqkDmmKSRFJinE0it5ZIFEIIkR+TCfq8D9Y049oIwJ8fgYeXce1EKUwmxX/7NMaqNZ/+eQyALzcfx6Tgtf6SSEqKsysbri2pQIQQIl8mD2MtEqvFGK0F8Mf7RiK5baI9kbzStwlWK3weaySSxbHHMSnFq/2aSCIpAcVZlMpbKVVNKWV2ZUBCCJEnD0/oPxMa3J2577e3YfWr9k2TSTGlfxPuaVndvu+zTcd4cckurHKx3eWcTiJKqRZKqV8wFqU6CnS07a+olPpZKdXNxTEKIUQmDy8YMAfq3pm5b91UWDvVvmkyKV4b0JQBLTITyaKNR3npu104M3O5KJiz0540A34Fosl28VxrfRbwBYa7LDohhMiNpxnunQ91umfuW/0K/PqWfdPDpJg6sCn9m1ez7/tkw1EmfrdbEokLOXsm8h/gJNAIeBb7+pZ2PwOtXRCXEELkz9MbBn0CtbOM9/l5Mqx/377pYVK8cc/N9GlW1b7v4z+O8NKS3aSkylK7ruBsEukEzLSte55bKj8KVM1lvxBCuJ6XD9y3yFjIKsPKF4yFrmw8TIpp99zM3TdnfjQt2HCEjq//wv/9vJ9LyWmlGfF1x9kk4gNcyqe8QjFiEUII55n9YMjnxtK6GX58BjbNsm96eph4+96buatpFfu+81dSmfbT37R77Wde/n4Ppy6llGbU1w1nk8hBoGU+5V2BPUUPRwghisDsD0MXQ402mfuWTYDN8+2bnh4m3h3UjMm9G1Et2Ne+Pzk1nVm/Habz1NX864vtHDibWJqRl3vOJpFFwAPZRmBpAKXUBOAOYIGLYhNCiMLzDoShX0K1LN9zl/4Tti2yb3p6mBjePoo1T3XhrXtvpl6lQHtZWrrmy83H6fbWOsZ8HMvmIwmlGX25pZwZpWC7J2QF0BnYB9QHdgIRQGXgJ6Cn1rpMzsXcqlUrHRsb6+4whBAlKeUifNwHTm2z7VDQfwY0vTdHVa01q/86y//WHGRTXM6k0bpWKONjoulSLwKlbtwbFZVSm7XWrXItc3aom1LKE3gUGAo0wBihtR9jyO+7WmtL8cItOZJEhLhBJF+A+b3hjG3pI2WCAbOg8YA8n7L5yAX+t+YQq/aeyVFWv3Ig42Ki6dW0Cp4eRb5Hu9xyWRJRSs0BPtJab8yjvDUwTmv9YG7l7iZJRIgbyJXzML8XnLVdplUecM88aNg736f9fSaRj9YeYsm2E1iy3eFeLdiXMZ1qMeiWSHzNHiUUeNnjyiRiBe7XWi/Ko3wQsEhrXSbfXUkiQtxgks7BvLsg/i9j2+QJ9y6A+j0LfOqJiynM/vUwn206SnK2e0pC/c0MbxfFsHY1CfG//md+yi+JuPq8zB+QQddCiLIhIAKGfwdhdYxtqwUWD4O/Vxb41GrBvrx0d0N+f6YrT3SrS2iWZHHhSipvr/qb9q/9wuSluzlx8cYdHlzgmYhSKhKIsm2uAV4GVuVSNRR4AfDVWjd2XYiuI2ciQtygLp+EuT0h4bCxbfKE5vdDp39BcI1CNZGSms7i2GPMWHcoR9LwNCl6N6vKuJho6mYZ8XW9KFZ3llJqIjCR3O9Qd6gKWIGRWusyOcxXkogQN7CLx2BeT7h4NHOfyQtaPAAdnyx0MklLt7Jsxyk+XHuQfadz3lPSrUFFxsVE0yoq1FWRu11xk8jNQDOMJDEHmAH8ka2axlg6d5PW+lixIy4hkkSEuMFdPArfjIMjvzvuN3lBi2HQaQIEVcv9udlorVnz9zn+t+Ygfx6+kKP8lqgQxsVEc2u9iuV+HRNXXlifCHyltd7lquBKkyQRIQRaw+G1sHoKHNvgWOZhhhbDodOTUKHw0wBuPpLAh2sP8tOenMOD61YK4KHO0fRuVhWvcjo82KX3iZRnkkSEEHZaw6HVRjI5/qdjmYcZWo4wurkqVMn16bk5cNYYHvztthOkpeccHjyqYy3ua10DP7OzK5O7lyQRG0kiQogctIaDv8CaKXB8k2OZhze0Ggkdn4DAyoVu8tQlY3jwoj9zDg8O9vNieLsohrePchjxVZZJErGRJCKEyJPWcOBnWPMqnNjsWObpA60ehA6PQ2ClQjd5MTmVBX8cYd76OM5fSXUo8/XyYNAtNRjdqRbVQ/xccQQlRpKIjSQRIUSBtIb9PxnJ5ORWxzJPH2g1Cjo+DgEVC93k1bR0vog9xkfrDnE8wXF4sIdJ0fvmqjwUU5v6lcvmahqSRGwkiQghCk1r2L8SVr+aZTJHG09fuGWUcWYSEFHoJi3pVpbtPMWHaw+x99TlHOVd6xvDg2+JCilTEz5KErGRJCKEcJrW8PePxjWTU9sdy7z84JbR0OGf4B/uRJOatX+f48O1B9lwKOfw4BaRwYzvUofb6peN4cGSRGwkiQghikxr+Gu5kUxO73Qs8/KD1mOg/T/BP8ypZrceNYYHr9xzhuwfxzdVDOChmGh631wVs6f7hgdLErGRJCKEKDatYd8yWPNa5lTzGbz8oc1YaP8Y+Dl3x/rBc0nMWHuIr7cezzE8uEqQD6M61mJw60j8vUt/eLAkERtJIkIIl7FaYd/3RjI5u9uxzBwAbR6Cdo84nUxOX7rKnN8Ps3DDEa5kGx4c5OvF8HY1Gd4+irAA7+IeQaFJErGRJCKEcDmrFfYttSWT3/gAlgAADudJREFUPY5l5kBbMnnY6WRyKTmNTzYeYe7vh4lPchwe7ONl4t5WNRjTqTY1Qkt+eLAkERtJIkKIEmO1wt4lsOZ1OLfXscy7ArQZB+3+Ab4hTjV7NS2dLzcfZ8a6Qxy9kOxQ5mFS9GpahXEx0TSoUnLDgyWJ2EgSEUKUOKsV9nxjJJOMxbAyeFeAtv+AtuPBN9ipZi3pVn7YdZoP1x5k98mcw4O71ItgXEw0bWqFunx4sCQRG0kiQohSY02H3d/A2tch/m/HMu8g46yk7XjwCXKqWa01v+6P58O1B1l/8HyO8uaRwYyLiaZ7g0ouGx4sScRGkogQotRZ02HX10YyOb/fscwnyLj43mYc+DjfHbX92EU+XHuQH3efzjE8ODrCn4diounbrFqxhwdLErGRJCKEcBtrOuz6ypZMDjiW+QTbkslDRUomh84lMfPXQ3y1+QSp6VaHssoVbMOD20QSUMThwZJEbCSJCCHcLt0Cu740ksmFQ45lviGZycTb+WV2z16+yuzfD7Nww1GSrlkcyir4ePL0HfW5v21Np9vNL4mUzxVShBCivPLwhJvvg4c3Qd//QUitzLKUBPjlv/BOU/j1LbiW5FTTFSv48NydDfj92a48fUc9wrPcS3L5qgV/bw9XHYWdnIkIIYQ7pVtgx+ewbiokxDmW+YZCh8fgljHgHeB001fT0vl6ywk+WncQS7pmzVNdirS6onRn2UgSEUKUWelpsP0zWPcGXDziWOYXZkzyeMtoMPs737RVczwhmZphzj8XJInYSRIRQpR56Wmw/VNbMjnqWOYfYSSTVqPAXHoLWck1ESGEKC88vKDFMHhkM9z9LgTVyCy7cg5WvgjvNoX170Nqct7tlBJJIkIIURZ5mqHlCHh0C/R6GypUzyy7cg5WvgDvNYM/pkNaSp7NlDRJIkIIUZZ5mo313R/bAne9BRWqZZYlnYEVz8G7zWDDh5B2tdTDkyQihBDlgae3sSTvY1uh55sQWDWzLOk0/PiMcWay8aNSTSaSRIQQojzx9DZWUXxsK9z5BgRWySxLPAU/PA3vNYc/Z4LlWomHI0lECCHKIy8fYxXFx7bBHa9DQOXMssSTsPxfRjLZNKtEk4kkESGEKM+8fKDtOPjnNrjjNQiolFl2+QQsmwDvtYBNs8GSmnc7ReTWJKKUmqSU0vk80rLVr6eU+lYplaCUuqKU+lUp1dVd8QshRJnh5WtMLf/P7XD7q+BfMbPs8nH4//buPEiuqgzj8O+FLGxFCJAEARUwiITVyiICSihBWQpQERcUZQn7UiiFiFBsYqQCKQIISoAQoMBCi0Apq4RAZBUCxEDESAgJFsQw7InFCEk+/zi3pafTM0x67p3bM/M+VVM9c+7S30mn++t7zrnn3PWT1NSVs+5f8b2tacD8OuU7AqcDf6oUSPoM8BiwHJgAvAscDdwnad+ImF58uGZmTa7/2mk53pFHwKwp8OikNCQYUl9KzkpNIhExB5hTWy7p6uzX66qKfwVsAIyMiNnZfjcCc4ErJX0u+tLt92ZmHRmwDux6Eow6IjVlvb0Qhm2X+9OUfSWyCknrAN8FXgXuzcrWBQ4EHqokEICIWCbpWuACYDTwZPdHbGbWxAasmyZxLEgzdqx/G1gfuD4iVmRlOwIDgcfr7P9E9ji6G2IzM7MqzZhEjgICmFJVVrmr5tU6+1fKNquzzczMCtRUSUTSNsDuwIyIeLlqU2W6ynqDnVtr9qk95zGSZkma1dLSkl+wZmbWXEmEdBUCcG1NeWWqyoGsaq2afdqIiMkRMSoiRg0ZMiSHEM3MrKJpkoikfsAPgbeA22s2v5Y91muyqpTVa+oyM7MCNU0SAQ4AhgE3RURts9VzpKasL9Y5bpfs0atNmZl1s2ZKIpWmrOtqN0TEMtKNh2Ml7VQpl7QeMA54EQ/vNTPrdk2xPK6kTYFXgKcj4gvt7DOclCg+BC4F3iPdsb4DsH9E3NeJ52kBFn3cfu3YGHijwWN7Kte5b3Cd+4au1PnTEVG3U7lZbjY8HFiTVTvU/y8i5kvaDbgI+BkwAHgG2KezU56094/QGZJmtbfGcG/lOvcNrnPfUFSdmyKJRMR4YHwn9nsBOKj4iMzMrDOaqU/EzMx6GCeRzptcdgAlcJ37Bte5byikzk3RsW5mZj2Tr0TMzKxhTiJmZtYwJxEzM2uYk0gDJA2VdL2kOZLektQqab6k67KbInsdSZtJOlPSTEmLszXu50q6WNJGZcdXFEnHSrpZ0j8krZDU4zsRJa0h6cdZnVol/UvSxGzxt14p+7/7B0kLJIWkhWXHVCRJn5V0gaQnJLVIWipptqSz8n6d3bHegGzK+imkRbIWAe8DWwNHkmYa3iUi/l5ehPmTdBxwGXAX8AiwFBhDulF0MTAmIv5dWoAFyT5sNgKeBbYENo8IlRpUF0m6DDiFNNHpPcC2wMnAw8BeEbGyxPAKkSX/t0g3KI8E3ouILUoNqkCSLgJOBP5IWrjvQ2BP0qJ/c0ifUe/n8mQR4Z+cfkirKwZwVdmxFFC37YBN6pSPy+p8SdkxFlTvLYA1st/vTG+Z8uPq4uu4Eritpvzk7HU8tOwYC6r3VlW/Pw8sLDumgus7ChhUp/zC7HU+Ka/ncnNWvirzcg0uNYoCRMTcqH+lcWv2uH13xtNdImJh9K5v5t8DBEyqKb+GtCbPD7o9om4QEQvKjqE7RcSsiHi3zqbc369NMe1JTyWpPzAI6A8MB87LNt1dVkwl2Dx7XFJqFNZZo0lXIm1mvY6IVkmzs+3We+X+fnUS6Zqvkaaor1gCnBYRN5UUTxnOzx5vKDUK66xNgTdi1TV7IC3stqukARHxQTfHZQWTtCZwDrAcuCWv8/bpJCJpA+DU1Tjk8oh4q+rvJ4C9gbWBEcB3gMGS+kXE8vwizU8Oda4+12nAIcDkiJiRR3xFyLPOvcA6pAXe6mmt2sdJpPeZRFrE7+cRMS+vk/bp0VmStgBeXo1Dto6I+R2cb1PSyIfbIuLYrkVXjLzqLGkcaS6eu4FvRMSHuQRYgBzrfCdp7ZoeOzpL0nPA0IgYVmfb70lfCgb25isRSc8D60UvHp1VS9IvgLNJX/hy/Wzq01ciEbGQ1MmY1/lekzQdOErSKe00GZQqjzpLOpKUQP4MHNzMCQTyf517uNeAEZIG1vn/uRmpqavXJpC+SNJ5pARyPXBc3uf36Kz8rU1aYGv9sgMpgqQjSCN5pgNfb8ZEaR16ivS+H1NdKGktYGdgVhlBWTEknQucC9wIjIsCmp6cRBogaZWmgKx8BPAVYEFEtHRvVMWTdDhp9ckHgYMiorXjI6wJ3Uq6T6C2j+hoUl/Izd0ekRVC0jmkEaM3AUcUNVS9TzdndcGZkvYm3b29kNRUsj1wGGm47wnlhVYMSQcC15HWtr8VOFhq00K0LCLuKCO2Ikk6ANgp+3N4VnZ29vc7EfHrUgJrUEQ8J+lK4CRJ00h9WtuS7mCfSY6jdpqJpMOAT2d/DgEGVL2Oi3rbiEpJJ5JGTr5CajU4tOb9uiQi7s/lufpyx3qjJO0FHE+aPmEoqfnqVdKb8JKImFtieIXI2lXP7WCXRb2xo1LSVOBH7WzukXXOhnqeChxDuiP/DdIXg3MiYlmJoRVG0kPAHu1snhkRY7svmuJ9zP9byLHOTiJmZtYw94mYmVnDnETMzKxhTiJmZtYwJxEzM2uYk4iZmTXMScTMzBrmJGJmZg1zEjHrBEljJUU29UvpJB0v6T1JG63mcTtLWimpvRvvzFaLk4hZlexD9rxs+vimJGkQaUqLSyPizdU5NiJmA3cAE1UzD4ZZI5xEzNramTS9yxY15X8hzdDcDHMsnQBsADQ6b9ck0pQ9++UWkfVZTiJmnRARKyOiNSJWlBmHpDVIc17d04WZoh8mTRya+9oS1vc4iZhlskkmr8/+fDDrAwlJU+v1iVSXSTpB0jxJrZKek7R/ts8Oku7N+i/elHS5pP51nntrSTdJWizpA0kLJV0sad2aXceQrpLurnOOT0qaImmRpP9Kel3SY5LaTMSXrSlxH7CPpPUa/xcz81TwZtWmAZ8gfdMfD7yQlb8EDOzguBOBwaS1VlpJ06rfIekQ0gJevyP1Q3wVOBl4HbiwcrCkkcAM4B3gatKM0Dtl59lN0h5Vq0dWOsSfrA5AUj/gftLqhFcB/wQGATsCXwJuqIn5ceBYYHfg3g7qZtYhJxGzTETMkfQ4KYncHxEPVbZJGtvBoZsCIyLi3WzfGcDfSEnpWxExLdvvt5KeJiWdC6uOnwIsBkZHxNKq53wgO8f3galZ8Yjs8aWaGEYA2wBnRMSETlS3cvx2OIlYF7g5y6zrplYSCKRkRFq867WqBFLxCLBJpRlJ0g6kq4VbgIGSNq78ZPv+h3QFUzEEWB4R79Wct/L8e0oa2omYK6O6OrOvWbucRMy6bkGdsreBl9spB6jc37Ft9ng+0FLz8zqwLlC9HHMAqh2eGxGLgF+SEs5iSU9LmiBpdDsxV473gkLWJW7OMuu69kZsdTSSSzWPE2m/Wentqt9bSCtprs9HVx8ARMTZkqYA+5P6QcYBp0uaEBFn1Jxzw6rzmTXMScSsre7+Zv5i9rgiIqZ3Yv/ns8etgVm1GyNiAXAFcIWktUijsH4qaWJEvF616/Ca85k1xM1ZZm1V1hjfsMO98vMs6YP8OElb1W6U1E9SdSwPZY+71Ow3qHbocES08tEIs8E1p94FWA482njoZr4SMav1FLASOEvSYFLHdr2+jVxEREg6jDTEd07WHDUXWId0tfBN4Ew+Gp31NKkPZj/a3rG+JzBZ0m3APFIyHElq0vprRMyr7Jj1p+wL3BsRyzDrAicRsyoR8YqkI4EzgN8A/Un3WEwt8DlnS/o8KVkcSLqTfCnprvKpwANV+4akq4HxkoZFxJJsU2VI8VjSkOA1gVdI97tMrHnKLwOfIk2fYtYlSjevmllPIWl9Ul/KNRFxdgPH305KIqPCHwDWRU4iZj2QpOOACcCWqzOTr6SdgWeAPSNiZlHxWd/hJGJmZg3z6CwzM2uYk4iZmTXMScTMzBrmJGJmZg1zEjEzs4Y5iZiZWcOcRMzMrGH/A9E6K1ka5FXPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "#set up the function to find T(t)\n",
    "def funcTemp(tempInit,K,time): #function accepts 4 arguments: initial temperature (tempInit), the constant K value (K), and the elapsed time (time)\n",
    "    temp=0\n",
    "    if time == -3:\n",
    "        temp = 55\n",
    "    elif time == -2:\n",
    "        temp = 58\n",
    "    elif time == -1:\n",
    "        temp = 60\n",
    "    elif time == 0:\n",
    "        temp = 65\n",
    "    elif time == 1:\n",
    "        temp = 66\n",
    "    elif time == 2:\n",
    "        temp = 67\n",
    "    currentTemp = temp + (tempInit-temp)*(math.exp((-K)*time)) #use the arguments to solve for Temperature as a function of time (currentTemp)\n",
    "    if currentTemp > 98.6:\n",
    "        currentTemp = 98.6\n",
    "    return currentTemp\n",
    "\n",
    "#function for euler formula\n",
    "def funcEuler(initialTemp,initialTime,finalTime,K,stepSize): #function accepts 6 arguments: ambient temperature (ambientTemp), initial temperaturte(tempInit), the constant K value (K), and the elapsed time (time)\n",
    "    currentTime = initialTime\n",
    "    currentTemp = initialTemp\n",
    "    ambientTemp=65\n",
    "\n",
    "    while currentTime <= finalTime: #increment the current time by the step size until it reaches the desired time\n",
    "        if currentTime == -3:\n",
    "            ambientTemp = 55\n",
    "        elif currentTime == -2:\n",
    "            ambientTemp = 58\n",
    "        elif currentTime == -1:\n",
    "            ambientTemp = 60\n",
    "        elif currentTime == 0:\n",
    "            ambientTemp = 65\n",
    "        elif currentTime == 1:\n",
    "            ambientTemp = 66\n",
    "        elif currentTime == 2:\n",
    "            ambientTemp = 67\n",
    "        currentTime += stepSize #increment the current time by the step size until it reaches the desired time\n",
    "        currentTemp += stepSize * (-K*(currentTemp-ambientTemp)) #Add the previous value of the formula to the current value to get the total temperature\n",
    "    return currentTemp\n",
    "\n",
    "initialTemp=85\n",
    "K=0.611\n",
    "startTime = 0\n",
    "endTime = 10\n",
    "timeStep = .1\n",
    "timeArray = np.zeros(6)\n",
    "tempAnalytical = np.zeros(6)\n",
    "tempEuler = np.zeros(6)\n",
    "for t in range (0,6,1):\n",
    "    timeArray[t] = t-3\n",
    "    tempAnalytical[t] = funcTemp(initialTemp,K,t-3)\n",
    "    tempEuler[t] = funcEuler(initialTemp,startTime,t-3,K,timeStep)\n",
    "\n",
    "plt.rcParams.update({'font.size': 18})\n",
    "plt.rcParams['lines.linewidth'] = 3\n",
    "plt.xlabel('time(s)')\n",
    "plt.ylabel('temperature(F)')\n",
    "plt.plot(timeArray,tempAnalytical,label='analytical solution')\n",
    "plt.plot(timeArray,tempEuler,label='Euler approximation')\n",
    "plt.legend()\n",
    "print ('The function looks correct, and although I am sure there is a better way to find the solution, I have no idea what it is.')\n",
    "print ('The Euler approximation does not work with times lower than 0 and I can not figure out why.')"
   ]
  },
  {
   "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
}