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 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": 451,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1306,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The K is 0.275\n"
     ]
    }
   ],
   "source": [
    "t_amb= 65 #ambient temperature\n",
    "t1= 85 #starting temperature\n",
    "t2= 74 #ending temperature\n",
    "dt=t1-t\n",
    "dt1=t1-t_amb #ending temerature- ambient temperature\n",
    "deltatime= 2 # 2 hours later\n",
    "K = dt/(2*dt1)\n",
    "print('The K is {:2,.3f}'.format(K))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Change your work from problem 1 to create a function that accepts the temperature at two times, ambient temperature, and the time elapsed to return $K$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1307,
   "metadata": {},
   "outputs": [],
   "source": [
    "def measure_K(Temp_t1,Temp_t2,Temp_ambient,delta_t):\n",
    "    ''' Determine the value of K based upon temperature of corpse \n",
    "    when discovered, Temp_t1\n",
    "    after time, delta_t, Temp_t2\n",
    "    with ambient temperature, Temp_ambient'''\n",
    "    \n",
    "    dtemp= Temp_t1- Temp_t2\n",
    "    dtemp1= Temp_t1- Temp_ambient\n",
    "    k= dtemp/(delta_t*dtemp1)\n",
    "    \n",
    "       \n",
    "    return k\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1308,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The K is 0.275\n"
     ]
    }
   ],
   "source": [
    "k=measure_K(85,74,65,2)\n",
    "print('The K is {:2,.3f}'.format(k))"
   ]
  },
  {
   "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": 1309,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_con(t):\n",
    "    timestep= np.linspace(0,2,t)\n",
    "    kconst=.275\n",
    "    t_ana=t_amb+(t1-t_amb)*np.exp(-kconst*timestep)\n",
    "    delta_time= np.diff(timestep)\n",
    "\n",
    "    tnum=np.empty(len(timestep))\n",
    "    for i in range(0,len(timestep)-1):\n",
    "        tnum[0]=85\n",
    "        tnum[i+1]= tnum[i] - (tnum[i]-t_amb)*kconst*delta_time[i]\n",
    "        \n",
    "    return t_ana, tnum, timestep\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1311,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 223 µs, sys: 0 ns, total: 223 µs\n",
      "Wall time: 223 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "o=20\n",
    "\n",
    "t_ana,tnum,timestep=time_con(o);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1312,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3gU5fbA8e9JgYQmUvRSBb2IiPRepEpTpIqIKB0EFNtPFGygVwUFr4iCCtiwIIIQQJEmYkPk0lEBEQVpCoKhSEg9vz9mE5eQsgnZTHZzPs+zT3annl2Gs7PvzHteUVWMMcbkHyFuB2CMMSZ3WeI3xph8xhK/McbkM5b4jTEmn7HEb4wx+YwlfmOMyWcs8RtjTD5jid/keSJy2uuRJCIxXq/7uh3fhRCR30WkudtxmPwlzO0AjMmMqhZJfi4ie4EhqrrKvYh8IyJhqpoQ6PswwcfO+E3AE5FQEXlMRH4RkT9F5D0RKe6Zd5WIJIjIYBE5KCLHRGSQiDQRke9FJFpE/uu1reEislpEXhORkyLyo4i08JpfQkRme87U94vIOBEJSbXuNBH5Cxjj2f8aETkuIkdF5G0RKepZfh5wCbDC8+vlbhHpKCI/p3p/Kb8KRGSiiLwvInNF5BRwS0bv35i0WOI3wWA00B5oDpQH4oEXvOaHAjWBy4GBwEvAA0BLz/SBItLIa/kWwFagJDARiBKRYp557wEnPNtqCHQDbk+17hagFPC8Z9qTwL+AGkBV4BEAVe0FHAHaq2oRVZ3q4/vtCbwNXAR85MP7N+YclvhNMLgDGKOqh1T1LPAE0FtExGuZJ1U1VlUXe17PVtVjqvobsBao47XsflWdrqrxqjobOAB0EJHLcBL7/ap6RlUPA1OBW7zW/UVVZ6pqoqrGqOpOVV2tqnGq+jswBecL50J8oapLVTVJVWN8fP/GpLA2fhPQPMmtArBURLwrDobgnLEDJKrqMa95McAfqV4X8Xp9INVu9gFlgcuACOCoV04NAbybZvaniq8s8CLQFCjqWf6wL+8tAyn78OH9/3mB+zJByM74TUBTp7zsQaCNqhb3ekSoanaTXvlUrysCh3AS7mngYq/9FFPVut4hpVp3EvA3cI2qFgOGAJLB8n8DhZJfiEg4UCLVMinr+On9myBnid8Eg1eBiSJSAUBELhGRGy9gexU8F2rDROQ2nMS/QlV/BdYBz4lIUREJEZEqmdyOWRTny+KkiFQE7k81/w+c6wXJdgAlRKStJ+k/Qeb/T3P6/ZsgZ4nfBIPngFXAas+dLmuBuhmvkqEvcdr8j+NciO2uqic88/oAxYGdnvlzgUsz2NbjOBddTwALcS7GensaeNpzd9FdnrP0e3AuIh8Afifz5pqcfv8myIkNxGLMP0RkOHCTql7ndizG+Iud8RtjTD5jid8YY/IZa+oxxph8xs74jTEmnwmIDlylSpXSSpUquR2GMcYElI0bN/6pqqVTTw+IxF+pUiU2bNjgdhjGGBNQRGRfWtOtqccYY/IZS/zGGJPPWOI3xph8JiDa+I0xuS8+Pp4DBw5w9uxZt0MxmYiIiKB8+fKEh4f7tLwlfmNMmg4cOEDRokWpVKkSVto/71JVjh07xoEDB6hcubJP6wRt4o/afJBJy3dxKDqGssUjGd2hKt3qlHM7LGMCxtmzZy3pBwARoWTJkhw9etTndYIy8UdtPsjYBduJiU8E4GB0DGMXbAew5G9MFljSDwxZ/XcKyou7k5bvIiY+kbryE0NDPwaUmPhEJi3f5XZoxhjjuqBM/IeiYwDoFvoNj4S/z/Phr1CQuJTpxpi8b//+/bRu3Zpq1apRvXp1XnzxxZR5x48fp127dlSpUoV27drx119/nbf+3r17iYyMpHbt2imP2bNnZ7jP8ePHM3ny5Bx9H9HR0UyfPj1Ht3mhgjLxly0eCcDjCQP4b/xN9Az9mg8KPEWNiyzxG+MvUZsP0mziaiqP+YRmE1cTtfngBW0vLCyM559/nh07drBu3TqmTZvGjz/+CMDEiRNp27Ytu3fvpm3btkycODHNbVxxxRVs2bIl5dGvX78Liim1hISETJexxJ9LRneoSmR4KCBMTezB8Lh7qSr7mSsPw8FNbodnTNBJvq52MDoG5Z/raheS/MuUKUPdus5AYkWLFqVatWocPOhsb9GiRfTv3x+A/v37ExUVlaVtFylSJOX5/PnzGTBgwHnL7Nmzh44dO1KvXj2uvfZadu7cCcCAAQO4//77ad26NQ899NA56/zwww80bNiQ2rVrU7NmTXbv3s2YMWPYs2cPtWvXZvTo0QBMmjSJBg0aULNmTcaNGwc4v1Cuuuoq+vfvT82aNbnppps4c+YMAGPGjOHqq6+mZs2aPPDAA1l6r2kJyou7yRdwk+/q2V6sJesat6DN5nvgzU7QdRrUuMnlKI0JHsnX1bwlX1fLiRsq9u7dy+bNm2nUqBEAf/zxB2XKlAGcL4gjR46kuV5ywk320ksvce211/q0z2HDhvHqq69SpUoVvvvuO0aOHMnq1asB+Omnn1i1ahWhoaHnrPPqq69yzz330LdvX+Li4khMTGTixIl8//33bNmyBYAVK1awe/du1q9fj6rSpUsXvvzySypWrMiuXbt4/fXXadasGYMGDWL69OkMGjSIhQsXsnPnTkSE6OjorH14aQjKxA9O8j/vgKv/OXzYDz4aDEd+hNaPQkhQ/ugxJleld/0sJ66rnT59mp49ezJlyhSKFSuWpXWTm3qys8+1a9fSq1evlGmxsbEpz3v16nVe0gdo0qQJTz/9NAcOHKBHjx5UqVLlvGVWrFjBihUrqFOnTsq+du/eTcWKFalQoQLNmjUD4LbbbmPq1Knce++9REREMGTIEG644QY6d+6c5feTWv7KekVKQ79FULcffPU8zL0NYk+5HZUxAS/5upqv030VHx9Pz5496du3Lz169EiZfumll3L48GEADh8+zCWXXJKl7Xrf/phWz+SkpCSKFy9+zvWBHTt2pMwvXLhwmtu99dZbWbx4MZGRkXTo0CHlF4I3VWXs2LEp2/35558ZPHjweXElvw4LC2P9+vX07NmTqKgoOnbsmKX3mpb8lfgBwgrAjVOh03Pw0zJ4vT38tdftqIwJaP9cV/tHZHgooztUzfY2VZXBgwdTrVo17r///nPmdenShbfffhuAt99+m65du2Zp25deeik7duwgKSmJhQsXnje/WLFiVK5cmXnz5qXEsnXr1ky3+8svv3D55Zdz991306VLF7Zt20bRokU5deqfE8wOHTrwxhtvcPr0aQAOHjyY0lT122+/8e233wIwZ84cmjdvzunTpzlx4gTXX389U6ZMydYvmNTyX+IHEIFGd8BtH8HJgzCjNez92u2ojAlY3eqUY0KPGpQrHokA5YpHMqFHjQtq3//mm2945513WL16dcrtmEuXLgWci50rV66kSpUqrFy5kjFjxqS5jeQ2/uTH1KlTAeeuoM6dO9OmTZuUawWpvffee7z++uvUqlWL6tWrs2jRokxjnjt3Ltdccw21a9dm586d9OvXj5IlS9KsWTOuueYaRo8eTfv27bn11ltp0qQJNWrU4Kabbkr5YqhWrRpvv/02NWvW5Pjx44wYMYJTp07RuXNnatasScuWLXnhhRey83GeIyDG3K1fv776bSCWY3tgzi1w/Be4fhLUH+Sf/RgTYHbs2EG1atXcDiPf2Lt3L507d+b777/P1vpp/XuJyEZVrZ962fx5xu+t5BUwZBVc3ho+vg8+eQAS492Oyhhj/MYSP0DERXDrXGg6Cv43E97tAWeOux2VMSYfqVSpUrbP9rMqaG/nzLKQUGj/FFxyNSy5h9PTWnBH/GjWnixl1T2NMUHFzvhTq30rXzR9m5jTJ3k19iFah2zKkV6IxhiTV1jiT8PD/4vgxtin2KuXMiv8ee4MjSImPsGqexpjgoIl/jQcio7hd0rSK24ci5KaMjr8Q14Jn8KJaGv3N8YEPkv8aUjubXiWgtwXP5L/xN9Gu5CNLIkcB3/+7HJ0xpiseuutt7jrrrsyXebQoUMpr4cMGZJSDTQr1qxZkyNlFfzJEn8azu2FKLyeeD1Dkh6hbPhpmNkGflruanzGmJyXOvHPmjWLq6++2sWI/McSfxrS6oXYrUcfCo78Ci6+DN7vDV9MgqQkt0M1Juh169aNevXqUb16dWbMmAE4ZZUfeeQRatWqRePGjfnjjz8AWLJkCY0aNaJOnTpcd911KdOTnTp1isqVKxMf7/TVOXnyJJUqVWLevHls2LCBvn37Urt2bWJiYmjVqhXJHUeXLVtG3bp1qVWrFm3btgVg/fr1NG3alDp16tC0aVN27Qqca4B2O2c60qzuCTBoOSy5Bz5/Cg5vge6vQsGiuR+gMbnp0zHw+/ac3ea/akCntAdQ8fbGG29QokQJYmJiaNCgAT179uTvv/+mcePGPP300zz44IPMnDmTRx99lObNm7Nu3TpEhFmzZvHcc8/x/PPPp2yraNGitGrVik8++YRu3brxwQcf0LNnT3r16sW0adOYPHky9euf29H16NGjDB06lC+//JLKlStz/Lhzre+qq67iyy+/JCwsjFWrVvHwww/z0Ucf5exn5CeW+LOqQCHoMQPK1oEVj8LMtnDL+1Dq325HZkxQmjp1akohtf3797N7924KFCiQ0o5er149Vq5cCcCBAwfo3bs3hw8fJi4ujsqVK5+3vSFDhvDcc8/RrVs33nzzTWbOnJnh/tetW0eLFi1StlWiRAkATpw4Qf/+/dm9ezcikvIrIhD4NfGLyH3AEECB7cBAVT3rmfcAMAkorap/+jOOHCcCTUbCv66BeQOcdv+eM+HKDm5HZox/+HBm7g9r1qxh1apVfPvttxQqVIhWrVpx9uxZwsPDU0oYh4aGpgyBOGrUKO6//366dOnCmjVrGD9+/HnbbNasGXv37uWLL74gMTGRa665JsMYVPW8cskAjz32GK1bt2bhwoXs3buXVq1aXfD7zS1+a+MXkXLA3UB9Vb0GCAVu8cyrALQDfvPX/nNF5RYwbI21+xvjJydOnODiiy+mUKFC7Ny5k3Xr1mW6fLlyThNtctnmtPTr148+ffowcODAlGmpyycna9KkCV988QW//vorQEpTj/e+3nrrrSy9L7f5++JuGBApImFAISD5kvkLwIM4vwQCW/GKMHgF1LzZaff/8HYb3MWYHNKxY0cSEhKoWbMmjz32GI0bN85w+fHjx9OrVy+uvfZaSpUqle5yffv25a+//qJPnz4p0wYMGMDw4cNTLu4mK126NDNmzKBHjx7UqlWL3r17A/Dggw8yduxYmjVrRmJi4nn7yMv8WpZZRO4BngZigBWq2ldEugBtVfUeEdmL84vgvKYeERkGDAOoWLFivX379vktzhyhCutecdr9S/7b2v1NwAvmsszz589n0aJFvPPOO26HkmOyUpbZb238InIx0BWoDEQD80SkH3An0D6z9VV1BjADnHr8/oozx6Rq949/tSUPy93MP3WNFXkzJg8ZNWoUn376acqgLvmRPy/uXgf8qqpHAURkATAQ54tgq+diSXlgk4g0VNXf/RhL7qncguXNPqDCiqE8ywT+FdqTl6O7MXaBcyucJX9j3PXSSy+5HYLr/NnG/xvQWEQKiZPl2wILVPUSVa2kqpWAA0DdoEn6Hk9+dZruseOJSmrG/4XPZ2b484THn7AibybgBMIIfSbr/05+S/yq+h0wH9iEcytnCJ6mm2B3KDqGWApwf/wIHosfQIuQbSwp8CgXndjhdmjG+CwiIoJjx45Z8s/jVJVjx44RERHh8zp+vY9fVccB4zKYX8mf+3dL2eKRHIyOAYR3EtvzfVJlphd4kYUFx8Pm4lCnr9shGpOp8uXLc+DAAY4ePep2KCYTERERlC9f3uflbbB1P4jafJCxC7YTE//PLV7lwk8TdenrlP7zO6g3ADo+C+G+f0MbY0xW2WDruSitIm+jezSn9Iil0Pw+2PgWvNkRogO7/5oxJjDZGb8bdn4CC4c74/z2nAX/vs7tiIwxQcjO+POSq25wSj0ULQvv3gRfPGelHowxucYSv1tKXgFDVnlKPTwNc3rDGRva0Rjjf5b43VSgEHR/DW54HvZ8DjNawuGtbkdljAlylvjdJgINhsCgZZCUCLPawabgqR9ijMl7bCCWvKJ8fbjjS/hoMCy+i71b1zDw917sPZFotX6MMTnKzvjzksKl4LYF7KoyjEr75vNizBjKyVEORscwdsF2ojYfdDtCY0wQsMSf14SEMmh/J4bE/R+V5A8+LvAwbUI2EROfaLV+jDE5whJ/HnQoOoZVSfXoHPc0B7Q0bxSYzJiw9zkSbQO8GGMunCX+PKhs8UgAftNL6Rk3nncSrmN42Md8FPkMnLDmHmPMhbHEnweN7lCVyPBQAGIpwGMJg/i/pLupFvIbvNocdq9yOUJjTCCzxJ8HpVXr59ruwwkf8SUULQPv9YTPnoTEBLdDNcYEIKvVE2jiY+DTB2HTbLisuVPrp1gZt6MyxuRBVqsnWIRHQpeXnB6/hzY5TT97Pnc7KmNMALHEH6hq3QJDP3fu/X+nO3w+wen5a4wxmbDEH8guuQqGroZafeCLic4XwOkjbkdljMnjLPEHugKFofsr0HUa7P/Oafr59Su3ozLG5GFWqydY1LkNytaBD/vD7C78WPUuhv3akoMnYq3WjzHmHHbGH0wurQ7D1rC/XCeu3jmVZ86M52JOWq0fY8w5LPEHm4JFuOXoYMbGD6ZRyE6WFhxL45AfrdaPMSaFJf4gdOjEWeYktqV73BP8rRG8H/4094XN54/o026HZozJAyzxB6HkWj8/aiVujHuaBUnXck/YAuZHToATB1yOzhjjNkv8Qci71s8ZInggfjgPJt1F9ZC9zl0/O5e6G6AxxlWW+INQWrV+mnYfSfjIr6F4RfigDyx9EOLPuh2qMcYFVqsnv0mIhVXjYd10+FcNuOlNKFXF7aiMMX5gtXqMI6wgdJwAfeY6tf1fawlb3ocAOAEwxuQMS/z5VdWOMOIbp9NX1AhYeAfE2ghfxuQHlvjzs2Jlof9iaPUwbJ8Hr7WAQ5vdjsoY42eW+PO7kFBo9RD0/9hp/5/VDr6dbk0/xgQxq9VjHJWawfCvYdGdsHwsv29ZxuDogfx4ooDV+jEmyNgZv/lHoRJwy/tsq/EwJX7/htfP3kfjkB+s1o8xQSbTxC8i9UXkPhGZJCJPisjNIlIiN4IzLhBhxO4GdI97kr81gvfCn2F02AfEx8darR9jgkS6iV9EBojIJmAsEAnsAo4AzYGVIvK2iFTMnTBNbjoUHcMPnnIPHya25M6wxXxUYDwFTvzidmjGmByQURt/YaCZqsakNVNEagNVgN/8EZhxT9nikRyMjuEMEYxJGMaapNpMDJ/J0oIPw6YCUOd2EHE7TGNMNqV7xq+q09JL+p75W1T1M/+EZdzkXesHYFlSQ7olTeJ0qdqweBR82A/OHHcxQmPMhcioqWeF1/Ox2dm459rADyLyvYjMEZEIz7WCnSKyTUQWikjx7Gzb+E9atX7u7dGK0iOXQbsnYden8Eoz+OULt0M1xmRDurV6RGSzqtbxPN+kqnWztGGRcsDXwNWqGiMiHwJLgUPAalVNEJFnAVT1oYy2ZbV68phDm+GjoXDsZ2g6Cto8BmEF3I7KGJNKdmr15EQPnjAgUkTCgELAIVVdoaoJnvnrgPI5sB+Tm8rWgTu+gHoDYO1UeP06OPqT21EZY3yUUeK/XEQWi8gSr+cpj8w2rKoHgck4F38PAydUdUWqxQYBn2Y3eOOiAoXhxinQ+z2I3u+Ue9jwpvX4NSYAZHRXT1ev55OzumERudizjcpANDBPRG5T1Xc98x8BEoD30ll/GDAMoGJFu2s0z6rWGcrVcwq9fXwv7F4JXV6CwiXdjswYkw6/1eMXkV5AR1Ud7HndD2isqiNFpD8wHGirqmcy25a18QeApCSnxv9nT0BkCej+ClzRxu2ojMnXstzGLyJLRORGEQlPY97lnl68gzLY529AYxEpJCICtAV2iEhH4CGgiy9J3wSIkBBoehcM+QwiLoJ3ujPnP7dTdUwUzSautnIPxuQhGTX1DAXuB6aIyHHgKBABVAL2AC+r6qL0VlbV70RkPrAJp0lnMzAD+AEoiNP7F2Cdqg6/8Ldi8oQyNVnS6H1OLRnLrSymVoGt3HtiJGMXxAFYoTdj8gCfmnpEpBJQBogBfsrtM3Vr6gkszSau5mB0DG1CNvFs+AyKcYZJCb1ZVqQHX4+9zu3wjMk30mvq8akss6ruBfbmcEwmSB2Kdjp8r06qS4fY55gQPotHw9/juphNEH2lM+C7McY1VpbZ5LiyxSNTnh+nGHfE38cD8XdQI2QvTG8Km9+z2z6NcZElfpPjUtf6AeGTkDZ8024xlKkJi0bC3Nvg7z9di9GY/MynxC8ikSJS1d/BmOCQVq2fCT1q0L5ZQ2eIx/ZPwe4VML0x7FrmdrjG5DuZXtwVkRtxOnAVUNXKnnLMT6pql9wIEOziblD64wdYcAf8sR3q9oMOz0DBom5HZUxQyU6tnmTjgYY4vW9R1S04t3Qak32XVoehn0Gze2HTO061z33fuh2VMfmCL4k/QVVP+D0Sk/+EFYR2T8BAT7mmNzvBynGQEOtuXMYEOV8S//cicisQKiJVROQlYK2f4zL5yWVNYMQ3UPd2+GYKzGzjNAUZY/zCl8Q/CqgOxALvAyeAe/0ZlMmHChZ1irv1+QBO/0Hiay2Z9tQorhizxEo+GJPDMuzAJSKhwBOqOhp4JHdCMvla1U4svXYh4Z/ex51Js6lfYD0PnLjDSj4Yk4MyPONX1USgXi7FYgwAT685ytDYe/m/uOFUk30sKzCGnknLmLxsh9uhGRMUfCnZsNkz8Mo84O/kiaq6wG9RmXzNKfkgfJTUgrWx1Xk2fAZPhb/JN2fWw19V4OLL3A7RmIDmSxt/CeAY0Aa40fPo7M+gTP7mXfLhMCXpFz+GMfFDqBO6B15paiN9GXOBMk38qjowjUdGdfiNuSBplXxYFNKOr9otgXJ1nZG+3u0BJw64FqMxgSzTph4ReZM0Bl635G/8JfkC7qTluzgUHUPZ4pGM7lCVDnXKQZNFsPENWPE4TG8CHSdA7b7gjO1gjPGBLyUbenq9jAC6A4dU9W5/BubNSjaY8xz/FRbdBfu+hirt4cYXoVhZt6MyJk9Jr2RDlsfcFZEQYJWq5tqAqpb4TZqSkmD9DFg1HsIKQKfnoGZvO/s3xuNCavWkVgWwkTSM+0JCoPFwp9dv6Wqw8A744FY49YfbkRmTp2Wa+EXklIicTH4AS3AGSzcmbyh5BQxcCu2fhj2rYXoj2D7f7vwxJh1ZbupxgzX1GJ8d/QmiRsDBDawJacwDZ/pTsPi/GN2hqvX6NflOtpt6ROQzX6YZkyeUvpKoem8yOelWmiRuYHnBB6l78jPGLthm9X6M8Ug38YtIhIiUAEqJyMUiUsLzqATY7RMmz5q04mdejuvMDXHPsF8v4aUCLzOVSby5zOr9GwMZ38d/B04VzrLARiD5VomTwDQ/x2VMtjklH+BnLU+PuCcYFPopD4R9SKOzo2DTGahzu935Y/K1dM/4VfVFVa0MPKCql6tqZc+jlqq+nIsxGpMl3iUfkghhVuINdIh7lp9DKsPiUfBON/hrr3sBGuMyX0o2vCQi14jIzSLSL/mRG8EZkx3nl3yAI2Hl+O3GuXDDf+HARpjeFL57zekLYEw+48vF3XHAS55Ha+A5INcGWjcmq7rVKceEHjUoVzwSAcoVj2RCjxp0q1sBGgyGO9fBZU3h0wed4R7/3O12yMbkKl9KNmwHagGbVbWWiFwKzFLVG3MjQLDbOY0fqMK2ufDpQxAfA63GQNO7IdSXSuXGBIYL6bkbo6pJQIKIFAOOAJfndIDG5CoRqHUL3LkeruwAnz0Bs9rA79vdjswYv/Ml8W8QkeLATJy7ezYB6/0alTG5peil0PsduHk2nDwMM1rB6qchIdbtyIzxmwybekREgPKqut/zuhJQTFW35Up0HtbUY3LFmeOw/GHYOgdKXwVdp0H5834lGxMwstXUo863QpTX6725nfSNyTWFSkD3V6HvfM6ciiZxVjtmPdqHthOWWq9fE1R8aepZJyIN/B6JMXlE1OmrafH3BOYktGZI2FLeirmbjxe8Y8nfBA1fEn9rnOS/R0S2ich2EbGzfhO0Ji3fxZ/xBXk0YTC9Yh8nlnBmhU4gcslwOH3U7fCMuWC+3LvWye9RGJOHJJd8APifXsX1cRMYGbaIkaGLYFoD6PAM1OpjZR9MwPKl5+4+oALQxvP8jC/rGROovEs+AMQRzpSEmxhQ8AUoVdUp+zy7Cxzb41KExlwYX3vuPgSM9UwKB971Z1DGuCmtkg+R4aHc3Ok6GPgpdH4BDm2BV5rCV/+FxHiXIjUme3w5c++OU6LhbwBVPQQU9WdQxrgp3ZIPdco5wz3WH+R0/KrS3un49VpLOGC3G5vA4Usbf5yqqogogIgU9nNMxriuW51yGY/YVayM0/Fr5yfwyQMw6zpodAe0eRQK2nmRydt8OeP/UEReA4qLyFBgFU4v3kyJyH0i8oOIfC8ic5IHdxGRlSKy2/P34gt5A8a46qob4M7voOFQp9rntMawa5nbURmTIV8u7k4G5gMfAVcCj6vqS5mtJyLlgLuB+qp6DRAK3AKMAT5T1SrAZ57XxgSuiGJw/SQYvNJ5Pqc3fNgfTv3udmTGpMnXUoTbgUhAPc+zsv1IEYkHCgGHcC4St/LMfxtYg3Px2JjAVqEBDPsC1k4lcc2znPlxBRPi+/Blket5oGM1G+zd5Bm+3NUzBKcoWw/gJpzOXIMyW09VDwKTgd+Aw8AJVV0BXKqqhz3LHAYuSWe/w0Rkg4hsOHrUOs2YABFWgKiifegc/yzfJ1bimfDXmRIzljcWfGI9f02e4Usb/2igjqoOUNX+QD18OEP3tN13BSrjjNtbWERu8zUwVZ2hqvVVtX7p0qV9Xc0Y101avosd8ZfSJ/4R/i9uOJXlMAtCHuLMx2Mh9rTb4RnjU+I/AJzyen0K2O/DetcBv6rqUVWNBxYATS/g9wsAABS8SURBVIE/RKQMgOfvkayFbEze9k/PX+GjpBa0jZ3Mh4ktuTVxEUxr5NwJZIyLfEn8B4HvRGS8pzPXOuBnEblfRO7PYL3fgMYiUshT3rktsANYDPT3LNMfWJT98I3Je1L3/I2mKA8nDGV4gWcg4iL44FaY0weif3MpQpPf+ZL49+CUZk4u3L8Ip82+KBl05FLV73DuBtqEc0E4BJgBTATaichuoJ3ntTFBI72evx07dYM7voB2/4Ff1jhn/19PsZ6/JtdlOuZuXmADsZhAE7X5IJOW7+JQdAxli0cyukPVc+/qid4Py8bAzo/hkqvhhv/CZU3cC9gEpfQGYvFlsPX6wCPAZXjd/qmqNXM6yPRY4jdBa+dS+PRBOLEf6twG1z0JhUu6HZUJEuklfl/u438P586e7UBSTgdmTL521fVweUv44ln4dprzRdD+P1DrVqcukDF+4MuRdVRVF6vqr6q6L/nh98iMyS8KFIZ2T8IdX0GpK2HRnfDWDXBkh9uRmSDlS1NPW6APTnmF2OTpqrrAv6H9w5p6TL6RlARb3iPu00cIiT/NzITrmV/4VkZ1rGU9f02WXUhTz0DgKpw6/MlNPYpzX74xJieFhBAlbXj27GTu03cZEbaErme/4dkF/UGH0K1uebcjNEHAl8RfS1Vr+D0SYwzg9Pw9HF+YB7mDuYmt+E/4W7wY+gL/+3g1VJgFpa90O0QT4Hxp418nIlf7PRJjDHDumL8btSo3xj3F4/H9qZq42xn1a+U4K/1gLogvib85sEVEdonINhHZLiLb/B2YMflV6p6/iYQyO7EDfSKmQc3e8M0UmNYQvl8AAdAPx+Q9viT+jkAVoD1wI9DZ89cY4wfp9fwd2rERdJvm1P0vVALmD4TZXeHoLpciNYHKl4FY9gEVgDae52d8Wc8Ykz0ZjvkLUKGhU/f/+slw2DPo+4rHIPZUhts1Jpkvt3OOA+oDVVX1ShEpC8xT1Wa5ESDY7ZzGpOvvP2HVONj8LhQtAx2ehuo9QMTtyEwekN7tnL6cuXcHugB/A6jqITIozmaMyUWFS0HXaTB4FRS5BOYPgtld4MhOtyMzeZgviT9OnZ8FCiAihf0bkjEmyyo0gKGfww3Pw+Ft8GozWPGoNf+YNPmS+D8UkdeA4iIyFFgFzPJvWMaYLAsJhQZDYNQm9lboRtLal/n9mRo88fTjRG064HZ0Jg/x5eLuZJy6+h8BVYHHVXWqvwMzxmRP1E9n6fRLL7rHPsHvejHj4l/kskXdWbN6uduhmTzCl8HWn1XVlao6WlUfUNWVIvJsbgRnjMm6Sct3EROfyFb9N93jnmR0/DDK8wctvuztFIA7baOd5ne+NPW0S2Nap5wOxBiTM7x7/iohzEtsRevY55mZcANsnQtT68I3UyEhzsUojZvSTfwiMkJEtgNVPT12kx+/AtZz15g8KnXPX4DTFGJ2kcEwch1c1hRWPgavNIGfVrgQoXFbRmf87+P00F3s+Zv8qKeqt+VCbMaYbEiv5+/oDlWh1L+h74fQd74z4/1e8F4v+HO3C5Eat9iYu8YEoUzH/AWnqWf9a/DFcxB/BhoNh5YPQsRF7gRtcly2x9zNCyzxG+NHp4/AZ086vX8Ll4K246B2Xxv6MQhcSM9dY0wwK3IJdH0Zhq6GiyvD4rtgZmv47Tu3IzN+YonfGOMoVxcGr4AeM4n56xC80Z6ox66n24QPidp80O3oTA6yxG+M+YcIUYnNaP73JF5O6EqnkPXMOXsnBxc+xpL/2QXgYGGJ3xhzjknLd3EsPpzJCb1pGzeJz5LqcmfIRzT5pB1sfs8ZEN4ENEv8xphzeHcAO6CXcFf83fSIHc+BpBKwaCTMbAV7v3YvQHPBLPEbY86RVgewTXold0U+Bz1mOmMAvHUDfNAXju1xIUJzoSzxG2POkV4HsAc6VoOaN8NdG6D1o7Dnc5jWCJY/AjHRLkVrssMSvzHmHJkO/VigELQcDaM2Qq3e8O00mFoH1s+ExARXYze+sQ5cxpgLc3irc9a/9ysoVRXaPwVV2tnwj3mAdeAyxvhHmVrQfwnc8j4kxTv1f97tAX/86HZkJh2W+I0xF04ErroBRn4HHSYQ99sGEl9pxnuP9qDzhAXWASyPscRvjMk5YQWIiuhKy7P/ZXZCO24OXcMHZ0ewf+E46wCWh1jiN8bkqEnLd3E4vhBPJPSnfdxzfJVUg1Eh85wOYBvfsgvAeYAlfmNMjvLuAParlmFE/H30jB3HvqRSsOQeeLUZ7FoGAXBjSbCyxG+MyVFpdQDbqFW5O/JZuPkdSIyHOb3h7Rvh4CYXIjSW+I0xOSrdEcA6XgVXd4E7v4PrJ8ORH53yz/MHw1973Qk2n/Jb4heRqiKyxetxUkTuFZHaIrLOM22DiDT0VwzGmNyXaQew0HBoOBTu3gLXPgA7P4GXGzh9Ac4cdzX2/CJXOnCJSChwEGgEzAReUNVPReR64EFVbZXR+taBy5ggduIgrHnGqfwZUcz5Mmg4DMIj3I4s4LndgastsEdV9wEKFPNMvwg4lEsxGGPyoovKQddpMOIbKN8AVj7m/ALY9qGVgPaT3DrjfwPYpKovi0g1YDkgOF88TT1fCKnXGQYMA6hYsWK9ffvOW8QYE4z2fE704rEUP7GD7UmVeL1gf1pd3/v8weJNplw74xeRAkAXYJ5n0gjgPlWtANwHvJ7Weqo6Q1Xrq2r90qVL+ztMY0weEXXySpoef5x740ZysZxmSvwTXBp1M2tWL3c7tKCRG009nXDO9v/wvO4PLPA8nwfYxV1jTIpJy3dxJl6JSmpOm9jneSL+dq5kH62+vBk+7Ad/Wg/gC5Ubib8PMMfr9SGgped5G8D+FY0xKbw7gMURzpuJnWgZ+wJTEnrAz585YwAsvhtO2uXB7PJr4heRQkA7/jnDBxgKPC8iW4Fn8LTjG2MMpN0B7DSFmFfkducW0IZDYcv7zhgAKx+HmL9ciDKw+TXxq+oZVS2pqie8pn2tqvVUtZaqNlLVjf6MwRgTWNLtANahKhQpDZ2ehVEb4Opu8M1UeLEWfPVfiDvjUsSBx3ruGmPylEw7gAFcXAl6vAbDv4YKjeGzJ5xfABvecEpCmAzZCFzGmMC3by2sGg/7v4MSV0CbR51fBCH5+9zW7Q5cxhjjP5c1hUHLoc8HEFoA5g906gDtWW1VQNNgid8YExxEoGonGPENG+tM4PDhg/BOdzb9pzlffvax29HlKZb4jTFBJWrr79y24XJanp3M+Ph+VEjcT4uv+vL79M5waIvb4eUJlviNMUFl0vJdxMQnEkc4byV2pEXsC0yMv4VCRzbDjJYw93Y4ssPtMF1lid8YE1S8O4ABxBDBq4ldaH52CrQcA3s+h+lN4KOhcGyPS1G6yxK/MSaopNUBDKBo8ZLQeizcuw2a3Q07ljhVQBffDScO5HKU7rLEb4wJKhl2AAMoVALaPQn3bIEGQ2DrHKcPwKcPwekjLkSc+yzxG2OCik8dwACK/guufw5GbYSavWH9TKcX8MpxQT8SmHXgMsYYcNr710yA7fOhYFFocic0HumMChagrAOXMcZkpOQV0HMWjFjLoRINYc0E/ppQjVeeupMl//vJ7ehylCV+Y4zxEnXoItoeGMqNsU+xJekKRiS8S7OP2/D9h09C3N9uh5cjLPEbY4yX5H4A2/VyBsY/RPfYJ9iWdDnX/Pg8TKkJ37wY8F8AlviNMcZL6n4Am7UKA+IfokfseChT0xkD4MVasPalgC0FbYnfGGO8pNcP4I+LasHtC51icJdWhxWPOl8A306D+Jg018mrLPEbY4yXTPsBVGwM/RbBwGVwSTVY/rDnC2B6wHwBWOI3xhgvPvcDuKwJ9F8MA5ZCqSth+Vh4sTasexXiz7oSu6/sPn5jjMkJv34FaybCvq+haBlofj/U7QfhEa6FlN59/Jb4jTEmJ/36JXw+AX5byxFKMj2+M18U6cQ9HWue/6vBz6wDlzHG5IbKLYiqPZMBiY+xN6kU48PfZu7Z4exa+AxL/rfb7egAS/zGGJPjJq34iTXx1bg5bhy9Yx9jd1I5Hgp5l+aftIGv/guxp1yNzxK/McbkMO++AN9pNfrGP0KP2PFsTawMnz0BL1wDa56FmGhX4rPEb4wxOSytvgCb9EoeKTwehq52Bodf8wxMqQGf/SfXq4Fa4jfGmByWYV+AcvWgzxy44yu4vBV8Ndn5Alj5OJw+mivxWeI3xpgc5lNfgDI1ofc7MHIdXNnRKQExpQYsGwsnD/s1Prud0xhj8oI/f4avnodtcyEkDOreDs3uheIVsr1Ju4/fGGMCwfFf4ev/krT5fRIVRsTdw45izRndoWqW+wHYffzGGBMISlQmqsIYrkt4kXcSrmN9UlUORscwdsF2ojYfzJFdWOI3xpg8ZtLyXfwSfzFPJvTjJEUAiIlPZNLyXTmyfUv8xhiTx6QeEyCz6Vllid8YY/KY9MYESG96VlniN8aYPCbTMQEuUFiObMUYY0yOSb57Z9LyXRyKjqFs8chs3dWTHkv8xhiTB3WrU85vZZytqccYY/IZS/zGGJPPWOI3xph8xhK/McbkM5b4jTEmnwmIIm0ichTYl83VSwF/5mA4OcXiyhqLK2ssrqzJq3HBhcV2maqWTj0xIBL/hRCRDWlVp3ObxZU1FlfWWFxZk1fjAv/EZk09xhiTz1jiN8aYfCY/JP4ZbgeQDosrayyurLG4siavxgV+iC3o2/iNMcacKz+c8RtjjPFiid8YY/KZgE78ItJRRHaJyM8iMiaN+SIiUz3zt4lIXV/X9XNcfT3xbBORtSJSy2veXhHZLiJbRCRHR5j3Ia5WInLCs+8tIvK4r+v6Oa7RXjF9LyKJIlLCM88vn5eIvCEiR0Tk+3Tmu3VsZRaXW8dWZnG5dWxlFleuH1uebVcQkc9FZIeI/CAi96SxjP+OMVUNyAcQCuwBLgcKAFuBq1Mtcz3wKSBAY+A7X9f1c1xNgYs9zzslx+V5vRco5dLn1Qr4ODvr+jOuVMvfCKzOhc+rBVAX+D6d+bl+bPkYV64fWz7GlevHli9xuXFsebZdBqjreV4U+Ck381cgn/E3BH5W1V9UNQ74AOiaapmuwGx1rAOKi0gZH9f1W1yqulZV//K8XAeUz6F9X1Bcflo3p7fdB5iTQ/tOl6p+CRzPYBE3jq1M43Lp2PLl80qPq59XKrlybAGo6mFV3eR5fgrYAaQuvu+3YyyQE385YL/X6wOc/8Glt4wv6/ozLm+Dcb7VkymwQkQ2isiwHIopK3E1EZGtIvKpiFTP4rr+jAsRKQR0BD7ymuyvzyszbhxbWZVbx5avcvvY8pmbx5aIVALqAN+lmuW3YyyQR+CSNKalvjc1vWV8WTe7fN62iLTG+c/Z3GtyM1U9JCKXACtFZKfnrCU34tqEU9vjtIhcD0QBVXxc159xJbsR+EZVvc/g/PV5ZcaNY8tnuXxs+cKNYysrXDm2RKQIzpfNvap6MvXsNFbJkWMskM/4DwAVvF6XBw75uIwv6/ozLkSkJjAL6Kqqx5Knq+ohz98jwEKcn3W5EpeqnlTV057nS4FwESnly7r+jMvLLaT6Ke7HzyszbhxbPnHh2MqUS8dWVuT6sSUi4ThJ/z1VXZDGIv47xvxx4SI3Hji/Vn4BKvPPBY7qqZa5gXMvjqz3dV0/x1UR+Blommp6YaCo1/O1QMdcjOtf/NOpryHwm+ezc/Xz8ix3EU5bbeHc+Lw826xE+hcrc/3Y8jGuXD+2fIwr148tX+Jy8dgSYDYwJYNl/HaMBWxTj6omiMhdwHKcq9xvqOoPIjLcM/9VYCnOlfGfgTPAwIzWzcW4HgdKAtNFBCBBnep7lwILPdPCgPdVdVkuxnUTMEJEEoAY4BZ1jjS3Py+A7sAKVf3ba3W/fV4iMgfnTpRSInIAGAeEe8WU68eWj3Hl+rHlY1y5fmz5GBfk8rHl0Qy4HdguIls80x7G+eL2+zFmJRuMMSafCeQ2fmOMMdlgid8YY/IZS/zGGJPPWOI3xph8xhK/McbkM5b4TUATkfEi8oCvy4jIABEpm8V9dBORq71ePyki12Uv4gz3IyKyWkSKiUil9CpKprHeXSIyMKfjMcHLEr/JbwYAWUr8QDcgJfGr6uOquiong/K4Htiq53fdT5eIhAFvAHf7IR4TpCzxm4AjIo94apGvAqp6Tb9CRJZ5imp9JSJXpVrvJqA+8J6nxnqkiDwuIv/z1GKfIZ4eO17rNAW6AJM861whIm95tpVcs/0ZEflWRDaISF0RWS4ie5I7oXmWG+3ZzzYReSKdt9YXWOT1OlREZopTr32FiER6trXGs88vgHtU9QywV0RyrQSDCWyW+E1AEZF6OHVV6gA9gAZes2cAo1S1HvAAMN17XVWdD2wA+qpqbVWNAV5W1Qaqeg0QCXROtc5aYDEw2rPOnjTC2q+qTYCvgLdweqk2Bp70xNwepyBZQ6A2UE9EWqSxnWbARq/XVYBpqlodiAZ6es0rrqotVfV5z+sNwLVpbNOY8wRsyQaTb10LLPSc5SIiiz1/i+AMQjLP66S9oA/bay0iDwKFgBLAD8CSLMa02PN3O1BEnfrqp0TkrIgUB9p7Hps9yxXBSeqpKz2W8Kyb7FdVTe7OvxGn5kyyuanWPQJchTE+sMRvAlFadUZCgGhVre3rRkQkAudXQX1V3S8i44GIbMQT6/mb5PU8+XUYTpGtCar6WibbSRCREFVNSrVdgEScXyTJvOvKgBN3TJaiNvmWNfWYQPMl0N3TPl8Up446nguiv4pIL0i5Q6ZWGuufwhnqDv5J8n96fjHclM4+vdfJjuXAIM8+EJFynhrvqe3CGU4vO64EfLoLyBhL/CagqDNc3VxgC04t86+8ZvcFBovIVpwmm7SGo3sLeNVTETEWmInTRBMF/C+d3X4AjBaRzSJyRTZiXgG8D3wrItuB+aT9RfIJTiXJ7GgG+ONOIxOErDqnMXmEOOOpzlbVdllcrw5wv6re7p/ITLCxM35j8ghVPQzMFJFiWVy1FPCYH0IyQcrO+I0xJp+xM35jjMlnLPEbY0w+Y4nfGGPyGUv8xhiTz1jiN8aYfOb/AShh4QpqNE9VAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(timestep,tnum,'o',label=str(o)+' Euler steps')\n",
    "plt.plot(timestep,t_ana,label='analytical')\n",
    "plt.title('Temperature')\n",
    "plt.xlabel('delta time (hr)')\n",
    "plt.ylabel('temperature (F)')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1314,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 318 µs, sys: 0 ns, total: 318 µs\n",
      "Wall time: 318 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "p=100\n",
    "\n",
    "t_ana,tnum,timestep=time_con(p);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1315,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZzN9f7A8df7bLaskZIkLcqSNfnRooS6N0uE5HZlaUdRQhIt4qLlKi2S1G0hZGsjCbfkVkqhUgnZKkuDNOZs798f58x0Os7MnDFzZubMeT8fj3mY8z3f7/d8zvGd93l/35/P9/MVVcUYY0zqcBR1A4wxxhQuC/zGGJNiLPAbY0yKscBvjDEpxgK/McakGAv8xhiTYizwG2NMirHAb4o9Efk94icoIukRj3sXdfvyQ0R+FpELirodJrW4iroBxuRGVY/L/F1EtgIDVHVZ0bUoPiLiUlV/sr+GKXks4zdJT0ScIjJaRH4Ukb0i8oqIVAo/d7aI+EWkv4jsFJF9ItJPRP5PRDaISJqIPBqxr5tFZLmIPCsiB0XkaxG5KOL5KiLyUjhT3y4iY0TEEbXtVBH5DRgRfv0VIrJfRPaIyIsiUj68/hzgBGBp+OxlsIhcLiI/RL2/rLMCEZkgIq+KyGwROQRck9P7NyYWC/ymJBgGtAcuAGoCPuCxiOedwLlAHaAv8ARwF3BxeHlfETk/Yv2LgC+B44EJwAIRqRB+7hXgQHhfLYAuwHVR264DqgKPhJc9AJwINATqAqMAVLU78CvQXlWPU9Upcb7fbsCLQEVgXhzv35i/sMBvSoKbgBGquktVjwD3Az1FRCLWeUBVM1R1UfjxS6q6T1V/AlYDTSLW3a6qT6mqT1VfAnYAHUTkVEKBfaiq/qGqu4EpwDUR2/6oqs+pakBV01X1W1VdrqpeVf0ZeJzQF05+rFTVt1U1qKrpcb5/Y7JYjd8ktXBwOwV4W0QiZxx0EMrYAQKqui/iuXTgl6jHx0U83hH1MtuAGsCpQGlgT0RMdQCRpZntUe2rAfwbaAWUD6+/O573loOs14jj/e/N52uZEsgyfpPUNDS97E7gUlWtFPFTWlWPNejVjHpcC9hFKOD+DlSOeJ0Kqto0sklR204CDgMNVLUCMACQHNY/DJTNfCAibqBK1DpZ2yTo/ZsSzgK/KQmeASaIyCkAInKCiHTMx/5OCXfUukTkH4QC/1JV3QKsASaKSHkRcYjImbkMxyxP6MvioIjUAoZGPf8Lof6CTN8AVUSkbTjo30/uf6cF/f5NCWeB35QEE4FlwPLwSJfVQNOcN8nRKkI1//2EOmKvUtUD4ed6AZWAb8PPzwaq57Cv+wh1uh4A5hPqjI00DhgXHl00MJyl306oE3kH8DO5l2sK+v2bEk7sRizG/ElEbgauVtXLirotxiSKZfzGGJNiLPAbY0yKsVKPMcakGMv4jTEmxSTFBVxVq1bV2rVrF3UzjDEmqaxdu3avqlaLXp4Ugb927dp89tlnRd0MY4xJKiKyLdZyK/UYY0yKscBvjDEpxgK/McakmKSo8RtjEsfn87Fjxw6OHDlS1E0xx6h06dLUrFkTt9sd1/oW+I1JcTt27KB8+fLUrl0bm8I/+agq+/btY8eOHZx22mlxbVNiA/+CL3YyackmdqWlU6NSGYZ1qEuXJicXdbOMKXaOHDliQT+JiQjHH388e/bsiXubEhn4F3yxk5FvrCfdFwBgZ1o6I99YD2DB35gYLOgnt7z+/5XIzt1JSzaR7gvQVL7jBuebgJLuCzBpyaaibpoxxhS5Ehn4d6WlA9DZ+RGj3K8y0TUNN352pqXTesJyFnyxs4hbaIyJ1K9fP0444QQaNGjwl+X79++nXbt2nHnmmbRr147ffvst67nx48dzxhlnULduXZYsWRJzv23atKFu3bo0btyYxo0bc/XVV+fYjq1btx7VhoIwc+ZMdu3aVeD7PVYlMvDXqFQGgLH+Pvzb35UerpW87HmYyhzMKvtY8Dfm2Cz4YietJyzntBFvFVgidf311/Puu+8etXzChAm0bduW77//nrZt2zJhwgQAvv76a2bNmsXGjRt59913ufXWWwkEAjH3/corr7Bu3TrWrVvH3Llz893WSNm9ZjQL/IVgWIe6lHE7URw85r+awd6BNJbNzPeM4XTZSbovwB2z11n2b0weZfaf7UxLR6HAEqmLLrqIKlWiby0MCxcupE+fPgD06dOHBQsWZC2/5pprKFWqFKeddhpnnHEGn3zySdyvd/311//lS+C44447ap1AIMCwYcM477zzOPfcc3n22WcBWLFiBZdccgnXXnstDRs2PGqb66+/ngYNGtCwYUMee+wx5s6dy2effUbv3r1p3Lgx6enprF27losvvphmzZrRoUMHdu/eDYTOUO644w5atWpFgwYNst7TypUrs85amjRpwqFDh+J+r7GUyM7dzA7cSUs2sTMtnUXBVmz3VmOa51Hme+5joG8wq4KNrNPXmDzK7D+LlNl/loi/oV9++YWTTjoJgJNOOolff/0VgJ07d9KyZcus9WrWrMnOnbG/fHr37k2ZMqEqQLt27Zg0aVJcr/38889TsWJFPv30UzIyMmjdujXt27cH4JNPPmHDhg1HDZ9ct24dO3fuZMOGDQCkpaVRqVIlnnzySSZPnkzz5s3x+XwMGjSIhQsXUq1aNWbPns2oUaOYMWMGAIcPH2b16tWsWrWKfv36sWHDBiZPnszUqVNp3bo1v//+O6VLl473I4ypRAZ+CAXyLk1OpvWE5exMS+cLPZPOGQ8y3fMIL7gn8qD/OmYGOiT0oDWmpMnsP4t3eaLEuo9IdiNbXnnlFZo3b57n11i6dClfffVV1pnBgQMH+P777/F4PLRo0SLmmPk6derw448/MmjQIP7+979nfVFE2rRpExs2bKBdu3ZA6Cwh88sNoFevXkDoLOjgwYOkpaXRunVrhg4dSu/evenatSs1a9bM8/uJVCJLPZEyyz4Au6jK1d4xvB9sylj3Szzseh6XdfoaE7fM/rN4l+dX9erVs8ogu3fv5oQTTgBCGf727duz1tuxYwc1atSIe78ul4tgMAiEvkS8Xu9R66gqTzzxRFb/wJYtW7ICebly5WLut3Llynz55Ze0adOGqVOnMmDAgJj7rV+/ftZ+169fz9KlS7Oej/4CExFGjBjB9OnTSU9Pp2XLlnz77bdxv9dYSnzg79LkZMZ3bcjJ4QPzD0pzk28IT/k7ca1rOf9xT6ASh6zT15g4RCZSmcq4nQzrUDchr9epUydefPFFAF588UU6d+6ctXzWrFlkZGSwZcsWvv/+e1q0aBH3fmvXrs3atWuBUH+Bz+c7ap0OHTrw9NNPZz333Xffcfjw4Rz3u3fvXoLBIN26dePBBx/k888/B6B8+fJZdfm6deuyZ88ePv74YyA0ZcbGjRuz9jF79mwAPvzwQypWrEjFihXZvHkzDRs2ZPjw4TRv3jzfgb/ElnoiZZZ9/rywCyb6r+G7YE3+5X6OhZ7R9PfdxQ++mlb2MSYHkf1nBXlVfK9evVixYgV79+6lZs2a3H///fTv358RI0bQo0cPnn/+eWrVqsWcOXMAqF+/Pj169KBevXq4XC6mTp2K0+mMue/IGn/VqlVZtmwZN9xwA507d6ZFixa0bds2ZgY/YMAAtm7dStOmTVFVqlWrltW5nJ2dO3fSt2/frLOJ8ePHA6HO5JtvvpkyZcrw8ccfM3fuXAYPHsyBAwfw+/3ccccd1K9fHwidNbRq1YqDBw9m1f0ff/xxPvjgA5xOJ/Xq1eOKK644hk/5T0lxz93mzZtrQd2IJXMqh53hmmQT+Z5pnkcpjZdBvoGsCDbhZJviwaSQb775hnPOOaeom2EIjerJ7ATOq1j/jyKyVlWP2lmJL/VE69LkZD4acWlW6ecLPZNOGQ+xTaszwz2ZG5xvsjPtDyv7GGNKrJQL/Jkia5W7OZ7u3vt4O9iCUe5Xmex+lqAvnTtf/7JAL1IxxpicrFix4piy/bxK2cAf3embTmkG+gbzqO9qrnau4jXPQ1TR/QV6kYoxxhQHKRv44eiyDwhTAl25yXsHZ8t2FpUaTUP5EcAmeTPGlBgpHfgzRQ9RWxJsQTfvWAI4mOO5n06OjwBsvL8xpkSwwM9fyz4COEX4Rk+lU8ZDfKmnM8UzlRGu13AQtLKPMSbpWeAPyyz7bJnwdx7p0Ygybif7qcA/vPfwsr8tN7sWM909mfL8YZO8GVNMzZw5k4EDB+a6TuRMmQMGDODrr7/O82utWLGCK6+8Ms/bFQcW+GOIPAPw4eJef39G+fpxoWM9CzyjqSOhg8ayf2OST3Tgnz59OvXq1SvCFhU+C/zZiO74fSVwGb2991BJfmeBZzRtHF8A1ulrTEHp0qULzZo1o379+kybNg0ITZc8atQoGjVqRMuWLfnll18AWLx4Meeffz5NmjThsssuy1qe6dChQ5x22mlZ0y0cPHiQ2rVrM2fOnKOmSG7Tpg2ZF4i+++67NG3alEaNGtG2bVsgNBNnq1ataNKkCa1atWLTpuT/e0+JKRvyY1iHuln37/1Ez6FTxkNM8zzKDPdkJvt78FSgU1anr13ta5LeOyPg5/UFu88TG8IVE3JdbcaMGVSpUoX09HTOO+88unXrxuHDh2nZsiXjxo3j7rvv5rnnnuPee+/lggsuYM2aNYgI06dPZ+LEiTzyyCNZ+ypfvjxt2rThrbfeokuXLsyaNYtu3brRvXt3pk6dGvPq2D179nDDDTewatUqTjvtNPbv3w/A2WefzapVq3C5XCxbtox77rmHefPmFexnVMgs8Ociem7/nVSjm3csE93TuNs9m3qOrQzz3cTONGxuf2PyYcqUKcyfPx+A7du3Z02BnFlHb9asGe+99x4Qmo2zZ8+e7N69G6/XG3OK5AEDBjBx4kS6dOnCCy+8wHPPPZfj669Zs4aLLrooa1+ZN4Y5cOAAffr04fvvv0dEYk7olmwSGvhFZAgwAFBgPdBXVY+En7sLmARUU9W9iWxHfh09yVspBvsGsiFYm+GuWZzu2c2NviFs91XnjtnrmLRkk2X/JjnFkZknwooVK1i2bBkff/wxZcuWpU2bNhw5cgS32501TbHT6cTv9wMwaNAghg4dSqdOnVixYgVjx449ap+tW7dm69atrFy5kkAgkOu9dFU15pz+o0eP5pJLLmH+/Pls3bqVNm3a5Pv9FrWE1fhF5GRgMNBcVRsATuCa8HOnAO2AnxL1+onw16t9hWmBjvT13U0N2csiz2haO0IZv3X6GpM3Bw4coHLlypQtW5Zvv/2WNWvW5Lr+ySeHEqvMaZtj+ec//0mvXr3o27dv1rLIKZIj/d///R8rV65ky5YtAFmlnsjXmjlzZp7eV3GV6M5dF1BGRFxAWSCzK/0x4G5CZwJJJbrTd1WwEZ28D/GrVuIl9wQGON8C1Dp9jcmDyy+/HL/fz7nnnsvo0aP/clvFWMaOHUv37t258MILqVq1arbr9e7dm99++y3rrlbw5xTJmZ27mapVq8a0adPo2rUrjRo1omfPngDcfffdjBw5ktatW8d9c/XiLqHTMovI7cA4IB1Yqqq9RaQT0FZVbxeRrYTOCHIs9RTktMwF5c+yT+hAKEc6k93PcIXzUxYGWjHcdwNHKGVTPJtiryRPyzx37lwWLlzIf/7zn6JuSsLlZVrmhNX4RaQy0Bk4DUgD5ojIP4HbgKNvRHn09jcCNwLUqlUrUc08ZtGdvocpwy2+O7g1uJC7XHM407OTG31D2JF2gnX6GlMEBg0axDvvvMPbb79d1E0pdhKW8YtId+ByVe0ffvxPoC9QH/gjvFpNQuWfFqr6c3b7Ko4Zf6To7L+NYx3/dj9JEAeDfIP4MNgQwLJ/UyyV5Iw/lRSXG7H8BLQUkbIS6ipvC7yhqieoam1VrQ3sAJrmFPSTQfQUzyuCjenkfYhftDIvuidwk3MxoNbpa4qtZLgTn8leXv//Ehb4VfV/wFzgc0JDOR3AtES9XlGL7vTdpifS1Xs/7wTPZ6T7NZ50T6EsR2yeH1PslC5dmn379lnwT1Kqyr59+yhdunTc26TcPXcTLbrsA8qNzjcZ7prFZq3Bjb6hbNWTACjjdjK+a0Mr/Zgi5fP52LFjB0eOHCnqpphjVLp0aWrWrInb7f7L8uxKPRb4EyD6hu4ArR3recL9BC4CDPHdyvvBZkCo7v/RiEuLqqnGmBLMbrZeiDLLPo/3bJx1g5ePgg3pmDGObVqd5z2PMNT1etb8/lb2McYUJgv8CRTd6buTalztHcsc/0UMdi1ghnsSFfndOn2NMYXKSj2F5K+1f6W3833GuF7kF63Czb4hbNTagA35NMYUHCv1FLHoeX5eCVxGT+99uCTAPM8YrnauBGyeH2NM4lngL0TRQz6/0DO5MmMcnwfPZLL7Wca5nseDz+b5McYklAX+IjCsQ92sTt99VOQ630ie9nekt+t9XvfcTw32WqevMSZhLPAXgehO3wBO/uXvxY3eIdSR3bxZ6h4ucKy3so8xJiGsc7eIRV/wdZrs5mn345wlO3jE352nAp1QHNbpa4zJM+vcLaais/8tehJXee9nYbAVw9yvM939CBVsyKcxpgBZxl+MtJ6wPOJqX+UfzmXc53qJX7QKt/huZ4PWAWzIpzEmPpbxJ4HITl8QXg60o4d3DA4JMs8zll7O97FZPo0x+WWBvxiJLvsArNMzuDJjHGuC9Rjvfp5H3M9Qmgwb8mmMOWZW6immojt9HQQZ5JzP7a43+E5rcqvvdn7UGghQw0o/xpgYrNSTZKKz/yAO/h3oxj99I6gmaSzy3MvfHWtQ7GpfY0zeWMafBKKz/xPZx1TPFJo5vmemvz0P+3vjJTQPt3X8GmMyWcafxCKzfwF+5nh6ekcz3X8F17uW8rrnfmrKHsCyf2NM7izjT0KRwz47OD5lkvsZgjgY6ruF5cGmWetZ9m9MarOMvwSJHPa5JHgeV3ofZodWY4ZnMiNcr+LCD1j2b4yJzQJ/Eoru+P1Jq9PNO5aX/W252fUmr3ke4kT2AdiwT2PMUSzwJ6no2ztm4OFef38GewdST7bxVql7uNjxJYDN9GmM+QsL/EkuOvtfFGxFR+84ftVKvOj5F8Ncs3ASsLKPMSaLde6WIJHDPkvhZazrRXq5PuCTYF0GewfyM8cD1ulrTKqwzt0UEJn9Z+BhpP8GbvfeSn3ZytulRtLGsQ6wTl9jUp1l/CVU5JDPOrKLqe5/c45jO0/7O/KIvzt+XIBl/8aUZJbxp5jIIZ8/ag26eB/kVf+l3OJazCzPQ9RgL2DZvzGpyAJ/CRXd6ZuBh3v8AxjkHUhd2c7bpUZymWMtYEM+jUk1FvhLsOghnwCLg6240juO7VqN6Z5HuM/1Eh58NuTTmBRigT8FRGf/2/REunnv5wV/B/q53mWuZyynys9W9jEmRVjnboqJnumzveNTJrqn4STIKF9/FgVbAdbpa0xJYJ27Bjg6+18aPI+/ZYznWz2FKZ4n+ZdrGmU4Ytm/MSWYZfwpLHLIp5MAd7jmcZtzIZu1BoN8g/hWawGW/RuTrCzjN0eJHPIZwMkj/h78wzeSinKYhZ7R/MP5HnZzd2NKnlwzfhFpDlwI1ADSgQ3AMlXdn/jmhVjGnzgLvtjJpCWbsjJ/gOM5wGT3M1zi/JJ3A+cx3HcDBzgOsOzfmGSS54xfRK4Xkc+BkUAZYBPwK3AB8J6IvCgitRLVYFM4Yg353EdF+vmG8aCvN5c6PuedUiNoId8AdsGXMSWBK4fnygGtVTU91pMi0hg4E/gpEQ0zhSszg8/M/hUHzwf+zifBc5jifoLXPA/xROAqnvBfRbovtJ5l/cYkJ+vcNUeJHvJZjnQecL9AN+eHfBo8izu8t7GTalb2MaaYO5ZSz9KI30cmqmGm+Ike8nmYMtzpu5XbvbdytmznnVIj+ZtjjZV9jElSOY3qqRbxe/dj2bmIDBGRjSKyQUReE5HSIjJJRL4Vka9EZL6IVDqWfZvEilX7Xxi8gL95H+ZHPYmnPFP4l2sa+A5zx+x1Nt2DMUkkp8CfrxqQiJwMDAaaq2oDwAlcA7wHNFDVc4HvCHUem2IqOvvfrtW52juGJ/2d6e5cyZueUdSXLZb9G5NEcgr8dURkkYgsjvg96yfO/buAMiLiAsoCu1R1qar6w8+vAWoee/NNYcjM/jODvx8Xk/09udY3irKSwXzPfdzoXMwRn8+yf2OSQLaduyJycU4bqurKXHcucjswjtD4/6Wq2jvq+cXAbFV9Oca2NwI3AtSqVavZtm3bcns5k2DRnb4AlTjEePd0rnB+yoeB+tzpu4VfqEIZt5PxXRtax68xRSi7zt2EjeoRkcrAPKAnkAbMAeZmBnkRGQU0B7pqLo2wUT3FR6wLvkDp6VzBGNdLZOBmhO8GlgTPA+yCL2OK0rGM6lksIh1FxB3juToi8oCI9MvhNS8DtqjqHlX1AW8ArcLb9wGuBHrnFvRN8RKr0xeE2YFLsub5f9bzGA+7nrPJ3owppnIq9ZwIDAW6AfuBPUBpoDawGXhSVRdmu2OR84EZwHmESj0zgc+A74FHgYtVdU88jbSMv3iKlf278TPENZebnYvZoidyu+82NmgdnCIEValhZwDGFJp8lXpEpDZwEqEA/p2q/hHni95PqNTjB74ABgAbgVLAvvBqa1T15pz2Y4G/eItV+2/p+JpH3U9RjQM85r+aZwIdCYZPMK3+b0zhKPQaf0GywF/8xcr+K/A7D7tncKVzDf8Lns2dvlvYoaHLQ06uVIaPRlxaVM01JiXYtMwmoWLV/g9yHAN9gxjivYV6so23PSPo4viQzKmebdinMUXDAr8pUJEXfAngFAfzgxdyhXc8m/QUHvc8xRPuJ6jA79bxa0wRibfGXwaopaqbEt+ko1mpJ3lF1v8dBLnZuYghrnnspSJ3+m5mdbABYMM+jUmEYy71iEhHYB3wbvhx4zxcuWtSXOQZQBAHTwW60NV7P39oKV71PMy9rv9QCq9l/8YUonjuwLUWuBRYoapNwsu+Cs+1Uygs4y8ZIu/xW5oMRrpepY/rPb4NnsJQ3y18rbUBy/6NKSj56dz1q+qBBLTJpJjIe/weoRRj/H3p4x1OFTnEAs9obnYuwkHQsn9jEiyewL9BRK4FnCJypog8AaxOcLtMCRQ90yfAymAj2mf8i2XBZoxwz2KW50FOkV9I9wVswjdjEiSewD8IqA9kAK8CB4A7EtkoU3LFGvaZRnlu9d3OEO8tnC0/8Y5nJD2cH5A57NOyf2MKVo41fhFxAhNUdVjhNeloVuMvmWJd9FWDvUx2P0Mr59e8F2jKSN8N7KUiYLV/Y/LqmGr8qhoAmiWsVSalxcr+d1GV3r57eMB3HRc51rOk1N10cHwCYNm/MQUknlE9jwBnEppW+XDmclV9I7FN+5Nl/CVfrOz/DNnBo+6nOdexhXmBC7jf14eDlAMs+zcmHsc8V4+IvBBjsapqTlMyFygL/KkjesI3F34GuRZwm3MBv1KJu3038WGwIWCTvRmTG5ukzSSNWNn/ubKZR91Pc4ZjFy/52zHe34t0SgOW/RuTnfxm/EetZBm/SbTo7L8UXu52zaa/6x22Bqtzp+9m1mpdwLJ/Y2LJzwVcbwJvhX/eByoAvxds84w5WvS4/ww8POi/jmu89+IkyBzPA4xwvUopvDbu35g8yHOpR0QcwDJVLbTJ1C3jN9HZfznSGeV6hWtdy/kueDJ3+m5hvdYBLPs3JlNBzsd/JlAr/00yJn7R2f9hynCPfwD/9A6nvKQz33MfQ12v48Zv2b8xuYinxn+Iv9b4fwZGquq8RDYskmX8JlJ09l+Bw4xxv0Q353/5JliLO303Z034Ztm/SWXHnPGranlVrRDxc1ZhBn1jokVn/wcpx52+WxjgvZPj5SALPaMZ7HwDVzj7n7SkSG4jYUyxFc98/O/Hs8yYwhTrqt9lwWa0z/gXbwfPZ6h7Lgs891FXfrLbPBoTJdvALyKlRaQKUFVEKotIlfBPbaBGYTXQmJxEZ/9plOd230Bu8g6huuxnsWcUA53z+SXtkE33YExYtjV+Ebmd0CycNYCdgISfOgg8p6pPFkoLsRq/iU907b8yB3nAPZOOzjWsD9bmLt/NbNJadsGXSRn5uYBrkKo+kbCWxcECv4lXrKt+L3d8woPuGVTkMFP8XXkm0BG3u5R1+poSL19TNohIA6AehK+RB1T1pQJtYQ4s8Ju8irzNI4Sy//vdL9LJ+TEbgrUZ5ruJb/RUy/5NiZafm62PAZ4I/1wCTAQ6FXgLjSlAkbd5BPiNCgz2DQrX/n9jkedehrjm8qvV/k0KiqfUsx5oBHyhqo1EpDowXVU7FkYDwTJ+c2xilX0AKnGIMe6XuMr5Ed8ETwll/5xOUJUadgZgSpD8XLmbrqpBwC8iFYBfgToF3UBjClqsIZ8QGvkzxHcb/b13Ull+Z4HnPoY6X8OD1272YlJCPIH/MxGpBDwHrAU+Bz5JaKuMKUCxbvIO8H6wGe0zJjIvcBG3uRbxlucemsp3NuWDKfFyu+euADVVdXv4cW2ggqp+VSitC7NSjyko0UM+M13o+Irx7unUYB8vBC5nsr876ZS2KR9MUsvPcM61qlqk9921wG8KUmbtf1daOg4RAuG/gXKkc7drFn1c7/FTsBoj/DewOtgAsJu9mOSUn8A/FZipqp8mqnG5scBvEiXWGUAL+YYJ7ueo4/iZ1/yXMN5/LQcpZ9m/STr5CfxfA3WBrYRuti6E7rl7bgLaGZMFfpNIsUb/lMLLENc8bnC+yV4qcq+vH+8FQ38/lv2bZJGfwH9qrOWquq2A2pYrC/ymMMTK/hvKj0x0T+Mcx0+8GWjJWF8f9lLRsn+TFPIzLfM24BTg0vDvf8SznTHJJtbon/Vah47eh5jk60E7x2csK3UX3RyrSPf5beSPSVrxZPxjgOZAXVU9S0RqAHNUtXVhNBAs4zeFL1b2f7rsZLx7Oi0cm1gVaMg9/v7s0BMs+zfFVn4u4LqK0BQNhwFUdRdQvmCbZ0zxEiv736wn09M7mnt9fWni+IGlnuH0d76F1+e17N8klXgCv1dDpwUKICLlEtskY4qHWDgG75YAABccSURBVFf+Kg5eDrSjfcZEVgfrMdr9CvM991FPttpVvyZpxBP4XxeRZ4FKInIDsIzQVbzGpIRY2f9ujmeA7y5u8w7mJNnHIs+9DHe9RtCXbtm/KfbinZa5HdA+/HCpqr4X185FhgADCJ0trAf6AmWB2UBtQkNEe6jqbzntx2r8priIVfuvyO/c43qVnq4VbA1W5x5/f1YHG1jt3xS5/M7HfyLQglAA/1RVf45jm5OBD4F6qpouIq8DbxOa13+/qk4QkRFAZVUdntO+LPCb4iS7WT//z7GRh13TOc3xC3MDFzHOdy2/UcHG/Zsik59x/AOA+4DlhC7euhh4QFVn5LLdycAaQlM6HwQWAFMIzevfRlV3i8hJwApVrZvTvizwm+IoVvZfCi+DXPO5yfkmBynLg77rWBBsjSAodvGXKVz5CfybgFaqui/8+HhgdW7BOrzu7cA4IJ1Qiai3iKSpaqWIdX5T1coxtr0RuBGgVq1azbZtK7TrxYyJW3bZf135ifHu6TR1/MCqQEPu9ffjJ60OYCUgU2jyE/jfB65QVW/4sQd4W1Uvy2W7ysA8oCeQBswB5gJPxhP4I1nGb4q7WNm/gyDXOt/nbtcsPPiZ4u/KtMDf8eMCLPs3iZefcfw7gf+JyNjwxVxrgB9EZKiIDM1hu8uALaq6R1V9wBtAK+CXcImH8L+/5vXNGFPcxBr5EwwP/bwsYzLLg0242z2bNz2jaCrfAdjwT1Nk4gn8mwnV5zNPDRYCuwldxJXThVw/AS1FpGx4Xv+2wDfAIqBPeJ0+4f0Zk/Syu+PXr1TmVt8d9PfeyXGSzhulxjLO9TwV+N1u+mKKRFyjeo555yL3Eyr1+IEvCA3tPA54HahF6Muhu6ruz2k/VuoxySay9i/8mTWV5QhDXHPp53yH/ZTnQd91LAq2AsRq/6bA5afG3xwYBZwK4eIkYNMyGxOfWB3A9WUr49zTaez4kVWBhoz292WbnghY7d8UnPyO6hlG6AKsYOZym5bZmLyJ7gB2EKS3cxnDXLMphZ8n/Z15NtARL27L/k2ByE/g/1BVL0hYy+Jggd+UFLGy/xP4jdHu/9DRuYbNwZO419+Pj4P1Acv+Tf7kJ/C3BXoB7wMZmctV9Y2CbmR2LPCbkibW8M+LHV/ygOsFTnX8yhuBC3jY19tu+mLyJT+B/2XgbGAjf5Z6VFX7FXgrs2GB35RE2d3y8TbXAm52LiadUkzy9+TVQFuCOCz7N3mWn8C/XlUbJqxlcbDAb0qy7G768qDrBVo5v2Zd8HRG+fqxUU+z7N/kSX4C/3PAY6r6daIalxsL/Kakiz31g9LZ8RH3ul+mCod4KdCeR/3d+UPKEVSlhp0BmFzkJ/B/A5wObCFU4xdCpR4bzmlMAYuV/VfgMHe6Xuc65zL2UYGHfL1ZGGyNjf03uclP4D811nIbzmlMYmQ38VsD+ZGH3C/Q2LGZjwP1GO2/nh+0JmCjf0xs+Z2P/wLgTFV9QUSqAcep6pYEtDMmC/wmFcXK/oUgvZwfcLdrFuU4wvOBK5ji78oflLbs3xwlPxn/GKA5UFdVzxKRGsAcVW2dmKYezQK/SVWZ2f+utHQcIgTCf69VOMhw1yx6ulawW6vwoO8fvB08HxDL/k2W/AT+dUAT4HNVbRJe9pXV+I0pXLHOAJrKdzzofoH6jm2sCjRkrL8PP2qNrPmB7EsgteVnWmavhr4dNLyjcgXdOGNM7mJN/fy5nkUn70Pc5+tDY8dm3vUMZ7jrNcpwBLCpn01s8WT8dwFnAu2A8UA/4DVVnZL45oVYxm/MX8XK/qtygOGu1+juWsUurcI43z94K1z+Acv+U1F+O3fbAe0JHUFLVPW9gm9i9izwG3O07Eb/RJZ/PgrUZ4y/T9boH+sATi35qfH/S1WH57YskSzwG5O9nG77OMw1m7JkMDPQgX/7u/I7ZQHL/lNFfmr87WIsuyL/TTLGFITo2r/w520fL8l4lLmBi+jvfIflpe7iKsd/AbXaf4rLNuMXkVuAW4E6hG6/mKk88JGq/iPxzQuxjN+Y+MUqAZ0rm3nAPZPGjs18FjyLMb4+bNTTAMv+S7I8l3pEpCJQmVCH7oiIpw7ldqvEgmaB35i8iy4BCUGudq5iuGsWVTjErMAlTPL34DcqWO2/hMpX525Rs8BvzLGJlf1X4DC3u96gj3MJhynNI/7uvBK4jABOy/5LGAv8xqSwWB3AZ8gOxrpe5ALnRr4NnsL9/n/ycbC+XfxVguSnc9cYk+RiXfz1g9bkH757uMk7hHIc4TXPOKa6H6cGewC7+Ksks4zfmBQTK/svhZcbnG9xm2shgvJsoCNP+ztyhFKAZf/JyjJ+YwwQO/vPwMOTgau4NOMRlgabc7vrDZaXupOOjtXY8M+SxzJ+Y1JYrOwf4Dz5ljHul2jg2MqnwbO43/dPNmgdwLL/ZGIZvzHmKLEu/gL4VM+mk/chhvtuoLb8zCLPaCa6nqUaaZb9lwCW8RtjssQa/lmePxjomk9f57t4cTPV34UZgcvJwGPZfzFnwzmNMXGLVQKqLbsZ5XqVds61/BSsxjh/b5YEz6OM22UXfxVTFviNMXmS3eyfFzjWc6/rZc52bGdN8Bwe9F3HRq1t2X8xZIHfGHNMYmX/TgL0ci5niGsulfmd1wMX84i/B2mOyhxX2kXaHz5q2BdBkbPOXWPMMYk1/DOAM2v2z+mBv9HV+V8+KDWUG2Q+f/xxGMUuACvOLOM3xsQtu+GftWU3I12v0cH5GTu0KhN917Ao+H/Y3b+KlmX8xph8i5X9A2zVk7jJN5Re3lEc0HJM8TzJG54xNJHvAcv+ixvL+I0xxyS77N9BkK7O/zLMNZvqksbiQEv+5b+GHXoCYNl/YbLOXWNMgcsc+bMrLZ2KZdwc9vrxBUIxpSxHuMn1Jjc638RBkBcClzPV34VDlLUZQAuJBX5jTMLFGgJanf3c5Xqdbs7/kkY5Hvd349VAW/y4ALsBfCJZ4DfGFJpYZaD6spVRrpdp5fyazcGTmODvxXvBZlgHcOJY564xptDE6gTeqLW51jeKft67UITnPI8y2/Mg50rolt7WAVx4Epbxi0hdYHbEojrAfcAK4BmgNOAHblXVT3Lal2X8xiSv7C4Au8b5AUNcc6kqB1kYaMUkf092aDXAsv+CUqSlHhFxAjuB84HngMdU9R0R+Rtwt6q2yWl7C/zGJLfI2n9mxy7AcfzBTa43GeB8GwdBXgx0YKq/Mwc4zmr/BaCoA397YIyqthaRJcAMVZ0tIr2Ajqp6bU7bW+A3puSI1QF8IvsY6prL1c5VHKQsT/q78J9AO5sBNJ+KOvDPAD5X1SdF5BxgCaEeHQfQSlW3xdjmRuBGgFq1ajXbtu2oVYwxSSxWCehs+YkRrtdo4/ySHVqVSb4eLAq2Ahw2/PMYFFngFxEPsAuor6q/iMgUYKWqzhORHsCNqnpZTvuwjN+Ykim7GUBbO9Yz0vUaDRxb2RCszcP+a1kdbADY8M+8KMrA3xm4TVXbhx8fACqpqoqIAAdUtUJO+7DAb0zJFiv7F4J0cqxmmPt1aspeVgbOZYK/F9/oqYBl//EoyuGcvYDXIh7vAi4O/34p8H0htMEYU4zFGv6pOFgYvIC2GZN50NebRo7NvOW5h0fdT1FT9tjwz3xIaMYvImWB7UAdVT0QXnYB8G/ABRwhNJxzbU77sYzfmNSR3RxAFTjMLa5F9HW+i6D8J9COqf7O/EYFy/6zYVfuGmOSRnbDPyE0AugO1zy6O1dymNI84+/IC4HLOUJp6wCOYoHfGJOUsusAPkN2MNw1m3bOtfyqlfi3vyuzA23w47IO4DAL/MaYpJZdCaiZbGK4exYtHJvYEqzOI/4evBU8H8WR8tm/zdVjjElq2d0EZq3WpYf3Pvp77yQDD096nmCR514udHzFzrQ/rAM4Bsv4jTFJJ6ebwHRxfMhQ91xqyl5WB+ox0X8N6/SMlMz+rdRjjClRcuoA9uDjWuf7DHQtoKocZEmgOZP9PfhBa6ZUB7AFfmNMiZVdB3A50unnfIcbXW9RliMsCF7AY/6r2aHVUqID2AK/MabEy64EVJmD3OJaTB/nUoQgrwbaMtXfhf1SmaAqNUroGYAFfmNMSsgu+4fQNQCDXW/Qw7kSL25mBjrwrP9KDnBcibwPsAV+Y0xKyS77BzhVfmaIay6dHB/zO2WY5v87LwQu5zChEUMlpQxkgd8Yk3Jy6gAGqCs/cadrDu2da9mrFXja34mXA5eRgQdI/uzfAr8xJqVlfgnsSkvHIUIgIvY1lh+4yzWbC5wb2a1VeNLfhdcDbfAl+VXAFviNMSYsuzJQS8fX3OV6neaO7/gpWI1/+7uxINiaAM6kzP4t8BtjTITsy0BKG8eXDHXN4VzHFjYHT+JxfzfeDLYk2e4EZoHfGGOyEXskkNLe8RlDXXM527Gdb4On8Ji/G0uC5wGSFCUgC/zGGJOL7O4EdqVjDXe45nG6YzcbgrV51H81y4NNACnW2b8FfmOMiUN21wE4CdDZ8RG3u97gVMevrAvW4TF/d1YGz0WQYlkCssBvjDF5kF0HsAs/XZ3/ZbBrPjVlL58Fz+Ixfzc+CjaguJWALPAbY0we5XQdgBs/3Z0rGeiaTw3Zz/+CZ/O4vxsfB+sDxSP7t8BvjDH5kF0JyIOPns4PuM21kBPlN9YEz+FxfzfWBOsV+TQQFviNMaYAZFcCKoWXa5wfcKtrIdUljY8D9Xjc343/6TlA0UwDYYHfGGMKSE4loFJ46eVczq2uRZwQ4wugMLN/C/zGGJMA2ZWASuGlt/N9bnYtPuoLoLBKQBb4jTEmgXIqAfVyLucW1yKqS1pEJ3A9Ej0KyAK/McYkWG4loGucH3CLaxEnym98EqzLFH9XPgwPA01E9m+B3xhjClFOJaCezg+4xbWYk2Q/nwfPYIq/KyuCjQr8QjAL/MYYUwSyKwF58NHduZJbXIuoKXv5MliHKf6reD/YlIIqAVngN8aYIpLbhWDdnKu41bmQWo49bAyeyhP+q1gSbI5DnPm6J7AFfmOMKQayKwG58NPZsZrbXAuo4/iZTcGaTPV34c1gS4I4jukMILvA78j/2zDGGBOvLk1O5qMRl/J4z8aUcTuzlvtxMS94EZd5JzPYOxBBmeJ5kr85/gdAui/ApCWbCqQNrgLZizHGmDzJzNyjS0BBHCwKtmKxtyVtHV+wItgoa5tdUWcJx8oCvzHGFJEuTU7O+gI4+p7ADpYFm/1l/RqVyhTI61rgN8aYYiD6SyB6JFAZt5NhHeoWyGtZ4DfGmGImsgy0Ky39mEf1ZMcCvzHGFEORZwAFzUb1GGNMirHAb4wxKcYCvzHGpBgL/MYYk2Is8BtjTIpJirl6RGQPsO0YN68K7C3A5hQUa1feWLvyxtqVN8W1XZC/tp2qqtWiFyZF4M8PEfks1iRFRc3alTfWrryxduVNcW0XJKZtVuoxxpgUY4HfGGNSTCoE/mlF3YBsWLvyxtqVN9auvCmu7YIEtK3E1/iNMcb8VSpk/MYYYyJY4DfGmBST1IFfRC4XkU0i8oOIjIjxvIjIlPDzX4lI03i3TXC7eofb85WIrBaRRhHPbRWR9SKyTkQK9EbDcbSrjYgcCL/2OhG5L95tE9yuYRFt2iAiARGpEn4uIZ+XiMwQkV9FZEM2zxfVsZVbu4rq2MqtXUV1bOXWrkI/tsL7PkVEPhCRb0Rko4jcHmOdxB1jqpqUP4AT2AzUATzAl0C9qHX+BrwDCNAS+F+82ya4Xa2AyuHfr8hsV/jxVqBqEX1ebYA3j2XbRLYrav2OwPJC+LwuApoCG7J5vtCPrTjbVejHVpztKvRjK552FcWxFd73SUDT8O/lge8KM34lc8bfAvhBVX9UVS8wC+gctU5n4CUNWQNUEpGT4tw2Ye1S1dWq+lv44RqgZgG9dr7alaBtC3rfvYDXCui1s6Wqq4D9OaxSFMdWru0qomMrns8rO0X6eUUplGMLQFV3q+rn4d8PAd8A0ZPvJ+wYS+bAfzKwPeLxDo7+4LJbJ55tE9muSP0JfatnUmCpiKwVkRsLqE15adf/iciXIvKOiNTP47aJbBciUha4HJgXsThRn1duiuLYyqvCOrbiVdjHVtyK8tgSkdpAE+B/UU8l7BhL5jtwSYxl0WNTs1snnm2PVdz7FpFLCP1xXhCxuLWq7hKRE4D3ROTbcNZSGO36nNDcHr+LyN+ABcCZcW6byHZl6gh8pKqRGVyiPq/cFMWxFbdCPrbiURTHVl4UybElIscR+rK5Q1UPRj8dY5MCOcaSOePfAZwS8bgmsCvOdeLZNpHtQkTOBaYDnVV1X+ZyVd0V/vdXYD6h07pCaZeqHlTV38O/vw24RaRqPNsmsl0RriHqVDyBn1duiuLYiksRHFu5KqJjKy8K/dgSETehoP+Kqr4RY5XEHWOJ6LgojB9CZys/AqfxZwdH/ah1/s5fO0c+iXfbBLerFvAD0CpqeTmgfMTvq4HLC7FdJ/LnRX0tgJ/Cn12Rfl7h9SoSqtWWK4zPK7zP2mTfWVnox1ac7Sr0YyvOdhX6sRVPu4rw2BLgJeDxHNZJ2DGWtKUeVfWLyEBgCaFe7hmqulFEbg4//wzwNqGe8R+AP4C+OW1biO26DzgeeEpEAPwamn2vOjA/vMwFvKqq7xZiu64GbhERP5AOXKOhI62oPy+Aq4Clqno4YvOEfV4i8hqhkShVRWQHMAZwR7Sp0I+tONtV6MdWnO0q9GMrznZBIR9bYa2B64D1IrIuvOweQl/cCT/GbMoGY4xJMclc4zfGGHMMLPAbY0yKscBvjDEpxgK/McakGAv8xhiTYizwm6QmImNF5K541xGR60WkRh5fo4uI1It4/ICIXHZsLc7xdURElotIBRGpnd2MkjG2GygifQu6PabkssBvUs31QJ4CP9AFyAr8qnqfqi4ryEaF/Q34Uo++dD9bIuICZgCDE9AeU0JZ4DdJR0RGheciXwbUjVh+uoi8G55U678icnbUdlcDzYFXwnOslxGR+0Tk0/Bc7NMkfMVOxDatgE7ApPA2p4vIzPC+Mudsf1hEPhaRz0SkqYgsEZHNmRehhdcbFn6dr0Tk/mzeWm9gYcRjp4g8J6H52peKSJnwvlaEX3MlcLuq/gFsFZFCm4LBJDcL/CapiEgzQvOqNAG6AudFPD0NGKSqzYC7gKcit1XVucBnQG9Vbayq6cCTqnqeqjYAygBXRm2zGlgEDAtvszlGs7ar6v8B/wVmErpKtSXwQLjN7QlNSNYCaAw0E5GLYuynNbA24vGZwFRVrQ+kAd0inqukqher6iPhx58BF8bYpzFHSdopG0zKuhCYH85yEZFF4X+PI3QTkjkRSXupOPZ3iYjcDZQFqgAbgcV5bNOi8L/rgeM0NL/6IRE5IiKVgPbhny/C6x1HKKhHz/RYJbxtpi2qmnk5/1pCc85kmh217a/A2RgTBwv8JhnFmmfEAaSpauN4dyIipQmdFTRX1e0iMhYofQztyQj/G4z4PfOxi9AkW+NV9dlc9uMXEYeqBqP2CxAgdEaSKXJeGQi1Oz1PrTYpy0o9JtmsAq4K1+fLE5pHnXCH6BYR6Q5ZI2Qaxdj+EKFb3cGfQX5v+Izh6mxeM3KbY7EE6Bd+DUTk5PAc79E2Ebqd3rE4C4hrFJAxFvhNUtHQ7epmA+sIzWX+34inewP9ReRLQiWbWLejmwk8E54RMQN4jlCJZgHwaTYvOwsYJiJfiMjpx9DmpcCrwMcish6YS+wvkrcIzSR5LFoDiRhpZEogm53TmGJCQvdTfUlV2+VxuybAUFW9LjEtMyWNZfzGFBOquht4TkQq5HHTqsDoBDTJlFCW8RtjTIqxjN8YY1KMBX5jjEkxFviNMSbFWOA3xpgUY4HfGGNSzP8DL7wItOOuOJIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(timestep,tnum,'o',label=str(p)+' Euler steps')\n",
    "plt.plot(timestep,t_ana,label='analytical')\n",
    "plt.title('Temperature')\n",
    "plt.xlabel('delta time (hr)')\n",
    "plt.ylabel('temperature (F)')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1317,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_inf(infinity,step):\n",
    "    timestep= np.linspace(0,infinity,step)\n",
    "    kconst=.398995\n",
    "    t_ana=t_amb+(t1-t_amb)*np.exp(-kconst*timestep)\n",
    "    delta_time= np.diff(timestep)\n",
    "\n",
    "    tnum=np.empty(len(timestep))\n",
    "    for i in range(0,len(timestep)-1):\n",
    "        tnum[0]=t1\n",
    "        tnum[i+1]= tnum[i] - ((tnum[i]-t_amb)*kconst*delta_time[i])\n",
    "  \n",
    "    \n",
    "    return t_ana, tnum\n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1394,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3xV9Znv8c+XJBACkZtoRZwhOh5bL4CaWi/ValGrjndH1NpTL62Mp15O6+l4qVPFzplOa72MWo8ziFan46Vo5aittV5a6+lY1OAlgkrxgiVgMYiNqOESeM4fayVsYO+wA1nZSfb3/Xrltfdee629nsXWPFm/y/NTRGBmZrahAaUOwMzMeicnCDMzy8sJwszM8nKCMDOzvJwgzMwsr8pSB9Cdtt566xg3blypwzAz6zNmz569NCJG53uvXyWIcePG0dDQUOowzMz6DEnvFHrPTUxmZpaXE4SZmeXlBGFmZnn1qz4IM+sfVq9eTVNTEytWrCh1KP1GdXU1Y8eOpaqqquhjnCDMrNdpamqitraWcePGIanU4fR5EcH7779PU1MTdXV1RR+XaROTpG9JmitpjqR7JFVLmippkaSX0p+jChx7hKR5kt6QdGlmQTbOgOt3h6nDk8fGGZmdysyKs2LFCkaNGuXk0E0kMWrUqC7fkWV2ByFpe+BCYNeIaJU0Azg1ffv6iLimk2MrgJuBw4Am4HlJD0XEq90aZOMMePhCWN2avG5ZmLwGGD+5W09lZl3j5NC9NuffM+tO6kpgsKRKoAZYXORx+wBvRMRbEbEKuBc4rtuje/J765JDu9WtyXYzszKXWYKIiEXANcCfgHeBloh4LH37fEmNkm6XNCLP4dsDC3NeN6XbNiJpiqQGSQ3Nzc1dC7KlqWvbzcy66I477uD888/f5D6LF6/7+/nrX/86r77a9QaTp556iqOPPrrLxxWSWYJIf/EfB9QBY4Ahkr4C3ALsBEwkSRzX5js8z7a8KxtFxLSIqI+I+tGj884WL2zY2K5tNzPLwIYJYvr06ey6664ljCiRZRPTocDbEdEcEauBB4D9I2JJRKyJiLXArSTNSRtqAnbIeT2W4punijfpCqgavP62qsHJdjMre8cffzx77703u+22G9OmTQNg6NChXH755UyYMIF9992XJUuWAPDwww/zuc99jj333JNDDz20Y3u75cuXU1dXx+rVqwH48MMPGTduHPfddx8NDQ2cfvrpTJw4kdbWVg4++OCOskGPPvooe+21FxMmTGDSpEkAPPfcc+y///7sueee7L///sybNy+T689ymOufgH0l1QCtwCSgQdJ2EfFuus8JwJw8xz4P7CypDlhE0rn95W6PMO2IXjLzO2wTS9GwsUlycAe1Wa9x1cNzeXXxh936mbuO2Yorj9ltk/vdfvvtjBw5ktbWVj772c9y0kkn8fHHH7Pvvvvyz//8z1x88cXceuut/OM//iOf//znmTVrFpKYPn06V199Nddeu66BpLa2loMPPphf/vKXHH/88dx7772cdNJJnHzyydx8881cc8011NfXr3f+5uZmzjnnHJ5++mnq6upYtmwZAJ/+9Kd5+umnqays5IknnuA73/kOP//5z7v13wgyTBAR8ayk+4EXgDbgRWAaMF3SRJImowXA3wNIGgNMj4ijIqJN0vnAr4EK4PaImJtJoOMnM+X/bc9W1ZX89Gufy+QUZtY33XjjjcycOROAhQsXMn/+fAYOHNjRzr/33nvz+OOPA8ncjVNOOYV3332XVatW5Z1v8PWvf52rr76a448/np/85CfceuutnZ5/1qxZHHTQQR2fNXLkSABaWlo444wzmD9/PpI67kq6W6YT5SLiSuDKDTb/9wL7LgaOynn9CPBIdtGtU1NVwYrVa3riVGbWRcX8pZ+Fp556iieeeII//OEP1NTUcPDBB7NixQqqqqo6hoxWVFTQ1tYGwAUXXMBFF13Esccey1NPPcXUqVM3+swDDjiABQsW8Lvf/Y41a9aw++67dxpDROQdnvrd736XQw45hJkzZ7JgwQIOPvjgLb7efFyLCagZWMEnq5wgzGydlpYWRowYQU1NDa+//jqzZs3a5P7bb58MtrzzzjsL7vfVr36V0047jbPOOqtjW21tLcuXL99o3/3224/f/e53vP322wAdTUy557rjjju6dF1d4QQBVA+soNUJwsxyHHHEEbS1tTF+/Hi++93vsu+++3a6/9SpUzn55JM58MAD2XrrrQvud/rpp/PBBx9w2mmndWw788wzOffcczs6qduNHj2aadOmceKJJzJhwgROOeUUAC6++GIuu+wyDjjgANasye53lyLyjh7tk+rr62NzFgz6h/te5vdvLOUPl03KICoz66rXXnuNz3zmM6UOIxP3338/Dz74ID/96U97/Nz5/l0lzY6I+nz7u1gfbmIys55xwQUX8Ktf/YpHHumR7tUt5gSBm5jMrGfcdNNNpQ6hS9wHAdRUVbJqzVra1qwtdShmZr2GEwRJExNAq4e6mpl1cIIABrcnCDczmZl1cIIABlclCcId1WZm6zhB4CYmM+t9HnroIX7wgx9s1rHjxo1j6dKlWxyDRzGxronJdxBmfVTjjGShr5ampFx/Hy+62dbWxrHHHsuxxx5b0jh8B8G6Jib3QZj1Qe1LB7csBGLd0sFbuL78ggUL+MxnPsM555zDbrvtxuGHH75RKe6lS5cybtw4ICl5cfzxx3PMMcdQV1fHj3/8Y6677jr23HNP9t13344yGW+++SZHHHEEe++9NwceeCCvv/46kMymvuiiizjkkEO45JJL1ltoaMmSJZxwwglMmDCBCRMm8MwzzwD5y5F3JycIoGZgciPlJiazPijDpYPnz5/Peeedx9y5cxk+fPgmS2rPmTOHu+++m+eee47LL7+cmpoaXnzxRfbbbz/+4z/+A4ApU6Zw0003MXv2bK655hq+8Y1vdBz/xz/+kSeeeGK9MuEAF154IV/4whd4+eWXeeGFF9htt6SA4e23387s2bNpaGjgxhtv5P3339/ia87lJiZym5jaShyJmXVZhksH19XVMXHiRCAp7b1gwYJO9z/kkEOora2ltraWYcOGccwxxwCwxx570NjYyEcffcQzzzzDySef3HHMypUrO56ffPLJVFRUbPS5v/nNbzoSTEVFBcOGDQPylyMfNWrU5l/wBpwg8DBXsz5t2Ni0eSnP9i00aNCgjucVFRW0trZSWVnJ2rXJpNoVK1YU3H/AgAEdrwcMGEBbWxtr165l+PDhvPTSS3nPN2TIkKJjK1SOvDu5iYlkPQhwJ7VZn9TDSwePGzeO2bNnA0nhva7YaqutqKur47777gOS9R5efvnlTR43adIkbrnlFgDWrFnDhx9+2OVy5Jsj0wQh6VuS5kqaI+keSdWSfiTpdUmNkmZKGl7g2AWSXpH0kqSul2jtgsEe5mrWd42fDMfcCMN2AJQ8HnNjZqOYvv3tb3PLLbew//77b9ZQ0rvuuovbbruNCRMmsNtuu/Hggw9u8pgbbriB3/72t+yxxx7svffezJ07t8vlyDdHZuW+JW0P/B7YNSJaJc0gWSFuMfCbdFnRHwJExCV5jl8A1EdE0d/A5pb7jgh2/M4jnHfw3/DtL+3S5ePNrHv153LfpdTVct9ZNzFVAoMlVQI1wOKIeCwi2nuDZwFb3lC4hSRRU+WS32ZmuTJLEBGxCLgG+BPwLtASEY9tsNvZwK8KfQTwmKTZkqYUOo+kKZIaJDU0NzdvdryDB1a6icnMLEdmCULSCOA4oA4YAwyR9JWc9y8H2oC7CnzEARGxF3AkcJ6kg/LtFBHTIqI+IupHjx692fHWDKyg1cNczXqN/rTaZW+wOf+eWTYxHQq8HRHNEbEaeADYH0DSGcDRwOlRIOqIWJw+vgfMBPbJMFYGu4nJrNeorq7m/fffd5LoJhHB+++/T3V1dZeOy3IexJ+AfSXVAK3AJKBB0hHAJcAXIuKTfAdKGgIMiIjl6fPDgS2fFtmJwQMr3MRk1kuMHTuWpqYmtqTZ2NZXXV3N2LFd6/LNLEFExLOS7gdeIGlKehGYBswFBgGPSwKYFRHnShoDTI+Io4BtgZnp+5XA3RHxaFax0jiDW5ddzqg1zXB93y/0ZdbXVVVVUVdXV+owyl5mw1xLYbOGubYX+sqt5VI1ONNx1GZmvUUph7n2fhkW+jIz68ucIDIs9GVm1pc5QRQq6NUNhb7MzPoyJ4geLvRlZtZXOEGkhb6WV2/H2hBrtxrrDmozM5wgEuMn8/Ahj7HjyrtY8rUGJwczM5wgOgwZlJT8/nily22YmYETRIehg5I5gx+v9GxqMzNwguhQM7A9QfgOwswMnCA6dNxBuGCfmRngBNGhxn0QZmbrcYJIrbuDcIIwMwMniA41A30HYWaWywkita6T2n0QZmbgBNGhYoAYXFXhOwgzs5QTRI4hgyrdB2Fmlso0QUj6lqS5kuZIukdStaSRkh6XND99HFHg2CMkzZP0hqRLs4yz3dBBFW5iMjNLZZYgJG0PXAjUR8TuQAVwKnAp8GRE7Aw8mb7e8NgK4GbgSGBX4DRJu2YVa7uagZVuYjIzS2XdxFQJDJZUCdQAi4HjgDvT9+8Ejs9z3D7AGxHxVkSsAu5Nj8vUUDcxmZl1yCxBRMQi4BrgT8C7QEtEPAZsGxHvpvu8C2yT5/DtgYU5r5vSbRuRNEVSg6SG5ubmLYq5xk1MZmYdsmxiGkHyV38dMAYYIukrxR6eZ1vk2zEipkVEfUTUjx49evOCTbmT2sxsnSybmA4F3o6I5ohYDTwA7A8skbQdQPr4Xp5jm4Adcl6PJWmeytSQgR7mambWLssE8SdgX0k1kgRMAl4DHgLOSPc5A3gwz7HPAztLqpM0kKRz+6EMYwWSO4hP3MRkZgYknciZiIhnJd0PvAC0AS8C04ChwAxJXyNJIicDSBoDTI+IoyKiTdL5wK9JRj/dHhFzs4q13dBBlXy0qo2IIMlpZmblK7MEARARVwJXbrB5JcndxIb7LgaOynn9CPBIlvFtaOigSiKSkt/txfvMzMqVZ1K3a5zBV/7wt7w16MtU/3gCNM4odURmZiXlP5MhSQYPX8iQ1a0gGLC8CR6+MHlv/OTSxmZmViK+gwB48nuwunX9batbk+1mZmXKCQKgpalr283MyoATBMCwsV3bbmZWBpwgACZdAVWD199WNTjZbmZWpjbZSS2pHjiQpFxGKzAHeCIilmUcW89JO6LXPnEVtCzi4+pPUfu333MHtZmVtYJ3EJLOlPQCcBkwGJhHUhbj88Djku6U9Fc9E2YPGD8ZvjmHnVbdxa17P+jkYGZlr7M7iCHAARHRmu9NSROBnUlmQ/cLAwaIoYMq+XCF6zGZmRVMEBFxc2cHRsRL3R9O6dUOqmS5E4SZWadNTI/lPL+sZ8IpvdrqKpavWF3qMMzMSq6zUUy5iyucnHUgvUVtdSUfueS3mVmnCSLvAj39XW21m5jMzKDzTuodJT1Esrpb+/MOEXFsppGVyNDqKt5e+nGpwzAzK7nOEsRxOc+vyTqQ3sJNTGZmic5GMf1uSz5Y0i7Az3I27QhcAewH7JJuGw78JSIm5jl+AbAcWAO0RUT9lsRTrNpqD3M1M4NOEoSkh0lWgHs0XVM6970dgTOBBRFxe77jI2IeMDHdvwJYBMyMiH/N+ZxrgZZO4jskIpYWdyndo3ZQJava1rKybQ2DKit68tRmZr1KZ01M5wAXAf8qaRnQDFQD44A3gR9HRL71pPOZBLwZEe+0b0jXqZ4MfHEz4s5MbXUVAB+taGPQUCcIMytfnTUx/Rm4GLhY0jhgO5JaTH+MiE+6eJ5TgXs22HYgsCQi5hcKAXhMUgD/HhHT8u0kaQowBeCv/mrLK3/UVif/JMtXtDFq6KAt/jwzs76qqBXlImIBsGBzTiBpIHAsSU2nXKexcdLIdUBELJa0DUntp9cj4uk8sU0jaQqjvr5+i4fmtt9BeKirmZW7nij3fSTwQkQsad8gqRI4kfU7sdcTEYvTx/eAmcA+GccJwFbpHURLq2dTm1l564kEke9O4VDg9YjIu2SbpCGSatufA4eTlBnP3PCagYAThJlZUQlC0uB02GqXSKoBDgMe2OCtjfokJI2R9Ej6clvg95JeBp4DfhkRj3b1/F3WOIOd7tqXtwZ9mYMfOQQaZ2R+SjOz3qqYBYOOIZkoNxCoS8t8f6+YmdRpZ/aoPNvPzLNtMXBU+vwtYMKmPr9bNc6Ahy+kcnUrCIaseBcevjB5z2tDmFkZKuYOYipJ+/9foKPM97jsQiqRJ78HqzdY+mJ1a7LdzKwMFZMg2iKis8ls/UNL3u6QwtvNzPq5YhLEHElfBiok7SzpJuCZjOPqecPGdm27mVk/V0yCuADYDVgJ3E1SGuObWQZVEpOugKrB62+rGpxsNzMrQ512Uqc1lK6KiH8ALu+ZkEqkvSP6ye8RLU28N2Brtj3m++6gNrOy1WmCiIg1kvbuqWBKbvxkGD+Z//Wzl3huwTJ+P75XlYkyM+tRxZTaeDFdLOg+oGMlnYjYcG5Dv7HV4CpaPvFEOTMrb8UkiJHA+6xfdTXYePJbvzG8porlK9tYszaoGKBSh2NmVhKbTBARcVZPBNKbDBucFOz7sHU1I4YMLHE0ZmalUcxM6p+Q3DGsJyLOziSiXqA9QfzFCcLMylgxTUy/yHleDZwALM4mnN6hPUG4YJ+ZlbNimph+nvta0j3AE5lF1AsMr3GCMDPbnHLfOwNbvnRbL9bRxPTJqhJHYmZWOsX0QSxn/T6IPwOXZBZRL7BVTie1mVm5KqaJqbYnAulNRqSLBn3guRBmVsY22cQk6clitvUnVRUDqK2uZNnHbmIys/JVMEFIqpY0Etha0ghJI9OfccCYTX2wpF0kvZTz86Gkb0qaKmlRzvajChx/hKR5kt6QdOnmXuDmGjlkoBOEmZW1zpqY/p6kausYYDbQPqX4Q+DmTX1wRMwDJkJH0b9FwEzgLOD6iLim0LHp/jeTLFfaBDwv6aGIeHVT5+0WjTO4r/VyRs9rhuvHJhVdXbTPzMpMwQQRETcAN0i6ICJu2sLzTALejIh3pKJKV+wDvJEuPYqke4HjgOwTRLr06DZr09XlWhZ66VEzK0ub7IOIiJsk7S5psqSvtv908TynAvfkvD5fUqOk2yWNyLP/9sDCnNdN6baNSJoiqUFSQ3NzcxfDysNLj5qZAcV1Ul8J3JT+HAJcDRxb7AkkDUz3vy/ddAuwE0nz07vAtfkOy7Nto3IfABExLSLqI6J+9OjRxYZVmJceNTMDipso93ckTUR/Tgv3TQAGdeEcRwIvRMQSgIhYEhFrImItcCtJc9KGmoAdcl6PpafKe3jpUTMzoLgE0Zr+Mm+TtBXwHrBjF85xGjnNS5K2y3nvBGBOnmOeB3aWVJfegZwKPNSFc24+Lz1qZgYUlyAaJA0n+Wt/NvAC8FwxHy6phmQkUu7aEVdLekVSI0mT1bfSfcdIegQgItqA84FfA68BMyJibnGXtIXGT4ZjbqS1ZgxrQ6wauj0cc6M7qM2s7Cgib9N+8mYy5GhsRCxMX48DtoqIxh6Jrovq6+ujoaGhWz5r9jsfcNItz3DHWZ/l4F226ZbPNDPrbSTNjoj6fO91egcRSfb4vzmvF/TW5NDdRqbrQHiynJmVq2KamGZJ+mzmkfQyI2ucIMysvBWzYNAhwLmSFgAfkwxBjYgYn2VgpbbV4EoqBogPXPLbzMpUMQniyMyj6IUkMaLG9ZjMrHwVM5P6HZI5CV9Mn39SzHH9waghA1n6kROEmZWnYmdSXwJclm6qAv4zy6B6i9G1g1j60cpSh2FmVhLF3AmcQFIq42OAiFgMlMUiQqNrB9G83AnCzMpTMQliVTrcNQAkDck2pF6icQZT3zyFpz85kbh+96TKq5lZGSkmQcyQ9O/AcEnnAE+QzKruv9KS38NW/ZkBCtRe8ttJwszKSDGd1NcA9wM/B/4bcEU3rA/Ru7nkt5lZUcNcAV4BBpM0M72SXTi9hEt+m5kVNYrp6yTF+U4kKf09S9LZWQdWUi75bWZWVB/EPwB7RsSZEXEGsDfJsNf+yyW/zcyKShBNwPKc18tZfznQ/ict+R3DdmBtiJaBn3LJbzMrO8X0QSwCnpX0IEkfxHHAc5IuAoiI6zKMr3TGT0bjJ7Pf95/kwJ235kfjJ5Q6IjOzHlVMgngz/Wn3YPpYPpPlPJvazMrQJhNERFy1OR8saRfgZzmbdgSuALYHjgFWkSSesyLiL3mOX0DSnLUGaCu0oEXWRtcO4r3lK0pxajOzkipmFFO9pJmSXpDU2P6zqeMiYl5ETIyIiSQd258AM4HHgd3TcuF/ZF2Np3wOST+jJMkBYPTQQSz50HcQZlZ+imliuotkJNMrwNrNPM8k4M20Guw7OdtnkQyd7bU+NayapR+tZPWatVRVlEURWzMzoLhRTM0R8VBEvB0R77T/dPE8pwL35Nl+NvCrAscE8Jik2ZKmFPpgSVMkNUhqaG5u7mJYm7bdsGoi4D0X7TOzMlPMHcSVkqYDTwIdvyUj4oFiTiBpIEk12Ms22H450EZyh5LPARGxWNI2wOOSXo+IpzfcKSKmAdMA6uvro5iYuuJTw6oB+HNLK9sPH7yJvc3M+o9iEsRZwKdJ1oFob2IKoKgEQbIi3QsRsaR9g6QzgKOBSWml2I2kZcWJiPckzQT2ATZKEJlqnMEBj03lrUGLWHHvGDjyKs+FMLOyUUyCmBARe2zBOU4jp3lJ0hEkM7G/EBGf5DsgLSk+ICKWp88PB3q2Ul5a0bVqdSsIaloXJxVdwUnCzMpCMX0QsyTtujkfLqkGOIz17zZ+TDKH4nFJL0n6t3TfMZIeSffZFvi9pJdJ6kD9MiIe3ZwYNpsruppZmSvmDuLzwBmS3ibpgxAQ6TDVTqV3CKM22PY3BfZdDByVPn8LKO3UZVd0NbMyV0yCOCLzKHqjYWOhJU/JKVd0NbMyUcyCQe8AOwBfTJ9/UsxxfZ4ruppZmStmJvWVJJ3K7cNUq4D/zDKoXiGt6MqwHQjEotiatUe7oquZlY9i7gROIJnH8DF09BWURaE+xk+Gb83hP494mQNW3sh7444tdURmZj2mmASxKp2rENAxBLWsjB2RNDU1fZB3VK6ZWb9UTIKYIenfgeGSzgGeAKZnG1bvssOIGgAWOkGYWRkpptz3NZIOAz4EdgGuiIjHM4+sF2m/g1i4rHUTe5qZ9R+bTBCSfhgRl5CU6d5wW//XOIPqJ7/HW9VNtMzaBkb/b3dUm1lZKKaJ6bA8247s7kB6pbTcBi0LGUAwYvWS5HXjjFJHZmaWuYIJQtL/kPQKsEvuQkHpjOpNLhjUL7jchpmVsc6amO4mWavhX4BLc7Yvj4hlmUbVW7jchpmVsYIJIiJagBaSaqzlyeU2zKyM9f+SGVsiT7mNNRUut2Fm5cEJojMblNtoWrs1s3a70qOYzKwsFFPNtbyNnwzjJxNrgy9e8ShnDPprDih1TGZmPSCzOwhJu6QLArX/fCjpm5JGSnpc0vz0cUSB44+QNE/SG5IuzbdPTxowQNSNGsJbzR+XOhQzsx6RWYKIiHkRMTEiJgJ7k5QJn0kyIurJiNgZeJL1R0gBIKkCuJlkvsWuwGmbu6pdt2mcwd0ff41bFxwG1+/uuRBm1u/1VB/EJODNdD2J44A70+13Asfn2X8f4I2IeCsiVgH3pseVRjphblTbewwgkpFNnjBnZv1cTyWIU4F70ufbRsS7AOnjNnn23x7IHV/alG7biKQpkhokNTQ3N3djyDk8Yc7MylDmCULSQJL1JO7rymF5tkW+HSNiWkTUR0T96NGjNyfETfOEOTMrQz1xB3Ek8EJELElfL5G0HUD6+F6eY5pIljltNxZYnGmUnSk0Mc4T5sysH+uJBHEa65qXAB4CzkifnwE8mOeY54GdJdWldyCnpseVhtenNrMylGmCkFRDUg32gZzNPwAOkzQ/fe8H6b5jJD0CEBFtwPnAr4HXgBkRMTfLWDuVM2FuLeK9AaOT154wZ2b9mJLVRPuH+vr6aGhoyPQcVz08l3ufW8jcq77EgAH5ukrMzPoOSbMjoj7fey610RWNM/j2qycxd8AprL1+Nw9zNbN+zaU2ipXOhRiyuhUEA5YvSuZCgJuazKxf8h1EsTwXwszKjBNEsTwXwszKjBNEsTwXwszKjBNEsfLMhYhKz4Uws/7LCaJY7XMhBo8kgAhYPaC61FGZmWXGo5i6qq01KRQlGLjqA49kMrN+y3cQXeGRTGZWRpwgusIjmcysjDhBdIVHMplZGXGC6Io8I5naKqo9ksnM+iUniK7YcCQTsIJBpY7KzCwTHsW0OdpHMgFD17R4JJOZ9Uu+g+gqj2QyszLhBNFVHslkZmUi0yYmScOB6cDuJE32ZwPfBHZJdxkO/CUiJuY5dgGwHFgDtBVa0KLHDRsLLQvzbzcz60ey7oO4AXg0Iv4uXVu6JiJOaX9T0rVASyfHHxIRSzOOsWsmXZH0OeQ0M60F+JvDfTtmZv1KZr/TJG0FHATcBhARqyLiLznvC5gM3JNVDJkYPxkmfBlYt9zoACBeutsrzJlZv5LlH707As3ATyS9KGm6pCE57x8ILImI+QWOD+AxSbMlTSl0EklTJDVIamhubu6+6Dsz/7E0vHUq1rij2sz6lywTRCWwF3BLROwJfAxcmvP+aXR+93BAROwFHAmcJ+mgfDtFxLSIqI+I+tGjR3dT6JtQoEM63FFtZv1IlgmiCWiKiGfT1/eTJAwkVQInAj8rdHBELE4f3wNmAvtkGGvXFOqQ1gA3M5lZv5FZgoiIPwMLJbWPWJoEvJo+PxR4PSLy/sktaYik2vbnwOHAnKxi7bI8JTcAFGuSDmwnCTPrB7IeeHMBcJekRmAi8P10+6ls0LwkaYykR9KX2wK/l/Qy8Bzwy4h4NONYi9deckMVG7/nSXNm1k8oIja9Vx9RX18fDQ0NPXfCqcPZsLN63Xudjd41M+sdJM0uNM/MQ/e3RMHJcXIzk5n1eU4QW2LSFeTOh1gn3MxkZn2eE8SWGD+Zgk1M+cpxmJn1IU4QW2rYDnk3h5uZzKyPc2qbtoYAAAilSURBVILYUgWamUTAry7p+XjMzLqJE8SW6qSZKVqX+S7CzPosJ4juUKCZSUDMPNdJwsz6JCeI7jDpioJvKdbAA1PgFxf1YEBmZlvOCaI7jJ8Mg0d2skNAw+2+kzCzPsUJorsc+cO89ZnWcae1mfUtThDdpbP6TKloXeamJjPrM5wgutP4yXDCv5F/dnXaad1wm5OEmfUJThDdbfxkqD+74NsdSWLqMPhhnfslzKzXcoLIwtHXddpp3XF/0boMHjgHvj/GicLMeh0niKwc+UMKNTVtZNXHSaKYOsx3FmbWa1SWOoB+a/xk+NMsaLit68e231k8cM7G72kAxNpkct6kK9KZ3GZm3S/TBCFpODAd2J2kHsXZwJeAc4DmdLfvRMQjeY49ArgBqACmR8QPsow1E0dflzw23E7Bqq9dFWuTx5aFhZOImZWfwSOTlotu/KMx6yamG4BHI+LTwATgtXT79RExMf3JlxwqgJuBI4FdgdMk7ZpxrNk4+jo4cdomJtKZmW2h1mXw4Hnd2jydWYKQtBVwEHAbQESsioi/FHn4PsAbEfFWRKwC7gWOyybSHjB+MlzyNpx4K1QNKXU0ZtZfrVnVrYuVZXkHsSNJM9JPJL0oabqk9t+O50tqlHS7pBF5jt0eyF1xpyndthFJUyQ1SGpobm7Ot0vvMX4yXL44SRS+ozCzLLQ0ddtHZZkgKoG9gFsiYk/gY+BS4BZgJ2Ai8C5wbZ5jC6zjmWdjxLSIqI+I+tGjR3dL4Jlrv6OY2uJkYWbda9jYbvuoLBNEE9AUEc+mr+8H9oqIJRGxJiLWAreSNCflOza3hvZYYHGGsZaOk4WZdZeKgZ1Wl+6qzEYxRcSfJS2UtEtEzAMmAa9K2i4i3k13OwGYk+fw54GdJdUBi4BTgS9nFWuvMX7ypkcgNM5Iiv61LuuZmMysb8hgFFPW8yAuAO6SNBB4CzgLuFHSRJImowXA3wNIGkMynPWoiGiTdD7wa5JhrrdHxNyMY+0bikkiZmbdQBHdND6/F6ivr4+GhoZSh2Fm1mdImh0R9fnec6kNMzPLywnCzMzycoIwM7O8nCDMzCyvftVJLakZeGczD98aWNqN4fQFvuby4GsuD5t7zX8dEXlnGferBLElJDUU6snvr3zN5cHXXB6yuGY3MZmZWV5OEGZmlpcTxDrTSh1ACfiay4OvuTx0+zW7D8LMzPLyHYSZmeXlBGFmZnmVfYKQdISkeZLekHRpqePJiqQFkl6R9JKkhnTbSEmPS5qfPuZb3a9PSVcpfE/SnJxtBa9T0mXpdz9P0pdKE/WWKXDNUyUtSr/vlyQdlfNen75mSTtI+q2k1yTNlfQ/0+39/XsudN3ZfdcRUbY/JKXE3yRZHnUg8DKwa6njyuhaFwBbb7DtauDS9PmlwA9LHWc3XOdBJCsZztnUdQK7pt/5IKAu/W+hotTX0E3XPBX4dp59+/w1A9uRLD4GUAv8Mb2u/v49F7ruzL7rcr+D2Ad4IyLeiohVwL3AcSWOqScdB9yZPr8TOL6EsXSLiHga2HA1pULXeRxwb0SsjIi3gTfIv8Jhr1bgmgvp89ccEe9GxAvp8+XAayRr1vf377nQdReyxddd7glie2BhzusmOv8H78sCeEzSbElT0m3bRrq6X/q4Tcmiy1ah6+zv3//5khrTJqj25pZ+dc2SxgF7As9SRt/zBtcNGX3X5Z4glGdbfx33e0BE7AUcCZwn6aBSB9QL9Ofv/xZgJ2Ai8C5wbbq931yzpKHAz4FvRsSHne2aZ1ufvGbIe92ZfdflniCagB1yXo8FFpcolkxFxOL08T1gJsmt5hJJ2wGkj++VLsJMFbrOfvv9R8SSiFgTEWuBW1nXtNAvrllSFckvybsi4oF0c7//nvNdd5bfdbkniOeBnSXVpetmnwo8VOKYup2kIZJq258DhwNzSK71jHS3M4AHSxNh5gpd50PAqZIGSaoDdgaeK0F83a79F2XqBJLvG/rBNUsScBvwWkRcl/NWv/6eC113pt91qXvmS/0DHEUyGuBN4PJSx5PRNe5IMprhZWBu+3UCo4Angfnp48hSx9oN13oPyW32apK/oL7W2XUCl6ff/TzgyFLH343X/FPgFaAx/UWxXX+5ZuDzJE0ljcBL6c9RZfA9F7ruzL5rl9owM7O8yr2JyczMCnCCMDOzvJwgzMwsLycIMzPLywnCzMzycoIwyyFpuKRvpM/HSLo/w3NNzK28adbbOEGYrW848A1IZp9HxN9leK6JJOPYzXolz4MwyyGpvaLvPJIJV5+JiN0lnUlSHbQC2J2k3s1A4L8DK4GjImKZpJ2Am4HRwCfAORHxuqSTgSuBNUALcChJdc3BwCLgX4BfADcBewCVwNSIeDA99wmsK9t8d0RclfE/hRmVpQ7ArJe5FNg9IiamFTN/kfPe7iQVNKtJfrlfEhF7Sroe+CrwryQLx58bEfMlfQ74P8AXgSuAL0XEIknDI2KVpCuA+og4H0DS94HfRMTZkoYDz0l6Ij33Pun5PwGel/TLiGjI8h/CzAnCrHi/jaQO/3JJLcDD6fZXgPFplc39gfuSsjlA8lc/wH8Bd0iaATxAfocDx0r6dvq6Gvir9PnjEfE+gKQHSMouOEFYppwgzIq3Muf52pzXa0n+XxoA/CUiJm54YEScm95R/C3wkqSN9iEpz3xSRMxbb2Ny3IZtwW4btsy5k9psfctJlnPsskhq87+d9jegxIT0+U4R8WxEXAEsJSnDvOG5fg1ckFbtRNKeOe8dlq65PJikL+S/NidGs65wgjDLkTbj/JekOcCPNuMjTge+Jqm9cm77ErY/kvRK+rlPk1TW/S2wa7rQ/CnAPwFVQGO63z/lfO7vSap2vgT83P0P1hM8ismsl0tHMXV0Zpv1FN9BmJlZXr6DMDOzvHwHYWZmeTlBmJlZXk4QZmaWlxOEmZnl5QRhZmZ5/X8yoOkSlK0gXwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ta, tn =time_inf(100,250)\n",
    "plt.plot(ta,'-',label='analytical')\n",
    "plt.plot(tn,'o',label='numerical')\n",
    "plt.xlabel('timestep')\n",
    "plt.ylabel('temperature (F)')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final temperature as t$\\rightarrow\\infty$ is 65$^{o}$F."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1520,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The time of death was at 907 AM\n"
     ]
    }
   ],
   "source": [
    "tnew=98.6 # the new temperature\n",
    "kconst=.275\n",
    "dtnew= tnew-t_amb \n",
    "hours= np.log(dt1/dtnew)/kconst # getting how many hours to add on to the existing time\n",
    "mins= 60\n",
    "hr=12\n",
    "minutes= hours*mins # using 24 hour clock 1100 would already have 12 hours mulitplied by 60 minutes\n",
    "time_mins=(hr*mins)+minutes #the answer for time_mins= 642 knowing that at 0900 there are 600 minutes, I just add 300 for the actual time.\n",
    "act_time= time_mins+300\n",
    "\n",
    "print('The time of death was at {:.0f} AM'.format(act_time))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Now that we have a working numerical model, we can look at the results if the 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) *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. Compare the new nonlinear Euler approximation to the linear analytical model. At what time was the corpse 98.6$^{o}$F? i.e. what was the time of death?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1383,
   "metadata": {},
   "outputs": [],
   "source": [
    "def temp_time(number):\n",
    "    \n",
    "    \n",
    "    if number==-3 or number < -2:\n",
    "        print(-3, 'hours =',8,'AM', 55,'F')\n",
    "    elif number== -2 or number < -1:\n",
    "        print(-2, 'hours =',9,'AM', 58,'F')\n",
    "    elif number== -1 or number < 0:\n",
    "        print(-1, 'hours =',10,'AM', 60,'F')\n",
    "    elif number==0 or number < 1:\n",
    "        print(0, 'hours =',11,'AM', 65,'F')\n",
    "    elif number == 1 or number < 2:\n",
    "        print(1, 'hour =',12,'PM', 66,'F')\n",
    "    else:\n",
    "        print(2, 'hours =',13,'PM', 67,'F');\n",
    "        \n",
    "    return\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1384,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-3 hours = 8 AM 55 F\n",
      "-2 hours = 9 AM 58 F\n",
      "-1 hours = 10 AM 60 F\n",
      "0 hours = 11 AM 65 F\n",
      "1 hour = 12 PM 66 F\n",
      "2 hours = 13 PM 67 F\n"
     ]
    }
   ],
   "source": [
    "temp_time(-3)\n",
    "temp_time(-2)\n",
    "temp_time(-1)\n",
    "temp_time(0)\n",
    "temp_time(1)\n",
    "temp_time(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It does look correct,  but I find it easier plotting with an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1416,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature= np.array([55,58,60,65,66,67])\n",
    "clock= np.array([8,9,10,11,12,13])\n",
    "other= np.array([-3,-2,-1,0,1,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1417,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5xN9f7H8dfH/R4at1waoplSKJdEIlKSzBxUunfEVE4pTnckUelwqtNNRyp1jvpFZdwiQqJcmlCiGZG7MJI7Yy7f3x97cybNTIPZe+09+/18POYxe69Z1vf9NebjO9+11neZcw4REYkcRbwOICIiwaXCLyISYVT4RUQijAq/iEiEUeEXEYkwKvwiIhFGhV8kCMyspJkdMLOzvM4iosIvYc9fUI99ZJnZ4Wzvb/Eo02Izu/XYe+dcmnOunHNumxd5RLIr5nUAkdPlnCt37LWZbQB6O+c+9y6RSGjTiF8KPTNrbWZLzGyvmW0zsxfNLMdBj5l9YWa9T9iWYmadzayomb1qZqn+Y31nZjE5HOOfQHNgrP+3jn+aWSkzc2ZWy7/P/5nZv8xstpkd9Ldb1cxeN7M9ZrbKzC7MdszaZjbZzHaZ2c9mdk/B/i1JJFHhl0iQDtwHVAbaANcBvXPZ933gpmNvzKyp/8/NAroATYFzgErAzcBvJx7AOfd34Bt8v3mU87/PyY3AQ0AUvt++FwPzgTOBT4F/+DMU9b//GjgL6AQ8YWZt89V7kROo8Euh55xb6pz7xjmX6ZxbB4wFciuaHwGtzKyG//3NwETnXAa+/0AqALG+w7pVzrmdpxFtonPuO+fcYWAysNc596FzLhOYAFzk3+8yoJRz7nnn3FHn3BrgHaDnabQtEUyFXwo9MzvfzGaY2Q4z2wc8iW+U/QfOud3AbOAGMyuCb1Q+3v/lGcBbwL+BHf5pmXI5HSefdmR7fTiH98eOfTYQ7Z8C2mNme4ABQPXTaFsimAq/RII3gWXAOc65CsDTgOWx/wf4pnvaAhn4plhwPi845y4CGgGNgQdyOUZBLnu7GUh2zlXM9lHeOfeXAmxDIogKv0SC8vimUQ6YWUOgz5/sPxloCAwEPnD+tcvNrKWZNfOfGD4IHAUycznGDqBegaSHhf72H/SfJC5mZo3M7OICOr5EGBV+iQT9gd5mdgB4Dfgwr52dc4eAKUAHfCd7j6kIjAP2AD8DG4GXcznMi8DtZvabmf3jdMI759KBzkArf5upwGj+NxUkclJMD2IREYksGvGLiEQYFX4RkQijwi8iEmFU+EVEIkxYLNIWFRXloqOjvY4hIhJWvv32213OuSonbg+Lwh8dHU1SUpLXMUREwoqZbcxpu6Z6REQijAq/iEiEUeEXEYkwKvwiIhFGhV9EJMIErPCb2dtmttPMfsi2bZiZfW9mK8xslpmdFaj2RUQkZ4Ec8Y/D94i47EY65xo555oA0/A9EENERIIoYIXfOfclsPuEbfuyvS1LwT6sQkSk0Dicfph+M/qx/cD2Aj920Of4zewZM9sM3EIeI34zSzCzJDNLSk1NDV5AEZEQ8PT8p3ll6SusTl1d4McOeuF3zg10ztXG9xzT+/LYb4xzrplzrlmVKn+441hEpNBasX0FI78eSa8mvWhft32BH9/Lq3reB7p72L6ISMjJyMqg95TeRJWJYtRVowLSRlALv5k1yPa2K5AczPZFRELdS4tf4ttfvuXVzq9SqXSlgLQRsEXazOwDoB0QZWZbgCFAZzOLAbLwPTv0nkC1LyISbtbtXseT854kLiaO7ucFbkIkYIXfOXdTDpvfClR7IiLhzDnH3dPupnjR4rzW+TXMLGBthcWyzCIihd27373LnPVzGH3taGpWqBnQtrRkg4iIx3Yc2MGAzwbQpk4bEpomBLw9FX4REY/1m9mPg+kHefO6NyligS/LmuoREfHQlJQpTFg1geFXDCcmKiYobWrELyLikX1p++g7vS8XVr2Qh1s/HLR2NeIXEfHIY58/xrb92/jkxk8oUbRE0NrViF9ExAMLNy1kdNJoHrjkAVrUbBHUtlX4RUSC7EjGEfpM7UN0xWiGtx8e9PY11SMiEmTPLniW5F3JzLxlJmVLlA16+xrxi4gE0codK3lu4XPc1ug2rq5/tScZVPhFRIIkMyuTPlP7ULFURV64+gXPcmiqR0QkSF5d+ipLti7h/W7vE1UmyrMcGvGLiATBhj0bGDh3IJ0bdKbnBT09zaLCLyISYM457pnmW4V+9LWjA7ryZn5oqkdEJMDGrxzPZ+s+4+VOL1PnjDpex9GIX0QkkFIPpvLgzAdpWaslfZv39ToOoMIvIhJQ/T/rz760fYy9bixFixT1Og6gwi8iEjAzfprB+JXjeaLNEzSs2tDrOMep8IuIBMD+tP3cM/0ezos6j8cve9zrOL+jk7siIgEwaO4gNu/dzMJeCylZrKTXcX5HI34RkQK2eMtiXln6Cn2b96VV7VZex/kDFX4RkQJ0NPMovaf0pmaFmjzX4Tmv4+RIUz0iIgVoxMIRrEpdxdSbplK+ZHmv4+RII34RkQKyOnU1w78cTs8LetLl3C5ex8mVCr+ISAHIcln0mdqH8iXL869O//I6Tp401SMiUgDeSHqDrzd/zbi4cVQtW9XrOHnSiF9E5DRt3ruZxz5/jI71OnJ749u9jvOnVPhFRE6Dc46+n/Yl02Xy7y7/9nzlzfzQVI+IyGmYsGoC09ZM459X/ZO6lep6HSdfNOIXETlFvx76lftn3E+zs5rR75J+XsfJN434RURO0UOzH+K3I78x+7rZFCsSPuU0YCN+M3vbzHaa2Q/Zto00s2Qz+97MJplZxUC1LyISSLPXzWbcinE80uoRGldv7HWckxLIqZ5xQKcTts0GLnDONQLWAKG1ZJ2ISD4cPHqQu6fdzblnnsvgtoO9jnPSAlb4nXNfArtP2DbLOZfhf7sYqBWo9kVEAmXIF0NYv2c9Y7qMoVSxUl7HOWlentztBczI7YtmlmBmSWaWlJqaGsRYIiK5S9qWxIuLXyTh4gTaRrf1Os4p8aTwm9lAIAMYn9s+zrkxzrlmzrlmVapUCV44EZFcpGemc9eUu6hWthr/6PgPr+OcsqCfhjazO4AuQAfnnAt2+yIip2rU16P4fsf3TLpxEmeUOsPrOKcsqIXfzDoBjwJtnXOHgtm2iMjpWPPrGobOH0r387oTHxvvdZzTEsjLOT8AFgExZrbFzO4CXgXKA7PNbIWZvRGo9kVECkqWyyJhagKli5fmlWte8TrOaQvYiN85d1MOm98KVHsiIoEydtlY5m+cz9jrxlKjfA2v45w2LdkgIpKHbfu38fDsh7ki+gp6XdTL6zgFQoVfRCQP9316H0czj4bNypv5ET6LS4iIBNknP37CpORJjOgwggZnNvA6ToHRiF9EJAd7juzhb5/+jSbVmzDg0gFexylQGvGLiOTg4VkPk3owlek3T6d40eJexylQGvGLiJxg3vp5jF0+lgGXDuDiGhd7HafAqfCLiGRzOP0wCdMSqFepHk+1e8rrOAGhqR4RkWyenv80a3ev5fPbPqdM8TJexwkIjfhFRPyW/7KckV+PpFeTXnSo18HrOAGjwi8iAmRkZdB7am+iykQx8qqRXscJKE31iIgALy1+iWW/LGNCjwlULl3Z6zgBpRG/iES8dbvX8eS8J+ka05Ue5/fwOk7AqfCLSERzzpEwLYHiRYvzeufXC82yDHnRVI+IRLRxK8Yxd/1cRl87mpoVanodJyg04heRiLX9wHb+PuvvXFbnMhKaJngdJ2hU+EUkYj0w8wEOph/kzevepIhFTjmMnJ6KiGQzJWUKE1ZNYPDlg4mNivU6TlCp8ItIxNl7ZC99p/flwqoX8kjrR7yOE3Q6uSsiEefxOY+zbf82PrnxE0oULeF1nKDTiF9EIsrCTQsZnTSaBy55gBY1W3gdxxMq/CISMY5kHKHP1D6cfcbZDGs/zOs4ntFUj4hEjGe+fIbkXcnMvGUm5UqU8zqOZzTiF5GIsHLHSkZ8NYLbGt3G1fWv9jqOp1T4RaTQy8zKpPfU3lQsVZEXrn7B6zie01SPiBR6ry59laVblzK+23iiykR5HcdzGvGLSKG2Yc8GBs4dSOcGnbnpgpu8jhMSVPhFpNByznHPtHsAGH3t6IhYeTM/NNUjIoXW+JXj+WzdZ7zc6WXqnFHH6zghQyN+ESmUUg+m8uDMB2lZqyV9m/f1Ok5IUeEXkUKp/2f92Ze2j7HXjaVokaJexwkpASv8Zva2me00sx+ybbvezFaZWZaZNQtU2yIS2Wb8NIPxK8fzRJsnaFi1oddxQk4gR/zjgE4nbPsB6AZ8GcB2RSSC7U/bz93T7ua8qPN4/LLHvY4TkgJ2ctc596WZRZ+w7UdAZ9ZFJCBSD6byyOePsGXfFhb2WkjJYiW9jhSSQvaqHjNLABIA6tTR2XgRydnPv/1MYnIiicmJfLX5K7JcFgNaDqBV7VZeRwtZIVv4nXNjgDEAzZo1cx7HEZEQ4Zxj+fblx4v9yp0rAWhcrTGDLx9MXEwcTao38ThlaAvZwi8ickx6ZjpfbvySxOREJqdMZvO+zRSxIrSp04YXr36RuJg46laq63XMsKHCLyIh6cDRA8xcO5PJKZOZtmYae47soXSx0lxd/2qGXTGMa8+9VuvunKKAFX4z+wBoB0SZ2RZgCLAbeAWoAkw3sxXOucheH1VEjttxYAdT10wlMTmRz3/+nLTMNM4sfSZ/if0LcTFxdDynI2WKl/E6ZtgL5FU9ua2GNClQbYpI+Pnp15988/UpiSzavAiHo27FuvRt3pf42Hha1W5FsSKanChI+tsUkaDKclkkbUticvJkElMSWZ26GoCLa1zM0HZDiY+N54KqF+iy7wBS4ReRgDuaeZQvNnxx/OTstv3bKGpFaRvdlnua3kNcbJwWUQsiFX4RCYh9afuY8dMMJqdMZvpP09mXto+yxcvSqX4n4mPj6dygM5VLV/Y6ZkRS4ReRAvPL/l+YkjKFxJRE5vw8h/SsdKqUqcIN599AXGwcHep2oHTx0l7HjHgq/CJyWpJ3JR+/mWrJ1iUAnFPpHB645AHiY+NpWaulVscMMSr8InJSslwWS7YsOT5fn/JrCgDNz2rOM+2fIS4mjvOrnK+TsyHsTwu/mZ0LjAaqOecuMLNGQFfn3PCApxORkJCWkcbc9XNJTE5kypopbD+wnWJFitG+bnv6XdKPrjFdqVWhltcxJZ/yM+J/E3gY+DeAc+57M3sfUOEXKcT2HNnDpz99SmJyIjPWzuDA0QOUK1GOzg06Ex8TzzUNrqFiqYpex5RTkJ/CX8Y5t/SEX9syApRHRDy0Zd8WJidPZnLKZOZtmEdGVgbVy1XnlgtvIS4mjvZ122up40IgP4V/l5mdAzgAM+sB/BLQVCISFM45VqeuPn7nbNK2JABizozh75f+nfjYeFrUbEER01NaC5P8FP6/4VseOdbMtgLrgVsDmkpEAiYzK5NFWxYdPzm7dvdaAFrWasmIDiOIi40jNirW45QSSH9a+J1zPwNXmllZoIhzbn/gY4lIQTqcfpg56+f4Ts6mTCH1UColipagQ90OPHTpQ3SN6UqN8jW8jilBkp+reioCtwPRQLFjc/3OuX4BTSYip2X34d1MXzOdxJREZq6dyaH0Q1QoWYFrG1xLfGw8nep3okLJCl7HFA/kZ6rnU2AxsBLICmwcETkdG/dsZHKK7+Ts/A3zyXSZnFX+LO5sfCfxsfG0jW5LiaIlvI4pHstP4S/lnBsQ8CQictKcc6zcufL4nbPLty8H4Pwq5/No60eJj42n6VlNdXJWfic/hf8/ZtYHmAakHdvonNsdsFQikquMrAy+2vTV8StxNuzZgGG0qt2KkR1HEhcTR4MzG3gdU0JYfgr/UWAkMBD/JZ3+z/UCFUpEfu9Q+iFmr5tNYkoiU1Om8uvhXylZtCRX1ruSgW0Gct2511GtXDWvY0qYyE/hHwDUd87tCnQYEfmfXYd2MW3NNBKTE5m1bhaHMw5TsVRFupzbhfiYeK6ufzXlSpTzOqaEofwU/lXAoUAHERH4+befj985u2DTArJcFrUr1Kb3xb2Ji4nj8rMvp3jR4l7HlDCXn8KfCawws3n8fo5fl3OKnCbnHMu3Lz/+GMLvd3wPQKNqjRjYZiDxsfFcVP0irXQpBSo/hT/R/yEiBSA9M50FmxYcv3N2095NFLEiXFbnMl646gXiYuOoV0mn0CRw8nPn7rvBCCJSmB04eoDP1n7G5JTJTFszjd+O/EapYqW4+pyreartU3Q5twtVylbxOqZEiFwLv5lNcM7dYGYr+d/VPMc451zjwEYTCW87D+5kaspUElMSmb1uNmmZaVQuXZmuMV2Jj42nY72OlC1R1uuYEoHyGvE/4P/8I771+I8x4B8BSyQSxtbuXnv8ZqqvN3+NwxFdMZp7m91LXGwcl9W5jGJF9OA78Vau/wKdc8eWXq7vnNuY/WtmpqX7RPCdnE3alsTklMkkJieyKnUVAE2qN2FI2yHEx8bTqFojnZyVkJLXVM+9QF+gnpl9n+1L5YGvAh1MJFQdzTzK/A3zj5+c3bp/K0WtKJeffTkJTRPoGtOV6IrRXscUyVVev3O+D8wAngMey7Z9v5ZrkEizL20fM9fOJDE5kU9/+pS9aXspU7wMnep3Ij4mnmvPvZbKpSt7HVMkX/Ka6tkL7AVuCl4ckdCx/cB2pqRMITE5kTnr53A08yhRZaLofl534mPjubLelZQuXtrrmCInTWeZRLJJ2ZVyfPGzJVuW4HDUq1SP+1vcT3xsPJfWupSiRYp6HVPktKjwS0TLclks3br0+Hx98q5kAJrWaMrTVzxNfGw8Das01MlZKVRU+CXipGWkMW/DvOPFfvuB7RQrUox20e24r/l9dI3pSu0zansdUyRgAlb4zextoAuw0zl3gX9bZeBDfI9x3ADc4Jz7LVAZRI7Ze2Qvn/70KYkpicz4aQb7j+6nXIlyXFP/GuJi4ujcoDOVSlfyOqZIUARyxD8OeBV4L9u2x4A5zrkRZvaY//2jAcwgEWzrvq2+k7MpicxbP4/0rHSqla1Gzwt6Eh8bT/u67SlVrJTXMUWCLmCF3zn3pZlFn7A5Dmjnf/0u8AUq/FLAdh3aRa/JvZi6ZioADSo3oH/L/sTHxnNJrUv0GEKJeMGe46927I5g59wvZlY1tx3NLAFIAKhTp06Q4km4S9qWRPcJ3dl+YDtD2g7hxoY3EhsVq5OzItmE7Mld59wYYAxAs2bNTlwkTuR3nHOMXTaW+2bcR/Vy1Vn414U0r9nc61giISnYhX+HmdXwj/ZrADuD3L4UQofTD3Pfp/fx9oq36VivI+93f5+oMlFexxIJWcGe7JwC3OF/fQcwOcjtSyGz/rf1tH67NW+veJtBbQYx45YZKvoifyKQl3N+gO9EbpSZbQGGACOACWZ2F7AJuD5Q7UvhN+OnGdzyyS1kuSym9JzCdTHXeR1JJCwE8qqe3Nb46RCoNiUyZLkshs0fxtD5Q7mw2oV8fMPH1K9c3+tYImEjZE/uiuRk9+Hd3PrJrcxYO4PbGt3GG13eoEzxMl7HEgkrKvwSNpb/spxuE7qxdd9WXu/8Ovc0u0eXaYqcAhV+CQvvLH+He6ffS5WyVVjw1wVcUusSryOJhC0VfglpRzKO0G9GP95c9ibt67bn/7r/H1XKVvE6lkhYU+GXkLVxz0Z6TOxB0rYkHmv9GMPaD9ODykUKgH6KJCTNWjeLmz6+iYysDCbdOIn42HivI4kUGlqtSkJKlsti+JfD6fTfTpxV/iyS+iSp6IsUMI34JWTsObKH2ybdxrQ107j5wpsZ02UMZUuU9TqWSKGjwi8h4bvt39FtQjc27d3EK9e8wt+a/02XaooEiAq/eO69797j7ml3U7l0ZebfOZ9WtVt5HUmkUFPhF8+kZaTR/7P+jE4aTduz2/Jhjw+pVq6a17FECj0VfvHE5r2b6TGxB0u3LuXhVg/zbIdndammSJDoJ02Cbs7Pc+j5cU+OZBzho+s/ovv53b2OJBJRdDmnBI1zjhELR3DVf6+iSpkqfNPnGxV9EQ9oxC9BsffIXu6cfCeJyYnc2PBGxnYdS7kS5byOJRKRVPgl4FbuWEn3Cd1Zv2c9L179Ig9c8oAu1RTxkAq/BNT7K9+nz9Q+VChZgbm3z6XN2W28jiQS8TTHLwFxNPMo/Wb045ZPbqFpjaYsS1imoi8SIjTilwK3dd9Wrp94PYu2LKJ/y/48f+XzFC9a3OtYIuKnwi8F6osNX3DjRzdy8OhBPuzxITc0vMHrSCJyAk31SIFwzjHq61Fc+d6VVCpViaV9lqroi4QojfjltO1L20evyb34+MeP6X5ed96Oe5sKJSt4HUtEcqHCL6dldepqun3YjbW71zKq4ygGXDpAl2qKhDgVfjllH/7wIXdNuYuyJcry+e2f0y66ndeRRCQfNMcvJy09M53+M/vT8+OeNK7emGUJy1T0RcKIRvxyUn7Z/ws3fHQDCzct5P4W9zPqqlGUKFrC61gichJU+CXfFmxcwA0f3cC+tH2M7zaemy+82etIInIKNNUjf8o5x4uLXuSKd6+gfInyLOm9REVfJIxpxC95OnD0AHdNuYsJqyYQHxvPuLhxnFHqDK9jichpUOGXXCXvSqbbh91I+TWFER1G8EjrR3Sppkgh4MlUj5k9YGY/mNkqM3vQiwySt49Xf0zzN5uz69AuZt06i0cve1RFX6SQCHrhN7MLgD5AC6Ax0MXMGgQ7h+QsIyuDh2c9TI+JPWhYpSHfJnxLh3odvI4lIgXIixH/ecBi59wh51wGMB/4iwc55ARb9m2h4386MmrRKPo268v8O+dT+4zaXscSkQLmxRz/D8AzZnYmcBjoDCSduJOZJQAJAHXq1AlqwEiz8+BORiwcwevfvE4RK8J78e9xW+PbvI4lIgES9MLvnPvRzJ4HZgMHgO+AjBz2GwOMAWjWrJkLasgIsefIHkZ9PYqXFr/E4YzD3N74doa0HUJ0xWivo4lIAHlyVY9z7i3gLQAzexbY4kWOSHXg6AFeXvIyI78eyZ4je7ih4Q0MbTeU2KhYr6OJSBB4UvjNrKpzbqeZ1QG6AZd6kSPSHMk4whtJb/DsgmdJPZRKl3O7MOyKYTSp3sTraCISRF5dx/+xf44/Hfibc+43j3JEhPTMdN5Z8Q7DvhzGln1baF+3Pc+0f4aWtVp6HU1EPODVVI+euh0EmVmZfPDDBzz1xVOs+20dLWu15N34d2lft73X0UTEQ7pztxByzjEpeRKD5w1mdepqGldrzNSbpnJtg2t1E5aIqPAXJs45Plv3GYPmDuLbX74l5swYPuzxIT3O70ER03p8IuKjwl9IfLnxSwbOHcjCTQuJrhjNO3HvcGujWylWRN9iEfk9VYUw983Wbxg0bxCz1s2iRrkavNb5NXpf3FsPRxGRXKnwh6kfdv7A4HmDSUxO5MzSZzKy40j6Nu9LmeJlvI4mIiFOhT/MrN29liFfDOGDlR9QvmR5hrYbyoMtH6RCyQpeRxORMKHCHyY2793M0/Of5p0V71CiaAkeaf0Ij7R+hMqlK3sdTUTCjAp/iNtxYAfPLniWN759A4C+zfvyRJsnqF6uusfJRCRcqfCHqN2HdzPyq5G8vPRl0jLSuLPJnTzZ9knqnKGVSkXk9Kjwh5j9aft5afFLjFo0iv1p++l5QU+GthtKgzP1rBoRKRgq/CHicPphXv/mdUZ8NYJdh3YRFxPHsCuGcWG1C72OJiKFjAq/x45mHuWtZW8xfMFwtu3fRsd6HRnefjgtarbwOpqIFFIq/B7JzMrkv9//l6Hzh7J+z3pa127N+93ep210W6+jiUghp8IfZFkui49Xf8yTXzxJ8q5kLq5xMa91fo1O9TtpATURCQoV/iBxzvHpT58yeN5glm9fznlR5/HR9R/R7bxuKvgiElQq/EEwb/08Bs0bxNebv6ZepXq8F/8eN194M0WLFPU6mohEIBX+AFqyZQkD5w5kzvo51CxfkzeufYNeF/WieNHiXkcTkQimwh8A3+/4nkFzBzF1zVSqlKnCC1e9wL3N76VUsVJeRxMRUeEvSCm7UhjyxRA+XPUhZ5Q8g+FXDOeBlg9QrkQ5r6OJiBynwl8ANu7ZyND5Q3n3u3cpXaw0T1z2BA+1eohKpSt5HU1E5A9U+E/DL/t/4dkFz/Lvb/9NEStCvxb9eLzN41QtW9XraCIiuVLhPwW/HvqV5796nleXvkp6Vjq9mvRicNvB1KpQy+toIiJ/SoX/JOxL28eLi17kn4v+yYGjB7il0S081fYpzql8jtfRRETyTYU/Hw6lH+K1pa/x/FfP8+vhX+l2Xjeebvc0Das29DqaiMhJU+HPQ1pGGmOXjWX4guFsP7CdTvU7MfyK4TQ9q6nX0URETpkKfw4ysjL4z3f/Yej8oWzcu5HLz76ciddP5LI6l3kdTUTktKnwZ5Plspi4aiJDvhhCyq8pNDurGWOuG0PHeh21no6IFBoq/PgWUJu2ZhqD5w3mux3fcUHVC5h04yTiYuJU8EWk0In4wj/n5zkMnDuQJVuXUL9yfcZ3G8+NDW/UAmoiUmhFbOFftHkRA+cOZN6GedSuUJs3r3uTOxrfoQXURKTQi7jCv2L7CgbNHcT0n6ZTtWxV/tXpXyQ0TdACaiISMTwp/GbWH+gNOGAl8Ffn3JFAtpm8K5kn5z3JxNUTqVSqEs91eI77W9xP2RJlA9msiEjICXrhN7OaQD/gfOfcYTObAPQExgWivfW/rWfo/KH85/v/UKZ4GQZfPpgBlw6gYqmKgWhORCTkeTXVUwwobWbpQBlgWyAaGf7lcJ6e/zRFixSlf8v+PNr6UaqUrRKIpkREwkbQC79zbquZjQI2AYeBWc65WSfuZ2YJQAJAnTp1Tqmt6IrR9L64NwPbDKRmhZqnkVpEpPAw51xwGzSrBHwM3AjsASYCHznn/pvbn2nWrJlLSkoKUkIRkcLBzL51zjU7cXsRD7JcCax3zqU659KBT4BWHuQQEYlIXhT+TUBLMytjvttiOwA/epBDRCQiBb3wO+eWAB8By/BdylkEGBPsHCIikcqTq8dBENYAAAdHSURBVHqcc0OAIV60LSIS6byY6hEREQ+p8IuIRBgVfhGRCKPCLyISYYJ+A9epMLNUYKPXOf5EFLDL6xAFoLD0A9SXUFRY+gHh0ZeznXN/WKcmLAp/ODCzpJzukAs3haUfoL6EosLSDwjvvmiqR0Qkwqjwi4hEGBX+glNY7j4uLP0A9SUUFZZ+QBj3RXP8IiIRRiN+EZEIo8IvIhJhVPhPgZltMLOVZrbCzJKybb/fzFLMbJWZ/cPLjPmVU1/MrImZLT62zcxaeJ3zz5hZRTP7yMySzexHM7vUzCqb2Wwz+8n/uZLXOfMjl76M9L//3swmmVlYPDQ6p75k+9pDZubMLMrLjPmVW1/C8ece55w+TvID2ABEnbDtCuBzoKT/fVWvc55GX2YB1/hfdwa+8DpnPvrxLtDb/7oEUBH4B/CYf9tjwPNe5zyNvlwFFPNvez6c++J/XRv4DN+NmVFe5SuA70tY/txrxF9w7gVGOOfSAJxzOz3OczocUMH/+gxgm4dZ/pSZVQAuB94CcM4ddc7tAeLw/bDi/xzvTcL8y60vzrlZzrkM/26LgVpeZcyvPL4vAC8Cj+D7txby8uhLWP7cq/CfGgfMMrNv/Q+FBzgXaGNmS8xsvpk19zDfycipLw8CI81sMzAKeNyzdPlTD0gF3jGz5WY21szKAtWcc78A+D9X9TJkPuXWl+x6ATOCH+2k5dgXM+sKbHXOfedxvpOR2/clLH/uVfhPTWvn3MXANcDfzOxyfA+1qQS0BB4GJvgfLRnqcurLvUB/51xtoD/+UU4IKwZcDIx2zl0EHMQ3tROO8uyLmQ0EMoDx3sQ7KTn15SlgIPCkh7lORW7fl7D8uVfhPwXOuW3+zzuBSUALYAvwifNZCmThW8QppOXSlzuAT/y7TPRvC2VbgC3O91hP8D3a82Jgh5nVAPB/Dodfw3PrC2Z2B9AFuMX5J5RDXG59qQt8Z2Yb8E1ZLTOz6t5EzLfc+hKWP/cq/CfJ/6tq+WOv8Z10+wFIBNr7t5+L7+RPSK/cl0dftgFt/bu1B37yJmH+OOe2A5vNLMa/qQOwGpiC7z8x/J8nexDvpOTWFzPrBDwKdHXOHfIs4EnIpS/LnHNVnXPRzrlofIXzYv++ISuPf2Nh93MPHj1zN8xVAyb5f5srBrzvnJtpZiWAt83sB+AocEcYjMpy68sB4F9mVgw4AiTkcYxQcT8w3v99+Bn4K76BzQQzuwvYBFzvYb6TkVNfvgFKArP936/Fzrl7vIuYbzn1JVzl1JeDhN/PvZZsEBGJNJrqERGJMCr8IiIRRoVfRCTCqPCLiEQYFX4RkQijwi8hz78qYl+vc+SHmT1oZmUCePwaZjbN/7qdme31r6K6wsw+92+/z8zC+bJJCTAVfgkHFYGQKPzmk9fPzYPASRV+//0S+TUAeDPb+wXOuSb+jyv9294G+p1MBoksKvwSDkYA5/hHtSMBzOxhM/vGvz79UP+2aP9a6WPN7AczG29mV5rZV+Zbk7+Ff7+nzOw/ZjbXv73PsYbyOO6PZvY6sAyobWajzfesglXZ9usHnAXMM7N5/m0Hsh27h5mN878eZ2Yv+Pd73n8X9dv+tpebWVwufxfdgZl5/WX57+zdYGHwHAXxhu7clXDwGHCBc64JgJldBTTAt4aQAVP8i8ttAurju0M3Ad/drjcDlwFdgSf439LMjfAtrFUWWG5m04EL8jhuDPBX51xff4aBzrndZlYUmGNmjZxzL5vZAOAK51x+bts/F7jSOZdpZs8Cc51zvcz3kJWlZva5c+7gsZ3NrC7w27ElgP3amNkK/+uJzrln/K+TgDbA0nzkkAijwi/h6Cr/x3L/+3L4CvYmYL1zbiWAma0C5jjnnJmtBKKzHWOyc+4wcNg/6m6B7z+I3I670Tm3ONufv8F8y1gXA2oA5wPfn2Q/JjrnMrP1qauZPeR/XwqoA/yYbf8a+JYGzm6Bc65LDsfeCcSeZB6JECr8Eo4MeM459+/fbTSLBrKPhrOyvc/i9//eT1yrxP3JcU8ceT8ENHfO/eafvimVS9bs7Zy4z8Fsrw3o7pxLyeU4AIfzaOdEpfz7i/yB5vglHOwHymd7/xnQy8zKAZhZTTM72YesxJlZKTM7E2iHb1oov8etgK9o7zWzavieZZBb1h1mdp7/hPBf8sjzGXC/+VdgM7OLcthnDb//rSUv5+JbaVXkDzTil5DnnPvVf4L2B2CGc+5hMzsPWOSvkweAW4HMvI5zgqXAdHzTKcP8zyXYlp/jOue+M7PlwCp8qzR+le3LY4AZZvaLc+4KfOcnpgGb8RXicrnkGQa8BHzvL/4b8K29n73dg2a2zszqO+fW/kn/WgND/2QfiVBanVMijpk9BRxwzo3yOsvJMrO/AE2dc4Py2OciYIBz7rbgJZNwohG/SBhxzk3yT0/lJQoYHIw8Ep404hcRiTA6uSsiEmFU+EVEIowKv4hIhFHhFxGJMCr8IiIR5v8B7cysBYtrpX0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.plot(temperature, clock, color = 'green', linestyle= '-' )\n",
    "plt.title('Ta vs time')\n",
    "plt.xlabel('temperature (F)')\n",
    "plt.ylabel('time');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1445,
   "metadata": {},
   "outputs": [],
   "source": [
    "'this is the nonlinear model'\n",
    "tp= np.linspace(65,67,100) #the ranging ambient temperature\n",
    "\n",
    "tp2=np.empty(len(tp)) #tp2 is then new temperature for the changing ambient temperature\n",
    "for i in range(0,len(tp)-1):\n",
    "    tp2[0]=85\n",
    "    tp2[i+1]= tp2[i]-(tp[i]-tp2[i])*-kconst\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1454,
   "metadata": {},
   "outputs": [],
   "source": [
    "'this is the linear analytical model'\n",
    "timestp=np.linspace(0,6,100)\n",
    "tnum3=np.empty(len(timestp))\n",
    "for i in range(0,len(timestp)-1):\n",
    "    tnum3[0]=85\n",
    "    tnum3[i+1]= tnum3[i] - (t_amb-tnum3[i])*-kconst\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1455,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU5fX48c/JRhJICEu0JBQBRaiswYACCiLK4oJIVaRqAauoX5fa/kSh/ip0+0rF1q1WxK36E3EFXKACCrgWWUWUpWxRkiBEMGELWc/vj3tnmIRJMklmMlnO+/XKKzN3PTPiPbnP89zziKpijDHGlBcR7gCMMcbUT5YgjDHG+GUJwhhjjF+WIIwxxvhlCcIYY4xfliCMMcb4ZQnCGAOAiJwpIrnhjsPUH5YgTJ0QkSM+P6Uiku/z/rpwx1cbIvK9iJwX7jiqq3zcqvpfVU0KZ0ymfokKdwCmaVDVFp7XIpIB3KSqH4QvosCISJSqFjf0cxhTE3YHYeoFEYkUkd+LyC4R+UFE5opIkruum4gUi8ivRCRLRA6IyI0iMkBEvhaRXBH5u8+xbhWR5SLytIgcEpHNIjLYZ31rEXnJ/Qt6j4hMF5GIcvs+KSI/AlPd868UkYMikiMiL4pIgrv9G8ApwFL3buguERkpIjvKfT7vX+siMlNEXhGR10TkMHBtZZ/fz3e1S0Qu8nnfTETyROQsEWkuIq+6seaKyBci0srPMfzF3U1Ein22WSUiM0RktbvNfBFpIyKvu9/rKhFp77N9D/e7+1FEtojImOr8GzD1jyUIU19MAYYD5wHtgSLgEZ/1kUAvoDMwCXgCuAcY4i6fJCLn+Gw/GNgItAFmAgtFJNFdNxfIc4/VHxgD3FBu3y+BtsDf3GV/BH4C9AS6AvcDqOrVwH5guKq2UNXHA/y8PwdeBFoCbwXw+X29Coz3eX8ZsFtVNwM34bQMpLrx3wEUlj9ANeIeB1wDdMD57J8BTwKtgW9xvwf3u10GPOee95fA8yJyRhXfg6nHLEGY+uIWYKqqZqvqceAPwDgREZ9t/qiqBar6jvv+JVU9oKrfAZ8DaT7b7lHVf6pqkaq+BGQCI0TkNJwE8FtVPaaqe4HHgWt99t2lqs+oaomq5qvqVlVdrqqFqvo98ChOYqqNj1R1saqWqmp+gJ/fYy4wVkRi3Pe/AF5xXxcBycDpqlqsqmtU9Wgt4nxWVTNU9SCwFNiiqh+5TWJvcuI7vxL4WlXnut/bGuBdnERoGijrgzBh514EfwosFhHf6pEROHcAACWqesBnXT6wr9z7Fj7vM8ud5lsgBTgNiAVyfK69EYBvk9CecvGlAI8BA4EEd/u9gXy2SnjPEcDn/8F3R1X9RkT2AKNEZAUwCrjbXf0czp3OmyLSAngJ+L2qltQwzvLfcUXf+WnA4HKjoKKAH2t4XlMPWIIwYaeqKiJZwFhVXVd+vYi0rcFh25d73wHIxrkwHwFaacWljMsvnwUcBXqo6o8ici3w50q2PwrEe96ISDROk4zfc1T1+SswD6eZqTWwRlX3uMcqAB4AHhCRzsAS4Bucu47yglnKeQ+wVFUvD+IxTZhZE5OpL2YDM0XkpwAicoqI1OZi81O3wzlKRK7HSRBLVXU3sAp4SEQSRCRCRLpI5cNUE3CSyiER6QD8ttz6fTj9GR5bgNYiMsxNDn+g6v/Xqvv55+H0PdzEieYlROQit7M6AjgEFAMV3T2Uj7s2FgJpIjJORKJFJEZEzhWRM4N0fBMGliBMffEQ8AGw3B3Z8znQtxbH+xinffwgTkfqlaqa564bDyQBW931rwGnVnKsB3A6j/OABTidyr7+AvzFHTV0h6r+APwa56/2TOB7yjUT+VGtz6+qGTid8P2AN3xWpQJvA4eBr4HFwOsVHKZM3FXEVylV/REYgTOAYC/O3dqfgejaHNeEl9iEQaaxEZFbgatU9aIqNzbGVMjuIIwxxvhlCcIYY4xf1sRkjDHGL7uDMMYY41ejeg6ibdu22rFjx3CHYYwxDca6det+UNVkf+saVYLo2LEja9euDXcYxhjTYIjItxWtsyYmY4wxflmCMMYY45clCGOMMX41qj4IYxqioqIiMjMzOX78eLhDMY1YbGws7du3Jzo68OonliCMCbPMzEwSEhLo2LEj/qd/MKZ2VJUDBw6QmZlJp06dAt4vpE1MIvIbEflGnGkh54lIrDuFYZaIfOn+XFLBviNFZJuI7BCRqaGKceGGLAbNXE6nqYsYNHM5CzdkhepUxvh1/Phx2rRpY8nBhIyI0KZNm2rfpYbsDkJEUoG7gLNUNV9EXufErF2PqOrDlewbiTOt4cU41TDXiMg77pSKQbNwQxbT5m8iv8iphpyVm8+0+ZsAGJOWGsxTGVMpSw4m1GrybyzUndRRQJyIROFMoJId4H79gR2quktVC3Hm4L0i2MHNWrKN/KIS7ohcwOCIjQDkF5Uwa8m2YJ/KGGManJAlCFXNAh4GvsOpD5+nqkvd1XeIyFci8ryItPKzeyplp33MdJedREQmi8haEVmbk5NTrRizc/MBuCXqPYZEfHXScmNM8Dz66KMcO3asxvuvXLmSzz//PIgRnXDBBRf4fcj2nXfeYebMmTU6Zm5uLv/85z+rva4+CVmCcC/8VwCdcOYCbu7O7PUUcDrQBydx/M3f7n6W+a0qqKpzVDVdVdOTk/0+LV6hlKQ4AA4TRwvyT1puTH3UUPvN6nOCqMjo0aOZOrVmXaCWICp3EbBbVXNUtQiYDwxU1X2qWqKqpcAzOM1J5WXiTOLu0Z7Am6cCNmVEV+KiIzmicSSI8w83LjqSKSO6BvtUxgSFp98sKzcf5US/WbCSxPvvv0/fvn3p3bs3w4YNA+DgwYOMGTOGXr16ce655/LVV87d9owZM7jxxhu54IIL6Ny5M48//jgAR48e5dJLL6V379706NGD1157jccff5zs7GyGDh3K0KFDAbjttttIT0+ne/fuTJ8+3RtDx44dmT59On379qVnz55s3bqVjIwMZs+ezSOPPEKfPn345JNPysS9evVqBg4cSFpaGgMHDmTbNqeZ+F//+hdjxozh8ssvp1OnTvzjH//g73//O2lpaZx77rkcPHjQe4yXX36ZgQMH0qNHD1avXu3d/447nMn2cnJy+PnPf06/fv3o168fn332WaXfw9SpU9m5cyd9+vRhypQpZeL1t27WrFn069ePXr16eb+PjIwMunXrxk033USPHj247rrr+OCDDxg0aBBdunTxxjljxgxuuOEGLrzwQrp06cIzzzxT838EvlQ1JD/AOTiTpcfj3BG8CNwJtPPZ5jfAq372jQJ24dx9xOBMrdi9qnOeffbZWl0L1mfqVzP66cf/d5AOfPBDXbA+s9rHMKY2Nm/eHPC2Ax/8UE+7772TfgY++GGt49i/f7+2b99ed+3apaqqBw4cUFXVO+64Q2fMmKGqqh9++KH27t1bVVWnT5+uAwYM0OPHj2tOTo62bt1aCwsL9c0339SbbrrJe9zc3FxVVT3ttNM0JyfHu9xz/OLiYh0yZIhu3LjRu93jjz+uqqpPPvmk/upXv/Keb9asWX5jz8vL06KiIlVVXbZsmY4dO1ZVVV944QU9/fTT9dChQ7p//35NTEzUp556SlVV7777bn3kkUdUVXXIkCHemD/66CPt3r27d//bb79dVVXHjx+vn3zyiaqqfvvtt9qtW7dKv4fdu3d7j1Ne+XVLlizRm2++WUtLS7WkpEQvvfRS/eijj3T37t0aGRmpX331lZaUlGjfvn110qRJWlpaqgsXLtQrrrjCG0OvXr302LFjmpOTo+3bt9esrKyTzuvv3xqwViu4poZsFJOqfiEibwLrcSZO3wDMAZ4VkT44TUYZwC0AIpICPKuql6hqsTtH7hIgEnheVb8JRZxj0lLh6/ZwPJfPbr4wFKcwJmgq6h8LRr/ZqlWrGDx4sHecfOvWrQH49NNPeestZxruCy+8kAMHDpCX50zvfemll9KsWTOaNWvGKaecwr59++jZsyf33HMP9913H5dddhnnn3++3/O9/vrrzJkzh+LiYvbu3cvmzZvp1asXAGPHjgXg7LPPZv78+VXGnpeXx4QJE9i+fTsiQlFRkXfd0KFDSUhIICEhgZYtW3L55ZcD0LNnT+/dEMD48eMBGDx4MIcOHSI3N7fMOT744AM2bz4xkPLQoUMcPny4wu+hOpYuXcrSpUtJS0sD4MiRI2zfvp0OHTrQqVMnevbsCUD37t0ZNmwYIkLPnj3JyMjwHuOKK64gLi6OuLg4hg4dyurVqxkzZky14igvpA/Kqep0YHq5xTdUsG02cInP+8U4E66HXmwi5O2pejtjwiwlKY4sP8kgGP1mqup3KKT6mVTMs12zZs28yyIjIykuLubMM89k3bp1LF68mGnTpjF8+HAeeOCBMvvv3r2bhx9+mDVr1tCqVSsmTpxYZoy+57ieY1bl97//PUOHDmXBggVkZGRwwQUXnHQsgIiICO/7iIiIMscu/9nLvy8tLeU///kPcXEnf9f+vofqUFWmTZvGLbfcUmZ5RkZG0OKvCavFBNAsAY4fCncUxlTJ02/mK1j9ZgMGDOCjjz5i9+7dAN72+cGDBzN37lzA6Shu27YtiYmJFR4nOzub+Ph4rr/+eu655x7Wr18PQEJCgvcv7kOHDtG8eXNatmzJvn37+Pe//11lfL77l5eXl0dqqjPQ8V//+ldgH7ic1157DXDumFq2bEnLli3LrB8+fDj/+Mc/vO+//PLLGsdbft2IESN4/vnnOXLkCABZWVns37+/WvG//fbbHD9+nAMHDrBy5Ur69etXrf39sVIbAM0SocD/f0hj6hPPA5yzlmwjOzeflKQ4pozoGpQHO5OTk5kzZw5jx46ltLSUU045hWXLljFjxgwmTZpEr169iI+P58UXX6z0OJs2bWLKlClEREQQHR3NU089BcDkyZMZNWoU7dq1Y8WKFaSlpdG9e3c6d+7MoEGDqozv8ssv56qrruLtt9/miSeeKNN0de+99zJhwgT+/ve/c+GFNWsqbtWqFQMHDuTQoUM8//zzJ61//PHHuf322+nVqxfFxcUMHjyY2bNnV3i8Nm3aMGjQIHr06MGoUaOYNWtWpeu2bNnCgAEDAGjRogUvv/wykZGRFR3+JP379+fSSy/lu+++4/e//z0pKSnV+PT+Nao5qdPT07VGEwat/Cus/F/4/QGItJxp6taWLVv42c9+Fu4wTAM2Y8YMWrRowT333FPpdv7+rYnIOlVN97e9NTGB08QEUGh3EcYY42F/LoPTSQ1OP0Scvwe7jTGm/poxY0ZIjmt3EHDiDsL6IYwxxssSBDid1AAFNpLJGGM8LEGAT4KwOwhjjPGwBAFl+yCMMcYAliAc3j4ISxDGAEycOJE333wzpOdYu3Ytd911V0jPEaiOHTvyww8/nLR89uzZvPTSSzU6ZkZGBq+88kq119UnliDA+iCMCYP09HRv5dP66tZbb+WXv/xljfa1BNFYRMeBRFofhGmSXnrpJXr16kXv3r254YYTpdI+/vhjBg4cSOfOnb13E0eOHGHYsGHeUtxvv/024Fzwfvazn3HzzTfTvXt3hg8fTn6+UzNqzZo19OrViwEDBjBlyhR69OgBOGU7LrvsMqDiktkAf/rTn+jWrRsXX3wx48eP5+GHT56t+N133+Wcc84hLS2Niy66yFssb8aMGUyYMIHhw4fTsWNH5s+fz7333kvPnj0ZOXJkmaJ+s2bNon///vTv358dO3Z49/ecb+fOnYwcOZKzzz6b888/n61btwLO3dZdd9110nc1depUPvnkE/r06cMjjzxSJt7y60pKSpgyZYq33PfTTz/t/Y6GDBnCNddcw5lnnsnUqVOZO3cu/fv3p2fPnuzcudMbw6233sr555/PmWeeyXvvvVfNfwX+2XMQACJOP4T1QZhw+/dU+H5TcI/5k54wyv+saN988w1/+ctf+Oyzz2jbtm2Z+RH27t3Lp59+ytatWxk9ejRXXXUVsbGxLFiwgMTERH744QfOPfdcRo8eDcD27duZN28ezzzzDNdccw1vvfUW119/PZMmTWLOnDkMHDiw0sl3tm7dyooVKzh8+DBdu3bltttuY+PGjbz11lts2LCB4uJi+vbty9lnn33Svueddx6rVq1CRHj22Wd56KGH+NvfnLnIdu7cyYoVK9i8eTMDBgzgrbfe4qGHHuLKK69k0aJF3oqniYmJrF69mpdeeom77777pIvs5MmTmT17Nl26dOGLL77gf/7nf1i+fHmF39XMmTN5+OGH/V6sy6+bM2cOLVu2ZM2aNRQUFDBo0CCGDx8OwMaNG9myZQutW7emc+fO3HTTTaxevZrHHnuMJ554gkcffRRwkvRHH33Ezp07GTp0KDt27CA2NrbC7zsQliA8miXYHYRpcpYvX85VV11F27ZtgRMlvgHGjBlDREQEZ511lvcvclXld7/7HR9//DERERFkZWV513Xq1Ik+ffoATpnujIwMcnNzOXz4MAMHDgTgF7/4RYV/3formf3pp596y1gD3lLd5WVmZjJu3Dj27t1LYWGht2Q5wKhRo4iOjqZnz56UlJQwcuRIgJPKZXvKfY8fP57f/OY3ZY5/5MgRPv/8c66++mrvsoKCgkq/q+pYunQpX331lffuIy8vj+3btxMTE0O/fv1o164dAKeffro3cfTs2ZMVK1Z4j3HNNdcQERFBly5d6Ny5M1u3bvX+96gpSxAezVpaH4QJvwr+0g+Vikp8Q9kS1p6abXPnziUnJ4d169YRHR1Nx44dvWW6y5e8zs/P91sqvCL+SmYHuv+dd97Jb3/7W0aPHs3KlSvLPFnsWx47Ojra+3krK5ftr9R3UlJShRVc/X1X1aGqPPHEE4wYMaLM8pUrV1q573rB7iBMEzRs2DBef/11Dhw4AFCmicmfvLw8TjnlFKKjo1mxYgXffvttpdu3atWKhIQEVq1aBcCrr75arfjOO+883n33XY4fP86RI0dYtGhRhXF5yn1XVW22Ip5y36+99pq3qqpHYmIinTp14o033gCcC/rGjRsrPV51y30/9dRT3j6R//73vxw9erRa8b/xxhuUlpayc+dOdu3aRdeutS8BH9I7CBH5DXATzuxxm4BJwJ+Ay4FCYCcwSVVz/eybARwGSoDiiqoNBk1sIhwK+rTXxtRr3bt35/7772fIkCFERkaSlpZW6XwK1113HZdffjnp6en06dOHbt26VXmO5557jptvvpnmzZtzwQUXnDTPQmX69evH6NGj6d27N6eddhrp6el+958xYwZXX301qampnHvuud45LaqjoKCAc845h9LSUubNm3fS+rlz53Lbbbfx5z//maKiIq699lp69+5d4fF69epFVFQUvXv3ZuLEiWWarcqv+/Wvf01GRgZ9+/ZFVUlOTmbhwoXVir9r164MGTKEffv2MXv27Fr3P0AIy32LSCrwKXCWquaLyOs4M8RlA8vdaUX/CqCq9/nZPwNIV9WTBydXoMblvgHeugky18KvK58ExJhga+zlvo8cOUKLFi0Ap3N27969PPbYY9Xe/9ixYwwePJg5c+bQt2/fUIXbIE2cOJHLLruMq666qtLtqlvuO9R9EFFAnIgUAfFAtqou9Vm/Cqj8E9UVmzTImJBYtGgRDz74IMXFxZx22mnVnvFt8uTJbN68mePHjzNhwgRLDnUoZAlCVbNE5GHgOyAfWFouOQDcCLxW0SGApSKiwNOqOsffRiIyGZgM0KFDh5oH3CzBOqmNCYFx48Yxbty4Gu/fEB4oC7eaTrNalZB1UotIK+AKoBOQAjQXket91t8PFANzKzjEIFXtC4wCbheRwf42UtU5qpququnJyck1Dzg2EUoKobig6m2NCbLGNLOjqZ9q8m8slKOYLgJ2q2qOqhYB84GBACIyAbgMuE4riFpVs93f+4EFQP8Qxnqi3IY9LGfqWGxsLAcOHLAkYUJGVTlw4EC1O65D2QfxHXCuiMTjNDENA9aKyEjgPmCIqh7zt6OINAciVPWw+3o48McQxlq2HlOLWtyJGFNN7du3JzMzk5ycnHCHYhqx2NhY2rdvX619QtkH8YWIvAmsx2lK2gDMAb4BmgHL3Ac5VqnqrSKSAjyrqpcApwIL3PVRwCuq+n6oYgWsoqsJm+jo6DJP/hpTX4R0FJOqTgeml1t8RgXbZgOXuK93ARUPMA6FWJs0yBhjfNmT1B6eOwjrgzDGGMASxAk27agxxpRhCcLDJg0yxpgyLEF4WCe1McaUYQnCIyoGomKtD8IYY1yWIHxZyW9jjPGyBOGrWaI1MRljjMsShC+7gzDGGC9LEL5iE60PwhhjXJYgfNmcEMYY42UJwpf1QRhjjJclCF82aZAxxnhZgvAV6zYxWV1+Y4yxBFFGswTQUig8Gu5IjDEm7CxB+LJ6TMYY42UJwpe3HpONZDLGmJAmCBH5jYh8IyJfi8g8EYkVkdYiskxEtru/W1Ww70gR2SYiO0Rkaijj9Ipt6fy2ZyGMMSZ0CUJEUoG7gHRV7QFEAtcCU4EPVbUL8KH7vvy+kcCTwCjgLGC8iJwVqli9rKKrMcZ4hbqJKQqIE5EoIB7IBq4AXnTXvwiM8bNff2CHqu5S1ULgVXe/0LI+CGOM8QpZglDVLOBh4DtgL5CnqkuBU1V1r7vNXuAUP7unAnt83me6y04iIpNFZK2IrM3Jyald0NYHYYwxXqFsYmqF81d/JyAFaC4i1we6u59lfh9OUNU5qpququnJyck1C9Yj1r2DsD4IY4wJaRPTRcBuVc1R1SJgPjAQ2Cci7QDc3/v97JsJ/NTnfXuc5qmQWrg5D4BHF69j0MzlLNyQFepTGmNMvRXKBPEdcK6IxIuIAMOALcA7wAR3mwnA2372XQN0EZFOIhKD07n9TghjZeGGLKYt2MxhjaMF+WTl5jNt/iZLEsaYJiuUfRBfAG8C64FN7rnmADOBi0VkO3Cx+x4RSRGRxe6+xcAdwBKcpPK6qn4TqlgBZi3ZRn5RCUeII5FjAOQXlTBrybZQntYYY+qtqFAeXFWnA9PLLS7AuZsov202cInP+8XA4lDG5ys7Nx+AXG1Bkhw5abkxxjQ19iS1KyUpDoCDmkArOXzScmOMaWosQbimjOhKXHQkP5JAa5wEERcdyZQRXcMcmTHGhEdIm5gakjFpzmMWBe8l0ar4a1KT4pgyoqt3uTHGNDWWIHyMSUuF3N7w0RI+u3cIRESGOyRjjAkba2IqL74NoJCfG+5IjDEmrCxBlBffxvl97EB44zDGmDCrsolJRNKB83HKZeQDXwMfqOrBEMcWHvGtnd+WIIwxTVyFdxAiMlFE1gPTgDhgG05ZjPOAZSLyooh0qJsw65DdQRhjDFD5HURzYJCq+n1STET6AF1wSmo0HpYgjDEGqCRBqOqTle2oql8GP5x6IM6amIwxBipvYlrq83pa3YRTD8TEQ3S8JQhjTJNX2Sgm38kVrg51IPVKfBs41jj74I0xJlCVJQi/E/Q0CfGt7Q7CGNPkVdZJ3VlE3sGZ3c3z2ktVR4c0snCKb2MJwhjT5FWWIK7wef1wqAOpV+LbwMHd4Y7CGGPCqrJRTB/V5sAi0hV4zWdRZ+ABYADgKZGaBOSqah8/+2cAh4ESoFhV02sTT7VYH4QxxlScIETkXZwZ4N5355T2XdcZmAhkqOrz/vZX1W1AH3f7SCALWKCqj/oc529AXiXxDVXVHwL7KEEU3wYK8qCkCCKj6/z0xhhTH1TWxHQz8FvgURE5COQAsUBHYCfwD1X1N5+0P8OAnar6rWeBO0/1NcCFNYg7tLzlNg5CwqnhjcUYY8Kksiam74F7gXtFpCPQDqcW039V9Vg1z3MtMK/csvOBfaq6vaIQgKUiosDTqjrH30YiMhmYDNChQ5Aqf/g+TW0JwhjTRAU0H4SqZgAZNTmBiMQAo3FqOvkaz8lJw9cgVc0WkVNwaj9tVdWP/cQ2B6cpjPT09OAMzbVyG8YYUyflvkcB61V1n2eBiEQBYynbiV2Gqma7v/cDC4D+IY7zBEsQxhhTJwnC353CRcBWVc30t4OINBeRBM9rYDhOmfG6YQnCGGMCSxAiEucOW60WEYkHLgbml1t1Up+EiKSIyGL37anApyKyEVgNLFLV96t7/hqL8+mkNsaYJiqQCYMux3lQLgbo5Jb5/mMgT1K7ndlt/Cyf6GdZNnCJ+3oX0Luq44dMVAw0S7Q7CGNMkxbIHcQMnPb/XPCW+e4YupDqCavHZIxp4gJJEMWqWtnDbI2T1WMyxjRxgSSIr0XkF0CkiHQRkSeAz0McV/hZgjDGNHGBJIg7ge5AAfAKTmmMu0MZVL1g9ZiMMU1cpZ3Ubg2lP6jqFOD+ugmpnrA7CGNME1fpHYSqlgBn11Es9Ut8ayg6CkX54Y7EGGPCIpBSGxvcyYLeAI56Fqpq+WcbGhfvw3IHoWVqeGMxxpgwCCRBtAYOULbqqnLyw2+Ni+/T1JYgjDFNUJUJQlUn1UUg9Y6V2zDGNHGBPEn9As4dQxmqemNIIqovLEEYY5q4QJqY3vN5HQtcCWSHJpx6xLcPwhhjmqAqn4NQ1bd8fubizALXI/ShhdfbW49SivDou/9h0MzlLNyQFe6QjDGmTtWk3HcXIEhTt9VPCzdkMXXhFvK0Oa04TFZuPtPmb7IkYYxpUqpMECJyWEQOeX6Ad4H7Qh9a+Mxaso38ohIOagJt5DAA+UUlzFqyLcyRGWNM3QlkFFNCXQRSn2TnOg/H/UBLkiX3pOXGGNMUBHIH8WEgyxqTlKQ4AL7X1pzKjyctN8aYpqDCBCEisSLSGmgrIq1EpLX70xFIqerAItJVRL70+TkkIneLyAwRyfJZfkkF+48UkW0iskNEptb0A9bElBFdiYuOZK+2pp0cBJS46EimjKj2pHrGGNNgVdbEdAtO1dYUYB0g7vJDwJNVHVhVtwF9wFv0LwtYAEwCHlHVhyva193+SZzpSjOBNSLyjqpuruq8wTAmzXlyOmPRKTQrLuKslsVMHpnmXW6MMU1BhQlCVR8DHhORO1X1iVqeZxiwU1W/FZEqN8aZwW6HO/UoIvIqcAVQJwkC3CTR7N3uK0gAABnJSURBVAJ4/XkWTzodfmLJwRjTtATSSf2EiPQAzsJ5UM6z/KVqnOdaYJ7P+ztE5JfAWuD/qOqP5bZPBfb4vM8EzvF3YBGZDEwG6NAhyKNvE92kcCgbftIzuMc2xph6LpBO6unAE+7PUOAhYHSgJxCRGHf7N9xFTwGn4zQ/7QX+5m83P8tOKvcBoKpzVDVdVdOTk5MDDSswCe2c34ca/4PjxhhTXiAPyl2F00T0vVu4rzfQrBrnGAWsV9V9AKq6T1VLVLUUeAanOam8TOCnPu/bE47yHi1OBYmwBGGMaZICSRD57sW8WEQSgf1A52qcYzw+zUsi0s5n3ZXA1372WQN0EZFO7h3ItcA71ThncERGQYufWIIwxjRJgRTrWysiSTh/7a8DjgCrAzm4iMTjjES6xWfxQyLSB6fJKMOzTkRSgGdV9RJVLRaRO4AlQCTwvKp+E9hHCrLEdnDYEoQxpumpak5qAR5U1Vxgtoi8DySq6leBHFxVjwFtyi27oYJts4FLfN4vBhYHcp6QSkyBH7aHOwpjjKlzVc1JrcBCn/cZgSaHRiMhxZqYjDFNUiB9EKtEpF/II6mvElOg4BAUHA53JMYYU6cC6YMYCtwqIhnAUZwhqKqqvUIZWL3hfRZiLyQ3ubqFxpgmLJAEMSrkUdRniZ5nIbIg+czwxmKMMXUokBnlvsV5JuFC9/WxQPZrNBLduoSH94Y3DmOMqWOBPkl9HzDNXRQNvBzKoOqVBDdBHLLZ5IwxTUsgdwJX4pTKOAre4ahNpzE+OhbiWttIJmNMkxNIgih0h7sqgIg0D21I9VBiqtNJbYwxTUggCeJ1EXkaSBKRm4EPcJ6qbjoSU6yJyRjT5ARS7vthEbkYZ6KgM4EHVHVZyCOrTxLbQda6cEdhjDF1KpBhrgCbgDicZqZNoQunnkpMhWM/QHEBRFWnkK0xxjRcgYxiugmnON9YnNLfq0TkxlAHVq945oWwoa7GmCYkkDuIKUCaqh4AEJE2wOfA86EMrF7xPAtxKBtadQxrKMYYU1cC6aTOBHwLER2m7HSgjZ/v1KPGGNNEBHIHkQV8ISJv4/RBXAGsFpHfAqjq30MYX/2QaFOPGmOankDuIHbilPz2zAn9Ns5c0gk0kQfmFm4+zDFief7fnzFo5nIWbrAhr8aYxi+QYa5/qMmBRaQr8JrPos7AA0AqcDlQiJN8JrkTEpXfPwOnOasEKFbV9JrEUVsLN2QxbcHX9JDWnCoHycrNZ9p8ZyDXmLTUcIRkjDF1IpBRTOkiskBE1ovIV56fqvZT1W2q2kdV+wBn4xT5WwAsA3q45cL/y4kaT/4MdY8RluQAMGvJNvKLStirrUmVHwDILyph1pJt4QrJGGPqRCB9EHNxRjJtAkpreJ5hwE63Guy3PstX4Qydrbeyc/MByNCfMDric5yWNvEuN8aYxiqQPogcVX1HVXer6reen2qe51pgnp/lNwL/rmAfBZaKyDoRmVzRgUVksoisFZG1OTk51QyrailJcQDs1na0lGO0dgd0eZYbY0xjFUiCmC4iz4rIeBEZ6/kJ9AQiEoNTDfaNcsvvB4px7lD8GaSqfXEmLLpdRAb720hV56hquqqmJycnBxpWwKaM6EpcdCS71BnJ1En2EhcdyZQRXYN+LmOMqU8CaWKaBHTDmQfC08SkwPwAzzEKWK+q+zwLRGQCcBkwzK0UexK3rDiqul9EFgD9gY8DPGfQeDqiX/n3ASiEvs0PcMOlPa2D2hjT6AWSIHqras9anGM8Ps1LIjISZwKiIap6zN8ObknxCFU97L4eDvyxFjHUypi0VMb0Ggd/+TX3nxMFlhyMMU1AIE1Mq0TkrJocXETigYspe7fxD5znJ5aJyJciMtvdNkVEFrvbnAp8KiIbcepALVLV92sSQ9BERkHrTnBgR1jDMMaYuhLIHcR5wAQR2Q0UAAKoO0y1Uu4dQptyy86oYNts4BL39S6gdwCx1a02Z8CBneGOwhhj6kQgCWJkyKNoKNqcATs+hNISiIgMdzTGGBNSVTYxuUNafwpc6L4+Fsh+jVKbM6CkAPIywx2JMcaEXCBPUk/H6VT2PPEcDbwcyqDqrTZu65j1QxhjmoBA7gSuxHmO4Sh4+wqaRJG+k1iCMMY0IYEkiEL3WQUF7xDUpqnFKdAs0RKEMaZJCCRBvC4iTwNJInIz8AHwbGjDqqdEoM3pliCMMU1CIOW+HxaRi4FDQFfgAVVdFvLI6qs2Z8CeL8IdhTHGhFwgndR/VdVlqjpFVe9R1WUi8te6CK5eanMG5O6BouPhjsQYY0IqkCami/0sGxXsQBqMNmcACgd3hTsSY4wJqQoThIjcJiKbgK6+EwW5T1RXOWFQo2UjmYwxTURlfRCv4MzV8CAw1Wf5YVU9GNKo6rH3MuO5DPjr3EW8k9CCKSO6WmVXY0yjVGGCUNU8IA+nGqvBnZ/63V30i0iik+y1+amNMY1a0yyZUUOe+al3azs6R+wFbH5qY0zjZQmiGjzzUP+3tD1dZQ/izp9k81MbYxojSxDV4JmH+ivtTILk01n2lllujDGNScgShIh0dScE8vwcEpG7RaS1iCwTke3u71YV7D9SRLaJyA4Rmepvm7rmmZ96Y+npAPSWnTY/tTGm0QpZglDVbaraR1X7AGfjlAlfgDMi6kNV7QJ8SNkRUgCISCTwJM7zFmcB42s6q10wjUlL5cGxPTme2Jmj2oyBcd/x4Fibn9oY0zgFMmFQMAwDdqrqtyJyBXCBu/xFYCVOOXFf/YEd7sxyiMirwBXA5jqJthJj0lKdhPBCOlcV77f5qY0xjVZd9UFcC8xzX5+qqnsB3N+n+Nk+Fdjj8z7TXXYSEZksImtFZG1OTk4QQ65CShp8vwmKC+vunMYYU4dCniBEJAZnPok3qrObn2Xqb0NVnaOq6aqanpycXJMQaya1rzO73P6w39QYY0xI1MUdxChgvaruc9/vE5F2AO7v/X72ycSZ5tSjPZAd0iirK6Wv8zt7fXjjMMaYEKmLBDGeE81LAO8AE9zXE4C3/eyzBugiIp3cO5Br3f3qj1YdIa41ZG8IdyTGGBMSIU0QIhKPUw12vs/imcDFIrLdXTfT3TZFRBYDqGoxcAewBNgCvK6q34Qy1moTcfohsixBGGMap5COYlLVY0CbcssO4IxqKr9tNnCJz/vFwOJQxldrKWnw6SNQeAxi4sMdjTHGBJU9SV0bqX1BS5zRTMYY08hYgqiF9w+mAPCHp/8fg2YuZ+GGrDBHZIwxwWMJooYWbsjiN//ex/fail4Ru7ylvy1JGGMaC0sQNeQp/f1VaWd6y07ASn8bYxoXSxA15Cnx/UVpNzpHfE8KP5RZbowxDZ0liBrylPheWdoHgAsiN5ZZbowxDZ0liBrylP7eqSnsKU3mgogvrfS3MaZRqatqro2Op8T3rCXbWHmkNz+P+pS/XtaV0Vbd1RjTSFiCqAVv6e9tBTDvA0a3/hboFO6wjDEmKKyJKRg6DYbIGNi+LNyRGGNM0FiCCIaY5nDaIEsQxphGxRJEsHS5GH7YBrnfhTsSY4wJCksQQfJBUW8A7p/1iJXdMMY0CpYggmDhhizuXHaY79zhrlZ2wxjTGFiCCAKn7EYpK0r7MCjiG+I5bmU3jDENniWIIPCU13i7ZBDxUsDlkf8ps9wYYxqiUM8olyQib4rIVhHZIiIDROQ1EfnS/ckQkS8r2DdDRDa5260NZZy15SmvsV67sK20PeMjPyyz3BhjGqJQ30E8Bryvqt2A3sAWVR2nqn1UtQ/wFmWnIy1vqLtteojjrBVP2Q0Q5pVcSJ+IXaRF77GyG8aYBi1kCUJEEoHBwHMAqlqoqrk+6wW4BpgXqhjqypi0VB4c25PUpDgWlpxHATH87fQvveU4jDGmIQplqY3OQA7wgoj0BtYBv1bVo+7684F9qrq9gv0VWCoiCjytqnP8bSQik4HJAB06dAhm/NXiLbsBMH85nbctgsKjzkN0xhjTAIWyiSkK6As8pappwFFgqs/68VR+9zBIVfsCo4DbRWSwv41UdY6qpqtqenJycpBCr52PEy+FgkNM+cMMeybCGNNghTJBZAKZqvqF+/5NnISBiEQBY4HXKtpZVbPd3/uBBUD/EMYaNAs3ZHHLR9FsL01lfORyeybCGNNghSxBqOr3wB4R8fTUDgM2u68vAraqaqa/fUWkuYgkeF4Dw4GvQxVrMHmeiZhbMoy+ETs4N2KzPRNhjGmQQj2K6U5groh8BfQB/tddfi3lmpdEJEVEFrtvTwU+FZGNwGpgkaq+H+JYg8Lz7MO8kgvJ0jb8LmouQqk9E2GMaXBCOh+Eqn4JnDREVVUn+lmWDVzivt6FMyy2wUlJiiMrN58CYphVNI5HY/7J6IjPWZt4cbhDM8bUd6WlUHgYjudBfq7z+3hu1a9j4mHyyqCHYxMGBdmUEV2ZNn8T+UUlvF06kF+VLua+6NdZe9FN4Q7NGFMXSop8LvC5J34fzy273F8SKDgEWlrJwQViWzo/cUkQmwSJ7SChXUg+iiWIIPOdijQ7N5/HIibwLH9g84KH+OsH45gyoqs9H2FMfaYKRccqv5B7XvtbV3S08uNHNnMv7i2dC3yLU6BtF+e173LfJOBZHpMAEXVXIckSRAh4nolYuCGLafMjWUZf/ifqbRbmDmLa/ELvNsaYECktcS/a5S/kuQFc+POgtKjy4zdLPHEhj0uC1p3Kvi/z2ndZS4huOCV4LEGEkDOiqYQ/y/W8F3M//4x5jHGFDzBryTZLEMZUpeh4AH/Bl2+Tdy/wBXmVH1siy128kyDptBOvyzfj+L5ulgiRTePS2TQ+ZZh4Ri59qz9hStEtzI55lPujXuYPuRPDG5gxdUEVCg5X0vZexV/zJQWVHz+6edkLeWJ7OKW7/4u697W7fUwLEKmb76EBswQRQp4RTQDvl/ZnTvGlTI5axIbSMxg0M876I0z9V77DtbIRNeWTwPG86ne4JvzE/wXeu8znIh8VU2dfQ1NlCSKEfEc0Afy1+Fp6R+zkwejnmJDX1vojTOiV6XCt4EJ+0gXf58IflA7XCtrm67jD1VSfJYgQ8h3RlJWbTwmR3FF4F6/G/In/F/MgtxfdxawlMZYgTOVKS5zhjwF3spZbF3CHq3vxbt25UXa4muoTVQ13DEGTnp6ua9fWz7mFOk1dhOebbsUhXoh5iB6SwX1Fk5lfOpiUJGtyatSKjgfQLONvXZ6THKjk/9OIqKo7Vq3D1VRARNZVNOeO/cuoI779ET+SyHWF9zM7+hH+FjObPsU7+GvutUybvwmwJqd6KeAO1wqacYqPV3786PiyF+8yHa5+/pr3vcjHNLcOVxMSdgdRR5xnIk70RwDEUMS9Ua9yY+T7fE8r7i/6FStK00i1u4nQKNPhmgfHf6z6oSffdYF2uPrrXC1zcU86eTildbiaMLE7iHqg/BPWChQSzZ+Lb2BRybnMjH6GF2Jmsar0Z8w5dCm/m3+8zH4G56/4wqOBP9Fafl21OlxbQvNkt8O1ovZ4n7/mmyVah6tpdOwOIkwGzVzubXICiKaYX0Yu5caof5MqB9hRmsL8kvPY2Pw8rh45jDF924cx2iDy/BVf5UW+gqaa0uLKj98s0c+494ou7uXeW4eraYIqu4OwBBEm/pqcAKIo5pKI1UyMep++ETsA2F16KqtKz+K72G70HzSMoeedD9GxdR+0Z8jk8UNOe3zBIffCfch97bmw+74u91PVX/GeDtcqL+4VdMpGRNbNd2FMI2EJop5auCHLOwTWn1M5yMWR6xgWsZ60iB0kyYmL6z5ak8WpZJW0pDCmFd3P6ES3ju2dv4Kj4yEy2ikn4LlgaqkzXLK0xHlCtbgASgqhKN/pQC06BoXHnCacoqNQcAQKjzjvPcmg4HAV7fA454xNPNHsUv6v9ZMu7L7bJTnxW4erMXXGEkQ9V9HdRFlKB9lPL9lFJ9lLB9lPh4j9JJNLKzlCKzlSuyCiYp3RMNHNndryMS2gWQv3dyI0S3DeN0t0EkCzRJ/x8z4JwUbUGNOghK2TWkSSgGeBHjgDuW8ERgA3AznuZr9T1cV+9h0JPAZEAs+q6sxQxhpO5R+o80/4Tk/lOz3V79pISoingFgKiJNCYigiklIiUASlhAj3J5JCjaKIKAqIJjY2niKJ5scjJbSMi0aOQ+6xIlKS4hjaLZkVW3PIzs131omzruzrPETyyD22q5JtQvs68FgtvoYea32PL5yxhuJZqpDeQYjIi8AnqvqsiMQA8cDdwBFVfbiS/SKB/wIXA5nAGmC8qm6uaB9ouHcQvgK7mzDGmJPFRUfy4Nie1UoSld1BhGxcnogkAoOB5wBUtVBVcwPcvT+wQ1V3qWoh8CpwRWgirV/GpKXy4NiepCY5I2qsscYYE6j8ohJmLdkWtOOFcuB2Z5xmpBdEZIOIPCsizd11d4jIVyLyvIi08rNvKrDH532mu+wkIjJZRNaKyNqcnBx/mzQ4Y9JS+WzqhWTMvJRHxvWxZGGMCVh2hc3U1RfKBBEF9AWeUtU04CgwFXgKOB3oA+wF/uZnX3/XQr9tYao6R1XTVTU9OTk5KIHXJ/6ShQBJcdG0io8GLHEYY05ISQre8zyh7KTOBDJV9Qv3/ZvAVFXd59lARJ4B3qtg35/6vG8PZIcq0IbCM5VpeZ7hsoF0hP14rAih0tJvxpgGKi46kikjugbteCFLEKr6vYjsEZGuqroNGAZsFpF2qrrX3exK4Gs/u68BuohIJyALuBb4RahibegqShwVCSSh2GiRphNfQ4q1vsfX2EYxhboW053AXHcE0y5gEvC4iPTB+SM2A7gFQERScIazXqKqxSJyB7AEZ5jr86r6TYhjbTKqm1CMMU2TPShnjDFNWFiGuRpjjGnYLEEYY4zxyxKEMcYYvyxBGGOM8atRdVKLSA7wbQ13bwv8EMRwwqmxfJbG8jnAPkt91Fg+B9Tus5ymqn6fMm5UCaI2RGRtRT35DU1j+SyN5XOAfZb6qLF8DgjdZ7EmJmOMMX5ZgjDGGOOXJYgT5oQ7gCBqLJ+lsXwOsM9SHzWWzwEh+izWB2GMMcYvu4MwxhjjlyUIY4wxfjX5BCEiI0Vkm4jsEJGp4Y6nptzZ+faLiL/y6Q2KiPxURFaIyBYR+UZEfh3umGpKRGJFZLWIbHQ/yx/CHVNtiEikO0Okv3lcGgwRyRCRTSLypYg06AqfIpIkIm+KyFb3/5kBQTt2U+6DEJFI4L/AxTiTFK0Bxqvq5rAGVgMiMhg4Arykqj3CHU9tiEg7oJ2qrheRBGAdMKaB/ncRoLmqHhGRaOBT4NequirModWIiPwWSAcSVfWycMdTUyKSAaSraoN/UE5EXgQ+UdVn3akV4lU1NxjHbup3EP2BHaq6S1ULgVeBK8IcU42o6sfAwXDHEQyquldV17uvDwNbqGBO8vpOHUfct9HuT4P8q0xE2gOXAs+GOxbjEJFEYDDwHICqFgYrOYAliFRgj8/7TBrohaixEpGOQBrwReVb1l9us8yXwH5gmc80vA3No8C9QGm4AwkCBZaKyDoRmRzuYGqhM5ADvOA2/T0rIs2DdfCmniDEz7IG+dddYyQiLYC3gLtV9VC446kpVS1R1T44c6v3F5EG1wQoIpcB+1V1XbhjCZJBqtoXGAXc7jbRNkRRQF/gKVVNA44CQetLbeoJIhP4qc/79kB2mGIxPtz2+reAuao6P9zxBIN7678SGBnmUGpiEDDabbt/FbhQRF4Ob0g1p6rZ7u/9wAKc5uaGKBPI9LkrfRMnYQRFU08Qa4AuItLJ7dy5FngnzDE1eW7H7nPAFlX9e7jjqQ0RSRaRJPd1HHARsDW8UVWfqk5T1faq2hHn/5Plqnp9mMOqERFp7g5+wG2OGQ40yNF/qvo9sEdEurqLhgFBG8wRFawDNUSqWiwidwBLgEjgeVX9Jsxh1YiIzAMuANqKSCYwXVWfC29UNTYIuAHY5LbdA/xOVReHMaaaage86I6YiwBeV9UGPUS0ETgVWOD8HUIU8Iqqvh/ekGrlTmCu+0fuLmBSsA7cpIe5GmOMqVhTb2IyxhhTAUsQxhhj/LIEYYwxxi9LEMYYY/yyBGGMMcYvSxDG+HArY/6Pz/sUEXkzROcaIyIPuK//JSJXBbBPjIh8LCJNeoi6qRuWIIwpKwnwJghVzVbVKi/cNXQv8M9ANxaRSLeo5IfAuBDFZIyXJQhjypoJnO7OEzBLRDp65tgQkYkislBE3hWR3SJyh4j81i2StkpEWrvbnS4i77uF4D4RkW7lTyIiZwIF5cpNDxaRz0Vkl+duQkQucOfGeAXY5G63ELgulF+CMdDEn6Q2xo+pQA+3uJ6nmqyvHjjVZWOBHcB9qpomIo8Av8SpeDoHuFVVt4vIOTh3CReWO84gYH25Ze2A84BuOCVfPE1b/d2Ydrvvvwb61eIzGhMQSxDGVM8Kd46KwyKSB7zrLt8E9HIr0A4E3nBLOQA083Ocdjhlmn0tVNVSYLOInOqzfLVPckBVS0SkUEQS3FiMCQlLEMZUT4HP61Kf96U4/z9FALmeO5BK5AMtKzm2byn6o372bwYcrzJaY2rB+iCMKeswkFDTnd15K3aLyNXgVKYVkd5+Nt0CnFGTc4hIGyBHVYtqGqcxgbAEYYwPVT0AfCYiX4vIrBoe5jrgVyKyEfgG/9PYfgykiU87VDUMBRpiZVvTwFg1V2PCREQeA95V1Q+qud98YJqqbgtNZMY47A7CmPD5XyC+Oju4Nf8XWnIwdcHuIIwxxvhldxDGGGP8sgRhjDHGL0sQxhhj/LIEYYwxxi9LEMYYY/z6/96xtez8YffyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(timestp,tnum3,'o',label=str()+' constant ambient temp')\n",
    "plt.plot(timestp,tp2,label='changing ambient temp')\n",
    "plt.title('Temperature vs time')\n",
    "plt.xlabel('time (hr)')\n",
    "plt.ylabel('temperature (F)')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1595,
   "metadata": {},
   "outputs": [],
   "source": [
    "def charttemp(r):\n",
    "    '''chart temp will print out a chart of the difference in time the ambient temperature and the final temperature\n",
    "       to see which ambient temperature will be the closest to 98.6'''\n",
    "    Amb_Tem= np.linspace(55,65,r)\n",
    "    difftime=np.linspace(-3,0,r)\n",
    "    NewTemp=np.empty(len(Amb_Tem))\n",
    "    print('Delta Time(hr)| Ambient Temperature (F) | New Temperature (F)')\n",
    "    print('--------------------------------------------------------------')\n",
    "    for i in range(0,len(Amb_Tem)):\n",
    "        NewTemp[i]=Amb_Tem[i]+(85-65)*np.exp(-kconst*difftime[i])\n",
    "    for j in range(0,len(Amb_Tem)):\n",
    "        print('{:10.6f} {:20.6f} {:30.6f}'.format(difftime[j],Amb_Tem[j],NewTemp[j]))\n",
    "\n",
    "    return \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1599,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Delta Time(hr)| Ambient Temperature (F) | New Temperature (F)\n",
      "--------------------------------------------------------------\n",
      " -3.000000            55.000000                     100.637615\n",
      " -2.842105            55.526316                      99.224704\n",
      " -2.684211            56.052632                      97.894194\n",
      " -2.526316            56.578947                      96.642584\n",
      " -2.368421            57.105263                      95.466522\n",
      " -2.210526            57.631579                      94.362797\n",
      " -2.052632            58.157895                      93.328335\n",
      " -1.894737            58.684211                      92.360193\n",
      " -1.736842            59.210526                      91.455554\n",
      " -1.578947            59.736842                      90.611719\n",
      " -1.421053            60.263158                      89.826104\n",
      " -1.263158            60.789474                      89.096235\n",
      " -1.105263            61.315789                      88.419744\n",
      " -0.947368            61.842105                      87.794362\n",
      " -0.789474            62.368421                      87.217919\n",
      " -0.631579            62.894737                      86.688333\n",
      " -0.473684            63.421053                      86.203615\n",
      " -0.315789            63.947368                      85.761857\n",
      " -0.157895            64.473684                      85.361235\n",
      "  0.000000            65.000000                      85.000000\n"
     ]
    }
   ],
   "source": [
    "charttemp(20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1519,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The time of death was at 810 AM\n"
     ]
    }
   ],
   "source": [
    "'''knowing the ambient temperature was around 55 F, we can now find the delta time and figure out the time of death'''\n",
    "tnew=98.6 # the new temperature\n",
    "kconst=.275\n",
    "T_am= 55\n",
    "dtnew= tnew-T_am\n",
    "hours= np.log(dt1/dtnew)/kconst # getting how many hours to add on to the existing time\n",
    "mins= 60\n",
    "hr=12\n",
    "minutes= hours*mins # using 24 hour clock 1100 would already have 12 hours mulitplied by 60 minutes\n",
    "time_mins=(hr*mins)+minutes \n",
    "act_time= time_mins+260\n",
    "\n",
    "print('The time of death was at {:.0f} AM'.format(act_time))\n",
    "\n"
   ]
  }
 ],
 "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
}