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": 352,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K is equal to 0.275\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "T2=74\n",
    "T1=85\n",
    "deltat=2\n",
    "Ta=65\n",
    "K = -(T2-T1)/((T1-Ta)*deltat)\n",
    "\n",
    "print(\"K is equal to %0.3f\" %K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": 331,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The value for K using the initial conditions mentioned above is 0.275\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "def measure_K(Temp_t1,Temp_t2,Temp_ambient,delta_t):\n",
    "    return -(Temp_t2-Temp_t1)/((Temp_t1-Temp_ambient)*delta_t)\n",
    "    \n",
    "print(\"The value for K using the initial conditions mentioned above is %0.3f\" %measure_K(85,74,65,2))    \n",
    "    \n",
    "    "
   ]
  },
  {
   "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": 346,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "answer to part a)\n",
      "Eueler integration with only 5 time steps\n",
      "[85.         73.76469985 68.84099817 66.68325981 65.73766335 65.32326989]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "print(\"answer to part a)\")\n",
    "\n",
    "print(\"Eueler integration with only 5 time steps\")\n",
    "t=np.linspace(0,15,6)\n",
    "T_num=np.linspace(85,85,6)\n",
    "for i in range(1,len(t)):\n",
    "    K=0.275\n",
    "    T_num[i]=65+(T_num[i-1]-65)*np.exp(-K*3)\n",
    "print(T_num)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 347,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "answer to part a)\n",
      "Eueler integration with 100 time steps\n",
      "[85.         83.92970296 82.91668271 81.95787408 81.05037596 80.19144246\n",
      " 79.37847467 78.60901272 77.88072842 77.19141815 76.53899621 75.92148853\n",
      " 75.33702669 74.78384224 74.26026137 73.76469985 73.29565823 72.85171731\n",
      " 72.43153382 72.03383639 71.65742167 71.30115074 70.96394559 70.64478592\n",
      " 70.34270604 70.05679192 69.78617844 69.53004681 69.28762203 69.05817057\n",
      " 68.84099817 68.63544772 68.44089728 68.25675817 68.08247324 67.91751514\n",
      " 67.76138475 67.61360965 67.47374272 67.34136074 67.21606317 67.09747087\n",
      " 66.98522503 66.87898601 66.77843235 66.68325981 66.59318041 66.50792159\n",
      " 66.42722539 66.35084764 66.27855722 66.21013542 66.14537521 66.08408062\n",
      " 66.02606621 65.97115643 65.91918513 65.86999508 65.82343742 65.77937129\n",
      " 65.73766335 65.6981874  65.66082401 65.62546011 65.5919887  65.56030852\n",
      " 65.53032369 65.50194349 65.47508206 65.44965812 65.42559473 65.40281909\n",
      " 65.38126229 65.36085909 65.34154777 65.32326989 65.30597015 65.2895962\n",
      " 65.27409851 65.25943016 65.2455468  65.2324064  65.2199692  65.20819758\n",
      " 65.19705592 65.1865105  65.17652942 65.16708247 65.15814108 65.14967818\n",
      " 65.14166818 65.13408683 65.12691119 65.12011956 65.11369138 65.1076072\n",
      " 65.10184862 65.0963982  65.09123947 65.0863568  65.08173543]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "print(\"answer to part a)\")\n",
    "\n",
    "print(\"Eueler integration with 100 time steps\")\n",
    "t=np.linspace(0,20,101)\n",
    "T_num=np.linspace(85,85,101)\n",
    "K=0.275\n",
    "for i in range(1,len(t)):\n",
    "    T_num[i]=65+(T_num[i-1]-65)*np.exp(-K*0.2)\n",
    "print(T_num)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "part a) The difference by increasing the time step from 5 to 100 time steps is roughly 0.25 degrees. It is evident that the temperature converges to 65 degrees.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 329,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "As the time approaches infinity, the Temperature approaches 65.00000\n"
     ]
    }
   ],
   "source": [
    "t=np.linspace(0,100,501)\n",
    "T_num=np.linspace(85,85,501)\n",
    "for i in range(1,len(t)):\n",
    "    K=0.275\n",
    "    T_num[i]=65+(T_num[i-1]-65)*np.exp(-K*0.2)\n",
    "\n",
    "\n",
    "print(\"As the time approaches infinity, the Temperature approaches %0.5f\" %T_num[500])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer to part b)\n",
    "\n",
    "As demonstrated in the numerical solution above, the temperature of the corpse is 65 degrees farenheit as t approaches infinity. This is also evident from the analytical equation. As t approaches infinity, e^-x approaches zero so the second term in the equation approaches 0. As a result, the final temperature approaches the ambient temperature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Part c) at what time was the time of death?\n",
    "\n",
    "By keeping the same analytical solution, but changing the initial temperature of time 0 to 98.6 degrees farenheit, the time taken for the body to reach 85 degrees from 98.6 can be solved by counting the time steps it takes to reach that temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K in terems of minutes is equal to 0.005\n",
      "the corpse reaches 85 degrees in roughly 1.88333 hours\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "#solving for a new constant K in minutes for a more accurate answer\n",
    "\n",
    "import numpy as np\n",
    "T=np.array([85,74])\n",
    "k = -(T[1]-T[0])/((T[0]-65)*2*60)\n",
    "\n",
    "print(\"K in terems of minutes is equal to %0.3f\" %k)\n",
    "\n",
    "# creating a while loop to determine when the time reaches 85 degrees\n",
    "\n",
    "T_ambient= 65.00\n",
    "corpse = 98.60\n",
    "minutes = 0.00\n",
    "temperature_loss=0\n",
    "\n",
    "while (corpse>85.00):\n",
    "    temperature_loss= k*(corpse-T_ambient)\n",
    "    corpse-= temperature_loss\n",
    "    minutes+=1.00\n",
    "    \n",
    "#returning the answer back into hours\n",
    "\n",
    "print(\"The corpse reaches 85 degrees in roughly %0.5f hours\" %(minutes/60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer to part c) \n",
    "If the corpse reached a temperature of 85 degrees at 11 am. Then the time of death was approximately 9:07 am."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Now that we have a working numerical model, we can look at the results if the\n",
    "ambient temperature is not constant i.e. T_a=f(t). We can use the weather to improve our estimate for time of death. Consider the following Temperature for the day in question. \n",
    "\n",
    "    |time| Temp ($^o$F)|\n",
    "    |---|---|\n",
    "    |8am|55|\n",
    "    |9am|58|\n",
    "    |10am|60|\n",
    "    |11am|65|\n",
    "    |noon|66|\n",
    "    |1pm|67|\n",
    "\n",
    "    a. Create a function that returns the current temperature based upon the time (0 hours=11am, 65$^{o}$F) \n",
    "    *Plot the function $T_a$ vs time. Does it look correct? Is there a better way to get $T_a(t)$?\n",
    "\n",
    "    b. Modify the Euler approximation solution to account for changes in temperature at each hour. \n",
    "    Compare the new nonlinear Euler approximation to the linear analytical model. \n",
    "    At what time was the corpse 98.6$^{o}$F? i.e. what was the time of death? \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 336,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "part a)\n",
      "Using the function, the ambient temperature at time -1  is 58.00\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Temperature')"
      ]
     },
     "execution_count": 336,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEaCAYAAAAsQ0GGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU5fX48c9JQhZCIJAQ9p2wuICBKC4sAipYxSIq4kpr1VorLlVc2qr9ys+vWHHX+lXrUveqIC5UqWwKiiiL4AIh7BAEspBAQvac3x93MiQwEzLJTGaSnPfrldfk3vvcOyeB5OTe53nOI6qKMcYY429hwQ7AGGNM02QJxhhjTEBYgjHGGBMQlmCMMcYEhCUYY4wxARER7ABCRWJiovbs2TPYYRhjTKOyatWqLFVt7+mYJRiXnj17snLlymCHYYwxjYqIbPd2zB6RGWOMCQhLMMYYYwLCEowxxpiAsARjjDEmICzBGGOMCQgbRWaMMc3U3DUZPDI/jd25hXSOj2H6uP5MTOnit+tbgjHGmGZo7poM7p6zjqLSCgAycgu5Z84PAH5LMvaIzBhjmglVZdO+fF79ait3zz6cXCoVlpbzyPw0v72f3cEYY0wTlp1fzFebs1m6MZNlm7L4Ja+oxva7cwv99t6WYIwxpgkpKi1n1fb9LE3PYml6Jj/tPuDT+Z3jY/wWiyUYY4xpxFSVDXsOsiw9iy/TM/luW85Rj76qiouK4LQ+CbSOjuDjdb9QXHa4bUyLcKaP6++32CzBGGNMI7PvQBFL07NYtimLpelZZOUXe20bHiac1C2eEcmJjEhOZHDXeCLCne734cntbRSZMcY0Z4dKylixNYdl6VksS88ibe/BGtv3SoxlRHIiw/smcmqfBFpHt/DYbmJKF78mlCNZgjHGmBBTUaH8uDvPuUtJz2LV9v2UlHt/7BXfsgVn9HHuUIYnJ9K1bcsGjNY7SzDGGBMCdu0/xLL0LJZuyuLrTVnsP1TqtW2LcGFoj7aMSG7PiOREju/chvAwacBoa8cSjDHGBMHBolKWb85m2SbnLmVLVkGN7ft1aMWI5PYMT05kWK92tIwM/V/foR+hMcY0AWXlFazdlet+7LVmZy7lFeq1fWKrKHc/yvDkRDq0jm7AaP3DEowxxgSAqrIt+xDL0jNZmp7F8s3ZHCwu89o+KiKMYb0TGOFKKAM6xiESeo+9fGEJxhhj/CT3UAlfbcpm2SYnqezaX/Os+OM7t3b3owzt0ZboFuENFGnDsARjjDF1VFJWwart+1m2KZNl6Vmsy8hDvT/1olObaNdIr/ac0SeBhFZRDRdsEFiCMcaYWlJV0vflu/pRMlmxNYdDJeVe28dGhnNanwRXP0p7+rSPbfSPvXxhCcYYY2qQebCYr1wz5pdtymTvAe+z5sMEBnWNZ6TrLiWlezwtwptv0XpLMMYYQ/XFt9rFRjKoaxv2HChm/S81F4vs3q6luwzLab0TadPS86z55sgSjDGmWauoUP7xxSae+DydMtew4eyCEhanZXps3zo6gtP7JDKiXyIj+rane0JozJoPRZZgjDHNzi95he75KF9tyiK7oMRr24gwYUj3tu4yLIO6xofkrPlQZAnGGNPkFRSX8c2WbHcF4k378mt97vf3n0OrKPtVWRf2XTPGNDnlFcoPGXks3ZjJ0k1ZrNmxn9Jy7+OHwwQ8TarvEh9jyaUe7DtnjGkSduYccq/i+PXmbPIKvReLjIwI4+SeTrHI4X0TSd9zkD/P/ZHC0sNDjv29+FZzZAnGGNMo5RU6xSKXpjtrzW/PPlRj+wEd4xjZz0koJ/dsR0zk4VnzJ3Rpg4RJQBffao4swRhjGoXS8gq+35nrvktZuzPX42OtSklxUQxPTmRkcnvO6JtI+7iaZ80HevGt5sgSjDEmJKkqW7IKnDVS0jP5ZksO+TUUi4xpEc6w3u3ctb2Sk1o1q1nzocgSjDEmZOQUlPCVa32UpemZ7M4r8tpWBE7s0obhfRMZkdyeIT3iiYpoWsUiGztLMMaYoCkuK2fVtv0s3eQklJ92H6ixWGSX+BjXrPn2nN4ngbaxkQ0XrPGZJRhjTINRVdL2HmRZehZfpmfx7dZsikq9rzUfFxXBqX0S3LW9eia0tMdejUjIJRgRiQGmAZcAyUAksBdYCTyhql9VafsqMLWGy6Wp6oDARWuMOZZ9B4rcywIv3ZRF5kHvxSLDw4STusUzvG8iI/slMrhrPBHNuFhkYxdSCUZEegH/BfoC+4AvgGKgJ/BrYC3wlYdTvwI2edj/S0ACNcZ4VVhSzoqt2a5+lCzS9h6ssX2vxFhXP0oip/ZJoHW0FYtsKkImwYhILPA50AeYAcxQ1dIqxxOABC+n/1NVXw14kMaYo1RUKD/tPsBS16JbK7ftp6Tc+2Ov+JYtOKOPU9dreN9EurWzYpFNVcgkGOCvOMnlNVW978iDqpoNZDd4VMaYo2TkFrrXmv9qUxb7D3mfNd8iXBjao617+PDxndtYschmwucE4+ojmQycBrQHvlDVp1zHegKdgZWq6r086dHXjASuc23O9DUmY4x/VV0bpXN8DNPG9CGhVbQ7qWzJKqjx/H4dWjG8b3tG9EtkWK92tIwMpb9lTUPx6V9dREYB7wBJgAAK5FVpMhz4F04Cmu3DpYfiPP7aqarrReR04HzXvj3AZ6q6vIbzR4vIIKAVzoCAZcDnqur9Pt0Y49HcNRncM+cHd12ujNxC7p7zY43nJLaKYnjfBIa7ant1bBPdEKGaEFfrBCMi/YFPgJbAa8CXwEtHNJuL0yk/Ed8SzImu13QvI8PuE5HZwFWqWujh/Ks97PtZRKao6g8+xGFMs/fI/LRqRR89iYoI45Re7dxzUgZ0jLPhw+YovtzB/AUnuVyhqu8AiEi1BKOq+SKyHkjxMY52rteRQDgwC/g/nD6XkcA/gIuAA8A1Vc77HlgFLAS2A62BIcCDwGBggYgMUdUMT28qItcD1wN0797dx5CNaZp253r6G87x+1G9GdG3Pak92xLdwmbNm5r5MsB8DLCuMrnUYCfQqY5xRAAvqep0Vd2sqrmq+hHOHZECU0Wkd+VJqvqEqj6tqj+raoGq/qKq84BTgG9wHuXd4+1NVfUFVU1V1dT27dv7GLIxTU9JWQUR4Z7vRLrEx3DPuQMZnpxoycXUii8Jpj2QVot2ZTh3Or6oOlD+xSMPqupKnDuVMODMY13MNcDgIdfmr3yMxZhma9Z/0zwuzGVro5i68CXB5AK1qWXdG8j0MY5tVT7f6qVN5f6OtbzmBter1d82phYWp+3jhS+3uLdbR0cgOHcuD0060UrZG5/50gezEhgjIj1VdZunBiIyGKfv4z0f41hd5fMEPCeoRNdrbRfTrpyUWfvFt41ppvbkFXH7u2vd22f2b8/LU08mzOarmHrw5Q7meSAKeFdEjuoRF5FOHB5V9rwvQbg64Ve4Nsd6uHZbnM57cBJdbUx2vX7nSyzGNDflFcot76whp8CZutahdRSPXjLYkoupt1onGFdn+8tAKrBJRL52HTpTRBYCm3GSwAuqurgOsTzoer1PRE6q3Cki0cBzQBucfpjlrv0nicj5IlKtt1FEIkTkT8DNrl2P1yEWY5qNpxams2JrDgBhAk9cmkJCq5pXfzSmNnyaaKmq14rIBpyRWae6dvd0feQDf1XV/61LIKr6sYjMAu4AVojICpxhyqfgVAfIAC5Tda8W0RP4AMgRkY3ALiAOZ05NZ6ACuEtV59clHmOag+Wbs3l6Ubp7e9qYZE7r463knzG+8bl+g6rOEpGncBJMb5x5KzuBZap6qD7BqOp0153RNJy5NC2BHcBjwExVrdo3sxZ4EicB9XC1V5xE8wrwrKquqk88xjRl2fnF3PLOGve69sN6tePmscnBDco0Kb7M5B8DlKnql65hwF+6PvxKVT/AuTM5VrutwK3+fn9jmoOKCuX299ayz7U2S7vYSJ6ckmJFKI1f+dLJvwA4qsqxMabxeXHpFpakHX4g8OjkwVY/zPidLwkmB6eQpDGmEVu9Yz+PzD88Z/r6kb0Z3T8piBGZpsqXBLMSsOWHjWnE8gpLufntNZS5Ol4Gd4vnjnNshr4JDF8SzKPASSIyJVDBGGMCR1W5e/Y6du13ilnGRUfwzGUpREbYmvcmMHwZRZaNU9X4TRG5CKcjfjvgsfSqqq72tN8YExxvrNjBpz/ucW8/fNEgW67YBJSvpWIUZ6GxSa4Pb9THaxtjAuin3XnM+ORn9/aVp3bnVyf6WvTcGN/4kgRW4yQOY0wjUlBcxrS31lBS5izwOqBjHH8977ggR2Wag1onGFVNDWQgxpjAuPfDH9mSVQA4ZfefuXyIrediGoT17hnThM1etYs5qw8v6Dpj4gn0TWoVxIhMc2IJxpgmanNmPvd++KN7e9KQLlw8tGsQIzLNjS+lYoYcu9VhNorMmOApKi3nj2+u5lBJOQC928cy49cnBDkq09zUZRRZbdgoMmOC6P/N+5kNe5yVyCMjwnjmsiHERtmPpGlY/hhFFoZTzbid6/haoLz+oRlj6uLTH37hjW92uLfvPW8gx3VuHcSITHPlt1FkIjIMeBGnXP7EesZljKmDnTmHuHP2Ovf2uSd05MpTewQxItOc+a2TX1VX4CSWsTiLhhljGlBpeQXT3l7DwaIyALq2jWHmRYMQsRL8Jjj8OopMVbcA3wG/8ed1jTHHNmt+Gt/vzAUgIkx4+rIU2sS0CHJUpjkLxDDlLJzljI0xDWRx2j6e/3KLe3v6uP6kdG8bxIiM8XOCEZGWwDAg35/XNcZ4tyeviNvfXevePrN/e64b0TuIERnj8GUeTLsaDrfCWSvmHqAz8G494zLG1EJ5hXLrv9eQU1ACQFJcFI9eMpgwW/rYhABfhilnHrsJgrPq5T11C8cY44unF6XzzZYcAMIEnpySQkKrqCBHZYzDlwSzH+8TLUuADGAh8Liq7qtvYMaYmi3fnM1TC9Pd29PGJHNan4QgRmRMdb7Mg0kMZCDGmNrLzi/mlnfW4Fr5mGG92nHz2OTgBmXMEazYpTGNTEWFcvt7a9l3sBiAdrGRPDklhXDrdzEhptYJRkQ+EpGba9HuJhH5qH5hGWO8+eeyLSxJO9wl+uglg+nYJjqIERnjmS99MOfjzHE5lsHAeXULxxhTkzU79vP3z9Lc29eP7M3oAUlBjMgY7wLxiCwSqAjAdY1p1vIKS5n29hrKXB0vg7vFc8c5/YMclTHeBSLBpADZAbiuMc2WqnL37HXs2l8IQFx0BM9clkJkhHWjmtBV4yMyD30pY2voX4kA+uOUifmg/qEZYyq9sWIHn/64x7398EWD6NauZRAjMubYjtUHc36VzxXo5vqoSRpwV32CMsYc9vPuA8z45Gf39hXDuvOrEzsFMSJjaudYCWaC61WAj4AFwJNe2pYAGar6s5fjxhgfFRSXcdPbqykpc7o1B3SM497zjwtyVMbUTo0JRlXnVX4uIquApVX3GWMC674Pf2JLZgEAMS3CeebyIUS3CA9yVMbUji8z+U8OZCDGmOpmr9rF7NW73NszJp5A36RWQYzIGN/YEBRjQtDmzHzu/fBH9/aklC5cPLRrECMyxne+TLQE3GX7zwP6Aa1x+meOpKp6Sz1jM6ZZKiot549vruZQSTkAvRNjmTHxhCBHZYzvfEowInId8DgQU3W361WrbCtgCcaYOnhw3no27DkIQGREGE9fnkJslM9/CxoTdL7UIhsN/J9r82ngO9fntwPPAb+4tp8FjlmzzBhztE9/+IXXv9nu3r73vIEc37lNECMypu586YO5zfU6XlVvBdYDqOrjqvpHnEdmbwNXAJ/4NUpjmoGdOYe4c/Y69/a5J3TkylN7BDEiY+rHlwRzCrBGVZd5Oqiqh4DfAcXA//ghNmOajdLyCqa9vYaDRWUAdImPYeZFgxCxEvym8fIlwbQFNlXZLgUQEXe9ClUtBpYBY/0SnTHNxKz5aXy/MxeAiDDh6ctTaBPTIshRGVM/viSYbKDqIPwc1+uR9/ARgK3bakwtLU7bx/NfbnFvTx/XnyHd2wYxImP8w5cEsx3oXmV7Hc6IsYsqd4hIW2AUsNMv0RnTxO09UMTt7651b5/Zvz3XjegdxIiM8R9fxj4uAu4SkS6qmgHMA/KAv4lIL2AXMAVoA7zi90iNaWLKK5Rb3llDTkEJAElxUTx6yWDCbOlj00T4kmD+DfTGGS2Woap5IvIH4F/Ab6u0W4918htzTE8vSuebLc6TZhF4YspJJLSKCnJUxviPL7XI1gGXHbHvHRH5DpgItAM2AO+6OvuNMV58syWbpxamu7enjUnm9D6JQYzIGP+rdYIRkTCcEjBadb+qbgYe9XdgxjRV2fnF3PLOGlwrHzOsVztuGZsc3KCMCQBfOvnLODx73xhTBxUVyu3vrWXvAecmv11sJE9OSSHc+l1ME+RLgsnHeQRmjKmjfy7bwpK0TPf2o5cMpmOb6CBGZEzg+JJg0gBbp9WYOlqzYz9//yzNvX3diF6MHpAUxIiMCSxfEszLwHARsbrhxvgor7CUaW+voczV8TK4WzzTxw0IclTGBFatE4yqPge8CSwUkT+KSOfAhWVM06Gq3DNnHbv2FwIQFx3BM5elEBlh6/2Zps2XUWQHXJ/GAk8BT4lIKVDiobmqqtUYNwZ4c8UO/vPDHvf2wxcNolu7ljWcYUzT4MtES0+LgUe6PowxHvy8+wAPfPKze/uKYd351YnWlWmaB18STFzAojCmCSooLuOmt1dTUlYBwICOcdx7/nFBjsqYhuPLTP6CQAZiTFNz34c/sSXT+bGJaRHOM5cPIbpFeJCjMqbhWC+jMQEwZ/UuZq/e5d6eMfEE+iZ5espsTNPlc4IRka4icr+IfCYiq0TkgSrHhojI5SJij9NMs7U5M5+/zv3RvT0ppQsXD+0axIiMCQ5f+mAQkSnAP4EYnLVgFFhbpUkP4HXgN65XY5qVotJy/vjmag6VlAPQOzGWGRNt6phpnmp9ByMiJwOv4SSWB3CWRT6ygNJ/cErK/NpfARrTmDw4bz0b9hwEIDIijKcvTyE2yqe/44xpMnz5n383TkI6T1WXAIhUzy+qWiwiG4Dj/RWgMY3Fpz/8wuvfbHdv//W8gRzf2aaDmebLlz6Y4cB3lcmlBruwmmWmmdmZc4g7Z69zb48/viNXndojiBEZE3y+JJh4YPsxWznXtMmXptkoLa9g2ttrOFhUBkCX+BgevnjQUXf4xjQ3viSYLKBnLdr1A36pUzTGNEKz5qfx/c5cACLChKcvT6FNTIsgR2VM8PmSYL4GhorIIG8NRGQUMBD4or6BGdMYLE7bx/NfbnFv3zGuP0O6tw1iRMaEDl8SzFNAODBHRE4/8qCIDMEp6V8BPOuf8IwJXXsPFHH7u4dH6Y/s157rR/QOYkTGhBZfyvUvBf4G9AaWishunHkw54vIZpzllHsB96vqqgDEakzIKK9Qbn3ne3IKnGLiSXFRPDZ5MGG29LExbj7N5FfVB4BLcJZO7ogzDyYRJ7FsBa5U1QfrE5CIxIjInSLynYjkisghEdkqIu+JyBlezrlcRJaKSJ6I5IvISteaNVYKxwTEM4s2sXxLNgAi8MSUk0hsFRXkqIwJLT7PAFPV2cBsEemOczcTDuxU1Y31DUZEegH/BfoC+3D6copxBhf8GqdqwFdHnPMscCNQBCwESnEmgT4DjBWRS1S1vL6xGVPpmy3ZPLnw8H/3aWOSOb1PYhAjMiY01XmKsaruAHb4KxARiQU+B/oAM4AZqlpa5XgCkHDEORfhJJc9wEhVTXft7wAsBi4EbgKe9FecpnnLzi/mlnfW4Fr5mFN6tePmMX2DG5QxIarOj5BEpLWIDBSRASLij+nKf8VJLq+p6n1VkwuAqmZ7uEu6x/V6V2VycbXdC/zBtXm3PSoz/lBRodzx3lr2HigGoG3LFjw1JYWIcPvvZYwndammfKWIrAJygB+Bn4BsV2Xlq+sShIhEAte5NmfW8pyuwFCcJZvfO/K4qn4BZOD0FZ1al7iMqeqlZVtZnJbp3n508mA6tokOYkTGhLZaPyITZ1rya8DlHC5ymed6bQOkAK+IyDiczn71IY6hOI+/dqrqetcw6PNd+/YAn6nq8iPOSXG9/qSqhV6u+x3QxdX2ax/iMcZt7poMHpy3nsz8Yve+60b0YsyADkGMypjQ58sdzI3AFTgz+m8G2qpqW1Vti1NGZhpOx/wU4I8+xnGi6zVdRF7F6ci/B7geuA/4WkTeF5GYKuf0cr3WVL6mso+oVw1tjPFq7poM7nx/XbXkIgL9O9iSR8Yciy8J5nqcEV1nquozqlp594KqHlDVZ4ExOKO4rvcxjnau15HA1cAsnJFkbXFGj2UAF1F9Amfl8oA1LeWc73q13wbGZz/vPsBds9dRUl5Rbb8qPL4g3ctZxphKvowi6wcsUdX13hq4Hm8tBkb5GEdloosA/qmq06sc+8g1qfNbYKqI/D9V3cLhx3S+PIqrRkSux5UMu3fvXtfLmCZmZ84hHvt8I3O/z8Dbg97dud6eyhpjKvmSYA4A+2vRLhc46GMcVdu/eORBVV3pGliQCpwJbKlyTk0LnVce8xiPqr4AvACQmppa50Rlmob9BSU8u3gTry3fftRdy5E6x8fUeNwY41uCWQKcLiIRqlrmqYGIRACn4Xuxy21VPt/qpc1WnATT8Yhzalp0o5uH6xtTTWFJOa98vZXnlmx2l9yvdFynODZnFlBcdjjhxLQIZ/q4/g0dpjGNji8J5q84o7JeFJGbVbXaXYGItMIpiBnvauuL1VU+TwAyPbSpnCpd2a+yxvV6vIjEeBlJdvIRbY1xKyuvYPbqXTz+eTp7DhRVO3ZSt3juOXcAw3onMHdNBo/MT2N3biGd42OYPq4/E1O6BClqYxoPXxLMBOB94Brg1yLyCYfvNnq6jrfBqah8vofllB/zdmFVzRCRFcAwnDIvG6oeF5G2wBDX5krXOTtFZLVr/yU4Q6irnjMK6IozzPnIIc6mGVNVFqzfx98/20D6vvxqx3onxnLn+P6MO76je8GwiSldLKEYUwdS2+kqIlKB06FeNXNUnuxpH1WOqaqGH+P6E4CPcIY6j1PV7137o4FXgUuBVcDJlXNsRORinEmWe4ARqrrJtT8Jp1TMccCtqnrMUjGpqam6cuXKYzUzjdyq7fuZ+el6vttWvTuxfVwUt56VzOTUbrSwmfnG1JqIrFLVVE/HfLmDeYx6jNg6FlX9WERmAXcAK1x3NNnAKUBnnKHKl1WdwKmq74vIczhlYX4QkQUcLnbZGpiLU/TSNHObM/P5+2cbmP/T3mr7W0VF8PuRvfndiF60jKxzaT5jjAe1/olS1TsCGYjrPaaLyNc4kzZTgJY4kyUfA2aq6lF9M6p6o4gsw5ncOQqnuvMGnEd1z6lqzcOBTJO270ARjy9I592VOymvOPz3UYtw4YphPZg2pi8JVmbfmIAIuT/ZVPUD4AMfz3kLeCswEZnG6GBRKc9/sYWXlm2lsLT6ag0XDO7M7ef0o0dCbJCiM6Z5CLkEY0x9FJeV8+Y3O3hm8Sb3apOVzuibwN3jB3JiV38U/zbGHIvPCUZEBgOjcfpFvJWSVVW9pT6BGeOLigrl43W7mfXfNHbmVB+xflyn1tx97gBGJCdy5OhGY0zg+FJNORp4A2cRL6g+cuxICliCMQ1iWXoWMz9bz48ZB6rt79o2hjvO6c8FgzsTFmaJxZiG5ssdzP8Ck3AmOr4LpHN40qMxDe7HjDwe/mwDS9Ozqu1v27IFN41J5spTuxMVUePoeGNMAPmSYC7FWf9lqKvYpDFBsTPnELP+m8aH3++utj+6RRi/G96L34/qQ+voFkGKzhhTyZcE0xb43JKLCZacghKeXpTOG99sp7T88JDjMIHJqd249ax+tsKkMSHElwSzlQBOtDTGm0MlZby8bCvPf7GFg8XVi1GefVwH7hrfn75JtuSPMaHGlwTzGnCPiCSoanagAjKmUll5Be+t2sXjn29k38HiaseG9mjLPecOILVnOy9nG2OCzZcEMwtnpvwCEblBVVcEKCbTzKkq//15L3//bAObM6svWNqnfSx3jR/A2cd1sCHHxoQ4X0rFlIvIFGAZ8LWIHAR2AZ5KsaiqDvZTjKYZWbkth4c+3cCq7dWLUXZoHcVtZ/Xj4qFdibBilMY0Cr7Mg+mCs5BYL5w5MK1xqhV7Yn01xifpew/y8GdpLFhfvRhlXFQEN5zZh2vO6EVMpA05NqYx8eUR2SNAb+Bb4GlgEzYPxtTTnrwiHv98I++t2kmVWpS0CBeuOrUnN43pS7vYyOAFaIypM18SzFnATmC0l9Ujjam1vMJSnv9iMy9/tZWi0upPWSee1Jnbz+lPt3YtgxSdMcYffEkw0cBiSy6mPorLynl9+XaeWbyJ3EOl1Y6NSE7krvEDOKGLFaM0pinwJcH8ACQEKhDTtFVUKB+uzWDW/I1k5Fb/G+WELq25e/xAhicnBik6Y0wg+JJgHgfeFpEhqro6UAGZpkVV+TI9i5mfbmD9L9WLUXZr5xSjnDDIilEa0xT5kmAWAY8CC0VkJjAf78OUUdWc+odnGrMfduXx0Kfr+Xpz9Xm57WIjmTamL1cM60FkhA05Nqap8iXBVC5XLDiVlf+3hrbq47VNE7I9u4BZ/93Ix2urF6OMaRHOtSN6cf3I3sRZMUpjmjxfksB+bH6LqUFWfjHPLNrEmyuqF6MMDxMuPbkbt45NJqm1FaM0prnwZSa/9cAajwqKy3hp2Vae/2IzBSXl1Y6NP74j08f3p0/7VkGKzhgTLPYYy/hk7poMHpmfxu7cQjrFR3N6n0SWpGWSlV+9GOXJPdty97kDGdqjbZAiNcYEmyUYU2tz12Rwz5wfKCx17lJ25xbx/qpd1dokJ7XirvEDGDswyYpRGtPM+ZxgRCQV+CNwGtAe+Leq3ug6Ntq1/0VVzfR+FdMYPTI/zZ1cjtSxdTR/OrsfFw3tSrgNOTbG4GOCEZHpOKPHKqsOKs4M/0otgBlANvC8PwI0oWN3rvciDkumn0l0CytGaYw5rNaTEETkHOBhIAv4LdAHZ8hyVQtwRptN8Pg4rf0AABrmSURBVFeAJjTMXrXL6xDCLvExllyMMUfx5Q7mT0AJMF5V1wJHPWNX1QoR2Qj081uEJqjKyiuY+ekG/rlsq8fjMS3CmT6ufwNHZYxpDHxJMCcDKyqTSw12ASfWPSQTKnIPlTDt7TUsTc9y7+sQFwUC+w4U0zk+hunj+jMxpUsQozTGhCpfEkwssKcW7Vpy9KMz08hs3HuQ615byfbsQ+59Zx/XgccvPYlWUTb40BhzbL78pthD7R59DcBZN8Y0UvN/2sOf/v19tUmTN49N5taxyVaU0hhTa75UGvwCGCQiI701EJGJOKteLqxvYKbhVVQoTy5I5/evr3Inl5aR4Tx3xRD+dHY/Sy7GGJ/4kmBmAeXAByJyuYi4lxsUkQgRmQS8CBQDT/o3TBNoBcVl3Pjmah5fsNG9r2vbGGb/4XTOPbFTECMzxjRWvtQi+0FEbsCZ3/I6UIYzD2YKcBkQiVO6/1pVTQ9ArCZAdmQf4rrXVpK296B73+l9Enjm8iG0i40MYmTGmMbMp8U4VPVlYDjwGU4yEZyJluHAEmCMqv7LzzGaAPpqUxYXPLusWnL57Rk9ee2aUyy5GGPqxefhQKq6AjhPRCKBzjjJ5RdVPVTzmSaUqCqvfLWNB/+znvIKZwplZHgY/+/CE5ic2i3I0RljmgKvCUZEXgaWue5ajqKqJcC2AMVlAqiotJy/zv2xWqHKpLgo/u+qoQzpHlrVj8vKysjJySEvL4+ysrJgh2NMkxUREUGbNm1o164dERH+mYpQ01V+43r1mGBM47T3QBG/f30V3+/Mde8b3C2eF64aSocQWwysoqKCnTt3EhUVRffu3YmMjLQKzcYEgKpSUlJCdnY2O3fupEePHoSF1X85c1sQvRlZs2M/E55eVi25XDSkK/++/tSQSy4A+/fvJyIigk6dOhEVFWXJxZgAERGioqLo1KkTERER7N+/3y/XtSnZzcR7K3fylw9+pKS8AnCWMf7zrwZyzRk9Q/YXd35+Pu3atQvZ+IxpakSE+Ph49u/fT0JCQr2vZwmmiSsrr+DB/6znla+2ufe1iWnBs5cPYXhyaK+CXVRURMuWLY/d0BjjNy1btmT37t1+uZYlmCZsf0EJN729mq82Zbv39evQihevTqVHQmwQI6udiooKvzwHNsbUXlhYGBUVFX651rESzMUicmYdrquq2qcO5xk/2bDnANe9tpKdOYcXCRt3fAcendy4ilXa4zFjGpY/f+aO9ZumlevDV97WpjIN4LMff+FP767lUJVilbed1Y9pY/paPTFjTIM5VoL5DGcVS9MIVFQoTyxM56mFhyv1xEaG89ilJzHu+I5BjMwY0xwdK8HsUdUvGiQSUy/5xWXc9u/v+fznve593du15MWrU+nfMS6IkRljmivrQW0CtmcXMOkfX1VLLsP7JvLRTWdYcmmCRMTnj9/85jfBDts0Q42nt9d4tDQ9k5veWkNeYal73++G9+KecwcQEW5/PzRFU6dOPWrfnj17mD9/PrGxsVx88cVHHR8+fHhDhNYsffLJJ0yYMIHzzjuPTz75JNjhhBRLMI2UqvLSsq3873/W46pVSWREGA9deCIXDe0a3OBMQL366qtH7VuyZAnz588nMTHR43FjgsESTCNUVFrOnz/4gTmrM9z7OrSO4vmrUjmpW3wQIzPGmMO8PkNR1TBVvaYhgzHHtieviEufX14tuaR0j+fjm4ZbcvGDuWsyOGPmInrdPY8zZi5i7pqMY5/UCK1bt46pU6fSo0cPoqKiaNeuHePGjWP+/Pke2ycmJiIiZGVl8f7773P66acTFxdH+/btueyyy9i1y6nMXVZWxsMPP8xxxx1HTEwMnTt35rbbbqOwsPCoa95xxx2ICLNmzWLjxo1ceumlJCUlER0dzaBBg3jmmWdqnPC3dOlSLrnkEjp37kxkZCRJSUlMmjSJb7/99qi2+fn5iAitWrVCVXnuuedITU0lLi4OEXFX6l67di1/+ctfOPXUU+nUqRORkZF06NCBCRMmsGjRoqOum5qayoQJEwCYN29etX6v888/3+P3z5PU1FREhJUrV3rd//nnnzNu3DgSEhIQERYsWFCt7UcffcR5551HUlISkZGRdOnShauuuooNGzZ4/R4Gmj2kb0RWbd/PhGeWsXZXnnvf5NSuvHP9qSSFYLHKxmbumgzumfMDGbmFKJCRW8g9c35ocknm5ZdfZujQobz22mvEx8dzwQUXcPzxx7N48WLGjx/Pww97n5nwyCOPMGXKFKKjoxk/fjzR0dG88847jBw5kry8PC644AIefPBB+vbty1lnnUVBQQFPPPEEV1xxhddrrl+/npNPPpnly5czZswYRo0aRVpaGtOmTeOqq67yeM4DDzzAyJEjmTNnDl27dmXixIn06tWLuXPncsYZZ/DWW295fb9rrrmGm2++mbi4OCZMmMBJJ53knlz40EMP8dBDD5Gfn09KSgoTJ06kS5cufPLJJ5x11lk8//zz1a41YcIExo4dC0DXrl2ZOnWq++Oss87yGoOvXn31VcaNG8e+ffs455xzGD16tLukvqpy7bXX8utf/5oFCxaQnJzMxIkTSUxM5I033iA1NZXFixf7LRZfiKrNiQRITU3VI/96CCXvfreTv86tXqzyvvOP4+rTejTZ2e7r169n4MCBXo/3vHteA0bjX9tmnufX6y1ZsoTRo0fTo0cPtm3b5rXdihUrOOOMM4iLi2POnDmMHj3afWzNmjWce+65ZGZmsnz5ck455RT3scTERLKzs2ndujULFizg5JNPBqCgoIAxY8bw7bffcsIJJxAWFsbnn39OUlISAJs2bSIlJYX8/HxWr15NSkqK+5p33HEHjz76KABXXXUV//znP4mMdFZR/emnnxg9ejSZmZn861//4uqrr3afN3v2bC6++GJ69uzJnDlzql1z0aJF7juHDRs20L17d8C5g4mLi3N/LQsWLGDw4MFHfX8WLlxIv3796Nat+qJ7X375Jeeeey6qyo4dO0hMPFzHrzad/JXfv8zMzGrnVkpNTWXVqlV89913pKamHrUf4I033vCYqGfNmsX06dMZMmQI7777Ln36HC6i8tZbb3HllVfSvn17tmzZQmxs7UpEHetnryoRWaWqqZ6O2R1MiCstr+D+D3/kztnr3MmlbcsWvH7NKUw9PXQrIZvQ9MADD1BeXs5TTz1VLbkApKSkMHPmTCoqKnj22Wc9nn/XXXe5kwtAbGws06ZNA+DHH3/kueeecycXgL59+zJ58mQAr39Ft2nThqefftqdXACOP/547r33XgAef/zxau3vu+8+AF577bVqyQVgzJgx3HnnnRQWFvLSSy95fL97773XY3IBGDt27FHJBWDkyJFce+21FBYWMm9ew/9hc+GFF3pMLsXFxTz00EOEh4fz/vvvV0suAJdffjlXX301+/bt4913322ocN2skz+E5RSUcOObq/hmS45734COcbx4dSrd2lmVYeOb4uJiFi5cSEREBBdeeKHHNqNGjQJg+fLlHo+PHz/+qH19+/YFIC4ujtNPP/2o48nJyQBeK/Sef/75tGnT5qj9V111FTfffDNr164lNzeX+Ph4tm/fzs8//0yHDh0YMWJEnb6GSZMmedxfKTc3l3nz5rFu3TpycnIoLXWmAKxfvx6AjRs31nh+IHiL+ZtvviEnJ4dhw4bRq1cvj21GjRrFv/71L5YvX85vf/vbQIZ5FEswIern3Qe4/vWV7Np/uHP03BM6MuuSwcQ2omKVgeTvx0yVfTCFpYdruMW0COehSScyMaWLX98rGHbv3k1xcTGA+3GRN5mZmR73d+169BD4Vq1aeT1W9XhRUZHH495+McbHx9OmTRvy8vLIyMggPj6eLVu2ALB3795j3r17+hpatGhB586dvZ7zzjvvcMMNN5CXl+e1zYEDB2p830Do0aOHx/2V348VK1bU6fsRaPabKgTNW/cLd7y3ttovutvP7sdNY/raI7EAqkwij8xPY3duIZ3jY5g+rn+TSC4A5eXO/6eoqCimTJlSY9voaM+DRmpaPiGQSytU/r+v/BoSEhKqjdLyxFPCi4yM9Brnpk2buPrqqykrK+Nvf/sbF198MT169CA2NhYR4bHHHuP2228nEP3WxyqPHxMT43F/5fejR48enHnmmTVe46STTqpTbPVhCSaEVFQojy/YyNOLNrn3tYqK4PFLT+Ls4zoEMbLmY2JKlyaTUI5UuRyuqvLCCy9U6/MIJm+DEvLy8sjLy0NE6NSpE4C7fyQuLs7vE0rnzp1LaWkpU6dO5f777z/q+KZNmzycVTuV3+v8/HyPnfzbt2+v03Urvx+9evUKyQm21skfIg4WlXL96yurJZeeCS354MbTLbkYv4iNjWXEiBGUlJTw0UcfBTsct48//tjjY6c33ngDgEGDBtG2bVsA+vfvT+/evdm2bdtRc0bqKyfH6ev01MlfUFDAhx9+6PG8yuRROZfGky5dnD9aPM1JWbFihfu9fTV8+HBatWrF8uXLycgIveH0lmBCwNasAi78x9csWL/PvW9EciIf/nE4yR2sWKXxn7/97W+Eh4dzww038MEHHxx1vLy8nPnz53ucVBgoeXl53HLLLe7OdHA61GfMmAHALbfcUq39Aw88AMDkyZNZsmTJUdcrKipizpw57uG9tTVgwADA6Yep+gu/qKiI3//+914HKVQmj7S0NK+Puirnyjz00EMUFBS492/ZsoVrr73Wpzirio2N5e6776a4uJgLLriA77///qg2Bw8e5PXXX3f31zQke0QWZF9szGTaW6s5UHT4r5/rR/bmznH9rVil8buRI0fy4osvcsMNNzBp0iR69+7NgAEDaNmyJbt27WLjxo3k5OQwY8YMxowZ0yAxXXPNNbz33nssWrSI0047jdzcXBYvXkxJSQmXXnrpUSOfrrjiCrZs2cL999/P6NGjGThwIP369SMsLIxdu3axYcMG9y/VoUOH1jqOyZMnM3PmTNavX0+fPn0YNWoUERERLFu2jKKiIm688Ub+8Y9/HHXecccdR//+/UlLS2Pw4MGkpKQQGRnJ4MGD3UO4b7vtNl599VW+/PJL+vfvz7Bhw8jKyuK7775jzJgxREREeEwOtfHnP/+ZHTt28MILLzBkyBAGDx5M7969KSsrY+fOnaxfv56ioiKWLl1K79696/QedWW/wYJEVXnhy8389pVv3cklMiKMxy8dzJ9/NdCSiwmY3/72t6xbt44//OEPREREsHjxYv7zn/+QmZnJsGHDePbZZ7nuuusaLJ6BAwfy7bffcvLJJ7Nw4UKWLFlCcnIyTzzxBG+++abHc+69915WrFjB1KlTKSws5LPPPmPBggXk5eVx1lln8corr7hLuNRWdHQ0X3/9NbfeeiuJiYl89tlnLF++nLPPPps1a9Z4nXgoInz88cdceOGF7N27lzfffJOXXnqpWtmdDh068PXXXzN58mSKioqYN28ee/fu5d5772Xu3LmEh4f7FOuR7//888+zYMECLr74YjIzM/nkk0/48ssvKSoqYtKkSfz73/+uNoGzodhMfpeGnMlfVFrO3bPXMff7w7fcHVtH88LVQxnU1eqJVfJlNrFpfCpn8j/yyCPccccdwQ7HVOGvmfz2iKyB7c4t5Pevr+KHjMPj7If2aMtzVw4hKc7qiRljmg5LMA1o5bYcbnhjFVn5Je59U07uxv/8+niiIup+i2yMMaHIEkwDefvbHdz34Y+UljuPJCPChPsnHMeVpzbdYpXGmObNEkyAlZZX8MDHP/P6N4cnUrWLjeTZy4dwWp+EIEZmTHDNmjWLWbNmBTsME0CWYAIoO7+YP7y5mm+3Hh5TP7BTa168eihd21qxSmNM0xYyY2FF5FUR0Ro+jpoCW5dzGspPu/O44JmvqiWX8wZ1YvYfTrPkYoxpFkLxDuYrwFPRn1/8fE7AfLx2N9PfX0tRqTOrVwTuOKc/N57Zx/pbfKSq9j0zpgH5c+pKKCaYf6rqqw1wjt+VVyiP/jeNfyzZ7N7XKiqCJ6ecxNiBVk/MV2FhYVRUVNRrEpoxxjcVFRV+q4wdigmmUTpQVMqt73zPog2H64n1SozlxauH0jfJ6onVRXR0NIcOHTrm2iXGGP85dOiQ1+UBfGUJph7mrslwrx0SHiaUVRy+tRzVrz1PXZZCm5gWQYywcWvVqhW5ubm0atXKHpMZ0wBUldzcXGJjY/1yvVBMMKNFZBDQCtgLLAM+V9WaVuSpyzn1cuTqh1WTyw2j+jB9XH/Cw+yXYn20bduWAwcO8Msvv5CQkEBkZKQlGmMCQFUpKSkhOzubsrIy9/II9RWKCeZqD/t+FpEpqvqDH8+pl0fmp1VbcbJS25YtuPvcAYF4y2YnLCyMbt26kZOTw44dO2pcb8MYUz8RERG0adOGpKSkJtkH8z2wClgIbAdaA0OAB4HBwAIRGaKqGfU8x01ErgeuB+jevbtPwe7OLfS4P/dQqcf9pm4iIiJISkoiKSkp2KEYY3wUMvNgVPUJVX1aVX9W1QJV/UVV5wGnAN8AScA99T3niPNfUNVUVU1t3769T/F2jvfcCeZtvzHGNDchk2C8UdUS4CHX5q8CdY6vpo/rT0yL6sNnY1qEM31c/0C8nTHGNDqh9IisJpUz8rsE+Jxam5jiXLZyFFnn+Bimj+vv3m+MMc1dY0kwlVUh8wN8jk8mpnSxhGKMMV6E/CMyl8mu1+8CfI4xxhg/CYkEIyInicj5IhJ+xP4IEfkTcLNr1+P1OccYY0zDCZVHZD2BD4AcEdkI7ALigBOBzkAFcJeqzq/nOcYYYxpIqCSYtcCTOMOLewApgOIkjVeAZ1V1lR/OMcYY00BCIsGo6lbg1kCfY4wxpuGIP2v/N2YikolTDcDUTiKQFewgTMDZv3PzUJ9/5x6q6nGmuiUYUycislJVU4Mdhwks+3duHgL17xwSo8iMMcY0PZZgjDHGBIQlGFNXLwQ7ANMg7N+5eQjIv7P1wRhjjAkIu4MxxhgTEJZgjDHGBIQlGFNvItJfRG4TkU9FZJOIFIlInogsF5FbRSQy2DEa34jI5SKy1PXvmC8iK0XkjyJivzMaMRFpISJjReRREflGRH4RkRIRyRCR90XkTL++n/XBmPoSkV046+4UAStxyvV0AE4DooE1wFmqmhO0IE2ticizwI04/54LgVJgLE6tvw+AS1S1PHgRmroSkbOAz12be3CWnC8AjgNOcO2foar3+eX9LMGY+hKRhcCbwLuqml9lf0/gE+B44DVVnRqUAE2tichFwPs4v3xGqmq6a38HYDEwELhVVZ8MXpSmrkRkDM4fD0+q6tIjjl2K83McDoxR1cX1fj9LMCaQRGQ4sBTnr+E2ruWsTYgSkZXAUGCqqr52xLFRwBKc5NNFVSsaPkITSCLyT+B3wMuq+rv6Xs+ep5pAW+N6jebwKqMmBIlIV5zkUgK8d+RxVf0CyAA6Aqc2bHSmgVT+vHb1x8UswZhAS3a9lgDWBxPaUlyvP6lqoZc23x3R1jQtlT+vv/jjYpZgTKDd7Xr9RFWLgxqJOZZerteaqorvOKKtaSJEpCPwG9fmbH9c0xKMCRgR+Q1wKXAI+HNwozG10Mr1WlBDm8pBHHEBjsU0IBGJAN4A2gALVfVjf1w3JBYcM8EjIn8HLqjDqWNVNaOG644FnsdZZfT3qppWxxBNwxHXq438aX7+D2co+k7gSn9d1BKM6Qz0r8N5LbwdcI0c+xCIBG5W1TfqGJtpWAddr61qaFN57GANbUwjIiJP4owc24Pzh+Mef13bHpE1c6p6papKHT62ebqeiJwO/AeIBe5S1acb8usx9bLN9dqjhjbdjmhrGjEReRS4GcjESS7p/ry+JRjjNyJyKvApzvP5v6rq34MckvFN5RDV40Ukxkubk49oaxop1+PxPwHZwNmq+rO/38MSjPELETkFmA+0Bv6mqg8GOSTjI1XdCazGebR5yZHHXRMtu+I8SlnesNEZfxKRmcB0YD9OclkbiPexBGPqTUSGAv/FSS4zVPV/ghySqbuHXK8Pi0jfyp0ikgT8w7U502bxN14iMgO4C8jFSS4Buxu1UjGm3kQkB2iL8x/2wxqa3qGqWQ0TlakrEfkH8Aec8j4LOFzssjUwF7jYil02TiJyAYd/RlcCP3lpukFVZ9b7/SzBmPoSkdr+J+rlbXCACS0icjnwR+BEnOKHG4CXgefs7qXxcs1Ne6UWTb9Q1TPr/X6WYIwxxgSC9cEYY4wJCEswxhhjAsISjDHGmICwBGOMMSYgLMEYY4wJCEswxhhjAsISjDHGmICwBGNMHYiI1uHjVde5Z7q2lwT3qzAmsGw9GGPq5l8e9nUExuGsCPm+h+PLAhqRMSHGZvIb4yciciawGNiuqj1raNcS6A4cUtUd3toZ09jZHYwxDUxVD+HU9jKmSbM+GGMamLc+GBHp6dq/TUTCRORPIvKTiBSKyC4Recx194OItBWRJ1xti0UkXUT+VMN7iohMEZH/ikiW65wdIvKiiPQM6Bdsmi1LMMaEpreAB4CtOGvtxAK3AbNFpB2wArgU+A6nb6cn8KiI/PnIC4lIC5w+obeB4cDPwEc4fUXXAqtFJDXAX49phuwRmTGhpwfOWiz9VHU3gIh0w1mmeDzwBbAWuEpVi1zHzwM+Ae4WkSdcj+EqzQAmAV8CV6jqrsoDInIT8DTwjogMUNWygH91ptmwOxhjQtPNlckF3MsZv+Ha7AH8oTK5uI7PA9YBcYD7bsR1t3MzkA9cUjW5uM57BpgH9AHODcyXYporSzDGhJ5SYJGH/Ztcryu9rAya7nrtXGXfaCAGZwGpfV7e7wvX62m+BmpMTewRmTGhZ4+XR1X5rtddHo5VPR5dZV9v1+t5tVh5tH0t4zOmVizBGBN6jrUksS9LFoe7XtOAb47RdoUP1zXmmCzBGNO07XS9/qCqvwlmIKb5sT4YY5q2BTh9OmeJSHywgzHNiyUYY5owVd0LPAvEAx+JyIAj27gmbV4rIh0aPEDTpNkjMmOavjtxRpZNBn4Uke9xJnBGA92AgUCk63VvsII0TY/dwRjTxKlqqapeCvwaZzJmZ9fnp+P8kfkWcCGwOWhBmibJqikbY4wJCLuDMcYYExCWYIwxxgSEJRhjjDEBYQnGGGNMQFiCMcYYExCWYIwxxgSEJRhjjDEBYQnGGGNMQFiCMcYYExD/H4wam+3c4XNsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "#where 0 is 11 am\n",
    "\n",
    "time = np.array([-3,-2,-1,0,1,2])\n",
    "T_a = np.array([55,58,60,65,66,67])\n",
    "\n",
    "print(\"part a)\")\n",
    "\n",
    "#function to return ambient temperature at a given time (11am = t=0)\n",
    "\n",
    "def Temp_a(x):\n",
    "    for i in range(len(time)):\n",
    "        if (time[i]==x):\n",
    "            return T_a[i]\n",
    "        \n",
    "print(\"Using the function, the ambient temperature at time -1  is %0.2f\" %Temp_a(-2))\n",
    "\n",
    "plt.rcParams.update({'font.size': 22})\n",
    "plt.rcParams['lines.linewidth'] = 3\n",
    "plt.plot(time,T_a,'o-',label='Temperature')\n",
    "plt.legend()\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Temperature')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The figure above looks correct based on the  given data. A better way to obtain T_a(t) would to be given more data at intermediate steps between each hour. In this figure, the temperature slope at each time is assued to be linear between each hour, however, in reality, it is not the case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 348,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "part b)\n",
      "Incorporating the new ambient temperatures as a function of time, the temperature from 11 am to 1 pm in half hour intervals is:\n",
      "[85.         82.49491982 80.37588923 78.5933141  77.10397147]\n"
     ]
    }
   ],
   "source": [
    "#### import numpy as np\n",
    "print(\"part b)\")\n",
    "\n",
    "\n",
    "#assuming at each half hour interval the temperature is half of the difference  of the temperatures\n",
    "ambient_temperature = np.array([65,65.5,66,66.5,67])\n",
    "#starting at time t=0 (11am)\n",
    "time=np.array([0,0.5,1,1.5,2])\n",
    "Temp_nonlinear=np.linspace(85,85,5)\n",
    "for i in range(1,len(time)):\n",
    "    K=0.275\n",
    "    Temp_nonlinear[i]=ambient_temperature[i]+(Temp_nonlinear[i-1]-ambient_temperature[i])*np.exp(-K*0.5)\n",
    "print(\"Incorporating the new ambient temperatures as a function of time, the temperature from 11 am to 1 pm in half hour intervals is:\")\n",
    "print(Temp_nonlinear)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 342,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After two hours, the temperature with changing ambient temperature is 77.10\n",
      "After two hours, the temperature without changing ambient temperature is 76.54\n"
     ]
    }
   ],
   "source": [
    "print(\"After two hours, the temperature with changing ambient temperature is %0.2f\" %Temp_nonlinear[4])\n",
    "\n",
    "#using previous code to find temperature after 2 hours with constant ambient temperature\n",
    "\n",
    "t=np.linspace(0,2,101)\n",
    "T_num=np.linspace(85,85,101)\n",
    "for i in range(1,len(t)):\n",
    "    K=0.275\n",
    "    T_num[i]=65+(T_num[i-1]-65)*np.exp(-K*0.02)\n",
    "print(\"After two hours, the temperature without changing ambient temperature is %0.2f\" %T_num[100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Part b) (continued)\n",
    "It is logical that with the new ambient temperature increasing, that the corpse will decrease in temperature slower since the difference in temperature between the corpse and the ambient temperature becomes smaller as time increases. This is displayed above by finding the two temperatures at time = 2 hours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 326,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 85.          88.31653838  92.49044445  97.4269931  103.23859914\n",
      " 110.12794839 118.25390203]\n",
      "from this data, after 1.5 hours, the corpse is 97.4 degrees\n",
      "Due to a lack of data for the abmbient temperature between each hour, the time of death is best estimated to be a few minutes before 9:30am\n"
     ]
    }
   ],
   "source": [
    "#Because the ambient temperature is changing, a while loop can not be used to solve this problem \n",
    "#To solve this problem, newtons law of cooling is reversed and solved for a law of heating\n",
    "#t=0 is set to 11 am and the order of the ambient temperature is reversed\n",
    "#by working backwards, the time it takes to reach 98.6 will be the time of death\n",
    "\n",
    "#the ambient temperature is decomposed into half hour intervals to obtaina more accurate model (assuming linear slopes)\n",
    "ambient_temperature = np.array([65,62.5,60,59,58,56.5,55])\n",
    "#starting at time t=0 (11am)\n",
    "time=np.linspace(0,3,7)\n",
    "Temp_nonlinear=np.linspace(85,85,7)\n",
    "for i in range(1,len(time)):\n",
    "    K=0.275\n",
    "    Temp_nonlinear[i]=ambient_temperature[i]+(Temp_nonlinear[i-1]-ambient_temperature[i])*np.exp(K*0.5)\n",
    "print(Temp_nonlinear)\n",
    "\n",
    "print(\"from this data, after 1.5 hours, the corpse is %0.1f degrees\" %Temp_nonlinear[3])\n",
    "print('Due to a lack of data for the abmbient temperature between each hour, the time of death is best estimated to be a few minutes before 9:30am')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}