Final Intial Commit of 01_Getting-started-project.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computational Mechanics Project #01 - Heat Transfer in Forensic Science\n",
    "\n",
    "We can use our current skillset for a macabre application. We can predict the time of death based upon the current temperature and change in temperature of a corpse. \n",
    "\n",
    "Forensic scientists use Newton's law of cooling to determine the time elapsed since the loss of life, \n",
    "\n",
    "$\\frac{dT}{dt} = -K(T-T_a)$,\n",
    "\n",
    "where $T$ is the current temperature, $T_a$ is the ambient temperature, $t$ is the elapsed time in hours, and $K$ is an empirical constant. \n",
    "\n",
    "Suppose the temperature of the corpse is 85$^o$F at 11:00 am. Then, 2 hours later the temperature is 74$^{o}$F. \n",
    "\n",
    "Assume ambient temperature is a constant 65$^{o}$F.\n",
    "\n",
    "1. Use Python to calculate $K$ using a finite difference approximation, $\\frac{dT}{dt} \\approx \\frac{T(t+\\Delta t)-T(t)}{\\Delta t}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The constant value K is 0.6111111111111112\n"
     ]
    }
   ],
   "source": [
    "T1 = 85\n",
    "T2 = 74\n",
    "delta_t = 2\n",
    "T_ambient = 65\n",
    "dT_dt = (T2-T1)/delta_t\n",
    "K=-dT_dt/(T2-T_ambient)\n",
    "print(\"The constant value K is\", K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Change your work from problem 1 to create a function that accepts the temperature at two times, ambient temperature, and the time elapsed to return $K$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6111111111111112"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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",
    "    dT_dt = (Temp_t2-Temp_t1)/delta_t\n",
    "    K = -dT_dt/(Temp_t2-Temp_ambient)\n",
    "    return K\n",
    "measure_K(85,74,65,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. A first-order thermal system has the following analytical solution, \n",
    "\n",
    "    $T(t) =T_a+(T(0)-T_a)e^{-Kt}$\n",
    "\n",
    "    where $T(0)$ is the temperature of the corpse at t=0 hours i.e. at the time of discovery and $T_a$ is a constant ambient temperature. \n",
    "\n",
    "    a. Show that an Euler integration converges to the analytical solution as the time step is decreased. Use the constant $K$ derived above and the initial temperature, T(0) = 85$^o$F. \n",
    "\n",
    "    b. What is the final temperature as t$\\rightarrow\\infty$?\n",
    "    \n",
    "    c. At what time was the corpse 98.6$^{o}$F? i.e. what was the time of death?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[85.         83.80276397 82.67719664 81.6190078  80.62416405 79.68887344\n",
      " 78.80957101 77.98290521 77.20572511 76.47506841 75.78815014 75.14235204\n",
      " 74.53521257 73.96441757 73.42779138 72.92328861 72.44898627 72.00307654\n",
      " 71.58385976 71.18973805 70.81920918 70.47086083 70.14336525 69.83547413\n",
      " 69.54601394 69.27388136 69.01803912 68.77751206 68.55138338 68.33879117\n",
      " 68.13892512 67.9510234  67.77436983 67.60829105 67.45215405 67.30536369\n",
      " 67.16736046 67.03761836 66.91564286 66.80096902 66.69315977 66.59180418\n",
      " 66.49651591 66.40693177 66.3227103  66.24353048 66.16909051 66.09910664\n",
      " 66.03331214 65.97145621 65.91330309 65.85863112 65.80723192 65.75890956\n",
      " 65.71347987 65.67076968 65.6306162  65.59286637 65.55737632 65.52401077\n",
      " 65.49264254 65.46315207 65.43542696 65.40936151 65.3848564  65.3618182\n",
      " 65.34015911 65.31979657 65.30065297 65.28265535 65.26573509 65.24982771\n",
      " 65.23487257 65.22081267 65.20759443 65.19516745 65.18348438 65.17250067\n",
      " 65.16217447 65.15246642 65.1433395  65.13475894 65.12669203 65.11910801\n",
      " 65.11197799 65.10527479 65.09897285 65.09304816 65.08747813 65.08224153\n",
      " 65.0773184  65.07268998 65.06833863 65.06424776 65.06040177 65.05678601\n",
      " 65.0533867  65.05019087 65.04718636 65.0443617 ]\n",
      "[85.         83.81437419 82.69903381 81.64981225 80.6627899  79.7342795\n",
      " 78.8608124  78.03912556 77.26614937 76.53899621 75.85494962 75.2114542\n",
      " 74.60610602 74.03664366 73.50093977 72.99699309 72.52292102 72.07695255\n",
      " 71.65742167 71.26276113 70.89149657 70.54224105 70.21368985 69.90461559\n",
      " 69.61386365 69.34034785 69.08304643 68.84099817 68.61329884 68.39909783\n",
      " 68.19759492 68.00803737 67.82971703 67.66196776 67.50416287 67.35571287\n",
      " 67.21606317 67.08469208 66.96110885 66.84485178 66.73548659 66.63260471\n",
      " 66.53582179 66.44477629 66.35912809 66.27855722 66.2027627  66.13146138\n",
      " 66.06438689 66.00128866 65.94193098 65.88609209 65.83356341 65.78414869\n",
      " 65.73766335 65.69393371 65.65279643 65.61409781 65.5776933  65.5434469\n",
      " 65.51123066 65.48092425 65.45241444 65.42559473 65.40036492 65.37663077\n",
      " 65.35430362 65.33330004 65.31354158 65.29495443 65.27746916 65.26102043\n",
      " 65.2455468  65.23099047 65.21729705 65.2044154  65.1922974  65.18089776\n",
      " 65.17017391 65.16008578 65.15059569 65.14166818 65.13326991 65.12536949\n",
      " 65.11793743 65.11094595 65.10436893 65.0981818  65.09236146 65.08688615\n",
      " 65.08173543 65.07689005 65.07233191 65.06804398 65.06401024 65.06021563\n",
      " 65.05664597 65.05328793 65.05012895 65.04715724]\n"
     ]
    }
   ],
   "source": [
    "#a\n",
    "import numpy as np\n",
    "\n",
    "T_0=85\n",
    "T_ambient = 65\n",
    "K = 0.61111111111112\n",
    "t = np.linspace(0,10,100)\n",
    "Temp_analytical = T_ambient + (T_0-T_ambient)*np.exp(-K*t)\n",
    "delta_t = 0.1\n",
    "Temp_numerical=np.zeros(len(t));\n",
    "for i in range(1,len(t)):\n",
    "    Temp_numerical[0] = T_0\n",
    "    Temp_numerical[i]=T_ambient+(Temp_numerical[i-1]-T_ambient)*np.exp(-K*delta_t)\n",
    "print(Temp_analytical)\n",
    "print(Temp_numerical)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f9758427dd8>"
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeXhU5fn4//c9S0jCFpSIJKjgBrLIFhFBEUVksQJ1BbFFraJtwa1fFetSl58VRbtYsRbR4sdaN0TQiiCIqKCoICiLIoggCQhBSNhCMsv9++OciUOYJJOQyXq/rmuumfOc7TkJzJ1nF1XFGGOMKclT0xkwxhhTO1mAMMYYE5MFCGOMMTFZgDDGGBOTBQhjjDExWYAwxhgTkwUIU61EpK2IqIj4ajovJnFE5Dci8lZN58McHgsQplwislFECkRkr4jsEpG3ReSYas7Davf+e0UkJCIHorb/WJ15qWoiskRErqzpfFSEiHxXxu/jVlV9VlUvrOl8msNjAcLE60JVbQK0BrYB/6jOm6tqJ1Vt4ubhI2BcZFtV/1ydeamI6igp1URpTFVPiPp9fA5cG/X7+Et158ckhgUIUyGqegCYDnSMpIlIcxH5PxHJFZFNInK3iHjcfV4ReUxEdojIBuCCqPMuFZFl0dcXkT+IyMzK5E1ErheRtSKy0y3lZLrpyW611g3uX7673Ty2F5HPRCRfRF6MfNGKyGARWS8i97vX2iAil0bdJ0VE/iYim0XkRxH5h4g0KnHuPSKyDfiniKSLyDvuz2eniMwSkdbu8Y8DpwFT3b++HxeRDiISLPFsxaUM9zkWiMhkEdkFTCjr+WP8nBaKyLUl0taKyFD39/Wkm9d8EflSRNpX4ndxg4jMr8zP3z3nlyLylYjkichHItKx9LuZRLEAYSpERFKBy4ElUcn/AJoDxwNnA78Grnb3XQf8AugOZAGXRJ33JtBORE6JSrsSeKES+RoJ3AxcCLQClgP/KXHYAKCrm8c/ufm+FGgH9AIujjq2LZAEHA2MBZ4XkXbuvr8CbYAuQHvgZNwv6ahz/cAxwI04/8+eBo517xW5Bqr6Bw7+C/wPcT5yP2AF0BJ4PM7nj/gvMCqyISI9gSOAd3F+Vz2BE4AWwBXArjjzVJ64fv4i0ht4Cuff0JE4/x5m1kRJqcFTVXvZq8wXsBHYC+QBQWAL0MXd5wUKgY5Rx18PLHQ/LwBuiNp3PqCAz93+J/CQ+7kTzpdRo3LysxDnCzU67X1gdNS2HwjgfFkmu/fsGbV/NXBT1PZkYKL7eTBwAEiO2v8mcBvgA4qAzKh95wBfR527D/CXkf/ewNao7SXAlVHbHYBgiXOKjwFuAL6N9/lj3P8IoABo7W4/Djzlfh7q/mx6ARLnv4+D8h+Vx/nu54r+/P8N3FXiepuA02v6/0JDe1kJwsRrhKqmAY2AccAHInI0zl+wSTj/gSM2AZHqjQxgc4l90Z4HrhARAX4FvKqqhZXI33HA026VRB6QixPM2kQdsy3qc0GM7SZR27nqVKdF5zvDffmB1VH3mgkcFXXsj6oaiGyISFMReU5EfhCR3Th/qbesxDNG21xiO57nB0BVdwLzgMvcqsDLgRfd3e8AzwL/AraJyFMi0qTkNSop3p//ccAfI8/iPk86P/+bMtXEAoSpEFUNqeoMIAScCezA+Uv1uKjDjgVy3M9bcapaovdFX28Jzl/kZ+FUZ1S4esm1GbhKVdOiXimquqzcM2NrKSLJUdvH4pSctuJ88Z4QdZ/mqnpk1LElp0iegPNFfZqqNsMpRUkZx+8DvJF2DdfRJY4peU5Fn/8lnGqms93n+RhAHX9R1e7AqThVQjeVco1E2QzcW+JZUt1/d6YaWYAwFSKO4Tj101+ragh4FXjI/Uv5OOBWfq7/fhW4UUTaiEgLDq6rj/g/4EmcapVFlcza08DdkQZVEWkhIheXc05Z/MA9IpIkIucCA4HX3ZLBc8DfRaSl+/M4RkQGlnGtpsB+IE9EWgJ3l9i/Daf9JmILTglgtNto/DvK/+u5os8/C6dK7y7gJVWnHkdEeotIllvfvw8neIfKuXdVmwKMd/MhItJERIa57V+mGlmAMPF6S0T2AruBh4Axqrra3Tce58tkA7AIpxH0OXffM8Bc4EvgCyDWX4EvAJ2pfOkBVX0JJ8jMcKtxVuB8qVfWRpy/rH/EeZarVXWDu+9mnC/xpUA+MAc4sYxrPYZTpfQTzs9ndon9fwV+Lc4Yk0fdoHstTkPuDpwSWJkloYo+v6rux2lXGYDz+4pIA6bhtDdtwKlae6Kse1c1VV2M07j/Lzcf3+KULm3xmmom7h8OxtQYEUkBtgM9VHVdLcjPYOBJVS3rS9+Yes9KEKY2+C3weW0IDsaYn1m/YlOjRGQjToPtiBrOijGmBKtiMsYYE5NVMRljjImpXlUxtWzZUtu2bVvT2TDGmDpj2bJlO1Q1Pda+ehUg2rZty9KlS2s6G8YYU2eISMnZDYpZFZMxxpiYLEAYY4yJyQKEMcaYmOpVG4QxJnECgQDZ2dkcOHCg/INNrZOcnEybNm3w+/1xn2MBwhgTl+zsbJo2bUrbtm1xZmc3dYWq8tNPP5GdnU27du3KP8GV0ComEblFnMXmV4nIS+7Sg/eJSI6IrHBfQ0s5d7C7DOJ6EYk1A2iVmLk8h74TF9Buwtv0nbiAmctzyj/JmAbowIEDHHnkkRYc6iAR4cgjj6xw6S9hJQh3PdwbcVYaKxCRV4GR7u6/qupjZZzrxVlhaiCQDXwuIm+q6pqqzOPM5TncOWMlBQFnNuOcvALunLESgBHdbW0SY0qy4FB3VeZ3l+hGah+Q4s4tn4ozRXI8egHrVXWDqhYBLwPDqzpzk+auJRQ4wFjvW5zpcQJDQSDEpLlrq/pWxhhT5yQsQKhqDs48+D/grMKVr6rvurvHichX7jKMLWKcnsnBSypmU8qCKSIyVkSWisjS3NzcCuVxS14BAbyM9b3NL70fHZRujKldrrnmGo466ig6d+58UPrOnTsZOHAgJ510EgMHDmTXrl3F+x5++GFOPPFE2rdvz9y5c2Net3///rRv355u3brRrVs3LrnkkjLzsXHjxkPyUBWmTZvGli3x/g1dPRIWINwv/uFAO5x1fBuLyJU4i9SfAHTDCRyPxzo9RlrMWQVVdYqqZqlqVnp6zNHipcpIS0HxsCTckb6e1cW3yEhLqdB1jDGHqur2vauuuoo5c+Yckj5x4kQGDBjAunXrGDBgABMnTgRgzZo1vPzyy6xevZo5c+bwu9/9jlAo9uJ4L774IitWrGDFihVMnz79sPJZUmn3LKlBBQjgPOB7Vc11l2mcAfRR1W3uusZhnNXGesU4N5uD1zFuQ/zVU3G7bVB7UvxeFoc7cbTs4njZSorfy22D2lf1rYxpUCLtezl5BSg/t+8dTpDo168fRxxxxCHps2bNYsyYMQCMGTOGmTNnFqePHDmSRo0a0a5dO0488UQ+++yzuO931VVXHRQsmjRpcsgxoVCI2267jdNOO41TTz2Vf/3rXwAsXLiQc845hyuuuIIuXboccs5VV11F586d6dKlC3/961+ZPn06S5cuZfTo0XTr1o2CggKWLVvG2WefTc+ePRk0aBBbt24FnBLPzTffTJ8+fejcuXPxM33wwQfFpaDu3buzZ8+euJ+1NIns5voD0NtdR7YAZ2nDpSLSWlW3usf8ElgV49zPgZNEpB2Qg9O4fUVVZzDSEP3C7FwIwIBG39Bp+FBroDamHPe/tZo1W3aXun/5D3kUhcIHpRUEQtw+/Ste+uyHmOd0zGjGny7sVOG8bNu2jdatWwPQunVrtm/fDkBOTg69e/cuPq5Nmzbk5MQOUKNHjyYlxak5GDhwIJMmTYrr3s8++yzNmzfn888/p7CwkL59+3L++ecD8Nlnn7Fq1apDupWuWLGCnJwcVq1yvvry8vJIS0vjySef5LHHHiMrK4tAIMD48eOZNWsW6enpvPLKK9x1110895yzku++ffv4+OOP+fDDD7nmmmtYtWoVjz32GJMnT6Zv377s3buX5OTkeH+EpUpYgFDVT0VkOs46xEFgOc5i5FNFpBtOfc5G4HoAEckApqrqUFUNisg4nLWMvcBzUesfV6kR3TMZ0e1KfnzgfganrqWnBQdjDlvJ4FBeeiLEWuumtJ48L774IllZWRW+x7vvvstXX31VXNLIz89n3bp1JCUl0atXr5hjDo4//ng2bNjA+PHjueCCC4oDSrS1a9eyatUqBg50lhUPhULFQRBg1KhRgFOq2r17N3l5efTt25dbb72V0aNHc9FFF9GmTZsKP09JCR0op6p/wll4PdqvSjl2CzA0ans2hy7unhgiZDfP4sS8jwiFQni93mq5rTF1VXl/6feduICcGJ09MtNSeOX6M6o0L61atWLr1q20bt2arVu3ctRRRwFOiWHz5p/7umRnZ5ORkRH3dX0+H+GwE9BUlaKiokOOUVX+8Y9/MGjQoIPSFy5cSOPGjWNet0WLFnz55ZfMnTuXyZMn8+qrrxaXDKKv26lTJz755JOY1ygZ6ESECRMmcMEFFzB79mx69+7N/Pnz6dChQ9zPG4vNxeSS488mjb18tzL2L8QYE79I+160RLXvDRs2jOeffx6A559/nuHDhxenv/zyyxQWFvL999+zbt06evWK1eQZW9u2bVm2bBngtGcEAoFDjhk0aBD//Oc/i/d9++237Nu3r8zr7tixg3A4zMUXX8yDDz7IF198AUDTpk2L2w3at29Pbm5ucYAIBAKsXv1zJcorr7wCwKJFi2jevDnNmzfnu+++o0uXLtxxxx1kZWXxzTffxP2spbGpNlxts4bAFxPYsXI+J3c7s6azY0ydFmnHmzR3LVvyCshIS+G2Qe0Pq31v1KhRLFy4kB07dtCmTRvuv/9+fvOb3zBhwgQuu+wynn32WY499lhee+01ADp16sRll11Gx44d8fl8TJ48udTageg2iJYtWzJ//nyuu+46hg8fTq9evRgwYEDMEsG1117Lxo0b6dGjB6pKenp6cSN5aXJycrj66quLSycPP/ww4DSK33DDDaSkpPDJJ58wffp0brzxRvLz8wkGg9x888106uSU3Fq0aEGfPn3YvXt3cenjb3/7G++//z5er5eOHTsyZMiQSvyUD1av1qTOysrSw1kwaPMDHfkpKZNuE+ZVYa6MqR++/vprTjnllJrORoPXv3//4sbsior1OxSRZaoa82JWxRRlbWoPTiz4ipMmzLJ5mYwxDZ4FCNfM5TnMzDuRJnKALrKhSvptG2NMVVu4cGGlSg+VYQHCNWnuWhYFOxBWoa/H6Z9s8zIZYxoyCxCuLXkF5NGU1XocZ3pXHZRujDENkQUIV2T+pY/Cp9JD1tGE/QelG2NMQ2MBwhXpt/1h+FT8EuIMzxqbl8kY06BZgHCN6J7Jwxd1YUvTLuzTRpzjX8XDF3WxeZmMqWemTZvGuHHjyj0membVa6+9ljVrKr5e2cKFC/nFL35R4fNqCxsoF2VE90xGdM9k5aQenLXvK9p0i39YvjGm/pg2bRqdO3cunppj6tSpNZyjmmEliBgCbftzDD/y3dqVNZ0VY0yUESNG0LNnTzp16sSUKVMAZxruu+66i65du9K7d2+2bdsGwFtvvcXpp59O9+7dOe+884rTI/bs2UO7du2Kp8nYvXs3bdu25bXXXjtk6u3+/fsTGYQ7Z84cevToQdeuXRkwYADgzNzap08funfvTp8+fVi7tn70frQSRAzHnHYhrH6YrV/M5sQOp9Z0doypfd6ZAD9W8R9QR3eBIRPLPOS5557jiCOOoKCggNNOO42LL76Yffv20bt3bx566CFuv/12nnnmGe6++27OPPNMlixZgogwdepUHn30UR5//Of1yZo2bUr//v15++23GTFiBC+//DIXX3wxl156KZMnT445Wjk3N5frrruODz/8kHbt2rFz504AOnTowIcffojP52P+/Pn88Y9/5PXXX6/an08NsAARQ/pxHflRjiL5hw+ACTWdHWOM64knnuCNN94AYPPmzcVTa0fq+Xv27Mm8ec5UOdnZ2Vx++eVs3bqVoqKimFNvX3vttTz66KOMGDGCf//73zzzzDNl3n/JkiX069ev+FqRBYzy8/MZM2YM69atQ0RiTuxXF1mAiEWEbxqfRs89CzhpwiyOSmt62BONGVOvlPOXfiIsXLiQ+fPn88knn5Camkr//v05cOAAfr+/ePprr9dLMBgEYPz48dx6660MGzaMhQsXct999x1yzb59+7Jx40Y++OADQqFQuWtNq2rMNSXuuecezjnnHN544w02btxI//79D/t5awNrg4hh5vIcpuedTFMpoJust2k3jKkF8vPzadGiBampqXzzzTcsWbKk3OMzM50/6iLTgcfy61//mlGjRnH11VcXp0VPvR3tjDPO4IMPPuD7778HKK5iir7XtGnTKvRctVlCA4SI3CIiq0VklYi8JCLJIjJJRL4Rka9E5A0RSSvl3I0islJEVohI5adorYRJc9fyYbAjIRXO8n4F2LQbxtS0wYMHEwwGOfXUU7nnnnsOWk40lvvuu49LL72Us846i5YtW5Z63OjRo9m1a1fxKm3w89TbkUbqiPT0dKZMmcJFF11E165dufzyywG4/fbbufPOO+nbty+hUOgwn7T2SNh03yKSCSwCOqpqgYi8irNC3BZggbus6CMAqnpHjPM3AlmquiPeex7udN8R7Sa8jQLTk+4jiQDDih5y8gR8P/GCw76+MXVRfZ3ue/r06cyaNYsXXnihprOScBWd7jvRbRA+IEVEAkAqsEVV343avwS4JMF5qLCMtBRy8gp4P9SN2/yvks4ucmlh024YU8+MHz+ed955h9mzq2d147omYVVMqpoDPAb8AGwF8ksEB4BrgHdKuwTwrogsE5Gxpd1HRMaKyFIRWZqbm1sVWS+eduP9cDcA+nu/tGk3jKmH/vGPf7B+/XpOPvnkms5KrZSwACEiLYDhQDsgA2gsIldG7b8LCAIvlnKJvqraAxgC/F5E+sU6SFWnqGqWqmalp6dXSd4j027kN+vAVj2C87wrbNoNY3B68Zi6qTK/u0Q2Up8HfK+quaoaAGYAfQBEZAzwC2C0lpJrVd3ivm8H3gDiX228CozonsniOwewvVU/+shKzjkpZlu6MQ1GcnIyP/30kwWJOkhV+emnn0hOTq7QeYlsg/gB6C0iqUABMABYKiKDgTuAs1V1f6wTRaQx4FHVPe7n84EHEpjXUjXucgFNt8/kk0/mcMbAi2siC8bUCm3atCE7O5uqqso11Ss5OZk2bdpU6JyEBQhV/VREpgNf4FQlLQemAKuBRsA8d8DJElW9QUQygKmqOhRoBbzh7vcB/1XVOYnKa1nanTaUwvf8FH79DliAMA2Y3++PORrZ1F8J6+ZaE6qqm2tJy/98LmkHcji36C9kpKXYqGpjTL1RVjdXG0ldjpnLc3iroDPtPD9ynGy1UdXGmAbDAkQ5Js1dy7yg0931XM8KwEZVG2MaBgsQ5diSV8BmbcX6cAbner44KN0YY+ozCxDliIyenh/uyemeb2jGvoPSjTGmvio1QIjIb6M+d6ie7NQ+kVHV74Z64pcQ/T0rbFS1MaZBKKsEcV3U5/8mOiO1VWRU9fZmXcjV5gz2LbNR1caYBiHeKqZDV8hoQEZ0z2TRnefxY+tz6Scr6Nu2SU1nyRhjEq6sAJEmIheKyHCgmYgMi35VVwZrk7Qev6SJHGD14rdqOivGGJNwZY2kXgxc5n7+GLg0ap8CbyYqU7VVm+6D2D87Gf1mNvziyvJPMMaYOqzUAKGqv6rOjNQF4k9mVWovOu9ZxPET3qJ1WmMbVW2Mqbesm2sFzFyew8t7TiVd8m2tamNMvWcBogImzV3L/EBXAurlfO8ywEZVG2PqLwsQFbAlr4DdNOaTcEcGeT7DaYqxUdXGmPqp3AAhIheJSFP38wQReVVEuiU+a7VPZPT0nHAv2nm2cYr8cFC6McbUJ/GUIO5zF+7pA1wIvAI8ndhs1U6RUdVzQ1mEVBjq/dRGVRtj6q14AkTIff8F8JSqvo6z4E+DExlVnZx2NEvCHbnA8ykPDOtovZiMMfVSPAFiq4hMBi4HZotIUpznISK3iMhqEVklIi+JSLKIHCEi80RknfveopRzB4vIWhFZLyIT4n+kxBrRPZPFE86lVe+RHO/ZSuvCDTWdJWOMSYh4vugvAz4ALlDVXUBLoNwvbBHJBG4EslS1M+AFRrrnvqeqJwHvxbqWiHiBycAQoCMwSkQ6xvVE1aTtWZcTwsPe5dNrOivGGJMQ5QYIVd0L/AD0cpMKcdaVjocPSBERH5AKbAGGA8+7+58HRsQ4rxewXlU3qGoR8LJ7Xq3ha9aKb5K6cFLufNpN+B99Jy6w8RDGmHolnl5MdwN/Au52k5KJY3ZXVc0BHsMJLluBfFV9F2ilqlvdY7YCR8U4PRPYHLWd7abFyt9YEVkqIktzc3PLy1aVmbk8h1cLenKCZysny2YbNGeMqXfiqWK6BBgKzko57hd/s/JOctsWhgPtgAygsYjEO4FRrNljNdaBqjpFVbNUNSs9PT3Oyx++SXPX8nbgtOLeTGCD5owx9Us8AaJQVRX3C1pEUuO89nnA96qaq6oBYAbQB9gmIq3da7UGtsc4Nxs4Jmq7DU71VK2xJa+AHTTns/ApXOD5FBs0Z4ypb+IJEDPcXkzNReRq4F3guTjO+wHoLSKpIiLAAOBrnFlgx7jHjAFmxTj3c+AkEWnn9poaSS2bPTYyOO6t8Bmc6NlCR9l0ULoxxtR18TRSPwL8D+cLuivwkKr+LY7zPgWmA18AK917TQEmAgNFZB0w0N1GRDJEZLZ7bhAYB8zFCSqvqmq8DePVIjJobnaoFwH1Mty72AbNGWPqFXFqj8o5SKQNcJKqvi8iyYBXVfclPHcVlJWVpUuXLq22+81cnsOkuWu5f98DdPJs4r3B73HlGcdX2/2NMeZwicgyVc2KtS+eXkzX4JQeprpJxxK7WqjBiQyaO+X839BadpK+84uazpIxxlSZeNogbgR6A7sBVPVbYndNbbAyT7+IApKRla/VdFaMMabKxBMgDriD1YDiUc6xuqE2XEmNWdn0THoVfET7CTNt0Jwxpl6IJ0AsFpHbgWQROQdnNtf/JTZbdcvM5Tk8k9eTNNnHWZ6vbNCcMaZeiCdA3A7sAb4BbsKZP+muRGaqrpk0dy3vBzrxkzZlhHcxYIPmjDF1n6+snW510nOqOgb4Z/Vkqe7ZkleA4uPtUG8u8y6kKfvZQ6oNmjPG1GllliBUNQS0FhF/NeWnTooMjns9dBbJEiieesMGzRlj6rJ4qpg2AB+JyJ0icmPkleiM1SWRQXNf6gmsC2dyifcDkv0eGzRnjKnT4gkQucA8nOm606NexhVZaS4zLZXpoX6c5vmW6ztjK80ZY+q0MtsgAFT1nurISF03onsmI7pnEs7vQOivr3B8zps46x0ZY0zdVG6AEJEZMZLzgaXAM9FjJAx4mmewMiWL0/LncMKEtzg6rTG3DWpvpQljTJ0TTxVTNhAEXnBfRcBO4FTgmcRlrW6auTyHZ/eeQYbspLdntY2JMMbUWfEEiK6qepmqvqGqbwCjgNNU9XrgtMRmr+6ZNHct7wS6k6+pXOL9ELAxEcaYuimeANHKnc01IoOfG6kLqz5LdduWvAIKSeLNUB+GeD6jmbMQn42JMMbUOfGOpP5EROaJyHzgE+AOEWkMvJjQ3NVBkbEPr4T6kywBhrsjq21MhDGmromnF9ObIjIP6IgzSd9qVY38OfxYaeeJSHuceZsijgfuBc4AIgME0oA8Ve0W4/yNOFN8hIBgafOV1za3DWrPnTNWsipwPCvDbbnC+x6veQbZmAhjTJ0TTy+mFJw5mNqq6g0icqKInKSq75R1nqquBbq51/ACOcAb0avRicjjOD2iSnOOqu6I4zlqjUhvpUlz1/LSngH82f8sY4/dab2YjDF1TrkBAmf96ZXAme72FuA1oMwAUcIA4DtV3RRJcNepvgw4twLXqRMiYyIoPI2CiS9yypYZqI7GeWRjjKkb4mmDOElV/wwEAFR1PxVfD2Ik8FKJtLOAbaq6rpRzFHhXRJaJyNgK3q92aNSUlS0G0j/wEafeOd3WiTDG1CnxBIgidx1qBRCRdjhjIeIiIknAMJxSR7RRHBo0ovVV1R44w5F/LyL9Srn+WBFZKiJLc3Nz481WtZi5PIdHc3uTIkUM9y62MRHGmDolngDxADAHaCMizwPvA3dW4B5DgC9UdVskQUR8wEUc3Ih9EFXd4r5vB94AepVy3BRVzVLVrPT02jVF1KS5a1kaaMfKcFtGe98D1MZEGGPqjHIDhKrOAS4FrsP9olbV9ypwj1glhfOAb1Q1O9YJItJYRJpGPgPnA6sqcM9aITL24aXQAE7x/EAPWXdQujHG1GZlBggR8YrIEOBKnG6qB4Cf4r24iKQCA4GS8zkd0iYhIhkiMtvdbAUsEpEvgc+At91AVadExj7MDPVlt6YyxvfuQenGGFOblRogRKQ1zl/td+EEhxOAu4GVInJ0PBdX1f2qeqSq5pdIv0pVny6RtkVVh7qfN6hqV/fVSVUfqthj1Q6RdSL2k8xrobMZ6vmUDG+ejYkwxtQJZZUg/gxMVdUzVXW8qo5T1TNxJuh7uHqyV7f9vE5ECi+EBuKXENc1/tDGRBhj6gRR1dg7RL5R1Q4V3VeTsrKydOnSpTWdjVKtmXQ+Lfeu5czCJ0hPa2rTgBtjapyILCttpoqyShBltaRaK2sFzVyew1939+coyWOw5zPr8mqMqfXKGkndXESGxUgXoFmC8lNvTZq7li2BLmyQo/m1713eLOpT3OXVShHGmNqorACxGKd7aywfJyAv9dqWvAIUDy+EBvIn/wt0kQ2s1OOty6sxptYqNUCo6q+qMyP1XUZaCjl5BbwWOptbfNO5zvc2NwbGW5dXY0ytFc9IalMFIl1e95LKS6FzGer5lOO8O6zLqzGm1rIAUU2iu7w+HxyMIoxt9C7DumbUdNaMMSamUru51kW1vZtrtOxnryTth3kM9f6Lzfv9ZKSlWLdXY7FSH4EAACAASURBVEy1q2w31+gL9BKRy0TkisirarPY8Kxt92uayAGGFM5Bwbq9GmNqnXIDhIhMA57EmWDvLPd1ZlnnmPLd+5mPxaFOXOWbi58ggM30aoypVeJZUa430FFVw4nOTEOyJa+AZzwXMM37KMO9i5keOrs43RhjaoN4qphWAy0TnZGGJiMthYXhrqwJH8dvvW/iIVycbowxtUE8AaI58LWIvC0iMyKvRGesvnO6vfqYHBzOCZ6tDPF8Rorfa91ejTG1RjxVTDZzawJEeis9PsfLdwWv8XvfLD7y9OWWV1Ywae5a69FkjKlx5QaICq4eZypgRPdMRnTPZPH039J31b1kFX3OAnoU92iKHGOMMTWhrAWDPnDfd4nIzqjXLhHZWd6FRaS9iKyIeu0WkZtF5D4RyYlKH1rK+YNFZK2IrBeRCZV/xNrvznWnkK0tGeebCTjjUqxHkzGmppXVBnGO+94SSI96RbbLpKprVbWbqnYDegL7cda0BvhrZJ+qzi55roh4gcnAEKAjMEpEOsb5THXO5vwATwcvpIdnPX08q4vTrUeTMaYmlRogIt1aVTUU61XB+wwAvlPVTXEe3wtY7y49WgS8DAyv4D3rjIy0FF4Lnc2P2oJbfNOJlCKsR5MxpiZV11xMI4GXorbHichXIvKciLSIcXwmsDlqO9tNO4SIjBWRpSKyNDc3t+pyXI1uG9Qejz+FJ4MjOM3zLf08X1mPJmNMjUt4gBCRJGAY8Jqb9E/gBKAbsBV4PNZpMdJiThqlqlNUNUtVs9LTy635qpUiE/ktajKEbG3Jrb7X8HvhlldW0HfiApt+wxhTI+Kdi6mNiJzjfm4kIo0rcI8hwBequg1AVbe51VRh4Bmc6qSSsoFjorbbAFsqcM86Z0T3TBbeOYjszr+nm2cDWUWf2xxNxpgaFc9cTNcAbwJT3aTjgFkVuMcooqqXRKR11L5fAqtinPM5cJKItHNLICPdPNR7t6/vzKbwUdwa1RZhPZqMMTUhnhLEjTjzMe0GUNVvgaPiubiIpAIDgeiR14+KyEoR+Qqnp9Qt7rEZIjLbvUcQGAfMBb4GXlXV1TQAm/MD/D14EZ09Gxns+bw43Xo0GWOqWzwjqQ+oapGI0yzgdkGN1UZwCFXdDxxZIi3mUqaqugUYGrU9GzikC2x9l5GWwsy8M7kh/Ba3+V5hflEPgvisR5MxptrFU4JYLCK3A8luO8QrwP8Sm62G67ZB7Wnk9/NocCQneLZyuXchKX6P9WgyxlS7eALE7cAe4BvgJuA94K5EZqohi/Ro+rppXz4Nd+Bm3+s0kULr0WSMqXZlBgi3Ouk5Vf2nqv5SVUe4n21tiAQa0T2TxXcOYH+/e0mXfEaF3rIeTcaYaldmgHBHTLcWEX815cdEuXtpCrNDvbje9xYtyQesR5MxpvrEU8W0AfhIRO4UkRsjr0RnzDg9lyYFLyeJoDsFx8/pxhiTaPEEiFxgHpDKwZP2mQTLSEvhe23Nf0LnMdK7gFNkU3G6McYkmqjGnMGiTsrKytKlS5fWdDaqzMzlOdw5YyX+QD4LG93Kt3oMI4vuBoTMtBRbVMgYc9hEZJmqZsXaV+44CBGZR4x5kFT1/CrImylD5Mt/0ty1PL7nMh7yP8cQz2e8Ez7dFhUyxiRcPAPl7o76nAxcDBQmJjumpMiqc2c9HOLrgvnc5X+RBYXdKSSpuMHaAoQxJhHKbYNQ1U+jXh+o6o3EnmDPJFB2fhEPBH9FG9nB9d6fxylag7UxJlHimayvWdQrTUQGAK3LO89UrYy0FD4Jd+J/od783jeLY2VbcboxxiRCPL2YVuPMuLoaWI4zivq6RGbKHOq2Qe1J8Xt5MHAlRfh4wDcNUHLyCmyEtTEmIeJpgzheVQPRCSISz3mmCkU3WP9lzyX8yf8CQ0OfMjvc2xqsjTEJEU8J4tMYaZ9VdUZM+UZ0z2TxhHN5r8lwVoXbcq//BZqwH7AR1saYqldqgBCRo0SkK5AiIl1E5FT3dSbOoDlTQzbnF3FX4BqOIo8/+F4rTrcGa2NMVSqrqugC4Bqc5T6fikrfA9yTyEyZsmWkpfBl3om8EDqPMd53+V+oN8u0vTVYG2OqVLkjqUXkMlV9tcIXFmmPs3ZExPHAvUAmcCFQBHwHXK2qeTHO34gTjEJAsLSRftHq20jq0kRGWHsCe5nb6A4K1c/QoocpJMlGWBtjKqSskdRxTbUhIoOATjgD5QBQ1T9XIANeIAc4HWgPLFDVoIg84l7rjhjnbASyVHVHvPdpKAECnCAxae5ajt/9KS8kTeSp4DAeDY4EIMXv5eGLuliQMMaUq6wAEc84iKeAMcCtQApwJXBiBfMwAPhOVTep6rvumtMAS3CqsEwFRRqsNzQ7nVeC/Rnr/R+dZQNgDdbGmKoRTy+mM1X1CuAnVb0HpxRQ0S/1kcBLMdKvAd4p5RwF3hWRZSIytrQLi8hYEVkqIktzc3MrmK26b0teAQ8FR/MTzXjM/y8aUVScbowxhyOeAHEg8i4iR7vbbeO9gYgkAcOA10qk3wUEgRdLObWvqvYAhgC/F5F+sQ5S1SmqmqWqWenpDW8W8oy0FHbTmDsCY+ng2Vzcq0nBBtAZYw5LPAFitoikAY8BK4CNwPQyzzjYEOALVd0WSRCRMcAvgNFaSiOIqm5x37cDb2DzP8UUGWG9MNyNF4Lnca13Nmd4VgO2RKkx5vCUtya1B3hHVfNU9TWgHdBFVf9YgXuMIqp6SUQGA3cAw1R1fyn3bSwiTSOfgfNxpvswJYzonsnDF3UhMy2FPwevYKO24jH/0zS1AXTGmMNU3prUYeDvUdsFqroz3ouLSCowEJgRlfwk0BSYJyIrRORp99gMEZntHtMKWCQiX+KM2n5bVefEe9+GJtJgfYBkbg38jlbs4gH/v4v3W3uEMaYy4qlimiciwytzcVXdr6pHqmp+VNqJqnqMqnZzXze46VtUdaj7eYOqdnVfnVT1ocrcv6HJSEthhZ7IE8GL+KV3MRd7PgSsPcIYUznxBIhxwBsiUiAiO0Vkl4jEXYow1SfSHvFkaASfhDryoP/fnCBOULD2CGNMRcUTIFoCfqAJkO5uN7zuQnVApD2idVpjbgr8ngKSeNL/RHHXV2uPMMZURDwryoWAS4E73M+tgW6JzpipnEh7RC4tuDXwO07xbOZe3wvF+609whgTr3hGUj8JnAP8yk3aDzydyEyZw5eRlsIH4a48HbyQ0b73rD3CGFNh8VQx9VHV63EHzLm9mJISmitz2CLtEZOCl/FxqCMP+Z+lk3wPWHuEMSY+8QSIgDseQgFE5EggnNBcmcMWaY84Oq0J4wI38hPN+FfSX0ljD2DtEcaY8sUTICYDrwPpInI/sAh4JKG5MlUi0h6xi2b8tuhm0snjCf+TeAkB1h5hjClbPI3U/wfcjTPVxk7gUlV9OdEZM1UnIy2Fr/QE7g5eQz/vSu72/Qew9ghjTNniKUEAeIEAziI/8Z5jaolIe8Rrof5MDQ7hat9cRnvnA9YeYYwpXTy9mO7CmUspA2ea7/+KyJ2JzpipOgfP1zSaBaFu3O+bRh+PM72VtUcYY2KJZ8nRr4GekYn13PmVlqnqKdWQvwppSCvKVVa7CW/TmP1MT7qf1vITlxTdxzp1lvew5UqNaXgOa0U5YBPgi9r2ARuqImOm+mWkpbCXVH5T9P8oJIlpSY/QCmfmFKtuMsZEiydA7AdWi8hUEXkGWAnkichfROQvic2eqWqR9ogc0rmq6HaasZ9pSY/Y9ODGmEP4yj+Et91XxJIE5cVUg0j10aS5a1mT15YbAjczzf8o//L/hasDt1NIknV/NcYAcbRB1CXWBlExfScuICevgOGeRfw96SnmhXry28BNBPFZe4QxDcRhtUGIyGAR+VxEttt03/VLpLppVvhM7g5czUDvMh7zP42HsLVHGGPiaoN4ErgeyKQC032LSHt3xbjIa7eI3CwiR4jIPBFZ5763KOX8wSKyVkTWi8iEijyUiU9099f/hAYyMTCSEd6P+f98zwFq7RHGNHDxdHNdCJzrLj9auZuIeIEc4HTg98BOVZ3ofvG3UNU7Yhz/Lc5ypdnA58AoVV1T1n2siqny2k14GwX+n+8Vxvlm8ULwPO4NXoXiQXB6P1mVkzH1T1lVTPE0Ut8OvOUGisJIoqo+UYE8DAC+U9VN7vKl/d3054GFwB0lju8FrFfVDe4DvAwMB8oMEKbyMtJSyMkr4LHgZfgIc4PvLQTlnuDVKJ7iKifAgoQxDUQ8VUz3AyEgDadqKfKqiJE4o7EBWqnqVgD3/agYx2cCm6O2s920Q4jIWBFZKiJLc3NzK5gtExFpjwBhYnAkTwWHcaXvPR7yPYe4k/dalZMxDUs8JYijVLVnZW8gIknAMKAi03NIjLSYdWGqOgWYAk4VU4UzaICDu79uySvg0eDlhBHG+WbRSIq4IzCWID5y8groO3GBVTcZ0wDEEyDeE5FzVXVBJe8xBPhCVbe529tEpLWqbhWR1sD2GOdkA8dEbbcBtlTy/iZOI7pnFn/p9524gMfyLqNAG3Gb/1WaUcC4wHgKSbLqJmMaiHiqmK4D5ovI3kp2cx3Fz9VLAG8CY9zPY4BZMc75HDhJRNq5JZCR7nmmmjhVTj4mh0ZwT+AqBnqX8W//ozSJGnF98ysrbLpwY+qxeAJES8APNKcC3VyheGK/gcCMqOSJwEARWefum+gemyEiswFUNQiMA+YCXwOvqurqeO5pqkZ0F9gXQudzU9Hv6OX5hleTHiyeuwls/iZj6rO4RlKLyEjgeFX9s4i0wWloXpbw3FWQdXNNjMiI636eL3nK/3d2k8pVRXfwrf5cC5iZlsLiCefWYC6NMZVxuCOpnwTOAX7lJu0Hnq667JnaLtLD6cNwVy4ruhcvYaYn3ceZnpXFx0Qar60kYUz9EU8VUx9VvR44AKCqO4GkhObK1CrR1U1rtC2/LHyAHG3JNP8jXOWdQ6SDmVU3GVO/xBMgAiLiwf0WEJEjgUqPqjZ104jumSyecC5/u7wbu/ytuKToPhaEu3Of//942DcVP0HAGq+NqU9KDRAiEukCOxl4HUgXkfuBRcAj1ZA3UwtFShNpaUdwfeAW/hEcwSjf+7yS9ABH81PxcVaaMKbuK7WRWkS+UNUe7udOwHk4A9jmq+qq6sti/KyRunpFGq+HepbwqH8KhfgZHxjPx+HOxcdY47UxtVtlG6mLRzOr6mpV/buq/q22BgdT/SKN17PDvRle9CA7tRkv+B9mvHcGHrcW0hqvjam7yhpJnS4it5a2U1VtudEGLnp6ju/yMhle9CAP+Z/lD/7p9PGs4ebA79jGETby2pg6qqwShBdoAjQt5WXMQY3X6m/MLYHf8YeiGzjV8x3vNJrAQI9T5WeN18bUPWWVILaq6gPVlhNTp0WXJl7P68fyohN5wv8kzyT9hdeC/Xgg+Gv2kGqlCWPqkLIaqZeravdqzs9hsUbq2iHSeO0nyI2+GfzOO4utHMntgbGHNGDbrLDG1KzKNlIPSFB+TD0XabwO4OPx4GVcXHQ/hernv0l/5hHfFJqxF7CusMbUdnHNxVRXWAmi9pi5PIdJc9eSk1cAQCOKuMk3g7He/7GLptwXGMPb4dOJdJaz0oQxNaOsEoQFCJNQM5fncOeMlRQEQgB0ko1M9E+hi2cjH4a68KfgVXyvrQEnVCgWLIypToc1WZ8xhyN6HieA1dqWEUUPcl/g13TzrGdO0h38P98rpHKgeMlAq3oypnawEoSpNiVLE+nkcaf/v1zkXcQ2TePRwEhmhM9Eo/5usdKEMYlVY1VMIpIGTAU649QeXAPcDLR3D0kD8lS1W4xzNwJ7gBAQLO0BolmAqP1Ktk0A9JBvucf/H7p71rMy3JZHgqNYFO5SvD/F7+Xhi7pYkDAmAWoyQDwPfKSqU92lQ1NVNS9q/+NAfqzxFm6AyFLVHfHezwJE3VGyNCGEudDzCbf7X6GN7GBRqBOPBkfylZ5QfI6VJoypejUSIESkGfAlzkp0h9xERAT4AThXVdfF2L8RCxD1WnRpItJAnUSA0d75jPPN5EjZw7xQD/4WvJjV2g6whmxjqlpNBYhuwBRgDdAVWAbcpKr73P39gL+UmjGR74FdON8H/1LVKeXd0wJE3VWy6qkJ+7nKO5frfG/TXPYzL9STJ4PD+VJPLD7HgoUxh6+mAkQWsAToq6qfisjfgd2qeo+7/5/AelV9vJTzM1R1i4gcBcwDxqvqhzGOGwuMBTj22GN7btq0KSHPY6pHyaqnpuznau8crvG9Q5rsY3GoE0+FhrE43JmoCYetncKYSqqpAHE0sERV27rbZwETVPUCdzGiHKCnqmbHca37gL2q+lhZx1kJon6I1ZDdmAKu8L7Htb7ZtJI81oSPY2pwCG+F+xBwpxTzihBWJcNKFMbErSYbqT8CrlXVte6XfGNVvU1EBgN3qurZpZzXGPCo6h738zzgAVWdU9b9LEDULyVLE+CMyB7hXcxvvLM52ZPDdk3jv6Fz+W9wANtpUXycVT8ZE5+aDBDdcLq5JgEbgKtVdZeITMMpXTwddWwGMFVVh4rI8cAb7i4f8F9Vfai8+1mAqH9iNWQ7lLM8K7nG+w5ne74ihIe54Sz+GxrAJ+GOB42lsGBhTOlsqg1TL5QWLI6VbVzpnc9l3oWkyT42hY/ilVB/ZoTO4keOPOgaFiyMOZgFCFPvRILFlrwCPCKEVGlEEYM9nzHSu5AzvGsIq7Ao3JkZobN4N5zFfpIPuoYFC2MsQJh6LlZbxXHyIxd5F3GR5yOO8eSyXxsxP9yDN0N9+DB8KkX4D7qGBQvTUFmAMPVeadVPQpjTZC3DvB8z1PspR8hedmsKC8LdeSd0Oh+ET+UAjQ66lgUL05BYgDANSmnBwkeQvp7VDPF8yiDvUlrIXgo0iUXhLswL92BBqAc7aH7QtSxYmPrOAoRpsMoKFqd7vmagZxnneb+gjTgzunwVbsf74W58EOrKl3oCIbzF14qcn5biRwTy9gdszIWp8yxAGEPZXWY7yib6e1ZwjncFPWQdXlF2ayofhzuxKNyZT8Id+U4ziB69HWGlDFOXWYAwpoTSgwU0Zy99Pas4y7OSs7wri0sX2zWNJeFT+Czcgc/CHVinmQeNtwArZZi6xwKEMWUoK1iAcpxs4wzPGs7wrOF0z9ccLbsAyNPGfBE+iWXhk/lCT2JluB17SY15D79HaJLss4Bhah0LEMbEqexgAaAcI9s53fMNWbKWHp51nOxxlkYNq7BeM/gyfAIrtR2rwu1Yo8cd0ksKrKRhag8LEMZUQvnBwtGMvXTzfEdX+c5593xHS9kNQEiF7zSDNXocX4ePY60ewzfhY/iRIyirPSM6cDS3IGISyAKEMYcpeuR25At71/5AqaWM1uyks+d7Onu+p6P8wCmeTcVtGeBUT32rbVgfznTeNZPvwhlsLSVwRLPSh6lKFiCMSZB4SxngNH53kM209/xAB9nMiZ4cTpZs0mRf8TH7tBEb9Wi+16P5XluzSVuxKdyKTdqK7aRRVvCw0oepDAsQxlSDWKWMyJf0vqIggVCs/2tKOvmc4NnC8bKVE2QL7WQrbeVHjpXt+CRcfOQB9ZOt6WzWdLI1nRxtSbams0WPJEdbkksa4RK9qkqKJ4ic0yGd97/JZUtegQWVBsAChDE1rCIljQgfQTJlB8fJNo6V7Rwr2zlGtnOM5NJGcg8qeQAE1Ms2WrBNW/CjtuBHPZJtmsY2bcF2WrBd08jV5uymMeVVY0UrL6iU9tmCS91gAcKYWqRi7Rmla0wBmbKDDNlBhuwkQ3bQWnbSip20lp0cLTtpLIWHnFeoPnbQnB3qvHZqU36iOT9pU3ZqM3bSlF3alF00YZc2YQ+ph4z3iFdlg4sFnepjAcKYOqCsKqrDCSKtZBetZBfp5JMueaRLPi0ln5Y470fIbo5kD40kEPMaIRXyaUyeNiGfJuRrY/JpTL42Zjep7NZU9pDKnuL3FHaTyj5NYS8p7CO53KqviogeU3I4wSYRn+tiAKvJFeXScFaU64zz7/oaYBBwHZDrHvZHVZ0d49zBwN8BL85KcxPLu58FCFPfVVXp41BKEwpoIXs4gj0cIXtowR5ayF7SZC/N2Uea7CWNvTSTfTRnH81lH83Yf1A7SWn2aSP2kcJeTWY/yewjmX2azH4asd99L6AR+7UR+2nEARpRoEkU0IgCkjhAIw5oEgdIcrY1iUL8HMB5r2wJJxGip16Jbs+prcGpJgPE88BHqjpVRJKAVOBmYK+qPlbGeV7gW2AgkA18DoxS1TVl3c8ChGmoElH6iI+SQiHN2E9T2e++F9CEAppIAU3ZTxMpoDEHaEwBTeQAqRygMYU0lgJSKSRVDpBKISkU0kiClcpFofooxE8hforwU6jOexE+512d9wA+Ct20gPoI4CWAjyJ8BPAdlBaMeg/iJaBegsXpXkKRd3Xeg25aEA8hPD9/VmdfGA9BPISL9zvHhPAQRqhIu1BpUvxeHr6oS4WCRFkBwnfYOSr9ps2AfsBVAKpaBBSJxPVD6AWsV9UN7rVeBoYDZQYIYxqqEd0zy/1SKC+IlOzFFF+DulBAMgUks02PcJIOIwp5CZFKIckUkiJF7uciUsR5L35JEY0IFH9OIkgjAjTC/SwBGhEgKfKSIKkcoBFBkgjgI0SSx/nsJ4SfIH6CcZWGEiWk4gaLnwNIuPhdCEdtqwphnOMVYQfNubzoXgoCISbNXVtlVVwJCxDA8TjVSP8Wka7AMuAmd984Efk1sBT4g6ruKnFuJrA5ajsbOD3WTURkLDAW4Nhjj6263BtTz8QTREqKJ6hUZYklhNdpxyD14JOrqalUCOMnhK84aITwEcQvQXyEi9O9hJ2AQgivRH0mXLzP4372yc/pkWMi+zyE8RHGKweneaM+e9CDPnsI4xE9aHuPphQ/w5a8gir7eSQyQPiAHsB4Vf1URP4OTACeBB7E+ZU/CDyO0zYRLVYxI+Y/EVWdAkwBp4qparJujIHKBZWIygaX8j6XPqbk8CkeivAcsiRtmQGqln3rZKSllH9QnBIZILKBbFX91N2eDkxQ1W2RA0TkGeB/pZx7TNR2G2BLojJqjKl6hxNcypKowHO4nxPXzhO/FL+X2wa1r7LrJSxAqOqPIrJZRNqr6lpgALBGRFqr6lb3sF8Cq2Kc/jlwkoi0A3KAkcAVicqrMabuSFTgqQrRwavkqPTa2oupLIksQQCMB150ezBtAK4GnhCRbjiBdiNwPYCIZOB0Zx2qqkERGQfMxenm+pyqrk5wXo0x5rDU5uBVGTZQzhhjGrCyurnWntElxhhjahULEMYYY2KyAGGMMSYmCxDGGGNiqleN1CKSC2yq5OktgR3lHlW/2DPXfw3tecGeuaKOU9X0WDvqVYA4HCKytLSW/PrKnrn+a2jPC/bMVcmqmIwxxsRkAcIYY0xMFiB+NqWmM1AD7Jnrv4b2vGDPXGWsDcIYY0xMVoIwxhgTkwUIY4wxMTX4ACEig0VkrYisF5EJNZ2fRBORY0TkfRH5WkRWi8hN5Z9VP4iIV0SWi0isNUjqHRFJE5HpIvKN+/s+o6bzlGgicov773qViLwkIsk1naeqJiLPich2EVkVlXaEiMwTkXXue4uquFeDDhAi4gUmA0OAjsAoEelYs7lKuCDOMq+nAL2B3zeAZ464Cfi6pjNRjf4OzFHVDkBX6vmzi0gmcCOQpaqdcZYKGFmzuUqIacDgEmkTgPdU9STgPXf7sDXoAAH0Atar6gZVLQJeBobXcJ4SSlW3quoX7uc9OF8a9WcC+1KISBvgAmBqTeelOohIM6Af8CyAqhapal7N5qpa+IAUEfEBqdTDlShV9UNgZ4nk4cDz7ufngRFVca+GHiAygc1R29k0gC/LCBFpC3QHPi37yHrhb8DtQLimM1JNjgdygX+71WpTRaRxTWcqkVQ1B3gM+AHYCuSr6rs1m6tq0yqyUqf7flRVXLShBwiJkdYg+v2KSBPgdeBmVd1d0/lJJBH5BbBdVZfVdF6qkQ/oAfxTVbsD+6iiaofayq13Hw60AzKAxiJyZc3mqm5r6AEiGzgmarsN9bBIWpKI+HGCw4uqOqOm81MN+gLDRGQjTjXiuSLyn5rNUsJlA9mqGikdTscJGPXZecD3qpqrqgFgBtCnhvNUXbaJSGsA9317VVy0oQeIz4GTRKSdu272SODNGs5TQomI4NRLf62qf6np/FQHVb1TVduoaluc3/ECVa3Xf1mq6o/AZhFp7yYNANbUYJaqww9AbxFJdf+dD6CeN8xHeRMY434eA8yqiov6quIidZWqBkVkHDAXp8fDc6q6uoazlWh9gV8BK0VkhZv2R1WdXYN5MokxHnjR/eNnA3B1DecnoVT1UxGZDnyB01tv+f/f3v27NhWFYRz/PoOCOAji5KD/gIJQEASFDm6OKg4uopOLW1cdVdykmyD+ooKbg6ObOHRTN5d2EHeh1UHxdbg3NsipGrmxof1+INyT5N7kZEge7g3nfdmGZTeSPAXmgQNJPgA3gFvAsyRX6ILy/CDvZakNSVLLTr/EJEnahAEhSWoyICRJTQaEJKnJgJAkNRkQ0gSSrP1y/1KSxa2ajzRNBoQ0A/rKwtJMMSCkgSQ5nORlkrf99lD/+IMk58b2W+u3831vjiW6hYt7k7xI8qbvZ3Bhiz6KBOzwldTSP9gztgIdYD8b5VkWgUdV9TDJZeAufy67fBw4UlUrSc4CH6vqDECSfQPPXZqIZxDSZL5U1bHRDbg+9twJYKkfPwZO/sXrLVfVSj9+B5xOcjvJqar6NNy0pckZENL0jOrYfKP/rvVF5HaP7bP+c+eq98AcXVDcTDIePtJ/Z0BIw3nNRovLi8CrfrxK98MPXb+CXa2DkxwEPlfVE7rGN9u9PLdmnP9BSMO5BtxPskDXzW1UPfUe8DzJ1t9gGQAAAEFJREFUMl2/4PVNjj8K3EnyHfgKXJ3yfKXfspqrJKnJS0ySpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnpB2mRU47mpIK1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(t,Temp_numerical,'o-',label='100'+' Euler steps')\n",
    "plt.plot(t,Temp_analytical,label='analytical')\n",
    "plt.title('Body Temperature vs Time')\n",
    "plt.xlabel('Hours')\n",
    "plt.ylabel('Temperature in Degrees F')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [],
   "source": [
    "#b\n",
    "import numpy as np\n",
    "\n",
    "T_0=85\n",
    "T_ambient = 65\n",
    "K = 0.61111111111112\n",
    "t = np.linspace(0,1000,1000)\n",
    "Temp_analytical = T_ambient + (T_0-T_ambient)*np.exp(-K*t)\n",
    "delta_t = 1\n",
    "Temp_numerical=np.zeros(len(t));\n",
    "for i in range(1,len(t)):\n",
    "    Temp_numerical[0] = T_0\n",
    "    Temp_numerical[i]=T_ambient+(Temp_numerical[i-1]-T_ambient)*np.exp(-K*delta_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f97583a3c88>"
      ]
     },
     "execution_count": 161,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZyVZf3/8debAQTcQBmNJYPSL4aaoCMhlGFoKJZQ5paWS2p+v6mZ34eGaWX1Ky21xbQFcauvuRGCJYqaoVluIKi4IIqoDISjCCiyDMPn98d9n+Ewc2bmPuOcGZh5Px+P45z7urfrOoPnM9dyX5ciAjMzs7o6tXUGzMxsy+QAYWZmBTlAmJlZQQ4QZmZWkAOEmZkV5ABhZmYFOUBYq5I0QFJI6tzWebHSkfR1SX9t63zYB+MAYU2StEjSGknvSXpH0t2SPtzKeXguvf97kmokrc3b/m5r5qWlSXpM0oltnY9iSHqlkd/HeRFxXUR8oa3zaR+MA4Rl9YWI2A7oAywDftOaN4+IvSJiuzQP/wTOym1HxE9bMy/FaI2aUlvUxiLiY3m/jyeB0/J+H79o7fxYaThAWFEiYi0wGRicS5O0o6Q/SqqS9JqkiyV1SveVSbpC0luSFgJH5J13tKTZ+deX9L+SpjYnb5K+IWm+pOVpLadfmt4tbdY6M/3Ld1Wax0GSnpC0UtLNuS9aSYdJelnSD9NrLZR0dN59ukv6laQ3JP1H0m8kbVPn3O9JWgb8TlK5pHvSz2e5pGmS+qTHXwkcAExK//q+UtKekjbUKVttLSMtx4OSrpH0DjChsfIX+JxmSjqtTtp8SWPT39fVaV5XSnpa0qBm/C7OlPRAcz7/9JwvSnpG0gpJ/5Q0uOG7Wak4QFhRJPUAjgUey0v+DbAj8FHgM8DXgFPSfacDnweGAhXAl/POuwsYKOnjeWknAn9qRr6OA84FvgDsCswB/q/OYaOBfdM8/iDN99HAQGAYcFTesQOArsCHgDOAmyQNTPf9EugP7AMMAv6L9Es679wuwIeBc0j+P/s9sFt6r9w1iIj/ZfO/wP83Y5EPAuYCvYErM5Y/58/A8bkNSfsDOwH3kfyu9gc+BvQCvgK8kzFPTcn0+UsaDvyW5N/QziT/Hqa2RU2pw4sIv/xq9AUsAt4DVgAbgCXAPum+MmAdMDjv+G8AM9P3DwJn5u37HBBA53T7d8BP0vd7kXwZbdNEfmaSfKHmp/0DOCFvuwtQTfJl2S295/55+58DvpW3fQ1wWfr+MGAt0C1v/13A+UBnYD3QL2/fwcALeeeuBro0kv/hwNK87ceAE/O29wQ21Dmn9hjgTOClrOUvcP+dgDVAn3T7SuC36fux6WczDFDGfx+b5T8vjw+k74v9/G8ALqpzvdeAT7b1/wsd7eUahGU1PiJ6AtsAZwEPSfoQyV+wXUn+B855Dcg1b/QF3qizL99NwFckCfgqcHtErGtG/j4C/D5tklgBVJEEs/55xyzLe7+mwPZ2edtVkTSn5ee7b/rqAjyXd6+pwC55x/4nIqpzG5K2l3S9pNclrSL5S713M8qY740621nKD0BELAfuB45JmwKPBW5Od98DXAf8AVgm6beStqt7jWbK+vl/BPhurixpecrZ9G/KWokDhBUlImoiYgpQA3wKeIvkL9WP5B22G1CZvl9K0tSSvy//eo+R/EX+aZLmjKKbl1JvACdHRM+8V/eImN3kmYX1ltQtb3s3kprTUpIv3o/l3WfHiNg579i6UyRPIPmiPiAidiCpRamR41cDZbl+jdSH6hxT95xiy38LSTPTZ9Ly/BsgEr+IiKHAJ0iahL7VwDVK5Q3g+3XK0iP9d2etyAHCiqLEOJL26Rcioga4HfhJ+pfyR4Dz2NT+fTtwjqT+knqxeVt9zh+Bq0maVR5pZtZ+D1yc61CV1EvSUU2c05guwPckdZX0WeBQ4C9pzeB64NeSeqefx4clHdrItbYH3gdWSOoNXFxn/zKS/pucJSQ1gBPSTuP/oem/nost/zSSJr2LgFsiknYcScMlVaTt/atJgndNE/duaROBs9N8SNJ2ko5M+7+sFTlAWFZ/lfQesAr4CXBSRDyX7jub5MtkIfAISSfo9em+a4EZwNPAU0ChvwL/BOxN82sPRMQtJEFmStqMM5fkS725FpH8Zf0fkrKcEhEL033nknyJzwJWAvcCuzdyrStImpTeJvl8ptfZ/0vga0qeMfl5GnRPI+nIfYukBtZoTajY8kfE+yT9KqNJfl85PYEbSfqbFpI0rV3V2L1bWkT8i6Rz/w9pPl4iqV168ZpWpvQPB7M2I6k78CawX0Qs2ALycxhwdUQ09qVv1u65BmFbgv8GntwSgoOZbeJxxdamJC0i6bAd38ZZMbM63MRkZmYFuYnJzMwKaldNTL17944BAwa0dTbMzLYas2fPfisiygvta1cBYsCAAcyaNauts2FmttWQVHd2g1puYjIzs4IcIMzMrCAHCDMzK6hd9UGYWelUV1ezePFi1q5d2/TBtsXp1q0b/fv3p0uXLpnPcYAws0wWL17M9ttvz4ABA0hmZ7etRUTw9ttvs3jxYgYOHNj0CamSNjFJ+raSxebnSbolXXrwEkmVkuamr7ENnHtYugziy5IKzQDaIqbOqWTkZQ8ycMLdjLzsQabOqWz6JLMOaO3atey8884ODlshSey8885F1/5KVoNI18M9h2SlsTWSbgeOS3f/MiKuaOTcMpIVpg4FFgNPSrorIp5vyTxOnVPJhVOeZU11Mptx5Yo1XDjlWQDGD/XaJGZ1OThsvZrzuyt1J3VnoHs6t3wPkimSsxgGvBwRCyNiPXArMK6lM3f5jPmsqa7h7LIpHNTpaQDWVNdw+Yz5LX0rM7OtTskCRERUksyD/zrJKlwrI+K+dPdZkp5Jl2HsVeD0fmy+pOJiGlgwRdIZkmZJmlVVVVVUHpesWAPA/3S+i5Gd5tVLN7Mtx6mnnsouu+zC3nvvvVn68uXLOfTQQ9ljjz049NBDeeedd2r3XXrppey+++4MGjSIGTNmFLzuqFGjGDRoEEOGDGHIkCF8+ctfbjQfixYtqpeHlnDjjTeyZEnWv6FbR8kCRPrFPw4YSLKO77aSTiRZpP5jwBCSwHFlodMLpBWcVTAiJkZERURUlJcXfFq8QX17dk8vrM1umEs3s+Zr6f69k08+mXvvvbde+mWXXcbo0aNZsGABo0eP5rLLLgPg+eef59Zbb+W5557j3nvv5X/+53+oqSm8ON7NN9/M3LlzmTt3LpMnT/5A+ayroXvW1aECBHAI8GpEVKXLNE4BRkTEsnRd440kq40NK3DuYjZfx7g/2ZunMjt/zCC6dykjAKXxp3uXMs4fM6ilb2XWoeT69ypXrCHY1L/3QYLEQQcdxE477VQvfdq0aZx00kkAnHTSSUydOrU2/bjjjmObbbZh4MCB7L777jzxxBOZ73fyySdvFiy22267esfU1NRw/vnnc8ABB/CJT3yCP/zhDwDMnDmTgw8+mK985Svss88+9c45+eST2Xvvvdlnn3345S9/yeTJk5k1axYnnHACQ4YMYc2aNcyePZvPfOYz7L///owZM4alS5cCSY3n3HPPZcSIEey99961ZXrooYdqa0FDhw7l3XffzVzWhpRymOvrwPB0Hdk1JEsbzpLUJyKWpsd8EZhX4NwngT0kDQQqSTq3v9LSGcx1RMfUpP7Qr2d3zh8zyB3UZk344V+f4/klqxrcP+f1Fayv2bhZ2prqGi6Y/Ay3PPF6wXMG992BH3xhr6LzsmzZMvr06QNAnz59ePPNNwGorKxk+PDhtcf179+fysrCAeqEE06ge/ek5eDQQw/l8ssvz3Tv6667jh133JEnn3ySdevWMXLkSD73uc8B8MQTTzBv3rx6w0rnzp1LZWUl8+YlX30rVqygZ8+eXH311VxxxRVUVFRQXV3N2WefzbRp0ygvL+e2227joosu4vrrk5V8V69ezb///W8efvhhTj31VObNm8cVV1zBNddcw8iRI3nvvffo1q1b1o+wQSULEBHxuKTJJOsQbwDmkCxGPknSEJImo0XANwAk9QUmRcTYiNgg6SyStYzLgOvz1j9uUeOH9mP1VNin3w78678/W4pbmHU4dYNDU+mlUGitm4ZG8tx8881UVFQUfY/77ruPZ555pramsXLlShYsWEDXrl0ZNmxYwWcOPvrRj7Jw4ULOPvtsjjjiiNqAkm/+/PnMmzePQw9NlhWvqampDYIAxx9/PJDUqlatWsWKFSsYOXIk5513HieccAJf+tKX6N+/f9HlqaukD8pFxA9IFl7P99UGjl0CjM3bnk79xd1LIhB44SSzzJr6S3/kZQ9SWWCwR7+e3bntGwe2aF523XVXli5dSp8+fVi6dCm77LILkNQY3nhj01iXxYsX07dv38zX7dy5Mxs3JgEtIli/fn29YyKC3/zmN4wZM2az9JkzZ7LtttsWvG6vXr14+umnmTFjBtdccw233357bc0g/7p77bUXjz76aMFr1A10kpgwYQJHHHEE06dPZ/jw4TzwwAPsueeemctbiOdiItf77QBh1lJy/Xv5StW/d+SRR3LTTTcBcNNNNzFu3Lja9FtvvZV169bx6quvsmDBAoYNK9TlWdiAAQOYPXs2kPRnVFdX1ztmzJgx/O53v6vd99JLL7F69epGr/vWW2+xceNGjjrqKH784x/z1FNPAbD99tvX9hsMGjSIqqqq2gBRXV3Nc89takS57bbbAHjkkUfYcccd2XHHHXnllVfYZ599+M53vkNFRQUvvvhi5rI2xFNtACHXIMxaUq4f7/IZ81myYg19W6B/7/jjj2fmzJm89dZb9O/fnx/+8Id8/etfZ8KECRxzzDFcd9117Lbbbtxxxx0A7LXXXhxzzDEMHjyYzp07c80111BWVlbw2vl9EL179+aBBx7g9NNPZ9y4cQwbNozRo0cXrBGcdtppLFq0iP3224+IoLy8vLaTvCGVlZWccsoptbWTSy+9FEg6xc8880y6d+/Oo48+yuTJkznnnHNYuXIlGzZs4Nxzz2WvvZKaW69evRgxYgSrVq2qrX386le/4h//+AdlZWUMHjyYww8/vBmf8uba1ZrUFRUV0ZwFg1Zd0pcXyg/nk9+8rgS5MmsfXnjhBT7+8Y+3dTY6vFGjRtV2Zher0O9Q0uyIKHgxNzGRNi61o0BpZtYS3MRE2kltZrYVmDlzZqvdyzWIWq5BmJnlc4DANQgzs0IcIAAQROs9wGNmtjVwgMCNS2ZmhThAkGticpgw6whuvPFGzjrrrCaPyZ9Z9bTTTuP554tfr2zmzJl8/vOfL/q8LYUDBJ5qw8w2VzdATJo0icGDB7dhjtqGA0RKrkGYbfHGjx/P/vvvz1577cXEiROBZBruiy66iH333Zfhw4ezbNkyAP7617/yyU9+kqFDh3LIIYfUpue8++67DBw4sHaajFWrVjFgwADuuOOOelNvjxo1itxDuPfeey/77bcf++67L6NHjwaSmVtHjBjB0KFDGTFiBPPnt49VKf0cBB7FZFa0eybAf55t2Wt+aB84/LJGD7n++uvZaaedWLNmDQcccABHHXUUq1evZvjw4fzkJz/hggsu4Nprr+Xiiy/mU5/6FI899hiSmDRpEj//+c+58spN65Ntv/32jBo1irvvvpvx48dz6623ctRRR3H00UdzzTXXFHxauaqqitNPP52HH36YgQMHsnz5cgD23HNPHn74YTp37swDDzzAd7/7Xf7yl7+07OfTBhwgctzEZLbFu+qqq7jzzjsBeOONN2qn1s618++///7cf//9QDJ767HHHsvSpUtZv359wam3TzvtNH7+858zfvx4brjhBq699tpG7//YY49x0EEH1V4rt4DRypUrOemkk1iwYAGSCk7stzVygMCd1GZFa+Iv/VKYOXMmDzzwAI8++ig9evRg1KhRrF27li5dutROf11WVsaGDRsAOPvssznvvPM48sgjmTlzJpdcckm9a44cOZJFixbx0EMPUVNT0+Ra0xFRcE2J733vexx88MHceeedLFq0iFGjRn3g8m4J3AeR4xqE2RZt5cqV9OrVix49evDiiy/y2GOPNXl8v37J7LG56cAL+drXvsbxxx/PKaecUpuWP/V2vgMPPJCHHnqIV199FaC2iSn/XjfeeGNR5dqSlTRASPq2pOckzZN0i6Ruki6X9KKkZyTdKalnA+cukvSspLmSip+itQjugzDb8h122GFs2LCBT3ziE3zve9/bbDnRQi655BKOPvpoPv3pT9O7d+8GjzvhhBN45513aldpg01Tb+c6qXPKy8uZOHEiX/rSl9h333059thjAbjgggu48MILGTlyJDU1NR+wpFuOkk33Lakf8AgwOCLWSLqdZIW4JcCD6bKiPwOIiO8UOH8RUBERb2W9Z3On+/7PJR/jjZ4HcMC5txZ9rllH0V6n+548eTLTpk3jT3/6U1tnpeSKne671H0QnYHukqqBHsCSiLgvb/9jwJdLnIcM/ByEWUd09tlnc8899zB9equsbrzVKVmAiIhKSVcArwNrgPvqBAeAU4HbGroEcJ+kAP4QERMLHSTpDOAMgN122615eW3WWWa2tfvNb37T1lnYopWsD0JSL2AcMBDoC2wr6cS8/RcBG4CbG7jEyIjYDzgc+KakgwodFBETI6IiIirKy8s/QI4dJsya0p5WoOxomvO7K2Un9SHAqxFRFRHVwBRgBICkk4DPAydEA7mOiCXpzzeBO4Hsq40XKeRhrmZN6datG2+//baDxFYoInj77bfp1q1bUeeVsg/idWC4pB4kTUyjgVmSDgO+A3wmIt4vdKKkbYFOEfFu+v5zwI9Kl1Uh/6M3a1T//v1ZvHgxVVVVbZ0Va4Zu3brRv3//os4pZR/E45ImA0+RNCXNASYCzwHbAPenD5w8FhFnSuoLTIqIscCuwJ3p/s7AnyPi3pLlNe+/ZlZYly5dCj6NbO1XSUcxRcQPgB/USd69gWOXAGPT9wuBfUuZt83u7ecgzMzq8ZPUgIe5mpnV5wBBUoPwdN9mZptzgDAzs4IcIGq5BmFmlq/BACHpv/Pe79k62WkbIfdBmJnV1VgN4vS8938udUbakvsgzMzqy9rE1M7HgfpJajOzuhp7DqKnpC+QBJEdJB2ZvzMi7ippzlpR1P7HzMxyGgsQ/wKOSd//Gzg6b18A7SZAuAZhZlZfgwEiIr7amhlpS36S2sysPg9zTbmT2sxscw4QgKfaMDOrzwECz+ZqZlZIkwFC0pckbZ++nyDpdklDSp+11pMsGGRmZvmy1CAuSRfuGQF8gWQN6d+XNlutzQ/KmZnVlSVA1KQ/Pw/8NiL+QrLgT7sR7oMwM6snS4BYKuka4FhguqSuGc9D0rclPSdpnqRbJHWTtJOk+yUtSH/2auDcwyTNl/SypAnZi9Q8rkGYmW0uyxf9McBDwBER8Q7QG2jyC1tSP+AcoCIi9gbKgOPSc/8eEXsAfy90LUllwDXA4cBg4HhJgzOVqBn8HISZWX1NBoiIeA94HRiWJq0jWVc6i85Ad0mdgR7AEmAccFO6/yZgfIHzhgEvR8TCiFgP3JqeV0KuQZiZ5csyiuliknWlL06TupFhdteIqASuIAkuS4GVEXEfsGtELE2PWQrsUuD0fsAbeduL07RC+TtD0ixJs6qqqprKVmGe7tvMrJ4sTUxfBsYCq6H2i3+Hpk5K+xbGAQOBvsC2kk7MmK9CbT4Fv8EjYmJEVERERXl5ecbL172w3MhkZlZHlgCxLiKC9AtaUo+M1z4EeDUiqiKiGpgCjACWSeqTXqsP8GaBcxcDH87b7k/SPFVCrkGYmeXLEiCmpKOYdpR0CnAfcH2G814HhkvqIUnAaOAFkllgT0qPOQmYVuDcJ4E9JA1MR00dRwlnjw3P5mpmVk9j030DEBE/k3Q4sB7YF/hJRNyT4bzHJU0GngI2AHOAicB2wO2Svk4SRI4GkNQXmBQRYyNig6SzgBkko5+uj4isHePNIOQ+CDOzzTQZIFLPAmsj4h/pswzbRsTqpk6KiB+QdHDnW0dSm6h77BKSvo7c9nRgesb8fSAe5mpmVl+WUUynkjTvTEqTdqNws9BWzjUIM7N8WfogzgGGA6sAIuIlCg9N3XrJczGZmdWVJUCsTR9WA2qfcm5XbTIODWZm9WUJEP+SdAHQTdLBJLO5/q202Wpt7qQ2M6srS4C4AHgXeBH4Fsn8SReVMlOtzcNczczqa3QUU9qcdH1EnAT8rnWy1AbkAGFmVlejNYiIqAH6SOrSSvlpE0E761QxM2sBWZ6DWAj8U9I00vmYACLiqpLlqtW5BmFmVleWAFEF3E8yXXfWeZi2Pu6kNjPbTJapNr7XGhlpW34OwsysriYDhKQpBZJXArOAa/OfkdhahdeDMDOrJ8sw18Ukk+39KX2tB5YDnwCuLV3WWpc7qc3MNpelD2LfiPhMbkPSVOChiDhI0vOly1rrCVyDMDOrK0sNYldJ/fO2+wK5pdvWtXyW2oL7IMzM6spSg7gAeFTSiyQtMf8FnCVpW+DmUmautYTcwGRmVleWUUx3SbofGEwSIJ6LiDXp7isaOk/SIJJ5m3I+CnwfOBAYlKb1BFZExJAC5y8imeKjBtgQERVNluYDcQ3CzCxfllFM3UnmYBoQEWdK2l3SHk2tKhcR84Eh6TXKgErgzoj4Vd61ryQZEdWQgyPirQzl+IDcxGRmVleWPojr0+M+lW4vAX5a5H1GA69ExGu5hHSd6mOAW4q8Vgm4icnMrK4sAWKPiPgpUA0QEe9T/DfqcdQPBJ8GlkXEggbOCeA+SbMlndHQhSWdIWmWpFlVVVVFZmvTjVyDMDPbXJYAsV5SN9JGekkDSZ6FyERSV+BI4I46u46n8drDyIjYDzgc+KakgwodFBETI6IiIirKy8sLHZIlkx7mamZWR5YA8SPgXqC/pJuAfwAXFnGPw4GnImJZLkFSZ+BLbN6JvZmIWJL+fBO4ExhWxD2L5D4IM7O6soxiulfSbGAESdPS+emXdlaFagqHAC9GxOJCJ6RDaDtFxLvp+8+RBKqSCPdBmJnV02gNQlKZpMOBE0mGqa4F3s56cUk9gEOBuvM51euTkNRX0vR0c1fgEUlPA08Ad0fEvVnvWyyHBzOz+hqsQUjqAzxIEhDmkHyPHgX8QtJnI+I/TV087dDeuUD6yQXSlgBj0/cLgX2zFeGDC7mJycysrsaamH4KTIqIK/MTJX0buBQ4pZQZa01ek9rMrL7GAsSBEVEvCETEL9NpN9oVeRSTmdlmGuuDWNPMfVshuR/CzKyOxmoQO0o6skC6gB1KlJ+2ITcxmZnV1ViA+BdwdAP7/l2CvLSZ8HMQZmb1NBggIuKrrZmRtuUGJjOzurI8Sd1BuAZhZpbPAQJArkOYmdXV4QPE1DmVLFu1nuoNNYy87EGmzqls6yyZmW0Rsiw5iqRhwID84yPizyXKU6uZOqeSC6c8y+WAFFSuWMOFU54FYPzQfm2bOTOzNpZlRbkbSZYbnUuy/CckDfZbfYC4fMZ81lTXEF02pa2pruHyGfMdIMysw8tSgxgODI6IjaXOTGtbsiJ53m8jnTYb5ppLNzPryLL0QTwH9C51RtpC357dgaQ61CkvQOTSzcw6siwBYkfgBUl3S5qSe5U6Y63h/DGD6N6lbLMaRPcuZZw/ZlAb58zMrO1laWK6tOS5aCO5fgbdKToR9OvZnfPHDHL/g5kZ2VaU+3trZKStjB/ajydnbkvnlfCvCZ9t6+yYmW0xGmxikvRQ+vMdScvzXu9IWt7UhSUNkjQ377VK0rmSLpFUmZc+toHzD5M0X9LLkiY0v4hNC3XyXExmZnU0VoM4OP3ZrA7qiJgPDIFk6VKgEriTZKGhX0bEFQ2dmx5/DclypYuBJyXdFRHPNycvWXSi3Q3SMjP7QBqsQeSGtUZETaFXkfcZDbwSEa9lPH4Y8HJELIyI9cCtwLgi75mdOvwD5WZm9bTWN+NxwC1522dJekbS9ZJ6FTi+H/BG3vbiNK0eSWdImiVpVlVVVfNyJ7kGYWZWR8kDhKSuwJHAHWnS74CPkTQ/LQWuLHRagbSCnQQRMTEiKiKiory8vFl5DNcgzMzqyfTNKKm/pIPT99tI2raIexwOPBURywAiYlnaTLURuJakOamuxcCH87b7A0uKuGeRXIMwM6uryQAh6VTgLmBSmvQRYFoR9zievOYlSX3y9n0RmFfgnCeBPSQNTGsgx6V5KA2PYjIzqydLDeIckvmYVgFExEvALlkuLqkHyUik/Cevfy7pWUnPkIyU+nZ6bF9J09N7bADOAmYALwC3R8RzmUrULF5y1MysrixPUq+NiPVS0i2QDkHNtL5ORLwP7FwnreBSphGxBBibtz0dmJ7lPh+YtNlcTGZmlq0G8S9JFwDd0n6I24C/lTZbrSvppHaAMDPLlyVAXAC8C7wIfAv4O3BRKTPV+kSncIAwM8vXaBNT2px0fUScRDI8tX1yJ7WZWT2N1iDSJ6b7SOrS2HFbPQcIM7N6snRSLwT+KWkasDqXGBFXlSxXbcDPQZiZbS5LgKgC7gd6pK92J1SWbViWmVkHkmU9iO+1RkbalECuQZiZbabJACHpfgqMAY2Iz5UkR21BnVyDMDOrI0sT08V577sBRwHrSpOdNqJO7oMwM6sjSxPT43WSHsqtNtdeCFEmj2IyM8uXpYlph7zNTsD+QJ8GDt8q5ab7jo0bUSdP/W1mBtmamJ4j6YMQsAF4FTi9lJlqdek8UxHhvggzs1SWAPHRiKjOT5CU5bytR64G4ek2zMxqZWlPqdsHAfBES2ekLeVqDRs3FrvUtplZ+9VgTUDSLiR9Dd0l7cOm79EdaG8PzKkMcIAwM8vXWFPREcCpJMt9/jYv/V2gXT08F3l9EGZmlmgwQETEDcANko6JiNuLvbCkQSRrR+R8FPg+0A/4ArAeeAU4JSJWFDh/EUkwqgE2RERFsXnIntlNo5jMzCyR5TmI2yWNAfYieVAul/7TJs6bDwyB2mnDK4E7gUHAhRGxQdLPgAuB7zRwmYMj4q0sBfkgVFuDcIAwM8vJ8hzEb4GewEHADSRPUj9W5H1GA9Y/OSUAABHGSURBVK9ExGvAa3npjwFfLvJaLS+tQWx0DcLMrFaWUUyfioivAG+nE/d9kqRfohjHAbcUSD8VuKeBcwK4T9JsSWc0dGFJZ0iaJWlWVVVVkdmqvQgAG90HYWZWK0uAWJv7KelD6faArDeQ1BU4ErijTvpFJA/e3dzAqSMjYj/gcOCbkg4qdFBETIyIioioKC8vz5qtOpl0H4SZWV1ZAsR0ST2BK4C5wCJgchH3OBx4KiKW5RIknQR8HjghGhg6FBFL0p9vkvRdDCvinkVKR/B6mKuZWa2m1qTuBNyTjjK6Q9LfgO4RsbyIexxPXvOSpMNIOqU/ExHvN3DfbYFOEfFu+v5zwI+KuGdx/CS1mVk9Ta1JvRH4dd72mmKCg6QewKHAlLzkq4HtgfslzZX0+/TYvpKmp8fsCjwi6WmSp7bvjoh7s963WLlRTH5QzsxskyxzKt0vaVxETCv24mkNYec6abs3cOwSYGz6fiGwb7H3a7bcKCbXIMzMamUJEGcBO0paB6whabCPiNippDlrTWmAwM9BmJnVyhIgepc8F20t96CcRzGZmdVqchRTRNQARwPfSd/3IX1Cur2QO6nNzOppMkBIuho4GPhqmvQ+8PtSZqrVuZPazKyeLE1MIyJiP0lzACJiefrwW7uh2j4I1yDMzHKyPChXnT4PEQCSdgbaVWN9eC4mM7N6sgSIa4C/AOWSfgg8AvyspLlqZXIntZlZPVk6qf8IXEwy1cZy4OiIuLXUGWstU+dUMu3ppQCcfctsps6pbOMcmZltGbLUIADKgGqSRX6ynrPFmzqnkgunPMt765Kaw/L31nLhlGcdJMzMyDaK6SKSuZT6kkzz/WdJF5Y6Y63h8hnzWVNdw8Z0sr5OBGuqa7h8xvw2zpmZWdvLMorpRGD/3MR6kn4CzAYuLWXGWsOSFWsAiDRAKOmHr003M+vIsjQXvcbmgaQzsLA02WldfXt2B6Am/RjK0sFZuXQzs44sS4B4H3hO0iRJ1wLPAisk/ULSL0qbvdI6f8wguncp2yxAdO9SxvljBrVxzszM2l6WJqa701dOsetRb7HGD+0HwEN/nQ0bYedtO/ONI/apTTcz68iaDBARcV1rZKStjB/ajwHLPwb/hO+PHcR/OTiYmQHZRjEdJulJSW9KWi7pHUlNLhokaVC6IFDutUrSuZJ2knS/pAXpz16N3He+pJclTWhO4bJSpzIANm7cUMrbmJltVbL0QVwNfAPoB5STTP9d3tRJETE/IoZExBBgf5K+jDuBCcDfI2IP4O/p9mYklZE8wX04MBg4XtLgTCVqhk5pgKDGk/WZmeVkCRCLgbkRUR0RNblXkfcZDbwSEa8B44Cb0vSbgPEFjh8GvBwRCyNiPXBrel5JbKpBOECYmeVk6aS+APirpJnAulxiRFxVxH2OI3nYDmDXiFiaXmOppF0KHN8PeCNvezHwyUIXlnQGcAbAbrvtVkSW8q7RKfkY3MRkZrZJlhrED4EaoCdJ01LulUk6NfiRwB1F5EsF0grOxR0REyOiIiIqysszZ2vzm5UlNYio8WR9ZmY5WWoQu0TE/h/gHocDT0XEsnR7maQ+ae2hD/BmgXMWAx/O2+4PLPkAeWhUp065FeXcxGRmlpOlBvF3SZ/9APc4nk3NSwB3ASel708CphU450lgD0kD0xrIcel5JZHrgwh3UpuZ1coSIE4HHpD0XjHDXAEk9QAOBabkJV8GHCppQbrvsvTYvpKmA0TEBuAsYAbwAnB7RDyXtVDF6lSWVKTCndRmZrWyNDH1bu7F0wn+dq6T9jbJqKa6xy4BxuZtTwemN/fexaitQbiT2sysVpYFg2qAo4HvpO/7AENKnbHWlBvFhPsgzMxqZXmS+mrgYOCradL7wO9LmanW1ikdxbTRo5jMzGplaWIaERH7SZoDEBHL047jdiP3JHXS9WFmZpCtk7paUifS5xAk7Qy0qz+1czUIT7VhZrZJgwFCUq52cQ3wF6Bc0g+BR4CftULeWo1HMZmZ1ddYE9MTwH4R8UdJs4FDSJ5wPjoi5rVK7lpJrpPaD8qZmW3SWICone4ifQahZM8htLXaJibXIMzMajUWIMolndfQzojYqpcbzVfmJiYzs3oaCxBlwHYUnjivXck9KOcahJnZJo0FiKUR8aNWy0kb+sdLb3M0MP2ZSs579UHOHzPI61KbWYfX2DDXdl9zAJg6p5LL718AQBkbqVyxhgunPMvUOZVtnDMzs7bVWICoN19Se3T5jPmsqU7eK328Y011DZfPmN+GuTIza3sNBoiIyDRj69ZuyYo1bEg/hs55z/8tWbGmrbJkZrZFyPIkdbvWt2d3NqRdMZ2p2SzdzKwj6/AB4vwxg+jcpQsAXZTMxdS9SxnnjxnUltkyM2tzWSbra9dyo5Wqp5bRhQ3069ndo5jMzChxgJDUE5gE7E0y2d+pwLlA7s/znsCKiKi3voSkRcC7QA2wISIqSpXP8UP7sWZqGfv3355vnvlBVlc1M2s/Sl2D+DVwb0R8OZ0ivEdEHJvbKelKYGUj5x8cEW+VOI8AVKsz1KxvjVuZmW0VShYgJO0AHAScDBAR64H1efsFHANsEX+y19AZbaxu62yYmW0xStlJ/VGgCrhB0hxJkyRtm7f/08CyiFjQwPkB3CdptqQzGrqJpDMkzZI0q6qqqtmZraYz8prUZma1ShkgOgP7Ab+LiKHAamBC3v7jgVsaOX9kROwHHA58U9JBhQ6KiIkRURERFeXl5c3ObA1lrkGYmeUpZYBYDCyOiMfT7ckkASO3GNGXgNsaOjkilqQ/3wTuBIaVMK/UyDUIM7N8JQsQEfEf4A1JuRFLo4Hn0/eHAC9GxOJC50raVtL2uffA54CSLlK0QZ3p5BqEmVmtUo9iOhu4OR3BtBA4JU0/jjrNS5L6ApMiYiywK3Bn0o9NZ+DPEXFvKTNao84oXIMwM8spaYCIiLlAvecXIuLkAmlLgLHp+4XAvqXMW76pcyoZWN2JFetXM/IyT/dtZgaeaoOpcyq5cMqzVJM8Se3pvs3MEh0+QCTTfddQHZ3pomSyPk/3bWbmAFE7rfc6urDNpuf4PN23mXV4HT5A5Kb1XkNXuucFCE/3bWYdXYcPEOePGUT3LmWsYRu6sw7wdN9mZuDpvmtHK62fug3dWe/pvs3MUh0+QEASJB57tDfdl63nXxO2iLkDzczaXIdvYoJkqOszy9bTjXWMvPQBD3E1M8MBovY5iJUbutBZG3lz5Wo/B2FmhgNE7XMQa+kKQHfW+jkIMzMcIGqfd1hFDwB20PubpZuZdVQdPkDknnd4J7YHoBfvAbBj9y5tliczsy1Bhw8Q548ZRJdOYkUki931VBIgVq/f4H4IM+vQOnyAGD+0H9t168w7bF6DqK4J90OYWYfm5yCAd96vZgO9APiQ3q5Nr3Q/hJl1YB2+BgEg4F16sDy2Yze92dbZMTPbIpS0BiGpJzAJ2BsI4FRgDHA6UJUe9t2ImF7g3MOAXwNlJCvNXVaqfEb685Xoy+BOr222b8CEu0t1WzOzFjXyYztx8+kHttj1Sl2D+DVwb0TsSbJC3Atp+i8jYkj6KhQcyoBrgMOBwcDxkgaXOK88vvHjfEILGailpb6VmVmL+9cryznh2kdb7Holq0FI2gE4CDgZICLWA+vTdaabMgx4OV16FEm3AuOA50uR1149uvDO+9XcsuGzfLXsfmZ0vYA36cXa6EqQKb/NEk0fYmbWpHfYnmPXfx9IgkRLKWUT00dJmpFukLQvMBv4VrrvLElfA2YB/xsR79Q5tx/wRt72YuCThW4i6QzgDIDddtutWRn9wRf24tzb5lJJOV9c/yOOKZtJuVawDdXNul4WcngwsxayKh2m39JKGSA6A/sBZ0fE45J+DUwArgZ+TPIH9I+BK0n6JvIV+rO94DdqREwEJgJUVFQ061t3/NB+nHvbXAAWRl8u2/CV5lzGzKxdKWUfxGJgcUQ8nm5PBvaLiGURURMRG4FrSZqTCp374bzt/sCSEuaVE4c3r/ZhZrYlGfmxnVrsWiULEBHxH+ANSbml2UYDz0vqk3fYF4F5BU5/EthD0kBJXYHjgLtKlVeA/zd+HwcJM9uqtfQoJkWUri1c0hCSYa5dgYXAKcBVwBCSJqNFwDciYqmkviTDWcem544FfkUyzPX6iPhJU/erqKiIWbNmlaIoZmbtkqTZEVFRcF8pA0Rrc4AwMytOYwHCT1KbmVlBDhBmZlaQA4SZmRXkAGFmZgW1q05qSVXAa00eWFhv4K0WzM7WwGVu/zpaecFlLtZHIqK80I52FSA+CEmzGurJb69c5vavo5UXXOaW5CYmMzMryAHCzMwKcoDYZGJbZ6ANuMztX0crL7jMLcZ9EGZmVpBrEGZmVpADhJmZFdThA4SkwyTNl/SypAltnZ+WIunDkv4h6QVJz0n6Vpq+k6T7JS1If/bKO+fC9HOYL2lM2+W++SSVSZoj6W/pdrsuL4CknpImS3ox/X0f2J7LLenb6b/peZJukdStPZZX0vWS3pQ0Ly+t6HJK2l/Ss+m+q5Rx3WcAIqLDvkimEn+FZHnUrsDTwOC2zlcLla0PyQJNANsDLwGDgZ8DE9L0CcDP0veD0/JvAwxMP5eyti5HM8p9HvBn4G/pdrsub1qWm4DT0vddgZ7ttdwkyxG/CnRPt28nWfe+3ZUXOIhkVc55eWlFlxN4AjiQZKXOe4DDs+aho9cghgEvR8TCiFgP3AqMa+M8tYiIWBoRT6Xv3wVeIPmfaxzJFwrpz/Hp+3HArRGxLiJeBV6m8Gp/WyxJ/YEjSNYgyWm35QWQtAPJF8l1ABGxPiJW0L7L3RnoLqkz0INktcl2V96IeBhYXie5qHKmC7TtEBGPRhIt/ph3TpM6eoDoB7yRt704TWtXJA0AhgKPA7tGxFJIggiwS3pYe/gsfgVcAGzMS2vP5YWk9lsF3JA2rU2StC3ttNwRUQlcAbwOLAVWRsR9tNPyFlBsOful7+umZ9LRA0Shtrh2Ne5X0nbAX4BzI2JVY4cWSNtqPgtJnwfejIjZWU8pkLbVlDdPZ5JmiN9FxFBgNUnTQ0O26nKnbe7jSJpR+gLbSjqxsVMKpG015S1CQ+X8QOXv6AFiMfDhvO3+JNXVdkFSF5LgcHNETEmTl+XWBU9/vpmmb+2fxUjgSEmLSJoKPyvp/2i/5c1ZDCyOiMfT7ckkAaO9lvsQ4NWIqIqIamAKMIL2W966ii3n4vR93fRMOnqAeBLYQ9JASV2B44C72jhPLSIdqXAd8EJE/CJv113ASen7k4BpeenHSdpG0kBgD5LOra1CRFwYEf0jYgDJ7/HBiDiRdlrenIj4D/CGpEFp0mjgedpvuV8Hhkvqkf4bH03Sv9Zey1tXUeVMm6HelTQ8/by+lndO09q6p76tX8BYkhE+rwAXtXV+WrBcnyKpSj4DzE1fY4Gdgb8DC9KfO+Wdc1H6OcyniJEOW9oLGMWmUUwdobxDgFnp73oq0Ks9lxv4IfAiMA/4E8nInXZXXuAWkn6WapKawNebU06gIv2sXgGuJp1BI8vLU22YmVlBHb2JyczMGuAAYWZmBTlAmJlZQQ4QZmZWkAOEmZkV5ABhVgRJ79XZPlnS1W2VH7NScoAw2wJIKmvrPJjV5QBh1kIkfUTS3yU9k/7cLU2/UdKX8457L/05SsmaHX8GnpW0raS7JT2drnVwbBsVxQxIJvoys+y6S5qbt70Tm6ZnuRr4Y0TcJOlU4Cqanlp5GLB3RLwq6ShgSUQcASBpxxbOu1lRXIMwK86aiBiSewHfz9t3IMliRZBMAfGpDNd7IpL5+wGeBQ6R9DNJn46IlS2XbbPiOUCYlU5uHpsNpP+vpROmdc07ZnXtwREvAfuTBIpLJeUHH7NW5wBh1nL+TTKTLMAJwCPp+0UkX/yQrGXQpdDJkvoC70fE/5EsirNfyXJqloH7IMxazjnA9ZLOJ1nl7ZQ0/VpgmqQnSGbgXN3A+fsAl0vaSDKD53+XOL9mjfJsrmZmVpCbmMzMrCAHCDMzK8gBwszMCnKAMDOzghwgzMysIAcIMzMryAHCzMwK+v8bxCwaWK1ZDwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(t,Temp_numerical,'o-',label='100'+' Euler steps')\n",
    "plt.plot(t,Temp_analytical,label='analytical')\n",
    "plt.title('Body Temperature vs Time')\n",
    "plt.xlabel('Hours')\n",
    "plt.ylabel('Temperature in Degrees F')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As t approaches infinity the body temperature reaches the ambient tempertaure because it cannot fall further than that, since the heat has nowhere to go."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the time of death is 0.8489352983033802 hours before the body was initially found at 85 degrees F at 11am\n",
      "Therefore the time of detah is about 10:09am\n"
     ]
    }
   ],
   "source": [
    "#c\n",
    "#Using the equation outlined in the problem 3 description\n",
    "#we assume the body was found at 85 degrees farenheit\n",
    "#and we assume T(0) to be the initial 98.6 degrees\n",
    "#from there we can solve for t and subtract that\n",
    "#from the time the body was actually found to find the time \n",
    "#of death\n",
    "T_found = 85\n",
    "T_0 = 98.6\n",
    "K = 0.61111111112\n",
    "T_ambient = 65\n",
    "t = np.log((T_found-T_ambient)/(T_0 - T_ambient))/(-K)\n",
    "print('the time of death is',t,'hours before the body was initially found at 85 degrees F at 11am')\n",
    "print('Therefore the time of detah is about 10:09am')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Now that we have a working numerical model, we can look at the results if the\n",
    "ambient temperature is not constant i.e. T_a=f(t). We can use the weather to improve our estimate for time of death. Consider the following Temperature for the day in question. \n",
    "\n",
    "    |time| Temp ($^o$F)|\n",
    "    |---|---|\n",
    "    |8am|55|\n",
    "    |9am|58|\n",
    "    |10am|60|\n",
    "    |11am|65|\n",
    "    |noon|66|\n",
    "    |1pm|67|\n",
    "\n",
    "    a. Create a function that returns the current temperature based upon the time (0 hours=11am, 65$^{o}$F) \n",
    "    *Plot the function $T_a$ vs time. Does it look correct? Is there a better way to get $T_a(t)$?\n",
    "\n",
    "    b. Modify the Euler approximation solution to account for changes in temperature at each hour. \n",
    "    Compare the new nonlinear Euler approximation to the linear analytical model. \n",
    "    At what time was the corpse 98.6$^{o}$F? i.e. what was the time of death? \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [],
   "source": [
    "#a\n",
    "def ambient_temp(t):\n",
    "    i=t+3\n",
    "    Temp_ambient = np.array([55,58,60,65,66,67])\n",
    "    print(Temp_ambient[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "65\n"
     ]
    }
   ],
   "source": [
    "ambient_temp(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "66\n"
     ]
    }
   ],
   "source": [
    "ambient_temp(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "67\n"
     ]
    }
   ],
   "source": [
    "ambient_temp(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "55\n"
     ]
    }
   ],
   "source": [
    "ambient_temp(-3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f975831e518>"
      ]
     },
     "execution_count": 168,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUZfbA8e9JKIFQAoTeAwFEagQEDAqigOjPLmBZsazormLXxbLu6q6IbV0Vd13FvoAgAroqRbFRBCQk9NBb6KEGSEg7vz/uTRxCyiRkMsnM+TxPnszcmXvveWeSM++8997ziqpijDEmeIT4OwBjjDFlyxK/McYEGUv8xhgTZCzxG2NMkLHEb4wxQcYSvzHGBBlL/KZIIrJNRC4p4LF+IrK+rGMy5Z+IbBaRPv6Ow5zJEn8AE5EfReSwiFT11T5Udb6qti+NbRXxAXOziBx3f1JFJNvj/vHS2L+/iMgQEdnk7ziKQ0TuLOT9SAZQ1Taq+ou/YzVnssQfoESkFdAPUOBKvwZTClR1oqrWUNUawGXA7pz77rJySURCRMSn/2ciUsmX28+Pqr7n8dpfA2zxeD8iyzoeUzyW+APXrcBi4ENgpOcDIvKhiPxLRGa5PbSFItJIRP7pfkNIFJHuebbXU0TWuo9/ICJh7rb6i0iSx7abiMjnInJARLaKyP0ej/1VRKaKyMcikiIia0Skh/vYJ0AL4H9uTI8Xt8Ei0lxEvhCRZBHZIiL3eDw2TkQmisgUd/sJItJaRP7iPn+biAzweP5iEfmbiMSJyFG3TbU9Hu8nIktE5IiILBeRC/Ks+5yILAFOAk1E5G73dU0RkU0icof73HrADCDKo8dcT0Q+FZGnPbZ52rcCEdkrIo+KyBrgWFHtz/M69ReRHSIiHstuFJGl7u0LRCReRI65+3mhuO+FR4yxJXz967p/J3tFZKf7PMtXpUVV7ScAf4BNwB+B84AMoKHHYx8Cye5jYcD3wFacD4tQ4O/ADx7P3wasBpoDdYGFwN/dx/oDSe7tECAOeAaoAkQBW4DB7uN/BdKAoe5+XgAW59nPJV60LXefHstCgVXAn9x9twN2ABe5j4/DScIDgErAFLfNj7r3RwPrPLa3GNgOdABqAP8DJriPtQIOApe4bR4KHADqeKy7BWgPVHa3fyXQGhB3vVTgXPf5Q4BNedrzKfC0x/3TngPsBX4FmgDVimp/nm0LsBPo57Hsf8CD7u144Ab3dk3g/CLejzPi94gxtoSv/yzgTaA60NiNaaS//68C5cc+QQOQ28tqCUxV1ThgM3BTnqfNUNU4VU3D6XGmqerHqpqF80+Zt8c/XlV3quoh4Hngxnx23ROor6rPqWq6qm4B3gVGeDxngap+4+7nE6DrWTY3RywQpqovuvveAHyQZ9/zVPUHVc0EpgG1gFfd+58CHUSkmsfzP1DVRFU9DvzFo80jgemq+p2qZqvqN8BaYJDHuhNUdb2qZqhqpqp+qapb1fEd8JMb89l4TVV3q2qql+0HQJ3MOiWnPSJSF+fDaIr7lAygnYjUU9UUVV1ylnHm8Or1F5GWwIXAw6p6UlX3AG/k1xZTMmU+NmjKxEhgrqomu/cnucte83jOPo/bqfnczztuvtPj9nacnmZeLXGGNY54LAsF5nvc3+tx+yQQJiKV3H/+s9ESaJXPvr/zuJ+3jQfcJJhzHyDc43beNld3h3taAjeKyA0ej1fm9NfEc11E5ErgaaAtzreE6pz+upSE5z68ab+nScAcdyjuBpwP5D3uYyNxvp1tcIeXnlHVOWcZK3j/+rfE+SZ6wGM0KgTnW6wpBZb4A4zbYx0GhIpITpKtCkSISFdVXVHCTTf3uN0C2J3Pc3YCW1U1uoT7OJtSsTuBRFXtfBbbyCtvm0+q6lER2YnTox9dyLq5bRGRcOAz4HpglqpmishsnCGX057r4QTOh0OORoXtg2K2X1WXi0jOcNVNOMN/OY+tA4aLSChOL3u6iNRR1XRvtl0KdgLHcYbOrHywD9hQT+C5GsgCOgLd3J9zcHqXt57Fdu8VkWbusMCT/DYs4GkpcExE/uR+ZQ8VkU4i0tPLfezDOS5QEgsARORBEQkTkUoi0kVEYkq4PYDbRKSdiNTA6QHntPkj4AYRGei2sZp7O7/kDM4YfGVgP5Dt9v77ezy+D2jg7idHAnCFiESISFOcMfDClKT9k3HG2HsC03MWisit7jBPFnAU5wMmu4j9lxpV3YpznOQlEakpzplR0TkHis3Zs8QfeEbijE3vUNW9OT/AeOBmKfmpf5OAuTgHLbfgHAA+jZso/g/nw2YrzgHkCUDtvM8twAvA0+6ZMo8WJzhVzcA5yNoXZ1jmAPBvzhyyKo5PcJLjLpzE94i7ry3AdcCzOG3cDjxAAf9P7pDbozgHUA/ifDh/4/GUFcCXwHa37XWB93GGNnYAX7lxFKiE7Z8EDMT5FnLUY/kVwHoRScF5T4aVwlBccd0IRACJwCGcD92GZRxDwBL7JmXMmURkMc4B7f/6OxZjSpv1+I0xJshY4jfGmCBjQz3GGBNkrMdvjDFBpkKcxx8ZGamtWrXydxjGGFOhxMXFJatq/bzLK0Tib9WqFcuWLfN3GMYYU6GIyPb8lttQjzHGBBlL/MYYE2Qs8RtjTJCpEGP8+cnIyCApKYm0tDR/h2IqqLCwMJo1a0blypX9HYoxZarCJv6kpCRq1qxJq1at8CjdaoxXVJWDBw+SlJRE69at/R2OMWWqwib+tLQ0S/qmxESEevXqceDAAX+HYky+Zsbv4uU569l9JJUmEdV4bHB7ru7etFS2XWETP2BJ35wV+/sx5dXM+F08MX0VqRlZAOw6ksoT01cBlEryt4O7xhhTzrw0OzE36edIzcji5TnrS2X7lvjP0owZMxAREhMTi73ubbfdxrRp085YvmzZMu6///4SxzR27Nh8l59//vl069aNFi1aUL9+fbp160a3bt3Ytm1biffla9OnTy/Ra2tMRaKqJO49xoT5W7j1/aXsPpr/SSu7j6Tmu7y4KvRQT3H4arxs8uTJxMbG8umnn/LXv/717AMFevToQY8ePUq8/tixY3nyySfPWL5kiTNn9ocffsiyZcsYP358ifdRmjIzM6lUKf8/xenTpxMSEkKHDh1KZXvGlBf7j6WxYFMyCzYmM39TMgdSTgEQVT+c8KqhnDiVdcY6TSKqlcq+g6LHnzNetutIKspv42Uz43ed1XaPHz/OwoULee+99/j0009zl//4449cdNFFDBs2jHbt2jFmzBgmTpxIr1696Ny5M5s3b8597nfffUe/fv1o164dX331Ve76V1xxBQAnTpzgjjvuoGfPnnTv3p0vvvgCcJL3tddey5AhQ4iOjubxxx8HYMyYMaSmptKtWzduvvlmr9sya9Ys+vTpQ0xMDMOHD+fEiRMANGvWjKeeeorevXvTs2dPli9fzqBBg2jTpg3vvvtubhsGDBjA1VdfTceOHbn33nvJqfpa2Hb/9re/ccEFFzBjxgzefvttevbsSdeuXbnhhhtITU1l/vz5fPPNNzz00EO530xiY2NJSEgAYO/evbRt2xaACRMmMGLECK644gouu+wyAMaNG0evXr3o0qULzz33XHHeWmNKXWp6Fj+u38/fv1rL4Nd+ptfYeTw8dQU/rN9P76h6vHRdFxaNuZjvH+nP81d3plrl0NPWr1Y5lMcGty+VWAKiW/Ts/9awdvexAh+P33GE9KzTpwxNzcji8Wkrmbx0R77rdGxSi7/837mF7nfmzJkMGTKEdu3aUbduXZYvX05MjDPF6YoVK1i3bh1169YlKiqK3//+9yxdupTXX3+dN998k3/+858AbNu2jZ9++onNmzczYMAANm3adNo+nn/+eS6++GLef/99jhw5Qq9evbjkkksASEhIID4+nqpVq9K+fXtGjx7NuHHjGD9+fG5y9Mb+/fsZN24c8+bNo3r16jz//PO8/vrrud8aWrVqxeLFixk9ejR33nknCxYs4Pjx43Tt2pW77roLcL5NrF27lubNm3PppZfyxRdf0Ldv30K3Gx4ezsKFCwE4ePAg99xzD+B8eH344Yf84Q9/YOjQoVx//fVcffXVRbbjl19+ISEhgTp16vDNN9+wY8cOlixZgqoydOhQFi1aRN++fb1+XYw5G9nZyprdx5i/6QALNiazbNth0rOyqRIaQs/Wdbi6ewf6RUfSsXEtQkJOP9EgZzTCzuo5C3mTflHLvTV58mQefPBBAEaMGMHkyZNzE3/Pnj1p3LgxAG3atGHQoEEAdO7cmR9++CF3G8OGDSMkJITo6GiioqLOGM+eO3cuX375Ja+88grgnMa6Y4fzYTVw4EBq13ams+3YsSPbt2+nefPmxW7HokWLWLt2bW5STE9PJzb2t3mtr7zyytzYMzMzCQ8PJzw8nJCQEI4fPw5A7969yamgOmLECBYsWABQ6HaHDx+ee3vlypU888wzHDlyhJSUlNxvPMUxaNAg6tSpAziv26xZs+jevTvgfDvbsGGDJX7jU7uOpLJg4wHmb0xm4aZkDp/MAKBDo5qM7NuS2Oj69GpVl2pVQovYkpP8SyvR5xUQib+onvkF475nVz4HRZpGVGPK3X1KtM+DBw/y/fffs3r1akSErKwsRISXXnoJgKpVq+Y+NyQkJPd+SEgImZm/zVud95TCvPdVlc8//5z27U//irdkyZLT9hEaGnradotDVRkyZAiffPJJvo97xp63XTn7zK8dRW03PDw89/att97KrFmz6NSpExMmTGDx4sX5rlOpUiWys50P7LxXbXtuT1V5+umnufPOO/PdjjGlISUtg8VbDuUm+y3JzlBm/ZpVGdChAf2iI7mgbSQNaob5OdLTBcUY/2OD25f6eNm0adO49dZb2b59O9u2bWPnzp20bt06t6frrc8++4zs7Gw2b97Mli1bzkjwgwcP5s0338wdM4+Pjy9ym5UrVyYjI8PrGPr27ctPP/3Eli1bAOe4wsaNG4vRCli8eDE7duwgKyuLqVOnEhsbW6ztnjhxgkaNGpGRkcGkSZNyl9esWZOUlJTc+61atSIuLg4g3zOicgwePJj33nsv95hCUlISycnJxWqTMXllZmUTt/0wr3+3kRveXkT3577lro+XMWXZTprXrc7Tl5/DnAcvZOmTA/nHsG5c071ZuUv6ECA9/qL4Yrxs8uTJjBkz5rRl1113HZMmTTptCKMo7du356KLLmLfvn28/fbbhIWd/kfy5z//mQcffJAuXbqgqrRq1Sr3IHBBRo0aRZcuXYiJiWHixIlFxtCwYUPee+89hg8fTnp6OuCcGRQdHe11O/r27csjjzzCmjVr6N+/P1deeSUi4vV2n3vuOXr16kWLFi3o1KlTbm/+xhtv5O677+bVV19l5syZPPbYYwwfPpwPPviAAQMGFBjP0KFDSUxMpHfv3oDzATJp0iQiIyO9bpMxqsr2gyeZvymZBRsPsGjzQVLSMhGBTk1qM+rCKGKjIzmvZR2qVip6+Ka8qBBz7vbo0UPzTsSybt06zjnnHD9FZDx99913jB8/npkzZ/o7lGKzvyOT15GT6SzafJD5G5NZsOkAOw85w8RNI6rRLzqS2OhI+raJpG54FT9HWjQRiVPVM84ND4oevzHGFCQ9M5vlOw7nnk+/KukI2Qo1qlaid1Q97uoXRWzbSFpHhgdMmQ9L/OasXXLJJbmnmBpT3qkqm/Yfd3v0ySzecpCT6VmEhghdm9Vm9MXR9IuOpGvzCCqHBuZh0Aqd+FU1YD6BTdmrCMOcpnQkHz/Fwk3JTrLfmMzeY84xpFb1qnNdTDNioyPp06YetcKCY26GCpv4w8LCOHjwIPXq1bPkb4otpx5/3oPpJjCkZWTx67ZDzvDNxmTW7nEu8KxdrTKxbZ1x+ti2kTSvW93PkfpHhU38zZo1IykpyeqpmxLLmYHLVHzZ2cq6vcdY4A7fLN16iFOZ2VQOFc5rWYfHBrcntm0knZrWJjTEOoo+TfwiEgFMADoBCtyhqr+IyGjgPiAT+FpVHy/utitXrmwzJxkTJPIrstg7qh7zNx5gwSbnKtnk484pw9ENanDz+S3pFx1Jr9Z1Ca9aYfu3PuPrV+R1YLaqXi8iVYDqIjIAuArooqqnRKSBj2MwxlRg+U1K8tCUBHKO0ETWqOIO39Qntm0kjWrb8F1RfJb4RaQWcCFwG4CqpgPpIvIHYJyqnnKX7/dVDMaYiisrW1m16yjPfLH6jElJFKgVVolPR/WhQ6OaZxQ5M4XzZY8/CjgAfCAiXYE44AGgHdBPRJ4H0oBHVfXXvCuLyChgFECLFi18GKYxprzYeehk7oVTCzcd5GhqwaVHUtIy6dikVhlGFzh8mfgrATHAaFVdIiKvA2Pc5XWA3kBPYKqIRGmec+tU9R3gHXCu3PVhnMYYPzmamsEvmw+ywC1dvO3gSQAa1QpjUMeGxEZH8sKsRPbmMyNVaU1KEox8mfiTgCRVXeLen4aT+JOA6W6iXyoi2UAkzrcDY0wAy8jKJmHnEfd8+gOsSDpKVrZSvUoovaPqMbJvK/pFR9Kmfo3c07RVOW2MH0p3UpJg5LPEr6p7RWSniLRX1fXAQGAtsBm4GPhRRNoBVQArm2hMAFJVtiafYL57Pv3iLQc5fiqTEIHOzSL4Y/82xLaNpHuLOlSplP9Vsr6elCQY+fqsntHARPeMni3A7cAJ4H0RWQ2kAyPzDvMYYyquQyfSWejOJbtgU3LuXBjN61bjym5N6NfWKXJWu7r3V8n6clKSYOTTxK+qCUB+s4bf4sv9GmPKzqnMLOK2HXZLFyezevdRVKFmWCUuaBPJH/q3oV90JC3rhRe9MVMm7MoGY0yxqCob9h1nvjvr1NKth0jNyKJSiBDTog4PXdKO2OhIujStTaUALXJW0VniN8YUaX9KmjN04w7f7E85BUBU/XCG92xObNtIerepRw27SrZCsHfJGHOG1PQslm47xPwNTkmExL3O9Jd1qlcmNro+/dxCZ3ZKZcVkid8YQ3a2snbPMX7e6JxPv2zbYdKzsqkSGkLP1nX405AO9IuOpGPjWnaVbACwxG9MkNp1JJUF7jj9os0HOXTCKXLWoVFNRvZtSWx0fXq1qku1KhVnLlnjHUv8xgSJlLQMFm855CT7TclsOXACgAY1q9K/fX36RUdyQdtIGtS0ImeBzhK/MQEqMyubFUlH3QOyB4jfcYTMbCWscgi9o+pxU68W9IuuT7uGNWwyoyBjid+YCia/2vQ5FzdtP3iCn91yCIs2HyQlLRMR6Ny0NqMujCI2OpLzWtahaiUbvglmlviNqUDyq03/+LQVTPl1B0lHUtl5yLlKtmlENS7v3JjY6EguaBNJnfAq/gzblDOW+I2pQF6es/6M2vTpWcriLYe4pGND7uoXRWzbSFpHhtvwjSmQJX5jKpDdbt2b/Lx7a37VUYw5k11PbUwFUr9m1XyX24VUpjgKTPzuFIk5tzuUTTjGmIKkZWSR37VTVpveFFdhPf67PG5P8nUgxpjCPfu/tew9doq7L4qiaUQ1BOcg7gvXdraSxaZYvB3jt6NExvjRVyt3M3npDu65qA1jLuvAE5ed4++QTAVWWOKPEJH/w/lWUEtErvR8UFW/9GlkxhgAdhw8yROfryKmRQSPDGrn73BMACgs8S8Ehrm3FwE3eDymgCV+Y3wsPTOb+yYvRwTeuLE7la2+vSkFBSZ+Vf1dWQZijDnTi7MTWZl0lLdviaFZner+DscECOs+GFNOzVu3j/cWbGVkn5YM6dTY3+GYAGKJ35hyaM/RVB75bAUdG9fiiaF2INeULkv8xpQzmVnZPDA5gfTMbMbf1J2wylZQzZSuIhO/iFwrIjXd22NEZKqIdPN9aMYEpzfmbWTptkM8f00nourX8Hc4JgB50+P/q6qmiEhf4P+AKcDbvg3LmOC0cFMyb/6wievPa8Y13Zv5OxwToLxJ/DmlAK8A/qWqnwP5FwwxxpTYgZRTPDglgajIcJ676lx/h2MCmDeJf4+IvAUMB74RkSperoeIRIjINBFJFJF1ItLH47FHRURFJLJkoRsTOLKzlYenJnAsNYO3bo6hehUrnGt8x5sEPgz4CbhcVQ8DkcAYL7f/OjBbVTsAXYF1ACLSHLgU2FHsiI0JQP/5eQvzNybzzP91pEOjWv4OxwS4IhO/qh7HSdC93EWngDVFrScitYALgffc7aSr6hH34deAx3GuADYmqMVtP8Qrc9dzeefG3NSrhb/DMUHAm7N6ngb+AjztLgrDu2qdUcAB4AMRiReRCSIS7tb82aWqK4rY7ygRWSYiyw4cOODF7oypeI6cTOf+yQk0iQjjhes626xZpkx4M9RzPTAUOAGgqrsAb76LVgJigH+rand3/b8CTwHPFLWyqr6jqj1UtUf9+vW92J0xFYuq8vi0lexPSWP8jTHUCqvs75BMkPAm8Z9SVcUdlhERbwuGJAFJqrrEvT8N54OgNbBCRLYBzYDlItKoWFEbEwA+/mU7c9fu409DOtC1eYS/wzFBxJvEP909q6e2iNwOzAXeL2olVd0L7BSRnKmBBgLLVbWBqrZS1VY4Hw4x7nONCRqrdx3l+a/XcXGHBtwZ29rf4ZggU+Q5Y6r6oohcBqTjnJnzvKrO8nL7o4GJ7imgW4DbSxypMQHi+KlMRk+Op254FV65oauN65sy5+3JwquANFX9QUTCRCRcVU8UtZKqJgA9Cnm8lZf7NyYgqCpPz1jF9oMnmHxXb+qGV/F3SCYIeXNWzx04k65McBe1AL7wZVDGBKrP4pKYmbCbBy9px/lR9fwdjglS3ozx3w/0Bo4BqOoGoIEvgzImEG3an8JfvlhD3zb1uHdAW3+HY4KYN4k/TVXTc+6ISCg2+boxxZKWkcW9E+OpXiWUfw7vRmiI/QsZ//Em8S8UkceBMBEZgFOd8yvfhmVMYHnuq7Ws35fCq8O60qBWmL/DMUHOm8T/OJACJAIPAPNwLsIyxnjhq5W7mbRkB3dfFEX/9jZKavyv0LN63GGd91V1JPDvsgnJmMCx4+BJnvh8Fd1bRPDooPZFr2BMGSi0x6+qWUBjEbFryY0ppvTMbO6bvBwRePPG7lQOtZlOTfngzXn8W4D5IvIFbr0eAFV9w2dRGRMAXpqdyMqko7x9SwzN6nhb6cQY3/Mm8R8AvgWquz/GmCLMW7ePCQu2cmuflgzp1Njf4RhzGm9KNvy5LAIxJlDsOZrKI5+toGPjWjw59Bx/h2PMGYpM/CIyPZ/FR4FlwLue5/gbE+wys7J5YHIC6ZnZjL+pO2GVQ/0dkjFn8OZoUxKQCXzi/qQDh4AuwLu+C82YiueNeRtZuu0Qz1/Tiaj6NfwdjjH58maMv6uqXpRzR0RmAj+p6oUistZ3oRlTsSzalMybP2zi+vOacU33Zv4Ox5gCedPjbyginn/FTYCcKbFOlX5IxlQ8ycdP8cCUBKIiw3nuqnP9HY4xhfKmx/848IuIJOLU6GkH3Cci4cBEXwZnTEWQna08NCWBo6kZfHxHL6pX8bbauTH+4c1ZPV+KyLdAR5zEv0ZVU92HX/FlcMZUBP/5eQvzNybz/DWdOKexN9NRG+Nf3tTjr4ZTo+cuVV0GNHVn5DIm6MVtP8wrc9dzeefG3NSrhb/DMcYr3ozxv+8+L9a9vxsY67OIjKkgjp7M4P7J8TSJCOOF6zrbFIqmwvAm8Uer6lggA0BVT2L1+E2QU1Ue/3wF+46l8eaNMdQKs3JWpuLwJvGni0gYoAAi0hrnXH5jgtbHv2xnzpp9/GlIB7o1j/B3OMYUizenHzwHzAaaichHwEXAnT6NyphybPWuozz/9Tou7tCAO2Nb+zscY4rNm7N6ZotIHNAXZ4jnMVXd7/PIjCmHjp/KZPTkeOqGV+GVG7oSYlMomgrIm4lYBgEd3EXrgIO+DsqY8khV+fPM1Ww/eILJd/WmbngVf4dkTIkUOMYvIo2B1TjTLEYBbYCngVUi0sibjYtIhIhME5FEEVknIn1E5GX3/koRmSEiNkBqKoRpcUnMiN/FAwPbcX5UPX+HY0yJFXZwdywwQVVjVXW0qt6nqrE4hdle8HL7rwOzVbUD0BXnG8O3QCdV7QJsAJ4oefjGlI1N+1N45os19Imqx30Xt/V3OMaclcKGevqo6u15F6rqa275hkKJSC3gQuA2d710nLOB5no8bTFwfXECNqaspWVkce/EeKpXCeWfI7oRauP6poIrrMefWsLHckThzN71gYjEi8gEt76PpzuAWV5syxi/ee6rtazfl8Krw7rSsFaYv8Mx5qwV1uOvLSJX5rNcAG8KklQCYoDRqrpERF4HxgB/BhCRp3Dq/Odb6E1ERgGjAFq0sEvhjX98vXIPk5bs4O6LoujfvoG/wzGmVBSW+BcCNxTw2CIvtp0EJKnqEvf+NJzEj4iMBK4ABqqq5reyqr4DvAPQo0ePfJ9jjC/tOHiSMZ+vpHuLCB4d1N7f4RhTagpM/Kr6u7PZsKruFZGdItJeVdcDA4G1IjIE+BNwkVv+wZhyJz0zm/smL0cE3hjRncqh3lzkbkzF4OvC4aOBiSJSBdgC3A78ClQFvnWLWi1W1Xt8HIcxxfLS7ERWJh3l7VtiaF63ur/DMaZU+TTxq2oC0CPPYjsXzpRr3yfuY8KCrdzapyVDOjX2dzjGlDr7/mqMhz1HU3lk6go6Nq7Fk0PP8Xc4xviEVz1+EekFtPJ8vqpO8lFMxvhFZlY2D0xO4FRmNuNv6k5Y5VB/h2SMTxSZ+EXkQ5xpFxOALHexApb4TUB5Y95Glm47xD+GdSWqfg1/h2OMz3jT4+8NdFTVbF8HY4y/LNqUzJs/bOL685pxbUwzf4djjE95M8a/Boj0dSDG+Evy8VM8MCWBqMhwnrvqXH+HY4zPedPjrw2sE5HFwKmchap6rc+iMqaMZGcrD09dwdHUDD6+oxfVq/j6DGdj/M+bv3JvK3EaU+G8M38LP284wN+v7sQ5jb2pRGJMxefNDFzzyiIQY8pa3PbDvDxnPUM7N+Lm860elAkeBSZ+EflJVS8SkcO4E63nPASoqtb1eXTG+MjRkxncPzmeJhFhvHBtF9yryI0JCoX1+Ae4v+3ArgkoqgwnuSsAABzZSURBVMrjn69g37E0pv2hL7WrVfZ3SMaUqcKKtGW7v7MKeo4xFdEni7czZ80+nhp6Dt2a28yfJvhYyQYTVNbsPsrfv1rHgPb1uTO2tb/DMcYvLPGboHH8VCb3TYqnTnhlXh3WjRCbQtEEKa8Sv4g0E5EB7u2q+UyhaEy5pqr8eeZqth88wesjulM3vIq/QzLGb4pM/CJyB/AlMMFd1BL4wpdBGVPapsUlMSN+Fw8MbEfvqHr+DscYv/Kmx38/Tr2eYwCqugGwyUdNhbFpfwrPfLGGPlH1uO9imw7CGG8Sf5qqpufcEZFQnHP5jSn30jKyuG9SPNWrhPLPEd0ItXF9Y7xK/AtF5HEgzB3nnwJ85duwjCkdz321lsS9Kbw6rCsNa4X5OxxjygVvEv/jQAqQCDwAzAOe8mVQxpSGr1fuYdKSHdx9URT929vopDE5Cq3V4w7rvK+qI4F/l01Ixpy9HQdPMubzlXRvEcGjg9r7OxxjypVCe/zuVbuNRcSuaTcVRnpmNqMnL0cE3hjRncqhdrmKMZ68Kcu8BZgvIl8AJ3IWquobPovKmLPw8pxEViQd5e1bYmhet7q/wzGm3PEm8R8AvgWquz/GlFvfJ+7j3flb+V3vlgzp1Njf4RhTLnlTj//PZRGIMWdrz9FUHpm6gnMa1+Kpy8/xdzjGlFtFJn4R+ZbT6/EDoKqDvFg3AueK307uNu4A1uOcEtoK2AYMU9XDxQnamLwys7J54NMETmVm89ZN3QmrHOrvkIwpt7wZ6nna43YYcB0ec+8W4XVgtqpeLyJVcIaKngTmqeo4ERkDjAH+VIyYjck1M34XL89Zz64jqQDcfH5zourX8HNUxpRv3gz1LMmz6CcR+amo9USkFnAhcJu7nXQgXUSuAvq7T/sI+BFL/KYEZsbv4onpq0jN+G3KiOnLd9OzVT2u7t7Uj5EZU755U6StlsdPhIgMBLw5ahaFc2D4AxGJF5EJblXPhqq6B8D9bVfWmBJ5cXbiaUkfIDUji5fnrPdTRMZUDN4M9azBGZ8XIBPYCtzl5bZjgNGqukREXscZ1vGKiIwCRgG0aGETYZvfpGVk8dGibew5mpbv47vdYR9jTP68SfxRqprhuUBEvFkvCUjyGCqahpP494lIY1XdIyKNgf35rayq7wDvAPTo0eOMg8sm+GRlKzPid/GPuevZfTSNqpVCOJWZfcbzmkRU80N0xlQc3lzSmHeMH2BpUSup6l5gp4jkXC8/EFiLU9t/pLtsJFbb3xRBVfkhcT+XvzGfRz9bQf2aVZl8V29evK4L1fKcvVOtciiPDbYSDcYUpsCeu4g0wBnLryYinfmtFHMtvL+QazQw0T2jZwtwO86HzVQRuRPYAdxQwthNEEjYeYQXvlnHkq2HaFWvOm/dFMPQzo0Q+a288stz1rP7SCpNIqrx2OD2dmDXmCKIav6jKCJyO855992ABI+HUoAPVPUz34fn6NGjhy5btqysdmfKga3JJ3h5TiLfrNpLZI0qPDAwmhG9WljdHWOKQUTiVLVH3uUF9vhV9QOcM3KGqepUn0ZnjOtAyilen7eBT5fupEqlEB68JJrf94uiRlVvDisZY7zhzXn8U0VkMHAuzgVcOcvH+jIwE1yOn8rknZ+3MGH+FtIzs7mxVwvuHxhN/ZpV/R2aMQHHm5IN/wIicC7G+gDnyt3FPo7LBIn0zGw+/XUHb8zbSPLxdC7v3JhHB7endWS4v0MzJmB58/05VlW7iMgKVf2ziLwEfO7rwExgU1W+XrWHl+esZ/vBk/SOqsuEkefQrXmEv0MzJuB5k/hzrpJJE5FGwEGcAmvGlMiizcmMm5XIyqSjdGhUkw9u70n/dvVPO1PHGOM73iT+b9wqm6/gnN2ThVNjx5hiWbfnGONmJfLThgM0qR3Gqzd05eruTQkNsYRvTFkqas7dEGCWqh4BPhORr4BqqnqoTKIzASHp8En+MXcDMxJ2USusMk8O7cCtfVpZ6WRj/KTQxK+q2W6Nnd7u/VTACqEYrxw+kc6/ftzER4u2g8CoC6P440VtqV3dpnA2xp+8Ger5VkSuUlUrrWC8kpaRxQcLt/GvHzdx4lQm18U046FL21kNHWPKCW8S/31AbRE5hdPbF0BVta5PIzMVTla28nlcEv/4dgN7j6UxsEMDHh/SgfaNavo7NGOMB28Sf6TPozAVmqoyb91+XpydyMb9x+nWPILXR3Tj/Kh6/g7NGJMPb67czRKRETjlmceKSDOgIRDn8+hMuRe3/TAvzkpk6bZDREWG8++bYxjSqZGdmmlMOebNlbvjgco4V+6OBU4CbwM9fRuaKc82HzjOy7PXM3vNXiJrVOXvV3dieM/mVkTNmArAm6GevqoaIyLxAKp6yC2zbILQ/mNp/HPeRqb8upOwSiE8fGk77oxtTbgVUTOmwvDmvzXDPZ9fAUSkHnDmtEcmoKWkZbhF1LaSmZ3N73q35L6L2xJZw4qoGVPReJP438KpzVNfRJ4FhgHP+jQqU26kZ2Yzccl23vx+E4dOpPN/XZvw6KB2tKxnRdSMqai8Obj7sYjEAZe4i25Q1dW+Dcv4W3a28r+Vu3ll7np2Hkqlb5t6jLmsA12aWRE1Yyo6bwdmQ4EMnOEeO3oX4BZsTGbc7HWs3nWMcxrX4qM7OnNhdKSdqWNMgPDmrJ6ngJuAGTgXb00SkYmq+oKvgzNla/Wuo7w4O5H5G5NpGlGN14Z35aquTQmxImrGBBRvevy3AOep6kkAEXke5xx+S/wBYuehk7w6dz0zE3YTUb0yT19+Dr/r05KqlayImjGByJvEvz3P8yoBW3wTjilLh06kM/77Tfx38XZCQuCP/dtwT/821AqzImrGBDJvEv9JYI2IzMEZ4x8ELBCRfwCo6sM+jM/4QGp6Fu8v3MrbP27mRHomw3o058FL2tGodljRKxtjKjxvEv/X7k8Om2+3gsrMyuazuCRe+3YD+1NOcWnHhjw+uD3RDa2ImjHBxJvTOd8ri0CM76gqc9fu46XZiWw+cILzWtbhrZtj6NnKCqwaE4y8OatnCPA3oKX7fK/LMovINiAFZ7rGTFXtISLdcGr9hAGZwB9VdWmJW2ByzYzfxctz1rP7SCpNIqrx2OD2NKtTjRdmJRK3/TBt6ofzn9+dx6CODe3UTGOCmDdDPeNxrtZdRclKNQxQ1WSP+y8Bz6rqLBEZ6t7vX4LtGg8z43fxxPRVpGZkAbDrSCoPT00gW6FBzaq8cG1nbjivGZWsiJoxQc+bxJ8EJKhqadXnUaCWe7s2sLuUthvUXp6zPjfp58hWqBVWiR8f60/1KlZEzRjj8CYbPA78T0R+BE7lLFTVN7xYV4G5IqLAf1T1HeBBYI6IvIJzFXDf/FYUkVHAKIAWLVp4savgtvtI/lMhp6RlWtI3xpzGm+/9z+KM0UcA9T1+vHGBqsYAlwH3isiFwB+Ah1S1OfAQkO/BY1V9R1V7qGqP+vW93V1wUlVqVcv/3Hub59YYk5c3XcEGqnpeSTauqrvd3/tFZAbQCxgJPOA+5TNgQkm2bRzpmdn85cvVHE3NIESc4Z0c1SqH8tjg9v4LzhhTLnnT458nIhcXd8MiEi4iNXNu41z4tRpnTP8i92kXAxuLu23jOJByipveXczkpTv5Y/82vHJ9V5pGVEOAphHVeOHazlzdvam/wzTGlDPe9PjvAh4VkZNAOt6fztkQmOGeNlgJmKSqs0XkOPC6iFQC0nDH8U3xrEw6wt2fxHH4ZDpv3NidK7s2AeDa85r5OTJjTHnnTeKPLMmGVXUL0DWf5QuAEg0dGcfM+F386fOVRNaoyrR7+tKpaW1/h2SMqUC8uXI3S0RGAFGqOlZEmuH05uN8Hp05TVa28uLsRN75eQu9WtXlX7fE2NSHxphi8+bK3fFAZeBCYCxO0ba3gZ6+Dc14Onoyg9GfxvPzhgPc0rsFz1xxLlUq2cVYxpji82aop6+qxohIPICqHhKRKj6Oy3jYtD+F33+0jF1HUhl7TWduOt+uazDGlJw3iT9DREJwLsZCROpRstINpgS+W7uPB6ckEFY5hEl39bbCasaYs1Zg4heRSqqaCbwFfA7UF5Fncer2PFtG8QUtVeWtHzbx6rcbOLdJLd75XQ+7GMsYUyoK6/EvBWJU9WMRiQMuwTmV8wZVXV0m0QWpk+mZPPbZSr5etYerujVh3LVdqFbFpkE0xpSOwhJ/bt1eVV0DrPF9OGbnoZPc9fEy1u9L4YnLOjDqwigroWyMKVWFJf76IlLgtIqq+g8fxBPUFm1O5t6Jy8nMVj64rSf92zfwd0jGmABUWOIPBWrg0fM3vqGqfPzLdp77ai2t6lXn3Vt7EFW/hr/DMsYEqMIS/x5Vfa7MIglSpzKzeGbmGqYs28nADg3454hu1AzLv9KmMcaUBq/G+I1v7D+Wxj3/jWP5jiPcN6AtD1/ajpAQe9mNMb5VWOIfWGZRBKEVO50ia0dTM3jrphgu79LY3yEZY4JEgYlfVQ+VZSDBZPryJMZMX0X9GlX5/A996dikVtErGWNMKbE5+cpQZlY242YlMmHBVnpH1eWtm2KoZ0XWjDFlzBJ/GTlyMp3Rk+OZvzGZkX1a8vQVHakcakXWjDFlzxJ/GdiwL4W7Pl7G7iOpjLu2MyN6WZE1Y4z/WOL3sblr9vLQlASqVanEp6N6c15LK7JmjPEvS/w+kp2tvPn9Jl77bgNdmtXmP787j8a1rciaMcb/LPH7wIlTmTwydQWz1+zlmu5NeeHazoRVtiJrxpjywRJ/Kdtx8CSjPlnGhn0pPH35OdwZ29qKrBljyhVL/KVo4aZk7p20HFX46I5e9Iuu7++QjDHmDJb4S4Gq8uGibfz963VERYbz7q09aBUZ7u+wjDEmX5b4z9KpzCyenrGaz+KSuLRjQ14b3o0aVe1lNcaUX5ahzsK+Y2nc/UkcCTuPcP/AaB4cGG1F1owx5Z5PE7+IbANSgCwgU1V7uMtHA/cBmcDXqvq4L+Pwhfgdh7n7kziOn8rk3zfHcFlnK7JmjKkYyqLHP0BVk3PuiMgA4Cqgi6qeEpEKN83UtLgknpy+ioa1q/LxnX3p0MiKrBljKg5/DPX8ARinqqcAVHW/H2IokcysbMZ+k8j7C7fSt0093rophjrhVfwdljHGFIuvq4QpMFdE4kRklLusHdBPRJaIyE8i0jO/FUVklIgsE5FlBw4c8HGYRTt8Ip2RHyzl/YVbuf2CVnx8Ry9L+saYCsnXPf4LVHW3O5zzrYgkuvusA/QGegJTRSRKVdVzRVV9B3gHoEePHoofrd/rFFnbezSNl67vwrAezf0ZjjHGnBWfJn5V3e3+3i8iM4BeQBIw3U30S0UkG4gE/N+tz8fs1Xt4eOoKalStxKd39yamRR1/h2SMMWfFZ0M9IhIuIjVzbgODgNXATOBid3k7oAqQXNB2/CU7W3nt2w3c89/lRDesyf9Gx1rSN8YEBF/2+BsCM9w6NZWASao6W0SqAO+LyGogHRiZd5jH346fyuThKQnMXbuP62Ka8fw1nazImjEmYPgs8avqFqBrPsvTgVt8td+ztf3gCe76eBmbD5zgmSs6cvsFrazImjEmoNiVux4WbHSKrAF8dHsvYqMj/RyRMcaUPkv8OEXW3luwlbHfrCO6QU3eufU8WtazImvGmMAU9Ik/LSOLJ2esYvryXQw+tyH/GNaNcCuyZowJYEGd4fYeTePu/8axYucRHrqkHaMvbmtF1owxAS9oE3/c9sPc8984Tp7K5D+/O4/B5zbyd0jGGFMmgjLxT122k6dnrKZR7TD+e+f5tG9U098hGWNMmQmqxJ+Rlc3zX6/jw0XbiG0byfibuhNR3ertGGOCS9Ak/kMn0rl34nJ+2XKQO2Nb88RlHagU6usadcYYU/4EbOKfGb+Ll+esZ/eRVOrXrEpGVjYn0rN49YauXHdeM3+HZ4wxfhOQiX9m/C6emL6K1IwsAPannALgoUujLekbY4JeQI51vDxnfW7S9zT11yQ/RGOMMeVLQCb+3UdSi7XcGGOCSUAm/iYR1Yq13BhjgklAJv7HBrenWp4yytUqh/LY4PZ+isgYY8qPgDy4e3X3pgC5Z/U0iajGY4Pb5y43xphgFpCJH5zkb4neGGPOFJBDPcYYYwpmid8YY4KMJX5jjAkylviNMSbIWOI3xpggI6rq7xiKJCIHgO0lXD0SSC7FcCoCa3NwsDYHh7Npc0tVrZ93YYVI/GdDRJapag9/x1GWrM3BwdocHHzRZhvqMcaYIGOJ3xhjgkwwJP53/B2AH1ibg4O1OTiUepsDfozfGGPM6YKhx2+MMcaDJX5jjAkyAZ34ReQhEVkjIqtFZLKIhPk7Jl8SkQfctq4RkQf9HY+viMj7IrJfRFZ7LKsrIt+KyEb3dx1/xljaCmjzDe57nS0iAXWKYwHtfVlEEkVkpYjMEJEIf8ZY2gpo89/c9iaIyFwRaVIa+wrYxC8iTYH7gR6q2gkIBUb4NyrfEZFOwF1AL6ArcIWIRPs3Kp/5EBiSZ9kYYJ6qRgPz3PuB5EPObPNq4Frg5zKPxvc+5Mz2fgt0UtUuwAbgibIOysc+5Mw2v6yqXVS1G/AV8Exp7ChgE7+rElBNRCoB1YHdfo7Hl84BFqvqSVXNBH4CrvFzTD6hqj8Dh/Isvgr4yL39EXB1mQblY/m1WVXXqep6P4XkUwW0d677tw2wGGhW5oH5UAFtPuZxNxwolbNxAjbxq+ou4BVgB7AHOKqqc/0blU+tBi4UkXoiUh0YCjT3c0xlqaGq7gFwfzfwczzGt+4AZvk7iLIgIs+LyE7gZqzHXzh3jPcqoDXQBAgXkVv8G5XvqOo64EWcr8OzgRVAZqErGVMBichTOH/bE/0dS1lQ1adUtTlOe+8rjW0GbOIHLgG2quoBVc0ApgN9/RyTT6nqe6oao6oX4nxl3OjvmMrQPhFpDOD+3u/neIwPiMhI4ArgZg2+i5AmAdeVxoYCOfHvAHqLSHUREWAgsM7PMfmUiDRwf7fAOeg32b8RlakvgZHu7ZHAF36MxfiAiAwB/gRcqaon/R1PWchzgsaVQGKpbDeQPzRF5FlgOM7Xwnjg96p6yr9R+Y6IzAfqARnAw6o6z88h+YSITAb645Sr3Qf8BZgJTAVa4Hzo36CqeQ8AV1gFtPkQ8CZQHzgCJKjqYH/FWJoKaO8TQFXgoPu0xap6j18C9IEC2jwUaA9k45Smv8c9fnl2+wrkxG+MMeZMgTzUY4wxJh+W+I0xJshY4jfGmCBjid8YY4KMJX5jjAkylvjLkIioiHzicb+SiBwQka/c+1eKyBj39l9F5FH39nMicol7+0G3JENpx/aWWwFwrYikurcTROR6z/37mojcVpIKhCJytYh09Ljv85gLey9EZKCILHdfwwUi0tYH+79HRG4txvNbeVZ+dJfl/p2VQjy3icj4s9zGj/lVGs27PL+2GO9V8ncAQeYE0ElEqqlqKnApkHtOrqp+iXMh0mlU1bM+x4PAfwGvL2ARkVBVzSrsOap6r/vcVsBXbjXAHNO83VcpuA2n7pDXBfXcInxX41QvXAtnvGa+Uth78W/gKlVdJyJ/BJ7GaVupUdW3S3N7viIilTyKq5U73vx/BBrr8Ze9WcDl7u0b8bi6tqAek4h86Pa878epO/SDiPzgPvZvEVnm1mV/1mOdbSLyjIgsAMaIyHKPx6JFJM7bgHP277HdsSLyi7vfGBGZIyKbReQej3UeE5Ff3Vriz+azzVB3u6tFZJU4cydcD/QAJro95WpuG351n/eOexV2Tg9wrIj8hHs1J/Cyu16bfGJ+1u2BrxKRDu7y+uLU7l8uIv8Rke0iEplPrGe8xvm9F3koUMu9XRv3g0xEeonIIhGJd3+3d5ffJiIzReR/IrJVRO4TkYfd5y0Wkbr5xOX5rfBHEXlRRJaKyAYR6efFW5t3e7m9ahGJFJFtHrFNF5HZ4sx38JLHOre7+/sJuMBj+Yci8g/3tXlRRMLFqTf/q9umq9znVRORT92/kylAtRLEHSYiH7jvbbyIDPCIe7zH874Skf7u7ePifCtcAvQRkXHifNtdKSKvFDeGisZ6/GXvU+AZcYZ3ugDvA179k6rqGyLyMDBAVZPdxU+p6iERCQXmiUgXVV3pPpamqrEAInKJiHRT1QTgdpza3yW1U1X7iMhr7nYuAMKANcDbIjIIiMaZG0CAL0XkQrfsbI5uQFN3rgREJEJVj4jIfcCjqrrMXT5eVZ9zb3+CU6flf+42IlT1IvexaJxvKtPc+3ljTlbVGHF6348Cv8e5MvJ7VX1BnHIAowpob36vcX7vhaffA9+ISCpwDOjtLk8ELlTVTHGGosbyW/2VTkB397XcBPxJVbu7r/OtwD8LiC9HJVXtJSJD3bblN9TVRkQSPO43wqliW5RubmyngPUi8ibOFfHPAucBR4EfcK6Qz9EOuERVs0RkLM5rfYc4E6gsFZHvgLuBk6raRUS6AMsp2ET39QSognM1K8C9AKra2f1Qnysi7YpoTziwWlWfcT9U3wM6qKpKgE3wkh/r8ZcxNym3wuntf1MKmxwmTm8+HjgX6Ojx2BSP2xOA293kNRyn4FNJ5QxHrQKWqGqKqh4A0tx/mkHuTzzOP3IHnA8CT1uAKBF50026x8jfABFZIiKrgItx2phjSgHr5Ge6+zsO5/UHiMX5IEZVZwOHC1i3sNe4IA8BQ1W1GfAB8A93eW3gM3HGp1/j9Pb84PFaHuW3D7hVHjEXJr825rVZVbvl/ADeDhfNU9WjqpqGM5zWEjgf+NEthJjOme/HZx5DKINwvnkmAD/ifLi1AC7EGS7L+d9YScFu9oh7qMfyWOATdxuJOKUNikr8WcDn7u1jQBowQUSupRjDqBWV9fj940ucXlZ/nNo6JSIirXF6rz1V9bCIfIjzD5XjhMftz3F7uECcqh6k5HLqHWV73M65Xwmnl/+Cqv6noA248XYFBuP02Ibh1FjPJc5Umf/CmUVtp4j8lYLb523MWfz2d3/G14K8vHiN81unPtBVVZe4i6bglMoG+BtOgr9GnOMpP+YTI5z+2ua8rkXJr43FkclvncG8bfSMzXP7hdV88Xx/BLgu78Qx7jezs60bU9D76NkeOL1NaTkfSu63r144hRxH4JQ+vvgsYyrXrMfvH+8Dz6nqqhKsmwLUdG/XwvnnOioiDYHLClrJ7anNwTno+EEJ9lscc4A7RKQGONNgils5NIc4Y+khqvo58Gcgxn3Is305/6jJ7rauL2Sfnut5awHOBw7u8FR+8/QW9hoXtM/DQG2P4YZL+a0ybG1+O6B/WzHj9bVtOMM2UPhrnWMJ0F+cyX8qAzcU8tw5wGiR3GM03d3lP+NMMJIzfWiXEsTtuY12ON8k1uO0p5uIhIhIc5yhxzO4f1u1VfUbnAP23fJ7XiCxHr8fqGoS8HoJV38HmCUie1R1gIjE44ytbwEWFrHuRJxyzT6diUxV54rIOcAv7v/5ceAWTq+R3xT4QERyOh8586d+iHOcIBXoA7yLM9SxDfi1kN1+CrwrzkFXb5IWOOPTk0VkOM5UlXtwkrlnW1YU8hqf9l54rJMpIncBn4tINs4HQc63mZeAj9zjA997GWdZeQWYKiK/w4vYVHWP+y3sF5zXbjnO3Nb5+RvOMYqVbvLfhnO85t84fwcrgQRgaQni/hfO38wqnF7+bap6SkQWAltx/n5WU/Dxg5rAF+43TMEZpgtoVp0ziIhzBkhtVf2zv2MpD0SkKpDlJuo+wL/znMZqTECyHn+QEJEZQBsCfOyymFrg9HBDgHTgLj/HY0yZsB6/McYEGTu4a4wxQcYSvzHGBBlL/MYYE2Qs8RtjTJCxxG+MMUHm/wFOE05eyPRdYQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Temp_ambient = np.array([55,58,60,65,66,67])\n",
    "t = np.array([8,9,10,11,12,13])\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(t,Temp_ambient,'o-',label='Ambient Temperature')\n",
    "plt.title('Ambient Temperature vs Time')\n",
    "plt.xlabel('Military Time starting at 8am in Hundred Hours')\n",
    "plt.ylabel('Temperature in Degrees F')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a better way to plot the ambient temperature with respect to time, which is to add more data points, preferrably every minute. This current data set assumes the single temperature reading at each hour persists for the enitre hour and then changes immediately to the next hour's ambient temeprature when the hour starts. It would be better to have temperature data for every minute so we could get a better idea of how the temperature changes throughout the day. We could assume that the temperature changes linearly from each data point, but real temperatures don't behave this predictably and we have no way of knowing how the temparatures in this model change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [],
   "source": [
    "#b\n",
    "import numpy as np\n",
    "\n",
    "T_0=85\n",
    "T_ambient = 65\n",
    "K = 0.61111111111112\n",
    "t = np.linspace(0,3,150)\n",
    "Temp_analytical = T_ambient + (T_0-T_ambient)*np.exp(-K*t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_t = 0.02\n",
    "Temp_numerical=np.zeros(len(t));\n",
    "Temp_numerical[0] = T_0\n",
    "for i in range(1,len(t)):\n",
    "    if 0<=i<= 50:\n",
    "        T_ambient = 65\n",
    "    elif 50<=i<=100:\n",
    "        T_ambient = 66\n",
    "    elif 100<=i<=150:\n",
    "        T_ambient = 67\n",
    "    Temp_numerical[i]=T_ambient+(Temp_numerical[i-1]-T_ambient)*np.exp(-K*delta_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f9758315ef0>"
      ]
     },
     "execution_count": 171,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3gUVffA8e9JsgsB6USkqGBD6UhEBAtWsGNBRfS1oz8FxQ42sKNgQcHCq4j6KoJIsYKiIjZUqjQREJAERERCDWSzOb8/ZjYMYZNskt0km5zP8+yTnZk7M3cyyZ69Ze4VVcUYY4zJK6GsM2CMMaZ8sgBhjDEmLAsQxhhjwrIAYYwxJiwLEMYYY8KyAGGMMSYsCxCmVIlIUxFREUkq67yY2BGR60Tko7LOhykZCxCmUCKyWkQyRWS7iGwWkU9E5MBSzsNi9/zbRSQoIrs8y/eVZl6iTURmicgVZZ2PohCRlQXcjztU9XVVPbes82lKxgKEidS5qrof0BDYALxYmidX1Zaqup+bh2+BvqFlVX2iNPNSFKVRUiqL0piqHuq5H78A13vux7OlnR8TGxYgTJGo6i5gAtAitE5EaonIWyKyUUTWiMgDIpLgbksUkWEi8o+I/AGc7dmvp4jM8R5fRO4UkcnFyZuI3Cgiy0TkX7eU09hdX9Wt1rrJ/ea71c1jcxH5WUS2iMg7oQ9aEekuIitE5GH3WH+ISE/PeZJF5HkRWSsif4nIiyJSJc++D4rIBuBlEUkRkc/c38+/IjJFRBq66Z8BjgFec799PyMiR4pIdp5ryy1luNfxlYiMFJHNwICCrj/M72mGiFyfZ90yETnLvV8j3LxuEZEFItK8GPfiJhGZXpzfv7vPBSLyq4hkiMi3ItIi/7OZWLEAYYpERKoBlwKzPKtfBGoBhwAnAf8BrnG33QCcA7QHUoGLPft9CDQTkaM8664A3i5Gvi4D+gPnAg2AecD/8iQ7FWjr5nGQm++eQDOgI3CRJ21TwA8cAPQB3hSRZu6254AmQGugOXAE7oe0Z18fcCBwK87/2SvAQe65QsdAVe9k72/gd0Z4yScC84H6wDMRXn/Iu0Cv0IKIdADqAp/j3KsOwKFAHeByYHOEeSpMRL9/EekEvITzN1QP5+9hclmUlCo9VbWXvQp8AauB7UAGkA2sA1q72xKB3UALT/obgRnu+6+AmzzbzgAUSHKXXwYed9+3xPkwqlJIfmbgfKB6130N9PYs+4AAzodlVfecHTzbFwO3eZZHAkPc992BXUBVz/YPgbuBJCALaOzZdjKw1LPvDsBXQP47Aes9y7OAKzzLRwLZefbJTQPcBPwe6fWHOX9dIBNo6C4/A7zkvj/L/d10BCTCv4+98u/J43T3fVF//28A9+c53hrg2LL+X6hsLytBmEj1UNXaQBWgL/CNiByA8w3Wj/MPHLIGCFVvNALW5tnm9SZwuYgIcCUwXlV3FyN/BwOvuFUSGcBGnGDWxJNmg+d9Zpjl/TzLG9WpTvPmu5H78gGLPeeaDOzvSfuXqgZCCyJSQ0RGi8ifIrIV55t6/WJco9faPMuRXD8Aqvov8AVwiVsVeCnwjrv5M+B14FVgg4i8JCL75T1GMUX6+z8YuC90Le71pLDnb8qUEgsQpkhUNaiqE4EgcDzwD8431YM9yQ4C0t3363GqWrzbvMebhfON/ASc6owiVy+51gJXq2ptzytZVecUumd49UWkqmf5IJyS03qcD95DPeeppar1PGnzDpE8AOeD+hhVrYlTipIC0u8AEkPtGq4D8qTJu09Rr38sTjXTSe71/ACgjmdVtT3QBqdK6LZ8jhEra4GH8lxLNffvzpQiCxCmSMRxPk799FJVDQLjgcfdb8oHA3ewp/57PHCriDQRkTrsXVcf8hYwAqda5btiZu0V4IFQg6qI1BGRiwrZpyA+4EER8YvIKcDpwAduyWA0MFxE6ru/jwNF5PQCjlUD2AlkiEh94IE82zfgtN+ErMMpAfR2G41vpvBvz0W9/ik4VXr3A2NVnXocEekkIqluff8OnOAdLOTc0TYK6OfmQ0RkPxE5z23/MqXIAoSJ1Ecish3YCjwOXKWqi91t/XA+TP4AvsNpBB3tbvsvMA1YAMwFwn0LfBtoRfFLD6jqWJwgM9GtxpmP86FeXKtxvln/hXMt16jqH+62/jgf4rOBLcBU4LACjjUMp0ppE87v59M8258D/iPOMyZPu0H3epyG3H9wSmAFloSKev2quhOnXeVUnPsVUhsYg9Pe9AdO1doLBZ072lT1e5zG/VfdfPyOU7q0yWtKmbhfHIwpMyKSDPwNHK2qy8tBfroDI1S1oA99Yyo8K0GY8uD/gF/KQ3Awxuxh/YpNmRKR1TgNtj3KOCvGmDysiskYY0xYVsVkjDEmrApVxVS/fn1t2rRpWWfDGGPixpw5c/5R1ZRw2ypUgGjatCmzZ88u62wYY0zcEJG8oxvksiomY4wxYVmAMMYYE5YFCGOMMWFVqDYIY0zsBAIB0tLS2LVrV+GJTblTtWpVmjRpgs/ni3gfCxDGmIikpaVRo0YNmjZtijM6u4kXqsqmTZtIS0ujWbNmhe/gimkVk4jcLs5k84tEZKw79eBgEUkXkfnu66x89u3uToO4QkTCjQAaFZPnpdNlyFc0G/AJXYZ8xeR56YXvZEwltGvXLurVq2fBIQ6JCPXq1Sty6S9mJQh3PtxbcWYayxSR8cBl7ubnVHVYAfsm4swwdTqQBvwiIh+q6pJo5nHyvHQGTlxIZsAZzTg9I5OBExcC0KO9zU1iTF4WHOJXce5drBupk4Bkd2z5ajhDJEeiI7BCVf9Q1SzgPeD8aGdu6LRl7A4EuDlxCm1kJQCZgSBDpy2L9qmMMSbuxCxAqGo6zjj4f+LMwrVFVT93N/cVkV/daRjrhNm9MXtPqZhGPhOmiEgfEZktIrM3btxYpDyuy8hkPzLpnTSd4b4RVMMpfqVnZFpVkzHlzLXXXsv+++9Pq1at9lo/ePBgGjduTLt27WjXrh2ffrpnuo0nn3ySww47jObNmzNt2rSwx+3atSvNmzfP3f/iiy8uMB+rV6/eJw/RMGbMGNati/Q7dOmIWYBwP/jPB5rhzONbXUSuwJmk/lCgHU7geCbc7mHWhR1VUFVHqWqqqqampIR9WjxfjWons5Xq3J51MwfL3wxOejN328CJCy1IGFMC0W7fu/rqq5k6dWrYbbfffjvz589n/vz5nHWW06y5ZMkS3nvvPRYvXszUqVO5+eabCQbDT473zjvv5O4/YcKEEuUzr/zOmVelChDAacAqVd3oTtM4EeisqhvceY1zcGYb6xhm3zT2nse4CZFXT0Xs7m7NSfYl8rMexYjg+VyS9A3nJPwIWFWTMSURat9Lz8hE2dO+V5IgceKJJ1K3bt2I00+ZMoXLLruMKlWq0KxZMw477DB+/vnniPe/+uqr9woW++233z5pgsEgd999N8cccwxt2rTh1VdfBWDGjBmcfPLJXH755bRu3Xqffa6++mpatWpF69atee6555gwYQKzZ8+md+/etGvXjszMTObMmcNJJ51Ehw4d6NatG+vXrwecEk///v3p3LkzrVq1yr2mb775JrcU1L59e7Zt2xbxteYnlt1c/wQ6ufPIZuJMbThbRBqq6no3zQXAojD7/gIcLiLNgHScxu3Lo53BUEN0/3HzeSH7Qo5PWMQTvteZt/sw0knJrWqyBmtj9vbwR4tZsm5rvtvn/ZlBVjBnr3WZgSD3TPiVsT//GXafFo1qMujclsXKz4gRI3jrrbdITU3lmWeeoU6dOqSnp9OpU6fcNE2aNCE9PXyA6t27N8nJyQCcfvrpDB06NKLzvv7669SqVYtffvmF3bt306VLF8444wwAfv75ZxYtWrRPt9L58+eTnp7OokXOR19GRga1a9dmxIgRDBs2jNTUVAKBAP369WPKlCmkpKQwbtw47r//fkaPdmby3bFjBz/88AMzZ87k2muvZdGiRQwbNoyRI0fSpUsXtm/fTtWqVYv2Swwjlm0QPwETcOYhXuieaxTwtIgsFJFfgZOB2wFEpJGIfOrumw30xZnLeCkw3jP/cVT1aN+YxrWTySaJ2wK3ICjP+0eS6M7TblVNxhRd3uBQ2PqS+L//+z9WrlzJ/PnzadiwIXfeeSfg9P3PK7+ePN4qpkiDA8Dnn3/OW2+9Rbt27Tj22GPZtGkTy5c7EyN27Ngx7DMHhxxyCH/88Qf9+vVj6tSp1KxZc580y5YtY9GiRZx++um0a9eOxx57jLS0tNztvXr1ApxS1datW8nIyKBLly7ccccdvPDCC2RkZJCUVPLv/zF9UE5VB+FMvO51ZT5p1wFneZY/Zd/J3WPi7m7NGThxIWsDDXggcA3D/S9xS+IUXghemFvVZKUIY/Yo7Jt+lyFfkZ6Ruc/6xrWTGXfjcVHNS4MGDXLf33DDDZxzzjmAU2JYu3ZPX5e0tDQaNWoU8XGTkpLIyXECmqqSlZW1TxpV5cUXX6Rbt257rZ8xYwbVq1cPe9w6deqwYMECpk2bxsiRIxk/fnxuycB73JYtW/Ljjz+GPUbeQCciDBgwgLPPPptPP/2UTp06MX36dI488siIrzccG4sJpxTx5IVOPeGUnOOZGDye25I+oIM4bRDWq8mYogm173kl+xK5u1vzqJ8rVDcPMGnSpNweRueddx7vvfceu3fvZtWqVSxfvpyOHcM1eYbXtGlT5syZAzjtGYFAYJ803bp14+WXX87d9vvvv7Njx44Cj/vPP/+Qk5PDRRddxKOPPsrcuXMBqFGjRm67QfPmzdm4cWNugAgEAixevKcSZdy4cQB899131KpVi1q1arFy5Upat27NvffeS2pqKr/99lvE15ofG2rD1aN9Y4ZOW0Z6RiYPBa4m1b+M4f6RnLl7CNuoZg/QGVMEof+TodOWsS4jk0a1k7m7W/MS/f/06tWLGTNm8M8//9CkSRMefvhhrrvuOu655x7mz5+PiNC0adPchuKWLVtyySWX0KJFC5KSkhg5ciSJiYlhj+1tg6hfvz7Tp0/nhhtu4Pzzz6djx46ceuqpYUsE119/PatXr+boo49GVUlJSWHy5MkFXkd6ejrXXHNNbunkySefBJxG8Ztuuonk5GR+/PFHJkyYwK233sqWLVvIzs6mf//+tGzplNzq1KlD586d2bp1a27p4/nnn+frr78mMTGRFi1acOaZZxbjt7y3CjUndWpqqpZkwiDvk9XtZTnv+x9mak5H+gb6AUKiCM9c0taChKmUli5dylFHHVXW2aj0unbtmtuYXVTh7qGIzFHVsAezKiYPb1XTPD2cZ7N7ck7iLC5P/AqAoKo1WhtjKg0LEHmEejUBvBw8l5nB1gxKeouj3Fn57PkIY0xZmjFjRrFKD8VhASKMUAObksDtgZvJoDojfC9QHadXhjVaG2MqAwsQYYSqmhJF2EQtbgv0pan8xWO+0YRG/LCqJmNMRWcBIh892jfmmUvakuxLZFZOC4ZnX8QFid/TM/EbwKqajDEVnwWIAngbrUcEe/B9sCWPJI3hcHGeaLSqJmNMRWYBohChRuscEugfuIXtJPOSbzjJ7tDgVtVkTHwZM2YMffv2LTSNd2TV66+/niVLij5f2YwZM3Kf7I5HFiAiEGq03khtbgvcwqGyjkd9YwC1qiZjKqC8AeK1116jRYsWZZijsmEBIgLeqqYfclrxYvACLk6cyWWJXwNW1WRMaenRowcdOnSgZcuWjBo1CnCG4b7//vtp27YtnTp1YsOGDQB89NFHHHvssbRv357TTjstd33Itm3baNasWe4wGVu3bqVp06a8//77+wy93bVrV0IP4U6dOpWjjz6atm3bcuqppwLOyK2dO3emffv2dO7cmWXLKsaXRhtqI0LeoTiGZ19Ie1nOw0ljWJTTlEV6iA3FYSqXzwbAXwuje8wDWsOZQwpMMnr0aOrWrUtmZibHHHMMF110ETt27KBTp048/vjj3HPPPfz3v//lgQce4Pjjj2fWrFmICK+99hpPP/00zzyzZ36yGjVq0LVrVz755BN69OjBe++9x0UXXUTPnj0ZOXJk2KeVN27cyA033MDMmTNp1qwZ//77LwBHHnkkM2fOJCkpienTp3PffffxwQcfRPf3UwasBFEEoaqmHBK4LXAL/1CLl33DqcV2MgNB7hy/wEoSxsTQCy+8kFtSWLt2LcuXL8fv9+fW83fo0IHVq1cDzuit3bp1o3Xr1gwdOnSvwe5Crr/+et544w0A3njjDa655poCzz9r1ixOPPHE3GG8QxMYbdmyhZ49e9KqVStuv/32sOeKR1aCKALvBEObqcnNWbcx3v8Iz/le4rrAXQQ1wUoSpnIo5Jt+LMyYMYPp06fz448/Uq1aNbp27cquXbvw+Xy5w18nJiaSnZ0NQL9+/bjjjjs477zzmDFjBoMHD97nmF26dGH16tV88803BIPBQueaVtWwc0o8+OCDnHzyyUyaNInVq1fTtWvXEl9veWAliCLyDsWxQA/jkez/cErifG5JnALY8xHGxMqWLVuoU6cO1apV47fffmPWrFmFpm/c2Pmi9uabb+ab7j//+Q+9evXaq/TgHXrb67jjjuObb75h1apVALlVTN5zjRkzpkjXVZ7lGyBE5HXP+ytKJzvxwTvW/TvBU/kgeDx3JE3g+ASn9GCN1sZEX/fu3cnOzqZNmzY8+OCDe00nGs7gwYPp2bMnJ5xwAvXr1883Xe/evdm8eXPuLG2wZ+jtUCN1SEpKCqNGjeLCCy+kbdu2XHrppQDcc889DBw4kC5duhAMBkt4peVHvsN9i8hcVT067/vyrKTDfRfF5Hnp3Dl+AUFVqrKbSf6HaCCbOWf3E6yjPsm+RJ68sLVVNZkKo6IO9z1hwgSmTJnC22+/XdZZiTkb7ruUeIfi2EUV/i/QnySCvOQfThWyrNHamDjQr18/BgwYwIMPPljWWSmXCmqkbiIizwLieZ9LVe+Iac7igLfRerU25K7ATYzyP8cjSWO4N/sGgoo1WhtTjr344otlnYVyraASxEBgMbDI8977KpSI3C4ii0VkkYiMFZGqIjJURH4TkV9FZJKI1M5n39UislBE5otI6dQbFYO30frznGMYnn0BlybN4IrE6YA1WpuKpSLNQFnZFOfe5VuCUNXX89sWCRFpDNwKtFDVTBEZD1wGfAEMVNVsEXkKJ/jcm89hTlbVf0qSj9Jwd7fmuVOVPp99ES1lNYOS3uL3nCb8rEflNlpbKcLEs6pVq7Jp0ybq1asXtqunKb9UlU2bNlG1atUi7Rfr5yCSgGQRCQDVgHWq+rln+yzg4hjnIeZCH/xOo3UCtwduYbL/QUb6h3Pe7sdZTz2rajJxr0mTJqSlpbFx48ayzoophqpVq9KkSZMi7ROzAKGq6SIyDPgTyAQ+zxMcAK4FxuV3COBzEVHgVVUdFau8RkPog3/gxIVsC1SjT+AOJvsf4lX/s/TMGkRmwM+d4xfsldaYeOLz+XKfIDaVQ8x6MYlIHeB8oBnQCKjufZ5CRO4HsoF38jlEF7dr7ZnALSJyYj7n6SMis0Vkdll/s/EO6rdSG3N74GbaJKziCd/rgBJUteHBjTFxo9AAISJPikhNEUkSkWkiskFELo/g2KcBq1R1o6oGgIlAZ/eYVwHnAL01n5YTVV3n/vwbmAR0zCfdKFVNVdXUlJSUCLIVW95G6+k5HXg2cDEXJX7L1YnTAGu0NsbEj0hKEGeq6lacD/S/gZbk36js9SfQSUSqidOidSqwVES6u/ufp6o7w+0oItVFpEboPXAGTm+quOB90vrFYA+mBVN5IOl/dLEnrY0xcSSSABFqpzgLGOv2Kiq0v5Sq/gRMAOYCC91zjQJGADWAL9wurK8AiEgjEfnU3b0B8J2ILAB+Bj5R1amRX1bZClU1JYqgJHBH4P9YoY15yTecZrIesJnojDHlX75DbeQmEBmK0w4QBFKBWjgf2MfGPntFU5pDbURi8rz03O6vTWQjk/0PskWrc0HWw2xlPxJFeOaSttZobYwpMyUaakNV7wZOATq4bQm7gAujm8WKydtonaYp3JTVnwPlb0b6XiCJbGu0NsaUa5E0UifjdEcNPZN+ANAmlpmqSLyN1rP1SO7Lvp4TEhfxYJIzMJg1WhtjyqtI2iBGu+lOcJfXAU/ELEcVkLfRekLwJF7NPpurkr7gisQvAGu0NsaUT5EEiMNV9QkgAOD2PLLn7IvA22gN8FR2L6YH2zM46U06Jzids6yqyRhT3kQSILJEpCpuzyURaQZkxTRXFZB3eHBnTuu+rNDGvOx7nkNknQ0PbowpdyIJEI8AU3GG/H4T+BpngD1TRN5G6x0kc33gLrJI4g3f09RlqzVaG2PKlUh6MU0FegI34D7RrKpfxjpjFZW30TpNU7gh6y4ayGZe8w/LnWio/7j5dBnylQUKY0yZinQsplOBVqo6GagiIh1imKcKz9toPV8P47bALbSTlTzrewkhB3Aarq00YYwpS5F0cx0BnAyEBtrbAbwSy0xVdHkbrafldOTx7Ms5O/Fn7k3aM7itdYE1xpSlSEoQnVX1RpwH5FDVfwF/THNVCXgbrQFeD57FW9mnc1PSR1yeuKcGz7rAGmPKSiTzQQREJIE9vZjqgVsPYkokNMTG0GnLSM/I5OHs/9BENvJI0hus03rMyGkH2LzWxpiyEUkJYiTwAZAiIg8D3wFPxTRXlUiP9o35fsApPH9pO/w+P30Dt7JMD2SE7wVaymoA6wJrjCkTkfRiegt4ABgGbAZ6qup7sc5YZRNql9hJVa7NupstVOcN/9M0kb8BrAusMabUFRggRCRRRBao6mJVHa6qz6tq3MzLEG9CXWA3UJersu7FT4C3fEOow1bAShLGmNJVYIBQ1SCwRESs8ruUhLrArtAmXJd1F41kE6P9w0h2+ghYScIYU2oiaYOojzMT3DQRmRh6xTpjlZW3C+wcbc6tgb60kZW5Q4SDdX81xpSOSHoxDYl5LsxeQr2VBk5cyOeBY3gw+1qe8L3O4zqae7NvACS3+6v1bDLGxEqhAcKG1SgboQ/+O8cv4N3gqTSQzdyWNJEN1ObZ7EsA6/5qjImtQgOEiGxm3zmotwCzgbtVdXUM8mXYuyTxXOAi9mcztyZN5h+txVvBbrmN1t60xhgTLZFUMb0IbADexZkH4jIgBVgBvIEzDIeJkdAHf/9x83kg+1rqyVYe8b3JVq3O5JzjcxutvWmNMSYaImmkPkNVR6rqZlX9V1VfAs5U1XeAujHOn2FP99cgifQL9OP7YEuG+V7h9ITZgHV/NcbERkSjuYrIhXneh2aUK3DIDRG5XUQWi8giERkrIlVFpK6IfCEiy92fdfLZt7uILBORFSIyINILqqhC3V9346dP4A4WaTNG+F7guITFgHV/NcZEXyQB4grgBhH5V0Q24cwLcaWIVAP657eT++zErUCqqrYCEnGqpwYAX6rq4cCX7nLefRNxhvg4E2gB9BKRFkW6sgrG2/11B8lclXUvq7Qhr/mG0U5WAFaSMMZEVyRDbaxQ1TNVta6q1nPf/66qO1X1m0J2TwKSRSQJqAasA84H3nS3vwn0CLNfR2CFqv6hqlnAe+5+lZp3BNgt7MeVWQP4R2sxxv8UR8hawEoSxpjoiWQ+iMPch+QWuMttRKTQKUdVNR1n/KY/gfXAFlX9HGigquvdNOuB/cPs3hhY61lOc9eFy18fEZktIrM3btxYWLbinrcksZE69A7cxy78/M//JAfJBsBKEsaY6Iikiuk14GH2tDcsZM/kQfly2xbOB5oBjYDqIlLofqHdw6zL29XWWak6SlVTVTU1JSUlwsPHN29JIk3354qsgSSRzTu+J2jAv4CVJIwxJRdJgKiuqj+EFlRVgUAE+50GrFLVjaoaACYCnYENItIQwP35d5h904ADPctNcKqnjMtbklihTbgqawC1ZTtv+5+0wf2MMVERSYDYJCLN2DNhUA/grwj2+xPoJCLVRERw5rVeCnwIXOWmuQqYEmbfX4DDRaSZiPhxGrc/jOCclYq3JLFQD+H6rLs4SP7mbf8QarIdsJKEMab4IgkQfYHXgSNFZA1Or6ObCttJVX8CJgBzcaqlEoBROGM7nS4iy4HT3WVEpJGIfOrum+2edxpOUBmvqouLdmmVg7ck8ZMexY2BOzhc0twgsQOwkoQxpnjEqTGKIKFILTd9RmyzVHypqak6e/bsss5GmZg8L52BExeSGQhySsJcXvE9x2JtxpVZA9hONQCSfYk8eWFre+LaGJNLROaoamq4bYVNGHSYiDwlIlOAt4ABInJYLDJpSsZbkvgq52j6Bm6llaxijP9pqpMJWEnCGFM0+QYIETkWZ/7pAE5weBsIAjNF5JjSyZ4pCm+bxOc5x9Av0I92soLR/qE24ZAxpsgKKkEMAnqr6gOq+oGqTlDV+4HeON1eTTnkLUlMzenIbYG+pMoyRvuGUZXdgJUkjDGRKShAHBZuLghV/Ro4NHZZMiXlLUl8ktOJ2wM30zFhKa/5hlGFLMBKEsaYwhUUILYVsG1HtDNiostbkvgwpwt3BW6ic8ISRvmezQ0SVpIwxhSkoPkgDhSRZ8OsF/IZ9sKUL94JhyYFTiApO8hQ3yhe5TluDNzObvw2n4QxJl8FBYiCxlu6L9oZMbHhnbr0/WBXBGVI0mu87hvKDYE7yaSqzUxnjAkr3wChqq+XZkZM7HhLEuMDJxPQJIb5XmGM/2muy7qL7VSzkoQxZh8RTRhk4p+3TWJSzgncGujH0bLcnrg2xuTLAkQlkrd3082B22gpq3jH/zi13T4J1rvJGBNiAaKS8ZYkvshJpU/gTg6XdMb6H6M+WwArSRhjHJFMGFRfRO4RkZdEZFToVRqZM7HhLUnMyGnHtYG7OVj+5j3/ozafhDEmVyQliClAA5xhN770vEwc85YkfshpxVVZ93KA/Ms4/6M0xpmZz0oSxlRukU4YdKeqvquq40KvmOfMxJy3JPGLHsmVWQOpK9sYV+VRmsl6wEoSxlRmkQSIz0TkjJjnxJQJb0linh5Or6wHqEoW7/sfpqWsAqwkYUxlFUmAuAmYKiLbReRfEdksIv/GOmOm9HhLEmW6WGkAACAASURBVIu1KT2zBrELP+/5H6NTwhLAShLGVEaRBIj6gA+oBaS4yymxzJQpfd6SxCptyMW7B7Fe6/Km7ynOSPgFsJKEMZVNQfNBHO6+bZnPy1Qw3pLEX9TjkqyHWKIH87LveXomzgCsJGFMZVLQWEwDgOuAkWG2KXBiTHJkypR37KYMrUHvrPt4xfccQ32jqMM2RgXPtbGbjKkkChqL6Tr35wmllx1THnjHbtoZqMp1gbt5jpe4zzeWurKdIdmXEVRs7CZjKriCShAlIiLNAW932EOAh4DjgObuutpAhqq2C7P/apw5KYJAdn6TapvY8JYkAprErYG+ZGh1bkr6iNps4/7s68gMYCUJYyowUdXYn0QkEUgHjlXVNZ71zwBbVPWRMPusBlJV9Z9Iz5OamqqzZ8+OQo5NyOR56QycuJDMQBBQ7kh6n1uTJvNVsB19A7eyk6oA1KnmY9C5LS1QGBNnRGROfl/AS2ssplOBlXmCgwCXAGNLKQ+mGLy9m0B4NvsS7g9cy0kJC3jP/ygpZACweWfAGq+NqWAiChAicpmI3O++P1BEOhTxPJexbyA4Adigqsvz2UeBz0Vkjoj0KSBvfURktojM3rhxYxGzZSLh7d0E8E7wNG4I3Mlhso6J/kEcKk5QyAwE6T9uPl2GfGWBwpgKIJLB+kYAJwNXuKt2AK9EegIR8QPnAe/n2dSLgksPXVT1aOBM4BYRCdtrSlVHqWqqqqampNjjGbGyd0kCvso5msuyHqCq7OYD/2COkd9y06ZnZFppwpgKIJISRGdVvRHYBaCq/wL+IpzjTGCuqm4IrRCRJOBC9m7E3ouqrnN//g1MAjoW4ZwmBvKWJH7VQ7kg6xE2aU3+53+CcxJ+zE1rD9UZE/8iCRABEUnAqfJBROoBOUU4R7iSwmnAb6qaFm4HEakuIjVC74EzgEVFOKeJkVBJonayD4A03Z+LsgazQA9lhP9Fbkj8GPdPxR6qMybORRIgRgIfACki8jDOsN9PRXJwEakGnA5MzLNpnzYJEWkkIp+6iw2A70RkAfAz8ImqTo3knCb2erRvzPxBZ/D8pe1oXDuZDGpwZdZAPg524n7fuzycNIYE9zuElSSMiV8RdXMVkZY43/oFmK6q5fLbvHVzLRuhrrC7AgEGJo2lT9InfBVsx62BvmynWm466wprTPlTUDfXAh+Uc59fmKuqbYHFsciciX/eh+qeyO7NGm3Aw0ljmOgfxHWBu1irDYA9XWG9+xhjyq8Cq5hUNQgsERH7bzYF8jZgvxM8jSsDA9lfMpjsf2ivHk5W5WRM/Ih0uO+lIjJNRCaGXrHOmIk/3q6wP+a0pEfWI2Tofrzjfzx3NFiwxmtj4kUkYzENiXkuTIXhHehvdaAhF2Q9zEjfCwz1jeIwSeep7F7kkGAjwhoTB0plLKbSYo3U5cfkeekM/nAxGZkBEgnyYNLbXJ30OdOD7bkt0JcdJOemtcZrY8pOicZiEpFtIrLVfe0Ukd0isjX62TQVibcrLJLE4OyreSBwDV0TFvCBfzBN5O/ctDaOkzHlU6EBQlVrqGpNVa0J7Af0BobHPGemQvA2Xv8veDpXBe6loWziI/8DHJ+wMDedNV4bU/4UaTRXVc1R1Qk4D78ZExFv4/X3Oa05N+txNmgd3vQN4cbEj7Anr40pnwptpBaR8zyLCUAqzgNzxkTM23j9Z6ABF2Y9zNO+UQz0jaV1wiruCfRhJ1Wt8dqYcqTQRmoReduzmA2sBl5V1b9imK9isUbq8s/beA1Kn8SPuTfpPZZrE24M3M4aPSA3rTVeGxN7BTVSRxIgOqnqrMLWlQcWIOLH5Hnp3Dl+AUFVuiQsZITvRRLI4bbALczIaZ+bLtmXyJMXtrYgYUyMlHRGuZfCrBtZsiyZys7beB1ql0jTFEb7htE3cRJig/0ZU+byLUGISEfgOOAuYKhnU03gElVtE/vsFY2VIOKPtyRRld084XudCxO/44vg0dwZuImt7Jeb1qqcjIm+4pYgquMMs5EEpHheWUDPaGfSVE7eksQuqnBH4P8YFLiKkxIW8GmV+2grK3LT2vMSxpSuSNogDlHVP0opPyViJYj4tXfjNbSVFYz0v8D+bObx7Ct4M3gGoc5ziSI8c0lbK0kYEwUlbaSuD9wJtASqhtar6hnRzGQ0WICIf94qp1ps5xnfy5yWOI9Pgh0ZEOjDNptfwpioKmkj9f9wurYegTOT3F/A/KjlzhgPb5XTFvbjhsCdPBHoRbeE2Xzkv5+Wsjo37eadAfqPm0/7Rz63aidjYiCSAJGiqq8CWar6JXAV0DG22TKVmXfeayWBUcFzuTTrQapIgIn+QVye+CWhp6/B2iaMiZVIAkTA/fmXiHQDWgEHxi5Lxuw92F+iCHO0OWfvfoKfco7kCd/rDPeNpAY7c9Nbd1hjoi+SNojzgG+Ag3Gef6gJPKyq5W7SIGuDqJhCc15nBoIIOdySOIX+SR+wTuvRP3ALc/WIvdJb24QxkSt2G4Q7J3VTVd2iqr+q6gmq2rY8BgdTceWtchoRvIBLsh5CgPH+R+ibOIkE98E6sLYJY6IlkjmpLyzOgUWkuYjM97y2ikh/ERksIume9Wfls393EVkmIitEZEBx8mAqDm+VU+1kH3P1CM7KepJPc47lLt/7jPU/RkM27bWPtU0YUzKRVDE9BtQA3gN2hNar6q8Rn8QpiaQDxwLXANtVdVgh6X/HGVY8DfgF6KWqSwo6j1UxVR57usPmcGHCtzziG0M2iQwI3MDUnL37UNhzE8bkr6TdXE8CjgaexmmDGAmMKGIeTgVWquqaCNN3BFao6h+qmoUTnM4v4jlNBbanO2wSE3NO5OysJ1ijDXjF/zxPJP2XZHblpg2qWpWTMcUQyYxyJ4R5nVjE81wGjPUs9xWRX0VktIjUCZO+MbDWs5zmrtuHiPQRkdkiMnvjxo1FzJaJZ962iTV6ABdnDebl7HO5LHEGH/vvp42s3Cu9tU0YUzSRzEmdIiKvisjH7nILEbk60hOIiB84D3jfXfUycCjQDlgPPBNutzDrwtaFqeooVU1V1dSUlJRIs2UqCG/bRPXkZJ7K7sUVgYEky24m+gdxe9L7JJG91z4WKIyJTCRVTGNwurmGnn1YjjP0RqTOBOaq6gYAVd2gqkFVzQH+S/iH7tLY+1mLJsC6IpzTVDLeQPGTtqb77qeYktOZ25ImMcn/EIdL2j77WCO2MQWLJEDsr6rvgtOPUFUDQLAI5+iFp3pJRBp6tl0ALAqzzy/A4SLSzC2BXAZ8WIRzmkoq1DYR8NXkzsDN3JjVn0ayiY/993N94id7dYcFe8DOmIJEEiB2iEhd3CoeETkG2BbJwUWkGk5PJO9zE0+LyEIR+RU4GbjdTdtIRD4FUNVsoC8wDVgKjFfVxZFdkqnsvG0T03I60m3308zMacMDvncY63+MJvL3XumtEduY8CLp5poKDMcZzXUBTmPxxapa7gbss26uJq89w4hncXHiTAYlvYWgPJp9JeOCXcnb3CU434Qa107m7m7NrWusqfBKNNy3ewA/cBTO/88St+tpuWMBwuQnFCiqZ65jmO9Vjktcwsxga+7Lvp40zb9zgw3bYSq6ks4HUQW4ETge58vVt8B/VXV3tDNaUhYgTGEmz0vnrvHz6JXwBQOSnKaxp7Iv4+3g6WgBNa4WKExFVVCASIpg/zeB3Tg9jsBpdH4Tp+HYmLgS+oAfODGJr3a350nfazzie5NzEmdxb6APq7Rh2P1CPZ68xzCmooskQLRQ1Tae5S9EZEGsMmRMrIU+4Ad/mMB/Mgdwcc5MHkx6m6n+ATybfTGvBc8iSOI++4V6PHmPYUxFFkkvpvluzyUARKQD8GPssmRM7O15bqI9P9bozmm7h/J1TjsG+sYy0T+I5vJn2P2sx5OpTCJpg1iE00C9yl3VDFiM8yyEqurRMc1hEVgbhCmJyXPT+P7D17lXX6MWO3g1eA4vZl/Abvz57mNtEybelbSR+tCCtqvqyoK2lyYLECYaPp21iODU+ziXb1iTsz8PZl/DzJy2Be5jgcLEq2h0c62JM9xFbptFUYb7Li0WIEw0fff5RBp9dx+HJKznw+BxPBq4ko3ULnAfCxQm3pS0BDEI6INTxRRKrMUY0TXmLECYaPtw9irWfPg4fWQSu/HzdPalvBM8tcAusWCBwsSPkgaIZUCb8vjcQ14WIEwsTJ6XzugpX3BvcBRdEhczL+cw7gtcx1I9uNB9LVCY8q6kAWIi0EdV/4lF5qLJAoSJpclz0/j5w1e4Q9+kNtsZE+zG8OyL2Ea1Qve1QGHKq5IGiA7AZOBXnAfmAFDVYs1VHUsWIExp+OSnxez8bBAX6XQ2UZMhgV5MzDm+0GonsEBhyp+SBohFwGhgIewZK1lVv4xmJqPBAoQpTV9//Tl1v7mPtixnbs5hPBS4mkV6SET7WqAw5UVJA8TM8tggHY4FCFPqcnKY89HLNJ33FHV0K+8FT2Zo9iVspmZEu1ugMGWtpAHiGWAnzoQ93iom6+ZqTMiuLax4/0Garnyb7ZrMM9k9eTd4atghO8KxQGHKSkkDxLdhVls3V2PC+XspG8ffRso/P7E050Aey76C73NaR7y7BQpT2kr8oFy8sABhygVVWDKFHZ/cR/Wd6XwZbM8T2ZezUiP/0LdAYUpLSUsQKcBjQGNVPUdEWgAdVXVM1HNaQhYgTLkS2AU/v0rg66eR7J38L/s0hmdfGHH7BFigMLFXUICIZDTXMcA3wIHu8nLgzuhkzZgKzFcVutyG7/YFJKVew1VJ05lZ9Q6uT/wEP4GIDrF5Z8BGjzVlJpIAsb+qvovbxVVVAzgjuRpjIlG9PpzzLHLzD9Q4rAsP+N7hy6r3cGbCT+wZvaZgFihMWYgkQOwQkbq4f8nu3BDbCttJRJqLyHzPa6uI9BeRoSLym4j8KiKTRCTs6GcislpEFrr7Wr2RiX/7HwVXTIArPuDAlLq87B/OR1UHcVzC4ogPYYHClKZI2iBSgeFAS2AB0Bi4WFXnR3wSkUQgHTgWaA58parZIvIUgKreG2af1UBqUYb4sDYIEzdygrBgLHz9BGxN53va8vjuS1miTYt0GGujMCVVrEZqEemkqrPc936cSYMEWKKqWUXMwBnAIFXtkmf9BTjBpneYfVZjAcJUdIFM+Pm/8O0zsCuDTzmeJ3dfxFptUKTDWKAwxVXcADE3WrPFichoYK6qjsiz/iNgnKr+L8w+q4DNOFVbr6rqqHyO3QdnOHIOOuigDmvWrIlGlo0pXZkZ8P1wmPUyOcEA4/Q0hu06n03UivgQgvPP0rh2Mnd3a27BwkSkTAOEW/pYB7RU1Q2e9fcDqcCFGiYTItJIVdeJyP7AF0A/VZ1Z0LmsBGHi3tb18M1TMPctshP8vBnszgu7zmQL+xX5UFaqMJEoboDIAPL9QFbV8yI8+fnALap6hmfdVcBNwKmqujOCYwwGtqvqsILSWYAwFcY/K2DGE7BoIoGkaozOPpORu7qxlepFPpQFClOQ4gaI5cD1+R1UVb+J8OTvAdNU9Q13uTvwLHCSqm7MZ5/qQIKqbnPffwE8oqpTCzqXBQhT4WxYAjOehKUfkuWryajAWby863R2kFzkQ1mgMOGUWRWTiFQD1gKHqOoWd90KoAqwyU02S1VvEpFGwGuqepaIHAJMcrcnAe+q6uOFnc8ChKmw1i+Ar5+E3z9jt782L2edw6u7TiGTqkU+VIJAjlpbhXEUN0BMLI+TAhXEAoSp8NLmOFVPK6azq0o9Xg+ezUvbT2InyRE+crcvK1lUbjZYnzEVzZ+zYMYQ+ONrqFqbpU2voM9vHVi7q0qxD2mBonKyAGFMRZU2B74dBss+BX8Nfj/4Um5c3olVuwqfJzs/FigqFwsQxlR0fy1yHrZbPAmSqrLyoIv5vz+68PuuyEeOzcvaKiqHkg73/QHOnNSfqWpOgYnLmAUIU+n9sxy+fRZ+HQcJiaxqcj5915zE4l11S3xoK1lUTCUNEKcB1wCdgPeBMar6W9RzGQUWIIxxbV7tPJk973+QE2Rt4+48+PcpzNjaMPeJ6+KyQFGxRKWKSURqAb2A+3G6rv4X+J87/He5YAHCmDy2roMfR8KcMZC1HZqdCJ1vY/K2Ixn80RIyMov/72tVUBVDiQOEiNQDrgCuxBk24x3geKC1qnaNXlZLxgKEMfnIzIC5b8KsV2DbOti/BRzXlyk5nRn08fISBYoQK1nEp5JWMU0EjgTexqleWu/ZNju/A5cFCxDGFCI7CxZ9AD+8CH8vhhoN4dgb+cTfnfs/WxuVQGEli/hS0gBxiqp+FZOcRZkFCGMipAorv4IfXoA/ZoB/Pzj6Kj6v0YOHv91OekZmidsqQixglG/FfZK6wKeoVXViFPIWVRYgjCmG9b/CjyOckoXmQPOzoGMfaHYik+evY/CHi6NSsgixqqjypbgB4g337f5AZyBUijgZmFEeh+GwAGFMCWxJh9mvOw3aOzdBylFwbB9ocymTF2dEPVBYyaJ8KGkV08fADaG2BxFpCIy0AGFMBRXY5ZQmfnoF/voVqtaC9lfCMdczeY2fodOWRbUKKqSaL4EqvkQydgZoZEGj1JQ0QCxS1Vae5QTgV++68sIChDFRpAprf4KfXoUlU5zqpyO6O6WKQ04GESbPS496ycLLShmxV9IAMQI4HBiL84XhMmCFqvaLdkZLygKEMTGydR3MHg2z34Cd/0C9w6HD1dDucqhWl8nz0mNWsvCygBF90XgO4gLgRHdxpqpOKih9WbEAYUyMBXY54z3NHg1pP0NiFWhxPqReAwcdByIAFjDiSDQCRAOgI859/llV/45uFqPDAoQxpeivRU6D9q/jYPdWqN/cKVW0vQyq7T32U6yrokIsYBRdSauYLgGGAjMAAU4A7lbVCVHOZ4lZgDCmDGTtcEsVb0D6bKdU0fICJ1gc1Cm3VAGlV7IIsYBRuJIGiAXA6aFSg4ikANNVtW3Uc1pCFiCMKWN/LXQCxa/jIWsbpBwJ7XpDm0uhRoN9kocCxrqMTGol+8jKDrIzELtBoy1g7KukAWKhqrb2LCcAC7zrygsLEMaUE7u3w+KJMPctSPsFJBEOP90JFkd0hyR/vruWZinDAkbJA8RQoA1OLyaAS3G6ud4b1VxGgQUIY8qhjctg/ruw4D3Y/hck14U2lzjBomGbQne3gBFb0WikvhBn9FbBejEZY4ojmO2M/zT/HWeK1GAWHNAa2l0BrXtC9XoRHcYCRnRFbcpREakPbNIIdhKR5sA4z6pDgIeAt9z1TYHVwCWqujnM/t2B4UAi8JqqDinsnBYgjIkTO/+FhROcYLF+PiT44LDToE1POOJM8Ec+p3ZZBIxEEYKqFSJwFHcspk7AEOBf4FGc4b7rAwnAf1R1ahEykAikA8cCtwD/quoQERkA1MlbXeWm/x04HUgDfgF6qeqSgs5jAcKYOPTXIlgw1hneY9t6Z2TZo851ShXNToLEpCIdrrR7SkF8DxNS3AAxG7gPqAWMAs5U1VkiciQwVlXbFyEDZwCDVLWLiCwDuqrqendcpxmq2jxP+uOAwarazV0eCKCqTxZ0HgsQxsSxnCCs/g4WjoclHzrPVlTfH1pd5JQsGh29V5fZSJVFwID4KW0UN0DMV9V27vulqnqUZ9u8IgaI0cBcVR0hIhmqWtuzbbOq1smT/mKgu6pe7y5fCRyrqn3DHLsP0AfgoIMO6rBmzZpIs2WMKa8Cu2D5NKe77PLPnfaKeoc5pYpWF0H9w4t96LIKGCHlrV2juAFirqoenfd9uOVCTu7Hmaa0papuiDBA9AS65QkQHQsb/8lKEMZUQJmbnRLFwvedEgYKDVpBix7OA3n1DyvR4ctLwCirkkZxA0QQ2IHTcykZ2BnaBFRVVV+EJz8fuEVVz3CXrYrJGFM8W9c5I8sungxrZznrGrSClj2gRcmDBZR9wAgprcARtV5MxTz5e8A0VX3DXR6K0xMq1EhdV1XvybNPEk4j9ak4jdu/AJer6uKCzmUBwphKZEs6LP3QGeZj7U/OugatoeX5UQsWsHfACH1Yl4fAEa2AUWYBQkSqAWuBQ1R1i7uuHjAeOAj4E+ipqv+KSCOc7qxnuenOAp7H6eY6WlUfL+x8FiCMqaS2pDsliyWT9wSL/VvCUefAkWfDAW2K1cBdkNIeJiQ/JQ0YZVqCKE0WIIwxbElz2ix++xj+/NGZ6KjWQU6gOOocOLBTkbvORqqsSxvJvkSevLB1kYKEBQhjTOW04x9Y9hn89onzFHdwtzPUR/Mz4chz4NCTwZcc82yUZrtG49rJfD/glIjTW4Awxpjd22Hll06w+H0q7NoCvmpw6Clw+BnOq2bDUslKLEsaAqwacnbk6QsIELEpZxljTHlTZT9n9rsW50Mw4HSZ/e1jWDbV+QnQsC0c3g2O6OY8mJeQEJOs9GjfOGw1UDQCR6Pa0SsRWQnCGFO5qcLfS+D3ac5DeWt/ctotqtV3hig//AynlJFcu/BjxUikVVTWBlEACxDGmBLb+S+s+NJ5knvFdOdBPUl05tw+4gynhJHSPOq9oooiXEnDejEVwgKEMSaqgtnONKq/T3Nef7uPYtVs7DRwH3oKNOsa8VDl5ZEFCGOMiYaMtU6pYuVXsOobp6EbcdouDj3FeR14bIEz5pU3FiCMMSbacoKwbp4TLFZ+5UytmpMNvurQ9Pg9JYz6R5RpdVRhrBeTMcZEW0IiNEl1XifdA7u2Oj2jQgFj+TQnXc3G0OxEaHoCNDsBah9UtvkuAgsQxhgTDVVrwpFnOS+AzWvgj6/dYPG5MykSQO2D9wSLpidArfI1P4SXVTEZY0ys5eTAxqVOCWPVTFjzvdM7CqDuIU6VVNMTnaBR44BSzZq1QRhjTHmSkwMbFsHqb52gsfp72L3F2VbvcGjaxelWe1Anp8QRwzYMCxDGGFOe5QThr19h1bdO0Pjzpz0Bo0YjJ1CEAkaDlk77R5RYI7UxxpRnCYnQqL3z6nKrU8L4e4kzGu2fs5yfiyc6aavUhAM77gkajTvEbMBBCxDGGFPeJCTAAa2cV8cbnHUZa92A4QaNrx5z0/qgyTFw9SdRHzvKAoQxxsSD2gc6rzaXOMs7/3WevfjzR+d9DAYWtABhjDHxqFpdZ9TZI7rF7BSxGcvWGGNM3LMAYYwxJiwLEMYYY8KKaRuEiNQGXgNa4cxxcS3QH2juJqkNZKhquzD7rga2AUEgO79+usYYY2Ij1o3Uw4GpqnqxiPiBaqp6aWijiDwDbClg/5NV9Z8Y59EYY0wYMQsQIlITOBG4GkBVs4Asz3YBLgFOiVUejDHGFF8s2yAOATYCb4jIPBF5TUSqe7afAGxQ1eX57K/A5yIyR0T65HcSEekjIrNFZPbGjRujl3tjjKnkYhkgkoCjgZdVtT2wAxjg2d4LGFvA/l1U9WjgTOAWETkxXCJVHaWqqaqampKSEqWsG2OMidlgfSJyADBLVZu6yycAA1T1bBFJAtKBDqqaFsGxBgPbVXVYIek2AmuKmeX6QEVp76go11JRrgPsWsqjinIdULJrOVhVw367jlkbhKr+JSJrRaS5qi4DTgWWuJtPA37LLzi4VVEJqrrNfX8G8EgE5yx2EUJEZleUnlIV5VoqynWAXUt5VFGuA2J3LbHuxdQPeMftwfQHcI27/jLyVC+JSCPgNVU9C2gATHLasUkC3lXVqTHOqzHGGI+YBghVnQ/sE9VU9eow69YBZ7nv/wDaxjJvxhhjCmZPUu8xqqwzEEUV5VoqynWAXUt5VFGuA2J0LRVqRjljjDHRYyUIY4wxYVmAMMYYE1alChAi0l1ElonIChEZEGa7iMgL7vZfReTosshnJCK4lq4iskVE5ruvh8oin4URkdEi8reILMpnezzdk8KuJV7uyYEi8rWILBWRxSJyW5g0cXFfIryWeLkvVUXkZxFZ4F7Lw2HSRPe+qGqleAGJwEqcIUD8wAKgRZ40ZwGfAQJ0An4q63yX4Fq6Ah+XdV4juJYTcZ64X5TP9ri4JxFeS7zck4bA0e77GsDvcfy/Esm1xMt9EWA/970P+AnoFMv7UplKEB2BFar6hzoDB74HnJ8nzfnAW+qYBdQWkYalndEIRHItcUFVZwL/FpAkXu5JJNcSF1R1varOdd9vA5YCjfMki4v7EuG1xAX3d73dXfS5r7y9jKJ6XypTgGgMrPUsp7HvH0okacqDSPN5nFsc/UxEWpZO1qIuXu5JpOLqnohIU6A9zrdVr7i7LwVcC8TJfRGRRBGZD/wNfKGqMb0vsX6SujyRMOvyRt9I0pQHkeRzLs4YK9tF5CxgMnB4zHMWffFyTyIRV/dERPYDPgD6q+rWvJvD7FJu70sh1xI390VVg0A7cSZjmyQirVTV2+YV1ftSmUoQacCBnuUmwLpipCkPCs2nqm4NFUdV9VPAJyL1Sy+LURMv96RQ8XRPRMSH84H6jqpODJMkbu5LYdcST/clRFUzgBlA9zybonpfKlOA+AU4XESauWNDXQZ8mCfNh8B/3J4AnYAtqrq+tDMagUKvRUQOEHcwKxHpiHOvN5V6TksuXu5JoeLlnrh5fB1YqqrP5pMsLu5LJNcSR/clxS05ICLJuIOe5kkW1ftSaaqYVDVbRPoC03B6AY1W1cUicpO7/RXgU5xeACuAnewZXLBcifBaLgb+T0SygUzgMnW7OZQnIjIWpxdJfRFJAwbhNL7F1T2BiK4lLu4J0AW4Eljo1ncD3AccBHF3XyK5lni5Lw2BN0UkESeIjVfVj2P5GWZDbRhjjAmrMlUxGWOMKQILEMYYY8KyAGGMMSYsCxDGGGPCsgBhjDEmLAsQlYyIbM+zfLWIjCjF83cSkZ/cUTOXishgd/15EmZU2iiet6uIdPYsjxGRi0twvNUistC9joUiUqSxsERksIjcVdzzR5OI3Orei3dieI6mks8ot/mklG67qwAABmtJREFUry0iNxewPezIuSLSU5yRTnNEZJ/pjk3RWIAwUeH2zY7Em0AfVW0HtALGA6jqh6o6JFb5w3k+oXNhiYroZPc6LgZeiPKxS9PNwFmq2rusM+JRGydf+RnDvk8RAywCLgRmxiBPlY4FCJNLRA4WkS/FGUf+SxE5yF2/17ftUCnE/Vb+tYi8i/MgUnUR+cQd9GyRiFwa5jT7A+vBGVdGVZe4x8otybjne0FEfhCRP/Kc+x73G/sCERnirjtURKaKyBwR+Vb+v71zDbGqCsPw80aZjpNTagoVOGJiBaagEJbXAilSSCqiNEwNKrqgUhBIJRFe0vpjYaChYAZmNpZFOXgp75j3S2YlWtgPQ8wotUj9+vF9e2Z33NsZihil9cDmrLP2uu21zlnfupzzLumGiueqBR4HJsSIf0DcGliSx3OSvox6OEdzv4B2wM+5+BPj+fdIGp/znyQ/w2MF0CNX9m25MN0lba0of6fMT1IvSZZrmwOSqiQNj5nZdkkrJHWO+4PUeM7BdklXVKT9Fi4b/5GkCZLaS1oaz75J0s0R7m8znni22rj2SZoTI/d6+b98kdQn2mkj8GRRxUmqjs/atoqZ2DSgW5R7RmW8MuVcM9tnZvuL8kr8A/6NVni6Lr4LOAPsyF0/AG/EvWXA6HCPBZaGez5wXy6N3+J1MHAC6Brv7wXm5MLVFOT/It6Z1gGPAa3D/5FcOeYDi/EBzE24tDnAXcAGoCret4/XlUD3cN8CrCrIdzLwbO59WR5D8QPgFfc+BgYWpHcI2I2PWE8Cw8K/T/i3BaqBvbiCaOZfhRuU77LyAKuB3uGeAjxdkN/eiPcULrUyEugCbIz7V9H4x9dHgddybXpbuKuBS0uepWO4ZwEvhft2YEdJ/e0BauM6nSv/e8CocO8CBoV7BgXnZOBqDu3C3THqRZFu4bkaubilYXCdor4t/X272K//jdRGooFT5ssigI/cgWytth8+PQdYALzajPQ2m9nBcO8GZkqajh/AsrYysJm9HGvdQ4GHgAdxQ1PJUjM7C3yVjYZx7Zl5ZnYy0jomV+m8FVgsNQhZXt6McpflMTSu7fG+Glf2LFqyGGJmRyV1A1ZK+hzoD9SZ2QkASR8AA3BjU5eVXVJeO2suMEbSROAB/LyPSjbgshEDcSNyJ96RZnV8HbBIrv3fCsjaZD3wetT5B2Z2uIk66Y8besxslaQOkmqaiHPQzDIZi61AbcS50sy+CP8FuIGvRMAUSQOBs7g0deeCcIkWIC0xJc5HpsNymvisyHvhVrkwJxoCm31D40h5qkqObjSzA2Y2G7gD6CWpQ0GwP3Ju5V4rtWEuAY6bWe/cdWOznq48j6m5tK43s7fPl4iZHQCO4DORIrnlhqAl/kvwznMYsNXMioTi1uKGpgvwIdAL78wzwzULn4H1JGZmUbZp+IyiDbCpcvmtgDK56IbPQNA6587X4xl8VlDUVkWMBK4G+sTA5UhF2okWJBmIRJ4NuDIs+Bd3XbgP4R0/+IlVlxVFlnQNcNLM3gFm4sdvVoa5W41D/e54h3K8meWrB8ZKqoq02ptr+x+UdH/4SVKvgri/4kdONsXyyKM60rtWUqfzRYj7XYHv8Q77ntgXaAuMwDv3NcAISW1iH2B4Ft/Mfo98ZwPzSrJZA4wCvo1ZzzFclG193K8Bfgz36FzZupnZbjObDmwBmjIQa/C2R9Jg4GjU8SGiPeXnHHc9XyLmctS/SOofXmUb4DXAT2b2p6QhuAGE5rdX4j8kGYhEnmfwpY5duAJmdsD7HGCQpM34Gv+Jkvg9gc1y1cxJwCsFYR4G9keYBcBI80NQmsTMPsPljLdE/GzTdCQwTtJOfK2+6Ceny/AOOr9JXZRHPfAusFHSbuB9yjuq1VGO1cDzZnbE/HjL+cBm/OSyuWa2PfwX4fs+S2hcGspYiI+460vKdSic2YxhHT5zyjbHJ+PLbGuBo7mo42NDeSeuVPpp2bPn0ukbn4FpNBqbJUD7eN4n8LOdm2IM8GZsUp8qCbMw8tuCt+PXADGLWh9lP2eTWq6cuxHoIemwpHHhP0KupNsP+ETS8maUM1FCUnNNJC4A4hdCNWb2QkuXJZHISJvUiUQLI6kO6Ib/aiiRuGBIM4hEIpFIFJL2IBKJRCJRSDIQiUQikSgkGYhEIpFIFJIMRCKRSCQKSQYikUgkEoX8BUwdDRg9wDMUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(t,Temp_numerical,'o-',label='150'+' Euler steps')\n",
    "plt.plot(t,Temp_analytical,label='analytical')\n",
    "plt.title('Body Temperature vs Time')\n",
    "plt.xlabel('Hours Since the Body was found at 11')\n",
    "plt.ylabel('Body Temperature in Degrees F')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thanks to our solution from problem 3 where the time of death was 10:09 with an ambient temperature of 65 degrees, we can assume that the new caluclated time of death will be later in the day if the ambient temperature from 10-11am is 60 degrees, because the lower ambient temperature will cause the corpse to cool down faster to hit 85 degrees at 11 am. With this assumption in mind and the equation from problem three we can calculate the following time of death."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "60\n"
     ]
    }
   ],
   "source": [
    "ambient_temp(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the time of death is 0.710797829884799 hours before the body was initially found at 85 degrees F at 11am\n",
      "Therefore the time of detah is about 10:17am\n"
     ]
    }
   ],
   "source": [
    "T_found = 85\n",
    "T_0 = 98.6\n",
    "K = 0.61111111112\n",
    "T_ambient = 60\n",
    "t = np.log((T_found-T_ambient)/(T_0 - T_ambient))/(-K)\n",
    "print('the time of death is',t,'hours before the body was initially found at 85 degrees F at 11am')\n",
    "print('Therefore the time of detah is about 10:17am')"
   ]
  }
 ],
 "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
}