pribec2021.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fd116332",
   "metadata": {},
   "source": [
    "### A discrete reactive collision scheme for the lattice Boltzmann method\n",
    "Implemented based on the paper of Ivan Pribec et al. (2021) published in the Journal of Molecular Liquids."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a433469",
   "metadata": {},
   "source": [
    "The code solves the reaction-diffusion-advection (RDA) equation using the lattice Boltzmann method (LBM) with D2Q5 lattice. It simulates two reactive components $A$ and $B$ which undergo a reaction creating another (inert) component $C$. It is a irreversible reaction.\n",
    "\n",
    "$A + B \\longrightarrow C$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b22293c",
   "metadata": {},
   "source": [
    "Diffusion coefficient $D_s = \\frac{1}{3} (\\tau_s - \\frac{1}{2})$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4d53132",
   "metadata": {},
   "source": [
    "Source term from the reaction component\n",
    "$R_s^\\prime(\\mathbf{x}, i, t) = \\sum_{j,k} K_{AB,s}(i,j,k) \\sigma_{AB}(g_{j,k}) g_{j,k} f_A(\\mathbf{x},j,t) f_B(\\mathbf{x}, k, t)$\n",
    "\n",
    "$i = 1\\ldots Q$, $j \\in 1\\ldots Q$, $k \\in 1\\ldots Q$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdc69871",
   "metadata": {},
   "source": [
    "Matrices $K_{AB,A} = \\{ \\}$ and $K_{AB,B} = \\{ \\}$, while matrix $K_{AB,C}$ is the most complicated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "0281c07b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0YAAAF0CAYAAAAHCSmXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAA9hAAAPYQGoP6dpAABiNklEQVR4nO3dd3xT1fsH8E+6Usooo3RBKUPZq4LQgsiSDSKIsmSLIqIMEQH1CzhAHIggQ5QpiijLgUyZSsEiBYGyZK+yactqaXt+fzy/NA1NS9omuRmf9+uV1705GfdJmubcJ/fc5+iUUgpERERERERuzEPrAIiIiIiIiLTGxIiIiIiIiNweEyMiIiIiInJ7TIyIiIiIiMjtMTEiIiIiIiK3x8SIiIiIiIjcHhMjIiIiIiJye0yMiIiIiIjI7TExIiIiIiIit8fEiLJYsGABdDpdxsXLywulS5dGv379cP78eatv786dOxg/fjy2bNmS68eeOHECQ4YMQcWKFVGgQAH4+fmhWrVqeOedd/IU6/jx46HT6UzamjRpgiZNmpi06XQ6jB8/Pt/P7azMvSe2fP6cPiOGz+upU6dsFg8RuRfD98ru3but8nwTJ07EqlWrsrTHxcVh/PjxZr+/+vbti7Jly1pl+1o5deoUdDodFixYYLfn37FjB8aPH4+bN29mub+t+y5yfl5aB0COa/78+ahcuTLu3r2Lbdu2YdKkSdi6dSv279+PggULWm07d+7cwYQJEwAgV19Yv/32G7p164aAgAAMGTIEERER0Ol02L9/P+bNm4fVq1cjNjbWanFmFh0djdKlS+fqMS+++CJat25tk3hczcyZM02u5/QZadeuHaKjoxESEmKv8IiIcmXixIno0qULnnnmGZP2uLg4TJgwAU2aNMmSBL377rsYOnSo/YJ0QiEhIYiOjkaFChUy2nbs2IEJEyagb9++KFq0qMn9H+xbiB7ExIiyVb16ddStWxcA0LRpU6SlpeH999/HqlWr0LNnT01jO3nyJLp164aKFSti8+bN8Pf3z7itWbNmeP3117Fy5UqbbT8yMjLXjyldunSukyl3VbVqVYvvW7JkSZQsWdKG0RAR2V/mnX0yT6/X56o/zk3fQu6JQ+nIYoYvn9OnTwMA7t27hzFjxqBcuXLw8fFBqVKl8Oqrr2Y5fL1p0yY0adIEJUqUQIECBVCmTBk8++yzuHPnDk6dOpWxUzthwoSM4Xt9+/bNMZYpU6bg9u3bmDlzpklSZKDT6dC5c2eTtnnz5qFWrVrw9fVF8eLF0alTJxw6dChP78WDQ+nu3LmDkSNHoly5chnPX7duXSxZsiTjPuaG0qWnp+Pjjz9G5cqVodfrERgYiN69e+PcuXMm92vSpAmqV6+OmJgYNGrUCH5+fihfvjw++ugjpKen5xhrREQEGjVqlKU9LS0NpUqVMnmfUlJS8MEHH2TEU7JkSfTr1w9Xrlx56Hty/fp1DB48GKVKlYKPjw/Kly+Pt99+G8nJyVle8/Tp01G7dm0UKFAARYsWRWRkJH755ReT12s4MvSwz0h2Q+k2btyI5s2bo0iRIvDz80PDhg3xxx9/mNznypUreOmllxAWFpbxehs2bIiNGzc+9PUSkXu7d+8e3njjDdSuXRv+/v4oXrw4oqKi8PPPP5vcT6fT4fbt21i4cGHG91eTJk2wYMECPPfccwDkx0fDbYZhYeaG0lny/QkAS5cuRVRUFAoWLIhChQqhVatWDx1BsW/fPuh0OsydOzfLbWvWrIFOpzPZzrFjx9CjRw8EBgZCr9ejSpUqmDFjhkXv3Z9//onmzZujcOHC8PPzQ4MGDbB69eos9zt//nzGd7SPjw9CQ0PRpUsXXLp0CUDWoXTjx4/Hm2++CQAoV65cxntqGIZtbiidpf1eTvsy5Dp4xIgs9t9//wGQX+iVUnjmmWfwxx9/YMyYMWjUqBH+/fdfjBs3DtHR0YiOjoZer8epU6fQrl07NGrUCPPmzUPRokVx/vx5rF27FikpKQgJCcHatWvRunVrDBgwAC+++GLGNnKyfv16BAUFWfxL0aRJkzB27Fh0794dkyZNwrVr1zB+/HhERUUhJiYGjz76aL7emxEjRuDbb7/FBx98gIiICNy+fRsHDhzAtWvXcnzcK6+8gjlz5mDIkCFo3749Tp06hXfffRdbtmzBnj17EBAQkHHf+Ph49OzZE2+88QbGjRuHlStXYsyYMQgNDUXv3r2z3Ua/fv0wdOhQHDt2zOR1rl+/HhcuXEC/fv0ASIfbsWNHbN++HaNGjUKDBg1w+vRpjBs3Dk2aNMHu3btRoEABs9u4d+8emjZtiuPHj2PChAmoWbMmtm/fjkmTJmHv3r0mHV7fvn2xePFiDBgwAO+99x58fHywZ8+ebM8RystnZPHixejduzc6duyIhQsXwtvbG1999RVatWqFdevWoXnz5gCAXr16Yc+ePfjwww9RsWJF3Lx5E3v27Hno342IKDk5GdevX8fIkSNRqlQppKSkYOPGjejcuTPmz5+f8b0cHR2NZs2aoWnTpnj33XcBAEWKFEHJkiUxceJEjB07FjNmzMBjjz0GIOcjRZZ8f06cOBHvvPMO+vXrh3feeQcpKSn45JNP0KhRI/z999/ZHjWpVasWIiIiMH/+fAwYMMDktgULFiAwMBBt27YFIEMAGzRogDJlyuCzzz5DcHAw1q1bh9dffx1Xr17FuHHjsn0NW7duRYsWLVCzZk3MnTsXer0eM2fORIcOHbBkyRJ07doVgCRFjz/+OO7fv4+xY8eiZs2auHbtGtatW4cbN24gKCgoy3O/+OKLuH79OqZPn44VK1ZkDLHO7jVb2u89bF/Gz88v29dLTkYRPWD+/PkKgNq5c6e6f/++SkpKUr/99psqWbKkKly4sIqPj1dr165VANTHH39s8tilS5cqAGrOnDlKKaWWLVumAKi9e/dmu70rV64oAGrcuHEWx+jr66siIyMtuu+NGzdUgQIFVNu2bU3az5w5o/R6verRo0dG27hx49SD/xaNGzdWjRs3Nml7MN7q1aurZ555Jsc4HnzuQ4cOKQBq8ODBJvfbtWuXAqDGjh1rEgMAtWvXLpP7Vq1aVbVq1SrH7V69elX5+PiYPJ9SSj3//PMqKChI3b9/Xyml1JIlSxQAtXz5cpP7xcTEKABq5syZJvFkfk9mz56tAKgff/zR5LGTJ09WANT69euVUkpt27ZNAVBvv/12jjE/+Pw5fUYMn9eTJ08qpZS6ffu2Kl68uOrQoYPJ/dLS0lStWrVUvXr1MtoKFSqkhg0blmMsROR+DN8rMTExFj8mNTVV3b9/Xw0YMEBFRESY3FawYEHVp0+fLI/56aefFAC1efPmLLf16dNHhYeHZ1y35PvzzJkzysvLS7322msm7UlJSSo4OFg9//zzOb6GadOmKQDqyJEjGW3Xr19Xer1evfHGGxltrVq1UqVLl1YJCQkmjx8yZIjy9fVV169fV0opdfLkSQVAzZ8/P+M+kZGRKjAwUCUlJWW0paamqurVq6vSpUur9PR0pZRS/fv3V97e3iouLi7beM09/yeffGLSJ2T2YN9iab9nyb4MuQYOpaNsRUZGwtvbG4ULF0b79u0RHByMNWvWICgoCJs2bQKALEPennvuORQsWDBjyFLt2rXh4+ODl156CQsXLsSJEydyFUNqaqrJRSmV69cRHR2Nu3fvZok1LCwMzZo1yzK8Ki/q1auHNWvWYPTo0diyZQvu3r370Mds3rwZQNb3sF69eqhSpUqWuIKDg1GvXj2Ttpo1a2YMbcxOiRIl0KFDByxcuDBj2N2NGzfw888/o3fv3vDykgPHv/32G4oWLYoOHTqYvOe1a9dGcHBwjlUDN23ahIIFC6JLly4m7YbXZngta9asAQC8+uqrOcacHzt27MD169fRp08fk9eRnp6O1q1bIyYmBrdv3wYg7/WCBQvwwQcfYOfOnbh//77N4iIi1/PTTz+hYcOGKFSoELy8vODt7Y25c+fmeZh2Tiz5/ly3bh1SU1PRu3dvk+8/X19fNG7c+KHVX3v27Am9Xm9S5W3JkiVITk7OGF1w7949/PHHH+jUqRP8/PxMttO2bVvcu3cPO3fuNPv8t2/fxq5du9ClSxcUKlQoo93T0xO9evXCuXPncOTIkYzX27RpU1SpUsWStydPLO338rsvQ86DiRFla9GiRYiJiUFsbCwuXLiAf//9Fw0bNgQAXLt2DV5eXlmGM+l0OgQHB2cMRapQoQI2btyIwMBAvPrqq6hQoQIqVKiAL7744qHbP3XqFLy9vU0uW7duBQCUKVMGJ0+etOh1GGIxV7UsNDTUKsOmpk2bhrfeegurVq1C06ZNUbx4cTzzzDM4duyY1eIqUaJElvvp9XqLkrD+/fvj/Pnz2LBhAwBjR5c5Kbt06RJu3rwJHx+fLO97fHw8rl69muNrCQ4OznIOVWBgILy8vDJey5UrV+Dp6Yng4OCHxpxXhrHnXbp0yfI6Jk+eDKUUrl+/DkDG4ffp0wfffPMNoqKiULx4cfTu3Rvx8fE2i4+IXMOKFSvw/PPPo1SpUli8eDGio6MRExOD/v374969e1bfniXfn4bvv8cffzzL99/SpUtz/B4HgOLFi+Ppp5/GokWLkJaWBkCG0dWrVw/VqlUDIN/3qampmD59epZtGIbaZbedGzduQCmVbb9neH7D67V1wSJL+7387MuQc+E5RpStKlWqZFSle1CJEiWQmpqKK1eumCRHSinEx8fj8ccfz2hr1KgRGjVqhLS0NOzevRvTp0/HsGHDEBQUhG7dumW7/dDQUMTExJi0VapUCQDQqlUrTJ8+HTt37nzoeUaGhOLixYtZbrtw4YLJeTx5VbBgQUyYMAETJkzApUuXMo4edejQAYcPH35oXA9++VsrLoNWrVohNDQU8+fPR6tWrTB//nzUr1/fZNx1QEAASpQogbVr15p9jsKFC2f7/CVKlMCuXbuglDJJji5fvozU1NSM11KyZEmkpaUhPj7eZuW1DduaPn16tp8Nw9j0gIAATJ06FVOnTsWZM2fwyy+/YPTo0bh8+XK27wMRESDnMpYrVw5Lly41+d57sOCMtVjy/Wn4/lu2bBnCw8PztJ1+/frhp59+woYNG1CmTBnExMRg1qxZGbcXK1Ys4whPdkevypUrZ7a9WLFi8PDwyLY/zvwaSpYsmaUQkbXlpt/L674MORceMaI8MZy8vnjxYpP25cuX4/bt2xm3Z+bp6Yn69etnVK3Zs2cPADnqASDLkQ8fHx/UrVvX5GL4kho+fDgKFiyIwYMHIyEhIcu2lFIZ5bqjoqJQoECBLLGeO3cOmzZtMhtrfgQFBaFv377o3r07jhw5km3FmmbNmgHI+h7GxMTg0KFDVo3L0ImtWrUK27dvx+7du9G/f3+T+7Rv3x7Xrl1DWlpalve9bt26GUmpOc2bN8etW7eyTGC4aNGijNsBoE2bNgBg0slaIrvPiDkNGzZE0aJFERcXZ/Z11K1bFz4+PlkeV6ZMGQwZMgQtWrTI+GwSEWVHp9PBx8fHJCmKj4/PUpUOyP7ofm6+2yz5/mzVqhW8vLxw/PjxbL//HqZly5YoVaoU5s+fj/nz58PX1xfdu3fPuN3Pzw9NmzZFbGwsatasaXYb5kY4APIjYv369bFixQqT15yeno7FixejdOnSqFixYsbr3bx5c8bQOkvl5j3NS7+X3b4MuQYeMaI8adGiBVq1aoW33noLiYmJaNiwYUZVuoiICPTq1QsAMHv2bGzatAnt2rVDmTJlcO/ePcybNw8A8NRTTwGQX2TCw8Px888/o3nz5ihevDgCAgJynPG7XLly+OGHH9C1a1fUrl07Y4JXQKrlzJs3D0opdOrUCUWLFsW7776LsWPHonfv3ujevTuuXbuGCRMmwNfXN8fqOZaqX78+2rdvj5o1a6JYsWI4dOgQvv32W0RFRWVbraZSpUp46aWXMH36dHh4eKBNmzYZVenCwsIwfPjwfMeVWf/+/TF58mT06NEDBQoUyKj8Y9CtWzd89913aNu2LYYOHYp69erB29sb586dw+bNm9GxY0d06tTJ7HP37t0bM2bMQJ8+fXDq1CnUqFEDf/75JyZOnIi2bdtm/K0bNWqEXr164YMPPsClS5fQvn176PV6xMbGws/PD6+99prZ58/NZ6RQoUKYPn06+vTpg+vXr6NLly4IDAzElStXsG/fPly5cgWzZs1CQkICmjZtih49eqBy5cooXLgwYmJisHbt2iyl3onIPW3atMlsxcy2bduiffv2WLFiBQYPHowuXbrg7NmzeP/99xESEpJlGHWNGjWwZcsW/PrrrwgJCUHhwoVRqVIlVK9eHQAwZ84cFC5cGL6+vihXrpzZxMKS78+yZcvivffew9tvv40TJ06gdevWKFasGC5duoS///47Y3RDTjw9PdG7d29MmTIFRYoUQefOnbNMi/HFF1/giSeeQKNGjfDKK6+gbNmySEpKwn///Ydff/014zxkcyZNmoQWLVqgadOmGDlyJHx8fDBz5kwcOHAAS5YsyUg033vvPaxZswZPPvkkxo4dixo1auDmzZtYu3YtRowYgcqVK5t9/ho1amTE2KdPH3h7e6NSpUpmRz1Y2u9Zsi9DLkLDwg/koCytxnP37l311ltvqfDwcOXt7a1CQkLUK6+8om7cuJFxn+joaNWpUycVHh6u9Hq9KlGihGrcuLH65ZdfTJ5r48aNKiIiQun1egXAbPUec44fP64GDx6sHnnkEaXX61WBAgVU1apV1YgRI7JUpPnmm29UzZo1lY+Pj/L391cdO3ZUBw8eNLlPXqvSjR49WtWtW1cVK1ZM6fV6Vb58eTV8+HB19erVHJ87LS1NTZ48WVWsWFF5e3urgIAA9cILL6izZ89miaFatWpZXv+DVYsepkGDBgqA6tmzp9nb79+/rz799FNVq1Yt5evrqwoVKqQqV66sXn75ZXXs2DGTeB58T65du6YGDRqkQkJClJeXlwoPD1djxoxR9+7dy/KaP//8c1W9evWMv0VUVJT69ddfc3z+7D4jD1alM9i6datq166dKl68uPL29lalSpVS7dq1Uz/99JNSSql79+6pQYMGqZo1a6oiRYqoAgUKqEqVKqlx48ap27dvW/yeEpHrMXyvZHcxfN989NFHqmzZskqv16sqVaqor7/+2ux3/d69e1XDhg2Vn5+fAmDy/TZ16lRVrlw55enpaVJhzdz3uyXfn0optWrVKtW0aVNVpEgRpdfrVXh4uOrSpYvauHGjRa//6NGjGa91w4YNZu9z8uRJ1b9/f1WqVCnl7e2tSpYsqRo0aKA++OADk/vggapxSim1fft21axZM1WwYEFVoEABFRkZmeU1KKXU2bNnVf/+/VVwcLDy9vZWoaGh6vnnn1eXLl3K8fnHjBmjQkNDlYeHh0nVP3N9iyX9nqX7MuT8dErlocwXERERERGRC+E5RkRERERE5PaYGBERERERkdtjYkRERERERG6PiREREREREbk9JkZEREREROT2mBgREREREZHbc7kJXtPT03HhwgUULlzYZDZqIiKyPaUUkpKSEBoaCg8P/vZmwL6JiEgbuemXXC4xunDhAsLCwrQOg4jIrZ09exalS5fWOgyHwb6JiEhblvRLLpcYFS5cGIC8+CJFimgcDRGRe0lMTERYWFjGdzEJ9k1ERNrITb/kcomRYYhCkSJF2PkQEWmEw8VMsW8iItKWJf0SB4ATEREREZHbY2JERERERERuj4kRERERERG5PSZGRERERETk9pgYERERERGR22NiREREREREbo+JERERERERuT0mRkRERERE5PaYGBERERERkdtjYkRERERERG6PiREREREREbk9JkZEjmTnTkCvB6ZN0zoSIiIi8cEH0jfFxWkdCZFNMTEiciRz5wIpKcDs2YBSWkdDRETuLj0d+Oor6Zvmz9c6GiKbYmJE5Chu3wa+/17WDx0C9uzRNh4iIqItW4Bz52R93jzg/n1NwyGyJSZGRI5i+XLgzh1Z9/KSDoiIiEhL8+ZJnwQA168Da9ZoGw+RDTExInIU33wDePz/v2RqKvDtt8C9e9rGRERE7ishAfjpJ+mTAMDTU/oqIhfFxIjIERw/DmzfLmO5DZKSgJ9/1i4mIiJyb0uXmg6dS0sDfv8duHRJu5iIbMimidG2bdvQoUMHhIaGQqfTYdWqVQ99zNatW1GnTh34+vqifPnymD17ti1DJHIMCxbIL3GZeXoCX3+tSThErox9E5GF5szJ2qYUsHix/WMhsgObJka3b99GrVq18OWXX1p0/5MnT6Jt27Zo1KgRYmNjMXbsWLz++utYvny5LcMk0lZamlSjS0vL2r5pE3DmjDZxEbko9k1EFoiLA/75J2uF1PR0SZhYOZVckJctn7xNmzZo06aNxfefPXs2ypQpg6lTpwIAqlSpgt27d+PTTz/Fs88+a6MoiTT2xx/AxYvmb9PpgEWLgHfesW9MRC6MfRORBebPl5ELD/5oBwBHjwIxMUC9evaPi8iGHOoco+joaLRs2dKkrVWrVti9ezfuZ1MeMjk5GYmJiSYXIqeyaFH2t6WnszodkcbYN5HbUUqGeJtLigxy6ruInJRNjxjlVnx8PIKCgkzagoKCkJqaiqtXryIkJCTLYyZNmoQJEybYK0Qi63vhBeMcEUlJxvmLGjeW5dNPaxMXEQFg30RuSKcDXn9dRjQAwNatsqxfH/D1lfWuXbWJjciGHCoxAgCdTmdyXf3/GNYH2w3GjBmDESNGZFxPTExEWFiY7QIksrbWreUCAAcOADVqAAMGsCQqkQNh30Ru59135QJIH7Vuncy3V6qUtnER2ZBDJUbBwcGIj483abt8+TK8vLxQokQJs4/R6/XQ6/X2CI+IiNwQ+yYiIvfgUOcYRUVFYcOGDSZt69evR926deHt7a1RVERE5M7YNxERuQebJka3bt3C3r17sXfvXgBS8nTv3r048//lh8eMGYPevXtn3H/QoEE4ffo0RowYgUOHDmHevHmYO3cuRo4cacswiYjIjbBvIiIic2w6lG737t1o2rRpxnXDeOs+ffpgwYIFuHjxYkZHBADlypXD77//juHDh2PGjBkIDQ3FtGnTWA6ViIishn0TERGZY9PEqEmTJhknqJqzYMGCLG2NGzfGHkNVLiIiIitj30REROY41DlGREREREREWmBiREREREREbo+JERERERERuT0mRkRERERE5PaYGBERERERkdtjYkRERERERG6PiREREREREbk9JkZEREREROT2mBgREREREZHbY2JERERERERuj4kRERERERG5PSZGRERERETk9pgYERERERGR22NiREREREREbo+JERERERERuT0mRkRERERE5PaYGBERERERkdtjYkRERERERG6PiREREREREbk9JkZEREREROT2mBgREREREZHbY2JERERERERuj4kRERERERG5PSZGRERERETk9pgYERERERGR22NiREREREREbo+JEZEjuX5dllevahsHERERkZthYkTOad06oG9fIClJ60is69dfZblvn7ZxWFt6OjB2LDBtmtaREBEREZnFxIic09KlwMKFQIMGwIULWkdjHUoBy5bJ+unTxqNHzu7uXeC554BJk4Avv9Q6GiIiIiKzmBiRc6paVZaHDgF16wIHDmgbjzXs3g2cOiXrSgHff69pOFZx9SrQtCmwapVcr1lT03CIiIiIssPEiJyTITFKSwMuXwYiI4E//tA2pvyaOxfw9DRenzNHu1is4fhxoF49SfjS0wEvL6BaNa2jIiIiIjKLiRE5J0NiBEhydPcu0KoVsGiRdjHlx507wOLF8loM9u8H9u7VLKR82bULePxx4MwZ42tKTTX9uxERERE5ECZG5JzKlAF8fY3X09NlB7xPH+D992UomjNZuRK4fdu0zcsLmDdPm3jyY9UqoHFjICHBNNEDmBgRERGRw2JiRM7JwwN49FHzt/3vf8CAAcD9+/aNKT+++UZeU2apqXIELDlZm5jyYvp0oHNnICVFktXMPDyAihW1iYuIiIjoIZgYkfOqVcv0nJzM5s8HBg2ybzx5deoUsGVL1kQCkKMuhhLejm7ePOD11+VonbkjdmXKAHq9/eMiIiIisgATI3Je2Q3L8vQEdDqgQgX7xpNXCxZkf5uHB/D113YLJV9CQmR4o06X9TadjhXpiIiIyKExMSLnVbVq1nNYAGDIEODkSZlQ1BnkNEltejpw86bdQsmXNm2AixeBTz81thmO6LEiHRERETk4JkbkvDIfMapQAWjYUNaPHgXCw7WJKS8mTgQOH5bL9OnSVrq0sW39em3jy42iRYGwMOP1yEhZ3r/PwgtERETk0Ly0DoAozx55BFi2DPDzk1LdOp0MPVuzBoiPB4KDtY7QMno9UKmSrMfFydLPz9jmTNLTgeefl/Xr14FixYA9e2Quo+ee0zY2IiIiohwwMSLnpdMBzz5r2rZmjQzpCglxvpLdruD112XZs6ckRQDw2GNyISIiInJgHEpHrqV1a+P6n39qF4c7unULmDFD1nMqKEFERETkgJgYkes5dUqWjRrxqJE9RUTIcu5cKbZARERE5ESYGJHrCQ83nvQ/ZYq2sbiLw4eB//6T9f79tY2FiIiIKA+YGJFr2rxZliNHAikp2sbiDqpUkWVsrLZxEBEREeUREyNyTb6+wP/+J+sdO2obi6v78UdZBgQAtWtrGgoRERFRXjExItc1frws166V8t1kfenpQNeusn70qLaxEBEREeUDEyNyXTodsG6drIeEaBuLq3rtNVm+8IKxPDcRERGRE2JiRK6tZUvj+vbt2sXhipKSgJkzZZ3luYmIiMjJMTEi13f6tCyffJLlu63JUJ573jzA01PbWIiIiIjyiYkRub4yZYAGDWT9s8+0jcVVHDoEHD8u6/36aRsLERERkRUwMSL38McfsnzzTZbvtoaqVWXJ8txERETkIpgYkXvw9QXGjZP1Dh20jcXZLV0qy5IlWZ6biIiIXAYTI3IfhsRo/Xrg4kVtY3FW6elAt26yzvLcRERE5EKYGJH70OkkKQKA0FBtY3FWQ4bIsndvoGhRTUMhIiIisiabJ0YzZ85EuXLl4Ovrizp16mB7DiWTt2zZAp1Ol+Vy+PBhW4dJ7qJFC+P6tm3axeGMkpKAWbNkfd48bWMhyif2TURE9CCbJkZLly7FsGHD8PbbbyM2NhaNGjVCmzZtcObMmRwfd+TIEVy8eDHj8uijj9oyTHI3hs9f48Ys350btWrJcv58lucmp8a+iYiIzLFpYjRlyhQMGDAAL774IqpUqYKpU6ciLCwMswy/OmcjMDAQwcHBGRdP7oSRNYWFAQ0byvqnn2obi7OIiwNOnpT1vn01DYUov9g3ERGROTZLjFJSUvDPP/+gZcuWJu0tW7bEjh07cnxsREQEQkJC0Lx5c2zevDnH+yYnJyMxMdHkQvRQGzfKctQolu+2RLVqsty3T9s4iPKJfRMREWXHZonR1atXkZaWhqCgIJP2oKAgxMfHm31MSEgI5syZg+XLl2PFihWoVKkSmjdvjm05nAsyadIk+Pv7Z1zCwsKs+jrIRfn6AhMmyHr79trG4uiWLJFlUBBQs6a2sRDlE/smIiLKjpetN6DT6UyuK6WytBlUqlQJlSpVyrgeFRWFs2fP4tNPP8WTTz5p9jFjxozBiBEjMq4nJiayAyLLvPuulPDesEHKd4eEaB2R40lPB3r0kHWeaE4uhH0TERE9yGZHjAICAuDp6ZnlF7jLly9n+aUuJ5GRkTh27Fi2t+v1ehQpUsTkQmQRnU6SIoDlu7MzeLAs+/ZleW5yCeybiIgoOzZLjHx8fFCnTh1sMOx4/r8NGzagQYMGFj9PbGwsQvhLPtnKU08Z17du1S4OR5SYCHz1lax/8422sRBZCfsmIiLKjk2H0o0YMQK9evVC3bp1ERUVhTlz5uDMmTMYNGgQABlqcP78eSxatAgAMHXqVJQtWxbVqlVDSkoKFi9ejOXLl2P58uW2DJPc3dmzUqmuSRMZOpbNcBq3YzifaOFClucml8K+iYiIzLFpYtS1a1dcu3YN7733Hi5evIjq1avj999/R3h4OADg4sWLJvNGpKSkYOTIkTh//jwKFCiAatWqYfXq1Wjbtq0twyR3V7o00KgRsH078PHHwFtvaR2R9g4eBE6flvXevbWNhcjK2DcREZE5OqVca4bLxMRE+Pv7IyEhgWO6yXLJyVKpDgDu3QP0em3iWLkS6NwZqFgROHJEmxgA41GzfftYiY5yhd/B5vF9IafWujWwbh1w7hxQqpTW0RDlSm6+f206wSuR09Drgfffl/V27bSNRWvffy/LkBAmRUREROQ2mBgRGbz9tiz/+AO4cEHbWLSSng707CnrcXHaxkJERERkR0yMiAx0OmDjRll316ECr7wiy379WJ6biIiI3AoTI6LMmjc3rm/ZolkYmkhMBObMkfWvv9Y2FiIiIiI7Y2JE9KCzZ2XZtClg79ok9+/LMjXVvtsFgOrVZbloEctzExERkdthYkT0oNKlgcaNZX3yZPtue8UKWV68aN/tHjxoTAh79bLvtomIiIgcABMjInPWrZPlmDFSytte1q+X5d279j1qZDha9O+/9tsmERERkQNhYkRkjl4PfPCBrNtrEscjR4AbN4zXp0yxz3a/+06WoaFAjRr22SYRERGRg2FiRJSdsWNluWkTcP687bc3Zozp9RkzbL/NtDTghRdkneW5iYiIyI0xMXJ1Fy4ACQlaR+GcdDqZ0wiQ845sKT0dWL3atO3MGeDkSdtud9AgWQ4YAPj723ZbruruXeDUKa2jICIionzy0joAsqGUFCAiArh8GShcGAgPBx55BChXDihb1ngJD+dOcXaaNTOub94slepsYckS+Xs9aOxYuc0WEhKAb76R9a++ss02XMHdu5KknjplevnvP0lcr12TJPrYMaBCBW1jJSIiojxjYuTKPD2Bzp2B2bOBpCTgwAGpPublJUOo0tON9zUkTt99B9SsqV3MjujcOTli1KyZvGc6nfW38dFH5tt//tn62zIwFFz49luW537QrVtAy5bA0aOS+BjodOb/fxo2BIoXt3+cREREZDUcSufKPD2BmTNNS04rJXPlZN6pA4yJ059/2jdGZ1CqFNCkiaxnl8Dkx82b8t6bc/cusHKl9bd54IAkfIDxHCMyOnIEiI42TYoA8/8/3boBGzcCxYrZN0YiIiKyKiZGrk6nA0aNAn74QX7p9sjmT67TAWXKyLkmlNXatbIcO9b65bvHjcv59g8/tO72AGP1uf37rf/crqBOHeCZZ+R/JidjxwLffy9VDImIiMipMTFyF127SnW1QoXMD5tSSo5c3Llj99Ccgl4PTJwo661bW/e5V63K+fa9e627vcWLZVm6tHE4HZlKSZH/E3NzSel08gPDV19J0mqLoZVERERkd0yM3EmjRsCuXTJfTeZfwg2JUmKinCfh4QHExmoToyMbPVqWW7ZYt3z3lClA7dpyKVnS2G5oM5QNt4a0NKBXL1k/eNB6z+sqLl0CHntMEuHly7Pe7ukpt/32G/DSS/aPj4iIiGyGiZG7qVwZiImRoVSGYXXp6XJUwrAjqJTsHOp0wNy5cp3k/di0SdatWb772WclEY2NBfr1kzZvb2Pbe+9Zb1uGnfmBA4EiRaz3vM7ur7/k7xscbPxR4P33gevXgRIl5LqXl5xH9NdfQJs22sVKRERENsHEyB0FBQHbtxt37gYMAGrVkgp2SkklLsORixdflASqb18pBODuMpfrNiRJziIhAZg3T9ZnzdI2FkeQni5H63Q64IknjO0bNsj/wTvvSCJkKF5Svjywe7f8aEBEREQuh4mRuypYUEpB//knMHWq6W2PPipzH926BTz3nLQtXAj4+cm8R6dP2ztax2IYRte8uXMdTataVZbffefe5bmTkqSwgqcn8MYb0laxInD2rPw9n3rK9P59+8oQ1F27pKQ9ERERuSQmRu7M01PmXylY0PztBQsCP/4o56VMmyZtp09LcqTTAWvW2C1UhxIaapz41VCQwdHt3w9cuCDrPXpoG4tWDh8GihaVIYSG+aFeegm4d0/Kc2c3PNLTE6hXTx5LRERELouJET2chwfw2mvya3p0tLG9bVtJkMaNM1+9y5UZksJ33rF++W5bMEzam918Sa5KKUnudTqgShUZTgjIEdD0dKksx1LbREREBCZGlFuRkbKzeemS/IoOSHEAb2+gcWM5Wd0d+PgAkybJeqtW2sbyMIsWybJMGaBaNW1jsZeUFGDoUEnqu3aVNh8fYN8++fz27s0y20RERGSCiRHlTWCgnHORkiITyALAtm1SwUunA/75R9v47OGtt2S5dStw7py2sWQnLQ3o00fW3WEy1/h4KSSi1xuHf7ZqJXN0JScbj5wRERERPYCJEeWPt7dU7VIKWLnS2F63riRIX3/tXAUKckOnAzZvlvWwMG1jyc7AgbJ86SXXLs+9bZv8PUJCgH//lbaJE2WI59q1gL+/tvERERGRw2NiRNbzzDOSBB07JiXBAdkh9/CQoUuuWO67SRPj+h9/aBaGWTdvAvPny7orludOTwc+/VQSosaNje1//CGfwzFj3Lv6HhEREeUKEyOyvkcekSFNt28D3bpJ27ffSrnvsDDg5Elt47M2Q/nup55yrKNjhvLc339vnMzXFSQlAU8/LUnPm29KW+XKMpxRKWPFQCIiIqJccKG9JXI4fn7AkiXyy/6MGdJ27pxMlKnTAatXaxuftYSGGue++fBDbWMx+Pdf4OJFWe/eXdtYrCUuDihcWIYE/vqrtL3yipTbPnQIKFVK2/iIiIjIqTExItvT6YDBg+XX/F27jO3t28ttb7/t/OW+DUneu+86RvnuWrVkefCgtnHkl1LADz/I56RaNZl0GAAWL5aEe+ZMltsmIiIiq2BiRPZVr57s7F65AkRFSdvEiVLE4YkngGvXtI0vr3x8gI8+kvUWLbSNZeFCWYaHG4fTOZuUFGDIEBkCaDji5ecn8zApBfTsyXLbREREZFVMjEgbAQHAjh2yAzx6tLT99Ze063TA7t3axpcXhrLl27cDZ89qE0NaGtC3r6wbqrM5kwsXgOrV5SiQYfhl27YyMevt2+4zDxMRERHZHRMj0pa3t0yUqhTwyy/G9scflwTpq68cq6BBTnQ6mdMIkMlUtTBggCwHDXKu8txbtsj7V6qUcfjfRx/JEMvVq53rtRAREZFTYmJEjqNDB0mCjh+XggaA7OB7eMjQqTt3tI3PEk8+aVzfuNG+27550ziMznC0xZGlpwMffywJUdOmxvbNm+Vz8NZbLLdNREREdsPEiBxP+fJSAvv2bUmIACk5XbCgJEwnTmgb38NcuCDLFi3se7SrcmVZLlni2OW5ExOBdu0k6XnrLWmrVk3+5kqZzg1FREREZCcOvPdEbs/Pz1h9bPZsabt4EahQQY4yGEo2O5qQEGMBhvfft8829+0DLl2SdcPcUY7m4EH5m/r7A7//Lm2vvipV/A4cMB4lJCIiItIAEyNyfDod8PLLcjTh77+N7U8/LbeNGeN45b5/+02W48bJPDu2Vru2LOPibL+t3FAK+O47+TtVrw7cvSvthvmtvvxSKvoRERERaYyJETmXxx+Xne2rV6W8NyAn6Xt7S/nvq1e1jc/Ax0fOnwFsX757wQJZlisHVKli221ZKjlZJl/18ABeeEHaChWSo0ZKyVEtltsmIiIiB8LEiJxTiRJSFvv+fZkgFgB27gRKlpQd7sxHlrQycqQs//zTduW709KAfv1kfd8+22wjN86fl+TM19c4/LFDBym3nZTkvPMqERERkctjYkTOzcsL+OADOQphGL4GAPXrS4I0Y4Z25b51OmDbNlm3Vfnu/v1l+corQOHCttmGJTZtktdbujRw+LC0ffKJJG6//MJy20REROTwmBiR62jXTpKgEyeAsDBpGzJEhnN16yZV7uytUSPj+oYN1n3uGzeARYtk/csvrfvclkhLkzmodDqgeXNj+9at8ncYOdKxq+MRERERZcK9FnI95coBZ87IvEe9eknb0qVyjktwMPDff/aNx1C+u2VL6x69MpTn/uEH+yYgCQlAq1ZytG7sWGmrWVNep1KmczkREREROQkmRuS6ChSQIyrp6cCcOdJ26RLw6KNylOPnn+0TR0iIJBIA8N571nnOvXuBy5dlvWtX6zznw+zfD+j1QNGiwPr10vbaa1JoYd8+eZ1EREREToqJEbk+nQ4YOFCOZuzebWx/5hm5bdQoKeJgS7/8Isvx461TvjsiQpaHDuX/uXKiFPDtt/I+1awJpKRI+9KlknBOm8Zy20REROQSmBiRe6lTR3b2r10DGjeWtk8+kZ37evWAK1dss10fH9kOYHo+Tl7Mny/LChWMw+msLTlZ5o7y8AB695Y2f39JxJQCnn+e5baJiIjIpTAxIvdUvDiwZYscKfrf/6QtJgYIDJQd/p07rb/NN96Q5Y4dcg5UXqSlGSvRxcZaJ67Mzp4FKlaUctuG4YfPPAMkJgI3b9ouESMiIiLSGBMjcm9eXsCECXIU5Pffje1RUZIgTZ8uQ8asQaeTuZcAIDw8b89hmLPo1VetW55740aJr0wZ4NgxaZsyRRKxlSu1LQVOREREZAdMjIgM2rSRBOnUKaBsWWl7/XXA0xN47jnrlPt+4gnjuqGAgaVu3JDzfQA5tye/0tJkDiidDmjRwti+fbu8D8OHs9w2ERERuQ3u9RA9KDwcOHlSyn337Stty5ZJue+SJY1HVPLq4kVZtmqVu/LdFSvK8scf85ew3LwpiZCXF/Duu9IWESFxKWWavBERERG5CSZGRNkpUEAKHaSnA3PnStvVq5Kg6HTAihV5e97gYDk6BUiVOkvExsq2ATl6lRf79gHe3kCxYjJ0DpCjQsnJwJ49EhcRERGRm2JiRPQwOp0UPFBKEghDNbZnn5X1kSNzX+571SpZvveeafnumzdl+eB5TY89JsvDh3O3HaWABQskztq1gdRUaV+2TLYxZQrLbRMRUc7i42VpmD+PyEUxMSLKjYgISSiuXweaNpW2zz6T5KJuXZlA1hI+PvI4wPg8APDvv7LMnBgZjlY98ghQqZJlz3/vHjBggAy5MxRsKF5cEiuljEkdERHRw5w4Icvly7WNg8jGmBgR5UWxYsCmTXIEZsIEafvnHxmOptNJSe6HGT5cljt3AqdPSzGEffukTSlpS00FXnxR2iwpz33mjMxvVKAAMG+etD37LJCUJHM3WZpYERERAcCBA9KHAMD33+fu3FgiJ8PEiCg/PD1lHiSlgLVrje0NG0qC9MUX2Zf71umAP/+U9bJl5byfu3eNty9caCz+MGSIFH/Izvr18nzh4cZf9qZOlWTLUDiCiIgot+bPN44wOHkS2LVL23iIbIiJEZG1GKrMnT4tR20AYNgwSZ46dwZu3cr6GEMCBQAffmh621dfAd99J+tffJH1sWlpco6STifbNvjrL4lj6FCW2yYiory7f18SI8NRIk9P42gEIhfEvSYiaytTBvjvPzn6M2CAtBkmSS1eHDhyxPT+hvLdhslfDS5ckOVPP5kmODduAM2aSbntceOkzXB+k1JAgwbWf01EROR+Vq+WPscgLU1+sLtzR7uYiGzI5onRzJkzUa5cOfj6+qJOnTrY/uDO3wO2bt2KOnXqwNfXF+XLl8fs2bNtHWLeXL4MLFokw6ic9QtCKTnZ/6OPgCVLtI4m727elIlKH38cWLdO62iMfH2Bb76RoXTz50vbjRtA5cpylGfZMvkbBAUBVatm/zxdushy715JkIoXBzZvlrY33gBSUoCYGCAw0KYvJ1dOnwaiomSC3JMntY4m73btAsaMAbZtM1b0czbp6cD06TIp8H//aR2Nw3DJvkkpKbDy2WfA11/zXBDKv7lz5ShRZnfu5H26CiKDxETg7bflvLVr17SOxkjZ0A8//KC8vb3V119/reLi4tTQoUNVwYIF1enTp83e/8SJE8rPz08NHTpUxcXFqa+//lp5e3urZcuWWbzNhIQEBUAlJCRY62WItDSlYmKUGj9eqcceU0qnU0q6HaW++86627KlW7eU+vlnpQYOVCo42PgabPtRsL4dO5Tq3l2pkiVNX0OZMlpHlrO9e5Xy8jKNefhwpapVM20zXLy8lJo+PWv7ihVav5KcdepkGm+hQkq1aqXUsmXyv+QsIiKMr6FwYaWee06phQuVunxZ68gsd/y46d+iXDn5zG3YoFRystU3Z7PvYCtyqb7p7l2l1qxRasgQpcLCTP/WV65Yd1vkXi5eVMrDI2v/4+Gh1JNPah0dObtZs4yfKZ1OqXr1lPrwQ9lPSk+36qZy8/1r073hevXqqUGDBpm0Va5cWY0ePdrs/UeNGqUqV65s0vbyyy+ryMhIi7dp1c7n5k2lfvpJqb59lSpRQv54np5ZvyQcPTE6dkypL75Q6qmnjDvlD+6cO3pidPeuUtOmKVW/vlLe3uaTCGdIjAxu3JC/R3avw9ylZEmljh7VOnLLPJgYPdipVq6s1NixSl26pHWkOcucGGX+/9fp5AeSCRPkBxNHTvYeTIwy//8XKKDU008r9fXXSp07Z5XNOUNi5PR90+nTslPRtq1Sen323+lMjCg/PvnEfGJkuJw4oXWE5MxmzTI9yGDYPwCUCgxU6sUXlVq5UqmkpHxvKjffvzqlbHOsPSUlBX5+fvjpp5/QqVOnjPahQ4di79692Lp1a5bHPPnkk4iIiMAXmU40X7lyJZ5//nncuXMH3t7eWR6TnJyM5OTkjOuJiYkICwtDQkICihQpkrfgv/tOTnzfsUPG03p5PXwYTU5DobSglJxzkpRk+eSjer1tY8qLtDTnHcJElvPwAMz8f2su03dLjjw95RyykiUd73Vcv26cnNEcT0/5PwOA6tWB558H3n03z5tLTEyEv79//r6Dbcip+6YxY+R8xSNHZDiuTpd91UtAqlQWLJi3bRH9958M1c5OyZJyIcqLuLicbzfse3t5AY0aSXXezp3ztKnc9EteedqCBa5evYq0tDQEBQWZtAcFBSE+m046Pj7e7P1TU1Nx9epVhISEZHnMpEmTMMEwj4y1LFhgeiK8JTvmD/sDOwNLdwKJrC093bk/f2lpcq7bzZtaR5J7hqQIkPlKzp3LV2Lk6Jy6b/r0U2N/ZPiNNSenT1t3+0SZXbkiFyJbMHzXpabKedXFi+c5McoNmyVGBjpDKeL/p5TK0vaw+5trNxgzZgxGjBiRcd3wq1y+rFolf4TffgN++UWqhnl4ZN8RzZsHdO2av23awvnzUozg999lMtLkZNNfhx/kiF9wqalSpOD774E9ex6+8xwcDGzZYpfQ8i01FZg5Uy4P07w5MHmy88xHNGSIzMuUE09PoEoV4JlngH79AAc8uoCmTSVZeJDh/0inA+rXB9q1A1q3lteTw/ebJk6elCNB5hh+kStcWF6D4XW4Aafsm86elfnSfvtNlrdv5zyi4fRpICAgf9skAoCOHeU7/dgxIDRU62jIFXzzjUxp8uB+deaj4WFhso/Qti3QpIldwrJZYhQQEABPT88sv8Bdvnw5yy9vBsHBwWbv7+XlhRIlSph9jF6vh97aQ8AKFgTat5fLrFnAwYOSXPzyCxAdLX+szJ2RXg/4+Vk3Bmt49FG5DBkiCcXWrfI6Vq2SDvPBZM9RO9AhQ+QCyE7q1KmyU3D+fNb7+vgAlSrZNbxcu34d6NRJKp0Z1K8vw2MePOIwcKBUl/rjDynJ7e8P7Nwple0cWeHC5tuLFQMaNwYGDwZatLBvTHmReYiUIRkqVgzo0EGSiJYtgaJFNQvPIgUKGNcz/88bktK2bYHISPlOcwNO3TcFB8ukz337yhDpv/6Scso//yw7rA8Or/Pzc8y+iZyPoTJdgQL8TJF1+PgY9z8zD5tr2BB4+mnpYytWtPuPjTYr1+3j44M6depgw4YNJu0bNmxAg2zmWYmKispy//Xr16Nu3bpmx3DbhU4nv7aOGgX8+Sdw9SqwdCnQvbvsIHl5ybw1jk6vl524qVOBU6dkJ/yzzyQDd5bXAMjf4ptvZLjP7dvA559LmW7D5yM8XNv4crJnj3yeSpQwJkWjRskY7pkzTZMiwxdBUJDs5CxcKNcTEoxHJX780XHL8VapIkudThLV0aPlyOv163KOhDMkRYD8sAAAtWsD77wj5buvXpW/x/PPO35SBMiRuGLFpHx827bA7Nly5CEuDpg4EXjiCbdJigAX6pu8veX7+5NPgKNHgRMnpCx7y5aywxEYKH9zIiJHVLas/FhXsiTQuzewfLmU7d6yBRgxQvYdtBiBke9SDzkwlESdO3euiouLU8OGDVMFCxZUp06dUkopNXr0aNWrV6+M+xtKog4fPlzFxcWpuXPnOk5JVHPS0pSKj7f9dmwtMVEuzi421vGqg6WnS8WvB6v5rFxpej9D+4svyjJz5b07d4z3+/dfYxUqw+X1121Sdjnf9u1zvL9HbqWmOn7lPEtcu6bUvXt22ZQzVKVz+b7p9m2pfElkLa1aSX9jpeqVREopKQlvh/0EhynXrZRSM2bMUOHh4crHx0c99thjauvWrRm39enTRzVu3Njk/lu2bFERERHKx8dHlS1bVs2aNStX23OGTpncwJ07SvXqZZrABAVJ6fQHGRKnypWVGjXKmBh9/rmsP/541sfcvGnsqAyXmjXlS4ZIQ87yHcy+iSgXmBiRE3OIct1acfRSseTiTp4EnnxShvoZdO0qs4ebK5ubmmocBpiUBLz/PvDxx9KWnCyHmQ3PW7Zs1senpcn9x441bd+2TcpbEtkZv4PN4/tCTq11aynmdO4cUKqU1tEQ5Upuvn9tdo4RkVtZvVrGwpYvb0yKZs6U84N++CH7uUR695bl0KFZK87pdHJyNQCUK2f+8Z6eMreJUlKgweDJJ+Xxn36a8zwnRERERASAiRFR3qWmygn5Op1UMDTYtUsSlVdeyfnEwevXgSVLZH3KFPP3adDAWA1ozZqc42nWTLZ77pyxat2bb8rjO3SQI1JEREREZBYTI6LcunZNhql5ewMffihtUVEyD5RSQL16lj1PhQqyXL7cOGTOnAsXZNm2rWVV6EqVAg4dAu7dk+QMkHlPihSRMtquMBkxERERkZUxMSKy1O7dcgQoIEBKtwNShjolBdixI3fzQP3zj7E898Nmcg4MlCM+gByhspRebxzOt3ixtN26BVSrJq/jhx8ct9w3ERERkZ0xMSLKiVLAV19JIvH448b2n3+W2yZNMp0E1FJ168ry6FHL7r98uSwnTgTu3s3dtnQ6oGdPiXf/fuOEn927y5Gq116T5I6IiIjIjTExIjLnzh1JJjw8gEGDpK1UKeD4cUkwnn467889Z44sq1QxTiD6MN7eMjkvIIUV8qp6dXltCQkyNA8AvvxSji5Vr24ctkdERETkZpgYEWV24oQkQAULAt9/L209ewK3b0tRg/Ll8/f8qanAyy/LekxM7h77+uuy3L0bOHUqf3EUKSKV9FJTgcmTpe3gQXntOp3MPE1ERETkRpgYEQHAr79KQlChgvGoyezZxvNz/Pyss51evWQ5fHj2Jbyzo9PJuUxA9uW7c8vTExg1So6Cbd5sbG/aVLb38ccs901ERERugYkRua/UVJkYVaczHRr399+SKLz8cs7ltnPr2jUpeADI/EJ5ERVlPKdp9WrrxGXQpIm87vPnpUADALz1liRP7doBiYnW3R4RERGRA2FiRO7n6lWZH8jbW4onAMATTxjLbWcusmBNhmF4K1bkXJ77YQwTyLZvb5uqcqGhwIEDUu771Vel7fffAX9/Ocp18KD1t0lERESkMSZG5D7+/luOAJUsCURHS9s77wD37wPbt+eu3HZu7d5tPOLSqVP+nisw0HiEa+zY/D1XTvR6KcyQnm6ciPbOHSnSoNPJOVgs901EREQugokRuTalZC4fnQ6oX9/Y/ttvctv77wNeXraPw3AU6tgx6zzfsmWy/Oij3Jfvzi2dDujWTd6vgwdlkljAWLVv8GAgOdm2MRARERHZGBMjck137hjn6TEMBwsLk6pzSsk5M/by1VeyrFYNeOQR6zyntzfwxRey/sQT1nlOS1StKke+EhKMk87OmgX4+kr58fPn7RcLERERkRUxMSLXcvw4EBIi58IYCh306iWJ0pkz1qvmZqnUVOM8SLt2Wfe5DeW79+wBTp607nM/TJEiwC+/AGlpwCefSNvhw0Dp0nKEadMm+8ZDRERElE9MjMg1/Pyz7JA/8ggQHy9tX30l58csWgQUKKBNXD17ynLEiNyX57aE4Vyp/M6vlFceHsDIkXIULvPcR82by9/jo48keSIiIiJycEyMyHndvy/lpHU64JlnjO27d8uO+ksvWbfcdm5duwb8+KOs57U898NERgI+PrL+22+22YalGjeW9/3CBaBmTWkbM0bO4WrdWobfERERETkoJkbkfK5ckUIKPj4yASkgO+XXrsmOeZ062sZnYBi2t2qVbRM0Q/nuDh0co0pcSAiwb58UZDAM91u3DihaVCrd7d+vaXhERERE5jAxIuexc6ckGIGBUnobAP73PzlytGULULy4puGZ2L0bSEqS9Y4dbbutkiWNR8zGjLHttnLDx0cKRKSnA0uXSltKihxN0umAb791jESOiIiICEyMyNEpJXPp6HRAVJSx/fff5bYJE+xTbju3rF2e+2EMQ/YmT5ZCE45EpwOef17+XocOyUSxANC7t5yj9PLLLPdNREREmmNiRI7p9m3ZmfbwAF57TdrKlpXqa0oBbdpoGl6OZs2SZfXq1ivP/TDe3sD06bLesKF9tpkXlSsDN29KyW/DkbQ5c6Tcd6VKxmGBRERERHbGxIgcy7FjMjSsUCHgp5+krV8/OQpy8qQkR44sNVUmPAWsX577YYYMkeXevTJfkyMrXFjOvUpLA6ZMkbajR2WuKZ0O+OMPTcMjIiIi98PEiBzDypWyQ1yxInD1qrTNnSvnp8ybp1257dzq0UOWI0cCfn723/7OnbKsUMH+284LDw9g+HA5Crh9u7H9qafk8/DBByz3TURERHbBxIi0c/8+8OabsgPcubO06XQyYalSQP/+2pbbzq2rV41HuQzV8uytfn0ZlgYAv/6qTQx59cQT8ne/eBGIiJC2d9+Vc8hatJAheEREREQ2wsSI7O/yZSlO4ONjnN+naVPg+nU5QmTYKXY2hmF+hslmtXLmjCyffto5q74FB0tynJwMDBsmbRs3AsWKyblU//6raXhERETkmpgYkf1ER0vCEBQk5awBYPx4OXK0aZPs+Dqrv/+WghGAJCRaKlnSeATurbe0jSU/fHyAzz+XZNlwJC41FahVSz5HCxY4Z+JHREREDomJEdlWejowbZrsyDZoYGxfu1Z2aseNc8xy27lVv74s//tP2zgMfvhBlp984njlu3NLpwO6dJHPy+HDxvmq+vWTc5RefBG4d0/bGImIiMjpMTEi27h1C3j2WcDTExg6VNoqVABOn5Yd3FattI3PmmbOlGXNmo5T9MDbW+Z/AkwTUmdXqRJw7ZpMnvvss9I2d64U53jkEeDsWW3jIyIiIqfFxIis6+hRoEQJKce8YoW0DRgA3L0rR1PKlNE2PmtLTQVefVXWo6O1jeVBhrj27QOOH9c2FmsrVAhYtkwq1k2dKm3Hj8vnS6cDNmzQNDwiIiJyPkyMKP+UApYvlx3SSpWkiAIAzJ8vQ+m++cZYKc3VdOsmyzff1KY898MYynfba6JZe/PwkCOSSgF//WVsb9lSPo/vvcdy30RERGQRJkaUd/fvAyNGyM5ply7S5ukJxMbKjmrfvs5Vbju3rl6VhBAAJk/WNpbs1K9vTNh++UXbWGytQQP53MXHA3XqSJvhHLZmzYAbN7SNj4iIiBwaEyPKvUuXgMceM1YNA2RCzhs3ZGhZ7dqahmc34eGy/OUXx04AT5+WZceO7lHFzVD1MCUFeOMNadu8WYo2eHgAe/dqGh4RERE5JiZGZLm//pIEIDhYjgoBwPvvSzK0YQNQtKim4dnVrl3Gam8dOmgby8MEBBgLFYwapW0s9uTtLfNkKWU8300pmSdLpwPmzXOPRJGIiIgswsSIcpaeLkeFdDrgiSeM7evXy07lO+/I8Dl3ExkpS2cpamAo3/3pp85fvjsvOnWSz+vRozLPEyBFQTw8pOw3y30TERG5PSZG7uyff2Qy0u3bs95265bsTHp6ynlEAPDoo8CZM7KD2aKFfWN1JDNmyLJ2baB8eU1DsZiXl7GsuCGpc0ePPgpcviyf7+eek7YFC6Tcd7lyxmGHmZ07J0fcDMklERERuSQmRu5q9Wo5AvTrr/LL+f370n74MFCsmJTbXrVK2gYOlHLbR48CYWGahewQ7t8HhgyR9cxV0JzBK6/Icv9+x5mIVisFCwI//igV66ZNk7ZTp4CyZeXo6Nq1xvuOGiVD8bp3lyIbHH5HRETkkpgYuaOvvpLzYpKT5fqxY0DPnrJDWKUKcPOmtC9cKEPp5sxx3XLbudW1qyxHjbJeee70dOMl80535nZr+ftvWT76qPWe05l5eACvvSbve+Z5qNq0kf+HXr2AJUuM7aNHS4KZmmr/WImIiMimmBi5k/R02bEbNEh2BDPvhP/0kyz1epkQVCmgd2/HrrZmb1euACtXyvpHH1nvedu0kSGLnp7AJ59I2/37xrbgYOtt6/HH5WgJYDwiSCIyUj73ly9LmXMAWLw46/3mzJEhqLdu2Tc+IiIisikmRu4iOdk4FCg7gwbJSeg1a9ovLmdSpowsf/3Vugnj1as53373rvW2BciQMcBYkIBMlSwpE+MuWmT+dqWk+MgTT8icSUREROQSmBi5gxs3gObNgWXLcr7fnDnAoUP2icnZ7NplrFzWvr11n/udd3K+3TB5rrUEBBif8803rfvcruL2bXlvskuA09KAAweAunX5P0NEROQimBi5ulOnZFjQzp3Zn6vi5SVDttLTgU2b7Bqe07Blee5OnaQqWnY+/ND62zScN/PZZ5IEkKm4OJnIWKeT+ZDMSUuTI0b16wNbt9o3PiIiIrI6L60DIBu6d0/OKbl6Nesv315eQGgoUKGCXMqWlUunTlpE6ti+/FKWERG2K8/dsaP5ctAVKsjfydq8vIBZs6SQQGSkVKojo8cfl3Owjh2THxdOnJBKfmfPms55lJ4OJCUBTz0F7NkD1KihVcRERESUT0yMXJmvL/DCC3KSuCHxMVyCg91zYtbcun9fqpYBti3P/dFH5hOjYcNst81BgyQxOnBAEgBWqjPVsWPWNqWAa9ckWcp8iY8HKla0b3xERERkVUyMXN3nn2sdgXMzTAI6enTOw93yKzxcijucOWNs8/QEBg+23TYBICZGjo5UrMhCDJbQ6eQcrYAAOb+IiIiIXAbPMSLKzpUrwM8/y/rEibbf3quvml5v1kzm2bGlunWBQoVk3VCKnIiIiMgNMTEiyk7p0rL87Tf7zOc0YoTpdmxRdMEcQ/nuzp2tO5ksERERkRNhYkRkzs6dQEqKrLdrZ59tenkBVavKuoeHDHGzhxIlgOefl/WRI+2zTSIiIiIHw8SIyJyoKFmeOGHf7XbvLstixey73e++k+Xnn7N8NxEREbklJkZED5o2TZZ16gDlytl324YjRiVK2He7Xl7AV1/Jer169t02ERERkQNgYkSU2f37wNChsr59u7ax2NtLL8kyLk7KdxMRERG5ESZGRJl16SLLMWNsW57bUcXEyJJz8hAREZGbYWJEZHD5MvDLL7Jur4pwjqZuXaBIEVlfsULbWIiIiIjsiIkRkYGhPPfq1fYpz+2oDAUnnn2W5buJiIjIbTAxIgKA6Gg5vwgA2rbVNhatlSgBdOsm62+8oW0sRERERHbCxIgIABo0kOXJk9rG4Si+/VaWU6cCt25pGgoRERGRPdgsMbpx4wZ69eoFf39/+Pv7o1evXrh582aOj+nbty90Op3JJTIy0lYhEokvvpBl3bpA2bKahuIwMpfvttdEs0R2wL6JiIiyY7PEqEePHti7dy/Wrl2LtWvXYu/evejVq9dDH9e6dWtcvHgx4/L777/bKkQiGT43bJisb9umaSgOx1C++/Bh4OhRbWMhshL2TURElB0vWzzpoUOHsHbtWuzcuRP169cHAHz99deIiorCkSNHUKlSpWwfq9frERwcbIuwiLLq3FmWb7/tnuW5H2b3bjmSVqkSoJTW0RDlC/smIiLKiU2OGEVHR8Pf3z+j4wGAyMhI+Pv7Y8eOHTk+dsuWLQgMDETFihUxcOBAXL58Ocf7JycnIzEx0eRCZJFLl4DffpP199/XNhZHVacOULSorC9frmkoRPnFvomIiHJik8QoPj4egYGBWdoDAwMRHx+f7ePatGmD7777Dps2bcJnn32GmJgYNGvWDMnJydk+ZtKkSRljxf39/REWFmaV10BuoFQpWf7+u3uX536Y48dl2aULy3eTU2PfREREOclVYjR+/PgsJ6A+eNm9ezcAQGdmR1MpZbbdoGvXrmjXrh2qV6+ODh06YM2aNTh69ChWr16d7WPGjBmDhISEjMvZs2dz85LIXf31F5CWJutt2mgbi6MrXhzo3l3Whw/XNhYiM9g3ERGRNeTqHKMhQ4agm2F+k2yULVsW//77Ly5dupTltitXriAoKMji7YWEhCA8PBzHjh3L9j56vR56vd7i5ySCUsATT8j6qVOahuI0Fi0CliwBpk0DPvwQKFRI64iIMrBvIiIia8hVYhQQEICAgICH3i8qKgoJCQn4+++/Ua9ePQDArl27kJCQgAaG+WIscO3aNZw9exYhISG5CZMoZ1OnyrJePSA8XNNQnIaXF/D118DAgVKM4fBhrSMiysC+iYiIrMEm5xhVqVIFrVu3xsCBA7Fz507s3LkTAwcORPv27U2q/lSuXBkrV64EANy6dQsjR45EdHQ0Tp06hS1btqBDhw4ICAhAp06dbBEmuaP794ERI2R9yxZNQ3E6L74oyyNH5ELkZNg3ERFRTmw2j9F3332HGjVqoGXLlmjZsiVq1qyJb7/91uQ+R44cQUJCAgDA09MT+/fvR8eOHVGxYkX06dMHFStWRHR0NAoXLmyrMMndPPOMLN95h+W582LPHllWrqxtHER5xL6JiIiyo1PKtSYnSUxMhL+/PxISElCkSBGtwyFHcukSYJiHJD3dMSvRrVwpcytVrOi4R2VKlACuXwd++kkq1RFlwu9g8/i+kFNr3RpYtw44d85Y0ZXISeTm+9dmR4yIHE5oqCzXrHHMpMhZ/PefLJ97juW7iYiIyGUwMSL38Oefxp341q21jcXZFSsG9Ogh68OGaRoKERERkbUwMSLXpxTQqJGsszy3dSxcKMvp04Fbt7SNhYiIiMgKmBiR6/v8c1lGRrI8t7V4eQHffCPrjz2mbSxEREREVsDEiFxbSgrwxhuyvnmztrG4mgEDZHnsGOc1IiIiIqfHxIhcm6E89//+B/j6ahqKS4qNlWWVKtrGQURERJRPTIzIdcXHSwU6ABg/XtNQXFbt2lK+GwB+/FHTUIiIiIjyg4kRua6QEFmuXcvy3LZ07Jgsu3Zl+W4iIiJyWkyMyDVt325cb9VKuzjcQbFiwAsvyPrrr2sbCxEREVEeMTEi16MU8OSTsn76tLaxuIsFC2Q5YwaQlKRpKERERER5wcSIXM+UKbKMigLKlNE2Fnfh6QnMmyfrLN9NREREToiJEbmWlBRg5EhZ37RJ21jcTb9+svzvP+DQIW1jISIiIsolJkbkWp5+WpbjxrE8txYM5burVtU2DiIiIqJcYmJEzuvWLaBuXaBpU+C334Dz54F16+S2ceO0jc1d1a4NlCwp60uXArt3Az17AuXLAwcOaBoaERERUU68tA6AKM8OHAD++Qfw8AC2bDG2L1/uXOW5Y2OBr76S9ePHZXnsGDBokKy3bAl07qxNbHkRFyfJUbduct3LC0hLA/74A6heXdvYiIiIiLLBxIicV1ycLB+cO6d3b+DPP4GhQ4HwcPvHlVsffwz88IMkEEoZ2+fOlYRizRrnSIwSEoCvvwY+/9y0PTUV8PY2/r2IiIiIHBCH0pHziouTHe4H3b4tO+cjRtg/prx49llZpqZKIgRIgpSaKkfDunbVLrbcmDsXePNN4MKFrLfdvw/s32//mIiIiIgsxMSInNfBg7LD/SBPT0CvB/r2tXtIedKhA+Dvb/62tDRjtTdH17o1EBIi7785Bw+aHhEjIiIiciBMjMh5mTsC4eUlSca2bZJwOAO9Xob/eT0wstXDQ4pLVKmiTVy5VbWqFFuoWtV8cpSYCFy5Yv+4iIiIiCzAxIic0+3bUoUuM09POacoJgaoV0+buPJqwAAZOpeZUsBLL2kTT16FhgJ//QU0a2a+AAbPMyIiIiIHxcSInNPhw6bXPTwkGfr7bykN7Wxq1QJq1jRt8/FxnvOLMitcGPj9d6B/f9N2nY6JERERETksJkbknB7cwX72WWDTJqB4cW3isYaBA43rOh3w3HNAkSLaxZMfXl5Soe7DD41tnp5MjIiIiMhhMTEi53TokHH9zTel3LWvr3bxWEOPHsZzc5TKesTF2eh0wNixwOLFkiilpnKSVyIiInJYTIzIOcXHy473jBkyD5CHC3yUixcHnnpK1gsWBBo31jYea+nZE1i/HihUyHwpbyIiIiIH4AJ7k+SW5swB9u0DBg/WOhLr6tRJlpUquUayZ9C0KbB3r1StIyIiInJAXg+/C5ED8vICatTQOgrri4qSpSu+tgoVtI6AiIiIKFsu9JM0kQswHCV6cE4jIiIiIrIpJkZEREREROT2mBgREREREZHbY2JERERERERuj4kRERERERG5PSZGRERERETk9pgYERERERGR22NiREREREREbo+JERERERERuT0mRkRERERE5PaYGBERERERkdtjYkRERERERG6PiREREREREbk9JkZEREREROT2mBgREREREZHbY2JERERERERuj4kRERERERG5PSZGRERERETk9pgYERERERGR22NiREREREREbo+JERERERERuT0mRkRERERE5PaYGBERERERkdtjYkRERERERG6PiREREREREbk9JkZEREREROT2mBgREREREZHbY2JERERERERuz2aJ0YcffogGDRrAz88PRYsWtegxSimMHz8eoaGhKFCgAJo0aYKDBw/aKkQiInIz7JuIiCg7NkuMUlJS8Nxzz+GVV16x+DEff/wxpkyZgi+//BIxMTEIDg5GixYtkJSUZKswiYjIjbBvIiKi7NgsMZowYQKGDx+OGjVqWHR/pRSmTp2Kt99+G507d0b16tWxcOFC3LlzB99//72twiQiIjfCvomIiLLjMOcYnTx5EvHx8WjZsmVGm16vR+PGjbFjxw4NIyOysbt3gX375HL4sLSdOWNsu3ZN2/iI3Bj7JnJbly8b+6Hz56Xt4EFjW0qKtvER2YCX1gEYxMfHAwCCgoJM2oOCgnD69OlsH5ecnIzk5OSM64mJibYJkMhWRo8Gpk0zbduwAahdW9Yfewz45x+7h0VE7JvITSkFPPEEcOyYaXurVsb1//0PmDDBvnER2ViujhiNHz8eOp0ux8vu3bvzFZBOpzO5rpTK0pbZpEmT4O/vn3EJCwvL1/aJ7K5kyexv8/AAQkLsFwuRE2LfRGRlOh0QGCjL7OTUdxE5qVwdMRoyZAi6deuW433Kli2bp0CCg4MByK9zIZl2BC9fvpzll7rMxowZgxEjRmRcT0xMZAdEzqVPH/nlTamst6WnAwMG2D8mIifCvonIBl58EfjrL/O3eXkBPXrYNx4iO8hVYhQQEICAgACbBFKuXDkEBwdjw4YNiIiIACDVg7Zu3YrJkydn+zi9Xg+9Xm+TmIjsIiwMaNYM2LIFSEszva1YMaBdO03CInIW7JuIbKBLF2DwYDkPNjMvL6BTJ6B4cW3iIrIhmxVfOHPmDPbu3YszZ84gLS0Ne/fuxd69e3Hr1q2M+1SuXBkrV64EIMMUhg0bhokTJ2LlypU4cOAA+vbtCz8/P/TgrxLk6gYOzJoUeXkBffsCPj6ahETkitg3EVmoUCGgWzfpizJLTQX699cmJiIbs1nxhf/9739YuHBhxnXDL22bN29GkyZNAABHjhxBQkJCxn1GjRqFu3fvYvDgwbhx4wbq16+P9evXo3DhwrYKk8gxdOwIFCkCZD5Bm50PkdWxbyLKhQEDgPnzTduCgoAWLbSJh8jGdEqZO7HBeSUmJsLf3x8JCQkoUqSI1uEQWe7VV4E5cyQh0umkKt2ePVpHRZQr/A42j+8LOSWlgEceAU6ckOuenlJJ9YMPtI2LKBdy8/3rMPMYEbm9AQMkKTIYOFC7WIiIiHQ66Ys8/n93MS1NhngTuSgmRkSOIiICqFZN1r29ge7dtY2HiIiod29j1dSGDeUIEpGLYmJE5CgMv8wBQOfOQNGimoZDRESE0FDjxK4vvqhtLEQ2ZrPiC0SUBy+8AJw/z/khiIjIcYwbB9SoISW8iVwYEyMiR1KiBPDxx1pHQUREZBQZKRciF8ehdERERERE5PaYGBERERERkdtjYkRERERERG6PiREREREREbk9JkZEREREROT2mBgREREREZHbY2JERERERERuj4kRERERERG5PSZGRERERETk9pgYERERERGR22NiREREREREbs9L6wCsTSkFAEhMTNQ4EiIi92P47jV8F5Ng30REpI3c9EsulxglJSUBAMLCwjSOhIjIfSUlJcHf31/rMBwG+yYiIm1Z0i/plIv9rJeeno4LFy6gcOHC0Ol0Nt9eYmIiwsLCcPbsWRQpUsTm23NWfJ8sx/fKMnyfLGfP90ophaSkJISGhsLDg6O1DezZN/F/w3J8ryzD98lyfK8s46j9kssdMfLw8EDp0qXtvt0iRYrwH8ACfJ8sx/fKMnyfLGev94pHirLSom/i/4bl+F5Zhu+T5fheWcbR+iX+nEdERERERG6PiREREREREbk9Jkb5pNfrMW7cOOj1eq1DcWh8nyzH98oyfJ8sx/fKvfDvbTm+V5bh+2Q5vleWcdT3yeWKLxAREREREeUWjxgREREREZHbY2JERERERERuj4kRERERERG5PSZGRERERETk9pgY5dKHH36IBg0awM/PD0WLFrXoMUopjB8/HqGhoShQoACaNGmCgwcP2jZQB3Djxg306tUL/v7+8Pf3R69evXDz5s0cH9O3b1/odDqTS2RkpH0CtpOZM2eiXLly8PX1RZ06dbB9+/Yc779161bUqVMHvr6+KF++PGbPnm2nSLWXm/dqy5YtWT47Op0Ohw8ftmPE9rdt2zZ06NABoaGh0Ol0WLVq1UMf486fKVfFvsky7Jeyx77JcuybHs5Z+yYmRrmUkpKC5557Dq+88orFj/n4448xZcoUfPnll4iJiUFwcDBatGiBpKQkG0aqvR49emDv3r1Yu3Yt1q5di71796JXr14PfVzr1q1x8eLFjMvvv/9uh2jtY+nSpRg2bBjefvttxMbGolGjRmjTpg3OnDlj9v4nT55E27Zt0ahRI8TGxmLs2LF4/fXXsXz5cjtHbn+5fa8Mjhw5YvL5efTRR+0UsTZu376NWrVq4csvv7To/u78mXJl7Jssw37JPPZNlmPfZBmn7ZsU5cn8+fOVv7//Q++Xnp6ugoOD1UcffZTRdu/ePeXv769mz55twwi1FRcXpwConTt3ZrRFR0crAOrw4cPZPq5Pnz6qY8eOdohQG/Xq1VODBg0yaatcubIaPXq02fuPGjVKVa5c2aTt5ZdfVpGRkTaL0VHk9r3avHmzAqBu3Lhhh+gcEwC1cuXKHO/jzp8pd8C+KXvsl7LHvsly7Jtyz5n6Jh4xsrGTJ08iPj4eLVu2zGjT6/Vo3LgxduzYoWFkthUdHQ1/f3/Ur18/oy0yMhL+/v4Pfd1btmxBYGAgKlasiIEDB+Ly5cu2DtcuUlJS8M8//5h8FgCgZcuW2b4n0dHRWe7fqlUr7N69G/fv37dZrFrLy3tlEBERgZCQEDRv3hybN2+2ZZhOyV0/U2TKHfsm9kvmsW+yHPsm23GUzxQTIxuLj48HAAQFBZm0BwUFZdzmiuLj4xEYGJilPTAwMMfX3aZNG3z33XfYtGkTPvvsM8TExKBZs2ZITk62Zbh2cfXqVaSlpeXqsxAfH2/2/qmpqbh69arNYtVaXt6rkJAQzJkzB8uXL8eKFStQqVIlNG/eHNu2bbNHyE7DXT9TZMod+yb2S+axb7Ic+ybbcZTPlJfdtuTAxo8fjwkTJuR4n5iYGNStWzfP29DpdCbXlVJZ2pyBpe8VkPU1Aw9/3V27ds1Yr169OurWrYvw8HCsXr0anTt3zmPUjiW3nwVz9zfX7opy815VqlQJlSpVyrgeFRWFs2fP4tNPP8WTTz5p0zidjTt/ppwJ+ybLsF+yDvZNlmPfZBuO8JliYgRgyJAh6NatW473KVu2bJ6eOzg4GIBkwiEhIRntly9fzpIZOwNL36t///0Xly5dynLblStXcvW6Q0JCEB4ejmPHjuU6VkcTEBAAT0/PLL8q5fRZCA4ONnt/Ly8vlChRwmaxai0v75U5kZGRWLx4sbXDc2ru+plyRuybLMN+KX/YN1mOfZPtOMpniokR5IMeEBBgk+cuV64cgoODsWHDBkRERACQMapbt27F5MmTbbJNW7L0vYqKikJCQgL+/vtv1KtXDwCwa9cuJCQkoEGDBhZv79q1azh79qxJx+2sfHx8UKdOHWzYsAGdOnXKaN+wYQM6duxo9jFRUVH49ddfTdrWr1+PunXrwtvb26bxaikv75U5sbGxLvHZsSZ3/Uw5I/ZNlmG/lD/smyzHvsl2HOYzZddSDy7g9OnTKjY2Vk2YMEEVKlRIxcbGqtjYWJWUlJRxn0qVKqkVK1ZkXP/oo4+Uv7+/WrFihdq/f7/q3r27CgkJUYmJiVq8BLtp3bq1qlmzpoqOjlbR0dGqRo0aqn379ib3yfxeJSUlqTfeeEPt2LFDnTx5Um3evFlFRUWpUqVKucx79cMPPyhvb281d+5cFRcXp4YNG6YKFiyoTp06pZRSavTo0apXr14Z9z9x4oTy8/NTw4cPV3FxcWru3LnK29tbLVu2TKuXYDe5fa8+//xztXLlSnX06FF14MABNXr0aAVALV++XKuXYBdJSUkZ30MA1JQpU1RsbKw6ffq0UoqfKXfBvsky7JfMY99kOfZNlnHWvomJUS716dNHAchy2bx5c8Z9AKj58+dnXE9PT1fjxo1TwcHBSq/XqyeffFLt37/f/sHb2bVr11TPnj1V4cKFVeHChVXPnj2zlKvM/F7duXNHtWzZUpUsWVJ5e3urMmXKqD59+qgzZ87YP3gbmjFjhgoPD1c+Pj7qscceU1u3bs24rU+fPqpx48Ym99+yZYuKiIhQPj4+qmzZsmrWrFl2jlg7uXmvJk+erCpUqKB8fX1VsWLF1BNPPKFWr16tQdT2ZSgF++ClT58+Sil+ptwF+ybLsF/KHvsmy7Fvejhn7Zt0Sv3/mU1ERERERERuiuW6iYiIiIjI7TExIiIiIiIit8fEiIiIiIiI3B4TIyIiIiIicntMjIiIiIiIyO0xMSIiIiIiIrfHxIiIiIiIiNweEyMiIiIiInJ7TIyIiIiIiMjtMTEiIiIiIiK3x8SIiIiIiIjcHhMjIiIiIiJye/8HcUWwTOASM2QAAAAASUVORK5CYII=",
      "text/plain": [
       "Figure(PyObject <Figure size 1000x400 with 2 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "D, Q = 2, 5\n",
    "# K(i,j,k)\n",
    "K_ABA = [(i ≠ j) ? 0.0 : -1.0 for i=1:Q, j=1:Q, k=1:Q]\n",
    "K_ABB = [(i ≠ k) ? 0.0 : -1.0 for i=1:Q, j=1:Q, k=1:Q]\n",
    "# e(i)\n",
    "e = [ # lattice velocities\n",
    " 0.  0.\n",
    " 1.  0.\n",
    " 0.  1.\n",
    "-1.  0.\n",
    " 0. -1.\n",
    "]\n",
    "ω = [ # weights\n",
    " 1/3\n",
    " 1/6\n",
    " 1/6\n",
    " 1/6\n",
    " 1/6\n",
    "]\n",
    "# g(j,k)\n",
    "g =[ # relative velocites\n",
    "    0. 1. 1. 1. 1.\n",
    "    1. 0. √2 2. √2\n",
    "    1. √2 0. √2 2.\n",
    "    1. 2. √2 0. √2\n",
    "    1. √2 2. √2 0.]\n",
    "\n",
    "mA, mB = 1.0, 2.0\n",
    "mC = mA + mB\n",
    "e′ = [mA/mC * e[j,d] .+ mB/mC * e[k,d] for j=1:Q, k=1:Q, d=1:2]\n",
    "using PyPlot\n",
    "figure(figsize=(10,4))\n",
    "subplot(1,2,1)\n",
    "title(\"Post-Collision velocities\")\n",
    "for j=1:Q,k=1:Q\n",
    "    arrow(0., 0., e′[j,k,1], e′[j,k,2], head_width=0.05, color=\"r\")\n",
    "end\n",
    "subplot(1,2,2)\n",
    "title(\"Lattice velocities\")\n",
    "for i=1:Q\n",
    "    arrow(0., 0., e[i,1], e[i,2], head_width=0.05, color=\"r\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "113887bf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×5×5 Array{Float64, 3}:\n",
       "[:, :, 1] =\n",
       " 0.0  1.0  1.0  1.0  1.0\n",
       "\n",
       "[:, :, 2] =\n",
       " 1.0  0.0  1.0  1.0  1.0\n",
       "\n",
       "[:, :, 3] =\n",
       " 1.0  1.0  0.0  1.0  1.0\n",
       "\n",
       "[:, :, 4] =\n",
       " 1.0  1.0  1.0  0.0  1.0\n",
       "\n",
       "[:, :, 5] =\n",
       " 1.0  1.0  1.0  1.0  0.0"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K_ABC = zeros(Q, Q, Q)\n",
    "for j=1:Q, k=1:Q\n",
    "    i0, i′, i″ = 1, j, k\n",
    "    if j == k\n",
    "        continue\n",
    "    elseif j == 1\n",
    "        i′ = (k == Q ? k - 1 : k + 1)\n",
    "    elseif k == 1\n",
    "        i″ = (j == Q ? j - 1 : j + 1)\n",
    "    elseif (j + k) % 2 == 0\n",
    "        i″ = (j + k) ÷ 2\n",
    "    elseif (j + k) % 2 == 1\n",
    "        nothing\n",
    "    end\n",
    "    \n",
    "    A = [1.00000 1.00000 1.00000\n",
    "         e[i0,1] e[i′,1] e[i″,1]\n",
    "         e[i0,2] e[i′,2] e[i″,2]]\n",
    "    b = [1.000; e′[j,k,1]; e′[j,k,2]]\n",
    "    \n",
    "    K_ABC[[i0, i′, i″], j, k] .= A \\ b\n",
    "end\n",
    "sum(K_ABC, dims=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c76f25f5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.3",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}