CSD_2D_model_v3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Some Packages\n",
    "import numpy as np\n",
    "import imageio\n",
    "import scipy as sp\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import os, glob, ntpath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Given Some Coefficients in Partial(K)/Partial(t)\n",
    "D = 0.75\n",
    "\n",
    "# Given Some Coefficients in F(K)\n",
    "A = -0.3\n",
    "Kr = 0.03\n",
    "Ktheta = 0.2\n",
    "Km = 1.0\n",
    "r = 8\n",
    "\n",
    "# Given Some Coefficients in Partial(r)/Partial(t)\n",
    "B = 0.0001\n",
    "C = 10.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0025000000000000005"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Timesteps\n",
    "maxIter = 2000\n",
    "\n",
    "# Given the range of the rectangle\n",
    "lenX = 200\n",
    "lenY = 200\n",
    "\n",
    "deltaX = 0.1\n",
    "deltaY = 0.1\n",
    "deltaT = (deltaX*deltaY)**2/(deltaX**2 + deltaY**2)/2/D\n",
    "D*deltaT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Given Initial values\n",
    "k_init = np.full([lenX, lenY], Kr, dtype=np.float64)\n",
    "kt_init = np.zeros([lenX, lenY], dtype=np.float64)\n",
    "r_init = np.zeros([lenX, lenY], dtype=np.float64)\n",
    "rt_init = np.zeros([lenX, lenY], dtype=np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   0.,  100.,    0.],\n",
       "       [ 100., -400.,  100.],\n",
       "       [   0.,  100.,    0.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Definie the filter of the laplacian operator\n",
    "filter_laplacian = np.array([[1/deltaX/deltaY, 1/deltaY**2, 1/deltaX/deltaY],\n",
    "                           [1/deltaX**2, -4/deltaX/deltaY - 2/deltaX**2 - 2/deltaY**2, 1/deltaX**2],\n",
    "                           [1/deltaX/deltaY, 1/deltaY**2, 1/deltaX/deltaY]])\n",
    "\n",
    "filter_laplacian = np.array([[0, 1/deltaY**2, 0],\n",
    "                           [1/deltaX**2, - 2/deltaX**2 - 2/deltaY**2, 1/deltaX**2],\n",
    "                           [0, 1/deltaY**2, 0]])\n",
    "\n",
    "filter_laplacian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Dirichlet_bd( maxIter, clamping_ticks, K, Kt, R, Rt, top, bottom, left, right):\n",
    "    \n",
    "    center_point_loc = np.zeros((2,1))\n",
    "    center_point_loc[0,0] = int(lenX/2)\n",
    "    center_point_loc[1,0] = int(lenY/2)\n",
    "    center_point_value = np.zeros([maxIter, 2])\n",
    "    R_value = np.zeros([maxIter, 2])\n",
    "    \n",
    "    # Fixed the boundary values\n",
    "    K[:,0] = left\n",
    "    K[:,-1] = right\n",
    "    K[0,:] = top\n",
    "    K[-1,:] = bottom\n",
    "    \n",
    "    for i in range(maxIter):\n",
    "        \n",
    "        if i == clamping_ticks:  \n",
    "            r2 = r**2\n",
    "            for a in range(lenX):\n",
    "                for b in range(lenY):\n",
    "                    dis2 = (a-center_point_loc[0,0])**2 + (b-center_point_loc[1,0])**2\n",
    "                    if dis2 < r2:\n",
    "                        K[a,b] = Km\n",
    "            \n",
    "\n",
    "        # Fick's second law\n",
    "        R = R + deltaT*Rt\n",
    "        Kt = D*sp.signal.convolve2d(K, filter_laplacian, mode = 'same', boundary='symm') + A*(K - Kr)*(K - Ktheta)*(K - Km)*(K + 0.1) - R*K\n",
    "        K = K + deltaT*Kt\n",
    "\n",
    "        # Fixed the boundary values\n",
    "        K[:,0] = left\n",
    "        K[:,-1] = right\n",
    "        K[0,:] = top\n",
    "        K[-1,:] = bottom\n",
    "\n",
    "        # plot the images\n",
    "        plt.imsave('./images/{}.png'.format(i), K, cmap = 'hot')\n",
    "#         plt.imsave('./images/{}.png'.format(i), K, cmap = 'hot', vmin = 0, vmax = 1)\n",
    "        \n",
    "#         im_K = K\n",
    "#         im_K[0,0] = 1\n",
    "#         im_K[lenX-1, lenY-1] = 0\n",
    "#         plt.imsave('./images/{}.png'.format(i), im_K, cmap = 'hot')\n",
    "        \n",
    "#         fig, ax = plt.subplots()\n",
    "#         cax = fig.add_axes()\n",
    "\n",
    "#         if i < clamping_ticks:\n",
    "            \n",
    "#             im_K = K\n",
    "#             im_K[0,0] = 1\n",
    "#             im_K[lenX-1, lenY-1] = 0\n",
    "#             im = ax.imshow(im_K, cmap='hot')\n",
    "#             fig.colorbar(im, cax=cax, orientation='vertical')\n",
    "#             plt.savefig('./images/{}.png'.format(i))\n",
    "            \n",
    "#         else:\n",
    "#             im = ax.imshow(K, cmap='hot')\n",
    "#             fig.colorbar(im, cax=cax, orientation='vertical')\n",
    "#             plt.savefig('./images/{}.png'.format(i))\n",
    "\n",
    "        # Fick's second law\n",
    "        Rt = B*((K-Kr) - C*R)\n",
    "        \n",
    "        center_point_value[i,0] = i\n",
    "        center_point_value[i,1] = K[int(center_point_loc[0,0]),int(center_point_loc[1,0])]\n",
    "        R_value[i,0] = i\n",
    "        R_value[i,1] = R[int(center_point_loc[0,0]),int(center_point_loc[1,0])]\n",
    "        \n",
    "    return center_point_value, R_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Neumann_bd( maxIter, clamping_ticks, K, Kt, R, Rt, top, bottom, left, right):\n",
    "    \n",
    "    center_point_loc = np.zeros((2,1))\n",
    "    center_point_loc[0,0] = int(lenX/2)\n",
    "    center_point_loc[1,0] = int(lenY/2)\n",
    "    center_point_value = np.zeros([maxIter, 2])\n",
    "    R_value = np.zeros([maxIter, 2])\n",
    "    \n",
    "    \n",
    "    for i in range(maxIter):\n",
    "\n",
    "        if i == clamping_ticks:  \n",
    "            r2 = r**2\n",
    "            for a in range(lenX):\n",
    "                for b in range(lenY):\n",
    "                    dis2 = (a-center_point_loc[0,0])**2 + (b-center_point_loc[1,0])**2\n",
    "                    if dis2 < r2:\n",
    "                        K[a,b] = Km\n",
    "        \n",
    "        # Fick's second law\n",
    "        R = R + deltaT*Rt\n",
    "\n",
    "        Kt = sp.signal.convolve2d(K, filter_laplacian, mode = 'same')\n",
    "        \n",
    "        # Fixed the differential boundary values\n",
    "        Kt[:,0] = left\n",
    "        Kt[:,-1] = right\n",
    "        Kt[0,:] = top\n",
    "        Kt[-1,:] = bottom\n",
    "        \n",
    "        Kt = D*Kt + A*(K - Kr)*(K - Ktheta)*(K - Km)*(K + 0.1) - R*K\n",
    "        K = K + deltaT*Kt\n",
    "        \n",
    "        # plot the images\n",
    "#         plt.imsave('./images/{}.png'.format(i), K, cmap = 'hot')\n",
    "        plt.imsave('./images/{}.png'.format(i), K, cmap='hot', vmin = 0, vmax = 1)\n",
    "        \n",
    "#         fig, ax = plt.subplots()\n",
    "#         cax = fig.add_axes()\n",
    "\n",
    "#         if i < clamping_ticks:\n",
    "            \n",
    "#             im_K = K\n",
    "#             im_K[0,0] = 1\n",
    "#             im_K[lenX-1, lenY-1] = 0\n",
    "#             im = ax.imshow(im_K, cmap='hot')\n",
    "#             fig.colorbar(im, cax=cax, orientation='vertical')\n",
    "#             plt.savefig('./images/{}.png'.format(i))\n",
    "            \n",
    "#         else:\n",
    "#             im = ax.imshow(K, cmap='hot')\n",
    "#             fig.colorbar(im, cax=cax, orientation='vertical')\n",
    "#             plt.savefig('./images/{}.png'.format(i))\n",
    "\n",
    "        # Fick's second law\n",
    "        Rt = B*((K-Kr) - C*R)\n",
    "\n",
    "        center_point_value[i,0] = i\n",
    "        center_point_value[i,1] = K[int(center_point_loc[0,0]),int(center_point_loc[1,0])]\n",
    "        R_value[i,0] = i\n",
    "        R_value[i,1] = R[int(center_point_loc[0,0]),int(center_point_loc[1,0])]\n",
    "        \n",
    "    return center_point_value, R_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Periodic_bd( maxIter, clamping_ticks, K, Kt, R, Rt):\n",
    "\n",
    "    center_point_loc = np.zeros((2,1))\n",
    "    center_point_loc[0,0] = int(lenX/2)\n",
    "    center_point_loc[1,0] = int(lenY/2)\n",
    "    center_point_value = np.zeros([maxIter, 2])\n",
    "    R_value = np.zeros([maxIter, 2])\n",
    "    \n",
    "    \n",
    "    for i in range(maxIter):\n",
    "\n",
    "        if i == clamping_ticks:  \n",
    "            r2 = r**2\n",
    "            for a in range(lenX):\n",
    "                for b in range(lenY):\n",
    "                    dis2 = (a-center_point_loc[0,0])**2 + (b-center_point_loc[1,0])**2\n",
    "                    if dis2 < r2:\n",
    "                        K[a,b] = Km\n",
    "        \n",
    "        # Fick's second law\n",
    "        R = R + deltaT*Rt\n",
    "        Kt = D*sp.signal.convolve2d(K, filter_laplacian, mode = 'same', boundary = 'wrap') + A*(K - Kr)*(K - Ktheta)*(K - Km)*(K + 0.1) - R*K\n",
    "        K = K + deltaT*Kt\n",
    "        \n",
    "        # plot the images\n",
    "#         plt.imsave('./images/{}.png'.format(i), K, cmap='hot')\n",
    "        plt.imsave('./images/{}.png'.format(i), K, cmap = 'hot', vmin = 0, vmax = 1)\n",
    "\n",
    "#         fig, ax = plt.subplots()\n",
    "#         cax = fig.add_axes()\n",
    "\n",
    "#         if i < clamping_ticks:\n",
    "            \n",
    "#             im_K = K\n",
    "#             im_K[0,0] = 1\n",
    "#             im_K[lenX-1, lenY-1] = 0\n",
    "#             im = ax.imshow(im_K, cmap='hot')\n",
    "#             fig.colorbar(im, cax=cax, orientation='vertical')\n",
    "#             plt.savefig('./images/{}.png'.format(i))\n",
    "            \n",
    "#         else:\n",
    "#             im = ax.imshow(K, cmap='hot')\n",
    "#             fig.colorbar(im, cax=cax, orientation='vertical')\n",
    "#             plt.savefig('./images/{}.png'.format(i))\n",
    "\n",
    "        # Fick's second law\n",
    "        Rt = B*((K-Kr) - C*R)\n",
    "\n",
    "        center_point_value[i,0] = i\n",
    "        center_point_value[i,1] = K[int(center_point_loc[0,0]),int(center_point_loc[1,0])]\n",
    "        R_value[i,0] = i\n",
    "        R_value[i,1] = R[int(center_point_loc[0,0]),int(center_point_loc[1,0])]\n",
    "        \n",
    "    return center_point_value, R_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_fick_law(boundary, clamping_ticks = 100, top = 0.0, bottom = 0.0, left = 0.0, right = 0.0):\n",
    "    '''N: The number of the random points you want to pick'''\n",
    "    '''clamping_ticks: '''\n",
    "    \n",
    "    # remove other images in the dir\n",
    "    for image in sorted(glob.glob('./images/*.png'), key=lambda k: int(ntpath.basename(k).replace('.png', ''))):\n",
    "        os.remove(image)\n",
    "    \n",
    "    # Fixed the initial values\n",
    "    K = k_init\n",
    "    Kt = kt_init\n",
    "    R = r_init\n",
    "    Rt = rt_init\n",
    "    \n",
    "    if boundary == 'Diri':\n",
    "        point_value, R_value = Dirichlet_bd( maxIter, clamping_ticks, K, Kt, R, Rt, top, bottom, left, right)\n",
    "        \n",
    "    elif boundary == 'Neum':\n",
    "        point_value, R_value = Neumann_bd( maxIter, clamping_ticks, K, Kt, R, Rt, top, bottom, left, right)\n",
    "        \n",
    "    elif boundary == 'Perid':\n",
    "        point_value, R_value = Periodic_bd( maxIter, clamping_ticks, K, Kt, R, Rt)\n",
    "        \n",
    "    return point_value, R_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Diri, Neum, Perid\n",
    "point_value, R_value = compute_fick_law('Diri', clamping_ticks = 200, top = Kr, bottom = Kr, left = Kr, right = Kr)\n",
    "# point_value, R_value = compute_fick_law('Neum', clamping_ticks = 200)\n",
    "# point_value, R_value = compute_fick_law('Perid', clamping_ticks = 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load all timesteps images\n",
    "images = []\n",
    "for image in sorted(glob.glob('./images/*.png'), key=lambda k: int(ntpath.basename(k).replace('.png', ''))):\n",
    "    images.append(imageio.imread(image))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1c30361c18>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD8CAYAAAB+fLH0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztfV2sLclV3rfOvmffM/fHc20PHjseExs0IGEeJmA5SCjICYEYRDI4CsR+wBZYXKPYSiLxgE1QQOTFIRhLKImjsTKykWIbK5HBsiyMYyWCBxw8A44ZYxyPzQDXM77DHXN/j885++5deehap79ee1X/7d5n/5z6pK3qrq6uru7d6+tvraqulhACMjIyMhQ7q25ARkbGeiGTQkZGRgWZFDIyMirIpJCRkVFBJoWMjIwKMilkZGRUsDRSEJHXisgXReRJEXn7so6TkZExLGQZ4xREZATg/wH4AQBXAHwGwBtCCH86+MEyMjIGxbKUwqsBPBlC+EoI4QjAhwA8vKRjZWRkDIgzS6r3pQD+itavAPi7qcIiEnJwIyNjuZgB10II39RUblmkIE5exU8RkcsALmvhvSU1JCMjo8A+8Bdtyi2LFK4AeBmtPwDgaS4QQngEwCMAMBLJL2BkZKwJlqXaPwPgQRF5hYiMAbwewEeXdKyMjIwBsRSlEEK4KyJvA/AJACMAj4YQPr+MY2VkZAyLpXRJdsVIJOSYQkbGcrEPPB5CeFVTuRz0z8jIqCCTQkZGRgWZFDIyMirIpJCRkVFBJoWMjIwKMilkZGRUkEkhIyOjgkwKGRkZFWRSyMjIqCCTQkZGRgXLeksyY8MwOoFjTE/gGBmLI5PCKcJJGH7f42fCWB9k9yEjI6OCrBS2GKtWBl1g25qVw+qQSWGLsEkk0IRMEqtDJoUNxjaRQBMySZwcMilsGE6KCIYINs0GqCMFvg6ZIIZFJoUNwTLIYNlR5rr6hyQMvTaZHIZB7/tCRF4mIv9LRL4gIp8XkX8V839JRL4qIp+Nvx8errmnDyMsTgg7id8qsYw2DXGtMhZTCncB/GwI4Y9E5CKAx0Xkk3Hbu0MIv7p4804nhrixV230fcHt7qsmcvxhMfQmhRDCMwCeicu3ROQLKL4MldETQyiCVR27Dboa5xAEARTnlomhPQZ5oIjIywH8HQD/J2a9TUQ+JyKPisjzhzjGNqOv7O0rvUfO7yTgHbftsRd1M7Jr0R4Lk4KIXADwPwD86xDCTQDvAfCtAB5CoSTeldjvsog8JiKPrX6S+dWgz43a1TD6GmHKgLv82mJRkuiCTA7NWOi7DyKyC+BjAD4RQvg1Z/vLAXwshPCddfWctu8+9FUFQ9e9SuPoIufblu3jYpwmt2Lp330QEQHwXwF8gQlBRF5CxV4H4Im+x8jIyDh5LNL78L0AfgLAn4jIZ2PezwN4g4g8hOIr008BeMtCLdwy9HEXhqhzGcdtQt2T22tP6qndtjdB29xFMeQg5DzyZ+NOCEMbZVN9bY63ym7LNobbZKxN27u6E9tODm3dhzyi8QTQJXi2SF1NxxlKdXSFZ2y2LZ4BNw1lbhrJ2LVLM6uGApkUloguBraIMkhtW1RtdKkLSBteG1ehiSTqCKLNMOcdp04Pech0JoWl4CTIoAsRDKEg2qALcaTapMZY95RPGW5b5ZDJoR6ZFAbGUK6CV89QRJA69jK7KD1jV7CReoogpSJS6mFocjhtxJBJYUAsSghdVEFbIhjiWItimqi3r9F7Ru0RQVNMootLcZqIIZPCABhCHbRVBm3IwJYZwgXpi1niOHVE4bkOdQqiiRw0vy7g2UQOp8mdyKSwIE5KHTQZetP2RdTGIujylLVE0UQQTeTQJeaQVUOJTAoLYBFC6KMMuhBBU1mvzLJiCkMYkkcQKQLpG3PIxFAgk0JPLBLR7/MUH7XclirXdMxFxzA0GYk1uCbDYqOvUxBMAF3dCq8Nbcc2bDMxbOpcHBkZGUtCVgo9cJIqoeuTP1WmThkMEXBsO0Jxkaer7stPfVYFfWINdW1qcie2VS1kUuiIZRFCW3dh5GxvKlMXT2jbU6HbpjVpqh7bnZgyNG+bjSV46EoOQxMDH2sbkEmhJdoE4dr2MLQNGDY99dsSQep4feIMOw2pos3Tua47UI0xZbC8rYtySPVgLNozsU2qIZNCC/QlhKYehi6G6+1XRxRdy6faoXmpcQUKa5xNXYtsREMblB7PI4A2qgGYb89pIoZMCg04SUJIGWjK4OsMv0vZNgphF/UDjnS73d/rOvSUQp0r0gZeHXVEZEmjzWCn00IMmRRqsGxCaBs3SBlxWyLoqixsexQeMTAheC6EpxhSSsFu8/KB7kaXqiMTg49MCgugiRDqeha6yPsuZFBHEE29GinloGVmKAiAsUv5nrFMKbWGaduTIoKuBtZEOnV5tt19iWGTkUkhgSaV0JUQ6tyFthJ/x1nW9SYSSZFCm8CkGsYYafB27+UlJo2Um5FSCsC8caZ6KtqoOz1eXcCzKeaxzd2VmRQcrBMhWDLQ7SNnnxQ57MBvn0cgXtuAeYVgjWLXLNuoPlAlA2/EItdpCaJNFyLv3xSz8JSIPafTSgwLk4KIPAXgForzvxtCeJWIvADAbwJ4OYrJW388hPA3ix7rJDAkITR1NbaJHVh1kNpmy3hl28QUUjGGVJ7C3vxMDGzoXM/UlPUMtKlngPepG9PQhRi4zW2IgctbbCIxLDxxaySFV4UQrlHerwD4egjhnSLydgDPDyH8XKqOdZq4NXXjtxmD0JcQ6kjCbmcD95RCSjGkyqaOa8+Hz6Np4BIwTwhTJ98rM02sc763zdbrlbXHTbURJt9uSxl5ihjWhRRWPXHrwwBeE5ffD+B/A0iSwrqgrT/qle9KCE0xAF23Rr1rtnnlNF99fC+mwNv42G1iDGdRGICm1rXQHomZWbf5R7SPxhRGtE3XbUDQUw99n8h1Qc1U70TXY22aWkg9ALsgAPhdEXlcRC7HvPvjB2j1Q7Qvsjvlz8ZlZKwnhlAK3xtCeFpEXgTgkyLyZ212CiE8AuARoHAfBmjHQqhTCW3HItiydcFDzatzI+zT3y7XuRUc9QeqPQcpFyOlGqam7MxJtRwA7MGX6NY1sP64KodJrGNG5bx4Q0o9pMDnkRqzYOMWqWBi18DjJqmFhUkhhPB0TJ8VkY8AeDWAqyLykhDCM/Ezcs8uepxloishpPZflBC8bbtm312Tz8tjVA08FYPwjlUXT9ijbR4hWGNIuQscB7CBSE/Gs0uh7oTCEpeC25+KdTQRg90/5V542AZiWIgUROQ8gJ0Qwq24/IMAfhnARwG8CcA7Y/rbizZ0FegaXLR59iZuIoSmJ7/Ns8Tg7TOm7TDLlnBSPRGalxoH4AUc7RPcBgQndHxdnsT1I6qPiUe3j2hZiceqibquwiZi4Py6UY4pI2/qqlx3LKoU7gfwkeJbszgD4AMhhN8Rkc8A+LCIvBnAXwL4sQWPszQMHVxsSwipQKFd12VPDbCa8MjCEgUTgyULLzjpnaeud+l5ODLrk3gczVdCsASg+exC6P6gvCNUSaMJ1mjr3AYv0GmJgfdrQtfyq8Cp/ZZkExmkehC89aa+/raEYMvVPf2ZKEB51oVIkYiXphSDtmmG+esCpF0I3WYVAqe6zK7FEe3H60dmH1uWuyi1TtsG2xXqdYPa9nvrnJdab1ILJ00Mq+6S3Gg0EYJX1iMEa1RNhFDnDuzCJwo42+r2tfV7ddlztn46b2vjOrDUH6OU/JrPioAVA69ryvXNnG0eUi6P3afNC1xdYgab6kacSlLo6jKk9m9SCKB8r4wlhCal4Bk/aJmVAasJz3WwLgkwTwi6zgHDOqWQ8vn3UI0XAKULoXlKDJqqChhRauv3nrSeS1M3qGgH8zEGW1eKGFLuU1ssuv+ycCpJoQ5d3AYPNi5g83ibJQQu7xGAZ/RIlK0jEc+t8FwJUHnbfqsQ1PiBwmC4W1K7GNXoVTHYmMLIpDtOCirjkUZXI+PzsP9tXdxEYWME26AWTh0p9FUJ3tOf8+u6GG1eihDq3AFLCLuYVwojyk8Rg63f66Fo60JYo/H8bSaEPZQxgFTvg00tUcAsH6EkGDuysivqiAGodyP6PvXXUS2cOlKoQ9tYQlNPg0VbQkgZuSUCNn6gJIs9s+9Zpy6tX5/qdps3DsJTCArvaZoiBDb2PQAHdCxLCIfmOh3QcTngyG7LUMTA4HNrG1/YdLVwqkihTiXUEUKb/UaokoS33kQI3pPcIwImAT32HubJw5ZtoxTYjWkiBs+FUNhRiU1KYYrC8PWpzwpB15VE+JqnCGIo1A1wsuBr0YUY1k0t1AVtMzIyTiFOjVLoG0uw8HoYrO+taiDV07BjyqZ6E1glqBLYM9tBebaMpx5A20Dp2Zh6IyCtm9QlngCUT/NDlG6D5z7ok1/dBZvq8XU/Dj6ykmgDr0syhbo3JXdMmb5YJ7VwakihDm1cBy9uYPOa+so9dyEVCGT5rzKfjd3GFM7SNo0nWFdC6+cyuq71eF2V3CWpvQdMDNxDoOdqX246QkkEZ1EQhB5HDZ/jCjso4woHtG6vJxMMd332Qcqw9X+texfCy9vU2EImhRp0jSXwujc2oYkQOMjHhOCt75nyY1TJQst4JGJVBZyUlYJwkMF5E2ocGeGeWCzMgHMoDfQAxTob8BjzSmGCkiD0WqmCYBVge3a8UZJdUaccvCHNmxQj6IpTQQp9A4xembZug3Ufdml7SjFYI7eEYJ/+bMB78InAKgfAJ4QRgNEOZWjF9mTVGviRboYwCoptYwA4iKMYZ/Ougx7mAIWSUPWgPQ3WdUi5bn0CY21enLJoevfBIwpPLaSOuy5kcipIoQ/aEkmTKoBJd1AlDe4FSMUUmBDY2PVY5+Py+bjvOZPyPto2JYKxNsLzJRS7tKM2sC6owMGEKBFGAEaHwHgKHM1KY9fqWT3sAtindATgDl1HXeZuSg+2m3SXlu3/lRqjYEc8cr3bqha2nhRSN433dKlTCXabdRu4jE2ZBDyXQYOB7PNr2RQhpJTCeacMKwagsFPZQVUypPwYBb8+mVIK9o0l9SF4sMIBMJ4V6V7cTwmAlQEb/Z14XkoivJwCj7Dkrk6NiTDURbH7e9wHKmvz2qoFJMqlyp40tp4U+iBl7MA8maTcBi/OMDLrTAjsPqid1hHCOVTdAc7TfVgxjAGM+ADnUSWEVHcEUCUEPYE6pcDdDUBpwcZ3kDtF9vkZMJ3MKwQNLu6gVAxajebVgZUIN03XdQzFiJrPLkIKfKrsDnjksWoD74OtJoWhVELbenk7P3AtAbBi4HXbxWgDiNb4z8Xyuk3zWDGo7cs55wDsX9RFHVl762uOTAzc1QAAF1FYpD7SbZRxD6X1o8gbTYC9/aJ6dQ/YZYjFAFT/F48Y1Njt6/j6zsXMrNeNhKyLOaSMftPVwlaTwpDwApJ1wUUub8crqKtgXYcxlfN6Eiwh6I3MBGDLyC4Vtv5Gqi8TqFoUD17g9541vYhq6J9iCTiH6gAF7o80QQXZBUb7wMVJqRRSMQUlRKBKDDo/Aqh6JQJLADx6UjtTOF7QNJZhE7sb2yCTAvwehS6xBAWThHUXuPeBg43sOnDPnxcXsIRwnsqfj+V02z078JkCmI9EsmKAk+pBvAkRVCmMUCUGVgqsGPZRKgWt30QV5Q5wfgLszIpq1eijt4HzSBOD9w4GCx1gfkJYOwuTt5x6Tdtut7EF21OxCe7E1pJCF9ehCd4+nkqAWbeqIRVk9LoYOY9jDEwIWp7JYg/AWHe+iHliAOYDDtaFYN8HKN+D1pcZrFLYQbUrQAlAg42e66CpnvA+St9hF8At4J5J4VbwNQeqJKGnZF0GVgwMJQgbClFD3jFlNfXuJ6sM+qiFdXQhepOCiHw7ik/DKb4FwL8FcAnATwP465j/8yGEj/duYUZGxomiNymEEL4I4CEAEJERgK8C+AiAnwTw7hDCrw7SwiWjretQF1z01IJNbS+fVQ22S5LdCtvN6D34VRSMz6IqJc5RAT2xi/FgXLH2SIAq5llT+HNQHKFTt+Eiyt4GDTQCxWNd9T4PZ1RFABRfIh3HlLsV7sRxFIfFYW6hxNRUy+tTzLsBGlfYNafBKVDthRhRCrPdwnMTUi4EsN5uxFDuw/cD+HII4S/izM5bgzYBRrtN97PkYgOL1o3gV6G5K1K3WWLQEMEcITwvVqBGfxElMWjjmBh2TApU4wlKCGp96jbci+o4ZK3/ENV4gloup2dRWjkTwi1zkW/G8zwsq79lquWm6TpQGh43m1M7doGN1wYaOd9DmwFNm4KhSOH1AD5I628TkTcCeAzAz3pfnI6fmLsMxGGxA6JPPKFubEJTF6RVCTbQyESg9zvHE9ilB+aHKttYIYsAgAhhL248i4IYOK6gxACUqiFFDF6gkS3LDl5SpWBHF3E8wSMEbY8lBE1vxvO4UxIDUA0wMgFY5cDxBQ6HMDGoGmDho/Va1aB1pVAXVKyLEaxbXGGIr06PATwN4JUhhKsicj+Aayi+MfnvALwkhPBTdXUMPcV7G1JIkYB9urMq4O2pqdR53AH35O1SnjfaUH9A+ZC36l/X92gdcBQCp5YYgKo7ocSgAUfrQgBVK7PTKNlhiEDVdQBKTb8f01sorECVwlFcnlB6EwUJcBqDjkeHxSFuxUPto1w/onXtrTiM6/qbmPQI5QtZPELbTiPvTT2vqJsu3q7DWa+bXHYonOQU7z8E4I9CCFcBQFMAEJH3AvjYAMdYCdoOXvJiCOxKsHJghcA9D6nBS2qvuyjsnIcdjDlA4RGCJQZlESWCi1ShRww8qtE+avXRqmMUWCno2AR9dJ9FtcdB3QWOKXgKwUtRDJOeTorDTGP1U5TrKmI4JDI1P6sYZlQOKBWEjS20mZdx012IIUjhDSDXQb8hGVdfB+CJAY6xMBZxHTzlkep25DEISORpnV5sQbcrAbAroa7DGDQOAai6DE0pkO6uVGLQOhXcV8pjFS6iGmzUfTiW4BECdTsen/BNNBOC1j8D7pkCs1mRpcarRMBEoc3mHlQNOmqzOTCo4DybKiwhWBdjCBdiFVj0W5LnAPwAgLdQ9q+IyEMo3IenzLalo8n/r0MdcdinvndMG0sYmzKsEHZonUkAKAmABy3ZtOIOaCzBUwyXYv4lOgAHHrUeVgqWGNgpP4uqhd2L8gUFVQoXUQ4isIOVVCFoHAFx+RKA65QyIVxHaf2KKXDusKoQNNWmsmxXklBe4+uvsQX+FoV6Q6wgPCNn1E3S0tfgV9FbsRAphBD2AbzQ5P3EQi06AbQhjraDnKzr4O1vux09haD5ui8HIDk9jzh0WUkAKNlCjdwjBA1EgPJGKIxaXYmUUjgHv/dB0+eh2vugBKBRP0sIN2J6nS5eihiux3rV6oHjqKLMipGPqVmj1ecfm+ZaxcDdldqrykOetYmeC5FCGyJYVSCxCVs7opHRZRRjm16HNq6D7a60ZSwh7FIZHaNgeyHOaRk7ztk+9ZUYmBBYKVwyZTnGAEpZKZwZA7gL4EKRnomDDs6cAXAbGJ8B7h6V+6lSuIV5haAkoI+T6zGvjhAuorSy8zgeLblzo6oQ2IVQkrWEYBUDnFSNv6sLofv2Gdq8Li7EqSCFOjSpBu4y98p7roPtsbDKYIR5QlDjB6qBR44tjEFvO/IbUZ47oEqBCeESHSBFDNoteYa6H85cQoUQjkch3Ufp7UgQAHA9Ft2fjyFwQFGVglUImlpCULIiX0HOAeP9khDYTWB3LPWhGb3+SiLcDcldlm1dCKCeBFbhDnTF1pDCIrEEoN87EUwI3mAlrwvUBhUtMdg8jS2MgGI+BGYItVvtWuS+S9sdqQQBFC7D2GzbAXBBW3xfTJVFgJIQLqAgghcDuI1jQsClmGrZ68CZCwCuFce6PStSVQW3ALwgFr+BeSJ4Hkoi0FT9ATN5y2gXxavXKLsV+RpbQvC+XakxBfufaZ4SfRcXgrGoCjhJV6OPLWRkZGwxtkYptMVQQUbrg9r6bTyBexA1tROucJckxxN2EEUBj3fmV6F5IJIXX1BX4V468L2UjgDco5rkEorbglXCvSjHDN4FcD8KtfBCFOrgXhQ+gh7gBo7VAvaK9MIB8I2jQh3ciO29Yaq/BOA5VNWBjoHQMRGI286hMm+8WT1WDXq99frqtyetarCxAo0f2Ce0jS/UYVODjaeOFNoiFTvw4gu2K5LHIcBs41GOPPx5j8pYQthDnFPxvCms7oMGCXmEohIExxP04EwIYwDj56G4Fe6jVAsqrPtg3QbrPigxXIsNvgbccxc4ullUewNVDuF4wtfhE4JaOQ+fjoOkZALszeZ7IoBqnu2J0HX7n/HYh1Rcwd4jbcckrDu2nhTajEq0sQAuYwci2f1TXZGpeIKm3ONgFQRQBMyUEEZAVTLwfAhASRZMDLZXgcc1zBHCHkpCuD8Wug8FCQClalBi0ADjJUpvx+1AYd1KCC8FcBVFDCLGF0DEAFQDihOzzoSgvSE8KSxJhNF+qRTOYn626FQ8wcYVgOr/2LZrku8FjyBYOaRIYx16ILaeFJaBum4sr6uSl9lFAJXnLkn7PtJY/Q9lCJ1Cjd+IYmLgcQz3oupKAA4hvBglIbBS8AKN6kK8EIWpqetwAXOBRjyAeWL4WpUYgJIApqh6KjzVm64DpVJQy49MMN4BDmZVFwIoiYBfeGKiPkL6P1Pj14Cj10W5bciBRoOmC1LnVqTK2S5K6zpwyjM0H6c2yKCEMKI8di94MBMTxSWUg5YqhLCHghCYIDRft51FYdy8j6Yvjdvvr6nn/uo+4+eV7efuUW4vv4uhpGjP31wX7wVP7zqzQrNkrnmMOjVZF6faRAPbxDbPoW135CLdltzlmKrXGv/I+dmbK+U6AOQ68OSqrBj2UCUEdS/YaNjQ9CTuGaMaQ3ghqi7EizFv5A/EbQ+gSgg2n0nkxXScvXgcOu494+JnYyHcfh2xpednz98oKM3Wy6bX0l5n7z9LxYZShMBYxJAWuS+Xga0ghaGRUgNdyngxiyZC4AE3x0FK63dYg2CDSX0iSg3twk4ci2BjCLruGfQZWn4A1ae/l88K4WyifiaKvaJNKUJgpXAeVdXADLBbvV783ol3nYF5lVenFIa4JxbFSZFHJoWI1B/aJciosDEGzuPRjTCplld5eyyBrQGwlNC8OmLQ7cfQbkc2UFUMHFOwRg7ah1NgnkSQqJf3VX+G2ldHCHw97HWg66N8qe5Y3fXm900A/z+z6HNvrJsaqEMmhRp0Ga9QF2S0N0pKMbCLcuxOmCdhxWnWg9o5Fe1EKcfzM96H+W5HTtlw2cj1JQUlgBealPNfaPZvOl7cpi4Et9tODXcO6enozfVhdywV5LX/Q0rdtSEKxqYb1aa3vxZduiOHPA7fQPbGS92IfBPPxRn0cWY/1GIH+APzBqUpv8tw3KvABjrCvOGqkYNSa+wvNtt5ual+UgraPttuS3jA/DnrdRlXN+8i7TKk/g/vP4Mpuwjq6h/yOItg1cffOHhKgOE9dZoUAq/bmMLcXe2F122eGok1MDZC3Iuia/ESbdMxCEDVyL2nfCq1++t4B3scHhhlXAg78xOfnz1XqxRi6sUUUtfdM9K2LmJdAHpTsW3nk5GRsSAyKSyANhHqEapDYvnpxOPt7VNqiji0GajKB54oQJ+cOiOS5mnK0xHNvcvAebrOw5v16W5dDVYUd53tXMar1x6X8uz0SXw+Os38jPLs9YjXSXbKa972uvP/1KfnYZuQSYHQ9WKwsdeNbkv5sbtmnW/I43iCZpgb//gtHp1rzH78ACgJozJByoVYSIck60jFS5ROad1ua0pT9fCcDHpcPQGduGVcNXo+HyU+njCBrwe/aIJqjKbNdefUwpJI11GMm2ZkrdorIo+KyLMi8gTlvUBEPikiX4rp82O+iMivi8iTIvI5EfmuZTX+JOE9dfjpMnO2M5g0vK8P2eDnDJi/e/muVmNRg7EGtEv5lQlSeKgyEwK/y2CVQhdi4OV7Mf+OBB9XBx5Qntd+O/Gi5tVcH/2KtGf49t0F+4YkavbxVJ9d3nS0JbH3AXityXs7gE+FEB4E8Km4DhRTvj8Yf5cBvGfxZq4XUi+7wOQr1MjrJKyVrrtAeWeDUpbMaiRKEvyGT2VKY37DkZfVMC0hsAHzm5EeidjttowlBktIvHxhfpZVO4eakoV1Hfg6zeaVgqfq+PrzpQbmjVz34fxtfO8BaEkKIYTfQ/FCK+NhAO+Py+8H8KOU/xuhwKcBXBKRlwzR2HVBGx+TZwzmG65ObUzNcoVN+FFljUSX+SnK7whXplBT49VXhqxRH2LeyA86pDbP1sdt4fpjW/idZqt89ALtmtQ+/k1MwXvd2XPpWF2k5l6sW98WLOLu3K/fd4jpi2L+SwH8FZW7EvO2BnVPCPv00Rut7itBfINWbt7Ud874fd/UJAH8pHWNkI2R0/O0fhvlK9G6z22zj91uy5xPHMcjkb2qMvDOy+Z5RDmbJ1lQUXv99VLr/2RVXQqnWil0hPdpyLlv04nIZRF5TEQeW+zDdSePuieE9VNTSoHrmtI+c0qBX/RnA1DD8CYJ4A+07ALzRr6HeeO2H330jJwJgMvchk8itj5LJnuYawsbu3deNo+JUlNHKQBVguDrr5faUwp1cxtkpTCPq+oWxPTZmH8FwMuo3AMovjVZQQjhkRDCq0IIr9q071SnnjaMJqVgb1RP7h7vzC/vWwOws47qV0zYeHAGpTFaQjiDqlGfxbyR73VIbZ6tj9vC9ce2WGI7QJUc7flzUIdmMvHcMUsQnKaUAqPN/74NWIQUPgrgTXH5TQB+m/LfGHshvgfADfqM3FbABps4ks35CrZtXbfRb76JgWjLzCaaqrHozCBWTfCT9lAbq8Z42yyfQfkk1/kUb9D6GcyTxhlU5l6c227L3DD13zHH5+Xb1XZ70yIp+fF14OsTL7bOx2Kvb+r686UG5g2eQzqKU60UROSDAP4AwLeLyBUReTOAdwL4ARH5EopPx70zFv84gK8AeBLAewH8i8FbvQKwwXpP+R1nOyPV++Dtc+xuTCiDUzUSdclZIfCTVfPnjNAaJhsk3tSmAAAfOUlEQVS5GrLmX6flNikvW0K4jnlCMm3x2s/nxXk114ef+nx9ragAfFJP7eOpDbu86Wg1HVsI4Q2JTd/vlA0A3rpIo1YF2y3VBCYE7SXz4ElVoJxQ1AtCTgDcw5lHtEFTZZlDlF991vQ8rd89Kr/kBKDoJlTf/3YsfB3Fy03XUby3cB3lSET+QssFk7Kx69ddzqCYhm2EYnpmXT8T1xGPO6Ll+IWpYyLDPCHoN92YEDTV6zOtZvMl5OvO63Vkbg9niaENVj3nYlds2mCrjIyMJSOTwgJoCjzN9Sag+tSyfi3vNwIQ7KPN+tHsPhxSnqYjFB93HQHlUxwop1G+btav0Tq7D7xNVcE1VNWBbucyXr32uJTH7eXzAMovW7OSsNcjXqdgAo1A83X34jowZVLr24ZMCh1h3VeLOn/T3ogTZ92q4TmZbI3Fy1OCuGPSipGz4es2NXaglPmgvDap3d+ShyUMlNtsez2is+eqqblOx9cP/nXm9RQx2zwG3web5h40YatJoclw68oschy+2ewNl/Jj2Qe2PHB8h1sjYRbRvP2YWgO7qxuA+af7tdgSa+TPoTRuTb+GYlTitZh+zWzn5ab6iaS0fbbd2mwmQXvOel2OqpsnmL+WqfhO3X8GU3YR1NU/5HEWwVaTwqJo8+fwHzyFrySabkQOfFmPYe5JyE9IPag1HGtQt7RifZKrUcOk1tj13QQ18qsxtUTB+c+Z/ZuOF7fNYju53XoeTHj78ZxTSuGgugrMX9MmguZtM1Rdii73xKYik0JEym/02LyttPRupqOYZ59gEyp/FMsfPxyVAOyTkeX0nbjTHZREwYRx7EIAxdP5LgpjZsN9DlUD/hpK41ZFcNXse5XK3jVlvXp5X1IKTGTcfj4vvh72OtD1OYi76ADIuutthzm0IYA+98YmxSEyKTho84d2CUZZ9WAVgid3j3vY7GOODUClRR0hTFE8gacoPgd/W+mGDVTXD+L6IeaN/C6Kwaq6fpDIZxI5TNR/l9YPijZxO+uIQZctUU6q18uLKaRciTp113a9bZlFcFLEshWk0PZiLXJR6wJKfHPZiLf9NREDMP8AxCGqRMAGoXlqLClCUBdihuLrz3MKgYmCDVqNWg3/CkpC8PI9YlECMIrhG0fFT12HFDHw+dnz5/XDeUGl17KOELz/CUj/Zx4WcRnWTUXkb0kaNA1gmprtuu4FG/WdfksW+sIff+NQU2D+9aHx1GScRWEsOkMRv3uk28YoDGw3piOU8cUXoPj68xgoDPfFKAz5/rh+X0wv0QEuAPgqioFOOqxZexKeQ3VgElANZrLLEMnk6GZ5Qa7Hi8LEMMG8UlBrZ9I4qKZe54y6EJzOUCUINn7NY6SIocmgNzG+kEmhB7iv2051sINynN7M7ANUX+4DyptxQmXU7tXWj2bAeEYZhzHVeMEYJSHcovVdFD1/I5TEANBn4JkY7kNh9PejeOrfF1Og+lVpJYTnMP9SlcJ2bxqX4ehmYZnaK6mEMIl5k9h+JgpdB0pCMArqaFZeP075+lqfPzVSGpgnc1tmEw2+DbaeFFJPfv2jedadKdIz69iXX7Ss3Sfln/KUAHaKAB3urDeo2v0YxY1deft5P2awUrCEwAqBlYJG93fQQAya6qSqX0N1ePM1SvlFKSYGHrhELgMTgkojJQSrFNTlYWLQC7of91f34qDqUajbAFRdByUIL/hYFweCKZOKKcycbV7gsilIvUpsPSn0RR1BAPMvNelgO1YMR6asnfqAXYhj46flygwEM0DumA2qFNTwLSHsoPoqwwtieSUETac340dn7VRr/K7EDZTvSPB7Dzo46Q6qUgQoezmuF3V+46i4EEoIWuwIpQtxPa7bWIiuA/Ouw51iFKMXcoEpal0HfmVi6qzzsiWEunjDusUJumArAo0ZGRnD4dQpBasAPDQFG7UMxxO4fp0ukcvq08lOJLRr8oGqctB0H8D5mclQacEqYQdVd0EDjJdQPpn1wX8xNvgSgOkRcOEuCnVwX0w10KjxBFUJz8X0KtLugxk1eXtWXAQOKupjWeMIvI3LsFoASrdhgmN5QF5ERRkAVXdBxy9oPijfiwGlAo5DDWJaR0WxNaTAMYI+aEME3jE5rjBz1rmsnS/EEgLfxJp3hJIHphNgNEGVGACfEG6isGlLDNoY/SzDRRQ2fhHAZBa/8vxs/LajGrl9bVrdBnYflBiAytBlHj/Bxj4pi7mEcDOepyUGoIwnHBRlp5OSCI4wHzOYOnleGcD3+5Uw6uIJTVg0VnCS5LE1pNAXTcrB+zO8LslUXAEoJ0maUFm9IZkg1AfmSYZswHFvHxAlBJ3AoUkhcApaVqVwMR6M52TAfvlNx+l+/KCMBhh1RmZDCHePyhPQoOAtSj1iAKrxhIOYHsJXDEDZXTkBwn6FH3BI186O95qYdGrWtcsSqAYNvXgCzHoqyGixjsrA4lSQQhcVwCSRIgxmfds1qTcGd03yPilXQgOPwHzwkYXBGMBIXy9Ww/MIQZUCE4J+6v25mDdB4UooIUxQfLZNJ2lRJXIOwOwoTh9/MzZ4vyp9Do7Ki2yHVtuBVDdiykrhFkpCUKXguRJabyQGHa+lAkoVw4SK2x4HdiVS3Yw8dsErY8c1MKaJ5SasQ88DcEpIwWKouALXleqaBEq7UcXAyoGVgt7ETBaWGO4AuDihngignIOACeF5mFcKeodeRPEVD6sUlBCYGIByGsdbKLo+LaPpXAcsdVLvMnhGzu6EpoeYJwYlmziOIUxKz8TGEnjcB7//4CkD+34EMG+gXV2HTY0nAC1IQUQeBfAjAJ4NIXxnzPsPAP4xiuv5ZQA/GUK4LiIvB/AFAF+Mu386hPAzS2h3Em0MPoU6IrB/chsXgkcoaruUCKxiOJ5jFcWNzfXsO+l5NRZuEBMCEwOnejJKCOcxTwh3UFUKPD7iJsrvSuhJ8nxlegKATwx2YBJQKoPDRKr7aL2RJLQ6m2psgV+L0HFOrBC8QU38zkQb18GLRVhs2hDoNkrhfQD+I4DfoLxPAnhHCOGuiPx7AO8A8HNx25dDCA8N2soBUGfw3pgEdiGAeaJp60LA5PGHjjzFAFRJgpXCDgrb2AGwMwPu0R10iIAlBCaGi9RojiGkFMIhyi9Zq1KYoTomggmBlQJQneiFlYJaLo87UNeBg4weQSjJTIBvzIrD30FJAOw6HFL1KYVg4wE2yNjVdeA6AZ9EUlgX1wFoQQohhN+LCoDzfpdWPw3gnw3brPVA3WhIj0Q45ZgCKwgmAE5HqCoFditUIShR6PooksJ4B82EoG4CUBLCEcrgoo0lqGIAStdBh13qhM/Wv+ETAKqxBetC2C5GG4BMKQYAR5OqMtDgolavxs9KgX8zk7IisDEFOCkS64p1MvA+GCKm8FMAfpPWXyEif4zib/yFEMLvD3CMwVGnAOqGRQPz3xr0pmpnSaoEoXlebMHGFNTGNFWFoG3meZRwGIOUasQpQmCloIQwQTnzM3/hTdUBUBKCHT6tCuE42Eht4hcQtFvSdk9qe7h3wSOG6G4cHZZVsstwhGpsQV2IVEzBxhJYyfH/x6og1cuQGs8As91iXeMJwIKkICL/BkX/1H+LWc8A+OYQwnMi8t0AfktEXhlCuOnsexnFV6nd78wtgrpeg74uRF2cQv94nqBVlYHm87dbPKVgH7Q7qLoM7OJX5kvR/ZQYpiiVwnlUCUGVQh0h6IF4GLW+V6HBEdttwtOwx7Ycn4g2OEUMegEtMdxBRSEcHZbu1X6sWl0HG6qww5x5UJNVCuzKWeP3XIw6DO06rIo4epOCiLwJRQDy++O3HhBCOH6FPYTwuIh8GcC3AXjM7h9CeATAIwAwElnbz0laIvH+KC/+wP6oug+2FyL1oOXJmXdQJQ1POdxCtHclBiWEGWJ3Ikp1AFQJQd0Bz3VgpaDEwB9/tTEFbrTnQqh1qgVzTIGJga2eCIG9DSYAVgYaZOSYgm7zYguaD1QN38YQYNatYoDZb5PR690HEXktisDiPwkh7FP+N4nIKC5/C4AHUXwtKiMjY0PQpkvygwBeA+A+EbkC4BdR9DacBfBJEQHKrsfvA/DLInIXBYn+TAjh60tq+yDw3ARVB216IXibbueAI1B1IyYou/rZfTiicqwYVBWouvfcCaB4ip4DCrVgZYq6EedinioI7nmwM7to7AAx7yaqPRAaaOSeBxto5Igfa3g7f6TnA0Q34mhSxhBU6HjuAq/rwCV+S9IbqwCTpnoPUqoATlle9uIN6xxLUEhU/ivFSCTsNRfrXm8i38qjkbO8U7POebq+g8JGdH2Xyu3S+i79xijsjPOA8tvNuu18TM/Ffc5RPmh5j7bds5PYGXF55KT6J9hUD6JuA8cVmsYqcGoJQS2E50Yg6/7GrAwmqjeBapHUrjgwh+beB35HQn9TWgaqBMKxCK+XAom8NqRwkvGEfeDxEMKrmsqdyhGNFm3Vgm4HfMLhP3inZtn2RCi0zgmq9Vv3HKj2GJ6nZRUEsxlw7hAQVQVTlEZ+hLJ7kad42o/rGk/QSpVMbM+DfcNLwV+wAqozJE3juh1dpJa7X4xUTBk7kCaEA5Nq+Rmtp3ocrFLwuiktrLIAuhHCuiKTQkvon6wGXvcH82Am0LISDLsKdvizHoNHDDNBWGIAqmOL1P7VYzg/AXZuAKLuAlBO7eS5DloZz+ykRABURzPyCYHq56ip7Zq0ikG3TYqXm7RzwlMDzCEpQmA1oNV7czN6PQ5s/N5gpW0PMCq2mhS6dE16ZdsOXvK2e9CxCDoFG5x1hT7IvePvO2X5pjyL+bDBeB8Y8dRO51EqBlYKVjEApVJgJrOBE25A6otW3Feo1h/3nU5KI1fO0F4ENn6tJqUQ+NVpO3hphmpsQddtzwTgqwCL1PbUaMc25erqPSlsNSn0hecy8FOfDbXuCWGDlDpOQbsotcyRs98B6qGTimj9OoRAFcLM5O+hHAF5DoBwVyQHF/nNKyURoEoMGmyc0DrMidgInhN01A/o8uxqnhvAcQIkynBX5BEdRsur62BHMYKWrdRPqQSGHRqNmvVNwdaTQp1aAObHIKTUgo0leLGFKapGb28K63oDpWtuex/2nLIWWv+xq4CyF2JCKY9ORjzW3gwY7cch0vyegxZSEgDKsdZtlQLLID0h4Di2cDSrPsl5+EIqHqnDHID5IQ+WECa0Xa+TpwhsANG6BVOTcn5TLMHbxlhXlQCcAlLoizoXgdUC/4kcoLR/OpOENwKSp2MD6glBb0qezJnDAVYhcDcoUCWJg1kkiAkw2kHVhWClwF0uKaXAF4O1PoCpIQJ+ktsBR5YgeB2oKoQJ5gmBexS0PBMCk4A3ihHUNj41byBTCutg3H1xKkihrYGnyloS8IY/dwk0tQlWpuAFufSG5WEETAhKEGyn1nPYi5XvAcB+lQ9EG6yTT+44KbGYugWGF9yYo7bfBghT3YgweUwME2edy7NCYJfA5sMpk3ITmkiiS0ByXYjkVJBCX7Qhk5RbYaEqYpG22PrtU80qAyYIHjfh9UYewHxnAqUHMZ1FPjgs62MPQgUDdz6wAc+Q9CSO85gY+Ok/pe38xXmPRGaYJwhtg3UZPCKwPRFenKHOTeBym4xMCqhXC7rMMQjOs0/9uheu+kKNnQ3Yu2Gt66DGOEZ1cJSNMWpIQYkBqM4L6w3gqgspsGwH6jsjtLwNBlpj5y5GVgheGSYXrj9FCF6MgJFyG5ryuqiEdUL+7kNGRkYFp0Yp1LkCXeCNdOwTW+jqTlhZb59Sqg7GqLoN3mfp7Lcn2G3Q6Ra1HLdXU+s+eK6NHkufsqwaOOgHzD/xvfiAjUF4Za1KsPWn3AUbfITJB+XZ5aGe/uvkcpwaUgD6Bxy7xhaGlolq3HxMO4RA3QY2Qm+OVW0jD004xPyn6+yU85raa2GJAagao01twBOoGqymh2adeytSRFA3EMm6CzxWweuS9AgCqP9v+7oO60QIwCkjhSa06YngcqnYgt1n0T+dZ2WyBKHt2UVVIXDKxADUKwVO+TxS8YS2cYU2SsFLrfEDPhGkxiFw/W0JgVEXN2g71mBTYgmKU0cKfd0IVgKATwy63d4UXl5XsAEzQeg2JgK98S0hqArgekaYVwp1hMDnk3IdrPFYpcBBUC3f1EPgjTtoIgQmkS4KAU5eEyH0/X/XTSUAp5AUmtCkFppIxZPToDxvpGMT9CZmNcAuhI0hcMqEoUaocQaOUXiEoNv52F5k2g5V4JGblhhsV6XNm6Le6AGfLFJqAE5eigg842/6r+z2TVcJwCklhb5qwe6fciNs/awibAwiVX9TW9nPZcNn0uBJYZkEOI8neuGYAr9/wXOotIkppJQCGyM/ybVMnXG3cQe8fFu/JQQ7WInhlfHW+2IdVQJwSkmhCW1jC1zWIwaFF3ysI4ZUP7HezGOnfVYh6MQvvO7FFIAqEehAJF3nNjUpBS+mYMcB2BiDLlu3Yuas2/etbNkmpeARAudzu5sIwcM2qATgFJOCjRFYeMTA5ZkolkEMXjvU6DUWoPXtmjw2et1fZ30Cqq9r8ysMwPxb0nzOvJ5yk4B5Azoy69aV0DzvKW8Jw8uzLoln6G0JoU7xpLonvfUmQlhXlQD0/2zcLwH4aQB/HYv9fAjh43HbOwC8GcV5/8sQwieW0O7B0NWV4PJtiEHLAfPE0OXYKvEtIWg9SgB6Y6vB8gjFVO+DTuTCM0sD88RQ58o0uRApxcDLtofCe/rzPnZ7ShF4ZTgfqLYt1X6YPLvcButMBoq+n40DgHeHEH6VM0TkOwC8HsArAfwtAP9TRL4thLAJ12IOKbmcIga7XW+k1DiGpl4JvnlZBdgbloOMKuM91WLbrAphhHki4HPibTq4aWbSQypfJ8E9qe4ZszXk1FQNmnqEY+tvQwh1asDmpf67TXUbFL0+G1eDhwF8KH7/4c9F5EkArwbwB71beAJoEzNos79VB56KANoRA9dl27ZjUi3DKoKDjEB1diePCFQpaF2sGA5MWUsI1si9rko2TM1PuRkTs0+bJ7+3DKc80J0QugQWu8Qc1hWLxBTeJiJvRPGhl58NIfwNgJei+Lak4krMW3t0JYZU+Tq3oa9i4Lq1Xj6OVQZa1xjzT3bQNi+omHIdWLp7hMDgyVO0vZxaEuBz8OIQnsLg82RlYLcB3QihbWwgReQpbAohAP1fiHoPgG8F8BCKT8W9K+Z7X4Bz55AXkcsi8piIPLb6SeYzMjIUvZRCCOGqLovIewF8LK5eAfAyKvoAgKcTdazdZ+OGDjracl6vRJ9jWSbnp7aW03EGPKU8aH/rnvCYBaseWCXUDXO27dW28bp1J2wZ70luVYS3LeWmILG9ro0w+XbbNgYXGb1IQUReEkJ4Jq6+DsATcfmjAD4gIr+GItD4IIA/XLiVawAOGjIWIYa2boMlKyuRR5RnyWFq1nl/PZ8UITAZ2DybdiEGzwjbxA3alEVNeZh9Um3hfK+dHjY9uMjo+9m414jIQyhcg6cAvAUAQgifF5EPA/hTFF+jfuum9Tw0qYWm+IIlBqAaY0gRgzeOQWHnd0wpBbvONzwbtO154Lbq8TjQqPDiCinUBe+8J7nNTykGJPIsIdQFCusIwZ5TEyE0kcFG3fwRW/3ZuEXQ5EZ4wRivG88rn+r+syMI7bYdZ1nXbU+CPY63j9cWr+11YxQUlhjhLKeMU7elnuhtXZCU0af2TbXPa/s2EEL+bNyCGFIxaHnAf7OyCZ6b4bkUqYFU7FrY9ngxBNsuLcNlGdwFatHGUDW/znBTZduQSNPxvf9h2wihCzIpLIA2xKBIGazN07JdicMeP+WqcNsBn6RSxs9jIezx7Adt+BhclvPrjLaJTLq4BG0JoS4uwtim+IGHTAo1aFILQPsxDE0k4I1vaAtr+E0xi6Z4htdDwj0UnpKoG4TFSLkO3rZUgLJN2VTgMOUuDEkIm6wSgEwKjVg2MXC+Z5iLqAY+jufaeF2knnHXEYKCA5Oeq2PbZPPbPM3blk+5BZkQ2iGTQgssQgyAL7d1HyD9RmXfMQ11beTjWQLwiAmYN7YU2gbt6sp2CUzavC51tmkX47QQApBJoTU8A7fwiMHb167XjWkAusUbrLrwYJWCzbftVqQGLjWNt2ijFFLb2hh0WwVQt62u/aeJEIBMCp3RpBpSxODt66kGIP2ORBvl4JVpQxAKTylw2bqUkTKktk/tJlVQV1cdGTTVY7HNvQwp9H33ISMjY0uRlUIPLEstePu3UQxth0un0OT28HFSaHP8IeIMXco1jT/w1pv2b7vvJiOTQk+0IQZFU8+EF69gArD7eNt0f49EPLQhJ4tF+ufbGGOXEYVtexbqjte3hyG177Ygk8ICaGNIQLeeCZj8ph4KJLb3Qdvz6VOvRZ+neJPSaNrelJ8JoUAmhQWxCDHo/nDqSAUAvQFOTQTRxb2oc336oE3AsS6/DXksogzq6vCw7YQAZFIYBG26K4F5w7Z1eMRg6/Xq8Lo4QeW8GzlFFH0HSdVhiCd2X7XR51hd9t9GZFIYEEOpBqB9zAFIqwdbjst3IYoh0OcpvShpNJ1LJgQfmRQGRhdiAOp7KYDmmAPXxfWlytryqX2WhT4R/UVdkK5t6FLXNiKTwhLQ1p0A+pED59ttTQTBqFMTFm3iDEMaW1PPydB1dq1rm5FJYYnoSg51hldXV2pbShHY/Riptg7ZHdmn/ro6hqi/S13bjkwKJ4CuLgXQXTnwttT2JpKwdSwTQ7xPMKQyaFPfaUEmhRNCF9UAtHcrUnW2GY1YZzRDdEsuyygzGSwXjf+9iDwqIs+KyBOU95si8tn4e0pEPhvzXy4i36Bt/2WZjc/IyBgevb4lGUL457osIu8CcIPKfzmE8NBQDdw2dB012KQYtE5F3duQFm2HaS8DXZ7OXQZeLbMdpwULfUtSRATAjwP4B8M2a7vR1ZUA2pED163o+xLTEMOd+xpc1/0yGQyLRWMKfw/A1RDClyjvFSLyxwBuAviFEMLvezuKyGUAlwH/W3OnAYuQA9DO7++qEOr2Wxb6HKuvkslk0IxFSeENAD5I688A+OYQwnMi8t0AfktEXhlCuGl3XMfPxq0KfcgBaNebUHc8xjJehGpz3LZY1J3JZNAevUlBRM4A+KcAvlvz4ifoD+Py4yLyZQDfhuLL1BkN6EsOiq4qwjv2OmGouMY6nts6Y5Gep38I4M9CCFc0Q0S+SURGcflbUHxL8iuLNfH0YUq/vpjRb5MwRLunGOYanla06ZL8IIA/APDtInJFRN4cN70eVdcBAL4PwOdE5P8C+O8AfiaE8PUhG3zaMMSNPUv8VolltCmTwDDI35LcMJyE7w+sZvBSX2QiaIf8LcktRZsxCUNg1UqiCZkIlodMChuMLmMSNh2ZBE4OmRS2CNtEEpkEVof83YeMjIwKslLYYmyScsjKYH2QSeEUocnwlkka2eg3B5kUMo6RDTcDyDGFjIwMg0wKGRkZFWRSyMjIqCCTQkZGRgWZFDIyMirIpJCRkVFBJoWMjIwKMilkZGRUkEkhIyOjgkwKGRkZFWRSyMjIqCCTQkZGRgWZFDIyMirIpJCRkVHBWrw6PQOu7QN3AFxbdVuWgPuwnecFbO+5bet5/e02hdZiincAEJHH2kw/vWnY1vMCtvfctvW82iK7DxkZGRVkUsjIyKhgnUjhkVU3YEnY1vMCtvfctvW8WmFtYgoZGRnrgXVSChkZGWuAlZOCiLxWRL4oIk+KyNtX3Z5FISJPicifiMhnReSxmPcCEfmkiHwpps9fdTubICKPisizIvIE5bnnIQV+Pf6HnxOR71pdy5uROLdfEpGvxv/tsyLyw7TtHfHcvigi/2g1rT45rJQURGQE4D8B+CEA3wHgDSLyHats00D4+yGEh6hb6+0APhVCeBDAp+L6uuN9AF5r8lLn8UMAHoy/ywDec0Jt7Iv3Yf7cAODd8X97KITwcQCI9+PrAbwy7vOf4327tVi1Ung1gCdDCF8JIRwB+BCAh1fcpmXgYQDvj8vvB/CjK2xLK4QQfg/A10126jweBvAbocCnAVwSkZecTEu7I3FuKTwM4EMhhMMQwp8DeBLFfbu1WDUpvBTAX9H6lZi3yQgAfldEHheRyzHv/hDCMwAQ0xetrHWLIXUe2/I/vi26P4+Si7ct59YaqyYFcfI2vTvke0MI34VCUr9VRL5v1Q06AWzD//geAN8K4CEAzwB4V8zfhnPrhFWTwhUAL6P1BwA8vaK2DIIQwtMxfRbAR1BIzasqp2P67OpauBBS57Hx/2MI4WoIYRpCmAF4L0oXYePPrStWTQqfAfCgiLxCRMYoAjofXXGbekNEzovIRV0G8IMAnkBxTm+Kxd4E4LdX08KFkTqPjwJ4Y+yF+B4AN9TN2BSYGMjrUPxvQHFurxeRsyLyChTB1D886fadJFb6lmQI4a6IvA3AJ1B89PjREMLnV9mmBXE/gI+ICFBc2w+EEH5HRD4D4MMi8mYAfwngx1bYxlYQkQ8CeA2A+0TkCoBfBPBO+OfxcQA/jCIItw/gJ0+8wR2QOLfXiMhDKFyDpwC8BQBCCJ8XkQ8D+FMAdwG8NYSw1d/izSMaMzIyKli1+5CRkbFmyKSQkZFRQSaFjIyMCjIpZGRkVJBJISMjo4JMChkZGRVkUsjIyKggk0JGRkYF/x9m6I6aUuofvwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x115f5b780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(images[1500])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_video(images, image_folder, video_name, outimg=None, fps=5, size=None,\n",
    "               is_color=True, format=\"XVID\"):\n",
    "    \n",
    "    from cv2 import VideoWriter, VideoWriter_fourcc, imread, resize\n",
    "    fourcc = VideoWriter_fourcc(*format)\n",
    "    vid = None\n",
    "    for image in images:\n",
    "        if not os.path.join(image_folder, image):\n",
    "            raise FileNotFoundError(image)\n",
    "        img = imread(os.path.join(image_folder, image))\n",
    "        if vid is None:\n",
    "            if size is None:\n",
    "                size = img.shape[1], img.shape[0]\n",
    "            vid = VideoWriter(video_name, fourcc, float(fps), size, is_color)\n",
    "        if size[0] != img.shape[1] and size[1] != img.shape[0]:\n",
    "            img = resize(img, size)\n",
    "        vid.write(img)\n",
    "    vid.release()\n",
    "    return vid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_folder = 'images'\n",
    "video_name = 'Diffusion.mp4'\n",
    "\n",
    "images = []\n",
    "for i in range(maxIter):\n",
    "    images.append( str(int(i)) + \".png\")\n",
    "\n",
    "make_video(images, image_folder, video_name, outimg=None, fps=50*int(1/deltaT), size=None, is_color=True, format=\"XVID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "K_R_value = np.zeros((N,2))\n",
    "K_R_value[:,0] = point_value[:N,1]\n",
    "K_R_value[:,1] = R_value[:N,1]\n",
    "plt.plot( K_R_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.max(R_value[:,1]), np.min(R_value[:,1]))\n",
    "np.argmax(R_value[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(R_value[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1\n",
      "50 2\n",
      "150 3\n",
      "950 4\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAERCAYAAACU1LsdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xu0JXV55vHvQwN2FNDB1jjSiBcaR9QQFdEVQ0ImRhrGEY2JAW+YMTrRIePETBxjIoloVpZOjDETNOKK1zHeGI2tomQcNEiMplsTxEbBFlRaFGzQDs2tofudP6pq8ztFXffZdc7ZtZ/PWmfVPlW169Tp3qee/f7eqtqKCMzMzAAOWO0dMDOztcOhYGZmEw4FMzObcCiYmdmEQ8HMzCYcCmZmNuFQMDOziYUOBUmHS/qopJslfUfSsxvWlaTXS7oh/3qDJK3k/tp4STpL0jZJt0t6V8Xye0p6i6RdknZLuri0/LGSLpa0R9J1kl5W83Oek69TfN0iKSQ9bqBfzebMgau9A6vsXGAv8JPATwOflHRpRGyvWPfFwNOB44AA/i9wFfBXK7SvNm7XAq8DTgZ+omL5eWR/r48AbiR7vQIgaQPwaeC3gfOBg4GNVT8kIt4HvC957guAVwNfmcHvYCOwsJWCpHsBzwReHRF7IuISYAvwvJqnnAm8MSJ2RsT3gDcCL6jZ9kmSdkp6haTrJX1f0tMlnSrpSkk3SnpVsv4J+bvEf83f5f3ZTH9ZW/Mi4iMR8bfADeVlkh4OPA14cUT8MCL2RcSXk1VeDlwYEe+LiNsj4qaI+HrHH30m8J6oubWBpM9Jep2kL+SVxccl3VfS+/LX61ZJD87XlaQ35a/53ZK+KulRPf4ZbA1Y2FAAjgH2RcSVybxLgUfWrP/IfHmXdQEeAKwHjgDOBt4OPBd4HHAicLakh+brvhl4c0QcBjwM+FC/X8VG7gnAd4DX5MNHl0l6ZrL8icCN+YH7+vzA/aC2jUo6Cvg54D0tq55O9mbpCLLX5z8C7wQOB74O/GG+3lPy7R0D3Af4NSpCzta2RQ6FQ4DdpXm7gUM7rr8bOKShr3AH8McRcQfwAWAD2YH/pnx4ajvwU8m6R0vakFctX+z/69iIbQQeRfaaeyBwFvBuSY9Ilp8JvAx4EHA18P4O230+8PmIuLplvXdGxLciYjfwKeBbEfGZiLgT+DDwmHy9O8j+fv4doIj4ekR8v+svaWvDIofCHuCw0rzDgJs6rn8YsKeu7AZuiIh9+eNb8+l1yfJbyYIG4IVk766+kZfjT+2w/7Y4biU74L4uIvZGxN8DnyV7Z14s/2hEbI2I24DXAD8j6d4t230+8O4OP7/8uq18HUfERcBfkvXqrpN0nqTy35itcYscClcCB0ralMw7juwdfJXt+fIu6/YSEd+MiDOA+wOvB87Pex5mAF/tsDx9c1I8rj07TtKTyKqO85e3a0tFxF9ExOPIhlaPAX53ltu34S1sKETEzcBHgHMk3Sv/IzkNeG/NU94DvFzSEZIeCPwO8K5Z7Iuk50q6X0TsB36cz97X9BwbF0kHSloPrAPWSVovqTg78GLgu8Dv5es9CTgJuDBf/k7gGZJ+WtJBZGcTXRIRP6bemcD/iYi6ynia3+Hxkp6Q78PNwG34dTx3FjYUci8lO/3verIx2JcUp6NKOlHSnmTdtwEfBy4DvgZ8Mp83C5uB7fnPezNwej4MYIvjD8iGYl5JdkLCrfk88r7UacCpZH2FtwPPj4hv5MsvAl5F9pq8HjgamFxzI2m7pOck368HnkW3oaM+Dsv37UdkjfEbgD+d8c+wgckfsmNmZoVFrxTMzCzhUDAzswmHgpmZTTgUzMxswqFgZmYTa+IuqUdKPgVqjbgmwrcDb+DX6trh1+owXCmYmdmEQ8HMzCYcCmZmNuFQMDOzCYeCmdkckvSO/EOVvlazXJL+QtKO/FPwHttluw4FM7P59C6ym2nWOQXYlH+9GHhrl406FMzM5lBEXAzc2LDKaeSfv51/muN9JP3btu2uiesUzMzGbPPmzbFr165ez/nyl7+8newzKQrnRcR5PTZxBHBN8v3OfF7jR6Q6FMzMBrZr1w/Ztu1LvZ4jHXRbRBy/jB9bdXFf68WXDgUzsxVx50r/wJ3Akcn3G4Fr257knoKZ2eCCLBT6fC3bFuD5+VlITwR2R0Tj0BG4UjAzWwFFKMyOpPeTfVb3Bkk7gT8EDgKIiL8CLiD7CNcdwC3Ar3fZrkPBzGxw+1naM16+iDijZXkA/6Xvdh0KZmaDm32lMBSHgpnZinAomJkZkFUK+1Z7JzpxKJiZDW72PYWhOBTMzFaEh4/MzAxwo9nMzBIOBTMzm3BPwczMlnClYGZmgIePzMws4eEjMzObcKVgZmZLOBTMzAxwpWBmZon9wO2rvROdOBTMzAbnSsHMzJZwKJiZGeBKwczMEr5OwczMJlwpmJnZEg4FMzMDPHxkZmYJDx+ZmdmEQ8HMzJZwKJjZGrIO2Fd6nM6zIbmnYGZrQPngvy6fX57ngBiah4/MbJU1HfTrljkYhjQf/7IOBbMRajro07DMwTAUVwpmtkrqDvIHA3vzdYrHxdTBMDT3FMxslTQFwsH5OntpDwabpfmpFA5Y7R0ws9lqCoS9ZO9Xy/MOpnqYyWalCIU+X+0kbZZ0haQdkl5ZsfxBkj4r6Z8lfVXSqW3bdKVgNiJtgbA+X+8W4J5UB4QrhiEEsx4+krQOOBf4JWAnsFXSloi4PFntD4APRcRbJR0LXAA8uGm7DgWzESm/2y8Hwi35evekOhgcCEMZZPjoBGBHRFwFIOkDwGlAGgoBHJY/vjdwbdtGPXw0Q+9d7R0wq7F/ymU2K4MMHx0BXJN8vzOfl/oj4LmSdpJVCb/VtlGHwoy8F3jeau+ELby64aOiSrgnS6uEW/Jl7isMbapQ2CBpW/L14tJGVfODUmcA74qIjcCpwHslNR73PXw0Iw4EWwuaho+KEAAPH628qXoKuyLi+IblO4Ejk+83cvfhoRcCmwEi4h8lrQc2ANfXbdSVgtmINDWab+PulUK50exKYSiDDB9tBTZJeoikg4HTgS2ldb4L/CKApEeQFYY/bNqoKwWzEWmqFIpg2E92ZPCZRytp9o3miLhT0lnAhWT/Ze+IiO2SzgG2RcQW4HeAt0v67XwnXhAR5SGmJRwKZiPTFgxUzHMgDG2Yi9ci4gKyBnI67+zk8eXAk/ps06FgNjJVt6xID/7QHggOhlkL4PbV3olOHApmI1N3k7vi4E/y2IGwUubnNhcOBbMRarv7qe+QutIcCma2yqa9dbYNxaFgZqtsXzJt+zhOB8KQfOtsM1tj9lU8dhCsFA8fmZlZKuYjgh0KZmZDC+CO1d6JbhwKZmZDC+ZmrM6hYGa2EubkHuUOBTOzoblSMDOzCfcUzMxsCVcKZmYGZJWCewpmZjbhSsHMzAD3FMzMLOGzj8zMbAn3FMzMDMgCYW/rWmuCQ8HMbCW4UjAzM8A9BTMzK3EomJkZ4FNSzcws4SuazcxsCQ8fmZkZ4EazmZkl3FMwM7MlXCmYmRngRrOZmZW4UjAzM2CuegoHrPYOmJmNXnH2UZ+vDiRtlnSFpB2SXlmzzrMkXS5pu6S/adumKwUzs5Uw456CpHXAucAvATuBrZK2RMTlyTqbgN8DnhQRP5J0/7btOhTMzIY2zK2zTwB2RMRVAJI+AJwGXJ6s8yLg3Ij4EUBEXN+2UQ8fmZmthP09v9odAVyTfL8zn5c6BjhG0j9I+qKkzW0bdaVgZja06a5o3iBpW/L9eRFxXvK9an5S6kBgE3ASsBH4vKRHRcSP636oQ8HMbCX07ynsiojjG5bvBI5Mvt8IXFuxzhcj4g7gaklXkIXE1rqNevjIzGxoQdZT6PPVbiuwSdJDJB0MnA5sKa3zt8AvAEjaQDacdFXTRl0pmJkNbYArmiPiTklnARcC64B3RMR2SecA2yJiS77sKZIuJxvA+t2IuKFpuw4FM7OVMMAVzRFxAXBBad7ZyeMAXp5/deJQMDMbmu99ZGZmE0VPYQ44FMzMVoIrBTMzA/zJa31dE1F1EYbZmuPXqk1lju6SuiZCwcxs1FwpmJnZEu4pmJkZ4ErBzMwS7imYmdkSc1IpjP6GeJI+J+k2SXvyrytKy58t6TuSbpb0t5IOT5YdLumj+bLvSHp2w8+RpNdLuiH/eoMkn6lilST9b0nfl/Svkq6U9BvJsoMlnS/p25JC0kkVz3+spIvz1/R1kl5WWv4ySVfnr92vSzqmZj8+lfxt7JG0V9JlM/+FF91AH8c5hNGHQu6siDgk/3p4MVPSI4G3Ac8DfhK4BXhL8rxzya5D/EngOcBb8+dUeTHwdOA44KeApwL/eda/iI3GnwAPjojDgKcBr5P0uGT5JcBzgR+Un5jf7fLTZK/d+wJHA3+XLP8N4IXAfwAOIXst7qraiYg4JfnbOAT4AvDh5f96djez/5CdQSxKKNR5DvDxiLg4IvYArwZ+WdKhku4FPBN4dUTsiYhLyG5L+7yabZ0JvDEidkbE94A3Ai+oWlHSSZJ2SnqFpOvzd4xPl3Rq/q7xRkmvStY/QdK2/F3ldZL+bHb/BLYaImJ7RNxefJt/PSxftjci/jx/zVW9Z3w5cGFEvC8ibo+ImyLi6wCSDgD+EPjtiLg8Mt+KiBvb9knSg4ETgffWLc8rl1+XdI2kH0n6TUmPl/RVST+W9JfJ+kdL+ntJuyXtkvTBrv8+o1P0FPp8rZJFCYU/yV+U/1AqxR8JXFp8ExHfIqsMjsm/9kXElcn6l+bPqbJkWy3rAjwAWE/28XlnA28ne2f4OLI/zLMlPTRf983Am/N3lQ8DPtSwXZsTkt4i6RbgG8D3Kd3tssETgRslfSF/U/FxSQ/Kl23Mvx6VH7ivlvSaPCzaPB/4fERc3bLeE8g+qOXXgD8Hfh94Mtnr/VmSfj5f77VkFcy/yffpf3X8/cbHw0dryv8AHkp28D0P+Likh+XLDgF2l9bfDRzasqxKef3dwCENfYU7gD/OPxHpA8AGsgP/TRGxHdhONgxVrHu0pA151fLF2t/W5kZEvJTs9XQi8BHg9uZnTGwkq0xfBjwIuBp4f7IM4CnAo8k+YOUMsuGkNs8H3tVhvddGxG0R8XfAzcD7I+L6vEL+PPCYfL07gKOAB+brX9Jh2+PlUFgbIuJL+YH29oh4N/APwKn54j3AYaWnHAbc1LKsSnn9w4A9+f3Mq9wQEcV//a359Lpk+a1kQQPZH/QxwDckbZX01Jpt2pyJiH35wXIj8JKOT7sV+GhEbI2I24DXAD8j6d7c9Vp6Q0T8OCK+TdZ7OLV6UxlJP0tWvZ7f4eeXX6d1r9tXkH2O8D9J2i7pP3XY9jjN0fDRIp6SGtz1gdfbyRrDAOTDNfcAriRr9RwoaVNEfDNf5bj8OVWKbf1Th3X77XD288/IhwB+GThf0n0j4uZZbN/WhAPJewodfJWlH9BePBZwBdkQaN2bkTpnAh/Je2szERE/AF4Ek9D5jKSLI2LHrH7G3Jiji9dGXSlIuo+kkyWtl3SgpOcAP0f2EXUA7wP+o6QT88byOWR/GDflB9yPAOdIupekJwGnUdOEA94DvFzSEZIeCPwO3UrxLr/HcyXdLyL2Az/OZ8/JS8zKJN1f0umSDpG0TtLJZEM8FyXr3EPS+vzbg/PXcPFm5p3AMyT9tKSDyE6QuCSvDG4BPgi8Ij9hYiPZgfkTDfvzE8CvMqPXa7LdX81/PsCPmKtD4wB89tGacBDwOuCHZKfk/Rbw9Ii4ArIzQIDfJAuH68nGd1+aPP+lwE/ky94PvCR/DnmQpO+q3gZ8HLgM+BrwyXzeLGwGtuc/783A6fmwgc2nIBsq2kl2sPxT4L9FxMeSda4gG4o5guxNzK1k4/NExEXAq8heY9eTnZKaXkNzFtlw5rXAPwJ/A7wDKl+3kJ1KvRv47Mx+w8zjgS/lP28L8LIOTexxmqNGs+qHvM3MbBaOv49i24n9nqNP8OWIOH6YPaq3iD0FM7OVNUcDZw4FM7OV4Ftnm5kZ4ErBzMwSvnV2P0dK7navEf4M4mZ+ra4dc/dadaVgZmZAVinMSU9h7NcpmM2tdQNNh9x2+jOsZE6uU3AomK2gPgfXffn0jhlPh9hmMS32uc/vuhDm6N5HDgWzgVUd6MvTpoPtQTXT9R2mt1VMq+al0y7brtunLoFRnqb/RqM1R1c0OxTMBtAUBMs98BcH75s7TO9VMa2al067bLspQPoExsIEhEPBbHG1BUGfAOhyoD+0YnpTw7RpWTqt2naXAJk2KKoCYjQ8fGS2mOoCoUsQVFUATQFQHLx3V0zvnU8Pr5hWzUun967ZZl1gNAVFVUXRFBCjDgbfJdVssUwTCE1B0BYAVQf+G0vTH1ZMq+al0/I22gKjKSjqAmLhgmGg4SNJmyVdIWmHpFc2rPcr+edrt95gz6FgNgNdA6FuaKgqCLoGQHpQv18+fUDF9Af5NH38g5p1021VBUZbULQFRFPlMNpgmHEoSFoHnAucAhxL9kFcx1asdyjwX4EvddlNh4LZMvUJhLQyqBoaaguCugAoH+S/VzE9Ip9+DziyNK/qOVXBURUUdQGRVhBVQ0xNlUP533TuDdNTOAHYERFXRcRess96P61ivdcCbyD7J2/lUDBbpuUEQrkymCYI0oP7kcn0mprpkcB3WtapC44uAdGncugaDHNvuuGjDZK2JV8vLm31CLL/ssLOfN6EpMcAR0ZE7Sfvlfk2F2bLtJxASHsFaSCkB900CNKDdBoE3yH7WLZ0+lDgqtL0m/k+b8ofb6pY56qKbR1FdvQp/+wjyIKivK/3IxtqKn6ne7P07Ka0ckivk6jqNYymYujfPN7V8iE7Vfd+mtybK/9M9zcBL+jzQ10pmC3DcgKhbqioa0VQvKsvAqAcBN+smG5iaSBUrZNuI912uZLoUkE0VQ59K4a5NkyjeSfZf0dhI9lHsBYOBR4FfE7St4EnAlvams2uFMyWYTmBkAZB2jSuqwzKFUF6oE+nDyf7gOe6KR3WKablbReVRVMFUa4c0rOZypVDl4phFJXCMLfO3gpskvQQsn/+00k+qzsidgMbiu8lfQ747xGxrWmjrhTMlqE4WC0nENJx+XJlkI7vNwXCw6kOhGMrplXz0mm6jXTbaWWR7ku5B1E+i6mp59ClYhhNb2HGlUJE3AmcBVwIfB34UERsl3SOpKdNu5uuFMymVDd01DcQpq0M6t7dHwtc3jAFeDRwWYd1u1QQ01QObRVDXW9hbg106+yIuAC4oDTv7Jp1T+qyTYeC2ZTqho6GDISuAVAc9KumNCxLp12DYiWCYe6HkALYu9o70Y2Hj8ym1DR0NItAKDeN2wLh0VQHwnGladW847h7MKTbTH9W1RBTuTld1Yxua0KXh5JG12yek9tcuFIwm1Ld0NG0gXAN7UNFbRXBccClDdPH5vv+lfzxVzo8p62CaBpS6loxlE9XLTed575SYH723aFgNqVypVAcxMrXIUwbCF0qg6YgKA765SkNy9qCYppgKFcOTcHQNoQ0r4ozUueBQ8FsSuVKodxLSC9Mm1UgtFUG5YP748nOW0ynVMxLp20BUVc5tPUaqoJhmt7CPBrmjNRhOBTMptRn6OhGZhcI0wTBVuAJ3HVHtOJxMV1uQAwRDFVDSPNsFdsEvTgUzKbUNnRU1UvoGghtQ0VtQVA+6KdTGpY1BUVdQHQZUmrrMTQFw1gqhXnZd4eC2ZTqho66BML3mH0g1B3cfwb4QjKlYl4xbQuIWQdDelZSW29hXg6qVRwKZgugXCmkvYSq5nLdaafTBEJbZVA+2P8scEk+JXmcTtsCoq5yWE4wtDWdx1QpuKdgNnLlSqE8dFRuLs86EPoEQTqlYVnfgGgLhqYeQzkYxj6ENC/77VAwm0LV6ahdmst1gfBwpguEtiD4eeDvS1Mq5qXTrgHRNxiarmMoB8PYhpAGusvFIBwKZlNIT0dtGzoqB8JR3HVw7HKWUd9AqDvYF1OAfw9c1GHdoYKh6srnumBwpbCyHApmU+gzdFR1tlHdlcp9A6GtMigO/umUinnptGvlsJxgaOotpLfAKAfDvBxYy9xTMBu5uiZz1dBR117CtIHQJQguAp4MfCbf/+JxMe0bEEMGwxgrBZ99ZDZydZVC29DRSgdC+eD/5Hz/y/OaAmItBcO8HFiruKdgNmJNTeaqoaO0lzDrQOgSBJ8BTib7NBaSx8W0b0AsNxjqegvpB/OU7546z8Gwn7m5c7ZDwWwaaaVQbjJXDR2lvYTy2UZDBkL54H9yvv/leU0BMatgcKUwHxwKZlNIK4W6JnOfoaNZBkLdwf5U7vqIruJxMW0LiGmDwUNIGfcUFkzx2bFjuOe7dVN1OmrRZK4aOurSSyhuHzHLQCgf/E/N9788rykgpg2Grr2FdAhpjIFQcKUwcunn8+6rmeeAGK+iUqhqMlcNHaXXJXTpJQwVCEWl8FTgEw3rTBsMXXsLi1gpuKcwQnUH/apKwQExbuVKoW3oKL0uoWsvobhYrE8g1B3kixB4ar7/nyjN6xIQbcEwq95C2mweQyCAr2gepS4H+7Z1HAzjUa4Uyk3mrkNHxWcTtPUSilNC+wRC+aBfTGlYVhcQdcEwq97CQxl3pQDzs98OhQ66HOzpsI6DYTzaKoV06Ogq7n7WUdXQUZ9ewrSB0FQpDBUMyx1COpT5DwQ3mkekz0G+67oOhvlXVymkTeams46aegnF0NGsA+E04GP5/hePi+msg2HaIaSxNpt9m4uR6HOQp8e6Dob5V1cppE3m4l5H0w4ddQmEk2kOhPLB/7R8/8vzmgKia/O5rbfQ5UK2ukrh5vzfep7/ZtxTmHN9Du7FO4Dy/VkcDONVrhTS01GrmszTDh2VewnFRWVVZxl1CYS6SqFvMAw9hPQAxlcpzMt+OxRqdAmCNASomNclIGw+lSuF9HTUH9DcZL6M6gvWljN01DUQulQKswqGuiGkaZvN81wpePhoBPoEQlWl4GAYt6JSSC9cK5+OWm4y112wNu3QUTre3zUQikrhGcBHO6w7bW+h6xDSIp2WOi/77VCo0DcQ2ioFB8P4FK+F9MK14nTUukqh7vqE4nOP+w4d1fUSmg7yz8j3/6O0B8OshpC6Npvbrmye90phXvbdoVCh6eBdd9CH7EVbvve7A2Gcllsp1DWZhxw6KkKA5PFQwTBNs7nqtNQxVQpuNM+xvpXC+vx5t+FgWBRNlULdhWtNlUJ56KhcKSx36CgNAVgaCF2DoW0IqW+zue9pqfNcKexnfnoKB6z2Dqw10wTCbSwNhPU167ad0mrzo6lS6HLPo75N5uUMHZUP/lWB8Iya56bbbuspFCFWhFq5UiiffVSEZLmnkH6uQnr2UfE3NK/29fzqQtJmSVdI2iHplRXLXy7pcklflfT/JB3Vtk2HgpmtmHk+qC9H0VOYZShIWgecC5xClq9nSDq2tNo/A8dHxE8B5wNvaNuuQ6GkTz8hrQzSx3XVgoePxqNt+Ci9RiEdPupyq+x0+Gi5/YSqiqBqXl21UHd1c1GxFBVM0QwvDx8VFVC50Vw+JTX9FLb0lNT009fmuZ8A2RBSn68OTgB2RMRVEbEX+AB3nXUMQER8NiJuyb/9IrCxbaPuKZT0GT5KQwDuHgzuK4xX2/BRejVzOnzU9xqFtjOPyremKB/Uqw7+VMyr6yvM+lqFRW00T3nr7A2StiXfnxcR5yXfFy+1wk6yf/I6LwQ+1fZDHQolfc88KkIAHAiLpGujuemW2V0qhVk0mssHfyrmudE8rClvnb0rIo5vWK6aH3X3FaXnAseT/Tc1cihUmCYYwIGwSLo2mruekjpko7kcDOBTUlfDAPu9kyxDCxuBa8srSXoy8PvAz0fE7W0bdShUmPbitbpAcDCMz5AXr1VVCssZQioHAwwbCL547e4Gus3FVmCTpIeQvQ85HXh2uoKkxwBvAzZHxPVdNupQqDDtxWsOhMUx5MVrQ97molwp+DYXK2OIK5oj4k5JZ5H9V6wD3hER2yWdA2yLiC3A/wQOAT4sCeC7EfG0pu06FGr0DQYq5jkQxmu5lULabE4rhXQIaYjbXBSnpvQJBN8QbzaGuKI5Ii7gro/eLuadnTx+ct9tOhRqtN3u+iCWBgTcPQi63kLb5k/x/9t06+wun7rWp9k8iyGkj+X7P3Qg+NbZS/neRyPQ5WCeBgT0CwIHwnzrcuvs9LTUultn1zWbp/mQna7BQMOyLoHgD9npz7fOHom+B/eu6zoQ5l+5Umj6OM5ypdDWbO4zhNTUWyiCoRwQUB0ETYHgj+NcHlcKI9L1IN9n3Xl5cVi9cqXQ9HGcdUNIaaUwiyGkumAoBwQ0B8FyA2G5Q0dVlcI8B0LBd0kdkVlXCjb/6iqFtjulTjOENMtg+ES+/3XLZh0IbUNHXZvMrhRWjkOho1lVCjYOdZVC2wVsfYeQ6noLXYOhKiComPeJ0nO6BEIxbQuErkNHYz0dFdxTGK19ybRPpZA+18ahrVKoaja3Xa/Q51PYmoKhLSComFcVBF0CoW8voe/Q0VgqBZif/XYoTKkcEOljB8H4lSuF8mmpXe+B1NRbaBtCamo+VwVDMaVhWZezjJbbS+h7JfMYAsHDRwtmX81jG69ypVA+LbXPEFLbXVOXGwzlgID2ymC5gdC1l7BIlYIbzWYjlt4Hq6gU2oaQ2noLVUNI0wZDU0BAexDMIhD69hLGHAjuKZiNXHo1e1EplJvN5SGkqgvZuvQWpgmGpoCA9iCYdSA0XbBWDB2lH8N5OEs/YGfeg8HDR2Yjl1YKdc3m8hBS04Vsyw2G8umqTQEB7UFQPu10loGwiJUCzM++OxTMppBWClXN5mmGkJYTDE2VQzkgoD0IqiqD1Q6E9L5i88bDR2YjV1UppM3mqiGkIYOhqXIoBwS0B0FVZbASgZAOHY2pUvDwkdnI1VUKTUNIxb2QytctLDcY2iqH8hSqA+Ci0nPLlcGsAqGql7AIlYLPPjIbsfQOuWmlUDSbm4aQ0usWlhsMXSqH8hSqA6DLUNEsAqHP0NEYmsyFedl3h4LZFIoLFeuazeUL2ep6C1VnI/UJhi6VQ3kGBgzeAAAGi0lEQVRKw7K2ymDaQOjTS6gbOpr3SsE9BbORaxtCSi9k69N07hMMbZVDVVDA3XsEdUFQrgz6BELfXsJYh47APQWzhVDXbD6U6t5C8YlsswqGtsohDYg0KODuAVAXBFWVwawCIb0uYexDR+CegtnolSuFdAipqbfQJRjaegxdKoeqgHhCvu/leXVBUFcZtA0ZdQmERRk6AlcKZguhXCnUDSE1BUNxNlI5GNqaz10qh7qggOYA6FIZLKeH0CcQirCd50AA9xTMFkJaKZSHkMrBkF63kAZDeppqn2BoqxzSgCgHBVQHQFMQlCuDvoHQdvppORDKQ0fzXinA/Oy7Q8FsSulpqeUhpKbeQtdgqDtdta1yqAuIYkrDsrogqKoM6k47raoQqgKhay9hDIHg6xTMFkDVaanlIaSq3kJdMFSdrtq1Cd0WEOl0XcOyuiDo2kyuuw6hSyAUIVo3dDTPn14YwN7V3omOHApmUypXCm1DSG3BUHcdQ1sTui4gqoLi0nzfy/OqAqCtMphVICzC0BG4UjAbvXKlMEQwPJTmIaWmgKgLCirmXVp6TtcgSIeKin0dKhDmvVKYl313KJhNqapSqOotLCcYyr2GcuXQJSDKUzqsUxcETZVBsa+zDoQxVAoOBbMFkL57LVcKfYIhPV31ASy98rnca2iqHKoCoqoXQcW8qudUBUFTZVDsaxEExe+S/o5dAiHtJYypUvApqWYLoGkIqWswTFs51AVEW1BQMa8pAOqCoG9lsKiBUHBPwWwB1A0hdQmGrkNKdZVDXUCkQ0xVUxqWVQ0N1QVBuTLoOlS0iIHg4SOzBVE3hNQlGNLrGJqGlOoqhzQgqoaYytPiIE/yuHzgL0+PKv2MNAjKlUGXoaL0tNO2QBhbMMzL/h+w2jtgNu+WEwx9h5TSnkNaQTQFRVpRHFV6XFUBVAVAVUWQ9gy6DhUtaiAUPYU+X11I2izpCkk7JL2yYvk9JH0wX/4lSQ9u26ZDwWyZioPWNMHQZUjpcKoDom9QFNOqeX0CoCoIin3sMlS0aIEAdw0f9flqI2kdcC5wCtn5AmdIOra02guBH0XE0cCbgNe3bdehYLZM6cGrLRjS6xjaKoe2gOgSFFWBUT7glw/8bQHQJwiqKoP0tNMugbCu+3/Fmra/51cHJwA7IuKqiNgLfAA4rbTOacC788fnA78oSU0bdU/BbAa6BkNd5TBNQDQFRVNgUDGv7sBfFQDLCYKmyqAuEMZUKcxYUdgVdnLX3dHvtk5E3ClpN3BfYFfdRh0KZjMyTTD0DYi0OV2elg/a6cE8ndKwrGobVT+r3DTuGgRNlcFYAwFgP1x4E2zo+bT1krYl358XEecl31e944/S913WWcKhYDZDdcHQJyCqhpiagqIcGE3BcVO+n3XLygf8qgN/UwBUDQ1NUxmMKRAAImLzAJvdSdYOKmwErq1ZZ6ekA7kr62u5p2A2Y/uSaVNAzCoo+gTHvVqWdTngLzcAulQGYwqEAW0FNkl6iKSDgdOBLaV1tgBn5o9/BbgoIlwpmK2WpoCYVVDUDUXVTemwTtWQzxAB4CCYXt4jOAu4kOyf8h0RsV3SOcC2iNgC/DXwXkk7yCqE09u2q5bQWBFHSqu/EwbANRGNZyYsupV6rdYdPGcxZcBtr+TQj1+rw3ClYLYGld89z3q6Uj/D5o97CmZmNrEmho/MzGxtcKVgZmYTDgUzM5twKJiZ2YRDwczMJhwKZmY24VAwM7MJh4KZmU04FMzMbMKhYGZmEw4FMzObcCiYmdmEQ8HMzCYcCmZmNuFQMDOzCYeCmdkMSXqRpH/Jv/Ynj/+sYt23SXpSad6e5PGpkr4p6UErse/gz1MwMxuEpCOAL0TEUQ3r/AvwuIjYl8zbExGHSPpF4DzgKRHxreH3OOOP4zQzG8ajgMvqFkp6BHBlGgjJshOBtwOnrmQggEPBzGwojwa+1rD8FODTFfPvAXwMOCkivjHEjjVxT8HMbBiNlQJwMtWhcAfwBeCFQ+xUG4eCmdkwaisFSfcE7hMR11Ys3g88C3i8pFcNuH+VPHxkZjZjkg4ANgF1wz+/AHy27vkRcYukpwKfl3RdRPz1ALtZyaFgZjZ7RwM7I+L2muWnAOc3bSAibpS0GbhY0q6I+Nisd7KKT0k1M1thkr4CPCEi7ljtfSlzKJiZ2YQbzWZmNuFQMDOzCYeCmZlNOBTMzGzCoWBmZhMOBTMzm3AomJnZhEPBzMwm/j9F5XFEaIAC3AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cbd9312e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# plate size, mm\n",
    "w = h = 10.\n",
    "# intervals in x-, y- directions, mm\n",
    "dx = dy = 0.1\n",
    "# Thermal diffusivity of steel, mm2.s-1\n",
    "D = 0.75\n",
    "\n",
    "Tcool, Thot = 0., 1.\n",
    "\n",
    "nx, ny = int(w/dx), int(h/dy)\n",
    "\n",
    "dx2, dy2 = dx*dx, dy*dy\n",
    "dt = dx2 * dy2 / (2 * D * (dx2 + dy2))\n",
    "\n",
    "r0 = np.zeros((nx-2, ny-2))\n",
    "u0 = Kr * np.ones((nx, ny))\n",
    "u = np.empty((nx, ny))\n",
    "\n",
    "# Initial conditions - ring of inner radius r, width dr centred at (cx,cy) (mm)\n",
    "r, cx, cy = 0.02, 5, 5\n",
    "r2 = r**2\n",
    "for i in range(nx):\n",
    "    for j in range(ny):\n",
    "        p2 = (i*dx-cx)**2 + (j*dy-cy)**2\n",
    "        if p2 < r2:\n",
    "            u0[i,j] = Thot\n",
    "\n",
    "def do_timestep(u0, u, r0):\n",
    "    # Propagate with forward-difference in time, central-difference in space\n",
    "    \n",
    "    u[1:-1, 1:-1] = u0[1:-1, 1:-1] + D * dt * (\n",
    "          (u0[2:, 1:-1] - 2*u0[1:-1, 1:-1] + u0[:-2, 1:-1])/dx2\n",
    "          + (u0[1:-1, 2:] - 2*u0[1:-1, 1:-1] + u0[1:-1, :-2])/dy2 ) + A*(\n",
    "        u0[1:-1, 1:-1]-Kr)*(u0[1:-1, 1:-1]-Ktheta)*(u0[1:-1, 1:-1]-Km)*(u0[1:-1, 1:-1]+0.1) - r0*u[1:-1, 1:-1]\n",
    "    \n",
    "    u[:,0:2] = Kr\n",
    "    u[:,ny-2:ny] = Kr\n",
    "    u[0:2,:] = Kr\n",
    "    u[nx-2:nx,:] = Kr\n",
    "    \n",
    "    r0 = r0 + dt*B*((u0[1:-1, 1:-1] - Kr) - C*r0)\n",
    "\n",
    "    u0 = u.copy()\n",
    "    \n",
    "    return u0, u, r0\n",
    "\n",
    "# Number of timesteps\n",
    "nsteps = 1000\n",
    "# Output 4 figures at these timesteps\n",
    "mfig = [0, 50, 150, 950]\n",
    "fignum = 0\n",
    "fig = plt.figure()\n",
    "for m in range(nsteps):\n",
    "    u0, u, r0 = do_timestep(u0, u, r0)\n",
    "    if m in mfig:\n",
    "        fignum += 1\n",
    "        print(m, fignum)\n",
    "        ax = fig.add_subplot(220 + fignum)\n",
    "        im = ax.imshow(u.copy(), cmap=plt.get_cmap('hot'), vmin=Tcool,vmax=1)\n",
    "        ax.set_axis_off()\n",
    "        ax.set_title('{:.1f} ms'.format(m*dt*1000))\n",
    "fig.subplots_adjust(right=0.85)\n",
    "cbar_ax = fig.add_axes([0.9, 0.15, 0.03, 0.7])\n",
    "cbar_ax.set_xlabel('$T$ / K', labelpad=20)\n",
    "fig.colorbar(im, cax=cbar_ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}