mirror of
https://git.gfz-potsdam.de/naaice/tug.git
synced 2025-12-15 18:38:23 +01:00
490 lines
82 KiB
Plaintext
490 lines
82 KiB
Plaintext
{
|
|
"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": [
|
|
"# 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": 12,
|
|
"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",
|
|
"# C_t[19,19] = 2000\n",
|
|
"\n",
|
|
"\n",
|
|
"plt.imshow(C_t, vmin=0, vmax=2000)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 11,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# calculate mean between two cells\n",
|
|
"def alpha_interblock(alpha1, alpha2, harmonic=True):\n",
|
|
" if not harmonic:\n",
|
|
" return 0.5 * (alpha1 + alpha2)\n",
|
|
" else:\n",
|
|
" return 2 / ((1/alpha1) + (1/alpha2))\n"
|
|
]
|
|
},
|
|
{
|
|
"attachments": {},
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Implementation for constant boundary conditions"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# simulate one time step\n",
|
|
"def simulate_constant_boundary(C_t): \n",
|
|
" C_t1 = np.zeros((grid_size['x'],grid_size['y']))\n",
|
|
"\n",
|
|
" # inner cells\n",
|
|
" for i in range(1, grid_size['y']-1): #rows\n",
|
|
" for j in range(1, grid_size['x']-1): #columns\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) \n",
|
|
" + 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]) \n",
|
|
" + 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]) \n",
|
|
" + 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]) \n",
|
|
" + 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]) \n",
|
|
" + 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]) \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
|
|
"\n",
|
|
" # corner top left\n",
|
|
" i = 0\n",
|
|
" j = i\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1] \n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + 2 * alpha_x[i,j]) * C_t[i,j]\n",
|
|
" + 2 * alpha_x[i,j] * bc_left) \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * C_t[1,j] \n",
|
|
" - (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) + 2 * alpha_y[0,j]) * C_t[0,j]\n",
|
|
" + 2 * alpha_y[0,j] * bc_top)\n",
|
|
" \n",
|
|
"\n",
|
|
" # corner top right\n",
|
|
" i = 0\n",
|
|
" j = grid_size['x']-1\n",
|
|
" n = j\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (2 * alpha_x[i,n] * bc_right \n",
|
|
" - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) + 2 * alpha_x[i,n]) * C_t[i,n]\n",
|
|
" + alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * C_t[i,n-1]) \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * C_t[1,j] \n",
|
|
" - (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) + 2 * alpha_y[0,j]) * C_t[0,j]\n",
|
|
" + 2 * alpha_y[0,j] * bc_top)\n",
|
|
"\n",
|
|
" # corner bottom left\n",
|
|
" i = grid_size['y']-1\n",
|
|
" m = i\n",
|
|
" j = 0\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1] \n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + 2 * alpha_x[i,j]) * C_t[i,j]\n",
|
|
" + 2 * alpha_x[i,j] * bc_left) \\\n",
|
|
" + time_step/delta_y**2 * (2 * alpha_y[m,j] * bc_bottom \n",
|
|
" - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) + 2 * alpha_y[m,j]) * C_t[m,j]\n",
|
|
" + alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * C_t[m-1,j])\n",
|
|
"\n",
|
|
" # corner bottom right\n",
|
|
" i = grid_size['y']-1\n",
|
|
" j = grid_size['x']-1\n",
|
|
" m = i \n",
|
|
" n = j\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (2 * alpha_x[i,n] * bc_right \n",
|
|
" - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) + 2 * alpha_x[i,n]) * C_t[i,n]\n",
|
|
" + alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * C_t[i,n-1]) \\\n",
|
|
" + time_step/delta_y**2 * (2 * alpha_y[m,j] * bc_bottom \n",
|
|
" - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) + 2 * alpha_y[m,j]) * C_t[m,j]\n",
|
|
" + alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * C_t[m-1,j])\n",
|
|
"\n",
|
|
" return C_t1"
|
|
]
|
|
},
|
|
{
|
|
"attachments": {},
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Implementation for closed boundary conditions"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 115,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# simulate one time step\n",
|
|
"def simulate_closed_boundary(C_t): \n",
|
|
" C_t1 = np.zeros((grid_size['x'],grid_size['y']))\n",
|
|
"\n",
|
|
" # inner cells\n",
|
|
" for i in range(1, grid_size['y']-1): #rows\n",
|
|
" for j in range(1, grid_size['x']-1): #columns\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j]) \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
|
|
" \n",
|
|
" # boundary conditions\n",
|
|
" # left without corners / looping over rows\n",
|
|
" for i in range(1, grid_size['y']-1):\n",
|
|
" j = 0\n",
|
|
" C_t1[i,0] = C_t[i,0] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * (C_t[i,j+1] - C_t[i,j]))\\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j]) \n",
|
|
"\n",
|
|
" # right without corners / looping over rows\n",
|
|
" n = grid_size['x']-1 # maximum index in x-direction (columns)\n",
|
|
" for i in range(1, grid_size['y']-1):\n",
|
|
" j = n\n",
|
|
" C_t1[i,n] = C_t[i,n] \\\n",
|
|
" + time_step/delta_x**2 * ( \n",
|
|
" - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * (C_t[i,n] - C_t[i,n-1])))\\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j])\n",
|
|
"\n",
|
|
" # top without corners / looping over columns\n",
|
|
" for j in range(1, grid_size['x']-1):\n",
|
|
" i = 0\n",
|
|
" C_t1[0,j] = C_t[0,j] \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * (C_t[1,j] - C_t[0,j])) \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
|
|
"\n",
|
|
" # bottom without corners / looping over columns\n",
|
|
" m = grid_size['y']-1 # maximum index in y-direction (rows)\n",
|
|
" for j in range(1, grid_size['x']-1):\n",
|
|
" i = m\n",
|
|
" C_t1[m,j] = C_t[m,j] \\\n",
|
|
" + time_step/delta_y**2 * ( \n",
|
|
" - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * (C_t[m,j] - C_t[m-1,j]))) \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
|
|
"\n",
|
|
" # corner top left\n",
|
|
" i = 0\n",
|
|
" j = i\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * (C_t[i,j+1] - C_t[i,j]))\\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * (C_t[1,j]-C_t[0,j]))\n",
|
|
" \n",
|
|
"\n",
|
|
" # corner top right\n",
|
|
" i = 0\n",
|
|
" j = grid_size['x']-1\n",
|
|
" n = j\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (- (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * (C_t[i,n] - C_t[i,n-1])))\\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * (C_t[1,j] - C_t[0,j])) \n",
|
|
" \n",
|
|
"\n",
|
|
" # corner bottom left\n",
|
|
" i = grid_size['y']-1\n",
|
|
" m = i\n",
|
|
" j = 0\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * (C_t[i,j+1] - C_t[i,j])) \\\n",
|
|
" + time_step/delta_y**2 * (- (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j] * (C_t[m,j] - C_t[m-1,j]))))\n",
|
|
"\n",
|
|
" # corner bottom right\n",
|
|
" i = grid_size['y']-1\n",
|
|
" j = grid_size['x']-1\n",
|
|
" m = i \n",
|
|
" n = j\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (- alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * (C_t[i,n] - C_t[i,n-1]))\\\n",
|
|
" + time_step/delta_y**2 * (- alpha_interblock(alpha_y[m,j], alpha_y[m-1,j] * (C_t[m,j] - C_t[m-1,j])))\n",
|
|
"\n",
|
|
" return C_t1"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 13,
|
|
"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(100):\n",
|
|
" C_t = simulate_constant_boundary(C_t)\n",
|
|
" records.append(C_t)\n",
|
|
"\n",
|
|
"# with np.printoptions(precision=2,floatmode='fixed',linewidth=200):\n",
|
|
"# print(C_t)\n",
|
|
"\n",
|
|
"plt.imshow(C_t, vmin=0, vmax=2000)\n",
|
|
"plt.colorbar()\n",
|
|
"plt.show()\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 14,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgIAAAGdCAYAAABgnRvHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA96UlEQVR4nO3df1xUZb4H8M8ZYGbMmGlNZSARyZf5e43QBAzLTIzS3NSVagPdtNab7WZsr4rMTdt7o3bLSPNH7qrktiF18de9sqv4SkRX8yZB25YZbSRksCy+VgY0ZmDmuX8Ysx6ZGRjnGRjnfN69zuvVnHnOl+8cBuc7z/Oc8yhCCAEiIiLSJF1vJ0BERES9h4UAERGRhrEQICIi0jAWAkRERBrGQoCIiEjDWAgQERFpGAsBIiIiDWMhQEREpGHhvZ2ALE6nE99++y0iIyOhKEpvp0NERD4QQqC5uRkxMTHQ6QL3HbW1tRV2u11KLL1eD6PRKCVWbwqZQuDbb79FbGxsb6dBRER+qK2txaBBgwISu7W1FfFxV6O+wSElnsViQXV19RVfDIRMIRAZGQkAmHzVXIQrEX7F0g24VkZKAACnqa+UOI4+8n5VQienx0Rxyrs7ta5Nzh+m7rt2KXEAQPmuVU6g72xy4gAQNjk5CXublDgAICT97oRDTpwLwZyS4vAO7D2lHW04jGLXv+WBYLfbUd/gQHV5HEyR/vU6WJudiE88BbvdzkIgWHQMB4QrEQhX9H7F0ukMMlICADjD5MRSwv0rbi4WlIWAU1IhECbvA07RSXp9suIAEIqcDzghcfhMKHKKL6HI7A6WVAiAhUCP+f5U98TQrilS53chEEpCphAgIiLqDodwwuFnjeeQ1esUBFgIEBGRpjgh4PSzt8ff44NJwPpG1q1bh/j4eBiNRiQmJuLQoUNe2x88eBCJiYkwGo24/vrrsWHDhkClRkREGuaU9F+oCEghUFhYiKVLl2LZsmWoqKhAamoq0tPTUVNT47Z9dXU17rrrLqSmpqKiogLPPvssfvGLX6CoqCgQ6REREdH3AlIIrFq1CgsXLsSiRYswcuRI5OXlITY2FuvXr3fbfsOGDRg8eDDy8vIwcuRILFq0CA899BBeeeWVQKRHREQa5hBCyhYqpBcCdrsd5eXlSEtLU+1PS0vDkSNH3B5z9OjRTu2nT5+O48ePo63N/Sxwm80Gq9Wq2oiIiLrSMUfA3y1USC8EGhsb4XA4EBUVpdofFRWF+vp6t8fU19e7bd/e3o7Gxka3x+Tm5sJsNrs23kyIiIjIdwGbLHjptaBCCK/Xh7pr725/h5ycHDQ1Nbm22tpaPzMmIiItcELA4ecWSj0C0i8f7N+/P8LCwjp9+29oaOj0rb+DxWJx2z48PBzXXuv+Ln8GgwEGg7wb/xARkTbw8kE16T0Cer0eiYmJKCkpUe0vKSlBSkqK22OSk5M7td+3bx/Gjx+PiAh5d9QjIiIitYAMDWRnZ+P3v/89Nm/ejBMnTuCJJ55ATU0NFi9eDOBCt35WVpar/eLFi3Hq1ClkZ2fjxIkT2Lx5MzZt2oQnn3wyEOkREZGG8aoBtYDcWTAjIwNnzpzBCy+8gLq6OowZMwbFxcWIi4sDANTV1anuKRAfH4/i4mI88cQTWLt2LWJiYrB69WrMmTMnEOkREZGGOeH/ahShczshQBEiNMoaq9UKs9mM2/ve7/+iQwP7S8oKcJolrT54VYgvOmSXtfqgxEWHzstafVBSHACiNRhXH5S06BBXH9S0dtGGUuxCU1MTTCZTQH5Gx+fE5yeiEOnnokPNzU6MGPmPgObbU7jWABERaUrHzH9/Y4SKkCsEdAOu9XsZ4bMToiVlA5wdKmcaRusAeR1RQi/nDazY5S0XqrfKOU9G97eduCxXNcg5533+aZcSBwAi/nleShxdU4uUOAAgzsvJScjsOWEvBXnhEJCw+qCcXIJByBUCRERE3nCOgFrAbihEREREwY89AkREpClOKHDAv6FNp5/HBxMWAkREpClOcWHzN0ao4NAAERGRhrFHgIiINMUhYWjA3+ODCQsBIiLSFBYCahwaICIi0jD2CBARkaY4hQKn8POqAT+PDyYsBIiISFM4NKDGoQEiIiINY48AERFpigM6OPz8HixxFYpex0KAiIg0RUiYIyA4R4CIiOjKxDkCapwjQEREpGHsESAiIk1xCB0cws85AiG01gALASIi0hQnFDj97BB3InQqAQ4NEBERaVjI9Qg4TX3hDDP4FePsUHn1kUi0Sokz6boaKXEAIMrQLCVOo/1qKXEA4G+N0VLiNH5zjZQ4AGA/JefPw6HXS4kDALLOuF7I+zajOJ1yAjkkXpAlKychKQ4AEUrXm13hOFlQLeQKASIiIm/kzBHg0AARERGFAPYIEBGRplyYLOjnokMcGiAiIroyOSXcYphXDRAREVFIYI8AERFpCicLqrEQICIiTXFCxxsKXYSFABERaYpDKHD4uXqgv8cHE84RICIi0jD2CBARkaY4JFw14AihoQH2CBARkaY4hU7K5qt169YhPj4eRqMRiYmJOHTokMe227dvx7Rp0zBgwACYTCYkJydj7969qjb5+flQFKXT1tra6lNeLASIiIgCrLCwEEuXLsWyZctQUVGB1NRUpKeno6bG/ToyZWVlmDZtGoqLi1FeXo4pU6Zg5syZqKioULUzmUyoq6tTbUaj0afcODRARESa0htDA6tWrcLChQuxaNEiAEBeXh727t2L9evXIzc3t1P7vLw81eMXX3wRu3btwv/8z/8gISHBtV9RFFgsFt9fwEXYI0BERJrixL+vHLjcrWNdSqvVqtpsNlunn2e321FeXo60tDTV/rS0NBw5cqR7OTudaG5uRr9+/VT7W1paEBcXh0GDBmHGjBmdegy6Q3ohkJubiwkTJiAyMhIDBw7Ej370I5w8edLrMaWlpW7HOT7//HPZ6REREUkTGxsLs9ns2tx9u29sbITD4UBUVJRqf1RUFOrr67v1c1599VWcO3cO8+bNc+0bMWIE8vPzsXv3bhQUFMBoNGLSpEmoqqry6TVIHxo4ePAglixZggkTJqC9vR3Lli1DWloaPvvsM/Tt29frsSdPnoTJZHI9HjBggOz0iIhI4+TcUOjC8bW1tarPLYPB4PEYRVHfe0AI0WmfOwUFBVixYgV27dqFgQMHuvYnJSUhKSnJ9XjSpEm46aabsGbNGqxevbrbr0V6IfDnP/9Z9XjLli0YOHAgysvLMXnyZK/HDhw4ENdcc43slIiIiFzk3GL4wvEmk0lVCLjTv39/hIWFdfr239DQ0KmX4FKFhYVYuHAh3nvvPdxxxx1e2+p0OkyYMKH3ewQu1dTUBACdxjXcSUhIQGtrK0aNGoXnnnsOU6ZM8djWZrOpxmKsVisAwNEnHEp4hF85tw5wdt2omyZd535GqK8eGnhYShwAuD7CKiXO1+1XS4kDAH82/FBKnD85RkmJAwDN566REkdvlTcCZ2iS8ycb3uzf38jFws5J+mckLExOHADQcfoTBQ+9Xo/ExESUlJTg3nvvde0vKSnBrFmzPB5XUFCAhx56CAUFBbj77ru7/DlCCFRWVmLs2LE+5RfQQkAIgezsbNxyyy0YM2aMx3bR0dHYuHEjEhMTYbPZ8Ic//AFTp05FaWmpx16E3NxcrFy5MlCpExFRiHJCgRP+3SLY1+Ozs7ORmZmJ8ePHIzk5GRs3bkRNTQ0WL14MAMjJycHp06exdetWABeKgKysLLz++utISkpy9Sb06dMHZrMZALBy5UokJSVh2LBhsFqtWL16NSorK7F27VqfcgtoIfDYY4/hr3/9Kw4f9v5tdvjw4Rg+fLjrcXJyMmpra/HKK694LARycnKQnZ3temy1WhEbGysncSIiClkyhwa6KyMjA2fOnMELL7yAuro6jBkzBsXFxYiLiwMA1NXVqe4p8Oabb6K9vR1LlizBkiVLXPvnz5+P/Px8AMDZs2fxyCOPoL6+HmazGQkJCSgrK8PNN9/sU24BKwR+/vOfY/fu3SgrK8OgQYN8Pj4pKQlvv/22x+cNBoPXSRlERETuyLmPgO/HP/roo3j00UfdPtfx4d6htLS0y3ivvfYaXnvtNZ/zuJT0QkAIgZ///OfYsWMHSktLER8ff1lxKioqEB0dLTk7IiIiupj0QmDJkiV45513sGvXLkRGRrrGNcxmM/r06QOg81hIXl4ehgwZgtGjR8Nut+Ptt99GUVERioqKZKdHREQa5xQKnH4uI+zv8cFEeiGwfv16AMBtt92m2r9lyxYsWLAAQOexELvdjieffBKnT59Gnz59MHr0aOzZswd33XWX7PSIiEjjnBKGBvy9D0EwCcjQQFcuHQt56qmn8NRTT8lOhYiIiLrARYeIiEhTLncZ4UtjhAoWAkREpCkOKHD4eR8Bf48PJqFT0hAREZHP2CNARESawqEBNRYCRESkKQ7437XvkJNKUAidkoaIiIh8xh4BIiLSFA4NqLEQICIiTemNRYeCGQsBIiLSFCFhGWLByweJiIgoFLBHgIiINIVDA2ohVwgInQKh87PLR9/1egndFWVolhLn+girlDgAMDj8akmRWiTFAaL1TVLiRBptUuIAQJOk94EzQkqYC7HCJXVHKhK7NWXFkpkTkRdcfVAtdEoaIiIi8lnI9QgQERF545CwDLG/xwcTFgJERKQpHBpQC52ShoiIiHzGHgEiItIUJ3Rw+vk92N/jgwkLASIi0hSHUODws2vf3+ODSeiUNEREROQz9ggQEZGmcLKgGgsBIiLSFCFh9UHBOwsSERFdmRxQ4PBz0SB/jw8moVPSEBERkc/YI0BERJriFP6P8TvlLUnT61gIEBGRpjglzBHw9/hgEjqvhIiIiHzGHgEiItIUJxQ4/Zzs5+/xwYSFABERaQrvLKjGoQEiIiINC7keAcUpoPg5nVOxy6v0Gu1XS4nzdbucOBe0SIkiM6c6u1lKnOZWg5Q4AKCzyXkf6NqkhPk+lqSpyk6nnDgAICTlJCsOURc4WVAt5AoBIiIib5yQcIvhEJojEDolDREREfmMPQJERKQpQsJVAyKEegRYCBARkaZw9UE1FgJERKQpnCyoJv2VrFixAoqiqDaLxeL1mIMHDyIxMRFGoxHXX389NmzYIDstIiIiciMgPQKjR4/G/v37XY/DwsI8tq2ursZdd92Fhx9+GG+//Tb+8pe/4NFHH8WAAQMwZ86cQKRHREQaxqEBtYAUAuHh4V32AnTYsGEDBg8ejLy8PADAyJEjcfz4cbzyyissBIiISDreYlgtIIMcVVVViImJQXx8PO677z589dVXHtsePXoUaWlpqn3Tp0/H8ePH0dbm+U4sNpsNVqtVtREREZFvpBcCEydOxNatW7F371787ne/Q319PVJSUnDmzBm37evr6xEVFaXaFxUVhfb2djQ2Nnr8Obm5uTCbza4tNjZW6usgIqLQ1DE04O8WKqQXAunp6ZgzZw7Gjh2LO+64A3v27AEAvPXWWx6PURT1CRXf32r00v0Xy8nJQVNTk2urra2VkD0REYU6FgJqAb98sG/fvhg7diyqqqrcPm+xWFBfX6/a19DQgPDwcFx77bUe4xoMBhgM8u4rT0REpEUBvxDSZrPhxIkTiI6Odvt8cnIySkpKVPv27duH8ePHIyIiItDpERGRxrBHQE16IfDkk0/i4MGDqK6uxrFjxzB37lxYrVbMnz8fwIUu/aysLFf7xYsX49SpU8jOzsaJEyewefNmbNq0CU8++aTs1IiIiFgIXEL60MA333yD+++/H42NjRgwYACSkpLwwQcfIC4uDgBQV1eHmpoaV/v4+HgUFxfjiSeewNq1axETE4PVq1fz0kEiIqIeIL0Q2LZtm9fn8/PzO+279dZb8dFHH8lOhYiIqBMB/+8DIOSkEhS41gAREWkK7yyoxkKAiIg0hYWAWsgVAro2B3ROh18x9FZ5cyj/1uj+aglf/dnwQylxACBa3yQlTp3dLCUOABw7M0RKnLNnrpYSBwD6NMl5H0S0yOtEDP/Ov/d2B8UmJw4AoF1SLIfEnJxOebGIJFm3bh1++9vfoq6uDqNHj0ZeXh5SU1Pdtt2+fTvWr1+PyspK2Gw2jB49GitWrMD06dNV7YqKirB8+XL8/e9/x9ChQ/Ff//VfuPfee33KK3TWUSQiIuqG3rhqoLCwEEuXLsWyZctQUVGB1NRUpKenqybPX6ysrAzTpk1DcXExysvLMWXKFMycORMVFRWuNkePHkVGRgYyMzPx8ccfIzMzE/PmzcOxY8d8yk0RHbfxu8JZrVaYzWZMuekZhIcZ/Yp1aoZJUlbA1eM93ybZF9Ou+1xKHCC0ewS+qhkoJQ4A9KnWS4lzdY28P7HIb+xS4uj/cU5KHADQnW2WEke0tEiJAwCi1SYnTnu7lDgAIGT1eITGP9mdtIs2lGIXmpqaYDLJ+zf4Yh2fE7fsXoLwvv7dkK79nA2H71nb7XwnTpyIm266CevXr3ftGzlyJH70ox8hNze3Wz9z9OjRyMjIwK9+9SsAQEZGBqxWK/70pz+52tx55534wQ9+gIKCgm6/FvYIEBERBZDdbkd5eXmnBfbS0tJw5MiRbsVwOp1obm5Gv379XPs8LdrX3ZgdQm6OABERkTdCKBB+TvbrOP7SlW/d3f6+sbERDofD7QJ7l95i35NXX30V586dw7x581z7PC3a192YHdgjQEREmuKEImUDgNjYWNVKuN66+d0tsOdtcb0OBQUFWLFiBQoLCzFwoHr483JjXow9AkRERJeptrZWNUfA3WJ4/fv3R1hYmNsF9i79Rn+pwsJCLFy4EO+99x7uuOMO1XOeFu3rKual2CNARESaIvOqAZPJpNrcFQJ6vR6JiYmdFtgrKSlBSkqKxzwLCgqwYMECvPPOO7j77rs7Pe9p0T5vMd1hjwAREWmKzDkC3ZWdnY3MzEyMHz8eycnJ2LhxI2pqarB48WIAFxbkO336NLZu3QrgQhGQlZWF119/HUlJSa5v/n369IHZfOGKrccffxyTJ0/Gyy+/jFmzZmHXrl3Yv38/Dh8+7FNu7BEgIiIKsIyMDOTl5eGFF17AjTfeiLKyMhQXF3tckO/NN99Ee3s7lixZgujoaNf2+OOPu9qkpKRg27Zt2LJlC374wx8iPz8fhYWFmDhxok+5sUeAiIg0pbduMfzoo4/i0UcfdfvcpQvylZaWdivm3LlzMXfuXJ9zuRgLASIi0pTeGBoIZiwEiIhIU4SEHoFQKgQ4R4CIiEjD2CNARESaIuD/kg2htOIDCwEiItIUJxQo8HOyoJ/HBxMODRAREWkYewSIiEhTeNWAWsgVArrv2qELa/MrhrFRUjIAGr+5RkqcPzlGSYkDAJFGOWu1N7f6t573xc6euVpKHH1dhJQ4AGA4IyeOsUnSOvQAws/5997uoNjsUuIAgGiTk5NwOKXEAS4svCIljjOURoKpg1MoUHrhPgLBikMDREREGhZyPQJERETeCCHhqoEQ6ixiIUBERJrCOQJqHBogIiLSMPYIEBGRprBHQI2FABERaQqvGlBjIUBERJrCyYJqnCNARESkYewRICIiTbnQI+DvHAFJyQQBFgJERKQpnCyoxqEBIiIiDWOPABERaYr4fvM3RqhgIUBERJrCoQE1Dg0QERFpmPRCYMiQIVAUpdO2ZMkSt+1LS0vdtv/8889lp0ZERPTvsQF/txAhfWjgww8/hMPx7/XX//a3v2HatGn48Y9/7PW4kydPwmQyuR4PGDBAdmpERESAhKEBhNDQgPRC4NIP8JdeeglDhw7Frbfe6vW4gQMH4pprrpGdDhERkQrvLKgW0DkCdrsdb7/9Nh566CEoivfqKSEhAdHR0Zg6dSoOHDgQyLSIiIjoewG9amDnzp04e/YsFixY4LFNdHQ0Nm7ciMTERNhsNvzhD3/A1KlTUVpaismTJ3s8zmazwWazuR5brVYAgPJdKxSdf6XaVQ1Ov46/mP2UnFPcfO4aKXEAoEkvp5TV2eR1jfVpklOTGs5ICQMA6Nvg6LpRNxj+1SYlDgCENdu6btQNyndy4gCAsEt6fW3yzhMccn53EPL+LQipr5BXOF41oBbQQmDTpk1IT09HTEyMxzbDhw/H8OHDXY+Tk5NRW1uLV155xWshkJubi5UrV0rNl4iINEAo/o/xh1AhELChgVOnTmH//v1YtGiRz8cmJSWhqqrKa5ucnBw0NTW5ttra2stNlYiISLMC1iOwZcsWDBw4EHfffbfPx1ZUVCA6OtprG4PBAIPBcLnpERGRRnGyoFpACgGn04ktW7Zg/vz5CA9X/4icnBycPn0aW7duBQDk5eVhyJAhGD16tGtyYVFREYqKigKRGhERaR3vMawSkEJg//79qKmpwUMPPdTpubq6OtTU1Lge2+12PPnkkzh9+jT69OmD0aNHY8+ePbjrrrsCkRoRERFdJCCFQFpaGoSHfpP8/HzV46eeegpPPfVUINIgIiLqhFcNqHHRISIi0p4Q6tr3FxcdIiIi0jD2CBARkaZwaECNhQAREWkLrxpQYSFAREQao3y/+RsjNHCOABERkYaxR4CIiLSFQwMqLASIiEhbWAiocGiAiIhIw9gjQERE2sJliFVYCBARkaZw9UG10CsEvrMBOv9+Q33+aZeUDODQ66XE0VvljeI4I+TE0bXJiQMAES1y/qqMTQ4pcQDA8C85LzD8X99JiQMASst5KXFEa6uUOAAg7HL+XoTDKSUOAAhnCP0rTRRgoVcIEBERecPJgiosBIiISFs4R0CFVw0QERFpGHsEiIhIUxRxYfM3RqhgIUBERNrCOQIqLASIiEhbOEdAhXMEiIiINIw9AkREpC0cGlBhjwAREWmLkLT5aN26dYiPj4fRaERiYiIOHTrksW1dXR0eeOABDB8+HDqdDkuXLu3UJj8/H4qidNpafbxhGAsBIiKiACssLMTSpUuxbNkyVFRUIDU1Fenp6aipqXHb3mazYcCAAVi2bBnGjRvnMa7JZEJdXZ1qMxqNPuXGQoCIiLSlF3oEVq1ahYULF2LRokUYOXIk8vLyEBsbi/Xr17ttP2TIELz++uvIysqC2Wz2GFdRFFgsFtXmKxYCRESkLR1XDfi7AbBararNZrN1+nF2ux3l5eVIS0tT7U9LS8ORI0f8eiktLS2Ii4vDoEGDMGPGDFRUVPgcg4UAERHRZYqNjYXZbHZtubm5ndo0NjbC4XAgKipKtT8qKgr19fWX/bNHjBiB/Px87N69GwUFBTAajZg0aRKqqqp8isOrBoiISFNk3lmwtrYWJpPJtd9gMHg+RlHfe0AI0WmfL5KSkpCUlOR6PGnSJNx0001Ys2YNVq9e3e04LASIiEhbJF4+aDKZVIWAO/3790dYWFinb/8NDQ2degn8odPpMGHCBJ97BDg0QEREFEB6vR6JiYkoKSlR7S8pKUFKSoq0nyOEQGVlJaKjo306jj0CREREAZadnY3MzEyMHz8eycnJ2LhxI2pqarB48WIAQE5ODk6fPo2tW7e6jqmsrARwYULgP//5T1RWVkKv12PUqFEAgJUrVyIpKQnDhg2D1WrF6tWrUVlZibVr1/qUGwsBIiLSFAUS5gj42D4jIwNnzpzBCy+8gLq6OowZMwbFxcWIi4sDcOEGQpfeUyAhIcH1/+Xl5XjnnXcQFxeHr7/+GgBw9uxZPPLII6ivr4fZbEZCQgLKyspw8803+/ZahBAhcaNEq9UKs9mMqddkIlzR+xVLDPKtW8WbtgFXyYlztbyazRkuZ7EMXZu8t074dw45cc61SYkDAGHNnS8DuhxKy3kpcQBAnPftjmEe4/h45zGvsex2OYEcct4DACBkxQqNfx6vCO2iDaXYhaampi7H3C9Xx+dE3Ev/BZ2PN925lLO1FaeeWRbQfHsK5wgQERFpGIcGiIhIW7jokAoLASIi0hYWAiocGiAiItIw9ggQEZGmyLyzYCjwuUegrKwMM2fORExMDBRFwc6dO1XPCyGwYsUKxMTEoE+fPrjtttvw6aefdhm3qKgIo0aNgsFgwKhRo7Bjxw5fUyMiIupaL6w+GMx8LgTOnTuHcePG4Y033nD7/G9+8xusWrUKb7zxBj788ENYLBZMmzYNzc3NHmMePXoUGRkZyMzMxMcff4zMzEzMmzcPx44d8zU9IiIi8oHPQwPp6elIT093+5wQAnl5eVi2bBlmz54NAHjrrbcQFRWFd955Bz/72c/cHpeXl4dp06YhJycHwIU7LB08eBB5eXkoKCjwNUUiIiLPOFlQRepkwerqatTX16vWXDYYDLj11lu9rrl89OjRTus0T58+3esxNput0zrQREREXemYI+DvFiqkFgIdKyv5uuZyfX29z8fk5uaq1oCOjY31I3MiIiJtCsjlg5ez5rKvx+Tk5KCpqcm11dbWXn7CRESkHUKRs4UIqZcPWiwWABe+4V+8DGJXay5bLBaf12k2GAwwGAx+ZkxERJrDOQIqUnsE4uPjYbFYVGsu2+12HDx40Ouay8nJyZ3Wad63b5/UdZqJiIgAzhG4lM89Ai0tLfjyyy9dj6urq1FZWYl+/fph8ODBWLp0KV588UUMGzYMw4YNw4svvoirrroKDzzwgOuYrKwsXHfddcjNzQUAPP7445g8eTJefvllzJo1C7t27cL+/ftx+PBhCS+RiIiIPPG5EDh+/DimTJniepydnQ0AmD9/PvLz8/HUU0/hu+++w6OPPop//etfmDhxIvbt24fIyEjXMTU1NdDp/t0ZkZKSgm3btuG5557D8uXLMXToUBQWFmLixIn+vDYiIqLOODSgoggRGgtud6wzPfWaTIQrer9iiUHRXTfqprYBV8mJc7W86RzOcDmTXHRt8t464d/JWT8+/FyblDgAENZskxJHaTkvJQ4AiPOtcuK0yokDAMJulxPIIec9AABCVqzQ+OfxitAu2lCKXWhqaoLJZArIz+j4nLh++YsIMxr9iuVobcVXv342oPn2FC46REREpGEht+iQsLdBdHGpYld0TS2SsgH0kr5RhDdHSIkDAPDz/Lg4nXLiAFBscr7BKTZJ304BKN/J6RGQ+u1b1utrk9dzIuubvHDy2zf1EA4NqIRcIUBEROQVCwEVDg0QERFpGHsEiIhIU2TcByCU7iPAHgEiIiINYyFARESkYRwaICIibeFkQRUWAkREpCmcI6DGQoCIiLQnhD7I/cU5AkRERBrGHgEiItIWzhFQYSFARESawjkCahwaICIi0jD2CBARkbZwaECFhQAREWkKhwbUODRARESkYewRICIibeHQgAoLASIi0hYWAiocGiAiItKwkOsREG0OCKXdvxjnz0vKBlCcTilxws5J/FUpipw4QmJJ3O6QEka0tUmJAwDCLieWsNulxAEA0e7fe9vFIed8A4BwSnofCDl/KxdihdDXNZKOkwXVQq4QICIi8opDAyosBIiISFtYCKhwjgAREZGGsUeAiIg0hXME1FgIEBGRtnBoQIVDA0RERBrGHgEiItIUDg2osRAgIiJt4dCACocGiIiINIw9AkREpC3sEVBhIUBERJqifL/5GyNUcGiAiIhIw9gjQERE2sKhARUWAkREpCm8fFDN56GBsrIyzJw5EzExMVAUBTt37nQ919bWhqeffhpjx45F3759ERMTg6ysLHz77bdeY+bn50NRlE5ba2urzy+IiIjIKyFpCxE+FwLnzp3DuHHj8MYbb3R67vz58/joo4+wfPlyfPTRR9i+fTu++OIL3HPPPV3GNZlMqKurU21Go9HX9IiIiMgHPhcC6enp+M///E/Mnj2703NmsxklJSWYN28ehg8fjqSkJKxZswbl5eWoqanxGldRFFgsFtVGREQUEL3QG7Bu3TrEx8fDaDQiMTERhw4d8ti2rq4ODzzwAIYPHw6dToelS5e6bVdUVIRRo0bBYDBg1KhR2LFjh895BXyOQFNTExRFwTXXXOO1XUtLC+Li4uBwOHDjjTfi17/+NRISEjy2t9lssNlsrsdWqxUAIBwOCMW/iyHEdxKHJBwOOXHCwuTEAQBF0oUvQmLfmKTzJBxOKXEAAG1tUsJIzUnIiSWcEn93knKS+n4i8qI35ggUFhZi6dKlWLduHSZNmoQ333wT6enp+OyzzzB48OBO7W02GwYMGIBly5bhtddecxvz6NGjyMjIwK9//Wvce++92LFjB+bNm4fDhw9j4sSJ3c4toJcPtra24plnnsEDDzwAk8nksd2IESOQn5+P3bt3o6CgAEajEZMmTUJVVZXHY3Jzc2E2m11bbGxsIF4CERGR31atWoWFCxdi0aJFGDlyJPLy8hAbG4v169e7bT9kyBC8/vrryMrKgtlsdtsmLy8P06ZNQ05ODkaMGIGcnBxMnToVeXl5PuUWsEKgra0N9913H5xOJ9atW+e1bVJSEh588EGMGzcOqampePfdd3HDDTdgzZo1Ho/JyclBU1OTa6utrZX9EoiIKBRJnCxotVpV28U91R3sdjvKy8uRlpam2p+WloYjR45c9ss4evRop5jTp0/3OWZACoG2tjbMmzcP1dXVKCkp8dob4DYpnQ4TJkzw2iNgMBhgMplUGxERUVc6hgb83QAgNjZW1Tudm5vb6ec1NjbC4XAgKipKtT8qKgr19fWX/Trq6+ulxJQ+R6CjCKiqqsKBAwdw7bXX+hxDCIHKykqMHTtWdnpERETS1NbWqr6IGgwGj22VS+ZnCSE67fOVjJg+FwItLS348ssvXY+rq6tRWVmJfv36ISYmBnPnzsVHH32E//3f/4XD4XBVJv369YNerwcAZGVl4brrrnNVTitXrkRSUhKGDRsGq9WK1atXo7KyEmvXrvU1PSIiIu8k3lmwOz3S/fv3R1hYWKdv6g0NDZ2+0fvCYrFIienz0MDx48eRkJDgmtGfnZ2NhIQE/OpXv8I333yD3bt345tvvsGNN96I6Oho13bxmEVNTQ3q6upcj8+ePYtHHnkEI0eORFpaGk6fPo2ysjLcfPPNvqZHRETklcyhge7Q6/VITExESUmJan9JSQlSUlIu+3UkJyd3irlv3z6fY/rcI3DbbbdBeLnMx9tzHUpLS1WPX3vtNY+XRxAREV3psrOzkZmZifHjxyM5ORkbN25ETU0NFi9eDODCBPjTp09j69atrmMqKysBXOiJ/+c//4nKykro9XqMGjUKAPD4449j8uTJePnllzFr1izs2rUL+/fvx+HDh33KjWsNEBGRtvTCokMZGRk4c+YMXnjhBdTV1WHMmDEoLi5GXFwcgAs3ELr0xnsX30unvLwc77zzDuLi4vD1118DAFJSUrBt2zY899xzWL58OYYOHYrCwkKf7iEAAIrozlf4K4DVaoXZbMZtutkIVyL8iqUzep7s4StF718uLryhULfwhkLdDcMbClFwaRdtKMUuNDU1BewqsI7PiR8ueBFhev9uYe+wt+Kv+c8GNN+ewh4BIiLSFK4+qBbQOwsSERFRcGOPABERaUsvzBEIZiwEiIhIUxQhoPg5J8Xf44MJhwaIiIg0jD0CRESkLRwaUGEhQEREmsKrBtQ4NEBERKRh7BEgIiJt4dCASugVAsIJwL87nYm2djm5AIBT0l3XdCHeeSPpPEm9Uaasux0G4138ZAqh2dOkDRwaUAvxTxciIiLyJvR6BIiIiLzh0IAKCwEiItIUDg2osRAgIiJtYY+ACucIEBERaRh7BIiISHNCqWvfXywEiIhIW4Tw/7LXELpslkMDREREGsYeASIi0hReNaDGQoCIiLSFVw2ocGiAiIhIw9gjQEREmqI4L2z+xggVLASIiEhbODSgwqEBIiIiDWOPABERaQqvGlBjIUBERNrCGwqpsBAgIiJNYY+AWugVAsL/WSDC4ZCTCwCIEJpaegUQTol/ncH4uwuhbyFEFBxCrxAgIiLyhlcNqLAQICIiTeHQgBovHyQiItIw9ggQEZG28KoBFRYCRESkKRwaUPN5aKCsrAwzZ85ETEwMFEXBzp07Vc8vWLAAiqKotqSkpC7jFhUVYdSoUTAYDBg1ahR27Njha2pERETkI58LgXPnzmHcuHF44403PLa58847UVdX59qKi4u9xjx69CgyMjKQmZmJjz/+GJmZmZg3bx6OHTvma3pERETeCUlbiPB5aCA9PR3p6ele2xgMBlgslm7HzMvLw7Rp05CTkwMAyMnJwcGDB5GXl4eCggJfUyQiIvKIQwNqAblqoLS0FAMHDsQNN9yAhx9+GA0NDV7bHz16FGlpaap906dPx5EjRzweY7PZYLVaVRsRERH5RnohkJ6ejj/+8Y94//338eqrr+LDDz/E7bffDpvN5vGY+vp6REVFqfZFRUWhvr7e4zG5ubkwm82uLTY2VtprICKiEOYUcrYQIf2qgYyMDNf/jxkzBuPHj0dcXBz27NmD2bNnezxOURTVYyFEp30Xy8nJQXZ2tuux1WplMUBERF3jnQVVAn75YHR0NOLi4lBVVeWxjcVi6fTtv6GhoVMvwcUMBgMMBoO0PImISBsUSJgjICWT4BDwOwueOXMGtbW1iI6O9tgmOTkZJSUlqn379u1DSkpKoNMjIiLSNJ97BFpaWvDll1+6HldXV6OyshL9+vVDv379sGLFCsyZMwfR0dH4+uuv8eyzz6J///649957XcdkZWXhuuuuQ25uLgDg8ccfx+TJk/Hyyy9j1qxZ2LVrF/bv34/Dhw9LeIlEREQX4Z0FVXwuBI4fP44pU6a4HneM08+fPx/r16/HJ598gq1bt+Ls2bOIjo7GlClTUFhYiMjISNcxNTU10On+3RmRkpKCbdu24bnnnsPy5csxdOhQFBYWYuLEif68NiIiok54+aCaz4XAbbfdBuGlEtq7d2+XMUpLSzvtmzt3LubOnetrOkREROQHrjVARETawqsGVFgIEBGRpihCQPFzjN/f44MJCwF3hFNeKIe0UNTTQugPnYjIExYCRESkLc7vN39jhAgWAkREpCkcGlAL+A2FiIiIKHixR4CIiLSFVw2osBAgIiJt4Z0FVTg0QEREmtJxZ0F/N1+tW7cO8fHxMBqNSExMxKFDh7y2P3jwIBITE2E0GnH99ddjw4YNqufz8/OhKEqnrbW11ae8WAgQEREFWGFhIZYuXYply5ahoqICqampSE9PR01Njdv21dXVuOuuu5CamoqKigo8++yz+MUvfoGioiJVO5PJhLq6OtVmNBp9yo1DA0REpC29MDSwatUqLFy4EIsWLQIA5OXlYe/evVi/fr1rAb6LbdiwAYMHD0ZeXh4AYOTIkTh+/DheeeUVzJkzx9VOURRYLJbLfx1gjwAREWmM4pSzdZfdbkd5eTnS0tJU+9PS0nDkyBG3xxw9erRT++nTp+P48eNoa2tz7WtpaUFcXBwGDRqEGTNmoKKiovuJfY+FABER0WWyWq2qzWazdWrT2NgIh8OBqKgo1f6oqCjU19e7jVtfX++2fXt7OxobGwEAI0aMQH5+Pnbv3o2CggIYjUZMmjQJVVVVPr0GFgJERKQtHUMD/m4AYmNjYTabXZu7bv4OiqJckobotK+r9hfvT0pKwoMPPohx48YhNTUV7777Lm644QasWbPGp9PBOQJERKQtEu8jUFtbC5PJ5NptMBg6Ne3fvz/CwsI6fftvaGjo9K2/g8Vicds+PDwc1157rdtjdDodJkyYwB4BIiKinmIymVSbu0JAr9cjMTERJSUlqv0lJSVISUlxGzc5OblT+3379mH8+PGIiIhwe4wQApWVlYiOjvbpNbAQICIiTelYa8DfzRfZ2dn4/e9/j82bN+PEiRN44oknUFNTg8WLFwMAcnJykJWV5Wq/ePFinDp1CtnZ2Thx4gQ2b96MTZs24cknn3S1WblyJfbu3YuvvvoKlZWVWLhwISorK10xu4tDA0REpC29cPlgRkYGzpw5gxdeeAF1dXUYM2YMiouLERcXBwCoq6tT3VMgPj4excXFeOKJJ7B27VrExMRg9erVqksHz549i0ceeQT19fUwm81ISEhAWVkZbr75Zp9yU4Tw92wEB6vVCrPZjNswC+GK+26TbvMyeYM0JDT+NIiuCO2iDaXYhaamJtWYu0wdnxNTEnMQHu7bTXcu1d7eigPluQHNt6ewR4CIiLRFAPDhPgAeY4QIFgLu8JsgEVHIupwxfncxQgULASIi0hYBCXMEpGQSFHjVABERkYaxR4CIiLSlF64aCGYsBIiISFucAPy9OMzfyYZBhEMDREREGsYeASIi0hReNaDGQoCIiLSFcwRUODRARESkYewRICIibWGPgAoLASIi0hYWAiocGiAiItIw9ggQEZG28D4CKiwEiIhIU3j5oBoLASIi0hbOEVDxeY5AWVkZZs6ciZiYGCiKgp07d6qeVxTF7fbb3/7WY8z8/Hy3x7S2tvr8goiIiKj7fO4ROHfuHMaNG4ef/vSnmDNnTqfn6+rqVI//9Kc/YeHChW7bXsxkMuHkyZOqfUaj0df0iIiIvHMKQPHzG70zdHoEfC4E0tPTkZ6e7vF5i8Wierxr1y5MmTIF119/vde4iqJ0OpaIiEg6Dg2oBPTywX/84x/Ys2cPFi5c2GXblpYWxMXFYdCgQZgxYwYqKiq8trfZbLBaraqNiIiIfBPQQuCtt95CZGQkZs+e7bXdiBEjkJ+fj927d6OgoABGoxGTJk1CVVWVx2Nyc3NhNptdW2xsrOz0iYgoJIl/9wpc7gb2CHTL5s2b8ZOf/KTLsf6kpCQ8+OCDGDduHFJTU/Huu+/ihhtuwJo1azwek5OTg6amJtdWW1srO30iIgpF/hYBMoYWgkjALh88dOgQTp48icLCQp+P1el0mDBhgtceAYPBAIPB4E+KREREmhewHoFNmzYhMTER48aN8/lYIQQqKysRHR0dgMyIiEjTnELOFiJ87hFoaWnBl19+6XpcXV2NyspK9OvXD4MHDwYAWK1WvPfee3j11VfdxsjKysJ1112H3NxcAMDKlSuRlJSEYcOGwWq1YvXq1aisrMTatWsv5zURERF5JpwXNn9jhAifC4Hjx49jypQprsfZ2dkAgPnz5yM/Px8AsG3bNgghcP/997uNUVNTA53u350RZ8+exSOPPIL6+nqYzWYkJCSgrKwMN998s6/pERERkQ8UIUJjxoPVaoXZbMZtmIVwJaK30yEiIh+0izaUYheamppgMpkC8jM6PifuiP0PhOv8m2PW7rRhf+36gObbU7jWABERaYtTwuV/Wp4jQEREdEXjnQVVAnofASIiIgpu7BEgIiJtEZDQIyAlk6DAQoCIiLSFQwMqHBogIiLSMPYIEBGRtjidAPy8IZBTwzcUIiIiuqJxaECFQwNEREQaxh4BIiLSFvYIqLAQICIibeGdBVU4NEBERKRh7BEgIiJNEcIJ4ecywv4eH0xYCBARkbYI4X/XPucIEBERXaGEhDkCIVQIcI4AERGRhrFHgIiItMXpBBQ/x/g5R4CIiOgKxaEBFQ4NEBERaRh7BIiISFOE0wnh59AALx8kIiK6UnFoQIVDA0RERBrGHgEiItIWpwAU9gh0YCFARETaIgQAfy8fDJ1CgEMDREREGsYeASIi0hThFBB+Dg2IEOoRYCFARETaIpzwf2iAlw8SERFdkdgjoMY5AkRERBoWMj0CHdVZO9r8vk8EERH1rHa0AeiZb9rtwuZ3135HvqEgZAqB5uZmAMBhFPdyJkREdLmam5thNpsDEluv18NiseBwvZzPCYvFAr1eLyVWb1JEiAx0OJ1OfPvtt4iMjISiKG7bWK1WxMbGora2FiaTqYczvHzMu+ddqbkz757FvOURQqC5uRkxMTHQ6QI3at3a2gq73S4lll6vh9FolBKrN4VMj4BOp8OgQYO61dZkMgXNm98XzLvnXam5M++exbzlCFRPwMWMRmNIfHjLxMmCREREGsZCgIiISMM0VQgYDAY8//zzMBgMvZ2KT5h3z7tSc2fePYt5UygImcmCRERE5DtN9QgQERGRGgsBIiIiDWMhQEREpGEsBIiIiDQs5AqBdevWIT4+HkajEYmJiTh06JDX9gcPHkRiYiKMRiOuv/56bNiwoYcyvSA3NxcTJkxAZGQkBg4ciB/96Ec4efKk12NKS0uhKEqn7fPPP++hrIEVK1Z0+vkWi8XrMb19rjsMGTLE7flbsmSJ2/a9db7Lysowc+ZMxMTEQFEU7Ny5U/W8EAIrVqxATEwM+vTpg9tuuw2ffvppl3GLioowatQoGAwGjBo1Cjt27OixvNva2vD0009j7Nix6Nu3L2JiYpCVlYVvv/3Wa8z8/Hy3v4PW1tYeyRsAFixY0OnnJyUldRm3N883ALfnTVEU/Pa3v/UYsyfONwWPkCoECgsLsXTpUixbtgwVFRVITU1Feno6ampq3Lavrq7GXXfdhdTUVFRUVODZZ5/FL37xCxQVFfVYzgcPHsSSJUvwwQcfoKSkBO3t7UhLS8O5c+e6PPbkyZOoq6tzbcOGDeuBjP9t9OjRqp//ySefeGwbDOe6w4cffqjKu6SkBADw4x//2OtxPX2+z507h3HjxuGNN95w+/xvfvMbrFq1Cm+88QY+/PBDWCwWTJs2zbXuhjtHjx5FRkYGMjMz8fHHHyMzMxPz5s3DsWPHeiTv8+fP46OPPsLy5cvx0UcfYfv27fjiiy9wzz33dBnXZDKpzn9dXZ3UO8R1db4B4M4771T9/OJi7/es7+3zDaDTOdu8eTMURcGcOXO8xg30+aYgIkLIzTffLBYvXqzaN2LECPHMM8+4bf/UU0+JESNGqPb97Gc/E0lJSQHLsSsNDQ0CgDh48KDHNgcOHBAAxL/+9a+eS+wSzz//vBg3bly32wfjue7w+OOPi6FDhwqn0+n2+WA43wDEjh07XI+dTqewWCzipZdecu1rbW0VZrNZbNiwwWOcefPmiTvvvFO1b/r06eK+++6TnrMQnfN25//+7/8EAHHq1CmPbbZs2SLMZrPc5Lxwl/f8+fPFrFmzfIoTjOd71qxZ4vbbb/fapqfPN/WukOkRsNvtKC8vR1pammp/Wloajhw54vaYo0ePdmo/ffp0HD9+HG1tvbPEZFNTEwCgX79+XbZNSEhAdHQ0pk6digMHDgQ6tU6qqqoQExOD+Ph43Hffffjqq688tg3Gcw1ceN+8/fbbeOihhzwuVtWht8/3xaqrq1FfX686pwaDAbfeeqvH9zvg+ffg7ZhAa2pqgqIouOaaa7y2a2lpQVxcHAYNGoQZM2agoqKiZxK8SGlpKQYOHIgbbrgBDz/8MBoaGry2D7bz/Y9//AN79uzBwoULu2wbDOebekbIFAKNjY1wOByIiopS7Y+KikJ9fb3bY+rr6922b29vR2NjY8By9UQIgezsbNxyyy0YM2aMx3bR0dHYuHEjioqKsH37dgwfPhxTp05FWVlZj+U6ceJEbN26FXv37sXvfvc71NfXIyUlBWfOnHHbPtjOdYedO3fi7NmzWLBggcc2wXC+L9Xxnvbl/d5xnK/HBFJrayueeeYZPPDAA14XvxkxYgTy8/Oxe/duFBQUwGg0YtKkSaiqquqxXNPT0/HHP/4R77//Pl599VV8+OGHuP3222Gz2TweE2zn+6233kJkZCRmz57ttV0wnG/qOSGz+mCHS7/VCSG8ftNz197d/p7w2GOP4a9//SsOHz7std3w4cMxfPhw1+Pk5GTU1tbilVdeweTJkwOdJoAL/yh2GDt2LJKTkzF06FC89dZbyM7OdntMMJ3rDps2bUJ6ejpiYmI8tgmG8+2Jr+/3yz0mENra2nDffffB6XRi3bp1XtsmJSWpJuZNmjQJN910E9asWYPVq1cHOlUAQEZGhuv/x4wZg/HjxyMuLg579uzx+sEaLOcbADZv3oyf/OQnXY71B8P5pp4TMj0C/fv3R1hYWKdKu6GhoVNF3sFisbhtHx4ejmuvvTZgubrz85//HLt378aBAwe6vZzyxZKSknq1Wu/bty/Gjh3rMYdgOtcdTp06hf3792PRokU+H9vb57vjCg1f3u8dx/l6TCC0tbVh3rx5qK6uRklJic9L4ep0OkyYMKFXfwfR0dGIi4vzmkOwnG8AOHToEE6ePHlZ7/dgON8UOCFTCOj1eiQmJrpmgHcoKSlBSkqK22OSk5M7td+3bx/Gjx+PiIiIgOV6MSEEHnvsMWzfvh3vv/8+4uPjLytORUUFoqOjJWfXfTabDSdOnPCYQzCc60tt2bIFAwcOxN133+3zsb19vuPj42GxWFTn1G634+DBgx7f74Dn34O3Y2TrKAKqqqqwf//+yyoEhRCorKzs1d/BmTNnUFtb6zWHYDjfHTZt2oTExESMGzfO52OD4XxTAPXWLMVA2LZtm4iIiBCbNm0Sn332mVi6dKno27ev+Prrr4UQQjzzzDMiMzPT1f6rr74SV111lXjiiSfEZ599JjZt2iQiIiLEf//3f/dYzv/xH/8hzGazKC0tFXV1da7t/PnzrjaX5v3aa6+JHTt2iC+++EL87W9/E88884wAIIqKinos71/+8peitLRUfPXVV+KDDz4QM2bMEJGRkUF9ri/mcDjE4MGDxdNPP93puWA5383NzaKiokJUVFQIAGLVqlWioqLCNbv+pZdeEmazWWzfvl188skn4v777xfR0dHCarW6YmRmZqqumvnLX/4iwsLCxEsvvSROnDghXnrpJREeHi4++OCDHsm7ra1N3HPPPWLQoEGisrJS9Z632Wwe816xYoX485//LP7+97+LiooK8dOf/lSEh4eLY8eO9Ujezc3N4pe//KU4cuSIqK6uFgcOHBDJycniuuuuC+rz3aGpqUlcddVVYv369W5j9Mb5puARUoWAEEKsXbtWxMXFCb1eL2666SbVZXjz588Xt956q6p9aWmpSEhIEHq9XgwZMsTjH0qgAHC7bdmyxWPeL7/8shg6dKgwGo3iBz/4gbjlllvEnj17ejTvjIwMER0dLSIiIkRMTIyYPXu2+PTTTz3mLETvn+uL7d27VwAQJ0+e7PRcsJzvjssWL93mz58vhLhwCeHzzz8vLBaLMBgMYvLkyeKTTz5Rxbj11ltd7Tu89957Yvjw4SIiIkKMGDFCekHjLe/q6mqP7/kDBw54zHvp0qVi8ODBQq/XiwEDBoi0tDRx5MiRHsv7/PnzIi0tTQwYMEBERESIwYMHi/nz54uamhpVjGA73x3efPNN0adPH3H27Fm3MXrjfFPw4DLEREREGhYycwSIiIjIdywEiIiINIyFABERkYaxECAiItIwFgJEREQaxkKAiIhIw1gIEBERaRgLASIiIg1jIUBERKRhLASIiIg0jIUAERGRhrEQICIi0rD/B9qt2KsVy0FaAAAAAElFTkSuQmCC",
|
|
"text/plain": [
|
|
"<Figure size 640x480 with 2 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"plt.imshow(C_t)\n",
|
|
"plt.colorbar()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"attachments": {},
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Widget"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 127,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "77835379be7343d0aff52d275a86be97",
|
|
"version_major": 2,
|
|
"version_minor": 0
|
|
},
|
|
"text/plain": [
|
|
"interactive(children=(IntSlider(value=0, description='w', max=999), Output()), _dom_classes=('widget-interact'…"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"<function __main__.update(w=1)>"
|
|
]
|
|
},
|
|
"execution_count": 127,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"def update(w = 1):\n",
|
|
" fig = plt.figure(figsize = (10,7))\n",
|
|
" y = records[w]\n",
|
|
" plt.imshow(y, vmin=0, vmax=50)\n",
|
|
" \n",
|
|
"interact(update, w = IntSlider(min=0, max = 999, step = 1, value = 0))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 35,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD4CAYAAADl7fPiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAREUlEQVR4nO3dfYxc1X3G8e/jtaGVMbGJwQHjEppaSFZU3MhyGpVWpiTURihOqjS1VbVOS7UkClIjNapoK0GU/kNVUaTECOIkFqRKgL45sRQLsGglgpQXFmRenEBxLUfs4thNTHFwiGDtX//Yu9Gc9cz63rn37NyZPh9pNXfuPXPuuTvLw7z8fI4iAjOzWYsGPQAzaxeHgpklHApmlnAomFnCoWBmicWDHkA3YxcujSUXLx/0MErL8wWOKgygQre5vmyq1G+ma6siR7+Zxlrht1XaW6+e4PSpU127bmUoLLl4OWvu+Hjj/Vb6jzfKPxVnSraNM+X7rNT2dIU/m1xtK/y+VGkM5ZtW6VdnSjas8DyU7hMqBUiOUJj83F09j/ntg5klaoWCpM2SXpR0SNKtXY6fL+mh4vh3Jb2zzvnMLL++Q0HSGHA3sAVYB2yXtG5Os5uAVyPi14C7gL/v93xmtjDqvFLYCByKiMMR8SbwILB1TputwP3F9r8C10nK8RbJzBpSJxRWAy933J8s9nVtExHTwGvA27t1Jmlc0oSkidMnT9UYlpnV0ZoPGiNiV0RsiIgNYxcuHfRwzP7fqhMKU8CajvuXF/u6tpG0GHgb8JMa5zSzzOqEwpPAWklXSjoP2AbsndNmL7Cj2P4I8B/hf6tt1mp9Fy9FxLSkW4BHgDFgd0QclPRZYCIi9gJfBv5J0iHgBDPBYWYtVquiMSL2Afvm7LutY/vnwB/UOce5x1ClcYUqwQptS/fbhnLkXEb52oZJA7/b1nzQaGbt4FAws4RDwcwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLNEKyduhQrly5lKlwdeiZujzDpr2/JNK8k1hpJtlal8u9JMQwv8x+hXCmaWcCiYWcKhYGYJh4KZJRwKZpZwKJhZwqFgZok6K0StkfSfkr4v6aCkv+jSZpOk1yQdKH5u69aXmbVHneKlaeAvI+JpScuApyTtj4jvz2n3rYi4scZ5zGwB9f1KISKORsTTxfZPgR9w9gpRZjZkGilzLlaT/g3gu10Ov0/SM8ArwKcj4mCPPsaBcYDFK99GnGn+445K1aJVSqLPlG1XpSS7dFMoef6Zjiu0rUBVypErjLdKmXG1kuTmZ+DOVrqc63fQQ+3/8iRdAPwb8KmIODnn8NPAFRFxNfB54Ou9+ulcNm7RMi8bZzYotUJB0hJmAuGrEfHvc49HxMmIeL3Y3gcskbSyzjnNLK863z6ImRWgfhAR/9ijzTtml56XtLE4n9eSNGuxOp8p/Bbwx8Bzkg4U+/4G+BWAiLiXmfUjPyFpGngD2Oa1JM3arc5akk9wjs9WImInsLPfc5jZwnNFo5klHApmlnAomFnCoWBmCYeCmSXaO5tz6YaZSocrtC1dvlypFDjTrMtV+s1VPp2pJLpSOXDJfrPN5jzo0vR5+vQrBTNLOBTMLOFQMLOEQ8HMEg4FM0s4FMws4VAws4RDwcwSDgUzS7S3ovF0pWkwy/VZpZIuR0VjpWrCNlT9ZZqMNVvbDM9vG6oUh23iVjMbLQ4FM0s0McX7EUnPFcvCTXQ5Lkmfk3RI0rOS3lP3nGaWT1OfKVwbET/ucWwLsLb4eS9wT3FrZi20EG8ftgJfiRnfAZZLunQBzmtmfWgiFAJ4VNJTxdJvc60GXu64P0mXNScljUuakDRx+uSpBoZlZv1o4u3DNRExJekSYL+kFyLi8aqdRMQuYBfA+b+62mtDmA1I7VcKETFV3B4H9gAb5zSZAtZ03L+82GdmLVR3LcmlkpbNbgPXA8/PabYX+JPiW4jfBF6LiKN1zmtm+dR9+7AK2FMsF7kY+FpEPCzp4/CLpeP2ATcAh4CfAX9a85xmllGtUIiIw8DVXfbf27EdwCerdaw8k6FWGUKOktlcpcsV+s1WulyhLD1b6XKG8ulcpcuVrqsClzmbWeMcCmaWcCiYWcKhYGYJh4KZJRwKZpZwKJhZwqFgZgmHgpklHApmlmjxbM4Z8irTDLqly4xzlE5TsWT2dJ5+K5UuV5mpe8CzRA98NumqbRvo068UzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEn2HgqSriqXiZn9OSvrUnDabJL3W0ea22iM2s6z6Ll6KiBeB9QCSxpiZtn1Pl6bfiogb+z2PmS2spt4+XAf8d0T8sKH+zGxAmipz3gY80OPY+yQ9A7wCfDoiDnZrVCw5Nw4wdtFyqFIKW1au0tKyM0q3YhbjCqXLVUqih6h0uUrbbM9Dtr/FCm17aGIp+vOADwL/0uXw08AVEXE18Hng6736iYhdEbEhIjaMLVtad1hm1qcm3j5sAZ6OiGNzD0TEyYh4vdjeByyRtLKBc5pZJk2EwnZ6vHWQ9A4Vy0dJ2lic7ycNnNPMMqn1mUKxfuQHgJs79nUuGfcR4BOSpoE3gG3FilFm1lJ1l407Bbx9zr7OJeN2AjvrnMPMFpYrGs0s4VAws4RDwcwSDgUzSzgUzCzRztmcgzxlzlXHUFLpUtgKZbBVZn7OVgo8bKXLFcabYzbnVpQ5N3B+v1Iws4RDwcwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEu0sc0bVSmzLGvAMuvlmaK7SNk/5dKXxVpklOtfs0znKnKuUxg+4zHk+fqVgZolSoSBpt6Tjkp7v2HeRpP2SXipuV/R47I6izUuSdjQ1cDPLo+wrhfuAzXP23Qo8FhFrgceK+wlJFwG3A+8FNgK39woPM2uHUqEQEY8DJ+bs3grcX2zfD3yoy0N/D9gfESci4lVgP2eHi5m1SJ3PFFZFxNFi+0fAqi5tVgMvd9yfLPaZWUs18kFjsZZDrc9IJY1LmpA0cfr115sYlpn1oU4oHJN0KUBxe7xLmylgTcf9y4t9Z0nWkrzgghrDMrM66oTCXmD224QdwDe6tHkEuF7SiuIDxuuLfWbWUmW/knwA+DZwlaRJSTcBdwAfkPQS8P7iPpI2SPoSQEScAP4OeLL4+Wyxz8xaqlRFY0Rs73Houi5tJ4A/77i/G9jd1+jMbMG1s8w5gLIlq5Fn1uccZajVyoYzzeZcpRS3QtlwlTLnRblmic5R5pypHDlX+XQTfbrM2cwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEg4FM0u0tsxZ082XL2ebQbdkqXWl8+cqXc40Q/MwlS5XaZttNudMz28T/ErBzBIOBTNLOBTMLOFQMLOEQ8HMEg4FM0s4FMwscc5Q6LGO5D9IekHSs5L2SFre47FHJD0n6YCkiQbHbWaZlHmlcB9nL/W2H3h3RPw68F/AX8/z+GsjYn1EbOhviGa2kM4ZCt3WkYyIRyNiurj7HWYWeTGzEdBEmfOfAQ/1OBbAo5IC+EJE7OrViaRxYBxg8fIVLMpQ5lxJjtmcq5w/Q8lu9bajWbpcpW2lsbagjL303+I87WqFgqS/BaaBr/Zock1ETEm6BNgv6YXilcdZisDYBXD+mjULXO1tZrP6/vZB0seAG4E/KhaYPUtETBW3x4E9wMZ+z2dmC6OvUJC0Gfgr4IMR8bMebZZKWja7zcw6ks93a2tm7VHmK8lu60juBJYx85bggKR7i7aXSdpXPHQV8ISkZ4DvAd+MiIezXIWZNeacnyn0WEfyyz3avgLcUGwfBq6uNTozW3CuaDSzhEPBzBIOBTNLOBTMLOFQMLNEK2dzVlQrLx24BkpLz5JtNufypcuVnoM2lGVnKJ8edJk1kKXkfr52fqVgZgmHgpklHApmlnAomFnCoWBmCYeCmSUcCmaWcCiYWcKhYGaJVlY00oaKxgyzRFabgLNC5WGuSroWVPMNum2+SWbL/zFUes4a4FcKZpZwKJhZot9l4z4jaaqYn/GApBt6PHazpBclHZJ0a5MDN7M8+l02DuCuYjm49RGxb+5BSWPA3cAWYB2wXdK6OoM1s/z6WjaupI3AoYg4HBFvAg8CW/vox8wWUJ3PFG4pVp3eLWlFl+OrgZc77k8W+7qSNC5pQtLE6VOnagzLzOroNxTuAd4FrAeOAnfWHUhE7IqIDRGxYWzp0rrdmVmf+gqFiDgWEacj4gzwRbovBzcFrOm4f3mxz8xarN9l4y7tuPthui8H9ySwVtKVks4DtgF7+zmfmS2cc1Y0FsvGbQJWSpoEbgc2SVrPTN3fEeDmou1lwJci4oaImJZ0C/AIMAbsjoiDOS7CzJqTbdm44v4+4KyvK88pYNFbJct8y1cDD5dKk7FW6DdX2XCm8Q56ktdqfeYpXc4yyasnbjWzshwKZpZwKJhZwqFgZgmHgpklHApmlnAomFnCoWBmCYeCmSUcCmaWaOdszoCmyzas0mn5phUmU86i2szPmfqtIEspbs5+S7bNNetylVmiF2WYUXq+6/crBTNLOBTMLOFQMLOEQ8HMEg4FM0s4FMws4VAws0SZORp3AzcCxyPi3cW+h4CriibLgf+NiPVdHnsE+ClwGpiOiA2NjNrMsilTvHQfsBP4yuyOiPjD2W1JdwKvzfP4ayPix/0O0MwWVpmJWx+X9M5uxyQJ+Cjwuw2Py8wGpG6Z828DxyLipR7HA3hUUgBfiIhdvTqSNA6MAyy5cAVjb5UbwJmxCqOtUro8RJ+2DLoUuDX9ZpjVOtesy1VKlxdNl/8lLCr5zwPmew7qhsJ24IF5jl8TEVOSLgH2S3qhWLD2LEVg7AL45UvXZKrQN7Nz6fv/h5IWA78PPNSrTURMFbfHgT10X17OzFqkzovk9wMvRMRkt4OSlkpaNrsNXE/35eXMrEXOGQrFsnHfBq6SNCnppuLQNua8dZB0maTZFaFWAU9Iegb4HvDNiHi4uaGbWQ79LhtHRHysy75fLBsXEYeBq2uOz8wW2BB9xm5mC8GhYGYJh4KZJRwKZpZwKJhZopWzOesMjL1RsvEvle/3TIWrrVKyWnrm5ypl1m2o6cw1hjZcW1ktKMkuW7oMsPjnJTueZ5Zqv1Iws4RDwcwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEopoX82ppP8Bfjhn90pgFNePGNXrgtG9tlG4risi4uJuB1oZCt1ImhjFFaZG9bpgdK9tVK9rlt8+mFnCoWBmiWEKhZ6rSw25Ub0uGN1rG9XrAoboMwUzWxjD9ErBzBaAQ8HMEkMRCpI2S3pR0iFJtw56PE2RdETSc5IOSJoY9HjqkLRb0nFJz3fsu0jSfkkvFbcrBjnGfvS4rs9ImiqetwOSbhjkGJvW+lCQNAbcDWwB1gHbJa0b7KgadW1ErB+B773vAzbP2Xcr8FhErAUeK+4Pm/s4+7oA7iqet/URsa/L8aHV+lBgZqXqQxFxOCLeBB4Etg54TDZHRDwOnJizeytwf7F9P/ChhRxTE3pc10gbhlBYDbzccX+y2DcKAnhU0lOSxgc9mAxWRcTRYvtHzCw6PCpukfRs8fZi6N4WzWcYQmGUXRMR72HmrdEnJf3OoAeUS8x89z0q33/fA7wLWA8cBe4c6GgaNgyhMAWs6bh/ebFv6EXEVHF7HNjDzFulUXJM0qUAxe3xAY+nERFxLCJOR8QZ4IuM2PM2DKHwJLBW0pWSzgO2AXsHPKbaJC2VtGx2G7geeH7+Rw2dvcCOYnsH8I0BjqUxs0FX+DAj9ry1coWoThExLekW4BFgDNgdEQcHPKwmrAL2SIKZ5+FrEfHwYIfUP0kPAJuAlZImgduBO4B/lnQTM/8U/qODG2F/elzXJknrmXk7dAS4eVDjy8FlzmaWGIa3D2a2gBwKZpZwKJhZwqFgZgmHgpklHApmlnAomFni/wCOWVT0CjszYQAAAABJRU5ErkJggg==",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"import matplotlib.animation as animation\n",
|
|
"\n",
|
|
"\n",
|
|
"frames = []\n",
|
|
"fig = plt.figure()\n",
|
|
"for i in records:\n",
|
|
" frames.append([plt.imshow(i, vmin=0, vmax=10, animated=True)])\n",
|
|
"\n",
|
|
"ani = animation.ArtistAnimation(fig, frames, interval=50, blit=True, repeat_delay=1000)\n",
|
|
"\n",
|
|
"ani.save(\"/Users/hannessigner/Desktop/test.gif\")"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.11.4"
|
|
},
|
|
"orig_nbformat": 4
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|