Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

compmech-project01/01_Getting-started-project.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
278 lines (278 sloc) 29.1 KB
Raw Blame
Open in GitHub Desktop
{
"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": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"K= -0.275 1/Hr\n"
]
}
],
"source": [
"dT = 11 #85F-74F\n",
"dt = 2 #HOURS\n",
"Ta = 65\n",
"T = 85\n",
"K = -(dT/dt)/(T-Ta)\n",
"print('K=', K, '1/Hr')"
]
},
{
"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": 44,
"metadata": {},
"outputs": [],
"source": [
"def measure_K(Temp_t1,Temp_t2,Temp_ambient,delta_t):\n",
" ''' Determine the value of K based upon temperature of corpse \n",
" when discovered, Temp_t1\n",
" after time, delta_t, Temp_t2\n",
" with ambient temperature, Temp_ambient\n",
" Arguments\n",
" ---------\n",
" your inputs...\n",
" Temp_t1 - initial body temp\n",
" Temp_t2 - final body temp\n",
" Temp_ambient - room temp\n",
" delta_t - change in time\n",
" \n",
" Returns\n",
" -------\n",
" your outputs...\n",
" K - empirical constant\n",
" '''\n",
" \n",
" dT = Temp_t1 - Temp_t2\n",
" dTa = Temp_t1 - Temp_ambient\n",
" dT_dt = dt/delta_t\n",
" K = -(dT/dt)/(dTa)\n",
" return K"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"K= -0.275 1/Hr\n"
]
}
],
"source": [
"T1 = 85\n",
"T2 = 74\n",
"Ta = 65\n",
"dt = 2\n",
"K = measure_K(85,74,65,2)\n",
"print('K=', K, '1/Hr')"
]
},
{
"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": 49,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3hUZfbA8e9JLwRSCC2A9F4ChICiotJZwYqCvSL7U1fdXVddd0Xdta91rYiFtYBYEDsIKggoVUqQ3kMJAUIC6eX8/rhDCDAJE8hkUs7nee4zc+/cciZoTu773vc9oqoYY4wxx/PzdQDGGGOqJksQxhhj3LIEYYwxxi1LEMYYY9yyBGGMMcatAF8HUJHq16+vLVq08HUYxhhTbSxdunSfqsa6+6xGJYgWLVqwZMkSX4dhjDHVhohsK+0za2IyxhjjliUIY4wxblmCMMYY41aN6oMwxlSM/Px8kpOTycnJ8XUopoKEhITQtGlTAgMDPT7GEoQx5gTJyclERETQokULRMTX4ZjTpKrs37+f5ORkWrZs6fFxXm1iEpF7RGS1iCSJyGQRCRGRh0Vkp4gsdy3DSzl2qIisE5GNInK/N+M0xhwrJyeHmJgYSw41hIgQExNT7jtCryUIEYkD/gQkqGoXwB8Y7fr4eVWNdy3fuDnWH3gFGAZ0AsaISCdvxWqMOZElh5rlVP49vd1JHQCEikgAEAbs8vC4RGCjqm5W1TxgCnCRNwLMyTrML+8/TNK8L7xxemOMqba8liBUdSfwH2A7sBtIV9WZro/vEJGVIvK2iES5OTwO2FFiPdm17QQiMlZElojIktTU1HLHGRAYRNuN75A3/7VyH2uM8Y79+/cTHx9PfHw8jRo1Ii4urng9Ly+vUmJYvHgxIsLs2bNP6zz/+Mc/eOGFF8rc57nnnjum+WfIkCEcOnSo3NeaOHEid999d7mPK403m5iicP7qbwk0AcJF5BrgNaA1EI+TOJ51d7ibbW4rG6nqBFVNUNWE2Fi3o8XLFBAYxJbGw+matZA9e5LLfbwxpuLFxMSwfPlyli9fzrhx47jnnnuK14OCgiolhsmTJ3P22WczefJkr1/r+AQxY8YMIiIivH7dk/FmE9NAYIuqpqpqPvAZcJaqpqhqoaoWAW/iNCcdLxloVmK9KZ43T5Vbs/NvJlAKWTvzHW9dwhhTQSZNmkRiYiLx8fH83//9H0VFRRQUFBAZGcm9995Lz549GTJkCAsXLqR///60atWKb75xujonTpzIJZdcwpAhQ2jfvj3//ve/3V6jqKiITz/9lEmTJvHtt98W37Vs3LiRLl26cPPNN9O5c2eGDRtW/Iv99ddfp3fv3nTv3p1Ro0aRnZ19zDnXrVtHYuLRX3dr1qwhMTGR559/nr1793LOOecwcOBAAJo2bcrBgwcBeOedd+jWrRvdu3fnxhtvBGD69On06dOHHj16MHjwYPbu3VuBP+GjvPmY63agr4iEAdnAAGCJiDRW1d2ufS4BktwcuxhoKyItgZ04ndtXeSvQRu0S2BrYmoZbPqew6B/4+1nnnDFHPPLlan7flVGh5+zUpC7jR3Qu93FJSUlMmzaNBQsWEBAQwNixY5kyZQpXXHEF6enpDB48mGeeeYYRI0bw8MMPM3v2bFasWMFtt93G8OHOA5OLFi0iKSmJoKAgevfuzYUXXkh8fPwx15k7dy4dOnSgVatW9OvXj++++46RI0cCzi/6yZMn07VrVy699FI+//xzRo8ezahRoxg3bhwA999/P++++y5//OMfi8/Zvn17QkJCSEpKokuXLrzzzjvceOON/PGPf+TZZ5/l559/JjIy8pg4VqxYwVNPPcWCBQuIjo7mwIEDAJx77rmMHDkSEeH111/n2Wef5amnnir3z/NkvNkHsRD4BFgGrHJdawLwtIisEpGVwPnAPQAi0kREvnEdWwDcAcwA1gBTVXW1t2IFyOowio66kSVLFnjzMsaY0zBr1iwWL15MQkIC8fHxzJkzh02bNgEQGhrKoEGDAOjatSvnnXceAQEBdO3ala1btxafY8iQIURFRREeHs7FF1/MvHnzTrjO5MmTGT3aeehy9OjRxzQztWnThq5duwLQq1ev4nOvXLmSc845h65duzJlyhRWrz7xV9bNN9/MO++8Q0FBAR9//DFjxowp8/v+8MMPXHnllURHRwMUv27fvp3BgwfTtWtXnnvuObfXqgheHSinquOB8cdtvraUfXcBw0usfwOc8Aist7QZcCMFq54mbcH/ILFfZV3WmCrvVP7S9xZV5aabbuJf//rXMdsLCgqO6Zvw8/MjODi4+H1BQUHxZ8c/7nn8en5+PtOmTeObb77hkUceoaioiIMHD5KZmQlQfF4Af3//4nNfd911fPvtt3Tp0oWJEyfy66+/nhD/qFGjePzxx+nXrx9nnnnmCXcM7r6vu8dTb7/9dv7+978zfPhwZs2axZNPPlnmeU6VzcXkEhTZiK31ziQ+bSZ7D2b6OhxjjBsDBw5k6tSp7Nu3D3Cedtq+fXu5zjFz5kwOHjxIVlYW06dPp1+/fid83rt3b3bs2MHWrVvZvn07I0aM4Isvyn4UPjMzk0aNGpGfn8+HH37odp+wsDAuuOAC7rjjjuL+BICIiAi3Ty0NHDiQKVOmFDctHXlNT08nLi4OVWXSpEnl+v7lYQmihIg+19BIDvDL7Gm+DsUY40bXrl0ZP348AwcOpFu3bgwePJiUlJRynePss8/mqquuokePHowZM+aE/ofJkydzySWXHLPtsssuK/WX/hGPPvooiYmJDBo0iE6dSh/Xe/XVVxMYGMiAAQOKt40dO5aBAwcWd1If0a1bN/72t79x7rnnEh8fz7333gvAww8/zCWXXEL//v1p2LChR9/7VIiq26dHq6WEhAQ9rYJB+TlkPt6KeX4JDHpwOn7WWW1qqTVr1tCxY0dfh1HhJk6cSFJS0knHJXjTk08+SW5uLuPHH9/67n3u/l1FZKmqJrjb3ybrKykwhNQzhnPOli9YtG4bfTu28HVExpgaZMSIEezYsYMffvjB16F4xBLEcRr3v4ngrR+z4acP6dvx774OxxhTgW655RafXv/LL7/06fXLy/ogjhPc4kz2Bzel7e4v2X8419fhGGOMz1iCOJ4I2m00ff1+Z+b8xb6OxhhjfMYShBv1z3KGamQt/ZCa1IlvjDHlYQnCnagWpEYncH7ObBZt3u/raIwxxicsQZSi3pnX0cpvD7/8PMPXoRhT6/h6uu+mTZty5ZVXFq9PmTKl0ju4d+zYcUwM5XHNNdfw+eefn3YMliBKEdT1EvIkmAabP+NgVuXMP2+McVSF6b4XLlzIunXrKuVaxysoKKBZs2Z89NFHPrn+EZYgShNSl6zWw/iDLGD60i2+jsYY41IZ030D/OUvf+Hxxx8/YfvxBYA6dOhAcnJy8VTgN910E507d+a6665jxowZnHXWWbRr144jg3gPHz7MDTfcQGJiIj169Ch+9HXixImMHj2aCy+8kGHDhrFx48biUd4FBQXcc889dOnShW7duvHqq68CMH78eHr37k2XLl0YN25chfeZ2jiIMkT2vRY2fs72Xz5Dz77favSa2unb+2HPqoo9Z6OuMKz8E8xV1nTfAGPGjOHll19myxbP/0Bct24dU6dOpUOHDvTs2ZPg4GAWLFjAp59+ypNPPsknn3zCo48+ytChQ3n33XdJS0ujT58+xbPQ/vLLLyxfvpyoqCg2btxYfN7XXnuNXbt2sWLFCvz9/YvnZLrrrrt45JFHUFWuuuoqvvvuO4YNG1bun2tp7A6iLK3OJys4lr6HZrJs+0FfR2NMrVdZ030DBAQE8Je//KVcM6W2adOGTp064efnR6dOnYrnVioZw8yZM3nssceIj4/n/PPPJycnp3jCwcGDBxMVdWIV5lmzZjFu3Dj8/f2Bo9N+z549m8TERLp3786cOXMqfNpvu4Moi58/AfGjOe/XV3lswQp6nXGeryMypvKdwl/63lIZ032XdMMNN/D000/Trl274m0BAQEUFRUVr5csFVpyKvDSYlBVPv/8c1q3bn3MtebOnUt4eHip3/v4OLOysrjjjjtYtmwZcXFx/OMf/zgmlopgdxAnEdTzKgKlkKA1n5GRk+/rcIyp1Spjuu+SgoKC+NOf/sSLL75YvK1FixYsXboUcJqrduzYUa7rDxkyhJdeeql4/bfffjvpMYMHD+a1116jsLAQcKb9zs7Oxs/Pj/r163Po0CE+/fTTcsXhCa8mCBG5R0RWi0iSiEwWkRAReUZE1orIShGZJiJuK2aIyFZX5bnlInIaU7SepoadyIrpwgjmMH2518piG2M8UBnTfR/v1ltvPebR2lGjRpGSkkKPHj146623aNWqVbmuP378eLKysujatSudO3fm4YcfPukxt912G40aNSquTT116lRiYmK4/vrr6dKlC5dccgl9+vQpVxye8Np03yISB8wDOqlqtohMxakQtwv4QVULROQpAFW9z83xW4EEVd3n6TVPe7rvUugvryIzHuD/6r7Mq392WxDPmBrFpvuumco73be3m5gCgFARCQDCgF2qOtNVcxrgV6Cpl2M4bdJ1FEUSQLcD37EqOd3X4RhjTKXwWie1qu4Ukf8A24FsYKaqzjxut5uA0kaCKDBTRBR4Q1UnuNtJRMYCYwGaN29eIbGfoE4sha0HcOmG+by4cDNdm/bwznWMMV7l6+m+qxuv3UGISBRwEdASaAKEi8g1JT5/ECgAPijlFP1UtScwDLhdRM51t5OqTlDVBFVNiI2NrdDvUFJgz6tpIGmkrphhndWmVrCJKmuWU/n39GYT00Bgi6qmqmo+8BlwFoCIXA9cCFytpUStqrtcr3uBaUCiF2M9uXZDKQiux3Cdw0eLyvfUgjHVTUhICPv377ckUUOoKvv37yckJKRcx3lzHMR2oK+IhOE0MQ0AlojIUOA+oL+qZrk7UETCAT9VPeR6Pxh41IuxnlxAMAFdL2fY0vcZMW8VN/RrQaC/PSVsaqamTZuSnJxMamqqr0MxFSQkJISmTcvX5evNPoiFIvIJsAynKek3YAKwGggGvncN/PhVVceJSBNgoqoOBxoC01yfBwAfqup33orVY71uIHjJW5ybOYMvV/Ti0p5Vvn/dmFMSGBhIy5YtfR2G8TGvPebqC956zLUkfXsYe3Zs5Oa6E/j67vNsfiZjTLXmy8dcaxzpO47Gupe41Ln8vMHjIRrGGFPtWIIor/Z/QOvGMTb4eybM3ezraIwxxmssQZSXfwCSeCu9dRWpm5aRtNMGzhljaiZLEKei5/VoQAi3BH3Pmz/bXYQxpmayBHEqwqKRbldwsf885q1cT3Ka26d1jTGmWrMEcaoSbyOwKJcr/H7knflbfR2NMcZUOEsQp6pRF2hxDreG/MDHi7aQnm3TbxhjahZLEKejz21EF6RwZsEiPli4zdfRGGNMhbIEcTraD4d6zbm7zg+8M38ruQWFvo7IGGMqjCWI0+HnD4m30jFvJfUPr2f6b1ZxzhhTc1iCOF09r0UDw7gr4kcm/LyZoqKaM3WJMaZ2swRxukKjkG5XMqhgDvv37uKn9Xt9HZExxlQISxAVoc9t+BflMjb8Z96YYwPnjDE1gyWIitCgI7Tsz7UBs1i6ZS8rdhz0dUTGGHPaLEFUlD7jqJObwkXBv9kkfsaYGsESREVpNwQiz+Duuj/wbdJutu+36TeMMdWbVxOEiNwjIqtFJElEJotIiIhEi8j3IrLB9RpVyrFDRWSdiGwUkfu9GWeF8POHxLE0O7SCrn5beWue3UUYY6q3MhOEiDQWkbtF5FMR+UVEfhCRl0RkiJyklJqIxAF/AhJUtQvgD4wG7gdmq2pbYLZr/fhj/YFXgGFAJ2CMiHQ6lS9YqXpcA4Hh/CP2Z6YuSSYtM8/XERljzCkrNUGIyJvA+659XgRuBP4MzAMuBuaLyNknOX8AECoiAUAYsAu4CJjk+nyS61zHSwQ2qupmVc0DpriOq9pCIyF+DL0yZhOWf4B3F2z1dUTGGHPKyrqDeFlVB6jqc6o6V1XXqupyVZ2qqn8ELgBKfehfVXcC/wG2A7uBdFWdCTRU1d2ufXYDDdwcHgfsKLGe7Np2AhEZKyJLRGRJampqGV+nkiSOxa8oj382Wcxb87bYXYQxptoqK0GUWXBZVXNUdX1pn7v6Fi4CWgJNgHARucbDuNw1X7kdoqyqE1Q1QVUTYmNjPTy9F8W2h9YXcGHu1+Tm5fDanE2+jsgYY05JWQniyyNvRGTqKZx7ILBFVVNVNR/4DDgLSBGRxq7zNsb9XUgy0KzEelOc5qnqoc84AjJTeKjVBiYt2Mqe9BxfR2SMMeVWVoIo+Vd821M493agr4iEuTq0BwBrgC+A6137XA9Md3PsYqCtiLQUkSCczu0vTiEG32gzCGLacmXOJ6gW8tIPG3wdkTHGlFtZCUJLee8RVV0IfAIsA1a5rjUBeBIYJCIbgEGudUSkiYh84zq2ALgDmIGTVKaq6uryxuAzfn7Q/28E7V/Do+22MHXxDrbuy/R1VMYYUy6i6v53v4gUAuk4dxIRQMaRjwBV1ehKibAcEhISdMmSJb4Ow1FUCK/0oUD86ZryEIM6NeGlMT18HZUxxhxDRJaqaoK7z8q6gwgCYoH6QLDr/ZH1KtAbXMX5+UP/+wjYt5YnOmzlixW7+H1XxsmPM8aYKqKsBBGsqoWlLQAiElZJcVZPXS6FmLaMSHuPeiF+PDtzna8jMsYYj5WVIL4SkadE5CwRCTmyUUSai8j1rv6CC70fYjXmuovw37eGpzptZfbavSzZesDXURljjEfKShADgPnAXcBGEUkXkQM4Hc8tgFtV9VQef61dulwK9dsxeO+7xIYH8vSMdZTW72OMMVVJqQlCHV+o6pWq2lRV66lqtKomquojrpHS5mRcdxF++9byTOetLNpygLkbyhyDaIwxVYJHs7mKyGgR+bvrfTMR6eXdsGqYzpdA/facu/stmkUG88yMtVa72hhT5Z00QYjIy8D5wLWuTZnA694Mqsbx84f+f8Mv1bmLSNqZwbdJe3wdlTHGlMmTO4izVPU2IAdAVQ/gPAJrysN1F9Fn+5u0iw3j2e/XUVBY5OuojDGmVJ4kiHwR8cM1mlpEYgD7zVZefv5w3n1I6lqe7rSFzamZfLbMunGMMVWXJwniFeBTIFZEHsGpB/GUV6OqqTpdDLEd6L75DXrE1eGFWevJyS/0dVTGGOPWSROEqv4P+AdObYc0YJSqTvF2YDWS64kmSV3LEx22sCs9hw8Xbvd1VMYY49bJSo76i8gKVV2tqi+q6guqmlRZwdVIrruIDutepV+rSF75cSOHcwt8HZUxxpygzAThmlLjd1d9aVMR/Pyg/32wbx3/breR/Zl5vD1vi6+jMsaYE3jSB1EfWCMiM0TksyOLtwOr0TpdDLEdaZn0CoM71ufNuZtJPZTr66iMMeYYAR7s86TXo6ht/PzgvPvg4xv4V48NnL0+mie+WcNzV8b7OjJjjCl20gShqrMrI5Bap+NF0KATDX97kdvOeZeXf9rKlb2b0adVjK8jM8YYwLOR1IdEJMO1ZIlIroictLCBiLQXkeUllgwRuVtEPiqxbauILC/l+K0issq1XxWpAlSBivsi1vOnhquJiwzln9OTyLfBc8aYKsKTx1wjVLWuqtYF6gBXAy96cNw6VY1X1XigF5AFTHNN/ndk+6dAWf0Z57v2dVvtqNrrOBIadCZo3tM8/Id2rE85zLvzt/o6KmOMATycrO8IVS1S1U9wakmXxwBgk6puO7JBRAS4AphcznPVHH5+cMGDsH8DA7O+5oIODXhh1np2p2f7OjJjjPGoiWlkieViEfk3Tl3q8hjNiYngHCBFVTeUcowCM0VkqYiMLSO+sSKyRESWpKamljOsKqD9cGh1HvLj4zw6sBEFRcq/v1rj66iMMcajO4hRJZaLgHzXq0dEJAgYCXx83EdjKPvuoZ+q9gSGAbeLyLnudlLVCaqaoKoJsbHVsFS2CAx9CnIP0XT5c9x+fhu+XrWbueurYbIzxtQoHs3FpKrXupYbVfURnIpynhoGLFPVlCMbRCQAuBT4qLSDVHWX63UvMA1ILMc1q5cGHSBxLCx5h9vaZdIiJozxX6wmt8DmaTLG+I4nCeJVN9teKcc13N0pDATWqmqyuwNEJFxEIo68BwYDNXuKj/Puh7Bogr9/gEdGdmbLvkzenLvZ11EZY2qxUsdBiEgicCbOLK5/KvFRXSDQk5OLSBhOh/Ztx310Qp+EiDQBJqrqcKAhMM3pxyYA+FBVv/PkmtVWaCQMeAi+vIv+vecyvGtL/vvDRi6Kj6NZdJivozPG1EJl3UGE40yzEQDElljycPojTkpVs1Q1RlXTj9t+g6q+fty2Xa7kgKpuVtXurqWzqj7m+VeqxnpcC427w8x/8tCQM/D3Ex75crWvozLG1FKl3kGo6o/AjyLyjqpaW0dl8POHYU/D20NotPJ17howhie+Xcus31MY2Kmhr6MzxtQynvRBZIjIEyLyhYjMPLJ4PbLaqnlf6HoFzH+JmzoLbRvU4eEvV5OdZx3WxpjK5UmCeB/YCrTDqSS3B3A7PYapIIMeAb8AAmf9k0cv6kJyWjav/rTR11EZY2oZTxJErKq+AeS5Ju67npr8yGlVULcJnPtXWPsVZ7KCS3rE8caczWxOPezryIwxtYgnCSLf9bpHRIYAXYBm3gvJAHDm7RDVEr69nweGtCY4wI/xX6xGVX0dmTGmlvAkQTwuIvWAv+LUpn4fuNerURkICIahT8C+dTRY8x5/GdyOnzfs46uVu30dmTGmljhpTWqghaqmq+pKVT3H9eipVZSrDO2GQpuB8NMTXNM1jG5N6/HQ9CT2ZuT4OjJjTC3gSU3qSyspFnM8ERjyBORnEfDTv3n+yniy8wu595OV1tRkjPE6T5qY5onIiyJypoh0O7J4PTLjiG0HfcbBsvdonbeeB4d3ZM76VN77ddvJjzXGmNPgSU3q/q7XniW2KeB2dlXjBf3vg5VT4dv7uOamGcxeu5fHvl7DWa3r06ZBHV9HZ4ypoTypKHeOm8WSQ2UKqQsDx0PyImTlFJ6+rBthQf7c89Fy8gqsRKkxxjs8KRgUKyJviMhXrvVOInKD1yMzx+p+FTTrA9/dTwMO8MSl3Vi1M52XZpdWb8kYY06PJ30Q7wJzODr2YQPwF28FZErh5wcXvwaF+fDFnQzt3JBRvZry6k8bWbrtgK+jM8bUQJ4kiAaq+iFQBKCq+YBNDOQLMa1h0KOwcRYsfZfxIzsTFxXKPR+t4HBuga+jM8bUMJ4kiEwRicbpmEZEegOHvBqVKV3CzdCyP8x4kDqZO3juiniS07J41KYFN8ZUME8SxF+BL4FWIjIHp9DPnV6NypTOzw8uesWZGnz67fRuHsm4/q2ZuiSZGav3+Do6Y0wN4slTTEuA83Eed70L6KSqJ53NVUTai8jyEkuGiNwtIg+LyM4S24eXcvxQEVknIhtF5P7yfrEaLbIZDH0Sts2Hha9x98B2dImrywOfrWLvIRtlbYypGJ48xRQMjAMeBB4Axrq2lUlV16lqvKrGA72ALGCa6+Pnj3ymqt+4uaY/Tt3rYUAnYIyIdPL0S9UK8VdBu2Ew6xGC0jbwwpXxZOYW8DcbZW2MqSCeNDFNwvkF/yYwEWfA3KRyXmcAsElVPR3+mwhsdJUezQOmABeV85o1mwiMeBGCwmHabbSJCeWBYR34aV0q79soa2NMBfAkQXRS1etV9XvXchPQsZzXGY3Td3HEHSKyUkTeFpEoN/vHATtKrCe7tp1ARMaKyBIRWZKamlrOsKq5iIZw4XOw6zeY9xzXndmCc9rW57Fv1rDJakcYY06TJwliuevJJQBEpBfwi6cXEJEgYCTwsWvTa0BrIB7YDTzr7jA329y2m6jqBFVNUNWE2NhYT8OqOTpfAl0uhzlP4Zeykv+M6k5IoD93T1lOboE9jWyMOXWeJIiewK+uzuKNwCLgLBH5TUSWeXD8MGCZqqYAqGqKqhaqahFOs5W76nTJHFuUqCmwy4Nr1U7Dn4GwGJg2joZhwpOuUdYPfW4Fhowxp86TyfpOt+1/DCWal0SksaoeqXpzCZDk5pjFQFsRaQnsxGmiuuo046i5wqJh5Mvw4Sj48XGGDnqE289vzSs/bqJL03pc2/cMX0dojKmGPHnMdROQCgQD4UcWVd3k+qxUIhIGDAJKFhh6WkRWichKnMdn73Ht20REvnFdswC4A5gBrAGmqqqNBCtLu8HQ8zpY8BJsX8ifB7Xngg4NeOSL1SzcvN/X0RljqiE5WROEiIwHxgJbONoPoFVxRteEhARdsmSJr8PwnZwMeK0f+AfAuHmkFwZxySvzSc/O58s7z6ZJZKivIzTGVDEislRVE9x95kkfxFVAK1U926b7ruJC6sLFr8KBzfD9eOqFBjLhul7kFhRx23tLycm3TmtjjOc8SRCrgQhvB2IqSMtzoM8fYfGbsPpz2jSI4IUr41m1M50HPltlndbGGI950kn9GPCbq88g98hGVbVa1VXVoEcgeTF8/n8Q256BnTry50HteO779XRuUpdbzmnl6wiNMdWAJwliEvA8sArXlN+migsIhivfgzf6w5Sr4dYfuOP8Nqzelc7j36yhQ6O6nN22vq+jNMZUcZ40MR1Q1edco6hnH1m8Hpk5PXWbwBWT4OA2mHYbfijPXhFPmwZ1uGPyMnYcyPJ1hMaYKs6TBLFYRP4lIr1FpNuRxeuRmdN3xlkw5AlY/x3MeYo6wQFMuDaBoiLl1v8tISvPigwZY0rnSYJIBM4DnsOZYfUV4GUvxmQqUuKtTj3rOU/Cum9pUT+c/17Vk/Uph7j3Y5v51RhTOk8Gyp3jZrHHXKsLEWdCv8bd4bOxsG8j/dvFct/QDny9ajevzSlzrKMxphbzpB5ErIi8ISJfudY7icgNXo/MVJzAULjyffAPhClXQe4hxp7bipHdm/DMjHXM+j3F1xEaY6ogT5qY3gXmcHTyvA3AX7wVkPGSyOZw+TuwfwN8/kcEeOqybnSNq8ftHy5j0ZYDvo7QGFPFeJIgGqjqh7gecVXVfHXgAxsAACAASURBVMCG5FZHrfrDoEdhzZcw7zlCg/x554bexEWFcvO7i0name7rCI0xVYgnCSJTRKJxzcPkqg1xyKtRGe858w7ochnM/hdsnEVMnWDev7kPESEBXP/2IjZboSFjjIsnCeKvwJdAKxGZgzN1951ejcp4jwiM/C807Ayf3AwHttAkMpT3bukDwLVvLWJ3eraPgzTGVAWlJggR6QugqktwpuXuD9yFU4J0eeWEZ7wiKNwZaY3CR9dAXiatY+sw6aZEMrLzuWbiQg5k5vk6SmOMj5V1B/HqkTeqmqeqK1R1uarab46aILoVXPY2pKyGT2+BwgK6xNVj4vUJJKdlc8M7iziUk+/rKI0xPuRJE5OpqdoOhGFPw7pv4Is7oaiIPq1iePXqnqzelcHY/9kU4cbUZmVN1tdKRL4o7UNVHVnWiUWkPfBRyfMBDwFxwAggD9gE3KiqB90cvxWnM7wQKCitoIU5TX3GQvYB+OkJCI2CIY8xoGNDnh3Vnbs/Ws6dk3/jtat7EuBvf0sYU9uUlSBSgWdP9cSqug6IBxARf5za0tOA9sADqlogIk8BDwD3lXKa81V136nGYDzU/z7IToNfX4GwKDj3Xi7uEUd6dj7jv1jNfZ+u4pnLu+HnJ76O1BhTicpKEIdUdU4FXWcAsElVtwHbSmz/Fbi8gq5hTpWIM6lfdhr88G/nTqL3LVx/VgvSs/N57vv11A0N4KELOyFiScKY2qKsBLG1Aq8zGufx2OPdxLHNUCUpMFNEFHhDVSe420lExuLUzKZ58+YVEGot5ecHF73i1LX++q8QEgldL+fOC9pwMCuft+dvoV5oIHcPbOfrSI0xlUS8PZuniAQBu4DOqppSYvuDQAJwqboJQkSaqOouEWkAfA/cqapzy7pWQkKCLlmypGK/QG2Tnw3vXwY7FsKYj6DtQIqKlHs/Wcmny5K5/fzW/HVwe7uTMKaGEJGlpfXxVkbP4zBg2XHJ4XrgQuBqd8kBQFV3uV734vRdJFZCrCYwFMZMhgadnDES2xfi5yc8fXk3Rvduxis/buKh6aspKrJpwo2p6SojQYyhRPOSiAzF6ZQeqapuy5qJSLiIRBx5DwwGkiohVgMQUg+u+cypSvfhKNiThL+f8MSlXbnt3Fa89+s2/jx1OfmFVoHWmJrMk+m+PxWRP4hIuZOJiIQBg4DPSmx+GYgAvheR5SLyumvfJiLyjWufhsA8EVkBLAK+VtXvynt9cxrqxMJ1n0NgOLx/KRzYjIhw/7AO3DukPZ8v38W492ychDE12Un7IERkIHAj0Bf4GHhXVddWQmzlZn0QXpC6Dt4eCiF14aYZENEIgPd+3cZD05NIbBHNxOsTiAgJ9HGgxphTcVp9EKo6S1WvBnriPNn0vYgsEJEbRcR+K9R0se3hmk8gcx+8dwkccrqSru17Bi9cGc/SbWlc9abN3WRMTeRRs5GIxAA3ALcAvwEv4iSM770Wmak64no5Hddp2+DtIZC2FYCL4uOYcF0v1qcc4oo3frFZYI2pYTzpg/gM+BkIA0ao6khV/UhV7wTqeDtAU0W0PBeu/8IZTPfWEEj5HYALOjRk0k2J7EnP4fLXfmHLvkwfB2qMqSie3EG8rKqdVPUJVd1d8gObH6mWaZoAN33njLx+ZxjsWAxA31YxTL61L1l5BYx6/Rd+35Xh40CNMRWh1E5qEbm0rANV9bOyPvcF66SuJGlbXf0Re2D0B9D6AgA27j3ENRMXkZVXwMTre5PYMtq3cRpjTupUO6lHuJabgbeAq13LROCaig7SVCNRLZwnmqJbwwdXwOppALRpEMEnfzyT+nWCuXrir3ywcFvZ5zHGVGmlJghVvVFVb8SZE6mTql6mqpcBnSstOlN11WkAN3zlNDt9fCMsfReAplFhTLu9H2e1rs+D05J4cNoq8gpsQJ0x1ZEnfRAtjut7SAFsxjYDoZHOiOu2g+DLu+Dn50CVeqGBvH1Db8b1b80HC7dz9cRfST2U6+tojTHl5EmC+ElEZojIDa45lL4GfvRyXKa6CAqD0R9C11Ew+xH4/p+gir+fM+r6xdHxrNqZzsiX57EqOd3X0RpjysGTgXJ3AK8D3XEKAE1wPeJqjMM/EC6ZAL1vhQX/hS/ugMICwBkr8cm4s/AT4fLXFzB9+U4fB2uM8ZSn8ystAH4AZgPzvReOqbb8/GD4M051ut/ehw+vgKwDAHSJq8f0O/rRvVkkd01ZzuPfrKHQZoM1psrzZKDcFTgT5l0OXAEsFBGrAmdOJALn/x1GvAhb5sKb58MeZxLe+nWC+eCWPlzb9wwmzN3Mje8uJj0r38cBG2PK4slkfSuAQa66DIhILDBLVbtXQnzlYuMgqpAdi2DqdZCTDiP/C12P/k0xedF2HpqeRFxkKG9el0DbhhE+DNSY2u10Cwb5HUkOLvs9PM7UZs0SYewcaNwdPr0ZZjxY3C8xJrE5k2/ty+HcQi56ZT4fLd6OtysbGmPKz5Nf9N+VeIrpBpynmL45yTHGQERDuO4Lp/P6l5fh/UucWWGBhBbRfHXn2cQ3i+S+T1dx23tLbUZYY6oYj2pSu6bdOBsQYK6qTvN2YKfCmpiqsN8+gK/ucQbYXfkeNOkBQFGR8vb8LTz93TrqhQXy9OXdOL99Ax8Ha0ztcdo1qVX1M1X9M/AY8LmHF23vqhh3ZMkQkbtFJFpEvheRDa7XqFKOHyoi60Rko4jc78k1TRXW42pnoj9VZzbY5R8C4Ocn3HJOK6bf0Y/osCBufGcxD01PIjvPKtUZ42ulJggR6SsiP4nIZyLSQ0SScOpCp7jqSpdJVdeparyqxgO9gCxgGnA/MFtV2+I8NnvCL38R8QdeAYYBnYAxItLpFL6fqUriesJtc5z+ic//CF//FQqcZqWOjesy/Y5+3Hx2S/73yzYu/O/PJO20gXXG+FJZdxAvA48Dk3HGQNyiqo2Ac4EnynmdAcAmVd0GXARMcm2fBFzsZv9EYKOqblbVPGCK6zhT3YXXh2s/hzPvgMVvwqQRxQWIQgL9+eeFnXj/5j5k5hZy8SvzefWnjTZmwhgfKStBBKjqTFX9GNijqr8CnGI96tE4iQag4ZG5nVyv7hqc44AdJdaTXdtOICJjRWSJiCxJTU09hdBMpfMPgCGPwWVvQcpqeK2fM9mfqz/s7Lb1+e7ucxjSuRFPf7eOMRN+ZceBLN/GbEwtVFaCKDkF5/G1JD3+k05EgoCRwMfliEvcbHN7TVWdoKoJqpoQGxtbjksYn+t6OfzfAqfp6cu74IPLIWMXAJFhQbx8VQ+eu6I7v+/OYPiLPzN50XaK7G7CmEpTVoLo7upYPgR0c70/st61HNcYBixT1RTXeoqINAZwve51c0wy0KzEelNgVzmuaaqLyOZw7XQY9gxsnQ+v9oWVU0EVEeHSnk359q5z6BxXlwc+W8Vlry+wvgljKklZ9SD8VbWuqkaoaoDr/ZH1wHJcYwxHm5cAvgCud72/Hpju5pjFQFsRaem6AxntOs7URH5+0GcsjJsH9dvBZ7c6o7BdYyaaRYcx+da+PH9ld3YcyGLky/N4+IvVZOTYVB3GeJNH4yBO+eQiYTh9Ca1UNd21LQaYCjQHtgOjVPWAiDQBJqrqcNd+w4EXAH/gbVV97GTXs3EQNUBRISx4CX58HELqOfM6dfhD8cfp2fk8O3Md7/26jfp1gvnHHzoysnsTRNy1ShpjTqascRBeTRCVzRJEDZKyGqbdBntWQfcxMPRJp0CRy8rkg/zj8yRWJqdzVusYHr2oC20a1PFhwMZUT5YgTPVUkAdzn4Gfn4WIRvCHZ6HdUGfWWKCwSJm8aDtPf7eW7PxCbj2nFXde0JbQIH8fB25M9XHaI6mN8YmAILjgQbj5ewiqA5NHO0867dsAgL+fcE3fM/jhr+cxsnscr/60iYHPzWHG6j02+Z8xFcDuIEz1UJgPiybAT09Cfhb0GecUJwqpW7zLws37+ef0JNanHCbhjCj+NrQDiS2jfRi0MVWfNTGZmuPwXpj9qFO1LjwWBo6H7lc5T0IB+YVFTF2yg5dmbyAlI5fz2sdy75D2dG5Sz8eBG1M1WYIwNc/OpfDtfZC8GOJ6wbCnoenR/8az8wqZ9MtWXvtpE+nZ+Yzo3oQ/D2pHy/rhvovZmCrIEoSpmYqKYNVU+P4hOJwC8VfDgPFOHQqX9Ox8JszdxNvztpJXWMQVCc24a0BbGtUL8WHgxlQdliBMzZZ7COb+B355BQJC4Ny/QuJYCAor3mXvoRxe+WEjHy7ajp8IN5zVgnH9WxMVHuTDwI3xPUsQpnbYvwm+ewA2zIDwBnD23dDrxmMSxY4DWTz//XqmLd9JnaAAbujXguvPakH9OsE+DNwY37EEYWqXbQucp522zDmaKBJugsDQ4l3W7sng+e/XM/P3FAL9/bi8V1NuPaeV9VGYWscShKmdts6HOU/ClrlQpyH0uxsSbjwmUWxKPczEnzfz6bKd5BcWMbhTQ27r35qezd0WOjSmxrEEYWq34xPF2fdArxuOSRSph3KZtGAr//tlKxk5BfRuEcXYc1szoEMD/PxsnidTc1mCMAZg6zyn6Wnrz0cTRc/rj+mjyMwt4KPFO3hr3hZ2HsymdWw4Y89txcU94ggOsCk8TM1jCcKYkkomitAoJ0n0vgUij5YgyS8s4ptVu3ljzmZ+351B/TpBjEpoxujezTgjxvopTM1hCcIYd7b9AgtfgzVfOusdLnSm8DjjrOIJAVWVeRv3MWnBNn5Ym0KRQr82MYxJbM7gTo0ICrDpzEz1ZgnCmLIc3AGLJ8KySZCdBg27Qp/bnJKoJfop9qTn8PGSHUxZvIOdB7OJCQ/i8l5NubJ3M1rF2lTjpnqyBGGMJ/KyYNXHsPAN2LsaQqOdzuzeN0O9psW7FRYpP29IZfKi7cxas5fCIqVvq2jGJDZnaJdG1ldhqhWfJQgRiQQmAl0ABW4C7gbau3aJBA6qarybY7cCh4BCoKC0L1CSJQhTIVSdfoqFr8O6bwBxqtr1uAZaXwD+Ryvu7s3I4eOlyUxZvJ0dB7KJCgtkRPcmjOjehF7No+wJKFPl+TJBTAJ+VtWJrtrSYap6sMTnzwLpqvqom2O3Agmqus/T61mCMBUubRssfhOWfwhZ+50ZZLteAd1HQ+NuxbsVFSnzN+1jyuIdzPo9hdyCIprUC+HC7k0Y0a0JXeLqWllUUyX5JEGISF1gBU496hMuIs7/LduBC1R1g5vPt2IJwlQVBXmwcRas+BDWfQdF+dCwi5Mouo5yKt65HM4tYNbvKXy5YhdzN6SSX6i0rB/OiG6NGdG9CW0bRvjwixhzLF8liHhgAvA70B1YCtylqpmuz88Fnis1MJEtQBpO09QbqjqhlP3GAmMBmjdv3mvbtm0V/VWMOVbWAUj6FFZMgZ1LQPycpqfuY5ymqBId2wez8vguaQ9frtzFL5v2U6TQoVEEI7o34cJuje2RWeNzvkoQCcCvQD9VXSgiLwIZqvpP1+evARtV9dlSjm+iqrtEpAHwPXCnqs4t65p2B2Eq3b4NsGIyrPgIMpIhuC60H+Y8MttmAAQdTQB7D+XwzcrdfLlyN0u3pQHQtkEdLujYgIEdG9KzeRT+1mdhKpmvEkQj4FdVbeFaPwe4X1X/ICIBwE6gl6ome3Cuh4HDqvqfsvazBGF8pqjIGXi38iOnYzs7DQJCnTuLjhdCu6EQdrT8aXJaFjNWp/DD2hQWbj5AQZESGRbI+e0bMKBjA85tF0vdkMAyLmhMxSgrQQR466KqukdEdohIe1VdBwzAaW4CGAisLS05iEg44Keqh1zvBwMndGQbU2X4+UGr/s5SmO/MKLv2K1j7Naz7GsQfWvSDDiOgwx9oGhXHzWe35OazW5KRk8/c9an8sGYvP67by7TfdhLgJyS2jGZAx4YM6NCAFjbLrPEBbz/FFI/zmGsQsBm4UVXTRORdnLuL10vs2wSYqKrDRaQVMM31UQDwoao+drLr2R2EqXJUYdcyWPOVkzD2rXe2N+np9Fe0vgAaxxfX1C4sUn7bnsasNXuZvSaFDXsPA9A8Oox+bWI4s3V9zmwVQ2yE1a8wFcMGyhlTVaSuh7VfOglj1zJnW2i0c+fR+gJodf4xc0Jt35/FD2tTmL9pP79u3s+hnAIA2jeM4MzWMfRrU58+raKtOcqcMksQxlRFh/fC5jmw6QfY/CMc2u1sj2kLrc93kkWLsyGkLgAFhUWs3pXB/E37+GXTfhZvPUBOfhF+Al2bRnJW6xjObBVDj+aRRFjCMB6yBGFMVacKqWth049Owtg2H/KzwC8AmvZ2EkXzvs77kHoA5BYU8tv2gyzYtJ8FG/exfMdBCooUEecOo0fzKHqd4SwtYsJsoJ5xyxKEMdVNQS7sWOgkjM0/wu6VoIWAOAP0mvc9urjmicrMLWDZ9jSWbTvI0u1p/LY9rbhJKjo8iJ7NI4uTRvemkYQG2ZxRxhKEMdVf7mFnUN72X50leTHkOR3Y1G16NFk0S4QGncA/kKIiZWPqYZZuS2PZtjSWbk9jc2omAP5+QtsGdegSV48uTerSJa4eHRvXJTzYaw82mirKEoQxNU1hgTPj7PZfYfsvzuuRPgz/YGjYGZrEO09INYmH2I4QEERaZh6/7Uhj6bY0knZmkLQznf2ZeYBTAqN1bJ3ihNG5ST06NalLvVDrz6jJLEEYU9OpwsHtzp3F7uWwa7nTLJWb7nzuH+QkjSMJo3E8NOiI+geRkpFL0s50knalk7Qzg9W70tmdnlN86mbRobRrEEG7RhG0a1iHdg0jaB1bh5BAa6KqCSxBGFMbFRVB2hbY9Zv7pCH+ENMaGnR0mqUadHTuNKJbsS+7kNW7nDuM33dnsCHlEJtTMykocn5f+AmcERNenDDaNnSSR4uYcEsc1YwlCGOM40jS2L0c9q5xLb/DgS0482Li3G3Ub+9KHB2cpBHThvx6zdmals+6lEOsTznMhpRDrEs5xLb9WRS6EocIxEWG0rJ+OC3rh9MiJpyWseG0qh9OXGQoAf5WorWqsQRhjClbXhbsW1ciabiWjBKz4YgfRDaH6NYQ08a1tCK3Xis250WxPjWLLfsyjy6pmRzKLSg+PNBfaBYdRsuYcM6ICadZdCjNosJoFh1G06hQ6yD3EZ/MxWSMqUaCwqBJD2cpKSfdGf19YBPs3wT7NzrLjoXFT1EFAx39g+gY1RKiWjhJpPkZaL1mHAxuwpbC+mzMCGDL/iy2pGaydX8m8zftIye/6JhLRYcH0SwqlKZRYTR1JY+mUaHERYbSqF6IDf7zAbuDMMaUnyocTjmaNI4kkLRtcHAb5GYcu39QhJM4os6AyOZovaYcCmrIHo1mW0E9NmXXYdvBApLTskhOy2ZnWjZ5hccmkDrBATSuF0KjeiE0rhdC43qhxetNIkNpWDeEuiEBNiCwnOwOwhhTsUScKnoRjZxZao+XfdBJFAe3H13StjnLlrlI3mHqAnWBdsAgBOo0gIjGEBeHdmzC4aBYUv1i2FMUyc68OmzJCWdLZjC7MnJZt+cQqYdzOf7v2+AAP2IjgmkQEex6DXG7Hh0eRFCA9YecjCUIY0zFC410lsbdT/xMFXIOQsZuyNgFGTudMRwZO531tC3ItvlE5BwkAmhV8ljxd+qCxzSg6IwGZAfGkBEQzX4i2VsYwZ7COuzMDWVbbhib94awcEsRB7Py3YYYERJA/TpOsogODyImPIiYOkFEhwcXv48KCyIqPIjI0EDCgvxr3d2JJQhjTOUSgdAoZ2nYqfT98jKdJHI4BTL3OpMbHk5xve7F73AK4Yd/JzxzL42LCtyfI6gO2iiG/OAosgOjyPSvx0Gpy0Gtw4HCMPYWhpKSF8qulGDW5ISwIzuIg0WhFHHi3UWQvx/1wgKJCgskMjSIyLBAIsMCiQoLol5YIPVCA6kbEkjd0EDqhgS4XgOpGxpAcED1fPTXEoQxpmoKCof6bZylLEVFzh3J4b2Qtd+17HO9HkAy9xGUtZ+grP3UO7SRJpn7oCC7lGuCIhQF1yU/sB65ARHk+IeTJeEcJowMDeVgUSgHMkPYnx5MSn4w63ODSCsI5TAhZGoomYSQSQiFHE0KwQF+xYkjIiSQiJAA6gS7lpAAIlyvdYIDi9fDgwMID/YnPCiAMNdraKA/fpVYltYShDGmevPzc8q5lijpelIFuU4/SXaak1yy04rXJecg/tkH8c9OIyTnIPVyMiA3BXIynM733EMUjxkB8Hctx1/CP4R8/zDy/ELJ8Qsji1AyC0PIPBzC4YwgDhcFkV4UTEZBIAcLA9mnIWRpMNkEk0UwWRpMDsHkEEi2BpNDENkE4RcUSmhQEOHB/oQFBRAe5E+DusG8enWv0/1JnsCrCUJEInEqynXB+YneBAwBbgVSXbv9XVW/cXPsUOBFnB/9RFV90puxGmNqkYBgiGjoLOVVVAR5h44mjJKJI++wM7Fi3mECcg8RkJdJaN5h6rm2OfukOVO552VCURZIXrl/ExcUBpKXHUxejpM4Mg7WB+aV/7uchLfvIF4EvlPVy0UkCAjDSRDPq+p/SjtIRPyBV4BBQDKwWES+UNXfSzvGGGMqhZ+fU5PDVZfjtBUWQH6mM1jxSOI48lqQA/nZR5cC5zXAtYS51hsFhlZMLMfxWoIQkbrAucANAKqaB+R5+BRAIrBRVTe7zjUFuAiwBGGMqVn8A8C/AhNOBfLmg8CtcJqR3hGR30RkooiEuz67Q0RWisjbIhLl5tg4YEeJ9WTXthOIyFgRWSIiS1JTU93tYowx5hR4M0EEAD2B11S1B5AJ3A+8BrQG4oHdwLNujnV3m+F2yLeqTlDVBFVNiI2NrZDAjTHGeDdBJAPJqrrQtf4J0FNVU1S1UFWLgDdxmpPcHdusxHpTYJcXYzXGGHMcryUIVd0D7BCR9q5NA4DfRaRxid0uAZLcHL4YaCsiLV2d26OBL7wVqzHGmBN5+ymmO4EPXL/kNwM3Ai+JSDxOk9FW4DYAEWmC8zjrcFUtEJE7gBk4j7m+raqrvRyrMcaYEmw2V2OMqcXKms3VpjM0xhjjliUIY4wxbtWoJiYRSQW2neLh9YF9FRhORavq8YHFWBGqenxQ9WOs6vFB1YrxDFV1O0agRiWI0yEiS0prh6sKqnp8YDFWhKoeH1T9GKt6fFA9YgRrYjLGGFMKSxDGGGPcsgRx1ARfB3ASVT0+sBgrQlWPD6p+jFU9PqgeMVofhDHGGPfsDsIYY4xbliCMMca4VesThIgMFZF1IrJRRO73dTzHE5FmIvKjiKwRkdUicpevY3JHRPxddT++8nUs7ohIpIh8IiJrXT/LM30d0/FE5B7Xv3GSiEwWkRAfx/O2iOwVkaQS26JF5HsR2eB6dVfPxdcxPuP6d14pItNcpY+rVIwlPvuriKiI1PdFbCdTqxNEidKmw4BOwBgR6eTbqE5QAPxFVTsCfYHbq2CMAHcBa3wdRBmOlL/tAHSnisUqInHAn4AEVe2CM0nlaN9GxbvA0OO23Q/MVtW2wGzXui+9y4kxfg90UdVuwHrggcoO6jjvcmKMiEgznLLK2ys7IE/V6gRBidKmrpKoR0qbVhmqultVl7neH8L5xea2up6viEhT4A/ARF/H4k6J8rdvgVP+VlUP+jYqtwKAUBEJwKnf7tMaKKo6Fzhw3OaLgEmu95OAiys1qOO4i1FVZ6pqgWv1V5x6Mj5Tys8R4Hngb5RSDK0qqO0JwuPSplWBiLQAegALy96z0r2A8x96ka8DKUVZ5W+rBFXdCfwH56/J3UC6qs70bVRuNVTV3eD88QI08HE8J3MT8K2vgzieiIwEdqrqCl/HUpbaniA8Lm3qayJSB/gUuFtVM3wdzxEiciGwV1WX+jqWMpRW/rbKcLXlXwS0BJoA4SJyjW+jqt5E5EGcJtoPfB1LSSISBjwIPOTrWE6mtieIalHaVEQCcZLDB6r6ma/jOU4/YKSIbMVportARN73bUgncFv+1ofxuDMQ2KKqqaqaD3wGnOXjmNxJOVIV0vW618fxuCUi1wMXAldr1Rvs9f/t3U+ITWEYx/HvT4rUsBAhNEXExt9EsyGULGQrSbKRIkpkIxuahUhJFnYWamKBnX9FQqYmjEwRIUpih8nGY/G+1xy3MzHXXOdqfp+6zXvP3HPPM3/envOec9/nnUU6EXiU+810oEfSlEqjKjHSE0TLL20qSaRr530RcbzqeOpFxMGImB4R7aTf382IaKkz38GWv60wpDJvgOWSxuW/+Wpa7EZ6dhnYmttbgUsVxlJK0jrgALAhIr5WHU+9iOiNiMkR0Z77zVtgcf4/bSkjOkHkG1m1pU37gK4WXNq0A9hCOjN/mB/rqw7qP1Rb/vYxsBA4WnE8v8ijmwtAD9BL6puVlmOQdB64B8yV9FbSdqATWCvpOekTOJ0tGOMpoA24lvvLmRaM8b/gUhtmZlZqRI8gzMxscE4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFWIGli4ePE7yW9Kzy/O4zH2SjpUG4flrSv7vuvBqvwKel61VVUbWQYXXUAZq0kIj6R5kkg6TDwOSKONeFQ+4ENQ9khT6ATcA7YCRxpQlxmP3kEYfaHJH3OX1dKuiWpS9IzSZ2SNkt6IKlX0qz8ukmSLkrqzo+OvH0O8C0iPv7BMdvz+hWnSZPoZpBmM29q2g9qlnkEYdaYBcA8Uhnnl8DZiFiWF3TaBewhrUFxIiLuSJpJmrE/jzQ7vqfu/fbWFeebVmjPBbZFxM7aBkljJE3MIx6zpnCCMGtMd63staQXQK00dy+wKrfXAPPTlSEAxktqA6aSyo8XnSheyspF3GpeR8T9utd/ICURJwhrGicIs8Z8K7S/F55/Z6BfjQJWRER/cUdJ/cCEIRzrS8m2sUB/yXazYeN7EGbNc5VUDBIASQtzsw+Y3eib5pvVU4BXfxOc2e84QZg1z25g9ZmgawAAAGdJREFUqaTHkp4CO/L228AiFa49DdES4H5hWU2zpnA1V7MKSDoJXImI6w3uezkibgx/ZGYDPIIwq8ZRYFyD+z5xcrB/wSMIMzMr5RGEmZmVcoIwM7NSThBmZlbKCcLMzEo5QZiZWakfOAYiX/XtkfoAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#a\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"T = 74\n",
"Ta= 65\n",
"T0 = 85\n",
"K = -0.275\n",
"t = np.linspace(0,15,30)\n",
"\n",
"T_analytical = Ta + (T0-Ta)*np.exp(K*t)\n",
"\n",
"T_numerical = np.zeros(len(t))\n",
"T_numerical[0] = T0\n",
"for i in range(1, len(t)):\n",
" T_numerical[i] = T_numerical[i-1] + K*(T_numerical[i-1]-Ta)*(t[i]-t[i-1])\n",
"\n",
"\n",
"plt.plot(t, T_analytical, label='Temp Analytical')\n",
"plt.plot(t, T_numerical, label='Temp Numerical')\n",
"plt.xlabel('Time(Hr)')\n",
"plt.ylabel('Body Temperature(F)')\n",
"plt.legend();\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"T_Infinite = 65.0 F\n"
]
}
],
"source": [
"import numpy as np\n",
"t= np.linspace(0,10000,150)\n",
"T = 85\n",
"Ta = 65\n",
"K = -0.275\n",
"for i in range(0, len(t)-1):\n",
" T_analytical = Ta + (T-Ta)*np.exp(K*t)\n",
"\n",
"print('T_Infinite = ', T_analytical[len(t)-1], \"F\")\n"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The Time of Death Was 9 : 11 am\n"
]
}
],
"source": [
"import numpy as np\n",
"dT = 85 - 65\n",
"dTd = 98.6 - 65\n",
"dt = 11 - np.log((dTd/dT)/-K)\n",
"\n",
"print('The Time of Death Was', int(dt),':',format(int(60*(dt-int(dt))),'01d'), 'am')"
]
},
{
"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": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"def TempNow(t):\n",
" T= np.array([55,58,60,65,66,67])\n",
" Hr= np.array()\n",
" "
]
}
],
"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
}