TugJulia/proto/FTCS.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib notebook\n",
    "%matplotlib inline\n",
    "from ipywidgets import interact, IntSlider"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up environment\n",
    "dim = 2\n",
    "grid_size = {'x': 20, 'y': 20} # to start out\n",
    "domain_size = {'x': 20, 'y': 20} # to start out\n",
    "\n",
    "# boundary_condition\n",
    "# constant boundary condition\n",
    "bc_left = 0\n",
    "bc_right = 0\n",
    "bc_top = 0\n",
    "bc_bottom = 0\n",
    "\n",
    "delta_x = domain_size['x'] / grid_size['x']\n",
    "delta_y = domain_size['y'] / grid_size['y'] \n",
    "\n",
    "alpha_x = np.ones((grid_size['x'],grid_size['y']))\n",
    "alpha_y = np.ones((grid_size['x'],grid_size['y']))\n",
    "\n",
    "# time dimension\n",
    "max_stable_time_step_x = delta_x**2 / (2 * np.max(alpha_x))\n",
    "max_stable_time_step_y = delta_y**2 / (2 * np.max(alpha_y))\n",
    "max_stable_time_step = max_stable_time_step_x if max_stable_time_step_x <= max_stable_time_step_y else max_stable_time_step_y\n",
    "\n",
    "time_step = 0.1\n",
    "if time_step > max_stable_time_step:\n",
    "    print('Warning: Time step more than maximum stable time step!')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAGdCAYAAABKG5eZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAht0lEQVR4nO3dfWxUZeK38e/BwoCmHa3SzgyUWkkVeQmLgIWqtehSLCsLglolgRI3KhFdsSFqVWLdbBjxtxqCRY2uggQFsinQ7oIrJdBWlkJAKboG2RKrZaWzXYjOQJVpK+f5w4dxx77gwIztPVyf5CSeM/d9vGcy8fJ0zrSWbdu2AAAwRJ+eXgAAAJEgXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMktDTC4iW06dP6+jRo0pMTJRlWT29HABABGzb1okTJ+TxeNSnT/fXVHETrqNHjyotLa2nlwEAOA9HjhzR4MGDux0TN+FKTEyUJN2oqUpQ3x5eDQAgEu1q005tCf23vDtxE64zPx5MUF8lWIQLAIzy/39r7s/5qIebMwAARolZuF555RVlZGSof//+Gjt2rD744INux1dXV2vs2LHq37+/rrrqKr322muxWhoAwGAxCdf69eu1cOFCPf3009q/f79uuukm5efnq7GxsdPxDQ0Nmjp1qm666Sbt379fTz31lH7/+9+rrKwsFssDABjMisXf48rKytJ1112nV199NXTs2muv1YwZM+T1ejuMf+KJJ1RRUaGDBw+Gjs2fP18HDhxQbW3tz/p3BgIBOZ1O5Wo6n3EBgGHa7TZVqVx+v19JSUndjo36FVdra6s+/PBD5eXlhR3Py8vTrl27Op1TW1vbYfyUKVO0b98+tbW1dTonGAwqEAiEbQCA+Bf1cB07dkzff/+9UlNTw46npqbK5/N1Osfn83U6vr29XceOHet0jtfrldPpDG18hwsALgwxuznjp7c02rbd7W2OnY3v7PgZxcXF8vv9oe3IkSPnuWIAgAmi/j2uK664QhdddFGHq6vm5uYOV1VnuFyuTscnJCTo8ssv73SOw+GQw+GIzqIBAMaI+hVXv379NHbsWFVWVoYdr6ysVHZ2dqdzJk6c2GH81q1bNW7cOPXty40WAIAfxeRHhUVFRfrzn/+st956SwcPHtRjjz2mxsZGzZ8/X9IPP+abO3duaPz8+fP15ZdfqqioSAcPHtRbb72lN998U4sWLYrF8gAABovJr3wqKCjQ8ePH9Yc//EFNTU0aOXKktmzZovT0dElSU1NT2He6MjIytGXLFj322GNasWKFPB6Pli9frlmzZsVieQAAg8Xke1w9ge9xAYC5evR7XAAAxBLhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYJaGnFwBcKN4/Wnfe55ji+dV5nwMwHVdcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGCUqIfL6/Vq/PjxSkxMVEpKimbMmKFDhw51O6eqqkqWZXXYPvvss2gvDwBguKiHq7q6WgsWLNDu3btVWVmp9vZ25eXlqaWl5axzDx06pKamptCWmZkZ7eUBAAwX9T9r8ve//z1sf+XKlUpJSdGHH36onJycbuempKTo0ksvjfaSAABxJOZ/j8vv90uSkpOTzzp2zJgxOnXqlIYPH65nnnlGkyZN6nJsMBhUMBgM7QcCgfNfLBBD/C0tIDpienOGbdsqKirSjTfeqJEjR3Y5zu126/XXX1dZWZk2bNiga665Rrfeeqtqamq6nOP1euV0OkNbWlpaLJ4CAKCXsWzbtmN18gULFmjz5s3auXOnBg8eHNHcadOmybIsVVRUdPp4Z1dcaWlpytV0JVh9z2vdAIBfVrvdpiqVy+/3KykpqduxMbvieuSRR1RRUaEdO3ZEHC1JmjBhgurr67t83OFwKCkpKWwDAMS/qH/GZdu2HnnkEW3cuFFVVVXKyMg4p/Ps379fbrc7yqsDAJgu6uFasGCB3n33XZWXlysxMVE+n0+S5HQ6NWDAAElScXGxvvrqK61evVqStGzZMl155ZUaMWKEWltbtWbNGpWVlamsrCzaywMAGC7q4Xr11VclSbm5uWHHV65cqXnz5kmSmpqa1NjYGHqstbVVixYt0ldffaUBAwZoxIgR2rx5s6ZOnRrt5QEADBfTmzN+SYFAQE6nk5szAMBAveLmDAAAYoFwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBglKiHq6SkRJZlhW0ul6vbOdXV1Ro7dqz69++vq666Sq+99lq0lwUAiBMJsTjpiBEjtG3bttD+RRdd1OXYhoYGTZ06Vffff7/WrFmjf/zjH3rooYc0cOBAzZo1KxbLAwAYLCbhSkhIOOtV1hmvvfaahgwZomXLlkmSrr32Wu3bt09/+tOfCBcAoIOYfMZVX18vj8ejjIwM3XPPPfr888+7HFtbW6u8vLywY1OmTNG+ffvU1tbW5bxgMKhAIBC2AQDiX9TDlZWVpdWrV+v999/XG2+8IZ/Pp+zsbB0/frzT8T6fT6mpqWHHUlNT1d7ermPHjnX57/F6vXI6naEtLS0tqs8DANA7RT1c+fn5mjVrlkaNGqVf//rX2rx5syTp7bff7nKOZVlh+7Ztd3r8fxUXF8vv94e2I0eORGH1AIDeLiafcf2vSy65RKNGjVJ9fX2nj7tcLvl8vrBjzc3NSkhI0OWXX97leR0OhxwOR1TXCgDo/WL+Pa5gMKiDBw/K7XZ3+vjEiRNVWVkZdmzr1q0aN26c+vbtG+vlAQAME/VwLVq0SNXV1WpoaNCePXt05513KhAIqLCwUNIPP+KbO3duaPz8+fP15ZdfqqioSAcPHtRbb72lN998U4sWLYr20gAAcSDqPyr897//rXvvvVfHjh3TwIEDNWHCBO3evVvp6emSpKamJjU2NobGZ2RkaMuWLXrssce0YsUKeTweLV++nFvhAQCdsuwzd0IYLhAIyOl0KlfTlWDxI0YAMEm73aYqlcvv9yspKanbsfyuQgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGIVwAAKMQLgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGIVwAAKMQLgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGIVwAAKMQLgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGiXq4rrzySlmW1WFbsGBBp+Orqqo6Hf/ZZ59Fe2kAgDiQEO0T7t27V99//31o/5///KcmT56su+66q9t5hw4dUlJSUmh/4MCB0V4aACAORD1cPw3O888/r6FDh+rmm2/udl5KSoouvfTSaC8HABBnYvoZV2trq9asWaP77rtPlmV1O3bMmDFyu9269dZbtWPHjrOeOxgMKhAIhG0AgPgX03Bt2rRJ33zzjebNm9flGLfbrddff11lZWXasGGDrrnmGt16662qqanp9txer1dOpzO0paWlRXn1AIDeyLJt247VyadMmaJ+/frpr3/9a0Tzpk2bJsuyVFFR0eWYYDCoYDAY2g8EAkpLS1OupivB6nvOawYA/PLa7TZVqVx+vz/sfofORP0zrjO+/PJLbdu2TRs2bIh47oQJE7RmzZpuxzgcDjkcjnNdHgDAUDH7UeHKlSuVkpKi3/zmNxHP3b9/v9xudwxWBQAwXUyuuE6fPq2VK1eqsLBQCQnh/4ri4mJ99dVXWr16tSRp2bJluvLKKzVixIjQzRxlZWUqKyuLxdIAAIaLSbi2bdumxsZG3XfffR0ea2pqUmNjY2i/tbVVixYt0ldffaUBAwZoxIgR2rx5s6ZOnRqLpQEADBfTmzN+SYFAQE6nk5szAMBAkdycwe8qBAAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGCUiMNVU1OjadOmyePxyLIsbdq0Kexx27ZVUlIij8ejAQMGKDc3V59++ulZz1tWVqbhw4fL4XBo+PDh2rhxY6RLAwBcACIOV0tLi0aPHq3S0tJOH3/hhRf00ksvqbS0VHv37pXL5dLkyZN14sSJLs9ZW1urgoICzZkzRwcOHNCcOXN09913a8+ePZEuDwAQ5yzbtu1znmxZ2rhxo2bMmCHph6stj8ejhQsX6oknnpAkBYNBpaamaunSpXrwwQc7PU9BQYECgYDee++90LHbbrtNl112mdauXfuz1hIIBOR0OpWr6Uqw+p7rUwIA9IB2u01VKpff71dSUlK3Y6P6GVdDQ4N8Pp/y8vJCxxwOh26++Wbt2rWry3m1tbVhcyRpypQp3c4JBoMKBAJhGwAg/kU1XD6fT5KUmpoadjw1NTX0WFfzIp3j9XrldDpDW1pa2nmsHABgipjcVWhZVti+bdsdjp3vnOLiYvn9/tB25MiRc18wAMAYCdE8mcvlkvTDFZTb7Q4db25u7nBF9dN5P726Otsch8Mhh8NxnisGAJgmqldcGRkZcrlcqqysDB1rbW1VdXW1srOzu5w3ceLEsDmStHXr1m7nAAAuTBFfcZ08eVKHDx8O7Tc0NKiurk7JyckaMmSIFi5cqCVLligzM1OZmZlasmSJLr74Ys2ePTs0Z+7cuRo0aJC8Xq8k6dFHH1VOTo6WLl2q6dOnq7y8XNu2bdPOnTuj8BQBAPEk4nDt27dPkyZNCu0XFRVJkgoLC7Vq1So9/vjj+u677/TQQw/p66+/VlZWlrZu3arExMTQnMbGRvXp8+PFXnZ2ttatW6dnnnlGixcv1tChQ7V+/XplZWWdz3MDAMSh8/oeV2/C97gAwFw99j0uAABijXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGCUiMNVU1OjadOmyePxyLIsbdq0KfRYW1ubnnjiCY0aNUqXXHKJPB6P5s6dq6NHj3Z7zlWrVsmyrA7bqVOnIn5CAID4FnG4WlpaNHr0aJWWlnZ47Ntvv9VHH32kxYsX66OPPtKGDRv0r3/9S7/97W/Pet6kpCQ1NTWFbf379490eQCAOJcQ6YT8/Hzl5+d3+pjT6VRlZWXYsZdfflnXX3+9GhsbNWTIkC7Pa1mWXC5XpMsBAFxgYv4Zl9/vl2VZuvTSS7sdd/LkSaWnp2vw4MG6/fbbtX///m7HB4NBBQKBsA0AEP9iGq5Tp07pySef1OzZs5WUlNTluGHDhmnVqlWqqKjQ2rVr1b9/f91www2qr6/vco7X65XT6QxtaWlpsXgKAIBexrJt2z7nyZaljRs3asaMGR0ea2tr01133aXGxkZVVVV1G66fOn36tK677jrl5ORo+fLlnY4JBoMKBoOh/UAgoLS0NOVquhKsvhE/FwBAz2m321Slcvn9/rP2IuLPuH6OtrY23X333WpoaND27dsjipYk9enTR+PHj+/2isvhcMjhcJzvUgEAhon6jwrPRKu+vl7btm3T5ZdfHvE5bNtWXV2d3G53tJcHADBcxFdcJ0+e1OHDh0P7DQ0NqqurU3Jysjwej+6880599NFH+tvf/qbvv/9ePp9PkpScnKx+/fpJkubOnatBgwbJ6/VKkp577jlNmDBBmZmZCgQCWr58uerq6rRixYpoPEcAQByJOFz79u3TpEmTQvtFRUWSpMLCQpWUlKiiokKS9Ktf/Sps3o4dO5SbmytJamxsVJ8+P17sffPNN3rggQfk8/nkdDo1ZswY1dTU6Prrr490eQCAOHdeN2f0JoFAQE6nk5szAMBAkdycwe8qBAAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGCUiMNVU1OjadOmyePxyLIsbdq0KezxefPmybKssG3ChAlnPW9ZWZmGDx8uh8Oh4cOHa+PGjZEuDQBwAYg4XC0tLRo9erRKS0u7HHPbbbepqakptG3ZsqXbc9bW1qqgoEBz5szRgQMHNGfOHN19993as2dPpMsDAMS5hEgn5OfnKz8/v9sxDodDLpfrZ59z2bJlmjx5soqLiyVJxcXFqq6u1rJly7R27dpIlwgAiGMx+YyrqqpKKSkpuvrqq3X//ferubm52/G1tbXKy8sLOzZlyhTt2rWryznBYFCBQCBsAwDEv6iHKz8/X++88462b9+uF198UXv37tUtt9yiYDDY5Ryfz6fU1NSwY6mpqfL5fF3O8Xq9cjqdoS0tLS1qzwEA0HtF/KPCsykoKAj988iRIzVu3Dilp6dr8+bNmjlzZpfzLMsK27dtu8Ox/1VcXKyioqLQfiAQIF4AcAGIerh+yu12Kz09XfX19V2OcblcHa6umpubO1yF/S+HwyGHwxG1dQIAzBDz73EdP35cR44ckdvt7nLMxIkTVVlZGXZs69atys7OjvXyAACGifiK6+TJkzp8+HBov6GhQXV1dUpOTlZycrJKSko0a9Ysud1uffHFF3rqqad0xRVX6I477gjNmTt3rgYNGiSv1ytJevTRR5WTk6OlS5dq+vTpKi8v17Zt27Rz584oPEUAQDyJOFz79u3TpEmTQvtnPmcqLCzUq6++qk8++USrV6/WN998I7fbrUmTJmn9+vVKTEwMzWlsbFSfPj9e7GVnZ2vdunV65plntHjxYg0dOlTr169XVlbW+Tw3AEAcsmzbtnt6EdEQCATkdDqVq+lKsPr29HIAABFot9tUpXL5/X4lJSV1O5bfVQgAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKBGHq6amRtOmTZPH45FlWdq0aVPY45Zldbr93//9X5fnXLVqVadzTp06FfETAgDEt4jD1dLSotGjR6u0tLTTx5uamsK2t956S5ZladasWd2eNykpqcPc/v37R7o8AECcS4h0Qn5+vvLz87t83OVyhe2Xl5dr0qRJuuqqq7o9r2VZHeYCAPBTMf2M6z//+Y82b96s3/3ud2cde/LkSaWnp2vw4MG6/fbbtX///m7HB4NBBQKBsA0AEP9iGq63335biYmJmjlzZrfjhg0bplWrVqmiokJr165V//79dcMNN6i+vr7LOV6vV06nM7SlpaVFe/kAgF7Ism3bPufJlqWNGzdqxowZnT4+bNgwTZ48WS+//HJE5z19+rSuu+465eTkaPny5Z2OCQaDCgaDof1AIKC0tDTlaroSrL4R/fsAAD2r3W5Tlcrl9/uVlJTU7diIP+P6uT744AMdOnRI69evj3hunz59NH78+G6vuBwOhxwOx/ksEQBgoJj9qPDNN9/U2LFjNXr06Ijn2raturo6ud3uGKwMAGCyiK+4Tp48qcOHD4f2GxoaVFdXp+TkZA0ZMkTSDz+2+8tf/qIXX3yx03PMnTtXgwYNktfrlSQ999xzmjBhgjIzMxUIBLR8+XLV1dVpxYoV5/KcAABxLOJw7du3T5MmTQrtFxUVSZIKCwu1atUqSdK6detk27buvffeTs/R2NioPn1+vNj75ptv9MADD8jn88npdGrMmDGqqanR9ddfH+nyAABx7rxuzuhNAoGAnE4nN2cAgIEiuTmD31UIADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCgRhcvr9Wr8+PFKTExUSkqKZsyYoUOHDoWNsW1bJSUl8ng8GjBggHJzc/Xpp5+e9dxlZWUaPny4HA6Hhg8fro0bN0b2TAAAF4SIwlVdXa0FCxZo9+7dqqysVHt7u/Ly8tTS0hIa88ILL+ill15SaWmp9u7dK5fLpcmTJ+vEiRNdnre2tlYFBQWaM2eODhw4oDlz5ujuu+/Wnj17zv2ZAQDikmXbtn2uk//73/8qJSVF1dXVysnJkW3b8ng8WrhwoZ544glJUjAYVGpqqpYuXaoHH3yw0/MUFBQoEAjovffeCx277bbbdNlll2nt2rU/ay2BQEBOp1O5mq4Eq++5PiUAQA9ot9tUpXL5/X4lJSV1O/a8PuPy+/2SpOTkZElSQ0ODfD6f8vLyQmMcDoduvvlm7dq1q8vz1NbWhs2RpClTpnQ7JxgMKhAIhG0AgPh3zuGybVtFRUW68cYbNXLkSEmSz+eTJKWmpoaNTU1NDT3WGZ/PF/Ecr9crp9MZ2tLS0s71qQAADHLO4Xr44Yf18ccfd/qjPMuywvZt2+5w7HznFBcXy+/3h7YjR45EsHoAgKkSzmXSI488ooqKCtXU1Gjw4MGh4y6XS9IPV1Butzt0vLm5ucMV1f9yuVwdrq7ONsfhcMjhcJzL8gEABovoisu2bT388MPasGGDtm/froyMjLDHMzIy5HK5VFlZGTrW2tqq6upqZWdnd3neiRMnhs2RpK1bt3Y7BwBwYYroimvBggV69913VV5ersTExNBVktPp1IABA2RZlhYuXKglS5YoMzNTmZmZWrJkiS6++GLNnj07dJ65c+dq0KBB8nq9kqRHH31UOTk5Wrp0qaZPn67y8nJt27ZNO3fujOJTBQDEg4jC9eqrr0qScnNzw46vXLlS8+bNkyQ9/vjj+u677/TQQw/p66+/VlZWlrZu3arExMTQ+MbGRvXp8+PFXnZ2ttatW6dnnnlGixcv1tChQ7V+/XplZWWd49MCAMSr8/oeV2/C97gAwFy/2Pe4AAD4pREuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoyT09AKi5cwfcm5XmxQXf9MZAC4c7WqT9ON/y7sTN+E6ceKEJGmntvTwSgAA5+rEiRNyOp3djrHsn5M3A5w+fVpHjx5VYmKiLMvq8HggEFBaWpqOHDmipKSkHlhhfOH1jC5ez+ji9YyuX+L1tG1bJ06ckMfjUZ8+3X+KFTdXXH369NHgwYPPOi4pKYk3chTxekYXr2d08XpGV6xfz7NdaZ3BzRkAAKMQLgCAUS6YcDkcDj377LNyOBw9vZS4wOsZXbye0cXrGV297fWMm5szAAAXhgvmigsAEB8IFwDAKIQLAGAUwgUAMMoFE65XXnlFGRkZ6t+/v8aOHasPPvigp5dkpJKSElmWFba5XK6eXpYxampqNG3aNHk8HlmWpU2bNoU9btu2SkpK5PF4NGDAAOXm5urTTz/tmcUa4Gyv57x58zq8XydMmNAzi+3lvF6vxo8fr8TERKWkpGjGjBk6dOhQ2Jje8v68IMK1fv16LVy4UE8//bT279+vm266Sfn5+WpsbOzppRlpxIgRampqCm2ffPJJTy/JGC0tLRo9erRKS0s7ffyFF17QSy+9pNLSUu3du1cul0uTJ08O/S5OhDvb6ylJt912W9j7dcsWfp9pZ6qrq7VgwQLt3r1blZWVam9vV15enlpaWkJjes37074AXH/99fb8+fPDjg0bNsx+8skne2hF5nr22Wft0aNH9/Qy4oIke+PGjaH906dP2y6Xy37++edDx06dOmU7nU77tdde64EVmuWnr6dt23ZhYaE9ffr0HlmP6Zqbm21JdnV1tW3bvev9GfdXXK2trfrwww+Vl5cXdjwvL0+7du3qoVWZrb6+Xh6PRxkZGbrnnnv0+eef9/SS4kJDQ4N8Pl/Ye9XhcOjmm2/mvXoeqqqqlJKSoquvvlr333+/mpube3pJRvD7/ZKk5ORkSb3r/Rn34Tp27Ji+//57paamhh1PTU2Vz+froVWZKysrS6tXr9b777+vN954Qz6fT9nZ2Tp+/HhPL814Z96PvFejJz8/X++88462b9+uF198UXv37tUtt9yiYDDY00vr1WzbVlFRkW688UaNHDlSUu96f8bNb4c/m5/+qRPbtjv98yfoXn5+fuifR40apYkTJ2ro0KF6++23VVRU1IMrix+8V6OnoKAg9M8jR47UuHHjlJ6ers2bN2vmzJk9uLLe7eGHH9bHH3+snTt3dnisN7w/4/6K64orrtBFF13U4f8ImpubO/yfAyJ3ySWXaNSoUaqvr+/ppRjvzN2ZvFdjx+12Kz09nfdrNx555BFVVFRox44dYX8qqje9P+M+XP369dPYsWNVWVkZdryyslLZ2dk9tKr4EQwGdfDgQbnd7p5eivEyMjLkcrnC3qutra2qrq7mvRolx48f15EjR3i/dsK2bT388MPasGGDtm/froyMjLDHe9P784L4UWFRUZHmzJmjcePGaeLEiXr99dfV2Nio+fPn9/TSjLNo0SJNmzZNQ4YMUXNzs/74xz8qEAiosLCwp5dmhJMnT+rw4cOh/YaGBtXV1Sk5OVlDhgzRwoULtWTJEmVmZiozM1NLlizRxRdfrNmzZ/fgqnuv7l7P5ORklZSUaNasWXK73friiy/01FNP6YorrtAdd9zRg6vunRYsWKB3331X5eXlSkxMDF1ZOZ1ODRgwQJZl9Z735y96D2MPWrFihZ2enm7369fPvu6660K3eCIyBQUFttvttvv27Wt7PB575syZ9qefftrTyzLGjh07bEkdtsLCQtu2f7jl+Nlnn7VdLpftcDjsnJwc+5NPPunZRfdi3b2e3377rZ2Xl2cPHDjQ7tu3rz1kyBC7sLDQbmxs7Oll90qdvY6S7JUrV4bG9Jb3J3/WBABglLj/jAsAEF8IFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMMr/AzeOm9E8htvqAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# initialize environment with start values\n",
    "C_t = np.zeros((grid_size['x']+2,grid_size['y']+2))\n",
    "C_t[2,10] = 2000\n",
    "\n",
    "\n",
    "plt.imshow(C_t, vmin=0, vmax=2000)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate mean between two cells\n",
    "def alpha_interblock(alpha1, alpha2, harmonic=False):\n",
    "    if not harmonic:\n",
    "        return 0.5 * (alpha1 + alpha2)\n",
    "    else:\n",
    "        return 2 / ((1/alpha1) + (1/alpha2))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulate one time step\n",
    "def simulate(C_t): \n",
    "        C_t1 = np.zeros((grid_size['x'],grid_size['y']))\n",
    "\n",
    "        # inner cells\n",
    "        for i in range(1, grid_size['x']-1):\n",
    "                for j in range(1, grid_size['y']-1):\n",
    "                        C_t1[i,j] = C_t[i,j] \\\n",
    "                                + time_step/delta_x**2 * (alpha_interblock(alpha_x[i+1,j], alpha_x[i,j]) * C_t[i+1,j]\n",
    "                                        - (alpha_interblock(alpha_x[i+1,j], alpha_x[i,j]) + alpha_interblock(alpha_x[i-1,j], alpha_x[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_x[i-1,j], alpha_x[i,j]) * C_t[i-1,j]) \\\n",
    "                                + time_step/delta_y**2 * (alpha_interblock(alpha_y[i,j+1], alpha_y[i,j]) * C_t[i,j+1]\n",
    "                                        - (alpha_interblock(alpha_y[i,j+1], alpha_y[i,j]) + alpha_interblock(alpha_y[i,j-1], alpha_y[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_y[i,j-1], alpha_y[i,j]) * C_t[i,j-1])\n",
    "                        \n",
    "        # boundary conditions\n",
    "        # left\n",
    "        for i in range(1, grid_size['y']-1):\n",
    "               C_t1[i,0] = C_t[i,0] \\\n",
    "                        +  time_step/delta_x**2 * (alpha_interblock(alpha_x[i,1], alpha_x[i,0]) * C_t[i,1] \n",
    "                                - (alpha_interblock(alpha_x[i,1], alpha_x[i,0]) + 2 * alpha_x[i,0]) * C_t[i,0]\n",
    "                                + 2 * alpha_x[i,0] * bc_left)\n",
    "\n",
    "        # right\n",
    "        n = grid_size['x'] - 1 # maximum index in x-direction (columns)\n",
    "        for i in range(1, grid_size['y']-1):\n",
    "               C_t1[i,n] = C_t[i,n] \\\n",
    "                        +  time_step/delta_x**2 * (2 * alpha_x[i,n] * bc_right \n",
    "                                - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) + 2 * alpha_x[i,n]) * C_t[i,n]\n",
    "                                + alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * C_t[i,n-1])\n",
    "\n",
    "        # top\n",
    "        for j in range(1, grid_size['x']-1):\n",
    "               C_t1[0,j] = C_t[0,j] \\\n",
    "                        +  time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * C_t[1,j] \n",
    "                                - (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) + 2 * alpha_y[0,j]) * C_t[0,j]\n",
    "                                + 2 * alpha_y[0,j] * bc_top)\n",
    "\n",
    "        # bottom\n",
    "        m = grid_size['y'] - 1 # maximum index in y-direction (rows)\n",
    "        for j in range(1, grid_size['x']-1):\n",
    "               C_t1[m,j] = C_t[m,j] \\\n",
    "                        +  time_step/delta_y**2 * (2 * alpha_y[m,j] * bc_bottom \n",
    "                                - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) + 2 * alpha_y[m,j]) * C_t[m,j]\n",
    "                                + alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * C_t[m-1,j])\n",
    "\n",
    "\n",
    "        return C_t1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAAGiCAYAAACcWg7FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCv0lEQVR4nO3df1hUZd4/8PeoMCgLk4gwM4lEXkolxCKYgJWYBdLij7RVsy9iuZhPZkvIZVpbsc/uE9au6ZWmmYvir9JnH0VtdVFYBSLQ/EXrr1jcSNAYKS+YEVsBmfv7h8tZR4YD48wAznm/nuu+4pzzOffcc3QfPt6/jkoIIUBEREQEoFd3N4CIiIh6DiYGREREJGFiQERERBImBkRERCRhYkBEREQSJgZEREQkYWJAREREEiYGREREJGFiQERERBImBkRERCRhYkBEROREmZmZGDlyJLy8vODn54fJkyejvLzcIkYIgYyMDOj1evTt2xexsbE4c+aMRUxjYyMWLFgAX19feHp6YuLEibh48aJFTF1dHZKSkqDRaKDRaJCUlIT6+nqb2svEgIiIyIkKCwsxf/58HD58GHl5ebhx4wbi4uJw7do1Keb999/HBx98gFWrVuHo0aPQarV46qmncPXqVSkmNTUVOTk52LZtG4qLi9HQ0IDExES0tLRIMTNnzkRZWRlyc3ORm5uLsrIyJCUl2dZgQURERF2mtrZWABCFhYVCCCHMZrPQarVi6dKlUsz169eFRqMRH3/8sRBCiPr6euHm5ia2bdsmxVy6dEn06tVL5ObmCiGEOHv2rAAgDh8+LMWUlpYKAOKbb77pdPv62JcH9Rxmsxnff/89vLy8oFKpurs5RERkAyEErl69Cr1ej169nNeZff36dTQ1NTmkLiFEm983arUaarVa9j6j0QgA8PHxAQBUVlbCYDAgLi7Oop4xY8agpKQEL730Eo4fP47m5maLGL1ej5CQEJSUlCA+Ph6lpaXQaDQYNWqUFBMVFQWNRoOSkhIEBwd36nu5TGLw/fffIyAgoLubQUREdqiursagQYOcUvf169cRFPgzGGpbOg7uhJ/97GdoaGiwOPfOO+8gIyOj3XuEEEhLS8Ojjz6KkJAQAIDBYAAA+Pv7W8T6+/vjwoULUoy7uzv69+/fJqb1foPBAD8/vzaf6efnJ8V0hsskBl5eXgCAR/E0+sCtm1tDRES2uIFmFGOf9P/LnaGpqQmG2hZUHg+Et5d9vRKmq2YERVxAdXU1vL29pfMd9Ra88sor+Pvf/47i4uI2127vfbDWI3G722OsxXemnlu5TGLQ+qX7wA19VEwMiIjuKuLmf7piKNjbq5fdiYFUl7e3RWIgZ8GCBdizZw+KioosekW0Wi2Am//i1+l00vna2lqpF0Gr1aKpqQl1dXUWvQa1tbWIiYmRYi5fvtzmc3/44Yc2vRFyuCqBiIgUpUWYHVI6SwiBV155BTt37sTBgwcRFBRkcT0oKAharRZ5eXnSuaamJhQWFkq/9CMiIuDm5mYRU1NTg9OnT0sx0dHRMBqN+Oqrr6SYI0eOwGg0SjGd4TI9BkRERJ1hhoC5tYvCjjo6a/78+fj000+xe/dueHl5SeP9Go0Gffv2hUqlQmpqKt59910MHToUQ4cOxbvvvot+/fph5syZUuycOXOwcOFCDBgwAD4+PkhPT0doaCiefPJJAMCDDz6I8ePHIyUlBWvXrgUAzJ07F4mJiZ2eeAg4scdg9erVCAoKgoeHByIiIvDFF1/IxhcWFiIiIgIeHh64//778fHHHzuraUREpGBmB/1fZ61ZswZGoxGxsbHQ6XRS2b59uxSzaNEipKam4uWXX0ZkZCQuXbqEAwcOWMy5WL58OSZPnoxp06Zh9OjR6NevHz7//HP07t1bitm6dStCQ0MRFxeHuLg4PPzww9i8ebNNz0clhLAvbbJi+/btSEpKwurVqzF69GisXbsWf/rTn3D27FkMHjy4TXxlZSVCQkKQkpKCl156CV9++SVefvllfPbZZ5g6dWqnPtNkMkGj0SAWkzjHgIjoLnNDNKMAu2E0Gjs9Zm+r1t8T35cPcsjkQ33wRae2t7s4JTEYNWoURowYgTVr1kjnHnzwQUyePBmZmZlt4l9//XXs2bMH586dk87NmzcPX3/9NUpLSzv1mUwMiIjuXl2ZGFR/c69DEoOABy65ZGLg8KGEpqYmHD9+3GITBgCIi4tDSUmJ1XtKS0vbxMfHx+PYsWNobm62ek9jYyNMJpNFISIi6kjrHAN7i6tyeGLw448/oqWlxepGDe1tsGAwGKzG37hxAz/++KPVezIzM6WXRGg0Gm5uRERE5ABOm3xo60YN1uKtnW+1ZMkSGI1GqVRXV9vZYiIiUgIzBFrsLK7cY+Dw5Yq+vr7o3bt3m96BWzdquJ1Wq7Ua36dPHwwYMMDqPZ3Zj5qIiOh2Xb1c8W7j8B4Dd3d3REREWGzCAAB5eXntbrAQHR3dJv7AgQOIjIyEmxsnEhIREXUVpwwlpKWl4U9/+hPWr1+Pc+fO4bXXXkNVVRXmzZsH4OYwwKxZs6T4efPm4cKFC0hLS8O5c+ewfv16ZGVlIT093RnNIyIiBWsRwiHFVTll58Pp06fjypUr+O///m/U1NQgJCQE+/btQ2BgIICb2zhWVVVJ8UFBQdi3bx9ee+01fPTRR9Dr9fjwww87vYcBERFRZ5n/Xeytw1U5ZR+D7sB9DIiI7l5duY/BN+f84WXnPgZXr5rxwIOXXXIfA74rgYiIFKV1ZYG9dbgqJgZERKQoLeJmsbcOV8XEgIiIFIVzDOQ5bYMjIiIiuvuwx4CIiBTFDBVa0P5OvJ2tw1UxMSAiIkUxi5vF3jpcFYcSiIiISMIeAyIiUpQWBwwl2Ht/T8bEgIiIFIWJgTwOJRAREZGEPQZERKQoZqGCWdi5KsHO+3syJgZERKQoHEqQx6EEIiIikrDHgIiIFKUFvdBi57+LWxzUlp6IiQERESmKcMAcA8E5BkRERK6BcwzkcY4BERERSdhjQEREitIieqFF2DnHwIXflcDEgIiIFMUMFcx2dpib4bqZAYcSiIiISMIeAyIiUhROPpTHxICIiBTFMXMMOJRARERECsAeAyIiUpSbkw/tfIkShxKIiIhcg9kBWyJzVQIREREpAhMDIiJSlNbJh/YWWxQVFWHChAnQ6/VQqVTYtWuXxXWVSmW1/OEPf5BiYmNj21yfMWOGRT11dXVISkqCRqOBRqNBUlIS6uvrbWorEwMiIlIUM3o5pNji2rVrCAsLw6pVq6xer6mpsSjr16+HSqXC1KlTLeJSUlIs4tauXWtxfebMmSgrK0Nubi5yc3NRVlaGpKQkm9rKOQZERKQoLUKFFjvfjmjr/QkJCUhISGj3ulartTjevXs3xo4di/vvv9/ifL9+/drEtjp37hxyc3Nx+PBhjBo1CgCwbt06REdHo7y8HMHBwZ1qK3sMiIiI7pDJZLIojY2Ndtd5+fJl7N27F3PmzGlzbevWrfD19cXw4cORnp6Oq1evStdKS0uh0WikpAAAoqKioNFoUFJS0unPZ48BEREpSosDViW0/HtVQkBAgMX5d955BxkZGXbVvXHjRnh5eWHKlCkW559//nkEBQVBq9Xi9OnTWLJkCb7++mvk5eUBAAwGA/z8/NrU5+fnB4PB0OnPZ2JARESKYha9YLZz50Pzv3c+rK6uhre3t3RerVbbVS8ArF+/Hs8//zw8PDwszqekpEg/h4SEYOjQoYiMjMSJEycwYsQIADcnMd5OCGH1fHs4lEBERHSHvL29LYq9icEXX3yB8vJy/OpXv+owdsSIEXBzc0NFRQWAm/MULl++3Cbuhx9+gL+/f6fbwMSAiIgUpXUowd7iDFlZWYiIiEBYWFiHsWfOnEFzczN0Oh0AIDo6GkajEV999ZUUc+TIERiNRsTExHS6DRxKICIiRTHD9lUF1uqwRUNDA86fPy8dV1ZWoqysDD4+Phg8eDCAmxMZ//znP2PZsmVt7v/nP/+JrVu34umnn4avry/Onj2LhQsXIjw8HKNHjwYAPPjggxg/fjxSUlKkZYxz585FYmJip1ckAE7oMcjMzMTIkSPh5eUFPz8/TJ48GeXl5bL3FBQUWN3Y4ZtvvnF084iIiLrcsWPHEB4ejvDwcABAWloawsPD8fbbb0sx27ZtgxACzz33XJv73d3d8be//Q3x8fEIDg7Gq6++iri4OOTn56N3795S3NatWxEaGoq4uDjExcXh4YcfxubNm21qq0oIx747cvz48ZgxYwZGjhyJGzdu4M0338SpU6dw9uxZeHp6Wr2noKAAY8eORXl5ucUkjoEDB1p8YTkmkwkajQaxmIQ+KjeHfBciIuoaN0QzCrAbRqPR4veAI7X+nlhzYiT6/sy+DvN/NdzAf4046tT2dheHDyXk5uZaHG/YsAF+fn44fvw4Hn/8cdl7/fz8cM899zi6SURERJI72dLYWh2uyunfzGg0AgB8fHw6jA0PD4dOp8O4ceNw6NAh2djGxsY2G0sQERGRfZyaGAghkJaWhkcffRQhISHtxul0OnzyySfYsWMHdu7cieDgYIwbNw5FRUXt3pOZmSm9JEKj0bTZZIKIiMgaM1QOKa7K4XMMbjV//nzs3bsXxcXFGDRokE33TpgwASqVCnv27LF6vbGx0WLrSZPJhICAAM4xICK6C3XlHIPlx2IcMsfgtcgSzjGwxYIFC7Bnzx4UFRXZnBQAN/d33rJlS7vX1Wq1Q3aYIiIiZXHMlsiuO8fA4YmBEAILFixATk4OCgoKEBQUdEf1nDx5Utq0gYiIiLqGwxOD+fPn49NPP8Xu3bvh5eUlvbhBo9Ggb9++AIAlS5bg0qVL2LRpEwBgxYoVuO+++zB8+HA0NTVhy5Yt2LFjB3bs2OHo5hERkcKZhQpmezc4svP+nszhicGaNWsAALGxsRbnN2zYgNmzZwMAampqUFVVJV1rampCeno6Ll26hL59+2L48OHYu3cvnn76aUc3j4iIFM7sgKEEM4cSOq8zcxmzs7MtjhctWoRFixY5uilERERkI74rgYiIFMUxr11mjwEREZFLaIEKLXbuQ2Dv/T2Z66Y8REREZDP2GBARkaJwKEEeEwMiIlKUFtg/FNDimKb0SK6b8hAREZHN2GNARESKwqEEeUwMiIhIUVpEL7TY+Yvd3vt7MiYGRESkKMIBr00WXK5IRERESsAeAyIiUhQOJchjYkBERIrCtyvKc92Uh4iIiGzGHgMiIlKUFge8dtne+3syJgZERKQoHEqQ57opDxEREdmMPQZERKQoZvSC2c5/F9t7f0/GxICIiBSlRajQYudQgL3392Sum/IQERGRzdhjQEREisLJh/KYGBARkaIIB7xdUXDnQyIiItfQAhVa7HwJkr3392Sum/IQERGRzdhjQEREimIW9s8RMAsHNaYHYmJARESKYnbAHAN77+/JXPebERERkc2YGBARkaKYoXJIsUVRUREmTJgAvV4PlUqFXbt2WVyfPXs2VCqVRYmKirKIaWxsxIIFC+Dr6wtPT09MnDgRFy9etIipq6tDUlISNBoNNBoNkpKSUF9fb1NbmRgQEZGitO58aG+xxbVr1xAWFoZVq1a1GzN+/HjU1NRIZd++fRbXU1NTkZOTg23btqG4uBgNDQ1ITExES0uLFDNz5kyUlZUhNzcXubm5KCsrQ1JSkk1t5RwDIiIiJ0tISEBCQoJsjFqthlartXrNaDQiKysLmzdvxpNPPgkA2LJlCwICApCfn4/4+HicO3cOubm5OHz4MEaNGgUAWLduHaKjo1FeXo7g4OBOtZU9BkREpCitkw/tLQBgMpksSmNj4x23q6CgAH5+fhg2bBhSUlJQW1srXTt+/Diam5sRFxcnndPr9QgJCUFJSQkAoLS0FBqNRkoKACAqKgoajUaK6QwmBkREpChmqKRtke+4/HuOQUBAgDSer9FokJmZeUdtSkhIwNatW3Hw4EEsW7YMR48exRNPPCElGgaDAe7u7ujfv7/Fff7+/jAYDFKMn59fm7r9/PykmM7gUAIREdEdqq6uhre3t3SsVqvvqJ7p06dLP4eEhCAyMhKBgYHYu3cvpkyZ0u59QgioVP+Z73Drz+3FdIQ9BkREpCjCASsSxL97DLy9vS3KnSYGt9PpdAgMDERFRQUAQKvVoqmpCXV1dRZxtbW18Pf3l2IuX77cpq4ffvhBiukMJgZERKQodg8jOODtjB25cuUKqqurodPpAAARERFwc3NDXl6eFFNTU4PTp08jJiYGABAdHQ2j0YivvvpKijly5AiMRqMU0xkcSiAiIkXpjp0PGxoacP78eem4srISZWVl8PHxgY+PDzIyMjB16lTodDp89913eOONN+Dr64tnnnkGAKDRaDBnzhwsXLgQAwYMgI+PD9LT0xEaGiqtUnjwwQcxfvx4pKSkYO3atQCAuXPnIjExsdMrEgAmBkRERE537NgxjB07VjpOS0sDACQnJ2PNmjU4deoUNm3ahPr6euh0OowdOxbbt2+Hl5eXdM/y5cvRp08fTJs2Df/6178wbtw4ZGdno3fv3lLM1q1b8eqrr0qrFyZOnCi7d4I1KiGEQ18FkZGRgd/+9rcW526dNWlNYWEh0tLScObMGej1eixatAjz5s2z6XNNJhM0Gg1iMQl9VG531HYiIuoeN0QzCrAbRqPRYjKfI7X+nph04EW4ebrbVVfztSbsjlvv1PZ2F6f0GAwfPhz5+fnS8a3ZzO0qKyvx9NNPIyUlBVu2bMGXX36Jl19+GQMHDsTUqVOd0TwiIlKwO9nS2FodrsopiUGfPn3a3b3pdh9//DEGDx6MFStWALg5RnLs2DH88Y9/lE0MGhsbLTaSMJlMdrWZiIiInLQqoaKiAnq9HkFBQZgxYwa+/fbbdmNLS0stdnICgPj4eBw7dgzNzc3t3peZmWmxqURAQIDD2k9ERK7rbliV0J0cnhiMGjUKmzZtwv79+7Fu3ToYDAbExMTgypUrVuMNBkOb9ZX+/v64ceMGfvzxx3Y/Z8mSJTAajVKprq526PcgIiLXxMRAnsOHEm59SURoaCiio6MxZMgQbNy4UZqFebvbd2RqnQ8pt1OTWq122EYSREREdJPTlyt6enoiNDRU2r3pdlqtts2KhdraWvTp0wcDBgxwdvOIiEhhHPEvflfuMXD6zoeNjY04d+6ctHvT7aKjoy12cgKAAwcOIDIyEm5uXHZIRESOxaEEeQ5PDNLT01FYWIjKykocOXIEzz77LEwmE5KTkwHcnBswa9YsKX7evHm4cOEC0tLScO7cOaxfvx5ZWVlIT093dNOIiIioAw4fSrh48SKee+45/Pjjjxg4cCCioqJw+PBhBAYGAri5t3NVVZUUHxQUhH379uG1117DRx99BL1ejw8//JB7GBARkVMI2L8PgUN3BuxhHJ4YbNu2TfZ6dnZ2m3NjxozBiRMnHN0UIiKiNjjHQB7flUBERIrCxEAeX7tMREREEvYYEBGRorDHQB4TAyIiUhQmBvI4lEBEREQS9hgQEZGiCKGCsPNf/Pbe35MxMSAiIkUxQ2X3Pgb23t+TcSiBiIiIJOwxICIiReHkQ3lMDIiISFE4x0AehxKIiIhIwh4DIiJSFA4lyGNiQEREisKhBHlMDIiISFGEA3oMXDkx4BwDIiIikrDHgIiIFEUAEML+OlwVEwMiIlIUM1RQcefDdnEogYiIiCTsMSAiIkXhqgR5TAyIiEhRzEIFFfcxaBeHEoiIiEjCHgMiIlIUIRywKsGFlyUwMSAiIkXhHAN5HEogIiIiCXsMiIhIUdhjII+JARERKQpXJcjjUAIRESlK6+RDe4stioqKMGHCBOj1eqhUKuzatUu61tzcjNdffx2hoaHw9PSEXq/HrFmz8P3331vUERsbC5VKZVFmzJhhEVNXV4ekpCRoNBpoNBokJSWhvr7eprYyMSAiInKya9euISwsDKtWrWpz7aeffsKJEyfw1ltv4cSJE9i5cyf+8Y9/YOLEiW1iU1JSUFNTI5W1a9daXJ85cybKysqQm5uL3NxclJWVISkpyaa2ciiBiIgU5ea/+O2dY3DzvyaTyeK8Wq2GWq1uE5+QkICEhASrdWk0GuTl5VmcW7lyJR555BFUVVVh8ODB0vl+/fpBq9VarefcuXPIzc3F4cOHMWrUKADAunXrEB0djfLycgQHB3fqu7HHgIiIFKV18qG9BQACAgKkbnuNRoPMzEyHtNFoNEKlUuGee+6xOL9161b4+vpi+PDhSE9Px9WrV6VrpaWl0Gg0UlIAAFFRUdBoNCgpKen0Z7PHgIiI6A5VV1fD29tbOrbWW2Cr69evY/HixZg5c6ZF3c8//zyCgoKg1Wpx+vRpLFmyBF9//bXU22AwGODn59emPj8/PxgMhk5/PhMDIiJSFPHvYm8dAODt7W3xy9tezc3NmDFjBsxmM1avXm1xLSUlRfo5JCQEQ4cORWRkJE6cOIERI0YAAFSqtkMkQgir59vDoQQiIlIURw4lOFJzczOmTZuGyspK5OXldZhwjBgxAm5ubqioqAAAaLVaXL58uU3cDz/8AH9//063g4kBERFRN2tNCioqKpCfn48BAwZ0eM+ZM2fQ3NwMnU4HAIiOjobRaMRXX30lxRw5cgRGoxExMTGdbguHEoiISFkcOZbQSQ0NDTh//rx0XFlZibKyMvj4+ECv1+PZZ5/FiRMn8Je//AUtLS3SnAAfHx+4u7vjn//8J7Zu3Yqnn34avr6+OHv2LBYuXIjw8HCMHj0aAPDggw9i/PjxSElJkZYxzp07F4mJiZ1ekQA4ocfgvvvua7MBg0qlwvz5863GFxQUWI3/5ptvHN00IiIiwBHDCDYOJRw7dgzh4eEIDw8HAKSlpSE8PBxvv/02Ll68iD179uDixYv4+c9/Dp1OJ5XW1QTu7u7429/+hvj4eAQHB+PVV19FXFwc8vPz0bt3b+lztm7ditDQUMTFxSEuLg4PP/wwNm/ebFNbHd5jcPToUbS0tEjHp0+fxlNPPYVf/vKXsveVl5dbjKcMHDjQ0U0jIiLqltcux8bGQsjcJHcNuLkssrCwsMPP8fHxwZYtW2xr3G0cnhjc/gt96dKlGDJkCMaMGSN7n5+fX5v1mkRERNS1nDr5sKmpCVu2bMGLL77Y4VKJ8PBw6HQ6jBs3DocOHeqw7sbGRphMJotCRETUkZ66KqGncGpisGvXLtTX12P27Nntxuh0OnzyySfYsWMHdu7cieDgYIwbNw5FRUWydWdmZlrsNhUQEODg1hMRkUtqnSNgb3FRKtHRwIYd4uPj4e7ujs8//9ym+yZMmACVSoU9e/a0G9PY2IjGxkbp2GQyISAgALGYhD4qtztuMxERdb0bohkF2A2j0ejQDYNuZTKZoNFocF/WW+jVz8Ouusw/Xcd3c37n1PZ2F6ctV7xw4QLy8/Oxc+dOm++NiorqcPJEey+qICIiktMdkw/vJk5LDDZs2AA/Pz/84he/sPnekydPShs2EBEROVQ37GNwN3FKYmA2m7FhwwYkJyejTx/Lj1iyZAkuXbqETZs2AQBWrFiB++67D8OHD5cmK+7YsQM7duxwRtOIiIhIhlMSg/z8fFRVVeHFF19sc62mpgZVVVXScVNTE9LT03Hp0iX07dsXw4cPx969e/H00087o2lERKRwjlhV4MqrEpySGMTFxbW7WUN2drbF8aJFi7Bo0SJnNIOIiMg6Fx4KsBdfokREREQSvkSJiIgUhUMJ8pgYEBGRsnBVgiwmBkREpDCqfxd763BNnGNAREREEvYYEBGRsnAoQRYTAyIiUhYmBrI4lEBEREQS9hgQEZGyOOK1yVyuSERE5Br4dkV5HEogIiIiCXsMiIhIWTj5UBYTAyIiUhbOMZDFoQQiIiKSsMeAiIgURSVuFnvrcFVMDIiISFk4x0AWEwMiIlIWzjGQxTkGREREJGGPARERKQuHEmQxMSAiImVhYiCLQwlEREQkYY8BEREpC3sMZDExICIiZeGqBFkcSiAiIiIJewyIiEhRuPOhPCYGRESkLJxjIItDCURERE5WVFSECRMmQK/XQ6VSYdeuXRbXhRDIyMiAXq9H3759ERsbizNnzljENDY2YsGCBfD19YWnpycmTpyIixcvWsTU1dUhKSkJGo0GGo0GSUlJqK+vt6mtTAyIiIic7Nq1awgLC8OqVausXn///ffxwQcfYNWqVTh69Ci0Wi2eeuopXL16VYpJTU1FTk4Otm3bhuLiYjQ0NCAxMREtLS1SzMyZM1FWVobc3Fzk5uairKwMSUlJNrWVQwlERKQoKjhgjsG//2symSzOq9VqqNXqNvEJCQlISEiwWpcQAitWrMCbb76JKVOmAAA2btwIf39/fPrpp3jppZdgNBqRlZWFzZs348knnwQAbNmyBQEBAcjPz0d8fDzOnTuH3NxcHD58GKNGjQIArFu3DtHR0SgvL0dwcHCnvht7DIiISFlalyvaWwAEBARI3fYajQaZmZk2N6eyshIGgwFxcXHSObVajTFjxqCkpAQAcPz4cTQ3N1vE6PV6hISESDGlpaXQaDRSUgAAUVFR0Gg0UkxnsMeAiIjoDlVXV8Pb21s6ttZb0BGDwQAA8Pf3tzjv7++PCxcuSDHu7u7o379/m5jW+w0GA/z8/NrU7+fnJ8V0BhMDIiJSFgeuSvD29rZIDOyhUllumiSEaHOuTTNui7EW35l6bsWhBCIiUhbhoOIgWq0WANr8q762tlbqRdBqtWhqakJdXZ1szOXLl9vU/8MPP7TpjZDDxICIiKgbBQUFQavVIi8vTzrX1NSEwsJCxMTEAAAiIiLg5uZmEVNTU4PTp09LMdHR0TAajfjqq6+kmCNHjsBoNEoxncGhBCIiUpTu2PmwoaEB58+fl44rKytRVlYGHx8fDB48GKmpqXj33XcxdOhQDB06FO+++y769euHmTNnAgA0Gg3mzJmDhQsXYsCAAfDx8UF6ejpCQ0OlVQoPPvggxo8fj5SUFKxduxYAMHfuXCQmJnZ6RQJwBz0GjtikwZodO3bgoYceglqtxkMPPYScnBxbm0ZERNSxbhhKOHbsGMLDwxEeHg4ASEtLQ3h4ON5++20AwKJFi5CamoqXX34ZkZGRuHTpEg4cOAAvLy+pjuXLl2Py5MmYNm0aRo8ejX79+uHzzz9H7969pZitW7ciNDQUcXFxiIuLw8MPP4zNmzfb1FaVEMKmr/fXv/4VX375JUaMGIGpU6ciJycHkydPlq6/9957+J//+R9kZ2dj2LBh+P3vf4+ioiKUl5dbfMFblZaW4rHHHsPvfvc7PPPMM8jJycHbb7+N4uJii2UXckwmEzQaDWIxCX1UbrZ8JSIi6mY3RDMKsBtGo9Fhk/lu1/p74r7f/w96eXjYVZf5+nV895s3ndre7mJzYmBxs0plkRgIIaDX65GamorXX38dwM0tHP39/fHee+/hpZdeslrP9OnTYTKZ8Ne//lU6N378ePTv3x+fffZZp9rCxICI6O7VpYnB7xyUGLzlmomBQycfdmaTBmtKS0st7gGA+Ph42XsaGxthMpksChERUUda5xjYW1yVQxMDuU0a5DZXMBgMNt+TmZlpsdtUQECAHS0nIiIiwEnLFe9kkwZb71myZAmMRqNUqqur77zBRESkHA7cEtkVOXS54q2bNOh0Oun8rRswtHef3MYO1rT3ogoiIiJZDtz50BU5tMegM5s0WBMdHW1xDwAcOHDApg0ZiIiIOoNzDOTZ3GNg7yYNADBr1izce++90luofv3rX+Pxxx/He++9h0mTJmH37t3Iz89HcXGxA74iERERdZbNicGxY8cwduxY6TgtLQ0AkJycjOzsbCxatAj/+te/8PLLL6Ourg6jRo1qs0lDVVUVevX6T2dFTEwMtm3bht/85jd46623MGTIEGzfvr3TexgQERF1GocSZNm1j0FPwn0MiIjuXl25j8H9b72L3nbuY9By/Tq+/d0b3MeAiIiIXBtfokRERMrCoQRZTAyIiEhZmBjI4lACERERSdhjQEREiuKIfQhceR8D9hgQERGRhIkBERERSTiUQEREysLJh7KYGBARkaJwjoE8JgZERKQ8LvyL3V6cY0BEREQS9hgQEZGycI6BLCYGRESkKJxjII9DCURERCRhjwERESkLhxJkMTEgIiJF4VCCPA4lEBERkYQ9BkREpCwcSpDFxICIiJSFiYEsDiUQERGRhD0GRESkKJx8KI+JARERKQuHEmQxMSAiImVhYiCLcwyIiIhIwh4DIiJSFM4xkMfEgIiIlIVDCbI4lEBERORE9913H1QqVZsyf/58AMDs2bPbXIuKirKoo7GxEQsWLICvry88PT0xceJEXLx40SntZWJARESK0jqUYG/prKNHj6KmpkYqeXl5AIBf/vKXUsz48eMtYvbt22dRR2pqKnJycrBt2zYUFxejoaEBiYmJaGlpccgzuRWHEoiISFm6eChh4MCBFsdLly7FkCFDMGbMGOmcWq2GVqu1er/RaERWVhY2b96MJ598EgCwZcsWBAQEID8/H/Hx8ba3XwZ7DIiIiO6QyWSyKI2NjbLxTU1N2LJlC1588UWoVCrpfEFBAfz8/DBs2DCkpKSgtrZWunb8+HE0NzcjLi5OOqfX6xESEoKSkhKHfycmBkREpCzCQQVAQEAANBqNVDIzM2U/eteuXaivr8fs2bOlcwkJCdi6dSsOHjyIZcuW4ejRo3jiiSekJMNgMMDd3R39+/e3qMvf3x8Gg8GeJ2EVhxKIiEhRVP8u9tYBANXV1fD29pbOq9Vq2fuysrKQkJAAvV4vnZs+fbr0c0hICCIjIxEYGIi9e/diypQp7dYlhLDodXAUJgZERER3yNvb2yIxkHPhwgXk5+dj586dsnE6nQ6BgYGoqKgAAGi1WjQ1NaGurs6i16C2thYxMTF33vh2cCiBiIiUxYFDCbbYsGED/Pz88Itf/EI27sqVK6iuroZOpwMAREREwM3NTVrNAAA1NTU4ffq0UxID9hgQEZGidMfOh2azGRs2bEBycjL69PnPr96GhgZkZGRg6tSp0Ol0+O677/DGG2/A19cXzzzzDABAo9Fgzpw5WLhwIQYMGAAfHx+kp6cjNDRUWqXgSDb3GBQVFWHChAnQ6/VQqVTYtWuXdK25uRmvv/46QkND4enpCb1ej1mzZuH777+XrTM7O9vq5g/Xr1+3+QsRERHJ6oYeg/z8fFRVVeHFF1+0ON+7d2+cOnUKkyZNwrBhw5CcnIxhw4ahtLQUXl5eUtzy5csxefJkTJs2DaNHj0a/fv3w+eefo3fv3nfwAOTZ3GNw7do1hIWF4YUXXsDUqVMtrv300084ceIE3nrrLYSFhaGurg6pqamYOHEijh07Jluvt7c3ysvLLc55eHjY2jwiIqIeJy4uDkK0zSb69u2L/fv3d3i/h4cHVq5ciZUrVzqjeRZsTgwSEhKQkJBg9ZpGo7EYAwGAlStX4pFHHkFVVRUGDx7cbr0qlardzR2IiIgcyoXfdWAvp08+NBqNUKlUuOeee2TjGhoaEBgYiEGDBiExMREnT56UjW9sbGyzsQQREVFHunpL5LuNUxOD69evY/HixZg5c6bsco4HHngA2dnZ2LNnDz777DN4eHhg9OjR0lINazIzMy02lQgICHDGVyAiIlIUpyUGzc3NmDFjBsxmM1avXi0bGxUVhf/3//4fwsLC8Nhjj+F///d/MWzYMNmxlCVLlsBoNEqlurra0V+BiIhcUTctV7xbOGW5YnNzM6ZNm4bKykocPHiw05s/tOrVqxdGjhwp22OgVqs73GGKiIjodt2xXPFu4vAeg9akoKKiAvn5+RgwYIDNdQghUFZWJm3uQERERF3D5h6DhoYGnD9/XjqurKxEWVkZfHx8oNfr8eyzz+LEiRP4y1/+gpaWFukFDz4+PnB3dwcAzJo1C/fee6/0sonf/va3iIqKwtChQ2EymfDhhx+irKwMH330kSO+IxER0X908WuX7zY2JwbHjh3D2LFjpeO0tDQAQHJyMjIyMrBnzx4AwM9//nOL+w4dOoTY2FgAQFVVFXr1+k9nRX19PebOnQuDwQCNRoPw8HAUFRXhkUcesbV5REREsjiUIM/mxCA2NtbqJg2t5K61KigosDhevnw5li9fbmtTiIiIyMH4rgQiIlIWDiXIYmJARETKwsRAFhMDIiJSFM4xkOf0LZGJiIjo7sEeAyIiUhYOJchiYkBERIqiEgKqTqyg66gOV8WhBCIiIpKwx4CIiJSFQwmymBgQEZGicFWCPA4lEBERkYQ9BkREpCwcSpDFxICIiBSFQwnyOJRAREREEvYYEBGRsnAoQRYTAyIiUhQOJchjYkBERMrCHgNZnGNAREREEvYYEBGR4rjyUIC9mBgQEZGyCHGz2FuHi+JQAhEREUnYY0BERIrCVQnymBgQEZGycFWCLA4lEBERkYQ9BkREpCgq881ibx2uiokBEREpC4cSZHEogYiIiCRMDIiISFFaVyXYWzorIyMDKpXKomi1Wum6EAIZGRnQ6/Xo27cvYmNjcebMGYs6GhsbsWDBAvj6+sLT0xMTJ07ExYsXHfVILDAxICIiZWnd4MjeYoPhw4ejpqZGKqdOnZKuvf/++/jggw+watUqHD16FFqtFk899RSuXr0qxaSmpiInJwfbtm1DcXExGhoakJiYiJaWFoc9llacY0BERIrSHfsY9OnTx6KXoJUQAitWrMCbb76JKVOmAAA2btwIf39/fPrpp3jppZdgNBqRlZWFzZs348knnwQAbNmyBQEBAcjPz0d8fLx9X+Y27DEgIiK6QyaTyaI0NjZajauoqIBer0dQUBBmzJiBb7/9FgBQWVkJg8GAuLg4KVatVmPMmDEoKSkBABw/fhzNzc0WMXq9HiEhIVKMIzExICIiZREOKgACAgKg0WikkpmZ2ebjRo0ahU2bNmH//v1Yt24dDAYDYmJicOXKFRgMBgCAv7+/xT3+/v7SNYPBAHd3d/Tv37/dGEfiUAIRESmKI4cSqqur4e3tLZ1Xq9VtYhMSEqSfQ0NDER0djSFDhmDjxo2Iioq6WZ9KZXGPEKLNudt1JuZOsMeAiIjoDnl7e1sUa4nB7Tw9PREaGoqKigpp3sHt//Kvra2VehG0Wi2amppQV1fXbowjMTEgIiJl6YZVCbdqbGzEuXPnoNPpEBQUBK1Wi7y8POl6U1MTCgsLERMTAwCIiIiAm5ubRUxNTQ1Onz4txTgShxKIiEhRunpVQnp6OiZMmIDBgwejtrYWv//972EymZCcnAyVSoXU1FS8++67GDp0KIYOHYp3330X/fr1w8yZMwEAGo0Gc+bMwcKFCzFgwAD4+PggPT0doaGh0ioFR7K5x6CoqAgTJkyAXq+HSqXCrl27LK7Pnj27zUYOrWMocnbs2IGHHnoIarUaDz30EHJycmxtGhERUY9z8eJFPPfccwgODsaUKVPg7u6Ow4cPIzAwEACwaNEipKam4uWXX0ZkZCQuXbqEAwcOwMvLS6pj+fLlmDx5MqZNm4bRo0ejX79++Pzzz9G7d2+Ht9fmHoNr164hLCwML7zwAqZOnWo1Zvz48diwYYN07O7uLltnaWkppk+fjt/97nd45plnkJOTg2nTpqG4uBijRo2ytYlERETt6+J3JWzbtk32ukqlQkZGBjIyMtqN8fDwwMqVK7Fy5crOf/AdsjkxSEhIsJhhaY1arba6kUN7VqxYgaeeegpLliwBACxZsgSFhYVYsWIFPvvsM1ubSERE1K7u2ODobuKUyYcFBQXw8/PDsGHDkJKSgtraWtn40tJSi40bACA+Pl5244bGxsY2G0sQERGRfRyeGCQkJGDr1q04ePAgli1bhqNHj+KJJ55odzco4OYyDbnNHazJzMy02FQiICDAYd+BiIhcmFk4prgoh69KmD59uvRzSEgIIiMjERgYiL1790r7QFtj6+YOS5YsQVpamnRsMpmYHBARUce6eI7B3cbpyxV1Oh0CAwNRUVHRboxWq5Xd3MEatVrdqY0kiIiIbqWCA+YYOKQlPZPTNzi6cuUKqqurodPp2o2Jjo622LgBAA4cOOCUjRuIiIiofTb3GDQ0NOD8+fPScWVlJcrKyuDj4wMfHx9kZGRg6tSp0Ol0+O677/DGG2/A19cXzzzzjHTPrFmzcO+990ovm/j1r3+Nxx9/HO+99x4mTZqE3bt3Iz8/H8XFxQ74ikRERLewc+dCqQ4XZXNicOzYMYwdO1Y6bh3nT05Oxpo1a3Dq1Cls2rQJ9fX10Ol0GDt2LLZv326xUUNVVRV69fpPZ0VMTAy2bduG3/zmN3jrrbcwZMgQbN++nXsYEBGRw3G5ojybE4PY2FgImUxp//79HdZRUFDQ5tyzzz6LZ5991tbmEBERkQPxXQlERKQsXJUgi4kBEREpikoIqOycI2Dv/T0ZX7tMREREEvYYEBGRspj/Xeytw0UxMSAiIkXhUII8DiUQERGRhD0GRESkLFyVIIuJARERKQt3PpTFxICIiBSFOx/K4xwDIiIikrDHgIiIlIVDCbKYGBARkaKozDeLvXW4Kg4lEBERkYQ9BkREpCwcSpDFxICIiJSF+xjI4lACERERSdhjQEREisJ3JchjYkBERMrCOQayOJRAREREEvYYEBGRsggA9u5D4LodBkwMiIhIWTjHQB4TAyIiUhYBB8wxcEhLeiTOMSAiIiIJewyIiEhZuCpBFhMDIiJSFjMAlQPqcFEcSiAiIiIJEwMiIlKU1lUJ9pbOyszMxMiRI+Hl5QU/Pz9MnjwZ5eXlFjGzZ8+GSqWyKFFRURYxjY2NWLBgAXx9feHp6YmJEyfi4sWLDnkmt2JiQEREytI6x8De0kmFhYWYP38+Dh8+jLy8PNy4cQNxcXG4du2aRdz48eNRU1MjlX379llcT01NRU5ODrZt24bi4mI0NDQgMTERLS0tDnksrTjHgIiIyIlyc3Mtjjds2AA/Pz8cP34cjz/+uHRerVZDq9VarcNoNCIrKwubN2/Gk08+CQDYsmULAgICkJ+fj/j4eIe1lz0GRESkLA7sMTCZTBalsbGxw483Go0AAB8fH4vzBQUF8PPzw7Bhw5CSkoLa2lrp2vHjx9Hc3Iy4uDjpnF6vR0hICEpKShzxVCRMDIiISFkcmBgEBARAo9FIJTMzs4OPFkhLS8Ojjz6KkJAQ6XxCQgK2bt2KgwcPYtmyZTh69CieeOIJKdEwGAxwd3dH//79Lerz9/eHwWBw6OPhUAIREdEdqq6uhre3t3SsVqtl41955RX8/e9/R3FxscX56dOnSz+HhIQgMjISgYGB2Lt3L6ZMmdJufUIIqFT2rr20xB4DIiJSFrODCgBvb2+LIpcYLFiwAHv27MGhQ4cwaNAg2SbqdDoEBgaioqICAKDVatHU1IS6ujqLuNraWvj7+9v09TvCxICIiBSlq5crCiHwyiuvYOfOnTh48CCCgoI6vOfKlSuorq6GTqcDAERERMDNzQ15eXlSTE1NDU6fPo2YmBjbH4IMDiUQEZGydPGWyPPnz8enn36K3bt3w8vLS5oToNFo0LdvXzQ0NCAjIwNTp06FTqfDd999hzfeeAO+vr545plnpNg5c+Zg4cKFGDBgAHx8fJCeno7Q0FBplYKjMDEgIiJyojVr1gAAYmNjLc5v2LABs2fPRu/evXHq1Cls2rQJ9fX10Ol0GDt2LLZv3w4vLy8pfvny5ejTpw+mTZuGf/3rXxg3bhyys7PRu3dvh7bX5qGEoqIiTJgwAXq9HiqVCrt27bK4fvvOTa3lD3/4Q7t1ZmdnW73n+vXrNn8hIiIiWWbhmNJJQgirZfbs2QCAvn37Yv/+/aitrUVTUxMuXLiA7OxsBAQEWNTj4eGBlStX4sqVK/jpp5/w+eeft4lxBJt7DK5du4awsDC88MILmDp1apvrNTU1Fsd//etfMWfOHKuxt/L29m6zRaSHh4etzSMiIpLHtyvKsjkxSEhIQEJCQrvXb9+1affu3Rg7dizuv/9+2XpVKlW7Oz5Z09jYaLGRhMlk6vS9REREZJ1TVyVcvnwZe/fuxZw5czqMbWhoQGBgIAYNGoTExEScPHlSNj4zM9NiUwlndKcQEZErcsTmRq7bY+DUxGDjxo3w8vKS3ZwBAB544AFkZ2djz549+Oyzz+Dh4YHRo0dL6zetWbJkCYxGo1Sqq6sd3XwiInJFXfwSpbuNU1clrF+/Hs8//3yHcwWioqIsXi85evRojBgxAitXrsSHH35o9R61Wt3hDlNERERkG6clBl988QXKy8uxfft2m+/t1asXRo4cKdtjQEREdEfMDhgKsGFVwt3GaUMJWVlZiIiIQFhYmM33CiFQVlYm7fhERETkMMLsmOKibO4xaGhowPnz56XjyspKlJWVwcfHB4MHDwZwc4XAn//8ZyxbtsxqHbNmzcK9994rvYXqt7/9LaKiojB06FCYTCZ8+OGHKCsrw0cffXQn34mIiIjukM2JwbFjxzB27FjpOC0tDQCQnJyM7OxsAMC2bdsghMBzzz1ntY6qqir06vWfzor6+nrMnTsXBoMBGo0G4eHhKCoqwiOPPGJr84iIiORxHwNZKiFc49uZTCZoNBrEYhL6qNy6uzlERGSDG6IZBdgNo9Fo8RpjR2r9PfHkvfPQp5d9k9dvmBuRf+ljp7a3u/BdCUREpCzsMZDF1y4TERGRhD0GRESkLAIO6DFwSEt6JCYGRESkLBxKkMWhBCIiIpKwx4CIiJTFbAZg5wZFZm5wRERE5Bo4lCCLQwlEREQkYY8BEREpC3sMZDExICIiZeHbFWVxKIGIiIgk7DEgIiJFEcIMYedrk+29vydjYkBERMoihP1DAZxjQERE5CKEA+YYuHBiwDkGREREJGGPARERKYvZDKjsnCPAOQZEREQugkMJsjiUQERERBL2GBARkaIIsxnCzqEELlckIiJyFRxKkMWhBCIiIpKwx4CIiJTFLAAVewzaw8SAiIiURQgA9i5XdN3EgEMJREREJGGPARERKYowCwg7hxKEC/cYMDEgIiJlEWbYP5TgussVOZRARESKIszCIcVWq1evRlBQEDw8PBAREYEvvvjCCd/OfkwMiIiInGz79u1ITU3Fm2++iZMnT+Kxxx5DQkICqqqqurtpbbjMUELreM8NNNu9bwUREXWtG2gG0DVj9zdEo91DAa3tNZlMFufVajXUanWb+A8++ABz5szBr371KwDAihUrsH//fqxZswaZmZl2tcXRXCYxuHr1KgCgGPu6uSVERHSnrl69Co1G45S63d3dodVqUWxwzO+Jn/3sZwgICLA498477yAjI8PiXFNTE44fP47FixdbnI+Li0NJSYlD2uJILpMY6PV6VFdXw8vLCyqVymqMyWRCQEAAqqur4e3t3cUtvHNsd9e7W9vOdncttttxhBC4evUq9Hq90z7Dw8MDlZWVaGpqckh9Qog2v2+s9Rb8+OOPaGlpgb+/v8V5f39/GAwGh7TFkVwmMejVqxcGDRrUqVhvb+8e8z8GW7DdXe9ubTvb3bXYbsdwVk/BrTw8PODh4eH0z7Hm9iTCWmLRE3DyIRERkRP5+vqid+/ebXoHamtr2/Qi9ARMDIiIiJzI3d0dERERyMvLszifl5eHmJiYbmpV+1xmKKEz1Go13nnnHatjQD0Z29317ta2s91di+2mzkpLS0NSUhIiIyMRHR2NTz75BFVVVZg3b153N60NlXDlfR2JiIh6iNWrV+P9999HTU0NQkJCsHz5cjz++OPd3aw2mBgQERGRhHMMiIiISMLEgIiIiCRMDIiIiEjCxICIiIgkLpcY2Ppay8LCQkRERMDDwwP3338/Pv744y5q6U2ZmZkYOXIkvLy84Ofnh8mTJ6O8vFz2noKCAqhUqjblm2++6aJWAxkZGW0+X6vVyt7T3c+61X333Wf1+c2fP99qfHc976KiIkyYMAF6vR4qlQq7du2yuC6EQEZGBvR6Pfr27YvY2FicOXOmw3p37NiBhx56CGq1Gg899BBycnK6rN3Nzc14/fXXERoaCk9PT+j1esyaNQvff/+9bJ3Z2dlW/wyuX7/eJe0GgNmzZ7f5/KioqA7r7c7nDcDqc1OpVPjDH/7Qbp1d8byp53KpxMDW11pWVlbi6aefxmOPPYaTJ0/ijTfewKuvvoodO3Z0WZsLCwsxf/58HD58GHl5ebhx4wbi4uJw7dq1Du8tLy9HTU2NVIYOHdoFLf6P4cOHW3z+qVOn2o3tCc+61dGjRy3a3brpyC9/+UvZ+7r6eV+7dg1hYWFYtWqV1evvv/8+PvjgA6xatQpHjx6FVqvFU089Jb1QzJrS0lJMnz4dSUlJ+Prrr5GUlIRp06bhyJEjXdLun376CSdOnMBbb72FEydOYOfOnfjHP/6BiRMndlivt7e3xfOvqalx6Na2HT1vABg/frzF5+/bJ/8ynu5+3gDaPLP169dDpVJh6tSpsvU6+3lTDyZcyCOPPCLmzZtnce6BBx4Qixcvthq/aNEi8cADD1ice+mll0RUVJTT2tiR2tpaAUAUFha2G3Po0CEBQNTV1XVdw27zzjvviLCwsE7H98Rn3erXv/61GDJkiDCbzVav94TnDUDk5ORIx2azWWi1WrF06VLp3PXr14VGoxEff/xxu/VMmzZNjB8/3uJcfHy8mDFjhsPbLETbdlvz1VdfCQDiwoUL7cZs2LBBaDQaxzZOhrV2Jycni0mTJtlUT0983pMmTRJPPPGEbExXP2/qWVymx6D1tZZxcXEW5+Vea1laWtomPj4+HseOHUNzc7PT2irHaDQCAHx8fDqMDQ8Ph06nw7hx43Do0CFnN62NiooK6PV6BAUFYcaMGfj222/bje2Jzxq4+fdmy5YtePHFFzt8mUl3P+9bVVZWwmAwWDxTtVqNMWPGyL7Gtb0/h+589avRaIRKpcI999wjG9fQ0IDAwEAMGjQIiYmJOHnyZNc08BYFBQXw8/PDsGHDkJKSgtraWtn4nva8L1++jL1792LOnDkdxvaE503dw2USgzt5raXBYLAaf+PGDfz4449Oa2t7hBBIS0vDo48+ipCQkHbjdDodPvnkE+zYsQM7d+5EcHAwxo0bh6Kioi5r66hRo7Bp0ybs378f69atg8FgQExMDK5cuWI1vqc961a7du1CfX09Zs+e3W5MT3jet2v9O23ra1zb+3Porle/Xr9+HYsXL8bMmTNl3/L3wAMPIDs7G3v27MFnn30GDw8PjB49GhUVFV3W1oSEBGzduhUHDx7EsmXLcPToUTzxxBNobGxs956e9rw3btwILy8vTJkyRTauJzxv6j4u964EW19raS3e2vmu8Morr+Dvf/87iouLZeOCg4MRHBwsHUdHR6O6uhp//OMfu2x7zYSEBOnn0NBQREdHY8iQIdi4cSPS0tKs3tOTnnWrrKwsJCQkyL4Dvic87/bcyWtce8qrX5ubmzFjxgyYzWasXr1aNjYqKspiot/o0aMxYsQIrFy5Eh9++KGzmwoAmD59uvRzSEgIIiMjERgYiL1798r+ou0pzxsA1q9fj+eff77DuQI94XlT93GZHoM7ea2lVqu1Gt+nTx8MGDDAaW21ZsGCBdizZw8OHTqEQYMG2Xx/VFRUt2bznp6eCA0NbbcNPelZt7pw4QLy8/Pxq1/9yuZ7u/t5t64AsfU1ru39OXT1q1+bm5sxbdo0VFZWIi8vT7a3wJpevXph5MiR3fpnoNPpEBgYKNuGnvK8AeCLL75AeXn5Hf197wnPm7qOyyQGd/Jay+jo6DbxBw4cQGRkJNzc3JzW1lsJIfDKK69g586dOHjwIIKCgu6onpMnT0Kn0zm4dZ3X2NiIc+fOtduGnvCsb7dhwwb4+fnhF7/4hc33dvfzDgoKglartXimTU1NKCwslH2Na3t/Dl356tfWpKCiogL5+fl3lBgKIVBWVtatfwZXrlxBdXW1bBt6wvNulZWVhYiICISFhdl8b0943tSFumvWozNs27ZNuLm5iaysLHH27FmRmpoqPD09xXfffSeEEGLx4sUiKSlJiv/2229Fv379xGuvvSbOnj0rsrKyhJubm/i///u/Lmvzf/3XfwmNRiMKCgpETU2NVH766Scp5vZ2L1++XOTk5Ih//OMf4vTp02Lx4sUCgNixY0eXtXvhwoWioKBAfPvtt+Lw4cMiMTFReHl59ehnfauWlhYxePBg8frrr7e51lOe99WrV8XJkyfFyZMnBQDxwQcfiJMnT0qz95cuXSo0Go3YuXOnOHXqlHjuueeETqcTJpNJqiMpKcliVc6XX34pevfuLZYuXSrOnTsnli5dKvr06SMOHz7cJe1ubm4WEydOFIMGDRJlZWUWf+cbGxvbbXdGRobIzc0V//znP8XJkyfFCy+8IPr06SOOHDnSJe2+evWqWLhwoSgpKRGVlZXi0KFDIjo6Wtx77709+nm3MhqNol+/fmLNmjVW6+iO5009l0slBkII8dFHH4nAwEDh7u4uRowYYbHsLzk5WYwZM8YivqCgQISHhwt3d3dx3333tfs/HGcBYLVs2LCh3Xa/9957YsiQIcLDw0P0799fPProo2Lv3r1d2u7p06cLnU4n3NzchF6vF1OmTBFnzpxpt81CdP+zvtX+/fsFAFFeXt7mWk953q3LJG8vycnJQoibSxbfeecdodVqhVqtFo8//rg4deqURR1jxoyR4lv9+c9/FsHBwcLNzU088MADDk9w5NpdWVnZ7t/5Q4cOtdvu1NRUMXjwYOHu7i4GDhwo4uLiRElJSZe1+6effhJxcXFi4MCBws3NTQwePFgkJyeLqqoqizp62vNutXbtWtG3b19RX19vtY7ueN7Uc/G1y0RERCRxmTkGREREZD8mBkRERCRhYkBEREQSJgZEREQkYWJAREREEiYGREREJGFiQERERBImBkRERCRhYkBEREQSJgZEREQkYWJAREREkv8PvJqzF8cRDSIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "records = []\n",
    "\n",
    "for i in range(1000):\n",
    "    C_t = simulate(C_t)\n",
    "    records.append(C_t)\n",
    "\n",
    "# with np.printoptions(precision=2,floatmode='fixed',linewidth=200):\n",
    "#     print(C_t)\n",
    "\n",
    "plt.imshow(C_t, vmin=0, vmax=2000)\n",
    "plt.colorbar()\n",
    "plt.show()\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "802dbc9c1b5843abaedfa3b0fd4ebbb8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, description='w', max=999), Output()), _dom_classes=('widget-interact'…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.update(w=1)>"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def update(w = 1):\n",
    "    fig = plt.figure(figsize = (10,7))\n",
    "    y = records[w]\n",
    "    plt.imshow(y)\n",
    "    \n",
    "interact(update, w = IntSlider(min=0, max = 999, step = 1, value = 0))"
   ]
  },
  {
   "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.11.4"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}