01_Getting-started-project.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "# 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 . \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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2 µs, sys: 3 µs, total: 5 µs\n",
      "Wall time: 8.11 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "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",
    "    your inputs...\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    your outputs...\n",
    "    \n",
    "    '''\n",
    "    differential=(Temp_t2-Temp_t1)/delta_t\n",
    "    K= -(1/(Temp_t1-Temp_ambient))*(differential)\n",
    "    '''differential is \"dT/dt\" '''\n",
    "    return K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.275"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "measure_K(85,74,65,2)\n"
   ]
  },
  {
   "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": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 5 µs, sys: 0 ns, total: 5 µs\n",
      "Wall time: 8.11 µs\n"
     ]
    }
   ],
   "source": [
    "def measure_K(Temp_t1,Temp_t2,Temp_ambient,t1,t2):\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",
    "    your inputs...\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    your outputs...\n",
    "    \n",
    "    '''\n",
    "    differential=(Temp_t2-Temp_t1)/(t2-t1)\n",
    "    K= -(1/(Temp_t1-Temp_ambient))*(differential)\n",
    "    '''differential is \"dT/dt\" '''\n",
    "    K\n",
    "    return K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.275"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "measure_K(85,74,65,11,13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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",
    "    your inputs...\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    your outputs...\n",
    "    \n",
    "    '''\n",
    "    K= -(1/(T-T_a))*(delta)"
   ]
  },
  {
   "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": 216,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "T_a= 65\n",
    "T_i= 85\n",
    "K=.275\n",
    "N=12\n",
    "T=[T_i]*(N+1)\n",
    "t=np.linspace(0,14,N)\n",
    "for i in range (0,len(t)):\n",
    "    slope= -K*(T[i]-T_a)\n",
    "    T[i+1]=T[i] + slope\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "def T_analytical(T_0,T_a,K,t):\n",
    "    T_t= T_a+(T_0-T_a)*np.exp(-K*t)\n",
    "    return T_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Temperature of Corpse vs. Time')"
      ]
     },
     "execution_count": 220,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3hUVfrA8e+bTkJIKAkdAkjvEBBp0lREESwgimtvP0Vsq6LY1l1ddRXLWhELrhULigIqSBELCEjvHUINLRAgkPL+/rg3MIRJMoFMJuX9PM88M3PvPXPfO2LeOefcc46oKsYYY0xOQYEOwBhjTPFkCcIYY4xXliCMMcZ4ZQnCGGOMV5YgjDHGeGUJwhhjjFeWIIwpRCLSQ0TWiUiqiPQNdDwllYj8Q0ReC3QcZZ0liDLA/WOV/cgSkSMe74cGOr4zISI7RKRroOPw8DTwvKqWV9UfvB0gIteJyF8ickhEtovI9yLSqYjjDBgRCffyb/Kwx/vLVfUJVR0W6FjLupBAB2D8T1XLZ78WkY3Azao6NXAR+UZEQlQ1o4Sdoy6wLI/zPQLcBdwGTAUygH7AAGB2QU5UFN+PP6jqUcDz3+QO4ApV/TVwURlvrAZhEJFgEXlMRNaLyG4R+VhEYt19TUQkQ0RuEpGtIrJHRG4UkXNEZKmI7BeRUR6fdbuITBORt0XkgIgsF5HuHvsriciH7i//LSLyhIgE5Sj7uojsA0a4558hIntFJFlExopItHv8F0A88JP7y3O4iPQVkbU5ru94LUNEnhWRT0TkcxE5CAzJ6/pz+b7udJuR9ojI1yJS1d2eBNTIjsdLucrA48CtqjpBVQ+r6jFV/UZVH3aPKede/3YRSRKR/4hIqLuvr4isdWPdCbzpse0f7ne0XkQGeZxzgIisFJGD7vc93GPfpSKy2P1vOEtEmuVyvR+IyL9ybPtRRO5wXz/mxntARFaISLfcvjtfuf+dxrivC/Rv0C1zm4iscr+TiSJS80xjKpNU1R5l6AFsBPrk2DYCmIXzxy0C+AB4393XBFDgFSAcuAQ4BHwFVAbqAPuAs93jb8f5VXwHEApcC+wFKrj7JwP/BSKB6sAC4LocZW8BgoFy7vl7AWFANZxf2c96xL4D6Orxvi+wNsf1HT8GeBY4ivOrPcg9R67X7+X76+d+Xiv32NHAlNziyVF2IHAEkDz++zzvxlIFqArMBUZ6XFsG8JT7fZTz2PZvd1sf4DBQzy2zB+jovq4MtHVfdwK2A+3d7/pWYDUQ4iWm8z2/U5ykfMSNsTWw3o1VgPrZ5y7Av8lTvjP3v9OY0/w3OARYATRy/w3+C5ge6P/3SuIj4AHYo4j/g3tPEBuALh7v67l/ZMTjf87KHvsPAQM83k8Ebndf3w5syPH5i4FBOM0vh4BQj303AJM9yq7OJ/4hwB8e708nQfzk6/V7Of/HwFMe72OBLKCat3hylL0J2JjP9W0Fenm8HwCs9Li2nN9fXyANiPDYNgF4wH290/2Oo3Oc533cxOOxbVP2H9kc24Pd68pONHcBk9zXzXESTU+8JBcf/036miB8/Tc4HRjqsS8USAeqFvX/byX9YU1MZZyICFAbmORW1ffj/KoPwvl1BpCpqns8ih3B+cPj+b68x/ukHKfZhPPrvC7Or+5kj3O9gvPrM9uWHPHVEJEv3KaFA8AYnF+uZ+L4OXy8fk813OsBQFX3AwcAX5ow9gBV3XOewt1ezfPz3deen71DVdNzFE1W1bQcZWq4rwcClwOb3ea7RHd7XeCR7Gt2rzvO23WoaiYwDrjK3XQ1TqJEVZfh1MCeBna5zXNVc35GISjIv8G6wFse15WMU8uq5Ye4SjVLEGWcOj+xsn+1xno8IlR192l+bM7/EesA23D+MKcCFT3OU0FV23mGlKPsf3B+LbZQ1QrAzTg1m9yOP4TTfAWA235fKccxx8ucxvVvw/kDlP35MUAF9zPy86sb+0Xedrqx7PD8fJzvzvOzvU2/XEVEInKU2eZ+5h+qejFOEv4J+NQ9ZgvweI5rjlTVr3OJ/VPgShE5C2gJfOMR91hV7YzTvBSB06QTSFuA63NcWzlVnR/guEocSxAG4C3gWRGpDSAi8SLS/ww+r7bb4RwiItfg/MH6SVU34PQhPC8i0SISJCINJe/bVKNxksoBEakD3Jdj/06cP0zZVgCVRKS3mxz+Qf7/zgty/Z8Ct4hIC/eP8nPANFXdkc85cBPOP4G3ReRit0M6VET6i8gzHp//hIhUFpF4YCTwUT4fHQo8JiJhItILOA/4SkSiRGSIiFTAaWI5CGS6ZUYDd4lIojjKi8glIhLp7QSq+gdOU9abwHeqesj9rpqJyLkiEo7zK/6IxzkC5S3gURFpDCAiFUXk8gDHVCJZgjDgdIxOBaaJc2fP70C7vIvk6RegLU7n9EjgUlVNcfddhdNuv9Ld/zknNzHl9DjQFUgBxuN0THp6GnjabU4Y5v4RvhunCSQJ5xd5fjUhn69fVb/H6RCegPMrvRrwt3w+37P80zjfyT/duDbjdBB/63G9y3FulV0I/ObGl5eNOE0oO4D3gBtUdb2770acJqcUnBsGrnPj+A0YDrwN7MfpoL4a7zWUbJ/idIJ/4rGtHPCiey3bcZp5Hgdw7zoq8l/tqvop8BrwtdssuRAnaZoCEqdWa0zhEJHbce5p7xPoWMoCcUZrv6aqZwU6FlP6WA3CGGOMV5YgjDHGeGVNTMYYY7yyGoQxxhivStVkfVWqVNGEhIRAh2GMMSXG/Pnzd6tqnLd9pSpBJCQkMG/evECHYYwxJYaIbMptnzUxGWOM8coShDHGGK8sQRhjjPGqVPVBGGNKh/T0dJKSkkhLS8v/YOOTiIgIatWqRWhoqM9lLEEYY4qdpKQkoqOjSUhIIJfZ0U0BqCp79uwhKSmJevXq+VzOr01MInKviCxzlwX8VEQiRORJd27/he6jXy5l+7pLBq4VkRF+C3LxOHipBTwZ6zwvHue3UxljfJOWlkblypUtORQSEaFy5coFrpH5rQbhrgE7HGimqkdEZBzOamAAL6nqC3mUDQZex5mBMQmYKyITVHV5oQa5eBx8NxzSjzjvU7Y47wFaDS7UUxljCsaSQ+E6ne/T353UIUA5EQnBWcRlm4/lOuIsG7leVY8Bn+EsvVi4fn7qRHLIln7E2W6MMWWc3xKEqm4FXsCZ7347kKKqP7m7h4nIYhF5T0Qqeilek5OXnkwilyUdReRWEZknIvOSk5MLFmRKzpUx89lujDFFZMKECTz77LOnVTYhIYHdu093QcgT/JYg3D/8A3AWgK8BRLmri70JNADa4CSOF70V97LN66yCqjpaVRNVNTEuzuto8dzF5LJEbW7bjTGmCGRkZHDJJZcwYoT/ul994c8mpj7ABlVNdhdZ/xrorKo7VTVTVbOAd3Cak3JKwllIPlstfG+e8l3vxyG03MnbQss5240xZdrGjRtp2rQpt9xyC82bN+f888/nyJEj9OjR4/iUPrt37yZ7/rcPPviAgQMH0r9/f+rVq8drr73GqFGjaNu2LZ06dWLv3r0ArFu3jr59+9K+fXu6devGypUrAbj++uu577776NmzJw899BAffPABw4YNA2Dnzp1ceumltG7dmtatW/P7778DMHDgQNq3b0/z5s0ZPXp0oX8H/rzNdTPQyV3j9gjQG5gnItVVdbt7zKXAUi9l5wINRaQezoLtQ3CWQyxc2R3RPz+FpmwBhU0t7yXBOqiNKTb+8d0ylm87UKif2axGBZ7o3zzf49asWcOnn37KO++8w+DBg/nqq5wr3p5s6dKlLFiwgLS0NM466yyee+45FixYwL333suHH37IPffcw6233spbb71Fw4YNmTNnDnfccQfTpk0DYPXq1UydOpXg4GA++OCD4587fPhwzj33XMaPH09mZiapqakAvPfee1SqVIkjR47QoUMHLr/8cipXrnz6X0wOfksQqjpHRL4E/sJZL3cBzkLpY0SkDU6T0UbgNgARqQGMUdV+qpohIsOAH4Fg4D1VXeaXQFsNhlaDObpvG/JKK9asXEjCJX45kzGmhKlXrx5t2rQBoH379mzcuDHP43v27El0dDTR0dHExMTQv39/AFq2bMnixYtJTU3l999/Z9CgQcfLHD169PjrQYMGERwcfMrnTps2jQ8//BCA4OBgYmJiAHj11VcZP348AFu2bGHNmjUlI0EAqOoTwBM5Nntd4F1VtwH9PN5PAib5L7qTRVSswaqal9A9aQLzlq4ksUWTojq1MSYPvvzS95fw8PDjr4ODgzly5AghISFkZWUBnDKuwPP4oKCg4++DgoLIyMggKyuL2NhYFi5c6PV8UVFRPsc2Y8YMpk6dyh9//EFkZCQ9evQo9JHnNheTh4T+DxEqGWz54aVAh2KMKaYSEhKYP38+AF9++WWBylaoUIF69erxxRdfAM4I50WLFuVbrnfv3rz55psAZGZmcuDAAVJSUqhYsSKRkZGsXLmS2bNnF/BK8mcJwkN4tcZsju9Nr4MT+HNVrlOkG2PKsL///e+8+eabdO7c+bRuJf3444959913ad26Nc2bN+fbb7/Nt8wrr7zC9OnTadmyJe3bt2fZsmX07duXjIwMWrVqxWOPPUanTp1O53LyVKrWpE5MTNQzXTDo6Ka5hL/fh48q3MI19+U62NsY40crVqygadOmgQ6j1PH2vYrIfFVN9Ha81SByCK/bgW0VO9A75Sv+XLM9/wLGGFNKWYLwosoFD1Bd9vLXxHcCHYoxxgSMJQgvwhqfz57yjei99zP+XH/mw9WNMaYksgThjQjRvf9Ow6CtzJr4v0BHY4wxAWEJIhdhrS7nYER1zk3+hLkb9wY6HGOMKXKWIHITHEJ4t7tJDFrN5InjAx2NMcYUOUsQeQjrcC1pobGcs+Mj5lktwhhzGjwn3cvrmG3bTsxHevPNN7N8ecHXR5sxYwYXX3xxgcvlxhJEXsKiCO50G+cF/8WXP0wJdDTGmNyU8KWDcyaIMWPG0KxZswBG5LAEkY/Qc24nPSiCxK0fMX+T1SKMKXaylw5O2QLoiaWDCyFJeJtOu3z58owcOZLWrVvTqVMndu7cCcB3333H2WefTdu2benTp8/x7dkOHjxIvXr1SE9PB+DAgQMkJCTwxRdfMG/ePIYOHUqbNm1OmVL8hx9+oF27drRu3ZrevXsD8Oeff9K5c2fatm1L586dWbVq1Rlfqzd+nayvVIisBO2uZeC8d7n/h99of1v/QEdkTNkyeQTsWJL7/qS5kHn05G3pR+DbYTB/rPcy1VrChfmv1uZtOu1Dhw7RqVMnnn76aR588EHeeecdHn30Ubp27crs2bMREcaMGcPzzz/Piy+eWA8tOjqaHj16MHHiRAYOHMhnn33G5ZdfzqBBg3j99dd54YUXSEw8eUBzcnIyt9xyC7/88gv16tU7vqZEkyZN+OWXXwgJCWHq1Kk88sgj+U5FfjosQfggtMswsua/S4stnzB/U1fa1/W2SqoxJiByJof8theAt+m0w8LCjrfzt2/fnilTnObnpKQkrrzySrZv386xY8eoV6/eKZ9388038/zzzzNw4EDef/993nkn78G4s2fPpnv37sc/q1KlSgCkpKRw3XXXsWbNGkTkeK2ksFmC8EXFumQ1u5Srl33P/VP+ov3NvQMdkTFlR36/9F9q4TYv5RBTG26YeNqnzW067dDQUEScVZGDg4PJyMgA4K677uK+++7jkksuYcaMGTz55JOnfGaXLl3YuHEjM2fOJDMzkxYtWuQZg6oeP5enxx57jJ49ezJ+/Hg2btxIjx49Tvs682J9ED4K6XYvUaTRYMNn/LV5X6DDMcZk89PSwQWdTjslJYWaNWsCMHZsLk1bwLXXXstVV13FDTfccHxbdHQ0Bw8ePOXYc845h5kzZ7JhwwaA401MnufyXHmusPk1QYjIvSKyTESWisinIhIhIv8RkZUislhExotIbC5lN4rIEhFZKCJnNkVrYajWgswGfbgp9Afe+MnbKqnGmIBoNRj6v+rUGBDnuf+rJ5YUPk0FnU77ySefZNCgQXTr1o0qVarketzQoUPZt28fV1111fFt119/PbfffvvxTupscXFxjB49mssuu4zWrVtz5ZVXAvDggw/y8MMP06VLFzIzM8/oOvPit+m+RaQm8CvQTFWPiMg4nBXitgHT3GVFnwNQ1Ye8lN8IJKqqz5MhFcZ033naMAvGXszI9Bu5/LbHaVfH+iKM8YfSPN33l19+ybfffsv//lf00/gUt+m+Q4ByIhICRALbVPUnVc1w988Gavk5hsKT0JXM6u24PXQSr05ZGehojDElzF133cWIESN47LHHAh2KT/yWIFR1K/ACsBnYDqSo6k85DrsRmJzbRwA/ich8Ebk1t/OIyK0iMk9E5iUnJxdG6LkTIbjbvdRmB1HrJrHA+iKMMQXw3//+l7Vr19KoUaNAh+ITvyUIEakIDADqATWAKBG5xmP/SCAD+DiXj+iiqu2AC4E7RaS7t4NUdbSqJqpqYlxcXKFeg1dNLiKrUgPuDPuOV6au9v/5jCmjStNql8XB6Xyf/mxi6gNsUNVkVU0HvgY6A4jIdcDFwFDNJWpV3eY+7wLGAx39GKvvgoIJ6jKcZmwgfe10Fm7ZH+iIjCl1IiIi2LNnjyWJQqKq7Nmzh4iIiAKV8+c4iM1AJxGJBI4AvYF5ItIXeAg4V1UPeysoIlFAkKoedF+fDzzlx1gLptUQsqY9wzAm8srUXrx/Q/HIXcaUFrVq1SIpKQm/NxuXIREREdSqVbAuX78lCFWdIyJfAn/hNCUtAEYDy4BwYIo7AGS2qt4uIjWAMaraD6gKjHf3hwCfqOoP/oq1wEIjCDrn/zhn6pP8a/WfLNrSiNa1vd6ta4w5DaGhoV5HIpui5bfbXAPB77e5ekpLQUc148djrRiX8BTvXd+haM5rjDGFKJC3uZZeETFI4o2cz2zWrlrC4iTrizDGlC75JghxtBSRC0Sku4hULorASoROdyDBIdwZ/gOvTF0T6GiMMaZQ5ZogRCRBRN4A1gEvAzcA9wG/iMhvIvI38TaLVFlSoTrS6kouD5rBwpVrWJKUEuiIjDGm0ORVg3ge+AI4S1V7q+oQVR2oqs2BK3A6kq8riiCLtS53E5x1jNsipvLKzzYuwhhTeuSaIFR1sKpOV9UsL/u2q+oLqvqBX6MrCao0RJpcxLUhU/h9xWaWbrVahDGmdMiriemfHq97FU04JVSXe4jIOMD1Eb/wsvVFGGNKibyamC7yeP2CvwMp0Wp3gLpd+L/wycxYsdVqEcaYUsFucy0sXe4h+uhOroyYwys/Wy3CGFPy5TWSOl5EhgPi8fo4VX3Vr5GVNA3Pg/hm3Jc6mcTl57B0awotasYEOipjjDltedUg3gfigCoerz0fxpMIdLmbyofXc1HEEl61WoQxpoTLtQahqiVjRYvipMXlMO1fPKw/0WV5a5ZtS6F5DatFGGNKprzuYionIjeKyE3ujKwmP8GhcM6d1DywgK4R660WYYwp0fJqYvoY2AMkA58UTTilQNu/QUQs/6g0hR+X7WT5tgOBjsgYY05LXgmiIrDcfVQqmnBKgfDy0PFWGuydSauIHVaLMMaUWHkliGuAW4HbgWuLJpxS4uzbIKQcz8TP4IdlO1ix3WoRxpiSJ6+pNraq6gOq+ndV3ViEMZV8UVWg7TU03z2ZBuEHrBZhjCmR8uqkHpLXbK3ubK+d/RNWKdB5GKKZPFfzVyYvtVqEMabkyauJqSawQERGi8htInKZiFwtIo+LyDScKcD35PXhInKviCwTkaUi8qmIRIhIJRGZIiJr3OeKuZTtKyKrRGStiIw4/UsMkIoJ0PxS2u/+hhrhR/nvNKtFGGNKlryamF4EEoHxQG2cuZk64ySFm9ypv1flVl5EagLDgURVbQEEA0OAEcDPqtoQ+Nl9n7NsMPA6cCHQDLhKRJqd1hUGUpe7kWOpPFd3HpOW7GDlDqtFGGNKjjznYlLVDFWdrKqPqupNqjpMVV9X1Q0+fn4IUE5EQoBIYBswABjr7h8LDPRSriOwVlXXq+ox4DO3XMlSvTU06EWXPeOoFJ5lfRHGmBLFb5P1qepWnFlgNwPbgRRV/Qmoqqrb3WO2A/FeitcEtni8T3K3nUJEbhWReSIyLzk5uTAvoXB0uYegQ8k822AZk5bsYNWOg4GOyBhjfOK3BOH2LQwA6gE1gCgRucbX4l62qbcDVXW0qiaqamJcXDGcIqped6jeht57P6NCeJDVIowxJYY/p/vuA2xQ1WRVTQe+xunD2Cki1QHc511eyibh9Htkq4XTPFXyiEDXewjev4F/Nt7IpKXbrRZhjCkR8k0QIhInIm+LyPfu+2Yicr0Pn70Z6CQike7tsr2BFcAETqxlfR3wrZeyc4GGIlJPRMJwOrcn+HDO4qnpJVCpPhcd+JyosGBetTuajDElgC81iA+AmZz4Rb8GuD+/Qqo6B/gS+AtY4p5rNPAscJ6IrAHOc98jIjVEZJJbNgMYBvyIk1TGqeoyn6+quAkKhs53EbJjAY82382kJdtZnLQ/0FEZY0yeRNVr0/6JA0TmqmoHEVmgqm3dbQtVtU2RRFgAiYmJOm/evECH4V36EXi5JenxLem69U6iI0L5/q6uRIQGBzoyY0wZJiLzVTXR2z5fahCHRKQSbiexiHQArBG9oELLQd0uhG6YxuxjV/DB/huZ/IktymeMKb7yWnI029+B74D6IjIT53bTK/waVWm0eBys/hEAQakVtJtK659hzdQYGva5McDBGWPMqfJNEKo6T0R6Ak1xbj9d7g5eMwXx81OQceSkTZFyjOjfnuFg178RHREaoMCMMcY7X+5iugwIV9VFQF/gIxEpdv0PxV5KktfN8Vm7+df3K4o4GGOMyZ8vfRBPqupBd+bW/sDnwFv+DasUiqnldfPB8Kp8Pm8LU5fvLOKAjDEmb74kiEz3+WLgDVX9Cgj3X0ilVO/HnY7qHKLOHU7T6hUY8fVi9qQeDUBgxhjjnS8JYruIvI4zWG2SO3DNnyOwS6dWg6H/qxBTGxCIrg7B4YSs/ZFRg1px4EgGI8cvJb/bjo0xpqj48od+MM5AuX6qug+ogpcpuo0PWg2Ge5fCk/vh/pVw4bOwYSZNt4/nvvMb8cOyHXyzcGugozTGGMCHBKGqqcBKoJeI/B9QRVUn+z2ysqDd9ZDQDX58lFtahZFYtyKPf7uMbfuP5FvUGGP8zZe7mEYCn+KMf6gFfCIiD/s7sDIhKAgu+S9oJsET7+XFQa3IzFIe+HIRWVnW1GSMCSxfmpiuATqo6khVHYmzmM+1/g2rDKlUD3o/AWunUDdpAo9e1Izf1u7hf7M3BToyY0wZ50uC2MTJA+pCgPX+CaeM6ngr1O4EP4zgqqah9Ggcx78nr2B9cmqgIzPGlGG+JIjDwDIRGSMi7+DMzLpfREaJyCj/hldGBAXBgNch4ygy8X6eu6wl4SHB3DduERmZWYGOzhhTRvmSICYCTwJ/ALOBp4BpwDL3YQpDlbOg50hYNZGqmyfyr4EtWLhlP2/NXBfoyIwxZVSeczGJSDDQXVWvy+s4U0jOuROWfwOTHqD/nX/yU+savDx1DT0ax9OiZkygozPGlDF51iBUNROoLiI2k1xRCAqGAW/AsVSY9Hf+OaA5laLCuG/cQtLSM/Mvb4wxhciXJqb1wCwReVhEhmc/8iskIo1FZKHH44CI3CMin3ts2ygiC3Mpv1FElrjHFdNVgPwgvgmc+xAs/4bYjZN5/opWrN6ZyqgpqwMdmTGmjPFlPYhkYAoQ6T58oqqrgDZwvKlqKzBeVV/OPkZEXgRS8viYnqq629dzlhpd7obl38LE++lx558MPbsO78xaT+8m8Zxdv3KgozPGlBG+rAfxGICIlHPfn84w397AOlU9fnO/iAjONB69TuPzSrfgUBj4BozuAT+M4JGL3uDXtbu5/4tF/HBPd8qH+5LXjTHmzPgykrqZiMwF1gBrRWSOiDQt4HmG4IzG9tQN2Kmqa3Ipo8BPIjJfRG7NI75bRWSeiMxLTk4uYFjFWLWW0O1+WPw5URun8uKg1mzbf4R/fb880JEZY8oIX/ogRgOPqGotVa0JjATe8fUE7uyvlwBf5Nh1FacmDU9dVLUdcCFwp4h093aQqo5W1URVTYyLi/M1rJKh298hvhl8fw+JVYO47dwGfDZ3Cz+vsLUjjDH+50uCiFbVKdlvVHUqEF2Ac1wI/KWqx/+qiUgIcBnO4kNeqeo293kXMB5nio+yJSTMGUCXugt+Gsk9fRrSpFo0D321hL2HbNVXY4x/+ZIgNrp3MNVyHyNwpt/wlbeaQh9gpap6XYdTRKJEJDr7NXA+sLQA5yw9araDLsNhwUeEb5zBS1e2IeXIMR79ZomtHWGM8StfEsSNQG1gkvuoBdzgy4eLSCRwHvB1jl2n9EmISA0RmeS+rQr8KiKLgD+Biar6gy/nLJXOHQFVGsF3d9O0knDfeY2ZtGQHExZtC3RkxphSTHL7FSoi4UB5Vd2TY3sV4KCqFrv1MRMTE3XevFI6ZGLLn/Du+ZB4I5n9XmTw23+wZudBfry3O9VjTl3K1BhjfCEi81U10du+vGoQr+D9FtR+gE3SV9Rqd4ROd8C8dwneNItRg1uTkaU8+OVia2oyxvhFXgmiu6rmvPMI4H9AD/+EY/LU61GoWA8m3EXdaHikX1NmrdnNR7Z2hDHGD/JKEOJtozo/V73uM34WFgkDXoN9G+HnfzL07Dp0bxTH05NWsGH3oUBHZ4wpZfJKELtFpH3OjSLSDtjrv5BMnhK6QodbYM5byJY5PH95K3ftiIW2doQxplDllSAeAL4SkUdF5EL38RjwlbvPBEqfJyGmNnx7J9UilX8ObMGCzft5+xdb6M8YU3hyTRCqOhvoBJQDbncf5YDOqvpH0YRnvAovD5e8CnvWwox/c0nrGlzcqjovTVnN0q15zX1ojDG+y289iB2qOlJVB7iPR1R1e1EFZ/LQoCe0uxZ+/y8kzeefA1pQKSqM+8ctsrUjjDGFwpeBcqa4Ov9fUL4afHsnFcOV565oxaqdB3nJ1o4wxhQCSxAlWUQM9JexyHEAACAASURBVH8FklfAL/+hZ+N4rj67DqNnrefPDXYfgTHmzOSaIETkA/d5WJFFYwqu0fnQ+iqYNQq2L2Jkv6bUrhjJ/V8sJPVoRqCjM8aUYHnVIDqKSE3gFhGJFpEKno+iCtD44IJnIKoKfHsnUSHKqMGtSdp3hKcnrgh0ZMaYEiyvBDEGmAE0AZbleJTNmVWLq8hKcNEo2LEEfn2ZxIRK3Na9AZ/+uZnpK3cFOjpjTAmV122uo1S1IfChqtZR1doejzpFGKPxRdOLofllMPM52Lmce89z1o548KvF7LO1I4wxpyHfTmpVvUVEWojI7e6jWVEEZk5Dv/9ARAX49g7CRRk1uA37Dx9jxNeLycqyCf2MMQXjy5rUdwLjgDru4wsRucPfgZnTEFXFSRLbFsAfr9GsRgUe6tuEH5ft5LFvl9qsr8aYAgnx4ZjbgI6qmgogIs8AvwNv+DMwc5qaXwZLv4bpz0DjftzUtSHJqUd5e+Z6yoeHMOLCJojYXIvGmPz5Mg5CgHSP9+n4MJuriDQWkYUejwMico+IPCkiWz2298ulfF8RWSUia91lTo0vRJwO69By8O2diGYxom8TrulUh7d/Wc9r09YGOkJjTAnhSw3if8BsEfnKfX8pMDa/Qqq6CmgDICLBwFZgPM5ypS+p6gu5lXWPfx1nudIkYK6ITFDV5T7Ea6KrwoXPwfjb4Pn6SFoK/4ypRZP61/HoFIgKD+HGrvUCHaUxppjLN0Go6vMiMh3ohlNzuF1V5xbwPL2Bdaq6ycfmjY7AWlVdDyAinwEDAEsQvpIg55G233mbsoWhh18krc49PPU9lA8PYXCH2gEO0hhTnPk01YaqznVve33xNJIDwBDgU4/3w0RksYi8JyIVvRxfE9ji8T7J3XYKEblVROaJyLzk5OTTCK2U+vkp0JPXh5D0I9x09H90bxTHQ18v5rtF2wIUnDGmJPD7XEwiEgZcAmQvX/om0ACn+Wk78KK3Yl62eb0FR1VHq2qiqibGxcUVQsSlREqS182SspW3r2lPh7qVuPfzhfy8YmcRB2aMKSmKYrK+C4G/VHUngKruVNVMVc0C3sFpTsopCfBs/6gF2M/dgoiplcv2mpQLC+bd6xNpVqMC//fxX/y+dnfRxmaMKRF8ShAiUktEerqvw0UkqgDnuAqP5iURqe6x71K8T9sxF2goIvXcGsgQYEIBzml6P+7cyZRT5UagSnREKGNv6EhC5Uhu/nAe8zftK/oYjTHFmi8D5W7E+eM8xt1UF/jWlw8XkUicO5G+9tj8vIgsEZHFQE/gXvfYGiIyCUBVM4BhwI/ACmCcqi7z6YqMo9Vg6P+qszQp4jw3OA/WT4MZ/wagYlQYH910NvHR4dzw/p8s22ar0RljTpD8RteKyEKcZqA5qtrW3bZYVVsVQXwFkpiYqPPmzQt0GMWXKkwYBgs+ggv+Dec4A+KT9h1m8Ft/cDQji89vO4ez4ssHOFBjTFERkfmqmuhtny9NTGmqeny2N3eMgg3FLYlEnFpF00vgx4edRAHUqhjJRzefjQhcM2YOW/YeDnCgxpjiwJcE8ZuIPAhEuP0QnwPf+zcs4zdBwXD5GGjQCybcBcudrp36ceX5301ncyQ9k6Fj5rDzQFqAAzXGBJovCeJB4CCwErgb+BkY6c+gjJ+FhMOVH0HNRPjqJlg3DYCm1Ssw9saO7Ek9yjVj5rDXpgk3pkzLM0G4zUnvqeqbqnqpqg50X2flVc6UAGFRMHQcVGkEnw2FLX8C0KZ2LGOu68DmvYe59r05HEhLz+eDjDGlVZ4JQlUzgeoiElpE8ZiiVK4iXPM1RFeDj6+AHc4dx+c0qMxb17Rn1Y6D3Pj+XA4fs7WtjSmLfGliWg/MEpGHRWR49sPfgZkiEl0V/vYNhEbB/y6FPesA6NkknleGtOWvzfu47X/zSUvPDHCgxpii5kuCSAamAJFAnMfDlBYV68K130BWBnw4EA44g9b7tazOc5e3Ytaa3dz16QLSM61l0ZiyJN9xECWJjYM4Q9sWwAf9oUINuGEyRFUGYOzvG3liwjIGtqnBqMFtCAqyu5yNKS3yGgeR73TfIjIFLxPlqer5hRCbKU5qtIWrP4OPLoePL4drJ0BEBa7rnEDq0Qz+8+MqIsNDeHpgC1uVzpgywJcFgx71eB0BXA4c9U84JuASusKgsfD5UPj0KrjmSwgtx509z+LQ0QzemLGOqLBgHunX1JKEMaWcLwsGzcmxaaaIzPRTPKY4aNwXBr4FX98CX1zvjJkIDuWBCxpz6GgG78zaQPnwUO7u0zDQkRpj/MiXJqYKHm+DgPZA9VwON6VFq0FwNAUm3g/f/B9cOhoJCuKJ/s05dCyTl6auJio8mJu71Q90pMYYP/GliWkZTh+EABnABuAWfwZliokON0NairM6XUQM9HuBoCDh2ctacvhYBv+auIKo8BCu6lgn0JEaY/zAlwRRX1VPGk4rIr6UM6VB1/vgyH74/VWIiIXejxESHMTLV7bl8LF5PDJ+CZFhwQxo43VFWGNMCebLOIicfRAAfxZ2IKaYEoHznoJ218GsF+C3VwEICwnirWva0zGhEveNW8SU5bZ0qTGlTa4JQkTiRaQ1UE5EWopIK/fRFWfQnCkrRODil6D5pTDlMZg/FoCI0GDevb4DLWrGcOcnf/GbLV1qTKmSV1PRRcCNOOtBv+Gx/SDwWH4fLCKNcaYGz1YfeByoCfQHjgHrgBtUdb+X8hvdc2UCGbkN5DBFJCgYLh0NRw/Cd3dDRAVofinlw0MYe0MHhoyezc1j5/HRzR1pX7dSoKM1xhQCX1aUG6yq487oJM6ssFuBs4HGwDRVzRCR5wBU9SEvZTYCiarq889SG0ldBI4dduZs2jofrvoMGvYBIPngUQa//Qe7U4/yv5vOpk3t2AAHaozxxRmtKKeq40TkAhG5T0QeyX4UMIbewDpV3aSqP7lrTgPMxqmhmJIiLBKu/hzim8Dn18Dm2QDERYfz0c1nE1MulMFv/8EnczZTmqZxMaYsyjdBiMgbwHXAfUA54BrgrAKeZwjwqZftNwKTcymjwE8iMl9Ebs0jvltFZJ6IzEtOTi5gWOa0lIuFa8ZDTE34eDBsXwxAzdhyTBjWlbPrVeKR8Uu4f9wimyrcmBLMlyamxaraSkQWqWprEYkGvvJ1LiYRCQO2Ac1VdafH9pFAInCZeglCRGqo6jYRiceZTfYuVf0lr3NZE1MR278F3usLmUfhhh+givO7ITNLeW3aWl7+eTUN48vzxtD2nBVfPsDBGmO8OaMmJiB7ceI0Eanmvk8owPkvBP7KkRyuAy4GhnpLDgCqus193gWMBzoW4JymKMTWdqYJV4UPB0BKEgDBQcLdfRry4Y0d2Z16jAGv/cqERdsCHKwxpqB8SRCTRCQWeAFYCGwEvizAOa7Co3lJRPoCDwGXqOphbwVEJMqtqSAiUcD5wNICnNMUlSoN4W9fw9EDMLoXjGoKT8bCSy3odmQ6E4d3pUn1Cgz/dAFPfLuUoxm28JAxJUV+a1IHAZNVdb+qfgHUA1qqqk+d1CISCZwHfO2x+TUgGpgiIgtF5C332BoiMsk9pirwq4gswhmUN1FVfyjIhZkiVL01nH0bHNrpLjakkLIFvhtO9U3f8dmtnbi5az3G/rGJwW/PJmmf198Fxphixpc+iNmq2qmI4jkj1gcRQC+1cJJCTjG14V6n8vfD0u088MVigoOFl65sQ8/G8UUcpDEmpzPtg5giIgMKOSZT2rj9D3lt79uiOt/d1ZXqMeW44f25vPDjKjKz7FZYY4orXxLEMGC8iBwRkb0isk9E9vo7MFPCxOQynCU00hlc50qoEsX4OzpzZWJtXpu+lr+9O4fkg7b+lDHFkS8JogoQCpQH4tz3cf4MypRAvR+H0HInbwsKgfRDMKY3JK8+vjkiNJjnrmjF81e0Yv6mfVz06iz+3GC/OYwpbnwZSZ0JDAIecl9XB9r4OzBTwrQaDP1fdfocEOd54JtwzVeQuhNG94DFX5xUZHBibb65s4uzpsQ7sxn9yzobfW1MMeJLJ/VrODWI7qraVEQqAT+qaoeiCLAgrJO6mDqwDb68ETb/Ae2vh77PnlTbOJiWzkNfLWbSkh2c36wq/xnUmphyoYGL15gy5Ew7qTur6m24A+ZUdS8QVojxmdKuQg247nvocg/M/wDGnAd71h3fHR0RyutXt+Pxi5sxbeUu+v/3V5ZuTQlcvMYYwLcEke6Oh1AAEakMZPk1KlP6BIfAef+Aq8fBgSR4+1xYemJ4jIhwY9d6fH7bOaRnZnHZm7/z6Z824Z8xgeRLgngd+AqIE5F/AL8Cz/k1KlN6NboAbpsF8U3hyxtg4t8h48RdTO3rVmTi8G6cXa8SD3+9hPu/sAn/jAmUfPsgAESkOdDHffuzqhbLaS+sD6IEyUyHqU/CH685I7EHjYVK9U7szlL+O20Nr/y8hkbx0bxxTTsaxNmEf8YUtjPtgwAIBtJxVoHztYwxuQsOhQuehiGfwL6NTpPT8gkndgcJ9/RpxIc3diQ59SiX/PdXvl9sE/4ZU5R8WQ9iJM5kezVwFvf5REQe9ndgpoxocpHT5FS5AYz7G0weARnHju/u1jCOicO70rhaNMM+WcCTE5ZxLMO6wIwpCr7c5roCaJ8986o7Ad98VW1aBPEViDUxlWAZx2DK4zDnTajZHq54HyrWPb47PTOLZyev5N1fN9CmdiyvD21HzdhyeXygMcYXZ9rEtAkI8XgfAqwvjMCMOS4kDC58FgZ/CLvXwNvdYOWk47tDg4N47OJmvDm0Het2pXLRq7OYsWpXAAM2pvTzJUEcBpaJyBgReQdYAuwXkVEiMsq/4Zkyp9kAuG0mVEyAz66CH0c6HdquC1tWZ8JdXalWIYIbPpjLiz+tsiYnY/zElyamm/Lar6rvFmpEZ8CamEqR9DT4aSTMHQO1OsKg90+aEDAtPZPHv13KuHlJJFSOZMSFTbigeTVEJIBBG1Py5NXE5NNtriWFJYhSaOlXMOFuZ6DdpaOh0clLoU9fuYtnJq1gza5UEutWZORFTWlbp2KAgjWm5DmjPggR6Ssic0Vkl033bYpci8vh1hlQoSZ8MsgZO5F5YuBczybxTL67G89c2pKNew5z6Ru/M+yTv9iy11atM+ZM+dLEtBYYjNP3cLyx153ZNa9yjYHPPTbVBx4HPnS3J+Csbz1YVfd5Kd8XeAVnDMYYVX02v4uxGkQpln4EfhjhzOVUpzM0Hwi//9dZkCimFvR+nNTGlzF65jpGz1pPVhZc17kuw3o2JCbSJv4zJjdn1MQkIjOAXqp62j2BIhIMbAXOBu4E9qrqsyIyAqioqg95OX41znrWScBc4CpVXZ7XeSxBlAGLx8E3d0LWsZO3h5ZzphtvNZgdKWmMmrKKL+YnUSEilOG9G/K3TnUJC7ExnsbkdKa3uT4IfCciD4jI8OxHAWPoDaxT1U3AAGCsu30sMNDL8R2Btaq6XlWPAZ+55UxZ12owRFY6dXv6Efj5KQCqxUTw/BWtmXhXN1rViuGf3y/nvJdmMmnJdpv8z5gC8CVB/APIBGJxVpLLfhTEEJzR2ABVVXU7gPvsbeX6msAWj/dJ7rZTiMitIjJPROYlJycXMCxTIqXu9L49x7rYzWpU4MMbO/LBDR2ICAnmjo//4vI3f2f+plNaNI0xXoTkfwjxqtr+dE8gImHAJUBBpufwdq+i159+qjoaGA1OE1OBAzQlT0wtSNly6vagYFg1GRr1Bfd2VxGhR+N4ujWM48v5W3jhp9Vc/ubvXNSyOg/2bUzdylFFHLwxJYcvNYifRaTXGZzjQuAvVc3+2bdTRKoDuM/ehsMmAbU93tcCbKY24/C2/nVwGJSrBJ8OgfcugE2/n7w7SLiyQx1m/L0H9/RpyLSVu+gzaib//H45+w/n6M8wxgC+JYhbgKkiknqat7lexYnmJYAJwHXu6+uAb72UmQs0FJF6bg1kiFvOGO/rXw94He5bDhe/DPs3w/sXwseDYMeSk4pGhYdwT59GzHygB5e3q8X7v22g+/PTGTNrPUcz8rwxz5gyx5e7mIK9bc/vNle3bCROX0J9VU1xt1UGxgF1gM3AIFXdKyI1cG5n7ece1w94Gec21/dU9en8zmd3MRkAjh2GP0fDr6Mg7QC0vAJ6jjxpvYlsK3cc4N+TVjJzdTK1K5Xjob5NuKhldRuRbcqMMx5JLSJDcP7IPyMitXA6mucXcpxnzBKEOcmRffDbqzD7TchKh/Y3QPcHILrqKYf+sjqZZyatYOWOg7StE8vIfk1JTPByt5QxpcyZjoN4DQgFuqtqUxGpBPyoqh0KP9QzYwnCeHVwB8x8Hv4a6/RVdLoDugyHiJiTDsvMUr76K4kXf1rFzgNHubBFNR7q24SEKtaRbUqvM00Qf6lqOxFZoKpt3W2LVLW1H2I9I5YgTJ72rIPpTzvzO5WrCF3vg463nNLhffhYBmNmbeCtmetIz8zimk51Gd6rIRWjwgIUuDH+c6YD5dJFJAj3NlO3D8HmVzYlT+UGcMV7cNsvzqJEUx6DV9vB/LEnze8UGRbC8N4NmfFAD65oX5uxv2+k+3+mM/qXdaSlW0e2KTtyrUGISIiqZojItcClQCLwHs68TP9Q1c+KLkzfWA3CFMiGWfDzPyBpLlRuCL0fg6aXHB9DkW31zoP8e9IKpq9KpmqFcK7sUIchHWpTw1a0M6XAaTUxZTctua+bA31wBrBNVdWl/gr2TFiCMAWmCqsmOdN0JK+EGm2hz5NQv8cph/62djfvzFrPzNXJCNCzcTxXn12HHo3jCQ6yu55MyXS6CeJ4n0NJYQnCnLasTFj8OUx/xhmlXb8H9H4CarY75dAtew/z+dwtfD5vC8kHj1I9JoIrO9RmSIc6VIuJKPLQjTkTp5sgkoBclxRV1WK33KglCHPGMo7CvPfgl//A4T3OEqi9HoNtC5xahsf04unNr+DnFTv5eM5mZq3ZTZBAryZVGXp2Hbo3irNahSkRTjdBbAfexPu8SKjqPwotwkJiCcIUmqMH4Y/XnTUnjh0CCQLPsaEe04sDbN5zmE/nbuaLeVvYnXqMmrHlGNKhNoM71KZqBatVmOLrjPsgSgpLEKbQHdoNr7SGY6mn7oupDfee3B13LCOLqSt28smczfy6djfBQUKfpvFcfXZdup1VhSCrVZhiJq8EkddsrvYv2ZioKk4NwpuULXBoD0RVPr4pLCSIfi2r069ldTbuPsSnczfz5bwkfly2k9qVyjGkQx0GJdYiPtpqFab4y6sGUUlVS9Ta01aDMH7xUgvv04sDBIVC477QZiic1QeCT13e9GhGJj8tc2oVf6zfQ0iQcH7zqlzdsS6dG1S2WoUJqDOei6mksARh/GLxOPhuuLNqXbbQctD9Qacje/HncCgZouKg5WBoczVUa+H1o9Ynp/Lpn5v5cn4S+w6nU7dy5PFaRZXy4UV0QcacYAnCmDO1eNwpdzFld1CTmQ5rp8LCj2HVD87EgNVaObWKloNOaoLKlpaeyY/LdvDJnM3M2bCX0GDh/ObVGNqxDuc0qGyzyZoiYwnCmKJyeC8s+dJJFtsXOk1QjS5wkkXD87w2Qa3d5dQqvvorif2H06lXJYqrOtbmiva1qWTzPxk/swRhTCDsXAYLP3FqH4d2QWQVp9bR5mqo1vKUw9PSM5m8dDufzNnM3I37CAsOolvDKvRsEk/PJvHUtKk9jB9YgjAmkDLTYe3PbhPUZLcJqqVHE1SVU4qs3nmQz/7cwpQVO9iy1+n7aFItmh6N4+nVJJ52dWIJCfZlrk1j8hawBCEiscAYoAXObLA3AvcAjd1DYoH9qtrGS9mNwEEgE8jI7QI8WYIwxd7hvc504ws/dkZnB4VAo77Q+ipoeD6EuE1Kbp+HpiSRUb4Gs+rcwTv7E5m7cS8ZWUqFiBC6N4qjZ+N4ejSOo7J1cJvTFMgEMRaYpapj3LWlI1V1v8f+F4EUVX3KS9mNQKKq7vb1fJYgTImyczks+gQWfe42QVV27oKKqgKzXjj1rqn+r3Kg0aX8tmY301buYvqqZHanHkUEWteKpadbu2heo4LdOmt8FpAEISIVgEU4S5WechJxbtPYDPRS1TVe9m/EEoQpCzIzYJ1HE1TmMe/H5Ri5nZWlLNt2wE0Wu1iUtB9ViIsOp0ejOHo1iadrwypER5zaMW5MtkAliDbAaGA50BqYD9ytqofc/d2BUbkGJrIB2IfTNPW2qo7O5bhbgVsB6tSp037Tpk2FfSnGFJ3De+H5ernvf3ADRHpfK3t36lFmrkpm+qpd/LI6mQNpGYQECR0SKtHL7ehuEBdlt9CakwQqQSQCs4EuqjpHRF4BDqjqY+7+N4G1qvpiLuVrqOo2EYkHpgB3qeoveZ3TahCmVMhr5DYCNdpAvXOdKcnrdDplyVSAjMws5m/ax/RVyUxfuYtVOw8CULtSOXo1jqdHk3jOqV+ZiNBgv12GKRkClSCqAbNVNcF93w0YoaoXiUgIsBVor6pJPnzWk0Cqqr6Q13GWIEypkNvI7c7DnVll18+EpD8hKwOCw6HO2U6yqNfDSR5Bp/7R37r/CNNX7mL6yl38tm43aelZRIQG0bmBcxttL7uNtswKZCf1LOBmVV3l/pGPUtUHRKQv8LCqnptLuSggSFUPuq+nAE+p6g95nc8ShCk18hq5DXA0FTb9DhtmwvoZsNPtm4iIgYRubsI4F6o0PGUJ1bT0TGav38OMVclMW7mLzXsPA9AwvjyJCRVpXSuW1rVjaRhf3m6lLQMCmSDa4NzmGgasB25Q1X0i8gFO7eItj2NrAGNUtZ+I1AfGu7tCgE9U9en8zmcJwpRZqclOstgwE9bNgJTNzvboGk6yqH+ukzAqVD+pmKqyLvkQM1bt4pc1u1m4eR8H0jIAKBcaTMuaMbSuHUOb2hVpXTuGmrHlrA+jlLGBcsaUJaqwb4PTFLV+Bmz4BY64EzPHNTnRf5HQxalxeNRWNKYWyR0f4vfIXizcsp9FSftZtu0AxzKyAKhSPux4DaN17Vha14ohNtKmAynJLEEYU5ZlZcHOJU6yWD/TaZrKOAISDLF1nA7xrIwTx+dYLe9YRhardhxkYdJ+Fm1xHmuTU8n+05FQOdJNFk7SaF6jgnV+lyCWIIwxJ2QchaS5TsL47RXv4y4iYuHqz6FqcwiPPmX3gbR0lialeCSNFHYcSAMgJEhoWr0CrWvH0LpWLG1qx1I/rryt0V1MWYIwxnj3ZCzOUKM8VKznzB3l+ahQ85TO7x0paSzKThhJ+1m8JYWDR52aSfnwELc/I5Y2tZ3nahUirD+jGDjdJUeNMaVdTC3vYy6iq8PFL8OOJU7z1I4lsGLCif0RsackjWpVGlOteTUuaF4NcEZ6r9996HjCWLhlP+/+up70TCchVY4Ko0FceerHRVE/Lsp9XZ7aFcvZ3VPFhNUgjCnLchtz4dEHcdzRg878UdkJY8cS532GWzYo1OkEr9bSWVGvWkuo2uKkkd9p6Zns+PVDqsx5jqi0HewOiuNlHcLHRzqdOH2wUKdS5PGE4SSPKOpXKU9FWx+j0FkTkzEmd/mNuchLVibsWXdy0tixFFJ3nDimQq0TSePoQZj/AWSkndgfWo5DF4xiZdyFrE9OZf3uQ6zb5Txv2nPoeI0DoFJUGPWrRLm1jvLHayB1KkUSarWO02IJwhhTtFKTT00au1eDZno/vlxFGDQWKiY4/RvBTut3RmYWSfuOsC45lfXJh1i/O5V1yYdYn5zK7tQTneshQU6t40RTlZNA6leJolJUmPV15MEShDEm8NLT4Olq5NspHhTi1GQqJkBsXee5ovscm+A0WYmQciTdqXFkJ45dzvPG3Yc5lpl1/ONiI0OpXyWK2pUiqRYTQY2Ycic9V44KK9PTo1sntTEm8EIj8u4Uv/Rt2L8J9m10H5tg5UQ4nGPG/7BoqFiXmIoJtI2tS9uKCVAvAdrVhdhmZAZHsNWtdaxLTiV69Xh6bXuLyruS2a6VeS59MBOyup74uOAgqsaEU71COarHRlgS8WA1CGNM0SlIp3i2o6lu4vBIHscTyaYTneTZylc7Ues4dgjW/HTSWA8NKceWrs+yIq4vO1LS2JZyhB0paWzfn8b2A85rz34PKN1JxJqYjDHFx5l0iuekCqm7Tq15ZCeR3KZNDwqGmh2gfLz7qOo8R8WTFRXPvqCKbM+IZmtqVoGSyNDI2dyR+QlVspJJCYtndr1h7Kk/gMpRYVSMDKNyeec5NjKs2AwctARhjCmb8hoImNDNSS6pOyFtv/djImKc5BHlmUjiyIqK52BIZXZlRbM1owKb06KI2fAdF67/N2F69HjxwxrG/7d377FV1nccx98fqAXalYKKTrkMZUwlDtEQL1O3OTQyZ7wkm3Fxm9NtzmReNzJRE2OWxZDppi4ajGMOkxE34yVz6pzGOdxFEUURFK/otFgoDKgCpS30uz9+vzOenj7n9EB7+pzD+b6S5lz6XD5tes63v+d5zu87t/v7vQ5pQfiM4diGesY27MN+jSPYt7GesY31oZDE233zvso1fYmfg3DO1aZC5zyaJ8J3H931eEcnbF0fisWW3G1b6BWeu9+6PNx2fcIwoDl+TQVA4V3fenrtpkFd3Na0iGuOP5R2a2TjzlFs2DmKtu5RrO0cQWtHHRu3dfPu+i1sfL+LTdu66ClQzxrqh/cpGid3PMOsj+6iqXMdGuhoLIUXCOfc3mvWDennPGbd0Hu5uhGhmDRP6H+bXdti4WjbNQLZ0gaL56UuPqyznfGL5zA+7ZsaHkYpI5vhwDHYyDF079NEx/Amtg1r4mMaaaeBTTsb2LBjFOu6R9LaOYKPPh7BmJbnmb3jLkYpnl9p/zD8rDBoRcILhHNu75V7oxyscx4A9Q1QPzmcCE96ZVH6aGX0eLjocdjeDh2bw+GsAve1vZ36dzJEaAAACExJREFU9hbqt2+muWMzB/V0F8+SfxqjuyP8rF4gnHOuBNPPG9TDLgUVGq2cemPfYlIKs7Ct7ZtjAWnvff+Ja9LXa++3i3PJylogJI0hdJQ7knCm6GLgdOAHwPq42HVm9njKurOB24HhhE5z6eM355yrBIM9WpHiaKUBRh/c9/vP3VHg/EoJh8lKVO4RxO3AE2b2dUn1QAOhQNxqZrcUWknScOBO4DSgBVgq6REze73MeZ1zbs8N1WgFSj+/MgBlm91K0mjgi8BvAcysy8wKXEvWx7HAO2a22sy6gD8AZ5cnqXPOVaHp54UPGDZPBBRui33gcA+UcwRxKOEw0u8kHQW8BFwZv3eZpO8ALwI/MbNNeeuOB5JjpxbguLSdSLoEuARg0qRJg5feOecqXZlHLOWcH7cOOAaYb2ZHA1uBucB8YAowA2gFfpmybtpHDFOvDjazu81sppnNHDdu3KAEd845V94C0QK0mNmS+PgB4BgzW2dmO82sB/gN4XBS2roTE48nAB+VMatzzrk8ZSsQZrYW+FDSYfGpWcDrkg5KLHYusDJl9aXAVEmHxJPb5wOPpCznnHOuTMp9FdPlwKL4Jr8auAj4taQZhENG7wM/BJB0MOFy1jPMbIeky4C/Ei5zvcfMXitzVueccwk+WZ9zztWwmpnNVdJ64D97uPr+wIZ+l6oM1ZQVqitvNWWF6spbTVmhuvIOJOtnzCz1Cp+9qkAMhKQXC1XRSlNNWaG68lZTVqiuvNWUFaorb7mylvMqJuecc1XMC4RzzrlUXiB2uTvrALuhmrJCdeWtpqxQXXmrKStUV96yZPVzEM4551L5CMI551wqLxDOOedS1XyBkDRb0puS3pE0N+s8xUiaKOkZSaskvSbpyv7Xypak4ZJelvRo/0tnS9IYSQ9IeiP+jk/IOlMhkq6OfwMrJd0naWTWmZIk3SOpTdLKxHP7SnpK0tvxdmyWGXMKZL05/h28Kunh2PysIqTlTXxvjiSTtP9g7KumC0SiMdFXgWnANyVNyzZVUTsI06MfARwP/KjC80KY4n1V1iFKlGtwdThwFBWaW9J44ApgppkdSZiO5vxsU/WxEJid99xc4Gkzmwo8HR9XgoX0zfoUcKSZTQfeAq4d6lBFLKRvXiRNJDRZ+2CwdlTTBYIqa0xkZq1mtize/4TwBjY+21SFSZoAfI3QdraiDbDBVRbqgFGS6gidGitqtmMzexbYmPf02cC98f69wDlDGqqAtKxm9qSZ7YgPnyfMKF0RCvxuAW4FfkqB1gh7otYLRFpjoop9w02SNBk4GlhSfMlM3Ub4g+3JOkgJkg2uXpa0QFJj1qHSmNka4BbCf4qtQLuZPZltqpIcaGatEP7ZAQ7IOE+pLgb+knWIYiSdBawxs+WDud1aLxAlNyaqJJI+BTwIXGVmH2edJ42kM4E2M3sp6ywlKtTgquLEY/dnA4cABwONkr6Vbaq9k6TrCYd2F2WdpRBJDcD1wOA1o45qvUBUXWMiSfsQisMiM3so6zxFnAicJel9wqG7r0j6fbaRikptcJVhnmJOBd4zs/Vm1g08BHwh40ylWJfrBxNv2zLOU5SkC4EzgQussj8wNoXwz8Ly+HqbACyT9OmBbrjWC0RVNSaSJMIx8lVm9qus8xRjZtea2QQzm0z4vf7NzCr2v9xCDa4yjFTMB8Dxkhri38QsKvSEep5HgAvj/QuBP2WYpShJs4FrgLPMbFvWeYoxsxVmdoCZTY6vtxZC9861A912TReIeBIq15hoFXB/hTcmOhH4NuG/8Vfi1xlZh9qL5BpcvUromX5TxnlSxVHOA8AyYAXhdVxR00JIug94DjhMUouk7wHzgNMkvU242mZelhlzCmS9A2gCnoqvs7syDZlQIG959lXZIyfnnHNZqekRhHPOucK8QDjnnEvlBcI551wqLxDOOedSeYFwzjmXyguEqwmS9ktcGrxW0prE43+XaZ9HS1oQ798oaU459lNg35+XtHCo9uf2TnVZB3BuKJjZfwmfbUDSjcAWM7ulzLu9Dvh5OXcgqS4xqdz/mdkKSRMkTTKzQZvd09UWH0G4midpS7z9sqTFku6X9JakeZIukPSCpBWSpsTlxkl6UNLS+HViyjabgOl5k6dNk/R3SaslXZFY9sexr8NKSVfF5ybn9SeYEwsbcRs3SVoMXCnpG3Hd5ZKeTezvz1TeNOCuivgIwrnejgKOIEynvBpYYGbHKjRnuhy4itA34lYz+6ekSYRP4h+Rt52ZQH5Dl8OBUwif0H1T0nxgOnARcBxh8sgl8Y1/Uz85x5jZlwAkrQBON7M1eY1tXiRMOPiLkn965xK8QDjX29LclNSS3gVy02ivILy5Q5gsb1qYBgmA0ZKaYo+OnIMI04cnPWZmnUCnpDbgQOAk4GEz2xr3+RBwMv3PCfbHxP1/AQsl3U+YuC+njTDbq3N7xAuEc711Ju73JB73sOv1Mgw4wcw6imynA8hvA5rc9s64vbQp5yFMMZ08BJy/ra25O2Z2qaTjCM2ZXpE0I55zGRlzOLdH/ByEc7vvScIkjwBImpGyzCrgsyVs61ngnDgzayNwLvAPYB1wQLz6agRh2ulUkqaY2RIzuwHYwK4p7D9H38NczpXMRxDO7b4rgDvjrK91hDf5S5MLmNkbkppTDj2Rt9yyeDnqC/GpBWb2MoCknxE6Br4HvFEkz82SphJGI08DuRPjpwCP7e4P51yOz+bqXJlIuhr4xMyGvCd3HHUsBk5KuwzWuVL4ISbnymc+vc87DKVJwFwvDm4gfAThnHMulY8gnHPOpfIC4ZxzLpUXCOecc6m8QDjnnEvlBcI551yq/wHi1LXpbI6kagAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,T[:-1],'-',label='numerical')\n",
    "#t=np.linspace(0,14)\n",
    "plt.plot(t,T_analytical(85,65,.275,t),'o-',label='analytical')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('Time (hours)')\n",
    "plt.ylabel('Temperature of Corpse (°F)')\n",
    "plt.title('Temperature of Corpse vs. Time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_a= 65\n",
    "T_i= 85\n",
    "K=.275\n",
    "T_f = 120\n",
    "T=[T_i]*(N+1)\n",
    "t=np.linspace(0,T_f, T_f)\n",
    "for i in range (0,len(t)):\n",
    "    slope= -K*(T[i]-T_a)\n",
    "    T[i+1]=T[i] + slope\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Temperature of Corpse vs. Time')"
      ]
     },
     "execution_count": 221,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXwV1fnH8c83C9khLGEn7OIGokRQQSqgVq11qUtxq1Xr0rp0r1vt8vvV/qxVW1utitRqrVL31lasuG8VZBFEBATCvgaQLRAg8Pz+mAlewk1yA7mZLM/79ZpX7p2ZM/NMCPe555yZc2RmOOecc5WlRB2Ac865hskThHPOubg8QTjnnIvLE4Rzzrm4PEE455yLyxOEc865uDxBOFeHJJ0gaYGkLZJOiTqexkrSLyXdF3UczZ0niGYg/LCqWHZL2hbz/qKo4zsQklZJGhZ1HDFuB+40s1wz+0+8HSRdKmmapFJJKyX9W9Ix9RxnZCRlxPmb3Brz/hwz+7mZXRd1rM1dWtQBuOQzs9yK15IWAd8ys9eiiygxktLMrLyRnaM7MKua890CXA9cDbwGlAOnAWcCE2tzovr4/SSDmW0HYv8mVwHnmtl70UXl4vEahENSqqTbJBVLWivpCUn54baDJZVLukLScknrJF0u6VhJn0jaIOmemGNdI+kNSQ9J2iTpU0nDY7a3kfTX8Jv/Ukk/l5RSqez9kj4HbgrP/5ak9ZJKJD0mKS/c/xmgPTAh/OZ5g6RTJM2vdH17ahmS7pD0pKSnJG0GRld3/VX8vq4Nm5HWSXpeUodw/TKgc0U8ccq1BX4GXGVmL5rZVjPbYWb/MLObw32ywutfKWmZpN9KSg+3nSJpfhjrauCBmHW/DH9HxZLOiznnmZLmSNoc/r5viNl2tqSPw3/DdyUdWsX1PirpV5XWvSLpO+Hr28J4N0maLen4qn53iQr/ncaGr2v1NxiWuVrS3PB38pKkLgcaU7NkZr40owVYBJxYad1NwLsEH26ZwKPAX8JtBwMG3AtkAGcApcBzQFugEPgcGBLufw3Bt+LvAOnAN4D1QMtw+8vAH4FsoBPwEXBppbJXAqlAVnj+kUALoCPBt+w7YmJfBQyLeX8KML/S9e3ZB7gD2E7wrT0lPEeV1x/n93daeLwB4b5jgFeriqdS2bOAbYCq+fe5M4ylHdABmAzcGnNt5cD/hL+PrJh1/xeuOxHYCvQMy6wDBoev2wJHhq+PAVYCg8Lf9VXAZ0BanJhOjv2dEiTlbWGMRwDFYawCelWcuxZ/k/v8zsJ/p7H7+Tc4GpgNHBT+Df4KeDPq/3uNcYk8AF/q+R88foJYCAyNed8z/JBRzH/OtjHbS4EzY96/BFwTvr4GWFjp+B8D5xE0v5QC6THbLgNejin7WQ3xjwY+iHm/PwliQqLXH+f8TwD/E/M+H9gNdIwXT6WyVwCLari+5cDImPdnAnNirq3y7+8UoAzIjFn3IvDj8PXq8HecV+k8fyFMPDHrFld8yFZanxpeV0WiuR4YH74+jCDRjCBOcknwbzLRBJHo3+CbwEUx29KBnUCH+v7/1tgXb2Jq5iQJ6AaMD6vqGwi+1acQfDsD2GVm62KKbSP44Il9nxvzflml0ywm+HbeneBbd0nMue4l+PZZYWml+DpLeiZsWtgEjCX45nog9pwjweuP1Tm8HgDMbAOwCUikCWMd0CE85z7C9R1jjx++jj32KjPbWaloiZmVVSrTOXx9FnAOsCRsvisK13cHbqm45vC6C+Jdh5ntAp4GLghXXUiQKDGzWQQ1sNuBNWHzXIfKx6gDtfkb7A48GHNdJQS1rK5JiKtJ8wTRzFnwFaviW2t+zJJpZmv387CV/yMWAisIPpi3AK1jztPSzI6KDalS2d8SfFs83MxaAt8iqNlUtX8pQfMVAGH7fZtK++wpsx/Xv4LgA6ji+K2AluExavJeGPtX4m0MY1kVe3yC313sseMNv9xOUmalMivCY35gZqcTJOEJwLhwn6XAzypdc7aZPV9F7OOAr0vqA/QH/hET92NmdhxB81ImQZNOlJYC36x0bVlmNjXiuBodTxAO4EHgDkndACS1l/TVAzhet7DDOU3SxQQfWBPMbCFBH8KdkvIkpUjqq+pvU80jSCqbJBUCP6i0fTXBB1OF2UAbSaPC5PBLav47r831jwOulHR4+KH8G+ANM1tVwzkIE87/Ag9JOj3skE6X9FVJv445/s8ltZXUHrgV+FsNh04HbpPUQtJI4CTgOUk5kkZLaknQxLIZ2BWWGQNcL6lIgVxJZ0jKjncCM/uAoCnrAeBfZlYa/q4OlfQlSRkE3+K3xZwjKg8CP5XUD0BSa0nnRBxTo+QJwkHQMfoa8IaCO3v+CxxVfZFqvQMcSdA5fStwtpltDLddQNBuPyfc/hR7NzFV9jNgGLAReIGgYzLW7cDtYXPCdeGH8HcJmkCWEXwjr6kmlPD1m9m/CTqEXyT4lt4RuKSG48eWv53gd/K/YVxLCDqI/xlzvZ8S3Co7HXg/jK86iwiaUFYBjwCXmVlxuO1ygianjQQ3DFwaxvE+cAPwELCBoIP6QuLXUCqMI+gEfzJmXRZwd3gtKwmaeX4GEN51VO/f2s1sHHAf8HzYLDmdIGm6WlJQq3Wubki6huCe9hOjjqU5UPC09n1m1ifqWFzT4zUI55xzcXmCcM45F5c3MTnnnIvLaxDOOefialKD9bVr18569OgRdRjOOddoTJ06da2ZFcTb1qQSRI8ePZgyZUrUYTjnXKMhaXFV27yJyTnnXFyeIJxzzsXlCcI551xcniCcc87F5QnCOedcXElNEJK+L2lWOC3gOEmZkn4Rju0/PVxOq6LsKeGUgfMl3ZTMOJ1zzu0raQkinAP2BqDIzA4nmJVqdLj5d2Y2MFzGxymbCtwPnAocClxQ1Xy5zjnnkiPZTUxpQJakNIJJXFYkWG4wwbSRxWa2A/g7wdSLda5s5y7GvLOAScXrat7ZOeeakaQlCDNbDtxFMN79SmCjmU0IN18n6WNJj0hqHad4F/aeenIZVUzpKOkqSVMkTSkpKal1nBI88t4i7n71s1qXdc65piyZTUytCb719ySYHzcnnF3sAaA3MJAgcdwdr3icdXFHFTSzMWZWZGZFBQVxnxavVkZaKtd8qRcfLlzPBwu8FuGccxWS2cR0IrDQzErCSdafB44zs9VmtsvMdgMPEzQnVbaMYCL5Cl1JvHmq1kYPLqR9XgZ/eH1esk7hnHONTjITxBLgGEnZkgSMAmZL6hSzz9nAJ3HKTgb6SuopqQVB5/aLyQo0Mz2Vq7/Umw+K1/HhwvXJOo1zzjUqyeyDmAQ8C0wDZobnGkMwYf1MSR8DI4DvA0jqLGl8WLYcuA54hWAS+qfNbFayYgW4cHAh7XIzuPd174twzjloYhMGFRUV2YGM5vrwO8XcPn42z15zLEU92tRhZM451zBJmmpmRfG2+ZPUMS46ppC2OS241/sinHPOE0Ss7BZpXDm8F+/OW8u0JZ9HHY5zzkXKE0QllxzTndbZ6X5Hk3Ou2fMEUUlORhrfOr4Xb80tYcbSDVGH45xzkfEEEcelx/Ug32sRzrlmzhNEHLkZaVwxtCevz1nDJ8s3Rh2Oc85FwhNEFS4d2oOWmWl+R5NzrtnyBFGFlpnpXD6sJ69+uppZK7wW4ZxrfjxBVOOyoT3Jy0jjj6/PjzoU55yrd54gqtEqK53LhvbgP7NWMWfVpqjDcc65euUJogaXD+tJrtcinHPNkCeIGuRnt+DS47oz/pOVfLZ6c9ThOOdcvfEEkYBvDetFVnoqf3zDaxHOuebDE0QCWue04BvH9uDfH69g/potUYfjnHP1whNEgq48vieZaanc94Y/F+Gcax48QSSobW4GlxzbnRdnrKC4xGsRzrmmzxNELVx5fC9apKVw35veF+Gca/qSmiAkfV/SLEmfSBonKVPSbyXNkfSxpBck5VdRdlE4Nel0Sfs/TVwdKsjL4KIh3fnn9BUsWlsadTjOOZdUSUsQkroANwBFZnY4kAqMBl4FDjezAcBnwM3VHGaEmQ2sajq8KFw9vBdpKeJ+r0U455q4ZDcxpQFZktKAbGCFmU0ws/Jw+0Sga5JjqFPtW2ZyweBCnv9oOUvWbY06HOecS5qkJQgzWw7cBSwBVgIbzWxCpd0uB16u6hDABElTJV1V1XkkXSVpiqQpJSUldRF6jb59Qm9SU8Sf3vJahHOu6UpmE1Nr4EygJ9AZyJF0ccz2W4Fy4IkqDjHUzI4CTgWulTQ83k5mNsbMisysqKCgoE6voSodWmYy+uhuPDt1GUvXey3COdc0JbOJ6URgoZmVmNlO4HngOABJlwKnAxeZmcUrbGYrwp9rgBeAwUmMtda+fUJvUiQeeHtB1KE451xSJDNBLAGOkZQtScAoYLakU4AbgTPMLO7Xb0k5kvIqXgMnA58kMdZa69Qqi/OKuvLMlKUs37At6nCcc67OJbMPYhLwLDANmBmeawxwH5AHvBrewvoggKTOksaHxTsA70maAXwIvGRm/0lWrPvrOyP6APDgW16LcM41PWnJPLiZ/Rz4eaXVfarYdwVwWvi6GDgimbHVhS75WZw7qCtPTV7Kd0b0plOrrKhDcs65OlNjDUKB/pK+LGm4pLb1EVhj8Z0T+rDbjIfeLo46FOecq1NVJghJPST9CVgA/B64DPgB8I6k9yVdEvYtNGvd2mTztaO68OSHS1i9qSzqcJxzrs5UV4O4E3gG6GNmo8xstJmdZWaHAecS9BNcWh9BNnTXjujDrt1ei3DONS1VJggzO9/M3jSz3XG2rTSzu8zs0aRG10h0b5vDWQO78MSkxazZ7LUI51zTUF0T0//GvB5ZP+E0XteN7MPOXbt5+B2vRTjnmobqmpi+EvP6rmQH0tj1bJfDmQO78PjExazdsj3qcJxz7oD5fBB16LqRfdhevpuH3/VahHOu8avuOYj2km4AFPN6DzP7Q1Ija4R6F+Ty1QGdefyDxVw9vDdtclpEHZJzzu236moQfwEKgHYxr2MXF8f1I/uwbecuxnotwjnXyFVZgzCz2+ozkKaib4c8Tuvficf+u4grj+9Fa69FOOcaqeruYsqSdLmkKyRl12dQjd0NI/tSumMXj7y/MOpQnHNuv1XXxPQEsA4oAZ6sn3Cahn4d8zj18I48+v4iNm7dGXU4zjm3X6pLEK2BT8OlTf2E03RcP7Ivm7eX82evRTjnGqnqEsTFwFXANcA36iecpuPQzi05+dAO/OX9hWzc5rUI51zjU91QG8vN7Mdm9iMzW1SPMTUZN4zqy+aych59f1HUoTjnXK1V10k9urrRWsPRXo9LTlhNw+FdWnHiIe3583vFbC7zWoRzrnGprompC/CRpDGSrpb0NUkXSvqZpDcIhgBfV93BJX1f0ixJn0gaJylTUhtJr0qaF/5sXUXZUyTNlTRf0k37f4nRumFUXzaVlfPYfxdFHYpzztVKdU1MdwNFwAtAN4KxmY4jSApXhEN/z62qvKQuwA1AkZkdDqQCo4GbgNfNrC/wevi+ctlU4H7gVOBQ4AJJh+7XFUZsQNd8RvQrYOx7C9myvTzqcJxzLmHVjsVkZuVm9rKZ/dTMrjCz68zsfjNL9NacNCBLUhqQDawAzgQeC7c/BpwVp9xgYL6ZFZvZDuDvYblG6YZRfdmwdSd//WBR1KE451zCkjZYn5ktJxgFdgmwEthoZhOADma2MtxnJdA+TvEuwNKY98vCdfuQdJWkKZKmlJSU1OUl1JkjC1sz/KACxr67kFKvRTjnGomkJYiwb+FMoCfQGciRdHGixeOss3g7mtkYMysys6KCgoY7RNR3R/VlfekO/jZxcdShOOdcQpI53PeJwEIzKzGzncDzBH0YqyV1Agh/rolTdhlBv0eFrgTNU43WoO6tGdanHWPeKWbbjl1Rh+OcczWqMUFIKpD0kKR/h+8PlfTNBI69BDhGUnZ4u+woYDbwIl/MZX0p8M84ZScDfSX1lNSCoHP7xQTO2aB998S+rCvdwROTvBbhnGv4EqlBPAq8zRff6OcBP6ypkJlNAp4FpgEzw3ONAe4ATpI0DzgpfI+kzpLGh2XLgeuAVwiSytNmNivhq2qgju7RhmN7teXBt70W4Zxr+GQWt2n/ix2kyWZ2tKSPzOzIcN10MxtYLxHWQlFRkU2ZMiXqMKo1edF6znvwA64Y1pPbTm+Ud+4655oQSVPNrCjetkRqEKWS2hB2Eks6Gthch/E1K0f3aMPFxxTyyPsLmVhc7XOGzjkXqUQSxI+AfwG9JL0NjAOuT2pUTdzNpx5CYZtsfvTMDH94zjnXYNWYIMxsCjAC+BLwXeBQM5ue7MCaspyMNO467wiWb9jG7S/Njjoc55yLK5G7mL4GZJjZDOAU4G+SGlz/Q2NzdI82XHV8L8Z9uIQ358a709c556KVSBPTL8xsczhy61eBp4AHkxtW8/D9kw7ioA653Pjsx2zYuiPqcJxzbi+JJIiK+zFPB/5kZs8BGckLqfnITE/l7vMGsr50Bz9/sdHfxeuca2ISSRArJd1P8LDa+PDBtWQ+gd2s9O/aiutG9uGf01cwfubKqMNxzrk9EvmgP5/gQbnTzOxzoB1xhuh2++/aEX3o36UVt74wk5LN26MOxznngMTuYtoCzAFGSvo20M7MXk56ZM1IemoK95x/BKU7dnHz8zOp6eFF55yrD4ncxXQrwbMPXQgGzXtS0s3JDqy56dshjx+f3I/XZq/muWnLow7HOecSamK6GDjazG41s1sJJvP5RnLDap4uH9aTwT3a8MsXZ7Fiw7aow3HONXOJJIjFBDPDVUgDipMTTvOWmiLuOu8Idpnxk2c/Zvdub2pyzkUnkQSxFZglaaykhwlGZt0g6R5J9yQ3vOansG02t37lEN6bv9aHBXfORSqt5l14KVwqTExSLC504eBCXpm1ml+Pn8PxfQvo0S4n6pCcc81QtQlCUiow3MwurW4/V7ckcec5Azj5d2/zw2dm8PTVx5KaEm8WVuecS55qm5jMbBfQSVJ6PcXjQh1bZfLLMw9j6uLPefhd7/JxztW/RJqYioF3Jf0TKK1YaWZ/qK6QpH4E4zZV6AX8DDgW6Beuywc2xJt8SNIignkndgHlVU1o0ZSdNbALr3yymnsmfMaIfu3p1zEv6pCcc81IIp3UJcCrQDZQELNUy8zmmtnA8MN/EEFn9wtm9vWY9c8Bz1dzmBHhvs0uOUDQ1HT72YeTl5nGD56ezo7y3VGH5JxrRmqsQZjZbQCSssL3+3OD/ihggZntuS1HkgiG8Ri5H8drNtrmZvDrr/Xn6senct+b8/nBSQdFHZJzrplI5EnqQyVNBuYB8yVNknRILc8zmuBp7FjHA6vNbF4VZQyYIGmqpKtqeb4m5cuHdeRrR3bh/jfnM2PphqjDcc41E4k0MY0BbjGzrmbWBbgVeDjRE4Sjv54BPFNp0wXsmzRiDTWzo4BTgWslDa/i+FdJmiJpSklJSaJhNTo/P+MwCnIz+OEzMyjbuavmAs45d4ASSRB5ZvZqxRszew2oTW/pqcA0M1tdsUJSGvA19u7E3ouZrQh/rgFeIBjiI95+Y8ysyMyKCgpq7BpptFplpXPnuQOYv2YLd70yN+pwnHPNQCIJYpGkmyV1DZebCIbfSFS8msKJwBwzWxavgKQcSXkVr4GTgU9qcc4mafhBBVx8TCF/fn8hk4rXRR2Oc66JSyRBXA50A8aHS1fgskQOLikbOIl971Tap09CUmdJ48O3HYD3JM0APgReMrP/JHLOpu7mUw+hW+tsfvTsDLZsL486HOdcE6aq5h6QlAHkmtm6SuvbAZvNrMHNbFNUVGRTpkyJOoykm7xoPec/9AEXDC7k12f3jzoc51wjJmlqVY8SVFeDuJf4t6CeBvggfRE6ukcbrjy+F09OWsJbc9dEHY5zromqLkEMN7PKdx4BPA6ckJxwXKJ+cNJB9G2fy43PfczGrTujDsc51wRVlyDijg5nQZuUjxwXscz0VO45fyDrtuzg5y82+/5751wSVJcg1koaVHmlpKOA9ckLySWqf9dWXDeyD/+YvoKXZ66MOhznXBNT3VAbPwaekzQWmBquKyK4q+nCZAfmEnPtiD68PnsNt/7jE4p6tKEgLyPqkJxzTUSVNQgzmwgcA2QB14RLFnCcmX1QP+G5mqSnpnD3+UewZXs5t74wk6ruSnPOudqqdrA+M1tFMLSGa8AO6pDHj04+iF+Pn8Pz05ZzzqCuUYfknGsCEnlQzjUCVwzrxdE9WvOLf81ixYb9GXDXOef25gmiiUhNEXeddwS7dhs3PvexNzU55w5YlQlC0qPhz+vqLRp3QLq3zeGW0w7h3Xlr+dukJVGH45xr5KqrQQyW1AW4UlKepJaxS30F6GrnoiGFHN+3Hb9+aTaL1pbWXMA556pQXYIYC7wFHAzMqrT4k1kNlCTuPHcAaaniR8/MYNdub2pyzu2f6m5zvcfM+gJ/NbNCM+sWsxTWY4yuljq1yuKXZxzGlMWfM/bd4qjDcc41UjV2UpvZlZIOl3RNuBxaH4G5A3P2kV348mEduHvCZ3y2enPU4TjnGqFE5qS+FngaKAyXZyR9J9mBuQMjidvP7k9eZho3jPuITWU+oJ9zrnYSuc31amCwmd1iZrcAQwieqnYNXLvcDH739YHMX7OFKx6dzLYdPpe1cy5xiSQIAbFfP3fio7k2GsMPKuD3owcydfHnXPO3qewo3x11SM65RiKRBPE4MFHSTyX9FPgv8FhNhST1kzQ9Ztkk6XuSfiFpecz606oof4qkuZLmh/Ngu/10+oDO/N/X+vP2ZyV876mPKN/lScI5V7Nqx2ICMLM7Jb0JHE9Qc7jGzCYnUG4uMBBAUiqwHHiBYD7r35nZXVWVDfe/n2A+62XAZEkvmtmnNV+Si+frRxeyuaycX700m5wWM/nNOQNISfGKoHOuajUmCIAwIdSYFKoxClhgZoulhD6UBgPzzawYQNLfgTMBTxAH4FvH92JTWTl/eH0eeZnp3Hb6IST47+Gca4bqayym0cC4mPfXSfpY0iOSWsfZvwuwNOb9snDdPiRdJWmKpCklJSV1F3ET9f0T+3LZ0B488v5Cfv/avKjDcc41YElPEJJaAGcAFfNbPwD0Jmh+WgncHa9YnHVxHwk2szFmVmRmRQUFBXUQcdMmidu+cijnDerKva/P8wfpnHNVSqiJSVJXoK+ZvSkpA0gzs0QH+jkVmGZmqwEqfobHfRj4d5wyy4BuMe+7AisSPJ+rQUqK+L+v9WfL9qBPomVmOucf3a3mgs65ZiWRB+UuB14kGJsJoDvwz1qc4wJimpckdYrZdjbxx3WaDPSV1DOsgYwOY3B1JC01hd+PHsjwgwq46fmPeeljn9PaObe3RJqYbiCYenQTgJl9BrRP5OCSsgnuRHo+ZvWdkmZK+hgYAXw/3LezpPHhOcqB64BXgNnA02Y2K6ErcgnLSEvloYsHMah7a7731Ee8OXdN1CE55xoQ1TSxjKSJZnaMpI/M7MjwFtTpZta/fkJMXFFRkU2ZMiXqMBqdTWU7uWDMROav2cJfLx/MkF5tow7JOVdPJE01s6J42xKpQbwv6SdApqQRwFPE7zdwjVTLzHT+evlgurbO4orHpjBz2caoQ3LONQCJJIifAJuBOcB3gdeBW5MZlKt/bXMz+Nu3htAqK51vPDKJeT4CrHPNXrUJImxOesTMHjCzs83srPC1j9XQBHVqlcUT3xpCWmoKF/95EkvXb406JOdchKpNEGa2C+gkKb2e4nER69Euh8evGEzZzt1cNHYSqzeVRR2Scy4iiTQxFQPvSrpZ0g0VS7IDc9E5uGNLHrt8MOu2bOeSP0/i89IdUYfknItAIgmiBHgVyAYKYhbXhA3sls/DlxaxaN1WLv3Lh2z2CYeca3ZqvM21MfHbXOvea5+u5pq/TWVQ99Y8dvlgMtNTow7JOVeHDug2V0mvSppQean7MF1DdOKhHbj7/CP4cNF6vvPENJ9wyLlmJJGxmH4a8zoTOAfYnpxwXEN05sAubNlezq0vfMIPnp7OvaOPJNXnknCuyUtkwqBJlVa9LentJMXjGqiLhnRnc1k5d7w8h7zMNH59dn+fS8K5Jq7GBCGpZczbFGAQ0KmK3V0Tds2XerO5bCf3v7mAvMx0bj71YE8SzjVhiTQxzSKYi0FAObAQuDKZQbmG60cn92NzWTlj3ikmLyON60f1jTok51ySJJIgepnZXvc4SkpoHgnX9EjiF189jC1l5dz96mfkZabxzaE9ow7LOZcEiTwHUbkPAuDDug7ENR4pKeLOcwdw8qEd+MW/PuXZqcuiDsk5lwRV1gQktSfoa8iS1J8vpgFtSfDQnGvG0lJT+OOFR3LFo1P4ybMzyM1I5ZTDvWvKuaakuhrEV4D7CKb7/BNwf7jcAtyW/NBcQ5eRlspDlwxiYLd8rh/3Ee98VhJ1SM65OpTIhEHnm9nTtT6w1I9g7ogKvYCfAV2ArwI7gAXAZWa2IU75RQTDjO8Cyqt60i+WP0kdjY1bdzL64YksWlvK41cMpqhHm6hDcs4lqLonqRMaakPSl4HDCB6UA8DMfl2LAFKB5cAQoB/whpmVS/pNeKwb45RZBBSZ2dpEz+MJIjolm7dz/kMfsHbLdv52xRCO6JYfdUjOuQQc6FAbfwIuBX4AZAEXA31qGcMoYIGZLTazCeGc0wATCZqwXCNXkBdMONQyM53zHvqApyYviTok59wBSuQupmFmdiGwzsxuI6gF1PZDfTQwLs76y4GXqyhjwARJUyVdVdWBJV0laYqkKSUl3gYepS75Wbx43VAG92jDjc/N5MfPzKBs566ow3LO7adEEkTFjDFlkjqG73skegJJLYAzgGcqrb+V4MG7J6ooOtTMjgJOBa6VNDzeTmY2xsyKzKyooMBHIY9a29wMHrt8MNeP7MMzU5dx9p/+y6K1pVGH5ZzbD4kkiPGS8oG7gOnAIuDZWpzjVGCama2uWCHpUuB04CKrohPEzFaEP9cALwCDa3FOF6HUFPHDk/vxyDeLWLFhG1+97z0mzFoVdVjOuVqqaU7qFOBlM9tgZs8APYH+ZnZLLc5xAUq2CgUAABgSSURBVDHNS5JOAW4EzjCzuJMeS8qRlFfxGjgZ+KQW53QNwMiDO/Dv64fRo20OVz0+lf97eTblu3y4cOcai5rmpN4N3BvzfpuZrU/04JKygZOA52NW3wfkAa9Kmi7pwXDfzpLGh/t0AN6TNIPgqe2XzOw/iZ7XNRzd2mTzzDXHcuGQQh56u5iL/zyJNZt9nmvnGoNEnoP4X2CKmf2zfkLaf36ba8P23NRl3PqPmbTMTOe+C49icE9/XsK5qB3Qba7AdcALkrZJWi/pc0kJ1yKcq3DOoK688J2hZLdI5YKHJ/LwO8U0pSlvnWtqEkkQ7YB0IBcoCN/77UJuvxzSqSUvXj+MEw9pz+3jZ/Ptv01jU9nOmgs65+pdjQnCzHYB5wE3hq87AQOTHZhrulpmpvPgxYO49bRDeHX2as68733mrNoUdVjOuUoSeZL6PmAEcEm4aivwYDKDck2fJK4c3osnvzWELdvLOev+93l+mg8b7lxDkkgT03FmdjXhA3PhXUwtkhqVazaG9GrLSzcMY0DXfH7w9AxueWGmP33tXAORSILYGT4PYQCS2gJ+M7urM+3zMnnyW0O4+ku9eHLSEs578AOWro/7iIxzrh4lkiDuB54DCiT9EngP+E1So3LNTlpqCjefeggPXTKIRWtLOf2P7/HmnDVRh+Vcs5ZIJ/VfgZ8SDLWxHjjPzP6e7MBc8/Tlwzryr+uH0alVJpc9Opl7Jsxl126/Fda5KCRSgwBIBXYSTPKTaBnn9kuPdjm88J2hnDuoK394Yz7f/MuHrNuyPeqwnGt2ErmL6VaCsZQ6Ewzz/aSkm5MdmGveslqk8ttzB3DH1/ozaeF6Tv/je0xb8nnUYTnXrCRSG7gYONrMfmpmtxKMqvqN5IblXHAr7OjBhTz/7eNITRFff+gDHvvvIn/62rl6kkiCWAykxbxPA4qTE45z+zq8Syteuv54hvct4OcvzuKGv0+ndHt5zQWdcwckkQSxFZglaaykh4GZwAZJ90i6J7nhORdolZ3Ow98o4sdf7sdLH6/gzPvfZ/6azVGH5VyTllbzLrwULhUmJikW56qVkiKuHdGHgd3yuWHcR5xx3/v85pwBfPWIzlGH5lyTVONw342JD/fdfKzaWMa1T05j6uLPGXVwe24+7RD6tM+NOiznGp0DGu5b0imSJkta48N9u4aiY6tM/n7VMdx4ysF8uHA9X/79O9z2j0/8dljn6lAiEwbNB84n6HvYM8RGOLJrg+I1iOZp3Zbt/P61eTz54RKy01O5dmQfvnlcDzLTU6MOzbkG70AnDFoGTDeznWa2q2JJ4KT9wilFK5ZNkr4nqY2kVyXNC3+2rqL8KZLmSpov6aYE4nTNVNvcDP73rMN55XvHM7hnG+54eQ6j7n6bF2es8FtinTsAidQgBgM/B94C9tTfzewPCZ9ESgWWA0OAa4H1ZnZH+MHf2sxujLP/ZwTzWS8DJgMXmNmn1Z3HaxAO4P35a/nVS7OZvXITA7vlc9vphzCou09v6lw8B1qD+CWwC8gnmEmuYqmNUcACM1sMnAk8Fq5/DDgrzv6DgflmVmxmO4C/h+Wcq9HQPu349/XDuPPcAazYsI1zHviA7zwxlcXrSqMOzblGJZHbXNub2aADPM9oguE6ADqY2UoAM1spqX2c/bsAS2PeLyOofexD0lXAVQCFhYUHGKZrKlJTxPlF3Th9QCfGvFPMQ28X8+qnq7n02B5cP7IvrbLTow7RuQYvkRrE65JG7u8JJLUAzgCeqU2xOOvitoWZ2RgzKzKzooICnyrb7S27RRrfO/Eg3vrxCZx9ZBf+/P5CvnTXmzzy3kJ2lPu0Js5VJ5EEcSXwmqQt+3mb66nANDNbHb5fLakTQPgz3qD/y4BuMe+7AitqcU7n9tKhZSZ3nnsEL11/PId3bsX//PtTvvz7d3hl1irvyHauCokkiHZAOtCKoO+hHbXrg7iAL5qXAF4ELg1fXwr8M06ZyUBfST3DGsjosJxzB+TQzi15/IrB/OWbR5OaIq5+fCpfHzORmcs2Rh2acw1OIhMG7QLOA24MX3cCBiZycEnZBHciPR+z+g7gJEnzwm13hPt2ljQ+PGc5cB3wCjAbeNrMZiV6Uc5VRxIjDm7Pf757PL8663AWrNnCV+97jx88NZ0VG7ZFHZ5zDUYit7neR1CDGG5mh0hqA7xiZkfXR4C14be5uv2xuWwnf3prAX9+byECrjy+F9ec0JvcjETu4XCucTvQ21yPM7OrgTIAM1sPtKjD+JyLVF5mOjeecjBv/PBLnHJ4R+57cz4n/PYtnpy0hPJd3pHtmq9EEsROSSmEdxFJakvMkBvONRVdW2dz7+gj+ce1Q+nRNptbXpjJaX94l7fmxruPwrmmr8oEIamifn0/8BxQIOmXwHvAb+ohNuciMbBbPs9ccywPXHQU28t3882/TOaSP09izqpNUYfmXL2qsg9C0jQzOyp8fRhwIsHzCa+Z2Sf1F2LivA/C1bUd5bv56weL+OMb89lctpOzj+zKpcd1Z0DX/KhDc65OVNcHUV2C+MjMjkxqZHXME4RLlg1bd/CH1+cz7sMlbNu5i/5dWnHRkEK+ekRncrwz2zVi+5sglgFVTilqZg1uulFPEC7ZNpXt5B8fLeeJiUuYu3ozuRlpnH1kFy4cUsghnVpGHZ5ztVZdgqjuq08qkEv8YS+ca5ZaZqbzjWN7cMkx3Zm25HOemLiEp6Ys5fGJizmqMJ+LhnTnKwM6+VwUrklIqA+isfAahIvC56U7eG7aMp6ctITitaW0ykrn3EFduXBIIb0LfBpU17B5H4Rz9cDM+KB4HU9MWsIrn6yifLdxTK82XDSkO18+rCMt0hK5q9y5+rW/CaJN+FBco+EJwjUUazaX8cyUZYz7cAnLPt9Gu9wWnDuoGxcOLqSwbXbU4Tm3x34liMbIE4RraHbvNt6ZV8ITk5bw+uzV7DYYflABFw0pZNTB7UlL9VqFi5YnCOcagJUbt/H3D5fy1OSlrNpURoeWGXz96EJGH92NzvlZUYfnmilPEM41IOW7dvPGnDU8MWkJ78wrQcDIgztw0ZBChh9UQGqK3zjo6s/+3ubqnEuCtNQUTj6sIycf1pGl67cy7sMlPD1lKa/NXk2X/CwuHFLIeUVdaZ+XGXWorpnzGoRzDcCO8t1M+HQVT0xcwgfF60hLEScf1oEzB3ZhWJ92/rS2SxqvQTjXwLVIS+H0AZ05fUBnFpRsYdykJTw7bRnjZ66iRWoKQ3q1YeTB7Rl5cHu6t82JOlzXTCS1BiEpHxgLHE4wXPjlwPeAfuEu+cAGM9tnhjpJi4DNwC6gvKoMF8trEK4p2VG+mymL1vPGnDW8MXcNxSWlAPQuyGHkwe0ZcXB7ju7RhnS/E8odgMg6qSU9BrxrZmPDuaWzzWxDzPa7gY1m9j9xyi4CisxsbaLn8wThmrJFa0t5Y84a3py7hknF69mxazd5GWkMP6iAEQe354R+BbTLzYg6TNfIRJIgJLUEZgC9LM5JJAlYAow0s3lxti/CE4RzcZVuL+e9+Wt5Y3aQMNZs3o4EA7rmMypsijqsc0uC/2bOVS2qBDEQGAN8ChwBTAW+a2al4fbhwD1VBiYtBD4naJp6yMzGVLHfVcBVAIWFhYMWL15c15fiXINmZsxasSloipqzhhnLNmAG7fMyGNEvaIoa1redz7Ht4ooqQRQBE4GhZjZJ0r3AJjO7Ldz+ADDfzO6uonxnM1shqT3wKnC9mb1T3Tm9BuEcrN2ynbfmlvDmnDW881kJm7eX7+noHtEvqF30aOcd3S4QVYLoCEw0sx7h++OBm8zsK+F0psuBQWa2LIFj/QLYYmZ3VbefJwjn9rZz126mLPqcN+eu4fXZq1kQdnT3KshhZJgsinq08YEEm7FIbnM1s1WSlkrqZ2ZzgVEEzU0QTF86p6rkICkHSDGzzeHrk4F9OrKdc9VLT03h2N5tObZ3W2457RCWrNvKG3NW88bcEv76wWLGvreQ3Iw0hh/UjhP6taeoe2t6tM0hxZ/mdiT/OYjrgSfCO5iKgcvC9aOBcbE7SuoMjDWz04AOwAthB1sa8KSZ/SfJsTrX5BW2zeabQ3vyzaE9Kd1ezvvz1/Lm3KDvYvzMVQC0zEzjiG75HNktnyPCxe+Oap78SWrnHGbGvDVbmL5kAx8t3cD0pRv4bPVmdu0OPh+6ts5iYLf8PcthnVuR1cJnzWsK/Elq51y1JHFQhzwO6pDH+Ud3A2DrjnI+Wb6JGWHC+GjJBv798UoAUlPEwR3zOCImafQuyPWBBpsYr0E45xK2ZnMZM5Zu3JM0ZizbwOaycgByM9Lo36UVAwvzOaJrPkcW5tOhpQ842NB5DcI5Vyfa52Vy0qGZnHRoByCYEKl4beleCePhd4opD5umOrbMZGDYjzGwWz4DurbygQcbEf+Xcs7tt5QU0ad9Ln3a53LOoK4AlO3cxacrNzF9yRdJ4z+zgg7wFEHf9nkc0a0V/bu0ondBLr3b59I+L8Of+m6APEE45+pUZnoqRxW25qjC1nvWrS/dwYxlG5i+JEgYEz5dzdNTvrjLPTcjjV4FOfQuyKVXuxx6t8+lV0EOPdrmkJnuneFR8T4I51y9MzNWbiyjuKSU4rVbWLBmC8VrS1mwZgsrNpbt2U8K7qAKEkcuvduHSaQgh4Jcr3XUBe+DcM41KJLonJ9F5/wshvVtt9e2rTvKKS4pZUHJlr1+TixeR9nO3Xv2y8tMo1dBLr3DmkfFz8K22WSkea2jLniCcM41KNkt0ji8SysO79Jqr/W7dxsrN5UFtY2SLSwIk8d/56/j+WnL9+yXIihsk70neQQ/c+nWJov2eZl+K24teIJwzjUKKSmiS34WXfKzGH5QwV7btmwvZ+Ge2sYXyeP9+WvZXv5FrSM1RbTPy6BTq0w6tcoKfuaHP1tl0jk/i3a5GZ5EQp4gnHONXm5GGv27tqJ/171rHbt2Gys2bGNByRaWb9jGyg1lrNi4jVUby/h05SZem716rwQCkJYiOrQMEkbHMGlUJJBOrbLolJ9Ju5yMZjFelScI51yTlZoiurXJplub7LjbzYwNW3eyYmOQPFZuKmPlhm2s3FjGyo3bmLl8IxM+Xc2OSkkkPTVIIp1bZdGxVSad8r94XfGzTU6LRl8T8QThnGu2JNE6pwWtc1pwWOdWcfcxM9aX7mDlxjJWbNjGqk1lrNgQJJCVG8v4aOnnvPxJGTt3WaVjQ6usdNrktKBtTgva5LSgTU4GbcPzfbGuBW1zW9A6u0WDu6XXE4RzzlVDEm1zM2ibm7FPx3mF3buNdaU7WLlxGys2lLFq4zbWb93J+tLtrC/dwbotO1i4tpSpiz9nfekOdlfxdEFOi1Ta5MYkkuwgeexJJJWSS25GWlJv9fUE4ZxzByglRRTkZVCQl8GArtXvu3u3sXHbTtZv3bEneawv3cHnWyteb2dd6Q5WbypjzspNrCvdsU8/SYUWqSm0yWlBtzZZPHPNcXV+XZ4gnHOuHqWkfNGs1bug5v3NjK07drG+dMeeZV3pjrB2EtRSktXX4QnCOecaMEnkZKSRk5FWZWd7svhEtM455+JKaoKQlC/pWUlzJM2WdKykX0haLml6uJxWRdlTJM2VNF/STcmM0znn3L6SXYO4F/iPmR0MHAHMDtf/zswGhsv4yoUkpQL3A6cChwIXSDo0ybE655yLkbQEIaklMBz4M4CZ7TCzDQkWHwzMN7NiM9sB/B04MzmROueciyeZNYheQAnwF0kfSRorKSfcdp2kjyU9Iql1nLJdgKUx75eF6/Yh6SpJUyRNKSkpqdMLcM655iyZCSINOAp4wMyOBEqBm4AHgN7AQGAlcHecsvHu2Yr7aImZjTGzIjMrKihI4J4x55xzCUlmglgGLDOzSeH7Z4GjzGy1me0ys93AwwTNSfHKdot53xVYkcRYnXPOVZK0BGFmq4ClkvqFq0YBn0rqFLPb2cAncYpPBvpK6impBTAaeDFZsTrnnNtXUqcclTQQGAu0AIqBy4A/EDQvGbAIuNrMVkrqDIw1s9PCsqcBvwdSgUfM7PYEzlcCLN7PcNsBa/ezbH1rTLFC44q3McUKjSvexhQrNK54DyTW7mYWt32+Sc1JfSAkTalqXtaGpjHFCo0r3sYUKzSueBtTrNC44k1WrP4ktXPOubg8QTjnnIvLE8QXxkQdQC00plihccXbmGKFxhVvY4oVGle8SYnV+yCcc87F5TUI55xzcXmCcM45F1ezTxCNaVhxSd0kvRkOnT5L0nejjqkmklLDsbj+HXUsNYk3PH3UMVVF0vfDv4FPJI2TlBl1TLHCcdbWSPokZl0bSa9Kmhf+jDcOW72rItbfhn8HH0t6QVJ+lDHGihdvzLYfSTJJ7eriXM06QTTCYcXLgR+a2SHAMcC1DTxegO/yxTDvDV1Vw9M3KJK6ADcARWZ2OMHDpKOjjWofjwKnVFp3E/C6mfUFXg/fNwSPsm+srwKHm9kA4DPg5voOqhqPsm+8SOoGnAQsqasTNesEQSMbVtzMVprZtPD1ZoIPsLij3DYEkroCXyF4mr5BO8Dh6aOQBmRJSgOyaWBjlZnZO8D6SqvPBB4LXz8GnFWvQVUhXqxmNsHMysO3EwnGg2sQqvjdAvwO+AlVDGy6P5p7gkh4WPGGRlIP4EhgUvV7Rur3BH+wu6MOJAHVDU/foJjZcuAugm+KK4GNZjYh2qgS0sHMVkLwZQdoH3E8iboceDnqIKoj6QxguZnNqMvjNvcEkfCw4g2JpFzgOeB7ZrYp6njikXQ6sMbMpkYdS4KqGp6+wQnb7s8EegKdgRxJF0cbVdMk6VaCpt0noo6lKpKygVuBn9X1sZt7gmh0w4pLSidIDk+Y2fNRx1ONocAZkhYRNN2NlPS3aEOqVtzh6SOMpzonAgvNrMTMdgLPA8dFHFMiVleM5hz+XBNxPNWSdClwOnCRNewHxnoTfFmYEf5/6wpMk9TxQA/c3BNEoxpWXJII2shnm9k9UcdTHTO72cy6mlkPgt/rG2bWYL/lVjU8fYQhVWcJcIyk7PBvYhQNtEO9kheBS8PXlwL/jDCWakk6BbgROMPMtkYdT3XMbKaZtTezHuH/t2UEc++sOtBjN+sEEXZCXQe8QvAf7GkzmxVtVNUaClxC8G18ericFnVQTcj1wBOSPiYYkv7XEccTV1jLeRaYBswk+H/coIaFkDQO+ADoJ2mZpCuAO4CTJM0juNvmjihjrFBFrPcBecCr4f+zByMNMkYV8SbnXA275uSccy4qzboG4ZxzrmqeIJxzzsXlCcI551xcniCcc87F5QnCOedcXJ4gXLMgqW3MrcGrJC2Pef/fJJ3zSEljw9e/kPSjZJyninP3l/RofZ3PNU1pUQfgXH0ws3UEzzYg6RfAFjO7K8mnvQX4VTJPICktZlC5PcxspqSukgrNrM5G93TNi9cgXLMnaUv48wRJb0t6WtJnku6QdJGkDyXNlNQ73K9A0nOSJofL0DjHzAMGVBo87VBJb0kqlnRDzL4/COd1+ETS98J1PSrNT/CjMLERHuPXkt4GvivpvLDsDEnvxJzvXzS8YcBdI+I1COf2dgRwCMFwysXAWDMbrGBypuuB7xHMG/E7M3tPUiHBk/iHVDpOEVB5QpeDgREET+jOlfQAMAC4DBhCMHjkpPCD//Ma4sw3sy8BSJoJfNnMllea2GYKwYCDdyZ89c7F8ATh3N4mVwxJLWkBUDGM9kyCD3cIBss7NBgGCYCWkvLCOToqdCIYPjzWS2a2HdguaQ3QARgGvGBmpeE5nweOp+YxwZ6Kef0+8KikpwkG7quwhmC0V+f2iycI5/a2Peb17pj3u/ni/0sKcKyZbavmONuAytOAxh57V3i8eEPOQzDEdGwTcOVjlVa8MLNrJA0hmJxpuqSBYZ9LZhiHc/vF+yCcq70JBIM8AiBpYJx9ZgN9EjjWO8BZ4cisOcDZwLvAaqB9ePdVBsGw03FJ6m1mk8zsZ8BavhjC/iD2beZyLmFeg3Cu9m4A7g9HfU0j+JC/JnYHM5sjqVWcpicq7TctvB31w3DVWDP7CEDS/xDMGLgQmFNNPL+V1JegNvI6UNExPgJ4qbYX51wFH83VuSSR9H1gs5nV+5zcYa3jbWBYvNtgnUuENzE5lzwPsHe/Q30qBG7y5OAOhNcgnHPOxeU1COecc3F5gnDOOReXJwjnnHNxeYJwzjkXlycI55xzcf0/tuFWdmSnNhoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t,T[:-1])\n",
    "plt.xlabel('Time (hours)')\n",
    "plt.ylabel('Temperature of Corpse (°F)')\n",
    "plt.title('Temperature of Corpse vs. Time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As t approaches infinity, T asymptotically approaches the ambient temperature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The time of death was 1.8756842373710174 hours before the time of measurement\n"
     ]
    }
   ],
   "source": [
    "K = 0.275\n",
    "t = -(np.log((98.5 - 65)/(85 - 65)))/-K\n",
    "print('The time of death was',t,'hours before the time of measurement')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": 224,
   "metadata": {},
   "outputs": [],
   "source": [
    "def current_temp(time):\n",
    "    if time <= -3:\n",
    "        temp = 55\n",
    "    else:\n",
    "        if time < -2:\n",
    "            temp = 55\n",
    "        else:\n",
    "            if time == -2:\n",
    "                temp = 58                \n",
    "            else:\n",
    "                if time < -1:\n",
    "                    temp = 58\n",
    "                else:\n",
    "                    if time == -1:\n",
    "                        temp = 60\n",
    "                    else:\n",
    "                        if time < 0:\n",
    "                            temp = 60\n",
    "                        else:\n",
    "                            if time == 0:\n",
    "                                temp = 65\n",
    "                            else:\n",
    "                                if time < 1:\n",
    "                                    temp = 65\n",
    "                                else:\n",
    "                                    if time == 1:\n",
    "                                        temp = 66\n",
    "                                    else:\n",
    "                                        if time < 2:\n",
    "                                            temp = 66\n",
    "                                        else: \n",
    "                                            if time >= 2:\n",
    "                                                temp = 67\n",
    "    return temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAfVklEQVR4nO3de3hddZ3v8fcHGBARRG1gFBkLKjA6h5sBxQtHwOvxAujgFQeYOQd1OCDOYRwQLzhHHAeUI4i0cMDKmaci2FIuA2KxI2GUaxJSKLbQ0oYS2zTBtkkIbZLufM8fa+10U5Jm7c1eSZr1eT3PfvZea6/fb31Xkue7v/uXtX5LEYGZmRXHTpMdgJmZTSwnfjOzgnHiNzMrGCd+M7OCceI3MyuYXSY7gCxmzJgRM2fOnOwwzMx2KC0tLc9GRMO263eIxD9z5kyam5snOwwzsx2KpKdHW++hHjOzgnHiNzMrGCd+M7OCceI3MysYJ34zs4Jx4jczm6LWDQxxUutyugaG6tqvE7+Z2RR1WXsnD/b0c1l7Z137deI3M5uC1g0McWPnegL4Ref6ulb9TvxmZlPQZe2dDKf3SxmOqGvV78RvZjbFlKv9wfQ+WYNR36rfid/MbIqprPbL6ln1O/GbmU0xzb39I9V+2WDAw739del/h5ikzcysSBYddUiu/bviNzMrGCd+M7OCceI3MysYJ34zs4Jx4jczKxgnfjOzgnHiNzMrGCd+M7OCceI3MyuYXBO/pL0lzZO0TNJSScek68+W9ISkxyVdkmcMZlZMpd5Buq5eTKlvsKb2fX19zJkzh76+vppjGBjooqXlswwMdNfcRx7yrvgvB+6KiEOAw4Clko4DTgQOjYi3Aj/IOQYzK6DeRasZbO+ld9Hqmto3NTWxevVqmpqaao5h1aor2djzMKvaf1xzH3nILfFL2gs4FrgOICIGI2Ij8GXg+xExkK7vyisGMyumUu8g/S3rIKC/eV3VVX9fXx9tbW1EBG1tbTVV/QMDXaztnAcEa9fOn1JVf54V/4FANzBH0iOSrpW0B3AQ8B5JD0pqknTUaI0lnSmpWVJzd/fU+YGZ2dTXu2g1lKc1jqi66m9qaiLS9hFRU9W/atWVRAynfZSmVNWfZ+LfBTgSmBURRwD9wPnp+lcB7wD+EbhJkrZtHBHXRERjRDQ2NDTkGKaZTScj1X4pTfylqKrqL1f7pVIpaV4qVV31l6v9iOTGKRFDU6rqzzPxdwAdEfFgujyP5IOgA7g5Eg8Bw8CMHOMwswJ5QbVfVkXVX1ntb21eXdVfWe1v7WPqVP25Jf6I6ASekXRwuuoE4A/ALcDxAJIOAnYFns0rDjMrlsHVvVur/bJSMPh0b6b2HR0dI9X+SPNSiY6Ojswx9PS2jlT7ZRFD9PS0Zu4jT9r2k62unUuHA9eSJPeVwBkkQz4/BQ4HBoHzIuI/ttdPY2NjNDc35xanmdl0JKklIhq3XZ/rHbgiog140U6BU/Pcr5mZjc1X7pqZFYwTv5lZwTjxm5kVjBO/mVnBOPGbmRWME7+ZWcE48ZuZFYwTv9k01f18N6ffdTrPbqrtwvihri7aT/0CW2qcJLG/Z4AFP2yhv2egpvbPbVjPjRedT//GDTW1t7E58ZtNU7MfnU3rulZmL55dU/tnr5rFppYWuq+aVVP75jtWsWZFD813ttfU/oH5N9Cx7HHun39DTe1tbE78ZtNQ9/Pd3LriVoLglhW3VF31D3V10bNgAUTQc/PNVVf9/T0DLL2/EwKW3re26qr/uQ3rWXLPIojg8Xt+46q/zpz4zaah2Y/OZjidHXI4hquu+p+9ahYxnM4lPzxcddXffMcqYjidz344qq76H5h/A8TW/bvqry8nfrNpplztDw0ns0MODQ9VVfWPVPtD6eySQ0NVVf3lan84nSFzuBRVVf3lar+0ZQsApS1bXPXXmRO/2TRTWe2XVVP1V1b7ZdVU/ZXV/tb22av+ymq/cv+u+uvHid9smlnctXik2i8bGh6irastU/tNbW1bq/2RDobY9Mgjmdp3ruwdqfbLhktB51M9mdqvWf7ESLVfVtqyhTVPLsvU3saX63z89eL5+M3MqjfWfPyu+M3MCsaJ38ysYJz4zcwKxonfzKxgnPjNzArGid/MrGCc+M3MCsaJ38ysYJz4zcwKJlPil/QGSe9LX+8uac+M7faWNE/SMklLJR1T8d55kkLSjNpCNzOzWoyb+CX9D2AecHW66vXALRn7vxy4KyIOAQ4DlqZ97g+8H1hdbcBmZvbSZKn4zwLeBfQCRMRyYJ/xGknaCzgWuC5tNxgRG9O3/w/wNWDqTxRkZjbNZEn8AxExWF6QtAvZEvaBQDcwR9Ijkq6VtIekjwN/jIjF22ss6UxJzZKau2u856eZmb1YlsTfJOnrwO6S3g/8Erg9Q7tdgCOBWRFxBNAPXARcCHxrvMYRcU1ENEZEY0NDQ4bdmZlZFlkS//kklftjwBeBO4FvZGjXAXRExIPp8jySD4IDgMWS2kn+X9Aq6c+rjNvMzGq0y/belLQzcH1EnAr832o6johOSc9IOjgingBOAFoj4oSK/tuBxoio7k7QZmZWs+0m/ogoSWqQtGvlOH8VzgbmStoVWAmcUUuQZmZWP9tN/Kl24PeSbiMZpwcgIi4br2FEtAEvuvtLxfszM+zfzMzqKEviX5M+dgIyXbhlZmZT17iJPyK+MxGBmJnZxBg38Uv6LaOctx8Rx+cSkZmZ5SrLUM95Fa9fBnwS2JJPOGZmlrcsQz0t26z6vaSmnOIxM7OcZRnqeXXF4k7A2wBfcGVmtoPKMtTTQjLGL5IhnlXA3+UZlJmZ5SdL4v/LiNhcuULSbjnFY2ZmOcsyV899o6y7v96BmJnZxBiz4k8nTtuPZFbOI0iGegD2Al4+AbGZmVkOtjfU80HgdJIZNCunZ+gDvp5jTGZmlqMxE39EXA9cL+mTETF/AmMyM7McZTmPf76kjwBvJbmAq7z+n/MMzMzM8pHlZuuzgU+TTLEs4BTgDTnHZWZmOclyVs87I+JvgA3phG3HAPvnG5aZmeUlS+Ivn8P/vKTXAUMkt080M7MdUJYLuG6XtDdwKdBKchVvVbdhNDOzqWO8e+7uBCyKiI3AfEn/DrwsInomJDozM6u77Q71RMQw8MOK5QEnfTOzHVuWMf6Fkj4pSeNvamZmU12WMf5/APYASpI2kZzSGRGxV66RmZlZLrJcwOUbrJuZTSNZLuCSpFMlfTNd3l/S0fmHZmZmecgyxn8VyUVbn0uXnwN+kqVzSXtLmidpmaSlko6RdGm6/KikBempomZmNkGyJP63R8RZpBdyRcQGYNeM/V8O3BURhwCHAUuBu4G/iohDgSeBC6qO2szMapYl8Q9J2pnkwi0kNQDD4zWStBdwLHAdQEQMRsTGiFgYEVvSzR4gmfbZzMwmSJbEfwWwANhX0sXA74DvZWh3INANzJH0iKRrJe2xzTZ/C/yqmoDNzOylGTfxR8Rc4GskyX4NcFJE/DJD37sARwKzIuIIoB84v/ympAtJbt4+d7TGks6U1Cypubu7O8PuzMwsiywVPyS3Wtw53X73jG06gI6IeDBdnkfyQYCk04CPAp+PiBitcURcExGNEdHY0NCQcZdmZjaeLKdzfgu4Hng1MINk6OYb47WLiE7gGUkHp6tOAP4g6UPAPwEfj4jna47czMxqkuXK3c8CR0TEZgBJ3yeZpfO7GdqeDcyVtCuwEjgDeBjYDbg7nQXigYj4Ug2xm5lZDbIk/naSWy6W5+XfDXgqS+cR0QY0brP6TVmDMzOz+suS+AeAxyXdTXJK5/uB30m6AiAizskxPjMzq7MsiX9B+ii7J59QzMxsImSZpO36iQjEzMwmRpazej6aXoC1XlKvpD5JvRMRnJmZ1V+WoZ4fAZ8AHhvrnHszM9txZLmA6xlgiZO+mdn0kKXi/xpwp6QmkjN8AIiIy3KLyszMcpMl8V9MMgf/y8g+HbOZmU1RWRL/qyPiA7lHYmZmEyLLGP9vJDnxm5lNE1kS/1nAXZI2+3ROM7MdX5YLuPaciEDMzGxiZLmAS5JOlfTNdHl/SUfnH5qZmeUhy1DPVcAxwOfS5eeAn+QWkZmZ5SrLWT1vj4gjJT0CEBEb0vn1zcxsB5Sl4h+StDPJlMxIagCGc43KzMxykyXxX0EyLfM+ki4Gfgf8S65RmZlZbrKc1TNXUgvJPXMFnBQRS3OPzMzMcjFu4pf0bxHxBWDZKOvMzGwHk2Wo562VC+l4/9vyCcfMzPI2ZuKXdIGkPuDQ9Ird3nS5C7h1wiI0M7O6GjPxR8S/pFftXhoRe6WPPSPiNRFxwQTGaGZmdTTuUI+TvJnZ9JJljN/MzKaRXBO/pL0lzZO0TNJSScdIerWkuyUtT59flWcMZmb2Qlkmafu3LOvGcDlwV0QcAhwGLAXOBxZFxJuBRemyWd0NDHTR0vJZBga6a2q/bmCIk1qX0zUwVFP7rt7NfOrq++nq21xTe/o6Yc6HoW9dbe3NxpDb6ZyS9gKOBa4DiIjBiNgInAhcn252PXBSNQGbZbVq1ZVs7HmYVe0/rqn9Ze2dPNjTz2XtnTW1v2LRch5uX88Vi1bU1J6mS2D1A9D0r7W1NxtDnqdzHgh0A3MkPSLpWkl7APtGxFqA9Hmfl34YZi80MNDF2s55QLB27fyqq/51A0Pc2LmeAH7Rub7qqr+rdzO/bOkgAuY1P1N91d/XCW1zIYaTZ1f9Vkd5ns65C3AkMCsijgD6qWJYR9KZkpolNXd31/ZV3Ypr1aoriUjmEowoVV31X9beyXAEAMMRVVf9VyxaPtK+FFF91d90SZL0IXl21W91lOl0Tkn7SXqnpGPLjwx9dwAdEfFgujyP5INgnaTXAqTPXWPs95qIaIyIxoaGhmxHY8bWaj8iqdIjhqqq+svV/mCStxmM6qr+crU/VEo6GCpFdVV/udovDSbLpUFX/VZXWf65+33g98A3gH9MH+eN1y4iOoFnJB2crjoB+ANwG3Bauu40fBWw1VlltV9WTdVfWe2XVVP1V1b7ZVVV/ZXVfpmrfqujLDdiORk4OCIGauj/bGBueuOWlcAZJB82N0n6O2A1cEoN/ZqNqae3daTaL4sYoqenNVP75t7+kWq/bDDg4d7+TO1bV28cqfbLhkpB69MbMrWn46Gt1X5ZaTBZb1YHim0qkxdtIP0KOCUinpuYkF6ssbExmpubJ2v3ZmY7JEktEdG47fosFf/zQJukRcBI1R8R59QxPjMzmyBZEv9t6cPMzKaBLHfgul7S7sBfRMQTExCTmZnlKMtZPR8D2oC70uXDJfkbgJnZDirLlA0XAUcDGwEiog04IMeYzMwsR1kS/5aI6Nlm3fZPBTIzsykryz93l0j6HLCzpDcD5wD35RuWmZnlJUvFfzbJDJ0DwA1AL3BunkGZmVl+spzV8zxwYfowM7Md3JiJX9KPIuJcSbczyph+RHw818jMzCwX26v4y3fZ+sFEBGJmZhNjzMQfES3pc1M6ydohJJX/ExExOFY7MzOb2sYd45f0EWA28BQg4ABJX4yIX+UdnJmZ1V+W0zl/CBwXESsAJL0RuANw4jcz2wFlOZ2zq5z0UysZ465ZZmY29W3vrJ5PpC8fl3QncBPJGP8pwMMTEJuZmeVge0M9H6t4vQ74r+nrbuBVuUVkZma52t5ZPWdMZCBmZjYxspzVcwDJtA0zK7f3BVxmZjumLGf13AJcB9wODOcbjpmZ5S1L4t8cEVfkHomZmU2ILIn/cknfBhbywputt+YWlZmZ5SZL4v8vwBeA49k61BPpspmZ7WCyJP6TgQM9P4+Z2fSQ5crdxcDetXQuqV3SY5LaJDWn6w6X9EB5naSja+nbxtffM8CCH7bQ3zMw/sajeG7Dem686Hz6N26oqX2pd5CuqxdT6qutZujr62POnDn09fXV1N7MRpcl8e8LLJP0a0m3lR9V7OO4iDg8IhrT5UuA70TE4cC30mXLQfMdq1izoofmO9trav/A/BvoWPY498+/oab2vYtWM9jeS++i1TW1b2pqYvXq1TQ1NdXU3sxGlyXxf5tkuOd7JBO2lR+1CmCv9PUrgTUvoS8bQ3/PAEvv74SApfetrbrqf27DepbcswgiePye31Rd9Zd6B+lvWQcB/c3rqq76+/r6aGtrIyJoa2tz1W9WR+Mm/ohoqnwAW4BPZew/gIWSWiSdma47F7hU0jMkN3m5YLSGks5Mh4Kau7u7M+7OyprvWEUMJzdOi+Gouup/YP4NEMNp++Gqq/7eRash0hu3RVRd9Tc1NRFp+4hw1W9WR1kq/vK4/CWS2oHvAksz9v+uiDgS+DBwlqRjgS8DX42I/YGvklwc9iIRcU1ENEZEY0NDQ8bdGWyt9odLSeIcLkVVVX+52i9t2QJAacuWqqr+kWo/3T+lqKrqL1f7pVIpaV4queo3q6MxE7+kgyR9S9JS4ErgGUARcVxEXJml84hYkz53AQuAo4HTgJvTTX6ZrrM6qqz2y6qp+iur/a3ts1f9L6j2RzrIXvVXVvtbm7vqN6uX7VX8y4ATgI9FxLsj4sdAKWvHkvaQtGf5NfABYAnJmH55ps/jgeW1BG5j61zZO1Ltlw2Xgs6nejK1X7P8iZFqv6y0ZQtrnlyWqf3g6t6t1f5IB8Hg072Z2nd0dIxU+yPNSyU6OjoytTez7dO2ldXIG9LJwGeAdwJ3Ab8Aro2IAzJ1LB1IUuVDcr3AzyPiYknvBi5P120G/r58f9+xNDY2RnNzc5bdmplZSlJLxRmVI7Y3LfMCYEFarZ9EMh6/r6RZwIKIWLi9HUbESuCwUdb/DnhblfGbmVmdZDmrpz8i5kbER4HXA23A+blHZmZmuch0Vk9ZRKyPiKsjwvP0mJntoKpK/GZmtuNz4jczKxgnfjOzgnHiNzMrGCd+M7OCceI3MysYJ34zs4Jx4jczKxgnfjOzgnHiNzMrGCd+M7OCceI3MysYJ34zs4Jx4jczKxgnfjOzgnHiNzMrGCd+M7OCceI3MysYJ34zs4Jx4jczKxgnfjOzgnHiNzMrmFwTv6R2SY9JapPUXLH+bElPSHpc0iV5xmBmZi+0ywTs47iIeLa8IOk44ETg0IgYkLTPBMRgZmapyRjq+TLw/YgYAIiIrkmIwcyssPJO/AEslNQi6cx03UHAeyQ9KKlJ0lGjNZR0pqRmSc3d3d05h2lmVhx5D/W8KyLWpMM5d0talu7zVcA7gKOAmyQdGBFR2TAirgGuAWhsbAzMzKwucq34I2JN+twFLACOBjqAmyPxEDAMzMgzDjMz2yq3xC9pD0l7ll8DHwCWALcAx6frDwJ2BZ4dqx8zM6uvPId69gUWSCrv5+cRcZekXYGfSloCDAKnbTvMY2Zm+ckt8UfESuCwUdYPAqfmtV8zM9s+X7lrZlYwTvxmZgXjxG9mVjBO/GZmBePEb2ZWME78ZmYF48RvZlYwTvxmZgUzrRN/V+9mPnX1/XT1ba6tg75OmPNh6FtXU/Pu57s5/a7TeXZTbTNSDHV10X7qF9ji2UnNrI6mdeK/YtFyHm5fzxWLVtTWQdMlsPoBaPrXmprPfnQ2retamb14dk3tn71qFptaWui+alZN7c3MRjNtE39X72Z+2dJBBMxrfqb6qr+vE9rmQgwnz1VW/d3Pd3PrilsJgltW3FJ11T/U1UXPggUQQc/NN7vqN7O6mbaJ/4pFyxlO534rRVRf9TddkiR9SJ6rrPpnPzqb4bT9cAxXXfU/e9UsYng43f2wq34zq5tpmfjL1f5QKUn8Q6WoruovV/ulwWS5NFhV1V+u9oeGh5L9Dw9VVfWPVPtDQ+mKIVf9ZlY30zLxV1b7ZVVV/ZXVflkVVX9ltV9WTdVfWe1v3b2rfjOrj2mZ+FtXbxyp9suGSkHr0xuyddDx0NZqv6w0mKzPYHHX4pFqf2T/w0O0dbVlar+prW1rtT/SwRCbHnkkU3szs+3RjnAPlMbGxmhubp7sMMzMdiiSWiKicdv107LiNzOzsTnxm5kVjBO/mVnBOPGbmRWME7+ZWcHsEGf1SOoGnq6x+QygtlnSdlw+5mLwMRfDSznmN0REw7Yrd4jE/1JIah7tdKbpzMdcDD7mYsjjmD3UY2ZWME78ZmYFU4TEf81kBzAJfMzF4GMuhrof87Qf4zczsxcqQsVvZmYVnPjNzAqmEIlf0v+W9KikNkkLJb1usmPKm6RLJS1Lj3uBpL0nO6a8STpF0uOShiVN21P+JH1I0hOSVkg6f7LjmQiSfiqpS9KSyY5lIkjaX9JvJS1N/6a/Us/+C5H4gUsj4tCIOBz4d+Bbkx3QBLgb+KuIOBR4ErhgkuOZCEuATwD3TnYgeZG0M/AT4MPAW4DPSnrL5EY1IX4GfGiyg5hAW4D/FRF/CbwDOKuev+dCJP6I6K1Y3AOY9v/RjoiFEbElXXwAeP1kxjMRImJpRDwx2XHk7GhgRUSsjIhB4BfAiZMcU+4i4l5g/WTHMVEiYm1EtKav+4ClwH716n+XenU01Um6GPgboAc4bpLDmWh/C9w42UFYXewHPFOx3AG8fZJisQkgaSZwBPBgvfqcNolf0m+APx/lrQsj4taIuBC4UNIFwP8Evj2hAeZgvGNOt7mQ5Gvj3ImMLS9Zjnma0yjrpv032KKS9ApgPnDuNiMXL8m0SfwR8b6Mm/4cuINpkPjHO2ZJpwEfBU6IaXLBRhW/5+mqA9i/Yvn1wJpJisVyJOnPSJL+3Ii4uZ59F2KMX9KbKxY/DiybrFgmiqQPAf8EfDwinp/seKxuHgbeLOkASbsCnwFum+SYrM4kCbgOWBoRl9W9/2lSCG6XpPnAwcAwyfTOX4qIP05uVPmStALYDfhTuuqBiPjSJIaUO0knAz8GGoCNQFtEfHByo6o/Sf8N+BGwM/DTiLh4kkPKnaQbgPeSTFG8Dvh2RFw3qUHlSNK7gf8EHiPJWwBfj4g769J/ERK/mZltVYihHjMz28qJ38ysYJz4zcwKxonfzKxgnPjNzArGid/GJek16cymbZI6Jf2xYvm+nPZ5hKRr09cXSTovj/28FJLOSWdPnCvpvZLeWef+fybpr9PX1071ydgkNUi6a7LjsPFNmyt3LT8R8SfgcEiSMPBcRPwg591+HfhunjuQtEvFRHa1+HvgwxGxqvxzATJ/EFaz/4j477WFWF+Sdo6I0mjvRUS3pLWS3hURv5/o2Cw7V/z2kkh6Ln1+r6QmSTdJelLS9yV9XtJDkh6T9MZ0uwZJ8yU9nD7eNUqfewKHRsTiitVvkXSPpJWSzqnY9h8kLUkf56brZlbO2y7pvDQxk/bxPUlNwFfSOfyXSFos6UXTOUt6haRFklrT4zgxXT8bOBC4TdJXgS8BX02/Bb1nrONMv71cI2kh8P+22ZckXSnpD5LuAPapeO8eSY2Sdk6/CSxJ4/lq+v6bJP0mPY5WSW9M+7u0YttPp9vemF4EVu77Z5I+mfZ9aRrvo5K+WPG7/a2knwOPKbm/xVcq2l9c8Tu5Bfj8KH8qNpVEhB9+ZH4AFwHnVSw/lz6/l+Rq2deSXDH8R+A76XtfAX6Uvv458O709V+QXJK+7T6OA+Zvs8/70n5nkFyN/GfA20iubNwDeAXwOMkshjOBJRXtzwMuSl/fA1xV8d5jwH7p671HiWUXYK/09QxgBVsvfGwHZozxcxn1ONPtWoDdR9nXJ0juo7Az8Lr05/nXFXE3psd8d0WbvdPnB4GT09cvA14OfLKiv32B1env52Tg+nTbXUlm+9wdOBP4Rrp+N6AZOCD93fYDB6TvzQRa09c7AU8Br0mX9wMem+y/Uz+2//BQj9XTwxGxFkDSU8DCdP1jbJ0K+30k1Xu5zV6S9oxkzvGy1wLd2/R9R0QMAAOSukgS2buBBRHRn+7zZuA9jD93TeUU1b8HfibpJmC0ibAEfE/SsSSXzu+X7rtznH2Mepzp69siYtMobY4FbohkKGWNpP8YZZuVwIGSfkwy2eDCtN/9ImIBQERshpHL/sv9rUu/5RwF/Aq4QtJuJDc3uTciNkn6AHBo+f8KwCuBNwODwEMRsSrtv13SnyQdkf4sHolkOBCgi+RDy6YwJ36rp4GK18MVy8Ns/VvbCThmjMRXtomkah2r71La32hTFEMyDXXlMOa2ffWXX0TElyS9HfgI0Cbp8IokBsmwRQPwtogYktQ+Sn+jGfU40w+C/lFbpCFtr9OI2CDpMOCDwFnAp4Bzx9h81J9PRGyWdE/ax6eBGyq2Pzsifr1NzO8dJeZrgdNJpsj+acX6l5H8/mwK8xi/TbSFJPdDAEDS4aNssxR4U4a+7gVOkvRySXuQDGH8J8kkXvsoORtpN5KpqUcl6Y0R8WBEfAt4lhdOeQxJ1duVJv3jgDeM0VUfsGfFcpbjHO14PpOOtb+WUW4YJGkGsFNEzAe+CRwZyTztHZJOSrfZTdLL0/4+nfbXQPKN4qG0q18AZ5B8Qyon+l8DX1YyHTCSDkp/rqNZQPJt4aiK9gAHkdwC06YwV/w20c4BfiLpUZK/v3tJ/jE6IiKWSXrlKENAbLNdq6SfsTWZXRsRjwBI+meSce9VbH8a7kuVTNstYBGweJv35wK3S2oG2rbT1+3AvPSfv2dnOc5RLACOJxkaexJoGmWb/YA5kspFW/leyl8Ark6Pewg4Je3vmPSYAvhaRJSHqMr/XL4tkls4QlLFzwRalXw16QZOGi3QiBiU9FtgY7zwLJ/jSIagbArz7Jw2JaVnq/RFxLWTHYu9WPrB0wqcEhHLK9bfC5wYERsmLTgbl4d6bKqaxQvH9W2KUHIh2Qpg0TZJvwG4zEl/6nPFb2ZWMK74zcwKxonfzKxgnPjNzArGid/MrGCc+M3MCub/A8yaG8DDReGgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in np.linspace(-3,2,20):\n",
    "    plt.plot(i,current_temp(i),'^')\n",
    "    plt.xlabel('Time (hours after discovery)')\n",
    "    plt.ylabel('Ambient temperature')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the given temperature values , my plot of ambient temperature appears to be correct. In order to better visualize T_a(t) would be to interpolate between our values at each hour. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Disregard these 3 cells below that are formatted as markdown"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "F = 10\n",
    "T_a=[55]*F + [58]*F + [60]*F + [65]*F + [66]*F + [67]*F\n",
    "T_i= 85\n",
    "K=.275\n",
    "N=len(T_a)\n",
    "T=[T_i]*(N+1)\n",
    "t=np.linspace(0,12,N)\n",
    "for i in range (0,len(t)):\n",
    "    slope= -K/F*(T[i]-T_a[i])\n",
    "    T[i+1]=T[i] + slope \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "plt.plot(t,T[:-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "F = 10\n",
    "T_a=[55]*F + [58]*F + [60]*F + [65]*F + [66]*F + [67]*F\n",
    "T_i= 85\n",
    "K=.275\n",
    "N=len(T_a)\n",
    "T=[T_i]*(N+1)\n",
    "t=np.linspace(0,12,N)\n",
    "for t in range (0,len(t)):\n",
    "    slope= K/F*(T[i-1]-T_a[i-1])\n",
    "    T[i-1]=T[i] + slope\n",
    "    if T==98.6:\n",
    "        break\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {},
   "outputs": [],
   "source": [
    "def varyT_a(N):\n",
    "    t = np.linspace(0,12,N)\n",
    "    T_num = np.zeros(len(t))\n",
    "    T_num[0] = 85\n",
    "    delta_t = t[1] - t[0] \n",
    "    K = 0.275\n",
    "    \n",
    "    for i in range (1,len(t)):\n",
    "        T_num[i] = ((-K*(T_num[i-1] - current_temp(t[i])))*(delta_t)) + T_num[i-1]\n",
    "    \n",
    "    T_ana = 65 + (T_num[0] - 65)*np.exp(-K*t)\n",
    "    \n",
    "    plt.plot(t,T_num, 'o-', label = 'T_Numerical with varying T_a')\n",
    "    plt.plot(t,T_ana, '-', label = 'T_Analytical')\n",
    "    plt.plot(t,[current_temp(t[i]) for i in range(0,len(t))],'o',label='ambient Temp')\n",
    "    plt.legend(loc='best')\n",
    "    plt.xlabel('Time (hours)')\n",
    "    plt.ylabel('Temperature of Corpse (F)')\n",
    "    plt.title('Temperature of Corpse vs. Time')\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeZzN9f7A8dd7Fox1bGHIRYuyDiYKSbKFkCXaqUiXSP3qcluobuXSnm6lrhYVSbbCFaGiwtgjSpYwiNHIMpjl/fvj+z3TmXFm5sxyZn0/H4/zcL7L5/v9fGfGeZ/PLqqKMcYYk1ZQfmfAGGNMwWQBwhhjjE8WIIwxxvhkAcIYY4xPFiCMMcb4ZAHCGGOMTxYgjMlFItJeRH4VkZMi0jW/81NYiciTIjI5v/NR3FmAKAbcDyvPK1lE4r22b83v/OWEiBwSkbb5nQ8vzwATVbWsqv7P1wkicqeIrBeRUyJyUES+EJEr8zif+UZESvr4mzzttd1XVcep6oj8zmtxF5LfGTCBp6plPe9FZA9wj6ouzb8c+UdEQlQ1sZDd42/A1gzu90/gfuBeYCmQCHQDegE/ZOVGefHzCQRVPQt4/00eAvqp6sr8y5XxxUoQBhEJFpHHRWSXiBwVkY9EJNw9dpmIJIrI3SJyQERiReQuEblKRH4UkTgRedHrWsNEZJmIvCUif4rINhFp53W8koh84H7z3yci40QkKE3a10XkD2CMe/8VInJMRI6IyPsiUs49/1PgAuBL95vnSBHpKiI70zxfSilDRCaIyMci8omInAAGZvT86fy8hrvVSLEiMltEqrn79wMRnvz4SFcZeAIYqqrzVfW0qp5T1bmqOtY9J8x9/oMisl9EJolIqHusq4jsdPN6GHjDa9+T7s9ol4j097pnLxHZLiIn3J/3SK9jN4rIZvd3+K2INEjned8TkX+l2bdYRP7uvn/cze+fIvKTiFyd3s/OX+7v6R33fZb+Bt0094rIDvdnskBEauY0T8WSqtqrGL2APUDHNPvGAN/ifLiVAt4D3nWPXQYo8ApQEugJnAI+AyoDtYE/gFbu+cNwvhX/HQgF7gCOAeXd44uA14DSQA1gA3BnmrRDgGAgzL1/B6AEUB3nW/YEr7wfAtp6bXcFdqZ5vpRzgAnAWZxv7UHuPdJ9fh8/v27u9Zq4504BlqSXnzRpewPxgGTw+5no5qUKUA1YCzzq9WyJwFPuzyPMa99z7r6OwGmgrpsmFmjpvq8MNHPfXwkcBFq4P+uhwM9AiI88dfb+meIE5Xg3j02BXW5eBajnuXcW/ibP+5m5v6d3svk3OBD4CbjU/Rv8F7A8v//vFcZXvmfAXnn8C/cdIHYDbby267ofMuL1n7Oy1/FTQC+v7QXAMPf9MGB3mutvBvrjVL+cAkK9jg0GFnml/TmT/A8Evvfazk6A+NLf5/dx/4+Ap7y2w4FkoLqv/KRJezewJ5PnOwB08NruBWz3era0P7+uwBmglNe++cDD7vvD7s+4XJr7vIsbeLz27fV8yKbZH+w+lyfQ3A8sdN83xAk01+IjuPj5N+lvgPD3b3A5cKvXsVAgAaiW1//fCvvLqpiKORER4EJgoVtUj8P5Vh+E8+0MIElVY72SxeN88Hhvl/Xa3p/mNntxvp3/Dedb9xGve72C8+3TY1+a/EWIyKdu1cKfwDs431xzIuUefj6/twj3eQBQ1TjgT8CfKoxYoJp7z/O4+6t7X999733tQ6qakCbpEVU9kyZNhPu+N9AX+M2tvoty9/8N+Kfnmd3nrurrOVQ1CZgJ3OzuugUnUKKqW3FKYM8Av7vVc9XSXiMXZOVv8G/Am17PdQSnlFUrAPkq0ixAFHPqfMXyfGsN93qVUtWj2bxs2v+ItYEYnA/mk0BFr/uUV9Xm3llKk3YSzrfFRqpaHrgHp2ST3vmncKqvAHDr7yulOSclTTaePwbnA8hz/QpAefcamVnp5r27r4NuXg55Xx/nZ+d9bV/TL1cRkVJp0sS41/xeVXvgBOEvgenuOfuAJ9I8c2lVnZ1O3qcDA0TkYqAxMNcr3++ramuc6qVSOFU6+WkfMCjNs4Wp6rp8zlehYwHCALwJTBCRCwFE5AIRuSEH17vQbXAOEZHbcD6wvlTV3ThtCBNFpJyIBInIJZJxN9VyOEHlTxGpDTyY5vhhnA8mj5+ASiJynRscniTzv/OsPP90YIiINHI/lP8NLFPVQ5ncAzfgPA28JSI93AbpUBG5QUSe9br+OBGpLCIXAI8CH2Zy6VDgcREpISIdgE7AZyJSRkQGikh5nCqWE0CSm2YKcL+IRImjrIj0FJHSvm6gqt/jVGW9AXyuqqfcn1UDEblGRErifIuP97pHfnkTeExE6gOISEUR6ZvPeSqULEAYcBpGlwLLxOnZ8x3QPOMkGfoGaIbTOP0ocKOqHneP3YxTb7/dPf4JqauY0noCaAscB+bgNEx6ewZ4xq1OGOF+CI/CqQLZj/ONPLOSkN/Pr6pf4DQIz8f5ll4duD2T63unfwbnZ/K0m6/fcBqI53k97zacrrIbgVVu/jKyB6cK5RAwFRisqrvcY3fhVDkdx+kwcKebj1XASOAtIA6ngfoWfJdQPKbjNIJ/7LUvDHjBfZaDONU8TwC4vY7y/Fu7qk4HJgOz3WrJjThB02SROKVaY3KHiAzD6dPeMb/zUhyIM1p7sqpenN95MUWPlSCMMcb4ZAHCGGOMT1bFZIwxxicrQRhjjPGpSE3WV6VKFa1Tp05+Z8MYYwqNdevWHVXVqr6OFakAUadOHaKjo/M7G8YYU2iIyN70jlkVkzHGGJ8sQBhjjPHJAoQxxhifilQbhDG5ISEhgf3793PmzJnMTzamkChVqhS1atUiNDTU7zQWIIxJY//+/ZQrV446deqQzszcxhQqqkpsbCz79++nbt26fqcLaBWTiIwWka3usoDTRaSUiIx35/bf6L66pZO2q7tk4E4RGROoPM7dcIA2E5ZRd8wC2kxYxtwN/szabIqyM2fOULlyZQsOpsgQESpXrpzlUnHAShDuGrAjgQaqGi8iM3FWAwN4SVWfzyBtMPA6zgyM+4G1IjJfVbflZh7nbjjA2NlbiE9wZic+EBfP2NlbAOjdzJawLc4sOJiiJjt/04FupA4BwkQkBGcRlxg/07XEWTZyl6qeA2bgLL2YqyYt3kFyQjxDgr+gpfwEQHxCEpMW78jtWxljTKETsAChqgeA53Hmuz8IHFfVL93DI0Rks4hMFZGKPpLXJPXSk/tJZ0lHERkqItEiEn3kyJEs5TEmLh5FuDtkEQ+EfJZqvzHGFHcBCxDuB38vnAXgI4Ay7upibwAXAZE4geMFX8l97PM5q6CqTlHVKFWNqlrV52jxdEWEh3GOUN5O7Ebr4G00k19S9hvjr9xux4qNjSUyMpLIyEiqV69OzZo1U7bPnTt33vkiwkMPPZSy/fzzzzN+/Pgc5SGroqOjGTlyZLbStm/fPlszIDzxxBMsXboUgJdffpnTp0+nHCtbtmx6yQJu/vz5TJgwIcfXGT58OJGRkTRo0ICwsLCUv4FZs2blQi79E8heTB2B3ap6BEBEZgOtVTVl+UQReRv4wkfa/TgLyXvUwv/qKb893KU+Y2dvYXrCdYwImcffQ+Yzkkd4uEv93L6VKaIC0Y5VuXJlNm7cCMD48eMpW7Ys//d//5fu+SVLlmT27NmMHTuWKlWqZOueOZGYmEhUVBRRUVF5et+nnnoq5f3LL7/MbbfdRunSPldMzXWJiYmEhPj++OzZsyc9e/bM8T1ef/11APbs2UOPHj1S/ibyUiADxG/Ale4at/HAdUC0iNRQ1YPuOTcCP/pIuxa4RETq4izYPhBnOcRc5fkPPGnxDt470YXRoZ/xZEuxBmqT4snPt7It5s90j2/4LY5zScmp9sUnJPHIrM1MX/ObzzQNIsoz7oaGuZbHkJAQhg4dyksvvcQzzzyT6tigQYPo0aMH/fr1A5xv1idPnmTFihWMGzeOatWqsXHjRvr06UPjxo155ZVXiI+PZ+7cuVx00UUcOXKEYcOG8dtvzrO8/PLLtGnThvHjxxMTE8OePXuoUqUKQ4cO5fnnn+eLL77g5MmT3H///URHRyMijBs3jr59+3Lfffexdu1a4uPj6devH08++WS6z7RmzRomTJjA7NmzmTdvHgMHDuT48eMkJyfToEEDdu3alfJsMTExxMTEcO2111KlShWWL18OwKOPPsoXX3xBWFgY8+bNo1q1v1a2TU5Opl69emzcuJHw8HAALr74YlatWsWaNWv417/+xblz56hcuTIfffQR1apVO++Z9+3bx2uvvUZkZCQAbdq04Y033mD9+vVER0czefJkBg0aRPny5YmOjubQoUNMnDiRfv36kZyczIgRI/j666+pW7cuycnJ3HXXXSm/p+x4++23mTJlCufOnePiiy9m2rRpOQ6YgWyDWA3MAtYDW9x7TcFZsH6LiGwGrgVGA4hIhIgsdNMmAiOAxTiL0M9U1a2ByGfvZjVZNaYDt496llNakrDVr1qXV+O3tMEhs/2BMnz4cD766COOHz+e+cmuTZs28corr7BlyxamTZvGzz//zJo1a7jnnnt47bXXABg1ahSjR49m7dq1fPbZZ9xzzz0p6detW8e8efP4+OOPU1336aefpkKFCmzZsoXNmzfToUMHAJ555hmio6PZvHkzX3/9NZs3b043b82bN2fDhg0AfPvttzRq1Ii1a9eyevVqWrVqlerckSNHEhERwfLly1OCw6lTp7jyyivZtGkT7dq14+23306VJigoiF69ejFnzhwAVq9eTZ06dahWrRpt27blhx9+YMOGDQwcOJCJE/9aEtz7me+55x7ee+89AH7++WfOnj1LkyZNznuWgwcPsnLlSr744gvGjHF67M+ePZs9e/awZcsW3nnnHb7//vt0fxb+6tOnD2vXrmXTpk1cfvnl/Pe//83xNQM6UE5VxwHj0uz2ucC7qsYA3by2FwILA5e71FYeSOJockcGBy1kovRlX1w16/JqMv2m32bCMg746NRQMzyMT+69KlDZOk/58uW54447ePXVVwkL868N7YorrqBGjRoAXHTRRXTu3BmAxo0bp3zQLl26lG3b/upd/ueff3LixAnAqUrxda+lS5cyY8aMlO2KFZ1+KDNnzmTKlCkkJiZy8OBBtm3b5vMDFZxS0cUXX8xPP/3EmjVrePDBB/nmm29ISkri6quvzvTZSpQoQY8ePQBo0aIFS5YsOe+cAQMG8NRTTzF48GBmzJjBgAEDAGeg5IABAzh48CDnzp1LNbDM+5n79+/P008/zaRJk5g6dSqDBg3ymZfevXsTFBREgwYNOHz4MAArV66kf//+BAUFUb16da699tpMnykzP/74I4899hhxcXGcPHmSLl265PiaNheTa9LiHUxJ6EYiwQwLdppFrMuryczDXeoTFhqcal9YaHC+tGM98MAD/Pe//+XUqVMp+0JCQkhOdkozqpqqkbtkyZIp74OCglK2g4KCSExMBJyqmO+//56NGzeyceNGDhw4QLly5QAoU6aMz3yo6nl97nfv3s3zzz/PV199xebNm+nevXumg7auvvpqFi1aRGhoKB07dmTlypWsXLmSdu3aZfqzCA0NTclDcHBwyvN4u+qqq9i5cydHjhxh7ty59OnTB4D777+fESNGsGXLFt56661U+fR+5tKlS9OpUyfmzZvHzJkzueUW37Xg3j9nzwqegVjJc9CgQUyePJktW7Ywbty4XJkqxgKEKyYunt+pyKdJ19A/eAU1iE3Zb0x6ejeryXN9GlMzPAzBKTk816dxvpQ6K1WqxE033ZSqaqFOnTqsW7cOgHnz5pGQkJCla3bu3JnJkyenbPvTUJo2zR9//MGff/5JmTJlqFChAocPH2bRokWZXqddu3a8/PLLXHXVVVStWpXY2Fi2b99Ow4bnl+rKlSuXUrLxl4hw44038uCDD3L55ZdTuXJlAI4fP07Nms7v7/3338/wGvfccw8jR47kiiuuoFKlSn7fu23btnz22WckJydz+PBhVqxYkaW8+3LixAlq1KhBQkICH330UY6vBxYgUni6tr6R2BMB7g35PNV+Y9LjacfaPaE7q8Z0yNcqyYceeoijR4+mbA8ZMoSvv/6ali1bsnr16nS/9afn1VdfJTo6miZNmtCgQQPefPPNTNM89thj/PHHHzRq1IimTZuyfPlymjZtSrNmzWjYsCF33XUXbdq0yfQ6rVq14vDhwyklhiZNmtCkSROfI4KHDh3K9ddfn+WqmgEDBvDhhx+mVC+B03Osf//+XH311Zn2CmvRogXly5dn8ODBWbpv3759qVWrFo0aNeLee++lVatWVKhQIUvXSOvpp5+mVatWdOrUicsuuyxH1/KQQBR18ktUVJRmd0U57+6KE0KmcGPwKjomv8pDfa6xNohi5qeffuLyyy/P72yYQiAmJob27duzfft2goKy9n375MmTlC1bltjYWFq2bMmqVauoXr16gHLq8PW3LSLrVNVnH2UrQbi8qwreSOpFMEncLZ8z+pON1qPJGHOeDz74gFatWvHMM89kOTgA9OjRg8jISK6++moef/zxgAeH7LAShA9zNxwgafa9dAtaTduzrxBLBcJCg/OtbtnkrYJYgoiNjeW66647b/9XX32VUnduiofhw4ezatWqVPtGjRrlVzVXVksQth6ED5MW76BkYi9uLLGSISELmJB4S0qPJgsQJj94j642xZtnhHVesComH2Li4tmlEcxPvoo7gpdQmeMp+40xpriwAOGDp+fSq4l9KMk57g35ItV+Y4wpDixA+OAZ/LRLI5ib3JY7gr/kwtDjNomfMaZYsQDhg3ePptcSbySEJIbIfOvRZIwpVixApMMz+OmBm7oyO7kdA2QpVfkjZTpnCxImULK6HgTAnDlzEBG2b9+eo3sPGjQo0/UGnn322VTbrVu3zta9xo8fz/PPp7vysCkALEBkYtLiHbya2JsgkhkeMhewOZpMYHl6LG3cuJFhw4YxevTolO0SJUr4TDN9+nTatm2bapK8QEkbIL777ruA39PkD+vmmglnWdIL+DTpGm4OXsbbST3Yr1WtR1NxsWgMHNqSu9es3hiuz/mKYx4nT55k1apVLF++nJ49e6asJrdixQrGjx9PlSpV+PHHH2nRogUffvghIsJTTz3F559/Tnx8PK1bt+att95KNYXFV199xeTJk1Omw16yZAlvvPEGl156KfHx8URGRtKwYUM++uijlDUmACZOnMi0adMICgri+uuvZ8KECQFZp8DkDStBZOKvHk03ogQxKvizVPuNyW9z586la9euXHrppVSqVIn169enHNuwYQMvv/wy27ZtY9euXSkDrEaMGMHatWv58ccfiY+P54svUi/s2KFDB3766Sc867y/++67DB48mAkTJhAWFsbGjRvPmxBu0aJFzJ07l9WrV7Np0yYeeeQRIDDrFJi8YSWITHiWJT2UUJkPkjpxV/Ai3pVeDO1yfX5nzeSFXPymHyjTp0/ngQceAGDgwIFMnz6d5s2bA9CyZUtq1aoFQGRkJHv27KFt27YsX76ciRMncvr0aY4dO0bDhg254YYbUq4pItx+++18+OGHDB48mO+//54PPvggw3wsXbqUwYMHp5QOPLObBmKdApM3AhogRGQ0cA+gOKvKDQaeBm4AzgG/AoNVNc5H2j3ACSAJSExvKHigeS9L+kZcT24OXsbIoJnc90kNJi3ewcNd6tvoapNvYmNjWbZsGT/++CMiQlJSEiKSsgqa91oEnnURzpw5w9///neio6O58MILGT9+vM+1AwYPHswNN9xAqVKl6N+/f7prMHv4WgcCnIbvuXPn0rRpU957771cmdra5I2AVTGJSE1gJBClqo2AYJy1pZcAjVS1CfAzMDaDy1yrqpH5FRw8PD2anhjQjv8mdaOrrKah7LIeTSbfzZo1izvuuIO9e/eyZ88e9u3bR926dVm5cmW6aTzBoEqVKpw8eTLdXksRERFERETwr3/9K9VqaaGhoT7XlejcuTNTp07l9OnTABw7dgwIzDoFJm8Eug0iBAgTkRCgNBCjql+6a04D/ADUCnAecs2kxTt4O7Ebf2hZHgn5BLAeTSZ/TZ8+nRtvvDHVvr59+563TrS38PBwhgwZQuPGjenduzdXXHFFuufeeuutXHjhhTRo0CBl39ChQ2nSpAm33nprqnO7du1Kz549iYqKIjIyMqULayDWKTB5I6CzuYrIKOAZIB74UlVvTXP8c+ATVf3QR9rdwB841VNvqeqUdO4xFBgKULt27RZ79+7N3YfwUnfMAhS4J3gBj4V+xC3n/sl3yY0QYPeE7gG7r8lbBXE21/wyYsQImjVrxt13353fWTG5oMCsByEiFYFeQF0gAigjIrd5HX8USATSK3O2UdXmwPXAcBHxuRCtqk5R1ShVjapatWquPkNanp5L05I6sV+rMCZkOkKy9WgyRVKLFi3YvHkzt912W+YnmyIpkI3UHYHdqnoEQERmA62BD0XkTqAHcJ2mU4RR1Rj3399FZA7QEvgmgPnNlKdHU3xCCV5M6MeLJd6kd+garukyLD+zZYqB/FgPwrOWtSm+AhkgfgOuFJHSOFVM1wHRItIV+Adwjaqe9pVQRMoAQap6wn3fGXgqgHn1i3ePpnlxbRmSvJDRQTPo+EmU9WgyAWXrQZj8ELAqJlVdDcwC1uN0cQ0CpgCTgXLAEhHZKCJvAohIhIgsdJNXA1aKyCZgDbBAVf8XqLxmhadH0wsDmjMp6WZqy+/cHPyV9WgyxhQ5AR0HoarjgHFpdl+czrkxQDf3/S6gaSDzllOTFu/gQFITvgtqwMiQ2cxOupoTCaVt1TljTJFhU21kkzMXk/BM4q1U5CTDQ+Z57TfGmMLPAkQ2eXoubdW6zEluy+DgRdSS361Hk8l37733HiNGjPB5LLtTc3uuGxMTc97+4cOHExkZSYMGDQgLC0uZmjyzacNNwWdzMWXTXz2akpiUcBPdSq5mTOgnJHaxiciKmwW7FvDK+lc4dOoQ1ctUZ1TzUXSvVzDHxeRkau733nuPRo0aERERkWr/66+/DsCePXvo0aOHNaYXIVaCyCbvVecOU5m3k7rTI+h7Ppg501adK0YW7FrA+O/Gc/DUQRTl4KmDjP9uPAt2LcjRdXv37k2LFi1o2LAhU6b8NUa0bNmy/OMf/6BFixZ07NiRNWvW0L59e+rVq8f8+fNTztu3bx9du3alfv36PPnkk6nSe0yaNIkrrriCJk2aMG6c01S4Z88eLr/8coYMGULDhg3p3Lkz8fHxzJo1i+joaG699VYiIyOJj/evKvWXX36hS5cutGjRgnbt2vHzzz8DcNtttzF8+HCuvfZaLrroIr755hvuvPNOLrvsspRBeYmJiYSHhzN69GiaN29Op06diI2Nzf4P1WSdqmb6AsoD9YHauKOvC+KrRYsWmh/mrN+vTR+drYefqK3rHm+uf/vH53rZY4t0zvr9+ZIfkzPbtm3z+9xOn3bSRu81Ou/V6dNOOcpDbGysqqqePn1aGzZsqEePHlVVVUAXLlyoqqq9e/fWTp066blz53Tjxo3atGlTVVV99913tXr16nr06NGU9GvXrlVV1TJlyqiq6uLFi3XIkCGanJysSUlJ2r17d/3666919+7dGhwcrBs2bFBV1f79++u0adNUVfWaa65JuY4vu3fv1oYNG6ba1759e925c6eqqq5cuVI7dXJ+Lrfeeqveeuutqqo6a9YsLV++vG7dulWTkpK0adOmumXLFk1ISFBAZ8yYoaqqjz/+uI4aNSonP9Ziz9ffNhCt6XymplvFJCLlgPuAW4CywFGgFFBZRFYC/1HVbwMavQqJSYt3EJdYgkl6E5NCp9Az6DvmJ7SxHk3FwKFTh7K031+vvvpqymI9+/bt45dffqFy5cqUKFGCrl27AtC4cWNKlixJaGgojRs3Zs+ePSnpO3XqlDKArk+fPqxcuZKoqL9mU/jyyy/58ssvadasGeAsOvTLL79Qu3Zt6tatS2RkJOCMpva+blbExcXxww8/0Ldv35R9iYmJKe8904s3btyYiIiIlPmeGjRowJ49e7jssssICQmhf//+gFPquOWWW7KVF5M9GbVBzMGZBuM6VU0p14kzn29L4HYRuURVpwY4jwWep+fSrKR23B68hLGh01lytgUx501iboqa6mWqc/DUQZ/7s2vFihUsXbqU77//ntKlS9O+ffuUGVhDQ0NTptQOCgpKmc47KCgo1Ydv2mm3026rKmPHjuXee+9NtX/Pnj3nTRHub3VSWqpKlSpV0m2T8M679z29nyWz5zCBlW4bhKp2VNV3vYODu19VdbWqjrDg4PD0XFKCeDLhDmrIMe4LmW89moqBUc1HUSq4VKp9pYJLMar5qGxf8/jx41SsWJHSpUuzfft2fvjhhyxfY8mSJRw7doz4+Hjmzp1LmzZtUh3v0qULU6dOTVkq9MCBA/z+++8ZXrNcuXKcOHHC7zxUrFiRGjVqpJSEkpOT2bRpU5aeIyEhgdmzZwPw8ccf07Zt2yylNzmTboAQkfu83tscvRl4uEt9wkKDAVin9Zmb1Jp7gxcgx3+zBusirnu97oxvPZ4aZWogCDXK1GB86/E56sXUtWtXEhMTadKkCY8//jhXXnlllq/Rtm1bbr/9diIjI+nbt2+q6iVw1m645ZZbuOqqq2jcuDH9+vXL9MN/0KBBDBs2LEuN1DNmzODNN9+kadOmNGzY8LylTTNToUIF1q9fT/PmzVm5ciWPPfZYltKbnEl3um8RWa/ObKqp3hdkUVFRGh0dnS/3nrvhgDO6Oi6e6sSyrOT/sSK5KX9PeICw0GCe69PY2iMKCZvuu2BITEykSpUqxMVZXW1uCdR031bxlwnPHE01w8M4RGXeSLyBbsFraBO0xRYVMsYUShkFiHARuUFEegHlRaSn9yuvMljYeBqspyT1YG/yBTwZ8j6hJNoUHMZkUUhIiJUe8llGAWIVcBPQD/gO6O/16hf4rBVOnobps5RgfOKdXBwUw+DgRdZgXcikV/VqTGGVnb/pdLu5qurtOcpNMeU9Bcfy5GYsSWrOqJDZzET6tlIAACAASURBVI9rTZsJy2zNiEKgVKlSxMbGUrlyZetWaYoEVSU2NpZSpUplfrKXjAbKDcRZL9pn2BGROkCEqmZ/cpciyHtRoQNx8TyVeDtLSzzCo6EfcX/cSMbO3pLqPFPw1KpVi/3793PkyJH8zooxuaZUqVLUqlUrS2kyGihXE9ggImuAdcARnJHUFwPtgT9xVoYzafRuVpPezWrSZsIy9sVV4z+JPRkd+hkzkq5lVUJjG2FdwIWGhlK3bt38zoYx+S6jgXIvAFE4I6ovBLrjrCkdC9ytqr1VNcOuOSIyWkS2isiPIjJdREqJSCURWSIiv7j/VkwnbVcR2SEiO0VkTLafMB95GqbfTLqB3cnVeDrkXUpyzhqsjTGFQobdXFU1UVUXqepjqnq3O3r6dVXdndmFRaQmMBKIUtVGQDAwEBgDfKWqlwBfudtp0wYDrwPXAw2Am0WkQVYfLr95N1g/nngX9YIO2QhrY0yhEejpvkOAMBEJAUoDMUAv4H33+PtAbx/pWgI7VXWXqp4DZrjpChXvEdYrkxszL6k19wXPp+TxX22EtTGmwAtYgFDVA8DzwG/AQeC4qn4JVFPVg+45B4ELfCSvCezz2t7v7juPiAwVkWgRiS5ojYrea0YA/CvhNs5SgqdD3uVA3GnGzt5iQcIYU2AFLEC4bQu9gLpABFBGRG7zN7mPfT57U6nqFFWNUtWoqlWrZi+zAeQ9wvoI4fw7cSBtgrfSN+hbG2FtjCnQMg0QIlJVRN4SkS/c7QYiMsiPa3cEdqvqEVVNAGbjNHIfFpEa7rVqAL6mkNyP0zDuUQuneqrQ8jRMf5zUgejkS3ks9EMqc9warI0xBZY/JYj3gK/56wP7F+AhP9L9BlwpIqXdNSSuA34C5gN3uufcCczzkXYtcImI1BWREjiN2/N9nFdoeE8JPibhHkpzhidCpxEkQt0xC6xNwhhT4PgTIC5Q1Y+BZAC3NJCUWSJVXQ3MAtYDW9x7TQEmAJ1E5Begk7uNiESIyEI3bSIwAliME1RmqurWrD1aweLdYL1Ta/F6Ym96BX/H1bIeBQ7ExVubhDGmQEl3uu+UE0RWAH2AparaXESuAF5U1avzIH9Zkp/TffvDMyV4TFw8JSSR+aH/pKzE0/nsRE7hlDBqhoexakyHfM6pMaa4yOl03/8HfA7UE5GvgenA/bmYv2LD02C9e0J3zmkIYxPuoQbH+EfIjJRzrE3CGFNQZDTVBgCqGi0i1wKX4/Qu2uaOTTA5EBEexvq4S5ma1JV7QhaxKLkl3yc3tEF0xpgCw59eTH2Akqq6CegKfCgikQHPWRHnaZN4PvEmdidX498hUwjjDAfi4q3B2hhTIPhTxTReVU+ISGvgBuAT4M3AZqvo8wyiqxweziMJ91JLjvJIyCeANVgbYwoGfwKEp8dSD+A/qvoZUDJwWSo+PG0SMRWa8X5SZwaHLKal/ARgg+iMMfnOnwBxUERexxmLsNAdlxDoOZyKlZi4eCYmDmBPcjWeD32TMsSn7DfGmPzizwf9TTgD5bqp6h9AFXzMwGqyLyI8jHhK8VDCMGrJUR4N+RBw5hax9ghjTH7JNECo6klgO9BBRO4DqqjqooDnrBjxNFiv0/pMSerBLSHLaR+0AbD2CGNM/vGnF9OjOGMfauLMifSxiIwNdMaKE+9ZX19M7Mf25Av5d+jbhHMCsPYIY0z+8KeK6TbgClV9VFUfxVmr4Y7AZqv48TRYJxDKgwn3UZETPBP6XzyT2Fp7hDEmr/kTIPaSekBdCLArMNkxEeFhbNM6vJjYn+7Ba+gX/A1g7RHGmLznT4A4DWwVkXdE5G2ciffiRORFEXkxsNkrfjztEVOSevBD8uWMD3mf2nIYsPYIY0ze8idALADGA98DPwBPAcuAre7L5CJPe0SN8DI8eO4+kgnipdD/EOwOR7H2CGNMXslwLiYRCQbaqeqdGZ1nclfvZjXp3awmdccs4NGEu3itxGRGhszhpcR+gLVHGGPyRoYlCFVNAmqISGge5cd4iQgP4/Pk1nyWdDUjgufQyh1lbe0Rxpi8kOlsrjgN0t+KyDzglGenqr6aUSIRqY8zb5NHPeAJ4CqgvrsvHIhT1fMm/xORPcAJnKk+EtObr7woe7hLfcbO3sLjCYNpVuIXXi7xOteffY44yqW0R4BT4jDGmNzmTxvEEWAJUBqo6vXKkKruUNVI98O/BU5j9xxVHeC1/zOctarTc617brELDvBXe0TF8IrcnzCSSvzJpNC38HR9tfYIY0wg+bMexOMAIhLmbmenAvw64FdV3evZ4a5TfRNgy6dl4K/2iHgmJN7MuNBpDE7+H+8mXQ9Ye4QxJnD8GUndQETWAr8AO0VktYhcnsX7DMQZje3tauCwqv6SThoFvhSRdSIyNIP8DRWRaBGJPnLkSBazVXhEhIfxblJXliS1YGzIxzSVnYC1RxhjAsefKqYpwD9VtZaq1gQeBd729wbu7K89gU/THLqZ84OGtzaq2hy4HhguIu18naSqU1Q1SlWjqlbNtOar0HLGR4Twfwn38jsVeb3Eq1TgJGDjI4wxgeFPgCinqks8G6q6FCiXhXtcD6xX1cOeHSISAvQhdSN2Kqoa4/77OzAHZ4qPYsvTHlE2vCrDz43kAv7g+dA3sfYIY0yg+BMg9ojIWBGp5b7G4Ey/4S9fJYWOwHZV3e8rgYiUEZFynvdAZ+DHLNyzSPLM17RZL+bZxFvpFLyeocFfpBy39ghjTG7yJ0DcBVwILHRftYDB/lxcREoDnTi/p9J5bRIiEiEiC93NasBKEdkErAEWqOr//LlncRARHsZ7SV1YkNSSR0I+4cqgbYC1Rxhjcpeoqu8DIiWBsqoam2Z/FeCEqp7Ng/xlSVRUlEZHR+d3NgJu7oYDjJ29haCEk8wr8TjhcpIeZ5/hEJUBCAsN5rk+jW18hDEmUyKyLr2hBBmVIF7BdxfUboBN0pePPO0R4eGVuDdhNKU4x39KvEIoiYC1RxhjckdGAaKdqqbteQQwDWgfmOwYf3naI3ZpTR5JGErzoJ08EfJByvEDcfFW3WSMyZGMAoT42qlOnZTPYybvRYSHsTD5St5M7MHtIUsZELw85Zh1fzXG5ERGAeKoiLRIu1NEmgPHApclkxWe9SMmJg7km6TGPB0ylebyc8pxq24yxmRXRgHiYeAzEXlMRK53X4/jzJ/0cN5kz2TGe/2IEQn3E6NVeLPEy1TziuFW3WSMyY50A4Sq/gBcCYQBw9xXGNBaVb/Pm+wZf3jaI8qFV2VIwkOU5gxTSrxISc6lnGPVTcaYrMpsPYhDqvqoqvZyX/9U1YN5lTmTNQ93qc/+kL8xOuHvNJbdqUZag1U3GWOyxp+BcqaQ8FQ3bSt/NRMTB3BD8A+MDJ6T6hyrbjLG+MufBYNMIeKZHrzNc8rFp2N4MHQWv2oEC5KvTDnHFhsyxvgj3RKEiLzn/jsiz3Jjcs3DXS/jKYayNvlSXgh9g2aSelZ1q24yxmQmoyqmliJSExgiIuVEpLz3K68yaLKnd7OaPNWnOePCHuWQVuLtEi9woRxOdY5VNxljMpJRgHgHWAFcBmxN8yr2M6sWBr2b1WTh2N6MKfU4wSTzXujElDUkPKx3kzEmPRl1c31RVS8BPlDV2qp6oderdh7m0eTQwOs7cH/y/1FLjpzX/RWsuskY41umvZhUdYiINBKRYe6rQV5kzOSe3s1q0q/PTTxbYiStgrbzQugbCMmpzrHqJmNMWv6sST0cmAnUdl+fisjfA50xk7t6N6vJk4+OY3LInfQIXs0TIdPwHiMBVt1kjEnNn3EQ9wIt3UFy/wRa4YyqNoVQrW6P8H5yNwaHLE61Gp2HVTcZYzz8CRACJHhtJ+DHbK4iUl9ENnq9/hSRB0RkvIgc8NrfLZ30XUVkh4jsdJc5Nbmgd/NaVOg9kaVBbfln6HT6BX993jlW3WSMAf8Gyk0DfhCRz9ztG4H3M0ukqjuASAARCQYOAHNwlit9SVWfTy+te/7rOMuV7gfWish8Vd3mR35NJno3vxCazGbNs52YEPI2x7UMS5JTLyhlg+mMMf40Uk8EhgKngXhgWEYf7um4DvhVVff6eX5LYKeq7lLVc8AMoFcW72kyElKSw9e/w1bqMTn0tZR1rb3FJyTx0MxN1B2zwEoUxhRDfs3FpKpr3W6vL6jq2mzcZyAw3Wt7hIhsFpGpIlLRx/k1gX1e2/vdfecRkaEiEi0i0UeOHMlG1oqvG664lP3dPiAmqDpvh75AE/n1vHOSVFGsAduY4ijgk/WJSAmgJ+BZvvQN4CKc6qeDwAu+kvnYpz72oapTVDVKVaOqVq2aCzkuXrq3akjd0V/yZ1B5Pigxgcvkt3TPtQZsY4qXvJjN9XpgvaoeBlDVw6qapKrJwNs41Ulp7Qcu9NquBcQEPKfFVfkItl43jTOUZFqJZ7lI0i8lWAO2McWHXwFCRGqJyLXu+5IiUiYL97gZr+olEanhdexGfE/bsRa4RETquiWQgcD8LNzTZFHntq348bppBIvwYYnnqCu/p3uuVTcZUzz4M1DuLpwP53fcXX8D5vlzcREpjdMTabbX7okiskVENgPXAqPdcyNEZCGAqiYCI4DFwE/ATFXd6tcTmWzr2K4tlYYtpEZpZUH4BC4OjU333PiEJB74ZKOVJowpwkTVZ9X+XyeIbMSpBlqtqs3cfZtVtUke5C9LoqKiNDo6Or+zUfgd3ATv9+RUUBluT3qC9cfLZXh6WGgwz/VpbN1hjSmERGSdqkb5OuZPFdMZt6up52LB+DFQzhRiNZrCHXMpk3yS2WHP0rzCyQxPt8ZrY4omfwLEKhF5BCjltkN8Apw/R4MpWiKawe1z4Ewc04KfzLC6Cazx2piiyJ8A8QhwAtgOjAK+Ah4NZKZMAVGzBdwxnzLJp5hf9hlalo/L8HRrvDamaMmwDcKtTpqqqnfmXZayz9ogAuTQFni/J4SUZOkVb3P/klPEJyRlmKRmeBgPd6lv7RLGFHDZboNQ1SSghoiEBiRnpnCo3hgGLYDkRDr+MIjXrwulZnhYhkmsNGFM4edPL6Y3cUY9zwNOefar6quBzVrWWQkiwI7uhA96wbkTcOtntPnoBAfi4jNNZqUJYwqunPZiOgIsAUoDVb1epripcjHctQjCKsEHvfh382OEhQZnmsxKE8YUTpmWIAoTK0HkkROHYNqNELuTNc3/zegtdfwqSYCVJowpaHJUghCRJSLyZdpX7mfTFBrlqsPghRDRnJZrH2TVdbt5eUCklSaMKWL8aYNo5bVZCugLnFXVhwOZseywEkQeO3caPh0EvyyG9mOZW/42Jn35s5UmjClEMipBZKuKSUS+VtVrcpyzXGYBIh8kJcD8kbDpY2h+B3R/ibmbDzN29pZMu8KCTdNhTH7LKEBkuuSoiJT32gwCWgA10jndFDfBodD7P1ChJnwzCU4cone/d4HGTFq8I9PShGeaDgsQxhQ8/vRi2oozJfdWYAPOKOohgcyUKWREoMNj0OMl2LkU3utO74uCWDWmg19tEzZNhzEFkz9tEKGqmpBmX4g7JXeBYlVMBcDPi+HTwRAWDrfMhOqNmLvhgF+lCcFZNtDaJozJOzkdB7Hax741OcuSKbIu7QJ3/Q80GaZ2gZ+/pHezmn6VJjxfVaynkzEFQ7oBQkQuEJGmQJiINBaRJu6rLc6gOWN8q9EEhiyDSnVh+gD4bjKo0rtZTZ7r0zjTaTrAFiQypiBIt4pJRAYDd+FMs7HR69AJ4F1V/TTDC4vUx5ka3KMe8ARQE7gBOAf8CgxW1fOmCRWRPe69koDE9IpA3qyKqYA5dwrmDIOf5kPkbdDjRQgpCUCbCcv87g5rPZ2MCZwcdXMVkZtUdWYOMxAMHABaAfWBZaqaKCL/BlDVf/hIsweIUtWj/t7HAkQBlJwMX0+Ar/8NF7aCm6ZBuWrM3XDA766wHtY2YUzuy1E3V1WdKSJdgIY4A+U8+5/NQh6uA35V1b3AXq/9PwD9snAdU9gEBcG1/4Sql8G84TDlGhjwIb2bOX+PnsZrTwN1RjxtE4AFCWPygD8liP8A4UA74F2ckdQ/qOpdft9EZCqwXlUnp9n/OfCJqn7oI81u4A+cz423VHVKOtceCgwFqF27dou9e/f6Os0UBId+hBm3wImD0P0FZ2Cdy9+eTh5WmjAmd+S0immzqjYRkU2q2lREygGfqWpnP29eAogBGqrqYa/9jwJRQB/1kQkRiVDVGBG5AGc22ftV9ZuM7mVVTIXA6WMw6y7YtRya3Q7dnofQlIJplqqerG3CmJzLaTfXM55/RaS6u10nC/e/Hqf04B0c7gR6ALf6Cg4Aqhrj/vs7MAdomYV7moKqdCW47TO4+v9gwzSY2hn+2JNy2Ho6GVNw+FOCGA+8DHQCXsPpVfS+qv7TrxuIzAAWq+q77nZX4EXgGlU9kk6aMkCQqp5w3y8BnlLV/2V0LytBFDI7/gdzhjrve/0HLu+R6nBWShOhQULZUiHEnU4gwqqfjPFbtquYRCQIuEJVV7vbYUCYqh7z88algX1APVU97u7bCZQEYt3TflDVYSISAbyjqt1EpB5OqQGchvSPVfWZzO5nAaIQOrYbZg2GmA1w5d+h45MQUiLlcFbbJjys+skY/+S0DeIHVb0yIDnLZRYgCqnEs7DkCVj9JkQ0h37/hUr1Up2SnW6xYI3ZxmQmpwHiaSBaVecFInO5yQJEIbdtPswf4Yyd6PESNOmf6nB2SxM2x5Mx6ctpgPgDqACcBeJx/7+paqXczmhOWYAoAuL2wewh8Nv30PRmuH4ilCqf6pTslibAqp6MSSunvZiqAKFAWaCqu10197JnjJfwC+HOL+CaMbD5E3izDez9PtUp3j2dBAgPCyU0WPy6vPV8MsZ/fq0oJyIDcRqanxWRWkA1VV0X8NxlkZUgipjfVju9nOJ+gzYPQPuxqRqwvWWn+smqnozJYQlCRCYD1wK3u7tOA2/mXvaMSUftVjBsJUTeAitfhLc7OKOxffB3SnFvNr24MRnzp4qptareiztgzu3i6vtrnDG5rWQ56PU6DJwOJw/DlPbwzfOQ5Hu9qrQD7fyreLKqJ2N88aeRejVwFU5PpuYiUhlYqqrN8iKDWWFVTEXcqVhY+BBsnQM1Ip3AUb1Rhkms6smYjOW0F9MdwI048yZNBW4CnlTVGbmd0ZyyAFFMbJ0LC/8P4v9wpuy4+sGUdSbSk92eTxYsTFGXowDhXqAh0NHd/EpVfVcE5zMLEMXI6WOw6B+wZaYzlfgNr0DtjMdzepcm/JlePC3rImuKotwIEE2Atjj/p1ap6ubczWLusABRDP38JSx4EI7vg6i74bonICw802TZHXQHVpowRUtOq5geBW7BmRtJgF7AR6r6XG5nNKcsQBRTZ0/Ciufgh/9A6SrQ9Tlo1Bck8yZqq3oyxV1OA8RPQAtVPe1ulwbWqerluZ7THLIAUczFbIQvRkPMeqh7jbPWRNVLM02W06onCxamMMtpgPgfcJOq/ululwemq2r3XM9pDlmAMCQnwbp3YelTkHAarrwPrnnE6S7rBwsWprjJaYCYDVwBLMb52+8MrAQOA6jqg7ma2xywAGFSnDwCX42HDR9C2erQ6Slo3N9ZI9tPOWmnAAsWpnDIaYC4O6PjqvrfHOQtV1mAMOfZH+10iY3ZADWj4Pp/Qy2f/xfSlZPJAT1sQSNTUOW4F1M2b1of+MRrVz3gCeADd38dYA9O9dUfPtJ3BV4BgnEWEpqQ2T0tQBifkpNh8wxYOt4Zjd24v9PbKby235fIadVTWla6MAVFTksQXYGngb/hrO6W5em+RSQYOAC0AoYDx1R1goiMASqq6j98nP8zzjKn+4G1wM2qui2j+1iAMBk6ewJWvgTfvw6qTvvE1Q9CqQpZuowFC1OU5DRA7MQZPb0FSPbsV1W/y9si0hkYp6ptRGQH0F5VD4pIDWCFqtZPc/5VwHhV7eJuj3XvmWHXWgsQxi/H98Oyf8Gm6RBW0RmNfcU9EFoqy5eyYGEKu5wGiBVAB1VNzvDEjK8xFVivqpNFJE5Vw72O/aGqFdOc3w/oqqr3uNu3A61UdYSPaw8FhgLUrl27xd69e7ObTVPcHNzkVDv9ugwqXOhMJ95kAASHZOtyFixMYZTTANESGAeswFlVDgBVfdXPm5cAYoCGqnrYzwDRH+iSJkC0VNX7M7qXlSBMtuxa4QSKmA1Q+RK49p/QoHeWejyl5QkWMXHxVAgL5dS5RBKSsh8yLFiYQMlpgFgEJHB+FdPjft68FzBcVTu721bFZAoeVfjpc1j+DBzZDtUaQ/t/wGU9/BqRnZncLF1YsDC5KacBYp2qtsjBzWcAi1X1XXd7EhDr1UhdSVUfSZMmBKeR+jqcxu21wC2qujWje1mAMDmWnARbZsHXE+DYLqjeGK75B9TvnqMShTcLFqYgyWmAmAj8T1WXZePGpYF9OMuVHnf3VQZmArWB34D+qnpMRCJwurN2c8/rBryM0811qqo+k9n9LECYXJOUCFs+hW8mOoHigobQ7iG36sm/Fev8EYhgER4Wigg25sL4JacB4g+gAs5So+fIRjfXvGIBwuS6pET48TP49nk4+jNUvthZH7vJgHTXx86u3G7k9rBShslITgOEz69LWenmmlcsQJiASU6Cn+bDty/AoS1QLgKuGg4t7vR7nqessGBh8kpurAcxEKea6FkRqQVUU9V1uZzPHLMAYQJOFX79Cr59CfauhJIVIGowtBoG5WsE5JaBDhZWJVW85bQEMRkIBdqp6uUiUgmn0fmK3M9qzliAMHlqfzR895pTspBgaNQHrvw7REQG7JaBChbebN6o4iWnAWK9qjYXkQ2q2szdt0lVmwYgrzliAcLki2O7YfWbzsyx507C39pAq3udnk/ZHHTnj7wIFmAljaIupwFiNXAVEO0GisrAUk+wKEgsQJh8deY4rJ8Ga96CuN+gfC244i5odgeUrRrQW6cdmCcCf5xOsMBhMpWtACEiIaqaKCJ3ADcCUcBUnHmZnlTVGYHKcHZZgDAFQnIS/LzYCRS7VkBQKDToBVfcDbWvypWBd/7Kq1KGN2sIL1yyGyDWq2pz931DoCPO736pqv4YqMzmhAUIU+Ac+Rmip8LGj+Hscah6GTS/E5oOhNJ521M8P4OFlTIKruwGiA0FsRopIxYgTIF17pQznmLd+3AgGoJLwuU9oNntzvrZuTRK21++qqTiTifkyrxR/rDAUXBkN0DsB15M76Kqmu6x/GIBwhQKh7Y4gWLLTKfdIrw2NL3FKVVUqpvfucuXkoaHr8BRwYJIQGU3QBwE3sD5nZ1HVZ/MtRzmEgsQplBJOAPbv4AN02DX14BC7dZOoGjQC8LCM71EoOVH43dmvLvhWvDIuRy3QRQWFiBMoXV8P2ya4SxiFLvTqYKq39VZHvXiTtlazCiQCmLgACuBZIe1QRhTWKhCzHrYPNOZVfb0UWe09uU9nIF4da+B4ND8zmW68rN6yl8WRFLLboCopKrHApqzXGYBwhQpSYmwewVs+cxZq+LcCQirBJffAA1vhDptC0WwKGilDH8UpyCS47mYCgsLEKbISjjjzAG1dQ7sWOSM2A6rCJd1h8t7Qr32EFIyv3Ppl8IcOLz5E0Suvawqy7cfISYuvsAGFQsQxhQlCfGw8ytnDqgdi+Dsn1CiLFzSyVkB7+KOBaKBO6sy6nqb191wAyWzoJIfpRQLEMYUVYlnYfc3Tm+o7Qvh1O8QFAJ/aw31u8GlXaBSvfzOZa4qKiUQf+RFQMm3ACEi4cA7QCOc57wLeADwrEEdDsSp6nnTX4rIHuAEkAQkpvcA3ixAmGItOdkZhLdjoRMsju5w9le+xAkUl3RypvooJFVRWZVZCaSoBhFvYaHBPNencZaCRH4GiPeBb1X1HREpAZRW1Tiv4y8Ax1X1KR9p9wBRqnrU3/tZgDDGy7Fd8POX8Mti2LMSks5BaBmodw1c1MF5VaqXp3ND5bfiEERqhoexakwHv8/PlwAhIuWBTTgLDZ13ExERnDWpO6jqLz6O78EChDG549wp2P0t7FwCvyyBuL3O/vC/wUXXOo3cda/J8/mhCqLCHkQE2D2hu//n51OAiASmANuApsA6YJSqnnKPtwNeTDdjIruBP3B+D2+p6pR0zhsKDAWoXbt2i7179+b2oxhT9MT+Cr8uc167v3W60CJQownUbecEi9pXBmQ51aLAnwZ1715MeTkupLCUIKKAH4A2qrpaRF4B/lTVx93jbwA7VfWFdNJHqGqMiFwALAHuV9VvMrqnlSCMyYakRGdw3q4VToP3vtVOdZQEO6vj1WnrLIJ0YatC2TuqoPC3l1ZOSim53QYRuOWuYD+wX1VXu9uzgDFuhkKAPkCL9BKraoz77+8iMgdoCWQYIAws2LWAV9a/wqFThyhfojwiwvGzx7P8vnqZ6rSr1Y5v9n+T42sF+n1Bz2tBz1+qvJatTrsOw/lm3zIOnf2D8noU2TeL4wdmUz4pGQkO4bgo5YPDkOBQjp87YT/LLOQ1qI5Q7uxxyrrnJWXwvlRQWeITktGgU0hyaQA06HS674OSKtK37pBc7RYb6Ebqb4F7VHWHiIwHyqjqwyLSFRirqtekk64MEKSqJ9z3S4CnVPV/Gd2vuJcgFuxawPjvxnMm6Ux+Z8UYkw9KBZdifOvxdK+XO20QgZ6E/n7gIxHZDEQCz7r7BwLTvU8UkQgRWehuVgNWisgmYA2wILPgYOCV9a9YcDCmGDuTdIZX1r+Sa9cLZBUTqroRZ6nStPsH+dgXA3Rz3+/Cadg2WXDo1KH8zoIxJp/l5udA3i5jZQKqepnq+Z0FY0w+y83PAQsQRcio5qMoFVyw1g0wxuSdUsGlGNV8VK5dzwJEEdK9XnfGtx5PjTI1EIQKJSoQXjI8W+9rlKnBgPoDcuVagX5f0PNa0POXr3kNdkb+9wAACzdJREFULU94aFkEqEAI4RqEqFIhMYnwpCTnfTKEE4wANULLM6DWddQoXb3A/MwKys+yRpkaWW6gzoxN1meMKVhUnRX2Dm+Fwz/Cke3w+09wZAckJ7gnibOWd9X6UOXSv/6tcqmNBs+i/BoHYYwxWScC4Rc6r/pd/9qflODML3VkO/y+3fn36M/Oet5JZ/86L6wSVLnEmaSw8kVQ+WLn34p1ocT/t3fnMVZWZxzHvz9BBQYGRHYQQaoooYqEuFRrF2u1akT/aKK1jdGmrUlVtDUWbdKYpjGm2sWmhsaOW6LRGpfU1rZqbN1qq1A3UBQRFAZBBpBtUBB5+sd5Ry6Xd5CZufe+c2d+n+TmXXiX54SZ+8x7znvOGVD78tQxJwgzqw999k1PCsMnw5SZO/fv+ATWL4M1b6URbNcuhjWL07hTL9+16zUGjcmSxQQYOjEljaET03b/A2pZmrrgBFFDpb2cq9UjdFTDKGZNn1XRekizbm2fPulLfuhEOOzru/7bRxvTU8faxbBuaVpf9zYsejTNnVGq3+A0eOEBB6flkINTNdaQ8elppheOS+U2iBqpZS/nzvSmNOt1tm6GD94p+SyFD95NI92uXwbby35X+w1JiWLweBg8FgaPg8aS5aBR3XqO8Pa4DaIbqGUv57belE4QZnuw/0AYNTV9yu3YAa0tsGF5ljCWZ+vLUyJ551nYuqHsJMHAkdA4Jn0GjYbG0alaq3F02h40CvZvrJs5OJwgaqTWvZzdq9qsC/bZBwaNTJ9x7Uxm+dFG2LgCNqyAjc1puek92PheqtJa+kxOEgH69k/XHTgqW46EgSPSsmFEtj4CGoYXPvufE0SNjGoYxcrWlTW9n5lVUb/G9BlxRPvHbGuFTatg08qy5SrY/H56fXfJk/BRTiKB1C7SMDz7DIMBw0rWD9y5HDAsPaVUmBNEjcyaPqumbRCV7E1pZp20X0P2qu2kPR/38YepSmvz6vRpXQ2bW9KytQVa10DLItjyHGxZx24zRQw4EK5aUvHwnSBqpK09wG8xmdlu9u2/842pz7LjE/jwg5Q0tqyFLWtSH5Eq8FtMZma9WJHzQZiZWZ1ygjAzs1xVbYOQNARoAqaSWlUuAk4Fvge0ZIddExF/yzn3NOAmoA/QFBHXVyPGSs3h7PYBM+tpqj0n9Z3AMxHRJGk/YABwObA5Im7cw3l9gEXAKUAzMBc4LyJe39P9OtoGUcQczu7lbGbdSSFtEJIagZOAWwEiYltErN/L048BFkfEkojYBtwLzPyMczqsiDmcKz1nrJlZtVSzDeIQUjXS7ZJektQkqSH7t0skvSrpNkl5QyiOBZaXbDdn+3Yj6fuS5kma19LSkndIu4rqbexezmZWD6qZIPoC04E5EXE00ArMBuYAk4BpwErgVznn5g1UklsXFhG3RMSMiJgxfPjwDgVYVG9j93I2s3pQzQTRDDRHxPPZ9v3A9Ih4PyI+iYgdwB9J1Ul55x5Usj0OeK/SARYxh7N7OZtZvajaW0wRsUrSckmTI+JN4GTgdUmjI6JtUKJzgAU5p88FDpU0EVgBnAt8q9Ixlvdu9ltMZmY7VXuojUuBu7M3mJYAFwK/kzSNVGX0DvADAEljSK+znh4R2yVdAjxKes31toh4rRoBnnHIGf7CNjPL4aE2zMx6MQ+1YWZmHeYEYWZmuZwgzMwslxOEmZnl6lGN1JJagHc7efowYE0FwylSTylLTykHuCzdUU8pB3StLAdHRG4v4x6VILpC0rz2WvLrTU8pS08pB7gs3VFPKQdUryyuYjIzs1xOEGZmlssJYqdbig6ggnpKWXpKOcBl6Y56SjmgSmVxG4SZmeXyE4SZmeVygjAzs1y9PkFIOk3Sm5IWS5pddDydJekgSf+StFDSa5LqetIJSX2ymQj/WnQsXSVpiKT7Jb2R/f8cX3RMnSHpiuxna4GkeyTVdjKVLshmr1wtaUHJvqGSHpf0VrbMm92y22mnLDdkP1+vSnpI0pBK3KtXJwhJfYCbgW8AU4DzJE0pNqpO2w78OCKOAI4DfljHZQGYBSwsOogKuQn4R0QcDhxFHZZL0ljgMmBGREwlDcN/brFRdcgdwGll+2YDT0TEocAT2XY9uIPdy/I4MDUijgQWAVdX4ka9OkGQZrNbHBFLImIbcC8ws+CYOiUiVkbEi9n6JtKXUO483t2dpHHAGUBT0bF0laRG4CTgVoCI2BYR64uNqtP6Av0l9QUGUIVZHqslIp4G1pXtngncma3fCZxd06A6Ka8sEfFYRGzPNv9LmoWzy3p7ghgLLC/ZbqZOv1RLSZoAHA08v+cju63fAlcBO4oOpAIOAVqA27MqsyZJDUUH1VERsQK4EVhGmkt+Q0Q8VmxUXTaybXbLbDmi4Hgq5SLg75W4UG9PEMrZV9fv/UoaCDwAXB4RG4uOp6MknQmsjoj/FR1LhfQFpgNzIuJooJX6qcr4VFY/PxOYCIwBGiR9u9iorJykn5Kqm++uxPV6e4JoBg4q2R5HHT02l5O0Lyk53B0RDxYdTyedAJwl6R1Sld9XJd1VbEhd0gw0R0Tb09z9pIRRb74GLI2Iloj4GHgQ+ELBMXXV+5JGA2TL1QXH0yWSLgDOBM6PCnVw6+0JYi5wqKSJ2bzZ5wIPFxxTp0gSqZ57YUT8uuh4Oisiro6IcRExgfT/8c+IqNu/VCNiFbBc0uRs18nA6wWG1FnLgOMkDch+1k6mDhvbyzwMXJCtXwD8ucBYukTSacBPgLMiYkulrturE0TWqHMJ8Cjph/2+iHit2Kg67QTgO6S/uF/OPqcXHZQBcClwt6RXgWnAdQXH02HZE9D9wIvAfNJ3R90MVSHpHuA/wGRJzZK+C1wPnCLpLeCUbLvba6csvwcGAY9nv/t/qMi9PNSGmZnl6dVPEGZm1j4nCDMzy+UEYWZmuZwgzMwslxOEmZnlcoKwXkHSgSWv/66StKJk+7kq3fNoSU3Z+rWSrqzGfdq59+cl3VGr+1nP1LfoAMxqISLWkvogIOlaYHNE3Fjl214D/KKaN5DUt2SQtk9FxHxJ4ySNj4hl1YzBei4/QVivJ2lztvyypKck3SdpkaTrJZ0v6QVJ8yVNyo4bLukBSXOzzwk51xwEHBkRr5TsniLpSUlLJF1WcuyPsjkWFki6PNs3oWy8/yuzxEZ2jeskPQXMkvTN7NxXJD1dcr+/UF9Dcls34ycIs10dBRxBGk55CdAUEcdkEzBdClxOmt/hNxHxrKTxpJ74R5RdZwawoGzf4cBXSD1e35Q0BzgSuBA4ljR45PPZF/8HnxHnkIj4EoCk+cCpEbGibKKYeaSBAX+516U3K+EEYbaruW1DQEt6G2gb0no+6csd0sB1U9KQRAA0ShqUzcPRZjRpmO9Sj0TEVmCrpNXASOBE4KGIaM3u+SDwRT57TLA/laz/G7hD0n2kQfTarCaNvGrWKU4QZrvaWrK+o2R7Bzt/X/YBjo+ID/dwnQ+B8ik5S6/9SXa9vCHnIQ3ZXFoFXH6t1raViLhY0rGkSZZeljQta3Ppl8Vh1ilugzDruMdIgzwCIGlazjELgc/txbWeBs7ORkltAM4BngHeB0Zkb1/tTxrGOZekSRHxfET8DFjDziHsD2P3ai6zveYnCLOOuwy4ORudtS/pS/7i0gMi4g1Jg3Oqnig77sXsddQXsl1NEfESgKSfk2YFXAq8sYd4bpB0KOlp5AmgrWH8K8AjHS2cWRuP5mpWJZKuADZFRM3n1s6eOp4CTsx7DdZsb7iKyax65rBru0MtjQdmOzlYV/gJwszMcvkJwszMcjlBmJlZLicIMzPL5QRhZma5nCDMzCzX/wEfTsjpGes+uQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "varyT_a(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 88.3\n",
      "0.6 92.9695\n",
      "1.2 98.73946749999999\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "98.73946749999999"
      ]
     },
     "execution_count": 236,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K=.275\n",
    "dt=.6\n",
    "T=85\n",
    "for i in range(0,3):\n",
    "    T= T+K*(T-current_temp(-dt*i))*dt\n",
    "    print(i*dt,T)\n",
    "T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It makes sense that the new euler approximation with varying T_a settles at a higher value than the analytical solution because it cant go below the highest value of T_a(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corpse was at 98.6 degrees at approximately 1.2 hours before initial time of discovery (9:48 am). The meaning of this is that the person was alive before this time."
   ]
  }
 ],
 "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
}