01_Getting-started-project.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computational Mechanics Project #01 - Heat Transfer in Forensic Science\n",
    "\n",
    "We can use our current skillset for a macabre application. We can predict the time of death based upon the current temperature and change in temperature of a corpse. \n",
    "\n",
    "Forensic scientists use Newton's law of cooling to determine the time elapsed since the loss of life, \n",
    "\n",
    "$\\frac{dT}{dt} = -K(T-T_a)$,\n",
    "\n",
    "where $T$ is the current temperature, $T_a$ is the ambient temperature, $t$ is the elapsed time in hours, and $K$ is an empirical constant. \n",
    "\n",
    "Suppose the temperature of the corpse is 85$^o$F at 11:00 am. Then, 2 hours later the temperature is 74$^{o}$F. \n",
    "\n",
    "Assume ambient temperature is a constant 65$^{o}$F.\n",
    "\n",
    "1. Use Python to calculate $K$ using a finite difference approximation, $\\frac{dT}{dt} \\approx \\frac{T(t+\\Delta t)-T(t)}{\\Delta t}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Answer 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.275\n"
     ]
    }
   ],
   "source": [
    "T = 85\n",
    "Ta = 65\n",
    "dT = (74-85)/2\n",
    "K = dT/(Ta-T)\n",
    "print(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": "markdown",
   "metadata": {},
   "source": [
    "# Answer 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "metadata": {},
   "outputs": [],
   "source": [
    "def measure_K(Temp_t1,Temp_t2,Temp_ambient,delta_t):\n",
    "    ''' Determine the value of K based upon temperature of corpse \n",
    "    when discovered, Temp_t1\n",
    "    after time, delta_t, Temp_t2\n",
    "    with ambient temperature, Temp_ambient\n",
    "    Arguments\n",
    "    ---------\n",
    "    Temp_t1 - initial temperature when corpse is discovered\n",
    "    Temp_t2 - final temperature\n",
    "    Temp_ambient - the ambient temperature surrounding the corpse\n",
    "    delta_t - the time elapsed between Temp_t1 and Temp_t2\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    K - constant coefficient of cooling\n",
    "    \n",
    "    '''\n",
    "    dT = (Temp_t2 - Temp_t1)/delta_t\n",
    "    K = dT/(Ta-T)\n",
    "    return K\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.275"
      ]
     },
     "execution_count": 241,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "measure_K(85,74,65,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. A first-order thermal system has the following analytical solution, \n",
    "\n",
    "    $T(t) =T_a+(T(0)-T_a)e^{-Kt}$\n",
    "\n",
    "    where $T(0)$ is the temperature of the corpse at t=0 hours i.e. at the time of discovery and $T_a$ is a constant ambient temperature. \n",
    "\n",
    "    a. Show that an Euler integration converges to the analytical solution as the time step is decreased. Use the constant $K$ derived above and the initial temperature, T(0) = 85$^o$F. \n",
    "\n",
    "    b. What is the final temperature as t$\\rightarrow\\infty$?\n",
    "    \n",
    "    c. At what time was the corpse 98.6$^{o}$F? i.e. what was the time of death?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {},
   "outputs": [],
   "source": [
    "def T_analytical(Temp_ambient,Temp_t1,K,t):\n",
    "    '''This function returns the analytical solution for a first order\n",
    "    Thermal system.\n",
    "        Arguments\n",
    "    ---------\n",
    "    Temp_ambient - the ambient temperature surrounding the corpse\n",
    "    Temp_t1 - initial temperature when corpse is discovered\n",
    "    K - constant coefficient of cooling\n",
    "    t - the time elapsed between Temp_t1 and Temp_t2\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    T_analytical - analytical solution for a first order thermal system\n",
    "    '''\n",
    "    T_analytical = Temp_ambient + (Temp_t1 - Temp_ambient)*np.exp(-K*t)\n",
    "    return T_analytical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[85.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]\n",
      "[85.         82.25       79.878125   77.83238281 76.06793018 74.54608978\n",
      " 73.23350243 72.10139585 71.12495392 70.28277275 69.5563915  68.92988767\n",
      " 68.38952812 67.923468   67.52149115 67.17478612 66.87575303 66.61783698\n",
      " 66.3953844  66.20351904 66.03803518 65.89530534 65.77220085 65.66602324\n",
      " 65.57444504 65.49545885 65.42733326 65.36857493 65.31789588 65.2741852\n",
      " 65.23648473 65.20396808 65.17592247 65.15173313 65.13086983 65.11287522]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "time_step = 0.5\n",
    "t = np.arange(0,18,time_step)\n",
    "Temp_ambient = 65\n",
    "Temp_t1 = 85\n",
    "K = 0.275\n",
    "T_numerical = np.zeros(len(t))\n",
    "T_numerical[0] = 85\n",
    "print(T_numerical)\n",
    "for i in range (1,len(t)):\n",
    "    T_numerical[i] = T_numerical[i-1] + (-K*(T_numerical[i-1] - Temp_ambient))*time_step\n",
    "print(T_numerical)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Answer 3a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Temperature (F)')"
      ]
     },
     "execution_count": 253,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3zU9f3A8df7sgkQVhghQYagAkkYARyAKIpABRWFgrhwoK2Iq1V+1SJFbW1FrRQXipa2iiKKE5UhKKAgAZkCsiHMMAMSQsb798ddMITvJZeQu8t4Px+Pe9zdd9z3nS9H3vlsUVWMMcaYwlzBDsAYY0z5ZAnCGGOMI0sQxhhjHFmCMMYY48gShDHGGEehwQ6gLNWrV0+bNm0a7DCMMabCWLp06X5VjXXaV6kSRNOmTUlNTQ12GMYYU2GIyDZv+6yKyRhjjCNLEMYYYxxZgjDGGOOoUrVBGGMqvuzsbNLS0jhx4kSwQ6lUIiMjiY+PJywszOdzLEEYY8qVtLQ0atSoQdOmTRGRYIdTKagqBw4cIC0tjWbNmvl8nl+rmETkQRFZIyKrRWSKiESKyBgR2Skiyz2Pvl7O7S0i60Vko4iM8luQK6fCC21hTC3388qpfruUMaZ4J06coG7dupYcypCIULdu3RKXyvxWghCRxsBIoLWqZorIVGCwZ/cLqjquiHNDgJeAK4E0YImIfKKqP5VpkCunwqcjITvT/f7IDvd7gKRBZXopY4zvLDmUvdLcU383UocCUSISClQDdvl4Xmdgo6puVtWTwLvANWUe3ZyxvyaHfNmZ7u3GGFPF+S1BqOpOYBywHdgNHFHVmZ7dI0RkpYi8KSK1HU5vDOwo8D7Ns+0MIjJcRFJFJDU9Pb1kQR5JK9l2Y0yVMX36dESEdevWlfozbrvtNqZNm1bkMX/9619Pe3/xxReX6lpjxoxh3DivFTOl4rcE4fnFfw3QDIgDokXkJuAVoAXQDnfieM7pdIdtjisbqepEVU1R1ZTYWMfR4t7FxJdsuzGmypgyZQpdu3bl3Xff9et1CieI7777zq/XKwl/VjFdAWxR1XRVzQY+BC5W1b2qmquqecDruKuTCksDEgq8j8f36inf9RwNYVGnbwuLcm83xlRZx44dY+HChUyaNOlUgpg3bx49evTghhtu4Pzzz2fo0KHkr8g5duxYOnXqRNu2bRk+fDiFV+qcM2cO11133an3s2bNYsCAAYwaNYrMzEzatWvH0KFDAahevfqp4/7xj3+QmJhIcnIyo0a5++q8/vrrdOrUieTkZK6//nqOHz/ut/vgz26u24ELRaQakAn0BFJFpJGq7vYccx2w2uHcJUBLEWkG7MTduH1jmUfoaYjW2WMgYydZEkVkv/HWQG1MOfGXT9fw066MMv3M1nE1eaJfmyKP+eijj+jduzetWrWiTp06LFu2DIAff/yRNWvWEBcXxyWXXMLChQvp2rUrI0aMYPRo9x+WN998M5999hn9+vU79XmXX3459957L+np6cTGxvLWW28xbNgw+vXrx4QJE1i+fPkZMXzxxRd89NFHLF68mGrVqnHw4EEABgwYwF133QXA448/zqRJk7jvvvvK5N4U5s82iMXANGAZsMpzrYnAP0RklYisBC4DHgQQkTgRmeE5NwcYAXwFrAWmquoavwSaNAh56CdW1+zGkbwosi4Y4JfLGGMqjilTpjB4sLvT5eDBg5kyZQoAnTt3Jj4+HpfLRbt27di6dSsAc+fOpUuXLiQmJvL111+zZs3pv65EhJtvvpn//e9/HD58mO+//54+ffoUGcPs2bMZNmwY1apVA6BOnToArF69mm7dupGYmMjbb799xrXKkl8HyqnqE8AThTbf7OXYXUDfAu9nADP8F93pQi+4mgaL55O65BtSLu4ZqMsaY4pQ3F/6/nDgwAG+/vprVq9ejYiQm5uLiNC3b18iIiJOHRcSEkJOTg4nTpzg97//PampqSQkJDBmzBjH8Qb5JYbIyEgGDhxIaGjRv35V1bFr6m233cZHH31EcnIy//73v5k3b95Z/8ze2FxMHs0vGUCuCoeXTQ92KMaYIJo2bRq33HIL27ZtY+vWrezYsYNmzZqxYMECx+Pzk0G9evU4duyY115LcXFxxMXF8dRTT3Hbbbed2h4WFkZ2dvYZx/fq1Ys333zzVBtDfhXT0aNHadSoEdnZ2bz99ttn86MWyxKER0TN+myJTuac/fPIzs0LdjjGmCCZMmXKaQ3KANdffz3vvPOO4/G1atXirrvuIjExkWuvvZZOnTp5/eyhQ4eSkJBA69atT20bPnw4SUlJpxqp8/Xu3Zv+/fuTkpJCu3btTnVhffLJJ+nSpQtXXnkl559/fml/TJ9I4db2iiwlJUXPZsGgddOf4fwVf2NJ/7l06tChDCMzxvhq7dq1XHDBBcEOwy9GjBhB+/btueOOO4Jyfad7KyJLVTXF6XgrQRTQ7JKBAOxL/TDIkRhjKpuOHTuycuVKbrrppmCH4jObzbWAiPot2BnejIa75pCb9yQhLpsPxhhTNpYuXRrsEErMShCFHG/Wi3a6luXrNgY7FGOMCSpLEIXEX3gDIaJsX/xRsEMxxpigsgRRSNQ5KRwKqUftHbPIy6s8DfjGGFNSliAKc7k4lNCTzrnLWbF1T7CjMcaYoLEE4aBBp+upJlls+P6zYIdijAkCEeHhhx8+9X7cuHGMGTMmoDGkpqYycuTIUp3bo0cPzqbLfz5LEA6iz+tBplSj2uYvz5iV0RhTzvhh2eCIiAg+/PBD9u/fXwYBllxOTg4pKSmMHz8+KNfPZwnCSWgE6Q270SXnB9akHQp2NMYYb/KXDT6yA9Bflw0+yyQRGhrK8OHDeeGFF87YV3gRoPzpuefNm8ell17KoEGDaNWqFaNGjeLtt9+mc+fOJCYmsmnTJgDS09O5/vrr6dSpE506dWLhwoWAe8Gf4cOH06tXL2655RbmzZvH1VdfDbinHx82bBiJiYkkJSXxwQcfAPC73/2OlJQU2rRpwxNPFJ727uzZOAgv6nQcQPXdXzFz8WzaJtj038YExRejYM8q7/vTlkBu1unbsjPh4xGwdLLzOQ0Toc8zxV763nvvJSkpiUceecTncFesWMHatWupU6cOzZs358477+SHH37gxRdf5F//+hf//Oc/uf/++3nwwQfp2rUr27dv56qrrmLt2rWAe6zEggULiIqKOm0SvieffJKYmBhWrXLfi0OH3H+4Pv3009SpU4fc3Fx69uzJypUrSUpK8jne4liC8KJ6m97kfBaCrJ+B6kBbRN2Y8qhwcihuewnUrFmTW265hfHjxxMVFVX8CUCnTp1o1KgRAC1atKBXr14AJCYmMnfuXMA9jfdPP/106pyMjAyOHj0KQP/+/R2vNXv27NNWtqtd271S89SpU5k4cSI5OTns3r2bn376yRJEQETV4kC9TnTZt4jlOw7TvonT0tnGGL8q7i/9F9p6qpcKiUmAYZ+f9eUfeOABOnTowLBhw05tCw0NJS/PPaGnqnLy5MlT+wpOB+5yuU69d7lc5OTkAJCXl8f333/vmAiio6Md43Ca+nvLli2MGzeOJUuWULt2bW677TbHacbPhrVBFCGm/bW0cO1mwfflZ41YY0wBfl42uE6dOgwaNIhJkyad2ta0adNT02Z8/PHHjlN1F6VXr15MmDDh1Hun1eSKO+fQoUNkZGQQHR1NTEwMe/fu5YsvvihRHL7wa4IQkQdFZI2IrBaRKSISKSLPisg6EVkpItNFpJaXc7d6Vp5bLiJn31+rFCLbuBuI8tZ9blOAG1MeJQ2CfuPdJQbE/VzGywY//PDDp/Vmuuuuu/jmm2/o3Lkzixcv9vpXvzfjx48nNTWVpKQkWrduzauvvlrsOY8//jiHDh2ibdu2JCcnM3fuXJKTk2nfvj1t2rTh9ttv55JLLinxz1Ycv033LSKNgQVAa1XNFJGpuFeI2wV8rao5IvJ3AFV91OH8rUCKqvrcz+xsp/t2kvHPi9hwMJvDQz6n5wUNyvSzjTFnqszTfQdbeZvuOxSIEpFQoBqwS1VnetacBlgExPs5hrMSnXwN7V0bmb2kiJ4UxhhTCfktQajqTmAcsB3YDRxR1ZmFDrsd8FZxpsBMEVkqIsO9XUdEhotIqoikpqenl0XopwkJDcOF8tdN15P3fJsyGYRjjDEVgd8ShIjUBq4BmgFxQLSI3FRg/2NADuBtUdVLVLUD0Ae4V0S6Ox2kqhNVNUVVU2JjY8v0Z2DlVPj2WU+84MpIK5NBOMaYotkMBmWvNPfUn1VMVwBbVDVdVbOBD4GLAUTkVuBqYKh6iVpVd3me9wHTgc5+jNXZnLHuQTcFZWe6txtj/CIyMpIDBw5YkihDqsqBAweIjIws0Xn+HAexHbhQRKoBmUBPIFVEegOPApeq6nGnE0UkGnCp6lHP615A4H8rH0kr2XZjzFmLj48nLS0Nf1QZV2WRkZHEx5esyddvCUJVF4vINGAZ7qqkH4GJwBogApjlGfixSFXvEZE44A1V7Qs0AKZ79ocC76jql/6K1auYeC+DcMp1u7oxFVpYWBjNmjULdhgGP4+kVtUngMIzSJ3r5dhdQF/P681Asj9j80nP0e42hwLVTFmEE1FGg3CMMaY8s5HURSk8CAeYkZPCuvq9gxuXMcYEgCWI4iQNggdXw5jD5DRsz3munUz/cWewozLGGL+zBFECoe0G09q1jVVLv7f1qo0xlZ4liJJoO4A8CaHribks2nIg2NEYY4xfWYIoier10eaXcV3IQj5e5tC7yRhjKhFLECUUkjyYRnKAPavncfxkTvEnGGNMBWUJoqTO70tuaDWuyv2Gz1fuDnY0xhjjN5YgSio8Glfrq+kf+gPTFm8MdjTGGOM3liBKQZJ+S3V+odbOeazbkxHscIwxxi8sQZRGsx7kVYvl+tCFvPuDNVYbYyonSxClERKKK2kgl4csZ9aydWSezA12RMYYU+YsQZRW0iBCNZtu2d8xY5U1VhtjKh9LEKXVqB1arxVDIr5nyg/bgx2NMcaUOUsQpSWCJA4iOW8Nu7f9zM97jwY7ImOMKVOWIM5G4g0AXBf2nTVWG2MqHUsQZ6NOM0i4kKFRi/hg6Q5OZFtjtTGm8vBrghCRB0VkjYisFpEpIhIpInVEZJaIbPA81/Zybm8RWS8iG0VklD/jPCtJg2h0chuNszby5eo9wY7GGGPKjN8ShIg0BkYCKaraFggBBgOjgDmq2hKY43lf+NwQ4CWgD9AaGCIirf0V61lpcx2Ki/cjnqT/x23ghbawcmqwozLGmLPm7yqmUCBKREKBasAu4Bpgsmf/ZOBah/M6AxtVdbOqngTe9ZxX/mycjYgQTSYu1L2G9acjLUkYYyo8vyUIVd0JjAO2A7uBI6o6E2igqrs9x+wG6juc3hgo2Oqb5tl2BhEZLiKpIpKanp5elj+Cb+aMBS3U9pCd6d5ujDEVmD+rmGrj/qu/GRAHRIvITb6e7rDNcQk3VZ2oqimqmhIbG1u6YM/GkbSSbTfGmArCn1VMVwBbVDVdVbOBD4GLgb0i0gjA87zP4dw0IKHA+3jc1VPlT0x8ybYbY0wF4c8EsR24UESqiYgAPYG1wCfArZ5jbgU+djh3CdBSRJqJSDjuxu1P/Bhr6fUcDWFRp206QQTac3SQAjLGmLLhzzaIxcA0YBmwynOticAzwJUisgG40vMeEYkTkRmec3OAEcBXuJPKVFVd469Yz0rSIOg3HmLcBZ48XIw6OYzvoy8PcmDGGHN2RNWxar9CSklJ0dTU1OAFsG4GvDuEP7j+wOFz+vDGrSnBi8UYY3wgIktV1fGXlY2kLkutroKYBO6r8Q1z1u1l6/5fgh2RMcaUmiWIsuQKgZRhnHNkCa1cu3hr4ZZgR2SMMaVmCaKstb8FQsJ5rMH3vL80jSOZ2cGOyBhjSsUSRFmrHgutr+WSozPh5C+8t8TWijDGVEyWIPyh052EZB/lgQbLmfzdNnJy84IdkTHGlFiRCUJEwkXkWhF5zjMb65si8pCInB+oACukhM7QMJEhfMXOw8f5co3N8mqMqXi8JggReRxYDFwGrMA9sd4nuCfge0FEvhSRtgGJsqIRgU53UuPIevrGbOfNBdZYbYypeEKL2LdKVZ/ysu8fnmkyErzsN4kDYeZoHq79LT23nsOP2w/Rvonj0hfGGFMuFVXFVOTUFqq6W1V/KON4Ko/waGh3I833zaZpxDEmWSnCGFPBFJUglua/EJF/BiCWyqfTHUheNk/EL+WL1XvYeTgz2BEZY4zPikoQBafc7u7vQCqlei2heQ+6HfmUEHJ5Y/7mYEdkjDE+KypBVJ5JmoKp052EHtvFoy22MeWH7ew/lhXsiIwxxideJ+sTkePAOtwlifM8r/G8V1XtEJAISyDok/U5yc2BZ1uQl50JOVkcjWhIzNVPumeBNcaYICtqsr6iejEl+imeqmXNh3DyF1x52SAQc3IP+slId/2dJQljTDnmNUGo6qZABlJpzRkLeafPxyQ5njWrLUEYY8qxogbKzRWR34lIXKHtoSLSXUQmicgw/4dYwXlZm1ptzWpjTDlXVBXTb4A7geki0hg4CER6HnOAl1TVa4W/iJwHvFdgU3NgNHAR7jYNgFrAYVVt53D+VuAokAvkeKsjK/di4uHIjjM2H41oQM0ghGOMMb4qqorpODAeGC8iEUB9IFNV9/vywaq6HmgHICIhwE5guqqeGlMhIs8BR4r4mMt8vV651XM0fDoSsn8dA5ElEfwtaxCPZ+UQHVFUjjbGmODxaTZXVc1S1R1n8cu6J7BJVbflbxARAQYBU0r5mRXDaWtWu4eWHGvWhyknLuR/i7YVfa4xxgRRoKb7HsyZiaAbsFdVN3g5R4GZIrJURIZ7+2ARGS4iqSKSmp6eXkbhlrGkQfDgahhzGJp2o276Ynq0qMnr8zeTeTI32NEZY4wjvycIEQkH+gPvF9o1hKJLD5d4xlr0Ae4VEcfR3Ko6UVVTVDUlNja2TGL2q+5/gKO7eSJhOfuPneRdW1DIGFNO+ZQgRCReRC7zvI4QkegSXKMPsExV9xb4vFBgAKc3Yp9GVXd5nvcB04HOJbhm+dXsUmicQrN1E7moaU1e+2YzWTlWijDGlD/FJggRuR33zK5veDadA3xcgms4lRSuANapqmNfTxGJFpEa+a+BXsDqElyz/BKB7n+Ew9v5S7O17Mk4wfup1uXVGFP++FKCGAlcCGQAqOrPuHs0FUtEqgFXAh8W2nVGm4SIxInIDM/bBsACEVkB/AB8rqpf+nLNCqHVVdAgkZY/T6RTk5r86+sN1hZhjCl3fEkQJ1T1ZP4bT5dVKeL4U1T1uKrWVdUjhbbfpqqvFtq2S1X7el5vVtVkz6ONqj7ty/UqDBHo/jByYAN/vWArezOyeOs7Wy/CGFO++JIgForII0Ckpx3iPeAz/4ZVBVzQH+q2pOW617ji/FhembeJQ7+cLP48Y4wJEF8SxCO4RzSvA+7HPYr6MX8GVSW4QqDbQ7B3FWMu2MkvWTm8PG9jsKMyxphTikwQnuqkN1X1FVW9TlWv9bzOC1B8lVviQKjVhPiVE7i+fWMmf7eNtEPHgx2VMcYAxSQIVc0FGolIWIDiqVpCwuCSB2BnKo+evw8EXpjlbdygMcYEli9VTJuB+SLyfyIyMv/h78CqjHZDoUYj6i0bz7CLm/Lhj2ms25MR7KiMMcanBJEOzAKqAbEFHqYshEXCxffB1vk8urofmyKGUu/1jrByarAjM8ZUccVOJaqqfw5EIFVahHvib9dx91xS9XL3kfvxSELAFhUyxgRNsQlCRGbhnjjvNKrayy8RVUXf/P2MTSG5meicsYglCGNMkPiyGMHjBV5HAtcDWf4Jp4rytrqcrTpnjAkiX6qYFhfa9I2IfOOneKomL6vO7ZV61M7JJSI0JAhBGWOqOl8m66tZ4FFLRHoCjQIQW9XRczSERZ22KTckkr9mDWTSApuCwxgTHL5UMa3B3QYhQA6wBbjLn0FVOfntDHPGeqqVlJDEAZzIuJ5/zdnIte0aE1crqsiPMMaYsuZLN9fmqtpEVRNUtZmqXg4s9HdgVU7+qnNPHIJzusLPXzH6ynjyVHl6xtpgR2eMqYJ8SRCF2yDAPQW38QcRuOopOH6A+NWv8Pse5/L5yt0s3Fja5cCNMaZ0vCYIEakvIslAlIgkikiS59EV96A54y9x7SFpMCx6hXuSQ2lSpxpPfLKGkzk2BZYxJnCKKkH8BpgAxAMvAy95Hn8Cih08JyLnicjyAo8MEXlARMaIyM4C2/t6Ob+3iKwXkY0iMqrkP1oF1/PPIELEN08x+urWbNx3jMnfbQ12VMaYKkRUzxgDd/oBIoNU9azmffDMCrsT6AIMA46p6rhijv8Z92p0acASYIiq/lTUdVJSUjQ1NfVsQi1f5jwJ88fBnV9z++w8Fm8+wNd/6EGDmpHBjswYU0mIyFJVTXHaV2wbhKpOFZGrROQhEflT/qOEMfQENqnqNh+P7wxs9KwsdxJ4F7imhNes+Lo+ANH14as/8cTVF5Cdq/zNGqyNMQHiyziIl4FbgYeAKOAm4NwSXqfwGtQjRGSliLwpIrUdjm8MFBw5lubZVrVE1IDL/gQ7FnHOvjncfWlzPlq+i8WbDwQ7MmNMFeBLL6auqnojcMAzcV8X3O0SPhGRcKA/8L5n0ytAC6AdsBt4zuk0h22OdWEiMlxEUkUkNT093dewKo72N0PsBTDrCX7frQmNa0Ux+uM1ZOdag7Uxxr98SRAn8p9FpKHnfdMSXKMPsExV9wKo6l5VzfWsSvc67uqkwtKAhALv44FdTh+uqhNVNUVVU2JjK+Es5CGh0OspOLSFqBfPZ8GJAUw6NIw5UycEOzJjTCXnS4KYISK1gHHAcmArMK0E1xhCgeolESk4Tcd1wGqHc5YALUWkmacEMhj4pATXrFwyD4K44MRhBCXetZ/u655k5/zJwY7MGFOJFbcmtQv4QlUPq+r7QDMgUVV9aqQWkWq4eyJ9WGDzP0RklYisBC4DHvQcGyciMwBUNQcYAXwFrAWmquqakv1olcicsVBoGfBqcpKwuU+SY1VNxhg/KXIuJlXNE5EXgQs97zOBTF8/XFWPA3ULbbvZy7G7gL4F3s8AZvh6rUrNy7Tf9XL389r8LfyuR4sAB2SMqQp8qWKaJSJVr4tpeRLj3CfgUFgsL8z+mY37jgU4IGNMVeBLghgBTBeRTBE5KCKHROSgvwMzBThMB44rjNBeY6gWHsIfp60gN6/oAY/GGFNSviSIekAYUB2I9byvhN2FyrGkQdBvPMQkAOJOFppHTLNOPNGvNT9uP8xbC23dCGNM2fJlJHUuMBB41PO6Ee4xDCaQ8qcDH3MY7l8JkTXhkxFcm9yInufX59mv1rNl/y/BjtIYU4n4MpJ6Au7eRvmNy8eBV/0ZlClG9frQ+xnYsRhZMomnr0skPNTFo9NWkmdVTcaYMuJLFdPFqno3ngFzqnoQCPdrVKZ4Sb+Fc6+A2WNoqPv489Wt+WHrQV6fvznYkRljKglfEkS2ZzyEAohIXcA63webCFz9gvv1pw8wsENj+rRtyLNfrWf5jsPBjc0YUyn4kiBeAj4AYkXkL8AC4O9+jcr4plYTuGIMbJqDrHyPZwYk0aBmJCOn/MjRE9nBjs4YU8H50kj9H+Bx3FNtHAQGquq7/g7M+KjTnZBwIXw5ipi8Q7w4uB07D2fy2PTVFLfWhzHGFMWXEgRACJANnCzBOSYQXC7o/y/IPg7vDCZlenc2ht/II2tvYPHH1pfAGFN6vvRiegz3ZHtxuGdVfUdE/s/fgZkSiG0F5/WFXUvhyI5TE/ol/ziavQv/E+zojDEVlC+lgZuATqr6uKo+hnt67lv8G5YpsZ1nLrUaJSdhzliycnKDEJAxpqLzJUFs4/RJ/UIB60tZ3hzZ6bg5Nnc/f5uxLsDBGGMqA18SxHFgjYi8ISKvA6uAwyLyvIg879/wjM+8TOiXEV6ff3+3la/W7AlwQMaYiq7I6b49Pvc88i3yUyzmbPQcDZ+OhOwCs7GHRBD9m7Ekzo/h4akraHFvNOfWrxG8GI0xFUqxCUJVJwUiEHOWkga5n+eMda8fIS4IjybsvKt4rWkk/Scs4K7/LOWjey8hJiosuLEaYyoEX3ox9RaRJSKyryTTfYvIeSKyvMAjQ0QeEJFnRWSdiKwUkeme5Uydzt/qWXluuYic2QJrzlRwQr87ZkHWUfjoXuJiInl5aEd2HDzO/e/+aFODG2N84ksbxATgbqAxJZjuW1XXq2o7VW0HdMTdljEdmAW0VdUk4GegqC6zl3k+I8WHOE1B8R3hyrGw/nNY9DKdm9VhTP82zFufznMz1wc7OmNMBeBLG0QasFxVz2b+pZ7AJlXdhrtXVL5FwA1n8bmmKBf+DrYugFmjIaELQ7t0ZM2uI7w8bxOt42pydVJcsCM0xpRjvpQgHgE+FZE/isjI/EcJrzMY92C7wm4HvvByjgIzRWSpiAz39sEiMlxEUkUkNT09vYRhVXIicO1LUDMO3h+GZB7iL/3b0vGc2vzx/ZX8tCsj2BEaY8oxKW6+HhH5Avc0G6soMIurqv7ZpwuIhAO7gDaqurfA9seAFGCAOgQhInGquktE6uOulrpPVb8t6lopKSmammrNFWdIWwpvXgX1W0PmQfRIGnuoy2thNzHygceoE22ztxtTVYnIUm/V+L5UMdVX1Y5ncf0+wLJCyeFW4Gqgp1NyAFDVXZ7nfSIyHfcI7iIThPEiviO0HQAr3wNAgEbs55GTL/P6xHDuHjGKyLCQ4MZojCl3fKlimiMil5/FNYZQoHpJRHoDjwL9VfW40wkiEi0iNfJfA72A1WcRg9n23RmbqslJrj/0Jg9NXW4r0RljzuBLgrgLmC0ix0rSzRVARKoBVwIfFtg8AagBzPJ0YX3Vc2yciMzwHNMAWCAiK4AfgM9V9Usffybj5Eia4+bGrgPMWLWHsZ/9ZNODG2NO40sVU73SfrinhFC30LZzvRy7C+jreb0ZSC7tdY2DmHg4suOMzRITzx0pzZi0YAsNYyK559IWQQjOGFMe+bJgUC4wEHjU87oR0M7fgY+QGfUAABr1SURBVJky1nM0hEWdvk1ccNmfeKzvBfRLjuOZL9bx4TLnkoYxpurxZST1BOAy4GbPpuOArURT0SQNgn7jISYBEIiqA5oHW+bjEhg3MImLmtflkWkrmb/BugsbY3xrg7hYVe8GTgCo6kHA+kVWRAWn4nh0C1w6Cla8A9/8nYjQEF67pSPn1q/OPf9dyuqdR4IdrTEmyHxJENki4sI9cA0RqUuB8RCmAusxCpJvhHl/g+XvUDMyjMm3d6ZWtXBuffMHNuw9GuwIjTFB5DVBiEh+A/ZLwAdArIj8BVgA/D0AsRl/E4F+L0Kz7vDJfbB5Hg1qRvLfOzrjcglDXl/Mxn3Hgh2lMSZIvI6kFpFlqtrB87oNcAXuMVazVbVcjkmwkdSldOIIvNkbDmyGqFpwbC/Z1eMY88sNzArtzrvDL6R5bPVgR2mM8YOiRlIXVcUk+S9UdY2qvqiq/yyvycGchcgY6HAr5J6AY3sAJezYTp4MmciVOd8w5PVFbNn/S7CjNMYEWFHjIGJF5CFvO1XVlhutTL6fcMYmV04mY6p/wBeZ3RkycRHv3X0h59SNDkJwxphgKKoEEQJUxz3q2elhKhMvI63Dju3i7Tu7kJWTy5CJi9h+wHF2FGNMJVRUCWK3qo4NWCQmuLyMtKZmHBc0qsn/7uzC0DcWM+T1Rfzvzi40q2clCWMqO5/aIEwV4DTSGsAVCpmHaBMXw//u6EJmdi43vPKdjZMwpgooKkH0DFgUJvgKj7SOSYCL7oOju+E/18Dxg7RtHMO0ey4iMiyEwRMX8d3G/cGO2hjjR8UuGFSRWDdXP9gwC94dCrGt4JZPoFod9hw5wS1vLmbr/uO8OLgdfRIbBTtKY0wplbabqzHQ8koY8g6k/wyT+8MvB2gYE8nUuy8iMT6G37+zjHcWbw92lMYYP/Blum9T1Z17BQyZAu/eCK92A1FqZezm/ZqNeT3uJv40HQ7+ksW9l52LiDVdGVNZWAnC+ObcntDld3B0J2TsAhRXRhrDM17kyeY/MW7mz4z6YBUnc2yaLmMqC78lCBE5z7NiXP4jQ0QeEJE6IjJLRDZ4nmt7Ob+3iKwXkY0iMspfcZoSWD3tjE2SnclNv0xmxGXn8l7qDm56YzEHjmUFIThjTFnzW4JQ1fWq2k5V2wEdca8jMR0YBcxR1ZbAHM/704hICO5JAvsArYEhItLaX7EaH3kZTCdH0vjDVefx4uB2rEg7TP8JC1m7OyPAwRljylqgqph6AptUdRtwDTDZs30ycK3D8Z2Bjaq6WVVPAu96zjPBFBPvvD3KXQi8pl1jpt59Edm5eVz/ynfMXLMngMEZY8paoBLEYGCK53UDVd0N4Hmu73B8Y6DgsN40z7YziMhwEUkVkdT0dFsJza+8LVuaeRDmPAl5eSQn1OLT+7rSsn517v7fUl6au5HK1JXamKrE7wlCRMKB/sD7JTnNYZvjbxlVnaiqKaqaEhsbW5oQja+cBtNd8zK0vxnmj4MP7oDsTBrUjOS9uy+iX1Icz361nt+/vYyME9nBjt4YU0KB6ObaB1imqns97/eKSCNV3S0ijYB9DuekAQkF3scDu/wcp/FF0iD3o6DkwVCvJcwaDbt+hNwsIjN282JMPP3b383dK4Q14xcw4cb2JMXXCk7cxpgSC0QV0xB+rV4C+AS41fP6VuBjh3OWAC1FpJmnBDLYc54pj0Tgkvuhyz1waMupbrByZAdXbHyaWVfsIcfTLvHWwi1W5WRMBeHXBCEi1YArgQ8LbH4GuFJENnj2PeM5Nk5EZgCoag4wAvgKWAtMVdU1/ozVlIF1n5+5LTuT5iue5/OR3ejeMpa/fPoT9/xvKUcyrcrJmPLO5mIyZWdMLZybigTGHEZVeWP+Fv7+5ToaxkTyryHtad/EcRiMMSZAbC4mExjeusG6QmDfOkSEu7o3Z+o9F6EKN7z6Pc9+tY6snNzAxmmM8YklCFN2nLrBhkRASCRM7AHL/guqdGhSmxn3d+O69o15ae4mrpmwkDW7bH0JY8obSxCm7Dh2g50AI5dCQif4ZAR8cCcsnUzMq+0Zt+ZS1tb5AykZs7lmwkLGz9lAdq7N5WRMeWFtECYw8nJhwfPw9VO4h7n8+r3T0Cj+U+8hntjahsTGMTw3KJlWDWzZc2MCwdogTPC5QqD7HyG6PoUbsiUnk1sz/8MrQzuw83Amvxk/n2e/WkfmSWubMCaYLEGYwPrFy3QoR9Lok9iImQ92p19yHC/N3cSVL3zDnLV7nY83xvidJQgTWN56OoVHQ9Yx6lWP4PlB7Xh3+IVEhYVwx+RUhv8nlZ2HMwMbpzHGEoQJMKeeTq4QOHkMXr4Qfp4JK6dy4ceXMjPjGlbFPEzNjdO54rlveHneRusSa0wA2ZKjJrDy53GaM9a9vkRMvDtp1DoHPh0J7wwECQHNRYAaWbv5R/gkzq9dg6e+zOWdxdv541Xn0S8pDpfLljc1xp+sF5MpP3JOwrMtIMthsaGYBOZfPZe/zljH2t0ZJMfH8Ke+F9Cled3Ax2lMJWK9mEzFEBoOWUed9x1Jo1vLWD67ryvjBiazNyOL305cxF3/SWVT+rHAxmlMFWEJwpQv3hqxQ8Nh90pCXMINHeOZ+4ce/PGq8/h+0wF6vfAtj0xbwfYDxwMbqzGVnCUIU744NmKHudslXusO0++BwzuIWvcB9y6/llXyW5ZGP0DOiqlc9tw8SxTGlCFrpDbli7dG7Ja9YMELsOgVWDnVvQZFXg4C1Mrey7iISVwSX5f/Wy58sGwnA9o3ZsTl53JO3eig/jjGVGTWSG0qlsM74KXOkO1QSohJYO8dqbz6zSbeWbydnDylf3Icd3ZrRpu4mMDHakwFUFQjtZUgTMVSKwGyvQyaO5JGg5qRPNGvDb+7tAWvfbuZKT9sZ/qPO7nk3Lrc2a05PVrFImLdY43xhb9XlKslItNEZJ2IrBWRi0TkPRFZ7nlsFZHlXs7dKiKrPMdZscD8yltDNgqfPQSHtlJ/6yf8ecMg1rgGsyrmYZrvnsGwt5bQ64Vvmbpkhw24M8YHfq1iEpHJwHxVfcOztnQ1VT1cYP9zwBFVHetw7lYgRVX3+3o9q2KqIlZOdQ+qK1iSCI2EximwYzHk5YC4QH9NAhoWxZLEMTyxpQ1rd2dQNzqcgSkJ3Ni5CU3qVgvCD2FM+VBUFZPfEoSI1ARWAM3V4SLiLudvBy5X1Q0O+7diCcJ4s3LqmQ3ZSYMgYxdM6AwnHcZTxCSgD6xi4cYD/Of7rcxZt4/cPKV7q1hu6tKEy8+vT2iIdewzVUuw2iCaA+nAWyKSDCwF7lfVXzz7uwF7nZKDhwIzRUSB11R1otNBIjIcGA7QpEmTsozflGdJg37t8VRQzTj3vE5OjuxAfkmna8v6dG1Zj0OL/kfI3Cepvn0vu7bW5cmwm6h14VBu6BhPQh0rVRjjzxJECrAIuERVF4vIi0CGqv7Zs/8VYKOqPufl/DhV3SUi9YFZwH2q+m1R17QShAHghbZwZIfzPlcYXHA11GkOi14+rZoqSyJ45OQdfJzblc5N63B9x8b0SWxEzciwAAVuTOAFa6qNNCBNVRd73k8DOngCCgUGAO95O1lVd3me9wHTgc5+jNVUJk6D7cKi4Mqx0Hk4bJoL8587ozdUhGbxXJ1P+ONV57H/lywe/WAVnZ6azcgpPzJv/T5bDtVUOX6rYlLVPSKyQ0TOU9X1QE/gJ8/uK4B1qprmdK6IRAMuVT3qed0LOKMh2xhH3gbb5W/v+Wd4uqHjqaFHd3LvZefy+x4t2DZvMrW++ys11+9j19q6jHbdSF7bgfwmqREXtahLmLVXmErO3+Mg7gPe9vRg2gwM82wfDEwpeKCIxAFvqGpfoAEw3dNfPRR4R1W/9HOspjLx1kYB7tJETIL3aqj3hyE1GtJ06VunShnxrv38RV7j8VXKLakXUataGFe1bkjfpEZc1Lwu4aGWLEzlYyOpTdXk1FU2JAKaXAR7V8HxA46n5dWMZ/ZVc5ixajez1+7jWFYONSJC6X5eLFdcUJ8erepTOzo8QD+EMWfPRlIbU1hR1VC52fBkLO6OdKdzZaTRq3kEvdq050R2Lhtmv0n8smeJ+Xkfu9bV5S+5g9iV0J+eF9TnsvPr07J+dRu5bSosK0EY46SonlASAudcDDUawtpPIefEqV3Zrkiej7yXVw52BKBBzQi6tYylW8t6XHJuPepVjwhE9Mb4LCgD5YLBEoQpM05VUGFRcPEDkHcS1n8B+35yPjcmgV3DljB/QzrHl77LVXteo6HuZ5fW450at5LTZiBdmtUhpWkdYqKsC60JLksQxpSGt9Ha+cbUwqkaCoBuD4Oqe3rynF+TzAki+FPOXXyYczEugQsa1aRLs7p0aV6Hzk3rWPuFCThLEMb4g7dqqJBwyMs9bS6ogvJqxrP4mm9ZvOUAIavfZ8ChN2mEu4Tx76hbOHzutXRoUpsO59SiZf0ahLisDcP4jyUIY/zBWzVUv/HQqjc80wSvJYxOd7pLGMvfPq0NI0si+At3807mhQBUjwilXUItEuNjSGocQ9vGMcTXjrKGb1NmrBeTMf5Q3IC8mHjnEkZoBKx413HOqAjN4uka0xh+7yiWbT/Esu2HqLlhOkN3TD5VyviT60bSEq4mKT6G1o1iuKBRDc6pG20lDVPmrARhjL8UVcJoez2MrYvXEkb1htAo2T131MaZkHvy1K6TEsE/o0Yw8XCKe9U81wIeDZtKIznAkbD6rDrvfsLaD6Zlg+rWa8oUy0oQxgRDaUsYkbWgxWWwe4VjT6lwzeKRvEk8cOvFpG9dQ4Mf3iQ0111NVTt7LymrxjBqWRqf5HWlTnQ459avzg2h39E3/XWiT+wht0ZjXFc8gSvZy0hzYzysBGFMsBRVwshPIkX1lCpCVlR9plz0Gev3n6Th1o8ZnvEiUfxaCsnUcF6sdh9b435Ds9houp+YS4cN4wn/ZTfExCOFe2yZSstKEMaUR8WVMMB7KaNmHFz3Gkzu5/jREZn7uO3rLu45p47vhQLJASBKTnLnyf8yaF93otZ/QHLI60SI55gjOzjx4Qgmf7OJ3ef0p3GtKBrXjqLN/q+I//FZXBk7EadYTaVjJQhjyrPiShneutpG1XFPbX5wE6x63/vnN0xE929EcjLP2LVPYumZ9xJHs3Lo71rAM2FvUE1+TTRZEsG0uEfYe04/GsZEkXxoJi1XP0/YsV1WCqlArARhTEVV7NTlo50TSJ+//3rM9kXOSSS8OtSIQ/ascrx0fU1nVcxD5FRvhGvvaly5p5dCIjSLK3a+zEWbz+dqWci1YW8QXqAUkvnhCF6Z9TPr6/emQc1IYqtH0DFjNu03jCfy+G5yazQm7/I/E95+8K8fWtzgRBNQVoIwpqIr7pdqaUshETXg/KshYyds8b6Yo7rCQPMQh4GBR101eSbifrZlRtEyew2PhL5PVIFSyHENZwx383305Vwb8h0jjo0nQrNO7c8JiWRNx6fIuuB6alcLI6ZaGHU2fUzo3CeL/nktyfjMBsoZU9UV9UvTl8Zyb0kkshak3A4Lni91aNkSwbIaPUg6+i1RemZV116txXVZYzlKNS53LeVvYZNOq+o6QQT/jX2I9fX7cNEvX9N/+zOE5f06+DA3JIqd3Z4ht+1AqkeEErNxOmHznkKKSiC+JN1KkoQsQRhjiuavUkiNRvDbt+H4fniniF+gRS3gVIACTsMBTxDBbNfF9Mj9nupy4oz9+7Umvz95PymudYwM/YhIyT7t3JdrjOTHWldSLTyEbplzGbT7WcILlmRckfzY7i8caH4t5+z8nFY/PEZI7q/3QkOjyP7NPwlr99tfR7mXRZIJQCIKWoIQkVrAG0Bb3P+2twNXAXcB6Z7D/qSqMxzO7Q28CITgXmnumeKuZwnCGD/yVykkJgEeXO19f7W60PMJyMqAmY97j69mY3d1WCnkIWRIDCcIp64eIIwzq8t+0Qim53blmpCF1HBIQge1Ok/k3Ia6wkkJ2cgQviCCXxPRScKZWut2VtTsQccTixmQ/vJpSSjbFcncVo+zvfFvCA91ce6eL+i8esypMS7gLg1tvfivHG01gLAQITzERXioi4jQEBrGRJbqZw9mgpgMzFfVNzzLjlYDHgCOqeq4Is4LAX4GrgTSgCXAEFX1Mr+ymyUIY4LobEsh/kwy0fVhwET477Xe4+84zH3tle867lYgN7IuIScOOJZiykKeCkeIJptQ6pBBqOSdcUxaXj26nhx/2rZ61cNJffzKUl0zKL2YRKQm0B24DUBVTwInfZxkrDOwUVU3ez7rXeAaoMgEYYwJoqLWAc/fD96TiC/jQrz12uo5uuj9Vz3tHp3urSorJgH6/dP9ettCx2MkJoHQopJQjUZwy8eQkwWvdcfrAMd+L8Kn9zvuElGqdfgteTknCVn5X8djGrsO8NZtncjKySM71/3w1zxc/uzm2hx3NdJbIpIMLAXy78oIEbkFSAUeVtVDhc5tDBT8F0gDujhdRESGA8MBmjRpUnbRG2PKni9JxJ9JprgE48sx3vZfORZiz3O/9zbAMSYBOt4G347zmoQi+nsa/LfN83JMPJedX9/5/pQxlx8/OxToALyiqu2BX4BRwCtAC6AdsBt4zuFcp3TomI5VdaKqpqhqSmxsbJkEbowpx5IGuauTxhx2PxdOKEXtTxrkrrKKSQDEU3IYX7JjfPmMnqPdSaOgwkmmqP2+HuNn/ixBpAFpqrrY834aMEpV9+YfICKvA595OTehwPt4YJe/AjXGVCHFlVJ8OSYQ1Wm+HONn/m6kng/cqarrRWQMEA08r6q7PfsfBLqo6uBC54XibqTuCezE3Uh9o6quKep61khtjDElE8ypNu4D3vb0YNoMDAPGi0g73FVGW4G7PUHG4e7O2ldVc0RkBPAV7m6ubxaXHIwxxpQtGyhnjDFVWFElCH82UhtjjKnALEEYY4xxZAnCGGOMo0rVBiEi6cC2Up5eD9hfhuH4i8VZ9ipKrBZn2aoocYJ/Yz1HVR0HkVWqBHE2RCTVW0NNeWJxlr2KEqvFWbYqSpwQvFitiskYY4wjSxDGGGMcWYL41cRgB+Aji7PsVZRYLc6yVVHihCDFam0QxhhjHFkJwhhjjCNLEMYYYxxVqQQhIr1FZL2IbBSRUQ77RUTGe/avFJEOQYozQUTmishaEVkjImcsPyUiPUTkiIgs9zwCN0n86XFsFZFVnhjOmAirPNxTETmvwH1aLiIZIvJAoWOCdj9F5E0R2SciqwtsqyMis0Rkg+e5tpdzi/xOByDOZ0VkneffdrpnHXqnc4v8ngQgzjEisrPAv29fL+cG7H4WEet7BeLcKiLLvZzr/3uqqlXigXtW2E24V7oLB1YArQsd0xf4AveCRRcCi4MUayOgg+d1DdxTnxeOtQfwWTm4r1uBekXsLxf3tND3YA/uwUHl4n7iXpq3A7C6wLZ/4F4/BdwLbf3dy89S5Hc6AHH2AkI9r//uFKcv35MAxDkG+IMP342A3U9vsRba/xwwOlj3tCqVIE6tc63u9bHz17ku6BrgP+q2CKglIo0CHaiq7lbVZZ7XR4G1uJdhrYjKxT0toCewSVVLO+K+zKnqt8DBQpuvASZ7Xk8GrnU41ZfvtF/jVNWZqprjebsI9+JeQeXlfvoioPcTio5VRAQYBEzxZwxFqUoJwmmd68K/dH05JqBEpCnQHljssPsiEVkhIl+ISJuABvYrBWaKyFLP+uCFlbd7Ohjv/+HKw/3M10A9C2t5np0WIS5v9/Z23KVFJ8V9TwJhhKcq7E0vVXbl7X52A/aq6gYv+/1+T6tSgvBlnWuf18IOBBGpDnwAPKCqGYV2L8NdTZIM/Av4KNDxeVyiqh2APsC9ItK90P5yc0/FvXBVf+B9h93l5X6WRHm6t48BOcDbXg4p7nvib68ALYB2wG7cVTeFlZv76TGEoksPfr+nVSlB+LLOdblZC1tEwnAnh7dV9cPC+1U1Q1WPeV7PAMJEpF6Aw0RVd3me9wHTcRfTCyo39xT3f6RlWmBd9Hzl5X4WsDe/Ks7zvM/hmHJxb0XkVuBqYKh6KscL8+F74lequldVc1U1D3jdy/XLxf2EU8suDwDe83ZMIO5pVUoQS4CWItLM85fkYOCTQsd8Atzi6XlzIXAkv5gfSJ66x0nAWlV93ssxDT3HISKdcf9bHghclCAi0SJSI/817gbL1YUOKxf31MPrX2Tl4X4W8glwq+f1rcDHDsf48p32KxHpDTwK9FfV416O8eV74leF2r2u83L9oN/PAq4A1qlqmtPOgN1Tf7aAl7cH7h41P+PuqfCYZ9s9wD2e1wK85Nm/CkgJUpxdcRdtVwLLPY++hWIdAazB3dNiEXBxEOJs7rn+Ck8s5fmeVsP9Cz+mwLZycT9xJ63dQDbuv2LvAOoCc4ANnuc6nmPjgBlFfacDHOdG3PX2+d/TVwvH6e17EuA4/+v5/q3E/Uu/UbDvp7dYPdv/nf/dLHBswO+pTbVhjDHGUVWqYjLGGFMCliCMMcY4sgRhjDHGkSUIY4wxjixBGGOMcWQJwhgHIlJLRH5f4H2ciEzz07WuzZ89VkT+LSI3OBwTKyJf+uP6xnhjCcIYZ7WAUwlCVXep6hm/uMvII8DLRR2gqunAbhG5xE8xGHMGSxDGOHsGaOGZa/9ZEWmaP2e/iNwmIh+JyKciskVERojIQyLyo4gsEpE6nuNaiMiXnsnU5ovI+YUvIiKtgCxV3V9gc3cR+U5ENhcqTXwEDPXjz2zMaSxBGONsFO5pwdup6h8d9rcFbsQ9/83TwHFVbQ98D9ziOWYicJ+qdgT+gHMp4RLcEwUW1Aj3aPqrcSeqfKm4Z/g0JiBCgx2AMRXUXHWv1XFURI4An3q2rwKSPDPxXgy875niCSDC4XMaAemFtn2k7knlfhKRBgW278M93YIxAWEJwpjSySrwOq/A+zzc/69cwGFVbVfM52QCMUV8dsEpqCM9xxsTEFbFZIyzo7iXey0Vda/fsUVEBsKptbmTHQ5dC5zr48e2IsCzoJqqzRKEMQ5U9QCwUERWi8izpfyYocAdIpI/46bT8pXfAu2lQD1UES4DPi9lLMaUmM3makyQiciLwKeqOruY474FrlHVQ4GJzFR1VoIwJvj+inu9Cq9EJBZ43pKDCSQrQRhjjHFkJQhjjDGOLEEYY4xxZAnCGGOMI0sQxhhjHFmCMMYY4+j/AXprH+Sc0yoLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "t=np.linspace(0,18)\n",
    "t_1=np.arange(0,18,time_step)\n",
    "Temp_ambient = 65\n",
    "Temp_t1 = 85\n",
    "plt.plot(t,T_analytical(Temp_ambient,Temp_t1,K,t),'-',label='Analytical')\n",
    "plt.plot(t_1,T_numerical,'o-',label='Numerical')\n",
    "plt.legend()\n",
    "plt.xlabel('time (h)')\n",
    "plt.ylabel('Temperature (F)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Answer 3b:\n",
    "As shown in the graph above, as time goes to infinity the solutions converge at 65 degrees"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Answer 3c:\n",
    "The time where the temperature was aproximately 98.6 degrees was -1 hours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Temperature (F)')"
      ]
     },
     "execution_count": 254,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU9b3/8ddnJntCEhISSAhLQHYQ0ACKCu6ireJ+sdTt+hN7q13vbUtvvWrXa/fqw9perbhUxVpbl7origgiCuLCKsgaQBIISyBkne/vj5nEAFmGJJMzybyfj8c8zpnvnJnzyTwC75zz/Z7vMeccIiIiAD6vCxARkeihUBARkQYKBRERaaBQEBGRBgoFERFpEOd1Ae3Rq1cvN3DgQK/LEBHpUpYtW7bLOZfT1GtdOhQGDhzI0qVLvS5DRKRLMbPNzb2m00ciItJAoSAiIg0UCiIi0qBL9ymISPdSU1NDcXExlZWVXpfSLSQlJVFQUEB8fHzY71EoiEjUKC4upkePHgwcOBAz87qcLs05x+7duykuLqawsDDs9+n0kYhEjcrKSrKzsxUIHcDMyM7OPuajLoWCiEQVBULHact3GZOhsG3vIX758hp27DvkdSkiIlElYqFgZnPMrMTMVjRqu8LMVppZwMyKjtj+h2a23szWmtl5kaoL4GBVLX+a/xnz15ZGcjci0kU9/fTTmBlr1qxp82dcd911PPXUUy1u84tf/OKw55MnT27Tvu644w5+85vftOm9R4rkkcJDwLQj2lYAlwILGjea2UhgBjAq9J57zcwfqcKG5KaRl5HEWwoFEWnC3LlzOfXUU3niiSciup8jQ+Gdd96J6P7CEbFQcM4tAMqOaFvtnFvbxObTgSecc1XOuY3AemBipGozM6YOzWHR+l3U1AUitRsR6YIOHDjAokWLeOCBBxpCYf78+Zx++ulcfvnlDB8+nJkzZ1J/18qf/OQnTJgwgdGjRzNr1iyOvJvlvHnzuOSSSxqev/baa1x66aXMnj2bQ4cOMW7cOGbOnAlAWlpaw3a/+tWvGDNmDGPHjmX27NkA3H///UyYMIGxY8dy2WWXUVFR0eE/f7QMSe0LvNvoeXGo7ShmNguYBdC/f/827/D0YTk88f5Wlm/Zy8TCrDZ/johExo//tZJV2/d36GeOzE/n9gtHtbjNM888w7Rp0xg6dChZWVl88MEHACxfvpyVK1eSn5/PKaecwqJFizj11FO55ZZbuO222wC4+uqref7557nwwgsbPu/MM8/k5ptvprS0lJycHB588EGuv/56LrzwQu655x4+/PDDo2p46aWXeOaZZ1iyZAkpKSmUlQX/vr700ku58cYbAbj11lt54IEH+MY3vtEh3029aOlobqqLvMmbRzvn7nPOFTnninJympzkLyyTj+uF32e89WlJmz9DRLqfuXPnMmPGDABmzJjB3LlzAZg4cSIFBQX4fD7GjRvHpk2bAHjzzTeZNGkSY8aM4Y033mDlypWHfZ6ZcfXVV/Poo4+yd+9eFi9ezPnnn99iDa+//jrXX389KSkpAGRlBf9wXbFiBaeddhpjxozhscceO2pfHSFajhSKgX6NnhcA2yO5w/SkeE7s35O3Pi3le+cNj+SuRKQNWvuLPhJ2797NG2+8wYoVKzAz6urqMDMuuOACEhMTG7bz+/3U1tZSWVnJ17/+dZYuXUq/fv244447mrwuoP7IICkpiSuuuIK4uJb/63XONTmc9LrrruOZZ55h7NixPPTQQ8yfP7/dP/ORouVI4TlghpklmlkhMAR4L9I7nToshxXb9lNSrkvqRQSeeuoprrnmGjZv3symTZvYunUrhYWFLFy4sMnt6wOgV69eHDhwoNnRRvn5+eTn5/Ozn/2M6667rqE9Pj6empqao7Y/99xzmTNnTkOfQf3po/LycvLy8qipqeGxxx5rz4/arEgOSZ0LLAaGmVmxmd1gZpeYWTFwMvCCmb0C4JxbCTwJrAJeBm52ztVFqrZ6U4cGTz+9/emuSO9KRLqAuXPnHtYpDHDZZZfx+OOPN7l9ZmYmN954I2PGjOHiiy9mwoQJzX72zJkz6devHyNHjmxomzVrFscff3xDR3O9adOmcdFFF1FUVMS4ceMahpv+9Kc/ZdKkSZxzzjkMHx6ZMxx2ZE95V1JUVOTac5OdQMAx8RevM3lwL+6+anwHViYibbF69WpGjBjhdRkRccsttzB+/HhuuOGGTt1vU9+pmS1zzhU1tX209Cl4wuczpgzJ4c21JdQFHH6fLq8XkY534oknkpqaym9/+1uvS2lVtPQpeGbqsBz2VNTwybZ9XpciIt3UsmXLWLBgwWGd1dEq5kPhtCE5mKGrm0WiRFc+pR1t2vJdxnwoZKUmcHxBpq5XEIkCSUlJ7N69W8HQAervp5CUlHRM74vpPoV6U4fmcM8b69hbUU1mSoLX5YjErIKCAoqLiykt1ZF7R6i/89qxUCgQDIW7561j4fpdfPn4fK/LEYlZ8fHxx3SXMOl4MX/6CGBsQQYZyfGaSltEYp5CAYjz+zh1SC/e+rRU5zJFJKYpFEKmDs2htLyK1TvKvS5FRMQzCoWQ00NTXrz1qU4hiUjsUiiE5KYnMSIvXUNTRSSmKRQamTo0h6Wb9nCgqtbrUkREPKFQaGTq0BxqA4531mvWVBGJTQqFRk4c0JO0xDjmrdYpJBGJTQqFRhLifJwxPJfXV++kLqChqSISexQKRzhvVG92H6xm6aYyr0sREel0CoUjnD4sl4Q4H6+s3Ol1KSIinU6hcIS0xDhOPa4Xr6z8XFc3i0jMUSg0YdqoPmzbe4iV2/d7XYqISKdSKDThrBG5+AxeXfm516WIiHQqhUITstMSmTAwS/0KIhJzFArNOG9UH9buLGfjroNelyIi0mkUCs04d1RvAF7RKSQRiSEKhWYU9ExhdN90hYKIxBSFQgvOG9mH5Vv2snN/pdeliIh0CoVCC6aN7gPAq6vU4SwisUGh0ILjctMY1CtVQ1NFJGYoFFpgZpw7qg+LP9vNvooar8sREYk4hUIrzhvVm9qAY94anUISke5PodCKsQWZ9E5P1CgkEYkJCoVW+HzGuSP78NanpRyqrvO6HBGRiFIohGHa6D5U1gRYsK7U61JERCJKoRCGiYVZZCTH88oKnUISke5NoRCGeL+P80b15tVVO6ms0SkkEem+IhYKZjbHzErMbEWjtiwze83M1oWWPRu99kMzW29ma83svEjV1VYXj+vLgapaXl+tUUgi0n1F8kjhIWDaEW2zgXnOuSHAvNBzzGwkMAMYFXrPvWbmj2Btx2zSoGx6pyfyzPJtXpciIhIxEQsF59wCoOyI5unAw6H1h4GLG7U/4Zyrcs5tBNYDEyNVW1v4fcb0cX2Zv7aUPQervS5HRCQiOrtPobdzbgdAaJkbau8LbG20XXGo7ShmNsvMlprZ0tLSzh0NNH1cPrUBxwuf7OjU/YqIdJZo6Wi2JtpcUxs65+5zzhU554pycnIiXNbhRualMyQ3TaeQRKTb6uxQ2GlmeQChZUmovRjo12i7AmB7J9fWKjPj4vF9Wbp5D1vLKrwuR0Skw3V2KDwHXBtavxZ4tlH7DDNLNLNCYAjwXifXFpbp4/IBeO6jqMssEZF2i+SQ1LnAYmCYmRWb2Q3AncA5ZrYOOCf0HOfcSuBJYBXwMnCzcy4qLwgo6JnChIE9eXr5Npxr8gyXiEiXFRepD3bOXdXMS2c1s/3PgZ9Hqp6OdPH4vvzo6RWs3L6f0X0zvC5HRKTDREtHc5fypTF5xPuNZz9Uh7OIdC8KhTbITElg6tBcnv1wO3UBnUISke5DodBGl4zvS0l5Fe9u2O11KSIiHUah0EZnjcglLTFO1yyISLeiUGijpHg/00b34aUVn2vmVBHpNhQK7XDJ+ODMqfNWl7S+sYhIF6BQaIeTBmWT2yORZzQKSUS6CYVCO/h9wWkv3lxTQml5ldfliIi0m0Khna4s6kdtwPH3ZVtb31hEJMq1eEWzmSUAFwCnAfnAIWAF8KJzbk3ky4t+x+WmMakwiyfe28rXpgzG52tqwlcRka6h2SMFM7sVWAKcAXxE8KY4zxEMkt+b2ctmNrpTqoxyX5nUny1lFSz6bJfXpYiItEtLRwqfOOd+1sxrvwpNfd2vmddjyrTRfeiZEs/jS7Zw2pDOvceDiEhHaqlP4bmW3uic2+Gci8rprTtbYpyfy08s4LVVOykpr/S6HBGRNmspFJbVr5jZHzqhli7tqon9gx3OS4u9LkVEpM1aCoXGPaZTIl1IVzcoJ42TB2XzxPtbCGiSPBHpoloKBf3PdoyumtSfrWWHeHu9OpxFpGtqqaN5uJl9QPCIYVhondBz55w7IeLVdTHnjepNVmoCc5dsYepQdTiLSNfTUiiM6bQquonEOD9XnFjAXxZupGR/JbnpSV6XJCJyTJo9feSc+6ylR2cW2ZXMmNifuoDjyaW6wllEup6WLl5708z+w8zyj2iPM7MpZvaAmV0f+RK7lsJeqUwenM3c97aqw1lEupyWOpq/BMQDT5tZsZl9bGafAhuA64E/Oece7Iwiu5qvTOrPtr2HWLCu1OtSRESOSbN9Cs65CuBu4G4zSwRygUPOOQ2tacW5I/uQnZrA40u2cPqwXK/LEREJW1izpDrnqpxzWxUI4UmI83F5UQHz1pSwY98hr8sREQmbps6OkK9OGoBzjofe2eR1KSIiYVMoREi/rBTOH5PH40u2cKCq1utyRETCElYomFmBmZ0RWk80s9TIltU9zDptEOWVtfztfQ1PFZGuodVQMLN/Jzhj6l9CTQOAZyNZVHcxtl8mEwdmMWfhRmrrAl6XIyLSqnCOFL4JnATsB3DOfUpwJJKE4cYpg9i29xAvrvjc61JERFoVTihUOueq65+YmZ/DZ1CVFpw1PJdBvVK5f8EGnNPFbCIS3cIJhUVm9n0gKdSv8Dfg+ciW1X34fMYNpxXyybZ9LNlY5nU5IiItCicUvg+UA2uAbwHzgB9Fsqju5rITCshKTeAvb2/wuhQRkRa1NEtq/amiOc65a4E/dU5J3U9SvJ9rTh7AH15fx/qSAxyXm+Z1SSIiTWrxSME5VwfkmVl8J9XTbV190gAS43w8sFBHCyISvVo8UgjZALxtZs8CB+sbnXN3R6yqbig7LZHLTizgqWXF/Oe5w+iVluh1SSIiRwmnT6EUeA1IAXIaPeQY3XBqIdW1AR5ZvNnrUkREmtTqkYJz7n86eqdm9i3gRoJDW+93zv3BzLIIjmwaCGwCrnTO7enofXtpcE4aZ4/ozaPvbuY/pg4mOcHvdUkiIocJ54rm18zs1SMfbd2hmY0mGAgTgbHAl81sCDAbmOecG0JwhNPstu4jms2aMoiyg9X87f0tXpciInKUcPoUbm20ngRcBlS1Y58jgHdD92vAzN4CLgGmA6eHtnkYmA/8oB37iUoTBvZkUmEW987/jBkT+5MUr6MFEYkerR4pOOeWNHq85Zz7JsG/8ttqBTDFzLLNLAW4AOgH9HbO7QjtcwfNTKVhZrPMbKmZLS0t7Xp3NjMzvnPOUErKq3hsiY4WRCS6hHP6KL3RI9PMzgLy2rpD59xq4JcEO69fBj4Cwp5b2jl3n3OuyDlXlJPTNfu7TxqUzSnHZfOn+eupqNa02iISPcIZfbSS4F/3K4HlBK9mvrE9O3XOPeCcO8E5NwUoA9YBO80sDyC0LGnPPqLdd84eyq4D1Tz6rkYiiUj0CCcUBjnn+jvn+jnnCp1zZwKL2rNTM8sNLfsDlwJzCU7PfW1ok2vp5tNzFw3MYsrQHP781gYO6iY8IhIlwgmFJU20vdfO/f7DzFYB/wJuDg09vRM4x8zWAeeEnndr3zl7CGUHq3l48SavSxERAVoYfRT6az4PSDazMXwxXXY6wQvZ2sw5d1oTbbuBs9rzuV3N+P49OXN4Lvct2MDVJw2gR5JmExERb7V0pPAl4B6gALgX+GPo8d9Ah1/QFqu+c/ZQ9lbU8NCiTV6XIiLS/JGCc+5B4EEzu9I592Qn1hRTxhRkcM7I3tz/9gaumTyQjGQdLYiId8K5TuFJMzvPzL5rZv9d/+iM4mLFt88ewv7KWuYs3Oh1KSIS48K5TuFegqOBvgskA18FjotwXTFlVH4G54/uw5yFG9lbUd36G0REIiSc0UenOue+AuwOTY43iWA/g3Sgb589lAPVtfz5Ld1vQUS8E04oVNYvzaxP6PnAiFUUo4b16cEl4/oyZ9FGtpZVeF2OiMSocELhRTPLBH4DfEhwWuunIllUrPretGH4zbjzpTVelyIiMarFUDAzH/CSc26vc+7vQCEwxjmnjuYIyMtI5mtTB/PCJztYsmG31+WISAxq7R7NAeCuRs8POefKIl5VDJs1ZRD5GUn85PlV1AWc1+WISIwJ5/TRa2Y2PeKVCADJCX5+cP5wVm7fzz8+KPa6HBGJMeGEwi3A02Z2yMzKzGyPmeloIYIuGpvPCf0z+fUrazmgyfJEpBOFEwq9gHggDcgJPe+aNzLoIsyM2y4cRWl5Ffe+ud7rckQkhoRzRXMdcAXwg9B6HjAu0oXFunH9Mrl0fF/+slBDVEWk84RzRfM9wBnA1aGmCuDPkSxKgjREVUQ6WzinjyY7524idBFbaPRRQkSrEkBDVEWk84UTCjWh6xUcgJllA4GIViUN6oeo3v7cSmrq9LWLSGSFEwp/BP4B5JjZj4GFwC8jWpU0SE7wc/tFo1jzeTn3v615kUQkspq9n0I959wjZrYMODvUdIVzbkVky5LGzhvVh/NH9+EPr69j2qg+DMpJ87okEemmwjlSAPADNUD1MbxHOtCPp48iKc7HD//5CQFd6SwiERLO6KMfAXOBfIJTZj9uZj+MdGFyuNweSfzoSyNYsrGMJ97f6nU5ItJNhfNX/1eBCc65W51zPwImAtdEtixpypVF/Zg8OJv/fXE1n++rbP0NIiLHKJxQ2MzhfQ9xgHo8PWBm/O+lY6iuC/A/z67AOZ1GEpGOFU4oVAArzewvZnY/8Amw18x+Z2a/i2x5cqQB2al895yhvLZqJy+t+NzrckSkm2l19BHwQuhR790I1SJhuuHUQp7/eAe3PbuSUwb3IiMl3uuSRKSbCGdI6gOdUYiEL87v487LxnDRPYv4+Yur+NXlY70uSUS6iXBGH00zs/fNrERTZ0ePUfkZ3DRlEE8uLeaVlTqNJCIdI5w+hXuAm4C+aOrsqPLts4cypm8GP/jHx+zYd8jrckSkGwgnFIqBD51zNc65uvpHpAuT1iXE+bj7qvFU1wb4zt8+1O07RaTdwgmF7wP/MrPvmdk36x+RLkzCU9grlR9fNIp3N5Tx57c+87ocEeniwgmFHwN1QCbB00b1D4kSl59YwIVj8/nda5+ybPMer8sRkS4snCGpuc65EyNeibSZmfHzS0azfMsevvXEcl781mmkJ2mYqogcu3COFOaZ2ZkRr0TaJT0pnrtmjGfHvkpufVpXO4tI24QTCjcCr5vZAQ1JjW4nDujJt88awnMfbecfH2zzuhwR6YLCCYVeQDyQgYakRr2vn3EckwqzuO3ZFawvKfe6HBHpYloNhdDw0yuAH4TW84Bx7dmpmX3HzFaa2Qozm2tmSWaWZWavmdm60LJne/YRq/w+4w8zxpGS4OfGR5axr6LG65JEpAsJ54rme4AzgKtDTRXAn9u6QzPrC3wTKHLOjSZ4A58ZwGxgnnNuCDAv9FzaIC8jmT999USK91TwjSeW6/oFEQlbOKePJjvnbgIqAZxzZUBCO/cbBySbWRyQAmwHpgMPh15/GLi4nfuIaRMGZvGT6aNZ8Gkpv3x5jdfliEgXEU4o1JiZD3AAZpYNBNq6Q+fcNuA3wBZgB7DPOfcq0Ns5tyO0zQ4gt6n3m9ksM1tqZktLS0vbWkZMuGpif645eQD3LdjAPz8o9rocEekCmg2F0F/xAH8E/gHkmNmPgYXAL9u6w1BfwXSgkOAtPlPN7Kvhvt85d59zrsg5V5STo/7u1vzPl0dy0qAsZv/zEz7autfrckQkyrV0pPAegHPuEeBWgn/d7wGucM490Y59ng1sdM6VOudqgH8Ck4GdZpYHEFqWtGMfEhLv93HvzBPJSUtk1l+XUrJft/EUkea1FApWv+KcW+mcu8s59wfn3Ip27nMLcJKZpZiZAWcBq4HngGtD21wLPNvO/UhIVmoC919TxP5DtXzt0WVU1Wo+QxFpWkvTXOSY2Xebe9E516ZbcTrnlpjZU8AHQC2wHLgPSAOeNLMbCAbHFW35fGnayPx0fnvlWL7+2Af855MfcfeM8fh81vobRSSmtBQKfoL/UXf4/xzOuduB249oriJ41CARcsGYPH54/nD+96U1ZKcmcMdFowgerImIBLUUCjuccz/ptEqkU9w0dTC7D1Zz34INZKUm8q2zh3hdkohEkZZCQX9CdlM/PH84uw9U8/vXPyUrLYGrTxrgdUkiEiVaCgWdyummzIxfXjaGvRXV3PbsCrJSEvjS8XlelyUiUaDZ0UehK5elm4rz+/jjzBMoGtCTb/9tOQvX7fK6JBGJAuFc0SzdVFK8n79cO4HBOWnc9NelurhNRBQKsS4jOZ5H/n0iWWkJfPWBJSzfott5isQyhYKQm57EE7NOJis1ga/+ZQnvbdSZQ5FYpVAQAPpmJvPkTSfTJyOJa+e8xzvr1ccgEosUCtKgd+iIoX9WCtc/9D7z12r6KZFYo1CQw+T0SGTurJMYnJPGrEeW8fqqnV6XJCKdSKEgR8lKTWDujScxIq8HX3t0GS98vMPrkkSkkygUpEkZKfE8+v8mMa5fJrfM/YAHFm7EOd3WU6S7UyhIs3okxfPXGyZx3sg+/PT5Vdzx3Epq69p80z0R6QIUCtKi5AQ/9848gVlTBvHw4s3M+usyDlbVel2WiESIQkFa5fMZ/33BCH568Wjmry3hyv9bzE7dwU2kW1IoSNiuPmkAD1w3gU27DnLxHxexesd+r0sSkQ6mUJBjcsawXP7+tck4B5f/6R2e/3i71yWJSAdSKMgxG5mfzjM3n8KwPj245fHl3PHcSqpr1QEt0h0oFKRN+mQEr37+91MKeeidTfzbfYvZvveQ12WJSDspFKTNEuJ83HbhSO6deQLrdh7gS3e/zVuflnpdloi0g0JB2u2CMXk8d8sp9E5P4roH3+P3r31KXUAXuol0RQoF6RCDctJ4+uuncMn4vtw1bx1X/t9iNu066HVZInKMFArSYZIT/Pz2irH84d/GsW5nOeff9TaPLN5EQEcNIl2GQkE6lJlx8fi+vPqdqUwozOK2Z1dyzZz31Akt0kUoFCQi+mQk8fD1E/j5JaP5YMsezvv9Ap5aVqxJ9USinEJBIsbMmDlpAC9/awoj8tL5r79/xLUPvs9G9TWIRC2FgkRc/+wU5s46idsvHMkHm4NHDb97dS2VNXVelyYiR1AoSKfw+4zrTynkjf+cyvlj+nD3G+s55/dvMW+17uwmEk0UCtKpctOTuGvGeB6/cRKJcX5ueHgp/+/hpWwtq/C6NBFBoSAemTy4Fy9+8zRmnz+cRet3cdZv3+In/1pF2cFqr0sTiWkKBfFMQpyPr00dzJv/dTqXntCXh97ZyNRfvck9b6yjolo38hHxgnXlIYJFRUVu6dKlXpchHWR9STm/enktr67aSU6PRL599hCuLOpHvF9/u4h0JDNb5pwravI1hYJEm2Wby7jzpTW8v2kPA7JT+I+pg7n0hAIS4hQOIh1BoSBdjnOOeatLuGveOj7Zto+8jCRuPG0QV03sT3KC3+vyRLo0hYJ0Wc453l63i3veXM97G8vISk3ghlMLufrkAaQnxXtdnkiXFFWhYGbDgL81ahoE3AY8EmofCGwCrnTO7WnpsxQKseX9TWX88c31zF9bSo/EOC4vKuCakwdS2CvV69JEupSoCoXDdm7mB7YBk4CbgTLn3J1mNhvo6Zz7QUvvVyjEphXb9nH/2xt48ZMd1NQ5pg7N4brJA5k6NAefz7wuTyTqRXMonAvc7pw7xczWAqc753aYWR4w3zk3rKX3KxRiW0l5JY8v2cJjS7ZQWl7FwOwUrj55IJed0JfMlASvyxOJWtEcCnOAD5xz95jZXudcZqPX9jjnejbxnlnALID+/fufuHnz5s4rWKJSdW2Al1d+zsPvbGLZ5j0k+H2cM7I3lxcVMGVIDn4dPYgcJipDwcwSgO3AKOfcznBDoTEdKciRVm3fz9+XbeXZD7dTdrCa3umJXDK+gCuKChick+Z1eSJRIVpDYTpws3Pu3NBznT6SDlNdG+CNNTt5alkxb64tpS7gOL4ggy+NyeOCMXn0y0rxukQRz0RrKDwBvOKcezD0/NfA7kYdzVnOue+39BkKBQlHSXklzy7fzvMfb+ej4n0AjO2XyZfH5HHB8Xn0zUz2uEKRzhV1oWBmKcBWYJBzbl+oLRt4EugPbAGucM6VtfQ5CgU5VlvLKnjhkx288PEOPtkWCoiCDM4a0Zszh+cyKj8dM/VBSPcWdaHQURQK0h6bdx/k+Y938NqqnXxUvBfnoE96EmeOyOXsEblMHtyLpHhdPS3dj0JBpBWl5VW8ubaEN1aX8Pa6Ug5W15EY52PCwCxOOa4Xpx7Xi5H56RrJJN2CQkHkGFTV1vHuhjIWfFrKovW7WPN5OQAZyfFMHpzN5ON6Makwi+Ny0nSxnHRJLYVCXGcXIxLtEuP8TB2aw9ShOUCwo3rxZ7tZuG4Xi9bv4qUVnwOQmRJP0YCeTBiYRdHALMb0zdBMrtLlKRREWpHbI4np4/oyfVxfnHNsKavgvY1lLN20h/c3lfH66hIAEuN8jO6bwfEFGYzrl8nYgkwGZKeo41q6FJ0+EmmnXQeqWLqpjPc37eGjrXtZsX0flTUBIHjK6fiCDEb3zWBkXjoj89MZmJ2qvgnxlPoURDpRbV2AT3ce4KPivXxcvJcPt+5j3c5yagPBf2vJ8X6G5/VgRF46I/r0YEjvHgzJTSM7LdHjyiVWKBREPFZVW8f6kgOs2r6fVTv2s3rHflZt38/+yi/uRZ2VmsCQ3DSG9E5jSG4PCnulUtgrlfzMZB1ZSIdSR7OIxxLj/IzKz2BUfkZDm3OOzxeG1rUAAAgASURBVPdXsm7nAdaVHGDdznLWlRzg2Q+3U94oLBLifAzISmkIiX5ZKfQPPfIzk9W5LR1KoSDiETMjLyOZvIxkpoRGOkEwLErLq9i46+BRj/lrS6muCzRs6zPIy0imoGcyfXsm0zcz+Hn5mUnB9cxk0hL1z1zCp98WkShjZuSmJ5GbnsSkQdmHvRYIOHaWV7K17BBbyirYUlZBcWj57me7+Xx/JYEjzgj3SIyjd0YSvdMT6Z2eRO/0JPqkJ5HbI5GcHon0SgsuUxUegkJBpEvx+b44uphYmHXU67V1AUrKq9i+9xDb9h5i+95Kdu6v5PN9lewsr+Tdz3ZTUl7V0OndWHK8PxQSCWSlJpKdmkDP1ASyUxPISk0gKy2BnikJ9EyJJzM5gR5Jcbp4rxtSKIh0I3F+H/mZyeRnJtNkLyLBo43dB6spKa9k14FqSsur2HWgqmG560AVxXsq+Lh4L2UHq5sMEAieuspIjqdnSgLpyfFkJMeHlnGkJ33xPC0xjh5JcfRIig8t40hLjCM1QaESjRQKIjHG5zNyQqeOWuOcY39lLWUHqyk7WMXeihr2VNSwt6I6tB5c7q8Mrm/efZD9lbXsO1RDXTNh0lhKgp/UxFBIJPpJTYgjNTGOlAR/6BHXsE1SfLAtOd7/xXqCn6Q4P8kJPhLjgu1J8T6S4v3E+9UB3xYKBRFplpmREToKKOyVGvb7nHMcrK6jvLKG8sra0OOL9QNVNRyoquNgVS0Hq2o5EFoerKqjpLySiqo6KqrrOFhdS0V1XVgBcyS/z0iMCwZEYpwv9PCTGO8jwe9rWCbE+UiI8zesJ8b5iPcb8aHn8f7gdvF+I65+Pc6I832xXZzfR7wv+Hqc34j3BZdx9W0+I85v+H3B9wWXdtgyWq58VyiISIczM9JCRwB5Ga1v3xLnHNV1ASqrA1TU1HKouo5DNXUcqg4GR2VNHZW1ASpr6qiqqaOyJsChmjqqauuoqglQVf9abSDYVhugujZAZU2A/Ydqqa4NUF0XoKqmjuo6R01d8PWaukCzp84iwe8z/BYMiKMeoXafD/xm+HzGmcNyufXLIzu8DoWCiEQ1Mwv+hR/nJ4P4Tt13IBAMpJq6ALX1gdFovSa0rA0E12vrHDWB4Ot1oba6QHCbuoCjJuAINHpeG3ANy0D90rmG99c5R12A4HoAAi64fZ1z5EXojoEKBRGRZvh8RpLPH1M3W1JPjIiINFAoiIhIA4WCiIg0UCiIiEgDhYKIiDRQKIiISAOFgoiINFAoiIhIgy59O04zKwU2e11HJ+kF7PK6iCim76d1+o5aFkvfzwDnXE5TL3TpUIglZra0uXuqir6fcOg7apm+nyCdPhIRkQYKBRERaaBQ6Dru87qAKKfvp3X6jlqm7wf1KYiISCM6UhARkQYKBRERaaBQiHJmNs3M1prZejOb7XU90cjMNpnZJ2b2oZkt9bqeaGBmc8ysxMxWNGrLMrPXzGxdaNnTyxq91Mz3c4eZbQv9Hn1oZhd4WaNXFApRzMz8wB+B84GRwFVm1vE3Ze0eznDOjdM48wYPAdOOaJsNzHPODQHmhZ7Hqoc4+vsB+H3o92icc+7FTq4pKigUottEYL1zboNzrhp4ApjucU3SBTjnFgBlRzRPBx4OrT8MXNypRUWRZr4fQaEQ7foCWxs9Lw61yeEc8KqZLTOzWV4XE8V6O+d2AISWuR7XE41uMbOPQ6eXYvL0mkIhulkTbRpDfLRTnHMnEDzNdrOZTfG6IOmS/gQMBsYBO4DfeluONxQK0a0Y6NfoeQGw3aNaopZzbntoWQI8TfC0mxxtp5nlAYSWJR7XE1Wcczudc3XOuQBwPzH6e6RQiG7vA0PMrNDMEoAZwHMe1xRVzCzVzHrUrwPnAitaflfMeg64NrR+LfCsh7VEnfrADLmEGP09ivO6AGmec67WzG4BXgH8wBzn3EqPy4o2vYGnzQyCv8+PO+de9rYk75nZXOB0oJeZFQO3A3cCT5rZDcAW4ArvKvRWM9/P6WY2juAp2k3ATZ4V6CFNcyEiIg10+khERBooFEREpIFCQUREGigURESkgUJBREQaKBREQsws08y+3uh5vpk9FaF9XWxmt4XWHzKzy5vYJsfMYn54rXQuhYLIFzKBhlBwzm13zh31n3UH+T5wb0sbOOdKgR1mdkqEahA5ikJB5At3AoNDc+n/2swG1s+3b2bXmdkzZvYvM9toZreY2XfNbLmZvWtmWaHtBpvZy6HJ+d42s+FH7sTMhgJVzrldjZqnmNk7ZrbhiKOGZ4CZEfyZRQ6jUBD5wmzgs9Bc+t9r4vXRwFcIzonzc6DCOTceWAxcE9rmPuAbzrkTgf+i6aOBU4APjmjLA04FvkwwnOotBU5r248jcuw0zYVI+N50zpUD5Wa2D/hXqP0T4HgzSwMmA38PTbsBkNjE5+QBpUe0PROaiG2VmfVu1F4C5HfUDyDSGoWCSPiqGq0HGj0PEPy35AP2OufGtfI5h4CMFj678ZTpSaHtRTqFTh+JfKEc6NHWNzvn9gMbzewKAAsa28Smq4HjwvzYocTobJ3iDYWCSIhzbjewyMxWmNmv2/gxM4EbzOwjYCVN3z51ATDeGp1jasEZwAttrEXkmGmWVBEPmNldwL+cc6+3st0CYLpzbk/nVCaxTkcKIt74BZDS0gZmlgP8ToEgnUlHCiIi0kBHCiIi0kChICIiDRQKIiLSQKEgIiINFAoiItLg/wOej5yDKOi4zQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t=np.linspace(-3,18)\n",
    "Temp_ambient = 65\n",
    "Temp_t1 = 85\n",
    "plt.plot(t,T_analytical(Temp_ambient,Temp_t1,K,t),'-',label='Analytical')\n",
    "plt.legend()\n",
    "plt.xlabel('time (h)')\n",
    "plt.ylabel('Temperature (F)')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Now that we have a working numerical model, we can look at the results if the\n",
    "ambient temperature is not constant i.e. T_a=f(t). We can use the weather to improve our estimate for time of death. Consider the following Temperature for the day in question. \n",
    "\n",
    "    |time| Temp ($^o$F)|\n",
    "    |---|---|\n",
    "    |8am|55|\n",
    "    |9am|58|\n",
    "    |10am|60|\n",
    "    |11am|65|\n",
    "    |noon|66|\n",
    "    |1pm|67|\n",
    "\n",
    "    a. Create a function that returns the current temperature based upon the time (0 hours=11am, 65$^{o}$F) \n",
    "    *Plot the function $T_a$ vs time. Does it look correct? Is there a better way to get $T_a(t)$?\n",
    "\n",
    "    b. Modify the Euler approximation solution to account for changes in temperature at each hour. \n",
    "    Compare the new nonlinear Euler approximation to the linear analytical model. \n",
    "    At what time was the corpse 98.6$^{o}$F? i.e. what was the time of death? \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Answer 4a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def current_temp(time):\n",
    "    '''Returns the temperature at the current time\n",
    "    '''\n",
    "    Ambient_Temp = [55, 58, 60, 65, 66, 67]\n",
    "    length = len(time)\n",
    "    current_temp = np.zeros(length)\n",
    "    current_time = np.zeros(length)\n",
    "    for i in range(0,length):\n",
    "        current_time[i] = time[i] + 11\n",
    "        if time[i] < -2:\n",
    "            current_temp[i] = 55\n",
    "        elif time[i] < -1:\n",
    "            current_temp[i] = 58\n",
    "        elif time[i] < 0:\n",
    "            current_temp[i] = 60\n",
    "        elif time[i] < 1:\n",
    "            current_temp[i] = 65\n",
    "        elif time[i] < 2:\n",
    "            current_temp[i] = 66\n",
    "        else:\n",
    "            current_temp[i] = 67\n",
    "    return current_temp, current_time\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[55. 55. 55. 58. 60. 65. 66. 67. 67. 67. 67.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Temperature (F)')"
      ]
     },
     "execution_count": 258,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAeLklEQVR4nO3df5xVdb3v8dcbBmT0aGP8SAEV8eSUAg04kIW/RabUEsrUtA7pMX8d7V57hGlyj9gP80SPU5qhx/BX92GZIaC3TECO4tUSHAQDQfShmQ5oDnrHSEd+fu4few0O4/zYM8yazd7r/Xw8eMxa373XWp8Fw3uv/V1rfZciAjMzy45ehS7AzMx6loPfzCxjHPxmZhnj4DczyxgHv5lZxpQVuoB8DBgwIIYNG1boMszMisqyZcs2RMTAlu1FEfzDhg2jtra20GWYmRUVSX9trd1dPWZmGePgNzPLGAe/mVnGFEUff2u2bNlCXV0d7733XqFLsSLRr18/hg4dSp8+fQpdillBFW3w19XVsffeezNs2DAkFboc281FBG+++SZ1dXUcfPDBhS7HrKCKNvjfe+89h77lTRL9+/envr6+0KUUrXnL1zFj/lrWNzQyuKKcqTWVTBo9pKS3Xar7XLTBDzj0rVP8+9J185av46o5K2ncsg2AdQ2NXDVnJUDqQViobZfyPvvkrpl1aMb8tTtCqEnjlm3MmL+2ZLddyvvs4N8Fr7/+OmeddRaHHHIIhx12GCeffDLPP/98j9bQ0NDAzJkzP9D+5ptvUlVVRVVVFfvttx9DhgzZMb958+YerbEzbr/9dl5//fVCl2EtrG9o7FR7KWy7lPc5M8E/b/k6xl//3xx85e8Zf/1/M2/5ul1aX0QwefJkjjvuOF588UVWr17Nddddx9/+9re817Ft286f6Fu3bu10HW0Ff//+/VmxYgUrVqzgoosu4vLLL98x37dv305vpzu13O/muhL8Xfl7s84ZXFHeqfZS2HYp73Mmgr+pv2xdQyPB+/1luxL+jzzyCH369OGiiy7a0VZVVcXRRx/No48+yqmnnrqj/dJLL+XOO+8EcsNPfPe73+Woo47it7/9Lccddxzf+c53OPbYY7nhhhuor6/ni1/8ImPHjmXs2LE88cQTAEyfPp3zzjuP4447juHDh3PjjTcCcOWVV/Liiy9SVVXF1KlT867/rrvuYty4cVRVVXHJJZewfft2tm7dSkVFBVOnTmXMmDHU1NSwZMkSjj32WIYPH86DDz4IwKxZs5g8eTI1NTVUVlby/e9/P6/1Tps2jXHjxrF06VKuueYaxo4dy4gRI7jooouICH7zm9+wYsUKzjzzzB3fTIYOHUpDQwMATz75JBMmTABg2rRpXHjhhZx00kmce+65bN26lW9+85uMGzeOUaNGMWvWrC78q1pbptZUUt6n905t5X16M7WmsmS3Xcr7XNQnd/PVXn9ZV0+UrFq1iiOOOKJLy/br14/HH38cgFtuuYWGhgYWL14MwNlnn83ll1/OUUcdxSuvvEJNTQ1r1qwB4LnnnuORRx5h48aNVFZWcvHFF3P99dezatUqVqxY0ana586dyx//+EfKysq44IILuOeeezjjjDN4++23mThxIjNmzOBzn/sc06dPZ9GiRTzzzDNceOGFnHzyyQAsXbqUVatW0bdvX8aOHcupp55KWVlZu+sdM2bMjg+JyspKrr32WiKCs88+m4ceeogzzzyTn/3sZ9x0001UVVV1uB/Lly/nscceo1+/fsycOZNBgwaxdOlSNm3axJFHHsnEiRM58MADO/VvY61r+n9SiCtcCrXtUt7nTAR/IfvqWnPmmWe2Of/www+zevXqHfN///vf2bhxIwCnnHIKe+yxB3vssQeDBg3qVLdScw8//DBPPfUU1dXVADQ2NnLAAQcAUF5ezkknnQTAyJEj+dCHPkRZWRkjR47k5Zdf3rGOmpoa9t13XwAmTZrE448/ztatW9tcb9++fZk8efKO5RctWsSMGTN477332LBhA0cccQSf/exnO7Ufp512Gv369QNgwYIFrFmzhnvuuQeAt99+mxdeeMHB340mjR7SY5cy7i7bLtV9zkTwD64oZ10rIb8r/WWHH344s2fPbvW1srIytm/fvmO+5d3Fe+21V5vz27dv509/+hPl5R+sbY899tgx3bt37y73bUcE5513Ht/73vd2at+6detO/f+9evXasc1evXrttL2Wl0ZKane95eXlO5Z59913ufTSS3n66acZMmQI06ZNa/MO7OZ/l+39PUYEM2fO5MQTT8zr78AsyzLRx59Gf9kJJ5zApk2b+MUvfrGj7amnnmLx4sUcdNBBrF69mk2bNvH222+zaNGivNc7ceJEbrrpph3zHXXh7L333ju+EeRrwoQJ3HvvvWzYsAHIXQH0yiuvdGodCxYsoKGhgXfffZf777+f8ePH573exsZGevXqxYABA9i4cSP33Xdfm/szbNgwli1bBrDT+1qqqalh5syZOz6c1q5dS2NjYb7Rme3uMhH8k0YP4YdfGMmQinIEDKko54dfGLlLX6MkMXfuXBYuXMghhxzC4YcfzvTp0xk8eDAHHHAAZ5xxBqNGjeKcc85h9OjRea/3xhtvpLa2llGjRnHYYYdxyy23tPv+/v37M378eEaMGJH3yd2RI0dyzTXXMGHCBEaNGsXEiRM73W101FFHcfbZZzN69Gi+/OUvU1VVlfd6+/fvz5QpUxgxYgSTJ0/mk5/85I7Xzj33XM4///wdJ3enT5/OJZdcwtFHH93u1UgXXnghH/3oR6mqqmLEiBFcfPHFvtrHrA2KiELX0KHq6upo+SCWNWvW8PGPf7xAFWXbrFmzWLVqFT/96U8LXUqn+ffGskTSsoiobtmeiSN+MzN7XyZO7lr3Ov/88wtdgpntgqI+4i+Gbirbffj3xSynaIO/X79+vPnmm/7PbHlpGo+/6bp/sywr2q6eoUOHUldX5/HVLW9NT+Ayy7pUg19SBTALGAEEcF5E/EnSZcClwFbg9xFxRWfX3adPHz9JyTKpkA8HsdKQ9hH/DcBDEXG6pL7AnpKOB04DRkXEJkmDUq7BrGQU8uEgVjpS6+OXtA9wDHAbQERsjogG4GLg+ojYlLS/kVYNZqWmkA8HsdKR5snd4UA9cIek5ZJmSdoLOBQ4WtISSYsljW1tYUkXSKqVVOt+fLOc3W3AQStOaQZ/GTAGuDkiRgPvAFcm7fsCRwJTgXvVysNQI+LWiKiOiOqBAwemWKZZ8Sjkw0GsdKQZ/HVAXUQsSeZnk/sgqAPmRM5SYDswIMU6zEpGIR8OYqUjteCPiNeBVyU1/UaeCKwG5gEnAEg6FOgLbEirDrNSksaAg5Y9aV/Vcxlwd3JFz0vAueS6fG6XtArYDEwJ34VllrdCPhzESkOqwR8RK4APjAwHfCXN7ZqZWduKdsgGMzPrGge/mVnGOPjNzDLGwW9mljEOfjOzjHHwm5lljIPfzCxjivZBLGaF5DHxrZg5+M06yWPiW7FzV49ZJ3lMfCt2Dn6zTvKY+FbsHPxmneQx8a3YOfjNOslj4lux88lds05qOoHrq3qsWDn4zbrAY+JbMXNXj5lZxjj4zcwyxsFvZpYxDn4zs4xx8JuZZYyD38wsYxz8ZmYZ4+A3M8sYB7+ZWcakGvySKiTNlvScpDWSPtXstW9JCkkD0qzBzMx2lvaQDTcAD0XE6ZL6AnsCSDoAOAl4JeXtm5lZC6kd8UvaBzgGuA0gIjZHREPy8k+AK4BIa/tmZta6NLt6hgP1wB2SlkuaJWkvSZ8H1kXEM+0tLOkCSbWSauvr61Ms08wsW9IM/jJgDHBzRIwG3gGmA1cD/97RwhFxa0RUR0T1wIEDUyzTzCxb0gz+OqAuIpYk87PJfRAcDDwj6WVgKPC0pP1SrMPMzJpJLfgj4nXgVUlNjyU6EXg6IgZFxLCIGEbuw2FM8l4zM+sBaV/Vcxlwd3JFz0vAuSlvz8zMOpBq8EfECqC6ndeHpbl9MzP7IN+5a2aWMQ5+M7OMcfCbmWWMg9/MLGMc/GZmGdPuVT3JZZgnA0cDg4FGYBXwYEQ8l355ZmbW3do84pc0DVgCHA88A9wFPEDuw+Inkh6SNKJHqjQzs27T3hH/yoj4fhuv/UjS/sABKdRkZmYpai/4H2hvwYh4DXite8sxM7O0tXdyd1nThKSf9kAtZmbWA9oLfjWbPibtQszMrGe0F/x+OpaZWQlqr4//Y5KeJnfkX5lMk8xHRIxJvTozM+t27QX/yB6rwszMekybwR8RL/ZkIWZm1jPau4HrEUkXSxrcor1M0jGSbpPkB6uYmRWZ9rp6TgHOB+ZKGgK8BfRL/iwCfh4RtemXaGZm3am9rp53gRuBGyXtAQwCGiNiQ08VZ2Zm3S+vRy9GxCbg1ZRrMTOzHuBhmc3MMsbBb2aWMXkFv6Shko5PpveQtFe6ZZmZWVo6DH5J55EbqXNW0nQQcH+aRZmZWXryOeL/BnAk8HeAiHie3BU+HZJUIWm2pOckrZH0KUkzkvk/S5orqaLr5ZuZWWflE/zvRcTmphlJvdl55M723AA8FBEfAz4BrAEWAiMiYhTwPHBV50o2M7NdkU/wPyHpCqBf0s//G+B3HS0kaR9ywznfBhARmyOiISIWRMTW5G1PAkO7VrqZmXVFPsF/BbAReA74H+Tu2r06j+WGA/XAHZKWS5rVyknh84A/dKJeMzPbRe0Gf9Ktc3tE3BwRkyNiUjK9PY91lwFjgJsjYjTwDnBls3VfDWwF7m5j2xdIqpVUW19fn+/+mJlZB9oN/ojYBuwvqU8X1l0H1EXEkmR+NrkPAiRNAU4FzomIVh/4EhG3RkR1RFQPHDiwC5s3M7PW5DNkw0vA/5V0P7mjdgAi4sb2FoqI1yW9KqkyItYCJwKrJX0G+DZwbDIekJmZ9aB8gr+e3JU4eyZ/OuMy4G5Jfcl9gJwLPAXsASyUBPBkRFzUyfWamVkXdRj8EfG/urryiFgBVLdo/ueurs/MzHZdh8EvaSGtPHg9IiamUpGZmaUqn66eac2m+wFfBDalU46ZmaUtn66eJS2aFktanFI9ZmaWsny6evZpNtsLOALYP7WKzMwsVfl09TxLro9f5G64+gvw9TSLMjOz9OQT/MMjYkvzBkl5PbLRzMx2P/mM1dOyjx9gaXcXYmZmPaPNI3dJg8j15ZdLGsn7QzHvQ+dv5DIzs91Ee102p5AbPXMoMLNZ+0agyzd1mZlZYbUZ/BFxB7khlc+IiHt7sCYzM0tRPtfx3yupBjic3A1cTe3XpVmYmZmlI5/r+GcCFeSepnUHuTt3n0y5LjMzS0k+V/UcFRFnA28mA7Z9Ej8u0cysaOX1sPWmn5L2S+aHpVaRmZmlKp8bsR6UVAH8GFgBbAPuSrUqMzNLTbvBL6kX8IeIaAB+K+l3QHlEvNUj1ZmZWbfr6Jm724Ebms03OvTNzIpbPn38CyWdlnolZmbWI/Lp478U+JCkTUAjuaEbIiI+nGplZmaWinyCf0DqVZiZWY/psKsnIrYBXwK+nUzvD1SlXZiZmaWjw+CXdBNwPPDVpOld4JY0izIzs/Tk09Xz6YgYI2k5QES8JalvynWZmVlK8rmqZ0tyPX8ASOoPbE+1KjMzS00+wf9z4D5goKRrgceB/8hn5ZIqJM2W9JykNZI+JenDkhZKeiH5ue8u1G9mZp2Uz7DMv5S0DJiQNH0pIlbluf4bgIci4vSke2hP4DvAooi4XtKVwJXAt7tQuxnzlq9jxvy1rG9oZHBFOVNrKpk0ekihyzLbreVzxA/QG9gCbM53GUn7kBvK+TaAiNicDP1wGu+P9XMXMKkzBZs1mbd8HVfNWcm6hkYCWNfQyFVzVjJv+bpCl2a2W8vnqp6rgV8Dg8kNx/wrSVflse7hQD25p3gtlzRL0l7ARyLiNYDk56AuV2+ZNmP+Whq3bNuprXHLNmbMX1ugisyKQz5H718BxkbEtIi4GhgH/Esey5UBY4CbI2I08A65bp28SLpAUq2k2vr6+nwXswxZ39DYqXYzy8kn+P/KzucCyoCX8liuDqiLiCXJ/GxyHwR/k7Q/QPLzjdYWjohbI6I6IqoHDhyYx+YsawZXlHeq3cxy8gn+d4Fnk66aXwArgQZJ/ynpP9taKCJeB16VVJk0nQisBh4ApiRtU4D7u1y9ZdrUmkrK+/Teqa28T2+m1lS2sYSZQX43cP0++dOkM8/bvQy4O7mi5yXgXHIfNvdK+lfgFXLDQZh1WtPVO76qx6xzFBGFrqFD1dXVUVtbW+gyzMyKiqRlEVHdsj2fq3o+I+kpSW9IekvS/5Pkh7GYmRWpfLp6bgLOINe376EazMyKXD7BXwesSB7DaGZmRS6f4L8C+D+SHgU2NTVGxI1pFWVmZunJJ/ivJTdcQwXu6jEzK3r5BP+giDgi9UrMzKxH5HMD1yJJJ6ReiZmZ9Yh8gv/rwMOS/uHLOc3Mil8+XT0DUq/CzMx6TIdH/BGxjdywCt9OpvcHqtIuzMzM0pHPnbs3AccDX02a3gVuSbMoMzNLTz5dPZ+OiDGSlgNExFvJoGtmZlaE8jm5u0VSLyAAJPXH1/ObmRWtNoNfUtO3gZ8D9wEDJV0LPA78Rw/UZmZmKWivq2cpMCYifilpGTABEPCliFjVI9WZmVm3ay/41TQREc8Cz6ZfjpmZpa294B8o6ZttvRgRbT520czMdl/tBX9v4J9oduRvZmbFr73gfy0ivttjlZiZWY9o73JOH+mbmZWg9oL/xB6rwszMekybwR8RHoHTzKwE5XPnrpmZlRAHv5lZxuQzSFuXSXoZ2AhsA7ZGRLWkKnKje/YDtgKXRMTSNOuwdM1bvo4Z89eyvqGRwRXlTK2pZNLoIYUuy8zakGrwJ46PiA3N5n8EXBsRf5B0cjJ/XA/UYSmYt3wdV81ZSeOWbQCsa2jkqjkrARz+ZrupQnT1BLBPMv0hYH0BarBuMmP+2h2h36RxyzZmzF9boIrMrCNpH/EHsEBSAP8VEbcC/xOYL+nH5D54Pt3agpIuAC4AOPDAA1Mu07pqfUNjp9rNrPDSPuIfHxFjgM8C/ybpGOBi4PKIOAC4HLittQUj4taIqI6I6oEDB6ZcpnXV4IryTrWbWeGlGvwRsT75+QYwFxgHTAHmJG/5bdJmRWpqTSXlfXrv1FbepzdTayoLVJGZdSS14Je0l6S9m6aBicAqcn36xyZvOwF4Ia0aLH2TRg/hh18YyZCKcgQMqSjnh18Y6RO7ZruxNPv4PwLMldS0nV9FxEOS/gHckDzh6z2SfnwrXpNGD3HQmxWR1II/Il4CPtFK++PAEWlt18zM2uc7d83MMsbBb2aWMQ5+M7OMcfCbmWWMg9/MLGMc/GZmGePgNzPLGAe/mVnGOPjNzDLGwW9mljEOfjOzjHHwm5lljIPfzCxjHPxmZhnj4DczyxgHv5lZxjj4zcwyxsFvZpYxDn4zs4xx8JuZZYyD38wsYxz8ZmYZ4+A3M8sYB7+ZWcakGvySXpa0UtIKSbXN2i+TtFbSs5J+lGYNZma2s7Ie2MbxEbGhaUbS8cBpwKiI2CRpUA/UYGZmiUJ09VwMXB8RmwAi4o0C1GBmlllpB38ACyQtk3RB0nYocLSkJZIWSxrb2oKSLpBUK6m2vr4+5TLNzLIj7a6e8RGxPunOWSjpuWSb+wJHAmOBeyUNj4hovmBE3ArcClBdXR2YmVm3SPWIPyLWJz/fAOYC44A6YE7kLAW2AwPSrMPMzN6XWvBL2kvS3k3TwERgFTAPOCFpPxToC2xoaz1mZta90uzq+QgwV1LTdn4VEQ9J6gvcLmkVsBmY0rKbx8zM0pNa8EfES8AnWmnfDHwlre2amVn7fOeumVnGOPjNzDLGwW9mljEOfjOzjHHwm5lljIPfzCxjHPxmZhnj4Dczy5ieGI+/IOYtX8eM+WtZ39DI4IpyptZUMmn0kJLediH32cyKR0kG/7zl67hqzkoat2wDYF1DI1fNWQmQehAWatuF3GczKy4l2dUzY/7aHQHYpHHLNmbMX1uy2y7kPptZcSnJ4F/f0Nip9lLYdiH32cyKS0kG/+CK8k61l8K2C7nPZlZcSjL4p9ZUUt6n905t5X16M7WmsmS3Xch9NrPiUpInd5tOZhbiCpdCbbuQ+2xmxUXF8AyU6urqqK2tLXQZZmZFRdKyiKhu2V6SXT1mZtY2B7+ZWcY4+M3MMsbBb2aWMQ5+M7OMKYqreiTVA3/t4uIDgA3dWE4x8D5ng/c5G3Zlnw+KiIEtG4si+HeFpNrWLmcqZd7nbPA+Z0Ma++yuHjOzjHHwm5llTBaC/9ZCF1AA3uds8D5nQ7fvc8n38ZuZ2c6ycMRvZmbNOPjNzDKmpINfUoWk2ZKek7RG0qcKXVOaJF0u6VlJqyT9WlK/QteUBkm3S3pD0qpmbR+WtFDSC8nPfQtZY3drY59nJL/bf5Y0V1JFIWvsTq3tb7PXviUpJA0oRG1paWufJV0maW3yf/tH3bGtkg5+4AbgoYj4GPAJYE2B60mNpCHAN4DqiBgB9AbOKmxVqbkT+EyLtiuBRRHxUWBRMl9K7uSD+7wQGBERo4Dngat6uqgU3ckH9xdJBwAnAa/0dEE94E5a7LOk44HTgFERcTjw4+7YUMkGv6R9gGOA2wAiYnNENBS2qtSVAeWSyoA9gfUFricVEfEY8FaL5tOAu5Lpu4BJPVpUylrb54hYEBFbk9kngaE9XlhK2vg3BvgJcAVQcleltLHPFwPXR8Sm5D1vdMe2Sjb4geFAPXCHpOWSZknaq9BFpSUi1pE7GngFeA14OyIWFLaqHvWRiHgNIPk5qMD19LTzgD8Uuog0Sfo8sC4inil0LT3oUOBoSUskLZY0tjtWWsrBXwaMAW6OiNHAO5Te1/8dkj7t04CDgcHAXpK+UtiqrCdIuhrYCtxd6FrSImlP4Grg3wtdSw8rA/YFjgSmAvdK0q6utJSDvw6oi4glyfxsch8EpWoC8JeIqI+ILcAc4NMFrqkn/U3S/gDJz275Sry7kzQFOBU4J0r7ppxDyB3UPCPpZXLdWk9L2q+gVaWvDpgTOUuB7eQGbdslJRv8EfE68KqkyqTpRGB1AUtK2yvAkZL2TI4ITqSET2a34gFgSjI9Bbi/gLX0CEmfAb4NfD4i3i10PWmKiJURMSgihkXEMHKBOCb5f17K5gEnAEg6FOhLN4xOWrLBn7gMuFvSn4Eq4LoC15Oa5JvNbOBpYCW5f9uSvL1d0q+BPwGVkuok/StwPXCSpBfIXfVxfSFr7G5t7PNNwN7AQkkrJN1S0CK7URv7W9La2OfbgeHJJZ73AFO645udh2wwM8uYUj/iNzOzFhz8ZmYZ4+A3M8sYB7+ZWcY4+M3MMsbBb0UhGY3xfzebL5NUL+l3yfznJV2ZTE+X9K1k+muSBqdQz9XJJZQrJG1rNv0NSRdJ+pfu3may3dmShifT/2jjPZdKOjeN7VtpKCt0AWZ5egcYIak8IhrJXau/runFiHiA3E1cLX0NWEUnBqyTVNZs8LNWRcQPgB8k7/9HRFTlu/6uknQ40DsiXurgrbcDTwB3pF2TFScf8Vsx+QNwSjL9ZeDXTS8kR/Y3NX+zpNOBanI38a2QVC7p3yU9lTyz4NamcU8kPSrpOkmLgasl/UVSn+S1fSS93DTfkRbfOB6V9BNJjyn3TIixkuYkzw34frNlviJpaVLnf0nq3cqqz6HFHcmSfiDpGUlPSvoIQHIX78uSxuVTr2WPg9+KyT3AWco9YGYUsKS9N0fEbKCW3Dg2Vck3hZsiYmzyzIJycuPcNKmIiGMj4lrgUd7/kDkLuC8ZA6krNkfEMcAt5IL734ARwNck9Zf0ceBMYHzyzWEbuZBvaTywrNn8XsCTEfEJ4DHg681eqwWO7mK9VuIc/FY0IuLPwDByR/sPdnE1xydD3K4kNwbK4c1e+02z6VlAUz/5uexat0lTF9RK4NmIeC0ZX/0l4ABy4yodATwlaUUyP7yV9exPbqjxJpuB3yXTy8j93TR5g9worWYf4D5+KzYPkHvuwHFA/84smHxTmEnuKWWvSpoONH885TtNExHxhKRhko4l16/+gUcAdsKm5Of2ZtNN82WAgLsioqMnaDW2qHdLs3FbtrHz/+d+yfvNPsBH/FZsbge+GxEr83z/RnIDmcH7oblB0j8Bp3ew7C/JnUdI+yTpIuB0SYNgx/ODD2rlfWuAf85znYeSO6lt9gEOfisqEVEXETd0YpE7gVuSLpRNwC/IdbnMA57qYNm7yT0E49cdvG+XRMRqYBqwIBlJdiG5bp2Wfk/um04+xgMPd0uBVnI8OqdZG5Krgk6LiK8WuhYASeXAI+ROAm9r532jgW/uLnXb7sfBb9YKST8DPgucHBHPF7qeJpJqgDUR8Uo77zkJeCEiXu6xwqyoOPjNzDLGffxmZhnj4DczyxgHv5lZxjj4zcwyxsFvZpYx/x81j53YlKs5WgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "some_time = np.linspace(-5,5,11)\n",
    "current_Temp,current_Time = current_temp(some_time)\n",
    "print(current_Temp)\n",
    "plt.plot(current_Time,current_Temp,'o',label='Current Temperature')\n",
    "plt.legend()\n",
    "plt.xlabel('Miltary Time (h)')\n",
    "plt.ylabel('Temperature (F)')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[85.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
      "  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]\n",
      "[85.         82.25       79.878125   77.96988281 76.32402393 75.04197064\n",
      " 73.93619967 72.98247222 72.15988229 71.45039847 70.83846868 70.31067924\n",
      " 69.85546084 69.46283498 69.12419517 68.83211833 68.58020206 68.36292428\n",
      " 68.17552219 68.01388789 67.8744783  67.75423754 67.65052988 67.56108202\n",
      " 67.48393324 67.41739242 67.36000096 67.31050083 67.26780697 67.23098351\n",
      " 67.19922328 67.17183008 67.14820344 67.12782547 67.11024947 67.09509016]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "time_step = 0.5\n",
    "t = np.arange(0,18,time_step)\n",
    "Temp_ambient = 65\n",
    "current_Temp,current_Time = current_temp(t)\n",
    "Temp_t1 = 85\n",
    "K = 0.275\n",
    "T_numerical = np.zeros(len(t))\n",
    "T_numerical[0] = 85\n",
    "print(T_numerical)\n",
    "for i in range (1,len(t)):\n",
    "    T_numerical[i] = T_numerical[i-1] + (-K*(T_numerical[i-1] - current_Temp[i-1]))*time_step\n",
    "print(T_numerical)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Answer 4b\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Temperature (F)')"
      ]
     },
     "execution_count": 257,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU1fn48c+TTHaSkEBAICCoCCKELSwqrqgFFVxQigsILtRWrNp+a2lrlVrrr4srarUqVtsqlqLgUlEBpQICEkBZlX0JRAhLCCF78vz+uDcYkplkskwmy/N+veY1c8+9d+4zlzDP3HPuOUdUFWOMMaaikGAHYIwxpnGyBGGMMcYrSxDGGGO8sgRhjDHGK0sQxhhjvPIEO4D61LZtW+3atWuwwzDGmCZj1apVB1U1ydu6ZpUgunbtSlpaWrDDMMaYJkNEdvlaZ1VMxhhjvLIEYYwxxitLEMYYY7xqVm0Qxpimr6ioiPT0dPLz84MdSrMSGRlJcnIyYWFhfu9jCcIY06ikp6cTGxtL165dEZFgh9MsqCqHDh0iPT2dbt26+b1fQKuYROR+EdkgIutFZKaIRIrINBHZKyJfuY8rfOw7QkS+FZGtIjI1YEGunQVP9YZprZ3ntbMCdihjTPXy8/Np06aNJYd6JCK0adOmxldlAbuCEJFOwE+BXqqaJyKzgHHu6qdU9fEq9g0FngcuA9KBlSLynqpurNcg186C938KRXnO8tE9zjJAyth6PZQxxn+WHOpfbc5poBupPUCUiHiAaGCfn/sNBraq6nZVLQTeAq6u9+gWPvJ9cihTlOeUG2NMCxewBKGqe4HHgd1ABnBUVT9xV08RkbUi8qqIJHjZvROwp9xyultWiYhMFpE0EUnLzMysWZBH02tWboxpMebMmYOI8M0339T6PSZOnMjs2bOr3Oaxxx47afncc8+t1bGmTZvG44/7rJiplYAlCPeL/2qgG9ARiBGRW4AXgNOBfjiJ4wlvu3sp8zqzkaq+pKqpqpqalOS1t7hv8ck1KzfGtBgzZ85k2LBhvPXWWwE9TsUE8cUXXwT0eDURyCqmS4EdqpqpqkXAO8C5qrpfVUtUtRR4Gac6qaJ0oHO55WT8r57y3/CHICzq5LKwKKfcGNNi5eTksHTpUmbMmHEiQSxatIiLLrqI66+/np49e3LzzTdTNiPnI488wqBBg+jduzeTJ0+m4kydCxcu5Nprrz2xPH/+fK677jqmTp1KXl4e/fr14+abbwagVatWJ7b785//TJ8+fejbty9Tpzr36rz88ssMGjSIvn37MmbMGHJzcwN2HgJ5m+tuYKiIRAN5wHAgTUQ6qGqGu821wHov+64EuotIN2AvTuP2TfUeodsQrQumQfZeCiSKyFHTrYHamEbid+9vYOO+7Hp9z14d43h41NlVbjN37lxGjBjBmWeeSWJiIqtXrwZgzZo1bNiwgY4dO3LeeeexdOlShg0bxpQpU3joIeeH5fjx4/nggw8YNWrUife75JJLuPvuu8nMzCQpKYm///3vTJo0iVGjRvHcc8/x1VdfVYph3rx5zJ07lxUrVhAdHc3hw4cBuO6667jzzjsBePDBB5kxYwb33HNPvZybigLZBrECmA2sBta5x3oJ+LOIrBORtcDFwP0AItJRRD509y0GpgAfA5uAWaq6ISCBpoxFfraR9bHnc7Q0ioKzrgvIYYwxTcfMmTMZN8656XLcuHHMnDkTgMGDB5OcnExISAj9+vVj586dAHz22WcMGTKEPn368Omnn7Jhw8lfVyLC+PHj+de//kVWVhbLli1j5MiRVcawYMECJk2aRHR0NACJiYkArF+/nvPPP58+ffrwxhtvVDpWfQpoRzlVfRh4uELxeB/b7gOuKLf8IfBh4KI7WWivK2m/YjFpK/9H6rnDG+qwxpgqVPdLPxAOHTrEp59+yvr16xERSkpKEBGuuOIKIiIiTmwXGhpKcXEx+fn5/OQnPyEtLY3OnTszbdo0r/0Nyq4YIiMjueGGG/B4qv76VVWvt6ZOnDiRuXPn0rdvX1577TUWLVpU58/si43F5Dr9vDGUqJC1Zm6wQzHGBNHs2bOZMGECu3btYufOnezZs4du3bqxZMkSr9uXJYO2bduSk5Pj866ljh070rFjRx599FEmTpx4ojwsLIyioqJK219++eW8+uqrJ9oYyqqYjh07RocOHSgqKuKNN96oy0etliUIV0RcO3ZEp9AlcxHFJaXBDscYEyQzZ848qUEZYMyYMbz55ptet2/dujV33nknffr04ZprrmHQoEE+3/vmm2+mc+fO9OrV60TZ5MmTSUlJOdFIXWbEiBGMHj2a1NRU+vXrd+IW1t///vcMGTKEyy67jJ49e9b2Y/pFKra2N2WpqalalwmDvnnn/9Fz7R9Ju2YRqf3612Nkxhh/bdq0ibPOOivYYQTElClT6N+/P7fffntQju/t3IrIKlVN9ba9XUGU0/W8GwDY/+U7QY7EGNPcDBw4kLVr13LLLbcEOxS/2Wiu5US2P4N94V1pn7GQktJHCA2x8WCMMfVj1apVwQ6hxuwKooKcrj+gX+kmvtq8I9ihGGNMUFmCqCD5nDF4pJTdy62ayRjTslmCqCD61EFkhbah9e75lbrLG2NMS2IJoqKQEA4nD2dwyRrW7twf7GiMMSZoLEF40X7QdcRIAd8u+yDYoRhjgkBE+PnPf35i+fHHH2fatGkNGkNaWho//elPa7XvRRddRF1u+S9jCcKLmJ6XkCdRRG3/yKqZjGnsAjBtcEREBO+88w4HDx6shwBrrri4mNTUVKZPnx6U45exBOGNJ4LM9sMYUrSSTfuOBjsaY4wvZdMGH90D6PfTBtcxSXg8HiZPnsxTTz1VaV3FSYDKhudetGgRF154IWPHjuXMM89k6tSpvPHGGwwePJg+ffqwbds2ADIzMxkzZgyDBg1i0KBBLF26FHAm/Jk8eTKXX345EyZMYNGiRVx11VWAM/z4pEmT6NOnDykpKbz99tsA/PjHPyY1NZWzzz6bhx+uOOxd3Vk/CB8SBlxD7IfzWbB8Ab3GXB/scIxpmeZNhe/W+V6fvhJKCk4uK8qDd6fAqte973NKHxj5x2oPfffdd5OSksIDDzzgd7hff/01mzZtIjExkdNOO4077riDL7/8kmeeeYZnn32Wp59+mnvvvZf777+fYcOGsXv3bn7wgx+wadMmwOkrsWTJEqKiok4ahO/3v/898fHxrFvnnIsjR44A8Ic//IHExERKSkoYPnw4a9euJSUlxe94q2MJwofY3ldQ8mEIfPMhqmNsEnVjGqOKyaG68hqIi4tjwoQJTJ8+naioqOp3AAYNGkSHDh0AOP3007n88ssB6NOnD5999hngDOO9cePGE/tkZ2dz7NgxAEaPHu31WAsWLDhpZruEBGem5lmzZvHSSy9RXFxMRkYGGzdutATRIKITOdhmEIMyl7M2/Sh9O7cOdkTGtDzV/dJ/qrdbvVRBfGeY9N86H/6+++5jwIABTJo06USZx+OhtNQZ0FNVKSwsPLGu/HDgISEhJ5ZDQkIoLi4GoLS0lGXLlnlNBDExMV7j8Db0944dO3j88cdZuXIlCQkJTJw40esw43VhbRBViOt3Nd1D9vK/ZcuDHYoxxpsATxucmJjI2LFjmTFjxomyrl27nhg249133/U6VHdVLr/8cp577rkTy95mk6tunyNHjpCdnU1MTAzx8fHs37+fefPm1SgOfwQ0QYjI/SKyQUTWi8hMEYkUkb+IyDcislZE5oiI15/mIrLTnXnuKxGp+/1atRDVx5kysPSbD2wIcGMao5SxMGq6c8WAOM/1PG3wz3/+85PuZrrzzjv53//+x+DBg1mxYoXPX/2+TJ8+nbS0NFJSUujVqxcvvvhitfs8+OCDHDlyhN69e9O3b18+++wz+vbtS//+/Tn77LO57bbbOO+882r82aoTsOG+RaQTsATopap5IjILZ4a4fcCnqlosIn8CUNVfetl/J5Cqqn7fZ1bX4b69yX5qCN8cUXJv/oCLerSr1/c2xlTWnIf7DrbGNty3B4gSEQ8QDexT1U/cOacBlgPJAY6hTmJSRjMwZDPzVwZu3ldjjGmMApYgVHUv8DiwG8gAjqrqJxU2uw3wVXGmwCciskpEJvs6johMFpE0EUnLzMysj9BPEuoJJxTl0S3XUPrk2fXSCccYY5qCgCUIEUkArga6AR2BGBG5pdz63wDFgK9JVc9T1QHASOBuEbnA20aq+pKqpqpqalJSUr1+BtbOgiVPuPFCSHZ6vXTCMcZUzUYwqH+1OaeBrGK6FNihqpmqWgS8A5wLICK3AlcBN6uPqFV1n/t8AJgDDA5grN4tfMTpdFNeUZ5TbowJiMjISA4dOmRJoh6pKocOHSIyMrJG+wWyH8RuYKiIRAN5wHAgTURGAL8ELlTVXG87ikgMEKKqx9zXlwMN/618NL1m5caYOktOTiY9PZ1AVBm3ZJGRkSQn16zJN2AJQlVXiMhsYDVOVdIa4CVgAxABzHc7fixX1btEpCPwiqpeAbQH5rjrPcCbqvpRoGL1KT7ZRyecRt2ubkyTFhYWRrdu3YIdhiHAPalV9WGg4ghSZ/jYdh9whft6O9A3kLH5ZfhDTptDuWqmAsKJqKdOOMYY05hZT+qqVOyEA3xQPIitp4wMblzGGNMALEFUJ2Us3L8epmVR1K4P3UMymLtmX7CjMsaYgLMEUQNh/caRErKdVau/pLTU7rAwxjRvliBqovcYFGHo8YWs2n0k2NEYY0xAWYKoibgOlHa9gOtClzJntd3qaoxp3ixB1FBov3F0lgPsWbuI/KKSYIdjjDEBYwmipnpeRUloBJcVf85H678LdjTGGBMwliBqKjKOkJ5XMtqznFkrtgc7GmOMCRhLELUgKT+kNceI3L2IbZk5wQ7HGGMCwhJEbZwxnNLIRK4LXcpbX+4OdjTGGBMQliBqIzSMkD7XcblnFfPSNlNQbI3VxpjmxxJEbaX8kHAtZEjBMj7esD/Y0RhjTL2zBFFbyYPQhK6Mi1xm1UzGmGbJEkRtiSB9xpJaupat27ay4+DxYEdkjDH1yhJEXaSMRVCu9izjrZV2FWGMaV4sQdRF2+7QsT/jo5czOy2dwuLSYEdkjDH1JqAJQkTuF5ENIrJeRGaKSKSIJIrIfBHZ4j4n+Nh3hIh8KyJbRWRqIOOsk5Qf0qVwKwm521mwyRqrjTHNR8AShIh0An4KpKpqbyAUGAdMBRaqandgobtccd9Q4HlgJNALuFFEegUq1jpxR3idG/EwI98+C57qDWtnBTsqY4yps0BXMXmAKBHxANHAPuBq4HV3/evANV72GwxsVdXtqloIvOXu1/hsX4SI0Io8BHXmsH7/p5YkjDFNXsAShKruBR4HdgMZwFFV/QRor6oZ7jYZQDsvu3cC9pRbTnfLKhGRySKSJiJpmZmZ9fkR/LPwEdAKbQ9FeU65McY0YYGsYkrA+dXfDegIxIjILf7u7qXM6xRuqvqSqqaqampSUlLtgq2Loz7mhfBVbowxTUQgq5guBXaoaqaqFgHvAOcC+0WkA4D7fMDLvulA53LLyTjVU41PfHLNyo0xpokIZILYDQwVkWgREWA4sAl4D7jV3eZW4F0v+64EuotINxEJx2ncfi+Asdbe8IcgLOqkonwinHJjjGnCAtkGsQKYDawG1rnHegn4I3CZiGwBLnOXEZGOIvKhu28xMAX4GCepzFLVDYGKtU5SxsKo6RDvXPAowm8LJ7C69WVBDswYY+pGVL1W7TdJqampmpaWFrwA9qyEGZfyGLex78zxPHfTgODFYowxfhCRVaqa6m2d9aSuT50HQfIg7or4mI/X72NvVl6wIzLGmFqzBFHfzrmbxIK9XCJp/GPZzmBHY4wxtWYJor71HAXxXfi/2AXMXLGb4wXFwY7IGGNqxRJEfQv1wNC76J6/jq4F3/LOausPYYxpmixBBEL/8Wh4LD+PW8CrS3dSWtp8bgQwxrQcVSYIEQkXkWtE5Al3NNZXReRnItKzoQJskiLjkIG3cn7hEvIP7mbRZm99AY0xpnHzmSBE5EFgBXAx8DXOwHrv4QzA95SIfCQivRskyqZoyI8QlLujFzJjyY5gR2OMMTXmqWLdOlV91Me6P7vDZHT2sd607oL0upobvpnPY1tHsykjm7M6xAU7KmOM8VtVVUxVDm2hqhmq+mU9x9O8nDOFiJIcbgr/nL8vtasIY0zTUlWCWFX2QkSeboBYmp/kVOg8hB9HfsJ7X6VzMKcg2BEZY4zfqkoQ5YfcviDQgTRb59xNm8J9XFi6kteW7gx2NMYY47eqEoTdm1kfel4FUW14Lvw5fvbFEEqfPNtmmzPGNAlVNVL3FJHVOFcSPdzXuMuqqjYSnT/Wvw0F2YRR5Jy57HRnSlJwRoI1xphGqqoE0afBomjOFj4CpUUnl5VNSWoJwhjTiPlMEKq6rSEDabZsSlJjTBNVVUe5z0TkxyLSsUK5R0QuEJEZIjIp8CE2cT6mHi2N69TAgRhjTM1UVcV0JXAHMEdEOgGHgUj3sRB4XlV9zs4jIj2Af5crOg14CDgH6OGWtQayVLWfl/13AseAEqDY14QWjd7wh5w2h6Lv54bI1zCWd/kJFwUvKmOMqVZVVUy5wHRguohEAO2APFU96M8bq+q3QD8AEQkF9gJzVPVEnwoReQI4WsXbXOzv8RqtsnaGhY+cqFba6zmVqZt78r/iEiI8oUEMzhhjfPNrNFdVLVDVPXX4sh4ObFPVXWUFIiLAWGBmLd+z6UgZC/evh2lZcMlvOL1kK+2ObWD2KmuHMMY0Xg013Pc4KieC84H9qrrFxz4KfCIiq0Rksq83FpHJIpImImmZmZn1FG4ADbkLjW7D71rN4YVF2ygqKQ12RMYY41XAE4SIhAOjgf9UWHUjVV89nOf2tRgJ3C0iXntzq+pLqpqqqqlJSUn1EnNARcQiw+6nf9FqOmatZu6avcGOyBhjvPIrQYhIsohc7L6OEJGYGhxjJLBaVfeXez8PcB0nN2KfRFX3uc8HgDnA4Bocs3FLvR1tdQoPxbzDXz/bSolNKGSMaYSqTRAichvOyK6vuEWnAu/W4BjerhQuBb5RVa+V8CISIyKxZa+By4H1NThm4xYejVzwf/Qu3kCnIyv4YO2+YEdkjDGV+HMF8VNgKJANoKqbce5oqpaIRAOXAe9UWFWpTUJEOorIh+5ie2CJiHwNfAn8V1U/8ueYTcaACWh8Mg9Gvc3T8zdbW4QxptGpqh9EmXxVLXRuOjpxy6pUvYvDvVW2jZfyiV7K9gFXuK+3A339OUaT5YlALvwlPd+7h9OOLObfK0/jlqGnBjsqY4w5wZ8riKUi8gAQ6bZD/Bv4ILBhtRB9b0ITT+O30XN4Zv63HC8oDnZExhhzgj8J4gGcHs3fAPfi9KL+TSCDajFCPchFv6Zr8XYG5S3mVZu72hjTiFSZINzqpFdV9QVVvVZVr3FfW4V5fel9HSSdxW9j5vLy51s5fLww2BEZYwxQTYJQ1RKgg4iENVA8LU9IKFz8azoU7WYRd5Dwl3bwVG+bVMgYE3T+NFJvBxaLyLvA8bJCVZ0esKhamuJ8QEiUHGf56B6bVMgYE3T+tEFkAvOBaCCp3MPUl4WPUGmG17JJhYwxJkiqvYJQ1d82RCAtmk0qZIxphKpNECIyn0o/b0FVLw9IRC1RfLJTreSt3BhjgsSfNogHy72OBMYABYEJp4XyMqlQoYays/fPODOIYRljWjZ/qphWVCj6n4j8L0DxtEwVJhVSTwQUF/H4xlheHK6EhPjVcd0YY+qVP1VMceUWQ4CBQIeARdRSpYw9kSjk6F7k2cFMOvgks1YOYtyQrsGNzRjTIvlzF9MGnJFUNwBrcHpR3xnIoFq8+E54Rj7GOaEb2f7Rs2TlWuc5Y0zD8ydBnKaqXVS1s6p2U9VLgKWBDqylkwETyEm+gHtL/8mM9xcFOxxjTAvkT4Ko2AYBzhDcJpBEaHX9X/GEhnLOhmmsT88KdkTGmBbGZ4IQkXYi0heIEpE+IpLiPobhdJozgda6MyWX/Z5zQzaw5N+PU2ozzxljGlBVjdRXArcBycBfy5UfA6zzXAOJHno73636DxMzX6Dgz/8mKj/T6R8x/CEbhsMYE1A+E4Sq/h34u4iMVdUajxwnIj04ec7p04CHgNY4jdyZbvmvVfXDCrsjIiOAZ4BQ4BVV/WNNY2gWRGg3cDTy8XIk/4BTZmM1GWMagKhWX20hIj8AzsbpKAeAqj7m90GcYcP3AkOASUCOqj5ezfabcaYrTQdWAjeq6saqjpOamqppaWn+htV0PNXbR0/rznB/85mq2xjT8ERklaqmeltXbSO1iPwVuBX4GRAF3AKcUcMYhgPbVHWXn9sPBraq6nZVLQTeAq6u4TGbDxuryRgTBP7cxTRMVW8CDrkD9w3BaZeoiXHAzHLLU0RkrYi8KiIJXrbvBJT/yZzullUiIpNFJE1E0jIzM71t0vT5GpPJxmoyxgSQPwkiv+xZRE5xl7v6ewARCQdGA/9xi14ATgf6ARnAE95281LmtS5MVV9S1VRVTU1KaqajkA9/CMKiTioqJcQpN8aYAPEnQXwoIq2Bx4GvgJ3A7BocYySwWlX3A6jqflUtcactfRmnOqmidKBzueVkYF8Njtm8pIyFUdMhvjOKcFxiCKGUo0ea6RWTMaZRqHIsJhEJAeapahbwHxH5AIhS1cM1OMaNlKteEpEOqprhLl6LM4xHRSuB7iLSDadxexxwUw2O2fy4YzUJcPDgMVY+O5phnz2Enj4ESfbavmSMMXVS3ZzUpTi3mpYt59UkOYhINM6dSO+UK/6ziKwTkbXAxcD97rYdReRD9zjFwBTgY2ATMEtVN/h73Obu1LaxpF/0FBmlCeS+cQvk1iRfG2OMf6q9zVVEfg+kqeq7DRNS7TXb21y9KClVfvXc6zx6+Odo1wuImPA2hPhTY2iMMd+r022uOL/k54hInogcFpEjImI/WYMsNET40Y3X82jJRCJ2fop+/pdgh2SMaWb8mVGubcCjMLVyelIrOg3/Me8s/IZrF/0/KMqF9W87/SNsOA5jTB1VewWhqiXADcAv3dcdcG5RNY3A7eefxlvt7mc/CejSp90e1/r9cBxrazxKijHGAP71pH4OpzF5vFuUC7wYyKCM/zyhITw6dgiiWrnzSFGeM42pMcbUgj9VTOeq6gARWQOgqofdzm+mkTizfSwqPuaLsOE4jDG15E8jdZHbH0IBRKQNUBrQqEzNxXsdicSG4zDG1Jo/CeJ54G0gSUR+BywB/hTQqEyNyfCHKfWcPBwHIWE2HIcxptaqrWJS1X+IyCrgUrfoBlW1MaYbm5SxhAC58x4mMncfJaERhJUWVBrDyRhj/OVvz6pQoAgorME+pqGljCX6l5v4dcpiUvJeJLttf5h9O+xaFuzIjDFNkD93Mf0GZyyljjiD5r0pIr8KdGCm9h4edTadktpwXdZPKY5LhpnjIPPbYIdljGli/LkauAUYpKoPqupvcEZfnRDYsExdRIWH8uyN/dmdH8XUqIfR0HD41xjIzqh+Z2OMcfmTIHZxcluFB9gemHBMfTmrQxwPXnkWs7d7mNvracg7Aq9cCk/2gmmtnWlMrROdMaYK/iSIXGCDiLwiIi8D64AsEXlSRJ4MbHimLsYPPZXLerXngS+E/T1vgex0yN6L9bQ2xvjDn45y/3UfZZYHKBZTz0SEP49J4crpiyld+3blDcp6Wtt4TcYYL/y5zXVGQwRiAiMhJpy/jU+l/csHvU/kaj2tjTE++HMX0wgRWSkiB2y476apT3I8+dGneF8Z56MHtjGmxfOniuk5YCxO24PfQ2yISA/g3+WKTgMeAjoBo3D6VGwDJrlTmlbcfydwDCgBin1NaGH8Ez3yEQrn3kN4af7JKyJiofA4hMcEJzBjTKPlTyN1OvCVqhapaknZo7qdVPVbVe2nqv2AgTiN3XOA+UBvVU0BNgNV9am42H0PSw51lTKW0NHTyQxtR6kKhTGdYOAkOPgtvHEDFBwLdoTGmEbGnyuIB4D3RWQRUFBWqKrTa3Cc4cA2Vd2Fc9tsmeXA9TV4H1MHof1+iKf7tVz4/BIKi0p5/8JhtOt2Prx9J/zzWrh5NkS1DnaYxphGwp8riN/hVPO0BpLKPWpiHE5v7IpuA+b52EeBT0RklYhM9vXGIjJZRNJEJC0zM7OGYbU8CTHhvDwhley8Yu761yoKel4DY1+HfV/Bi+dbPwljzAmiqlVv4ExoPbDWB3DmjtgHnK2q+8uV/wZIBa5TL0GISEdV3Sci7XCqpe5R1c+rOlZqaqqmpaXVNtQW5cN1GfzkjdVcPzCZv1yfgnzyW1j27MkbhUXBqOl2G6wxzZj7He+1Gt+fK4iFInJJHY4/ElhdITncClwF3OwtOQCo6j73+QBO28XgOsRgKriiTwfuHd6d2avSeWr+Ztg4t/JGNiOdMS2aPwniTmCBiOTU8jbXGylXvSQiI4BfAqNVNdfbDiISIyKxZa+BywEbYrye3Xdpd36Y2pnpn25FffWHsH4SxrRY/jRSt63tm4tINHAZ8KNyxc8BEcB8EQFYrqp3iUhH4BVVvQJoD8xx13uAN1X1o9rGYbwTEf5wbW8OHMtn7842JMvByhtFJTR8YMaYRqHaKwj3ltYbgF+6rzsA/fx5c1XNVdU2qnq0XNkZqtq57BZYVb3LLd/nJgdUdbuq9nUfZ6vqH2rz4Uz1PKEhPH/zAP4dN4k8rTDVuIRA3mGnmqnUZpk1pqXxpyf1c8DFwHi3KBd4MZBBmYYVHe5h4o9+weMRd7OPtigC8Z3h6udhwK2w+AmYNd7pUGeMaTH8qWI6V1UHiMgaAFU97N6ZZJqRNq0iuPVHDzD6hfOIDAvlndvPpV1cJPS9EdqdBR//Gl79AfS7BZY957RNxCc7c17bXU7GNEv+NFIXiUgITr8ERKQNNRhywzQdXdpE8/eJgzl8vJAJr35JVm4hiMDQH8NNsyBzC3z0S2eocBsy3Jhmz2eCEJGyq4vngbeBJBH5HbAE+FMDxGaCoE9yPH8bP5DtB49zy4wVHM0tclZ0v8x7L2u7FdaYZod6iuYAABqdSURBVKuqK4gvAVT1H8CDwOPAEeAGVX2rAWIzQXJ+9yReGj+Qzd/lMP7VFRzNc5NEzn7vO9itsMY0S1UliBOzB6jqBlV9RlWfVlXrj9ACXNSjHS+OH8CmjGwmvPol2flFTpuDNzE1HXnFGNMU+BxqQ0TSAZ9Tiqpqo5tu1IbaqH/zN+7nJ2+sok+neN4cuofIefc51UonuL8jLvkNDPsZhIQGJU5jTO3UdqiNUKAVEOvjYVqAy3q159kbB7A2/Si3rOhC/sinnVtgy26FHfU09B4Dnz4Kr4+C5S86A/3ZgH/GNHlVXUGsVtUBDRxPndgVRODMW5fBlJlrGNglgRkTU4mNDPt+pSp8/Ra8fy+UFJy8ow34Z0yjVtsrCG8zGJsWamSfDjwzrh+rdx/hxpeXczCnXCIQgX43QnRi5R3tLidjmqyqEsTwBovCNAlXpXTk5VtT2XoghxteXMaewxXGWjz2nfcd7S4nY5oknwlCVWsyYqtpIS7u0Y437hjK4eOFjHnhC779rtxUpb7ucvJEwJGdDRKfMab++NOT2piTDDw1gf/cdQ4icMOLX5C20/0tMfwhp82hvJAwp43i+aGw+EkoKXIarq0h25hGzxKEqZUz28cy+65zadMqgltmrODTb/Y7DdGjpp98l9M1f4WfroEzhsPC38HTKfDuFBuuw5gmoNopR5sSu4up4R3MKWDi379kU8Yxfn91b24a0sX3xt/Og7duBi2pvC6+M9xvfTCNaWh1nXK0tgftISJflXtki8h9IpIoIvNFZIv77HVGGhEZISLfishWEZkaqDhN3bRtFcHMO4cy7Iy2/HrOOqa9t4HiEh9jOfYYCepjnTVkG9PoBCxBqOq3ZZMCAQNx5pGYA0wFFqpqd2Chu3wSEQnFGSRwJNALuFFEegUqVlM3sZFhvDpxEHcM68ZrX+xk0msrvx/kryJfDdlRrZ32CWNMo9FQbRDDgW2qugu4GnjdLX8duMbL9oOBre7McoXAW+5+ppEKDREevKoXfx6TwvLth7j2r0vZnplTeUNvDdkSAnlH4PkhsGGONWIb00g0VIIYB8x0X7dX1QwA97mdl+07AXvKLae7ZaaRGzuoM2/eOZSsvCKueX4pi7dknryBt4bsa/8G42ZCaDj8ZyLMmWyN2MY0AgFvpHZnn9sHnK2q+0UkS1Vbl1t/RFUTKuxzA/ADVb3DXR4PDFbVe7y8/2RgMkCXLl0G7tq1K4Cfxvhrz+Fc7vxHGlsO5PDLET248/zTEKmmc35pCfzlDGce7IqsEduYgAhKI3U5I4HVqlo2mcB+EengBtYBOOBln3Sgc7nlZJwkU4mqvqSqqaqampRkw043Fp0To5n943O57Kz2PPbhN9zxehpHjhdWvVNIqFPV5M3RPU4CMcY0mIZIEDfyffUSwHvAre7rW4F3veyzEuguIt3cK5Bx7n6mCWkV4eGFWwYwbVQvPt+SyZXTF7Nql48EUMZXIzbAswNh5QxnfCdrpzAm4AJaxSQi0ThtCaep6lG3rA0wC+gC7MaZoe6wiHQEXlHVK9ztrgCexhl2/FVV/UN1x7N+EI3X2vQs7n5zNRlZ+fziB06VU0iIlyqntbOcNofyc06ERcGAW2HPl7BvNYTHQnE+lBadvI2NGmtMjVVVxWQd5UyDOZpXxNS31zJv/Xdc0rMdT9zQl4SY8Mobrp3ljAB7NN25ohj+kPPFrwo7l8Ab1zsJoiJrpzCmxixBmEZDVfnn8l08+sEmWkeH8acxKVzc09uNbFWY1hrw8Xf7q70Q0cp3kjHGnCTYjdTGnCAiTDinK+/85FwSosOZ9NpKHpj9tTPntb+qaqd4oge8Ngres/GejKkrSxAmKHp3iue9e87jJxedzuxV6Yx46nOWbDno387eOtuFRcGFU6HX1bBzMRRXmNnOJi4ypsYsQZigifCE8sCInrz943OJDA/llhkr+O3c9RwvKK56R2+d7UZNh4t/5Ywe68vRPZB/9PtluxPKmCpZG4RpFPKLSvjLx9/y6tIdJCdE8cjo3jVvmyjzVG+3esmL0HDofjnEdYDV/4LiCndL2Z1QpoWxNgjT6EWGhfLbq3rx1p1DCQ8NYdJrK7nrn6vYl5VX/c4V+aqCuuhXMOgOSE+DL18+OTmAVUMZU4FdQZhGp7C4lJcXb+fZT7cQIsJ9l3Zn0nndCAutwe+Zqu5iKi2BR9rg806oKWnQ5gxY9x+7E8o0e3abq2mS9hzO5Xfvb2DBpgP0aB/Lo9f2ZlDXxPp586qqoQBikpxhP0rLtYdYFZRphqyKyTRJnROjeeXWQbw0fiA5BcXc8OIy7n5zNbsP5db9zX1VQ434I1z5BBQcOzk5gFMF9clvoKRcuTV0m2bMriBMk5BbWMyL/9vOy59vp7i0lPFDu3LPJWd474ntr6qqoarqjBcRB6eeCxGxsOm9k2+ptasM08RYFZNpNvZn5/PU/M3MSttDqwgPUy45gwnndCUyLLR+D+SrCiq6DZw1CnYshsPbvO8b1wl+ttF5bT26TSNnCcI0O99+d4z/N28Ti77NJDkhinuHd+fa/p3w1KQhuyq+Bg0sf3VQ1VVGu7MhKhHSl588lWrF97AEYoLM2iBMs9PjlFhemzSYf90+hNbRYfxi9lqGP/k/Zq9Kp7iktO4H8NUZr/yXt68hPyLiIPYU2LWk8jzbRXkw7wHI+Bq+etNJQjYkiGmk7ArCNHmqyoJNB3h6wWY27Mvm1DbR3HNJd67p17H+rii8qe4qo6orjKpUHJXWrjJMAFkVk2kRvCWKH194Otf071T/bRRlqvry9tWOEXsK/OAxmH2b7/ftPASSekBRPmx8F0qqaAi3BGLqwBKEaVHKEsUzCzezfm82bVuFM+Gcrtwy9FQS63LXU01Vd4XhK4GEx0CH/pD5DeT6GMAwIhYu/wNk7YJlz588P4a3O6ksiRgfgpYgRKQ18ArQG+da+zbgPqCHu0lrIEtV+3nZdydwDCgBin19gPIsQZjyVJVl2w7x0uLtLPo2k8iwEK4fmMztw06jW9uYhgmiqi/mujaEVyUyHq56ClqfCvvWwPzfVn0cSyAtVjATxOvAYlV9xZ1bOlpVs8qtfwI4qqqVBsBxE0Sqqvo5BrQlCOPb5v3HeGXxduau2UdRaSnDe7bnlqFduKB7kvepTxtKdV/Mvq4y4pNh0jx4OoVaJRCAyNYw+lk4sBGWPFX3qxBLMk1SUBKEiMQBX+PMR13pICIiOHNSX6KqW7ys34klCFPPDhzL5x9f7OKtlbs5mFNI58Qobhp8KjekJtO2VUSww6usttVUcZ3g5tmQtRtm/rB2xw6LhsF3Qkw7OLwd1vzLd1uIP1dDZZ/HkkyjEqwE0Q94CdgI9AVWAfeq6nF3/QXAkz4DE9kBHMH5efQ3VX3Jx3aTgckAXbp0Gbhr1676/iimGSosLuXjDd/xr+W7WLHjMGGhwsjeHbhpSBeGdEvE+f3SSNS1mspnEukIN74Ff7sQn1choeFQUug7thAPdOwP363zPk94TBLc/B+ISoBti+Djqb5jbagkY0noJMFKEKnAcuA8VV0hIs8A2ar6W3f9C8BWVX3Cx/4dVXWfiLQD5gP3qOrnVR3TriBMbWw9cIx/Ld/N26vTOZZfTOfEKK7rn8x1AzpxapsGaquoC3++EGtzFRLfGe5b50yy9Keu+Ewip10M2z+rffyeSGcmwG/+C4U5ldfHJMFNs5yG+W2fwfyHfM/jUd1nbcgrnYZ4j3oQrARxCrBcVbu6y+cDU1X1ShHxAHuBgaqa7sd7TQNyVPXxqrazBGHqIq+whI82ZPD2qr0s3XYQVRjUNYHrBiRzZUoH4iLDgh1i7QXqKqSsz4av9THtYNTTkJcF7/7Ed3ytuzjVYbUV4oEOfeG79SdXg5WJiIPzfwZLnob8rMrrY9rBLW87n3vrAlgwzXebjD/nqz4Slb/JrI6C2Ui9GLhDVb91v+RjVPUXIjIC+JWqXuhjvxggRFWPua/nA4+o6kdVHc8ShKkv+7LymPvVXt5elc62zOOEe0K46MwkrkzpwPCz2tMqwhPsEOtXXa9CAppkkpzG9IIceOcO35/h9OGwbWHNPneNCEQnOslOSyqvDo2Abhc41XLbPq08IRU4V0FD7oIVf4OC7MrroxJh5J8gJBQ+/AXkHqq8TcWOlHUUzATRD+c213BgOzBJVY+IyGs4Vxcvltu2I/CKql4hIqcBc9xVHuBNVf1DdcezBGHqm6qyNv0oc9bsZd76DPZnFzT/ZOFLfdT9ByvJxCc7E0E9OxCy91ZeH93WudIpLoC3b/d9DlJvh7QZvtd3HOAMr7J/ne9tEGp951nZ/tO8XAXV9t2so5wxdVdaqqzafYT/rs04KVkMO6Mtw89qx/Ce7TklPjLYYTZuwU4yAU1C5X7ZV/seZzufsaLYDjDxv06S+cdoyNlf9XHqgSUIY+pZWbL4cF0GCzbtZ89h5wund6c4hvdsz6Vntad3p7jGdTdUcxHou5gaov3A2iAaniUIEwyqypYDOSzYtJ8FG/ezZk8WqpAUG8H5Z7Tl/DPbct4ZbWkXa1cXTYbdxeSsswRhTP06mFPAZ98c4PMtB1myJZMjuc6Q32d1iOOC7k6yGHhqAjEtpe3CNGqWIIwJktJSZcO+bBZvzWTx5oOs2nWEwpJSQkOEPp3iGXJaIkO6JZLaNbFp30ZrmixLEMY0ErmFxazadYQV2w+zfPshvk7PoqhECRHo1TGOAV0SGHhqAgO6JJCcEGVtGCbgLEEY00jlFZawZvcRlu84zModh/k6PYvcQuce+7atIhjQpTUDTk0gJTme3p3i7SrD1LuqEoRVghoTRFHhoZx7RlvOPaMtAMUlpXy7/xird2exZtcRVu8+wicbv7/VsVvbGHp3iqdPpzj6dGpNr45xxEdZ0jCBYVcQxjRyh48Xsn7vUdbtPcq6dOd5b9b3tz52ah1Fz1NiOatDHGd1iKNnh1i6tokhNJjDmJsmw64gjGnCEmPCueDMJC44M+lE2aGcAtbtPcqmjGNsysjmm++yWbQ5k5JS5wdfuCeE05NacWb7VnRv14oz2sXSvX0rTk2MDuw83aZZsQRhTBPUplUEF/Vox0U92p0oyy8qYeuBHDZlZLPlQA5b9h9j1a4jvPvVvhPbeEKELm2iOa1tDN3axtCtbSv3OYZ2sRHBnTzJNDqWIIxpJiLDQundyWnMLu94QTHbM4+zef8xtmbmsPPgcXYcPM7iLQcpKC49sV2EJ4QuidF0SYymc2I0p7ZxXndKiKJT6yhirYG8xbEEYUwzFxPhoU9yPH2ST04cpaVKRnY+OzKPs+NgDrsP57L7cC67DuWybPuhE3dTlYmPCqNT66gTCaNj60g6xEfRIT6SU+IjaR8XSZhVXzUrliCMaaFCQsT5wm8dxbDubU9ap6ocOl7I7sO57D2Sx96sPNKPOK93HTrOF1sPcrxCAhGBpFYRtI+LpH1cBEmxkbSLjaBdXATtYyNJio2gbWwEbVuFE+EJbciPamrJEoQxphIRoW2rCLcvRkKl9arKsYJiMrLyyTiax3dH88k46rzen13A3qx81uzO4tBx79OVxkV63GQRQVKrCBJiwkiMiaBNTDiJMeG0iQknwX0dHxVGZJgllGCwBGGMqTERIS4yjLhTwuhxSqzP7QqLSzmYU8CBYwUcyM7nYE4hB3MKvn8cK2TTd9kcPl5IljtmlTdRYaEkRIcRHx1OQnQYraPDiI8KIy7KeS57xEU6ZbGRHmIjPcRFWnKpC0sQxpiACfeE0LF1FB1bR1W7bXFJKVl5RRw+XsihnEInaeQ5iePI8UKy8orIyi3kSG4Rm/fncDSviKN5RRSWa2j3GkNoCLGRHlpFemgV4Txi3dcx7nJ0uIeYiFBiTpSFEhXmITo8lJiIUKLCPUSHhRIVHkqEJ6TFDIES0AQhIq1xZpTrjTOF0m3AD4A7gUx3s1+r6ode9h0BPAOE4sw098dAxmqMCS5PaMiJai3a+79fflEJR/OKyMotIju/iGP5RRzLLyY7v5jsPOf1sfwijhcUk1NQzLH8YjKO5pNTUExOfjHHC4vJL6o6yZQXIs4VTVS4h6jwEOd1WCiRJx4hJy1HhIUQ6fH+HOEJJdwTQoQn5MRzhCeE8FCnPNwTQlioOK9DGz4xBfoK4hngI1W9XkTCgWicBPGUqj7uaycRCQWeBy4D0oGVIvKeqm4McLzGmCam7Iu4fVzt59soLiklt6iE4wXFHC9wnwuLySssIbewhNzCYve5hLzCEvKK3EeF5azcQvKKSsgvKiXfLSsoKqWwxP8EVJXwUCdhhHlCCAt1kka4J4SkVhHMuuucejlGeQFLECISB1wATARQ1UKg0M8MOBjYqqrb3fd6C7gasARhjKl3ntAQ4kJDAjYYYmmpUljiJI2C4lIKikrJLy6hsLiUgmK3zC0vKC6hqEQpLC6lsOx1ibO+qKSUIve5sKSUwmKlqKSUmIjAtLME8griNJxqpL+LSF9gFXCvu26KiEwA0oCfq+qRCvt2AspP6JoODPF2EBGZDEwG6NKlS/1Fb4wx9SQkRIgMCW1yDeaB7NXiAQYAL6hqf+A4MBV4ATgd6AdkAE942dfbZYbXUQVV9SVVTVXV1KSkJG+bGGOMqYVAJoh0IF1VV7jLs4EBqrpfVUtUtRR4Gac6ydu+ncstJwP7vGxnjDEmQAKWIFT1O2CPiPRwi4YDG0WkQ7nNrgXWe9l9JdBdRLq5jdvjgPcCFasxxpjKAn0X0z3AG+6X/HZgEjBdRPrhVBntBH4EICIdcW5nvUJVi0VkCvAxzm2ur6rqhgDHaowxphybMMgYY1qwqiYMsqEXjTHGeGUJwhhjjFeWIIwxxnjVrNogRCQT2FXL3dsCB+sxnECxOOtfU4nV4qxfTSVOCGysp6qq105kzSpB1IWIpPlqqGlMLM7611RitTjrV1OJE4IXq1UxGWOM8coShDHGGK8sQXzvpWAH4CeLs/41lVgtzvrVVOKEIMVqbRDGGGO8sisIY4wxXlmCMMYY41WLShAiMkJEvhWRrSIy1ct6EZHp7vq1IjIgSHF2FpHPRGSTiGwQkXu9bHORiBwVka/cx0NBinWniKxzY6g0EFZjOKci0qPcefpKRLJF5L4K2wTtfIrIqyJyQETWlytLFJH5IrLFfU7wsW+Vf9MNEOdfROQb9992jjsPvbd9q/w7aYA4p4nI3nL/vlf42LfBzmcVsf67XJw7ReQrH/sG/pyqaot44IwKuw1nprtw4GugV4VtrgDm4UxYNBRYEaRYO+DMnQEQC2z2EutFwAeN4LzuBNpWsb5RnNMKfwff4XQOahTnE2dq3gHA+nJlfwamuq+nAn/y8Vmq/JtugDgvBzzu6z95i9Ofv5MGiHMa8H9+/G002Pn0FWuF9U8ADwXrnLakK4gT81yrMz922TzX5V0N/EMdy4HWFeavaBCqmqGqq93Xx4BNONOwNkWN4pyWMxzYpqq17XFf71T1c+BwheKrgdfd168D13jZ1Z+/6YDGqaqfqGqxu7gcZ3KvoPJxPv3RoOcTqo5VRAQYC8wMZAxVaUkJwts81xW/dP3ZpkGJSFegP7DCy+pzRORrEZknImc3aGDfU+ATEVnlzg9eUWM7p+Pw/R+uMZzPMu1VNQOcHwxAOy/bNLZzexvO1aI31f2dNIQpblXYqz6q7Brb+Twf2K+qW3ysD/g5bUkJwp95rv2eC7shiEgr4G3gPlXNrrB6NU41SV/gWWBuQ8fnOk9VBwAjgbtF5IIK6xvNORVn4qrRwH+8rG4s57MmGtO5/Q1QDLzhY5Pq/k4C7QXgdKAfkIFTdVNRozmfrhup+uoh4Oe0JSUIf+a5bjRzYYtIGE5yeENV36m4XlWzVTXHff0hECYibRs4TFR1n/t8AJhD5TnGG805xfmPtFpV91dc0VjOZzn7y6ri3OcDXrZpFOdWRG4FrgJuVrdyvCI//k4CSlX3q2qJqpYCL/s4fqM4nwAi4gGuA/7ta5uGOKctKUH4M8/1e8AE986bocDRssv8huTWPc4ANqnqkz62OcXdDhEZjPNveajhogQRiRGR2LLXOA2WFecYbxTn1OXzF1ljOJ8VvAfc6r6+FXjXyzZBn7tdREYAvwRGq2quj238+TsJqArtXtf6OH7Qz2c5lwLfqGq6t5UNdk4D2QLe2B44d9RsxrlT4Tdu2V3AXe5rAZ53168DUoMU5zCcS9u1wFfu44oKsU4BNuDcabEcODcIcZ7mHv9rN5bGfE6jcb7w48uVNYrziZO0MoAinF+xtwNtgIXAFvc50d22I/BhVX/TDRznVpx6+7K/0xcrxunr76SB4/yn+/e3FudLv0Owz6evWN3y18r+Nstt2+Dn1IbaMMYY41VLqmIyxhhTA5YgjDHGeGUJwhhjjFeWIIwxxnhlCcIYY4xXliCM8UJEWovIT8otdxSR2QE61jVlo8eKyGsicr2XbZJE5KNAHN8YXyxBGONda+BEglDVfapa6Yu7njwA/LWqDVQ1E8gQkfMCFIMxlViCMMa7PwKnu2Pt/0VEupaN2S8iE0Vkroi8LyI7RGSKiPxMRNaIyHIRSXS3O11EPnIHU1ssIj0rHkREzgQKVPVgueILROQLEdle4WpiLnBzAD+zMSexBGGMd1NxhgXvp6q/8LK+N3ATzvg3fwByVbU/sAyY4G7zEnCPqg4E/g/vVwnn4QwUWF4HnN70V+EkqjJpOCN8GtMgPMEOwJgm6jN15uo4JiJHgffd8nVAijsS77nAf9whngAivLxPByCzQtlcdQaV2ygi7cuVH8AZbsGYBmEJwpjaKSj3urTccinO/6sQIEtV+1XzPnlAfBXvXX4I6kh3e2MahFUxGePdMZzpXmtFnfk7dojIDXBibu6+XjbdBJzh59ueSQOPgmpaNksQxnihqoeApSKyXkT+Usu3uRm4XUTKRtz0Nn3l50B/KVcPVYWLgf/WMhZjasxGczUmyETkGeB9VV1QzXafA1er6pGGicy0dHYFYUzwPYYzX4VPIpIEPGnJwTQku4IwxhjjlV1BGGOM8coShDHGGK8sQRhjjPHKEoQxxhivLEEYY4zx6v8D5LL7CcvdzkkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "t=np.linspace(0,18)\n",
    "t_1=np.arange(0,18,time_step)\n",
    "Temp_ambient = 65\n",
    "Temp_t1 = 85\n",
    "plt.plot(t,T_analytical(Temp_ambient,Temp_t1,K,t),'-',label='Analytical')\n",
    "plt.plot(t_1,T_numerical,'o-',label='Numerical')\n",
    "plt.legend()\n",
    "plt.xlabel('time (h)')\n",
    "plt.ylabel('Temperature (F)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[106.60180545 101.57301578  97.27298091  93.49273047  90.2903125\n",
      "  87.475       85.          82.525       80.3903125   78.54914453\n",
      "  76.96113716  75.5914808   74.41015219  73.39125626  72.51245853\n",
      "  71.75449548  71.10075235  70.5368989   70.0505753   69.6311212\n",
      "  69.26934203  68.9573075   68.68817772  68.45605329  68.25584596\n",
      "  68.08316714  67.93423166  67.80577481  67.69498077  67.59942091\n",
      "  67.51700054  67.44591296  67.38459993  67.33171744  67.28610629\n",
      "  67.24676668]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "time_step = 0.5\n",
    "t = np.arange(0,18,time_step)\n",
    "Temp_ambient = 65\n",
    "Temp_t1 = 85\n",
    "K = 0.275\n",
    "current_Temp,current_Time = current_temp(t)\n",
    "T_numerical = np.zeros(len(t))\n",
    "begin = -3\n",
    "if begin >= 0:\n",
    "    T_numerical[0] = 85\n",
    "    for i in range (1,len(t)):\n",
    "        T_numerical[i] = T_numerical[i-1] + (-K*(T_numerical[i-1] - current_Temp[i-1]))*time_step\n",
    "else:\n",
    "    T_numerical[abs(begin)*2] = 85\n",
    "    for i in range (abs(begin*2)+1,len(t)):\n",
    "        T_numerical[i] = T_numerical[i-1] + (-K*(T_numerical[i-1] - current_Temp[i-1]))*time_step\n",
    "    for i in range(abs(begin*2),0,-1):\n",
    "        T_numerical[i-1] = T_numerical[i] + (K*(T_numerical[i] - current_Temp[i]))*time_step\n",
    "print(T_numerical)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Answer 4b\n",
    "Shown in the figure below is the temperature was 98.6 around -3 hours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 276,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Temperature (F)')"
      ]
     },
     "execution_count": 276,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU1fn48c8zk8kKhGxAICCgIAJhDaCAiqKAK4jCFzcEF6p1b/tV+qtV1Npq64pW/WqxalUsbqhVrIKAgrKERQQCIiCQECAEyEL2yfn9cSdDCDPJkGSWJM/79ZrXnbvM3Cc3yTxzzrnnHDHGoJRSSgHYgh2AUkqp0KFJQSmllJsmBaWUUm6aFJRSSrlpUlBKKeUWFuwAGiIxMdF07do12GEopVSTsmbNmoPGmCRP+5p0UujatSvp6enBDkMppZoUEdnlbZ9WHymllHLTpKCUUspNk4JSSim3Jt2moJRqHsrLy8nMzKSkpCTYoTQrkZGRpKSk4HA4fH6NJgWlVNBlZmbSunVrunbtiogEO5xmwRhDbm4umZmZdOvWzefXafWRUiroSkpKSEhI0ITQiESEhISEky59aVJQSoUETQiNrz7XtEUmhawjxfxlQQYH8rX+UimlqvNbUhCR10TkgIhsrLZtkohsEpFKEUmrcfzvReRnEdkqImP9FRfA0dIK/m/pDv67aZ8/T6OUamI++ugjRIQtW7bU+z2mTZvG+++/X+sxf/7zn49bHz58eL3ONWvWLJ588sl6vdYbf5YUXgfG1di2EZgIfFN9o4j0BqYAfVyveVFE7P4KrEe7VpyaFMOCjZoUlFLHzJ07l5EjR/Luu+/69Tw1k8J3333n1/OdDL8lBWPMN8ChGtsyjDFbPRw+HnjXGFNqjNkJ/AwM9VdsIsJFfZNZufMQh46W+es0SqkmpLCwkOXLlzNnzhx3UliyZAmjRo3iqquuolevXlx77bVUzVb5yCOPMGTIEPr27cuMGTOoOYvlokWLuOKKK9zrX331FRMnTmTmzJkUFxczYMAArr32WgBatWrlPu6vf/0rqamp9O/fn5kzZwLw6quvMmTIEPr378+VV15JUVGR365DqNyS2glYUW0907XtBCIyA5gB0KVLl3qfcFzfDryw+Ge+2ryP/xlS//dRSjWuhz/dxOa9+Y36nr07tuGhy/rUesz8+fMZN24cPXv2JD4+nrVr1wKwbt06Nm3aRMeOHRkxYgTLly9n5MiR3HHHHTz44IMAXH/99fznP//hsssuc7/f+eefz+23305OTg5JSUn885//ZPr06Vx22WW88MILrF+//oQYFixYwPz581m5ciXR0dEcOmR9r544cSK33HILAA888ABz5szhzjvvbJRrU1OoNDR7aiL3OHm0MeYVY0yaMSYtKcnjIH8+6dOxDZ3jo/j8R61CUkpZVUdTpkwBYMqUKcydOxeAoUOHkpKSgs1mY8CAAfzyyy8ALF68mGHDhpGamsrXX3/Npk2bjns/EeH666/nrbfe4siRI3z//fdcdNFFtcawcOFCpk+fTnR0NADx8fEAbNy4kbPPPpvU1FTefvvtE87VmEKlpJAJdK62ngLs9ecJq6qQ/rl8J3nF5cRG+d7jTynlP3V9o/eH3Nxcvv76azZu3IiI4HQ6EREuvvhiIiIi3MfZ7XYqKiooKSnh17/+Nenp6XTu3JlZs2Z57A9QVTKIjIxk0qRJhIXV/pFrjPF4G+m0adOYP38+/fv35/XXX2fJkiUN/pm9CZWSwifAFBGJEJFuQA9glb9POq5vB8qdhkUZ+/19KqVUCHv//feZOnUqu3bt4pdffmHPnj1069aNZcuWeTy+KgEkJiZSWFjo9W6jjh070rFjR/70pz8xbdo093aHw0F5efkJx48ZM4bXXnvN3WZQVX1UUFBAcnIy5eXlvP322w35Uevkz1tS5wLfA6eLSKaI3CQiV4hIJnAW8JmI/BfAGLMJmAdsBr4AbjfGOP0VW5UBKW1Jjo3Uu5CUauHmzp17XKMwwJVXXsk777zj8fi2bdtyyy23kJqayoQJExgyZIjX97722mvp3LkzvXv3dm+bMWMG/fr1czc0Vxk3bhyXX345aWlpDBgwwH276aOPPsqwYcO48MIL6dWrV31/TJ9IzRbzpiQtLc00dJKdWZ9sYu6q3az944XERIRKbZpSLUtGRgZnnHFGsMPwizvuuIOBAwdy0003BeX8nq6tiKwxxqR5Oj5Uqo+C5qK+HSitqGTx1gPBDkUp1cwMHjyYDRs2cN111wU7FJ+1+K/GaV3jSWwVzoKN+7i0X8dgh6OUakbWrFkT7BBOWosvKdhtwpg+HVi85QAl5X5vxlBKqZDW4pMCWFVIRWVOlv6UE+xQlFIqqDQpAGd2TyA2ysEXeheSUqqF06QAOOw2xvRuz8KM/ZRVVAY7HKWUCpqWmRQ2zINn+sKsttZywzwuSu1AQUkFy7cfDHZ0SqkgEBF++9vfuteffPJJZs2aFdAY0tPTueuuu+r12lGjRtHQW/ShJSaFDfPg07sgbw9grOWnd3F2yWJaR4TxhY6FpFTo8/DFrqEiIiL48MMPOXgwOF8MKyoqSEtLY/bs2UE5f5WWlxQWPQLlxcdvKy/GsfhPnH9GO77cvI8Kp1YhKRWyvHyxa2hiCAsLY8aMGTzzzDMn7Ks5cU7VUNdLlizh3HPPZfLkyfTs2ZOZM2fy9ttvM3ToUFJTU9m+fTsAOTk5XHnllQwZMoQhQ4awfPlywJokZ8aMGYwZM4apU6eyZMkSLr30UsAaynv69OmkpqbSr18/PvjgAwBuu+020tLS6NOnDw899FCDfmaP16HR3zHU5WV63X7RhR34eP1eVu48xIjTEgMbl1LKsmAm7PvR+/7M1eAsPX5beTF8fAesecPzazqkwkWP13nq22+/nX79+nHffff5HO4PP/xARkYG8fHxdO/enZtvvplVq1bx3HPP8fzzz/Pss89y9913c++99zJy5Eh2797N2LFjycjIAKy+DMuWLSMqKuq4ge4effRRYmNj+fFH61ocPnwYgMcee4z4+HicTiejR49mw4YN9OvXz+d469LykkJsiusbxonbz+3ZjuhwO//ZsFeTglKhqmZCqGv7SWjTpg1Tp05l9uzZREVF+fSaIUOGkJycDMCpp57KmDFjAEhNTWXx4sWANST25s2b3a/Jz8+noKAAgMsvv9zjuRYuXHjcDHBxcXEAzJs3j1deeYWKigqys7PZvHmzJoUGGf2gVdSsXoXkiILRDxIVbmdsnw58tiGbWZf3ISLMbzOCKqW8qesb/TN9vXyx6wzTP2vw6e+55x4GDRrE9OnT3dvCwsKorLSqlY0xlJUdm7Gx+tDaNpvNvW6z2aioqACgsrKS77//3uOHf0xMjMc4PA2jvXPnTp588klWr15NXFwc06ZN8zhkd0O0vDaFfpPhstnWHxCAPcJa7zcZgAkDO5FfUsHiLdqRTamQNPpB64tcda4vdo0hPj6eyZMnM2fOHPe2rl27uoes+Pjjjz0Oe12bMWPG8MILL7jXPc26VtdrDh8+TH5+PjExMcTGxrJ//34WLFhwUnH4ouUlBbASwL0bYfRDVpGz02D3rhGnJpDYKoL567KCGKBSyqvjvtiJtaz2xa4x/Pa3vz3uLqRbbrmFpUuXMnToUFauXOn12703s2fPJj09nX79+tG7d29efvnlOl/zwAMPcPjwYfr27Uv//v1ZvHgx/fv3Z+DAgfTp04cbb7yRESNGnPTPVpeWPXR2fjY80xtG3AMXHGvFf+TTzby1YherH7hAZ2RTKgCa89DZwaZDZ5+MNsnQYwysfwecFe7NEwZ2pMxZyYIfs4MYnFJKBV7LTgoAA6+Dwn3w80L3ptROsXRPiuEjrUJSSrUwmhR6joOYJFj3L/cmEeGKAZ1YufMQWUeKa3mxUqqxNOWq7FBVn2uqScHugP5T4KcvoPDY7GvjB3QC4OP1WlpQyt8iIyPJzc3VxNCIjDHk5uYSGRl5Uq9ref0UPBl4PXz3PPzwLoywBqPqkhBN2ilxfLQ2i9vOPfWE+4WVUo0nJSWFzMxMcnL0VvDGFBkZSUpKykm9RpMCQNLpkDIU1r0Fw+8EVwKYMLATD8zfyObsfPp0jA1ykEo1Xw6Hg27dugU7DIVWHx0z6Ho4uNUaV8XlktRkHHbRPgtKqRZDk0KVPleAIwbWvuneFBcTzqjT2/Hx+r04K7WuUynV/GlSqBLR2koMmz6C0kL35isGduJAQSnfb88NYnBKKRUYmhSqG3Q9lBXC5vnuTef3akfriDDts6CUahE0KVTXeRgk9IC1x/osRDrsXJyazBcbsykucwYxOKWU8j9NCtWJWD2c96yAg9vcmycM7MTRMidfZewPYnBKKeV/mhRq6n81iP24Hs7DusXTMTaSD9d6mbVNKaWaCU0KNbVuD+37Wp3ZXJOC2za+x5WDU1j6U44Oe6GUatY0KdS0YR7kZICppPqk4NNaW/0X5q32MOOTUko1E5oUalr0CDjLjt9WXkzCisc5p0cS89L3UOGsDE5sSinlZ5oUasrz0m6Ql8nVQzuTnVfC0p90fBalVPOkSaGmWC+DR8WmMPqM9iS2imDuqt2BjUkppQJEk0JNniYFD7MmBXfYbUxOS+HrLQfIztMGZ6VU8+O3pCAir4nIARHZWG1bvIh8JSLbXMu4avt+LyI/i8hWERnrr7jqVHNScIAeF7onBZ8ypAuVBt5L19tTlVLNjz9LCq8D42psmwksMsb0ABa51hGR3sAUoI/rNS+KiN2PsdWu32S4dyPMOgI9xsKu76DcKhl0SYhm5GmJ/Hv1Hh0kTynV7PgtKRhjvgEO1dg8HnjD9fwNYEK17e8aY0qNMTuBn4Gh/ortpAy/A4oOwoZ/uzddPbQLWUeK+WabNjgrpZqXQLcptDfGZAO4lu1c2zsB1TsAZLq2BV/XsyG5P3z3AlRat6Je2Ls9CTHhzF2pDc5KqeYlVBqaPc116bFuRkRmiEi6iKQHZOo+ETjrTsjdBtu+BCA8zMZVg1NYtOUAB/JL/B+DUkoFSKCTwn4RSQZwLQ+4tmcCnasdlwLs9fQGxphXjDFpxpi0pKQkvwbr1mcCtEmxhr5w+Z8hnXFWGt5bow3OSqnmI9BJ4RPgBtfzG4CPq22fIiIRItIN6AGsCnBs3tkdcOatsGsZ7F0HQPekVpzVPYG5q3ZTqQ3OSqlmwp+3pM4FvgdOF5FMEbkJeBy4UES2ARe61jHGbALmAZuBL4DbjTGhNXnBoBsgoo3VtuBy9bAuZB4uZtnPB4MYmFJKNZ4wf72xMeZqL7tGezn+MeAxf8XTYJFtYNBUWPESXDAL2nZmbJ/2xEU7eHf1bs7pGaCqLKWU8qNQaWhuGobdai1XvgxARJidqwan8N9N+7WHs1KqWdCkcDLadoa+E2HNG1CSB8DUs7pijOHN73cFOTillGo4TQon66w7oKzASgxA5/hoxvTuwDsrd1NUVhHk4JRSqmE0KZysjgOsDm0rXwZnOQA3nd2NvOJyPlybFeTglFKqYTQp1EfyAMjPgkcT4Zm+pOUvJLVTLK8t36m3pyqlmjRNCidrwzxI/8ex9bw9yKd38eApm9iRc5SlOh6SUqoJ06RwshY94h4x1a28mLTts2nXOoLXlu0MTlxKKdUINCmcLC/TdUpeFjcM78q32w7y0/6CAAellFKNQ5PCyaplus6rh3YhIsympQWlVJNVa1IQkXARmSAiT4nIXNdsar8RkV6BCjDkeJquE4Hz/0h8TDgTB6Xw4boscgtLgxKeUko1hNekICIPACuB84AfsCbF+QRraIxnROQLEekbkChDSc3pOqPiAQOOSABuHNGVsopK3tG5FpRSTZAY4/kWShEZb4z52ONO3ENfdzbGBG0007S0NJOenh6s01sqnfDimSB2uG052OxMfW0VGdn5LLv/PCLCgjerqFJKeSIia4wxaZ721VZ99Eltb2qMyQ5mQggZNjuMmgk5GbDpIwBuGtmNnIJSPtuQHeTglFLq5NSWFNZUPRGRZwMQS9PV+wpo1wcW/xmcFZzTI5HT2rVizrKdeCuJKaVUKKotKVSfIvMcfwfSpNlscN7/g0PbYcO/ERFuHNGNTXvzWf5zbrCjU0opn9WWFPQr7snodYk1/MXSJ8BZzsRBnWjfJoIXFm8LdmRKKeWz2pJCLxFZKyLrqj1fKyLrRGRtoAJsMkTgvD/AkV2w7i0iHXZmnHMqK3YcYvUvh4IdnVJK+aS2mddSAxZFc9HjQkgZAt/8DfpfzTVDu/Di4p+ZvWgb/7ppWLCjU0qpOnktKRhjttf2CGSQTUZVaSE/C9a+QVS4nZvP7s632w6yfs+RYEenlFJ1qq3z2mIRuU1EOtbYHiYi54jIHBGZ7v8Qm5juo+CUkfDtU1BWxPVnnUJslIMXvta2BaVU6KutTeESwAF8JCKZIrJBRH4CdgDTgZeMMf8MRJBNigic/wco3A9P96LVXxL5NvxOord+yKa9ecGOTimlauW1TcEYUwTMBmaLSATQDig2xhwMVHBNVl4miM09j3Ob0n084fgHcz9uS5/b7g9ycEop5Z1Po6QaY0qNMXs0Ifho0SNgKo/bFCVljNn3Ctt0WG2lVAjTobP9wcucCx3J5e+Lfw5wMEop5TtNCv7gZc6Fgoj2fPLDXn45eDTAASmllG98SgoikiIi57meR4hIjH/DauI8zblgj4ALHsRht/HiEi0tKKVCU51JQURuxBoxtWq2+lMAr0NqK06cc8EWBpFtiB14JVcP7cKHa7PYc6go2FEqpdQJfCkp3AWcCeQDGGN+wroTSdWm32S4dyPMOgLXzIOjOfDdbH51bndsNmH2Iu23oJQKPb4khRJjTFnViojYOX4EVVWX00ZD7/Hw7VMkV+7nhrNO4YO1mWzdp3ciKaVCiy9JYbmI3AdEutoV/g38x79hNUNj/2LNzrZgJr8edRox4WH87b9bgh2VUkodx5ekcB9QAGwB7gYWAX/wZ1DNUmwnGHU//LSAuMxF3DrqVBZmHNARVJVSIaXWpOCqKnrNGPOSMeYKY8wE1/PK2l6nvDjz15DUCxbcx41D29OudQSPL9iis7MppUJGrUnBGOMEkkXEEaB4mje7Ay5+Eo7sJmrlc9x7YU/W7DrMV5v3BzsypZQCfKs+2gF8KyK/F5G7qh7+DqzZ6nY2pE6G5c8xqWsp3ZNi+Ot/t1Lh1MKXUir4fEkKOcBXQDSQVO2h6mvMnwAbYa+MYFHBBF4/ciOrP/2/YEellFK1zrwGgDHmj419UhG5G7gF69bWV40xz4pIPNadTV2BX4DJxpjDjX3ukLBzKRgnOMsRIMV2kIT1D1F2ShzhA6cEOzqlVAvmS4/mr0Tky5qP+p5QRPpiJYShQH/gUhHpAcwEFhljemDd4TSzvucIeYsegcry4zZFUUrpF7OCE49SSrnUWVIAHqj2PBK4EihtwDnPAFa45mtARJYCVwDjgVGuY94AlgDNc/IBL6OoxpTuI6+onNhobddXSgVHnSUFY8zKao+lxpi7sL7l19dG4BwRSRCRaOBioDPQ3hiT7TpnNl6G0hCRGSKSLiLpOTk5DQgjiLyMorrXJPC8TtuplAoiX6qP2lR7tBWR0UByfU9ojMkAnsBqvP4C+AGoOInXv2KMSTPGpCUlNdH2bk+jqAKrOk3l9e9+0Yl4lFJB48vdR5uwvt1vAtZh9Wa+pSEnNcbMMcYMMsacAxwCtgH7RSQZwLU80JBzhLSao6i2ag9i55KYLUSH25j16Sbt0KaUCgpf2hS6G2OOaxUVEV9e55WItDPGHBCRLsBE4CygG3AD8Lhr2byH5+432XpUWf4cEV89yIt9h3Ndenc+/3Efl/Srd4FMKaXqxZeSwkoP21Y18LwfiMhm4FPgdtetp48DF4rINuBC13rLcdYd0GU4I7b9lXPbl/CnzzZTVOZzrZpSSjUKr0lBRNqJSH8gSkRSRaSf6zESqyNbvRljzjbG9DbG9DfGLHJtyzXGjDbG9HAtW9ZIcTY7THgRMZU8H/0P9uUV8cLXOkObUiqwaqsGugS4EUgBXqy2vQBo9A5tCojvBmP/TJtP7+LpUwZz37d2rhqcQvekVsGOTCnVQngtKRhj/mmMORu4yfXNvupxsTHmvQDG2LIMmgo9xjIh91XOCMvm4U83a6OzUipgxJcPHBEZC/TB6rwGgDHmz36MyydpaWkmPT092GE0voL98OKZ5FdGUlBcSkfbISQ2xbqVtXrjtFJK1YOIrDHGpHna50s/hRex7gb6DRAFXAec1qgRquO1bg+pk2hTupdOtlwEA3l74NO7YMO8YEenlGrGfLn7aKQx5hog1zU43jCsdgblT1s/P3FbebE1bpJSSvmJL0mhpGopIh1c6139FpGyeBkfyet2pZRqBL4khc9FpC3wJLAea1jr9/0ZlMLr+EgmtlOAA1FKtSR1zdFsAxYYY4647jjqBqQaY/5fQKJryTyMj1Ru7CzqeGuQAlJKtQR1zdFcCTxXbb24xXUqC5Ya4yMZRzQOcTJvYx7bcwqDHZ1SqpnypfroKxEZ7/dI1In6TYZ7N8KsI8h9Oyhvl8qT9hd56t0vcFZq3wWlVOPzJSncAXwkIsUickhEDouIlhYCzRGF4+q3iHTYuTPnYd76NiPYESmlmiFfkkIi4ABaAUmu9SY6kUETF9cVx+TXON22h7aL/pddB7UaSSnVuHyZec0JTALudz1PBgb4OzDlmfS4kKNn/S/jbctY9K/HqNRqJKVUI/KlR/MLwHnA9a5NRcDL/gxK1a71hb8nq925TD3yEqWPnwqz2sIzfbW3s1KqwXypPhpujPkVrk5srruPwv0alaqdzUbHoVdgE4gqywUdBkMp1Uh8SQrlrv4KBkBEEoBKv0al6iTfPoWNGlVHOgyGUqqBfEkKfwc+AJJE5GFgGfCEX6NSddNhMJRSflDnXMvGmDdFZA1wgWvTJGPMRv+GpeoUm2JVGXnarpRS9eRLSQHADpQDZSfxGuVPHobBMEBBv2lBCUcp1Tz4cvfRH4C5QEesIbPfEZHf+zswVYcaw2CUx3SgwERR+P0/qSzOC3Z0Sqkmqs6Z10QkAxhsjClyrUcDa4wxZwQgvlo125nX6mnh5x8wauXN7E46l+63fwQiwQ5JKRWCGjTzGrCL49sewoAdjRGYalyjL5rI/KRf0f3gYrI++0uww1FKNUG+JIUiYJOI/ENEXgV+BI6IyNMi8rR/w1MnQ0S4cPojLLSNpEP63yjMWBTskJRSTUyddx8Bn7keVVb4KRbVCGJjwkm45mV2vDmWDu9Nx3nnMuxxXYIdllKqifDlltQ5gQhENZ6Bp3Xm43NepMO3U8h/aSxxUTbIy7JuVx39oNVIrZRSHvhy99E4EVktIgd06Oym4/LR57I+4VLiyva6OrTpUBhKqbr50qbwAvAroBM6dHaTISKMdH5/4g4dCkMpVQtf2hQygfWuqTlVEyJ5WZ536FAYSikvfEkK9wGfisgSoLRqozFmtr+CUo1Eh8JQSp0kX6qPHgacQFusaqOqhwp1XobCIO3GoISjlAp9vpQU2hljBvs9EtX4qu4yWvQI5GWS70jAWVZM5LK/E9V3IsR1DWp4SqnQ40tJYZGInO/3SJR/9JsM926EWUeImrmNxzs8Q2lJMUWvjYfCnGBHp5QKMb4khVuAhSJSqLekNm0Ou40Hpk/k0TazkPy9FL0+EUoLgh2WUiqE+FJ9lOj3KFTAtI50cP8tU3nw+cP85eDjlP3faMKdR7Vzm1IK8KGkYIxxApOA+13Pk4EB/g5M+U+7NpHcOuN2PuJ8wg9t1c5tSik3X3o0vwCcB1zv2lQEvNyQk4rIvSKySUQ2ishcEYkUkXgR+UpEtrmWcQ05h6rdqUmtuDwm48Qd2rlNqRbNlzaF4caYXwElAMaYQ0B4fU8oIp2Au4A0Y0xfrFndpgAzgUXGmB7AIte68qOIo3s979DObUq1WL4khXIRseG6xV1EEoCG9m4OA6JEJAyIBvYC44E3XPvfACY08ByqLt46scV2CmwcSqmQ4TUpuD6wAf4OfAAkicjDwDLgifqe0BiTBTwJ7AaygTxjzJdAe2NMtuuYbKCdl7hmiEi6iKTn5OgtlQ3ioXMbgIlOgopSDy9QSjV3tZUUVgEYY94EHsD6ID8MTDLGvFvfE7raCsYD3bDmfY4Rket8fb0x5hVjTJoxJi0pSTtWN0iNeZ7zIzrwQcUIJHsdZu41UFYU7AiVUgFW2y2p7gl+jTGbgE2NdM4LgJ3GmBwAEfkQGA7sF5FkY0y2iCQDBxrpfKo2/Sa7b0FtbQw7v/yJ+755lce3/wPevhK5+t8Q2SbIQSqlAqW2pJAkIr/xttMYU9+pOHcDZ4pINFAMjAbSgaPADcDjruXH9Xx/VU8iwm/H9OSJypu5+9sInt31ErY3xyODrodvn7YaoLUvg1LNWm1JwQ60olqJoTEYY1aKyPvAWqACWAe84jrXPBG5CStxTGrM8yrfiAj3jzudv5jrmLEsglf2Pott7zrEus/gWF8G0MSgVDMkxhjPO0TWGmMGBTiek5KWlmbS09ODHUazZIzh0f9kcFv6WJIk/8QDYjtbYyoppZocEVljjEnztK+2huZGLSGopkVE+OOlZ5AoXsZG0r4MSjVLtSWF0QGLQoUkEamlL4NO1KNUc+Q1Kbh6LqsWTrz0ZSD+VHCWBz4gpZRf+dKjWbVk1foyGIQsk8gKexrsXAL/ugKO5gY7QqVUI/Jl6GzV0rn6MgiwY1sOv/rXGq6J/I4/7HkZefU8GDwN0l/TW1aVaga0pKBOytk9knjnljP5oGIk03iIiqOHYNHD1q2qOvy2Uk2eJgV10gZ0bst7tw7nl8jeHCxznHiADr+tVJOlSUHVy2ntWjH/1yNoJ4c9H6C3rCrVJGlSUPUWFxOOeBtmOyoOvHSMVEqFLk0KqkFk9EMn3LJqxAbFh+D9G6HYS0lCKRWS9O4j1TBVdxktegSTl0k2CfzdTOHm/uF0+3E27FkFV7wMBdlWO4PeoaRUSPM69lFToGMfhZ6dB4/yq3+ls+1AIY+mlXFt1qPIoe1gC4PKimMHOutA2R0AABSiSURBVKKs/g+aGJQKuPqOfaTUSeuWGMPHt49k0uAUHlgdzlTHU1Q6Yo5PCKB3KCkVojQpqEYXFW7nr1f15+nJ/UnfWwrlXmZw0zuUlAo52qag/GbioBT6pcSS81Ii7Y2H+bSj4q07lEQH5FUqVGhJQfnVae1aE3fZnyiTiBp7BIpz4c3xcHCb1QP6mb4wq6211B7RSgWFlhSU34UPmgJhNooWPERkcTbZJoGf+tzNuadEYfv6UXhhKNhsx9oddHY3pYJGSwoqMPpNJvr+DA7+dh8PdX+X6WtP5X/W9WHPdd+AI1IbopUKEZoUVEC1ax3Jq1PTeHJSf7bsK2DMK1sw5cWeD9aGaKUCTpOCCjgR4arBKXx57zmkdY0jqzLB84Fh4ZC1JrDBKdXCaVJQQZMcG8WbNw5l98DfUUz48TttDpAwePV8ePdaWPasNkQrFQCaFFRQiQjDr7iN8ouf41BYeyqNsI8kfhzyF/jdVhj1/2DbV7DwIZ2zQakA0KSgQkKbodcQ/8BPrL5hO9e3mcNlSzty87tbyex/J8R4qF7Shmil/EKTggopw7on8PndZzPzol4s//kgo59aisnP9nxw3h7Y/rUO0a1UI9J+CirkOOw2bj33VC7r35En/7uVrE0JpNgOnnig2OBfV0D7VBh+J5hKWPyYjsSqVAPoKKkq5O1e8jrtl95HhCl1bzOOKOSSp61E8N3zkLMFEKDa37OOxKqURzpKqmrSuoyaRvgVL1Ac3ZFKhMzKRGZH3cH3rcfCwOvgtu8hOpHjEgJou4NS9aDVR6pJkH6Tieo3mQpnJUvT9/D2wm088+oKzuwez92je3JWUa7nF+btgeXPQf+rYccSnehHqTpo9ZFqkkrKncxdtZuXlmznQEEpq2PuIcl54MQD7eHgLAPEaoMwzmP7tHpJtVBafaSanUiHnekjuvHNfefx0GW9mW2upsjU6ADniILxf4fbV0F4q+MTAmj1klIeaFJQTVpVcvjD7x9iVeosskmk0gj7JYkVfR6itPeVkHQ6lBV6foO8PfDVg7B3nXVrqw7hrVo4rT5SzUpZRSWf/LCXV7/Zwdb9BbRrHcENw7ty27rx2PI9DLAX5hqhtbLCaqwuOaJzSatmT6uPVIsRHmbjqsEpfHHP2bx541BO79Cav/13K/cdmXDiRD+OKLj8efjdNmtZVuB5CO+FDx9b15KEaua0pKCavYzsfOYs24nZMI975V062nIpjuqAY8wswgdOOXbgrLaccFtrlR5jIDoBNn0EFSXHtmtJQjVBtZUUAp4UROR04N/VNnUHHgTedG3vCvwCTDbGHK7tvTQpqJNxpKiM99dk8s7K3ew4eJS20Q4mDU7h6qFd6J7Uyvrmn7fnxBeGt4KYJDi80/Mbx3aGezceW98wT299VSEtpJLCcScXsQNZwDDgduCQMeZxEZkJxBlj7q/t9ZoUVH0YY/h+ey5vrdzFl5v2U1FpGNilLb/p8AMjNz+CVFSb9KeqJJA6CR6Ow2tJYvhd0PVsKMiGL+63qp1qvocmBhUiQjkpjAEeMsaMEJGtwChjTLaIJANLjDGn1/Z6TQqqoQ7klzB/fRYfrMli6/4CJoZ9xx8i5xFfkQOxnZDRDx37MPdWkrBHAMbVH8KL6qUJLUmoIAvlpPAasNYY84KIHDHGtK2277AxJq6212tSUI3FGMOmvfl8sDaTT9bvJfdoGfEx4Yzt055LUjtyZvd4wja9b83j4KkU0OtSyFwNb17u/SRj/wwlebB8NngqjWhiUAESkklBRMKBvUAfY8x+X5OCiMwAZgB06dJl8K5duwIWs2oZyp2VLN2awyc/7GVRxn6OljndCeKGVqs5fdPTSF6W52/53koTYj+x81x1bVLgN5uOrWtpQvlRqCaF8cDtxpgxrnWtPlIhp6TcyZKtOXz+YzYLM/ZTVOYkLtrBeb3accEZ7TmnZxKtIqoNIbZhnvfSRLdz4ame3k+W1As6pEKlE7Z8Bs7SE9+jKjFo0lANUFtSCOaAeFcDc6utfwLcADzuWn4cjKCUqi7SYWdc3w6M69vBnSC+2JjNoowDfLg2C4ddOLN7Ahec0Z7RZ7QjpeqD2dsHdmxnzyWJiDYQ3x12r/C8v7wYPv9fiI6Hg9tg0cPHEk/V9KRwfGLQxKHqISglBRGJBvYA3Y0xea5tCcA8oAuwG5hkjDlU2/toSUEFS4WzkjW7DrMwYz+LMg6w4+BRAE5r14qRpyVyTs9EhnVLICaixveu2koSVR/YtfWXqE1MO7jtO4hJhB/fq/s8mjRarJCsPmoMmhRUqNiRU8jXWw7wzbaDrNyRS2lFJQ67MPiUOM7ukcTwUxPo2ykWh91W94ext3aJNh1h4qvw+iW1BxPe2upgV1nu4T06wW82+5acNGk0W5oUlAqgknIn6b8c5tttOXy77SCbs/MBiA63M/iUOM7snsCwbvH0S2lLeJiHkWbq+sD2ljSiE+Gc38GhnbDq/7wHGNkWyo56Thqtk+HONVabRl1JoypWTRxNjiYFpYLoYGEpK3ccYuXOXFbuOMTW/QUARDpsDOjclsGnxDGoSxwDu8QRH+Ma/ru2D1tfvuV7SxyRsdD3KkifU0fUNaY2rRIVD1fNsXp4714BX/2x4aUNTSwBp0lBqRBy6GgZq3bmsmLHIdbsOszm7Hycldb/YbfEGAZ2acvALnH0T4nl9A6tiQizn/gmdX2Q1re0ERUPw++0GrLryxEDw2bAkT2Q8cnxnfrCIuHSZ2HA1b7F6cvP6uv10MTjpklBqRBWXOZkQ+YR1u4+wtrdh1m3+zAHC60PUodd6Nm+Nf1SYunbKZbUTrH0bN+aSIeHRFFTQ0ob3pJG6w5w1T+h8AC8d4P3c9vCThxxtjpHDES1td7HUzVWZCxc8DDs3whr36yRWKLg0qdhwDW+/Sy+JJ66rpcv+xvrPQJAk4JSTYgxhszDxWzMymNDVp61zMwjr9j68LTbhO6JMfRKbsMZya05o0MbzkhuQ/s2EYiI7yfyVxVVbGe458fax4o66w4oPgLr3/I93ppsDohoBSX5njsGOmJg0FTrHKUFJ+6PToRJr4MjGnYuhaVPeB8B19cSTUOTU4CShiYFpZq4qkTxY1YeGdn5ZGQXkJGdT9aRYx8wbSLD6NG+NT3ateK0dq3cz5NjI08uWVTxVxVV9XGgvN5p1QluXghP98ZrYhlxj9VgvvpV7z9DRBsoza/zR61VRBtr5j5TeeI+ewR0P9eaC3z7ouOvRfXXj7gLvnveGuakpugEaz6PXSusGwRq67TYSDQpKNVM5RWXs3WflSC2HShg2/5Cfj5QSO7RY9Ut0eF2uibE0C0phlMTrWW3xFZ0S4ghNtrRsAAaWtrwZ2KpOsbb/ph2VqN5eTG8U8uH7pm/hhUvet+fPMCq3jqw2fsxDVFzaPZGEKo9mpVSDRQb5WBot3iGdos/bntuYSk/HyjkpwOF7Mw5ys6DhWzKyuOLjfvcjdpglS66JERzSnwMneOj6eJ6dIqLomPbSM+N3NX1m+z9W2xdvbt9OWb0g56TxugHj63XdYy3/WMfg27nuC6kl57msZ1h3F8g41Pv+3+11HpeW3K6cy3MHgD5WSfub9UerpkHr4zCY6koz8M0sn6kSUGpZiihVQQJrSIY1j3huO1lFZXsPlTEjpxCdh8qYlduEbsPFZGRnc+Xm/dR7jz+QympdQSd2kaREhdlJYrYKDrERpIcG0mH2EgSYyKw2WqpmqotafhyTGMkFl/eo76JxdfkFBYOF8zyvH/Mn6DjACsuj0klxfO18ROtPlJKAeCsNOzLL2F3bhFZR4rJOlxM1pEiMg8Xk3WkmL1Hik9IGg670L5NJO3bRNKudYT1qHreJpKkVhEktg4nISYCe23JIxQE++4jX++SagTapqCUarDKSkPu0TL25ZWQnVfMvvwSsvNKyD5SzIGCUvbnl3CgoJSCkhNvRRWB+OhwEl1JIrFVBPEx4cRHhxPfKpyEmHDiYyKIj3HQNjqctlEOwuweens3d3r3UcNoUlAq9BSXOTlQUML+/FIOFroeBaXkFJYdWy8s5fDRcgpLvfdlaB0ZRttoB3HR4bSNDic2ykFsVJhreezRJtJBmygHrSPDaB1pLR0tMaGcBG1oVkoFTFS4nVMSYjglIabOY0vKnRwuKuPQ0WOPI0XlHC6ylkeKyjjsWu7OPUpecTn5JRXHNZZ7jMFhp3VkGK0iw2gVUe3hWo+JCCMm3G4tI8KICQ8jJsJOdHgY0eF2osLtxISHERVuJzrc3qKSjCYFpVTQRDrsJMdGkRwb5fNrjDEUllaQV1xOXnE5BSUV5LuWBSVW0igosUohBSUVFJZWUFhSwe6jRdbz0gqKSp2UOT30O/DCYRciHXaiHFbCqFpGhtmJdNjczyNcx0Q6bESEVS1tRFTbFm63EVH1PMzaHx5ms7aH2XDYXethNsJsUr8+Jg2gSUEp1aSIiKuayEFKrbO4166sopKjpRUcLavgaKmTo2UVFJc5OVpaQXG5k6OlTorKKigqc1Jc7qS4zElJ+bHnxeXWek5hOSXllZS41queV9RRmvHtZ8VKEq5E4bCLe/38Xu144NLeDT5HTZoUlFItkvVtPJy4qpFpG1mFs5IyZyUl5ZWUVjgpLa+kpMJJWUUlpRWVruWx9apt5U5rWVZhvb7MWUl5haHc6drn2p/c1vfS1cnQpKCUUn4QZrcRZrcR7Z+c4zctp/VEKaVUnTQpKKWUctOkoJRSyk2TglJKKTdNCkoppdw0KSillHLTpKCUUspNk4JSSim3Jj1KqojkALv8eIpE4KAf37+xNJU4oenEqnE2rqYSJzSdWBsS5ynGmCRPO5p0UvA3EUn3NrxsKGkqcULTiVXjbFxNJU5oOrH6K06tPlJKKeWmSUEppZSbJoXavRLsAHzUVOKEphOrxtm4mkqc0HRi9Uuc2qaglFLKTUsKSiml3DQpKKWUctOkUI2I/E1EtojIBhH5SETaejnuFxH5UUTWi0h6AOMbJyJbReRnEZnpYb+IyGzX/g0iMihQsVWLobOILBaRDBHZJCJ3ezhmlIjkua7fehF5MNBxVoul1t9liFzT06tdq/Uiki8i99Q4JijXVEReE5EDIrKx2rZ4EflKRLa5lh4nzazr7zlAsYbc/7yXOGeJSFa13+/FXl7b8GtqjNGH6wGMAcJcz58AnvBy3C9AYoBjswPbge5AOPAD0LvGMRcDCwABzgRWBuEaJgODXM9bAz95iHMU8J9g/759+V2GwjX18HewD6vzUdCvKXAOMAjYWG3bX4GZruczPf0f+fL3HKBYQ+5/3kucs4Df+fC30eBrqiWFaowxXxpjKlyrK4CUYMZTw1DgZ2PMDmNMGfAuML7GMeOBN41lBdBWRJIDGaQxJtsYs9b1vADIADoFMoZGFvRrWsNoYLsxxp89+X1mjPkGOFRj83jgDdfzN4AJHl7qy99zo/IUayj+z3u5pr5olGuqScG7G7G+IXpigC9FZI2IzAhQPJ2APdXWMznxw9aXYwJGRLoCA4GVHnafJSI/iMgCEekT0MCOV9fvMqSuKTAFmOtlX6hc0/bGmGywviQA7TwcE2rXFULvf76mO1zVXK95qZJrlGsaVt/omioRWQh08LDrD8aYj13H/AGoAN728jYjjDF7RaQd8JWIbHFld38SD9tq3k/syzEBISKtgA+Ae4wx+TV2r8Wq/ih01Y3OB3oEOkaXun6XoXRNw4HLgd972B1K19QXIXNdIWT/56t7CXgU6xo9CjyFlcSqa5Rr2uJKCsaYC4wxfT08qhLCDcClwLXGVVHn4T32upYHgI+wim3+lgl0rraeAuytxzF+JyIOrITwtjHmw5r7jTH5xphC1/PPAYeIJAY4zKpY6vpdhsQ1dbkIWGuM2V9zRyhdU2B/VRWba3nAwzEhc11D+H+++vn3G2OcxphK4FUv52+Ua9rikkJtRGQccD9wuTGmyMsxMSLSuuo5VkPVRk/HNrLVQA8R6eb6xjgF+KTGMZ8AU113zJwJ5FUV4wNFRASYA2QYY572ckwH13GIyFCsv8PcwEXpjsOX32XQr2k1V+Ol6ihUrqnLJ8ANruc3AB97OMaXv2e/C/H/+eoxVG/HusLL+RvnmgaiNb2pPICfserk1rseL7u2dwQ+dz3vjtWq/wOwCavaKVDxXYx1N8/2qvMCtwK3up4L8HfX/h+BtCBcw5FYRdYN1a7jxTXivMN17X7AatwbHqTft8ffZahdU1cc0Vgf8rHVtgX9mmIlqWygHOub6k1AArAI2OZaxruOdf8feft7DkKsIfc/7yXOf7n+/jZgfdAn++ua6jAXSiml3LT6SCmllJsmBaWUUm6aFJRSSrlpUlBKKeWmSUEppZSbJgWlXESkrYj8utp6RxF530/nmlA1kqmIvC4iV3k4JklEvvDH+ZXyRpOCUse0BdxJwRiz1xhzwod1I7kPeLG2A4wxOUC2iIzwUwxKnUCTglLHPA6c6hqv/m8i0rVqTHsRmSYi80XkUxHZKSJ3iMhvRGSdiKwQkXjXcaeKyBeugdO+FZFeNU8iIj2BUmPMwWqbzxGR70RkR41Sw3zgWj/+zEodR5OCUsfMxBqWeoAx5n897O8LXIM17sxjQJExZiDwPTDVdcwrwJ3GmMHA7/BcGhiBNYBddclYvcEvxUpOVdKBs+v34yh18lrcKKlKNcBiY80RUSAiecCnru0/Av1cI8MOB95zDUMEEOHhfZKBnBrb5htrsLPNItK+2vYDWEMZKBUQmhSU8l1pteeV1dYrsf6XbMARY8yAOt6nGIit5b2rD4Ec6TpeqYDQ6iOljinAmkK0Xow1b8ROEZkE7vmd+3s4NAM4zce37UmAR+RULZsmBaVcjDG5wHIR2Sgif6vn21wL3CQiVSNqepoO8RtgoFSrY6rFecBn9YxFqZOmo6QqFQQi8hzwqTFmYR3HfQOMN8YcDkxkqqXTkoJSwfFnrDkSvBKRJOBpTQgqkLSkoJRSyk1LCkoppdw0KSillHLTpKCUUspNk4JSSik3TQpKKaXc/j8qndigxHbpSQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "t=np.linspace(-3,15)\n",
    "t_1=np.arange(-3,15,time_step)\n",
    "Temp_ambient = 65\n",
    "Temp_t1 = 85\n",
    "plt.plot(t,T_analytical(Temp_ambient,Temp_t1,K,t),'-',label='Analytical')\n",
    "plt.plot(t_1,T_numerical,'o-',label='Numerical')\n",
    "plt.legend()\n",
    "plt.xlabel('time (h)')\n",
    "plt.ylabel('Temperature (F)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}