path: root/lbm_codegen.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Characteristic constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = [Matrix(x) for x in [(-1, 1), ( 0, 1), ( 1, 1), (-1, 0), ( 0, 0), ( 1, 0), (-1,-1), ( 0, -1), ( 1, -1)]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left [ \\left[\\begin{matrix}-1\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}-1\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}-1\\\\-1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\-1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\-1\\end{matrix}\\right]\\right ]$$"
      ],
      "text/plain": [
       "⎡⎡-1⎤  ⎡0⎤  ⎡1⎤  ⎡-1⎤  ⎡0⎤  ⎡1⎤  ⎡-1⎤  ⎡0 ⎤  ⎡1 ⎤⎤\n",
       "⎢⎢  ⎥, ⎢ ⎥, ⎢ ⎥, ⎢  ⎥, ⎢ ⎥, ⎢ ⎥, ⎢  ⎥, ⎢  ⎥, ⎢  ⎥⎥\n",
       "⎣⎣1 ⎦  ⎣1⎦  ⎣1⎦  ⎣0 ⎦  ⎣0⎦  ⎣0⎦  ⎣-1⎦  ⎣-1⎦  ⎣-1⎦⎦"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = [Rational(*x) for x in [(1,36), (1,9), (1,36), (1,9), (4,9), (1,9), (1,36), (1,9), (1,36)]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAAVCAYAAABi32uwAAAABHNCSVQICAgIfAhkiAAABzlJREFUeJztnWuoFVUUx3++wqKHWZEEkhZC4Y1I0jJMt0VShqFSXwIhpPoSFIhF9KBjlFlRaFpEVJr5IeghfQgrE063J0pZUWkvtRdpXtOyNM1HH9Ye7pxhzjkzc/acM/vM+sFh7t2zZ85ea/9n1r77dUFRFEVRFEVxggGOhj6bO1oaRVEURVGU4nAqte2ko8GJgTGZ3wUWAMsi6dcCS4H3gL/sTVYlLMDpwGHgiVDaw8A64GdgP/AHsBG4Dzilwb0uBV4FfgMO2OPbwPSEZWkHLn01AJgLfAzsBfYhfroVGOSuyE5waXdWfUDxNeL6WfJJI2Hm0P9CurFJXp/tLmt9a8xIjsaMdHa71Ack18g+pH20APix3s2MNaRS5/xn9vxeYBPpKvtmm39qKO0gUtnPA4sQh26w+X4FRsbc5x57fiewHFgIPGOveyRhWdqBS1+ttL/vAJ4FlgBf2bRXkAenKLi0O4s+wA+NuH6WfNJIwEhgD+KDJI0pn+0ua31rzEiOxox0drvSB2TXSJVQz1QYQ+PG1FRgDFIRQd6klb0G6KO2VTy0Tt4H7b2fiqRfZ9PXAifEXDckYVnagStfzbTXbkG6FwOGAKvtuRtcFNgRLjWSVh/gj0Zc+sk3jYDY/Q7wA/AoyRpTPttd1vrWmJEcjRnp7HahD2hNI1UyNqbi8iYx+iSk22x5grwA59NvXMBARBz/AKclvE9RMGT3VfAXxi0xeXvsuU9aL2IuGPLRSJw+wF+NGFrzk48auQ04AkxG3jfNGlPdYjeUs75BY0YaDBozGuFCH9C6RqqEGlODM9wgLVcDxwCvJcw/wx6/CKVdAoxGuih323v2AP8C64GPnJS080R9NcIet8TkDdLGAcOQIRNfSaOROH1AOTQS5yffNHIu0gW/BOgFLktwTTfYnQW1OxkaM8oVM1zoAxxrpB2NqVlIyy/aKgyYDxyPtDYvBCYhRi8K5RlvjzuAT4HzIvfoRSay7XRT5I4R9VWfPY6OyXtW6OdzkLFiX2mkkST6gHJoJM5PPmlkMPAi8BNwV4rrfLc7K2p3PBoz+iljzHChD8hRIwb3w3xDkYllLzfIs53aZYZrkJn6YR6y5w4B3wGXI84aC7xpz1UTlLsTGLL76np77ffA8FD6YGTlQeCzqxyV1SUGNxpJog/wVyOG1vzkk0buR1bgTAylVWg8zNcNdocxlKe+wxg0ZiTFoDGjHq70Aa1rpEpomC9uawSXTEMKt7pBnhHI5LMRwGyk9bwR6YoMCCaZDUBaiuuAv5EVCrOAX4Ap1L6kfSPOVy8hQjgb+BpZYbAYWQExHREASIDylWYaSaIP6H6N1POTLxqZgPRGPUa67nPf7c6K2l0fjRlCGWOGK32AY43k3ZiahSxXfCNB3h2Ig6Yhe0KsDJ3bbY9bgM8j1+0H3rI/T8hc0s4T56sjwDVIt+V2ZF+euUglTwJ22Xy/t6+YzkmqkUb6gO7XSD0/+aCRYHjvW+DelNf6bHcrqN3N0ZhRvpjhSh+Qo0YMbof5BiHjt2vSFgRpRR6lf2nnbPv7hjr5g+XVd2b4rrwx5OOrY5EK30exlvgGGPLTSFQf4K9GDPn5qSgaGUZtt3ujz+LQdb7bHYeh++s7DoPGjKQYNGbE4VIf0LpGqrRpNd9kpDXYqDuuHmfYY9AV2YuMa45BZvEfjOTvscdtGb6rCGTx1Rxk/PgF4L88CtUGsmokqg/obo1k9VNRNHIAeK7OuXHABcD7wDfUDgH6bndW1O70aMxojuqjTTHD4LZnailS8LiJX+fQv4QzzED6N9j6IHJulU1/IJJ+BdK1uQf5CzjMCjq/SZmhNV+dGJM2Htkqfy+1KzQCVuC33Vn0AX5qxNCaPsBfjUDjCejdaLeh/fW9gs7XtUFjRlIMGjPicK0PyKaRgCoZe6Zm2g/0F3oiUgkg3W/zI/k/RMYto1yJdKH1Irsg70IcNAWp6O3ATZFr5gEXAXcjLdT1wJnIGOphmz+6b0YwJ+xQc/Oc4tJXa5Gu2S+RB2EsMpHwANJNGbefiO92Z9EH+KMRl/oAvzSShm6xu9P17fv7ADRmaMxwrw/IppGmGBr3TFVoPNdhWyjveJs2r869eoAnkRUGfUgF/omMXVaoXdIZZjjwOLAV6ZLbBbwOXFwn/0bknyeeXOd8XlRw56vbkR1r9yAPw1bgaWBUg+/33e6s+gA/NFLBnT7AL41EqRDfM9VNdlfobH37/j4AjRnbQnnLGDPy0gek10hAlVDPVBhD8mG+Ziy094rbOKxdDENalkX6Z5ZxuPZVWe3Ogg++ysNPandx0fdB5yirr3ywuwj6iFKlSWMq+Gxu4Us2IS3ETjID2RY+bhy1SLj2VVntzoIPvsrDT2p3cdH3Qecoq698sLsI+gBZDRjtPQNks6qAUdROPusDluVfNkVRFEVRlMJzHHBHJK3SgXIoiqIoiqIoiqIoiqIoiuV/+NxET1RKZ7IAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left [ \\frac{1}{36}, \\quad \\frac{1}{9}, \\quad \\frac{1}{36}, \\quad \\frac{1}{9}, \\quad \\frac{4}{9}, \\quad \\frac{1}{9}, \\quad \\frac{1}{36}, \\quad \\frac{1}{9}, \\quad \\frac{1}{36}\\right ]$$"
      ],
      "text/plain": [
       "[1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAA0AAAASCAYAAACAa1QyAAAABHNCSVQICAgIfAhkiAAAAHZJREFUKJFjYKACCGFgYJjMwMBwmIGB4RMDA8N/BgaGJYQ0XYAq/MzAwHCdWE2ODAwMqgwMDIwMDAwOuDSxoPH3EzKVgYGBgYkYRaOaBlwTeuQGQDEDAwODBJS2ZGBgWABlv2FgYChBN6SBAZJ0cOEH5LiMzgAA6XoX52TB9a4AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$1$$"
      ],
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADEAAAAdCAYAAAAD8oRRAAAABHNCSVQICAgIfAhkiAAAAi9JREFUWIXt1t2LTVEYx/GPGTTFNJkhk5SoKeRCopAkaoSkhD/BtXArLrwkF8pbjQvuKcqUuZIkbpTxFlFKTSF5yUthzMTF2jNzzmnvffY5Z585Ls63dmu113p+63n2evazFk3+D6ZUOL8d3xq4fi7swuZGLJwn5zC90U6U0lLB3NboGa6TL1VTSRBrcb9ejkwWJ9DVaCfiqGQnuvCpXo7UQtYgevAyYewkbmEIP/EZgzgs+87NxSjOVKObtU4fQD9exYwN4yGe4wNmYDVW4m3UHyqjvxd92IjbOeoW0Zcy1pbw/hj+4kIG/QF8FKpfxbqF6dSCo9hSYtQp/V/4lfD+StT2pNhCh7AD/UJKVaxbGEQv1mB3idFW3CzjSBzbo/ZJmXnbhAP0Wl66y4X8KwzuvGwF4CCO4DTuClv+GHPK2F3FD8npU5XuG6yL+tMUV4w03kcLjD0DQtVJow3fhUCq1o37wjewI+pvwJ0yjozRLVS7buzEIqEkrkix6cVMXM9Z1yYTpfRUtEg1LMBvPEuZczma05Gzrqn4giWylcc0BoUUmB0z1iqU1YFadePSaUSoRofwtIoFCpkXtaMxY+uFkzctlarRHWePEOn8MmKLhVwtpcXEoXQvwfZs5ETcz1+L7jjtsl279+GPcMe5KNx0L+F1tNA7LE2wHRJKZt66RczKMGeZcI48EvJ7BF/xQKjtnQl2qyJn9uesO6kcF4JY2GhHauGF8JWbNGlSB/4BwBaSc+fGZlwAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{\\sqrt{3}}{3}$$"
      ],
      "text/plain": [
       "√3\n",
       "──\n",
       "3 "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_s = sqrt(Rational(1,3))\n",
    "c_s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_x, u_y, rho, tau = symbols('u_x u_y rho tau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([f_next_0, f_next_1, f_next_2, f_next_3, f_next_4, f_next_5,\n",
       "       f_next_6, f_next_7, f_next_8], dtype=object)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_next = symarray('f_next', 9)\n",
    "f_next"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([f_curr_0, f_curr_1, f_curr_2, f_curr_3, f_curr_4, f_curr_5,\n",
       "       f_curr_6, f_curr_7, f_curr_8], dtype=object)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_curr = symarray('f_curr', 9)\n",
    "f_curr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}u_{x}\\\\u_{y}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡uₓ ⎤\n",
       "⎢   ⎥\n",
       "⎣u_y⎦"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u = Matrix([u_x, u_y])\n",
    "u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.codegen.ast import Assignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAACmUAAAAaCAYAAABMp8pfAAAABHNCSVQICAgIfAhkiAAAIABJREFUeJztnXu4blO5wH9sti0bWxQ5nI2EitouJRV9FCpU0v1UZhedDqXL6SLdVhe3qFzilFPZOaUbXXmE0iKVREjZLmHtsN0VKmzsdf54x3y++c1vjjnHvI4x13p/z7OetdYYY8455jvHO9533EFRFEVRFEVRFEVRFEVRFEVRFEVRFEVRFEVRFEVRFEVpnXnm9wCYTvxck3PNqcCdwBqt5ix81gKOB6aA5YjcPuwzQ4HQB7lsj+Trbb4zkoPqmdCH8tQ1fZFJ6HqmOjakL2WqS1Qm46hMxlGZjKMyUZpG/QmlTUIvX6BlLEbtyzh9kUmRnq3HaH/kdEf5cuEi4ATgSGC+57zYeD8is9f7zsgMok69G3+PNzSao37RB9saAqq7zaL+0sygi+/YF/8pNFRuSmgU+Rvqj7ih/kizqD/SPuorhIvKTQmNqr7C65F+wLOAV2dct3rZjLwTWGD+HpiHTgITwLss1+wArEAM9Wznp4jMzgI+i8jtqT4zFAh9kcsPgdsIs2Nd9WxIX8pTl/RJJqHqmerYKH0qU12hMhlHZTKOymQclYnSBupPhM2rkElLvwLuR+qAbzpeuz7wGNJp54tQyxdoGUui9mWcPskkT88eh+R9Aum8D2lS5jSwyCHdfOAWk/5HreZonNPMc7fo+Ll9pWjSZN1695vm/ltVvD401gXejujwX4AHgfuQCctvA1a2XBeybU2iutsP2tZbZchRwC+AmxF9vxe4HPgkUh+UYQJYXCJ9V9+xT/5TSKjc/LMR8HVgGfAw4jcfC6zjMU9t8iaGC7bebklT5G/0wR/x6YuA+iNlUH9EUF9BsaFy80ud/vm+0ravMAm8NyP8NcAmrpl8BaOGY4BkeKLgunOBv1NhBugMYytEXj/znZHA6JNcno3k9VDfGclA9UzoU3nqir7JJFQ9Ux0b0rcy1QUqk3FUJuOoTMZRmShtof5E2FyBfJ8HgCWU6/R5h0m/aztZcyLU8gVaxmLUvozTN5m46tkk/ZyUeYxJuwL4a6s5GuffkfKwUsfP7SvxpMktLfF1692rgX9gn6zYN96JyGsZ8C3gCGRCyN9N+Olkl72QbWsS1d1+0LbeKkOWAxcjen4kMrD7e0T+twIbF1y/MPH3BKMTLTYmv27s4jv2zX8KBZWbf54M3MFw4t6RwPkMT70sOxEqdDZG6oMHyJ9oUeRv9MEf8emLgPojZVB/RFBfQclC5eafOv3zfaQLX2GS7EmZKwNfY3giuZW5wDmpsAHFkzK3QByDk4seMAs4EJHXAb4z0iAR8k6DGvfom1yWAEuBOb4zkkD1bEjfypMLEfX0rI8yCU3PVMdG6WOZyiNi9tmyIiJUJllEzL76uIgIlUmaiPr6ExoR/Xwn9SfCZVfgKUgn/oBynT5nA3fj/7uGVr5Ay1iSmWhf6tJHmbjo2ST9m5T5NGRg6myGneDrtZwvpTpLkE77rMG+uvXuGsjuzxdVvD5EdgP2YVxeGyATB6aB/SzXhmhbk6ju9oc29VYZxTaweBiiIyflXLsp8p2OBdZmONFiHvBxZMee51uu7eo79tF/CoHZKreIcPpOzkHy8u5U+BdM+Jcr3jcinHeMWQn4OXADcDT5Ey2g2N8I2R9RX6RfqD8iqK+gZDFb5RYRjh2t0z+fR0Q47xjTla8wSfakTICXAh8qyujBwCGpsAHFkzKPNGleaInf1cQfA2yHrNi5FznW5IdIhw2Io3EacKeJOxNZjdEH9mO4DWr6p+9H00RUV6q+yuWTSB739J2RBHl6Nht0DPpbnlyIqKZnfZZJaHqmtkzoc5nKI2L22bIiIlQmWUTMvvq4iAiVSZqI6voTqj2JqNeQ9vVe6k/0gwHunT5rI8eunZII0/I1RMvYzLYvVemzTFz0bJL+Tco8H3gEOZbrVHPN7jnpdwZ+gHQeP4To5iXA4RXSxfXA5zKeMw/p+73aXL8UWfk/BxkI+2Mq/e7mXkcAWyP1xh3Iro+/AXZMpW+rDtqP4YT95cD1iXzXye9R2HXnjSaNS3/c0chOCj827zsNPN2keZ75/3jz3t8Ebs/JE8hgQgRcANyDHAF4ObB/RtoydGEjDjXPOMESH6JtTVJGd1311jVtU7pbVg+Sz26ybLjqbdk819VbGJbTfTPiFpq4H1iu7QNd+YPPNM85ryDdRsiuMbciRyf+AfgLIuM8P6XoO9alrv9UpRwV2Q0Xu+KbJvzOsrJzkUtXsosIo+9kM3OfmxifDLYmUn/+E1kcUpaIMN4xyXuQiVe7IPMVpsmfaFHkb4Tsj/hsR4DdH2m7HZF8tg9/pOt2BPj1R9RXcEN9hWqorxCeHYWwJmX22VeYxD4pcw7wJ2BBHJA1Y/89yAuX5UXIit+LLfHbmd9bIKuCH0Mq1r8ix6V/HVldewlyLvs3gOuAvRCHow/cBnwKuAt41Pz9KeSDX+8vW97pq1x+bX7nObpdk6dns0HHoL/lqU36LJPQ9ExtmdDnMtUWKpNxVCbjqEzGUZlkM1Ptia/3Un9i5rEXcopHsmNNy9cQLWNqX7Los0xC1LO6vAHpYD4JWVl/lQnfzpL+UOBCYHvgF8juRj8FVgNeXCFd/JzLU89ZA/glMtD4L+A48/8nkJ1I5mdcs635vQVy/Fxcb1wA7AT8BBn4Tz+7qTpoDvBt5DjszYHvI3Jdgey68vWa+b3MxIMMtn4q8XOBCXfpj9saGTxcAXwFGTS4xsRtb34vBC5FFh+cmpOn1ZHj3U5BOvK/Yf7eAPlOdY7a7MJGPGJ+P2qJD1nny+iuqz6WSduU7pbVg+SzmygbZfW2bJ7r6m3yfS/LiNvB/E5/hz7RlT+4j/mdnoiT5hbgbchufjsh3/t44JUM68osir5jXer6T1XKUZHdcLErvmnC7ywrOxe59EF20Jx+7mZ+n4u8a5IHEHv7OOA5jb9BMU3XQU9FJl4dh9hTF4r8jVD9Ed/tiOSzkjrYRTsi+Wwf/kjX7Yjk+/rwR9RXcEN9hWqor1CP2dC3O1N9hcfMddYFtdsis/3TkzUHyAzQCct1ayDKdJUlHuBb5h63A89IhK8J/M1cfxtSycbMRYS+Aodz1wNhDuKMFBmWvhFRb6ZzH+WyNvLOl/jOiKFIz2aLjkE/y5MLEdX1rK8yCUnP1JaN0tcylUfE7LNlRUSoTLKImH31cRERKpM0EdVlEqo9iahXJ/h6L/Un+sEA95W430f6JpLvquVL0DI2ZKbalzr0VSYuejZp0oTCNPadMtdEdve4G1jHhO1prvluRvr1Eb38FaKLadYrmQ6G9cCWqTTfNOEfR3ZijHkhw10r3pO65jsm/E7GB4PPMHE7Zzy7qTroSwx3rFklEb4q0tE9jeycUDW/AO8w4e9gHNf+uAewT3pYbNLc4ZinH5qw9OTLJyC7Zj7EsGyVpW0bsQoiq7zdJkKzrTFldLeMPvrQ3Sp60GTZKKu3VfJcR29Bdpa7yxJ3hLn33jnXh05buv4BZJzwi0iZngauROqnPDYCvsro7lc3IN/WtluSy3dsgjr+U5VyVGQ3XOxKCNT1O8vKzkUuXckuIoy+k/hYzv+2PCeui/+rQh4jwnhHEDtyKXAtsnAF3Ha/KvI3QvRHQmhHQLY/0kU7IvlsH/5I1+0I8OuPqK/gjvoK1VBfIQw7mmRAODtl9tlXmMS+UybAfyL1WCYHIzP50wzIn5S5hYk/N+fBS0ya9GoLkBnM0wy3ck5ynolbN+feIbENkt/FnvPRNBH1lKqvcnkQqQhCoEjPZouOQX/LUxER1fWszzIJRc/Ulo3S5zJlI2J22rI8IlQmWUTMzvo4jwiVSZqI6jIJ1Z5E1KsTfL6X+hPhM8Ct02ce0jH2/VS4li9By9iQmWpf6tBnmRTp2STybqEwjX1S5hdM/EGJsA1N2HUZ6V9g4r5W8EzXdCD1wAOMLrx/jrnedkLSUrIHGq814XtlXPNZE5fcNaDJOmhHpLPdlud4EPQtNfIL8GUTvsPYFe79cW+1xIMMQk0DL3HI097m/+9Z7vUNE1/1eL62bcQxJt1ZBelCsq0xZXS3jD760N0qetBU2aiit1XyXEdvH2/if2aJj99pQ0t8H2hL129nOPlmGjkOdv2CvGyKlOvjkIHNCcRXmYfssHY/8PyM61z83iao6j9VLUdFdsPFroRAHb+ziuxc5NKV7CLC6Ds5mfzJBoeZ+I9UyGNEGO8I8Glk16nkpIwJiidaQLG/EZo/EkI7Asb9ka7aEfGzffkjXbcjfPsj6iu4o75CNdRXCMOOJhng1j/vQkQ479i1rzBJ/qTMHYGHkd2ox/gyshI3zQDJ8ITlpjuZ+KxVGiCz1B8DbrTE34us8l01I+46pKJNciAyM/ohZLvatJH3yf6ILA72nI+DkA6/+83Pb8l2IrKYYtRYFv0sdrhnKHIBeBLSgXkXokRXIQ5pFrdiP2qna/L0rGkd2wXZtn6ZeearKuS3TUIpT3XkNEWzehaKTD6CbOt/P6JjP0W2AM8jFD3r2pZVkVWXhFKmoJrdn2Jm27KYQ5E8fckh7RQzVyYTjOfftUNtKuPamVAfT5Gd/xNrXNt3mcwBPsOwPrkJ6VBbJe8iwxTNyaRpewJh1JNNvlcVH2u2+hN12n1dM8Ct0+dlJt0bEmG+/bBQyhd0X8ZiyvgcXRGKfYHq/VRTlKuLi/QnJJmUtbtFejaJvFsoTJM9KfPpyNHNf0JkkORuZGBwrVT4esDfzT1/AryW7F0QXdPF9UB6Vf7/mWuflfVCyE4kKxg9km++uddSRnfEiTnN3HOz1LObqoPiHXlOQ3z+9M/pJj4eUCmb35hLgOVkd5i79MfdwfjJUzHzkDJxgyU+naefmP+3t6SPB+v3sMTn0ZaNiDnY5G0JMoiWR0i2Fcrrrqs+lknblO5W0YMmy0ZZva2a56p6C3LM5TRil7K4B9HrLpmiOZ+gbV0HmVyxLzKBZRn2Y3VjFib+nmC0fbkx2XVo0XeEZuRW1X+qUo6K7EZRfJ3+5CnC8TvLys7F3halmWD8HV36Dqcyrgul76RoUubhJv6QgjxNEe47PhvxFz6XSjdB/rvHFPkbIfkjIbQjINsf6aIdkXy2D3+k63YE+PVH1FeYub4ClJsLk2SK2e0rTJH9nkVjarbrQrCjaQa4fbs0U4T7jj58hUnyJ2VuYp69FYx3hG6IvGBZHjS/bdugPhMp2OdZMrQO8APE2UgyH3gyw7PaQRyF45AO74uQbdfPRraX/mv5rDdObFAu7+h5qzIuN4BbEEf7esR52B9ZDbI9xVsEHwssSIUtAl6OVOBTqbgrHPLZpVxsMlkVUfBfI2VnL8QYbYZsRZ7F6gzLt2/y9KxJHQOR05XAKci24KERgp7F5amqnJrWsxBkAuIknYR0Dq2ErEz4OVJH2+xLKHrWpS0DcXzKyqpLQqm3X0k1u993Wwb2uicOew5wAO5b/89kmYA07geJ8Mcc7zlT6+NnMdqBtzVSf6V3nstipsrkY8gEtv2RjohnII3Dh5FJI3k0KZOm7UnV9lHT37nJ96riY81Wf6JOuy9U9kU6tJM7bPn2w0IpX9B9GYPyPkdXhGJf6vRT3YBMWHRlWUF8KDJZFfgg5exuSHpWhxOR/tb3Me6PXoXUP4uACxPhdyO7fnwSeCmwj7n2POCjyCBnmXRxPRD/H7MH0ql9qSXvT0LsyQOJsEXmXuciHcpptgPuQybdJp/dVB0UTzx8vSXPMbGulc0vyPfaBrgaKZ9pXPrjzkIGorN4pnmGbfeWdJ5egAwwXGZJ/yTzu0o/eBs2IuYgpC68GtnFs6hfIzSdL6u7rvpYJm1TultFD5osG2X1tkqe6+gtDCc9Z+nZZsikYtuuPG3RpE/Qpq7H3IFs8PIHZGD0VPInBy7NibvZEl70HaEZuVX1n6qUoyK7URQ/oHp/ckh+Z1nZudrbojRV+g5D7ju5z/xe2/KstVLpbIT6jqsgkwGvQ47MrkKRvxGSPxJCOwKy/ZEu2hHJZ/vwR7puR4Bff0R9BTt99xUWUG4uTJLZ7itUHVML1Y42SajvGKqvEPteGwLXpCN/AZyQcdEAMUATlpvGW2dfZIl/l4k/ICNuPxOXtX36zibu2ETY74D/TaW7HjjC8uyu+RXjq0FiNkI6oG9HKrQ/Mdzmendk4Ck523dj5P2fnLh+Gtkx5ELEwL85JzzNvcj59VWIzDMGFa+vIpc2ZHI47hXRyibPtpXsXZOnZ03qWJppwtsp04eeFelYE3KKqK5nodY98eqyfSz5DknPurRlWRTJqmvaKlNly1OTdj+ie1sG7dQ9ayN6sxuyIqfqrlURM0MmE+b+TREx8+rjY4G/kL3S14WI/stkBeMr9L4BnFnhnaC6TJq2J6HUk23ZSRcfS/2JUeq0+9pkQPFK3DnIYMHZqXCfflhI5Qu6L2NN+RxtEIp9uYlw+qnyZALd9r+cibvdddGzSbIHx3wxzfhOmW/EbceAvNX0c5HFjt8zae8me9eVvHRxPRAl0s8zYbbBkK1M/LdT4fHOh+/MuGZNpO78ZSKsyToozvMFljxnUTa/IBOGp4GvW+7p0h+XZ3cPxL47QzpPa5q0V1ruNQcZ+LqDan59W/b0vSb+KuCJDvkIzbbW1V1XvS1K25TuVtGDpspGFb2tkuc6egvwHRO/MCMulsXhqfBXIzYuec1xiM1/guU5vuiizZEkPl5wvZLXFVH0HZuiqk9ZpRwV2Q0Xu5LEZ39yVblBedm5yKUozQTN9R1GhNF38nYT9hXLs84x8S+skM8I/++4ADf7bKu3ivyNkPyRUNoRMO6PdNWOSD7bhz/SdTsCyteFtwDvT6VbhPgnT7M8w4b6CuXok69QZi5M2/TNV0hTZ0wtwr8dzWJg0vg+vrzvvsIk+TZ5FfPMveN/kiwnfya5jduQWdZbWuLjWdBZM5nzZjlva37HqzHmmvTHpNKdCzzXKaftshIyqze9GgSkM/u3yLu8CpmN/nxGV7Bezehs30Um/sbE/wAfQI4Oux74G8N3T4fHzAFegzTQflP15WpQVS4707xMDkEG+E5DnM5lwFeRFUjTqbxtafLuslNYF+TpWVM61gd86NnzMsKSOuabUOseEEdrZewrdUPSs65smY0iWXVJm2UqPlLVpTw9gDRoQ7D7odU9JyNHbJwPfKLWm1UnJJkcjKyOuxXxZy9GHPapWm9YnlDr47lI51581GGXhCST1yIda1shq9Oehkwy6nrSTJP2JKT2kU9/VP0JwXe7rwl2AdZFVtIn8emHhVS+oPsyFoLPkUUo9uUm5Jje9A54ofmr0H3/y0a4293Q9KwKawFHIzKMj/dLswkyKL5tRlzMcmTHq58jAxjPR459S++KmJcurgeS+vyo+bEdU/ihjGtI5DVrV5xtkbozWac0WQfFAx9lBu7K5heGZds20OzSH2fbNSiZxvbeyTzFC4ls3ylCdiI6imp+fRv29MPAkYj+7o4M7BcRks43obuueluUtindraIHTZWNKnqbvI9rnuvoLYhteoTxHZlWYzhImy7vpyPl/WPI4OEHgNchdvAuy3N80XXbbEPz2/XEEleKvmMT1PEpq5SjIrvhYleS+OpPriM3KC+7MvY2L00IfYdN6mc8QWwPhpMGYtZE+lEfRN61S5p6x4eBr+U8Y1tkIta1SJlLU+RvhOKPhNSOgHF/pKt2RPLZPvyRrtsRUL4u/C3jR8h/EVkkerXlGTbUV3Cnb77CK3CfC9MmffUVYnyOqc2GeT8z3VeI51wuz4r8FvDdjPABUtgmcm58ukmzeUbcFYhA5mbExSt1sgzjYhMXb2Ucz3bfJZXuE4ggbddHOfluki3N876TEXc2Mit9Zcu132J8B4FPMDqz/6PAv5CGAw7h2wD/QJylvyFbklclovpM56pyaUMmD5mfIxAlfAsio4Mynv8Wk+93ZcQtptuyFWPTs6Z0LItp7LsTxddHOdc3jQ89s5WnJHlyciWimp6FVvck+S7SIJljiQ9Nz7qwZTZssoqvjwqub5I2y1SZ8lTW7hcR0b0tg+brngMQZzQui5P42SkzJJm8BFk9tQ3S0PwF0iB/fME72IiYWfXxaxB/cENLvAsR/ZfJyogPuAJp0E8Dny16gRwiqsmkSXsSUj3Zlp108bFmuz/h0u6Lr4/yMt4yA4pX4p6AdNCunwr35YdBeOULuitjZXyO+PooJ99NEop9aboerkOeTKD7/pcydjdPz2Imye4IX4wfPZxmdKfMY03YUTnX7GDSXJUI25bhbqRJNkfq9KWILF3TgdQDDzK+GP/P5vnp3ZIOYrhzQDruSux1yvvMNW9IhDVdX19pwl+ZcQ3IwE2y7i6b32T4/pZnQLX+uJg/mDSrZsRl5eliE/aiVNoXIvb+JrKPKV1MsS40/X0+bsIvpVzbKyTbWkV3y+ijD92togdNlo2yelslz3X0FuASE7dFImwNxB7GMs36bnsgNu0Q4H6kbKRZjH+/u2ld3wrYICP9ysBh5pq2dmTK+45NUMenrFKOiuyGi11JUtT33hZ15AblZecil6I0TfYdRoTTdxJf8+5UeDyJ5MsV8ghhvWMWEyZt1k7kMUVtDFt8nI/IIR9NEFI7ArL9kS7aEfGzffkjXbcjoHxd+H5k176Y/ZBJ+eum7ruY7tsF6itk48NXKDMXpk366CskqTumFhGmHR1Q3D/vSkSY7xgzQXu+wiT5O2XG/cU7wHgD/wZgx5yL8zgDqfz3ZNQgrIashv8j2TNBt0MMftYq2u0Qw78kFZ7uhF0pIwyGivxoXsYbJGs1K8i2ui9GZLuCbLZl/LirRYweV7MIOebpxox0WeHXmrgFyLc5FVGKJo/3dKGqXNqQSbySJd7u9nKkMj8IWSGQZA9kQPDHGXnrumzFZOlZGzrmig85+NAzW3kKhdDqnpijkcHJnbGvfgpNz7q0ZUnyZDXT9KxKeXK1+20SSt2zJbIl/85YVth0SCgygdFjbq9CVjzdiDQKvmB7gRYItT5+GyKjZbaMt0hIMnkV8B9IB9qfTZrjkIFt22q6pmnLnviuJ336o6D+hEu7z1f75RXmB4adszsxnPx1N7LLUDL9b5BjWWN8+mEQXvmCbspYWZ9jtve/+K6HwS4T8NP/Usbu5ulZEb70MMk2SL/SzcCnc9L9CXnPpwKrIzp3MDKIdwkipzuBTYGXmWveinwz13RxPXAF4zI5Atl95yxkEsftyGDkU5DdTLdktPysZvJ6Fdn1QHr3gjbq6w+a/J6B7ObzR+Sb/5t5/qrAv1fML6n/D0M69v+JyPj7iTRV+uNABhi2RnTpkYz4rDx9FBmEOBM5VnIZssPInoj8dgfuy7hXkS40/X32R8r7Y8guSwdnXDfF+IRvCMe2VtVdV32kRNqmdLeKHjRdNsrobdU8V9XbmHOQnaYuQHZon49MJvkjMknrcWS3r88Ffo8sMtib7F12fNulNuriFyM+84XIGOI9yCKqFyCLMm4n++jBJsj7jk1Qx6csW46Kvo2LXUni0vfeFnXkBuVk5yIXlzQh9B22oZ8HIm3o4xEZLkHkvytwHeJXdInv/qEkRW0MW3wffJE22hFg90fabkckn+3DH+m6HRFT1o5cDHwemUz+T+BzSLm5J3XfrtsFoL5CSL5CmbkwbdJHXyGJrzG1NvSzbP9828wEX6GIBeZ35rHoewG3ZIQPkM7kiZwbz0Uq1N+lwrc3134l45qFJu6MjLjVkM6q5LbqcxED8upU2hMRhUxzObJi0baldtMcTfZK5pcj+bbN9l7dxO+aCr+ZUQN1HeJgp7GFp/k5sj1xFSKqz3SuIpe2ZLKUcRm8CXFekqyNKPSPMu4B3ZetmCw9a1LHspjGvjuRDzn40DMXHcuTkysR1fQsxLrn88hg+tOsuQ5Tz7qwZWmKZDXT9KxMeSpr94uI6NaWQfN1T2Ty8WjiZxppvDyKlLkyxPcblLwOwpGJjV8C/+OYNk3EzKmPFyKNhpdbc+1GRP9lcjPwnlTYx7A0jByIKC+Tpu1JKPVkm3ayyMdSf2KcrHafr/bLBMPVy1k/U4m0zzJh70/dw6cfFmL5gm7KWEQ5n2O29r80XQ/XwSYT8NP/4mp3i/QsZpLsia6+9HCa4U6ZF5r/93W4bolJGy+KfwUywHkt8h7Lkbrxq8ggJyXTxfWAbVekdyPfYDlyROGJSP1wD+ODJXl1SvwuDzAsV23V189Cdj+53aS7Gxkw/QqjO/KUzW+SdyGyfcjc47BUfNn+uLp5egFwPrKjyT+Q9/0UchSpjSJdaPr7TJBv46cRvU0Tkm2tqruu+kiJtE3pbpUy14buuupt1TxDNb2NmYcsFFiG7PZ8KXIM4gLE15m05GU3ZBxhBcOB5jS+7FJMG99za6TMXYF8y0eRyeG/R+qCqqeUuJD3HZugqk8J5ctRUVl3sSsxLn3vbVJHblBOdk3YWxtV+w4jwuo72Rg4BZmkshwZBz2OeroZEdY7ppkgf/erIn8jL74Pvgg0346AfH+kzXZE0TVt+yNdtyNiytqR1ZCd5fZAdu2+huzd+LtuF4D6CjZ8+Aquc2Haps++QhNjahHh2NEJ3PvnyxARzjtmMUF7vsIk+TtlvpqcE5TmIwr5hFT4gOJJmSAzrqcZnuXeBr8DTk6FXYes0kiyAFGWz7WYF1deilQOa1nit0bk9qRE2C4m7Nnm/zWQ99kpda0tPIvzaWYr2qbIk0tbMjkNWb2d5DPA1amwd5tn7ZxxD99lqws9S2IbCPcthzRt6ZmrjhVNGPCBr7rneNw6hULVsy51rEhWM03PqpQnV7vvi67rngXmnsmf3yP2bWtkVybf+K6PQRpstyHHboaAT19wApFFeod83/iQyT2MHzPwEapPygyF0OvJuhT5WOpPjJNu94XmT9g4HJHZph09z8VnDbV8QftlrIzP4VsWSXzYlz7Uwz76X1ztbp6eJZlkfFL0fHdFAAAFGklEQVSmz7I3zejx5X3mTcj7dLkzQh/xaduLCKkeLiJk29o3VHeLaVJvnwn8Hdlt7EeM7roXo2W4HXzUv0U+pU9c+959ELLc0oTWd6g0S1Ebwxav9Xg51Bdxo0k7djFwAjLpcu+M+NlahtVXGOI6F8YXocotyQRhjqkpzVLVV4DiSZmHA8fmPfwUxlf4DxidLXuN5dp5yOzrn+Y9oCavRVZhvB3ZSvpYZMXwwlS6fZBVCRvgn3WBe5FK8OnAVshuA/FZ9xsglU+8EmZHRMaPITsWgHRwP4Z0eCexhR+JFJBNkO3PjzDPeEkD79MUeXJpQyYgq2IeQbbs3xwp6/chWzbHrI7Muj/dkm/fZasLPZuPDC4sYuhQL2L0aBnfckjTlp7llScXOfnER91zEtIY2M3cP/6Zn0oXsp51oWPgJquZpmdV6m1Xu+8LH3VPmkngS1VfoAV8yOQYZDeZTc39zkT0a6aXk7xwkNWGSxG/MDR8yGQxcirAXoiPvC9wF7LDRJ8JvZ6sgquPpf6EW7svNH/CxhJkRX0XuPhhIZcv6K6MJZkk2+fwLYskPuxLH+phH/0viym2u0V6th7jK/iT+Cx7fZuUOQd4Ykb4i5Dy+lfG2+3KKD7qXVdCqofzCN22hojqbj2a0tuFwK0MjwF+BmI70wN0WobbwUf9W+RT+sK1790XocoNwu87VJqjyN/Ii9d6fBz1RerTpB37IuKDnGOJn61lWH2FIS5zYXwSqtxiQh5TU5qjjq8AxZMyL0LajFY2AX6QETaR+EmvdE+yC/BJ3CYRVOVAZBvVh4HLzDND53nArxEH5W/AuYw6MYcgFdCtwKmIDJck4t9J9mRYW/hipMJ4GLgTOcJuzzov0BJ5cmlaJjF7AVciTsl1wMGM7vLxVKScb1LmRTqmbT0bkL198eKWntcUbehZXnkaEL6cuq57bFtfT6TSha5nXdgyV1mFRp0yVbXeDt3ud133pJkkrEmZ0L1MvoM4yMvNPc8gvB0Duq6PQY4ymQa2qJPxFulaJmsiE2WWIscO3IisWJtX5yUCIfR6siwD3Hws9Sf60+4LDRc/LPTyBd2UsSSThOdzZOHD5vahHu66/8XF7hbp2eMY7Y+csKTzQd8mZW6D9IP9GPgCsrvKr5D3uIswd38Mka7r3ZlGH2xraKju1qeu3j4esYnp4/W+hwy2Kd3go/4t8il90If+5BDlBv3oO1SaocjfUH+kHOqLNENTdmx/5Ojnp9fO0cxDfYUhRXNhfBOq3CD8MTWlGer6CpPYJ2VuAXzbJRNHoQVNURRFURRFURRFURRFURQlTd8mZW6JLMK/FRmY+RfwZ6QPeH2P+VIUJR/VXUVRFEVRfKK+SFicB5zoOxOKoiiznEnskzKPRnapLWQu8BnCmjGtKIqiKIqiKIqiKIqiKIrim75NylQURVEURVEURVH6x8rIBNhDgduAdfxmR1EUZdYzSfakzB2A15W50XroCgdFURRFURRFURRFURRFUZQkFwMnI0e0r+k5L4qiKIqiKIqiKMrMZACsAJYAz/WbFUVRlFnNm5F+wHOB12bEb91tdhRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFUfrH/wPjymG5Nz1Z1QAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left ( \\left [ \\left ( m_{0}, \\quad f_{curr 3} + f_{curr 6}\\right ), \\quad \\left ( m_{1}, \\quad f_{curr 1} + f_{curr 2}\\right ), \\quad \\left ( m_{2}, \\quad f_{curr 0} + f_{curr 4} + f_{curr 5} + f_{curr 7} + f_{curr 8} + m_{0} + m_{1}\\right ), \\quad \\left ( m_{3}, \\quad \\frac{1}{m_{2}}\\right ), \\quad \\left ( m_{4}, \\quad f_{curr 0} - f_{curr 8}\\right )\\right ], \\quad \\left [ Assignment(rho, m2), \\quad Assignment(u_x, -m3*(-f_curr_2 - f_curr_5 + m0 + m4)), \\quad Assignment(u_y, m3*(-f_curr_6 - f_curr_7 + m1 + m4))\\right ]\\right )$$"
      ],
      "text/plain": [
       "⎛⎡                                                                            \n",
       "⎜⎢(m₀, f_curr_3 + f_curr_6), (m₁, f_curr_1 + f_curr_2), (m₂, f_curr_0 + f_curr\n",
       "⎝⎣                                                                            \n",
       "\n",
       "                                                ⎛    1 ⎞                      \n",
       "_4 + f_curr_5 + f_curr_7 + f_curr_8 + m₀ + m₁), ⎜m₃, ──⎟, (m₄, f_curr_0 - f_cu\n",
       "                                                ⎝    m₂⎠                      \n",
       "\n",
       "     ⎤                                                                        \n",
       "rr_8)⎥, [ρ := m₂, uₓ := -m₃⋅(-f_curr_2 - f_curr_5 + m₀ + m₄), u_y := m₃⋅(-f_cu\n",
       "     ⎦                                                                        \n",
       "\n",
       "                           ⎞\n",
       "rr_6 - f_curr_7 + m₁ + m₄)]⎟\n",
       "                           ⎠"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "moments = cse([\n",
    "    Assignment(rho, sum(f_curr)),\n",
    "    Assignment(u_x, sum([ (c_i*f_curr[i])[0] for i, c_i in enumerate(c) ]) / sum(f_curr)),\n",
    "    Assignment(u_y, sum([ (c_i*f_curr[i])[1] for i, c_i in enumerate(c) ]) / sum(f_curr))\n",
    "], optimizations='basic', symbols=numbered_symbols(prefix='m'))\n",
    "moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "const float m0 = f_curr_3 + f_curr_6;\n",
      "const float m1 = f_curr_1 + f_curr_2;\n",
      "const float m2 = f_curr_0 + f_curr_4 + f_curr_5 + f_curr_7 + f_curr_8 + m0 + m1;\n",
      "const float m3 = 1.0/m2;\n",
      "const float m4 = f_curr_0 - f_curr_8;\n"
     ]
    }
   ],
   "source": [
    "for expr in moments[0]:\n",
    "    print(\"const float %s = %s;\" % (expr[0], ccode(expr[1])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rho = m2;\n",
      "u_x = -m3*(-f_curr_2 - f_curr_5 + m0 + m4);\n",
      "u_y = m3*(-f_curr_6 - f_curr_7 + m1 + m4);\n"
     ]
    }
   ],
   "source": [
    "for i, expr in enumerate(moments[1]):\n",
    "    print(ccode(expr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAADnQAAAAfCAYAAABk38K4AAAABHNCSVQICAgIfAhkiAAAIABJREFUeJztnXv0Ndd4xz+JlyTiEgmNpqygREhISKREpRMiWq02IWpVSwZRligaFEG9S90abYmURbSVUm1d0rBSlyS1EEGIiKpKisTPLUEiLkGEyNs/9pyc85t3zjl79uw9e+/Z389a7zrvby7n7NnzPN/nefbMngEhhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCDE6hwCXAB8E3gPUnr9/R8/fNxYvAC4AfgRcCZwJ7B+1RSJVZCtCxEU+aMg13uZGzv0sX3FHfZcHOftnLHLtM/mksEW2IoQ78h9DrrEyN2L1s+x8uujcGqRh45BrP8tPhqH+M+Rq/7GJ0W+y2emic2uQHrmRa7/J7oUtshUh4iIfNOQab3Mj136WnwxD/Zc+ufpmbHLtN/mksEF2IsQw5EOG1GPlUzFzM8/GzNXcb9XGFbANODhAQ/YBHh/ge8fgLOAJGAO/F3AG8G1g95iNEkkiWxEiLvJBw32Aw2M3YuLknNeAfGUI6rs8kA72I2dNk08KW2QrQrgj/zEovwhPzJxEdj5ddG4N0rDwqK4qF/WfQTrTn1i6IZudLjq3BulRf5THiBKQrQgRF/mgQXlKeJTXlIv6L32kgf2RpompIzsRYhjyIUMuOUaFxVxNq40cuCXwn8DOnr83FrcAfgk8InZDRPLIVsrircB3gV1jN2TCHISJU0+y3L5kH3wdpqgX/plaXgNl+8pQ1HfpIh20Y2qaJp8UtshWykK12mpUZ9mj/CIcqeUkJdv51Cn53ErDwpGahg2lZD/xQcn9J52xJyXdKNlmp07J51Z6ZE9KeuSDku1e9EO2UhYaHw6PxpjtUZ4SDuU1YhH1X5pIA+2RpokSkZ2Uheq09ajOsieHHKMi4oTOU4EjPH/njL8GPgR8A7gWuBq4CHgJsEeg3/xVTD89cMU2e2Ic4nXN33sAx2FmPn+laesPgfMwThbrVa8x+i8WOdhKqnYC5djKDsATgfOBa4CfYo7zGcBNVux3MHADcELoBgrOAK7AJB7rmJIPHgOcAnwM83r0bcC/rNj+dsBHsW9zO24N8fkHAadjztN1zefZwMMt25I6U8trYL2vKK9Zjo3OhKKU2Az9NRD66WDbxsE9J8iNqWmaarU8SdFWZCfpoVotLKqzxqmzwN1vS66zUowTISlJ21M9t9IwQ1vDcq6RHofpn22Y8+kb1VXp+koO/Zdq34F0JiTLdCNVmw2F8p7451Z6ZJAehSFFu1cOkx66NjqMkmxF48Pp43OMORcfjHENW/fyLKe0vAbyyedT7b9QlBSfVauFozRNyyX2l+TfKdoJyFZSxEWbVafZozprOnVWRaQJnQcBX/D4fW1+jhGAfwJehTl5F2CO41vAHQP85jswJ2hVAvinTRtmr299avP35cDbgVc2bf5Bs/zdGEEbmxj9F4scbCVVO4FybOWtmGP6DvAPwMnA/7K+/8/GnKddRmhj6RyCOR8nWmw7JR/8XNOGa4CLsUtQ3gT8seX3t+OWq8+/qNnmSuAtwCswAwwXACdZtiVlppjXwHpfUV6zHBudCUUpsRncNBDsdbBt4+CeE+TEFDVNtVqepGgrspP0UK0WFtVZ49RZ4Oa3pddZKcaJkJSk7ameW2mYoa1hudZId8Scr2sIM6FTdZUhVV/Jof9S7TuQzoRilW6karOhUN4T/9xKjwzSozCkaPfKYdJD10aHUZKtaHw4fXyOMefigzGuYetenm5KzGsgn3w+1f4LRUnxWbVaGErUtFxif0n+naKdgGwlRVy0WXWaPaqzplNnVc22o0/oPAP4C4/f12bZ68RfjjmWN3j+vVdjZs3edc12HwCuYu4QD8a8vrY96/f2wNcxbX2Uv2ZaM3b/xSQHWwlpJ3Wzf+W4f062UuN2rEc1+10G3HZh+U0xWrat+e42+2CeFHFqz9/LgZphdhOKi4GvsTpxn5oPHg7cDZMcVdglKIcCl2JseB3tuOXi849u1p0D3LJjvU07UmdqeQ3Y+cpYeU1NXrHKVmdWUeN+zDnFZhh2rC4aCPY62LZx15wgN6amaanVajV5aVpMUrSVMWv6GsWCdUytVqtRnQXx/QfGr7Ogv9+qzkozToQkN20fQqrnVhpmWOyHXGukHYD/whzzqzHt9D2hU3WVQdfADNKZ8nTGhVW6karNrqKmnJp2CKmeW+mRQXoUhhTtXjlMepR+bXQoOdlKjXssmNr48JTxNcacy7W8sa9hg+7lWUaJeQ3kk8+n2n+hKCU+g2q1UJSoabo/Iz1StBNQrRaCmnFrNdVp/cmtzoK87hWGceqsCk8TOn8T+Ajm9aXXY4zjRXS/nvQOmNeR3nnVjwbiAOYd1ObEZt3RHev2btb9R8e6v8XMHr/nmt++NeZ1qW+xbOusPadYbj8GIfovVZYd65DjDGErQ+2kJswNoyH6byg1bsc6e0rE8R3r9m/WXdix7lXNuof0/L0+xOrPmmF2E6rdL2n2fdiS9VP3wQr7BOWbwGPWbNOnH5b5/I6YpPwnmFeX54RtbjO1vAbsfGXMvKZm3FgF4fPCddT4P+ZU87gaP8daYa+BsF4Hu2zcNSeIjWq1tGq1mnw0LVVi1Wpj1/Q1igXrCFWrqc7ajOqseHUWdPut6qzVxBzTi0Gq2h6CXPLFkjQs1xrpmZiLwocBWzHtXDehU3VVWnVVKHLpP+mMIWWdgXF0Q+Obc5T3GKRH0qMucs9jpjbmFwLpo0G2sp4UbaXGPRakei9Prn4Xst0+xphzvZZXEf4a9iqmdi9P7nkNKJ+3QWP8c3KMz20qVKstI3dNm1r+XaOxrCHkEvsW26NabTU149Zqqd7HMxTVWZup8aO1FdOps6rmu26cq9mVCKzjZcC5wG7AG4E3A7cC/gp4Ycf2jwF+CHzV4beG8ojm8/Md6+7bfHYlb7MOuqi1/HXAn2Bm/H5xzW//LnAz7J3uF83n9Zbbj4Hv/kuZZcfqepyhbCVFOwH//ReT2zefl3Wsmy27L0YDFzkCU0SdH6hds9+FvPoTwrX7483nQzvWleaD6/gs8EdrtunTD8t8/lDMQML7ge833/k8zI1zD7BtbAT65DZTymvA3leU1wzLC2NQUh5nwzod7LJx15wgJqrVVKvB9Hw8Vq0mO0nPTkLVajn2BajOSgHfdRZ0+63qrNXEHNOLwdS0fRW55IslaViONdI9MBeFT8Zokg2qq1RXQVr9l2Lf2VCKzsB4uqHxzTnKewzSIzukR3nlMRrzW4/00SBbWc/UbCXVe3ly7EsI224fY8xT8EEbXK5hr2JK9/JMIa8B5fM2aIx/Tq4xxRXVanlpmvLv9ZTk37nEPpCtjIGLNus+nv6ozrInqzpri+2GDX+OSR5OAp6PmR0K8CaMcZ2AeaXoDQv7HMF4DvMc4BaYWbMHY55o8XnMBfs2B2Fek/r1jnVdDvMGjKEfBVzNXHx+3PxrczRmtm3X7Pk2W4DHN///oMX2oQjZf6lhe6wuxxnKVlKxEwjbf7G5qvnserrNXRb+vy/zRGJX4EDM66x/Eq5pWfYnhGv3Bc3nYa3lJfhgXz6HSRB2ZHOMXmRVP9j6/P2az+9gEqJ7tdafCxwDXNmv+UHpm9tMJa+Bfr6ivGZYXjgGJeVxLqzTwS4bd8kJYqJaTbXajNx9PJVaTXaSnp2EqtVy7AtQnZUCQ+sssPNb1VmbSSVOjMXUtX2RHPPF0jQstxppC/A2jJ2caLmP6irVVTNS6b9U+s6FEnQGwupGSjY7Bsp7pEehkB7llcdozG97pI+6NmrL1G0l1Xt5cuxLCNtuH2PMOfqgCy7XsBeZ6r08U8lrQPl8F6nkhWMx9fg8BNVqeWma8u/tKcm/c4x9IFsZi77arPt43FCdZU/2dVZF6zWeDXthTuon6H6z58XNfnu3ll+NefrxGHy7acPs3weAPTu2271Zv8y4zmnW77WwbNuSf1s79t8ZuAZ4l2W7/6b5rvdZbh+KkP2XGjbH6nqcoWzFh53U+Hl1ccj+80WN27E+ttnvK5hjmLEFOJ35Mf/Owrp9mmVnuzXVipj9WeNuN6HbfS3GHhcpwQeh3yvE/7DZdv8l69f1g218eGWz/nrgy8BDMEnNfhgb2AZ8xKK9Y+GS20wlrwF7Xxk7r6kZL1ZB+LzQhprhx5xLHlfj5/xW2GsgrNbBZTbukhPEQrVaurVaTfqalhop1GoxavoaxYJ1hKjVVGd1ozprnDoL7PxWddZmUogTY5KLtvsgx3yxNA3LqUYCeCnm6b6LT0DdimnjcR3bq65Kt67yRY79J50xpKozoXUjFZu1oaacmtYHqZxb6ZFBejQOKdRyymHK1keQrczIdUw0xXt5cvW7Mdo9ZIw512t5MM417EWmeC/PVPIaUD6/jBTywjGZenxuU6FabZGpaNoU8++a4Xafg3/7IsfYB6rV+lAzXq2W8n08Q1CdtT01fnKMiunUWVWzvD1X02qjZzfLj1my32ea9XdYWHZb1ifEGyw3pK5/NidiT8zM2f8DLmf++toZRzTf9bIl+38PM1PWld9vvv+xFts+o9n2YjaL2DI28N9fbWL334wN4h7rGMdpayt97QT6999pvVpuiN1/Mzbwd6w7Yl59vA0jyKcCrwW+gAmEX2rWHbmwzwOaZe/w2Ma2Pafqd6et+b7Q7f4Ww17lnbMPVthr30ObbY9ast62H9bFh5Oa7/klcEBr3S7AN5r11q8SD0zf3MYmrwH/8St2XE4przmtV8sNsfuvzQZhjzml490g3LFW9Mv/VungMht3yQlioVrNHmmaajUbQtoJKBa4HmuIWi1VnzltzfepzgrnPxX2uuWrzoLVfqs6q5tUxqQgflyEdPUsx2OVhhm6+iGnGukQTCw7qbV8K6aNXRM6VVfZk1JdJZ3ZzAbSmTEJpRttYttsmw3C2lkqx7uB9GiG9KhcPdpg3DwmFZuHcey+r27E1oxFNpj2sSrfTXNMNMT4cN825nIfzzq7GaPdQ8aYU/LB03q1fJxr2F1M6V6eqeQ1kE68SjGfB43xTyU+t6mw7yPVavloWip6BmlqWq7+naOWpZQn5th/Mzbw6xN9tTmV+3jA73kvuc6C/O8V7iJEnVU1y26cq7nFoiEz/gDzytH3L1n/a5hXj17eWgbwoxXfeynwsx7tuHz9JnwHOAPz6tIvAW9l8wzbg5rPCzv2vQvG8Ia8GvZo4Oesn5F8POYpGl/EzMi92uK7Q/RXm9j9NyP2sY5xnDa24mInYILhbq1lB2J8+Z8xwrnI5yy/d5HY/TfD57HegBHrZwKPa/79AvOUnGOBvwfuBnx3YZ9rm8+dV3zvUHseqz99203odu/CvP9dyN0HbZnF4WVP9rCNW+viw/ebz8uA/27tey1wFvAkzE11n7RpeGD65jY2eQ34j1+x43LIvCZ2rIJx4xWEP+aUjjcXHVxm4y45QSxUq9kjTVOtZkNIOwHFgpRqNdVZ3ajOssNXnQWr/VZ1VjepjElB/LgI5eQAkEYeUKqG5VIjbQHehrGfF/fYT3WVPboGJp2xYco6MyOUbrSJbbNtSqlppUcG6VHZejR2HpOCzYO73YfWjdiascjUj1X5bppjoiHGh3O5j8e33YzR7iFjzLlfy+uDyzXsLqZ0L89U8hpII16lnM9rjH8a8XkIqtXy0bQU9AzS1bRU/LsELVOtlmYs6KvNqdzHA37Pe8l1FuSTY2RRZ1W0Zn1iLoD/jOWzgu/X7HNma/mhzfIn9WmAZy5q2nDbhWX/3izbu2P7pzfrXuH4ezcBrsK8WnUVz2p+53+AX3H8rTEYu/9i0j7W0MdpYyu+7aRuvq/y8F1txu6/ddT4P9ZZIPwpcNOF5Xs1v3Wex99qE7M/a9z7MmS7d8Qkgpc67p+7D1bYP3Fi32bb53ass41bbbriwyObZRcs2efVzfrn9/ytELjkNsprluPTV2rGi1UQP15BuGNO8Xhr/BxrRb+n7izTQVcNXJYTxEC1mj3SNNVqNsSq6WsUC4bgWqupztoe1Vlx6yzY3m9VZ60ntTGp0KSo7aFIMV+UhnWTUo20G/ZPf31ts4/qKnt0DUw6I50xxNKNVPOAmnJq2lBIj/pRIT2aMaU8JrbNg3KYFElRH0G24kpNmJwh1r08ufpd6HYPGWPO/VpeRdxr2JD3vTxTymsgfrxKOZ9vozH+6cTnCtVqM6akabH1DPLRtNT8OyQpxj5QreZKzXi1Wsr38QxBddb21Pixq4rp1FlVs+zGuZo7WjZgP2An4GZL9nl283lqa/kOzec2y98JwWxm7S8Xlu2Lmfn9tda2OwFPaf7/WcffOwzYAzMbdxnPA16DmVl8OOk8DaQLH/33aOA6NgvUycCXgdt5a+lw2sca0k5gva3kZCfgp/++CZzQWnYgxn7u6aeZg3gc5okQ78Qc24wrgCuBuwf87Vz9rm+7+9jA3TFxxvUpDVPzwVXc0Hzu0LHOJm510RUfzsW8zv1umJyhzezJFBs9fysELrmN8ppucvKVsfsvNr6ON/X4bMMyHXTVwGU5QQxUq9kjTTOkmDOuY8xaTXZiyC0WuNZqufqM6qw0CFFnwfZ+qzprPWOP6cWm5BwA4uYB0rDlpFQjXQf845J/FzXbnNf8PXsaquoqe5QvG6Qz65myzkA83dD4ptvx5pD7SI/CIT0y5JDHKIexR/ooW7HFR/+lPjYM8e7lydXvQl8jGDLGPDUfXIfva9iQ9708U8prQPl8HzTGP834vA7VaoYcNE35tz2q1WQrtpQSC7q02fd9PJCGz6jOSofs66yK1qxP4InMn1j88Nb2xzfLu2aqHtise1bfRvRgX+D2Hct3BF7e/P7HW+