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": 199,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K = 0.00458333333333 per min\n"
     ]
    }
   ],
   "source": [
    "Temp_t1 = 85\n",
    "Temp_t2 = 74\n",
    "Temp_ambient = 65\n",
    "delta_t = 120\n",
    "\n",
    "def measure_K(Temp_t1,Temp_t2, Temp_ambient,delta_t):\n",
    "    ''' Determine the value of K based upon temperature of corpse \n",
    "    when discovered\n",
    "    \n",
    "    Arguments\n",
    "    ----------\n",
    "    Temp_t1 = Temperature at point 1\n",
    "    Temp_t2 = Temperature at point 2\n",
    "    Temp_ambient = Constant ambient temperature\n",
    "    delta_t = Time in hours that have passed\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dff1 = dT/dt\n",
    "    K = an empirical constant\n",
    "    '''\n",
    "\n",
    "    dff1 = (Temp_t2 - Temp_t1)/(delta_t)\n",
    "   \n",
    "    K= -(dff1/(Temp_t1-Temp_ambient))\n",
    "\n",
    "print('K =', K, 'per min')\n"
   ]
  },
  {
   "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": 200,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Temperature at 3 hours = 66.5750000000054\n"
     ]
    }
   ],
   "source": [
    "Temp_t1 = 85\n",
    "Temp_t2 = 74\n",
    "Temp_ambient = 65\n",
    "delta_t1 = 120\n",
    "delta_t2 = 180\n",
    "\n",
    "def measure_K(Temp_t1,Temp_t2, Temp_ambient,delta_t,t):\n",
    "    ''' Determine the value of K based upon temperature of corpse \n",
    "    when discovered\n",
    "    \n",
    "    Arguments\n",
    "    ----------\n",
    "    Temp_t1 = Temperature at point 1\n",
    "    Temp_t2 = Temperature at point 2\n",
    "    Temp_ambient = Constant ambient temperature\n",
    "    delta_t1 = First amount of hours that have passed\n",
    "    delta_t2 = Second amount of hours that have passed\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    Temp_t3 = Temperature at 3 hours\n",
    "    dff1and2 = dT/dt\n",
    "    K = an empirical constant\n",
    "    '''\n",
    "\n",
    "Temp_t3 = Temp_t2 + (-K*delta_t2*(Temp_t2-Temp_ambient))\n",
    "   \n",
    "\n",
    "print('Temperature at 3 hours =', Temp_t3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K = 0.0048506944444436945 per min\n"
     ]
    }
   ],
   "source": [
    "\n",
    "dff1 = (Temp_t2 - Temp_t1)/(delta_t1)\n",
    "dff2 = (Temp_t3 - Temp_t1)/(delta_t2)   \n",
    "K= -((dff1+dff2)/2)/(Temp_t1-Temp_ambient)\n",
    "\n",
    "print('K =', K, 'per min')"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "A:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [],
   "source": [
    " \n",
    "Temp_t1 = 85\n",
    "Temp_ambient = 65\n",
    "K = 0.00458333333333\n",
    "e = 2.718281828459045\n",
    "\n",
    "def t_analytical(Temp_t1,Temp_ambient,t,K,e):\n",
    "    ''' Determine the value of T based upon time passed\n",
    "    \n",
    "    Arguments\n",
    "    ----------\n",
    "    Temp_t1 = Temperature at point 1\n",
    "    Temp_ambient = Constant ambient temperature\n",
    "    e = exp\n",
    "    K = empirical constant\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "  T = Temperature at time t\n",
    "   \n",
    "    '''\n",
    "\n",
    "    T=Temp_ambient+((Temp_t1 - Temp_ambient)*(e**(-K*t)))\n",
    "    \n",
    "    return T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    0 seconds after death, the Temperature of the body is 85 degrees F\n",
      "    1 seconds after death, the Temperature of the body is 85 degrees F\n",
      "    2 seconds after death, the Temperature of the body is 85 degrees F\n",
      "    3 seconds after death, the Temperature of the body is 85 degrees F\n",
      "    4 seconds after death, the Temperature of the body is 85 degrees F\n",
      "    5 seconds after death, the Temperature of the body is 85 degrees F\n",
      "    6 seconds after death, the Temperature of the body is 84 degrees F\n",
      "    7 seconds after death, the Temperature of the body is 84 degrees F\n",
      "    8 seconds after death, the Temperature of the body is 84 degrees F\n",
      "    9 seconds after death, the Temperature of the body is 84 degrees F\n",
      "   10 seconds after death, the Temperature of the body is 84 degrees F\n",
      "   11 seconds after death, the Temperature of the body is 84 degrees F\n",
      "   12 seconds after death, the Temperature of the body is 84 degrees F\n",
      "   13 seconds after death, the Temperature of the body is 84 degrees F\n",
      "   14 seconds after death, the Temperature of the body is 84 degrees F\n",
      "   15 seconds after death, the Temperature of the body is 84 degrees F\n",
      "   16 seconds after death, the Temperature of the body is 84 degrees F\n",
      "   17 seconds after death, the Temperature of the body is 84 degrees F\n",
      "   18 seconds after death, the Temperature of the body is 83 degrees F\n",
      "   19 seconds after death, the Temperature of the body is 83 degrees F\n"
     ]
    }
   ],
   "source": [
    "for t in range(0,20,1):\n",
    "    print('{:5.0f} seconds after death, the Temperature of the body is {:1.0f} degrees F'.format(t,t_analytical(Temp_t1,Temp_ambient,t,K,e)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 278,
   "metadata": {},
   "outputs": [],
   "source": [
    "t=('inf')\n",
    "\n",
    "\n",
    " \n",
    "Temp_t1 = 85\n",
    "Temp_ambient = 65\n",
    "K = 0.00458333333333\n",
    "e = 2.718281828459045\n",
    "\n",
    "def t_analytical(Temp_t1,Temp_ambient,t,K,e):\n",
    "    ''' Determine the value of T based upon time = infinity\n",
    "    \n",
    "    Arguments\n",
    "    ----------\n",
    "    Temp_t1 = Temperature at point 1\n",
    "    Temp_ambient = Constant ambient temperature\n",
    "    K = empirical constant\n",
    "    t = infinity \n",
    "    e = exp\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "  T = Temperature at time equals infinity\n",
    "   \n",
    "    '''\n",
    "\n",
    "    T=Temp_ambient+((Temp_t1 - Temp_ambient)*(e**-(K*t)))\n",
    "    \n",
    "    return T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "When t = inf Temperature = (55, 58, 60, 65, 66, 67)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "    print('When t =', t, 'Temperature =' ,T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hours passed since death =  1.8865228851474354\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "math.log(10)\n",
    "\n",
    "Temp_t1 = 85\n",
    "Temp_ambient = 65\n",
    "K = 0.00458333333333\n",
    "Temp_0 = 98.6\n",
    "t_1 = 11\n",
    "\n",
    "''' Determine the time of death, T = 98.6\n",
    "    \n",
    "    Arguments\n",
    "    ----------\n",
    "    Temp_t1 = Temperature at point 1\n",
    "    Temp_ambient = Constant ambient temperature\n",
    "    K = empirical constant\n",
    "    Temp_0 = Temp at time of death\n",
    "   t_1 = Time first given\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    " t = time in seconds since death\n",
    " t_hours = time in hours since death\n",
    " t_death = time in hours (decimal) since death\n",
    "   \n",
    "    '''\n",
    "\n",
    "t = math.log((Temp_0 - Temp_ambient)/(Temp_t1-Temp_ambient))/K\n",
    "t_hours = t/60\n",
    "\n",
    "\n",
    "print('Hours passed since death = ', t_hours)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time of Death = 9.113477114852564\n"
     ]
    }
   ],
   "source": [
    "t_death = t_1 - t_hours\n",
    "print('Time of Death =',t_death)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Answer__: The time of Death is 9:11 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": "markdown",
   "metadata": {},
   "source": [
    "A:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 227,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time= (800, 900, 1000, 1100, 1200, 1300)\n",
      "Temp = (55, 58, 60, 65, 66, 67)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "t = np.array1=(800, 900, 1000, 1100, 1200, 1300)\n",
    "T = np.array2=(55,58,60,65,66,67)\n",
    "print('time=', t)\n",
    "print('Temp =', T)\n",
    "'''Created arrays for both time of day in military time and the temperature at that time\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "If the time is 1100 hours, the Temperature is 65 degrees F\n"
     ]
    }
   ],
   "source": [
    "'''Change i to correspond with the number in the time array, starting at 0, to determine what temperature it was at that time'''\n",
    "\n",
    "i = 3\n",
    "print('If the time is', t[i], 'hours, the Temperature is', T[i], 'degrees F')\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAE0CAYAAAAR5SXAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU5fX48c9JQsK+JSyybwFEZbe4gCJUxWqr1g2pWlurte4bWlv112+trVqtG2qta8VdqXvFKm6AgrKLAiaEfU/CDiHb+f3x3Ewmw8wkk8xkJpnzfr3mdbn3PvfOSe4wJ/c+m6gqxhhjTLSlxDsAY4wxjZMlGGOMMTFhCcYYY0xMWIIxxhgTE5ZgjDHGxIQlGGOMMTFhCcYkBBHpJSJax1eveP8c9U1E3qrj7+xBv3Pt8LZdF8+fyTQelmCMMcbERFq8AzDGswZoFWLfGOC/3r8vB14MUW5vtINqACYS+v/xl8ARwLfAMSHKFMciKGPAEoxJEOqGlNgTbJ+I7PdbPaCqQcslI1UtCrVPRMq9f5bX5Hemqm2jFpgx2CMyY4wxMWIJxjQ6IjJMRP6fiMwUkW0iUuJVYC8Qkb+JSOdanneSX+X4odWUPcKv7DkB+w4RkXtEZKGI7BKRYhHZLCLfisgzIjKxNvHVVbhKfr/GBG9562O8bZtFZJ+ILBWRG0Skid8xWSJyl4gs98psE5FXRKRvDWLpICJ/FpF5IrJdRIpEZI2IPC8iw6P7k5uYUVV72SuhX8BYQL3XxdWUPd6vbKhXAXBMLeJoDuz2znFXNWXv8crtAJr6bR8JbK9BjGlR+L0t8s61qIbld3jlrwuy7y1v31vAVUBZiLhfAwQYCKwOUSYfyA4Tx2nArjC/m3Lgxnh/Lu1V/cvuYExjUw7MAm7ANQ7oB2QCg4Bf4Sq82wPTRCSiOgdV3Qf8x1udJCISrJy3/XxvdZpWrSd5FmgLrAcuAbK9+PoD44A7gZWRxFXPhgEPAu/hGg5k4X63L3v7zwEmUfl7ugDoCnQGLgOKcD/vQ8FOLiLH45JYK2A+cC7QwztmFPAKLoHdJyJnR/dHM1EX7wxnL3tV9yKCO5ganKspsNg71/W1OP7HfrGMqUG8Y/22d/fbfnw9/N5icQejwEuABOxPAZZ4+0uArUCXIOe5icq7kEMC9qUCed7+WUCTEHE+4pVZTRTu9OwVu5fdwZikou5u4lVv9cRanOITYKP371+EKHOBt1wHfO633b/V5oZavHciKMU9nqoykZSqluMej4H7Oe9R1Y2BB+PuQMDdhRwZsO9nQG/v3xerakmIGO7w4ugJHBtZ+KY+WYIxjY6IpIjIRBH5j4isFpG9/r3Xgbu8ogMiPbf3RfqSt3qOiKQHvHcGcJa3+lLAF/EaXP0LwJMiclik758AlqjqphD7/B/tfRisgKqup7LvTWBjix97yzxgs4i0DPbC3SGt9sqOjPgnMPXGEoxpVESkDe6u4WXgTNxfuc1DFG9Ty7d5wVu2B34SsO80XB2LfznAl5xu8FbHAktFJNdrOfYrEeley3jqU7C7kgr+/ZVCJSH/cs0Ctlck/D64xhThXv28sh2qD9nEiyUY09j8ExiNe8b/BHAK7gsrE1dx3Aq41Stbq47GqroY11gAKh+HVah4bLZIVZcGOfY5YALwBa4eoS+u8cEzwBoR+UhEBtcmrnpSFsVygY0kapPwm9biGFNPrCe/aTREJAvX6gjgj6p6d4hygX8518YLuKbIp4lIG1Xd6bVK+4nf/qBU9UPgQxFpj2uJdSwu6QzFPSb6UkRGqep3UYizIakYbeBjVa1N/ZhJMHYHYxqTwVR+pl8OU+6IKLzXi7i7pAygornsOd56eTXvD4CqFqrqe6p6q6oOwz1eKwFaAJOjEGNDk+cth4RqAm4aFkswpjHJ8Pt3arAC3l3OhLq+kapuAD7zVi8IWH4SogVVded8H/jKWw07UkAj9T9v2QGXbE0DZwnGNCZ5fv8+PXCniKTg6mii8YgMYKq3PF5EjsF17PTfHvj+nb1GCEF5w6z08FYLohRjQ/IfKluHPSoiPcKURUT6251OYrMEYxoNVV2B60QJcJeI3Cwi2d6YWOOBj3BNiL+P0ltOw7WIEtwjMwH8e/sHOgpYLyL/FpFzvS/IdiLSTUROwk1J0MsrW+0jtsbG6/dyEe4xYXdggYjcKiKDvd9TR2+cuctE5ANgGSHuVE1isEp+09hcAnyKay12j/fy9yQuwTxQ1zdS1d0i8jZuTpZe3ua3NfzQ+C1xX6IXhSnzmKoGvQtq7FR1poicgkuwHYC/eq9gSnAt8UyCsjsY06io6nxgBPA8ri9GCbAF1/HvHFW9LMpvGZgIQrYeA6bj6n/+DszGdbw8gBufayXuLmisql4Z5RgbFFWdgWu+fSPuj4VtuJ77+4Bc4A1cgu6oqjVtNm3iQAJGfDDGGGOiwu5gjDHGxIQlGGOMMTFhCcYYY0xMWIIxxhgTE9ZM2ZOVlaW9evWKdxjGGNOgzJ8/P19Vg45qbQnG06tXL+bNmxfvMIwxpkERkTWh9tkjMmOMMTFhCcYYY0xMWIIxxhgTE5ZgjDHGxIQlGGOMMTFhCcYYY0xMWIIxxpgkVVJWzterCnni85UxOb/1gzHGmCShqqzctpdZOduYmZPPnLwC9ha7GQ8mHN6Znpktovp+lmCMMaYRK9hzgNkrC5j5wzZm5eazaWdR0HIzc/ItwRhjjAmtqKSM+Wu2MzMnn5k52/hu466w5bu2bcZx/bM49JBWUY/FEowxxjRgqsryzbuZlZPPFznb+GZ1IUUl5SHLt8pI4+i+mYzJzmJ0dgd6ZTZHRGISmyUYY4xpYLbuKmJmTj6zcvOZmZNP/p4DIcumpghDu7dlTHYWY7KzGNKtLWmp9dO+yxKMMcYkuH3FpcxdVcisnHxm5eSzYsvusOV7Z7Vwdyj9sjiqbyatmzapp0irsgRjjDEJprxcWbpxp7tLycln/prtFJeFfuzVtnkTju2b5T32yqJbu+b1GG1olmCMMSYBrN++j1k5+czMzefL3Hy27ysJWbZJqjCiZzvGZHdgTHYWh3VpQ2pKbOpR6sISjDHGxMHuohK+WlnArFx3l5KXvzds+f6dWjImuwOjs7MY1bs9zdMT/+s78SM0xphGoLSsnMXrd/geey1ct4Oycg1ZPqtlhq8eZXR2Fp1aN63HaKPDEowxxsSAqrK6YJ+v1/xXKwvYfaA0ZPmMtBRG9clkjJdQBnZuFbPmw/XFEowxxkTJjn3FzM4tYFauSyrrt+8PW/6wLq199SgjerajaZPUeoq0fliCMcaYWiouLWf+mu3Myt3GrJx8lmzYiYZ+6sUhbZr6Ojge2zeTzJYZ9RdsHFiCMcaYGlJVcrbu8epRtjF3VSH7vMEig2mRnsrRfTO9epQO9O3QosE/9oqEJRhjjAlj2+4DzPZ6zM/K3caWXaF7zacIDO7WluO8u5RhPdrSpJ56zSciSzDGGOOnqKSMr1cV+oZhWbYp/GCRPdo39w3DcnSfLNo0j0+v+URkCcYYk9TKy5Vlm3f5mg9/vbqQ4tLQveZbN03jmL5ZjOmfxZh+HeiRmRi95hNRRAlGRFoAY4BhQCegLbAd2AosAGapavjeQsYYE2ebdu73JZTZufkU7C0OWTYtRRjeo51vGJbB3domZK/5RFRtghFXI3UacAUwzu8Y/99wRbuJUhGZATwGvK8arj2FMcbUj70HSpmTV+AbgTh3656w5ft2aOFrPjyqTyYtM+xhT22E/a2JyJnAvUAfXELZBXwNfA8UeuutgUxgEPAjYAJwMrBSRG5W1bdiFr0xxgRRVq58u2EnM3/YxszcfBau3U5JWei/d9u3SOfYflm+nvNd2jarx2gbr5AJRkQ+AY4HdgBPAP8Gvg53V+Ld7YwCLgbOBaaJyOeqOi6aQRtjTKB1hft8szh+ubKAnftDDxaZnpbCkb3cYJGj+2Ux6JDWpNhjr6gLdwczBPgD8LCqhu+O6vGSzxxgjohcD1wL3FTnKI0xJsDO/W6wyJk5bq75NQX7wpYf2LkVx/V3CeXIXu1plt64es0nonAJpreqhm+fF4aXlO4Wkcdqew5jjKlQUlbOonU7fHcpi9ftIMxYkXRslcHo7CyOy+7Asf2y6NCqcfeaT0QhE0xdkksszmOMSS6qSl7+XjdHSs425uQVsifMYJHNmqQyqk97X+V8dseWSdVrPhFZ0whjTMIo3FvMbG9+lJk529i4syhkWRE4omsbRvfLYkx2B4b3bEtGmj32SiThKvmfB75Q1afqMR5jTBI5UFrG/NXbmZnrEsp3G3eFHSyya9tmXq/5DhzTN5N2LdLrL1gTsXB3MBcApcBBCUZEdgFTVfXKWAVmjGl8VJUVW3YzKyefL3Ly+XpVAUUloXvNt8pI46i+mb6xvXplNrfHXg1IbR+RtQRi0lBcRJoBVwPnANlAOrAFmAc8qKqz/co+B/wyzOlWqOrAWMRpjKmZrbuKfNMCz8zNZ9vu0INFpqYIQ7u3ZXS/LI7rn8WQbm1JS+LBIhu6hKqDEZHewP+AfrjhZz4HDgC9gNOBxcDsIIfOBnKDbN8Uk0CNMSHtLy5j7qoCrx4lnxVbdoct3zurhVePksVRfTNp3dQGi2wsEibBeOOcfQT0Be4E7lTVEr/9mbgRA4J5SlWfi3mQxpiDlJcr323cxUxv0q15q7dTXBb6sVfb5k04tq8b12t0vyy6t7fBIhurhEkwwG245PK8qt4RuFNVC4CCeo/KGHOQDTv2++aan52bz/Z9oXvNN0kVRvRs52s+fFiXNjZYZJJIiAQjIunApd7q3fGMxRhzsN1FJczJK/Qllbz88IOm9+/UktH9OjCmfxajereneXpCfNWYelbdVR8lIg/XYp+q6rURxDEC9/hrnaouE5FjcCM4ZwKbgemq+lWY408QkcG4xgdbgFnAR6oa+j7dGBNW7tbdvL9kM7Nyt7Fw7Q5Kw3Sbz2qZweh+mYz2xvbq3KZpPUZqElV1CeZQ7xXpPsWNQ1ZTR3jLnBAtw+4QkWnAhSHGRbsoyLbvRWSiqn4bQRzGGGD60s1c9dKCkEklIy2FH/Vu7+uTMrBzK2s+bA4SLsHcX29RQHtveRyQCtwH/BNX53Icbn6Zs3DTA/za77hFwHxgBrAGN3XAcOAu3GCdH4vIcFXdEOxNReQy4DKAHj16RPcnMqaBWle4j8mvLz4ouRzWpTWjs90sjiN7taNpE+s1b8KTRJgTTET+gEsK4FqEXRqwfyRuHhoFslU1r5rzpeOaOB8FPKqqV1UXw8iRI3XevHm1Cd+YRqO4tJxznviKxet2AK7n/M0TBnBsvyyyWtpgkeZgIjJfVUcG25coPZj8G8o/GbhTVefh7lRSgLHVnUxVi4G/eas/iUJ8xiSF+/63wpdc0lKEKZOGcfrQrpZcTK0kSoJZ7ffvVSHKVGzvXMNzLveWXWsTkDHJ5tMVW/nXF5UPByafPIBhPdrFMSLT0IVMMCIyKRpvUMPzLPD7d6jOlFneMvxk2gefp6bljUlam3cWceNri33rYwd04NIxfeIYkWkMwt3BvCAii0TkLBGJqBG7iKSJyNkishiYWl15rxJ+rrc6Psj52uEq78GNSVYT53rLb2pY3pikVFauXPvKQgr3FgPQqXUG958zxKYQNnUWLsFcinsc9RqwUUQeFJGfikiHYIVFpKOI/MzrG7MReBXoCPymhrFUVPLfISJD/c7bFHgcaIOrh/nK2z5URE4TkSpNWbzkdgNwjbfpgRq+vzFJ6eEZOcxdVQhAisCD5w0j0+pcTBSEm9HyaRF5DbgV15T3Gtwox4hIIa4J8S5c0+BMKpsaC5AP/BW4V1XDj3RX+X7vish9wE3AXBGZ673Hj4AuwAbgfK1s9tYLeBMoFJEfgPVAK1yfmi5AOXCLqn5Yk/c3Jhl9tbKARz7J8a1fPS6bo/uGekptTGTCPvryksMfROTPwPnAGcBoXEIJ/BQW4JoGvwW85rXkioiqThaRL3GJbBjQHFgL/AO4W1W3+RVfDDyES0A9vfKKSzTP4ponz480BmOSRcGeA1z7ykLfvPajerfnmvHZ8Q3KNCo1qltR1SLcl/az4rrr9sQ9/moD7AC2quqaaASkqm/i7kyqK7cKuC4a72lMsikvV258fTFbvblZ2rdI56GJw2wQShNVEY9A5z2iWk3VpsXGmAbkyZl5fLai8oHA/ecOsfHDTNQlSj8YY0w9WbB2O3//cIVv/bLj+nDCgI5xjMg0VhHfwXhTGp8LHA10AD5X1Ye9fb1wFezzalMHY4yJrZ37S7jm5YW+ccaGdG/LTScNiHNUprGKtH/L8cAruPoXwVWq7/QrMhr4Ny4BTYtSjMaYKFBVfj9tCeu3uwHJWzVNY8r5w0hPswcZJjZq/MkSkQHAe7jk8jxwCS7J+HsLOIBrbWaMSSAvzF3LB0s3+9bvOWuwTVdsYiqSO5g/4poN/0JVXwEQkaf9C6jqHhFZhmsybIxJEN9t3Mmd733vW7/gqB785IhD4hiRSQaR3BuPA5ZUJJcw1gH2yTUmQew9UMrVLy2kuNRN8DqwcytuO3VQnKMyySCSBNMBWFFtKSjF3ekYYxLA7W8vJS9/LwDNmqQyZdJwmyzM1ItIEswOajb0fR9gW7WljDExN23+ev6zoHJC1zvPOJx+HVvGMSKTTCJJMPOAkV5T5KBEZAhuquIv6xaWMaauVm7bw+1vL/Wt/3x4V84e0S2OEZlkE0mCeQLIAF4TkYMmsBeRQ4Cn/coaY+KkqKSMK19cwL7iMgD6dGjBnacfHueoTLKpcYJR1XeAZ4CRQK43KCXAWBGZAazEzdnyL1X9NOqRGmNq7C/vf8/yzW4g8/S0FKacP5wWGRH3qzamTiLqYaWqvwFuBnYDR3mbewEn4Cr3b1PV30UzQGNMZD74dhMvzFnrW7/91EMZ1KV1HCMyyao2g13e500qdhSuQj8V1zR5lqrui3J8xpgIrCvcx83TlvjWTzm8Mxcc1TOOEZlkVuMEIyLjgFJV/cIbZ+wL72WMSQAlZeVc/fJCdheVAtCtXTPuPmswboYNY+pfJI/IPgbuiFUgxpi6ue/DFSxatwOAtBThkfOH0aZZkzhHZZJZJAmmENgSq0CMMbX36YqtPPFFnm998skDGNajXRwjMibyfjADYxWIMaZ2Nu8s4sbXFvvWxw7owKVj+sQxImOcSBLM/cBQEZkYq2CMMZEpK1eue3UhhXvd9EsdW2Vw/zlDSLGpj00CiKQVWQHwGPCiiJwFvAmsAfYHK6yqC+oenjEmnEc+yWFOXiEAKQIPTRxGZsuMOEdljBNJgpmHm2BMgJ97r1A0wnMbYyL01coCHp6R41u/elw2R/fNjGNExlQVSRJYgEscxpg4K9hzgGtfWYg38zGjerfnmvHZ8Q3KmAA1TjCqOjKWgRhjaqa8XLnx9cVs3X0AgPYt0nlo4jBSrd7FJBibjNuYBuapWXl8tqJyRoz7zxlC5zZN4xiRMcFZgjGmAVm4djv3Tq+c9++y4/pwwsCOcYzImNAiGSpmeCQntlZkxkTXzv0lXP3yQkq9ipch3dty00kD4hyVMaHVphVZTVgrMmOiSFX5/bQlrN/uegW0aprGlPOHkZ5mDyFM4opGK7IUoCfQ3tu/GCire2jGmAovzF3LB0s3+9bvOWsw3ds3j2NExlQvaq3IRGQU8CSwHjijjnEZYzzfb9zFne9971v/xage/OSIQ+IYkTE1E7X7a1Wdi0ss44GbonVeY5LZ3gOlXPXyAopLywEY2LkVt582KM5RGVMzUX2Aq6p5wDfAxdE8rzHJ6o63vyNv214AmjVJZcqk4TRtkhrnqIypmVjUEObjplE2xtTBtPnrmbZgvW/9zjMOp1/HlnGMyJjIRDXBiEhzYBSwJ5rnNSbZrNy2h9vfXupb//mwrpw9olscIzImcpH0g2kfZndL3FwxtwJdgNfqGJcxSauopIwrX1zAvmLXGLNPVgvuPOPwOEdlTOQiaaa8rfoiCG7Wy1trF44x5q73l7F8824A0tNSeGTSMFpkWLcy0/BE8qndTuiOlsXABmAG8ICqbq1rYMYkow++3cTUOWt867efeiiHdWkTx4iMqb1I+sFkxTIQY5LdusJ93DxtiW/9lMM7c8FRPeMYkTF1Y+NMGJMASsrKufrlhewuKgWga9tm3H3WYERsCH7TcNU4wYjIOyJyTQ3KXSUi79QtLGOSy30frmDRuh0ApKUIj0waRptmTeIclTF1E0kdzGm4Pi7VGQKcWrtwjEk+n67YyhNf5PnWJ588gOE92sUxImOiIxaPyNKB8hic15hGZ8uuIm58bbFvfeyADlw6pk8cIzImemKRYIYBBTE4rzGNSlm5cu0rCyncWwxAx1YZ3H/OEFJs6mPTSIR9RBakLmV8mPqVNGAAbpiYN+semjGN2yOf5DAnrxAAEXhw4lAyW2bEOSpjoqe6OpjT/P6tQHfvFc4K4Ja6BGVMYzcnr4CHZ+T41q8el80xfa0ngGlcqkswP/WWArwDfAw8FKJsMbBBVb8Psd8YAxTsOcC1ryzEm/mYUb3bc+347PgGZUwMhE0wqvp+xb9FZD4w03+bMSYy5eXKja8vZsuuAwC0b5HOQxOHkWr1LqYRiqQn/5GxDMSYZPDUrDw+W1E5rN/95wyhc5umcYzImNixnvzG1JOFa7dz7/QVvvVLx/TmhIEd4xiRMbEV8RCt3rD9pwL9gda4+plAqqrX1jE2YxqNnftLuPrlhZR6FS9Durdl8skD4xyVMbEVUYIRkUuBB4Bm/pu9pfqtK2AJxhhAVbn1P0tYv30/AK2apjHl/GGkp9kDBNO4RTIW2QnAP73VR4BvvH/fCDwObPLWHwWqHbPMmGTx4ty1/Pfbzb71e84aTPf2zeMYkTH1I5I/oa73lhNU9TpgGYCqPqCqV+Iemb0M/AJ4L6pRGtNAfb9xF39+r7Ll/i9G9eAnRxwSx4iMqT+RJJgfAQtVdVawnaq6D7gEOAD8XxRiM6ZB23uglKteXkBxqRuab2DnVtx+2qA4R2VM/YkkwbQDcv3WSwBExHevr6oHgFnA+KhEZ0wDdsfb35G3bS8AzZqkMmXScJo2SY1zVMbUn0gSTAHQ0m+90FsGTrmXBmTWJShjGrr/LFjPtAXrfet3nnE4/Tq2DHOEMY1PJAlmDdDDb30JrsXYWRUbRKQdcDywLirRGdMArdy2h9veWupb//mwrpw9olscIzImPiJppvwJcIuIdFXVDcD7wE7gTyLSG1gPTATaAM9GPVJjGoCikjKufHEB+4rLAOiT1YI7zzg8zlEZEx+RJJhXgT641mIbVHWniPwO+DfwK79yy7BKfpOk7np/Gcs37wYgPS2FRyYNo0VGxP2ZjWkUIhmLbAlwfsC2V0TkG+AMoD2wHHjNq+w3Jql88O0mps5Z41u/7dRDOaxLmzhGZEx81TjBiEgKbggY9d+uqiuB+6MdmDENybrCfdw8bYlvfcJhnbnwqMD2L8Ykl0gq+Uup7L1vjPGUlJVz9csL2V1UCkDXts245+zBiNgQ/Ca5RZJg9uAegRlj/Nz34QoWrdsBQFqK8MikYbRp1iTOURkTf5EkmBWAjXFhjJ9PV2zliS/yfOs3nTyA4T3axTEiYxJHJAnmGWC0iFibS2OALbuKuPG1xb714/p34LIxfeIYkTGJpcYJRlUfB14EZojIlSLSJXZhGZPYysqV615ZROHeYgA6tsrgH+cOIcWmPjbGJ5Lh+ncBZwNZwMPAOhEpEpFdQV47axuQiDQTkZtF5BsR2SEi+0RklYi8LiLHhjhmkojMFJGdIrJHROZ5SdAm3DAxMeWTXL7KKwBABB6cOJSslhlxjsqYxBJJD7BgAymle6+o8EYE+B/QD9gKfI4bnbkXcDqwGJgdcMyjwBVAETADNwjneGAKMF5EzlHVsmjFaMycvAIemvGDb/3qcdkc0zcrjhEZk5giSTCtYhYFICItgI+AvsCdwJ2qWuK3P5OAQTRF5CxcctkMHKeqOd72TsCnwJnAVcBDsYzdJI+CPQe49pWFeDMf86Pe7blmXL/4BmVMgoqkJ//eWAYC3IZLLs+r6h1B3r8AN6Kzv1u95S0VycUru8UbxuYz4Pci8oiqlscmbJMsysuVm15fzJZdbqCKds2b8PDEYaSl2pNYY4JJiP8ZIpIOXOqt3l3DY7oBI4Bi4PXA/ar6ObAB6AwcFZ1ITTJ7etYqPl2xzbd+/7lD6NymaRwjMiaxRTwKn/fFfglwNNABeL/ijkNEhgMDgXdVdXcEpx2Be/y1TlWXicgxwGnets3AdFX9KuCYYd7yO1XdH+K83wBdvbJfRhCPMVUsWreDe6ZX9jO+dExvxg3sFMeIjEl8ESUYEZkIPAU0w80Fo7iK9wo9ganAxd6ypo7wljki8hzwy4D9d4jINOBCv2TS21uuIbS1AWWNidiXuflMfmMJpV7Fy5DubZl88sA4R2VM4oukmfKRwPO4xPJnXEutwEb//8UNKXN6hHG095bHARcB9+FakrXzzrUBN7HZo37HVLRqC1c3tMdbxrSBgmmcvt+4i4ue+ZpJT81lww73d02rjDQemTiM9LSEeLpsTEKL5A7m97iEdKqqfgYcNJifqh4QkeXAYRHGUfG/NQ14SlUn++17R0Q2Al8DvxSRv6hqHpXJrcrozpEQkcuAywB69OhRTWmTLNYV7uMfH/3AW4s24D92eEZaCg+cN5Qemc3jF5wxDUgkCWY08E1FcgljPe7uJhL+9TVPBu5U1XkiMh8YCYwF8vyOCTfRecW+oPVBqvov4F8AI0eOrHWiMo3D9r3FPPppLs9/tYbisspGhykCZw3vxvUn9qdL22ZxjNCYhiWSBNOW8PUdFVKIvPPlar9/rwpRZhUuwXQOOCbcpBvdg5zfmCr2F5fx7JerePyzlb4h9yuMH9iRmycMZEBne8pqTKQiSTD5uB711ekPbIowjgV+/84EtgUpU9FVuqJeZaG3PExEmoVoSXZkQFljfErLypm2YD0PfJTD5l1FVfYN7d6WW08ZyKg+mSGONsZUJ5IE8yVwhogM9qZPPoiIHA8cCjwXSbH962oAAB9RSURBVBCqukFE5gKjcI/Xqsw7IyLtgOHe6jzvmHUissDbfg6uAUJgLN1wzZwDmzibJKaqfLxsK/dOX07O1j1V9vXJasHNEwZw8mGdbcIwY+ookgTzMK4l139E5CJVrdKvxOsD8wxQTtXWXjV1F/AOrknybFVd5J23KfA40AaYT9Vk8TdcJ8t7RORLVc31jukIPOaVudt68ZsK89ds5+4PlvHN6u1VtndolcF1P87m3JHdaWI9842JikiGipkpIn8C/gTMFJEtuBZcp4nIStzjMwFuV9X5kQaiqu+KyH3ATcBc746mAPgR0AXXVPl81cp2Par6hog8DvwO+FZEPqZysMvWwFu4QS9Nklu5bQ/3Tl/Oh99tqbK9ZUYavz2uD5eM6U3z9Ij7HRtjwojof5Sq/llEvsP1gznU25zlvfKAO1T1pdoGo6qTReRL4Gpc7/vmuM6S/8DdiRxUN6OqV4jILOBK4HggFfeI7Rngcbt7SW5bdxXxwMc5vDZvHWXllQ0Fm6QKvxjVk6vH9SPThtk3JiZEtXatc0WkB9AH94W+TlV/qOaQhDZy5EidN29evMMwUbK7qIQnPs/j6Vmr2F9SdbaGnw3pwo0n9adnZos4RWdM4yEi81V1ZLB9tX4moKprqRyKxZiEcKC0jBfnrGXKp7m+2SYrHNsvk99POJQjurWJU3TGJJdaJxgRaY0bSFKBTapa61ksjamr8nLl3SUbue9/K1hXWLXF+qBDWvP7UwYyJjvLWoYZU49qM5ryBcD1wBD8hmsRkcXAQ6r6fMiDjYmBWTn53D19GUs37KqyvVu7Ztx00gB+NqQLKSmWWIypbzVOMOL+9HsemERlYqm4a2mDq5R/VkROBi7Q2lbuGFNDSzfs5J7py5mZk19le7vmTbhqXDYXHNWDjLTUOEVnjInkDuYK4Be4XvZ3AlMrHot5j8suxM1KORHXV8WaB5uYWFe4j/v+t4K3F22ssr1pkxQuGd2b3x7fl9ZNm8QpOmNMhUgSzGXAAWCsqi7z36Gqu4BHReQT3LAsl2EJxkRZ4d5iHvkkhxfmrKGkrPIGOUXg3JHdue7H/W2GSWMSSCQJpj/wWWBy8efNRvkprj+KMVGxr7iUZ2at4onP89h9oOpglCcO6sQtEwbQr6MNRmlMookkwewCtldbCnYQYnh8YyJRWlbO6/PX88BHP7B194Eq+0b0bMetpwxkZK/2IY42xsRbJAnmM+AYEUlT1dJgBUQkDTga+DwKsZkkpar87/st3Dt9OSu3VZ2wtG+HFtwyYSAnDupkTY6NSXCRJJjbgG+AJ0XkGlWtcpciIi1xA2K29coaE7F5qwv52wfLmb+m6s1yp9YZXP/j/pw9ohtpNhilMQ1CJAnmp8AbwK+B00XkPSonB+vl7W+DGwPstCDTKf+jrsGaxitny27umb6Cj5dVHYyyVUYal4/ty6+P7U2zdGtybExDEkmCuQ/Xax/cXcoFfuv+2eTXAceJV84SjDnI5p1FPPDRD7w+fx1+Y1HSJFW48KheXDWuH+1bRDpBqjEmEUSSYP5BZUIxpk527i/hic9X8szsVRSVVB3w+oyhXbjxpAF0b988TtEZY6IhkvlgboplICY5HCgtY+pXa5jyaS479pVU2TcmO4tbJgzk8K42GKUxjYHNsGTqRXm58vbiDdz34Q9s2FF1MMrDu7bm9xMOZXR2VpyiM8bEgiUYE1Oqyhc5+dz9wXKWbao6GGX39m4wyp8OtsEojWmMajOa8hDgBNw0xqHG5VBVvbYugZmG79v1O/nbB8v4cmVBle3tW6Rz9bh+/GJUT9LTrMmxMY1VJKMpNwVeAM6s2BSmuAKWYJLUmoK93Pe/H3h3cdXBKJs1SeU3Y3pz2XF9aGWDURrT6EVyB/NX4OfAHuA1IMf7tzEA5O85wJRPcnlxbtXBKFNThPOO7M5147Pp2NoGozQmWUSSYM7Dzf8yQlXzYhSPaYD2Hijl6VmreOLzlewtLquyb8JhnZk8YQB9O7SMU3TGmHiJJMG0Az6y5GIqlJSV8+o363jw4xzy91QdjPLIXu34/SmHMqJnuzhFZ4yJt0gSzCqso6XBtQybvnQzf/9wBXn5VQejzO7YklsmDGT8oR1tMEpjklwkCeZ54FYRyVTVgmpLm0Zp664irn55IXNXFVbZ3rl1U244sT9njehGqjU5NsYQ+VhkxwMfi8jlqjo3RjGZBLV43Q5+O3U+m3cV+ba1aprGFWP78atje9G0iQ1GaYypFMlQMWUiMhGYBXwpIruB9UB58OI6JEoxmgQwbf56bn3zW4pL3eVOEfj1sb258oR+tLPBKI0xQUTSD6YrbiKx3rg+MK2BQSGKW11NI1FaVs7dHyznqVmrfNtaN01jyqThHNe/QxwjM8Ykukgekf0d6AN8DTwC5GL9YBq1HfuKufrlhczMyfdty+7YkicvGkmvrBZxjMwY0xBEkmB+DKwDTlDV/dUVNg3bD1t2c+nz81hTsM+37cRBnXjgvKG0zLAh7Iwx1Yvkm6Ip8Kkll8bvw+82c8Ori6p0mrxmfDbXjc+2QSmNMTUWSYL5FsiMVSAm/srLlUc+yeWBj3/wbWuensr95wzhlCMOiWNkxpiGKJIE8wDwsogMV9UFsQrIxMfeA6Xc+Npipn+32betW7tmPHnRSA49pHUcIzPGNFSRJJhPgPuBGSJyN/AhoZspo6qFwbabxLO2YB+XPj+PFVt2+7Yd0zeTKZOG096aIBtjaimSBLPNWwpuZOW/himrEZ7bxMns3HyufGlBlemLf3VsL/74k0NJS7W5WowxtRdJEtiO9W9pNFSVZ2ev5q7/LqOs3F3W9NQU/nLm4Zw7snucozPGNAaR9OS3CdMbiaKSMm57aylvzF/v29axVQb/vHAEw3vY6MfGmOiwx1hJZsuuIn47dT6L1u3wbRvSvS3/unAEnWwyMGNMFFmCSSIL127nt1Pns3V35dwtZw3vxl1nHm4DVRpjoi7iWlwRGSkiz4rIchEpEJHH/PadICJ/EBEbpCrBvD5vHec9MceXXFJThNtPG8R95wy25GKMiYmI7mBEZDKu9VjFN5LievhXaALcCRQAT0QjQFM3pWXl3PXfZTw7e7VvW5tmTXh00nBGZ1u1mjEmdmp8ByMiJwH3APnAr4C+uCbL/j7GtTb7abQCNLW3fW8xv3z26yrJpX+nlrxz1bGWXIwxMRfJHcwNQDEwQVUXAwdNiauq5SLyA9A/ahGaWlm+eReXPj+PdYWVQ8edfFgn7j/XBqs0xtSPSL5pjgTmViSXMNYDR9Q+JFNX05du4obXFrPPb7DK63/cn6vH9bPBKo0x9SaSBNMC2FxtKWjOwY/OTD0oL1cenJHDwzNyfNtapKfyj/OGcvJhneMYmTEmGUWSYDZTs0dfA3Hzxph6tOdAKde/uoiPvt/i29ajfXOevGgkAzq3imNkxphkFUkz5c+BwSJyXKgCInIGbtbLGXUNzNTcmoK9/Pyx2VWSy+h+Wbxz1bGWXIwxcRNJgrkPKAPeFJFJItK8YoeIpInIz4EngQPAQ9EN04QyM2cbP5symx+2VM5efcno3jz3qyNp29xGQjbGxE8kY5F9KyKX4/q3TAVKcf1gJgLnA+m4oft/o6o5IU9kokJVeXrWKv7632V4Y1WSnpbC3848grNGdItvcMYYQ4QdLVX1GRH5DrgDGIerzG+KSzafAX9S1ZnRDtJUVVRSxh/e/Jb/LNjg29apdQZPXDiSod3bxjEyY4ypFHGHCFWdC5wqIulAF1yv/k2qui/awZmDbd5ZxG+nzmPx+p2+bcN6tOWJC0bQ0QarNMYkkJAJRkSeAWap6jPB9qtqMbA6RnGZIOav2c7lL8xnm99gleeO7MadZxxORpqNJ2aMSSzh7mAu9pZBE4ypX699s47b3lpKcZmboTo1RbjjtEFcdHTPg0ZUMMaYRGBjhiS4krJy/vLe9/z7qzW+be2au8Eqj+ln44kZYxKXJZgEVri3mCtenM+cvELftoGdW/HkRSPp3r55mCONMSb+LMEkqO837uKyqfNYv71ysMpTDu/MfecMoYUNVmmMaQDsmyoBvb9kEze9vpj9JZWDVd54Yn+uGtfP6luMMQ1GdQnmbBEZW4vzqqr2rcVxSa28XHng4x945JNc37aWGWk8cN5QThzUKY6RGWNM5KpLMC29V6S0Fscktd1FJVz/6iI+XrbVt61XphusMruTjSdmjGl4qksw03GzWJoYWpW/l0ufn0fu1srxxMZkZzHl/OG0ad4kjpEZY0ztVZdgNqvq5/USSZL6/IdtXP3SAnYVlfq2XXZcH24+eQBpqZGMRWqMMYnFKvnjRFV5cmYed3+wvMpglfecdQRnDrPBKo0xDZ8lmDgoKinj99OW8Naijb5tnVs35V8XjWBwNxus0hjTOFiCqWcbd+znt1Pn8+2GysEqR/Rsx+MXDKdjKxus0hjTeFiCqUfzVhdy+Qvzyd9T7Ns28cju/N/ph9lglcaYRidkglFVq2GOope/Xssdby+lpMxVuKSlCP/vp4O44CgbrNIY0zjZHUyMlZSV8+d3v2fqnMrBKtu3SOfRScM5um9mHCMzxpjYsgQTQwV7DvC7Fxfw9arKwSoPPaQ1T140gm7tbLBKY0zjljCPwUTkORHRMK/l0Timvny3cSc/mzK7SnI5dfAhTPvd0ZZcjDFJIRHvYGYDuUG2b4ryMTHz7uKNTH5jMUUlbnIwEbjppAFcMbav1bcYY5JGIiaYp1T1uXo4JurKypX7/7eCxz5b6dvWMiONhyYOZfyhNlilMSa5JGKCaZB2FZVw3SuL+GR55WCVvbNa8ORFI+jX0QarNMYkH0swUbBy2x4ufX4eedv2+rYd378DD58/jDbNbLBKY0xySsQEc4KIDMZNE7AFmAV8pKrlUT4mKj5dsZVrXl7Ibr/BKi8/vi+TTx5AaorVtxhjklciJpiLgmz7XkQmquq3UTymzt6Yv57JbyxGvcEqM9JSuPfswZw+tGus3tIYYxqMhGmmDCwCrgEOw92JdAFOAxYDg4CPRSTwm7s2x/iIyGUiMk9E5m3bti3igEf1bk+75ukAdGnTlDcuP8aSizHGeEQ1sSefFJF04HPgKOBRVb0qFseMHDlS582bF3F8X60s4OEZOTwyaRhZLTMiPt4YYxoyEZmvqiOD7UvER2RVqGqxiPwNeBv4SayOqa2j+2ZyVJ/21r/FGGMCJNIjsnAqeuRH8vypNsfUiiUXY4w5WENJMBWjQu4JW6ruxxhjjImShpJgzvWW38T4GGOMMVGSEAlGRIaKyGkikhqwPU1EbsC1FAN4oC7HGGOMqT+JUsnfC3gTKBSRH4D1QCvgCFzT43LgFlX9sI7HGGOMqSeJkmAWAw8BPwJ6AsMAxSWNZ3FNjedH4RhjjDH1JCESjKquAq6L9THGGGPqT8J3tKwvIrINWFNtweCygPwohmMSk13n5GDXOTI9VbVDsB2WYKJAROaF6slqGg+7zsnBrnP0JEQrMmOMMY2PJRhjjDExYQkmOv4V7wBMvbDrnBzsOkeJ1cEYY4yJCbuDMcYYExOWYIwxxsRE0icYEekmIo+IyAoR2S8iRSKSIyL/FJE+YY6bJCIzRWSniOzxZsa8UkTC/k5re5wJT0QGiMi1IvKCiCwXkXIRURE5uwbH1uu1tM9A7UR6jUWkiYiMF5H7RWSOiGwSkWIR2SAib4jI2Bq8p13julDVpH3hhpfZjhtiZh3wlvda723bDRwT5LhHvf37gfdwY6Lt8rb9B0gN8X61Os5eNbqWD3q/x8DX2dUcV6/X0j4D9XeNgR/7ldnk/b5fBb712/7nRPlsNMZX3AOI6w8PX3oX/V9AE7/tTYCnvX2LA445y+8Dm+23vRPwvbfv2iDvVavj7FXja/kb4F7cNA19gc+qSzD1fS3tM1C/1xgYB7wBjAmy7zyg1Dv+BLvGMbpm8Q4gbj84NKXyr5jOQfZ38dvf3G/7PG/bRUGOOd7vw5USsK9Wx9mr1te3JgmmXq+lfQbq/xpXc/xT3vFPx/uz0VhfcQ8gbj+4u0sp8S74IUH2H+Lt20Nlc+5u3rYDQLMQ5614vHaM37ZaHWevOl3fsF8+9X0t7TNQ/9e4Bsdf6R3/YTw/G435lVwVTn5UtQSY4a3+n4g0qdjn/fsv3urT6n0ycHU2AN+p6v4Qp/4moGxdjjOxU9/X0j4DiSfbW24K2G7XOEoSYrj+OLoCmA5cCpwiIvO87UcC7XDzzUz2K9/bW4YbdXltQNm6HGdip76vpX0GEoiIdAYu9lanBey2axwlSZ1gVDVPRI4BngdOwd3iVpgHfOHd6VRo6S33hjntHm/ZKgrHmdip72tpn4EEISJpwAtAG2CGqr4bUMSucZQk7SMyAC+5LAX6Aafj5oHoAJyBu4OZJiJ3+B/iLZXI1PY4Ezv1fS3tM5A4/gmMx3VNuCDIfrvGUZK0dzAi0hbX56UFrsItz2/32yLyHbAEuF1EXlbVHFy/GKj8SyWYin27/bbV9jgTO/V9Le0zkABE5CHgEmAzMF5VNwcpZtc4SpL5DuZU3N3KnIDkAoCq5gJzcUl4rLd5tbfsGea83QPK1uU4EzurvWV9XcvaHmeiRETuB64BtuGSS06Ioqu9pV3jOkrmBNPDW+4MU2aHt2zvLRd6y8NEpFmIY44MKFuX40zs1Pe1tM9AHInIvcANQAFwoqp+H6a4XeMoSeYEs9FbjvBvolzB2zbCW10FoKrrgAVAOnBOkGOOxzUU2Ax8VbG9tseZ2Knva2mfgfgRkbtxrUG345LL4nDl7RpHUbw74sTrBXTEtfZQYAqQ4bcvA3jc21cItPHbdzaVvXH7BZzvO0IPIVGr4+xV6+v7GdX35K/Xa2mfgbhc4zu9MtuBERGc265xFF5JPeGYiPwSN+ZYKu6OZj6uJcgIXE/+A8BEVX0r4LjHgN8BRcDHuBEBxgOtcQ0HzlbVsiDvV6vjTPVEZDjwmN+mQbimoDm4PxIAUNWjAo6r12tpn4Hai/Qai8jPgLe9zfNwX/DBLFfVu4O8n13juop3hov3CxiO6wezCveBKAJW4sYpGhTmuEnAbNwoqXtxyelKqhljqLbH2ava6ziW4CPtVnklwrW0z0D9XGNcR8pqywOf2TWOzSup72CMMcbETjJX8htjjIkhSzDGGGNiwhKMMcaYmLAEY4wxJiYswRhjjIkJSzDGGGNiwhKMMcaYmLAEk8REZLWIqPf6WzVlX/Qr+1mQ/SoiB3WqEpHPvH1jA7b/ydv+pzr+GA2aiFzs/R6ei3csNSEij4tImYgMCthe8VnqFZ/I4kNE3hGR3SJySLxjSUSWYEyFi0QkNdgOEWkNnFkfQYjI2FBJrCEKlXgbIhEZgpte/AUNPxpxMvkjbk6pv8Y7kERkCcaAG6epC3BiiP0TgWbAN2HOcaj3qqkpXvkpERzTGL2J+z3cGu9AauDvuO+MP8c7kEShqt/iruEvRWRwvONJNJZgDMBz3vLiEPsvBsqAqaFOoKrLVXV5Td9QVfO9Y/JrekxjpKo7vd/DpnjHEo6IHIb7A+QzVV0Z73gSzDO4QXKvjncgicYSjAE3c+f3wOneVNI+IjIAOBr4EDcEeVCRPgoKVgfjPRb71Fs93q/Op8ojMxHpKSK3isinIrJORA6ISKG3PinE+/kevYlIcxH5i4gsF5H9IrJIRI7z9i8LE3OWiBR5x7QPVc7/5/Nb9/9Z/LcHrYMJiLepiNwpIrnee+eJyG0VjzRFpLuIPC0iG7z4vhWRYHPNV5y7iYhcLiIzRWS7d0yOiPxDRDqEOOwKb/nvcD+3d/4TRWSGiOwUkX0iMscb2ThU+SwRucfveuzyjrlCRA6a1r26+rswv1PfdhHJFJGHRWSViBSLyFt+5U4SkfdFZKuIlHifreUi8oy4EZ0DTQe2AJMC//8ku4MunklazwH3Aufj5sKpcLG3fLYeYpiOG836ZNx/2Ol++/zvji7EzfOx0ts+GzeR0xhgrIiMUtVrQ7xHU9w8IocCXwCLgXRV/UJEFgNDRGScqn4S5Njf4OYKek5VC4Ps97cI92X8S2+92i/mENKBj4DDvLhzgONwP39XEbkP9/PvA2YCXYHRwFQRUVV90f9kXn3a+16ZnbhRfnfgRhW/HjhLRI5X1dUBcZzuLT+uJt5LcPUS3wD/BQYAo4C3RORcVX0jIJ5+wCe4qYQ3A+8CzYETgEeBM0XkNFU9UM37RiLLi68N7nc2DzfTJSJyMe6zXo77w2sN0NKL72LgB9ykYj6qWub9AXQe8GOgys+Y1OI9nLO94vfCzQuuwEigM1AKzPXbnwpswP3nS6dyMqXPgpwr6HD4VE4KNTZg+5+87X8K2D421Hv4lTkSOCzI9mxgrXf8qBDnVdx0tZ2CHH+Jt/+NIPtScFM6KHBkBL/jkNMEePsv9so8FybemVSd9G4IUIx7bPk98CCQ6rf/Su+43CDv94q373WgXcC1vifY7977vSqwtgafpQPAhIB9t3n7coIc97W37zWgqd/27sAKb9/favLZqcHvtGK74u7IWwU5Ns/bf0yQfd0IMYUHcK133GPR+v/ZGF72iMwAoKqbcXcMPxKRisr6k3CV/y+panHcggugqt+o6kGTR6lqDu4ve3DJMJQrVXVLkO0v4SauOl1EugTsOxXoBXyjquEaO0RbOXCZqu6s2KBuyt//4pJec+BmrTqB1RO4n6OviPSo2CiuafF5uL/KL1LV7X7nLMM1NFiCezx5hN/5hnrLkI8P/TyiqtMDtt2Lu1vqFxDPGNwfC7uBy1W1yC+edcB13uqVItK0Bu9dUyXAb1V1d5B9nYAdqvpl4A5VXa+hW89VbB8WpRgbBUswxt9z3vLigOVzJBivXuJ0ry7lCe+5+nNUJpb+IQ7dEuzLA0BV9+MmmksDLgvYXVEH8WgdQ4/UGlUN9sWe6y0/CUz+qlqKu9sC9wdChVO85Xvez1qFqpYDs7zVo/12dfSWBTWI970g5y3G3RkExnO8t3xXgzxyVNUPcPV+rXCzzEbLAj34EWCFr4G2IvK8iAwTEanhOSvi71Tn6BoRq4Mx/t7BfYlcKCJ/xz13/1ZV58c3rKpE5GjcI5VuYYq1DrF9TTWnfxS4EbhURP6iqqUi0hdXL1QAvBppvHW0PsT2PTXc7/+Xfx9veaWIXFnN+/pX9rfxlruqOQbcI8pgKo71j6ert1xFaHm46cu7hikTqXCfgStwdVQXeq+dIvI1rh5sqnenH0zFz2eV/H4swRgfVS0WkZdwzS2fxVVo10flfo2JSHNcv4NOwNO4Bgm5wG5VLReRk3DP10P95XnQX+7+VHWtiLyD61h6Bq7C9nfe+Z7xf4xTT8rruN9fRUfa+cDSasr6P4Lc4S1DJe3axlNxjcK1PqzpHYS/6p7MhPwMqOoycS0nTwbGAcfiGhycCPxJRM4K8ggQKn8324PsS1qWYEyg53AJ5jRcpf+LYUvXv+NwyWW+qv4myP5+UXiPR3AJ5goReQ/4Fe6L859ROHc8rfOWn6rq5AiO2+otM6McT8XdV58wZXp7yw1+2yoeCbYMcUzPugSlqiW4R33vAYhIO+D/4Srynyb43VTF72ZrkH1Jy+pgTBWqugD3HL4AeF1V6/s/TMWXR6g/fir6n6wLsT9oP5hIqOqnuL/wTwD+z3vPD1Q1L+yBwZUABOvPEQcfeMszIoynolnuoLClIve5t/yp9yVehYicjHs8tgd311WhItkMDHKMABOiGaTXGGIy7o+MLiH6ClX8bhYE2Ze0LMGYg6jqGFXNUtU6f1nXQsWXR78QX4IV/WHGiYjvC0ZEUkTkDtwjjWioGMLmZm/5WC3PU/HzRDKMTkx4fzy8hbvLe01EDqrDEpFDROQ6/9+9l1jXAt2CHVOHeGbi+qO0Ah4VkQy/OLriml8DTAl4NPkp7sv+FBEZ7XdMKnAX8KPaxCOuA+4NIRLIqbjvy11UPjL0V9Eo4tMg+5KWJRiTUFR1DV4/FWCJiEwVkadEZLK3fwGuM15rYJGIfCAir+A6IN6OaxIbDVOpfJ6eR9VOn5F401vOEJFXvJ/lqTpHV3u/xN05nAnkeD3mX/V+j0txj60e4OA7yIqe7j+OcjyTvPc8H8jzYnkX1wdmIDAD1+/FR1XX4ureUoBPvFED3sRdp6uAh2sZSzpwP7BJRBaKyGveNfuGyut4i/cIzcdLxmNxnV2r64iaVCzBmET0c1wrsfa4L55LcH9BVjgb+D2ucn8sMB5XKT2aysdAdaKq+4CK5syPe014a+OPwD9wj3l+jvtZLql7hLWjqrtwv6+LcCMZ9PXiGoGrc/sncHKQxgwVozv8kihS1Vxc35G/435Hp+Ou6Xe4ZHGKBu/Ffw1wCy6pjPZec3GdhhfWMpw9uAYdb+AGdz0Z+BmuZdhLwFGqGqwebgKuKfdLqhrs7iZpiWqjGEncmKgSkY64ep4yoFuwfhrJRkSm4zrfZqsNeOkjIv/BtTgcqqpL4h1PIrE7GGOC+yPukcm/Lbn43Iyr+7gj3oEkCm/EgzNwnxNLLgHsDsYYj4gcA/wa99hoLG54k0GqujGecSUSEXkM+C0wONhwPclGRN7G9Zfprwk+5UI8WIIxxuM3ku4+XHPTyao6J65BGdOAWYIxxhgTE1YHY4wxJiYswRhjjIkJSzDGGGNiwhKMMcaYmLAEY4wxJib+P61Z8sQyQUisAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams.update({'font.size': 22})\n",
    "plt.rcParams['lines.linewidth'] = 3\n",
    "\n",
    "# x axis values \n",
    "x = np.array1\n",
    "# corresponding y axis values \n",
    "y = np.array2\n",
    "  \n",
    "# plotting the points  \n",
    "plt.plot(x, y) \n",
    "  \n",
    "# naming the x axis \n",
    "plt.xlabel('Military time(hours)') \n",
    "\n",
    "# naming the y axis \n",
    "plt.ylabel('Temperature(F)') \n",
    "  \n",
    "# giving a title to my graph \n",
    "plt.title('Ta vs Time') \n",
    "  \n",
    "# function to show the plot \n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Answer__ : This graph does look correct as the temperature increases with the time. I am not sure of a better way to get Ta(t)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hours passed since death =  1.6498199340873931\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "math.log(10)\n",
    "\n",
    "Temp_t1 = 85\n",
    "Temp_a1 = 65\n",
    "Temp_a2 = 60\n",
    "Temp_a3 = 58\n",
    "K = 0.00458333333333\n",
    "Temp_0 = 98.6\n",
    "t_1 = 11\n",
    "\n",
    "''' Determine the time of death, T = 98.6, considering different Temperatures at different times of the day\n",
    "    \n",
    "    Arguments\n",
    "    ----------\n",
    "    Temp_t1 = Temperature at point 1\n",
    "    Temp_a1 = ambient temperature at 11am\n",
    "    Temp_a2 = ambient temperature at 10am\n",
    "    Temp_a3 = ambient temperature at 9am\n",
    "    K = empirical constant\n",
    "    Temp_0 = Temp at time of death\n",
    "   t_1 = Time first given\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "t1 = time in seconds since death calculated with ambient temperature 1\n",
    "t2 = time in seconds since death calculated with ambient temperature 2\n",
    "t3 = time in seconds since death calculated with ambient temperature 3\n",
    " t_hours = time in hours since death calculated using average of time in seconds for each ambient temperature\n",
    " t_death = time in hours (decimal) since death\n",
    "   \n",
    "    '''\n",
    "\n",
    "t1 = math.log((Temp_0 - Temp_ambient)/(Temp_t1-Temp_ambient))/K\n",
    "t2 = math.log((Temp_0 - Temp_a2)/(Temp_t1-Temp_a2))/K\n",
    "t3 = math.log((Temp_0 - Temp_a3)/(Temp_t1-Temp_a3))/K\n",
    "t_hours = ((t1/60)+(t2/60)+(t3/60))/3\n",
    "\n",
    "\n",
    "print('Hours passed since death = ', t_hours)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time of Death = 9.350180065912607\n"
     ]
    }
   ],
   "source": [
    "t_death = t_1 - t_hours\n",
    "print('Time of Death =',t_death)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Answer__: Time of Death = 9:35 am"
   ]
  },
  {
   "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
}