TugJulia/proto/FTCS.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 109,
   "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": 124,
   "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.random.random((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": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAGdCAYAAABKG5eZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkE0lEQVR4nO3df3RU9Z3/8dfFhAlyklGUJDMQQuQglB+HQkACKj+kBENFWKmg7oGwtrZsqRVTTiG2HnH/aLCtLgdBWbv8kLWLnG4IZDdsJRzzQ0tgQRLXWsS4RpKVpBw4JQO4DAE+3z/8ZuqYmcGBmSSf+Hycc8/x3vv+fPKeTyKv3Jk7E8cYYwQAgCV6dXUDAABEg+ACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFgloasbiJUrV67oxIkTSk5OluM4Xd0OACAKxhidPXtWXq9XvXpFvqbqMcF14sQJZWRkdHUbAIDr0NTUpIEDB0as6THBlZycLEm6S7OVoMQu7gYAEI1LatPb2hP4tzySHhNc7U8PJihRCQ7BBQBW+f+fmvtVXurh5gwAgFUILgCAVeIWXC+99JKysrKUlJSk7OxsvfXWWxHrq6qqlJ2draSkJN12223auHFjvFoDAFgsLsG1Y8cOLV++XD/72c9UW1uru+++W3l5eWpsbAxZ39DQoNmzZ+vuu+9WbW2tnnrqKf34xz9WcXFxPNoDAFjMiccfkpw4caLGjRunl19+OXDsG9/4hubNm6eioqIO9StXrlRpaamOHj0aOLZ06VK9++67qqmp+Upf0+fzye12a5rmcnMGAFjmkmlTpXartbVVKSkpEWtjfsV18eJFvfPOO8rNzQ06npubq/3794ccU1NT06F+1qxZOnz4sNra2kKO8fv98vl8QRsAoOeLeXCdOnVKly9fVlpaWtDxtLQ0tbS0hBzT0tISsv7SpUs6depUyDFFRUVyu92BjTcfA8DXQ9xuzvjyvfjGmIj354eqD3W8XWFhoVpbWwNbU1PTdXYMALBBzN+AfOutt+qGG27ocHV18uTJDldV7dLT00PWJyQk6JZbbgk5xuVyyeVyxaZpAIA1Yn7F1bt3b2VnZ6u8vDzoeHl5uSZPnhxyzKRJkzrU7927V+PHj1diIjdaAAD+Ki5PFRYUFOif//mftXnzZh09elRPPvmkGhsbtXTpUkmfP823ePHiQP3SpUt1/PhxFRQU6OjRo9q8ebM2bdqkFStWxKM9AIDF4vJZhQsXLtTp06f1D//wD2pubtaoUaO0Z88eZWZmSpKam5uD3tOVlZWlPXv26Mknn9SGDRvk9Xq1bt06zZ8/Px7tAQAsFpf3cXUF3scFAPbq0vdxAQAQTwQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqMQ+uoqIiTZgwQcnJyUpNTdW8efN07NixiGMqKyvlOE6H7YMPPoh1ewAAy8U8uKqqqrRs2TIdOHBA5eXlunTpknJzc3X+/Pmrjj127Jiam5sD29ChQ2PdHgDAcgmxnvD3v/990P6WLVuUmpqqd955R1OmTIk4NjU1VTfddFOsWwIA9CBxf42rtbVVktSvX7+r1o4dO1Yej0czZsxQRUVFxFq/3y+fzxe0AQB6vrgGlzFGBQUFuuuuuzRq1KiwdR6PR6+88oqKi4u1c+dODRs2TDNmzFB1dXXYMUVFRXK73YEtIyMjHg8BANDNOMYYE6/Jly1bprKyMr399tsaOHBgVGPnzJkjx3FUWloa8rzf75ff7w/s+3w+ZWRkaJrmKsFJvK6+AQCd65JpU6V2q7W1VSkpKRFr43bF9fjjj6u0tFQVFRVRh5Yk5eTkqL6+Pux5l8ullJSUoA0A0PPF/OYMY4wef/xxlZSUqLKyUllZWdc0T21trTweT4y7AwDYLubBtWzZMv3rv/6rdu/ereTkZLW0tEiS3G63+vTpI0kqLCzUp59+qm3btkmS1q5dq8GDB2vkyJG6ePGiXnvtNRUXF6u4uDjW7QEALBfz4Hr55ZclSdOmTQs6vmXLFi1ZskSS1NzcrMbGxsC5ixcvasWKFfr000/Vp08fjRw5UmVlZZo9e3as2wMAWC6uN2d0Jp/PJ7fbzc0ZAGChbnFzBgAA8UBwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsEvPgWr16tRzHCdrS09MjjqmqqlJ2draSkpJ02223aePGjbFuCwDQQyTEY9KRI0dq3759gf0bbrghbG1DQ4Nmz56txx57TK+99pr+8Ic/6Ic//KH69++v+fPnx6M9AIDF4hJcCQkJV73Kardx40YNGjRIa9eulSR94xvf0OHDh/XrX/+a4AIAdBCX17jq6+vl9XqVlZWlhx56SB9//HHY2pqaGuXm5gYdmzVrlg4fPqy2traw4/x+v3w+X9AGAOj5Yh5cEydO1LZt2/TGG2/oN7/5jVpaWjR58mSdPn06ZH1LS4vS0tKCjqWlpenSpUs6depU2K9TVFQkt9sd2DIyMmL6OAAA3VPMgysvL0/z58/X6NGj9a1vfUtlZWWSpFdffTXsGMdxgvaNMSGPf1FhYaFaW1sDW1NTUwy6BwB0d3F5jeuL+vbtq9GjR6u+vj7k+fT0dLW0tAQdO3nypBISEnTLLbeEndflcsnlcsW0VwBA9xf393H5/X4dPXpUHo8n5PlJkyapvLw86NjevXs1fvx4JSYmxrs9AIBlYh5cK1asUFVVlRoaGnTw4EF95zvfkc/nU35+vqTPn+JbvHhxoH7p0qU6fvy4CgoKdPToUW3evFmbNm3SihUrYt0aAKAHiPlThf/7v/+rhx9+WKdOnVL//v2Vk5OjAwcOKDMzU5LU3NysxsbGQH1WVpb27NmjJ598Uhs2bJDX69W6deu4FR4AEJJj2u+EsJzP55Pb7dY0zVWCw1OMAGCTS6ZNldqt1tZWpaSkRKzlswoBAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAVol5cA0ePFiO43TYli1bFrK+srIyZP0HH3wQ69YAAD1AQqwnPHTokC5fvhzY/+Mf/6iZM2fqwQcfjDju2LFjSklJCez3798/1q0BAHqAmAfXlwNnzZo1GjJkiKZOnRpxXGpqqm666aZYtwMA6GHi+hrXxYsX9dprr+nRRx+V4zgRa8eOHSuPx6MZM2aooqIinm0BACwW8yuuL9q1a5fOnDmjJUuWhK3xeDx65ZVXlJ2dLb/fr3/5l3/RjBkzVFlZqSlTpoQd5/f75ff7A/s+ny+WrQMAuinHGGPiNfmsWbPUu3dv/fu//3tU4+bMmSPHcVRaWhq2ZvXq1Xr22Wc7HJ+muUpwEqPuFQDQdS6ZNlVqt1pbW4Pudwglbk8VHj9+XPv27dP3vve9qMfm5OSovr4+Yk1hYaFaW1sDW1NT07W2CgCwSNyeKtyyZYtSU1P17W9/O+qxtbW18ng8EWtcLpdcLte1tgcAsFRcguvKlSvasmWL8vPzlZAQ/CUKCwv16aefatu2bZKktWvXavDgwRo5cmTgZo7i4mIVFxfHozUAgOXiElz79u1TY2OjHn300Q7nmpub1djYGNi/ePGiVqxYoU8//VR9+vTRyJEjVVZWptmzZ8ejNQCA5eJ6c0Zn8vl8crvd3JwBABbqFjdnAAAQDwQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCpRB1d1dbXmzJkjr9crx3G0a9euoPPGGK1evVper1d9+vTRtGnT9P7771913uLiYo0YMUIul0sjRoxQSUlJtK0BAL4Gog6u8+fPa8yYMVq/fn3I87/85S/1wgsvaP369Tp06JDS09M1c+ZMnT17NuycNTU1WrhwoRYtWqR3331XixYt0oIFC3Tw4MFo2wMA9HCOMcZc82DHUUlJiebNmyfp86str9er5cuXa+XKlZIkv9+vtLQ0Pffcc/rBD34Qcp6FCxfK5/PpP//zPwPH7r33Xt18883avn37V+rF5/PJ7XZrmuYqwUm81ocEAOgCl0ybKrVbra2tSklJiVgb09e4Ghoa1NLSotzc3MAxl8ulqVOnav/+/WHH1dTUBI2RpFmzZkUc4/f75fP5gjYAQM8X0+BqaWmRJKWlpQUdT0tLC5wLNy7aMUVFRXK73YEtIyPjOjoHANgiLncVOo4TtG+M6XDsescUFhaqtbU1sDU1NV17wwAAayTEcrL09HRJn19BeTyewPGTJ092uKL68rgvX11dbYzL5ZLL5brOjgEAtonpFVdWVpbS09NVXl4eOHbx4kVVVVVp8uTJYcdNmjQpaIwk7d27N+IYAMDXU9RXXOfOndNHH30U2G9oaFBdXZ369eunQYMGafny5frFL36hoUOHaujQofrFL36hG2+8UY888khgzOLFizVgwAAVFRVJkp544glNmTJFzz33nObOnavdu3dr3759evvtt2PwEAEAPUnUwXX48GFNnz49sF9QUCBJys/P19atW/XTn/5U//d//6cf/vCH+stf/qKJEydq7969Sk5ODoxpbGxUr15/vdibPHmyXn/9df385z/X008/rSFDhmjHjh2aOHHi9Tw2AEAPdF3v4+pOeB8XANiry97HBQBAvBFcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrRB1c1dXVmjNnjrxerxzH0a5duwLn2tratHLlSo0ePVp9+/aV1+vV4sWLdeLEiYhzbt26VY7jdNguXLgQ9QMCAPRsUQfX+fPnNWbMGK1fv77Duc8++0xHjhzR008/rSNHjmjnzp368MMPdf/991913pSUFDU3NwdtSUlJ0bYHAOjhEqIdkJeXp7y8vJDn3G63ysvLg469+OKLuuOOO9TY2KhBgwaFnddxHKWnp0fbDgDgaybur3G1trbKcRzddNNNEevOnTunzMxMDRw4UPfdd59qa2sj1vv9fvl8vqANANDzxTW4Lly4oFWrVumRRx5RSkpK2Lrhw4dr69atKi0t1fbt25WUlKQ777xT9fX1YccUFRXJ7XYHtoyMjHg8BABAN+MYY8w1D3YclZSUaN68eR3OtbW16cEHH1RjY6MqKysjBteXXblyRePGjdOUKVO0bt26kDV+v19+vz+w7/P5lJGRoWmaqwQnMerHAgDoOpdMmyq1W62trVfNi6hf4/oq2tratGDBAjU0NOjNN9+MKrQkqVevXpowYULEKy6XyyWXy3W9rQIALBPzpwrbQ6u+vl779u3TLbfcEvUcxhjV1dXJ4/HEuj0AgOWivuI6d+6cPvroo8B+Q0OD6urq1K9fP3m9Xn3nO9/RkSNH9B//8R+6fPmyWlpaJEn9+vVT7969JUmLFy/WgAEDVFRUJEl69tlnlZOTo6FDh8rn82ndunWqq6vThg0bYvEYAQA9SNTBdfjwYU2fPj2wX1BQIEnKz8/X6tWrVVpaKkn65je/GTSuoqJC06ZNkyQ1NjaqV6+/XuydOXNG3//+99XS0iK3262xY8equrpad9xxR7TtAQB6uOu6OaM78fl8crvd3JwBABaK5uYMPqsQAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYJWog6u6ulpz5syR1+uV4zjatWtX0PklS5bIcZygLScn56rzFhcXa8SIEXK5XBoxYoRKSkqibQ0A8DUQdXCdP39eY8aM0fr168PW3HvvvWpubg5se/bsiThnTU2NFi5cqEWLFundd9/VokWLtGDBAh08eDDa9gAAPVxCtAPy8vKUl5cXscblcik9Pf0rz7l27VrNnDlThYWFkqTCwkJVVVVp7dq12r59e7QtAgB6sLi8xlVZWanU1FTdfvvteuyxx3Ty5MmI9TU1NcrNzQ06NmvWLO3fvz/sGL/fL5/PF7QBAHq+mAdXXl6efvvb3+rNN9/U888/r0OHDumee+6R3+8PO6alpUVpaWlBx9LS0tTS0hJ2TFFRkdxud2DLyMiI2WMAAHRfUT9VeDULFy4M/PeoUaM0fvx4ZWZmqqysTA888EDYcY7jBO0bYzoc+6LCwkIVFBQE9n0+H+EFAF8DMQ+uL/N4PMrMzFR9fX3YmvT09A5XVydPnuxwFfZFLpdLLpcrZn0CAOwQ9/dxnT59Wk1NTfJ4PGFrJk2apPLy8qBje/fu1eTJk+PdHgDAMlFfcZ07d04fffRRYL+hoUF1dXXq16+f+vXrp9WrV2v+/PnyeDz65JNP9NRTT+nWW2/V3/zN3wTGLF68WAMGDFBRUZEk6YknntCUKVP03HPPae7cudq9e7f27dunt99+OwYPEQDQk0QdXIcPH9b06dMD++2vM+Xn5+vll1/We++9p23btunMmTPyeDyaPn26duzYoeTk5MCYxsZG9er114u9yZMn6/XXX9fPf/5zPf300xoyZIh27NihiRMnXs9jAwD0QI4xxnR1E7Hg8/nkdrs1TXOV4CR2dTsAgChcMm2q1G61trYqJSUlYi2fVQgAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwStTBVV1drTlz5sjr9cpxHO3atSvovOM4Ibdf/epXYefcunVryDEXLlyI+gEBAHq2qIPr/PnzGjNmjNavXx/yfHNzc9C2efNmOY6j+fPnR5w3JSWlw9ikpKRo2wMA9HAJ0Q7Iy8tTXl5e2PPp6elB+7t379b06dN12223RZzXcZwOYwEA+LK4vsb15z//WWVlZfrud7971dpz584pMzNTAwcO1H333afa2tqI9X6/Xz6fL2gDAPR8cQ2uV199VcnJyXrggQci1g0fPlxbt25VaWmptm/frqSkJN15552qr68PO6aoqEhutzuwZWRkxLp9AEA35BhjzDUPdhyVlJRo3rx5Ic8PHz5cM2fO1IsvvhjVvFeuXNG4ceM0ZcoUrVu3LmSN3++X3+8P7Pt8PmVkZGia5irBSYzq6wEAutYl06ZK7VZra6tSUlIi1kb9GtdX9dZbb+nYsWPasWNH1GN79eqlCRMmRLzicrlccrlc19MiAMBCcXuqcNOmTcrOztaYMWOiHmuMUV1dnTweTxw6AwDYLOorrnPnzumjjz4K7Dc0NKiurk79+vXToEGDJH3+tN3vfvc7Pf/88yHnWLx4sQYMGKCioiJJ0rPPPqucnBwNHTpUPp9P69atU11dnTZs2HAtjwkA0INFHVyHDx/W9OnTA/sFBQWSpPz8fG3dulWS9Prrr8sYo4cffjjkHI2NjerV668Xe2fOnNH3v/99tbS0yO12a+zYsaqurtYdd9wRbXsAgB7uum7O6E58Pp/cbjc3ZwCAhaK5OYPPKgQAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYJW5/j6urlHz4nlKSry+PZ3m/GZtmAOBr4I0Tddc9h+/sFd18+1er5YoLAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYJUe8xeQjTGSJN+5K9c91yXTdt1zAMDXhe/s9f+72/5vd/u/5ZH0mOA6e/asJClz3CcxmO3jGMwBAF8PN98eu7nOnj0rt9sdscYxXyXeLHDlyhWdOHFCycnJchwnZI3P51NGRoaampqUkpLSyR1eO/rufLb2Tt+di75jxxijs2fPyuv1qlevyK9i9Zgrrl69emngwIFfqTYlJaXbfLOiQd+dz9be6btz0XdsXO1Kqx03ZwAArEJwAQCs8rUKLpfLpWeeeUYul6urW4kKfXc+W3un785F312jx9ycAQD4evhaXXEBAOxHcAEArEJwAQCsQnABAKzS44LrpZdeUlZWlpKSkpSdna233norYn1VVZWys7OVlJSk2267TRs3buykTj9XVFSkCRMmKDk5WampqZo3b56OHTsWcUxlZaUcx+mwffDBB53UtbR69eoOXz89PT3imK5e63aDBw8OuX7Lli0LWd9V611dXa05c+bI6/XKcRzt2rUr6LwxRqtXr5bX61WfPn00bdo0vf/++1edt7i4WCNGjJDL5dKIESNUUlLSaX23tbVp5cqVGj16tPr27Suv16vFixfrxIkTEefcunVryO/BhQsXOqVvSVqyZEmHr5+Tk3PVebtyvSWFXDfHcfSrX/0q7Jydsd7Xo0cF144dO7R8+XL97Gc/U21tre6++27l5eWpsbExZH1DQ4Nmz56tu+++W7W1tXrqqaf04x//WMXFxZ3Wc1VVlZYtW6YDBw6ovLxcly5dUm5urs6fP3/VsceOHVNzc3NgGzp0aCd0/FcjR44M+vrvvfde2NrusNbtDh06FNR3eXm5JOnBBx+MOK6z1/v8+fMaM2aM1q9fH/L8L3/5S73wwgtav369Dh06pPT0dM2cOTPwuZ2h1NTUaOHChVq0aJHeffddLVq0SAsWLNDBgwc7pe/PPvtMR44c0dNPP60jR45o586d+vDDD3X//fdfdd6UlJSg9W9ublZSUlKn9N3u3nvvDfr6e/bsiThnV6+3pA5rtnnzZjmOo/nz50ecN97rfV1MD3LHHXeYpUuXBh0bPny4WbVqVcj6n/70p2b48OFBx37wgx+YnJycuPV4NSdPnjSSTFVVVdiaiooKI8n85S9/6bzGvuSZZ54xY8aM+cr13XGt2z3xxBNmyJAh5sqVKyHPd4f1lmRKSkoC+1euXDHp6elmzZo1gWMXLlwwbrfbbNy4Mew8CxYsMPfee2/QsVmzZpmHHnoo5j0b07HvUP7rv/7LSDLHjx8PW7Nlyxbjdrtj21wEofrOz883c+fOjWqe7rjec+fONffcc0/Ems5e72j1mCuuixcv6p133lFubm7Q8dzcXO3fvz/kmJqamg71s2bN0uHDh9XW1jV/2qS1tVWS1K9fv6vWjh07Vh6PRzNmzFBFRUW8W+ugvr5eXq9XWVlZeuihh/Txx+E/Vb87rrX0+c/Na6+9pkcffTTshzO36+r1/qKGhga1tLQEranL5dLUqVPD/rxL4b8PkcbEW2trqxzH0U033RSx7ty5c8rMzNTAgQN13333qba2tnMa/ILKykqlpqbq9ttv12OPPaaTJ09GrO9u6/3nP/9ZZWVl+u53v3vV2u6w3uH0mOA6deqULl++rLS0tKDjaWlpamlpCTmmpaUlZP2lS5d06tSpuPUajjFGBQUFuuuuuzRq1KiwdR6PR6+88oqKi4u1c+dODRs2TDNmzFB1dXWn9Tpx4kRt27ZNb7zxhn7zm9+opaVFkydP1unTp0PWd7e1brdr1y6dOXNGS5YsCVvTHdb7y9p/pqP5eW8fF+2YeLpw4YJWrVqlRx55JOKHvQ4fPlxbt25VaWmptm/frqSkJN15552qr6/vtF7z8vL029/+Vm+++aaef/55HTp0SPfcc4/8fn/YMd1tvV999VUlJyfrgQceiFjXHdY7kh7z6fDtvvxbszEm4m/SoepDHe8MP/rRj/Tf//3fevvttyPWDRs2TMOGDQvsT5o0SU1NTfr1r3+tKVOmxLtNSZ//T9xu9OjRmjRpkoYMGaJXX31VBQUFIcd0p7Vut2nTJuXl5cnr9Yat6Q7rHU60P+/XOiYe2tra9NBDD+nKlSt66aWXItbm5OQE3Qhx5513aty4cXrxxRe1bt26eLcqSVq4cGHgv0eNGqXx48crMzNTZWVlEYOgu6y3JG3evFl/+7d/e9XXqrrDekfSY664br31Vt1www0dfpM5efJkh9942qWnp4esT0hI0C233BK3XkN5/PHHVVpaqoqKiq/851m+KCcnp0t/G+rbt69Gjx4dtofutNbtjh8/rn379ul73/te1GO7er3b7+CM5ue9fVy0Y+Khra1NCxYsUENDg8rLy6P+0xq9evXShAkTuvR74PF4lJmZGbGH7rLekvTWW2/p2LFj1/Tz3h3W+4t6THD17t1b2dnZgTvE2pWXl2vy5Mkhx0yaNKlD/d69ezV+/HglJibGrdcvMsboRz/6kXbu3Kk333xTWVlZ1zRPbW2tPB5PjLv76vx+v44ePRq2h+6w1l+2ZcsWpaam6tvf/nbUY7t6vbOyspSenh60phcvXlRVVVXYn3cp/Pch0phYaw+t+vp67du375p+cTHGqK6urku/B6dPn1ZTU1PEHrrDerfbtGmTsrOzNWbMmKjHdof1DtJVd4XEw+uvv24SExPNpk2bzJ/+9CezfPly07dvX/PJJ58YY4xZtWqVWbRoUaD+448/NjfeeKN58sknzZ/+9CezadMmk5iYaP7t3/6t03r++7//e+N2u01lZaVpbm4ObJ999lmg5st9/+M//qMpKSkxH374ofnjH/9oVq1aZSSZ4uLiTuv7Jz/5iamsrDQff/yxOXDggLnvvvtMcnJyt17rL7p8+bIZNGiQWblyZYdz3WW9z549a2pra01tba2RZF544QVTW1sbuPtuzZo1xu12m507d5r33nvPPPzww8bj8RifzxeYY9GiRUF31f7hD38wN9xwg1mzZo05evSoWbNmjUlISDAHDhzolL7b2trM/fffbwYOHGjq6uqCfub9fn/YvlevXm1+//vfm//5n/8xtbW15u/+7u9MQkKCOXjwYKf0ffbsWfOTn/zE7N+/3zQ0NJiKigozadIkM2DAgG693u1aW1vNjTfeaF5++eWQc3TFel+PHhVcxhizYcMGk5mZaXr37m3GjRsXdFt5fn6+mTp1alB9ZWWlGTt2rOndu7cZPHhw2G9svEgKuW3ZsiVs388995wZMmSISUpKMjfffLO56667TFlZWaf2vXDhQuPxeExiYqLxer3mgQceMO+//37Yno3p+rX+ojfeeMNIMseOHetwrrusd/tt+F/e8vPzjTGf3xL/zDPPmPT0dONyucyUKVPMe++9FzTH1KlTA/Xtfve735lhw4aZxMREM3z48JgHcKS+Gxoawv7MV1RUhO17+fLlZtCgQaZ3796mf//+Jjc31+zfv7/T+v7ss89Mbm6u6d+/v0lMTDSDBg0y+fn5prGxMWiO7rbe7f7pn/7J9OnTx5w5cybkHF2x3teDP2sCALBKj3mNCwDw9UBwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKzy/wBVnTUJwLaO3QAAAABJRU5ErkJggg==",
      "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'],grid_size['y']))\n",
    "C_t[19,0] = 2000\n",
    "C_t[19,19] = 2000\n",
    "\n",
    "\n",
    "plt.imshow(C_t, vmin=0, vmax=2000)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "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"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation for constant boundary conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulate one time step\n",
    "def simulate_constant_boundary(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['y']-1): #rows\n",
    "                for j in range(1, grid_size['x']-1): #columns\n",
    "                        C_t1[i,j] = C_t[i,j] \\\n",
    "                                + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
    "                                        - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j]) \\\n",
    "                                + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
    "                                        - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
    "                        \n",
    "        # boundary conditions\n",
    "        # left without corners / looping over rows\n",
    "        for i in range(1, grid_size['y']-1):\n",
    "               j = 0\n",
    "               C_t1[i,0] = C_t[i,0] \\\n",
    "                        + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1] \n",
    "                                - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + 2 * alpha_x[i,j]) * C_t[i,j]\n",
    "                                + 2 * alpha_x[i,j] * bc_left) \\\n",
    "                        + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
    "                                        - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j])        \n",
    "\n",
    "        # right without corners / looping over rows\n",
    "        n = grid_size['x']-1 # maximum index in x-direction (columns)\n",
    "        for i in range(1, grid_size['y']-1):\n",
    "               j = n\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",
    "                        + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
    "                                        - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j])\n",
    "\n",
    "        # top without corners / looping over columns\n",
    "        for j in range(1, grid_size['x']-1):\n",
    "               i = 0\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",
    "                        + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
    "                                        - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
    "\n",
    "        # bottom without corners / looping over columns\n",
    "        m = grid_size['y']-1 # maximum index in y-direction (rows)\n",
    "        for j in range(1, grid_size['x']-1):\n",
    "               i = m\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",
    "                        + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
    "                                        - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
    "\n",
    "        # corner top left\n",
    "        i = 0\n",
    "        j = i\n",
    "        C_t1[i,j] = C_t[i,j] \\\n",
    "                        + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1] \n",
    "                                - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + 2 * alpha_x[i,j]) * C_t[i,j]\n",
    "                                + 2 * alpha_x[i,j] * bc_left) \\\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",
    "\n",
    "        # corner top right\n",
    "        i = 0\n",
    "        j = grid_size['x']-1\n",
    "        n = j\n",
    "        C_t1[i,j] = C_t[i,j] \\\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",
    "                        + 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",
    "        # corner bottom left\n",
    "        i = grid_size['y']-1\n",
    "        m = i\n",
    "        j = 0\n",
    "        C_t1[i,j] = C_t[i,j] \\\n",
    "                        + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1] \n",
    "                                - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + 2 * alpha_x[i,j]) * C_t[i,j]\n",
    "                                + 2 * alpha_x[i,j] * bc_left) \\\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",
    "        # corner bottom right\n",
    "        i = grid_size['y']-1\n",
    "        j = grid_size['x']-1\n",
    "        m = i \n",
    "        n = j\n",
    "        C_t1[i,j] = C_t[i,j] \\\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",
    "                        + 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",
    "        return C_t1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation for closed boundary conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulate one time step\n",
    "def simulate_closed_boundary(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['y']-1): #rows\n",
    "                for j in range(1, grid_size['x']-1): #columns\n",
    "                        C_t1[i,j] = C_t[i,j] \\\n",
    "                                + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
    "                                        - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j]) \\\n",
    "                                + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
    "                                        - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
    "                        \n",
    "         # boundary conditions\n",
    "        # left without corners / looping over rows\n",
    "        for i in range(1, grid_size['y']-1):\n",
    "               j = 0\n",
    "               C_t1[i,0] = C_t[i,0] \\\n",
    "                        + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * (C_t[i,j+1] - C_t[i,j]))\\\n",
    "                        + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
    "                                        - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j])        \n",
    "\n",
    "        # right without corners / looping over rows\n",
    "        n = grid_size['x']-1 # maximum index in x-direction (columns)\n",
    "        for i in range(1, grid_size['y']-1):\n",
    "               j = n\n",
    "               C_t1[i,n] = C_t[i,n] \\\n",
    "                        + time_step/delta_x**2 * ( \n",
    "                                - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * (C_t[i,n] - C_t[i,n-1])))\\\n",
    "                        + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
    "                                        - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j])\n",
    "\n",
    "        # top without corners / looping over columns\n",
    "        for j in range(1, grid_size['x']-1):\n",
    "               i = 0\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] - C_t[0,j])) \\\n",
    "                        + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
    "                                        - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
    "\n",
    "        # bottom without corners / looping over columns\n",
    "        m = grid_size['y']-1 # maximum index in y-direction (rows)\n",
    "        for j in range(1, grid_size['x']-1):\n",
    "               i = m\n",
    "               C_t1[m,j] = C_t[m,j] \\\n",
    "                        + time_step/delta_y**2 * ( \n",
    "                                - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * (C_t[m,j] - C_t[m-1,j]))) \\\n",
    "                        + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
    "                                        - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
    "                                        + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
    "\n",
    "        # corner top left\n",
    "        i = 0\n",
    "        j = i\n",
    "        C_t1[i,j] = C_t[i,j] \\\n",
    "                        + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * (C_t[i,j+1] - C_t[i,j]))\\\n",
    "                        + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * (C_t[1,j]-C_t[0,j]))\n",
    "                        \n",
    "\n",
    "        # corner top right\n",
    "        i = 0\n",
    "        j = grid_size['x']-1\n",
    "        n = j\n",
    "        C_t1[i,j] = C_t[i,j] \\\n",
    "                        + time_step/delta_x**2 * (- (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * (C_t[i,n] - C_t[i,n-1])))\\\n",
    "                        + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * (C_t[1,j] - C_t[0,j])) \n",
    "        \n",
    "\n",
    "        # corner bottom left\n",
    "        i = grid_size['y']-1\n",
    "        m = i\n",
    "        j = 0\n",
    "        C_t1[i,j] = C_t[i,j] \\\n",
    "                        + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * (C_t[i,j+1] - C_t[i,j])) \\\n",
    "                        + time_step/delta_y**2 * (- (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j] * (C_t[m,j] - C_t[m-1,j]))))\n",
    "\n",
    "        # corner bottom right\n",
    "        i = grid_size['y']-1\n",
    "        j = grid_size['x']-1\n",
    "        m = i \n",
    "        n = j\n",
    "        C_t1[i,j] = C_t[i,j] \\\n",
    "                        + time_step/delta_x**2 * (- alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * (C_t[i,n] - C_t[i,n-1]))\\\n",
    "                        + time_step/delta_y**2 * (- alpha_interblock(alpha_y[m,j], alpha_y[m-1,j] * (C_t[m,j] - C_t[m-1,j])))\n",
    "\n",
    "        return C_t1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAAGiCAYAAACcWg7FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABC7UlEQVR4nO3df1hUZd4/8PeoMCgLk4gwM4lEXkolxCKYgJWYBVJqpq2afRHLxXwyW0Iu09qKfdonrF3TK00zF8UfFD77KGqri8IqEIHmL1p/xeJGgsZIecmM2ArI3N8/XM46MhwYZwZwzvu1130tc87n3HPP0ZwP96+jEkIIEBEREQHo1d0NICIiop6DiQERERFJmBgQERGRhIkBERERSZgYEBERkYSJAREREUmYGBAREZGEiQERERFJmBgQERGRhIkBERERSZgYEBEROVFGRgZGjhwJLy8v+Pn5YfLkyaioqLCIEUIgPT0der0effv2RWxsLE6dOmUR09jYiAULFsDX1xeenp6YNGkSzp8/bxFz+fJlJCYmQqPRQKPRIDExEfX19Ta1l4kBERGRExUVFWH+/Pk4ePAg8vPzcf36dcTFxeHq1atSzAcffIAPP/wQq1atwuHDh6HVavHEE0/gypUrUkxKSgpyc3ORk5ODkpISNDQ0YMKECWhpaZFiZs6cifLycuTl5SEvLw/l5eVITEy0rcGCiIiIukxdXZ0AIIqKioQQQpjNZqHVasXSpUulmGvXrgmNRiM++eQTIYQQ9fX1ws3NTeTk5EgxFy5cEL169RJ5eXlCCCFOnz4tAIiDBw9KMWVlZQKA+Pbbbzvdvj725UE9h9lsxg8//AAvLy+oVKrubg4REdlACIErV65Ar9ejVy/ndWZfu3YNTU1NDqlLCNHm+0atVkOtVsteZzQaAQA+Pj4AgKqqKhgMBsTFxVnUM2bMGJSWluKll17C0aNH0dzcbBGj1+sREhKC0tJSxMfHo6ysDBqNBqNGjZJioqKioNFoUFpaiuDg4E59LpdJDH744QcEBAR0dzOIiMgONTU1GDRokFPqvnbtGoICfwFDXUvHwZ3wi1/8Ag0NDRbH3nnnHaSnp7d7jRACqampePjhhxESEgIAMBgMAAB/f3+LWH9/f5w7d06KcXd3R//+/dvEtF5vMBjg5+fX5j39/PykmM5wmcTAy8sLAPAwnkIflVs3t4aIiGxxXTSjBLulf8udoampCYa6FlQdDYS3l329EqYrZgRFnENNTQ28vb2l4x31Frzyyiv4+9//jpKSkjbnbu19sNYjcatbY6zFd6aem7lMYtD6ofuo3JgYEBHdiYT1LzZH8/bqZXdiINXl7W2RGMhZsGABdu3aheLiYoteEa1WC+DGb/w6nU46XldXJ/UiaLVaNDU14fLlyxa9BnV1dYiJiZFiLl682OZ9f/zxxza9EXK4KoGIiBSlRZgdUjpLCIFXXnkF27dvx/79+xEUFGRxPigoCFqtFvn5+dKxpqYmFBUVSV/6ERERcHNzs4ipra3FyZMnpZjo6GgYjUZ8/fXXUsyhQ4dgNBqlmM5wmR4DIiKizjBDwAxhdx2dNX/+fHz22WfYuXMnvLy8pPF+jUaDvn37QqVSISUlBe+99x6GDh2KoUOH4r333kO/fv0wc+ZMKXbOnDlYuHAhBgwYAB8fH6SlpSE0NBSPP/44AOD+++/H+PHjkZycjLVr1wIA5s6diwkTJnR64iHgxB6D1atXIygoCB4eHoiIiMCXX34pG19UVISIiAh4eHjg3nvvxSeffOKsphERkYKZHfS/zlqzZg2MRiNiY2Oh0+mksnXrVilm0aJFSElJwcsvv4zIyEhcuHAB+/bts5hzsXz5ckyePBnTpk3D6NGj0a9fP3zxxRfo3bu3FJOdnY3Q0FDExcUhLi4ODz74IDZv3mzT/VEJIexLm6zYunUrEhMTsXr1aowePRpr167Fn/70J5w+fRqDBw9uE19VVYWQkBAkJyfjpZdewldffYWXX34Zn3/+OaZOndqp9zSZTNBoNIhVTeYcAyKiO8x10YxCsQNGo7HTY/a2av2e+KFikEMmH+qDzzu1vd3FKYnBqFGjMGLECKxZs0Y6dv/992Py5MnIyMhoE//6669j165dOHPmjHRs3rx5+Oabb1BWVtap92RiQER05+rKxKDm27sdkhgE3HfBJRMDhw8lNDU14ejRoxabMABAXFwcSktLrV5TVlbWJj4+Ph5HjhxBc3Oz1WsaGxthMpksChERUUda5xjYW1yVwxODn376CS0tLVY3amhvgwWDwWA1/vr16/jpp5+sXpORkSE9JEKj0XBzIyIiIgdw2uRDWzdqsBZv7XirJUuWwGg0SqWmpsbOFhMRkRKYIdBiZ3HlHgOHL1f09fVF79692/QO3LxRw620Wq3V+D59+mDAgAFWr+nMftRERES36urlincah/cYuLu7IyIiwmITBgDIz89vd4OF6OjoNvH79u1DZGQk3Nw4kZCIiKirOGUoITU1FX/605+wfv16nDlzBq+99hqqq6sxb948ADeGAWbNmiXFz5s3D+fOnUNqairOnDmD9evXIzMzE2lpac5oHhERKViLEA4prsopOx9Onz4dly5dwn//93+jtrYWISEh2LNnDwIDAwHc2Maxurpaig8KCsKePXvw2muv4eOPP4Zer8dHH33U6T0MiIiIOsv872JvHa7KKfsYdAfuY0BEdOfqyn0Mvj3jDy879zG4csWM++6/6JL7GPBZCUREpCitKwvsrcNVMTEgIiJFaRE3ir11uComBkREpCicYyDPaRscERER0Z2HPQZERKQoZqjQgvZ34u1sHa6KiQERESmKWdwo9tbhqjiUQERERBL2GBARkaK0OGAowd7rezImBkREpChMDORxKIGIiIgk7DEgIiJFMQsVzMLOVQl2Xt+TMTEgIiJF4VCCPA4lEBERkYQ9BkREpCgt6IUWO38vbnFQW3oiJgZERKQowgFzDATnGBAREbkGzjGQxzkGREREJGGPARERKUqL6IUWYeccAxd+VgITAyIiUhQzVDDb2WFuhutmBhxKICIiIgl7DIiISFE4+VAeEwMiIlIUx8wx4FACERERKQB7DIiISFFuTD608yFKHEogIiJyDWYHbInMVQlERESkCEwMiIhIUVonH9pbbFFcXIyJEydCr9dDpVJhx44dFudVKpXV8oc//EGKiY2NbXN+xowZFvVcvnwZiYmJ0Gg00Gg0SExMRH19vU1tZWJARESKYkYvhxRbXL16FWFhYVi1apXV87W1tRZl/fr1UKlUmDp1qkVccnKyRdzatWstzs+cORPl5eXIy8tDXl4eysvLkZiYaFNbOceAiIgUpUWo0GLn0xFtvT4hIQEJCQntntdqtRavd+7cibFjx+Lee++1ON6vX782sa3OnDmDvLw8HDx4EKNGjQIArFu3DtHR0aioqEBwcHCn2soeAyIiottkMpksSmNjo911Xrx4Ebt378acOXPanMvOzoavry+GDx+OtLQ0XLlyRTpXVlYGjUYjJQUAEBUVBY1Gg9LS0k6/P3sMiIhIUVocsCqh5d+rEgICAiyOv/POO0hPT7er7o0bN8LLywtTpkyxOP78888jKCgIWq0WJ0+exJIlS/DNN98gPz8fAGAwGODn59emPj8/PxgMhk6/PxMDIiJSFLPoBbOdOx+a/73zYU1NDby9vaXjarXarnoBYP369Xj++efh4eFhcTw5OVn6OSQkBEOHDkVkZCSOHTuGESNGALgxifFWQgirx9vDoQQiIqLb5O3tbVHsTQy+/PJLVFRU4Ne//nWHsSNGjICbmxsqKysB3JincPHixTZxP/74I/z9/TvdBiYGRESkKK1DCfYWZ8jMzERERATCwsI6jD116hSam5uh0+kAANHR0TAajfj666+lmEOHDsFoNCImJqbTbeBQAhERKYoZtq8qsFaHLRoaGnD27FnpdVVVFcrLy+Hj44PBgwcDuDGR8c9//jOWLVvW5vp//vOfyM7OxpNPPglfX1+cPn0aCxcuRHh4OEaPHg0AuP/++zF+/HgkJydLyxjnzp2LCRMmdHpFAuCEHoOMjAyMHDkSXl5e8PPzw+TJk1FRUSF7TWFhodWNHb799ltHN4+IiKjLHTlyBOHh4QgPDwcApKamIjw8HG+//bYUk5OTAyEEnnvuuTbXu7u7429/+xvi4+MRHByMV199FXFxcSgoKEDv3r2luOzsbISGhiIuLg5xcXF48MEHsXnzZpvaqhLCsc+OHD9+PGbMmIGRI0fi+vXrePPNN3HixAmcPn0anp6eVq8pLCzE2LFjUVFRYTGJY+DAgRYfWI7JZIJGo0GsajL6qNwc8lmIiKhrXBfNKBQ7YDQaLb4HHKn1e2LNsZHo+wv7Osz/1XAd/zXisFPb210cPpSQl5dn8XrDhg3w8/PD0aNH8eijj8pe6+fnh7vuusvRTSIiIpLczpbG1upwVU7/ZEajEQDg4+PTYWx4eDh0Oh3GjRuHAwcOyMY2Nja22ViCiIiI7OPUxEAIgdTUVDz88MMICQlpN06n0+HTTz/Ftm3bsH37dgQHB2PcuHEoLi5u95qMjAzpIREajabNJhNERETWmKFySHFVDp9jcLP58+dj9+7dKCkpwaBBg2y6duLEiVCpVNi1a5fV842NjRZbT5pMJgQEBHCOARHRHagr5xgsPxLjkDkGr0WWco6BLRYsWIBdu3ahuLjY5qQAuLG/85YtW9o9r1arHbLDFBERKYtjtkR23TkGDk8MhBBYsGABcnNzUVhYiKCgoNuq5/jx49KmDURERNQ1HJ4YzJ8/H5999hl27twJLy8v6cENGo0Gffv2BQAsWbIEFy5cwKZNmwAAK1aswD333IPhw4ejqakJW7ZswbZt27Bt2zZHN4+IiBTOLFQw27vBkZ3X92QOTwzWrFkDAIiNjbU4vmHDBsyePRsAUFtbi+rqaulcU1MT0tLScOHCBfTt2xfDhw/H7t278eSTTzq6eUREpHBmBwwlmDmU0HmdmcuYlZVl8XrRokVYtGiRo5tCRERENuKzEoiISFEc89hl9hgQERG5hBao0GLnPgT2Xt+TuW7KQ0RERDZjjwERESkKhxLkMTEgIiJFaYH9QwEtjmlKj+S6KQ8RERHZjD0GRESkKBxKkMfEgIiIFKVF9EKLnV/s9l7fkzExICIiRREOeGyy4HJFIiIiUgL2GBARkaJwKEEeEwMiIlIUPl1RnuumPERERGQz9hgQEZGitDjgscv2Xt+TMTEgIiJF4VCCPNdNeYiIiMhm7DEgIiJFMaMXzHb+Xmzv9T0ZEwMiIlKUFqFCi51DAfZe35O5bspDRERENmOPARERKQonH8pjYkBERIoiHPB0RcGdD4mIiFxDC1RosfMhSPZe35O5bspDRERENmOPARERKYpZ2D9HwCwc1JgeiIkBEREpitkBcwzsvb4nc91PRkRERDZjYkBERIpihsohxRbFxcWYOHEi9Ho9VCoVduzYYXF+9uzZUKlUFiUqKsoiprGxEQsWLICvry88PT0xadIknD9/3iLm8uXLSExMhEajgUajQWJiIurr621qKxMDIiJSlNadD+0ttrh69SrCwsKwatWqdmPGjx+P2tpaqezZs8fifEpKCnJzc5GTk4OSkhI0NDRgwoQJaGlpkWJmzpyJ8vJy5OXlIS8vD+Xl5UhMTLSprZxjQERE5GQJCQlISEiQjVGr1dBqtVbPGY1GZGZmYvPmzXj88ccBAFu2bEFAQAAKCgoQHx+PM2fOIC8vDwcPHsSoUaMAAOvWrUN0dDQqKioQHBzcqbayx4CIiBSldfKhvQUATCaTRWlsbLztdhUWFsLPzw/Dhg1DcnIy6urqpHNHjx5Fc3Mz4uLipGN6vR4hISEoLS0FAJSVlUGj0UhJAQBERUVBo9FIMZ3BxICIiBTFDJW0LfJtl3/PMQgICJDG8zUaDTIyMm6rTQkJCcjOzsb+/fuxbNkyHD58GI899piUaBgMBri7u6N///4W1/n7+8NgMEgxfn5+ber28/OTYjqDQwlERES3qaamBt7e3tJrtVp9W/VMnz5d+jkkJASRkZEIDAzE7t27MWXKlHavE0JApfrPfIebf24vpiPsMSAiIkURDliRIP7dY+Dt7W1RbjcxuJVOp0NgYCAqKysBAFqtFk1NTbh8+bJFXF1dHfz9/aWYixcvtqnrxx9/lGI6g4kBEREpit3DCA54OmNHLl26hJqaGuh0OgBAREQE3NzckJ+fL8XU1tbi5MmTiImJAQBER0fDaDTi66+/lmIOHToEo9EoxXQGhxKIiEhRumPnw4aGBpw9e1Z6XVVVhfLycvj4+MDHxwfp6emYOnUqdDodvv/+e7zxxhvw9fXFM888AwDQaDSYM2cOFi5ciAEDBsDHxwdpaWkIDQ2VVincf//9GD9+PJKTk7F27VoAwNy5czFhwoROr0gAmBgQERE53ZEjRzB27FjpdWpqKgAgKSkJa9aswYkTJ7Bp0ybU19dDp9Nh7Nix2Lp1K7y8vKRrli9fjj59+mDatGn417/+hXHjxiErKwu9e/eWYrKzs/Hqq69KqxcmTZoku3eCNSohhEMfBZGeno7f/e53FsdunjVpTVFREVJTU3Hq1Cno9XosWrQI8+bNs+l9TSYTNBoNYlWT0UfldlttJyKi7nFdNKNQ7IDRaLSYzOdIrd8TT+97EW6e7nbV1Xy1CTvj1ju1vd3FKT0Gw4cPR0FBgfT65mzmVlVVVXjyySeRnJyMLVu24KuvvsLLL7+MgQMHYurUqc5oHhERKdjtbGlsrQ5X5ZTEoE+fPu3u3nSrTz75BIMHD8aKFSsA3BgjOXLkCP74xz/KJgaNjY0WG0mYTCa72kxEREROWpVQWVkJvV6PoKAgzJgxA9999127sWVlZRY7OQFAfHw8jhw5gubm5navy8jIsNhUIiAgwGHtJyIi13UnrEroTg5PDEaNGoVNmzZh7969WLduHQwGA2JiYnDp0iWr8QaDoc36Sn9/f1y/fh0//fRTu++zZMkSGI1GqdTU1Dj0cxARkWtiYiDP4UMJNz8kIjQ0FNHR0RgyZAg2btwozcK81a07MrXOh5TbqUmtVjtsIwkiIiK6wenLFT09PREaGirt3nQrrVbbZsVCXV0d+vTpgwEDBji7eUREpDCO+I3flXsMnL7zYWNjI86cOSPt3nSr6Ohoi52cAGDfvn2IjIyEmxuXHRIRkWNxKEGewxODtLQ0FBUVoaqqCocOHcKzzz4Lk8mEpKQkADfmBsyaNUuKnzdvHs6dO4fU1FScOXMG69evR2ZmJtLS0hzdNCIiIuqAw4cSzp8/j+eeew4//fQTBg4ciKioKBw8eBCBgYEAbuztXF1dLcUHBQVhz549eO211/Dxxx9Dr9fjo48+4h4GRETkFAL270Pg0J0BexiHJwY5OTmy57OystocGzNmDI4dO+bophAREbXBOQby+KwEIiJSFCYG8vjYZSIiIpKwx4CIiBSFPQbymBgQEZGiMDGQx6EEIiIikrDHgIiIFEUIFYSdv/Hbe31PxsSAiIgUxQyV3fsY2Ht9T8ahBCIiIpKwx4CIiBSFkw/lMTEgIiJF4RwDeRxKICIiIgl7DIiISFE4lCCPiQERESkKhxLkMTEgIiJFEQ7oMXDlxIBzDIiIiEjCHgMiIlIUAUAI++twVUwMiIhIUcxQQcWdD9vFoQQiIiKSsMeAiIgUhasS5DExICIiRTELFVTcx6BdHEogIiIiCXsMiIhIUYRwwKoEF16WwMSAiIgUhXMM5HEogYiIiCTsMSAiIkVhj4E8JgZERKQoXJUgj0MJRESkKK2TD+0ttiguLsbEiROh1+uhUqmwY8cO6VxzczNef/11hIaGwtPTE3q9HrNmzcIPP/xgUUdsbCxUKpVFmTFjhkXM5cuXkZiYCI1GA41Gg8TERNTX19vUViYGRERETnb16lWEhYVh1apVbc79/PPPOHbsGN566y0cO3YM27dvxz/+8Q9MmjSpTWxycjJqa2ulsnbtWovzM2fORHl5OfLy8pCXl4fy8nIkJiba1FYOJRARkaLc+I3f3jkGN/7fZDJZHFer1VCr1W3iExISkJCQYLUujUaD/Px8i2MrV67EQw89hOrqagwePFg63q9fP2i1Wqv1nDlzBnl5eTh48CBGjRoFAFi3bh2io6NRUVGB4ODgTn029hgQEZGitE4+tLcAQEBAgNRtr9FokJGR4ZA2Go1GqFQq3HXXXRbHs7Oz4evri+HDhyMtLQ1XrlyRzpWVlUGj0UhJAQBERUVBo9GgtLS00+/NHgMiIqLbVFNTA29vb+m1td4CW127dg2LFy/GzJkzLep+/vnnERQUBK1Wi5MnT2LJkiX45ptvpN4Gg8EAPz+/NvX5+fnBYDB0+v2ZGBARkaKIfxd76wAAb29viy9vezU3N2PGjBkwm81YvXq1xbnk5GTp55CQEAwdOhSRkZE4duwYRowYAQBQqdoOkQghrB5vD4cSiIhIURw5lOBIzc3NmDZtGqqqqpCfn99hwjFixAi4ubmhsrISAKDVanHx4sU2cT/++CP8/f073Q4mBkRERN2sNSmorKxEQUEBBgwY0OE1p06dQnNzM3Q6HQAgOjoaRqMRX3/9tRRz6NAhGI1GxMTEdLotHEogIiJlceRYQic1NDTg7Nmz0uuqqiqUl5fDx8cHer0ezz77LI4dO4a//OUvaGlpkeYE+Pj4wN3dHf/85z+RnZ2NJ598Er6+vjh9+jQWLlyI8PBwjB49GgBw//33Y/z48UhOTpaWMc6dOxcTJkzo9IoEwAk9Bvfcc0+bDRhUKhXmz59vNb6wsNBq/LfffuvophEREQGOGEawcSjhyJEjCA8PR3h4OAAgNTUV4eHhePvtt3H+/Hns2rUL58+fxy9/+UvodDqptK4mcHd3x9/+9jfEx8cjODgYr776KuLi4lBQUIDevXtL75OdnY3Q0FDExcUhLi4ODz74IDZv3mxTWx3eY3D48GG0tLRIr0+ePIknnngCv/rVr2Svq6iosBhPGThwoKObRkRE1C2PXY6NjYWQuUjuHHBjWWRRUVGH7+Pj44MtW7bY1rhbODwxuPULfenSpRgyZAjGjBkje52fn1+b9ZpERETUtZw6+bCpqQlbtmzBiy++2OFSifDwcOh0OowbNw4HDhzosO7GxkaYTCaLQkRE1JGeuiqhp3BqYrBjxw7U19dj9uzZ7cbodDp8+umn2LZtG7Zv347g4GCMGzcOxcXFsnVnZGRY7DYVEBDg4NYTEZFLap0jYG9xUSrR0cCGHeLj4+Hu7o4vvvjCpusmTpwIlUqFXbt2tRvT2NiIxsZG6bXJZEJAQABiVZPRR+V2220mIqKud100o1DsgNFodOiGQTczmUzQaDS4J/Mt9OrnYVdd5p+v4fs57zq1vd3FacsVz507h4KCAmzfvt3ma6OiojqcPNHegyqIiIjkdMfkwzuJ0xKDDRs2wM/PD0899ZTN1x4/flzasIGIiMihumEfgzuJUxIDs9mMDRs2ICkpCX36WL7FkiVLcOHCBWzatAkAsGLFCtxzzz0YPny4NFlx27Zt2LZtmzOaRkRERDKckhgUFBSguroaL774YptztbW1qK6ull43NTUhLS0NFy5cQN++fTF8+HDs3r0bTz75pDOaRkRECueIVQWuvCrBKYlBXFxcu5s1ZGVlWbxetGgRFi1a5IxmEBERWefCQwH24kOUiIiISMKHKBERkaJwKEEeEwMiIlIWrkqQxcSAiIgURvXvYm8drolzDIiIiEjCHgMiIlIWDiXIYmJARETKwsRAFocSiIiISMIeAyIiUhZHPDaZyxWJiIhcA5+uKI9DCURERCRhjwERESkLJx/KYmJARETKwjkGsjiUQERERBL2GBARkaKoxI1ibx2uiokBEREpC+cYyGJiQEREysI5BrI4x4CIiIgk7DEgIiJl4VCCLCYGRESkLEwMZHEogYiIiCTsMSAiImVhj4EsJgZERKQsXJUgi0MJREREJGGPARERKQp3PpTHxICIiJSFcwxkcSiBiIjIyYqLizFx4kTo9XqoVCrs2LHD4rwQAunp6dDr9ejbty9iY2Nx6tQpi5jGxkYsWLAAvr6+8PT0xKRJk3D+/HmLmMuXLyMxMREajQYajQaJiYmor6+3qa1MDIiIiJzs6tWrCAsLw6pVq6ye/+CDD/Dhhx9i1apVOHz4MLRaLZ544glcuXJFiklJSUFubi5ycnJQUlKChoYGTJgwAS0tLVLMzJkzUV5ejry8POTl5aG8vByJiYk2tZVDCUREpCgqOGCOwb//32QyWRxXq9VQq9Vt4hMSEpCQkGC1LiEEVqxYgTfffBNTpkwBAGzcuBH+/v747LPP8NJLL8FoNCIzMxObN2/G448/DgDYsmULAgICUFBQgPj4eJw5cwZ5eXk4ePAgRo0aBQBYt24doqOjUVFRgeDg4E59NvYYEBGRsrQuV7S3AAgICJC67TUaDTIyMmxuTlVVFQwGA+Li4qRjarUaY8aMQWlpKQDg6NGjaG5utojR6/UICQmRYsrKyqDRaKSkAACioqKg0WikmM5gjwEREdFtqqmpgbe3t/TaWm9BRwwGAwDA39/f4ri/vz/OnTsnxbi7u6N///5tYlqvNxgM8PPza1O/n5+fFNMZTAyIiEhZHLgqwdvb2yIxsIdKZblpkhCizbE2zbglxlp8Z+q5GYcSiIhIWYSDioNotVoAaPNbfV1dndSLoNVq0dTUhMuXL8vGXLx4sU39P/74Y5veCDlMDIiIiLpRUFAQtFot8vPzpWNNTU0oKipCTEwMACAiIgJubm4WMbW1tTh58qQUEx0dDaPRiK+//lqKOXToEIxGoxTTGRxKICIiRemOnQ8bGhpw9uxZ6XVVVRXKy8vh4+ODwYMHIyUlBe+99x6GDh2KoUOH4r333kO/fv0wc+ZMAIBGo8GcOXOwcOFCDBgwAD4+PkhLS0NoaKi0SuH+++/H+PHjkZycjLVr1wIA5s6diwkTJnR6RQJwGz0GjtikwZpt27bhgQcegFqtxgMPPIDc3Fxbm0ZERNSxbhhKOHLkCMLDwxEeHg4ASE1NRXh4ON5++20AwKJFi5CSkoKXX34ZkZGRuHDhAvbt2wcvLy+pjuXLl2Py5MmYNm0aRo8ejX79+uGLL75A7969pZjs7GyEhoYiLi4OcXFxePDBB7F582ab2qoSQtj08f7617/iq6++wogRIzB16lTk5uZi8uTJ0vn3338f//M//4OsrCwMGzYMv//971FcXIyKigqLD3izsrIyPPLII3j33XfxzDPPIDc3F2+//TZKSkosll3IMZlM0Gg0iFVNRh+Vmy0fiYiIutl10YxCsQNGo9Fhk/lu1fo9cc/v/we9PDzsqst87Rq+/+2bTm1vd7E5MbC4WKWySAyEENDr9UhJScHrr78O4MYWjv7+/nj//ffx0ksvWa1n+vTpMJlM+Otf/yodGz9+PPr374/PP/+8U21hYkBEdOfq0sTgXQclBm+5ZmLg0MmHndmkwZqysjKLawAgPj5e9prGxkaYTCaLQkRE1JHWOQb2Flfl0MRAbpMGuc0VDAaDzddkZGRY7DYVEBBgR8uJiIgIcNJyxdvZpMHWa5YsWQKj0SiVmpqa228wEREphwO3RHZFDl2uePMmDTqdTjp+8wYM7V0nt7GDNe09qIKIiEiWA3c+dEUO7THozCYN1kRHR1tcAwD79u2zaUMGIiKizuAcA3k29xjYu0kDAMyaNQt333239BSq3/zmN3j00Ufx/vvv4+mnn8bOnTtRUFCAkpISB3xEIiIi6iybE4MjR45g7Nix0uvU1FQAQFJSErKysrBo0SL861//wssvv4zLly9j1KhRbTZpqK6uRq9e/+msiImJQU5ODn7729/irbfewpAhQ7B169ZO72FARETUaRxKkGXXPgY9CfcxICK6c3XlPgb3vvUeetu5j0HLtWv47t03uI8BERERuTbXe4iSqteNYg9hdkxbiIgczd5/33qsXl3XPc+hBFmulxgQERHJYWIgy1VTTyIiIroN7DEgIiJFccQ+BK68jwF7DIiIiEjCxICIiIgkHEogIiJl4eRDWUwMiIhIUTjHQB4TAyIiUh4X/mK3F+cYEBERkYQ9BkREpCycYyCLiQERESkK5xjI41ACERERSdhjQEREysKhBFlMDIiISFE4lCCPQwlEREQkYY8BEREpC4cSZDExICIiZWFiIItDCURERCRhj4E1KuZLBECYu7sFdwb+90J3GE4+lMfEgIiIlIVDCbKYGBARkbIwMZDFPkAiIiKSsMeAiIgUhXMM5DExICIiZeFQgiwOJRARETnRPffcA5VK1abMnz8fADB79uw256KioizqaGxsxIIFC+Dr6wtPT09MmjQJ58+fd0p7mRgQEZGitA4l2Fs66/Dhw6itrZVKfn4+AOBXv/qVFDN+/HiLmD179ljUkZKSgtzcXOTk5KCkpAQNDQ2YMGECWlpaHHJPbsahBCIiUpYuHkoYOHCgxeulS5diyJAhGDNmjHRMrVZDq9Vavd5oNCIzMxObN2/G448/DgDYsmULAgICUFBQgPj4eNvbL4M9BkRERLfJZDJZlMbGRtn4pqYmbNmyBS+++CJUKpV0vLCwEH5+fhg2bBiSk5NRV1cnnTt69Ciam5sRFxcnHdPr9QgJCUFpaanDPxMTAyIiUhbhoAIgICAAGo1GKhkZGbJvvWPHDtTX12P27NnSsYSEBGRnZ2P//v1YtmwZDh8+jMcee0xKMgwGA9zd3dG/f3+Luvz9/WEwGOy5E1ZxKIGIiBRF9e9ibx0AUFNTA29vb+m4Wq2WvS4zMxMJCQnQ6/XSsenTp0s/h4SEIDIyEoGBgdi9ezemTJnSbl1CCIteB0dhYkBERHSbvL29LRIDOefOnUNBQQG2b98uG6fT6RAYGIjKykoAgFarRVNTEy5fvmzRa1BXV4eYmJjbb3w7OJRARETK4sChBFts2LABfn5+eOqpp2TjLl26hJqaGuh0OgBAREQE3NzcpNUMAFBbW4uTJ086JTFgjwERESlKd+x8aDabsWHDBiQlJaFPn/989TY0NCA9PR1Tp06FTqfD999/jzfeeAO+vr545plnAAAajQZz5szBwoULMWDAAPj4+CAtLQ2hoaHSKgVHsrnHoLi4GBMnToRer4dKpcKOHTukc83NzXj99dcRGhoKT09P6PV6zJo1Cz/88INsnVlZWVY3f7h27ZrNH4iIiEhWN/QYFBQUoLq6Gi+++KLF8d69e+PEiRN4+umnMWzYMCQlJWHYsGEoKyuDl5eXFLd8+XJMnjwZ06ZNw+jRo9GvXz988cUX6N27923cAHk29xhcvXoVYWFheOGFFzB16lSLcz///DOOHTuGt956C2FhYbh8+TJSUlIwadIkHDlyRLZeb29vVFRUWBzz8PCwtXlEREQ9TlxcHIRom0307dsXe/fu7fB6Dw8PrFy5EitXrnRG8yzYnBgkJCQgISHB6jmNRmMxBgIAK1euxEMPPYTq6moMHjy43XpVKlW7mzsQERE5lAs/68BeTp9jYDQaoVKpcNddd8nGNTQ0IDAwEC0tLfjlL3+Jd999F+Hh4e3GNzY2WmwkYTKZHNVkohtUnJtL5Ir4dEV5Tv2X79q1a1i8eDFmzpwpu5zjvvvuQ1ZWFnbt2oXPP/8cHh4eGD16tLRUw5qMjAyLTSUCAgKc8RGIiIgUxWmJQXNzM2bMmAGz2YzVq1fLxkZFReH//b//h7CwMDzyyCP43//9XwwbNkx2LGXJkiUwGo1SqampcfRHICIiV9RNyxXvFE4ZSmhubsa0adNQVVWF/fv3d3rzh1a9evXCyJEjZXsM1Gp1hztMERER3YpDCfIc3mPQmhRUVlaioKAAAwYMsLkOIQTKy8ulzR2IiIioa9jcY9DQ0ICzZ89Kr6uqqlBeXg4fHx/o9Xo8++yzOHbsGP7yl7+gpaVFesCDj48P3N3dAQCzZs3C3XffLT1s4ne/+x2ioqIwdOhQmEwmfPTRRygvL8fHH3/siM9IRET0H1382OU7jc2JwZEjRzB27FjpdWpqKgAgKSkJ6enp2LVrFwDgl7/8pcV1Bw4cQGxsLACguroavXr9p7Oivr4ec+fOhcFggEajQXh4OIqLi/HQQw/Z2jwiIiJZHEqQZ3NiEBsba3WThlZy51oVFhZavF6+fDmWL19ua1OIiIjIwfisBCIiUhYOJchiYkBERMrCxEAWEwMiIlIUzjGQxz1fiYiISMIeAyIiUhYOJchiYkBERIqiEgKqTqyg66gOV8WhBCIiIpKwx4CIiJSFQwmymBgQEZGicFWCPA4lEBERkYQ9BkREpCwcSpDFxICIiBSFQwnyOJRAREREEvYYEBGRsnAoQRYTAyIiUhQOJchjYkBERMrCHgNZnGNAREREEvYYEBGR4rjyUIC9mBgQEZGyCHGj2FuHi+JQAhEREUnYY0BERIrCVQnymBgQEZGycFWCLA4lEBERkYQ9BkREpCgq841ibx2uiokBEREpC4cSZHEogYiIiCRMDIiISFFaVyXYWzorPT0dKpXKomi1Wum8EALp6enQ6/Xo27cvYmNjcerUKYs6GhsbsWDBAvj6+sLT0xOTJk3C+fPnHXVLLDAxICIiZWnd4MjeYoPhw4ejtrZWKidOnJDOffDBB/jwww+xatUqHD58GFqtFk888QSuXLkixaSkpCA3Nxc5OTkoKSlBQ0MDJkyYgJaWFofdllacY0BERIrSHfsY9OnTx6KXoJUQAitWrMCbb76JKVOmAAA2btwIf39/fPbZZ3jppZdgNBqRmZmJzZs34/HHHwcAbNmyBQEBASgoKEB8fLx9H+YW7DEgIiK6TSaTyaI0NjZajausrIRer0dQUBBmzJiB7777DgBQVVUFg8GAuLg4KVatVmPMmDEoLS0FABw9ehTNzc0WMXq9HiEhIVKMIzExICIiZREOKgACAgKg0WikkpGR0ebtRo0ahU2bNmHv3r1Yt24dDAYDYmJicOnSJRgMBgCAv7+/xTX+/v7SOYPBAHd3d/Tv37/dGEfiUAIRESmKI4cSampq4O3tLR1Xq9VtYhMSEqSfQ0NDER0djSFDhmDjxo2Iioq6UZ9KZXGNEKLNsVt1JuZ2sMeAiIjoNnl7e1sUa4nBrTw9PREaGorKykpp3sGtv/nX1dVJvQharRZNTU24fPlyuzGOxMSAiIiUpRtWJdyssbERZ86cgU6nQ1BQELRaLfLz86XzTU1NKCoqQkxMDAAgIiICbm5uFjG1tbU4efKkFONIHEogIiJF6epVCWlpaZg4cSIGDx6Muro6/P73v4fJZEJSUhJUKhVSUlLw3nvvYejQoRg6dCjee+899OvXDzNnzgQAaDQazJkzBwsXLsSAAQPg4+ODtLQ0hIaGSqsUHMnmHoPi4mJMnDgRer0eKpUKO3bssDg/e/bsNhs5tI6hyNm2bRseeOABqNVqPPDAA8jNzbW1aURERD3O+fPn8dxzzyE4OBhTpkyBu7s7Dh48iMDAQADAokWLkJKSgpdffhmRkZG4cOEC9u3bBy8vL6mO5cuXY/LkyZg2bRpGjx6Nfv364YsvvkDv3r0d3l6bewyuXr2KsLAwvPDCC5g6darVmPHjx2PDhg3Sa3d3d9k6y8rKMH36dLz77rt45plnkJubi2nTpqGkpASjRo2ytYlERETt6+JnJeTk5MieV6lUSE9PR3p6ersxHh4eWLlyJVauXNn5N75NNicGCQkJFjMsrVGr1VY3cmjPihUr8MQTT2DJkiUAgCVLlqCoqAgrVqzA559/bmsTiYiI2tUdGxzdSZwy+bCwsBB+fn4YNmwYkpOTUVdXJxtfVlZmsXEDAMTHx8tu3NDY2NhmYwkiIiKyj8MTg4SEBGRnZ2P//v1YtmwZDh8+jMcee6zd3aCAG8s05DZ3sCYjI8NiU4mAgACHfQYiInJhZuGY4qIcviph+vTp0s8hISGIjIxEYGAgdu/eLe0DbY2tmzssWbIEqamp0muTycTkgIiIOtbFcwzuNE5frqjT6RAYGIjKysp2Y7RarezmDtao1epObSRBRER0MxUcMMfAIS3pmZy+wdGlS5dQU1MDnU7Xbkx0dLTFxg0AsG/fPqds3EBERETts7nHoKGhAWfPnpVeV1VVoby8HD4+PvDx8UF6ejqmTp0KnU6H77//Hm+88QZ8fX3xzDPPSNfMmjULd999t/Swid/85jd49NFH8f777+Ppp5/Gzp07UVBQgJKSEgd8RCIiopvYuXOhVIeLsjkxOHLkCMaOHSu9bh3nT0pKwpo1a3DixAls2rQJ9fX10Ol0GDt2LLZu3WqxUUN1dTV69fpPZ0VMTAxycnLw29/+Fm+99RaGDBmCrVu3cg8DIiJyOC5XlGdzYhAbGwshkynt3bu3wzoKCwvbHHv22Wfx7LPP2tocIiIiciA+K4GIiJSFqxJkMTEgIiJFUQkBlZ1zBOy9vifjY5eJiIhIwh4DIiJSFvO/i711uCgmBkREpCgcSpDHoQQiIiKSsMeAiIiUhasSZDExICIiZeHOh7KYGBARkaJw50N5nGNAREREEvYYEBGRsnAoQRYTAyIiUhSV+Uaxtw5XxaEEIiIikrDHgIiIlIVDCbKYGBARkbJwHwNZHEogIiIiCXsMiIhIUfisBHlMDIiISFk4x0AWhxKIiIhIwh4DIiJSFgHA3n0IXLfDgIkBEREpC+cYyGNiQEREyiLggDkGDmlJj8Q5BkRERCRhjwERESkLVyXIYmJARETKYgagckAdLopDCURERCRhYkBERIrSuirB3tJZGRkZGDlyJLy8vODn54fJkyejoqLCImb27NlQqVQWJSoqyiKmsbERCxYsgK+vLzw9PTFp0iScP3/eIffkZkwMiIhIWVrnGNhbOqmoqAjz58/HwYMHkZ+fj+vXryMuLg5Xr161iBs/fjxqa2ulsmfPHovzKSkpyM3NRU5ODkpKStDQ0IAJEyagpaXFIbelFecYEBEROVFeXp7F6w0bNsDPzw9Hjx7Fo48+Kh1Xq9XQarVW6zAajcjMzMTmzZvx+OOPAwC2bNmCgIAAFBQUID4+3mHtZY8BEREpiwN7DEwmk0VpbGzs8O2NRiMAwMfHx+J4YWEh/Pz8MGzYMCQnJ6Ourk46d/ToUTQ3NyMuLk46ptfrERISgtLSUkfcFQkTAyIiUhYHJgYBAQHQaDRSycjI6OCtBVJTU/Hwww8jJCREOp6QkIDs7Gzs378fy5Ytw+HDh/HYY49JiYbBYIC7uzv69+9vUZ+/vz8MBoNDbw+HEoiIiG5TTU0NvL29pddqtVo2/pVXXsHf//53lJSUWByfPn269HNISAgiIyMRGBiI3bt3Y8qUKe3WJ4SASmXv2ktL7DEgIiJlMTuoAPD29rYoconBggULsGvXLhw4cACDBg2SbaJOp0NgYCAqKysBAFqtFk1NTbh8+bJFXF1dHfz9/W36+B1hYkBERIrS1csVhRB45ZVXsH37duzfvx9BQUEdXnPp0iXU1NRAp9MBACIiIuDm5ob8/Hwppra2FidPnkRMTIztN0EGhxKIiEhZunhL5Pnz5+Ozzz7Dzp074eXlJc0J0Gg06Nu3LxoaGpCeno6pU6dCp9Ph+++/xxtvvAFfX18888wzUuycOXOwcOFCDBgwAD4+PkhLS0NoaKi0SsFRmBgQERE50Zo1awAAsbGxFsc3bNiA2bNno3fv3jhx4gQ2bdqE+vp66HQ6jB07Flu3boWXl5cUv3z5cvTp0wfTpk3Dv/71L4wbNw5ZWVno3bu3Q9tr81BCcXExJk6cCL1eD5VKhR07dlicv3Xnptbyhz/8od06s7KyrF5z7do1mz8QERGRLLNwTOkkIYTVMnv2bABA3759sXfvXtTV1aGpqQnnzp1DVlYWAgICLOrx8PDAypUrcenSJfz888/44osv2sQ4gs09BlevXkVYWBheeOEFTJ06tc352tpai9d//etfMWfOHKuxN/P29m6zRaSHh4etzSMiIpLHpyvKsjkxSEhIQEJCQrvnb921aefOnRg7dizuvfde2XpVKlW7Oz5Z09jYaLGRhMlk6vS1REREZJ1TVyVcvHgRu3fvxpw5czqMbWhoQGBgIAYNGoQJEybg+PHjsvEZGRkWm0o4ozuFiIhckSM2N3LdHgOnJgYbN26El5eX7OYMAHDfffchKysLu3btwueffw4PDw+MHj1aWr9pzZIlS2A0GqVSU1Pj6OYTEZEr6uKHKN1pnLoqYf369Xj++ec7nCsQFRVl8XjJ0aNHY8SIEVi5ciU++ugjq9eo1eoOd5giIiIi2zgtMfjyyy9RUVGBrVu32nxtr169MHLkSNkeAyIiottidsBQgA2rEu40ThtKyMzMREREBMLCwmy+VgiB8vJyaccnIiIihxFmxxQXZXOPQUNDA86ePSu9rqqqQnl5OXx8fDB48GAAN1YI/PnPf8ayZcus1jFr1izcfffd0lOofve73yEqKgpDhw6FyWTCRx99hPLycnz88ce385mIiIjoNtmcGBw5cgRjx46VXqempgIAkpKSkJWVBQDIycmBEALPPfec1Tqqq6vRq9d/Oivq6+sxd+5cGAwGaDQahIeHo7i4GA899JCtzSMiIpLHfQxkqYRwjU9nMpmg0WgQ22sK+qjcurs5RERkg+uiGYXm7TAajRaPMXak1u+Jx++ehz697Ju8ft3ciIILnzi1vd2Fz0ogIiJlYY+BLD52mYiIiCTsMSAiImURcECPgUNa0iMxMSAiImXhUIIsDiUQERGRhD0GRESkLGYzADs3KDJzgyMiIiLXwKEEWRxKICIiIgl7DIiISFnYYyCLiQERESkLn64oi0MJREREJGGPARERKYoQZgg7H5ts7/U9GRMDIiJSFiHsHwrgHAMiIiIXIRwwx8CFEwPOMSAiIiIJewyIiEhZzGZAZeccAc4xICIichEcSpDFoQQiIiKSsMeAiIgURZjNEHYOJXC5IhERkavgUIIsDiUQERGRhD0GRESkLGYBqNhj0B4mBkREpCxCALB3uaLrJgYcSiAiIiIJewyIiEhRhFlA2DmUIFy4x4CJARERKYsww/6hBNddrsihBCIiUhRhFg4ptlq9ejWCgoLg4eGBiIgIfPnll074dPZjYkBERORkW7duRUpKCt58800cP34cjzzyCBISElBdXd3dTWvDZYYSWsd7rovmbm4JERHZqvXf7q4Yu78uGu0eCriOG+01mUwWx9VqNdRqdZv4Dz/8EHPmzMGvf/1rAMCKFSuwd+9erFmzBhkZGXa1xdFcJjG4cuUKAKBEfGH3hlZERNQ9rly5Ao1G45S63d3dodVqUWLY45D6fvGLXyAgIMDi2DvvvIP09HSLY01NTTh69CgWL15scTwuLg6lpaUOaYsjuUxioNfrUVNTAy8vL6hUKqsxJpMJAQEBqKmpgbe3dxe38Pax3V3vTm0729212G7HEULgypUr0Ov1TnsPDw8PVFVVoampySH1CSHafN9Y6y346aef0NLSAn9/f4vj/v7+MBgMDmmLI7lMYtCrVy8MGjSoU7He3t495j8GW7DdXe9ObTvb3bXYbsdwVk/BzTw8PODh4eH097Hm1iTCWmLRE3DyIRERkRP5+vqid+/ebXoH6urq2vQi9ARMDIiIiJzI3d0dERERyM/Ptzien5+PmJiYbmpV+1xmKKEz1Go13nnnHatjQD0Z29317tS2s91di+2mzkpNTUViYiIiIyMRHR2NTz/9FNXV1Zg3b153N60NlXDlfR2JiIh6iNWrV+ODDz5AbW0tQkJCsHz5cjz66KPd3aw2mBgQERGRhHMMiIiISMLEgIiIiCRMDIiIiEjCxICIiIgkLpcY2PpYy6KiIkRERMDDwwP33nsvPvnkky5q6Q0ZGRkYOXIkvLy84Ofnh8mTJ6OiokL2msLCQqhUqjbl22+/7aJWA+np6W3eX6vVyl7T3fe61T333GP1/s2fP99qfHfd7+LiYkycOBF6vR4qlQo7duywOC+EQHp6OvR6Pfr27YvY2FicOnWqw3q3bduGBx54AGq1Gg888AByc3O7rN3Nzc14/fXXERoaCk9PT+j1esyaNQs//PCDbJ1ZWVlW/wyuXbvWJe0GgNmzZ7d5/6ioqA7r7c77DcDqfVOpVPjDH/7Qbp1dcb+p53KpxMDWx1pWVVXhySefxCOPPILjx4/jjTfewKuvvopt27Z1WZuLioowf/58HDx4EPn5+bh+/Tri4uJw9erVDq+tqKhAbW2tVIYOHdoFLf6P4cOHW7z/iRMn2o3tCfe61eHDhy3a3brpyK9+9SvZ67r6fl+9ehVhYWFYtWqV1fMffPABPvzwQ6xatQqHDx+GVqvFE088IT1QzJqysjJMnz4diYmJ+Oabb5CYmIhp06bh0KFDXdLun3/+GceOHcNbb72FY8eOYfv27fjHP/6BSZMmdVivt7e3xf2vra116Na2Hd1vABg/frzF++/ZI/8wnu6+3wDa3LP169dDpVJh6tSpsvU6+35TDyZcyEMPPSTmzZtncey+++4Tixcvthq/aNEicd9991kce+mll0RUVJTT2tiRuro6AUAUFRW1G3PgwAEBQFy+fLnrGnaLd955R4SFhXU6vife61a/+c1vxJAhQ4TZbLZ6vifcbwAiNzdXem02m4VWqxVLly6Vjl27dk1oNBrxySeftFvPtGnTxPjx4y2OxcfHixkzZji8zUK0bbc1X3/9tQAgzp07127Mhg0bhEajcWzjZFhrd1JSknj66adtqqcn3u+nn35aPPbYY7IxXX2/qWdxmR6D1sdaxsXFWRyXe6xlWVlZm/j4+HgcOXIEzc3NTmurHKPRCADw8fHpMDY8PBw6nQ7jxo3DgQMHnN20NiorK6HX6xEUFIQZM2bgu+++aze2J95r4Mbfmy1btuDFF1/s8GEm3X2/b1ZVVQWDwWBxT9VqNcaMGSP7GNf2/hy689GvRqMRKpUKd911l2xcQ0MDAgMDMWjQIEyYMAHHjx/vmgbepLCwEH5+fhg2bBiSk5NRV1cnG9/T7vfFixexe/duzJkzp8PYnnC/qXu4TGJwO4+1NBgMVuOvX7+On376yWltbY8QAqmpqXj44YcREhLSbpxOp8Onn36Kbdu2Yfv27QgODsa4ceNQXFzcZW0dNWoUNm3ahL1792LdunUwGAyIiYnBpUuXrMb3tHvdaseOHaivr8fs2bPbjekJ9/tWrX+nbX2Ma3t/Dt316Ndr165h8eLFmDlzpuxT/u677z5kZWVh165d+Pzzz+Hh4YHRo0ejsrKyy9qakJCA7Oxs7N+/H8uWLcPhw4fx2GOPobGxsd1retr93rhxI7y8vDBlyhTZuJ5wv6n7uNyzEmx9rKW1eGvHu8Irr7yCv//97ygpKZGNCw4ORnBwsPQ6OjoaNTU1+OMf/9hl22smJCRIP4eGhiI6OhpDhgzBxo0bkZqaavWannSvW2VmZiIhIUH2GfA94X6353Ye49pTHv3a3NyMGTNmwGw2Y/Xq1bKxUVFRFhP9Ro8ejREjRmDlypX46KOPnN1UAMD06dOln0NCQhAZGYnAwEDs3r1b9ou2p9xvAFi/fj2ef/75DucK9IT7Td3HZXoMbuexllqt1mp8nz59MGDAAKe11ZoFCxZg165dOHDgAAYNGmTz9VFRUd2azXt6eiI0NLTdNvSke93q3LlzKCgowK9//Wubr+3u+926AsTWx7i29+fQ1Y9+bW5uxrRp01BVVYX8/HzZ3gJrevXqhZEjR3brn4FOp0NgYKBsG3rK/QaAL7/8EhUVFbf1970n3G/qOi6TGNzOYy2jo6PbxO/btw+RkZFwc3NzWltvJoTAK6+8gu3bt2P//v0ICgq6rXqOHz8OnU7n4NZ1XmNjI86cOdNuG3rCvb7Vhg0b4Ofnh6eeesrma7v7fgcFBUGr1Vrc06amJhQVFck+xrW9P4eufPRra1JQWVmJgoKC20oMhRAoLy/v1j+DS5cuoaamRrYNPeF+t8rMzERERATCwsJsvrYn3G/qQt0169EZcnJyhJubm8jMzBSnT58WKSkpwtPTU3z//fdCCCEWL14sEhMTpfjvvvtO9OvXT7z22mvi9OnTIjMzU7i5uYn/+7//67I2/9d//ZfQaDSisLBQ1NbWSuXnn3+WYm5t9/Lly0Vubq74xz/+IU6ePCkWL14sAIht27Z1WbsXLlwoCgsLxXfffScOHjwoJkyYILy8vHr0vb5ZS0uLGDx4sHj99dfbnOsp9/vKlSvi+PHj4vjx4wKA+PDDD8Xx48el2ftLly4VGo1GbN++XZw4cUI899xzQqfTCZPJJNWRmJhosSrnq6++Er179xZLly4VZ86cEUuXLhV9+vQRBw8e7JJ2Nzc3i0mTJolBgwaJ8vJyi7/zjY2N7bY7PT1d5OXliX/+85/i+PHj4oUXXhB9+vQRhw4d6pJ2X7lyRSxcuFCUlpaKqqoqceDAAREdHS3uvvvuHn2/WxmNRtGvXz+xZs0aq3V0x/2mnsulEgMhhPj4449FYGCgcHd3FyNGjLBY9peUlCTGjBljEV9YWCjCw8OFu7u7uOeee9r9D8dZAFgtGzZsaLfd77//vhgyZIjw8PAQ/fv3Fw8//LDYvXt3l7Z7+vTpQqfTCTc3N6HX68WUKVPEqVOn2m2zEN1/r2+2d+9eAUBUVFS0OddT7nfrMslbS1JSkhDixpLFd955R2i1WqFWq8Wjjz4qTpw4YVHHmDFjpPhWf/7zn0VwcLBwc3MT9913n8MTHLl2V1VVtft3/sCBA+22OyUlRQwePFi4u7uLgQMHiri4OFFaWtpl7f75559FXFycGDhwoHBzcxODBw8WSUlJorq62qKOnna/W61du1b07dtX1NfXW62jO+439Vx87DIRERFJXGaOAREREdmPiQERERFJmBgQERGRhIkBERERSZgYEBERkYSJAREREUmYGBAREZGEiQERERFJmBgQERGRhIkBERERSZgYEBERkeT/AygWtIhHw0BIAAAAAElFTkSuQmCC",
      "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_closed_boundary(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": 127,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "77835379be7343d0aff52d275a86be97",
       "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": 127,
     "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, vmin=0, vmax=50)\n",
    "    \n",
    "interact(update, w = IntSlider(min=0, max = 999, step = 1, value = 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD4CAYAAADl7fPiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAREUlEQVR4nO3dfYxc1X3G8e/jtaGVMbGJwQHjEppaSFZU3MhyGpVWpiTURihOqjS1VbVOS7UkClIjNapoK0GU/kNVUaTECOIkFqRKgL45sRQLsGglgpQXFmRenEBxLUfs4thNTHFwiGDtX//Yu9Gc9cz63rn37NyZPh9pNXfuPXPuuTvLw7z8fI4iAjOzWYsGPQAzaxeHgpklHApmlnAomFnCoWBmicWDHkA3YxcujSUXLx/0MErL8wWOKgygQre5vmyq1G+ma6siR7+Zxlrht1XaW6+e4PSpU127bmUoLLl4OWvu+Hjj/Vb6jzfKPxVnSraNM+X7rNT2dIU/m1xtK/y+VGkM5ZtW6VdnSjas8DyU7hMqBUiOUJj83F09j/ntg5klaoWCpM2SXpR0SNKtXY6fL+mh4vh3Jb2zzvnMLL++Q0HSGHA3sAVYB2yXtG5Os5uAVyPi14C7gL/v93xmtjDqvFLYCByKiMMR8SbwILB1TputwP3F9r8C10nK8RbJzBpSJxRWAy933J8s9nVtExHTwGvA27t1Jmlc0oSkidMnT9UYlpnV0ZoPGiNiV0RsiIgNYxcuHfRwzP7fqhMKU8CajvuXF/u6tpG0GHgb8JMa5zSzzOqEwpPAWklXSjoP2AbsndNmL7Cj2P4I8B/hf6tt1mp9Fy9FxLSkW4BHgDFgd0QclPRZYCIi9gJfBv5J0iHgBDPBYWYtVquiMSL2Afvm7LutY/vnwB/UOce5x1ClcYUqwQptS/fbhnLkXEb52oZJA7/b1nzQaGbt4FAws4RDwcwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLNEKyduhQrly5lKlwdeiZujzDpr2/JNK8k1hpJtlal8u9JMQwv8x+hXCmaWcCiYWcKhYGYJh4KZJRwKZpZwKJhZwqFgZok6K0StkfSfkr4v6aCkv+jSZpOk1yQdKH5u69aXmbVHneKlaeAvI+JpScuApyTtj4jvz2n3rYi4scZ5zGwB9f1KISKORsTTxfZPgR9w9gpRZjZkGilzLlaT/g3gu10Ov0/SM8ArwKcj4mCPPsaBcYDFK99GnGn+445K1aJVSqLPlG1XpSS7dFMoef6Zjiu0rUBVypErjLdKmXG1kuTmZ+DOVrqc63fQQ+3/8iRdAPwb8KmIODnn8NPAFRFxNfB54Ou9+ulcNm7RMi8bZzYotUJB0hJmAuGrEfHvc49HxMmIeL3Y3gcskbSyzjnNLK863z6ImRWgfhAR/9ijzTtml56XtLE4n9eSNGuxOp8p/Bbwx8Bzkg4U+/4G+BWAiLiXmfUjPyFpGngD2Oa1JM3arc5akk9wjs9WImInsLPfc5jZwnNFo5klHApmlnAomFnCoWBmCYeCmSXaO5tz6YaZSocrtC1dvlypFDjTrMtV+s1VPp2pJLpSOXDJfrPN5jzo0vR5+vQrBTNLOBTMLOFQMLOEQ8HMEg4FM0s4FMws4VAws4RDwcwSDgUzS7S3ovF0pWkwy/VZpZIuR0VjpWrCNlT9ZZqMNVvbDM9vG6oUh23iVjMbLQ4FM0s0McX7EUnPFcvCTXQ5Lkmfk3RI0rOS3lP3nGaWT1OfKVwbET/ucWwLsLb4eS9wT3FrZi20EG8ftgJfiRnfAZZLunQBzmtmfWgiFAJ4VNJTxdJvc60GXu64P0mXNScljUuakDRx+uSpBoZlZv1o4u3DNRExJekSYL+kFyLi8aqdRMQuYBfA+b+62mtDmA1I7VcKETFV3B4H9gAb5zSZAtZ03L+82GdmLVR3LcmlkpbNbgPXA8/PabYX+JPiW4jfBF6LiKN1zmtm+dR9+7AK2FMsF7kY+FpEPCzp4/CLpeP2ATcAh4CfAX9a85xmllGtUIiIw8DVXfbf27EdwCerdaw8k6FWGUKOktlcpcsV+s1WulyhLD1b6XKG8ulcpcuVrqsClzmbWeMcCmaWcCiYWcKhYGYJh4KZJRwKZpZwKJhZwqFgZgmHgpklHApmlmjxbM4Z8irTDLqly4xzlE5TsWT2dJ5+K5UuV5mpe8CzRA98NumqbRvo068UzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEn2HgqSriqXiZn9OSvrUnDabJL3W0ea22iM2s6z6Ll6KiBeB9QCSxpiZtn1Pl6bfiogb+z2PmS2spt4+XAf8d0T8sKH+zGxAmipz3gY80OPY+yQ9A7wCfDoiDnZrVCw5Nw4wdtFyqFIKW1au0tKyM0q3YhbjCqXLVUqih6h0uUrbbM9Dtr/FCm17aGIp+vOADwL/0uXw08AVEXE18Hng6736iYhdEbEhIjaMLVtad1hm1qcm3j5sAZ6OiGNzD0TEyYh4vdjeByyRtLKBc5pZJk2EwnZ6vHWQ9A4Vy0dJ2lic7ycNnNPMMqn1mUKxfuQHgJs79nUuGfcR4BOSpoE3gG3FilFm1lJ1l407Bbx9zr7OJeN2AjvrnMPMFpYrGs0s4VAws4RDwcwSDgUzSzgUzCzRztmcgzxlzlXHUFLpUtgKZbBVZn7OVgo8bKXLFcabYzbnVpQ5N3B+v1Iws4RDwcwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEu0sc0bVSmzLGvAMuvlmaK7SNk/5dKXxVpklOtfs0znKnKuUxg+4zHk+fqVgZolSoSBpt6Tjkp7v2HeRpP2SXipuV/R47I6izUuSdjQ1cDPLo+wrhfuAzXP23Qo8FhFrgceK+wlJFwG3A+8FNgK39woPM2uHUqEQEY8DJ+bs3grcX2zfD3yoy0N/D9gfESci4lVgP2eHi5m1SJ3PFFZFxNFi+0fAqi5tVgMvd9yfLPaZWUs18kFjsZZDrc9IJY1LmpA0cfr115sYlpn1oU4oHJN0KUBxe7xLmylgTcf9y4t9Z0nWkrzgghrDMrM66oTCXmD224QdwDe6tHkEuF7SiuIDxuuLfWbWUmW/knwA+DZwlaRJSTcBdwAfkPQS8P7iPpI2SPoSQEScAP4OeLL4+Wyxz8xaqlRFY0Rs73Houi5tJ4A/77i/G9jd1+jMbMG1s8w5gLIlq5Fn1uccZajVyoYzzeZcpRS3QtlwlTLnRblmic5R5pypHDlX+XQTfbrM2cwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEg4FM0u0tsxZ082XL2ebQbdkqXWl8+cqXc40Q/MwlS5XaZttNudMz28T/ErBzBIOBTNLOBTMLOFQMLOEQ8HMEg4FM0s4FMwscc5Q6LGO5D9IekHSs5L2SFre47FHJD0n6YCkiQbHbWaZlHmlcB9nL/W2H3h3RPw68F/AX8/z+GsjYn1EbOhviGa2kM4ZCt3WkYyIRyNiurj7HWYWeTGzEdBEmfOfAQ/1OBbAo5IC+EJE7OrViaRxYBxg8fIVLMpQ5lxJjtmcq5w/Q8lu9bajWbpcpW2lsbagjL303+I87WqFgqS/BaaBr/Zock1ETEm6BNgv6YXilcdZisDYBXD+mjULXO1tZrP6/vZB0seAG4E/KhaYPUtETBW3x4E9wMZ+z2dmC6OvUJC0Gfgr4IMR8bMebZZKWja7zcw6ks93a2tm7VHmK8lu60juBJYx85bggKR7i7aXSdpXPHQV8ISkZ4DvAd+MiIezXIWZNeacnyn0WEfyyz3avgLcUGwfBq6uNTozW3CuaDSzhEPBzBIOBTNLOBTMLOFQMLNEK2dzVlQrLx24BkpLz5JtNufypcuVnoM2lGVnKJ8edJk1kKXkfr52fqVgZgmHgpklHApmlnAomFnCoWBmCYeCmSUcCmaWcCiYWcKhYGaJVlY00oaKxgyzRFabgLNC5WGuSroWVPMNum2+SWbL/zFUes4a4FcKZpZwKJhZot9l4z4jaaqYn/GApBt6PHazpBclHZJ0a5MDN7M8+l02DuCuYjm49RGxb+5BSWPA3cAWYB2wXdK6OoM1s/z6WjaupI3AoYg4HBFvAg8CW/vox8wWUJ3PFG4pVp3eLWlFl+OrgZc77k8W+7qSNC5pQtLE6VOnagzLzOroNxTuAd4FrAeOAnfWHUhE7IqIDRGxYWzp0rrdmVmf+gqFiDgWEacj4gzwRbovBzcFrOm4f3mxz8xarN9l4y7tuPthui8H9ySwVtKVks4DtgF7+zmfmS2cc1Y0FsvGbQJWSpoEbgc2SVrPTN3fEeDmou1lwJci4oaImJZ0C/AIMAbsjoiDOS7CzJqTbdm44v4+4KyvK88pYNFbJct8y1cDD5dKk7FW6DdX2XCm8Q56ktdqfeYpXc4yyasnbjWzshwKZpZwKJhZwqFgZgmHgpklHApmlnAomFnCoWBmCYeCmSUcCmaWaOdszoCmyzas0mn5phUmU86i2szPmfqtIEspbs5+S7bNNetylVmiF2WYUXq+6/crBTNLOBTMLOFQMLOEQ8HMEg4FM0s4FMws4VAws0SZORp3AzcCxyPi3cW+h4CriibLgf+NiPVdHnsE+ClwGpiOiA2NjNrMsilTvHQfsBP4yuyOiPjD2W1JdwKvzfP4ayPix/0O0MwWVpmJWx+X9M5uxyQJ+Cjwuw2Py8wGpG6Z828DxyLipR7HA3hUUgBfiIhdvTqSNA6MAyy5cAVjb5UbwJmxCqOtUro8RJ+2DLoUuDX9ZpjVOtesy1VKlxdNl/8lLCr5zwPmew7qhsJ24IF5jl8TEVOSLgH2S3qhWLD2LEVg7AL45UvXZKrQN7Nz6fv/h5IWA78PPNSrTURMFbfHgT10X17OzFqkzovk9wMvRMRkt4OSlkpaNrsNXE/35eXMrEXOGQrFsnHfBq6SNCnppuLQNua8dZB0maTZFaFWAU9Iegb4HvDNiHi4uaGbWQ79LhtHRHysy75fLBsXEYeBq2uOz8wW2BB9xm5mC8GhYGYJh4KZJRwKZpZwKJhZopWzOesMjL1RsvEvle/3TIWrrVKyWnrm5ypl1m2o6cw1hjZcW1ktKMkuW7oMsPjnJTueZ5Zqv1Iws4RDwcwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEopoX82ppP8Bfjhn90pgFNePGNXrgtG9tlG4risi4uJuB1oZCt1ImhjFFaZG9bpgdK9tVK9rlt8+mFnCoWBmiWEKhZ6rSw25Ub0uGN1rG9XrAoboMwUzWxjD9ErBzBaAQ8HMEkMRCpI2S3pR0iFJtw56PE2RdETSc5IOSJoY9HjqkLRb0nFJz3fsu0jSfkkvFbcrBjnGfvS4rs9ImiqetwOSbhjkGJvW+lCQNAbcDWwB1gHbJa0b7KgadW1ErB+B773vAzbP2Xcr8FhErAUeK+4Pm/s4+7oA7iqet/URsa/L8aHV+lBgZqXqQxFxOCLeBB4Etg54TDZHRDwOnJizeytwf7F9P/ChhRxTE3pc10gbhlBYDbzccX+y2DcKAnhU0lOSxgc9mAxWRcTRYvtHzCw6PCpukfRs8fZi6N4WzWcYQmGUXRMR72HmrdEnJf3OoAeUS8x89z0q33/fA7wLWA8cBe4c6GgaNgyhMAWs6bh/ebFv6EXEVHF7HNjDzFulUXJM0qUAxe3xAY+nERFxLCJOR8QZ4IuM2PM2DKHwJLBW0pWSzgO2AXsHPKbaJC2VtGx2G7geeH7+Rw2dvcCOYnsH8I0BjqUxs0FX+DAj9ry1coWoThExLekW4BFgDNgdEQcHPKwmrAL2SIKZ5+FrEfHwYIfUP0kPAJuAlZImgduBO4B/lnQTM/8U/qODG2F/elzXJknrmXk7dAS4eVDjy8FlzmaWGIa3D2a2gBwKZpZwKJhZwqFgZgmHgpklHApmlnAomFni/wCOWVT0CjszYQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.animation as animation\n",
    "\n",
    "\n",
    "frames = []\n",
    "fig = plt.figure()\n",
    "for i in records:\n",
    "    frames.append([plt.imshow(i, vmin=0, vmax=10, animated=True)])\n",
    "\n",
    "ani = animation.ArtistAnimation(fig, frames, interval=50, blit=True, repeat_delay=1000)\n",
    "\n",
    "ani.save(\"/Users/hannessigner/Desktop/test.gif\")"
   ]
  }
 ],
 "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.9.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}