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-project.ipynb
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
404 lines (404 sloc)
89.4 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": 88, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"dT/dt = -5.5 degrees per hour\n", | |
"K = 0.275 hours^-1\n" | |
] | |
} | |
], | |
"source": [ | |
"T2 = 74\n", | |
"T1 = 85\n", | |
"delta_t = 2\n", | |
"T_ambient = 65\n", | |
"coeff_k = -(T1-T_ambient)\n", | |
"\n", | |
"\n", | |
"dT_dt = (T2-T1)/delta_t\n", | |
"print(\"dT/dt = \", dT_dt, \"degrees per hour\") # rate of temperature loss per hour\n", | |
"\n", | |
"K = dT_dt/coeff_k\n", | |
"print(\"K = \", K, \"hours^-1\")" | |
] | |
}, | |
{ | |
"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": 89, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"K = 0.275 hours^-1\n" | |
] | |
} | |
], | |
"source": [ | |
"T1 = 85\n", | |
"T2 = 74\n", | |
"T_ambient = 65\n", | |
"delta_t = 2 #elapsed time\n", | |
"\n", | |
"def measure_K(T1, T2, T_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: temp_t1,temp_t2,temp_ambient,delta_t \n", | |
" Returns: K'''\n", | |
" dT_dt = (T2-T1)/delta_t\n", | |
" K = -dT_dt/(T1-T_ambient)\n", | |
" return K\n", | |
"\n", | |
"result = measure_K(T1, T2, T_ambient, delta_t)\n", | |
"\n", | |
"print(\"K = \", result, \"hours^-1\")" | |
] | |
}, | |
{ | |
"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": 113, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"t = -1.8865228851460631 hours\n", | |
"Time of Death is 9:06.78 AM\n" | |
] | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAg4AAAEdCAYAAABkEe+RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeZxN9f/A8dd7doyxjrHvS5ZsYxsSlax9KyKyy5ItCWlRqFREaZGUfoW00MK33bdECYVkTYqiMbJmG7uZz++Pz5nrzjUrY85w38/H4z5m7uece8777O/zOZ9zjhhjUEoppZTKiAC3A1BKKaXUlUMTB6WUUkplmCYOSimllMowTRyUUkoplWGaOCillFIqwzRxUEoppVSG+X3iICK9ReSc23Eod4nIUhF5IwuGM15EtmVFTF7DzJLYLmH8RkS6uzX+7CAiO0Tk0SwcXllnvl2XVcNUVx4R2SMio9yOI6tlKHEQkUIi8qyIbBWRUyKyT0S+F5GeIhJ0uYO8zOYBJdwOwmtHk/Q5JyKxIjJdRPJmcliPisiOyxRqthKRUBE5ICInRaSw2/EkEZHrnOVU1qfTFKBR9kfkHif5NiLys4gE+HRzNenJhPrA1OwamTNfTDqfstkVjxtEpJ+InMoBcVyTxjIYmsFhzBWRr1LodC0wPWsjTnH8LZx4i17ucQGke9AXkZLAcuAcMBb4BTgLNAZGARuAdZcxxstCRAQIMsacBE66HY+X24BV2GVTA/g/IBTo62ZQLroD2AnsBXpjD8w5ljEmHoh3Ow4XGKA60BOY5W4oGSciIcaYM8aY/dk86g5AiNf3XcBI7IlMkuyOKUskzdNsHmewMebsJQ6mNbDep+zopQzQhfUqexhj0vwAnwJ7gHwpdAsG8nj9PxGIA84AvwJdffo3wL3YjeM48DfQEcgHvAMcA/4E7vD6TVnndz2AxdiD/F9AN59hPwVsAU4AscAM75ixB51zwA3Y5OcMcEtSuVd/EcBbzjSfdob1vM80Z2Q6BwNvO9MUC4xOZz4nTed1PuXPAet9yqKB/2EPUPuBj4EyXtNpfD7jgX7ArhTGN9errA/2AC3O9yjsQWC/Mx3Lget9YqkIfAQcBg45cV2bwnxvAqx1ls9qIDq9dc/5/XfAMOBOYGsK3ZcCbwCPOcvsXyfmPF791AW+BPY582w10Dql4XjNh8NAbp9+xmHXvXIpzOOlTj/jgW0+v2sBLHOm/YgzTRUuJrY05lNG1/80lwN2+9gAnHL+3uBMX/c0xp007KTtIk9qsac0LcCjwA6v77OAb7D7il3OfHkDu+0NxCaSh4DXgRCfYd0L/ObE/wcwBnuCkNR9BzABexZ4EFjtVf6oV39B2BOl7dj9QBzwslf3+7AnTPHY9e59oFh623Ma8/Ac0DuVbj28lslfwLNALq/uPzrTMwk44MybcUAg8KSzbu0FxvkMd4/T3yzs9r3f6V+8+glx1q2d2H3vJqCPV/cwZzoHAfOd4cx1uk12lsUJ7L7+ZSDc6daaC7ehGV7TM80n1gnAb17f3wc+wyZbO4FEZ5kJcD/wuzO/tgKjgcA05v01zvjrpdFPKPCSsx6cBnYDs51uE1OYli5e83iUzzwfi12fjzrfBwC5gFex+51YoL/P+JNO0o87454LFPGJ3/vzVSbWnxuAldh1+Sj2+HhDmutrOitzQSABrw0qjX4nYzfETkBl4BFnYd7k1Y9xZlQv7AFnurNSfYnd+VR0Vq7jQCGfDXA30A2o4qxEid4LGrvzaer0fxN2hZ3ts3NLxO4sbwTKA5FcmDi8hM06GwKlsTUr/S9iOvcC/YEK2AOfSWthkMKOxolxM/CqV1k1ZwE/7qww1wIfYDeUMOwKONFZ+Yo6n3DOH+yqOMPpi92h7PYa9tvA+87/ubBJ0UdAPWfZjMFuNFWdfqKc5fmqE0cVZ/kdBCJ95vv3zvK5BptcbMdrh57GBn0aKIzdcP8Fmvv0sxS7sU11+m/tfH/cq5/m2HWumrPMJmCTvso+w3nDa9oPAb28ugdgDy5jsDvkW535Wd+ZxwWd/sbjlThgk4YE4AWglhNjX+Cai4ktjXmV0fU/1eUAFMdue2858dyM3eFkNHGIwK73j6cWe0rTQsqJwxFgNlDVmdengC+AOU5st2APZIO8fjceexBpj13f22IPWE969bMDu3Mc78zval7l3onDbOz20QO7DTcC7vfqfp+zbMsBMcAK4Lu0tud0ll+KiQM2UToAdMXuD27Abpczvfr50ZlfTwKVsAdx48yvCc50DsBnH4Tddo9ik+7K2IT5JDDYq5/3sYnmTc60dsUmB92c7kmJw34n1gpARafbWGyiWhZoCWwDXnO6hQAjnOWatJ+K8JqejCQOx7DJSk3stiXYfd+fzjpTDvgP9tgxJp39THqJwyPYg+712ONCA+Bep1s48CHwrde0hHnNY9/E4V9sglvRWWYJzrJKKhvvrA8VvX43AnvcKufM01XAIqdbIPZ4ZLD74aJAgYysP9j96jFnvlV01oM7gJg019d0VuYGTjAd0ukvN3YHP9infAHwrdd3A7zg9T3SKfPO5As4Zbf4bIBP+gx7BV5nyynE1N6JKcBr52aApint9Ly+/xeYlQXT+ZJPP78Bz6QRb9J0nsAmBqec79/iZOleO9X3fX4b6vzudud7sh2xV39/JcWOreF5HLvjSNp57gIGeM2XXfgc3J14XnD+Hw/86NNdsAej4T7zva5XP43wSmLSmCdTgQVe36cD7/r0sxTY4FM2A1iZzrDX47Uz4cID3EvAD17fW2Ev0RVzvl/nTENZn+GOJ3nisAz4LK1YMhtbBoeR2vqf6nLA7px3kvwM/RYymDg4/9+DTT5KpjJfL5gW3/UVu47vw6s2AfgcuwMM9dlWP/TaNk9wYW1NT+Cw1/cdwOIUpmEHTuKA3YEaoGMm5ncd5zclfLbni04csNvSPymUt8Qmgbmd7z8CP/n0sx2nNsWrbCswwev7HuBrn36eT1p/OX9ALefTz9M42z3nE4dXMjCNdwHHvL73A06l0F9GE4cDJD9zzodd55v7/HYAsCeNuJKmM2nf6/2p4/TzGvYEV1IZxly8zvJ95rFv4vC+1/cg7L7+A6+yQOw21C+NmGOcmJNOsFs434tmZv0Bijm/a5TRdd0Yk27jSHH+mnT6q4jNIL/3Kf8Oe93Tm+cakrHXfxKwZzVJZYewZ1xFfH630uf7cuyZhw1UpIPTYHO3iMRjD4wh2OzL2+p0pmU60FFENonIiyLSxqvBV2am07fdRxz2DD09fYDa2Cy6FTab/cQrhvpAexGJT/pgz/DDsGcbaVmCzVrBZp6LsAe2G0WkCraR6Lde4ykKHPYZV1Ov8dQHon26H8PuNL1jMSS/dhjn/E11fohIGHanP9ureBbQQUQK+fSe5rwWkUinkelvInLYibM6UCa18WN3FE1EJGkd6w98boz5J43fpCTpslKKLjK2lIaTkfU/veVQDVhljPG+y+iHzMSBrYL9E3twuRRbTPLr5Huwl6pO+5Ql7SeqY2uKPvJZH18D8olIpNfvVqUz7rrO37SWW3MRWeQ0YD7G+fmUqeWWjpLY5TfdZ5oWYPfNFbz69d0G9uC1X/Uqy8h+tbyz/dV3yjb6jH8EF+5rLpinItJZRH4QkX+c370JhItIwdQmOJM2GttGLUlN7Dr/uU+8LwJRGWhk3hW77/X+/Op0ewN7Iv27s722F5Hgi4zb+xh4Drv/9j4GJmCTIs+ycho/fu21vn3jdEprfUt3/XH2Z3OBpSLyuYiMFpGK6U1Aeo0j/8BmJtWdkaXHN8GQFMpSasDiW2ZI/44P8fwj0hBbXf8M8AC2mrkR9qDj3QApwRiTZiteY8wiESmNPWg3x87UjSJyk098vrH4lvk2DsrINAHEGWOSbuf73VlJVjixfOsM421s1ZKvg+kM+1vgRRGpDuTFbuzfYqshE4BYr3EHYK+Zt09hOCe8+lkMpNTy+IjX/4nOxpAkaV6lNT86Yi+VfWjbsXoEYqv2n/cqS29ez8JWL47G1rqcxJ6xhJAKY8xmEfkB6CciE7FVn7enEW9a0kq8Mx2br0ys/+kth5TW4/ROGpL3bEyCiIwEvhKRl1LoJRGvbdeR0g44pX1CWvuJpL+dsJftfP3r9f/xFLpnmLN/+AK7HT6B3cmXxO7MM7zcMiBpmgZi9wG+Yr3+z+z8So33sgng/OU432El+nxPNk9F5HrgXez8GYG9fNgMp11KOjFkdB3xXY5J03YrtuYsvf597fLa/yVjjFnt3OXSEnvS9QowXkQaG2Myuz5lalk5B/LPsMnLOOx+vgK2Fi6teZmh9ccY00NEJmOn7WZggogMMMbMSm3AaSYOxph/ReRLYKiIvGyM8T4Y4GRcIdhrV6exK8Zmr16u9/l+KRphN9YkMdgDG9hq4wPGGM992CLS8WJHZIz5F3gPeE9E3sJm5dXInun0lXT2l9v5uwabWW83Tr1TCs5gD7C+FmMPxvcD3xtjzonIt9jr9omcr21IGk9P4KgxZl8q41mDraaO88n8s8I92IPqcz7lPbBn/8/7/iAN12Mbp34CICJ5sNf7NqXzu9ewbRP+xZ6ted9ulZSspDSfvf2MTUJfzuLYvGXV+r8Z6CEigV4JRqafQ2CM+Z9za9rzXHiA2YdtS+GtLpduM7bKt7wx5ov0ek7HWudvS+y1a1/1sbUbw5PWexGJvsRxpiQWO78qG2PmXIbhw4W3DscAfxljTonIGuwBvIQx5psLf5qmptgD8eNJBXLhs0BS209d7DqyAXsALmeMWZyJWDPEGHMM2+brIxGZgk30GwNfk/q0ZIWG2MRpeFJtoIg08eknpf1RhtcfY8wG7PybIiKzsPvYWan1n5Ez4MHYhfGziHQVkWoiUtFZCdYAlYwxJ7DXhJ8UkU4iUklEHsHeWnipVZZJ+jrjrywiT2BX8BecbluBSBHpKyLlRaSnE3emichTTrVvFRGphG2QGQ/8nU3TWVBEiopIcRGJwe5893E+Y3wa22Bsrog0EJFyInKD2Msq5Z1+/gKKikiMiBQWkdwATrXUVuwZe1KSsA67c7+V5InDO85wPheRlmKfM9FQRB4WkaQz72nYFXWhiDR1+rnOmYeNL3YGOJcHrgPeNMZs8v5gD+bXOGc0GbUV6CYi14pIbWxSmJGNPOmg8Rjwf8YY74NgUkvutiJSRETypTKMJ4E2IvKCiNR01qveYi8NXUpsvtOXFev/q9h2R6+LSFWnlu2pixgO2NbuMdjqXW/fAC1E5E5nP/IQ9iBzSYy9DfZp4GkRGerM5+oi0kVEJmVyWNuw6/90EekuIhVEpL6I3Of08gf2jHCks/3djm0ImKWc9e1RYJSIPOTse69x9k+vZNFoGorIGGdf1gvbsHKqM/7N2FqDWc6+t4KI1Bb7/IWR6Qx3K1BCRHo46+Td2DYN3v4CgkSkrbOfyuOUf4PdZjo4cT3G+csmqXIuc0/GHvwGOseKGk7sGVmPCzn7Xu9PXgBnv3eXswzKYy8pn8WeTCZNS3VnuyksIllZ8/Q79lh9v7O+3QE87NPPDudvO2d/FJGR9ccpe1pEmohIGbEJSQznL9GkLCMNIbA7k+c4f4vLPux1/e6cb5Gd0dsUu/uUpdQo6BROwxCS34651Om2A+jh85snsS26j2NrJu7Cq/EaPo0gvX6XrBx7kNiETRaSbp3zvtPhYqfzG1JpdOkznUmfROxZ7gKghk+/12Ibhh3CVm1vw1YBFvSK8V3smbIBxnv99hWnrI5X2UdOWUmf8RTCHkySpjXOicf7t2WwO9n92NqYndjLO+VSm+/Yal2DTyMmr+4vOuNKrSHSas7f8rWU9BvbXYtNvE46685g3+WR0nCc8qnYyzglU+g22okzgbRvx2yFrbU66axTS7BnxpcU26Wu/yktB+xlq43OstyEbRNzwfqc1jbkVT7N+a1348hgbMK/D1t9/Qq2Ott7ec0CvvEZ1htJ89irbAZeDVidsr7YZPgUdvv4ieR3XuwghbvEfMudOJ90ys9gGwp7N+wegj2jO4lt35B0e2Fzn+05K27H7OhMx0lsY+a1wMNe3VNqTPgDzu2Nqa1HnL8d823s/u4ANvkK8OonCFsj+bszH/Zj19+khthJjSM7+oxLsLeH7nfWyU+w+3BD8gZ8051+DOdvxwxx1ov9zjJ8AbvPveB2zFTm10Ds2fNp7D5wJWk3NEzpdsakzxSnn3uxtykecz4/AW29hlEE22bsKOnfjjnKZ/y7gIfSWR/vd/o76SzHpEbLjbz6eRR7B0kiyW/HTHX9AUoBCzl/m2kcdp+fN631Nel+/RxL7HWlv7B3Q2S2oZZSl0RE5mNbbv/H7ViUykoisgd7YMzRD1VTOc+V/rhopS4LESmArUJvj20wpJRSCk0clErNL9hLNc8aY5a6HItSSuUYOf5ShVJKKaVyDr9/rbZSSimlMk4vVbiscOHCpmzZsm6HoZRSV5Sff/75gDEmMv0+VVbTxMFlZcuWZc2aNW6HoZRSVxQRSenpkCob6KUKpZRSSmWYJg5KKaWUyjBNHJRSSimVYZo4KKWUUirDtHGkUipFiYmJ7Nq1i+PHL+kt1EpdlODgYIoUKUJERITboSgfmjgopVJ04MABRIQqVaoQEKCVkyr7GGM4efIkcXFxAJo85DC6N7hSrVwJo0dDQoLbkair1OHDh4mKitKkQWU7ESF37tyUKFGCffv2uR2O8qE1Dlei2Fho3x727oWNG+G99yB/frejUleZhIQEgoOD3Q5D+bFcuXJx9uxZt8NQPvRU4kr0yis2aQD46ito1Aj++MPdmNRVSUTcDkH5MV3/ciZNHK5ETz8NY8ac/751KzRoAF9/7V5MSiml/IImDleigACYMMFeoggLs2WHD0ObNvDSS6BvPFXqoo0fP57u3btf1G+XLVtGlSpVLjmGsmXL8s0331zycNITHh7On3/+ednHo64umjhcybp0gWXLoEQJ+z0hAe67DwYMgDNn3I1NqWzQvHlzChQowOnTp10Zv4iwbds2z/emTZuydevWyzrO3r17ExISQnh4uOdTq1atdH/XvHlz3njjjWRl8fHxlC9fPstjnDVrFtddd12WD1flDJo4XOnq1YPVq+2liiRvvAEtWsD+/e7FpdRltmPHDpYtW4aI8Mknn7gdTrYaPXo08fHxns/69evdDkn5EU0crgbFisF334F39eqyZVC/PmzY4F5cSl1Gc+bMoVGjRvTu3ZvZs2cn69a7d2+GDBlCu3btyJs3Lw0bNmT79u2e7vfddx+lSpUiIiKC6Oholi1bluI42rVrx8svv5ysrGbNmixcuJDrr78egFq1ahEeHs68efNYunQpJUuW9PQbGxtLhw4diIyMpFChQgwdOhSA7du3c+ONN1KoUCEKFy5Mt27dOHz48CXPk1OnTtG9e3cKFSpE/vz5qV+/Pnv37mXMmDEsW7aMoUOHEh4e7onDu8akd+/eDB48mDZt2hAeHk6TJk3Ys2cPw4cPp0CBAlxzzTX88ssvnnFNnDiRChUqkDdvXqpVq8aCBQsA2LJlCwMHDmTlypWEh4eT37nj6/Tp04waNYrSpUsTFRXFwIEDOXny5CVPs8p+mjhcLcLCYM4cmDQJkloi79wJjRuDs0ErdUlEsu+TAXPmzKFbt25069aNRYsWsTfpTiPHe++9x7hx4zh06BAVK1ZkjFeD4vr167Nu3Tr+/fdfunbtSqdOnTh16tQF4+jVqxdz5871fF+/fj1xcXG0bduW77//3lMWHx9P586dk/02ISGBW265hTJlyrBjxw7i4uLo0qULYB9w9PDDD7N79262bNlCbGws48ePz9B0p2X27NkcOXKE2NhYDh48yIwZM8iVKxdPPfUUTZs2Zdq0acTHxzNt2rQUfz9//nwmTJjAgQMHCA0NJSYmhrp163LgwAE6duzIiBEjPP1WqFCBZcuWceTIEcaNG0f37t35559/qFq1KjNmzCAmJob4+HhPQvTggw/y+++/s27dOrZt20ZcXBxPPPHEJU+zyn6aOFxNROxDoT75BPLmtWXHj0OHDrYxpTaaVFeJH374gZ07d3LnnXcSHR1NhQoVePfdd5P106FDBxo0aEBQUBDdunVj3bp1nm5JZ+VBQUGMHDmS06dPp9g24bbbbuOPP/7gD+d257fffpvOnTsTEhKSboyrVq1i9+7dTJ48mTx58hAWFua57l+xYkVuvvlmQkNDiYyMZMSIEXz33XcZnv4pU6aQP39+z6dXr16AfUzzwYMH2bZtG4GBgURHR2fqqYvt27cnOjqasLAw2rdvT1hYGD179iQwMJDOnTsnq3Ho1KkTxYsXJyAggM6dO1OpUiVWrVqV4nCNMcycOZOpU6dSsGBB8ubNyyOPPML777+f4dhUzqGJw9XollvskyW9Gz099hjcdRecOOFeXEplkdmzZ9OyZUsKFy4MQNeuXS+4XFG0aFHP/7lz5yY+Pt7z/bnnnqNq1arky5eP/Pnzc+TIEQ4cOHDBeEJDQ7nzzjuZO3cuiYmJvPfee/To0SNDMcbGxlKmTBmCgi58zt6+ffvo0qULJUqUICIigu7du6c4/tSMGjWKw4cPez5J096jRw9atWpFly5dKF68OKNHj87UA5SioqI8/+fKleuC797zcM6cOdSuXduTvGzatCnVadi/fz8nTpwgOjra03/r1q3Zr+2wrkiaOFytqleHVavghhvOl82bB9dfD7t2uReXunIZk32fNJw8eZL58+fz3XffUbRoUYoWLcrUqVNZv359hhoJLlu2jEmTJjF//nwOHTrE4cOHyZcvHyaV8fbq1Yt33nmHxYsXkzt3bmJiYjI0u0qVKsXff//NuXPnLuj28MMPIyJs2LCBo0ePMnfu3FTHnxnBwcGMGzeOX3/9lRUrVvDZZ58xZ84cIGsfprRz50769+/PtGnTOHjwIIcPH6ZGjRqeafAdV+HChcmVKxebN2/2JDtHjhxJloioK4cmDlezQoVg0SIYPPh82c8/Q716nFmRcmMwpXK6hQsXEhgYyK+//sq6detYt24dW7ZsoWnTpp6DZFqOHTtGUFAQkZGRnDt3jieeeIKjR4+m2n9MTAwBAQGMHDnygtqGqKioVJ+D0KBBA4oVK8ZDDz3E8ePHOXXqFMuXL/fEkNRwMC4ujsmTJ2diDqRuyZIlbNy4kYSEBCIiIggODiYwMDDdWDPr+PHjiAiRkZEAvPXWW2zatMnTPSoqil27dnHGuS08ICCA/v37c//993vePREXF8eiRYuyJB6VvTRxuNoFB9tHVL/6KjhVpstC91L5g+v58fWxLgenVObNnj2bPn36ULp0aU+NQ9GiRRk6dCjvvPNOimf43lq1akWbNm2oXLkyZcqUISwsjFKlSqX5m549e7Jx48YLHgw1fvx4evXqRf78+Zk/f36yboGBgXz66ads27aN0qVLU7JkSebNmwfAuHHjWLt2Lfny5aNdu3Z06NAhU/Pg2WefTfYch6RLNnv27KFjx45ERERQtWpVmjVr5on5vvvu48MPP6RAgQIMGzYsU+PzVa1aNUaOHElMTAxRUVFs3LiRJk2aeLrfeOONVK9enaJFi3pimzRpEhUrVqRRo0ZERETQokWLy/7MC3V5SFZUj6mLV69ePbNmzZrsGdmSJfzduz31Oh9hfx4IOQevnb6Z3k99bhMMpbxs2bKFqlWruh1GjjBnzhxef/11fvjhB7dD8TuprYci8rMxpp4LIfk9rXHwJzfcwK55M0kIslWXZ4KgT56vGT64HOf2/uNycErlTCdOnGD69OkMGDDA7VCUyhE0cfAzjRt1YvWgtdQ4ef4WrRdLxtFqbDkO/viti5EplfMsWrSIyMhIoqKi6Nq1q9vhKJUjaOLgh8qXqsnK8bG053z137fFT1N/Xgs2vjnRxciUyllatWrF8ePH+e9//5vibZVK+SNNHPxUeFgEH47dxOPFzp9F/ZXfELP9YT4e/R9Ip4GZUkop/6SJgx8LkADGDniHBddNI/ysXRWOh8AdeT5j/IBKJO7f53KESimlchpNHBS33zSElXcvp/zp3J6yx8vs4I4Hy3Js9XIXI1NKKZXTaOKgAKhRvhGrx+ykBecfU72wzEliZjdl++wXXIxMKaVUTqKJQyaJSFkR+UJEDonIHhGZJiJBPv2MExEjIi3civNiFMxTmC8f28r9kbd6yjZHGupvuZ9vHrpT2z0opZTSxOEiTAf2AcWA2kAzwPNMZxGpAHQErsgHIwQFBPH84P8yq8EzhCbY580fygWtQj/ghb7VMZl4EY9S/mrZsmVUqVLlkodTtmxZvvnmmyyIKG3h4eFZ9jhqdfXTxCHzygHzjTGnjDF7gK+A6l7dpwEPAmfcCC6r9GrzEN91/ZpiZ0IBSAyA+8v/Tp/7y3FqbcqvzlUqu5QtW5aoqCiOHz/uKXvjjTdo3ry5e0F5adq06WV/nHLv3r0JCQlJ9ujpWrVqpfu75s2b88YbbyQri4+Pp7z323SzyKxZszyvEldXD00cMu9FoIuI5BaREkAbbPKAiHQCzhhjvkhrACIyQETWiMianPxa2YbX3MSa0dtoaEp4ymZXjKfZazHsnvuqi5EpBefOnePFF190O4wLpPeujKw0evRo4uPjPZ+MvB1UqUuliUPmfYetYTgK7ALWAAtFJBx4Ghie3gCMMa8bY+oZY+olvV0upyqeryRLH91G7wI3espWFU8ket1glo/uAmfPuhid8mcPPPAAU6ZM4fDhwxd027FjByKS7CDufaY9a9YsmjRpwv3330/+/PkpX748K1asYNasWZQqVYoiRYowe/Zsz29Pnz7NqFGjKF26NFFRUQwcOJCTJ08CsHTpUkqWLMmkSZMoWrQoffr08ZQliY2NpUOHDkRGRlKoUCGGDh0KwPbt27nxxhspVKgQhQsXplu3bilOT2adOnWK7t27U6hQIfLnz0/9+vXZu3cvY8aMYdmyZQwdOpTw8HBPHCLCtm3bAFuTMXjwYNq0aUN4eDhNmjRhz549DB8+nAIFCnDNNdfwyy+/eMY1ceJEKlSoQN68ealWrRoLFiwA7DsmBg4cyMqVKz1vAnbnT60AACAASURBVE1vXqorgyYOmSAiAcAi4GMgD1AYKABMAh4H3jbG/OVehJdHWFAYb977DS/WfpjARFu2Jy80D5vHK72qYf65IptzqEySxyXbPhlRr149mjdvzpQpUy5qen766Sdq1qzJwYMH6dq1K126dGH16tVs27aNuXPnMnToUOLj4wF48MEH+f3331m3bh3btm0jLi6OJ554wjOsPXv28O+//7Jz505ef/31ZONJSEjglltuoUyZMuzYsYO4uDi6dOkCgDGGhx9+mN27d7NlyxZiY2MZP378RU2Pt9mzZ3PkyBFiY2M5ePAgM2bMIFeuXDz11FM0bdqUadOmER8fz7Rp01L8/fz585kwYQIHDhwgNDSUmJgY6taty4EDB+jYsSMjRozw9FuhQgWWLVvGkSNHGDduHN27d+eff/6hatWqzJgxg5iYGOLj4z0JUXrzUuV8mjhkTkGgFDDNGHPaGHMQeAtoC9wEDHPutNjj9DdfRB50L9ysIyIMu+1pFt2xkELn7Js0zwXC0Crb6DOyAie/1/dcqOz3xBNP8PLLL3Mxl/zKlStHnz59CAwMpHPnzsTGxjJ27FhCQ0Np2bIlISEhbNu2DWMMM2fOZOrUqRQsWJC8efPyyCOP8P7773uGFRAQwOOPP05oaCi5cuVKNp5Vq1axe/duJk+eTJ48eQgLC/Nc969YsSI333wzoaGhREZGMmLECL777rsMT8OUKVPInz+/59OrVy8AgoODOXjwINu2bSMwMJDo6GgiIiLSGdp57du3Jzo6mrCwMNq3b09YWBg9e/b0zCvvGodOnTpRvHhxAgIC6Ny5M5UqVWLVqpTbQWVkXqqcTx++ngnGmAMi8hcwSESmAOFAL2A9MAzwfjf1amAE8GW2B3oZ3VTzNn4uuZUOL1/HWnYDMLvKSTa+exMf/zyOMsPHgWTsjFGpS1WjRg1uueUWJk6cmOlXgEdFRXn+TzrY+5bFx8ezf/9+Tpw4QXR0tKebMYaEhATP98jISMLCwlIcT2xsLGXKlEnxXRf79u1j2LBhLFu2jGPHjpGYmEiBAgUyPA2jRo1iwoQJF5T36NGD2NhYunTpwuHDh+nevTtPPfUUwcHBKQzlQr7zIaX5kmTOnDk8//zz7NixA7ANLQ+kcvdVRualyvk0cci8DsAL2DsnEoAlwP1O7YOHiCQAh4wx8RcO4spWpmA5fnhkG4Nn3sas/V8DsLYYRO95nPf7f0+LFz+FPHlcjlJlNTPOuB1Cih5//HHq1q3LyJEjPWV5nPXvxIkTnjPtPXv2XNTwCxcuTK5cudi8eTMlSpRIsR9JI1kuVaoUf//9N+fOnbsgeXj44YcRETZs2EChQoVYuHChp93BpQgODmbcuHGMGzeOHTt20LZtW6pUqULfvn3TjDWzdu7cSf/+/Vm8eDExMTEEBgZSu3ZtjLHriu+4MjIvVc6nlyoyyRizzhjT3BhTwBhT2BjTyRhzwUsdjDFljTGX/wZsl+QKzsWbgxYxvdGTBCfancPB3NCqxBKe7VkBc5lvRVMqScWKFencuTMvvfSSpywyMpISJUowd+5cEhISePPNN9m+fftFDT8gIID+/ftz//33s2+f3dTj4uJYtGhRhn7foEEDihUrxkMPPcTx48c5deoUy5fbR7kfO3bM03AwLi6OyZMnX1SMvpYsWcLGjRtJSEggIiKC4OBgAgMDAVubkFXPbDh+/DgiQlIj77feeotNmzZ5ukdFRbFr1y7OnLF3p1/qvFQ5gyYO6qKJCINaPcrSHosplmDfc5EYAA/W3EunJ2pw7MN3XY5Q+YuxY8cme6YDwMyZM5k8eTKFChVi8+bNNG7c+KKHP2nSJCpWrEijRo2IiIigRYsWGX5OQ2BgIJ9++inbtm2jdOnSlCxZknnz5gEwbtw41q5dS758+WjXrh0dOnTIVFzPPvtssuc4FC5cGLC1Kx07diQiIoKqVavSrFkzunfvDsB9993Hhx9+SIECBRg2bFimxuerWrVqjBw5kpiYGKKioti4cSNNmjTxdL/xxhupXr06RYsW9cR2KfNS5QySVKWk3FGvXj2zZs0at8O4ZP8c+4dOr1zP8tPbPGVV98OCvP2o8uSrkML1XZWzbdmyJdPtBpTKaqmthyLyszGmngsh+T2tcVBZoljeYnz7wGaGlu3sKdsSCQ14g0/uqgv79BXdSil1NdDEQWWZkMAQXu71PrNbvEJYol21jobBbTU2MrZ/RRJ+XOFyhEoppS6VJg4qy/VsMpgVA1ZR1uT3lD1Z9xj/eeU6Dk1/DvTymFJKXbE0cVCXRZ0S0awZvY2WEXU9ZV9WNNTbOooNA26DEydcjE4ppdTF8pvEQUQCRKS+iPQTkVHO3/rOY6TVZVAodyG+uG8VD9cY5Cn7syA0ivqUtztVBm1JneNp42nlpsTERLdDUCm46g+aIpJPRCYAccCHQCegkfP3IyBORCaISP40BqMuUmBAIE/fMZ2Pbn+X8ER7Z8XJYOjZII6BD9Xg1LtzXI5QpSYsLIyDBw9q8qCynTGGM2fOEBcX53mYl8o5rvrbMUXkN2Au8K4x5oKnnohIOaAb0M0Yk+33nl0tt2NmxJZ9v9Jh5k38du78E/yid8OHwd0o+9z/QWioi9EpX2fPnmXXrl2cOnXK7VCUHwoKCiJfvnwULlyYgIALz3H1dkz3+EPikMsYk+47WzPaX1bzp8QB4NjpY/R/uxPz4s4/Ka7ASXh7Q0XaTVsE5cu7GJ1S6kqhiYN7rvpLFcaYkyISIiKVUuouIpVEJMSNpMEf5Q3Ny3t9v+TlGyZ7HlV9KBfc0nAbY4ZWI2Hhxy5HqJRSKi1XfeLgeAC4L5Vu9wIjU+mmLgMRYej1o/i+33JKBZx/E+DTDU/TcsEd7B01EM6edTFCpZRSqfGXxOEu4NlUuk3GtnFQ2axRqRjWjvydloUbesq+LQ91zWv8cFsdiI11MTqllFIp8ZfEoaQx5u+UOhhjYoFS2RyPchTOXZgvBi1nXIPRiNPcZncENK+/mef7XIP56it3A1RKKZWMvyQOp0QkxZe/i0hJQJuNuygwIJDxbSbxZbcvKIR9y2ZCAIxseoJOb7Xh6GOjISHB5SiVUkqB/yQOXwFPpNJtnNNduaxVpTasHb6FhhHVPGUfVYN6Ryaz8fYY2LMnjV8rpZTKDv6SOIwBWorIKhF5UER6On9/AloBj7gcn3KUzlea74f9wtAad3vK/igEDWutZnbna2DpUveCU0op5R+JgzEmDqgLfA90ASY4f78Hop3uKocICQzh5Tv+j/fav0MeggH7tMneNx6h/9QbOPXkOL10oZRSLrnqHwCV0/nbA6Aya8v+LdzxVmu2nDzftrXWHpgf24jKMz+GYsVcjE4p5RZ9AJR7/KLGQV25qkZWZdXwzXSt0N5Ttr4oRNf8kfc6VoFFi9L4tVJKqaymiYPK8cJDwpnb7SOmt36ZUBMIQHwodG15jAGvtObkQyP1gVFKKZVNNHFQVwQRYVDDoay8ZzUVw87fWTszGhoefZ7f2jaAnTtdjFAppfyDJg7qilKnWB1+Hv4rXSre7inbGAX1GqxjbpdqsGCBi9EppdTVTxMHQETap9+XyikiQiN4t+vHvNbuVUIJAuB4CPRofYK+szpwYug9oK+CVkqpy0ITB2uy2wGozBERBtQbyE/3rKFy7tKe8jfrQgPzOr+2rAO//+5ihEopdXXSxAEwxlR0OwZ1cWoVrcWaYZvoWqWjp2xzEajf9Ddm97gW3n7bxeiUUurqo4mDuuLlDc3L3M7zmXnL64Q5D4w6EQK9256h98c9Od6nG8THuxylUkpdHfzqAVAi8jWQ0gSfBnYBHxtjvs7OmPQBUFlr496N3Dn3Nn6L/8tTVnU/fLCqLNXf+C/UrOlidEqprKIPgHKPv9U4rAFqAOuBz5y/1YHfsMnDPBEZ7l546lJdG3Utq+/dQI9qd3nKtkRC/ZY7eLNfNOaVV8CPkmWllMpq/lbjsAK4xxiz0ausBvC6MaaxiDQG5mRnmwetcbg8jDHMWjeLIZ8O5KQ54ynvtgGmJ7YhYuYcKFzYxQiVUpdCaxzc4281DtWBrT5lfzjlGGNWAEWzOyiV9USEPnX6sHrgWqpGVPCUv1MT6pb+ktU3XQPffONihEopdWXyt8RhBfC6iJQUqyQwHVgJICLVgL1uBqiyVvUi1Vk9ZD29a/TwlG0vCI1vO8jkcTeTOPoBOHMmjSEopZTy5m+JQ08gP7ATOAfsAAoCvZzuAvRI8ZfqipUnJA9v3TGHdzq8Q96AXACcC4TRLaH1ninsuaG+PvNBKaUyyK8SB2PMfmPM7UBuoDyQxxjT3hiz1+m+2blcoa5CXa/tyrqhm2hQpI6n7OsKULPJBr7sUBPefFMbTiqlVDr8KnEAEJHywDBgmDHmtIhUFJGqbselskf5AuX5YcBPPNT4QQQBYH8eaNvpNCPn9eV05zvg0CGXo1RKqZzLrxIHEbkNWI1tDNnfKS4EvOhaUCrbBQcG88zNE/lfj/9RNPT8nRXPN4bGBRbwe9Nq8P33LkaolFI5l18lDsDTQEtjTG8gwSn7Baid0QGISFkR+UJEDonIHhGZJiJBItJIRL4WkX9FZL+IfCAixS7DNKgs0qJ8Czbc+ytty7fylK0tDnVv28Oc+5pjHnsUzp51MUKllMp5/C1xKGaM+dn5P+lidiI4ddYZMx3YBxTDJhzNgMFAAeB1oCxQBjgGvHXpIavLKTJPJJ91/5IXWr1AiJx/02av2w3df32Kozc0hj//dDlKpZTKOfwtcfhFRLr4lHXAPlEyo8oB840xp4wxe4CvgOrGmC+NMR8YY44aY04A04AmWRO2upxEhPsa3ceP/VdROV95T/m7NaFO9BpWtb4W5s51MUKllMo5/C1xuA94XkQWAblF5L/AFGBEJobxItBFRHKLSAmgDTZ58HU9sPlSA1bZp06xOvw8eD131+rjKfuzIDTpcoJnZ/QgsVtXOHzYxQiVUsp9fpU4GGM2AVWAd4CngAVATWPMlkwM5jts48qj2BdjrQEWevcgIjWBscADKQ1ARAaIyBoRWbN///5MT4e6fMJDwvm/29/kvTveIyIoD2Cf+fDgzdAy6D3iGlWDJUtcjlIppdzjV++quFQiEoB9aNRr2JqKcOBNYKsxZrTTT0VscvGQMebt9Iap76rIuf469BddP+jMj/+s9pQVPAEzP4UO7UbBhAkQGupihEr5L31XhXuu+sQhjVdpJ2OMaZmBYRUG9gP5jTFHnLLbgQnGmBoiUgabNEw0xszISHyaOORsZxPOMm7pOCb+MBHjtRrdvRZeiK1G3tnvw7XXuhihUv5JEwf3+MOlig+Bj5zPGqAmsBH7Wu0NwLXYZzukyxhzAPgLGOTcgpkf+7jq9U57h2+BVzKaNKicLzgwmKdvepolvZZQKk9xT/mbdaFO01/58da68NxzkJjoYpRKKZV9rvoaB28ishwYbIxZ71VWE3jVGJOhOyBEpDbwAlAL+yyIJcAQYBAwHjju3b8xJjyt4WmNw5Xj8KnDDP58MO9tes9TFpgIj30HYwKaETT7bShVysUIlfIfWuPgHn9LHI4AkcaYM15lYcAeY0x+N2LSxOHK886Gdxj82UCOno33lMXEwtyv81L+mRnQtauL0SnlHzRxcI8/XKrwthx4Q0RKOa/VLg28ivNabaUyolvNbqwfvJGmpa7zlK0sBbW6HWPW5G6Yu7ro+y6UUlctf0scegER2DsjzgF/Yl+z3du9kNSVqGz+sizpvZSnbnyKIAkEID4U+twOdybM49/6NWDxYpejVEqprOdXiUN6r9VWKjMCAwJ5pOkjrOz3I5ULVPSUf1gdat62m2/7t4ARI+DUKRejVEqprHXVJw4iUtC3zBhz2hiz0xhzOq3+lMqIesXrsXbgOgbUHeApi4uAm3rBqM1TOd2gLqxb52KESimVda76xAFYLiLPi0idlDqKSG0ReR74IZvjUleRPCF5eO0/r7Gw80IKh53PQZ9rDA2bbGFz2/rwzDNw7pyLUSql1KXzh8ShLrAbmC8ih0VktYh86/w9BLwPxDr9KXVJbrvmNjYO2UzrCq09ZeuLQnTfc0z95BESr2sCW7e6GKFSSl0af7sdsxL2gU/5gUPABmPMdjdj0tsxr07GGKatmsYD/xvF6UTP3b80/wtmfRlKmUcmwb33QoA/5O5KZT29HdM9fpU45ESaOFzdNu/bTPePu7Fur+eZY+Q9DS9+Cb3zN0PefAvKlXMxQqWuTJo4uEdPd5S6jKoXqc5P/VcxpukYApzN7Vgo3H073F70O/Y2qgEzZ4Im8EqpK4QmDkpdZiGBIUy4cQLL+y6nktdtm59cA9f2OsGC5wdAu3YQF+dilEoplTGaOCiVTRqVbMQvA9cxpP4QT9n+PNChC/QK/ZIjdavDO+9o7YNSKkfz28RBn9ug3JAnJA/T2k5jUfdFlAg//7bNObXh2m5HWPxYd+jUCfbvdzFKpZRKnV8lDiISISJvishx7C2YiMgtIvKYy6EpP9OyQks2Dt5Et2u7ecpi80GLXnDf8Y84UasaLFzoYoRKKZUyv0ocgFecvzWApHvkVgE93AlH+bMCuQowt8NcPuj0AYW8Hhr1UiOoe8cBVg9pDz17wuHDLkaplFLJ+Vvi0AoYaIz5CzAAxph9QJSrUSm/1rFaRzYN2cwtlW/xlG0tDDH9YFzs25y9tjp88YWLESql1Hn+ljgcAwp4F4hISUBfcqVcVTS8KJ90+YQ3/vMG4cHhACQEwBPNoVHb3Wzq0w5699bXdSulXOdvicMs7KOnYwBx3l/xJjDT1aiUAkSEvnX7sn7QepqWbuopX1scogfAM9tnc65GNfj0UxejVEr5O39LHJ4CvgTeAcKBj4FvgSluBqWUt/IFyrOk1xIm3zyZkIAQAM4EwSMtoHG7Pfza91bo3h0OHnQ5UqWUP/KbxEFEAoEuwFRjTHljTLAxppwxZqLR526rHCYwIJBRjUex9p611Ct+/qm6q0tAnXtg4s53bO3DggUuRqmU8kd+kzgYYxKAV40xp92ORamMql6kOiv7ruTpG59OVvvwcAtocss+fr2nA3Tpos99UEplG79JHBxfikgrt4NQKjOCAoJ4uOnD/HzPz0QXi/aUryoJde+BSbvm2dqHDz5wMUqllL/wt8QhAVgoIl+JyEwReT3p43ZgSqWnRpEa/NjvR5668SmCA4IBOB0ED90M1916gN8G3wkdO8JevUlIKXX5+Fvi8DfwAvALcAA46PVRKscLCgjikaaPsPaetclqH34qCbUHwuR/PiKhRjV47z1954VS6rIQbRfornr16pk1a9a4HYa6Ap1NOMuk5ZN44rsnOJt41lMeEwtvLYQq190Gr74KxYq5GKVSl4eI/GyMqZd+nyqr+VWNg4g0Tu3jdmxKZVZwYDCPXv8oawasoU7ROp7ylaVs7cNz+/5rax9mzdLaB6VUlvGrGgcR+cenqCAgwAFjTPEUfnLZaY2DygpnE84y8YeJPPn9k8lqHxr/DW/+F6rUvRleew3KlXMxSqWyjtY4uMevahyMMcW8P9iHQE0BnnY5NKUuSXBgMI81e4w1A9ZQu2htT/mK0lBrEEw8+TVna1aHqVMhIcHFSJVSVzq/Shx8GWPOAmOBMW7HolRWqBlVk1X9VvF488cJCggC7J0XD7eAht1O8suzI6BxY9i40eVIlVJXKr9OHBxN0+9FqStHcGAwY5uNZU3/NcnuvPilGNQfAGPCV3Gqfh147DE4rc9DU0pljl8lDiLyh4j87vXZBXyOrXVQ6qpSq2gtfuz3I5NaTCIsKAywb9x8+nqo0y+BFbMnQO3asHy5y5Eqpa4k/tY40vepkceBLcYY157joI0jVXb4/eDv9PukH8v+XuYpEwNDV8HTiyG8/xB45hnIm9fFKJXKOG0c6R6/qnEAKhljFnl9fjDGHBSRwW4HptTlVLlQZZb2Xsr0ttMJDwkHwAi83BBqDIavv3wFqlWDzz93OVKlVE7nb4lDandPTMjWKJRyQYAEMKj+IDYP3kzriq095TvzQ8uecHf0Lg51vAW6dtWXZimlUuUXiYPXQ54CRCTG5+FP3YF4t2NUKruUzleaL7p+wZzb51AwV0FP+Vt1oNoQWPDLe1C1Krz9tj44Sil1Ab9o4+D14KciwD6vTgbYAzxtjPkw2wND2zgod+2N38uwr4Yxf/P8ZOUdN8PLX0LRJq1g+nQoX96lCJVKmbZxcI9f1Dh4PfDpI5+HQBU3xtR1K2lQym1R4VHM6ziPBZ0XUDS8qKf8w+q29mH2nkWY6tVg4kQ4ezaNISml/IVfJA5JjDF3uh2DUjnR7dfczq+Df+Xu2nd7yg7lgt7tocWdp9k2+WGoWxdWrHAxSqVUTuBXiYOI5BGRp0VkuYhs9X6mQyaGUVZEvhCRQyKyR0SmiUiQ0+0mEflNRE6IyBIRKXP5pkaprFUgVwH+77b/4+seX1M2f1lP+bfl4dpB8Ez+TZxt2gQGDYLDh90LVCnlKr9KHIBXgObAS0Bx4DHgEDAzE8OYjm0nUQyoDTQDBotIYeBjZ5gFgTXAvKwKXKns0qJ8CzYN2sTImJEEiN1FnAqGR1pA9D3w42cz4JprYN48bTyplB/yt8ShDXC7MWYekOD87QRk5hJGOWC+MeaUMWYP8BVQHegAbDbGfGCMOQWMB2qJyDVZOgVKZYM8IXmY0nIKq/qtom6xup7yjVHQuC/cW3cvR3t1gbZt4a+/XIxUKZXd/C1xCAKSnhIZLyIRwC6gSiaG8SLQRURyi0gJbDKSlDysT+rJGHMc2O6UJyMiA0RkjYis2a/3y6scLLp4ND/1+4nnWj5H7uDcgH1w1LSGtvHkwh1fQfXq8Oyz2nhSKT/hb4nDBs6/1GoF8ILz2Z6JYXyHTQaOYpOONcBC7Cu6j/j0ewS44Bm+xpjXjTH1jDH1IiMjMzUBSmW3oIAgRsSMYPPgzbSp2MZTHhcB7btAh1tPEjfhQYiOhh9/dDFSpVR28LfEYSCQ9EyH+4AQoAzQOyM/FpEAYBG2LUMeoDBQAJiEfYhUhM9PIoBjlxq0UjlB2fxl+bzr57x/x/sUyVPEU76gKlQdCtPDNpLYOAYGD4Yjvjm0Uupq4TeJg4gEAncAOwCMMf8YY7obY24zxqxP88fnFQRKAdOMMaedl2O9BbQFNgO1vMaXB6jglCt1VRAROtfozJYhW+hXp5+n/FgoDGkH190Nmz561TaenD9fG08qdRXym8TBGJMAjALOXMIwDgB/AYNEJEhE8gO9sG0bFgA1ROQOEQnDvqp7gzHmt0uPXqmcpWCugsy8dSZLey2lcqHKnvKVpaDOPfBotT2c6tYZWreGP/5wMVKlVFbzm8TB8S7Q5xKH0QFoDewHtgHngPuNMfuxNRpPYW/xbAh0ucRxKZWjNSvbjPUD1zP2+rEEBwQDcC4QnrrePvvhf9v/BzVqwLhxcPKky9EqpbKCX7yrIomILMY2jvwTiMW+qwIAY0xLN2LSd1Woq8Wv+39lwKcDWB67PFl5503w/CIoHlkepk2DNm1SGYJSGafvqnCPvyUO96TWzRjzWnbGkkQTB3U1STSJzPx5Jg9+8yBHTp9vIJn3NDy1GAavhsD2HeCFF6BUKRcjVVc6TRzc41eJQ06kiYO6Gu2N38uor0cxd8PcZOV1d8OMz6D+kTz28sXw4RAc7FKU6kqmiYN7/K2NAyLSw3nXxGrnexMRae92XEpdTaLCo3i7/dss7rk4WePJtcWhYX8Y0uw4h8eOhtq14fvvXYxUKZVZfpU4iMhYYDgwH6jkFO8FxrgWlFJXsRvL3ciGgRt48oYnCQsKA+yTJ6c3gGuGwruBv2KaNYNevWDvXpejVUplhF8lDkBfoK0xZhbnG0ZuB8q7FpFSV7nQoFAevf5RNg3aROuKrT3le8Oh2x1wc0/Y+vkc++yHV1+FhAQXo1VKpcffEocQIOl9wEmJQx7ghDvhKOU/KhSswBddv+CDTh9QPG9xT/ni8lBzEIytc5iTwwZDo0awerWLkSql0uJvicP/gEnOUySTjAW+cCkepfyKiNCxWke2DNnC8IbDPa/tPhMETzaDawfDokNroGFDuOceOHDA5YiVUr78LXEYDlTGvnwqQkQOAdcCo12NSik/ExEawdTWU1nTfw0NSzT0lG8vCK17wJ0dDbvefx0qV4bp0/XyhVI5iF8lDsaYQ8aYtkA14AagjjGmjTHmcDo/VUpdBnWK1WFF3xXMaDeD/GH5PeUfVLeNJydVO8SZYUOgXj1YvjyNISmlsotfJQ4AIhIO1AdqA9EicsFrr5VS2SdAArin3j1sHbqVHjV7eMqPh8BDN9v2D18fXQfXXQc9e8I//6QxNKXU5eZXiYOINAV2Yts13Aw8BuxwypVSLiqSpwhz2s9haa+l1ChSw1O+tTC07AmdOkHsf9+GKlVgyhQ4c9Hvq1NKXQK/ShyAV4HhxphaxphbjTG1gfuAGS7HpZRyNCvbjLUD1jK11VQiQiM85R86ly+eqXWM0w89ALVqwddfuxipUv7J3xKHUtg3ZHp7DyjpQixKqVQEBwYzvNHwCy5fnAiBR1rYuy++OvcbtGwJHTvCzp0uRquUf/G3xOE9oJ9PWV8uTCaUUjlA0fCizGk/h2V9llEzqqan/I9C0KY7dOgMO7/5CKpWhSefhFOnXIxWKf/gVy+5cl6rfT3wN7ALW9NQCvgeSEzqLztfsa0vuVIqY84lnuPV1a/y6JJHOXr6qKc811l4ZBmMWgFhpcvD1Knwn/+AiIvRqstNX3LlHn9LHFJ9rba37HzFtiYOSmXO3vi9PPjNg8xePztZeYV/4aUvoe0f2EsYU6dCtWruBKkuO00c3ONXiUNOpImDUhdnRewKhnwxhHV71iUrv/U3mLoIyh8NhCFD7Ou7CxZ0KUp1uWji6hHC8AAAFpJJREFU4B6/SxxEpD5QBwj3LjfGPO9GPJo4KHXxEhITmLFmBo8ueZTDp84/xy3kHIxcaS9hhIcXtO0fBgyAoCAXo1VZSRMH9/hV40gReQ5YBNwKxHh9GrkZl1Lq4gQGBDKkwRC2Dt3K3bXv9pSfCYJnmkKVoTC35L+YIUOgTh1YvNjFaJW6OvhVjYOI/AvUMsbEuh1LEq1xUCrr/LTrJ4Z9NYxVcauSlcfE2vYP9XYDt99uHyBVoYI7QaosoTUO7vGrGgcgDjjmdhBKqcujYcmGrOy7klm3zaJoeFFP+cpS0KA/9L0V9n690DaafOghOKa7A6Uyy98Sh/7ATBG5RUQaeH/cDkwplTUCJIBetXuxdehWRjceTXBAMABG4M26UGkYTKl3hjNTJkGlSvDWW5CYmM5QlVJJ/O1SRR/gZefrSa9OxhhTxIWQ9FKFUpfZHwf/YMT/RvDZ758lK698wN590fYPIDoaXnwRmjRxJ0iVaXqpwj3+VuPwLNDFGBNujIn0+riSNCilLr9KhSrx6V2f8mW3L6lSqIqn/PfC0K4btOsKv+/42b5986674O+/XYxWqZzP3xKH09i7KpRSfqZ1xdZsHLSR51s+n+zlWV9UhhqD4YGb4eiC9+3bN8eM0fYPSqXC3xKHscCzIpLf7UCUUtkvODCY+2Pu5497/6BfnX4I9rHUZwNhShOodC/8X9VTJDzztG3/8PrrcO6cy1ErlbP4WxuHs0AgYICEpGJsG4cQN2LSNg7/3969h0dVnXsc/74kQIAEwQqIIBdBBEEugoiVAEUONCKiBWuFw7WIqSQBL9iCxyM9IlivxYI3igcsoEdb2oOKFC13tIhSQBHCVcQLeC0GFBFY/WNvhkkawgxMZk8yv8/z7IeZtdas/WZnkbzZe+21RYLz9sdvM3rhaFbtXlWovM0eePiv0H0n0KoVPPSQt4y1JAzNcQhOsp1xaAW0AC4ELvK3Vv6/IpJk2p/TnhXDVjDnJ3Ool1EvVL7+bLhiCPT9GWz55F3o1QuysmDjxgCjFUkMSZU4OOfyj23AZ0Xei0gSMjMGXDSA/Jx87u56N1VSq4Tq5jeHljfDLb3gy2ULoXVryM6GvXsDjFgkWEmVOJhZdTN72swOALv9sqvM7K6AQxORgFWrVI0J3SawJXcLg9sMDpUfToHfXubNf3j0kqN8P/1Jb/7DfffBwYMBRiwSjKRKHIBp/r+tgEP+6zeBQcGEIyKJpn71+sy6ZhZrblxD5wadQ+VfVoXRWXDRzfBS3QLcuHHeHRjPPgtJNFdMJNkSh15AtnNuJ94ESZxznwJ1Ao1KRBJOh3M6sHzocl647gUa12gcKs8/C/oMgJ6DYMN3H8CAAdCpE6xaVUJvIuVHsiUOBUDN8AIzqw/ogqWI/Bszo/+F/Xlv1Hvc3+P+Qus/vNYE2t0EI/vA3o1vegtI9e8PW7cGGLFI6Uu2xGEm8LyZXQaYmbUDngamBxqViCS0tNQ0xl4+lq25W8lun00F8350Hq0A09tD0zyY3BkO/v+fvAdo5eTAp58GHLVI6Ui2xOFe4BVgDpAOzAMWAw8GGZSIlA21q9Xm8aseZ332eno2Ob6uw/7KML4HNM+BOS0Oc/Sxad5ju++5Bw4cCDBikdhLigWgzOwG59yzQcdRHC0AJVI2OedYuG0hty26jU2fbypUd/HH8MCr/gJSdevCr38Nw4ZBamowwZZDWgAqOMlyxuHJWHRiZvuLbEfM7Hdh9T81s01mVmBm75nZNbHYr4gkHjMj6/ws1mevZ2rWVM6qelaobu053gJSvQfAxsOfwMiR3hoQ8+frDgwp85IlcbBYdOI/VTPdOZeOdyfGt8ALAGZWD5gN3ApUB8YCc81MT94UKccqplRkVMdRbMvdxrjO40hLTQvVLWgGrX8BN/aBT3Zvgr59oWtXWL06wIhFTk+yXKr4BuhNCQmEc25xlH0OAe4GmjjnnJldCrwY/ohuM/sMuNo598aJ+tGlCpHy5cOvP+SuJXcxa90sHMd/vlY9BLe/7m0Zh4DrroNJk6Bp0+CCLcN0qSI4yZI4HAF2ceLEwTnnzouyz8XAcufcBP99Ct5Ey4eAl4E+wFTgAufcgSKfHQmMBGjQoEH7Xbt2RbNrESkD1u9Zzx2v3cGi7YsKldfZDxOWwoi1kFoh1VvC+q67oLZOTkZDiUNwkiVx+No5V/3kLSPurwGwE2jqLyZ1rPznwBQgDW9lyuuccy+X1JfOOIiUb4u2L2Lsq2PZsHdDofLmn8FvXoM++WAZGXDHHTBmDKSnBxRp2aLEITjJMsch1gYDK4skDT2A+4FuQCWgK/B7M2sbSIQikhB6NunJ2pFrmdl3JvWr1w+Vb64FfW+AbkNhTUaBd9ahaVOYOhUOHTpxhyIBS5bEISaTI8MMBmYVKWuLd+niLefcUefcGmA10CPG+xaRMialQgpD2g5hS84WJnWfREaljFDd8kbQcST8rD9s+34v5OZC8+YwezYcORJc0CInkBSJg3Mu4+StImNmPwTq4d9NEWYNkHnsDIO/KmUmsAEREaBKxSqMyxzH9rzt5HbM9eY4+P6vFbQYBaOuhD2f7YRBg6BdO3jpJd3CKQklKRKHGBsCzHPOFYQXOueWAROAP5pZAfAnYJJzbtG/dyEiyaxWtVo8mvUo7938Hv1a9AuVH06BxzpCkzz4r+6wb8s70KcPZGbCypUBRixyXFJMjkxkmhwpIn//8O/86rVfsWzXskLlP/gGxq+Am9dA2mGgd2+4915o0yaYQBOIJkcGR2ccREQC1ql+J5YMWcKCAQtoU+d4UvBFVbitFzTLhZlt4ciCl73LFwMHwvbtAUYsyUyJg4hIAji2hPXam9Yy+9rZNK7ROFS3+wwYdg20+QXMb+Zwc+d6EyhHjYI9ewKMWpKREgcRkQRSwSowsPVANuds5tEfP0qtqrVCdRtre7dwdh4OK845DI895j2Fc/x4+OqrAKOWZKLEQUQkAVVKqUTupblsz9vOhK4TSK90fGGo1xtAl+Fw1QB4J/0bmDwZGjeGiROhoKCEXkVOnxIHEZEEllE5g7u73c32vO3kdcyjYoWKobqXm3mXLwZfCzsq7PMWkWrcGB58EL79NsCopTxT4iAiUgbUrlabKVlTyM/JZ1DrQZi/rp0z+EMbuCAHftEbPjr0BYwd613CmDYNvvsu4MilvFHiICJShjSu2Zhnrn2Gddnr6H1+71D54RR44hJomge394TP930COTnQrBnMmAGHDwcYtZQnShxERMqg1nVa89KAl1g+dDmZDTJD5QcrwkM/hMaj4b9/BPv2fgAjRkCLFjB3rpaxltOmxEFEpAzLbJjJsqHLWDhwIe3rtg+V768M93SFxmPgvs5wYNc2b/2HNm1g3jwtYy2nTImDiEgZZ2b0atqLNTeuYd5P59GyVstQ3VdVYFwPaDIaftcRvtu8Efr1gw4d4JVXlEBI1JQ4iIiUE2bGtS2uZX32emZfO5smNZuE6vamQ96V3iqUM9rB4XVr4coroXNn+NvflEBIxJQ4iIiUMykVUhjYeiCbRm3iqaueon71+qG6D2rAiL5w4Sh4rhUcfeN16NEDunaFxYuVQMhJKXEQESmnKqZU5Mb2N7I1dyuP9Hqk0CqUW38AN/SHdjfBX5qDW7ECrrgCunWDJUuCC1oSnhIHEZFyLi01jTGdxrBj9A7u7X4vNdJqhOo2nA3X/gza3wQvNgO3fDl07+4lEEuXBhazJC4lDiIiSSK9UjrjM8ezc/RO7sy8k2oVq4Xq/lEXrh4AHW+El88Ht2wZ/OhH3rZ8eYBRS6JR4iAikmRqpNVgYveJ7Bi9g9svu50qqVVCdW/Vg6sGQqcRsLApuKVLvfkP3bvDihXBBS0JQ4mDiEiSql2tNg/0fIAdo3dwS6dbSEtNC9W9WR+y/hMu/zm8eh64JUugSxdvHsTKlQFGLUFT4iAikuTOTj+bh3s9zI68HeR1zKNySuVQ3RvnQs/B0GUYLGmEd+dFZqZ3J8aqVYHFLMFR4iAiIgDUzajLlKwpbM/bzqhLRlEppVKobmVD6D4Uug2FZQ3x1n7o3Bl69oTNm4MKWQKgxEFERAqpV70eU6+cyrbcbWS3zy70KO9ljaDbMLhiMKxsgHfrZlraibqSckiJg4iIFOvcM87l8aseZ0vuFka0G0FqhdRQ3eLzIHM49Li9NsvtgwCjlHhT4iAiIiVqVKMR06+eTn5OPsPaDiPFUkJ1f0v7mK4zu9JtZjeW7FyC08qT5Z4SBxERich5Nc/j6b5PszlnM4PbDKaCHf8VsmzXMro/050uM7vw6vZXlUCUY0ocREQkKk3PbMqsa2aRn5PP8LbDC13CWPnBSnrO7sllMy5jwdYFSiDKISUOIiJySpqe2ZQZfWewJWcLIy8eWWgS5eqPVtN7bm8umX4J8/PnK4EoR5Q4iIjIaWlcszFP9nmSbXnbuLnDzYVu43z7k7fp+1xfLn7qYuZtmsdRdzTASCUWlDiIiEhMNDijAdN6TwstJBW+EuW6Pevo93w/2jzRhuc3Ps+Ro0cCjFROhxIHERGJqXrV6zElawo78nZwa6dbCz0L491P3+X6P17PRY9fxNx35iqBKIOUOIiISKmom1GXh3o9xPtj3ueXl/+y0NM4N32+iYHzBjJpxaQAI5RTocRBRERKVe1qtbmvx328P+Z97sy8k4xKGQBUTqnMiItHBBydRCv15E1ERERO31lVz2Ji94ncdtltTFk9hcNHD1M3o27QYUmUlDiIiEhc1axSkwndJgQdhpwiXaoQERGRiClxEBERkYgpcRAREZGIKXEQERGRiClxEBERkYgpcRAREZGIKXEQERGRiJkedRosM/sM2HWKHz8L+DyG4cSK4oqO4opeosamuKJzOnE1dM7VimUwEhklDmWYmb3lnOsQdBxFKa7oKK7oJWpsiis6iRqXlEyXKkRERCRiShxEREQkYkocyrangg7gBBRXdBRX9BI1NsUVnUSNS0qgOQ4iIiISMZ1xEBERkYgpcRAREZGIKXEQERGRiClxSDBmdqaZ/dnMDpjZLjMbUELbW8xsj5ntM7OnzazyqfQTy7jMbIiZvW1mX5vZh2Z2v5mlhtUvNbODZrbf3/LjFNdQMzsStt/9ZtYt2n5KIa4nisT0nZkVhNXH+njlmNlb/n5mnqRtPMdXRHEFML4ijSve4yvSuOI2vsysspnN8L++AjP7h5llldA+buNLYkuJQ+KZBhwC6gADgcfNrGXRRmbWC/gVcAXQCDgP+HW0/cQ6LqAqMAZvRbhL/fhuL9ImxzmX7m8XnEZM0cQF8EbYftOdc0tPsZ+YxeWcyw6PCXgWeKFIs1ger4+BicDTJTUKYHxFFBfxH1+RxgXxHV8RxRXn8ZUK7Aa6AmcAdwHPm1mjog0DGF8SS845bQmyAdXw/rM0Cyv7A3BfMW3nApPC3l8B7Im2n1jHVcxnbwVeDHu/FBgRwPEaCqyM9dcXy/78zxUAXUvjeBXZ10RgZgn1cRtf0cQVz/EV5fGK2/g61eMVz/EV1v8GoF+ijC9tsdl0xiGxNAOOOOe2hJWtB4rLtFv6deHt6pjZD6LsJ9ZxFdUF2FikbLKZfW5mq8JP58Yhrnb+freY2V1hp7gT5Xj1Az4Dlhcpj9XxikY8x9fpKM3xFa14ja9TFdfxZWZ18L72ot8fKDvjS4qhxCGxpAP7ipTtAzIiaHvsdUaU/cQ6rhAzGwZ0AB4MK/4l3mnJeniLv7xoZk3iENdyoBVQG+8H6A3A2FPoJ9ZxhRsCPOP8P7F8sTxe0Yjn+DolcRhf0Yjn+DpVcRtfZlYRmAPMcs5tLqZJwo8vOTElDollP1C9SFl1vNOLJ2t77HVBlP3EOi4AzOwa4D4gyzkXevqdc261c67AOfedc24WsAq4srTjcs7tcM7tdM4ddc69A/wP0D/afmId1zFmdi7eteFnisQdy+MVjXiOr6jFaXxFLM7jK2rxHF9mVgHv0sIhIOcEzRJ6fEnJlDgkli1AqpmdH1bWhuJP9W3068Lb7XXOfRFlP7GOCzP7MTAd6OP/EC2JAywecZWw30CPl28w8LpzbsdJ+j6d4xWNeI6vqMRxfJ2O0hxfpyIu48vMDJiBN6mxn3Pu+xM0TdjxJREIepKFtsIb8BzezOdqwOV4p+haFtPux8Ae4EKgJrCYsMlDkfZTCnF1B74AuhRTVwPoBaThzcAeCBwALohDXFlAHf91c+Bd4O6gj1dY+3xgeByOV6rf32S8vwrTgNQEGF+RxhXv8RVpXPEeXxHFFcD4egL4O5B+knZxHV/aYrsFHoC2It8QOBP4i/8f+ANggF/eAO8UXoOwtrcCe4Gvgf8FKp+sn9KOC1gCHPbLjm2v+HW1gDV4pxz/6f+A+Y84xfWgf6wOADvwTiVXDPp4+WWX+e0yivRRGsdrAt5fleHbhAQYXxHFFcD4ijSueI+vaL6PcRlfQEM/joNFvj8Dgx5f2mK76SFXIiIiEjHNcRAREZGIKXEQERGRiClxEBERkYgpcRAREZGIKXEQERGRiClxEBERkYgpcRBJcGY23sx+H8f9rTKzdv7rCWY2O177LiGmh80sO+g4RMRbLUxEAmRm+8PeVgW+A474729yzk2KYyx9gALn3D/itc8IPQC8aWZPO+cOBR2MSDLTGQeRgDnn0o9teKvk9QkrmxPncLLxljCOq7BHUBfLOfcJsBm4Oj4RiciJKHEQSXDhlwvMrJGZOTMbZma7zewrM8s2s0vMbIOZ/dPMphb5/HAz2+S3/auZNTzBfirhPQtiWZGqSmb2jJkVmNlGM+sQ9pkWZrbU3+9GM7s6rG6pmY0Iez/UzFaGvXdmNsrMtgJbzfOImX1qZvv8r6dVWBxLgd5RH0ARiSklDiJl06XA+cD1wG+BO4EeQEvgp2bWFUKPnx4P/ATv2QQr8B4eVJzzgaPOuQ+LlF+N99ChGsB8YKrfd0XgRWARUBvIBeaY2QVRfB3X+F/LhUBPoAvQzN/X9XgPtDpmE4WfqCgiAVDiIFI23eOcO+icW4T3IKBnnXOfOuc+wksO2vntbgImO+c2OecOA5OAtic461AD74FHRa10zi1wzh3Bu4xx7Jd3JyAd76mGh5xzi4GXgBui+DomO+e+dM59C3wPZOA9XdL8mD8Ja1vgxygiAVLiIFI27Q17/W0x79P91w2BKf6lhH8CXwIG1Cumz6/wfnEXtSfs9TdAmj8n4Rxgt3PuaFj9rhP0fSK7j73wE4+pwDRgr5k9ZWbVw9pm4D3FUUQCpMRBpHzbjXdnRo2wrYpz7vVi2m4FzMwi/cX/MXCumYX/HGkAfOS/PoB3l8gxZxfTR6HH8zrnHnXOtce75NIMGBtW3QJYH2FsIlJKlDiIlG9PAOPMrCWAmZ1hZtcV19A59z3wGtA1wr5X4yUHd5hZRTPrBvTBmw8BsA74iZlVNbOmwM9L6syf4HmpP3fiAHCQ47el4sf1SoSxiUgpUeIgUo455/4M/AZ4zsy+Bt4Fskr4yJPAoAj7PoQ3cTIL+Bx4DBjsnNvsN3kEOIR3GWUWcLJbS6sD0/EumezCmxj5IICZ1cWbQPmXSGITkdJjzrmTtxKRpOHfMpmbSItAmdlDwHbn3GNBxyKS7JQ4iIiISMR0qUJEREQipsRBREREIqbEQURERCKmxEFEREQipsRBREREIqbEQURERCKmxEFEREQi9i/DLMU+N3K2tQAAAABJRU5ErkJggg==\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"import numpy as np\n", | |
"step = .1\n", | |
"tmax = 2 #adjusting this variable will show temp. of corpse after two hours \n", | |
"current_temp = 85\n", | |
"analytical_array = np.zeros(int(tmax/step))\n", | |
"numerical_array = np.zeros(int(tmax/step))\n", | |
"xaxis = np.zeros(int(tmax/step))\n", | |
"xaxis_em = np.zeros(int(tmax/step))\n", | |
"j=0\n", | |
"\n", | |
"def analytical_sol(t): #function allows you to get population at any year\n", | |
" current_temp = T_ambient + (85 - T_ambient)* np.exp(-0.275*(t))\n", | |
" return current_temp\n", | |
"\n", | |
"for i in np.arange(0, tmax, step):\n", | |
" # x axis\n", | |
" xaxis[j] = i\n", | |
" #Analytical Estimate\n", | |
" analytical_array[j] = analytical_sol(i)\n", | |
" #Estimate based on Euler's Method\n", | |
" slope = -0.275 * (current_temp-65)\n", | |
" #print(\"slope = \", slope)\n", | |
" numerical_array[j] = slope * step + current_temp\n", | |
" current_temp = numerical_array[j]\n", | |
" xaxis_em[j] = j*step + step\n", | |
" j = j + 1\n", | |
" \n", | |
"import matplotlib.pyplot as plt\n", | |
"%matplotlib inline\n", | |
"\n", | |
"plt.rcParams.update({'font.size': 12}) #font size\n", | |
"plt.rcParams['lines.linewidth'] = 3 #line width\n", | |
"\n", | |
"#Plot \n", | |
"plt.plot(xaxis, analytical_array, color='red', linestyle='-', label='Analytical Estimate')\n", | |
"plt.plot(xaxis_em, numerical_array, color='green', linestyle='-', label='Numerical Estimate')\n", | |
"\n", | |
"plt.xlabel('Time (hours)')\n", | |
"plt.ylabel('Temperature (deg. C)')\n", | |
"plt.title('Comparison Between Analytical and Numerical Temperature Estimates')\n", | |
"plt.legend(loc='best'); \n", | |
"\n", | |
"\n", | |
"'''By increasing the step variable in the initial conditions, you will see the \n", | |
"two curves diverge. When the step variable is decreased, they converge.\n", | |
"\n", | |
"As t->infinity, the final temperature will approach Ta (65 deg. C)\n", | |
"By adjusting the variable t max in the initial conditions of this program,\n", | |
"you will see that both solutions approach this value of temp\n", | |
"'''\n", | |
"\n", | |
"#Solving for time at which t = 98.6 deg. C\n", | |
"time_since_death = (np.log((98.6-65)/(85-65)))/(-0.275)\n", | |
"time_of_death = '9:06.78 AM'\n", | |
"print(\"t = \", time_since_death, \"hours\")\n", | |
"print(\"Time of Death is \", time_of_death)\n" | |
] | |
}, | |
{ | |
"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": 206, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"Based on your current t value of -2 , the ambient temperature is 58 deg. C\n" | |
] | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEdCAYAAAAM1BBYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZgU1dXH8e9h30RARlBAcAfjghEFTcQNxegbNbhAJAqJRom7SVxxA3fFLUFFFBWDQnDDRBNFXIKKiBh3QVEi+yooyDbAnPePW033DLM0zHTXTPfv8zz90HWruuoUDXW67r11r7k7IiIiteIOQEREqgclBBERAZQQREQkooQgIiKAEoKIiESUEEREBFBCEJEKmNk4M3sm7jgk85QQZKuZmVfw+jbuGDPNzEaZ2csxxzA5je+idZwxSs1QJ+4ApEbbIeX9QcAL0Z9zorKNWY+oiphZPXcvrCHHPA6oF72vDcwFBhC+j4TFlQxP8oDuEGSrufvCxAtYFhUvSSlfAuFCZ2Y3m9ksM1tjZp+Z2W8T+zGzBtGv2AFm9qyZrTKzb83sRDNrYWZjzexHM/vazH6Z8rmO0ef6mNmbZrY22ubU1DjNbMfol/xSM1thZm+Z2SEp64+N9tPTzN41s3XAmWZWYGajzWxOFPd0M7sw5XO3AX2Bnim/xPuknM8pJeJ428yGpSwvNLPrzWy4mS0DXovKm5rZ/Wa2IPq7mJp63qV8D8tSvodFUfEPqd+PuxeZWS0zuy76HgrN7Csz+32JGFtHVUSrzWy+mV1V8njR9/KOmS2PXhPMbN+U9S+Y2dhSPvdfM7u3rPOQasDd9dKr0i/g54ADHUpZNwb4L3AUsDNwOrAS6ButbxB9dj7hArs7MCLa5mXgN8BuwHDgB2Db6HMdo8/NBXoDewJ3EO5M9om2aQLMAEYDP432PQhYA+wabXNstJ/PCb+2dwZ2BHYC/gzsH5X1jz7365R9PwO8DrSOXg1SzueUEn8PbwPDUpYXAiuAgVFcnQg/0t4BJgCHALsC5wPrgZ+n8T3UiY7dp5R1A6O/0zOi410a7ffklG3eAD4DDgX2BcZFMT6Tss2vgZOi72Tf6O92HtA4Wv8LYB1QkPKZA6K4fhL3v1W9yvn3E3cAeuXGq6yEkHLR3rlE+S3A5Oh94gJ6W8r6dlHZnSllO0RlPUrse2CJfX8APBK9HwDMBGqV2GZS4ngpCeHUNM7zIeCfKcujgJdLbLMlCeGlEtscC6xKXFxTyp8CxqQRX6kJATBCMr26RPljwH+j912izx6Ysr4JsDw1IZRyzPpAIXB8tFwL+B/w55RtHgTeifvfqV7lv9SGIJl2YPTnp2aWWl6HcOFL9XHK+4XRn5+UUrZ9ic+9W2J5EuEXaeL4OwErShy/PrCgxOempC6YWR3gSuA0oE30mXrAdKrOlBLLBwINgUUl4q0HfFqJ47QDmgITS5T/B+hj4WB7ES7sUxMr3f1HM0v9XjCzjoS7rAOBloQEUBdoH32myMweAc4GhphZI8Jd4cWViF+yQAlBMq0W0a9OQvVEqqISy6nrvWSZu3t0kayo7Sv1SloL+AjoU8p2JRNSyeWrCNUqlxIS00pCgvh5BcdPnJeVKK+bRgy1CA3ApR1jXQXHTUfJ4Y2tjPelMrPawHjgC+AcQlXRekIyr5ey6QjgBjPrDuwSHXezdgWpXpQQJNOmEi40bdx9QoaO0Y1Qj59wMCEJJI7fC1jm7stKfrAC3QnVQyMTBWa2e4ltCgk9ezZx90Iz+4HQDpH4XCNCG8eHFRxzKuEOyNx9xhbGW545hLaAwwhtFAndgS+iZPs54aLeBXg/irsJsB9RgzfQgXC30dvd34222QtolHowd19oZi8AvyckhCfdfXUVno9kgBKCZJS7f25mTwGPm9nlwHvANoSLzrbuflcVHOYPZvY1IQn8DugMJHoxjQQuAl4ys2uBrwmNvz2AD939pXL2+yXwKzM7lPCr/axo36lVTf8DjjWzTsASYIWHrqMTgAvMbDKwGrie9Hr1/ZvQ1vAPM7uCUE20HeGO4Xt3fzyNfWwmuuDfDlxpZnOAycDxhAb7PtE2U81sIjDCzP5ASCCDKJ7wFhDulP5gZksI7Tp3EhJjSQ8BLxHujM7bmrglu9TtVLKhH6FR8QZgGvAqoTfRN1W0/8uBCwnVOqcBp7v7JxDqwAkX08+AvwFfEXoGdQZmV7Df6wkJ7F+EX9X1gGEltnmIcNF+j5AQekXllxCSzwTgRcKF/hMq4O5FhF46LwF/JSSlF4FjCI3jlXE7cBdwM6FH1XnAee7+bMo2fYBvo7hfJTTQb7r7in7l9yZUAX5G+F4HEZJESRMIdyZT3P3jUtZLNWPumjFNaqaocXMaoVfM1Iq2l+yKqsnmEXobjYg7HqmYqoxEpEpFDc+tgD8SGsKfjDciSZcSgohUtX0IjedzgH7uvjbmeCRNqjISERFAjcoiIhKpsVVGLVu29A4dOsQdhohIjfLBBx8sdfeC0tbV2ITQoUMHpk5VxxIRkS1hZrPKWqcqIxERAZQQREQkooQgIiKAEoKIiERqbKNyeVasWMHixYtZv77kaMuSCxo3bkzbtm2pVUu/Z0SqUs4lhBUrVrBo0SLatGlDw4YNKTHJiNRwRUVFzJs3j6VLl7L99iXnyRGRysi5n1iLFy+mTZs2NGrUSMkgB9WqVYtWrVrxww8/xB2KSPa5w2efwT33wPz5Vb77nLtDWL9+PQ0bNow7DMmgunXrsmHDhrjDEMmOhQthwgR49dXwWhBNx9G0KZx1VpUeKucSAqA7gxyn71dy2urV8NZb4eI/fjx8WsZU2q++qoQgIpJTiorgo4+SdwBvvw3rypk+u1kzOOoo+OUvqzwUJQQRkWybMyeZACZMgKVLy962Th045BA4+mg45hg44ACoXbvs7StBCSGLmjRpsun96tWrqV+/PrWjL/ahhx6ib9++cYVWaa1bt+aZZ57h5z//edyhiFQ/K1fCm28mk8D06eVv36lTMgEcdhikXDsySQkhi3788cdN7zt06MAjjzxCjx49YowoPRs2bKBOncz+U8nGMUSyZsMGmDo1mQDefTeUlaWgAHr0CAmgRw9o2zZ7sabIuW6nNdnGjRu58cYb2WWXXWjZsiV9+/bl+++/B2D69OnUqVOHESNG0KZNG7bbbjseffRR3n33Xfbee2+aNWvGH//4x037GjZsGEceeSTnnnsuTZs2Za+99mLixImb1i9btowzzzyT1q1b065dOwYNGkRRUVGxz55//vk0b96c2267jenTp3P44YfTokULCgoK6NevHytXhnnVTz31VBYvXswxxxxDkyZN+Mtf/sLLL7/MbrvtVuz8Wrduzdtvvw3AlVdeyemnn07v3r3ZZpttGDNmTLnnL1LtzZwJw4bBySeHC/zBB8N114UG4pLJoH79cOG//Xb48MPQk+ipp6B//9iSAQDuXiNfBxxwgJfmiy++KF4Qeu5m57UF2rdv76+++mqxsltvvdV//vOf+7x583zNmjXer18/79+/v7u7T5s2zQG/6KKLfO3atf7CCy94o0aNvFevXr506VKfNWuWN2vWzCdPnuzu7g8++KDXrl3b77//fi8sLPSRI0d68+bN/YcffnB392OPPdYvuOACX7Vqlc+fP987d+7sjz/+eLHPDh8+3Dds2OCrV6/2adOm+Wuvvebr1q3zBQsWeNeuXf2KK67YFHurVq38rbfe2rT873//23fddddi55e6zRVXXOH16tXzl156yTdu3OirV68u9/wr/J5Fsm3ZMvdnnnE/91z3XXap+Pqw337uf/6z+/jx7qtXxxY2MNXLuK7GfmHf2lcuJoQOHTr422+/vWl55syZ3rBhQy8qKtqUEJYuXbppfaNGjXzcuHGblo877jh/8MEH3T1c1Dt06FBs//vss4+PHTvWv/32W2/UqJEXFhZuWvfoo4/6scceu+mzu+++e7nxjx492rt167ZpeWsSwtFHH532+ZekhCBZV1joPnGi+zXXuHft6l6rVvnXgx13dO/Xz33UKPeFC+OOfpPyEoIqbasJd2fOnDkcd9xxxfrZFxUV8d133wFQu3Zttttuu03rGjZsSKtWrYotp7ZTtC1x69m+fXvmz5/PrFmzWLt2LQUFyUmTioqKilXxtGvXrthn58+fz8UXX8ykSZNYuXIlRUVF7LDDDpU659RjVHT+LVu2rNSxRLaYO3z5ZXgW4NVXQ6Nwyv+vzTRuHBqAE43BnTpBDXtmJvcTgnvcEaTFzGjTpg3PPfccBxxwwGbrl5bXLa0Mc+fOLbY8e/ZsdtxxR9q1a0eTJk1Yvnx5mQ95lSy/7LLLaNy4MZ999hnNmzdnzJgxXHPNNWVu37hxY1avXr1pef369SxbtqzMY1R0/iJZsWRJ8aeCS/wfKsYMunQJF/+jjw5tBvXqZS/WDFCjcjUyYMAArrzySubMmQOEcZn++c9/bvX+5syZw7Bhw9iwYQOjRo1i9uzZHHPMMey8885069aNyy+/fNOv/RkzZmxq8C3NypUradKkCU2bNmX27Nncfffdxda3atWKmTNnblru1KkTy5Yt47XXXmP9+vVcf/31mxqts3X+IhVauxZeew2uuAJ++lPYfns4/XR47LHSk0GHDvD738PTT4dnB6ZMgZtuCncGNTwZgBJCtXL55ZfTo0cPjjzySLbZZhsOOeQQ/vvf/271/rp3786HH35IixYtuPnmm3n++efZdtttARg9ejTff/89HTt2pEWLFvTu3ZtFixaVua/Bgwfz9ttvs+222/KrX/2Kk08+udj6gQMHMnDgQJo1a8bQoUNp2bIl9913H3379qVt27a0bt26wmqfqj5/kc24wyefwJAh0LMnNG8eevvccUfo7VNS06Zw0knwwAMwY0boSTR8OJxyCrRokf34M8y8hlSplNSlSxefOnXqZuXTpk2jU6dOMURUvQwbNoxnnnmGCRMmxB1KRuh7lrTNnx+qgcaPD3+W88OH2rWhW7dkO8CBB4YnhXOImX3g7l1KW5dbZyoismoVTJyYbAz+/PPyt99jj2QCOPzwcFeQp7KaEMysD3A9sBOwEOjv7m+ZWSNgCHAaUBf42N27ZzM2EamhNm4M1T2JBDBpEhQWlr39dtuFweESjcE77ZS9WKu5rCUEMzsauB3oDUwBUvssDo9i6QQsAzpnK65cNWDAAAYMGBB3GCKZMWtWcnjo116DEj3YiqlXD372s+RdwP77g6ZfLVU27xAGAYPdfXK0PA/AzPYETgDauvuKaN0HWYxLRKq7FSvgjTeSdwEzZpS//d57JxPAoYeGZwSkQllJCGZWG+gC/MPMvgYaAOOAy4CuwCxgkJmdASwAbnD3Z0vZzznAOQA7lXObV1RUpAnYc1hN7QghW2DDhtClM5EA3nsvVA2VpVWrZALo0QMq+dBkvsrWHUIrQtvAKcChwHrgBeAaYBWwN/AssCNwMPCSmX3h7tNSd+LuwwnVS3Tp0qXUq0Ljxo2ZN28erVq1om7duppdK8e4O9999x0NGjSIOxSpSu7w9dfJaqA33gh3BWVp2BC6dw9J4OijYZ99atxTwdVRthLCmujPv7r7AgAzu5uQEJ4gJIib3H0D8B8zewM4BphW2s7K07ZtW5YuXcqsWbM0726OatCgwWbDckgNtGxZqP9PJIFZs8rf/qc/Td4FHHII6EdBlctKQnD35WY2FyjtV/0nVXmsWrVqsf3227P99ttX5W5FpLIKC0MPoMSwEFOnlj+0TLt2yTuAo44KQ0pLRmWzUfkx4EIze5lwR3AJ8CIwEZgNXGVmtxLaFA4ntC+ISE3lDl98kUwAb74ZJpAvS5MmcMQRybuAPfZQNVCWZTMh3Ai0BL4C1gJjgZvdfb2ZnQg8AlxJaGA+090rmGNORKqdRYuKDw43f37Z29aqBQcdlEwAXbtC3brZi1U2k7WE4O7rgfOiV8l1nxMak0WkJlmzJswIlkgAH39c/va77pqsBjrySGjWLDtxSlo0dIWIpK+oKFz0Ewngrbdg3bqyt2/WLNT/J5LALrtkL1bZYkoIIlK+uXOTCWDChDBnQFnq1Ak9gBIJoEuXMGCc1AhKCCJS3I8/hgbgRBKYVkHv706dkgngsMNgm22yEqZUPSUEkXy3cWPoAppIAJMmhSeFy1JQEJ4GTiQBPROSM5QQRPLRzJnJBPD667B8ednb1q8fxgNKJID99tPgcDlKCUEkH3z/fbjwJ5LAN9+Uv/1++yUTwKGHhqEiJOcpIYjkovXrYfLkZAKYMiX0ECrLjjsmE0CPHmGwOMk7SggiuWLJEhgzJvlU8MqVZW/bqFGYHSyRBPbaS08FixKCSE745pvQ3XPx4tLXm4UuoIkEcPDBoW1AJIUSgkhNV1gIffpsngzat09OE3nkkWHqSJFyKCGI1HRXXhm6jUIYC+jOO+G442C33VQNJFtECUGkJnvxRbjnnuTybbfBxRfHF4/UaBUmBDOrBRwA7Ac0A74HPgY+cPdyui2ISEbNnQv9+iWXjz8eLr00vnikxiszIZjZtoQ5Cc4CCoHpwEpgG+A6oK6ZjQCGuPv3WYhVRBI2bIDTTw+zjgG0aQOPP64qIqmU8u4Q3gNGAT9z95klV5rZzkBf4F2gU2bCE5FSDRoURhqF8NTw6NHQsmW8MUmNV15C2N/d15S10t3/B9xkZndVfVgiUqbXXoObb04u33BDeJpYpJLKHJDE3deYWT0z27209Wa2u5nVKy9piEgVW7QIfvOb5FzERx4JV18db0ySMyoaoeoyoKwuCxcCf6racESkTEVFcOaZsHBhWC4ogFGjNN+AVJmKEsKvgTvKWHcnoQ1BRLLhzjth/Pjk8t/+BjvsEF88knMqSght3X12aSvcfQ7QrupDEpHNTJoEAwcml6+4Anr2jC8eyUkVJYS1ZtamtBVm1hZYW/UhiUgxy5bBr38dJrKBMA7RjTfGG5PkpIoSwsvA4DLWXR+tF5FMcYezzoLZ0Y16s2ahi2nduvHGJTmpoieVBwKTzWwK8CywANgB6BX9eXBmwxPJc/ffD+PGJZdHjAiD1olkQLl3CO4+D/gpMBHoA9wU/TkROCBanzYz62Nm08xslZl9Y2aHllh/vZm5mfXYorMQyUUffgh/SunId8EF0KtXfPFIzqtwLCN3XwL8ubIHMrOjgduB3sAUwh1G6vpdgVMIdyEi+W3lSjjttDC0NUDnzqGXkUgGZXOm7EHAYHef7O5F7j6vxB3GUOAKwrhJIvnLHf7wB/j667DcuDH8/e/QoEG8cUnOy0pCMLPaQBegwMy+NrO5ZjbUzBpG608FCt39XxXs5xwzm2pmU5csWZKFyEVi8Pjj8OSTyeVhw2CPPWILR/JHtu4QWgF1CVVChwKdgf2Ba8ysCXALcElFO3H34e7exd27FBQUZDJekXhMmxbaChL69w9DVYhkQbYmyEmMd/RXd18AYGZ3A9cADYC/RYPlieSvNWtCu8Hq1WG5Y0cYOjTemCSvVOoOwcx+lc527r4cmAt4KauPAi4ys4VmtpDw9PNYM7uiMrGJ1DiXXgqffRbeN2gAY8eG9gORLKlsldGWdHt4DLjQzLY3s+aEKqIXCQlhb0I1UmdgPnAucH8lYxOpOcaOhYceSi7fey/ss0988UheqlSVkbvvtgWb3wi0BL4iDHkxFrjZ3YsNf2FmG4Hl7v5jZWITqTFmzoTf/z65fOqpcM458cUjeStbbQi4+3rgvOhV3nYdshKQSHVQWAi9e8OKFWF5553h4Yc1FabEIq2EYGavUnr9/zpC28Bz7v5qVQYmkheuugqmTg3v69YNzxtsu228MUneSrcNYSqhnv9jQr3/x8BPgOmEpPB3M6uw26iIpHjxRbj77uTybbfBgQfGF4/kvXSrjA4Derr7p4kCMxsJDHf3Q8xsLPAEcG8GYhTJPXPnQr9+yeXjjw+9jERilO4dwk+AL0uUzYjKcfdJQOsqjEskd23YAKefHuY5AGjTJjydrHYDiVm6CWESMNzM2lrQFngAeBfAzPYCFmUoRpHcMngwvPVWeF+rFjz1FLRsGW9MIqSfEM4EmgGzgA3At0ALIHHPa8AZVR2cSM55/XW46abk8g03QPfusYUjkiqtNoRoCOyTzKw+oWpoobuvS1n/eYbiE8kdixdD375hNFOAI4+Eq6+ONyaRFGk/qWxmuwAXARe5+zoz283MOmUuNJEcUlQEZ5wBCxeG5YICGDUKateONy6RFGklBDM7EXif0IiceKRyO+C+DMUlklvuvBPGj08uP/EE7LBD2duLxCDdO4RbgGPcvT+wMSr7kDD2kIiU5913YeDA5PLll8Oxx8YXj0gZ0k0IO7j7B9H7xBPLRYTGZBEpy/Ll0KcPbIx+R3XrVrxRWaQaSTchfGhmfUqU9SI8wSwipXGHs86C2bPDcrNmMGZMGKJCpBpK90nli4HxZvZboJGZvUCY8axnxiITqenuvx+efz65PGIEtG8fXzwiFUi32+lnZrYn8CvgbWAO0M/dv89kcCI11ocfwp/+lFw+/3zo1Su+eETSkPbw1+6+kjBekYiUZ+XKMKR1YWFY7twZhgyJNyaRNJSZEMoZ8roYdz+mSiMSqcnc4bzzYMaMsNy4cRjSukGDeOMSSUN5dwjPpLxvD/wOeJIwfMVOQF/g0cyFJlIDjRwZHjhLGDYM9tgjvnhEtkCZCcHdN03wambvEIa//jil7AngwcyGJ1KDTJsW2goS+veH3/wmtnBEtlS63U73BqaVKPuKaPhrkby3Zg2cdhqsXh2WO3aEoUPjjUlkC6WbEN4BHjGzdtHw1zsR7g7ezVxoIjXIpZfCZ5+F9/Xrh3aDxo3jjUlkC6WbEPoBTQnDXm8AZhKGw+6fkahEapKnn4aHHkou33cf7LtvfPGIbKUqGf5aJG/NnAlnn51cPvVUOOec+OIRqYQy7xDMrEXJMndf5+6zUpNBaduVs88+ZjbNzFaZ2TdmdqiZdTOzV81smZktMbOnzUzDQEr1V1gYnjdYsSIs77wzPPywpsKUGqu8KqN3zOxuM9u/tJVm1tnM7iY8uVwhMzsauB34LbAN0J1Q9dQcGA50IHRvXQk8lu4JiMTmqqtgajScV506YZyibbeNNyaRSiivyuinwPnAWDMrAGYQLtbbALsR5lB+KNouHYOAwe4+OVqeV+JPAMxsKPCfNPcpEo+XXoK7704u33YbHHRQfPGIVIHynkNYAwwBhpjZ7sA+hIbk5cAn7v5Nugcxs9pAF+AfZvY10AAYB1wWHSdVd0BTckr1NXcu9OuXXD7++NDLSKSGS7dReQbhDmFrtQLqAqcAhwLrgReAa4BNM4eY2b7AdcCJpe3EzM4BzgHYaaedKhGOyFbasAFOPx2++y4st2kDjz8OtdKejVak2srWv+LEXcBf3X2Buy8F7gaOS2xgZrsB/wYudve3StuJuw939y7u3qWgoCDjQYtsZvBgeCv651mrFjz1FLRsGW9MIlUkKwnB3ZcDcyljsDwzaw9MAG50979lIyaRLfb668VnO7v+eujePb54RKpYNu9zHwMuNLPtzaw5cAnwopm1AV4H7nf3YVmMRyR9ixdD375hNFOAI44oPk+ySA7Y4oSwJc8dlHAj8D5hDKRpwIfAzcDZwC7A9Wb2Y+K1lccQqXpFRXDmmbBwYVguKAgjmtauHW9cIlUsrYRgZk3N7FEzW0WYLQ0z+z8zuzbdA7n7enc/z92buXtrd7/I3de6+yB3N3dvkvrayvMRqXpDhsArrySXn3gCdtwxvnhEMiTdO4T7oz/3BqJpoJgCnFHlEYlUJ+++C1dfnVy+/HI49tj44hHJoHSn0OwJtHX3QjNzAHdfbGatMheaSMyWL4c+fWDjxrDcrVvxRmWRHJPuHcJKwhATm5hZW8LTyiK5xx3OOgtmzw7LzZrB6NFQt268cYlkULoJ4XHCEBYHAxaNb/Qo8HCmAhOJ1QMPwPPPJ5dHjIAOHWILRyQb0q0yupnwdPGTQBPgOcI4RkMyFJdIfD76CP74x+Ty+edDr17xxSOSJRUmhGgcot7APe5+W+ZDEonRypVhKszCqO9E586hl5FIHqiwysjdNwIPakIcyXnucN55MCMatqtx4zAVZoMG8cYlkiXptiH828x6ZjQSkbiNHBkeOEt48EHYY4/44hHJsnTbEDYC48zsP4QH0zaNSeTumi9Qar5p00JbQUL//nCGHrOR/JJuQpgN3JvJQERis2ZNmApz9eqw3LEjDB0ab0wiMUh3PoSrMh2ISGwuvRQ+/TS8r18/tBs0bhxvTCIxSCshmNkhZa1z90lVF45Ilj39NDz0UHL53nth333ji0ckRulWGT1bYrkFYMBSQKN8Sc00cyacfXZy+ZRT4Nxz44tHJGbpVhntkLpsZnWBQcD8TAQlknGFhWGcohUrwvLOO8PDD4NZvHGJxGirJshx9/WEuY81Q4jUTFdfDe+/H97XqQNjxoTxikTyWGVmTDu0yqIQyaaXXoK77kou33YbHHRQfPGIVBPpNirPoPh8yI0I7QgXZyIokYyZOxf69UsuH3dc6GUkImk3Kl9QYnkVMM3dv6vieEQyZ8OGMC/yd9E/2x13DE8n18rm1OIi1Ve6CWF3d9/sSR0zO8/dH6jimEQy48YbYeLE8L5WLXjqKWjZMt6YRKqRdH8a3VJGuaaPkprh9ddDQki4/no47LD44hGphsq9Q0h5IK1WYnKclNW7AD9mKjCRKrN4cagq8qgZ7IgjYKA6yImUVFGVUeKBtIaESXESHFgI/HGzT4hUJ0VFcOaZsHBhWC4oCCOa1q4db1wi1VC5CSHxQJqZjXX307ITkkgVGjIEXnklufzEE6ExWUQ2k1YbQlUlAzPrY2bTzGyVmX1jZodG5UeZ2XQzW21mb5hZ+6o4nuS5yZOLVw1dfjkce2x88YhUc2klBDNrbGa3mNk7ZvalmX2VeKV7IDM7Grgd+C2wDdAdmGlmLQnVUdcSnm2YCvx9S09EpJjly8PQFBs2hOVu3eAm9YEQKU+6vYzuBw4H/kIYzO5aYDnw8BYcaxAw2N0nu3uRu89z93lAL+Bzd3/a3dcCNwD7mVnHLdi3SCI+DvMAABUOSURBVNKMGXDyyTBrVlhu1gxGj4a6deONS6SaSzch/AI4yd3/DmyM/jwVSKsqycxqA12AAjP72szmmtlQM2sI/AT4OLGtu68CvonKRdK3aFGY9WyvveCNN5LlI0ZAhw6xhSVSU6T7YFodIPFU8o9m1hSYC+yZ5udbAXWBUwhjIK0HXgCuAZoAS0ps/wOhWqkYMzsHOAdgp512SvPQkvNWrgxjEw0ZAqtWFV83cCD06hVPXCI1TLp3CJ+QHMxuEmE6zXsJv+TTsSb686/uvsDdlwJ3A8cRnmVoWmL7psDKkjtx9+Hu3sXduxQUFKR5aMlZ69fD/ffDbrvBoEHFk0H37qFRWe0GImlLNyEMABZE7y8G6gHtgf7pfNjdlxPuKLyU1Z8D+yUWzKwxsGtULrI5dxg7NlQNXXBBePAsYe+94cUX4c03oWvX2EIUqYkqTAhR/f/JwLcA0S/837j7ie7+cbkfLu4x4EIz297MmgOXAC8CzwN7m9nJZtaAMM/CJ+4+fQvPRfLBG2+EC33v3vD118nytm3hscfgo4/g+OM10Y3IVqgwIbj7RuDPQGElj3Uj8D7wFTAN+BC42d2XEBLOzYSeS12BPpU8luSaTz4JQ1UfeWRyYhsIPYjuuAO++gr699cTyCKVkG6j8lOE5wce3doDRbOsnRe9Sq6bAKibqWxu9my47rrwhLGn1DjWrw8XXQRXXQXNm8cXn0gOSTchdALOMbPLgTmktAW4+zGZCEzy3LJlcOut8Ne/wrp1yXKzMMHN4MHQrl188YnkoHQTwtjoJZJZa9aEJHDrrfD998XXHX98KN9nn3hiE8lxaSUEd38o04FIntu4MVQLXXddmOYy1UEHhXYCzV8gklFpzx1oZmeY2b/M7P1o+Wdm9qvMhSZ5wT10E91vP/jd74ong913h6efDs8TKBmIZFy6g9tdR+gmOhbYPSpeBGiWEdl6kyfD4YfDL38Jn6c8dtKqFTzwQCg75RR1IRXJknTbEM4CDnL3RWZ2T1T2DWHWNJEt8+WXYUiJZ58tXt6kCVx2Gfzxj+G9iGRVugmhHpBo4Uv0MGoMrK7yiCR3LVgQegc9/HBoM0ioUwcGDIBrr4Xtt48vPpE8l24bwnjg9uip5YTrgH9VfUiSc1asCI3Fu+0Gw4YVTwa9e8O0aaFnkZKBSKzSvUO4BHiSMAppAzNbDkwGfp2pwCQHFBbCQw/BjTfCkhID2h5xBNx+Oxx4YDyxichm0u12uhw4zsx2IgxqN8fdv81kYFKDFRWFwecGDoSZM4uv23ffkAh69lRjsUg1k+4dAmbWBDiQMGPa9mb2nbtvNkS15LnXXoMrroAPPihevtNOYSjqvn2hVtq9nUUki9JKCGZ2KDCOMIT1LGAnYLiZneTub2UwPqkpPvoIrrwSXnmleHmLFuFO4bzzoEGDeGITkbSke4fwIHCJu/8tUWBmvwGGoaku89u334beQU8+WXzwuQYN4JJLwt1Cs2axhSci6Uv33r0dYcTTVKOBtlUbjtQY330XnhfYc08YNSqZDGrVgrPOChPd33qrkoFIDZJuQhgNnF2i7Cw2TxKS61avDhf6XXaBe+4JPYkSTjghzFvwyCNhwhoRqVHSrTLaHTgrGv56LuHOoB0w0czGJzbSUNg5bMMGePxxuP56mD+/+LqDDw49hw49tNSPikjNoOGvpXzu8I9/hIlopk0rvm7PPcPdwkknqQupSA7Q8NdStkmT4PLL4Z13ipe3bg2DBoXRSeuk3XNZRKq5LXkO4UBgf6DYqGPufndVByUxmz493BGMG1e8fJttQq+hSy6Bxo3jiU1EMibd5xDuIsypPAlYk7LKS/+E1Ejz58MNN8CIEeFp44S6dcNzBAMHQkFBbOGJSGale4fwW2A/d5+TyWAkJj/8EGYku+eeMIVlqtNPD2MR7aKRzkVyXboJYR6gYSpyzbp18OCDYUiJ774rvq5Hj9Bz6Kc/jSc2Ecm6dBPC74GHzWwksDh1hbtPqfKoJLOKimD0aLjmmvCkcar99w+J4OijYwlNROKTbkLoBPwiepVsQ0hrEHszexPoBmyIiua5+57RutOAQYTnG+YAV7v7uNL2I5U0fnxoGP7oo+LlHTrAzTdDnz4afE4kT6WbEO4A+rj7i5U83gXu/khqgZm1AUYBJwIvA8cBT5tZB3dfXMo+ZGv8978hEUyYULx8u+3CWEQDBkD9+vHEJiLVQroJYR3wSoVbbZ22wPfu/u9o+SUzWwXsSonqKdkKhYVhzKH77y9e3rBhKL/sMth223hiE5FqJd26geuAO8yssiOV3WpmS83sHTM7PCqbCkwzsxPMrLaZnURIQJ+U/LCZnWNmU81s6pKSM3DJ5hYtgqOOKp4MateGc86Br78OjclKBiISMfeKHyUws/VAbUKbQWJCXAPc3euldSCzrsAXQCHQBxgKdHb3b8zsLOA+oEG0/lR3f6m8/XXp0sWnTp2azqHz0wcfhCEl5s5Nlh1/PAwZAh07xheXiMTKzD5w9y6lrUu3ymjvygbh7u+lLI40s18TpuWcRmijOBz4L3AA8A8z+4W7f7T5nqRCTz0VhqBeuzYsm8Ftt4XqIY05JCJlSHcsoy8T782shbsvq4JjO+EuozMw0d0TP/ffN7P3gB6AEsKW2LgRrr46PGSWsO22oYvpL34RX1wiUiOk1YZgZk3N7NGosXdOVPZ/ZnZtmp9vZmY9zayBmdUxs75Ad0JD9fvAoWbWOdp2f+BQSmlDkHIsXw7/93/Fk8Gee8KUKUoGIpKWdBuVE62SexPq+AGmAGek+fm6wE3AEmApcCFwkrt/6e7/AW4AnjGzlcCzwC3uPr6snUkJ06ZB167w8svJsuOPh/fegz32iC8uEalR0m1D6Am0dfdCM3MAd19sZq3S+bC7LwEOLGf9UEIjs2ypF18M4w2tTBlZ5OqrYfDg0KNIRCRN6d4hrASapxaYWVtgUZVHJOlxh1tuCdNWJpJBw4YwZkx44ljJQES2ULp3CI8DY83sSsCiev7bgYczFZiUY9Uq+O1v4emnk2Xt24f5Czp3ji8uEanR0k0INwPrgScJE+Q8BzwEDMlQXFKWb78Nzxd8/HGy7LDDQnLQXAUiUgnlVhlFzwrg7kXufpu77+Ludd1952hZE+Rk05tvwoEHFk8G558Pr76qZCAilVZRG4LmUq4O3MPwEz16wNKloaxuXRg+HIYODe9FRCqpoiojPdYat3Xr4IIL4JGUQWJbtYJnn4Wf/Sy+uEQk51SUEGqb2RGUkxjc/fWqDUk2WbgQTj4ZJk1Klh1wADz/PLRrF19cIpKTKkoI9YERlJ0QHNBku5nw/vvwq1/BvHnJst/8JlQTNWwYX1wikrMqSgir3F0X/GwbNQrOPjtUF0GYweyOO8L8BRqcTkQyJN1up5INGzbAlVfCXXcly5o1Cw+b9ewZX1wikhfUqFxdLF8e5jMenzKEU6dO8MILsPvu8cUlInmj3G6n7r5NtgLJa198AQcdVDwZnHACTJ6sZCAiWZPuWEaSKS+8EEYq/frrZNm114aeRE2bxheXiOQdtSHEpagoDEJ33XXJskaNYORIOOWU+OISkbylhBCHH3+E/v3Dw2UJHTqEu4V9940rKhHJc0oI2fa//8GJJ8KnnybLjjgCxo6Fli3ji0tE8p7aELLp9dehS5fiyeDCC+GVV5QMRCR2SgjZ4A5//SsccwwsWxbK6tWDESPgL3/R4HQiUi2oyijT1q2DP/wBHnssWda6NTz3HBx8cHxxiYiUoISQSQsWQK9e4XmChAMPDF1K27SJLy4RkVKoyihTpkwJ7QWpyeDMM2HiRCUDEamWlBAyYeRI6N4d5s8Py7VqwT33wOOPQ4MGsYYmIlIWVRlVpQ0b4LLL4N57k2XNm8Pf/w5HHx1fXCIiacjaHYKZvWlma83sx+j1Zcq6Rmb2gJktNbMfzGxituKqMt99B8ceWzwZ/OQnYV4DJQMRqQGyfYdwgbs/Ukr58CiWTsAyoHNWo6qszz4LD5vNnJksO+kkeOIJ2EbjA4pIzRB7G4KZ7QmcAJzj7kvcfaO7fxB3XGl7/nno1q14Mrj++jAshZKBiNQg2U4It0bVQu+Y2eFRWVdgFjAoWvepmZ1c2ofN7Bwzm2pmU5csWZKtmEtXVASDBoVupatWhbLGjcPzBTfcEBqSRURqkGxeta4gzL/chlBF9E8z2xVoC+wN/ADsCFwAjDSzTiV34O7D3b2Lu3cpKCjIXuQlrVwJJ58cLvwJO+8M774b5kEWEamBspYQ3P09d1/p7uvcfSTwDnAcsAZYD9zk7oXu/h/gDeCYbMW2Rb75JjxhPG5csuyoo0Lj8T77xBeXiEglxVmv4YQpOj+JMYYtM2FCeNL488+TZZdcAi+/DNttF19cIiJVICsJwcyamVlPM2tgZnXMrC/QHXgFmAjMBq6K1v0MODxaVz24hwfLevYMcx9DGJzuscdCeR09ziEiNV+2rmR1gZuAjsBGYDpwkrt/CWBmJwKPAFcSGpjPdPfpWYqtfGvXwrnnhi6kCTvsEHoXde0aX1wiIlUsKwnB3ZcAB5az/nOg+g39OW9e6EU0ZUqyrGvX0JNoxx3ji0tEJAPUN7Is774bBqdLTQb9+8ObbyoZiEhOUkIozWOPweGHw8KFYbl2bbjvPnj0UQ1OJyI5S62hqdavhz//OcxiltCiBTz9NBx5ZHxxiYhkgRJCwtKlcNpp8MYbybJ99gnPG+yyS3xxiYhkiaqMAD75JDxfkJoMTj4ZJk1SMhCRvKGE8Oyz4cnjb79Nlg0eDGPHQpMmsYUlIpJt+VtlVFQUxiK68cZkWZMmMGpUGMpaRCTP5GdCWLECzjgD/vGPZNmuu8ILL4RJbURE8lD+JYSvvw53AF98kSw7+mgYMyb0KBIRyVP51YaQaDxOTQZ/+hP8619KBiKS9/IrIXTsCHvtFd7Xrx/GJxoyRIPTiYiQbwmhXr3Qq6hbN3jrrdCOICIiQD62IbRuHZ4vMIs7EhGRaiW/7hASlAxERDaTnwlBREQ2o4QgIiKAEoKIiESUEEREBFBCEBGRiBKCiIgAYO4edwxbxcyWALMqsYuWwNIqCqemyLdzzrfzBZ1zvqjMObd394LSVtTYhFBZZjbV3bvEHUc25ds559v5gs45X2TqnFVlJCIigBKCiIhE8jkhDI87gBjk2znn2/mCzjlfZOSc87YNQUREisvnOwQREUmhhCAiIoASgoiIRPI6IZjZKDNbYGYrzOwrMzs77pgyyczqm9kIM5tlZivN7EMz+0XccWWamV1gZlPNbJ2ZPR53PJlgZi3M7HkzWxV9v6fHHVOm5cP3miob/3/zb8a04m4FznL3dWbWEXjTzD509w/iDixD6gBzgMOA2cBxwFgz28fdv40zsAybD9wE9AQaxhxLptwPFAKtgM7AS2b2sbt/Hm9YGZUP32uqjP//zes7BHf/3N3XJRaj164xhpRR7r7K3W9w92/dvcjdXwT+BxwQd2yZ5O7Pufs44Lu4Y8kEM2sMnAxc6+4/uvvbwD+AnJ40PNe/15Ky8f83rxMCgJk9YGargenAAuBfMYeUNWbWCtgDyOVfkflgD2Cju3+VUvYx8JOY4pEsyMT/37xPCO5+HrANcCjwHLCu/E/kBjOrCzwJjHT36XHHI5XSBPihRNkPhH/XkoMy9f83ZxOCmb1pZl7G6+3Ubd19Y3Sb3Rb4QzwRV16652xmtYC/EeqcL4gt4CqwJd9zDvsRaFqirCmwMoZYJMMy+f83ZxuV3f3wrfhYHWpwG0I652xmBowgND4e5+7rMx1XJm3l95xrvgLqmNnu7j4jKtsPVQXmnEz//83ZO4SKmNn2ZtbHzJqYWW0z6wn8Gng97tgy7EGgE/BLd18TdzDZYGZ1zKwBUBuobWYNzCxnfgy5+ypCdedgM2tsZj8DTiT8isxZuf69liGz/3/dPS9fQAHwH+B7YAXwKfD7uOPK8Dm3J/SkWkuoZki8+sYdW4bP+waSvcgSrxvijquKz7EFMA5YReiSeHrcMel7rfLzzfj/Xw1uJyIiQB5XGYmISHFKCCIiAighiIhIRAlBREQAJQQREYkoIYiICKCEIDnMzK42s0eyeLx3zGz/6P0NZjYqW8cuJ6a7zWxA3HFIzZDrT/VJDjOzH1MWGxEGJtwYLZ/r7rdkMZZfAivd/cNsHTNNdwJTzOxRdy+MOxip3nSHIDWWuzdJvAhP5/4ypezJLIczgBiGiqhoqAZ3X0AY2v2E7EQkNZkSguSs1GobM+sQjYD6WzObY2bLzWyAmR1oZp+Y2fdmNrTE539nZtOibV8xs/ZlHKcecCRhKJRU9czsiWi6w8/NrEvKZzpFI7V+H607IWXdm6nTuZpZ/xKj1bqZnW9mM4AZFtxjZovN7IfofPZOieNN4Pgt/guUvKOEIPmmK7A70Bu4FxgI9CBMJnOamR0GYGYnAVcDvQjjXr0FjC5jn7sDRe4+t0T5CcAYoBlhBrOh0b7rAv8ExgPbAxcCT5rZnltwHidF57IXcAzQnTBZSrPo3FJnEZtGGP1UpFxKCJJvbnT3te4+njAQ3Gh3X+zu8wgX/f2j7c4FbnX3ae6+AbgF6FzGXUIzSp974G13/5e7byRUJyUuyt0Ik9rc5u6F7v468CJhtN103eruyzyMeLmeMBlOR8CimBekbLsyilGkXEoIkm8WpbxfU8pyk+h9e+C+qErne2AZYECbUva5nNJnJ1uY8n41kBieeUdgjrsXpayfVca+yzIn8SZKKEOB+4FFZjbczFInzNmGMKqvSLmUEERKN4fQU6lZyquhu08qZdsZhLlL0r2gzwfaRTNfJewEzIveryL0mkpoXco+ig1T7O5/cfcDCFVfewCXpazuRJhjWaRcSggipRsGXGVmPwEws23N7NTSNvQwa9UE4LA09/0e4aJ/uZnVNbPDgV8S2hsAPgJ6mVkjM9sNOKu8nUUN412jtolVhPHyN6Zschjw7zRjkzymhCBSCnd/HrgdGGNmK4DPgF+U85GHgDPS3HchocH5F8BS4AHgTE9Oln4PYb7cRcBIwmTq5WkKPEyouppFaFAeAmBmOxAanselE5vkN02QI1JFoq6hF1anh9PM7C7gG3d/IO5YpPpTQhAREUBVRiIiElFCEBERQAlBREQiSggiIgIoIYiISEQJQUREACUEERGJ/D+T0YIE2lPqNgAAAABJRU5ErkJggg==\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"import numpy as np\n", | |
"t = -2 #input value of t here to get ambient temperature at specific time point\n", | |
"\n", | |
"def ambient_temp(t):\n", | |
" if t >= -3 and t < -2:\n", | |
" # temp_at_t = 55\n", | |
" # linear interpolation (1-u)*x0 + u*x1, u is in [0,1]\n", | |
" temp_at_t = (-2-t)*55 + (3+t)*58\n", | |
" elif t >= -2 and t < -1:\n", | |
" #temp_at_t = 58\n", | |
" temp_at_t = (-1-t)*58 + (2+t)*60\n", | |
" elif t >= -1 and t < 0:\n", | |
" #temp_at_t = 60\n", | |
" temp_at_t = (-t)*60 + (1+t)*65\n", | |
" elif t >= 0 and t < 1:\n", | |
" #temp_at_t = 65\n", | |
" temp_at_t = (1-t)*65 + t*66\n", | |
" elif t >= 1 and t <= 2:\n", | |
" #temp_at_t = 66\n", | |
" temp_at_t = (2-t)*66 + (t-1)*67\n", | |
" else:\n", | |
" print(\"No data at this time point\")\n", | |
" temp_at_t = 0\n", | |
" return(temp_at_t)\n", | |
"\n", | |
"print(\"Based on your current t value of \", t, \", the ambient temperature is \",ambient_temp(t), \"deg. C\") \n", | |
" \n", | |
"xaxis = np.arange(-3, 3, 1)\n", | |
"yaxis = np.array([55, 58, 60, 65, 66,67])\n", | |
"\n", | |
"#Plot\n", | |
"plt.plot(xaxis, yaxis, color='red', linestyle='-', label='Temperature')\n", | |
"plt.xlabel('Time (hours)')\n", | |
"plt.ylabel('Temperature (deg. C)')\n", | |
"plt.title('Temperature Today')\n", | |
"plt.legend(loc='best'); \n", | |
"\n", | |
"'''The linear interpolation is flawed in this case. The temperature change between measured time \n", | |
"points is not linear in reality. A better way to do this would be to fit a polynomial curve to \n", | |
"given time points (or some higher-order curve). Numpy.polyfit is a python function that could accomplish this'''\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 227, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAg4AAAEdCAYAAABkEe+RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3gU1frA8e+bQgKEUEKkFykioCJVQhMriuiVoiDdH4pcRVHgYhdQ9IqA2EXgKs2GBa8dOyKgFKWKBbxgCCK9BEiAcH5/nNlkstkkCymTbN7P8+yTzJnZmXdmZ2fePXPmjBhjUEoppZQKRpjXASillFKq+NDEQSmllFJB08RBKaWUUkHTxEEppZRSQdPEQSmllFJB08RBKaWUUkEr8YmDiAwWkRNex6G8JSLfiMjMfJjPOBHZlB8xueaZL7HlYflGRPp7tfzCICJbROSBfJxfXWe7dciveariR0R2iMhor+PIb0ElDiISJyJPiMivIpIiIjtF5FsRGSgiEQUdZAF7E6jhdRCuA43vdUJEEkXkBREpd4rzekBEthRQqIVKRKJEZLeIHBWRyl7H4yMiHZzPqa7fqMlA28KPyDtO8m1EZJWIhPmN8zTpOQWtgamFtTBnu5hcXnULKx4viMhNIpJSBOI4O4fPYHiQ85gnIp8GGHUu8EL+Rhxw+Zc68VYt6GUB5HrSF5GawBLgBPAQ8BNwHGgHjAbWAqsLMMYCISICRBhjjgJHvY7H5R/Acuxncw7wHyAKGOJlUB7qCWwF/gYGY0/MRZYxJhlI9joODxigKTAQmOVtKMETkVLGmGPGmF2FvOgeQCnX8DZgFPaHjE9hx5QvfNu0kJcZaYw5nsfZXAGs8Ss7mJcZerBfFQ5jTI4v4ANgB1A+wLhIoKzr/8eBJOAY8DPQ1296A9yO/XIcBv4EegHlgVeBQ8AfQE/Xe+o67xsAfIk9yf8P6Oc370eBjcARIBGY5o4Ze9I5AVyETX6OAd185a7pYoFXnHVOdeb1pN86B7OetwJznXVKBMbksp1969nBr3wKsMavrCXwGfYEtQt4F6jjWk/j9xoH3ARsC7C8ea6yG7EnaHGGq2BPAruc9VgCdPKLpQHwDrAf2OfEdW6A7d4e+NH5fFYALXPb95z3LwLuAK4Hfg0w/htgJvCg85ntdWIu65qmBfAJsNPZZiuAKwLNx7Ud9gNl/KYZi933zgywjb9xphkHbPJ736XAYmfdDzjrVP90YsthOwW7/+f4OWC/H2uBFOfvRc769c9h2b55+74XZbOLPdC6AA8AW1zDs4AvsMeKbc52mYn97g3DJpL7gOlAKb953Q784sT/O3A/9geCb/wWYAL2V+AeYIWr/AHXdBHYH0qbsceBJOBZ1/gR2B9Mydj97g2gWm7f5xy24QlgcDbjBrg+k/8BTwClXeO/d9ZnIrDb2TZjgXDgEWff+hsY6zffHc50s7Df713O9OKappSzb23FHnvXAze6xkc76/lPYL4zn3nOuEnOZ3EEe6x/Fohxxl1B1u/QNNf6POcX6wTgF9fwG8CH2GRrK3DS+cwEuAv4zdlevwJjgPActv3ZzvJb5TBNFPCMsx+kAtuB2c64xwOsSx/XNh7tt80fwu7PB53hoUBp4EXscScRuNlv+b4f6YedZc8DzvCL3/369BT2n4uAZdh9+SD2/HhRjvtrLjtzJSAN1xcqh2knYb+I1wFnAfc5H+YlrmmMs6EGYU84Lzg71SfYg08DZ+c6DMT5fQG3A/2ARs5OdNL9QWMPPh2d6S/B7rCz/Q5uJ7EHy4uBekA8WROHZ7BZ5wVAbWzNys2nsZ5/AzcD9bEnPpPTh0GAA40T4wbgRVdZE+cDHu/sMOcCb2G/KNHYHfBxZ+er6rxiyDjZNXLmMwR7QNnumvdc4A3n/9LYpOgdoJXz2dyP/dI0dqap4nyeLzpxNHI+vz1AvN92/9b5fM7GJhebcR3Qc/hCpwKVsV/cvUBnv2m+wX7ZpjrTX+EMj3dN0xm7zzVxPrMJ2KTvLL/5zHSt+z5gkGt8GPbkcj/2gHyNsz1bO9u4kjPdOFyJAzZpSAOeApo5MQ4Bzj6d2HLYVsHu/9l+DkB17HfvFSeey7AHnGATh1jsfj8+u9gDrQuBE4cDwGygsbOtU4CPgTlObN2wJ7J/ut43DnsS6Y7d37tiT1iPuKbZgj04jnO2dxNXuTtxmI39fgzAfofbAne5xo9wPtszgQRgKbAop+9zLp9fwMQBmyjtBvpijwcXYb+XM1zTfO9sr0eAhtiTuHG21wRnPYfidwzCfncPYpPus7AJ81HgVtc0b2ATzUucde2LTQ76OeN9icMuJ9b6QANn3EPYRLUucDmwCXjJGVcKGOl8rr7jVKxrfYJJHA5hk5XzsN8twR77/nD2mTOBq7HnjvtzOc7kljjchz3pdsKeF9oAtzvjYoC3ga9c6xLt2sb+icNebILbwPnM0pzPylc2ztkfGrjeNxJ73jrT2abLgYXOuHDs+chgj8NVgYrB7D/Y4+ohZ7s1cPaDnkBCjvtrLjtzGyeYHrlMVwZ7gL/Vr3wB8JVr2ABPuYbjnTJ3Jl/RKevm9wV8xG/eS3H9Wg4QU3cnpjDXwc0AHQMd9FzD/wVm5cN6PuM3zS/Av3OI17eeR7CJQYoz/BVOlu46qL7h994o533XOsOZDsSu6f7nix1bwzMee+DwHTy3AUNd22Ubfid3J56nnP/HAd/7jRfsyehOv+3ewjVNW1xJTA7bZCqwwDX8AvCa3zTfAGv9yqYBy3KZ9xpcBxOynuCeAb5zDXfBXqKr5gx3cNahrt98x5E5cVgMfJhTLKcaW5DzyG7/z/ZzwB6ct5L5F3o3gkwcnP9vwSYfNbPZrlnWxX9/xe7jO3HVJgAfYQ+AUX7f1bdd380jZK2tGQjsdw1vAb4MsA5bcBIH7AHUAL1OYXs3d95Tw+/7fNqJA/a79FeA8suxSWAZZ/h74Ae/aTbj1Ka4yn4FJriGdwCf+03zpG//JeOEeqbfNI/hfO/JSByeD2IdbwAOuYZvAlICTBds4rCbzL+cy2P3+c5+7x0K7MghLt96+o697ldzZ5qXsD9wJZt5zMP1K99vG/snDm+4hiOwx/q3XGXh2O/QTTnEnODE7PuBfakzXPVU9h+gmvO+tsHu68aYXBtHivPX5DJdA2wG+a1f+SLsdU+39GtIxl7/ScP+qvGV7cP+4jrD733L/IaXYH952EBFejgNNreLSDL2xFgKm325rchlXV4AeonIehF5WkSudDX4OpX19G/3kYT9hZ6bG4HzsVl0F2w2+74rhtZAdxFJ9r2wv/Cjsb82cvI1NmsFm3kuxJ7YLhaRRthGol+5llMV2O+3rI6u5bQGWvqNP4Q9aLpjMWS+dpjk/M12e4hINPagP9tVPAvoISJxfpPnuK1FJN5pZPqLiOx34mwK1Mlu+dgDRXsR8e1jNwMfGWP+yuE9gfguKwV0mrEFmk8w+39un0MTYLkxxn2X0XenEge2CvYP7MklLzaazNfJd2AvVaX6lfmOE02xNUXv+O2PLwHlRSTe9b7luSy7hfM3p8+ts4gsdBowHyJjO53S55aLmtjP7wW/dVqAPTbXd03r/x3Ygeu46ioL5rhaz/n+tXbK1vktfyRZjzVZtqmI9BaR70TkL+d9LwMxIlIpuxU+ReuMbaPmcx52n//IL96ngSpBNDLviz32ul8/O+NmYn9I/+Z8X7uLSORpxu0+B57AHr/d58A0bFKU/lk5jR8/d+1vXzijctrfct1/nOPZPOAbEflIRMaISIPcViC3xpG/YzOTps7CcuOfYEiAskANWPzLDLnf8SHp/4hcgK2u/zfwL2w1c1vsScfdACnNGJNjK15jzEIRqY09aXfGbtR1InKJX3z+sfiX+TcOCmadAJKMMb7b+X5zdpKlTixfOfOYi61a8rcnl3l/BTwtIk2Bctgv+1fYasg0ING17DDsNfPuAeZzxDXNl0CglscHXP+fdL4MPr5tldP26IW9VPa2bceaLhxbtf+kqyy3bT0LW704BlvrchT7i6UU2TDGbBCR74CbRORxbNXntTnEm5OcEu9Tjs3fKez/uX0Ogfbj3H40ZJ7YmDQRGQV8KiLPBJjkJK7vriPQATjQMSGn44Tv73XYy3b+9rr+PxxgfNCc48PH2O/hw9iDfE3swTzozy0IvnUahj0G+Et0/X+q2ys77s8mjIzLcf7zOuk3nGmbikgn4DXs9hmJvXx4IU67lFxiCHYf8f8cfet2DbbmLLfp/W1zHf8yMcascO5yuRz7o+t5YJyItDPGnOr+dEqflXMi/xCbvIzFHufrY2vhctqWQe0/xpgBIjIJu26XARNEZKgxZlZ2M84xcTDG7BWRT4DhIvKsMcZ9MsDJuEphr12lYneMDa5JOvkN50Vb7JfVJwF7YgNbbbzbGJN+H7aI9DrdBRlj9gKvA6+LyCvYrLwJhbOe/ny//so4f1diM+vNxql3CuAY9gTr70vsyfgu4FtjzAkR+Qp73f4kGbUNvuUMBA4aY3Zms5yV2GrqJL/MPz/cgj2pTvErH4D99f+k/xty0AnbOPV9ABEpi73etz6X972EbZuwF/trzX27lS9ZCbSd3VZhk9Bn8zk2t/za/zcAA0Qk3JVgnHI/BMaYz5xb054k6wlmJ7YthVsL8m4Dtsq3njHm49wmzsWPzt/Lsdeu/bXG1m7c6dvvRaRlHpcZSCJ2e51ljJlTAPOHrLcOJwD/M8akiMhK7Am8hjHmi6xvzVFH7Il4vK9AsvYFkt1x6nT3kbXYE/CZxpgvTyHWoBhjDmHbfL0jIpOxiX474HOyX5f8cAE2cbrTVxsoIu39pgl0PAp6/zHGrMVuv8kiMgt7jJ2V3fTB/AK+FfthrBKRviLSREQaODvBSqChMeYI9prwIyJynYg0FJH7sLcW5rXK0meIs/yzRORh7A7+lDPuVyBeRIaISD0RGejEfcpE5FGn2reRiDTENshMBv4spPWsJCJVRaS6iCRgD747ycgYH8M2GJsnIm1E5EwRuUjsZZV6zjT/A6qKSIKIVBaRMgBOtdSv2F/sviRhNfbgfg2ZE4dXnfl8JCKXi+1n4gIRuVdEfL+8n8PuqO+JSEdnmg7ONmx3uhvAuTzQAXjZGLPe/cKezM92ftEE61egn4icKyLnY5PCYL7kvpPGg8B/jDHuk6CvJXdXETlDRMpnM49HgCtF5CkROc/ZrwaLvTSUl9j81y8/9v8Xse2OpotIY6eW7dHTmA/Y1u4J2Opdty+AS0Xkeuc4cg/2JJMnxt4G+xjwmIgMd7ZzUxHpIyITT3Fem7D7/wsi0l9E6otIaxEZ4UzyO/YX4Sjn+3cttiFgvnL2tweA0SJyj3PsPds5Pj2fT4u5QETud45lg7ANK6c6y9+ArTWY5Rx764vI+WL7XxiVy3x/BWqIyABnn/w/bJsGt/8BESLS1TlOlXXKv8B+Z3o4cT1IxmWTbDmXuSdhT37DnHPFOU7swezHcc6x1/0qB+Ac925wPoN62EvKx7E/Jn3r0tT53lQWkfysefoNe66+y9nfegL3+k2zxfl7lXM8ig1m/3HKHhOR9iJSR2xCkkDGJZrAgmkIgT2YTCHjFped2Ov6/clokR3sbYr9/coCNQpKwWkYQubbMb9xxm0BBvi95xFsi+7D2JqJG3A1XsOvEaTrfZnKsSeJ9dhkwXfrnPtOh9Ndzy/IptGl33r6Xiexv3IXAOf4TXsutmHYPmzV9iZsFWAlV4yvYX8pG2Cc673PO2XNXWXvOGU1/ZYThz2Z+NY1yYnH/d462IPsLmxtzFbs5Z0zs9vu2Gpdg18jJtf4p51lZdcQaQUZt3x9Q+6N7c7FJl5HnX3nVv/PI9B8nPKp2Ms4NQOMG+PEmUbOt2N2wdZaHXX2qa+xv4zzFFte9/9AnwP2stU657Ncj20Tk2V/zuk75Cp/znmvu3FkJDbh34mtvn4eW53t/rxmAV/4zWumbxu7yqbhasDqlA3BJsMp2O/HD2S+82ILAe4S8y934nzEKT+GbSjsbth9G/YX3VFs+wbf7YWd/b7P+XE7Zi9nPY5iGzP/CNzrGh+oMeF3OLc3ZrcfkXE75lzs8W43NvkKc00Tga2R/M3ZDruw+6+vIbavcWQvv2UJ9vbQXc4++T72GG7I3IDvBWcaQ8btmKWc/WKX8xk+hT3mZrkdM5vtNQz76zkVewxcRs4NDQPdzuh7TXamuR17m+Ih5/UD0NU1jzOwbcYOkvvtmKP9lr8NuCeX/fEuZ7qjzufoa7Tc1jXNA9g7SE6S+XbMbPcfoBbwHhm3mSZhj/nlctpffffrF1liryv9D3s3xKk21FIqT0RkPrbl9tVex6JUfhKRHdgTY5HuVE0VPcW9u2ilCoSIVMRWoXfHNhhSSimFJg5KZecn7KWaJ4wx33gci1JKFRlF/lKFUkoppYqOEv9YbaWUUkoFTy9VeKxy5cqmbt26XoehlFLFyqpVq3YbY+Jzn1LlN00cPFa3bl1WrlzpdRhKKVWsiEig3iFVIdBLFUoppZQKmiYOSimllAqaJg5KKaWUCpomDkoppZQKmjaOVEoFdPLkSbZt28bhw3l6CrVSpyUyMpIzzjiD2NhYr0NRfjRxUEoFtHv3bkSERo0aERamlZOq8BhjOHr0KElJSQCaPBQxejQorrZtgxtvBP01qArI/v37qVKliiYNqtCJCGXKlKFGjRrs3LnT63CUH61xKI4+/xz69oXduyEtDWbPBhGvo1IhJi0tjcjISK/DUCVY6dKlOX78uNdhKD/6U6I4SkqySQPA3Lnwn/94G48KWaIJqfKQ7n9FkyYOxdHgwfYyhc/w4bB6tWfhKKWUKjk0cSiunnsOzj3X/p+aCtddBwcOeBuTUiFg3Lhx9O/f/7Teu3jxYho1apTnGOrWrcsXX3yR5/nkJiYmhj/++KPAl6NCiyYOxVWZMvDWWxATY4c3bYIhQ0Afk65KkM6dO1OxYkVSU1M9Wb6IsGnTpvThjh078uuvvxboMgcPHkypUqWIiYlJfzVr1izX93Xu3JmZM2dmKktOTqZevXr5HuOsWbPo0KFDvs9XFQ2aOBRnjRqB+0Dwzjvw7LPexaNUIdqyZQuLFy9GRHj//fe9DqdQjRkzhuTk5PTXmjVrvA5JlSCaOBR3vXvDbbdlDI8eDT/84F08ShWSOXPm0LZtWwYPHszs2bMzjRs8eDC33XYbV111FeXKleOCCy5g8+bN6eNHjBhBrVq1iI2NpWXLlixevDjgMq666iqe9UvGzzvvPN577z06deoEQLNmzYiJieHNN9/km2++oWbNmunTJiYm0qNHD+Lj44mLi2P48OEAbN68mYsvvpi4uDgqV65Mv3792L9/f563SUpKCv379ycuLo4KFSrQunVr/v77b+6//34WL17M8OHDiYmJSY/DXWMyePBgbr31Vq688kpiYmJo3749O3bs4M4776RixYqcffbZ/PTTT+nLevzxx6lfvz7lypWjSZMmLFiwAICNGzcybNgwli1bRkxMDBUqVAAgNTWV0aNHU7t2bapUqcKwYcM4evRontdZFT5NHELBlCnQsqX9//hxuP562LPH25hU6BEpvFcQ5syZQ79+/ejXrx8LFy7k77//zjT+9ddfZ+zYsezbt48GDRpw//33p49r3bo1q1evZu/evfTt25frrruOlJSULMsYNGgQ8+bNSx9es2YNSUlJdO3alW+//Ta9LDk5md69e2d6b1paGt26daNOnTps2bKFpKQk+vTpA9gOju699162b9/Oxo0bSUxMZNy4cUGtd05mz57NgQMHSExMZM+ePUybNo3SpUvz6KOP0rFjR5577jmSk5N57rnnAr5//vz5TJgwgd27dxMVFUVCQgItWrRg9+7d9OrVi5EjR6ZPW79+fRYvXsyBAwcYO3Ys/fv356+//qJx48ZMmzaNhIQEkpOT0xOiu+++m99++43Vq1ezadMmkpKSePjhh/O8zqrwaeIQCqKibHsHJ7Pnzz9h4EA4edLbuJQqIN999x1bt27l+uuvp2XLltSvX5/XXnst0zQ9evSgTZs2RERE0K9fP1a77jzy/SqPiIhg1KhRpKamBmyb8I9//IPff/+d33//HYC5c+fSu3dvSpUqlWuMy5cvZ/v27UyaNImyZcsSHR2dft2/QYMGXHbZZURFRREfH8/IkSNZtGhR0Os/efJkKlSokP4aNGgQYLtp3rNnD5s2bSI8PJyWLVueUq+L3bt3p2XLlkRHR9O9e3eio6MZOHAg4eHh9O7dO1ONw3XXXUf16tUJCwujd+/eNGzYkOXLlwecrzGGGTNmMHXqVCpVqkS5cuW47777eOONN4KOTRUdmjiEijPPtB1B+Xz8Mfz7397Fo1QBmj17NpdffjmVK1cGoG/fvlkuV1StWjX9/zJlypCcnJw+PGXKFBo3bkz58uWpUKECBw4cYLevbxSXqKgorr/+eubNm8fJkyd5/fXXGTBgQFAxJiYmUqdOHSIisvazt3PnTvr06UONGjWIjY2lf//+AZefndGjR7N///70l2/dBwwYQJcuXejTpw/Vq1dnzJgxp9SBUpUqVdL/L126dJZh9zacM2cO559/fnrysn79+mzXYdeuXRw5coSWLVumT3/FFVewa9euoGNTRYcmDqHkmmtsGwefBx+ETz7xLh4VWowpvFcOjh49yvz581m0aBFVq1alatWqTJ06lTVr1gTVSHDx4sVMnDiR+fPns2/fPvbv30/58uUx2Sx30KBBvPrqq3z55ZeUKVOGhISEoDZXrVq1+PPPPzlx4kSWcffeey8iwtq1azl48CDz5s3LdvmnIjIykrFjx/Lzzz+zdOlSPvzwQ+bMmQPkb2dKW7du5eabb+a5555jz5497N+/n3POOSd9HfyXVblyZUqXLs2GDRvSk50DBw5kSkRU8aGJgx8RGS4iK0UkVURmucrbisjnIrJXRHaJyFsiUs01XkRkoojscV5PiBfdnj32GFx4of3fGNs1tatRmFLF3XvvvUd4eDg///wzq1evZvXq1WzcuJGOHTumnyRzcujQISIiIoiPj+fEiRM8/PDDHDx4MNvpExISCAsLY9SoUVlqG6pUqZJtPwht2rShWrVq3HPPPRw+fJiUlBSWLFmSHoOv4WBSUhKTJk06hS2Qva+//pp169aRlpZGbGwskZGRhIeH5xrrqTp8+DAiQnx8PACvvPIK69evTx9fpUoVtm3bxrFjxwAICwvj5ptv5q677kp/9kRSUhILFy7Ml3hU4dLEIavtwATgZb/yisB0oC5QBzgEvOIaPxS4FmgGnAd0A24p4FizioyE+fPB17J7/37o3l0fhqVCxuzZs7nxxhupXbt2eo1D1apVGT58OK+++mrAX/huXbp04corr+Sss86iTp06REdHU6tWrRzfM3DgQNatW5elY6hx48YxaNAgKlSowPz58zONCw8P54MPPmDTpk3Url2bmjVr8uabbwIwduxYfvzxR8qXL89VV11Fjx49TmkbPPHEE5n6cfBdstmxYwe9evUiNjaWxo0bc+GFF6bHPGLECN5++20qVqzIHXfccUrL89ekSRNGjRpFQkICVapUYd26dbRv3z59/MUXX0zTpk2pWrVqemwTJ06kQYMGtG3bltjYWC699NIC7/NCFQzJj+qxUCQiE4CaxpjB2YxvASwyxpRzhpcCs4wx053hIcDNxpi2OS2nVatWZuXKlfkaOwDLl0PHjuBk/PTpA6+9pg/DUkHbuHEjjRs39jqMImHOnDlMnz6d7777zutQSpzs9kMRWWWMaeVBSCWe1jicvk7ABtdwU8B9gXWNU5aFiAx1LoesLLDGQW3awAsvZAy/8QZMnVowy1IqhB05coQXXniBoUOHeh2KUkWCJg6nQUTOAx4C/uUqjgHcD4s4AMQEaudgjJlujGlljGnlu0ZYIIYMgVtcV0vGjIGvvy645SkVYhYuXEh8fDxVqlShb9++XoejVJGQ9T4hlSMRaQB8Aowwxri7m0sG3DdMxwLJxutrQU8/DWvWwPffQ1qa7Rxq1SqoXdvTsJQqDrp06cJhbR+kVCZa43AKRKQO8AXwiDFmrt/oDdiGkT7NyHwpwxtRUfYZFr77sXfvhh49QLt6VUopdRo0cfAjIhEiEg2EA+EiEu2U1QC+Ap43xkwL8NY5wEgRqSEi1YFRwKxCCzwn1avD22+DryOaVav0SZpKKaVOiyYOWT0AHAXuAfo7/z8A3ATUA8aKSLLv5XrfS8AHwDpgPfCRU1Y0dOhgL1v4vP669iyplFLqlGni4McYM84YI36vccaY8c7/Me6X633GGDPGGFPJeY3xvH2Dv3/+E4YNyxi+/35wnminlFJKBUMTh5JEBJ55Bjp3zigbMMA2nlRKKaWCoIlDSRMZads71Ktnhw8fts+4cLqBVUrl3eLFi2nUqFGe51O3bl2++OKLfIgoZzExMfnWHbUKfZo4lERxcfDBB+B73O6ff9o7LVJTvY1LqSDVrVuXKlWqZLpVcubMmXR216Z5qGPHjgXenfLgwYMpVapUpq6nmzVrluv7OnfuzMyZMzOVJScnU8/3YyIfzZo1K/1R4ip0aOJQUjVpYnuTDHN2gSVLbPuHItYsQ6nsnDhxgqfdDX6LiNyelZGfxowZQ3JycvormKeDKpVXmjiUZFdeCe6n8s2aBVOmeBaOUqfiX//6F5MnT2b//v1Zxm3ZsgURyXQSd//SnjVrFu3bt+euu+6iQoUK1KtXj6VLlzJr1ixq1arFGWecwezZs9Pfm5qayujRo6lduzZVqlRh2LBhHHX6Qvnmm2+oWbMmEydOpGrVqtx4443pZT6JiYn06NGD+Ph44uLiGD58OACbN2/m4osvJi4ujsqVK9OvX7+A63OqUlJS6N+/P3FxcVSoUIHWrVvz999/c//997N48WKGDx9OTExMehwiwqZNmwBbk3Hrrbdy5ZVXEhMTQ/v27dmxYwd33nknFStW5Oyzz+ann35KX9bjjz9O/fr1KVeuHE2aNGGB0+B648aNDBs2jGXLlqU/CTS3bamKB+05sqS76y5Yvx5ecR70OWaMbf9wik/rU6FPxhfeA9LM2FJt/8AAACAASURBVNxrvlq1akXnzp2ZPHkyEyZMOOVl/PDDD9x0003s2bOHsWPH0qdPH66++mo2bdrEokWL6NmzJz179iQmJoa7776bP/74g9WrVxMZGUnfvn15+OGH+bdzS/OOHTvYu3cvW7du5eTJk/zwww/py0lLS6Nbt25cfPHFzJ07l/DwcHwPtjPGcO+999KpUycOHjxIz549GTduHE899dQpr4/b7NmzOXDgAImJiURFRbF69WpKly7No48+ypIlS+jfvz833XRTtu+fP38+CxcupGnTpnTt2pWEhATGjx/PlClTGDt2LCNHjuRrp/v6+vXrs3jxYqpWrcpbb71F//792bRpE40bN2batGnMnDkz08PBctuWqujTGoeSTgRefNH28wD2UkW/fuA68ClVVD388MM8++yznM7D4s4880xuvPFGwsPD6d27N4mJiTz00ENERUVx+eWXU6pUKTZt2oQxhhkzZjB16lQqVapEuXLluO+++3jjjTfS5xUWFsb48eOJioqidOnSmZazfPlytm/fzqRJkyhbtizR0dHp1/0bNGjAZZddRlRUFPHx8YwcOZJFixYFvQ6TJ0+mQoUK6a9BgwYBEBkZyZ49e9i0aRPh4eG0bNmS2NjYXOaWoXv37rRs2ZLo6Gi6d+9OdHQ0AwcOTN9W7hqH6667jurVqxMWFkbv3r1p2LAhy5cvDzjfYLalKvq0xkHZbqkXLIB27eD33yElBa6+2j7fogAaTCmVX8455xy6devG448/fsqPAK/i64Yd0k/2/mXJycns2rWLI0eO0LJly/RxxhjS0tLSh+Pj44mOjg64nMTEROrUqUNERNbD7c6dO7njjjtYvHgxhw4d4uTJk1SsWDHodRg9enTA2pYBAwaQmJhInz592L9/P/379+fRRx8lMjIyqPn6b4dA28Vnzpw5PPnkk2zZsgWwDS13794dcL7BbEtV9GnioKzKleHjj6FtW9izB3btgq5dYelSqFTJ6+hUERDM5QMvjB8/nhYtWjBq1Kj0srJlywL2kdi+X9o7duw4rflXrlyZ0qVLs2HDBmrUqBFwmgAPwU1Xq1Yt/vzzT06cOJElebj33nsREdauXUtcXBzvvfdeeruDvIiMjGTs2LGMHTuWLVu20LVrVxo1asSQIUNyjPVUbd26lZtvvpkvv/yShIQEwsPDOf/88/H1fee/rGC2pSr69FKFytCgAbz/vq2BAPj1V+jeXW/TVEVagwYN6N27N88880x6WXx8PDVq1GDevHmkpaXx8ssvs3nz5tOaf1hYGDfffDN33XUXO53+TpKSkli4cGFQ72/Tpg3VqlXjnnvu4fDhw6SkpLBkyRIADh06lN5wMCkpiUnuxsp58PXXX7Nu3TrS0tKIjY0lMjKS8PBwwNYm5FefDYcPH0ZEiI+PB+CVV15h/fr16eOrVKnCtm3bOHbsGJD3bamKBk0cVGbt2sGcORnD334L//d/epumKtIeeuihLI+/njFjBpMmTSIuLo4NGzbQrl27057/xIkTadCgAW3btiU2NpZLL7006H4awsPD+eCDD9i0aRO1a9emZs2avPnmmwCMHTuWH3/8kfLly3PVVVfR4xQbJT/xxBOZ+nGoXLkyYGtXevXqRWxsLI0bN+bCCy+kf//+AIwYMYK3336bihUrcscdd5zS8vw1adKEUaNGkZCQQJUqVVi3bh3t27dPH3/xxRfTtGlTqlatmh5bXralKhqkqD1OoaRp1aqV8bWwLlKeeALuvjtj+IEH4JFHvItHFbqNGzeecrsBpfJbdvuhiKwyxrTyIKQST2scVGD/+hfcckvG8IQJ4NfbnFJKqZJHEwcVmAg89xxccUVG2S232DYQSimlSixNHFT2IiJg/nxo3twOnzwJvXvb7qmVUkqVSJo4qJyVK2dv0/T155CSAt26wYYN3sallFLKEyGTOIhImIi0FpGbRGS087e1iITMOnqmalVYuBDOOMMO798PXbrYp2qqkKaNp5WXTp486XUIKoBif1IVkfIiMgFIAt4GrgPaOn/fAZJEZIKIVPAwzOKvQQP45BOIibHDSUk2edizx9u4VIGJjo5mz549mjyoQmeM4dixYyQlJaV35qWKjlDoOfIHYB7Q3hiTpVcTETkT6AcsA/Tesrxo0QLee88+VfP4cfjlF3vZ4osvQL/cIadmzZps27bttJ4DoVReRUREUL58+fT+H1TRUez7cRCR0saYXJ/JGux0ha3I9uOQk/nzoU+fjE6huna1CUWQ/eArpVReaT8O3in2lyqMMUdFpJSINAw0XkQaikipopg0FFvXXw9PP50x/PHHMGgQ6INqlFIq5BX7xMHxL2BENuNuB0ZlM06drttvh/vvzxh+/XUYNky7plZKqRAXKonDDcAT2YybhG3joPLbI4/AbbdlDM+cCaNGafKglFIhLFQSh5rGmID3BhpjEoFahRxPySACzzwDAwdmlE2dCuPHexeTUkqpAhUqiUOKiAR8uLuI1ARSCjmekiMsDP7zH+jZM6Ns/HiYMsW7mJRSShWYUEkcPgUezmbcWGe8KigREfDqq5mfazF6NEyf7l1MSimlCkSoJA73A5eLyHIRuVtEBjp/fwC6APcFOyMRGS4iK0UkVURmucpLicjbIrJFRIyIdPZ7n4jIRBHZ47yeEBHJp/Ur+qKi4J13oFOnjLJhw2xCoZRSKmSEROJgjEkCWgDfAn2ACc7fb4GWzvhgbXfe/3KAcd8B/YEdAcYNBa4FmgHnAd2AWwJMF7rKlIEPPoDWre2wMbb9w5tvehuXUkqpfBMKPUcCYIzZBYzOh/m8CyAirYCarvJjwFPOuEAdFgwCphhjtjnTTAFuBqblNaZiJTbWdk3duTOsX2+fqNmvn20Lcd11XkenlFIqj0KixqGIaAqscQ2vccqyEJGhzuWQlSHZnW9cnO2GurHTw3daGtxwg72UoZRSqljTxCH/xAAHXMMHgJhA7RyMMdONMa2MMa3i4+MLLcBCVaUKfPUVnH22HU5Ls91UL1jgbVxKKaXyRBOH/JMMxLqGY4FkU9wfBpIXVava5KFRIzt84oTtrvq///U2LqWUUqdNE4f8swHbMNKnmVNWslWrZpOHhs6jRE6csG0d3n/f27iUUkqdlhKROIhI91OYNkJEooFwIFxEokUkwhkX5YwDKOWM812KmAOMFJEaIlId+3yMWfm3FsVY9erw9dfQoIEdPn4cevWyd2AopZQqVkpE4oB9XkWwHgCOAvdgb7086pQB/OoM1wAWOv/Xcca9BHwArAPWAx85ZQqgRg2bPNSvb4ePH4cePeCtt7yNSyml1CmRknwJviho1aqVWblypddhFJ7ERLjoIti82Q6HhcHs2dC/v7dxKaWKFRFZZYxp5XUcJVFJqXFQRUWtWrBoUUaDyZMnbSdRM2Z4G5dSSqmghEwHUAAi8jkQqAolFdgGvGuM+bxwo1JZ1Khhk4fLLoN162wPk0OHwtGjcMcdXkenlFIqB6FW47ASOAfb+dKHZHTC9As2eXhTRO70LjyVrkoV2+ahZcuMshEjYOJE72JSSimVq5CqcQAuBLoYY9b5CkRkNjDdGNNOROZj7354yqsAlUtcHHz5JXTtCkuX2rJ77oEjR2DcOChBzwhTSqniItRqHJpi73xw+90pxxizFKha2EGpHJQvDwsX2mdb+Dz8MIwcads/KKWUKlJCLXFYCkwXkZrOY65rAi8AywBEpAnwt5cBqgBiYuDjj+GKKzLKnnoKBg2yt20qpZQqMkItcRgIVAC2AieALUAl7JMrAQQY4ElkKmelS8N770HPnhll8+ZB9+720oVSSqkiIaQSB2PMLmPMtUAZoB5Q1hjT3RjztzN+g3O5QhVFUVHw5pv2Dgufjz6CLl1g/37v4lJKKZUupBIHABGpB9wB3GGMSRWRBiLS2Ou4VJDCw2HaNLjvvoyy776DCy+Ev/7yLi6llFJAiCUOIvIPYAW2MeTNTnEc8LRnQalTJwKPPgpTpmSUrV0LHTpk9DiplFLKEyGVOACPAZcbYwYDaU7ZT8D5nkWkTt/IkbY76vBwO/zHH9CuHSxf7m1cSilVgoVa4lDNGLPK+d/Xg+RJbKNIVRwNHAgLFkC081DSnTvtZYt33vE2LqWUKqFCLXH4SUT6+JX1wPYoqYqrq6+Gzz+HSpXscEqKfSz3xIm2u2qllFKFJtQShxHAkyKyECgjIv8FJgMjvQ1L5VmHDvD999CwYUbZPffYGomUFO/iUkqpEiakEgdjzHqgEfAq8CiwADjPGLPR08BU/mjYEJYtg06dMsrmzbOXLrZv9y4upZQqQUIqcQAwxhwyxswxxjxijJlljNEOAEJJXJy9bDFkSEbZ8uXQqpU2mlRKqUJQ7B9ylcOjtDMxxlxeCOGowlCqFMyYAc2awV13QVqa7eOhUydbPkA7B1VKqYISCjUObwPvOK+VwHnAOuxjtdcC52L7dlChRARuv90+IKtiRVuWmmrbPIwZY5MJpZRS+a7Y1zgYY17y/S8iS7CP1V7jKpsDvOhFbKoQXHIJrFgB11wDP/9syyZNgvXr4bXXoEIFb+NTSqkQEwo1Dm7nAP4NIX/Deay2ClH169tGk1dfnVH2ySe23cPatd7FpZRSISjUEoclwEwRqeU8Vrs2trZhmcdxqYIWG2ufrul+xsXmzdC2Lbz6qndxKaVUiAm1xGEQEIt9nPYJ4A/sY7YHexeSKjRhYfYZF2+9BTExtuzoUejf37aHOHbM2/iUUioEhFTikNtjtVUJ0auXvTWzUaOMsueeg4su0v4elFIqj4p94iAilfzLjDGpxpitxpjUnKZTIaxxY5s89OyZUbZ0KbRoAYsWeReXUkoVc8U+cQCWiMiTItI80EgROV9EngS+K+S4lNdiY+1liyeesJcxAP7+296JMXkynDzpbXxKKVUMhULi0ALYDswXkf0iskJEvnL+7gPeABKd6VRJIwL/+pftbTI+3palpdmya66BPXu8jU8ppYqZYp84GGOOGmMmG2MaAq2BfwPzgMeAVsaYs40xU40xQT0JSUSGi8hKEUkVkVl+4y4RkV9E5IiIfC0idVzjREQmisge5/WEiOjjvIuKiy+GVavgggsyyj76CM4/H5Ys8S4upZQqZop94uBmjPndGPOuMeZlY8wCY8zm05jNdmAC8LK7UEQqA+8CDwKVsL1UvumaZChwLdAM23tlN+CW01i+Kii1asG339puqn22bbMPyXr8cb10oZRSQQipxCE/OInHe4B/HXYPYIMx5i2n9mIc0ExEznbGDwKmGGO2GWOSgCnobaBFT6lS8OST8N//ZnRVnZYG994LV10Fu3Z5G59SShVxmjgErymQ3pW1MeYwsJmMXikzjXf+D9hjpYgMdS6HrNylJypvXHMNrF4NCQkZZZ9+ai9dfPutd3EppVQRp4lD8GKAA35lB4By2Yw/AMQEaudgjJlujGlljGkV72uwpwpf7dr21swxYzLKtm+3/T2MHw8nTngXm1JKFVEhmzgUQL8NydheKd1igUPZjI8Fko0xuT7yW3koMhImTrQNJePibNnJkzBunH1M9x9/eBqeUkoVNSGVOIhIrIi8LCKHsbdgIiLdROTBfJj9BmzDR9+yygL1nfIs453/N6CKh65d7aWLjh0zypYts5cu5swBzf+UUgoIscQBeN75ew7gezDBcmBAsDMQkQgRiQbCgXARiRaRCGABcI6I9HTGPwSsNcb84rx1DjBSRGqISHVgFDArz2ukCk/NmvD11/Z5FxHOE+cPHYJBg6BPH9i3z9v4lFKqCAi1xKELMMwY8z/AABhjdgJVTmEeDwBHgXuA/s7/DxhjdgE9gUeBfcAFQB/X+14CPgDWAeuBj5wyVZyEh9snbC5dCg0bZpTPnw/nnQfffONZaEopVRRIKF2CF5HNQDtjzN8istcYU0lEagJfGWPO8jq+QFq1amVWrlzpdRgqkORk2+fDzJkZZb6eKB95xN7aqZTyhIisMsa08jqOkijUahxmYbueTsB25tgc25HTDE+jUsVTTAzMmAHvvpvRcNIY++yL1q1tmwillCphQi1xeBT4BHgVe3vku8BXwGQvg1LFXPfusHYtXHZZRtnatTZ5ePhhOH7cu9iUUqqQhUziICLh2DYHU40x9YwxkcaYM40xj+stkSrPqle3HUQ99RSULm3LTpyAsWOhbVtYv97b+JRSqpCETOJgjEkDXjTGpHodiwpRYWEwYkTWHid//BFatoR//1s7jVJKhbyQSRwcn4hIF6+DUCHurLNg8WKYNAmiomzZsWP2boz27WHjRm/jU0qpAhRqiUMa8J6IfCoiM0Rkuu/ldWAqxISHw+jRtrahdeuM8uXLoXlz+7RNbfuglApBoZY4/Ak8BfwE7MY+4dL3Uir/NWli+3x49FHbfTVAaqp92mbr1qC32iqlQkxI9eNQHGk/DiFk3ToYPNjWQviEhdm+IMaPh7JlPQtNqVCj/Th4J6RqHESkXXYvr2NTJcC558IPP9h+HqKjbdnJkzBlih33+efexqeUUvkgpGocROQvv6JKgAC7jTHVPQgpV1rjEKI2b4ZbboEvv8xcPnAgPPlkRodSSqnTojUO3gmpGgdjTDX3C9sJ1GTgMY9DUyVN/fq2huGVV6BixYzyOXOgcWN94qZSqtgKqcTBnzHmOPYplvd7HYsqgURsm4eNG+3TNX127bJP3LzwQu04SilV7IR04uDo6HUAqoSrUgVefx0+/BBq1cooX7wYzj/f3tZ56JB38Sml1CkIqcRBRH4Xkd9cr23Yx1s/5HVsSnHVVbb24e67ISLClqWl2caTZ59tH92tly+UUkVcqDWO9O818jCw0RhTZPtx0MaRJdTPP8Ott8KiRZnLL7sMnnvO9k6plMqWNo70TkjVOAANjTELXa/vjDF7RORWrwNTKpMmTeDrr2HePHspw+fzz+2tm/feC8nJ3sWnlFLZCLXEIbu7JyYUahRKBUME+vWDX36B22+3nUWBfe7F449Do0bw6qt6+UIpVaSEROLg6uQpTEQS/Dp/6g/oTzdVdFWoAM88AytW2Ed0+2zfDv37Q4cOsGqVd/EppZRLSCQOwDvOqzTwrmv4bWCk81KqaGvRApYsgdmzoWrVjPKlS+1zL4YOtbdyKqWUh0IicXB1+PSOXydQ1Y0xLYwxb3sdo1JBCQuzvUv+9huMGZPx4CxjYMYMaNgQnn5an7yplPJMSCQOPsaY672OQal8Ua4cTJxoO4jq2jWj/MABuPNO2/+Df3fWSilVCEIqcRCRsiLymIgsEZFf3X06eB2bUqflrLPgo49s51ENGmSU//wzXHop9OwJW7Z4Fp5SquQJqcQBeB7oDDwDVAceBPYBMzyMSam8u+oqW/swcSLExGSUv/uuffbF2LFw5Ih38SmlSoxQSxyuBK41xrwJpDl/rwP0EoYq/qKibLuHX3+FAQMyylNS4OGHtfdJpVShCLXEIQLw9RKZLCKxwDagkXchKZXPqle3T9dcuhRatswoT0yE3r3hootg7Vrv4lNKhbRQSxzWkvFQq6XAU85rs2cRKVVQEhJg+XKYORPi4zPKFy2C5s3htttgT5HtbV0pVUyFWuIwDPjL+X8EUAqoAwzOrwWISGMR+UpEDojIJhHp7hp3iYj8IiJHRORrEamTX8tVKqCwMBgyxN6+edddEB5uy0+ehBdesI0rX3wRTpzwNk6lVMgImcRBRMKBnsAWAGPMX8aY/saYfxhj1uTTMiKA/wIfApWAocA8ETlLRCpjO5960Bm3EngzP5arVK4qVIAnn7SXKC69NKN87177MK2mTeHNN21CoZRSeRAyiYMxJg0YDRwrwMWcjb1bY6oxJs0Y8xWwBBgA9AA2GGPeMsakAOOAZiJydgHGo1RmTZrAZ5/BggVw5pkZ5b/9Bn362N4pP/pIG1AqpU5byCQOjteAGwtw/pJN2TlAUyC9ZsMYcxjbtqJpljeIDBWRlSKycpd2Iazymwhce63t6+Gxx6B8+Yxxa9ZAt272+Rf+j/RWSqkghFri0BiY5rQz+FxEPvO98mn+vwA7gX+JSKSIXA5cCJQBYoADftMfAMr5z8QYM90Y08oY0yre3ahNqfwUHW0fz/2//9m/ZcpkjFu6FDp3hssvh+++8yxEpVTxE2qJw3zgdmAq9gFX77heeWaMOQ5cC1wF7ABGOcvchn0CZ6zfW2KBQ/mxbKVOW8WKtuZh82YYPjzj+RcAn38OHTvaWzi/+kovYSilciVGDxR5IiJLgdmAAQYZY9o75WWBXUALY8wv2b2/VatWZuXKlYUSq1KA7aJ6/HjbF4R/Y8l27eDBB6FLF3vJQ6kiSkRWGWNaeR1HSRRqNQ6IyAAR+VhEVjjD7d23TObD/M8TkWgRKSMio4FqwCxgAXCOiPQUkWjgIWBtTkmDUp6oWxdeeQU2boTBgzNu4QR7CePKK6FNG9vAMi3NqyiVUkVUSCUOIvIQcCf28kFDp/hv4P58XMwAbF8RO4FLgMuMManGmF3Y20EfxT4f4wKgTz4uV6n8ddZZNoH4/XcYOjTzJYyVK6FHD/scjGnT4OhR7+JUShUpIXWpQkS2Am2MMX+LyD5jTEUREWCPMaaS1/EFopcqVJGRmAhPPAEzZkBqauZxlSvbnihvuy1zL5VKeUQvVXgnpGocsD1F7nf+92VEZQF9bKBSualVC5591t6Fcc89mW/j3L3btouoXRv++U/bL4RSqkQKtcThM2Ci04ukz0PAxx7Fo1TxU60a/PvftgZi6lSbLPikpNhLF40a2QaU77+v7SCUKmFCLXG4EzgL239CrIjsA84FxngalVLFUblycOed9jbO11+3vU66ffYZ/OMfUL8+PP44aGdmSpUIIdXGwUdEamMfbpVojNnicTg50jYOqtgwxvY2+eST8OGHWft8iIqC66+37SDatNHbOVWB0jYO3gm1GgdEJAZoDZwPtBSRLD03KqVOg4jtbfL9920txJgxEBeXMT41FebOhbZtoVkzePpp2zZCKRVSQqrGQUQ6Au9he3LcCtQGagHXGmMWexlbdrTGQRVrKSn2qZvPPw8rVmQdX6qUvZzxf/8Hl12Wuc8IpfJAaxy8E2o1Di8CdxpjmhljrjHGnA+MAKZ5HJdSoSk6GgYNguXL7WvQIChdOmP8sWPw1lu2U6m6dW2vlJs2eRauUirvQq3G4QBQyXnEtq8sHNhrjCmf/Tu9ozUOKuQcOGBrIV5+GX74IfA0rVtDv362TUS1aoUbnwoJWuPgnVCrcXgduMmvbAj2cdtKqcJQvrztifL772H9ehg5MmunUStW2Ds2ata0lzBeecUmHEqpIi/Uahy+BDoBf2LbOdTEtnH4Fkh/mo8x5nJPAgzgdGscDqYepFypcoi2XFfFwbFj8NFHMHs2fPwxHD+edZqoKHtJo2dPuPrqzB1QKeVHaxy8E2qJwy3BTGeMeamgYwnW6SYOl829jORjydzX4T66ndVNEwhVfOzbB++8A6++am/vDHQMioyESy+1z8v4xz+0m2uVhSYO3gmpxKE4Op3EYUXSCtrMbJM+fO4Z53JPh3u4vun1RIRF5HeIShWcbdtse4hXX4Wffgo8TVgYXHihTSK6d4caNQo3RlUkaeLgnZBLHESkNdAciHGXG2Oe9CainJ1O4vDSype449M7OJZ2LFN5vYr1uLv93QxqNoioiKj8DFOpgvfrr7Ym4t13YdWq7Kdr2xauuca+mjTRjqZKKE0cvBNSiYOITAFuBJYC7ucAG2PM9d5ElbPTvVSRdDCJJ5c9yUurXuLw8cOZxlWLqcaohFHc0uoWYkrFZDMHpYqwLVtgwQKbSCxdGvhyBsCZZ9oE4uqroVOnzI8GVyFNEwfvhFrisBdoZoxJ9DqWYOX1dsw9R/bw7PJneeaHZ9iXsi/TuIrRFbnjgju4vc3txJWJy2YOShVxf/0F771nk4hvvsn+oVrly8MVV9hE4soroWLFQg1TFS5NHLwTaonDOqCjMWZ/rhMXEfnVj8Oh1ENMXzWdKcum8FfyX5nGlY4ozeDzB3NX27toGNcwz8tSyjN79sAnn9hurz/9FA4dCjxdeDh07AjdutkkonFjvaQRYjRx8E6oJQ5tgVHAbGCne5wxZrknQeUivzuASjmRwuzVs3li6RP8se+PTOME4ZpG1zAqYRQdanfQOzFU8XbsmL0r4/337evPP7OftnZtWxtxxRVwySUQG1t4caoCoYmDd0ItcbgReNYZ9G/jcIYHIeWqoHqOPHHyBPM3zGfS0kms3rE663Krt2JUwih6Nemld2Ko4s8YWLfOJhAffGC7v85ORAS0b2+TiCuvhPPO09qIYkgTB++EWuKwC7jRGPOh17EEq6C7nDbG8PWWr5mybAof//5xlvG1y9dmxAUjuKnFTcRG6a8wFSL++ss++vvTT+GLL+DgweynrVYNunSxScQll2R+4qcqsjRx8E6oJQ7bgDONMQG6pSuaCvNZFT/v+pmpy6Yyd+1cUtNSM42LKRXDoGaDuK31bTSOb1wo8ShVKI4fh2XLbBLx6afZ9xcBtubh/PNtAnHJJbadRNmyhRerCpomDt4JtcTh/4BzgfHFpYGkFw+52nl4Jy+seIHnVzzP7iO7s4y/tN6lDG89nG5ndSM8TB+DrELMjh3w2We2keVnn8HevdlPGxlp+43wJRIXXKC3fBYRmjh4J9QSh+NAOGAA3z1bgm3jUMqzwHLg5dMxjx4/yty1c3n6h6f5edfPWcbXKV+HW1vfypDmQ/R2ThWa0tJg5UpbE7FwoW0bkd3tnmBrHzp1ykgkzj3X3sGhCp0mDt4JtcShUXbjjDG/FmYswSoKj9X2tYN4bvlz/PfX/3LSnMw0Pjoimr7n9OW2NrfRoloLj6JUqhAcPAjffmvbRXz5pX26Z04qVIAOHWwy0akTtGihNRKFRBMH74RU4uAmIpWMGmt3awAAFxlJREFUMTnUQRYNRSFxcNu6fyvTVk5jxo8z2HN0T5bxraq3YmiLofQ5pw/losp5EKFShejvv+Grr2wS8eWXtkfLnJQtC+3aZSQSbdpAdHShhFrSaOLgnZBKHEQkFngK6A1gjCkrIt2A5saYRzwNLhtFLXHwSTmRwhvr3+DZ5c/y418/ZhkfUyqGvuf0ZWjLobSs3tKDCJXywB9/ZCQR33xjE4uclCpl20X4Eom2bbUPiXyiiYN3Qi1xmAscBx4BfjTGVBSRM4DvjDFneRtdYEU1cfAxxvD9tu95fsXzvP3z21nuxgBoXrU5Q1sOpe+5ffWWTlVyGAObNtlOqL791r62bs35PSJwzjm2ViIhwb4aNtR+JE6DJg7eCbXEYSdQ0xhzTET2GmMqOeUHjDHl82kZdYEXgAQgFXgbuNMYc0JELgGeB2oDPwCDjTE5HkmKeuLgtvfoXuaumctLq15i4+6NWcaXiSzDDefcwJDmQ2hbs632TKlKnq1bYfHijGTit99yf09cXEYS0a4dtG6tt4AGQRMH74Ra4rAZaGeM+duXOIhITeCr/KpxEJGPsd1ZDwMqAJ8DM4DXgM3ATcAH2FqPjsaYtjnNrzglDj7GGJYmLmX6j9OZv2E+KSdSskxzVtxZDGo2iAHnDaBW+VoeRKlUEbBjh00kvv3W/l23Dk6ezPk94eG2N8uEBHuZo3VraNQIwsIKJ+ZiQhMH74Ra4vAgcClwD/AxcDEwEfjcGDMpn5axERhljPnYGZ4ExAKrsDUM7ZzyssBubPuKX7KbX3FMHNz2Ht3LvLXzmL5qOht2bcgyXhAurXcpg5oNonvj7pSJLONBlEoVEYcOwYoVtkOqpUvh++9z7kfCp1w5e8dG69YZr7p1S/QlDk0cvBNqiUMYMAYYCtQCtgEvARNNPq2oiAwD2mFrHCoCC4EHgc5AKWPMP13TrgfGGmPeyW5+xT1x8DHGsGzbMl7+6WXe3PAmyceSs0xTrlQ5rm96PYPPH0z7Wu31UoZSxtjLGUuX2mRi2TLYsMGW5yYuDlq1ypxMVKtW8DEXEZo4eCckEgcRucEY83ohLasxMA9ohu1sajZwIzAT2GWMucc17RJghjFmlt88hmKTG2rXrt1ya24NqoqZw8cOs+CXBcxeM5sv//gSQ9Z9rH7F+vQ7tx83nHsDZ1c+24MolSqiDhyAH36wtRErVthXbndv+NSoAc2b226zfa8zzwzJyxyaOHgnVBKHg8aYAm/O79RobMHWYkwGYoCXgV+BKCDSGHOra/p1wLiSUOOQnT8P/MncNXOZvWY2v+/9PeA0zas2p++5fendtLe2h1DKnzGQlJSRRKxYYXu73B9kr/qxsdCsWUYi0bw5NGkCUVEFG3cB08TBO6GSOBwyxhR4b0QiUhnYBVQwxhxwyq4FJgDPAIOMMe2d8rLOtC1CuY1DsHyXMmavns0bG97gYGrgpxV2qtOJG865gV5NelG5TOVCjlKpYsIY2Lw5cyKxahUcORLc+yMibPLgSyaaNbO3iZ5xRsHGnY80cfBOqCQOR4CrsM+lCMgY81U+LesPYDoZNQ6vAEeAO4H/b+/eg+M6yzuOfx/d77JlyfIltuxYtmVycRgbcFyIPdgh40CcQGiBzgQoMIQWmOGfdlo6bUPpAB3aQpmQtuEyhUICxYnB4IQGG+wkJBBSHDt25atsWfFFsmXdvLpLT/84K3mlSPZK1u6Rdn+fmTN79uw5Z5+T9Ua/Pec973sc+AiwE/gcsCEV76q4Xl19Xew8tpPHXn2Mncd20jvQ+7p1sjKyeMeyd/CBmz/APSvuoTRvSu6mFUldAwNw5Ai88kow7dsXTM2v7wF2XBUVQYCInW66CUqn3/dPwSE8qRIcBoB6xg8O7u43TtF73UbQO+VqgoG0fgV80t2bzGwz8DBQxZV+HE5dbX/pGBxitXa3sr12O48ffJzdJ3e/bpwMgJzMHO688U7e+4b3snXlVsryy0KoVGQGcoezZ68EiaFQceLExPazaNHrA8WqVZCfn5i646DgEJ5UCQ5JaeOQCOkeHGKdv3yeHx36EY8dfIzfvPabMdfJyshi09JN3L/qfu6ruY+KwookVymSAtra4MCBK0Hi4MHgbo5IJP59mEFVFdTUjJxWroTKyoTfKqrgEB4Fh5ApOIytrqWOHx78Idtqt405VgZAhmWwcclG7l91P++ueTfzi9PnVjSRKTc4GPR8efDgyKm2Fvr6Jrav0tKxA8WyZcH4HVNAwSE8qRIcktI4MhEUHK6trqWOJ/7vCbbVbuOlMy+NuY5hrLthHfeuvJetK7dSU16jfiJEpkJfXzAmR2yYePXV4HLHtXrBHC0zMwgPy5fDihXB47veFVwKmSAFh/CkRHCYyRQcJuZ022merH2SJ2qf4Nenfz1mHxEA1WXVbF2xlXtr7mX9ovVkZWQluVKRFNfTEwSKI0fg8OGRU0dH/Pv5xS9g8+YJv72CQ3gUHEKm4DB5ZzvOsr12O9tqt/Fs/bNjNqwEKMsv453L38nWlVu5a9ldFOfOyJNTIjODezBGx+gwceTI2KOH1tfD4sUTfhsFh/AoOIRMwWFqNHc289Sxp/jJkZ/w8+M/J9I3diOvnMwc3r707dxdfTdblm+huqw6yZWKpLFIBI4duzIdPw7f/OakerZUcAiPgkPIFBymXnd/N3tO7WHHkR3sOLKDMx1nxl23uqyaLdVb2FK9hY1LNpKfHd7tZSISPwWH8Cg4hEzBIbHcnd+f+30QIo7u4JXzr4y7bl5WHhuXbBwOEsvnLE9ipSIyEQoO4VFwCJmCQ3LVt9bz1LGnePr40+w+uZvOvvG76F02exlbqrdw57I72bhkIyW5M/KOX5GUpOAQHgWHkCk4hKenv4fnTj83HCQOXxx3SBEyLZM3L3wzm2/czOYbN7PuhnXkZE7N/egiMnEKDuFRcAiZgsP0car1FE8fezqusxEF2QVsqNowHCRunnszGZZ6QxeLTFcKDuFRcAiZgsP0NHQ24pkTz7D75G72nds3bp8RABUFFWy6cROblm5iQ9UGqsuq1QGVSAIpOIRHwSFkCg4zw8XOi/zq5K/YVbeLXSd3UddSd9X15xfN546qO9hQtYENSzawqnyVgoTIFFJwCI+CQ8gUHGamupY6dtftZvfJYLrYefGq61cUVHBH1R3DYeKWylt0aUPkOig4hEfBIWQKDjPfoA9yoPEAu+p2sefUHp4//TxtPW1X3WZ23mzeVvU23rroraxftJ41C9aQl5WXpIpFZj4Fh/AoOIRMwSH1DAwOcKDxAHvr97K3fi/P1j/Lpa5LV90mOyObNQvWsP6G9axftJ7bF93OguIFSapYZOZRcAiPgkPIFBxS36APcqjp0HCI2Fu/l6ZI0zW3WzJrCbffcDvrFwVh4tbKWzVYl0iUgkN4FBxCpuCQftydI81HeLb+WV587UVeaHiBo81Hr7ldYXYhaxas4U0L3sTaBWtZu2Aty2YvU6NLSUsKDuFRcAiZgoNAcNfGiw1BiHjxtRd56cxLdPV3XXO7WXmzWDN/zXCQWLtgLVWlVQoTkvIUHMKj4BAyBQcZS99AH/sb9/NCwwvDU0N7Q1zbzsmfw23zbmN15WpurbyV1fNWs6p8FblZuQmuWiR5FBzCo+AQMgUHideZ9jP87uzvePnsy8NTc1dzXNtmZWRRU15zJUxEH+cVzdPZCZmRFBzCo+AQMgUHmSx3p76tnt+diYaJc0GYaO9pj3sfZfllrCpfFUwVVx4Xly5WPxMyrSk4hEfBIWQKDjKVBn2QupY6DjQeYP/5/RxoCh5Ptp6c0H7ys/JZWb5yOFTUlNdQXVbNsrJlGiVUpgUFh/AoOIRMwUGSob2nnYNNB4Mw0XiA/Y37ebXpVS73Xp7wvioKKqguqw6CxOxlw4GiuqyaOflzdOlDkkLBITwKDiFTcJCwDPogr7W/Ru2FWmov1lJ7oZbDzYepvVDLhc4Lk9pnaW4pS2YtYXHp4hFTVWkVi0sXM69oHpkZmVN8JJKOFBzCo+AQMgUHmY6aO5uHw0TtxVqONh/lRMsJ6lrq6B3onfR+szKyuKHkhuEgsaB4AfOL5jO/eP6I+YLsgik8GklFCg7hUXAImYKDzCQDgwOc6TjD8UvHOXHpRPDYcuVxMpc+xlKSWzIiSMwvCqZ5RfOYWziXuYVzqSyqpLygXL1ppikFh/AoOIRMwUFShbtzofMC9a31nG47zem209S3jZy/1iiikzEnf85wmJhbOJfKwsoRz2OnktwStcFIEQoO4VFwmAAzG/1zKh94xN0/HX19E/B1YDHwW+DD7l5/tX0qOEg66ezrpKGtYThMnO04y7nL54Kp4xxnO85y/vJ5+gb7EvL+uZm5zC2cS0VhBXPy5zCnYE7wGDtfMIfygvLh+eKcYoWNaUjBITw6xzcB7l40NG9mhUAj8KPo83LgSeBjwE+BzwM/BNYlv1KR6akgu4CV5StZWb5y3HUGfZBLXZc41xEEirMdZ4fnmyJNNEWaaIw00hRpormzGSf+Hz89Az00tDfE3QsnBCOXluWXjQgWc/KvhIuy/DJm5c1iVt4sZufPHp4vzS1VQ1BJSQoOk/deoAl4Lvr8PcAhdx8KEg8BF82sxt0Ph1OiyMyTYRmUF5RTXlDOLZW3XHXd/sF+mjubhwPF6GAxeor0RSZcT99gH42RRhojjRPetjineESYGA4YeeMvK80rpSS3hOKcYrIzsyf8niKJpuAweR8CvutXrvXcBOwfetHdI2Z2Irp8RHAws48DHwdYvHhxcqoVSUFZGVlUFlVSWVQZ1/qR3ghNkSYudl6kuauZ5s7mEfPNXc0j5zubJxU2hnT0dtDR28HpttOT2j4vK4/inGKKc4uHw0RJbknwPKdk7OXR57HzRTlF5GTm6JKLTAkFh0kws8XABuCjMYuLgNE3v7cBxaO3d/dHgUchaOOQoDJFZJTCnEKW5ixl6eylcW/T3d89IkjEPl7svEhLdwut3a0jppauFtp62q673u7+brr7uyfdr0asTMukMKeQwuxCCrILhudHLIt5Xphz7WUF2QXkZeWRn5VPXlaegkmaUHCYnA8Cz7t7bD++l4HRffGWAB1Jq0pEplxeVh4LSxaysGThhLYbGBygo7djOEiMCBZjhY3osrbuNjp6O2jvaWfQB6fsOAZ8gPae9gmNZTJRuZm55GfnDweJofn87PzhgDFiPiufz6z7DFWzqhJWk0w9BYfJ+SDwpVHLDhFcvgCGG08uiy4XkTSTmZE53H5hyawlE97e3enq76K9p52OniBIDAWKjp6OEfMjXhtjncu9l+kf7J/6gxylZ6CHnoEeWmmNe5sHVj+g4DDDKDhMkJmtBxYSvZsixnbgy2Z2P7AT+FvggBpGishkmBkF2QUUZBcwr2jede+vb6CPSF+ESG9kxGNnX+frlkV6o8tjl42xvKu/i66+Lrr7u+kZ6JlUXflZ+dd9bJJcCg4T9yHgSXcfcQnC3S9EQ8PDwPcI+nF4fwj1iYi8TnZmNrMygzMgiTDog8NtMrr6uujq7xpzvru/e0TgmIpQJMml4DBB7v7gVV7bBdQksRwRkWkhwzKGz5CgkwgpLSPsAkRERGTmUHAQERGRuCk4iIiISNwUHERERCRuCg4iIiISNwUHERERiZuCg4iIiMTNrgzuKGEwswtA/SQ3LwcuTmE5M4GOOT3omNPD9RxzlbtXTGUxEh8FhxnMzF5297Vh15FMOub0oGNOD+l4zKlAlypEREQkbgoOIiIiEjcFh5nt0bALCIGOOT3omNNDOh7zjKc2DiIiIhI3nXEQERGRuCk4iIiISNwUHERERCRuCg4pwMy+Z2bnzKzdzI6a2cfCrimRzCzXzL5lZvVm1mFm+8xsS9h1JZqZfcrMXjazHjP7z7DrSQQzKzOz7WYWiX6+fxx2TYmWDp9rrHT9/qaSrLALkCnxReCj7t5jZjXAHjPb5+7/G3ZhCZIFNAAbgNPA3cB/m9kt7n4qzMIS7CzwD8BdQH7ItSTK14FeoBK4DdhpZvvd/VC4ZSVUOnyusdL1+5sydMYhBbj7IXfvGXoanZaFWFJCuXvE3R9y91PuPujuPwNOAmvCri2R3P1Jd/8x0Bx2LYlgZoXA/cDfuPtld38e2AE8EG5liZXqn+to6fr9TSUKDinCzB4xs07gMHAOeCrkkpLGzCqBFUAq/ypNByuAAXc/GrNsP3BTSPVIEuj7O/MoOKQId/8zoBh4G/Ak0HP1LVKDmWUD3we+4+6Hw65HrksR0DZqWRvBv2tJQfr+zkwKDtOcme0xMx9nej52XXcfiJ7evQH403Aqvn7xHrOZZQD/RXBN/FOhFTwFJvI5p7DLQMmoZSVARwi1SIKl0vc33ahx5DTn7hsnsVkWM7iNQzzHbGYGfIugEd3d7t6X6LoSaZKfc6o5CmSZ2XJ3PxZdthqdwk45qfb9TTc64zDDmdlcM3u/mRWZWaaZ3QV8APhl2LUl2L8Bq4B73L0r7GKSwcyyzCwPyAQyzSzPzFIm/Lt7hOAy29+bWaGZ/QFwL8Gv0pSV6p/rONLu+5tKNFbFDGdmFcA2gl9mGUA98DV3/0aohSWQmVUBpwjacfTHvPSgu38/lKKSwMweAv5u1OLPuftDya8mMcysDPg2cCfBXQZ/6e6PhVtVYqXD5xorXb+/qUTBQUREROKmSxUiIiISNwUHERERiZuCg4iIiMRNwUFERETipuAgIiIicVNwEBERkbgpOIhMc2b2WTP7ZhLf79dm9sbo/ENm9r1kvfdVavoXM/tE2HWIiLqcFgmdmV2OeVpA0DHOQPT5g+7+hSTWcg/Q4e77kvWecfoy8JKZfdvde8MuRiSd6YyDSMjcvWhoAk4TdMM7tCzZPel9ghC6eL5WF8vufo5gyPityalIRMaj4CAyzcVeLjCzJdERM//EzBrMrMXMPmFmbzKzA2bWamYPj9r+I2ZWG133f6Jd/o71PjnA24G9o17KMbPvmlmHmR0ys7Ux26yKjuzZGn1ta8xre8zsYzHPPzxqdFM3s0+a2THgmAW+YmZNZtYWPZ6bY+rYA7xzwv8BRWRKKTiIzExvAZYD7wO+Cvw1sBm4CfgjM9sAYGb3AZ8F3gNUAM8Bj4+zz+XAoLu/Nmr5VuAHwCxgB/BwdN/ZwE+BZ4C5wKeB75vZygkcx33RY3kD8A7gDmBF9L3eRzBexZBagjFZRCRECg4iM9Pn3b3b3Z8BIsDj7t7k7mcIwsEbo+s9CHzR3WvdvR/4AnDbOGcdZgEdYyx/3t2fcvcBgssYQ3+81wFFwJfcvdfdfwn8jGB01nh90d0vRUdI7AOKgRqCcXRqo5cohnREaxSRECk4iMxMjTHzXWM8L4rOVwH/Gr2U0ApcAgxYOMY+Wwj+cI92Pma+Exga9nkB0ODugzGv14+z7/E0DM1Eg8fDwNeBRjN71MxKYtYtBlonsG8RSQAFB5HU1kBwZ8asmCnf3V8YY91jgJlZvH/4zwKLzCz2/yOLgTPR+QjBXSJD5o2xjxHD87r719x9DcEllxXAn8e8vArYH2dtIpIgCg4iqe3fgb8ys5sAzKzUzP5wrBXdvQ/YBWyIc9+/JQgHf2Fm2Wa2EbiHoD0EwCvAe8yswMyqgY9ebWfRBp5vibadiADdXLktlWhdT8dZm4gkiIKDSApz9+3APwI/MLN24CCw5Sqb/AfwQJz77iVoOLkFuAg8AnzQ3Q9HV/kK0EtwGeU7wLVuLS0BvkFwyaSeoGHkPwGY2XyCBpQ/jqc2EUkcc/drryUiaSN6y+Snp1MnUGb2z8AJd38k7FpE0p2Cg4iIiMRNlypEREQkbgoOIiIiEjcFBxEREYmbgoOIiIjETcFBRERE4qbgICIiInFTcBAREZG4/T/N+jdoJRVmVAAAAABJRU5ErkJggg==\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"#revised code calls variable ambient_temp rather than variable T_ambient\n", | |
"import numpy as np\n", | |
"step = 0.01\n", | |
"tmax = 2 #adjusting this variable will show temp. of corpse after two hours \n", | |
"tmin =-3\n", | |
"current_temp = 85 # at t = 0\n", | |
"analytical_array = np.zeros(int((tmax-tmin)/step))\n", | |
"numerical_array = np.zeros(int((tmax-tmin)/step))\n", | |
"xaxis = np.zeros(int((tmax-tmin)/step))\n", | |
"xaxis_em = np.zeros(int((tmax-tmin)/step))\n", | |
"j=0\n", | |
"\n", | |
"def analytical_sol(t): #function allows you to get population at any year\n", | |
" #print(\"t = \", t)\n", | |
" current_temp = ambient_temp(t) + (85 - ambient_temp(t))* np.exp(-0.275*(t))\n", | |
" return current_temp\n", | |
"\n", | |
"for i in np.arange(tmin, tmax, step):\n", | |
" # x axis\n", | |
" xaxis[j] = i\n", | |
" #Analytical Estimate\n", | |
" analytical_array[j] = analytical_sol(i)\n", | |
" #Estimate based on Euler's Method\n", | |
" slope = -0.275 * (current_temp-ambient_temp(i))\n", | |
" # print(\"i = \", i)\n", | |
"# print(\"slope = \", slope)\n", | |
"# print(\"ambient_temp = \", ambient_temp(i))\n", | |
" #print(\"slope = \", slope)\n", | |
"# print(\"analytical_sol = \",analytical_array[j])\n", | |
" numerical_array[j] = slope * step + current_temp\n", | |
" current_temp = numerical_array[j]\n", | |
" xaxis_em[j] = j*step + tmin\n", | |
" j = j + 1\n", | |
" \n", | |
"#print(\"xaxis = \",xaxis) \n", | |
"#print(\"xaxis_em = \",xaxis_em) \n", | |
"#print(\"analytical_sol = \",analytical_array) \n", | |
"\n", | |
"\n", | |
"import matplotlib.pyplot as plt\n", | |
"%matplotlib inline\n", | |
"\n", | |
"plt.rcParams.update({'font.size': 12}) #font size\n", | |
"plt.rcParams['lines.linewidth'] = 3 #line width\n", | |
"\n", | |
"#Plot \n", | |
"plt.plot(xaxis, analytical_array, color='red', linestyle='-', label='Analytical Estimate')\n", | |
"plt.plot(xaxis_em, numerical_array, color='green', linestyle='-', label='Numerical Estimate')\n", | |
"\n", | |
"plt.xlabel('Time (hours)')\n", | |
"plt.ylabel('Temperature (deg. C)')\n", | |
"plt.title('Comparison Between Analytical and Numerical Temperature Estimates')\n", | |
"plt.legend(loc='best'); \n", | |
"\n", | |
"'''Using the updated model, the new time of death is at t = -1.53. This was found by \n", | |
"printing both analytical_array[j] and looking for t at 98.6 deg. C, which is represented by i\"" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "Python 3", | |
"language": "python", | |
"name": "python3" | |
}, | |
"language_info": { | |
"codemirror_mode": { | |
"name": "ipython", | |
"version": 3 | |
}, | |
"file_extension": ".py", | |
"mimetype": "text/x-python", | |
"name": "python", | |
"nbconvert_exporter": "python", | |
"pygments_lexer": "ipython3", | |
"version": "3.7.3" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 4 | |
} |