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": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6111111111111112"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "def k_calculation(temp_current,temp_initial,temp_ambient,time):\n",
    "    delta_T=temp_current-temp_initial\n",
    "    delta_time=time\n",
    "\n",
    "    k=-(delta_T/delta_time)/(temp_current-temp_ambient)\n",
    "\n",
    "    return k\n",
    "\n",
    "k_calculation(74,85,65,2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Change your work from problem 1 to create a function that accepts the temperature at two times, ambient temperature, and the time elapsed to return $K$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6111111111111112"
      ]
     },
     "execution_count": 192,
     "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",
    "    #change in temperature is dT\n",
    "    dT=Temp_t2-Temp_t1\n",
    "    #find K using formula in the problem statement\n",
    "    K=-(dT/delta_t)/(Temp_t2-Temp_ambient)\n",
    "    \n",
    "    return K\n",
    "\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",
    "\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": 221,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " The analytical solution is  70.89149787542375 \n",
      "\n",
      " The Euler solution for a time step of 60 minutes is  72.78 \n",
      " The Euler solution for a time step of 12 minutes is 71.1885682409906 \n",
      " The Euler solution for a time step of 1 minute is 70.92096677548729 \n",
      " The Euler solution for a time step of 0.1 minutes is 70.89560860759411\n",
      "\n",
      "\n",
      " As you can see, the solutions converge with decreasing time steps.\n"
     ]
    }
   ],
   "source": [
    "#a. \n",
    "\n",
    "import numpy as np\n",
    "\n",
    "#solve the problem analytically with the given data and formulas, and calculated K value \n",
    "#given the original temp, ambient temp, and time elapsed, it returns the body temperature\n",
    "def analytical_solution(temp_original,temp_ambient,hours):\n",
    "    exponent=-0.611111*hours\n",
    "    analytical_solution=temp_ambient+(temp_original-temp_ambient)*np.exp(exponent)\n",
    "    \n",
    "    return analytical_solution    \n",
    "    \n",
    "    \n",
    "    \n",
    "#the euler strategy recieves the time elapsed and the number of time steps requested, and it returns\n",
    "#the body temperature. it creates an empty list then fills the values with the euler approximations of \n",
    "#the body temperature using the slope and time passed since last approximation \n",
    "# the function only returns the last value in the list, the result\n",
    "\n",
    "def euler_solution(hours,n):\n",
    " \n",
    "    t=np.linspace(0,hours,n)\n",
    "    temp=np.zeros(len(t))\n",
    "    temp[0]=85\n",
    "\n",
    "    for i in range(0,len(t)-1):\n",
    "        temp[i+1]=temp[i]-(0.611*(temp[i]-65))*hours/n\n",
    "            \n",
    "    return temp[-1]\n",
    " \n",
    "print(' The analytical solution is ', analytical_solution(85,65,2),\n",
    "      '\\n\\n','The Euler solution for a time step of 60 minutes is ', euler_solution(2,2),\n",
    "      '\\n','The Euler solution for a time step of 12 minutes is', euler_solution(2,10),\n",
    "      '\\n','The Euler solution for a time step of 1 minute is', euler_solution(2,100),\n",
    "      '\\n','The Euler solution for a time step of 0.1 minutes is', euler_solution(2,1000)) \n",
    "\n",
    "print('\\n\\n','As you can see, the solutions converge with decreasing time steps.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The final temperature of the body can be seen using the analytical model with increasing elapsed time inputs: \n",
      " Body temp after 1 hr =  75.85495082939002 \n",
      " Body temp after 5 hr =  65.94193149832891 \n",
      " Body temp after 100 hr =  65.0 \n",
      " Body temp after 1000 hr =  65.0\n"
     ]
    }
   ],
   "source": [
    "#b\n",
    "\n",
    "print('The final temperature of the body can be seen using the analytical model with increasing elapsed time inputs:',\n",
    "      '\\n','Body temp after 1 hr = ',analytical_solution(85,65,1),\n",
    "      '\\n','Body temp after 5 hr = ',analytical_solution(85,65,5),\n",
    "      '\\n','Body temp after 100 hr = ',analytical_solution(85,65,100),\n",
    "      '\\n','Body temp after 1000 hr = ',analytical_solution(85,65,1000))\n",
    "      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At 10:09AM, 0.85 hours before the body was found, the temperature would be  98.62186581295349 according to the analytical model. Thus we can guess the time of death was around 10:09am\n"
     ]
    }
   ],
   "source": [
    "\n",
    "#c.\n",
    "    #using the analytical solution and plugging in values for the time elapsed, we see that the body temperature\n",
    "    #would have been around 98.6 about 0.85 hours before the body was found at 11am\n",
    "    \n",
    "import numpy as np\n",
    "def analytical_solution(temp_original,temp_ambient,hours):\n",
    "    exponent=-0.611111*hours\n",
    "    analytical_solution=temp_ambient+(temp_original-temp_ambient)*np.exp(exponent)\n",
    "    \n",
    "    return analytical_solution    \n",
    "    \n",
    "\n",
    "print('At 10:09AM, 0.85 hours before the body was found, the temperature would be ',analytical_solution(85,65,-0.85),\n",
    "      'according to the analytical model. Thus we can guess the time of death was around 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": 273,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The temperature at 8:30am is 56.5 \n",
      " The temperature at 9:30am is 59.0 \n",
      " The temperature at 10:30am is 62.5 \n",
      " The temperature at 11:30am is 65.5 \n",
      " The temperature at 12:45pm is 66.75 \n",
      " Yep! That all looks right to me.\n"
     ]
    }
   ],
   "source": [
    "#a\n",
    "#the function returns an approximation of the weather given a time of day relative to 11am\n",
    "\n",
    "import numpy as np\n",
    "from scipy import interpolate\n",
    "\n",
    "#interpolate from the given data values the temperature at any give time using a linear interpolation\n",
    "\n",
    "def amb_temp(hours):\n",
    "    time_values=[-3,-2,-1,0,1,2]\n",
    "    temp_values=[55,58,60,65,66,67]\n",
    "    \n",
    "    return np.interp(hours,time_values,temp_values)\n",
    "\n",
    "print('The temperature at 8:30am is', amb_temp(-2.5),\n",
    "     '\\n','The temperature at 9:30am is', amb_temp(-1.5),\n",
    "      '\\n','The temperature at 10:30am is', amb_temp(-0.5),\n",
    "      '\\n','The temperature at 11:30am is', amb_temp(0.5),\n",
    "      '\\n','The temperature at 12:45pm is', amb_temp(1.75),\n",
    "      '\\n', 'Yep! That all looks right to me.')\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 310,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The improved Euler model predicts a time of death of 10:15am. This is about 6 minutes different from the original linear analytical model.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXwV5b348c/3nJOdBJIQIBAgLBFkMSgBEdCiKC5VXNBK614r7a1b7b21eOtt6fqzvd5WrdvVul4t1uK+F3GpIKJhEcNW1kDYEpKQhew5398fcxICBjgkOWeSnO/79ZrXzDwzc+Z7IuabeZ55nkdUFWOMMQbA43YAxhhjOg9LCsYYY5pZUjDGGNPMkoIxxphmlhSMMcY087kdQHv07t1bMzMz3Q7DGGO6lOXLl+9T1bTWjnXppJCZmUlubq7bYRhjTJciIvlHOmbVR8YYY5pZUjDGGNPMkoIxxphmXbpNwRjTvdTX11NQUEBNTY3boXQLsbGxZGRkEBUVFfQ1lhSMMZ1GQUEBiYmJZGZmIiJuh9OlqSrFxcUUFBQwZMiQoK+z6iNjTKdRU1NDamqqJYQOICKkpqYe91OXJQVjTKdiCaHjtOVnGZFJYef+au59bwM7SqrcDsUYYzqVkCUFEXlSRApFJK9FWYqILBSRjYF1cqA8U0SqRWRVYHk0VHEBVNY08OCHm/hiW0kob2OM6YK8Xi/jxo1rXu65556jnv/0009zyy23hCm60AtlQ/PTwIPAsy3K5gKLVPUeEZkb2P9p4NhmVR0XwniaDe/Tg/hoL6sLyrjslIxw3NIY00XExcWxatWqkH1+Q0MDPl/nfccnZE8KqvpP4PA/xS8GnglsPwNcEqr7H43XI4wZ0JNVO/a7cXtjTBeUmZnJvn37AMjNzWXatGlfO6eoqIhZs2YxYcIEJkyYwJIlSwCYN28ec+bMYcaMGVx77bXhDPu4hTtd9VXV3QCqultE+rQ4NkREVgLlwN2q+kkoA8nO6MkzS/Opa/AT7YvIphVjOrVfvrGGtbvKO/QzR/VP4hcXjT7qOdXV1Ywbd7DS4q677uLKK68M6vNvv/127rjjDqZOncr27ds599xzWbduHQDLly9n8eLFxMXFtf0LhEFneYbZDQxS1WIRGQ+8KiKjVfVr/yJEZA4wB2DQoEFtvmH2wF7UfbKVDXsqGJvRs82fY4zpXtpTffT++++zdu3a5v3y8nIqKioAmDlzZqdPCBD+pLBXRNIDTwnpQCGAqtYCtYHt5SKyGTgB+NoQqKr6GPAYQE5OjrY1kOyMXgB8WbDfkoIxndCx/qIPN5/Ph9/vBzjiu/9+v5+lS5e2+ss/ISEhpPF1lHDXm7wOXBfYvg54DUBE0kTEG9geCmQBW0IZSEZyHCkJ0awusHYFY8yxZWZmsnz5cgBeeumlVs+ZMWMGDz74YPN+KBusQyWUr6TOB5YCI0SkQERuBO4BzhGRjcA5gX2AM4DVIvIlsAD4gaqG9H1REeGkjJ58uaMslLcxxnQxTW0KTcvcuXMB+MUvfsHtt9/O6aefjtfrbfXaBx54gNzcXE466SRGjRrFo4+G9O36kBDVNtfAuC4nJ0fbM8nOHxf+iwc/2MhX884lIaazNK8YE7nWrVvHiSee6HYY3UprP1MRWa6qOa2dH9Gv3Ywb2BO/Qt5Oe1owxhiI8KRwUqCxeXWBJQVjjIEITwq9e8QwoFccX1pjszHGABGeFACyB/a0pGCMMQGWFDJ6saOkmpIDdW6HYowxrov4pHBSi05sxhgT6SwpZPTEI7ByuyUFY4zjlVdeQURYv359mz/j+uuvZ8GCBUc953e/+90h+5MnT27TvebNm8e9997bpmsPF/FJISHGx4npSSzPt7kVjDGO+fPnM3XqVF544YWQ3ufwpPDpp5+G9H7BiPikADB+cDKrtu+nodHvdijGGJdVVlayZMkSnnjiieak8NFHHzFt2jQuv/xyRo4cyVVXXUVTx99f/epXTJgwgTFjxjBnzhwO7xC8aNEiLr300ub9hQsXctlllzF37tzm3tNXXXUVAD169Gg+7w9/+ANjx44lOzu7uVf1448/zoQJE8jOzmbWrFlUVXX87JGR2Y131yp49Ydw6SOQns34wck8uzSfDXsrGN3fBsczplN4Zy7s+apjP7PfWDj/6DOpvfrqq5x33nmccMIJpKSksGLFCgBWrlzJmjVr6N+/P1OmTGHJkiVMnTqVW265hZ///OcAXHPNNbz55ptcdNFFzZ931llncfPNN1NUVERaWhpPPfUUN9xwAxdddBEPPvhgq+MjvfPOO7z66qssW7aM+Ph4SkqcmozLLruMm266CYC7776bJ554gltvvbVDfjRNIvNJIak/FK6Bjf8A4JRByQCsyC91MypjTCcwf/58Zs+eDcDs2bOZP38+ABMnTiQjIwOPx8O4cePYtm0bAB9++CGnnnoqY8eO5YMPPmDNmjWHfJ6IcM011/Dcc8+xf/9+li5dyvnnn3/UGN5//31uuOEG4uPjAUhJSQEgLy+P008/nbFjx/L8889/7V4dITKfFHr0gf4nw8aFcMZPyEiOo29SDLn5pVxzWqbb0Rlj4Jh/0YdCcXExH3zwAXl5eYgIjY2NiAgXXHABMTExzed5vV4aGhqoqanhhz/8Ibm5uQwcOJB58+a1Oqx205NBbGwsV1xxxTGn41RVRORr5ddffz2vvvoq2dnZPP3003z00Uft/s6Hi8wnBYCsGVDwBVSVICKMH5zMcntSMCaiLViwgGuvvZb8/Hy2bdvGjh07GDJkCIsXL271/KYE0Lt3byorK4/4tlH//v3p378/v/nNb7j++uuby6Oioqivr//a+TNmzODJJ59sbjNoqj6qqKggPT2d+vp6nn/++fZ81SOK7KSgftj8AeBUIRWUVrO3vPXJM4wx3d/8+fMPaRQGmDVrFn/9619bPb9Xr17cdNNNjB07lksuuYQJEyYc8bOvuuoqBg4cyKhRo5rL5syZw0knndTc0NzkvPPOY+bMmeTk5DBu3Ljm101//etfc+qpp3LOOecwcuTItn7No4rcobP9jXBvFgw/Gy57jJXbS7n04U95+KpTuGBsescGaowJSnceOvuWW27h5JNP5sYbbwzrfW3o7GB5vE5C2PQ++BsZ3b8nMT6PVSEZYzrc+PHjWb16NVdffbXboRxTKGdee1JECkUkr0VZiogsFJGNgXVyi2N3icgmEdkgIueGKq5DZM2AqmLYtZJon4fsjF6WFIwxHW758uX885//PKSxurMK5ZPC08B5h5XNBRapahawKLCPiIwCZgOjA9c83DRnc0gNOwvEc/DV1MHJrNlVRk19Y8hvbYxpXVeu0u5s2vKzDFlSUNV/AoePHXEx8Exg+xngkhblL6hqrapuBTYBE0MVW7P4FMiY0JwUxg9Opr5RbdIdY1wSGxtLcXGxJYYOoKoUFxcTGxt7XNeFu59CX1XdDaCqu0WkT6B8APBZi/MKAmVfIyJzgDkAgwYNan9EWefAB7+BykLGD3Zqs77YVsLEISnt/2xjzHHJyMigoKCAoqIit0PpFmJjY8nIyDiuazpL57Wv99KAVv9UUNXHgMfAefuo3XfOOtdJChv/QcrJV5PVpwfLtpZw85nt/mRjzHGKiopiyJAhbocR0cL99tFeEUkHCKwLA+UFwMAW52UAu8ISUb+xkDQANrwDwKShqeRuK6HeBsczxkSgcCeF14HrAtvXAa+1KJ8tIjEiMgTIAj4PS0QiMOJ8pxNbfTWThqZSVddI3k5rVzDGRJ5QvpI6H1gKjBCRAhG5EbgHOEdENgLnBPZR1TXAi8Ba4F3gZlUN3ytAIy6A+irY8jGnDnXaEj7bYvMrGGMiT8jaFFT120c4NP0I5/8W+G2o4jmqzKkQnQgb3qL3iPMY3qcHn20p5t+mDXMlHGOMcUvk9mhuyRcDWWfDhnfB72fS0BRyt5XYpDvGmIhjSaHJiG/CgULYuZxJQ1M5UNdI3q5yt6MyxpiwsqTQJOtsEC9seItTh6QC8NmWYpeDMsaY8LKk0CQuGTKnwIZ3SEuMYVhagiUFY0zEsaTQ0ohvQtF6KN7MpKGpfLHV2hWMMZHFkkJLIwLzpq5/q7ldYY21KxhjIoglhZaSB0O/k2Dd6839FZZaFZIxJoJYUjjcqIuh4Av6+Is5oW8Plmza53ZExhgTNpYUDjcqMJr3ujc4PSuNZVtLbH4FY0zEsKRwuN7Doc8oWPsaU7N6U9fg5/OtNuSFMSYyWFJozaiLYftSTk2rJ9rr4ZONNra7MSYyWFJozaiLASV+8zvkZCbzyUZrVzDGRAZLCq1JGwmpWbD2NU7PSmP9ngoKy2vcjsoYY0LOkkJrRJynhW1LmJbhTAq32N5CMsZEAEsKRzLqYtBGRuz/J6kJ0VaFZIyJCK4kBRG5XUTyRGSNiPwoUDZPRHaKyKrAcoEbsTXrNxZShuJZ8zJTs3rzycZ9+P3tnxLaGGM6s7AnBREZA9wETASygQtFJCtw+E+qOi6wvB3u2A4hAmNmwbZPOHsg7KusZf2eCldDMsaYUHPjSeFE4DNVrVLVBuBj4FIX4ji2MZeD+vlG/WIAezXVGNPtuZEU8oAzRCRVROKBC4CBgWO3iMhqEXlSRJJbu1hE5ohIrojkFhWF+Jd0n5HQdwxJm15jZL9EPtxQGNr7GWOMy8KeFFR1HfB7YCHwLvAl0AA8AgwDxgG7gf85wvWPqWqOquakpaWFPuAxs6DgCy7NbOCLbaWUVdeH/p7GGOMSVxqaVfUJVT1FVc8ASoCNqrpXVRtV1Q88jtPm4L4xswC40PspjX7l439ZFZIxpvty6+2jPoH1IOAyYL6IpLc45VKcaib3JQ+GjIn0L3ib1IRoFq3b63ZExhgTMj6X7vuSiKQC9cDNqloqIv8nIuMABbYB33cptq8bMwt596dcOayK5zcU0dDox+e1Lh7GmO7Hreqj01V1lKpmq+qiQNk1qjpWVU9S1ZmqutuN2Fo1+lIQD5dFfUpZdT3L80vdjsgYY0LC/twNRmJfGDqNobveItqrfLDe3kIyxnRPlhSClf1tPOU7uK7/Lt63dgVjTDdlSSFYIy+E6B5cEbWYzUUH2LbvgNsRGWNMh7OkEKzoeBh1CcOL3ieOGhZZFZIxphuypHA8xn0bT/0BrkvO4x9r9rgdjTHGdDhLCsdj0GToOYjZMUv4YlsJRRW1bkdkjDEdypLC8fB4IHs2g/d/TpqW8J49LRhjuhlLCscrezaCcmPS57yT13m6UhhjTEewpHC8UofBwEnM8nzM0s37KK60KiRjTPdhSaEtTrmW1Jp8xrOBf6y1PgvGmO4jqKQgIhkicmZgO0ZEEkIbVic3+hI0JokbEz7h7a+sCskY030cMymIyHeB14G/BIoGA6+FMqhOLzoBGXs50/2f8tXm7ZQeqHM7ImOM6RDBPCncBkwCygFU9V9An1AG1SWcci1R/loulCUstCokY0w3EUxSqFHV5j+FRcQLSOhC6iLSx6H9xnJt9Me8aVVIxphuIpiksERE7gRiA+0KfwPeDG1YXYAIcsp1nKBb2L9pGYUVNW5HZIwx7RZMUrgTqADWA7cDi4CfteemInK7iOSJyBoR+VGgLEVEForIxsA6uT33CIuxV+D3xXKl50NeX7XL7WiMMabdjpoUAlVFT6rqI6p6qapeEtj2t/WGIjIGuAlnDuZs4EIRyQLmAotUNQsn8cxt6z3CJq4XntGXcanvU95dsdHtaIwxpt2OmhRUtRFIF5GoDrznicBnqlqlqg3AxzhzMl8MPBM45xngkg68Z+hM/B7xVDOq8C3+tbfC7WiMMaZdgqk+2gJ8IiJ3ichtTUs77pkHnCEiqSISD1wADAT6Nk3BGVh3jTecBoynvt/JXON7n5eXF7gdjTHGtEswSaEIWAjEA2ktljZR1XXA7wOf+S7wJdAQ7PUiMkdEckUkt6ioqK1hdKioU28iS3ayY+V7+P3qdjjGGNNmouruLzER+R1QgNOIPU1Vd4tIOvCRqo442rU5OTmam5sbjjCPrr6auj+M4P2aEfS6bj6Th/d2OyJjjDkiEVmuqjmtHQumR/NCEfnH4Us7A+oTWA8CLgPm4/Savi5wynV0pV7TUXF4xl/LDE8uiz5f6XY0xhjTZr4gzrm7xXYsMAto79CgL4lIKlAP3KyqpSJyD/CiiNwIbAeuaOc9wsp36vfwf/YgqRvmU1EzjcTYjmybN8aY8DhmUlDVZYcVfSwiH7fnpqp6eitlxcD09nyuq5IzqRh4JldsX8gby7fynSknuB2RMcYct2Cqj5JaLL1EZDqQHobYupykabeRJuXsXfIcbrfVGGNMWwRTfbQGUJzxjhqArTidz8xhZOg0SntkcV75y6zecQfZgzp/p2xjjGkpmFdSh6rqIFUdqKpDVPUsYEmoA+uSRIg941ZO9Ozgiw9fcTsaY4w5bsEkhcPbFAA+7+hAuou4U2ZT4Usha8uzVNTUux2OMcYclyMmBRHpIyLZQJyIjBWRkwLLVJyObKY1vhgqT7qeb8hKPly82O1ojDHmuBytTeGbwHeBDODhFuUVwH+FMqiurt9ZP6RuxZ/xfv4oOn0aIjb9hDGmazhiUlDVp4CnRORbqvpiGGPq8qRHGjsGzmT69tdZvmYDOWNGuh2SMcYE5ZhtCqr6ooicKyI/FpH/bFrCEVxXlvHNnxItDexdeJ/boRhjTNCC6afwMM6wEz8G4oCrgeEhjqvLi+k3go2p0zl9/6tsLbAJeIwxXUMwbx9NVdXvAMWq+l/AqTjtDOYY+pw/lySpZuNbD7gdijHGBCWYpNA0+XCNiPQL7GeGLKJuJHn4BDb0mMDJu/7K/rIyt8MxxphjCiYpvC0ivYB7gVXANmBBKIPqTuLO/AlpUsaqNx8+9snGGOOyY83R7AHeUdX9qvp3YAgwVlWtoTlIg06ZwcboEzlh45PU1NQc+wJjjHHRseZo9gP3t9ivVtWSkEfVnYjQMPnH9KeQFW8+6nY0xhhzVMFUHy0UkYtDHkk3NvKMy9nsy2Jw3sPU1trTgjGm8womKdwCvCIi1SJSIiKlImJPC8dBPB5qpvyEAexlxRv2tGCM6byCSQq9gSigB5AW2E9rz01F5A4RWSMieSIyX0RiRWSeiOwUkVWB5YL23KOzGfWNK9jky2LQmoepq23vxHXGGBMawfRobsSZGvOnge10YFxbbygiA4DbgBxVHQN4gdmBw39S1XGB5e223qMzEo+Hqsk/YYDuZaW1LRhjOqlgejQ/CJwJXBMoqgLa+1vNhzP6qg9nxNWI6PI7dtoVbPRlMTDvIeqsbcEY0wkFU300WVW/T6ATW+Dto+i23lBVd+L0edgO7AbKVPUfgcO3iMhqEXlSRFqdtkxE5ohIrojkFhUVtTUMV4jHQ/WUn9Jf97LilfuPfYExxoRZMEmhPtBfQQFEJBXwt/WGgV/2F+P0eegPJIjI1cAjwDCcqqndwP+0dr2qPqaqOaqak5bWrqYNV4z9xizWRY9l+PqHqaywXs7GmM4lmKTwEPASkCYivwQWA79vxz3PBraqapGq1gMv4zyN7FXVxkDfiMeBie24R6clHg++Gb+kN/tZ9ff/53Y4xhhziGAamp8F7sap8ikBrlDVF9pxz+3AJBGJF2f2menAOhFJb3HOpUBeO+7RqWXlTOfLhClk5z9NUWFENKcYY7qIYJ4UwHlDqB6oO45rWqWqy3DGTloBfBX4vMeAP4jIVyKyGqdh+4723KezS535GxKoYcPff+V2KMYY0yyYt49+BszHqf/PAP4qIne156aq+gtVHamqY1T1GlWtDazHqupJqjpTVXe35x6dXcaIU1iVch4TChewbfN6t8MxxhgguL/6rwYmqOrdqvoznLr+a0MbVmTIvOK3KLD7pbluh2KMMUBwSSGfQ+dy9gFbQhNOZEnpP4x1Q67jtKoPyV38ntvhGGNMUEmhClgjIn8Rkcdx2gH2i8gfReSPoQ2v+xv9rV+wT5KJ++C/qKtvdDscY0yECyYpvAXMA5YCnwG/Aj4A1gQW0w7R8UkUTbyT0f4NfPKKDX9hjHGX71gnqOoT4Qgkkp147g/IX/kEo9b8D4Ul36FPSquduY0xJuSCefvoPBH5QkQKbejsEPF4iP7mH0iXYnKf/4Xb0RhjIlgw1UcPAt8HBtBBQ2ebr0vPns76tHOZvu+vfL481+1wjDERKpikUACsUtX6wDAUjYEhtE0HG/KdP9IgPvxv/YSauga3wzHGRKBgksKdwBsi8hMRua1pCXVgkSgmOYPC8T9mkn8F77xkTTnGmPALJin8EmgEeuFUGzUtJgSGXHAHu2OGMGH9H/jXjkK3wzHGRJhjvn0E9FHV8SGPxDi8UcRfch89/3YxL/11LkP/4y/4vO0absoYY4IWzG+bRSJyVsgjMc16njiN7ZlXcHHVK7z01ltuh2OMiSDBJIWbgPdFpNJeSQ2fQVf+N1W+noxZfjcbdpW6HY4xJkIEkxR6A1FAT+yV1PCJS0YuuJfRso0lz/2K+sY2T3ZnjDFBC2aSnUbgCuCnge10nCkzTYglnjKLvf2n850D/8f/vfWh2+EYYyJAMD2aH8SZ9OaaQFEV0K5BekTkDhFZIyJ5IjJfRGJFJEVEForIxsDaxnoQoe/sB1FvDGNz7+KLLUVuR2SM6eaCqT6arKrfB2oAVLUEiG7rDUVkAHAbkKOqY3BmdZsNzAUWqWoWsCiwb5L6Ixf8ngmeDXz2199QVlXvdkTGmG4smKRQLyIeQAFEJBVobwW3D4gTER8QD+wCLgaeCRx/BriknffoNmLHX0XZoBnMqX+eB/72BqrqdkjGmG7qiEkh8Asb4CHgJSBNRH4JLAZ+39YbqupO4F5gO7AbKFPVfwB9m6bgDKz7HCGuOSKSKyK5RUURUp0iQs9vPYQ/Kp6ZW3/N35bZHEfGmNA42pPC5wCq+ixwN84v8lLgClV9oa03DLQVXAwMwZn3OUFErg72elV9TFVzVDUnLS2CXoLq0YfoSx4g27OFfW//hrydZW5HZIzpho6WFKRpQ1XXqOr9qnqfqua1855nA1tVtUhV64GXgcnAXhFJBwisbYyHw3jHXELt6G/xb55XefTZZ619wRjT4Y42zEWaiPz4SAdVta1TcW4HJolIPFANTAdygQPAdcA9gfVrbfz8bi1m5h+p2fE5Pyv7I/81fxT33XAWHo8c+0JjjAnC0Z4UvEAPIPEIS5uo6jJgAbACZ75nD/AYTjI4R0Q2AucE9s3hYhKJvfJJ+nrKOX/b/+PPiza6HZExphuRI73JIiIrVPWUMMdzXHJycjQ3NzInpNHF9yPv/5yf1X+XqbPv5Pyx6W6HZIzpIkRkuarmtHYsqDYF0/nI5FtpHHoW86Ke5S8vvmINz8aYDnG0pDA9bFGY4+fx4J31FzyJffmz90/8+zMfUlhe43ZUxpgu7ohJIdBz2XRmCal4r3yWfp5S5tbez/ee/pwDtTaNpzGm7Wz2lq4uIwfPub/jTFnBtMJn+cFzy6lrsBFVjTFtY0mhO5h4E4z9Fj/2/Z2Yze8x96XVNhSGMaZNLCl0ByIw8wHofzIPxz3CV6uW8bu311liMMYcN0sK3UVUHFz5PFGxPfhb4n38/ZPV3Pe+9WEwxhwfSwrdSc8ByOznSW7cx4KUR3l40Toe+Wiz21EZY7oQSwrdzcCJyMwHGV61kuf7PMfv313HE4u3uh2VMaaLONrYR6aryr4S9ucz8cPf8lB6P25+U6hr8PNv04a5HZkxppOzpNBdnfETKM3nm6ueZX9mOj97F6rrG7nj7CxErLO6MaZ1lhS6KxG46D4o38l3tt5LXdav+eUiqK5r4K7zT7SRVY0xrbI2he7MGwVX/h+Sns31u37Jz8eU8PgnW/nxi6usg5sxplWWFLq7mES4agGSnMkNO/6T309WXl21i+8+/QUVNTZJjzHmUJYUIkFCKlzzChLbkyvX38Zj58bz2ZZirnh0KQWlVW5HZ4zpRMKeFERkhIisarGUi8iPRGSeiOxsUX5BuGPr1noOgGtfA280M3JvYv6lyezcX80lDy0hd5uNfWiMcYQ9KajqBlUdp6rjgPFAFfBK4PCfmo6p6tvhjq3bSx0G170B4mHCx9fx1nf6khgbxbcf/4wXv9jhdnTGmE7A7eqj6cBmVc13OY7I0TsLrn0d/I0Mev1KXr8yjUlDU7nzpdX8dMFqauob3Y7QGOMit5PCbGB+i/1bRGS1iDwpIsluBdXt9RkJ178JKInzZ/L0+XHccuZw/pa7g0sf/pSt+w64HaExxiWuJQURiQZmAn8PFD0CDAPGAbuB/znCdXNEJFdEcouKisISa7fU50S44R3wxeJ99kL+Y3QlT10/gd1l1Vz058W8tLzARlk1JgK5+aRwPrBCVfcCqOpeVW1UVT/wODCxtYtU9TFVzVHVnLS0tDCG2w2lDoMb3oa4ZHh2JmdG5fHWbaczKj2Jf//7l9wyfyX7q+rcjtIYE0ZuJoVv06LqSETSWxy7FMgLe0SRKHkw3PAuJGfC899iwPY3mT9nEneeN4L38vZw7n3/ZNG6vW5HaYwJE1eSgojEA+cAL7co/oOIfCUiq4EzgTvciC0iJaXD9W/BwFPh5e/h/ewhfviNYbzywyn0iovmxmdyuf2FlRRX1rodqTEmxKQr1xvn5ORobm6u22F0H/U18MocWPsaTPgenPd76tTDQx9u4uGPNpEYG8Xc80dy+SkZNnaSMV2YiCxX1ZzWjrn99pHpTKJi4fKnYcrt8MVfYP6VRDdUcMc5J/DGrVPJTI3nzgWrmfXop+TtLHM7WmNMCFhSMIfyeOCcX8FFD8CWj+CJc6F4MyP7JbHgB5P578tPYntxFRc9uJg7F3zJnrIatyM2xnQgqz4yR7blY/j7daB+mPUEZJ0DQFl1PX9etJFnl+bj8cD3pg5lzjeGkhQb5XLAxphgHK36yJKCObrSbfDC1bA3D866G6b+2HmaALYXV/GH99bz5urd9IyL4qbTh3D9lCH0iLFpOozpzCwpmPapq4I3boOv/g5Z58Klj0J8SvPhvJ1l3Pf+v3h/XSHJ8VFcP3kI1542mOSEaILS5jMAABJOSURBVBeDNsYciSUF036qTuPze/8JCX3giqdg4KH9C7/csZ8HFm1k0fpC4qK8zJ44kO9OGcLAlHiXgjbGtMaSguk4u1bCi9dBWQGceVegOsl7yCkb9lTwv//czOurdtGoyvSRfbj2tExOz+pt80Mb0wlYUjAdq3o/vHkHrHkZBk6Cy/7X6RF9mD1lNTy/LJ/5n29nX2Udo9KTuP3sLGaM6mvJwRgXWVIwHU8VVr8Ib/+Hs33e7+Dka6CVX/a1DY28vmoXD324iW3FVYxKT+LWs4Zz7uh+1gnOGBdYUjChU5oPr/4Q8hfDsOkw8wHomdHqqQ2Nfl5btYs/f7CRbcVVDE1L4AffGMYl4wYQ7bMuM8aEiyUFE1p+P+Q+AQt/AeKBc+bB+O82v7p6uIZGP2/n7eGRjzazbnc5/ZJiuXbyYL4zcRC94u2NJWNCzZKCCY/SbfD6bbD1Y8iYABfdD31HH/F0VeWjfxXxxCdbWbxpH3FRXi47ZQDXnpbJiH6J4YvbmAhjScGET1Nbw3t3OQ3Sk/4NvvFTiE066mXrdpfz5OKtvPblLuoa/EwcksI1kwYzY3RfYnzeo15rjDk+lhRM+FWVwPvzYMWzkJAGZ8+D7G8fsUqpScmBOl7M3cFzn+VTUFpNSkI0l508gNkTBzK8jz09GNMRLCkY9+xcAe/cCQVfQP+TYcZvIHPqMS9r9CuLN+3jhc+3s3DtXhr8SnZGT2aNz+Cik/pbb2lj2sGSgnGX3+8MkbHoV1BeACMucJ4c0kYEdXlRRS2vrdrJSyt2sm53OVFe4YysNGaO68/ZJ/YlwcZaMua4dKqkICIjgL+1KBoK/Bx4NlCeCWwDvqWqpUf7LEsKXUx9NXz2CHzyR6g/ACdd6bQ3pAwJ+iPW7irntVU7eePLXewqqyE2ysO0E/pw/th+nDWyD4k2Uqsxx9SpksIhNxfxAjuBU4GbgRJVvUdE5gLJqvrTo11vSaGLOlAMS/4Enz8O/gY4+WqYekervaKPxO9XcvNLeXP1Lt7N20NhRS3RXg+nDUvlnFF9OWdUX/omxYbuOxjThXXmpDAD+IWqThGRDcA0Vd0tIunAR6p61PoFSwpdXPlu+ORepzHa3wjZs+H0f4fUYcf1MX6/smJ7Ke/m7WHhur3kF1cBMGZAEmeO6MO0EX3IzuiJz2sd5IyBzp0UngRWqOqDIrJfVXu1OFaqqsmtXDMHmAMwaNCg8fn5+eEL2IRG+S5Ycj8sfxoaauHEi2DKjyBj/HF/lKqysbCShWv38tGGQpbnl+JXSIr1MWV4b6Zm9Wbq8N4MSom38ZdMxOqUSUFEooFdwGhV3RtsUmjJnhS6mYq98Pn/OkN015TBoMkw6Qcw4pvgbVtj8v6qOj7ZuI/FG/fxz41F7A5MHzqgVxyThqYyaWgKE4ekWJIwEaWzJoWLgZtVdUZg36qPjKO2wqlSWvYo7N8OPQfBhBudAfcSUtv8sarK5qIDLN28j6Vbilm6uZjSqnoA+ibFkJOZwvhByYwfnMyo/klEWXWT6aY6a1J4AXhPVZ8K7P83UNyioTlFVe882mdYUujm/I2w4W347FFnwD1vNIy6BHJugEGntToi63F9vF/ZVFTJ51tLWLa1hBX5pezcXw1AjM/D6P5JjBuYTPbAnowZ0JMhqQk2qqvpFjpdUhCReGAHMFRVywJlqcCLwCBgO3CFqpYc7XMsKUSQwnWQ+xR8OR9qyyFlmPPWUva3ISm9w26zu6ya5fmlrNq+ny8L9vPVzjJq6v0A9IjxMap/EqP7JzEqPYlR/ZMY3qeHDcNhupxOlxQ6iiWFCFR3ANa+Biufg/wlzqisQ77hvLk08kKI6dGht2to9LOxsJKvCspYvXM/a3aVs353BdX1jQB4PcLQ3gmc0C+RE/okckLfHmT17cHg1ASrfjKdliUF0z0Vb3aeHFb/zWl7iIqHE86F0ZdC1gyIigvJbRv9ytZ9B1i/x0kQ6/eUs2FvBTtKqpvP8XmEQanxDEvrwdDeCQzpnUBm7wQyUxPokxhj1VDGVZYUTPemCjuWOaOzrn0NqvZBVAJkneO83po145ijtHaEqroGNhVWsnFvJVv2VbKpsJLNRQfYXlxFXaO/+bzYKA+DUxIYmBJHRnI8A1PiyUiOc5Ze8STF+exNKBNSlhRM5GhscBql17wC69+GA4VOA3Xm6TDifDjhPOg1MLwh+ZVd+6vZuu8A+SVV5O87wLbiKgpKq9hRUsWBusZDzu8R4yO9Zyz9e8WR3jOWfj1j6ZcUS9+esfRNdPaT46MscZg2s6RgIpO/EXZ8DuvfhA3vQMlmpzztRBg+HYaf7bzFFOXecBiqSsmBOnbur2ZnaTUFpdXsKqtm1/5qdu6vZk9ZLfsqa792XZRXSOsRQ1qis/QObKcmRJPaw1mn9IgmJSGa5Phoa98wh7CkYAzAvk3wr3dh00LI/xQa68AXC4MmwdBpMOQM6Jfd5o5yoVLX4GdveQ2FFTXsLa9lT1kNRZW1FJbXUlhRw77KOvZV1lJcWYv/CP87J8b6SI6PJjk+il4t1j3jog5ZkuKiSIrzkRgbRWKsjx7RPmv/6IYsKRhzuLoDsPUT2PKRM31o4VqnPDoRBp8Ggyc7TxH9TwZfjKuhBqvRr5RV11NcWUtRZS2lB+opOVBLyYF6SqvqAks9pQfqKKuuZ39VHeU1DUf9TBHoEe1zEkSsjx4xPhJiDl97iY92tuOivSRE+4iP9hIX7SUuykt8tJfYwDou2kusz2uJxmWWFIw5loq9TlvEtsVOsije6JR7Y5zEMHCCM+90xgRI6u9urB2o0a9U1NRTVu0sFTUNlLfYrqipp7ymgcraBiqb1oHlQNNS10jjkR5RjiDG5yE2yktsVGDtc7ZjfF5iDlk3LV6ifR6ivR5n3cp2lNdDlFeIarHv8wrRgbXP4zm47RWiPAfLnbVETDuNJQVjjldlkfNG0/alTrvE7lVOdRNAj34w4BTofwqkZztLYl9343WRqlLb4KeqrpGqugYO1Drr6rpGquoaqa4PLHUH1zX1TYufmgZnu7bB31xW1+CntiGw3Xhwv67Bf8Qqso7g8whej5MgfF4P3hb7LdfO4sHrAa8cLPO0su2s+VpZU7lHnGTkkYPniNB8btO2J7AWQEQY0S+RC8a2rePm0ZJC56o8Naaz6JEGJ17oLOCM3rp7NexcDrtWONOMbni7xfl9oe8Y6DcG+o6FvqMgNQt83X/aUBEJ/NXvJSUM06Q2NB5MFE6y8NPg1+b9er+fhkZt3q4PHK9v9FPfqDQ0+qn3O+uGRqXBf7CsMXBtfaPiV6XB76fRrzQ0Ko2qB7f9B/cPWdS5b4NfUW06xxlSpVH14FoVvx9nHThHm7cVVedYozrbqgS2tTkpXnhSepuTwtFYUjAmGL4Ypwpp4ISDZTXlsDcPdq2CPV/B3q9g6cPgdwbZw+OD1OHOtKNpI6H3CdA7yymLTnDne3QDPq8Hn9dDfPfPt0ekgWQRCpYUjGmr2CSnQXrw5INlDXVOe0ThusCyFvbkwbo3QA92YCNpAKQMdSYUShnmTEmaPMSZfa6Dh+ow3Y8EqphCwZKCMR3JFw19RztLS/U1Tj+JfRudpLFvk7O/7g2oKj703Pje0GsQJA+GngOdpddAJ5H0zIC45HaPEGvMkVhSMCYcomJbTxYA1fuhdCuUboOSrbA/H0rzYfeXsP6tgw3czZ8VD4npzltQienOKLGJ6U67RmI/Z92jrz1xmDaxpGCM2+J6QdzJzquvh/P74UARlO2AsgIo3wllO6FilzPH9Y7PoGLP1xMHOMkjoTck9IGENGeCooQ050kkPjWwpDhPHvEpENMTPNbzOdJZUjCmM/N4nNddE/tCRqtvEDqvplSXQsVuqNwLlYVOojhQ5CyVhU5C2b0KDuw72BB+OPFAbE8nScT2cpJV87qns8QkBdaJznZM4qGLNyp0PwsTFq4kBRHpBfwFGAMo8F3gXOAmoChw2n+q6tutf4IxppmI85d+fErr1VMtqTrzX1eXQFWJ055RVXJwv2a/k2CqS53zSvOdspryIyeTlrwxTrVVdA8nSUQnOEtUfCvreGc026i4QxdfnFPd1rwOLFFxzuCG1p4SUm49KdwPvKuql4tINBCPkxT+pKr3uhSTMd2fSKC6qpfz9lOwVKGhxkkUNeXO7Hc1Zc582nWVzrq20imvq3SGEak7ECivcJ5c6iqhrgrqq6H+QNu/gy/WST6+wOKNDiSO6EB5YO2NDmxHO08w3hbbnqZ9n7P2RAXKfS2O+5x9T1P5Yfseb2Dta7HvBfEeWiaeFtveg2WdNLmFPSmISBJwBnA9gKrWAXWR0r3cmC5J5OBf8on92v95qoHkUA31VYGl2kk89VXO21oN1YF1jdN5sKE6sK5xyhtrnVeAG2qcNpWG2oNltRXO2l8fKK+DxvrAEij3H33cp5ATT4sk0ZQ8PC22W5Q3H2uxZM2Ac3/b4WG58aQwFKeK6CkRyQaWA7cHjt0iItcCucC/q2rp4ReLyBxgDsCgQYPCE7ExpmOJONVH0fFAqjsx+P1OcmisD6wbWuw3OEvTMX9ji/L6wLUNB4/5G5y1Nh681t948Jg2tjjuP2y/0enDov6DZYdsc/A89OC5SQNC8mMJ+9hHIpIDfAZMUdVlInI/UA48COzD+RH8GkhX1e8e7bNs7CNjjDl+Rxv7yI33zwqAAlVdFthfAJyiqntVtVFV/cDjwEQXYjPGmIgW9qSgqnuAHSIyIlA0HVgrIi1HdroUyAt3bMYYE+ncevvoVuD5wJtHW4AbgAdEZBxO9dE24PsuxWaMMRHLlaSgqquAw+uzrnEjFmOMMQdZn3ZjjDHNLCkYY4xpZknBGGNMM0sKxhhjmoW981pHEpEiIN/tONqgN05HvUgTid87Er8zROb37krfebCqprV2oEsnha5KRHKP1JuwO4vE7x2J3xki83t3l+9s1UfGGGOaWVIwxhjTzJKCOx5zOwCXROL3jsTvDJH5vbvFd7Y2BWOMMc3sScEYY0wzSwrGGGOaWVIIMxE5T0Q2iMgmEZnrdjyhJiIDReRDEVknImtE5PZjX9V9iIhXRFaKyJtuxxIOItJLRBaIyPrAf/PT3I4pHETkjsC/7zwRmS8isW7H1FaWFMJIRLzAQ8D5wCjg2yIyyt2oQq4BZ2rVE4FJwM0R8J1buh1Y53YQYXQ/8K6qjgSyiYDvLiIDgNuAHFUdA3iB2e5G1XaWFMJrIrBJVbeoah3wAnCxyzGFlKruVtUVge0KnF8SoZlctpMRkQzgm8Bf3I4lHEQkCTgDeAJAVetUdb+7UYWND4gTER8QD+xyOZ42s6QQXgOAHS32C4iQX5AAIpIJnAwsO/qZ3cZ9wJ2A3+1AwmQoUAQ8Fagy+4uIJLgdVKip6k7gXmA7sBsoU9V/uBtV21lSCC9ppSwi3gkWkR7AS8CPVLXc7XhCTUQuBApVdbnbsYSRDzgFeERVTwYOAJHQbpaM88Q/BOgPJIjI1e5G1XaWFMKrABjYYj+DLvyYGSwRicJJCM+r6stuxxMmU4CZIrINp5rwLBF5zt2QQq4AKFDVpifBBThJors7G9iqqkWqWg+8DEx2OaY2s6QQXl8AWSIyJDA/9WzgdZdjCikREZw65nWq+ke34wkXVb1LVTNUNRPnv/MHqtpl/3oMhqruAXaIyIhA0XRgrYshhct2YJKIxAf+vU+nCzewuzJHc6RS1QYRuQV4D+cNhSdVdY3LYYXaFJz5t78SkVWBsv9U1bddjMmEzq3A84E/erYAN7gcT8ip6jIRWQCswHnbbiVdeMgLG+bCGGNMM6s+MsYY08ySgjHGmGaWFIwxxjSzpGCMMaaZJQVjjDHN7JVUY4IgIqnAosBuP6ARZ0gHgCpV7bKdlYxpyV5JNeY4icg8oFJV73U7FmM6mlUfGdNOIlIZWE8TkY9F5EUR+ZeI3CMiV4nI5yLylYgMC5yXJiIvicgXgWWKu9/AmIMsKRjTsbJx5lAYi9OT+wRVnYgzfPatgXPuB/6kqhOAWUTI0Nqma7A2BWM61hequhtARDYDTUMofwWcGdg+GxjlDJMDQJKIJAbmmzDGVZYUjOlYtS22/S32/Rz8/80DnKaq1eEMzJhgWPWRMeH3D+CWph0RGediLMYcwpKCMeF3G5AjIqtFZC3wA7cDMqaJvZJqjDGmmT0pGGOMaWZJwRhjTDNLCsYYY5pZUjDGGNPMkoIxxphmlhSMMcY0s6RgjDGm2f8H0D3ocarYMFgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#b \n",
    "    #modifying the original euler approximation to get a more accurate model using the weather\n",
    "\n",
    "import numpy as np\n",
    "from scipy import interpolate\n",
    "\n",
    "#interpolate from the given data values the temperature at any give time using a linear interpolation\n",
    "\n",
    "def amb_temp(hours):\n",
    "    time_values=[-3,-2,-1,0,1,2]\n",
    "    temp_values=[55,58,60,65,66,67]\n",
    "    \n",
    "    return np.interp(hours,plot_time,plot_temp)\n",
    "\n",
    "\n",
    "\n",
    "def euler_improved_solution(hours,n):\n",
    " \n",
    "    t=np.linspace(0,hours,n)\n",
    "    temp=np.zeros(len(t))\n",
    "    temp[0]=85\n",
    "\n",
    "    #simply replace the original \"65\" with the function amb_temp\n",
    "    for i in range(0,len(t)-1):\n",
    "        temp[i+1]=temp[i]-(0.611*(temp[i]-amb_temp(hours)))*hours/n\n",
    "            \n",
    "    return temp[-1]\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "#solve the problem analytically with the given data and formulas, and calculated K value \n",
    "#given the original temp, ambient temp, and time elapsed, it returns the body temperature\n",
    "def analytical_solution(temp_original,temp_ambient,hours):\n",
    "    exponent=-0.611111*hours\n",
    "    analytical_solution=temp_ambient+(temp_original-temp_ambient)*np.exp(exponent)\n",
    "    \n",
    "    return analytical_solution  \n",
    "\n",
    "\n",
    "#here i create two lists of values from the improved euler approximation and the original analytical approximation\n",
    "#i also make a time list so that the values can be plotted together as seen below\n",
    "euler=[]\n",
    "analytical=[]\n",
    "time=[]\n",
    "\n",
    "#i multiply the range by 10 and divide i values by 10 to generate more values and thus a smoother graph\n",
    "for i in range(-10,90):\n",
    "    euler.append(euler_improved_solution(i/10,100))\n",
    "    analytical.append(analytical_solution(85,65,i/10))\n",
    "    time.append(i/10)\n",
    "    # the below print command is just to see what is produced, and ensure correct values are generated\n",
    "#print(euler,'\\n\\n',analytical,'\\n\\n',time,'\\n')\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(time,euler,label='Euler')\n",
    "plt.plot(time,analytical,label='Analytical')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Temperature')\n",
    "plt.legend();\n",
    "\n",
    "#the graph shows the temperature prediction based on the different models in hours time relative to 11am.\n",
    "\n",
    "\n",
    "#a simple guess and check method using the updated euler function produces a new time of death value\n",
    "euler_improved_solution(-0.743,1000)\n",
    "\n",
    "print('The improved Euler model predicts a time of death of 10:15am. This is about 6 minutes different', \n",
    "'from the original linear analytical model.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}