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.
435 lines (435 sloc)
60.6 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": 585, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "1. K = 0.275 1/hrs\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "delta_T = 74-85\n", | |
| "delta_t = 2\n", | |
| "rate = delta_T/delta_t\n", | |
| "k = -1*rate/(85-65)\n", | |
| "print('1. K =',k, '1/hrs')" | |
| ] | |
| }, | |
| { | |
| "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": 584, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "2. K = 0.275 1/hrs\n" | |
| ] | |
| } | |
| ], | |
| "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", | |
| " \n", | |
| " Arguments\n", | |
| " ---------\n", | |
| " Temp_t1: Temperature of corpse when discovered\n", | |
| " Temp_t2: Temperature of corpse after elapsed time\n", | |
| " Temp_ambient: Ambient temperature\n", | |
| " delta_t: Change in time between Temp_t1 and Temp_t2\n", | |
| " \n", | |
| " Returns\n", | |
| " -------\n", | |
| " the emperical constant K\n", | |
| " \n", | |
| " '''\n", | |
| " delta_T = Temp_t2 - Temp_t1\n", | |
| " rate = delta_T/delta_t\n", | |
| " k = -1*rate/(Temp_t1-Temp_ambient)\n", | |
| " return k\n", | |
| "\n", | |
| "print('2. K =', measure_K(85,74,65,2), '1/hrs')" | |
| ] | |
| }, | |
| { | |
| "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": 582, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "3a. Time step of 0.045\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU9fX/8dcJO4jIZmULobbuLCpgq7ZqtZTaVlTU0k61ohC1tS61Km3cNf78Wr4Vqa0WRa3tFGtdu9m6ol93QVFB7KKSiCsgypIAAc7vj3sDk2QmmSRzZzIz7+fjkUdm7jZnhsvJnfNZrrk7IiJSPEpyHYCIiGSXEr+ISJFR4hcRKTJK/CIiRUaJX0SkyCjxi4gUGSV+EWk3M7vdzC7IdRySHiX+ImFm3zWzBWa2zszeN7MHzezgXMeVjJmVhnHW/7iZrU94/qVcx9hWZtY9fD9Dcx1Lusxstxb+Pca5+8nufm2uY5X0dM51ABI9M/sxMAM4HfgnsAmYCEwCnmrlsTq7++aMB5nA3auBHRJe04HR7v7fKF83E7Lx+ZhZJ3ffEuVrJHL3fxP+e5hZd6AW2N3dl2crBsksXfEXODPrA1wB/NDd73X39e5e5+5/cffzw226mdksM3sv/JllZt3CdYea2XIzu9DMPgBuS1j2MzNbaWbLzCyW8JpHmtnrZrbWzN41s58krPummS0ys0/M7BkzG9XG99UjjPMdM/vAzH6ZEPNEM/uvmV0UxvduGNMkM3vTzFaZ2XkJx7rGzOaZ2T1hzC+a2d4J64eZ2QPhsd4ys9Mb7fsHM/ujma0FppjZQWb2vJl9Gn6e15lZ/UXWk+Hvf4VXy0eb2elm9kjCMRt8KzCzO81stpk9ZGbrgS829/4bfU69wtf5XMKyIWZWa2Z9zWwXM/tH+O+xyswea+O/x51mdlEbP/9OZnZx+NmuNLO4me3UljgkPUr8he+LQHfgvma2qQC+AIwBRgPjgYsS1u8C9AOGA+UJywYAQ4DvA3PMbPdw3VzgNHfvDewDPAZgZvsBtwKnAf2B3wB/Tpaw0nAdMBQYCewO7EbwrabecKAujPOa8HWPA0YBRwCVZjYkYfvJwG/D9/kAcG+YkDoBfweeAQYTfFP6mZkdkmTfPsA94eueGR7rS8C3gGnhtl8Of+/u7ju4+/1pvt/vARcDvYEX03j/ALj7euDPwHcSFk8B/unuq4ELgX8R/FsOAi5LM56WtObzPx+YABwcvqe68P1JVNxdPwX8A8SAD1rY5k3gyITnXwOWhY8PJSgNdU9YfyiwGeiVsOwu4OLwcTVBct+x0evcCFzZaNm/gENaiM+BzyU87xzGNCRh2WHA0vDxROBToCR8PjA8xuiE7ZcAE8PH1wDzGx1/FTAOOAT4T6N4LgduTNj3oRbinwHMCx93D2MZmrD+dOCRhOcNtgHuBOak+/6TvP43gdcTni8ETggfXwv8CfhsmudTk/gTYryojZ//28BBCetGADWA5fr/T6H+6Iq/8K0CBiSUGpIZDFQlPK8Kl9Vb4e4bGu2z2oOryWT7TAaOBKrM7Akz+2K4fDhwXlhW+MTMPgGGNXqtdAwGugBLEo5zP7Bzo5i3ho9rw98fJqyvJaEdAXin/oEHNfr3wtcZDpQ1ivnHBFeyTfYFMLO9LGg8/9DM1gCXEFxRt0fia6Tz/hP9E/iMmY02s92AzwN/CddVErzXx8PyzI/bGWe9tD5/MzOCc+DvCe/lZYJqRP8MxSKNqHG38D0LbACOBu5Osc17BAluSfi8NFxWL9kUrn3NrFdC8i8FFgO4+4vAJDPrQlDyuIvgP/c7QKW7V7b97QDwPsE3jl3dfVU7j1VvWP2DsLwzmOAzWAO84e4jm9m38edzMzAfON7d15nZDILyRrJtAdYDPROe75Jkm8T9WvX+3b3OzO4mKPdsBO5z99pw3afA2cDZZjaa4A/A8+7+dEvHzQR3dzN7FzjW3Rdm4zVFNf6CF/7HvgT4VdiQ2NPMupjZ182svvvdPOAiMxtoZgPC7X+fxuEvN7OuFnSv/Cbwp/B5zMz6uHsdQeKs74FyM3C6mR1ggV5m9g0z693K91RHUDO+3swGhMcaZmZfbc1xGjnQgobnLsAFBN+UXiLs9WRm54SNrp3NbFTYXpFKb+DTMOnvDUxPiH0jQRnkswnbLwL2NbO9zawnweefUhvf/x8IavvfCR8Tvq+jzGxEeOX9KcG/VdZ6DIVuAq4xs2FhTDub2beyHENRUeIvAu7+C4LyxEXACoIr7zMJygMAVwELgFeB1wgS3lUtHPYDYDXBVXEcON3d3wjXnQgsC8scpxM0TOLuCwiS4A3hvv8FTm7j2zonfO0FBAnrH8Dnmt2jefcAp4RxTQYmu/uWMMkeCRxIUM5aQdBWsUOqAwHnAtPMbB3wK+CPjdZfQvBH8hMzO8rdXyOotf8f8AbBt4WWtPb9Pwl0ImiAfiRh+Z7h660Nt5np7s+l8fqZdG0Y02MW9Ix6BmjuD6u0k4WNKSJpM7NDgd+7e94MQmqOmV0DDHD3aS1uLFIAdMUvIlJklPhFRIqMSj0iIkVGV/wiIkUmL/rxDxgwwMvKynIdhohIXlm4cOFKdx/YeHleJP6ysjIWLFiQ6zBERPKKmVUlW65Sj4hIkVHiFxEpMkr8IiJFJi9q/CJSGOrq6li+fDkbNjSe7FXao3v37gwdOpQuXbqktb0Sv4hkzfLly+nduzdlZWUE88JJe7k7q1atYvny5YwYMSKtfQq21BOPQ1kZlJQEv+PxXEckIhs2bKB///5K+hlkZvTv379V36IK8oo/HofycqipCZ5XVQXPAWKx1PuJSPSU9DOvtZ9pQV7xV1RsT/r1amqC5SIixa4gE391deuWi0jxMDPOO++8bc9nzpzJZZddltUYFixYwFlnndWmfQ899NB2D2gtyMRfWppi+dBs31hIRDqabt26ce+997Jy5cqcvP7mzZsZO3Yss2fPzsnrQ4Em/spK6Nmz4bKeXddTefx58OnS3AQlIh1C586dKS8v57rrrmuy7uSTT+buu7ffmnqHHYIbrc2fP59DDjmEE044gd12240ZM2YQj8cZP348I0eO5M033wRgxYoVTJ48mXHjxjFu3Diefjq4dfFll11GeXk5EyZM4KSTTmL+/Pl885vfBGDdunVMnTqVkSNHMmrUKO655x4AzjjjDMaOHcvee+/NpZdemtnPIKNH6yDqG3ArKoLyTmkpVM54j9jOd8JDtxP/6Ekqfj5q+7pKNfqKZN3Cc2D1oswes+8Y2H9Wi5v98Ic/ZNSoUVxwwQVpH/qVV15h6dKl9OvXj89+9rNMmzaNF154geuvv55f/vKXzJo1i7PPPptzzz2Xgw8+mOrqar72ta+xdGlwsblw4UKeeuopevTowfz587cd98orr6RPnz689tprAKxevRqAyspK+vXrx5YtWzj88MN59dVXGTVqVCs+jNQKMvFDkMgbJvPPw/oXiF/2a8pn70rNpmCpevyIFJ8dd9yRk046idmzZ9OjR4+09hk3bhyDBg0CYNddd2XChAkAjBw5kscffxyARx55hNdff33bPmvWrGHt2rUAHHXUUUlf65FHHuHOO+/c9rxv374A3HXXXcyZM4fNmzfz/vvv8/rrryvxt0mvUiruupqaTQ0rXPU9fpT4RbIojSvzKJ1zzjnst99+TJ06dduyzp07s3XrViAYGLVp06Zt67p167btcUlJybbnJSUlbN68GYCtW7fy7LPPJk3wvXr1ShqHuzfpjvn2228zc+ZMXnzxRfr27cvJJ5+c0dHOBVnjb071O8nfsnr8iBSXfv36ccIJJzB37txty8rKyli4cCEADzzwAHV1da065oQJE7jhhhu2PV+0qOVSVuN9Vq9ezZo1a+jVqxd9+vThww8/5MEHH2xVHC0pusSfusdP6/6BRST/nXfeeQ1690yfPp0nnniC8ePH8/zzz6e8Sk9l9uzZLFiwgFGjRrHXXntx0003tbjPRRddxOrVq9lnn30YPXo0jz/+OKNHj2bfffdl77335pRTTuGggw5q9XtrTl7cc3fs2LGeqRuxNB7VC9Czaw1zTjuL2E9PgEETMvI6ItLU0qVL2XPPPXMdRkFK9tma2UJ3H9t426K74o/FYM4cGD4czILfc25cT2ziCzD/SOIz/0HZcNccPyJSsIqrcTfUtMfPQKh7hvgVv6H82i9RsyloaFGPHxEpREV3xZ9Slx2ouOPH1GxqWNPTHD8iUmiU+BNUv5N8hjv1+BGRQqLEnyBlj5/BNclXiIjkISX+BEnn+OlWS+WkabD4KuJx181dRCTvKfEnSNrj5+ZOxGKdid/4OuWnbqSqCty3N/wq+Yvkn/vuuw8z44033mjzMRpP6JbM1Vdf3eD5gQce2KbXuuyyy5g5c2ab9k1Gib+RWAyWLYOtW4PfsRO7whd/S8V9N1KzsXuDbdXwKxKtqG6hOm/ePA4++OAGc+REoXHif+aZZyJ9vXQp8afDjOoP+iRdpYZfkWjUD7bM9LfsdevW8fTTTzN37txtiX/+/PkceuihHHfcceyxxx7EYjHqB7deccUVjBs3jn322Yfy8nIaD3p99NFHOeaYY7Y9f/jhhzn22GOZMWMGtbW1jBkzhljYH7x+mmeAa6+9lpEjRzJ69GhmzJgBwM0338y4ceMYPXo0kydPpqbxrQQzRIk/TSkbfod1/JHPIvkoqluo3n///UycOJHddtuNfv368dJLLwHw8ssvM2vWLF5//XXeeuutbXPpn3nmmbz44ossXryY2tpa/vrXvzY43le+8hWWLl3KihUrALjtttuYOnUq11xzDT169GDRokXEG/21evDBB7n//vt5/vnneeWVV7ZND33sscfy4osv8sorr7Dnnns2mEcok5T405Ty5i6nzIUtG3MTlEgBi+oWqvPmzWPKlCkATJkyhXnz5gEwfvx4hg4dSklJCWPGjGHZsmUAPP744xxwwAGMHDmSxx57jCVLljQ4nplx4okn8vvf/55PPvmEZ599lq9//evNxvDII48wdepUeoZJpV+/fgAsXryYL33pS4wcOZJ4PN7ktTIl0pG7ZnYuMA1w4DVgqrtvCNf9BPg5MNDdc3MPtFZocnOXYU7l6X8jNnw6PHob8Xf/RsXlO+nmLiIZUloalHeSLW+rVatW8dhjj7F48WLMjC1btmBmHHnkkQ2mXe7UqRObN29mw4YN/OAHP2DBggUMGzaMyy67LOn0yFOnTuVb3/oW3bt35/jjj6dz5+ZTa7KpmCFoML7//vsZPXo0t99+e4MbtmRSZFf8ZjYEOAsY6+77AJ2AKeG6YcBXgbyqkDdo+K2yYFK3g+8i/pfdKD+jq3r8iGRQ0m/ZPYPlbXX33Xdz0kknUVVVxbJly3jnnXcYMWIETz31VNLt65P8gAEDWLduXcpePIMHD2bw4MFcddVVnHzyyduWd+nSJenUzhMmTODWW2/dVsP/+OOPAVi7di2DBg2irq6uSXkok6Iu9XQGephZZ6An8F64/DrgAoJvAvmt9Hgq7v8NNRsbnqHq8SPSPkm7V89p3zfpefPmNWiIBZg8eTJ/+MMfkm6/0047MX36dEaOHMnRRx/NuHHjmok3xrBhw9hrr722LSsvL2fUqFHbGnfrTZw4kaOOOoqxY8cyZsyYbV01r7zySg444AC++tWvsscee7T1bbYo0mmZzexsoBKoBR5y95iZHQUc7u5nm9kygm8ETUo9ZlYOlAOUlpbuX5XsO18HUVISXOk3ZhZ8OxCRQCFPy3zmmWey7777cuqpp+bk9TvEtMxm1heYBIwABgO9zOwkoAK4pKX93X2Ou49197EDBw6MKsyMSH1zly3ZDUREcmL//ffn1Vdf5Xvf+16uQ0lLlKWeI4C33X2Fu9cB9wJTCf4QvBJe7Q8FXjKzXSKMI3Ipe/xMPhtWLYhsEIqIdAwLFy7kySefbNBA3JFF2aunGviCmfUkKPUcDtzr7ofVb9BcqSefNOnxUwqVP1tB7DN/JX7pWspvnUtNbfBRa45/KXaperRI27W2ZB91jf9y4NvAZuBlYJq7b0xYv4w0En8mb72YVRtWUjZ8M1UfNf1CM3x40ENIpJi8/fbb9O7dm/79+yv5Z4i7s2rVKtauXcuIESMarEtV4y+6e+5mW0mJ4970BFfDrxSjuro6li9fnrQvvLRd9+7dGTp0KF26dGmwPFXiL8pbL2ZTaallfBCKSL7q0qVLk6tSyT5N2RCxlA2/P3o4eR9QEZGIKfFHrMkglNItzDn/V8QGTYDnTyH+u03q8SMiWaVSTxbEYok9eDqB/wReqyH+m39TPncLNWFzt3r8iEg2qHE3h8qG1lD1bs8my9XjR0QyIesjd6Vl1e81Tfqgm7uISLSU+HMo9c1dNNWDiERHiT+HUvb4Oe4nsPpVTfUgIpFQ424OJZ3qYcZyYp/5I/FLVlM+9xZN9SAiGafG3Y6o9gPKytBUDyLSLhq5m0967EL1iuR/kNXwKyLtpRp/B1VamnwCK031ICLtpcTfQaVs+D3tHtiqXj8i0nZK/B1U06ketjLnp3cQG34czD+S+O1r1ONHRNpENf4OrOFUDyXgp8ObnYjPeoryWzprqgcRaRP16slDZcM2UrW86S3e1ONHRBJpyoYCUv1u8vt6qsePiKRDiT8PpZzqYcim7AYiInlJiT8PJe3x062GyqOnE79uvhp9RaRZSvx5qEmPn+Ew56Y66L075TPGUVUV3NyrvtFXyV9EEqlxt4CUDXeqqpsO/FKjr0hxUuNuEah+J/loXzX6ikgiJf4CkrLR9zMrYXNNdoMRkQ5Lib+AJG307b6JymPPhn8eQPyWd9XwKyJK/IUkaaPvLV2JVZxE/KEDKf/hTmr4FRE17haLstLNVL3TdIYONfyKFC417ha56uXJp2VSw69I8VHiLxIpG34HrQtqPyJSNJT4i0Ty0b4bqDy6HJ49kfgdtWr4FSkSSvxFImXD7xl7EY875eWuhl+RIqHGXaFs2AaqlndvslwNvyL5TY27klL1u02TPqjhV6RQNXsHLjPrChwJfAkYDNQCi4G/u/sb0Ycn2VBaGpR3miwfsgFI/kdBRPJXyit+M7sIeB44DHgF+C3wZ4I/FteZ2T/MbJ+sRCmRSjnN81GnwmtXEv/9VjX8ihSQ5q74X3P3q1Ksu9bMBgHDIohJsqz+Pr0VFUF5p7QUKq+A2K4Qv2kp5XM3UbMxuPLX/X1F8l/Kxl0zM+8gLb9q3M0Rd8qGrqfqvR2arFLDr0jH15bG3YUJO8+KJCrp2Myofr9p0gc1/Irks+YSf+Lk7l+OOhDpmFKO+B2q+/uK5KvmEn+7yzxmdq6ZLTGzxWY2z8y6m9nPzewNM3vVzO4zs53a+zoSnaQNv11rqJxUDv/5DfG4q+FXJM80l/j3MLOXzOzlhMcvmdnLZvZSSwc2syHAWcBYd98H6ARMAR4G9nH3UcC/gZ+2/21IVJLf33cTsWPeIz7rCcpP3agRvyJ5prnG3V2b29Hd32z2wEHifw4YDawB7gdmu/tDCdscAxzn7s32D1HjbgfkWykbspaq9/s0WaWGX5GOIVXjbsrunC0l9pa4+7tmNhOoJhj49VBi0g+dAvwx2f5mVg6UA5SmKjRL7lgJ1R80Tfqghl+Rjq65AVyPm9kZZja40fLOZvZlM5trZlOb2b8vMAkYQTDqt5eZfS9hfQWwGUhaGHD3Oe4+1t3HDhw4sHXvSrIidcNvXXYDEZFWaa7G/w2gC3CfmS0PG2P/DbwFTAVudPfbmtn/COBtd1/h7nXAvcCBAGb2feCbQKyjjBWQ1mu24ffNW9XwK9JBNVfqqQFmA7PNrBuwM1Dr7ivTPHY18AUz60lQ6jkcWGBmE4ELgUPC15A8lXTE7yW1xEqXEf/FRsrnflcjfkU6oEinZTazy4FvE5R0XgamAUuAbsCqcLPn3P305o6jxt08s3ULZUPXqeFXJMda3bibCe5+KXBpo8Wfi/I1pQMo6aSGX5EOTPPxSyQ04lek40or8ZvZUDM7LHzczcx6RRuW5LvUDb/T4d+/UsOvSA61mPjN7BSCefhvCRcNBx6IMijJfylH/B67gvj1T2vEr0gOtdi4a2aLgPHA8+6+b7js1XDKhaxQ424Bcads6Fqq3tuxySo1/IpkVnvuubvB3bcVZs2sEw1n7hRJnxnV7zdN+qCGX5FsSSfxP21mFwDdwzr/H4G/RhuWFLKUDb9DNmQ3EJEilU7ivwBYC7wBnA08ClREGZQUtmbv8bvoZ8R/t1kNvyIRarYff1jWudXdvw/cmJ2QpNAlv8evEdu9B/Gbl1E+t46ajcGpqRG/IpmXTuPuQ8A3wvl2ckKNu8WjbGgNVe/2bLJcDb8irdeekbtvAf9nZg8A6+sXuvvsDMYnAkD1e02TPqjhVyST0qnxryC4a1ZPYGDCj0jGpWz4HbQO3InHUf1fpJ0inaQtU1TqKR7xeFDTr0mYt7Vntw3MOfUU6Lc/5b/4MTU123sT9+wZDBRT/V+kqTaXeszsYZLceN3dJ2QoNpFtkjb8VnYltt9Yyr48uUHSh+APREWFEr9Ia6TTuHtAwtPuwGRgo7ufH2VgiXTFLwAlJY5707GDZrB1aw4CEung2nzF7+7PN1r0hJk9kbHIRNJUWmpUVSVbnv1YRPJZOpO07Zjws5OZHQ4MykJsIg0kn/FzPZWnzIW6tWr4FUlTOt05lxDU+I3gTlpvA9OjDEokmSb1/2FbqSx/gNiI6cQvfI3ym/6XmtpOgAZ+iTQnnRp/l8aDt8yss7tvjjSyBKrxS7M+eoqyfUZQtWJIk1Ua+CXFrD2zczau8QO80P6QRDJk54OpXjk46SoN/BJpKmWpx8x2Jqjl9zCzkWyfinlHgsFcIh1GyobfYVvRHUZFGmruf8Q3gBuAocCvgV+FPz8DLo4+NJH0pWz4Pf4CWPm8Gn5FEqS84nf324DbzOwEd78rizGJtFrSgV/nv0VsyF3EL32f8ltvo2ZDV0ANvyJpTdlgZl8D9iYYwAWAu18dYVwNqHFX2qxuDWWlm6j6YECTVWr4lULXnikbfg3sBHwZuI1g5O5zGY9QJApddqT6w+Sr1PArxSqdVq+D3f27wCp3vxg4gKDuL5IXWrrVo+r/UmzSutl6/W8z2yV8XhZZRCIZ1tytHuOVccrLnaoqcN9e/1fyl0KWTuL/u5ntBMwEFgHLgLujDEokk2KxYOrm4cODCd2GD4c5c4zY93ekYtZBKWf8FClUzTbumlkJMK5+ojYz6wH0cPePsxQfoMZdiY5m/JRC1qaRu+6+Fbg+4XlttpO+SJRKS5sm/WB5lgMRyaJ0Sj0Pm9mkyCMRyYGUA7+m3kz8tzVq9JWClE7iPxO4z8xqzexjM1ttZrrql4LQpP5fupU5l/wZPppP+Wmo0VcKUjqzc3ZKttzdt0QSURKq8Uu2lQ3bSNXybk2Wa9CX5JM2z84ZJvjjgQvDx4OAMZkPUaTjqH63adIHqK5ueaS7SEeXzh24bgAOA04MF9UAN0UZlEiupRz0tfNHsP4dDfqSvJZOjf9Adz+NcCBX2Kuna6RRieRY0kbf7nVUHj+D+IxLKZ9Wp/q/5K10br1YF/bndwAz6w+oh7MUtKSzfVZ2ITbpEsp27UnNhi4Ntq8f9KXZPiUfpJP4fwXcAww0s8uBE4DLI41KpAOIxZIl8hFUr0he59ekb5Iv0mncvQO4iGDKho+B4939znQObmbnmtkSM1tsZvPMrLuZ9TOzh83sP+Hvvu17CyLZlXLQlyZ9kzyR7j3pOgF1wKZ09zGzIcBZwFh33yc8xhRgBvCou38eeDR8LpI3mp307crbNembdHjp9OqpAOYBgwmmY/6Dmf00zeN3Jrhnb2eC+/S+B0wCfhuu/y1wdGuDFsmllJO+TR1AxexDNembdHjpDOBaCuzv7jXh857AQnffs8WDm50NVAK1wEPuHjOzT9x9p4RtVrt7s+UeDeCSfKFJ36QjafMALqCKho3AnYG30njBvgRX9yMIvi30MrPvpRcumFm5mS0wswUrVqxIdzeRnEpZ/x+mgV/ScaST+GuAJWZ2i5ndDLwGfGJmvzCzXzSz3xHA2+6+wt3rgHuBA4EPzWwQQPj7o2Q7u/scdx/r7mMHDhzYmvckkjMpJ337zuXwyWI1/EqHkE53zr+FP/XSvd9uNfCFsDRUCxwOLADWA98Hrgl/P5B2tCIdXJP+/8OcynMWEiu9gXjFm5TfOpeaDcH4x/qG38T9RLKhxRp/uw4e9Pv/NrAZeBmYBuwA3AWUEvxxOL6lOf5V45e8t2ElZWVQ9eGAJqs08ZtEJVWNP53G3YnAlcBwgm8IBri794si0GSU+KUQlJQEXTwbU8OvRKU9jbs3AKcBQ4CBwIDwt4i0QsqJ3watBd+q+r9kTTqJfzmwyN3r3H1L/U/UgYkUmuQDvzZQefRpxH96BeXTt2rgl2RFOo27FwB/MbP5wMb6he4+O6qgRApR8onfuhH74gTKDjicmtqG12Ga+E2ikk6N/0GC6RpeI2FWTne/ONrQtlONXwqdBn5JFNpT49/Z3Y9y9wp3v7j+J4IYRYpWyoFfQzcDmvhNMiudxP+omX0l8khEiljKgV+TphO/9i+a+E0yKp3EPx14xMzWmdnHZrbazJrtdy8irZN04rcba4gdt4qKn4/UxG+SUenU+DslW57Nnj2q8UsxU/1f2qrNNf4wwR8PXBg+HgSMyXyIIpJMSxO/qf4vrZXOfPw3AIcBJ4aLaoCbogxKRLZLWf//9iXE57xFeTmq/0urpFPjP9DdTwM2AITz6nSNNCoR2aZJ/b/UmfM/LxP7ws1UVJRQU9Nwe9X/pSXpDOCqM7MSwAHMrD8J/flFJHoNb/xuwMGw6V9Ur9ox6fa68bs0J+UVf3i7RIBfAfcAA8PZNp8C/icLsYlIc7r2SV3/TzEvkAg0X+p5AcDd7wAuAmYCqwmmUb4zC7GJSAtS1v9PvR02rVbDryTVXKln26WEuy8BlkQfjoi0RtMbv2yhctqfiH3uVOLnP0f5nF9Ss6ELoBu/yHYp+/Gb2XIg5a0V3b252y5mlPrxi7TS6kWU7bELVa3m30YAAAvkSURBVB/t0mSVbvxSPFL142/uir8Twd2ykhcRRaTj6juG6hXJL+rU8CvNJf733f2KrEUiIhlVWmpUVSVZPmgNbO1FfF6nRlNEqwRULJpr3NWVvkgeS33jl9OJX3gR5dO3aOBXkWou8R+etShEJOOSTvw2txuxC4+l4rc/pKa24TRcGvhVPFqcpK0jUOOuSGZp4rfi0J4bsYhIgUk58GtwDbir/3+BU+IXKULJ6/+1VE6aRrzif3Tj9wKnxC9ShJLW/2/pQuxH46m4+bspb/wuhUE1fhFpQPX/wqEav4ikJWX9f8gGQDd+KQRK/CLSQPL6fw2VR51C/JLrVf8vAEr8ItJA0vr/zZ2Inb4HFTceo/p/AVCNX0TSpvp/flGNX0TaTfX/wqDELyJpU/2/MCjxi0jaVP8vDEr8ItIqsVhwI5etW4PfsRO7wchLqF41LOn21dUqAXU0SvwikhGp6v/9+myivByVgDoQJX4RyYhU9X82f0pNTcPlKgHllhK/iGRE8vp/Fz5ePyDp9roFZO4o8YtIxjSt/3dJ3QV0WNDxX/X/7FPiF5FIJS0BdV1P5bFnEb/uScrLXfX/LIss8ZvZ7ma2KOFnjZmdY2ZjzOy5cNkCMxsfVQwikntJS0Cz3iH2tReouLqUmpqG3whU/49eVqZsMLNOwLvAAcDNwHXu/qCZHQlc4O6HNre/pmwQKUC+lZJOpikgIpTrKRsOB9509yrAgR3D5X2A97IUg4h0JFaSuv4/dDOg+n9UspX4pwDzwsfnAD83s3eAmcBPk+1gZuVhKWjBihUrshSmiGRTyvr/pGnE/9+fVP+PSOSJ38y6AkcBfwoXnQGc6+7DgHOBucn2c/c57j7W3ccOHDgw6jBFJAeS1v9vXE9syiYqfjFO9f+IRF7jN7NJwA/dfUL4/FNgJ3d3MzPgU3ffsbljqMYvUnw0BXT75bLG/x22l3kgqOkfEj7+CvCfLMQgInlGU0BHJ9LEb2Y9ga8C9yYsng78r5m9AlwNlEcZg4jkp2angL7o55oCuh06R3lwd68B+jda9hSwf5SvKyL5LxYLfldUBNM7lJZC5ZWdie0/mrJDvp1yCuj6/SQ1jdwVkQ6r6RQQXWGvC5udAhpUBmqJEr+I5J2U9f9dPiX+uw2aBroFSvwikneS1/83UHnMGVT8eIWmgW6BEr+I5J2k/f/ndid26Q+oXjUk6T6aBno7JX4RyUtN6v8xYOeDKR2WahroYMyS6v9K/CJSYCqvtuTTQHznSuI3/Uv1f5T4RaTANCkDlTpzrn6a2Pgbqbi4m+r/ZGla5vbSlA0i0m6b11PStWdRTQOR62mZRURyq3MvTQMdUuIXkaKRchroo04hftUdRTMNtBK/iBSN5NNA1xI7sSsV1x9SNNNAK/GLSFFp0g30lAFwwC1UrypNun0hTgOhxC8iQkvTQGwsqG6gSvwiIrQ0DcRHBdUNVIlfRITimgZCiV9EJFQs00Ao8YuItCDlNBDH/Zj4dU/mXTdQJX4RkRYknQZi5hvEjphPxdWledcNVIlfRCQNDcpAVUbsR/vDxJeoXjU86fYduRuoEr+ISFuZpe4GOng98fjWDtkNVIlfRKQdkncDraVy0nQqzn2/Q3YDVeIXEWmH5N1AuxG7YBLVKwcl3SfX3UCV+EVE2qlpN9ASGP7tDtsNVIlfRCQiKbuBHn8+8VlP5awbqBK/iEhEknYD/fliYl95iIrKoTnrBqrELyISoSbdQM86AL6+KKfdQJX4RUSyzUpSdwMdtJb47+oi7QaqxC8ikgMpu4EefRoV530YaTdQJX4RkRxIORvoxadQvXJw0n0y1Q20c2YOIyIirRWLhTOAbmPAEZSWBuWdxkqT3ySs1XTFLyLSwSQtA/UMlmeCEr+ISAeTtAw0p/G3g7ZTqUdEpANqWgbKHF3xi4gUGSV+EZEio8QvIlJklPhFRIqMEr+ISJFR4hcRKTLm7rmOoUVmtgJIMo4tLQOAlRkMJ1MUV+sortZRXK3TUeOC9sU23N0HNl6YF4m/PcxsgbuPzXUcjSmu1lFcraO4WqejxgXRxKZSj4hIkVHiFxEpMsWQ+OfkOoAUFFfrKK7WUVyt01HjgghiK/gav4iINFQMV/wiIpJAiV9EpMjkdeI3s4lm9i8z+6+ZzUiy3sxsdrj+VTPbL919I44rFsbzqpk9Y2ajE9YtM7PXzGyRmS3IclyHmtmn4WsvMrNL0t034rjOT4hpsZltMbN+4bpIPi8zu9XMPjKzxSnW5+rcaimuXJ1bLcWVq3Orpbiyfm6Fxx5mZo+b2VIzW2JmZyfZJrpzzN3z8gfoBLwJfBboCrwC7NVomyOBBwnuZ/YF4Pl09404rgOBvuHjr9fHFT5fBgzI0ed1KPDXtuwbZVyNtv8W8FgWPq8vA/sBi1Osz/q5lWZcWT+30owr6+dWOnHl4twKjz0I2C983Bv4dzbzVz5f8Y8H/uvub7n7JuBOYFKjbSYBd3jgOWAnMxuU5r6RxeXuz7j76vDpc8DQDL12u+KKaN9MH/s7wLwMvXZK7v4k8HEzm+Ti3GoxrhydW+l8Xqnk9PNqJCvnFoC7v+/uL4WP1wJLgSGNNovsHMvnxD8EeCfh+XKafnCptkln3yjjSnQqwV/1eg48ZGYLzaw8QzG1Jq4vmtkrZvagme3dyn2jjAsz6wlMBO5JWBzV59WSXJxbrZWtcytd2T630pbLc8vMyoB9gecbrYrsHMvnWy9akmWN+6am2iadfdsq7WOb2WEE/zkPTlh8kLu/Z2Y7Aw+b2RvhVUs24nqJYG6PdWZ2JHA/8Pk0940yrnrfAp5298QruKg+r5bk4txKW5bPrXTk4txqjZycW2a2A8Efm3PcfU3j1Ul2ycg5ls9X/MuBYQnPhwLvpblNOvtGGRdmNgq4BZjk7qvql7v7e+Hvj4D7CL7WZSUud1/j7uvCx38HupjZgHT2jTKuBFNo9FU8ws+rJbk4t9KSg3OrRTk6t1oj6+eWmXUhSPpxd783ySbRnWNRNFxk44fg28pbwAi2N3Ds3Wibb9CwceSFdPeNOK5S4L/AgY2W9wJ6Jzx+BpiYxbh2YfugvvFAdfjZ5fTzCrfrQ1Cr7ZWNzys8ZhmpGyuzfm6lGVfWz60048r6uZVOXDk8twy4A5jVzDaRnWN5W+px981mdibwT4JW7lvdfYmZnR6uvwn4O0HL+H+BGmBqc/tmMa5LgP7Ar80MYLMHs+99BrgvXNYZ+IO7/yOLcR0HnGFmm4FaYIoHZ1quPy+AY4CH3H19wu6RfV5mNo+gJ8oAM1sOXAp0SYgp6+dWmnFl/dxKM66sn1tpxgVZPrdCBwEnAq+Z2aJw2c8I/nBHfo5pygYRkSKTzzV+ERFpAyV+EZEio8QvIlJklPhFRIqMEr+ISJFR4peCZ2b9E2Zg/MDM3k14/kxEr7mvmd0SPr7MzH6S5n6PmFnfKGISqZe3/fhF0uXB6NUxECRhYJ27z4z4ZX8GXJXuxhZ0GDfgd8APgMqI4hLRFb8UNzNbF/4+1MyeMLO7zOzfZnaNBXPbvxDOyb5ruN1AM7vHzF4Mfw5KcszewCh3fyVh8V5mNt/M3jKzs8LtysL52H9NMJfNMODPBLNEikRGV/wi240G9iQYvv8WcIu7jw9vkvEj4BzgeuA6d3/KzEoJRk/u2eg4Y4HGN/7YAziMYO71f5nZjeHy3YGp7v6D+g3NrJuZ9feEeXZEMkmJX2S7F939fQAzexN4KFz+GkHSBjiC4Oq9fp8dzay3B3Oq1xsErGh07L+5+0Zgo5l9RDAlAECVB3OtJ/oIGAwo8UsklPhFttuY8HhrwvOtbP+/UgJ80d1rmzlOLdC9mWNvSTjeeprqHh5DJBKq8Yu0zkPAmfVPzGxMkm2WAp9ry8HDRt5dCG77JxIJJX6R1jkLGBve/Pp14PTGG7j7G0CfsJG3tfYHnnP3ze2MUyQlzc4pEgEzOxdY6+63tHK/64E/u/uj0UQmoit+kajcSMO6froWK+lL1HTFLyJSZHTFLyJSZJT4RUSKjBK/iEiRUeIXESkySvwiIkXm/wO/O86N+Ec+2wAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "import numpy as np\n", | |
| "N = 45\n", | |
| "hours = 2\n", | |
| "t = np.linspace(0,hours,N)\n", | |
| "T_numerical = np.empty(N)\n", | |
| "T_numerical[0]=85\n", | |
| "T_analytical = np.empty(N)\n", | |
| "exponent = np.exp(-0.275*t)\n", | |
| "delta_t = np.diff(t)\n", | |
| "\n", | |
| "for i in range(1, len(t)):\n", | |
| " T_numerical[i] = -0.275*(T_numerical[i-1]-65)*delta_t[i-1] + T_numerical[i-1] # Calculated numerical solution\n", | |
| " T_analytical = (85-65)*exponent + 65 # Calculated analytical solution\n", | |
| "\n", | |
| "print('3a. Time step of', delta_t[0].round(3))\n", | |
| "\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "plt.plot(t, T_numerical, '-',color = 'orange', label = 'Numerical')\n", | |
| "plt.plot(t, T_analytical,'o', color = 'Blue', label ='Analytical')\n", | |
| "plt.xlabel('Time (hr)')\n", | |
| "plt.ylabel('Temperature (F)')\n", | |
| "plt.title('Corpse Temperature vs Time')\n", | |
| "plt.legend();\n", | |
| "\n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 588, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "3b. The final temperature as t approaches infinity is 65.0 degrees Fahrenheit\n", | |
| "3c. Time of death is 9:07 am\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "N = 45\n", | |
| "hours = 5000\n", | |
| "t = np.linspace(0,hours,N)\n", | |
| "T_analytical = np.empty(N)\n", | |
| "exponent = np.exp(-0.275*t)\n", | |
| "\n", | |
| "for i in range(0, len(t)-1):\n", | |
| " T_analytical = (85-65)*exponent + 65 # Analytical solution for corpse temperature\n", | |
| "\n", | |
| "print('3b. The final temperature as t approaches infinity is', T_analytical[len(t)-1], 'degrees Fahrenheit')\n", | |
| "t_0 = 85-65\n", | |
| "num = 33.6/t_0\n", | |
| "time_death = np.log(num)/-0.275 # Time of death calulation \n", | |
| "print('3c. Time of death is 9:07 am')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "4. Now that we have a working numerical model, we can look at the results if the\n", | |
| "ambient temperature is not constant i.e. T_a=f(t). We can use the weather to improve our estimate for time of death. Consider the following Temperature for the day in question. \n", | |
| "\n", | |
| "|time| Temp ($^o$F)|\n", | |
| "|---|---|\n", | |
| "|8am|55|\n", | |
| "|9am|58|\n", | |
| "|10am|60|\n", | |
| "|11am|65|\n", | |
| "|noon|66|\n", | |
| "|1pm|67|\n", | |
| "\n", | |
| " a. Create a function that returns the current temperature based upon the time (0 hours=11am, 65$^{o}$F) \n", | |
| " *Plot the function $T_a$ vs time. Does it look correct? Is there a better way to get $T_a(t)$?\n", | |
| "\n", | |
| " 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": 570, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd5wV5d3+8c/F0pGiNEVBRCl2RSQqasCe8otpGNsjMSYmMZroo7FEUzSamGiMJqZZosaSaCyPKTbsgtIUVBQwAiooZRFFkLq7398fMxuOCMth2dnZPed6v17ntWfmTLlmF77nPveZuUcRgZmZlY8WeQcwM7PG5cJvZlZmXPjNzMqMC7+ZWZlx4TczKzMu/GZmZcaF3zZK0huSDtvAawdJmtHYmazpkzRT0v5557CPc+EvYZKelPSepDZZ7SMinomIgQ2xrY28wZwgaVn6WCGppmB6WUPsv6mS9DdJF+Wdo5CkU+r4eywCiIgdI+K5vLPax7nwlyhJfYGDgAA+l2uYBhARt0fEFhGxBfAp4J3a6XRekySphaRM/59Japnl9tcnIm4s+N1/AZhV8Pfo1th5bNO48Jeuk4BxwM3AqMIXJN0s6feSHkxbaGMlbS3p6vQTwnRJe6+zvX0lvZq+fpOktum2hkuaW7DtXpLukVQpabak7xa89hNJd0n6i6Slkl6RNCR97VagD/DPNNO5m3rAknpLul/SIkmzJH2r4LXLJd0u6c50+1Mk7SDpx+nyb0gaUbD8OEk/lfS8pCXpMXUueP0gSeMlvS/pBUnD1ln3EknjgeVAL0nfTH+vSyW9LulrBct/S9KjBdNtJYWk7dLf35eAH6a5/54uM1/SOZJeAT7Y2PGv83saLuktSSqYd5ykCenzYZImS/og3c/PN/VvUZDxwHr+/rdK/53MlzQnXc71qqFEhB8l+ABeB04D9gHWAD0LXrsZWJS+1hZ4HJhN8mZRAVwKPFGw/BvAVKA3sBUwFrg0fW04MDd93gJ4HvgR0BroB8wCjkxf/wmwEvh0up+fA+PW2c9hRRzbf/dZMK8CeBk4L933AOAt4JPp65eTFOERQEvgzvSYz0mnzwCmFWxvHPAmMAjYAvgncEP6Wl/gXeCw9Jg/DVQCWxasOwsYCLRKt/85YAdA6XorgF3T5b8FPFqw77Ykn9S2S6f/Bly0zvHOByYCvYB2Gzv+ddYVMAc4qGDeP4Ez0+eTgZHp847AJzby9zgKeH098+cDB9bz9/8g8FugPbBNmmlU3v+vSuXhd9ASlLaytgfuiojngZnA8essdl9EPB8RK4H7gJUR8ZeIqCb5T7lui//aiJgTEYuBy4Dj1rPrfYHuEXFJRKyOiFnA9cCxBcuMiYgH0v3cCuy5mYdb60CgbUT8It33a8BN6+z7sYh4IiKqgLuBTsCv0um/AYMktStY/qaImB4Ry4AfFxzzKODeiHg0Imoi4gHgVeCIgnVviIgZEbEmIqoi4h8RMTsSjwJPpZk3x68j4p2IWFHk8QMQSWW9s/Z4JG1F8mZ0Z7rIGmCApK4RsTQixm9mzlpF/f4lbQ8cDPxvRCyPiHnAb9Z3LFY/jd43aI1iFPBIRCxKp+9I5/26YJkFBc9XrGd63X7zOQXP3yRpaa5re5JujfcL5lUAzxRMzy94vhxoK6ll+p9/c2wP9F3Pvh8tmF73GCvTIlg7DdCh4Pm6x9w+7e7ZHjhO0siC11vx0d9J4bpI+hxwEbATyaeE9nz091Ifhfso5vgL3QE8nHYljSR5Q56XvjaK5NPZa5JeB34UEQ9vZlYo/ve/PcmnnsqC3qgWJJ9irQG48JeYtMV6DFAhqbbItgG6SNozIl6s56Z7FzzvA7yznmXmALMjon8997E5Q8XOAaZHxO6bsY11rXvMyyNiiaQ5JC36M+pY97/HIqkD8Hfgy8CDEVEl6SGSLheAD0neCGptvaFt1TF/k44/Il6QVNtddTxJ91/ta9OAr0iqIGll3ytpy4hYXcy2G8AcYBlJ15mHD86Au3pKz+eBamAXYK/0sTNJ6/Kkzdjud9IvG7cCfsDaboFCE4APJJ2XfmSvkLSbpH2L3McCku8F6mMMgKQz0y9HW0raQ9Lgem4P4KuSBkjagqQFXHvMtwAjJR2aHmO79Pm6BbtWO5JPBAuBmrT1P7zg9SnA3pJ2ldSe5DuSQsX8Xupz/H8l6WPfF7i3dqakk9JunmpgCckbTM1G9t9gImI2yfckv5TUUcmZUf1rvyi2zefCX3pGkfRNvxUR82sfwLXACar/qX93AI+QfGk5i+QL4I9IC8X/I3mzmU3yBfINQOd1l92AnwMXpWfKnLMp4SJiDcmXrAeQdMtUAn/g411Wm+JWkuL4NknhOzvd1yySM20uJjnGN4HvsYH/T2mX2zkkX6C+S/Lm/EDB6y8DvyR5c54OPLnOJq4jOavqfUl/28A+6nP8dwCHknwKWVIw/7PADElLSf4mxzRAV9ymOg7oQvL7WEzyptuzkTOULPmTlNnHSRpH8oX2bXlnMWtobvGbmZUZF34zszLjrh4zszLjFr+ZWZlpFufxd+vWLfr27Zt3DDOzZuX5559fFBHd153fLAp/3759mTRpUt4xzMyaFUlvrm++u3rMzMqMC7+ZWZlx4TczKzMu/GZmZcaF38yszLjwm5mVGRd+M7My48JvZtbUVK+C+Y/DlAtg+fruebR5msUFXGZmJS0ClkyFeaNh/iOw8GmoXgFqCd2HQfv13em0/lz4zczysGIezH80LfajYWV6p9ROA2HHr8PWh0PP4dCqY4Pv2oXfzKwxVC1PWvLzR8O8R5IWPkCbrkmRr3106F33dhqAC7+ZWRaiBt6bvLZFXzkGalZDi9bQ/SDY4cSk0G+5F6hxv2514TczaygfvpW26EfDgkdh1bvJ/C67w4AzkkLf4yBo2T7XmC78Zmb1teYDWPDk2u6bpa8l89tuDb0+k3bfHAbtts415rpc+M3MilVTBe9OTAr9/NGwaBxEFVS0gx6fhP7fSop9511ByjvtBrnwm5ltSAQsm1nQffM4rFkCCLYaDDt/H7Y5HLodABVt8k5bNBd+M7NCqxYnBb622H84O5nfvg/0GZmeZnkItO2Wb87N4MJvZuWtejUsem5t983iSckZOS07Qs8RsPPZSbHv2L9Jd99sChd+MysvEfDBtLWnWS58Eqo+BFVA16Gw6w+T7puuQ6FFq7zTZsKF38xK38qFyVWytd03K95O5m+xE+wwKu2+GQGtO+ebs5G48JtZ6alakVwwVdt9896UZH7rLZPTK2uvkt2ib64x8+LCb2bNX9TA+y8VXCX7DFSvTLpqug2DPS9Lr5IdDC0q8k6bu0wLv6QuwA3AbkAAX4uI5ySdAZwOVAH/johzs8xhZiVo+dsfvUp25cJkfuddYKf0fPoeB0OrLfLN2QRl3eK/BngoIr4sqTXQXtII4Ghgj4hYJalHxhnMrBSsWQYLn1rbfbPk1WR+2x4Fg5wdBu23zTdnM5BZ4ZfUCTgY+CpARKwGVkv6NnB5RKxK5y/MKoOZNWM11bD4+YKrZJ+FmjVQ0TYZ5KzfyUmx77J7ow9y1txl2eLvB1QCN0naE3ge+B4wADhI0mXASuCciJi47sqSTgVOBejTp0+GMc2syVg2u6D75jFY/V4yf8u9YOCZsM0R0P3ApPhbvWVZ+FsCg4EzImK8pGuA89P5WwL7AfsCd0nqFxFRuHJEXAdcBzBkyJCPvGZmJWL1+7DgibXFftnryfx228J2n0+7bw5NunOswWRZ+OcCcyNifDp9N0nhnwvcmxb6CZJqgG4knw7MrJTVrIFF49d237w7AaIaWnaAHsNhYDp0cadBJXOVbFOUWeGPiPmS5kgaGBEzgEOBV4GZwCHAk5IGAK2BRVnlMLMcRcDS/yRDFs8fnbTuq5YmffJbDYFdLkivkt0PKlrnnbZsZH1WzxnA7ekZPbOAk4EPgT9LmgqsBkat281jZs3YykVJ/3xt983yt5L5HXaAvsen3TeHJBdTWS4yLfwRMQUYsp6XTsxyv2bWiKpXQeXYgkHOXgACWnVO+ud3vSAd5GzHvJNaylfumtmmiYAlr6ztvln4NFQvB7WEbvvD7hcn3TdbDYEWLjFNkf8qZrZxK+Z/9CrZFfOS+Z0Gwo6npIOcDYdWHXONacVx4Tezj6taDgufgflpq/79l5P5bboWXCV7OHTonW9OqxcXfjNLBjl7b8ra7pvKMVCzGlq0Tq6S3evydJCzvXyVbAlw4TcrVx++9dGrZFelZ1V32R0GpOfT9zgIWrbPN6c1OBd+s3Kx5gNY8OTas28+mJHMb7cN9Pr02kHO2m2da0zLngu/WamqqYJ3JxYMcjYOogoq2iVXye70zaTYd97VV8mWGRd+s1KydGbyhey80bDgcVizBBBstQ/s/P3kNMtuB0BFm7yTWo5c+M2as1WLkwJf21f/4exkfvs+0Gfk2kHO2nTNN6c1KS78Zs1NBLz2W3jjdlg8KTkjp2XHZBiEnc9Or5Lt7+4b2yAXfrPmZtbN8Pz3kitjd/1hOsjZ0OT+smZFcOE3a06WvAqTToeeh8CIR3zjcKsXX4lh1lxUrYAxX0nGrj/gNhd9qze3+M2aixfOhCVTYfhDybn3ZvXkFr9Zc/DmXfD6dbDLedDryLzTWDPnwm/W1C2bBRO+kQx5vMdP805jJcCF36wpq16d9OvTAob91WfuWINwH79ZUzbl/ORc/YPugQ7b553GSoRb/GZN1dv/ghm/hgGnQ+8v5p3GSogLv1lTtHwuPDcqGf9+7yvyTmMlxoXfrKmpqYKxx0PNKhh2J1S0zTuRlRj38Zs1NVMvgcpnYP9bodOAvNNYCXKL36wpmf8YTL0U+n0Vdjgx7zRWolz4zZqKFQvg2ROh00AYcm3eaayEZVr4JXWRdLek6ZKmSdq/4LVzJIWkbllmMGsWogaeOwnWvA8H3pWMx2OWkaz7+K8BHoqIL0tqDbQHkNQbOBx4K+P9mzUP065I7py17x+Tm52bZSizFr+kTsDBwI0AEbE6It5PX/41cC4QWe3frNmofBZevDC5Y9ZOp+adxspAll09/YBK4CZJkyXdIKmDpM8Bb0fEi3WtLOlUSZMkTaqsrMwwplmOVi2Gscclt0ocer3vmmWNIsvC3xIYDPwhIvYGPgR+AlwI/GhjK0fEdRExJCKGdO/ePcOYZjmJgPGnwMp5cOCd0Lpz3omsTGRZ+OcCcyNifDp9N8kbwQ7Ai5LeALYDXpC0dYY5zJqm134Hc/8P9rwcuu6bdxorI5kV/oiYD8yRNDCddSjwQkT0iIi+EdGX5M1hcLqsWflYPBkmnw29PgODzso7jZWZrM/qOQO4PT2jZxZwcsb7M2v61iyFsV+BNt1hv5vdr2+NLtPCHxFTgCF1vN43y/2bNTkRMPHbsGwmHPoEtPVlLNb4PFaPWWOadTO8cTvsfgn0ODjvNFamPGSDWWNZMg0mnQ49D4Fdf5B3GitjLvxmjaFqBYw5JhmK4YDboEVF3omsjLmrx6wxvHAWLJkKwx+EdtvkncbKnFv8Zll78y54/U+w87nQ66i805jV3eJPT8P8NHAQ0AtYAUwFHoiI6dnHM2vmls2CCd+ArvvBnpfmncYMqKPFL+kiYDwwAngRuAX4B8mbxa8lPSRpt0ZJadYcVa+GMV8BWsCBf4MWrfJOZAbU3eJ/OSI21ET5paRtgN4ZZDIrDS9eAIsnwUH3QIft805j9l91Ff5/1LViRMwD5jVsHLMS8fa/YPpV0P870PuLeacx+4i6vtx9vvaJpKsbIYtZaVg+F54bBVvuBYOvzDuN2cfUVfgLBxDxJYZmxaipgrHHQ80qGHYnVLTNO5HZx9TV1eO7Y5ltqqmXQOUzsP+t0GlA3mnM1quuwj9I0gskLf+B6XPS6YiIwZmnM2tO5j8OUy+Ffl+FHU7MO43ZBtVV+H3HZ7NirVwIz54AnQbCkGvzTmNWpw0W/oiY2ZhBzJqtqIFn/wdWvwcjHk7G4zFrwuq6gOsJSd+W1Gud+S0lHSzpRkm+sYrZtCtg/iOwzzWw5R55pzHbqLq6ej4DfB24T9K2wGKgbfp4DPhdREzKPqJZE1b5HLx4IfQZCTudmncas6LU1dWzHPgN8BtJbYAewIqIWNRY4cyatNXvwdhjoX0fGHq9b6FozUZRwzJHxCpgTsZZzJqPCBh3Cqx4Bw4fC607553IrGgej9+sPl77Hcy9D/a+EroNzTuN2SbxePxmm2rxZJh8NvT6DAw6K+80ZpusqMIvaTtJI9LnbST5fDUrT2uWwtivQJvusN/NILedrPnZ6L9aSV8jGanzhnTW9sD9WYYya5IiYOJpsGwmDLsD2nbLO5FZvRTTXPkusB/wAUBEvEZyhs9GSeoi6W5J0yVNk7S/pCvS6Zck3SepS/3jmzWi2bfAG7fBbj+GHh630JqvYgr/yohYXTshqYKPjtxZl2uAhyJiELAnMA0YDewWEXsArwEXbFpksxwsmQYTvwM9R8CuF+adxmyzFFP4x0o6F2ib9vPfCfxrYytJ6kQynPONABGxOiLej4hHIqIqXWwcsF39ops1kqoVMOaYZCiG/W+DFhV5JzLbLMUU/nOBpcB04HskV+0W0+TpB1QCN0maLOmG9Xwp/DXgwU3Ia9b4XjgLlkyF/f8C7XttfHmzJq7Owp926/w5Iv4QEV+IiM+nz2uK2HZLYDDwh4jYG/gQOL9g2xcCVcDtG9j3qZImSZpUWVlZ7PGYNay3/g6v/wl2Phd6HZV3GrMGUWfhj4hqYBtJreqx7bnA3IgYn07fTfJGgKRRwGeBEyJivTd8iYjrImJIRAzp3r17PXZvtpmWzYLxX4eu+8Gel+adxqzBFHPl7izgGUn3k7TaAYiI39S1UkTMlzRH0sCImAEcCrwq6SjgPOCT6XhAZk1P9WoY8xWgBQz7K7SoT9vHrGkqpvBXkpyJ0z59bIozgNsltSZ5AzkZmAi0AUYrGdRqXER8axO3a5atFy+AxZPgoHtgi755pzFrUBst/BHxw/puPCKmAEPWmb1Tfbdn1ije/jdMvwr6fwd6fzHvNGYNbqOFX9Jo1nPj9Yg4IpNEZnlaPhfGjYIt94LBV+adxiwTxXT1XFTwvC3wJWBVNnHMclRTBWOPh+qVMOxOqGibdyKzTBTT1TN+nVlPSXoqozxm+Zl6CVQ+k5yv32lA3mnMMlNMV0+ngskWwD7ANpklMsvD/Mdh6qXQ76uww//kncYsU8V09bxC0scvkguuZgPfyDKUWaNauRCePQE6DYQh1+adxixzxRT+fhGxpnCGJN+5y0pD1MBzJyX3zx3xcDIej1mJK2asnnX7+AEmNHQQs1xMuxLmPQz7XA1b7pF3GrNGscGWu6QeJH357STtztqhmDux6RdymTU9lc/Biz+A3l+Gnb6ZdxqzRlNXl81nSEbP3A74fcH8pUC9L+oyaxJWvwdjj4X2feAT14OKvcWEWfO3wcIfETeRDKl8TETc1YiZzLIVAeNOgRXvwOFjobVvAmflpZjz+O+SdCSwK8kFXLXzf5ZlMLPM/Of3MPc+2PtK6DY07zRmja6Y8/h/D3QhuZvWTSRX7o7LOJdZNt6bAi/8L/T6NAw6K+80Zrko5qyeAyPieODddMC2T+DbJVpztGZpcgvFNt1gv1tAxfzzNys9xZyPv7L2p6StgXeBvpklMstCBEw8DZbNhEMeh7bd8k5klptiCv8DkroAVwJTgGrglkxTmTW02bfAG7fB7hdDz0/mncYsV3UWfkktgAcj4n3g75L+BbSLiMWNks6sISyZBhO/Az1HwK4X5p3GLHcbu+duDXBNwfQKF31rVqpWwNivJEMx7H8btKjIO5FZ7or5dmu0pKMzT2KWhRfOgvdfToZabt8r7zRmTUIxffynA50lrQJWkAzdEBGxVabJzDbXW3+H1/8EO58LvY7KO41Zk1FM4ffpD9b8LJsF478OXfeDPS/NO41Zk7LRrp6IqAZGAuelz7cB9so6mFm9Va+GMccCLWDYX6FFq7wTmTUpGy38kq4FRgC1tyVaDvwxy1Bmm+XFH8DiibDfjbBF37zTmDU5xXT1HBARgyVNBoiIxZJaZ5zLrH7e/jdM/xX0Pw16fzHvNGZNUjFn9axJz+cPAEldgZpMU5nVx/K5MG4UdNkTBv8q7zRmTVYxhf93wD1Ad0kXA2OAXxSzcUldJN0tabqkaZL2l7SVpNGS/pP+3HIz8pslaqqS++ZWr4QD74KKthtfx6xMFfPl7l+Ai0iGbFgMjIyIvxW5/WuAhyJiELAnMA04H3gsIvoDj6XTZptn6k9h4dOw7x+g04C805g1acUOT1gBrAFWF7uOpE4kQznfCBARq9OhH45m7Vg/twCf35TAZh+z4Imk8O8wCnb4n40vb1bmijmr50Lgr0AvkuGY75B0QRHb7gdUktzFa7KkGyR1AHpGxDyA9GePeqe38rb6PZh8LjzxqaSVP+TavBOZNQvFnNVzIrBPRCwHkHQZ8Dzw8yK2PRg4IyLGS7qGTejWkXQqcCpAnz59il3NykH1SpjxW3jlZ7BmCexwEux5GbTaIu9kZs1CMd02b/LRN4iWwKwi1psLzI2I8en03SRvBAskbQOQ/ly4vpUj4rqIGBIRQ7p3717E7qzk1VTDrFvgnwNgyrnQ/QD41BTY/2Zov23e6cyajWJa/MuBVyQ9THJK5xHAGElXAUTE/65vpYiYL2mOpIERMQM4FHg1fYwCLk9/3r/5h2ElLQLeeRBePD8ZcG2rfZNB13oOzzuZWbNUTOH/d/qotSn32z0DuD294GsWcDLJp4y7JJ0CvEUyHITZ+i2akLTuFz4FW+yUnKrZ+8sg5Z3MrNnaaOGPiBvru/GImAIMWc9Lh9Z3m1YmPvhPMvTCnLuhbQ8Y8jvY6Rsed8esAWy08Es6CvgpsH26vIdltuysWABTL4bXr4eKNrD7T2DQ/0KrjnknMysZxXT1XAscA7yMh2qwrKxZCtOuTMbZqV4FO50Ku/0I2vXMO5lZySmm8M8FpqS3YTRrWNWrYeb1MPUSWLkQ+oyEPS6DTv3zTmZWsoop/OcC/5T0JLCqdmZE/CarUFYGIpI7ZL34A1g2E3oMh4P/Cd2G5p3MrOQVU/gvJhmuoQvu6rGGsOCJ5IrbxZOgy+4w/AHY5iifqWPWSIop/D0iYp/Mk1jpe+8lmHIezHsI2veG/W6BvidAi4q8k5mVlWIK/2OSDomIxzNPY6XpwzfhxR/CG7dB6y6w9xUw4HQPnWyWk2IK/zeAcyQtJxmd06dzWnFWvQuv/Bxe+y0g2Pn7sOv50Nq3YDDLUzGFv1vmKay0VK2A136TFP2qpclwybtfDB16553MzCjuyt1qSccC/SLiZ5K2A3qSjNBptlZNNcy+BV76Eax4G3p9Fvb6OXTZLe9kZlagmCt3rwVakdxU5Wckg7b9Edg322jWbETA2/9KBlFb8ip0/QQMuwN6HJx3MjNbj2K6eg6IiMGSJgNExOJ00DUzqHwuOVOn8hnoOAAOvBt6f9GnZpo1YcUU/jWSWpAMyYykrvh8fvtgRjqI2r3Qtmdyr9sdT/EgambNwAYLv6SWEVEF/A64B+gu6WKScXsubqR81tSsmAcvXwwzb4CKdrD7JTDoLN/9yqwZqavFPwEYHBF/kfQ8cBjJqZwjI2Jqo6SzpmPNB/DqFTD9Kog10P802O2iZMhkM2tW6ir8/+2kjYhXgFeyj2NNTvVqeP2PMPWnsGoRbH8s7HEpdNwx72RmVk91Ff7uktZ7W0WAiLgqgzzWVEQNvHknvHghfDgbeh4Ce/0Cuq7vvjpm1pzUVfgrgC0oaPlbmZj/KEw+D957AbrsCcMfgm2O8Jk6ZiWirsI/LyIuabQklr/Fk2HK+TD/EeiwPex/K/Q9HtQi72Rm1oCK6uO3ErdsNrz0Q3jjdmi9FQy+KvnytqJN3snMLAN1FX7fEL3UrVwEr1wG//k9qAJ2uQB2OQ9ad847mZllaIOFPyIWN2YQa0RVy2HG1fDqL6BqGfT7WnJT8/bb5p3MzBpBMVfuWqmoqYJZN8HLP04uxNruaNjzZ9B5l7yTmVkjcuEvBxEw93548QL4YDp0OwCG3QU9Dsw7mZnlINPCL+kNYClQDVRFxBBJe5GM7tkWqAJOi4gJWeYoa5Vjk/vbLnoWOg2Cg+5LWvo+NdOsbDVGi39ERCwqmP4lcHFEPCjp0+n08EbIUV6WTEta+HPvh3bbwNDroN/J0MIf8szKXR5VIIBO6fPOwDs5ZChdK+YlN0KZ9Weo6JAMrzDoTGjZIe9kZtZEZF34A3hEUgB/iojrgDOBhyVdCbQADljfipJOBU4F6NOnT8YxS8SCp2DMl2HNEhhwBux6EbT1nTPN7KOyLvzDIuIdST2A0ZKmA18GzoqIeyQdA9xIMvLnR6RvEtcBDBkyJDLO2bxFJAOpTfoubNEPDnsaOu+cdyoza6IyvRY/It5Jfy4E7gOGAqOAe9NF/p7Os/qqXg0TvgkTT0vG0zlyvIu+mdUps8IvqYOkjrXPgSOAqSR9+p9MFzsE+E9WGUreigXw+CEw8/rkqtuD/wGtu+SdysyauCy7enoC9yk5bbAlcEdEPCRpGXCNpJbAStJ+fNtE706CZ74Aq96FA/4KfY/NO5GZNROZFf6ImAXsuZ75Y4B9stpvWZh9O0z4OrTpAYePha32zjuRmTUjHm+3Oamphsnfh+dOhK32haMmuuib2Sbz1TzNxer3YOxxMO9h6P9tGHw1VLTOO5WZNUMu/M3Bkmnw1Odg+Zsw9E+wk78WMbP6c+Fv6ub+E549AVq2g0Me98BqZrbZ3MffVEXA1Mvg6aOhY384cpKLvpk1CLf4m6KqD2HcyfDW32H74+ET10PL9nmnMrMS4cLf1CybDU9/Ht5/Gfb6Jex8jodQNrMG5cLflCx4AsaMTO6UNfwB6HVU3onMrAS5j78piIAZ18Ljh0Ob7nDkBBd9M8uMW/x5q14Fk74DM2+EXp+FYbdDq04bX8/MrJ5c+PO0Yh488yVY9BzseiHscQnIH8LMLFsu/Hl5dyI8/YXkitwD74I+I/NOZGZlws3LPMy+FUYflJ9Gtz4AAAi8SURBVNz/9ohnXfTNrFG58Demmip44Wx47iTotj8cORG2/NgApmZmmXJXT2NZtRjGHgvzR8OA02HwVdCiVd6pzKwMufA3hvdfSYZeWP4WDL0edvp63onMrIy58Gdt7v3w7InQsgMc+iR0PyDvRGZW5tzHn5WogZcvSYZf6DQIjprkom9mTYJb/FlYswzGjYI590LfE2HodcmwymZmTYALf0NbNitp5S95Bfb+FQw6y4OsmVmT4sLfkOY/BmOOAQKGPwTbHJ53IjOzj3Eff0OIgBm/gSeOhHZbJ4OsueibWRPlFv/mql4FE78Ns26C7Y6G/W+FVh3zTmVmtkEu/Jtj+TvwzBfh3fGw249g9x97kDUza/IyLfyS3gCWAtVAVUQMSeefAZwOVAH/johzs8yRiUXj4ZkvwJoP4MC7oc+X8k5kZlaUxmjxj4iIRbUTkkYARwN7RMQqST0aIUPDmnULTDgV2m0LRzwMXXbPO5GZWdHy6Or5NnB5RKwCiIiFOWSon5oqmPx9mHE19DwkGU65Tde8U5mZbZKsO6QDeETS85JOTecNAA6SNF7SU5L2Xd+Kkk6VNEnSpMrKyoxjFmHVu/DEUUnRH/g9GPGwi76ZNUtZt/iHRcQ7aXfOaEnT031uCewH7AvcJalfREThihFxHXAdwJAhQ4I8vT81HWRtLnziz7DjybnGMTPbHJm2+CPinfTnQuA+YCgwF7g3EhOAGqBbljk2y5x74ZH9oHoFHPaUi76ZNXuZFX5JHSR1rH0OHAFMBf4POCSdPwBoDSza0HZyEzXw0k+Se+J23hWOnATd9ss7lZnZZsuyq6cncJ+ScWpaAndExEOSWgN/ljQVWA2MWrebJ3drliZ3yZr7f7DDKBj6R6hom3cqM7MGkVnhj4hZwMfuKxgRq4ETs9rvZls6M+nP/2A6DL4aBn7Xg6yZWUnxlbuF5j+aDrIGjHgItj4s3zxmZhnw+AKQDLI2/dfpIGvbJjdBd9E3sxLlFn/1SpjwTZj9F9juC7D/X6DVFnmnMjPLTHkX/uVvp4OsTYDdL4bdLvIga2ZW8sq38Fc+lxT9qmVw0H3Q+/N5JzIzaxTl2bydeRM8NhxatocjnnPRN7OyUl4t/po18MLZ8Npvky9vh90JbbbKO5WZWaMqn8K/chGMPQYWPAEDz4K9fwktyufwzcxqlUfle++l5KKsFfNgv1ug30l5JzIzy03pF/637obnRkHrLnDY09BtaN6JzMxyVdpf7r7yMxgzErrsAUdNctE3M6PUC3/H/rDjKXDYk9Bum7zTmJk1CaXd1dNnZPIwM7P/Ku0Wv5mZfYwLv5lZmXHhNzMrMy78ZmZlxoXfzKzMuPCbmZUZF34zszLjwm9mVmYUEXln2ChJlcCb9Vy9G7CoAeM0Bz7m8uBjLg+bc8zbR0T3dWc2i8K/OSRNiogheedoTD7m8uBjLg9ZHLO7eszMyowLv5lZmSmHwn9d3gFy4GMuDz7m8tDgx1zyffxmZvZR5dDiNzOzAi78ZmZlpiwKv6SfSnpJ0hRJj0jqlXemrEm6QtL09Ljvk9Ql70xZkzRS0iuSaiSV7Cl/ko6SNEPS65LOzztPY5D0Z0kLJU3NO0tjkNRb0hOSpqX/pr/XkNsvi8IPXBERe0TEXsC/gB/lHagRjAZ2i4g9gNeAC3LO0ximAl8Ens47SFYkVQC/Az4F7AIcJ2mXfFM1ipuBo/IO0YiqgLMjYmdgP+A7Dfl3LovCHxEfFEx2AEr+G+2IeCQiqtLJccB2eeZpDBExLSJm5J0jY0OB1yNiVkSsBv4GHJ1zpsxFxNPA4rxzNJaImBcRL6TPlwLTgG0bavulfc/dApIuA04ClgAjco7T2L4G3Jl3CGsQ2wJzCqbnAp/IKYs1Akl9gb2B8Q21zZIp/JIeBbZez0sXRsT9EXEhcKGkC4DTgR83asAMbOyY02UuJPnYeHtjZstKMcdc4rSeeSX/CbZcSdoCuAc4c52ei81SMoU/Ig4rctE7gH9TAoV/Y8csaRTwWeDQKJELNjbh71yq5gK9C6a3A97JKYtlSFIrkqJ/e0Tc25DbLos+fkn9CyY/B0zPK0tjkXQUcB7wuYhYnnceazATgf6SdpDUGjgW+EfOmayBSRJwIzAtIq5q8O2XSEOwTpLuAQYCNSTDO38rIt7ON1W2JL0OtAHeTWeNi4hv5Rgpc5K+APwW6A68D0yJiCPzTdXwJH0auBqoAP4cEZflHClzkv4KDCcZongB8OOIuDHXUBmSdCDwDPAySd0C+EFEPNAg2y+Hwm9mZmuVRVePmZmt5cJvZlZmXPjNzMqMC7+ZWZlx4TczKzMu/FbyJHVNR2adImm+pLcLpp/NaJ97S7ohff4TSecUud6jkrbMIpNZrZK5ctdsQyLiXWAvSIowsCwirsx4tz8ALi124fSCHQG3AqcBJX9uvuXHLX4ra5KWpT+HS3pK0l2SXpN0uaQTJE2Q9LKkHdPluku6R9LE9DFsPdvsCOwRES8WzN5F0pOSZkn6brpc33S89d8DL5AMxfAP4Lisj9vKm1v8ZmvtCexMMvzvLOCGiBia3gTjDOBM4Brg1xExRlIf4OF0nUJDSO4NUGgQyaiwHYEZkv6Qzh8InBwRp9UuKKmNpK7pJxWzBufCb7bWxIiYByBpJvBIOv9l1g7lfRhJ6712nU6SOqZjptfaBqhcZ9v/johVwCpJC4Ge6fw3I2LcOssuBHqxdrgNswblwm+21qqC5zUF0zWs/b/SAtg/IlbUsZ0VQNs6tl1dsL0P17N+23QbZplwH7/ZpnmE5H4OAEjaaz3LTAN2qs/G0y95twbeqM/6ZsVw4TfbNN8FhqQ3sX8V+NiIpxExHeicfsm7qfYhGUm1aqNLmtWTR+c0y4Cks4ClEXHDJq53DfCPiHgsm2RmbvGbZeUPfLRfv1hTXfQta27xm5mVGbf4zczKjAu/mVmZceE3MyszLvxmZmXGhd/MrMz8f1mcmVqSTikfAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "def temp(time):\n", | |
| " ''' Returns the ambient temperature based on the time in hours between -3 and 2 where time = 0 is 11 am\n", | |
| " \n", | |
| " Argument\n", | |
| " ---------\n", | |
| " time: time in hours between -3 and 2 where time = 0 is 11 am\n", | |
| " Returns\n", | |
| " -------\n", | |
| " temperature: the ambient temperature at the specified time '''\n", | |
| " \n", | |
| " import numpy as np\n", | |
| " temperature = np.array([55, 58, 60, 65, 66, 67])\n", | |
| " if time == -3:\n", | |
| " return temperature[0]\n", | |
| " elif time == -2:\n", | |
| " return temperature[1]\n", | |
| " elif time == -1:\n", | |
| " return temperature[2] \n", | |
| " elif time == 0:\n", | |
| " return temperature[3]\n", | |
| " elif time == 1:\n", | |
| " return temperature[4]\n", | |
| " elif time == 2:\n", | |
| " return temperature[5]\n", | |
| " else:\n", | |
| " print('Enter hours between -3 and 2')\n", | |
| "\n", | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "temperature = np.array([55, 58, 60, 65, 66, 67])\n", | |
| "t = np.array([-3,-2,-1,0,1,2])\n", | |
| "plt.plot(t,temperature,color = 'Orange')\n", | |
| "plt.xlabel('Time (hr)')\n", | |
| "plt.ylabel('Temperature (F)')\n", | |
| "plt.title('Ambient Temperautre vs Time');" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "4a. The grpah does not look correct because it does not account for the time in between each temperaure change. A better way to get Ta(t) is by providing more points to potentially have a continuous function instead of discrete points." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 583, | |
| "metadata": { | |
| "scrolled": true | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "4b.\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUZfbA8e9JCB1RmoKUgKL0ZqSDICquujbE1R8KLmrsrMJaY0FXdl0XRdS1oCiWLFgQsdOkuqKAoFKFFQgBpDcJLeT8/rg3wyTMTCbJ3Jkkcz7Pkyczd+bee24COfOWe15RVYwxxhiAhFgHYIwxpuSwpGCMMcbHkoIxxhgfSwrGGGN8LCkYY4zxsaRgjDHGx5KCMcZTIjJORO6LdRwmPJYUDCLyfyKyUER+F5HNIvKliHSPdVyBiEhDN87cLxWR/X7Pe8Q6xqISkYru9dSPdSzhEpEzCvh9nK2qN6jq07GO1YSnXKwDMLElIkOBB4BbgSnAYeBC4DJgXiGPVU5VsyMepB9VzQCq+p1TgbaqusbL80ZCNH4+IpKoqke9PIc/Vf0F9/chIhWBA8CZqpoZrRhMZFlLIY6JSHXgCeAOVf1IVfer6hFV/VRV73XfU0FEnhORTe7XcyJSwX2tl4hkisj9IvIb8KbftodEZLuIrBORAX7nvEhElovIPhHZKCJ/9XvtEhFZIiK7ReS/ItKmiNdVyY1zg4j8JiIv+MV8oYisEZGH3fg2ujFdJiL/E5EdIjLM71hPich4EZnoxrxARFr6vd5ARCa7x/pVRG7Nt+9/ROQ9EdkHXCMi3UTkOxHZ4/48R4lI7oezOe73Ve6n7MtF5FYRme53zDytCRGZICLPi8hUEdkPdAl1/fl+TlXc85zut+1UETkgIieJyCki8pX7+9ghIl8X8fcxQUQeLuLPP1FEHnF/tttFJF1ETixKHCY8lhTiWxegIjApxHvSgM5AO6At0BF42O/1U4AaQCMg1W9bLeBUYBAwRkTOdF8bC9yiqtWAVsDXACLSAXgDuAWoCbwKfBLoj1kYRgH1gdbAmcAZOK2hXI2AI26cT7nnvQpoA5wHjBCRU/3e3w94y73OycBH7h+rROAL4L9APZwW1kMick6AfasDE93z3ukeqwfwR+Am97093e9nqmpVVf04zOu9DngEqAYsCOP6AVDV/cAnwLV+m68BpqjqLuB+YBXO77IuMDzMeApSmJ//vcAFQHf3mo6412e8oqr2FadfwADgtwLe8z/gIr/nfYF17uNeON1NFf1e7wVkA1X8tr0PPOI+zsD5w39CvvO8DPwt37ZVwDkFxKfA6X7Py7kxneq3rTewwn18IbAHSHCf13aP0dbv/cuAC93HTwGz8h1/B3A2cA6wOl88jwMv++07tYD4HwDGu48rurHU93v9VmC63/M87wEmAGPCvf4A578EWO73fBFwtfv4aeADoEmY/56Oi98vxoeL+PNfC3Tze60xkAVIrP//lNUvaynEtx1ALb/ui0DqAev9nq93t+XapqoH8+2zS51PoYH26QdcBKwXkdki0sXd3ggY5nZV7BaR3UCDfOcKRz0gCVjmd5yPgTr5Ys5xHx9wv2/xe/0AfuMWwIbcB+qMCWxyz9MISM4X81CcT8DH7QsgIi3EGcjfIiJ7gUdxPokXh/85wrl+f1OAk0WkrYicATQFPnVfG4FzrTPdLp+hxYwzV1g/fxERnH8DX/hdy2KcHo6aEYrF5GMDzfHtW+AgcDnwYZD3bML547fMfd7Q3ZYrUJndk0Skil9iaAgsBVDVBcBlIpKE043yPs5//A3ACFUdUfTLAWAzTkvlNFXdUcxj5WqQ+8DtMqqH8zPYC6xU1dYh9s3/83kNmAX0V9XfReQBnC6TQO8F2A9U9nt+SoD3+O9XqOtX1SMi8iFOF9IhYJKqHnBf2wP8BfiLiLTFSQ7fqeo3BR03ElRVRWQjcKWqLorGOY2NKcQ19z/9o8C/3UHNyiKSJCJ/EJHcKYTjgYdFpLaI1HLf/24Yh39cRMqLM0X0EuAD9/kAEamuqkdw/qjmzpR5DbhVRDqJo4qIXCwi1Qp5TUdw+qhHi0gt91gNROT8whwnn67iDIInAffhtLB+wJ2dJSJ3uwPA5USkjTs+Ekw1YI+bEFoCN/vFfgina6WJ3/uXAO1FpKWIVMb5+QdVxOv/D85YwrXuY9zrulREGruf2Pfg/K6iNrPJ9QrwlIg0cGOqIyJ/jHIMccWSQpxT1WdxujweBrbhfGK/E6fLAeBJYCHwE/Azzh/DJws47G/ALpxP0+nAraq60n3temCd23VyK84gKaq6EOcP5IvuvmuAG4p4WXe7516I88fsK+D0kHuENhEY7MbVD+inqkfdP8AXAV1xusi24YyNVA12IOAe4CYR+R34N/BevtcfxUmgu0XkUlX9Gadvfy6wEqeVUZDCXv8cIBFnMHy63/bm7vn2ue8Zqarzwzh/JD3txvS1ODO4/guESrqmmMQdvDEmIkSkF/CuqpaaG7BCEZGngFqqelOBbzamDLCWgjHGGB9LCsYYY3ys+8gYY4yPtRSMMcb4eHafgoi8gTMVcauqtnK3/Qvntv7DOHfK/llVd7uvPQjciDPlbYiqTinoHLVq1dLk5GRvLsAYY8qoRYsWbVfV2oFe86z7SER6Ar8Db/slhQuAr1U1W0T+CaCq94tIC5z58B1xbgyaDpyhBVR7TElJ0YULF3oSvzHGlFUiskhVUwK95ln3karOAXbm2zZVj5UOno9T4AqcMs0TVPWQqq7FmaPe0avYjDHGBBbLMYXBwJfu41PJW78l0912HBFJFWdBmIXbtm3zOERjjIkvMUkKIpKGU58lPXdTgLcF7NdS1TGqmqKqKbVrB+wSM8YYU0RRL4gnIoNwBqD76LEBjUz8io7hdCttyr+vMeaYI0eOkJmZycGD+YvUGuOoWLEi9evXJykpKex9opoURORCnIU7zlHVLL+XPgH+IyLP4gw0NwW+j2ZsxpQ2mZmZVKtWjeTkZJyadcYco6rs2LGDzMxMGjduHPZ+nnUfich4nNLMZ4qzPOONOMXOqgHTxFl28RUAVV2GU0J5OU7xrjsKmnlkTLw7ePAgNWvWtIRgAhIRatasWeiWpGctBVW9NsDmsSHePwJnUQ9jTJgsIZhQivLvIz7vaD68BxYOcb4bY4zxic+ksHclrH4ZvhsMVvvJmCIREYYNG+Z7PnLkSIYPH16kY61bt45WrVoBsHDhQoYMGRKJEMM2fPhwKleuzNatW33bqlYNtSyGN2666SaWL19e6P1mzZrFJZdcEpEY4jMp1OoE7Z6CDR/BqudjHY0xpVKFChX46KOP2L59e0SPm5KSwvPPe/v/8ujR44csa9WqxTPPPOPpeUM5evQor7/+Oi1atIhZDBCvSQGg2VA49VJYci9s/y7W0RhT6pQrV47U1FRGjRp13Gvr16+nT58+tGnThj59+pCRkQHADTfcwJAhQ+jatStNmjThww+PXxrc/1Pv8OHDGTx4ML169aJJkyZ5ksW7775Lx44dadeuHbfccovvD/1tt91GSkoKLVu25LHHHvO9Pzk5mSeeeILu3bvzwQcfHHfewYMH895777FzZ55CDHlaMZC3RdSrVy/uueceevbsSfPmzVmwYAFXXnklTZs25eGHHy4w1qpVq/Loo4/SqVMnvv32W3r16kVu6Z6vvvqKDh060LZtW/r06QPA999/T9euXWnfvj1du3Zl1apVwX49RRb1+xRKDBHoMg6+bA/zroY/LIYKNWIdlTFFs+hu2LUkssc8qR2c9VzIt9xxxx20adOG++67L8/2O++8k4EDBzJo0CDeeOMNhgwZwscfOyu8bt68mXnz5rFy5UouvfRSrrrqqpDnWLlyJTNnzmTfvn2ceeaZ3HbbbaxZs4b33nuPb775hqSkJG6//XbS09MZOHAgI0aMoEaNGhw9epQ+ffrw008/0aZNG8CZtz9v3ryA56latSqDBw9m9OjRPP744+H+lChfvjxz5sxh9OjRXHbZZSxatIgaNWpw2mmncc8997B169agse7fv59WrVrxxBNP5Dnmtm3buPnmm5kzZw6NGzf2JapmzZoxZ84cypUrx/Tp03nooYeYOHFi2LGGI36TAkD5k6Db+zC9O3w7CM6ZDBK/jSdjCuuEE05g4MCBPP/881SqVMm3/dtvv+Wjjz4C4Prrr8+TNC6//HISEhJo0aIFW7ZsKfAcF198MRUqVKBChQrUqVOHLVu2MGPGDBYtWsTZZ58NwIEDB6hTpw4A77//PmPGjCE7O5vNmzezfPlyX1L405/+FPJcQ4YMoV27dnnGSgpy6aWXAtC6dWtatmxJ3bp1AWjSpAkbNmxg3rx5QWNNTEykX79+xx1z/vz59OzZ03d/QY0azgfWPXv2MGjQIFavXo2IcOTIkbDjDFd8JwWAWh2h/UhY9BdY8Qy0uDfWERlTeAV8ovfS3XffTYcOHfjzn/8c9D3+UyMrVKjgexxOlWb/9ycmJpKdnY2qMmjQIP7xj3/kee/atWsZOXIkCxYs4KSTTuKGG27IM0+/SpUqIc914okn8n//93+89NJLvm3lypUjJyfH9zz/vP/c+BISEvLEmpCQEDJWcFouiYmJx21X1YDTSR955BF69+7NpEmTWLduHb169Qp5PUVhH4sBzrgLGvSDHx+ErXNjHY0xpUqNGjW4+uqrGTv22G1IXbt2ZcKECQCkp6fTvXv3iJ6zT58+fPjhh77ZQjt37mT9+vXs3buXKlWqUL16dbZs2cKXX35ZwJGON3ToUF599VWys52CzieffDJbt25lx44dHDp0iM8++ywisYbSpUsXZs+ezdq1a337gNNSOPVUp1bouHHjChVHuCwpgDO+0GksVG3ijC8c+C3WERlTqgwbNizPLKTnn3+eN998kzZt2vDOO+8wevToiJ6vRYsWPPnkk1xwwQW0adOG888/n82bN9O2bVvat29Py5YtGTx4MN26dSv0sWvVqsUVV1zBoUOHAEhKSvINBl9yySU0a9YsIrGGUrt2bcaMGcOVV15J27Ztfd1e9913Hw8++CDdunULOIMqEkr1Gs0RX2Rn988wpRPUPBvOnQEJ1rtmSq4VK1bQvHnzWIdhSrhA/05isshOqXRia+j4KmydAz+mxToaY4yJOksK+TW+Hk6/FVY8DRs+jnU0xhgTVZYUAjnrOaiRAvMHwb41sY7GGGOixpJCIIkVoMeHIOVg7pWQvT/WERljTFRYUgimSiPoNh52L4XvbrbCecaYuGBJIZS6F0DbJ2H9eCucZ4yJC5YUCtLiAah/OSweBltmxzoaY0qMQKWlX3nlFd5+++2ox7Jt2zaSkpJ49dVXi3WccePGceeddxb4nk2bji0hXxLKXUeSJYWCSAJ0eQuqngbfXA1ZG2MdkTFFkp4OycmQkOB8T0+P/DluvfVWBg4cGPkDu1Q1T8mJXB988AGdO3dm/Pjxnp07V/6kUBLKXUeSJYVwJJ0APSdBdhbMvQqOHop1RMYUSno6pKbC+vXO8Nj69c7zSCeG4cOHM3LkSMApK33//ffTsWNHzjjjDObOdUrIHD16lHvvvZezzz6bNm3a+D7d//777/Tp04cOHTrQunVrJk+eDDilq5s3b87tt99Ohw4d2LBhw3HnHT9+PM888wyZmZls3Hjsg1vVqlVJS0ujbdu2dO7c2VeA79NPP6VTp060b9+e884777jCfPv27aNx48a+gnN79+4lOTmZDz74gIULFzJgwADatWvHgQMHSkS560iypBCu6i2g8zjYMd8pnmdMKZKWBllZebdlZTnbvZSdnc3333/Pc8895ytHPXbsWKpXr86CBQtYsGABr732GmvXrqVixYpMmjSJH374gZkzZzJs2DBfwbxVq1YxcOBAFi9eTKNGjfKcY8OGDfz222907NiRq6++mvfee8/32v79++ncuTM//vgjPXv25LXXXgOge/fuzJ8/n8WLF3PNNdfw9NNP5zlmtWrV6NWrF59//jkAEyZMoF+/fvTv35+UlBTS09NZsmRJnsqwueWuJ06cyI8//uhbsyG33PXixYt54okneOihhyL8U44sq+NQGA37QYsHYfk/4KT20PSWWEdkTFjcNW7C3h4pV155JQBnnXUW69atA2Dq1Kn89NNPvgV29uzZw+rVq6lfvz4PPfQQc+bMISEhgY0bN/o+wTdq1IjOnTsHPMeECRO4+uqrAbjmmmu48cYbGTp0KOCsdZDbb3/WWWcxbdo0ADIzM/nTn/7E5s2bOXz4sK9Etb+bbrqJp59+mssvv5w333zTl1CCiWW560iypFBYbf7mLGay6C6o3hLqRLb6ozFeaNjQ6TIKtN1LuaWkc0tegzMu8MILL9C3b9887x03bhzbtm1j0aJFJCUlkZyc7CtTHark9fjx49myZQvpbl/Ypk2bWL16NU2bNiUpKclXgto/hrvuuouhQ4dy6aWXMmvWrIBrS3fr1o1169Yxe/Zsjh49mmf1tUBiWe46kqz7qLASEqHbf6BKMsy7CrIyYx2RMQUaMQIqV867rXJlZ3u09e3bl5dfftn3ifmXX35h//797Nmzhzp16pCUlMTMmTMLLC8NTrfS/v372bhxI+vWrWPdunU8+OCDvrLdwfiXoH7rrbeCvm/gwIFce+21edaKqFatGvv27TvuvbEsdx1JlhSKovyJ0PNj507nOVfC0YMF72NMDA0YAGPGQKNGTqX4Ro2c5wMGFP2YWVlZ1K9f3/f17LPPhrXfTTfdRIsWLejQoQOtWrXilltuITs7mwEDBrBw4UJfn304JarHjx/PFVdckWdbv379CpyFNHz4cPr370+PHj2oVatW0PcNGDCAXbt2ce211/q23XDDDdx6662+geZcsSx3HUlWOrs4NnwMc6+AxoOg85vO/zZjosRKZ3vvww8/ZPLkybzzzjuxDqXICls628YUiqPB5dDqMVj6uLPIebO7Yx2RMSZC7rrrLr788ku++OKLWIcSVZYUiqv1o7D7J+eO5xOaQ72+Be9jjCnxXnjhhViHEBM2plBckgBd3obqreCbP8HeX2IdkYkjpbn713ivKP8+LClEQlJV6DkZEpJg9h/h8O5YR2TiQMWKFdmxY4clBhOQqrJjxw4qVqxYqP2s+yhSqiZDj49gxrnwzTVwzufO9FVjPFK/fn0yMzPZtm1brEMxJVTFihWpX79+ofaxpBBJdXrA2S/B96mw+F44K7wpesYURVJSUsA7cY0pDksKkXb6zbD7Z1g1yqmXdPpNsY7IGGPCZmMKXujwLNTtCwtugy0zYx2NMcaEzZKCFxLKQbf3oFpTmNsP9q6OdUTGGBMWSwpeKV8dzvnUmbI6+xI4vCvWERljTIHiMilEYwUqAKqdBj0mwf61MLc/5JTskrnGGBN3SSFaK1D51OkBHcfAlhmw8E7npMYYU0LFXVKIyQpUTW5wFudZMwZWPuPhiYwxpng8Swoi8oaIbBWRpX7b+ovIMhHJEZGUfO9/UETWiMgqEfGsgFCsVqCi7ZPQsD8svg82TPL4ZMYYUzRethTGARfm27YUuBKY479RRFoA1wAt3X1eEhFPbgcOttKU1ytQIQnQ+S2o2Qn+OwB2LPD4hMYYU3ieJQVVnQPszLdthaquCvD2y4AJqnpIVdcCa4COXsQV0xWoylWCcyZDxVOcGkn7C15ZyhhjoqmkjCmcCmzwe57pbjuOiKSKyEIRWViUmi9erEBVKBXrQK/PndXaZl0Ch/dE6cTGGFOwkpIUAi1ZFnCajqqOUdUUVU2pXbt2kU42YACsWwc5Oc73qCWEXNWbO8Xz9q50bm47ejjKARhjTGAlJSlkAg38ntcHNsUikKjdw3DKudBprDNV9fubbaqqMaZEKClJ4RPgGhGpICKNgabA99EOIur3MDQZCK2fgLVvw8/DPTqJMcaEz8spqeOBb4EzRSRTRG4UkStEJBPoAnwuIlMAVHUZ8D6wHPgKuENVj3oVWzAxuYeh1cPQZDAsfQL+N9bDExljTMGkNK/alJKSogsXLozY8RISAvfiiDjjD57JOeIMOm+Z4SzOY+s8G2M8JCKLVDUl0GslpfuoRIjZPQwJSdDjA2ed53lXwc4fPD6hMcYEZknBT0zvYUg6AXp9ARVqwqw/wO+/RuGkxhiTlyUFPzG/h6FyPej1FeRkw8wL4aCtvWuMiS5LCvkUdA+D51NWqzdz1mHI2uCsw5C9P8InMMaY4CwpFELUpqzW7grdJsDOhTDvT07LwRhjosCSQiFEdcpq/csg5SXY9Dl8n2o3txljoqJcrAMoTaJedrvpLXBgMyx93KmZ1O4pj05kjDEOSwqF0LCh02UUaLtnWj8Gh7bC8n9ChTrQfKiHJzPGxDvrPiqEmExZFYGzXoAGV8HiYbD2HQ9PZoyJd5YUCiFmU1YTEqHru3DyuTD/z7DxC49PaIyJV5YUCilmU1YTK0DPj+Gkds5dz9u+idCBjTHmGEsKEeT5lNWkas5dz5UbwKyLYdeSCB3YGGMclhQiKCpTVivWgXOnOQliZl/Y+0sED26MiXeWFCIoalNWqzSE3tNAc+Dr82H/hoL3McaYMFhSiKCoVlmt3gx6T4Eju2HmBVYnyRgTESGTgoiUF5HLReQZERkvIm+IyFARaRatAEuTqE9ZrdHBqZO0f51TQO/wHo9OZIyJF0GTgog8DHwH9AZ+BN7CWTazHDBKRL4SkVZRibKUiMmU1To9oftE2P0TzL7YCugZY4ol6MprInKZqk4OuqNIXaCBqkZ9LeVckV55LRrS052B54wMp1tpxIgIJY2MD+Cba5x7Gc75FBIrRuCgxpiyKNTKa6HKXHwS6qCquhnYXJzA4k3ulNXcGUq5U1YhAomhYX/IzoL5N8C8q6HHRGdFN2OMKYRQYwqLch+IyHNRiKXM83zKapNBkPJv2Pgp/Pd6yDkaoQMbY+JFqJaC+D3u6XUg8SAqU1bPuB2yf4cl9zsthc5vQoLVPTTGhCfUXwsr4B9hUauy2uI+yDkMPz0COYega7p1JRljwhIqKTQTkR9wWgxnuo9xn6uqdvA8ujJmxIi8Ywrg4ZTVVg9DYiVY/FcnQXR7z6mfZIwxIYQaU2gN9Aeu8nuc+7y/96GVPeFMWY1oQb3mwyDlRcicDHMuh+wDxbwCY0xZF3RKamlQGqekhpJ/dhI4LYli3+vwv7Hw3c1wci/o+QkkVS1uqMaYUizUlNRQN6/NFJHbRKRevu3lRKSniIwVkT9HOth45tnspNNuhC5vw9Y5Tq2kw7uKeUBjTFkVqvvoYiAJmCQimSLyk4j8AvwK/Bl4WVXfjEaQ8cLT2UmNr4Pu78OuRTDjXDi4NQIHNcaUNUGTgqpmqerzqtoJOA0nSXRV1Yaq+mdVLTv9NiWE5wX1GlwJPT+Fvatg+jmQlRmhAxtjyoqwqqSq6iFV3aCq270OKJ5FpaBevb5OddUDm2BaD/j91wge3BhT2lnp7BKkoNlJEZuZVKcH9PkasvfBtO6w++cIXYExprSz2UelhCczk/Ysdwaes7OcZT5rd4lIrMaYkq1Is4/yHaC+iPR2H1cQkSqRDNAUzJOZSdVbwPnfQIVa8PV5sGlKsWI0xpR+BSYFERmMUzH1dXdTIyBoSW3jDc9mJlVNhvPnwQlnwJw/wvr3i3lAY0xpFk5LYQjQGdgLoKq/AHW8DMocz9OZSZVOhj4zoWYnZ02G1a9E4KDGmNIonKRwUFUP5z4RkUTyVlA1UeD5zKTyJzqzkupdBAtug58fh1I83mSMKZpwksI3InIfUNEdV3gP+MzbsEx+UambVK4y9JwEjQfBz8Nh4R22JoMxcabA2UduyyAVuACnhTAFeFVVc7wPL7R4mn1UkIjOTlKFJQ/AiqehQT/o+q4t72lMGRJq9lHIpOAmhDdUdZBXwRWHJYVjkpMDr9XQqBGsW1fEg654FhYPg5N7Q49JUL56MSI0xpQURZ6SqqpHgboiUugVWkTkDRHZKiJL/bbVEJFpIrLa/X6S32sPisgaEVklIn0Le75458nspOZDocs7sHUuTO8JWRuLcTBjTGkQzpjCr8Bc94/2kNyvMPYbB1yYb9sDwAxVbQrMcJ8jIi2Aa4CW7j4vua0UEybPZic1vg56fe6Uw5jaxbnhzRhTZoWTFLYB04DKQG2/r5BUdQ6wM9/my4C33MdvAZf7bZ/g1lhaC6wBOoYRm3F5Ojup7gVw3hzIOQJTuzktB2NMmVRgUlDVRwJ9FfF8J6vqZve4mzl2v8OpwAa/92W620yYPJ+dVKM9XPCtc0/D1+dDxocRvgJjTEkQao1mAERkGnDcaLSqXhDBOALd9xBwBFxEUnFmQ9Ew4ivel24DBgSfaZR/dtL69c7z3P3CUjXZKYsx+1KYdzW0HwnN7nGykDGmTAin++hh4BH3awSwEvixiOfbIiJ1AdzvuSu9ZAIN/N5XH9gU6ACqOkZVU1Q1pXbtAnuxjCtitZMq1IRzpztrMyweBouG2L0MxpQh4XQffef3NVtVh1D0/v5PgNzprYM4VkPpE+Aat9heY6Ap8H0Rz2ECiOjspHKVnFXcmg2DX16EuVdC9v5ixWeMKRnCKYh3gt/XiSLSB6gbxn7jgW+BM93lPG8EngLOF5HVwPnuc1R1GfA+sBz4CrjDnQ5rIiTis5MkATqMhJQXYdNnML0XHPitqOEZY0qIcO5o3oDTvy9ANrAWeFxVZ3sfXmh281r4PFmPIVfmp04hvYq14ZzP4MRWxTygMcZLxV1PoYm7LnMDVW2squcC30Q2ROM1T2cn1f8jnD8Hcg7D1K62LoMxpVg4LYUfVLVDQdtiwVoKkRORlkRWJsy6BPYshZQXoOltnsRqjCmeIrUURKSOiLQFKolIaxFp4351x7mRzZQhEZmdVLk+nD8X6l4IC26HRffYzCRjSplQ9ylcDAzGmR76kt/2fTjTU00ZErHZSUnVoOdk+GEorHoO9q2Gbv+BpBOKHaMxxntBk4Kqvgm8KSJXq6qt0VjGNWwYuMpqkWYnJSRCymg44UznPoapXeGcT6Fq42LHaYzxVjj3KbwvIn1FZKiIPJT7FY3gTPSEUzup0APRZ9wOveFSBjQAABdQSURBVL9yqqtO6Qhb50U4amNMpIVzn8JLODeaDQUqAdcBp3scl4mygmYn5Q5Er1/vrMGTWyajwMRwynnQdz6UPwm+Phd+Hef1pRhjiiGc2Uc/qWobEflRVduKSDVgYoRrHxWJzT6KnmIv4nNop1MvacsMOPMeaP80JBRYessY44Hi3qdwMPe7iJziPk+OUGymlCj2QHSFGtD7SzjjLlg1CmZdBId3RSw+Y0xkhJMUvhCRE4GRwBJgHWB1k+NMRMpkJCRByvPQ8TXYOgu+6gh7VkQiPGNMhIRMCiKSAHypqrtV9QOgMdBaVW2gOc5EdBGf02+Cc7+G7L0wpRNs/CwiMRpjiq+gNZpzgNF+zw+oav7V1EwcCGcgulAzk+p0h74LoFpTZ32GpU+C5nh8FcaYgoQz0Pw3YKGqTg75xhiwgeaSoVglMrIPwPepsO5dqH8FdHnLuQHOGOOZUAPN4SSFXUB14BBwAKdaqqpqjUgHWliWFEqGYs9MUnXufl58r3PDW4+P4YSmEY7SGJOruLOPagFJQFWgtvvcljwzPsWemSTiLOvZeyoc3AJTzoaNn0csPmNM+MK5o/ko0B+4331cF2jndWCm9IjYAj6nnAsXLoKqTWD2JfDTY1ZQz5goC+eO5heB3sD17qYs4BUvgzKlS0RLZFRpBOd/A01ugKVPOMnhkM1tMCZawuk+6qqqt+DexObOPirvaVSmVIl4iYxylaDTG3D2K84d0F+lwM7FUbseY+JZOAPN3wFdcGYgdRCRmsB0VW0fjQBDsYHm0qFYA9Hbv4O5/eDwDkj5N5w22IMIjYkvxR1o/jcwEagtIo8D84B/RjA+U8YVayC6Vif4ww9Qqxt8dyPMHwzZWQXvZ4wpknAGmt8GHsYpc7ET6K+qE7wOzJQd4QxEhxxzqFgHek+Blg/Dr2/C1C6wb42HERsTv8JpKQAkAkeAw4XYxxig4IHosMYcEhKh7d+g1xfOWtBfnQUbPoraNRgTL8KZfZQGjAfq4SzN+R8RedDrwEzZUdBAdKHWh673B6c7qdqZzljDwr/A0cOeX4Mx8SKcgeYVwFmqmuU+rwwsUtXmUYgvJBtoLhsSEpwWQn4ikBOsHNLRQ7Dkflg1GmqcDd3fh6rJXoZpTJlR3IHm9eRdy7kc8GskAjMGinjzW2IFOOs56DER9v0CX7aHDR97Ep8x8SScpJAFLBOR10XkNeBnYLeIPCsiz3obnokHxbr5rcGVbnfSaTD3Clhwh81OMqYYwlkP8XP3K9d8j2Ixccp/bCEjw2khjBhx/M1vueMOuQPRvn2rNnHugv7xIVj5LPw2Hbq8A7U6Rv1ajCntChxTKMlsTCE+FOrmt9++hvk3wIFN0OoRaPmQs+KbMcanWGMKInKhiCwQka0islNEdomIFaMxUVOom99OORcu+gkaXQs/D4ep3WDPSi/DM6ZMCWdM4UXgFuBUrHS2iYFCD0SXPxG6vuPMSNr/K3zVHlaOspXdjAlDOEkhE1iiqkdU9Wjul9eBGZOryAPRDfvDRUvh5PPgh6Ewozf8vjaKkRtT+oQz0Hwf8KmIzMJZfQ0AVX3eq6CM8Ve8gehT4JxPYO1bsOgv8EVraP8MnJ7q3AhhjMkjnJvXvsQpcfEz4Gt/q+oj3oZWMBtoNlCIgej9GU5BvS0z4JTzoNPrzvoNxsSZUAPN4bQU6qjqWRGOyZiICXsgukpDOHcarHnVWQ/681bQfqS1GozxE86YwgwROdfzSIwpooIGovOMNzQW0r+/FS76GWp2ggW3wswL4Pd10QrXmBItnKRwMzBdRH63KammJAo1EB20AuvkZKfV0PFV2D4fvmgFq160GUom7oWTFGoBSUB1bEqqKYFCVWENWYFVxOk6ungZ1O4Bi+6C6efA3lUxuQ5jSoKw7mgWkWuAJqr6dxGpD5ysqos8j64ANtBsChJ2BVZVWPsO/HC3UzupzePQbBgkhDPsZkzpUtw7ml8EegPXu5uygFeKGdBfRGSpiCwTkbvdbTVEZJqIrHa/n1SccxgDhVj1rbGQcPpAkoduJX35U7DkAZjSEXb+EJ1AjSkhwuk+6qqqtwAHAVR1J1C+qCcUkVY44xQdgbbAJSLSFHgAmKGqTYEZ7nNjiqXQq75tKEfqM3eTvv07OPgbTDnbmalklVdNnAgnKRwRkQRAAUSkJn73KxRBc2C+qmapajYwG7gCuAx4y33PW8DlxTiHMUAxVn17tiNcvBxOuwlWjHSmr26eFv0LMCbKgiYFEcntTP03MBGoLSKPA/OAfxbjnEuBniJS013F7SKgAc44xWYA93udYpzDGJ8BA5yb2HJynO+5CQEKuMeh/InO7KTzZjuVVmdeAP+9Dg5ujULUxsRGqJbC9wCq+jbwMDAS2AX0V9UJRT2hqq7ASSrTgK+AH4HscPcXkVQRWSgiC7dt21bUMIwBwhxzmNaT5DtXkHBdDsmXjSD9/jRY87pNXzVlUqik4LvFU1WXqepoVX1OVZcW96SqOlZVO6hqT2AnsBrYIiJ1AdzvAT+OqeoYVU1R1ZTatW1mrCmesMccMhJQFdZvb0Tqqy+QPuprmN4Tdi+LftDGeCjUfLvaIjI02IuqWuSlOEWkjqpuFZGGwJVAF6AxMAh4yv0+uajHNyZcBRXbCzjmcKgiaR+PYUDvRvBlO2g2FFo/CuWqRDd4YzwQ9D4FEdkMvIxfi8Gfqj5e5JOKzAVq4hTaG6qqM9wB7PeBhkAGTjdVyDun7T4F47WQ9zlkbYcl98Ovb0DlBpDyAtS/LPpBGlNIRS2It1lVn/AiIFXtEWDbDqCPF+czpqgaNgxcgbVhQ6BiLeg8lvTF95J29wlkbDuFhidvZcTfYcBgmydhSqewxhSMiVdhjTn8tRnrt9VDSWD9ljqk3laF9KcmwtGD0Q/YmGIKlRTsU7uJe0W6z+FwFdKeOQs+bwkbP49+0MYUQ1i1j0oqG1MwsRZ8zEHJ+aQl7F0Bp/4ROoyCaqdFP0BjAihW7SNjTHDB73MQ+MMSaPc06R/WIrlJORISlORGOc760caUUJYUjCmGkGMOieVJX3wvqWNfZ/32Rs59DhkJpN6cTXp66W2hm7LNkoIxxRDWmMOBvP/Nsg6UI23oFqvAakokG1MwxkPBxxxyyHm3HJw2GNqMgEonRz84E7dsTMGYGAk65tBAodk9pL99iOSGh9zxBrXxBhNzlhSM8VDQMYe/J5K+4hlS33yL9dsbuuMNYuMNJuYsKRjjoQLXjw423rD9+xhFbOKdJQVjPBZsPYegazlsqwNTO8F/ryP9jW0kJztjE8nJWPeS8ZwlBWNiJPh4A9DyIdLHJ5J6W+VjS4Wud8p4W2IwXrKkYEyMBB9vSIC2I0ibPJasw3nLcWdlOd1OxnjFkoIxMVLQPQ4ZmYGLGGdkqK36ZjxjScGYGAq1fnTQ7qWa62FKR/htBunp2JiDiShLCsaUUIG7l5QRaRvg4DbS/zaW1BsP2piDiSi7o9mYEiw9PchSoUcPktzgMOs3n3DcPo0aOa0OY4KxO5qNKaWCdi8lViTjt+MTAkBGRg7sWgpg3Uum0CwpGFNKBR9zyIAvW5Oe9k9Sb8627iVTKJYUjCmlgk5p/Wc1aD2ctDH/R9aBvDOYbEqrKYglBWNKqaBTWgfXhNaPkbGjfsD9MjKOjSNa95LJL/BEaGNMqTBgQN5prP4aNhTWrw+wveZ6+Ko/6Sv+RerD55CVJcCx7qXc45r4ZC0FY8qooFNa718F2ftJ+3uyLyHksu4lY0nBmDIqcPeSMOCvfeHiZWTsaBRwv9zuJetaik/WfWRMGRa0e0mEhg0J0r20gfQnZ5H69+t8pb2tayl+WEvBmDgVsHup0lFG3DaZtNHnHL/Wg3UtxQVLCsbEqYDdS68lMuCJu8jYEfgmCKcYn3UvlWXWfWRMHAvWvRRy5tIXF5O+/GlSH/uDdS+VQdZSMMYcJ+jMpQd/hXInkPZ0S+teKqMsKRhjjhN05tLQc6HvtwXOXALrXiqtrPvIGBNQkW+M++Iy0pf93bqXSilrKRhjCi1o99IDv0BiJeteKsUsKRhjCi1o99KwC6Dv/NDdSzlHAeteKqms+8gYUyRF7l6a3JX0ZU+T+vdryTqQCFj3UkliLQVjTMQF7V56ZAvU7ETa6B6+hJDLupdKBksKxpiIC9q9NKQT9JwU+ua4/RmAdS/FinUfGWM8UfTupcak//QYqaPTyDqYBFj3UjRZS8EYE3VBV417qhK0HUHaG6m+hJDLv3vJWhHeiUlSEJF7RGSZiCwVkfEiUlFEaojINBFZ7X4/KRaxGWO8F3TVuBtPhpYPkbG9XsD9MjJySH9mOqk359ja0x4RVS34XZE8ocipwDyghaoeEJH3gS+AFsBOVX1KRB4ATlLV+0MdKyUlRRcuXOh90MaYqEpODlzWu1HtjaBHWL89+fjXGsG6dV5HVjaIyCJVTQn0Wqy6j8oBlUSkHFAZ2ARcBrzlvv4WcHmMYjPGxFjQ7qVR9ULfA7Hvf4B1LxVH1JOCqm4ERgIZwGZgj6pOBU5W1c3uezYDdaIdmzGmZAjavTRAaNhQAu7TsOZ6+PR00v/6V1JvOmzdS0UU9aTgjhVcBjQG6gFVROS6QuyfKiILRWThtm3bvArTGBNjAwY43UE5Oc733FlHQVsR/6wEZ40m7Z17yTpYPs/rWVmQ9uBR33NrSQQXi+6j84C1qrpNVY8AHwFdgS0iUhfA/b410M6qOkZVU1Q1pXbt2lEL2hhTMgRtRQw+Gc4cQsa2kwPul7FBYEpn0p+aSOrNR60lEUQskkIG0FlEKouIAH2AFcAnwCD3PYOAyTGIzRhTCgRrRQA0DHxfHA3r7gWUtGfOKvBu6nhuScRiTOE74EPgB+BnN4YxwFPA+SKyGjjffW6MMYUStHvpXydC3+9CDFTnwI+PkP7ySlJTNW5bElGfkhpJNiXVGBNIerrzyT8jw2k5jBhxrDURdLprnc2se64+yUP+V+anvJbEKanGGOOZUN1LQVsSz9aFfttDT3nd8BEc3lWmu5csKRhj4krw6a5A+ZNCTHndAHP7kX7PnaTeeKDMdi9ZUjDGxJ0itSRG1YPz5pA26d9kHaqU5/WsLEi7dw/sWED6uzmluhVhScEYY/wEbUlcVw7q9CDjtxMD7pexuRrpj44i9caD+VoRWqoSgyUFY4zJp0hTXhsoaZNfI+tw3mZGVpaQNmwrrHoRdi8lPb1ktyQsKRhjTCEE7V76RyIZm6oE3Cdjay1YdBfpaX8ndXDJbklYUjDGmEIINVAdtBXRMAEu/ZW0j8cEb0mseMYdkzga05aE3adgjDERkp7uzETKyjq2rXLlY0kjIcFpIeQnkkPOu4mkf3Mtqa+/nidxVK6szlKmEVxxzu5TMMaYKAg53ZUCWhJXbAo+JjF0M/zwV9gwifRxezxtSVhLwRhjoqRYLYn/VCJ9bj9SX3+NrMNVAu4fLmspGGNMCVCslkT/vaR98kaehADHF/MrLksKxhgTRUW6cW4EkFiBjI0VAx4zIyNy8VlSMMaYEqLoLYnIxWBJwRhjSpAityQixJKCMcaUEgW1JCKhXOQOZYwxxmsDBkQ2CeRnLQVjjDE+lhSMMcb4WFIwxhjjY0nBGGOMjyUFY4wxPqW69pGIbAPWF+MQtYDtEQqnNIi36wW75nhh11w4jVS1dqAXSnVSKC4RWRisKFRZFG/XC3bN8cKuOXKs+8gYY4yPJQVjjDE+8Z4UxsQ6gCiLt+sFu+Z4YdccIXE9pmCMMSaveG8pGGOM8WNJwRhjjE9cJwUR+ZuI/CQiS0RkqojUi3VMXhORf4nISve6J4nIibGOyWsi0l9ElolIjoiU6WmLInKhiKwSkTUi8kCs4/GaiLwhIltFZGmsY4kGEWkgIjNFZIX7b/ovkT5HXCcF4F+q2kZV2wGfAY/GOqAomAa0UtU2wC/AgzGOJxqWAlcCc2IdiJdEJBH4N/AHoAVwrYi0iG1UnhsHXBjrIKIoGximqs2BzsAdkf4dx3VSUNW9fk+rAGV+1F1Vp6pqtvt0PlA/lvFEg6quUNVVsY4jCjoCa1T1V1U9DEwALotxTJ5S1TnAzljHES2qullVf3Af7wNWAKdG8hxxv8iOiIwABgJ7gN4xDifaBgPvxToIEzGnAhv8nmcCnWIUi/GYiCQD7YHvInncMp8URGQ6cEqAl9JUdbKqpgFpIvIgcCfwWFQD9EBB1+y+Jw2nKZoezdi8Es41xwEJsK3Mt37jkYhUBSYCd+fr8Si2Mp8UVPW8MN/6H+BzykBSKOiaRWQQcAnQR8vIjSqF+D2XZZlAA7/n9YFNMYrFeEREknASQrqqfhTp48f1mIKINPV7eimwMlaxRIuIXAjcD1yqqlmxjsdE1AKgqYg0FpHywDXAJzGOyUSQiAgwFlihqs96co4y8kGxSERkInAmkINTgvtWVd0Y26i8JSJrgArADnfTfFW9NYYheU5ErgBeAGoDu4Elqto3tlF5Q0QuAp4DEoE3VHVEjEPylIiMB3rhlJHeAjymqmNjGpSHRKQ7MBf4GefvFsBDqvpFxM4Rz0nBGGNMXnHdfWSMMSYvSwrGGGN8LCkYY4zxsaRgjDHGx5KCMcYYH0sKxnMiUtOtRLtERH4TkY1+z//r0Tnbi8jr7uPhIvLXMPebLiInhfG+u0WkcnHjdI/1oFvVdJWIBJwqKyKvh1v4TER6isgPIpItIlfle+0rEdktIp9FInZT9lhSMJ5T1R2q2s6tRvsKMCr3uap29ei0D+HcmxAWcSQA7wC3h7HL3UCxk4L7h/4aoCVOtc+X3GqneajqTaq6PMzDZgA34Nyln9+/gOuLFq2JB5YUTEyJyO/u914iMltE3heRX0TkKREZICLfi8jPInKa+77aIjJRRBa4X90CHLMa0EZVf/Tb3EJEZonIryIyxH1fsluX/iXgB5wSEZ8A1xYQ8xCgHjBTRGYW80dwGTBBVQ+p6lpgDU610/znnCUiKSKSKCLjRGSp+3O5J/97VXWdqv7EsZub/F+bAewLFZCI3Oz+bH90f9aV3e3jRORlt57/ryJyjruewQoRGVe0yzcljSUFU5K0Bf4CtMb5NHuGqnYEXgfuct8zGqelcTbQz30tvxScNRT8NQP64vzBfcytHwPOHe1vq2p7VV2vqruACiJSM1iQqvo8Tk2h3qp6XGVdERnl1z3m/xVo0ZtAlU1DlUJuB5yqqq1UtTXwZoj3FtVHqnq2qrbFKc18o99rJwHnAvcAnwKjcFo5rUWknQexmCgr8wXxTKmyQFU3A4jI/4Cp7vafOVbW/DycT/25+5wgItXc2vK56gLb8h37c1U9BBwSka3Aye729ao6P997t+K0BHZQBKp63Kf3EApb2fRXoImIvIBTwHFqiPcWVSsReRI4EagKTPF77VNVVRH5Gdiiqj8DiMgyIBlY4kE8JoosKZiS5JDf4xy/5zkc+7eaAHRR1QMhjnMAqBji2Ef9jrc/wP4V3WMUiYiMIvDaHBNU9al82wpV2VRVd4lIW5xWzx3A1TjrYkTSOOByVf1RRG7AqS2Uy/93kv/3ZX9PygDrPjKlzVScdS8ACNJlsQI4vSgHd6tQngKsc5/PEJFA3Tn7gGqBjqGq9/gNpPt/5U8I4IxhXCMiFUSkMdAU+D5EfLWABFWdCDwCdCjM9YWpGrDZ7WIb4MHxTQlmScGUNkOAFBH5SUSWA8dVeFXVlUB1d8C5sM7CqRyb7c5GOp3Ayz2OAb4s7kCzqi4D3geWA18Bd6jq0RC7nArMEpElOJ/oj1tjW0TOFpFMoD/wqtu1k/vaXOADoI+IZAaZAvsIzmpe04iDcvImL6uSasokd1bOPlUNNBAdar/RwCeqOkNEWgGDVXWoJ0EaUwJZS8GUVS+Tt887XEvdaZuo6lJLCCbeWEvBGGOMj7UUjDHG+FhSMMYY42NJwRhjjI8lBWOMMT6WFIwxxvj8P4sGRQ63vnfjAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "import numpy as np\n", | |
| "N = 40\n", | |
| "hours = 2\n", | |
| "t = np.linspace(0,hours,N)\n", | |
| "t_BD = np.linspace(0,-3,N)\n", | |
| "t_analytical = np.linspace(-3,hours, N)\n", | |
| "T_numerical = np.empty(N)\n", | |
| "T_numerical[0]=85\n", | |
| "T_numerical2 = np.empty(N)\n", | |
| "T_numerical2[0] = 85\n", | |
| "exponent = np.exp(-0.275*t_analytical)\n", | |
| "delta_t = np.diff(t)\n", | |
| "delta_t1 = np.diff(t_BD)\n", | |
| "\n", | |
| "# Caluating a piecewise function for the numerical solution\n", | |
| " \n", | |
| "for i in range(0, len(t)-1):\n", | |
| " time1 = int(t[i])\n", | |
| " temp_a = temp(time1)\n", | |
| " T_numerical[i+1] = -0.275*(T_numerical[i]-temp_a)*delta_t[i] + T_numerical[i] # Numerical solution from t = 0 to 2s\n", | |
| "\n", | |
| "\n", | |
| "for n in range(0, len(t_BD)-1):\n", | |
| " if t_BD[n] == 0:\n", | |
| " time2 = int(t_BD[n])\n", | |
| " temp_a = temp(time2)\n", | |
| " T_numerical2[n+1] = -0.275*(T_numerical2[n]-temp_a)*delta_t1[n] + T_numerical2[n] #Numerical solution from t = -3 to 0s \n", | |
| " \n", | |
| " elif t_BD[n] == -1:\n", | |
| " time2 = int(t_BD[n])\n", | |
| " temp_a = temp(time2)\n", | |
| " T_numerical2[n+1] = -0.275*(T_numerical2[n]-temp_a)*delta_t1[n] + T_numerical2[n] \n", | |
| " \n", | |
| " elif t_BD[n] == -2:\n", | |
| " time2 = int(t_BD[n])\n", | |
| " temp_a = temp(time2)\n", | |
| " T_numerical2[n+1] = -0.275*(T_numerical2[n]-temp_a)*delta_t1[n] + T_numerical2[n] \n", | |
| " \n", | |
| " elif t_BD[n] == -3:\n", | |
| " time2 = int(t_BD[n])\n", | |
| " temp_a = temp(time2)\n", | |
| " T_numerical2[n+1] = -0.275*(T_numerical2[n]-temp_a)*delta_t1[n] + T_numerical2[n] \n", | |
| " \n", | |
| " else:\n", | |
| " time2 = int(t_BD[n])\n", | |
| " temp_a = temp(time2-1)\n", | |
| " T_numerical2[n+1] = -0.275*(T_numerical2[n]-temp_a)*delta_t1[n] + T_numerical2[n] \n", | |
| " \n", | |
| "\n", | |
| " \n", | |
| "time1 = int(t_analytical[i])\n", | |
| "temp_a = temp(time1)\n", | |
| "T_analytical = (85-65)*exponent + 65 #Linear analytical solution \n", | |
| "\n", | |
| "print('4b.')\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "plt.plot(t, T_numerical, '-',color = 'orange', label = 'Nonlinear Numerical')\n", | |
| "plt.plot(t_BD, T_numerical2, '-',color = 'orange')\n", | |
| "plt.plot(t_analytical, T_analytical,'o', color = 'Blue',label = 'Linear Analytical')\n", | |
| "plt.xlabel('Time (hr), t = 0 is 11 am')\n", | |
| "plt.ylabel('Temperature (F)')\n", | |
| "plt.title('Corpse Temperature vs Time')\n", | |
| "plt.legend();\n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The numerical solution is now nonlinear or follows more of an exponential shape since it accounts for the change in temperature. This difference is seen in the graph above when graphing the analytical and numerical solution. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 573, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "4b. The time of death is -1.5490549054905491 hours before 11 am which is approximately 9:27 am for a body temperature of 98.602 F\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "for j in range(len(t_BD)):\n", | |
| " body_temp = T_numerical2[j].round(3)\n", | |
| " #print(body_temp)\n", | |
| " if body_temp >= 98.6 and body_temp <= 98.603:\n", | |
| " print('4b. The time of death is', t_BD[j], 'hours before 11 am which is approximately 9:27 am for a body temperature of',body_temp,'F')\n", | |
| " \n", | |
| " " | |
| ] | |
| }, | |
| { | |
| "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 | |
| } |