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": 72,
   "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']+2,grid_size['y']+2))\n",
    "alpha_y = np.ones((grid_size['x']+2,grid_size['y']+2))\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": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAGdCAYAAABKG5eZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhtElEQVR4nO3de2xUZeL/8c/BwoCmHa3QzgyUWkm9cAmL3KtikaVQVhYUpcoulLhRicCKDRHrJeLmG0bc1bAIQnQFJN7YTYGyC66UQFtZCgGl6BJkS6y2C53tQmQGUKYtnN8f/hx37AUHZmyf4f1KTuI58zzHZyYT357OmdaybdsWAACG6NTeCwAAIBKECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBREtp7AdFy/vx5HTt2TImJibIsq72XAwCIgG3bOnXqlDwejzp1avuaKm7CdezYMaWlpbX3MgAAl6C2tla9evVqc0zchCsxMVGSdJsmKEGd23k1AIBINKlRO7Ul9N/ytsRNuL778WCCOivBIlwAYJT//1tzf8xHPdycAQAwSszC9eqrryojI0Ndu3bV4MGD9eGHH7Y5vqysTIMHD1bXrl11/fXXa+XKlbFaGgDAYDEJ17p16zRv3jw9/fTT2r9/v26//Xbl5uaqpqamxfHV1dWaMGGCbr/9du3fv19PPfWUfvvb36qoqCgWywMAGMyKxd/jGj58uG655RatWLEidOzmm2/W5MmT5fV6m41fsGCBNm3apEOHDoWOzZo1SwcOHFBFRcWP+ncGAgE5nU5laxKfcQGAYZrsRpWqWH6/X0lJSW2OjfoVV0NDgz766CPl5OSEHc/JydGuXbtanFNRUdFs/Lhx47Rv3z41Nja2OCcYDCoQCIRtAID4F/VwHT9+XOfOnVNqamrY8dTUVPl8vhbn+Hy+Fsc3NTXp+PHjLc7xer1yOp2hje9wAcDlIWY3Z/zwlkbbttu8zbGl8S0d/05hYaH8fn9oq62tvcQVAwBMEPXvcXXv3l1XXHFFs6ur+vr6ZldV33G5XC2OT0hI0LXXXtviHIfDIYfDEZ1FAwCMEfUrri5dumjw4MEqKSkJO15SUqKsrKwW54wcObLZ+K1bt2rIkCHq3JkbLQAA34vJjwoLCgr0pz/9SatWrdKhQ4f0+OOPq6amRrNmzZL07Y/5ZsyYERo/a9YsffnllyooKNChQ4e0atUqvfHGG5o/f34slgcAMFhMfuVTXl6eTpw4od/97neqq6tT//79tWXLFqWnp0uS6urqwr7TlZGRoS1btujxxx/X8uXL5fF4tHTpUk2ZMiUWywMAGCwm3+NqD3yPCwDM1a7f4wIAIJYIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARol6uLxer4YOHarExESlpKRo8uTJOnz4cJtzSktLZVlWs+2zzz6L9vIAAIaLerjKyso0e/Zs7d69WyUlJWpqalJOTo7OnDlzwbmHDx9WXV1daMvMzIz28gAAhkuI9gn//ve/h+2vXr1aKSkp+uijjzRq1Kg256akpOjqq6+O9pIAAHEk5p9x+f1+SVJycvIFxw4aNEhut1tjxozRjh072hwbDAYVCATCNgBA/ItpuGzbVkFBgW677Tb179+/1XFut1uvvfaaioqKtH79et14440aM2aMysvLW53j9XrldDpDW1paWiyeAgCgg7Fs27ZjdfLZs2dr8+bN2rlzp3r16hXR3IkTJ8qyLG3atKnFx4PBoILBYGg/EAgoLS1N2ZqkBKvzJa0bAPDTarIbVapi+f1+JSUltTk2Zldcc+fO1aZNm7Rjx46IoyVJI0aMUFVVVauPOxwOJSUlhW0AgPgX9ZszbNvW3LlztWHDBpWWliojI+OizrN//3653e4orw4AYLqoh2v27Nl65513VFxcrMTERPl8PkmS0+lUt27dJEmFhYU6evSo1q5dK0lasmSJrrvuOvXr108NDQ166623VFRUpKKiomgvDwBguKiHa8WKFZKk7OzssOOrV6/WzJkzJUl1dXWqqakJPdbQ0KD58+fr6NGj6tatm/r166fNmzdrwoQJ0V4eAMBwMb0546cUCATkdDq5OQMADNQhbs4AACAWCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEaJergWLlwoy7LCNpfL1eacsrIyDR48WF27dtX111+vlStXRntZAIA4kRCLk/br10/btm0L7V9xxRWtjq2urtaECRP00EMP6a233tI//vEPPfroo+rRo4emTJkSi+UBAAwWk3AlJCRc8CrrOytXrlTv3r21ZMkSSdLNN9+sffv26Q9/+APhAgA0E5PPuKqqquTxeJSRkaH7779fn3/+eatjKyoqlJOTE3Zs3Lhx2rdvnxobG1udFwwGFQgEwjYAQPyLeriGDx+utWvX6oMPPtDrr78un8+nrKwsnThxosXxPp9PqampYcdSU1PV1NSk48ePt/rv8Xq9cjqdoS0tLS2qzwMA0DFFPVy5ubmaMmWKBgwYoJ///OfavHmzJOnNN99sdY5lWWH7tm23ePx/FRYWyu/3h7ba2toorB4A0NHF5DOu/3XVVVdpwIABqqqqavFxl8sln88Xdqy+vl4JCQm69tprWz2vw+GQw+GI6loBAB1fzL/HFQwGdejQIbnd7hYfHzlypEpKSsKObd26VUOGDFHnzp1jvTwAgGGiHq758+errKxM1dXV2rNnj+69914FAgHl5+dL+vZHfDNmzAiNnzVrlr788ksVFBTo0KFDWrVqld544w3Nnz8/2ksDAMSBqP+o8N///rceeOABHT9+XD169NCIESO0e/dupaenS5Lq6upUU1MTGp+RkaEtW7bo8ccf1/Lly+XxeLR06VJuhQcAtMiyv7sTwnCBQEBOp1PZmqQEix8xAoBJmuxGlapYfr9fSUlJbY7ldxUCAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMErUw3XdddfJsqxm2+zZs1scX1pa2uL4zz77LNpLAwDEgYRon3Dv3r06d+5caP+f//ynxo4dq/vuu6/NeYcPH1ZSUlJov0ePHtFeGgAgDkQ9XD8MzgsvvKA+ffrojjvuaHNeSkqKrr766mgvBwAQZ2L6GVdDQ4PeeustPfjgg7Isq82xgwYNktvt1pgxY7Rjx44LnjsYDCoQCIRtAID4F9Nwbdy4USdPntTMmTNbHeN2u/Xaa6+pqKhI69ev14033qgxY8aovLy8zXN7vV45nc7QlpaWFuXVAwA6Isu2bTtWJx83bpy6dOmiv/71rxHNmzhxoizL0qZNm1odEwwGFQwGQ/uBQEBpaWnK1iQlWJ0ves0AgJ9ek92oUhXL7/eH3e/Qkqh/xvWdL7/8Utu2bdP69esjnjtixAi99dZbbY5xOBxyOBwXuzwAgKFi9qPC1atXKyUlRb/4xS8inrt//3653e4YrAoAYLqYXHGdP39eq1evVn5+vhISwv8VhYWFOnr0qNauXStJWrJkia677jr169cvdDNHUVGRioqKYrE0AIDhYhKubdu2qaamRg8++GCzx+rq6lRTUxPab2ho0Pz583X06FF169ZN/fr10+bNmzVhwoRYLA0AYLiY3pzxUwoEAnI6ndycAQAGiuTmDH5XIQDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjRByu8vJyTZw4UR6PR5ZlaePGjWGP27athQsXyuPxqFu3bsrOztbBgwcveN6ioiL17dtXDodDffv21YYNGyJdGgDgMhBxuM6cOaOBAwdq2bJlLT7+4osv6uWXX9ayZcu0d+9euVwujR07VqdOnWr1nBUVFcrLy9P06dN14MABTZ8+XVOnTtWePXsiXR4AIM5Ztm3bFz3ZsrRhwwZNnjxZ0rdXWx6PR/PmzdOCBQskScFgUKmpqVq8eLEeeeSRFs+Tl5enQCCg999/P3Rs/Pjxuuaaa/Tuu+/+qLUEAgE5nU5la5ISrM4X+5QAAO2gyW5UqYrl9/uVlJTU5tiofsZVXV0tn8+nnJyc0DGHw6E77rhDu3btanVeRUVF2BxJGjduXJtzgsGgAoFA2AYAiH9RDZfP55Mkpaamhh1PTU0NPdbavEjneL1eOZ3O0JaWlnYJKwcAmCImdxValhW2b9t2s2OXOqewsFB+vz+01dbWXvyCAQDGSIjmyVwul6Rvr6DcbnfoeH19fbMrqh/O++HV1YXmOBwOORyOS1wxAMA0Ub3iysjIkMvlUklJSehYQ0ODysrKlJWV1eq8kSNHhs2RpK1bt7Y5BwBweYr4iuv06dM6cuRIaL+6ulqVlZVKTk5W7969NW/ePC1atEiZmZnKzMzUokWLdOWVV2ratGmhOTNmzFDPnj3l9XolSY899phGjRqlxYsXa9KkSSouLta2bdu0c+fOKDxFAEA8iThc+/bt0+jRo0P7BQUFkqT8/HytWbNGTzzxhL755hs9+uij+uqrrzR8+HBt3bpViYmJoTk1NTXq1On7i72srCy99957euaZZ/Tss8+qT58+WrdunYYPH34pzw0AEIcu6XtcHQnf4wIAc7Xb97gAAIg1wgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFEiDld5ebkmTpwoj8cjy7K0cePG0GONjY1asGCBBgwYoKuuukoej0czZszQsWPH2jznmjVrZFlWs+3s2bMRPyEAQHyLOFxnzpzRwIEDtWzZsmaPff311/r444/17LPP6uOPP9b69ev1r3/9S7/85S8veN6kpCTV1dWFbV27do10eQCAOJcQ6YTc3Fzl5ua2+JjT6VRJSUnYsVdeeUXDhg1TTU2Nevfu3ep5LcuSy+WKdDkAgMtMzD/j8vv9sixLV199dZvjTp8+rfT0dPXq1Ut33XWX9u/f3+b4YDCoQCAQtgEA4l9Mw3X27Fk9+eSTmjZtmpKSklodd9NNN2nNmjXatGmT3n33XXXt2lW33nqrqqqqWp3j9XrldDpDW1paWiyeAgCgg7Fs27YverJlacOGDZo8eXKzxxobG3XfffeppqZGpaWlbYbrh86fP69bbrlFo0aN0tKlS1scEwwGFQwGQ/uBQEBpaWnK1iQlWJ0jfi4AgPbTZDeqVMXy+/0X7EXEn3H9GI2NjZo6daqqq6u1ffv2iKIlSZ06ddLQoUPbvOJyOBxyOByXulQAgGGi/qPC76JVVVWlbdu26dprr434HLZtq7KyUm63O9rLAwAYLuIrrtOnT+vIkSOh/erqalVWVio5OVkej0f33nuvPv74Y/3tb3/TuXPn5PP5JEnJycnq0qWLJGnGjBnq2bOnvF6vJOn555/XiBEjlJmZqUAgoKVLl6qyslLLly+PxnMEAMSRiMO1b98+jR49OrRfUFAgScrPz9fChQu1adMmSdLPfvazsHk7duxQdna2JKmmpkadOn1/sXfy5Ek9/PDD8vl8cjqdGjRokMrLyzVs2LBIlwcAiHOXdHNGRxIIBOR0Ork5AwAMFMnNGfyuQgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGIVwAAKMQLgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGIVwAAKMQLgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGIVwAAKMQLgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGiThc5eXlmjhxojwejyzL0saNG8MenzlzpizLCttGjBhxwfMWFRWpb9++cjgc6tu3rzZs2BDp0gAAl4GIw3XmzBkNHDhQy5Yta3XM+PHjVVdXF9q2bNnS5jkrKiqUl5en6dOn68CBA5o+fbqmTp2qPXv2RLo8AECcS4h0Qm5urnJzc9sc43A45HK5fvQ5lyxZorFjx6qwsFCSVFhYqLKyMi1ZskTvvvtupEsEAMSxmHzGVVpaqpSUFN1www166KGHVF9f3+b4iooK5eTkhB0bN26cdu3a1eqcYDCoQCAQtgEA4l/Uw5Wbm6u3335b27dv10svvaS9e/fqzjvvVDAYbHWOz+dTampq2LHU1FT5fL5W53i9XjmdztCWlpYWtecAAOi4Iv5R4YXk5eWF/rl///4aMmSI0tPTtXnzZt1zzz2tzrMsK2zftu1mx/5XYWGhCgoKQvuBQIB4AcBlIOrh+iG326309HRVVVW1OsblcjW7uqqvr292Ffa/HA6HHA5H1NYJADBDzL/HdeLECdXW1srtdrc6ZuTIkSopKQk7tnXrVmVlZcV6eQAAw0R8xXX69GkdOXIktF9dXa3KykolJycrOTlZCxcu1JQpU+R2u/XFF1/oqaeeUvfu3XX33XeH5syYMUM9e/aU1+uVJD322GMaNWqUFi9erEmTJqm4uFjbtm3Tzp07o/AUAQDxJOJw7du3T6NHjw7tf/c5U35+vlasWKFPP/1Ua9eu1cmTJ+V2uzV69GitW7dOiYmJoTk1NTXq1On7i72srCy99957euaZZ/Tss8+qT58+WrdunYYPH34pzw0AEIcs27bt9l5ENAQCATmdTmVrkhKszu29HABABJrsRpWqWH6/X0lJSW2O5XcVAgCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCiECwBgFMIFADBKxOEqLy/XxIkT5fF4ZFmWNm7cGPa4ZVktbr///e9bPeeaNWtanHP27NmInxAAIL5FHK4zZ85o4MCBWrZsWYuP19XVhW2rVq2SZVmaMmVKm+dNSkpqNrdr166RLg8AEOcSIp2Qm5ur3NzcVh93uVxh+8XFxRo9erSuv/76Ns9rWVazuQAA/FBMP+P6z3/+o82bN+s3v/nNBceePn1a6enp6tWrl+666y7t37+/zfHBYFCBQCBsAwDEv4ivuCLx5ptvKjExUffcc0+b42666SatWbNGAwYMUCAQ0B//+EfdeuutOnDggDIzM1uc4/V69fzzz8di2UBMfHCs8pLPMc7zs0s+B2C6mF5xrVq1Sr/61a8u+FnViBEj9Otf/1oDBw7U7bffrj//+c+64YYb9Morr7Q6p7CwUH6/P7TV1tZGe/kAgA4oZldcH374oQ4fPqx169ZFPLdTp04aOnSoqqqqWh3jcDjkcDguZYkAAAPF7IrrjTfe0ODBgzVw4MCI59q2rcrKSrnd7hisDABgsoivuE6fPq0jR46E9qurq1VZWank5GT17t1bkhQIBPSXv/xFL730UovnmDFjhnr27Cmv1ytJev755zVixAhlZmYqEAho6dKlqqys1PLlyy/mOQEA4ljE4dq3b59Gjx4d2i8oKJAk5efna82aNZKk9957T7Zt64EHHmjxHDU1NerU6fuLvZMnT+rhhx+Wz+eT0+nUoEGDVF5ermHDhkW6PABAnLNs27bbexHREAgE5HQ6la1JSrA6t/dygGa4qxBoXZPdqFIVy+/3Kykpqc2x/K5CAIBRCBcAwCiECwBgFMIFADAK4QIAGIVwAQCMQrgAAEYhXAAAoxAuAIBRCBcAwCgx/UOSAL7Hr2sCooMrLgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGIVwAAKMQLgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGIVwAAKMQLgCAUQgXAMAohAsAYBTCBQAwSkTh8nq9Gjp0qBITE5WSkqLJkyfr8OHDYWNs29bChQvl8XjUrVs3ZWdn6+DBgxc8d1FRkfr27SuHw6G+fftqw4YNkT0TAMBlIaJwlZWVafbs2dq9e7dKSkrU1NSknJwcnTlzJjTmxRdf1Msvv6xly5Zp7969crlcGjt2rE6dOtXqeSsqKpSXl6fp06frwIEDmj59uqZOnao9e/Zc/DMDAMQly7Zt+2In//e//1VKSorKyso0atQo2bYtj8ejefPmacGCBZKkYDCo1NRULV68WI888kiL58nLy1MgEND7778fOjZ+/Hhdc801evfdd3/UWgKBgJxOp7I1SQlW54t9SgCAdtBkN6pUxfL7/UpKSmpz7CV9xuX3+yVJycnJkqTq6mr5fD7l5OSExjgcDt1xxx3atWtXq+epqKgImyNJ48aNa3NOMBhUIBAI2wAA8e+iw2XbtgoKCnTbbbepf//+kiSfzydJSk1NDRubmpoaeqwlPp8v4jler1dOpzO0paWlXexTAQAY5KLDNWfOHH3yySct/ijPsqywfdu2mx271DmFhYXy+/2hrba2NoLVAwBMlXAxk+bOnatNmzapvLxcvXr1Ch13uVySvr2CcrvdoeP19fXNrqj+l8vlanZ1daE5DodDDofjYpYPADBYRFdctm1rzpw5Wr9+vbZv366MjIywxzMyMuRyuVRSUhI61tDQoLKyMmVlZbV63pEjR4bNkaStW7e2OQcAcHmK6Ipr9uzZeuedd1RcXKzExMTQVZLT6VS3bt1kWZbmzZunRYsWKTMzU5mZmVq0aJGuvPJKTZs2LXSeGTNmqGfPnvJ6vZKkxx57TKNGjdLixYs1adIkFRcXa9u2bdq5c2cUnyoAIB5EFK4VK1ZIkrKzs8OOr169WjNnzpQkPfHEE/rmm2/06KOP6quvvtLw4cO1detWJSYmhsbX1NSoU6fvL/aysrL03nvv6ZlnntGzzz6rPn36aN26dRo+fPhFPi0AQLy6pO9xdSR8jwsAzPWTfY8LAICfGuECABiFcAEAjEK4AABGIVwAAKMQLgCAUQgXAMAohAsAYBTCBQAwCuECABiFcAEAjEK4AABGIVwAAKMQLgCAUQgXAMAohAsAYBTCBQAwSkJ7LyBavvtDzk1qlOLibzoDwOWjSY2Svv9veVviJlynTp2SJO3UlnZeCQDgYp06dUpOp7PNMZb9Y/JmgPPnz+vYsWNKTEyUZVnNHg8EAkpLS1Ntba2SkpLaYYXxhdczung9o4vXM7p+itfTtm2dOnVKHo9HnTq1/SlW3FxxderUSb169brguKSkJN7IUcTrGV28ntHF6xldsX49L3Sl9R1uzgAAGIVwAQCMctmEy+Fw6LnnnpPD4WjvpcQFXs/o4vWMLl7P6Opor2fc3JwBALg8XDZXXACA+EC4AABGIVwAAKMQLgCAUS6bcL366qvKyMhQ165dNXjwYH344YftvSQjLVy4UJZlhW0ul6u9l2WM8vJyTZw4UR6PR5ZlaePGjWGP27athQsXyuPxqFu3bsrOztbBgwfbZ7EGuNDrOXPmzGbv1xEjRrTPYjs4r9eroUOHKjExUSkpKZo8ebIOHz4cNqajvD8vi3CtW7dO8+bN09NPP639+/fr9ttvV25urmpqatp7aUbq16+f6urqQtunn37a3ksyxpkzZzRw4EAtW7asxcdffPFFvfzyy1q2bJn27t0rl8ulsWPHhn4XJ8Jd6PWUpPHjx4e9X7ds4feZtqSsrEyzZ8/W7t27VVJSoqamJuXk5OjMmTOhMR3m/WlfBoYNG2bPmjUr7NhNN91kP/nkk+20InM999xz9sCBA9t7GXFBkr1hw4bQ/vnz522Xy2W/8MILoWNnz561nU6nvXLlynZYoVl++Hratm3n5+fbkyZNapf1mK6+vt6WZJeVldm23bHen3F/xdXQ0KCPPvpIOTk5YcdzcnK0a9eudlqV2aqqquTxeJSRkaH7779fn3/+eXsvKS5UV1fL5/OFvVcdDofuuOMO3quXoLS0VCkpKbrhhhv00EMPqb6+vr2XZAS/3y9JSk5OltSx3p9xH67jx4/r3LlzSk1NDTuempoqn8/XTqsy1/Dhw7V27Vp98MEHev311+Xz+ZSVlaUTJ06099KM9937kfdq9OTm5urtt9/W9u3b9dJLL2nv3r268847FQwG23tpHZpt2yooKNBtt92m/v37S+pY78+4+e3wF/LDP3Vi23aLf/4EbcvNzQ3984ABAzRy5Ej16dNHb775pgoKCtpxZfGD92r05OXlhf65f//+GjJkiNLT07V582bdc8897biyjm3OnDn65JNPtHPnzmaPdYT3Z9xfcXXv3l1XXHFFs/8jqK+vb/Z/DojcVVddpQEDBqiqqqq9l2K87+7O5L0aO263W+np6bxf2zB37lxt2rRJO3bsCPtTUR3p/Rn34erSpYsGDx6skpKSsOMlJSXKyspqp1XFj2AwqEOHDsntdrf3UoyXkZEhl8sV9l5taGhQWVkZ79UoOXHihGpra3m/tsC2bc2ZM0fr16/X9u3blZGREfZ4R3p/XhY/KiwoKND06dM1ZMgQjRw5Uq+99ppqamo0a9as9l6acebPn6+JEyeqd+/eqq+v1//93/8pEAgoPz+/vZdmhNOnT+vIkSOh/erqalVWVio5OVm9e/fWvHnztGjRImVmZiozM1OLFi3SlVdeqWnTprXjqjuutl7P5ORkLVy4UFOmTJHb7dYXX3yhp556St27d9fdd9/djqvumGbPnq133nlHxcXFSkxMDF1ZOZ1OdevWTZZldZz35096D2M7Wr58uZ2enm536dLFvuWWW0K3eCIyeXl5ttvttjt37mx7PB77nnvusQ8ePNjeyzLGjh07bEnNtvz8fNu2v73l+LnnnrNdLpftcDjsUaNG2Z9++mn7LroDa+v1/Prrr+2cnBy7R48edufOne3evXvb+fn5dk1NTXsvu0Nq6XWUZK9evTo0pqO8P/mzJgAAo8T9Z1wAgPhCuAAARiFcAACjEC4AgFEIFwDAKIQLAGAUwgUAMArhAgAYhXABAIxCuAAARiFcAACjEC4AgFH+H9cAmb2OLAPMAAAAAElFTkSuQmCC",
      "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[18,10] = 2000\n",
    "\n",
    "\n",
    "plt.imshow(C_t, vmin=0, vmax=2000)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "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": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulate one time step\n",
    "def simulate(C_t): \n",
    "        C_t1 = np.zeros((grid_size['x']+2,grid_size['y']+2))\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'] # 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'] # 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": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAAGiCAYAAACcWg7FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAWklEQVR4nO3df3RU1b3//9dAYECaTAkxmYyEGF1glXApBCSJFoNoIAqoYMGiAb7SqFfApoGFovWa3k8vUVuFK/H3RSKCQrv42UKRUEyQ8kN+RQEpxTaaoBmjFCaEQhKS8/2DcsqQZMgwEwNzno+19pI5Z++dPdOszjvv/ePYDMMwBAAAIKldWw8AAABcOggMAACAicAAAACYCAwAAICJwAAAAJgIDAAAgInAAAAAmAgMAACAicAAAACYCAwAAICJwAAAgFaUl5engQMHKjw8XNHR0br77rt18OBBrzqGYSg3N1cul0udO3dWWlqa9u/f71WnpqZG06ZNU1RUlLp06aJRo0bp8OHDXnWOHj2qzMxMORwOORwOZWZm6tixY36Nl8AAAIBWVFxcrClTpmjbtm0qLCzU6dOnlZ6erhMnTph1nn/+eb344ovKz8/Xjh075HQ6dfvtt+v48eNmnezsbK1YsUJLlizR5s2bVV1drREjRqi+vt6sM378eJWUlGjdunVat26dSkpKlJmZ6d+ADQAA8J2prKw0JBnFxcWGYRhGQ0OD4XQ6jWeffdasc+rUKcPhcBivvfaaYRiGcezYMaNDhw7GkiVLzDpffvml0a5dO2PdunWGYRjGp59+akgytm3bZtbZunWrIcn4y1/+0uLxhQUWB106Ghoa9NVXXyk8PFw2m62thwMA8INhGDp+/LhcLpfatWu9ZPapU6dUW1sblL4Mw2j0fWO322W3232283g8kqTIyEhJUmlpqdxut9LT0736ueWWW7RlyxY9/PDD2rVrl+rq6rzquFwuJSYmasuWLRo2bJi2bt0qh8OhQYMGmXWSk5PlcDi0ZcsWXXfddS16XyETGHz11VeKi4tr62EAAAJQXl6u7t27t0rfp06dUkL89+SurL9w5Rb43ve+p+rqaq9rzzzzjHJzc5ttYxiGcnJydPPNNysxMVGS5Ha7JUkxMTFedWNiYvTFF1+YdTp27KiuXbs2qnO2vdvtVnR0dKOfGR0dbdZpiZAJDMLDwyVJN+sOhalDG48GAOCP06rTZq01/7+8NdTW1spdWa/SXfGKCA8sK1F1vEEJSV+ovLxcERER5vULZQumTp2qTz75RJs3b2507/zsQ1MZifOdX6ep+i3p51whExicfdNh6qAwG4EBAFxWjDP/+S6mgiPC2wUcGJh9RUR4BQa+TJs2TatXr9amTZu8siJOp1PSmb/4Y2NjzeuVlZVmFsHpdKq2tlZHjx71yhpUVlYqNTXVrPP11183+rnffPNNo2yEL+xKAABYSr3REJTSUoZhaOrUqVq+fLk2btyohIQEr/sJCQlyOp0qLCw0r9XW1qq4uNj80k9KSlKHDh286lRUVGjfvn1mnZSUFHk8Hn300Udmne3bt8vj8Zh1WqLVAoNXXnlFCQkJ6tSpk5KSkvThhx/6rF9cXKykpCR16tRJ11xzjV577bXWGhoAwMIaZASltNSUKVO0aNEivfvuuwoPD5fb7Zbb7dbJkyclncmSZGdna/bs2VqxYoX27dunSZMm6YorrtD48eMlSQ6HQ5MnT9b06dP1pz/9SXv27NEDDzygPn366LbbbpMkXX/99Ro+fLiysrK0bds2bdu2TVlZWRoxYkSLFx5KrTSVsHTpUmVnZ+uVV17RTTfdpNdff10ZGRn69NNP1aNHj0b1S0tLdccddygrK0uLFi3Sn//8Zz366KO68sorNWbMmNYYIgDAohrUoJb/vd98Hy316quvSpLS0tK8ri9YsECTJk2SJM2cOVMnT57Uo48+qqNHj2rQoEFav36915qLOXPmKCwsTGPHjtXJkyc1dOhQFRQUqH379madxYsX67HHHjN3L4waNUr5+fl+vTebYRgtD3taaNCgQerfv7/5YUhnIpm7775beXl5jeo//vjjWr16tQ4cOGBee+SRR/Txxx9r69atLfqZVVVVcjgcStNdrDEAgMvMaaNORVolj8fT4jl7f539nvjqYPegLD50XXe4VcfbVoI+lVBbW6tdu3Z57bWUpPT0dG3ZsqXJNlu3bm1Uf9iwYdq5c6fq6uqabFNTU6OqqiqvAgDAhdQbRlBKqAp6YPDtt9+qvr6+yf2Yze2jdLvdTdY/ffq0vv322ybb5OXlmWdBOxwOzjAAALTId73G4HLTaosP/d2P2VT9pq6fNWvWLHk8HrOUl5cHOGIAABD0xYdRUVFq3759o+zAufsxz+d0OpusHxYWpm7dujXZpiXHTgIAcL4GGaoP8C9+MgZ+6Nixo5KSkrz2WkpSYWFhs/soU1JSGtVfv369BgwYoA4dWEgIAAgephJ8a5WphJycHP3f//2f3nrrLR04cEA///nPVVZWpkceeUTSmWmACRMmmPUfeeQRffHFF8rJydGBAwf01ltvaf78+ZoxY0ZrDA8AADSjVc4xGDdunI4cOaL//u//VkVFhRITE7V27VrFx8dLOnNaU1lZmVk/ISFBa9eu1c9//nO9/PLLcrlceumllzjDAAAQdMHYVRDKuxJa5RyDtsA5BgBw+fouzzH4y4EYhQd4jsHx4w36wfVfc44BAAAIbSHzdEUAAFqiPgi7EgJtfykjMAAAWEq9caYE2keoIjAAAFhKw79KoH2EKtYYAAAAExkDAIClNMimejV/RH9L+whVBAYAAEtpMM6UQPsIVUwlAAAAExkDAICl1AdhKiHQ9pcyAgMAgKUQGPjGVAIAADCRMQAAWEqDYVODEeCuhADbX8oIDAAAlsJUgm9MJQAAABMZAwCApdSrneoD/Lu4PkhjuRQRGAAALMUIwhoDgzUGAACEBtYY+MYaAwAAYCJjAACwlHqjneqNANcYhPCzEggMAACW0iCbGgJMmDcodCMDphIAAICJjAEAwFJYfOgbgQEAwFKCs8aAqQQAAGABZAwAAJZyZvFhgA9RYioBAIDQ0BCEI5HZlQAAACyBwAAAYClnFx8GWvyxadMmjRw5Ui6XSzabTStXrvS6b7PZmiy//vWvzTppaWmN7t93331e/Rw9elSZmZlyOBxyOBzKzMzUsWPH/BorgQEAwFIa1C4oxR8nTpxQ3759lZ+f3+T9iooKr/LWW2/JZrNpzJgxXvWysrK86r3++ute98ePH6+SkhKtW7dO69atU0lJiTIzM/0aK2sMAACWUm/YVB/g0xH9bZ+RkaGMjIxm7zudTq/Xq1at0pAhQ3TNNdd4Xb/iiisa1T3rwIEDWrdunbZt26ZBgwZJkt58802lpKTo4MGDuu6661o0VjIGAABcpKqqKq9SU1MTcJ9ff/211qxZo8mTJze6t3jxYkVFRal3796aMWOGjh8/bt7bunWrHA6HGRRIUnJyshwOh7Zs2dLinx/0wCAvL08DBw5UeHi4oqOjdffdd+vgwYM+2xQVFTU5t/KXv/wl2MMDAFhc/b92JQRaJCkuLs6cz3c4HMrLywt4fG+//bbCw8M1evRor+v333+/3nvvPRUVFenpp5/WsmXLvOq43W5FR0c36i86Olput7vFPz/oUwnFxcWaMmWKBg4cqNOnT+upp55Senq6Pv30U3Xp0sVn24MHDyoiIsJ8feWVVwZ7eAAAi2sw2qkhwJMPG/518mF5ebnX95bdbg+oX0l66623dP/996tTp05e17Oyssx/JyYmqmfPnhowYIB2796t/v37SzqziPF8hmE0eb05QQ8M1q1b5/V6wYIFio6O1q5duzR48GCfbaOjo/X9738/2EMCAKBVREREeAUGgfrwww918OBBLV269IJ1+/fvrw4dOujQoUPq37+/nE6nvv7660b1vvnmG8XExLR4DK2+xsDj8UiSIiMjL1i3X79+io2N1dChQ/XBBx/4rFtTU9NobgcAgAsJ5lRCsM2fP19JSUnq27fvBevu379fdXV1io2NlSSlpKTI4/Hoo48+Muts375dHo9HqampLR5Dq+5KMAxDOTk5uvnmm5WYmNhsvdjYWL3xxhtKSkpSTU2N3nnnHQ0dOlRFRUXNZhny8vL0y1/+srWGDgAIUQ3yf1dBU334o7q6Wp999pn5urS0VCUlJYqMjFSPHj0knVnI+Lvf/U4vvPBCo/Z/+9vftHjxYt1xxx2KiorSp59+qunTp6tfv3666aabJEnXX3+9hg8frqysLHMb40MPPaQRI0a0eEeCJNkMo/UeETVlyhStWbNGmzdvVvfu3f1qO3LkSNlsNq1evbrJ+zU1NV6rP6uqqhQXF6c03aUwW4eAxg0A+G6dNupUpFXyeDxBTc2fq6qqSg6HQ6/vTlLn7wX2d/HJ6tN6uP+uFo+3qKhIQ4YMaXR94sSJKigokCS98cYbys7OVkVFhRwOh1e98vJyPfDAA9q3b5+qq6sVFxenO++8U88884xXRv4f//iHHnvsMfO7c9SoUcrPz/drmr7VMgbTpk3T6tWrtWnTJr+DAunMFotFixY1e99utwdlkQcAwFou5oCipvrwR1pami70d/hDDz2khx56qMl7cXFxKi4uvuDPiYyM9Pnd2RJBDwwMw9C0adO0YsUKFRUVKSEh4aL62bNnjzlvAgBAsFzMkcZN9RGqgh4YTJkyRe+++65WrVql8PBwc++kw+FQ586dJUmzZs3Sl19+qYULF0qS5s6dq6uvvlq9e/dWbW2tFi1apGXLlmnZsmXBHh4AAPAh6IHBq6++KulM2uRcCxYs0KRJkySdORO6rKzMvFdbW6sZM2boyy+/VOfOndW7d2+tWbNGd9xxR7CHBwCwuAbZ1KBAFx8G1v5S1ipTCRdydqHFWTNnztTMmTODPRQAABphKsE3HqIEALCUYJxD0FrnGFwKQvedAQAAv5ExAABYSoNhU0OgBxwF2P5SRmAAALCUhiBMJQR6DsKlLHTfGQAA8BsZAwCApQTnscuh+3c1gQEAwFLqZVN9gOcQBNr+Uha6IQ8AAPAbGQMAgKUwleAbgQEAwFLqFfhUQH1whnJJCt2QBwAA+I2MAQDAUphK8I3AAABgKTxEyTcCAwCApRhBeOyywXZFAABgBWQMAACWwlSCbwQGAABL4emKvoVuyAMAAPxGxgAAYCn1QXjscqDtL2UEBgAAS2EqwbfQDXkAAIDfyBgAACylQe3UEODfxYG2v5QRGAAALKXesKk+wKmAQNtfykI35AEAAH4jYwAAsBQWH/pGYAAAsBQjCE9XNDj5EACA0FAvm+oDfAhSoO0vZaEb8gAAAL+RMQAAWEqDEfgagQYjSIO5BBEYAAAspSEIawwCbX8pC913BgAA/Bb0wCA3N1c2m82rOJ1On22Ki4uVlJSkTp066ZprrtFrr70W7GEBACBJapAtKMUfmzZt0siRI+VyuWSz2bRy5Uqv+5MmTWr03ZmcnOxVp6amRtOmTVNUVJS6dOmiUaNG6fDhw151jh49qszMTDkcDjkcDmVmZurYsWN+jbVVMga9e/dWRUWFWfbu3dts3dLSUt1xxx360Y9+pD179ujJJ5/UY489pmXLlrXG0AAAFnf25MNAiz9OnDihvn37Kj8/v9k6w4cP9/ruXLt2rdf97OxsrVixQkuWLNHmzZtVXV2tESNGqL6+3qwzfvx4lZSUaN26dVq3bp1KSkqUmZnp11hbZY1BWFjYBbMEZ7322mvq0aOH5s6dK0m6/vrrtXPnTv3mN7/RmDFjWmN4AAB8pzIyMpSRkeGzjt1ub/a70+PxaP78+XrnnXd02223SZIWLVqkuLg4bdiwQcOGDdOBAwe0bt06bdu2TYMGDZIkvfnmm0pJSdHBgwd13XXXtWisrZIxOHTokFwulxISEnTffffp73//e7N1t27dqvT0dK9rw4YN086dO1VXV9dsu5qaGlVVVXkVAAAu5Oziw0CLpEbfQzU1NRc9rqKiIkVHR6tXr17KyspSZWWleW/Xrl2qq6vz+r50uVxKTEzUli1bJJ35PnU4HGZQIEnJyclyOBxmnZYIemAwaNAgLVy4UO+//77efPNNud1upaam6siRI03Wd7vdiomJ8boWExOj06dP69tvv2325+Tl5ZlzKA6HQ3FxcUF9HwCA0NQgm3ks8kWXf60xiIuL8/ouysvLu6gxZWRkaPHixdq4caNeeOEF7dixQ7feeqsZaLjdbnXs2FFdu3b1ahcTEyO3223WiY6ObtR3dHS0Waclgj6VcG6qpE+fPkpJSdG1116rt99+Wzk5OU22sdm852oMw2jy+rlmzZrl1V9VVRXBAQDgO1VeXq6IiAjztd1uv6h+xo0bZ/47MTFRAwYMUHx8vNasWaPRo0c3284wDK/vyqa+N8+vcyGtfo5Bly5d1KdPHx06dKjJ+06ns1EkU1lZqbCwMHXr1q3Zfu12+0X/DwAAsC7jInYVNNWHJEVERHgFBsESGxur+Ph487vT6XSqtrZWR48e9coaVFZWKjU11azz9ddfN+rrm2++aZSZ96XVzzGoqanRgQMHFBsb2+T9lJQUFRYWel1bv369BgwYoA4dOrT28AAAFhPwNEIQns54IUeOHFF5ebn53ZmUlKQOHTp4fV9WVFRo3759ZmCQkpIij8ejjz76yKyzfft2eTwes05LBD0wmDFjhoqLi1VaWqrt27fr3nvvVVVVlSZOnCjpzBTAhAkTzPqPPPKIvvjiC+Xk5OjAgQN66623NH/+fM2YMSPYQwMAIKiLD1uqurpaJSUlKikpkXRmq35JSYnKyspUXV2tGTNmaOvWrfr8889VVFSkkSNHKioqSvfcc48kyeFwaPLkyZo+fbr+9Kc/ac+ePXrggQfUp08fc5fC9ddfr+HDhysrK0vbtm3Ttm3blJWVpREjRrR4R4LUClMJhw8f1k9+8hN9++23uvLKK5WcnKxt27YpPj5e0pkIp6yszKyfkJCgtWvX6uc//7lefvlluVwuvfTSS2xVBACEjJ07d2rIkCHm67Nr5CZOnKhXX31Ve/fu1cKFC3Xs2DHFxsZqyJAhWrp0qcLDw802c+bMUVhYmMaOHauTJ09q6NChKigoUPv27c06ixcv1mOPPWbuXhg1apTPsxOaYjPOrvS7zFVVVcnhcChNdynMxhQEAFxOTht1KtIqeTyeVpmzl/79PXHX+gfVoUvHgPqqO1GrVelvtep42woPUQIAWMrFHGncVB+hiocoAQAAExkDAIClBGNXQWvvSmhLBAYAAEshMPCNqQQAAGAiYwAAsBQyBr4RGAAALIXAwDemEgAAgImMAQDAUgwFfg5BSJwM2AwCAwCApTCV4BuBAQDAUggMfGONAQAAMJExAABYChkD3wgMAACWQmDgG1MJAADARMYAAGAphmGTEeBf/IG2v5QRGAAALKVBtoDPMQi0/aWMqQQAAGAiYwAAsBQWH/pGYAAAsBTWGPjGVAIAADCRMQAAWApTCb4RGAAALIWpBN8IDAAAlmIEIWMQyoEBawwAAICJjAEAwFIMSYYReB+hisAAAGApDbLJxsmHzWIqAQAAmMgYAAAshV0JvhEYAAAspcGwycY5Bs1iKgEAAJiCHhhcffXVstlsjcqUKVOarF9UVNRk/b/85S/BHhoAADKM4JRQFfSphB07dqi+vt58vW/fPt1+++368Y9/7LPdwYMHFRERYb6+8sorgz00AABYY3ABQQ8Mzv9Cf/bZZ3Xttdfqlltu8dkuOjpa3//+94M9HAAA4IdWXWNQW1urRYsW6cEHH5TN5ju66tevn2JjYzV06FB98MEHF+y7pqZGVVVVXgUAgAs5mzEItISqVg0MVq5cqWPHjmnSpEnN1omNjdUbb7yhZcuWafny5bruuus0dOhQbdq0yWffeXl5cjgcZomLiwvy6AEAoejs0xUDLaGqVQOD+fPnKyMjQy6Xq9k61113nbKystS/f3+lpKTolVde0Z133qnf/OY3PvueNWuWPB6PWcrLy4M9fABACGqLxYebNm3SyJEj5XK5ZLPZtHLlSvNeXV2dHn/8cfXp00ddunSRy+XShAkT9NVXX3n1kZaW1mih/n333edV5+jRo8rMzDT/aM7MzNSxY8f8GmurBQZffPGFNmzYoJ/+9Kd+t01OTtahQ4d81rHb7YqIiPAqAABcik6cOKG+ffsqPz+/0b1//vOf2r17t55++mnt3r1by5cv11//+leNGjWqUd2srCxVVFSY5fXXX/e6P378eJWUlGjdunVat26dSkpKlJmZ6ddYW+2AowULFig6Olp33nmn32337Nmj2NjYVhgVAMDqzvzFH+iuhDP/PX99m91ul91ub1Q/IyNDGRkZTfblcDhUWFjodW3evHm68cYbVVZWph49epjXr7jiCjmdzib7OXDggNatW6dt27Zp0KBBkqQ333xTKSkpOnjwoK677roWvbdWyRg0NDRowYIFmjhxosLCvGOPWbNmacKECebruXPnauXKlTp06JD279+vWbNmadmyZZo6dWprDA0AYHHBXHwYFxfntd4tLy8vKGP0eDyy2WyNdustXrxYUVFR6t27t2bMmKHjx4+b97Zu3SqHw2EGBdKZDLzD4dCWLVta/LNbJWOwYcMGlZWV6cEHH2x0r6KiQmVlZebr2tpazZgxQ19++aU6d+6s3r17a82aNbrjjjtaY2gAAARNeXm511R2U9kCf506dUpPPPGExo8f79X3/fffr4SEBDmdTu3bt0+zZs3Sxx9/bGYb3G63oqOjG/UXHR0tt9vd4p/fKoFBenq6jGZWZhQUFHi9njlzpmbOnNkawwAAoBHjXyXQPiQFfY1bXV2d7rvvPjU0NOiVV17xupeVlWX+OzExUT179tSAAQO0e/du9e/fX5KaPBrAMIwLHhlwLp6VAACwlEv1HIO6ujqNHTtWpaWlKiwsvGDA0b9/f3Xo0MFcrO90OvX11183qvfNN98oJiamxeMgMAAAoI2dDQoOHTqkDRs2qFu3bhdss3//ftXV1ZmL9VNSUuTxePTRRx+ZdbZv3y6Px6PU1NQWj4XHLgMArCWYcwktVF1drc8++8x8XVpaqpKSEkVGRsrlcunee+/V7t279Yc//EH19fXmmoDIyEh17NhRf/vb37R48WLdcccdioqK0qeffqrp06erX79+uummmyRJ119/vYYPH66srCxzG+NDDz2kESNGtHhHgkRgAACwmmBMBfjZfufOnRoyZIj5OicnR5I0ceJE5ebmavXq1ZKkH/7wh17tPvjgA6Wlpaljx47605/+pP/93/9VdXW14uLidOedd+qZZ55R+/btzfqLFy/WY489pvT0dEnSqFGjmjw7wRcCAwCApQTjscn+tk9LS2t2Uf6Z/nx3GBcXp+Li4gv+nMjISC1atMi/wZ2HNQYAAMBExgAAYCnB2FUQyk9XJDAAAFiLYfN7jUCTfYQophIAAICJjAEAwFLaYvHh5YTAAABgLW1wjsHlhKkEAABgImMAALAUdiX4RmAAALCeEJ4KCBRTCQAAwETGAABgKUwl+EZgAACwFnYl+ERgAACwGNu/SqB9hCbWGAAAABMZAwCAtTCV4BOBAQDAWggMfGIqAQAAmMgYAACshccu+0RgAACwFJ6u6BtTCQAAwETGAABgLSw+9InAAABgLawx8ImpBAAAYCJjAACwFJtxpgTaR6giMAAAWAtrDHwiMAAAWAtrDHzye43Bpk2bNHLkSLlcLtlsNq1cudLrvmEYys3NlcvlUufOnZWWlqb9+/dfsN9ly5bphhtukN1u1w033KAVK1b4OzQAABAgvwODEydOqG/fvsrPz2/y/vPPP68XX3xR+fn52rFjh5xOp26//XYdP3682T63bt2qcePGKTMzUx9//LEyMzM1duxYbd++3d/hAQDgmxGkEqJshnHx5zfZbDatWLFCd999t6Qz2QKXy6Xs7Gw9/vjjkqSamhrFxMToueee08MPP9xkP+PGjVNVVZX++Mc/mteGDx+url276r333mvRWKqqquRwOJSmuxRm63CxbwkA0AZOG3Uq0ip5PB5FRES0ys84+z0R98L/U7vOnQLqq+HkKZVPf7pVx9tWgrpdsbS0VG63W+np6eY1u92uW265RVu2bGm23datW73aSNKwYcN8tqmpqVFVVZVXAQAAgQlqYOB2uyVJMTExXtdjYmLMe82187dNXl6eHA6HWeLi4gIYOQDAMphK8KlVDjiy2bxXaxqG0ehaoG1mzZolj8djlvLy8osfMADAOs7uSgi0hKigbld0Op2SzmQAYmNjzeuVlZWNMgLntzs/O3ChNna7XXa7PcARAwCAcwU1Y5CQkCCn06nCwkLzWm1trYqLi5Wamtpsu5SUFK82krR+/XqfbQAAuBhnTz4MtIQqvzMG1dXV+uyzz8zXpaWlKikpUWRkpHr06KHs7GzNnj1bPXv2VM+ePTV79mxdccUVGj9+vNlmwoQJuuqqq5SXlydJ+tnPfqbBgwfrueee01133aVVq1Zpw4YN2rx5cxDeIgAA5+DkQ5/8zhjs3LlT/fr1U79+/SRJOTk56tevn/7rv/5LkjRz5kxlZ2fr0Ucf1YABA/Tll19q/fr1Cg8PN/soKytTRUWF+To1NVVLlizRggUL9B//8R8qKCjQ0qVLNWjQoEDfHwAAbS4YhwPW1NRo2rRpioqKUpcuXTRq1CgdPnzYq87Ro0eVmZlpLszPzMzUsWPH/BprQOcYXEo4xwAALl/f5TkGPZ77VVDOMSh7/BctHu8f//hH/fnPf1b//v01ZswYrzOAJOm5557T//zP/6igoEC9evXSr371K23atEkHDx40/7D+z//8T/3+979XQUGBunXrpunTp+sf//iHdu3apfbt20uSMjIydPjwYb3xxhuSpIceekhXX321fv/737f4vfGsBACApdgUhKcr/uu/55+h09zC+IyMDGVkZDTZl2EYmjt3rp566imNHj1akvT2228rJiZG7777rh5++GF5PB7Nnz9f77zzjm677TZJ0qJFixQXF6cNGzZo2LBhOnDggNatW6dt27aZGfc333xTKSkpOnjwoK677roWvbdW2a4IAMAlK4jbFePi4rzO1Dm7ds4fLTkccNeuXaqrq/Oq43K5lJiYaNbZunWrHA6H1zR8cnKyHA6HzwMDz0fGAACAi1ReXu41lXAx2+h9HQ74xRdfmHU6duyorl27Nqpztr3b7VZ0dHSj/qOjo30eGHg+AgMAgLUEcVdCRERE0NZEXMzhgOfXaap+S/o5F1MJAABrucSORD73cMBznXvQn9PpVG1trY4ePeqzztdff92o/2+++cbngYHnIzAAAKANteRwwKSkJHXo0MGrTkVFhfbt22fWSUlJkcfj0UcffWTW2b59uzwej18HBjKVAACwlGCcXOhv+0APB3Q4HJo8ebKmT5+ubt26KTIyUjNmzFCfPn3MXQrXX3+9hg8frqysLL3++uuSzmxXHDFiRIt3JEgEBgAAq2mDkw937typIUOGmK9zcnIkSRMnTlRBQYFmzpypkydP6tFHH9XRo0c1aNCgRocDzpkzR2FhYRo7dqxOnjypoUOHqqCgwDzDQJIWL16sxx57zNy9MGrUKOXn5/s1Vg44AgC0ue/ygKOrf/U/atcpwAOOTp3S5794qlXH21bIGAAArIVnJfhEYAAAsJS2WGNwOWFXAgAAMJExAABYyzlHGgfUR4giMAAAWAtrDHwiMAAAWAprDHxjjQEAADCRMQAAWAtTCT4RGAAArCUIUwmhHBgwlQAAAExkDAAA1sJUgk8EBgAAayEw8ImpBAAAYCJjAACwFM4x8I2MAQAAMBEYAAAAE1MJAABrYfGhTwQGAABLYY2BbwQGAADrCeEv9kCxxgAAAJjIGAAArIU1Bj4RGAAALIU1Br75PZWwadMmjRw5Ui6XSzabTStXrjTv1dXV6fHHH1efPn3UpUsXuVwuTZgwQV999ZXPPgsKCmSz2RqVU6dO+f2GAADAxfM7MDhx4oT69u2r/Pz8Rvf++c9/avfu3Xr66ae1e/duLV++XH/96181atSoC/YbERGhiooKr9KpUyd/hwcAgG9GkEqI8nsqISMjQxkZGU3eczgcKiws9Lo2b9483XjjjSorK1OPHj2a7ddms8npdPo7HAAA/MJUgm+tvivB4/HIZrPp+9//vs961dXVio+PV/fu3TVixAjt2bPHZ/2amhpVVVV5FQAAEJhWDQxOnTqlJ554QuPHj1dERESz9X7wgx+ooKBAq1ev1nvvvadOnTrppptu0qFDh5ptk5eXJ4fDYZa4uLjWeAsAgFDDVIJPrRYY1NXV6b777lNDQ4NeeeUVn3WTk5P1wAMPqG/fvvrRj36k3/72t+rVq5fmzZvXbJtZs2bJ4/GYpby8PNhvAQAQiggMfGqV7Yp1dXUaO3asSktLtXHjRp/Zgqa0a9dOAwcO9JkxsNvtstvtgQ4VAACcI+gZg7NBwaFDh7RhwwZ169bN7z4Mw1BJSYliY2ODPTwAgMWdXXwYaAlVfmcMqqur9dlnn5mvS0tLVVJSosjISLlcLt17773avXu3/vCHP6i+vl5ut1uSFBkZqY4dO0qSJkyYoKuuukp5eXmSpF/+8pdKTk5Wz549VVVVpZdeekklJSV6+eWXg/EeAQD4N04+9MnvwGDnzp0aMmSI+TonJ0eSNHHiROXm5mr16tWSpB/+8Ide7T744AOlpaVJksrKytSu3b+TFceOHdNDDz0kt9sth8Ohfv36adOmTbrxxhv9HR4AAL4RGPjkd2CQlpYmw2j+E/F176yioiKv13PmzNGcOXP8HQoAAAgynpUAALAUDjjyjcAAAGAtTCX41OonHwIAYGVXX311kw8KnDJliiRp0qRJje4lJyd79VFTU6Np06YpKipKXbp00ahRo3T48OFWGS+BAQDAUr7r7Yo7duzwekDg2WcK/fjHPzbrDB8+3KvO2rVrvfrIzs7WihUrtGTJEm3evFnV1dUaMWKE6uvrg/KZnIupBACAtXzHUwlXXnml1+tnn31W1157rW655Rbzmt1ub/ZBgh6PR/Pnz9c777yj2267TZK0aNEixcXFacOGDRo2bJj/4/eBjAEAABfp/If51dTU+KxfW1urRYsW6cEHH5TNZjOvFxUVKTo6Wr169VJWVpYqKyvNe7t27VJdXZ3S09PNay6XS4mJidqyZUvQ3xOBAQDAWoL4rIS4uDivB/qdPbivOStXrtSxY8c0adIk81pGRoYWL16sjRs36oUXXtCOHTt06623mkGG2+1Wx44d1bVrV6++YmJizEMEg4mpBACApdj+VQLtQ5LKy8u9ngd0oWf4zJ8/XxkZGXK5XOa1cePGmf9OTEzUgAEDFB8frzVr1mj06NHN9mUYhlfWIVgIDAAAuEgREREtflDgF198oQ0bNmj58uU+68XGxio+Pt58kKDT6VRtba2OHj3qlTWorKxUamrqxQ++GUwlAACspY0eu7xgwQJFR0frzjvv9FnvyJEjKi8vNx8kmJSUpA4dOpi7GSSpoqJC+/bta5XAgIwBAMBS2uLkw4aGBi1YsEATJ05UWNi/v3qrq6uVm5urMWPGKDY2Vp9//rmefPJJRUVF6Z577pEkORwOTZ48WdOnT1e3bt0UGRmpGTNmqE+fPuYuhWAiMAAAWEsbnHy4YcMGlZWV6cEHH/S63r59e+3du1cLFy7UsWPHFBsbqyFDhmjp0qUKDw83682ZM0dhYWEaO3asTp48qaFDh6qgoEDt27cP8I00RmAAAEArS09Pb/Ihg507d9b7779/wfadOnXSvHnzNG/evNYYnhcCAwCA9YTwsw4CRWAAALAUnq7oG7sSAACAiYwBAMBaeOyyTwQGAABLYSrBN6YSAACAiYwBAMBamErwicAAAGApTCX4xlQCAAAwkTEAAFgLUwk+ERgAAKyFwMAnAgMAgKWwxsA31hgAAAATGQMAgLUwleATgQEAwFJshiFbE49A9rePUMVUAgAAMJExAABYC1MJPhEYAAAshV0Jvvk9lbBp0yaNHDlSLpdLNptNK1eu9Lo/adIk2Ww2r5KcnHzBfpctW6YbbrhBdrtdN9xwg1asWOHv0AAAQID8DgxOnDihvn37Kj8/v9k6w4cPV0VFhVnWrl3rs8+tW7dq3LhxyszM1Mcff6zMzEyNHTtW27dv93d4AAD4ZgSphCi/pxIyMjKUkZHhs47dbpfT6Wxxn3PnztXtt9+uWbNmSZJmzZql4uJizZ07V++9956/QwQAoFlMJfjWKrsSioqKFB0drV69eikrK0uVlZU+62/dulXp6ele14YNG6YtW7Y026ampkZVVVVeBQAABCbogUFGRoYWL16sjRs36oUXXtCOHTt06623qqamptk2brdbMTExXtdiYmLkdrubbZOXlyeHw2GWuLi4oL0HAEAIYyrBp6DvShg3bpz578TERA0YMEDx8fFas2aNRo8e3Ww7m83m9dowjEbXzjVr1izl5OSYr6uqqggOAAAXxFSCb62+XTE2Nlbx8fE6dOhQs3WcTmej7EBlZWWjLMK57Ha77HZ70MYJALAIzjHwqdVPPjxy5IjKy8sVGxvbbJ2UlBQVFhZ6XVu/fr1SU1Nbe3gAAOAcfmcMqqur9dlnn5mvS0tLVVJSosjISEVGRio3N1djxoxRbGysPv/8cz355JOKiorSPffcY7aZMGGCrrrqKuXl5UmSfvazn2nw4MF67rnndNddd2nVqlXasGGDNm/eHIS3CACAt1CeCgiU34HBzp07NWTIEPP12Xn+iRMn6tVXX9XevXu1cOFCHTt2TLGxsRoyZIiWLl2q8PBws01ZWZnatft3siI1NVVLlizRL37xCz399NO69tprtXTpUg0aNCiQ9wYAQGOGcaYE2keI8jswSEtLk+HjA3n//fcv2EdRUVGja/fee6/uvfdef4cDAACCiGclAAAshV0JvhEYAACshV0JPrX6rgQAAHD5IGMAALAUW8OZEmgfoYrAAABgLUwl+MRUAgAAMBEYAAAs5eyuhEBLS+Xm5spms3kVp9Np3jcMQ7m5uXK5XOrcubPS0tK0f/9+rz5qamo0bdo0RUVFqUuXLho1apQOHz4crI/EC4EBAMBazh5wFGjxQ+/evVVRUWGWvXv3mveef/55vfjii8rPz9eOHTvkdDp1++236/jx42ad7OxsrVixQkuWLNHmzZtVXV2tESNGqL6+Pmgfy1msMQAAWEpbnGMQFhbmlSU4yzAMzZ07V0899ZT5BOK3335bMTExevfdd/Xwww/L4/Fo/vz5euedd3TbbbdJkhYtWqS4uDht2LBBw4YNC+zNnIeMAQAAF6mqqsqr1NTUNFnv0KFDcrlcSkhI0H333ae///3vks48b8jtdis9Pd2sa7fbdcstt2jLli2SpF27dqmurs6rjsvlUmJiolknmAgMAADWYgSpSIqLi5PD4TDL2YcDnmvQoEFauHCh3n//fb355ptyu91KTU3VkSNH5Ha7JUkxMTFebWJiYsx7brdbHTt2VNeuXZutE0xMJQAALCWYUwnl5eWKiIgwr9vt9kZ1MzIyzH/36dNHKSkpuvbaa/X2228rOTn5TH82m1cbwzAaXTtfS+pcDDIGAABcpIiICK/SVGBwvi5duqhPnz46dOiQue7g/L/8KysrzSyC0+lUbW2tjh492mydYCIwAABYSxvsSjhXTU2NDhw4oNjYWCUkJMjpdKqwsNC8X1tbq+LiYqWmpkqSkpKS1KFDB686FRUV2rdvn1knmJhKAABYyne9K2HGjBkaOXKkevToocrKSv3qV79SVVWVJk6cKJvNpuzsbM2ePVs9e/ZUz549NXv2bF1xxRUaP368JMnhcGjy5MmaPn26unXrpsjISM2YMUN9+vQxdykEE4EBAACt6PDhw/rJT36ib7/9VldeeaWSk5O1bds2xcfHS5JmzpypkydP6tFHH9XRo0c1aNAgrV+/XuHh4WYfc+bMUVhYmMaOHauTJ09q6NChKigoUPv27YM+XpthBJAPuYRUVVXJ4XAoTXcpzNahrYcDAPDDaaNORVolj8fjtZgvmM5+T6QM/2+FdegUUF+n605p67r/atXxthUyBgAAS2mLA44uJyw+BAAAJjIGAABraTDOlED7CFEEBgAAaznn5MKA+ghRBAYAAEuxKQhrDIIykksTawwAAICJjAEAwFoCPLnQ7CNEERgAACyF7Yq+MZUAAABMZAwAANbCrgSfCAwAAJZiMwzZAlwjEGj7SxlTCQAAwETGAABgLQ3/KoH2EaIIDAAAlsJUgm9+TyVs2rRJI0eOlMvlks1m08qVK73u22y2Jsuvf/3rZvssKChoss2pU6f8fkMAAODi+R0YnDhxQn379lV+fn6T9ysqKrzKW2+9JZvNpjFjxvjsNyIiolHbTp0Ce142AACNGEEqIcrvqYSMjAxlZGQ0e9/pdHq9XrVqlYYMGaJrrrnGZ782m61RWwAAgo6TD31q1V0JX3/9tdasWaPJkydfsG51dbXi4+PVvXt3jRgxQnv27PFZv6amRlVVVV4FAIALOXvyYaAlVLVqYPD2228rPDxco0eP9lnvBz/4gQoKCrR69Wq999576tSpk2666SYdOnSo2TZ5eXlyOBxmiYuLC/bwAQCwnFYNDN566y3df//9F1wrkJycrAceeEB9+/bVj370I/32t79Vr169NG/evGbbzJo1Sx6Pxyzl5eXBHj4AIBSdnUoItISoVtuu+OGHH+rgwYNaunSp323btWungQMH+swY2O122e32QIYIALAgW8OZEmgfoarVMgbz589XUlKS+vbt63dbwzBUUlKi2NjYVhgZAABojt8Zg+rqan322Wfm69LSUpWUlCgyMlI9evSQJFVVVel3v/udXnjhhSb7mDBhgq666irl5eVJkn75y18qOTlZPXv2VFVVlV566SWVlJTo5Zdfvpj3BABA89iV4JPfgcHOnTs1ZMgQ83VOTo4kaeLEiSooKJAkLVmyRIZh6Cc/+UmTfZSVlaldu38nK44dO6aHHnpIbrdbDodD/fr106ZNm3TjjTf6OzwAAHzj6Yo+2QwjNMKeqqoqORwOpekuhdk6tPVwAAB+OG3UqUir5PF4FBER0So/w/yeGPiUwsICO0Dv9OlTKtrxP6063rbCsxIAAJbCsxJ8IzAAAFgLawx8atVzDAAAwOWFjAEAwFoMSYGeQxC6CQMCAwCAtbDGwDcCAwCAtRgKwhqDoIzkksQaAwAAYCJjAACwFnYl+ERgAACwlgZJtiD0EaKYSgAAACYCAwCApZzdlRBoaam8vDwNHDhQ4eHhio6O1t13362DBw961Zk0aZJsNptXSU5O9qpTU1OjadOmKSoqSl26dNGoUaN0+PDhoHwm5yIwAABYy9k1BoGWFiouLtaUKVO0bds2FRYW6vTp00pPT9eJEye86g0fPlwVFRVmWbt2rdf97OxsrVixQkuWLNHmzZtVXV2tESNGqL6+Pigfy1msMQAAoBWtW7fO6/WCBQsUHR2tXbt2afDgweZ1u90up9PZZB8ej0fz58/XO++8o9tuu02StGjRIsXFxWnDhg0aNmxY0MZLxgAAYC1BzBhUVVV5lZqamgv+eI/HI0mKjIz0ul5UVKTo6Gj16tVLWVlZqqysNO/t2rVLdXV1Sk9PN6+5XC4lJiZqy5YtwfhUTAQGAABrCWJgEBcXJ4fDYZa8vLwL/GhDOTk5uvnmm5WYmGhez8jI0OLFi7Vx40a98MIL2rFjh2699VYz0HC73erYsaO6du3q1V9MTIzcbndQPx6mEgAAuEjl5eWKiIgwX9vtdp/1p06dqk8++USbN2/2uj5u3Djz34mJiRowYIDi4+O1Zs0ajR49utn+DMOQzRbo3ktvZAwAANbSEKQiKSIiwqv4CgymTZum1atX64MPPlD37t19DjE2Nlbx8fE6dOiQJMnpdKq2tlZHjx71qldZWamYmBi/3v6FEBgAACzlu96uaBiGpk6dquXLl2vjxo1KSEi4YJsjR46ovLxcsbGxkqSkpCR16NBBhYWFZp2Kigrt27dPqamp/n8IPjCVAACwlu/4SOQpU6bo3Xff1apVqxQeHm6uCXA4HOrcubOqq6uVm5urMWPGKDY2Vp9//rmefPJJRUVF6Z577jHrTp48WdOnT1e3bt0UGRmpGTNmqE+fPuYuhWAhMAAAoBW9+uqrkqS0tDSv6wsWLNCkSZPUvn177d27VwsXLtSxY8cUGxurIUOGaOnSpQoPDzfrz5kzR2FhYRo7dqxOnjypoUOHqqCgQO3btw/qeAkMAADW0mBItgAzBg3+TSX40rlzZ73//vsX7KdTp06aN2+e5s2b1+KffTEIDAAA1sLTFX1i8SEAADCRMQAAWEwQMgYK3YwBgQEAwFqYSvCJqQQAAGAiYwAAsJYGQwFPBfixK+FyQ2AAALAWo+FMCbSPEMVUAgAAMJExAABYC4sPfSIwAABYC2sMfPJrKiEvL08DBw5UeHi4oqOjdffdd+vgwYNedQzDUG5urlwulzp37qy0tDTt37//gn0vW7ZMN9xwg+x2u2644QatWLHCv3cCAEBLnM0YBFpClF+BQXFxsaZMmaJt27apsLBQp0+fVnp6uk6cOGHWef755/Xiiy8qPz9fO3bskNPp1O23367jx4832+/WrVs1btw4ZWZm6uOPP1ZmZqbGjh2r7du3X/w7AwAAfrMZF3q6gw/ffPONoqOjVVxcrMGDB8swDLlcLmVnZ+vxxx+XJNXU1CgmJkbPPfecHn744Sb7GTdunKqqqvTHP/7RvDZ8+HB17dpV7733XovGUlVVJYfDoTTdpTBbh4t9SwCANnDaqFORVsnj8SgiIqJVfsbZ74nbYh9WWLuOAfV1uqFWGypeb9XxtpWAdiV4PB5JUmRkpCSptLRUbrdb6enpZh273a5bbrlFW7ZsabafrVu3erWRpGHDhvlsU1NTo6qqKq8CAMAFMZXg00UHBoZhKCcnRzfffLMSExMlSW63W5IUExPjVTcmJsa81xS32+13m7y8PDkcDrPExcVd7FsBAAD/ctGBwdSpU/XJJ580meq32Wxerw3DaHQt0DazZs2Sx+MxS3l5uR+jBwBYVkNDcEqIuqjtitOmTdPq1au1adMmde/e3bzudDolnckAxMbGmtcrKysbZQTO5XQ6G2UHLtTGbrfLbrdfzPABAFbGOQY++ZUxMAxDU6dO1fLly7Vx40YlJCR43U9ISJDT6VRhYaF5rba2VsXFxUpNTW2235SUFK82krR+/XqfbQAAQPD5lTGYMmWK3n33Xa1atUrh4eHmX/kOh0OdO3eWzWZTdna2Zs+erZ49e6pnz56aPXu2rrjiCo0fP97sZ8KECbrqqquUl5cnSfrZz36mwYMH67nnntNdd92lVatWacOGDdq8eXMQ3yoAACJjcAF+BQavvvqqJCktLc3r+oIFCzRp0iRJ0syZM3Xy5Ek9+uijOnr0qAYNGqT169crPDzcrF9WVqZ27f6drEhNTdWSJUv0i1/8Qk8//bSuvfZaLV26VIMGDbrItwUAQDM4+dCngM4xuJRwjgEAXL6+03MMIv+/4Jxj8I8FIXmOAc9KAABYimE0yAjwscmBtr+UERgAAKzFMAKfCgiNZHuTCAwAANZiBGGNQQgHBgEdiQwAAEILGQMAgLU0NEi2ANcIsMYAAIAQwVSCT0wlAAAAExkDAIClGA0NMgKcSmC7IgAAoYKpBJ+YSgAAACYyBgAAa2kwJBsZg+YQGAAArMUwJAW6XTF0AwOmEgAAgImMAQDAUowGQ0aAUwkh8mDiJhEYAACsxWhQ4FMJobtdkakEAIClGA1GUIq/XnnlFSUkJKhTp05KSkrShx9+2ArvLnAEBgAAtLKlS5cqOztbTz31lPbs2aMf/ehHysjIUFlZWVsPrZGQmUo4O99zWnUBn1sBAPhunVadpO9m7v60URPwVMDZ8VZVVXldt9vtstvtjeq/+OKLmjx5sn76059KkubOnav3339fr776qvLy8gIaS7CFTGBw/PhxSdJmrW3jkQAALtbx48flcDhape+OHTvK6XRqszs43xPf+973FBcX53XtmWeeUW5urte12tpa7dq1S0888YTX9fT0dG3ZsiUoYwmmkAkMXC6XysvLFR4eLpvN1uh+VVWV4uLiVF5eroiIiDYYYWjh8wwuPs/g4vMMru/i8zQMQ8ePH5fL5WqV/iWpU6dOKi0tVW1tbVD6Mwyj0fdNU9mCb7/9VvX19YqJifG6HhMTI7fbHZSxBFPIBAbt2rVT9+7dL1gvIiKC/6MIIj7P4OLzDC4+z+Bq7c+ztTIF5+rUqZM6derU6j+nKecHEU0FFpcCFh8CANCKoqKi1L59+0bZgcrKykZZhEsBgQEAAK2oY8eOSkpKUmFhodf1wsJCpaamttGomhcyUwkXYrfb9cwzzzQ5/wP/8XkGF59ncPF5BhefZ+BycnKUmZmpAQMGKCUlRW+88YbKysr0yCOPtPXQGrEZoXyuIwAAl4hXXnlFzz//vCoqKpSYmKg5c+Zo8ODBbT2sRggMAACAiTUGAADARGAAAABMBAYAAMBEYAAAAEyWCQwul8ddXupyc3Nls9m8itPpbOthXTY2bdqkkSNHyuVyyWazaeXKlV73DcNQbm6uXC6XOnfurLS0NO3fv79tBnsZuNDnOWnSpEa/r8nJyW0z2EtcXl6eBg4cqPDwcEVHR+vuu+/WwYMHverw+2kNlggMLqfHXV4OevfurYqKCrPs3bu3rYd02Thx4oT69u2r/Pz8Ju8///zzevHFF5Wfn68dO3bI6XTq9ttvNx8SBm8X+jwlafjw4V6/r2vX8qC1phQXF2vKlCnatm2bCgsLdfr0aaWnp+vEiRNmHX4/LcKwgBtvvNF45JFHvK794Ac/MJ544ok2GtHl65lnnjH69u3b1sMICZKMFStWmK8bGhoMp9NpPPvss+a1U6dOGQ6Hw3jttdfaYISXl/M/T8MwjIkTJxp33XVXm4zncldZWWlIMoqLiw3D4PfTSkI+Y3D2cZfp6ele1y/Vx11eDg4dOiSXy6WEhATdd999+vvf/97WQwoJpaWlcrvdXr+rdrtdt9xyC7+rASgqKlJ0dLR69eqlrKwsVVZWtvWQLgsej0eSFBkZKYnfTysJ+cDgcnvc5aVu0KBBWrhwod5//329+eabcrvdSk1N1ZEjR9p6aJe9s7+P/K4GT0ZGhhYvXqyNGzfqhRde0I4dO3TrrbeqpqamrYd2STMMQzk5Obr55puVmJgoid9PK7HMsxIul8ddXuoyMjLMf/fp00cpKSm69tpr9fbbbysnJ6cNRxY6+F0NnnHjxpn/TkxM1IABAxQfH681a9Zo9OjRbTiyS9vUqVP1ySefaPPmzY3u8fsZ+kI+Y3C5Pe7yctOlSxf16dNHhw4dauuhXPbO7u7gd7X1xMbGKj4+nt9XH6ZNm6bVq1frgw8+UPfu3c3r/H5aR8gHBpfb4y4vNzU1NTpw4IBiY2PbeiiXvYSEBDmdTq/f1draWhUXF/O7GiRHjhxReXk5v69NMAxDU6dO1fLly7Vx40YlJCR43ef30zosMZVwOT3u8lI3Y8YMjRw5Uj169FBlZaV+9atfqaqqShMnTmzroV0Wqqur9dlnn5mvS0tLVVJSosjISPXo0UPZ2dmaPXu2evbsqZ49e2r27Nm64oorNH78+DYc9aXL1+cZGRmp3NxcjRkzRrGxsfr888/15JNPKioqSvfcc08bjvrSNGXKFL377rtatWqVwsPDzcyAw+FQ586dZbPZ+P20ijbdE/Edevnll434+HijY8eORv/+/c0tOPDPuHHjjNjYWKNDhw6Gy+UyRo8ebezfv7+th3XZ+OCDDwxJjcrEiRMNwzizJeyZZ54xnE6nYbfbjcGDBxt79+5t20Ffwnx9nv/85z+N9PR048orrzQ6dOhg9OjRw5g4caJRVlbW1sO+JDX1OUoyFixYYNbh99MaeOwyAAAwhfwaAwAA0HIEBgAAwERgAAAATAQGAADARGAAAABMBAYAAMBEYAAAAEwEBgAAwERgAAAATAQGAADARGAAAABM/z8dxt+i2NKqwQAAAABJRU5ErkJggg==",
      "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": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2855e9596fdf4770b7021cfff78c612d",
       "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": 86,
     "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
}