01_Getting-started-project.ipynb

{
 "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": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K: 0.275 1/hr\n"
     ]
    }
   ],
   "source": [
    "T_a = 65 #deg f\n",
    "T_o = 85 #deg f @ 11am\n",
    "T_f = 74 #deg f @ 1pm (2 hours later)\n",
    "dt = 2\n",
    "dT_dt = (T_f - T_o) / 2\n",
    "K = -dT_dt / (T_o - T_a)\n",
    "print('K:' , K , '1/hr')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Change your work from problem 1 to create a function that accepts the temperature at two times, ambient temperature, and the time elapsed to return $K$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.275"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def measure_K(Temp_t1,Temp_t2,Temp_ambient,delta_t):\n",
    "    ''' Determine the value of K based upon temperature of corpse \n",
    "    when discovered, Temp_t1\n",
    "    after time, delta_t, Temp_t2\n",
    "    with ambient temperature, Temp_ambient\n",
    "    Arguments\n",
    "    ---------\n",
    "    Temp_t1: initial temperature\n",
    "    Temp_t2: final tempertature\n",
    "    Temp_ambient: ambient temperature\n",
    "    delta_t: time elapsed\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    Calculates K, the empirical constant using the finite difference approach.\n",
    "    \n",
    "    '''    \n",
    "    return -((Temp_t2 - Temp_t1) / delta_t) / (Temp_t1 - Temp_ambient)\n",
    "\n",
    "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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytical_T(t):\n",
    "    T_0 = 85\n",
    "    T_a = 65\n",
    "    K = 0.275\n",
    "    return T_a + (T_0 - T_a)*np.exp(-K*t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Part A:\n",
    "T_0 = 85\n",
    "T_a = 65\n",
    "K = 0.275\n",
    "\n",
    "N1 = 5\n",
    "times1 = np.linspace(0,2,N1)\n",
    "dt = np.diff(times1)\n",
    "\n",
    "analytical = np.empty(len(times1))\n",
    "for i in range(len(times1)):\n",
    "    analytical[i] = T_a + (T_0 - T_a)*np.exp(-K*times1[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Numerical Solution for N = 5\n",
    "numerical_T1 = np.empty(len(times1))\n",
    "numerical_T1[0] = 85\n",
    "\n",
    "for i in range(1, len(times1)):\n",
    "    numerical_T1[i] = numerical_T1[i-1] - K*(numerical_T1[i-1] - T_a)*dt[i-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Numerical Solution for N = 50\n",
    "T_0 = 85\n",
    "T_a = 65\n",
    "K = 0.275\n",
    "\n",
    "N2 = 50\n",
    "times2 = np.linspace(0,2,N2)\n",
    "dt = np.diff(times2)\n",
    "numerical_T2 = np.empty(len(times2))\n",
    "numerical_T2[0] = 85\n",
    "\n",
    "for i in range(1, len(times2)):\n",
    "    numerical_T2[i] = numerical_T2[i-1] - K*(numerical_T2[i-1] - T_a)*dt[i-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3gUVffA8e8hCRAQCFWpEiwgSC8qJfoqzUKRoqDSi6KIWBB+8qrYFfBVsSvNgigdBJEiAlIEQhGkClID0qskEOD8/tgJbMIm2ZAtKefzPPuwM3tn5uw43tw9c+deUVWMMcZkHzmCHYAxxpjAsorfGGOyGav4jTEmm7GK3xhjshmr+I0xJpuxit8YY7IZq/iNyYDEZZSIHBWR5QE+9iAR+dZ5X0ZETolISCBjMP5lFb9JkYjsEJH9IpLXbV13EZnv4+PMF5HuPtpXWRHZcQXbjRaR130Rgw/UBxoBpVS1jqcCIlJcREaIyD4ROSkim0TkFff/VumlqrtU9SpVPe+rfZrgs4rfeCMUeCrYQWQz1wI7VPVfTx+KSCFgKRAO3Kaq+XD9oYgArgtYlCZTsorfeGMI8JyIRHj6UEQqiMgcETkiIptF5AFnfaSIHBORHM7ycBE54LbdtyLSV0TeABoAHzlphY+cz+uKyAoROe78W9dt2/ki8pqILHZau7NFpEgy8fUXkRin3GYRuSutJ0BEPhCR3SJyQkRWikgDZ31uEYlNOLaI/FdEzolIfmf5dRF5P5l9lhCRac552yoiPZz13YDhwG3O+XjFw+bPACeBR1R1B4Cq7lbVp1R1rRfnz+OxPcRYVkRUREKd5RTPu4h0FJGdInJYRF50fjE2TNvZNv5mFb/xRjQwH3gu6QdOWmEO8B1QDGgPfCIilVR1O3ACqO4UbwCcEpGbnOUoYIGqDgR+A3o7aYXeTot2BjAMKAz8D5ghIoXdDv8Q0MU5bs6E+FR1h6qWdeIrD/QGajut4ibAjis4ByuAakAh57uOF5HcqhrnfHa723faCdRz/47J7HMssAcoAbQB3hSRu1R1BPAYsNQ5Hy972LYhMElVL3jasRfnz+OxUz8NQDLnXUQqAp8ADwPFgQJASS/3aQLIKn7jrZeAJ0WkaJL19+FKSYxS1XOqugqYiKsyAVeld7uIXOMsT3CWI4H8wB/JHO9e4C9V/cbZ71hgE9DMrcwoVd2iqrHAOFwVc1LngVxARREJc/4obEvTNwdU9VtVPezE8q6zz/JJvmMoUAVXZXu7iOQGauP6o5aIiJTGlcfvr6pxqroGVyu/g5chFQb2pfB5sufPB8dO7ry3AX5U1UWqehbXNWODgWVAVvEbr6jqn8B0YECSj64FbnFSOsdE5BiuFl9CRb8AuANXy3chrl8Otzuv35JrseJqie5Msm4niVuQ/7i9Pw1c5SHurUBfYBBwQES+F5ESyX7RZIjIsyKy0UmbHMPVmk1IcSR8xxrAOly/gG4HbgW2quohD7ssARxR1ZMpfL+UHMbVqk5OSucvvcdO7ryXAHYnfKCqp504TQZjFb9Ji5eBHiSuIHbjStdEuL2uUtVezucLcKV47nDeL8KVBrmdxCmQpC3Dvbj+qLgrA8SkNWhV/U5V6zv7U+CdtGzv5PP7Aw8ABVU1AjgOiFNkCa7W//24zsUGJ9Z7ST7NsxcoJCL53Nal5fvNBe5PuH+SzP6TO3/pPXZy9gGlEhZEJBzXLxOTwVjFb7zmtJ5/APq4rZ4O3CgiHUQkzHnVTsjjq+pfQCzwCLBQVU8A+4HWJK4U9wPl3JZ/cvb7kIiEisiDQEXneF4TkfIicqeI5ALinFhS6poY4tywTXjlBPIB54CDQKiIvIQrTZVwXk4DK4En3L7TEuBRkqn4VXW3U+Yt5zhVgG7AGC+/2v+cGL4SkWud71pSRP7n7CvZ8+eDYydnAq5UUl3nvL3CpT+OJgOxit+k1avAxX7iTrqgMdAOV0vyH1wt6lxu2ywADqvqLrdlAVa7lfkAaCOuB5aGqephXPcPnsWVLngeuC+ZtElKcgFvA4ec2IoBL6RQfgCuPw4Jr3nALGAmsAVXSiQOt5SG23cKA5a7LefDld5KTnugLK7zNhl4WVXnePOlVPUIUBeIB5aJyEngF1y/RLZ6cf6u+NgpxLQeeBL4Hlfr/yRwADiTnv0a3xObiMUY4w8ichVwDLjB6eFlMghr8RtjfEZEmolIHqeb71BcN7t3BDcqk5RV/MYYX2qBK320F7gBaKeWVshwLNVjjDHZjLX4jTEmmwkNdgDeKFKkiJYtWzbYYRhjTKaycuXKQ6qa9Gn7zFHxly1blujo6GCHYYwxmYqIJH16G7BUjzHGZDtW8RtjTDZjFb8xxmQzmSLHb4xJXXx8PHv27CEuLi7YoZgAy507N6VKlSIsLMyr8lbxG5NF7Nmzh3z58lG2bFlEbGy07EJVOXz4MHv27CEyMtKrbbJmxb/ofRadLkP/VRHsPRZLiYhw3qlxjPp5dkH9vsGOzhi/iIuLs0o/GxIRChcuzMGDB73eJkvm+BedLkPFxX0ocyIaBcqciKbi4j4sOl0m2KEZ41dW6WdPaf3vniVb/P1XRVAmvg89r/qYvDkr8vbJdTwR34ddqyJY3DjY0RljTHBlyRb/3mOxLL1QiU/zXseyonvoVKwsS4lk77HYYIdmTJYlIjz77LMXl4cOHcqgQYMCGkN0dDR9+vRJvaAHd9xxR6oPinbu3JmSJUty5oxrioFDhw7hi1EFRo8eTdGiRalWrRrVqlVj+PDh6d5nSrJkxV8iIpzbcqznkyNbuPVAafbmPUahcu9SrOg/qW9sTDYxZXUM9d6eR+SAGdR7ex5TVqdv5sVcuXIxadIkDh1K61w5vnHu3Dlq1arFsGHD/HqckJAQRo4c6fP9Pvjgg6xZs4Y1a9bQvXt3n+/fXZas+N+pcYyPw4bxZHwf5hx+gjK77qEIJzhT+ANGLXyRC+7ze29fCIveD16wxgTBlNUx/N+kdcQci0WBmGOx/N+kdemq/ENDQ+nZsyfvvffeZZ917tyZCRMmXFy+6irX/Ozz58/n9ttv54EHHuDGG29kwIABjBkzhjp16lC5cmW2bdsGwMGDB2ndujW1a9emdu3aLF68GIBBgwbRs2dPGjduTMeOHZk/fz733XcfAKdOnaJLly5UrlyZKlWqMHHiRAB69epFrVq1qFSpEi+//HKav2ffvn157733OHfuXJq3zSiyZMVfP88uNtQbxq78tRDgUM6mPFdqIHfmLML/tk/hyekPczTuqKvSH98ZStYIdsjGBNSQWZuJjU889XBs/HmGzNqcrv0+8cQTjBkzhuPHj3u9zR9//MEHH3zAunXr+Oabb9iyZQvLly+ne/fufPjhhwA89dRTPP3006xYsYKJEycmahGvXLmSqVOn8t133yXa72uvvUaBAgVYt24da9eu5c477wTgjTfeIDo6mrVr17JgwQLWrl2bpu9YpkwZ6tevzzfffJNiuQYNGlxM3bi/5s6d67H8xIkTqVKlCm3atGH37qQze/pWlry5S/2+1IfLbuQ20of4fsnrDPnrB9pOaMKQg0eo3nY0REYFI0pjgia5+13pvQ+WP39+OnbsyLBhwwgPD/dqm9q1a1O8eHEArrvuOho3dv2PW7lyZX799VcA5s6dy4YNGy5uc+LECU6ePAlA8+bNPR5r7ty5fP/99xeXCxYsCMC4ceP44osvOHfuHPv27WPDhg1UqVIlTd/zhRdeoHnz5tx7773Jlvntt9+83l+zZs1o3749uXLl4rPPPqNTp07MmzcvTTGlRdas+JMhIrSv9yJVTx7juT0z6FIoD71PbaGr1ieHZMkfP8Z4VCIinBgPlXyJCO8q65T07duXGjVq0KVLl4vrQkNDuXDBlWJVVc6ePXvxs1y5cl18nyNHjovLOXLkuJhOuXDhAkuXLvVYwefNm9djHKp6WTfH7du3M3ToUFasWEHBggXp3LnzFT3pfP3111OtWjXGjRuXbJkGDRpc/OPkbujQoTRs2DDRusKFC19836NHD/r375/mmNIi+9V22xdSce1kxpVtR8O4eD5Y9QGP//I4R+KOBDsyYwKmX5PyhIeFJFoXHhZCvybl073vQoUK8cADDzBixIiL68qWLcvKlSsBmDp1KvHx8WnaZ+PGjfnoo48uLq9ZsybN2xw9epQTJ06QN29eChQowP79+5k5c6bHbTt27Mjy5ctT3P/AgQMZOnRosp//9ttvF2/Wur+SVvoA+/btu/h+2rRp3HTTTal9vXTJXhV/Qk6/7WiuajiIIY0+48UTZ1ixdxltJ9zNyjUjLy9vN35NFtSyeknealWZkhHhCFAyIpy3WlWmZfWSPtn/s88+m6h3T48ePViwYAF16tRh2bJlybbSkzNs2DCio6OpUqUKFStW5LPPPkt1m//+978cPXqUm2++mapVq/Lrr79StWpVqlevTqVKlejatSv16tXzuO3atWsvpp+SU6lSJWrU8M39wWHDhlGpUiWqVq3KsGHDGD16tE/2m5xMMedurVq11CcTsSx633Uj1z2nv30hm/6ew3OHl7D79H56RzanW9Tr5Nix6OIfCbsHYDKDjRs3+r2lmB2cOHGCbt26MX78+GCHkiae/vuLyEpVrZW0bLbK8XscpycyigqRUfwQ/y+vzOnNsB0/Ev3PCt7ctY3CVukbk+3kz58/01X6aZW9Uj0pyBuWl3fuHsnLBWuz8vRe2l5TmBVe9kowxpjMxCp+N7LjN9psms+Y4veQ9+xpus/qxmd/fMb5C+dT39gYYzIJq/gTuN34Ld90CD/c/h73xMbz8ZqPefTHBzkUeyhxWbvpa4zJpKziTxCzKtGN3DzXN+bNxp/xap4KrDm6ibZTWrBs3zJ72tcYk+llr5u7KfFw41fK3c795W7n5j+/59llr9FjdnceO3WWR9uMJMRu+hpjMilr8Xvhhpvb8f21D9Ds5Ck+vSonPf/6hoOnvZ/txpjswIZlvnIpDcv81VdfccMNN3DDDTfw1VdfpftYYBW/d7YvJM+qr3ijYndeO36GtQdW0+bHNizduzTYkRlzZRa970pbukvnvSsbljl9PA3LfOTIEV555RWWLVvG8uXLeeWVVzh69Gi6j2UVf2rcbvpy50BaNhvO9weOUzD+DI/OeZQPV3/IuQvnLpW1m74mMyhZw3VdJ1T+Prh3ZcMy+96sWbNo1KgRhQoVomDBgjRq1Iiff/453fu1ij81SW76EhnFda1G8V2+WrSIPcsXa7+gx+weHNg0zW76mswjMsp1XY/vDPPe8NlT6jYs8yW+GJY5JiaG0qVLXyxTqlQpYmLSN2EO+Pnmrog8DXQHFFgHdFHVOOez54AhQFFVDc5vQ28k87RvnsgoXtu+kDo/duc1WU3bfSt4867nqWc3fU1mERkFtbrBwsEQ9bxPnlK3YZkv8cWwzJ6G1EnrxOqe+K3FLyIlgT5ALVW9GQgB2jmflQYaAbv8dfyAiIyi2c2d+H73bgrliuCxPz/mg1UfXEr9GJORbV8I0SNclX70iMtz/leob9++jBgxgn///ffiOl8Ny5yQA4+JiSFfvnzAlQ3L/Msvv7B27Vruvfdevw7L7G2Lv3Dhwhe/d48ePS6OZFqqVKlEk7Ls2bOHEiVKpDnepPyd6gkFwkUkFMgD7HXWvwc8j+uXQObl/I9T7rZn+G7PXlpfU5fh64bTbVY3/vnX5vc1GViSe1cX0z4+qPxtWGYXXwzL3KRJE2bPns3Ro0c5evQos2fPpkmTJinG5Q2/VfyqGgMMxdWq3wccV9XZItIciFHVP1LaXkR6iki0iEQfPJgBu04m+R8nvO1oBv0xh7fzVmTT4fW0/bEtv+357VJZu+lrMhIP965oO9q13gdsWOa0SW5Y5kKFCvHiiy9evKn90ksvUahQoXQfz2/DMotIQWAi8CBwDBgPTAKeABqr6nER2YErFZRijt9nwzL7UjJDPLNuItv/ms5z197Aln9j6Fq6Eb2jpxBmI30aP7NhmX0jOwzL7M9UT0Ngu6oeVNV4XJV+FyAS+MOp9EsBq0TkGj/G4R/1+15ekUdGQfMPiGw1ijF/b6Ft3usYuXsOXW+ozD/FbgxOnMaYNLFhmdNnF3CriOQR1x2Wu4BJqlpMVcuqallgD1BDVbNWQjwyity1uvHSn78yuPBtbDn9D21+bMOC3QuCHZkxxvg1x78MmACswtWVMwfwhb+Ol6G49Za4e8NcxlXvR3ENofe83rwb/S7xF+IvlbPcvzEmwPzaq0dVX1bVCqp6s6p2UNUzST4vm6H78F8JD70lrp3+PN+Wvp8HT8czev1oOv/cmb0bJ9sDX8aYoLAnd30tmd4SuYD/Nv2Cocfi2HZoA21/f5Ff73zGbvgaYwLOKn5fS+6mr7O+SeXOjN+1i5K5CtFn/ecMXjGY+PNp69NsjDHpYRV/IDm5/9J1n+HbPXt4qMTtfLPhGzr93ImYU+kff8OYjGDy5MmICJs2bbrifSQd1M2TN998M9Fy3bp1r+hYgwYNSvFBrKzIKv5ASZL7z9l2NP+3+if+d1N3dhzZTNup9/PLrl8Sl7cbvyYTGjt2LPXr1080To4/JK34lyxZ4tfjZSVW8QdKMrn/RmeVH2oOpEzsKfr+2pe3l7/N2W2/2I1fkymdOnWKxYsXM2LEiIsV//z587njjjto06YNFSpU4OGHH744+Nirr75K7dq1ufnmm+nZs+dlg5L98ssv3H///ReX58yZQ6tWrRgwYACxsbFUq1aNhx9+GLg01DPA4MGDqVy5MlWrVmXAgAEAfPnll9SuXZuqVavSunVrTp8+7ddzkZHZ1IuBkswon0RGURr4OlcB3pv1BN9uHMPq+K8Zet9QStuNX3OF3ln+DpuOXHmqxZMKhSrQv07/FMtMmTKFpk2bcuONN1KoUCFWrXINAbF69WrWr19PiRIlqFevHosXL6Z+/fr07t2bl156CYAOHTowffp0mjVrdnF/d955J0888QQHDx6kaNGijBo1ii5dutCsWTM++ugjj2P2zJw5kylTprBs2TLy5MnDkSNHAGjVqhU9evQAXMM5jBgxgieffNIn5yazsRZ/BpHzurvoX6ED7+8/yO6cuXlgzRDm7JwT7LCMSZOxY8fSrl07ANq1a8fYsWMBqFOnDqVKlSJHjhxUq1aNHTt2APDrr79yyy23ULlyZebNm8f69esT7U9E6NChA99++y3Hjh1j6dKl3H333SnGMHfuXLp06UKePHkALo5t8+eff9KgQQMqV67MmDFjLjtWdmIt/ozCufF7V60nqbBqJP3K3cQz85+hfYnbee7O/5EzJOelcjGrPP+CMMaRWsvcHw4fPsy8efP4888/ERHOnz+PiHDPPfckGno5JCSEc+fOERcXx+OPP050dDSlS5dm0KBBHodITmjh586dm7Zt2xIamnK15Wk4ZnDdMJ4yZQpVq1Zl9OjRzJ8/P93fObOyFn9GkOTGb8nWo/hq6wY6RlRh7N4FPDKlJbtO7PLJ9HjG+MuECRPo2LEjO3fuZMeOHezevZvIyEgWLVrksXxCJV+kSBFOnTqVbC+eEiVKUKJECV5//XU6d+58cX1YWJjH4Z0bN27MyJEjL+bwE1I9J0+epHjx4sTHxzNmzJj0fNVMzyr+jMDDjd+wtqPpV7g2wyo9SsyJXTwwpQU/T+vmk+nxjPGHsWPHJroRC9C6devLpkRMEBERQY8ePahcuTItW7akdu3aye774YcfpnTp0lSsWPHiup49e1KlSpWLN3cTNG3alObNm1OrVi2qVat2savma6+9xi233EKjRo2oUKHClX7NLMFvwzL7UoYcljmA9s5+gX7bJ7A2dy4eLP8g/Wr3I1dIrtQ3NNlKVh6WuXfv3lSvXp1u3boFO5QMK6MMy2x8YftCSqz5ntHXd6DLv2f5YfMPPPLTI+w8sTPYkRkTEDVr1mTt2rU88sgjwQ4ly7CKPyNzy/2H3fUiz9z9JR8fjWPf0a08MK01P/39U+Ky9sCXyYJWrlzJwoULE90gNuljFX9G5iH3H9ViBBPy16F87Gn6/9afV5a+QtzWuXbT1wBc9gCUyR7S+t/dcvyZVPy2X/l41mOMyJuTG+PPM/S2V4is1DbYYZkg2r59O/ny5aNw4cIeuzOarElVOXz4MCdPniQyMjLRZ8nl+K0ffyYVdt1/6FuhIzVXDGNgidI8uGYIL4WHc1+5+4IdmgmSUqVKsWfPHg4ePBjsUEyA5c6dm1KlSnld3lr8mVVC/r9WN/avGsnz11Vm1YlttLqmLgPuep/w0PBL5eyBL2OyJevVk5UkeeDr6tajGPHXWnoUrM7kfYt5aHJL/j72tz3wZYzxyFr8mdGi912VufuDXE7Lfkl4OP+39kNiQ0L577HTNG823B74Miabshx/VpLCSJ91gfFHdtN/2w8MLJCb5TGzeaFULfKE5Ql4mMaYjMlSPVnN9oUUW/0dX97QkUdPnWXatmk8NOMhth7dGuzIjDEZhFX8WYlb7j/0rhfpfc+XfH70NEdP7KL99AeY/NfkS/197YEvY7Itq/izEg8PfN3WYiQT8tehSlwcLy15iYGLBnJ662y76WtMNmY3d7OJ83/P54uZj/LpVTkpe+4CQ2/5Lzfe3D7YYRlj/OiKunOKSHER6SsiE0VkqYjME5FhItJE7NHATCWk3B30uqkjX+7bz8mceXhozbtM+muSPeJvTDaUbMUvIl8C3zplPgC6AM8Ai4CWwGIRqR+III0PODN83VLnKcbvP0z1fJG8vORl/m/S/fz71+zLy1r+35gsK6XunB+p6h8e1q8BxolIbqCMf8IyPuX+wFdkFEUiG/DZ+M4Mr/Mgn+ycwfqFTzM0biDlKz+UuKwxJkuyHH92kMIDXyuur0f/X5/m+JljDChUmzabFiA2y5cxWUJyOf5kK34RuQ8oraqfOsuLgaLOxwNUdZK/gk3KKn7/Ohx7mBemtWNJ3D80DS/DwOZjiMgdEeywjDHpdCU3dwcAbjN9cBXQAGgKPOHb8EwwFf5nPZ/u3MZTBaow9/ROWky6m593/Gw3fo3JolKq+HOpqvv8fktUdb+q/g3Y8/9ZhZPTX1JtKKM2daPszmYUPXWMfgv60XdGBw6ePpi4rN30NSbTS6niL+i+oKq93BaL+SccE3Axq1hUbSg9fstDzLFYVsfWJ3RHJ+48XJjFh9bQYvK9rid+/15gD30Zk0WklOMfC8xR1ZFJ1ncDGqvqgwGID7Acv7/Ve3seMcdiL1t/V8GVnCv0PatyhnDbmfO83OANSt50fxAiNMZciSsZnfNpYKqItAdWOetqAvmBFr4P0QTLXg+VPsC8ozXZVisH49Z8yntFi3H/qrd4itO0r9CeHGKjfRiTWaX0f29OVb0FGAL847wGq2odVd0XkOhMQJSICPe4/r58W8mxciTtqj/O5AMnqJGvLG8vf5vOU+5n+/Htlwpa7t+YTCWlin+y8+9zqvqe85qdQvnLiMjTIrJeRP4UkbEikltEhojIJhFZKyKTRcT6DQZZvyblCQ8LSbTu9rCNDJX3Ls7yVaL1KD7dFM0bV/+Hbce20WZqK4avG078tl8t929MJpNSjn8NMB54DFerPxFVHZbijkVK4hreoaKqxorIOFzdQ/cC81T1nIi84+yrf0r7shy//01ZHcOQWZvZeyyWEhHhfF5uETfXvsPjQ1+Hil7Pm/P7MSd3KDfFn+fVWwZSofLDQYvdGOPZleT42wOtnDJFUyiXklAgXETicXUB3ZvkV8PvQJsr3LfxoZbVS9Kyekm3NXdeXsiZ5asI8L+YP5gb/SGvFy9N+9VD6RJ/mEerPkqukFyBCtkYc4WSrfhVdSPwhoisVdUf07pjVY0RkaHALiAWmO0hVdQV+CGt+zZB5gz41rDWk9ReOYIhle/iy3VfMnfrNF69YyjVilW7VC5mleepIo0xQeNN14xlIvK5iEwHEJGKItI5tY1EpCCu3j+RQAkgr4g84vb5QOAcMCaZ7XuKSLSIRB88eNBTERMM7oO43TmQAm1G8/raX/isxD3EndpHx5kdeWf5O5zeOsdy/8ZkUN5U/KOABUBpZ/kv4FkvtmsIbFfVg6oaD0wC6gKISCfgPuBhTeYmg6p+oaq1VLVW0aJXmmkyPudhli/ajqZeeHEmN3iPB0+f5duN39JqYV+WNhxgg70ZkwF5U/EXU9XvgAsATiV+3ovtdgG3ikgeZ9KWu4CNItIU6A80V9XTVxi3CZb6fS+vzCOjoH5f8t7QmIHlOzB6735CcxWg57oPeXnJy5w4eyI4sRpjPPKm4v9XRAoBCiAitYGTqW2kqsuACbge/lrnHOsL4CMgHzBHRNaIyGdXGLvJaJzcf81bnmJCzD90Ld2IqVun0nJcI+at+PDystb335igSHU8fhGphWsGrkrAH0BJoI2qrvF/eC7WnTMTSDLZS8Ly+qav8tKWb9nybwxNi9ZkwH/epfA/6xOXNcb4RZrH40+ycU7gJkCADap61vchJs8q/kwghcle4m97gpG/vcTnO34kb0gu+h89wb33DUfK3R68eI3JBq5kIpbmKe1QVaf5KLZUWcWfNWyb9Twv7ZjC2ty5iCoVxYu3vsg1ea8JdljGZFlXMhFLW+fVC/gGV5/7bsDXzr/GeG/7Qq77YwJfX9+B50+eYcXe37l/6v2M3zKeC3oh2NEZk60kW/GragdV7QDE4xp2oaWqtsCV6z8XqABNFuCW/w+560U63DuciQeOUUnDeHXpq3Sf3Z3dJ3ZfKms3fY3xK2969ZRT1Ri35b1AeT/FY7IiD33/S7caxZdXVWPQiTNsPLiOVtNa8dVvL3PeHvoyxu9SGqsnwUIRmQGMxdWlsx2w0K9RmazF05ANkVFIZBStty+k/sQuvF76eob+PYlZ11fk1YgSXB/4KI3JNrxp8T8BjAZuAW7FleO3ydaNb0RGcXWNrgzb+DuDC9/GnrMnaDu9LZ/+8Snx5+ODHZ0xWVKqFb+6jFfVJ53X+OSGWTAmzZyHviTqee7eMJcp1fvTKE8ZPlnzCQ/OeJD1h9ZfKme5f2N8ItmKX0R+FZFeIlIiyfpQEYkSkREi0sX/IZosK/yQwwUAACAASURBVMmAb7QdTaGpTzK4eCM+PBrH8X8P8NBPD/G/X54mznL/xvhMSi3+e4EwYLKI7HFmzPoL+BvoAnyqqqMCEaTJopIZ8I0L57ijxQgm79nL/XnKMmrPXNqUjSQ6PE8wozUmy/D2yd1cQDEgVlUP+T2qJOwBrmxq3huwcDDL6nTi5dgtxJyK4cHyD/J0zafJG5Y32NEZk+FdyQNcF6nqGVXdHYxK32RTTu6fqOe55c/pTKr6LI8UqMS4zeNoObUli2IWXSpnuX9j0sSrit+YgPKQ+88z6VH6F6vPN0dOk0eh19xeDJzZneMTOlvu35g0sorfZDwp5P6rthzJ+O3b6Jm/Ij/t/53mJa9htsQFM1pjMh1vc/ylgBtU9Vcn3x+qqv/6PTqH5fhNIk7uf/OtPXjp3G42HN5AwzINeeGWFyiax2ZrMybBFef4RaQrMA0Y7qy6Fpjq2/CM8ZJb7r/82kmMqfgYTxeqzcLd82kxtQVTtk5BVS33b0wKvEn19MH1xO4JAFXdgquHjzGB5SH3HzqhG10L12DiwRPckLsYLy5+kV4/PsjeiV0s929MMryp+OPcJ14RkRBcE7IYE1gp5P7LthrFqL/W8kJETVYf3kDLqyP4Lm6PDflsjAfeVPyLReR5ILeI/Af4AZju37CM8SCFid6JjCJHrW60Xz2ZySWbU+Oa2ry1/C26/NyF7ce3BydeYzIobyr+53FNrr4JeAr4BRjoz6CMSTO33H+JNd/z6XXteb3e62w99CdtprZixLoRnLtw7lJZy/+bbCzFit9J64xU1U9V9X5nMpZPVe33s8lAPOT+ZUIXWoQUZGqdQUTFxvH+qvd5aMZDbF73naus5f9NNpZixa+q54HiIhIWoHiMSbvkcv8xqyhSvhnvNf6c/x2L48Cx7bRb+SYf1m7N2TK3BjNiY4Iq1X78IvIZUA1XF86LffdVdZh/Q7vE+vGbdJv3BscXDWXwTfWZdnoH5QqU49V6r1K1aNVgR2aM36RnrJ6DwBwgD1DU7WVM5uDk/wvUf443tq/nk5sf5/S/++nwUwfeWf4Op+NPXypnuX+TDaQ69aKqvhiIQIzxC/f8f2QURDagwfjOTKnbm/fWfcG3G79l/u75DIpsxS1z33KVMyaL8ybVMwfXXLuJqGpjfwWVlKV6zBVb9L7rRq57N9DtC133BUrWIHpKVwYVK8bOcydpfU1dnv3PUPLlzBe8eI3xoeRSPd5U/Le4LeYGWgNnVLWfb0NMnlX8xm/mvUHcb0P4pHIjvjq1hSI5cvNihY7cUcttWumEPxSeJo03JgO74hy/qi5zey1Q1T5AHb9EaUwgObn/3A368czWlXxX7TkichfkyfWf8fxPXTgSd+RSqsi6f5osJNUcv4jkd1vMAdQEivstImMCwUPuv9L4znzfejgjds/h8x0/8vv4xgw4epy724xCkj4xbEwm5k2vnvXAn86/q3E9tdvDn0EZ43fJ9P0P27eWx+54i/HF76H06RP0L5CbPtsnsP/f/UEN1xhf8ibHH6aq8UnWharqOb9G5sZy/CagnF8D52t2YczGb/iwwFWEAs9e15bWtw1ARC6Vs9y/ycDS049/mYd1y9MfkjEZkFsKKOSuF+l473AmHThOxVxFeeWv7+g+tQ27T+y23L/J1JLN8YtIMVy5/HARqcyloZjz43qYy5isx0MKqHSrUQzfs5KJIXG8u3kMraY048kTp3m4zUhCLPdvMqGUbu7eC3QFSgGfuK0/CdhDXSZr8pS2iYxCIqNoA9Q/dYrXt/3AkHzhzNr0Ja8WKs11EdcFPExj0iPZVI+qjlLVBkA3VW3g9rpHVccHMEZjMobtC7lmzVg+vLEzbx+PY9fRbbSZ2oohc5/ixNkTicrZ0A8mI/NmyIZxItIEqITrAa6E9W/6MzBjMhS33L9ERnFvuShundCZD6+vyTd7fmHa+GU8UbMvbXJeTeiEbjb0g8nQvOnV8wkQAUQBo3A9ufu7qnZNdeciTwPdcQ35sA7oguv+wA9AWWAH8ICqHk1pP9arxwRdCkM/bCpQlCG/v8HynCFcd+4Cz1XrTf2avYIXqzGO9AzZsFZVq4jIH6paVUTyARNTG6tHREoCi4CKqhorIuOAn4CKwBFVfVtEBgAFVbV/Svuyit9kdPrL6/y68iPeLRnJrnOnqF+yPv1q9aNcRLlgh2aysfR054xL+FdErnGWy3p53FBcvYJCcbX09wItgK+cz78CWnq5L2Mypu0LkZUjubPWk0yJ2c9z5Vrxx4E/aDW1JW/MeYKjcUcTlbX8vwm2VHP8wE8iEgEMBdYA57lUcSdLVWNEZCiwC4gFZqvqbBG5WlX3OWX2Od1GLyMiPYGeAGXKlPHqyxgTKFNWxzBk1mbKnIjmk5wfsr7uB9S/sxVhkQ3oNL4zzVt+xMc7ZjA+ZgEzJjbhsWq9aZ+7FGETu1v+3wRdiqkeEckB1FbVZc5yOBCuqkdS3bFIQWAi8CBwDBgPTAA+UtUIt3JHVbVgSvuyVI/JSKasjuH/Jq0jNv48j4b8yFotx5qQKrzVqjItq5dM9ETvtj9/YMiSQSzOFcq15y7wTJXH+E+t3pee/jXGj64o1eNMqv6B23KsN5W+oyGwXVUPOkM+TALqAvtFpLgTVHHggJf7MyZDGDJrM7Hx5wH4/Hwzll6oRGz8eYbM2uwqEBl18XmA625+kM+uf4RP/jlASO4IntrwBT1m92Dzkc3BCt8Yr3L8c0SkxRXsexdwq4jkEVfz5i5gIzAN6OSU6YRrLl9jMo29x2K9X+8M/dygdh8mxOzjhesfYPOBNbT9sS2DlgziUOyhS+Us928CxJscf2+ggIicwZWrF0BVtVBKG6nqMhGZAKwCzuEa2fML4CpgnIh0w/XHoW064jcm4EpEhBPjoZIvERGeeEWSoZ/DIhvQfnxn7rntCT5f9yVjt07m5x0/071kQzos+45clvs3AeJNi78IEIarwi7qLHs12bqqvqyqFVT1ZlXtoKpnVPWwqt6lqjc4/3qbOjImQ+jXpDzhYSGJ1oWHhdCvSfnEBZMZ+rmAhPD83V8w5eAp6oRG8MGOqbQoU4ZZEktq3auN8YVU+/EDiEg7oJyqvikipYCrVXWl36Nz2M1dk9Ek9OrZeyyWEhHh9GtS3nVjNy3mvQELB7OsTicGn4thy9Et1Mh/Hc83eINKRSq5ytjQzyYd0vMA10e4WvxRqnqTiBQCZqlqbf+Eejmr+E2Wk5AGqtUNokdwvs1IpvyzlGGbxnAkRGh+XXP6FLmFq6f1TfyrwZg0SM8DXHVV9VGcB7mc1ExOH8dnTPbhnvu/c6Br7P8JXWldvB4zGrxLt3/P8vO26TT7/QU+rdOW2NIBa2OZbMKbij/e6c+vACJSGLjg16iMycqSyf0Ts4qrbmhC3wodmbp7Nw3ylOGTnT9x3+T7+HHbj1xQ+9/O+IY3Ff/HuB7EKioir+Aaf+cdv0ZlTFZWv+/lqZuEvv9O989SdZ/l3R2b+Krq0xQ9f4EXFr3AwzMeZs2BNa7y1v3TpIM3wzJ/LSIrcT2QBdBWVf/0b1jGZENJun8S2YAa4zvzXb2nmb7yQz7IuZsOMzvQtGhNnt7wGyVajwp2xCaT8qbFDxACxANn07CNMSYtkkkB5dDzNG82nB9376VX/krMPxBNs2L5+eDoav6N/zeoIZvMyZtePQOBh4DJuB7eagGMUdW3/B+ei/XqMYaL3T//qdubD8IvMP3v6RQOCadPuVa0uLUfITmcZwusC6hxpKdXzyO4Bmr7r6oOBOoAHX0doDEmBU7un6jnuWbNWN4qdTff3fMdpfMW5+W/xtBu0n0s37f8UrqoZI1gR2wyMG8q/p0kvhcQCvztn3CMMZfx0P2T8Z2pfOooX7ecwpAKXTlxYhfdZnfjqdmPsuu+wdbv36TIm4r/NLBeRIaLyJe4plA8JiL/E5H/+Tc8Y0xK3T9FhKa3PM3UMg/w1JFj/B6emxYr32DoiqGJJ4A3xo03Of5uKX2uqiN8GpEHluM3JgVuTwEfWjWSDys3ZPI/S4gIu4rHazxFmxvbEJoj1HL/2dAVD9mQEVjFb0wyknYBdZY31e7I4K3jWJEzhOsKXEe/0k2p98tgG/4hm7nim7si0lREVojIARE5IiJHRcRG1DQmI0gmDVQhLIIRd33K+8fiOPvvfh7782N6VajN3wVLBTVckzF4k+rZCjyAK7d/8ZlxVT3v39AusRa/MVdo3hucXTiYsdWa8/npbZyOP8UDJaJ4vMHrROR2ZkC1FFCWlZ7unHuANaoar6rnE16+D9EY41NOF9CcUc/TactiptccSJviDRgXM597Jjbm6/VfE79tnnX/zIa8afHXAV4G5gNnEtar6jC/RubGWvzGpFEyuX/ajmbrv3sZuuSVixPAP1u1F3fUfMImgM+C0tPifwU4D0Tgmnkr4WWMyahS6AJ6/c3tEk0A32f95zYBfDbjTYt/parWDFA8HlmL3xgfcuv+GR89gvF1O/PJ9mmcuHCWVje0onf13hQJL2K5/ywgPS3+X0TkTj/EZIwJtCRPAYe1Hc1DS0Yzo2x7Hvn3LFO3Tua+yfcxfMF/OWO5/yzLm4q/BzBXRE5Zd05jMrlUJoCffPAUtUNsAviszptUT4in9dad05gsyBkB9Pc6HRlybq9NAJ/JXXGqx6ng2wL9nffFgWq+D9EYE1RuI4De+ucMxt38JC/f8BA7jm6l3Yx2DFw0kP2bpln3zyzAmyd3PwL+A3RwVp0GPvNnUMaYAEtmAvg2NgF8luRNjr+uqj4KxAGo6hEgp1+jMsYElk0An614U/HHi0gOQAFEpDBuQzcYY7KAtE4AH17UNQH85OaXJoAHmwQ+k0i24heRhMlXPgYmAkVF5BVgEfBOAGIzxgSbhxRQjdmv813Fx3ijfEcOHNtOh5kd6LegH3s3Trb8fyaRbK8eEVmlqjWc95WAhrjm3J2rqn8GLkTr1WNM0Cx631WRu/8acOvVc3rrHEbP6s2ofOFcuHCOjmWa0L3Bq+QNyxu8mM1FaR6PX0RWq2p1v0fmBav4jcnA5r3BP4vf5YMKdZl+eqdNAJ+BJFfxh3oq7CgqIs8k96Gq2rSLxmR3Tv7/mnrP8lb0CB5qMojBf0/k5b/GMHbvAvrVe4U6cXGX0kUmQ0ip4g8BrsKV3jHGmMSSjgAa2YDK4zvzdZtRzDq4ivfWD6fb7G7cGXeOZ+8bTBmb+SvD8CrHH2yW6jEmA0ol/x839xW+XTecLwsX4azAw/kq0LNSJ/Lf0NRjeeN7V/LkrrX0jTHJS6ULaO5VX9G96mPMOHCC5sXq8PXx9dy76DlGLPwvp+NPX/rFYL2AAi6lFn8h52GtoLMWvzGZSHITwN/9Gu9vn8bioxuIyJGTjsdP0b7JMK66oUmwI86y0tzizyiVvjEmk0luAvjjB/ms+Q+MKdaIKqeOMyxfTppEv8rnf3zOybMngxpydpPq6JxXvGOR8sAPbqvKAS/hmsLxMyA3cA54XFWXp7Qva/Ebk0W4TQKzfs0oPruhDvOPrCNfaDgdKnXh4YoPkz9nfsv9+0ia+/H7+OAhQAxwC/Al8J6qzhSRe4DnVfWOlLa3it+YLCCZFNDGWh357K8fmJc7lKvCruLh4g3oED2BAm1GX34PwaRJembg8oW7gG2quhPXmD/5nfUFgL0BisEYE0zJpIBuyhnBB40/Z8LhWG4LjeDzXTNpUrwww46u4VjcsaCGnFUFqsU/Elilqh+JyE3ALFy9hnLgGv1zp4dtegI9AcqUKVNz587LihhjshJnEpgtt/bki3Bl9o5ZhIfkov1Nj9CxUkcK5S5kKaA0ClqqR0Ry4mrVV1LV/SIyDFigqhNF5AGgp6o2TGkfluoxJnObsjqGIbM2s/dYLCUiwunXpDwtq5e8VMAt90/0CGg7mq3/7uWLxa/wc+4wcofmpt019ei0aiqFLQXktWCmeu7G1drf7yx3AiY578cDdQIQgzEmSKasjuH/Jq0j5lgsCsQci+X/Jq1jyuoYVwEPI4AyvjPX5y3B4MafM+XQv9yZsxhf7Z5D02siGHJoGYdiDwXvC2UBgaj42wNj3Zb3Arc77+8E/gpADMaYIBkyazOx8Ymn6I6NP8+QWZtdCylMAkNkFOWqd+Xt9b8xtfi9NI68m283fkvT8Q15Z24fDp4+eGmnNheA1/ya6hGRPMBuoJyqHnfW1Qc+wDVOUByu7pwrU9qPpXqMybwiB8zAUy0jwPa37015Yw8poF2Fy/LFkteY/s9SQnKE0ab8A3QtUImrp/VN/AfEXNHonOmmqqeBwknWLQJq+vO4xpiMo0REODHHYj2uT5GHQeAY35kybUfzetMveXTDJIYvfIFxm8YyXi/QqtrddC92I9f45VtkLYHqzmmMyab6NSlPeFhIonXhYSH0a1I+5Q1TSgEBpSu24pUbOzB99x5aXHU9E/cv4e4JTXhtdi/2nnLrJW4poMsEpDtnelmqx5jMLdVePVciSRpob7N3Gbl7DpP2LkRzhNLi+pZ0j6hCqenPZdsUUFCf3E0vq/iNMYkk8xQwbUfzT9xRRszvz8Q8uVA9T7Nr6tKj7ouUzl86yEEHXrCf3DXGGN9JIQ10zU0tGFi+AzN37ebBfOX56dAqmk2+l4E/d2fnCbcHQbNxCsha/MaYrCVJCuhg8w8YtWcu4/f8wlnJwb3l7qNHoepEzuif5VNAQenVY4wxAeWhJ1DR8Z15vu1oupZuyOhf+jFu+0/M2DaNpjffzqMFS1Eu2DEHgaV6jDFZRwopoCLlm/FchQ7M3LGTTvlv4tej62k5tQX9furM1qNbL+0jG6SALNVjjMkekqSAjrb8iK/3/MJ3O2dxOofQ6NpGPFakNjf+NDDLpIAs1WOMyb48pIAKju/MU21H06nUXXw952m+2/Urc3bOoWGlejxaoBgVgh2zH1mqxxiT9aWQAoq48R763NSRWTu281j+Siw7vo22P7alz7w+bDi8IZhR+42leowx2VuSFNCJ+z9hzMaxfHNyIyfPxXJ7qdvpVbUXlU4eznRzAVg/fmOMScrDkND5Jz9OryK1mbX3EE+Wbcaag2toN6Mdj895lLX5CgY7Yp+wit8Yk30llwK6cI58bUbTc9n3zCrSkKdOnmVd3gI8vGYoj815jDUH1gQz6nSzVI8xxiTHmQ6SqOc53eBpvt/8PV+t/pQjF+K4tfit9KraixpX18iwU0JaqscYY9Ji+0LXHABRz0P0CPLsiabrzV2ZedubPHfyDFsOrafTz53oNqU1K6Z0hZI1gh2x16ziN8aYpJKZDpLtC8lzfSM63Tucn/f8w/MR1fn7yCa6Fgqny+ZRLN+3nMyQRbGK3xhjkkplLgAiowiv1Y0Oq6cys8wDDKgzgF2HNtBtdjc6/9yZpXuXuv4AZNCngC3Hb4wxaeVhSsgzF+KZ9FMvRhS5mv1nj1Etfzke27mBus2HI+VuT3WX/mA5fmOM8YVk0kC5coTR/r7h/LQ7hhcL1uKfo1t5rGBuHln+CgtXfpI4BRTkXwJW8RtjTFqklAaKjCJnrW48sGoSP137IC/d9hKHOM8Tf37KQ5ObsWD3AvTvBa4/HEG8GWypHmOM8RUPKaD4Mrfx47J3+WLjN8SE5uCm+PM8VrUX/6nVGxHxaziW6jHGGH9KJgUUtmspreoO4Mey7Xjt4GFOhRfgqQ1f0Pb7O5i7YhgX9ELifQQgBWQVvzHG+EJKKaDtCwlbOYqWNZ5gWsx+3ijfkbiQUJ7e8CVtJt7NrB2zuPD3/IClgCzVY4wx/pTMxPDn24xk5oFoPl8/ih2hObj+3AUevbk7jWr3ISRHiE8ObakeY4wJhmR+CYTsXcN9tz7HlLLtGHzgEJo7gn6bRtJqWitmzHyS83/PT7wfH6aBrOI3xhh/qt/38tm8IqNc67cvJGTlKO6u2ZtJMfsYclNXckgOBhyYT8tfn+DH34dw7sK5S78afJQGshm4jDEmGJKkgHJENqDp+M40bjOSX0Li+WzFu7yw+WvCDm2l6YY5Pp0O0lr8xhgTDMmkgHLsXUOjaxsxvvVPfFwkiobR37u6h/pwDmBr8RtjTDB4GsI5MupiBZ9jxyKi1s+8ODookQ2sxW+MMVlWCqOD+oJV/MYYk9GkNjpoOlmqxxhjMppU0kDpZS1+Y4zJZqziN8aYbMZSPcYYkwFNWR3DkFmb2XsslhIR4fRrUp6W1Uv6ZN9+a/GLSHkRWeP2OiEifZ3PnhSRzSKyXkQG+ysGY4zJjKasjuH/Jq0j5lgsCsQci+X/Jq1jyuoYn+zfby1+Vd0MVAMQkRAgBpgsIv8BWgBVVPWMiBTzVwzGGJMZDZm1mdj484nWxcafZ8iszT5p9Qcqx38XsE1VdwK9gLdV9QyAqh4IUAzGGJMp7D0Wm6b1aRWoir8dMNZ5fyPQQESWicgCEantaQMR6Ski0SISffDgwQCFaYwxwVciIjxN69PK7xW/iOQEmgPjnVWhQEHgVqAfME48zD+mql+oai1VrVW0aFF/h2mMMRlGvyblCQ9LPCZ/eFgI/ZqU98n+A9Gr525glarud5b3AJPUNQPMchG5ABQBrFlvjDFwMY/vr149gaj423MpzQMwBbgTmC8iNwI5gUMBiMMYYzKNltVL+qyiT8qvqR4RyQM0Aia5rR4JlBORP4HvgU6aGeZ/NMaYLMKvLX5VPQ0UTrLuLPCIP49rjDEmeTZkgzHGZDNW8RtjTDZjFb8xxmQzkhnuq4rIQWDnFW5ehIzZa8jiShuLK20srrTJqHFB+mK7VlUvexAqU1T86SEi0apaK9hxJGVxpY3FlTYWV9pk1LjAP7FZqscYY7IZq/iNMSabyQ4V/xfBDiAZFlfaWFxpY3GlTUaNC/wQW5bP8RtjjEksO7T4jTHGuLGK3xhjsplMXfGLSFNn7t6tIjLAw+ciIsOcz9eKSA1vt/VzXA878awVkSUiUtXtsx0iss6Zpzg6wHHdISLH3eZJfsnbbf0cVz+3mP4UkfMiUsj5zC/nS0RGisgBZzBBT58H69pKLa5gXVupxRWsayu1uAJ+bTn7Li0iv4rIRnHNPf6UhzL+u8ZUNVO+gBBgG1AO19DOfwAVk5S5B5gJCK6JX5Z5u62f46oLFHTe350Ql7O8AygSpPN1BzD9Srb1Z1xJyjcD5gXgfEUBNYA/k/k84NeWl3EF/NryMq6AX1vexBWMa8vZd3GghvM+H7AlkPVXZm7x1wG2qurf6hrx83tck7i7awF8rS6/AxEiUtzLbf0Wl6ouUdWjzuLvQCkfHTtdcflpW1/vO+n8Dn6hqguBIykUCca1lWpcQbq2vDlfyQnq+UoiINcWgKruU9VVzvuTwEYg6eD7frvGMnPFXxLY7ba8h8tPXHJlvNnWn3G564brr3oCBWaLyEoR6emjmNIS120i8oeIzBSRSmnc1p9xJczv0BSY6LbaX+crNcG4ttIqUNeWtwJ9bXktmNeWiJQFqgPLknzkt2ssEDNw+ctl8/Ti+g/lTRlvtr1SXu9bRP6D63/O+m6r66nqXhEpBswRkU1OqyUQca3CNbbHKRG5B9dsaTd4ua0/40rQDFisqu4tOH+dr9QE49ryWoCvLW8E49pKi6BcWyJyFa4/Nn1V9UTSjz1s4pNrLDO3+PcApd2WSwF7vSzjzbb+jAsRqQIMB1qo6uGE9aq61/n3ADAZ18+6gMSlqidU9ZTz/icgTESKeLOtP+Ny044kP8X9eL5SE4xryytBuLZSFaRrKy0Cfm2JSBiuSn+Mqk7yUMR/15g/blwE4oXr18rfQCSXbnBUSlLmXhLfHFnu7bZ+jqsMsBWom2R9XiCf2/slQNMAxnUNlx7qqwPscs5dUM+XU64Arlxt3kCcL2efZUn+ZmXAry0v4wr4teVlXAG/tryJK4jXlgBfA++nUMZv11imTfWo6jkR6Q3MwnWXe6SqrheRx5zPPwN+wnVnfCtwGuiS0rYBjOslXFNSfiIiAOfUNfre1cBkZ10o8J2q/hzAuNoAvUTkHBALtFPXlRbs8wVwPzBbVf9129xv50tExuLqiVJERPYALwNhbjEF/NryMq6AX1texhXwa8vLuCDA15ajHtABWCcia5x1L+D6w+33a8yGbDDGmGwmM+f4jTHGXAGr+I0xJpuxit8YY7IZq/iNMSabsYrfGGOyGav4jc+JSGG3EQ//EZEYt+UlfjpmdREZ7ryvICJLReSMiDyXpFyKozUGm4h0EpG/nFenZMrMF5GATQwuIr1FpEugjmf8z7pzGr8SkUHAKVUd6ufjjAdeV9U/nEfsrwVaAkfdjy0iUcApXINf3ezPmNLKGQ44GqiF6xH8lUBNvTToWkK5+cBzqurroYJDVPW8h/V5cA1nUN2XxzPBYy1+E1Aicsr59w4RWSAi40Rki4i8La6x5JeLawz065xyRUVkooiscF71POwzH1BFVf8A1yP2qroCiE9aVr0YrVFEmonIMhFZLSJzReRqZ/0gEflKRGaLa6z2ViIy2In3Z+cR/PRoAsxR1SNOZT8H18BhnrR1ztUWEWngxJdbREY58ax2xutBRDqLyEdu32+6iNzhvD8lIq+KyDJcg6i9LSIbxDX++1AAVT0N7BCRgA3xYPzLKn4TTFWBp4DKuJ5ivFFV6+AaZ+ZJp8wHwHuqWhto7XyWVC3Al6mbRcCtTgv3e+B5t8+uw/UofQvgW+BXVa2M62nUe5PuSBJP9OH+GubhuGkZdTHUOVd9cT2NCvAEgBNPe+ArEcmdynfNi2s4g1uADbieYq2kqlWA193KRQMNUtmXySQy7ZANJktYoar7AERkGzDbWb8O+I/zviFQ0Xl0HiC/iORT1xjmCYoDB30YVyngB3GNfZ4T2O722UxVjReRdbgel094jH8drjFhElHVIcAQL4+bllEXEwb1Wul23PrAh85xN4nITuDG/40gHQAAAb5JREFUVI55nktDEZ8A4oDhIjIDmO5W7gBQIZV9mUzCWvwmmM64vb/gtnyBS42SHMBtqlrNeZVMUumDq7WdWss2LT4EPnJazo8m2fcZAFW9AMTrpZtk7jFflMYWf1pGXUw4V+fdjuvpDwfAORL/v+7+feIS8vqqeg7XAGoTcd0f+TnJNrHJ7N9kMlbxm4xuNtA7YUFEqnkosxG43ofHLADEOO899qzxlqoOcfuj5f7q46H4LKCxiBQUkYJAY2edtxYCDwOIyI24BvzajGsKwWoikkNESpPM8MLiGhu+gLqGTe4LuJ/rG/FtOs0EkVX8JqPrA9RybjZuAB5LWkBVNwEFnJu8iMg1zkiMzwD/FZE9IpLf+WwssBQo76zv5uGYg4DxIvIbcMgv38oDdU0C8hqwwnm9qoknBknNJ0CIk4b6AeisqmeAxbjSVeuAobgmRfEkHzBdRNYCC4Cn3T6rB8xNQywmA7PunCZLEJGngZOq6unmr0kHEakOPKOqHYIdi/ENa/GbrOJTEt8zML5TBHgx2EEY37EWvzHGZDPW4jfGmGzGKn5jjMlmrOI3xphsxip+Y4zJZqziN8aYbOb/AVK1sFA6BmxcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"Newton's Law of Cooling\")\n",
    "plt.ylabel(\"Temperature (defG)\")\n",
    "plt.xlabel(\"Time (11 am = 0 hours)\")\n",
    "plt.plot(times1, numerical_T1,'o', label = 'Numerical, N = 5')\n",
    "plt.plot(times2, numerical_T2, 'x', label = 'Numerical, N = 50')\n",
    "plt.plot(times1, analytical, label = 'Analytical')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`b.` The final temperature at $T \\rightarrow \\infty $ will be the ambient temperature $T_{a} = 65^{o} F$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`c.` The time of death is around 9:07 am."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 am Subtract: 53.19137310876379 min\n"
     ]
    }
   ],
   "source": [
    "T = 98.6\n",
    "Ta = 65\n",
    "To= 85\n",
    "\n",
    "t = np.log((T-Ta)/(To-Ta))/-K\n",
    "print((int(t)+11), 'am', 'Subtract:', -(t-int(t))*60, 'min')"
   ]
  },
  {
   "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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5xcdX3/8debXFgCSaDhYiCEBAigYkBcA4pYMIKiCIUCSgsN2DZWLYIX2qW09VL1F8X+Wrw3PwKGFtogcvFGIPIQKYqBhEaIEG4h1k2AbO5kN8nk8vn9cc4sk2R39sxmzs7Mzvv5eMxj5pw5c85nNvDZ737P9/v5KiIwM7PmsVetAzAzs4HlxG9m1mSc+M3MmowTv5lZk3HiNzNrMkNrHUAWBx54YEyYMKHWYZiZNZSFCxeuioiDdt3fEIl/woQJLFiwoNZhmJk1FEm/62m/u3rMzJqME7+ZWZNx4jczazIN0cffk61bt9Le3s7mzZtrHcqAaGlpYdy4cQwbNqzWoZhZg2vYxN/e3s7IkSOZMGECkmodTq4igtWrV9Pe3s7EiRNrHY6ZNbiG7erZvHkzY8aMGfRJH0ASY8aMaZq/bswsXw2b+IGmSPpFzfRdzSxfDZ34zcwGq45Xt/DVuUt4oWNj1c/txN9P69at49vf/nZVzrVmzRrOPPNMJk2axJlnnsnatWurcl4za1xLXt7Atx98gVc2VL+L14m/n/qT+AuFAp2dnbvtnzFjBlOnTuW5555j6tSpzJgxo1phmlmDenFVkiuOOmi/qp/bib+f2traeOGFFzjxxBO55ppryh779NNP8+lPf5pjjz2WZ599drf377nnHqZNmwbAtGnTuPvuu3OJ2cwax9KOTvYdPoSDR+5d9XM37HDOUp//0W95asWGqp7zDYeO4rMfeGOv78+YMYPFixezaNGiHt/v7Ozk9ttvZ9asWUQEV1xxBU888QQjR47c7dhXXnmFsWPHAjB27FhWrlxZnS9hZg1r6apOJh60by4DOwZF4q9HY8eOZfLkydx4440cd9xxtQ7HzBrM0o6NnDT+gFzOPSgSf7mWea3ccccdzJo1i/PPP59LLrmEadOmccQRR/R47CGHHMJLL73E2LFjeemllzj44IMHOFozqyebt25n+bpNXPiWcbmc3338/TRy5EheffXVXt8/66yzmDNnDg8//DCjR4/mvPPO493vfjfLli3b7dhzzz2X2bNnAzB79mzOO++8vMI2swbwu9VdRMDEA/fN5fxO/P00ZswYTj31VI4//viyN3fHjBnDVVddxaJFi/jyl7/MkCFDdjumra2NefPmMWnSJObNm0dbW1ueoZtZnVuajt3PY0QPDJKunlq57bbbKjp+ypQpPe4fM2YMDzzwQDVCMrNBYGk6lHOCW/xmZs1haUcnh4zam/32zqdt7sRvZlZnlq7ayJEH5tPNA078ZmZ158V0DH9enPjNzOrIms4C67q2cmRO/fvgxG9mVlfyHtEDOSd+SftLukPSEklPS3pbuv9KSc9I+q2kr+YZg5lZIymO6MlrDD/kP5zzBmBuRFwoaTgwQtIZwHnA5IjYIqkhp6muW7eO2267jY997GN7fK41a9bwwQ9+kGXLljFhwgRuv/12Djggn6naZo1sx46gsH1HrcPI1fMrNzJsiBh3wD65XSO3xC9pFPBO4HKAiCgABUkfBWZExJZ0f0NWJCuWZe4t8a9duzZz8i6WZW5ra2PGjBnMmDGDr3zlK9UM12xQ+LObHuXh51fVOozcTTp4P4YOya9DJs8W/5FAB3CzpBOAhcBVwDHAaZK+BGwGPhMRj+36YUnTgekA48ePzzHM/ikty3zmmWdy/fXX7/R+a2srJ598Mn/xF3/BGWecUbbC3j333MODDz4IJGWZTz/9dCd+s11s3xE8umwNbz9qDKdNOqjW4eTqrRPy/Ys/z8Q/FDgJuDIi5ku6AWhL9x8AnAK8Fbhd0pEREaUfjoiZwEyA1tbWnd7bzb1t8PKT1Y3+dW+Cs3tfEKWvsszPPvss9957L9/85jf5+Mc/zmWXXcbll1/OoYceutuxLsts1rcV6zZR2LaDc084lA9Nqb/GYCPJ8+ZuO9AeEfPT7TtIfhG0A3dG4lFgB3BgjnHUxJAhQzjnnHO48847eeihh1i6dCnjx4/n0UcfrXVoZg2puPbskTmOdmkWubX4I+JlSb+XdGxEPANMBZ4CXgDeBTwo6RhgOLBnnXZlWua1tH79eubMmcPNN9/MsGHDmDVrFpMnT97tOJdlNuvbiwMw2qVZ5D2q50rg1nREz1LgCqATuEnSYqAATNu1m6cR9FWW+dJLL+WRRx7hoosu4pZbbmHSpEm9Hlssy9zW1uayzGa9WNrRyciWoRy43/Bah9Lwck38EbEIaO3hrUvzvO5AKC3LfPbZZ+92c/fiiy/me9/7HkOH9v0jbmtr4+KLL2bWrFmMHz+e73//+3mFbdawkvo1+SxF2GxclnkPlCvLfO6552Y+j8sym/XtxY5OTj5yTK3DGBRcssHM6l5XYRsr1m/OtX5NM3HiN7O6131jN8eKlc2koRN/A94T7rdm+q5muyom/jxr1DeThk38LS0trF69uikSYkSwevVqWlpaah2KWU0s7SguRTiixpEMDg17c3fcuHG0t7fT0dFR61AGREtLC+PGjat1GGY1sbRjI4eObmHE8IZNWXWl7E8xHX//PuA04FBgE7AY+GlELMk/vN4NGzaMiRMn1jIEMxsgL67q9IzdKuq1q0fS3wPzgTOA3wCzgR+S/LL4F0lzJR0/IFGaWdOKCJZ2dHrGbhWVa/E/GRFf7OW9r0oaCxyeQ0xmZt06Nm7h1S3bONIjeqqmXOL/YbkPRsRLwEvVDcesn15+EhbdBk1wsx/g0WVreGG/k3h2/9NqHUruVm0sAK7RU03lEv9CkmqaSPrXiLh6YEIy64dHZ8Lj/w57j6p1JAPi9Vu2svCV7dxBc9znGnfAPkwet3+twxg0yiX+0oIY78w7ELM9smUjjDkarlxQ60gGxEjgo+nDrFLlxvE3x9/MNjgUOmG4x3ibZVGuxX+cpMdJWv7Hpq9JtyMiTso9OrOstnbBcA/3M8uiXOJ/04BFYbanChth38G9DqtZtfSa+CPihYEMxGyPFDph/yNqHYVZQ+hz/rOkteze378eWABcExHLcojLrDIFd/WYZZWl8MU3gFeA20j69z8EHAQ8D9xMMrPXrLYKG31z1yyjLIn/rIg4pWT725J+HRGnSPqbvAIzq8jWLhjuCT5mWWQqyyzpgl1eF8f478gjKLOKbCvA9gIMc+I3yyJL4r8U+EtJayStBv4SuEzSCMCzea32tia12t3iN8umz66eiHgeOLuXt39R3XDM+qHQlTw78Ztl0meLX9LRku6T9Jt0e7Kka/MPzSyjglv8ZpXI0tVzI/B5XuvPf5Kk+8esPrirx6wiWRL/vhHxq+JGJIvcbs1yckn7S7pD0hJJT0t6W8l7n5EUkg6sPGyzEsUW/zAP5zTLIstwztWSJpJO4pL0R8DLGc9/AzA3Ii5Ml3EckZ7jcOBM4H8rD9lsF919/J7AZZZFlsT/18AskqJtvyNZfOWSvj4kaRRJOefLASKiABTSt/8F+BvgnspDNttFYWPy7AlcZplkHdXzLkmjAUXEuoznPhLoAG6WdALJwi5XAVOB5RHxG0m9fljSdGA6wPjx4zNe0pqSb+6aVaTXxC/pE73sByAivp7h3CcBV0bEfEk3AJ8j+SvgrL4Ci4iZwEyA1tZWrw1gvdvqrh6zSpRr8Rdr3E4CpgA/SrfPIdv4/XagPSLmp9t3kCT+iUCxtT8OeFzSlIjIet/AbGfFrh7f3DXLpFxZ5n8AkHQfcGJEbEi3/wGY09eJI+JlSb+XdGxEPEPSxfN4REwtHiNpGdAaEav27GtYUyt0gYbA0L1rHYlZQ8hyc/cIYHPJ9hbIvMLzlcCt6YiepcAVlYVnlkGhM+nfL3PPyMxekyXx3wbMl/QDkiGdFwD/keXkEbEIaC3z/oQs5zErq7DRN3bNKpBlVM8XJN1LclMW4K8i4rF8wzKrwNYu9++bVaDcqJ4REdEFkCb63ZJ96TFmNVPs6jGzTMqVbPixpK9IerukluJOSeMlTZP0U5IRPma1Vej0UE6zCpTr6pkKfIBk0tWpkvYjKdT2PPAT4C8jYnn+IZr1odAJI/6g1lGYNYxywzkD+GH6MKtfhU7Y//BaR2HWMDItvWhW17Z2edlFswo48Vvj83BOs4o48VvjK3S5MqdZBTIlfkmnSPqz9PUYSS6XafVh+zbYvsWjeswq0OcELkl/D5wKHAXcArSQzOZ9R76hmWXgZRfNKpalxX8h8D6gEyAdwjkqz6DMMvOyi2YVy5L4t6RDO4tLL/r/MKsf3YuwuKvHLKssif9OSd8CRku6ArgfuCnfsMwy6k78bo+YZZWlSNtXJJ1Nsl7uCcCXIuLe3CMzy8LLLppVrGzilzQE+GlEvAdwsrf6U1x20RO4zDIr29UTEduBgiTfzLX6VFx20S1+s8yyLMSykWSN3PtJR/YARMSncovKLCt39ZhVLEvi/1n6MKs/hbSrx4nfLLMsN3dnDUQgZv3irh6zimWZufsc6Rj+UhFxTC4RmVViaxcgGNrS56FmlsjS1VNamqEFuAgYnU84ZhUqrr4l1ToSs4aRpavnlV12fU3SwznFY1aZwkZP3jKrUJaunsklm3sBrbjFb/Wi0OX+fbMKZenq+VbJ623Ai8AHs5xc0v7AjcDxJPcJPgxcQLKWbwF4AbgiItZVELPZawqdTvxmFcqS+C+NiN+V7qigHv8NwNyIuFDScGAEMA+4NiK2SfoKcC3wt5UEbdZta6dn7ZpVKEuRtrt62Hd3Xx9KZ/u+E5gFEBGFiFgXEfdHxLb0sF8D47IGa7Ybt/jNKtZri1/SMcDrSapynlvy1iiS0T19ORLoAG6WdAKwELgqIjpLjvkwMKfiqM2KCl0w6tBaR2HWUMq1+N9IsgjL/iRDOIuPtwMfyXDuocBJwHci4s0k5R7aim9Kuo7knsGtPX1Y0nRJCyQt6OjoyHA5a0rF4ZxmllmvLf6IuAu4S9I7IqI/wzfbgfaImJ9u30Ga+CVNA84BpqaLvPR0/ZnATIDW1tYejzGjsNGrb5lVKMvN3cckfYTkL4DuLp6ImF7uQxHxsqTfSzo2Ip4BpgJPSXovyc3cP4yIrj2I3SyZues+frOKZLm5ewswgaSFPp9k0fXNGc9/JXCrpCeAE4EvA98ERgLzJC2S9N1KgzYDYMd22LbZXT1mFcrS4j8mIj4o6f0RMUvSLcB9WU4eEYtIJnyVOrrSIM165GUXzfolS4t/a/q8TtLrSVrrR+QXkllGrsVv1i9ZWvyzJB0AfJakpT8C+MdcozLLopj4PYHLrCJZ1txdFRFrgZ8DWWfsmuVvq1v8Zv2RZc3dqwcoFrPKuKvHrF+y9PHfJ+lqSWMljSo+co/MrC9edtGsX7L08Rdn6X6apMKm0md3+1htedlFs37JshDL4QMRiFnFtqYtfs/cNatIn109kvaR1CbpO+n20ZLOzj80sz509/F7ApdZJbL08d+UHndaur2CZAauWW11d/W4xW9WiSyJf1JEfJl0IldaX8crW1vtFdL/FIfuU+tIzBpKlsRfkNRCckMXSRNJlk00q63iIix7ZfnP2MyKsozq+QIwFxgnaTbwh8Cf5xqVWRZbO31j16wfsozqmStpIckCLAKuiYiVuUdm1hcvu2jWL1la/ABvA04l6e7ZDvwot4jMsnLiN+uXLMM5vwFcBTwHPA98It1nVltO/Gb9kqXF/y7g+OISiZJuAp7INSqzLAqdsPfIWkdh1nCyDId4FhhXsj0WWJxPOGYV8LKLZv2SpcU/Gnha0q/T7ZOBX0m6EyAiLsgrOLOyChs9a9esH7Ik/i/lHoVZfxQ6PWvXrB+yDOd8AEDSiNLjI2JDjnGZ9a3grh6z/ugz8Uv6c+CLJMM4d+CyzFYPdmyHbZu87KJZP2Tp6mkDTvCkLasrW70Ii1l/ZRnVsxRwt47Vl+6SzO7jN6tU1hb/L9NRPVuKOyPiU7lFZdYX1+I367csif+7wC+BJ0n6+DOTtD9wI3A8yX2BDwPPAHOACcAy4OKIWFvJec280LpZ/2VJ/Dsi4hP9PP8NwNyIuFDScGAE8HfAAxExQ1IbyV8Uf9vP81uzKiZ+V+c0q1iWxP+ApA+TFGYr7eop2+8vaRTwTuDy9PgCSW3/84DT08NmAw/ixJ+PVc9D4dVaR5GPV9LJ4+7qMatYlsQ/LX3+fMm+LMM5jwQ6gJslnQAsJCn2dkhEvAQQES9JOriykC2Tl5+E776j1lHkb8SYWkdg1nCyTOA6fA/OfRJwZUTMl3QDSbdOJpKmA9MBxo/3lIGKrXo2eX7/P8Oow2obS15a9ocDj651FGYNJ8sErn1IWupHRMRHJR1Nsg7vvX18tB1oj4j56fYdJIn/FUlj09b+WKDH+QERMROYCdDa2hrZvo51W788eX7TRdAyuraxmFldyTKO/6b0uNPS7RXAl/v6UES8DPxe0rHprqnAU8APea37aBpwTyUBW0YblsPwkU76ZrabLH38kyLiEkkXAURElyRlPP+VwK3piJ6lwBUkv0RuT0tB/C9wUT/itr6sb4fRg7SLx8z2SJbEX5DUQnJDF0kTgUKWk0fEIqC1h7emZo7Q+mfD8sHbt29meyRLV88XgLnAOEmzgZ8D1+Yale259cvd4jezHvXa4pc0PiL+NyLmSloIvJ2kMuc1LthW57Ztgc6VMGpc38eaWdMp19VzN8lwTCKiA9+EbRwbViTPbvGbWQ/KdfVkvYFr9WZDOpTTffxm1oNyLf7DJH29tzf3oH6P5W29E7+Z9a5c4t9EUmbBGk2xxe+uHjPrQbnEvzoiZg9YJFY9G5Yn5QxcstjMelCujz/TWH2rQ+uXw2iP6DGznvWa+CPilIEMxKpoQ7v7982sV1kmcFmj8eQtMyvDiX+wKXTBpjVu8ZtZr/pM/JL+Pcs+qxPdk7fcx29mPcvS4n9j6YakIcBb8gnH9tiG9uTZLX4z60WviV/StZJeBSZL2pA+XiVZOMXlG+rVeo/hN7Pyyo3q+T8RMRK4PiJGpY+RETEmIlyds165XIOZ9SHLmrvXSjoMOKL0+Ih4KM/ArJ/Wt8O+B8HQvWsdiZnVqSxr7s4APkSybOL2dHcATvz1yAuwmFkfsqzAdT5wbERsyTsYq4L1y2HMUbWOwszqWJZRPUuBYXkHYlXiFr+Z9SFLi78LWCTpAaC71e+yzHVo8wbYssEjesysrCyJ/4fpw+qdR/SYWQZZRvXMlrQPMD4inhmAmKy/usfwe9aumfUuS8mGDwCLgLnp9omS/BdAPfKsXTPLIMvN3c8BU4B1ABGxCJiYY0zWX+uXg/aCkWNrHYmZ1bEsffzbImK9tNPa65Hl5JKWAa+SjP/fFhGtkk4Evgu0ANuAj0XEoxVFbT3bsBz2ex0MyfLPambNKkuGWCzpT4AhkiYBnwB+VcE1zoiIVSXbXwU+HxH3Snpfun16Beez3qxv94geM+tTlsR/JXAdyVDO/wTuA/5pD64ZwKj09WhgxR6cqzIbVsDqFwbscgNuzYtw2Em1jsLM6lyWUT1dJIn/un6cP4D7JQXwbxExE7gauE/S10juMby9pw9Kmg5MBxg/fnw/Lt2D/7gQVv62OueqV5MvrnUEZlbnek38kv41Iq6W9CN66NOPiHMznP/UiFgh6WBgnqQlwIXAJyPiB5IuBmYB7+7h/DOBmQCtra2Z7in0aUM7vP5cmDK9KqerO9rLLX4z61O5Fn9xla2v9ffkEbEifV4p6S6S0UHTgKvSQ74P3Njf81dk+zbYvB4OeSNMPG1ALmlmVo/K1eNfmD7/AngEWAusAR5J95UlaV9JI4uvgbOAxSR9+n+YHvYu4Lk9+QKZbVqbPI8YMyCXMzOrV1nKMr+fZPjlC4CAiZI+EhH39vHRQ4C70mGgQ4HbImKupI3ADZKGAptJ+/Fzt2lN8rzPAQNyOTOzepVlVM8/kwzJfB5A0lHAT4CyiT8ilgIn9LD/YWqxZm/X6uR5xB8M+KXNzOpJlpm7K4tJP7WUZN3dxtKVtvjd1WNmTa7cqJ4L0pe/lfRT4HaS0T0XAY8NQGzV1d3V4xa/mTW3cl09Hyh5/Qqv3ZDtABqvo7y7xe/Eb2bNrdfEHxFXDGQguetaDUP2hmEjah2JmVlNZRnVM5GkbMOE0uMzTuCqH5vWJP37OxebMzNrOllG9dxNMrv2R8COfMPJUddad/OYmZEt8W+OiK/nHkneulZ7DL+ZGdkS/w2SPgvcz86LrT+eW1R52LQGDn5DraMwM6u5LIn/TcBlJOUVil09kW43jq417uoxMyNb4j8fODIiCnkHk5sdO5IWv8fwm5llmrn7G2D/vAPJ1Zb1EDs8a9fMjGwt/kOAJZIeY+c+/sYZzunJW2Zm3bIk/s/mHkXeulyuwcysKMvSizvV3pd0KvAnQJ81+evGJrf4zcyKsrT4kXQiSbK/GHgR+EGeQVWdu3rMzLqVq855DPAh4BJgNTAHUEScMUCxVU+xFr+7eszMyrb4lwD/DXygZBGWTw5IVNW2aQ1oCLSMrnUkZmY1V2445x8DLwM/l/T/JE0lWXqx8RQnb7lAm5lZ2cXW74qIDwLHAQ8CnwQOkfQdSWcNUHzV0bXa3TxmZqk+J3BFRGdE3BoR5wDjgEVAW+6RVdMmV+Y0MyvKMnO3W0SsiYh/i4gGrNPjWbtmZlBh4m9Ym9a4JLOZWWrwJ/6IpI/fXT1mZkAzJP5CJ2wvuKvHzCyVa+KXtEzSk5IWSVpQsv9KSc9I+q2kr+YZQ3e5Bo/qMTMDMpZs2ENnRMSq4oakM4DzgMkRsUXSwblevThr1109ZmZAbbp6PgrMiIgtABGxMteruTKnmdlO8k78AdwvaaGk6em+Y4DTJM2X9AtJb+3pg5KmS1ogaUFHR0f/I9i0Nnl2H7+ZGZB/V8+pEbEi7c6ZJ2lJes0DgFOAtwK3SzoyIqL0gxExE5gJ0NraGvSXu3rMzHaSa4s/IlakzyuBu4ApQDtwZyQeJVnA/cDcgih29bQ09uqRZmbVklvil7SvpJHF18BZwGLgbuBd6f5jgOHAqt7Os8c2rUmS/pCBuI9tZlb/8syGhwB3KamIORS4LSLmShoO3CRpMVAApu3azVNVnrxlZraT3BJ/RCwFTuhhfwG4NK/r7qZrjUf0mJmVGPwzdze5QJuZWanBn/iLi7CYmRnQLInfXT1mZt0Gd+Lfuhm2dsIIl2Q2Mysa3Im/WKDNffxmZt0Gd+J3nR4zs90M7sTf3eJ34jczKxrcib/LXT1mZrsa5Ik/LdDmrh4zs26DO/G7q8fMbDeDO/F3rYXh+8HQvWsdiZlZ3Rjcif+gY+GN59c6CjOzujK4axW/ZVryMDOzboO7xW9mZrtx4jczazJO/GZmTcaJ38ysyTjxm5k1GSd+M7Mm48RvZtZknPjNzJqMIqLWMfRJUgfwu35+/EBgVRXDaQT+zs3B37k57Ml3PiIiDtp1Z0Mk/j0haUFEtNY6joHk79wc/J2bQx7f2V09ZmZNxonfzKzJNEPin1nrAGrA37k5+Ds3h6p/50Hfx29mZjtrhha/mZmVcOI3M2syTZH4Jf2TpCckLZJ0v6RDax1T3iRdL2lJ+r3vkrR/rWPKm6SLJP1W0g5Jg3bIn6T3SnpG0vOS2modz0CQdJOklZIW1zqWgSDpcEk/l/R0+t/0VdU8f1MkfuD6iJgcEScCPwb+sdYBDYB5wPERMRl4Fri2xvEMhMXABcBDtQ4kL5KGAN8CzgbeAFwi6Q21jWpAfA94b62DGEDbgE9HxOuBU4CPV/PfuSkSf0RsKNncFxj0d7Qj4v6I2JZu/hoYV8t4BkJEPB0Rz9Q6jpxNAZ6PiKURUQD+CzivxjHlLiIeAtbUOo6BEhEvRcTj6etXgaeBw6p1/sG95m4JSV8C/gxYD5xR43AG2oeBObUOwqriMOD3JdvtwMk1isUGgKQJwJuB+dU656BJ/JJ+Bryuh7eui4h7IuI64DpJ1wJ/DXx2QAPMQV/fOT3mOpI/G28dyNjykuU7D3LqYd+g/wu2WUnaD/gBcPUuPRd7ZNAk/oh4d8ZDbwN+wiBI/H19Z0nTgHOAqTFIJmxU8O88WLUDh5dsjwNW1CgWy5GkYSRJ/9aIuLOa526KPn5Jk0o2zwWW1CqWgSLpvcDfAudGRFet47GqeQyYJGmipOHAh4Af1jgmqzJJAmYBT0fE/636+QdJQ7AsST8AjgV2kJR3/quIWF7bqPIl6Xlgb2B1uuvXEfFXNQwpd5LOB74BHASsAxZFxHtqG1X1SXof8K/AEOCmiPhSjUPKnaT/BE4nKVH8CvDZiJhV06ByJOkdwH8DT5LkLYC/i4ifVuX8zZD4zczsNU3R1WNmZq9x4jczazJO/GZmTcaJ38ysyTjxm5k1GSd+qzpJY9JKqIskvSxpecn2r3K65psl3Zi+Pk7SI5K2SPrMLsfVdZVHSdMkPZc+pvVyzIMDWX1U0l9LumKgrmf583BOy5WkzwEbI+JrOV/n+8AXI+I3kg4GjgD+CFhbem1J7wQ2ArdExPF5xlQpSX8ALABaScowLATeEhFrdznuQeAzEbGgytcfEhHbe9g/AvhlRLy5mtez2nGL3waUpI3p8+mSfiHpdknPSpoh6U8lPSrpSUlHpccdJOkHkh5LH6f2cM6RwOSI+A1ARKyMiMeArbsem6XKo6QPSJov6X8k/UzSIen+z0mana7psEzSBZK+msY7N51ivyfeA8yLiDVpsp9H76WIL0p/Vs9KOi2Nr0XSzWk8/yPpjHT/5ZK+WfL9fizp9PT1RklfkDQfeFv67/BUuo7D1wDSmd/LJE3Zw+9ndcKJ32rpBOAq4E3AZcAxETEFuBG4Mj3mBuBfIuKtwB+n7+2qlaQWf7U8DJyStnD/C/ibkveOAt5PUgr5P4CfR8SbgE3p/p1Iuqakm6v08fUerttT5c3eSvEOTX9WV/Na3amPA6TxXALMltTSx3fdF1gcEScDTwHnAyWxHVgAAAIuSURBVG9M13H4YslxC4DT+jiXNYhBU6TNGtJjEfESgKQXgPvT/U/yWunsdwNvSEqXADBK0si0RnnRWKCjinGNA+ZIGgsMB14see/eiNgq6UmSkglzS2KesOuJIuJ64PqM162k8maxaNfCkuu+g6RkBRGxRNLvgGP6uOZ2kkJgABuAzcCNkn5CsmhR0UrguD7OZQ3CLX6rpS0lr3eUbO/gtUbJXsDbIuLE9HHYLkkfktZ2Xy3bSnwD+Gbacv7ILufeAhARO4CtJVVPS2PuVmGLv5LKm8Wf1faS6/b0iwOSstyl/6+Xfp/NxX79dOGeKSS/CP6I136pFT+zqZfzW4Nx4rd6dz/J+gkASDqxh2OeBo6u4jVHA8Uifj2OrMkqIq4v+aVV+vhED4ffB5wl6QBJBwBnpfuyegj4UwBJxwDjgWeAZcCJkvaSdDhJct+Nktrvo9NCYFcDpT/rY6hud5rVkBO/1btPAK3pzcangN0qjEbEEmB0epMXSa+T1A58Cvh7Se2SRqXv/SfwCHBsuv/Pe7jm54DvS/pvYFUu36oHEbEG+CeS0suPAV9I92X1bWBI2g01B7g8IrYAvyTprnoS+BrweC+fHwn8WNITwC+AT5a8dyrwswpisTrm4Zw2KEj6JPBqRPR089f2gKQ3A5+KiMtqHYtVh1v8Nlh8h53vGVj1HAj8Q62DsOpxi9/MrMm4xW9m1mSc+M3MmowTv5lZk3HiNzNrMk78ZmZN5v8DXCSCFkh8uYMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Part A.\n",
    "def ambient_T(t, spec):\n",
    "    \"\"\"Accepts times in hours where 11 am = 0 hrs and returns the ambient temperature.\n",
    "    \n",
    "    spec is a parameter that specifies a time before 11 am or after,\n",
    "    given as 'b' for before\n",
    "    and 'a' for after\n",
    "    \"\"\"\n",
    "    if spec == 'b':\n",
    "        B = {-3 : 55, -2 : 58, -1 : 60, 0 : 65}\n",
    "        return B.get(int(t))\n",
    "    else:\n",
    "        A = {0 : 65, 1 : 66, 2 : 67}\n",
    "        return A.get(int(t))\n",
    "    \n",
    "t_A = np.linspace(0,2)    \n",
    "t_B = np.linspace(-3,0)\n",
    "ambient_A = [ambient_T(i,'a') for i in t_A]\n",
    "ambient_B = [ambient_T(i,'b') for i in t_B]\n",
    "\n",
    "plt.xlabel('Time (11 am = 0 hours)')\n",
    "plt.ylabel('Ambient Temperature (degF)')\n",
    "plt.plot(t_A, ambient_A, label = 't < 0')\n",
    "plt.plot(t_B, ambient_B, label = 't > 0')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Part a continued.`The changes in temperatures between the hours are not this abrupt so this graph does not accurately reflect how the ambient temperature is changing. Perhaps having more data points between each hour rather than just at the start and end of the hour would produce a more accurate temperature distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Part B:\n",
    "N = 100\n",
    "t_A = np.linspace(0,2,N)\n",
    "dt_A = np.diff(t_A)\n",
    "t_B = np.linspace(0,-3,N)\n",
    "dt_B = np.diff(t_B)\n",
    "ambient_A = [ambient_T(i,'a') for i in t_A]\n",
    "ambient_B = [ambient_T(i,'b') for i in t_B]\n",
    "\n",
    "numerical_A = np.empty(len(t_A))\n",
    "numerical_A[0] = 85\n",
    "numerical_B = np.empty(len(t_B))\n",
    "numerical_B[0] = 85\n",
    "\n",
    "for i in range(1, len(t_A)):\n",
    "    numerical_A[i] = numerical_A[i-1] - K*(numerical_A[i-1] - ambient_A[int(i-1)])*dt_A[i-1]\n",
    "        \n",
    "\n",
    "for i in range(1, len(t_B)):\n",
    "    numerical_B[i] = numerical_B[i-1] - K*(numerical_B[i-1] - ambient_B[-int(i-1)])*dt_B[i-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "times1 = np.linspace(0,2,10)\n",
    "times2 = np.linspace(0,-3,20)\n",
    "analytical1 = np.empty(len(times1))\n",
    "analytical2 = np.empty(len(times2))\n",
    "\n",
    "for i in range(len(times1)):\n",
    "    analytical1[i] = T_a + (T_0 - T_a)*np.exp(-K*times1[i])\n",
    "    \n",
    "for i in range(len(times2)):\n",
    "    analytical2[i] = T_a + (T_0 - T_a)*np.exp(-K*times2[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUZfbA8e8JoYTewgJCKCuIdCQUKYJKiYhBUBAWVFggysICIosg68JPkdVdFBRUlia6VKUrLlWKgAqhKAgquksTpCsgkZbz+2Nu4gCTZJJMS3I+zzPPzL1z575nIubk3vd9zyuqijHGGAMQFuwAjDHGhA5LCsYYY5JZUjDGGJPMkoIxxphklhSMMcYkCw92AJlRsmRJrVixYrDDMMaYLGX79u2nVDXS03tZOilUrFiR+Pj4YIdhjDFZiogcTOk9u31kjDEmmSUFY4wxySwpGGOMSZal+xSMMdnDlStXOHLkCL/++muwQ8lW8uXLR7ly5cidO7fXn7GkYIwJuiNHjlCoUCEqVqyIiAQ7nGxBVTl9+jRHjhyhUqVKXn/Obh8ZY4Lu119/pUSJEpYQfEhEKFGiRLqvviwpGGNCgiUE38vIzzRHJoULFy4wcOBAfvrpp2CHYowxISVHJoUvv/ySyZMn06lTJy5fvhzscIwxIUBEePrpp5O3x40bx+jRowMaQ3x8PAMHDszQZ1u2bOmTybw5Mik0adKEGTNmsG7dOnr37o0tNGSMyZs3L4sWLeLUqVNBaf/q1atER0fz+uuvB6X9JDkyKQD06NGDMWPGMGvWLJ577rlgh2OMCbLw8HDi4uIYP378Te/17NmTBQsWJG8XLFgQgPXr19OiRQu6dOlC1apVGT58OLNnz6Zhw4bUqlWL77//HoCTJ0/y0EMP0aBBAxo0aMDmzZsBGD16NHFxcbRp04bHHnuM9evX0759e8B1m7tXr17UqlWL2rVrs3DhQgD69etHdHQ0NWrUYNSoUb7/Ofj8jFnIs88+y6FDh3jxxRcpV64cTz75ZLBDMibHGzx4MLt27fLpOevWrcuECRPSPK5///7Url2bYcOGeX3uL774gn379lG8eHEqV65Mnz592Lp1K6+99hoTJ05kwoQJDBo0iKeeeopmzZpx6NAh2rZty759+wDYvn07mzZtIiIigvXr1yef94UXXqBIkSLs3r0bgLNnzwLw4osvUrx4ca5du8a9997Ll19+Se3atdPx00id35KCiMwA2gMnVLWms68zMBq4HWioqvHO/orAPuAb5+Ofqarff0OLCG+88QY//PAD/fv3p2zZssTGxvq7WWNMiCpcuDCPPfYYr7/+OhEREV59pkGDBpQpUwaA3//+97Rp0waAWrVqsW7dOgDWrFnD3r17kz9z7tw5zp8/D0BsbKzHttasWcO8efOSt4sVKwbAe++9x5QpU7h69SrHjh1j7969WSMpADOBScC7bvv2AJ2Af3k4/ntVrevHeDwKDw9n/vz53H333XTt2pW1a9dy5513BjoMY4zDm7/o/Wnw4MHccccd9OrVK3lfeHg4iYmJgGtSmPsAlbx58ya/DgsLS94OCwvj6tWrACQmJvLpp596/OVfoEABj3Go6k1DSv/3v/8xbtw4tm3bRrFixejZs6fPZ4H7rU9BVTcCZ27Yt09Vv0nhI0FToEABli9fzi233EL79u35+uuvgx2SMSZIihcvTpcuXZg+fXryvooVK7J9+3YAli5dypUrV9J1zjZt2jBp0qTkbW9uj934mbNnz3Lu3DkKFChAkSJFOH78OP/5z3/SFYc3QqmjuZKI7BSRDSLSPKWDRCROROJFJP7kyZM+azwyMpKVK1cSHh5OTEwMR48e9dm5jTFZy9NPP33dKKS+ffuyYcMGGjZsyOeff57iX/cpef3114mPj6d27dpUr16dyZMnp/mZv/71r5w9e5aaNWtSp04d1q1bR506dahXrx41atTgj3/8I02bNk33d0uL+HM4ptNX8GFSn4Lb/vXAULc+hbxAQVU9LSL1gSVADVU9l9r5o6Oj1deL7Gzfvp2WLVtSqVIlNm7cSNGiRX16fmPMzfbt28ftt98e7DCyJU8/WxHZrqrRno4PiSsFVb2kqqed19uB74GqwYilfv36LFq0iK+//poHH3zQqjYaY3KUkEgKIhIpIrmc15WBKsB/gxVP69ateeedd9iwYQN/+MMfuHbtWrBCMcaYgPJbUhCRucCnwG0ickREeotIRxE5AtwJLBeRlc7hdwFfisgXwALgSVU94/nMgdGtWzdee+01Fi9eTL9+/WzWszEmR/DbkFRV7ZbCW4s9HLsQWOivWDJq4MCBHD9+nLFjxxIZGcmLL74Y7JCMMcavcvSMZm+MGTOGkydPMnbsWEqWLMlTTz0V7JCMMcZvLCmkQUR46623OHPmDEOGDKF48eI8/vjjwQ7LGGP8IiQ6mkNdrly5mD17Nq1ataJ3794sW7Ys2CEZY/xg8eLFiEimJrDeWDzPk7Fjx1633aRJkwy1NXr0aMaNG5ehz6bEkoKX8ubNy+LFi4mOjqZLly7JNU2MMcEwG6iI61dYRWc78+bOnUuzZs2uqznkDzcmhS1btvi1vfSwpJAOBQsWZPny5dx6663ExsaydevWYIdkTA40G4gDDgLqPMeR2cRw4cIFNm/ezPTp05OTwvr162nZsiUPP/ww1apVo3v37skjEZ9//nkaNGhAzZo1iYuLu2mE4tq1a+nYsWPy9urVq+nUqRPDhw8nISGBunXr0r17d+C3UtwA//jHP6hVqxZ16tRh+PDhAEydOpUGDRpQp04dHnroIS5evJip75oaSwrpVKJECVatWkWpUqWIiYlJLmtrjAmUkcCNvxQvOvszbsmSJcTExFC1alWKFy/Ojh07ANi5cycTJkxg7969/Pe//01eC2HAgAFs27aNPXv2kJCQwIcffnjd+e655x727dtHUjmet99+m169evHSSy8RERHBrl27mD37+kT2n//8hyVLlvD555/zxRdfJJfw7tSpE9u2beOLL77g9ttvv64uk69ZUsiAsmXLsnr1aiIiImjdujX79+8PdkjG5CCH0rnfO3PnzqVr164AdO3alblz5wLQsGFDypUrR1hYGHXr1uXAgQMArFu3jkaNGlGrVi0+/vhjvvrqq+vOJyI8+uijzJo1i59++olPP/2U++67L9UY1qxZQ69evcifPz/gKs4HsGfPHpo3b06tWrWYPXv2TW35ko0+yqDKlSuzevVqWrRoQatWrfjkk0+IiooKdljG5ABRuG4ZedqfMadPn+bjjz9mz549iAjXrl1DRGjXrt11pbFz5crF1atX+fXXX/nTn/5EfHw85cuXZ/To0R5L4vTq1YsHHniAfPny0blzZ8LDU/+V66lcNrg6r5csWUKdOnWYOXPmdYvx+JpdKWRC9erVWbVqFT///DOtWrXixx9/DHZIxuQALwL5b9iX39mfMQsWLOCxxx7j4MGDHDhwgMOHD1OpUiU2bdrk8fikBFCyZEkuXLiQ4mijsmXLUrZsWcaMGUPPnj2T9+fOndtj+e02bdowY8aM5D6DM2dchR3Onz9PmTJluHLlyk23nHzNkkIm1atXj48++oijR4/SunVrTp8+HeyQjMnmugNTgAqAOM9TnP0ZM3fu3Os6hQEeeugh5syZ4/H4okWL0rdvX2rVqsWDDz5IgwYNUo62e3fKly9P9erVk/fFxcVRu3bt5I7mJDExMcTGxhIdHU3dunWTh5u+8MILNGrUiNatW1OtWrWMfk2v+LV0tr/5o3R2Rn388ce0a9eOGjVqsHbtWiu5bUw6ZOfS2QMGDKBevXr07t07KO1nydLZ2cE999zDokWL2L17N+3atUtef9UYk3PVr1+fL7/8kh49egQ7FK9ZUvChdu3aMW/ePLZu3Ur79u39OpbYGBP6tm/fzsaNG6/rrA51OTQp+Gc2JLjGE8+aNYtNmzbRoUMHEhISfHZuY7KzrHwrO1Rl5GeaA5OCf2ZDuuvatSszZsxg7dq1dOrUiUuXLvns3MZkR/ny5eP06dOWGHxIVTl9+jT58uVL1+dy4DyF1GZDZnz0wo0ef/xxrly5Qt++fXn44YdZuHAhefLk8dn5jclOypUrx5EjR5Jn/xrfyJcvH+XKlUvXZ/yWFERkBtAeOKGqNZ19nYHRwO1AQ1WNdzt+BNAbuAYMVNWVN53UJ/wzG9KTPn36cPXqVfr160eXLl147733LDEY40Hu3LmpVKlSsMMw+Pf20Uwg5oZ9e4BOwEb3nSJSHegK1HA+82bSms2+l9KsR//MRn7yySeZOHEiS5cupWvXrh4nrBhjTKjwW1JQ1Y3AmRv27VPVbzwc3gGYp6qXVPV/wHdAQ/9E5vvZkGkZMGBA8nrPlhiMMaEsVDqabwEOu20fcfbdRETiRCReROIzdv/R97MhvTFw4EAmTJjAokWL6NatmyUGY0xICpWkcHMFKNfQoJt3qk5R1WhVjY6MjMxgc92BA0Ci8+zfhJBk0KBBjB8/noULF9oVgzEmJIVKUjgClHfbLgccDVIsfjV48GDGjx/PokWLeOSRR7h8+XKwQzLGmGShkhSWAV1FJK+IVAKqANl2WbPBgwcn9zF07tzZ5jEYY0KG35KCiMwFPgVuE5EjItJbRDqKyBHgTmC5iKwEUNWvgPeAvcAKoL+qXvNXbN7x36xncPUxvPHGGyxbtoxOnTp5rMVujDGBZlVSPUqa9ew+yS0//uiQnjJlCk888QRt27Zl8eLFRERE+PT8xhhzI6uSmm7+WQPWk7i4OGbMmMGqVato3749v/zyi8/bMMYYb1lS8Chws57BtWTfu+++y/r167nvvvs4d+6cX9oxxpi0WFLwKLCzngF69OjBvHnz+PTTT2ndujVnz571W1vGGJMSSwoeBX7WM0Dnzp1ZsGABu3bt4p577rHiYMaYgLOk4FFwZj0DdOjQgWXLlvH111/TokULjh7NltM1jDEhypJCioIz6xmgbdu2rFixgsOHD9O8eXMOHDgQsLaNMTmbJYUQ1aJFC9asWcOZM2do3rw533zjqY6gMcb4liWFTPPfJLdGjRqxYcMGLl++TPPmzdm5c6fPzm2MMZ5YUsgU/y/tWbt2bT755BPy5ctHy5Yt2bRpk8/ObYwxN7KkkCmBmeRWtWpVNm3aROnSpWnTpg0rVqzw6fmNMSaJJYVMCdwkt6ioKD755BOqVatGbGws8+fP93kbxhhjSSFTAjvJrVSpUqxbt47GjRvTrVs3Jk+e7Jd2jDE5lyWFTAn8JLciRYqwcuVK7r//fvr168eYMWPIykUNjTGhJTytA0SkLtAcKAskAHuAtar6s59jywKS5i6MxHXLKApXQvDvnIaIiAgWLVpE7969ee655zhx4gQTJkwgLMxyvDEmc1JMCiLSA3gK16po23ENrckHtAKeE5EdwChVPRKIQENXdwI5sS1J7ty5mTlzJpGRkbz66qucOnWKmTNnkidPnoDHYozJPlK7UigB3KWqHms5i0g0cDuupGGCICwsjHHjxlGqVCmGDx/OqVOnWLhwIYUKFQp2aMaYLCrF+w2q+lpKCcF5P15VV6f0vojMEJETIrLHbV9xEVktIvud52LO/ooikiAiu5xHNutB9d8ENxHhmWee4e233+bjjz/m7rvv5vjx4z47vzEmZ0kxKYjIf9xeD8vAuWcCMTfsG46rP6IKsNbZTvK9qtZ1Hk9moL0Q5f8JbgA9e/Zk6dKl7Nu3j6ZNm/Ldd9/59PzGmJwhtZ7J0m6vu6b3xKq6EThzw+4OwDvO63eAB9N73qwncKu43X///Xz88cf89NNPNGnShK1bt/q8DWNM9pZaUvDHOMffqeoxAOe5lNt7lURkp4hsEJHmKZ1AROJEJF5E4rPGegOBXcWtUaNGbNmyhYIFC3L33XezcuVKv7RjjMmeUksKlUVkkYgsdnud/PBxHMeAKFWtBwwB5ohIYU8HquoUVY1W1ejIyEgfh+EPgV/FrWrVqmzZsoUqVarQoUMHK4thjPFaaqOPHnJ7PclH7R0XkTKqekxEygAnAFT1EnDJeb1dRL4HqgLxPmo3iF7E1YfgfgvJ/6u4lS5dmrVr19K6dWs6dOjA4sWLadeunV/bNMZkfamNPlqb2iOD7S0DHndePw4sBRCRSBHJ5byuDFQB/pvBNkJM8FZxK1GiBGvWrKFmzZp07NiRDz74wO9tGmOyNkmrRIKI7OTm/oWfcf0V/3dVvbEzOelzc4GWQEngODAKWAK8h+veySGgs6qeEZGHgOeBq8A1XJPi0vwNFh0drfHx2eBiws/Onj1L27Zt2bVrF/Pnz6djx47BDskYE0Qisl1Voz2+50VS+AeuP3HnOLu64vrFfQForKqxPow1XbJfUpiNv0pm/Pzzz8TExLBt2zbmzJlDly5dfHJeY0zWk1pSSLP2EdBEVZu5be8UkU2q2kxEdvsmRPPbfIakvoek+Qzgi8TgXkivW7duXL58mR49emT6vMaY7MWbCmqFRKR+0oaI3AEkjQy66peociT/z2coXLgwK1asoGXLljz22GPMmDHDZ+c2xmQP3lwpPAH8W0Ry47qNdBnoLSIFgH/4M7icJTDzGQoUKMCHH35Ip06d6N27NwkJCfTv39+nbRhjsq40k4KqfgZUF5ESuPogTrm9PddvkeU4UbhuGXna71sREREsWbKERx55hAEDBnDx4kX+8pe/+LwdY0zWk+btI2e46L+Ad1T1lIhUF5Ge/g8tpwnsgj158+bl/fffp0uXLgwbNoxRo0bZYj3GGK/6FGYCG4DyzvZ+4Gl/BZRzBX4+Q+7cuZkzZw69evXi+eefZ+jQoZYYjMnhvOlTKKWqc0TkLwCqekVErvk5rhwq8Av25MqVi2nTplGwYEFeffVVzp8/z1tvvUWuXLkCGocxJjR4c6Xwi4gUx5nAJiINgPN+jcqkwbfrM4SFhfHaa68xcuRIpk6dSo8ePbhy5UqmozTGZD3eXCkMBT7AVRRvA3AL8LBfozKp8M98BhFhzJgxFC5cmGeeeYbz58/z/vvvExERkblwjTFZSpozmgFEJA+upTcF2Kuql/0dmDey34xmb1TE8yilCsABn7Twr3/9i379+tG8eXOWLVtGkSJFfHJeY0xoyNCMZhFJqXxFlIigqst8Ep1JJ//PZ3jiiScoWrQoPXr04O6772bFihWUKlUq7Q8aY7K81G4fdXaeSwJNgHW4rhRa4BqNZEkhKAIzn+GRRx6hSJEidOrUiWbNmrF69WoqVKjg0zaMMaEntdLZj6rqo8AVoLqqPqiqHYAaWHmLIArcfIaYmBjWrFnDyZMnadKkCV999ZXP2zDGhBZvRh9VVtUf3LaPArf5KR6TpsDOZ2jSpAkbN25EVWnevDlbtmzxSzvGmNDgTVLYKCLLRaSHiHTHddtoo5/jMqnqjqtTOdF59u/chlq1arF582ZKlChBq1atWL58uV/bM8YEjzdJoT+uWc2NgMbAu84+E9J8O5ehUqVKbN68merVq9OhQwdmzpyZ2QCNMSEozaSgLu+r6p+dx/vqxThWEZkhIidEZI/bvuIislpE9jvPxdzeGyEi34nINyLSNuNfyfw2l+EgrjmHSXMZMpcYSpUqxbp167jnnnvo1asXL730kpXFMCabSTEpiMg6EeknImVv2B8uIneJyHQR6ZXKuWcCMTfsGw6sVdUqwFpnGxGpjmtFtxrOZ95MWrPZZIT/1mYoVKgQH374Id26dWPEiBEMHDiQa9es6okx2UVqQ1LvB/oAi0XkFuAMEAHkxfUL/Q1VTXHmmKpuFJGKN+zugGvdZoB3gPXAM87+eap6CfifiHwHNAQ+Td/XMS7+ncuQJ08eZs2aRdmyZXnllVc4evQos2bNstnPxmQDKSYFVb0IvA68LiJ5gVJAwg3rKaTX71T1mHP+YyKSNCPqFuAzt+OOOPtuIiJxOHUdoqJ8v9ZA9uD/uQxhYWGMGzeOcuXKMWTIENq0acPSpUspXry4z9owxgSeNx3NqOolVT2cyYSQGvHUbAqxTFHVaFWNjoyM9FM4WV3g5jIMHjyYefPmsXXrVpo2bcqBAwd83oYxJnC8Sgo+dFxEygA4zyec/Uf4bb0GgHK45kOYDAnsXIYuXbqwatUqfvzxRxo3bsz27dv90o4xxv8CnRSWAY87rx8Hlrrt7yoieUWkElAF2Brg2LIZb+cy+GboaosWLdi8eTP58uXjrrvusrkMxmRRXiUFESknInc7r/OKSAEvPjMXV0fxbSJyRER6Ay8BrUVkP9Da2UZVvwLeA/YCK4D+qmpDWvzOt0NXq1evzmeffUa1atWIjY3lrbfe8lmkxpjASLN0toj8ERgAFFHV34tIVeBNVW0ViABTkzNLZ/tSRfxRhvvChQt07dqV5cuXM3ToUF5++WXCwgJ9UWqMSUlqpbO9+T91IK6ZzOcAVPVbXCORTJbnn6GrBQsWZMmSJfTv359x48bRuXNnLl68cd6EMSYUeZMUfnVfVMeZVOZptJDJclIaopr5oavh4eFMnDiR8ePHs3jxYlq2bMmPP/6Y6fMaY/zLm6SwWUSGAfmcfoX5wIf+DcsEhn+HrooIgwcPZsmSJezdu5dGjRrx5Zdf+uTcxhj/8CYpDAPOA18Dg3DNZs58vQQTAgIzdDU2NpaNGzdy9epVmjZtaiOTjAlhqSYF51bRDFV9S1U7OgvtvKWqiQGKz/idN0NXMz9s9Y477mDr1q1UrVqV2NhYxo8fb8X0jAlBqSYFZ1hoGRHJHaB4TMjx3bDVW265hY0bN9KhQweGDBnCE088weXLl9P+oDEmYLy5ffRf4BOntPXApIe/AzOhwrcVVwsUKMCCBQsYMWIEU6dOpW3btpw+fTqzQRpjfMSbpHASWI2rBzLS7WFyBN8PWw0LC2Ps2LH8+9//5tNPP6Vhw4bs3bs3w+czxvhOmpPXQplNXguEivhjgluSzz77jAcffJCLFy8yZ84c2rdvn+lzGmNSl6nJa84KaatufPg+TBOa/DtstXHjxmzbto0qVaoQGxvLyy+/bB3QxgRRaovsJPmr2+t8wEPAJf+EY0JP0mikkbhuGUXhSgi+G7Zavnx5PvnkE/74xz8yfPhwvvjiC6ZNm0b+/DcmI2OMv3mzRvPnbo8NqjoQ16poJsfwf8XV/PnzM3fuXMaOHcu8efNo1qwZhw75ZqU4Y4z3vLl9VNjtUVRE7gXKBCA2k6VkfuiqiDBixAg++OADvv/+e6Kjo9mwYYNfojXGeObN6KOvgD3O805c9xH6+jMokxX5bujq/fffz9atWylevDitWrVi4sSJ1s9gTIB4kxQqq2qUqpZX1Uqqeg+w2d+BmazGt0NXb7vtNj7//HPuu+8+Bg4cSM+ePUlISMh4eMYYr3iTFD73sC9Tq6KJyCAR2SMiX4nIYGffaBH5QUR2OY92mWnDBJrvK64WKVKEJUuWMHr0aN59912aN2/OwYOehscaY3wlxaQgIqVEpA4QISK1RKS282jGzWMUvSYiNXHdfmoI1AHai0gV5+3xqlrXeXyU0TZMMPhn6GpYWBijRo1i2bJl7N+/n/r167N69epMndMYk7LUrhTuByYB5YA3gTecx7PAc5lo83bgM1W9qKpXgQ1Ax0ycz4QE/1ZcfeCBB4iPj6dMmTLExMQwduxYEhOtLqMxvubNcpxdVPU9nzUocjuwFLgTSMBVijseOA30xLXCWzzwtKqe9fD5OFzDWoiKiqpvtxOyotlkdN7DL7/8Qp8+fZg3bx6xsbG88847FC1a1H+hGpMNpTaj2asyFyLSFqiBa/IaAKo6NhMB9Qb6AxeAvbiSw0vAKVzjGV8AyqjqH1M7j5W5yIqShq66j1TKT3quKlSV119/naFDh1KhQgUWLVpE7dq1fR+qMdlUZstcvAk8DgwBIoAewK2ZCUhVp6vqHap6F3AG2K+qx1X1mrNWw1Rsglw2lfmhqyLCoEGDWL9+PQkJCTRq1IiZM2f6MEZjci5vRh81U9U/AKdV9TmgEa5+hgwTkVLOcxTQCZgrIu4T4jrimhthsh3fDV1t2rQpO3bsoEmTJvTq1Ys+ffrYsFVjMsmbpPBr0rOIlHa2K2ay3YUishf4AOjv9B38Q0R2i8iXwN3AU5lsw4Qk3w5d/d3vfseqVasYOXIk06dPp3Hjxnz//fcZD8+YHM6bpPCRiBQFxgG7cBW/WZCZRlW1uapWV9U6qrrW2feoqtZS1dqqGquqxzLThglVvh+6mitXLsaMGcPy5cs5fPgw0dHRrFixIjNBGpNjpbVGcxjwH1X9SVXfByoBtVT12YBEZ7Ihb4eupr+4Xrt27YiPjycqKop27drx97//3cpjGJNOaa3RnAi85radoKpn/B6VyebSqrqa8eJ6lStXZsuWLXTt2pVnn32Whx9+mHPnzvkudGOyOW9uH60WkQ5+j8SYZJkboVSgQAFmz57Nq6++ytKlS2nUqBH79u3zdZDGZEveJIUBwGIRSRCRMyJyVkTsasH4UeZHKIkITz31FGvWrOH06dM0bNiQhQsX+iY8Y7Ixb5JCSSA3UBCIdLYj/RmUyel8N0KpZcuW7Nixg5o1a/Lwww8zdOhQrl69mrnwjMnGvFl57RrQGXjGeV0GqOvvwExO5tsRSuXKlWPDhg3079+fV155hXvvvZdjx2xwmzGeeDOjeRKueQOPOrsuApP9GZTJ6dJTXM+7UUp58uRh0qRJzJo1i23btnHHHXfYqm7GeODN7aMmqvoEziQ2Z/RRHr9GZYxX60Knf5RS9+7d2bp1K4ULF+bee+/l5ZdftmqrxrjxJilcceYrKICIlMD1f6oxQZaxUUo1a9Zk27ZtdOrUieHDh9OhQwfOnLGxE8aAd0nhDWAhECki/wdsAl72a1TGeCXjo5QKFy7M/PnzmThxIitXrqRevXp8/rmnRQaNyVm86Wh+F/grrjIXZ4DOqjrP34EZk7bMjVISEQYMGMCmTZsQEZo1a8b48eNtFrTJ0by5UgDIBVwBLqfjM8b4mW9GKTVs2JCdO3fSvn17hgwZQseOHe12ksmxvBl9NBKYC5TFVTJ7joiM8HdgxqTNd3WUihUrxqJFixg/fqqk1b4AABZbSURBVDwfffQR9erVY8uWLf4L3ZgQ5c1ynPuA+qp60dnOD2xX1dsDEF+qbOU1k7b0r/QWHx/PI488wsGDBxkzZgzDhg0jLMwukE32kamV13CN8wt32w4H/uuLwIzxv/SPUIqOjmbHjh089NBDjBgxgpiYGH788Ud/BmlMyPAmKVwEvhKRaSIyFdgN/CQir4rIq/4Nz5jMytgIpSJFijBv3jymTJnCpk2bqFOnjq3RYHIEb5LCcmA08CnwGfA88DHwlfNINxEZJCJ7ROQrERns7CsuIqtFZL/zXCwj5zbmehkfoSQi9O3bl/j4eEqVKsV9993H008/zaVLl3wbojEhxJshqdNTe6S3QRGpCfQFGgJ1gPYiUgUYDqxV1SrAWmfbmExKzwglzx3S1atXZ+vWrfTv359XX32Vxo0b8/XXX/stYmOCyZvRRzEisk1ETviodPbtwGeqelFVrwIbgI5AB+Ad55h3gAcz0YYxjvSMUEq5ZEZERASTJk1i6dKlHD58mPr16zNlyhSb02CyHW9GH30HdMHVl5Bc3sKpmJr+BkVuB5YCdwIJuK4K4oFHVbWo23FnVfWmW0giEofr/1aioqLqHzx4MCNhGHODirgSwY0q4Kq99JujR4/Ss2dPVq9ezYMPPsjUqVMpWbKk/0M0xkcyO/roCLBLVa+o6rWkR0aDUdV9uMpkrAZWAF8AXhe4V9UpqhqtqtGRkbasg/EV7zuky5Yty4oVK3jllVf46KOPqFWrFitXrvRveMYEiDdJYRjwgYj8RUQGJj0y06jTH3GHqt6Fq3TGfuC4iJQBcJ5PZKYNY9InfR3SYWFhDBkyhK1bt1KiRAliYmL485//TEJCgv9CNCYAvEkK/wdcA4riWnEt6ZFhIlLKeY4COuGaMb0MeNw55HFct5iMCZCMdUjXqdOBHTueZtCgQUyaNIk77riD7du3+zlWY/zHmz6F7apa36eNinwClMBVT2mIqq51SnK/h+tPs0O4Cu+l2qFtM5qNb83GNantEK5/hi+Scof0zTOkV68uRc+ePTlx4gSjRo1i+PDhhIeHY0yoSa1PwZuk8A9ghap+7I/gMsOSggm8iqTWIX3mzBn69+/PvHnzaNSoEe+++y5Vq1YNbIjGpCGzHc19gTUicsFHQ1KNycJS75AuXrw4c+fOZd68eXz77bfUrVuXiRMn2upuJsvwJimUBHIDRXD1JZQkk30KxmRd3nVIP/LII+zZs4eWLVsycOBAWrVqhQ2fNlmBNzOarwGdgWec12WAuv4OzJjQ5G2H9GzKlm3C8uUrOHeuOBUrbqFWrVpMnTrVJryZkObNjOZJwN3Ao86ui8BkfwZlTOjyZob0b7OjRZRChc4wbVoYw4aVIy4ujpiYGA4fPhzwyI3xhje3j5qo6hPArwDOiKA8fo3KmJDWHdcs50Tn+cYRSjeX6w4LS2DkyF+YNGkSmzZtokaNGkybNs2uGkzI8SYpXBGRMFwFYXCGjlqvmTEp8twZLXKY/v37s3v3burXr0/fvn1p27at9TWYkJJiUhCRpAHWbwALgUgR+T9gE64yFcYYj1LvjK5cuTJr167ljTfeYMuWLdSsWZM333zTRiiZkJDiPAUR2aGqdzivawCtcN1EXaOqewIXYspsnoIJTd4vAXry5AQSE0cQGfkrJ07kRfVFypR5OnChmhwptXkKqU23lKQXqprhBXWMyXmSfvGnPTs6MnIkTncdpUtf4pdfhrJs2TbatZtls6FNUKT2ry5SRIak9Kaq2lKcxqSoOzcngRvd3CFdoADUrj2fhg2/Zfr06dSrV89fARrjUWodzbmAgkChFB7GmEzx3CFdoYJw7NgxGjRowLBhw7h48aLH44zxh9SuFI6p6vMBi8SYHCcKT3WURKLYu3cnw4YN45///CcLFy5k8uTJtG7dOvAhmhwntSsFSeU9Y0ympTw7ulixYkydOpV169YRHh7O22+34eTJAqhev360Mb6W2pXCvQGLwpgcKe0O6ZYtW7J79wggjjx5km4jHUQ1DhFIu9/CmPRJs3R2KLMhqSZnqIin20yXL5chT56jAY/GZH2ZLZ1tjAkqzx3S4eHHGD58OL/88kuA4zHZWVCSgog8JSJficgeEZkrIvlEZLSI/CAiu5xHu2DEZkzo8TxD+syZArz88svUqFGDVatWAX/CdUdYnOc/BS5Ek20EPCmIyC3AQCBaVWviGvra1Xl7vKrWdR4fBTo2Y0KT5w7pkiX/xcaNG4mIiODbb9ui+hau5dRxnt/CEoNJr2BNmQwHIkTkCq5/7Udx3Tg1xtwk5Q7p5s1hx44d5MlTABFP/YNTgDcDFKfJDgJ+paCqPwDjcP3rPgb8rKqrnLcHiMiXIjJDRIp5+ryIxIlIvIjEnzx5MkBRGxNsKZfrjoiIIFcuzwNGVK9x7do1j+8Z40kwbh8VAzoAlYCyQAER6YHrWvf3uFZ1Owa84unzqjpFVaNVNToy0lYFNcYll8e9165BlSpVmDBhAufOnQtwTCYrCkZHcyvgf6p6UlWvAItwLeRzXFWvqWoiMBVoGITYjMmi4m7aowoHDrSlbNmyPPXUU5QrV45Bgwbx3XffBSE+k1UEIykcAhqLSH4REVyT5PaJSBm3YzoCIVGe25is4U2gH79dMeRCpB+33rqCTZs2sW3bNmJjY3nrrbeoWrUqDzzwAKtXr3ZWfpuNq0vPZkubIE1ecxbreQS4CuwE+gDTcN06Ulw3TZ9Q1WOpnccmrxmTPseOHWPy5MlMnjyZEydO8PTTZfj730+RO/cVt6M8r/1gso/UJq/ZjGZjcqBLly4xf/58WrfuS5kylz0cUQHX32YmO7IZzcaY6+TNm5fHHnuM0qWveHxf9RBXr14NcFQmFFhSMCYHE/E8W/rgQaVSpUq88MILHDuW6l1ck81YUjAmR7t5trRqfk6eHEz16tX529/+RlRUFJ07d2bt2rUkJiY6R1nndHZlScGYHK07rk7lCrhqJlVAZAoNGoxn5cqV7N+/n8GDB7Nu3TpatWpFtWrV+OCDbqj2xVW5VZ3nOCwxZA/W0WyMSdOvv/7KggULmDx5MrNmbaZiRU9HWed0VmGjj4wxPqMa5rHOkqogkujhEybU2OgjY4zPpNQ5feiQ0rFjR5YvX24jl7IwSwrGmHS6uXM6MTEfmze3Y8uWLbRv354KFSowcuRIK6mRBVlSMMak082d02Fh0/jDH5Zz+PBhFi1aRN26dXnppZeoUqUKd911F2+//TYXLlzARi2FPutTMMb4xQ8//MC7777L22+/zf79++nZMw+TJ18jb173Ut5WUiMYrKPZGBM0qsqWLVu47ba2lCzpaT1pG7UUaNbRbIwJGhGhadOmlCx50eP7iYkHmTp1Kj/99FOAIzOeWFIwxgSI51FLx46FExcXR+nSpencuTPLli3j8mX3In3WDxFIlhSMMQFy86glyE/ZsjPZunUrcXFxbNiwgQ4dOlCmTBn+9Kc/8c03o1CNw2ZPB471KRhjAmg2MBLXWltRuBLFb53MV65cYfXq1fz73/9m6dKl7N2bYLOn/cA6mo0xWc758+cpWLAwIje/Z7OnMyfkOppF5CkR+UpE9ojIXBHJJyLFRWS1iOx3nosFIzZjTGgoVKgQIhU8vnfwoHLnnXfy2muvcfTo0QBHlr0FPCmIyC3AQCBaVWviWlS2KzAcWKuqVYC1zrYxJkfzNHs6gi++6EJCQgKDBw+mXLlytGzZkokTJ3L48GG3I62DOiOC1dEcDkSISDiu/+JHgQ7AO8777wAPBik2Y0zI8DR7eiodOsxn165d7Nu3j1GjRnHixAkGDhxIVFQUDRo0YPHiziQm9sE6qNMvKH0KIjII158ACcAqVe0uIj+palG3Y86q6k23kEQkDtd/XaKiouofPHgwUGEbY0LYN998w5IlS1iyZAlz535mHdSpCKmOZqevYCHwCPAT8D6wAJjkTVJwZx3NxhhPUi7vDQsWvEfbtm0pXLhwECILDaHW0dwK+J+qnlTVK8AioAlwXETKADjPJ4IQmzEmG0ipvPeRI2F06dKFkiVL0rp1a15//XW+//57tyOsHyIYSeEQ0FhE8ouIAPcC+4BlwOPOMY8DS4MQmzEmW0h5otzGjRsZPHgwR44cYdCgQdx6661Ur16dOXPu59q13uT0fohg9Sn8H67bR1eBnUAfoCDwHq4ZLYeAzqp6JrXz2O0jY0zKUp8oB/Ddd9/x4Ycfsnz5cqZNW0MFjyNgs18/REj1KfiSJQVjjK+k1A+RmAh/+9tI7rvvPho1akR4eHgQovOtUOtTMMaYkJNSP8SJE3l56aWXaNasGZGRkXTu3Jnp06dz5MgRt6OyT1+EJQVjjAFS6ocoXXo6p06d4r333qNTp05s2bKFPn36UL58eWrWrMmsWe2yVV+E3T4yxphkafdDqCp79uxh5cqVrFy5MsW+CNUoREJzHpX1KRhjjJ+k1hfRo0c3WrVqRatWrYiK8nx7KhisT8EYY/wkpb6I06fzs27dOnr37k2FChWoUqUK/fr1Y8GCBZw+fZpQ7YewpGCMMZniuS8iMnIKR48eZffu3UyYMIFq1aoxa9YsOnfuzMCBJUlIeJRQ7Iew20fGGJNpafdFgGsRofj4eG6//T6KFv35pvfPni3M1q3v0bRpUwoWLOi3aK1PwRhjQkoYriuE6yUmQq5cEB4eToMGDWjRogUtW7b0eZKwPgVjjAkpKXU6l2flypUMHToUVWXcuHHExMRQtGhRGjduzLBhw9ixYyjXrpXHX30RWX9qnjHGZDkv4upDuOi2Lz9hYX+nTZs2tGnTBoALFy6wZcsWNmzYwIYNGzh27BVuuy2RXLmSPpPUFwGebldlhN0+MsaYoPCuH8JdYmIUYWGHPbyTvvpMqd0+sisFY4wJiu6k96/7sLAjKbxzKNPRJLfhszMZY4zxs5T6Inw3Mc6SgjHGZBme50S49vuGJQVjjMkyugNTcPUhiPM8BV91MkMQ+hRE5DZgvtuuysDfgKJAX+Cks/9ZVf0owOEZY0yIS39fRHoEPCmo6jdAXQARyQX8ACwGegHjVXVcoGMyxhjjEuzbR/cC36tqaNaXNcaYHCbYSaErMNdte4CIfCkiM0SkWLCCMsaYnCpoSUFE8gCxwPvOrreA3+O6tXQMeCWFz8WJSLyIxJ88edLTIcYYYzIomFcK9wE7VPU4gKoeV9VrqpoITAUaevqQqk5R1WhVjY6MjAxguMYYk/0Fc0ZzN9xuHYlIGVU95mx2BPakdYLt27efksytd1cSOJWJz2c1Oe37gn3nnMK+c/p4WEDUJSi1j0QkP3AYqKyqPzv7/o3r1pHiKuLxhFuS8Fcc8SnV/8iOctr3BfvOOYV9Z98JypWCql4EStyw79FgxGKMMeY3wR59ZIwxJoTk9KQwJdgBBFhO+75g3zmnsO/sI1l6PQVjjDG+ldOvFIwxxrixpGCMMSZZjk4KIvKCU1Zjl4isEpGywY7J30TknyLytfO9F4tI0WDH5G8i0llEvhKRRBHJ1sMWRSRGRL4Rke9EZHiw4/E3pyTOCRFJc15TdiAi5UVknYjsc/5ND/J1Gzk6KQD/VNXaqloX+BBXCe/sbjVQU1VrA98CI4IcTyDsAToBG4MdiD85VYffwFUtoDrQTUSqBzcqv5sJxAQ7iAC6CjytqrcDjYH+vv5vnKOTgqqec9ssgGviXLamqqtU9aqz+RlQLpjxBIKq7nNKtmd3DYHvVPW/qnoZmAd0CHJMfqWqG4EzwY4jUFT1mKrucF6fB/YBt/iyjWCWuQgJIvIi8BjwM3B3kMMJtD9y/YJHJmu7BVelgCRHgEZBisX4mYhUBOoBn/vyvNk+KYjIGqC0h7dGqupSVR0JjBSREcAAYFRAA/SDtL6zc8xIXJeiswMZm794851zAPGwL9tf/eZEIlIQWAgMvuGOR6Zl+6Sgqq28PHQOsJxskBTS+s4i8jjQHrhXs8lElXT8d87OjgDl3bbLAUeDFIvxExHJjSshzFbVRb4+f47uUxCRKm6bscDXwYolUEQkBngGiHVqUJnsYxtQRUQqOeuVdAWWBTkm40MiIsB0YJ+qvuqXNrLJH4oZIiILgduAROAg8KSq/hDcqPxLRL4D8gKnnV2fqeqTQQzJ70SkIzARiAR+AnapatvgRuUfItIOmADkAmao6otBDsmvRGQu0BJXGenjwChVnR7UoPxIRJoBnwC7cf3eAnhWVT/yWRs5OSkYY4y5Xo6+fWSMMeZ6lhSMMcYks6RgjDEmmSUFY4wxySwpGGOMSWZJwQSMiJRwKtLuEpEfReQHt+0tfmqznohMc15XE5FPReSSiAy94biQrrYpIo+LyH7n8XgKx6wPZBVYERkgIr0C1Z4JDBuSaoJCREYDF1R1nJ/beR8Yo6pfiEgpoALwIHDWvW0RuQu4ALyrqjX9GVN6iUhxIB6IxlW2YjtQX1XP3nDcemCoqsb7uP1cqnrNw/78wGZVrefL9kxw2ZWCCQkicsF5bikiG0TkPRH5VkReEpHuIrJVRHaLyO+d4yJFZKGIbHMeTT2csxBQW1W/AFDVE6q6Dbhy47HeVNsUkQdE5HMR2Skia0Tkd87+0SLyjrMmxwER6SQi/3DiXeGUJciMtsBqVT3jJILVpFwuurPzs/pWRJo78eUTkbedeHaKyN3O/p4iMsnt+30oIi2d1xdE5HkR+Ry40/nvsNdZh2McgDMj/oCINMzk9zMhxJKCCUV1gEFALeBRoKqqNgSmAX92jnkNGK+qDYCHnPduFI1rLQVf2QQ0dv4yngcMc3vv98D9uEpVzwLWqWotIMHZfx0R+YvbrTP3x+se2vVU/TSlcsnhzs9qML/V8eoP4MTTDXhHRPKl8V0LAHtUtRGwF+gI1HDW4Rjjdlw80DyNc5ksJNsXxDNZ0jZVPQYgIt8Dq5z9u/mtvHkroLqrFAwAhUWkkFNjPkkZ4KQP4yoHzBeRMkAe4H9u7/1HVa+IyG5cJSZWuMVc8cYTqeo/gX962W56qp8mFUjb7tZuM1xlPlDVr0XkIFA1jTav4Sq6BnAO+BWYJiLLcS1IleQEUC2Nc5ksxK4UTCi65PY60W07kd/+kAkD7lTVus7jlhsSArj+Sk/rL+L0mAhMcv7ifuKGc18CUNVE4Ipb9Vn3mJOl80ohPdVPk35W19za9ZRUwFU63f13gPv3+TWpH8FZlKkhriTxIL8lvKTPJKRwfpMFWVIwWdUqXOtfACAidT0csw+41YdtFgGSCiZ6HAHkLVX9p1tCc38M9HD4SqCNiBQTkWJAG2eftzYC3QFEpCoQBXwDHADqikiYiJTH9Yv/JuKq3V/EKbo2GHD/WVfFt7foTJBZUjBZ1UAg2un43AvcVOlVVb8GijgdzohIaRE5AgwB/ioiR0SksPPeXOBT4DZnf28PbY4G3heRT4BTfvlWHqjqGeAFXKWxtwHPO/u89SaQy7m1NR/oqaqXgM24boHtBsYBO1L4fCHgQxH5EtgAPOX2XlNgTTpiMSHOhqSabE1EngLOq6qnjmiTCSJSDxiiqo8GOxbjO3alYLK7t7i+j8L4TknguWAHYXzLrhSMMcYksysFY4wxySwpGGOMSWZJwRhjTDJLCsYYY5JZUjDGGJPs/wFUeOj5bI2o/QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xlabel('Time (11 am = 0 hours)')\n",
    "plt.ylabel('Temperature (degF)')\n",
    "plt.plot(t_A, numerical_A, color = 'black', label = 'Numerical')\n",
    "plt.plot(t_B, numerical_B, color = 'black')\n",
    "plt.plot(times1, analytical1, 'o', color = 'yellow')\n",
    "plt.plot(times2, analytical2, 'o', color = 'yellow', label = 'Analytical')\n",
    "plt.legend();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the non-linear Euler approximation varies when the time is before 11 am and the person has a higher temperature for the non-linear case and the non-linear approximation actually accounts for the variations in temperature which are physically real. Thus, the non-linear model is a better representation of what actually happens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 am Subtract: 30.909090909090907 min\n"
     ]
    }
   ],
   "source": [
    "for i in range(len(t_B)):\n",
    "    if numerical_B[i] <= 98.7 and numerical_B[i] >= 98.4:\n",
    "        print((int(t_B[i])+11), 'am', 'Subtract:', -(t_B[i]-int(t_B[i]))*60, 'min')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Part b  continued.` The time of death is approximately 9:29 am."
   ]
  }
 ],
 "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
}