321 lines
47 KiB
Plaintext
321 lines
47 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 87,
|
|
"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": 88,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# set up environment\n",
|
|
"dim = 2\n",
|
|
"grid_size = {'x': 20, 'y': 20} # to start out\n",
|
|
"domain_size = {'x': 20, 'y': 20} # to start out\n",
|
|
"\n",
|
|
"# boundary_condition\n",
|
|
"# constant boundary condition\n",
|
|
"bc_left = 0\n",
|
|
"bc_right = 0\n",
|
|
"bc_top = 0\n",
|
|
"bc_bottom = 0\n",
|
|
"\n",
|
|
"delta_x = domain_size['x'] / grid_size['x']\n",
|
|
"delta_y = domain_size['y'] / grid_size['y'] \n",
|
|
"\n",
|
|
"alpha_x = np.ones((grid_size['x'],grid_size['y']))\n",
|
|
"alpha_y = np.ones((grid_size['x'],grid_size['y']))\n",
|
|
"\n",
|
|
"# time dimension\n",
|
|
"max_stable_time_step_x = delta_x**2 / (2 * np.max(alpha_x))\n",
|
|
"max_stable_time_step_y = delta_y**2 / (2 * np.max(alpha_y))\n",
|
|
"max_stable_time_step = max_stable_time_step_x if max_stable_time_step_x <= max_stable_time_step_y else max_stable_time_step_y\n",
|
|
"\n",
|
|
"time_step = 0.1\n",
|
|
"if time_step > max_stable_time_step:\n",
|
|
" print('Warning: Time step more than maximum stable time step!')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 102,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAGdCAYAAABKG5eZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkA0lEQVR4nO3dfXSU5Z3/8c/EhAmyyViUJDMQQuAglIdDISABlQcpwVARKhXUPRDW1pYttWLKKcTWI+4fDbbV5SAoa8uD1C5yuiGQ3bCVcMyDlMCCJK61iHFNSVaScuBIBnAZErh+f/jL6JiZgYGZJFd8v865z/G+7+915TtXIp/cM/dMHMYYIwAALBHX1Q0AABAJggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYJX4rm4gWq5cuaKTJ08qKSlJDoejq9sBAETAGKNz587J4/EoLi78NVWPCa6TJ08qPT29q9sAANyAxsZGDRgwIGxNjwmupKQkSdKJo4OU/Hc39gzot28fHY2WAADXqE2t2q89/n/Lw+kxwdX+9GDy38UpOenGgivekRCNlgAA1+r/f2rutbzUw80ZAACrEFwAAKvELLheeuklZWZmKjExUVlZWXrrrbfC1ldWViorK0uJiYkaPHiwNm7cGKvWAAAWi0lw7dixQ8uXL9fPfvYz1dTU6O6771Zubq4aGhqC1tfX12v27Nm6++67VVNTo6eeeko//vGPVVRUFIv2AAAWc8TiD0lOnDhR48aN08svv+w/9vWvf13z5s1TYWFhh/qVK1eqpKREx44d8x9bunSp3nnnHVVXV1/T1/R6vXK5XPrkg8E3fHPGLM83bmg8ACAybaZVFdqtlpYWJScnh62N+hXXpUuX9PbbbysnJyfgeE5Ojg4cOBB0THV1dYf6WbNm6ciRI2ptbQ06xufzyev1BmwAgJ4v6sF1+vRpXb58WampqQHHU1NT1dzcHHRMc3Nz0Pq2tjadPn066JjCwkK5XC7/xpuPAeCrIWY3Z3z5XnxjTNj784PVBzverqCgQC0tLf6tsbHxBjsGANgg6m9Avu2223TTTTd1uLo6depUh6uqdmlpaUHr4+PjdeuttwYd43Q65XQ6o9M0AMAaUb/i6tWrl7KyslRWVhZwvKysTJMnTw46ZtKkSR3q9+7dq/HjxyshgU+xAAB8LiZPFebn5+u3v/2tNm/erGPHjunJJ59UQ0ODli5dKumzp/kWL17sr1+6dKlOnDih/Px8HTt2TJs3b9amTZu0YsWKWLQHALBYTD6rcOHChTpz5oz+6Z/+SU1NTRo1apT27NmjjIwMSVJTU1PAe7oyMzO1Z88ePfnkk9qwYYM8Ho/WrVun+fPnx6I9AIDFYvI+rq7A+7gAwF5d+j4uAABiqcf8WZN23759NH+WBAB6MK64AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAVol6cBUWFmrChAlKSkpSSkqK5s2bp+PHj4cdU1FRIYfD0WF7//33o90eAMByUQ+uyspKLVu2TAcPHlRZWZna2tqUk5OjCxcuXHXs8ePH1dTU5N+GDh0a7fYAAJaLj/aEf/zjHwP2t2zZopSUFL399tuaMmVK2LEpKSm65ZZbot0SAKAHiflrXC0tLZKkvn37XrV27NixcrvdmjFjhsrLy8PW+nw+eb3egA0A0PPFNLiMMcrPz9ddd92lUaNGhaxzu9165ZVXVFRUpJ07d2rYsGGaMWOGqqqqQo4pLCyUy+Xyb+np6bF4CACAbsZhjDGxmnzZsmUqLS3V/v37NWDAgIjGzpkzRw6HQyUlJUHP+3w++Xw+/77X61V6erqmaa7iHQk31DcAoHO1mVZVaLdaWlqUnJwctjZmV1yPP/64SkpKVF5eHnFoSVJ2drbq6upCnnc6nUpOTg7YAAA9X9RvzjDG6PHHH1dxcbEqKiqUmZl5XfPU1NTI7XZHuTsAgO2iHlzLli3Tv/7rv2r37t1KSkpSc3OzJMnlcql3796SpIKCAn388cfatm2bJGnt2rUaNGiQRo4cqUuXLum1115TUVGRioqKot0eAMByUQ+ul19+WZI0bdq0gONbtmzRkiVLJElNTU1qaGjwn7t06ZJWrFihjz/+WL1799bIkSNVWlqq2bNnR7s9AIDlYnpzRmfyer1yuVzcnAEAFuoWN2cAABALBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCpRD67Vq1fL4XAEbGlpaWHHVFZWKisrS4mJiRo8eLA2btwY7bYAAD1EfCwmHTlypPbt2+ffv+mmm0LW1tfXa/bs2Xrsscf02muv6U9/+pN++MMfql+/fpo/f34s2gMAWCwmwRUfH3/Vq6x2Gzdu1MCBA7V27VpJ0te//nUdOXJEv/71rwkuAEAHMXmNq66uTh6PR5mZmXrooYf00Ucfhaytrq5WTk5OwLFZs2bpyJEjam1tDTnO5/PJ6/UGbACAni/qwTVx4kRt27ZNb7zxhn7zm9+oublZkydP1pkzZ4LWNzc3KzU1NeBYamqq2tradPr06ZBfp7CwUC6Xy7+lp6dH9XEAALqnqAdXbm6u5s+fr9GjR+ub3/ymSktLJUmvvvpqyDEOhyNg3xgT9PgXFRQUqKWlxb81NjZGoXsAQHcXk9e4vqhPnz4aPXq06urqgp5PS0tTc3NzwLFTp04pPj5et956a8h5nU6nnE5nVHsFAHR/MX8fl8/n07Fjx+R2u4OenzRpksrKygKO7d27V+PHj1dCQkKs2wMAWCbqwbVixQpVVlaqvr5ehw4d0ne+8x15vV7l5eVJ+uwpvsWLF/vrly5dqhMnTig/P1/Hjh3T5s2btWnTJq1YsSLarQEAeoCoP1X4v//7v3r44Yd1+vRp9evXT9nZ2Tp48KAyMjIkSU1NTWpoaPDXZ2Zmas+ePXryySe1YcMGeTwerVu3jlvhAQBBOUz7nRCW83q9crlcmqa5infwFCMA2KTNtKpCu9XS0qLk5OSwtXxWIQDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqUQ+uQYMGyeFwdNiWLVsWtL6ioiJo/fvvvx/t1gAAPUB8tCc8fPiwLl++7N//85//rJkzZ+rBBx8MO+748eNKTk727/fr1y/arQEAeoCoB9eXA2fNmjUaMmSIpk6dGnZcSkqKbrnllmi3AwDoYWL6GtelS5f02muv6dFHH5XD4QhbO3bsWLndbs2YMUPl5eWxbAsAYLGoX3F90a5du3T27FktWbIkZI3b7dYrr7yirKws+Xw+/e53v9OMGTNUUVGhKVOmhBzn8/nk8/n8+16vN5qtAwC6KYcxxsRq8lmzZqlXr17693//94jGzZkzRw6HQyUlJSFrVq9erWeffbbD8Wmaq3hHQsS9AgC6TptpVYV2q6WlJeB+h2Bi9lThiRMntG/fPn3ve9+LeGx2drbq6urC1hQUFKilpcW/NTY2Xm+rAACLxOypwi1btiglJUXf+ta3Ih5bU1Mjt9sdtsbpdMrpdF5vewAAS8UkuK5cuaItW7YoLy9P8fGBX6KgoEAff/yxtm3bJklau3atBg0apJEjR/pv5igqKlJRUVEsWgMAWC4mwbVv3z41NDTo0Ucf7XCuqalJDQ0N/v1Lly5pxYoV+vjjj9W7d2+NHDlSpaWlmj17dixaAwBYLqY3Z3Qmr9crl8vFzRkAYKFucXMGAACxQHABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArBJxcFVVVWnOnDnyeDxyOBzatWtXwHljjFavXi2Px6PevXtr2rRpeu+99646b1FRkUaMGCGn06kRI0aouLg40tYAAF8BEQfXhQsXNGbMGK1fvz7o+V/+8pd64YUXtH79eh0+fFhpaWmaOXOmzp07F3LO6upqLVy4UIsWLdI777yjRYsWacGCBTp06FCk7QEAejiHMcZc92CHQ8XFxZo3b56kz662PB6Pli9frpUrV0qSfD6fUlNT9dxzz+kHP/hB0HkWLlwor9er//zP//Qfu/fee/W1r31N27dvv6ZevF6vXC6Xpmmu4h0J1/uQAABdoM20qkK71dLSouTk5LC1UX2Nq76+Xs3NzcrJyfEfczqdmjp1qg4cOBByXHV1dcAYSZo1a1bYMT6fT16vN2ADAPR8UQ2u5uZmSVJqamrA8dTUVP+5UOMiHVNYWCiXy+Xf0tPTb6BzAIAtYnJXocPhCNg3xnQ4dqNjCgoK1NLS4t8aGxuvv2EAgDXiozlZWlqapM+uoNxut//4qVOnOlxRfXncl6+urjbG6XTK6XTeYMcAANtE9YorMzNTaWlpKisr8x+7dOmSKisrNXny5JDjJk2aFDBGkvbu3Rt2DADgqyniK67z58/rww8/9O/X19ertrZWffv21cCBA7V8+XL94he/0NChQzV06FD94he/0M0336xHHnnEP2bx4sXq37+/CgsLJUlPPPGEpkyZoueee05z587V7t27tW/fPu3fvz8KDxEA0JNEHFxHjhzR9OnT/fv5+fmSpLy8PG3dulU//elP9X//93/64Q9/qE8++UQTJ07U3r17lZSU5B/T0NCguLjPL/YmT56s119/XT//+c/19NNPa8iQIdqxY4cmTpx4I48NANAD3dD7uLoT3scFAPbqsvdxAQAQawQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqEQdXVVWV5syZI4/HI4fDoV27dvnPtba2auXKlRo9erT69Okjj8ejxYsX6+TJk2Hn3Lp1qxwOR4ft4sWLET8gAEDPFnFwXbhwQWPGjNH69es7nPv000919OhRPf300zp69Kh27typDz74QPfff/9V501OTlZTU1PAlpiYGGl7AIAeLj7SAbm5ucrNzQ16zuVyqaysLODYiy++qDvuuEMNDQ0aOHBgyHkdDofS0tIibQcA8BUT89e4Wlpa5HA4dMstt4StO3/+vDIyMjRgwADdd999qqmpCVvv8/nk9XoDNgBAzxfT4Lp48aJWrVqlRx55RMnJySHrhg8frq1bt6qkpETbt29XYmKi7rzzTtXV1YUcU1hYKJfL5d/S09Nj8RAAAN2Mwxhjrnuww6Hi4mLNmzevw7nW1lY9+OCDamhoUEVFRdjg+rIrV65o3LhxmjJlitatWxe0xufzyefz+fe9Xq/S09M1TXMV70iI+LEAALpOm2lVhXarpaXlqnkR8Wtc16K1tVULFixQfX293nzzzYhCS5Li4uI0YcKEsFdcTqdTTqfzRlsFAFgm6k8VtodWXV2d9u3bp1tvvTXiOYwxqq2tldvtjnZ7AADLRXzFdf78eX344Yf+/fr6etXW1qpv377yeDz6zne+o6NHj+o//uM/dPnyZTU3N0uS+vbtq169ekmSFi9erP79+6uwsFCS9Oyzzyo7O1tDhw6V1+vVunXrVFtbqw0bNkTjMQIAepCIg+vIkSOaPn26fz8/P1+SlJeXp9WrV6ukpESS9I1vfCNgXHl5uaZNmyZJamhoUFzc5xd7Z8+e1fe//301NzfL5XJp7Nixqqqq0h133BFpewCAHu6Gbs7oTrxer1wuFzdnAICFIrk5g88qBABYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFgl4uCqqqrSnDlz5PF45HA4tGvXroDzS5YskcPhCNiys7OvOm9RUZFGjBghp9OpESNGqLi4ONLWAABfAREH14ULFzRmzBitX78+ZM29996rpqYm/7Znz56wc1ZXV2vhwoVatGiR3nnnHS1atEgLFizQoUOHIm0PANDDxUc6IDc3V7m5uWFrnE6n0tLSrnnOtWvXaubMmSooKJAkFRQUqLKyUmvXrtX27dsjbREA0IPF5DWuiooKpaSk6Pbbb9djjz2mU6dOha2vrq5WTk5OwLFZs2bpwIEDIcf4fD55vd6ADQDQ80U9uHJzc/X73/9eb775pp5//nkdPnxY99xzj3w+X8gxzc3NSk1NDTiWmpqq5ubmkGMKCwvlcrn8W3p6etQeAwCg+4r4qcKrWbhwof+/R40apfHjxysjI0OlpaV64IEHQo5zOBwB+8aYDse+qKCgQPn5+f59r9dLeAHAV0DUg+vL3G63MjIyVFdXF7ImLS2tw9XVqVOnOlyFfZHT6ZTT6YxanwAAO8T8fVxnzpxRY2Oj3G53yJpJkyaprKws4NjevXs1efLkWLcHALBMxFdc58+f14cffujfr6+vV21trfr27au+fftq9erVmj9/vtxut/7617/qqaee0m233aZvf/vb/jGLFy9W//79VVhYKEl64oknNGXKFD333HOaO3eudu/erX379mn//v1ReIgAgJ4k4uA6cuSIpk+f7t9vf50pLy9PL7/8st59911t27ZNZ8+eldvt1vTp07Vjxw4lJSX5xzQ0NCgu7vOLvcmTJ+v111/Xz3/+cz399NMaMmSIduzYoYkTJ97IYwMA9EAOY4zp6iaiwev1yuVyaZrmKt6R0NXtAAAi0GZaVaHdamlpUXJycthaPqsQAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYJWIg6uqqkpz5syRx+ORw+HQrl27As47HI6g269+9auQc27dujXomIsXL0b8gAAAPVvEwXXhwgWNGTNG69evD3q+qakpYNu8ebMcDofmz58fdt7k5OQOYxMTEyNtDwDQw8VHOiA3N1e5ubkhz6elpQXs7969W9OnT9fgwYPDzutwODqMBQDgy2L6Gtff/vY3lZaW6rvf/e5Va8+fP6+MjAwNGDBA9913n2pqasLW+3w+eb3egA0A0PPFNLheffVVJSUl6YEHHghbN3z4cG3dulUlJSXavn27EhMTdeedd6quri7kmMLCQrlcLv+Wnp4e7fYBAN2Qwxhjrnuww6Hi4mLNmzcv6Pnhw4dr5syZevHFFyOa98qVKxo3bpymTJmidevWBa3x+Xzy+Xz+fa/Xq/T0dE3TXMU7EiL6egCArtVmWlWh3WppaVFycnLY2ohf47pWb731lo4fP64dO3ZEPDYuLk4TJkwIe8XldDrldDpvpEUAgIVi9lThpk2blJWVpTFjxkQ81hij2tpaud3uGHQGALBZxFdc58+f14cffujfr6+vV21trfr27auBAwdK+uxpuz/84Q96/vnng86xePFi9e/fX4WFhZKkZ599VtnZ2Ro6dKi8Xq/WrVun2tpabdiw4XoeEwCgB4s4uI4cOaLp06f79/Pz8yVJeXl52rp1qyTp9ddflzFGDz/8cNA5GhoaFBf3+cXe2bNn9f3vf1/Nzc1yuVwaO3asqqqqdMcdd0TaHgCgh7uhmzO6E6/XK5fLxc0ZAGChSG7O4LMKAQBWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFaJ7+oGosUYI0lqU6tkurgZAEBE2tQq6fN/y8PpMcF17tw5SdJ+7eniTgAA1+vcuXNyuVxhaxzmWuLNAleuXNHJkyeVlJQkh8MRtMbr9So9PV2NjY1KTk7u5A6vH313Plt7p+/ORd/RY4zRuXPn5PF4FBcX/lWsHnPFFRcXpwEDBlxTbXJycrf5ZkWCvjufrb3Td+ei7+i42pVWO27OAABYheACAFjlKxVcTqdTzzzzjJxOZ1e3EhH67ny29k7fnYu+u0aPuTkDAPDV8JW64gIA2I/gAgBYheACAFiF4AIAWKXHBddLL72kzMxMJSYmKisrS2+99VbY+srKSmVlZSkxMVGDBw/Wxo0bO6nTzxQWFmrChAlKSkpSSkqK5s2bp+PHj4cdU1FRIYfD0WF7//33O6lrafXq1R2+flpaWtgxXb3W7QYNGhR0/ZYtWxa0vqvWu6qqSnPmzJHH45HD4dCuXbsCzhtjtHr1ank8HvXu3VvTpk3Te++9d9V5i4qKNGLECDmdTo0YMULFxcWd1ndra6tWrlyp0aNHq0+fPvJ4PFq8eLFOnjwZds6tW7cG/R5cvHixU/qWpCVLlnT4+tnZ2VedtyvXW1LQdXM4HPrVr34Vcs7OWO8b0aOCa8eOHVq+fLl+9rOfqaamRnfffbdyc3PV0NAQtL6+vl6zZ8/W3XffrZqaGj311FP68Y9/rKKiok7rubKyUsuWLdPBgwdVVlamtrY25eTk6MKFC1cde/z4cTU1Nfm3oUOHdkLHnxs5cmTA13/33XdD1naHtW53+PDhgL7LysokSQ8++GDYcZ293hcuXNCYMWO0fv36oOd/+ctf6oUXXtD69et1+PBhpaWlaebMmf7P7QymurpaCxcu1KJFi/TOO+9o0aJFWrBggQ4dOtQpfX/66ac6evSonn76aR09elQ7d+7UBx98oPvvv/+q8yYnJwesf1NTkxITEzul73b33ntvwNffsyf8Z6N29XpL6rBmmzdvlsPh0Pz588POG+v1viGmB7njjjvM0qVLA44NHz7crFq1Kmj9T3/6UzN8+PCAYz/4wQ9MdnZ2zHq8mlOnThlJprKyMmRNeXm5kWQ++eSTzmvsS5555hkzZsyYa67vjmvd7oknnjBDhgwxV65cCXq+O6y3JFNcXOzfv3LliklLSzNr1qzxH7t48aJxuVxm48aNIedZsGCBuffeewOOzZo1yzz00ENR79mYjn0H81//9V9Gkjlx4kTImi1bthiXyxXd5sII1ndeXp6ZO3duRPN0x/WeO3euueeee8LWdPZ6R6rHXHFdunRJb7/9tnJycgKO5+Tk6MCBA0HHVFdXd6ifNWuWjhw5otbW1pj1Gk5LS4skqW/fvletHTt2rNxut2bMmKHy8vJYt9ZBXV2dPB6PMjMz9dBDD+mjjz4KWdsd11r67Ofmtdde06OPPhryw5nbdfV6f1F9fb2am5sD1tTpdGrq1Kkhf96l0N+HcGNiraWlRQ6HQ7fcckvYuvPnzysjI0MDBgzQfffdp5qams5p8AsqKiqUkpKi22+/XY899phOnToVtr67rfff/vY3lZaW6rvf/e5Va7vDeofSY4Lr9OnTunz5slJTUwOOp6amqrm5OeiY5ubmoPVtbW06ffp0zHoNxRij/Px83XXXXRo1alTIOrfbrVdeeUVFRUXauXOnhg0bphkzZqiqqqrTep04caK2bdumN954Q7/5zW/U3NysyZMn68yZM0Hru9tat9u1a5fOnj2rJUuWhKzpDuv9Ze0/05H8vLePi3RMLF28eFGrVq3SI488EvbDXocPH66tW7eqpKRE27dvV2Jiou68807V1dV1Wq+5ubn6/e9/rzfffFPPP/+8Dh8+rHvuuUc+ny/kmO623q+++qqSkpL0wAMPhK3rDusdTo/5dPh2X/6t2RgT9jfpYPXBjneGH/3oR/rv//5v7d+/P2zdsGHDNGzYMP/+pEmT1NjYqF//+teaMmVKrNuU9Nn/xO1Gjx6tSZMmaciQIXr11VeVn58fdEx3Wut2mzZtUm5urjweT8ia7rDeoUT68369Y2KhtbVVDz30kK5cuaKXXnopbG12dnbAjRB33nmnxo0bpxdffFHr1q2LdauSpIULF/r/e9SoURo/frwyMjJUWloaNgi6y3pL0ubNm/X3f//3V32tqjusdzg95orrtttu00033dThN5lTp051+I2nXVpaWtD6+Ph43XrrrTHrNZjHH39cJSUlKi8vv+Y/z/JF2dnZXfrbUJ8+fTR69OiQPXSntW534sQJ7du3T9/73vciHtvV691+B2ckP+/t4yIdEwutra1asGCB6uvrVVZWFvGf1oiLi9OECRO69HvgdruVkZERtofust6S9NZbb+n48ePX9fPeHdb7i3pMcPXq1UtZWVn+O8TalZWVafLkyUHHTJo0qUP93r17NX78eCUkJMSs1y8yxuhHP/qRdu7cqTfffFOZmZnXNU9NTY3cbneUu7t2Pp9Px44dC9lDd1jrL9uyZYtSUlL0rW99K+KxXb3emZmZSktLC1jTS5cuqbKyMuTPuxT6+xBuTLS1h1ZdXZ327dt3Xb+4GGNUW1vbpd+DM2fOqLGxMWwP3WG9223atElZWVkaM2ZMxGO7w3oH6Kq7QmLh9ddfNwkJCWbTpk3mL3/5i1m+fLnp06eP+etf/2qMMWbVqlVm0aJF/vqPPvrI3HzzzebJJ580f/nLX8ymTZtMQkKC+bd/+7dO6/kf//EfjcvlMhUVFaapqcm/ffrpp/6aL/f9z//8z6a4uNh88MEH5s9//rNZtWqVkWSKioo6re+f/OQnpqKiwnz00Ufm4MGD5r777jNJSUndeq2/6PLly2bgwIFm5cqVHc51l/U+d+6cqampMTU1NUaSeeGFF0xNTY3/7rs1a9YYl8tldu7cad59913z8MMPG7fbbbxer3+ORYsWBdxV+6c//cncdNNNZs2aNebYsWNmzZo1Jj4+3hw8eLBT+m5tbTX333+/GTBggKmtrQ34mff5fCH7Xr16tfnjH/9o/ud//sfU1NSYf/iHfzDx8fHm0KFDndL3uXPnzE9+8hNz4MABU19fb8rLy82kSZNM//79u/V6t2tpaTE333yzefnll4PO0RXrfSN6VHAZY8yGDRtMRkaG6dWrlxk3blzAbeV5eXlm6tSpAfUVFRVm7NixplevXmbQoEEhv7GxIinotmXLlpB9P/fcc2bIkCEmMTHRfO1rXzN33XWXKS0t7dS+Fy5caNxut0lISDAej8c88MAD5r333gvZszFdv9Zf9MYbbxhJ5vjx4x3OdZf1br8N/8tbXl6eMeazW+KfeeYZk5aWZpxOp5kyZYp59913A+aYOnWqv77dH/7wBzNs2DCTkJBghg8fHvUADtd3fX19yJ/58vLykH0vX77cDBw40PTq1cv069fP5OTkmAMHDnRa359++qnJyckx/fr1MwkJCWbgwIEmLy/PNDQ0BMzR3da73b/8y7+Y3r17m7NnzwadoyvW+0bwZ00AAFbpMa9xAQC+GgguAIBVCC4AgFUILgCAVQguAIBVCC4AgFUILgCAVQguAIBVCC4AgFUILgCAVQguAIBVCC4AgFX+H4dDK9QkBqc7AAAAAElFTkSuQmCC",
|
|
"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[0,0] = 2000\n",
|
|
"\n",
|
|
"\n",
|
|
"plt.imshow(C_t, vmin=0, vmax=2000)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 90,
|
|
"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": 99,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# simulate one time step\n",
|
|
"def simulate(C_t): \n",
|
|
" C_t1 = np.zeros((grid_size['x'],grid_size['y']))\n",
|
|
"\n",
|
|
" # inner cells\n",
|
|
" for i in range(1, grid_size['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"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 103,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAAGiCAYAAACcWg7FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCv0lEQVR4nO3df1hUZd4/8PeoMCgLk4gwM4lEXkolxCKYgJWYBdLij7RVsy9iuZhPZkvIZVpbsc/uE9au6ZWmmYvir9JnH0VtdVFYBSLQ/EXrr1jcSNAYKS+YEVsBmfv7h8tZR4YD48wAznm/nuu+4pzzOffcc3QfPt6/jkoIIUBEREQEoFd3N4CIiIh6DiYGREREJGFiQERERBImBkRERCRhYkBEREQSJgZEREQkYWJAREREEiYGREREJGFiQERERBImBkRERCRhYkBEROREmZmZGDlyJLy8vODn54fJkyejvLzcIkYIgYyMDOj1evTt2xexsbE4c+aMRUxjYyMWLFgAX19feHp6YuLEibh48aJFTF1dHZKSkqDRaKDRaJCUlIT6+nqb2svEgIiIyIkKCwsxf/58HD58GHl5ebhx4wbi4uJw7do1Keb999/HBx98gFWrVuHo0aPQarV46qmncPXqVSkmNTUVOTk52LZtG4qLi9HQ0IDExES0tLRIMTNnzkRZWRlyc3ORm5uLsrIyJCUl2dZgQURERF2mtrZWABCFhYVCCCHMZrPQarVi6dKlUsz169eFRqMRH3/8sRBCiPr6euHm5ia2bdsmxVy6dEn06tVL5ObmCiGEOHv2rAAgDh8+LMWUlpYKAOKbb77pdPv62JcH9Rxmsxnff/89vLy8oFKpurs5RERkAyEErl69Cr1ej169nNeZff36dTQ1NTmkLiFEm983arUaarVa9j6j0QgA8PHxAQBUVlbCYDAgLi7Oop4xY8agpKQEL730Eo4fP47m5maLGL1ej5CQEJSUlCA+Ph6lpaXQaDQYNWqUFBMVFQWNRoOSkhIEBwd36nu5TGLw/fffIyAgoLubQUREdqiursagQYOcUvf169cRFPgzGGpbOg7uhJ/97GdoaGiwOPfOO+8gIyOj3XuEEEhLS8Ojjz6KkJAQAIDBYAAA+Pv7W8T6+/vjwoULUoy7uzv69+/fJqb1foPBAD8/vzaf6efnJ8V0hsskBl5eXgCAR/E0+sCtm1tDRES2uIFmFGOf9P/LnaGpqQmG2hZUHg+Et5d9vRKmq2YERVxAdXU1vL29pfMd9Ra88sor+Pvf/47i4uI2127vfbDWI3G722OsxXemnlu5TGLQ+qX7wA19VEwMiIjuKuLmf7piKNjbq5fdiYFUl7e3RWIgZ8GCBdizZw+KioosekW0Wi2Am//i1+l00vna2lqpF0Gr1aKpqQl1dXUWvQa1tbWIiYmRYi5fvtzmc3/44Yc2vRFyuCqBiIgUpUWYHVI6SwiBV155BTt37sTBgwcRFBRkcT0oKAharRZ5eXnSuaamJhQWFkq/9CMiIuDm5mYRU1NTg9OnT0sx0dHRMBqN+Oqrr6SYI0eOwGg0SjGd4TI9BkRERJ1hhoC5tYvCjjo6a/78+fj000+xe/dueHl5SeP9Go0Gffv2hUqlQmpqKt59910MHToUQ4cOxbvvvot+/fph5syZUuycOXOwcOFCDBgwAD4+PkhPT0doaCiefPJJAMCDDz6I8ePHIyUlBWvXrgUAzJ07F4mJiZ2eeAg4scdg9erVCAoKgoeHByIiIvDFF1/IxhcWFiIiIgIeHh64//778fHHHzuraUREpGBmB/1fZ61ZswZGoxGxsbHQ6XRS2b59uxSzaNEipKam4uWXX0ZkZCQuXbqEAwcOWMy5WL58OSZPnoxp06Zh9OjR6NevHz7//HP07t1bitm6dStCQ0MRFxeHuLg4PPzww9i8ebNNz0clhLAvbbJi+/btSEpKwurVqzF69GisXbsWf/rTn3D27FkMHjy4TXxlZSVCQkKQkpKCl156CV9++SVefvllfPbZZ5g6dWqnPtNkMkGj0SAWkzjHgIjoLnNDNKMAu2E0Gjs9Zm+r1t8T35cPcsjkQ33wRae2t7s4JTEYNWoURowYgTVr1kjnHnzwQUyePBmZmZlt4l9//XXs2bMH586dk87NmzcPX3/9NUpLSzv1mUwMiIjuXl2ZGFR/c69DEoOABy65ZGLg8KGEpqYmHD9+3GITBgCIi4tDSUmJ1XtKS0vbxMfHx+PYsWNobm62ek9jYyNMJpNFISIi6kjrHAN7i6tyeGLw448/oqWlxepGDe1tsGAwGKzG37hxAz/++KPVezIzM6WXRGg0Gm5uRERE5ABOm3xo60YN1uKtnW+1ZMkSGI1GqVRXV9vZYiIiUgIzBFrsLK7cY+Dw5Yq+vr7o3bt3m96BWzdquJ1Wq7Ua36dPHwwYMMDqPZ3Zj5qIiOh2Xb1c8W7j8B4Dd3d3REREWGzCAAB5eXntbrAQHR3dJv7AgQOIjIyEmxsnEhIREXUVpwwlpKWl4U9/+hPWr1+Pc+fO4bXXXkNVVRXmzZsH4OYwwKxZs6T4efPm4cKFC0hLS8O5c+ewfv16ZGVlIT093RnNIyIiBWsRwiHFVTll58Pp06fjypUr+O///m/U1NQgJCQE+/btQ2BgIICb2zhWVVVJ8UFBQdi3bx9ee+01fPTRR9Dr9fjwww87vYcBERFRZ5n/Xeytw1U5ZR+D7sB9DIiI7l5duY/BN+f84WXnPgZXr5rxwIOXXXIfA74rgYiIFKV1ZYG9dbgqJgZERKQoLeJmsbcOV8XEgIiIFIVzDOQ5bYMjIiIiuvuwx4CIiBTFDBVa0P5OvJ2tw1UxMSAiIkUxi5vF3jpcFYcSiIiISMIeAyIiUpQWBwwl2Ht/T8bEgIiIFIWJgTwOJRAREZGEPQZERKQoZqGCWdi5KsHO+3syJgZERKQoHEqQx6EEIiIikrDHgIiIFKUFvdBi57+LWxzUlp6IiQERESmKcMAcA8E5BkRERK6BcwzkcY4BERERSdhjQEREitIieqFF2DnHwIXflcDEgIiIFMUMFcx2dpib4bqZAYcSiIiISMIeAyIiUhROPpTHxICIiBTFMXMMOJRARERECsAeAyIiUpSbkw/tfIkShxKIiIhcg9kBWyJzVQIREREpAhMDIiJSlNbJh/YWWxQVFWHChAnQ6/VQqVTYtWuXxXWVSmW1/OEPf5BiYmNj21yfMWOGRT11dXVISkqCRqOBRqNBUlIS6uvrbWorEwMiIlIUM3o5pNji2rVrCAsLw6pVq6xer6mpsSjr16+HSqXC1KlTLeJSUlIs4tauXWtxfebMmSgrK0Nubi5yc3NRVlaGpKQkm9rKOQZERKQoLUKFFjvfjmjr/QkJCUhISGj3ulartTjevXs3xo4di/vvv9/ifL9+/drEtjp37hxyc3Nx+PBhjBo1CgCwbt06REdHo7y8HMHBwZ1qK3sMiIiI7pDJZLIojY2Ndtd5+fJl7N27F3PmzGlzbevWrfD19cXw4cORnp6Oq1evStdKS0uh0WikpAAAoqKioNFoUFJS0unPZ48BEREpSosDViW0/HtVQkBAgMX5d955BxkZGXbVvXHjRnh5eWHKlCkW559//nkEBQVBq9Xi9OnTWLJkCb7++mvk5eUBAAwGA/z8/NrU5+fnB4PB0OnPZ2JARESKYha9YLZz50Pzv3c+rK6uhre3t3RerVbbVS8ArF+/Hs8//zw8PDwszqekpEg/h4SEYOjQoYiMjMSJEycwYsQIADcnMd5OCGH1fHs4lEBERHSHvL29LYq9icEXX3yB8vJy/OpXv+owdsSIEXBzc0NFRQWAm/MULl++3Cbuhx9+gL+/f6fbwMSAiIgUpXUowd7iDFlZWYiIiEBYWFiHsWfOnEFzczN0Oh0AIDo6GkajEV999ZUUc+TIERiNRsTExHS6DRxKICIiRTHD9lUF1uqwRUNDA86fPy8dV1ZWoqysDD4+Phg8eDCAmxMZ//znP2PZsmVt7v/nP/+JrVu34umnn4avry/Onj2LhQsXIjw8HKNHjwYAPPjggxg/fjxSUlKkZYxz585FYmJip1ckAE7oMcjMzMTIkSPh5eUFPz8/TJ48GeXl5bL3FBQUWN3Y4ZtvvnF084iIiLrcsWPHEB4ejvDwcABAWloawsPD8fbbb0sx27ZtgxACzz33XJv73d3d8be//Q3x8fEIDg7Gq6++iri4OOTn56N3795S3NatWxEaGoq4uDjExcXh4YcfxubNm21qq0oIx747cvz48ZgxYwZGjhyJGzdu4M0338SpU6dw9uxZeHp6Wr2noKAAY8eORXl5ucUkjoEDB1p8YTkmkwkajQaxmIQ+KjeHfBciIuoaN0QzCrAbRqPR4veAI7X+nlhzYiT6/sy+DvN/NdzAf4046tT2dheHDyXk5uZaHG/YsAF+fn44fvw4Hn/8cdl7/fz8cM899zi6SURERJI72dLYWh2uyunfzGg0AgB8fHw6jA0PD4dOp8O4ceNw6NAh2djGxsY2G0sQERGRfZyaGAghkJaWhkcffRQhISHtxul0OnzyySfYsWMHdu7cieDgYIwbNw5FRUXt3pOZmSm9JEKj0bTZZIKIiMgaM1QOKa7K4XMMbjV//nzs3bsXxcXFGDRokE33TpgwASqVCnv27LF6vbGx0WLrSZPJhICAAM4xICK6C3XlHIPlx2IcMsfgtcgSzjGwxYIFC7Bnzx4UFRXZnBQAN/d33rJlS7vX1Wq1Q3aYIiIiZXHMlsiuO8fA4YmBEAILFixATk4OCgoKEBQUdEf1nDx5Utq0gYiIiLqGwxOD+fPn49NPP8Xu3bvh5eUlvbhBo9Ggb9++AIAlS5bg0qVL2LRpEwBgxYoVuO+++zB8+HA0NTVhy5Yt2LFjB3bs2OHo5hERkcKZhQpmezc4svP+nszhicGaNWsAALGxsRbnN2zYgNmzZwMAampqUFVVJV1rampCeno6Ll26hL59+2L48OHYu3cvnn76aUc3j4iIFM7sgKEEM4cSOq8zcxmzs7MtjhctWoRFixY5uilERERkI74rgYiIFMUxr11mjwEREZFLaIEKLXbuQ2Dv/T2Z66Y8REREZDP2GBARkaJwKEEeEwMiIlKUFtg/FNDimKb0SK6b8hAREZHN2GNARESKwqEEeUwMiIhIUVpEL7TY+Yvd3vt7MiYGRESkKMIBr00WXK5IRERESsAeAyIiUhQOJchjYkBERIrCtyvKc92Uh4iIiGzGHgMiIlKUFge8dtne+3syJgZERKQoHEqQ57opDxEREdmMPQZERKQoZvSC2c5/F9t7f0/GxICIiBSlRajQYudQgL3392Sum/IQERGRzdhjQEREisLJh/KYGBARkaIIB7xdUXDnQyIiItfQAhVa7HwJkr3392Sum/IQERGRzdhjQEREimIW9s8RMAsHNaYHYmJARESKYnbAHAN77+/JXPebERERkc2YGBARkaKYoXJIsUVRUREmTJgAvV4PlUqFXbt2WVyfPXs2VCqVRYmKirKIaWxsxIIFC+Dr6wtPT09MnDgRFy9etIipq6tDUlISNBoNNBoNkpKSUF9fb1NbmRgQEZGitO58aG+xxbVr1xAWFoZVq1a1GzN+/HjU1NRIZd++fRbXU1NTkZOTg23btqG4uBgNDQ1ITExES0uLFDNz5kyUlZUhNzcXubm5KCsrQ1JSkk1t5RwDIiIiJ0tISEBCQoJsjFqthlartXrNaDQiKysLmzdvxpNPPgkA2LJlCwICApCfn4/4+HicO3cOubm5OHz4MEaNGgUAWLduHaKjo1FeXo7g4OBOtZU9BkREpCitkw/tLQBgMpksSmNj4x23q6CgAH5+fhg2bBhSUlJQW1srXTt+/Diam5sRFxcnndPr9QgJCUFJSQkAoLS0FBqNRkoKACAqKgoajUaK6QwmBkREpChmqKRtke+4/HuOQUBAgDSer9FokJmZeUdtSkhIwNatW3Hw4EEsW7YMR48exRNPPCElGgaDAe7u7ujfv7/Fff7+/jAYDFKMn59fm7r9/PykmM7gUAIREdEdqq6uhre3t3SsVqvvqJ7p06dLP4eEhCAyMhKBgYHYu3cvpkyZ0u59QgioVP+Z73Drz+3FdIQ9BkREpCjCASsSxL97DLy9vS3KnSYGt9PpdAgMDERFRQUAQKvVoqmpCXV1dRZxtbW18Pf3l2IuX77cpq4ffvhBiukMJgZERKQodg8jOODtjB25cuUKqqurodPpAAARERFwc3NDXl6eFFNTU4PTp08jJiYGABAdHQ2j0YivvvpKijly5AiMRqMU0xkcSiAiIkXpjp0PGxoacP78eem4srISZWVl8PHxgY+PDzIyMjB16lTodDp89913eOONN+Dr64tnnnkGAKDRaDBnzhwsXLgQAwYMgI+PD9LT0xEaGiqtUnjwwQcxfvx4pKSkYO3atQCAuXPnIjExsdMrEgAmBkRERE537NgxjB07VjpOS0sDACQnJ2PNmjU4deoUNm3ahPr6euh0OowdOxbbt2+Hl5eXdM/y5cvRp08fTJs2Df/6178wbtw4ZGdno3fv3lLM1q1b8eqrr0qrFyZOnCi7d4I1KiGEQ18FkZGRgd/+9rcW526dNWlNYWEh0tLScObMGej1eixatAjz5s2z6XNNJhM0Gg1iMQl9VG531HYiIuoeN0QzCrAbRqPRYjKfI7X+nph04EW4ebrbVVfztSbsjlvv1PZ2F6f0GAwfPhz5+fnS8a3ZzO0qKyvx9NNPIyUlBVu2bMGXX36Jl19+GQMHDsTUqVOd0TwiIlKwO9nS2FodrsopiUGfPn3a3b3pdh9//DEGDx6MFStWALg5RnLs2DH88Y9/lE0MGhsbLTaSMJlMdrWZiIiInLQqoaKiAnq9HkFBQZgxYwa+/fbbdmNLS0stdnICgPj4eBw7dgzNzc3t3peZmWmxqURAQIDD2k9ERK7rbliV0J0cnhiMGjUKmzZtwv79+7Fu3ToYDAbExMTgypUrVuMNBkOb9ZX+/v64ceMGfvzxx3Y/Z8mSJTAajVKprq526PcgIiLXxMRAnsOHEm59SURoaCiio6MxZMgQbNy4UZqFebvbd2RqnQ8pt1OTWq122EYSREREdJPTlyt6enoiNDRU2r3pdlqtts2KhdraWvTp0wcDBgxwdvOIiEhhHPEvflfuMXD6zoeNjY04d+6ctHvT7aKjoy12cgKAAwcOIDIyEm5uXHZIRESOxaEEeQ5PDNLT01FYWIjKykocOXIEzz77LEwmE5KTkwHcnBswa9YsKX7evHm4cOEC0tLScO7cOaxfvx5ZWVlIT093dNOIiIioAw4fSrh48SKee+45/Pjjjxg4cCCioqJw+PBhBAYGAri5t3NVVZUUHxQUhH379uG1117DRx99BL1ejw8//JB7GBARkVMI2L8PgUN3BuxhHJ4YbNu2TfZ6dnZ2m3NjxozBiRMnHN0UIiKiNjjHQB7flUBERIrCxEAeX7tMREREEvYYEBGRorDHQB4TAyIiUhQmBvI4lEBEREQS9hgQEZGiCKGCsPNf/Pbe35MxMSAiIkUxQ2X3Pgb23t+TcSiBiIiIJOwxICIiReHkQ3lMDIiISFE4x0AehxKIiIhIwh4DIiJSFA4lyGNiQEREisKhBHlMDIiISFGEA3oMXDkx4BwDIiIikrDHgIiIFEUAEML+OlwVEwMiIlIUM1RQcefDdnEogYiIiCTsMSAiIkXhqgR5TAyIiEhRzEIFFfcxaBeHEoiIiEjCHgMiIlIUIRywKsGFlyUwMSAiIkXhHAN5HEogIiIiCXsMiIhIUdhjII+JARERKQpXJcjjUAIRESlK6+RDe4stioqKMGHCBOj1eqhUKuzatUu61tzcjNdffx2hoaHw9PSEXq/HrFmz8P3331vUERsbC5VKZVFmzJhhEVNXV4ekpCRoNBpoNBokJSWhvr7eprYyMSAiInKya9euISwsDKtWrWpz7aeffsKJEyfw1ltv4cSJE9i5cyf+8Y9/YOLEiW1iU1JSUFNTI5W1a9daXJ85cybKysqQm5uL3NxclJWVISkpyaa2ciiBiIgU5ea/+O2dY3DzvyaTyeK8Wq2GWq1uE5+QkICEhASrdWk0GuTl5VmcW7lyJR555BFUVVVh8ODB0vl+/fpBq9VarefcuXPIzc3F4cOHMWrUKADAunXrEB0djfLycgQHB3fqu7HHgIiIFKV18qG9BQACAgKkbnuNRoPMzEyHtNFoNEKlUuGee+6xOL9161b4+vpi+PDhSE9Px9WrV6VrpaWl0Gg0UlIAAFFRUdBoNCgpKen0Z7PHgIiI6A5VV1fD29tbOrbWW2Cr69evY/HixZg5c6ZF3c8//zyCgoKg1Wpx+vRpLFmyBF9//bXU22AwGODn59emPj8/PxgMhk5/PhMDIiJSFPHvYm8dAODt7W3xy9tezc3NmDFjBsxmM1avXm1xLSUlRfo5JCQEQ4cORWRkJE6cOIERI0YAAFSqtkMkQgir59vDoQQiIlIURw4lOFJzczOmTZuGyspK5OXldZhwjBgxAm5ubqioqAAAaLVaXL58uU3cDz/8AH9//063g4kBERFRN2tNCioqKpCfn48BAwZ0eM+ZM2fQ3NwMnU4HAIiOjobRaMRXX30lxRw5cgRGoxExMTGdbguHEoiISFkcOZbQSQ0NDTh//rx0XFlZibKyMvj4+ECv1+PZZ5/FiRMn8Je//AUtLS3SnAAfHx+4u7vjn//8J7Zu3Yqnn34avr6+OHv2LBYuXIjw8HCMHj0aAPDggw9i/PjxSElJkZYxzp07F4mJiZ1ekQA4ocfgvvvua7MBg0qlwvz5863GFxQUWI3/5ptvHN00IiIiwBHDCDYOJRw7dgzh4eEIDw8HAKSlpSE8PBxvv/02Ll68iD179uDixYv4+c9/Dp1OJ5XW1QTu7u7429/+hvj4eAQHB+PVV19FXFwc8vPz0bt3b+lztm7ditDQUMTFxSEuLg4PP/wwNm/ebFNbHd5jcPToUbS0tEjHp0+fxlNPPYVf/vKXsveVl5dbjKcMHDjQ0U0jIiLqltcux8bGQsjcJHcNuLkssrCwsMPP8fHxwZYtW2xr3G0cnhjc/gt96dKlGDJkCMaMGSN7n5+fX5v1mkRERNS1nDr5sKmpCVu2bMGLL77Y4VKJ8PBw6HQ6jBs3DocOHeqw7sbGRphMJotCRETUkZ66KqGncGpisGvXLtTX12P27Nntxuh0OnzyySfYsWMHdu7cieDgYIwbNw5FRUWydWdmZlrsNhUQEODg1hMRkUtqnSNgb3FRKtHRwIYd4uPj4e7ujs8//9ym+yZMmACVSoU9e/a0G9PY2IjGxkbp2GQyISAgALGYhD4qtztuMxERdb0bohkF2A2j0ejQDYNuZTKZoNFocF/WW+jVz8Ouusw/Xcd3c37n1PZ2F6ctV7xw4QLy8/Oxc+dOm++NiorqcPJEey+qICIiktMdkw/vJk5LDDZs2AA/Pz/84he/sPnekydPShs2EBEROVQ37GNwN3FKYmA2m7FhwwYkJyejTx/Lj1iyZAkuXbqETZs2AQBWrFiB++67D8OHD5cmK+7YsQM7duxwRtOIiIhIhlMSg/z8fFRVVeHFF19sc62mpgZVVVXScVNTE9LT03Hp0iX07dsXw4cPx969e/H00087o2lERKRwjlhV4MqrEpySGMTFxbW7WUN2drbF8aJFi7Bo0SJnNIOIiMg6Fx4KsBdfokREREQSvkSJiIgUhUMJ8pgYEBGRsnBVgiwmBkREpDCqfxd763BNnGNAREREEvYYEBGRsnAoQRYTAyIiUhYmBrI4lEBEREQS9hgQEZGyOOK1yVyuSERE5Br4dkV5HEogIiIiCXsMiIhIWTj5UBYTAyIiUhbOMZDFoQQiIiKSsMeAiIgURSVuFnvrcFVMDIiISFk4x0AWEwMiIlIWzjGQxTkGREREJGGPARERKQuHEmQxMSAiImVhYiCLQwlEREQkYY8BEREpC3sMZDExICIiZeGqBFkcSiAiIiIJewyIiEhRuPOhPCYGRESkLJxjIItDCURERE5WVFSECRMmQK/XQ6VSYdeuXRbXhRDIyMiAXq9H3759ERsbizNnzljENDY2YsGCBfD19YWnpycmTpyIixcvWsTU1dUhKSkJGo0GGo0GSUlJqK+vt6mtTAyIiIic7Nq1awgLC8OqVausXn///ffxwQcfYNWqVTh69Ci0Wi2eeuopXL16VYpJTU1FTk4Otm3bhuLiYjQ0NCAxMREtLS1SzMyZM1FWVobc3Fzk5uairKwMSUlJNrWVQwlERKQoKjhgjsG//2symSzOq9VqqNXqNvEJCQlISEiwWpcQAitWrMCbb76JKVOmAAA2btwIf39/fPrpp3jppZdgNBqRlZWFzZs348knnwQAbNmyBQEBAcjPz0d8fDzOnTuH3NxcHD58GKNGjQIArFu3DtHR0SgvL0dwcHCnvht7DIiISFlalyvaWwAEBARI3fYajQaZmZk2N6eyshIGgwFxcXHSObVajTFjxqCkpAQAcPz4cTQ3N1vE6PV6hISESDGlpaXQaDRSUgAAUVFR0Gg0UkxnsMeAiIjoDlVXV8Pb21s6ttZb0BGDwQAA8Pf3tzjv7++PCxcuSDHu7u7o379/m5jW+w0GA/z8/NrU7+fnJ8V0BhMDIiJSFgeuSvD29rZIDOyhUllumiSEaHOuTTNui7EW35l6bsWhBCIiUhbhoOIgWq0WANr8q762tlbqRdBqtWhqakJdXZ1szOXLl9vU/8MPP7TpjZDDxICIiKgbBQUFQavVIi8vTzrX1NSEwsJCxMTEAAAiIiLg5uZmEVNTU4PTp09LMdHR0TAajfjqq6+kmCNHjsBoNEoxncGhBCIiUpTu2PmwoaEB58+fl44rKytRVlYGHx8fDB48GKmpqXj33XcxdOhQDB06FO+++y769euHmTNnAgA0Gg3mzJmDhQsXYsCAAfDx8UF6ejpCQ0OlVQoPPvggxo8fj5SUFKxduxYAMHfuXCQmJnZ6RQJwBz0GjtikwZodO3bgoYceglqtxkMPPYScnBxbm0ZERNSxbhhKOHbsGMLDwxEeHg4ASEtLQ3h4ON5++20AwKJFi5CamoqXX34ZkZGRuHTpEg4cOAAvLy+pjuXLl2Py5MmYNm0aRo8ejX79+uHzzz9H7969pZitW7ciNDQUcXFxiIuLw8MPP4zNmzfb1FaVEMKmr/fXv/4VX375JUaMGIGpU6ciJycHkydPlq6/9957+J//+R9kZ2dj2LBh+P3vf4+ioiKUl5dbfMFblZaW4rHHHsPvfvc7PPPMM8jJycHbb7+N4uJii2UXckwmEzQaDWIxCX1UbrZ8JSIi6mY3RDMKsBtGo9Fhk/lu1/p74r7f/w96eXjYVZf5+nV895s3ndre7mJzYmBxs0plkRgIIaDX65GamorXX38dwM0tHP39/fHee+/hpZdeslrP9OnTYTKZ8Ne//lU6N378ePTv3x+fffZZp9rCxICI6O7VpYnB7xyUGLzlmomBQycfdmaTBmtKS0st7gGA+Ph42XsaGxthMpksChERUUda5xjYW1yVQxMDuU0a5DZXMBgMNt+TmZlpsdtUQECAHS0nIiIiwEnLFe9kkwZb71myZAmMRqNUqqur77zBRESkHA7cEtkVOXS54q2bNOh0Oun8rRswtHef3MYO1rT3ogoiIiJZDtz50BU5tMegM5s0WBMdHW1xDwAcOHDApg0ZiIiIOoNzDOTZ3GNg7yYNADBr1izce++90luofv3rX+Pxxx/He++9h0mTJmH37t3Iz89HcXGxA74iERERdZbNicGxY8cwduxY6TgtLQ0AkJycjOzsbCxatAj/+te/8PLLL6Ourg6jRo1qs0lDVVUVevX6T2dFTEwMtm3bht/85jd46623MGTIEGzfvr3TexgQERF1GocSZNm1j0FPwn0MiIjuXl25j8H9b72L3nbuY9By/Tq+/d0b3MeAiIiIXBtfokRERMrCoQRZTAyIiEhZmBjI4lACERERSdhjQEREiuKIfQhceR8D9hgQERGRhIkBERERSTiUQEREysLJh7KYGBARkaJwjoE8JgZERKQ8LvyL3V6cY0BEREQS9hgQEZGycI6BLCYGRESkKJxjII9DCURERCRhjwERESkLhxJkMTEgIiJF4VCCPA4lEBERkYQ9BkREpCwcSpDFxICIiJSFiYEsDiUQERGRhD0GRESkKJx8KI+JARERKQuHEmQxMSAiImVhYiCLcwyIiIhIwh4DIiJSFM4xkMfEgIiIlIVDCbI4lEBERORE9913H1QqVZsyf/58AMDs2bPbXIuKirKoo7GxEQsWLICvry88PT0xceJEXLx40SntZWJARESK0jqUYG/prKNHj6KmpkYqeXl5AIBf/vKXUsz48eMtYvbt22dRR2pqKnJycrBt2zYUFxejoaEBiYmJaGlpccgzuRWHEoiISFm6eChh4MCBFsdLly7FkCFDMGbMGOmcWq2GVqu1er/RaERWVhY2b96MJ598EgCwZcsWBAQEID8/H/Hx8ba3XwZ7DIiIiO6QyWSyKI2NjbLxTU1N2LJlC1588UWoVCrpfEFBAfz8/DBs2DCkpKSgtrZWunb8+HE0NzcjLi5OOqfX6xESEoKSkhKHfycmBkREpCzCQQVAQEAANBqNVDIzM2U/eteuXaivr8fs2bOlcwkJCdi6dSsOHjyIZcuW4ejRo3jiiSekJMNgMMDd3R39+/e3qMvf3x8Gg8GeJ2EVhxKIiEhRVP8u9tYBANXV1fD29pbOq9Vq2fuysrKQkJAAvV4vnZs+fbr0c0hICCIjIxEYGIi9e/diypQp7dYlhLDodXAUJgZERER3yNvb2yIxkHPhwgXk5+dj586dsnE6nQ6BgYGoqKgAAGi1WjQ1NaGurs6i16C2thYxMTF33vh2cCiBiIiUxYFDCbbYsGED/Pz88Itf/EI27sqVK6iuroZOpwMAREREwM3NTVrNAAA1NTU4ffq0UxID9hgQEZGidMfOh2azGRs2bEBycjL69PnPr96GhgZkZGRg6tSp0Ol0+O677/DGG2/A19cXzzzzDABAo9Fgzpw5WLhwIQYMGAAfHx+kp6cjNDRUWqXgSDb3GBQVFWHChAnQ6/VQqVTYtWuXdK25uRmvv/46QkND4enpCb1ej1mzZuH777+XrTM7O9vq5g/Xr1+3+QsRERHJ6oYeg/z8fFRVVeHFF1+0ON+7d2+cOnUKkyZNwrBhw5CcnIxhw4ahtLQUXl5eUtzy5csxefJkTJs2DaNHj0a/fv3w+eefo3fv3nfwAOTZ3GNw7do1hIWF4YUXXsDUqVMtrv300084ceIE3nrrLYSFhaGurg6pqamYOHEijh07Jluvt7c3ysvLLc55eHjY2jwiIqIeJy4uDkK0zSb69u2L/fv3d3i/h4cHVq5ciZUrVzqjeRZsTgwSEhKQkJBg9ZpGo7EYAwGAlStX4pFHHkFVVRUGDx7cbr0qlardzR2IiIgcyoXfdWAvp08+NBqNUKlUuOeee2TjGhoaEBgYiEGDBiExMREnT56UjW9sbGyzsQQREVFHunpL5LuNUxOD69evY/HixZg5c6bsco4HHngA2dnZ2LNnDz777DN4eHhg9OjR0lINazIzMy02lQgICHDGVyAiIlIUpyUGzc3NmDFjBsxmM1avXi0bGxUVhf/3//4fwsLC8Nhjj+F///d/MWzYMNmxlCVLlsBoNEqlurra0V+BiIhcUTctV7xbOGW5YnNzM6ZNm4bKykocPHiw05s/tOrVqxdGjhwp22OgVqs73GGKiIjodt2xXPFu4vAeg9akoKKiAvn5+RgwYIDNdQghUFZWJm3uQERERF3D5h6DhoYGnD9/XjqurKxEWVkZfHx8oNfr8eyzz+LEiRP4y1/+gpaWFukFDz4+PnB3dwcAzJo1C/fee6/0sonf/va3iIqKwtChQ2EymfDhhx+irKwMH330kSO+IxER0X908WuX7zY2JwbHjh3D2LFjpeO0tDQAQHJyMjIyMrBnzx4AwM9//nOL+w4dOoTY2FgAQFVVFXr1+k9nRX19PebOnQuDwQCNRoPw8HAUFRXhkUcesbV5REREsjiUIM/mxCA2NtbqJg2t5K61KigosDhevnw5li9fbmtTiIiIyMH4rgQiIlIWDiXIYmJARETKwsRAFhMDIiJSFM4xkOf0LZGJiIjo7sEeAyIiUhYOJchiYkBERIqiEgKqTqyg66gOV8WhBCIiIpKwx4CIiJSFQwmymBgQEZGicFWCPA4lEBERkYQ9BkREpCwcSpDFxICIiBSFQwnyOJRAREREEvYYEBGRsnAoQRYTAyIiUhQOJchjYkBERMrCHgNZnGNAREREEvYYEBGR4rjyUIC9mBgQEZGyCHGz2FuHi+JQAhEREUnYY0BERIrCVQnymBgQEZGycFWCLA4lEBERkYQ9BkREpCgq881ibx2uiokBEREpC4cSZHEogYiIiCRMDIiISFFaVyXYWzorIyMDKpXKomi1Wum6EAIZGRnQ6/Xo27cvYmNjcebMGYs6GhsbsWDBAvj6+sLT0xMTJ07ExYsXHfVILDAxICIiZWnd4MjeYoPhw4ejpqZGKqdOnZKuvf/++/jggw+watUqHD16FFqtFk899RSuXr0qxaSmpiInJwfbtm1DcXExGhoakJiYiJaWFoc9llacY0BERIrSHfsY9OnTx6KXoJUQAitWrMCbb76JKVOmAAA2btwIf39/fPrpp3jppZdgNBqRlZWFzZs348knnwQAbNmyBQEBAcjPz0d8fLx9X+Y27DEgIiK6QyaTyaI0NjZajauoqIBer0dQUBBmzJiBb7/9FgBQWVkJg8GAuLg4KVatVmPMmDEoKSkBABw/fhzNzc0WMXq9HiEhIVKMIzExICIiZREOKgACAgKg0WikkpmZ2ebjRo0ahU2bNmH//v1Yt24dDAYDYmJicOXKFRgMBgCAv7+/xT3+/v7SNYPBAHd3d/Tv37/dGEfiUAIRESmKI4cSqqur4e3tLZ1Xq9VtYhMSEqSfQ0NDER0djSFDhmDjxo2Iioq6WZ9KZXGPEKLNudt1JuZOsMeAiIjoDnl7e1sUa4nB7Tw9PREaGoqKigpp3sHt//Kvra2VehG0Wi2amppQV1fXbowjMTEgIiJl6YZVCbdqbGzEuXPnoNPpEBQUBK1Wi7y8POl6U1MTCgsLERMTAwCIiIiAm5ubRUxNTQ1Onz4txTgShxKIiEhRunpVQnp6OiZMmIDBgwejtrYWv//972EymZCcnAyVSoXU1FS8++67GDp0KIYOHYp3330X/fr1w8yZMwEAGo0Gc+bMwcKFCzFgwAD4+PggPT0doaGh0ioFR7K5x6CoqAgTJkyAXq+HSqXCrl27LK7Pnj27zUYOrWMocnbs2IGHHnoIarUaDz30EHJycmxtGhERUY9z8eJFPPfccwgODsaUKVPg7u6Ow4cPIzAwEACwaNEipKam4uWXX0ZkZCQuXbqEAwcOwMvLS6pj+fLlmDx5MqZNm4bRo0ejX79++Pzzz9G7d2+Ht9fmHoNr164hLCwML7zwAqZOnWo1Zvz48diwYYN07O7uLltnaWkppk+fjt/97nd45plnkJOTg2nTpqG4uBijRo2ytYlERETt6+J3JWzbtk32ukqlQkZGBjIyMtqN8fDwwMqVK7Fy5crOf/AdsjkxSEhIsJhhaY1arba6kUN7VqxYgaeeegpLliwBACxZsgSFhYVYsWIFPvvsM1ubSERE1K7u2ODobuKUyYcFBQXw8/PDsGHDkJKSgtraWtn40tJSi40bACA+Pl5244bGxsY2G0sQERGRfRyeGCQkJGDr1q04ePAgli1bhqNHj+KJJ55odzco4OYyDbnNHazJzMy02FQiICDAYd+BiIhcmFk4prgoh69KmD59uvRzSEgIIiMjERgYiL1790r7QFtj6+YOS5YsQVpamnRsMpmYHBARUce6eI7B3cbpyxV1Oh0CAwNRUVHRboxWq5Xd3MEatVrdqY0kiIiIbqWCA+YYOKQlPZPTNzi6cuUKqqurodPp2o2Jjo622LgBAA4cOOCUjRuIiIiofTb3GDQ0NOD8+fPScWVlJcrKyuDj4wMfHx9kZGRg6tSp0Ol0+O677/DGG2/A19cXzzzzjHTPrFmzcO+990ovm/j1r3+Nxx9/HO+99x4mTZqE3bt3Iz8/H8XFxQ74ikRERLewc+dCqQ4XZXNicOzYMYwdO1Y6bh3nT05Oxpo1a3Dq1Cls2rQJ9fX10Ol0GDt2LLZv326xUUNVVRV69fpPZ0VMTAy2bduG3/zmN3jrrbcwZMgQbN++nXsYEBGRw3G5ojybE4PY2FgImUxp//79HdZRUFDQ5tyzzz6LZ5991tbmEBERkQPxXQlERKQsXJUgi4kBEREpikoIqOycI2Dv/T0ZX7tMREREEvYYEBGRspj/Xeytw0UxMSAiIkXhUII8DiUQERGRhD0GRESkLFyVIIuJARERKQt3PpTFxICIiBSFOx/K4xwDIiIikrDHgIiIlIVDCbKYGBARkaKozDeLvXW4Kg4lEBERkYQ9BkREpCwcSpDFxICIiJSF+xjI4lACERERSdhjQEREisJ3JchjYkBERMrCOQayOJRAREREEvYYEBGRsggA9u5D4LodBkwMiIhIWTjHQB4TAyIiUhYBB8wxcEhLeiTOMSAiIiIJewyIiEhZuCpBFhMDIiJSFjMAlQPqcFEcSiAiIiIJEwMiIlKU1lUJ9pbOyszMxMiRI+Hl5QU/Pz9MnjwZ5eXlFjGzZ8+GSqWyKFFRURYxjY2NWLBgAXx9feHp6YmJEyfi4sWLDnkmt2JiQEREytI6x8De0kmFhYWYP38+Dh8+jLy8PNy4cQNxcXG4du2aRdz48eNRU1MjlX379llcT01NRU5ODrZt24bi4mI0NDQgMTERLS0tDnksrTjHgIiIyIlyc3Mtjjds2AA/Pz8cP34cjz/+uHRerVZDq9VarcNoNCIrKwubN2/Gk08+CQDYsmULAgICkJ+fj/j4eIe1lz0GRESkLA7sMTCZTBalsbGxw483Go0AAB8fH4vzBQUF8PPzw7Bhw5CSkoLa2lrp2vHjx9Hc3Iy4uDjpnF6vR0hICEpKShzxVCRMDIiISFkcmBgEBARAo9FIJTMzs4OPFkhLS8Ojjz6KkJAQ6XxCQgK2bt2KgwcPYtmyZTh69CieeOIJKdEwGAxwd3dH//79Lerz9/eHwWBw6OPhUAIREdEdqq6uhre3t3SsVqtl41955RX8/e9/R3FxscX56dOnSz+HhIQgMjISgYGB2Lt3L6ZMmdJufUIIqFT2rr20xB4DIiJSFrODCgBvb2+LIpcYLFiwAHv27MGhQ4cwaNAg2SbqdDoEBgaioqICAKDVatHU1IS6ujqLuNraWvj7+9v09TvCxICIiBSlq5crCiHwyiuvYOfOnTh48CCCgoI6vOfKlSuorq6GTqcDAERERMDNzQ15eXlSTE1NDU6fPo2YmBjbH4IMDiUQEZGydPGWyPPnz8enn36K3bt3w8vLS5oToNFo0LdvXzQ0NCAjIwNTp06FTqfDd999hzfeeAO+vr545plnpNg5c+Zg4cKFGDBgAHx8fJCeno7Q0FBplYKjMDEgIiJyojVr1gAAYmNjLc5v2LABs2fPRu/evXHq1Cls2rQJ9fX10Ol0GDt2LLZv3w4vLy8pfvny5ejTpw+mTZuGf/3rXxg3bhyys7PRu3dvh7bX5qGEoqIiTJgwAXq9HiqVCrt27bK4fvvOTa3lD3/4Q7t1ZmdnW73n+vXrNn8hIiIiWWbhmNJJQgirZfbs2QCAvn37Yv/+/aitrUVTUxMuXLiA7OxsBAQEWNTj4eGBlStX4sqVK/jpp5/w+eeft4lxBJt7DK5du4awsDC88MILmDp1apvrNTU1Fsd//etfMWfOHKuxt/L29m6zRaSHh4etzSMiIpLHtyvKsjkxSEhIQEJCQrvXb9+1affu3Rg7dizuv/9+2XpVKlW7Oz5Z09jYaLGRhMlk6vS9REREZJ1TVyVcvnwZe/fuxZw5czqMbWhoQGBgIAYNGoTExEScPHlSNj4zM9NiUwlndKcQEZErcsTmRq7bY+DUxGDjxo3w8vKS3ZwBAB544AFkZ2djz549+Oyzz+Dh4YHRo0dL6zetWbJkCYxGo1Sqq6sd3XwiInJFXfwSpbuNU1clrF+/Hs8//3yHcwWioqIsXi85evRojBgxAitXrsSHH35o9R61Wt3hDlNERERkG6clBl988QXKy8uxfft2m+/t1asXRo4cKdtjQEREdEfMDhgKsGFVwt3GaUMJWVlZiIiIQFhYmM33CiFQVlYm7fhERETkMMLsmOKibO4xaGhowPnz56XjyspKlJWVwcfHB4MHDwZwc4XAn//8ZyxbtsxqHbNmzcK9994rvYXqt7/9LaKiojB06FCYTCZ8+OGHKCsrw0cffXQn34mIiIjukM2JwbFjxzB27FjpOC0tDQCQnJyM7OxsAMC2bdsghMBzzz1ntY6qqir06vWfzor6+nrMnTsXBoMBGo0G4eHhKCoqwiOPPGJr84iIiORxHwNZKiFc49uZTCZoNBrEYhL6qNy6uzlERGSDG6IZBdgNo9Fo8RpjR2r9PfHkvfPQp5d9k9dvmBuRf+ljp7a3u/BdCUREpCzsMZDF1y4TERGRhD0GRESkLAIO6DFwSEt6JCYGRESkLBxKkMWhBCIiIpKwx4CIiJTFbAZg5wZFZm5wRERE5Bo4lCCLQwlEREQkYY8BEREpC3sMZDExICIiZeHbFWVxKIGIiIgk7DEgIiJFEcIMYedrk+29vydjYkBERMoihP1DAZxjQERE5CKEA+YYuHBiwDkGREREJGGPARERKYvZDKjsnCPAOQZEREQugkMJsjiUQERERBL2GBARkaIIsxnCzqEELlckIiJyFRxKkMWhBCIiIpKwx4CIiJTFLAAVewzaw8SAiIiURQgA9i5XdN3EgEMJREREJGGPARERKYowCwg7hxKEC/cYMDEgIiJlEWbYP5TgussVOZRARESKIszCIcVWq1evRlBQEDw8PBAREYEvvvjCCd/OfkwMiIiInGz79u1ITU3Fm2++iZMnT+Kxxx5DQkICqqqqurtpbbjMUELreM8NNNu9bwUREXWtG2gG0DVj9zdEo91DAa3tNZlMFufVajXUanWb+A8++ABz5szBr371KwDAihUrsH//fqxZswaZmZl2tcXRXCYxuHr1KgCgGPu6uSVERHSnrl69Co1G45S63d3dodVqUWxwzO+Jn/3sZwgICLA498477yAjI8PiXFNTE44fP47FixdbnI+Li0NJSYlD2uJILpMY6PV6VFdXw8vLCyqVymqMyWRCQEAAqqur4e3t3cUtvHNsd9e7W9vOdncttttxhBC4evUq9Hq90z7Dw8MDlZWVaGpqckh9Qog2v2+s9Rb8+OOPaGlpgb+/v8V5f39/GAwGh7TFkVwmMejVqxcGDRrUqVhvb+8e8z8GW7DdXe9ubTvb3bXYbsdwVk/BrTw8PODh4eH0z7Hm9iTCWmLRE3DyIRERkRP5+vqid+/ebXoHamtr2/Qi9ARMDIiIiJzI3d0dERERyMvLszifl5eHmJiYbmpV+1xmKKEz1Go13nnnHatjQD0Z29317ta2s91di+2mzkpLS0NSUhIiIyMRHR2NTz75BFVVVZg3b153N60NlXDlfR2JiIh6iNWrV+P9999HTU0NQkJCsHz5cjz++OPd3aw2mBgQERGRhHMMiIiISMLEgIiIiCRMDIiIiEjCxICIiIgkLpcY2Ppay8LCQkRERMDDwwP3338/Pv744y5q6U2ZmZkYOXIkvLy84Ofnh8mTJ6O8vFz2noKCAqhUqjblm2++6aJWAxkZGW0+X6vVyt7T3c+61X333Wf1+c2fP99qfHc976KiIkyYMAF6vR4qlQq7du2yuC6EQEZGBvR6Pfr27YvY2FicOXOmw3p37NiBhx56CGq1Gg899BBycnK6rN3Nzc14/fXXERoaCk9PT+j1esyaNQvff/+9bJ3Z2dlW/wyuX7/eJe0GgNmzZ7f5/KioqA7r7c7nDcDqc1OpVPjDH/7Qbp1d8byp53KpxMDW11pWVlbi6aefxmOPPYaTJ0/ijTfewKuvvoodO3Z0WZsLCwsxf/58HD58GHl5ebhx4wbi4uJw7dq1Du8tLy9HTU2NVIYOHdoFLf6P4cOHW3z+qVOn2o3tCc+61dGjRy3a3brpyC9/+UvZ+7r6eV+7dg1hYWFYtWqV1evvv/8+PvjgA6xatQpHjx6FVqvFU089Jb1QzJrS0lJMnz4dSUlJ+Prrr5GUlIRp06bhyJEjXdLun376CSdOnMBbb72FEydOYOfOnfjHP/6BiRMndlivt7e3xfOvqalx6Na2HT1vABg/frzF5+/bJ/8ynu5+3gDaPLP169dDpVJh6tSpsvU6+3lTDyZcyCOPPCLmzZtnce6BBx4Qixcvthq/aNEi8cADD1ice+mll0RUVJTT2tiR2tpaAUAUFha2G3Po0CEBQNTV1XVdw27zzjvviLCwsE7H98Rn3erXv/61GDJkiDCbzVav94TnDUDk5ORIx2azWWi1WrF06VLp3PXr14VGoxEff/xxu/VMmzZNjB8/3uJcfHy8mDFjhsPbLETbdlvz1VdfCQDiwoUL7cZs2LBBaDQaxzZOhrV2Jycni0mTJtlUT0983pMmTRJPPPGEbExXP2/qWVymx6D1tZZxcXEW5+Vea1laWtomPj4+HseOHUNzc7PT2irHaDQCAHx8fDqMDQ8Ph06nw7hx43Do0CFnN62NiooK6PV6BAUFYcaMGfj222/bje2Jzxq4+fdmy5YtePHFFzt8mUl3P+9bVVZWwmAwWDxTtVqNMWPGyL7Gtb0/h+589avRaIRKpcI999wjG9fQ0IDAwEAMGjQIiYmJOHnyZNc08BYFBQXw8/PDsGHDkJKSgtraWtn4nva8L1++jL1792LOnDkdxvaE503dw2USgzt5raXBYLAaf+PGDfz4449Oa2t7hBBIS0vDo48+ipCQkHbjdDodPvnkE+zYsQM7d+5EcHAwxo0bh6Kioi5r66hRo7Bp0ybs378f69atg8FgQExMDK5cuWI1vqc961a7du1CfX09Zs+e3W5MT3jet2v9O23ra1zb+3Porle/Xr9+HYsXL8bMmTNl3/L3wAMPIDs7G3v27MFnn30GDw8PjB49GhUVFV3W1oSEBGzduhUHDx7EsmXLcPToUTzxxBNobGxs956e9rw3btwILy8vTJkyRTauJzxv6j4u964EW19raS3e2vmu8Morr+Dvf/87iouLZeOCg4MRHBwsHUdHR6O6uhp//OMfu2x7zYSEBOnn0NBQREdHY8iQIdi4cSPS0tKs3tOTnnWrrKwsJCQkyL4Dvic87/bcyWtce8qrX5ubmzFjxgyYzWasXr1aNjYqKspiot/o0aMxYsQIrFy5Eh9++KGzmwoAmD59uvRzSEgIIiMjERgYiL1798r+ou0pzxsA1q9fj+eff77DuQI94XlT93GZHoM7ea2lVqu1Gt+nTx8MGDDAaW21ZsGCBdizZw8OHTqEQYMG2Xx/VFRUt2bznp6eCA0NbbcNPelZt7pw4QLy8/Pxq1/9yuZ7u/t5t64AsfU1ru39OXT1q1+bm5sxbdo0VFZWIi8vT7a3wJpevXph5MiR3fpnoNPpEBgYKNuGnvK8AeCLL75AeXn5Hf197wnPm7qOyyQGd/Jay+jo6DbxBw4cQGRkJNzc3JzW1lsJIfDKK69g586dOHjwIIKCgu6onpMnT0Kn0zm4dZ3X2NiIc+fOtduGnvCsb7dhwwb4+fnhF7/4hc33dvfzDgoKglartXimTU1NKCwslH2Na3t/Dl356tfWpKCiogL5+fl3lBgKIVBWVtatfwZXrlxBdXW1bBt6wvNulZWVhYiICISFhdl8b0943tSFumvWozNs27ZNuLm5iaysLHH27FmRmpoqPD09xXfffSeEEGLx4sUiKSlJiv/2229Fv379xGuvvSbOnj0rsrKyhJubm/i///u/Lmvzf/3XfwmNRiMKCgpETU2NVH766Scp5vZ2L1++XOTk5Ih//OMf4vTp02Lx4sUCgNixY0eXtXvhwoWioKBAfPvtt+Lw4cMiMTFReHl59ehnfauWlhYxePBg8frrr7e51lOe99WrV8XJkyfFyZMnBQDxwQcfiJMnT0qz95cuXSo0Go3YuXOnOHXqlHjuueeETqcTJpNJqiMpKcliVc6XX34pevfuLZYuXSrOnTsnli5dKvr06SMOHz7cJe1ubm4WEydOFIMGDRJlZWUWf+cbGxvbbXdGRobIzc0V//znP8XJkyfFCy+8IPr06SOOHDnSJe2+evWqWLhwoSgpKRGVlZXi0KFDIjo6Wtx77709+nm3MhqNol+/fmLNmjVW6+iO5009l0slBkII8dFHH4nAwEDh7u4uRowYYbHsLzk5WYwZM8YivqCgQISHhwt3d3dx3333tfs/HGcBYLVs2LCh3Xa/9957YsiQIcLDw0P0799fPProo2Lv3r1d2u7p06cLnU4n3NzchF6vF1OmTBFnzpxpt81CdP+zvtX+/fsFAFFeXt7mWk953q3LJG8vycnJQoibSxbfeecdodVqhVqtFo8//rg4deqURR1jxoyR4lv9+c9/FsHBwcLNzU088MADDk9w5NpdWVnZ7t/5Q4cOtdvu1NRUMXjwYOHu7i4GDhwo4uLiRElJSZe1+6effhJxcXFi4MCBws3NTQwePFgkJyeLqqoqizp62vNutXbtWtG3b19RX19vtY7ueN7Uc/G1y0RERCRxmTkGREREZD8mBkRERCRhYkBEREQSJgZEREQkYWJAREREEiYGREREJGFiQERERBImBkRERCRhYkBEREQSJgZEREQkYWJAREREkv8PvJqzF8cRDSIAAAAASUVORK5CYII=",
|
|
"text/plain": [
|
|
"<Figure size 640x480 with 2 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"records = []\n",
|
|
"\n",
|
|
"for i in range(1000):\n",
|
|
" C_t = simulate(C_t)\n",
|
|
" records.append(C_t)\n",
|
|
"\n",
|
|
"# with np.printoptions(precision=2,floatmode='fixed',linewidth=200):\n",
|
|
"# print(C_t)\n",
|
|
"\n",
|
|
"plt.imshow(C_t, vmin=0, vmax=2000)\n",
|
|
"plt.colorbar()\n",
|
|
"plt.show()\n"
|
|
]
|
|
},
|
|
{
|
|
"attachments": {},
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Widget"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 105,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "2f9cde999c5849f9b782f7e68bf424c9",
|
|
"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": 105,
|
|
"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=0.1)\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
|
|
}
|