{ "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//AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAbT0lEQVR4Ae2d/23cOLeG7Qv/HQQbIAWMO3BuKojTQZwOctNBgFQQOB04qWCx6WB3KzDWHfgrIMAaQRrwfU9WvFfWcnQ4PygeZh4CtEY8ks7Lh9IZmtJQx0dHR2fKfynn0pf7+/uLnIEyCEAAAhDYjMDx8fGt9ljl9lKsPT4ZGT7qs208Tv8Zr/AZAhCAAAR2InCZ2fulyl5Z+TggXylCE4AztCiCAAQgsA8CirGfpsdRr9mKfgTk/5oavXXtfGbdbuXH3rb7trf07dUlsrbetUdmizbv7MrbI3PLK35YWkv/uIf80ONoTc4t+H5WvlP+b+XsGIjK955a+vYqE1lb79ojs0Wbd3bl7ZG55RU/LF1Kv93Uu1deqTt95GVt927Y/rG37b7tLX17dYmsrXftkdmizY8ZufMvMrec3mnZPvXrWP+jLBf3RxsPWWhHEgQgAAEIVCBAQK4AlUNCAAIQ2IYAAXkbauwDAQhAoAIBAnIFqBwSAhCAwDYECMjbUGMfCEAAAhUIEJArQOWQEIAABLYhQEDehhr7QAACEKhAgIBcASqHhAAEILANAQLyNtTYBwIQgEAFAgTkClA5JAQgAIFtCBTNZWEH1u+4f9PC5rSwuSws/aWyGy2v9ZM/m7qzWmrp26tUZG29a4/MFm3e2ZW3R+aWV/ywtLZ+m/ctTVB/qsDK9JsP+bMGAQhAoCoBBXmby8KmPz5myKIqag4OAQhAoJwAAbmcFVtCAAIQqEqAgFwVLweHAAQgUE6AgFzOii0hAAEIVCVAQK6Kl4NDAAIQKCdAQC5nxZYQgAAEqhIgIFfFy8EhAAEIlBMgIJezYksIQAACVQkQkKvi5eAQgAAEygkQkMtZsSUEIACBqgQ2Csj6id/lkN9peaW8qqpudPCWvkcysh8ja8sK3rBQ9TtTvlW2uUwWT63971Lhltpb+t6Fme0bXXstfZtMLvSXOH3Q762/DMDs4rQJhl7WngNDPpr5trrOpcja5nR7NtXL2vez8p2yTSi12JevfNkF2dS/adg2tdTe0ve2vNJ+0bUvpc8mF7pXXimwHuWybDb5xe3UprJL5d+n5ftcb+nbq0dkbZ72Teyq5ztlO0ceb7LfvrZt7X+XerTU3tL3Lsxs3+ja96lPx7L4qmrfH5UOWVxoB5tqc5quVXA+fHNMbftab+nbq0NkbZ527BCAQDACpQH5XLrtX9dp+jYUmL1Waunbq1NkbZ527BCAQDACbkAu7P3+UqNeLX179YmszdOOHQIQiEnADciSnYJt6g2Pa5J6zXYDpkZq6durT2RtnnbsEIBAQAIlAblE9pOSjSpt09K3V6XI2jzt2CEAgYUJnMjfo8FnWk4lpF7wtNzWUy/x75xxD2UtfXvyI2tLj439qUps8t/Lhe703ngV/9ntw3AU7DZsaLhtCOz/N3+aPlpAnk26QL8JtG2Tu7BTWZV38bX0PQtFxsjaTLvp0+KZfSZtRgB2m/FKW8Mtkdh+aUMW34fd0zJ3tD9UuMoYUg/Z7LVSS99enSJr87RjhwAEYhD4mmSUjiH/ph3s11rTZD2wm+GbcWrb13pL314dImvztGOHAASCESgKyAq4n6T7TkMXr5L+YbzotdbfpLIay5a+vfpE1uZpxw4BCMQj4I4hjyRbb/i9AnEauniu9RcKSkvcBGrpe4Qg+zGytqzg0kK1tf0HYPcJ0n9HNneJtfe12v1j6XG23a61/211234ttbf0vQuz1txKtNdma3frzpRt8p5TXWRVbs6VVJRtIAABCBwiAQV5m8viSvH3uGjI4hAhUWcIQAACSxMgIC9NHH8QgAAE1hAgIK8BQzEEIACBpQkQkJcmjj8IQAACawgQkNeAoRgCEIDA0gQIyEsTxx8EIACBNQQsINujbm+V5ybLWbM7xRCAAAQgsCMBm4LBYvCPVzittLxSTvNSWDkJAhCAAASWIWBvHrIYXPxOvWVk4QUCEIDAARNgDPmAG5+qQwACsQhsFJD1E7/LIb/T8krZhjsWS/J3pnyrnOZhXsz3ro6iapeupm3qce1AX9hzMuo557W52Q9Ve/HkQgJk81180O+tvwzALCjaZDMva86BoeObn8/KdtPRJrlZ9EtA/rZO0bW3atNSoFH1RW7XyNq8dkf7P4RscqF75ZUCq71p4l9ZNpv84nZqU9ml8u/T8lrr8vVO2bQ+ruWj1nGjaZeeEG26jnd0fUl3tHZNumwZWdtYZ+7zIWlXXe1aFIb74pt6F9ohN83mtcrPh282fSR1RCB6m0bX11FTI7UXAqVjyPZYRu455W9DRc1O6otA9DaNrq+v1kZtFwTcgFzY++UZ5i6a+x+R0ds0ur6OmhqpnRFwA7Lqk4Jt6g2Pq5h6zXbjjdQPgehtGl1fPy2N0q4IlATkkgo9KdmIbboiEL1No+vrqrERG4PAiWQ8GqSk5VRZ6gVPy2099WT+zhl7Lxv+df5T9djkP4AL3S29CV736G0aXV/w5m0jr+frpbH2p6nFLCB/H1bSMtl+LBVcvkmsfc4FpVT2U76Lz+quettLTH+qFL1No+v7qU6GPVam5+ulsfavqRlKhyxsNqJV2mm0TD1ks5P6IhC9TaPr66u1UdsFgdKAbK+DT6+CH1fMeo83w7fLuJzP8QlEb9Po+uK3MAq7I1AUkBVwP6lmdxq6eJVqOIy5vNb6m1TGsh8C0ds0ur5+WhqlPRGwMeTSZL3h9wrEaejiudZf6MKpfgNLPq23ZOPVqZduc2iY32v5/6hl2BRce7M2LWywsPoit2tkbV67H7p2u1t3pmwTB50quP2UN+e8kwA7BCAAgVYE9CVkc1lcKf4eFw1ZtBKKXwhAAAKHRICAfEitTV0hAIHQBAjIoZsHcRCAwCERICAfUmtTVwhAIDQBAnLo5kEcBCBwSAQIyIfU2tQVAhAITcACsj3q9lZ5bkKX0JVAHAQgAIGOCdg0ARaDf7zCaaXllXKal8LKSRCAAAQgsAwBezuOxeDid+otIwsvEIAABA6YAGPIB9z4VB0CEIhFYKOArJ/4XQ75nZZXyjbcsUhq6durYGRtnnazS/+Z8q1ymt+6ZLfFtomsD23bnQaRuZXUqJb+4smFJMDmu/ig31t/McFat4vXJvl5WXsOjJa+ra5zKbI2R7e132dlu5lrkzYt9uUqX24azq+Q+tDmNl92g8jcsoInhUvpt8mF7pVXCqxHuSybTX5xO7Wp7FL592n5Ptdb+vbqEVmbp31sVz3eKds58HhcHuVzZH1oy8cM79yJzM3TbvZ96texLL7qsPfFN/UutENums1rlZ8P3xz6WCW19O1VKLI2Tzt2CEAgGIHSMWR7LCP3nPK3oT5mr5Va+vbqFFmbpx07BCAQjIAbkAt7v1WeYW7p22unyNo87dghAIGYBNyALNkp2Kbe8LgmqddsN4hqpJa+vfpE1uZpxw4BCAQkUBKQS2Q/Kdmo0jYtfXtViqzN044dAhBYmIAF5EeDz7ScSki94Gm5rade4t854x7KWvr25EfW5mnHDgEIxCHwNEmxgPx9WEnLZPux1KMYaagiNyyRyqq8i6+l7wcQMiuRtWXkUgQBCMQl8DVJKx2ysNmIVmmn0TL1kM1eK7X07dUpsjZPO3YIQCAYgdKA/Jt026+5psle034z6i1O7ftYb+nb0x9Zm6cdOwQgEIxAUUBWwP0k3Xd61OtV0j889vVa629SWY1lS99efSJr87RjhwAE4hE42UCS9YbfKxCnoYvnWn+hoHSzwTG23bSlb09zZG2z2tWW1sO3+wDpvx+bm8Ta81rt+nF25wWMkfWhbbsTIDK3khrV1n8sEWfKNnHQqS7CKjfnSirKNhCAAAQOkYCCvM1lcaX4e1w0ZHGIkKgzBCAAgaUJEJCXJo4/CEAAAmsIEJDXgKEYAhCAwNIECMhLE8cfBCAAgTUECMhrwFAMAQhAYGkCBOSlieMPAhCAwBoCFpDtUbe3ynOT5azZnWIIQAACENiRgE3BYDH4xyucVlpeKad5KaycBAEIQAACyxCwNw9ZDC5+p94ysvACAQhA4IAJMIZ8wI1P1SEAgVgENgrI+onf5ZDfaXmlbMMdiyX5O1O+VU7zMC/m23MUWducdulu2qZz2szWgT7OSa8RM/Zer5dUlVr6iycXkgCb7+KDfm/9xURp3YKiTUbzsuYcGIOfz/JlNx1tEpxFvwTkb22KrG2t6JFB+pu06UjC7Meo+iK3O9pmT6mdjEuxtcmF7pVXCqxHuSybTX5xO7Wp7FL592l5rXX5eqdsWh/X8rHtcSNry9VJekO0aU6blUXXl3RHbne05eNZartdlvtkO5zrknNffFPvQjvlptm8Vvn58M2hj6SOCERv0+j6OmpqpPZCoHQM2R7LyD2n/G2oqNlJfRGI3qbR9fXV2qjtgoAbkAt7vzzD3EVz/yMyeptG19dRUyO1MwJuQFZ9UrBNveFxFVOvOdxTD2ORfP4XgehtGl3fv4BSAIF9ECgJyCV+npRsxDZdEYjeptH1ddXYiI1BwALyo0FKWk6VpV7wtNzWU0/m75yRsrAEordpdH1hGxZhXRJ4mlRbQP4+rKRlsv1Y6lGMNFSRG5ZIZbyL7wG12CvR2zS6vtiti7oOCXxNmkuHLGw2olXaabRMPWSzk/oiEL1No+vrq7VR2wWB0oBsr4tPr4ofV+yZVm5GPZqxjc+xCURv0+j6Yrcu6rokUBSQFXA/qXZ3ehzpVarl8GjSa62/SWUs+yEQvU2j6+unpVHaE4GTDcRab/i9AnEauniu9Re6cG42OMZWm8qn9ZZsvDr10m0ODfN7Lf8ftzronnaKrK2gis3atECbbRJWX+R2R1vh2bXFZrXZHkvTmbJNMnOq4MbNuS0aiV0gAAEIbEtAQd7mlblS/D0uGrLY1hH7QQACEIBAOQECcjkrtoQABCBQlQABuSpeDg4BCECgnAABuZwVW0IAAhCoSoCAXBUvB4cABCBQToCAXM6KLSEAAQhUJWAB2R51e6s8N6FLVREcHAIQgMABE7BpAiwG/3iF00rLK+U0L4WVkyAAAQhAYBkC9nYci8HF79RbRhZeIAABCBwwAcaQD7jxqToEIBCLwCZzWRzpJ36Xg3ybkP5U+XKpn1u39O01WWRtvWvvma2xl36bmsDmYnmmayXNLW6m6qmlb69ykbV52j37LnUrDshyYvNdfNBJ9cUEad0m+7FJfl7WDsotfVtd51JkbXO6zRZde3R96/hKt10bn5XtRrlNiGX3aRZJLX17FYyszdPu2fdVt6IhCzmzyS8ep2Bs4oZvewvOPwajraxGaunbq09kbb1r75mtXRvKF8p25/xXry32aW/p26tHZG2eds++r7oVBWSJuVDOTbN5rfJzXTzWI6iVWvr26hRZW+/ae2brsccOgSyB0oBsj2XknlNOY2Jmr5Va+vbqFFlb79p7Zuuxxw6BLAE3IBf2fqs8w9zSd5bWqDCytpHM7Mfo2qPry0KlEAJ7IOAGZPlIwTb1hsduU6+51pBFS9/jeuY+R9aW0zsui649ur4xSz5DYG8ESgJyibMnJRtV2qalb69KkbX1rr1nth577AdK4ET1fjTUPS2nKFIveFpu66knY88l10gtfXv1iaytd+1N2Q5DJn8K4ib/+dlTFTceeOwxCTRu86eJigXk78NKWibbj6U9ziGx9jl3cqayKu/ia+n7AYTMSmRtGbkPiqJrb63P/AvYswfQWPmpCTRu868JbumQhc1GtEo7jZaph2z2Wqmlb69OkbX1rr1nth577BDIEigNyPbTT/vF0TRZL+Jm+HaZ2va13tK3V4fI2nrX3jNbjz12CGQJFAVkBdxP2vtOQxev0lGGMZfXWn+TymosW/r26hNZW+/ae2brsccOgXUEbAy5NFlv+L0CcRq6eK71F7pwlriR0dK3xyeytt61d8tW14n18O0eS/rP0uZ9sWvlWtfMR69hdrG39O3pjqzN0+7Z91E3u1t3pmwTB53qRKlyc86rCHYIQAACh0pAgdzmCrpS/D0uGrI4VFDUGwIQgMCSBAjIS9LGFwQgAIEZAgTkGTiYIAABCCxJgIC8JG18QQACEJghQECegYMJAhCAwJIELCA/Ghym5ZL+8QUBCEDg0An831wWFpDTHBZpeehwqD8EIACBJQlsPJfFkuLwBQEIQOAgCTCGfJDNTqUhAIGIBAjIEVsFTRCAwEES2GQuiyP9xO9yoGQT0p8qXy71c+uWvr0zI7K23rVHZyt9NvWAzVvxTNeCzaMcJkXW5kE6VO3FAVmAbL6LDzrpvhhMrdvEKTZhysvaQbmlb6vrXIqsbU632aJrj6pPuuzc/6x8p2yTB62UQ6TI2jxAaD86KhqyECib/OJxCsYGdugNWHC+svVaqaVvr06RtfWuPTJbO/eV7ZVNb8X5V4/1kvbI2jwOaC8MyAJ5oWxTB07TtQrOdfFYj6FWaunbq1Nkbb1r75mtxx47BLIEinrI2vNc2f49m6Y0Zmb2Wqmlb69OkbX1rr1nth577BDIEnADcmHv95fs0XcsbOnbkx5ZW+/ae2brsccOgTkCbkDWzinYpt7w+Hip11xryKKl73E9c58ja8vpHZdF1x5d35glnyGwNwIlAbnE2ZOSjSpt09K3V6XI2nrX3jNbjz32AyVwonrba5vsbnHq7U5RrCu37VJPxp5LrpFa+vbqE1lb79p7Zuux/2ntw1DTn6rgJv8x29MqN62hNNb+h+pvMfjIAvJK2R5ds8J/DUvYoygSK1MWcgJf5V18LX1bhedSZG1zus0WXXt0fR7fQ7Vbu6nuz3qsf2PtdgPbYvCn0iELC9YWuKcp9ZDNXiu19O3VKbK23rX3zNZjjx0CWQKlAdl+GppeZz4+kH0b3gzfLuPyfX5u6durR2RtvWvvma3HHjsEsgSKArIC7iftfaehi1fpKMOYy2utv0llNZYtfXv1iaytd+09s/XYY4fAOgI2hlyarDf8XoE4DV081/oLXThLDMi39O3xiaytd+1h2eo6sB683UNJ/znavC52LVzrmvjoga9pj6zNq/eha7e7dWfKNnHQqU6kKjfnvEbADgEIQOBQCehLyOYKulL8PS4asjhUUNQbAhCAwJIECMhL0sYXBCAAgRkCBOQZOJggAAEILEmAgLwkbXxBAAIQmCFAQJ6BgwkCEIDAkgQsID8aHKblkv7xBQEIQODQCTxNACwgfx9W0jLZWEIAAhCAQH0CX5MLhiwSCZYQgAAEGhMgIDduANxDAAIQSAQIyIkESwhAAAKNCWwyl8WRfuJ3Oei1CelPlS+X+rl1S99eG0XW1rv2ntkae+m3qQls3otnulb+Nd+4bVMrtfS9a52ia6+lrzggS4DNd/FBJ9UXg611m1jFJlR5WTsot/RtdZ1LkbXN6TZbdO3R9a3jK912bXxWvlO2yYdWyouklr53rWB07UvoKxqykBCb/OJxCsYGfvi2t+BsM91XSy19e5WKrK137T2ztWtD2V5NZK/l+dVri33aW/retR7RtS+hryggC/SFsk0tOE3XKjjXxWM9glqppW+vTpG19a69Z7Yee+wQyBIoDcj2zif792ua0piY2Wullr69OkXW1rv2ntl67LFDIEvADciFvd9fskffsbClb096ZG29a++ZrcceOwTmCLgBWTunYJt6w+PjpV5zrSGLlr7H9cx9jqwtp3dcFl17dH1jlnyGwN4IlATkEmdPSjaqtE1L316VImvrXXvPbD322A+UwInuHNrNOnuV07qUesE5e+rJ2HPJNVJL3159ImvrXXtTtsOQyZ+CuMl/fvZUxY0H/me2w2271tV5Yy+Rtnx0IohnWq59p5496qFtbNvcyZnKqryLr6Vvq/BciqxtTrfZomtvrc/8C9MzjyP2hwTg9pBH6Zri68bv1PtDB19lHKQestlrpZa+vTpF1ta79p7ZeuyxQyBLoHQM2X76mV53Pj6Q9SJuhm/Gcfk+P7f07dUjsrbetffM1mOPHQJZAkUBWQHXxjfu1LV+lY6izzZc8Vr5TSqrsWzp26tPZG29a++ZrcceOwTWEThZZ8iUW2/4vQJxGrp4rvUXunCWuJHR0ncGxYOiyNoeCM2sRNceXV8G6T9Fuk6sh2+dlvSfpc37YtfKta6Zj/9sVedvS9+71ii69tr67G7d7E29XQGzPwQgAAEIrCegIL/xTb31R8MCAQhAAAJ7IVA0hrwXTxwEAhCAAARmCRCQZ/FghAAEILAcAQLycqzxBAEIQGCWAAF5Fg9GCEAAAssRICAvxxpPEIAABGYJEJBn8WCEAAQgsBwBAvJyrPEEAQhAYJYAAXkWD0YIQAACyxEgIC/HGk8QgAAEZglsMpfFkX7idzkczSakP1W+1O/yq8yFPFXd0vdUy3Q9srap1ul6dO0d6LOpB2zeime6Fmwe5TBJ7NBWqTVqsS0OyBJgk9h/0En3xeqodZs4xSZMeVk7KLf0bXWdS5G1zek2W3TtUfVJl537n5XtzSY2edBKOURCW71mWIJt0ZCFhNjkF49TMLYqD70BC85Xtl4rtfTt1Smytt61R2Zr576yvbLprTj/6rFe0o62erSXYFsUkFXFC+XcNJvXKj/XxWM9hlqppW+vTpG19a69Z7Yee+wQyBIoDcjn2tv+PZumNGZm9lqppW+vTpG19a69Z7Yee+wQyBJwA3Jh7ze9Wy/rZNvClr49zZG19a69Z7Yee+wQmCPgBmTtnIJt6g2Pj5d6zbWGLFr6Htcz9zmytpzecVl07dH1jVnyGQJ7I1ASkEucPSnZqNI2LX17VYqsrXftPbP12GM/UAIlATn1gnOIUk/GnkuukVr69uoTWVvv2ntm67HHDoG1BNyAbI96DHvnhiVSWZUfh7T0vZbYYIisrXftPbP12GOHwBwBNyAPO/+h5SpzoNRDNnut1NK3V6fI2nrX3jNbjz12CGQJlAZk+2loep35+EDPtHIz6tGMbfv63NK3V4fI2nrX3jNbjz12CGQJFAVkBdxP2vtOjyO9SkcZHk16rfU3qazGsqVvrz6RtfWuvWe2HnvsEFhH4GSdIVNuveH3CsRp6OK51l/owrnJbLvvopa+vbpE1ta79rBsdR1YD97uoaT/HG1eF7sWrnVNfPTA17SjrR7d2myPJf1M2SYOOtWJVOXmXD08HBkCEIBA3wQU5G2uoCvF3+OiIYu+q4t6CEAAAn0QICD30U6ohAAEDoAAAfkAGpkqQgACfRAgIPfRTqiEAAQOgAAB+QAamSpCAAJ9EBg/9naru31T1V90588mCidBAAIQgMCOBBRjb3WI1brDWEC2R93sVTS5xGNwOSqUQQACENiOQHpRdHbv/wW+16yOj3UhPgAAAABJRU5ErkJggg==\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": "iVBORw0KGgoAAAANSUhEUgAAAucAAAAVCAYAAADlwkDzAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJVUlEQVR4Ae2dgXHdNgyG61wHyDUbeIS02SDdIOkGTTdIZ8gI7grpBk5HaDZwNmjsDVz8CunQFCkAJCSSFnWnSiJA4sdHUuFT9fwu7u/vfwi3i4uLSyr7EpaF55w99G1x3oO+HjTk2PegrQcNOT4W5T3k14MGC5a5NnrIrwcNOT4o70FfDxq2GNXYesitBw1bDHvQ14OGLUatbT3w6UHDnv3A5ZeyPwsFkcN7un4ZloXnnD30bXh+6XQ2kTAAo8ln/5ExGU/G+xPgI8xxyDOq8Zh8eXqTEc+otcfso/17gGO8sj8szmlR+Yb0vaCn5n+ndG7ZyfbB7e/peEX7ZaoNizIuFun/hDjkh3wO3VzMJENOt6VQivWS9hvan8ftTj4xkbLrERmXZfq41lbejz3rrijO5j2l13GMrCejh/tv8l5YNzK+1T6Q8XD30iPHIBer53laOw6PGoO1Okfuo7MwTs4TKkTfYRH3L85T+5Yd9Wh/4+s53xs64vWYZHul5ZpYzvd5aSxtPZd3kqFGtzau93fxP9LxysVDx2bzdz5Zu2/X6ggtiJlqz2k5ZAyl4kvLXA5DMpbmmPLT5p1qQ1OmGQ/Ot/k4noy+3+sdi+Rc14yD2Pcoxto45I9/A081BkdlFI8p7bU2b237lv5areR/unFcy9uacdgHy+KZCrCge5cTmrOjDu03cT0q+0D7dVxec62N5fyvamJq6lK8JEOtbk3MnC/FxOtJ3OIcfXdKPjlumvKRGGvy4nwleXNtbNmpfdU9xfk3H8dhTpNR+l4YMqo935ux1yeJc8Yx6PngOCqjMIeSc0neJe3uUUeilXy6WBOE+Ut0h/4tzyVaOcah3S/Ob7eSogpJO5Vf0/4xrktleKVkc3EY1+GuS2JRnVvaD3migVipHEp0p9rRlFFMdnGO9qCZ9tPx0bDM+Y7EOJdDSbk075K2UYc29T2F6jQfx2G+k1H6Xhgyqj3fm7HXJ41ztjHo+eA4KqMwh5Jzad4lbVvXkWolv1PdSy05WzH2ffCM3ul5TRdbf51lyw7bV9rj7c4VwG61lcRCXr9ZCci1wzAs0Z0LZV0++VgTXbfXA+O1qn5LSubL2RhPRv2N37ONwZIemIxKqB1bZ/bR/rw5xosdXwh9S/vyJcqMpqSdFqSrLxsm6v+UKFMXVcRCXr+qA+orNGekl7zUmHwKwSmqNWWs0NncddR5fiS4yehI2qpYc57zuCYjnlFrj9lH+/cAx3ixY3H+C+34Amduy9n9wts/JQ/r+6fpkgV8WC93Xhpr+WJqrlHD8h4YlaQz+ZRQ09VpzVintq33qPP8SGqT0ZG05bHmPOdZTUY8o9Yes4/27wGO8WLH4hwLaL+YTsni7Kk6vuyFPzngmIqFvKw+IGylMAqjOIfJJyZifz0CY/us92tx1Hm+H5F1y5PRmsneJXOe84QnI55Ra4/ZR/v3AMd4sf9IOvAk5m5DT86OBnKbf7rzX85BWV4aC+/u7PY314McemAUyBGfnoqPeyXgH6Kj+cD2lr408llMdO3YmjH+7jbyPTrvNQm+ZNR5zmdm5zEso4HGYUlvNZ/nJaIPrtOc0UhjsJHW5n105JjslPHSB1icF220YLmjxFA3tdDxZQhSvVXEyi2aqzVJGqjQLWnewudUfNAfBO1nC3CKNpoyhs5GeSsQfXOtmC/NGauTLawwMqNRxmFh15xmDBbyQbXmjEYag420Nu+jivGlrtop46UP8FoLFtB+MZ1KbsuOF9dTT6bRODbYrbaSWMhr60mTlbZeGGnzmXy0xPT+PTDWq25XY9R5fiSxyehI2rJYc57znCYjnlFrj9lH+/cAx3ixY3GOxWtqge0lbtnxa4n4MmS84enkZ/epZLG5/30Q+2muxbGCRvEhgX16b6CtF0ZB6qLTyUeEqcqpB8ZVCWgqG8ylUee5GFOvjAx0gcHWvVDMqNbRKBeNjKHmeQM+YDkUI03np3wbMU5J0ZQN1UdPlPHSB1icf6b91UbvZe20+P6L6n0lQG98fQcLf1v896jslmz4ediiTRorahwfEvCjJtnN6a3SRo13wSibZN4w+eTZWFmaM7ZKhGvHYi6NOs85Nt7eKyMLXS7H7L3QM9j7aJiLRuow87wRH7AchpGm41O+DRmn5GjKhumjJ8z4oQ/wgxY39I/i8muh8ZFsnB2P4D/Qjl+lxI4nXy8T7eDPw2Av/kVK1KWdjeVju3iX/jp3dH7F2qh+F4xIB9jjwwh+5Qs/3YacUPY+lbuzn4ZPioG2zPEcjrE2z9hfmrcbU8VzCXFpG3Wei+Zfr4yMdG3eC+Nxpb0mjYcwlsbx+h27Hu6lh/Bx81QUqzdGXk/pUTo23JiouheWavT1pFoDf+g91Tj2uZcerRm7cXO5LMjdxWpB7cVydu/HHakdPGEvXpxz7Yd2ioNXdbIfOkJfnNNWpQ2xaM8yjOPlrmt15NqNyynO5JP5QBqzKr0elXFpvr7eUWMY8UZl3CsjC13Uhsm90I+n0qNFLpLYFGfeS5l76aiMJP2/5XPUGNzSILWN2kdPiXHYB3itBRueRv+xnKX/w9nTtdalr2ig3K2Ldyn5k1qFbulWq200RpOPdGSU+43KuDzjbzVr55Im/qiMe2VkocvqXqgZBylfi1xS7cZlcwzGRNbXozJaZ6IrOWoM6lSlvUfto6fE+KEPlsU5LZjx7vglvcODJwCrjbOvKiQK3PtBVn/3PBHhe5HLA/9rBnmxm4W2kRhNPuyQqHYYlXFt4hZzSaphVMa9MrLSZXEvlI6BnJ9VLrn2ffkcg55E/jgqo3xGMstRY1CmZttr1D56SoxXfUA3Uv9qC97zvPbX8ZFsm/bYP76m+u/isr2ukQf0Stu30oaYiC2NG/tZ6Yjbja+hEVrj8ty1lS7EROxcHK7cSgcXx8KOPJGvtC2r3BATsaVxrf2s8pDoQp7IV+ILHyttiEl7MWMrHZK8oRN6hb5m92jERGxJ3D18KLZZLlv6kCNy3fIJbVa6EBOxw7Y151Y6JDGhE3olvvCx0oaYtBczkurN+VnlkWvfshycwEvaplVuiInY0rixn5WOuN09rpEn8s21Hdsv4Og3t3J/TWXJJ86c3bfT8kga31H8T5QD+ycU99DZO6PJZ49ef9zmZPyYxx5XkzFPdTLiGdV4TL48vcmIZ9TaY/bR/j3AMU7Z/wd4PfpmzgVlQgAAAABJRU5ErkJggg==\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": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAF5CAYAAABtIcr0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXRlVYGo8W8nNc9AioJiLkAUlEFaRZBB1MaIrx1BbURQVJwQRX1Ka7e+9xwaFWgFRcSxFXyNc2sThxYUROQpKs4IMoNVRYoaoaYk+/1xc6uSW7nJHc/4/dZiQW7OSbauU1+du8++54QYI5KkfOtJewCSpPYZc0kqgGm1L4SBwQA8EzgeWAiEOvtGYB1wPfDD2N830q1BSo0IA4PTgX7gaGAuHrvKiTAwOAs4GXgKMJvJj901wHWxv+/asd/YIebAPwOnNTGOFwNXj+4npekS4OlNbP9i4KvAe7ozHGlqYWCwF7gCeHITu50aBgY/F/v7Lqi+MG6aJQwMLgZe1sJ4TgkDg7u3sJ/UEWFg8FCaC3nVKWFgcGmnxyM14ak0F/KqM8LA4KLqF7Vz5k+Y4LVGBOCwFvaTOuXwNvb12FWaDm1xv17g8dUvasM9s+XhtLev1K4Zbezrsas0tXP8bTvuJ5oz39F5z96LO343i97eytc77TrE5355VxsDkLrvq5cs4rqvLuT+O2bw1Oes5/zPLE97SFJTfnDVfP7j4l1YtXw6C3cZ4tx/W84TT9g40aaNxRzgrPeu5HmvXduxQUrdtsvuQ5x67ipuuW4uWzbVWx0gZdPPvzeHL31oMe/41IM8/qhNPPTgpL1uPOZS3pz44g0A/OU3s1j1N4915ctVH+njlDev4tBjNgGwZK+hyTZv/GLnlR/p49QD9ufNz9ybX/5odnujlCTVNTwEd/1xFmtX9XLmEftx2iHLuPjNu7Lp0brvMBuL+Sv/5SE+f8udfPn3d3LSaWv4wKv25L7bp3ds4JKk7VYt72V4CG66Zj4f+e69XHrdPdz1x1l88YO71NulsZg/4ehNzF0QmTErcvIr1/GYwzfy84G5HRu4JGm7mbMrd0A8+ZWrWbzHMDvtOszzz36YX/+4bndbvDdLiHi3RUnqjoW7jLDTrkOExq/bTx3zdQ/3cNM1c9i8MTC0Fb73pfncdsscnvysR9oZq9R1Q1th88bAyDCMDLPtGJby4OkvWst3P7+IVct7Wbuqh/+8YieOPHFDvc2nvsI/tDXwpQv6+PDrZ9LTE9l9vy2864oH2Pdg/1Qo2774gV34+ie2zzHe+N0FvOiNqzjrfatSHJXUmDPes4p1q3t57VP3Y/qMyFH963nF+Q/X23zqmO+8ZJhP/uTejg5SSsJZ7zPcyq/pM+C8S1Zy3iUrG9nc+5lLUgHUxrydq5peEVVeeewqTR05/mpjPuFn/hv0aDsDkdrUzgV5j12lqSPdrY35r4BWLmwOA79sY0BSu25ucb9h4BedHIjUpFaP3Y3ArdUvxsU89vetBz7awg+9KPb3rWlxQFLbYn/fncAXW9j1Yo9dpSn29/0G+HqzuwEXxP6+bWf1IU7w4Z8wMLiMJp4BGvv77mhyIFJXhIHBQ2jiGaAeu8qKMDB4GHAUUz8DdDXw49jfd8+4/SeKuSQpX1yaqMILIXw9hHBJ2uOQuskzcxVaCGF34C4qFzp3jTF6GwoVkmfmKrpXU5lnHAFOTXksUtd4Zq7CCiH0AMuBxaMv/THGeEiKQ5K6JrEz8xDCziGE74cQzgghzErq96rUngWMPdb2DSEcltZgVB4hhHkhhDeGEL6TVO+SnGbZCBwDXA6sDCFcEELYI8Hfr/J5KzB/zNczgXNSGotKIISwfwjhk8AK4GLgicCWJH53YjGPMW4EPkVl/eR84FzgjtG/uY4JoYm7sEtTGL3weULNy73Ay0IIPiVLHRMq/j6EcB3weyrXaeZQ+TT9BTHGkSTGkfQF0I9TuRAFlbOkWcDJwPeB25yCUQdVL3zW8kKoOqI6lQLcS+UTnCdQaVr1+cgB+EJi40n6AmgIYQA4iYk/4bSByh/Ay4CPxxgfSHJsKoYJLnzW8kKoWhZC2B94G3AGlV5N9E5vCPhcjPHspMaVxtLED1H/DnfzGD8F812nYNSC2guftbwQqqZMMpVSb8puK3BhUuODdM7MA3A7sH8Dm0cqt3h8AHhajPGhbo5NxRBC+B6Vd3/1DANfiDG+OqEhKcdCCAcA1wGLqJxwNuKnMcZjuzeqHSV+Zh4rf3t8kMqUylQClZvO7IwfcFID6lz4rOWFUDUjULnpYKMh30ClcYlKK5BfaXC7EWAV8JQY44oujkfFUe/CZy0vhKohMcbbgWOp3GmzEeuoLOpIVCoxH12meDmTr7+shvyoGOOdiQxMuTZ64fMcJp8vr5oHvL27I1JRxBhvBY5j6qA/SoLLEcdK7eP8IYS9gduY+A/eCLAGeJIhV6NCCH8P/BcwrcFdRoAjYoy/7d6oVCSjF85/Tv0Tho3AbjHGRs/iO6bRg77jYoz3jl4ZfjbjlymOUHnHsDOwOY2xKbf+BFw9wev/OPrvq2peH6ayRlhq1BD1Qz4EfCmNkEPKN9oKIRxH5UyqemGhOrVyHJU/mAB7ut5c7QghRGBFjHG3tMei/AohHEJlWSJUPqb/Y2DBmE02AYfFGP+S8NCA9FeI3EDlHgYwfo78z2z/2+9+7+EiKU01Ie+NMf6aHefQf5lWyCHlmI9ZpriVmoudMcbNGHRJKZsg5COww0XRzaSwHHGstM/MobJM8fNMsGrFoEtKU72QV40J+qdJYTniWLl4OEUIYSaV+ShwDl1Ncs5crZgq5FmThTPzKXmGLilJeQs55CTmYNAlJSOPIYccxRwMuqTuymvIIWcxB4MuqTvyHHLIYczBoEvqrLyHHHIaczDokjqjCCGHHMccDLqk9hQl5JDzmINBl9SaIoUcChBzMOiSmlO0kENBYg4GXVJjihhyKFDMwaBLmlxRQw4FizkYdEkTK3LIoYAxB4MuabyihxwKGnMw6JIqyhByKHDMwaBLZVeWkEPBYw4GXSqrMoUcShBzMOhS2ZQt5FCSmINBl8qijCGHEsUcDLpUdGUNOZQs5mDQpaIqc8ihhDEHgy4VTdlDDiWNORh0qSgMeUVpYw4GXco7Q75dqWMOBl3KK0M+XuljDgZdyhtDviNjPsqgS/lgyCdmzMcw6FK2GfL6jHkNgy5lkyGfnDGfgEGXssWQT82Y12HQpWww5I0x5pMw6FK6DHnjjPkUDLqUDkPeHGPeAIMuJcuQN8+YN8igS8kw5K0x5k0w6FJ3GfLWGfMmGXSpOwx5e4x5Cwy61FmGvH3GvEUGXeoMQ94ZxrwNBl1qjyHvHGPeJoMutcaQd5Yx7wCDLjXHkHeeMe8Qgy41xpB3hzHvIIMuTc6Qd48x7zCDLk3MkHeXMe8Cgy6NZ8i7z5h3iUGXKgx5Mox5Fxl0lZ0hT44x7zKDrrIy5Mky5gkw6CobQ548Y54Qg66yMOTpMOYJMugqOkOeHmOeMIOuojLk6TLmKTDoKhpDnj5jnhKDrqIw5NlgzFNk0JV3hjw7jHnKDLryypBnizHPAIOuvDHk2WPMM8KgKy8MeTYZ8wwx6Mo6Q55dxjxjDLqyypBnmzHPIIOurDHk2WfMM8qgKysMeT4Y8wwz6EqbIc8PY55xBl1pMeT5YsxzwKAraYY8f4x5Thh0JcWQ55MxzxGDrm4z5PllzHPGoKtbDHm+GfMcMujqNEOef8Y8pwy6OsWQF4MxzzGDrnYZ8uIw5jln0NUqQ14sxrwADLqaZciLx5gXhEFXowx5MRnzAjHomoohLy5jXjAGXfUY8mIz5gVk0FXLkBefMS8og64qQ14OxrzADLoMeXkY84Iz6OVlyMvFmJeAQS8fQ14+xrwkDHp5GPJyMuYlYtCLz5CXlzEvGYNeXIa83Ix5CRn04jHkMuYlZdCLw5ALjHmpGfT8M+SqMuYlZ9Dzy5BrLGMug55Dhly1jLkAg54nhlwTMebaxqBnnyFXPcZc4xj07DLkmowx1w4MevYYck3FmGtCBj07DLkaYcxVl0FPnyFXo4y5JmXQ02PI1QxjrikZ9OQZcjXLmKshBj05hlytMOZqmEHvPkOuVhlzNcWgd48hVzuMuZpm0DvPkKtdxlwtMeidY8jVCcZcLTPo7TPk6hRjrrYY9NYZcnWSMVfbDHrzDLk6zZirIwx64wy5usGYq2MM+tQMubrFmKujDHp9hlzdZMzVcQZ9R4Zc3WbM1RUGfTtDriQYc3WNQTfkSo4xV1eVOeiGXEky5uq6MgbdkCtpxlyJKFPQDbnSYMyVmDIE3ZArLcZciSpy0A250mTMlbgiBt2QK23GXKkoUtANubLAmCs1RQi6IVdWGHOlKs9BN+TKEmOu1OUx6IZcWWPMlQl5CrohVxYZc2VGHoJuyJVVxlyZkuWgG3JlmTFX5mQx6IZcWWfMlUlZCrohVx4Yc2VWFoJuyJUXxlyZlmbQDbnyxJgr89IIuiFX3hhz5UKSQTfkyiNjrtxIIuiGXHllzJUr3Qy6IVeeTat9IQwMTgdeCDwNWACEOvtGYD1wI/CN2N+3uVuDlMaKMW4OIcwCNlEJ+p4xxgfqHrtvuxSmzVgYBgb/vfojgNXAdcB3Yn/fiCFXmsLA4Gwqx+7RwDzqd3cEWAX8d+zvGxj7jR1iDnwcOLGJcTxrdPvXNLGP1JaJgs41D70beMYOGx90JPT0TgeeUvOdfuCJIYSrMeRKSRgYDMCngKOa2O25YWDwkNjf99HqC+OmWcLA4IE0F/Kq48LA4MEt7Ce1bNyUy9Jl9zO09aSmf8iWTS9n3iJDrjQdQXMhr3pFGBicU/2ids788W0M6Alt7DtOCKHeWwxpnG1B3+exsPyeZQxtnejd5sS2bJ7ByvuXsfdBYMjVhA43qtXuzgQOrH5RG/MZLQ+nvX0JFU8LIXwNWBVC2K+dn6fyiDFu5m2feDnQeNC3bJ7Byvv2BeBD33imIVejQghHAqtDCF8IITyxAz+yI91tbDXL+8/cjZc9dn9euM8BvPLI/fj2pxe28cvHCSHsHEJ4C3A3MEDlIsB0YOdO/Q6VwJx5Q+yx7HZgfNDXDvbwyXfCG47r5fQnLOP7X54/LuR77P8Xps+IaQ1bubQb0Au8HLghhPDnEMJrQgjzu/Lb7vnzdP5hjwN5/5m7TbZZY29JX3rew+z9mBXMmBW56w8zOP+Fe3Hg4Zs4+MktrWAZfYtyDPAW4GQqV2jnjNlkuJWfq5ILPZE9lt3OA3ceyPJ7lrHbPnfysfMWM206XPS9YdYO/o33n7EnO+0aWLqsEnJn9NSaYSpBnwMcBFwEfGz0YvrHY4y/6thv+sQ7lrDskE1TbdbYmfkBh25hxqzK2UsIkRDggTubfmtQ5yx8FuNDLrWuGnSAe/68jF/8cD7Pey3MmgOHHDXMoccGfv49Q65OmwfMZsez9Xlt/dQfXDWfOQuGecLRj061aeMfGrronF15/p4H8obj92PR4iGOOXlDI7vVzIU/AHwA2JvJ11JKrasGfcW90NMDS/YGIqy8b1/2OgCW3/OIIVeX1J6trxydWz+i6Z+0YU0PX7mwj9d98KFGNm/8yv95l6zk3ItX8tsbZ3PrDbOZPnPKecbRs/C3Upn/nkvj8e4BjgshLG54fCq3d1x2GEc9e/w7vN5pDzF7XuUYGhnpBWBh3xo2bpjNxg3btx349yeF57y3rQv4KpWjaKxl1bPylwOnhBDuAy6MMV7R0G/57P/q48RT17LbPkONbN54zAF6p8ERx2/kR1cv4JuXLeLUc9dMscfFVObDm71twHwqf6tJjbn5+3Dg4eNfGx6CjY+Mf23tqkXMnA2rlu+57bVf/PBfExihyqt6tn4gcDkwdcxvu2Umv/vZHD55/d2N/pLmYl41MgR/u7uRM5nHAm8Czhz9utH5o7XAM2KMt7QwOpVQGBh8CfC/x704f6eZjAzvw4r7YMleldf+dtcG9nncFvY8YHDbdh/65utjf9+1CQ5XORZCOBm4Emh0Vd96KhdML6MS86n9+vo5DD44nVccuj8Amzf2MDICrz92JpfdcM9Eu0x9xrxqeS8/uGo+j64PDA/BTdfM4WfXLODwY6eckI8x3hZjPAdYDLwW+BWwEdja0P8gqVVbNs9g/ep9OOIE+M8rYOuWYe64FW65bh7POPWRKfeX2rOFyq0mrgdOBxbHGP8pxjhhiHfwvNeu4TM338ml193NpdfdzTNfsobDj32ED3zt/nq7TH1mHgIMfHERl797CXEEdtl9iDPfs5LjX9jQBVCAGOMm4CvAV0IIB9Ha2brUmLHryM/7+B1c8LoDeOtJvcxbOMRp75jGnPl7MbT1TqZNb2guUmrCuLPwhuNda/bcyOy525doz5o7wvSZI+y8pO6y7aljvvOSYS7+/n0tDWgCMcbbgHNCCO8AXgC8HXjc6Fimd+r3qKRqPxAUArzxAujpHWbpfncSR8K4degGXe3bQuXa4C+AC4H/ijF29rg6632rptoktfuZxxg3xRi/EmM8ksqNZi4HNoz+M2vSnaWJTBTyWmPXoTd7Lxdpu1lUzsLXUFms8dgY43Exxm93POQNqo15Ox9rbnnfCebWv03lg0VSY3574x5ThryqNuh/uLkvgRGqOH4HfIftc+HntzydUtGR7tbGvJ0LQw3Podcz5mz9JTHGKd9WSDD6hKCvXXoh0PgnO8cG/Yp//o8kHhKtYogx3htjPKWDZ+HtdHfbvrUxv5nK3E+zInBTGwOSWrLtCUG33QJ77P+npj7ZGXoiS5f9hjv/AF1+SLQ0iZ+1uN8gcFv1i3Exj/19g8CHW/ihF8X+vhUtDkhqybhHva1f3UsIH2nyRwzT0/NehrZ09SHR0mRif9+9wKVN7rYV+JfY37ft5DvEuON0TRgYXMr25yhOZh1wY+zve6DJgUhtqffMztFj9xjGfqDjY2+9gmnT1/HGD79tzI94GLh+9ASGEMJMKuuCAfaMMXpMK1FhYHBvtj8DtJ5I5Yz8+tjft3rc/hPFXMqyZh++HEKIwIoY46T3gzboyrPUliZKrWg25M0Y90xRp1yUM8ZcudHNkFcZdOWVMVcuJBHyKoOuPDLmyrwkQ15l0JU3xlyZlkbIqwy68sSYK7PSDHmVQVdeGHNlUhZCXmXQlQfGXJmTpZBXGXRlnTFXpmQx5FUGXVlmzJUZWQ55lUFXVhlzZUIeQl5l0JVFxlypy1PIqwy6ssaYK1V5DHmVQVeWGHOlJs8hrzLoygpjrlQUIeRVBl1ZYMyVuCKFvMqgK23GXIkqYsirDLrSZMyVmCKHvMqgKy3GXIkoQ8irDLrSYMzVdWUKeZVBV9KMubqqjCGvMuhKkjFX15Q55FUGXUkx5uoKQ76dQVcSjLk6zpDvyKCr24y5OsqQ12fQ1U3GXB1jyKdm0NUtxlwdYcgbZ9DVDcZcbTPkzTPo6jRjrrYY8tYZdHWSMVfLDHn7DLo6xZirJYa8cwy6OsGYq2mGvPMMutplzNUUQ949Bl3tMOZqmCHvPoOuVhlzNcSQJ8egqxXGXFMy5Mkz6GqWMdekDHl6DLqaYcxVlyFPn0FXo4y5JmTIs8OgqxHGXDsw5Nlj0DUVY65xDHl2GXRNxphrG0OefQZd9RhzAYY8Twy6JmLMZchzyKCrljEvOUOeXwZdYxnzEjPk+WfQVWXMS8qQF4dBFxjzUjLkxWPQZcxLxpAXl0EvN2NeIoa8+Ax6eRnzkjDk5WHQy8mYl4AhLx+DXj7GvOAMeXkZ9HIx5gVmyGXQy8OYF5QhV5VBLwdjXkCGXLUMevEZ84Ix5KrHoBebMS8QQ66pGPTiMuYFYcjVKINeTMa8AAy5mmXQi8eY55whV6sMerEY8xwz5GqXQS8OY55ThlydYtCLwZjnkCFXpxn0/DPmOWPI1S0GPd+MeY4YcnWbQc8vY54ThlxJMej5ZMxzwJAraQY9f4x5xhlypcWg54sxzzBDrrQZ9Pww5hllyJUVBj0fjHkGGXJljUHPPmOeMYZcWWXQs82YZ4ghV9YZ9Owy5hlhyJUXBj2bjHkGGHLljUHPHmOeMkOuvDLo2WLMU2TIlXcGPTuMeUoMuYrCoGeDMU+BIVfRGPT0GfOEGXIVlUFPlzFPkCFX0Rn09BjzhBhylYVBT4cxT4AhV9kY9OQZ8y4z5Corg54sY95FhlxlZ9CTY8y7xJBLFQY9Gca8Cwy5NJ5B7z5j3mGGXJqYQe8uY95BhlyanEHvHmPeIYZcaoxB7w5j3gGGXGqOQe88Y94mQy61xqB3ljFvgyGX2mPQO8eYt8iQS51h0DvDmLfAkEudZdDbZ8ybZMil7jDo7THmTTDkUncZ9NYZ8wYZcikZBr01xrwBhlxKlkFvnjGfgiGX0mHQm2PMJ2HIpXQZ9MYZ8zoMuZQNBr0xxnwChlzKFoM+NWNew5BL2WTQJ2fMxzDkUrYZ9PqM+ShDLuWDQZ+YMceQS3lj0HdU+pgbcimfDPp4pY65IZfyzaBvV9qYG3KpGAx6RSljbsilYjHoJYy5IZeKqexBL1XMDblUbGUOemlibsilcihr0EsRc0MulUsZg174mBtyqZzKFvRCx9yQS+VWpqAXNuaGXBKUJ+iFjLkhlzRWGYJeuJgbckkTKXrQCxVzQy5pMkUOemFibsglNaKoQS9EzA25pGYUMei5j7khl9SKogU91zE35JLaUaSg5zbmhlxSJxQl6LmMuSGX1ElFCHruYm7IJXVD3oOeq5gbckndlOeg5ybmhlxSEvIa9FzE3JBLSlIeg556zEMIu4cQ7gshvKDO9w25pMQ1GvQQwtkhhL+GEOYlN7odpR5z4E3AEuDKEMLzx37DkEtK01RBDyGcDVwELAVOT3h444QYY3q/PISZwEpgwehLG4F/jDF+y5CrU0IIEVgRY9wt7bEon0ZbtWn0yz1jjA+MCfmc0dfvBfaNKUU17TPzU2rGMBu4KoTwPzHkkjJigjP0dzM+5AA7AycmPbaq1M7MQwgB+BNw0CSbGXI1LIRwNPAdINR8a6fRf6+ueX0IODbGeFu3x6ZiqDlDn8i1McZnJDWesdI8M38KsOck398I/ENCY1ExrAXmUYn32H+qal9fBKxJeIzKtzOBRyf5/tEhhP0SGss4acb8XYx/i1KrOuXy/Em2kbaJMf4B+HMTu/x3jHFFt8ajYplgjnwivcBbkxnReKnEPISwFDiJHd8O1zLoatZHgQ0NbLeByh9MaUoNhhxgOvCqEMLc7o9qvLTOzN/I1CGvmg18PYSwfxfHo+L4Go0dW48A13Z5LCqAEMJRwGVMHfKxEl+mmHjMRy8gvAmY2cDmG4D1VM62HuzmuFQMMcaNwJeoXNysZyPwcS+uq0G3AZ+iMlf+SAPbzwXOH13kkZg0zsxrlyPWilT+D7sdOAfYNcb4ztE/pFIjLgW2TvL9HuCzCY1FORdjXB1jfAOwG/BO4AGmnspLfJlioksTp1iOuJlKyH8EfAj4WVqL75V/IYRbgUPrfHsgxvicJMej4ggh9ADPBP4JeDIwjcpcea1ElykmfWY+0XLE6lTKx4ADY4zPjTHeaMjVpnoXQr3wqbbEGEdijD+IMZ5A5YThM0w8BZPoMsWkz8y/xfa1449SmQf/IPB/Y4yTLcSXmhJCmA08RGX+cqwVwFLny9VJIYT5wBlUllwvpPJ5h63Ap2KMb05iDImdmYcQdgP+BzAMXENlaeJBMcYvGHJ1Wp0LoV74VFfEGNfHGC8F9gZeBPyEyprzV4cQmlkF07LEzsxDCL3Aa4DvxhjvT+SXqtRGb9b2CyrLW6FyXWYfPyikJIQQDgSOBT6fxLRxqndNlLqt5kKoFz5VWGnfNVHqtuqFUC98qtA8M1ehjV4IHaSyYsoLnyosYy5JBeA0iyQVwLTaF8LA4GzgNOBpbH+cWz3rgZ8CV8b+vsnu8St1XRgYnE/l2D2ayjrfyawBbgCuiv19m7s9NmkyYWBwMfBy4O/YvvqqngmP3R1iTuXuYE9tYhxHAcePDkRKRRgYrN5v5bAmdjuGyqeSX9eVQUkNCAOD84ArgX2a2O0YKu09u/rCuGmWMDB4MM2FvOpJYWCw3n0wpCT8Hc2FvOrpYWAwlSfDSKOeRXMhrzohDAwuq35RO2f+2DYG9Lg29pXa1c7x57GrNHXk2K2N+UR3/mrUjDb2ldrlsau8aufY3bbvRHPm471g7wPHfb1lc+BZL13DWz62so0BSN23ZVPg4jfvyu9vmsuGdb0s2WsLp79rkGOe28gDBqR0PXDnNC552xLuuHU206ZHnnLSes65cCXTJm7/1DH/5r23b/vvRzcETjv4AI57/vqODVjqlqEh6Fs6xL9++15232eIG787l4++cSn7HnwXeyyb7ElEUvouedsSFu4yzJV/+CvrV/dw/gv34hufXMSp566ZaPPm1plf99X5zN95iCOO96k/yr458yJnvW8VeywboqcXjn3eIyxeuoXbbpmV9tCkKT10/3SOe956Zs6O9C0d5vDjHuHe2+o+brO5mF/71QUc//x1BD9rpBxa9bdelt87g/0O3pL2UKQpnfyq1fzkm/PZ+EhgxX3T+M31cznyGXWnCBuv8t/unsafb5nDSaev68hApSRt3QIfes3uHPe8dex3iDFX9h1+7Ebuu30mp+x/IK984jKWPX4TJ7yg7rNHG4/597+8gMccvpE995/sQblS9owMwwfP2p1p0yPn/pv3Mlf2jQzDv7x0T4569nq+cfftXPXHO9iwtofLzl9cb5fGY/6Tby7kxFPWdmSgUlLiCHz4dbuxdnAa77vyQaa7ClE5sHZVLw+vmMYL37CGGbMiixaP8KyXruPXP6l9DOI2jcX81htmsXrlNE481VUsypcLz1nC/X+dwfuvvp9Zc7xFqPJhp12H6dtjK9+6fBFDW2Hdwz386OoF7H1Q3fsITb00EeCHX1nIk565nrkL/MOg/Hjwrmlce/VCps2InHbIAdteP/sDy3n26Z6YKNve/dkHufzdu/LtT+9MT0/kcU/eyBsuqPv5nsZi/vZPOs+o/Fm63xDXPPSXtIchteSgIzdz0ffua3Rz17OX18EAAADhSURBVBhKUgHUxrydaRSnYJRXHrtKU0e6WxvzduYRXX+uNLVz/Dl/rjS1c/xt27c25jcBwy38wBHgZ20MSGrXT1vcbwtwcycHIjXphhb32wL8v+oX42Ie+/vWAO+luaAPA/8n9vcNtjggqW2xv+9B4AKaO3aHgHfH/j7voqjUxP6+XwJfbHK3IeA9sb9v2ydCQ4w7TteEgcFFVB5JtAAI9cZA5RT/57G/b3WTA5G6IgwM7kLl2J3L5MfuWuCm2N/n9KAyIQwMLqXyxKxZtHDsThhzSVK+uDRRkgrAmEtSAfx/3Sy8hI22nQYAAAAASUVORK5CYII=\n", "text/plain": [ "