{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from lbmpy.session import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial 03: Defining LB methods in *lbmpy*\n", "\n", "\n", "## A) General Form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The lattice Boltzmann equation in its most general form is:\n", "\n", "$$f_q(\\mathbf{x} + \\mathbf{c}_q \\delta t, t+\\delta t) = K\\left( f_q(\\mathbf{x}, t) \\right)$$\n", "\n", "with a discrete velocity set $\\mathbf{c}_q$ (stencil) and a generic collision operator $K$.\n", "\n", "So a lattice Boltzmann method can be fully defined by picking a stencil and a collision operator. \n", "The collision operator $K$ has the following structure:\n", "- Transformation of particle distribution function $f$ into collision space. This transformation has to be invertible and may be nonlinear.\n", "- The collision operation is an convex combination of the pdf representation in collision space $c$ and some equilibrium vector $c^{(eq)}$. This equilibrium can also be defined in physical space, then $c^{(eq)} = C( f^{(eq)} )$. The convex combination is done elementwise using a diagonal matrix $S$ where the diagonal entries are the relaxation rates.\n", "- After collision, the collided state $c'$ is transformed back into physical space\n", "\n", "![](../img/collision.svg)\n", "\n", "\n", "The full collision operator is:\n", "\n", "$$K(f) = C^{-1}\\left( (I-S)C(f) + SC(f^{(eq}) \\right)$$\n", "\n", "or\n", "\n", "$$K(f) = C^{-1}\\left( C(f) - S (C(f) - C(f^{(eq})) \\right)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## B) Moment-based relaxation\n", "\n", "The most commonly used LBM collision operator is the multi relaxation time (MRT) collision.\n", "In MRT methods the collision space is spanned by moments of the distribution function. This is a very natural approach, since the pdf moments are the quantities that go into the Chapman Enskog analysis that is used to show that LB methods can solve the Navier Stokes equations. Also the lower order moments correspond to the macroscopic quantities of interest (density/pressure, velocity, shear rates, heat flux). Furthermore the transformation to collision space is linear in this case, simplifying the collision equations:\n", "\n", "$$K(f) = C^{-1}\\left( C(f) - S (C(f) - C(f^{(eq})) \\right)$$\n", "\n", "$$K(f) = f - \\underbrace{ C^{-1}SC}_{A}(f - f^{(eq)})$$\n", "\n", "in *lbmpy* the following formulation is used, since it is more natural to define the equilibrium in moment-space instead of physical space:\n", "\n", "$$K(f) = f - C^{-1}S(Cf - c^{(eq)})$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Use a pre-defined method\n", "\n", "Lets create a moment-based method in *lbmpy* and see how the moment transformation $C$ and the relaxation rates that comprise the diagonal matrix $S$ can be defined. We start with a function that creates a basic MRT model.\n", "Don't use this for real simulations, there orthogonalized MRT methods should be used, as discussed in the next section." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", "
MomentEq. Value Relaxation Rate
$1$$\\rho$$\\omega_{0}$
$x$$u_{0}$$\\omega_{1}$
$y$$u_{1}$$\\omega_{2}$
$x^{2}$$\\frac{\\rho}{3} + u_{0}^{2}$$\\omega_{3}$
$y^{2}$$\\frac{\\rho}{3} + u_{1}^{2}$$\\omega_{4}$
$x y$$u_{0} u_{1}$$\\omega_{5}$
$x^{2} y$$\\frac{u_{1}}{3}$$\\omega_{6}$
$x y^{2}$$\\frac{u_{0}}{3}$$\\omega_{7}$
$x^{2} y^{2}$$\\frac{\\rho}{9} + \\frac{u_{0}^{2}}{3} + \\frac{u_{1}^{2}}{3}$$\\omega_{8}$
\n", " " ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from lbmpy.creationfunctions import create_lb_method\n", "method = create_lb_method(stencil='D2Q9', method='mrt_raw')\n", "# check also method='srt', 'trt', 'mrt'\n", "method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first column labeled \"Moment\" defines the collision space and thus the transformation matrix $C$.\n", "The remaining columns specify the equilibrium vector in moment space $c^{(eq)}$ and the corresponding relaxation rate.\n", "\n", "Each row of the \"Moment\" column defines one row of $C$. In the next cells this matrix and the discrete velocity set (stencil) of our method are shown. Check for example the second last row of the table $x^2 y$: In the corresponding second last row of the moment matrix $C$ where each column stands for a lattice velocity (for ordering visualized stencil below) and each entry is the expression $x^2 y$ where $x$ and $y$ are the components of the lattice velocity.\n", "\n", "In general the transformation matrix $C_{iq}$ is defined as;\n", "\n", "$$c_i = C_{iq} f_q = \\sum_q m_i(c_q)$$\n", "\n", "where $m_i(c_q)$ is the $i$'th moment polynomial where $x$ and $y$ are substituted with the components of the $q$'th lattice velocity" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADhCAYAAAD/Ec//AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAb1ElEQVR4Ae2d0W3cyJaGpYWeDeMO4ACkDKzrCMbKwHIGXmcwrz1vA28GtiMYjDPwTASGlcHcAAysYDiC/f+erkWLYrOqu1msU9JXAF1kFcnzn6/aR+wiefp0tVo9Pzk5+aplrHz69ddfr8c6aIMABCAAgf0IKJ7+rSPOx45S3+nZVsf/aN07b5f/bG+wDgEIQAACRxF4N3L0ldpeuX07IL9XhCYAj9CiCQIQgMAcBBRjPwzPozY33QvIw/1Gt3Wwpzj+0HKp9e+jO1VqbGk751Jkbb1rj8wWbblP13h/ZG7jiu+21tK/fYV81+LWlow/1eZHLbda/q1ldA5E7bOXlrZzzkTW1rv2yGzRlvt0jfdH5jau+G7rEvpLA7KvhNc39yTqF637KnmRInvNbOccjKytd+2R2aIt9+ka74/MbVzx3dYl9P/XXZNsQQACEIBAKwIE5FbksQsBCEBgQICAPADCJgQgAIFWBAjIrchjFwIQgMCAAAF5AIRNCEAAAq0IEJBbkccuBCAAgQEBAvIACJsQgAAEWhEgILcij10IQAACAwIE5AEQNiEAAQi0InBIQP5pI/ZfDUS3tJ1zN7K23rVHZou23KdrvD8yt3HFd1ur6D/dyod8oVcDd2Z7U58TCrm81OLcFjdavP9n9d3LYKT22UpL2zknImvrXXtktmjLfbrG+yNzG1d8t7WGfp3zv2XF2TZPiwPyXVlsQQACEIDAHAS2A/IhUxZzaOAcEIAABCAwIEBAHgBhEwIQgEArAgTkVuSxCwEIQGBAgIA8AMImBCAAgVYECMityGMXAhCAwIAAAXkAhE0IQAACrQgQkFuRxy4EIACBAQEC8gAImxCAAARaESAgtyKPXQhAAAIDAgTkARA2IQABCLQicLaPYb3i926z//+qvtDyTm0781/sc+7cvi1t96wtp72kX+yfaz/nMrnU+veSY+bcp7X9Y3xpqb2l7WOY+djo2mvpKw7IEvBVnH5T/WkDzAmGvmr7SkvVoNzStn2dKpG1TenO9ckvj+9HLbda/q3lXMtipbX9Yxxtqb2l7WOY+djo2pfQVxSQJcTZiJ6qXgfjDbzvm+332r5yW43S0nbOn8jactpz/fLNV8LX3k/rv6jyVfJipbX9Yxxtqb2l7WOY+djo2pfQVzqH7P+YTrc5LF/U8FJCfTVVq7S0nfMpsracdvohAIFgBEoDsnMg+6vrsKSpCvfXKi1t53yKrC2nnX4IQCAYgWxALrz6rfLrIS1t58YpsracdvohAIGYBLIBWbJTsPWc4q5Sa8qipe1dvqb2yNqSRmoIQKAjAiUBucSd9PtSJfvOvU9L2zlfImvLaacfAhBYmMCZ7D3Z2Ez1UMLY3HHaJ10l+rnkGqWl7Zw/kbWdbKZU/pIT+3x7udZxNznHH3o/7A4bYbgdxk1HPUtHOiBPFkH2423eZ+w/dmpLN/cmz7VvZ0vbOa2RtVm79am6zPlB/30CsLvPpKQFbiWUpvfxlMWPzS6pHjviTzWej3SkK2T31yotbed8iqwtp51+CEAgBoFvSUbpHLJfm/XbWsPiK7CbzV/GYd9c2y1t53yIrC2nnX4IQCAYgaKArID7QbpvVb9K+rXu6YrXWt6kthp1S9s5fyJry2nfsz/dnEzfiPY8/OjdW9s/xoGW2lvaPoaZj42uvYq+09Vq9VzOO0/FhQLMzrngTQB2ciHPTfom3gstzm1R/SZQS9vycbJE1jYpvKBTvvkbgItfgPEfYI+1PyOf1ec/0lVLa/vHONdSe0vbxzDzsdG119Cnczo1xXvVp8UB+VjQHA8BCEAAAvcJbAfkoimL+6egBQIQgAAE5iZAQJ6bKOeDAAQgcCABAvKB4DgMAhCAwNwECMhzE+V8EIAABA4kQEA+EByHQQACEJibAAF5bqKcDwIQgMCBBByQ/VzpWy1TyXIOPD2HQQACEIBAhoBTMDgGnzggn2vx7+K1egtLpikQgAAEHi0Bv3jlGLwOyI+WAo5DAAIQiESAOeRIo4EWCEDgURPI5kPepqNX/JzLwsW5LC60vFPbzvwX3nHOIlvOu+H8Cpdad06NbkpU7dLVdExzA9iBvrCfyaifudyYu/+xai8OyALkBEROJvRpA8zJZr5q+0pLtaCsc9vORy2+6egUoJ7z7qJE1y59Tca0dPCi6os8rpG15cYd7ScnRQFZoJyN6KnqdTA2WK37l0S87cnoK7fVKLaj81773Fr/RZWvSLookbVLW7MxLRm8yPqCjyv/X0o+YDPvM9dnonQO2QFxLM3mF7W/lBhfxVL6IhB9TKPr62u0UdsFgdKA7Mcyxp5TTlMV7qf0RSD6mEbX19doo7YLAtmAXHj1yzPMXQz3PyKjj2l0fR0NNVI7I5ANyPInBVvPTe0qTFnsIhOzPfqYRtcXc1RR1T2BkoBc4mT6famSfdmnDwLRxzS6vj5GGZWhCJxJzZONolQPBY7NHad90pWMn0t+cGXz1fkvObbPN4BrHXcTHEb0MY2uL/jwtpHX8/+XxtqfpRFzQP6x2Uh16lvXEurH27w+FpRSW7q5tz7mofxj3+XL5UPxJ/kRfUyj60scqe8S6Pn/S2Pt3xLJ0ikLZyM6Twdt1ekK2f2UvghEH9Po+voabdR2QaA0IPt1Zb8lNyy+erzZ/HUZ9rEdm0D0MY2uL/booq5LAkUBWQH3g7y7Vf0qeal1T1e81vImtS1Qpxs56cp8AZOzmQilPdCYjgKOrm9LdKhx3dLl1cjaBlLvbT5K7aer1eq5UDinwYX+E+ycC94EYCei8byqb+K90OLcFtVvYMmGr5Zc/LKA/xDYprV+Vp//WIQtkbVLm1k2GdOSAYusL/i48v+l5AM28z6HfiZ0nNMYvFd9WhyQZ9bO6SAAAQhAQAS2A3LRlAXUIAABCECgPgECcn3GWIAABCBQRICAXISJnSAAAQjUJ0BArs8YCxCAAASKCBCQizCxEwQgAIH6BAjI9RljAQIQgEARAQdkP8/7VstUQpeik7ETBCAAAQjsTcBpAhyDTxyQz7X4d/F6fPtNsikQgAAEuibgF94cg9cBuWtPEA8BCEDgoRBgDvmhjCR+QAAC3RNwPuTiolf8nPfAxbksLrS8U9vO/Bfeca7S0nbOh8jactrdL/3OZ+L8B5dad66SUCWyPrQd9lGJzK3Eo1r6iwOyBDgBkZMJfbJg1U5M81X1lZaqQVnnb2bbvk6VyNoyuj1+H7X4Zq5Tq/peQpgirmH1oe2wj0lkbiUeLaG/KCBLiLMRPVW9DsYWr3X/koi3PRl95bYaRTaa2c75E1lbgXZfCV97P/nxiypfJYcp0hRWH9oO+5hE5lbi0RL6S+eQ/R93LM3mF7W/lFBfzdQqLW3nfIqsLaedfghAIBiB0oDsxzLGnlNOUxXur1Va2s75FFlbTjv9EIBAMALZgFx49VvlGeaWtnPjFFlbTjv9EIBATALZgCzZKdh6Tm9XqTVl0dL2Ll9Te2RtSSM1BCDQEYGSgFziTvr9q5J9596npe2cL5G15bTTDwEILEzAAfnJxmaqhxLG5o7TPukq0c8l1ygtbef8iawtp51+CEAgDoFnSYoD8o/NRqpT37rWXGmaqhiblkht6ebenWOP3WhpO6c9sracdvohAIFQBL4lNaVTFs5GdJ4O2qrTFbL7a5WWtnM+RdaW004/BCAQjEBpQPZrtX6ba1gu1XCzdbU47J9ju6XtnP7I2nLa6YcABIIRKArICrgfpPtW9aukX+uernit5U1qq1G3tJ3zJ7K2nPZBf7r5mL7xDLqbb0bWh7bDPh6RuZV4VEX/6Wq1ei7rzhVxoQCzcy54E4CdXMhzyr6J90KLc1vcqK5aWtrOORZZW4F2X+G7+AUX/4H1WPoz8Fl++Y9w0yINYfWh7bCPRmRuJR7V0K9zOj3Ee9WnxQG5RCz7QAACEIDAfgS2A3LRlMV+p2dvCEAAAhA4hAAB+RBqHAMBCECgAgECcgWonBICEIDAIQQIyIdQ4xgIQAACFQgQkCtA5ZQQgAAEDiFAQD6EGsdAAAIQqEDAAdnPnb7VMpUsp4JpTgkBCEAAAiLgFAyOwScOyOda/Lt4Ud/SkjQKBCAAgQdLwC9mOQavA/KD9RLHIAABCPREgDnknkYLrRCAwIMmcLaPd3rFz7ksXJzL4kLLO7XtzH/hHecssuW8G85vcKl159QIUyJrm4Ik3U3HdEqb+zrQx2cyN4gj/b3+f0mu1NJfHJAlwAmInEzok0WpdjKar6qvtFQLyhs7H2XLNx2dAtRz3iFKZG0lgKS/yZiWaPM+UfVFHne0lX669t9vCbZFAVlCnI3oqep1MLYrWv++2fZk9JXbahTb0XmvfW6t/6LKVyQhSmRtOUDS3mxMc9rcH1lf5HFHW8mn67B9lmBbOofsgDiWZvOL2l9KqK+WKX0RiD6m0fX1Ndqo7YJAaUD2YxljzymnqQr3U/oiEH1Mo+vra7RR2wWBbEAuvPrlGeYuhvsfkdHHNLq+joYaqZ0RyAZk+ZOC7dRTDUxZ9DXw0cc0ur6+Rhu13RAoCcglzqTflyrZl336IBB9TKPr62OUURmKgAPyk42iVA8Fjs0dp33SlYyfS6b0QyD6mEbX189Io7QHAs+SSAfkH5uNVKe+da35vDRVMTYtkdrSzb07x7IRk0D0MY2uL+aooqpjAt+S9tIpC2cjOk8HbdXpCtn9lL4IRB/T6Pr6Gm3UdkGgNCD7dWW/JTcsl2q42bqiGfazHZdA9DGNri/uyKKsWwJFAVkB94M8vFX9KnmqdU9XvNbyJrUtUKcbOenKfAGTxSYia7vnRKAxvafNDdH1bYmOPO5o2xqomVersD1drVbPJdQ5DS70n2DnXLD6HICdiOa7Ft/Ee6HFuS1uVFctsuGrJRe/LGAdtmmtn9XnPxbNSmRtOSjS3mxMc9rcH1lf5HFHW8mn67B9arDVOZ3G4L3q0+KAfJh8joIABCAAgSkC2wG5aMpi6mT0QQACEIDAPAQIyPNw5CwQgAAEjiZAQD4aISeAAAQgMA8BAvI8HDkLBCAAgaMJEJCPRsgJIAABCMxDgIA8D0fOAgEIQOBoAg7Ifp73rZaphC5HG+IEEIAABCAwSsBpAhyDTxyQz7X4d/Eivv0mWRQIQAACD5qAX3hzDF4H5AftKc5BAAIQ6IUAc8i9jBQ6IQCBB0/gbB8P9Yqfc1m4OJfFhZZ3atuZ/8I7zlVa2s75EFlb79p7Zmv20u9cMc7Fcql154FZrLS0nXMysrac9lz/Mb4VB2QZcQIiJxP6ZEGqnZjmq+orLVWDss7fzLZ9nSqRtU3pdl907dH17eIr3f6/8VGLb5Q7ba3v0yxSWtrOORhZW057rn8u34qmLGTM2Yieql4HY4vTuv/ae3s9Ge22GqWl7Zw/kbX1rr1ztt+l/1qL75z/nhuLOftls5ntnB+RteW05/rn8q0oIEvMtZaxNJtf1P5SYnxFUKu0tJ3zKbK23rX3zDbHnn4IjBIoDch+LGPsOeU0VeH+WqWl7ZxPkbX1rr1ntjn29ENglEA2IBde/VZ5hrml7VFaW42RtW3JHF2Nrj26vlGoNEJgBgLZgCwbKdh6znhXqTVl0dL2Ll9Te2RtSeOuOrr26Pp2caUdAkcRKAnIJQbS70uV7Dv3Pi1t53yJrK137T2zzbGn/5ESOJPfTza+p3qIYmzuOO2TrmT8XHKN0tJ2zp/I2nrX3pTtZsrkL0Hc55ufn6q4yYGnPyaBxmP+LFFxQP6x2Uh16lvXEurHaLw+9uFMbenm3vqYuf5paTvnQ2RtvWtvzdb2xfAyx5H+h0Og8Zh/SyRLpyycjeg8HbRVpytk99cqLW3nfIqsrXftPbPNsacfAqMESgOyX/30G0fD4quIm81fl2HfXNstbed8iKytd+09s82xpx8CowSKArIC7gcdfav6VTqL1j1d8VrLm9RWo25pO+dPZG29a++Z7YB9uvmYvk0OuqtutrSdcyyytpz2XP/BvnkOubRcakcnE3qh2jfxXP+s7SVuZLS0LTcnS2Rtk8LVGV17dH07+er/ha/wXV7+U538oTbfa/ms2hc41UpL2zmnImvLac/1z+Hb6Wq1ei5DTt5zoRNWuTmXc4R+CEAAAo+VgOKucwW9V31aNGXxWEHhNwQgAIElCRCQl6SNLQhAAAITBAjIE3DoggAEILAkAQLykrSxBQEIQGCCAAF5Ag5dEIAABJYk4ID8ZGMw1UvaxxYEIACBx07g/3NZOCCnHBapfuxw8B8CEIDAkgT2zmWxpDhsQQACEHiUBJhDfpTDjtMQgEBEAgTkiKOCJghA4FES2CeXxYle7Xu3oeRcFhdanNtikdetW9rOfTIia+tde3S20ufUA85bcal151EOUyJry0F6rNqLA7IAOd/Fb6o/GaZqZ3v7qvpKS9WgrPM3s21fp0pkbVO63Rdde1R90uXP/kctt1qclvZcS4gSWVsOENpPToqmLATKyS+eql4HY4PVuq8GvP3e27VKS9s5nyJr6117ZLb+7GvxTza9Feffc6yX7I+sLccB7YUBWSCvtYyl2fyi9pcC6SuGWqWl7ZxPkbX1rr1ntjn29ENglEDRFbKOdE5Xfz0bljRVkXK+Dvvn2G5pO6c/srbetffMNseefgiMEsgG5MKr33+Nnv3Ixpa2c9Ija+tde89sc+zph8AUgWxA1sEp2HrOeFepNWXR0vYuX1N7ZG1J4646uvbo+nZxpR0CRxEoCcglBn4q2anSPi1t51yKrK137T2zzbGn/5ESOJPfngf23eKxOWJj2dXuvnQl4+eSa5SWtnP+RNbWu/ae2ebYP9j+zVTTX3Jwn2/MflrlpjWUxtr/lP+OwScOyOda/OiaG+9NS0ioH/FR1yjkBD7d3PN+s5WWtnNORNbWu/ae2ebYP+R+j5v8u+zRx8bafQPbMfhD6ZSFg7UD97CkK2T31yotbed8iqytd+09s82xpx8CowRKA7JfDfUbScPiv4Y3m78uw765tlvazvkQWVvv2ntmm2NPPwRGCRQFZAXcDzr6VvWrdBate7ritZY3qa1G3dJ2zp/I2nrX3hHbdHMxfVvMoV+yP7K2HIdHqd1zyKXFV8NOJvRCtW/iuf5Z20tMyLe0LTcnS2Rtk8LVGV17WH363PsK3sXzfy5/qM33Uj6r9gVMsxJZWw7KY9d+ulqtnguSk/dcbD5QOWb0QwACEIDATAQUd50r6L3q06Ipi5nschoIQAACEJggQECegEMXBCAAgSUJEJCXpI0tCEAAAhMECMgTcOiCAAQgsCQBAvKStLEFAQhAYIKAA/KTTX+qJ3anCwIQgAAEZibwLJ3PAfnHZiPVqY8aAhCAAATqE/iWTDBlkUhQQwACEGhMgIDceAAwDwEIQCARICAnEtQQgAAEGhPYJ5fFiV7te7fR61wWF1qc26JKLuQhl5a2h1qG25G1DbUOt6Nrj65vyHO4Lf1OTeC8F5dav5dvfLj/nNstbR/rR3TttfQVB2QJcL6L31R/MmzVzvb2VfWVlqpBWedvZtu+TpXI2qZ0uy+69uj6dvGVbv/f+KjlVovT1p5rWaS0tH2sg9G1L6GvaMpCQpz84qnqdTA2eK37r723nem+WmlpO+dUZG29a++crX9lxz9N5J/l+T03FnP2y2Yz28f6EV37EvqKArJAX2sZS7P5Re0vJdRXBLVKS9s5nyJr6117z2xz7OmHwCiB0oDsnK/++jUsaaoi5YQd9s+x3dJ2Tn9kbb1r75ltjj39EBglkA3IhVe//xo9+5GNLW3npEfW1rv2ntnm2NMPgSkC2YCsg1Ow9ZzxrlJryqKl7V2+pvbI2pLGXXV07dH17eJKOwSOIlASkEsM/FSyU6V9WtrOuRRZW+/ae2abY0//IyVwpq+Hvll3OuH/2Nxx2j1dyfi55Bqlpe2cP5G19a69KdvNlMlfgrjPNz8/VXGTA/+Q++F22OiKm3+Dcf07jA7Iz7Xh53xHf1NP/X6MRt2jH870gU0397zfbKWl7ZwTkbX1rr01W9sXw8scR/rvEoDbXR6lW+K292/q/amTn48YSFfI7q9VWtrO+RRZW+/ae2abY08/BEYJlM4h+9VPv3E0LL6KuNn8ZRz2zbXd0nbOh8jaetfeM9sce/ohMEqgKCAr4Hp+41b1q3QWrXu64rWWN6mtRt3Sds6fyNp6194z2wH7dPMxfZscdFfdbGn7WMeia6+i72wPapfa18mEXqj2TTzXP2t7iRsZLW3LzckSWdukcHVG1x5d306++n/hK3wXv+Di8ofafK/ls+r1DZx1a4V/Wto+1p3o2mvrO12tVpM39Y4FzPEQgAAEILCbgIL83jf1dp+NHghAAAIQmIVA0RzyLJY4CQQgAAEITBIgIE/ioRMCEIDAcgQIyMuxxhIEIACBSQIE5Ek8dEIAAhBYjgABeTnWWIIABCAwSYCAPImHTghAAALLESAgL8caSxCAAAQmCRCQJ/HQCQEIQGA5AgTk5VhjCQIQgMAkgX1yWZzoFb93m7M5l8WFFue2qJILeai6pe2hluF2ZG1DrcPt6No70OfUA85bcSmtzqMcpkgP2iqNRi22xQFZApzE/jfVn+yjamd7+6r6SkvVoKzzN7NtX6dKZG1Tut0XXXtUfdLlz/5HLf5lE6elPdcSoqCt3jAswbZoykJCnPziqep1MLbLWvfVgLffe7tWaWk751Nkbb1rj8zWn30t/smmt+L8e471kv1oq0d7CbZFAVkuXmsZS7P5Re0vJdRXDLVKS9s5nyJr6117z2xz7OmHwCiB0oDsnK7+ejYsaaoi5Xwd9s+x3dJ2Tn9kbb1r75ltjj39EBglkA3IhVe/VX4NoaXtUVpbjZG1bckcXY2uPbq+Uag0QmAGAtmALBsp2HrOeFepNWXR0vYuX1N7ZG1J4646uvbo+nZxpR0CRxEoCcglBn4q2anSPi1t51yKrK137T2zzbGn/5ESKAnIY3PHCVe6kvFzyTVKS9s5fyJr6117z2xz7OmHwE4C2YCs+bw0VTE2LZHa0s29nYYO6WhpO6c3srbetffMNseefghMEcgG5M3Bf6o+HzlRukJ2f63S0nbOp8jaetfeM9sce/ohMEqgNCD71VC/kTQsl2q42bqiGfbPsd3Sdk5/ZG29a++ZbY49/RAYJVAUkBVwP+joW9Wv0lm07umK11repLYadUvbOX8ia+tde0ds083F9G0xh37JfrTVo12F7dkeen017GRCL1T7Jp7rn7V9o7p2aWk751tkbb1rD8tWn3tfwbv4BRaXP9TmeymfVfsCpllBWz30tdmerlar55Lv5D0Xmw9UPW84MwQgAAEI3CGguOtcQe9VnxZNWdw5mg0IQAACEKhCgIBcBSsnhQAEILA/AQLy/sw4AgIQgEAVAgTkKlg5KQQgAIH9CRCQ92fGERCAAASqENh+7O1v3eUbGvmkNicKp0AAAhCAwJEEFE//1inOd53GAdnPTvqnaMZKlRwVY4ZogwAEIPAICKQfih519f8AJvvxikUCJ4QAAAAASUVORK5CYII=\n", "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\\\0 & 0 & 0 & -1 & 1 & -1 & 1 & -1 & 1\\\\0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1\\\\0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1\\\\0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1\\\\0 & 0 & 0 & 0 & 0 & -1 & 1 & 1 & -1\\\\0 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1\\\\0 & 0 & 0 & 0 & 0 & -1 & 1 & -1 & 1\\\\0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\\end{matrix}\\right]$" ], "text/plain": [ "⎡1 1 1 1 1 1 1 1 1 ⎤\n", "⎢ ⎥\n", "⎢0 0 0 -1 1 -1 1 -1 1 ⎥\n", "⎢ ⎥\n", "⎢0 1 -1 0 0 1 1 -1 -1⎥\n", "⎢ ⎥\n", "⎢0 0 0 1 1 1 1 1 1 ⎥\n", "⎢ ⎥\n", "⎢0 1 1 0 0 1 1 1 1 ⎥\n", "⎢ ⎥\n", "⎢0 0 0 0 0 -1 1 1 -1⎥\n", "⎢ ⎥\n", "⎢0 0 0 0 0 1 1 -1 -1⎥\n", "⎢ ⎥\n", "⎢0 0 0 0 0 -1 1 -1 1 ⎥\n", "⎢ ⎥\n", "⎣0 0 0 0 0 1 1 1 1 ⎦" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Transformation matrix C\n", "method.moment_matrix" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAucAAAAVCAYAAADlwkDzAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAIzUlEQVR4Ae2d63EcNwyAZY8LUJIOLh1IUQW2OpDcQewO4p/Sv4zSgeMKMnYHiSvIRB04HdhRBxdgb1fH45EEH+ADWuzMekkCJIAP5Iq3XulOttvtiXne3NxszLpdpuS2fuv6CP6N4IOP+wi+jeCDjw9H+wjxjeADB0vfGCPEN4IPPj7YPoJ/I/gQYlQiGyG2EXwIMRzBvxF8CDHqLRuBzwg+1MwDFZ9L/vzEOG5vb3+B6pnRdFCk5AfK/Sqb2c8uHghgpHzqzwxlrIzrE6At6DykGZVoKF+anjKiGfXW0BzVzwDF+Ej+DD8t4AGbyiu4XMD13dRg/ROSg+xuVv8K1x/hvIO2f60hWKoxtkAHP2Q8wPV3FqORg4A9L8MYvyPNkGpgCz9gfYTzHMoPdgdoWzUfm0dOXSLjnDjtPlTctn5uHeyQ95QR5zHGq4zonye582Lp15CxuHtpyzkYY2vUdbrMpdxrqzmY65/Zj/J11BxRfpsx9i5TvlKMbfm0OYfGUwjsM1zPXQGG5CD7B/r8CtdP2HfWxbZLKLNu0GG8aFuz7ku4Hm1O0U/uA+x4Gab4nevXbP8D9P8G509w4g+V73zxzz6thg+wKD4kMy4JPjXuElvYd56bUfeUWbf7PFZG+6zPLLw/T/aaaaVWjFPtgD7+XFrVHJTKKG3GHWunxn08QruWVF/XOI9Ls8HN2MzBi9k5fEr1PuCoUw4DvYE+p3CdNubYH8r4xBrrON4ltnEcMGaqLbSPfr/lsB8xRldGyB18vEY/oRx8PWmOZVV85piLLlIZFwUNnTPizjYJtqSu89T1tzpG2QHPHVvNwww7o9xLW87BVFtDMJIyB0v9xP46jzkohseowPhxnSzvnL8GI6FXQHxy3AzeO9z/G9pewZinDlluU5KtOR70m9OHkO8jMAr5dyBTPgc4qlQGYlwlvkqDSl3nlXA4h1VGTix9GnWd09yVEc2ot4bmqH4GKMam/DlU8D1p7+snhPwV9MXXKOxjGQ/lXEeOLfTjNZcDvnEGYuRz0deufHxk+NpHYMwXTf2RpK7z+mT2FpTRnsUoJV3ndCaUEc2ot4bmqH4GKMaTHJ+c46snfwX8ccphQxrzRPr7wLjRogJbGBfbqzUBh7szCvgWEimfEB0eWVfGPCG0GUXqOm9DZ2dFGbWknWRL1zmNSxnRjHpraI7qZ4BiPMlxc46/PPgl4I9Pvmy88d0z3xGzgff1NdtzbWFcG3OgSuURGOWEpnxyqKX16c04zdu+2lLXeUtqyqgl7Xhbus5pVsqIZtRbQ3NUPwMU40mOm3PcQLteTVlcpOSLnuv6g6uxUpvLFsbVYnMuhZGNXvnYRPjrEhjzR11vRKnrvB6R45GV0TGT2i26zmnCyohm1FtDc1Q/AxTjSf4C/MAnMQ8Bf3xyHMB3LE93vvoUEttzbeG7O7hxrn2MwCgnxlXxmV8J+AygUubENfS7z4E79+nN+KRT3DnIpK7znFhz+4hlJGge5uSm+zrPcbpxn+6MJM3BTr52z1HLOTko4ykHuDnPOiAo/JOJ2Ne10Vna0EjxUWDLt2ku9ilmgAK/Y4bn0FkVH8wHQDvnAJcwRlfG6GenuBMQ7VQL1kt3xsnBZnaQzEjKPMxMzWrmYCYf7NadkaQ52MnX7jkqmF/JXQdlPOUAX2vBJzHLZtoVXEiOL667XhvBwfFAOdeRYwvjCj1p4vJtFEap8SifVGLp+iMwTve6Xw+p67wlMWXUknacLV3nNCdlRDPqraE5qp8BivEkx805Pt12bbAXF0Ny/Ip4/GVI+8Cnk/fzp5JJBmU0WHJE2zKM4IcE8uk9g2+jMDJCjyoqnyhMRUojMC4KIKUzw1qSus6jMY3KiMEvZBC6F0YzKlVkiiXFDVHrvAMfZCmKUUryXbqdGLtcSWkTlaMnynjKAW7O7+G8CGTPKwcw+MVF3+B6tfSfYeHfFv/ZavsPZPg1x1lHrC1r8OlDgtV2UJ39LfINBhyCkRHY8gthmOTQsVY+ISaxMjGMYwOK1PPGzbGWpK5zi504Rhy5mxl474UWo9JqVcaGc147hg4Wz+HE2L0HE+NYvl6/mfww4/TaMpWgPBojy73kqjfuCoyTnbM6eH219EbLkdfvJ8x4ygFuzv+A88xKkFml5DjQJYC6gxO/Nv4DnC+h/HijgvIDtOETlVMolzxBJ22BDfPAL+v402ywy0y+DcEIYvmIJ8SIX4GOx1SHtqW+a93/uyo++7DzS4IZ5wcNPWPiBp1Vr3PJjBhzR90LRczDmFxagQxxL43xmyvXMbZGZGT5lFyNiZuLcbJzVocYX60uq5vHVvzJ1QqMdznYbrcnNzc3X+A8w7LrpOSuPq42GOcKzlOXjLsN7GzQ79hxS31DW3B6GbbyI8GO8vHM91iGlB7MB5GMqbgoeelaosY35VIZj8qIwy8Yg+VeaOY5p8wRS4xdsCNynbfigwylMorJf0inJeOQHzEyqTl6SozNHOCTczzu4Hw7ldz/UHJ3r+PWi/kT5bGEv+UdDIl+xx6lvkljpHxiZ0a+nlTG+RHvepaupRT7UhmPyojDL657Yco8cOlyxOIa127TOWgTOa5LZXQcSVpLqzmY5pVbW2qOnhLjxxxMm3PYMOO74xu4blw5o+SuPnYbjIGvs3D93XN7+IP6HAfGg3GRB4dvsy0vQ9IJUODwI9IO5ln5xMDK1IFcimScGe5jt1ZzGA1KZTwqIy6/YJzgz5PHyVKxwBUL5SLYEbnOW/FBflIZUbmn5C0ZU75Qcqk5ekqM7Rw8N5J2DeX3Rt0uUnJb366/AeO/2Y2V6hhH6H8CbLNcvklhpHzsGcBfl8q4lATXWorxQyrjURlx+lV6L4zJf0iHM5aQHZ2DITo7mVRGdGRhjVZzMOxFnFRqjp4S44McPMN3kZZj3rlf+TbRlHwZp+cVfMRfSv0EV/JPKNbwc3RGyqdG1g/HVMaHPGrUlDFNVRnRjEo0lC9NTxnRjHpraI7qZ4Bi7JL/D0hX4JBsdFfbAAAAAElFTkSuQmCC\n", "text/latex": [ "$\\displaystyle \\left( \\left( 0, \\ 0\\right), \\ \\left( 0, \\ 1\\right), \\ \\left( 0, \\ -1\\right), \\ \\left( -1, \\ 0\\right), \\ \\left( 1, \\ 0\\right), \\ \\left( -1, \\ 1\\right), \\ \\left( 1, \\ 1\\right), \\ \\left( -1, \\ -1\\right), \\ \\left( 1, \\ -1\\right)\\right)$" ], "text/plain": [ "((0, 0), (0, 1), (0, -1), (-1, 0), (1, 0), (-1, 1), (1, 1), (-1, -1), (1, -1))" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAF5CAYAAABtIcr0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAg5ElEQVR4nO3deXRlVYGo8W8nNc9AiqEoQApKVJRBWkGQQRww4mucQG2ccVZEQV83jd363lNpVKAFBBHHtsHXOLc2QW1BQUSeouKMIMhoUaSoEWpKst8fN7cquZWb3PGM328tlubmnGTrOvk4d599zwkxRiRJ+daT9gAkSe0z5pJUANNqXwgDgwF4DnAssBAIdfaNwDrgBuD7sb9vpFuDlBoRBganA/3AkcBcPHaVE2FgcBZwInA4MJvJj901wPWxv++6sd/YIebAPwGnNjGOlwFXj+4npeli4FlNbP8y4CvA+7szHGlqYWCwF7gCeHoTu50SBgY/F/v7zqu+MG6aJQwMLgZe2cJ4Tg4Dg3u0sJ/UEWFg8CCaC3nVyWFgcEmnxyM14Rk0F/Kq14aBwUXVL2rnzJ8ywWuNCMDBLewndcohbezrsas0HdTifr3Ak6tf1IZ7ZsvDaW9fqV0z2tjXY1dpauf423bcTzRnvqMzn78Xd/5mFr29la932nWIz/387jYGIHXfVy5exPVfWcj9d87gGS9Yz9mfWZH2kKSmfO+q+fzHhbuwasV0Fu4yxBn/uoKnHrdxok0biznAaR9YyUlvXtuxQUrdtsseQ5xyxipuvX4uWzbVWx0gZdNPr53Dl85dzPs+9SBPPmITDz84aa8bj7mUN8e/bAMAf/rVLFb91WNd+XLVx/o4+V2rOOioTQDsttfQZJs3frHzyo/1ccr++/Gu5+zNz38wu71RSpLqGh6Cu38/i7Wrenndofty6oHLuPBdu7LpsbrvMBuL+ev/+WE+f+td/Ptv7+KEU9fw4Tcs5b47pnds4JKk7Vat6GV4CG6+Zj4f+869XHL9Pdz9+1l88SO71NulsZg/5chNzF0QmTErcuLr1/H4Qzby04G5HRu4JGm7mbMrd0A88fWrWbznMDvtOsyL3vIIv/xh3e62eG+WEPFui5LUHQt3GWGnXYcIY2dVJr+GP3XM1z3Sw83XzGHzxsDQVrj2S/O5/dY5PP15j7Y5XKm7hrbC5o2BkWEYGWbbMSzlwbNeupbvfH4Rq1b0snZVD//56Z047PgN9Taf+gr/0NbAl87r46Nvm0lPT2SPfbfwD1c8wOOe6F+Fsu2LH96Fr31y+xzjTd9ZwEvfsYrTPrgqxVFJjXnt+1exbnUvb37GvkyfETmifz2vOfuReptPHfOddxvm0h/d29FBSkk47YOGW/k1fQacefFKzrx4ZSObez9zSSqA2pi3c1XTK6LKK49dpakjx19tzCf8zH+DHmtnIFKb2rkg77GrNHWku7Ux/wXQyoXNYeDnbQxIatctLe43DPyskwORmtTqsbsRuK36xbiYx/6+9cDHW/ihF8T+vjUtDkhqW+zvuwv4Ygu7XuixqzTF/r5fAV9rdjfgvNjft+2sPsQJPvwTBgaX0cQzQGN/351NDkTqijAweCBNPAPUY1dZEQYGDwaOYOpngK4Gfhj7++4Zt/9EMZck5YtLE1V4IYSvhRAuTnscUjd5Zq5CCyHsAdxN5ULnrjFGb0OhQvLMXEX3RirzjCPAKSmPReoaz8xVWCGEHmAFsHj0pd/HGA9McUhS1yR2Zh5C2DmE8N0QwmtDCLOS+r0qtecCY4+1x4UQDk5rMCqPEMK8EMI7QgjfTqp3SU6zbASOAi4HVoYQzgsh7Jng71f5vAeYP+brmcDpKY1FJRBC2C+EcCnwEHAh8FRgSxK/O7GYxxg3Ap+isn5yPnAGcOfov7mOCiH49HR1zOiFz+NqXu4FXhlC8ClZ6phQ8bwQwvXAb6lcp5lD5dP058UYR5IYR9IXQC+iciEKKmdJs4ATge8CtzsFow6qXvis5YVQdUR1KgW4l8onOI+j0rTq85ED8IXExpP0BdAQwgBwAhN/wmkDlT/Ay4CLYowPJDk2FcMEFz5reSFULQsh7AecBbyWSq8meqc3BHwuxviWpMaVxtLEc6l/h7t5jJ+C+Y5TMGpB7YXPWl4IVVMmmUqpN2W3FTg/qfFBOmfmAbgD2K+BzSOVWzw+ADwzxvhwN8emYgghXEvl3V89w8AXYoxvTGhIyrEQwv7A9cAiKiecjfhxjPHorg1qAomfmcfKvz0+QmVKZSqByk1ndsYPOKkBdS581vJCqJoRqNx0sNGQb6DSuESlFcgvN7jdCLAKODzG+FAXx6PiqHfhs5YXQtWQGOMdwNFU7rTZiHVUFnUkKpWYjy5TvJzJ119WQ35EjPGuRAamXBu98Hk6k8+XV80D3tvdEakoYoy3AccwddAfI8HliGOl9nH+EMLewO1M/Ic3AqwBnmbI1agQwvOA/wKmNbjLCHBojPHX3RuVimT0wvlPqX/CsBHYPcbY6Fl8xzR60HdcjPHe0SvDz2f8MsURKu8YdgY2pzE25dYfgKsneP3vRv/zqprXh6msEZYaNUT9kA8BX0oj5JDyjbZCCMdQOZOqXlioTq0cQ+UPE2Cp683VjhBCBB6KMe6e9liUXyGEA6ksS4TKx/R/CCwYs8km4OAY458SHhqQ/gqRG6ncwwDGz5H/ke3/9rvfe7hISlNNyHtjjL9kxzn0n6cVckg55mOWKW6l5mJnjHEzBl1SyiYI+QjscFF0MyksRxwr7TNzqCxT/DwTrFox6JLSVC/kVWOC/mlSWI44Vi4eThFCmEllPgqcQ1eTnDNXK6YKedZk4cx8Sp6hS0pS3kIOOYk5GHRJychjyCFHMQeDLqm78hpyyFnMwaBL6o48hxxyGHMw6JI6K+8hh5zGHAy6pM4oQsghxzEHgy6pPUUJOeQ85mDQJbWmSCGHAsQcDLqk5hQt5FCQmINBl9SYIoYcChRzMOiSJlfUkEPBYg4GXdLEihxyKGDMwaBLGq/oIYeCxhwMuqSKMoQcChxzMOhS2ZUl5FDwmINBl8qqTCGHEsQcDLpUNmULOZQk5mDQpbIoY8ihRDEHgy4VXVlDDiWLORh0qajKHHIoYczBoEtFU/aQQ0ljDgZdKgpDXlHamINBl/LOkG9X6piDQZfyypCPV/qYg0GX8saQ78iYjzLoUj4Y8okZ8zEMupRthrw+Y17DoEvZZMgnZ8wnYNClbDHkUzPmdRh0KRsMeWOM+SQMupQuQ944Yz4Fgy6lw5A3x5g3wKBLyTLkzTPmDTLoUjIMeWuMeRMMutRdhrx1xrxJBl3qDkPeHmPeAoMudZYhb58xb5FBlzrDkHeGMW+DQZfaY8g7x5i3yaBLrTHknWXMO8CgS80x5J1nzDvEoEuNMeTdYcw7yKBLkzPk3WPMO8ygSxMz5N1lzLvAoEvjGfLuM+ZdYtClCkOeDGPeRQZdZWfIk2PMu8ygq6wMebKMeQIMusrGkCfPmCfEoKssDHk6jHmCDLqKzpCnx5gnzKCrqAx5uox5Cgy6isaQp8+Yp8SgqygMeTYY8xQZdOWdIc8OY54yg668MuTZYswzwKArbwx59hjzjDDoygtDnk3GPEMMurLOkGeXMc8Yg66sMuTZZswzyKArawx59hnzjDLoygpDng/GPMMMutJmyPPDmGecQVdaDHm+GPMcMOhKmiHPH2OeEwZdSTHk+WTMc8Sgq9sMeX4Z85wx6OoWQ55vxjyHDLo6zZDnnzHPKYOuTjHkxWDMc8ygq12GvDiMec4ZdLXKkBeLMS8Ag65mGfLiMeYFYdDVKENeTMa8QAy6pmLIi8uYF4xBVz2GvNiMeQEZdNUy5MVnzAvKoKvKkJeDMS8wgy5DXh7GvOAMenkZ8nIx5iVg0MvHkJePMS8Jg14ehrycjHmJGPTiM+TlZcxLxqAXlyEvN2NeQga9eAy5jHlJGfTiMOQCY15qBj3/DLmqjHnJGfT8MuQay5jLoOeQIVctYy7AoOeJIddEjLm2MejZZ8hVjzHXOAY9uwy5JmPMtQODnj2GXFMx5pqQQc8OQ65GGHPVZdDTZ8jVKGOuSRn09BhyNcOYa0oGPXmGXM0y5mqIQU+OIVcrjLkaZtC7z5CrVcZcTTHo3WPI1Q5jrqYZ9M4z5GqXMVdLDHrnGHJ1gjFXywx6+wy5OsWYqy0GvXWGXJ1kzNU2g948Q65OM+bqCIPeOEOubjDm6hiDPjVDrm4x5uoog16fIVc3GXN1nEHfkSFXtxlzdYVB386QKwnGXF1j0A25kmPM1VVlDrohV5KMubqujEE35EqaMVciyhR0Q640GHMlpgxBN+RKizFXooocdEOuNBlzJa6IQTfkSpsxVyqKFHRDriww5kpNEYJuyJUVxlypynPQDbmyxJgrdXkMuiFX1hhzZUKegm7IlUXGXJmRh6AbcmWVMVemZDnohlxZZsyVOVkMuiFX1hlzZVKWgm7IlQfGXJmVhaAbcuWFMVempRl0Q648MebKvDSCbsiVN8ZcuZBk0A258siYKzeSCLohV14Zc+VKN4NuyJVn02pfCAOD04GXAM8EFgChzr4RWA/cBHw99vdt7tYgpbFijJtDCLOATVSCvjTG+EDdY/esS2DajIVhYPDfqj8CWA1cD3w79veNGHKlKQwMzqZy7B4JzKN+d0eAVcB/x/6+gbHf2CHmwEXA8U2M47mj27+piX2ktkwUdK55+Bzg2TtsfMBh0NM7HTi85jv9wFNDCFdjyJWSMDAYgE8BRzSx2wvDwOCBsb/v49UXxk2zhIHB5TQX8qpjwsDgk1rYT2rZuCmXJcvuZ2jrCU3/kC2bXsW8RYZcaTqU5kJe9ZowMDin+kXtnPmT2xjQU9rYd5wQQr23GNI424K+zxNgxT3LGNo60bvNiW3ZPIOV9y9j7wPAkKsJHW5Uq92dCSyvflEb8xktD6e9fQkVzwwhfBVYFULYt52fp/KIMW7mrE++Cmg86Fs2z2DlfY8D4NyvP8eQq1EhhMOA1SGEL4QQntqBH9mR7ja2muVDr9udVz5hP16yz/68/rB9+danF7bxy8cJIewcQng38BdggMpFgOnAzp36HSqBOfOG2HPZHcD4oK8d7OHSv4e3H9PLq5+yjO/++/xxId9zvz8xfUZMa9jKpd2BXuBVwI0hhD+GEN4UQpjfld92zx+n87d7LudDr9t9ss0ae0v6ijMfYe/HP8SMWZG7fzeDs1+yF8sP2cSTnt7SCpbRtyhHAe8GTqRyhXbOmE2GW/m5KrnQE9lz2R08cNdyVtyzjN33uYtPnLmYadPhgmuHWTv4Vz702qXstGtgybJKyJ3RU2uGqQR9DnAAcAHwidGL6RfFGH/Rsd/0yfftxrIDN021WWNn5vsftIUZsypnLyFEQoAH7mr6rUGds/BZjA+51Lpq0AHu+eMyfvb9+Zz0Zpg1Bw48YpiDjg789FpDrk6bB8xmx7P1eW391O9dNZ85C4Z5ypGPTbVp4x8auuD0XXnR0uW8/dh9WbR4iKNO3NDIbjVz4Q8AHwb2ZvK1lFLrqkF/6F7o6YXd9gYirLzvcey1P6y451FDri6pPVtfOTq3fmjTP2nDmh6+fH4fbz334UY2b/zK/5kXr+SMC1fy65tmc9uNs5k+c8p5xtGz8PdQmf+eS+Px7gGOCSEsbnh8Krf3XXYwRzx//Du83mkPM3tu5RgaGekFYGHfGjZumM3GDdu3Hfi3p4UXfKCtC/gqlSNorGXVs/JXASeHEO4Dzo8xXtHQb/ns/+rj+FPWsvveQ41s3njMAXqnwaHHbuQHVy/gG5ct4pQz1kyxx4VU5sObvW3AfCr/VpMac8t3Yfkh418bHoKNj45/be2qRcycDatWLN322s++/y/dH6BKrHq2vhy4HJg65rffOpPf/GQOl97wl0Z/SXMxrxoZgr/+pZEzmScA7wReN/p1o/NHa4FnxxhvbWF0KqEwMPhy4H+Pe3H+TjMZGd6Hh+6D3faqvPbXuzewzxO3sHT/wW3bnfuNt8X+vusSHK5yLIRwInAl0OiqvvVULpheRiXmU/vlDXMYfHA6rzloPwA2b+xhZATedvRMLrvxnol2mfqMedWKXr531XweWx8YHoKbr5nDT65ZwCHHPDrVrjHG22OMpwOLgTcDvwA2Alsb+h8ktWrL5hmsX70Phx4H/3kFbN0yzJ23wa3Xz+PZp0x57Ept2kLlVhM3AK8GFscY/zHGOGGId3DSm9fwmVvu4pLr/8Il1/+F57x8DYcc/Sgf/ur99XaZ+sw8BBj44iIuP2c34gjssscQr3v/So59ccN/EDHGTcCXgS+HEA6gtbN1qTFj15GfedGdnPfW/XnPCb3MWzjEqe+bxpz5ezG09S6mTW9oLlJqwriz8IbjXWv23MjsuduXaM+aO8L0mSPsvFvdZdtTx3zn3Ya58Lv3tTSgCcQYbwdODyG8D3gx8F7giaNjmd6p36OSqv1AUAjwjvOgp3eYJfveRRwJ49ahG3S1bwuVa4M/A84H/ivG2Nnj6rQPrppqk9TuZx5j3BRj/HKM8TAqN5q5HNgw+s+sSXeWJjJRyGuNXYfe7L1cpO1mUTkLX0NlscYTYozHxBi/1fGQN6g25u18rLnlfSeYW/8WlQ8WSY359U17Thnyqtqg/+6WvgRGqOL4DfBtts+Fn93ydEpFR7pbG/N2Lgw19CGiyYw5W395jHHKtxUSjD4h6KuXnA80/snOsUG/4p/+I4mHRKsYYoz3xhhP7uBZeDvd3bZvbcxvoTL306wI3NzGgKSWbHtC0O23wp77/aGpT3aGnsiSZb/irt9Blx8SLU3iJy3uNwjcXv1iXMxjf98g8NEWfugFsb/voRYHJLVk3KPe1q/uJYSPNfkjhunp+QBDW7r6kGhpMrG/717gkiZ32wr8c+zv23byHWLccbomDAwuYftzFCezDrgp9vc90ORApLbUe2bn6LF7FGM/0PGJ91zBtOnreMdHzxrzIx4Bbhg9gSGEMJPKumCApTFGj2klKgwM7s32Z4DWE6mckd8Q+/tWj9t/ophLWdbsw5dDCBF4KMY46f2gDbryLLWliVIrmg15M8Y9U9QpF+WMMVdudDPkVQZdeWXMlQtJhLzKoCuPjLkyL8mQVxl05Y0xV6alEfIqg648MebKrDRDXmXQlRfGXJmUhZBXGXTlgTFX5mQp5FUGXVlnzJUpWQx5lUFXlhlzZUaWQ15l0JVVxlyZkIeQVxl0ZZExV+ryFPIqg66sMeZKVR5DXmXQlSXGXKnJc8irDLqywpgrFUUIeZVBVxYYcyWuSCGvMuhKmzFXoooY8iqDrjQZcyWmyCGvMuhKizFXIsoQ8iqDrjQYc3VdmUJeZdCVNGOuripjyKsMupJkzNU1ZQ55lUFXUoy5usKQb2fQlQRjro4z5Dsy6Oo2Y66OMuT1GXR1kzFXxxjyqRl0dYsxV0cY8sYZdHWDMVfbDHnzDLo6zZirLYa8dQZdnWTM1TJD3j6Drk4x5mqJIe8cg65OMOZqmiHvPIOudhlzNcWQd49BVzuMuRpmyLvPoKtVxlwNMeTJMehqhTHXlAx58gy6mmXMNSlDnh6DrmYYc9VlyNNn0NUoY64JGfLsMOhqhDHXDgx59hh0TcWYaxxDnl0GXZMx5trGkGefQVc9xlyAIc8Tg66JGHMZ8hwy6KplzEvOkOeXQddYxrzEDHn+GXRVGfOSMuTFYdAFxryUDHnxGHQZ85Ix5MVl0MvNmJeIIS8+g15exrwkDHl5GPRyMuYlYMjLx6CXjzEvOENeXga9XIx5gRlyGfTyMOYFZchVZdDLwZgXkCFXLYNefMa8YAy56jHoxWbMC8SQayoGvbiMeUEYcjXKoBeTMS8AQ65mGfTiMeY5Z8jVKoNeLMY8xwy52mXQi8OY55QhV6cY9GIw5jlkyNVpBj3/jHnOGHJ1i0HPN2OeI4Zc3WbQ88uY54QhV1IMej4Z8xww5EqaQc8fY55xhlxpMej5YswzzJArbQY9P4x5RhlyZYVBzwdjnkGGXFlj0LPPmGeMIVdWGfRsM+YZYsiVdQY9u4x5Rhhy5YVBzyZjngGGXHlj0LPHmKfMkCuvDHq2GPMUGXLlnUHPDmOeEkOuojDo2WDMU2DIVTQGPX3GPGGGXEVl0NNlzBNkyFV0Bj09xjwhhlxlYdDTYcwTYMhVNgY9eca8ywy5ysqgJ8uYd5EhV9kZ9OQY8y4x5FKFQU+GMe8CQy6NZ9C7z5h3mCGXJmbQu8uYd5AhlyZn0LvHmHeIIZcaY9C7w5h3gCGXmmPQO8+Yt8mQS60x6J1lzNtgyKX2GPTOMeYtMuRSZxj0zjDmLTDkUmcZ9PYZ8yYZcqk7DHp7jHkTDLnUXQa9dca8QYZcSoZBb40xb4Ahl5Jl0JtnzKdgyKV0GPTmGPNJGHIpXQa9cca8DkMuZYNBb4wxn4Ahl7LFoE/NmNcw5FI2GfTJGfMxDLmUbQa9PmM+ypBL+WDQJ2bMMeRS3hj0HZU+5oZcyieDPl6pY27IpXwz6NuVNuaGXCoGg15RypgbcqlYDHoJY27IpWIqe9BLFXNDLhVbmYNempgbcqkcyhr0UsTckEvlUsagFz7mhlwqp7IFvdAxN+RSuZUp6IWNuSGXBOUJeiFjbsgljVWGoBcu5oZc0kSKHvRCxdyQS5pMkYNemJgbckmNKGrQCxFzQy6pGUUMeu5jbsgltaJoQc91zA25pHYUKei5jbkhl9QJRQl6LmNuyCV1UhGCnruYG3JJ3ZD3oOcq5oZcUjflOei5ibkhl5SEvAY9FzE35JKSlMegpx7zEMIeIYT7QggvrvN9Qy4pcY0GPYTwlhDCn0MI85Ib3Y5SjznwTmA34MoQwovGfsOQS0rTVEEPIbwFuABYArw64eGNE2KM6f3yEGYCK4EFoy9tBP4uxvhNQ65OCSFE4KEY4+5pj0X5NNqqTaNfLo0xPjAm5HNGX78XeFxMKappn5mfXDOG2cBVIYT/iSGXlBETnKGfw/iQA+wMHJ/02KpSOzMPIQTgD8ABk2xmyNWwEMKRwLeBUPOtnUb/c3XN60PA0THG27s9NhVDzRn6RK6LMT47qfGMleaZ+eHA0km+vxH424TGomJYC8yjEu+x/1TVvr4IWJPoCJV3rwMem+T7R4YQ9k1oLOOkGfN/YPxblFrVKZcXJTMc5V2M8XfAH5vY5b9jjA91azwqlgnmyCfSC7wnmRGNl0rMQwhLgBPY8e1wLYOuZn0c2NDAdhuo/GFKU2ow5ADTgTeEEOZ2f1TjpXVm/g6mDnnVbOBrIYT9ujgeFcdXaezYehS4rstjUQGEEI4ALmPqkI+V+DLFxGM+egHhncDMBjbfAKyncrb1YDfHpWKIMW4EvkTl4mY9G4GLvLiuBt0OfIrKXPmjDWw/Fzh7dJFHYtI4M69djlgrUvk/7A7gdGDXGOPfj/6RSo24BNg6yfd7gM8mNBblXIxxdYzx7cDuwN8DDzD1VF7iyxQTXZo4xXLEzVRC/gPgXOAnaS2+V/6FEG4DDqrz7YEY4wuSHI+KI4TQAzwH+Efg6cA0KnPltRJdppj0mflEyxGrUymfAJbHGF8YY7zJkKtN9S6EeuFTbYkxjsQYvxdjPI7KCcNnmHgKJtFlikmfmX+T7WvHH6MyD/4R4P/GGCdbiC81JYQwG3iYyvzlWA8BS5wvVyeFEOYDr6Wy5Hohlc87bAU+FWN8VxJjSOzMPISwO/A/gGHgGipLEw+IMX7BkKvT6lwI9cKnuiLGuD7GeAmwN/BS4EdU1py/MYTQzCqYliV2Zh5C6AXeBHwnxnh/Ir9UpTZ6s7afUVneCpXrMvv4QSElIYSwHDga+HwS08ap3jVR6raaC6Fe+FRhpX3XRKnbqhdCvfCpQvPMXIU2eiF0kMqKKS98qrCMuSQVgNMsklQA02pfCAODs4FTgWey/XFu9awHfgxcGfv7JrvHr9R1YWBwPpVj90gq63wnswa4Ebgq9vdt7vLQpEmFgcHFwKuAv2H76qt61jDBsbtDzKncHewZTYzjCODY0YFIqQgDg9X7rRzcxG5HUflU8lu7MiipAWFgcB5wJbBPE7sdRaW9b6m+MG6aJQwMPonmQl71tDAwWO8+GFIS/obmQl71rDAwmMqTYaRRz6W5kFcdFwYGl1W/qJ0zf0IbA3piG/tK7Wrn+PPYVZo6cuzWxnyiO381akYb+0rt8thVXrVz7G7bd6I58/FevPfycV9v2Rx47ivW8O5PrGxjAFL3bdkUuPBdu/Lbm+eyYV0vuy7dwmvOHuSoFzbygAEpXQ/cNY2Lz9qNO2+bzbTpkcNPWM/p569k2sTtnzrm37j3jm3//bH1gVMP3J9jXrS+YwOWumVoCPqWDPEv37qXPfYZ4qbvzOXj71jC4550N3sum+xJRFL6Lj5rNxbuMsyVv/sz61f3cPZL9uLrly7ilDPWTLR5c+vMr//qfObvPMShx/rUH2XfnHmR0z64ij2XDdHTC0ef9CiLl2zh9ltnpT00aUoP3z+dY05az8zZkb4lwxxyzKPce3vdx202F/PrvrKAY1+0juBnjZRDq/7ay4p7Z7Dvk7akPRRpSi98w2p+9I35bHw08NB90/jVDXM57Nl1pwgbr/KDd0/jj7fO4YRXr+vIQKUkbd0C575pD445aR37HmjMlX0HH72R++6Yycn7Lef1T13Gsidv4rgX1332aOMx/96VC3j8IRtZut9kD8qVsmdkGD5y2h5Mmx4541+9l7myb2QYPvDKpRzx/PV8/S93cNXv72TD2h4uO3txvV0aj/mPvrGQ409e25GBSkmJI/DRt+7O2sFpfPDKB5nuKkTlwNpVvaxaMY2XvH0NM2ZFFi0e4bmvWMcvf1T7GMRtGov5bTfOYvXKaRx/iqtYlC/nn74b9/95Bh+6+n5mzfEWocqHnXYdpm/PrXzz8kUMbYV1j/Twg6sXsPfj695HaOqliQDf//JCnvac9cxd4B+D8uPBu6dx3dULmTYjcuqB+297/S0fXsHzX+2JibLtnM8+yOXn7Mq3Pr0zPT2RJx3+GG8/r+7nexqL+XsvdZ5R+bNk3yGuefhPaQ9DaskBh23mgmvva3Rz1xhKUgHUxrydaRSnYJRXHrtKU0e6WxvzduYRXX+uNLVz/Dl/rjS1c/xt27c25jcDwy38wBHgJ20MSGrXj1vcbwtwSycHIjXpxhb32wL8v+oX42Ie+/vWAB+guaAPA/8n9vcNtjggqW2xv+9B4DyaO3aHgHNif593UVRqYn/fz4EvNrnbEPD+2N+37ROhIcYdp2vCwOAiKo8kWgCEemOgcor/09jft7rJgUhdEQYGd6Fy7M5l8mN3LXBz7O9zelCZEAYGl1B5YtYsWjh2J4y5JClfXJooSQVgzCWpAP4/Vte8hUm3Po0AAAAASUVORK5CYII=\n", "text/plain": [ "