diff --git a/doc/notebooks/demo_central_moments.ipynb b/doc/notebooks/demo_central_moments.ipynb index c5997429f0eb6903528875a5ef4a2f3eb5d68726..1d59c9f947aee462c621fd10eec0ec556da76812 100644 --- a/doc/notebooks/demo_central_moments.ipynb +++ b/doc/notebooks/demo_central_moments.ipynb @@ -12,11 +12,9 @@ "\n", "from lbmpy.forcemodels import Guo\n", "from lbmpy.methods.abstractlbmethod import LbmCollisionRule\n", - "\n", - "from collections import OrderedDict\n", "from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments\n", "\n", - "from functions import get_shift_matrix, get_central_moments" + "from collections import OrderedDict" ] }, { @@ -257,7 +255,7 @@ " " ], "text/plain": [ - "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f9b5d1abfd0>" + "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f7cb10fb1c0>" ] }, "execution_count": 3, @@ -349,7 +347,7 @@ " " ], "text/plain": [ - "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f9b5d2683d0>" + "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f7cb0f67fd0>" ] }, "execution_count": 4, @@ -376,11 +374,12 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 5, "metadata": {}, "outputs": [ { "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAADhCAYAAACEL7KHAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2dUbLUNre2NxTXKQJVGUCYAYQRADOA/CMgzCAprvjuqDADwgg4YQaQEfDBDJIBpIp9qEzg/O/jbXUst91tt2W3Lb2qcsuSbVnrWcvqZUm2r7148eLuxcXFJy1d4d1//vOfJ10bnGcCJmACJmACJmACJnCcgHypP7XX9117atu1G40Nr7TOzs3wVzPhdRMwARMwARMwARMwgdEEfu044pHyHpPfdMZeyzuz89VBy1kmYAImYAImYAImcCoB+Ve/tY9VHll7zlh7v7OkVTmGTX/Xck/rX89SiUxOKn7BE/8ike5o+VV5drgT69echwEthVMpch7Tujn0EzKbmE0uPKbIcT1Gcp6UBLip5Xctr1WDN1o6x1XPU7ttnlUsmQf4UfEvWhiC/kXLe62bbUKVmvMwmKVwKkXOY1o3h35CZhOzyYXHVDnW4ox9lSBPtDyTmt7GqnJqLAFx/EnH4OC+C8dqnV5G0ji8DgkImPMwiKVwKkXOY1o3h35CZhOzyYVHCjlW4YzF6nEqAQGegP3cUc5H5T2U4dzs2Oas8QTMeRizUjiVIucxrZtDPyGzidnkwmOyHHbGYsPIJfVQglx2CBPmi7HdYToBcx7GsBROpch5TOvm0E/IbGI2ufCYLIedsdgwNp8a2Ot1a/OCnlkAcx6mgFI4lSLnMa2bQz8hs4nZ5MIjlRx2xmL7yCEVHK1DT6J6mHK6ps15GMNSOJUi5zGtm0M/IbOJ2eTCI4kcdsZi4ygldbsUQc8spzkPU0ApnEqR85jWzaGfkNnEbHLhcVQOO2Ox4nNIdc0VC3IFD573jjlMI2DOw/iVwqkUOY9p3Rz6CZlNzCYXHknksDMWG8fmUxq/DsOTXUORIS9M5N+8vOcSwJyHkS+FUylyHtO6OfQTMpuYTS48UslhZyy2j1xSHyTI9x3ChJ4xtjtMJ2DOwxiWwqkUOY9p3Rz6CZlNzCYXHpPlsDMWG0YuKT4n9UOHMPeU97nhyXfs4qwRBMx5GKxSOJUi5zGtm0M/IbOJ2eTCY7Ica3TGwkS30IsTq86powTkbPFB0kvF1QdIOUDrDFH+qOUpaYfpBMx5GMNSOJUi5zGtm0M/IbOJ2eTCI4UcN2I050tJGDxLAi9PI/CtSuY28T3Fva+dV3v45xABesH4MPh9xUzYJ36g9GfFDukImPMwlqVwKkXOY1o3h35CZhOzyYXHJDmuvXjx4q648FHpO7XzE2NyygRMwARMwARMwARMICkB+Vx8R/q14mtrHKZMKqwLMwETMAETMAETMIE1E7AztmbtuG4mYAImYAImYALZE7Azlr2KLaAJmIAJmIAJmMCaCdgZW7N2XDcTMAETMAETMIHsCdgZy17FFtAETMAETMAETGDNBOyMrVk7rpsJmIAJmIAJmED2BOyMZa9iC2gCJmACJmACJrBmAnbG1qwd180ETMAETMAETCB7AnbGslexBTQBEzABEzABE1gzATtja9aO62YCJmACJmACJpA9gdV8mzJ70mcSUJ9Z4HNXfPfznta/nqkaWZ9WXH+tBeQboHe08E1Qvqvq0CBQCqdS5GyotnPVHDqxVJlmE7PJhccUOdwzFttEFikZxE0tfGj9tQR6o+X7LARboRBizHddPyr+Rcsrrf+ihY/bm3lDX6VwKkXOhmo7V82hE0uVaTYxm1x4TJXDzlhsF1mkZBRftTzR8kwCvc1CqBUKIb585BXH912ontbpfSSNI+wgAqVwKkXOY0ZtDv2EzCZmkwuPFHLYGYttwykTGEPgiXb+3HHAR+U91AV6s2NbiVmlcCpFzmM2bA79hMwmZpMLj8ly2BmLDcMpExhD4KF2vuw4IMwXY7uDHFNBKIFTKXIes2lz6CdkNjGbXHhMluNGzOW8qbongaEfJkF/Uvq3Zo2UZiL6U8WeiN4E4/XFCdS2euy8t47tkPv2UjiVIucxezWHfkJmE7PJhUcqOdbWM/ZcgjEJ+r2W8IRapUHl46Q9VmxHrCLinzMTCI7WIXv0MOXFRSmcSpHz2GVnDv2EzCZmkwuPJHKspmdMTtZd6Ym5NoRHWtrDGuTtzc/RccFp47UCt5XmaTYHE1gDgdtrqMQG6lAKp1LkPGZy5tBPyGxiNrnwOCoHztg3tewhjlGMSMkRoifgDy1jegR46g8n66865ow/annJSiMwJhvlaX96y74opjftQjGTpnmtAI6bgwnMSeDyQOHhTokbhNJDKZxKkfOYPZtDPyGzidnkwmOKHN8FJEl7xuQEMWRzLxQ+Jq6PvVCM04Uzt5svprzv67wPrTLpFXsQ8rTfBy3VO54Uh0nUYbNjE0hGQPbF60Mor+vGI+QVb4OlcCpFzmMXkDn0EzKbmE0uPFLJwZyxf2pEIY6JLZ/iEVF6yXDsQnjMivI+hwytBwet/YfHcQx5OpjA3AS4OcAO2yH0jLVvHtr7lZIuhVMpch6zW3PoJ2Q2MZtceJwqx98BB87Y2gJ/bjunq64cw44Ii0P2kxZ6HkLvA9nNQJdh+DNs5nvdBFIT4OneHzoKpXf4s+y0eUPRsVsxWaVwKkXOY4ZrDv2EzCZmkwuPyXKs0RmLerr0h0YvF394wUG7c+RPDkesz1GLzaCMVJg4aAc1sb5lhwylXyquem4pXuvYHnMen5J2qJgUwcn2cGXt5tB/1ZtNzCYXHinkSDpnLMZ8coqnId9IOOaDMQEa54yeBj6+/LPit1oIfb0O/BlGDl21d2E/YoWnTmAOHoFvVcKFOXW7+XjVFv9MIRBs874KwV6JH4hxuHmYUnZOx5bCqRQ5j9mmOfQTMpuYTS48Jslx7cWLF/Q8fdJCj9OmnBjV939V7+iPT3n/t0VZVGcHEzABEzABEzCBQgjIX+GNEK8VX1vjMOUYNfCqi9DzcyGBWH+neFNO5RiBva8JmIAJmIAJmEBeBNY4TDmYsJyuV1oYvgxDmnxGyXN1BhP0jiZgAiZgAiZgAucmsGlnDHhyxPzG/XNbkc9vAiZgAiZgAiZwMoGtD1OeLLgPNAETMAETMAETMIE1ELAztgYtuA4mYAImYAImYALFErAzVqzqLbgJmIAJmIAJmMAaCOCM8eThMy28ud7BBEzABEzABEzABExgfgJ8WQj/6wJn7Hstr7X4De2C4GACJmACJmACJmACCxDgdVz4X5UztsD5fAoTMAETMAETMAETMIEuAp4z1kXFeSZgAiZgAiZgAiawEIHNv2dsIU6bPE39MlzqzjcTeSEuL8j11wkgkjiIK58V43ug97Te993UxGfdVnGl2GMpch6zPnPoJ2Q2MZtceEyRwz1jsU1kk5JR8L3Rj4p/0fJK67wcl4+EM0fQIQEBsbyphQ+wM+b/RovZ9nAVoyLssRQ5e9S8yzaHHYq9FbOJkeTCY6ocdsZiu8giJaPg46M4Cu+CQFqnt4Z0NVkw5Ds+nQBMtTzRwtMwb08vKe8jxacIeyxFzmPWag79hMwmZpMLjxRy2BmLbSOX1BMJ8rlDmI/KeyjDudmxzVkmMBeBUuyxFDmP2Yk59BMym5hNLjwmy2FnLDaMXFI8Ltv13rgwX4ztDiawFIFS7LEUOY/ZjTn0EzKbmE0uPCbLcSPmct5U3WPDkAaTzT8p/VuzRkozQfqpYk+QboJprNcMGzmdq36nXCcWZ6YmUIo9liLnMfswh35CZhOzyYVHKjnW1jP2XIIx2fy9ll+bqlM+TtpjxXbEmmD214OjdYiThyn3uTlnHgKl2GMpch6zEnPoJ2Q2MZtceCSRYzU9Y3Ky7kpPzGkiPNLSHmYjb28elI7DsfhRyzOt31PscJzA7eO7eA8TWIxAKfZYipzHDMcc+gmZTcwmFx5H5cAZ+6aWPcQxihGp2jH6Q4eM6XnhaTScrL/qmDPiXL1kpREYk43ytD8OHPmEMee8OiLP38sDYgUPnveOOZjAEgRKscdS5DxmM+bQT8hsYja58Jgix3cBCc7YP3UixGHb6FjO0VcddG/0gTqgPpYY5wrHajdfTHnf13l8VHMXlP9Zic+KH+8yC18RC163AIUu5zTkhYn8hdOy+HMTKMUeS5HzmL2YQz8hs4nZ5MJjohx/ByprmzNGvXhElF4yHLsQKmdLeThfDscJ4LTiwLZD6BmLnNr2Tk6bQGICpdhjKXIeMw9z6CdkNjGbXHhMlmONzhhORNvpYr4Ywl7IIftJS+jhIcthnwBPnf6wn131WtKT2HR0O3ZzlgkkJVCKPZYi5zHjMId+QmYTs8mFx2Q51uiMRUNochzuSnc4FsFBu2NnIrbmdkp8GOK9VLwbvtU6Dixz8Z6293c6CYEwQTP0PiYpNIdCSrHHUuQ8ZpPm0E/IbGI2ufBIIQdzxtYW+IbiGwnHqy2YaI5zxjw0PnL9s2J/dkYQBoTA7L72hSPxAzEMTu2AIrzLMQLiyR0RgbmOBL5Vic3yHdDdvMdqS9k/pdhjKXIes2Zz6CdkNjGbXHhMkuPaixcv6HniI770OEW9UjGv9aZUb3qAcNZ4WayDCZiACZiACZiACayagHwW3p/6WvG1NQ5TrhqeK2cCJmACJmACJmACKQmscZhysHzyJpns/0wLQ0TfK10NbSrmLf4OJmACJmACJmACJrB6Alt3xhhWZY6ZgwmYgAmYgAmYgAlskoCHKTepNlfaBEzABEzABEwgFwJ2xnLRpOUwARMwARMwARPYJAE7Y5tUmyttAiZgAiZgAiaQCwGcMeZdMQn+MhehLIcJmIAJmIAJmIAJrJwAXxbC/7pgAv/3Wl5rIXN1n8nRk5G8Of65FgJ1JTxV/urqelU1/5qACZiACZiACZjAUQK8CQL/67ctPE3Jy1wrzxGxtE7Fq5fUknYwARMwARMwARMwgS0T2MKcMT4MjvcYAu8S451ifDnAwQRMwARMwARMwAQ2TWALPWP0iv1305QXrLycVIZyw7vX+MA6cwF/Ub6/STmTHsSWGwQC3wDlk1z05m7y02IIMWcQF26i+J7nPa1nOdXA9nBlQebQfyWZTcwmFx5T5NjctylrYR8r9ncoY3u+EJNq/p/iR2GT1nEU+MD6I60zL9AhIQExZcj8peJ3FKuYOY7kwdsO2b9M3miVGwNuEHDIvhWf7JwxyWR7kHLNQRB6gtnEYHLhcYocOmab36ZUxWnE+Sg4X0d32CeA47WbX8dmMaOXjD89eiMcEhKoL6SbiitHjKK1DmvSzG10EAGYaHmiBdt8mysUyUfDWrw9mEO/hZtNzCYXHink2MKcsUp7EpZeH5yNbIc3YjM9KcXcuj/Fit6ZZqBHjD8JGDqkI/BERXUN/35U/sMOPaQ7s0taIwHbw5VWzKHfOs0mZpMLj8lybMIZq50I5j0x9MNdNhP47VjERk0Kp+svGO1vqnLaTlrPbs4eSADnl6G3dgjDk80HT9r7OJ0fAdvDlU7Nod+2zSZmkwuPyXLciLnMn5KjgENAdz5zvj4p/VvzrEoznLZ7j5jSOF0M+eCMhScoGe4Ik9S16gAB8cE77woVN22venEUj9JBV4Gl59UMj2G4dWwHb8+DgO3hSo/m0G/PZhOzyYVHKjnO0TP2XJV/JbW818Kw4y4oHyeNyfnNnh0mxOJ1EoeF110099Emhy4C4oQjhkPbdF7H6qCr6NLzgqN1yA5xeh3KIGB7uNKzOfTbu9nEbHLhkUSORXvGaseA+TQEnvhrD/GQF83B0THfsrPDyQToaXwnjjjAF6fooD4uOM5flL6tcprOHbs47BO4vZ/lnIIJ2B6ulG8O/ReB2cRscuFxVA6csW9q2UMco0ibYj5TcLZ+VNEvW8XTA9bOa+2Sb1Js6En5Q8uYHhWeUgtMIzjKZ3gX5s3hy9E60PH0WH5RHBw6Jqe/17J7hUZ04jISlwfEDHdKOK4OZRCwPVzp2Rz67d1sYja58Jgix3cBCc7YP3UixGFb8lh/3tWQjmKcLhyO3Xwx5TGURt7eu7C0rYhemZrPPTGYHFQWDtQtxZHDVJ/jQvEYHcD/QaiUjv2gBWeMBynCZPWwuYhYcvMgCbJ2Oc4hr0g2RRhAS0jbwxUQc2gZRiNpNg0YWs2Fx0Q5/g5UzjFnjHPTU9N+6o/3h6Ggz8QhKB16ZZjAX801U8x8M4ceAuIDyzuKdz1iWm8/gTpIBxynsnAu2o7FV+Xd1VJy4MYBPu0Qesb2bizaOzqdFQHbw5U6zaHfrM0mZpMLj8lynMsZ4w8scrqUpgcHgS7kADBBP/Qu0CtT5dfbWGeYrOtPkF2KDuKCg3RfcXtOFw5aszt1qA6CHtpcKSs4He1tpaSZj/dDh7D0bn6WDnBYHcohYHu40rU59Nu82cRscuExWY7rMZfFUlEvS+1A8KcWHDR6dar3iSkPZyDaX2n+5ErvlRGCOIgZDhZGwQteXzcX5T2DaeOIiKm2wXNPB43926s4Yuim2CBmDLNfKsbRrYLWYcJ8yKdXOf5tEQgTWbNz5G0PV5o2h5bFN5Jm04Ch1Vx4pJCDOWPnCPTavJEA9HoxyRnHgN4EPrDMdxTfaiH0/dm7V+aKT/uX4VscMoZ22yE4uiF/qA6aDlw4lhjdRA5dc2NB68Fu70tmbJn4gey4zbsgJPuiigc3CQTmKhJ+Vx72w9zD3dzRasu2f2wPV/ozh347NpuYTS48Jsmx6g+Fq5Gmt+aTluijwsr/X+Xxcebq6b5Yr06lJlDzjhwM5f2fzkMPph2y1MBdngmYgAmYQPYE9P9JxwmjWNfONUw5FLJ7ZYaSmnc/XjcSejQuZDis8+4yO2LzcnfpJmACJmACBRBYtTNW/9njkDH01g4eBmoTmSktPdADyYtewzAyT2J6TtRMvF2sCZiACZhAWQTONWdsDOXQK1M5X3II3Cszhl6ifcW9/XRmopJdjAmYgAmYgAmUTWD1zhi9MnWPTJjszwfG3StTtt1aehMwARMwARPIhsDqnTFIu1cmG3uzICZgAiZgAiZgAi0Cq54z1qqrkyZgAiZgAiZgAiaQHQGcMZ6Ie6blMjvpLJAJmIAJmIAJmIAJrJMAXxTC/7pgmJInFV9rIbPvVRLadJ6gIUpeLvq8Pnt4qvKp8ldX1/MQ8llNwARMwARMwAQ2SIAHEvG/ftvCnDFep1B5joDWOhXnRbBM5HcwARMwARMwARMwgU0T2MKcsZ/kgOE9hsBTld8rj7fzO5iACZiACZiACZjApglsoWeMXrH/bprymSovhxXHlcA3E+lJpJfRb82HyIQghgyXh/eu8XF15lv+ony/iPgA15LsUbJys8j3OO9pvcgpFSXp+4DZd24ymxhLLjymyLHqb1PG6rpK1cI+Vuxhyi5AdZ74MJTL9zvfkaWYuXfkPdK6HTKgnBDErppjqfhROFzrOL184B62zL10aBEQl+ztUTJyjb3RgnOOk45DFn1XV+kiQgn6PlWRZhOTy4XHKXLomM18mzLSmipO4/ZYC19Hd+ghUCv4puLKEWM3rXN3Tpo5dw6nE8Dx2s1hpBixpZcMvvSEOLQIiA8NTvb2yDWm5YkW7ONtC0MxyVL0fYpCzSamlguPFHJsYc5YpT0JS48Ef4TFdvvHZnwwxbcju4bMPir/oVhyB+9wGgHmL/7ZwZAeMRwO7NQhJmB7jHnknrK++zVsNjGbXHhMlmMTzlj9B8ecHIaBuPtkAr//9GKjbqZwGBgqaYcwPNl8IKK9j9OHCeB0/YUd9uxmR3cfjO1xn0nOOdZ3v3bNJmaTC4/JctyIucyf0p8Yf1YMWzDn65PSvzXPqjRDPbv3iCmN08XQGs5YeIKSYYAwgVqrDoGAGA1xBm6F/R2PIyC+3AF1hco2tb3qkaz1MNjOuwrMIc/2mIMWh8tgffezMpuYTS48Uslxjp6x56r8K6nlvRaGHXdB+fx5MTm/2evAxF+8TuKw8LqL5j7a5FATCI7WIT5DHDYDHUhAtogjxk1D8wZhrJ0PPNvmdrM9bk5lkypsfffjM5uYTS48ksixaM9Y/afFvCUCT6O1h9LIi+Y66Zhv2dkhKYHbSUtzYfTmvpOtcpNxcYqd18fhJP+o5ZnKKOkhFdsjBlBOsL77dW02MZtceByVA2fsm1r2EMco0qaYaxOcLf50XraKpwesndfaxckjBC4PbA8e/JcD+2S9SfaHw/OHljG9gzwhF+w24qN8htCx6+bw5Wg71/F3VQ72TxhTt6sj1vtre1yvbuaomfXdT9VsYja58Jgix3cBCc7YP3UixGFb8lh/OF8pVDF/Ovzh7OaLKY9hHvL23tOkbeSX2GMgsccFseIBBw6CWTuEvDCRv709+zR8JOS9FIKqLIbVbynevXOMcutzEA+2c+37WYd+VsyrW7IJksf2mI02jwtiffczMpuYTS48Jsrxd6CCM3aOQC9C+4m06k9IgvGntAtK59pjsJNxhhUcWpzbdrhVZ+w5vO0dnT5MQHaJvd5RvOsR03rFXHFwdgfb+eGzbX6r7XHzKhwlgPXdj8tsYja58Jgsx/WYy2Ip/rQip0tpehcQ6EJ/ZkzQr3pxFNNbwFyc8AfHLg6HCTCH6YeOXegRgie9Qw4nEhA/bhDuK25O2Kc0HLRml/VgO+fgjIPtMWPldohmfXdAqbPMJmaTC4/JcpzLGYscq/rPDechOGj0ONhhiI12cErsGP69VIxzUAWt49wy1Pv0Kse/pxAQRxwsLjxe8Pq6uSiPifdNu7WdC4qYlGiPYcJu6I0WiTJCofoepFyziTHlwiOFHOcapqRH4Y0E4NUWTCbnT4teGz5k/bPiYj8lItlThcDzvgqEMfED8Q0Ob6rzlFYOr2TBIWO+WDu02drO/yVUhD3q+sJRJzBfkPC78mjf3ivezZGttuT9U4S+T1Sh2cTgcuExSY7NfChcDRm9PDhr/kB4bMhOZUTAdp6RMi2KCZiACRwgoPaem3pGWK6da5jyQPW8yQRMwARMwARMwATKIXCuYcpyCFtSExhAQHdGDH3ymS+Gt/j2ajWEr7h6keyAIryLCZiACZjARgnYGduo4lztvAjI6WJeUfvpzLyEtDQmYAImYAKdBFbvjLnHoFNvzjQBEzABEzABE8iEwBacMfcYZGJsFsMETMAETMAETGCfgCfw7zNxjgmYgAmYgAmYgAksRgBnjJ4nJg433xy+WAV8IhMwARMwARMwARMokABfHcL/usAZ4ymu11qKe1O0ZHYwARMwARMwARMwgXMQ4Ol5/K/KGTtHBXxOEzABEzABEzABEzABEfCcMZuBCZiACZiACZiACZyRwOqfpjwjG5/aBEwgIYH6RbaUyLdS+awZnzdjzmpWoRQ5hyhNLO5qP77XeU/rX4ccU8I+tpFYy7nwmCKHe8Zim3DKBExgBgJqpD6p2I+Kf9HCVwV4wS0fz2bOajahFDkPKUwMbmrhA+nMhXmjJSsdH5J9yDbbSEwpFx5T5bAzFtuFUyZgAokJqJHiY7j8Qb8LRWudXhLS1eTVkL/luBQ5j+kI3Wp5ooWnxN4e27+k7baRWNu58Eghh52x2DacMgETSE/giYr83FHsR+U9VEN2s2PbFrNKkXOLullLnW0jsSZy4TFZDjtjsWE4ZQImkJ4Aj293vccwzBdjew6hFDlz0NW5ZLCNxORz4TFZjhsxl9NS9Z0tQxFMyv2k9G/NkpRmAudTxYMncM5RZrNOpayb4zKaNuduzjWX7o3/5m7+HYelyPmvyrw2loBtJCaWC49UcqTqGXuuCjEp972WX5vIlY+T9ljxYEesPn6OMptVK2XdHJfRtDl3cw6O1qHrP4dhylLk7Nayc4cQsI3ElHLhkUSOyT1jcrLuii9zPwiPtLSHI8jbmy+i44LTxmPut5Xm6aoqTCiTRv1HLc9Uxr2r0sr9NcdldD+Bc+81sEzNV3OW26upybwVKUXOeSnmXbptJNZvLjyOyoEz9k0te4hjFMdTf+nPKDhbOEIvW4cwlhrlaX96y74opjftQjGTeHnMHceNcEqZOIWci5DDnXYliJggyx8jZeJJJnRijhXF4z9n4HzsGjhe6W3scXmgmuGOkhuyrYdS5Ny6ns5Zf9tITD8XHlPk+C4gwRn7p06EOGwbFOtP7Cs7KsYRwnHYzRdT3vd1Hh/DbAZ6BB6EDO33QUv1ziHFOBCjy9Qxn1XeZ8WPQ7k5xDWLe6fIYo7DqS3NWTU7eA0Mr/m694SrFipJ29AOIS9M5G9v30y6FDk3o5AVVtQ2EislFx4T5fg7UEk1Z4zyeLRz50jVJ6gcI1UWR6kKWg8OWrsBxgGjd6sZBpXZPMDrnQTMsRNL8sxBnEdeA8kreYYCuRnjum+H0DPWvllr77eVdClybkUfa6ynbSTWSi48JstBz1iqQGO7c7rqQhl2pJIX+gNiWOZ/tIS7YbKbga6+0DiH/EFlqmwcOYd+AubYzybllkGcdcIx10DK+p2rLJ6mDvPjmnWgx5fe7Fyu31LkbOrQ6+MI2EZiXrnwmCzH9ZjLpFTU06UGll6uH7QEB+3OkUYXR6z9JzW1zEkCZXSwOS6jzKmcu66BZWo+41l03TN14VLxbgqB1rnWmWP6dMZTL1p0KXKOhBomLrdvtEcWk8futpFYj7nwSCFHyp4xnoZ8o0pxB8yEXP6YuPPlY8A/K36rhdB3F0zjHP2ZKT20TMp16Cdgjv1sUm4ZynnMNZCyfucsK7QF91UJ2gfiB2obws3aOeuW8tylyHmQmfRKTwGBucQEvlVJ+87c4N284mpLeT+2kVjnufCYJEcyZ0wXGH8wzJlphyiPC7Let2tIJ2qYh5bZPqHTMQFzjHnMlRrKWfsNvgbmquvS5dZs+FZh1qEUOY8pURyidv/Y/iVtt43E2s6Fx1Q5Ug5TxoQPp3jVRbhjupAQrL9T3O4ZO1yKt5rAdgn4Gtiu7lxzEzABE0hK4NqLFy+Y2/VJC3O6FnOGdK4woZchCz6j9Ivy+oZvtPlw0LH0tHHnjWOHTLzDbPcuM607DCBgjro08u0AABoCSURBVAMgJdol9TWQqFouxgRMwARMYAEC+g/gwcbXiq+dzRlbQE6fwgRMwARMwARMwARWSaDpjJ1rmHKVYFwpEzABEzABEzABE1iagJ2xpYn7fCZgAiZgAiZgAibQIIAz9k2dDnFjs1dNwARMwARMwARMwARmIJDu25QzVO5C46i8c+x5XTYT8wlPlX/yBP+rIvxrAssQsA0vw9lnMQETMIENE9h9mzLZe8YSw+BFsbt3Emn9tcqvnvhMfB4XZwJzEbANz0XW5ZqACZhAZgTWOmfsJzlgu/eQiTmvwfheebyywsEEtkDANrwFLbmOJmACJrACAmt1xugV++8K+LgKJnAqAdvwqeR8nAmYgAkURmCVw5TqAWt/u4w/Nj4hE30uKRddSS7mxfFdQwIfV7/Uwktws5QXIXMPpdlw7vocI590336hNUPWi71Qe0xdU+4rGd2OjQAqXoz08A3Pe1ovej605M/impkixyqdsaY91wb7GINt5ueyLvlowHgD76MgU63QT+Rp+RDyHW+TgHRIo5utDW9TK/PUWrpmbutLxe84g2IeRgrXcrYOmeR0O4bCj4TaHt5oN264ufGGW9FBTLK4ZqbKcX3NVlBf4HjMOd85IN/uYQX0IbnpJeNOibsmhw0TKMSGN6yhdFWXrvm0yU3FlSNGyVrnOibNQ0g5B7djA7SLPWh5ooU2/+2AQ7LeRRyyuGZSyLFaZ0zCccfAUB29QxgwE/hzvIvgQYU/JRt30M1AjxgNe44yN+XMdr3WXQk2nK0ORwr2RPt3TS34qPyHHdf4yOJXvbvbsVWrZ7WVy+WamSzHKocp6z8x7iT5IwtPUHInEeZVrdayTqgYTtddyckddFdoO2ld+zhvZQQKs+GV0T9bdXBI2vNdqUwYnmT7rteMDRkFt2MZKXNBUXK5ZibLcWMJ6PpjwqGgO/KOFuZPRA2W0gzHNV/qyhgyxxDvgvbbDeedUOaunDWtSA486q5QOaHaXt1pzyXvXOV2CbTlvBM4HbXhLfM4pe4nMDzlNGc5ppbt2LlvHdthq9sl/1nbsa1yK7neuVwzqeS4vpAxPFeFX+lc77Uwt2AXlI+T9ljxrmdI699qudZedgddrYwqs3XsqpOSG0eM4clmT+Bc8s5V7qoZn1C5UZwG2vAJ1dj0IaMYbkzS4Gjt2rGO+nODWUxYuB0rhmtGguZyzSSRY/aesfqCZM4EgScGeYqkGcjrmmfR3CdaP7VMHUdj+KOWZ1pf89OZ9BS+Ux1xYC8U45yNZqjjguP7RcffVrrp3E0pdyscwTc5nMp/8okzKuBUhjouJ1u7nZFKh4iySDs2pCLeZ7MEcrlmjsqBM8Z8Bob/2k7STnt1g/iHMsbc2fHECE5W8/1gOEIvdwVfrTDW2s5r7bKXHF1m/WfAuQhj5Lg6YqFf1ZO5csjX7PY/RV56HL+onODQMYH4vZbdKzTq8wRHeJBudDyO4eo5qo4pw2j+KU++lrKke66bxdqBjdlab/spZuHOmZuiVYaJut2TSeUt2Y7tnd8ZmyCw6WumQXiKHMy1rKZf4YwxHMaFQ+ZXLXtBFxb5J/Uk1cdeKOYPnMZ8N19MeZybPM4dBW3r7dXRtqqeY8rUvjgdnxXzvqdVBtUNB+qW4qbDdKH0aHlVDvweBEFVxgctOGM8lVpNKFY8ulwds3qOQeZU8SmcOLeO67XhVHVbspyaw2LtgM63GVuDjRbUQXvWDiEvTORvbz97eopu25VXWYu2Y+3zO70NAlu/ZgLliXLgF+F//bbUnDHqTU8PPQyVA0CGQuUYKY9GdxeUDr06PE1ZzTVTzHyzdhhcZvvAtaUlHyzuKN71iGm9/TqPQfJynMriD6Dd+MOenq12GFRu+6AC04M5SQdDbbg0jIMZbhAMN5Vce+0Qesb2bjrbO249Lbs/Zzu2dXwl1j+Xa2ayHEs6YzRSkdOlND1AVQPFn5eWcAdJj8Ku4VI+6wyztRu6MWWqiHUGyYWDdF9xNKdLeTRsl41aD5U3cGwcWq1SVvhjaG4bWm7zmBLXx3AaasOlcRzDcGtsmCP1Q0el6U2kV56boWzDCtqxbNlmLFgu18xkOa4vqOSol6a+cGm4goNGr1D1clfl4UxE+ytNQ9bu1Yn26StTx602qM78OaFIXvD6urkojwcNmg34VHlxxGDbDlPLbZeXa3oQp1qnQ204V1Z9cg1i2HfwmvOld6ZgXCrmJqoKWscOmI/59Conz9/a5s/djm0NbpjU3XWDvDVZTqpvLtdMCjmYM7ZUoNfnjSpNjwETWWmUuWPkI7o/K36rhdDlLJDf1asztEyOX2tg+BWHjGGtdgiOasgfKu/XcEAr7nIQ2GVoua3iiksO5TTGhkuDOJThVrmENu2+BKCdI36gNq59LW9Vvr56r6Ed66vbqvJlCzitBOYLEX5XHv+HzOndzamutpTxk8s1M0mOxZwxGRkOAvNF2qErr70P6b1enQRldp1n0TzJwItwB4Wh8mq/MDeva0ho709haLmDKpnxTgk47dlwxrg6RUvAsLPctWTW8lVPR62lTkvUQ3KfvR1bQs4U5xCrof95KU63+jJyuWamyrHkMOVQoxjbqzO03NL243Uh4c7rQobCOu8ui4aJSoOykLy24YVA+zTZE3A7lr2KLSAEFusZG4obZ0ELf2aDenVGlEt53LHilPCUYjVcqrh6D9fQcrayH3IhY5BT9ebOdfK8FZVXFMdT9C1Gs9jwKXXZ8jG2tS1rL03dZQOztGNpaudSTCAdgdU5Y7Vo4W6oGlLTBTm5V4c/SJXNfJVigmROLm+JHE80mOQ2fGI9NnuYbW2zqkta8TnasaQVdGEmkIDAKp0xXXy+G0qgXBdxPgK24fOx95lNwARMYGsEVumMAdF3Q1szJde3TcA23CbitAmYgAmYQBcBJvB/U28Icdd+zjMBEzABEzABEzABE0hH4LtQFD1j/9SJEIdtWcXqpbgpgZ7XQjEJnfBU+X1Pvl3tkdmvOWSmUItjAiZgAiawVQJ/h4qvdpgyVDBhzJOFPE1ZBa3zcc5PWga/H+fqyM3/msPmVWgBTMAETMAEciLAMGUp4Sc5YA8bwvJqC15xcbeRV8KqOZSgZctoAiZgAiawGQIlOWP0iv13M5qZr6LmMB9bl2wCJmACJmACowkUM0ypHrD2N79wSng5597ngUZT3NABJXOQ7MwVDO9e4yP1l6RLs4Fzmas40xtN4JuNTA9gyNxfhIBIBsHX13Al+lrYZyUmjFLx3c57Wt/kXO4pei3GGWuqvlb6Y5TezC9tvSQOkhVH7LXiR0HPWsc5+ESelg8h33F6AuLL/MyXit9RumIeqAns7ZABZcNB+vT1NVB/YuVroWZVtwNvlOTGmBtk7GiTYaper29S6gmVrhsN/oQ3631PEH93aIEc0PnuAQ5AiAG9ZNyBcTfmMBMBcf5JRd9UXDlinEbrcCfNgzQO2yfg62uADn0txJBoB7Q80ULb/Dbeup1UCr0W5YwJGF43w1L0hGAETODfrCd+qqnWMpfGgYc3/pTs9Mg0Az1iOArF2UETwszrT1R+13SAj8p/2KGTmavj4mcg4OtrGFRfC8M4bW2vyXotZpiy/rPlLhwnJDxBiTce5hBtTfkn1bdgDjhddyU/PTJdoe2kde3jvNMI8EfdnrNJSWF4ku27XjM2OGyOgK+vYSrztTCM09b2mqzXG0tIrD9A/ugYqmDSLvNEooZZaYaJRr2A9YQyGaenHsS7oHKioavdho2smMMwRYkTdy5doXLMtb3quTmBZ1eZUd4cZUYnWHGilv1YDW8d22Fr20vTueQ92/W1Fdso9VrYin5OrWcqvV4/tQIjj3uuCr/SMe+1MLdgF5SPk/ZYcV+PxW7f1sqoMlX+t1qutZdWmVtMmsOJWpMt4IgxPNnsHR3Fc+Cp5yhz4KnPvltwtA5d39wk5RZK1nmlywWvr63YTqnXwlb0c2o9k+h19p6x+oJkbgiBJ9l4aqIZyNubT6LjgtP2RdtvK737wzy1zOZJc1g3h8lapEf2nThyo3BxKk8dl9xWVSYOyo9anmk996d+b8M/lzCHHW2UTarry9fCRg3A1d4RONrG3VDDgSN0bXdIx0r9x/CHNo25g+UJCcpuvsuLP5eXrVMw1hrl6Th6y74oDn+STPJ9ryW8lmB0ma1z5pI0hxM1KVti/iD8msMro3nq+OS2qjLpseO6IIy55q6OWNdv++arWbtwR8kN16qCdAD3NbV5s/KZKO9e3VRequvL18IeXWesjMDJbZyuE6ZsVdO2cMYwduZR3dF6mFAbyar8r8o46e68PvZCMX8uNHC7+WLKY4iIvPY7nuhpeKClCtrvgxacMZ5+5A+T+lwoHlMm+3OubHobJI85VBYy7kfccKBuKQ7OfVXAiTyT26rqwU3MZ8W8C2/TAaZakIFrrx1CXme70955yTT11vlW0+bNLfsUedt1U1nJri+V5WuhDdjpVRHg2tFCnUJ71qxfyOts4+pr5bXia0vNGaNy9EDsHKm6ttWfjSrCBVcFrQcHrV15Gkccx2YYVCYHqFyOpZEATgCk1SyCOQxUo+wAm+PGY9cjpvX2K04G8eQ4lYUtJbXVgaJsaTdutmDVDrfqjPbNWHu/rabnsKNVs0h5fa1a0NMrV+q1cDqxbRw5Wa83FpSTxnjndNXnpWeiaoh1EeMo/Y+WPkfpUttC463VKgwqU2XjuWZzhxWEb8Tm0IDRtyobwCG/r3g3/7DeFwdt12Or9UE8td8stlrXKaeIuUNhXl1TLnqe6AH82szMaH0OO1otntTXV6Z2Ueq1sFq7TVSxyXq9nqgiQ4qJeg/qC/cHHRgcNHorDjXKOGLtP7+pZQ6p9xb2MYcjWpJt8cfIBcMLXukW3i3KY5J80/am8rStNvQhtji6l4qrnnA2aZ1rmSkDT0lnGuawo1Wikj6XvL5WyWBIpcSp1GthCJ4wyb3d6TLk2LPuk0KvS/aM0RvxRpXmDpkJuzRU3BnzseCfFb/VQmj+KV7lXP3SeEeNm9JDy2yWk+O6ORzXKq9V4Q+DHth2CDcEIX8oT9tqIHY8Dtf6fe3K9U/8QNd+m/3xkrazxxx2tFbp57i+1irr1HqVeC30MlMbwE0y4eFVdPG78vivZ554c8Si3rzaaJJeF3PGBJU/LuZQtEOUhxLqffnjbDfUUXpome0T5pY2h+MaFSNeODwoDOWp/Wyrg4hWPWFc/5t+wfJAUXe7zWFHu8JXtjLH9bUyEZNVp7aLoq6FQ/DEI/IBDu275m1T9Xp9pcLxqovgJV9ISNZ5H1S7Z2yl1Xe1CiJgWy1I2TOKajuaEa6LNoG1E1isZ2wMCDldr7QwfBmGNOnVyHluyRg83ndFBOawVZVJrzB3ztyE8KRndR1wrhWJ7qokJDCHHSWs3tmK8rVwNvQ+8cIEVumMwUAXIfMtHExg9QRS26rKowfY9r96zaetYGo7Slu785Tma+E83H3W5Qms1hlLjUIXtXsbBNUcUluWyzMBEzABEzCBaQRKcsbc2yBbkTNmDtOuGR9tAiZgAiZgAkkJrHUCf1IhXZgJmIAJmIAJmIAJrJXAWXrG1DvDO8Oe11AYPiQ8VX7fe5uu9jjwO0eZB07nTT0ErIcYjHnEPNqpUviUImdbv06bgAkMI3AWZ0xV40lJnhargtZfa6X6WHmddUo0R5mn1KP0Y6yH2ALMI+bRTpXCpxQ52/p12gRMYACBcw1T/iQH7GGjfjy6zyP8fDvw1DBHmafWpeTjrIdY++YR82inSuFTipxt/TptAiYwgMC5nDF6xf47oH5jdpmjzDHn975XBKyH2BLMI+bRTpXCpxQ52/p12gRMYACBswxTqges/b0pGio+LRN97mhA/Xe7zFHmrvAzr0g25tWF907xcfVL0lN4zSXSlvUwB+ct85jLRprllsJnS3KqroxQPFNMu+yQgIBYbqYNTyDuoCJyZVJfP3xv857WB8+DP4sz1tRUXfHHVLyZP2V9jjKn1GfKsZKFi/i14kehHK0zrPuJPC0fQv7aYtWNRj2pbueSUXWdnfOWeMzF+VC5pfDZgJz8kZx8Y3xIxyVuW6Jt2RrX3JhIHh5KfKOFjhI6TPg/GRWuj9o78c61QnAsRnmQh6oxR5mHzrfANvhEd6iSkV4yPO7wtfsFqjHuFBvUw6ycN8hjnMIn7l0Kn7XLqfr9PFGVPnyfwKxty/7pNpGTFRNdN1+1PNHCf/XbUzRwNmdMlcZzZKiN3h0EYQL/aG+yKXR9fNIym+WfaZ0HHf6UbHjezUCP2M2pzJoFplrfqB5m47xRHqnM4Wg5pfBZu5yqHz3Z3OSxOKQjMFvbkq6Ki5dkJi3kN1rpRZJ1o8TrLHCcwhOUeJRhXtToesxR5uhKzHMATtddydfXQLadtHlqMbDUDethFs4b5jFQ49N2K4XPRuT8f6onbXLUEz9Nwz5aBGZpWzZO1kxaCkzijOnixSH4ScsdLcxliiboK81wWvOlrrxTjGOId6HZCMxR5u5EG1oRhyc91a2cWG2v5nacwKun2Dj7hHKP6jY+wzpSM3I+yuMExuuA1lGLE2Q5yqfjNGfPyk1OycPwJDfIneEEeTvLaWbOUWaz/LWsS86ztuFr4dCsh5k0aVytJ3HGVNRzweWOisnaTGLbOWPKw0l73ISv9W+VdyzMUeaxc25iu/jhiFXDvI0Kj+LVOO7Y6qhyB+r22DlXsT0F54E8RjFeBZz+SoySZSCf/rOdb0s2ckoHtCVMFeG7tX1hlLx9hbTy5yizdYp1JsV6yTZ8nRBatSqdyWRnrAb4sebKE388TdAM5I16MmeOMpsVymCdnsZ34vQKWU7lpeOYREn4ouW20tEw8anlViXm8ZOE8yEUpzLWcfQs/6iFVxAkexL5UF2PbTtVlmPlrm17hnJiQ9G132R+qrw6rrd9mVDm6uy+yWrEepK2RRxz4QG6VEx67W6EfhbfdbIzpho33w/Gn8PLlhQPO/Jau+wl5yhz7yRLZNQXyx86FxfN0MBTGZ0OrPIZSoBPs+t7NC8dT4/lF8XBoXuo9fdadq/QqM8T6pFKt0MZjNpP9Ybv6jgPEOIU3XFXzXVFGGNXV0fM9ztalvmqMmvJ2cip64Z2oHd4sqY4Wt663EPtyyllnsXu19q2qF5n4YFNrJjJkP+12qzXFU12xqSUr7Vy+HPgj6E5REn3N3lM1ouCjuv1XieUybnW1lsAnyQ9F+KCod1S3HSYLk7kBf8HWqqgMj5owRnjqdZquEIxdb9QPFa3i+uhrusaOcMvqa2rPBzkz4qZFrCaoPpsxl6mQJsgZ68dTKnPqcdKjqp9VnxoePLiRHkPti+nlKljzmL3dV1X17aciwf2tlYmqtpBuzv1WlniuMnOWKOS9NRwt1M1yHV+9WdRG81uV6WHeq9jyjzbXcJOsBlXxAyWdxTvesS0TmN6oTg0poN41cfhMIXjKIaA7uDYzh9ULgWo7E3rQfVPxrnmkdzWKXfloRR7GSPnUDtYUrW0H/dl8wwPNQPXMDdl5NOmhyHMQfJq/3AT3m5HutqXQWU2K7fVdXFJ2rZslUOz3imZjLS7ZjVWsX49YS24AMOQViiWHpyqV0yg+FAuDgAB73XXW6Z81h/WMNkewuAydSw9BQy5tRuAUNZmY8lF40ijGRrFIAsX92VIKB7KK+ihcWi1Slm32plKDy33QnXcrB5m4AzK5LbeoZ+1ZRVhL4I+WM4RdrCYLmXv9IYzJSJaVAGcprCt2eYMlXdM+zK0zMW4zHGimdqWOaq6WJkzMBljd4vJOfREKXvGcIJ2f+Q16B+UF4Yt6dWpXu6qPKC1nSYaAJyOZv6gMnVMtkHMaKy4Q6VxbM/twIGt5nzVAAbx0jH17nsR+usy6EHl7pW2oYyZOKO7Em09e3upTXuQnLVtDbWDNVw11PXkdmBk+zKI4RqgnFqHOdqWU+uyluPmYDLS7uZEcbsunP/Tr0NPlNIZ4w7qjYDQE8DTeVxkjLP/qryfFb/VQui6yMnv6pUZWibH5xreSzD+1BnmaId2T+RQXn0Ggm7QWzsMLbd93JbSc3Au1dZLsBdse6icY+zgbNeM2mlu9mhrCNzocRP4UXG44Rsq75j2ZWiZVaU2+jNH27JRFLtqz8FkjN3tKpJqpb5eKO5hXebvyuP/lLnYoVOq3rQfJXPGdDJAMP7fDl157X1I7/XKJCiz6zybyhMDXqQ7KAzlpf3C3D4a3rZD105fDC13UCVXupNkTM75gKhZ23oJ9oJuE8i5ZwcHbGb2TZLn4Jv3h8qr/Qa3L0PLnF34GU8gGZdsW2aUJF3RczAZY3fpJPm3JJ1/qK/z70GNteuN9aVWz+q9LiXkBs7DK0iCB38hQ2Kdd5d19YxtQJxVVtG2vkq1LF6pEu3A7cviZuYTisBm7S5Zz9hQM+DPXguN06BemaHler9xBKSDV1oYQg7Dyty9PR1Xivc+REBsZ7F1lcu1Qy8GDjRPvVU6VByGkw5Vy9sWJjCXHSwsxqjTYYvYJYsOZNrK5PZFZdnuG1owjwaMelVMktvd/lnmyVncGavFCN5rNSQmgO6VmUe/B0sVd+ZrOMxLILmtS2/0Xlp38+otdenJ7SB1BVOXl7p9sd3HGjKPmEdIpba7UO7c8bUXL17c1Un4WC9POy42RKVzccdECHdNfNuyrzv/as8Dvzq2edeETPQS7N4AfeBQb0pIwHrYh5na1vfPsN2ckuzFdrBdO3XNTWAOAmoTeDDvteJrZ3PG5hDMZZqACZiACZiACZjAFgg0nbHmMOWf2tCuPxO6Jz0h0C7QaRMwARMwARMwARMoiYB8qT8lLyN4nQFnjKFJJgN3hcWGLbtO7jwTMAETMAETMAETyIBAmJrVKcr/B/cwelIWOd2xAAAAAElFTkSuQmCC\n", "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\- u_{0} & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\- u_{1} & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\\\u_{0}^{2} & - 2 u_{0} & 0 & 1 & 0 & 0 & 0 & 0 & 0\\\\u_{1}^{2} & 0 & - 2 u_{1} & 0 & 1 & 0 & 0 & 0 & 0\\\\u_{0} u_{1} & - u_{1} & - u_{0} & 0 & 0 & 1 & 0 & 0 & 0\\\\- u_{0}^{2} u_{1} & 2 u_{0} u_{1} & u_{0}^{2} & - u_{1} & 0 & - 2 u_{0} & 1 & 0 & 0\\\\- u_{0} u_{1}^{2} & u_{1}^{2} & 2 u_{0} u_{1} & 0 & - u_{0} & - 2 u_{1} & 0 & 1 & 0\\\\u_{0}^{2} u_{1}^{2} & - 2 u_{0} u_{1}^{2} & - 2 u_{0}^{2} u_{1} & u_{1}^{2} & u_{0}^{2} & 4 u_{0} u_{1} & - 2 u_{1} & - 2 u_{0} & 1\\end{matrix}\\right]$" ], @@ -409,7 +408,7 @@ "⎣u₀ ⋅u₁ -2⋅u₀⋅u₁ -2⋅u₀ ⋅u₁ u₁ u₀ 4⋅u₀⋅u₁ -2⋅u₁ -2⋅u₀ 1⎦" ] }, - "execution_count": 28, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" } @@ -421,7 +420,7 @@ "S = sp.diag(*method.relaxation_rates)\n", "c_eq = sp.Matrix(method.moment_equilibrium_values)\n", "\n", - "N = get_shift_matrix(raw_moments, stencil)\n", + "N = shift_matrix(raw_moments, stencil)\n", "\n", "N" ] @@ -483,10 +482,10 @@ { "data": { "text/html": [ - "<table style=\"border:none\"><tr><th>Name</th><th>Runtime</th><th>Adds</th><th>Muls</th><th>Divs</th><th>Total</th></tr><tr><td>OriginalTerm</td><td>-</td> <td>899</td> <td>1204</td> <td>3</td> <td>2106</td> </tr><tr><td>sympy_cse</td><td>195.11 ms</td> <td>393</td> <td>246</td> <td>1</td> <td>640</td> </tr></table>" + "<table style=\"border:none\"><tr><th>Name</th><th>Runtime</th><th>Adds</th><th>Muls</th><th>Divs</th><th>Total</th></tr><tr><td>OriginalTerm</td><td>-</td> <td>899</td> <td>1204</td> <td>3</td> <td>2106</td> </tr><tr><td>sympy_cse</td><td>321.34 ms</td> <td>393</td> <td>246</td> <td>1</td> <td>640</td> </tr></table>" ], "text/plain": [ - "<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f9b78f46610>" + "<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f7cb0f678b0>" ] }, "execution_count": 8, @@ -550,7 +549,7 @@ { "data": { "text/plain": [ - "<matplotlib.image.AxesImage at 0x7f9b78828280>" + "<matplotlib.image.AxesImage at 0x7f7ccc65cac0>" ] }, "execution_count": 11, @@ -593,7 +592,7 @@ { "data": { "text/plain": [ - "<matplotlib.image.AxesImage at 0x7f9b73f22d00>" + "<matplotlib.image.AxesImage at 0x7f7cc7f6d190>" ] }, "execution_count": 12, @@ -628,7 +627,7 @@ { "data": { "text/plain": [ - "<matplotlib.legend.Legend at 0x7f9b73ed9370>" + "<matplotlib.legend.Legend at 0x7f7cc7f18b80>" ] }, "execution_count": 13, @@ -677,7 +676,7 @@ { "data": { "text/plain": [ - "[<matplotlib.lines.Line2D at 0x7f9b73e689a0>]" + "[<matplotlib.lines.Line2D at 0x7f7cc7f45730>]" ] }, "execution_count": 14, @@ -732,7 +731,7 @@ "metadata": {}, "outputs": [], "source": [ - "cm = get_central_moments(u = u)\n", + "cm = get_central_moments(stencil, u = u)\n", "mtrr = OrderedDict(zip(cm, ω))" ] }, @@ -801,7 +800,7 @@ " " ], "text/plain": [ - "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f9b73948f70>" + "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f7ccc63c850>" ] }, "execution_count": 16, @@ -834,6 +833,7 @@ "outputs": [ { "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABd8AAADtCAYAAABOHAzIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2d65EcN7KFhwz93tAjYg0gPSC1FojygJQHEj0QQ79G/xSUB5Qs2Et5IMmCvaQHogEbQS5jHbj3ZE9VT3VNPRJAApUATkX01BNA5ncOqnrQ1dX3rq+vH11dXb3Ba2n67ccff3y2tIPbSIAESIAESIAESIAESIAESIAESIAESIAESIAESIAESKBXAhg7/wu5P1jKH/vufTLZ8TOW5eDp9G66wmUSIAESIAESIAESIAESIAESIAESIAESIAESIAESIAESIIETgZcLHL7GtqeyfTr4/gqj8RxsX6DFTSRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiQwJYDx9F+m67KMbTI7Db7flyVOJEACJEACJEACJEACJEACJEACJEACJEACJEACJEACJEACdgSmd77b1ZpQEz4ZkGfQv8brMZY/JlTFogcQoH4HQO+oSfqrI7Edpkr/ORSloZDor4bE7CgV+rYjsZmqigD7hAoTDyKBMwH2mTMKLpCAGQH2KzOUZhW5GHyHMT5FRr/i9QGvL/FafEg9tnNySID6ORSloZDor4bErDAV+q9C0SoKmf6qSCyGeiZA355RcIEETgTYJ2gEEggjwD4TxotHk4CGAPuVhtJxx3gZfJc73J8JBhjme8zk7ndOlRCAZtSvEq1qDJP+qlG1dmKm/9rR0mMm9JdHVRjTHgH6do8Q9/dGgH2iN8WZbyoB9plUgixPAncJsF/dZeJpC5/57kkNxkICJEACJEACJEACJEACJEACJEACJEACJEACJEACJNAEAQ6+NyEjkyABEiABEiABEiABEiABEiABEiABEiABEiABEiABEvBEgIPvntRgLCRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAk0Q4OB7EzIyCRIgARIgARIgARIgARIgARIgARIgARIgARIgARIgAU8EOPjuSQ3GQgIkQAIkQAIkQAIkQAIkQAIkQAIkQAIkQAIkQAIk0AQBDr43ISOTIAESIAESIAESIAESIAESIAESIAESIAESIAESIAES8ESAg++e1GAsJEACJEACJEACJEACJEACJEACJEACJEACJEACJEACTRDg4HsTMjIJEiABEiABEiABEiABEiABEiABEiABEiABEiABEiABTwQ4+O5JDcZCAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiTQBAGPg+9fDGQ/b4Jwf0lQv/40L5kx/VWSNtuaE6D/5kS4bkmA/rKkybpKEaBvS5FmO7UQYJ+oRSnG6YUA+4wXJRhHSwTYr5ypee/6+voRYnqD18Mff/zx3VHxoe3XQ9tPMP8Ur7d4STy/Y98vmHNyTID6ORangdDorwZErDgF+q9i8SoInf6qQCSGeIcAfXsHCTd0ToB9onMDMP1gAuwzwchYgAR2CbBf7SIqegD0+A4NvsL8npvB96IE2BgJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJGBOYDr57fOyMcbqsjgRIgARIgARIgARIgARIgARIgARIgARIgARIgARIgATKEuDge1nebI0ESIAESIAESIAESIAESIAESIAESIAESIAESIAESKADAhx870BkpkgCJEACJEACJEACJEACJEACJEACJEACJEACJEACJFCWAAffy/JmayRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAh0Q4OB7ByIzRRIgARIgARIgARIgARIgARIgARIgARIgARIgARIggbIEOPheljdbIwESIAESIAESIAESIAESIAESIAESIAESIAESIAES6IAAB987EJkpkgAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJlCXAwfeyvNkaCZAACZAACZAACZAACZAACZAACZAACZAACZAACZBABwQ4+N6ByEyRBEiABEiABEiABEiABEiABEiABEiABEiABEiABEigLAEOvpflzdZIgARIgARIgARIgARIgARIgARIgARIgARIgARIgAQ6IPBJBzkyRRIgARIgARKomsCPP/74CAm8xusxlj9WnUylwVODSoVj2E0TYL9sWt5sydE32dCyYhKojgDPB9VJ1kTA9F0TMgYlwcH3IFw8mARIgARIgATKEMCbsk/R0q94fcDrS7we4MWpIAFqUBA2myIBJQH2SyUoHnZBgL65wMEVEuiaAM8HXct/WPL03WHoXTTMwXcXMjAIEiABEiABErgkgDdocof7M9mK5e8xk7vfORUkQA0KwmZTJKAkwH6pBMXDLgjQNxc4uEICXRPg+aBr+Q9Lnr47DL2LhvnMdxcyMAgSIAESIAESIAESIAESIAESIAESIAESIAESIAESIIGWCHDwvSU1mQsJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkIALAq4eO4OvYcjzbb/D6yFeb7D+y5QS1uXH5r4dvq4x3cVlBwSonwMRGg6B/mpY3Iyp0TcZ4bLq4gTo5+LIm22QXmpWWiZWCQH2wUqEYpirBOjhVTTcQQJVEGAfLiuTtzvff4ABfgaC3/F6OUWB7TIo/xRzeQYuJ58EqJ9PXVqJiv5qRcmyedA3ZXmztbwE6Oe8fHuqnV7qSW3m6pEA+6BHVRhTCAF6OIQWjyUBfwTYhwtqIne+/21ob5wXbP62KQyqyw/J/WvY8jXmH273npZk29vZtiuUGwfp32PfF1h/MT+G6/kJUL/8jHtugf7qWf343OmbeHYs6Y8A/exPk1ojopdqVY5xt0KAfbAVJfvNgx7uV3tm3gYB9uFiOv59bMn0sTMQUB4b8ydeMtdOz1BOBtXfDXMp9w1eP8nCZHqC5YttOF7uhn+Pudwtf4X5E7x+x0sG6jkFEgA36hfIjIfrCdBfelY88pZALb5JjPM2YS41QyCDJ/g+qRl3pCeS6K9iXkqMMx0UayABAwIZfFysDxqkzyoaJZDo62IeToyzUfWYVm8EMvSDYn24N63W8pXB9/8OO8f52rG722EIeSTM490DFw4Yyl5hLoPsMgh8ft47tj0Ytv0xKyp3vX81bsNxf+Alg+8P8Ho3budcRwDMqJ8OFY+KIEB/RUBjkatafJMSJ2Vuk4C1J4b6rjDn+6Q2LROUVYq/SnopJc4gIDyYBDISsPZxyT6YEQurrpxAiq9LejglzsolYvgkcCZg3Q9K9uFzEn0u/HtM29sz3yWuZ3jJpzAyEDxOT2UB286PncHyOCA/H2SXcvIIG07HEKB+x3DvpVX6qxelbfOkb2x5srZjCdDPx/JvqXV6qSU1mUuNBNgHa1SNMU8J0MNTGlwmgfoIsA8X0kzufPc2yaD6eZB9CE4eI3O66x2D7vKomf/BS+6OX5o+YOPnSzu4rQgB6lcEc7eN0F/dSp+UOH2ThI+FnRGgn50JUnE49FLF4jH0JgiwDzYhY9dJ0MNdy8/kGyDAPlxIxPuF2glp5uJOdgy2y13sX+I1Dsg/xLbpXfHzumXgfW1gfn4s1+0JUD97pqzxlgD9dcuCS3oCLfjmiyFdfris1936SC8atOBna21YXxyBFrzkpV/GKcBSRxHw4psW+uBRGrJdHwRa8LCX84EPRRlFKQJefNdCHy6lWVI7Hu98f4GMfsUAuzzP/T1eYobHeL3Etu8x/ydeMq0NwMvA+4WBTkfzTykC1K8U6T7bob/61D0162p9g+ve6yF5ec63TK+xTa5x8vsm599GOe3hnywEHGpQrZ+zCMRKUwhU6yWH/TJFB5YtRMChb6rtg4UkYzP+CVTrYYfnA/9qM8JkAg59V20fThajcAXuBt9hRhlUl+cOzaeLbThufC780tckxrvk53VwPTMB6pcZcOfV01+dGyAy/Zp9g9gvrn2RCFgsgYA3DWr2c4IMLJqBQM1e8tYvM8jDKjMQ8OabmvtgBnlYZYUEavawt/NBhfIz5AgC3nxXcx+OwH9oEY+PnQkB8hMOHu8GvIJxZPk3zHnnewjF446lfsex76Fl+qsHle1zpG/smbLG4wjQz8exb61leqk1RZlPbQTYB2tTjPHOCdDDcyJcJ4G6CLAPJ+h17/r6Wp6p/gYveZZ6dYPWiFkeTyOTPKLmIV4vsG3tkTRyHCdHBKifIzEaDIX+alDUAinRNwUgs4liBOjnYqibb4heal5iJuicAPugc4EY3i4BengXEQ8gAdcE2IfD5AGv71DiFeb3qh98D0udR5MACZAACZAACZAACZAACZAACZAACZAACZAACZAACZBAHgLTwffaHzuThxBrJQESIAESIAESIAESIAESIAESIAESIAESIAESIAESIIEEAhx8T4DHoiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiSwREAG3+U578/x+rB0ALeRAAmQAAmQAAmQAAmQAAmQAAmQAAmQAAmQAAmQAAmQAAmoCPyBo2S8/UoG3x/g9Qqvz/HiRAIkQAIkQAIkQAIkQAIkQAIkQAIkQAIkQAIkQAIkQAIkEEfgCYrJePtp8D2uCpYiARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARJYJMBnvi9i4UYSIAESIAESIAESIAESIAESIAESIAESIAESIAESIAESiCfAwfd4dixJAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAiRAAosEPlncyo3NEvjxxx8fIbnXeD3G8sdmE2ViJGBMgH3HGCirIwESIAESKEKA168imNlIRQTYJyoSi6GSQGYCPB9kBszquyPAPrUsOQffl7k0tRXm/xQJ/YrXB7y+xEt+ZJcTCZDADgH2nR1A3E0CJEACJOCSAK9fLmVhUAcSYJ84ED6bJgFnBHg+cCYIw6meAPvUvoQcfN9nVP0R6Ahyh/szSQTL32Mmd79zIgES2CHAvrMDiLtJgARIgARcEuD1y6UsDOpAAuwTB8Jn0yTgjADPB84EYTjVE2Cf2peQz3zfZ8QjSIAESIAESIAESIAESIAESIAESIAESIAESIAESIAESCCIAAffg3DxYBIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARLYJ+DqsTP4qoI8m/w7vB7i9Qbrv0xTwLr8UOi3w1capru4XBEB6lyRWAahUm8DiKzicAL08eESNB0A/dW0vEwOBOhx2oAELgmwT1zy4Fp9BOjh+jRjxL4JsE/51ic1Om93vv8Aw/2MpH7H6+U0OWyXQfmnmMvzyznVTYA6161faPTUO5QYj/dIgD72qEo7MdFf7WjJTJYJ0OPLXLi1XwLsE/1q30rm9HArSjIPLwTYp7wokSEOufP9b0O94zxDM/tVYlBdfgT0X8ORX2P+YVZKtr2dbbtCOblb/hu8nmP58Xw/130RoM6+9MgdDfXOTZj1lyBAH5eg3G8b9Fe/2veSOT3ei9LMU0uAfUJLisd5JUAPe1WGcdVKgH2qVuV24/77eIQMvv93WBnn477gOQwjA+F/4iVz7fQM5WRQ/d0wl3IymP6TLEymJ1i+2DYYVLbLFNLmTQn+3SWQqOlS/dR5iYrjbYkeKKZ3YpyOFWBoFgQS/VHMxxa5so7yBOiv8szZYj4CiX5eCozn0CUq3FYNAfaJaqRioAEEEn1d7LyeGGcAER5KAmkEEr3KPpWG32vpf4+ByeC72QSzfURlUXefD2WvMJfBdBlIPz/vHdseDNv+wPw8YbsM2r/F/Ol5IxdMCYBttKZLgQz1XWFOnZcAOdyW4oGSeqfE6RA7QzImkOKPkj42TpvVFSJAfxUCzWaKEEjx81KAPIcuUeG2mgiwT9SkFmPVEkjxdcnzekqcWhY8jgQsCKR4lX3KQgHfdXh75rvQeoaXfOojg77jdBpcx7Y7j50ZD+C8OgLUuTrJkgKm3kn4WNgJAfrYiRCNhkF/NSos0zoToMfPKLhAAicC7BM0Qu0E6OHaFWT83giwT3lTxCgej4Pvcpf7fJBdnvd+uusdA/Df4SV3xnOqmwB1rlu/0OipdygxHu+RAH3sUZV2YqK/2tGSmSwToMeXuXBrvwTYJ/rVvpXM6eFWlGQeXgiwT3lRwjiO+8b1WVT3bloJBtrlh1i/xGsckH+IbdO74qeHc3mfwBfDIZ/vH5r1COqcFa+7ylvQ20vfcSduRwG14OOO5KouVfqrOsmqCdjL9Yser8YyzQfKPtG8xEywEIEWzutezgeFJGMzzgmwTzkXKDY802e+xwYxK/cC679igP0l5u/xEvM9xusltn2P+T/x4hRIAOxeD0XkWesyvcY2Yfs75ufn65/2lPlDnctw9tJKtXo77DteNO0xjmp93KNYFeZMf1UomueQHV6/6HHPhukgNvaJDkRmiqUJVHted3g+KK0d2/NJgH3Kpy7JUbkbfMdJUO5ql+cczaelbfNjuL5CAFxd8aPOK0I1urlmvb31nUYtUkVaNfu4CsCdB0l/dW6ADOl7u37R4xlEZpVBBNgngnDxYBLYJVDzed3b+WAXNg/oggD7VLsye3zsTLu0mRkJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkEAXBO5dX1/LM9Xf4CXPUr94vpB3AohXfozgOV7yKBXJ42e83mO7zDk1QoA6NyKkMg3qrQTFw1wToI9dy1N9cPRX9RIygR0C9PgOIO7ujgD7RHeSN5cwPdycpEzoYALsUwcLoGgeGn2Hw15hfq/qwXdFrjyEBEiABEiABEiABEiABEiABEiABEiABEiABEiABEiABIoQmA6+87EzRZCzERIgARIgARIgARIgARIgARIgARIgARIgARIgARIggZ4IcPC9J7WZKwmQAAmQAAmQAAmQAAmQAAmQAAmQAAmQAAmQAAmQQBECn6AVec67PDf9Q5EW2QgJkAAJkAAJkEA2Avh626eo/IehAfltFJm+xfaPN4v8SwLrBOifdTbcUy8B+rpe7Rh5HAF6Po4bS5HAFgH2qy063EcCdwmwz1z9ASoy3n4lg+/yj/krvGQj/zEHBE4kQAIkQAIkUDGBl3ijc7rISw5Ylmv86YfVK86JoZcjQP+UY82WyhGgr8uxZks+CNDzPnRgFG0RYL9qS09mk59A733mCRDL/+K/8LEz+c3GFkiABEiABEigJIHvMOAuF/pxeomFB9j2aNzAOQlsEKB/NuBwV7UE6OtqpWPgkQTo+UhwLEYCGwTYrzbgcBcJLBBgnxmgyJ3vricMFvDr864VWg+O2q2z4Z50AvRXOkPW0CwBuev9f5vNjonlJkD/5CbM+o8gQF8fQZ1tHkmAnj+SPttulQD7VavKMq9cBNhnBrLuB98RZ+9fU8jVCUrUS+1KUO63DfqrX+2Z+QYBfDD1y2y3vOl5h+1vZ9u5SgJ3CNA/d5BwQwME6OsGRGQKQQTo+SBcPJgEVATYr1SYeBAJnAmwz5xRXNXw2Bl+TeFWr9qWqF1titUVL/1Vl16r0eKiPP4o6Oox3HFDIJQVjpdHzTzF6zEZLhMIZbpcSz1bQ/Klf47TNUSn46L003IIL/raj26WkYR4wLLdo+oKyZeeP0qlY9sN8cixkfpoPZQX+5UP3bxFEeojb/GHxhOSb+99pobBd35NIbQH+Dme2vnRosVI6K8GVMVF+HukwWeR67WUZ7cLs91peDMkz3t/jGX+oPoCsU79p/IQ/bNgmEKbOvVlKl36OpVgxeU77TP0fMWezR16p30iFauqT0kj4Cs3DvE9dirxxsp32u9U/YZ95urq3vX1tQx6vMHrIYC88+5/xCgnuaeYP/QeK+O7JEDtLnlwzZYA/WXLs0Rt0EzuyP4H5i9KtFdLG+Ah1+XXeC0OmmO/DL5/xHz+eJlzitgn/xS8wFw+pLoa1mXu/jp/TiLzAlhs+g/7N3XIHF7W6pHbpoewn/7JqsB65Xu+XC/Z/p69Pklft++BpQz3+syeb5bqrGUbPV+LUmXj3OsTZaPx1dre+WCvT0k2OIbvkXzJ6iIa+IL/V6z8b9pzn0Hu38GgrzC/V9XgOwJu9h9hF2eMjEFQu4xwWfUV/VWfCaCZ/Jj2n5jzcSgAMfD4FYsf8PoSL7nefYbti3esY7t8aP7V0n5sk38KXuE1/VBDBuFlMH6xPuzragKHRf8N29U61AwNuS56CNvpn4OEXfPlQeG4aDa0T9LXLmQrFsRanwn1TbGAMzREz2eAWnGVa32i4pSSQw89H6z1KQkE+/geKVmR9ioYPHbn/9phO/+v6Pj/UnjgPPheww+unnrncKLjV3sqPFdRuwpFqyhk+qsisS5DlfO5DBBzAgH4WAbFnwkMLGsexSPshOHpznYpN5lkUFUGl2V+nlDv0rHn/Z0tLPovQoeasa15iP45TtVFXx4XzvEtR/RJ+vp42UpGsNhnInxTMmbrtuh5a6J117fYJ+pOKS36iPPBWp+SQPgeKU2OVksv9rsI79XMZ63fsM8MqlYx+A7Tjl/t+VriHtZlzq/PO++eg1ZytyW1c65VjeHRXzWqdo75G+jHweAzjrAFsPsFr//gdedudmz7LKy2Lo/u3n9rHqJ/Du0P3fsylT59nUqwuvLd9xl6vjrP5g64+z6RCnitT0m92Mf32KmA2yzffb9b6zfsM7eGdz/4DrHOX+3BsnwNX6bT1+dvFvnXKwFq51WZNuKiv+rVEdrJM/H44Wm6hMLwG7xWn/2e3kR7NdB/F5rSQxc4jluhL03Z09emOH1Wxj5zoQs9f4GjzxX2CVPd2adMcbZbGfvdhbbsNxc4Lldk8P1vw6ZxfnmE8RrMKV+Fl+feyA+mvsH6xaAB1uVH5r7FXL6CLxO/pnDD4fC/1O5wCZoOgP5qWt55cvJNmD/mGy3WI3xk0exRdQhDYXlxHT0qmIrazeY/SwaFvEwPWYqWVlc2XxbyUlr2tqXpa1ueXmvL1mesEy7QB+l5a9HqrC9bnyjgYW/E2ae8KeI3nmz9zjLlQn2Y/eauaH8fN8ng+3+HlXE+7ss1/wHCy9fk5c5H+fGB86ABtsmg/FPMT8+9lQCwzK/2CAgfE7XzoUOrUdBfrSp7Ny/5QVF5LlyOKchHOQIoWOdfaOtJwfZaaSqn/ywZlfAyPWSpWFpdOX1Zwktp2duWpq9teXqtLWefsc45dx+k560Vq7O+nH0it4e9EWef8qaI33hy9jvLrEv0Yfabu4r9e9xU9LEzGEiXx8b8a2hcPiH6MAYy2fZ2to2rDghQOwciNBwC/dWwuMupyTeg5uf/5SMDtsb6COXkR3Jkeo/XF1h/cVrz/0cYyqPZOIURyOK/sBC2jy7oZXpoW4qSe7P4sqCXSrLaa4u+3iPUxv4sfcYaTaE+SM9bC1dnfVn6RCEPeyPOPuVNEb/xZOl3lukW7MPsNxvCFR18RxzvIPzbIR55Tu1Ps9jkDr75ttkhXD2IALU7CHwnzdJfnQg9pPk55h+XUsY1Qt7A/ImXzLXTs+HaEuwjlJNvXL3H/GdpDPMneP2Ol3xAvDphf0qcq/UG7pDn6oVwCqy+2cNX/ZeSsbEnSnmZHkoR3bbsqi8TvVXKS1eJcVrSpK8tafqta7XPpIScwccl+iA9nyJ6O2VX+0Sir0t4+KRCYpyWSrJPWdJsu67VfpeStnFfKNWH2W82RC86+A4DfZRYMJdBdhkwmD5yRu7ek23ynKCLCcfXelfiRR41r1C7mtXzHzv95V+jUhEOXngc016kj+T68tXYHur4Ay8ZfH+Al7yBWJywT65nUXEuVhi3McubvbhQWMrSE0NdV5iHvF+K8TI9VIF1U7xV0EviVw/nRVGUvq7A115DtPZxoT5Iz3s1lJO4UnxdyMMnUilxGqNmnzIGyurCCFj2hYJ9mP1mQ+b7G/ty7pJnusunL/ImfZyeygK2jXfGn7ZjfbwrUZ4TL3cmyqDI76ed/HMEAWp3BPV+2qS/+tBavpImH7bmmlQ+wrVk/NB3Psgu16ZHuYIzrFcYCktOYQRy+y8smu2jc3uZHtrmX3Jvbl/m9lJJVntt0dd7hNrYn7vPWFPK2QfpeWu16qwvd5/I6WFvxNmnvCniN57c/c4y89x9mP1mQ62jBt9lwOPtLK7zrwTLgDteIpxMcifX+W54bJdleSyA1MGpPAFqV555Ty3SX32oLYPdOc/hWh+N15k5dXkTJZ/ce58kxvkHB95j9hBfbv9Z5pjby/SQpVppdeX2ZW4vpWVvW5q+tuXptbbcfcY675x9kJ63VqvO+nL3iZwe9kacfcqbIn7jyd3vLDPP3YfZbzbUur+xL+cuMeh5wkD6I6zIrwSPA/IPse3jMMAugyMXx2P9I15ShlN5AhdaULvyAjTeIv3VuMBDenKu/0fGVFU+2mhf3jisDcxvFDPb9cVQk8SxNcljb8br5tZx3HdJQOs/rQ6Xtduu5fYyPWSrV0ptWl/GtpHbS7FxhZTT9kn6OoRqvcdq+4zWN7lJ5OyD9Hxu9eqoX9snYrPJ6eHYmELLac8H7FOhZPs9XtvvtN7LSTJ3H2a/2VCv6DPfJ3G8wPKvGLiVu9rf4yUmEKFeYtv3mP8TL5nWBj8+YN/eoMSpAv4xJ0DtzJGywgkB+msCo+FFOce/zpif1kfyQe7SJNeeizcnSwdZb8P1b2Qiz/mW6TW2SRzyuLXzb6Sc9tz8keOeT9a5qCOw6b8IHXStxh2V28v0UJwuOUpt+tKgwdxeMghxuYqIPklfL6Nsbetmn4nwTW4+OfsgPZ9bvTrq3+wTBink9LBBeOtVRJwP2KfWcXLPJYHNfhfhvcvabddy92H2mw29PtnYl20XDCgDHvK8ofm0tG1+jKwffVfiUkxdbKN2Xch8WJL012HoizYMnd/idYXXI1m2bhx1qq4xOG787ZGlr+CZx7WXJ+LRXgOvcKzELPPzY9n26uf+GwJgtuk/7FfrkJspYsnmZdRND+UWMKD+PV8GVLV4aE4vLTZouBGxq/skfW0I3nlVe30mxDclUkU8Wc7n9JKjitAAACAASURBVHwJ9epoY69PpGaRy8OpcWnKI3ZeRzSgeEwwgb1+F+K94MYDC+Tsw6ib/1fs6HF/Z//Ru+VNytL0KTYWvytxKRBuWyVA7VbRcIcBAfrLAOLBVbxE+x7u2v4Jccin9KcJbxxk+TfMvV9j5M4FYcgpjoAX/8VFv1wq1Mv00DLHI7d68WWol45kNm+bvp4TaXvdS5+xphzSB+l5a/p11+elT4R42Btx9ilviviPx0u/syQV2ofZb3boux58HwY/ZJDt9CnKLJfidyXO2ufqBgFqtwGHu5IJ0F/JCA+vABrKY1QeYL50fi8WH9r/GY19gfn42DO5M+bbYgFENDQwE3ZLj6KJqLG/IgO7w/1nSR45qb1MD1mSt6vLiy9DvGSXfXpN9HU6w9pq8NJnrLlp+yA9b02+/vq89Amth70RZ5/ypkgd8Xjpd5a0Qvow+42O/CGPndGFdj5q/MTlNNgOYWu5K/GcQMcL1K5j8QukTn8VgJy5CRnoluecf525nc3qcV2RT+prml4h2Oc1Bew0Vhf+s2QT4GV6yBK8bV0ufBngJdvs02qjr9P41VraRZ+xhqfsg/S8Nfg26nPRJ5Qe9kacfcqbIvXE46LfWeIK6MPsNwrw966vr+URLt/g9T+Au/YoB0VV+Q5BXPI1Dpnkx1kf4vXCa6wSJKdbAtTulgWX7AnQX/ZMS9cIDeXO96eYy127nHYIgJP8KHkNj8XZycTH7h79Rw/58N5WFD36couHZh99raHU7jE99hl6vl0/W2TWY59I5cY+lUqQ5Xvsd+w3274fPPEE819k8P0RDn+D10Ns8P6M2+3MuJcESIAESIAESIAESIAESIAESIAESIAESIAESIAESIAEDiKAMfbv0PQrzO+5fub7QXzYLAmQAAmQAAmQAAmQAAmQAAmQAAmQAAmQAAmQAAmQAAkkEajhme9JCbJwnwTwyZI8TumHIfvxBx2/xXaXj1bqU6V2s6b/2tWWmZEACZAACfgjwOuuP00Y0bEE2CeO5c/W2yTAftWmrszqOAI99SkOvh/nM7acl8BLdOTzDxJiWX4E4vR4pbzNsnYSOBGg/2gEEiABEiABEihHgNfdcqzZUh0E2Cfq0IlR1kWA/aouvRitfwLd9Ck+dsa/GRlhHIHvMOD+ZFJUfrT3AbbJbxxwIoHcBOi/3IRZPwmQAAmQAAncEuB195YFl0hACLBP0AckYE+A/cqeKWvsm0A3fcr9ne8YLOXjQ/rujLHZy13v/xtbmOVIIJEA/ZcI0HtxXpu8K1R3fPRX3fox+kMI8Lp7CHY26pgA+4RjcRhatQTYr6qVjoE7JdBNn6rhznf5GsKL4fUMhvmAlzw+hNOEAPiMzzWfbG13cS9f7P8Fr+nz3aVTv8O2t+1Sictsj2Vcre2W0vDCMfRfuxYYM+O1aSTBeQ4C9FcOqg3VqbkWNZTu1V6+vO62pHZcLnseiavVb6m9fNkn/GpXKrI9j5SKo5Z2NLzYr2pRM1+cGp/ka72umjWseupTNQy+d/M1hNiuBMN+j7K9PU5FHiEjee9OOE7YPMXr8e7BnR3QqXdSVVZ7Txqi/1Jxuy3Pa5NbaZoIjP5qQsY8SXR67VZfe3ndzeM7z7WyT2yrwz6xzafFvZ32iVQp1dcZaYj9KhV3feXZr4I1Y5+aIKth8L2bryFMdFEv4gQgg8pfYP7bUiG5KOD1F17y+J5mJuTzhySD+XdbSWG/fCNAnvf+GMvTO+G3inWxDzw2vdMFhJUkwWa132CfyntSNY6l/1YYN7CZ16YGRHScAv3lWJwjQ8N1ZfPajf2r168j405tG3mprr287qaSrq88+wT/F6rPtXkj3usTeVv3XTvYrF4jsU91nZEMcSz/x/MttXl00Hzz/Zd5g5VUCC7sU0qt7l1fX8tdwfIYl4cA905Z7rDDEKMMpD7F/OFhQThpGAxkQP1PzC/u6B62/4p98oieL/ESjT/D9uYGn5GTePerpdywTS6K8sgiGcS4GtZl7t7nEm/OCQwWvZOzTe91D0zU/QbHr3pPcsV++s+76IbxQW9emwx5sqpLAvTXJY9e1+CDxWv3sF19/aqZH3JdvfZiH6+7NYsbETv7xOn9JvtEhHdaLbLWJ1rNV5PXwER9jcTxq31K2sN+Xms04Bs6ZvDQnXG3hlIMSmXgwT6loAZWcrPwK8zvuf/B1Wk+CJiPD5kCubmj+9XlptMFQQbZ5fn4V2DW+iNpJH8Z9DoNsEvOMiFvuSjKPhl8F9/IJMe8OC3xjzC7452escAnof1m0XvCkP7ry0nDOUbuhrj4ILQvCsw2FwH6KxfZKutdvHbDI6HXryqTH4JevPbyuluzpEmxs0/cvJ/n/0JJNmqq8GKfaCrDwGQirpGL1xlplteaQPjtHM5+NdGSfWoCI2DxPo7923D8OA8oXu7Q4UQnpufjQ26xfwMuv9yu9rc05C8c5G6w6SSfWD/BS+bjS56hK/+gcrq66t47qSbY8J5UTf+lAq6kPHwgH/Tx2lSJXrWFSX/Vplj2eLu/dm9ce3ndzW4/lw2wT9z8L8j/hVza85Cguu8TqdQ3rjNSNa81qYDrLM9+laBb533q7yM6ufP9v8PKOB/3uZkP/3zKHcxfS1DDusy7fXwIcpe7LLvNf2ZO4fANXucPIsDns9kxXB0I0DumVrjjPamd/jNl7LYy6Dx+7ZTXJrcq1RsY/VWvdjkihx/4vu8W7J1rL6+7t3B6WWKfuFCafeICR58r7BOmut/pU1I7rzWmjKuojP3KTKZe+9S/R4LFHzsD88odyvLcG3lm+xusnwdMJSisv8bsW8xPdyhjzseHCJi7kwz2nH4U5O4uP1ugX5DekZELB+Fx4aXIunooltU7hTT3ohO950WJxDhCfctrUyLwzorTX50JnifdrNduy5BD/R7RNq+9EdAaLMI+cSsq+8Qti56XsvWJAud1b7qxT3lT5Lh42K9s2Hffp+SxM6WnH3Dy/hmN/o6XfFX/PGG7DMrLj6lOHw3Cr/acCV0syA+p/nWxxedKqN4xWQgH+ZCGk45Abu+U0FyXaf6j6L38jEu1EOpbXptKKdNGO/RXGzoemUXua7dlbqF+D22b195QYm0ezz5xqyv7xC2Lnpdy9onc53VvurFPeVPkuHjYr2zYd9+nit75jkF1+eHLfw3aySdIH2Y6yra3020ow8eHTIHcLssd5XN+t3sdLMXoLWGjnOQmj5F5jmXNDxgKBw6+A4JyyuadgporU81+GL2XHXH+BmJ8izK8NuWXpokW6K8mZPSQRLZrt2VyMX6X9lEu5L0fr72WotVbF/vErXbsE7csel7K0icKnde96cY+5U2R4+Jhv7Jh332fKjr4Ds3e4eT9dtBOBld/mun4ZGHb7BCuDgQ+x/yjNY3hn58/Ua+cZLTTs4mu0zLBeg8Xd/GBTNoY5PlR2mNPFXf+Z9M7iR4oonlijJby03uWNI+rK9i3x4XKliskQH9VKJrDkDev3bHxZrieBvs94r0fr72xgrdVjn3iVk/2iVsWPS+t9onEc32J8/pJt8Q4LbVnn7KkWXdd7Fc2+nXfp4oOvuNkehosxlwGV2Ww9PyMbmyTO5dlmzwL6GIaTsIhd0JflOeKnsCgkeZu891Kh7quMFfrjWPlw5m3mD/dbeD2gNUT4u0hXNISAHvpp1EeGMpeYZ5V85QYtRyUx9F7SlCeD4vxreSDcnLN4rXJs7gOYqO/HIjAEFYJDP6MuuYvVRrjd5QJfe/Ha+8SfG4zIcA+YYKRlTgjkOLrQuf1E7GUOI2R8zpjDLTF6lL82mG/6r5PFR18n3S4Z1iWT1BlkG+cToOt2DbeGX/ajvVHWJCBPJlkoIPTDYEPmNXCQ613pLjCQXhw0hEo4Z3cmusyzX8UvZefcckW1L7ltamkLM20RX81I+UhiZS4dlsmpvZ7RKO89kZAa7AI+8StqOwTtyx6XsrdJ3Ke173pxj7lTZHj4mG/smHffZ864gdXRTq5y/3tTMPzrwhjUOM7vEScK8zlLmj5gVb5mgKnWwLCo5bnnKv1vk0vaEk+RaM/9MhKeCe35vps8x5J7+XlW7p2tW95bSotTRPt0V9NyHhYEiWu3ZbJqf0e0SivvRHQGizCPnErKvvELYuel3L3iZzndW+6sU95U+S4eNivbNh336fu23AMrkUMfJ4wiPEIK/IrwuOA/ENs+3g+gAtLBITVP5Z2zLZ9MayL2Y+acustX5cevXNUjjW1q/VOSk65NU+JTVNW22/oPQ3Neo6p3bf1kO4zUvqrT92tstZeu7XXL6u41urJ6Xdee9eo97WdfeJWb/aJWxY9L2n7RCyjnOf12JhCy2mvkexToWTbPZ79altb9qltPue9Rz125gUi+BUD7C8xf4+XnMjlBPcS277H/J94cdomIIxerx0CjuO+8ZE9r7FNOP+O+flZ+2vljbfn1ltyfG4cc8vVbXrHKPHcmhuFeVlNRL+h9y4R1r5WpW9rh95R/PRXR2JnSHXz2h1x/coQ4kWVOf3Oa+8F6m5X2CdupWefuGXR89JmnzAAk/O8bhDeehUR10j2qXWcve1hv1pQnH1qAcrOpnvX19dy1/kbvORu84tPM3fKFt2N2OSZ8DI4/7Bow44bA4u/EN4zzOXTuKYmrd44Tr7+Jh8o0BcBDgAvd97Rah6QZtZD6b2seKupvDbfVgOWgZ4I0F80wpQA/ODu2j2NL3VZ43ccw/d9qaAbKs8+cXpEK/tEQ55OTcVbn9Cc11Nzti7P64w10frrY79K07DnPoXcvwO9V5jfO+qxM2nqsfRI4CUWer/jWz6BFw6cwgjQO2G8lo6m95aocBsJkAAJkEAuArx2X13x2pvLXXXWyz7BPlGnc/NFzT6RzpbXmXSGrdXAfpWmKPsU+Mngu9ztLgO4H/DiVBEBfHoij495MHySVFHk66FKLnjJye2HIbfxUUR3Cg15y/GlH6NzJ5baNgzMXHhHdMRLpbkXzhIzYqH3vAjCOEiABEigAwKert2WuLXvA3jttaTeRl3sE3w/2oaT7bLw0ie053W7zG1q4nXGhmNrtbBfxSvKPnX1B+idbpj+BAsyiPQKL9nIHzkFhMqmZ4hXnu/+dWVxL4aLzikfBsknY5pJfNv7nf8aTmvHuPBOoOZruZTeTu+VJs72SIAESIAEhICLa7elFAHvA3jttQTfTl3sE+1oyUxsCBzeJwLO6zYZ29XC64wdy9ZqYr+KU7T3PiW/HyEMfnH/zHecuOXDARlglaDl+fQ/4/Ue22XOCQQGRk97YoJc5Yd5f8Pc7e8U1GDOHr2Tqgu9l0qwjfJD3+G1qQ053WVBf7mTxFVAPV67ee11ZUF3wbBPuJOEAR1MoMc+kYqc15lUgu2XZ78K05h96jRWe37mu/vB9zB5eTQJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJHEMAH0CcB9/lme+cSIAESIAESIAESIAESIAESIAESIAESIAESIAESIAESIAEDAlw8N0QJqsiARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgASHAwXf6gARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgASMCXDw3Rhoy9UNPzDRcooXufWW70XyGVbIUw+VrPSstEeSqZbU+Ue89QV45C6B3vxXc741x75rxIUDLPK1qGMhtGY3tcirxZy2DJiab2r5rdg87msx3xZzyuUdC1YWdeTKz2O9LfJqMact76Tmm1p+K7YW95XgxcH3Fp2TISeY8XtU+yhD1Z6rfDDk7TnGKmLr1D8p2tB7KfRmZem/GZD9Vfpvn5H6iE79V6WHqJXa1ucDO2V2zj9yocr+sZZrpx6I1pC81pxUz/ZONUwRKLq/SKPkHYU+iXlUixkLdeqBaA075ZXqwGje2oY5+K4l1fFx6LxPkf4XmP+2hAHbH+H1F16fLu2vdRvy+UNix1x+oZhTJAHw2/RPZLXVFwOX1X6DffSekcL03zJI+m+Zi/XWPf9t6WAdS8n6kFd157A9rYRfi3qlaKVhVtJ3ntra8koKc085Six7Htji4C2XkHhiNSSv+v+n2tMwxEctHbvV12P7i/Ah73WX5GK+3uIxe/Y8sMXhmIhtWkVeUe+l93jZRFdnLVteieUdQuLe9fW13M38Bq+HaPBdSGEe2z4BeEIG1P/E/PE022H7r9j2Aa8v8RIffYbtHzFvakJO0j++ajG33EINPrnjn9zteq1/4KHuN/RempL03yU/+u+SR+61Nf+F6pA7zpz1I9cqrp9rWgmbXvQK1WqLWU5Pea471CuhzL3lvuaBUA7e8gqJJ0RD8jqdT6u4Jqx5YE3DteNb3x7a10P6i7Aj77sOys38bovHblnzQCiHY7NIax25qs+ba7zSIqi7dKhXQnhryKA+uZH3Feb3PtEU4DFdE3gpZpkTgHlkkP2ZbMdy64+kkfyFw3PJl1MQgUX/BNXQ0MER/YbeS9Of/pvwo/8mMMosLvovQocy0eZppZZz2KJWgqQjvUK1WmWWx0r+a43wSihzbxAWPRDBwVteIfGEaEheN/9TCoda/6da1DDEMC0dG9HXQ/qLoCLvmWEKMJ+1ePjqogciOByeSEIAIf1mkVdC29UXjfBKCO8gPvLYmb8NJcZ5UAU8uHkC38CwvzSf5UaCQ/7CoanH6mykbLmre/+kwKT3UuidytJ/CQjpvwR4N0W7919FHqJWN+/1hIP2vU73zFLPEBX1j7VUu/dAoIbkFX6eWfPeUdu71zAFfGB/kabIOwU4ykYwT2zRvHj3HgjUsHteqQ4M5K1p7u/jQTL4/t9hZZyP+zjvnACMJ8/q5qOIbnwgHL7p3BJB6dM/Qbi2Dqb3tuis7KP/VsCEb6b/wpld0X8X0Fx7iFqFa0VmF8xSV1z3j7Xk6IELMrsaklcYr4ujnaxQQzMhdvuLtETeZrylIhVz0xYNKqMHLiDuakheF7xSV3Z5BzTw7/FYk8fOQGi5S0aeZfMQrzdYv7hTGuuvsf1bzOVRJaopR52qhis4qCCbr4Hj9EMPnrEU4iEchMeFtz1zWYutEC9pPpt/CuawhrHk9ma8J9AKakf/2biU/ovjmM1/ceEslyrUH717qAqtRMECemm1ysasQI7LneG4rVrmx0W43HI2Dyw3F7e1kJ80GlbBSygXYKbhFSd43lLZNCzAPC+ZsNq1+pN3GNeto7XMt+o4Yl82D1gnU6APazTMyqtAjtaypNSn4R1cv9z5bjH9ADF+RkW/4yXPGTpP2C6D8k8xVw+8D4Vz1HmOq/KFUmzkh1T/qoBVCR7C4UEFLDQhluAlceT0T6kcNDxzH9OS94RVKe3oPxtn0n9xHHP6Ly6i5VIl+qN3D9WilSiYWy+tVjmZ5c5xuScct1XL/LgIl1vO6YHlFuO2lvCTRsNaeAnl3Mw0vOLUzlsqp4a5meclE1a7Vn/yDuO6dbSW+VYdR+zL6QHrfHL3YY2GuXnlztFak5T6NLyD60++8x2D6o/Q6r+GluXTlg+zKGTb29m2K5QbB+nfY98XWH8xHpNQp9yBL48GeY46Ho/1tTRPYLPKe4OP8JzruXF4+V0JPEK9IhyqH3xP4OXGPwk5hGpe3tDLLTbhPUmtsHZZzl8JOcT0oWVHlN1K/918uy/0vUUW/1lKX9DL3j3kXivRPUGvkGufVqsszBJyrPX8KtJqmcuxnqYsHrBMsKCfNBq65yXsE5jlOM9Y2sGiriwaJjCv9byn6S+ilzfeIR638JtlHVrmlm1a1JXFAxaBTeso1Ic1GmbjlZBjrf1Gw3tqA9Vy8uA7WnkHMd4Orck/pz/NWn4y34bj5W7495jL3fJXmD/B63e8ZKBeppg65UMAaUsmEdnthDwlvj/xConzGcoJ5xg2e7zXWH2OHR/XdsZuT8x/3mwMjxivyHOfQvSax2m2nsgvhpe5fw7IIUbzq8Q4rTR34z1JKJFJjP+itEOoq+evA3KI6kOJcdJ/twRSrp/m/rsNK3zJ2BMx/THGy67OYQvUV88VC8cGbXKgV6h/tVqtMkvMuZQnrxLjDPLBzsFa5jvVFN+96oGUSIx1KeUnjYY18BLpYpjlOs+kWClH2VUNE30bwzzmWuzlvKfpL6KfJ96hHj/5L9EXlh7WMrds06KuVQ+kVJ5BlxJ9WKPhJq/EvGNyDO43iTGm2GJeVsN7XmZ3PXnwHYBOg7OYy8C3DE6en4mNbXKnsGyTZ+ZMJ/mk9qtxA477Ay8ZfH+AlwgbXCfKyMD0W8yfjvV6nQ/5Rd2ZH8MGHDZ5l+aUkv881hgeKBPjlc2T2TyunOsp/GJ4IRdz/5TOIVLzq5Q4DT3gxnuSUwqToazUob5e4NiY/rqJf4jD/Tk4Jc5NAGE76T9H7y0sPTHUdYW5uj/COjHXA1ceCrN/2tFH64X2Q8+fyVql5DyUvcI8tyelDflfI+o6kOaKO6WTmd+pseINlroMdYnWuf10mIaWvMQ2Q30yVzPDscXPM94sPnCLOp8MZa8wVzNH/jHXYmnDw3kvub+k5DGUFRZq3jg21OMniw5tRfnC2OPJzI3jObQ6a12G+q4wV3sKAEL7cLKGKXnH5Igywf0mJUZjUyXzXorn/tLGyG3PUO48cD7UcRoIH8CfNmF5HJCXTxOmk1wMHk03YFlV56xML6sqNoG85+w+YMOn841O11U8EmIXDsKjlUnFy7l/VDk0IFhr3hNJSmiX+/ylyiGxD3mwL/0Xp0Ju/8VFtVwqt5e9e6gmrURBlV7LUu9u1WqVm5kqxwbOryKIlvmueIUPyO0By3Ry+0mjYU28hL2KWaRIGl6RVWctlltDFfMGznta/V3wzuqocpVrmZeLSNdSbg/ootAflbMPazQswUuVox6Z2yM1vIOD/yS4xHoBGVR/O9t9/sVdXCjk61H/g5cksjSJWeQThumkqhN1f5wW6mRZxQYsQnjP0ckHJNJODZOKR4JXxJvzD4xq4LIWo4oXCnv2jyqHBM3X2JXe3pr3hF8J7XKfv1Q5INeUPlTaa0vt0X9LVPa35fbffgT6I3J72buHatJKVFXpFXnt02qVm5kqR7Co/fwqemqZy7GeptwesMw1t580GtbES9irmGU+z1h6wKKu3BqqmCOR2s97mv4ierngHelxC79Z1qFlbtmmRV25PWAR47SOnH1Yo2EJXqocG+g3Gt5T7VXL91VH6Q4Ssc8TgD/Civzi7jgg/3BHBElwfjFJrfMcT4MLqWyWeM8xiXb/mG9cWP9i2CZ1HjWl8tiLW74yNnp579ga9qfysvRPLK/UHGLbtSqn7TeteU/4ldBOe/6K1TM1B00fio1NU47+GyhFvF/R8NX6T6uDps3YY3J72fs5TKuV8G1Bry2faLUKYbbV3tq+3J5ca9dyu9YrWuaWsVnUpfWAloNFTGt15PaTRsOaeAnHVGZrWsh2Da+t8kft02oYG18q89beV3rnHauzZTnt+bX1PqflYMl+qa6cfVijYe4+Izmn5rjEreQ2rVc0vIPjtrzz/QVa/xX/xMrzi97jJcJI0C+x7XvM/4mXTB9vZnf+ysD7hZhY19Z5p7IONmjZhPCeYxPNXs83juvQddwnz7eS6TW2iYby/P7zs/9Pe/L/0fKIjURyfB5b2GE5La9s/jFgos3BoCm7KiL6TWveE5gltNs8fxkoqs0hpQ8ZhHlZBf134qHV7hJe2Nqm/yJ0CGs97Ggtj1gvez+HbWolKCvVK8wFN0drtdplFtP4pExuT06asl2M8IqWuW2g6bVteiCCQ3pE6zXk9pNGw5p4CUkts3Xq63s0vNZLH7dnU0ODsLTMY6/FBiHerSKir2v198L7btIHb8nI/ODM7jS/6YEIDncaMN6Qsw9r+s0mL6NctTkaNWdTTYRXNLyDg7t3fX39CKXe4CV3ps8Hv4Mr1BRAO//BcV9hLp/OnCYs/x8WkmNAPfKceRnwf3hTM/+CRTRvlP0LBJ9hftaqFaLISeUVHCdfr5EPFLr0FPJuxj9azb14vHfvTXWI0Q5lXJy/EEd0H5oyKL1M/90Sr9l/t1mkL4V6uRYPIU4X54p0hZZr0Pg3VCsvzBBHledXUSqU+bK6x2314gFLAqF+CtGwRV5T9shv9/+qEF7Tur0se9EQcVR53gvV3wvv0X+IZ9fj47Fe5qHMvcQ9xuHNA2NcqXPkpe7DIRp65IWYquo3Ibw1PkB98vj1V5jfk8fO/G0oNM41daQe8xMqkE8TThMCkeXfMC8y+D8029MshfdLgGrpju8Y3eUTPuHQ60T/HKd8795LJe/l/JXSh1IZpJSn/1Lo3Vw3Wrt+hnq5Fg95OVekOS6tdKhWXpiFejKNkm3pUOa2rafX5sUD6Znc1hDqpxANW+R1S063FMJLV2PZo7xoGOrTspTWWwvV3wvv9Yz87wll7i2jVj0Q0odDNGyVV0lfhvDWxPX38aBD7nyXxjHQLsaQSR5RI3cUv8C2ta9RyXGbE8rK3cnP8ZKBfLmb/2e83mO7zLufUnij7O8A+BzzJj4cQR5qrwzHyidV8uPB3U7IP7q/ouzh/hl0rOr8MMRM7wX016UO6sF/EhfiiO5DS3nl3kb/3RAeOESfO1D+8POftVe0Xq7NQ41qpXq/E6uVF2aIo6rzq/TJWObW/Tm1Pi8eSM1jWl7rpxgNG+WV9Twz1cbDshcNtT71wExiQLzik+D/azzwHmKPfi94lAaxzI+Kd61dDx5Yiy1lO/Lafe8So6EXXkPsVfWbGN57HkCd5zvfP9k7ONd+BCGfKJhNqE8Ghk3rNAvOQUWJvJ8hBXm+exMD0IFeeYW85aTR9VS7fwI196I1vQclDLRzcf5K7ENHeJL+a8h/lgYK8HJtHnJxrjDWSvveOFYrF8wCPGmJN7WuWOap7VqXd+EBy6QC/BSjYYu8cp9nLOW1qMuFhgE+tcjZoo6YzwnOtQAAIABJREFU/iLtHs4brLUet+BkWUcsc8sYLOo63AMWSczrUPbhGA1d8Kq038Twnku7un5/dQ93kMBAAB1HvpEgd77LD+d2Mw35NnPH/1HC9eqfFN70Xgq9y7L03yUPzRr9p6GkO6ZX/9XoIWoV/u3GXpnpev/6UTX2j7VsevVArIbkFX6eWfPeUdt71TCFd2x/kTbJO458CvO4FvOV6tUDsRr2yivVgbG8Q9o97LEzIUHyWBIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARLwTgCD+ufHzvDOd+9qMT4SIAESIAESIAESIAESIAESIAESIAESIAESIAESIIHqCBz2zPfqSDFgEiABEiABEuiQAD6x/xRp/zCkLj+YJdO3w9cab9b4lwRWCNA/K2C4mQQGAuwjdVuB+tWtH6MvT4B9pjxztlg3AfaZuvUbo+fg+0iCcxIgARIgARIggSUCL/Gm7/zD01iWH6N5g9fDpYO5jQRmBOifGRCuksCMAPvIDEhlq9SvMsEY7uEE2GcOl4ABVEaAfaYywZbC5WNnlqhwGwmQAAmQAAmQwEjgOwy4PxlXMH+J1wNsezTZxkUSWCNA/6yR4XYSuCHAPlK3E6hf3fox+vIE2GfKM2eLdRNgn6lbv1P0Lu98xz/0/Ip7gLnIKwAWD71DgP65g4QbSIAELgnIXe//e7mJaySgJkD/qFHxwE4JsI/ULTz1q1s/Rl+eAPtMeeZssW4C7DN163eK3uud7/K1ihfD6xki/YCXfMWd0zIBE17gPT7Ld7mVxrb2lu+GfMn+IcsNugu7yGsBiuEm8g2DuccL+3/B6+OkVnkD+A7b3k62cXEgsMezNVB7+WK/a//sxd+bXnv59sZrj8fefg0v731knqMmp3mZmtf38rXSb6+dmhkuxV5TvjXFusS69LY9Xthv8r5gr53SeXtvryZeNcVqoftevtjPPmMBOrCOPV0Cq7vyOvjOr1WEKZnMC8b6Hk329ggBeWyC5N37lOSfTr2T6hl6L5XgSnn6cQXM9ma1H8FXrhNP8Xq8XWWfezv1X7X+oV5h/bRTXmGQ7h6t7h9SFIxdn2M79YBaw1j9yPVux/GypVNtUvGzz6QSzFNerUue5nW1dtrn1NrwOqPzkdFRal007XkdfOfXKjTq3R6TxAsdWAZSvsD8t9sqb5ekg+P1F17yOKBmJuTzhySD+XfNJBWXSLR/wG7TO3HhtFEKbFb7DfbRexlkph/XoVr4EXXIt6Pkee+PsTy9E3694Y727PlvS4OaMSEv1fkMx7nyD/UKe++zx6tmD6fGDjYm13tvfWTOZc8DWxzmddW0jryynuPINexcVNI7e9qUjMVbW1v9/eg+441VyXgsdCkZ77ytvT63ld+8rprW2WeOU2vLU1pdtNHfu76+lgHVb/D6H1Tu8p9pxCX/7D/F/KE2sZ6PC+GFY0X/PzG/uItx2P4r9skjf77ES+7E+QzbXXoEsUVPyEkeafRVi7nFQAEHVX/DcYveiWmzlTIDE3W/wfH0npH49ONdkJZ+RF0ycCqPg5MP666GdZm/u9tyf1vAYfF8OGxXnxNqJodcV89n2OfKP9Tr1IdX9Zr7cI3X/Lie1gcm6r6N4zd5Y7+rPjLXcs0DoRzm9da0vqVhrH7kGnYuKumXNW1KxuCtrYGJyXnPus94Y1UyHktdSsY9b2utz4XmN6+3pnXkuvpeAfui3iesca2Ji3WsoZ7a0mUvtkG3J5j/Ij+4KiK+wks+1Xc3sIogXX/9EsxcTRG8ZKBV9L+YUI94QZ63f4Xl1h9JI/kLh9OgkuTc6xTon0Xv9MpO8ga/0H5D79kZhn6csbTyI+oZ3yfI4Ltck2WS8+WL0xL/CIFF/0VoUDPNxfOZU/9Qr5v3fsJB895nkVfNZk2NPaJvL/YPicNpH5kjWvRABId5vTWtL2qYqB+5hp2LSvplUZuSAXhrK6K/F+sz3liVjMdKl5Ixr7S12Oci8lupvorN7DMFZIrw1KIuylCf4Dgp/8t9ZYFDDgMU+YdfOiG/4q5QIJLXNyj3i6L6Zg8Z8hcOcudit1OEf7r3TqpZ6L1Ughfl6ccLHOErG36UuzDkjYPMx5f8VoS7D+zDszYr0b3/KvMP9bp57yccNO99uueVeqbY6B9SdQ3n2O49sKFhin7kGnYuSu2KIeW71yYE1tKx7DNLVI7ftqHL0cF13+c2tOF15kB3bugSFJXc+e5yQoIy8C532X0tAQ7rMudX3BcUG/gE8UIZeV43ed7wFA7y+KUuP4gI9Q+9c2Mao79de8+CIf1oQfFcxx0/gu9n571cuEOA/rtA4t4/1Gtbr4u9WCGvOZGk9Tv9Q2oDY9fnWHrgQvM7GsbqR67bXC/2Fl6hNqbA2WdMcZpVdkcXs5ojKmKfu4B2RxteZy74HLVyR5fQQIoMvsMscleN/KilPLP9DdYvBjix/hrbv8X8dBcd5rtfcccxQXWi/mqm0Nw0vFaSlw82Tj8itLLfxeZQHpFBCwfhceHNyLoOLRbKK9I/Wb0TmsOhwNMbb8Z7cxQFdczmx4I5zPEdtd6MHwtql81/1iYowKQG/1CvW2Np9MrGq4AfbzP1saTh7SPSyyiyeeCymfS1Ap6y1JBcbyW35Hpba/xSNm0KeDQ+6zwlLbWlLnYaWepiEVU2bS2Cm9ZRoA9bapONawEOU+welpN1KfXYmR8gzs8g9jte8hiZ84TtMigvP6Z6Gngfdmi+VhFa57nNChZCc9PwWkpbfkj1r6UdzraF8ogJXzjIhz4tTKG8YvyT2zuhOdSsW0vem+tQSsecfiyVw5zdUest+bGUdjn9Z+2D3Exq8A/1unWVRq+cvHL78TZTH0sa3j4ivYwipwcuW0pfy+0pSw3J9VZvS663tcYv5dQmt0fjs85T0lJb6mKnkaUuFlHl1NYivmkdufuwpTY5uebmMGXuYTlZl+x3vmNQXX6c7V8DLfnk5cOMnGx7O92GMptfv4ypU+pHOblbXh4t8hzLj2WbtykmN5TZ5LWRo/CY67FxePldMTwkSpQL1Vo4VD/4HsMLZWL8k807MTlEai7FPExNeG8OsrCOWfyYkMP4IfN7cPkC9dT0w6RN+DFBu9Brh1g/i//mfSp1PYFJiJ9r8E/reoV4WKNXFl6F/JjabazLa3hbt2lRXxYPWAQ2raOQpyw1bJ2r9bloKnfu5SzaFPJobjah9bvvMwm6hHg8lFvu4y11sYg1S5+zCGxaR4JXjnovnYVrAoeu+4xq8B1wBdKfeMlcOz1Dubc4+N0wl3Iy8P2TLEymJ1ieb5vsXlwMrnMwiLQlU0geNyXK/Q3OLSG0z1H2Y0L5xaKJfpnXGcwjUmt5hpNnX8y5rK0H81qraGf7pncSPRCcQ4zmiTHu4Ana7dZ7iYyK6DiQXvXjATnIt7neo135ttcV5k/w+h0v+aB5dcJ+Of/EXmdX643Y4caPiUxc+C+C/7lIYv7neiYLMUxC/ezGP5O854ur54v5gSHrTvR6hJhD3utq9FrllZhzCT+eJEyMM8QGe8dqeO/VccT+VQ+kBJNBlxKestSwZa45zkUpdgstu6pNom9LePSUa2Kcoby2ji/SZxLzjdEl1OMt67Klv3bfap/TVrB0XKIvlqqM8cqR76VXuSayieHQfZ/RDr7LAG3UneIQ9TS4i7n8MyCDC+dnamOb3Gks2+T5ORcT9q1+OhRTJ8rIBwFvMX960ZCzlZjcJAWUW+VVOsUhhyi/zGON4YEyMVqvnpjmMXlej+El+aCcqX+GOKI8EJMDygRrnhKjsQfcei+FUSkd97QonQPikb701RgX2v8DLxl8f4CX/AOyOGFf9HV2scL4jW78mMJkKHuFufq9B44NPo/EY94vOeQQdR5dqn2o7wpzNRPUE+pnN/5ZYpBzmwe9EEOoh5P0Ssl5KHuFeU4/niRPidPYM0m8jWM5vDprXYb6rjDP6Sn3GnrgihiKnotKmjmF71D2CvOcHj3hSInTmGeRPpOS71D2CnO1Ljg21ONd6mLspeDqBm35XnqBXAqboewV5uwzC2zXNt1f25Fh+zPUKZ+QnAbjh/pPA+HYJiev84T18dOhF1g+PSsec3le/HxS1zkvWMG6OrcAXvO0P2CDfPhRw6TmEZmMcBAerUxqXpH+KeEddQ6Vi9aa9+ZylNAxtx9VOaAvjR8ozwfZ5br3aA7G6XprflRpl6hFbv8lhnenuIpJpJ9r8E+Tet1RWbdBo1duXjn9qKNQ7igN73LR6FvK7QF9JLojc3rKUsMmueokunOUJdc7lUdsyK1NTo9GpJu1iKW2LnTJSqtc5Za6WESdW1uLGKd15OzDltrk5qriMAVX8XKyLiUH32VQ4u0M9vnXd2UAEC9JSCa50+p8Nzy2y7J8dV/qmE4hdU7L1bAckpuW1zxvGSCaM50f42U9hEdMzPKp/HzALKYeL2VCeMX4p4R3QnLwwj0mjta8N2dQQsfcftTmMF7D5gzkjY/oXMPUmh+12qVok9t/KbEtldUyifFzDf5pVa8lrfe2afTKzSunH/fyL71fw7t0TJr2cntAE0PIMTk9Zalhq1xDtBqPteQ61pkyz61NTo+m5J2jrKW2XnTJwal0nZa6WMSeW1uLGKd15OzDltrk5qrlMGVX63KyLvcLZi7CnycMpD/Civz67jgg/xDbPg4D7PIP38XxWP+Il5SZThfHrNU5LVDRsiq3QF7z9IX9P+YbF9a/GLaJ4Y6aVDwSgpOvI41eTKjGTVEVrwT/aL2TAkSVQ0oDmctq+01r3ptjLaFjbj+m5iDnzrWBzDmvXOu9+jFVO40eWv9pNdC0mXJMKpMtP9dwPutNry2vaPTS8tpqZ2tfTj9utWu5T9u3Nbwt47KqS+sBLQeruNbqyekpSw1747qml2y35LrVjnafVhttffPjcnp03laudW1/t9TWuy65WIfUe4QuIfGtHavVVpvfWjtW23P2YfYZK5V09Wg9lazLJ4jnb0NM41wXYvhRL1DkVwz2yV227/ESw0oCL7Hte8z/iZdMawMUH7BvPvirrfNUcWV/tLmF8JojEOav5xvHdegy7pNnOcn0GttEN3l+8fnZ/ac9+f9oecRGIjk+jy3ssJyWV6x/Nr1jxEObg1FzNtVE9JvWvDcHWULH3H7U5vBxnvywLv3s4k3aynHmm+nHK612Kew3/RehQUosmrJaJjF+ruF81qpeGu3nx2j02uQ1rzBiPacfI8LRF4no2xre+gDKHbnpgQgOuSPP6SlLDVvlGqOvJdeY9udlNrWZHxyxntOjEeHoi0T0d0ttveiiB1boyIN1schyU9uI/Cxi2qojZx9mn9kib7QvwlOxuvx9DPne9fX1I6y8wUvuPD9kcGAMRuaIYYznMyyf//HD8n+w+yfM5Rnw0RPKy3PmZcD/YXQljgoijyReKP8X0nmGuXza2NSEnFRa4zj5uox8oNCEJ0JERM7R/kFZd95BTCrNQxjlPBbxduu9La4xOnrxI+KQa9VXmJ/PqVj+P2xzcY3d4U4/AhD0Cj6PoIy78+GW1tp9yEvtZxxbjX9a1WvUFfntejhELy+8EIfajyMLL/MQ3l5insbhxQPTmCyWQzyVQ8NWuY7aID/Tc9FYb4m5F20QB897E8G96DKGpPH4eKyXOWJ2+X7Nm7ZWeiEvdR/OoY03rohn97pgxd6qnhRdUFZ+z/QV5vfuWwVkWM95wH1W52F3Ds7i8LaayuslEmrpju8YfeSTS+HQ45TiH3on3TE9ey+d3mUNXvz4E8KST8ZPEy60svwb5od/uD2EtDWjH7fobO/z4r/tKMP3hvi5Jv+0qleIwiF6eeEV4scQFiWODeFdIp7QNrx4IDTuveNDPJVDw1a57nGf7s/BdVp/7LIXbUI8GptrrnI5tPWiSy5mJerNoYtF3K1qG9KHc2jTKlcLz2nrMNHF3Z3vkj0GKdSfDmlpoU75hO85XjIQInf7yh3077E96U561HH4lMoL5X9HEs8xr2FwaJc38lBrPRwrn0TJj/92OSH36P6Gsi68M+hYVf8eYu7ae/MOl6qjIz/KmxyZ5BFr8o2aF4ht7YMuOe7wiX48vfdQXzuWBPPiv6XYUrYhr10/1+ifFvUadNi9Fsbo5YUX4tj1Y4rfc5SN4Z0jjtQ6vXggNY95eY2ncmrYIteBV5Zz0Vy/nOtetEEcPO9NhPagi9bjk7BdLA5xu/3/04O2OYRCXrt9OKc2HrgO+e1eF3LwT6kzVReUP9/5/klKIBnL/oS6ZZD89LV9BJx85yDqkIFl+cSixSmV1zNAkee7NzEAHaj1K+QtJ4GepxT/uPBOoOZetKb3ZkoY6OjFjzVea7r3Yyv+m3Wr5FVw0fi5Rv+4OF8kCzSpIMDDMXq54KX044SKi8UY3i4CnwXhwgOzmJJXlZ7KqWFzXDOfi5I1D6jAhTZKjwakVeTQpvtMgMeLwA5oJKcuAWGsHuqiz61GF7lD2YdzanM4V/aZqyuXd76LpyHO7qdDkd5vslgqL5SXO/6eYl79NwG0AiNX+aHfWh4HoU0r6rgU//TonSjIk0L03gSG8SL9GA6UfgxntlaiR//V7B/qtebk5e098lomod9ac/9YyrJHD5TQkFyX3OZjW4/apJJnn0klmKd8CV0sIu+xz5XQpkeuqX600AV1nO98dzv4ngqK5UmABEiABEiABEiABEiABEiABEiABEiABEiABEiABEigJIHp4Pv9kg2zLRIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARLogYDXZ773wJ45VkQAn1h9inB/GEKWR/TI9C22u/4RxZsw+bd2AvRf7QoyfhIgARKogwCvN3XoxChJoGYCPM/UrB5jr5EA+1yNqjHmIwnk6DMcfD9SUbZdE4GX6IDnH2bFsvwgxhu8HtaUBGOtlgD9V610DJwESIAEqiLA601VcjFYEqiSAM8zVcrGoCsmwD5XsXgM/RAC5n2Gj505REc2WiGB7zDg/mQSt/wg8ANsezTZxkUSyEWA/stFlvWSAAmQAAlMCfB6M6XBZRIggRwEeJ7JQZV1ksA6Afa5dTbcQwJLBMz7TDd3vmOQlI8NWbIUt2kJyF3v/6s9mMeRgDGBZv3Hc7OxU1hdEAH6LwgXD+6DQLPXmz7kY5YkUAUBnmeqkIlBNkSAfa4hMZlKEQLmfaabwXfIY/61gSKSsxEXBDBA88ssEOmM77D97Ww7V0nAnEDj/uO52dwxrDCAAP0XAIuHtk+g8etN+wIyQxKogADPMxWIxBCbIsA+15ScTKYAgRx9Rh478w4vGUj8UCCHI5sw/9rAkcnkbhtmG39UNHdTLuoPyRfHyqNmnuL12EXwDoMI4ekw/OIhhfBq0H8m5+YQhsUFdtYgWV0IQv9d4NCttOShlnLRqBeSL46Nfr8T0o4m7paPqY1VbfGmesciX4s6UvMoWT4kXxwbfZ4plVNIPqVi8tyOBS+LOjwzsowtlBWOZ5+zFMCgrlANl5q0qGOp3la3hfBK7DN/gKGMt1/J4LsMssqPR36OV8uT+dcGWoUFc32P3OSk3NMkz2+XvDenoZPK894fY/nj5sGd7uzUP6lq9+y/5HMzPRdsP5XfgmutswD9F6dbEx7q9Nyh0g5s5P+DqPc7nXKN60k3pVSapDRgVbZTbZP0IbN196WcZ9Zrtd3TqX6pENlnUgmGlVfzZp8LA1vwaLWGSzHxPLVEZXebirlBn5HfjZTx9qt719fXMsj6Bq+HqFjugu9iQq7yD8VTzB92kbAySfCQO7r/gfmLpSLYLn55jVdzg8/ITQbfP2I+f8TMCQW2yz+iLzA/fXI1rF9h3k2/OYHY+AMWm/7ZKNr8LrDZ7DvYT//BBeAQdG7G8fTcQu8BlyS/LVTZxSZr/+3pUDNU5LZ5zvKeG+LfPHf0rB1yj36/s8fVuy9yxbfnJ+x335/2tN3LMRfbEvXG6kNmdf9ftadfCe95bWOvv+fqM1555I7LgjfqiL62585vrB8x8r3ZyljUyGg+32M2P76n9dR+Y9FnUMd3YP4K83tdDr7vidCTIae5gov8KO2fmF88TmXY/iv2yaOJvsRLBnU+w/bm7vxGTvJB1Ffz3LAuFyv5xGr6oYQMwstgfHMckFfwBA6L/gmuqKECAxN138HxXfsP+W8OGM+tQc9dErHy22Wt/axZ+S9Uh5oJI9fFc5b3nAaN+H7H+P3OGlfvfsgVX+i5wHN/WtM2NMdcrEvUG6oPmZ1uqFi8RoCN+/+r1vQr4TWvbQxMkv+vWcuPzC/JWPJGXdX2uVAOlxTrWkOui+fMtSwGNnfez64d38P2UL+sMcd2kz6Des6D7z394OrJawPEqK/RdmBW4XL6SsQ0VzCTweVnsg3LcmeODJC1Okn+wuF0d/skSTkRyuCyzM8TeMyPO+/rcGHRPx1yOKcc0Xe69R9YyQUu9NxMz53ddjo/h56r1/w2qbWPRUv/RfT7miHX6qHFcwe1O1kx5f3OIteaDZ4Se4SfPPenRW0jckxBenTZUH3I7Ob/SuEw/38p5TxTygeL+pVq3GM7Ef3dpM94ZFEiJmPe1fa5CA4l5MnVBvtMItkIv6wxN+8zXQ2+QwgZ3JE7lb8WTYd1mfOxITcm/wYs5m+ObvZ08hf5/4LXf/C6uKMd6591giAlze79kwJPysJnXfoPeceem+m5BNOt+S2hyiqL0n/xslXsoe7PHWvaYXvK+53uucb3pvX3ACl1GpbtXtu1PrPBmMzqfl/bvX4b3lbtYp9RYTI7aIt34rXdLMadirrvc1sarrDrntkKF/XmNeY5+sx9dVSVHwh4Mrgjn2rI83YeyQvL8ggReZRK9xN4yPO1+CHEjROEwzfdmyIAAP0TAGv/0K78B+9EnZvpuX0jKY/oym9zJvTfnEjUelUe4rnjQmMz7cj1gmvKipkmKUFMy1LbKY3T/0q7/yOQWTizixIHr1A/UwFU5zQyN2Ou4m3WmlFF1P8CpEpDMrtglrqiYp7aSJE732EMeVyHPOtGftz0DdYvftAS6/IDnt9iLl+ZV00RdZp/bUAVaOJBEXnGtijfBvgjtnDJcgWYCAfhceHTkjlatVWA1RhqNv8UzGHM5eh51f6L0Cv23EzP2Ti1ar/NEbTgv3lOKesRPGKaq81D2c4dMfC2yhTQz1K7bFwLcNiSofQ+S02sYs+mrVWAYz0FvKLVpwpmBXiJNFpmo4we5tn0K8TcA8MxBq3+ZD4SS5treae1Yl86m/7WoRbow1oNszErkKO1LKn1aZkntVPqzvcfIODPiPR3vOT5aecJ22VQ/inm6oH3oXBQnahffiD03vx1DsTvQlCeCWnID6n+lVC+ZNHcTISD3I3bwpSb1cgop39K5TDmcvS8dv8F6ZVwbqbnbJxau9/mFFrw3zynlPUgHpEN1eahnOeOSISrxXLrZ6ldTq65OawKcMAOS02sws+prVWMYz25vaLVpxZmuXmJLlpmo4Ye5jn1K8HcA8MxBq3+ZD4SS5treae1Yl86p/7W0ebuw1oNczLLnaO1Jqn1aZkntZP9zncMrMjjXf41RCmfzswf8yLb3g77zzOUGwfp32PjF1iXR8Scptg6x/K1zGPz3GK3kbt8O2Guzcbhx+xKYCL5yddEn6OOxzvRC4fqB98TWK32vQ1uWfyTkEOI3htpHbKrWv/F6hVJ2ZvnYvpNZOqmxar125xCC/6b55SyHssD5UK9XJuHspw7UrRaKpugX8j1z1K7LFwTOIT6eEmGI7ZZamIVfxZtrYIb6ynkFa0+7pkl8Ao5x4g8WmajlB7mWfRLYF7r+SxEf2/MQ33uwbchvL3EO8aRRf+xcqt5oT6sPWdmYZaQY619RuyhZZ5kJdXgOwQQkH/iJXPt9Azl3uLgd8Ncysng50+yMJmeYPliG46Xu+HfYy53y19h/gSv3/GSgXqZguu8KVb+L2L2xm4NwufY8XFtZ+z2xPyXmg3WHjHIB0DiM5k0HpZnPmmOO1WY808ivxhWe31vLd1V/xyQQ6jep5wS41zjErPdjf8igg/2XEQbYxFPnovqN048V7PfRi+Mcxf+G4OJmRt7IpgH2o/xcm0eWj13xGg2ljHWTqqN0S/0+mep3SrXRDYxHGJ8fJUY52iF1LmlJqmxjOVXtR0PiJln4F3CK1p9amAWwyv0HCPW0DKLsVGuMqv6Jfo2hnnN57MQ/T0xj/E5ryFpvXFV/5RqE/vrUtMl+rD2nLnKLDHvmBxr7jMh56klT6i3aQffZVB2727hxUYh/GlAF3MZ/JQBzfNztLFN7i6WbfKMnekkn+5+NW7AcX/gJYPvD/ASM8TUeYVy0pb27uex+aT5EKsbdknJRBROyX+puaG+K8zVfsKx8iHQW8yfLtW5sG31RLZwbNZNiFm8XrV/SucQofdJw5Q4jU3gxn+heQ0MrzBX909pA8ebnptRXxX9JiXOUG02jq/Wb/OcBp5XmB/qv3lcIeuWnojksfn+ayWXZjy0kp9qs6V20mCMfijj8v1OCpsYDsAX4+OrlDhVJtEd1E1/suZdyCuH6mPJLIYXyoSeY8T1hzLTdTv9USkaxDBHZDWfzwRssv6lmUf6/ColTr0Dd49M5r3bQkUHWGsy1Cdah/yvEdqHkzVMyTsmR5SJuTZcpcRpbMNk5pp47msOMjrmGeo5D5wPdZ4GQwexTpuwPA7Iyyc+00kGUx5NN2BZVaeUQb1SVj45lkEeedU0qfIMZDfP/wM21MRFxWSepHJdOAiPViYVK+f+UeXQiGAt+E+tV8K5Ofc5S5VDYr/xYNkW/DbnqNJOCjn23zynlHUVjwQv1+ah3OeOFK2Wyqr0Wyqo2GapXW6uKg4JPlbgKnKIpSZWAefW1irOsZ6cXtHqUxMzFa8RbsRcyyyi6mxFcuunYt7A+UwE0urvgnk2R5WrWMu7XES6lnLrr4tCf1Tf4w8nAAAbAklEQVTOPqzVMDczVY56ZK6P1DJPSuKTpNJhhWVQ/e2syPkXenFxkYHx/8FLEl+axFzyicR0UtWJuj/iJW2H3P08befoZVWeCDKE3Twn+bBD2qllUjER7SMSEp/NP/yJqMZNERUrROvZP6ocIvV2I9QQSAv+U+uVcG7Ofc5S5ZDYbzx4rwW/zTmqtJPzhWP/zXNKWVfxQAOx14DaPJT73JGi1VJZlX6R1z9L7XJzVXFI8PES+yO2WWpiFX9uba3iHOvJ6RWtPjUxU/GKPMeIJlpmo34e5rn1UzEHiNjrsgeGYwxa/V0wT/D5mO/Rcy3vo+Oct59b/3l7qes5+7BWw9zMVDk20GfEC1rmSb65n1Q6rLCY4zxBpEdYkV/oHQfkH+4IJ0DmF6DUOs/xOF9IzXOJ3Txl0eEf840L618M26TOI6dUJluxy2NeRl9uHVfLvlRWlv6JZZaaQ2y7luW0facF/5XQS3vOitUwNQdNv4mNTVOuJ7/NeaRqN69vaV3rP60OS21YbUvlsefl2s5ZNWknHkjVb8tHltppuW7Fs7UvlcOej7fattinPRdYamIRt9Sh1Vabo1Vca/Xk9IpWn5qYpfJa02HcrmU2Hu9hrtUvNtZU5kefzyRvbX/X6u+deazWVuWseVvFZVWPVn8tB6u41urJ2YfZZ9aoh2/X+kXLPDyCSYmSd76/QLu/YoBdnnn0Hi8xrCT5Etu+x/yfeMn08WZ2568MvF+YHOvaOu9UVtkGbZ4h7OYIhP/r+cZxHRqN++T5VjK9xjbRQ57Ff36O/2lPmT9aJjHRSI7PYwo6LaNllc0/Bly0ORg0ZVtFRN9pwX8l9No8ZxmoqM0hpd8YhHlZRad+u4RQ5r3Bpv8idJjnYLme28u1nbNq0k58oNUvxjOW2m1yjQluVkbLgefkGTiD1U1tnZ3vJN2cXtH2mZqYaXnFWknLLLb+HOU29TNoUMvc1flM8o7o71r9vTA3kNeuioy87YK0qWlT/wgONlGt15KzD7PPrHNX7Ynwi5a5qv21g+5dX18/ws43eMmd5/PB7bVyWbcjjv+gga8wl0/AThOW/w8LyTGiHnnOvAz4P7ypua2/yCuaHcr+BRrPMD9zb4kO8trVHsfI12vkA4Um/bGnJ/Juxj8avfd4lN7fs/9i9EIZF+csxBHdb0p7bNpez36bcpBlsNi9PiyUceG/eVwp66FertVDiLs57aa6a/ycQzsvXBEHz8lTQxgse9HWIJWLKkK8gmOD/kdolZkARG6qa2YoswtxDl7xoh/iqPJ8Nvik6j6j9fnBVj03X3N/G/zS5HuzkD4cqiGOd8UM8aiuDWfTOlgIZR4aMuqXx6u/wvxeycfOhMT5Ew6WTx9OEwKV5d8wd/HhwBCW11kKu5dIqqU7vmM0kk8xhUOvE/1zrPK9+y+UvpdzVkq/Cc3Z8nj6LY2mF/+lZXFZOtTLtXqoRe0uldxfy6GdF66hPt6nVeaIHJpYRe5FW6t8xnpCvBKqT6vMRnaaeSgzTZ2ljvGiX4hHS7HRthOqvxfm2vy8HRfK21v8reof0odDNWyVWUlvhjKPjs3lne+SDQbaxUgyySNq5C7kF9i29tUrOU41oY7qPo1RJTY5KIUdyv6Oqp5j3swHHchFPnV/jpd8iCPf9PgZr/fYLvPzNBwnn0rJDwF3OyH/6L6Hsof7Z9BxV29vAg9xd+s/5B91bvbgOfES4ojuN0d4EfHKebFbv82Zg0fV/pvnk7Ku9XLtHkL8h1+vUnRaKjtosnv9y6mdF66Ig+fkJZMkbPOibUIKi0U1XsExUdfM1pgNHHbPMQI6ltmiSAdt9KIf4qjqfJaivwfmg3dVPj/ImneaHWKu/n29B/3vwDXYgLx2+3Cshh6YDbFX1WdE1ljmIZZAG+c73z8JKVjyWAQpn0BwiiCQyO4ZmpTnuzczAA0e8kGCxk+vcJycNLqeavdPgN7edKb/4hRxcc5K7DdxmaeVot/S+I2lXfhvDMZiHuDl2j3UonYe3u+44BrgY4tuY1FHDf3JhbYWsKd1KL0Sq09TzMBKe44RxLHMpvIcvexCP6VHj2Y1bz9W/8OZB/p8nvdR67G8j4p3rd3D9V8LLGW7sg/Hang4s0r7jEgayzzKDm7vfI/KZqMQDKG6+3mjim52DayeYn5xZ3jLAJCr/OgvH21kIHKP/knF1rP/Br8kfVJOz4U5sGe/zUnRf3MiuvVWPNTjuaOEdj1y1fWc5aNKaLLccvjWHrVN1YfMwn3mqUSP+qXyZ59JJRhWPpV3WGv5j+6xz6Vq2COzVCemMte2j3bOd753M/iuhcPjSIAESIAESIAESIAESIAESIAESIAESIAESIAESIAESCCGwHTw3e1jZ2ISYxkSiCWATvEpyv4wlJdvScj0LbYn/87ATVX8SwLrBOi/dTbck58A/Zefccst0D/1qkvtfGpHXXzqIlFRmzzakGserl5qpb5elLiMg7pc8vC0Rm08qXEbS6ouHHy/Zcmlvgm8RGc6P+8dy/L8pzd4yY/9ciKB3ATov9yEWf8WAfpviw737RGgf/YI+d1P7XxqQ1186iJRUZs82pBrHq5eaqW+XpS4jIO6XPLwtEZtPKlxG0uSLvdv6+ESCXRN4DsMuD+ZEJBfpH6AbY8m27hIArkI0H+5yLJeDQH6T0OJx6wRoH/WyPjfTu18akRdfOoiUVGbPNqQax6uXmqlvl6UuIyDulzy8LRGbTypcRtLki6H3PmOAU3zR3zkqPOWsZ+lXvI8gLjc9f6/B7RbtEn6pyjukMbovxBaPNaaAP1nTbSv+rrwT6OSUjufwlIXn7pIVNQmjzbkmoerl1qprxclLuOgLpc8PK1RG09q3MaSpMshg++IPel2/dvcL5Zy1HnRgJOVXvIsihuD0r/MGpSO9Q7b3862175K/zhUkP7j452OtCX9R/+l+K8j/6RgclmW2rmU5Yq6+NRFoqI2ebQh1zxcvdRKfb0ocRkHdbnk4WmN2nhS4zaWVF2OeuxM0u36t+lfLOWo86IBJysmecI444+KOkkrbxgh+eJYedTMU7we543qkNrpnwOw039n6PTfGUW5BfrvzJr+O6PQL7Tkn5Bc9IT8HhmSL46Nfu8T0o5fWuUiC+GVootFRiGxWrR3dB0h+aZoE9LO0Uws2g/JN4Vr6Vgt2qu9jhBtJddYfUPbqZ1ravyhvGJ1SY1zLB8a71iu1nlIvrHahLRRK0fLuEN5xehy1OB70u36K5Bz1LnS1KGbk/OEUb5HBr09y1ye3y55b05Dp5PnvT/G8sfNg+vcSf8coxv9d8Od/qP/jiFA/6Vwb+L8xfc+6xZIee/TKdd1mLo9VfSpTrXNrg25rneSlHPReq36PZ1qowe0fKSqz0jRWH2pyzL4na3ZddlpX727U31V+rDPqG1kcaBKE2koVpd719fXMgj7Bq+HqOSdVFZ6Qrsy2PkU84dWbeeo0yo2y3pC88Txckf3PzB/sRQHtosfXuPV3OAzcpPB94+Yzx8xc0KB7fJtgBeYywDh1bAu80P6xSmozH+QW1Dfw/Gb/skcruvqwWaz72A//TdTkP6bAUlYpf/C4Vn7b0+D8Aj9lEBuVZ+/EP/mtYvaxb332ePqx8HlI9nzFPa77lN72u7lV564XYs5tSFXv/+H7Wlj57D6atrr73t9RjLGMVH/Z6Pc5vW7Ppp2EYNN0v+eKbpYZbGn716OVnEcUQ9yy/I+YI/pEbl6aXPPT3uaSB44JuhchuO/Q7FXmN87fPAdQWyeNGKEylFnTBy5y4TmiePlh27/xPzicSrD9l+x7wNeX+IlmnyG7c3d+Y2c5IOmr+a5YV060Su8ph9KyCC8DMY3xwF5XSGvoL6H4xf9I3X1Og1M1H0Hx9N/g1nov/ReQ//FM7TyX6gG8REfXxK5Vnn+GjTiex/j9z5rXI936nERhJ4PvPapNW1D8ztOifSWc2hDrqf/PVxeR9a0SXdSvTWE9ve1PiMEsC/q/+whhjvX73qppkc+MEn+3zNFl/QsbmpY0zc0R6t4jqgHuZqeE9eYHpGblzZD/bSmieSDfcHnMpQ5D74f9YOrJy2G4E0f8ZGjTi/GmcYRmaewlgHmiwl1yeDyM9mIZfkETgZlW50kf+Fwurt9kqSc+GRwWebnCTzmx5331byAvOTEEdr3Fv1TM4fU2MExtO/Qf4BO/6U676Y8/RfH0dJ/ERrEBe2jVK3nr8VrF7U7mSrlvc8iVx9WPSaKCE957VOL2kbkd4wQNq3m0IZcb/4PFQ7z/69SzkUWii9qY1FxrXVE9Pe1PiMIYvWlLjMDOdFlFlX06qK+ETlGB+Cg4Fq/YZ8xEifCT2uaSESxupyyOWzwHRBk8E/uKv5aIhnWZR79iI+hDtM6T5Sc/UnI8xuUnb/ZcZZd3nCQ/y94/Qevizvasf5Z3pb91I5cY/te9/5JVRHs6T/6L9VG0eXpv/N7jZj3Cd2f/yr2D7XLc+3pnmv0yXgo6LhPda9tJm3INc+5KLUrSvnutUmFuNZnpF7si/0/m7okCpNJl8SozsW713dNH/aZs0eKL6xpIoEk6HLK45DBdwR9vl0fy+Nd1jIoPH3kxylA7Z8cdWrbLnlcbJ4oJ89Li/5go2SOBdoSDt/gtfjs9wLtH9YE/XMY+mnD9N/NB6/qcz/PX1P7JC/Tf/Rfiomq8g/PHRdSm2lHrhdcU1fMdEkNRMpT2wuKZtqQax6uF7VGrlCbSHDLxdhnlrkcvdVMF6tE2O8uSJroQ6YXTFNXTDSZB2Ey+A6h5XEd8iwb+cHUN1i/GNTEuvyA57eYyyMaZNq9XT9HnTdN+/pbME/5hsEfvrJfjiaCyXJF61uFg/C48On64X73RLDa7Xsr2WbzT0QOKyFWs5n+u7kGnAWDB/a+kUP/nWklL9B/jvyXrOasggLn09r8k+3cMUOfvFqZdtm4FuCQrJVxBd76VDZtjbldFfCKpTbkemsAS663tcYvZdOmgEfjs85T0lJb6mKnkaUuVlFl09cqQKmnUB+20icr00IsLOVLqctKk4sYTAbfUeMPEEO+xi13V8sPQJwHNbFNBuWfYn56pjiWr7Cs+epRjjqleW9TqTzlh1Tl+UU1TEFMIhL6C2WeRJTzWCSIlbLvLeWZ0z9BOSwFV9k2+i9cMPovnNlaCfpvjcz69pz+W281bk/u82lt/qF2tz6y1C4n19weviXiY8lSF4uMcmprEd+0jtxesdSGXG+Vs+R6W2v8Uk5tcns0Pus8JS21pS52GlnqYhVVTn2tYpR6SvRhK31yMy3BwlK7lLqsNLmIIXnwHYN58uiAfw21yqctHy5auLnD+O1s2+Zqjjo3GzxoZ+E85dsJc20Oyny92VgmKCf5yaNknmP58XoLpz3CQR59VPUUyyoy6Sz+ic0hUO/IlLMVo//C0Xrzn/xAkEzv8foCfox+ZNqplrJ/6L9w3ln8Fx7GdomE82mIn2vzT+vaSX5HvPfJwrWQh7c7Uvm93vpUFm2tsRbyiqU2rXM96lxkYa0s2hTyqEX+lnW47zMJuoR43JKpRV2WuljEI3Vk6XdWwUk9CV4JeV8tTVnpk41pAota+42VJqLveUoefEdN7yDG26FG+Qfgp3PtNwtyh/F82+yQO6s56rzTiIMNJfP8HPl+tM4Z2kuH+hMvmWunZxPPzMsEMxlOBuOd7Jo45BlOmuPmsXlbD2aVkMCqfxI9EJxDhN6ntBPjTEB3pyj9dwfJ7gZP/pNvc72Hn36WqDF/gtfveMmHz6sT9lufK1fb2tlB/+0AWti96r+FY9WbMngi5nwa6ufa/NOydo9gtqPe+6xyTfR1CQ+f+mhinOp+rjjQW59a1VaRy+ohGXiX8IqlNi1zPfJctOq5gB2r2iT6toRHT2kmxhmAavfQIn0mMd8YXUI93rIuuyZQHrDa75TlFw9L9Ma8zhivhL6vljat+s0m00Q2MSyC+01ijHP9UtatNLmIIXnwHYA+So2Yyz8AMrgwfeSM3F0s2+SZORcTjl/9RCihTmlLewfQRTxHrCTkucqudB5DDnt3mqvDimGCMvLhz1vMnyob2jwxKes4/LAYVhI0ypn6Z4gjygMxOaBMqN4nrVLiNBa7d/+ZnqdTdI3xH7wg/eer0ROo4w+8ZPD9AV5yoV6csO8jdkT1k8UK4zfSf07eJ1h7YqjvCvOQ92Ohfm7CP/Hd56akB+0QQ+i1sIh2KWyGsleY5/TwScSUOFP9MytfRJdZm8VXrXkP9V1hntMr7rXxwBUxuDwXWZg8he9Q9grznB49pZkSpwWnSR1F+kxKvkPZK8zVuuDYUI93qcvEB4ctDvqa/L811HWFudorSDz0fbWwYr8ZHGOp31Bl7CyLJvdjo1koJ890l09EZIBhnE6DodgmJ6zzhPXxEyF5TrzcPSgDF7+fD7hdCKnzEYpJvTKwI6+appA8tezm+X/Ahpq4qJnME1WsCwfh0cqkZhXQ9+ZscvtHncM8sArXe/Zf7Hnahf/Qf8YPlOeD7B/hQ8mthon+C3+fkNt/1r5RnU8j/Vybf5rULtIwltrl5prTw5H4shWz1MUiyNzaWsQ4rSOnVyy1aZLrVIiAZUuuAc2uHppbm5weXU3qoB2W2rrQ5SCO1s1a6mIVW259reKUenL3YSt9SjBVsbCEf1BdVppchG85+C6DEm8var953vvprncZ9MNLkpBJPhE63w2P7bIsX92XOqaTuk6UlTufZSB/Pigyrc/rsjpPJKBlN89VuMz5zo/xtB7CJDRu+SSrRp+s5RnCyqt/QnJY41DL9m79l3Cezn3+0vpvvIbNvSZvdkTXGib6L/z8n9t/1r7J6efa/NOqdjGesdQuN9ecHo5hl7OMpS4WcebW1iLGaR05vWKpTatcp1poly25atvcOi63Njk9upXXEfsstfWiyxEcrdu01MUqttz6WsUp9eTuw1b6lGCqZWHJ/4i6rDS5iP3+xVraioh9njDI8ggr8ou744D8Q2z7iJcIJgMYF8dj/SNeUmY6XRyzVue0QKXLqjwD2c1RiA7/mG9cWP9i2CaGO3JSMYkMUL6KNPoysgpXxVSsCvknFowqh9jKC5XT9p0u/Zeogfb8FdtMqv/kfLk2MB8bU2g5+g/EMr1P0PpPq0GotqHH5/Rzbeev3rTb8oqldlquW/Fs7cvp4a12LfdpzweWuljEr9VWm59FTFt15PSKpTa9cd3SzJLrVjvafVpttPXNj8vp0Xlbuda1/d1SW++65GIdUu8RuoTEt3WsVl9tjlttpe7L3Yet+o2WaQqPVBYpbVuU1frJSpOLmD+5WEtbeYHiv+IfX7mz9j1eIowE/RLbvsf8n3jJtDZA8QH75gO+2jpPFVf8R5tnCLs5DuH/er5xXIdG4z55ppVMr7FNNJRHAp2f43/aU+aPlklMNJLj85iCTstoWWXzjwEXbQ4GTdlWEdF3evVfCvjN81dKxUNZrf8+rrQlfevizcjKceab6b8rrXYp7Df9F6FBSiyaslomMX6u7fzVqnYaH8yPsdRuk+u84Yj1nB6OCEdfJOJ8YKmLPtD1Ize1jchvvSWbPTm9YqlNq1xjVLTkGtP+vMymNvODI9ZzejQiHH2RiP5uqa0XXfTACh15sC5WWW7qG5GjVVxL9eTuw1b9ZpPpUmIR27QsIqrOVyTCT1aaXCR17/r6Wu42f4OX3JmeffAAbYztfYbl8z9+WP4PYvgJc3l0TPSE8vKceRnwfxhdidOCyCmJHcr/hdSeYS6fijU3Ia9d7XGMfPNCPlBozh97giLnpvyj0XuPSen9PftvyjpGO5Rxcf5CHHKt+grz83kUy/+HbUWuoVOOocuIs9vz35QVOOxeK6bHyzLKuPDfPK7UdeSl9jOOrdI/rWo3ao/8dv2cQzsvXBGH2sMjMy/zHLpY5OZFW4tcpnWEeCWHNq1yHRkjv0PORWP7KXMv2iAOns8mQnrRZQxJ4/HxWC9zxOz2vZs3fS00Q05BfdhaH49MEdPutcGCvVUdGTSR3+x8hXrv3bcKMqCe84D7rMxhdw7O4vC8msruJZJr6Y7vGK3k0zrh0ONE/xyves/+S6Xv5fz1ExKRT8NPEy6ksvwb5tk/vB6aTJnRf/H0vPgvPoPlkiF+rtU/rWq3rOjy1hzaeeEa4uFlOsdtzaGLRTZetLXIZVpHiFdyaNMq1ynjveUcXPfa1Oz3ok2IRzV5lTwmh7ZedCnJ0bqtHLpYxdiivqF92FqfFpla+U1bj7Um53aLD74PAxQyCCifws2n852E8x1cP915J4M70ezAXh4f8wDzJfbVIpZ88JITzQ9DfuOjji5yGvKWY494jM5FLEesIO8m/CM64rWr9xGMt9qUuLG/W/9tsdHsAz8X5y/EId/O+kI8iJc8Uk1+9f1bTQ5HHoNY6b8EAcDPhf8SUlgsirxUfq7ZPw1rp7oW5tLOC1fEofLwYgc4cGMuXSxS8qKtRS7TOrReyaVNw1wPPRdNNY5d9qKN1qOxeeYq13qfkfzw4v+exgYC0+beWyMn9XsS8RWQmo4NeGIq+eFVVb/Jocm023wyXSm4PH4idBpsR5I13TlYENNiU6nsZKBInu/+9WLtFW6Ef2RQWT6h2pte4YDnewc1vr96/wTo7U1K+i9dERfnL3hQc75Jz9a2BvovnacL/6WncVmD0s+1+6c57QKuhTm1c8FV6eFL4x+/llMXi+xcaGuRyLQOpVdyatMcVyfnoqnMscsutFF6NDbHXOWa7jMBHs/FN7benLrExjQv56LfzYNKWQ/ow7n0ccG00n6TS5OTpYo/8300MsSQT0Fkkh9nledvv8C2tcdiyHGbE8rKJ0cysCoD+fJsa/nU6T22y7ypKZXdwOppi2zWhEaucodqLY+GWEvDZDv9Y4IxqBL67wYXOCSfp4c6ujp/BZlt4WD6j/5bsIV6Uyv+6fHcUUK7HrmqO8/KgSV0WWk6aHOP2pbQhlyDbFj04B61SQXMPpNKME/5ErpYRd5jv8utT49MU/2YSxPUe37m+2GD76lwWJ4ESIAESIAESIAESIAESIAESIAESIAESIAESIAESIAEPBGYDr5PHzvzF3bM45Q7heVrC5xIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgAQGAhg7/wuL8k3/xUkG3+V52WvPwZZ9nEiABEiABEiABEiABEiABEiABEiABEiABEiABEiABEiABC4JjI9Wv9w6rP0/IxjgA/F2g/wAAAAASUVORK5CYII=\n", "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\\\- u_{0} & - u_{0} & - u_{0} & - u_{0} - 1 & 1 - u_{0} & - u_{0} - 1 & 1 - u_{0} & - u_{0} - 1 & 1 - u_{0}\\\\- u_{1} & 1 - u_{1} & - u_{1} - 1 & - u_{1} & - u_{1} & 1 - u_{1} & 1 - u_{1} & - u_{1} - 1 & - u_{1} - 1\\\\u_{0}^{2} & u_{0}^{2} & u_{0}^{2} & \\left(- u_{0} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2} & \\left(- u_{0} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2} & \\left(- u_{0} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2}\\\\u_{1}^{2} & \\left(1 - u_{1}\\right)^{2} & \\left(- u_{1} - 1\\right)^{2} & u_{1}^{2} & u_{1}^{2} & \\left(1 - u_{1}\\right)^{2} & \\left(1 - u_{1}\\right)^{2} & \\left(- u_{1} - 1\\right)^{2} & \\left(- u_{1} - 1\\right)^{2}\\\\u_{0} u_{1} & - u_{0} \\left(1 - u_{1}\\right) & - u_{0} \\left(- u_{1} - 1\\right) & - u_{1} \\left(- u_{0} - 1\\right) & - u_{1} \\left(1 - u_{0}\\right) & \\left(1 - u_{1}\\right) \\left(- u_{0} - 1\\right) & \\left(1 - u_{0}\\right) \\left(1 - u_{1}\\right) & \\left(- u_{0} - 1\\right) \\left(- u_{1} - 1\\right) & \\left(1 - u_{0}\\right) \\left(- u_{1} - 1\\right)\\\\- u_{0}^{2} u_{1} & u_{0}^{2} \\left(1 - u_{1}\\right) & u_{0}^{2} \\left(- u_{1} - 1\\right) & - u_{1} \\left(- u_{0} - 1\\right)^{2} & - u_{1} \\left(1 - u_{0}\\right)^{2} & \\left(1 - u_{1}\\right) \\left(- u_{0} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2} \\left(1 - u_{1}\\right) & \\left(- u_{0} - 1\\right)^{2} \\left(- u_{1} - 1\\right) & \\left(1 - u_{0}\\right)^{2} \\left(- u_{1} - 1\\right)\\\\- u_{0} u_{1}^{2} & - u_{0} \\left(1 - u_{1}\\right)^{2} & - u_{0} \\left(- u_{1} - 1\\right)^{2} & u_{1}^{2} \\left(- u_{0} - 1\\right) & u_{1}^{2} \\left(1 - u_{0}\\right) & \\left(1 - u_{1}\\right)^{2} \\left(- u_{0} - 1\\right) & \\left(1 - u_{0}\\right) \\left(1 - u_{1}\\right)^{2} & \\left(- u_{0} - 1\\right) \\left(- u_{1} - 1\\right)^{2} & \\left(1 - u_{0}\\right) \\left(- u_{1} - 1\\right)^{2}\\\\u_{0}^{2} u_{1}^{2} & u_{0}^{2} \\left(1 - u_{1}\\right)^{2} & u_{0}^{2} \\left(- u_{1} - 1\\right)^{2} & u_{1}^{2} \\left(- u_{0} - 1\\right)^{2} & u_{1}^{2} \\left(1 - u_{0}\\right)^{2} & \\left(1 - u_{1}\\right)^{2} \\left(- u_{0} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2} \\left(1 - u_{1}\\right)^{2} & \\left(- u_{0} - 1\\right)^{2} \\left(- u_{1} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2} \\left(- u_{1} - 1\\right)^{2}\\end{matrix}\\right]$" ], @@ -941,10 +941,10 @@ { "data": { "text/html": [ - "<table style=\"border:none\"><tr><th>Name</th><th>Runtime</th><th>Adds</th><th>Muls</th><th>Divs</th><th>Total</th></tr><tr><td>OriginalTerm</td><td>-</td> <td>1543</td> <td>1917</td> <td>2</td> <td>3462</td> </tr><tr><td>sympy_cse</td><td>155.96 ms</td> <td>248</td> <td>179</td> <td>1</td> <td>428</td> </tr></table>" + "<table style=\"border:none\"><tr><th>Name</th><th>Runtime</th><th>Adds</th><th>Muls</th><th>Divs</th><th>Total</th></tr><tr><td>OriginalTerm</td><td>-</td> <td>1543</td> <td>1917</td> <td>2</td> <td>3462</td> </tr><tr><td>sympy_cse</td><td>193.52 ms</td> <td>248</td> <td>179</td> <td>1</td> <td>428</td> </tr></table>" ], "text/plain": [ - "<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f9b738c1100>" + "<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f7cc7993190>" ] }, "execution_count": 19, diff --git a/lbmpy/moments.py b/lbmpy/moments.py index b1e5af5833c808a71dec6c4e1347926b62a752dc..cfb868a73f0255e3f3442b950b4f307f22c9e462 100644 --- a/lbmpy/moments.py +++ b/lbmpy/moments.py @@ -379,14 +379,13 @@ def shift_matrix(moments, stencil): directions = asarray(stencil) - u = sp.symbols("u_:{dim}".format(dim=dim)) - f = sp.symbols("f_:{dim}".format(dim=nr_directions)) - m = sp.symbols("m_:{dim}".format(dim=nr_directions)) + u = sp.symbols(f"u_:{dim}") + f = sp.symbols(f"f_:{nr_directions}") + m = sp.symbols(f"m_:{nr_directions}") Mf = moment_matrix(moments, stencil=stencil) * sp.Matrix(f) shift_matrix_list = [] - shift_equations = [] for nr_moment in range(len(moments)): exponent_central_moment = moments[nr_moment].as_powers_dict() diff --git a/lbmpy_tests/test_moments.py b/lbmpy_tests/test_moments.py index f6ae3cf27c8280386721e437d45469588ae4fb2f..565f436d2d104aaaaf08b97dd2a122e0b0224a9e 100644 --- a/lbmpy_tests/test_moments.py +++ b/lbmpy_tests/test_moments.py @@ -1,4 +1,5 @@ from lbmpy.moments import * +from lbmpy.session import create_lb_method from lbmpy.stencils import get_stencil @@ -100,3 +101,24 @@ def test_is_shear_moment(): assert is_shear_moment(x * y, 2) assert is_shear_moment(x * y - 1, 2) assert is_shear_moment(x * y - x, 2) + + +def test_shift_matrix(): + stencils = (get_stencil("D2Q9", ordering='walberla'), + get_stencil("D3Q27", ordering='walberla'), + get_stencil("D3Q19", ordering='walberla')) + for stencil in stencils: + method = create_lb_method(stencil=stencil, method='mrt_raw', compressible=True, equilibrium_order=6) + moments = method.moments + + sm = shift_matrix(moments, stencil) + + m_eq = sp.Matrix(method.moment_equilibrium_values) + m_eq_rel = sp.simplify(sm * m_eq) + + m_eq_sol = sp.Matrix(method.moment_equilibrium_values) + + for i in range(0, len(stencil)): + m_eq_sol[i] = m_eq_sol[i].subs({sp.Symbol("u_0"):0, sp.Symbol("u_1"):0, sp.Symbol("u_2"):0}) + + assert m_eq_rel == m_eq_sol