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?
Go to file
 
 
Cannot retrieve contributors at this time
{
"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": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.275 1/hours\n"
]
}
],
"source": [
"import numpy as np\n",
"\n",
"T_a = 65 #degrees F ambient temperature\n",
"T_i = 85 #temperature at 11 am in degrees F\n",
"T_f = 74 #temperature 2 hours later in degrees F \n",
"k = -(T_f-T_i)/2 *1/(T_i - T_a)\n",
"print(k, '1/hours')"
]
},
{
"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": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.275, '1/hours')"
]
},
"execution_count": 9,
"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",
" Temp_t1 -- the first known temperature\n",
" Temp_t2 -- the second known temperature\n",
" Temp_ambient -- the ambient temperature ie. the temperature of the surrounds (assumed constant)\n",
" delta_t -- time elapsed between the first and second temperature \n",
" \n",
" Returns\n",
" -------\n",
" K -- the empirical constant\n",
" \n",
" '''\n",
" delta_T = Temp_t2 - Temp_t1\n",
" K = -(delta_T/delta_t)*(1/(Temp_t1-Temp_ambient))\n",
" return K, \"1/hours\"\n",
"\n",
"measure_K(85, 74, 65, 2)\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": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f53d567ac18>"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd1yV9RfA8c+XoWAOXOXMUU4EceDIvW24s7ScmTY0rX6VluUurSzLrExztZypWVq5IvdAxW2audAyHDgBGef3xwM3UMYF7mWe9+t1X3Cf+4zD9Xp4+D7nOV8jIiillMo9XDI7AKWUUhlLE79SSuUymviVUiqX0cSvlFK5jCZ+pZTKZdwyOwB7FCtWTMqXL5/ZYSilVLaya9euCyJS/Pbl2SLxly9fnsDAwMwOQymlshVjzKnElutQj1JK5TKa+JVSKpfRxK+UUrlMthjjV7lDZGQkwcHBhIeHZ3YoSmUrHh4elClTBnd3d7vW18Svsozg4GAKFChA+fLlMcZkdjhKZQsiwsWLFwkODqZChQp2bZNjh3p2rviCf8bcT8zoQvwz5n52rvgis0NSKQgPD6do0aKa9JVKBWMMRYsWTdVfyjnyjH/nii+osetNPM0tMFCCEArtepOdgH/HZzI7PJUMTfpKpV5q/9/kyDP+srvft5J+PJ7mFmV3v59JESmlVNaRIxP/3RKSxPILGRyJyo6WLVuGMYYjR46keR/9+vVjyZIlya7zzjvvJHj+wAMPpOlYY8aMYfLkyWnaVuVOOTLx/2vuuEMZgEumYAZHorKj+fPn07hxYxYsWODU49ye+Lds2eLU4ykVJ0cm/jO1XyVM8iRYFiPgJVf4bf5jSExMJkWmsrrr16+zefNmZs2aZUv8AQEBNG/enEcffZSqVavy5JNPEjdz3bhx4/D396dGjRoMGjSI22e0W7duHV26dLE9X7NmDV27dmXEiBGEhYXh5+fHk08+CUD+/Plt67333nv4+PhQs2ZNRowYAcDMmTPx9/enZs2adOvWjZs3bzr1vVA5V468uOvf8Rl2Yo313y0X+NcU46TP05iTM2nxx69s/bASVQeso3Dh8pkdqkrCi7+8SNA/QQ7dp18JPz5q/1Gy6yxfvpz27dtTuXJlihQpwu7duwHYs2cPBw8epFSpUjRq1IjNmzfTuHFjhgwZwqhRowDo3bs3P/30Ex06dLDtr2XLlgwePJiQkBCKFy/OnDlz6N+/Px06dGDatGkEBd35M/78888sX76c7du3ky9fPi5dugRA165dGThwIABvvvkms2bN4oUXXnDIe6Nylxx5xg9W8i8x5k9cxoZSYsyfNHh0BPVePsbGqu2pez2Ea1P92Lt7bmaHqbKY+fPn06NHDwB69OjB/PnzAahXrx5lypTBxcUFPz8/Tp48CcBvv/1G/fr18fHxYf369Rw8eDDB/owx9O7dm2+++YbQ0FC2bt3Kgw8+mGwMa9eupX///uTLlw+AIkWKAHDgwAGaNGmCj48P33777R3HUspeOfKMPynGxYUmPRZyeM9XFPpxGNVWDGXNkZ9o1WMhLi6umR2eiielM3NnuHjxIuvXr+fAgQMYY4iOjsYYw0MPPUTevHlt67m6uhIVFUV4eDjPP/88gYGBlC1bljFjxiRaSx13hu/h4UH37t1xc0v+v52IJFqe169fP5YvX07NmjWZO3cuAQEB6f6ZVe6UY8/4k1OtVh/yv7CH/fnvps3RNWz+4H5CLh3P7LBUJluyZAl9+vTh1KlTnDx5kjNnzlChQgU2bdqU6PpxSb5YsWJcv349ySqeUqVKUapUKSZMmEC/fv1sy93d3YmMjLxj/bZt2zJ79mzbGH7cUM+1a9coWbIkkZGRfPvtt+n5UVUulysTP0DBwuWp/fIfbK32EA1vXOT6J3WYtXQ6jSatp8KIlTSatJ7le85mdpgqA82fPz/BhViAbt268d133yW6vpeXFwMHDsTHx4fOnTvj7++f5L6ffPJJypYtS/Xq1W3LBg0ahK+vr+3ibpz27dvTsWNH6tati5+fn61Uc/z48dSvX582bdpQtWrVtP6YSmFur0LIiurWrSvOnIjlz6BvyP/D/ygcE8X4qN58E90aMHi6uzKxqw+da5V22rHVfw4fPky1atUyOwynGDJkCLVq1WLAgAGZHYrKoRL7/2OM2SUidW9fN9ee8cd3v18v+rl/ztYYbya4z+ET90/Iz03CIqN5/9c/Mjs8lc3VqVOHffv20atXr8wORSkgl13cTc4fV/PSn1d5NuYn/ue2CO88JxkSOZTDoeUzOzSVze3atSuzQ1AqAT3jj1XKyxPBhc+jO9Lz1pvkMxEsyzOaPnl/JjwyLLPDU0oph9HEH+vVdlXwdLdKOndKVR6OeIftUpWx5mt+++A+jp7bnckRKqWUY2jij9W5VmkmdvWhtJcnBvDwKsGlzvM54tudtuE3YEZzlgVMuOOWfKWUym50jD+ezrVK31nBU/tLQqo8TNGlT9M+4D1mHvuFx3v9RCFPr8wJUiml0knP+O1Q3LsLXsMO8E/h8gw6u5/1Uyqx48RvmR2WcgJXV1f8/Pxsj0mTJiW7/ty5cxkyZEi6jztlyhQ8PDy4cuVKuveVHk8//TSHDh1K935OnjyZ4P6HwMBAhg4dmu792iP+zxC/A+rJkyepUaNGitvPnTuX4sWL2z4DX375pe21efPmUalSJSpVqsS8efPsjikgICBLdV/VM347uRYsSYUX9nDmp2F02v01f8zryIx6Ayh2z//4YM0xzoWGUcrLk1fbVdG6/wyyfM9Z3v/1D4e+956enok2TnOUqKioRFs2zJ8/H39/f5YtW5bg7t70iI6OxtU1da1I4ie59IhL/E888QQAdevWpW7dO8rJnSL+z/DOO+/wxhtvpHofjz/+ONOmTUuw7NKlS4wdO5bAwECMMdSpU4eOHTtSuHDhFPcXEBBA/vz50zzngqPpGX9quLhStuM0bvaYT0k3T3ptn8uWpR9xNjQMAc6GhvH60v16x28GWL7nLK8v3Z9h73358uW5cMGayCcwMJDmzZvfsU5ISAjdunXD398ff39/Nm/eDFgTpQwaNIi2bdvSp0+fO7Y7fvw4169fZ8KECbamcGCdeXbq1In27dtTpUoVxo4dC1hJtWrVqvTt2xdfX18effRRW3uH8uXLM27cOBo3bszixYsJCgqiQYMG+Pr60qVLFy5fvkxUVBT+/v62Xj+vv/46I0eOBKB58+bE3SyZP39+hg8fTp06dWjdujU7duygefPmVKxYkRUrVthiadKkCbVr16Z27dq2s9oRI0awceNG/Pz8mDJlCgEBATzyyCOAlUA7d+6Mr68vDRo0YN++fbb36amnnrIdY+rUqXe8V4sWLeLll18G4OOPP6ZixYq297Bx48YJfobEWl9HR0czcOBAvL29adu2LWFh9lfs/frrr7Rp04YiRYpQuHBh2rRpwy+//HLHelOnTqV69er4+vrSo0cPTp48yfTp05kyZQp+fn5s3Lgx2c9K7969admyJZUqVWLmzJkA/P333zRt2hQ/Pz9q1KjBxo0b7Y47USLitAfwEnAQOADMBzzivfYKIECxlPZTp04dyWpirv4t20fXFxldUBa/+YhUHb5Eyg3/ScoN/0kemLgus8PLlg4dOmT3ug9MXGd7v+M/0vveu7i4SM2aNW2PBQsWiIhIuXLlJCQkREREdu7cKc2aNRMRkTlz5sjgwYNFRKRnz56yceNGERE5deqUVK1aVURERo8eLbVr15abN28meszx48fLuHHjJDo6WsqVKyfnz5+37btEiRJy4cIFuXnzpnh7e8vOnTvlxIkTAsimTZtERKR///7y/vvv2+J89913bfv28fGRgIAAERF56623ZNiwYSIicuDAAalataqsXr1a/Pz8JCIiQkREmjVrJjt37hQREUBWrVolIiKdO3eWNm3ayK1btyQoKEhq1qwpIiI3btyQsLAwERE5evSoxP1f/e233+Thhx+2xRH/+ZAhQ2TMmDEiIrJu3TrbvkaPHi0NGzaU8PBwCQkJkSJFisitW7cSvFd///231K1bV0REunXrJnXr1pXg4GCZO3eujBgx4o6f4a677rJte+LECXF1dZU9e/aIiEj37t3l66+/vuPfI+599/HxkW7dusnp06dFROT999+X8ePH29YbN26c7X2Pr2TJkhIeHi4iIpcvX7b9bPHXTe6z4uvrKzdv3pSQkBApU6aMnD17ViZPniwTJkwQEZGoqCi5evXqHcdN7P8PECiJ5FSnDfUYY0oDQ4HqIhJmjFkE9ADmGmPKAm2A0846vrOZAiXoGT6SoW5LecF1Gb55jvN85DD+lDKcC9W6f2dL6j1O73ufnqGetWvXJhgfv3r1KteuXQOgY8eOeHp6JrrdggULWLZsGS4uLnTt2pXFixczePBgANq0aUPRokUBqx//pk2b6Ny5M2XLlqVRo0YA9OrVi6lTp/LKK68A1jAFwJUrVwgNDaVZs2YA9O3bl+7duwPg7e1N79696dChA1u3biVPnoQTFwHkyZOH9u3bA+Dj40PevHlxd3fHx8fH1pY6MjKSIUOGEBQUhKurK0ePHk3xfdq0aRPff/89YM1XcPHiRdu1jYcffpi8efOSN29e7r77bs6fP0+ZMmVs25YoUYLr169z7do1zpw5wxNPPMGGDRvYuHEjXbt2TfHYFSpUwM/PD7DuqI77OeLr0KEDPXv2JG/evEyfPp2+ffuyfv36RCv6EuuiGtd/qXPnznTu3DnROJL7rHTq1AlPT088PT1p0aIFO3bswN/fn6eeeorIyEg6d+5s+xnSytlDPW6ApzHGDcgHnItdPgV4DeuMP9sq4XUXU6IepU/kCIqYa6zI8xZdXDZSrIC2eHa2Ul6JJ9GklqeXm5sbMbEztyXWehkgJiaGrVu3EhQURFBQEGfPnqVAgQIA3HXXXYlus2/fPo4dO0abNm0oX748CxYsSDDcc3tiiXue1PLkjnW7/fv34+Xlxfnz5xN93d3d3bZfFxcXW2tqFxcXoqKiAOui9D333MPevXsJDAzk1q1bKR43uQSaWPvr2zVs2JA5c+ZQpUoVmjRpwsaNG9m6davtF2Fy7Nl/0aJFbesNHDjQdud1mTJlOHPmjG294OBgSpUqdcf2K1euZPDgwezatYs6deokeozkPiuJ/ds2bdqUDRs2ULp0aXr37s1XX32V4s+aHKclfhE5C0zGOqv/G7giIquNMR2BsyKyN7ntjTGDjDGBxpjAkJDEJ0/PbHE3fW2K8eGhiInslwpMyfM5Q28N5qONE4kRneLRWeLfcBfH092VV9tVccrxypcvb0sAcWert2vbtm2CC4L2/OUwf/58xowZw8mTJzl58iTnzp3j7NmznDp1CrCmarx06RJhYWEsX77cltxOnz7N1q1bbfuIG9+Or1ChQhQuXNg2Hvz111/bzv6XLl3KxYsX2bBhA0OHDiU0NNTetyKBK1euULJkSVxcXPj666+Jjo4GoECBArYz2Ns1bdrU1lY6ICCAYsWKUbCg/fNhN23alMmTJ9O0aVNq1arFb7/9Rt68eSlUqNAd6ybV+jo5f//9t+37FStW2BqftWvXjtWrV3P58mUuX77M6tWradeuXYJtY2JiOHPmDC1atOC9994jNDSU69ev3/F+JPdZ+eGHHwgPD+fixYsEBATg7+/PqVOnuPvuuxk4cCADBgywzQyXVk5L/MaYwkAnoAJQCrjLGNMHGAmMSml7EZkhInVFpG7x4olPnp7Z4t/0FUJhXvWcwP77+tPbXKDVurcZMKsJZ6/qhV5nuP2Gu9Jeng7ppBp3MTDuETff7ejRoxk2bBhNmjRJslJm6tSpBAYG4uvrS/Xq1Zk+fXqKx1uwYMEdraC7dOlim++3cePG9O7dGz8/P7p162arjKlWrRrz5s3D19eXS5cu8dxzzyW6/3nz5vHqq6/i6+tLUFAQo0aN4sKFC4wYMYJZs2ZRuXJlhgwZwrBhw+x+j+J7/vnnmTdvHg0aNODo0aO2vzZ8fX1xc3OjZs2aTJkyJcE2Y8aMsb1PI0aMSFVZJECTJk04c+YMTZs2xdXVlbJlyyb6iw+Sbn2dnKlTp+Lt7U3NmjWZOnUqc+fOBayZ0N566y3bBdlRo0bZZkeLEx0dTa9evfDx8aFWrVq89NJLeHl50aFDB5YtW2a7uJvcZ6VevXo8/PDDNGjQgLfeeotSpUoREBCAn58ftWrV4vvvv0/zv1ccp7VlNsZ0B9qLyIDY532A/oA3EDdLdBms4Z96IvJPUvtydltmR5Njawlf3IfoW9d5OY8r7bvMomu1lMcfc7uc3JY5LebOnUtgYOAdZYUnT57kkUce4cCBA5kUmXKWMWPGkD9/ftv1mtTIKm2ZTwMNjDH5jDVo1QpYKiJ3i0h5ESkPBAO1k0v62ZGp1BrPwTtxKVmLGbdiuLDwSZ5b3p/rt65ndmhKKeXciViMMWOBx4EoYA/wtIhExHv9JFBXRC4kt5/sdsZvEx1F9LpxuG75mCCieaVQMdr7zmPZzmi94SsResavVNpllTN+RGS0iFQVkRoi0jt+0o99vXxKST9bc3XDte04eGIx3nm9WH7lKifWrdIbvpRSmUrv3M0Ildvi/vw2Tpj7+TTPdCa4zSIvVtmbzvKllMpomvgzSqEydA0byfSoR+jlto7v84yhnLEubegNX0qpjKSJPwPd7VWASVFPMODW/yhjQvgpz0gectmGm/s1Qm5kzXsVlFI5jyb+DBR309G6mDo8HPEOf0ppPsszlVdd3qLOZzVYeXRlZoeY62lbZm3LPH36dHx8fPDz86Nx48apej8++ugjW8O8LC2xBj5Z7ZEVm7Sl1bLdwfLAxHVSfvhP0vSdX+TYVy+IjC4oB8YXkwqjjTz747NyPeJ6ZoeZKVLTpE1ERPYuFPnQW2R0Ievr3oXpjiF+Uy97xG/SZo/IyMhEl/v7+0vjxo1lzpw5qTp+cqKiohy2r9S6vUlbZrm9SZu3t3eK21y5csX2/Q8//CDt2rWz+3jxm/lltNQ0adMz/gzWuVZpNo9oyYlJD/P76+24v/dU6PEd1d08OehalH8DZ+H3hR/bgrdldqhZ275F8ONQuHIGEOvrj0Ot5U6gbZlzT1vm+O0jbty4kWgjths3bvDwww9Ts2ZNatSowcKFC5k6dSrnzp2jRYsWtGjRAoDVq1fTsGFDateuTffu3bl+/brt32n48OHUq1ePevXq8eeffwKwePFiatSoQc2aNWnatOkdx3WYxH4bZLVHTjrjT9KlEyJfNBcZXVBmvV1M8o5xkTfXvSkRURGZHVmGSXDGsmq4yOyHkn6MKy4yuuCdj3HFk95m1fAUY9C2zNqWWURk2rRpUrFiRSlTpowcPXr0jteXLFkiTz/9tO15aGio7f2P+5yEhIRIkyZN5Pp16y/4SZMmydixY23rxbVZnjdvnu29qVGjhgQHB4vIfy2d7aVn/NlR4fLw1K9Q/zmeunWLQ55l+XrD2zSc1ZBPN2yn0aT1VBixkkaT1mvdP0B0ROqW2ymuLXPcI67FsT3Wrl3LkCFD8PPzo2PHjqlqy9yjR48EbZnjxLVl9vT0tLVlBu5oyxy3HJJvy7xhwwYgYVvm2bNn29WWuVmzZom2ZR44cCA+Pj50797drvHwTZs20bt3byDptszFihWztWWOL7m2zE2aNEnx2Pa0ZQYYPHgwx48f591332XChAl3vO7j48PatWsZPnw4GzduTLRB3LZt2zh06BCNGjXCz8+PefPm2ZrvAfTs2dP2Na7hXqNGjejXrx8zZ860NbxzBp16MStxywMPToJyDan4wxCOut/DwH/y8O5fwbjgAfx30xeQs+/4fTD5i6pMqRE7zHObQmWhv+MvkqemLXNiCd6etswAt27domLFirZ+/NmpLXNMTAweHh4pHlcc3JZ59uzZbN26lQ8++CDFY9++/5Rm4OrRo0eiDfAqV67Mrl27WLVqFa+//jpt27Zl1KiEvSdFhDZt2iQYvosv/r9Z3PfTp09n+/btrFy5Ej8/P4KCgmzzMTiSnvFnRdU7wTO/k6dYJebFHOIttyW48d9/AL3pC2g1CtxvS7DuntZyJ9C2zInLiW2Zjx07Zvt+5cqVVKpU6Y51zp07R758+ejVqxevvPKKrU1y/J+7QYMGbN682TZ+f/PmzQQT1SxcuND2tWHDhoB1raJ+/fqMGzeOYsWKJej/70ia+LOqIhVhwBq+imrDQLdVLMozjlL8190i19/05fsYdJhqneFjrK8dplrL00HbMqdOTmzLPG3aNLy9vfHz8+PDDz9MNL79+/dTr149/Pz8ePvtt3nzzTdtx3vwwQdp0aIFxYsXZ+7cufTs2dN2IfvIkSO2fURERFC/fn0+/vhj23v06quv4uPjQ40aNWjatCk1a9ZMzVtjN6c2aXOUbNukzQEaTVqP39XfmOQ+k2hceDnyOdbH1Mbd/RobhzenRP4SmR2iw2iTtoS0LXPOVb58eQIDAylWrJjD9pllmrSp9Hu1XRXWuzaiw60JnJVizM4zmTfcvyXUZSben3mz6KBzyheVUjmXnvFnA8v3nOX9X//gQugVJt01ny7Rv3KzhC+PyXVW/hvE496P8+lDn1I0n+MvAmUkPeNXKu1Sc8avVT3ZQOdapeNV8HSB/UvI9+MwfnTNw0Lfp+hz4Gt+P/U7Mx6ZQYcqHTI11vQSkURvmFFKJS21J/A61JMd+TwKgwIwBUvRY+8Szvg+Tcl8xem4oCP9lvcjNDxtFRqZzcPDg4sXL6b6Q6xUbiYiXLx40a5S2jg61JOdRYbBz8Nh9zxi7m3A+yWrMXLnNEoWKMnAal/wa5BHtprpKzIykuDg4CTr5JVSifPw8KBMmTK4u7snWJ7UUI8m/pxg3yL48UVw9+CPpq/wyIbNRFzqarvpC8DT3ZWJXX2yfPJXSjmOVvXkZL6PwaAAyH8PVX55g5ERHriT8De/3vSllIqjiT+nKF4Znl4HtZ6kX/QSvs3zDsW5nGCVXH/Tl1IK0MSfs+TJB50+Zbz7UHzNX6zK+zqNXPbbXhaXi6z7a10mBqiUygo08edAPg89y2Mx73BZCvC1+yRecluCh1sMrgVX0frr1jz707Nci0i8j4pSKufTOv4cyLqA+yDP/VKK529+xjC3pTxR7Cx3PTmDUTvL8OHWD/n5z5+Z1XEWrSu2zuxwlVIZTKt6coM938DKVyBvAej2JVvc3Xnqh6f44+IfDKw9kPfbvE8hjzs7Gyqlsjct58ztzh+CxX3h4p/QbARhDQczZsN4Jm+dTKkCpXi62heszmZ1/0qp5GniVxBxHX56CfYvgorNoeuXbA/9i17fTdW6f6VyIK3jV5A3P3SdYfWtP70NpjemflQkxaOfSpD0Qev+lcrJNPHnNsZAnb5WzX/e/DCvA12vL8AQc8eqWvevVM6kiT+3KlHDutvXuwuvuC9irvt7FOFqglWiXS7w3f7vtGmaUjmMJv7cLG8B6DaLoJqjaeBymFV5X8ffWFPD5XUzeBX7nSeXPknHBR0JvhqcycEqpRxFE39uZwx+XV5mS/P5RLp4MD/PBIbnX8W7XX3Y8+JMPmz7Iev+Wkf1T6vzReAXxMidQ0JKqexFq3rUf8KvwooX4NByqNQWOk+Hu4py/NJxBv44kN9O/kazcs2Y2WEmlYpWyuxolVIp0KoelTKPgtB9Ljw0Gf4KgC+awOlt3FfkPtb1WcfMDjMJ+icI3+m+vLf5PaJiojI7YqVUGjj1jN8Y8xLwNCDAfqA/MB7oANwCjgP9RSTZKaP0jD8TnNsDi/tB6BloPRoavgAuLpy7do4hq4aw7MgyqufvS77wHly4Fq03fSmVBaXpBi5jTB7gIaAJUAoIAw4Aq0TkSAoHLA1sAqqLSJgxZhGwCjgHrBeRKGPMuwAiMjy5fWnizyThV+CHIXB4BVRuD50/h3xFEBFGrlzOt5sEQ17b6nrTl1JZS6qHeowxbwLbgRbAXmAesAKrsdsUY8wvxpgaKRzXDfA0xrgB+YBzIrJaROLGCLYBZVL906iM4VEIHvsKHnwP/lwHXzSFMzsxxvD7gUIJkj7oTV9KZRfJdefcLyITknjtPWNMSaBsUhuLyFljzGTgNNZfCqtFZPVtqz0FLExse2PMIGAQwL333ptMmMqpjIH6z0CZutbQz5z20Hos50LLA+aO1c+G3uRS2CWKeBbJ6EiVUnZK7uLuiuQ2FJG/RWRHUq8bYwoDnYAKWMNEdxljesV7fSQQBXybxP5niEhdEalbvHjx5EJRGaF0HXhmA1RqB6tHMi/fxxTk+h2rRZkQqn1ajYUHFuqNX0plUckl/l1x3xhjPkrDvlsDJ0QkREQigaXAA7H76ws8Ajwpmh2yD8/C0ONbaPcOjWU3P+cdSU3z538vu7vyarvK3FvoXnp834MO8ztwKvRUJgaslEpMcok//t/xTdOw79NAA2NMPmOMAVoBh40x7YHhQEcRuZmG/arMZAw0HIzLU79QOJ8bS/KOo7/rz5Qu5MHErj4Ma96IbQO2MaXdFAJOBuD9mTcfbfuI6JjozI5cKRUrucSfrjNxEdkOLAF2Y5VyugAzgGlAAWCNMSbIGDM9PcdRmaSsP/le2IJ75TaMdv+azRXn0LnqXQC4urjyYoMXOfj8QZqXb85Lv75E/S/rs+fvPZkctFIKkinnNMbcBI5gnflXif2e2OciIrUzJEK0nDNLE4Etn8DaMVCojHUDWOna8V4WFh9azNCfhxJyM4Su5SYSfK42/1yJ0Np/pZws1XX8xpj7ktuhiBx3UGwp0sSfDZzeDkv6w40QaPs21BtoDQvFuhx2mT7ffcLeY9464YtSGSTVdfwicjy5h3PDVdnOvfXh2U3WzF4/v2qVfoZfsb1c2LMwl0Ia64QvSmUByd3A9Zsx5jljTKnblrsZY5oaY2YZY/o7P0SVbeQrAj0XQusxcPhH+KIZnAuyvZzUxC5nQ2/qxV+lMlByF3cfBtyBZcaYYGPMPmPMUeAvrJ47n4vInIwIUmUjLi7Q+CXotxKiImBWG9j5JYhQyssz0U2iTAgNZzXUi79KZZDkhnpuishUEakP3If1i+/yw0EAACAASURBVOABEblXRPqLiA66q6SVawjPboTyTWDl/2DJU7zesgye7q4JVvN0d+GJhvk5deUUdWfW5eVfX+ZaxLVMClqp3MGutswiEiEiZ0TkgrMDUjnIXcXgySXQ8i04tJxHtvXg01bulPbyxAClvTyZ2NWX9zp258jgIwysPZAp26ZQ/bPqLDu8TO/8VcpJdCIWlTFOboIlAyA8FB58F2r3TVD1E2frma08u/JZ9p3fxyOVH+GTBz+hvFf5jI9XqRwgTW2ZswpN/DnE9RBYOhD++g18HoNHpkDe/HesFhUTxcfbPmZ0wGjcIx6gtMtz3AjLo3X/SqVSumbgMsaUMca0iP0+rzHmLkcHqHKB/MWh1/fQYiQcWAIzmsP5g3es5ubixv8e+B8fN9+C163nuR6WBwHOhobx+tL9LN9zNsNDVyonSTHxG2OewurU+WXsonLAD84MSuVgLq7Q7DXo84NV5z+zFez+2roD+DazN1xAxD3BsrDIaCb9ciijolUqR7LnjH8o0AC4CiAiR4G7nRmUygUqNLVu+CrrDyuGwPLn4NaNBKskVff/95VwPt3xqdb+K5VG9iT+cBG5FffEGONKYjNwKJVaBe6B3suh2QjYuwBmtIB/D9teTqru393tOkN+HkL9L+uz8+zOjIpWqRzDnsS/2RjzGuARO86/EPjJuWGpXMPFFVq8Dr2XQdglmNkSgr4D4NV2VRKp+3fl/a5NWNBtAeeunaP+l/V59qdnuXjzYmZEr1S2ZE/ifw24htWdcxiwDhjpzKBULnRfC2vop3Qda9hn+WA6exdmYlef2+r+fehSuwyP13icI0OO8GKDF/ly95dUmVaFWbtnESMxmf2TKJXlJVvOGTusM1tE+mZcSHfScs5cJDoKfp8EGyZD8arw2DwoXiXZTfad38fgVYPZdHoT9UvXp8d9k1m8/RbnQsO0BFTlamkq5xSRaKCkMcY9ufWUchhXN2j5plX2eeNfa9x/78JkN/G9x5cN/TYwr/M8Tv9zNx/+ep6zoWFaAqpUEuwZ6vkL2GiMed0YMzTu4ezAVC53fytr6KdkTVg2CFa8AJGJV/kAGGPoU7MP5VyHautnpVJgT+IPAdYA+YDi8R5KOVfBUtD3R2j8Muz+yqr5v3As2U3+uRKR6PKkSkOVyo3cUlpBRN7KiECUSpSrG7QeDeUegKWDrLt9O3wMPo8munopL0/OJpLkI82/PPPjM7zd6m2K5Svm5KCVytrsuXN3jTFm9e2PjAhOKZtKbayhn3tqwPcD4McXITL8jtUSKwH1cHehQdV/mbVnFpU/qcznOz/Xm79UrmbPUM+bwFuxj7exyjr3OjMopRJVqDT0+wkaDYNdc+DL1nAx4SygnWuVvqMEdFJXX5b0HU7Qs0H4lfDj+VXPU3dmXTaf3pw5P4dSmSxN3TmNMb+LSDMnxJMoLedUd/jjF1j+rFX+2XEq1Ohq12YiwuJDi/nf6v8RfDWYXr69eLf1u5QqUCrljZXKZtLcltkYUzDeUxegDta0i5UdG2LSNPGrRIWegSVPQfAO8H8a2r4N7h4pbwfcuHWDiZsm8v6W93F3cefxipM5fKISf18J19p/lWOkJ/GfAQSrP08UcAIYKyK/OyPQxGjiV0mKjoS1Y2DrNKv0s/tcKFLR7s2PXzpOvwWfcup0owRloJ7urkzs6qPJX2Vr6enHXzF2nt2yIlJBRFoCOjiqsgZXd2j3NvSYD5dPwhfN4JD9XcPvK3IfMVcf0dp/lavYk/i3J7Jsh6MDUSpdqj4Ez2yEYpVgUR9Y9RpEJV7Tf7ukavzPht7kasRVR0apVJaQZOI3xtxtjKkJeBpjfIwxvrGPxlg3cymVtRQuB/1/gfrPwY4vYHY766+AFCTV/jnKhFD5k8rM2TNHm7+pHCW5M/6HgWlAGeAz4NPYxxtYpZ1KZT1ueeDBSfD4N3DxL5jeFA7/mOwmSbV/frnt/VQsXJGnVjxFvZn12HJmizMjVyrD2HNx9zERWZRB8SRKL+6qNLl0Ahb3g7+DoMHz0Hqs9YshEcv3nOX9X/+4o6OniDD/wHxeW/MaZ6+dpWeNnrzb+l3KFiqbsT+LUmmQ5qqe2I3bAd7w3xUwEXnHoREmQxO/SrOoCFj9JuyYYfX6f3SONSSUSjdu3WDSpklM3joZg6FbhUn8eao6/1yJ0PJPlWWlp5zzM8ALaArMAboB20TkKWcEmhhN/CrdDi63OnwaA52nWxeD0+BU6CmeWvgZx07U0/JPleWlp5yzsYg8AVyMbdhWH2vcX6nsw7szDAoAr3KwoCf8OtK6ByCVynmVI/xyOy3/VNmaXZOtx301xpSIfV7enp0bY14yxhw0xhwwxsw3xngYY4rENn47Fvu1cBpjVyp1it4HA9ZA3QHWDV9zHrTu/k2l5Mo/g68GpzdKpZzOnsS/yhjjBUwGgoCTwJKUNjLGlAaGAnVFpAbgCvQARgDrRKQS1vy9I9IWulJp4O4Bj3wIj86Gf4/AF03g6K+p2kVK5Z+jfxvNjVs3HBGtUk6RbOI3xrgAP4tIqIgsBioAPiLyhp37d8O6D8ANq/b/HNAJmBf7+jygc5oiVyo9anSDZ36HgmXgu8dgzSi7h36SKv8c/UgtOlTpwLgN46g8rTLzguZp/b/KklKaczcG+Dje8zARuWTPjkXkLNZfCaeBv4ErIrIauEdE/o5d52/g7sS2N8YMMsYEGmMCQ0JC7PphlEqVovfB02ugTj/Y/DHMfQSupDw3b2Ktnyd29WFgIz8WPrqQzU9tpkzBMvT7oR/+M/35/WSGtbVSyi72VPWMBwJFxP4GKNZ2hYHvgceBUGAx1hDRNBHxirfeZRFJdpxfq3qU0+1bDD8OA7e80HWGNfFLOsRIDPP3z+f1da9z5uoZulTtwntt3uP+Ivc7KGClUpaecs7LQCEgAgjD6tIpIlIkhe26A+1FZEDs8z5AA6AV0FxE/jbGlAQCRKRKcvvSxK8yxIVjsKgv/HvQmue3xUhr6sd0uBl5kylbpzBx00Rcw+tTyjxHeISn1v6rDJGecs5igDuQH2uS9WLYN9n6aaCBMSafMcZgJfzDwAqgb+w6fYFU/SWhlNMUqwRPr4VavWHThzCvA1w9l65d5nPPx8imI5nWcjtFo4YSFuGJAGdDwxixdB/L96Q8tKSUo6WY+EUkGugODI/9viTgZ8d227GGdnYD+2OPNQOYBLQxxhwD2sQ+VypryJMPOk2DLl9YrR6mN4E/16V7tzMDzhMTk/Cvh/DIGMb+tJu0zIKnVHrYM9n6NKAF0Dt20U1guj07F5HRIlJVRGqISG8RiRCRiyLSSkQqxX6162KxUhmqZg/rhq+7isM33WD9BEjHBO1J1f5fugFN5jRhW/C2NO9bqdSyZ6jnARF5htgbuWITdeKdrpTKSYpXgYHrwO8J2PA+fNUJrv2Tpl0lVfvvlS+G45eP03BWQx5b/BjHLx1PdD2lHMmexB8ZW88vAMaYooAWJ6vcIc9d0Pkz6PQZBAfC9MbwV0Cqd5NU7f/YDnU59sIxxjQbw8pjK6n2aTVe/OVFLt686KAfQKk72VPV0wfoAtQFZgOPYc25u8D54Vm0qkdlCf8etqp+LhyFah3g3G6r7r9QGWg1CnwfS3bzpFo/x/n72t+MDhjNrD2zKJCnAG80eYOh9Yfi4WbfBPJK3S69bZm9gdaxT9eJyAEHx5csTfwqy4i4bo35n7ltTN7dEzpMTTH52+PgvwcZvnY4K4+t5F73rnhF9eXqTVctAVWplp5yTrD67EQCt1KxjVI5T978cDWREszIMFg3ziGH8L7bm5+e+Im3G/6KudaLKzddbSWgry/dryWgKt3sqeoZCcwHSmG1Y/7OGPO6swNTKsu6kkQHzivB4MDSzJW73UAS1lGERUbz9qr9DjuGyp3sOXvvBfiLyJsiMhKoB/RxblhKZWGFkpqOQmDuw9ZFYAdIqgT032uRPLn0SU5cPuGQ46jcx57Efwqry2YcN+Av54SjVDbQapQ1ph+fu6dV9nnhKHzZyroIfDF9pZlJlYDe5XGLZYeXUfXTqrz0y0tcuHkhXcdRuY89if8mcNAY86UxZibWXbihxpgPjTEfOjc8pbIg38esC7mFygLG+tphKnT+HIbugWYj4Nhq+LQerHoNbqQtMSdVAvpOpwYce+EYvX17M3XHVO6beh9vb3hb5wBQdrOnnHNAcq+LyCyHRpQIrepR2c61fyBgEuz+CtzzQeNh0GCw1RIiFVIqAT0Ucog31r3BD3/8QIn8JRjdbDQDag3A3dXd0T+RyobSVc6Z2TTxq2wr5A9YOxb+WAkFSkKLN6DmE+nu+nm7LWe2MHztcHb/5UrxmAEQXZhSXp681q6qln/mYmku5zTGtDfG7DTG/GuMuWSMuWyM0f46StmjeBXo+R30/8W6KLziBZjeCP74xaEVQA+UfYCXa86nVMwrEF0EMJwLDee174NYtlvnAVYJ2TPGPw14BihN6toyK6XilGtoTfT+2FcQfQvmP27N+BW8y2GHmLz6KJHRJsGyW1Hw6tINbA/e7rDjqOzPnsQfDASJSKSIRMc9nB2YUjmOMVC9EwzeAQ9NhpAj8GVLWNwPLqW/UC6p8s/IqII0mNWArgu7cjjkcLqPo7I/exL/a8CPxphXjTFD4x7ODkypHMvVHeoNtCqAmr4GR3+FafXg5+FwI+3N2ZIq/yxVyIOxzcey9q+11Pi8Bv1/6M+p0FNpPo7K/uxJ/GOBaMALa4gn7qGUSg+PgtByJLyw27oHYMcMmOoHGz+AWzdTvbukyj+Ht6/GqGaj+GvYX7zU4CXm759P5WmVGfbzMM5fP++on0ZlI/aUc+4SkToZFE+itKpH5Qr/HoF1Y+GPVVCglFUB5PcEuLimvG2slMo/AYKvBjPu93HM3jObvG556XTvBE6c8eGfKxHaCC6HSc9k6+8Bv4jIemcFlxJN/CpXObkZ1rwFZ3fB3dWh9Vio1Ma6RuBARy8e5fklczh6oi4u/Nf62dPdlYldfTT55wDp6c45EFhrjLmu5ZxKZYDyjeDpddB9HkSFw3fdrYnfz+526GEqF61M2OU2CZI+WI3g3vvliEOPpbIWexJ/McAdKISWcyqVMYwB787w/HZ48H1rEpiZLWDJU3DJcc3ZkqoEOnsljJm7ZhIZHemwY6msI8XEH1u62R0YHvt9ScDP2YEppQC3PFB/UGwF0KtwZBVM84efR6SrAihOUpVArm5XGPTTIKp9Wo1v9n1DdDommldZjz137k4DWgC9YxfdBKY7Myil1G08CkLLN61fAH49YccXsRVAH1qTwKRRUpVAH3Rtxo89fyR/nvz0XtYbn899WHJoCTGi023nBPYM9TwgIs8A4QAicgnIk/wmSimnKFgSOn4Cz22Fco2sKqBP6sCebyANZ+Wda5VmYlcfSnt5YoDSXp5M7OpDl9pleKTyI+x+ZjeLuy8GoPvi7tSZUYcVf6wgO/T4Ukmzp6pnO9AQCBSR2saYosBaEamVEQGCVvUolaSTm2D1W9bE73dXhzbj4P7WDq8Aio6JZv6B+Yz9fSx/XvqTGgX6436zG5evoyWgWViqyzmNMW4iEmWM6QN0AeoCs4HHgLEissCZAceniV+pZIjAoeVWF9DLJ6BCU+sXQCnHn5tFxUTx6opFLNvmAeS1Lfd0d2FiV19N/llMWso5dwCIyFfAm8Bk4DLQPSOTvlIqBcaAdxerB9CD78H5gzCjOSwZAJdPOvRQbi5u7DhSgvhJHyAsMobxK/c69FjKeZJL/La/FUXkoIh8LCIficiBDIhLKZVabnmg/jMwNAiavAJHVsIndeGX1+Gm4269SaoE9ML1aFp91YpNpzc57FjKOZKbDaK4MeblpF4UEZ12UamsyKMgtHoL/AfAb+/A9umw51to8hLUf/bO+YJTqZSXJ2cTSf4FPaM58O8BmsxpQpuKbRjbfCwNyzZM17GUcyR3xu8K5AcKJPFQSmVlBUtBp2nw7Ga4twGsHWNVAAV9l6YKoDhJlYCO7+jPiWEnmNxmMkH/BPHA7Ado/017tp7Zms4fRDlachd3d4tI7QyOJ1F6cVcpBzix0eoBdG4P3FPD6gF0f6s0VQCl1Azuxq0bfLbzM97b8h4Xbl6gQdEhRIY+xMXrMVoFlIHSUtWzJyNLNpOjiV8pB4mJgUPLrAqg0FNQoVlsBZBzbsa/fus6Ly77hjV77sYkqALSRnAZIS1VPa2cGI9SKjO4uECNbjAkENq/C//shxnN4PuBcNnxk7Pkz5OfwycqJ0j6YDWCm7Byn8OPp+yTZOKPvUNXKZUTueWBBs/CsCBo/DIcXgHT6sKvIx1aAQRJVwGFXI+i9Vet2Xhqo0OPp1JmT8uGNDHGVDHGBMV7XDXGvGiM8TPGbItdFmiMqeesGJRSKfAoBK1HW7OA+T4GWz+1egBt/hgiwx1yiKQawRX0jGb/v/tpOrcpLee1JOBkgEOOp1LmtMQvIn+IiJ+I+AF1sJq7LQPew7rz1w8YFftcKZWZCpWGTp/Cc5uhbH1YMyq2Amh+uiqAIOUqoA/bfsjhC4dpMa8FTec0Ze1fa7UXkJM5LfHfphVwXEROAQIUjF1eCDiXQTEopVJyjzc8uRj6/gh3FYPlz8IXzeDPdWneZVKN4DrXKk0+93y81PAl/hr6F588+AknQk/Q5us2PDD7AVYdW6W/AJwkxSZtDjmIMbOB3SIyzRhTDfgV685gF6zun3dcVTLGDAIGAdx77711Tp1y/IUnpVQyYmLg4FJYN86qAKrY3KoAKlnTaYeMiIpgbtBcJm6aSMjF8twjA4mJKkQpL09ea1dVq4BSKc1z7jrgwHmwzuq9ReS8MWYq8LuIfG+MeQwYJCKtk9uHlnMqlYmiImDnLNjwHoRdBt/HrbkBvO512iG/33WKEUv3ERn936CEu6vwbteadK1T1mnHzWnSM+duej2IdbZ/PvZ5X2Bp7PeLAb24q1RW5pYXGj5v9QBq/BIc+sEa/3dCBVCcD9ccT5D0ASKjDS8vDeCrvV/plJDplBGJvycwP97zc0Cz2O9bAscyIAalVHp5ekHrMfDCLvCJXwE01WEVQHGSKgElujB9l/el8rTKTA+cTkRUhEOPm1s4NfEbY/IBbfjvDB9gIPCBMWYv8A6x4/hKqWyiUBno/Ck8uwnK1LPaQEyrC3sXWNcFHCCpEtDSXp6s6LGCe+66h+dWPkfFqRWZsnUKN27dcMhxcwunJn4RuSkiRUXkSrxlm0SkjojUFJH6IrLLmTEopZykRA3otQT6rIB8RWDZMzCjKRxfn+5dJ1UC+mq7qnSo0oGtA7aypvcaKhetzMurX6bcR+WYsGECoeGh6T52bpAhVT3ppRd3lcribBVAYyH0NNzX0moCV9I3zbtMqRFcnC1ntvDOxndYeWwlRWnH3TKQsHAPbQZHJlb1OIImfqWyiagI2Pkl/P4ehF+JrQAa6dQKoDgfB2zm49UXiIn5b5qRvG6Gd7vVzLXJPzOrepRSuYVbXmg42OoB1GgoHFxmzQK2+i2rFNSJFm2LSJD0ASKihOHLNnMo5JBTj53daOJXSjmeZ2HrZq8XdlndQLd8Ah/7WV8dXAEUJ6lKoPBbnnh/5k3nBZ3ZHrzdKcfObjTxK6Wcx6ssdPkcnt0IpevA6jdhmj/sW+SwCqA4SVUClSzkwaimo9hwagMNZjWg5byWrD6+Ole3g9DEr5RyvhI+0Hsp9F5u3Q+wdKA1D8Dx3xx2iKQqgUa0r87YFmM5/dJpPmj7AX9c/IN237Sjzow6LDywkOh0NqHLjvTirlIqY8XEwIElsG48XDkN97WCNmOtXw7pZE8lUERUBN/u/5Z3N79L8PmSFI8ZANGFKVXIg9faV8tRF4K1qkcplbVEhlsVQBvetyqAavaAFiOt4aEMsHTXGYYv3Utk9H9zDru5xjC2UxWerFc5Q2JwNq3qUUplLe4e8MAQqwLogRfgwFKrB9CaURDm/BuxPlhzLEHSB4iKdmHEsi28svoVgq8GOz2GzKKJXymVuTwLQ9vxsRVAXa3eP1P9YMs0674AJ0mqCshVivHRto+o8HEF+i3vx4F/DzgthsyiiV8plTV4lYUu0+GZDVCqFqweafUA2rfY4RVAkFw/oHz8OfRPBvsPZvGhxfh87sPD3z1MwMmAHFMJpIlfKZW1lPSF3sush0chWPo0zGwOfwU49DBJ9wOqQnmv8nzU/iNOv3iacc3HsfPsTlrMa0G9L+ux6OAiomKiHBpLRtOLu0qprCsmBvYvhvXj4coZuL+11QOoRA2H7N7efkBhkWF8tfcrJm+dzLl/S1M8egDEFKZkIQ+GZ+FKIK3qUUplX5HhsHNmbAXQVfB7Alq8YbWIzkCJVQK5ukQzquP99G1QPUNjsYdW9Silsi93D6vyZ2iQVQm0f7FVAbR2TIZUAMVJrBIoOsaVkT9s5+kVT2ebnkCa+JVS2Ue+ItB2glUBVL0zbPrIqgDa+plTK4DiJFUJ5CbF+Xb/t3h/5s1D3z7Eur/WZekLwZr4lVLZj9e90PULeOZ3KOkHv75u9QDav8QpFUBxkqsEOv3iacY2H0vguUBaf92a2jNq8/Xer7kVfctp8aSVjvErpbK/P9fBmtFwfr/1i6DteKjQ1OGHWb7nLK8v3U9Y5H/9fTzdXZnY1cd2gTc8Kpxv933Lh9s+5FDIIUoXKE27Um9y6MT9/HMlIkMniNGLu0qpnC0m2ur6uX4CXA2GSm2tyeHv8XboYeytBBIRfvnzF8b+spJzZ1vggofttdt/WTiLJn6lVO4QGQ47voANH0DEVfB7MrYCKHNKLhtNWs/ZRK4NFM1vCBz5IMaYRLZyDK3qUUrlDu4e0GiY1QOo4WDYvwg+qQ1rx1rN4DJYUheEL1yPpvaM2swLmkdEBlyYjk8Tv1IqZ8pXBNq9DUMCoVpH2PShNQvYts8hKuMuuCZ1QbhQvhhuRd+i3w/9KP9xecb/Pp6QGyEZEpMmfqVUzla4HHSbCYN+t3r+/zICPvWHA987tQIoTlKtIcZ1qMuB5w7wa69f8Svhx6iAUZSdUpanVzzN/vP7nRqTjvErpXIPETgeVwF0wGoG12Y8VGji1MPac0H4cMhhpm6fyry98wiLCsO/8HNEXnmI0BsmzZVAenFXKaXixETDvoWw/u3YCqB2sRVAmd924VLYJV79YRFrg+7GkNe2PC2VQHpxVyml4ri4Wv1+Xgi0mr6d3gbTG8EPQ+DquUwNrYhnEY6crJwg6QOERUbz/q9/OOQYmviVUrmXuyc0ftGqAKr/nPVXwNTasG5cplQAxUmqEiip5amliV8ppfIVgfbvwJCdUO0R2PgBTK0F27/I0AqgOElVAiW1PLU08SulVJzC5aHblzAowLrj9+fX4NN61nzAGXg9NLlJYhxBE79SSt2uVC3oswKe/B7c88GS/vBlKzi5KUMO37lWaSZ29aG0lycGKO3l6dAWD1rVo5RSyYmJhr0L4Le34epZqNzeqgC6u1pmR5YirepRSqm0cHGFWk9acwC0HgOntsDnD8CKF+Dq35kdXZpo4ldKKXu4e0Ljl6xZwOo/C0HzrQvA68Zb00FmI05L/MaYKsaYoHiPq8aYF2Nfe8EY84cx5qAx5j1nxaCUUg53V1FoP9GqAKr6MGycbM0Ctn1GplQApUWGjPEbY1yBs0B9oCIwEnhYRCKMMXeLyL/Jba9j/EqpLOvsblgzCk5uhCIVodUoa1pIJ7Zbtldmj/G3Ao6LyCngOWCSiEQApJT0lVIqSytdG/r+CE8sBte8sLgffNkaTm7O7MiSlFGJvwcwP/b7ykATY8x2Y8zvxhj/xDYwxgwyxgQaYwJDQjKmValSSqWJMVC5LTy3GTp9arV9mPsQzO8JIY5ps+BITh/qMcbkAc4B3iJy3hhzAFgPDAP8gYVARUkmEB3qUUplK7duwvbPYdNHcOs61OoNzV+HgiUzNIzMHOp5ENgtIudjnwcDS8WyA4gBimVAHEoplTHy5IMm/7MqgOo9A0HfWbOArZ+QJSqAMiLx9+S/YR6A5UBLAGNMZSAPcCED4lBKqYx1V1F4cBIM2WHd+LXhfasEdMdMiI7MtLCcmviNMfmANsDSeItnAxVjh3wWAH2TG+ZRSqlsr0hF6D4HBq6H4lVh1SvwaX049EOG9gCKoy0blFIqI4nAsdXWLGAhh6GMvzULWLmGDj9UZpdzKqWUgtgKoHZWBVDHaXAlGOa0h/lPQMjRDAlBE79SSmUGF1eo3Rte2A0t34ITG+CzBvDji3DtH9i3CKbUgDFe1td9ixx2aB3qUUqprODGBevi784vAQMIxET997q7J3SYCr6P2b1LHepRSqms7K5i8OC7MHgHuLolTPoAkWHWlJAOoIlfKaWykqL3QWR44q9dCXbIITTxK6VUVlOoTOqWp5ImfqWUympajbLG9ONz97SWO4AmfqWUymp8H7Mu5BYqCxjrayov7CbHzSF7UUop5Vi+jzks0d9Oz/iVUiqX0cSvlFK5jCZ+pZTKZTTxK6VULqOJXymlcpls0avHGBMCnErj5sXImhO9aFypo3GljsaVOlk1LkhfbOVEpPjtC7NF4k8PY0xgYk2KMpvGlToaV+poXKmTVeMC58SmQz1KKZXLaOJXSqlcJjck/hmZHUASNK7U0bhSR+NKnawaFzghthw/xq+UUiqh3HDGr5RSKh5N/Eoplctk68RvjGlvjPnDGPOnMWZEIq8bY8zU2Nf3GWNq27utk+N6MjaefcaYLcaYmvFeO2mM2W+MCTLGOHSiYTviam6MuRJ77CBjzCh7t3VyXK/Gi+mAMSbaGFMk9jWnvF/GmNnGmH+NMQeSeD2zPlspxZVZn62U4sqsz1ZKcWX4Zyt232WNMb8ZYw4bYw4aY4Ylso7zPmMiki0fgCtwHKgI5AH2AtVvW+ch4GesmYsbANvte8JrnQAAB3pJREFU3dbJcT0AFI79/sG4uGKfnwSKZdL71Rz4KS3bOjOu29bvAKzPgPerKVAbOJDE6xn+2bIzrgz/bNkZV4Z/tuyJKzM+W7H7LgnUjv2+AHA0I/NXdj7jrwf8KSJ/icgtYAHQ6bZ1OgFfiWUb4GWMKWnntk6LS0S2iMjl2KfbAMfMp5bOuJy0raP33ROY76BjJ0lENgCXklklMz5bKcaVSZ8te96vpGTq+3WbDPlsAYjI3yKyO/b7a8BhoPRtqzntM5adE39p4Ey858Hc+cYltY492zozrvgGYP1WjyPAamPMLmPMIAfFlJq4Ghpj9hpjfjbGeKdyW2fGhTEmH9Ae+D7eYme9XynJjM9WamXUZ8teGf3ZsltmfraMMeWBWsD2215y2mcsO8/AZRJZdnttalLr2LNtWtm9b2NMC6z/nI3jLW4kIueMMXcDa4wxR2LPWjIirt1YvT2uG2MeApYDlezc1plxxekAbBaR+Gdwznq/UpIZny27ZfBnyx6Z8dlKjUz5bBlj8mP9snlRRK7e/nIimzjkM5adz/iDgbLxnpcBztm5jj3bOjMujDG+wJdAJxG5GLdcRM7Ffv0XWIb1Z12GxCUiV0Xkeuz3qwB3Y0wxe7Z1Zlzx9OC2P8X/397ZhkhVhXH89y9fIhMFiwok1kypjFIyP2hCpVFUhFFg2IsEfdCyyLB36OVLKH5IK0TCJBfUSjGzMJU0X9DSKE0Nbd3MYiHSIEFqtXb36cM5047TnZ2Zde4dl3l+MHDn3Oee88yZZ54595x7/zfF/ipFLWKrLGoQWyWpUWxVQuaxJak3IekvNbNVCSbpxVgaCxdZvAhnK4eBIXQucIwosLmT0xdHdpV7bMp+XQY0A2MLyvsB/fO2dwC3Z+jXJXTe1DcG+CX2XU37K9oNIMzV9suiv2KdDRRfrMw8tsr0K/PYKtOvzGOrHL9qGFsCGoF5XdikFmM9dqrHzNokzQDWE1a5F5vZ95Kmxf0LgbWElfFm4C/gka6OzdCvl4FBwAJJAG0W1PcuBj6KZb2AZWa2LkO/7gOmS2oDWoH7LURarfsL4B5gg5n9mXd4av0laTnhSpQLJbUArwC983zKPLbK9Cvz2CrTr8xjq0y/IOPYiowDHgL2SdoTy14k/HGnHmMu2eA4jlNn9OQ5fsdxHKcbeOJ3HMepMzzxO47j1Bme+B3HceoMT/yO4zh1hid+pypIGijpsQqPGSJpp6RDkj6Q1KeCYxskTeli/zpJxyV9WlA+IyoaWryBKHUkTcqpUcbtq4vYTZP0cAX13lT4+dJEUh9JWyX12MvAnYAnfqdaDAQqSvzAHOANMxsG/EGQGCiXBqBo4gfmEq6TLmQ7MBH4uYK2zpRngQVxexKQmPjNbKGZNWbmVRGKJXYLgmAbgcnZeuRUG0/8TrWYDQyN2uVzSxkr3BlzC7AyFi0hJMVCuwZJ2yR9G19j89obH9ubWXicmW0ETiSU7zazIyV8S2wzjrC3SPpQUpOk2Qr697sUdNuHJtQ1HDhlZr/Heu4G5ka/hxbYvippVtzeLGlOrLtJ0vgi7l4gaaWkg5KWxn5F0gRJu6NfiyX1jeVHcmc6kkZL2pzX9juSNgCNkkbEtvcoaMEPi+2tBh7oqv+csx8/ZXOqxfPANWY2EkBSf2BbEdspwFHguJm1xbJiCoNHgVvN7GRMPsuB0bG9WWZ2VxU/Q6k2Aa4DriLc4n8YWGRmYxQepPEE8FRBXeMIAmWY2Q5Jawi69CspTa9Y9x2EO04nJtiMAkYQtFq2A+MUHhryHjDBzJokNQLTgXkl2rseuNHMWiW9Bcw3s6VxCu7caLMfuKEM352zGE/8TipY0BgfWWy/pIuSDkso6w28LWkk0A4Mr46HXdJVm1+b2a8Akn4ENsTyfcDNCXVdChzrph854a5vCFNbSewys5boz55odwL4ycyaos0S4HFKJ/41ZtYat78EXpI0GFhlZocAzKxd0t+S+sfv2OmBeOJ3UqGMEf8BwoMlesVRfzGFwZnAb4SR9jnAyRTcraTNU3nbHXnvO0j+PbUSRMC6Q67u9iJ1F/qTs0uS7c3RRucU73kF+/7TqjGzZZJ2EoTC1kt61Mw2xd19yeZ7cFLCE79TLU4QHiEHlB7xA0j6giDe9T4wFfg4wWwA0GJmHZKm0jnlcFp7VaZYm93hAPBg3vs0/c5xEGiQdIWZNRMWubfEfUcIUzqfAfcWq0DS5cBhM3szbl8LbJI0CDhmZv+k+QGcdPHFXacqWNB9367wwOqSi7uR54CnJTUTFCXfTbBZAEyV9BVhyiU3Kt1LUHb8LmlxV9I2YAUwQVKLpNti+ZMKKo2Dgb2SFlXQZnfYCozKLboS/uSeiQuv/1sMrgZmdpKg5LhC0j7C2UhOifI1YH7sn/YuqpkM7I/TR1cSJIQhTGetTcNvJztcndNxUkbSfOATM/u81r6cKZJWAS+Y2Q+19sXpPj7id5z0eR04v9ZOnCnx6p7VnvR7Pj7idxzHqTN8xO84jlNneOJ3HMepMzzxO47j1Bme+B3HceoMT/yO4zh1xr/9AJQ5X9V/TgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"N = 100\n",
"t = np.linspace(0, 2, N)\n",
"dt= t[1]-t[0]\n",
"T_analytical = np.zeros(len(t))\n",
"T_0 = 85\n",
"T_a = 65\n",
"K = 0.275\n",
"\n",
"for i in range(0, len(t)):\n",
" T_analytical[i] = T_a + (T_0-T_a)*np.exp(-K*t[i])\n",
"\n",
" \n",
"plt.plot(t, T_analytical, label=\"Analytical\", color=\"green\")\n",
"plt.xlabel('t=0 at 11 am (t in hours)')\n",
"plt.ylabel('Temperature (F)')\n",
"\n",
"N = 25\n",
"t = np.linspace(0, 2, N)\n",
"dt = dt= t[1]-t[0]\n",
"T_E_10 = np.zeros(N)\n",
"T_E_10[0] = 85\n",
"for i in range(1, len(t)):\n",
" T_E_10[i]= T_E_10[i-1] - K*(T_E_10[i-1]-T_a)*dt\n",
"\n",
"plt.plot(t, T_E_10, 'o', label= \"Euler Approximation with 25 steps\")\n",
" \n",
"\n",
"N=3\n",
"t = np.linspace(0, 2, N)\n",
"dt= t[1]-t[0]\n",
"T_E_2 = np.zeros(N)\n",
"T_E_2[0] = 85\n",
"for i in range(1, len(t)):\n",
" T_E_2[i]= T_E_2[i-1] - K*(T_E_2[i-1]-T_a)*dt\n",
"\n",
"plt.plot(t, T_E_2, 'o-', label=\"Euler Approximation with 3 steps\")\n",
"\n",
"plt.legend(loc='best') "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-1.8865228851460631 hours\n"
]
}
],
"source": [
"T = 98.6\n",
"t = -np.log((T-T_a)/(T_0-T_a))/K\n",
"print(t, \"hours\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"3b. As t approaches infinity the temperature approaches the ambient temperature, 65 Farenheit.\n",
"\n",
"\n",
"\n",
"3c. Time of death is 9:07 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",
"\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": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'Ambient Temperature (F)')"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU5dn/8c/FEgJhJ2GTNRA2ZUdEDO61Vttad6VaH/URbRX3Wq396VOfWrVad7vwuFdEcYHWDbFWK0gFCftO2MOasIYlCclcvz/mYCOFMEAmJ5n5vl+vvDJzMnPOd6JcOXPPfa7b3B0REUketcIOICIiVUuFX0Qkyajwi4gkGRV+EZEko8IvIpJk6oQdIBbp6eneqVOnsGOIiNQoOTk5Be6esf/2GlH4O3XqxPTp08OOISJSo5jZqgNt11CPiEiSUeEXEUkyKvwiIklGhV9EJMmo8IuIJBkVfhGRJKPCLyKSZFT4RUSqmeLSMqbkFvDIhEVs3FFU6fuvERdwiYgkMndn8cZCJi8t4IulBUxbsZmivRHq1DIGdWxGq8aplXo8FX4RkRBs2lHE5NwCJi8tYFJuAfmFxQBkZqRx2fEdyO6azpAuLWhYr/LLtAq/iEgV2FNSxtQVm6OFfmkBizcWAtCsQV2yszIY1jWd7Kx02jatH/csKvwiInEQiTjz1+1gUm4+k5cWMH3lVkrKIqTUrsXxnZvxo/49GJaVTq82jalVy6o0mwq/iEglWbttD5OX5jNpaQFf5hawdfdeAHq0bsRVQzuSnZXB4E7NqZ9SO9ScKvwiIkeosGgvXy3f8k2xX16wC4CMRvU4rUdLhmWlc1LXdFo2qtwPZ4+WCr+ISIxKyyLMztvO5KUFTM7NZ+bqbZRGnNS6tTihcwuGn9CBYVkZdGvVELOqHb45HCr8IiIH4e6s2rybSbkFTF6az5RlmyksKsUMjmvbhBEnZ5Kdlc7Ajs2oVyfc4ZvDocIvIlLOtt0lTFm2mUnBWf2aLXsAOKZpfc7t3YbsrHSGdkmneVpKyEmPnAq/iCS1ktIIM1Zv/WY+/dy8bUQcGtarw5DMFlw3LJPsrul0Tk+r1sM3h0OFX0SSiruTu2lncEZfwFfLN7O7pIzatYy+7Zow8vQshmWl07d9U+rWTsyuNir8IpLwCnYW82Vu9MKpyUsL2BD0v+nUogEXDmhHdlY6J3ZpQePUuiEnrRoq/CKScIr2lvH1yi3fXCW7YP0OAJrUr0t2cIVsdtd02jdvEHLScKjwi0iNF4k4CzfsCKZZFjBtxRaKSyPUrW0M7NiMn3+3O9ld0znumCbUruKrZKujuBZ+M2sKPA8cBzhwjbv/y8xGAjcBpcAH7n5XPHOISOLZsL2ISUvzmZwbvUq2YGcJAFktG/LjEzoyLCudwZ2bkxaHJmc1Xbx/I08BE9z9IjNLARqY2WnAeUAfdy82s5ZxziAiCWBXcSlTV2z+Zpx+6aadAKQ3TAmGbzLI7ppO6ybV6yrZ6ihuhd/MGgMnA/8F4O4lQImZ/RR42N2Lg+2b4pVBRGqusogzd+32b9ohzFi9lb1lTr06tRjcuTkXD2pHdtcMerRuVOVNzmq6eJ7xZwL5wEtm1hfIAW4BugHDzOxBoAi4092/3v/JZjYCGAHQoUOHOMYUkepizZbd31w49WXuZrbviTY569WmMdec1JlhWRkM6tSM1Lo15yrZ6iiehb8OMAAY6e5Tzewp4O5gezNgCHA8MNbMMt3dyz/Z3UcBowAGDRr0rZ+JSGLYvmcv/1q2mclB6+KVm3cD0LpxKmf1akV20OQsvWG9kJMmlngW/jwgz92nBvffJlr484B3g0I/zcwiQDrRdwciksD2lkWYtWZbME6fz+y87ZRFnAYptRmS2YKrhnZiWFY6XTKqd5Ozmi5uhd/dN5jZGjPr7u6LgTOABcAy4HTgczPrBqQABfHKISLhcXdWFOxiUjCf/qvlm9lZXEotg97tmvKzU7uQ3TWd/h2akVInMa+SrY7iPatnJDA6mNGzHLga2AW8aGbzgBLgqv2HeUSk5tqyq4Qvg7VkJ+cWsHZbtMlZ++b1+WG/tgzrGm1y1qRBclwlWx3FtfC7+yxg0AF+dEU8jysiVae4tIyclVuD1sUFzFu3HXdolFqHk7qk89NTuzAsK52OLdLCjioBXdkgIofF3VmycSeTgmmW01ZsYc/eMurUMgZ0aMZtZ3YjOyudPsc0oU6CNjmr6VT4ReSQNhUWRYduguGbTYXFAGRmpHHp8e3J7prOkC4taKirZGsE/VcSkf+wp6SMaSu3MGlJtCXCog2FADRrUJfsrAyGBY3O2jatH3JSORIq/CJCJOIsWL+DL5ZG59NPX7mVkrIIKbVrcXznZvzi7B4My0qnV5vGuko2AajwiySptdv2fNMOYcqyzWzZFW1y1qN1I64a2pHsrAwGd2pO/RRdJZtoVPhFkkRh0V6+Wr4lWuxzC1ievwuAlo3qcWr3DIYFV8m2bKQmZ4lOhV8kQZWWRZidtz34QDafmau3URpxUuvWYkhmC4YP7sCwrAy6tdJVsslGhV8kgazavIsvgnYIU5ZtprCoFDPofUwTRpycSXZWOgM7NqNeHQ3fJDMVfpEabNvuEqYs2/xNR8s1W6JXyR7TtD7n9m4TbXLWJZ1maSkhJ5XqRIVfpIZxd16espLxs9YxN28bEYeG9epwYpcWXDcsk+yu6XROT9PwjRyUCr9IDfNWTh6/fm8Bfdo1YeTpWQzLSqdv+6bU1VWyEiMVfpEaZOnGQu7/63yGdmnBX649QQuHyxHRKYJIDVG0t4ybXp9Jg5TaPHlpPxV9OWI64xepIX793gIWbyzklWsG07Kx5trLkdMZv0gN8P6cdYyZtpobTunCKd0ywo4jNVyFZ/xm1ga4FBgGtAX2APOAD4CJWkBFJP5Wb97NPe/MZUCHptxxVrew40gCOOgZv5n9H/Ba8JiniK6edTswGfgR8KWZZVdFSJFkVVIa4aYxMzCDpy/vr5k7UikqOuN/1t1nH2D7LGCsmaUCHeITS0QAHpmwiDl52/nTFQNo16xB2HEkQVRU+CtcAN3di4AllRtHRPb5dOFGXpi8gqtO7MjZx7UJO44kkIreN76374aZja2CLCISWL99D3e8NZtebRpzzzk9w44jCaaiwl9+knBWvIOISFRpWYRbxsyipDTCs8P7k1pXDdWkclU01OMHuS0icfT0p0uZtnILT1zal8yMhmHHkQRUUeHva2ZbiJ75NwpuE9x3d28e93QiSebL3AKe+SyXiwa24/z+7cKOIwmqosKvPq4iVSi/sJhb35xFZnoaD5x3bNhxJIFVNMZfz93LDvYFYGYVzi8zs6Zm9raZLTKzhWZ2Yrmf3WlmbmbplfRaRGqsSMS5fewsduzZy3M/HkCDFHVTkfipqPC/b2aPmNnQYM4+AGbWwcyuMrMPge8fYv9PARPcvQfQF1gY7KM98B1g9dHFF0kMf/5iOZOWFnDfD3rRo3XjsONIgquo8J8BfAncAuSa2bZgnP9toBNwnbsfdJqnmTUGTgZeAHD3EnffFvz4CeAu9KGxCDmrtvDYxMWc27sNwwfrmkiJv4O+nwz68Pwt+DoSmUA+8JKZ9QVyiP4ROQNY6+6zK1ohyMxGACMAOnTQPwZJTNt2l3DzmFm0bZrKQxf21qpZUiXi2fijDjAA+KO79wd2Af8D3Avcd6gnu/sodx/k7oMyMtSNUBKPu3PX23PYVFjEs5cPoHFq3bAjSZKIZ+HPA/LcfWpw/22ifwg6A7PNbCXQDphhZq3jmEOkWnr1X6uYuGAjvzi7B33bNw07jiSRuBV+d98ArDGz7sGmM4AZ7t7S3Tu5eyeifxwGBI8VSRrz1m7nwQ8WcnqPllyb3TnsOJJkYpozZmZDgG7u/qqZtQDS3D2WGTkjgdFmlgIsJ9raWSSp7SwuZeSYmTRPS+Gxi/tqXF+q3CELv5n9CjgJ6AK8CqQCrwOH7MXv7rOAQRX8vFOsQUUSgbvzq3FzWbV5F2OuG0LzNF0nKVUvlqGei4BziH44i7uvBTTRWOQIvJWTx/hZ67j1zG6ckNki7DiSpGIp/MXB1E6HQ1+tKyIHlrupkPv/Op+hXVpw42ldw44jSSyWwv+umT0HNDGzq4GJwIvxjSWSWIr2lnHj6Jk0SKnNk5f2o3YtjetLeA45xu/uj5jZ94ASom0XHnT3j+KeTCSBPPD+AhZvLOTlq4+nZePUQz9BJI4qLPxmVhv40N2/C6jYixyB9+es4/Wpq7n+lExO7d4y7DgiFQ/1BF04S4K+OyJymFZv3s0978ylf4em3HlW90M/QaQKxDKPfyfRK20nEszsAXD32+OWSiQBlJRGuGnMDMzgmcv7U7d2PC+UF4ldLIX/78GXiByG301YxJy87fzpigG0a6bJcFJ9xPLh7gtVEUQkkXy6cCPPT17BT07syNnHtQk7jsi3xHLl7lIO0Dff3bvFJZFIDbd++x7ueGs2vdo05pfn9Aw7jsh/iGWop3xrhlTgYqBJfOKI1GylZRFuGTOLktIIzw7vT2rd2mFHEvkPsQz1bNxv02NmNjlOeURqtKc/Xcq0lVt44tK+ZGY0DDuOyAHFMtTTp9zdWkSbrumMX2Q/U3ILeOazXC4a2I7z+7cLO47IQcUy1PNcudulwArg0vjEEamZCnYWc8ubs8hMT+OB844NO45IhWIp/Fe4+6ryG8xMi+CKBCIR57Y3Z7F9z15evWYwDVJiWuZCJDSxXFEy7gDbxld2EJGa6s9fLGfS0gLu/0EverbRRe5S/R301MTMugE9iXbl/GG5HzUmOrtHJOnlrNrKYxMXc27vNgwfrDfCUjNU9J70WOACoCnRKZz7FALXxzOUSE2wffdebh4zk7ZNU3nowt5aQlFqjIMWfncfB4wzs2x31/RNkXLcnbvemc3GHUW8/dOhNE6tG3YkkZjF8inU12Z2PdF3AN8M8bj7iLilEqnmXv3XKj6ev5F7z+lJv/ZNw44jclhi+XD3VaAT8H1gKtFF14vimEmkWpu3djsPfrCQ03u05NrszmHHETlssRT+bu5+D7AzaNh2NnBcfGOJVE87i0sZOWYmzdNSeOzivtTSEopSA8VS+PcG37eZWU+gEdAxfpFEqid35/+Nn8eqzbt46rJ+NE9LCTuSyBGJZYz/BTNrBtwPfAw0AO6LZedm1hR4nug7BAeuITpT6AdE1/BdBlzt7tsOP7pI1Xo7J49xM9dy25ndOCGzRdhxRI5YhWf8wZq7Be6+1d0/c/cO7p7u7n+Icf9PARPcvQfRhdoXAp8Ax7l7H2AJcM9R5BepErmbCrnvr/M5MbMFN53eNew4IkclljV3bz2SHQfr9J4MvBDsq8Tdt7n7RHcvDR72FaBuVlKtFe0t48bRM2mQUpsnL+tHbY3rSw0Xyxj/x2Z2q5m1MbPG+75ieF4mkA+8ZGYzzex5M0vb7zHXAB8dbmiRqvTA+wtYvLGQ31/Sl1aNddG61HyxFP7rgTuAacA8YH7w/VDqAAOAP7p7f6ILtd+974dmdi/Rbp+jD/RkMxthZtPNbHp+fn4MhxOpfB/MWc/rU1dz/SmZnNq9ZdhxRCpFLAuxtD/CfecBee4+Nbj/NkHhN7OriF4XcIa7/8eyjsFxRwGjAAYNGnTAx4jE0+rNu7n7nTn079CUO8/qHnYckUpzyDN+M6tvZneb2R+D+13N7HuHep67bwDWmNm+fzFnAAvM7GzgF8AP3X33UWQXiZuS0gg3jZmBGTx9WX/q1o7lzbFIzRDLdM4XgbnAsOD+OuAtYhubHwmMNrMUYDlwNfA1UA/4JGhq9ZW733CYuUXi6ncTFjEnbzt/umIA7Zs3CDuOSKWKpfBnufvlZnYxgLvvthjbELr7LKJLNZanuXBSrf1j0Uaen7yCn5zYkbOPaxN2HJFKF8v71xIzSyV6ARZm1pnoxVciCWf99j3cMXY2vdo05pfn9Aw7jkhcxHLG/wAwAWhnZq8ApwDXxjWVSAhKyyLcMmYWxaURnh3en9S6tcOOJBIXsczqmWBmOcBQwICfu/umuCcTqWJPf7qUaSu38PglfcnMaBh2HJG4iXVV6BOBk4gO95QB78UtkUgIpuQW8MxnuVw0sB0XDNDF5JLYYpnO+QxwC7AUyAVuDraJJISCncXc8uYsMtPTeOC8Y8OOIxJ3sZzxn060qdq+D3dfBObENZVIFYlEnNvHzmb7nr28es1gGqTE+iZYpOaKZVbPEr7dSK0NsbVsEKn2Rk1azhdL8rnv+73o2SaWFlQiNV8spzdNgIVm9lVw/wRgipm9C+DuF8QrnEg85azayqMfL+ac3q358Qkdwo4jUmViKfwPxj2FSBXbvnsvN4+ZSdumqTx0QR9ivCZRJCHEMp3zUwAza1D+8e6+I465ROLG3bnrndls3FHE2z8dSpP6dcOOJFKlDln4zexa4DdEp3FGiM7ld0DvjaVG+stXq/h4/kbuPacn/do3DTuOSJWLZajnbqCvLtqSRDB/3XZ+8/5CTuuewbXZncOOIxKKWGb1LAc0rCM13s7iUm56fSbN0ury+0v6UUtLKEqSivWM/8tgVk/xvo3ufnvcUolUMnfn/42fx6rNu3j9uiE0T0sJO5JIaGIp/H8CviTakz8S3zgi8fF2Th7jZq7ltjO7MSSzRdhxREIVS+GPuPvNcU8iEie5mwq576/zOTGzBTedruUgRGIZ4//UzK4xswwza7zvK+7JRCpB0d4ybnp9Jg1SavPkZf2orXF9kZjO+K8Kvv+63DZN55Qa4YH3F7BoQyEvX308rRqnhh1HpFqI5QKu9lURRKSyfTBnPa9PXc31p2RyaveWYccRqTZiactc38zuNrM/Bve7mtn34h9N5Mit3rybu9+ZQ/8OTbnzrO5hxxGpVmIZ438xeNyw4P464LdxSyRylEpKI4wcMwMzePqy/tStHcv/5iLJI5Z/EVnu/ltgL4C77ybatkGkWnr040XMztvO7y7qQ/vmDcKOI1LtxFL4S8wslegHuphZZ6AkrqlEjtA/Fm3k/yat4MohHTn7uDZhxxGplmKZ1fMAMAFoZ2avAKcA18Y1lcgRWL99D3eMnU3PNo2599yeYccRqbYOWvjNrIO7r3b3CWaWAwwlOsTz81gbtplZU+B54Dii7xiuARYDbwKdgJXAJe6+9WhehEhpWYRb3phFcWmE54b3J7Vu7bAjiVRbFQ31jN93w93z3f2v7j7+MLt0PgVMcPceQF9gIdHeP5+6exbwaXBf5Kg8/Y9cpq3Ywm9+dByZGQ3DjiNSrVVU+I/qA9zg6t6TgRcA3L3E3bcB5wGvBA97BfjR0RxHZMqyAp75x1IuHNCOCwa0O/QTRJJcRWP8x5jZ0wf7YQz9ezKBfOAlM+sL5AC3AK3cfX2wj/Vmpitr5Ihs372XP3yey0tTVtI5PY0Hzjs27EgiNUJFhX8P0WJ9NPseAIx096lm9hSHMaxjZiOAEQAdOqg7hPxb0d4yXpmykuc+y6WwuJQL+rfj59/tTlq9WOYqiEhF/1I2u/srFfz8UPKAPHefGtx/m2jh32hmbYKz/TbAAT8zcPdRwCiAQYMG+VHkkARRFnHGzVzL4xMXs257Ead1z+Cus3vQs416BoocjooK/1HN1Xf3DWa2xsy6u/ti4AxgQfB1FfBw8P2vR3McSXzuzueL83lkwiIWbSikb7sm/P6SfpzYRX31RY7EQQu/uw+phP2PBEabWQrRJRyvJvqB8thgEffVwMWVcBxJULPWbOOhDxcydcUWOrVowHPDB3BO79aY6eJxkSMV10FRd58FDDrAj86I53Gl5ltRsItHP17Eh3M3kN4whf8971guG9xBfXdEKoE+DZNqJb+wmKc+XcIb09aQUqcWt56ZxX8Py6ShPrgVqTSH/NdkZn9x9ysPtU3kaOwsLmXUF8t5ftJySkojXD64AzefkUVGo3phRxNJOLGcRn1rcrSZ1QYGxieOJJuS0ghvfL2apz9dSsHOEs7t3YY7v9udzulpYUcTSVgV9eq5B/glUN/MduzbTHS2z6gqyCYJzN35YO56Hv14Mas272ZIZnOev6on/do3DTuaSMKraFbPQ8BDZvaQu99ThZkkwU1ZVsDDHy1iTt52erRuxEtXH8+p3TI0U0ekisSy5u49ZnYM0LH84939i3gGk8SzcP0OHv5oEf9ckk/bJqn8/uK+/Kj/MdSupYIvUpVi+XD3YeAyohdelQWbHVDhl5jkbd3N4xOXMG7WWhqn1uWX5/TgJyd2UutkkZDE8uHu+UB3dy+OdxhJLFt3lfCHz3N5ZcoqMBhxciY/O6UrTRrUDTuaSFKLpfAvB+oCKvwSk6K9Zbz05Ur+8Hkuu4pLuXBAO277TjfaNq0fdjQRIbbCvxuYZWafUq74x9CWWZJMWcR5JyePxz9ZwoYdRZzRoyV3nd2D7q0bhR1NRMqJpfD/LfgSOSB359OFm3hkwiKWbtpJv/ZNeeqyfpyQqSZqItVRLLN6XjGz+kCHoMumyDdyVm3lkY8WMW3lFjLT0/jjjwdw9nFqoiZSncUyq+cHwGNACtDZzPoBD7j7D+MdTqqvZfk7eXTCYibM30B6w3r85kfHcenx7dVETaQGiGWo53+AwcDnEO24aWad45hJqrFNO4p48tOlvPn1GlLr1OL273Tj2uzOWv1KpAaJ5V9rqbtv3++tu1bESjKFRXuDJmorKI1EuHJIR246vSvpDdVETaSmiaXwzzOz4UBtM8sCbgamxDeWVBclpRFGT13FM//IZcuuEn7Qty13ntWNji3URE2kpoql8I8E7iU6lXMM8DHwv/EMJeGLRJz35qzjsYmLWbNlD0O7tODu7/WgTzs1UROp6WKZ1bObaOG/N/5xpDqYvLSAhycsZN7aHfRs05hXrunNyVnpmqkjkiAqasv8pLvfambvcYAxfc3qSTzz1m7nkQmLmLS0gGOa1ueJS/tyXt9jqKUmaiIJpaIz/r8E3x+riiASnjVbdvP7iYsZP2sdTRvU5Vfn9uTKEztSr46aqIkkoor68ecE3/9pZilAD6Jn/ovdvaSK8kkcbdlVwrP/yOW1r1ZRqxb87NQu3HBqFxqnqomaSCKL5QKuc4E/AcuIrsDV2cyud/eP4h1O4mNPSRkvfrmCP32+jF0lpVwyqD23ntmN1k1Sw44mIlUgllk9vwdOc/dcADPrAnwAqPDXMKVlEd7KyeOJT5awqbCY7/RqxV3f7U5WKzVRE0kmsRT+TfuKfmA5sClOeSQO3J2JCzbyuwmLWJa/i4Edm/HcjwdwfKfmYUcTkRBUNKvnguDmfDP7EBhLdIz/YuDrWHZuZiuBQqIrd5W6+6Cg18+fgFSgFPiZu0874lcgFZq+cgsPfbSInFVb6ZKRxp+vHMhZvVppaqZIEqvojP8H5W5vBE4JbucDzQ7jGKe5e0G5+78Dfu3uH5nZOcH9Uw9jfxKD3E2FPDJhMZ8s2EjLRvV46ILeXDywHXXURE0k6VU0q+fqOB3TgcbB7SbAujgdJylt2lHE458sYez0NTRIqcOdZ3XjmuzONEhREzURiYplVk9nom0bOpV/fIwXcDkw0cwc+LO7jwJuBT42s8eAWsDQgxx3BDACoEOHDjEcSr5avpmfjZ5BYdFerhraiZGnZ9E8LSXsWCJSzcRyGjgeeAF4D4gc5v5Pcvd1ZtYS+MTMFgEXAbe5+ztmdkmw7zP3f2LwR2IUwKBBg9QNtALuzmtTV/Prv82nQ/MGjL1+CF1baqaOiBxYLIW/yN2fPpKdu/u64PsmMxtHtK//VcAtwUPeAp4/kn1LVElphPv/No8x09ZwWvcMnrysP03q6wIsETm4WAr/U2Z2PzCRby+2PqOiJ5lZGlDL3QuD22cBDxAd0z+F6MIupwNLjyy65BcW89PXcpi+ais/O7ULd5zVndrqqyMihxBL4e8NXEm0SO8b6vHgfkVaAeOCaYN1gNfdfYKZ7ST6x6QOUEQwji+HZ07eNq7/Sw5bd5fw9OX9+WHftmFHEpEaIpbCfz6Qebj9edx9OdD3ANsnAwMPZ1/ybeNnruUX78whvWE93r5hKMcd0yTsSCJSg8RS+GcDTdHVuqEriziPTFjEqC+WM7hTc/5wxQAtfSgihy2Wwt8KWGRmX/PtMX71469C23fvZeQbM/liST5XDOnAfd8/lpQ6uhhLRA5fLIX//rinkArlbirkv1+Zztpte/jt+b0ZfoKuaxCRIxfL0ov/LH/fzE4ChgP/PPAzpDL9fcFGbn1zFql1a/H6dUPUWE1EjlpM1/EHjdWGA5cAK4B34hlKohdlPfdZLr//ZAnHtm3MqCsH0bZp/bBjiUgCqKg7ZzfgMuByYDPwJmDufloVZUtau0tK+flbc/hg7nrO69eWhy/oQ/0ULYMoIpWjojP+RcAk4AflFmG5rUpSJbE1W3Zz3avTWbyxkHu+14MRJ2eqhbKIVKqKCv+FRM/4PzOzCcAbRJdelDiZsqyAG0fPoDTivPRfx3Nq95ZhRxKRBHTQ+YDuPs7dLyW6yPrnwG1AKzP7o5mdVUX5koK788qUlVz5wjSap6Xw1xtPUtEXkbiJZVbPLmA0MNrMmhNdgetuor175CgVl5Zx3/j5vDl9DWf0aMmTl/WjUaqarIlI/BzW6hzuvgX4c/AlR2nTjiJueC2HGau3cdNpXbn9O92opSZrIhJnWpYpJLPXRJusbd+zl+eGD+DcPm3CjiQiSUKFPwTvzsjj7nfnktGwHu/8dCi92jY+9JNERCqJCn8VKi2L8PBHi3h+8gqGZDbnueEDaKEmayJSxVT4q8i23SWMHDOTSUsLuOrEjvzq+72oW1tN1kSk6qnwV4ElGwu57tXprNu2h4cv6M1lg9VkTUTCo8IfZxPnb+C2N2dRP6UOb4wYwsCOarImIuFS4Y+TSMR55h+5PPH3JfRp14Q/XzmQNk3UZE1EwqfCHwe7iku5Y+xsJszfwPn9j+GhC3qTWldN1kSkelDhr2SrN+9mxF+ms2RjIb86tyfXZndWkzURqVZU+CvRl7kF3Pj6DNzhlWsGMywrI+xIIiL/QYW/Erg7L09ZyW8+WEhmehr/95NBdEpPCzuWiA9uwPYAAAqrSURBVMgBqfAfpeLSMn41bh5v5eTxnV6teOLSfjSsp1+riFRfqlBHYeOOIq7/Sw6z1mzj5jOyuPWMLDVZE5FqL66F38xWAoVAGVDq7oOC7SOBm4BS4AN3vyueOeJh5uqtXP+XHHYWl/LHHw/ge73VZE1EaoaqOOM/zd0L9t0xs9OA84A+7l5sZjVuxZG3c/L45btzadWkHq9eO5QerdVkTURqjjCGen4KPOzuxQDuvimEDEektCzCbz9cxItfrmBolxY8N3wAzdJSwo4lInJY4t0lzIGJZpZjZiOCbd2AYWY21cz+aWbHH+iJZjbCzKab2fT8/Pw4xzy0rbtKuOqlabz45QquPqkTr14zWEVfRGqkeJ/xn+Tu64LhnE/MbFFwzGbAEOB4YKyZZbq7l3+iu48CRgEMGjTICdHiDdEmaxu2F/G7i/pwyaD2YcYRETkqcS387r4u+L7JzMYBg4E84N2g0E8zswiQDoR/Wn8AE+at5/axs2lYrw5vXD+EAR2ahR1JROSoxG2ox8zSzKzRvtvAWcA8YDxwerC9G5ACFBxsP2GJRJwnPlnCDa/NIKtVI94bma2iLyIJIZ5n/K2AcUGfmjrA6+4+wcxSgBfNbB5QAly1/zBP2HYWl3L7m7OYuGAjFw5ox4PnH6cmayKSMOJW+N19OdD3ANtLgCviddyjtWrzLq57dTrL8ndx3/d7cfVJndRkTUQSiq7cLWfy0miTNYBXrh5MdlZ6yIlERCqfCj/RJmsvTF7Bbz9cSFbLRoz6yUA6tlCTNRFJTElf+Iv2lvHLcXN5d8ZavntsKx6/pB9parImIgksqSvchu1FXP9aDrPXbOO2M7sx8vSuarImIgkvaQt/zqqt3PBaDruLS/nzlQP57rGtw44kIlIlkrLwj52+hl+Nm0frJqm8du0JdG/dKOxIIiJVJqkK/96yCA9+sJCXp6wku2s6zw7vT9MG6rcjIsklaQr/ll0l3Dh6Bv9avplrsztzz/d6UKd2vHvUiYhUP0lR+Beu38F1r05nU2Exv7+4LxcObBd2JBGR0CR84f9w7nruGDubxvXrMPb6E+nXvmnYkUREQpXQhf+5z3J59OPF9O/QlD9fMZCWjVPDjiQiErqELvydWqRx6aD2PPCjY6lXR03WREQgwQv/uX3acG4fLYIuIlKeprWIiCQZFX4RkSSjwi8ikmRU+EVEkowKv4hIklHhFxFJMir8IiJJRoVfRCTJmLuHneGQzCwfWHWET08HCioxTk2g15wc9JqTw9G85o7unrH/xhpR+I+GmU1390Fh56hKes3JQa85OcTjNWuoR0Qkyajwi4gkmWQo/KPCDhACvebkoNecHCr9NSf8GL+IiHxbMpzxi4hIOSr8IiJJJikKv5n9r5nNMbNZZjbRzNqGnSnezOxRM1sUvO5xZpbwiw2b2cVmNt/MImaWsFP+zOxsM1tsZrlmdnfYeaqCmb1oZpvMbF7YWaqCmbU3s8/MbGHw//Qtlbn/pCj8wKPu3sfd+wHvA/eFHagKfAIc5+59gCXAPSHnqQrzgAuAL8IOEi9mVht4Dvge0Au43Mx6hZuqSrwMnB12iCpUCtzh7j2BIcCNlfnfOSkKv7vvKHc3DUj4T7TdfaK7lwZ3vwLahZmnKrj7QndfHHaOOBsM5Lr7cncvAd4Azgs5U9y5+xfAlrBzVBV3X+/uM4LbhcBC4JjK2n9Cr7lbnpk9CPwE2A6cFnKcqnYN8GbYIaRSHAOsKXc/DzghpCxSBcysE9AfmFpZ+0yYwm9mfwdaH+BH97r7X939XuBeM7sHuAm4v0oDxsGhXnPwmHuJvm0cXZXZ4iWW15zg7ADbEv4dbLIys4bAO8Ct+41cHJWEKfzufmaMD30d+IAEKPyHes1mdhXwfeAMT5ALNg7jv3OiygPal7vfDlgXUhaJIzOrS7Toj3b3dytz30kxxm9mWeXu/hBYFFaWqmJmZwO/AH7o7rvDziOV5msgy8w6m1kKcBnwt5AzSSUzMwNeABa6++OVvv8EORGskJm9A3QHIkTbO9/g7mvDTRVfZpYL1AM2B5u+cvcbQowUd2Z2PvAMkAFsA2a5+3fDTVX5zOwc4EmgNvCiuz8YcqS4M7MxwKlEWxRvBO539xdCDRVHZpYNTALmEq1bAL909w8rZf/JUPhFROTfkmKoR0RE/k2FX0Qkyajwi4gkGRV+EZEko8IvIpJkVPil0sXSPdLMXjazi6o625Ews5uC1+Jmll5uew8z+5eZFZvZnXE8fhszez+4/V9m9my8jnWAY2eY2YSqOp5UDRV+qVRhd48Mjl/ZvgTOJHoNSHlbgJuBx+JwzPJuB/4vngcwswNexe/u+cB6MzspnseXqqXCL5XtcLpHnmxmU8xs+b6zf4t61MzmmdlcM7s02H7qvrPe4P6zZvZfwe2VZnafmU0GLjazm81sQbAWwRtH+4Lcfaa7rzzA9k3u/jWwd/+fmdl4M8sJeqmPKLd9p5k9Evzs72Y22Mw+D34HPzxIhAuB8mfdbc1sgpktNbPfldv35cHvbJ6ZPVL+mOVuX2RmLwe3Xzazx83sM+ARMzslWLNilpnNNLNGwdPGAz+O4VclNUTC9OqRauNwuke2AbKBHkTbDrxNtJ9+P6Av0as0vzazWPrrF7l7NoCZrQM6u3vxgRagMbPuHLxb6anuvi2G4x3KNe6+xczqE30N77j7ZqJtwT9391+Y2TjgN8B3iL47eoX92i+YWWdgq7sXl9vcj2i3xmJgsZk9A5QBjwADga3ARDP7kbuPP0TObsCZ7l5mZu8BN7r7l0FzsKLgMdODnJIgVPilsh1O98jx7h4BFphZq2BbNjDG3cuAjWb2T+B44FCdCcsX8jnAaDMbT/Rs9dthoj37+x1if0fr5qCFBESbqmURbZ9Rwr/P3ucCxe6+18zmAp0OsJ82QP5+2z519+0AZrYA6Ai0IPoHJT/YPho4mQO8/v28FfyuITqk9Xjw3HfdPS/YvglI+FXrkomGeqSyHU73yPJnsbbf9/2V8u3/X1P3+/mucrfPJfo5w0AgZ//xazPrXm5IY/+vo16i0sxOJfqZwInu3heYWS7v3nKdUiMEv4PgD+CBTsT28J+vtfzvrSx43sF+b/DtP7wH/b25+8PAfwP1ga/MrEe55+ypYP9Sw6jwS2U72u6RXwCXmlltM8sgetY6jegHq73MrJ6ZNQHOONCTzawW0N7dPwPuApoCDcs/xt0Xu3u/g3xVxjBPE6LDM7uD4jnkKPa1hAO/E9jfVOAUM0sPPuC+HPhn8LONZtYz+N2cf7AdmFkXd5/r7o8QHd7ZV/i7EV3WUhKEhnqkUrl7qZndBHzMv7tHzj+MXYwDTgRmEz1TvcvdNwCY2ViiwzhLiZ5FH0ht4LXgj4MBTxxtMTezm4n+EWkNzDGzD939v82sNdEC2RiImNmtRMfqJwA3mNkcYDHRpS+PiLvvMrNlZtbV3XMreNx6iy4y9BnR1/1huYVp7ia61vQaogW84YH3wq1mdhrRdxELgI+C7acRXcNCEoS6c4pUc8FnBQPd/VchHf8L4Dx33xrG8aXy6YxfpJpz93Fm1iKMYwfDbY+r6CcWnfGLiCQZfbgrIpJkVPhFRJKMCr+ISJJR4RcRSTIq/CIiSeb/A/viQVobeU1jAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"time = np.array([-3, -2, -1, 0, 1, 2])\n",
"T_a = np.array([55, 58, 60, 65, 66, 67])\n",
"\n",
"plt.plot(time, T_a)\n",
"plt.xlabel('0 hours = 11am (hours)')\n",
"plt.ylabel(\"Ambient Temperature (F)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"2a. No, this does not look correct. For a better answer get more values for T_a so that it is a continuous function."
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f53d55e7048>"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUVfrA8e87k95IQkIIJDSliDSRqlioUhRBiqDSFHFddbGh2Nayu66r2F07rIoIUgT9odJBBBEBBQRBOhJCCem9nt8fM8QQkjBJZjIp7+d55iFz58657w0w79x7znmPGGNQSimlACzuDkAppVT1oUlBKaVUIU0KSimlCmlSUEopVUiTglJKqUIe7g6gMsLCwkyzZs3cHYZSStUo27ZtO2OMCS/ptRqdFJo1a8bWrVvdHYZSStUoInK0tNf09pFSSqlCmhSUUkoVcllSEJFZInJaRHYV2faSiOwVkZ0islhEgou89piIHBCR30XkOlfFpZRSqnSu7FP4CHgL+KTItpXAY8aYPBH5D/AY8KiItAXGAJcCjYBVItLKGJPvwviUUlUsNzeXmJgYsrKy3B1KneDj40NUVBSenp4Ov8dlScEYs15EmhXbtqLI0x+BkfafbwTmGWOygcMicgDoBmxyVXxKqaoXExNDYGAgzZo1Q0TcHU6tZowhPj6emJgYmjdv7vD73NmncDvwrf3nxsCxIq/F2LcppWqRrKws6tevrwmhCogI9evXL/dVmVuSgog8AeQBc85uKmG3Esu3isgUEdkqIlvj4uJcFaJSykU0IVSdivyuqzwpiMgE4HrgVvNn3e4YILrIblFAbEnvN8a8b4zpYozpEh5e4tyLC0rKSmLqt1NJzkqu0PuVUqq2qtKkICIDgUeBocaYjCIvfQWMERFvEWkOtAR+clUc++P3898t/+Xeb+911SGUUtWUiPDQQw8VPp8xYwbPPPOMU9qeOHEiCxcurFQbV1xxRYXe98wzzzBjxoxKHRtcOyR1LraO4tYiEiMid2AbjRQIrBSR7SLyLoAxZjcwH/gNWAbc48qRR10bd+Wpq5/i052fMn/3fFcdRilVDXl7e/PFF19w5swZd4dyjvx820feDz/84NY4XJYUjDFjjTGRxhhPY0yUMWamMeZiY0y0MaaT/fGXIvv/yxhzkTGmtTHm27LadoYnrn6C7o27c9fSu4hJiXH14ZRS1YSHhwdTpkzh1VdfPe+14t/0AwICAFi3bh3XXHMNo0ePplWrVkyfPp05c+bQrVs32rdvz8GDBwvfs2rVKq666ipatWrF0qVLAdsH/rRp0+jatSsdOnTgvffeK2y3d+/e3HLLLbRv3/6cYwK8+OKLtG/fno4dOzJ9+nQAPvjgA7p27UrHjh0ZMWIEGRlFb7o44ffj1NZqEA+LB7OHz+ay9y5jwpIJrBy3EovoBG+lqsr9y+5n+8ntTm2zU8NOvDbwtQvud88999ChQwceeeQRh9vesWMHe/bsITQ0lBYtWjB58mR++uknXn/9dd58801ee8123CNHjvDdd99x8OBBevfuzYEDB/jkk0+oV68eW7ZsITs7myuvvJIBAwYA8NNPP7Fr167zho1+++23LFmyhM2bN+Pn50dCQgIAN910E3feeScATz75JDNnzuS+++5z+DwupE5/Cras35LXBr7GmsNreGXTK+4ORylVRYKCghg/fjxvvPGGw+/p2rUrkZGReHt7c9FFFxV+qLdv354jR44U7jd69GgsFgstW7akRYsW7N27lxUrVvDJJ5/QqVMnunfvTnx8PPv37wegW7duJc4jWLVqFZMmTcLPzw+A0NBQAHbt2sVVV11F+/btmTNnDrt3767or6FEdfZK4aw7LruDb/Z/w+OrH6dfi350atjJ3SEpVSc48o3ele6//346d+7MpEmTCrd5eHhQUFAA2CZ/5eTkFL7m7e1d+LPFYil8brFYyMvLK3yt+DBQEcEYw5tvvsl1151bwWfdunX4+/uXGJ8xpsQhpRMnTmTJkiV07NiRjz76iHXr1jl4xo6p01cKYPsL++CGDwj3D2fsorFk5Dr3/pxSqnoKDQ1l9OjRzJw5s3Bbs2bN2LZtGwBffvklubm55W53wYIFFBQUcPDgQQ4dOkTr1q257rrreOeddwrb27dvH+np6WW2M2DAAGbNmlXYZ3D29lFqaiqRkZHk5uYyZ86cspqokDqfFADq+9Xn42Efs/fMXh5a/tCF36CUqhUeeuihc0Yh3XnnnXz33Xd069aNzZs3l/otviytW7fmmmuuYdCgQbz77rv4+PgwefJk2rZtS+fOnWnXrh133XXXOVcXJRk4cCBDhw6lS5cudOrUqXC46T/+8Q+6d+9O//79adOmTbnjuxD5c/5YzdOlSxfjzEV2pq2YxoxNM1h882KGtRnmtHaVUjZ79uzhkksucXcYdUpJv3MR2WaM6VLS/nqlUMS/+v6LzpGdueOrO3SYqlKqTtKkUISX1Yu5I+aSnZfNuMXjyC/Qyt1KqbpFk0Ixreq34q3Bb7HuyDpe2PCCu8NRSqkqpUmhBBM6TmBsu7E8ve5pNv6x0d3hKKVUldGkUAIR4d3r36VpcFNu+eIWEjMT3R2SUkpVCU0KpQjyDmLeiHnEpsZyx1d3UJNHaSmllKM0KZSha+OuvND3BRbvXczbW952dzhKKSfQ0tll06RwAQ/0fIAhLYfw4IoH+eXEL+4ORylVSVo6u2yaFC7AIhY+GvYR4X7hjF44mpTsFHeHpJSqBC2dfYHfj1Nbq6XC/MKYO2IuvT/uzV1L7+Kzmz7TdWaVqiQtna2ls2u0q5pexXO9n2Pernm8v+19d4ejlKoELZ1dOr1SKIfpvaaz/uh6pi6bSrfG3bgs8jJ3h6RUjaWls7V0do1nEQuzh88mzC+M0QtHk5yV7O6QlFIVpKWzS6ZJoZzC/cP5fOTnHE48rPMXlKrhtHT2+bR0dgXN+GEG01ZO47XrXmNqj6luiUGpmkZLZ1c9LZ1dRR7q+RDD2gzj4ZUP88Mx944rVkopZ9GkUEEiwv9u/B9N6jVh9ILRxKXHuTskpZSqNE0KlRDsE8yi0YuIz4xn7KKxuv6CUqrG06RQSZ0aduLtwW+z+vBqnlr7lLvDUUqpStGk4ASTLpvE5Msm8+8N/+bLvV+6OxyllKowTQpO8ubgN7k88nLGLxnPvvh97g5HKaUqRJOCk/h4+LBo9CI8LZ7c9PlNpOWkuTskpVQpFi9ejIiwd+/eCrfhSJns559//pzn7i6L7QhNCk7UNLgpc0fMZc+ZPTqxTSknWJWWwJiY3+hzdAdjYn5jVVqCU9qdO3cuvXr1Yt68eU5przTFk4K7y2I7QpOCk/W/qD/P93me+bvn8/Kml90djlI11qq0BGYkxHAqPxcDnMrPZUZCTKUTQ1paGhs3bmTmzJmFSWHdunVce+21jBw5kjZt2nDrrbcWfql77rnn6Nq1K+3atWPKlCnnfdlbvXo1w4cPL3y+cuVKbrrpJqZPn05mZiadOnXi1ltvBdxfFtsRWhDPBR658hG2xG7h0VWP0qlhJ/q16OfukJSqdt5KOM6BnMxSX/8tOwNbOvhTtjG8GB/D0lISw8Vevtwb2rjM4y5ZsoSBAwfSqlUrQkND+fnnnwH45Zdf2L17N40aNeLKK69k48aN9OrVi3vvvZe///3vAIwbN46lS5dyww03FLbXp08f7rnnHuLi4ggPD+d///sfkyZN4oYbbuCtt95i+/bzy4O7qyy2I/RKwQXOTmxrE9aGMQvHcCTpiLtDUqrGKZ4QLrTdUXPnzmXMmDEAjBkzhrlz5wK2EtZRUVFYLBY6depUWA577dq1dO/enfbt27NmzZrzSlWLCOPGjePTTz8lKSmJTZs2MWjQoDJjcFdZbEfUySuFVWkJfJh0ktP5uTSwejI5uCH9AkKdeoxA70CW3LyErh90Zdi8Yfxwxw/4efo59RhK1WQX+kY/JuY3TuWfX6U0wurJaw0vrtAx4+PjWbNmDbt27UJEyM/PR0QYPHjwOaWxrVYreXl5ZGVl8de//pWtW7cSHR3NM888Q1ZW1nntnr0y8PHxYdSoUXh4lP3R6q6y2I6oc1cKrrpPWZKW9Vsyd8Rcdp7aye1f3q4dz0qVw+TghngX++D0FmFycMMKt7lw4ULGjx/P0aNHOXLkCMeOHaN58+Zs2LChxP3PJoCwsDDS0tJKHW3UqFEjGjVqxD//+U8mTpxYuN3T07PE8tvuKovtiDqXFD5MOkm2Of8+5YdJJ11yvEEtB/F83+f5fPfnvLjxRZccQ6naqF9AKA+HRhFh9USwXSE8HBpVqav6uXPnntMpDDBixAg+++yzEvcPDg7mzjvvpH379gwbNoyuXbuW2vatt95KdHQ0bdu2Ldw2ZcoUOnToUNjRfJa7ymI7wmWls0VkFnA9cNoY086+bRTwDHAJ0M0Ys7XI/o8BdwD5wN+MMcsvdIyKlM7uc3RHiXckBVjTtGO52nKUMYaxi8Yyf/d8lt6ylMEtB7vkOEpVd7W5dPa9997LZZddxh133OHuUM5RnUpnfwQMLLZtF3ATsL7oRhFpC4wBLrW/520RsboiqAZWz3JtdwYRYebQmXRs2JGxi8ay90zFJ8wopaqfyy+/nJ07d3Lbbbe5O5RKc1lSMMasBxKKbdtjjPm9hN1vBOYZY7KNMYeBA0A3V8RV0n1KAW4NauCKwxXy9/LnyzFf4m31ZujcoSRmJrr0eEqpqrNt2zbWr19/Tmd1TVVd+hQaA8eKPI+xbzuPiEwRka0isjUurvxrGBS/TxlssSLA/6UnkJpf9vJ4ldWkXhO+uPkLjiQdYfTC0eQVuPZ4SlVHOuCi6lTkd11dksL5Y7MoeTCyMeZ9Y0wXY0yX8PDwCh2sX0Ao86LasqZpRxZHt+NfDZpzJCeLaacPkebiNRF6NenFu9e/y6pDq3hw+YMuPZZS1Y2Pjw/x8fGaGKqAMYb4+Hh8fHzK9b7qMk8hBogu8jwKiK2qg/fwDeLZ8Gb8Pe4I004d4qWIFgRYXNKlAcDtl93OrtO7ePXHV7k0/FLu6nKXy46lVHUSFRVFTEwMFbnKV+Xn4+NDVFRUud5TXZLCV8BnIvIK0AhoCfxUlQH09Avi6fCmPBN3hEdPHeLFiBb4uzAxvNT/Jfae2cu9395Ly/ot6dO8j8uOpVR14enpSfPmzd0dhiqDy24fichcYBPQWkRiROQOERkuIjFAT+BrEVkOYIzZDcwHfgOWAfcYY6p8bctefvV4OrwZv+dk8OjpQ2S48FaS1WJl3sh5tKrfipHzR+oaDEqpasFl8xSqQkXmKTjiu/QknjtzlEu9/flPg+b4uvCK4VDiIbp/2J0QnxB+nPwjob7OLbehlFLFuWueQo11jX8wT4U1ZXd2OtNPHybThVcMLUJasPjmxRxNPsqI+SPIyc9x2bGUUupCNCmU4lr/YJ4Ia8Ku7HQeP32YrIIClx2rV5NezBw6k3VH1nH30rt1ZIZSym00KZShj38Ij4U1YWd2Oo/HuTYx3NbhNp66+ilmbZ/FCxtecNlxlFKqLJoULqCffwjT6zdhe1YaT8YdJtuFieHZa59lbLuxPL7mcebvnu+y4yilVGk0KTigf0AIj9aP5md7YsgxrkkMIsKsG2dxZfSVjF88nh+OVf/1XJVStYsmBQddFxDKtPrRbMtK46nTR1yWGHw8fFgyZgnR9aK5cd6NHEg44JLjKKVUSTQplMOggFAeCo3ip6xU/h7nusQQ5hfGN7d8gzGGwXMGcybjjEuOo5RSxWlSKKchgfV5MDSKzZmpPBt3lFwXJYaW9Vvy1div+CP5D4bOHUpmbukLnCullLNoUqiAGwLr80BoY37ITHFpYrgi+grm3DSHH2N+5LbFt5Hv4mJ9SimlSaGChgaG8bfQxmzMTOEfcX+Q56K5BSPajuCV617hiz1f8ODyB3UOg1LKpapLQbwaaXhgGMYY3kyM5d4T+0ksyCMuP5cGVk8mBzes1FqyRd3f436OJh3ltc2vEV0vmoeveNgp7SqlVHGaFCrppqBw9mZnsDIjqXDbqfxcZiTEADgtMbx83cvEpsUybeU0GgU24pb2tzilXaWUKkpvHznBzuz087ZlG8OHSSeddgyLWPh42Mdc0/QaJi6ZyKpDq5zWtlJKnaVJwQlO5+eWa3tFnZ3D0CasDcM/H8622G1ObV8ppcpMCiLiJSLDRORlEZkrIrNE5EERaVNVAdYEDayeJW4Pszr/7lywTzDLbltGfd/6DP5ssE5uU0o5ValJQUSeBDYDvYEdwMfYVkjzAF4VkWUi0q5KoqzmJgc3xFvOX2Y6t8AQm5vt9OM1CmzE8tuWk1+Qz3WfXseJ1BNOP4ZSqm4q60rhV2PMZcaYqcaYT4wxy4wxS4wxLxpjBgGTAL8qirNa6xcQysOhUURYPREgwurJhKAGFAjce/IA+7IznH7M1mGt+ebWbziVdopBcwaRlJV04TcppdQFlLrymoiIqeaD4l218pqz/JGbxSOnDpFSkM9z4c3o4hvo9GOsOLiC6z+7nh5RPVh+23J8PX2dfgylVO1S0ZXXCnsxReQ1p0dVBzTx9OGthi2J9PBi+ulDrExLdPoxBlw0gNnDZ7Phjw2MXjiaXCd3biul6paykkLRm+RXuzqQ2irMw5PXG15MO29/no//g/kpcU4/xs3tbua/g//L0n1Luf2r2ylwUdkNpVTtV9bwmGp966gmCbBYeTGiBc+f+YN3EmOJz8/lruBILCV0TlfU3V3vJiEzgSfXPkmwdzBvDHoDcWL7Sqm6oayk0EZEfsZ2xdDa/jP258YY09nl0dUiXmLhqbCmhCYeZ35KHPF5uTwaFo2nOG+qyONXPU5CZgKv/PgKIb4hPNf7Oae1rZSqG8pKCu2rLIo6wirCfSGNCbN68kHSSZJO5/FceDP8LFantC8izBgwg+TsZP6x/h8EeQdpnSSlVLmUmhSMMQerMpC6QkS4pV4EoVZPXoo/xv2nDvJCg+aEljIBriLtv3f9e6TmpDJt5TQCvQK5q8tdTmlbKVX7lTV5ba2I3C0ijYpt9xCRq0VkpohMcn2ItdPAgFCeb9CcY7nZ3HfyADFOnORmtViZPXw2Q1oO4e6v7+bTnZ86rW2lVO1W1g3tIYAnsFhEYkRkp4jsAw5hm7j2jjHmf1URZG3V3TeIVyIuIr0gn/tOHmCvEye5eVm9WDBqAdc2u5aJSybyxZ4vnNa2Uqr2KnXy2jk7iXgDDYBMY0y1WTC4uk9ec9Sx3GweOXWIpII8ng1vSjffIKe1nZaTxoDZA9gau5XFNy9mSKshTmtbKVUzVXTyWiFjTLYx5lh1Sgi1SbSnN/+NvJgoDy8eP32YFWkJTms7wCuAb279hg4RHRgxfwQrD650WttKqdpHS2dXE6FWT15reDEdfQL4d/wx5iafdtrSm8E+wSy/bTmt6rfixnk38t2R75zSrlKq9tGkUI34W6z8u0Fz+vgF837SCf6bGEuBkxJDfb/6rBq/imbBzRjy2RA2/rHRKe0qpWoXR/sUooCWxpi19v4FD2PM+cuNVbHa0qdQXIExvJMYy8LUM1zi6Uu8E9d+PpF6gms/vpbY1FhW3LaCntE9nRi5UqomqFSfgojcjm0dhQ/tm5oCXzovPFWcRYS/hjSir2899uRmcjo/F8Ofaz+vqkSfQ2RgJGvGryHCP4KBcwby0/GfnBe4UqrGc+T20d+AHkAKgDFmH7aRSMqFRIRdOecPUXXG2s+NgxqzdsJawvzCGDB7AFuOb6lUe0qp2sORpJBljMk5+0RErJxbQVW5iCvXfo6uF83aCWsJ9Q2l/+z+bI2tfbfhlFLl50hS2CgijwA+ItIb+BxYeqE32ddzPi0iu4psCxWRlSKy3/5nSJHXHhORAyLyu4hcV5GTqW1KW/vZXyxOGZnUpF6TwsTQ75N+esWglHIoKTwCpAJ7ganAauAJB973ETCw2LbpwGpjTEt7O9MBRKQtMAa41P6et+1XJHVaSWs/W4A0U8CzZ46SWZBf6WM0DW7KuonrCq8YNsdsrnSbSqmaq8ykYP9gnmWMeccYM9wYM8z+8wVXcTHGrAeK94jeCHxs//ljYFiR7fPsk+QOAweAbuU5kdqopLWfp4dGc3dIJN9nJPPXkwc47oSaSU3qNeG7id/Z+hg+HcCmY5sqH7xSqkYqMykYY/KBSBFxTglPiDDGnLC3fYI/O6wbA8eK7Bdj33YeEZkiIltFZGtcnPNXMatu+gWEMi+qLWuadmReVFv6B4YyOqgBLzZoQUJ+Ln85uZ/NmSmVPk50vWjWTVxHhH8EAz4dwPdHv3dC9EqpmsaR20eHgO/t9/z/dvbh5DhK6rgu8aa5MeZ9Y0wXY0yX8PBwJ4dRc1zuG8h7ka1oaPXksdOH+TT5VKX7GaKColg3cR1RQVEMnDOQtYfXOilapVRN4UhSiANWAn5AeJFHRZwSkUgA+5+n7dtjgOgi+0UBsRU8Rp3R0MOLNxu2pI9fMDOTTvL0maNkVLKfoVFgI9ZNWEfz4OYM/mwwyw8sd1K0SqmawKEZzRVuXKQZsNQY087+/CUg3hjzgohMB0KNMY+IyKXAZ9j6ERph64Ruab99VaraOqO5vIwxLEw9w7uJsUR7evPP8OZEeXpXqs0zGWfoP7s/v8X9xoJRCxjaeqiTolVKuVtlZzSvFJEVxR8OvG8usAnb+s4xInIH8ALQX0T2A/3tzzHG7AbmA78By4B7LpQQ1J9EhFFB4bzUoAWJ+Xn85cQ+NmVUrp8hzC+MNePX0KlhJ0bMH8Hnuz53UrRKqersglcKItK9yFMfYASQbYyZ5srAHKFXCuc7mZfD3+OOcCAnk0n1GnJrvQZYpOJzDVOyU7j+s+vZeGwjH97wIZMu08X2lKrpyrpSKHWN5rOMMcUHrn8nIlp7uZpq6OHFmxEX83JCDLOST7I/J5PpYdH4WSo27SPIO4hlty1j+OfDuf2r20nNSeVv3Z09zkApVV04cvsoqMgjWET6ApFVEJuqIG+LhcfqR3NPSCM2Zibz15P7+SM3q8Lt+Xn68dWYrxjeZjhTl03ln+v/6bS1HpRS1Ysjo492A7vsf/6CbTbzna4MSlWeiDAyKJyXIy4iKT+Pv57YX6l+Bm8Pb+aPms+4DuN4au1TTFs5TRODUrXQBW8fAS2MMedUYBMRR96nqoFOPgG8F9mKp+KO8HjcYSbWi2BcvYgK9TN4WDz4aNhH1POux8ubXiYxM5H3bngPD4v+c1CqtnDkf/NmoHOxbT+VsE1VUxH2foZXEmL4KPkU+3My6eEbxKfJpzhdzsV7LGLhjUFvEOobynPrnyMpO4k5N83Bx8OnCs5EKeVqpSYFEWmAre/AV0Ta8+es4yBsE9lUDeJtsTC9fjStvXx5MzGWHzJTCqeMn128B3AoMYgIz/Z+lhDfEB5Y/gBDPhvCkpuXEOgd6MIzUEpVhbL6FIYAb2GbXfw28F/743HgKdeHppxNRLgpKJwQi8d5NUQqsnjP/T3u55Nhn/Ddke/o/XFvTqefvvCblFLVWqlXCsaY/wH/E5HRxpj5VRiTcrGkgrwSt1dk8Z5xHccR6hvKqAWj6DWrF8tvW07zkOaVDVEp5SYXHH1kjJkvIteJyIMi8vjZR1UEp1yjtMV7QirYYTyk1RBWjV/FmYwzXDHrCnac3FGZ8JRSbuTIPIW3gQnAg4AvcBtwsYvjUi5U0uI9AAkFecxOOkV+BYaaXhF9BRtu34CHxYOrP7paK6wqVUM5Mk+hlzHmFmyF7J4CumPrZ1A1VEmL9zwUGkVfv2BmJZ/kgVMHOZmXc8F2imsb3pYfbv+hsPS21ktSquZx5H7B2amwWSLSEIgHmrksIlUl+gWEnjfS6PrA+nRPC+L1hBgmx/7O1NAo+geElNJCyaLrRbNh0gaGfT6MMYvGEJMSw4M9H0QqUX9JKVV1HLlS+EZEgoEZwHbgCLDQlUEp9+kfEMKHjVrT3MuH5+P/4J9xR0kr5xoNIb4hLL9tOSPbjuThlQ/zwPIHyHfCetJKKde70BrNFuBbY0ySMWYB0Bxob4zRjuZarKGHF69FXMzt9RqyNiOJybG/szMrrVxt+Hj48PnIz7m/+/28vvl1Ri8cTWZuposiVko5y4XWaC4AXi/yPNMYk+DyqJTbWUUYFxzBmw0vxirCA6cOMjPxBHnl6IS2iIVXB77Kq9e9yuI9i+n7SV/i0mv/utpK1WSO3D5aKSI3ujwSVS219fbng8hWXOcfyqcpp7nv5H5icrPL1cb9Pe5n4eiF/HLyF3rO7Mm++H0uilYpVVmOJIV7gcUikikiCSKSKCJ6tVCH+FmsPBIWzTNhTTmel8OdJ/bxdWp8uaqk3nTJTaydsJaU7BR6zuzJ+qPrXRixUqqiHEkKYYAnEACE25+HuzIoVT1d4x/MzMhWXOLtx4yEGJ4+c5Tk/JJnR5ekR1QPfpz8Iw38G9Dvk37M3jHbhdEqpSrCkRnN+cAo4FH7z5FAJ1cHpqqncA8vZjRowV+CI9mUkcLkE7+zLTPV4fe3CGnBD7f/QK8mvRi/ZDxPrXmKAlPgwoiVUuXhyIzmt4DewDj7pgzgXVcGpao3iwg312vA25EX4ydWHj59iHcSY8lx8MM9xDeEZbct447L7uCf3/+TsYvG6sgkpaoJRyavXWGM6SwivwAYYxJExMvFcakaoKWXH+9FtuLdxFjmp8SxLTOVvn7BfJkWf8F1GrysXnxwwwe0rt+aR1c9yuHEw3w55ksiA3WlV6XcyZE+hVz7fAUDICL1Ab3eVwD4WCzcXz+K58ObE5ubzfvJJzmVn4vhz3UaVqWVPC5BRJh25TQW37yY3+J+o+sHXfn5xM9VewJKqXM4khT+CywCwkXkWWAD8B+XRqVqnJ5+Qfhbredtd2Sdhhvb3MjG2zdiEQu9ZvViwe4FrgpTKXUBjnQ0fwI8ia3MRQIwyhgzz9WBqZonvpSRSI6s09CxYUe23LmFTg07MXrhaJ5e+7R2QCvlBo5cKQBYgVwgpxzvUXVMaes0WMChMhkRARGsnbCWiZ0m8tz65xi1YBRpOeUrr6GUqhxHRh89AcwFGmErmf2ZiDzm6vbJlKMAABu0SURBVMBUzVPSOg2eCP5iYeqpg7wcf4zUC8xr8PbwZtbQWbwy4BWW7F1Cz5k9OZR4yJVhK6WKkAvNShWRPcDlxpgM+3M/YJsx5pIqiK9MXbp0MVu3bnV3GKqIVWkJfJh08pzRR1f61eOj5FMsTImjnsWDe0Mb0dsv+ILltFceXMnNC28G4PORn9P/ov5VcQpK1Xoiss0Y06XE1xxICsuA0caYFPvzIGCuMWaI0yMtJ00KNcu+7AxeTohhX04m3X0DuT80ioYeZY9uPphwkGGfD+O3uN/4T7//8FDPh3RtBqUqqayk4Ej/QAawW0Q+FJEPgF+BJBF5RURecWagqnZr5e3H2w1bck9II3ZkpTMp9nfmp8SVufznRaEXsemOTQxvM5xpK6cxdtFY0nPSqzBqpeoWR64U7ijrdWPMTKdGVA56pVBzncrL4bWEGH7MTKWlly8Ph0bRytuv1P2NMbyw4QWeWPME7Rq0Y/HNi7ko9KIqjFip2qNSt4+qM00KNZsxhu8yknkz8ThJ+XmMCAxnUnAEvpbz5zucteLgCsYuGkuBKeDT4Z8ypJXb72IqVeNU6vaRiAwUkS0iclpLZytnEhGu9Q/m40ZtGBJQnwWpcUyK/Z1NGSmlvmfARQPYNmUbzYObc/3c63l67dO61KdSTuRIn8JbwF1AY7R0tnKBAIuVB+tH8UbERfhYLDwed5hn446QUMqkt2bBzdh4+0YmdZrEc+ufY/BngzmTcaaKo1aqdnKkT2Ed0Me+NGe1orePap8cU8C85Dg+TT6Ft8XClOBIvBFmJZ88r8ieMYYPf/6Qe7+9lwj/CBaMWkD3qO7uPgWlqr3KDkntBjwNrAMK12E0xrxRiYCmAncCAnxgjHlNREKBz4FmwBFsw2ATy2pHk0LtdSw3m1fij7E9Ox3BXo3RzluEh0OjCquvbovdxsgFIzmecpwZA2ZwX7f7dNiqUmWo7JDUZ4F8IBjbbaOzj4oG0w5bQugGdASuF5GWwHRgtTGmJbDa/lzVUdGe3rwScRFBFivFv7YUL7J3eaPL+XnKzwxqOYipy6YyeuFokrOSqzZgpWoJR9ZTaGCMudyJx7wE+LHIDOnvgOHAjcC19n0+xnZl8qgTj6tqGBEhtZRO5OJF9kJ8Q1hy8xJm/DCDx1Y/xi8nfmHBqAVcFnlZVYSqVK3hyJXCahHp48Rj7gKuFpH69pIZg4FoIMIYcwLA/mcDJx5T1VClFdmzIvxSrMje2fUZ1k1cR1ZeFj1m9uDtLW9Tk4ddK1XVHEkKdwKrRCTNGUNSjTF7sK3HsBJYBuwAHF79XUSmiMhWEdkaFxdX0TBUDVFakT0/ER48dZAnTh/mj9ysc17v1aQX2/+ynb7N+3LPN/cwasEokrKSqjJspWosRzqaS5xJZIxxyuBwEXkeiAGmAtcaY06ISCSwzhjTuqz3akdz3VBSkb2r/IJZlBrHnOTTZJsChgaGMaFeBPWsf94RLTAFvLLpFR5b/RiNAxszb+Q8ekT1cOOZKFU9VHpGs4iMAVoYY54XkShst3q2VSKgBsaY0yLSBFgB9AQeB+KNMS+IyHQg1BjzSFntaFJQifm5fJR0iqVp8fhZrIyr14BhgWF4yZ8XwZtjNjNm0RhiUmL4R+9/8MiVj2ARXRZE1V2VHZL6FuAJXG2MucQ+dHS5MaZrJQL6HqiPbeGeB40xq+1rP88HmgB/YFvhrczbVJoU1FmHc7J4LzGWzVmpNPLwYkpwJFf71SscmpqUlcRdS+9i/u759G3el0+Gf0KjwEZujlop96hsUvjZGNNZRH4xxlxm37bDGNPRBbGWiyYFVdxPmSm8m3iCw7lZtPP2468hjbnEXmjPGMPMX2YyddlU/Dz9mDV0Fje0vsHNEStV9So7TyFXRCzY5w/Zv9FXu9nNSgF08w3ig8hWPBQaxfHcHP56cj//jDvKqbwcRITJnSezbco2ooKiGDpvKPd8fQ8ZuRnuDlupaqPUpCAiZ3vs/gssAsJF5FlgA7bRQ0pVS1YRrg+sz6eN23BrUAO+z0xm3PG9fJB4gvSCfNqEteHHO37koZ4P8fbWt+nyfhe2n9zu7rCVqhZKvX109raR/edLgX7YylKsMsbsqroQS6e3j5QjTuXl8GHSCValJxFi8WBicEOGBISyNj2RN88cIRkhO/0U3bPjeLnTLVjLKN2tVG1QoT6Fon0I1ZUmBVUee7MzeDsxll+z06lvsZJSUEBukSIa+XmZWH5fwMdX/I2mwU3dGKlSrlVWUiirzEW4iDxY2ovGGF2KU9Uobbz9eD3iIr7PTOa5uKMUn2hj9fAlu9kAOrzbgTcGvsH4juO1sJ6qc8rqaLYCAUBgKQ+lahwR4Wq/4FJHSnj7N6RjREcmfjmREfNHEJeus+ZV3VLWlcIJY8xzVRaJUlWogdWTUyUs4uMtwhtjv+HLbe/w5NonufTtS3n/hvcZ1maYG6JUquqVdaWg182q1iqpppIV21yGu04dJLXlKL6cvJWooCiGfz6ccYvHkZhZ5vIeStUKZSWFvlUWhVJVrF9AKA+HRhFh9USACKsn0+tHszDqUsbXi+DnrFT+k13A4GELeaj/K8z9dS7t3mnH1/u+dnfoSrmUQ7WPqisdfaRcJTU/j0WpZ1iYEke6KaCdGL7//kl2HPiaCR0n8Op1rxLiG+LuMJWqkMrOaFaqzgm02uYzzItqy8R6ERzGSkivf3HTqK9ZHLOFS9++lK9+/8rdYSrldHqloJQD0gryWZQSxwL7lUPOqW3s3PwSQxt35vWBrxPmF+buEJVymF4pKFVJARYrE4pcOYQ27EqXofP4tUF3Os8Zyrxd83SFN1Ur6JWCUhWQVpDPFylnmJd8kkzgzB/riE78jfevfYq9Fr/zFgXqFxDq7pCVKlTpRXaqK00Kyt3SCvJZlHyaOUnHybV4kp7wO4EhLSkosoiPtwgPh0ZpYlDVht4+UspFAixWJoREsrjpZQz39sK/WEIAyDaGD5NOuilCpcpHk4JSTuBvsfK3hpcgpSzzebqE2dNKVUeaFJRyogZWzxK3FxTk85cd8zmWfqaKI1KqfDQpKOVEJZbPMAbJOM3vwa259cR+hm+bzXendrspQqXKpklBKScqsXxGWBPWXjqQhz2hXuoREkIu4enMHPps+5S39y0nLz/P3WErVUhHHylVxfalnuSFwxs44BOB1TuI7MQDdCvI4LHWg4jw10lwyvV0SKpS1VBqbjYvH/6OdXkG8W9IdvppIpL38UCT7lzduMT/r0o5hSYFpaqxAmOYf+JXZiccIyMwivzcTPJPbGJUUAOmtB2On6efu0NUtYwmBaVqiB2pcbx6/GeOeNUHi5WU45u4LC+Zq1oOYVmBh86SVk6hSUGpGiY+L4c3Y37h+7x8Cjz9MabgnDkQOktaVYbOaFaqhqnv4cUzzbrz7UU9CBA5b1JctjG8nXDMTdGp2kyTglLVmJdYSC/laj6hwHDl/93Du9s+JCU7pYojU7WVJgWlqrnSZklbMHh1mMLsgIvotvZ5bl32IOuPrtcS3qpSPNwdgFKqbJODGzIjIYbsIh/23iI8GNqUYKsHs+Py2dVmNLFiZeqpn8n/+WNGhF3M7R1uI7petBsjVzWRJgWlqrmzncmlrdHQrWlnEvNz+b/kUyzKv5SUiM6sykllzqZ3aJJxnIkX9+emS27Soa3KITr6SKlaxBjDjux0Pk/4g83ZWRiLlZQzu0k49A1X+/gzsd0Yrml2DZZSqrmquqGs0Ud6paBULSIidPIJoFOjtqTm57EiPYH5Jp/TYZdyODeDuw8sJ/+75xkR1YVxHW7jhF+krhKnzqFXCkrVcsYY9uRk8FXKaVanJ5EnFtITD5ASv4eGzQcgVu/CfXX+Q92gVwpK1WEiQltvf9qGN+dv9fNZk57EEg9PDoZcfN6+2cbwXmKsJoU6TG8sKlWH+FmsXB9Ynw8bty11n7j8PK6fewNzds4hNTu1CqNT1YFbkoKIPCAiu0Vkl4jMFREfEQkVkZUist/+Z4g7YlOqrogoZf6DiJDY/m6ePLyBJu93YeT8kSzYvYCM3IwqjlC5Q5X3KYhIY2AD0NYYkyki84FvgLZAgjHmBRGZDoQYYx4tqy3tU1Cq4lalJZQ4/6GvXzAn8nLYnp2GQchM3E/sgf8j5dh3DIzuwai2oxjUcpAOca3BqmOfggfgKyK5gB8QCzwGXGt//WNgHVBmUlBKVdyF5j+cyctlXUYSq7x88Q1pCV0e4FDcTu7bMZf0b+5jYLOrGHnJSAa1HESAV4A7T0U5kVtGH4nIVOBfQCawwhhzq4gkGWOCi+yTaIw57xaSiEwBpgA0adLk8qNHj1ZV2ErVWcdzs1mdnsSq9ASO5eWAySf1xBaO7f+StBM/MqDZNYy4ZAQ3tLqBej713B2uuoBqVTrb3lewCLgZSAIWAAuBtxxJCkXp7SOlqpYxhoO5WaxKT2RNehJx+blIQS7Jxzdy9PcvSD25hd5Nr6ZzxzvZH9yK+IICnf9QDVW320f9gMPGmDgAEfkCuAI4JSKRxpgTIhIJnHZDbEqpMogIF3v5crGXL1OCI9mVnc7q9CTWefQlKPparPnZJCfsY6NfNJaCAgBO5efyUrytzLcmhurPHUnhD6CHiPhhu33UF9gKpAMTgBfsf37phtiUUg6yiNDBJ4AOPgHcF9qYrVmprElPYmWRyXBn5QDPHfuZ1cdXM7T1ULo17obVYq36oNUFVXlSMMZsFpGFwM9AHvAL8D4QAMwXkTuwJY5RVR2bUqpiPETo4RtED98gVqUnUtJNaatvGP87vo0XN79GmHcQQ1oOYWjrofRr0U87qqsRLXOhlHKqMTG/cSo/97ztAhjAagrwSD7Ivt+/IObwMiy5GVzb7Fqub3U9Q1oOoXlI8yqPua6pVh3NzqRJQanqp7T5Dw+GRNHA04sNGclsyEi2JQ5j8M88Tezh5ezdu5DM1D9oG96WwRcPZnDLwfRq0gvPUibZqYrTpKCUqlKr0hLKrL56dhTTRnuCOJCbBUBgbjppJzax/dfZJMbtJNArgH4t+jHo4kEMajmIqKAod51SraJJQSlVrZ3My2FjRjIbM1LYkZ1GAeBv8vCI38uevfPZf+gbTEEenTveSVi7CeR4+NPAw5M7gyN1RFMFaFJQStUYKfl5/JiZwobMFLZkppJlCvABfLKTSPIMAEuR8TH5OXRO2c/tjTrSNrwtIuK2uGuS6jZPQSmlShVk9WBAQCgDAkLJLijg56w0NmQm821JO1u9+MEzhFfeaUejwEb0b9Gf/i3607dFXxoGNKzq0GsFvVJQStUIfY7uKHGoK0Cz7DMkHN/I5p2zOZV0AIB2DdrRt3lf+rXox9VNrybIO6jqgq3m9PaRUqrGK22oq7cIPmIhuSAfgEYYvFMOc/Twcjbt/ozMnFSsYqVr4670adaHPs37cEX0Ffh6+lb1KVQbmhSUUjVeaUNdHw6Noo9/CPtzMtmSlcqWzFR2Z6eTD/iIEF2QTd6Z3ezdv5if9n9NvsnHy+pFj6ge9G7Wm97NetM9qjs+Hj7uO7kqpklBKVUrXGio61npBflsz0pjS2YqW7JSic3LAaCB1YOGOcmkn9zKL3vn88uxHwhvPpAWl/8NH/+GeOdlMthawJ1RnWv1lYQmBaVUnXY8N5ut9quIX7LSyDAFWIAIi5VT+XkUFBm1lJ+XyYEf/kWL3ESubno1VzW5iiuir6hVJcE1KSillF2eMezOTmdLZiqfp8SRV0L3tUdeBmnrHmJr7BbyCvKwiIUOER24qslVXNXkKq5sciWNAhu5IXrn0KSglFIlKGtEk79YaOvlQ2BmHIknt7D9wFI2H/uhcK3q5sHNubLJlVwZbXtc2uBSLOKWZe/LTecpKKVUCRpYPUsc0RRksXK1Xz12ZKWzxeIPja7Fp3Ef7vbyJTw3mdTT2/nt0DJWHlzJpzs/BaCedz16RPWgZ1RProi+gu5R3WvkMFi9UlBK1VlljWg624GdkJ/Ljqx0dmSlsSM7nSP2Ok3eIrT18qeJySEnfg8Hjq7mxz++Z9fpXRgMDZoPonXXB7H6huFfkMNIHx/GR7arFlcTevtIKaVK4eiIprOS8/PYkZ1WmCgO5WZhAE+ES7z9aOPhybHUGDbjQ0GRkhz5eZn88dNLXJybTPfG3eke1Z1ujbu5Zea1JgWllHKR1Pw8fs1OZ0e2LUnsz8mkoJR9PXJSSVg7lZ2ndpJXkAdAdFA03Rp3o2ujrnRt3JXLIy93+UgnTQpKKVVF0gvyuf7YrlJf7+0XTEsPTyTlKDGxm9kWu5mfjv/EocRDhfu0rt+aLo26FD46Nezk1NXpNCkopVQVKqskR7DFo/A1T4SLvXxp6+1HNAVkJfzOvtgf2Ra7lS2xW4hNjQVAENqEteHyRpdzeaTt0alhJwK9AysUn44+UkqpKjQ5uGGZHdhn8nLZk5PBb9np/JadwdK0eNu+nuGENB9Ouza3Mtrbj/D8TFLidvHria1sO7GN1YdWsyI1nhaBrfE5cZAID68L9oGUl14pKKWUC5SnAzvPGA7nZvJbdga/ZWewJzuDY3nZAFiAZp4+tPX2AwPL0xMoeg1SfLSUI/T2kVJK1TAp+Xnszcn4M1HkZJBmrwRbXITVk3lRbR1uW28fKaVUDRNk9aCbbxDdfG0T4AqMoe8fO0vc93QJ/RcV5f5ZFEoppS7IIkKE1bPE1xqUsr1Cx3FaS0oppVxqcnBDvIutQ+0twuRg502A09tHSilVQ5ztTC7PDOzy0qSglFI1SL+AUKcmgeL09pFSSqlCmhSUUkoV0qSglFKqkCYFpZRShTQpKKWUKlSjy1yISBxwtBJNhAFnnBSOO9WW8wA9l+qotpwH6Lmc1dQYE17SCzU6KVSWiGwtrf5HTVJbzgP0XKqj2nIeoOfiCL19pJRSqpAmBaWUUoXqelJ4390BOEltOQ/Qc6mOast5gJ7LBdXpPgWllFLnqutXCkoppYrQpKCUUqpQnU4KIvIPEdkpIttFZIWINHJ3TBUlIi+JyF77+SwWkWB3x1RRIjJKRHaLSIGI1LjhgyIyUER+F5EDIjLd3fFUlIjMEpHTIrLL3bFUlohEi8haEdlj/7c11d0xVYSI+IjITyKyw34ezzr9GHW5T0FEgowxKfaf/wa0Ncb8xc1hVYiIDADWGGPyROQ/AMaYR90cVoWIyCVAAfAe8LAxpsYsxC0iVmAf0B+IAbYAY40xv7k1sAoQkauBNOATY0w7d8dTGSISCUQaY34WkUBgGzCspv29iIgA/saYNBHxBDYAU40xPzrrGHX6SuFsQrDzB2pshjTGrDDG5Nmf/ghEuTOeyjDG7DHG/O7uOCqoG3DAGHPIGJMDzANudHNMFWKMWQ8kuDsOZzDGnDDG/Gz/ORXYAzR2b1TlZ2zS7E897Q+nfm7V6aQAICL/EpFjwK3A390dj5PcDnzr7iDqqMbAsSLPY6iBHz61mYg0Ay4DNrs3kooREauIbAdOAyuNMU49j1qfFERklYjsKuFxI4Ax5gljTDQwB7jXvdGW7ULnYt/nCSAP2/lUW46cSw0lJWyrsVegtY2IBACLgPuL3SmoMYwx+caYTtjuBnQTEafe2qv1y3EaY/o5uOtnwNfA0y4Mp1IudC4iMgG4HuhrqnlnUTn+XmqaGCC6yPMoINZNsagi7PfgFwFzjDFfuDueyjLGJInIOmAg4LTBALX+SqEsItKyyNOhwF53xVJZIjIQeBQYaozJcHc8ddgWoKWINBcRL2AM8JWbY6rz7B20M4E9xphX3B1PRYlI+NmRhSLiC/TDyZ9bdX300SKgNbaRLkeBvxhjjrs3qooRkQOANxBv3/RjDR5JNRx4EwgHkoDtxpjr3BuV40RkMPAaYAVmGWP+5eaQKkRE5gLXYivRfAp42hgz061BVZCI9AK+B37F9v8d4HFjzDfui6r8RKQD8DG2f1sWYL4x5jmnHqMuJwWllFLnqtO3j5RSSp1Lk4JSSqlCmhSUUkoV0qSglFKqkCYFpZRShTQpKJcSkWburrIpIgtFpIX957QL7V/VxOYNe1XVnSLSuZT9vnG0+m1plWZFpL69WmiaiLzlxHMIF5FlzmpPuY8mBVUjiYhDs/FF5FLAaow55MJYrJVsYhDQ0v6YArxT0k7GmMHGmCQH29wF3ASsL7Y9C3gKeLhioZbMGBMHnBCRK53Zrqp6mhRUVbCKyAf2b64r7DMxEZFOIvJjkTUgQuzb1539disiYSJyxP7zRBFZICL/B6wQkUgRWS+29TB2ichVJRz7VuDLohvsRRB32I8dYd/WVERW22NZLSJN7Ns/EpGRRd6bZv/zWvs37s+wTYiqjBuxlac29hLIwfZSz+cQkSP234e/iHxtP4ddInJz8X1LqzRrjEk3xmzAlhyKt/+OiGyVYnX67cd9XkQ22V/vLCLLReSgiBSdILkE2+9b1WCaFFRVaAn81xhzKbYZyiPs2z8BHjXGdMD2wepI3amewARjTB/gFmC5vThYR2B7Cftfia12/ln+2GZ7d8T2LfpO+/a3sH0wd8BWTPANB2LpBjxhjGlb/AUR+dyerIo/xpfQTnkrqw4EYo0xHe3rHDjrts0TxpguQAfgGvvs2bOOGWN6YpsV/BEwEugBFJ1NuxUoKTGrGqTWF8RT1cJhY8zZD+xtQDMRqQcEG2O+s2//GFjgQFsrjTFna/xvAWbZC50tKXKMoiKBuCLPc4ClRWLpb/+5J7bbLQCzgRcdiOUnY8zhkl4wxpz37b0M5a2s+iswQ2yLKS01xnxfjmOVZbSITMH2uRAJtAV22l87W7/pVyDAviZBqohkiUiw/bbWaaDGrl6obPRKQVWF7CI/53PhLyN5/Plv06fYa+lnf7AvAnM1cByYXcq38MxibeQWqSBbVixn9ymMxV5UzaukWIor55VCuSqrGmP2AZdj+4D+t4hUeh0QEWmOrZ+hr/1q6WvO/b2d/Tss4Ny/zwL+/B36YPt9qxpMk4JyC2NMMpBYpB9gHHD2quEItg89sN2mKJGINAVOG2M+wFYBs6RRO3uAix0I6QdsFU3Bdl98Qwmx3IhtpasLMsbcbIzpVMLjkxJ2/woYbx+F1ANINsacKK1tsa0lnmGM+RSYQcnnXV5B2JJcsr2fZVAF2miFE0s4K/fQ20fKnSYA74qIH3AImGTfPgOYLyLjgDVlvP9aYJqI5GJbS7ikb+Ff2/dbdYFY/obtVtQ0bLebzsbyAfCliPwErKaMq4NK+AYYDBwAMoocuzTtgZdEpADIBe4uvoOcW2n2axEprDRr77gPArxEZBgwwBizQ0R+AXZj+7vYWIHz6I3t961qMK2Sqmo1+0intcCVxph8d8dTm4nIeuBGY0yiu2NRFadJQdV6InIdtsVV/nB3LLWViIRjS7xL3B2LqhxNCkoppQppR7NSSqlCmhSUUkoV0qSglFKqkCYFpZRShTQpKKWUKvT/DAORb7kBDgQAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"##BD means Before Death\n",
"N = 100\n",
"K = .275\n",
"time_BD = np.array([0, -1, -2, -3])\n",
"timing_BD = np.linspace(0, -3, N)\n",
"dt_BD = timing_BD[0]-timing_BD[1]\n",
"T_a_BD = np.array([65, 60, 58, 55])\n",
"T_BD = np.zeros(N)\n",
"T_BD[0] = 85 \n",
"for i in range(1, len(timing_BD)):\n",
" for j in range(0,3):\n",
" T_BD[i]= T_BD[i-1] + K*(T_BD[i-1]-T_a_BD[j])*dt_BD\n",
"\n",
"\n",
"##AD means After Death\n",
"timing_AD = np.linspace(0, 3, N)\n",
"time_AD = np.array([0, 1, 2])\n",
"dt_AD = timing_AD[0]-timing_AD[1]\n",
"T_a_AD = np.array([65, 66, 67])\n",
"T_AD = np.zeros(N)\n",
"T_AD[0] = 85\n",
"for i in range(1, len(timing_AD)):\n",
" for j in range(0,2):\n",
" T_AD[i]= T_AD[i-1] + K*(T_AD[i-1]-T_a_AD[j])*dt_AD\n",
"plt.plot(timing_BD, T_BD, color = \"green\", label=\"Numberical Before Death\")\n",
"plt.plot(timing_AD, T_AD, color=\"green\", label=\"Numberical After Death\")\n",
"plt.xlabel('hours (hour = 0 is 11am)')\n",
"plt.ylabel('Temperature (F)')\n",
"\n",
"N = 15\n",
"t = np.linspace(-3, 3, N)\n",
"dt= t[1]-t[0]\n",
"T_analytical = np.zeros(len(t))\n",
"T_0 = 85\n",
"T_a = 65\n",
"K = 0.275\n",
"\n",
"for i in range(0, len(t)):\n",
" T_analytical[i] = T_a + (T_0-T_a)*np.exp(-K*t[i])\n",
"plt.plot(t, T_analytical, '-o', color='turquoise', label= \"Analytical\")\n",
"\n",
"plt.legend(loc=\"best\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The nonlinear Euler Approximation makes more logical sense because the ambient temperature in reality would not be constant as first assumed. This graph has an exponential decay curvature which is what is expected. However, before death does not really make sense because the further you go back before his death the hotter the person gets.\n",
"\n",
"The time of death is 9:30 am. "
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Temperature is 98.637 -1.494949494949495 hours before 11am\n"
]
}
],
"source": [
"for i in range(len(timing_BD)):\n",
" temp = T_BD[i].round(3)\n",
" if temp> 98.4 and temp<98.8:\n",
" print(\"Temperature is\", temp, timing_BD[i], \"hours before 11am\")"
]
}
],
"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
}