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.
630 lines (630 sloc)
63.8 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": 18, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"approximation of DT/dt is -5.5\n", | |
"calculated value of K is 0.275\n" | |
] | |
} | |
], | |
"source": [ | |
"import numpy as np\n", | |
"import math\n", | |
"\n", | |
"Temp_ambient=65\n", | |
"Temp_t1=85 #temp 1\n", | |
"t1=11 #time 1\n", | |
"Temp_t2=74 #temp 2\n", | |
"t2=13 #time 2\n", | |
"delta_t=t2-t1\n", | |
"dTdt = (Temp_t2-Temp_t1)/delta_t #approximation of DT/dt\n", | |
"\n", | |
"print('approximation of DT/dt is',dTdt)\n", | |
"\n", | |
"k=-dTdt/(Temp_t1-Temp_ambient)\n", | |
"print('calculated value of K is',k)" | |
] | |
}, | |
{ | |
"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": 10, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def measure_K(Temp_t1,Temp_t2,Temp_ambient,delta_t):\n", | |
" ''' Determine the value of K based upon temperature of corpse \n", | |
" when discovered, Temp_t1\n", | |
" after time, delta_t, Temp_t2\n", | |
" with ambient temperature, Temp_ambient\n", | |
" Arguments\n", | |
" ---------\n", | |
" your inputs...\n", | |
" \n", | |
" Returns\n", | |
" -------\n", | |
" your outputs...\n", | |
" \n", | |
" '''\n", | |
" #k=-dTdt/(T1-Tamb)\n", | |
" #dTdt=(T2-T1)/(t2-t1)\n", | |
" #k=-(T2-T1)/((t2-t1)*(T1-Tamb))\n", | |
" k=-(Temp_t2-Temp_t1)/(delta_t*(Temp_t1-Temp_ambient))\n", | |
" return k\n", | |
" " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 11, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"0.275\n" | |
] | |
} | |
], | |
"source": [ | |
"print(measure_K(85,74,65,2))" | |
] | |
}, | |
{ | |
"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": 219, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def temperature(N):\n", | |
" ''' \n", | |
" computes the analytical and numerical solution for the temperature of a corpse\n", | |
" N is number of timesteps\n", | |
" temp_analytical is the 64-bit floating point \"true\" solution\n", | |
" temp_numerical is the 32-bit approximation of the velocity'''\n", | |
" \n", | |
" T0=85\n", | |
" Ta=65\n", | |
" t=np.linspace(0,N,N+1)\n", | |
" temp_analytical = Ta+((T0-Ta)*np.exp(-k*t))\n", | |
" temp_numerical=np.empty(len(t),dtype=np.float32)\n", | |
" delta_time =np.diff(t)\n", | |
" temp_numerical[0]=85\n", | |
" for i in range(0,N-1):\n", | |
" temp_numerical[i+1]=temp_numerical[i]+(k*(Ta-temp_numerical[i]))*delta_time[i]\n", | |
" \n", | |
" return temp_analytical,temp_numerical,t\n", | |
" " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 136, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"(array([85. , 80.19144246, 76.53899621, 73.76469985, 71.65742167,\n", | |
" 70.05679192, 68.84099817, 67.91751514, 67.21606317, 66.68325981,\n", | |
" 66.27855722, 65.97115643, 65.73766335, 65.56030852, 65.42559473,\n", | |
" 65.32326989, 65.2455468 , 65.1865105 , 65.14166818, 65.1076072 ,\n", | |
" 65.08173543]),\n", | |
" array([8.5000000e+01, 7.9500000e+01, 7.5512497e+01, 7.2621559e+01,\n", | |
" 7.0525627e+01, 6.9006081e+01, 6.7904411e+01, 6.7105698e+01,\n", | |
" 6.6526634e+01, 6.6106812e+01, 6.5802437e+01, 6.5581764e+01,\n", | |
" 6.5421776e+01, 6.5305786e+01, 6.5221695e+01, 6.5160728e+01,\n", | |
" 6.5116531e+01, 6.5084488e+01, 6.5061256e+01, 6.5044411e+01,\n", | |
" 1.3452465e-43], dtype=float32),\n", | |
" array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,\n", | |
" 13., 14., 15., 16., 17., 18., 19., 20.]))" | |
] | |
}, | |
"execution_count": 136, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"temperature(20)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 137, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"(array([85. , 80.19144246, 76.53899621, 73.76469985, 71.65742167,\n", | |
" 70.05679192, 68.84099817, 67.91751514, 67.21606317, 66.68325981,\n", | |
" 66.27855722, 65.97115643, 65.73766335, 65.56030852, 65.42559473,\n", | |
" 65.32326989]),\n", | |
" array([85. , 79.5 , 75.5125 , 72.62156 , 70.52563 , 69.00608 ,\n", | |
" 67.90441 , 67.1057 , 66.526634, 66.10681 , 65.80244 , 65.581764,\n", | |
" 65.421776, 65.305786, 65.221695, 2.8125 ], dtype=float32),\n", | |
" array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,\n", | |
" 13., 14., 15.]))" | |
] | |
}, | |
"execution_count": 137, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"temperature(15)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 138, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"(array([85. , 80.19144246, 76.53899621, 73.76469985, 71.65742167,\n", | |
" 70.05679192, 68.84099817, 67.91751514, 67.21606317, 66.68325981,\n", | |
" 66.27855722]),\n", | |
" array([85. , 79.5 , 75.5125 , 72.62156 , 70.52563 , 69.00608 ,\n", | |
" 67.90441 , 67.1057 , 66.526634, 66.10681 , 0. ],\n", | |
" dtype=float32),\n", | |
" array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]))" | |
] | |
}, | |
"execution_count": 138, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"temperature(10)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 133, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"import matplotlib.pyplot as plt\n", | |
"\n", | |
"plt.rcParams.update({'font.size': 22})\n", | |
"plt.rcParams['lines.linewidth'] = 3" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 154, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"n=16\n", | |
"\n", | |
"temp_analytical,temp_numerical,t=temperature(n);\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 155, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"<matplotlib.legend.Legend at 0x7f4aa24ea5c0>" | |
] | |
}, | |
"execution_count": 155, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAE0CAYAAACfP2y8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeXwU9f348dc7JwFCwi3hSFAQBZVTUSnlUEEFK4eotSKI1Ape3/LDKt9qvb5VKlbUFq1aNSj2ULRYRaVeoOAJggcoigWVQ+4QjgA53r8/PrPJJtndZDebbDa8n4/HPGZ35vOZ+cy1752Zz3xGVBVjjDEmHiTEugDGGGNMdVnQMsYYEzcsaBljjIkbFrSMMcbEDQtaxhhj4oYFLWOMMXEjboOWiEwUEQ3QFYtInoisFJH7ReS4WJc13ojzSxF5z1uXJd66XRDGNHzbY2ItFjVuiMgGb33cFuuyxAsROUZEcr11d8hvn+oV67KZ2InboBVCApAB9AKuBz4VkV/W5gwb4A/SH4BHgdNw61JiW5zY8H4wVUQWx7osRxoR6Qh8BEwAsoGU2JbI1BcNJWidC6R7XQvgZOBeoAi3s/9FRE6JXfHih4g0xQV7gBeA43CBKx24OFblMkecKbhjuQAYDbSl7Bj/LIblMjGWFOsCREmBqu7z+74cWC4im4H7cMH518DPY1G4OHM8Zf9q71TVtbEsjDli9fT6/1HVal+WNg1fQznTCubPwH7v809iWZA40tjvc17MSmGOdL790PZBU06DDlqqWgh8631tEyydiGSJyK9E5EXv/tRBETkgIv8VkadEpH+QfLkiorhr7gC3BqgYMjFAviQRmSQir4nIjyJyWES2i8h/ROQXIlLje0giMlBE/iYi33vLs1tEPhSRm7xLgBXT3+Yty2K/wev9l6WG5ZkgIku9chwQkU9F5DciUuW9ChHpJyKPi8g6EdkvIntF5DMRuUtEWoXIly4iY73t9KWX97CIbBGRl0XkwkDr2lfJB3c/BWBQgO2aG/HKCF7e1t4yfSoi+SJSICLfishjItKjirzZIjLbWy97veXcLCKrRORRERkdJN/pIjLP29cLvG3znYi8LyIzReTkGixPuojM8Pa73d5++L2IPCMiA4Lk8e1rg71BEyqs99vCLMMJIvK/IrJYRLaKSKG3bj8TkT96985qTEQyRORGEXlHRLaJqziy0dvn/1dEjg6SL0FExnu/Bdu87bbN+x7yt8D/N0ZEEkVkqogsE5Gd3vD/8dLl+KUdLCJNRORWEVntHRO7ROR1ERlZxTJGtI95eaP3m6eqcdkBEwH1usEh0n3upfkxRJrdftMK1JUAMwLky60inwITK+TpBKyqIs8rQJMI14sAs6uY/ndA9wr5bqtqWcIshy/fFcBzIaa7CEgIMo0E4I/e+g+WfzvQP0j+f1Vj+7wIpITYt4J1uWGujw1evtuCjB+CO6sINr8i4H+C5B0E7KuivPsC5JtejeV8OcL98ERgYxXTngVIkP0mWBdw/QUpQ89qTC8fOKuGv0VDvf0w1HwWBMiXCSypIt+bQLMqjrFfAW8HyPs/Xrocv2GjKftNDNTdG819rDZ+8yLeULHuqEbQApJxlwcVeCvEtJZ7B9AwoAfQytvQw4D5fvM5u0K+VKApLggocJf33b9L8kvfDFjnpd0BTMNVdMgEugI3AQe88fMiXC//61fed4EzvOU5BpjhN/0fgOZ++VK88p7jl7+7/7KEWQ7fNL7F/eDOBE4AmuN+TF7wS3NlkGnM8kszFxjoLUtbYJTfwbcNaBcg/1+BJ4GxQG+gHZAF9Mfd6/Sti7sr5Evylnme33qsuF1Tw1wfGwjyo+tte98Pwk5cJYROuKsD5wFf+K2HcRXyJvjtf+uA8d62bgEc623/mcDqCvmO9baLAm8Aw715tsQFnJ956+6ZCPbBlsBmb9oHvP2ui7fthgJL/Zbnhgp5fev3XW/8vArrPSWMcpwEfAD8xtt3unrrpRtwCfCxN4/dwFERHm+nAIe86ezCHcM9vP28IzASeKziesT9uVzktx7+6u2jLbz+437jAv5x8Bu/ESjGHS8n+k2jl5cuxy/teq+8v/PbJmfhfgN9aa6Ixj5WW795YW+k+tJRvaA1zS/N2BrM6w/eNJYEGb+BavwLBB6k7N9dtyBpzvQrc78wy9nW7wBaHOgAx9W09E3/vgDjB/uNz6nBOvP/F/WLAOMTgJXe+A8CjO9L2RlWpbNcL01TYI2X5s8RlNEXoPcC6QHG5/rWZRT216D7CLDAG1cA9AwwPhNY66XZ4r9dcT9SvvVcKW+I8lzr5dkaaD+p4bLe7027BBgeYHwKZUHpINAmQJrF3vjcaJatwjySKAugt0eQX/z2v+1A11DzqvB9lN92uytInnv80vwswHj/Y2xKiHnnVEh7SYA0TSj7E7gdvz9lke5jXt6o/+bVys5QFx3lg9Y5lP0Ta+794N0DFBLkxznMeXX3plMENA4wPugPUoWdwnfWN72K+b3ppZsdZjn9L/ecFCLdi5T9w0ysMG6w3zRyarDOfNN4L0SaX3tpDgU4qOd64z6jwiWkCunGU3aGEjRdiPzbvPzDAozLpZaDFu5sqriq/RR3Wce3Tsf4De/tN7xVGOXxrfsvarpsFaabSNllzhdCpPMv97QA4xdTy0HLm8/UqvbTEHmH+y3D5WHmfZkAf0IqpEnF/alQ4N8BxvvmvaaKeeX4pa30B9Ev3Qi/dBdEYR+rld+8hlIR4xXcv+W9uFP05cANuB/Dc1R1WlUTEHez/y/eTcY94lrW8N0UXu0lS8SdFkfidMpqRC0RkabBOsqeQ+kX5jx8NSS/VtVQz7I85/UzcZfsatOrIcb5qtOn4P5s+DvT678NNAmxrtZ46VoAlW52i0gHEblTXOseO70b8f6VS1p7SbtFtHQ1dzplFaKeC5HuZdyZGLhLXT5r/YbPFZFjqznflV6/h1fhomU181XlRNxzfRBieVR1JWWVpAYGSxcNIjJaRJ4VV6llf4XtP8dLFsn2P8PrHwL+HkZ5BPBVRPm3qh4OlE5VDwEveV9D1X5+pbrzxt3nDeY13KU6/MoHke9jtfKb11Ce0wqmCTBbRFap6o/BEonI/+HuBVWnBktG1UkC8j8oPqpmntZVJykn2+uvCZmqLAj78nwa5nzCsTnEuAN+n0ur2ns7cZb39Tqvq47WlP0QIiLnU3ZPpCqRbteayvb7HHS7qWqhiHyDu0+T7Tf8gIjMwF2SOxc4V0TW4i6/vQO8oapbAkxvsbhmuUYBNwLTRWQ57nLZO8Cbqrq/Yr5oLY9nNe5PYHYV6SIiIo1x906HVyN5JNvf9wd2raoeDCNfM9wfRqj+sdpcRJqpan6ANP8NY95fBRuhqsXePtaTKOxj1NJvXkM50xqiqqKqgtshTsXdJwB30+9vwapUishFwG9xAetd4Be4G6mtvWml4/49+kQa6CM5KBqFmT7d6+8LmcqdkVbMU1uKq5nOf/tEGkBK15eI5AD/wAWsDbjA1w9XGcPXwkc67iY2xO4PnP/6r+52K7fNVPUB3OXDD3CXWLoBk4GngI0i8lKQf8cX4q5I/Bd3FaE/8P9wl4+3icifRKRZeItT8+WJoj9SFrDm4ipF+Cof+Lb/VG98YgTT962bvSFTVRbJOqqYz9+BIMMDqWp+vvHR2Mdq5TevwZ1pqepe4EMRGYOr+TcGV514PG4FV3S1138PV6GjpGICEUmOQtH8d5Z0Ld+CR7T4dvCqziz8x4d7wNUF/3VzjarOCZoyuEm4AyAfOFVVtwZKFMGPcrT5r/+mwJ4QaX3brdI2U9dqxAIRaY27LPMT3L/i7rgf6wEi0kdVN/jlKcQ1d3aviHTBtTX5Uy/9UcA1wKkicpqqFkW4PKEEXZ6a8s6yJnpf/6CqNwVJF+4fQ3+RBt1I1lHFfJGKeJtEsI/Vym9eQznTqkTd3b2rcD9aALcHCT6+FqOfDRSwPCcGGR4O/1P42mqleoPX715FOv8HVTcESxQrqroHV7EC3E3gSPjW8VshAlYnyv4tx8oGv89Bt5uIJOGqF1fMU46qblfVF1X1BlXtgWu6rAR3z/B/QuRbp6pPq+ovcVW1H/RG9cP9IFWXf9mqux9uCJUoQsdR9q891P2mmhzb67x+tzCDXz6uEhRUfx3tCnJpMFxB33ohIom4aujgqrgHFMY+Viu/eQ02aIFbubhLBOBq0IwPkCzV64e6PBAon7/CakxjCe6GLcDlVUwvUku9/rEiEupgvMDr5+GeAaqPXvf6YyI8G6qr7VpT7+MOeHDPkwUzAkjzPi8Nka4cVf0HZdv4+GrmKcI9bO5TrXyeLyg7Wwy6PCLSE3epDsJYnjCk+n0OuP28e6ejajAP3z6aShiNSXt/qJd5X88LdiVHXGsx53lflwVKE4GgrVYAZ1N2b7na8wuxj9XKb16DDlqe+yn7V3OT92/Cn+/fwM8C3fcSkQmU1WQLxndWkBUsgfcv6XHv60QRCfUDhYg0E5F2Vcy3onmArybSA96/84rTPZuyA/WJEGeXsTbb6zcHHq/qEm2A6+m+7TogUM04ETkB92BjKFVu15pS1W2U1RCb6pWrHBHJwD0rCPAjriahb1x7CdAsl9/4NMrKv9NveFcRCXX8+9eS3Rk0VQWqWox7VADcH45Kx463LX1ncgeBp6s7/TCs9/t8foAyiFeGirVWw/EGZRUp7hGRoDWLAxyLvt+CdrgHfQO5HffsJbgHlKOhv4hcUnGgiDTBPSAM7iHgl/zGRbSP1dpvXjjPFtSnjmo24+SlvcUv7SUVxt3mN+7vuGe8WuJqad2HezZrdah54arN+p4VGoa7AZnkdeKXLgNXe0dx/64fxzWP0hZ38HTF/Tv9K+7f6gURrBf/FjEW4+7ntcRVB7+RsucmNuLXIoZf/sF++XNqsH1805gYIk3IeeEOIt/4T3CVZDrjal619/Lf4m2fl0NMeznuqf82Xv7rcY9GbPG2WcBn7IBxftO42svv264Bm54KsawbQszHv0WM7bhmeTriKgONpHyzOxcGOA7yvH1mDO7sxdcawwjKtz5xnl++XNwP+93euunorddjcP+Kv/fy7CNAayNVLKt/ixj7cC1SHO0NH4KrceYr02+CTGMxNXxOy28ah4FbcZUHWuLu2/mekyo9tiOch3+LGDu9ZT3eW5cdcM+QPkzgFjFe81sPj+AuobXA1d57zG/cwkiPMS9djl/a9V55b/G2dUtv+3/sl6ZiixgR7WNe3qj/5kW0M9SHjvCCVjPcj5TiTmP9A0ljXHVMDdJ94e2YoYLWCX47bsVuYoW07Sps5FBdpafgq7FehLIWCYJ1ldoe9Ms/2C9dTg22TzSCluD+VBRVsTwKPB8g/5wQ6XfhbiRvIHgwScNVoQ+UPzfM9RF0Pt746rQ9eH0Vx0GormJTVbnVyHMAGBXh9q9O24P3EuSBcKITtI7DnTUEm/8/cBV2lAiDljefMyn7fQnWRdr24FtU3fZg0GPMS5fjl3Y05ZsFq9j9MVr7mF/+qP7mRbSR6kNHGEHLS3+zX/oxFcY1xv04fom7XJGH+2f/v964nKrmhasu/DywCffPLugOhfsxHgU8iwsgBbigt9nbSW8EutRw/QzEnTn+4E07D/gQd0ksaDuC1KOg5ZeuK+6sd5W3HEVefxUuMJ1BhRY1/PJOwN032of7Ef4Gd1ko2xu/gdDBJAv3T/kbbzvVStDy0rTBtV/5Ka72VgEuaP4VOCFInsa4f7uzveX83tuHD+AeCn0SV3uyYr62wKXe+JW4y46FuEoCn+BalOlUw30w3TuGPvK21yGvfH8DBlSRd3Ek6znAdLJxZy0bveNyG671hUu98RN927SG82mOO3v5EBfAfMv6rnfMZQfJl4C7t7rIK1uh11/kbZ9QrcFEErQG42oI3om7tLkfd/vkdWBkNPexCtOI2m+eeBM0xhjTAHnPLK73vg5R1cUxK0wUHAkVMYwxxjQQFrSMMcbEDQtaxhhj4oYFLWOMMXHDKmJ4WrVqpTk5ObEuhjHGRNWhQ4f44gvXYMWxxx5Lenp02ydesWLFDlUN940UEWtwDeZGKicnh+XLl8e6GMYYE1dEJGg7hbXBLg8aY4yJGxa0jDHGxA0LWsYYY+KGBS1jjDFxw4KWMcaYuGG1B2tgwcpNzFq0ls15BWRlpnHD8G6M6t0+1sUyxpgGy860IrRg5SZmvPA5XfI/4IrEhWzL28uMFz5nwcpNsS6aMcY0WBa0IjRr0VpKCgu4PSmXm5Of4bWUG+lXvJJZi9bGumjGGNNg2eXBCG3OK+DSxMXkJGwF4JiELTydMpP/7O8LuzpDi84xLqExxjQ8dqYVoazMNP5RPJQ7C3/BXk0rHT4scQXM6Q9v/R8cPhDDEhpjTMNjQStCNwzvRlJyKo8Xj2DooT/yXNFPy0YWH4J3ZsGfT4YvXgBr39EYY6LCglaERvVuz91jTqR9Zho7yOT+pr9mycC/Q1afskT5G2H+5TD3PNi6OnaFNcaYBsJaeff069dPo9JgbkkJrHoG3rgNDuwoGy6JcPJkGDID0prXfD7GGFMPiMgKVe1XV/OzM61oS0iAPuPh2hVw6lQXrAC0GD56BP7UF1bkQklxTItpjDHxyIJWbUnLhLPvhinLoLPf/a4DO+Gl6+GxofD9h7ErnzHGxCELWrWtzfFw2b/hwqcgo2PZ8C2r4Ilh8MKvYO+PsSufMcbEEQtadUEEup8PV38Eg26CpEZl4z77h7tkuOwBKDocuzIaY0wcsKBVl1Iau4oYV38Ex59XNvzwPnj9d/DwafDNG7ErnzHG1HMWtGKheTZcNA/GL4BW3cqG71wHz4zlk9tP54NX57maiMYYY0pZ0IqlY4bAlGV8fsJN5VrV6KOrOfXDq9n7x97w8ePWsoYxxngsaMVaYjJXrevPkEP38Y+iwRRp2SZJ378BFk6D2T1cs1B7t8aunMYYUw9Y0KoHNucVsIMMbiq6kp8eup9HikaQ73fmRcEu1yzU/SfAgqnWuoYx5ohlQaseyMosC1CbacXdRb/gtEN/5oGkSZDZqSxh8WHX2sbDp8NTo1ylDWvRxBhzBLGgVQ/cMLwbacmJ5YaVJDcle8R0uHYljJsLHU4un+m/b8MzY+GhU+GTp6DwYB2W2BhjYqNeBC0RGSwiWs2uk1++3CrSfhXL5aou/8Z3BWifmcbdY05kVO/2kJgEPUbB5Dfgitfd817it9m2fwX/vtZdOlz8B9i/I+h8jDEm3tWLBnNF5DjgphBJTgGOB74FuqpXaBHJBSYAy4B1AfJtUdUZ1SlD1BrMrQu7N8AHf4GVT7tnvPwlNYKeF7t2D1t3C5jdGGOipa4bzK0XQasqIrIa6A78VlXv8hueiwtal6tqbk3mEVdBy6cgz10a/PAvkL+p8viuw6DPZdDlLEhuVHm8McbUUF0HraS6mlGkROQ0XMAqBubGuDj1S1omDLgOTp0Ca16E9/7k2jT0+eY/8M1/KExqQnL3kdBjDBwzFJJSYldmY4ypgXpxT6sKk7z+a6oa4HTCkJgMJ14AVy7m3Z88xRt6MiUqpaOTi/bDZ/+Ev18E93aBBVfDujeguDB2ZTbGmAjU6zMtEWkMXOR9fTxE0iEichLQFNgKLAVeV9Ujqx0kEW5ans6mQ78mR7YwNvFdRia8T+cEv4eSD+6BVfNcl9bCtYF4whjI/omr9GGMMfVYvb6nJSITgFxgG9BBVQsrjM/F3dMKZA1wsap+Xp15xeU9rQA637SQ8ltU6SEbOC/xA65q9SnkfR84Y5PWrmZij9HQ6TRISAyczhhj/Nibi8vzXRp8qmLA8qwCrgN64M6ysoCRwKe4+2BviEj7YBMXkStFZLmILN++fXt0Sx4j/g8qO8Jq7czTTSfB9Z/B5LfgtGugWYXVsn87fPxXyB0B93WHV290L6m0RnuNMfVIvT3TEpEuwDfe1+6q+mUYeVOAJcCpwBxVvaaqPA3lTGvByk3MeOFzCgqLS4elJSeWPfflU1ICGz+CL16ANQtgX5B2DZt1cM+J9RgDWb0hob7/zzHG1CWr8u4RkbuAGcD7qnp6BPl/BrwIrFfVo6tK31CCFrjANWvRWjbnFZCVmcYNw7uVD1gVlRTDd+/B6n+5WogHgjyg3LgVdB4InX8KnQdBi6PdCy6NMUcsC1qAiCQC3wHtgcmqGqoSRrBpHAusBQ6rampV6RtS0KqR4iLY8C6sfgG+fAkKdgdP26y9F8C8LqND3ZXTGFMvWNACRORcYCGwHzhKVfdVkSXQNE4D3gN2qWrLqtJb0KrsxRUbWLjg75yl73FG4ie0kCo2Q4ujywJYzkBo2qZuCmqMiRl7uNi5wuv/M5KA5bnQ638chfIcke55/b9sKjyJ/3ASUlTCcfIDpyesZkjqV/wk+Ws4lF8+w67/um5FrvveprsLXp1/CjkDIK15nS+DMaZhqXdBS0Ra4WoAQohns0SkF9ABeFVVi/2GJ+FqFF7nDZpdS0Vt8DbnFZR+VhL4UrP5sjibJw6cy/rfD4ctn8L6JbD+Hfj+AygqKD+BbWtc99EjgEC7nmX3wzr0tSBmjAlbvQtawHggBfhKVd8LkS4H+BewS0S+BjYC6cCJuKrvJcCNqrqodovbcGVlprEpryDgcBKTXODp0BcGToOiQ7BxuQtg69+BjR9Dif9TCuqamNqyCt570A3KzHaBrF1PaNfL9Zu2rpuFM8bEpXp3T0tEPsMFnt+o6qwQ6ToD1+NagM8GWgKKC17v4qq6r6jufO2eVmXVrj4fyOH97uxr/TuuYsfmlVCdBkrSs/wCmdc1y7JaisbUU1YRI0YsaAUWdvX5YA7u4YO3X2Lt+wvpzZd0kx9IlaLq5W3SunIgy8y2QGZMPWBBK0YsaNW+ATPfKr3cmEwRXWUjPRI2cGqjHxjbbif8+Hnl+2LBNMooC2Ctj4OWXVzXuKUFM2PqkNUeNA2Wf8WOQpJYozmsKc5h/n4YO3mEe0Zs5zeugkdp9xkc3lt5Ygf3lN0/89cooyyAtewCLY9x/RbHQGrTWl5CY0xts6Bl6kzIih3gKne0Od51PS92w0pKYPd6rxKHF8g2r4KDeYFncnAPbFrhuorS21UIaF7XPNu93sUYU+9Z0DJ15obh3QJW7LhheLfgmRISvLOlY+CEsW6YKmPv/iet9n1J94Tv6SxbvO5HmsrB4NPau8V1G94tP1wSoXmOC2CZHV1LH83aQ4bXb5YFSVU2qmKMqQMWtEyd8VXgqHHFDhE+yU9HOYVFJaf4jVDakMdHV+XAznVe9y3s+MadrZUEqfihxbDrW9cF06RNWRDL6FAWzHyf09vZ+8iMqQN2lJk6Nap3+8hqH1YQ+FKjkJyZBTk/cZ2/4iLI+84FsdKA5gW1/I1Vz3D/NtdtXhl4vCRA06PKB7Ymrf26VmWfkxtFtMzGGAtaJk6FfakxMansMiPDyo0aeverpORvIEd+JEt20s7rcpJ2c0LTfbDvx6qfMdMS2LvZdVW1HJaSXj6IBf3cGhq3sBdyGuPHgpaJS1G71Ais31OC0omvtFO54VII62/1ajXu+xH2bHJnZXs2Qf4m2LPR629yZ2HVdXiv63avrzKpIuSRzu6SJhQkptOqdRvatmkLjTIhLdPVliz97H33fU5tZu8/Mw2OBS0Tt2r3UmOFWo0ZHbxXr/QPPJGiw7B3M1MffonkfZtpJ7toKfm0lHxasYejkvZxbNOD7g3Rwe6tBSAozcmneUI+6BbY9jVUOz4KNGpWGsy2F6Xxxa4EdhxOhtR0Tjo6i26d2rkzv5Qm7pGAFK+r+DkxxZ5/M/WCBS1zxIuoVmNFSSnQPIdX849GqfzO0dKzNlVXXX//DhfASrsdlT7v2bGZDCJ9yQGAukcADu4BoDUwBNxRX4x7L/g3QTOXl5DkBTAX4HYVpfDtHiGvMAmS0zi2Q2uy27Zy9+uSG0NyGiSluX5yY294oGGNeXlNHve8+R0/7Dlcs1ZXzBHBgpY54kXzUmOVZ20irnX7tObQqmvIafW6aSGJFNGcvWTIfppxgAzZTyb7mX1+jgtGBXkuCAb6HOih7EiVFHnTds/HtfA6EnFNU3/vdREY6XVFqQkcLkjm8IIkCl5rTFpaGiSmQlIj96cgMdU9epCU6s78yg13n7/afojF3+azo0BpkpbGkB7t6ZXd2qVPTPL6KS4IJ6a45/MSkyEhufz3xBReWb2D2W9v4Ic9hbTJaMK04d0Z1SeyF51GrTm0KE8rHlkzTh5rxslEQ40aGa7Av9krf+0z01h209CqJ1Bc5N55VrCb8+9dSDMv8DWRAppykMYcpKkc5KpT27gGjg/vg0N7/T7vc/3D+8K6pNmQKQlIYpILcglJrpJMQpIXBJPKPpeOS2ZnQQnrdx3kcEkCxSRQTCIkJNKtXQbtMpuUpZXE8v0Aw77efoC31u7kUIlQrAmUICQkJjL8hCx6tG/uarEmJLq+f1c6LNHvu5Qf1vGUiF7cas04GRPHonnWVuPLlolJrvZh4xbsyDiBT4MEwKtGVhEAVd2rZw7vh8N7OfueV72AV0Aah0nlMGlymDQOc9s5naGwoKwrKij/vcKwHXl7SOMQaRwmQer/H2ihBIoPu66aWgItBXdW6u9HrwvDscCxCUDF+jVrvK4mLnkOjh1WdboYs6BlTJRFq4JIvQmAIt79p0bQpCV7M7rxVV6BexGQn/aZadz2k2qcAfo5v/RsUkmimFQKSaGQThlJvPirvq6CS/Eh1y86WOHzYRdMiw5B8SFmLfyMFCkilUKSKCK5tCtmXO+2UFzoupJCL/B434sPuzNJv2Hb9+wliWKSKSKJYhIpIVmKq1yeuBYnNU0taBlTjzW4ABhiWkUkUUQSJclNufzsE6FFeGVb8E7wy6njxoYXTEcFvDSrdMpI5Z3pA12QKymCkmIXBMt9L/KCo/v+q6c+ZPfeAySJC36JXhBs3SSZmaO7e3mKXcss/v2SIvf8X+mwIh5Z/A37Cg6RKCUkUEICSgJKRmoCvzilg0tfmsf7rN7nkorf/dOUuFZf4oAFLWOOEPUxAMZDMC2bVhLTzu4edosm55xzVOD7nOeeCN3DWxAQR7kAACAASURBVM62TTZxf6BpjTwRjpDKGFYRw2MVMYyJf/W1ll59nVY02EsgY8SCljHGhK+ug1Z83HkzxhhjsKBljDEmjkRUEUNE2gFtgUxgN7BNVbdEs2DGGGNMRdUKWiKSBIwFzsc1X1apbqSIbAUWAwuAF1TVHqE3xhgTVSGDlog0BqYDU3CBytfMcz6wy+s3wz30fRRwMXARsE1EHgL+qKoHaqfoxhhjjjRB72mJyGXA18BtQArwKDAO6KCqmap6tKr28voZQEfgQuCvQCpwO/CNNx1jjDGmxkJVxMgF8oBLgHaqOkVVn1fVzYESq+omVZ2vqr8C2gGXevmfjHKZjTHGHKFCXR6cAMzTCB7kUtVDwN9E5O/ALyItnDHGGOMvaNBS1adrOnEv4M2r6XSMMcYYsOe0jDHGxJFQFTGGishxdVkYY4wxJpRQZ1pvADcEGiEiD4rIJbVTJGOMMSawqi4PSpDh1wBnRrksxhhjTEh2T8sYY0zcsKBljDEmbljQMsYYEzcsaBljjIkbVbXy3lJE+kQwDlX9JPJiGWOMMZVVFbRGel1FGmKcb3xE7+oyxhhjggkVWHbhgo8xxhhTL4Rqe7BVXRbEGGOMqYpVxDDGGBM3LGgZY4yJG6EazE2PxgyqOx0RyRURDdF9FSLvJSLyrojsEZF9IrJcRK4WEQvKxhjTgISqiPFfEbkbmOO91DEsIpIKXAvcCLQOI+syYF2A4VuCzGcOMBU4CLwJFAJnAH8GzhCRcapaHMb8jTHG1FOhgtaXwL3ATSLyDPCUqq6saoLes1sTgJ8DrYB3wizTX1U1tzoJRWQsLmD9CPxUVb/xhrcF3gZG4xr3fSDMMhhjjKmHQtUe/KmIXAz8AbgeuE5EdgMf4ALaTiAfaAa0BLoDpwKZuNbhNwDXqOqztVj+GV7/Rl/A8sq+VUSmAItxQfdPqlpSi+UwxhhTB0I+AKyq/xCR54ALcGc0PwHO9Tr/Z7h8rzApwZ3hPAS8WJuX5USkA9AXOAw8F6DsS0RkE9AeF0zfq62yGGOMqRtVtlrhBZ5/Av8UkebAIKAX0AbIAPKAbcAnwBJVza9hmYaIyElAU2ArsBR4PcCZUm+vv1pVC4JM62Nc0OqNBS1jjIl7YTW1pKq7gQVeV1suCzBsjYhcrKqf+w3r7PW/CzGt7yukNcYYE8fqU5XwVcB1QA/cWVYWrm3DT3H3y94QkfZ+6Zt6/f0hprnP60el+r4xxpjYqjeN2qrq/RUG7QcWisjrwBLcfakZuNqAUHYfLeL2EUXkSuBKgE6dOkU6GWOMMXWkPp1pBaSqh4G7va/n+o3a6/WbEpxv3N5AI1X1UVXtp6r9WrcO51EyY4wxsVDvg5bH1xqG/+XBDV4/O0S+jhXSGmOMiWPxErRaev19fsN8Dzr3EJG0IPlOrpDWGGNMHIuXoHWh1//YN0BVf8BVs08BxlXMICKDgA641jLer4MyGmOMqWX1ImiJSC8RGSkiiRWGJ4nINFytQoDZFbL67nX9QUS6+OVrg3vAGWCmtYZhjDENQ32pPZgD/AvYJSJfAxtx1dRPxFV9L8E11bTIP5OqzheRh4EpwOci8gZlDeY2wz1P9ue6WghjjDG1K+KgJSLNcBUjFNiiqntqUI5PcY3anoKrWNHbm+5G4ElcS/MrAmVU1akishS4GtdaRyKu4sYTwMN2lmWMMQ1H2EFLRC4Ffg30xO9ZKRH5FHhAVZ8Kd5qquh74n3Dz+eX/G/C3SPMbY4yJD9W+pyXO08Bc3JlQAq6V93zvc2/gSRF5RkQk+JSMMcaYyIRTEWMq8AtgB65iRHNVba6qzXGvI7kW13DuxbhLdcYYY0xUhRO0rgQOAYNV9c/+97BUNV9V5wBDcRUhroxuMY0xxpjwgtaxwGJV/TJYAm/c20CXYGmMMcaYSIUTtPKB3dVIl0eQtv6MMcaYmggnaC0GTheRoDUOvXGn4VplN8YYY6IqnKB1M67CxWMiUun9VCLSFHjUS3NzdIpnjDHGlAnnOa3zgPnAJOB8EXkZWO+Ny/HGZ+Ae6h1Zsda7qt5X08IaY4w5soUTtO6l7IWLmcClft/9I9SkCvnES2dByxhjTI2EE7TuowZvCTbGGGNqqtpBS1Wn12ZBjDHGmKrUi1eTGGOMMdVhQcsYY0zciKSV957AENx7rhoFSaaqen1NCmaMMcZUVO2gJSKNgHnAaN+gEMkVsKBljDEmqsI507oLGAPsA54FvvE+G2OMMXUinKB1EbAH6Kuq/62l8hhjjDFBhVMRoznwjgUsY4wxsRJO0FqPPVxsjDEmhsIJWk8Bg0WkZW0VxhhjjAklnKB1L/Ae8IaI9K+l8hhjjDFBhdOMU7GIXAwsBd4Tkb3ARqAkcHLtGaUyGmOMMUB4z2m1x73csTPuGa1mQPcgye3elzHGmKgLp8r7LOBo4CPgT8A67DktY4wxdSicoHUm8AMwRFULaqk8xhhjTFDhVMRoBHxgAcsYY0yshBO0PgesursxxpiYCSdozQYGiUif2iqMMcYYE0o497TeAv4IvCkiM4FFBK/yjqruqnnxjDHGmDLhBK3tXl9wLb7fFSKthjltY4wxpkrhBJbd2PNXxhhjYiicFjFa1WZBjDHGmKqEUxHDGGOMiSkLWsYYY+JG2EFLRPqJyJMi8pWI7BSRh/zGDRGR/xWR1tEtpjHGGBNmDT8RuQFXazDRG6S4ljJ8koE7gZ3AI9EooDHGGONT7TMtERkG/AHYAVwOHIOr/u7vDVwtw/OiVUBjjDHGJ5wzrWnAYeBsVf0UQKR8zFLVEhH5Gjg2aiWsp1SVvXv3kp+fz4EDByguLo51kYwxfpKSksjIyKBFixYkJdljow1FOFvyZOBDX8AKYSNwYuRFqv9UlW3btrF//35atGjBUUcdRWJiYqUgboyJDVXl8OHD7Ny5kx9++IHs7GwSEqzeWUMQzlZsAvxYjXSNqXzZsEHZu3cv+/fvJzs7m8zMTJKSkixgGVOPiAipqam0a9eOpKQkdu/eHesimSgJJ2j9SPUu+x2He+9Wg5Wfn0+LFi1ITEysOrExJmZEhMzMTPbv3x/ropgoCSdoLQFOEpGfBksgIqNwbzd+s6YFq88OHDhA06ZNY10MY0w1NG7cmIICew1gQxFO0LoXKAb+JSKXiEhj3wgRSRKRMcBjwCHggXAKISLJInKGiPxRRD4QkS0iclhENonIfBEZHCRfrohoiO6rcMpRXcXFxXaWZUycSEhIoKQk4MsoTBwKp+3Bz0XkKtzzV08DRbjntC4Gfg6k4F5TMllVvwmzHIOA173PPwIrgP1Ad2AsMFZE7lTV3wXJvwxYF2D4ljDLUW12D8uY+GDHasMSVj1QVX1CRFYDvwOG4ipcNMIFsMXAbar6bgTlKAGeBx6omF9ELgKeAW4RkbdV9e0A+f+qqrkRzNcYY0wcCfvhBVX9EBghIilAFq51jC2qeiDSQqjqW7iXTAYa908ROQu4ArgUCBS0jDHGHAGC3tMSkSdEZFKw8ap6WFU3qOq3NQlY1bTS63eo5fkYY4ypx0JVxJgI/KSOylGVrl4/2D2qISJyn4g8KiJ3ishwEbEnCWNo7dq1PPDAA1x66aUcd9xxJCQkICLMnz+/WvkLCgq45557OPnkk8nMzKRx48Z07tyZcePGsWzZsrDKkpOTg4hU2S1evDiCJa3MNz1jTPTV+7ZNROQoXAAFd98rkMsCDFsjIher6ue1UjAT0sMPP8wDD4RVibTU+vXrGTZsGOvWraNNmzYMGjSI1NRUNmzYwIsvvkjPnj0ZMGBA2NMdPnw4Rx11VNDxocYdKSZOnMjcuXN58sknmThxYqyLY0wl9TpoiUgSMA/IAN5U1ZcqJFmFq2n4JvAd0AzoA/we6Am8ISJ9VHVTkOlfCVwJ0KlTp1pZhiPVCSecwA033EC/fv3o27cvV1xxBUuWLKky3/79+znrrLP49ttvueWWW7jllltITk4uHb9z50527twZUZluuukmBg8eHFFeY0z9UK+DFvAX4AxcCxuXVhypqvdXGLQfWCgir+Mehj4VmAFcE2jiqvoo8ChAv379NHrFNpMnT44o3//93//x7bffctlll3HHHXdUGt+yZUtatmxZ0+IZY+JUvb3vIyIP4GoM/gicoarVafcQcJVEgLu9r+fWQvHq3IKVmxgw8y0637SQATPfYsHKgCePce3w4cM89thjgDsriqXFixcjIkHPzDZs2ICIkJOTE9Z0CwsL+ctf/sLAgQNp3rw5jRo1omvXrkybNo3t27dXSp+bm4uIMHHiRHbu3Ml1111H586dSUlJYdSoUdWa53/+8x9GjBhBmzZtSE5OpkWLFhx33HFMmjSJTz75pNzyzJ07F4DLL7+83P2+3NzcctPcuXMnN998MyeeeCJNmzalSZMm9OnTh9mzZ1NYWFipDBMnTiydzqpVqxg1ahStWrWicePG9O3blyeffDJg2Q8ePMjMmTPp06cPTZs2LW1P8LTTTuPmm2/m4MGD1VoHpuGo6kzrgmCtUVRBVfWYCPIBICJ/BK4DtuMCVrgPKwP4WsNoH2k56osFKzcx44XPKSh0rz/ZlFfAjBfcrbpRveN+8UqtWLGCnTt30rFjR44//njee+89Xn75ZXbu3MlRRx3F2WefzWmnnRbrYkYsPz+fESNGsHTpUjIyMujbty+ZmZl88sknzJ49m+eff54lS5YEDIQ7duzg5JNPZs+ePQwcOJB+/fpV64wzNzeXyy+/nISEBPr37092djb79u3jhx9+IDc3l2OPPbY0IEyYMIGlS5fy7bffMmDAALp06VI6Hf/Pn3/+OWeffTabN2+mQ4cODB48mJKSEj788EOmTZvGwoULeeWVV0hJSalUng8//JApU6bQvn17zjrrLLZt28aSJUuYNGkSK1eu5MEHHyxNW1JSwogRI3jrrbfIyMhg0KBBZGRksHXrVtauXcvvf/97rrnmGrsXeaRR1YAd7oHfSLviYNOtqgPuwbW0sQPoWYPpnOZNZ2d10vft21era82aNdVOGw2n3/2mZt/4cqXu9LvfrNNy1MSgQYMU0Oeeey5omkceeUQBHTp0qE6YMEG97VeuGzt2rB44cCCseWdnZyugb7/9drXzvP322wrooEGDAo5fv369ApqdnV1pnK+sFV100UUK6AUXXKC7du0qHV5UVKS/+c1vAs7vySefLJ3esGHDND8/v9rLoKrauXNnBXTZsmWVxv3www+6evXqcsN86/3JJ58MOL0DBw6UTvOuu+7SwsLC0nE7d+7UM888UwG99dZbA04X0Ouuu06LiopKx33wwQeanp6ugC5cuLB0+JIlSxTQPn366L59+8pNr6SkRJcuXar79++v1nqo62P2SAIs1wh/pyPpqro8+BowJIJuaHUCZkUiMhO4Aff247O06nd3hXKh1/+4BtOoFzbnBW7sM9jweLVr1y4A3nnnHZ566immT5/OunXr2L17Ny+++CLt27fn+eef5+qrr45o+kOGDAla3T0zMzOai1LJmjVr+Oc//0l2djZPPfUUzZs3Lx2XmJjI3XffzUknncSSJUv4/PPKFV6Tk5N55JFHSE9PD2u+W7duJTMzk9NPP73SuA4dOtC9e/ewppebm8v69eu58MILmTFjRrmXK7Zo0YK5c+eSnJzMnDlzfH8ey8nKyuKee+4p13Zn//79+fWvfw3A7Nmzy5UdYODAgTRp0qTcdESEAQMG0LhxY8yRparLgz+qatVVvqJARO4EbgTycAFrZRXpe+EeNn5VVYv9hifhLi1e5w2aHSB7XMnKTGNTgACVlZkWg9LUHl+jpkVFRUyePJlZs2aVjvvZz35GVlYWp5xyCnPnzuXmm2/m6KOPDmv6oaq81/aP36uvvgrAyJEjSUurvN0SEhL4yU9+wmeffcb777/PiSeWf49qnz59wr5/BnDKKaewePFiLrvsMn7961/Tq1evGj1D9sorrwAwbty4gOOzsrLo2rUra9as4ZtvvuHYY8u/zWjcuHGkpqZWyjd+/HjuuOMOli5dSlFREUlJSfTp04fExEQef/xxjj32WMaOHUvbtm0jLrtpGOpF7UER+Rlws/d1HXBtkAPrK1Wd6X3OAf4F7BKRr3FvTE7HvTU5C3eZ8kZVXVSLRa8TNwzvVu6eFkBaciI3DO8Ww1JFn/9ZxC9/+ctK433V55cvX87ixYvDDlqxrPL+3//+F4A5c+YwZ86ckGkDVcjIzs6OaL4PPfQQI0aM4Omnn+bpp58mIyODU045hbPOOovx48eHfT/ItxzBgpa/7du3VwpanTt3Dpi2U6dOJCQkcPDgQXbu3Enbtm055phjmD17NtOnT+fqq6/m6quv5uijj+b000/n/PPPZ/To0fa2hSNQvQhaQAu/z/28LpAlgC9ofYp7BcopQDbQG3fNfCPwJDBHVVfUSmnrmK+yxaxFa9mcV0BWZho3DO/WoCphAOXOJIL9uHXu3Jnly5fz44/VrkxaK8J91UVxsfvD0bdvX0444YSQaXv06FFpWKCzs+o4/vjjWbt2LYsWLeKtt95i2bJlvP3227z++uvcdtttPP/885x99tnVnp5vOUaMGEGrVq1Cpo300QT/P6zXXnst48aNY8GCBSxdupSlS5cyb9485s2bR69evViyZAnNmjWLaD4mPtWLoKWuhfbcMPOsB/6nNspTH43q3b7BBamK+vTpU/p5586dtG7dulKaHTt2ANT6Szh9Nd/27dsXcPx3330X1vQ6duwIuPtq/pc960JycjIjR45k5MiRAOzevZvbb7+dBx54gCuuuIJNm6r/+ETHjh1Zu3YtU6ZMYcSIEWGXZcOGDQGHf//995SUlNCoUSNatGhRbtxRRx3FVVddxVVXXQXAp59+yvjx41m1ahUzZ87krrvuCrscJn7V2+e0zJGnffv29O/fH4A336z88uvdu3eXPlfUr1+wk/HolQXg22+/Dfjcke/eTnWdc845ACxYsICioqKaF7AGmjdvzqxZs0hISGDz5s3lLkf6gnWwMvqW47nnnoto3s899xyHDx+uNPyZZ54BYMCAAeUqdwTSs2dPrr/+esAFMHNkCRq0VDVBVYO28m5Mbfjtb38LwB133MGqVatKhx88eJApU6awZ88e+vbtW+vPa2VnZ3PMMceQl5fHvffeW27cggULyj1PVB19+vRh1KhRrFu3jgsvvJCNGzdWSrNlyxbuv//+qAW1AwcOcN999wW8R7Zw4UJKSkpo1qxZuZqTvmD95ZdfBpzmlVdeSceOHZk7dy633norBw5UfsHDF198EfRh4U2bNnHTTTeVu7z68ccfc9999wGUBiOAt956i1deeaXS+iguLi790xDpvT4Tx+qyfn197urzc1rxaMWKFdq/f//SzvccTteuXcsND2T69OkKaEpKig4cOFBHjRqlWVlZCmj79u3166+/Dqssvue0hg8frhMmTAjaLVq0qFy+Z599VkWk9FmhCy64QHv27Kkior/97W/Dfk5rz549pc+rNWrUSPv3768XXnihnn322dqjRw9NSEhQQAsKCkrz+J7TmjBhQljLrKq6e/duBTQxMVF79eql48aN04suuki9JssU0IcffrhcnpUrV2pCQoImJCTosGHDdNKkSXrFFVeUe87rs88+006dOimgLVq00MGDB+vFF1+sQ4YM0ZycHAUqbVvfc1pXXXWVpqamapcuXfTiiy/WM844Q5OSkhTQqVOnlssze/ZsBTQjI0OHDBmil1xyiY4aNUrbtWungB511FG6YcOGaq0LO2ZrD3X8nFbMg0V96SxoRZfv4dyqumBeeOEFHTJkiGZmZmpKSop26dJFp02bptu2bQu7LL6gVVU3e/bsSnlffPFFPfXUUzUtLU3T09N10KBB+tprr0X0cLGqe5D4qaee0mHDhmmrVq00KSlJW7durT179tSpU6dWCpw1CVqFhYX68MMP60UXXaTdunXTZs2aaVpamnbp0kUvueQS/eCDDwLme+6557R///7atGnT0mWp+LBxXl6e3nXXXdq/f39t1qyZpqSkaPv27fXUU0/VW265RT/99NNy6f0fWl6xYoWOHDlSmzdvro0aNdLevXvrY489piUlJeXyrFu3Tm+99VYdMmSIduzYUVNTU7Vly5bau3dvvf3228PaF+yYrT11HbTEzdP069dPly9fXq20X375Jccff3wtl8iYhiPWrzyxY7b2iMgKVa3dm8x+rCKGMcaYuGFByxhjTNywoGWMMSZu1IuHi40xDVtubm6ld3IZEwk70zLGGBM3LGgZY4yJGxa0jDHGxA0LWsYYY+KGBS1jjDFxw4KWMcaYuGFByxhjTNywoGWMMSZuWNAyxhgTNyxoGWOMiRsWtMwRKScnBxFhw4YNdTI/EUFE6mReFW3YsAERIScnJybzNyaaLGgZU0MTJ05ERKxtPWPqgDWYa0wd+PLLL2NdBGMaBAtaxtSB4447LtZFMKZBsMuDplZ8+OGH3HDDDfTr14+2bduSkpJCVlYWF1xwAR988EGl9Lfddhsiwm233cbWrVv51a9+RYcOHUhNTaVz587cdNNNHDx4sFK+vXv38uijjzJq1Ci6dOlC48aNadq0Kb179+b3v/89BQUF1SqvqtK1a1dEJGD5fMaMGYOI8NBDD5XeK5o7dy4Al19+eem9q4qXC0Pd0yosLOTRRx9lyJAhtGjRgtTUVDp16sTIkSN55plnyqX97rvvuPvuuxkyZAgdO3YkNTWVFi1aMGTIEP72t79Va1mNiWuqap0qffv21epas2ZNtdMeqc444wxNTEzUk046SUeOHKljx47VE044QQFNTEzUZ599tlz6W2+9VQGdNGmStm/fXrOysvSCCy7QYcOGaePGjRXQ8847r9J83n33XQW0TZs2OnDgQL3ooov0zDPP1PT0dAX0lFNO0YKCgkr5srOzFdD169eXDps9e7YCOn78+IDLtHHjRk1KStL09HTNz8/X7du364QJE/SYY45RQAcMGKATJkwo7d59993SvIC6w628Xbt26WmnnaaApqam6tChQ/Xiiy/Wn/70p5qZmanZ2dnl0t95550K6DHHHKNnnnmmXnTRRTpgwABNTExUQK+77rpK81i/fr0ClaZ1JLFjtvYAy7UOf6vt8mBtuC0j1iWI3G17ojKZ6dOn88wzz9C2bdtyw1966SXGjh3LVVddxYgRI2jcuHG58U888QSTJ09mzpw5pKSkAO5+0CmnnMJLL73EsmXLGDBgQGn6nJwc3nzzTQYPHkxCQtmFg7y8PH7+85/z2muv8cADD3DjjTdWWebLL7+cW265hWeffZb77ruPVq1alRv/yCOPUFRUxGWXXUZ6ejrp6enk5uYyceJEvv32WyZPnszEiRPDWk8TJ07k/fff57TTTmP+/PlkZWWVjjt48CBvv/12ufTDhw9n9OjR9OjRo9zwb775hjPOOIMHH3yQSy65hP79+4dVDmPihV0eNLXi7LPPrhSwAM477zzGjRvHrl27Kv0gA3Ts2JEHH3ywNGABHH/88YwfPx6AN998s1z6Dh06MHTo0HIBCyAzM5MHH3wQgPnz51erzBkZGYwfP55Dhw7xxBNPlBtXWFjIY489BsDUqVOrNb2qrFq1in//+980bdqUF198sVzAAmjUqBHnnHNOuWEnn3xypYAF0LVrV2655Rag+strTDyyMy1Ta3bs2MHLL7/MF198QV5eHkVFRQB88cUXAHz99deMGDGiXJ6hQ4eSlpZWaVq+igybN2+uNE5VWbZsGe+88w4bN26koKCg9FKCbz7Vdc011/Dwww/zl7/8henTp5cGwxdeeIEff/yRwYMH071792pPL5TXXnsNgPPPP5/WrVtXO9/BgwdZtGgRH3/8Mdu3b+fQoUMAbNmyBQhveY2JNxa0akOULrHFs0ceeYRp06Zx4MCBoGny8/MrDevUqVPAtM2aNQOoVBlj69atjBkzhvfeey+s+QTTvXt3zjzzTN544w1ee+01zj33XAAeeughAK6++upqT6sq3333HRBezcL333+fCy+8kI0bNwZNE87yGhNv7PKgibrly5czZcoUCgsLmTVrFl999RX79u2jpKQEVWXGjBkApWdC/ipe5qvK5MmTee+99xgwYACvv/4627Zt4/Dhw6hq6RlIuK699lqgLFCtXr2ad955h6ysLEaNGhXRNKPhwIEDjB49mo0bN3LFFVewfPly8vLyKC4uRlVZtGgREHi9GtNQ2JmWibr58+ejqlx33XVMnz690vh169ZFZT779+/nlVdeITExkZdffpnMzMyozGfkyJF07tyZV199lQ0bNjBnzhwArrzySpKSonfIZGdnA7B27dpqpX/nnXfYunUrffv25a9//Wul8dFar8bUZ3amZaJu165dgKtUUdH27dt5/fXXozKfPXv2UFJSQnp6eqWABVR6xqm6EhISmDp1KiUlJcyaNYt58+aRlJTElVdeGTC9r9KI755ddQ0fPhyAF198kR07dlSZPtR6Bew5LXNEsKBlos53j+app55i3759pcP37t3LpEmTyMvLi8p82rZtS/PmzcnLy6v0g/3aa69x3333RTztK664gsaNG/PQQw+xd+9eRo8eTbt27QKmbd++PRB+U029e/fmvPPOK52+ryKFz8GDB3n11VdLv/vW61tvvcVXX31VOrykpIQ77riDZcuWhTV/Y+KRBS0TdZdffjkdO3bkk08+4eijj2bMmDGMHj2anJwcli9fzqRJk6Iyn8TERH77298C8Itf/ILTTz+99Bmlc845h2nTpkU87ebNm3PppZeWfg9VAeP8888nISGB+++/n+HDh3PFFVeU3murSm5uLieffDJLly7l6KOP5qyzzuKSSy5h8ODBtGvXjilTppSm7dOnD+eddx75+fn06tWLc845h4svvpiuXbty55138pvf/Cbi5TUmUTQ8JgAADc5JREFUXljQMlHXvHlzli9fzpVXXknTpk1ZuHAhy5cvZ8yYMXzyySdBL29F4v/9v//H/PnzOfXUU1m9ejUvv/wyiYmJzJs3j9///vc1mvZZZ50FQI8ePRg0aFDQdL169eKf//wnJ598Mu+99x5PPPEEjz/+eLWqnrdo0YJ3332XP/3pT/Tp04ePPvqIF154gfXr1zNw4EBmzpxZLv38+fOZOXMmXbp0YfHixbz55pv06NGDpUuXVnqmy5iGSKymkdOvXz9dvnx5tdJ++eWXHH/88bVcIhNro0ePZsGCBTz00EPlznhM/LFjtvaIyApV7VdX87MzLWMCWLFiBf/+979p2bIll112WayLY4zxWJV3Y/xMnjyZffv28corr5RWcGjSpEmsi2WM8VjQMsbP448/TkJCAtnZ2fzud7+LWjuDxpjosKBljB+7x2tM/Wb3tIwxxsSNBhO0ROQSEXlXRPaIyD4RWS4iV4tIg1lGY4w50jWIH3QRmQM8A/QD3gVeB44F/gzMF5HEaM/TLiMZEx/sWG1Y4j5oichYYCrwI3CSqo5U1dFAV+BLYDRwTTTnmZiYSHFxcTQnaYypJSUlJWG/PcDUXw1hS87w+jeq6je+gaq6FfA9EXpTNC8TNm7cuFybesaY+uvAgQMBXyxq4lNc1x4UkQ5AX+Aw8FzF8aq6REQ2Ae2BU4GqG4OrhmbNmrFjxw7S09NJTIz6lUdjTJSoKnl5efasXRALVm5i1qK1bM4rICszjRuGd2NU7/axLlZI8X6m1dvrr1bVgiBpPq6QtsbS09Np0qQJ3333Xelr5O26uTH1h+8loFu2bKGoqIjmzZvHukj1zoKVm5jxwudsyitAgU15Bcx44XMWrNwU66KFFNdnWkBnr/9diDTfV0hbYyJCmzZt2Lt3L/n5+Wzbts3ucRlTzyQlJZGRkUGbNm3snlYAsxatpaCw/O9WQWExsxatrddnW/EetJp6/f0h0vhuPqVXHCEiVwJXAnTq1CmsGYsIzZo1o1mzZmHlM8aY+mBzXuCLU8GG1xfx/vdDvH5E1+ZU9VFV7aeq/Vq3bh3FYhljTP2WlRm4ckqw4fVFvAetvV6/aYg0vnF7Q6Qxxpgjyg3Du5GWXL4iWVpyIjcM7xajElVPvF8e3OD1s0Ok8b1xcEOINMYYc0Tx3beKt9qD8R60Vnr9HiKSFqQG4ckV0hpjjMEFrvoepCqK68uDqvoD8AmQAoyrOF5EBgEdcK1lvF+3pTPGGBNtcR20PHd7/T+ISBffQBFpAzzkfZ2pqiV1XjJjjDFRFe+XB1HV+SLyMK7Jps9F5A2gEDgDaAYswDWca4wxJs7FfdACUNWpIrIUuBoYBCQCXwFPAA/bWZYxxjQMDSJoAajq34C/xbocxhhjao9Ym3mOiGwndHNQobQCdkSxOCY8tv5jx9Z97NSXdZ+tqnXWOoMFrSgQkeWq2i/W5ThS2fqPHVv3sXOkrvuGUHvQGGPMEcKCljHGmLhhQSs6Ho11AY5wtv5jx9Z97ByR697uaRljjIkbdqZljDEmbljQMsYYEzcsaNWAiFwiIu+KyB4R2Sciy0XkahGx9VpLRCRXRDRE91WsyxjPRKSbiFwvIvNE5CsRKfHW6wXVyGvHQw1Esu6PxOOhwbSIUddEZA4wFTgIvElZe4d/Bs4QkXGqWhzDIjZ0y4B1AYZvqeuCNDBTgOvDzWTHQ1REtO49R8zxYEErAiIyFneA/gj8VFW/8Ya3Bd4GRgPXAA/ErJAN319VNTfWhWiAvgBmAcuBFcDjuPY8g7LjIWrCXvd+jpjjwYJWZGZ4/Rt9ByiA6v9v7/6DpazqOI6/PwIGikZikigJ5Uwh5Y8wsmwQzGkM+cMpMCeNHHUqrXEmpkmsP4oZKmwybURkzHEqf0wz5I9UZhpg6pKWlYgQYpqapJhSWPFbofz2xzkL67J77+69yz53dz+vmWfO7Dnnee53L5w9d59znnNis6QrgB5gnqQbvVivtZOIuLX8taR6TnN7aIJ+/u67ju81N0jS8cBkYA+wtLI8IlYBLwHvAM5obXRmreX2YK3mb1qNOy2nGyJid406jwLH5bq/a0lU3We6pJOBkcBm4GFghf+Sbzm3h8Gha9qDO63GTchpbyvCv1BR15pvTpW8JyVdGBHrWx5N93J7GBy6pj349mDjRuZ0Zy91duT0iIMcSzdaC1wFTCL9W4wFZgLrgJOAlZKOKy68ruP2UKyuaw/+ptW40uio178qQETcUJG1E1gmaQWwijRucg1ptpodfG4PBerG9uBvWo3bntORvdQplW3vpY41UUTsAb6bX84oMpYu4/YwCHVye3Cn1biNOT2hlzrjKupaa5Se/u+o2yGD3Macuj0MPh3ZHtxpNe7xnE6SNKJGnQ9W1LXWGJ3THb3WsmZyexi8OrI9uNNqUES8CKwBDgVmV5ZLOgs4nrQ6wCOtja7rXZDTRwuNoou4PQxqHdke3Gn1T+le8bWSTixlSjoGWJxfLuzEZySKJOlUSTMlDanIHyppLmkWFcD1rY+uq7k9FKBb24M3gewnSYtJC1y+Bqxk/wKhRwL3AbO8QGhzSTofuBf4F/AXYBNpGvX7SVN93wCuiYjvFRZkm5P0AfZ3NJCmTR8BPEP6vQMQEWdUnOf2MECN/u67tT240xoASZ8BvkT6TzKENPB5G3Cz/6psPkkTSKtgTyEN/I8mTbXeBDwE3BQRjxUXYfuTNI20yG2vIuKAhfHcHgam0d99t7YHd1pmZtY2PKZlZmZtw52WmZm1DXdaZmbWNtxpmZlZ23CnZWZmbcOdlpmZtQ13WmZm1jbcaVnXkTRNUkjqKTqWgZJ0dX4v51bk9+T8aQWFVopjrKRdkn5eZBzWOdxpWceRtDF/YI8vOpaDSdKxwDeA30TEL4uOp5qI+DtpaaJPFd2BWmdwp2Xd6I/ARGBO0YEM0HzSWnPziw6kDwuB14HvFx2ItT93WtZ1ImJXRDwVES8UHUt/SRoNfBb4K3WsV1ekiNgCPABMlnRm0fFYe3OnZR1D0iWSgv276D6fbxNG+e3CWmNaksbn/I2SDpE0V9IGSbslbZL0A0mH5bpvk3RDrvu6pGfydhC1YpOkCyUtl7Qln/OCpB/18zbmpcBw4KfRxwKikiZLul/Sq/m9rJN0WY26+8bCJE2VtCzH+0ZeVRxJwyXNk7RG0o78Xl6W9IikBZKGV7n0T3J6ZT/eq9k+Q4sOwKyJniV9OM4CDgfu5s27tjayg+tdwEygJ193KvAVYKKki4Dfk27NPQwclcuvkzQ8Ir5TfiFJw4CfAZ8EdgOrgc3A+4DLSeM9H4+I1Q3Ed35OV/ZR71xgLvA0sBx4J/AR4FZJoyLiuhrnzQa+CDwJrACOBvZKOgRYBpwNbAVW5XQM8B7SGNsi0qaP5X4N/A+YKWmItymxfosIHz466gA2krZoGF+jfFou76nIH5/zg7StxtiysnHAlly2HlgKDC8rPy+XbQMOq7juwly2Cji+ouzLuexZYGid7+8wYE8+hteo01P2Xi6tKLs452+tEmv5eZ+vct2puewx4PCKMgFnVl6zrPzxfO6Uov+P+Gjfw7cHzaq7KtLMN2DftvJ35JcnAFdExGtl5cuAP5G+fZ1eypd0FGkH2R3A7IjYVP5DImIR6ZvLu4FP1BnbJGAY8Hx5DDXcHRG3VfzMO4A/kzZoPL3qWbAiIm6pkj8mpw9FxM6K60ZE/DYidtW45pM5Pa2PmM1qcqdldqC9wK+q5D+b09WRJhdUeianY8vypgMjgFUR8Y8aP29VTj9cZ3zH5PTVOuo+WCP/qZyOrVF+T438NaTbfJdJulLSmBr1qintvtvIOWZv4k7L7ECvRMR/q+SXxsQ2VSkrLy+fiPCunJ5XMSlk3wGUtkN/e53xvTWn2+qoW2uGZOncapMmAP5WLTMiniON7R0K3AS8Iuk5SbdLmiVpSC+xlH7mqD5iNqvJEzHMDtTX1vCNbB1f+hB/mjR5ozd/qPOa/8npkXXU7e8297trFUTEjZKWkiaDfDQfF+djraSzIqJah1qK99/9jMnMnZbZQfZiTtdHxCVNumbpNuPoJl2vYRHxCrAkH0g6BbgdOBWYB3y9ymmleGvdJjXrk28PWifak9PB8EfZStIY2TmSmnVbbANphYkJkkY06ZoDEhHrgB/ml6fUqHZSTtcc/IisU7nTsk70Uk4nFhoFEBGbSWM/o4D7Jb23sk5+UPnyeic1RMRu0q3EYcDkZsbbF0lnS5ohaWhF/hBgRn55wHiYpJGk59K24U7LBmAw/CVq1mz3kp7FulPScvaPAV0dEfXMuGu2r5Fm6V0APCFpLfA8aRLEOFLnemhON9d5zftIz0ydQ3rAuVVOBq4HtkpaA7xMem7sQ8CxpIeKr61y3nTS+N6D4QeLbQDcaVknWkQa9L+ItKrFW3L+AuqbJt5UEbEX+LSkO0nLL00hffhvJ33o3wX8Aniugcv+GPg2MEfS/IjodSmnJnqA9K1xKnAiaXWNHaRZikuAmyPin1XO+1xOF7ciSOtcat3/dTNrJklLgC8AH4uIas+VDQqSjiY9JvBERNR6mNmsLh7TMmtf3yJ9W/tmwXH0ZR7p2+5Xiw7E2p87LbM2laedLwCmVu5cPFhIGkta2f2eiOgpOBzrAL49aGZmbcPftMzMrG240zIzs7bhTsvMzNqGOy0zM2sb7rTMzKxtuNMyM7O28X9EuGg+L8v+XgAAAABJRU5ErkJggg==\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"plt.plot(t,temp_numerical,'o',label=str(n)+' Euler steps')\n", | |
"plt.plot(t,temp_analytical,label='analytical')\n", | |
"plt.title('Rate of heat loss of a corpse')\n", | |
"plt.xlabel('time (hrs)')\n", | |
"plt.ylabel('Temp (F)')\n", | |
"plt.legend()\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 142, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"As t goes to infinity, the temperature approaches the ambient temp\n", | |
"This can be seen in the graph above, as the lines approach the limit of Ta\n" | |
] | |
} | |
], | |
"source": [ | |
"print('As t goes to infinity, the temperature approaches the ambient temp')\n", | |
"print('This can be seen in the graph above, as the lines approach the limit of Ta')" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 167, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"1.6799999999999997\n", | |
"Time of death is 1.8865228851460631 hrs before T0\n" | |
] | |
} | |
], | |
"source": [ | |
"Ta=65\n", | |
"T0=85\n", | |
"dT=(98.6-Ta)/(T0-Ta)\n", | |
"print(dT)\n", | |
"tdeath=-(np.log(dT))/k\n", | |
"\n", | |
"print('Time of death is',-(tdeath),'hrs before T0')" | |
] | |
}, | |
{ | |
"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": 168, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"time=[8,9,10,11,12,13]\n", | |
"Tempa=[55,58,60,65,66,67]\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 206, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"x=[8,9,10,11,12,13]\n", | |
"y=[55,58,60,65,66,67]\n", | |
"\n", | |
"def tempattime(t):\n", | |
" \n", | |
" indexvalue=x.index(t)\n", | |
" tempa=y[indexvalue]\n", | |
" \n", | |
" return tempa\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 208, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"65" | |
] | |
}, | |
"execution_count": 208, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"tempattime(11)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 210, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"Text(0, 0.5, 'temp F')" | |
] | |
}, | |
"execution_count": 210, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAE0CAYAAAAR5SXAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3gc9bX/8feRZbl3uWDcbbnQXBOqTQstcBMInQAhJEASagKENLj8QnJjEgjNJCGEhEtPwAFuCr3aBkNcMM02cm+4SO5Fliyd3x8zu1rJuyutrNWuVp/X8+wzmpnvzJ7V7O7ZmfnOGXN3REREGltepgMQEZHcpAQjIiJpoQQjIiJpoQQjIiJpoQQjIiJpoQQjIiJpoQSTAWb2mJl5+Ph7puOJZWbLwrhu3Yd1RF7bJY0XmUh2MLNjYt7jgzIdTzZTgmliZtYROD1m0qlm1j1T8bREZvZw+OXwZiOsa58TsmSHxnxfSEAJpumdBXSIGS8AzstQLCIiaZOf6QBaoIvD4RJgNzAqnPa7jEXUyNzdMh2DSLq4+5uA3uP1oD2YJmRm/YGjw9FHwwfAoWZWlJmoRETSxN31aKIH8BPAw8cwoD9QFY7fVseyD4ft3gzHvwD8DVgD7ALmh+tvG7NMV+AW4BNgB1AK/B04IMnzLAuf59Zw/JvADGBjuI7ZwLVAfpJ1RF7jJUnatAOuA94CNgDlwFrgOeDLSZa7NVz3snB8EPD7MO7dwDrgGWBsnGUviYkt0ePhem7Lh+uxrr1eP8FRg0uBF8PXWx6+/peBrwOW4PmisYfjfYF7gMXh9l8OPAD0jVkmL3yud4BNwHZgOvBfKbzPJoXbZC1QBiwC7gC6N+D93yp8vzrwQD3avx62nVFrugHnA/8GPg//h1vC2F4BfggMSCGulN4XwDEx0wcle/+H2+AKYCawmeBz9AZwQq1lxgFPAitj/s+3Ae3qEf8E4KFwmR3ANuBD4H+AwlS3U2M+MvbELfFBkAQceCdmWuRDtDTRl0vYLvrBB74BVCT4ILxA8CU2CFiQoM1WYHSC51kWtrkVeCzJB+4toEOCdSRNMMAhMc+T6PFnoFWcZW8N5y8j2BvcnGD5MuBLtZZN6Yukjm35cD3WdUmtZQYAH9SxzL/j/V9rxT4GWJ9g+WXA/kAbgsSQ6Hkurcf77AqgMsHyq4ERDfgM3BkuvxEoSNJu/5jn/m7M9FZ1vK7I44YUYkrpfUH9E8y3gX8lWF8VcFHY/mKCJBmv3evE+RyEy+WF/8+qBMs6wY+XQzP2nZepJ25pD+CLMRs99gPzzZjpRydZPvLBX03wS/0Vgl+XPYDhwIMx67kCeJ/g1+I3gX5AL+AiguTixCS5Ws+zjOqE5wRJZnz4POOBx2Oe59EE64j7BRvOG0CwJ+UEv7ovB4YC3YADgduBPeH8X8RZ/tZw3iagBJgHnAH0AXqHH9ZNYZsVxOxpESTejlQnzmnheOyjTT23Z5uw/fJwXf8TZ12xz92Z4Bemh3H/ABhJsJdZBPwI2Bn5n8d5vkti/q9LgGLg7PB170+wVxn50fEIcHf4PvlvYATQHTgW+DRssw3okeR9torgS+994ASgkGCv+7+p/jL8jJg95nr+38bGvI4zkrS7MWxTHhsn8K2Y5X8HHArsF/4fvhBu/38B16YQU0rvC+qfYJaE8d8asw2Oo/qH5pbwf1tOsEc7keBzVgT8MWY9lyeI+zcxbf43XL6Q4HNwOvBROG89sF9Gvvcy8aQt8QFMCTf27lofmE4Eu7UOPJRk+Ydj3kz/JP6v++nh/IrwzTs0Tptvx6xnZJz5y2LmP5gglj/HtBkfZ36yBPN/VH+B9Uyw/siXyG5iDvmE826NWf9c4v/a/1pMm5OT/C/fbITtGvl/3VpHu3vDdltJ8Msf+FJM3BNqzbskZt4K4hz6AH4Rzt9D8Kv2rDhthlO9Z3BFHe+zD4D2cdpcHNPm+gb8zz4Ol52apE1kT++5WtOnhtOf3ddt19D3BfVPMA6cm2AbRPY6KoDngbw47aaFbd6NM298zDp+nCDOjlT/oJjS2P+v+jx0kr8JmFlr4Nxw9N/uXhqZ5+7bCN5gAGeZWbt6rPL77l4ZZ/pT4TAfuNfdF8dp81eCNxwEe1WJlBEcy47nRoIvfwj2kOrFzIYAp4Wj17n7hgRN/0xwbqGA4Fd6Ije5+444058j2IuB4FdtRplZB4KkCfBzd18Yr527v0pwSASC8zGJ/NzdS+JMj2z/VgTnLZ6J8xyfESRmCH79J/Mjd98ZZx2PALPC0Xpv/xiPhcNTzaxr7ZlmdiAwulbbiEjP19UNeN6mNsPd/1p7YrgN5oSj+QRJuirO8pFlx5lZ7R6/1xCci/oImBzvyd19O/CrcPR8M2vynm9KME3jVIJdV6juORbrkXDYGfhqHeta4u7FCebFJpSX4jUIE1rki71Pkud50903xZsRJsi3wtEjk6yjtuMJPhRVwDtm1jHeg+A6oXnhMhMSrGs3wcnSePFVERyOguSvsakcAbQP/34r0esOX/uHYbtErxsSbFvqsf1rtUv2v9lBcBg2kUgFigPiJYk6RA6ztiG4Lqy2C8PhFoK99ViR5HipmV1kZm1SfO6mVJ9tsNjdF9XRpoDg8FqsL4XDN4AOSd5Pn4btugNDUox/nynBNI3ItS+bCI4P1/YKQS+d2LaJrEkyb1fM35/Xo12yvaUFdcQxPxwOrKNdrBHhMI/gF+i2JI+vhW17JljXBnevSPJckV/e7ZO0aSojYv5+n+Sv+7qwXaLXDQneA+7emNu/OMFeckRk+xvBebV6c/eVVP9AuTB2Xvgr+4Jw9Bl3L6u1+F0E573aEfwwKzWzl8zsZjObaGatUoklzerzWa3PdoKYbRUmjr7h6DUkfz/NillHsvdUWijBpFlYBubUcHQGMMrMxsQ+gIOBt8M2J5pZ7ySrTPahT7Vdsl3m7XUsG5nfqX7hANAlhbYRbRNMr+//IRsuiGvM100dX/wRTbX9IbX3QERkT35SeH1YxESqE9Zee/vuvoXg0O59BD/YOgAnAj8n+AytNLNrzSwbvtsaaztBzW3VkPcTJHlPpUs2bIRcdx7BLi4E5x/mJnicE7ZpRfLj702lYz3nb0thnZEvpRJ3t3o+jkk58uxT48u4nq97UKaCDdV3+0Nq74GIZwjO8xk13++RPZoVVP/oqsHd17v7NQS/yCcAVxMcsttJ0KPsboLuu7kq9v10VQqfpTebOlAlmPSr65BXPBc1ehSpG1nH/FHhcHkK61wSDgvNrF/qITVbS2L+HpOxKFJTVMfhpsj2j/RqS4m7byXoUQhhgjGzAqrPyTzhYVeoJOuodPfZ7j7F3c8kuHB5ejj7ajPrkWpczUG4FxfpKDQ2k7HURQkmjcLyL5GeOr+q6xcGQe8sgDFmdnBmoo46xsy6xZsRHvY7OhydkcI6Y08aN6T3UWOJnLtpjOP19VnXWzSg112GdSC4RiORyDmyT919cwOfI9JD7CAzG01wKLlbrXn15u4bgd+Go60IugOnojHfF+kW+Sx9zcw6ZzSSJJRg0usbMX8/WY/2TxH0sIKG7fk0prYEFz3GcwdBDyCAv9R3he6+gOBKdYAfm9kRydqbWa9ESW4fRX799U3aqpHWFf5afygcvcTMzky2QjPrbGb7NUJs+2qyme3VScLMLqa6l1u9t38cLxJcdArBobHIobK57v5JvAXMrK4966Exf5cmbBVfY74v0u2ucNgNeCi8FCIhM0s12TYKJZg0CXvDRD4wn7j7R3Ut4+6rCC6uAvh6hnvELAMuM7NHzWycmXUPh49T/Sv8MXefneJ6v0NwZXE74A0zu9PMDjOzQjPrYWajzOwCM3uS4PDb0KRra5hIzEPM7MowkeWHj1Q/E5F1nW5mJ5pZl5h1xZ6Y/QmwkOAz97SZPWRmR5tZbzPrZmZFZnammf2JoB5VKt2/02ENcADwppl9Kdw2Q83sFuBPYZti4P6GPkHYCzByrceFVHeGSbb38oKZzTSz683scDPrE743DzSznxJcbAowJ7zeJBWN+b5IK3d/n+ofgGcB75nZ181ssJl1NbP9Lbgx2s1m9gnVe3ZNHqgeaXhQ82rfn6Sw3BUxy50UM/1h6rjKmDquMI5pt4wEV59TsxZZbFmY2o99qUU2kuorjJM9qqhVM41axS6TvMY3w3YPx5nXjuAag3jPuVf7Op7nIILDX/HWdUmttvtRXW2hrsdXai17SWReHfEk/d/X9V6iZi2y79DItcjiPN+htda7hyRlTai7hp0TlDka3oBY6vW+IIVilw3ZBik8j4Wfhz0JYo59JKyakM5HVmXlHBN7or4+h8cinqb6WHCmD5NdCFxGUAl2C0EvnbkE12oc7/Gvoq+TB4fKDiE4hPgPgl/L5QRf1CsJDp1cA/R393mJ1tNQHlwvMhH4A8EFmbWvtUhlXR8T1IT7O8HrSHhtjrt/Hj7vGQTbeUX43OUE10O8QVCTrMjd/y/RepqKu/+B4IK+fxLsde4m+AL+LXCwJ6hIkOJzvEdQ0yzi9fD/lMhJBO+N5wiuxdlE8AVbQtDr7HrgIE9976VR3xdNwQO3EnS4uIvg4uQtBD8KtoTjvyPYhucmWE1aWZgJRUQws4cJEv9bnhtdxCWDtAcjIiJpoQQjIiJpoQQjIiJpoQQjIiJpoZP8ocLCQh80aFCmwxDJqGXLllFaWkrHjh0ZMWJE3QtIizd79uwSd49bqbn2TWxarEGDBjFr1qy6G4qISJSZJaxHqENkIiKSFkowIiKSFkowIiKSFkowIiKSFkowIiKSFkowIiKSFkowIiItVEVlFe8v3cgDby1Oy/p1HYyISAvh7izesIPpxRuYVlzCzCWl7CivBODkg/owsEeHRn0+JRgRkRxWun03MxaXMu2zDUxfVMLnW+Lf5mZacYkSjIiIJFZWUcns5ZuYVlzCtOINfLJma9L2+3dtx6ThhYzar1Ojx6IEIyLSjLk7C9ZuY3pxCW8Xb+A/yzZSVlGVsH2nNvkcPrQHE4sKOaqoJ4N6tMfM0hKbEoyISDOzfmsZ04pLmL6ohGnFJZRs352wbas8Y0z/rkwsKmRiUSGj+3Ulv1XT9O9SghERyXI7y/fw3tKNTC8uYXpxCQvXbUvafnBhh2APZVghhw3tQee2rZso0pqUYEREskxVlfPxmi3BXkpxCbOXb6K8MvFhr67tW3Pk0MLwsFch/bq1b8JoE1OCERHJAqs27WR6cQnTFpXwzqISNu2sSNi2dStj/MBuTCzqycSiQg7s24VWeek5j7IvlGBERDJgW1kF7y4uZfqiYC9lScmOpO2H9+7IxKKeHFVUyKGDu9O+IPu/vrM/QhGRHLCnsop5qzZHD3vNXbmZyqrEdxQu7Ngmeh7lqKJCendu24TRNg4lGBGRNHB3lpXujF41/+7iUrbt3pOwfZv8PA4d0oOJYUIZ2adT2roPNxUlGBGRRrJ5ZzkzFpUyfVGQVFZt2pW0/YF9O0fPo4wf2I22rVs1UaRNQwlGRKSByvdUMXv5JqYv2sD04hI+XL0FT3zUi/26tI1e4Hjk0B706Nim6YLNACUYEZF6cneK128Pz6Ns4L2lG9kZFouMp0NBKw4f2iM8j9KToT07NPvDXqlQghERSWLDtt3MCK+Yn75oA+u2Jr5qPs/gkH5dmRTupYwd0JXWTXTVfDZSghERiVFWUcn7SzdGy7DM/zx5scgB3dtHy7AcPqSQLu0zc9V8NlKCEZEWrarKmb92a7T78PvLNlK+J/FV853b5nPE0EImDi9k4rCeDOiRHVfNZyMlGBFpcT7fsiuaUGYsKqF0R3nCtvl5xrgB3aJlWA7p1zUrr5rPRkowIpLzduzew8wlpdEKxIvWb0/afmjPDtHuw4cO6UHHNvqqbAj910Qk51RWOR+t3sK0zzYwbVEJc1dsoqIycf/h7h0KOHJYYfTK+b5d2zVhtLlLCUZEcsLKjTujd3F8Z3EpW3YlLhZZkJ/HFwYFxSKPGlbIAft1Jk+HvRqdEoyINEtbdgXFIqcVB/eaX166M2n7kX06MWl4kFC+MKg77Qpy66r5bKQEIyLNQkVlFR+s3BzdS5m3cjNJakXSq1MbjioqZFJRT44cVkjPTrl91Xw2UoIRkazk7iwp2RHcI6V4AzOXbGR7kmKR7Vq34tAh3aMn54t6dWxRV81nIyUYEckaG3eUMyO8P8q04g2s2VKWsK0ZHLx/F44aVsjEop6MG9iVNvk67JVNlGBEJGN276lk9rJNTFsUJJRP1mxNWixy/67twqvme3LE0B5061DQdMFKypRgRKTJuDsL121jenEJbxeX8P7SUsoqEl8136lNPocN7RGt7TWoR3sd9mpGsi7BmFk74GrgbKAIKADWAbOAu919Rkzbh4FvJFndQncfmb5oRaQu67eWRW8LPG1RCRu2JS4W2SrPGNO/K0cNK2TS8EJG9+tKfgsuFtncZVWCMbPBwMvAMGA98BawGxgEfBWYB8yIs+gMYFGc6Z+nJVARSWhXeSXvLS0Nz6OUsHDdtqTtBxd2CM+jFHLY0B50bqtikbkiaxKMmXUAXgGGArcBt7l7Rcz8HkCPBIv/yd0fTnuQIrKXqirnkzVbmRbedGvWsk2UVyY+7NW1fWuOHBrU9TpqWCH9u6tYZK7KmgQD/IwguTzi7rfUnunupUBpk0clIntZvXlX9F7zMxaVsGln4qvmW7cyxg/sFu0+fGDfLioW2UJkRYIxswLgsnB0ciZjEZG9bSurYOaSjdGksqRkR9L2w3t35KhhPZk4vJBDB3enfUFWfNVIE8uWrT6e4PDXSnefb2ZHAKeF09YCL7r7u0mWP9bMDgE6EnQImA684u6J99NFJKlF67fxrw/XMn3RBuau2MyeJJfNF3Zsw1HDenBUWNurT5e2TRipZKtsSTAHh8PiBD3DbjGzqcBF7r4rzvIXx5n2qZmd5+4fNWKcIi3Cix+v5aon5iRMKm3y8/ji4O7Ra1JG9umk7sOyl2xJMN3D4SSgFXAH8AeCcy6TgN8BZwJbgUtjlvsAmA28BiwHOgPjgF8Co4FXzWycu6+O96RmdjlwOcCAAQMa9xWJNFMrN+7kxqfn7ZVcDuzbmaOKgrs4ThjUjbatddW8JGee7LLZpgrC7CcESQGCHmGX1Zo/AXgfcKDI3ZfUsb4Cgi7OhwH3u/tVdcUwYcIEnzVrVkPCF8kZ5XuqOPuBd5m3cjMQXDn/w5NHcOSwQgo7qlik7M3MZrv7hHjzsuUKptiO8g/Wnunuswj2VPKAY+pambuXA78KR7/cCPGJtAh3vLwwmlzy84wpF4zlq2P2V3KRBsmWBLMs5u+lCdpEpvep5zoXhMP9GxKQSEvzxsL1/PHt6oMDN540grEDumUwImnusiXBzIn5O9HFlIXhMPnNtPdeT33bi7RYa7eUcf3f5kXHjxnRk8smDslgRJILsiLBhCfh3wtHj68938y6EZy8h6AmWX2cEw7/s2/RieS2yirn2qfmsnFHOQC9O7fhzrNH6xbCss+yIsGEIif5bzGzMZGJZtYW+D3QheA8zLvh9DFmdpqZ1ejKYmb5ZvYD4Jpw0l1pj1ykGbv3tWLeW7oRgDyDu88dSw+dc5FGkC3dlHH3f5jZHcANwHtm9h5BN+UvAn2B1cD5Xt3tbRDwLLDRzD4DVgGdCK6p6QtUATe5+0tN+kJEmpF3F5dy3+vF0fGrjyvi8KGJjlKLpCZrEgyAu99oZu8QlOsfC7QHVgC/BSa7+4aY5vOAewgS0MCwvRMkmr8QdE+e3YThizQrpdt3c+1Tc6P3tT90cHeuOb4os0FJTsmqBAPg7s8S7JnU1W4pcF36IxLJPVVVzvVPz2N9eG+W7h0KuOe8sSpCKY0qm87BiEgTeXDaEt5cWH1A4M5zRqt+mDQ6JRiRFmbOik385qWF0fHLJw3h2BG9MhiR5ColGJEWZMuuCq55cm60ztjo/l254cQRGY5KcpUSjEgL4e78aOqHrNoUFCTv1DafKeePpSBfXwOSHnpnibQQj723ghc+Xhsdv/3MQ3S7YkkrJRiRFuCTNVu47Z+fRscvPGwAXz54vwxGJC2BEoxIjtuxew9XPzGX8j3BDV5H9unEz049IMNRSUugBCOS425+/mOWlOwAoF3rVky5YJxuFiZNQglGJIdNnb2Kv8+pvqHrbacfxLBeHTMYkbQkSjAiOWrxhu3c/PzH0fGvjdufs8b3y2BE0tIowYjkoLKKSq58fA47yysBGNKzA7d99aAMRyUtjRKMSA76xb8+ZcHa4E7kBfl5TDl/HB3aZF3pQclxSjAiOeaFjz7nsZkrouM3nzqKA/p2zmBE0lIpwYjkkJUbd/LDqR9Gx085qA8XHjYwgxFJS6YEI5IjKiqruPrJuWwr2wNAv27tmHzmIZipBL9khhKMSI6446WFfLByMwD5ecZ954+lS7vWGY5KWjIlGJEc8MbC9Tzw9pLo+I0njWDsgG4ZjEhECUak2Vu7pYzr/zYvOn7MiJ5cNnFIBiMSCSjBiDRjlVXOdX+dy8Yd5QD06tSGO88eTZ5ufSxZQAlGpBm77/ViZi7ZCECewT3njaVHxzYZjkokkDDBmNkcM/vvBPMGmFn39IUlInV5d3Ep975WHB2/+rgiDh/aI4MRidSUbA9mDJCoA/1S4DeNH46I1Efp9t1c+9Rcwjsfc+jg7lxzfFFmgxKppaGHyCx8iEgTq6pyrn96Huu37Qage4cC7jlvLK103kWyjM7BiDQzf5q+hDcXboiO33n2aPp0aZvBiETiU4IRaUbmrtjEr19cGB2/fNIQjh3ZK4MRiSSmBCPSTGzZVcHVT85lT3jiZXT/rtxw4ogMRyWSmBKMSDPg7vxo6oes2rQLgE5t85ly/lgK8vURluxV1w0iTjaz1xswz939+H2IS0RiPPbeCl74eG10/PYzD6F/9/YZjEikbnUlmD7hI9V53uCIRKSGT9ds5bZ/fhod//qhA/jywftlMCKR+kmWYL7ZZFGISFw7du/hqifnUL6nCoCRfTpx82kHZDgqkfpJmGDc/X+bMhAR2dstz3/Ckg07AGjXuhVTLhhH29atMhyVSP3oDKFIlpo6exVT56yKjt92+kEM69UxgxGJpEYJRiQLLd6wnZuf/zg6/rWx+3PW+H4ZjEgkdUowIlmmrKKSKx+fw87ySgCGFHbgttMPynBUIqlTghHJMr/813wWrN0GQEF+HvddMJYOberq8CmSfZRgRLLICx99zqMzl0fHbz51FAf27ZLBiEQaTglGJEus3LiTH079MDp+ykF9uPCwRHfMEMl+SjAiWaCisoqrn5zLtrI9AOzftR2TzzwEM5Xgl+ZLCUYkC9zx0kI+WLkZgPw8474LxtKlXesMRyWybxp85tDM9gP6EZSFWePuaxotKpEW5I2F63ng7SXR8RtPGsG4Ad0yGJFI40h5D8bMvmNmC4FVwEzgPWClmS00s+81doAiuWzd1jKu/9u86PgxI3py2cQhGYxIpPHUO8GYWSszmwrcDxQR7rmEDw+n3Wdmz5qZalmI1KGyyrn2qbls3FEOQK9Obbjz7NHk6dbHkiNS2YO5FjiDIKFcCrRz9/7u3h9oR1AcczXwlbCtiCRx3+vFzFyyEQAzuPu8MfTo2CbDUYk0nlQSzKVAGXCMuz/s7hWRGe5eERbHPBbYDXyrccMUyS0zl5Ry72vF0fGrjyviiKGFGYxIpPGlkmCGAm+4++JEDcJ5rwM6iCySQOn23Vz71FzCOx9z6ODuXHt8UWaDEkmDVBLMFmBrPdptD9uKSC1VVc71T89j3dbdAHTvUMA9542llc67SA5KJcG8ChxlZgWJGoTzjiTYixGRWv40fQlvLtwQHb/z7NH06dI2gxGJpE8qCeZnBCfzHzWzvQ4Wm1kP4BGgLfCTxglPJHfMXbGJX7+4MDp+2cTBHDuyVwYjEkmvVC60vBj4Zzg81cxeBpaG8wYBJwLtgUeBi2uVuHB3v22foxVpprbsquDqJ+eyJzzxMrp/V248aWSGoxJJr1QSzK0E17tAkEhOT9DuIiCSXTz82wElGGmR3J0f//1DVm3aBUCntvlMOX8sBfmq1CS5LZUE83OqE4yI1NPj763g3x+tjY7ffuYh9O/ePoMRiTSNeicYd781jXGI5KRP12zl5//8NDr+9UMH8OWD98tgRCJNR/voImmyY/cernpyDuV7qgAY2acTN592QIajEmk6SjAiaXLL85+wZMMOANq1bsWUC8bRtrXK9EnLkXK5fjM7kqAkTF+CLsnxuLurXIy0WH+fs4qpc1ZFx287/SCG9eqYwYhEml69E4yZdQSeJuiODNU9xeJxVI9MWqjFG7bzs+c+jo5/bez+nDW+XwYjEsmMVPZgJgMnARuBx4BigrIwIhIqq6jkysfnsLO8EoAhhR247fSDMhyVSGakkmDOBDYBY9x9VV2NRVqiX/5rPgvWbgOgID+P+y4YS4c2Db5xrEizlspJ/s7A20ouIvG98NHnPDpzeXT8Z6eO4sC+XTIYkUhmpZJgFtOATgEiLcHKjTv54dQPo+MnH9iHiw4bmMGIRDIvlQTzF+AYM+uTrmBEmqOKyiqufnIu28r2ALB/13bcftYh1KrHJ9LipJJg7gZeBF43s+PNTNfQiAB3vLSQD1ZuBiA/z7jvgrF0adc6w1GJZF4qpWLczK4A3gReBirMbC1QlaD50MYJUSR7vbFwPQ+8vSQ6fsNJIxg3oFsGIxLJHqlcBzMIeBvYn+AamAJgQILmKoopOW/d1jKu/9u86Pik4T25fKLuFi4SkcpJ+98A/YC3CA6XLULXwUgLVVnlXPfUB2zcUQ5Ar05t+O05o8nTrY9FolJJMMcS3GDsJHcvT1M8mFk74GrgbKCIYE9pHTALuNvdZ8RZ5gLgu8AhQCtgAUGnhN+7e7xDeCL7ZMrri3h3SSkAZnD3eWMo7Ngmw1GJZJdUEkxr4D9pTi6DCc7vDAPWE+wt7Sa4Y+ZXgXnAjFrL3A98DygDXgMqgOOBKcDxZna2u1emK2ZpeWYuKeWe1z6Ljl99XBFHDN3rLuIiLV4qCWYe0DtdgZhZB+AVYCjB3S9vc/eKmPk9gN1waHgAABm+SURBVB61ljmTILmsBSa5e3E4vTfwBnAGcBVwT7rilpaldPturn1qLuGdj/ni4O5cc9ywzAYlkqVS6Wr8G2CimR2eplh+RpBcHnH3W2KTC4C7l7r7Z7WW+XE4vCmSXMK26wgOmQH8SF2qpTFUVTk3PD2PdVt3A9CtfWvuPW8s+a309hKJJ9U9mDuBV8zsLuAlYBXxuynj7ivqu2IzKwAuC0cn13OZfsB4oJygynPt53/LzFYT9Ho7DHinvvGIxPPQ9KW8sXBDdPzOc0bTp0uiO1aISCoJZmk4NOAn4SMRT3Hd4wkOf6109/lmdgRwWjhtLfCiu79ba5mx4fATd9+VYL3/IUgwY1GCkX3wwcrN3P7iguj4ZRMHc9zItB0xFskJqSSBlaTv+paDw2GxmT0MfKPW/FvMbCpwUUwyGRwOl5NYZC9qcJI2Ikm9s6iEG5/5kD3hiZfR/bty40kjMxyVSPZL5Ur+QWmMo3s4nETQzfgO4A9AaTjtdwS3C9gKXBq2jdwecEeS9Uau0+nUmMFKy/Dpmq1MfnEBb39WfVisU5t87jtvLAX5Ou8iUpdsqY4c+bTmA39y9xtj5v2fma0B3ge+YWa/cPclVN9Rs8F7VWZ2OXA5wIABiYoSSEuzcuNOfvvKZzz3wWo85t3VJj+Pu84dw4Ae7TMXnEgzki0JZlvM3w/Wnunus8xsNjABOAZYErNMshudR+ZtizfT3f8I/BFgwoQJKm/Twm3aUc79byzikXeXU15Z3Xclz+DMcf34/gnD6du1XQYjFGleUk4wZnY0wbUlhwM9gcfc/VvhvJOAo4F73X1tCqtdFvP30gRtlhIkmMjtAiLLJLvpRv846xepYVd5JX95Zym/f3NxtOR+xPEje/HDk0cyoo+OsoqkKqUEY2a3AjdTfXiKWn9XADcBq4H7U1j1nJi/ewAb4rSJXCodOa8yNxweaGbtEvQk+0KttiJReyqrmDpnFXe9UszarWU15o3p35UfnzKSQ4f0SLC0iNQllWrK/wXcQtCb7AcElZXX1Wr2BsGJ+dNIIcG4+2ozew84lKDMy4LY+WbWDRgXjs4Kl1lpZnPC6WcDj9Ra5miC4pxrgdpdnKUFc3denb+eX7+4gOL1Neu1DinswA9PHsFJB/bRDcNE9lEqezDXENQFO9nd5wN7fQDDe8Z8RlBLLFW/BP6PoEvyDHf/IHyOtsDvgS7AbGomi18RXGR5u5m94+6LwmV6EfQ8A5isgpcSMXv5Jia/MJ//LNtUY3rPTm247ktFnDOhP611Zb5Io0j1YsiZkeSSxCpgdKqBuPs/zOwO4AbgvXCPphT4ItCX4LDb+e7V/Xrc/Rkz+z1BWZiPzOxVqotddgaeIyh6KS3c4g3b+fWLC3jpk5o73R3b5HPFpCF8a+Jg2hdkS58XkdyQyieqHfHPjdTW4LOh7n6jmb1DUK5/LNCe4GLJ3xLsiez1/O7+PTObDlxJ0MEgUq7/z6hcf4u3fmsZd71azN9mraSyqrqjYOtWxtcPHcjVxw2jh8rsi6RFKgnmc6A+ly8fQPKr65Ny92eBZ1Nc5gngiYY+p+SebWUVPPDWEh6avpRdFTXv1vCV0X25/sThDOzRIUPRibQMqSSYN4BLzOxEd385XgMzO5eg27DK40tG7N5TyeMzVzDljUXRu01GHDmsBz86eRQH9+uSoehEWpZUb5n8deBpM7sRmBqZYWbtgbOAe4Gd4VCkyVRVOf/4cA13vLyQlRtr9lg/YL/O/OiUkUwsKlTPMJEmlEotsgVmdgnwMEGvrt8RlGm5kOrilHsIClImulhSpNFNLy5h8ovz+Xj11hrT+3Vrxw0njuAro/uSl6fEItLUUuo24+5PmdknBDcHO4mgp1Y+sAt4Ffi5u89u9ChF4vh49RZuf3EB04pLakzv1r41Vx1XxIWHDaBNfqsMRSciKffLdPePgHMtONbQg6DXVonuey9NZeXGndzx8kKe/2BNjeltW+fxraMGc8XRQ+nctnWGohORiFSu5B8AbHf3jRBcVAmUxGnXDeiUyh0tRepj445y7nu9mMdmLqeisrrLcZ7BORP6c92XhusOkyJZJNU7Wj4MfKuOdr8GvpniukUS2lm+hz9PX8oDby1h2+6axShPOKA3N508gmG9VIxSJNukkgSMmoUt62orsk/2VFbx9OxV3PXKZ6zftrvGvPEDu/HjU0YyYVD3BEuLSKalYy+jK0HNMpEGcXde/nQdv35xAYs31Lxh6dCeHbjp5JGccEBvdTkWyXJJE0x43iVWxzjTYtc1CjiRxPd0EUlq1rKN/OqFBcxeXrMYZe/Obfj+l4Zz1vh+5KsYpUizUNcezDJq3pL4zPCRjAGP70NM0gIVr9vG7S8u5NX5NYtRdmqTz3eOGcqlRw6mXYG6HIs0J3UlmBVUJ5gBBFfp79VzLFROUPH4WVTBWOpp7ZYy7nrlM56evZKYWpS0bmVcdNggrjpuGN07FGQuQBFpsKQJxt0HRf42syrgaXe/NN1BSe7bsquCB95azJ9nLKWsombB69PH9OX6E0fQv3v7DEUnIo0hlZP83wQWpSsQaRl276nk0XeXM+WNRWzeWVFj3sSiQm46eSQH7a9ilCK5IJVaZP+bzkAkt1VVOc/PW80dL33G6s01i1EetH9nfnTyKI4qKsxQdCKSDroYUtLK3Xm7uITJLyxg/uc1i1H27x4Uo/yvQ1SMUiQXKcFI2ny0agu/emE+7ywurTG9e4cCrj5uGF8/dCAF+epyLJKrlGCk0S0v3cEdL3/GP+bVLEbZrnUrvj1xMJdPGkInFaMUyXlKMNJoSrbvZsrri3j8vZrFKFvlGed+oT/XHV9Er84qRinSUijByD7bsXsPD01fygNvLWZHec27Npx8YB9uPHkEQ3t2zFB0IpIpSjDSYBWVVfz1Pyu5+9ViSrbXLD/3hUHd+NEpoxg/sFuGohORTFOCkZS5Oy9+vJbfvLSQJSU1i1EW9erITSeP5PhRvVSMUqSFU4KRlKzfWsbVT87lvaUba0zv07ktPzhhOGeO70crdTkWEZRgJAXzVm7mikdns3ZrWXRap7b5fO+YYXzzyEG0ba1ilCJSTQlG6mXq7FX8+NmPKN8T1A3LM7j0yMFceewwuqkYpYjEoQQjSe2prGLyCwv40/TqW/x0bpvPlAvGMWl4zwxGJiLZTglGEtq8s5yrn5zLtOLqOzQU9erIgxdPYFBhhwxGJiLNgRKMxPXZum1c9sgslpfujE474YDe3HXuGDq20dtGROqmbwrZy0ufrOUHf/2gxkWT1xxfxHXHF6kopYjUmxKMRFVVOfe9voi7Xv0sOq19QSvuPHs0pxy8XwYjE5HmSAlGgKDcy/V/m8eLn6yNTuvXrR0PXjyBUft1zmBkItJcKcEIK0p3ctkjs1i4blt02hFDezDlgnF0VxdkEWkgJZgWbsaiEq58Yk6N2xd/88hB/PTLo8hvpXu1iEjDKcG0UO7OX2Ys45f/nk9lVVBav6BVHr844yDOmdA/w9GJSC5QgmmByioq+dlzH/PM7FXRab06teEPF41n3ABVPxaRxqEE08Ks21rGFY/O5oOVm6PTRvfvyh8vGk9v3QxMRBqREkwLMnfFJq54dDbrt1Xfu+XMcf345RkHqVCliDQ6JZgW4ulZK/npsx9TXhkUq2yVZ/zky6O49MhBum+LiKSFEkyO21NZxS//PZ+/zFgWndalXWvuv2AcRxUVZi4wEcl5SjA5bNOOcq56cg4zFpVGpw3vHRSrHNhDxSpFJL2UYHLUgrVbueyRWazcuCs67aQDe3PnOSpWKSJNQ980OejFjz/nB3+bx86YYpXf/9Jwrj5umIpVikiTUYLJIVVVzt2vFXPva8XRaR0KWvHbc8dw0oF9MhiZiLRESjA5YvvuPXz/rx/wyqfrotMGdG/PgxdPYESfThmMTERaKiWYHLC8dAeXPTKLz9Ztj047alghUy4YS9f2KlYpIpmhBNPMTSvewFVPzGXLrupild86ajA/PmWkilWKSEYpwTRT7s5D05fyP/+eT1irkoL8PH51xsGcOb5fZoMTEUEJplkqq6jkJ89+xN/nrI5O6925DQ9cNIEx/btmMDIRkWpKMM3M2i1lXPHoLOat2hKdNnZAVx64cDy9VKxSRLKIEkwzMnv5Jr7z2Gw2xBSrPGdCP247/SDa5KtYpYhkFyWYZuJv/1nJz56rWazyltMO4OLDB6pYpYhkJSWYLFdRWcUv/vkp//vu8ui0bu2DYpVHDFOxShHJXkowWWzjjnK+9/hsZi7ZGJ02sk8nHrx4Av27t89gZCIidVOCyVKfrtnK5Y/OYtWm6mKVpxzUhzvOHk0HFasUkWZA31RZ6F8ffs4NT89jV0V1scrrTxjOVccN0/kWEWk2lGCySFWVc9ern3Hf64ui0zq2yeeuc8dwwgG9MxiZiEjqlGCyxLayCr7/1w94df766LRBPYJilUW9VaxSRJofJZgssLQkKFa5aH11scqJRYVMOX8cXdq3zmBkIiINpwSTYW99toGrn5jD1rI90WmXTxrCD08aoWKVItKsKcFkiLvz4LQlTH5hQY1ilbefeTBnjFWxShFp/pRgMqCsopIfTf2Q5z5YE53Wp3Nb/njxeA7pp2KVIpIblGCa2JrNu7ji0dl8tLq6WOX4gd34/YXj6NVJxSpFJHcowTShWcs28p3HZlOyvTw67bwv9Of/ffVAFasUkZyjBNNEnnx/Bbc8/zEVlcEJl/w847//6wAuPEzFKkUkNynBpFlFZRU//8enPDqzulhl9w4F3H/BOA4f2iODkYmIpJcSTBqVbt/Ndx+fw/tLq4tVjtqvMw9ePJ5+3VSsUkRyW9ZcaGFmD5uZJ3ksaIxlmsona7bwlSkzaiSXUw/Zj6nfPVzJRURahGzcg5kBLIoz/fNGXiZt/jFvDTc+M4+yiuDmYGZww4kj+N4xQ3W+RURajGxMMH9y94ebYJlGV1nl3PnyQn735uLotI5t8rnnvDEcP0rFKkWkZcnGBNMsbS2r4LqnPuD1BdXFKgcXduDBi8czrJeKVYpIy6ME0wgWb9jOZY/MYsmGHdFpRw/vyb3nj6VLOxWrFJGWKRsTzLFmdgjQEVgHTAdecfeqRl6mUbyxcD3XPDmXbTHFKr9z9FBuPGkErfJ0vkVEWq5sTDAXx5n2qZmd5+4fNeIy++yZ2au48Zl5eFissk1+Hr8+6xC+Omb/dD2liEizkTXdlIEPgGuAAwn2RPoCpwHzgAOAV82s9jd3Q5aJMrPLzWyWmc3asGFDygEfOrg73doXANC3S1ue+c4RSi4iIiHzyM/vLGVmBcBbwGHA/e5+VTqWmTBhgs+aNSvl+N5dXMq9rxVz3wVjKezYJuXlRUSaMzOb7e4T4s3LxkNkNbh7uZn9Cnge+HK6lmmow4f24LAh3XV9i4hILdl0iCyZyBX5qRx/asgyDaLkIiKyt+aSYCJVIbcnbbXvy4iISCNpLgnmnHD4nzQvIyIijSQrEoyZjTGz08ysVa3p+Wb2A4KeYgB37csyIiLSdLLlJP8g4Flgo5l9BqwCOgEHE3Q9rgJucveX9nEZERFpItmSYOYB9wBfBAYCYwEnSBp/IehqPLsRlhERkSaSFQnG3ZcC16V7GRERaTpZf6FlUzGzDcDyOhvGVwiUNGI4kp20nVsGbefUDHT3nvFmKME0AjOblehKVskd2s4tg7Zz48mKXmQiIpJ7lGBERCQtlGAaxx8zHYA0CW3nlkHbuZHoHIyIiKSF9mBERCQtlGBERCQtlGAayMz6mdl9ZrbQzHaZWZmZFZvZH8xsSKbjk/oxsxFmdq2ZPWZmC8ysyszczM6qx7IXmNk0M9tiZtvDu6NeaWb6XGWRVLexmbU2s+PN7E4zm2lmn5tZuZmtNrNnzOyYJn4JzZbOwTSAmY0FXge6EpSmiZSkmUBw/5ntwEnu/k5mIpT6MrO7gWvjzDrb3Z9Jstz9wPeAMuA1oAI4nqAe3rPh8pWNH7GkKtVtbGZfAl4JR9cSfL53ENyG/aBw+m3ufksaws0p+qXVMPcTJJcHgSHufrq7nw4MBv4MdAR+n8H4pP4+Bn4DnAsMI7jVdlJmdiZBclkLHOLup7n7GUARMB84A6jzNt3SZFLdxlXAVGCSu+8Xbt9z3f1g4DygErjZzI5NZ9C5QHswKTKztsCucHQ/d19ba35fYHU42sHddzZlfLJvzOxN4GiS7MGY2SxgPPANd3+k1ryjgTcJks/+7l6V1oAlZfXZxnUs/yfgW8Cf3f1bjRxeTtEeTOoqgT3h3/HulRzJ2DuoTkSSI8ysH0FyKQeerj3f3d8i+IHRBzisaaOTJjI3HPbLaBTNgBJMity9guCYO8D/M7PWkXnh378IRx9y7R7morHh8BN3T/QD4j+12kpuKQqHn2c0imYgK8r1N0PfA14ELgNOCQ+ZAHwB6EZwn5obMxSbpNfgcJis8vaKWm0lR5hZH+CScHRqBkNpFpRgGsDdl5jZEcAjwCnU3FWeBbwd7ulI7ukYDnckabM9HHZKcyzShMwsH3gM6AK85u7/yHBIWU+HyBogTC4fE/RI+SrB/SN6AqcT7MFMNTN1YcxNkfNuOvzZ8vyBoCv6SuDCDMfSLGgPJkVm1hV4DugAHOHuS2JmP29mnwAfEnRjfNLdizMRp6TNtnDYMUmbyLxtSdpIM2Jm9xD0HFsLHF+796jEpz2Y1J1KsLcys1ZyAcDdFwHvESTvY5o2NGkCy8LhwCRt+tdqK82Ymd0JXANsIEgu+tFYT0owqRsQDrckabM5HHZPcyzS9CJdVA80s3YJ2nyhVltppszs18APgFLgBHf/NMMhNStKMKlbEw7Hx3ZRjginjQ9HlzZZVNIk3H0lMAcoAM6uPT+80LIfwaGUd5s2OmlMZjaZoDfoJoLkMi/DITU7SjCpewHYSbAnc5eZtYnMCP++l+AQySbgpYxEKOn2q3B4u5kNi0w0s17A78LRybqKv/kys9uAmwiORpzg7tobbQCVimkAM/sG8BDQimCPZjZB76LxwH7AbuA8d38uY0FKvZjZOKqTAgQFDTsBxcDGyER3P6zWcr8DvktQ7PJVqotddiboBHKWil1mh1S3sZl9BXg+nDwL+CTBqhe4++RGDziHKME0UPimvQ6YSJBUICgR8gbwWx2rbR7C0utv1NXO3fcqC2RmFwBXAgcT/NhYQFDs9Pfae8keqW5jM7sE+Es9Vv2Wux+zL7HlOiUYERFJC52DERGRtFCCERGRtFCCERGRtFCCERGRtFCCERGRtFCCERGRtFCCERGRtFCCEYnDzI4xMzezNzMdS0M09/glNyjBSItkZsvCL+BBmY5FJFfphmMi8b0PjCIobCoiDaAEIxKHu+8kqC0mIg2kQ2TSopjZJWbmVN+Rcml4qMxjD5klOodhZoPC6cvMLM/MfmBmn5jZLjNbZWa/NbP2YdtuZnZ32Ha3mRWb2Q+SxGZmdp6ZvWxmJeEyK8zswX05lGdmrc3sp2a2wMzKzGy9mT1mZgPitI2+bjNrb2a/CJfbZWYfxLT7opk9bWarzazCzLaY2SIze8LMjmtorJJbtAcjLc0i4H+Bs4AOwFRge8z87fEWSuAJ4DTgzXC9k4DvA6PM7OvATIKy8NMJ7m46CbjTzNq6+//Erii8Ud1TwNeAXQRl4tcBBwHfBs40sxPdfVYqLxZoTXAPo0OBt4D5wOHA14FJZnaIu2+Os1zb8HWNAt4G5hHcZA0zOwH4V7juD4AZ4d/9CP6vW4HXU4xTcpG766FHi3sAywAHBiWYf0w4/81a0weF053gEFrfmHn9gZJw3kfA00DbmPmnhvO2Au1rrXdyOO8toF+teVeF8xYB+fV8fcfExPkfoFfMvC4E9zBy4KdJlpsL9I6z7tfD+efHmdcDGJ/p7atHdjx0iEyk4a5x98gttPHgdsqPhaMDge+6e1nM/H8BHxLs1UyITDez7sA1BHtPZ7v7qtgncfcpBHsMQ4FTUozRgUvdfX3M+rYAt4ejxydZ9kp3Xxdneu9w+MJeT+Ze6u6zU4xRcpQSjEjDVBD/MNCicDjL3UvizC8Oh31jph0LtCO4gdX6vRcBgj0bCA5vpWKFu38UZ3qkA0PfOPMA1rn7OwnmvR8OnzCzI82sVYoxSQuhczAiDbPW3ffEmR45h7MqzrzY+W1jpg0Jh6eGHRCS6VnP+CJWJJi+NU4csZYnWeePgTEEe1OnADvMbDZBwn3U3ZekGKPkKCUYkYap65bIqdwyObIHsJCgY0Ay76Ww3lTjiLUr0Qx3X2tm4wnO15wAHEnQiWAS8DMzu8Ld/9zA55UcogQjknkrw+FH7n5JJgOpL3evIthjeR3AzDoQdEaYDNxvZs+4+9Ykq5AWQOdgpKUqD4fZ8CPrVYJzOl8ys66ZDqYh3H2Hu99OcGiwLTAiwyFJFlCCkZZqdTgcldEogLCn1v1AV+D/zGxk7TbhRZvfNrPee62giZnZDWbWP870CcB+BIflEp2DkhYkG369iWTCswTnEB43s5eByMWGN7l7aQbi+SFBj65zgI/Dq+aXEuwN9CdIhAXhMF7X4ab0M+A3Zjaf4MLN3QQxHkHwo3Wyu3+ewfgkSyjBSEs1BehMcEX7aUCbcPovgCZPMO5eAZxrZo8DlwJfBA4BtgGfE1QNeB5Y3NSxxXElwcn9CVR3sf4c+AfwO3d/OYOxSRYx97p6RYqIiKRO52BERCQtlGBERCQtlGBERCQtlGBERCQtlGBERCQtlGBERCQtlGBERCQtlGBERCQtlGBERCQt/j//ji2I4U9y+AAAAABJRU5ErkJggg==\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"plt.plot(time,Tempa)\n", | |
"plt.title('Ambient temp vs time')\n", | |
"plt.xlabel('time hrs')\n", | |
"plt.ylabel('temp F')" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 212, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"This looks fairly correct, but the lines are jagged due to the linear interpolation\n", | |
"A better way to get Ta would be to would be using the analytical method above to find the exact solution\n" | |
] | |
} | |
], | |
"source": [ | |
"print('This looks fairly correct, but the lines are jagged due to the linear interpolation')\n", | |
"print('A better way to get Ta would be to would be using the analytical method above to find the exact solution')" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 216, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def temperature2(N,t):\n", | |
" ''' \n", | |
" computes the analytical and numerical solution for the temperature of a corpse\n", | |
" N is number of timesteps between 0 and 2 sec\n", | |
" temp_analytical is the 64-bit floating point \"true\" solution\n", | |
" temp_numerical is the 32-bit approximation of the velocity\n", | |
" this functions call tempattime within it, letting the ambient temp change depedning on time'''\n", | |
" \n", | |
" T0=85\n", | |
" Ta=tempattime(t)\n", | |
" t1=np.linspace(0,N,N+1)\n", | |
" temp_analytical = Ta+((T0-Ta)*np.exp(-k*t1))\n", | |
" temp_numerical=np.empty(len(t1),dtype=np.float32)\n", | |
" delta_time =np.diff(t1)\n", | |
" temp_numerical[0]=85\n", | |
" for i in range(0,N-1):\n", | |
" temp_numerical[i+1]=temp_numerical[i]+(k*(Ta-temp_numerical[i]))*delta_time[i]\n", | |
" \n", | |
" return temp_analytical,temp_numerical,t1\n", | |
" " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 217, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"(array([85. , 78.50844733, 73.57764488, 69.8323448 , 66.98751926,\n", | |
" 64.82666909, 63.18534753, 61.93864544, 60.99168528, 60.27240074,\n", | |
" 59.72605225, 59.31106117, 58.99584552, 58.7564165 , 58.57455288,\n", | |
" 58.43641435, 58.33148818, 58.25178918, 58.19125204, 58.14526972,\n", | |
" 58.11034283]),\n", | |
" array([8.5000000e+01, 7.7574997e+01, 7.2191872e+01, 6.8289108e+01,\n", | |
" 6.5459602e+01, 6.3408211e+01, 6.1920952e+01, 6.0842690e+01,\n", | |
" 6.0060951e+01, 5.9494190e+01, 5.9083286e+01, 5.8785381e+01,\n", | |
" 5.8569401e+01, 5.8412815e+01, 5.8299290e+01, 5.8216984e+01,\n", | |
" 5.8157314e+01, 5.8114052e+01, 5.8082687e+01, 5.8059948e+01,\n", | |
" 1.3452465e-43], dtype=float32),\n", | |
" array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,\n", | |
" 13., 14., 15., 16., 17., 18., 19., 20.]))" | |
] | |
}, | |
"execution_count": 217, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"temperature2(20,9)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 220, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"(array([85. , 80.19144246, 76.53899621, 73.76469985, 71.65742167,\n", | |
" 70.05679192, 68.84099817, 67.91751514, 67.21606317, 66.68325981,\n", | |
" 66.27855722, 65.97115643, 65.73766335, 65.56030852, 65.42559473,\n", | |
" 65.32326989, 65.2455468 , 65.1865105 , 65.14166818, 65.1076072 ,\n", | |
" 65.08173543]),\n", | |
" array([8.5000000e+01, 7.9500000e+01, 7.5512497e+01, 7.2621559e+01,\n", | |
" 7.0525627e+01, 6.9006081e+01, 6.7904411e+01, 6.7105698e+01,\n", | |
" 6.6526634e+01, 6.6106812e+01, 6.5802437e+01, 6.5581764e+01,\n", | |
" 6.5421776e+01, 6.5305786e+01, 6.5221695e+01, 6.5160728e+01,\n", | |
" 6.5116531e+01, 6.5084488e+01, 6.5061256e+01, 6.5044411e+01,\n", | |
" 1.3452465e-43], dtype=float32),\n", | |
" array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,\n", | |
" 13., 14., 15., 16., 17., 18., 19., 20.]))" | |
] | |
}, | |
"execution_count": 220, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"temperature(20)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 227, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"making Ta a function of time changes the results of the euler approximation by altering the smoothness of the exponential decay. It is know two functions added together, which makes the temperature decay less smooth\n" | |
] | |
} | |
], | |
"source": [ | |
"print('making Ta a function of time changes the results of the euler approximation by altering the smoothness of the exponential decay. It is know two functions added together, which makes the temperature decay less smooth')" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 228, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"1.4533333333333331\n", | |
"Time of death is 1.3594900679739386 hrs before T0\n" | |
] | |
} | |
], | |
"source": [ | |
"Ta=tempattime(8)\n", | |
"T0=85\n", | |
"dT=(98.6-Ta)/(T0-Ta)\n", | |
"print(dT)\n", | |
"tdeath=-(np.log(dT))/k\n", | |
"\n", | |
"print('Time of death is',-(tdeath),'hrs before T0')" | |
] | |
}, | |
{ | |
"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 | |
} |