tug/proto/FTCS.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# environment\n",
    "dim = 2\n",
    "grid_size = {'x': 20, 'y': 20} # to start out\n",
    "domain_size = {'x': 20, 'y': 20} # to start out\n",
    "\n",
    "# boundary_condition\n",
    "# constant boundary condition\n",
    "bc_left = 0\n",
    "bc_right = 0\n",
    "bc_top = 0\n",
    "bc_bottom = 0\n",
    "\n",
    "delta_x = domain_size['x'] / grid_size['x']\n",
    "delta_y = domain_size['y'] / grid_size['y'] \n",
    "\n",
    "alpha_x = np.ones((grid_size['x'],grid_size['y']))\n",
    "alpha_y = np.ones((grid_size['x'],grid_size['y']))\n",
    "\n",
    "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",
    "max_stable_time_step = 0.1\n",
    "\n",
    "C_t = np.zeros((grid_size['x'],grid_size['y']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1078a0550>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD4CAYAAADl7fPiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAAsTAAALEwEAmpwYAAANiklEQVR4nO3df6xk5V3H8ffH5UcTBMuKbPllS+qGhDa6NpvFRjQgLQVC3NY0dYlRVAzYlMQmJgY1KU39p8Yg0dC02dYN1LRAo67dpMuPzWpCm7SUhSwFWpCV0LC3lLXdyhZbi9t+/eOea+5zd+7udc7MnbnT9yu5mXOe55lznpNJPnvOmdnzTVUhSQt+YtITkDRdDAVJDUNBUsNQkNQwFCQ1Tpr0BAY5JafWazht0tOQZtZ/81+8Wj/IoL6pDIXXcBqX5IpJT0OaWQ/X3mX7vHyQ1OgVCkmuSvJMkgNJbhnQf2qSe7v+h5O8oc/+JI3f0KGQZB3wEeBq4GLguiQXLxl2A/Cdqvo54HbgL4fdn6TV0edMYQtwoKqeq6pXgXuArUvGbAXu6pb/AbgiycCbG5KmQ59QOA94YdH6wa5t4JiqOgq8DPz0oI0luTHJviT7/ocf9JiWpD6m5kZjVW2vqs1VtflkTp30dKQfW31CYQ64YNH6+V3bwDFJTgJ+Cvh2j31KGrM+ofAIsDHJhUlOAbYBu5aM2QVc3y2/G/iX8v9qS1Nt6B8vVdXRJDcDDwDrgB1V9VSSDwH7qmoX8HfA3yc5ABxmPjgkTbFM4z/cZ2R9+YtGaXwerr0cqcMDvwmcmhuNkqaDoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCp0adC1AVJ/jXJV5M8leSPBoy5LMnLSfZ3fx/oN11J49an6vRR4I+r6rEkpwOPJtlTVV9dMu7zVXVtj/1IWkVDnylU1YtV9Vi3/F3gaxxbIUrSGjOSewpdNelfBB4e0P3WJI8nuS/Jm46zDcvGSVOgz+UDAEl+EvhH4P1VdWRJ92PA66vqlSTXAP8MbBy0naraDmyH+Ue8952XpOH0OlNIcjLzgfCpqvqnpf1VdaSqXumWdwMnJzmrzz4ljVefbx/CfAWor1XVXy8z5nULpeeTbOn2Zy1JaYr1uXz4ZeC3gSeS7O/a/gz4WYCq+hjz9SPfm+Qo8H1gm7UkpenWp5bkF4CBZacWjbkDuGPYfUhaff6iUVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDU6B0KSZ5P8kRXFm7fgP4k+dskB5J8Jclb+u5T0vj0rvvQubyqvrVM39XM13rYCFwCfLR7lTSFVuPyYSvwyZr3JeC1Sc5Zhf1KGsIoQqGAB5M8muTGAf3nAS8sWj/IgJqTlo2TpsMoLh8uraq5JGcDe5I8XVUP/X83Ytk4aTr0PlOoqrnu9RCwE9iyZMgccMGi9fO7NklTqG8tydOSnL6wDFwJPLlk2C7gd7pvIX4JeLmqXuyzX0nj0/fyYQOwsysXeRLw6aq6P8kfwv+VjtsNXAMcAL4H/F7PfUoao16hUFXPAb8woP1ji5YLeF+f/UhaPf6iUVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUGDoUklzUlYpb+DuS5P1LxlyW5OVFYz7Qe8aSxmroZzRW1TPAJoAk65h/bPvOAUM/X1XXDrsfSatrVJcPVwD/XlVfH9H2JE3IqEJhG3D3Mn1vTfJ4kvuSvGm5DVg2TpoOmX8Ce48NJKcA3wDeVFUvLek7A/hRVb2S5Brgb6pq44m2eUbW1yW5ote8JC3v4drLkTqcQX2jOFO4GnhsaSAAVNWRqnqlW94NnJzkrBHsU9KYjCIUrmOZS4ckr0tXPirJlm5/3x7BPiWNSa8KUV39yLcDNy1qW1wy7t3Ae5McBb4PbKu+1yuSxqr3PYVx8J6CNF7jvqcgaYYYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqdHryUvSgge+sX/FY99x7qaxzUP9eaYgqbGiUEiyI8mhJE8ualufZE+SZ7vXM5d57/XdmGeTXD+qiUsaj5WeKdwJXLWk7RZgb1fHYW+33kiyHrgVuATYAty6XHhImg4rCoWqegg4vKR5K3BXt3wX8M4Bb30HsKeqDlfVd4A9HBsukqZIn3sKG6rqxW75m8CGAWPOA15YtH6wa5M0pUZyo7Gr5dDrWfHWkpSmQ59QeCnJOQDd66EBY+aACxatn9+1HaOqtlfV5qrafDKn9piWpD76hMIuYOHbhOuBzw4Y8wBwZZIzuxuMV3ZtkqbUSr+SvBv4InBRkoNJbgA+DLw9ybPA27p1kmxO8gmAqjoM/AXwSPf3oa5N0pRa0S8aq+q6ZbqOqe1WVfuAP1i0vgPYMdTsJK06f+askfCny7PDnzlLahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhonDIVl6kj+VZKnk3wlyc4kr13mvc8neSLJ/iT7RjhvSWOykjOFOzm21Nse4M1V9fPAvwF/epz3X15Vm6pq83BTlLSaThgKg+pIVtWDVXW0W/0S80VeJM2AUdxT+H3gvmX6CngwyaNJbjzeRiwbJ02HXo94T/LnwFHgU8sMubSq5pKcDexJ8nR35nGMqtoObAc4I+t71aWUNLyhzxSS/C5wLfBbXYHZY1TVXPd6CNgJbBl2f5JWx1ChkOQq4E+AX6+q7y0z5rQkpy8sM19H8slBYyVNj5V8JTmojuQdwOnMXxLsT/Kxbuy5SXZ3b90AfCHJ48CXgc9V1f1jOQpJI5Nlzvwn6oysr0tyTJlKSSPycO3lSB3OoD5/0SipYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIagxbNu6DSea65zPuT3LNMu+9KskzSQ4kuWWUE5c0HsOWjQO4vSsHt6mqdi/tTLIO+AhwNXAxcF2Si/tMVtL4DVU2boW2AAeq6rmqehW4B9g6xHYkraI+9xRu7qpO70hy5oD+84AXFq0f7NoGsmycNB2GDYWPAm8ENgEvArf1nUhVba+qzVW1+WRO7bs5SUMaKhSq6qWq+mFV/Qj4OIPLwc0BFyxaP79rkzTFhi0bd86i1XcxuBzcI8DGJBcmOQXYBuwaZn+SVs8Jq053ZeMuA85KchC4FbgsySbmS80/D9zUjT0X+ERVXVNVR5PcDDwArAN2VNVT4zgISaNj2Tjpx5Bl4yStmKEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqbGSZzTuAK4FDlXVm7u2e4GLuiGvBf6zqjYNeO/zwHeBHwJHq2rzSGYtaWxOGArMl427A/jkQkNV/ebCcpLbgJeP8/7Lq+pbw05Q0uo6YShU1UNJ3jCoL0mA9wC/NuJ5SZqQvvcUfgV4qaqeXaa/gAeTPJrkxuNtyLJx0nRYyeXD8VwH3H2c/kurai7J2cCeJE93BWuPUVXbge0w/4j3nvOSNKShzxSSnAT8BnDvcmOqaq57PQTsZHB5OUlTpM/lw9uAp6vq4KDOJKclOX1hGbiSweXlJE2RE4ZCVzbui8BFSQ4muaHr2saSS4ck5ybZ3a1uAL6Q5HHgy8Dnqur+0U1d0jhYNk76MWTZOEkrZihIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqTGVD5kJcl/AF9f0nwWMIv1I2b1uGB2j20Wjuv1VfUzgzqmMhQGSbJvFitMzepxwewe26we1wIvHyQ1DAVJjbUUCtsnPYExmdXjgtk9tlk9LmAN3VOQtDrW0pmCpFVgKEhqrIlQSHJVkmeSHEhyy6TnMypJnk/yRJL9SfZNej59JNmR5FCSJxe1rU+yJ8mz3euZk5zjMJY5rg8mmes+t/1JrpnkHEdt6kMhyTrgI8DVwMXAdUkunuysRuryqto0A9973wlctaTtFmBvVW0E9nbra82dHHtcALd3n9umqto9oH/NmvpQYL5S9YGqeq6qXgXuAbZOeE5aoqoeAg4vad4K3NUt3wW8czXnNArLHNdMWwuhcB7wwqL1g13bLCjgwSSPJrlx0pMZgw1V9WK3/E3miw7PipuTfKW7vFhzl0XHsxZCYZZdWlVvYf7S6H1JfnXSExqXmv/ue1a+//4o8EZgE/AicNtEZzNiayEU5oALFq2f37WteVU1170eAnYyf6k0S15Kcg5A93powvMZiap6qap+WFU/Aj7OjH1uayEUHgE2JrkwySnANmDXhOfUW5LTkpy+sAxcCTx5/HetObuA67vl64HPTnAuI7MQdJ13MWOf20mTnsCJVNXRJDcDDwDrgB1V9dSEpzUKG4CdSWD+c/h0Vd0/2SkNL8ndwGXAWUkOArcCHwY+k+QG5v8r/HsmN8PhLHNclyXZxPzl0PPATZOa3zj4M2dJjbVw+SBpFRkKkhqGgqSGoSCpYShIahgKkhqGgqTG/wLOTMif1Yv+LQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# initialization\n",
    "C_t[10,10] = 2000\n",
    "plt.imshow(C_t, vmin=0, vmax=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "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": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulating one time step\n",
    "def simulate(C_t): \n",
    "        C_t1 = np.zeros((grid_size['x'],grid_size['y']))\n",
    "\n",
    "        # inner cells\n",
    "        for i in range(1, grid_size['x']-1):\n",
    "                for j in range(1, grid_size['y']-1):\n",
    "                        C_t1[i,j] = C_t[i,j] \\\n",
    "                                + max_stable_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",
    "                                + max_stable_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(0, grid_size['y']):\n",
    "               C_t1[i,0] = C_t[i,0] \\\n",
    "                        +  max_stable_time_step/delta_x**2 * (alpha_interblock(alpha_x[i,1], alpha_x[i,0]) * C_t[i,1] \n",
    "                                - (alpha_interblock(alpha_x[i,1], alpha_x[i,0]) + 2 * alpha_x[i,0]) * C_t[i,0]\n",
    "                                + 2 * alpha_x[i,0] * bc_left)\n",
    "\n",
    "        # right\n",
    "        n = grid_size['x'] - 1 # maximum index in x-direction (columns)\n",
    "        for i in range(0, grid_size['y']):\n",
    "               C_t1[i,n] = C_t[i,n] \\\n",
    "                        +  max_stable_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(0, grid_size['x']):\n",
    "               C_t1[0,j] = C_t[0,j] \\\n",
    "                        +  max_stable_time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_x[0,j]) * C_t[1,j] \n",
    "                                - (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) + 2 * alpha_y[0,j]) * C_t[0,j]\n",
    "                                + 2 * alpha_y[0,j] * bc_top)\n",
    "\n",
    "        # bottom\n",
    "        m = grid_size['y'] - 1 # maximum index in y-direction (rows)\n",
    "        for j in range(0, grid_size['x']):\n",
    "               C_t1[m,j] = C_t[m,j] \\\n",
    "                        +  max_stable_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": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT4AAAD8CAYAAADub8g7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcgklEQVR4nO3df6xc5X3n8ffHPwDhQIF44/CrJGotJIqK01qmUeiuCYEYiw1Jxab2VikkVA5RkILUqktbCSLyT6ouQdslG3IDlqEiDrtJnKDGAVwayUFKCMYyiflVHESEDbFlTPiRX3Dv/e4f59xoMp6593nOnLlz5p7Pyzq6M+c8c84zc+d+/ZznOc/5KiIwM2uTRaOugJnZfHPgM7PWceAzs9Zx4DOz1nHgM7PWceAzs9Zx4DOzkZJ0pqTvSHpC0uOSPlWuP0XSDknPlD9P7vP6K8syz0i6MumYvo7PzEZJ0qnAqRGxW9IJwKPAB4GrgCMR8VlJ1wMnR8T/6HrtKcAuYDUQ5Wv/OCJenu2YbvGZ2UhFxIsRsbt8/BrwJHA6cDlwZ1nsTopg2O39wI6IOFIGux3AurmOuaSGetfuGB0bx7Fs1NUYLY26Agtcy090fsXPeSN+PdC37P0XLouXjkwllX30h79+HPhVx6qJiJjoLifpHcC7gIeBFRHxYrnpp8CKHrs+HXi+4/n+ct2sGhn4jmMZ5+ui+nesIUUT1d9w1iJHvmGK6RFHvpjOKFt/XR+OBwfex+EjUzx8/xlJZZee+uNfRcTq2cpIegvwNeC6iHhVHX+vERGSavsgfKprZhUFUzGdtMxF0lKKoHd3RHy9XH2w7P+b6Qc81OOlB4AzO56fUa6b1UCBT9I6SU9L2ld2PnZvP1bSPeX2h8tmrJktAAFME0nLbFQ07e4AnoyIz3VsuheYGaW9Evhmj5ffD1wi6eRy1PeSct2sKgc+SYuBzwOXAucAGyWd01XsauDliPh94BbgH6sez8yaZzrx3xzeA3wEeK+kPeWyHvgscLGkZ4D3lc+RtFrS7QARcQT4DPBIudxUrpvVIH18a4B9EfFsWZmvUIzCPNFR5nLg0+XjrwK3SlL4GhqzsRcEb+b0VfbbT8RD9B/OO6qzPyJ2AX/V8XwzsDnnmIOc6qaMpvymTERMAq8Ab+21M0mbJO2StOtNfj1AtcxsPgQwRSQtTdOYUd1yaHsC4ESd0rxPysyOMlf/XVMNEvhSRlNmyuyXtAT4HeClAY5pZg0RwNSY9loNcqr7CLBS0jslHQNsoBiF6dQ5KnMF8O/u3zNbOKYTl6ap3OKLiElJ11IMHS8GNkfE45JuAnZFxL0UQ9T/ImkfcIQiOJrZAhAN7b9LMVAfX0RsB7Z3rbuh4/GvgP82yDHmtGhxclEtTi9LxsyJ5P3mzBxZlN4Y17BmpDRh9kjGDIuckwlNJ7ZDck5QUvcJxFRO2bRpYUUdMsoOKALeHM+415zBDTMbN2JqTCeVO/CZWSVBVoO8URz4zKwyt/jMrFWKC5gd+MysRQJ4M8bzBk8OfGZWSSCmxvTOdg58ZlbZdPhU18xaxH18ZtZCYsp9fGbWJsUdmB346pU4DStnGpoWZ/ySli6tf79ZU+aG9L5yDGsqXI6cKWM5ZVOngeVMLZucTC6rRcOZWhapNwat4cLjCPFGZHynG6S5gc/MGm/afXxm1ibF4EY9ZxuSNgOXAYci4txy3T3A2WWRk4CfRcSqHq99DngNmAIm50pjCQ58ZlZZrYMbW4BbgbtmVkTEn//mSNLNFKkr+rkwIg6nHsyBz8wqqXNwIyJ29ks/W6af/DDw3loOhhOKm9kApkJJy4D+FDgYEc/02R7AA5IelbQpZYdu8ZlZJYF4M5JDyHJJuzqeT5QJxlJsBLbOsv2CiDgg6W3ADklPRcTO2XbowGdmlWQObhxOGXToViYp+zPgj/vWI+JA+fOQpG0UOb9nDXyVT3UlnSnpO5KekPS4pE/1KLNW0isd2dFv6LUvMxs/Qdpp7oCnuu8DnoqI/b02Slom6YSZx8AlwN65djpIi28S+OuI2F0e+FFJOyLiia5y342IywY4jpk1VF2DG5K2AmspTon3AzdGxB0UCcq2dpU9Dbg9ItYDK4BtZd6ZJcCXI+K+uY43SJa1F4EXy8evSXoSOB3oDnxmtgBFUNvlLBGxsc/6q3qsewFYXz5+Fjgv93i19PGVw9DvAh7usfndkh4DXgD+JiIe77OPTcAmgOM4Pv3gOZnAcqahLcn4aI5J268yMqeRc/ycKWs509CGNWVtSNPQlDG9LDUjWt40tPTPK95ILpqeEQ6Iqfm7UKMY3GjplDVJbwG+BlwXEa92bd4NnBURr0taD3wDWNlrP+UIzwTAiTplTFOYmLXLuN6IdKBaS1pKEfTujoivd2+PiFcj4vXy8XZgqaTlgxzTzJohENORtjRN5RZfeTX1HcCTEfG5PmXeTnHhYUhaQxFoX6p6TDNrlnFt8Q1yqvse4CPAjyTtKdf9PfC7ABFxG3AF8AlJk8AvgQ2Rk+7ezBqryKvbssAXEQ/B7PekiYhbKSYem9mCI9963szapUgv2dJRXTNrpwi171TXzMzJhsysVYr78bmPz8xaxekl66e0D3RoWdYSp6EBKHUq3JL0ukZG2bzsbQ2YspZjCNPQADSZluVMGZ9BznVamk4vHakZ4eZZcTlLA74jFTQ38JlZo7V6rq6ZtZcTiptZqxS3pfKprpm1jPv4zKxViruz+FTXzFqkmLLmwGdmrTK+Lb7xrLWZNcI0SlrmImmzpEOS9nas+7SkAx1ZGtf3ee06SU9L2ifp+pR6O/CZWSUzo7o1pZfcAqzrsf6WiFhVLtu7N0paDHweuBQ4B9go6Zy5Djb+p7o5swtyZnlkJQZK228szfi4M8rGkoy6ZryvGNLMDeXcizZnNsZUxmyIxPeW8wnkvK+s2RgZv4fkhEcZE2Jm3U19WdZ2lknLcq0B9pXZ1pD0FeBy5sj26BafmVWSmXNjuaRdHcumxMNcK+mH5anwyT22nw483/F8f7luVuPf4jOzkQhgMr3FdzgiVmce4gvAZ8pDfQa4GfhY5j56qiO95HPAa8AUMNn95sqkRP+LIgHwL4CrImL3oMc1s9Eb5qhuRByceSzpS8C/9ih2ADiz4/kZ5bpZ1dXiuzAiDvfZdilFLt2VwPkUUfz8mo5rZqMy5NSRkk6NiBfLpx8C9vYo9giwUtI7KQLeBuC/z7Xv+TjVvRy4q8yu9n1JJ3W9ITMbQ3XeiFTSVmAtRV/gfuBGYK2kVeWhngM+XpY9Dbg9ItZHxKSka4H7gcXA5oh4fK7j1RH4AnhAUgBfjIiJru39Oh9/K/CVnZ2bAI7j+BqqZWbDVleLLyI29lh9R5+yL1B0nc083w4cdanLbOoIfBdExAFJbwN2SHoqInbm7qQMmBMAJ+oU5941a7hW34g0Ig6UPw9J2kZxXU1n4KvU+WhmzRaIyenxvCJuoFpLWibphJnHwCUc3QF5L/CXKvwJ8Ir798wWhrqmrM23QVt8K4BtZV6CJcCXI+I+SdcARMRtFOfe64F9FJezfHTAY5pZE0RLT3XLaSLn9Vh/W8fjAD6ZtWNlTL3JmVq2KCM/wJKcKWOJ+82ZhrY0JzHRsKaspe82R15Snoz3NpkxvS2xXGRMQ1PG9Domc5JkZXwX0mswsFb38ZlZeznwmVmrBGJqTAc3HPjMrLImDlykcOAzs0qirYMbZtZu4cBnZu0y3JsUDJMDn5lV5hafmbVKBExNO/CZWct4VNfMWiXwqe7IKCcD1eKMiy2zyiZmWcuYWpZVNqOuWdPbhpRljawsa0OahJVYB01nTHOcTM+clvO9be492jy4YWYtlPN/WJM48JlZZT7VNbNWKUZ165mrK2kzcBlwKCLOLdf9E/BfgTeAHwMfjYif9Xjtc8yS6bGX8ZxhbGaNEJG2JNgCrOtatwM4NyL+EPgP4O9mef2FEbEqNXevA5+ZVRahpGXu/cRO4EjXugciYrJ8+n2KtBW1cOAzs0qCtKBXBr7lknZ1LJsyD/cx4Nt9q1Jkenw0db/u4zOzyjIGdQ+nnoZ2k/QPwCRwd58i2ZkeK7f4JJ0taU/H8qqk67rKrJX0SkeZG6oez8waJiCmlbRUJekqikGPv4g+eQA6Mz0CM5keZ1W5xRcRTwOrysotpkgZua1H0e9GxGVVj2NmzTXMy1kkrQP+FvgvEfGLPmWWAYsi4rWOTI83zbXvuvr4LgJ+HBE/qWl/ZjYG6hrVlbQV+B5wtqT9kq4GbgVOoDh93SPptrLsaZK2ly9dATwk6THgB8C3IuK+uY5XVx/fBmBrn23vLiv1AvA3EfF4r0Jlp+QmgOM4vqZqDSBnutZQMsINaXpbal0BFg/pf/OcUx8NaXpb4uebM7Us6zszjO/XPKtzrm5EbOyx+o4+ZV+gSFnbN9PjXAZu8Uk6BvgA8P96bN4NnBUR5wH/G/hGv/1ExERErI6I1Ut17KDVMrNhC4ocpClLw9RxqnspsDsiDnZviIhXI+L18vF2YKmk5TUc08waoMYLmOdVHae6G+lzmivp7cDBiAhJaygC7Us1HNPMRm6wEdtRGijwlaMoFwMf71h3DUBE3AZcAXxC0iTwS2BDvyFpMxtDY/rXPFDgi4ifA2/tWndbx+NbKUZmzGyhCd+dxczaqI0tPjNrO7f4zKxtpkddgWoc+Mysmpnr+MaQA5+ZVTau12iMf+DLmc4z4ulHkZNZK+c/0py6ZkxDy6lvDi1K/2vJ+sPKyDanqcRztJzfWcZ3UTnTF3Modb81/W4d+MysdXyqa2Ztk3MPiSZx4DOzakJ5d9ppEAc+M6vOLT4zax0HPjNrHQc+M2uVMb6A2Xl1zawyRdoy536kzZIOSdrbse4USTskPVP+PLnPa68syzwj6cqUejvwmVl1kbjMbQuwrmvd9cCDEbESeLB8/lsknQLcCJxPkVbyxn4BspMDn5lVVleLr0wAfqRr9eXAneXjO4EP9njp+4EdEXEkIl4GdnB0AD2K+/jMZpEzbW88e7sGlN7Ht1zSro7nExExMcdrVkTEi+Xjn1Kkkux2OvB8x/P95bpZOfCZWTXpp7EAhyNideVDFXl7ahtDTjrVne+ORzMbE/X18fVyUNKpAOXPQz3KHADO7Hh+RrluVql9fFuYx45HMxsPmk5bKroXmGksXQl8s0eZ+4FLJJ1cxpZLynWzSgp8893xaGZjoqYWn6StwPeAsyXtl3Q18FngYknPAO8rnyNptaTbASLiCPAZ4JFyualcN6tB+viG1vFoZs2XOmKbIiI29tl0UY+yu4C/6ni+Gdicc7xaBjfq6HiUtAnYBHAcx9dRLTMbthbO3Ki14zEiJiJidUSsXqpjB6iWmc2b4Q5uDM0ggW9oHY9mNh7quoB5vqVezjKvHY9mNgZi6KO6Q5PUxzffHY9mNiYa2JpL4Zkb/Qwhb54y9pl19Jy6ZtwqPCcbWpaMFkDOaVLO55tah6x9NkGkfrh1DcfWs5v55sBnZpU1sf8uhe/OYmat4xafmVU3pi0+Bz4zqyaaOWKbwoHPzKpzi8/M2kSM7+CGA5+ZVefAZ2at0tDpaCkc+MysOg9umFnbuMU3KtMZn3zO9KOcslOJ/+1Np//3qOmMa8tzPoOMb+qwZmtl/bFkvDdNZZRNfXM5v7PU70HmfhvNgc/MWqWh99pL4SlrZlZZHffjk3S2pD0dy6uSrusqs1bSKx1lbhik3m7xmVl1NbT4IuJpYBWApMUUd2nf1qPodyPissGP6MBnZgMYwpS1i4AfR8RPat9zB5/qmlk1qfk2ilbhckm7OpZNffa6AdjaZ9u7JT0m6duS/mCQqrvFZ2aVqFwSHY6I1bPuTzoG+ADwdz027wbOiojXJa0HvgGsTD/8b3OLz8yqqzfL2qXA7og4eNRhIl6NiNfLx9uBpZKWV632nIFP0mZJhyTt7Vj3T5KekvRDSdskndTntc9J+lE5CrOraiXNrJlqzrK2kT6nuZLeLknl4zUUseulqvVOafFtAdZ1rdsBnBsRfwj8B72bpjMujIhVczVzzWwM1dTik7QMuBj4ese6ayRdUz69Atgr6THgn4ENEdUvsZ+zjy8idkp6R9e6Bzqefr+slJm1SY03Io2InwNv7Vp3W8fjW4Fb6zlaPYMbHwPu6bMtgAckBfDFiJjot5NylGcTwHEcn3zwrKCfUXYY049yplQxOaQpTTnT2xYPpws4J3PZ0D6z1LI534OcLHo5U9ZyfmfzrcFVm81AgU/SPwCTwN19ilwQEQckvQ3YIempiNjZq2AZFCcATlx0yph+nGbtMq43Kaj8X7qkq4DLgL/od64dEQfKn4corsReU/V4ZtZA9Y7qzptKgU/SOuBvgQ9ExC/6lFkm6YSZx8AlwN5eZc1sPNU8qjtvUi5n2Qp8Dzhb0n5JV1N0Mp5Acfq6R9JtZdnTJG0vX7oCeKgchfkB8K2IuG8o78LM5l9Q3Ig0ZWmYlFHdjT1W39Gn7AvA+vLxs8B5A9XOzBrLyYbMrJ0c+MysbXIuTWoSBz4zq6ahI7YpHPjMrDL38ZlZ6wzhRqTzopmBLyASp+koZ+rP1FR62ZzsWpNp+w2l370s4z5neenQFqVfupk1bS9Hxm6z+pAypqzpzcm0com/2+L4Od+vnGx3DW5WNbhqs2lm4DOz5mvoxckpHPjMrDoHPjNrE1/AbGatpCbfMmsWDnxmVo2v4zOzNhrXy1mcZc3Mqqsv58asiclU+GdJ+8okZ380SLXd4jOzymoe3LgwIg732XYpRR7dlcD5wBfKn5W4xWdm1QTFxfMpy+AuB+6KwveBkySdWnVn49/iy/lQM2YixGTalf0ASpyRkTMbI+dqfU0vTi+bMXuEjLJZs1JyfmdZM3Pqn21D4gwPADK+M1llM2Ycpc54qmtQIqOPb3nXKexEV/KxuRKTnQ483/F8f7nuxbwaF8Y/8JnZSGRex3d4jtzayYnJ6uBTXTOrJvU0N6GFn5CY7ABwZsfzM8p1laTk3Ngs6ZCkvR3rPi3pQDkCs0fS+j6vXSfp6XIk5vqqlTSzZqoj2VBiYrJ7gb8sR3f/BHglIiqd5kLaqe4WiuRCd3WtvyUi/me/F0laDHweuJjifPwRSfdGxBMV62pmTVNPX+EKYFvZ/7wE+HJE3CfpGoCIuA3YTpHPZx/wC+CjgxwwJdnQTknvqLDvNcC+MukQkr5CMTLjwGe2QNRxOUu/xGRlwJt5HMAnBz9aYZA+vmvLCwk3Szq5x/Z+ozA9SdokaZekXW/y6wGqZWbzIoCpSFsapmrg+wLwe8AqiuHkmwetSERMRMTqiFi9lGMH3Z2ZzYNxTShe6XKWiDg481jSl4B/7VGs1lEYM2ugJt8dehaVWnxdV0x/iKNHYAAeAVZKeqekY4ANFCMzZrZALNgWn6StwFqKK6/3AzcCayWtojjLfw74eFn2NOD2iFgfEZOSrgXuBxYDmyPi8WG8CTMbgYV8W6qI2Nhj9R19yr5AMeQ883w7xTD08GRMacqahrYoY7pW6j6zpqFlTNXKSXKTMw0t4zPISo6UISvhUVZiosTPLOM7k/P9Glbiq/kkQA0cuEjhKWtmVlnWvOsGceAzs2oW8qmumVlvtd1yat458JlZZU0csU3hwGdm1bnFZ2atEh7VNbM2Gs+458BnZtX5chYzax8HPjNrlQCaOalkTs0NfJH2iUZOZq1FGdmq3kguihIzW0XONKXJ0WdO06IGpGTJmZKYM7UrNRvZsLKh5UyFy5q2N3+RSIRPdc2shRo6j3guDfgv3czG0sypbsoyC0lnSvqOpCckPS7pUz3KrJX0SkeCsxsGqbpbfGZWWU2nupPAX0fE7jLb2qOSdvRITPbdiLisjgM68JlZdTUEvjJN5Ivl49ckPUmRn2doicl8qmtmFdWXUHxGmdHxXcDDPTa/W9Jjkr4t6Q8GqblbfGZWzUyWtTTLJe3qeD4REROdBSS9BfgacF1EvNr1+t3AWRHxuqT1wDeAlZXqjQOfmQ0go4/vcESs7rsfaSlF0Ls7Ir7evb0zEEbEdkn/R9LyiDicW2dIy7mxGbgMOBQR55br7gHOLoucBPwsIlb1eO1zwGvAFDA52xs3szFUQx+figtR7wCejIjP9SnzduBgRISkNRTddC9VPWZKi28LcCtw18yKiPjzjgrdDLwyy+svrBqVzazBgvQLwWf3HuAjwI8k7SnX/T3wuwARcRtwBfAJSZPAL4ENEdWjbkqyoZ1lh+NRykj9YeC9VStgZuOqnjswR8RDzJGvKiJupWiA1WLQPr4/pWh+PtNnewAPSArgi92dmZ0kbQI2ARzH8ckfaNY0sAw5Wc6S65AzXWxx+pS1rK9eRua0sZPR+khuLAwpG1rONLSs7/h8TyFr6ZS1jcDWWbZfEBEHJL0N2CHpqYjY2atgGRQnAE7UKeP5aZq1SQA584gbpPJ1fJKWAH8G3NOvTEQcKH8eArYBa6oez8yaJoqbIqQsDTPIBczvA56KiP29NkpaVk4/QdIy4BJg7wDHM7OmqfkC5vkyZ+CTtBX4HnC2pP2Sri43baDrNFfSaZK2l09XAA9Jegz4AfCtiLivvqqb2UjNjOqmLA2TMqq7sc/6q3qsewFYXz5+FjhvwPqZWZM1sDWXwjM3zKw6Bz4za5WIvMt9GsSBz8yqc4vPzFrHgc/M2qWZI7Ypxj/wTWdktsq4kDKm6r9HqzKmizXi66Qh3ae2gRe09hPD+sPO+Qya2qqKvL+pJhn/wGdmozOmU9Yc+MysmoixTS/pwGdm1TX1NHwODnxmVlm4xWdm7dLMGxCkcOAzs2rqu/X8vHPgM7NKguHdAX3YnFDczKqJ+m5EKmmdpKcl7ZN0fY/tx0q6p9z+cL88QKkc+MysspiOpGU2khYDnwcuBc4BNko6p6vY1cDLEfH7wC3APw5Sbwc+M6uunhbfGmBfRDwbEW8AXwEu7ypzOXBn+firwEVllsdKGtnH9xovH/63+OpPulYvBwbLzzvqftje3SGDv6/mWqjvbSG8r7MG3cFrvHz/v8VXlycWP07Sro7nEx1ZF08Hnu/Yth84v+v1vykTEZOSXgHeSsXfQyMDX0T8p+51knZFxOpR1GeYFur7goX73hbq+8oVEetGXYeqfKprZqN2ADiz4/kZ5bqeZcoMj78DvFT1gA58ZjZqjwArJb1T0jEUiczu7SpzL3Bl+fgK4N8jOSv80Rp5qtvHxNxFxtJCfV+wcN/bQn1fI1H22V0L3A8sBjZHxOOSbgJ2RcS9wB3Av0jaBxyhCI6VaYCgaWY2lnyqa2at48BnZq0zFoFvruks40rSc5J+JGlP1zVOY0fSZkmHJO3tWHeKpB2Snil/njzKOlbR5319WtKB8ve2R9L6UdbR8jU+8CVOZxlnF0bEqgVwXdgWoPu6ruuBByNiJfBg+XzcbOHo9wVwS/l7WxUR2+e5Tjagxgc+0qaz2IhFxE6K0bZOndOM7gQ+OJ91qkOf92VjbhwCX6/pLKePqC51C+ABSY9K2jTqygzBioh4sXz8U2DFKCtTs2sl/bA8FR67U/i2G4fAt5BdEBF/RHEa/0lJ/3nUFRqW8mLThXLt1BeA3wNWAS8CN4+0NpZtHAJfynSWsRQRB8qfh4BtFKf1C8lBSacClD8Pjbg+tYiIgxExFUVS2S+x8H5vC944BL6U6SxjR9IySSfMPAYuAfbO/qqx0znN6ErgmyOsS21mgnnpQyy839uC1/gpa/2ms4y4WnVYAWwrbym2BPhyRNw32ipVJ2krsBZYLmk/cCPwWeD/Sroa+Anw4dHVsJo+72utpFUUp+7PAR8fVf2sGk9ZM7PWGYdTXTOzWjnwmVnrOPCZWes48JlZ6zjwmVnrOPCZWes48JlZ6/x/bLJC8XAYli8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "records = []\n",
    "\n",
    "for i in range(10):\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=20)\n",
    "plt.colorbar()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update(w = 1):\n",
    "    fig = plt.figure(figsize = (15,10))\n",
    "    y = records[w]\n",
    "    plt.imshow(y)\n",
    "    \n",
    "interact(update, w = IntSlider(min=0, max = 99, 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.10.1"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}