diff --git a/Tests/sgfrd_annihilation_SCK.ipynb b/Tests/sgfrd_annihilation_SCK.ipynb new file mode 100644 index 0000000..402cbc6 --- /dev/null +++ b/Tests/sgfrd_annihilation_SCK.ipynb @@ -0,0 +1,232 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from ecell4 import *\n", + "from ecell4_base.core import *\n", + "from ecell4_base import sgfrd" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "SCK theory (long time assumption)\n", + "done by W.X. Chew et al." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# SCK theory (relative error becomes less than 1% when 100R^2/D <= t)\n", + "def calc_k(ka, D, R, t):\n", + " if t == 0:\n", + " return ka\n", + "\n", + " kappa = ka / D\n", + " gamma = 0.5772156\n", + " # Cc = 4 * D / (R**2) * np.exp(4 * np.pi / kappa - 2 * gamma)\n", + " # to avoid numerical error\n", + " lnCct = 4 * np.pi / kappa - 2 * gamma + np.log(4 * D / (R**2)) + np.log(t)\n", + " return 4 * np.pi * D * (1.0 / lnCct - gamma / (lnCct**2) - 1.311 / (lnCct**3) + 0.25 / (lnCct**4))\n", + "\n", + "def calc_S(ka, D, R, concB, delta_t, total_t):\n", + " steps = int(total_t / delta_t) + 1\n", + " ks = [calc_k(ka, D, R, (tau+0.5) * delta_t) for tau in range(0, steps+1)]\n", + " return ([tau * delta_t for tau in range(0, steps)],\n", + " [np.exp(-concB * sum(k * delta_t for k in ks[0:tau])) for tau in range(0, steps)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "parameter setting" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# system size and number of particles\n", + "L = np.sqrt(10.0)\n", + "el = Real3(L, L, L)\n", + "NA = 100\n", + "NB = 100\n", + "concA = NA / (L * L)\n", + "concB = NA / (L * L)\n", + "\n", + "# physical parameters of species\n", + "rA = 0.005\n", + "rB = 0.005\n", + "DA = 1.0\n", + "DB = 0.0\n", + "R = rA + rB\n", + "D = DA + DB\n", + "\n", + "ka = np.pi * 4.0\n", + "\n", + "duration = 0.2\n", + "steps = 50\n", + "ts = np.linspace(0, duration, steps+1)\n", + "delta_t = duration / steps\n", + "\n", + "N = 30 # a number of samples" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# initial condition\n", + "y0 = {'A': NA, 'B': NB}" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# model\n", + "with species_attributes():\n", + " A | {'radius': rA, 'D': DA}\n", + " B | {'radius': rB, 'D': DB}\n", + "\n", + "with reaction_rules():\n", + " A + B > B | ka\n", + "\n", + "m = get_model()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "rng = GSLRandomNumberGenerator()\n", + "rng.seed(0)\n", + "\n", + "# obs = run_simulation(np.linspace(0, duration, steps+1), y0, volume=el, model=m, return_type='observer', solver='sgfrd')\n", + "obs = ensemble_simulations(ts, y0, volume=el, model=m, return_type='observer', solver='sgfrd', n=N)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXzU1b3/8dchBAKyBRKULYIICpKQSAQV0aAoaLWsgguKQqto8dp7r9utFlMqbW3rrW31XoqtYrUiKohxu9YFlx8FJYEQdglbCAKyhAjIFji/P85MmIRJmCQzmSXv5+ORR2a+3+985/Cd8MnJ53vO5xhrLSIiEv0ahbsBIiISHAroIiIxQgFdRCRGKKCLiMQIBXQRkRjROFxvnJSUZLt27RqutxcRiUp5eXm7rbXJ/vaFLaB37dqV3NzccL29iEhUMsZsqWqfUi4iIjFCAV1EJEYooIuIxIiw5dBFJHDHjh2juLiYw4cPh7spUk8SEhLo3Lkz8fHxAb9GAV0kChQXF9OyZUu6du2KMSbczZEQs9ayZ88eiouL6datW8CvU8pFJAocPnyYdu3aKZg3EMYY2rVrV+O/yBTQRaKEgnnDUpvPWwFdRCRGnDagG2OeN8Z8a4xZWcV+Y4z5kzGm0BhTYIy5MPjNPFVBAWRnw8SJ7ntBQX28q0jDNX36dC644ALS0tJIT0/nyy+/BNwN20ceeYQePXpw4YUXcskll/D+++8DbgLh7t27AcjLy6Nbt24sW7aswnnz8/N57733yp9nZ2fz+9//vp7+VbElkB76LGBYNfuvBXp4vu4C/rfuzapeQQH8/vdQUgKdO7vvv/+9grpIqCxatIh33nmHpUuXUlBQwEcffUSXLl0A+PnPf8727dtZuXIlS5cuZf78+ezfv7/C6wsKChgzZgxz5swhIyOjwr7KAT1Ujh8/HvL3CLfTBnRr7efA3moOGQ783TqLgTbGmA7BaqA/8+ZBYiLs2gX/93/ucWKi2y4iwbd9+3aSkpJo2rQpAElJSXTs2JHvv/+e5557jj//+c/l+84880zGjh1b/to1a9YwYsQIXnrpJfr371/hvEePHmXq1KnMmTOH9PR05syZA8Dq1avJysrinHPO4U9/+lP58S+//DL9+/cnPT2du+++uzxIz549m9TUVPr06cPDDz9cfnyLFi34z//8T/r27cv06dMZMWJE+b4PP/yQkSNHBvlKhVcwhi12Arb6PC/2bNte+UBjzF24XjwpKSm1fsOiItczX7MGliyBSy6B1q3ddpFY99OfQn5+cM+Zng5PP131/muuuYZp06bRs2dPhgwZwrhx47jiiisoLCwkJSWFVq1aVfna4cOH8/LLL3PZZZedsq9JkyZMmzaN3NxcnnnmGcClXNauXcuCBQvYv38/5513Hvfccw+FhYXMmTOHhQsXEh8fz7333ss//vEPhgwZwsMPP0xeXh6JiYlcc801zJ8/nxEjRnDw4EEGDBjAU089hbWWXr16sWvXLpKTk3nhhReYOHFina9dJKnXm6LW2pnW2kxrbWZyst9iYQFJSYHSUujd2z1fvdo9r8PvCBGpRosWLcjLy2PmzJkkJyczbtw4Zs2aFdBrhwwZwl//+tcapTx+8IMf0LRpU5KSkmjfvj07d+7k448/Ji8vj4suuoj09HQ+/vhjNm7cyJIlS8jKyiI5OZnGjRtz66238vnnnwMQFxfH6NGjATdq5LbbbuPll19m3759LFq0iGuvvbbG1yKSBaOHvg3o4vO8s2dbyIwa5XLmiYnQoQOsWOG+T5oUyncViQzV9aRDKS4ujqysLLKyskhNTeXFF19k7NixFBUV8d1331XZS3/mmWeYPHky9957L3/5y18Cei9v+sb7vmVlZVhrmTBhAr/+9a8rHPvWW29VeZ6EhATi4uLKn995553ccMMNJCQkcOONN9K4cWzNrQxGDz0HuN0z2uVioNRae0q6JZjS0uCBB1xA79gRdu6EW29120Uk+NatW8f69evLn+fn53P22WfTvHlzJk2axP3338/Ro0cB2LVrF6+//nr5sY0aNeKVV15h7dq1TJ069ZRzt2zZ8pSbqP5cddVVvPHGG3z77bcA7N27ly1bttC/f38+++wzdu/ezfHjx5k9ezZXXHGF33N07NiRjh078sQTT3DnnXfW6BpEg0CGLc4GFgHnGWOKjTGTjDGTjTGTPYe8B2wECoHngHtD1lofaWluuOLs2e75mjX18a4iDdOBAweYMGECvXv3Ji0tjdWrV5OdnQ3AE088QXJyMr1796ZPnz5cf/31p/TWExISyMnJIScnh2effbbCvsGDB7N69eoKN0X96d27N0888QTXXHMNaWlpXH311Wzfvp0OHTrwm9/8hsGDB9O3b1/69evH8OHDqzzPrbfeSpcuXejVq1ftL0iEMtbasLxxZmamDdYCF337QqtW8MUXQTmdSMRZs2ZNTAagcJgyZQoZGRlMioIcrb/P3RiTZ63N9Hd8TMwUHTMGFi6E7SFN9IhItOvXrx8FBQWMHz8+3E0JiZgI6KNHg7Xw5pvhbomIRLK8vDw+//zzCjddY0lMBPTevaFXL3jjjXC3REQkfGIioIPrpX/2mZs9KiLSEMVMQB8zBk6cgPnzw90SEZHwiJmAnpYG3bsr7SIiDVfMBHRjXC/9k09gb3WlxESkVqoqn1tWVsbPfvYzevToQXp6Ounp6UyfPr38dXFxceXb09PT2bx5M59++imtW7cmPT2d888/nwceeKD8+FmzZpGcnExGRgY9evRg6NCh/Otf/wq4nbm5ufzbv/1bUP7Ns2bNYsqUKdUe8+mnn1Zo34wZM/j73/8elPevqZia9zp6NDz5JOTkwB13hLs1IrHDt3xu06ZN2b17d/nM0Mcee4wdO3awYsUKEhIS2L9/P0899VT5a5s1a0Z+pWpimzdvZtCgQbzzzjscOnSIjIwMRo4cycCBAwEYN25cebGuBQsWMGrUKBYsWBDQWPzMzEwyM/0O0w6JTz/9lBYtWnDppZcCMHny5NO8InRipocOkJnpCnTNnRvuloiEV7AXgAmkfG5CQgLgpvJ7Z5EGolmzZqSnp7Ntm/8SUIMHD+auu+5i5syZp+x7/fXX6dOnD3379uXyyy8HXIC9/vrrAVe5ccKECQwaNIizzz6befPm8dBDD5GamsqwYcM4duwYUHEhjtzcXLKysk55r7fffpsBAwaQkZHBkCFD2LlzJ5s3b2bGjBn84Q9/ID09nS+++KLCAh35+flcfPHFpKWlMXLkSEpKSgDIysri4Ycfpn///vTs2ZMvgjQrMqYCujGul/7Pf8J334W7NSLhEYoFYK655hq2bt1Kz549uffee/nss88AysvntmzZssrXHjp0qDzd4q/+eElJCevXry8PyP5ceOGFrF279pTt06ZN44MPPmD58uXk5OT4fe2GDRv45JNPyMnJYfz48QwePJgVK1bQrFkz3n333dP908tddtllLF68mGXLlnHTTTfx29/+lq5duzJ58mT+/d//nfz8fAYNGlThNbfffjtPPvkkBQUFpKam8otf/KJ8X1lZGV999RVPP/10he11EVMBHVwe/ehReOedcLdEJDy8C8AkJkKjRsFZACbQ8rkvvPAC6enpdOnSha1b3TIJ3pRLfn4+b/rM/vviiy/o27cvnTp1YujQoZx11llVvn9VJUoGDhzIHXfcwXPPPVdled5rr72W+Ph4UlNTOX78OMOGuQXYUlNT2bx5c4BXAIqLixk6dCipqan87ne/Y9WqVdUeX1payr59+8oLhU2YMKG8rC/AqFGjADd7tSbtqE7MBfSLL3YVGDXaRRqqoiK34IuvYCwA4y2f+4tf/IJnnnmGuXPncu6551JUVFReLfHOO+8kPz+f1q1bn7b++aBBg1i+fDmrVq3ib3/72yl5dl/Lli3zmz+fMWMGTzzxBFu3bqVfv37s2bPnlGO8aaJGjRoRHx+PMab8eVlZGQCNGzfmxIkTABw+fNhvG+677z6mTJnCihUr+Mtf/lLlcYHytstbHjgYYi6gr1zp/szMyYGHHtI6o9LweBeA8VXXBWBOVz53ypQp5QHu+PHj5TdMA9GtWzceeeQRnnzySb/7P/vsM2bOnMmPf/zjU/Zt2LCBAQMGMG3aNJKTk8v/Kqiprl27kpeXB8DcKm7ClZaW0qlTJwBefPHF8u1Vlf9t3bo1iYmJ5fnxl156qcqyvsESUwHdmzs891w4fhwWLdLi0dLwjBrl8uYlJW6ynfex5y/8WqmufO706dPp0KEDffr0ISMjg0GDBjFhwgQ6duwY8PknT57M559/Xp568K4x2rNnT371q18xd+5cvz30Bx98sHwt0UsvvZS+ffvW6t/3+OOPc//995OZmVlhQQxf2dnZ3HjjjfTr14+kpKTy7TfccANvvvlm+U1RXy+++CIPPvggaWlp5Ofn+60HH0wxUT7XKzvb/eAmJsKsWa5XMn48tGvn9olEq5qWzy0ocDnzoiLXMx81SgvARKOals+NqXHo3sWjAS66yOXRd++GgwfD2y6R+paWpgDeEMVUysU3d3j++dCypUu7aPFoEWkIYiqg++YOjXFldbdsgX79wt0ykboLV3pUwqM2n3dMBXTfxaOLi2HgQGjcGBYsCHfLROomISGBPXv2KKg3ENZa9uzZUz77NlAxlUOHU3OH33wDzz8Pv/wlnHFG+NolUhedO3emuLiYXSr432AkJCTQ2XtTMEAxF9Ar+8lP4NVX4ZVXwM8wVpGoEB8fT7du3cLdDIlwMZVy8WfgQOjbF555xq07KiISq2I+oBsDU6a4cbn/7/+FuzUiIqET8wEd4JZboE0bePbZcLdERCR0GkRAb97c1YWeO9fdJBURiUUNIqAD3HMPlJWBnxr5IiIxocEE9HPPhWuvhRkzoI5VL0VEIlKDCegFBZCcDDt3uhmlqsAoIrGmQQR0b1nd1q1d8a4vvnCLSSuoi0gsaRAB3bskV9u2cMUVcOAA7NhRtyW5REQiTYMI6L5LcnXv7paoW7YMgrSMn4hIRGgQAd23rK4xcPnlriJj5WW6RESiWUAB3RgzzBizzhhTaIx5xM/+FGPMAmPMMmNMgTHmuuA3tfYqL8mVnOxSMMuWuaXqRERiwWkDujEmDngWuBboDdxsjOld6bDHgNestRnATcD/BLuhdVG5rG7btvCzn7la6a+/Hu7WiYgERyDVFvsDhdbajQDGmFeB4cBqn2Ms0MrzuDUQcfMxK5fVPXHCldWdPh3GjoVGDSL5JCKxLJAw1gnY6vO82LPNVzYw3hhTDLwH3OfvRMaYu4wxucaY3HDXdW7UCB59FFauhLfeCmtTRESCIlj90puBWdbazsB1wEvGmFPOba2daa3NtNZmJicnB+mta2/cODeD9Je/VGldEYl+gQT0bUAXn+edPdt8TQJeA7DWLgISgKRgNDCUGjeG225zN0evuQayszXZSESiVyABfQnQwxjTzRjTBHfTM6fSMUXAVQDGmF64gB7xa2UVFMDXX0PLlrB6Nezd62aUKqiLSDQ6bUC31pYBU4APgDW40SyrjDHTjDE/9Bz2n8CPjTHLgdnAHTYKVrOdNw/atYOsLFdWd8cONxJGM0hFJBoFtKaotfY93M1O321TfR6vBgYGt2mhV1Tkarukp8PixfDxxzB5stsuIhJtGvRgPe8M0kaNYMgQl3JZuNBtFxGJNg06oPvOIO3e3fXWFy+Gq68Od8tERGquQQd03xmk27a5BTCOHIH33jv9a0VEIk1AOfRYVnkG6cGD8N//7Zas69w5fO0SEampBt1D92f6dFcWYOrU0x8rIhJJFNAr6doV7rsPZs3SeHQRiS4K6H48+qibbDR6NEycqBmkIhIdFND92LoVevSAwkI4etSNgtEMUhGJdArofsybBwMGQJs28NFHbvk6zSAVkUingO5HUZErCXD11bBzJyxZ4oK6ZpCKSCRTQPfDO4O0Vy834WjBAlfrRTNIRSSSKaD74Z1Bum8fDB0Kx465Oi+jRoW7ZSIiVVNA98N3BumhQ3DZZbB5M+zeHe6WiYhUrcHPFK2K7wzSQ4egTx/4yU9g+XJo0iS8bRMR8Uc99AA0awZ//jOsXevKAoiIRCIF9ABddx2MHAnTpsGWLeFujYjIqRTQa+Dpp8EYuP/+cLdERORUCug1kJICP/4xvPWWG6OukgAiEkkU0GugoAB27XKjX5YscZOOVBJARCKFAnoNzJsHSUkwYoSbeLRkiUoCiEjkUECvgaIiVwIgJQUuvhjy8mDPHpUEEJHIoIBeA96SAABXXunqveTkwJlnhrddIiKggF4jvotKx8W5oL5/P6xfH+6WiYgooNeIb0mA4mLo3RsmTIC5c+HDD8PdOhFp6DT1v4YqLyp96BAsXgy33w533OFGvqSkuN6873EiIqGmHnodNWsGjz3mAvnbb0PnzlrhSETCQwE9CAoLISMDVq2CTZtcSkbDGUWkvimgB0FRkaubnpTkZpF+/71WOBKR+qeAHgQpKXDwoMubHzzohjLu26cVjkSkfimgB4F3OGNCAlx1FaxbB19+qRWORKR+KaAHge9wxk6doEcPyM8Ha8PdMhFpSAIatmiMGQb8EYgD/mqt/Y2fY8YC2YAFlltrbwliOyOe73DGb7+Fvn3hppsgNxfOOCO8bRORhuG0Ad0YEwc8C1wNFANLjDE51trVPsf0AP4LGGitLTHGtA9Vg6NB+/bw8suuxO7tt0NqqrtBqvHpIhJKgaRc+gOF1tqN1tqjwKvA8ErH/Bh41lpbAmCt/Ta4zYw+V10FEye6oYtffqnx6SISeoEE9E7AVp/nxZ5tvnoCPY0xC40xiz0pmlMYY+4yxuQaY3J37dpVuxZHkQ4d4KyzYMEC+O47jU8XkdAK1k3RxkAPIAu4GXjOGNOm8kHW2pnW2kxrbWZycnKQ3jpybdsGY8a4x2+8AWVlGp8uIqETSEDfBnTxed7Zs81XMZBjrT1mrd0EfI0L8A1aSgo0agTDh7vg/v77rvyuxqeLSCgEEtCXAD2MMd2MMU2Am4CcSsfMx/XOMcYk4VIwG4PYzqjkHZ9+1llw6aWwdCl89ZXGp4tIaJw2oFtry4ApwAfAGuA1a+0qY8w0Y8wPPYd9AOwxxqwGFgAPWmv3hKrR0cJ3fHqPHtC9u1vl6PDhcLdMRGKRsWGa/ZKZmWlzc3PD8t7hsmcPZGa6XHpenhveKCJSE8aYPGttpr99qodej9q1gzffhEsugR/8AIYNc7l1jU8XkWDQ1P96lp4OP/+5m0H6zjsany4iwaOAHgZHj7rSAPn5roa6xqeLSDAo5RIGRUUu5bJvn6uf3qaNK+ql8ekiUhfqoYdBSgocOABjx0KrVvDqqydrvYiI1JYCehh4x6cfOQI33wzHj7ubpUOGhLtlIhLNFNDDwHd8+qFDLqh//z08/rjLr4uI1IZy6GHiWz8d4MorXandsWPdDdOtWzWcUURqRj30CHHbbXD33e4m6YcfajijiNScAnoEOfNMOO88WLQIVq/WcEYRqRkF9AiydatLsaSkuJukGzao3K6IBE4BPYKkpMDBg+4maXIyzJkD69ZpOKOIBEYBPYJ4hzMeOgS33ALNm8P8+bopKiKB0SiXCOIdzjhvnkuzTJgAs2fDffdBs2ZubVItNi0iVVH53Ai3YgUMHOgejxvn1iktLXU9+QceUFAXaWiqK5+rlEuES02FG290aZi333YTjzT6RUT8UUCPAta6CUe7dsE//uFKBmj0i4hUpoAeBVJSICkJRo92C2L84x8uuGv0i4j4UkCPAt7RLx06uMfFxTB3LgwdGu6WiUgkUUCPAr7FvFq1cj31vXvhoYdg//5wt05EIoWGLUaJysW8XnvNjVW/4gq3NumOHRrOKNLQqYcepcaOhV/9yi1j9+KL0L69inmJNHQK6FHs++9dHn37dnjlFTf5SMMZRRouBfQoVlQEF10EY8a40S8vvghxcRrOKNJQKaBHsZQUN2u0d29X0Gv3bnjhBTdGXUQaHgX0KOYdzlhSAuecAyNGuGqNs2fD2rXhbp2I1DcF9CjmO5yxuBguuMBNOrIWLr3UrYA0cSJkZ+tGqUhDoOJcMejtt10K5tgxuOkmaNNGxbxEYoWKczUweXluWGOrVq7Hvm2bRr+INAQK6DGoqMgtMj1xovs+dy6sXAlbtoS7ZSISSpopGoNSUlyKJTERxo+HnBz45BO48EJYutQ910IZIrFHPfQY5Dv6pVEjyMqCPn1cML/hBti50/XcNbNUJLYEFNCNMcOMMeuMMYXGmEeqOW60McYaY/wm7KV+VB790raty6Vff72bVfrmm3DggBbKEIk1p025GGPigGeBq4FiYIkxJsdau7rScS2B+4EvQ9FQqZnKxbwAkpPdqJe5c+G559yN006dNLNUJFYE0kPvDxRaazdaa48CrwLD/Rz3S+BJ4HAQ2ydBlJLigvqkSRAf70oFLFqkhTJEYkUgAb0TsNXnebFnWzljzIVAF2vtu9WdyBhzlzEm1xiTu2vXrho3VurGm1uPj3dBvWNH+OgjWL/ejVkXkehW51EuxphGwH8Dd5zuWGvtTGAmuIlFdX1vqRlvbn3ePJdmueMO2LQJXn7ZBfWsLFcPRqNfRKJTIAF9G9DF53lnzzavlkAf4FNjDMBZQI4x5ofWWk0FjTD+cuvnnw9Tp8K6dTBu3MnRL5pZKhJdAkm5LAF6GGO6GWOaADcBOd6d1tpSa22StbartbYrsBhQMI8ix465ErzWumqNmza5cgEa/SISXU7bQ7fWlhljpgAfAHHA89baVcaYaUCutTan+jNIpCsqgvPOc6mWefNcLRh/PXkRiWwB5dCtte8B71XaNrWKY7Pq3iypT74zS2+9FT7/HD77zI1Znz/fLXOnmaUikU8zRaXCzFKAvn1hyBA4cgRuvNENbdTMUpHIp4Aup8wsTUyEp56CH/3ILT79z3+6NMwZZ2hmqUgkU3EuAfznzEtK4M47XQrmiy9c2mXECNi/3/XSvcMflYoRiQzqoUuVUlJc8L7ySjdm/fhxNwpm5Ur47W9dwFcqRiRyKKBLlXxz6126uDowZ58NS5a4crzGuGqOKvIlEhkU0KVKlXPrZ50Fb70Fl18Oe/fC//6vGwFjLbRurSJfIuGmHLpUy19uffBgOOcc10t/6y1YtQoGDVKRL5FwUw9damzUKJdPHz4chg6FzZtdPRhjYPlyyM52y99lZyuvLlKfFNClxrypmLZtXcXGe+5xY9ezs+G6606uaaqbpSL1SykXqZXKqZgTJ9yKSB995FZHGjIEMj3rVs2bpyGNIvVBAV2ColEjd9P07rvh3Xfhvfdc+uW66zRuXaS+KOUiQZOS4gL7bbfByJEu5fLXv0JeHvz61xq3LhJqCugSNN5x6/v2QZ8+MH48dO/uAve778K332rcukgoKaBL0FQet96hg1uQ+rrrICEBXn0VZs92Y9i949YLCjQqRiRYjLXhWQkuMzPT5uZqDYyGIDvbLW23bp2rC3P8OKSnu18AZWXuF0Dr1lBa6nr4WilJpGrGmDxrbaa/feqhS8iNGgXffQe9e8O990LPni6vPmcO7NjhVkdSKkak7hTQJeR8UzGlpTB2LLz0EjRvDh9+CM8/71I0oFSMSF0o5SJh8/jj8K9/uQU0Dh50PfjMTGjWDL7/XqkYEX+UcpGINHq0u3F6222u4Nf69a7nvnAhNG3qArpSMSKB08QiCRtvKmbePLfc3UUXuXTLG2/A11/DZZfBgAHQpEnFVIwmKIn4p5SLRJwpU+D992HjRrfs3WWXufHsx48rFSNSXcpFPXSJOHfd5UbF9OvnFtP44AN3A7VfPxe4ExPdcd7vqhUj4iigS8TxTcW0aOF65gUFbl3TFStcvj0tDeLilIoR8aWUi0QFa93N0w8/dCUE2rSBgQPdknhKxUhDopSLRD1j4KGH3Pd9+2DZMlcfplkzN+s0PV2pGBEFdIkaaWnw4IMuUCcluZ75mjVuHPvy5XDxxW6kTEKCUjHSMCnlIlFv4kS3vumWLW6IY0YG9OrlxrIrFSOxRikXiWk//akr8lVWBitXupExX33lhjpmZvpPxXi/q+cusUQ9dIkJvqmV1q1dT/yVV+DYMeja1U1Q6tnz5LGtWqnnLtFJPXSJeZXXOAW3JN7ChS6Az5njRsZccIEbMXP22af23P/nf9xr1GuXaKVaLhKzbrkFunWD22+HMWPcrNOFC+HLL92Y9p07Tx57+LBb4FrL5Ek0Uw9dYpbvBKWDB+FHP3JL4/361y5QL1sGXbq4GajFxdCunYY+SnQLKKAbY4YBfwTigL9aa39Taf9/AD8CyoBdwERr7ZYgt1WkxvylYnr2hOnT4ZtvYPVqmD/fVXXs3dtNWmrf3h2noY8SbU57U9QYEwd8DVwNFANLgJuttat9jhkMfGmt/d4Ycw+QZa0dV915dVNUwskbpLdscTn1hQtdMbATJ1yvPSMDOnZ0+zT0USJJXW+K9gcKrbUbPSd7FRgOlAd0a+0Cn+MXA+Nr31yR0Kvccy8ogCeegO3bXa89J8fViune3aVpunVzs1Q19FEiWSABvROw1ed5MTCgmuMnAe/722GMuQu4CyAlJSXAJoqEXloaPPaYC9DnnutSMKWlLrB//bWbuJSeDqmpbrRMfr7r0ScmVryJqp67hFNQb4oaY8YDmcAV/vZba2cCM8GlXIL53iJ15S/f/uijkJsLhYXw6afu66yzXBXIrCxNWpLIEkhA3wZ08Xne2bOtAmPMEOBR4Apr7ZHgNE8kvMaNg61b4fzz3fO8PJeS2bEDNmxwKZnUVDjvPJdjV89dwimQm6KNcTdFr8IF8iXALdbaVT7HZABvAMOstesDeWPdFJVo4W+Uy4wZrrxAYaFLzcTFuclKTZvCFVfAmWeefH1JiQvwo0ap5y51V91N0YCm/htjrgOexg1bfN5aO90YMw3ItdbmGGM+AlKB7Z6XFFlrf1jdORXQJZoVFLied5s2bnWl/HyXaz90CBo3hh493DDIHj0gPl7lBiR46hzQQ0EBXaJd5Z77iBHwxz+6CUsbNsCBA67n3rmz67lffjl06nTy9aBole4AAAwNSURBVOq5S20ooIvUE2/PvXVr2L/fPV+/3gV3cAH7vPPcV2Kieu5ScwroIvWocs995EiXc1+2DDZvPllDpk0bN1rmootc0bC4OLddPXepjgK6SJh5e+6JiW72qbfnvmOHm53atKkbMXPuuXDOObBpU9U9d1Cgb8gU0EUigL/RMq++6rZ/803F1MwZZ7ihkhdc4I6Ni3MB/ciRqksRgAJ9Q6CALhKhfHvurVq5Mexr17qgvHev673Hx7shkd26udozQ4dC27Ynz6FA37AooItEMH8993nzXK593z43YmbTJti92x3fvLlLy3Tr5lZjat0a3n7bjX/3zliF0wd6BfXopBWLRCKYv5ID4Hru7du7seylpa5me0mJS89s3uzWTwV3YzU+3vXuu3VzwdsYF8BzcioGepUoiG3qoYtEKH89dzg5oenYMZee2bTJ1XE/dMjtb9nSHZ+U5G66jhnjJjt5nThR/XBJUKCPZEq5iMQQf4HeWpg61Y19//Zb14M/eNAd36SJOy4lxU1yat4c1qyBvn1rlqIBBfpIoIAu0gD4BvouXeCSS2DpUpg1C/bscTdZwaVjWrZ0qZwuXVyQb9vW/VKonKIBBfpIo4Au0oB5A/16T9m8Fi3cgthbt7q0DUBCgsvXHzniAnqXLu44cCma2gR6BfXQ0E1RkQbM303XggL47W9dr7y01KVovvnGjap57TV3TKtW0KHDyQJk3pmsXrrpGnnUQxdpoPzl4g8dgscfd73ukhLYts1992rZ0gX5M890jzdtgtGjKwZ73XQNLaVcRCRglQP90KGwbh08+6y76bp3rxsT7w0dTZq4VZzOPNN9NWvm1ma98ELl4kNBAV1E6sw30Hfo4BbPXrfOlS/47jt34/Xo0ZPHt2njAnz79pCc7IZRLloEgwcr0NeFArqIhIw30G/Z4tIs3bvD/PknJ0Lt2XOyNw/Qrp0L8klJ7qtdO1i8GK68UoE+EAroIlKvfGvUnHGGC/ZFRS44l5a6m69791YM9K1bVwzyiYlu2OWQIapd40sBXUTqXXUzXb2BvqjIDZ88dMilbUpLXY/eO5wS3CzXdu1OBvm2bd1C3VlZbgUoY9xxtQ30/toZyb8AFNBFJGKcLtC3auWGUG7dCocPuxuxBw+6Hn1JiRtF4xUffzLIt2njXtOvnxtH36aN+2VQXaD/4Q/d0Mto6ukroItIxDtdoG/d2gXb4mLXgz9wwAX8kpKTX2VlFc/ZsqUL7Pv3u5mxZ53lnrdpA8ePw6pV0VcCQQFdRKJWIIG+tNT14I1x6ZuyMpen37cPdu1y9W3Kyirm7MGNn/cG+Vat3LlatoT8fDdhqmPHyEvpKKCLSMwJNNCXlLh8fVyc+/LelN2xw92sjY93KZ3SUtdr9xUX5wJ8q1bu+7ffQq9eJydWtWzpflGUldU8pVPboK6ALiINRk0CvW/AbdXKBflvvnGB/fBhl9rZv9/dsN23z333FzLj4tw5vAt/t2jhtu3aBb17u56+t269dxHw7Oza/ftUy0VEGoyqFgx54IGKgX7SJHdcz54Vt993nzve3y+A5s3dTdlGjVwO/7vvXNBev97dmN2/3600deDAycC/erX7PngwXH65O19RUWj+7QroItIgVBXoa/ILAE4G+vbtXaBv0wbS010JBO/NVWtdT3/ZMleeGFyuHtxrUlKC/+8DBXQREb9qE+jhZI/+yBF45JGKOXTvaBzva4JNAV1EpAbqmtLxbg8FBXQRkSCoaUonFBrVz9uIiEioKaCLiMQIBXQRkRgRUEA3xgwzxqwzxhQaYx7xs7+pMWaOZ/+XxpiuwW6oiIhU77QB3RgTBzwLXAv0Bm42xvSudNgkoMRaey7wB+DJYDdURESqF0gPvT9QaK3daK09CrwKDK90zHDgRc/jN4CrjPGWtBERkfoQSEDvBGz1eV7s2eb3GGttGVAKtKt8ImPMXcaYXGNM7q5du2rXYhER8ateb4paa2daazOttZnJycn1+dYiIjEvkIlF24AuPs87e7b5O6bYGNMYaA3sqe6keXl5u40xW2rQ1qokAbuDcJ5gi8R2qU2BicQ2QWS2S20KXLDadXZVOwIJ6EuAHsaYbrjAfRNwS6VjcoAJwCJgDPCJPU1dXmttULroxpjcqkpJhlMktkttCkwktgkis11qU+Dqo12nDejW2jJjzBTgAyAOeN5au8oYMw3ItdbmAH8DXjLGFAJ7cUFfRETqUUC1XKy17wHvVdo21efxYeDG4DZNRERqIhZmis4MdwOqEIntUpsCE4ltgshsl9oUuJC3K2xL0ImISHDFQg9dRERQQBcRiRkRF9DrUgjMGPNfnu3rjDFDAz1nqNpkjLnaGJNnjFnh+X6lz2s+9Zwz3/PVvp7a1NUYc8jnfWf4vKafp62Fxpg/1aZ8Qx3adatPm/KNMSeMMemefaG+VpcbY5YaY8qMMWMq7ZtgjFnv+Zrgs71O16q2bTLGpBtjFhljVhljCowx43z2zTLGbPK5Tun10SbPvuM+75vjs72b53Mu9HzuTWrSprq0yxgzuNLP1GFjzAjPvlBfq/8wxqz2fEYfG2PO9tkXkp8pAKy1EfOFGxa5ATgHaAIsB3pXOuZeYIbn8U3AHM/j3p7jmwLdPOeJC+ScIWxTBtDR87gPsM3nNZ8CmWG4Tl2BlVWc9yvgYsAA7wPX1le7Kh2TCmyox2vVFUgD/g6M8dneFtjo+Z7oeZxY12tVxzb1BHp4HncEtgNtPM9n+R5bX9fJs+9AFed9DbjJ83gGcE99tqvSZ7kXaF5P12qwz3vdw8n/fyH5mfJ+RVoPvS6FwIYDr1prj1hrNwGFnvMFcs6QtMlau8xa+41n+yqgmTGmaQ3eO+htquqExpgOQCtr7WLrfrr+DowIU7tu9rw2GE7bJmvtZmttAXCi0muHAh9aa/daa0uAD4FhQbhWtW6TtfZra+16z+NvgG+BYEzSq8t18svzuV6J+5zBfe5B/5kKsF1jgPettd/X8P1r26YFPu+1GDfDHkL3MwVEXsqlLoXAqnptIOcMVZt8jQaWWmuP+Gx7wfPn3s9r+OdVXdvUzRizzBjzmTFmkM/xxac5Z6jb5TUOmF1pWyivVU1fW9drVdefSQCMMf1xPcQNPpune/7M/0MNOw91bVOCcYX3FnvTGrjPdZ/nc67NOYPRLq+bOPVnqr6u1SRcj7u61wbj/1/EBfSYZIy5AFcj/m6fzbdaa1OBQZ6v2+qpOduBFGttBvAfwCvGmFb19N6nZYwZAHxvrV3pszlc1ypieXp0LwF3Wmu9PdP/As4HLsL9Sf9wPTbpbOumtd8CPG2M6V6P710tz7VKxc1296qXa2WMGQ9kAr8Lxfkri7SAXpNCYJiKhcCqem0g5wxVmzDGdAbeBG631pb3pKy12zzf9wOv4P6MC3mbPCmpPZ73zsP17np6ju/s8/qaXqc6tctn/yk9qXq4VjV9bV2vVZ1+Jj2/gN8FHrXWLvZut9Zut84R4AXq7zr5fkYbcfc8MnCfaxvP51zjcwajXR5jgTettcd82hvya2WMGQI8CvzQ5y/zUP1MObW5KRCqL1wpgo24m5remw0XVDrmJ1S8qfaa5/EFVLwpuhF38+K05wxhm9p4jh/l55xJnsfxuBzj5HpqUzIQ53l8jueHpq31f1Pmuvr6/DzPG3nac059XiufY2dx6k3RTbibV4mex3W+VnVsUxPgY+Cnfo7t4PlugKeB39RTmxKBpp7HScB6PDcJgdepeFP03mD/TFXVLp/ti4HB9XmtcL/QNuC5gR3qn6ny89f0BaH+Aq4DvvZcjEc926bhfssBJHh+SAo9F8D3P/+jntetw+cOsb9z1kebgMeAg0C+z1d74AwgDyjA3Sz9I54gWw9tGu15z3xgKXCDzzkzgZWecz6DZyZxPX5+WcDiSuerj2t1ES5neRDXq1zl89qJnrYW4tIbQblWtW0TMB44VulnKt2z7xNghaddLwMt6qlNl3red7nn+ySfc57j+ZwLPZ970xD8TFX3+XXFdRIaVTpnqK/VR8BOn88oJ9Q/U9ZaTf0XEYkVkZZDFxGRWlJAFxGJEQroIiIxQgFdRCRGKKCLiMQIBXQRkRihgC4iEiP+P98m4qT5fQR0AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# viz.plot_number_observer(obs)\n", + "t_nume = [time_frame[0] for time_frame in obs.data()]\n", + "S_nume = [time_frame[1] / NA for time_frame in obs.data()]\n", + "t_theo, S_theo = calc_S(ka, D, R, concB, delta_t, duration)\n", + "\n", + "plt.plot (t_theo, S_theo, label=\"SCK theory\", color=\"blue\")\n", + "plt.scatter(t_nume, S_nume, label=\"SGFRD simulation\", color=\"blue\", alpha=0.5)\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd3yUVfb/3ze99wAJIaFIT4XQBYIggtgVsYNrQ3+uZe1rR1zdVb/r6rqr2ECxAhZEXEUERQEhQKRHCIQQQkIgvbf7++POTCbJtPSY3PfrldfMPPXOZOY85zn3nM8RUko0Go1G0/1x6uwBaDQajaZj0AZfo9Foegja4Gs0Gk0PQRt8jUaj6SFog6/RaDQ9BJfOHkBjQkJCZP/+/TvkXKmpqQAMHTq0Q86n0Wg07cWOHTtOSylDbW3T5Qx+//79SU5O7pBzJSUlAbBx48YOOZ9Go9G0F0KIY/a20SEdjUaj6SF0OQ+/I3nsscc6ewgajUbTYfRogz9jxozOHoJGo9F0GD3a4KekpAAQHx/fySPRGKmuriYzM5OKiorOHopG0yXx8PAgIiICV1fXZu/bow3+PffcA+hJ265EZmYmvr6+9O/fHyFEZw9Ho+lSSCk5c+YMmZmZDBgwoNn760lbTZeioqKC4OBgbew1GgsIIQgODm7xHbA2+Jouhzb2Go11WvP76D4Gv7wANj4PJ3Z09kg0Go2mS9J9DD7Axufg2ObOHoWmh5CUlGS3SPDll1+mrKzM9Pr888+noKCgvYfWqTz66KP069cPHx+fBstfffVVoqOjOf/886mqqgLg559/5t5777V6rIkTJ7brWNuDL774gv3799vdbunSpWRlZZle33zzzQ7t1xq6j8H38AdXbyg66fAuf/vb3/jb3/7WjoPS/JGRUlJXV9eqYzQ2+GvXriUgIKC1Q3OI2traBq9ramoc2s/R7fLy8iwuv/DCC9m2bVuT5R988AG7d+9m4sSJfPvtt0gpeeaZZ3j88cetnmPz5j+eA9dSg//WW28xYsSI9hxaNzL4QoBfGBSdcHiXiRMn/iE9CE37kZ6eztChQ7nhhhuIjo7m+PHjfPfdd0yYMIFRo0Yxd+5cSkpKmux3++23k5iYyMiRI3nyyScBeOWVV8jKymLatGlMmzYNUNIhp0+f5uGHH+a1114z7f/UU0/x4osvAvDCCy8wZswYYmNjTcdqjLUx9e/fn4ceeohRo0axYsUKkpKSuOeee0hMTORf//oX6enpnHPOOcTGxjJ9+nQyMjIAWLBgAQsXLmTcuHE8+OCDVj8fKSU//PAD11xzDWPGjLG4zfjx4wkLC7O4b3V1NWVlZbi6urJ8+XJmz55NUFCQ1fMZ7xI2btxIUlISV1xxBcOGDePaa6/FUre+N998kzFjxhAXF8fll19uutiuWLGC6Oho4uLimDJlCgBTpkwxpWYDnH322fz222889dRTzJ8/n8mTJxMVFcVnn33Ggw8+SExMDLNmzaK6utr0WRuXjx07lsOHD7N582ZWr17NAw88QHx8PGlpaaSkpDB+/HhiY2O59NJLyc/PZ+XKlSQnJ3PttdcSHx9PeXl5gzvG//3vf4waNYq4uDimT59u9fNpNlLKLvU3evRo2WKWXiDlmzMc3vyXX36Rv/zyS8vPp2lz9u/f3+D11KlTm/y99tprUkopS0tLLa5/9913pZRS5ubmNllnj6NHj0ohhNyyZYvpGJMnT5YlJSVSSimff/55+fTTT5vGtn37dimllGfOnJFSSllTUyOnTp0qf/vtNymllFFRUTI3N9d0fOPrnTt3yilTppiWDx8+XGZkZMhvv/1W3nLLLbKurk7W1tbKOXPmyB9//LHBGG2NKSoqSv79739v8PndfvvtptcXXHCBXLp0qZRSyrfffltefPHFUkop58+fL+fMmSNramosfi4nTpyQzz77rBw2bJi87LLL5Jo1a6xua8Tb27vB6/fee0/Gx8fLa6+9VhYVFclp06bJqqoqh46xYcMG6efnJ48fPy5ra2vl+PHj5aZNm5psf/r0adPzRx99VL7yyitSSimjo6NlZmamlFLK/Px8KaWUS5culXfffbeUUsrU1FRptD1PPvmknDRpkqyqqpIpKSnS09NTrl27Vkop5SWXXCI///xzKaX6rBcvXiyllHLZsmVyzpw5Ukr1Wa5YscI0jpiYGLlx40YppZSPP/646Zzm3x/z16dOnZIRERHyyJEjUsr675Y5jX8nUkoJJEs79rX7ePgAfn2hKMv+dgb++te/8te//rUdB6T5IxIVFcX48eMB2Lp1K/v372fSpEnEx8ezbNkyjh1rqlH16aefMmrUKBISEti3b5/dW/qEhAROnTpFVlYWv/32G4GBgfTr14/vvvuO7777joSEBEaNGsXBgwc5dOhQg33tjWnevHkNtjd/vWXLFq655hoArr/+en7++WfTurlz5+Ls7NxkrNu2bSMyMpKMjAw2bdrEqlWrmDNnjsVtbXH99deza9culi9fzj//+U/uuusuvvnmG6644gruvfdeu+GzsWPHEhERgZOTE/Hx8aSnpzfZZu/evUyePJmYmBg++OAD9u3bB8CkSZNYsGABb775pinUNXfuXNasWUN1dTXvvPMOCxYsMB1n9uzZuLq6EhMTQ21tLbNmzQIgJiamwXmvvvpq0+OWLVuajKewsJCCggKmTp0KwPz58/npp59svs+tW7cyZcoUU569rTug5tK9Cq/8wqEkG+pqwal5X0ZN18RWUZyXl5fN9SEhIS0qqvP29jY9l1Jy7rnn8tFHH1nd/ujRo7z44ots376dwMBAFixY4FCe9Ny5c1m5ciXZ2dkmoyyl5JFHHuG2226zup+9MZmP39Jra1jbLjY2lrfffpu3336biy++mAULFjBv3jz8/PwcOm5jsrKy2LZtG0888QRTp07lhx9+YPHixaxfv55zzz3X6n7u7u6m587OzhbnGhYsWMAXX3xBXFwcS5cuNf3/X3/9dX799Ve+/vprRo8ezY4dOwgODubcc8/lyy+/5NNPP2XHjvoMP+O5nJyccHV1NaVCOjk5NTiveYrkHyGduJt5+OFQVwOluZ09Ek03Yfz48fzyyy8cPnwYgNLSUn7//fcG2xQVFeHt7Y2/vz85OTl88803pnW+vr4UFxdbPPa8efP4+OOPWblyJXPnzgXgvPPO45133jHF5E+cOMGpU6eaPSZrTJw4kY8//hhQk6iTJ0+2u4+Hh4fJM126dClpaWkkJCRw3XXXOXTOxjz++OMsWrQIgPLycoQQODk5NZjcbinFxcWEhYVRXV3NBx98YFqelpbGuHHjWLRoEaGhoRw/fhxQmTF33XUXY8aMITAwsNnn++STT0yPEyZMABr+z/39/QkMDGTTpk0AvP/++yZv39p3Y/z48fz0008cPXoUsD453hK6l4fvG64ei06Ab5/OHYumWxAaGsrSpUu5+uqrqaysBGDx4sUMGTLEtE1cXBwJCQkMGzaMfv36MWnSJNO6W2+9lVmzZhEeHs6GDRsaHHvkyJEUFxfTt29f0yTnzJkzOXDggMl4+Pj4sHz5cnr16tWsMVnj1Vdf5cYbb+SFF14gNDSUd999t1mfx+DBg3n++edZvHgxX3/9tcVtHnzwQT788EPKysqIiIjg5ptv5qmnngJg165dAIwaNQqAa665hpiYGPr162dzsthRnnnmGcaNG0doaCjjxo0zGdQHHniAQ4cOIaVk+vTpxMXFATB69Gj8/Py48cYbW3S+/Px8YmNjcXd3N91xXXXVVdxyyy288sorrFy5kmXLlrFw4ULKysoYOHCg6TM3TpR7eno2CAeFhoayZMkSLrvsMurq6ujVqxfr1q1rzcdiQkgLM92dSWJiomxxA5SsFFgyFeYth+EX2t1cN0Dpehw4cIDhw4d39jA0PYSsrCySkpI4ePAgTk7NC3gYmzWFhIS00+isY+l3IoTYIaVMtLVf9/Lw/fqqRwdz8V9++eV2HIxGo+nKvPfeezz66KP83//9X7ON/R+V7mXwvYLB2c3hXHwti6zR9FxuuOEGbrjhhhbvbylLqKvTvS5rTk4qdu9gaub333/P999/386D0mg0mq5B9/LwoVm5+IsXLwZ05yuNRtMzcMjDF0LMEkKkCiEOCyEetrD+L0KI/UKI3UKI9UKIKLN184UQhwx/89ty8BbxC4dix4uvNBqNpqdg1+ALIZyB14DZwAjgaiFEY4WfXUCilDIWWAn8w7BvEPAkMA4YCzwphGh+smtz8AtXHn4Xyz7SaDSazsYRD38scFhKeURKWQV8DFxsvoGUcoOU0lg1sRWIMDw/D1gnpcyTUuYD64BZbTN0K/iGQ00FlOe362k03Zdnn32WkSNHEhsbS3x8PL/++iug+u0+/PDDDB48mFGjRjFhwgRTkZVRFA1gx44dDBgwwJRzbiQlJYW1a9eaXpsLpmk0HYEjMfy+wHGz15koj90aNwHGUkNL+/ZtvIMQ4lbgVoDIyEgHhmQDP7PiK6+206DQ9Ay2bNnCmjVr2LlzJ+7u7pw+fdqk3f74449z8uRJ9u7di7u7Ozk5Ofz4448N9t+9ezdXXHEFn3zyCQkJCQ3WpaSkkJyczPnnn9+u76G2trbZOjeankGbZukIIa4DEoEXmrOflHKJlDJRSpkYGhraukE0Ixf/jTfe4I033mjd+TTdipMnTxISEmLSUgkJCSE8PJyysjLefPNNXn31VdO63r17c+WVV5r2PXDgAJdccgnvv/8+Y8eObXDcqqoqnnjiCT755BPi4+NNJfn79+8nKSmJgQMH8sorr5i2X758OWPHjiU+Pp7bbrvNJPj10UcfERMTQ3R0NA899JBpex8fH+677z7i4uJ49tlnueSSS0zr1q1bx6WXXtrGn5Tmj4gjHv4JoJ/Z6wjDsgYIIWYAjwJTpZSVZvsmNdp3Y0sG6jDmHr4dhg4d2q5D0bSOp7/ax/6sojY95ohwP568cKTV9TNnzmTRokUMGTKEGTNmMG/ePKZOncrhw4eJjIy0KRh28cUXs3z5cs4+++wm69zc3Fi0aBHJycn8+9//BlRI5+DBg2zYsIHi4mKGDh3K7bffzuHDh/nkk0/45ZdfcHV15Y477uCDDz5gxowZPPTQQ+zYsYPAwEBmzpzJF198wSWXXEJpaSnjxo3jpZdeQkrJ8OHDyc3NNckn/OlPf2r9h6f5w+OIh78dGCyEGCCEcAOuAlabbyCESADeAC6SUporPX0LzBRCBBoma2calrUfPr1BODmUmvnVV1/x1VdftetwNH8sfHx82LFjB0uWLCE0NJR58+axdOlSh/adMWMGb731VpNOU7aYM2cO7u7uhISE0KtXL3Jycli/fj07duxgzJgxxMfHs379eo4cOcL27dtJSkoiNDQUFxcXrr32WpPUrrOzM5dffjmgVBuvv/56li9fTkFBAVu2bGH27NnN/iw03Q+7Hr6UskYIcSfKUDsD70gp9wkhFqEE91ejQjg+wAqDRGiGlPIiKWWeEOIZ1EUDYJGUsu2k3yzh7KKMvgMG/6WXXgJUSzZN18OWJ96eODs7k5SURFJSEjExMSxbtowrr7ySjIwMioqKrHr5//73v1m4cCF33HGHw6FCS5K/Ukrmz5/Pc88912DbL7/80upxPDw8GsTtb7zxRi688EI8PDyYO3cuLi7dr+RG03wciuFLKddKKYdIKQdJKZ81LHvCYOyRUs6QUvaWUsYb/i4y2/cdKeVZhr/mSfO1FJ2Lr2khqampDRqOpKSkEBUVhZeXFzfddBN33323aRI3NzeXFStWmLZ1cnLiww8/5ODBgzzxxBNNjm1LKtmc6dOns3LlSpMscl5eHseOHWPs2LH8+OOPnD59mtraWj766COT1G5jwsPDCQ8PZ/HixS1WgtR0P7qXtIIRYy6+RtNMSkpKmD9/PiNGjCA2Npb9+/ebpH0XL15MaGgoI0aMIDo6mgsuuKCJt+/h4cHq1atZvXp1g561ANOmTWP//v0NJm0tMWLECBYvXszMmTOJjY3l3HPP5eTJk4SFhfH8888zbdo04uLiGD16NBdffLHV41x77bX069dPq49qTHQveWQjax+E3z6CR47b3EzLI3c9tDxy23HnnXeSkJDATTfd1NlD0bQxWh7ZHL9wqCyCymJw9+3s0Wg0Hc7o0aPx9vY2zVNpNNBtDb5ZLn6odYP//vvvd9CANJqOxbw/q0ZjpJsafLNc/FDrbd/69etndZ1Go9F0N7rppK3qD2pv4vaTTz6xOXmm0Wg03Ynu6eGbmpnbNvj//e9/AZg3b157j0ij0Wg6ne7p4bt6qHaHOhdfo9FoTHRPgw86F1/TYqzJI9fU1PDXv/6VwYMHEx8fT3x8PM8++6xpP2dnZ9Py+Ph40tPT2bhxI/7+/sTHxzNs2DDuv/9+0/ZLly4lNDSUhIQEBg8ezHnnncfmzZsdHmdycjJ33XVXm7znpUuXcuedd9rcZuPGjQ3G9/rrr/Pee++1yfk1HUP3DOmAodWhY83MNRojtuSRH3vsMbKzs9mzZw8eHh4UFxc3SHv09PQkJSWlwfHS09OZPHkya9asoby8nISEBC699FImTZoEqHCiUUxtw4YNXHbZZWzYsMGhWoTExEQSE22mXbcpGzduxMfHh4kTJwKwcOHCDju3pm3ovh6+b5j28DXNxhF5ZA8PD0BJJRircB3B09OT+Ph4Tpyw7IhMmzaNW2+9lSVLljRZt2LFCqKjo4mLi2PKlCmAMsAXXHABoJQ358+fz+TJk4mKiuKzzz7jwQcfJCYmhlmzZlFdXQ00bNSSnJxsKj4056uvvmLcuHEkJCQwY8YMcnJySE9P5/XXX+ef//wn8fHxbNq0qUEDl5SUFMaPH09sbCyXXnop+fmqAVFSUhIPPfQQY8eOZciQIWzatMnhz0vT9nRvD7/sDFRXqJi+BVauXNnBg9I0i28ehuw9bXvMPjEw+3mrq+3JI/v6Wq/rKC8vJz4+HoABAwbw+eefN1ifn5/PoUOHTAbbEqNGjbIovLZo0SK+/fZb+vbtS0FBgcV909LS2LBhA/v372fChAmsWrWKf/zjH1x66aV8/fXXDTTybXH22WezdetWhBC89dZb/OMf/+Cll15i4cKF+Pj4mMJS69evN+1zww038OqrrzJ16lSeeOIJnn76aV5++WVAhcK2bdvG2rVrefrpp/n+++8dGoem7em+Hr4xF7/YeiOUkJAQQkJCOmhAmj8Cjsojv/vuu8THx9OvXz+OH1cSHsaQTkpKSgNjv2nTJuLi4ujbty/nnXceffr0sXp+a1InkyZNYsGCBbz55ptW5Zdnz56Nq6srMTEx1NbWMmuW6iYaExNDenq6g58AZGZmct555xETE8MLL7zAvn37bG5fWFhIQUGBScht/vz5JtlmgMsuuwxQ1b/NGYem7enGHr5ZambQAIubGH/ICxYs6JgxaZqHDU+8PbElj1xcXIyvry833ngjN954I9HR0Xb1740x/KNHjzJ+/HiuvPJK051AY3bt2mUxfv/666/z66+/8vXXXzN69GiLlbTGMJSTkxOurq4YpMpxcnKipqYGABcXF+rq6gCoqKiwOIY///nP/OUvf+Giiy5i48aNzQpbWcI4LqP8s6bz6P4evo04/tKlSx1ubqHpGdiTR77zzjtNhrK2ttY0oesIAwYM4OGHH+bvf/+7xfU//vgjS5Ys4ZZbbmmyLi0tjXHjxrFo0SJCQ0NNdxXNpX///qaLxapVqyxuU1hYSN++Sp5k2bJlpuXW5J39/f0JDAw0xefff/99q7LNms6l+xt8nYuvaQa25JGfffZZwsLCiI6OJiEhgcmTJzN//nzCw8MdPv7ChQv56aefTKENY4/bIUOG8Le//Y1Vq1ZZ9PAfeOABUy/biRMnEhcX16L39+STT3L33XeTmJhotdH5U089xdy5cxk9enSDkOeFF17I559/bpq0NWfZsmU88MADxMbGkpKSYrEfgKbz6Z7yyEae6wfx18Bsyx6Vlkfuemh5ZI3GPi2VR3bIwxdCzBJCpAohDgshHrawfooQYqcQokYIcUWjdf8QQuwTQhwQQrwijIHFjsAvXOfiazQajQG7Bl8I4Qy8BswGRgBXCyFGNNosA1gAfNho34nAJCAWiAbGAB0X3NO5+BqNRmPCkSydscBhKeURACHEx8DFwH7jBlLKdMO6ukb7SsADcAME4ArktHrUjuLXF9JSra5eu3Zthw1F4zhSSjryRlCj+SPRmjC8IyGdvoB5SkCmYZldpJRbgA3AScPft1LKA423E0LcKoRIFkIk5+bmOnJox/ALh5JsqLWcCubl5YWXl1fbnU/Tajw8PDhz5kyrvtQaTXdFSsmZM2dM1d7NpV3z8IUQZwHDgQjDonVCiMlSygZT/FLKJcASUJO2bTYAv3CQdVCSA/5Nr1H/+c9/ALjjjjva7JSa1hEREUFmZiZteuHXaLoRHh4eRERE2N/QAo4Y/BOAeWuoCMMyR7gU2CqlLAEQQnwDTAA6RlDDPBffgsH/9NNPAW3wuxKurq4MGGC5UE6j0bQOR0I624HBQogBQgg34CpgtYPHzwCmCiFchBCuqAnbJiGddkPn4ms0Go0JuwZfSlkD3Al8izLWn0op9wkhFgkhLgIQQowRQmQCc4E3hBBG8Y2VQBqwB/gN+E1K+VU7vA/LmJqZa4Ov0Wg0DsXwpZRrgbWNlj1h9nw79XF6821qgdtaOcaW4xkILh46F1+j0WjoztIKAELoXHyNRqMx0H3VMo349YUiyxLJWlJBo9H0JLq3hw9aXkGj0WgM9AyDX3wS6hoXAcOLL75oatGm0Wg03Z2eYfBrq1S7w0asWbOGNWvWdMKgNBqNpuPpGQYfdC6+RqPp8fQcg68zdTQaTQ+nBxh8Y/GVnrjVaDQ9m25p8Esra3hw5W+s3XMSvENBOFv08D09PfH09OyEEWo0Gk3H0y3z8L3cnNl06DRF5TWcHxOmiq8Km3r433zzTSeMTqPRaDqHbunhCyGYNqwXmw7lUlVTBwGRUHjc/o4ajUbTjemWBh/gnKG9KK2qZdvRPAiMgvxjTbZ55plneOaZZzphdBqNRtPxdFuDP+msENxdnPjh4Cnl4RdnQU1Vg23Wr1/P+vXrO2mEGo1G07F0W4Pv6ebMhEHB/HAwBwKiVOeroszOHpZGo9F0Gt3W4AOcM6wX6WfKyBKhaoGFsI5Go9H0FLq1wZ82tBcAP50yNCovyOjE0Wg0Gk3n0q0Nfr8gL4b09mFNugAnFyho6OEHBwcTHBzcSaPTaDSajsUhgy+EmCWESBVCHBZCPGxh/RQhxE4hRI0Q4opG6yKFEN8JIQ4IIfYLIfq3zdAd45xhvdmaXkidb98mHv6qVatYtWpVRw5Ho9FoOg27Bl8I4Qy8BswGRgBXCyFGNNosA1gAfGjhEO8BL0gphwNjgVOtGXBzOWdYL2rqJHluYTqGr9FoejSOePhjgcNSyiNSyirgY+Bi8w2klOlSyt1AA9F5w4XBRUq5zrBdiZSyrG2G7hijIgPw93QlrTqoSUjnkUce4ZFHHunI4Wg0Gk2n4Yi0Ql/AvEw1Exjn4PGHAAVCiM+AAcD3wMOG5uYdgouzE1OHhLLzkB/j6nKguhxclX7Oli1bOmoYGo1G0+m096StCzAZuB8YAwxEhX4aIIS4VQiRLIRIzs3NbfNBnDOsF6kVQepFgZZY0Gg0PRNHDP4JoJ/Z6wjDMkfIBFIM4aAa4AtgVOONpJRLpJSJUsrE0NBQBw/tOFOHhHICw3F1aqZGo+mhOGLwtwODhRADhBBuwFXAagePvx0IEMJY+cQ5wP7mD7N1BHq7ERA+WL0oSO/o02s0Gk2XwK7BN3jmdwLfAgeAT6WU+4QQi4QQFwEIIcYIITKBucAbQoh9hn1rUeGc9UKIPYAA3myft2Kb+BFDqZSulOYcMS2LiIggIiKiM4aj0Wg0HY6QUnb2GBqQmJgok5OT2/y4B7OLcP3PWFzDY4hcuKLNj6/RaDSdiRBih5Qy0dY23brS1pyhvX057dybOp2Lr9Foeig9xuALIRBBUfhVZFFZo7JC77nnHu65555OHplGo9F0DD3G4AMERQwmSBSzPVWlZqakpJCSktLJo9JoNJqOoUcZ/MgBwwDYvXd3J49Eo9FoOp4eZfDdQgYCsG/fHk4VV3TyaDQajaZj6VEGn4BIAPrIHF78NrWTB6PRaDQdiyNaOt0H7xBw9WJm70qu2pHJuOFj8Kst7OxRaTQaTYfQswy+EBAQxSi/IoK83KiNupQ3bpvQ2aPSaDRdhU/ng5sPXPJaZ4+kXehZIR2AgEhci45z/3lD2Z6ez9d7Tnb2iDQaTVfh5G+Qua2zR9Fu9DyDHxgFBRlcmdgP35pC7n//FyqqO0ytWaPRdGVKTimBxS6mQNBW9DyDHxAJlYU4VxbgnbqWCmcvlvx0xP5+Go2me1NZDNWlUFMBJTmdPZp2oQca/Cj1WJCBR3EmXmdS+e/GNLILdZqmRtOjKTHrvtpNJVh6nsEPNBh8wz80MONHaqXk7/872ImD0mg0nY65V1+gDX73wJCLb/yHulYWcsvkAXy+6wQ7M/I7cWAajaZTMTf42sPvJngGgrs/FGQQHx9PfHw8dySdRS9fd55avc8krKbRaHoYxpCOi0e3bZTU8ww+KC8//xgvv/wyL7/8Mt7uLjx90Uh2ZxZy/4rd1NVZmaGvq4V9n6tHjUbTvSjJAeEMvaO1h9+tMKRmmjM7JoyHZg3jq9+y+Pu3VuL5e1bAigVwaF37j1Gj0XQsJTng0wuCBugYfrciIBIKjnHddddy3XXXmRYvnDqQ68ZH8saPR3hvS3rT/Xa+px5PdXhbXo1G096UnAKf3iqTr/AE1NZ09ojaHIcMvhBilhAiVQhxWAjxsIX1U4QQO4UQNUKIKyys9xNCZAoh/t0Wg241AVFQXUZJTjqZmZmmxUIInr4omhnDe/PU6n18ty+7fp/Th+HYL+p5rhZe02i6HSU5yuAHRoGshaJM+/v8wbBr8IUQzsBrwGxgBHC1EGJEo80ygAXAh1YO8wzwU8uH2cYYUjP7eFQ1WeXsJHj16gRiIgK46+Nd9R2jMHwAACAASURBVJk7u94zxPdiIFencGo03Y6SUyqkE9Awdbs74YiHPxY4LKU8IqWsAj4GLjbfQEqZLqXcDdQ13lkIMRroDXzXBuNtGwypmWEWDD6Ap5szb89PpJevBzcvS+ZoTgGkfAhDZ8OAyXD6d6hr8lY1Gk1zyD8GtdWdPQpFXW19SMdYq9MN4/iOGPy+wHGz15mGZXYRQjgBLwH329nuViFEshAiOTc315FDtw6jLr4Vgw8Q4uPOsj+NRUrJ2+/8F0pzYdQNEDoUqsugMMPqvg5RWazihBpNT6SiEF4bC9ve7OyRKMryVBjHpzf4Rai7+fz0zh5Vm9Pek7Z3AGullDaDYVLKJVLKRCllYmhoaDsPCXD3Bc8g4vsHMGGCdXnkASHevL1gDOdWfkcOQWx3SYBQ1Sax1XH8bx+F/4yHoqzWHQdgz0o4k9b642g0HUXOPqVZc/zXzh6Jwlh05dMLnF3Av2+PDemcAPqZvY4wLHOECcCdQoh04EXgBiHE880aYXsRGMW4IWE899xzNjcb5V/GFJHCOtfpXPN2Mp9n+qgVrY3jH/0JKovg6/tbp8xXcgpW3QTfPNi68Wg0HUnOPvWY3UX6S5sMfm/1GBDVY0M624HBQogBQgg34CpgtSMHl1JeK6WMlFL2R4V13pNSNsny6RQc/YemfIiQdVx040OMGxDMvauPUewaQt2pVhj84hzIPwohQyD1a9j/ZcuPdfj7+kft5Wv+KGTvUY95R1R4p7MxVtn69FKPgVE908OXUtYAdwLfAgeAT6WU+4QQi4QQFwEIIcYIITKBucAbQoh97TnoNiEgkurcNK64/DLr29TVqeycAVPxCx/M0hvHsGBif36r6E36gR0UVbRwwun4VvV44SsQFgdrH4DyFur4HFqn5CKcXGD72y07hkbT0eTsA2d39dxo/DuTJh5+fyg9BVVlnTak9sChGL6Ucq2UcoiUcpCU8lnDsieklKsNz7dLKSOklN5SymAp5UgLx1gqpbyzbYffCgKjcHWSyGIbutdHf1QVuaNuAMDF2YmnLhpJyMA4elemc9lrv/B7TnHzz52xVel19B0NF70KZWfgu8ebf5zaGkhbD0PnwPALIWU5VJU2/zgaTUdSV6eKF4fNUa9PdoGwTskp1drQ3RCyDayXUe9O9MxKWzDl2oZ5VFrfZtf7ynsedkGDxcOiE/EWlbiVZjHnlU28/P3vVNU0I00zY6sy9i5uysOf+Gd1riM/Nu89nEhWt8ODz4Wxt6rne1Y07xgaTUeTf1Rlug06B3z6qLaCnY1RVsFIQPdMzezxBt9qamZZHhz4CmKvAlePhut6DQfg40sDOT8mjJe/P8SFr/7MLkfklatK1Rc8cnz9sqSHIXAAfHU3VJc7/h4OfafSxwYmQeQE6DUStr3VbduzaboJxgnb3iOVw9MVJm6NVbZGAvurx24Wx+/BBl8lHoV5WjH4uz+B2ioYdX3TdYbUTL/iNP51VQLvLEikqKKay/67mWfW7KesyoYGx4kdKt+3n5nBd/WEi15Rns/GZiQxHfpOXTg8A0AIGHsz5OzpOqluGo0lcvaBcFKOU1isynjr7Fh5Yw/fpxe4eGoPv9vg6kkx3owaGNx0nZRKKK3vaOWFNMYrCLxDIfcAAOcM6813907h2nGRvP3zUc57+Se+35+DtORpZ2wFBPQb03D5gCmQcD1sftWxW9yik2qy66wZ9ctirlRa/12lmEWjsUTOXgg+Szk6YXEg6zpfkLCxhy+EQUY9vdOG1B70XIMP+EaMJCG4Grb+Vxnan/8JP72giqJO7TdN1lokdFiD4itfD1cWXxLDp7dNwNXZiZvfS+by/25mS9qZhvtlbFWejWdg02POfAa8gmH1n+1r7hvTMQfPrF/m7gPx16g0T1uT0RpNZ5Kzt96R6hOrHjszjl9doea/zD18MMioaw+/+xCeAKdT4X8Pw3ePwfdPwQ+LYetr4NcXRtpI2Qwdqgx+Iy9+7IAgvr1nCs9dFkNWQQVXv7mV69/+ld2ZBcqIH9/WMH5vjmcgzHpOffn3fmZ77Ie+A9/wpncgY26GumrYucz++9doOprKYuU1G7+3AZHgEdC5Br/UmIPfu+HygCjI715ZOi6dPYDOZParu/FzGcMnn65UeexOLuDkrCZCnexcC0OHqUrZ4pPgF95glauzE1ePjeTShL4s33qM1zYc5qJ//8LNg0t5rKpYTbBaY+RlsOkldacRfZkaT2Nqq+HIRhh5ibr1NCfkLJX9kPwunH0vOLs69mFoNB3BKRUGpXe0ehRChXUcMfhG56rxd761mIqu+jRcHhgFlYWqRsbSHfkfkB7t4ZeXV5BTXKMmPd19VDaOs6t9Yw9mmjrWK249XJ25efJAfnpwGvfMGAwZWwD4y1YPfvw913KM38kJpjyg7jysVeAe/1VdbMzDOeaMuQWKs+Dg1/bfhzkndsC/4qGkDQTs/vdXWHVz64+j6V7k7FWP5nemYbEqhGpPOfP7J2HpnLYfk7mOjjndUCa5Rxv8VmE0+A5ILPh6uHLPjCE8ODyfErdebDrlyfx3tnHuP3/iw18zqKhuFK8fcTGEDFVeviUZ5kPr1N3IgKmWTzjkPPCPhO1vNe89HVyrMoUyNjdvP0v8/o3SC+qu/P6tusvSNI/svSqxwN9MnissXmXE2dKnqq2BXcsNWW5tnHbcuMrWSDeUSdYGv6V4h4BnULNE1NyytuMzeBI/P3wO/3dlHO4uTvz18z1MeG49//jfQdJPG6pknZxhyv3K60m14KUfWqfCQh5+lk/k5Axj/gTpm+pvoR0hc7t6bG08tbIE8o6qH1J1ReuO1VVZ94Tqb1xe0Nkj+WORs0959+ZhmbA49Wir4vbYz6oivaai5TIk1ig5BQj1mzZHe/gaE0I0ydSxScFx1TItcgLuLs5cNiqCNX8+m49vHU9i/yD++2MaSS9u5MrXt/Bp8nFKB18EQYPgx7839GgKM+HUPuvhHCMJNyitEke9/LpaOLFTPc9KcWwfa5zaDxjGXNiF2sSV5KpGNm1B0UlleH55uW2O1xOQst7gmxM0CFy9bTsa5uHN4mzr27WEkhyVHdd4vsszADz8u1VqZo82+BdccAEXXHCB/Q2t0WuY8vAducU0FkP1G2daJIRg/MBg3rwhkc0Pn8MD5w0lt6SSB1fuZsxzG/nY80rI3kPdwW/qj2NKxzzX9vm8g1VoJ/V/jo0vNxWqisHdT/3wWnPbbIzTQusbxbQlO5fCF7eri29rqCpVk3nO7iqlty16GvQECjLUd6yxwXdygj4x1g1+Xa2qeveLUK+L2/jzLs5pGs4x0hqZ5Kxd9tOrO5gebfDvv/9+7r/fZjMu24QOg4qC+ll+W2RsUeJMxuyERoT5e/L/pp3FD/dNZeXCCVwUF87zmTEcq+vFwU8f48kv9rD1yBnk79+p+KdxDsEWA6aouwpHvrDGcE78NVB2unVGLHuvqqSE1hvXtsQoH513pHXHMXqYk+9TRUMbbfdU0BgwSSpY+A2ExalCQktzVsc2q45zY29Rr9vDw288YWukpTLJeUdgSVKX07bq0Qa/1YQOVY+5DsTJM36FiETVTccGQggS+wfx/OWxbH70PPJG/ZkRMo2sHV9x/ZKfKUtdzzaXUWxOO0N1rR3Btv5nq8f0n+2PL3ObSj0z1h6cbEVYJ2evqlIWzl1LbbDNDP5J9Rg5TtU97Fru0OR9j8do8A1aVA0Ii4XqUsiz0NNh/5dK5sBYCFl0sm3HZexla4mAKPUdbm4P69OH1KMxTNpF6NEGPykpiaSkpJYfwNF2hxWFygjayr+3gJebCwkXLAT/SF7vt573ZtTgTQVLTw3hmrd+ZdQz6/h/H+5k5Y5MTpdYUP0MHaZik+m/2D9ZZjJEjFG31sKp5RO3dXWQs19lXviFQ2EX8vCNxsSSUWkORoPjGwaT71fx5/WL7O9XY0OZtSeQs0eJBBoliM0xTdw2+t7V1cKB1SqE6RWknJK2DOlIqTx8XysGP7A/1FbWZ/I4ivGuwDy82QXo0YVXrcant5rUsZepk7kdkA3i9w7j4gaT78V5zb1M4N/g7MZL9/+Zi9LL+OHgKTak5vL17pMIAbERAUwbGsrUIaHE9PXHxdkJoibZ9/DLC9R7iL4c3LxUSmhLDX7BMRWn7ROtJm+7SkinvEBleYDKIGoNxWYG38MPzr5bVWhnbLVcRV1drpRQD6yBu3aCb5+m2/QEcvap74UlQoeBs5v63sVcUb/8+K/K2I64WL32DW/bkE5FoTLo1jx8o2pmwTHwC3P8uMaJ3uy96qLS1sViLaRHe/itRggIHW7fw8/YqsIbEYktO0/8tUrqIWsnRE3E08efWdFh/OOKOH59ZDpr/nw2984YgpOAf60/xKX/2UzConXcvGw7W2qHQ2EG0lamQZbhttM4vvD4lmfqmAprotVcQ1cJ6Ri9eheP1reCLD6pvHp3X/V6/B3KYKx7sulkd1EWvDtbqa9WlyppjZ5IVZn63K3MYeHsCr1GNHU09n+p/mdDzlOvffu07SR5iRVZBSMtTc00zptVFnapu1yHDL4QYpYQIlUIcVgI0aQnrRBiihBipxCiRghxhdnyeCHEFiHEPiHEbiHEvLYcfJcgdKh9Dz9jq/JsjAaiubi4w6R71PNG6ZhOToLovv7cNX0wn98xiR2Pncu/r0nggrhwDp0q4ak9qiT8qVff5M4Pd/LelnT2ZxVRW2dmmDKTAaHi7qBur0uyW+ZJZe9Vx+o1XElQF2fZr6DsCM4Y4vb9z1bFZc2NyZpTfFIZHqPX5uatehoc3wqpZhlVx7epibvTh2DuUlUsl7Wr5ef9I5N7AJCW1WeNGLXxjRfNujpl8M+aUf/b8QtrWw/fWpWtkYBI9djcTJ38dPAy5PVnd52wjt2QjhDCGXgNOBfIBLYLIVZLKc31TDOABahG5eaUATdIKQ8JIcKBHUKIb6WU3adaJXSYEioryQWf0Kbra6uVQR09v3XnGT0fasqVt2+DIG83LogN54JYpe+TmTeGyv8+x/leadxzLJ81u1U4wtfdhVFRgYwdEMQ1h3/BP2QoTh7+6iDm8dTmhh9y9kLwIGUEAyJVFktRVn3VYmeRlwYIGDRdpbYWnwT/vi07VnF2E/0kEm6ALa8pAb7BM+G3j+Drv6g7sxu+VBfATf/Xcw2+edMTa4TFqt9S4XH13cncrv5PIy6p38Y3XImd1dbYTYBwCGtVtkZcPZTGTnM8fCnV9iMuhpQP1G9i2PmtH2sb4MgnNhY4LKU8AiCE+Bi4GDAZfCllumFdA7dJSvm72fMsIcQpIBToEgb/yiuvbP1BTJk6By0b/JO7laFuSfzeHBd3mHR3s3eLCPKBQZMZl72HzQ+fw4mCcran57E9PZ/k9Dxe+PYg17gns6IukXdf/omEyAASwwK5DIHMSsHJeCvtKDl71YQt1JfPF2R0vsE/kwb+EfUZInlpLTf4RVnQb2zDZc4uMP1J+PR6pfdyfCsMnAZXvKMmG0Gpsx5Y3aViuh1Gzj4VBgvob30b4/fm5G/K4O//UsX1zb+Dvn2UE1F6qulFtyWYQjpWPHxovkxyWV59vUHQgK7RpN2AIwa/L2AehMoEmm29hBBjATegSQBVCHErcCtAZGRkcw/dYu64447WH8RcRG3A5Kbrj29Vj9YkkTuC/mfDwTWIwkwiAvsREejFpQmqiKXw+AH83y4hcMhEQqvc+Xr3ST7aVkO8Wx/SN6zjjdQpRIf7ExPhR3S4PwNDfXB2smKsKorUrWz8deq18Xa4K8Qw89IgaKC6+wCVmjlgSvOPI6Xy8H0tTOANvxD6Jqr/+YQ7YcbTDb3Q8HjlwRYcq58M7Clk74XeI2wLE/Yeqea6Tu6GoXOUwR80vaGEiNHIFzVVqW0RJdnqouIRYH2bgChD4yIHKUhXj4FRas6iC2XqdEiWjhAiDHgfmC+lbBI8lVIuAZYAJCYmdlhD1rIy1VbNy8ur5QfxC1fVqZYmbqVUjckDItvmy9lSjPn4x36BgKsarPI/oyZnZ868kJm9R1BXJzl6phTXLxIYfWoHNbV1fLjtGBW/qH+bp6szw8N8GRHux/AwP4b18WNoH1983F3qdXuMmRj+hsrI9sjUkRI+nKdy4SffZ3/bM4ch+goVYnF2a3kufnm+yuqwZPCFgHnvq3NZupiEJ6jHrF09y+BLqYzeyEtsb+fqCSFDlIeftVMVDU5/vOE2xhBjcRvl4htz8G3dcQVGwd6VKjzriNy4MUEisL9Kcz7wldKXspSO2sE4YvBPAGbSdkQYljmEEMIP+Bp4VErZjMtk+3P++SqutnHjxpYfRAjLE7dVpSoV79C3MPHPLT9+W9BrpPJg0n+GuIYGn8zt4OZrCk05OQkGhfrAyIlwYi2fzR9KjUcQR06XsvdEIXtPFLH3RCFf7spi+db6DJyoYC9u9dzAtcD3eaH0yy5mQIg3bj592idTZ/8X6rMtybFv8MvyVPpd8CAlLBfYv+WZOkZDYy1Fzy/c+sW91wh1scnaBSMvtX+uvaugdwyEDmnZWLsKRVmqIt1aho45YXFKhXT/F+DkCkNmNVzva/hs28zg26iyNRIQpcJIhZkqRGMPY7w/wODhI1WKcuMwYCfgiMHfDgwWQgxAGfqrgGscObgQwg34HHhPSrmyxaPs6oQOVXK5RnJ/V7Hc3FQ45zE4245Bam+cbOTjZ26HiNFNG62YTdy6nDWdIb19GdLbl8tGqcVSSk4UlHPgZDEHTxZxILsI36MHKZRe3PxlNpCDi5Ngtac/zqn7WLvud87q5cOgUB8Ghnrj4WqhsYuj1FbD+mfUc6OOui3Py5iSGTSo/rGlufjmRVfNxcVdhS0cmbgtz1f9BPqNhz99Y3/7tmb1XWoSdUwb9DSwJanQmLBY2P2xErkbNE0JmJnjHaLCPm3p4QfYmV8yl0l2yOCnqwwdd5/6u93sPX8Mgy+lrBFC3Al8CzgD70gp9wkhFgHJUsrVQogxKMMeCFwohHhaSjkSuBKYAgQLIRYYDrlAStlKOcYuRugwVV5fegbSf4Iv71Q/7us/V1/arkD/SUpqufBE/WRlVan6MU7+S9PtTb1GU+Cs6U1WCyGICPQiItCLc0cYMhzeyqPOKZ6vZ03mUE4Jv+cUU7wnnL5lB3jlh0MNGhb1DfBkUKi6AAwI8WJAiA/9Q7wI9/fEydocgZFdy5URj75C3Wrnplov6IF6b94Yvw8aqLzIujrHmt2YU9wKgw9qYnLvZ/Ynbo9uUl5lxmaV3tmRxkJKVTew9zMltWGcdG4pptqMEfa3NToaZWcaZucYcXI25OK3oYcfMcb2Ns3Nxc9Pr79I+PdTxZldJI7vUAxfSrkWWNto2RNmz7ejQj2N91sOLG/lGLs+xonbz29VKX8RY2DuspZngbQH5nH8WEN2UtYukLWWv/CeAaoM3tGKW4OkglPCdYwM92dkuCHF0zkOtvzMgadnkp5XTtqpUtJyS0x/247mUW7WAMbNxYn+wV70D/amf4g3kUFeRAV7ERXkTXiABy61FbDxeeX5Tn1IGfyTv9k2+HlpSi7C+MMNHqgyp0ospFfaw2TwW1gtG54AO95VcwjGC5AljmxQYntOLvDLv+CqD1p2vpZQnq9056mAza/CjCdbd7ycvaohjzHt1xZ9YtSjkwsMnW15G98+bePh19ZA6WnrKZlG/PqquwpHZZILjtXXtAih7my6SC6+llZoC4wG//D3MPY2mLlYSSJ0JXpHq05D6T/XG3yjQmZfKxXAYXGO543nH1WVpI0Nb0A/qKvGoyKXYX3CGdanYdMWKSWniis5kltK+plSjp5Wf0dOl7Lx91yqaurn+F2cBPd7f8PC6mxe7/04cq/kFmcv8g9tx2nwFQR5uyEsec1n0tTEufF/EjRQPeYdaZnB9wxSd3AtwXzi1pbBT9sA/Serz/OnF1XxVsjglp2zuRgrWb1C4Nc3YML/a9ocpDlY0sC3hoe/+j0FRFm/s/ANUxPjraXsNCDtx/CdXVQCgiOpmbU1KknBKEII6re3a3nL7ijbmB5t8BcsWNA2B/KPUBOzfUc7NhnXGTg5Q9TEhnH8zGQVz/YOtrxPeLyaPHOkibOlXqVQ71UXHLdoXIUQ9PbzoLefBxMGNRxHXZ0kp7iCY2fKOHamlOycbG7Y+RnbXcfwxtFe5O87zBi3CGr3bmHezu/xcnOmb4AnfQM96RvgSUSgF30DPZme/TsufgNwqZMqpdRo8M+k1d/5OIqloqvm0Gu40tHP2tVQM8ac/HR1AR23UOkb/fIKbH4FLnq15edtDkaDf+4iWH2navIyc3HLjlVdoS5Ww5rRd+K6z1TGjjV8w1TIq7XYK7oyJ7C/YyGdohPqrtk8C6tPtHKG8o/avsh3ANrgtwVCtPwH0ZH0n6R6zRYZpAEyt6viIGuYt54baKV/rpGcfSps0qtRnNa8+CqyeeUbTk6CMH9Pwvw9GT8wGNa9BXWljLn1ZXb1iaaksoaq1d/il7qSx88dxomCSk4UlJGZX07K8QIKyqoByR73ND7KCWfxY9/Q28+DCH9XPsCF7cnJ/F45lT7+HoT5e9DH34MQb3fbcwhFWa0TP3N2VQbAVqjM2Ct30DRVzJdwrfIQpz3aMcJrRjXKAVMgdh5sewsm/Nm6oqQtTqcqA2gr5NYYe6FQvzClUVNVqiq6W4o9HR1zAqNUz2d7mFIyzSaCjZPVxir0TqRHG/zTp08DEBLSitvVPxLmcfx+Yw0TVjYE3UyVjyn2DX72Xgg+q6lnFmAw+K3tfFWUpcILsVeajIePuwsMSoR9y7hpBBDc8GJTUllDTlYGvsvKiY5J4LaAgWQVVJBVUM4J0ZvCEwd58ui+Bvu4Ogt6+Srj38dw59Hbz50+/up5YuFJZK9oHMjGtk54Avz2ifVb/LQNKv0wxJCOOeFOSH4Xfn0dZjzVmjM7RlEWINTFZcoDsPtT+PmfMPv55h/LGBI0JgG0BcYJ8+Ls1hlQoyaPvZAOKEdm53uGgi8bE/bGsI+5h99ruHKGsvfWq352Ej3a4F9xhbqlblUe/h+JPrGqSMw8rGMrQ8ErSE22OTJxm7PH8lyAm7fS5G9t8dXG55U2+rS/NlxurvvT6Mfv4+6Cj5O6bU8cNYbEwWZdwj6IJrIok+3XzSC7sIKTheVkF1VwsrCCbMPfgZNFbEg9RVmVmlR2ppbf3U/x2o5Sluz6H739PAj1daeXnwe9fN3p5etOqOGvl69aF+Dp2vSOITxB9RrOS2sal6+rhaM/wtDz67N4ggfBiItg+ztw9l+sN69vK4qylNfr7KrOHX81JL+jwpbNTUQ4sUPVgBjDaG2ByeCfbJ3Bb05Ix9z5sWXw89PVBK+fWQ6LqycED+4SmTo92uD3OJycVROW9J+V5KyLp/3JtLBY+wa/olCFbEZZEYhrrUzy6UMqpDHm5qYVqqHDVYFO9m6IvqzpvsbJvcaGIXgQIn0ToT5uhPq6ExNhPYOkuKKanKIK8k4ew/lzSfyI4Vzp249TRZWcKq7gt+MFnCquoKK6qQKni5Mg2MeNEB93QnzUxWC4UwA3Acmb11MxPJBgHzeCfdwI8nLDJWe3mjMZmNTwQJPuVlIDO5e1fyFfUVZDozblQfjtY/j5/2DOS8071omdhu5nbagdZC6v0BpKTqlJYlcP+9v2iQGEkg23lj0EKs7vH9FU2K1PNBzf3qrhtgXa4Pc0+k9SFap11dB3lP1S8fB4OLhG6eRY8yxzDDp6xpS6xgRE2peQtsUPz6gL1JQHmq5zcVPN5E/utrxvXppK8fNvpNEUNBCqywyTsLZz6n09XPH1cIUq1bFq6uhYpg5teKGUUlJaVUtucaXp71RxBbnFlZwuqeR0SRWnSyo5lFPMmhLJNS5upGzbyOLN9UXsQsC9Hl9zF3DLzz647N6hLgTe7gR7BzMnZBy+P/+bI5FXE+znQ4CXG24u7ZD1UZTV8AIZGAUJ18OOZUqmO6Cf9X3NqSpVhXFD21gpsq3kFUpsNC9vjLuPQfbBTgmReQ6+Ob2jVeV0eUHTYrIORBv8noYxjp+f7lg80Xgrm71HXSwsYd70xBIBkXBoXfNVIqvK4JsHlGc79WHLaqSgwjqp/7N8/DNp6q6gscdlrJjMO+J4JyMbOfhCCBVCcndhQIjtiUQpJbVvxnGtzCN25gTOlFRyurSKMyWVzEk5SGbVQAqdA8k7VcLWI5UUlFcjJfzglMQyt7/z1n9eYFWd0urxdXchyMeNAC83grxcCfR2I9DLjSBvNwK8XAn0qn80Prdb5Vyc1TR7acr9Sur3pxfgolcc+7xO7lbFY8ac9LbC3U8pb7ba4NvoZWuJ8HiljWWLgmOW7wCMzlDO3uZnhrUh2uD3NPrEKe2cqmL7FYZgFiNPsW7ws/eoOK21dEX/fqrQqfS0daPdmNxU+HS+ujOYfL9l795InzgV8im2oKCYd6ReUsEc47K8I9bfV2NMBr91QnhCCFwiRuGyazljo/zrZS2qymDLHhh7K5+eV9//uKa2jvyyavJLJ1P68Rc8VbeehAkLySurIa+0irzSKvLLqsgtqeT3nBLyy6pM8w6W8HR1JtDLFX8vNwK9XAnwcsXfUz0PdqvhpopCUsv9yEs7Y1jnir9nH7xGzUfseBfOvtcxiYETO9Rj31Gt+biaIkTbFF+V5Cgj7ihh8aoCuTjbcrZUZQmU5lqWajA6Q9na4Hcat99+e2cPoeNxdlFSzYfXOWbwfXopA2crjp+zz9D83Ir3bpJJznDM4KcYmoe4esF1qyxKOzQgzCgDsbuhwZfSugyyfz8V6mlOQ/Pik2pCrjVFSEbC42HbG2p+opdhMjljC9RWNUmVdXF2Mk0GM+0++PxWrgv+HSZY71VQUV1LQVk1+WXqYmB8XlBWTX5pFQXl1RQYXv+eU2J63k9mcZM7vL6zjM+TG2odacDVJQAAHbpJREFU9nWOZ72rExv+cx9Lgu4zXAhc8fNwrX/u6WJaNuLwVnx8IyhxCsDXWAPRVviFt00Mv7kePhji+LOarjfOU1lSQvXto5IXcjpXG79HG/x587pfx0WHGHOz+sE4mtMdFmfd4NfVqjittQlbqI/5Fhy3fXtvDOHsWq7E3i5/27FwS+9oQKiJW/MfYvFJFae3lCHi7KJ+mM2RSTZ6do2F5lqCseL2ZEq9wT+yQalpRk20vl/0ZWpOY/OrDRuDNMLD1Zk+/s708XdgQtKAlJKK3zfAR3DXZVOZGziGwrJqCsqrKTT87fn9PKblf88nbn8hr7SKo6dLKSyvpqi8GvOumQA/uW3lZzmQOxetA1TWlJ+HC74e6sLg6+Ha5LWvR/2jn9lzH3cXvN1c6jOefPu0rj9wVam6y3UkJdNIn1hAqP+ZJYNvLovcmC4isdCjDf7x4ypVsF8/ByehugtDZ1n+wlojPF5N9FoqdMk7qoyqrcIa8+Ira5SeUZ2ijCGcpEccb2Hn7qNqABpflBqLpjUmaGDzDH5ri67MCRmi7mCydtVLVqdtVJ3R3Gz0Z3B2hdELlNHPO9Km6Y5CCDwrVDHSgAFDGBBs4U4mah588hVLZ9AgNCGlpKSyhqKKGgrLqinNzyZyRS4FI+bzRPgIiiqqKSqvoaiimmLD81PFFRw+VaNeV9Q07LNscXzg4+aCr4cL98gaLqnK4pa3f8XH0xVfw/yJj+Ei4ePujI+7Kz6Gi4XpomHYzrk5RVdGjBO3WVYmbi3l4JvTJwa2vdl27RlbQI82+Ndffz3Qg/LwW0pYnJp8O7YFBs9ouM6apII5ngFKx8dW56s9n6pG19d8atNztT7G2KZpb41lkRsTNBDSf3F8Mrm1hT7mODk31CoqyVW3++c8bns/gLirYcOzSkL4nMfaZjxGigytLqypgQ6YokJhh9c3MPhCCFM2U98ATyhR4aDYsecQ299+vF9KSXl1LcUV9RcA4/MSs+fFleo52WG4nammtvQMBwq8TNuYC/HZYoLrYT5yhge+yebApk14u6kLgZe7Cz7uzni71V8cvN1d8DYsi/MdTkDmFjJyS9T2bs71dx756Woy2cuKVEnvaNU858zh+ru6DqZHG3yNg0SMUQb7g8vVD37UfKWN4uqhDL5wVvnwtgjoZ7v4KnWtOkZLjD2o2+29q1SzE6Po1pk0FSLxbyLkqggapDROSk45JhtgKXulNYTFw46lyuM7asj+cERO27+vav2X8qG6E2qLEJORoiw1AW/tLsPDHyLGQtp62yqaJ3ao6lLjpL8dhBB4ubng5eZCbz8HwlD7jsCK/7D8ysgGd5c1tXWUVNZQUllDaWUtJZXVFFeo58UV1ablEdnpcBiC+0TSS3hQUllDdlEFpZU1lFTWUlpp+eLxJ2dfnnDN4aqXviSX+vRKT1dn3nBJpi8h3PnKz3i7OePl7qIe3dSFoX+NNzcBP/y4npwoL7wM68IDPOrVZdsZbfA19vEOgf+3VaXl7XwfVt2kjELcVaqwJmSw/eIVW8VX5fnK025Bk3YTxonb7D31MhBn0pTEszWDaFLNTLNv8KvKVIGZoymcjhCeAL/+F07/ruQUPALq02DtkXAtrFig4v5nzbC7ucM40iv2rHPgh8WGSU8rMfATO9QFvL3a+pnLK5gZfBdnJwK8VJqqTbZthMPw8NypVhMJauskpVU1lFXWUlJZQ1lVDSLDCb57n1eTBEcCYyitrFHbVNUydPcZclwiiAj0pKyqhqLyarILyymtrKW0qobqKsH1zs6k/raFv+/obzrPrJF9eP36Nk5dtYI2+BrH8AtXqZFn36eavOx8T5Xb11apRiT2CIhUGj6WOLxeCWy1pkCnj5nEgtHg56Wp2L41zHPxbU2UQusbn1jCJJW8UxnugVMd99aHnq8UTHctb2ODf8K+wR80XRn8tA0QZyHxQUpl8IfNabtxNcZk8LNatn9JjroztdHcxdlJ4OehMo5MhE6G7wTjPY4zfpxZMZ+UsCuH3vGzeXOWdX0q+Z9h3OxdxiWXnENZVS1llbV4urXhHZodtMHXNA8nJ1X2PzBJhU8OrFYZNfYI6AeVRZYrDVPXgndo6wp0vIOVfkm2oeK2rk5NKA8+18aYogypmQ5M3BqFttpSrTL4LNXkZM8KZWgH3u/4vi7uSsky+Z2GYazWUnzSesW0kbB4FadOW2/Z4OenQ3le2xdcmWP8P7Q0NbMkR33nmhsOc/dRd7SNK25Lc1Xygp3m9KJPDK5HNhLmb0P+uR1xqC5bCDFLCJEqhDgshHjYwvopQoidQogaIcQVjdbNF0IcMvzZyN3reO677z7uu6+T+83+kfEKUhkjjjTmsJapU1MFh75Xzapb2xwiLLZeYqEoU02QWZuwBZUpERDpWEPzNiq6aoCTU33TbrAtVW2JhOvUHdaeNmoXXVutwjR+dgTSnJxg0DmQ9oO6sDbGVHDVjgbfxV1ddFpafGUrHGWPsPimjYHMG5fbok+06rRWerpl524ldn9hQghn4DVgNjACuFoI0bg5ZQawAPiw0b5BwJPAOGAs8KQQwk4njY7jwgsv5MILL+zsYfQMTMVXjSZuMzYrbfO20FvpEwtnDqn0UXspmUYcTc1sbWtDaxjDOoH9HateNadPjLpg7Hq/bcZSnA1Ixxq8DJquvNpsCxpGJ3YqYb5edibyW4tveCsMfjN0dBoTHq/OW5xTv8xWDr45porbzinAcsSlGgscllIekVJWAR8DDURYpJTpUsrdQOPL/XnAOillnpQyH1gHNCMBvH1JTU0lNTW1s4fRMzAa/MaZOgfXKmG0gUmtP4cxfTRnn/2UTCNBg1ToR9rOAac4W+XNO9KXtTkYDf7ApJbtn3C9MrqO9h62hbHTlUMG/xz1mLa+6boTO9T/wp4wX2tpjbxCc6tszTGXSjZSkK4eAyKbbN4Ac02dTsARg98XMP+VZhqWOYJD+wohbhVCJAshknNzcx08dOu57bbbuO222zrsfD0ar2Dl9ZmHdKSE1G9UKMNWsZGjmCQWfoMzR9T57E2yBg1UFZeldr53xqKrtpT5BSVz4eLZ8sYY0Zerlom72qDJeXEzDL5vb+gdA4d/aLi8tlp9/u0ZzjHiF9ayGH5lsfLwW3q3FmaouDUvwMpPVxcQe99j7xAVMstMbtm5W0nndtQ1IOX/b+/cg+Sqqzz+OZMweQcySYCEAJnhESGCEx0CmzURMSDiQlxFEzArKD7QpbaybtUuD9eyUCnd1ZXdUiRsCdGgksiutXENUgbX1RUCeRIMMTAkBPJAkhBIQgjkcfaP87vJpe2Zvt19+97u6fOpmpo79/nt371z7q/P7/zO0btUtUtVu0aPTphcy2ksRKz3E6989ce19ndv+cXLYfhJVmD8hTXWw2/rKD0uEC9o3ht7XkjXfx9x7Di46fmjPeZyGdxm0TBrFlj92GqIevhJI5FOvwieX2oGNOLFdZYoL+2EacUYNsZe1IcOlHfc2p/C4YOVz/kYMCzM7I4b/E2l/fcR46daTYpi4x81JonB3wLEcw+MC+uSUM2xTl+jcPLV+gfs95kpeflEjg7c7nwGRiZIOxD5+EsN3O5JMa1CIdW6PibNhv0vW7RTNezeat82ShWsjzh9uhnOjb85ui6LAduIYWMAPVq5KimrfmgpEpIkD+yJsZ0FPfxNpf33Ee1TYd8Om1meMUkM/jLgDBFpF5FWYBawKOH5HwQuEZERYbD2krDOaUYKJ1+tX2xlESspjt0TY95mydx2PVvafw/2rUP69d7DV+05JW490HGhhaSuure680SVrpK6rU6+wFIJdMf8+FtW2LespMavGuKTr5Ky42n7VjJpdnXuubGTrBOw90X7hrF7c/HCJ8WIsrfGX5QZUdLgq+pB4AbMUK8DFqrqWhG5VUSuABCR80RkM/BhYK6IrA3HvgR8GXtpLANuDeucZuS4Uyw++41Xzfe6dWV67pyIE8+1UMXDB5Llvel3jOnqzeC/tgsO7k/m286Dln7QebWFSb6yufLz7N5aOiQzTv9W663GB25rUdKwJ6JZz7vLmHy16l57wZ87q7prRwO3W1db5JkeTv6SO+4U23fjb6vTUAGJJl6p6mJgccG6L8aWl2HummLH3g3cXYXGmvGFL6SceMrpnXikznOP2HLaszHjuVuS9PAhhGb24tKpxaSrtOm8Gn7zT1ZL4F29FIvpjT1breZxOZz2HnjqF+YSG3qCuSnO+ovKrl8u5fbwDx2Ex39svvtqv1WOiaVKjlxySX34YL38tf9l6cXTzIVUgroYtM2L6dOnM316itPSnd6JT75a/4D1ckannDWw7TRzM0DvaRXedExH76GZUfRKLQZt06Kt3QYDV1cYrXP4sH3rKjd1RFScpvshi86pRUnDnhg8ymZKJ02v0L3E/P2TZld/7Wjgduvq5DH4ccZPs/knaYTTlkFTG/zVq1ezenWJosROekQ9/O1/sNmlEy5L/6t/S4vFOrcOTT6TcuRplvZh387i2xuhhw+WRmLXRnNBlcu+neYGK8elA/ayHDHe3DrRgO3YDCJ0wO710BOT9/BXzbd0Cmdcks71x3ZaD//lTdByTHkuv/ap9vvZbN06TW3w58yZw5w5c/KW0TwMPcHSFa+ab2kP0vbfR7zjWpj8qeQvkyg0c8dTxbfXInFaLThSp3dj+cdGefDLzQYqYm6djb+F5x+1l3rSusVpMHxMMh/+qzvM9XTuzPQmhI3ptHbbvNwi0MpxzQw7EUZNyHzgtqkNvpMxLS3Wg9zxlM1YLddfnJTOq2D6l5Lvf1KXTV7qKSfN7m0WqlgqBXTeJJ1TUIxyZtkWcvp7rK7A+sXZuXMiks62XbPAQkjTcOdERDVuN/2uPP99RPtUKypU7jyCKnCD72RL5NY545LaT71PypCRcM6HbUCvmDtkzwv137uHoz7kSgx+NeMU46eaLz1L/33EsLGlXTqqVsfhpK508/ucGGZ2lxOhE6d9mr0ot6xMT1MJ3OA72RIVNK+VO6dSzv+MpbddWSQR2Z6tjWHwWwebAay0hy/9KssgOXC41eKF7A3+8DE2/vL63p732brSoofS7N2Dfe4oMCBpDH6c8cGPn6Fbxw2+ky3HT7QB1TSLdqTBmHMtr/9j/26hcnEapYcP5Rdmj4gidCoNETzrCnPTJSxpmBpJQjNX3WsziN/6wfSvH8XjV9LDH9xm+Yii8pYZ0NQG/7bbbuO2227LW0ZzMflT8Der0s86mQbnX2+5feIpCg4dtFC+NEsb1pKRlRr8LdV9xsmfhr9dC61DKj9HJZSqfHXgNXjiP+DsK2rzzEV+/Ep8+GBunecfqz4PUkKa2uBPmTKFKVNKlLZz0qXfMZUXnqg1Ey6zuQJL7zy67tXt5qOt95DMiLYO07x/d3nH7UlQy7Y3WlosNj1rSvXw1/23xbun7c6JeNtVVki+0m827dMsYm3zY+nq6oGmNvgPP/wwDz/8cN4ynHqhX3/7BrLp/44WqGiESVdxokidXWWEZqrCK1sa5zPGKZVeYdV8632f+s7aXH/IKLjwxspdYadOsbGTjPz4TW3wb775Zm6++ea8ZTj1xNs/ZoVOHg29/EaZdBURGfwkZRsjXt9t0SL1miuoNwYMszGhYj38px40//ik2dWXz6wVA4ebWyijvDp12gqOkxODRtjknDU/sck6UYx3oxjDEaFMYjl+/N0N9hkLGTbmzT78w4fgoS/Djz5ig6Jd1+WnLQnt02DL8t4jjVLCDb7jFHL+9eZXXTHPjKG02JT8RmDAUJvRXM5s2yOzbBvV4MfSK+zdDvM/AL/9hpV//OQvbZ5FPdM+zSaFPbe05pdyg+84hRz/Fiu7uOx7lvp26AmZZjSsmrbTyuvhN0rqiJ4YPtZezJsegblTLeplxndgxrfhmEF5qyvNyRdYLp5na+/Hd4PvOMW44LPmJnhyUeMZwnJj8cstbVhvDDvRCpDMe78Z+E8uqV1UTi1oHWzVtzIYuE2UD7+vcvvtt+ctwalXTr/4qOFsNEPY1g57X7BCM0ni4ndvtVTD9Z4rqCdGjLfQ2bMut559Pc7xKEX7NKtn8NrLMOi4ml2mqXv4nZ2ddHZ25i3DqUdaWmDyZ2y5USZdRRxJopbQjx+VNmxUOmfDtYvhI/Mb09iDJVLTw7CptmHiiQy+iFwqIutFpFtEbiyyfYCILAjbHxWR8WH9MSLyfRF5QkTWichN6cqvjiVLlrBkyZK8ZTj1SufVFps+psE6BUcMfsLQzHJLG9Yb/Vth/J9nU1axVow7D/oPrLlbp6RLR0T6Ad8BLgY2A8tEZJGqPhnb7Tpgl6qeLiKzgK8DM7EatwNU9RwRGQw8KSI/VtVn0/4glfCVr3wFwKteOcUZONzSBdRrDHdPlJsmec9WGNdVOz1OafoPgFMuyN/gA5OBblXdACAi9wEzgLjBnwF8KSzfD3xbRARQYIiI9AcGAW8AZc75dpwcaTRjD/aiGjI6mcE/sN+qXTVyD7+vcOnXau6SSvI0nwQ8H/t7c1hXdB9VPQi8AozEjP+rwDbgOeAbqvpS4QVE5NMislxElm/fvr3sD+E4TgFRnd5SHJlY1sA+/L7C8WfVfC5Erbsvk4FDwFigHfg7Eeko3ElV71LVLlXtGj26QSa4OE49kzQ0s5pKV07DkcTgbwFOjv09Lqwruk9w3xwL7ASuBn6hqgdU9UXgd4A7Cx2n1rR12AzaA6/1vt+RHr67dJqBJAZ/GXCGiLSLSCswC1hUsM8i4JqwfCXwK1VVzI1zEYCIDAEuAP6QhvA0mDt3LnPnzs1bhuOkz5Gsmc/2vl+UVqHR5ho4FVFy0FZVD4rIDcCDQD/gblVdKyK3AstVdRHwPWC+iHQDL2EvBbDonntEZC0gwD2quqYWH6QSJkyYkLcEx6kNbSGJ2s5neq/junsrtA6zgV6nz5Nopq2qLgYWF6z7Ymx5PxaCWXjc3mLr64Wf/exnAFx++eU5K3GclEkamtnok66csmjq1Arf/OY3ATf4Th9k0AgY1JbQ4PuAbbPQgEHGjuMkIkmkzp5tjVnpyqkIN/iO01cpFYt/+JDlkfceftPgBt9x+iptHZbP/+DrxbfveBr0EBw7LltdTm64wXecvkpbB6Cwa1Px7avmQ0t/mHBZprKc/GjqQdv58+fnLcFxakc8a+boM9+87cB+WP1DeMv7YdgJ2WtzcqGpDf7JJ59ceifHaVRGnma/iw3crlsEr+2Cd3w8W01OrjS1S2fBggUsWLAgbxmOUxsGjbDsi8UM/vK77RtA+7uy1+XkRlP38L/73e8CMHPmzJyVOE4NECkemvniOnjuEbj41sZM/+xUjN9tx+nLFDP4y++Bfq3Q+dF8NDm54QbfcfoybR3w8nNw8A37+4198Ph9cNYVMGRUvtqczHGD7zh9mbYOK479SqhhtPY/4fVXoOsT+epycsENvuP0ZdpCpM7OUNB8+T0wagKcOiU/TU5uNPWg7f3335+3BMepLfGsmdvWwJblVjtVJF9dTi40tcEfNcp9mE4fZ8goy3f/0gbYsR76D4S3zSp9nNMnaWqDP2/ePACuvfbaXHU4Ts0QsWIoL6yBF56AiX9p8flOU9LUPvx58+YdMfqO02dp67C4+zf2+mBtk5PI4IvIpSKyXkS6ReTGItsHiMiCsP1RERkf23auiDwiImtF5AkRGZiefMdxShL58Y+fCOPOy1eLkyslDb6I9MNq074POBu4SkTOLtjtOmCXqp4OfAv4eji2P3AvcL2qTgQuBA6kpt5xnNJEOXW6Pu6DtU1Okh7+ZKBbVTeo6hvAfcCMgn1mAN8Py/cD7xERAS4B1qjq4wCqulNVD6Uj3XGcRJz5PvizG6Dz6ryVODmTxOCfBDwf+3tzWFd0H1U9CLwCjATOBFREHhSRlSLy98UuICKfFpHlIrJ8+/bt5X4Gx3F6Y8hIeO9XoXVI3kqcnKl1lE5/4J3AecA+4CERWaGqD8V3UtW7gLsAurq6tMaajrB48eKsLuU4jpM7SXr4W4B44vhxYV3RfYLf/lhgJ/Zt4DequkNV9wGLgbdXKzotBg8ezODBg/OW4TiOkwlJDP4y4AwRaReRVmAWsKhgn0XANWH5SuBXqqrAg8A5IjI4vAjeBTyZjvTqueOOO7jjjjvyluE4jpMJJQ1+8MnfgBnvdcBCVV0rIreKyBVht+8BI0WkG/g8cGM4dhfwL9hLYzWwUlV/nv7HqIyFCxeycOHCvGU4juNkQiIfvqouxtwx8XVfjC3vBz7cw7H3YqGZjuM4To409Uxbx3GcZsINvuM4TpPgBt9xHKdJEAumqR9EZDuwqczDRgE7aiCnWupVF7i2SnFt5VOvuqBvaTtVVUf3tkPdGfxKEJHlqtqVt45C6lUXuLZKcW3lU6+6oPm0uUvHcRynSXCD7ziO0yT0FYN/V94CeqBedYFrqxTXVj71qguaTFuf8OE7juM4pekrPXzHcRynBG7wHcdxmoS6M/hV1s+9KaxfLyLvTXrOWmsTkYtFZEWo6btCRC6KHfPrcM7V4ef4jLWNF5HXYte/M3bMO4LmbhH5t1DFLCtdH41pWi0ih0WkM2zLqs2mhcI9B0XkyoJt14jI0+Hnmtj6qtusGm0i0ilHa0ivEZGZsW3zRGRjrN06s9QWth2KXX9RbH17uP/d4XlozVKbiLy74HnbLyIfCNuyarfPi8iT4b49JCKnxral87ypat38AP2AZ4AOoBV4HDi7YJ/PAXeG5VnAgrB8dth/ANAeztMvyTkz0DYJGBuW3wpsiR3za6Arx3YbD/y+h/M+BlwACPAA8L6sdBXscw7wTA5tNh44F/gBcGVsfRuwIfweEZZHpNFmKWg7EzgjLI8FtgHHhb/nxffNut3Ctr09nHchMCss3wl8NmttBff3JWBwxu327tg1P8vR/9HUnrd66+FXUz93BnCfqr6uqhuB7nC+JOesqTZVXaWqW8P6tcAgERlQgYbUtfV0QhEZAwxX1aVqT9YPgA/kpOuqcGyalNSmqs+q6hrgcMGx7wV+qaovqaUA/yVwaUptVpU2VX1KVZ8Oy1uBF4FeZ19mpa0nwv2+CLv/YM9Dpu1WwJXAA2pFm9Iiibb/iV1zKVZsClJ83urN4FdTP7enY5Ocs9ba4nwIqwvwemzdPeGr4j9W6AKoVlu7iKwSkf8Vkamx/TeXOGetdUXMBH5csC6LNiv32DTarFptRxCRyVhv8pnY6q8Gl8G3Kux0VKttoFj96qWRywS73y+H+1/JOdPSFjGLP33esm6367Aee2/Hlv281ZvB79OIyETg68BnYqs/qqrnAFPDz19lLGsbcIqqTsKK1/xIRIZnrKFHROR8YJ+q/j62Ou82q3tC728+8HFVjXqzNwFvwWpMtwH/kIO0U9XSBVwN3C4ip+WgoUdCu52DFXyKyLTdRGQ20AX8c9rnrjeDX0393J6OTXLOWmtDRMYBPwU+pqpHelyquiX83gP8CPvql5m24ALbGTSswHqDZ4b9x8WOr6TdqmqzwJ/0tjJss3KPTaPNqtVGeGH/HLhFVZdG61V1mxqvA/eQfbvF790GbCxmEna/jwv3v+xzpqUt8BHgp6p6IKY5s3YTkenALcAVMS9Aes9bNQMRaf9gFbg2YIOu0cDGxIJ9/po3D/ItDMsTefOg7QZsoKTkOTPQdlzY/4NFzjkqLB+D+TCvz1jbaKBfWO4ID0ybFh8QuiwrXeHvlqCnI482i+07jz8dtN2IDaCNCMuptFkK2lqBh4A5RfYdE34LcDvwtYy1jQAGhOVRwNOEgUvgJ7x50PZzWWqLrV8KvDuPdsNefs8QBt1r8byVJTqLH+Ay4KnwwW8J627F3ngAA8PD0R0+bNwY3BKOW09stLrYObPUBnwBeBWr6xv9HA8MAVYAa7DB3H8lGN8MtX0oXHs1sBK4PHbOLuD34ZzfJszMzvB+XggsLThflm12HuYXfRXrha6NHfuJoLkbc5uk1mbVaANmAwcKnrXOsO1XwBNB373A0Iy1TQnXfzz8vi52zo5w/7vD8zAgh3s6HutgtBScM6t2WwL8MXbfFqX9vHlqBcdxnCah3nz4juM4To1wg+84jtMkuMF3HMdpEtzgO47jNAlu8B3HcZoEN/iO4zhNght8x3GcJuH/Aaz3vWbpFuPEAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "A_nume = obs.data()\n", + "t_nume = [(A_nume[t][0] + A_nume[t+1][0]) * 0.5 for t in range(len(A_nume)-1)]\n", + "k_nume = [(A_nume[t][1] - A_nume[t+1][1]) / delta_t / (0.5 * (A_nume[t][1] + A_nume[t+1][1])) / concB / ka for t in range(len(A_nume)-1)]\n", + "\n", + "# since the theoretical formula is a long time asymptotic, it becomes irregular at t <= R^2/D.\n", + "# Therefore it skips t[0] here.\n", + "k_theo = [calc_k(ka, D, R, t) / ka for t in ts[1:-1]]\n", + "\n", + "plt.axvline(x = 100 * R**2 / D, color='k', linestyle='--', label=\"relative error < 1% in asymptotic\")\n", + "plt.plot(ts[1:-1], k_theo, label=\"SCK theory\")\n", + "plt.plot(t_nume, k_nume, label=\"SGFRD simulation\")\n", + "plt.legend()\n" + ] + } + ], + "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.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}