Project_01.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 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": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.275\n"
     ]
    }
   ],
   "source": [
    "Temp_t2 = 74\n",
    "Temp_t1 = 85\n",
    "Temp_ambient = 65\n",
    "delta_t = 2\n",
    "K = ((Temp_t2-Temp_t1)/delta_t) / (Temp_t1-Temp_ambient) #Finds K via a rearrange of Newton's Cooling\n",
    "print(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": 3,
   "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 - First temperature measurement of specimen\n",
    "    Temp_t2 - Second temperature measurement of specimen\n",
    "    Temp_ambient - Temperature of the surroundings\n",
    "    delta_t - The change in time between the first and second temperature measurements \n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    K - The cooling constant associated with the enviornment established via the inputs above\n",
    "    \n",
    "    '''\n",
    "    K = ((Temp_t2-Temp_t1)/delta_t) / (Temp_t1-Temp_ambient) #Finds K via a rearrange of Newton's Cooling\n",
    "    return K"
   ]
  },
  {
   "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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeZzN9f7A8debwRiMsYwMY80WhsEgS6gkLUKLSwiVSeVqvV33UolUt1Q3kRJabhIV2gtFRfZ9j8bSMKaxL2Mw4/3745yZ3zHOzJwzc2bOLO/n43Eec77nfD/f7/s7w/mc72d5f0RVMcYYY7JSzN8BGGOMKRiswjDGGOMRqzCMMcZ4xCoMY4wxHrEKwxhjjEeswjDGGOMRqzCM8TMR6SIif4jIaRHpnofnDRQRFZFw5/b7IvJUXp0/u0RkmIgs8nccRZFVGMYt54dX6uOiiJx12e7v7/hyQkQOiUhHf8fhYjzwsqqWVdXv3e0gIoNEZJ2InBGROBH5WkSu9mUQqjpYVV/25TEBRKSRs2JK/fdzSETeEJHivj6XyV1WYRi3nB9eZVW1LLAf6OHy2kx/x5cREQkogOeoBWzN5Hz/Bl4CxgChzv2nAT19HEduSnH599QC6AoM9XNMxktWYZhsEZHiIvK0iMSIyGERmSkiIc73GolIsojcJyIHROSIiNwrIu1EZIuIHBeR11yONUxEfhKRd0TkpIhsE5FOLu9XFJEPnd9M/xSRZ0WkWLqyk0XkGDDSef4lInJURBJE5AMRKefc/1OgCrDA+W13hIh0F5Hd6a4v7S5ERF4SkY9FZLaInAL6Znb9Gfy+HnY2Ox0RkbkicoXz9VigWmo8bspVAp4BolX1S1VNVNXzqjpfVf/l3Ke08/rjRCRWRF4RkRJZndvNuT4RkdHO591FZLeI/Nv5OzzgemcpIlVE5Dvn32uF83fkUTORqsYBPwKNXY4XISK/Ov9tbBKRm9Kd61vnuZbjqDBT35suIuPTXcdCERnmSSzGO1ZhmOz6B9AN6AiEAxeA113eLw40A+oCQ4A3gSeBzs7Xh4hIW5f9OwEbgUo4vk3PF5Fg53szgRPOY7UBegED05XdAFQGXnW+NhaoCkQADYFRAKp6F/AX0M35jXeih9d7B/ABUB743IPrTyMiNwNPA72B6sBh4CNnPOGu8bgpfg2gwNeZxPYcjt9pBNAK6AI8ldW5PVALEBwV2nDgbRFJjXEqkABcAUQDgzw8JuLoM7kBWOHcDsRxffNx3EH9A/hUROq4nOuo81wPAve6HO4D4G4REeexqgEdgDmexmO8oKr2sEemD2Av0DXda3uADi7bdYBEHB8wjXB8yFVyef8M0NNl+xtgmPP5MGBPuuNvAu7C8aF1Bijh8t4Q4DuXsr9nEX9fYLnL9iGgo8t2d2B3ujJp++CowBZ4ev1uzj8TGOuyHQJcBKq6iydd2fuAvVlc3wHgOpftnsCOrM4NBDr/TuHO9z4BRrv8Tk4AxVzKngQineUuArVc3psALMogvtR/D8edDwWWAGWc798A7HP93QHzgJEu56rt8t5rqedy/nuLAa5xbj8JzPX3/5nC+rA7DOM157e5GsC3ziaE48B6HHeslZy7pajqEZdiZ4H4dNuu36hj051mH45vtrVwfGgkuJzrDRzfNlP9mS6+aiLyqbMZ5SSO9v7K2bhUV2nn8PD6XVVzXg8Aqnocx4dvdQ/OewS4IvUbdHrO16u6Ht/5PPXYOTl3gqpedNlOxPE3q4rjg9r1b3bJ38CNFFUNUdUQ5zE28f93TdWA/er8xE93Dann+jPde6nXo8CHwADnSwOA/2V9aSY7rMIwXnP+J039Vhvi8ghU1cPZPGx4uu2awEEcHxSngQou5wlW1ZauIaUr+wqOu5KmqhoM3I/jQyej/c8AQakbzvb/iun2SSuTjes/yKXt7uWBYOcxsrLUGfst7t50xnLI9fg4fnepx87JuTNyCMfvw7XSqeFpYVU9g6MpqYuzieugM2ZXqdeQeq4a6d5z9SFwp4i0cu73jaexGO9YhWGy623gJRGpAWkdkz1ycLwazg7sABEZgONDYYGq7sHR1v2yiJQTkWIiUl8yHxZbDkclc1JEagKPp3s/Hkd/SKrtQEURud5ZWTxH1v83vLn+WcBQEWnqbK//D/CTqh7K4hw4K6BxwDsicquzg7uEiPQQkRdcjv+siFQSkSo4+ms+cnkvW+fOJKYk4CvgOXHM5WgK3O1peWccA4B9qnoa+BUoJiKPOv/+N+DoH/o03blKi0gz4JJh3aoaA2wD3gNmq+r57F6byZxVGCa7XgYWAT85Rw79BrTMvEimfsEx3PIojg+83qp6wvlePxxt7zuc78/m0iap9J7B0Rl9Akdb+Ofp3h8PjHc2Jw13fig/gqO9PxbHt9qs7pQ8vn5V/Rp4EfgSx7fpqlzaaZ8pVR2P43cyzhnXfhwdzV+4XO82HENzNwDLnPHl+NyZeABHU1ICjia/WcC5TPYvLs55GEAc0BzH4IXUCuhW4E4cTXCvAX9T1T9cznUFjor+HRwVQ3of4Oj0t+aoXCSXNhsak/ecQyDvVNWu/o7FZI+IvAEEquoDfjp/N+AtVa3nj/MXFbk+yckYU/g4m6EUx51NO+AeHHeC/oilJDACx/Bbk4usScoYkx3lcfQtnMHRX/K8ZpDWJDeJSCRwDEe/1eS8Pn9RY01SxhhjPGJ3GMYYYzxSqPswKleurLVr1/Z3GMYYU6CsXbv2sKqGpn+9UFcYtWvXZs2aNf4OwxhjChQR2efudWuSMsYY4xGrMIwxxnjEKgxjjDEeKdR9GKbwunDhArGxsSQlJfk7FGMKrMDAQMLDwylRokTWO2MVhimgYmNjKVeuHLVr1yaDzN/GmEyoKkeOHCE2NpY6depkXYA8bpISkcdEZKs4lumc5cx0Oca5bsEG5+PmDMp2F5GdzmUjR+ZmnFuWb2BZeAjbVm7KzdOYHEhKSqJSpUpWWRiTTSJCpUqVvLpLz7MKQ0Sq48j3EqWqTXEs4dnX+fbrqhrpfHzrpmxxHNP+b8KxDnA/EWmcfj9f2fLQAK4+eIJNwzzO2Gz8wCoLY3LG2/9Ded3pHQCUFpEAHAvWHPSwXBscS2jGOHPdf4JjGUqfOltCQIR6CReZcnUf+m7YCiKO140xpojLswpDVQ/gWPd3P458+CdUdYHz7eEisklEZohIBTfFq3PpEo2xZLDEpIhEi8gaEVmTkJDgVYx//LKe+Y1q8UP91kzodA8rq9djXqPa7Fm60avjmKJBRHjiiSfStidMmMCYMWPyNIY1a9YwYsSIbJXt0qXLZRNbe/fuTWRkJPXq1aN8+fJERkYSGRnJb7/95otwc8VPP/3EihUr8uRcy5Yto23btkRGRnLVVVcxbty4PI8hvWnTphEaGpr2t3rvPXfLhfhGXjZJVcBxV1AHx8IrZZwrq00BrsSxuHwc8Kq74m5ec5s1UVWnqmqUqkaFhl42sz1TTdtFklSqDA8v/5RKZ44xofP9JAWWoXHbZl4dx+RPcXHQuTMcyvZac5cqVaoUc+fO5fDh7K5KmzPJyclERUUxceJEnx1z3rx5bNiwgWnTpnHNNdewYcMGNmzYQPv27X12juxITk7O8L3sfFhndrzMDBo0iOnTp7Nhwwa2bNnCHXfcke0YfKl///5pf6shQ4bk2nnyskmqK7BHVRNU9QIwF2ivqvGqmuJcbP5dHM1P6cVy6Zq+4XjenOWV4NPH+O6qunSvE8TqGk2JKVcr60KmQBg3DpYuhbFjfXO8gIAAoqOjef311y97b/DgwXz22Wdp22XLlgVgyZIldO7cmT59+tCgQQNGjhzJzJkzadOmDREREfzxh2ORuYSEBO644w5at25N69atWbZsGQBjxowhOjqabt26cc8997BkyRJuvfVWAE6fPs2QIUOIiIigWbNmfP65Y6HBBx98kKioKJo0acKzzz6b7etdvXo1nTt3plWrVtx0003Ex8cD0LFjRx5//HGuueYaGjduzJo1a+jduzf169dPu+PavXs3TZo0YeDAgURERNCnTx/Onj2b5XFHjRpFp06dmDRpEl988QVt27alRYsWdOvWjb/++os//viDadOm8corr6TdCQ0YMID58+df9rtftGgRXbt2pW/fvrRo0QKADz74gDZt2hAZGclDDz3ExYsXM/0dJCQkULVqVQCKFy9O48aN3cYQHx/P7bffTlRUFG3atEmrTEaPHs2gQYO49tprqV+/PjNmzMj238MvVDVPHkBbHEtIBuG4Y/gA+DsQ5rLPY8AnbsoGADE47k5KAhuBJlmds1WrVppdF5JTtNtrP2unl3/SpAvJ2T6OyR3btm3zeN/AQFW4/BEYmLMYypQpoydOnNBatWrp8ePH9ZVXXtFnn31WVVUHDRqkn3766SX7qqouXrxYy5cvrwcPHtSkpCStVq2aPvPMM6qq+t///lcfeeQRVVXt16+f/vrrr6qqum/fPm3UqJGqqj777LPasmVLTUxMTDveLbfcoqqqTz31VFp5VdWjR4+qquqRI0dUVTU5OVk7d+6sGzduVFXVzp076+rVq91em+txVVWTkpK0Xbt2mpCQoKqqH330kQ4dOlRVVTt06KD//ve/VVV1woQJWr16dT106JCePXtWw8LC9NixY7pr1y4FdPny5aqqOnDgQH399dezPO7w4cMvuZ6LFy+qquqUKVP0qaeeUlXVUaNG6euvv562X//+/XXevHmX/e4XLlyoZcqU0X379qmq6ubNm7Vnz5564cIFVVUdOnSozpw50+3vI9UzzzyjISEh2rt3b506daomJSW5jaFPnz5p17pnzx5t0qRJ2n4tWrTQs2fPanx8fNrvKr127dpp8+bNL3v89NNPl+377rvvalhYmEZEROhdd92lsbGxmV5Deu7+LwFr1M1nap7Nw1DVlSLyGbAOSAbW41gha5pzERQF9uJYvxcRqQZMU9WbVTVZRIYDP+AYXTVDVbfmZrwBxYsx6paruGfGKv63fB/3X1M3N09nclFMDDz5JMyfD4mJEBQEvXvDhAk5P3ZwcDD33HMPEydOpHTp0h6Vad26NWFhYQBceeWVdOvWDYCIiAgWL14MOL4Nb9u2La3MyZMnOXXqFAC33Xab23MtWrSITz75JG27QgVHd+CcOXOYOnUqycnJxMXFsW3bNpo1866Zdfv27WzdupWuXR2r6KakpBAeHp72/m233ZZ2DREREVxxhWPJ9dq1axMbG0tgYCB16tTh6quvBmDAgAFMnTqVLl26ZHrcvn37pj3fv38/ffr04dChQ5w7d44GDRp4dQ0A7dq1o2bNmoDj97V69WqioqIAOHv2LDVq1MisOM899xwDBw5kwYIFfPjhh8yePZtFixZdtt+iRYvYuXNn2vaxY8fS7qh69epFYGAggYGBdOrUidWrV6fdJabyps+oV69eDBw4kFKlSjF58mSGDBnCggULsi6YDXk6cU9VnwXS3xO7XZBeVQ8CN7tsfwtcNuQ2N3VqEEqXhqG88eMubm8ZTsUyJfPy9MZHwsIgOBiSkiAw0PEzOBicLQs59uijj9KyZctL2o4DAgLSmjdUlfPnz6e9V6pUqbTnxYoVS9suVqxYWtv6xYsXWb58uduKoUyZMm7jUNXLhknu2bOHCRMmsHr1aipUqMDgwYOzNTteVWnWrBm//vqr2/ddryH99aVeU/rYRCTL47pe68MPP8y///1vbr75ZhYtWsRLL73ktozr7z4lJeWS/grX46kq9957b1rHtafq1atHvXr1GDp0KJUqVeLEiROX7aOqrFq1ipIlL//McPd7SK99+/YkJiZe9vrrr7/Otddee8lrlStXTnseHR3N6NGjPb4Wb1kuqSyMuvkqEs+n8NrCnVnvbPKt+HgYNgxWrHD89FXHN0DFihXp06cP06dPT3utdu3arF27FoAvvviCCxcueHXMbt26MWnSpLTtDRs2eF3m2LFjnDx5kjJlylC+fHni4+P57rvvvIojVePGjTlw4ACrVq0C4Pz582zd6t1N/p49e1i9ejUAs2bNomPHjl4d98SJE1SvXh1V5YMPPkh7vVy5cml3X3Dp737evHmkpKS4PV7Xrl2ZM2dO2qCFI0eOsH//fsDRibxu3brLynzzzTepzeT8/vvvlCpVinLlyl0WQ9euXZk8+f9XjHX9+82fP59z585x+PBhfv3117Q7HFe//fZbWie26yN9ZQEQFxd3ybGbNGni9np9wSqMLNS/ohwDr67Fxyv3sz3upL/DMdk0dy5MngzNmzt+zp3r2+M/8cQTl4yWGjp0KD///DNt2rRh5cqVGd4VZGTixImsWbOGZs2a0bhxY95+++0sy4wePZpjx47RtGlTmjdvzuLFi2nevDktWrSgSZMm3HvvvXTo0MHrawPHHcRnn33G448/nnbMlStXenWMJk2a8O6779KsWTPOnDlDdHS0V8cdM2YMvXv3pnPnzmlNXgA9e/Zkzpw5tGjRgt9++40HHniAhQsX0qZNGzZs2HDJHY+riIgInn32Wbp27UqzZs3o1q1bWof7pk2b0jq3Xb3//vs0bNiQyMhIBg8ezMcff0yxYsUui2Hy5MksW7Ys7e/37rvvph2jdevW3HTTTbRr147nnnvukmvJjtdeey3tbz5lypRLvrj4WqFe0zsqKkp9sYDSicQLdJmwmEZVg/l4aFubYZwPbN++nauuusrfYRgP7d69mzvvvNOjOyV/O3bsGA8++OAl/UG+Mnr0aCpXrsyjjz7q82Nnl7v/SyKyVlUvu/WxOwwPlA8qwePdGrI85gg/bPVhW4YxJt+pUKFCrlQWhYFlq/VQv9Y1mLliH89/s50uDasQWKK4v0MypsCoV69egbi7yG3PP/+8v0PIEbvD8FBA8WI806MxscfO8s7PMf4Oxxhj8pxVGF5of2Vlbm0WxltLdvPn0cuHvBljTGFmFYaXRt1yFcVEGPf1tqx3NsaYQsQqDC+FlS/N36+vx4Jt8SzZ+Ze/wzHGmDxjFUY23NexDnUql+G5r7ZxLtn9pCBT+Fl68/yhqKc3X7x4MS1atCAgIOCSpIsA06dPp379+tSvX5+PPvoox+eyCiMbSgUU59kejdlz+AzTvtnk25zZJtfEnYqj8/udOXTaN38rS2+edyy9ecZq167Nhx9+SJ8+fS55/fDhw7zwwgusXr2aFStW8PTTT7tNY+INqzCyqUvDKtzUtCoTf/uTPzfv8l3ObJNrxv0yjqX7lzL2Z9/8rSy9uaU3zw/pzevUqUNERATFil36cf7dd9/RvXt3QkJCqFSpEtddd13OkxK6S2FbWB45SW+epcBAPVCusl712Kc65I5n9KKvcmYbj3iV3vz5QGUMlz0Cn8/Z38rSm1t68/yQ3jyj637xxRf1xRdfvCR21xhT5cv05oVOTAzVnnySx1Z9yviO97CgaWdubB7um5zZxqdiRsTw5IInmb9jPonJiQQFBNH7qt5M6Jbzv5WlNyctptRrsPTmeZvePCPqJntxTtMaWYWRXc6c2YNXTuezRp157prBdGQtZXyVM9v4TFi5MIJLBZOUkkRgQCBJKUkElwqmalnf/K0svbmlN09P8zC9eUbCw8Mv6VeJjY2ladOmHpXNiPVh5ER8PCWihzK+bxQHg0P5b3J1f0dkMhB/Jp5hrYax4r4VDGs1zGcd32DpzT1l6c0dfJ3ePCPdu3fnu+++4/jx4xw5coQff/wx7W42u/K0whCRx0Rkq4hsEZFZIhIoIq+IyA4R2SQi80QkJIOye0Vks4hsEJGcp6D1BWfO7KiubejXpgYzKjdjy4GcjUIwuWPu3+Yy+ZbJNK/anMm3TGbu33yb39zSm2fN0ps7+Dq9+fLlywkPD2fevHncf//9ac2NoaGh/Otf/yIqKoq2bdsyduxYypcvn6Nz5Vl6cxGpDiwFGqvqWRGZg2MFvYPAT+pYhvU/AKr6Tzfl9wJRqurxGEZfpTf3xInEC1z/2s9UCwlk3kMdKF7MUqDnJktvXrBYenMHS2/unQCgtIgEAEHAQVVdoKqpjYwrgPAMS+dj5YNK8EyPxmyKPcEHv+31dzjGmGyy9OYZy7NOb1U9ICITgP3AWWCBqqYfFHwvMDujQwALRESBd1R1qrudRCQaiAbSRkPklR7Nwvh8bSyvLthJ96ZVqRbi2agZYwo7S2/uYOnNPSQiFYCeQB2gGlBGRAa4vD8KSAZmZnCIDqraErgJeFhEOrnbSVWnqmqUqkaFhob69BqyIiI836spFxVGz99CXjX3GWNMXsjLJqmuwB5VTVDVC8BcoD2AiAwCbgX6awafsqp60PnzL2Ae0CZPovZSjYpBPNGtAT/t+IsvNx70dzjGGOMzeVlh7AeuFpEgcQw8vh7YLiLdgX8Ct6mq20UmRKSMiJRLfQ50A7bkUdxeG9KhDs1rhPDcV9s4euZ81gWMMaYAyLMKQ1VXAp8B64DNznNPBSYB5YCFziGzbwOISDUR+dZZ/ApgqYhsBFYB36jq93kVu7eKFxNevqMZp5Iu2LoZxphCI09HSanqs6raSFWbqupAVT2nqvVUtYaqRjofw5z7HlTVm53PY1S1ufPRRFXH52Xc2dGwajke7FKPeesPsHiHrZtRGFl68/yhqKc3nzRpEhEREURGRnLNNdewY8eOtPd8nd7c7wkCc/ORq8kHPZB0IVlveG2Jth2/SI8nnvdrLIWNN8kH0xw8qNqpk2pcnE9iKFWqlNauXTstcZ5r8sG8kJo0L7u8ST7ob5lda/rEfzk9XmauvPJK3bx5s6o6kjlu3bo12zH4yokTJ9Kef/7552l/t4SEBK1bt64eO3ZMDx8+rLVr19bjx49fVt6b5IOWGiQXlQoozoS7mpNw+hzjv7GmKb8bNw6WLvVZKnpLb27pzfNDevPg4OC052fOnEnLTWXpzQvYHUaq/3y3XWv982tdfNsgn327Leq8usMIDFRNTT/v+shhKnpLb27pzfNLevM33nhD69atqzVq1NDdu3erqqU3L7Ae6VqfhT9u4F/Vu/DD2BcIfst3K6QZD8TEwJNPwvz5kJgIQUHQu7dPUtFbenPSYkq9BktvnvfpzUeMGMGIESP48MMPeeGFF5g+fTqqlt684CldmlJJSbwS1oDbB7zCuC3FeEUEAgPB+Q/I5DJnKnqSkhy/96Qkx7aPUtFbenNLb56eqn/Sm99999088sgjTJ8+3dKbF0gxMXD33USeiGXYys/5tNkNLLrvH7Bnj78jK1ri42HYMFixwvHTh2uwW3pzz1h6cwdfpzfftWtX2vOvvvqKhg0bAoUgvXmR5PLt9pG1c2n01x5GVmzL0XIV/R1Z0eJMRU/z5o6fcy29uTcsvblDfkxv/t///pcmTZoQGRnJpEmTeO+994ACnt7cH/IyvXmmbr/dUXFER7N92ixuC2pPt2bhTLq7RY7bFIsqS29esFh6c4eCnt7c+jDygsu32avebM6ji3fzyg876bbxCnpG2ip9xuQnlt48Y1Zh+MEDnery4/Z4Rs/fQlTtilS3NOimkLP05g6W3tx4LaB4MV7/WyQXLypPztnIxYuFt1nQGFN4WIXhJ7UqleHZHk1YHnOE6UttxJQxJv+zCsOP7ooKp1vjK3jlh51sO3jS3+EYY0ymrMLwIxHhpTuaUT6oBI98sp6z592PFzfGmPzAKgw/q1imJK/e1Zxdf53meUtQWODMmzcPEbkkpXR2pE9W6M4LL7xwyXb79u2zda4xY8YwIV1alPHjx6elMi9evHja84kT828am5iYGBvNlMeswsgHOjUIJbpTXWau3M/3W3w3A9nkvtQZy3nxwZW+wvDlGhWjRo1Km01cunTptOfZXWvDV1zTeqSX3Qojo5nfJmt5WmGIyGMislVEtojILBEJFJGKIrJQRHY5f1bIoGx3EdkpIrtFZGRexp0XnuzWkIjq5Rk5dxNxJyzHVEFw+vRpli1bxvTp0y/54FqyZAldunThzjvvpFGjRvTv3z8tncTYsWNp3bo1TZs2JTo6mvQTZ3/88Ud69+6dtr1w4UJuv/12Ro4cydmzZ4mMjKR///7A/6ftBnj55ZeJiIigefPmjBzp+O/x7rvv0rp1a5o3b84dd9zhNjeRJzJL1T148GC6detG7dq1mT9/Pk888QRNmzbllltuSfuwDw8PZ+TIkbRp04a2bdsSExOT5XEfeOABbrjhBoYMGcIff/zBNddcQ4sWLWjVqlXaTPCRI0eyePHitDuhadOmXTIhrnv37ixdupTk5GRCQkIYPXo0bdq0YdWqVRmmVDeZy7N5GCJSHRgBNFbVsyIyB+gLNAZ+VNWXnBXBSBxrfLuWLQ5MBm4AYoHVIvKlqhaaNpySAcWY2K8Ft0z8lUc+2cDH97cloLjdAHriua+2+nzQQONqwTzbo0mm+8yfP5/u3bvToEEDKlasyLp162jZsiUA69evZ+vWrVSrVo0OHTqwbNkyOnbsyPDhw3nmmWcAGDhwIF9//TU9evRIO+Z1113Hww8/TEJCAqGhobz33nsMGTKEHj16MGnSJLdzGb777jvmz5/PypUrCQoK4ujRowDcfvvtDB06FHB8CE+fPp2///3vXv8uRowYwVNPPcXVV1/N3r17ufXWW9myZQvgyA/1448/snHjRq655hq++OILXn31VXr06MH333+floW1QoUKrFq1ihkzZvD4448zf/78TI+7fv16fvnlFwIDA0lMTGThwoUEBgayY8cOBg0axMqVK3nppZeYNGlS2toX06ZNy/AaTpw4QcuWLXn++ec5d+4c1157LV9++SWVK1dm5syZPP3000ydOtXr301Rk9cT9wKA0iJyAQgCDgL/Aro43/8AWEK6CgNoA+xW1RgAEfkE6AkUmgoDoE7lMozr2ZQnPt3IxHtG8/irI3yWUdX43qxZs9K+0fbt25dZs2alVRht2rRJS9MdGRnJ3r176dixI4sXL+bll18mMTGRo0eP0qRJk0sqDBFh4MCBfPTRRwwZMoTly5fz4YcfZhrHokWLGDJkCEFBQYAjGSLAli1bGD16NMePH+f06dPceOON2brOzFJ133zzzQQEBBAREQHADTfcADjyNO3duzetTNI77rsAACAASURBVL9+/QBHUr/UO6DMjtuzZ08CAwMBOHfuHMOHD2fjxo0EBASkLTLljZIlS6bduWWVqt1kLM8qDFU9ICITgP3AWWCBqi4QkStUNc65T5yIVHFTvDrwp8t2LNDW3XlEJBqIBtLy3hckd7QK57eZX/NmeDuuHj+Z9m96l3q5KMrqTiA3HDlyhJ9++oktW7YgIqSkpCAivPzyy8ClKcyLFy9OcnIySUlJPPTQQ6xZs4YaNWowZswYt6nGU+8oAgMDueuuuwgIyPy/qbu05uDoSJ8/fz7Nmzfn/fffZ8mSJdm61sxSdbumNXd93zWtObhP4Z3ZcV2TNb766qvUqFGDjz76iAsXLlzSFOfKNa05cMnvtnTp0mkxZJVS3WQsz9o8nH0TPYE6QDWgjIgM8LS4m9fcTo9W1amqGqWqUaGhodkL1l9KlwYRxk56jLpHD/CINuBwmRDH6yZf+eyzz7jnnnvYt28fe/fu5c8//6ROnTosXbo0wzKpH2CVK1fm9OnTGY6KqlatGtWqVeP5559n8ODBaa+XKFHCbZr0bt26MWPGjLQ+itQmqVOnThEWFsaFCxeYOXNmdi8101Tdnpo9ezbguCtLzZjr6XFPnDhBWFgYIsIHH3yQ1u/jLq35+vXrUVX27t2bluI8PV+kai+q8rKRvCuwR1UTVPUCMBdoD8SLSBiA8+dfbsrGAq5LYYXjaM4qXJxrZ5QpUYxJX/yHk4FleWzoBFL+iPF3ZCadWbNmXdI5DXDHHXfw8ccfZ1gmJCSEoUOHEhERQa9evWjdunWG+/bv358aNWrQuHHjtNeio6Np1qxZWqd3qu7du3PbbbcRFRVFZGRk2pDZcePG0bZtW2644QYaNWqUncsEyDRVt6cSExNp06YNU6ZM4dVXX/XquMOHD2fatGlcffXV7Nu3L+2upkWLFqSkpNC8eXMmTpxI586dqV69OhEREYwcOZLIyEi3x/NFqvaiKs/Sm4tIW2AG0BpHk9T7wBqgJnDEpdO7oqo+la5sAPA7cD1wAFgN3K2qmX4tyDfpzb3x4IMwdSqULMknDTszsvvfebRrfR7t6v1ylIVZYU9vPnz4cFq0aMF9993n71ByLDw8nC1bthASEuLvUIwb3qQ3z/IOwznsNatHlv8SVHUl8BmwDtjsPPdU4CXgBhHZhWMU1EvO81YTkW+dZZOB4cAPwHZgTlaVRYHlsjLc3zpcye0ndvHGj7v4dVeCvyMzeaRVq1Zs2rSJAQM8bbE1Jm9keYchIkk4mn8yW+mnuKrmux7mAnmHkU7i+WR6TV7G4dPn+WZER8LKW38GFP47DGPyik/vMIDtqlpXVetk9ACO+Ch2k05QyQDe6t+KcxdSGP7xes4nX8y6UBFRmFeLNCYvePt/yJMKo52P9jHZVK9KWf5zZzPW7jvGeMs3BUBgYCBHjhyxSsOYbFJVjhw5kjbfxROezMN4BxiUxYkvH0xucizuVBx9P+/L7Dtnc2uzamzYf5xpS/cQWTOE3i2K9kSj8PBwYmNjSUiwvh1jsiswMNCrSYueVBjNUp+IyAJV7ZadwIz3xv0yjqX7lzL257G8dctbjLypEZsPnOBfczfT8IpgGlcL9neIflOiRAnq1Knj7zCMKVI86fRep6otnc/Xq2qLPInMBwpqp3fp8aVJSr78pi0wIJD9I47T482llAwoxpfDOxASdPksWWOMyYmcdHpXFZHBItKCzEdKGR+JGRHD3U3vJijAkRsoKCCI/hH92fPIHkLLleKtAS05dCKJv89aT3KKdYIbY/KGJxXGGCAK+C8QLiKbReQTEXlaRO7I1eiKqLByYQSXCiYpJYnAgECSUpIILhVM1bKORIQta1bg+V5N+XXXYf7zfc4W7jHGGE9l2Yehqpfk/BWRcBz9GhFAL+Dz3AmtaIs/E8+wVsOIbhXN1LVTiTsdd8n7fVrXYOvuON79dQ9Nyii9uuR9Aj5jTNGSZ6lB/KGg9mF46sJDD9P/RE02hF/FnOGdiKxhqReMMTmXk9Qg40TkUxF5X0Qa5k54xivOrLYlprzFlHkvUOVEAtEvfcmhytX9HZkxphDzpA8jRFXvwrHGhH8X+DUOzqy2BAVR6exJpn37CmfKBBM96iPOnrf1io0xucOTCuO8c4SUAmWy2tnkgbAwCA6GpCQIDKTRgV28kbSBzX8l8o/PNtrsZ2NMrvCkwhiFYy2LqcDs3A3HeMwlqy3DhtE1dhP/7N6IrzfF8fqiXf6OzhhTCHkySioJeCUPYjHemDv3/587Vy17QJWYhNNM/HEXtSsFcXvLop0+xBjjW550eq/zxT4m94kIz/eKoF3dSvzz802sjLEkwsYY3/GkSeoqEdmUyWMzUDm3AzWeKRlQjLcHtKJGxSAe+GgtMQmn/R2SMaaQ8CT5oCeLAWc5NMc5JNe1D6Qu8AyO1Oipw3VDgOOqetlivCKyFzjlPFeyuzHCxqF8UAneG9ya29/6jcHvrWbuQ+2pXLaUv8MyxhRwfpm4JyLFcazN3VZV97m8/ipwQlXHuimzF4hS1cOenqewT9zLyvr9x+j37goaVg1m1tC2BJX05PuBMaaoy0nywdxwPfBHuspCgD7ALD/FVOi0qFmBiX1bsCn2OCNmbSDlog23NcZkn78qjL5cXjFcA8SrakZjQhVYICJrRSQ6owOLSLSIrBGRNba4DnRrUpUxPZqwaHs8T3+xxeZoGGOyzeMKQxwGiMgzzu2aItLG2xOKSEngNuDTdG/1I/O7iw7OdTluAh4WkU7udlLVqaoapapRoaGh3oZXKA1qX5sHo67g45X7mfjFen+HY4wpoLy5w3gLRwd1P+f2KWByNs55E7BOVeNTXxCRAOB2MpkYqKoHnT//AuYBXldWRdlTP83gji0/8vqKOD5eud/f4RhjCiBvekHbqmpLEVkPoKrHnHcL3nJ3J9EV2KGqse4KiEgZoJiqnnI+7wZc1jFu3ChdGpKSEOClYsU5WjqY0Z+nUPG+e+i+ZYm/ozPGFCDe3GFccI5uUgARCQW8Wu5NRIKAG4C56d66rE9DRKqJyLfOzSuApSKyEVgFfKOq33tz7iLLJVFhiYspTF4wkeYXjjGi51Ms2+3xgDNjjPGqwpiIoymoioiMB5YCL3hzMlVNVNVKqnoi3euDVfXtdK8dVNWbnc9jVLW589FEVcd7c94iLV2iwqDTJ3jv3DrqVC5L9Idr2PjncX9HaIwpIDyqMJxDXn8BngJeBOKAXqqavuPa5EfpEhWGHPqTD+9rQ8WyJRn83ip2xZ/yd4TGmALA44l7zokcrXI5Hp8q6hP3srLvyBnufHs5Asx5oB21K1v2emOMbyburRCR1j6MyfhZrUplmHl/Wy6kXKT/tJUcOH7W3yEZY/IxbyqMa4HlIvJHatJBEdmUW4GZvNHginL87762nDx7gQHTVvLXqSR/h2SMyae8qTBuAq4ErgN6ALc6f5oCrmn18rx/b2viTybR/92VHD59zt8hGWPyIY8rDFXd5+6Rm8GZvNOqVkVmDG7Nn8cSGTBtJUfPnPd3SMaYfMbjiXupKUHSc5dZ1hRMV9etxPRb63Lv59sZ8PZSPn6wIyFB2ZmbaYwpjLxpkjrj8kjB0URVOxdiMn7U4YM3mPr58+z+6zT9p63keKLdaRhjHLK9HoaIlAK+VNUbfRuS79iwWi84U4ikWlKnJdG3j6be0QPMfDOaCmXsTsOYoiI31sMIwrFqnikMXFKIAHSJ38G7x5axO6wud1ufhjEG79Kbb3ZZx3srsBNHuhBTGKRLIUJSEp1LJTJtUGtiEk7Tb+oKEk7Z6CljijJv7jBSh9H2wJEttpqqvpkrURn/SJdChEOH6NQglBmDW7P/aCJ/m7qcQydsnoYxRZU3qUH+o6r/zOq1/MT6MHxn1Z6j3Pv+aiqWKcnHQ9sSXiHI3yEZY3KJL/owbnDz2k3ZD8kUJG3qVOR/97XheOJ5+ry9nD8STvs7JGNMHsuywhCRB0VkM9DQpQ9jk4jsATbnfogmv2hRswKfRLfjfMpF+ry9nK0HT2RdyBhTaHhyh/Exjn6LL/n/PoweQCtV7Z+LsZl8qHG1YOY80I5SAcXoO3UFa/Ye9XdIxpg8kmWFoaonVHWvqvYDTuJY/a4W0FREOnl6IhFpKCIbXB4nReRRERkjIgdcXr85g/LdRWSniOwWkZGentf4Xt3Qsnz6YHtCA4szYPIv/LR8p79DMsbkAW+G1d6PYxGlH4DnnD/HeFpeVXeqaqSqRgKtgEQcK/gBvJ76nqp+m76sc2nYyTj6TBoD/USksafnNr5XPaQ0c/78mvoJexk6/3c+X+t2OXZjTCHiTaf3I0BrYJ+qXgu0ABKyed7rgT+8SF7YBtjtXKr1PPAJ0DOb5zY5Vbo0iFD5rTeYNevfXL1vE098upF32/fxd2TGmFzkTYWRpKpJ4EgLoqo7gIbZPG9fYJbL9nBnR/oMEangZv/qwJ8u27HO1y4jItEiskZE1iQkZLc+M5lymRVe9vxZZnzzMrck7md8p0GM/WobFy9mL92MMSZ/86bCiBWREGA+sFBEvgAOentCESkJ3Aakrgc+Bcc6G5E41gp/1V0xN6+5/VRS1amqGqWqUaGhod6GZ9yIOxVH5/c7c+j0IccL6WaFl0o8zZvJmxnSoTYzlu3h75+sJ+lCin+DNsb4nEcVhogIMEJVj6vqGOBpYDrQKxvnvAlYp6rxAKoar6opqnoReBdH81N6sUANl+1wslFZmewZ98s4lu5fytifXTLZp5sVXuzQIZ65tTGjbr6KbzbFcc+MVZbp1phCxpuZ3mtVtVWOTyjyCfCDqr7n3A5T1Tjn88eAtqraN12ZAOB3HH0fB4DVwN2qujWzc9lM75wpPb40ScmXpwIJDAjk7KiM1//+cuNBnpyzkfCKpXlvcGtqVSqTm2EaY3zMFzO9V4hI6xwGEYRjxvhcl5dfdlkf/FrgMee+1UTkWwBVTQaG4xiZtR2Yk1VlYXIuZkQMdze9m6AARxqQoIAg+kf0Z88jezItd1vzanx0f1uOnjlP77d+Y+2+Y3kRrjEml3lTYVyLo9L4w9lBnfoh7zFVTVTVSqp6wuW1gaoaoarNVPW21LsNVT2oqje77PetqjZQ1StVdbw35zXZE1YujOBSwSSlJBEYEEhSShLBpYKpWrZqlmXb1KnI3AfbExwYQL93V/DFhgN5ELExJjd5vEQrljeqSIo/E8+wVsOIbhXN1LVTiTsd53HZuqFlmftQB4b9by2PfLKBPxLO8Oj19SlWzN0YBmNMfudNH4YA/YG6qjpWRGoCVVV1VW4GmBPWh5E/nPvzAKOe/R+fVYnglmZhTLizOaVLFvd3WMaYDPiiD+MtoB3Qz7l9Csfsa2MyVerF8bzywShGJm3n281x3PXObxw8nnGnuTEmf/Kmwmirqg8DSQCqegywhZ5NxpwzwpkyBbl4kWFv/INpnz7H3phD3DZpqSUuNKaA8abCuODM6aQAIhIKXMyVqEzhkG6dcIKCuP7qBswfGkXZUo7O8Jkr9+Fps6gxxr+8qTAm4kgWeIWIjAeWAi/kSlSmcHCzTjjBwdRrXIcvHu5Ih3qVGTVvCyM/32wzw40pADyuMFR1JvAUjkriINBLVT/NvJQp8tysEw5QPqgE0we1Zvi19Zi95k/+9s5yDli/hjH5mjejpAKBh4BrcDRFLQWmpCYkzI9slFTB8P2WQzz56UZKFBfe6NuCTg0sB5gx/uSLUVIfAk1wNE1NAq4C/ueb8ExR1r1pVb4c3oEq5QIZ9N4qJv64yzLeGpMPeVNhNFTV+1R1sfMRDTTIrcBM0VI3tCzzHm5PzwYVeG3h7wx6eylHTp/zd1jGGBfeVBjrReTq1A0RaQss831IpqgKKhnA66s/YvyCyazcd4ybJ/7Kqj029NaY/MKb1CBtgXtEZL9zuyawXUQ2A6qqzXwenSk6SpeGpCRS0wlEHtjBwz3/Rb8piTx641U8dG09iltKEWP8yps7jO5AHaCz81EHuBm4Fejh+9BMkZJuzkaT0/F8dXYZNzeqzKsLf2fAtJXEn8y34yuMKRK8qTD24xghNci5FrcCV6jqPi/W5jbGPTdzNsqVC2LikHa8fGczNvx5nO7//YVF2+L9HakxRZblkjL5h5s5GyJCn6gafPX3jlQtX5r7P1zD6PmbOXveJvoZk9e8mYexTlVbish6VW3hfG2jqjbP1QhzwOZhFC7nklOY8MNO3v11D1eGluGNvi1oWr28v8MyptDxxTyMHOWSEpGGIrLB5XFSRB4VkVdEZIdzUaZ5IhKSQfm9zkWbNoiI1QJFUKmA4oy6pTEf3deW0+eS6TV5GZN+2kVyiqU0MyYvZCeXVJXs5JJS1Z2qGqmqkUArINF5vIVAU+coq9+Bf2VymGudx7is5jNFR8f6lfmhX0O6/7WNCQt+p887y9lz+Iy/wzKm0MtOLqkXgThylkvqeuAPZ4f5Auea3QArgPBsHtMUISGvvMik9/7JG4nr2f3XaW564xfeW7bHZogbk4s87sPw6UlFZgDrVHVSute/Amar6kduyuwBjuFoEntHVadmcOxoIBqgZs2arfbtswFchYpzvoarQ2UrMfKWR1hSuyVt61TklTubU7NSkJ8CNKbgy3EfhohEOfsY1jn7GzaLyKZsBFISuA34NN3ro4BkYGYGRTuoaksca4s/LCKd3O2kqlNVNUpVo0JDLYldoeNmjY2qPbvz3uv38/Idzdh28CQ3/vcXZizdQ4rdbRjjU97M9J4J/APYTM4WTroJx91F2oB6ERmEYwLg9ZrBLY+qHnT+/EtE5gFtgF9yEIcpiDJYY0PCwugT5ujfGDVvM2O/3sY3m+P4zx3NqFelrL+jNqZQ8KbTO0FVv1TVPamT9bI5Ya8fMCt1Q0S6A/8EblPVRHcFRKSMiJRLfQ50A7Zk49ymMMhgjQ2AaiGlmTG4Na/1ac7uv05z8xu/8saiXZxPtpFUxuSUN/MwrsfxYf8jkJZGVFXnenwykSDgT6Cuqp5wvrYbKAUcce62QlWHiUg1YJqq3iwidXGMqALHXdHHqjo+q/PZPIyiLeHUOZ77aitfb4qjfpWyvHh7BFG1K/o7LGPyvYz6MLypMD4CGgFb+f8mKVXVe30WpY9ZhVG4xJ2Ko+/nfZl952yqlq3qcbmfftvB6E/Xc7BUMP3a1OCf3RsRElQyFyM1pmDzxcS95s7O5EGqOsT5yLeVhSl8xv0yjqX7lzL257Felbvuo4ksnHI/Q8/FMGdNLNe/+jOfr43FHyMEjSnIvLnDeBd4XVW35W5IvmN3GIVD6fGlSUq+PFNtYEAgZ0dlsg64myG420LrMOqm4awPa0ib2hUZ26sJjaoG+zpkYwo0X9xhdAQ2iMjOnAyrNcZbMSNiuLvp3QQFOIbSBgUE0T+iP3se2ZNFwcuH4Dbu1p7PX+rHS7dHsOuvU9wycSljv9rGyaQLuXwVxhR83gyr7Z5rUZgCLS4O+vaF2bOhquddCx4LKxdGcKlgklKSCAwIJCklieBSwVn3Y2QwBLdYWBh9w+DGJlV5ZcFO3vttD19uPMBTNzbizlbhFLOFmoxxK8frYeRKVKZAGTcOli6Fsd51LXgl/kw8w1oNY8V9KxjWahiHTh/KuhBkOgS3QpmSvNA7gq+Gd6RWpTI89fkmek5exuq9tiysMe5404cxBcfoqOtU9SoRqQAsUNXWuRlgTlgfRu5y00UAOL7Mn82kayE/UlW+WLKNl77czKFS5bilWRgjuzeiRkVLMWKKHl/0YbRV1YeBJABVPQbY2MQizE0XAf37w54suhbyIxGh16eT+WnKvTxy7nd+3B7P9a/9zIvfbbf+DWOcvOnDyNF6GKbwyaCLIFf6MXKVy61SEPDYfx/nb+UqM6HLYKamdGHO6j955Pr63N22FiUDvPmOZUzhkmfrYZjCKZMugoLDza1Stdtu5LWpT/DV8I5cFRbMmK+20fW1n/ly40FLoW6KrCzvMEQkQFWTVXWmiKzFsZaF4FgPY3uuR2jytbkuiWEmF9QV3jO5VWoKzLy/LUt+T+A/3+1gxKz1TP3lD/5xYyM61a+MiI2oMkWHJ01Sq4CWAKq6A9iRqxEZ4w+pt0rR0TB1qmOssJOIcG3DKnSuH8oXGw/w6oLfGTRjFW3qVOSpGxtafipTZGQ5SkpE1qtqizyKx6dslJTJDeeSU5i9aCsTF2zncMmydG4QymM3NCCyhtvl6I0pcDIaJeXJHUaoiDye0Zuq+lqOIjOmgCkVUJx7vnybO2e8z4fDxvJObAl6TV5G16uq8Mj1DYgIL+/vEI3JFZ5UGMWBsjj6LYwp2tKNqBr2xj8YULI0H7TtzdTSg+gxaSnXNarCI9fXp7ndcZhCxpMmqXXOpVELHGuSMj4XFwdPPgnz50NiomNkVe/eMGECp0Iq8cFve5m2dA/HEy/QuUEow6+rR2vr4zAFTE4m7tmdhTGpMhlRVS6wBMOvq8/Sf17HP25syOYDJ7jr7eX87Z3l/PJ7gqVTNwWeJxXG9b44kYg0FJENLo+TIvKoiFQUkYUissv5s0IG5bs7M+XuFpGRvojJmGzJYvJJ2VIBPHxtPZb+81qe7hTO3p37uGfGKm59cylfbzpIis3jMAWUx7mkfHpSx4zxA0Bb4GHgqKq+5KwIKqjqP93s/ztwAxALrAb6ZbU2hzVJGb976CHOTZvO/AfH8E6NdsQcPkOtSkHc37EOd7aqQemSxf0doTGXyfESrT4OphvwrKp2EJGdQBdVjRORMGCJqjZMt387YIyq3ujc/heAqr6Y2XmswjB+4yYzY4oUY2Hja3jnofGs33+cCkElGHh1LQa0q0WVcoF+CtSYy/ki+aAv9QVmOZ9foapxAM6fVdzsXx3402U71vnaZUQkWkTWiMiahIQEH4ZsjGNd8c7vd846vbqbdCPF7+5H90WfMPfB9nw6rB2talXkzcW76fjSYp78dCPb407m/gUYkwN5XmGISEngNuBTb4q5ec3trZGqTnWuPR4VGhqanRCNyZDH64pn0jkuIrSuXZFpg6L46Yku9G1Tg282xXHTG7/Sd+pyvt9yyPo5TL7kTbZaX7kJWKeq8c7teBEJc2mS+stNmVighst2OHAwl+M0Jk36dcWnrJnClDVTMl9XPJN0I6nqVC7D2J5NeeKGhnzy01Y+/HEbw2KOUj2kNP2vrsnfompQqWyp3LosY7zijyapfvx/cxTAl8Ag5/NBwBduyqwG6otIHecdSl9nOWPyRLbWFZ8715GRsXlzx0/XTI3plA8qwQPfTuXnNwfxduJaalYM4uXvd9LupZ94fPYG1u47ZsNyjd/l6R2GiAThGOn0gMvLLwFzROQ+HMvA3uXctxowTVVvVtVkERkO/IBj5vkMVd2al7Gboi3b64p7wqWDPADo/uazdOdZfq9Wj/9Nnse89QeYu/4AV4UF079tTXpGVqNcYImcn9cYL/lllFResVFSxpdun307YWXDiG4VzdS1U4k7Hcfcv2V81+CxTGaPU7Uqp88l8+WGg3y0Yh/b4k5SukRxejQPo2+bmrSoEWIp1o3P5ST5oDEGLqkcJt/iw8U/sli6sGypAO5uW5N+bWqwMfYEnyzZwZcrY5izJpYGV5SlT1QNereobn0dJtfZepPG5AceLF0oIkTWCOGl5R+yctJAXjy7iaCSATz/zXaufvFHHvjfGhZui+dCiq2cbHKHNUkZU1C4mQwI8Hu1esx55wvmbzjA4dPnqVSmJLdFVuOOluE0qRZsTVbGa/lqpndesQrDFCpZ9HVcSLnIL78n8NnaWH7c/hfnUy5Sv0pZerWozm3Nq1GjYpC/r8AUENaHYUxBl0VfR4nixbj+qiu4/qorOJF4gW82xzF3RQyv/LCTV37YSVStCvSMrMZNEWFUtv4Okw3Wh2FMQeJBXwc45nXc3bYmn+2cw6/v3M8/knZwMukCT3+xlTbjFzFw+kpmr97P8cTzeXwBpiCzJiljCqMM+jt2Vq/Pl9O/5KuNcew/mkhAMaFDvcrcEhHGDY2voEKZkn4I1uQ31odhTFGSRX+HqrLlwEm+3nyQbzbFEXvsLMWLCe2vrMSNTarSrfEVVAm2DLpFlfVhGFOUZNHfISJEhJcnIrw8I7s3YsvGGL6b8B7fBXdh9PzDPP3FFlrUCOHGJlW5ofEV1A0t6+cLMvmB9WEYUwB5lGbdw/4OESFi6qs8NetFfto/lwWPdeKxrg04l3yRF7/bwXWv/sz1ry7hxe+2s3rvUcukW4RZk5QxBdBD3zzEO2vf4YFWD/DWLW9l/0AZ9HUQGAhnzxJ7LJFF2+JZuD2elTFHSb6oVAgqQZeGVbiuURU6NQilfGnLa1XYWB+GMYVA+jTrqTJNs56ZLPo6XJ1MusAvvyfw47p9LNkcy7ESQRQvJrSqWYHODUPp3CDUJgoWEvltxT1jTDZkK816ZrLo63AVHFiCW5tV4/VV/2PNG3fz+ZllPNj5Ss6cT+aVH3Zy65tLafPCjzw+ZwNfbDjAkdPncnKpJh+yTm9jCpBcSbPuwUJPwCXNV8WBVpNepNWkF3kyMJC/4o/x8+8J/LLrMIt3/MXcdQcAaFItmI71KtOxfmWialWkdMni2Y/T+J01SRlTwORamvWseNh8lXJR2XLgBL+ui+HX71eyrkItLlxUShYvRstaIbS/sjLtr6xEs/AQSgZYI0d+ZH0Yxpice/BBx11IyZJw/jw88AC8lUGn+0MPwTvvkPjAQ6x85GmW/3GEZbsPsy3uJKpQukRxompX4Oq6lbi6bkUiqlsFkl/kiwpDREKAaUBTQIF7gUeBhs5dQoDjqhrppuxe4BSQAiS7u5j0rMIwxsduv93R7+HafJV+6dksRl4dO3OelXuOsPyPIyyPOcLv8acdxVLO07J2JVo3rEqbOhVpUaOCNWH5SX6pMD4AflXVac61uYNUBB86VQAADl1JREFU9bjL+68CJ1R1rJuye4EoVT3s6fmswjDGD7wYeQVw5PQ5Vv37P6zcfpBVkZ3YHlAeVQgoJjStXp6oWhWIql2BlrUqUKWczT7PC36f6S0iwUAnYDCAqp4Hzru8L0Af4Lq8iskYkzNxp+Lo+3lfZt85+/873r0YeUXp0lRKSuIm4CaARe9wolQZ1tZpzprXp7Nm7zE+XLGPaUsdo8BqVgyiZc0QWtaqQIsaFWgUVo4Sxa0ZK6/k5SipukAC8J6INAfWAo+o6hnn+9cA8aq6K4PyCiwQEQXeUdWp7nYSkWggGqBmzZq+jN8Yk864X8axdP9Sxv489tIJhJ6OvIqJuexupHzvXlw3YQLXOSuYc8kpbDlwknX7jrFm50GWrfqd+RscqUoCSxSjabXyRNYIIbJmCM3DQwivUNrmguSSPGuSEpEoYAXQQVVXisgbwElVfdr5/hRgt6q+mkH5aqp6UESqAAuBv6vqL5md05qkjMkdPp1A6GVHur7zDgeGPcq66CfYsP84G/48xpaDJzmf7FiatlKZkjQLL09EeAjNqpenWXh5S6ToJb/3YYhIVWCFqtZ2bl8DjFTVW0QkADgAtFLVWA+ONQY4raoTMtvPKgyTH8XFQd++MHu2+1aagiDuVBxPLniS+Tvmk5icSFBAEL2v6s2EbhO8nxPig47088kX2XnoFBtij7Pxz+Nsjj3Brr9OkZr2KrRcKSKql6dptWCaVC9Pk2rBVA+xO5GM+L0PQ1UPicifItJQVXcC1wPbnG93BXZkVFmISBmgmKqecj7vBlzWMW5MQTBuHCxdCmPHZvxFOr/z6QRC18ph8mT3+7hpukrrSAdK/l979x6kVX3fcfz9weXyLCzXZbPACgsBxUsqAYYxsRcbnMaYjNa2jgZic21KVDSdZJo0nUnTcZKxM04NMdbGkARtiPESmhrHMWZS1DQJVJM4gCAXlzu7uGjluoTbt3+cs7AsD8tZ9rns7vN5zTxzLjznPN+zwH6f8/v9zu9bNeDk7Lu3XjkJgEO3LWTt08tZ/Rd/zeqZf8aanXt5fv0bJ5PIiNxALh03nMvGD+eScclrat0wD+3tQqmf9F4ILE1HSDUBH0/33wI82vGNksYDiyPiOuAdwH+m3waqgB9ExLMli9qsADp/SX7wweSVfknuc3Yf3M2CWQtOe4CwaLrZkc7hw1QDs4HZi74Ai74AQ4bQtvcA61r28equfazdtZe1zfv5j19v4ffHkyxSNUBMrRvGxfU1TK8fzvT6Gi6ur2HciCG+G8EP7pmVTDdHm1pnWZquoNs/6GO33c6WJ55m7Yc/xbrrbuK15n2sb9nPrr2nsnvN4Couqq/hopoBTPvJ41z0uQVMmz6RuprB/TKRlL1JyqzSdedLsuXRRdPVacN7s/6g0zuRKmAqMPX+L3P9/V8+ecu399BR1u/ez/rd+9nQkiyfXdXMo5PnwrKNwEZqhlQxtW4YU8cOS5Z1w3jn2GE0jMpR1Q+H+zphmJVQ1tGm1j1nDO/N8oM+R7/IiOqBzJk8mjmTR59MLgHsqR7JxtqJbKydyKa6RjbeOI/l61t54jenumAHDhCTDrQy5V1TmdIwhim1Q5k8diiTa4cyZuigPntX4iYpM+uzejy8N+uQ3gzNXHsPHWVT6wFebz1A09JlvL5tD5unXMbWQcM5evzU79mawVU01g5l0phqGscMZVLVURrv+xqTFv0LY6c09Ipk4iYpM+t3mu5sOuvw3kyy3vJlaOYaUT2QWdPHM6vT8N9jGsCOuolsfmElW/YcZPOeg2x58xCrduzlmdXNyaity+fDt1eRG/gqE0dXc+Ho6nSZO7ndMCpH9aDy/sp2wjCzPqvHw3uzDOltd57NXFU33kjjvffSWF93appVgFyOI0eOsXNEHVtH1rN11Hi2jqxn++jxbHvvXH65aQ9tR4+fdvoxRw/S0DCWhroRNIzKMWFUjgkjTy1rhhS3XK4Thpn1aSUb3pshuTQPg5XNy7nhcBs618iGpiYGff7zTP7xj5m8eRfsfi1t5voM1NcTEbx58Ajb3jrE9rcOsWPJD9mxfivbNZO1xxv52drdHDl+4rRT1gypYsLIHOOHiNt/8m/M+s59BR1V4YRhZn1ax+JRD3zwHHcJRXb3i3dzTUszL157GX/ytaU9auaSRO2wwdSOHcnMjs1c6RNoJ4bkaN39FjvfbmPX223s/L90+fZhdq3ZwPHVawr+dKg7vc3Meui8O9+zPFvSnedKzjGFSlZn6/TufwOFzcxKrOnOJuZdPo/qqmoAqquqmf+u+Wy+a3PXBy5bljRvXXFFssz3IGJ3HuBpaoJ58ziRywEQuRzMnw+bzxFHRk4YZmY9VNC5tfJp73BfsSJZtrScJZD25NJGWxXE4baCPh3qPgwz69X6yuy+Re18zziaK/fVHN9//jAts+Ch2fDpl6F++YN85Kvf6/6083m4D8PMerXbboNvfavrMhmWKNS08+7DMLM+JZcDKZnR98SJZCkl+y2/YjeNOWGYWa+U9t9SnfQjU11d0P7bfqu9aWzFJ1ewYNYCWg6cpb/jPLgPw8x6Jc/ue36W3bws6fe5GR577IGC/rx8h2FmvVbWwUF2uo5VHQuppJ3ekkYCi4HLgQA+Abwf+BugNX3blyLimTzHXgssAi4gqcR3z7k+z53eZlZJCvTcXq/p9F4EPBsR04ErgHXp/vsiYkb6ypcsLgAeAD4AXAp8WNKlpQrazKwvKHa/T8kShqThwB8D3wGIiCMR8XbGw+cAmyKiKSKOAD8EbihOpGZmfVOx+31KeYcxhaTZ6XuSfidpsaSh6Z/dIWmVpO9KGpXn2AnA9g7bO9J9Z5D0aUkvS3q5tbU131vMzPqtYvb7lKwPQ9JsYAVwVUSslLQI2Ad8E9hD0qdxNzAuIj7R6dibgPdHxKfS7VuBORGxsKvPdB+GmVn39YY+jB3AjohYmW4/CcyMiN0RcTwiTgDfJml+ynfshR22G4BdRY3WzMxOU7KEEREtwHZJ7TWn5gJrJY3r8LYbgTV5Dn8JmCZpsqRBwC3AU0UN2MzMTlPqB/cWAkvTX/pNwMeBb0iaQdIktQX4WwBJ40mGz14XEcck3QH8lGRY7Xcj4tUSx25mVtE8+aCZmZ2mN/RhmJlZH+aEYWZmmfTrJilJrcDW8zy8lmS4byXxNVcGX3Nl6Mk1T4qIsZ139uuE0ROSXs7Xhtef+Zorg6+5MhTjmt0kZWZmmThhmJlZJk4YZ/dQuQMoA19zZfA1V4aCX7P7MMzMLBPfYZiZWSZOGGZmlokTRh6SrpW0XtImSV8sdzzFJulCScslrZP0qqS7yh1TKUi6IK3N8nS5YykFSSMlPSnptfTv+j3ljqnYJP1d+m96jaRHJQ0pd0yFltYRekPSmg77Rkv6maSN6TJfnaFuc8LopELLwR4DPhcRlwBXArdXwDUD3MWpMsGV4GwlkvslSROAO4HZEXE5ycSlt5Q3qqJYAlzbad8XgZ9HxDTg5+l2jzlhnKniysFGRHNE/DZd30/yiyRvRcP+QlID8EFgcbljKYUelkjuy6qAnKQqoJp+WEcnIl4E3uq0+wbg4XT9YeDPC/FZThhnylwOtj+S1Ai8G1jZ9Tv7vK8Dfw+cKHcgJdJVieR+KSJ2AvcC24BmYG9EPFfeqErmHRHRDMkXQqCuECd1wjiT8uyriLHHkoYBPwI+GxH7yh1PsUj6EPBGRPym3LGUUBUwE3gwIt4NHKRAzRS9VdpufwMwGRgPDJX0kfJG1bc5YZypIsvBShpIkiyWRsSycsdTZFcB10vaQtLk+D5J3y9vSEWXt0RyGeMphWuAzRHRGhFHgWXAe8scU6nsbq9mmi7fKMRJnTDOVHHlYCWJpG17XUT8a7njKbaI+IeIaIiIRpK/3/+OiH79zfNsJZLLGFIpbAOulFSd/hufSz/v6O/gKeCj6fpHgf8qxElLXaK116vQcrBXAbcCqyW9ku77UkQ8U8aYrPDylUjutyJipaQngd+SjAT8Hf1wihBJjwJXA7WSdgD/BNwDPC7pkySJ86aCfJanBjEzsyzcJGVmZpk4YZiZWSZOGGZmlokThpmZZeKEYWZmmThhmJlZJk4YVpEkjZH0SvpqkbSzw/avivB5H5PUKmlxun1152nVJS2R9FeF/uwO58+l13dEUm2xPsf6Lz+4ZxUpIt4EZgBI+gpwICLuLfLHPhYRdxT5M5BUFRHHOu+PiDZgRjolilm3+Q7DrBNJB9Ll1ZJekPS4pA2S7pE0X9L/Slot6Z3p+8ZK+pGkl9LXVQWIYW46q+zqtEDO4HT/lva7A0mzJT2frn9F0kOSngMekXRZGucrklZJmtbTmMx8h2HWtSuAS0jqDTQBiyNiTlqVcCHwWZLCRPdFxP9ImkgyrcwlGc79Rx2mYgGYCDydVoVbAsyNiA2SHgE+QzIle1dmAX8YEW2S7gcWRUT7VCAXZL1gs7NxwjDr2kvtdQUkvQ6011NYDfxpun4NcGkyvx0AwyXVpMWouvKLiPhQ+4akJenqxSSzrG5Itx8GbufcCeOptNkJ4NfAP6aFopZFxMZzHGt2Tm6SMuva7zusn+iwfYJTX7gGAO+JiBnpa0KGZNGVfDVZ2h3j1P/bzvWpD7avRMQPgOuBNuCnkt7Xg3jMACcMs0J4DjjZmS1pRg/P9xrQKGlqun0r8EK6voWk6QngL892AklTgKaI+AbJVNd/0MOYzJwwzArgTmB22rm8FljQk5NFxGGSqcefkLSa5G7m39M//mdgkaRfAMe7OM3NwJq0j2Q68EhPYjIDT29uVhKSPgbMLsWw2gyxbElj2VPuWKxv8R2GWWm0AR9of3CvHNof3AMGkty1mHWL7zDMzCwT32GYmVkmThhmZpaJE4aZmWXihGFmZpn8P/vtejc2d6ZzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Temp_initial = 85\n",
    "Temp_dead = 98.6\n",
    "Temp_ambient = 65\n",
    "t_ana = np.linspace(0,10,100) #create array of t values from 11 am to 9 pm with 100 steps\n",
    "Temp_ana = Temp_ambient + (Temp_initial - Temp_ambient)*np.exp(K*t_ana) #Analytical via integration of Newton's cooling\n",
    "N = 0\n",
    "num = []\n",
    "for N in [5,10,30]:\n",
    "    t = np.linspace(0,10,N) #create array of t values from 11 am to 9 pm with N steps\n",
    "    Temp_num = np.empty(len(t)) #create empty array with length of t\n",
    "    delta_t = np.diff(t) #establish size of time steps from t\n",
    "    for i in range(0,len(t)-1):\n",
    "        Temp_num[0] = Temp_initial #Establish first component in array as initial condition\n",
    "        Temp_num[i+1] = Temp_num[i] + (K*(Temp_num[i] - Temp_ambient)*delta_t[i]); #Numerical solution approach\n",
    "        num.append(Temp_num)\n",
    "\n",
    "    if N == 5:\n",
    "        plt.plot(t,Temp_num,'b*',label='Numerical Temperature, Step = {:.0f}'.format(N))\n",
    "    if N == 10:\n",
    "        plt.plot(t,Temp_num,'g*',label='Numerical Temperature, Step = {:.0f}'.format(N))\n",
    "    if N == 30:\n",
    "        plt.plot(t,Temp_num,'r*',label='Numerical Temperature, Step = {:.0f}'.format(N))\n",
    "plt.plot(t_ana,Temp_ana,label='Analytical Temperature')\n",
    "plt.title('Temperature of Cooling Body')\n",
    "plt.xlabel('Time [Hours]')\n",
    "plt.ylabel('Temeperature [$^oF$]')\n",
    "plt.legend(loc='best');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd5ycZbn/8c93d9M22eymb8gm2QRSCElACIFQBEEUjjTroSMiHlT0iAcVlJ/9WM6xHY8FEanSVRA5dKQoBkiCpBFKCAnZ9BDSe/b6/fE8m0yWLbPZnZ0t3/frNa+ZuZ8y1z0b5uIuz/0oIjAzM9tXBfkOwMzM2jcnEjMzaxYnEjMzaxYnEjMzaxYnEjMzaxYnEjMzaxYnEmtVko6XVNWM46+R9P9aMiZrOZJC0gH7eOy5kh5p6Zgs95xIrMkkLZS0RdJGScsl3SipVw4+5+OS/p5ZFhGXRsR3cvBZmXWqefyipT8ny1hGSKqW9Kt8fH5rkFSZJp2imrKIuDUi3pfPuGzfOJHYvjotInoBhwDvAq7Kczwt4bSI6JXxuKyunTJ//Boqa0gj+18AvA2cJalbU85rlg9OJNYsEbEceJgkoQAgqZukH0l6U9KKtDuqR13HS7pS0uuSNkh6SdIH0/IDgWuAKWnrYG1afqOk76av50k6NeNcRZJWSzo0fX+kpH9IWitppqTj96WOacvoGUk/lbQG+GY9ZQWSrpa0SNJKSTdLKk3PUfN/4BdLehP4awMfeQFwNbADOK1WLAdJelTSmvS7/WpaXijpqxnf5QxJQ+v6P39JT0r6ZB11WytpgaSj0vLFaT0urOvYjOP3ajVmbPuApH9KWp+e65sZm59On9emf98ptc+VxjFN0rr0+ahacXwnjX2DpEck9W/gO7UcciKxZpFUAZwCzM8o/iEwmiS5HAAMAb5ezyleB44FSoFvAb+XNDgi5gGXAlPT1kFZHcfeDpyd8f79wOqIeEHSEOD/gO8CfYErgD9KGrBvNeUIYAEwEPjPeso+nj7eA4wEegG1u8eOAw5MY30HSccCFcAdwF0kSaVmWwnwGPAQsB/Jd/t4uvmLJN/FvwC9gU8Am5tQt1lAP+C29LMPT89/HvCLfey63JTGXwZ8APi0pDPTbe9On8vSv+/UzAMl9SX5+/08jesnwP9J6pex2znARSTff1eSv7HlgROJ7at7JW0AFgMrgW8ASBJwCXB5RKyJiA3A94Cz6jpJRNwdEUsjojoi7gReAyZnGcNtwOmSitP356RlkPwAPhARD6TnfhSYTvJD21Cd1mY8LsnYtjQi/jcidkbElnrKzgV+EhELImIjSXffWbW6sb4ZEZsyzlHbhcCDEfF2WpdTJA1Mt50KLI+IH0fE1ojYEBHPpds+CVwdEa9EYmZEvNVAXTO9ERE3RMQu4E5gKPDtiNgWEY8A20mSSpNExJMRMTv9/meRJP7jsjz8A8BrEXFL+v3eDrzM3i20GyLi1fS7vIuMVrG1LicS21dnRkQJcDwwFqjpVhgAFAMzan6QSf4Pus6WgKQLJL2Yse/4jHM1KCLmA/OA09Jkcjp7Eslw4KOZiQE4BhjcSJ3KMh6/zdi2uI79a5ftByzKeL8IKAIGNXIeANLuv48Ct6b1mwq8SZIgIfmBf72ewxva1pgVGa+3pJ9du6zJLRJJR0h6QtIqSetIWpjZdj/V/i5J3w/JeL884/XmfYnRWoYTiTVLRDwF3Aj8KC1aTfLDc1DGD3JpOjC/F0nDgd8ClwH90u6rOYBqTp9FCDXdW2cAL6XJBZIf7FtqJYaeEfGDfatpnbHULltKksBqDAN2svcPdUN1+iBJt9SvlMyGW07yw1nTvbUY2L+eY+vbtil9Ls4oK28ghsZsasK5bgPuA4ZGRCnJmFe2f9va3yUk3+eS7EO11uJEYi3hZ8BJkg6JiGqS5PDTmi4ZSUMk1TUm0JPkB2VVut9FJC2SGiuACkldG/jsO4D3AZ9mT2sE4PckLZX3pwPR3ZVcw1Kxj3XMxu3A5Uqm7/Yi6dK7MyJ2Znn8hcD1wASSbppDgKOBQyRNAO4HyiV9QcmEhhJJR6THXgd8R9IoJSZK6hcRq0h+fM9Lv4dPUH8yysaLwIckFSu5XuTiBvYtAdZExFZJk9nTsoLkb15NMpZUlweA0ZLOUTKJ4l+BcSTfgbUxTiTWbOmP1c1AzYWCXyEZfH9W0nqSAeIxdRz3EvBjYCpJ0pgAPJOxy1+BucBySavr+exl6fFHkfTv15QvJmmlfJXkR2sx8CUa/jf/F+19Hck9jVS9tuuBW0hmJL0BbAU+l82B6eSAE4GfRcTyjMcMkq7BC9PxppNIxgmWk4wnvSc9xU9IxgkeAdYDvwNqZspdQlL3t4CDgH80sV6ZfkoyZrICuIm0G64enwG+nY6lfT2ND4CI2EwyQeGZtOvxyMwD0/GdU4H/SOP+MnBqRNT578DyS76xlZmZNYdbJGZm1ixOJGZm1ixOJGZm1ixOJGZm1ixNWmiuverfv39UVlbmOwwzs3ZlxowZqyOi0WWFOkUiqaysZPr06fkOw8ysXZFUe3WBOrlry8zMmsWJxMzMmsWJxMzMmsWJxMzMmsWJxMzMmsWJxMzMmsWJxMzMmsWJpAGPz1vBr56c3/iOZmadmBNJA55+dRW/fnJf715qZtY5OJE0oKR7FzZu20l1te/ZYmZWHyeSBpR0LyICNm3P9k6pZmadjxNJA0q6dwFgw1YnEjOz+jiRNKCke7KmpROJmVn9nEgaUJNINm7bkedIzMzarpwmEkknS3pF0nxJV9axfaykqZK2Sbqi1rbrJa2UNKeec18hKST1z1X8vXskXVvrt7hFYmZWn5wlEkmFwC+BU4BxwNmSxtXabQ3weeBHdZziRuDkes49FDgJeLOl4q1LaZpI1m1xi8TMrD65bJFMBuZHxIKI2A7cAZyRuUNErIyIacA7fqkj4mmSRFOXnwJfBnI6L7csTSRrN2/P5ceYmbVruUwkQ4DFGe+r0rJmkXQ6sCQiZjay36ckTZc0fdWqVfv0Wb13t0jctWVmVp9cJhLVUdasFoSkYuBrwNcb2zciro2ISRExacCARm85XKcuhQX07Frori0zswbkMpFUAUMz3lcAS5t5zv2BEcBMSQvTc74gqbyZ561XWXFXJxIzswYU5fDc04BRkkYAS4CzgHOac8KImA0MrHmfJpNJEbG6OedtSO8eXVi3xWMkZmb1yVmLJCJ2ApcBDwPzgLsiYq6kSyVdCiCpXFIV8EXgaklVknqn224HpgJj0vKLcxVrQ0p7FLlFYmbWgFy2SIiIB4AHapVdk/F6OUn3VF3Hnp3F+SubGWKjynp0ZcHqjbn+GDOzdstXtjeitEcXt0jMzBrgRNKIsuIurN3sRGJmVh8nkkaUFndh285qtmzfle9QzMzaJCeSRgwq6Q7Ayg1b8xyJmVnb5ETSiEG9k0SyfJ0TiZlZXZxIGjGodzcAVmzYludIzMzaJieSRgwqTVokK9wiMTOrkxNJI0q6FdGjSyEr1juRmJnVxYmkEZIoL+3OcicSM7M6OZFkYWBJN1au9xiJmVldnEiy4BaJmVn9nEiyMKh3d1as30pETm/IaGbWLjmRZGFQ7+5s21ntNbfMzOrgRJKF3deSeJzEzOwdnEiyUF5zdbvHSczM3sGJJAs1y6T4WhIzs3dyIsnCwJquLV/dbmb2Dk4kWehWVEif4i6s8ArAZmbv4ESSpSF9evDmmi35DsPMrM1xIsnS6EElvLJ8fb7DMDNrc5xIsjS2vIQV67fx9qbt+Q7FzKxNcSLJ0pjy3gC8vHxDniMxM2tbnEiydGB5CYC7t8zManEiydKAkm70Ke7CKyvcIjEzy+REkiVJjCkvYd4yJxIzs0xOJE0wtrw3r67YQHW1VwE2M6vhRNIEY8tL2Lx9F1Vv+3oSM7MaTiRNMCYdcJ/nAXczs92cSJpg9KCamVseJzEzq+FE0gQ9uxUxvF+xE4mZWQYnkiYaM6jEXVtmZhmcSJpobHkJC1dvYuuOXfkOxcysTXAiaaKxg3tTHTB/5cZ8h2Jm1iY4kTTR7plby9y9ZWYGTiRNVtmvJ92KCjzgbmaWciJposICJfcm8ZpbZmaAE8k+8ZpbZmZ75DSRSDpZ0iuS5ku6so7tYyVNlbRN0hW1tl0vaaWkObXK/1vSy5JmSbpHUlku61CXseUlrN64jdUbt7X2R5uZtTk5SySSCoFfAqcA44CzJY2rtdsa4PPAj+o4xY3AyXWUPwqMj4iJwKvAVS0Vc7beNSzJXf94/a3W/mgzszYnly2SycD8iFgQEduBO4AzMneIiJURMQ3YUfvgiHiaJNHULn8kInamb58FKlo88ka8a2gfBpR046E5y1r7o83M2pxcJpIhwOKM91VpWUv6BPBgXRskfUrSdEnTV61a1aIfWlAg3n/QIJ54eRVbtvvCRDPr3HKZSFRHWYvdyEPS14CdwK11bY+IayNiUkRMGjBgQEt97G6njB/Mlh27eOrVlk1SZmbtTS4TSRUwNON9BbC0JU4s6ULgVODciMjLXaYmj+hLWXEXHp67PB8fb2bWZuQykUwDRkkaIakrcBZwX3NPKulk4CvA6RGxubnn21ddCgs46cBBPDZvBdt3VucrDDOzvMtZIkkHxC8DHgbmAXdFxFxJl0q6FEBSuaQq4IvA1ZKqJPVOt90OTAXGpOUXp6f+BVACPCrpRUnX5KoOjTllQjkbtu7kmddX5ysEM7O8K8rlySPiAeCBWmXXZLxeTj2zriLi7HrKD2jJGJvj6AP606tbEQ/NXs57xgzMdzhmZnnhK9uboVtRISeMHcij81awc5e7t8ysc3IiaaZTxpezZtN2nl/4jktezMw6BSeSZjpuzAC6dyngoTmevWVmnZMTSTMVdy3i+NEDeWjOcqqr8zIT2cwsr5xIWsDJ48tZuWEb/1y8Nt+hmJm1ugYTiaRCSZe3VjDt1QkHDqRLobz2lpl1Sg0mkojYRa2FFu2denfvwjEH9OfBOcvJ04X2ZmZ5k03X1jOSfiHpWEmH1jxyHlk7c8r4wVS9vYW5S30vdzPrXLK5IPGo9PnbGWUBnNDy4bRf7x03iMJ7xINzljF+SGm+wzEzazWNJpKIeE9rBNLe9e3ZlSNG9OWhOcv50vvH5jscM7NW02jXlqRSST+pubeHpB9L8v9y1+GU8eW8vmoTr63w/dzNrPPIZozkemAD8LH0sR64IZdBtVfvP6gcCR70xYlm1olkk0j2j4hvpLfMXRAR3wJG5jqw9mhg7+4cNqyPE4mZdSrZJJItko6peSPpaGBL7kJq304eX868ZetZ9NamfIdiZtYqskkklwK/lLRQ0kKS+4H8W06jasfef1A5gNfeMrNOo7Er2wuAMRFxMDARmBgR74qIWa0SXTs0tG8xEytKuXPaYt850cw6hcaubK8mucshEbE+Iny1XRYuP2k0C1Zv4rd/W5DvUMzMci6brq1HJV0haaikvjWPnEfWjr1nzEBOGV/Ozx9/jcVr8nZbeTOzVpFNIvkE8FngaWBG+piey6A6gq+fNo6iAvGN++Z6/S0z69CyGSM5LyJG1Hp4+m8jBpf24PKTRvPXl1fy8NwV+Q7HzCxnshkj+VErxdLhfPyoSsaWl/Ctv8xl07ad+Q7HzCwnsunaekTShyUp59F0MEWFBfznByewbN1WfvbYq/kOx8wsJ7JJJF8E7ga2SVovaYMkz97K0mHD+3D25KFc/8xC5i3z12ZmHU+jiSQiSiKiICK6RkTv9H3v1giuo/jKyWMp7dGFq++d4/u6m1mHU28ikXRexuuja227LJdBdTRlxV256pSxzFj0NnfPWJzvcMzMWlRDLZIvZrz+31rbPpGDWDq0jxxWweTKvnz/wZdZs2l7vsMxM2sxDSUS1fO6rvfWCEl894Pj2bh1J99/YF6+wzEzazENJZKo53Vd7y0LoweV8MljR3L3jCqmLVyT73DMzFpEQ4lkrKRZkmZnvK55P6aV4utwPn/iAQwp68HX7pnNjl1e1NHM2r+G7tl+YKtF0YkUdy3im6cfxCU3T+d3f3+DS4/bP98hmZk1S72JJCIWtWYgnclJ4wbx3gMH8T+PvcapEwdT0ac43yGZme2zbC5ItBz45unjAPjWX17KcyRmZs3jRJInFX2K+ff3juLRl1bw2Ete1NHM2q+sEomkHpI8wN7CLj5mBKMH9eIb981l83Yv6mhm7VOjiUTSacCLwEPp+0Mk3ZfrwDqDLoUFfPfMCSxZu4WfPz4/3+GYme2TbFok3wQmA2sBIuJFoDJ3IXUuk0f05SOHVXDd3xbw6ooN+Q7HzKzJskkkOyNiXc4j6cSuOmUsvboXcfW9c3w3RTNrd7JJJHMknQMUShol6X+Bf2RzckknS3pF0nxJV9axfaykqZK2Sbqi1rbrJa2UNKdWeV9Jj0p6LX3uk00sbVm/Xt248uSxPP/GGv74wpJ8h2Nm1iTZJJLPAQcB24DbgHXAFxo7SFIh8EvgFGAccLakcbV2WwN8nrrvwngjcHId5VcCj0fEKODx9H2797FJQzl0WBnfe2Aeb3tRRzNrR7JJJGMi4msRcXj6uDoitmZx3GRgfkQsiIjtwB3AGZk7RMTKiJgG7Kh9cEQ8TZJoajsDuCl9fRNwZhaxtHkFBeI/PziBdVt28F8Pv5zvcMzMspZNIvmJpJclfUfSQU049xAg8+YbVWlZcw2KiGUA6fPAFjhnm3Dg4N5cdFQltz+/mBmL3s53OGZmWcnmDonvAY4HVgHXSpot6eoszl3XUvOtNpIs6VOSpkuavmrVqtb62Gb7wkmjKe/dna/dM5sNW9/RUDMza3OyuiAxIpZHxM+BS0muKfl6FodVAUMz3lcAS5sc4TutkDQYIH1eWddOEXFtREyKiEkDBgxogY9tHb26FfG9D43ntZUb+fCv/8Gbb23Od0hmZg3K5oLEAyV9M5099QuSGVsVWZx7GjBK0ghJXYGzgJa4kPE+4ML09YXAn1vgnG3KCWMHcfMnJrNi/TbO/NUzPP+G711iZm1XNi2SG4C3gfdFxHER8euIqLMVkCkidgKXAQ8D84C7ImKupEslXQogqVxSFcltfa+WVCWpd7rtdmAqMCYtvzg99Q+AkyS9BpyUvu9wjj6gP/d+9mjKenTh3Oue5a5pvte7mbVN6gwXwE2aNCmmT5+e7zD2ybrNO/jsbS/w9/mrueTYEVx5yoEUFvhOx2aWe5JmRMSkxvart0Ui6a70eXbG3RFn1bxvyWCtfqXFXbjxosO5cMpwfvu3N7jk5ukehDezNqWhOyT+e/p8amsEYvUrKizgW2eM54BBJXzzvrl8+Nf/4HcXHs7Qvr4hlpnlX70tkpprNYDPRMSizAfwmdYJzzKdf+RwbrpoMsvXbeWMX3oQ3szahmwG20+qo+yUlg7EsnPMqFqD8NM9CG9m+dXQGMmnJc0mmTWVOUbyBuAxkjwaOaAX93zmaI4Y0Y8v/2EW33tgHruqO/6kCTNrmxoaI7kNeBD4PnsvjLghItynkmelxV244aLD+c79L3Ht0wt4feVGfnbWIZR075Lv0Mysk2lojGRdRCyMiLPTcZEtJEuc9JI0rNUitHp1KSzg22eM5ztnHMSTr67iI7+eyuI1vhLezFpXVrfaTS/+ewN4ClhI0lKxNuL8KZXcdNFklq3bwhm/fIZpC91gNLPWk81g+3eBI4FXI2IEcCLwTE6jsiarGYQv7dGFc377LHd7EN7MWkk2iWRHRLwFFEgqiIgngENyHJftg5EDenFvOgj/pT/M4vsehDezVpBNIlkrqRfwNHCrpP8BduY2LNtXNYPwF0wZzm+eXsCnbp7Oxm3+c5lZ7mSTSM4gGWi/HHgIeB04LZdBWfPUHoT/8K/+4UF4M8uZbG5stSkidkXEzoi4KSJ+nnZ1WRuXOQh/pgfhzSxHGrogcYOk9RmPDZnPrRmk7buaQfjePbpw7m+f4w8zqvIdkpl1MA1dR1ISEb0zHiWZz60ZpDVPzSD84SP6cMXdMz0Ib2YtKqtb7Uo6RtJF6ev+kkbkNixracly9JM5/8hkEP7fbvEgvJm1jGwuSPwG8BXgqrSoK/D7XAZludGlsIDvnDmeb59xEE+8soqP/NqD8GbWfNm0SD4InA5sAoiIpUBJLoOy3LpgSiU3XnQ4S9cmg/DTPQhvZs2QTSLZHsn9eANAUs/chmSt4dhRA7gnHYQ/x4PwZtYM2SSSuyT9BiiTdAnwGHBdbsOy1rD/gF7c85mj9gzCP+hBeDNrumyuI/kR8Afgj8AY4OsR8fNcB2ato6y4KzdeNJnzjhzGb55awL/dMsOD8GbWJFnN2oqIRyPiSxFxBfBXSefmOC5rRV0KC/jumRPSQfiVfOTX/6DqbQ/Cm1l2GrogsbekqyT9QtL7lLgMWAB8rPVCtNZSMwi/JB2En7HIg/Bm1riGWiS3kHRlzQY+CTwCfBQ4IyLOaIXYLA+OHTWAez5zNL26FXH2tc/xRw/Cm1kjGrrV7siImAAg6TpgNTAsIja0SmSWNwcM7MW9nz2az9z6Av9x90xeW7mRL7x3FN27FOY7NDNrgxpqkeyoeRERu4A3nEQ6j7Lirtz0icmce8QwrnnqdY76wV/5wYMv+wJGM3sHJZeI1LFB2kV6ESIgoAewOX0d7Wm9rUmTJsX06dPzHUa7FBFMff0tbpq6kEdfWgHACWMHccGU4RxzQH8KCpTfAM0sZyTNiIhJje1Xb9dWRLgfw5DEUQf056gD+rNk7RZue24Rdzy/mMfmrWBE/56cd+RwPnJYBaU9uuQ7VDPLk3pbJB2JWyQta9vOXTw4ezk3T13IC2+upUeXQs58136cf2Ql4/ZrNw1VM2tEti0SJxJrljlL1nHL1EX8eeYStu6oZtLwPlxwVCUnH1RO16KsLlMyszbKiSSDE0nurdu8g7tnLOaWZxex6K3N9O/VjXMmD+WcI4ZTXto93+GZ2T5wIsngRNJ6qquDp19bxc1TF/HEKyspkHjfuEGcP2U4U0b2Q/LgvFl70ezBdrN9UVAgjh8zkOPHDOTNtzZz63OLuHP6Yh6cs5xRA3tx/pThfOjQCnp18z89s47CLRLLua07dvGXmUu55dlFzKpaR8+uhXz4sArOP3I4owb51jZmbZW7tjI4kbQdLy5ey81TF3L/zGVs31XNlJH9uGDKcE4aN4iiQg/Om7UlTiQZnEjanrc2buPO6Yu59dk3WbJ2C+W9u3POEcM4a/JQBpZ4cN6sLXAiyeBE0nbtqg7++vJKbp66kL+9tpouheLk8YO5YMpwJg3v48F5szzyYLu1C4UF4qRxgzhp3CAWrNrI7599k7tnLOYvM5dy4ODeXDBlOGccsh/FXf1P1aytymmntKSTJb0iab6kK+vYPlbSVEnbJF2RzbGSDpH0rKQXJU2XNDmXdbDWM3JAL75+2jie++qJfO+DE4gIrvrTbI743uN8+y8v8cbqTY2fxMxaXc66tiQVAq8CJwFVwDTg7Ih4KWOfgcBw4Ezg7fS2vg0eK+kR4KcR8aCkfwG+HBHHNxSLu7bap4hg+qK3uXnqIh6cvYyd1cGxo/pzwZRKThg7kEIvGGmWU22ha2syMD8iFqQB3QGcAexOJBGxElgp6QNNODaAmgWdSoGlOayD5ZEkDq/sy+GVfVl56oHc8fxibnvuTS65eTpDynpw3pHD+dfDh9K3Z9d8h2rWqeUykQwBFme8rwKOaIFjvwA8LOlHJF1zR9V1AkmfAj4FMGzYsOyjtjZpYEl3Pn/iKD59/P489tIKbp66iB8+9DI/fexVTp04mAumVHLI0LJ8h2nWKeUykdTV75BtP1pDx34auDwi/ijpY8DvgPe+Y+eIa4FrIenayvJzrY3rUljAKRMGc8qEwby6YgO3TF3En16o4k8vLGFiRSnnHzmc0w7ez3dzNGtFuRxsrwKGZryvIPtuqIaOvRD4U/r6bpJuMOuERg8q4TtnjufZr57It884iM3bd/GlP8ziyO8/zvcfnOe7OZq1klwmkmnAKEkjJHUFzgLua4FjlwLHpa9PAF5rwZitHSrp3oULplTy6OXv5rZLjmDKyH5c97c3ePd/P8HFN07jyVdWUl3tRqlZruSsaysidkq6DHgYKASuj4i5ki5Nt18jqRyYTjJ4Xi3pC8C4iFhf17HpqS8B/kdSEbCVdBzETBJH7d+fo/bvz7J1W7j9uTe57fnFPH7DNCr7FXPekcP56GFDKS323RzNWpKvbLcObfvOah6cs4xbpi5i+qK36d6lgDMOHsL5U4YzfkhpvsMza9O8REoGJxIDeGnpem55diH3/nMpW3bs4tBhZVx4VCUnjy+nW5EH581qcyLJ4ERimdZt2cEfZlTx+2cX8cbqTfTv1ZWzDh/GOUcMY7+yHvkOz6zNcCLJ4ERidamuDv42fzW3TF3I4y8nd3N874EDuWBKJUft77s5mrWFK9vN2rSCAnHc6AEcN3oAi9ds5tbn3uTOaW/y8NwV7D+gJxdMqeRDhw6hpLsH580a4haJWYatO3bxf7OWcfOzi5i5eC09uxbywUOHcMGUSkb7bo7WybhrK4MTie2LWVVruXnqIu6buZTtO6s5YkRfLphSyfsOGkQX383ROgEnkgxOJNYcazZt567pi/n9s4uoensLg3p345zJwzl78lAG9vbdHK3jciLJ4ERiLWFXdfDkKyu5eeoinnp1FUUF4uTx5VwwpZLDK303R+t4PNhu1sIKC8SJBw7ixAMHsXD1Jn7/7CLumr6Y+2ctY2x5CedPGc6ZhwyhZzf/Z2Wdi1skZs2wZfsu/vziEm6euoiXlq2npFsRpx+yH5Mq+zCxoowR/XpS4BtwWTvlrq0MTiSWaxHBC28md3N8ZO4KtuzYBUBJ9yImDCllYkUZB1eUMnFoGfuVdnc3mLUL7toya0WSOGx4Xw4b3pedu6qZv2ojsxavY2bVWmZVreN3f1/Ajl3J/7T179WViRVlTBhSysFDkyTTv1e3PNfAbN85kZi1sKLCAsaW92ZseW8+dnhyW51tO3cxb9kGZletZWbVOmZVreWJV1ZS0yEwpKwHEyv2tFzGV5TS2xdCWjvhRGLWCnQG5g4AAA+/SURBVLoVFXLI0DIOGVrG+WnZpm07mbNkHbOq9rRcHpyzfPcxIwf05OCKst0J5qD9evvOj9YmOZGY5UnPbkUcMbIfR4zst7vs7U3bmbVkHbMWJy2XZ+av5p5/LgGgqECMHlSyuztswpBSxpSX+OJIyzsPtpu1ccvXbU1bLEmrZVbVOtZt2QFAt6ICxu3Xe6+Wy8j+nilmLcOztjI4kVhHEhG8uWZzMtayOEkuc5auY/P2dKZYtyLGDyll4tDS3QlmSFkPzxSzJvOsLbMOShLD+/VkeL+enH7wfkBy1f38lRv3arlc//c3ds8U69ez657BfM8UsxbmRGLWARQWiDHlJYwpL+Fjk/bMFHt52QZmZcwUe+rVVVSnnRD7lXZnYkXZ7pbLBM8Us33kRGLWQXUrKuTgoWUcXGum2Nyl6/dKLg/NzZgp1r/nXi2Xg/Yr9Uwxa5QTiVkn0rNbEZNH9GXyiL67y9Zu3p4O4ifJZeqCt7j3xaVA0tIZPagkuSo/HW/xTDGrzYPtZvYOK9ZvZWY6kD+zai2zl6xj7ea9Z4pNHLKn5TKyfy/PFOuAPGsrgxOJWfPUzBTLbLnMWbJnplivbkWMH1IzDTlpuVT08Uyx9s6ztsysxWTOFDstY6bY66s27m65zKpayw3PLGT7rmoA+mbOFKsoZUJFKQNLfCOwjsiJxMz2Sc34yehBJXw0Y6bYK8s37HWNy9OvvrZ7ptjg0u4ZySWZKVbawzPF2jsnEjNrMd2KCtOurTI4cjgAm7fvZM6S9RlX5q/l4bkrdh8zInOmWEUyU6xHV88Ua0+cSMwsp4q7vnOm2LrNO5i1JB3MX7yW5xas4c8ZM8VGDeyVjLek17h4pljb5sF2M2sTVq7fuvvalprnmpliXYsKOHBw793TkA+uKGXkgF4UeqZYTnnWVgYnErP2JyJYvGbLXsu+zFmyjk3pTLGeXQsZP6SUg4cms8QOrijzTLEW5llbZtauSWJYv2KG9Svea6bYglUb92q53FhrptiEIaV7LqAc6plircEtEjNr17bvrE5niu1puby6YkOdM8UmVpQycUgZpcWeKZYNt0jMrFPoWlTAhPQ6FdgzU2zu0vV7XeOSOVOssl/x7sRy8NDk7pPFXf1zuK/8zZlZh1PctYjDK/tyeOXeM8VmL1m3u+UybeEa7puZzBQrEIweVLLXNS5jykvoWuSZYtlw15aZdVorN2xl1uK9Z4q9XTNTrLCAA/frvdeClft3splinrWVwYnEzLIREVS9XTNTLLnGpfZMsYMyBvMPrihjaN+OO1PMYyRmZk0kiaF9ixnat5hTJ+49Uyxzwcqbpi5i+843AOhT3IUJ6bUtNde4DOzduWaKuUViZtZE23dW8+qKdKbY4mTc5bWVG9mVThUr791990B+e54p1iZaJJJOBv4HKASui4gf1No+FrgBOBT4WkT8KJtjJX0OuAzYCfxfRHw5l/UwM8vUtaiA8UNKGT+klHOPSMq2bN/F3KXrdo+1zKpaxyMv7ZkpNjydKVbTchk/pOPMFMtZLSQVAr8ETgKqgGmS7ouIlzJ2WwN8Hjgz22MlvQc4A5gYEdskDcxVHczMstWjayGTKvsyKXOm2JYdzK65OVjVOmYsXMNfMmaKjRqYzhQbmiSYseW92+VMsVymw8nA/IhYACDpDpIEsDuRRMRKYKWkDzTh2E8DP4iIbRnnMDNrc0p7dOGYUf05ZlT/3WUrN2xNk0vScnn85ZXcPaMKSGeKDS7Z6xqX9jBTLJeJZAiwOON9FXBECxw7GjhW0n8CW4ErImJa7RNI+hTwKYBhw4Y1LXIzsxwZWNKdEw/szokHDgL2zBTbM5i/lnv+uYRbnl0EQHG6ptjEIXtaLsP6FrepmWK5TCR11TLbkf2Gji0C+gBHAocDd0kaGbVmDUTEtcC1kAy2Z/m5ZmatKnOm2AcmDgagujpYsHojMzOucbn52UVs/3syU6ysuEu6ptielsugPM4Uy2UiqQKGZryvAJa2wLFVwJ/SxPG8pGqgP7CqeeGambUNBQXigIElHDCwhA8fVgHAjl3JmmKZ05B//dTru2eKDerdba/B/IkVpZQVd22VeHOZSKYBoySNAJYAZwHntMCx9wInAE9KGg10BVa3ZOBmZm1Nl8I9M8XOOSLprt+yfRcvLVu3u+Uyq2odj9aaKfb9D03gqP3713faFpGzRBIROyVdBjxMMoX3+oiYK+nSdPs1ksqB6UBvoFrSF4BxEbG+rmPTU18PXC9pDrAduLB2t5aZWWfQo2shhw3vy2HD954pNqdmTbHF61plGX1fkGhmZnXK9oLE9jdh2czM2hQnEjMzaxYnEjMzaxYnEjMzaxYnEjMzaxYnEjMzaxYnEjMzaxYnEjMza5ZOcUGipFXAon08vD+dbwkW17lzcJ07h+bUeXhEDGhsp06RSJpD0vRsruzsSFznzsF17hxao87u2jIzs2ZxIjEzs2ZxImnctfkOIA9c587Bde4ccl5nj5GYmVmzuEViZmbN4kRiZmbN4kRSD0knS3pF0nxJV+Y7nqaSNFTSE5LmSZor6d/T8r6SHpX0WvrcJ+OYq9L6viLp/Rnlh0manW77uSSl5d0k3ZmWPyepsrXrWZukQkn/lHR/+r6j17dM0h8kvZz+rad0gjpfnv6bniPpdkndO2KdJV0vaWV6N9iaslapp6QL0894TdKFjQYbEX7UepDc3vd1YCTJPeFnktwCOO+xNaEOg4FD09clwKvAOOC/gCvT8iuBH6avx6X17AaMSOtfmG57HpgCCHgQOCUt/wxwTfr6LODONlDvLwK3Afen7zt6fW8CPpm+7gqUdeQ6A0OAN4Ae6fu7gI93xDoD7wYOBeZklOW8nkBfYEH63Cd93afBWPP9H0JbfKRf+sMZ768Crsp3XM2s05+Bk4BXgMFp2WDglbrqCDycfg+DgZczys8GfpO5T/q6iOTqWeWxjhXA48AJ7EkkHbm+vUl+VFWrvCPXeQiwOP2RKwLuB97XUesMVLJ3Isl5PTP3Sbf9Bji7oTjdtVW3mn+sNarSsnYpbbK+C3gOGBQRywDS54HpbvXVeUj6unb5XsdExE5gHdAvF3XI0s+ALwPVGWUdub4jgVXADWl33nWSetKB6xwRS4AfAW8Cy4B1EfEIHbjOtbRGPZv8++dEUjfVUdYu50lL6gX8EfhCRKxvaNc6yqKB8oaOaXWSTgVWRsSMbA+po6zd1DdVRNL18euIeBewiaS7oz7tvs7pmMAZJN03+wE9JZ3X0CF1lLWrOmepJevZ5Po7kdStChia8b4CWJqnWPaZpC4kSeTWiPhTWrxC0uB0+2BgZVpeX52r0te1y/c6RlIRUAqsafmaZOVo4HRJC4E7gBMk/Z6OW9+aeKoi4rn0/R9IEktHrvN7gTciYlVE7AD+BBxFx65zptaoZ5N//5xI6jYNGCVphKSuJANR9+U5piZJZ2b8DpgXET/J2HQfUDML40KSsZOa8rPSmRwjgFHA82nzeYOkI9NzXlDrmJpzfQT4a6Sdqq0tIq6KiIqIqCT5e/01Is6jg9YXICKWA4sljUmLTgReogPXmaRL60hJxWmsJwLz6Nh1ztQa9XwYeJ+kPmkL8H1pWf3yMYDUHh7Av5DMdHod+Fq+49mH+I8haY7OAl5MH/9C0gf6OPBa+tw345ivpfV9hXRmR1o+CZiTbvsFe1ZE6A7cDcwnmRkyMt/1TuM6nj2D7R26vsAhwPT073wvySybjl7nbwEvp/HeQjJTqcPVGbidZBxoB0kr4eLWqifwibR8PnBRY7F6iRQzM2sWd22ZmVmzOJGYmVmzOJGYmVmzOJGYmVmzOJGYmVmzOJFYmyEpJP044/0Vkr7ZQue+UdJHWuJcjXzOR5WswvtERtkESS+mjzWS3khfPyZpP0l/yEEcgyTdL2mmpJckPZCWV0o6p6U/r4E4Pi6pWtLEjLI5+VpR13LDicTakm3AhyT1z3cgmSQVNmH3i4HPRMR7agoiYnZEHBIRh5BcBPal9P17I2JpROQiwX0beDQiDo6IcexZOqUSaLVEkqoiucbBOignEmtLdpLcX/ry2htqtygkbUyfj5f0lKS7JL0q6QeSzpX0fHoPhv0zTvNeSX9L9zs1Pb5Q0n9LmiZplqR/yzjvE5JuA2bXEc/Z6fnnSPphWvZ1kgtBr5H039lUOG0hzElff1zSvZL+krZaLpP0RSULMj4rqW+63/6SHpI0I63P2DpOPZiMxfoiYlb68gfAsWmL6PJG6v+0pHvSFs01kgrS/W9M6z1b0jv+VnW4Hzgo4wp862CK8h2AWS2/BGZJ+q8mHHMwcCDJOkELgOsiYrKSm3l9DvhCul8lcBywP/CEpANIloxYFxGHS+oGPCPpkXT/ycD4iHgj88Mk7Qf8EDgMeBt4RNKZEfFtSScAV0TE9CbXPDGeZKXm7iRXFX8lIt4l6adprD8jSbaXRsRrko4AfkWydH6mXwJ3SroMeAy4ISKWkrRMroiImkT6qUbqPw5YBDwEfIhk2fohETE+Pb4sizpVk9xH46vsWZLDOhC3SKxNiWSF4puBzzfhsGkRsSwitpEsA1HzQzibJHnUuCsiqiPiNZKEM5ZkHaELJL1Issx+P5J1iiBZq2ivJJI6HHgykoUDdwK3ktyEqCU8EREbImIVybLef8msi5LVnI8C7k5j/g1J62MvEfEwyTLzv03r+U9JA+r4vMbqvyAidpEs13EMyfc2UtL/SjoZaGhF6Uy3kayRNSLL/a0dcYvE2qKfAS8AN2SU7ST9H5908bmuGdu2Zbyuznhfzd7/xmuvB1SzZPbn0h/e3SQdT7Ise13qWma7pTRWlwJgbTre0qCIWEPyA36bklsPvxt4q9ZuDdX/Hd9XRLwt6WDg/cBngY+RrMvUWCw704kUX2lsX2t/3CKxNif9AbyLZOC6xkKSriRI7kfRZR9O/dG0n39/kv9bf4VkVdNPK1lyH0mjldwcqiHPAcdJ6p8OxJ8NPLUP8TRZ2mJ7Q9JHIUmq6Q/7XiSdIKk4fV1C0p33JrCB5NbLNRqq/2QlK2AXAP8K/D2dCFEQEX8E/h/JsvWk4zmXNRL+jSTLwNfVMrJ2zInE2qofA5mzt35L8uP9PHAE9bcWGvIKyQ/+gyRjDFuB60iWXn8hHfT+DY201CNZmvsq4AmS+2S/EBF/buiYFnYucLGkmcBcksRa22HAdEmzgKkk40bTSFYJ3qlkWvDlNFz/qSSD83NIxkbuIblT3pNpV9iNJN8DJN1ntVs7e4mI7cDP2XNXP+sgvPqvmb1D2rW1e1A+i/3vBz6UJgvrZDxGYmbNlm3CsY7JLRIzM2sWj5GYmVmzOJGYmVmzOJGYmVmzOJGYmVmzOJGYmVmz/H9yhqdVwMqliQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "error = [] #Create empty array\n",
    "N = 30\n",
    "x = np.logspace(2,5,N-20) #Set x-scale to log space to better see the convergence curve\n",
    "for i in range(20,N): \n",
    "    error.append(np.abs((Temp_num[i]-Temp_ana[i])/Temp_ana[i])) #add each step size of error to error array\n",
    "plt.plot(x,error);\n",
    "plt.title('Relative Error Accumulation');\n",
    "plt.ylabel('Relative Error')\n",
    "plt.xlabel('Number of Time Steps, N');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Temperature approaches 65 degF as time approaches infinity\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xV9f3H8dc7A0jCDAmbsBWiIELEhXtUcYBWrRtXLW2dP1vrarXWVqtWa1sc1F0HLlRcKEVxi4QhG9kQZth75vP745zg9TbjBnJzE/J5Ph73cc/6nvM5d33u+Z5zvl+ZGc4551yskhIdgHPOuZrFE4dzzrkK8cThnHOuQjxxOOecqxBPHM455yrEE4dzzrkK8cRRTUk6StLMSljPfEknVkZMlbVeSTmSNkpKruy49gWSjpVUkOg4KkqSSepcyryLJH1U1TG5+Ki1iSP84VsuKSNi2lWSRicwrN3M7HMz2z+e25DURtIbklZKWidpsqTL4rCdHyUZM1toZvXNbFdlbyvcniTNlTQtHuvfl1TVa2VmL5rZyZWxrnIS1G3hn5KNkrZK2hUxPrUytp8okgZJ+m+i44BanDhCKcD1iQ4imqSUKtrUf4BFQDugKXApsLyKth1PRwPNgI6SDonXRqrwfYqnKnmtqoqZ/SX8U1IfGAR8XTxuZgckOr7SVMVnqTK3UdsTxwPAbyQ1jp4hqX34zyYlYtpoSVeFw5dJ+lLSw5LWhv/ajginL5K0QtLAiLJ1JT0oaWF4pPO4pLRw3rGSCiT9TtIy4Jno6gpJbSUNk1QoaZWkf4XTO0n6OJy2UtKLJe1PKQ4BnjWzTWa208wmmNkHEds8U9LUcP9GS+pW0kokPSvpnojx3bFL+g+QA7wT/uu7Ofq1ldRK0nBJqyXNlvTziHXdJelVSc9L2hDGk1fOfg0E3gbeD4cjYx0t6V5J34ZHWW9LygznFcd1taQlkpZKuikqltclvSBpPXBZ+L7+PVx+SThcN1y+iaR3w/dsTTjcJmJ9mZKeCcutkfRWVKw3hZ+jpZIuj5he1mcpK9zO2vD1/FxSWd/z8l6reyR9Fb5370hqGn7G1ksaK6l91Pr6hd+FlZIeKN52+L34ImLdXSWNDGOcKem8iHnPShos6b3wPR8jqVM477Nwse/CmH5Wxr6VSNKB4XdmjaTpkgZEzBsq6ZEwtk3ha9BM0qPhazpVUveI5ZeFn+kZ4b4MKX7/w/lnSZoUlv1cUm5U2d8oOBJaH077g6R54X5PkXRaOP1g4O/AseF+LwunfyPp4oh17j4qkVQv/Dz/UtIcYEp5+x8zM6uVD2A+cCIwDLgnnHYVMDocbg8YkBJRZjRwVTh8GbATuBxIBu4BFgKDgbrAycAGoH64/N+B4UAm0AB4B7g3nHdsuK6/hmXTwmkF4fxk4DvgYSADqAf0Ded1Bk4Ky2UDnwF/j97PUl6D/wJfAucDOVHz9gM2hetOBW4GZgN1otcLPFv8GkbsT0FpMUS/tsCnwKPhfvUECoETwnl3AVuBfuHrcC/wTRnvazrBl7Af8FNgZXHMEe/hYuDA8LV8A3ghKq6Xw3ndw1hOjIhlBzCA4E9XGnA38A3Bv/Zs4CvgT+HyTcMY0sP3/DXgrYhY3gNeAZqEr/ExUZ+Hu8Pp/YDNQJMYPkv3Ao+H5VKBowDtxWs1G+gENAKmAd8TfG9SgOeBZyKWN+CTMK6ccNnI78sX4XAGwZHu5eF6eoXbPiDi87Qa6BPOfxEYGrWdzjF8x3dvM2JaQ2ApcBHB5+mQcFudw/lDgWXAQeH7+wUwF/hZuPwDwAcR61sGTABahe//WOCOcN5h4bZ6h2WvDl+TlIiyY8OyaeG0nwEtCT5flxD8hmSF8wYB/43an2+AiyPGdy9D8H0ygs9Z43B/ytz/mH8/E/0DnqgHPySOA4F14Zte0cQxK2Je93D55hHTVhH8EIrgR7hTxLzDgXkRPxTbgXoR84/lh8RxOMEPWEoM+zUAmBC9n6Us2wS4D5gK7AImAoeE834PvBqxbBLBD+6x0etlLxIH0DbcdoOI+fcSHAlB8GP934h5ucCWMvb/4uLXiiCZrgXOinoP74ta3/bwS1QcV9eI+fcDT0XE8lnU9uYA/SLGfwLMLyW2nsCacLglUESYDKKWOxbYwo8/eysIfojK+yzdTXAEEcsPayyv1e0R43/jxz+aZwATI8YNOCVi/FfAqIjvS3Hi+BnweVQsTwB3RnyenoyY1w+YEbWdPU0cA4GRUdOeA34XDg8F/hkx77f8+Pt0CLAsYnwZcFnE+NnA1HD4mcjXL5y2ADg0ouyF5ezDDOAn4fCeJo4jYt3/WB+1vaoKM5sCvAvcsgfFI88HbAnXFz2tPkFSSgfGhYesa4ER4fRihWa2tZTttAUWmNnO6BnhYfRQSYvD6pMXgKxYgjezNWZ2iwV1v80JEsdbkkTwL2hBxLJFBP8SW8ey7gpoBaw2sw0R0xZEbWdZxPBmoJ5Kr68dSJDwdprZNoIjyoFRyyyK2lYqP37Noue3KmVecfwLIsZ3Ly8pXdITkhaE781nQGMFV5O1JdjvNaXsx6qo93szsX2WHiA4SvgorDIq63Mdy2sV/Xku6fMdqazXrlg74NDi+MN9uAhoEbFM9HsevZ091Q44OmrbPyVI5MUqa5/bAbdFbSubH3+2f/R5knRlRNXWWoIahZi+z2WI3EYs+1+ufeHkXmW4ExhP8I+q2KbwufhwHn78wa6IlQQfuAPMbHEpy1gZ5RcBOZJSSkge94Zle5jZqrC+8l8VDdDMVkp6kOCHIxNYQnAUBQRX3xD82JUU/yaC16lY9OtU1r4tATIlNYhIHjmlbKdMCs4fHA/0kfTTcHI6QaLJMrOV4bS2EcVyCKqfVkZMb0vwT694/pIy9mUJwZdxagnL3wTsT/APc5mkngTVGiJ4TzMlNTaztRXYzTI/S+FreBNwk6QDgE8kjTWzUZHLVeC1qqi2lPxaRFoEfGpmJ+3hNvbGIuAjMzujEtcZ/Xkq3udFwHtm9rf/LbLb7s+TpP2AfxK8L9+aWZGkGQSflx8tG6G87150uUrZ/1p/xAFgZrMJ6pqvi5hWSPDjdbGkZElXENT17sn6i4B/Aw9LagYgqbWkn8S4im8J6iXvk5QRnvQ6MpzXANgIrJXUmuDQOiaS/hqeKEuR1AD4JTDbzFYBrwKnSTpBUirBj9E2gjr8aBMJTopmSmoB3BA1fznQsaQYzGxRuM57w/3qAVxJUK9dUZcQ1CHvT1At1JPgXE0BcEHEchdLypWUTlC187r9+NLg34dHCwcQ1MO/UsY2XwbukJQtKQv4A8FRHwTvzRaC9yaT4A9K8X4vBT4AHlVwEj1V0tHl7WB5nyVJp0vqHCb69QTVgCVd9hzra1VRvw33py3BFYslvXbvAvtJuiTc71RJh6iUiy9KUOrnKQZvAQdL+lm43TqSDgt/tPfUdZJahu//Lfywz0OAayXlKVBfwQUn6aWspz5B9WUhkCRpEMERR7HlQNvw+1hsInBO+N3pSlA9V5ZK2X9PHD+4m+CkXaSfE/wQrwIOoOQfzVj9jqAK4Zuw2uK/BF/acoU/amcQfIgWEny5i68m+SPBycV1BCfBhlUgpnTgTYK67bkE/5zPDLc5k6AO/J8E/3LPAM4ws+0lrOc/BCfv5wMf8b8/FvcS/LiulfSbEspfQHB+YUkYz51mNrIC+1FsIPComS2LfBCcLI6sgvkPQT36MoJ64Oui1vMpwXs1CnjQzMq6ce0eIB+YBEwmOHItvsLs7wQnJFcS1EWPiCp7CcHRzgyCcxjRCbc0ZX2WuoTjG4GvCV6P0SWsI9bXqqLeBsYR/KC9BzwVvUB4VHQywUUZSwjeh+ILQ2JxF/Bc+Hk6r7yFo7a9huA81OUEf8aWELxfqWWVK8dQgosCZhF8Bu4Pt/UlwWfrCYLv2PfAhZRyBG5m4wle//wwtg7hcLERBN+xFfrhisv7CWqOCgkS1QuUobL2X+HJEedqBQU3eL5gZk+WMK89MA9ILel8knPRFFwWe46ZfVHuwvsQP+JwzjlXIZ44nHPOVYhXVTnnnKsQP+JwzjlXIbXiPo6srCxr3759osNwzrkaZdy4cSvNLDt6eq1IHO3btyc/P7/8BZ1zzu0maUFJ072qyjnnXIV44nDOOVchnjicc85ViCcO55xzFeKJwznnXIV44nDOOVchnjicc85ViCeOMoyZu4pHR89OdBjOOVeteOIow8hpy3ngw5lMWbwu0aE451y14YmjDNee0IUm6XW4+91peGOQzjkX8MRRhkZpqdx08n58O281H0xZluhwnHOuWvDEUY6f5bWla4sG/OX96WzdUVLXzc45V7t44ihHSnISfzg9l4I1W3jqi3mJDsc55xLOE0cMjuicxcm5zRn8yWyWr9+a6HCccy6hPHHE6PbTurFzl3H/iJmJDsU55xLKE0eM2jXN4PK+7XljfAHfLVqb6HCccy5hPHFUwDXHdSarfl2/PNc5V6t54qiABvVS+e1P9mPcgjUM/25JosNxzrmE8MRRQef0bssBrRry1w9msGW7X57rnKt94po4JJ0iaaak2ZJuKWH+RZImhY+vJB1UXllJmZJGSpoVPjeJ5z5ES04Sd55xAEvWbWXIZ3OrctPOOVctxC1xSEoGBgOnArnABZJyoxabBxxjZj2APwFDYih7CzDKzLoAo8LxKtWnQyandW/J45/OYem6LVW9eeecS6h4HnH0AWab2Vwz2w4MBfpHLmBmX5nZmnD0G6BNDGX7A8+Fw88BA+K4D6W65dSu7DLjrx/MSMTmnXMuYeKZOFoDiyLGC8JppbkS+CCGss3NbClA+NyspJVJulpSvqT8wsLCPQi/bG0z07n6qI68NXEJ4xeuKb+Ac87tI+KZOFTCtBKvYZV0HEHi+F1Fy5bGzIaYWZ6Z5WVnZ1ekaMx+eWwnmjWoyx/fmUZRkV+e65yrHeKZOAqAthHjbYD/uYZVUg/gSaC/ma2KoexySS3Dsi2BFZUcd8wy6qZw8yld+W7RWt6auDhRYTjnXJWKZ+IYC3SR1EFSHeB8YHjkApJygGHAJWb2fYxlhwMDw+GBwNtx3IdynX1waw5q04i/jpjBpm07ExmKc85VibglDjPbCVwDfAhMB141s6mSBkkaFC72B6Ap8KikiZLyyyoblrkPOEnSLOCkcDxhkpLEH87IZfn6bTw2ek4iQ3HOuSqh2tB0Rl5enuXn58d1GzcMncD7k5fx/vVH0blZ/bhuyznnqoKkcWaWFz3d7xyvJLeflku91CRuf3Oyt2PlnNuneeKoJNkN6nJbv26Mmbea1/ILEh2Oc87FjSeOSnReXlv6tM/kz+9PZ+XGbYkOxznn4sITRyVKShJ/OftANm/fyT3vTkt0OM45FxeeOCpZ52YN+OWxnXlr4hI+n1X5d6w751yieeKIg18d24mOWRnc/uYUb3rdObfP8cQRB/VSk/nzWd1ZuHozj4yalehwnHOuUnniiJPDOzXl3N5t+Pfnc5m+dH2iw3HOuUrjiSOObuvXjUZpqdw6bDK7vBFE59w+whNHHDXJqMPvT+/GxEVreXHMgkSH45xzlcITR5wN6Nmao7pkcf+ImSxbtzXR4Tjn3F7zxBFnkrhnwIHs2FXEXcOnll/AOeeqOU8cVaBd0wyuP7ELI6YuY+S05YkOxznn9kpKaTMkDYuh/Gozu6oS49ln/fyojgyfuIQ/vD2Fwzs1pX7dUl9655yr1sr69eoODCpjvoBHKjecfVdqchJ/Obs7P33sK/720UzuPOOARIfknHN7pKzEcaeZjSqrsKQ/V3I8+7ReOU24+NB2PPfVfAb0bM1BbRsnOiTnnKuwss5xvFxeYTN7qaz5kk6RNFPSbEm3lDC/q6SvJW2T9JuI6fuHPQIWP9ZLuiGcd5ekxRHz+pUXZ3Xy21P2J7tBXW5+fRJbd3hzJM65mqesxDGueEDS3yu6YknJwGDgVCAXuEBSbtRiq4HrgAcjJ5rZTDPraWY9gd7AZuDNiEUeLp5vZu9XNLZEalgvlft+2oOZyzfw8Mjvyy/gnHPVTFmJQxHDR+/BuvsAs81srpltB4YC/SMXMLMVZjYW2FHGek4A5pjZPnMH3XH7N+PCQ3MY8vlcxsxdlehwnHOuQspKHHvbRkZrYFHEeEE4raLO53+rza6RNEnS05KalFRI0tWS8iXlFxZWv+bNb+/XjZzMdG567Ts2btuZ6HCccy5mZSWOrpLGS5oQMTxe0gRJ42NYt0qYVqFkJKkOcCbwWsTkx4BOQE9gKfC3ksqa2RAzyzOzvOzs7Ipstkpk1E3hofMOYsnaLfzpHe/0yTlXc5R3Oe7eKADaRoy3AZZUcB2nAuPNbPddc5HDkv4NvLs3QSZS73aZDDqmE4+OnsOJuc05Kbd5okNyzrlylZo4zGzOXq57LNBFUgdgMUGV04UVXMcFRFVTSWppZkvD0bOAKXsZZ0LdcOJ+fDKzkFuHTaJXztE0rV830SE551yZym1yRNIaSaujHvMkvSapfWnlzGwncA3wITAdeNXMpkoaJGlQuO4WkgqA/wPukFQgqWE4Lx04CYi+g/1+SZMlTQKOA26s8F5XI3VSkvj7z3qyfstObntzMmbe/LpzrnpTeT9Uku4GlgMvEZy3OB/IBmYDV5nZcfEOcm/l5eVZfn5+osMo05DP5vCX92fw4LkHcU7vNokOxznnkDTOzPKip8fSyOHJZjbYzNaY2WozexQ41cxeBDIrPdJa6sq+HenTPpM/Dp9KwZrNiQ7HOedKFVPruJLOjhouvmKqKB5B1UbJSeJv5x1EkRm/fW0SRd5joHOumoolcVwM/Dw8t7EK+DlwSXgO4oa4RlfLtM1M584zDuDruat4+st5iQ7HOedKVG7b3mY2m+Cy2JJ8WrnhuHPz2vDRtGXc/+FMjtkvmy7NGyQ6JOec+5FYrqrqLOlDSd+F4z0k3Rr/0GonSdx7dg8a1E3hxlcnsn2n1wY656qXWKqqngT+yA/nMyYTVF+5OMluUJc/n9WdKYvX88+PZyU6HOec+5FYEkeGmX1VPGLB9btlNUroKsEpB7bgp73aMPiT2YxfuCbR4Tjn3G6xJI5V4d3fBiBpALAsrlE5AO48M5eWjdK49qUJrNvsudo5Vz3EkjiuAZ4iaOhwAXAL8Mu4RuWAoO+OwRf1YsWGrdz02nd+V7lzrlooN3GY2WwzOx5oCRxkZoeZmV8rWkV6tm3Mbf268d/py3nyc3/ZnXOJV+rluJKuK2U6AGb2jzjF5KJcdkR7vp23mvtGzODgnMbktfcb9p1ziVPWEUd2+DiC4Ea/TuHjeoK+MFwVkcRfz+lBmyZpXPPSBFZt3JbokJxztVipicPMfm9mvweaAD3N7Hozux44mKDaylWhhvVSGXxhL1Zv3s6Nr37nTZI45xImlpPj7YCtEePbgA7xCceV5cDWjbjzjFw++76Qxz7d2+5SnHNuz5Tb5AhBc+pjJL1BcEnu2cALcY3KlerCPjmMmbuav300k145TTi8U9NEh+Scq2ViuarqbuBqYAvBkccgM7sn3oG5kkniL2d3p31WBtcNnUDhBj/f4ZyrWjE1q25mY83sb+FjbKwrl3SKpJmSZku6pYT5XSV9LWmbpN9EzZsf9vQ3UVJ+xPRMSSMlzQqfm8Qaz76ift0UHr2oFxu27uD6oRPY5ec7nHNVqNTEIenb8gqXtYykZGAwQcu6ucAFknKjFlsNXAc8WMpqjjOznlE9UN0CjDKzLsCocLzW6dqiIXf3P5Cv5qzikVHenpVzruqUdY6ju6TxZcwXUFYFex9gtpnNBZA0FOgPTCtewMxWACsknRZ7yPQHjg2HnwNGA7+rQPl9xnl5bfl23mr++fEs8to14ej9shMdknOuFigrcRwYQ/mdZcxrDSyKGC8ADo0lqJABH0ky4AkzGxJOb25mSwHMbKmkZiUVlnQ1wbkZcnJyKrDZmuVP/Q9kcsE6bnxlIu9ddxQtGtVLdEjOuX1cWfdxzInhsaCMdauEaRWpjD/SzHoRVHX9WtLRFSiLmQ0xszwzy8vO3nf/iafVSWbwRb3YsmMX1748np27vP8O51x8xXRyfA8VAG0jxtsAS2ItbGZLwucVwJsEVV8AyyW1BAifV1RKtDVY52b1uffs7oydv4Y/vz890eE45/Zx8UwcY4EukjpIqgOcDwyPpaCkDEkNioeBk4Ep4ezhwMBweCDwdqVGXUP179may49szzNfzmfotwsTHY5zbh8Wyw2ASGoDdDGzTyTVBVLMbFNZZcxsp6RrgA+BZOBpM5sqaVA4/3FJLYB8oCFQJOkGgiuwsoA3wwYVU4CXzGxEuOr7gFclXQksBM6t2C7vu27v1405hZv4/dtT6JCVwaEd/eZA51zlU3l9PEi6gqBPjkZm1knSfsCjZnZiVQRYGfLy8iw/P7/8BfcB67bs4KxHv2TNpu0Mv6YvbTPTEx2Sc66GkjQu6nYIILaqquuAw4D1AGb2PVDilUwu8RqlpfLUwEMoMrjyubFs2Oo9BzrnKlcsiWOrmW0vHglv7CvpiilXTXTIyuDRi3oxp3ATNwyd6HeWO+cqVSyJ40tJNwP1JB0HvAK8G9+w3N46snMWd52Ry6gZK7h/xIxEh+Oc24fEkjhuBjYAMwg6cRoF3B7PoFzluOTw9lxyWDue+Gwur48rSHQ4zrl9RJlXVYXVUk+b2UDgsaoJyVWmP5yRy9yVG7lt2GQ6ZKXTu513O+uc2ztlHnGY2S6gpaTUKorHVbLU5CQGX9iLVo3r8Yv/jKNgzeZEh+Scq+FiqaqaC3wu6VZJ1xU/4h2YqzyN0+vw5MBD2LaziKuey2fTtrKaGHPOubLFkjgKgZFAOpAd8XA1SOdm9fnXhb34fvkGbnxlovdZ7pzbY+XeOW5mv6+KQFz8HbNfNneclsvd707joZHf85uf7J/okJxzNVC5iUPSSEpo1dbMTo5LRC6uLj+yPd8v38C/PplNTtN0zstrW34h55yLEEtbVXdEDNcDfgp4R9c1lCTu7n8gi9du4dZhk8lMr8OJuc0THZZzrgYp9xyHmY2JeHxqZtfxQxPnrgaqk5LEYxf35oBWDfn1S+MZO391okNyztUg5SYOSQ0jHo0lnQC0rILYXBzVr5vCM5cdQuvGaVz57FhmLFuf6JCcczVELFdVTSXoC2MqMIHgrvGfxzMoVzWa1q/Lc1f0Ia1OMgOf/tbv8XDOxSSWxNHRzHLMrK2ZdTCz44Ev4x2YqxptM9N57oo+bNm+i0uf+pZVG/30lXOubLEkjjElTPu2sgNxidO1RUOeuuwQFq/dwhXPjvUbBJ1zZSo1cUhqJukgIE1Sd0k9wkdfgpsByyXpFEkzJc2WdEsJ87tK+lrSNkm/iZjeVtInkqZLmirp+oh5d0laLGli+OhXsV12JTmkfSaDL+zFlCXrGfTCOLbvLEp0SM65aqqsy3FPA64A2gCPRkzfAJR7U2DYQOJg4CSgABgrabiZTYtYbDVBR1EDoorvBG4ys/Fh3+PjJI2MKPuwmT1YXgyuYk7Mbc69Z3fn5tcn8ZvXvuPvP+tJUpJ3veKc+7FSE4eZPQM8I+k8M3t1D9bdB5htZnMBJA0F+gO7E4eZrQBWSDotattLgaXh8AZJ04HWkWVdfJyX15ZVG7fz1xEzyMyow51n5BL2/e6cc0BsTY68KuknwAEENwAWT/9LOUVbA4sixguAQysaoKT2wMH8+FzLNZIuBfIJjkzWlFDuauBqgJycnIputlYbdExHVm7cxlNfzCO7QV1+fVznRIfknKtGYrmP41FgIPB/QBpwMRDLL0lJf1Mr1LKepPrAG8ANZlZ8o8FjQCegJ8FRyd9KKmtmQ8wsz8zysrO9TcaKkMTt/bpx1sGteeDDmQz9dmGiQ3LOVSOxXFXV18wuBFaFDR4eSnDeozwFQGRDSG2AJbEGFvYB8gbwopkNK55uZsvNbJeZFQH/xu9ij4ukJHH/OT04Zr9sbntzMm9PXJzokJxz1UQsiWNr8bOkFuF4+xjKjQW6SOogqQ5wPjA8lqAUVKo/BUw3s4ei5kXetX4Wwc2JLg5Sk5N47OJeHNI+kxtfmcjw72LO+865fVgsjRy+L6kx8CAwEdgFPFdeITPbKeka4EOguAvaqZIGhfMfDxNRPtAQKJJ0A5AL9AAuASZLmhiu8jYzex+4X1JPgmqv+cAvYt5bV2HpdVJ45vJDuOyZsdwwdAIAZx7UKsFROecSSWaln3aQlAQcYmZjwvE0IM3MalSreHl5eZafn5/oMGq0zdt3ctkzY8mfv5pHzj+YMzx5OLfPkzTOzPKip5fX53gR8EjE+JaaljRc5Uivk8Kzlx9CXvtMrh86gXe82sq5WiuWcxwjJfWPeySu2kuvE7So68nDudotlsRxDfCmpC2SVktaI8mPOmqpjLo/JI8bXpnoycO5WiiWxJEFpAL1gexw3G+MqMWKk0fvdk08eThXC8XSA+Au4Fzgd+FwS4Kb71wttjt55ATJ491Jnjycqy1iuXP8X8BxBJfHAmwGHo9nUK5myKgbXKrbO6cJ1w/15OFcbRFLVdURZvYLwhsBw6uq6sQ1KldjePJwrvaJJXHsCO/nMABJTQHvrMHtFpk8rnt5Aq+OXVR+IedcjRVL4hhM0GZUtqQ/Al8Af41rVK7GyaibwrNXHELfLtnc/MYkHh09m7JuLnXO1VyxNKv+vKRxwInhpHPNzNuHcv8jvU4KT16ax29f/477R8xk1cbt3N6vm3cG5dw+Jpa2qiBoa2oHQXVVLEcprpaqk5LEw+f1JDOjDk99MY/Vm7Zz/zk9SE32j41z+4pYrqq6HXgZaEXQNPpLkm6Nd2Cu5kpKEn84PZff/mR/3pywmJ8/n8/m7TsTHZZzrpLE8jfwYoKGDu8ws9sJ+r+4NL5huZpOEr8+rjP3nd2dz74v5KInx7Bm0/ZEh+WcqwSxJI4F/LhKKwWYG59w3L7m/D45PHZxb6YuWc+5T3zNkrVbEh2Sc24vxZI4NgNTJT0p6d/AZGCtpIckPVROWef4yQEteP6KPixft5VzHvuK2Ss2JDok59xeiCVxvAfcBXwNfAPcDXwMTN2L2LMAABx5SURBVA0fzpXrsI5NGfqLw9i+yzjn8a+ZsHBNokNyzu2hWNqqegr4D/CpmT0V/SirrKRTJM2UNFvSLSXM7yrpa0nbJP0mlrKSMiWNlDQrfG4S++66RDqgVSOG/fIIGqWlcuG/xzB65opEh+Sc2wOxXFV1GkH11MhwvKekN2Mol0xw8+CpBN3BXiApN2qx1cB1BN3Sxlr2FmCUmXUBRoXjrobIaZrO64OOoGN2Blc+l89/vp6f6JCccxUUS1XV3cChwFoAM5sIdI6hXB9gtpnNNbPtwFDgRx1CmdkKMxtLcI9IrGX780Of588BA2KIxVUj2Q3q8sovDue4/bP5/dtTueOtyezY5a3YOFdTxNRWlZmtjZoWS1sSrYHIRosKwmmxKKtsczNbChA+NytpBZKulpQvKb+wsDDGzbqqUr9uCk9ckscvjunIC98s5LJnvmXtZr9c17maIJbEMV3SeUCSpA6S/k5wkrw8JbUzEWvjRXtTNljYbIiZ5ZlZXna29ztVHSUniVtP7caD5x7E2HlrGDD4S2av2JjosJxz5Yi169jeBC3ivglsA26IoVwB0DZivA0Qa5vbZZVdLqklQPjsZ1hruHN6t+Hlqw9lw9adnPXol3z2vR8hOledxXJV1SYz+52ZHWxmPcPhzTGseyzQJTxKqQOcDwyPMa6yyg4HBobDA4G3Y1ynq8Z6t8vk7WuOpHXjNC575lue/XKet67rXDVVbiOHknoRXLnUPnJ5M+tVVjkz2ynpGuBDgkYSnzazqZIGhfMfl9QCyAcaAkWSbgByzWx9SWXDVd8HvCrpSmAhQbe2bh/Qpkk6b/zyCG54ZSJ3vTON71ds5I9nHuANJDpXzai8f3WSZgC3EVySu/vSFzObE9/QKk9eXp7l5+cnOgwXo6Ii44GPZvLY6Dkc1jGTxy7qTZMM73TSuaomaZyZ5UVPj+Wv3CozG2Zms8xsTvEjDjE6BwSt6/7ulK48dN5BjF+wlgGPfsms5d5MiXPVRSyJ44+SnpB0rqQzix9xj8zVemf3asPLVx/Gpm07GTD4S96euDjRITnniC1xXERwQ94AgvMJ5wLnxDMo54r1bteEd67tS7eWDbl+6ERuf3MyW3fsSnRYztVqsfQA2NvMDox7JM6VomWjNF6++jAe/GgmT3w6l4mL1vLoRb1o1zQj0aE5VyvFcsQxRtL+cY/EuTKkJidx66ndePLSPArWbOH0f3zBiClLEx2Wc7VSLImjDzBJ0lRJ4yVNkDQ+3oE5V5ITc5vz7rV96ZidwaAXxnP3O9PYvtPbuXKuKsVSVeWNCLpqpW1mOq8NOoK/vD+dp7+cx/iFaxh8US9aN05LdGjO1Qqx3Dk+B8gGjgyH1/K/rdk6V6XqpCRx15kH8OhFvZi9YiOn/eNzPp6xPNFhOVcrxNIfxx3AncAd4aR6wEvxDMq5WPXr3pJ3r+1Lq0ZpXPFsPn8dMYOd3kS7c3EVyzmOc4B+wCYAM1tM0ESIc9VC+6wMhv3qCC7ok8Njo+dw3hNfM3/lpkSH5dw+K5bEsc2CdkkMQFJ6fENyruLqpSZz79nd+ccFBzN7xUZOfeRzXvhmgTeU6FwcxJI4hkkaDDSSdDnwEfB0fMNybs+ceVArPrzxaPLaN+GOt6Zw+bNjWbF+a6LDcm6fUm4jhwCSTgVOJuhg6UMz+yDegVUmb+Sw9ikqMl4Ys4C/vD+deqnJ3DPgQE7v0SrRYTlXo5TWyGGpiUPSR2Z2ctwjqwKeOGqvOYUb+b9Xv+O7RWvp37MVd595II3SUxMdlnM1wp60juv9rboar1N2fd4YdDj/d9J+vDdpKT/5+2d8Pst7GHRub5SVOBpJOru0Rywrl3SKpJmSZku6pYT5kvSPcP6ksNMoJO0vaWLEY33YyROS7pK0OGJevz3ac1drpCQncd0JXRj2qyPIqJvMJU99y51vT2HLdm8s0bk9Udad442A0wnOa0QzYFhZK5aUDAwGTiLoQ3yspOFmNi1isVOBLuHjUOAx4FAzmwn0jFjPYoL+zos9bGYPlrV956L1aNOY9647ivtHzOTpL+fx+ayVPHBuD3q3y0x0aM7VKGUljgVmdsVerLsPMNvM5gJIGgr0ByITR3/g+fBy328kNZbU0swiW687AZhjZgv2IhbngOCy3T+ckcuJ3Zrx29cncc7jX3NhnxxuPqUrjdL83IdzsSirqqqkI42KaA0sihgvCKdVdJnzgZejpl0TVm09LalJSRuXdLWkfEn5hYVep+1+7IjOWXx049FcfkQHXv52ISc+9CnvTVrq9304F4OyEscle7nu0qq4Yl5GUh3gTOC1iPmPAZ0IqrKWAn8raeNmNsTM8swsLzvbz/O7/5VRN4U/nJHL27/uS7MGdfn1S+O58rl8CtZsTnRozlVrpSYOM5uyl+suANpGjLcBllRwmVOB8Wa2u/U6M1tuZrvMrAj4N0GVmHN7rHubRrz96yO547RufD1nFSc99BlPfj7X27xyrhSx3Dm+p8YCXSR1CI8czgeGRy0zHLg0vLrqMGBd1PmNC4iqppLUMmL0LGBvE5xzpCQncdVRHRn5f0dzeKem3PPedAY8+iWTC9YlOjTnqp24JQ4z2wlcA3wITAdeNbOpkgZJGhQu9j4wF5hNcPTwq+LyYZtYJ/G/V2/dL2mypEnAccCN8doHV/u0aZLOUwPzGHxhL5av30b/wV9w9zvT2LhtZ6JDc67aKLfJEUldgHuBXIIm1QEws47xDa3y+J3jbk+s27KD+0fM4MUxC2nVqB63n5ZLv+4tkPb2uhHnaoY9uXO82DMEJ6R3EvzDfx74T+WG51z10ygtlT+f1Z03fnk4DdNS+fVL4/nZE98wZbFXX7naLZbEkWZmowiOThaY2V3A8fENy7nqo3e7TN69ti9/PutAZhdu5Ix/fcHNr3/Hig3e6q6rnWLpc3yrpCRglqRrCO7ibhbfsJyrXlKSk7jo0Hac3qMV//p4Fs9+NZ/3Ji3l18d35oojO1AvNTnRITpXZWI54rgBSAeuA3oDFwOXxjMo56qrRmmp3H5aLh/deAyHd8ri/hEzOenhT/lgst886GqPWBJHezPbaGYFZna5mf0UyIl3YM5VZx2yMnhyYB4vXHko6akp/PLF8Zw/xM9/uNohlsRxa4zTnKt1+nbJ4r3r+nLPgAOZtSI4/3HLG5NYts7Pf7h9V6nnOMJe//oBrSX9I2JWQ4IrrJxzBOc/Lj6sHWcc1Ip/jprFc1/P580Ji7n08Hb88tjOZGbUSXSIzlWqso44lgDjgK3hc/FjOPCT+IfmXM3SKC2VO07P5eObjuX0Hq146ot5HPXXj3lo5Pes37oj0eE5V2liuQEwJbwLvMbyGwBdIsxesYGHRn7P+5OX0Tg9lUHHdGLg4e1Jq+NXYLmaYU/6HJ/M/7Zmu5uZ9ai88OLLE4dLpCmL1/HgRzMZPbOQ7AZ1ufb4zpx/SA51UuLZVJxze29PEke7slZYkzpW8sThqoOx81fzwIcz+Xbealo3TuP6E7tw9sGtSUn2BOKqpwonjqjC7YAuZvZfSWlAipltiEOcceGJw1UXZhZ0WfvhTCYvXkfH7Ax+dWxn+vdsRaonEFfN7HFbVZJ+DrwOPBFOagO8VbnhOVc7SOLo/bIZfs2RPH5xb+okJ/Gb177j2AdG89xX89myfVeiQ3SuXLGcHJ9I0FnSGDM7OJw22cy6V0F8lcKPOFx1ZWaMnlnI4E9mk79gDU0z6nBF3w5ccng7GtbzPtBdYpV2xBFLW1XbzGx7cVPSklIo46S5cy52kjiuazOO69qMb+etZvAns3ngw5k8PnoOlxzejiv6diCrft1Eh+ncj8SSOD6VdBuQJukkgs6W3olvWM7VPn06ZNKnQx+mLF7HY6Pn8Ninc3jqi3mcf0hbfn50R9o0SU90iM4BsVVVJQFXAicDIujR70mL4ay6pFOAR4DksMx9UfMVzu8HbAYuM7Px4bz5wAZgF7Cz+HBJUibwCtAemA+cZ2ZryorDq6pcTTSncCNPfDqHYeMXA3DmQa24om8HDmzdKMGRudpib6+qygYws8IKbDAZ+J6g+9cCgj7ILzCzaRHL9AOuJUgchwKPmNmh4bz5QJ6ZrYxa7/3AajO7T9ItQBMz+11ZsXjicDXZkrVb+Pfnc3ll7CI2b9/FIe2bcPmRHTg5t7lfyuviqsJXVSlwl6SVwAxgpqRCSX+IcZt9gNlmNtfMtgNDgf5Ry/QHnrfAN0BjSS3LWW9/4Llw+DlgQIzxOFcjtWqcxp1nHMDXt57AHad1Y9n6rfzqxfEcff8nPDZ6Dms3b090iK6WKevvyg3AkcAhZtbUzDIJjgqOlHRjDOtuDSyKGC8Ip8W6jAEfSRon6eqIZZqb2VKA8LnETqUkXS0pX1J+YWHMB0rOVVuN0lK56qiOjP7NcQy5pDftmmbw1xEzOOzeUdw6bDIzl9WYW6tcDVfWyfFLgZMiq4rMbK6ki4GPgIfLWbdKmBZdL1bWMkea2RJJzYCRkmaY2WflbPOHlZgNAYZAUFUVaznnqrvkJHHyAS04+YAWzFi2nme/nM+w8QW8/O1CjuzclMuP6MBxXZuRnFTS18u5vVfWEUdq9PkF2H2eI5YLzAuAthHjbQha3I1pGTMrfl4BvElQ9QWwvLg6K3xeEUMszu2TurZoyH0/7cE3t57Azafsz9zCTVz1fD7HPPAJ/xw1y/sFcXFRVuIoq+I0lkrVsUAXSR0k1QHOJ2iSPdJw4NLwfMphwDozWyopQ1IDAEkZBFd0TYkoMzAcHgi8HUMszu3TmmTU4VfHduazm4/jXxceTE5mOn8b+T1H3DeKq54by8hpy9m5qyjRYbp9RFlVVQdJWl/CdAH1yluxme2UdA3B5bvJwNNmNlXSoHD+48D7BFdUzSa4HPfysHhz4M3wpsMU4CUzGxHOuw94VdKVwELg3PJica62SE1O4vQerTi9Ryvmr9zEq/mLeG1cAf+dnk/zhnU5t3dbfnZIW9pm+j0hbs/FdDluTeeX47rabMeuIj6esYJXxi5i9MwVFBn07ZzF+X3aclJuc+qmeP8grmR7dR9HTeeJw7nAkrVbeH1cAa+MXcTitVvIzKjDgJ6tGXBwK7q3bkRx00LOgScOTxzORdhVZHwxeyWvjF3If6etYPuuIjpmZ3BWz9b079manKZeleU8cXjicK4U6zbv4IMpS3lzwmLGzFsNQO92TRhwcGtO696SzIw6CY7QJYonDk8czpVr8dotDJ+4hLcmLGbm8g2kJIlj98+mf8/WnJTbnHqpfj6kNvHE4YnDuQqZvnQ9b01YzNsTl7Bs/Vbq103hpNzmnHJgC47ZL9uTSC3gicMTh3N7ZFeRMWbeKt6asJiPpi1n7eYdpNdJ5riuzTj1wBYct38zMurG0kODq2k8cXjicG6v7dhVxJi5q3l/ylI+mrqMlRu3UzcliWP2y+bU7i04oVtz77lwH+KJwxOHc5VqV5GRP381H0xZxogpy1i2fiupyaJv5yxOPbAlx3dr5r0X1nCeODxxOBc3RUXGhEVrGTFlKe9PXsbitVuQoGfbxpzQtRnHd21Ot5YN/D6RGsYThycO56qEmTF1yXpGTV/BxzOW813BOgBaNqrH8V2bcXzXZhzRKYu0On5yvbrzxOGJw7mEWLFhK6NnFDJqxnI+n7WSzdt3UTcliSM7Z+1OJK0apyU6TFcCTxyeOJxLuG07dzFm7mo+nrGCUTOWs2j1FgA6N6tP385ZHNUli8M6NvWrtKoJTxyeOJyrVsyMOYUb+XjGCj6ftZJv561m284iUpJEr5wm9O2SRd8uWfRo3cj7Vk8QTxyeOJyr1rbu2MW4BWv4fNZKvphdyJTFQa8ODeqlcESnpvTtks1RnbNo1zTdT7JXkdIShx8POueqhXqpyRzZOYsjO2cBXVm9aTtfzl7JF7NW8sXslXw4dTkALRrW49COmRzaoSmHdsykY1aGJ5Iq5kcczrlqz8yYt3ITX85ZxZi5qxgzbzWFG7YBkFW/Lod2yNydTLo0q0+S97deKRJyxCHpFOARgh4AnzSz+6LmK5zfj6AHwMvMbLyktsDzQAugCBhiZo+EZe4Cfg4Uhqu5zczej+d+OOcSSxIds+vTMbs+lxzWbnciGTNv9e5E8t7kpQA0SU+lT4dM+nRoSq+cxhzQqhF1UvwcSWWKW+KQlAwMBk4CCoCxkoab2bSIxU4FuoSPQ4HHwuedwE1hEmkAjJM0MqLsw2b2YLxid85Vb5GJ5II+OZgZBWu28E2YRMbMW7W7aqtOShI9WjeiV7sm9MppTK+cJjRrWG7v164M8Tzi6APMNrO5AJKGAv2ByMTRH3jegvqybyQ1ltTSzJYCSwHMbIOk6UDrqLLOOQcEiaRtZjptM9M5N68tAMvWbWX8wjWMX7CG8QvX8OyX8xnyWREArRun/SiR5LZqSKpfuRWzeCaO1sCiiPECgqOJ8pZpTZg0ACS1Bw4GxkQsd42kS4F8giOTNdEbl3Q1cDVATk7Onu6Dc66GatGoHv26t6Rf95ZAcA/J1CXrGb9gDRMWriV//mre+W4JEByV5LZsSI82jejeuhE92jSmU3aGXwZcingmjpLOTkWfiS9zGUn1gTeAG8xsfTj5MeBP4XJ/Av4GXPE/KzEbAgyB4OR4RYN3zu1b6qYk0yunCb1ymuyetnTdFsYvWMvERWuYvHgdw8Yv5vmvFwCQlprMAa0a0r1NozChNKZjVoafeCe+iaMAaBsx3gZYEusyklIJksaLZjaseAEzW148LOnfwLuVG7ZzrrZo2SiN03qkcVqP4KikqMiYt2oTkwvWMalgHZMXr2Xot4t45sv5AGTUSSa3VUNyWzakW/jYv0WDWtepVTwTx1igi6QOwGLgfODCqGWGE1Q7DSWoxlpnZkvDq62eAqab2UORBSLOgQCcBUyJ4z4452qRpCTRKbs+nbLrM+Dg1kDQfPycwo1MKljHpIK1TF2yntfHFbBp+66gjKBDVga5rRrRrWUDurUMEkuzBnX32ftL4pY4zGynpGuADwkux33azKZKGhTOfxx4n+BS3NkEl+NeHhY/ErgEmCxpYjit+LLb+yX1JKiqmg/8Il774JxzyUliv+YN2K95A87p3QYIjkwWrdnMtCXrmb50PdOWbmD8gjW7z5kAZGbUYf/mDdiveX26hOX3a16fxul1ErUrlcZvAHTOuUqybvMOpi8Lksn0peuZuXwjs5dv2H10ApDdoG6QTJr9kEy6NG9Ao7Tq13OiNzninHNx1ig9lcM6NuWwjk13TysqMpas28Ks5Rv5fvkGvl++kVkrNvDK2EVs2fHjhNIxK4OO2fXplJ1Bx+wMOmbVp02TtGp3dZcnDueci6OkJNGmSTptmqRzXNdmu6cXFRmL127ZnUzmFm5k7spNjJiylDWbd+xerk5yEjlN03cnlY7ZGXTIyqBd03Sy6yfmPIonDuecS4CkpB9uWjyhW/MfzVu9aXuQSAo3MWdl+Fy4kU9mrmDHrh9OL6TXSSYnM512TdNp3zSDnOLnzHRaNU4jOU6XDnvicM65aiYzow6ZGZnktc/80fSdu4pYtGYLC1ZtYsGqzeFjE3MKN/HJzEK27yzavWxqsmjbJJ0/n9Wdwzs1jd7EXvHE4ZxzNURKchIdsoKqqmhFRcay9VuZv2oTC1dtZsHqIKlkZlT+VVyeOJxzbh+QlCRaNU6jVeM0jugU523Fd/XOOef2NZ44nHPOVYgnDueccxXiicM551yFeOJwzjlXIZ44nHPOVYgnDueccxXiicM551yF1Ipm1SUVAgsiJmUBKxMUTqL5vtdOvu+1197sfzszy46eWCsSRzRJ+SW1MV8b+L77vtc2tXnfIT7771VVzjnnKsQTh3POuQqprYljSKIDSCDf99rJ9732qvT9r5XnOJxzzu252nrE4Zxzbg954nDOOVchtS5xSDpF0kxJsyXdkuh4qoqktpI+kTRd0lRJ1yc6pqomKVnSBEnvJjqWqiSpsaTXJc0I3//DEx1TVZF0Y/h5nyLpZUn1Eh1TvEh6WtIKSVMipmVKGilpVvjcpDK2VasSh6RkYDBwKpALXCApN7FRVZmdwE1m1g04DPh1Ldr3YtcD0xMdRAI8Aowws67AQdSS10BSa+A6IM/MDgSSgfMTG1VcPQucEjXtFmCUmXUBRoXje61WJQ6gDzDbzOaa2XZgKNA/wTFVCTNbambjw+ENBD8erRMbVdWR1AY4DXgy0bFUJUkNgaOBpwDMbLuZrU1sVFUqBUiTlAKkA0sSHE/cmNlnwOqoyf2B58Lh54ABlbGt2pY4WgOLIsYLqEU/nsUktQcOBsYkNpIq9XfgZqAo0YFUsY5AIfBMWE33pKSMRAdVFcxsMfAgsBBYCqwzs48SG1WVa25mSyH48wg0q4yV1rbEoRKm1arrkSXVB94AbjCz9YmOpypIOh1YYWbjEh1LAqQAvYDHzOxgYBOVVF1R3YX1+f2BDkArIEPSxYmNat9Q2xJHAdA2YrwN+/ChazRJqQRJ40UzG5boeKrQkcCZkuYTVE8eL+mFxIZUZQqAAjMrPrp8nSCR1AYnAvPMrNDMdgDDgCMSHFNVWy6pJUD4vKIyVlrbEsdYoIukDpLqEJwoG57gmKqEJBHUc083s4cSHU9VMrNbzayNmbUneM8/NrNa8c/TzJYBiyTtH046AZiWwJCq0kLgMEnp4ef/BGrJhQERhgMDw+GBwNuVsdKUylhJTWFmOyVdA3xIcIXF02Y2NcFhVZUjgUuAyZImhtNuM7P3ExiTqxrXAi+Gf5bmApcnOJ4qYWZjJL0OjCe4qnAC+3DzI5JeBo4FsiQVAHcC9wGvSrqSIJGeWynb8iZHnHPOVURtq6pyzjm3lzxxOOecqxBPHM455yrEE4dzzrkK8cThnHOuQjxxOOecqxBPHK5Wk9RU0sTwsUzS4ojxr+KwvcskFUp6Mhw/NrqZd0nPSjqnsrcdsf60cP+2S8qK13bcvqtW3QDoXDQzWwX0BJB0F7DRzB6M82ZfMbNr4rwNJKWY2c7o6Wa2BegZNsHiXIX5EYdzpZC0MXw+VtKnkl6V9L2k+yRdJOlbSZMldQqXy5b0hqSx4ePISojhhLBV28lhRz11w+nzi48WJOVJGh0O3yVpiKSPgOclHRDGOVHSJEld9jYm5/yIw7nYHAR0I+jvYC7wpJn1CXtSvBa4gaDDpIfN7AtJOQRN23SLYd1HRTQDA5ADvBv2VvcscIKZfS/peeCXBE3El6U30NfMtkj6J/CImRU3OZIc6w47VxpPHM7FZmxxvwaS5gDF/TpMBo4Lh08EcoP29ABoKKlB2HFWWT43s9OLRyQ9Gw7uT9C66/fh+HPAryk/cQwPq6MAvgZuDzuyGmZms8op61y5vKrKudhsixguihgv4oc/YEnA4WbWM3y0jiFplKWk/mOK7eSH7290P9qbigfM7CXgTGAL8KGk4/ciHucATxzOVaaPgN0nvSX13Mv1zQDaS+ocjl8CfBoOzyeokgL4aWkrkNQRmGtm/yBoYrvHXsbknCcO5yrRdUBeeBJ6GjBob1ZmZlsJmkB/TdJkgqObx8PZfwQekfQ5sKuM1fwMmBKeQ+kKPL83MTkH3qy6c1VK0mVAXlVcjhtDLPPDWFYmOhZXs/gRh3NVawtwavENgIlQfAMgkEpwFONchfgRh3POuQrxIw7nnHMV4onDOedchXjicM45VyGeOJxzzlXI/wPSWKmmf9R6zwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "value = [] #Create empty array for delta temp between ambient and body to show convergence as t->infinity\n",
    "x = t[1:] #use t values as x-axis to show t->infinity\n",
    "N = 30\n",
    "for i in range(1,N): \n",
    "    value.append(np.abs((Temp_num[i]-Temp_ambient)/Temp_ana[i])) #add delta temp to empty array\n",
    "plt.plot(x,value);\n",
    "plt.title('Numerical Solution Approaches Ambient Temperature');\n",
    "plt.xlabel('Time [Hours]')\n",
    "plt.ylabel('Delta Temeperature [degF]');\n",
    "print('Temperature approaches',Temp_ambient, 'degF as time approaches infinity')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_t = -np.log((Temp_dead - Temp_ambient)/(Temp_initial - Temp_ambient))/K\n",
    "#Rearranged Newton's cooling to find the time between known dead temperature and temp at 11 am\n",
    "time_death = 11-d_t #Subtract delta time value from initial time, 11 am\n",
    "print('Time of Death:',int(time_death),':',format(int(60*(time_death-int(time_death))),'02d'), 'am') #Nice format time"
   ]
  },
  {
   "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?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_time = np.array([[8,55],[9,58],[10,60],[11,65],[12,66],[13,67]]) #Matrix of time to ambient temp\n",
    "def get_Temp(time,am_or_pm):\n",
    "    ''' Determine the temperature of the ambient air given a time of day\n",
    "    Arguments\n",
    "    ---------\n",
    "    time - The time at which you would like to recieve the ambient temperature\n",
    "    am_or_pm - Either 'am' or 'pm' inputs to establish the time of day\n",
    "    Returns\n",
    "    -------\n",
    "    T_amb - the ambient temperature at the time of day \n",
    "    '''\n",
    "    if am_or_pm == 'am': #Establishes am or pm to avoid inputs above 12\n",
    "        time = time #am responce\n",
    "    else:\n",
    "        time = time+12 # pm responce\n",
    "    for i in range(len(temp_time[:,0])): #Iterate through temp-time matrix\n",
    "        if time == temp_time[[i-1],0]: #Triggered if time input matches value\n",
    "             T_amb = temp_time[[i-1],1] #Sets temp as the cooresponding ambient temp\n",
    "    return T_amb, time\n",
    "amb_temp,time = get_Temp(11,'am') # Test case with time and associated temp as output\n",
    "print('Ambient Temperature at',time,':',int(amb_temp), 'degF')\n",
    "\n",
    "plt.plot(temp_time[:,0],temp_time[:,1])\n",
    "plt.xlabel('Time of day')\n",
    "plt.ylabel('Temperature [degF]')\n",
    "plt.title('Time vs. Ambient Temperature');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This ambient temperature vs time graph looks correct, except for the jumps between established values. As each step is taken, the path between the two x-values is linear. This could be improved by creating a trend line to hit each point. This would then allow the time input for the 'get_Temp' function to be at any time in between the extablished values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Temp_ambient = np.array([55,58,60,65,66,67])\n",
    "Temp_initial = 85\n",
    "Temp_dead = 98.6\n",
    "N = 20\n",
    "t = np.linspace(0,10,N)\n",
    "\n",
    "Temp_ana = Temp_ambient[3] + (Temp_initial - Temp_ambient[3])*np.exp(K*t) #Analytical with constant ambient temp\n",
    "\n",
    "Temp_num = np.empty(len(t)) #Establishes empty array for numerical solution of temperatures\n",
    "delta_t = np.diff(t) #establish size of time steps from t\n",
    "for i in range(-1,len(t)-1):\n",
    "    Temp_num[0] = Temp_initial #Sets first array entry as the initial condition\n",
    "    if i < 0: #Uses ambient temp at 10 am for time between 10-11\n",
    "        Temp_num[i+1] = Temp_num[i] + (K*(Temp_num[i] - Temp_ambient[2])*delta_t[i])\n",
    "    if i < 1:#Uses ambient temp at 11 am for time between 11-12\n",
    "        Temp_num[i+1] = Temp_num[i] + (K*(Temp_num[i] - Temp_ambient[3])*delta_t[i])\n",
    "    if 2 > i <= 1: #Uses ambient temp at 12 pm for time between 12-1\n",
    "        Temp_num[i+1] = Temp_num[i] + (K*(Temp_num[i] - Temp_ambient[4])*delta_t[i])\n",
    "    if 3 > i <= 2: #Uses ambient temp at 1 pm for time between 1-2 \n",
    "        Temp_num[i+1] = Temp_num[i] + (K*(Temp_num[i] - Temp_ambient[5])*delta_t[i])\n",
    "    if i >= 3: #Uses ambient temp at 1 pm for time after 2, due to lack of future temperature information\n",
    "        Temp_num[i+1] = Temp_num[i] + (K*(Temp_num[i] - Temp_ambient[5])*delta_t[i])\n",
    "\n",
    "plt.plot(t,Temp_num,'ro',label='Non-Linear Numerical')\n",
    "plt.plot(t,Temp_ana,label='Analytical Temperature')\n",
    "plt.title('Temperature of Cooling Body; Non-Linear Ambient Temperature')\n",
    "plt.xlabel('Time [Hours]')\n",
    "plt.ylabel('Temeperature [degF]')\n",
    "plt.legend(loc='best');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Temp_10 = Temp_ambient[2]+(Temp_initial-Temp_ambient[2])*np.exp(K*-1) #Finds the temp at 10 am with cooresponding T_a \n",
    "d_t = -np.log((Temp_dead - Temp_ambient[1])/(Temp_10 - Temp_ambient[1]))/K #Finds time before 10 person had died\n",
    "time_death = 10-d_t #Subtracts value of time from 10 to get the time of death\n",
    "print('Time of Death:',int(time_death),':',format(int(60*(time_death-int(time_death))),'02d'), 'am') #Nice fomat time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nonlinear Euler approximation varies only in final values, not trend shape. This is a result of the same equations being carried out with only different ambient temperatures. As the graph displays, the body temperature comes to a rest at the ambient temperature at 1 pm, 67 $^oF$ . Additonally, the body's temperature comes to a rest faster than that of the linear analytical solution. This is again due to the higher ambient temperature."
   ]
  }
 ],
 "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
}