Permalink
Cannot retrieve contributors at this time
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-project2.ipynb
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
419 lines (419 sloc)
64.9 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "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": 178, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "K= 0.275 1/hr\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "# 1\n", | |
| "Ti=85 #initial recorded temperature\n", | |
| "Tf=74 #final recorded temperature\n", | |
| "Ta=65 #ambient temperature\n", | |
| "k=-((Tf-Ti)/2)*(1/(Ti-Ta)) #emperical constant\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": 179, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "K= 0.275 1/hr\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "# 2\n", | |
| "import numpy as np\n", | |
| "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=Initial temperature (F)\n", | |
| " Temp_t2=Final temperature (F)\n", | |
| " Temp_ambient=Ambient temperature (F)\n", | |
| " delta_t=time difference between initial temperature and final temperature (hr)\n", | |
| " \n", | |
| " Returns\n", | |
| " -------\n", | |
| " k=empirical constant (1/hr)\n", | |
| " '''\n", | |
| " delta_T=Temp_t2-Temp_t1\n", | |
| " k=-(delta_T/delta_t)/(Temp_t1-Temp_ambient)\n", | |
| " return k\n", | |
| " \n", | |
| "print('K= ',measure_K(85,74,65,2),'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": 182, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x7fd02da83fd0>" | |
| ] | |
| }, | |
| "execution_count": 182, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUZfbA8e+hV+kiRQGVIhAIEIoUAVEEFBQVgVWqwlpQcV1d/K1r72LXXRuoq4AVlFUsiBQRRUIRQUSKgBQREQQEpOT8/jgzIYRJMimTm3I+zzNPMnfm3nsmYk7edl5RVZxzzrnUigQdgHPOubzJE4RzzrmIPEE455yLyBOEc865iDxBOOeci6hY0AHkpKpVq2rdunWDDsM55/KNhQsX/qqq1SK9VqASRN26dUlMTAw6DOecyzdEZH1ar3kXk3POuYg8QTjnnIvIE4RzzrmICtQYhHN50cGDB9m4cSP79+8POhRXiJUqVYratWtTvHjxqM/xBOFcjG3cuJHy5ctTt25dRCTocFwhpKps376djRs3Uq9evajPi2kXk4jcICLLRWSZiEwSkVIicoeIbBKRJaFHrzTO7SEiK0VktYiMiWWczsXS/v37qVKliicHFxgRoUqVKpluxcYsQYhILeA6IEFVmwJFgQGhlx9T1fjQY1qEc4sCzwA9gcbAQBFpHKtYnYs1Tw4uaFn5NxjrQepiQGkRKQaUATZHeV4bYLWqrlXVA8DrwPkxifDgQe7u/jnTHlnBoUMxuYNzzuVLMUsQqroJGAtsALYAv6vqJ6GXR4nIUhEZLyKVIpxeC/gpxfONoWPHEJGRIpIoIonbtm3LdJx7tv/JMzMacu7fT6NOHeWWW+CHHzJ9GefyvClTpiAifP/991m+xtChQ3n77bfTfc9999131PP27dtn6V533HEHY8eOzdK5LmfEsoupEvZXfz2gJlBWRC4D/gOcAsRjieORSKdHOBZxZyNVfV5VE1Q1oVq1iKvF01XuhHJsmDCXyfSlVeUfefhhaNgQOnWC8eNh9+5MX9K5PGnSpEl07NiR119/Pab3SZ0g5s2bF9P7udiJZRfTWcCPqrpNVQ8Ck4H2qrpVVQ+rahLwAtadlNpG4MQUz2sTffdUppXo35e+PfYzdX08PyVu5cEH4ddf4fLLoUYNGD4cPv8cfPM9l1/t2bOHL774gnHjxiUniFmzZtGlSxcuvvhiGjVqxKWXXkp4h8m77rqL1q1b07RpU0aOHEnqnSdnzJhB3759k59Pnz6dCy+8kDFjxrBv3z7i4+O59NJLAShXrlzy+x566CHi4uJo3rw5Y8bY3JMXXniB1q1b07x5cy666CL27t0b05+Fi14sp7luANqJSBlgH9ANSBSRGqq6JfSevsCyCOcuAOqLSD1gEza4/ZeYRSoCTz0FTZtS46EbuHniRG66Cb76yloRb7wBL70E9evDsGEweDDUitjh5VwGRo+GJUty9prx8fD44+m+5d1336VHjx40aNCAypUrs2jRIgAWL17M8uXLqVmzJh06dOCLL76gY8eOjBo1ittuuw2AQYMG8f7779O7d+/k65155plcc801bNu2jWrVqvHSSy8xbNgwevfuzdNPP82SCJ/xww8/5N1332X+/PmUKVOG3377DYALL7yQESNGAHDrrbcybtw4rr322hz50bjsieUYxHzgbWAR8G3oXs8DD4nItyKyFOgK3AAgIjVFZFro3EPAKOBjYAXwpqouj1WsAJx6KowZA5MmwYwZiMDpp8MLL8CWLfDKK1CzJvzf/8FJJ0GvXvD22/DnnzGNyrkcMWnSJAYMsEmEAwYMYNKkSQC0adOG2rVrU6RIEeLj41m3bh0AM2fOpG3btsTFxfHZZ5+xfPnR//uJCIMGDeK1115j586dfPnll/Ts2TPdGD799FOGDRtGmTJlAKhcuTIAy5Yto1OnTsTFxTFhwoRj7uWCE9OFcqp6O3B7qsOD0njvZqBXiufTgGOmwMbUP/4Br70G11wDS5dCiRIAlC1rrYbBg2HNGnj5ZXv06weVK8Nll1k3VPPmuRqty48y+Es/FrZv385nn33GsmXLEBEOHz6MiNCrVy9KliyZ/L6iRYty6NAh9u/fz9VXX01iYiInnngid9xxR8T58+EWQ6lSpejXrx/FiqX/60RVI061HDp0KO+++y7Nmzfn5ZdfZtasWdn+zC5neC2mlEqXtq6mlSvhkUhj53DKKXD33bBuHXz0EZx9Njz7rLXyW7aEp5+GUMvZuTzh7bffZvDgwaxfv55169bx008/Ua9ePebOnRvx/eFkULVqVfbs2ZPmrKWaNWtSs2ZN7rnnHoYOHZp8vHjx4hw8ePCY93fv3p3x48cnjzGEu5h2795NjRo1OHjwIBMmTMjOR3U5zBNEaj17woUXHskCaShaFM45B15/3bqgnnrKhjKuvdYGtvv3h48/hsOHcy905yKZNGnSUQPKABdddBETJ06M+P6KFSsyYsQI4uLiuOCCC2jdunWa17700ks58cQTadz4yDrWkSNH0qxZs+RB6rAePXrQp08fEhISiI+PT57Cevfdd9O2bVvOPvtsGjVqlNWP6WJAUs9OyM8SEhI0RzYM+uknOO006NYN3nsvU6cuWWID2hMmwPbtULs2DBlig9unnJL90Fz+s2LFCk477bSgw4iJUaNG0aJFCy6//PKgQ3FRiPRvUUQWqmpCpPd7CyKSE0+E22+HqVPhf//L1Knx8fDEE7BpE7z1FsTFwf332xh4ly422P3HH7EJ27nc1KpVK5YuXcpll10WdCguRjxBpGX0aGjcGK67DrIwL7tkSbj4Ypg2DTZsgPvus6QxdKh1QY0YAfPm+doKl38tXLiQOXPmHDXQ7QoWTxBpKV4c/v1vG4dItTI0s2rVIrmEx+efW+KYNAk6dLCerIcesnEM55zLSzxBpKdzZxg0CB5+2GY2ZZMIdOxoi++2bIFx46BaNZtde+KJ0Lu3dUvt25cDsTvnXDZ5gsjIww/b9NdRo3K0P6h8+SMlPFauhJtugkWL4JJL4IQTrMzHZ5/5LCjnXHA8QWSkenW491749FN4882Y3KJBAxvI3rDBbnPhhdaS6NYN6tSBm2+2dXvOOZebPEFE48orbRXcDTfArl0xu03RopYUXnoJtm61GlAtW8Jjj9kq7bg4ePBBm4XrXGaICDfeeGPy87Fjx3LHHXfkagyJiYlcd911WTq3S5cupDWFfdu2bRQvXpznnnsuO+Fl29SpU3nggQdy5FqPP/74UUULe/Xqxc6dO3Pk2pnhCSIaRYvCf/4DP/8MufQ/VenS1t00daqNV/z733DccVYu6qSTbMrsiy9CAP9mXD5UsmRJJk+ezK+//hrI/Q8dOkRCQgJPPvlkjl/7rbfeol27dsn1pXLCoSzsHtanT5/kCrXZlTpBTJs2jYoVK+bItTPDE0S02rSBv/4VnnwSvvkmV29dtSpcdRV88QWsXg133WVJY8QI6wG76CKYMsULB7q0FStWjJEjR/LYY48d81rqTYDC5blnzZpF586dueSSS2jQoAFjxoxhwoQJtGnThri4ONasWQPYX/AXXXQRrVu3pnXr1nzxxReAbfgzcuRIunfvzuDBg5k1axbnnXceYOXHhw0bRlxcHM2aNeOdd94B4KqrriIhIYEmTZpw++2py7hFNmnSJB555BE2btzIpk2bjvocN954Iy1btqRbt26ENxTr0qULo0ePpn379jRt2pSvv/46Yrz79+9PjrFFixbMnDkTgEcffZThw4cD8O2339K0aVP27t3Lyy+/zKhRo5J/pldddRVdu3bl5JNPZvbs2QwfPpzTTjvtqLIkkT7vk08+yebNm+natStdu3YFoG7dusnJ/dFHH6Vp06Y0bdqUx0O1vdatW8dpp53GiBEjaNKkCd27d2dfTsx2UdUC82jVqpXG1PbtqlWrqrZvr3r4cGzvlYGkJNUFC1Svv171+ONVQbViRdURI1Rnzw48PJfCd999l/z99derdu6cs4/rr884hrJly+rvv/+uderU0Z07d+rDDz+st99+u6qqDhkyRN96662j3quqOnPmTK1QoYJu3rxZ9+/frzVr1tTbbrtNVVUff/xxvT5044EDB+rnn3+uqqrr16/XRo0aqarq7bffri1bttS9e/cmX+/cc89VVdWbb745+XxV1d9++01VVbdv366qqocOHdLOnTvrN998o6qqnTt31gULFhzzuTZs2KCnnnqqqqrecsst+sgjjyS/Buhrr72mqqp33nmnXnPNNcnXuuKKK1RVdfbs2dqkSZOI8Y4dO1aHDh2qqqorVqzQE088Ufft26eHDx/WTp066eTJk7VVq1Y6d+5cVVV96aWXku8xZMgQ7d+/vyYlJem7776r5cuX16VLl+rhw4e1ZcuWunjx4nQ/b506dXTbtm3JnyX8PDExUZs2bap79uzR3bt3a+PGjXXRokX6448/atGiRZOv269fP3311VeP+Xml/LeY4ueUqGn8TvUWRGZUrmyzmubNsyXRARKBhAQrDrppkxUOPO88mDjRZufWq2elyb1ysgs77rjjGDx4cKa6eVq3bk2NGjUoWbIkp5xyCt27dwcgLi4uuTT4p59+yqhRo4iPj6dPnz7s2rWL3aGtGPv06UPp0qWPue6nn37KNddck/y8UiXbefjNN9+kZcuWtGjRguXLl/Pdd9+lG9/rr7/OJZdcAhxdxhygSJEi9O/fH4DLLrvsqOKEAwcOBOCMM85g165dyf37KeOdO3cugwZZ8elGjRpRp04dfvjhB4oUKcLLL7/MoEGD6Ny5Mx06dIgYW+/evRER4uLiqF69OnFxcRQpUoQmTZok/+wy+3nnzp1L3759KVu2LOXKlePCCy/k888/B6BevXrEx8cDtsp9XTq15KIV03LfBdLgwbaA4aaboE8fqFIl6IgoVswKB55zjpXxeO89q1r+0EM2Oyo+3kqSDxxoe1q44ARQ7fsoo0ePpmXLlgwbNiz5WLFixUhKSgKsR+HAgQPJr6VcJV2kSJHk50WKFEnup09KSuLLL7+MmAjKli0bMQ6NUPr7xx9/ZOzYsSxYsIBKlSoxdOjQiGXGU5o0aRJbt25NrgK7efNmVq1aRf369Y95b8r7pb53+HnKeDWdae2rVq2iXLlybN6c9kaXKX9WqX+Ohw4dytLnTS+m1KXbc6KLyVsQmVWkiI0Y79xpf6LnMWXLwl/+YiU+Nm2yulAlSsDf/26FA886y/ayiOFkLJeHVa5cmUsuuYRx48YlH6tbty4LFy4E4L333otYqjs93bt35+mnn05+Hmk3uYzO2bFjB7t27aJs2bJUqFCBrVu38uGHH6Z7jZUrV/LHH3+wadMm1q1bx7p167jllluSt1RNSkpKHluZOHEiHTt2TD73jTfeAOwv8goVKlChQoVjrn/GGWckJ54ffviBDRs20LBhQ37//Xeuv/565syZw/bt29Msh56R9D5v+fLlk1thqWN699132bt3L3/88QdTpkyhU6dOWbp/NDxBZEVcHFx/vW03N39+0NGkqXp1KyU1f74txvvXv6xyyLBh9lr//laLMMUfjK4QuPHGG4+azTRixAhmz55NmzZtmD9/fpp/9aflySefJDExkWbNmtG4cWOeffbZDM+59dZb2bFjB02bNqV58+bMnDmT5s2b06JFC5o0acLw4cPT7LoJS6uMebibqWzZsixfvpxWrVrx2WefJW+hCtal1b59e6688sqjkmVKV199NYcPHyYuLo7+/fvz8ssvU7JkSW644QauvvpqGjRowLhx4xgzZgy//PJLhp85tfQ+78iRI+nZs2fyIHVYy5YtGTp0KG3atKFt27ZcccUVtGjRItP3jlZMy32LyA3AFYBi244OA+4GegMHgDXAMFU9ZrKmiKwDdgOHgUOaRjnalHKs3Hc0du+GRo3sN+2CBTYVNh9QtYTx2mu2l8X27dZLdskl1g11+uk2vuFyTkEu952XlStXjj179hxzvEuXLowdO5aEhAx/pRQ4eabct4jUAq4DElS1KVAUGABMB5qqajPgB+CWdC7TVVXjo0kOua58eetQXrzY1kjkEyLQrp3tfLdlC7z/vu2K99JLVjzw1FOtpbFsWdCROueCFusupmJAaREpBpQBNqvqJ6oaXoXyFVA7xjHEzsUX22/Xf/7TFtHlM8WLw7nnWmXZrVttbOKUU6x4bVwcNGkCd94JK1YEHalzmRep9QC2vqMwth6yImYJQlU3AWOBDcAW4HdV/STV24YDaY1EKfCJiCwUkZGxijNbROxP8f37bVZTPnbccbbz3Sef2OD2M89Ypdk777RtMeLibBfWHChqWyjFsivXuWhk5d9gLLuYKgHnA/WAmkBZEbksxev/BA4Bae1S3kFVWwI9gWtE5Iw07jNSRBJFJDG8UjJXNWhg9bpfew1mzcr9+8fACSfA1Vfbx9m0yRaPV6wIt91mwy7x8dbKWLUq6Ejzh1KlSrF9+3ZPEi4wqsr27dspVapUps6L2SC1iPQDeqjq5aHng4F2qnq1iAwBrgS6qWqG27WJyB3AHlUdm977cnWQOqV9+6w/plQp25S6RIncjyEXbNoEb79tRW3nzbNjLVrYAPcll8DJJwcbX1518OBBNm7cmOEcd+diqVSpUtSuXZvixYsfdTy9QepYJoi2wHigNbAPeBlIBFYBjwKdVTXin/wiUhYooqq7Q99PB+5S1Y/Su2dgCQJstLd3b3jgAWtRFHA//XQkWXz1lR1LSLBE0a8f1K0baHjOuSgFkiBCN74T6I91JS3GprwuB0oC20Nv+0pVrxSRmsCLqtpLRE4GpoReLwZMVNV7M7pfoAkCoG9f68RfscJKrhYS69fb/hVvvmkzfsFqG4aTRSH6UTiX7wSWIHJb4Ali/Xob0T3nHJg8Obg4AvTjj0eSRWhxLu3aHUkWtfPvnDXnCqRA1kEUSnXq2CKCKVOs1kUhVK+e7YCXmGilye+7zyZ5/e1vtu92x4426J1OCRvnXB7hLYicduCAbf924ICtNotQwKww+uGHIy2LpUtthnDHjtayuPhimznlnMt93oLITSVKWDG/tWttwNoBNhv4n/+0vZZWrLD1FTt2wLXXWoXZrl1tQfrWrUFH6pwL8xZErFx6qU3zWbYMIpQedmb5cmtZvPEGfP+9Fcvt0sVaFhdeaIv1nHOx44PUQdiyxVaVtWtnu/l4Bbx0qVqyePNNSxY//GDJolMnmxx2wQU2xOOcy1nexRSEGjXgnnts2mtov12XNhFo2tT22/7+e+uK+uc/rdrs6NG2rqJlSyv3sWyZJRTnXGx5CyKWDh2yBQG//GId7+XLBx1RvrR6tU0MmzLFFuWpWtXZcMuiXTtrbTjnMs9bEEEpVswGrDdtslFZlyWnnmq1EOfNsx/ls89a1dnHH7cS5bVqwZVXwscf++ZHzuUkb0HkhpEjYfx4q9PUtGnQ0RQYO3facpMpU+DDD20/7goVrIT5BRdAz55QrlzQUTqXt/kgddC2b4eGDeG002DOHB+wjoH9++HTTy1ZTJ0Kv/4KJUvadh19+0KfPlC1atBROpf3eBdT0KpUgQcfhLlz4b//DTqaAqlUKTjvPBg3ziaQzZpl3U5Ll8Lll9vOsF26wBNPWEUU51zGvAWRW5KSbOnw6tW2606lSkFHVCioWs9eeJA7vJVqixbWsujb1yq1e6POFVbegsgLihSxpcLbt9v8TZcrRCwZ3HUXfPutbXL00EPW4rj9dtspr0EDqx/15ZeWx51zxlsQuW30aKtWN38+tG4ddDSF2pYtNl4xZQp89hkcPGg1oc4/31oWXbsW2L2fnEvmg9R5ya5dtsK6Zk1LEkWLBh2RI/0ZUeefbxXcK1QIOkrncp53MeUlxx0Hjz5qmyU8/3zQ0biQihXhL3+xulC//gr/+x9cdJEthO/f32ZAdesGjz1mZUCcKwy8BREEVZt/mZhoA9bVqwcdkUvD4cO2evv99+0RHuSuX99mTZ13ns098K4ol195F1NetHKljZAOHAivvBJ0NC5K69bBBx9YsvjsM1u5fdxx1gV13nm2OM8r0Lr8JLAuJhG5QUSWi8gyEZkkIqVEpLKITBeRVaGvEed7ikgPEVkpIqtFZEws4wxEw4ZWP+K//7XFcy5fqFsXrrnGxim2b4d337XS5HPnwpAh1hhs39520lu61IsKuvwtwxaEiMQDnYCawD5gGTBDVX/P4LxawFygsaruE5E3gWlAY+A3VX0g9Iu/kqr+I9W5RYEfgLOBjcACYKCqfpfePfNVCwJg717bw7pcOVi8GIoXDzoil0VJSbbeItwVtWCBHa9d+0hX1Jln+gaDLu/JUgtCRC4TkYXAnUAlYD2wCzgLmCUi40Qkoy3oiwGlRaQYUAbYDJwPhPtUXgEuiHBeG2C1qq5V1QPA66HzCpYyZWzK6/LltsTX5VtFilg58ttug6+/tim048bZTOZXX7UEUaUK9O4Nzz0HGzcGHbFzGUuzBSEi1wMvquofabyegP31Pz3Ni9s17sVaHp+o6qUislNVK6Z4zw5VrZTqvIuBHqp6Rej5IKCtqo6KcI+RwEiAk046qdX6/FhHoU8f69BesQJOPDHoaFwO+/NPmD37SOvixx/teHz8kdZF69ZestwFI0stCFV9Iq3kEHo9MYPkUAn7q78e1j1VVkQuizbmSLdMI47nVTVBVROq5dfRwSeesD6KG24IOhIXAyVLQvfu1lhcswa++85Wcx93HNx/v+1nccIJMHSo7VK7a1fQETtn0uti+jDF9zdn4dpnAT+q6jZVPQhMBtoDW0WkRui6NYBfIpy7EUj5p3RtrHuqYKpXD2691Xaemzw56GhcDIlYUd+bbrJWxS+/wMSJNut56lTo18/WXJx1lu13sXp10BG7wiy9Ru0JKb4fkIVrbwDaiUgZERGgG7ACmAoMCb1nCPBehHMXAPVFpJ6IlAjdf2oWYsg/brzR+hkGD7aiQa5QqFzZZjpPmGDJYs4ca0hu2WJf69e3hfd//zvMmGHdVc7llvTGIBapasvU32fq4iJ3Av2BQ8Bi4AqgHPAmcBKWRPqp6m8iUhMb8+gVOrcX8DhQFBivqvdmdL98N4sptc2bLUmUKGEjnfm1y8zliLVrbc3FBx/AzJm25qJMGasRdc459qhf3yvRuuzJ0kI5EdkJfIaNB3QNfZ9MVS/M4TizLd8nCLD5kWecYYni0099ia4DYM8eSxIff2yPcNdTvXpHksWZZ9q4hnOZkdUE0S29i6rqjByILUcViAQB1il96aVwxRVWr8n/RHSprFlzJFl89pklkGLFbJFeOGG0aOEzo1zGvNRGfvTPf9py3CeegOuuCzoal4cdOADz5h1JGIsX2/Fq1Wz21Dnn2Fcv+eUiyVaCEJHFHDvF9HcgEbhfVX/LkShzQIFKEElJcOGFVlb0o49smotzUdi61arQfvyxfd22zY7Hx0OPHpYw2rf33ktnspsgHsLGISaGDg0ADgN7gHaq2icHY82WApUgAHbvtv+TN260vSMaNAg6IpfPJCVZiyLcupg3Dw4dsuouZ555pDvqlFOCjtQFJbsJYq6qdox0TES+VdW4HIw1WwpcggBbdtumjc2HnD/fNi5wLot27bLB7o8+soQRXtV9yilHWhddu1oCcYVDdqu5lheRViku1hIIz5U4lAPxufTUq2cL6NauhQED7M8/57LouONsh7z//McGun/4AZ56ytZavPSSVX2pXNmSxIMPWgHCAjRM6TIpmhZEO2A8UBzrajoAXA4sBfqo6qRYBxmtAtmCCHvxRRgxwlZPPfpo0NG4AujPP+GLL460LpYutePVq9sgd48eNhTmy3MKlhyZxSQiVULv/zUng8tJBTpBAFx/vRX0GTcOhg8POhpXwG3ZYoPcH30E06fb/hci0Ly5bb/arRt06uTdUflddscgqgH3ALVU9TwRaQy0UdWXczzSbCrwCeLQIejVC2bNssnvHTtmeIpzOeHwYVi0yFoWM2bYYPeBA7aFSdu2RxJG27Y+Oyq/yW6C+ACYAPxDVZuLSHFgUV4anA4r8AkCYMcO+79w505bdV2nTtARuUJo717rjpoxwx4LF9pYRdmy1qoIJ4zmzX2xXl6X3QSxQFVbi8hiVW0ROrZEVeNjEGu2FIoEAbafddu2lhy++MLb+C5wO3ZYwzacML7/3o5XqWID3uGEceqpXhggr0kvQRSL4vw/RKQyocVyItIa2J2D8bnMatgQ3njDupsGD7ZNBPzPNBegSpWgb197AGzaZL2g4YTx9tt2/MQTjySLM8+EmjWDi9llLJoWRALwBNAE+AaoBVysqktiH17mFJoWRNjjj9uspltvhbvvDjoa5yJShVWrjiSLmTPht1D9hUaNjiSMLl0s0bjcle1ZTKE9GU7Dprl+F9onOs8pdAlC1Qr6jR8Pr78O/fsHHZFzGUpKgm++OZIw5syxMY3wvt7hhNGhg5U3d7GV1Wqu6ZbQUNU8t4FPoUsQYJPXu3WzUcLPP4eEiP+dncuzDhywIgHhhPHVVzZhr0QJqzQTThitW1vFWpezspogXg19WxXbKnQm1oLoDMxW1fNjEGu2FMoEAbYVWevWNhdxwQKoUSPoiJzLsj177G+dcMJYEurMLl8eOnc+Mn7RtKkPveWE7M5imgpcpaqbQs9rAU+q6kU5Hmk2FdoEAdZm79ABmjSx6SSlSwcdkXM54tdfj54htWqVHa9SxabUdu5sj2bNoGjRQEPNl7KbIJapatMUzwX4NuWxNM5rCLyR4tDJwG3A6UDD0LGKwM5IU2ZFZB02W+owcCitD5BSoU4QAFOmWInwSy+FV1/1+YSuQNqwwQa6Z8+2x9q1drxCBVs7Gk4YLVt6l1Q0spsg/g3UASZhU10HAD+p6tWZCKAosAloq6rrUxx/BPhdVe+KcM46ICEzpT0KfYIAuOce+Ne/4IEH4B//CDoa52Ju40ZLFHPm2NeVK+14uXLWqO7c+cguvr7K+1jZTRACXAycETo0B3hbM7EVnYh0B25X1Q6prrsBOFNVV0U4Zx2eIDJPFQYOhDffhPfeg969g47IuVz1889HksWcObBsmR0vXRpOP+xnSjEAACAASURBVP1IC6NtWyhVKthY84LAtxwVkfFYeY6nUxw7A3g0zcBEfgR2YK2W51T1+TTeNxIYCXDSSSe1Wr9+faS3FS5799qfTCtXwpdf2miec4XUr7/aoHe4S+qbb+zvqBIlLEmEE8bpp1upkMImq7OYZgJvAu+p6uYUx4ths5qGAHNV9aUMbl4C2Aw0UdWtKY7/B1itqo+kcV5NVd0sIscD04FrVXVOevfyFkQKmzbZlNfSpeHrr6Fq1aAjci5P2LHDKtSEE8aiRTYBsFgx+18mnDA6dLD9Mwq6rCaIMsAVwKXY6unfgNJASWAG8IyqZvjbWETOB65R1e4pjhXDxiRaqerGKK5xB7BHVcem9z5PEKnMn2//0tu1s7rN3gHr3DF27z6SMObMsZniBw8eWbgXThgdOxbMld45sZK6JHA8sC+z+0GIyOvAxylbGiLSA7hFVTuncU5ZoIiq7g59Px24S1U/Su9eniAieO01GDQI/vpX20bMZzY5l669e61nNtzCmD/f1qOK2FTacMI444yC0TAPbAwi1Ar5CThZVX9Pcfxl4CtVfTbFsZrAi6raS0ROBqaEXioGTFTVezO6nyeINIwZY/tHPv00XHNN0NE4l6/s32+9tOGEMW8e7NtnrzVubF1R4ccpp+S/v8ECH6TOLZ4g0nD4sJXZnDbNdnzp1i3oiJzLtw4csMo24S6pefPg99Cfv9WrW3mQcMJo2TLv9+x6gnDW0Xr66bB5s7WZ69cPOiLnCoSkJPjuOxvHCD/Ci/dKlbL1F+GE0b49VK4cbLyp5cQYRG2gvqrODI1HFFPVP3I4zmzzBJGBtWuhTRvbdf6rr2zpqXMux23ZYi2LcMJYtMgKEAKcdtrR3VJBb6KU3YVyw4FRQAVVPUVEGgD/VtWzcj7U7PEEEYVZs+Dss+Gss+D99714jXO5YO9emx0VThjz5tmuwQDHH39st1TJkrkXW3YTxBKgDTA/xZajS1W1WY5Hmk2eIKL03HNw5ZVw440wNt2Zw865GEhKghUrju6WWrPGXitZ8thuqSpVYhdLdhPEV6raLrwndaiu0hJVjYtFsNnhCSITRo2CZ56Bl16CoUODjsa5Qu/nn4/tljp40F5r1Ojobqn69XOuWyq7CeIRYCswDLgauAZYpaq35Ex4OccTRCYcPAg9e1oNgpkz7c8U51yesW/fsd1SO3bYa1WrHt0tlZCQ9W6p7CaIolito+7YhkEfY7WRkrIWTux4gsik336zYjS7dtm/xJNOCjoi51wakpLg+++P7pZavdpeq1wZtm3L2gZKWU4QoeQwXlWHZP62uc8TRBasWGGlOOrVs39xhbFamXP51Nat1rL4+We46qqsXSO9BJFuvlHVw0ANESmetVu7PO+00+D11+Hbb2HIEPszxTmXL1Svbmtgs5ocMhJNg2Qt8LmI3CIi14UfsQnHBaJnT3joIXjnHbjrmL2bnHOFVDQb8m3DiuWVCT1cQfS3v9nOKnfeafta9+sXdETOuYBlmCBU9V+5EYgLmAg8+yz88IN1NVWoAN27Z3yec67AyjBBiMh0bFe3o6Tc38EVECVLwpQpttL63HNh/HgrFe6cK5Si6WK6NcX3pYCLgD9jE44L3PHHW4nKvn1h8GAr7nfzzfmvhrFzLtui6WKan+rQbBGZHaN4XF5QoQJ8+KGtsB4zxrYvfewxr9vkXCETTRdTyl1ZiwCtgBoxi8jlDSVLwoQJULMmPPqolad89VWrX+ycKxSi6WJajo1BCHAI+BEYEcugXB5RpAg88gjUqmWF/bZuhffeK5gb8zrnjhFNgjhZVQ+mPCAi0ZznCoq//c1aEkOG2M7tH30EJ54YdFTOuRiLZqFc6jEIgK8zOklEGorIkhSPXSIyWkTuEJFNKY73SuP8HiKyUkRWi8iYKOJ0sTRggCWGjRttZ7pvvw06IudcjKWZIETkeBFpDpQWkTgRaRZ6dCSKBXOqulJV41U1Hhu32AtMCb38WPg1VZ0W4d5FgWeAnkBjYKCINM78x3M5qmtXq/6qCp062eZDzrkCK72uonOB4UBt4N8pju8GMrt4rhuwRlXXS3TTJdsAq1V1LYCIvA6cD3yXyfu6nNasGXz5JfToAeecYwPXl1wSdFTOuRhIswWhqi+paifgclXtlOLRS1XfyuR9BgCTUjwfJSJLRWS8iEQa8awF/JTi+cbQsWOIyEgRSRSRxG3btmUyLJclJ50Ec+fatlcDBsATTwQdkXMuBjIcg1DVN0XkHBH5m4j8X/gR7Q1EpATQBwgnlf8ApwDxwBbgkUinRQoljfieV9UEVU2oVq1atGG57KpcGaZPhwsugNGj4aabvBKscwVMhglCRP4NDAH+BpQGLgNOzcQ9egKLVHUrgKpuVdXDoQ2HXsC6k1LbCKScJlMb2JyJe7rcULo0vPUWXH217W09aBAcOBB0VM65HBLNLKaOqvoXYHuocF9b7Bd2tAaSontJRFIususLLItwzgKgvojUC7VABgBTM3FPl1uKFoWnn4b77oOJE6FXL9uhzjmX70WTIPaHv4rICaHndaO5uIiUAc4GJqc4/JCIfCsiS4GuwA2h99YUkWkAqnoIGIVtb7oCeFNVl0dzTxcAEbjlFnj5ZZg9G844w1ZeO+fytWgWvE0TkYrAWGAJcBh4JZqLq+peoEqqYxHLg6rqZqBXiufTgGOmwLo8bMgQ2+Lq4ottrcRHH0GjRkFH5ZzLonRbECJSBPhQVXeGZi7VA+JUNepBalfI9OhhrYh9+6BDB9sw1zmXL2W0J3US8ESK5/tU9beYR+Xyt1atbK1E5crQrZvVb3LO5TvRjEFMF5HzYx6JK1hOPtlaD82awYUXwnPPBR2Rcy6TohmDGAVUEJE/gX3YGgVV1coxjczlf9WqwWefQf/+cOWVVsfprrt88yHn8oloEkTVmEfhCq6yZeHddy1B3HOPbT703HNQvHjQkTnnMhDNjnKHRWQAVvb7PhGpDVQHFsY8OlcwFCsGL7wAtWvDnXfCzz/Dm29CuXJBR+acS0c0K6mfxtYrhKen7gWejWVQrgASgTvusNbDxx9bZdhffgk6KudcOqIZpG6vqn8ltGAuNIupREyjcgXXyJEwZQosX27TYNesCToi51waokkQB0PrIRRARKoAXpXNZV2fPjBjBuzYYQvqEhODjsg5F0E0CeIZ4B2gmojcCcwFHoxpVK7gO/10+OILG8Tu0gU+/DDoiJxzqURT7vu/wK1YqY3fgH6q+nqsA3OFQMOGtlaifn3o3dtqOTnn8oxoWhAARYGDwIFMnONcxmrUsNIcXbvCsGFw7722palzLnDRzGL6J1auuyZW5nuiiNwS68BcIXLccfDBB3DZZXDrrXDNNXD4cNBROVfoRbNQ7jKgVagyKyJyL7YG4v5YBuYKmRIl4JVXoFYtePBBKxc+caJtSuScC0Q03UXrOTqRFAPWxiYcV6gVKQIPPABPPmkF/rp2hdWrg47KuUIrmgSxF1guIi+KyAvAt8BOEXlURB6NbXiuULr2WtvK9PvvrdjfY495l5NzAYimi+mD0CPsqxjF4twRF10E7dpZDae//c0SxvjxvgGRc7komlpM47JyYRFpCLyR4tDJwG1ALaA3NiNqDTBMVXdGOH8dsBvbwe6QqiZkJQ6Xj9WqBVOn2ljEdddBfLyV6/j7362+k3MupqKZxdRDRBaIyC8i8puI7BCRDDcNUtWVqhqvqvFAK6yragowHWiqqs2AH4D0ZkR1DV3Dk0NhJQKXXmqlOc491/a+btcOvv026MicK/CiGYN4Gvgr9pd/Naz8d7VM3qcbsEZV16vqJ6p6KHT8K2zqrHPpO+EEeOcd62rasMF2rbvzTjhwIOjInCuwokkQG4ElqnpQVQ+HH5m8zwBsLUVqw4G0aiwo8ImILBSRkWldWERGikiiiCRu27Ytk2G5fOfii+G776BfP+tuat0aFnrleediQTSDVasi0ga4HZgF/Bk+rqpPRnUDkRLAZqCJqm5NcfyfQAJwoUYIQkRqqupmETke65a6VlXnpHevhIQETfTCb4XH1Kk2iP3LL3DzzXDbbVCqVNBROZeviMjCtLrxo2lB3IkNFFfEupbCj2j1BBalSg5DgPOASyMlBwBV3Rz6+gs2dtEmE/d0hUGfPtaaGDwY7r8fWrSAL78MOirnCoxopoIcr6qtsnGPgaToXhKRHsA/gM7h1dmpiUhZoIiq7g593x24KxsxuIKqYkWb/tq/P4wYYXtMjB5t25uWKRN0dM7la9G0IGaIyJlZubiIlAHOBianOPw0UB6YLiJLROTZ0Htrisi00HuqA3NF5Bvga+ADVf0oKzG4QuKcc2DZMutyeuwxW2A3e3bQUTmXr0UzBrEDqIBNUz0ACKCqWjn24WWOj0E4AGbOhCuugLVrrfDf/fdD+fJBR+VcnpTdMYiqQHEsSWR1mqtzuadrV1i61Lqa/v1viIuD6dODjsq5fCeaDYMOA/2Af4S+rwHExzow57KlbFnravr8cyhZErp3tzGK338POjLn8o1oVlI/DXQFBoUO7QWejWVQzuWYDh1gyRL4xz9sMLtJE9t7wjmXoWi6mNqr6l+B/QCq+htQIqZROZeTSpe2MuJffQWVKsF558GgQbB9e9CROZenRZMgDopIEWxlMyJSBUiKaVTOxULr1pCYaAvqXn/dWhOTJ2d8nnOFVJoJQkTCaySeAd4BqonIncBc4MFciM25nFeypNVwWrAAata0suKXXGKrsZ1zR0mvBfE1gKr+F7gVGAvsAPqp6uu5EJtzsRMfD/Pnw7332u51jRtbWfEMpn07V5iklyAk/I2qLlfVJ1T1cVVdlgtxORd7xYvD//0fLF4Mp55qZcUvuAA2bw46MufyhPRKbVQTkb+l9aKq+najrmBo3Bi++AIefxxuvdWeP/YYDB1q+1E4V0il14IoCpTDymJEejhXcBQtCjfeaAvsmjWD4cOhZ0/be8K5Qiq9FsQWVfUCea5wqV8fZs2yFdhjxthMp3vugb/+1UuJu0InqjEI5wqVIkVg1Cjb1rRdOyvZUa8ePPoo/PFH0NE5l2vSSxDdci0K5/KievXgk0/gs89sXOLGG6FOHZv55CU7XCGQZoIIrZh2rnATseJ/M2bYQHbbtjaQXacO/Otf8OuvQUfoXMxEs5LaOQfQvr3VcVq4ELp1s7GJunXhppvg55+Djs65HOcJwrnMatkS3nnHNig6/3wbm6hXD669Fn76KejonMsxniCcy6omTWDCBPj+e/jLX+DZZ+GUU6ys+Jo1QUfnXLbFLEGISMPQlqLhxy4RGS0ilUVkuoisCn2tlMb5PURkpYisFpExsYrTuWyrXx/GjYPVqy05vPoqNGhgFWNXrAg6OueyLGYJQlVXqmq8qsYDrbB9JKYAY4AZqlofmBF6fhQRKYoVCewJNAYGikjjWMXqXI6oUweeeQZ+/BFuuMEqxTZpAv362Z4UzuUzudXF1A1Yo6rrgfOBV0LHXwEuiPD+NsBqVV2rqgeA10PnOZf31agBY8fC+vVW6+mTT6BFC+jd2woEOpdP5FaCGABMCn1fXVW3AIS+Hh/h/bWAlKN9G0PHnMs/qla1mU7r18Pdd8O8ebbw7uyzYfZsrxzr8ryYJwgRKQH0Ad7KzGkRjkX8v0lERopIoogkbtu2LSshOhdbFSva2on16+Hhh22FdpcucMYZ8PHHnihcnpUbLYiewCJV3Rp6vlVEagCEvkbaqWUjcGKK57WBiDWYVfV5VU1Q1YRq1arlYNjO5bBy5eDvf7cxiqeegnXroEcPaNPG9qRI8o0aXd6SGwliIEe6lwCmAkNC3w8B3otwzgKgvojUC7VABoTOcy7/K13aaj2tWQMvvAC//Wb7UMTHwxtvwOHDQUfoHBDjBCEiZYCzgZQb/z4AnC0iq0KvPRB6b00RmQagqoeAUcDHwArgTVVdHstYnct1JUrAFVfAypU2NfbgQRgwwOo+vfKKPXcuQKIFqP8zISFBExMTgw7DuaxJSrKpsffcA998Y2U8xoyxjYtKlgw6OldAichCVU2I9JqvpHYuryhSBC6+2LZA/d//oHp1uPJKOPlkeOIJLzXucp0nCOfyGhE47zz48kuYPt1Wao8ebQlj0CCb+XToUNBRukLAE4RzeZUInHWW7XD3xRdW7+n9923mU61acP31sGCBT5N1MeMJwrn8oH17eP55Kys+eTJ06mTFAdu0gYYN4a67vECgy3GeIJzLT0qWhL594e23YetWePFFqF0b7rgDTj3VVmo/9RT8Eml5kXOZ4wnCufyqYkW4/HLbEnXDBnjoIdi/H667DmrWhF69YOJEH9x2WeYJwrmCoHZt29luyRIr5XHTTbah0aWX2uD2ZZfBhx/64LbLFE8QzhU0TZvC/fdbKY/Zsy1JTJtmLYpatayF8fXXPrjtMuQJwrmCqkgRKwj43HOwZQtMmWLPn38e2ra1TY3uuANWrQo6UpdHeYJwrjAoWdLqPb31lg1ujxsHJ51ks58aNLCE8eST9ppzIZ4gnCtsKlSA4cNhxgz46ScrQX7ggK2rqFULeva0vbb37Ak6UhcwTxDOFWa1alkJ8sWLbVD75pttH+3LLrPB7fD4hRcOLJQ8QTjnTJMmcN99sHYtzJljZT0+/BDOPdcSybXXwldf+eB2IeLVXJ1zaTtwwJLEhAkwdSr8+SccfzyceaaVAenWzarOunwrvWquniCcc9H5/Xfb+W76dPj0Uyv7AVZttls3Sxhdu4Lv7JiveIJwzuUsVRurmDHDksWsWbBrl73WvPmR1kWnTrbVqsuzPEE452Lr0CFYuNCSxYwZVn32wAEoVszqQ4UTRtu2ULx40NG6FDxBOOdy1759liTCCWPhQmt1lC0LnTtbsujWDeLibEGfC0x6CaJYjG9cEXgRaAooMBwYDTQMvaUisFNV4yOcuw7YDRwGDqX1AZxzeVDp0tZqOOsse75jh3VDhRPGtGl2vFo1G/AOJ4yTTw4sZHesmCYI4AngI1W9WERKAGVUtX/4RRF5BPg9nfO7quqvMY7RORdrlSpZmfK+fe35xo1WhTacMN54w47Xq3ckWZx5ps2YcoGJWReTiBwHfAOcrBFuIiICbADOVNVjisGEWhAJmUkQ3sXkXD6kCitXHkkWM2fajCmAZs2OzJA64wwf8I6BQMYgRCQeeB74DmgOLASuV9U/Qq+fATyaZmAiPwI7sK6p51T1+TTeNxIYCXDSSSe1Wr9+fU5/FOdcbjp8GBYtOpIw5s619RfFitkgd7iF0aqVjWm4bAkqQSQAXwEdVHW+iDwB7FLVf4Ve/w+wWlUfSeP8mqq6WUSOB6YD16rqnPTu6S0I5wqgfftg3jxLFjNmQGIiJCXZnt0NG0LLltCixZGvlSoFHXG+EtQg9UZgo6rODz1/GxgTCqgYcCHQKq2TVXVz6OsvIjIFaAOkmyCccwVQ6dJHWg0AO3daq2LRInt8/rntnBdWt+6xSaNGjUBCz+9iliBU9WcR+UlEGqrqSqAb1t0EcBbwvapujHSuiJQFiqjq7tD33YG7YhWrcy4fqVgRzjvPHmG//moFBxcvPpI4Jk8+8voJJxydMFq2tEQikuvh5yexnsV0LTAhNINpLTAsdHwAMCnlG0WkJvCiqvYCqgNTbBybYsBEVf0oxrE65/KrqlXh7LPtEbZrF3zzzdFJ45NPbIwDLNGkThoNGkDRosF8hjzIF8o55wqP/fttz+5w0li82JLIn3/a62XKWKmQlEmjSRMoUSLYuGPIV1I751xaDh6E778/OmksXgy7d9vrxYvbPt8pWxvNmxeYGVSeIJxzLjOSkmDNmqOTxqJFNtYBR2ZQxcdbt1T9+nDqqfa1cuV8NbbhCcI557JL1VaAp0waS5fC+vVHb6JUseLRCSP8NY8mD08QzjkXK3/+CT/+CKtWwerVR3+NNnmceipUqRJI8gisWJ9zzhV4JUtCo0b2SC2t5PHll1Z/KinpyHsrVozc6ggweXiCcM65WIkmeaRudXz1VXTJI/w1hsnDE4RzzgUho+Sxbt2xLY9IyaNCBStqOHt2jicKTxDOOZfXlCxps6QaNjz2tQMHju22OnAgJq0ITxDOOZeflCiRdvLIYb7Xn3POuYg8QTjnnIvIE4RzzrmIPEE455yLyBOEc865iDxBOOeci8gThHPOuYg8QTjnnIuoQFVzFZFtwPosnl4V+DUHw8kP/DMXfIXt84J/5syqo6rVIr1QoBJEdohIYlolbwsq/8wFX2H7vOCfOSd5F5NzzrmIPEE455yLyBPEEc8HHUAA/DMXfIXt84J/5hzjYxDOOeci8haEc865iDxBOOeci6jQJwgR6SEiK0VktYiMCTqeWBORE0VkpoisEJHlInJ90DHlFhEpKiKLReT9oGPJDSJSUUTeFpHvQ/+9Tw86plgTkRtC/66XicgkESkVdEw5TUTGi8gvIrIsxbHKIjJdRFaFvlbKiXsV6gQhIkWBZ4CeQGNgoIg0DjaqmDsE3KiqpwHtgGsKwWcOux5YEXQQuegJ4CNVbQQ0p4B/dhGpBVwHJKhqU6AoMCDYqGLiZaBHqmNjgBmqWh+YEXqebYU6QQBtgNWqulZVDwCvA+cHHFNMqeoWVV0U+n439kujVrBRxZ6I1AbOBV4MOpbcICLHAWcA4wBU9YCq7gw2qlxRDCgtIsWAMsDmgOPJcao6B/gt1eHzgVdC378CXJAT9yrsCaIW8FOK5xspBL8sw0SkLtACmB9sJLniceBmICnoQHLJycA24KVQt9qLIlI26KBiSVU3AWOBDcAW4HdV/STYqHJNdVXdAvZHIHB8Tly0sCcIiXCsUMz7FZFywDvAaFXdFXQ8sSQi5wG/qOrCoGPJRcWAlsB/VLUF8Ac51O2QV4X63c8H6gE1gbIiclmwUeVvhT1BbAROTPG8NgWwSZqaiBTHksMEVZ0cdDy5oAPQR0TWYd2IZ4rIa8GGFHMbgY2qGm4dvo0ljILsLOBHVd2mqgeByUD7gGPKLVtFpAZA6OsvOXHRwp4gFgD1RaSeiJTABrSmBhxTTImIYP3SK1T10aDjyQ2qeouq1lbVuth/489UtUD/ZamqPwM/iUjD0KFuwHcBhpQbNgDtRKRM6N95Nwr4wHwKU4Ehoe+HAO/lxEWL5cRF8itVPSQio4CPsRkP41V1ecBhxVoHYBDwrYgsCR37P1WdFmBMLjauBSaE/vhZCwwLOJ6YUtX5IvI2sAibrbeYAlh2Q0QmAV2AqiKyEbgdeAB4U0QuxxJlvxy5l5facM45F0lh72JyzjmXBk8QzjnnIvIE4ZxzLiJPEM455yLyBOGccy4iTxAuTZGqRoaO9wtVzEwSkXQ3Sg9V19wvIhXSeY+KyKspnhcTkW05UXU1miqXIlI3wme8Q0T+nt37Z4eIXCkig4OMIael/FmLSIKIPBl0TC5tniBcel7m2KqRAMuAC4E5UVxjILYgsW867/kDaCoipUPPzwY2RR9mumJS5TJaoaJxWaKqz6rqf3MynrRkJ84Mrls0rddUNVFVr4vFfV3O8ATh0pRG1UhUdYWqrszofBE5BSgH3IolivR8iFVbJfTeSZmLNk3ZrnIpIvEi8pWILBWRKeFWiIjMCregRKRqqJQHIjJURN4Skf8Bn4hIDRGZIyJLQvsUdIpwjwdE5LvQPcaGjiW3YkL3elBEvhaRH8LXCO1xMVZEvg2de23oeCsRmS0iC0Xk43AZhlT3fFlEHhWRmcCDIlI21GpcECrwd34G9+gWet+3ofNKho6vE5HbRGQu0C8Uyzci8iVwTYr7dwm3EkOfdXzoc64VketSvO9fYntaTBfb4yHQll1hUqhXUruYC/+i/xxoKCLHq2paNWJeB24L/cJoBowHIv0iLR+6XiR/UdXU5SSOqnIpImlVuTwlxcpygBOwyqAA/wWuVdXZInIXtnJ1dBrXCTsdaKaqv4nIjcDHqnpv6C/qMqk+U2WshdVIVVVEKqZxzWKq2kZEeoViOAsYiRWnaxGqDFBZrNbWU8D5qrpNRPoD9wLDI1yzAXCWqh4WkfuwMiTDQzF8LSKfAoMj3KMU1sLspqo/iMh/gauwqrkA+1W1Y+jzLU3x83s4nZ9ZI6ArUB5YKSL/wfaxuAirOlwMWyVdmIouBsoThIulAUBfVU0SkcnY8v9nIr1RVZeKlR8fCKRZ9iO0h0V8zofKGlVNvq6I3BH6WgGoqKqzQy+9ArwVxfWmq2q49bUAGB/6xf2uqi5J9d5dwH7gRRH5AEhr7CVcWHEhUDf0/VnAs6p6CCCUkJoCTYHpIgJWRmZLGtd8S1UPh77vjhU1DP+FXgo4KY17NMcK4/0Qeu8rWOsgnCDegIg/v1exDboi+UBV/wT+FJFfgOpAR+A9Vd0Xut7/0jjXxYAnCBcTItIMqM+RX1LhekARE0TIVOyv9i5AlTSum9kWxFYRqRFqPeRYlcuQQxzppk29teUf4W9UdY6InIF1ob0qIg+nHFsI/VXeBisuNwAYBZwZ4X5/hr4e5sj/u8KxJeoFWK6q0Wwx+keK7wW4KHX3odh/wEj3iOa6kc5Ny58pvg9/xozu42LIxyBcrAwE7lDVuqFHTaCWiNRJ55zxwF2q+m1ab1DV3aoan8YjUrXSbFW5VNXfgR0pxg0GAeG/htcBrULfX5zWNUKf+RdVfQGrpNsy1evlgAqhgomjyVwL6RPgSgkNMoe6q1YC1SS0B7WIFBeRJlFc62Pg2lBCQERapHOP74G6InJq6D0pfy7JQrvY/S4iHUOHLs3EZwOYC/QWkVKhn9O5GZ3gco4nCJcmsaqRX2LjBxvFKkUiIn3FqkieDnwgIh9HOH0AMCXVsSnAABGpKSLHdCOp6kZVfSJCHAkiktWtQh8AzhaRVdjsqAeycI0hwMOhvvR44K7Q8bHAVSIyD6iazvldgCUishjrT0/9GcsDYXxV8AAAAKtJREFU74euPxu4IROxvYhV71wqIt9gragDWMJ6MHRsCdHti3A3UDx0rWWh52ndYz9WHfYtEfkW26nv2TSuOwx4JjRIvS8Tnw1VXYAl+W+wLrZE4PfMXMNlnVdzdc7laSJSTlX3iEgZbGr1yPC+6i62fAzCOZfXPS8ijbFxnlc8OeQeb0E455yLyMcgnHPOReQJwjnnXESeIJxzzkXkCcI551xEniCcc85F9P9fXn03QH6FugAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "# 3a\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "N=10 \n", | |
| "Ti=85 \n", | |
| "Ta=65\n", | |
| "K=0.275\n", | |
| "t=np.linspace(0,10,N)\n", | |
| "T_analytical=np.zeros(len(t))\n", | |
| " \n", | |
| "for i in range(0, len(t)):\n", | |
| " T_analytical[i]=Ta+(Ti-Ta)*np.exp(-K*t[i])\n", | |
| "Tn=np.zeros(N) #spaceholder\n", | |
| "Tn[0]=Ti\n", | |
| "for i in range(1,len(t)):\n", | |
| " Tn[i]=Tn[i-1]-K*(Tn[i-1]-Ta)\n", | |
| " \n", | |
| "plt.plot(t,T_analytical,label=\"Analytical\",color=\"red\")\n", | |
| "plt.plot(t,T_E,'-',label=\"Numerical Approximation\",color=\"blue\")\n", | |
| "plt.xlabel('11 A.M.= 0 Hours since recording')\n", | |
| "plt.ylabel('Temperature (degF)')\n", | |
| "plt.legend(loc='best')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 183, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Death was 1.8865228851460631 hours ago\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "T=98.6\n", | |
| "K=0.275\n", | |
| "Ti=85\n", | |
| "Ta=65\n", | |
| "dT=np.log((T-Ta)/(Ti-Ta))/-K #change in time\n", | |
| "print(\"Death was\", -dT, \"hours ago\")" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## 3b. \n", | |
| "As t approaches infinity, the temperature will reach to the ambient temperature of 65 degrees F" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## 3c. \n", | |
| "Since the first temperature recorded on the body was at 11 am, subtracting 1.88 hours yields the time of death to be at 9:07 A.M. " | |
| ] | |
| }, | |
| { | |
| "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": 42, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "Text(0, 0.5, 'Temperature(F)')" | |
| ] | |
| }, | |
| "execution_count": 42, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZgU5dX+8e8NqIDiEsG4IlER44IEkRd3cSXRX4xG3OIrKkpMoibuC3GPeU00GhONirvRuIsm7ijuIopiFEWjQFQEFQQVlUXg/P54apx2nBl6YGpqZvr+XFdfdFV1VZ3qHk5XP/XUeRQRmJlZ5WhTdABmZta0nPjNzCqME7+ZWYVx4jczqzBO/GZmFcaJ38yswjjxm7VikiZI2qKJ99lG0k2SPpH0ZAPXXUbS55JWzys+c+JvtrI//qrHQkmzS6Z/VnR8S0LSB5K2LjqOPEnaQFLU+Byfz3mft0j6bem8iFg3Ikblud9a7AhsAawWEduWLpB0Vsn7MUfS/JLpFyNibkQsFxFTmjjmitKu6ACsdhGxXNVzSf8FDouIR4qLqDyS2kXE/Ja+j0ayoPRzrCBrAxMjYnbNBRFxBnAGgKQjgL0jYqcmjq/i+Yy/hZLUVtJpkiZKmp79tF4xW7ZBdiY1WNL7kj6WdKikLSSNy36CX1iyrSMkjZR0haTPJL0uaduS5d+RdEN2pv6epDMktamx7qWSZgInZ/t/XNIMSdMkXS+pU/b624FVgIezs7yjJQ2Q9HaN4/v6V4Gk8yT9Q9KtkmYB+9V3/LW8V10kPZDFMkPSPZJWK1l+uKT/SpqVbW9gHdvZStJoSZ9KmiLpIkkNPnnKjueqkukNJM0vmX4ue4+fyz6P+yWtVLJ8+2zZp5LelXSApKOBnwKnZe/r7bW8jx2yz2mqpMmSzpe0VLZsgKS3JZ2avU/v1/fLUlLXLK4Zkv4jaVA2/5fAJcD2WRynNvC9aa/0S2nNbPoWSRdLGiHpi+zvahVJf8v+jl+TtEnJ+mtln+/07LM8oiH7rxgR4UczfwD/BXaqMe9k4ClgdaA9cB1wbbZsAyCAi4FlgB8DXwB3AisDXYGZwP9krz8CmA/8ElgKOAiYASyfLX8A+CvQEVgNGAsMqrHu4UBboEO2/x2ApYFVgeeA80pi/wDYumR6APB2jeP7+jXAecBc4Eekk5UO9R1/Le/fd4E9svVWAO4BbsmWrQR8AqybTa8BfL+O7fQFNs+Oc13gbeCIOl67ATC/jmXnAVfV9drs/Xoz28eywLPAmdmy9YDPSUm+HdAF2DRbdgvw23rexz9m71nn7D15ARha8hl8BQzN/gb2BGYBy9VxDKOBi7K/rz7Z38tWJX8Tj5Txd/2t12WfZQBrlhzTB8Cm2ef3NDAR2Df7HM4HHshe2xZ4FTgp+9tbH3gX2K7o/8PN7VF4AH6U8SHVnvgnVf1Hy6a/B3wJiOrEv3LJ8i+APUqm76tKWtl/wEk1tv8KMJD0s/0LYKmSZYeU/Gc7AvjPIuLfDxhVMr04if/hco+/jPezHzA1e16V+PcA2jfwczkZuLmOZVWfwScljyNLjmdRif/4kuljgbuz52fVs89FJf73gR1Klu0BvFHyGXwKtClZ/hnQq5b9dAfmAB1K5l0EXF7yN9GYif+vJctPAMaWTG8OfJA93w54q8b2zgIua4z/h63p4Tb+FkiSgLWA+yWVVtlrQzqjh9S+/HHJstnAhzWmS9ufJ9fYzTuks+m1Sf8Zp6Xdfr2f0qaZ92rEtzrp18aWQKfs9VPLObZ6fL2PMo5/eo14OmXx7ARUNQd1AIiImVmTxrHA9Uq9UI6NiG80PWXb2RD4E9A7W78d8Ew9MS+IiFqbn8rwQcnzL6n+rNYCJjR0Y9l7tirpc63yDukXTpVpEbGwjv2WWj17bWkb/juki7p5qPl3W9ff8dpAN0mflCxvCzT7a2NNzW38LVCkU5mqs7cVSx7tI2L6otavw5o1prsCU0gJ93NgpZL9LB8RvUtDqrHu+aRfCRtHxPLAYaRfInW9/gtSMxIAWbvzd2q85ut1FuP4T86Ob/Msnl1K44mI+yJiR1JCexe4rJZtAFwJvERqFloeOLvGcZXrG8dLSsjleo/UBFSbOkvtZu/ZB6TkWKUr6X1sqClAF0kdGmFbjek90i+Y0r+JThGxZ8FxNTtO/C3X5cB5ktYCyC54/b8l2N5a2YXadpIOJP1HfjgiJpGaHv4oqZNSH+3uqr87ZifSl8VnkrqSzqZLfQisUzI9HviOpB2zpH8Wi/7bbMjxdyKdvX4iqTPwdZdHSWtI2k1SR9J1hM+BBfVs59OI+FzSRqTrGovjZaB/tu+VSG3S5boB2F3SntkF7i6SembLar6vNd0MnCFpZUmrkNrzb1yM+N8mNQX+TqnffW9gEHDTYmyrMT0NIOk32UXidpJ6ZvFZCSf+luuPpJ+wI5V6ujxLaoJYXE8CPyBdpBsK7BkRn2bL9ic1kbyRLb+VdHGwLqcDW5PajIeTLiqXOhc4N+uVcWR2lv5rUuKYTDozXdQvl4Yc/wWkC5ofk5LD/SXL2gKnZPv8mNRmfFQd2zkGOEzS58ClpPdhcdwH3Au8TvpSvbvcFSNiAqlt/lTSBfoxwEbZ4mHA5tn7ekstq5+e7fM10pfPM6T3sUGyXw/7ABuS3rdbgRMi4qmGbqsxRcRXpA4AW5KanqaRfr1VYpfaeim7AGIVTO5PbVZRfMZvZlZhnPjNzCqMm3rMzCqMz/jNzCpMi7iBq3PnztGtW7eiwzAza1FefPHF6RHRpeb8FpH4u3XrxpgxY4oOw8ysRZH0Tm3z3dRjZlZhnPjNzCqME7+ZWYVx4jczqzBO/GZmFcaJ38yswjjxm5lVGCd+M7NmZu5cGDkSTjkFpkxp/O23iBu4zMxaswgYNw5GjICHH4Ynn4TZs6FdO9hqK1h99cbdnxO/mVkBpk6FRx5JyX7ECPggG2W5Rw847DDYeWfYfnvo1Knx9+3Eb2bWBL78Mp3JV53VjxuX5q+8ckryVY+11so/Fid+M7McLFwIY8dWn9E//TTMmwdLLw3bbAMHHpgSfa9e0KaJr7Y68ZuZNZJ3361O9I88Ah9/nOZvsgkcdVRK9NtsAx07FhunE7+Z2WL67DN4/PHq5pv//CfNX3VV2G23lOh32ilNNydO/GZmZZo/H154ofqs/rnn0rwOHWC77eCII1Ky32gjkIqOtm5O/GZmdYiACROqE/3IkfDppymp9+4NJ5yQEv2WW8IyyxQdbfmc+M3MSsyYkRJ8VbKfNCnN79oVBg5MiX6HHaBz52LjXBJO/GZW0ebNg1GjqhP9mDGpR06nTtC/Pxx3XEr23bs37+abhnDiN7OKEgHjx1cn+scfhy++gLZtoW9fOO20lOj79oWllio62nw48ZtZq/fRR9+8S/b999P89daDQYNSou/fH1ZYodg4m4oTv5m1OrNnpxumqhL9yy+n+SutlLpXVt0l261boWEWxonfzFq8hQvhlVeqE/1TT8GcOampZqut4NxzU6Lv3Ts16VS6XBO/pBWBq4CNgQAOjYhRko4CjgTmA/dFxIl5xmFmrc/773/zLtmPPkrzN9ywuj/9ttvCcssVG2dzlPcZ/8XAgxGxt6SlgY6S+gN7AD0jYq6kVXKOwcxagc8/hyeeqE72r7+e5q+ySnXTzU47wRprFBtnS5Bb4pe0PLAtcDBARMwD5kn6BXBeRMzN5n+UVwxm1nItWAAvvlid6J99Fr76Ctq3T/VuDjkkJftNNmn6ImctXZ5n/OsA04BrJW0KvAj8Glgf2EbSucAc4PiIeKHmypKGAEMAunbtmmOYZtZcTJpUnegffRRmzkzze/WC3/wGdtkFtt46JX9bfHkm/nZAb+CoiBgt6WLg5Gz+SkA/YHPgNknrRESUrhwRw4BhAH369PnGMjNrHT75BB57rDrZv/12mr/GGvCTn6Qz+h13TM051njyTPyTgckRMTqbvoOU+CcDd2WJ/nlJC4HOpF8HZtaKffUVjB5dneiffz416Sy7bBptqqp08QYbtJ67ZJuj3BJ/RHwg6T1JPSLiTWBH4HVgArAD8Lik9YGlgel5xWFmxYmAt95KJYtHjEhn97NmpTb5Pn3SYOI77wz9+qUBSqxp5N2r5yjgpqxHz0TgEOAL4BpJ44B5wKCazTxm1nJNn57a56vO6t99N83/3vfggAOqi5yttFKxcVayXBN/RLwM9Kll0YF57tfMms7cufDMM9WJ/qWX0pn+Ciuk9vmqs/p11y06UqviO3fNrEEi4LXXqptvnnwyDSTerh1ssQWcdVZK9H36pHnW/PhjMbNF+uCDb94lO3Vqmt+jBwwenBL99tunUsbW/Dnxm9m3fPllqndTdVb/6qtp/sorV98lu/POsNZaxcZpi8eJ38xYuDBVsKxK9E8/nQYoWXrpdJfseeelRN+rl++SbQ2c+M0q1LvvfvMu2elZp+pNNqnuT7/NNtCxY7FxWuNz4jerEJ99lkabqkr2b76Z5q+2GvzoR9VFzlZdtdAwrQk48Zu1UvPnwwsvVCf6555L8zp0SBdif/7zlOw32sh3yVYaJ36zVmTChOp2+pEj4dNPU1LfbDM44YSU6LfcEpZZpuhIrUhO/GYt2IwZKcFXndVPmpTmd+0KAwdWFzlbeeVi47TmxYnfrIWJgL/+FW66CcaMST1yOnVKZRCOOy4l++7d3XxjdXPiN2thrrsOfv3rdGfsaaelRN+3bxpf1qwcTvxmLcjrr8ORR6az+4cf9sDhtnh8K4ZZCzF7Nuy7b6pdf+ONTvq2+HzGb9ZC/OY3MG4cPPhg6ntvtrh8xm/WAtx2GwwbBiedBLvuWnQ01tI58Zs1cxMnwuGHp5LH55xTdDTWGjjxmzVj8+aldv02beDmm91zxxqH2/jNmrGTT0599e+8E9Zeu+horLXwGb9ZM3XvvXDRRan75l57FR2NtSZO/GbN0OTJMGhQqn9//vlFR2OtjRO/WTMzfz4ccEAaxPzWW6F9+6IjstbGbfxmzczZZ6dhD//+d1h//aKjsdbIZ/xmzcijj8LvfgcHHwwHHlh0NNZaOfGbNRMffpiSfY8ecMklRUdjrVmuiV/SipLukPSGpPGStihZdrykkNQ5zxjMWoKFC+Ggg+CTT9JdussuW3RE1prl3cZ/MfBgROwtaWmgI4CktYCdgXdz3r9Zi3D++ana5uWXp8HOzfKU2xm/pOWBbYGrASJiXkR8ki2+CDgRiLz2b9ZSPPssDB2aRswaMqToaKwS5NnUsw4wDbhW0lhJV0laVtKPgfcj4t/1rSxpiKQxksZMmzYtxzDNijNjBuy/fxoq8corPWqWNY08E387oDdwWUT8APgCOBMYCpy+qJUjYlhE9ImIPl26dMkxTLNiRMDgwTB1auqvv8IKRUdklSLPxD8ZmBwRo7PpO0hfBN8D/i3pv8CawEuSVs0xDrNm6dJL4e674bzzYPPNi47GKkluiT8iPgDek9Qjm7Uj8FJErBIR3SKiG+nLoXf2WrOKMXZsGhh9t93gmGOKjsYqTd69eo4Cbsp69EwEDsl5f2bN3qxZqdRyly5p4HS361tTyzXxR8TLQJ96lnfLc/9mzU0E/OIXMGECPPYYdPZdLFYA1+oxa0LXXQc33ZTq8Wy7bdHRWKVyyQazJjJ+fKqtv8MOcOqpRUdjlcyJ36wJzJ4N++yTSjHceCO0bVt0RFbJ3NRj1gSOOQbGjYMHHoDVVis6Gqt0PuM3y9ltt8EVV8CJJ8KAAUVHY+bEb5ariRPh8MOhX79UZ9+sOVhkU4+kXsA2wOrAbGAc8GhEfJpzbGYt2rx5qb9+mzZwyy2w1FJFR2SW1HnGL+lASS8CZwErAe8AnwE7AY9LulrSmk0TplnLc8opMGYMXH01rL120dGYVavvjH9lYNuI+KK2hZL6AN8nlV0wsxL33gsXXgi/+hXstVfR0Zh9U52JPyIurm/FiBjT+OGYtXyTJ8OgQdCrF1xwQdHRmH1bfU09D5Q8P7FpwjFr2ebPhwMOgLlzU6nl9u2Ljsjs2+rr1VNaKnm/vAMxaw3OPhueeioNobj++kVHY1a7+hK/h0U0a4CRI1OXzYMPhgMPLDoas7rVd3F3HUl3ASp5/rWI8CUrs8xHH8HPfgY9esAllxQdjVn96kv8Py157j9lszosXAj/+78wcyY89FCqx2PWnNXXq+fRpgzErKU6/3x4+OHUrt+zZ9HRmC1afb167pb0Q0nf+nKQtLak0yUdmm94Zs3bqFEwdCgMHAhDhhQdjVl56mvq+RVwHHCppA+BaUB7YB3gXeDSiLgz/xDNmqeZM2G//aBrV7jySg+haC1HfU097wPHAsdKWg9YjVSr582ImNVE8Zk1SxEweDBMmQLPPAMrrFB0RGblK7c65xygXXa37jxJvnxlFe3SS2H4cDjvPOjbt+hozBpmkYk/a8f/J3BVNmtt4J48gzJrzsaOheOOg912SwOsmLU05ZzxHw30I1XmJCL+A6ySZ1BmzdWsWanUcpcuaeD0Nh7RwlqgcoZenBMR85RduZLUlnRTl1lFiYBf/hImTIDHHoPOnYuOyGzxlHO+8kxWpK29pP7ArcC95Wxc0oqS7pD0hqTxkraQdH42/Yqk4ZJWXJIDMGsq11+fBko/4wzYdtuiozFbfOUk/hOBWcAbwK+BR4GhZW7/YuDBiNgA2BQYD4wANo6InsB/gFMaGrRZUxs/PtXW798/9ds3a8nqberJmnWuiYhBwGUN2bCk5YFtgYMBImIeMA94uORlzwF7N2S7Zk1t9mzYZ59UiuHGG6Ft26IjMlsy9Z7xR8QCYDVJizNa6Dqkm76ulTRW0lW1dAM9FHjg26uaNR/HHAPjxsENN8DqqxcdjdmSK6epZyLwlKRTJB1d9ShjvXZAb+CyiPgB8AVwctVCSUOB+cBNta0saYikMZLGTJs2rYzdmTW+22+HK66AE0+EAQOKjsascZST+KeR2uU7Al1KHosyGZgcEaOz6TtIXwRIGgTsDvwsImqt+x8RwyKiT0T06dKlnN2ZNa6JE+Gww6Bfv1Rn36y1WGR3zog4bXE2HBEfSHpPUo+IeBPYEXhd0gDgJGC7iPhycbZtlrd581J//TZt4OabYanFaew0a6YWmfgljaCW0bgiYpcytn8UcJOkpUlNRocALwDLACOyewOei4gjGhK0Wd5OOQXGjIE774Ru3YqOxqxxlXMD129LnrcnDdAyt5yNR8TLQJ8as9crLzSzYtx3H1x4Yeq+uZfHmbNWqJymntE1Zj0h6Ymc4jEr1OTJMGgQ9OoFF1xQdDRm+SinqWf5ksk2wGakEs1mrcr8+XDAATBnDtx6K7RvX3REZvkop6nnNVIbv0jdLycBh+cZlFkRzj4bnnoq9ddff/2iozHLTzmJf52I+Kp0Rm3DMZq1ZCNHpi6bBx+cBk43a83K6cdfs40f4PnGDsSsKB99BD/7GfToAZdcUnQ0Zvmr88xd0iqktvwOkjahuhTz8qSbucxavIUL4aCD0vi5Dz2U6vGYtXb1NdnsRqqlsybwt5L5s4DFuqnLrLm54IKU8C+7DHr2LDoas6ZR32Dr15IKrO0TEbc1YUxmTWLUKDj1VNh7b/j5z4uOxqzplNOP/zZJuwIbkW7gqpr/+zwDM8vTzJmw337QtStceSXIY8pZBSmnH//fgBVJtfWvJd25+1zOcZnlJgIGD4YpU+CZZ2BFjwFnFaacXj1bR8QBwMdZwbb/IbX7m7VIf/sbDB8O550HffsWHY1Z0ysn8c+p+lfSqtl0t9wiMsvRyy/DscfCj36UBlgxq0Tl3Ih1fzYg+gXAy8AC4PpcozLLwaxZaQjFzp3TwOltyjntMWuFFjXmbhvggYj4BLhd0r1Ah4iY0STRmTWSCPjlL2HChHSXbufORUdkVpxFjbm7ELi4ZHq2k761RNdfnwZKP+MM2G67oqMxK1Y5P3ZHSNoj90jMcjJ+fKqt378/DB1adDRmxSunjf9IYAVJc4HZpNINERHfyTUys0Ywe3YaQnHZZdMZf9u2RUdkVrxyEr9bQ63FOuYYePVVeOABWH31oqMxax4W2dQTEQuAgcBJ2fPVgF55B2a2pG6/Ha64Ak48EQYMKDoas+ZjkYlf0iVAf6CqSvmXwOV5BmW2pCZOhMMOg379Up19M6tWTlPPlhHRW9JYgIiYIWnpnOMyW2zz5qU6PG3awM03w1JLFR2RWfNSTuL/KuvPHwCSVgYW5hqV2RI49VR44QW4807o1q3oaMyan3K6c14K3Al0kXQW8DTwh1yjMltM990Hf/pTullrr72KjsaseSqnLPMNkl4EdspmDYyIcfmGZdZwkyfDoEGw6aYp+ZtZ7cqtVtIW+AqY14B1kLSipDskvSFpvKQtJH1H0ghJb2X/rrQ4gZuVmj8/jZs7Zw7cdhu0b7/odcwqVTm9eoYCNwOrk8ox/0PSKWVu/2LgwYjYANgUGA+cDDwaEd2BR7NpsyVyzjnw5JNpCMX11y86GrPmTRFR/wuk8cBmEfFlNt0ReDEivr+I9ZYH/g2sEyU7kfQmsH1ETJW0GvB4RPSob1t9+vSJMWPGlHVAVnkeewx23DENmn7ddUVHY9Z8SHoxIvrUnF9Os807fPNaQDtgYhnrrQNMI43bO1bSVZKWBb4bEVMBsn9XKWNbZt8yc2a6OeuHP0xn+ZdcUnREZi1DOYn/S+C1LHFfCbwKfCLpQkkX1rNeO6A3cFlE/AD4ggY060gaImmMpDHTpk0rdzWrAHPmwPnnwzrrwAUXpD77jz4Kyy1XdGRmLUM5/fjvyx5Vyh1vdzIwOSJGZ9N3kBL/h5JWK2nq+ai2lSNiGDAMUlNPmfu0VmzBglRo7bTT4L330iha//d/0LNn0ZGZtSzldOe8enE2HBEfSHpPUo+IeBPYEXg9ewwCzsv+vWdxtm+VIyIVWTv55FRwbfPN4YYbYPvti47MrGVaZOKXNAA4B1g7e31DyjIfBdyUlXiYCBxCal66TdJg4F1SATizWj3/fGrHf+IJWG+91FVz771BKjoys5arnKaeS4B9SG37DSrVEBEvA9+6okw6+zer01tvpdILd9wBq6wCl14Khx/uujtmjaGcxD8ZeDkbhtEsVx9+CGedBVdeCcssA2eeCcceC506FR2ZWetRTuI/EfiXpMeBuVUzI+IveQVllWfWrNRD509/grlzYcgQOP10+O53i47MrPUpJ/GfRSrXsCKuymmNbN68dHZ/9tnw0UcwcCCcey507150ZGatVzmJf5WI2Cz3SKyiRKQRsk49FSZMSD10/vUv6Nu36MjMWr9ybuB6VNIOuUdiFeOxx1KC33df6NgR7r8fRo500jdrKuUk/sOBRyR9LmmGpJmSZuQdmLU+r7ySyivssEO6iHv99TB2bJrn7plmTaecpp7OuUdhrdo776S7bW+8EVZcMZVbOPJIl042K8oiz/gjYgHpJquTsuerAb3yDsxavo8/huOPTwXUbrsNTjghtecff7yTvlmRyqnHfwnQH/jfbNaXwOV5BmUt2+zZ8Ic/wLrrwkUXpQFS3norzVvJw+6YFa6cpp4tI6K3pLEAETEjK8Fg9g0LFqR2+9NPh/ffh913T0XUNt646MjMrFQ5F3e/ktQGCABJK+P+/FYiInXF7NkTBg+GNddMtXX+9S8nfbPmqM7EL6nq18ClwJ1AF0lnAU8Df2iC2KwFGDUKttsOfvzjNO7tHXekedtuW3RkZlaX+pp6ngd6R8QNkl4EdiJV5hwYEeOaJDprtt58M918ddddqazCZZels30XUTNr/upL/F/3rI6I14DX8g/HmrupU1MRtauugg4dUqmFY47x6FdmLUl9ib+LpGPrWhgR9Q27aK3MZ5+l/vcXXghffQW//CX89repZLKZtSz1Jf62wHKUnPlb5Zk3Dy6/HM45B6ZPT+Pb/u53qaummbVM9SX+qRFxdpNFYs3KwoVw660wdChMmpTKLPzhD9CntmF1zKxFqa87p8/0K9Qjj6RxbQ84AJZfHh58MM1z0jdrHepL/B4escKMHQu77go775zKLfz97/DSS2mei6iZtR51Jv6IcAXOCjFpEhx4IPTuDWPGpAu4b76Z5rUp5xY/M2tRyinZYK3U9OlptKu//Q3atoVTToGTToIVVig6MjPLkxN/BfryS/jzn9PF2s8/h0MPTYOar7FG0ZGZWVNw4q8g8+fDtdfCGWekG7H22AN+/3vYcMOiIzOzpuQW3AoQAXffDZtsAkOGwPe+B089leY56ZtVnlwTv6T/SnpV0suSxmTzekl6rmqeJI+0mqNnnoGtt4Y990zTw4fD00+neWZWmZrijL9/RPSKiKpe4H8EzoqIXsDp2bQ1svHj4Sc/SQl+0iQYNgxefTXNc9dMs8pWRFNPAMtnz1cAphQQQ6s1dSocfniqgz9yZCqv8NZbaV47X9ExM/K/uBvAw5ICuCIihgG/AR6SdAHpi2fL2laUNAQYAtC1a9ecw2wdnngC9t4bPv0UjjoqFVHr3LnoqMysuck78W8VEVMkrQKMkPQGsDdwTETcKWkf4GpSrf9vyL4khgH06dMnco6zRYtIhdSOPhrWWQeefBK+//2iozKz5irXpp6ImJL9+xEwHOgLDALuyl5yezbPFtO8efDzn6cyybvsAqNHO+mbWf1yS/ySlpXUqeo5sAswjtSmv132sh2At/KKobX78MNUNfPKK9Ndt//8J6y4YtFRmVlzl2dTz3eB4UpdSNoB/4iIByV9Dlycjek7h6wd3xpmzJjURfPjj+Hmm1OdfDOzcuSW+CNiIrBpLfOfBjbLa7+V4Kab4LDD0uhXzzwDP/hB0RGZWUviO3dbkAUL4IQTUtXMzTeHF15w0jezhnPP7hZi5kzYf3946CH4xS9SkbWlly46KjNriZz4W4Dx4+HHP4Z33oErrkj1dszMFpcTfzP3r3/Bz34GHTqkO3FdY8fMlpTb+JupiDRIyh57QPfuqRePk76ZNQaf8TdDX3wBhxwCt9+eBjy/8kro2LHoqMystXDib2YmTUoVNF99Ff74Rzj+eFfTNLPG5cTfjDz2GAwcmEbKuv9+GDCg6IjMrDVyG38zEKhdUtwAAAz3SURBVAGXXAI77wxdusDzzzvpm1l+nPgLNnduqpV/1FHwwx+mImvrr190VGbWmjnxF2jqVOjfH66+GoYOhXvugeWXX/R6ZmZLwm38BXnhhVRkbeZMuO221LZvZtYUfMZfgL//HbbZJg2F+OyzTvpm1rSc+JvQ/Plw3HFw0EGwxRbprH/Tb9UvNTPLl5t6msiMGalm/ogRcOSRcOGFsNRSRUdlZpXIib8JvPZaKr3w7rvpLtzDDis6IjOrZE78ObvnnlQ/f9ll4fHHYcsti47IzCqd2/hzsnAhnH12Kr+wwQapyJqTvpk1Bz7jz8Hnn8OgQXDXXelsf9iwVFbZzKw5cOJvZBMnprP8116DP/0JjjnGRdbMrHlx4m9Ejz4K++yTau88+GCqvWNm1ty4jb8RRMBf/gK77gqrrpqKrDnpm1lz5cS/hObOhcGD4de/ht13h+eeg/XWKzoqM7O6OfEvgSlTYLvt4Npr4fTT08XcTp2KjsrMrH65tvFL+i8wC1gAzI+IPtn8o4AjgfnAfRFxYp5x5GH06FRk7bPP4I474Kc/LToiM7PyNMXF3f4RMb1qQlJ/YA+gZ0TMlbRKE8TQqK6/HoYMgTXWgIcegk02KToiM7PyFdHU8wvgvIiYCxARHxUQw2KZPz91zzz4YNh661RkzUnfzFqavBN/AA9LelHSkGze+sA2kkZLekLS5rWtKGmIpDGSxkybNi3nMBft44/TcIh//nO6kPvQQ7DyykVHZWbWcHk39WwVEVOy5pwRkt7I9rkS0A/YHLhN0joREaUrRsQwYBhAnz59ggKNG5eKrE2eDNdcA4ccUmQ0ZmZLJtcz/oiYkv37ETAc6AtMBu6K5HlgIdA5zziWxF13Qb9+MHs2PPGEk76ZtXy5JX5Jy0rqVPUc2AUYB9wN7JDNXx9YGphe13aKsnAhnHlm6q2z0UapyFq/fkVHZWa25PJs6vkuMFypUE074B8R8aCkpYFrJI0D5gGDajbzFG3WrDRK1t13p2Jrl18O7dsXHZWZWePILfFHxETgWwMLRsQ84MC89rukJkxI7flvvJEu5B59tIusmVnr4iJtJR55JBVZg1Rkbaedio3HzCwPLtlAKrJ20UWpyNoaa6T++U76ZtZaVXzinzMn3ZB17LGpiWfUKFh33aKjMjPLT0Un/vffT0XWbrgBzjor1dxZbrmiozIzy1fFtvGPGgV77ZWGSRw+PI2aZWZWCSryjP/aa2H77aFjx/QF4KRvZpWkohL/V1+l7pmHHgrbbpsu4m68cdFRmZk1rYpJ/NOnp147f/1rqrD5wAPwne8UHZWZWdOriDb+V15JPXamTk219A86qOiIzMyK0+rP+O+4A7bYAubNgyefdNI3M2vVif/3v4eBA6Fnz1RkrW/foiMyMyteq0783bvD4MHw+OOw2mpFR2Nm1jy06jb+gQPTw8zMqrXqM34zM/s2J34zswrjxG9mVmGc+M3MKowTv5lZhXHiNzOrME78ZmYVxonfzKzCKCKKjmGRJE0D3lnM1TsD0xsxnJbAx1wZfMyVYUmOee2I6FJzZotI/EtC0piI6FN0HE3Jx1wZfMyVIY9jdlOPmVmFceI3M6swlZD4hxUdQAF8zJXBx1wZGv2YW30bv5mZfVMlnPGbmVkJJ34zswpTEYlf0jmSXpH0sqSHJa1edEx5k3S+pDey4x4uacWiY8qbpIGSXpO0UFKr7fInaYCkNyW9LenkouNpCpKukfSRpHFFx9IUJK0l6TFJ47O/6V835vYrIvED50dEz4joBdwLnF50QE1gBLBxRPQE/gOcUnA8TWEcsBfwZNGB5EVSW+BS4IfAhsD+kjYsNqomcR0woOggmtB84LiI+D7QD/hVY37OFZH4I+KzksllgVZ/RTsiHo6I+dnkc8CaRcbTFCJifES8WXQcOesLvB0REyNiHnALsEfBMeUuIp4EZhQdR1OJiKkR8VL2fBYwHlijsbbfqsfcLSXpXOAg4FOgf8HhNLVDgVuLDsIaxRrAeyXTk4H/KSgWawKSugE/AEY31jZbTeKX9Aiwai2LhkbEPRExFBgq6RTgSOCMJg0wB4s65uw1Q0k/G29qytjyUs4xt3KqZV6r/wVbqSQtB9wJ/KZGy8USaTWJPyJ2KvOl/wDuoxUk/kUds6RBwO7AjtFKbthowOfcWk0G1iqZXhOYUlAsliNJS5GS/k0RcVdjbrsi2vgldS+Z/DHwRlGxNBVJA4CTgB9HxJdFx2ON5gWgu6TvSVoa2A/4Z8ExWSOTJOBqYHxEXNjo228lJ4L1knQn0ANYSCrvfEREvF9sVPmS9DawDPBxNuu5iDiiwJByJ2lP4K9AF+AT4OWI2LXYqBqfpB8BfwbaAtdExLkFh5Q7STcD25NKFH8InBERVxcaVI4kbQ08BbxKylsAp0bE/Y2y/UpI/GZmVq0imnrMzKyaE7+ZWYVx4jczqzBO/GZmFcaJ38yswjjxV5C6Khw2pKqlpGMkzZG0Qj2vCUl/L5luJ2mapHsb4Ri+J2m0pLck3Zr1Za/5mu1r7kvSdZL2XtL9NyDOQyW9mlVHHSdpj2z+2ZJyuwktO85JWSXaf0vaMa991bH/MyUdnz3P9Vht8TnxV5brqL3CYUOqWu5Puoloz3pe8wWwsaQO2fTOQGPdN/EH4KKI6A7MBAY30nbLImmRd7tLWhMYCmydVUftB7wCEBGnR8Qj+UbJCVkl2t8Al+e1k6xSaJ2a6FhtMTjxV5C6KhyWW9VS0rrAcsBvSV8A9XkA2C17vj9wc8OirXX/AnYA7shmXQ/8ZDG2s6OksdkZ+TWSlsnm/1dS5+x5H0mPZ8/PlDRM0sPADZI2kvR8dlb9So07wwFWAWYBnwNExOcRMSnb1te/PLL9nSXppSyWDbL5y0m6tuQXw0+z+btIGpW9/vasjkt9RlFS0VHSZpKekPSipIckrZbNX0/SI9kvhJckravk/OzXyquS9s1eu71Snfh/kG4uQtJQpfEBHiHdKFm1v3KOtYukEdn8KyS9U/UZWH6c+K0hqhL4U0APSavU89pbgP0ktQd6UkdlQUk9sgRa26Pm4DErA5+UlJueTN2larcp3RapVAdZPNcB+0bEJqR6Vb9Y1IEDmwF7RMQBwBHAxdlZdZ8sjlL/Jt1dOilL4P+vnu1Oj4jewGXA8dm804BPI2KT7BfDyCwZ/hbYKXv9GODYRcQ8ALg7O+6lSHc17x0RmwHXAFV3/N4EXBoRmwJbAlNJvwB7AZsCOwHnV31RkEpDD42IDSVtRiob8YNsnc0beKxnACOz+cOBros4JmsEraZImzWJ/YA9I2KhpLuAgaRBQb4lIl5RKie7P1DnbebZL41eZe6/IZUpn4qI3b9eUboue9oDmBQR/8mmrwd+RSqBUJ9/RsTs7PkoUqXXNYG7IuKtbwQUsUCpVtLmwI7ARZI2i4gza9luVfGtF0mJE1Ki3a9kezMl7U4aeOWZ9MOHpbM4anO+pD+Sfnn0KznujYER2fptgamSOgFrRMTwbF9z4OuSATdHxALgQ0lPZMfzGfB81S8YYBtgeFU9KEn11Q2q7Vi3Jms2jIgHJc2sZ31rJE78VhZJPYHuVCeOpYGJ1JH4M/8ELiDVWFm5ju32oO6xAraPiE9KpqcDK0pql531L05lytq+PKrMp/pXcPsay76oehIR/5A0mtSU9ZCkwyJiZOmLs2qozwPPSxoBXAucWcs+52b/LqD6/6P49heagBERsagmNoATSEn2aNIX22bZ+q9FxBbf2Ki0fB3bqO99+qLGdLl1X+o6Vmtibuqxcu0PnBkR3bLH6sAaktauZ51rgLMj4tW6XhARb0ZErzoen9R4bQCPAVW9cwYBDa3B/wbQTdJ62fT/Ak9kz/9LSpIAP61rA5LWASZGxF9IX249ayxfXVLvklm9SMUBy/UwacyIqu2tRBpFbauquCV1lLR+XRuIiIXAxUAbSbsCbwJdJG2Rrb+UpI2yGu+TJf0km7+MpI6kC/37SmorqQuwLemLrKYngT0ldch+PdTXrFWbp4F9sn3vAqzUwPVtMTjxVxClCoejSO3zkyUNzubvKWkysAVwn6SHall9P1IbbKnhpHb81SV9qzknIiZHxMWNexScBByrVH10ZVLp2rJlTRmHALdLqqp8WNXz5SzgYklPkc5K67IvMC67drABcEON5UsBFygNdv9y9vqGDJb9O2Cl7MLqv4H+ETENOBi4WdIrpC+CDerbSPZF+TvgxGyYxr2BP2TbfJnUng/py+/obLvPkga6GU7qifRvYGS2jQ9q2cdLpF9sL5Nqxz/VgOOE9J7vIukl0jjCU0kXxi1Hrs5pZoVR6lG1ICLmZ79GLssumluO3MZvZkXqCtwmqQ0wDzi84Hgqgs/4zcwqjNv4zcwqjBO/mVmFceI3M6swTvxmZhXGid/MrML8f1i3448O+KlnAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "# 4a.\n", | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "time=np.array([-3,-2,-1,0,1,2])\n", | |
| "Temp=np.array([55,58,60,65,66,67])\n", | |
| "plt.plot(time,Temp,'-',color=\"blue\")\n", | |
| "plt.title('Temperature as a Function of Time')\n", | |
| "plt.xlabel('11 A.M = 0 Hours Since Recording')\n", | |
| "plt.ylabel('Temperature(F)')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "This plot does not look correct since temperature would not be linear. For a better way to get $T_a(t)$ would require more values recorded for a continous function. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 172, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stderr", | |
| "output_type": "stream", | |
| "text": [ | |
| "No handles with labels found to put in legend.\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x7fd02dc1f780>" | |
| ] | |
| }, | |
| "execution_count": 172, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3debzV49rH8c+3QXNK7UiDQkpSqa2kEIkMDzqkMhyHDOfg4DjGJxzzMT2O+TgZcyRzZioZMpTsUO1CotJWKk2aVft6/rh/e1u2PdZe+7f2Xtf79VqvtX7ztfawrnUPv/uWmeGcc84BVIs7AOecc6nDk4Jzzrl8nhScc87l86TgnHMunycF55xz+TwpOOecy+dJwbltIKmvpG8lrZE0IMnX6iBpczKvERdJj0u6LO44nCeFKin6gMp75Epan7B8ctzxbQtJP0rqE3ccCW4CbjOz+mb2VsGNUbzrop/9ckmvSNo5hjiRNFBSlqS1kn6S9ISk5hVw3T0K/E1aFEPe8r5m9iczuy3ZsbiSeVKogqIPqPpmVh/4HvifhHWj4o6vKJJqVMJr7ALMLGGfw6LfRQtgNXBnOcdQoujLwGPAbUAToDNQHZgoqWE5X+s3P2Mzm53w99g0Wt0+4W/y0/K8vts2nhTSkKTqkq6W9F30jXGUpEbRtg6SNksaJukHScsknSGpl6RsSSsl3Zlwrj9LekfSfyT9LGmWpAMTtu8QfSP9UdICSf+QVK3AsfdLWgFcEV3/vehb9VJJIyU1iPZ/DmgGjIu+YV4gaYCkOQXeX35pQtItkp6S9Iyk1cCQ4t5/ET+v86IqomWSXpS0Y7Q+B9g5L56Sfu5mth54AehY4OfzVPRe50q6TJKibTUk3R1ddw7QP+G4UyV9VCDO4ZKeLiT+GsAdwDVm9qyZbTCzhcAfo13Ol1Qv+pnunnBci6iU2ThaHihpevQ38IGkxPfxo6RLJM0Efi7pZ1FIjE9Luip6PUDSHElXRb+fHyQdKenYhN/D3xOOLdPv0xXPk0J6uhQ4DOgDtAQ2Af9K2F6d8E1yV+B04F7gEuCgaP3pknom7H8gMI3wDfQW4KWEb5+jgFXRuXoAxwGnFjj2C8I3yP+L1l0P7ATsDbQHhgOY2SBgCdE3bzO7p5Tv93hgJLA94UO5pPefT9KRwNXAQMI3/Z+AJ6N4WibGU1IQkuoDg4DJCasfBGoCbQkf+n8BToq2nQ8cQvg59AIGJxz3IrC3pF0T1p0M/LeQS3ci/DyfS1xpZlui8/Q3s7XAK8DQhF2GAGPNbIWk/YAHCH8PTaLrvFSgVDA4eg9Niv4plNouhN/LToS/qUeBEwh/f4cCN0lqEe1b6t+nKwUz80cVfgDzgEMLrJsL9E5YbgusAwR0AAxokrB9LXBswvLrwJ+j138G5hY4/3TCh98u0bE1E7adDryZcOzsEuIfAkxKWP4R6JOwPACYU+CY/H0IHyjjSvv+C7n+KOD6hOVGQC6wU2HxFHL8j4Qqo5XAZmABsGe0rRawBdg1Yf8Lgbei1x8Df0rYdgywOWH5MeDq6HUmIUHVKCSGQ6OYqxWy7SJgRvT6aGBWwrapwIkJ1xpe4Nj5QM+E93lSKf4ea0d/Xy0LrH8auCrhd7oqL14gIzqmS8L+M4EBZf19+qPkR9LrcF1qiaomWgFvSEocDbEav37D22JmyxK2rQcWF1hO/GacU+Ay8wnVKrsQPgSWRjUieddJrO5ZUCC+nYG7gf2BBtH+i0rz3oqRf41SvP+fChy7M/BO3oKZrZT0M6HU8GMpr3+EmX0YfaseRKjH70D4GVYjtPvkmR+dO+/aCwpsSzSSUNK4ATgFGG1mhfVO+omQ8Hfk9z/L5vz6nscCIyV1IfyO2wGvRtt2AU6UdGnCsdslxEqBWLfVUjPLjV6vj55/9ze4Fb9PVwKvPkozFr5K/QAcYmaNEh61zWxr/4FaFlhuDSwkfEisARonXKehmXVLDKnAsbcTShedzKwhcCbhA62o/dcCdfMWJNUEdiiwT/4xW/H+FxI+EPPOvz3QMDpHmZjZZjMbTUiUvQhJJZfw88rTOuHciwgfeInbEr0P1I6qdoZQeNURQDbhA3VQ4kpJ1QnVYhOi+DYBzxOqkE4GxlhoB4Hwu7ymwM+srpm9mPgWi3v/yZCkv+e05kkhPT0I3CKpFYCkZpL+ZxvO1ypqNK4h6RTCh9c4M5tLqD+/TVIDSdUktVPxXUobEBLJz5JaAxcX2L6Y0D6R50tgB0n9ooRwHSX/XZfl/Y8GzpLUSVJt4FbgHTMrbSkhX/T+BwF1gK/MbCMwBrg5aujdjVB99GR0yLPA3yQ1l9QU+E0//ugD8b/ACGC5mWUVdt2o9HA5cIOkQZJqRSWykUAN4L6E3Z8iJJih0es8I4C/SspUUF/SMZLqEr/y/ntOa54U0tNtwNvAOwo9cj4GuhV/SLEmAvsAywmNwgPNbFW0bSihHv6raPszhGqMolxDaDBcRfjAfKHA9psIjYwrJZ0ffRu8kFD3n0P49l3SN8RSv38zew34J6ERdiGh4fPUwvYtRl7vpFWERuuTzCyvCu2c6Hk+oZrq4ei9QPiw/oBQf/4JIUkUNJLQEF1UKSHvfYwEhgFXACuAGdGmPma2MmHXiYSOBtsTfkZ5x38EXAD8h9A+MpvQIJ4KE7KU999zWlP4suHc1pH0Z+AEMzs07ljSkUJ33cVABzP7vqT9nSuJlxScq9z+CrznCcGVF+995FwlJelHQtfLY+KOxVUdXn3knHMun1cfOeecy1epq4+aNm1qbdq0iTsM55yrVKZOnfqTmWUUtq1SJ4U2bdqQlVVo12znnHNFkFTw7vh8Xn3knHMunycF55xz+TwpOOecy1ep2xSccy7dbdq0iZycHDZs2PC7bbVr16Zly5bUrFmz1OfzpOCcc5VYTk4ODRo0oE2bNiQMUY+ZsWzZMnJycmjbtm2pz+fVR845V4lt2LCBJk2a/CYhAEiiSZMmhZYgiuNJwTnnKrmCCaGk9cVJz6SwfDlcdBGsXFnyvs45l0bSs01h7ly4915Yvx7+85+4o3HOuZSRniWF7t3h4othxAh4//24o3HOuW1S1MCmWzPgaXomBYDrroPddoMzzwwlBuecq4Rq167NsmXLfpcA8nof1a5du0znS8/qI4C6dUNJoV8/uPZauPXWuCNyzrkya9myJTk5OSxduvR32/LuUyiL9E0KAIccAsOGwf/9HwweDN18WlfnXOVSs2bNMt2HUJL0rT7Kc/vtkJERksOmTXFH45xzsfKk0Lgx3H8/fPFFKDE451wa86QA8Ic/hMe118Ls2XFH45xzsfGkkOe++6BOHTjrLMjNjTsa55yLhSeFPM2bh+qjiRPhoYfijsY552KRtKQg6VFJSyRlJ6y7QdJ0SV9IGidp54RtV0qaI+lrSYcnK65inX566KJ62WXwww+xhOCcc3FKZknhcWBAgXW3m1lnM+sKvAZcAyCpIzAE2Cs65gFJ1ZMYW+GkcO/Cpk1w7rmwFXcDOudcZZa0pGBmE4HlBdb9nLBYD8j71D0WeNrMNprZXGAO0CNZsRVr113hhhvglVfguediCcE55+JS4W0Kkm6StAA4maikALQAFiTslhOti8eFF0JmJvz1r7BsWWxhOOdcRavwpGBmw82sFTAKOD9aXdig34XW3Ug6W1KWpKzCbusuFzVqwCOPhCG2//735FzDOedSUJy9j54Cjo9e5wCtEra1BBYWdpCZjTCzTDPLzMjISF50nTvD5ZfDyJEwblzyruOccymkQpOCpHYJi8cAX0WvXwGGSKolqS3QDphSkbEV6qqroH17OOccWLMm7miccy7pktkldTQwCWgvKUfSMOAWSdmSpgOHARcCmNlM4FlgFvAWcJ6ZbUlWbKVWuzY8/DDMmwdXXx13NM45l3TamkkYUkVmZqZlZWUl/0LnnQf//jdMmgQ9eyb/es45l0SSpppZZmHb/I7m0vjnP6FFizAhzy+/xB2Nc84ljSeF0mjYEB58ELKz4ZZb4o7GOeeSxpNCaR11FAwdCjfeCLNmxR2Nc84lhSeFsrj77lBqGDYMtsTfDu6cc+XNk0JZZGTAXXfB5MlhYh7nnKtiPCmU1cknwxFHwP/+L8yfH3c0zjlXrjwplJUUuqdCuKmtEnfpdc65gjwpbI1ddgndVMeO9Ql5nHNViieFrXXeedC/P1x0kfdGcs5VGZ4Utla1amGwvHr1QlfVDRvijsg557aZJ4Vt0bx5SAzTp4cRVZ1zrpLzpLCtjjwyTMpzzz3w2mtxR+Occ9vEk0J5uPVW6NIFTj8dFi2KOxrnnNtqnhTKQ61aMHo0rF0Lf/wj5ObGHZFzzm0VTwrlZc89wzAYb78Nd9wRdzTOObdVPCmUpzPPhBNOgOHDYUr8E8c551xZeVIoTxKMGAE77wwnnQSrV8cdkXPOlYknhfLWuDGMGgVz54Yb3JxzrhLxpJAMffqEOZ3/+9+QIJxzrpLwpJAsV10VksNf/gLffht3NM45VyqeFJKlRo1QSqhePbQvbNoUd0TOOVciTwrJ1Lp1GEV1yhS45pq4o3HOuRJ5Uki2E04IXVVvvRXeeSfuaJxzrlhJSwqSHpW0RFJ2wrrbJX0labqkMZIaJWy7UtIcSV9LOjxZccXirrugfXs45RT46ae4o3HOuSIls6TwODCgwLrxQCcz6wzMBq4EkNQRGALsFR3zgKTqSYytYtWrF4bBWLYMzjjDZ2tzzqWspCUFM5sILC+wbpyZbY4WJwMto9fHAk+b2UYzmwvMAXokK7ZYdO0Kt90Gr74K998fdzTOOVeoONsUzgDejF63ABYkbMuJ1v2OpLMlZUnKWrp0aZJDLGcXXBCG2r7kkjAHg3POpZhYkoKk4cBmIO/OLhWyW6F1LGY2wswyzSwzIyMjWSEmhwSPPQaNGoXZ2tatizsi55z7jQpPCpJOA44GTjbLr1zPAVol7NYSWFjRsVWIZs3giSfCvM5//3vc0Tjn3G9UaFKQNAC4HDjGzBK/Jr8CDJFUS1JboB1QdYcZPewwuPRSePBBePHFuKNxzrl8yeySOhqYBLSXlCNpGHAf0AAYL+kLSQ8CmNlM4FlgFvAWcJ6ZbUlWbCnhxhshMzPcwzB3btzROOccALJK3D0yMzPTsrKy4g5j682ZA/vuG+58/ugjqF8/7oicc2lA0lQzyyxsm9/RHKfdd4enn4bsbPjTn/z+Bedc7DwpxO3ww8P9Cy+8ADfdFHc0zrk050khFVx8cRgC4+qr4ZVX4o7GOZfGPCmkgrxpPLt3D8lh1qy4I3LOpSlPCqmiTh146SWoWxeOPRZWrIg7IudcGvKkkEpatgxtC/Pnw5AhsHlzycc451w58qSQanr3DgPmjRsHV14ZdzTOuTRTI+4AXCHOOgumTYM77oAuXUI7g3POVQAvKaSqf/0LDjoo3PFcmW/Qc85VKp4UUlXNmvDcc7DTTnDccfDjj3FH5JxLA54UUllGRuiRtHw5HH88bNwYd0TOuSrOk0Kq69oVHn8cPv4Yzj/fh8JwziWVNzRXBieeGBqeb74Z9tkHzj037oicc1WUlxQqixtugKOPhgsvhPfeizsa51wV5UmhsqhWDZ58MoysOmgQzJsXd0TOuSrIk0Jlsv328PLLsGlT6JG0dm3cETnnqhhPCpXNHnuEORhmzIDTT/eGZ+dcufKkUBkNGAD//Ge4j+GWW+KOxjlXhXhSqKwuvRSGDoXhw+G11+KOxjlXRXhSqKwkePjh0EX1pJPg88/jjsg5VwV4UqjM6tYNDc+NGsERR8C338YdkXOukvOkUNm1bBmG2d68GQ47zMdIcs5tE08KVUGHDvD66yEhHHEErFoVd0TOuUoqaUlB0qOSlkjKTlg3SNJMSbmSMgvsf6WkOZK+lnR4suKqsnr2hBdfhOzscA/Dhg1xR+Scq4SSWVJ4HBhQYF028AdgYuJKSR2BIcBe0TEPSKqexNiqpsMPD4PnvfcenHwybNkSd0TOuUomaUnBzCYCywus+9LMvi5k92OBp81so5nNBeYAPZIVW5V28slw112h1HDuuX5zm3OuTFJllNQWwOSE5Zxo3e9IOhs4G6B169bJj6wyuvBCWLw43OC2445w/fVxR+ScqyRSJSmokHWFfsU1sxHACIDMzEz/GlyUm26CJUvC6KrNmoW5GJxzrgSpkhRygFYJyy2BhTHFUjVI8OCD8NNPcMEFYRa3wYPjjso5l+JSpUvqK8AQSbUktQXaAVNijqnyq1EDRo+GPn3g1FNh/Pi4I3LOpbhiSwqStgOOBA4AdgbWE3oQvWFmX5Vw7GigL9BUUg7wD0LD871ABvC6pC/M7HAzmynpWWAWsBk4z8y860x5qFMHXnkFDjoIBg4MPZMyM0s8zDmXnmRF9E6RdBVwPKH76FRgCVAb2AM4mNAOcImZZRd6ggqQmZlpWVlZcV2+clm0CPbfH9asgY8+CkNwO+fSkqSpZlbot8PiSgozzOzGIrbdJqk5v20HcKmsefMwHEbv3mE4jI8/hp13jjsq51yKKa5N4ZXiDjSzRWbm9f6VSbt28OabsGxZuNFtxYq4I3LOpZjiksLUvBeS7qqAWFxF6N4dXnoJZs+GY46B9evjjsg5l0KKSwqJ9w4cmOxAXAXq1w+efDK0LQweHEZYdc45ik8KfmNYVTZoENx3H7z6Kpx9tg+H4ZwDim9o7iDpM0KJoX30mmjZzKxb0qNzyXXuueGu5+uuCze33XJLuOnNOZe2iksKe1dYFC4+//hHSAy33QbVqsHNN3ticC6NFZkUzMzndkwHUqhGMgslhV9+gTvu8MTgXJoqMilIehd4FnjZzBYmrK8B7A+cBnxoZo8lPUqXXNWqwQMPQM2acOedsGkT3H23Jwbn0lBx1UdHAWcCYyS1IAxRUQeoBUwA7jczv524qpBCIkhMDPffHxKGcy5tFFd9tA64B7hHUi2gGbDezH6qqOBcBZNC1dF22/1alTRiBFT3SfCcSxelHTq7PdAHMEkfmtmMJMbk4iSFxubttguT82zaBI895onBuTRRYlKQNBw4CXgpWjVa0igz+2dSI3PxkUI31Zo14eqrQ2L473/DUNzOuSqtNP/lpwDdo+okJN1EGALDk0JVd9VVITFccUW46/mpp8Kyc67KKk1SmF9gvxrAd8kJx6Wcyy8PVUkXXxxKDM88A7VqxR2Vcy5JSpMU1gEzJY0lDH1xGPChpDsBzOziJMbnUsHf/hYSw/nnw/HHw/PPQ+3acUflnEuC0iSF16NHnslJisWlsvPOC1VH55wDxx4bRlqtUyfuqJxz5azEpGBmj1REIK4SOPvs0Nh85plw9NFhms969eKOyjlXjkq8M0nSAEmfSloiabmkFZKWV0RwLgWdcQY88USY6/nII2H16rgjcs6Vo9JUH90HnAjMAHKTG46rFE45JZQYTjkFBgwIs7k1bBh3VM65clCaMQxygC/MbJOZbcl7JDswl+KGDIGnn4YpU6B/f1i5Mu6InHPloDQlhcuAVyW9B2zMW2lm9yQrKFdJnHBCaHweNCjM5jZ+POywQ9xROee2QWlKCtcBW4BGQEbCo1iSHo3aIbIT1u0gabykb6LnxgnbrpQ0R9LXkg4v+1txscjriTRzJvTtCzk5cUfknNsGpUkKzczsGDMbbmZX5z1KcdzjwIAC664AJphZO8JIq1cASOoIDAH2io55QJIPtlNZHHkkvPYazJsH++0H06bFHZFzbiuVJilMkHRIWU9sZhMJw20nOhYYGb0eCRyXsP5pM9toZnOBOUCPsl7TxejQQ+HDD8PrAw6AsWPjjcc5t1VKkxTOAt6WtKYcuqTuaGaLAKLnZtH6FsCChP1yonW/I+lsSVmSspYuXbqVYbik6NwZPvkEdt0VjjoKHvFbXJyrbEqTFJoCNYHtCW0JTSlFm0IZFTbFlxW2o5mNMLNMM8vMyCjvMNw2a9ECJk4MJYczzwyD6lmhv0rnXAoqMSlE3U8HAZdHr5sDXbfyeoslNQeInpdE63OAVgn7tQQW4iqnhg3h1Vdh2DC46SY49VTYuLHk45xzsSvNHc33AQcDp0ar1gEPbuX1XiHM7Uz0/HLC+iGSaklqC7QDpmzlNVwqqFkTHnoIbrwRRo0KN7mtWBF3VM65EpSm+mh/MzsH2ABgZsuB7Uo6SNJoYBLQXlKOpGHALUB/Sd8A/aNlzGwm8CwwC3gLOM9vkKsCJBg+HJ58Ej76CHr3Dj2UnHMpqzQ3r22SVI2ojl9SE0ox3IWZDS1iU78i9r8JuKkU8bjK5uSTQ1vDwIGhy+prr0FmZtxROecKUWRJQVJewrgfeAHIkHQd8CFwawXE5qqSvn1DaaF2bTjooJAYnHMpp7jqoykAZvYEcBVwB7ACGGRmT1dAbK6q6dgRJk+GPfcMd0I/8EDcETnnCiiu+ii/m2hU5z8z+eG4Km+nneD992Ho0DBxz9y5cOutUK00zVvOuWQrLilkSCpyqk0zuzMJ8bh0UK8ejBkDF1wAd9wB8+eHORp8ik/nYldcUqgO1KfwG8uc2zbVq8N990HbtnDppfDDD/Dyy9C0adyROZfWiksKi8zs+gqLxKUfCS65BHbZJdzgtv/+8MYbsPvucUfmXNoqriLXSwiuYgwaBBMmwPLl0KNHSAzOuVgUlxQKvZ/AuaTo3TvM4ta6dRhM75prYIvfv+hcRSsyKUR3LjtXcXbdFSZNgj/9CW64IczTsGxZ3FE5l1a8H6BLLXXqwKOPwogR8N570K0bfPpp3FE5lzY8KbjUI8FZZ4U7oAH69AlJwofgdi7pPCm41JWZCZ99BgcfDOecA6efDuvWxR2Vc1WaJwWX2po0gddfDw3PI0eGbqvffht3VM5VWZ4UXOqrXh2uuy4kh++/h+7dwyQ+zrly50nBVR5HHglTp8Juu8Exx4S5GrzbqnPlypOCq1zatg0N0MOGwc03hxndli6NOyrnqgxPCq7yqV0bHn44PD74IFQnffJJ3FE5VyV4UnCV17Bh8PHHoc3hgAPC/AzebdW5beJJwVVu3bqFdob+/cP8DH/8o3dbdW4beFJwld8OO4TeSNdfD6NGheqkqVPjjsq5SsmTgqsaqlWDq6+GceNg9WrYb78wftLmzXFH5lyl4knBVS2HHgozZoThuK+5JgyRMXt23FE5V2l4UnBVT+PG8NRT8PTTISHssw/8+9/eCO1cKcSSFCRdKClb0kxJF0XrdpA0XtI30XPjOGJzVcjgwaHU0KcPnHtuuPlt4cK4o3IupVV4UpDUCTgL6AF0AY6W1A64AphgZu2ACdGyc9umRQt46y24/354/33Ye2949tm4o3IuZcVRUtgTmGxm68xsM/A+MBA4FhgZ7TMSOC6G2FxVJIWSwuefh/mfBw+GU06BlSvjjsy5lBNHUsgGDpTURFJd4EigFbCjmS0CiJ6bFXawpLMlZUnKWurDG7iyaN8+DJFx3XWhvWHvvcPc0M65fBWeFMzsS+BWYDzwFjANKHW/QTMbYWaZZpaZkZGRpChdlVWjRuiVNHky1K8feitdeCGsXx93ZM6lhFgams3sETPrZmYHAsuBb4DFkpoDRM9L4ojNpYm8CXwuuADuuSfcGZ2VFXdUzsUurt5HzaLn1sAfgNHAK8Bp0S6nAS/HEZtLI3XqwN13w/jxsGYN9OrlN7y5tBfXfQovSJoFvAqcZ2YrgFuA/pK+AfpHy84l36GHwvTpcOKJoWqpd2/46qu4o3IuFnFVHx1gZh3NrIuZTYjWLTOzfmbWLnpeHkdsLk01bhzGTXrmGZgzB7p0CQnC2xpcmvE7mp1LdOKJMGtW6LZ6ww2hh9LYsXFH5VyF8aTgXEE77ghPPBG6q9aoEWZ3GzzY74Z2acGTgnNFOeQQmDYtlBhefhk6dIB77/V5oV2V5knBueLUqgVXXQXZ2aF30gUXQM+e3n3VVVmeFJwrjd13D2MoPfNMqEbq0QP++ldYtSruyJwrV54UnCstKTREf/klnH9+mBO6Q4cwZIYPy+2qCE8KzpXV9tuHu6A/+SSMwjp0KBx+OHzzTdyRObfNPCk4t7UyM0NiuO++8Lz33mGe6I0b447Mua3mScG5bVG9Opx3XrgDeuBA+Mc/oHNnlj3zetyRObdVPCk4Vx6aN4fRo2HsWNZs2o7dh+zP4S0/4I0nppC7JTfu6JwrNU8KzpWnww7DPniHS/Z9kew1u3PUaT3Ys/U87rv6XVYv855KLvV5UnCunDVokcHwKcOYt7gpT909icYN1/PXGw+mZSu4+OR3+G7ad3GH6FyRPCk4lyQ1a9Vk6AW9mPzlXkwe+yVH9ZnFvc8cwO77tOG4Ph/z7nOfYLleteRSiycF5ypAz8P25KlxvZj31UquHPYxH07vwCEn9qTrbrN55KYJrP95ddwhOgd4UnCuQrXYPYObHurDgoX1efiWTzBqcOZV/WjdciPDzxjHvOx5cYfo0pwnBediUKf+dgy7vCfTvtudd8Z8Se+u8/jn44eya+fWHL7vZzz/0DR+2eh3SbuK50nBuRhJcPBxe/LSxEzmfbWUa86ZyKy5OzLo7C602mkZlw37lNkzlsUdpksjnhScSxGt99iRa//dl3kLm/L6I++wf8ds7nx8H9p3bkLffbIZde9nbFjn80e75PKk4FyKqb5dLY484xDGfNSXBbO+4+bz32DBovqcckE3dt5xNRee9AEzPv427jBdFSWrxKM7ZmZmWpaPa+/SQO7mTbz7/Kc89NAWxkzswS+ba9GzQzbDTlnOoLO60KjZ9nGH6CoRSVPNLLOwbV5ScK4SqFajJv2G7M/TEw7gh7k/c+fl7/Hz2jqcfdWB7NiiNsf3ncKYR6ezcYPf9+C2jZcUnKukLNfImvAlox79kdFvdmLJqmY0qreKQQPmcMqZO9PnsOZU8699rhApV1KQ9DdJMyVlSxotqbakHSSNl/RN9Nw4jticqyxUTezbvyN3jT6EH35swJuPv8fRvT5l1GsdOOiI5rRtvpgr/5LNzC/WxR2qq0QqvKQgqQXwIdDRzNZLehZ4A+gILDezWyRdATQ2s8uLO5eXFJz7vTWLF/DyY1mMer4x4z7vw5bcGnTZfS6nnCXNcOEAABDNSURBVPgzQ89uR4td6sYdootZypUUgBpAHUk1gLrAQuBYYGS0fSRwXEyxOVep1d+xFSdfMZA3Pj2QH6Zlcfffn6eWVnDpzV1o1bY2B3WdxT3XfUHOvPVxh+pSUCxtCpIuBG4C1gPjzOxkSSvNrFHCPivMrNgqJC8pOFdKuVv45pMsRj36E8+P3ZWZC/YEYL+9ZnP8cb9w/J/a0Xb3WjEH6SpKcSWFOKqPGgMvAIOBlcBzwPPAfaVJCpLOBs4GaN26dff58+dXSNzOVRm5m/n64yxeeOpHnn9rVz6f2xmAbu3ncvyxGzj+tN1o33G7mIN0yZRqSWEQMMDMhkXLfwT2A/oBfc1skaTmwHtm1r64c3lJwbltlLuZ76ZM5sWnFvHCW22Y/M2+AHTadQHHH7ue4//Ylk5daiLFHKcrV6mWFHoCjwL7EqqPHgeygNbAsoSG5h3M7LLizuVJwblylLuJnM8+4sWnFvLCW7vwwVe9MKtGu1Y/cswRqznqhJb06VuHmjXjDtRtq5RKCgCSriNUH20GPgfOBOoDzxKSw/fAIDNbXtx5PCk4lyRbfmHxjIm8NDqHF8e24r2Zffhlcy0a1lvL4Qcs5OiBjThiYAYZGXEH6rZGyiWF8uJJwbkKkLuJNfMmM+GlObw2th6vT+nDopU7I+XSs1MORx9tHHVCS7rsU92rmSoJTwrOuXKTu3I2n0/4lNdf3cRrH+zFp9+FdoiWzZZx1KErOer4neg3oB51/XaIlOVJwTmXHL+s5MdpE3nzpcW8NmEnxk3ry5oNDdiu5i/07vYj/Q+vS///aco++0D16nEH6/J4UnDOJV/uFjYu/ISJr33J2HHVGT91H6Z/3wWAHbZfyyF9VtP/qMYcengtdt015ljTnCcF51zFWzOXxdPeYcKbSxn/YXPGT+/HDytaArBry+X0P+QX+h/dlIP71WCHHWKONc14UnDOxWvzOmzJR3w9eRrjxxtvT9mDd2cdzOoNDZFyyey0hEMOqcaB/TPo3Uds79NDJJUnBedcatm4jE057zPlnW95+91ajJ/ajSnf9mDTlu2oVi2Xrnsu58C+23Fgv4YccAA0bRp3wFWLJwXnXGpb+z3r5r7HJ+8u4P0PazExuxuTvunFhk11ANhrj1Uc2LcWBx1SmwMPhObNY463kvOk4JyrPMxg1Uw2fv8uWe8vYOLHdXl/1n58NLs3azY0AKBdm5854ADY/8AG9NpfdOiATyhUBp4UnHOVV+4mWPYpm394jy8mL+b9jxswcea+fPD1AaxYG1qoGzXcQM/MjezXpwG99q9Gz57QqFEJ501jnhScc1WH5cKqmeQu/pBvPvuWSZOMSTP3YNI3vcjO6YRZKDLs2W4NvXrXolfvmuy3H3Ts6KWJPJ4UnHNV29oFsPQjVs/7lCmT1jDps2ZMmrMfk+fsx/I1TQBo2GAT+3bfQmbP2mRmQvfu0KYNaTk0hycF51x6+WUV/DQZW/Ih30ybz6RPajF5djc++bYn2Tmd2LQ5zBexQ6Nf6N4dMntsR/fuIVHsskvVTxSeFJxz6S13Eyz/HJZNZuOiz5nx+WqyZjRl6tzuZM3NJDunE5u3hDHBm+ywie7dRea+NejeHXr0gJYtY46/nHlScM65gn5ZBcunwvJP2bDwC6Z/tpapX7Ug67tMps7tTnZOJ7bk1uCPx33FyId+hsZdoHrVmLLUk4JzzpXGhiWwLAuWf8r6hdOY/vl66vADnVvPgGo1oVFnaNwtJIhGXaBxZ6jZMO6oy8yTgnPObQ0zWLcAln0Kyz8Nzys+h19W/LpPvbZRkuj8a7Ko3xaUul2diksKNSo6GOecqzQkqNc6PFofH9aZwbocWDkNVk6HFdPC6x9eCd1lAWrUh0Z7R6WJKFE06gQ1G8T3XkrJk4JzzpWFBPVahUeLo39dv3kdrJr5a5JYMQ3mj4Y5D/66T/3dQpLYfm/Yfk/YviM02COl2io8KTjnXHmoURea7Bseecxg3fchQSQmiwVjgKjqXtVCsmgYJYmGe4aE0XBPqFm/4t9GhV/ROefShQT1dgmPlsf8un7zelg9G1bNgp+//PV50Zuh+2yeuq0SEkVCwqjVJGkhe1JwzrmKVqNOqEZq3OW363M3wepvf5soVn0JSybClvW/7le7GbQ5FbrdUf6hlfsZnXPObZ1qNWH7DuHRauCv6y0X1n7/a6L4+Uuom5w76jwpOOdcqlM1qN8mPFocmdRLVXhHWkntJX2R8PhZ0kWSdpA0XtI30XPjio7NOefSXYUnBTP72sy6mllXoDuwDhgDXAFMMLN2wIRo2TnnXAWK+5a7fsC3ZjYfOBYYGa0fCRwXW1TOOZem4k4KQ4DR0esdzWwRQPTcrLADJJ0tKUtS1tKlSysoTOecSw+xJQVJ2wHHAM+V5TgzG2FmmWaWmZGRkZzgnHMuTcVZUjgC+MzMFkfLiyU1B4iel8QWmXPOpak4k8JQfq06AngFOC16fRrwcoVH5JxzaS6WpCCpLtAfeDFh9S1Af0nfRNtuiSM255xLZ5V6PgVJS4H523CKpsBP5RROefK4ysbjKhuPq2yqYly7mFmhjbKVOilsK0lZRU00ESePq2w8rrLxuMom3eKKu0uqc865FOJJwTnnXL50Twoj4g6gCB5X2XhcZeNxlU1axZXWbQrOOed+K91LCs455xJ4UnDOOZfPkwIg6RJJJqlp3LEASLpB0vRovolxknaOOyYASbdL+iqKbYykRnHHBCBpkKSZknIlxd51UNIASV9LmiMpZYaAl/SopCWSsuOOJY+kVpLelfRl9Du8MO6YACTVljRF0rQoruvijimRpOqSPpf0WnmfO+2TgqRWhDuov487lgS3m1nnaM6J14Br4g4oMh7oZGadgdnAlTHHkycb+AMwMe5AJFUH7ieM7dURGCqpY7xR5XscGBB3EAVsBv5uZnsC+wHnpcjPayNwiJl1AboCAyTtF3NMiS4EvkzGidM+KQD/Ai4DUqbF3cx+TlisR4rEZmbjzGxztDgZSM4ksWVkZl+a2ddxxxHpAcwxs+/M7BfgacJcIbEzs4nA8rjjSGRmi8zss+j1asIHXYt4owIL1kSLNaNHSvwfSmoJHAU8nIzzp3VSkHQM8IOZTYs7loIk3SRpAXAyqVNSSHQG8GbcQaSgFsCChOUcUuBDrjKQ1AbYB/gk3kiCqIrmC8KIzePNLCXiAu4ifJHNTcbJayTjpKlE0tvAToVsGg78L3BYxUYUFBeXmb1sZsOB4ZKuBM4H/pEKcUX7DCcU+0dVREyljStFqJB1KfENM5VJqg+8AFxUoKQcGzPbAnSN2s7GSOpkZrG2x0g6GlhiZlMl9U3GNap8UjCzQwtbL2lvoC0wTRKEqpDPJPUwsx/jiqsQTwGvU0FJoaS4JJ0GHA30swq8yaUMP6+45QCtEpZbAgtjiqVSkFSTkBBGmdmLJe1f0cxspaT3CO0xcTfS9waOkXQkUBtoKOlJMzulvC6QttVHZjbDzJqZWRsza0P4Z+5WEQmhJJLaJSweA3wVVyyJJA0ALgeOMbN1cceToj4F2klqG80uOIQwV4grhMI3skeAL83szrjjySMpI693naQ6wKGkwP+hmV1pZi2jz6whwDvlmRAgjZNCirtFUrak6YTqrZTopgfcBzQAxkfdZR+MOyAASQMl5QC9gNcljY0rlqgh/nxgLKHR9FkzmxlXPIkkjQYmAe0l5UgaFndMhG++pwKHRH9TX0TfguPWHHg3+h/8lNCmUO7dP1ORD3PhnHMun5cUnHPO5fOk4JxzLp8nBeecc/k8KTjnnMvnScE551w+TwoOKHoEzbKMQCrpb5I2SNq+mH1M0n8TlmtIWloeoz1G9wZ8IukbSc9E9wkU3KdvwWtJelzSCdt6/TLEeYakGdFos9mSjo3WXy8paTfpRe9zbtTtc5qkfsm6VhHXv1bSJdHrpL5Xt/U8Kbg8j1P4CJplGYF0KKFP98Bi9lkLdIpuCIIwQu0PpQ+zWLcC/zKzdsAKoEL74UsqcYSAaDCz4UCfaLTZ/YDpAGZ2jZm9ndwouTQaffciIGn3mUSjxRapgt6r2wqeFBxQ9AiapR2BVNJuQH3gKkJyKM6bhFEeifYdXbZoC72+gEOA56NVI4HjtuI8/aJx6mdEpada0fp5iubbkJQZDXuQ9+13hKRxwBOS9orG4f8iKgm0K3CJZsBqYA2Ama0xs7nRufJLLNH1rpP0WRRLh2h9fUmPJZQ0jo/WHyZpUrT/c9FYQsWZRMJAfZK6S3pf0lRJYyU1j9bvLuntqGTxmaTdFNwelXJmSBoc7dtXYW6Ep4AZ0brhCnNLvA20T7head5rhqTx0fr/SJqvFJnzpCrzpODKS96H+weEO2abFbPv08AQSbWBzhQxKqak9vr1LteCj4IT/DQBViYM7V3c6KQHJJ6LMJQIUTyPA4PNbG/C2GB/KemNA92BY83sJODPwN3Rt/HMKI5E04DFwNzow/1/ijnvT2bWDfg3cEm07mpglZntHZU03ok+KK8CDo32zwIuLiHmAcBL0fuuCdwLnGBm3YFHgZui/UYB90fzCuwPLCKUHLsCXQjDP9yel0QIQ4cPN7OOkroThmLYJzpm3zK+138QhnHoBowBWpfwnlw5qPID4rkKMwQYaGa5kl4EBhEmm/kdM5uuMEzyUOCNok4YlVC6lvL6ZRmd9AMzOzr/QOnx6GV7YK6ZzY6WRwLnEYYqLs4rZrY+ej2JMLptS+BFM/vmNwGZbVEYQ2pfoB/wL0ndzezaQs6bNzjcVMKHKoQP4SEJ51uhMHJmR+CjUGBiuyiOwtwu6TZCiSVv0pj2QCfC8CUA1YFFkhoALcxsTHStDQCS+gCjo1FEF0t6P3o/PwNT8ko+wAHAmLxxsiQVNwZUYe+1D1FVpJm9JWlFMce7cuJJwW0zSZ2Bdvz6obId8B1FJIXIK8AdQF/Ct/zCztseeKaI4/ua2cqE5Z+ARpJqRKWFrRmdtLDEkmczv5asaxfYtjbvhZk9JekTQvXYWElnmtk7iTtHo8tOAaZIGg88BlxbyDU3Rs9b+PV/Vfw+2YkwNk9J1XYAlxI+gC8gJL3u0fEzzazXb04qNSziHMX9nNYWWC7tODpFvVdXwbz6yJWHocC1eSPOmtnOQAtJuxRzzKPA9WY2o6gdzOxrM+taxGNlgX0NeBfI60V0GlDWeRa+AtpI2j1aPhV4P3o9j/ABCnB8USeQtCvwnZndQ0h8nQts31lSt4RVXYH5ZYhxHGHAvbzzNSbMgtc7L25JdSXtUdQJzCwXuBuoJulw4GsgQ1Kv6PiakvaK5jXIkXRctL6WpLqETgeDFSahyQAOJCS5giYCAyXViUodxVWVFeZD4MTo2ocBjct4vNsKnhQcUPQImirdCKRDCHW+icYQ2g12lvS7KiIzyzGzu8v3XXA5cLGkOYTSxyNlOTiqHjkdeE7SDMLMVnk9dK4D7pb0AeHbbFEGA9lRW0UH4IkC22sCd0j6KtpnMGUbBfdGoHHUyDsNONjMlgJ/AkYrjOo5Obp2kaIkeiNwWTRt6AnArdE5vyC0H0BIjBdE5/2YMNHRGEKPqWnAO9E5fjfkfDTN5jPR+V4gtDeVxXXAYZI+I8x5vYjQSO+SyEdJdc6lJIWeX1vMbHNUivl31IDvksjbFJxzqao18KykasAvwFkxx5MWvKTgnHMun7cpOOecy+dJwTnnXD5PCs455/J5UnDOOZfPk4Jzzrl8/w+YIOanUAi8xAAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "# 4b. \n", | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "Ti=85\n", | |
| "Tf=98.6\n", | |
| "Ta=65\n", | |
| "K=0.275\n", | |
| "N=10\n", | |
| "\n", | |
| "#Analytical Model\n", | |
| "t=np.linspace(0,4,N)\n", | |
| "Tanalytical=Ta+(Ti-Ta)*np.exp(-K*t)\n", | |
| "\n", | |
| "# Times & Temps Before Death\n", | |
| "#0,-1,-2,-3 : 65 60 58 55 \n", | |
| "BDt=np.array([0,-1,-2,-3]) #before time death\n", | |
| "numericalBD=np.linspace(0,-len(BDt),N)\n", | |
| "dtBD=numericalBD[0]-numericalBD[1] #dt before death \n", | |
| "BDT=np.array([65,60,58,55]) #Before death temp\n", | |
| "TBD=np.zeros(N) #Placeholder zeros\n", | |
| "TBD[0]=Ti \n", | |
| "for i in range(1,len(numericalBD)):\n", | |
| " for j in range(0,3):\n", | |
| " TBD[i]=TBD[i-1]+K*(TBD[i-1]-BDT[j])*dtBD\n", | |
| " \n", | |
| "#Times & Temps After Death\n", | |
| "#0,1,2 : 65,66,67\n", | |
| "ADt=np.array([0,1,2]) #time after death \n", | |
| "numericalAD=np.linspace(0,len(ADt),N)\n", | |
| "dtAD=numericalAD[0]-numericalAD[1] #dt after death\n", | |
| "ADT=np.array([65,66,67]) #temperature after death\n", | |
| "TAD=np.zeros(N) #placeholder zeros\n", | |
| "TAD[0]=Ti\n", | |
| "for i in range(1,len(numericalAD)):\n", | |
| " for j in range (0,2):\n", | |
| " TAD[i]=TAD[i-1]+K*(TAD[i-1]-ADT[j])*dtAD\n", | |
| "\n", | |
| "plt.plot(t,Tanalytical,'-',color='orange')\n", | |
| "plt.plot(numericalBD,TBD,'-',color='red')\n", | |
| "plt.plot(numericalAD,TAD,'-',color='blue')\n", | |
| "plt.title('Temperature of Body Over Time')\n", | |
| "plt.xlabel('11 A.M = 0 Hours Since Recording')\n", | |
| "plt.ylabel('Temp(F)')\n", | |
| "plt.legend(loc='best')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "In comparison to the previous linear model, this non-linear one makes a lot more sense as it follows the Newtons Law of Cooling exponential decay model. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 177, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Death was 1.8865228851460631 hours ago\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "T=98.6\n", | |
| "K=0.275\n", | |
| "Ti=85\n", | |
| "Ta=65\n", | |
| "dt=np.log((T-Ta)/(Ti-Ta))/-K\n", | |
| "print('Death was',-dt,\"hours ago\")" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Since the first temperature recorded on the body was at 11 am, subtracting 1.88 hours yields the time of death to be at 9:07 A.M. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 3", | |
| "language": "python", | |
| "name": "python3" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython3", | |
| "version": "3.7.3" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 4 | |
| } |