{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 毕达哥拉斯三体问题" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 问题的陈述" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "“勾三股四弦五”是毕达哥拉斯定理的一个特例。人们在研究三体问题的时候,尝试把三个星体摆放到边长分别是 3、4、5 的直角三角形的顶点之上,三个星体的质量分别是对边的长度,初始速度都为零,然后看系统轨道的演化,出乎人们的意料,轨道呈现出非常复杂的图案。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 定义系统的演化函数" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import wq.core.physics.unit.au as au\n", "\n", "from math import sqrt\n", "from wq.core.math.ode import verlet as solver\n", "from wq.core.physics.nbody.body3p import accelerationOf\n", "\n", "accel = accelerationOf(au, 5.0, 3.0, 4.0)\n", "step = solver(accel)\n", "\n", "time = 0\n", "x = np.array([0.0, 0.0, 0.0, 4.0, 3.0, 0.0])\n", "v = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0])\n", "\n", "cx = 4.0 * 3.0 / (5.0 + 3.0 + 4.0)\n", "cy = 3.0 * 4.0 / (5.0 + 3.0 + 4.0)\n", "\n", "K = 0.0\n", "U = - au.G * (5.0 * 4.0 / 3.0 + 3.0 * 5.0 / 4.0 + 3.0 * 4.0 / 5.0)\n", "E = K + U\n", "\n", "def evolve(tao):\n", " global time, x, v, K, U, E\n", "\n", " time, x, v = step(time, x, v, tao)\n", " x1 = x[0]\n", " y1 = x[1]\n", " x2 = x[2]\n", " y2 = x[3]\n", " x3 = x[4]\n", " y3 = x[5]\n", " vx1 = v[0]\n", " vy1 = v[1]\n", " vx2 = v[2]\n", " vy2 = v[3]\n", " vx3 = v[4]\n", " vy3 = v[5]\n", "\n", " r12 = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))\n", " r13 = sqrt((x1 - x3) * (x1 - x3) + (y1 - y3) * (y1 - y3))\n", " r23 = sqrt((x2 - x3) * (x2 - x3) + (y2 - y3) * (y2 - y3))\n", "\n", " v1 = sqrt(vx1 * vx1 + vy1 * vy1)\n", " v2 = sqrt(vx2 * vx2 + vy2 * vy2)\n", " v3 = sqrt(vx3 * vx3 + vy3 * vy3)\n", "\n", " K = (5.0 * (vx1 * vx1 + vy1 * vy1) + 3.0 * (vx2 * vx2 + vy2 * vy2) + 4.0 * (vx3 * vx3 + vy3 * vy3)) /2\n", " U = - au.G * (5.0 * 3.0 / r12 + 5.0 * 4.0 / r13 + 3.0 * 4.0 / r23)\n", " E = K + U\n", "\n", " dx = (5.0 * x1 + 3.0 * x2 + 4.0 * x3) / (5.0 + 3.0 + 4.0) - cx\n", " dy = (5.0 * y1 + 3.0 * y2 + 4.0 * y3) / (5.0 + 3.0 + 4.0) - cy\n", "\n", " x[0] = x[0] - dx\n", " x[1] = x[1] - dy\n", " x[2] = x[2] - dx\n", " x[3] = x[3] - dy\n", " x[4] = x[4] - dx\n", " x[5] = x[5] - dy\n", "\n", " return min(r12, r13, r23), max(v1, v2, v3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 绘制系统演化的轨迹图" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXecFEX2wL9vYZccRSQKBlDPjIqIoKuiIqfgKeepZw5n\nAgOcEb1lPVHPcGLWE1E4T0yY+BlRQTGhKAgKqKgIKCAoeUFg9/3+6N7dCT0zPTPdPT2z9f185rPT\n3dWvamdeval+9eqVqCoGg8FgKFyKct0Ag8FgMPiLMfQGg8FQ4BhDbzAYDAWOMfQGg8FQ4BhDbzAY\nDAWOMfQGg8FQ4ARm6EVkoYjMFpGZIvJJUPUaDJkgIp1FZIqIfCUiX4rIpQnK3SMi34rIFyKyb9Dt\nNBjcUD/AuhQoVdXfAqzTYMiULcAVqjpLRJoCn4nIZFWdV11ARAYAO6tqNxE5EHgQ6JWj9hoMCQna\ndSMB12cwZISqLlPVWfb79cA8oENMsYHAOLvMdKCliGwXaEMNBhcEaegVeEtEZojI+QHWazBkhYh0\nBfYFpsdc6ggsjjheAnQKplUGg3uCdN0crKpLRWRbYLKIzFfVaQHWbzCkje22eQ64zB7ZxxWJOTY5\nRQyhIzBDr6pL7b8rROQFoCdQY+hFxHQQg6+oalquQxEpBiYCT6jqiw5FfgI6Rxx3ss9FyjB6bfCd\nVLodiOtGRBqLSDP7fRPgKGBObDlV9fRVVlbmuUwjN//aqpq+rRURAR4F5qrq6ATFXgbOsMv3Alar\n6nK/9TofP3+jg7nV7aBG9NsBL1h9h/rA/1T1zYDqNhgy4WDgNGC2iMy0z10HbA+gqg+r6qsiMkBE\nFgAbgLNz09TsWLRyESc+eyIzfpkBU6G8vDxh2YY0pG+Xvjz0x4fYcdsdg2ukISsCMfSq+gOwTxB1\nGQxeoKrv4+KJV1WHBNAcT7h00qXc+/m9WcnYxCYm/ziZnR7YKep8McW8dMpLHNP9mKzkR9GiBaxd\na713OXI1OBPkZGzglJaWGrk+yc2nthYibj6n0v+U8u7Sd9MT3DWj5rCFLQyYMKDm+M1T3+TIbkfW\ntiVRe7t0gUWLUlcgEmfs800Hc6nb4tbH4zciomFpi6HwEBE0zclYj+oNVK+bljdlAxvSvm/nZjsz\n9rix9O3WN2m5dZvWcfmrl/P0V0+zoSp1PTceciM3HHZD9MkePWDmTOcb3PLgg3DhhdnJKBDc6HZg\nhl5E6gEzgCWqepzDdWPoDb5RyIb+7a/fpt9T/VyVfeCoB7jooIs8b8PspbM58JED2aSbHK9rmcKw\nYXDXXZ7XXdfdOmEz9MOA/YBmqjrQ4box9AbfKERDf/7z5zNmzpikZZ49/lkG7z3Yl/oT8fyc5znx\n+ROjTyosvBW6/O5ww++/Q0lJYoHi8muro/bDjW4HFV7ZCRgAjMGkQTDkASIyVkSWi0hcGLB9vVRE\n1thJ+maKyPVBtW3JyiVIuSQ08l+e9yVapmiZBm7kAU7Y8wS0TBmxz+W1y8cEul6T4IYGDSxjXv16\n6aXkFZx3nvP56vsNcQQyoheRZ4GbgebA343rJh4ptxRUy+ruZ+An6Y7oRaQvsB4Yr6p7OlwvBYY5\nPZ3GlPNUr6v1xInQ6M7UqXDYYTWH8g9qh5RVoDemIeu116B//3gDXv2ZJjLs++yT/TxAnhCKEb2I\nHAv8oqozMaN5RyI7b7KObAgOtdJzrEpRLNAvK5FuVI/eQ4FIlJEH27BHjOyPOg3YbTeoVy+1vGOO\nST5KV3V22cyaZd1XWem25QVNEOGVvYGBdkrXhkBzERmvqmfEFhw5cmTN+9LSUhNqZ8iYqVOnMnXq\nVD+rUKC3iHyBlfbg76o614+KqqqqqPfPeKMYGuNeTRKDfO1rcMsAQGDyzsB/E3xURx0Fkye7q0sV\nxo6Fc89NXK6+beK++gr+8IfUcguUQMMrReRQjOvGkcjR2usnvs7Rexydw9YUHplMxtpZKyclcN00\nAypVtUJEjgHuVtXuDuWy1uvYkfxOxTux4LoFWcn0nFgj36CBNclajWrU/+HqR2qvvWCO4xRJdjRu\nDBvSD0ENK250OxcLpuquNXdJ/4n90T3MxxRmVHVdxPvXROQBEWmtDhvrZPOkGmvk1w1fR9OmTTNo\nsY/Etuexx+Css7KfGJ092/rr9QRrRUWtzC1bakf9eUImT6tmwVSISHvEY3CNDyP67bDmnlREegLP\nqGpXh3IZ63WskQ+lTmzdCsXFtcdjx8LZdsqfGAMtI2vfp/xf0jXuTZrAeocs0v36wdtvp74/j21P\nKCZjDZlhJmVzi4hMAD4EdhGRxSJyjohcICIX2EUGA3NEZBYwGjg5V23NKZFGvri41sgDnH56dNlk\n86JnnBEdYpkuGzZY93XtGn3+rbfg449T3y8CoxMlKc1/zIg+REyaMYmBr9RG64VyBJen5NuCqbwY\nzUO0UXb6PyOu14RZKmjiBJnORMr+7TfYZpvk5XffHb78svZ4y5bki7KqadkSVqUKtgoXoRnRi0hD\nEZkuIrNEZK6I3BJEvfnGcftHz1GbUb0BYO2wtbluQmImT4aiIhg61Pl6pIGuVme3v1nVoZOxPyCt\nW6e+96uvrB+Ztm2t4+LiqFBLoQJhq/1awxQOsC6sXl2Qi64CMfSqugk4TFX3AfYCDhORPkHUnW+0\nr9c+100whIxmzZrlugmJ6dfPMqD33JO4jCod/tW+xtA3c9qQEeC//01s3B1kRnHfffDNN/HlVqyo\ndQcNHQqPPopQiRXpXc9+NedwpiNUIVSyO5+i9fJrgjYVgfnoVbXCfluC9enGRScY4Ofrf446bv9P\nY/gN+c/SjUtr3q+9U6MNevXrtNMyr2DIEOjWzZLz88/OZR54gKvOfQ/rF0eIfrSoPlfEXPanqGoL\nIlUccEDmTQoTgRl6ESmyJ66WA1P8WlxSaCyrWpbrJhgMWdHh9g7+CI4d1Ve7XNq3r/3xaNw4qkj7\nqPFlFaowYkSU0GphQBEzZtQ+ELz3nrfND5LAJ2NFpAXwBnCNqk6NOK9lZWU15er6ythI/3xzmrOm\nbE0OW5N/xMYal5eXp5vrZizwR6wQyrjwSrvMPcAxQAVwlp3mI7ZM1pOxoZ2IdcGi1YvocneXmuOt\n12+lnpvUB2656Sa4ISLf/ezZsKfD1/XRR9C7NwBCFdUj+nfekdiMDRwtL/Amg6gd5UfTuzd88IFX\n/0D2hCpNcVSlIjcAG1X1johzdT7qJpJut3Vjwcba1Y/53NnDgA9JzQYAQ1R1gIgciLUytpdDucKO\nuknC4jWL2X709jXH+7bdl88v+tz7ihIlPMumuAj/ZRBnMBHL8ZFcdU44ASZOTNlSXwhT1E0bEWlp\nv28EHAnUjdRyGfLtVd9GHZsInGBxkdRsIDDOLjsdaGkvovKEI7c7MnWhELN+0/ooIw/4Y+QhsQsn\nw+KzZoFQyRm8gDWdmLrvPf98rYvnjLgsXrknKB99e+Ad20c/HWu1oYvlanWba3tem+smGBLTEVgc\ncbwE6OSV8DcvfDPqON9+6Fv9q1XU8e/XOe044iFvvRV9nCJmfuPG6OMWLay/gwfDvvtC7She7VcV\nxcXWAtxU/Pe/lsFv1Sp12aAIKrxyjqr2UNV9VHUvVb09iHrznZuPuTnqON86ex0g9gvx1L8yqveo\n6Mry5Ps/bOxhbGVrzfFDAx6ipNjFYqVsOOKI6InXLVvg8MMTFm/Y0MrWUM3atdCmTaz7pYrHOQml\nCC1uyObNVpYFN9GfUBuSv8suaf83nmNWxoac5auX0+7udjXH+eirDQM+5Lp5CJiqqk/Zx/OBQ1V1\neUy5rIIMnIx7mHXgj//9I69+/2rUuUDbG+uHGTYM7rwzYfHLL4e773a+phojL4F9cru+atEi6NzZ\nXdlkZBJoYAx9HlAIE3O5xgdDHzkZ2wsY7eVkbJSMPDH2fcb04YOfosNR2jZuy/Irlye4wydiLe/p\np8P48QmLn3ACvPBC9Lmar8yFoU9UrRN+ZEgOTdSNiHQGxgNtsR5v/6Oq98SUMYY+CZGdffN1mymO\nTCZlSEkGUTcTgEOBNlhrP8qAYgBVfdgucx/QH9gAnK2qcbONXul1st2lwkCXO7uwaP2iuPMV11bQ\nqKRR8A2Ktbrt2ydeSOVQfJttYOVK0jL0AJs3W6n4U+GlqQuToW8HtFPVWSLSFPgMOF5V50WUMYY+\nCWZUnx35ltTMid3v3J25653XGeZSHxqUN2Azmx2v5VRPnYbYabhfGjaEjZvSM/TVvPQSHH988jJe\nmbvQhFeq6jJVnWW/Xw/MA3xaLleYxHaYjb9vTFDSUKh8NfyrhIZTygUpF857/rxA2yTlktDIv3Xa\nW47nAyNRNs0b3e1OvmmTIlSkLujAoEFW9R07Ji4TZO60XKyM7Qq8C+xuG/3q82ZEnwIzqs+cQhjR\nRzJj4QwOGJc6EcuGKzfQOCYNgFfE6mN96kdF24RGP7fd1vbDxDBmTM1+s5FGt7QUpk5VasMrt6IL\nf4YuXeJluCTV/ubZEBrXTU1llttmKnCTqr4Yc82kQHBBoSyN95tsUyB4hd8DmIqKCprc7iK4G7hg\n3wt4aOBDWdf52IzHOOeVc6LOXbbfZdz9WW34ypLLl9CxRZLhbND8/HPS4XVtWoStKMUMZhwTOZ1a\nYy9ZG+Trr4dRo5yvZSM7VIZeRIqB/wNeU9W4rVzMiN4drkf1S5dCB5fesYMOgg8/zLJl4abQRvRO\nnPPcOTz21WNJy2Q7ONj2lm1ZuTl6dPz1RV+zy4O1weItiluw+rrVWdXjG336OCaqqTX0VShWLp4r\nuJHRXE/kcgkvvspEo/tMZYfGRy8iAjwKzHUy8gb3JO2oN95Yuw7brZEHK+FT9X1ulv4ZQsnYwWPR\nMq15tZAWUdeP3CG7tApSLnFGXsuU3R7cLepcaI08wPvvWxb1H/9IUKDWXt7V/QnuvTdm31sPhgqq\nsNtu8ef99NkHFXXTB3gPmE3t6sFrVfX1iDKhHNF/ufRLDnvsMFZucfDxpUH3Vt159s/Pslf7vbJu\nU+So/qUTX2LgnoOylulICL+PTMkwjr4/1n6w9YAxqvqvmOulwEvA9/apiap6U0yZUOp1Oixdu5QO\nd8UPHLRM6XZXNxasrU2+9+NlP7J9y+3jyoYdkUqsr1mJVZPJk+Goo6LLexcxE38uXdmhct2kIgwd\n4trXruXWT24NpK7OTTuzaHh83LEbznzuTMZ/ZS8A2QrRpiWGVJ/puedGrwWP5fjj41eT5CEZxNHX\nA74G+gE/AZ8Cp8SEBJcCw1R1oKMQwqHX2dDl9i4sqojW0x5tevDZJZ8x6IlBvPzdyzXnB+08iBf/\n+mKsiLxg0CB42f5XDjsM3nkn+vqiRfFzsRs3WiGY2ZKtsTeGPgWVlZXUvym3W4YVUURlWWXqgtX8\n/js0bIiUUTNPFLfRcqaf48svWxrvREj0JFMyMPQHAWWq2t8+vgZAVW+NKFMKDFfV4xyFkN+G3mmR\n1o9DfmT7bbZn8FODmfh1bWKYpvWbsm7EuiCb5zmp1kY5LYZavBg6ZZnKbscd4Ycf4s+7VRs3ul1Y\nGyO6JN3kUKWdSply7pTMK1y1itf2as0Jp8CmGBd4FVVIuaSeJDvxRCsXajXV0V/VPP00nHRS5m0E\nGDjQ0q6lS7mkw2M8wNXUTONIbIXOeKH4IcEpO+WBMWUU6C0iX2CN+v9eKDunJUu7sP8D+/PZis+i\nruW7kY9l1izYZ5/ocyUlVveI/EHo3Bk++wx69Mi8ru+/9z+mPhBD72a3nkDa4cLAjzx4JGX9ylKW\nS4vWrTkG2Hg7sGkTWlJC0Y3R8+COxn7NGmjZMrlswdPR9sSJMHhwe+C6jO6PTNp09dVwazCeMD9w\n86F+DnRW1QoROQZ4EegeW2jkyJE178MeNtzz/p58uvLTqHMllPB7mZVmuKS8hC1sibpeKGG+991n\nbT0LVqriRN0q1tjvtx9MnQqHHup7E4H40GE3BDUZm3S3HruMr4+4yYy874oa+3O9eTMUF3P7e7dz\n1ZSrotvRqBFs2pRa5D+wBtuxrpvu3a2R/+67e9LU2nzcVk5uO91LWtxyC1xzTUbN8YwMXDe9gJER\nrptrgarYCdmYe34A9lPV3yLO5Y3rxqmPTBg0gZP3OTnh9UIx8tVE6v8bb8RPwiYqC9aesn37Zl8v\nWBPA/fq5vTck4ZUuduvxFScFvemgm2rC0AKnpAREuPLQqyzbCbYrRlwZeSCxF+Wbb2CPPWrDJWNf\n22wDf/sb/PZb3K2RW29CFaM438rFTT2U+ujr79Tk4nZ6JVo4eO21wS739ogZQDcR6SoiJcBfgJcj\nC4jIdnboMCLSE2vgFP/B5gGJjPjJ+5zM0zOfjrteTHHBGXmw0hZXc/TRycvG/n4fcgjMn+9NO9wa\nebcEuWCqKwlSvtrXfRv5hCJ1QAJLl3RSNQE37A83/ZHq9R2ou9QdKRnP0ZzJazWCqxeO1JDG99Os\nmbVJQ4a3e06G4ZXHUBte+aiq3iIiF4CVwVJELgEuArZibRA+TFU/jpER+hF9spG607UjOh/BW+fk\nOI+Nj0R21R12sHzobsuDNVZzk8EykYxp06x1Xe7vzbPJWD98mduVR2/jmbNRSHVnb9XK2nrGZsiH\ncF9vGDw7pnz37vD113Fi1m5Yy0131C6EOWOvM0DHWQdbtliTtpMmZdTEM3iDM2uOBKGSCoppRFXa\nVnqdPTeXq5F8Jn7MWFT1NeC1mHMPR7y/H7g/q0pyTCIjn8jVuXr4alo0beF4rVD44gvYe2/r/Q8/\nWK8ddkhcPtZn37Bhet3l2pgdQ9Mx8m4xI/qw4CLv9QlPnsAL30bHtKf9v7z9NowcCbNnW9Y4pq4x\nDOZ8niHSN5TpZgmxRj7fRvQe1RvqEX06EWih7Dc+0bhx9L6ybr7CTPS9a1f48cf074uuN0Rx9GEy\n9JGERnkjtaRlS1hVO6Vx8csX8+DMB+Nu8avt06dDr7i9kqBePdi6Nf58JP/+Nwwf7nzNGPrw4cbQ\nr7tyHU0bNw2gNeEiXcNdVWX1ETf3PPUUnHJK/PlMVCU0hj5it55tgF+Af6jqYzFlchZ1E0kLWrC6\nLAe5OkpKLNcL1vxs/5Ng8h8SFw/iBypZAFB16OS6ddC8eWpZubZ1xtA7k6xfTDphEsfueWyArQkX\n69dbc03VFBVBZYq1jU8+CX/9a/S5RGnxY/EzqVmdWxmb7mKpRNxRegfDD00wdHXBuo3ruGXaLTw5\n+0l+3GA/u0X++wmaeVe/u7j84MudL/qA0yglHdav9ydP2t//bu357OXqQT8Iu6GHPHFt5ogrroDR\nEWkYmze3lrcko0MHK3lsNU5ZRLx0axpDn4Sqqirq/TMLCxYE9sdRtBl+HA2dZs6HXXZJfk+WbNkC\nF18M48bVPGCkzZ57WlMAfpJuR/EjqZld5h7gGKyom7NUdWbM9dAbekNyOneGJUtqj5s2rQ02SESQ\n81PG0KfJ6/Ne55hnjslpGwDYCqd/CuPfSFHu6KPh9deTl1G14rUmTrSc70uWwDbbUH7CF0QEOfnK\nhx9aKe+9JN1HX5+Smg0AhqjqABE5ELhbVXvFyMm5Xhuyp7g4fn4q9eAivfKZEipD7yLla150iMEP\nD+aVZa+wCZcLmxLQtKgpPdr3YGTpSA7b+TDnQr16WcbZY4QtWF9D+p6MCy6AhyI2KZo61cr2ly5/\n/Ss88UT690WSTkfyKanZQ8AUVX3aPp4PHKqqyyPK5IVeG1KT7uDiu+9g553dlc0GV7qtqr6/sKzK\nAqAr1hr6WcBuMWXUkIBOnZItSE37BVWeiWvcWLVvX9W337aauvPO3jX10ENVlyyJ/ijOPz9x+WTY\n+pWOzg4GHok4Pg24N6bMJKB3xPFbWCkQjF4XKMXF8Xq3YUPi8unoaKa40e2gFkz1BBao6kIAEXkK\nGATMS3ZTnWPJkuisYL6xFqgOlclufrKiwvIMHXFE1o2K4913c5oJ0+34K/YDjLsvn5KaGZKzeTPs\numv0WsYmTSzXzubN8eVVvV80GOakZoOBo1X1fPv4NOBAVR0aUUaDaEso2X57K79vJvTu7bgHZjqo\nWqGSL7wACxfCr79a0Tb5RCrV8SOpme26maqqT9nHxnVTR0i0dcNNN8GIEdHnXKyFzIrQ+OhF5ESg\nfypDX1ZWVnNPnRr5uPnJv/xyuOsu/9uSJnfeaaV3XbjQW7nFxcmjfpw6VCSxo57y8vJ0DX19rMnY\nI4CfgU9IPhnbCxitZjK2TpGs61Z/7XXJ0LsZHdXtDuG3NoSM//s/OC7Bvky56gwO9yRNamaXuQ/o\nD2wAzlbVz2Nk1G29rgP072+lNHZDoRt6N6Mj0yEMvmEWTBn8xs2Dea4MfVD56LcCQ4A3gLnA05FG\n3mAwGPKd6tiaJ59MfD1XmAVThjqBGdEbCpXQjOgNBoPBkDuMoTcYDIYCx3dDLyJ/FpGvRKRSRHr4\nXV8k2e4wZOQGK9NPuekgIq1FZLKIfCMib4pIywTlForIbBGZKSKfBNnGfPv8jQ7mVreDGNHPAf4E\nvBdAXVHk2xeWT3Lzqa0ZcA0wWVW7A2/bx04oUKqq+6pqz8BaR/59/kYHC9zQq+p8Vf3G73oMBg8Z\nCNgb8TIOOD5J2RztimswuMf46A2GeLaLSGOwHNguQTkF3hKRGSJyfjBNMxjSx5PwShGZDLRzuHSd\nqk6yy0wBhseuHIyQYWLQDL4SGYKWRGdHAONUtVVE2d9UtXVsQRFpr6pLRWRbYDIwVFWnxZQxem3w\nnVThlZ5kr1TVIz2QYR6BDYGRTGdFZLmItFPVZSLSHmufYycZS+2/K0TkBawsrdNiyhi9NuScoF03\nRukN+cDLwJn2+zOBF2MLiEhjEWlmv28CHIUVeGAwhA7fV8aKyJ+Ae4A2wBpgpqqGYL8+g8EZEWkN\nPANsDywETlLV1SLSAWszkj+KyI7A8/Yt9YH/qeotOWmwwZCC0KRAMBgMBoM/hCbqRkT+KSJfiMgs\nEXlbRDzZaklEbheRebbs50WkhUdyPVsIJiL9RWS+iHwrIld70T5b7ljb3+yZS0FEOovIFPt//1JE\nLvVIbkMRmW5//3NFxLPRsYjUsxc1TfJKZpr1541ue73A0Q/d9kOvbbme67afem3Ld6fbqfYaDOoF\nNIt4PxRrA3Ev5B4JFNnvbwVu9UjurkB3YArQIws5KffTzUJ2X2BfYI6H31M7YB/7fVOs9NNetbex\n/bc+8DHQxyO5w4D/AS979TmkWX/e6LZXem3L8kW3/dBrW64vuu2XXtsyXel2aEb0qrou4rApsNIj\nuZNVtXpjvOmAJ7uQqncLwWr201XVLUD1frpZo1ao3yovZEXIXKaqs+z367H2/e3gkewK+20JlpH4\nLVuZItIJGACMIUfBAPmk2x7qNfik237otS3XF932Q68hPd0OjaEHEJFRIrIIK9LhVh+qOAd41Qe5\n2dARiNwwdol9LvSISFeskdV0j+QVicgsrEVKU1R1rgdi7wKuBHK6C67RbaCO6rZPeg1p6HZght72\nJa21X3NiXscBqOoIVd0eeBzrn3Are7KDzBq5dpkRwGZVTbAtQGZyPSAvZ8NFpCnwHHCZPfrJGlWt\nUtV9sEamh4hIaTbyRORY4BdVnYmPo/l80+2A9BqMbgPe6zWkr9ueLJhyyWXAJCx/5cAUZZ8kjdGJ\npliwJSJnYT3iHOFWphu5HvETEDk51xlr5BNaRKQYmAg8oapxMebZoqprROQVYH9gahaiegMDxdrE\nuyHQXETGq+oZHjQzkrzS7YD0GoxuR+GhXkOauh3IiN6NL0lEukUcDgJmelR3f6zHm0GquskLmU7V\nZHHvDKCbiHQVkRLgL1gLdkKJiAjwKDBXVUd7KLeN2OmARaQR1kRjVjqgqtepamdV3QE4GXjHayNf\n4Lqd7VNQnddtP/Qa0tftoFw3bnxJt9iPj7OAUmC4R3XfizUBNtkOQ3rAC6Ei8icRWQz0Al4Rkdcy\nkaM+7qcrIhOAD4HuIrJYRM72QOzBwGnAYfbnOdM2ONnSHnjH/v6nA5NU9W0P5EbihyuhoHTbK70G\n/3TbJ70Gf3Q7CL2GFLodxMrYY4FjVPUS2zc1XFXjfIFikj8ZfEY9zjvjRreNXhuCIJVuBzGir/Yl\n/QBMAA4XkfFOBb2KLU30KisrK4g6zP+S/iuXup1Pn1k+yc2ntvop1w1BbDziu5/UYMgFRrcN+UIu\n4ujNo6yhUDG6bQglQYZXoqrvAu8GWWckpaWlBVFHUPUU0v/iN0Hrtl+fWT7Jzae2+inXDaHJXiki\nGpa2GAoPEUFzsAmI0WuD37jR7VClQDAYDAaD9xhDbwg1IrUvg6EQadcOKipSl8uGQH30BkM6xBp3\nETBeEEOhEKnfTZr4q9tBpUDwNfm+obD44APnEXwYjbzRbUO65OIJNZARvapuEpHDVLVCROoD74tI\nH1V9P4j6DflDog4QRiMPRrcN7tmyBUpKnK/5rd+B+ejVp+T7hsKgbdv8M/LVGN02pELE2cj/5S/B\n6HeQ+ej9Sr5vyHNEYMWK+POq4TfyYHTbkJhhw5IPYJ56Kph2BDmi9zz5viG/SeSr/Ne/8sPAV2N0\n2+CECNzlsMVMLgYwgUfdaJLk+yNHjqx5X1paWhArJg3xDBvm3AHAuw4wdepUpk6d6o0wlyTSbaPX\ndYsDDoAZM+LPH388vPBC9vIz0e1AVsaKSBtgq6qutpPvvwGUa0ReZrOCsG6QKz+8XytjU+m20eu6\nRS70O0wrY4NKvm8IKYncNCedlF9uGgeMbhvo3t1Zv7/8Mhz6bXLdGHzlvvtg6FDna0F+3SbXjcEv\nch0t5ka3zcpYg2/kugMYDH7SrRssWBB/ftMmaNAg+PYkw+S6MXhOIjfNKacYI28oDEScjbxq+Iw8\nGENv8JChQ5OP4p98Mtj2GAxek8gXH/Y1H8Z1Y8iajRuhcWPna2FWfoMhHZwMfHExbN4cfFvSJaik\nZp1FZIqIfCUiX4rIpUHUa/AfEWcjf//9dcPIG90ufBK5IlXzw8hDcHH07YB2qjpLRJoCnwHHq+q8\niDImOiHh9cNVAAAgAElEQVSPaNoUNmyIP9+yJaxaFXx7UuFjHH1S3TZ6nd84GfjGjZ11P1eEJo5e\nVZep6iz7/XpgHtAhiLoN3nLRRZbyOym6ajiNvJ8Y3S5MWrZMPIoPk5F3S+CTsSLSFdgXa3GJIU/4\n6CNL8R96KP7a8uV1w02TCqPbPnPYYTBtmq9V/Pabpedr1kSf79Ejv3U80MlY+9H2OeAye/RjCDlb\nt1oTTk6MGwdnnBFse8KK0W0feOUVOPbY6HOHHFL7fvVqaNHCs+qKiy19jyWfDXw1gRl6ESkGJgJP\nqOqLTmXqZPKnRPGIK1ZAmzbBtiWGRE074gh4661g25IuQSY1S6XbdVKvs8XNFkwtW1p/s7TEn3wC\nBx4Yf37YMLjzzqxE+0KYk5oJMA74VVWvSFCm7kxarVsHzZunLpejzyPRyKZVK+vRNh/xcTI2qW7X\nKb32gmbNYH0GD0QZfsaFsHrbjW4HZej7AO8Bs4HqCq9V1dcjyuR1h1ixbgU73bUT63Sd4/VuLbrx\n3EnPsVfHvd0LDfjz2HFH+OGHUDTFc3w09El1O9/12gt6P9ybj5Z9lLTMsJ7DuHPAv50vOn1+WW4q\n/NJLVtrgWF59FY45xrWYUBAaQ++GfOwQQycO5b4v73N/Q8S/13wtrKnOyR77f7dqBePHw3HHZd1G\nN5x8Mjz9tPO1PPtKEpJPSc1Wrl3JlZOv5LETH/OpVf6z7ahtWbl1pfsbtPbvB/dC71WkVr5YY9+k\niaungaKieNH16jk/xeYDxtD7iJR7ZzO0LDf/95NPwl//6nwtjImZsiGfDH2kbuVKNzJl8eLFbD92\n+/RuUkCij1tIC1aXrU59b6yxT+JYnz4devWKP79oEXTu7Lq1ocMYep9IZOQ7N+jMomsWxV+oqOCh\nA5pw5TGwvgmWUseICLJDr1hhbcbtxPz5sMsugTUlMPLF0MfqVr4Z+kR947sLv2PH7XZ0vPZjI6Hr\nFVhbq0NN33D9v8ca+2++sVJLRtCmDfz6a3SxFi2swJ18xxh6n0i7M8YqomqcjL4d+vLe+e950byk\nOD22guUpOv1036vPGcbQ+09s2wd0HMAr573i4sba++QGoJ71Pq3/PbaPrVtnLd92uASwbBlst517\n8WEmNCtjC5mMOmJlJVqm7N2qdmJ22s/+LgTp29dS+Fibc/rp1rlCNvL5xMKzFkYdS7l46iYMCi1T\nd0Y+ArmBWoukWAo7fLjLCmMUu1kz/u/hH+OMfKNGVtFCMfJuCSqp2VgRWS4ic4KoL0galTdKXShW\nCevXBxG+WPmFP42K4Oefrf7y/vvR57t2tZo1frzvTShovNbtLl26ONdjG/xnZz7rRTWhYfbPs5GR\nIGVYI3nBMvKb7AL//ndtVrHqV6NGMHp0vLCIfnYij3PchdFzBU89BRUV/vwfYSeoEf1jQP+A6vKd\nhjSseb+JTZz43xNT3xQx67miGOQfRI1edCSWEpeUwH77wYQJWbdzzz2hY8f486qJwygNaeO5bmuZ\ncvUeVzteO+nlk2qMvpQLW7Zs8bJqz3j2i8Q/SJe9ellN+/d+xH6qrR55K2yzEvRfSYRv2gRXXBH/\nA2BvQLwf03ieM6j91ahEFf7yF2/+t3wkMB+9nQdkkqrumeB63vjoId4fuW2Dbfnlml+S3vPTmp/o\ndEen2pELgMIh8+DdZzJoRLNm0K+fNbrZPnr04uSX/Ogj56iDuoCfPvpkup2tXve9vy/vr3w/dcEI\nTu5+MhNOyX6gkC6xfWLpFUs56dmTmLbEvVsyyhW6di0ceijMmuX6/qOZyJv8iWojX5/1bMFenJhH\n9iUdQjUZW2iGHpwjDGJ99lVVVezxwB7M+3VeXFkUnhkHf17ocbuoJPJhrX17y4VTl8lXQx8ly0Nf\nfRFFdGvejRN2O4FLe11Ku5btUt6jqrz//fs8Nfsp3l/8Pt+v/p71HqT1uf3w2/l737+7v+Gzz6w1\nJkuXRp3+meZ0ZDXVRr4VS/mNmEfaf/4Trr/efV2PPw5nnw19+vieUC1TjKEPAKfOd9ZuZzFu3jiU\nxP/PXm334ouLEvjoX3oJbr3VGsls2uRcJlF72EpN2ALKeE7k9E6fwuLFackpNArB0MfJzsNJWoDe\nHXrz7jnvUr+eh6m2Nm5EGjfAGuAosBmNcLHGUVlphaAl4M9/hueeA6ii9vHb+rvjjtbTcaIQ5aBx\no9uh2kowH5M/PXficwyeODjq3OPzHk9Y/s4j72RY72HJhQ4aZL1SsXUr3HgjPPIILFtGGUOJdPz/\nTHPasx6W4BxyU8AEmdQsFX7ptVPE18XPX8yDcx70RL6XvHnqmxzZ7Uj/KmgUGRQh6LDrIEFGBcBa\nCgtxGTDfe8/yFtUS/2Pw/ffRUTs9e1qLsYIitEnNIL9H9I98+gi3TLuFH9ZlNoN53j7n8cigRzxu\nVTyRfvk9eZ/Z9I0vFNLP2G8KcUSfKVu3buXGt29kwpwJLNqwiM242w+viCKa1WtGlxZdOKLrEVx+\n8OVs39p5FeyEzydw6qRTo861KWnDimtXZN1+J379NTrZa81H3rMnfPppagFt28LSpUi9WMMe+90l\nV6G774ZLA95MMjSuGxGZABwKbAP8AvxDVR+LKZN2hzjz+TMZPyd5fGARRdSzXRlVVFFJZVp1eEEx\nxWwu839zSYd1WVb+7kjfYsiMTlD4mNQsqW6H0dAHiZN7qfKGSoqSuE0yYfPm6JQdUR95ZMeoqEi4\nk72wAai+VsXZ3MxYbogTeuONUFaWvD377AMzZ7puflaExtC7IducILmia/OuPH/S8+zbcd+ac2s3\nrqXFbfEbIvznmP9wfs/zfWvLd9/BzjtHn3uC4/gr/1d7IiTfd9Dky8rYQsSpn+7Xbj9mXDDD23oi\nqrn/frj4YocLkd/FTjtZfpjqYjVBDIo6RZ47fI+VlZbXKFGU66BB8KLj7hveUdCG/v4p9zPkvSE+\ntqiW3Vrvxs1H3Mzxf3DIa5qEIx4/gnd+fCfuvJ/L2uNTHFgHHfiWxV8XUdR9Z8f7Cp26bOi3bt1K\n8agE24Q5cHjHw3n6pKdp09y7jW8e+PgBLnnjkrjzXvaF2Cfac49fxpgX28dUmKC+nj2RTz/CCmSI\nMfSXXea8QCuGESPg5pudr/mpAgVt6L0YzRdTzPG7HM/jgx6ncSPnx7lsUVWKbowfHfhm7P/4R258\ndSfKGI1j9jSbCy+EB8M3Z+cbddXQ+/nU26a4DcMPHs41h16TVXua0Yy1ZWs9aVO0sVdgK1WUWL0g\nxffQoUNtxGavXlZkTSaMGwdnnRV/vnFjfzYWN4beI+pRj51a7MSpe53K8IOH07RB07RlSLnUpmNV\nOPI9eHMK1obH78SP+tPiscfgnHPiTt/AVdzELSQz+NUUFVlpRW67LbumhBVj6IOniCK6terGefuc\nx5CDhtCw2Ap3vH7y9Yz6cFRtQTvlgV71K7RunXmFn38O++2H8DtQTOSqxKKiSiorUwcZVv9Q3Hsv\nDMnSYbDbblY22Fi8Vgdj6MNEbDjuZtDYxzy3///rr7vbBmfUKLjuOgDOPx/GjHHZ1hgOOgjeeMNa\niJuvGEMfzdgBYzn7gLOjzi1btYwznz+Td5a8w1YC2oUjYvCj5fa5tm1h+XJ39/fvbylnDBdyBw8z\nDKcBTrIV4pWV1t4lHu45nu1mWC7kh8jQi0h/YDSWE2yManQ2i0w7RLoG/9tLvmXnNpaf+rFPHqP8\nvXJ+3PBj2vVmhcLpk2D85wmuN2gAv/+eufzZs61ENwnYsgV694YZHs2FFRdbez0MHeqNPD/wMerG\nF732kkR9JFv34cMfPcyoD0axeIMHi/EiDb1XqCbdgtblhlSe4KexD42hF5F6wNdAP+An4FPgFFWd\nF1Em6w6RrtGvT30WXraQji0dMn/FMOW7KYz+eDQzfp7BiooVbCGDZFJqvdrMhxWZ5LZJxsMPw9/+\nlvHtS5dakZgLFnjYpiwJetSTgcxA9NorcpHrfuFvC/nXtH8xcf5EVmxyiKG3+wTfgz7hQYVHHAFv\nvRV3ulWr5JuMBLHxfevWsGpV9DkvVCNMhv4goExV+9vH1wCo6q0RZbLuEFO/m8phTxyW8f2tS1qz\n4NIFtGrSKqt2pGS77eCX5AnQUpJuzo4M2bgRjjwSPvjA96riWL/eGnV5gU+GPhC99pqHP3mYfjv2\nY6c2O+W6KbU0bmwpWybMnWs5xF0wdCjc52Kb59des7xCXuO41iVrmeEx9IOBo1X1fPv4NOBAVR0a\nUcbzDpGtH3/fbfflkws+8TYnRzL8dub5wObNVnbYl1/2vql5MKLPiV7XCQLoC8ncOpF4nRQw8l/r\n0AF++ilbeeHJdZMTTY98NH3/h/fpO94hJUASZq6YSfFNtfHHA3YYwCtnpLdrTlqowu67WyOUakKe\no6akxP8FISEmvF9MvqMab+x79YKPP/asinXrat/Xr29NxDqxdGl0U0pKsptCiySorLJBGfqfgMh9\n1jtjpdqKws+kZn126BNl+D9d8im9Hu1FFVWuZbz6w6tRTwn3HnYvQw7xeNHWV185j2YMaRFQUrOc\n63VBE2vsfcwctjUiyChV99u82bnMv/9t7YeSjNj7und3175IQpvUTETqY01aHQH8DHxCyCatNm7e\nSI8HezB/tUPgayLsiSS90T5OMBGUEX448+owPrluQq/XBUEO+8Ly5dAudar+jAnKRx/IVoKquhUY\nArwBzAWejuwMYaBRSSPmXTYPLdOa16PHPkr9ZA89seuQ3n6buO3NXsnQ1XPnnZndZwiMfNDrgsCr\nGfkM2G47yxhXv7xy2UCwY7e8XTCVK3a5axe+WftN7QKoStCb0hSydKm7YUKiZEyGtKmrC6YKhpD3\nhR13dL8P85o10Ly5d3WHJurGDXnbITZtitn0IEMid7zJw+ibsGMMfZ5QrfvVyWauv95a4R2J+Tyj\nMIY+F3z9Ney6q7cyC+FzyTHG0OcBboMQzOcZRWh89HWKXXaJduqpZpck5l//Sl3GYCgE3CQ0M0Y+\nI4yhD4K1a+ONf6LE1dXMmWOVu+qqYNpoMOSaX39NbMir+40hI4zrxlAnMK4bQ6ESCteNiPxZRL4S\nkUoR6eF3fckIYAFNIHUEVU8h/S9+kEvd9uszyye5+dRWP+W6IQjXzRzgT8B7AdSVFGMcw1dHkPX4\nQM50O9+MkTH0udVz31MgqOp8sB4vDIZCwui2IV8wk7EGg8FQ4HgyGSsikwGnpZ7Xqeoku8wUYLiq\nOu6rJCJmxsrgK5lMxmar20avDUEQSJpiVT3SAxnm+dcQOrLVbaPXhjAQtOvGKL2hUDG6bQgtQYRX\n/klEFgO9gFdE5DW/6zQYgsDotiFfCM2CKYPBYDD4Q2iibkTkdhGZJyJfiMjzItLCp3p8W+QiIv1F\nZL6IfCsiV3spO6KOsSKyXETm+CHfrqOziEyxP6cvReRSH+poKCLTRWSWiMwVkVu8riOirnoiMlNE\nJvlVR4r6/2nr9SwReVtEOqe+y5Vcz/uM1/3Djz7hVx/wQ+/91nPXuq2qoXgBRwJF9vtbgVt9qmdX\noDswBejhodx6wAKgK1AMzAJ286H9fYF9gTk+fhftgH3s902xdlHy439pbP+tD3wM9PHp/xkG/A94\n2a/PLEX9zSLeDwXGeCTX8z7jZf/wq0/41Qf80ns/9dytbodmRK+qk1W1egPX6UAnn+qZr6rf+CC6\nJ7BAVReq6hbgKWCQ15Wo6jRglddyY+pYpqqz7PfrgXlABx/qqbDflmAZhd+8rkNEOgEDgDHkaMJU\nVSO2oaYpsNIjuZ73GY/7hy99wq8+4Jfe+6Xn6eh2aAx9DOcAr+a6EWnSEVgccbzEPpfXiEhXrNGT\n5zszi0iRiMwClgNTVHWu13UAdwFXQhq7wPuAiIwSkUXAmVijb68JY5/J2z7hpd77qOeudTsrQ+/G\npyUipSKyxvYjrRWRZSIyJ+Z1XET5EcBmVX0yi3ZNdqgjqh4fKLhZbRFpCjwHXGaPcDxFVatUdR+s\nkeghIlLqpXwRORb4RVVnksVo3sknLCKtbT37RkTetPtBQp1T1RGquj3wOFYHdVt3Sl1Ot88E2D/y\nsk94rfd+6Hm6up3tgqktwBWqOsv+cD4Tkckav0Hyu6o6MJUwETkL61HkiGwapR4s4MqAn4DISbbO\nWCOYvEREioGJwBOq+qKfdanqGhF5BdgfmOqh6N7AQBEZADQEmovIeFU9I005jwH3AuMjzl0DTFbV\n2+xJxlaqeo0LWU+Sxsg7lS5n0mcC7B951yf81HuP9Twt3c5qRJ+GTyvlL46I9Md6DBmkqpuyaVca\neOmznQF0E5GuIlIC/AV42UP5gSEiAjwKzFXV0T7V0UZEWtrvG2FNLM70sg5VvU5VO6vqDsDJwDsZ\nGPlEPuGBwDj7/Tjg+ET3i0i3iMNBePR/BtBnsu0fedUn/NB7v/Q8Xd32zEefxKelQG87BOxVEflD\nAhH3Yk1UTbbdPA941baYdvqyyEVVtwJDgDeAucDTDk82WSMiE4APge4islhEzva6DuBg4DTgMPu7\nmGkbFS9pD7xj+y6nA5NU9W2P64jFS1fCdqq63H6/HNguSdlbbNfILKAUGO5RGzzvM172D7/6hI99\nwA+9D0rPk+q2V0nNmmI9itwU+7gjIs2ASlWtEJFjgLtVtbuDjLz05xnyB80i74w9kJmkqnvax6tU\ntVXE9d9UNW7TU6PXhiBIpdtZj+hT+bRUdV11eJGqvgYUi4jjLsBexZaWlZUZWUZW1MsHlotIO7sP\ntAd+SVQwjJ+HkVUYslTd6Xa2UTcpfVoisp1dDhHpifUU4Xm8tMEQMC9jhUpi//V1wtpgyIZso26q\nfVqzRaR6guE6YHsAVX0YGAxcJCJbgQqsiQODIW+wfcKHAm1s//U/sGLhnxGRc4GFwEm5a6HBkJys\nDL2qvk+KpwJVvR+4P5t60qW0tNTIMrI8Q1VPSXCpX5DtCOtna2TlTpZbspqMFSs503igLdas739U\n9R6HcvcAx2CN6M9SK8g/toxm0xaDIRkiguZgExCj1wa/caPbvi+YsgP6d1bVbiJyIPAgVuiWwWAw\nGAIgiAVTNQtLVHU60FJEksUcGwx5g4hcK1YKkDki8qSINMh1mwyGWIJYMOWU2MiXzJQGQ5DYOn8+\nVjrfPbEyE5pgA0Po8GRzcBdJgGL9R6F0Wj71+VOMmzmOr1d9zfrf17NFt1BcVEzzhs3p1KITh3c5\nnIt7Xkyb5m1y3VRDOFiL5b5sLCKVQGOs/C4GQ6jI2tC7SAIUm9ioEwk6w8iRI2vel5aWejY7/cC0\nB7jknUsyu7kSVmxZwXfrvuPdJe9S9kFZ0uIHtj+QqWdPpWFxw8zqM3jC1KlTmTp1qq91qOpvInIn\nsAjYCLyhqm/5Wqkh75Dy2nGuluVmjJtt1I1g+d9/VdUrEpQZAAxR1QEi0gsYrapxk7FeRSds3ryZ\nBreEx03asUlHvrvsOxoUh6dNdRE/om5EZCdgEtaOR2uAZ4HnVPV/EWW0rKx2cODlAMYQbk5/+nSe\nmP9E1DkvDH3sIKa8vDylbmdr6PsA7wGzqXXHxC6YQkTuA/oDG4CzVfVzB1kZG/o1a9bQcnTLjO6N\npYQS2jdsT8fmHdmm0TaW/M1rWL5hOcvXL2dN1Ro0Q89Tn059mHbutMQFnngCxo2DevXg4othoJ3Z\nWSK+QxOqlxE+Gfq/AEeq6nn28elAL1W9JKKMCa+sg0SO4iPxY0TvRrc9SWrmBZl0iDFTx3D+u+e7\nKtuuXjuWXr807vyXP3/Jno/smVa9kTz6x0c5Z/9zao7/+uxfeXJu8v0fTt71ZCb0uQM6ZTAnvWkT\nNDBPB+nik6HfG2u/zgOATVibinxiLxKsLmMMfR3iwPsO5JNfP4k7//LAlzluX3/2PSpoQ5/oF7Oa\nZL+cJ084mae/edp1XW654oAr+PeAf0edu/j/LubBzx6MaRywBfTmLCusqoLVq6G1Q464hx+Gv/0t\nywoKB78WTInIVVi5bqqAz4Hz1Noftfq6MfR1hCBH8VH1BmHoRWQs8Eesba3ihsZibZv1EvC9fWqi\nqt7kUM51h8j0A0314+AVe26zJ7OHzI46t2LtCtre1dYy8oL1twL09pibmzWzDPiGDd426tlnYfBg\nb2XmEWZlrMEvOt7ckZ+3/Bx3fvXlq2nRooXv9Qdl6PsC64HxSQz9ME2xlaDbDjHyjZGUf1wedS6V\ngW9Q3oDNbE54fY+WezDnsjkJryfigx8+oM/4PgmvjzluDOf2OLf2hAjyd6AJNcZeRyZou/hkk7p3\nh6+/9kd2iDGG3uAHuRrFR7UhKNdN7KYMMddKgeGqmtRB5bZDDHh8AK/9WLvpTaaj+HP3OJcxJ45J\nWZ9bvv/1e3a6b6e4883qNWPt9WvtxlhtkeupCWw9cecTee6vz0Xf5JeRj2Tx4szmCPIUY+gNXtK4\nvDEb2Rh3Phfhk25027OVsUlwu5WgK14961Wa0IQGNMjIyC+4cAFapp4aeYAdt9kRLVPWXrU26vy6\nynVWO0pKas5t+0Pt9YlzJlqGXewyboz8qada0TeqsMsuzmUuvLC2jBOdOwfzg1LgiEhLEXlOROaJ\nyFw7hNhQwEi5xBn5EkpyFiPvhiBG9K63EvQy3jjWyDegAZvKgtpzHPZ7cD8+X/55jYtGtkCVPfkq\nI4Bi633Pr2D6sy4Ennoq/O9/ia83awbrYxYlFxVBZWXtcdeu8OOP8fcW4Igzk1jjTBCRccC7qjpW\nROoDTVR1TcR1M6IvEMLgpnEiFK4bh7I/APvF7jLlZYeI/UK+PudruneO+23xjyuvhDvuoOxAuLE/\nNcZ+8Ofw7CSQstpzWp5cVA2NGkFFRepyTqP02M/VTZkCw6fwyhbATFXdMUkZY+jznK1bt1I8qjju\nfJv6bVgxYkUOWhRNEGmK3TRiO6yIHA1iK8HHpz8edfzJaZ8Ea+R79oRPPwWgfDq81xWm7goIPLcf\nyH4RZavSkLtxY62B7tkTpsfmjrNRhS5dYNGi2nMi0YZc1VqUVVWVuIzBDTsAK0TkMWBv4DOsfE8u\nfpEN+UBYR/Hp4kXUTc02a8ByoAzbMaGqD4vIJcBFQPVWgsNU9WMHOZ6MfGK/mEC/kO+/h51iJmQ3\nb0ZuLnEsXvWPKqTaeMeOshcvtvzoqXjuOTjxxPjz990HQ4fWHB7O00zhJBo1gv33hzFjoPux3eHb\nb2vvadIk3v1TIPg0ot8f+AjoraqfishoYK2q/iOijEmBkIesXLWSbe/ZNu58t6bd+Gb4NzloUS2B\np0DwEj8MfeC/urHG2v5/Vm1YRes7ohc1Hd/teF449YWk99WwZUvUZG5CZs6EffapPR4zBjn/TKwH\nt0R6oPZrKyMYyU0//M3y5RcYPhn6dsBHqrqDfdwHuEZVj40oY1w3eUa+jeIDiboRkbEislxEEgai\ni8g9IvKtHXmzb7Z15gURnTvWyHdp3iWxkd9tt3hZxcW1ETQz43ZhrGXffWsjeIqLGTt2JcmNPPa1\nIqCEUYxCduhSI+KZZ5LcZkBVlwGLRaTaN9gP+CqHTTJkwbyf5zka+f7b9w+tkXdLEAumIrNXHogV\ndeNb9srQjOhV+f3332l4a3y64pp2ZTspOnSo5aJJ1BzWAU3toy1cygj+NOU2Hn4YXn/dyp4QUTHJ\nfhD69oX33nPftLDhYwqEvYExQAnwHVbSPhN1k2fk2yg+kkBG9Ko6DViVpEjOthIMKuWBE0c9fpSj\nka/Bi8iXe++tHelfdlnc5VGcT21S0XrcrbdRWgoTJsCqVbW3qsIPp15HfVZhzRDHt2PatNqHhUMO\nSa+ZhYyqfqGqB6jq3qp6QqSRN4SfKd9McbQTf9vrb3lh5N0SxIKpQLcSfObYaH/Dzjfv7FdVCZEb\nYPLCyc4XFX/CG0ePjrbcd93FdRKZuK2IlkkyOXf93y1soTVKPZQiRo9OXDbS6L/wQnbNNhhyhZQL\nh084PO68likP/+nhHLTIP4JYMDUJuFVVP7CP3wKuis1J72V0gtMvdBC/zseM2pPXt3xZ6wFRoBL0\nphSx8z4+2g8cCJMm1R7ff7+V6t6R886DRx+1LHhE6KVTMJETYfJQBLVgKhXGdRM+Ej3p33LwLVzT\n75qAW5M9oVgwJSIPAVNV9Sn7eD5wqKoujynnaYcI2tjH1adwzGx49QX4pinsMhzL0FeC/tMu889/\nwvXX+9ammrbFNO3LL2H33TOT9e67kOr3d80aaN48M/l+4WeuGxGpB8wAlsTmdDKGPlzksy8+EWHJ\ndfMycIbdoF7A6lgj7wdBfXHtbm3n/KNSbhl5gF2uoGY0f8NbWJE1qoEYeYgfae+xByxcmJmsQw+t\n9Q7tsYdzmRYtrB+Xl17KrI485DJgLiHd9N4A4z4Z59hPR/UdlddG3i2+L5iyy/i6lWAyisqLUJRO\nTTuxePji1DekgZPi7LXNXnwx5Iua49vev42r37665jiXShU7so8Nu8+UigprrVUiJk+Gfv2yrycb\nfIy66YS1s9QorMWAZkQfMgpxFB9JQe8wlUualjdlA/EbgzgpTqSSXdjjQh487sG4MkESa+wnTICT\nT/ZPfiS5/Hp9NPTPAjcDzYG/G0MfHv7z4X+4YPIFcecnDZ7Esbsf63BHfhJIrhsR6Q+MBuoBY1T1\nXzHXS3Gxw1S+4DQ66NmmJ9Mvic89E1s210YeLGMbaYxPOQVmzIA77vBOPjgb/OpzhWL3RORYrDxO\nM209d2TkyJE1700KhGAo5FF8bKCBG7Ia0duTUF9jrQj8CfgUOEVV50WUKcXDHaZyidsJ3oqNFTS5\nrUnKcrkk1hAny5PmZT3V/PILbBufSsQ3fEqBcDNwOlYep4ZYo/qJqnpGRJnQ63Uh8fgnj3P2a2fH\nnX/7lLc5vHt8KGUhEMSIviewQFUX2hU+BQwC5sWUy/sdLtwa+aMeP4rJP0bH0G++LvE2hrkidmT/\nyV+NVOQAAAz3SURBVCdWpMzatYnvybQeiDf4bdtGX89HVPU64DoAETkUy3VzRvK7DH5RyKP4bMk2\n6sZpMVTHmDKe7jCVC2IVaPQho+OUZ23FWqRc4oz81uu3Ulwcn8s6DMQmqly3LrmPPRtUnZNxFtgm\nV8ai5ICX5rzkaOSnnDrFGHmbbEf0bj7Fz4HOETtMvQgEmCDeW7478zt27Fq7z8TjMx7n7FfiHxUh\n/COJL790Pu9XavrqFPmxxr0QUuGr6rvAu7luR13DjOLdka2h/wmIHKd1xhrV16Cq6yLevyYiD4hI\na6fNR/Jh0mqncamXiB7c6WDeP/f9AFqTHQcemPian8Y31m3kR32ZTFgZ8ocPF3zIwf87OO78Gye/\nwVG7HJWDFoWbbCdj62NNxh4B/Ax8QvxkbOwOU8+oalcHWaGdtEonOVq+jSQiDe7220dvTAX+jrRT\npeH3ti7/VsamqDe0ep2vmFF8NL6vjFXVrcAQ4A2slYFPq+o8EblARKoDWAcDc0RkFlYYpodR2+Hg\n4HYHo2Wal4p211217xctgssvj77upw/dzVa2YUdEOovIFBH5SkS+FJFLc92mQqWiosLRyD953JN5\n2feCxCyYSpOKigoaN26c62Z4SqSBvfNOK4HZ8cdHl6mshCIfEmbsuWf0XIFfKuDjgql2QDtVnSUi\nTbH2jT2++qk2X/Q67LS7pR3LN8dnTjEG3qyMNbjk5JPh6YiMxqqwbBm0bx9dbv365KkOMiEo901Q\nrhsReRG4V1Xfto+NXmeJ0yi+S70uLLx+YfCNCSFBbSXYX0Tm21sFXp2gTN3bSjCPeOqp6OOSEmjX\nLt7oNm3q3d7hH3yQn66aZNhZXPcFfFh6VqDstlvt5gbVL5v7P7o/4foVY+TTI6uoG3tl7H1ErIwV\nkZdjJmMHADurajd7K8EHgbitBA25Zdo0a7tAsPYiX7AAdt45PkKmWbP03DjpGPN8HvjabpvngMtU\nNernMB+iyXKCKsyfH39ehAZ/h81No0/3aNuDzy76LJi2hZhcpEA4CChT1f728TUAqnprRJmHgCmq\n+rR9HEg+ekP6JHOjZOJiqayE+i6HEn5/9T7noy8G/g94TVVHx1wzep0IVccRg9yA5WuwU3sjxhef\njCBcN25Wxga6laAhc5JFwWQSIVOvnrs689kOiogAjwJzY428IQUiVvrUjh1BlUk9mlk7sdWj1shv\nBh2Z01YWBEGsjIX4XDeO95lH3Nwzfz7sumvtceRCpsrKaOPtZpFTrox4gAumDgZOA2aLyEz73LWq\n+noQlecdCUYIMkLgOIjchvOi9+DBKRUIJSCVQD222QZWrgyqsYVDtq6bXsDICNfNtUBVZKriXG0l\naMicrl3hxx9rj9u2heX2tzVhApx6anT5fPjaCnHBVGVlJXdOu5OrSq/yRb6nOBj4FpfA2tbU+hWq\nR/EKeiMIK4HWJMuJmA+65ze+h1e6XBk7ABiiqgPsH4bRqho3GWsMfbiI7Zc77WRN0Dpd+8Mf4Kuv\ngmlXphSioU9nxXZTmtKvaz+uPuhqenXPLhaiYmMFHy/5mPcWvcec5XOYs3wO3679NvENVTjbaqfn\n/K2go6ovLwPaJrg54rY6bjYCiaO3E5VVbzzyqKreUr0qNgxbCRoyJ9agN29ubfztdC3sX52PC6ZS\nbbwTCkMfWhS4+wFYfSrQGGsX0sqIAhLxckfr1vDOO7D33l42NLyYBVOGrHFyqSbKMR/mr8+njUfc\nbLzjm14vXbOUDqM7+CI7axL9y9Xnq4CbVmPt1QLpGPJs+fOf4ZlnAqvOd3w19CLSGnga6AIsBE5S\n1dUO5RYCa7F+preoas8E8oyhDymJjP2MGXDAAdHnwopPht5NeHFO9Xr9+vWc99J5TF4wmd+ISxib\nESWU0LykOW0bt2XHljuyV7u96NOpD/3/0B9JI0+GsIVs4kHuuguuuCLj21PSqBGMGGG9wozfhv42\nYKWq3maviG2lqtc4lPsB2M8pLXFMOWPoQ4zbhU9h/Qp9MvSDgaNV9Xz7+DTgQFUdGlGm7uh1KiVp\n1crafWb27NpbqKR2NO+NL/6UU+JXe/tJrr9ev+PoBwLj7PfjgOOTlC0AZ2LdJtfKHFLMp1JNiiRI\nTRmHrFqBzJ6JoDWv2pVRqU2Em3UZYEWGVa/PiH3de68/yfnCTjZx9NtFhEguB7ZLUE6Bt0SkEnhY\nVR/Jok5DDlG1Im+6dct1S0JDyo13oIDXh6SR32IDp5PteO/QQ7O6HYAhQ6xXKm6/HW64AX7/PXm5\ndu2yb1O6eJ4CQUQmA07/yghgnKq2iij7m6q2dpDRXlWXisi2wGRgqKpOcyinZWVlNccF1SEKkGST\ntGEgtjOUl5f74bpxE15cmK6bNDPSScIYS3fsuWeUx8cQgd8++vlAqaouE5H2WPlsdk1xTxmwXlXv\ndLhWmB3CEAp8DK+MCy+OuV54ej1pEgwcmLSIsAEoId4t485Ns2WL+zxJdR2/ffQvA2fa78/E2vQ7\ntgGNRaSZ/b4JcBQwJ4s6DYZQoaqvqeouqrpzrJEvKC6+mJo0wimN/FagEZZnuB6Wmal+ufutjd3p\nzJAd2YZXPgNsT0R4pYh0AB5R1T+KyI7A8/Yt9YH/JeoMBTnyMYSGQlwZGyhpuGqsSJrsZjwL4SML\nCrNgymCwMYY+C77+OjrTXTU9esCHHyING0BUmGQ17le05vtHlEt8dd2IyJ/tDZErRaRHknIpd6Ay\nGPIREbldRObZO6c9LyItct0mX9hlF+fzn38ODRtiBdbFumjcu2kM/pPN89Uc4E/Ae4kKROxA1R/4\nA3CKiOyWRZ2u8DI9rZFVGLJ84k1gd1XdG/gGuNbPynL62SYdck+lJu1kBgwdWvs+rPoTVlluydjQ\nq+p8Vf0mRbGewAJVXaiqW4CngEGZ1umWsH4pRlbuZPmBqk5W1Sr7cDo+b6iT8882cuXRpEk1p0t4\nHtiMHnhQwoVKyV733JNluxJQF2S5xe81Ym52oDIYCoFzgFdz3YjAOPbYGkt9bdk2qDaAjz/OdasM\nCUgaqZpkwdR1qjrJ4XwsZorFkNe46QMiMgLYrKpPBto4g8ElXuSjnwIMT5BjPuUOVBFlzY+CwVd8\nWjB1FnA+cISqbnK4bvTa4DupdNurtWeJKpkBdBORrlhLxP8CnOJUMBehbwZDNtibjlyJtTVmnJEH\no9eGcJBNeOWfRGQx0At4RURes893EJFXAFR1KzAEeAOYCzwdmQfEYMhz7gWaApNFZKaIPJDrBhkM\nToRmwZTBYDAY/CE0mZlF5J/2wpNZIvK2iHROfVdCWZ4tZHG7MCyFDM8WjYnIWBFZLiJZ5wwSkc4i\nMsX+/74UkUszlNNQRKbb391cEck654uI1LNHyW4m/VPJWigis215n2QrL4P6jW6nlhM6vbZlhVa3\n09JrVQ3FC2gW8X4o1kbLmco6Eiiy398K3JqFrF2B7sAUoEcG99cDFgBdsXY+ngXslkV7+gL7AnM8\n+MzbAfvY75tipdzNqG1AY/tvfeBjoE+WbRsG/A942YP/8wegdbZysqjf6HZqWaHUa1tGKHU7Hb0O\nzYheVddFHDYFVmYhy7OFLOpuYVgyPF00plYu/1VZtCdS1jJVnWW/Xw/MAzLabVpVK+y3JVgGIOMN\nSkWkEzAAGIN36+hzNilqdNtVW0Kp17aMMOu2KxmhMfQAIjJKRBZhpT2+NVV5l+R6IUteLBqzI6P2\nxTIemdxfJCKzsHYbm6Kqc7Nozl1Y0SxVqQq6pHqXsxkicr5HMtPC6HZuyFavbRlh1W3Xeh1oav9U\ni09UdQQwQkSuwfpAzs5Ull3G1UIWDxaGJSP0s90i0hR4DrjMHgGljT3K3Mf2Gb8hIqWqOjWDthwL\n/KKqM0WkNJO2OHCwRuxyJiLz1WGXs2wwuh0+vNBrCLVuu9brQA29qh7psuiTpBippJJlL2QZgLXN\nm1ftygRX+4rmChEpBiYCT6hq3OYx6aKqa+zw2v2xsl2lS29goIgMABoCzUVkvKqekUWbltp/V4jI\nC1guB08NvdFtIES67bVeQ/h0Ox29Do3rRkQit5weBMzMQlb1QpZBmmAhS6aiM7inZtGYiJRgLRp7\n2cM2ZYyICPAoMFdVR2chp42ItLTfN8KaMMzo+1PV61S1s6ruAJwMvJONkZcQ7HJmdDtYvNJrW1Yo\ndTttvc5m1tfLF9Yj1hysmfuJQNssZH0L/Ij1hcwEHshC1p+w/JAbgWXAaxnIOAZr5n8BcG2Wn9ME\nrFXGv9vtOjsLWX2wfIWzIj6r/hnI2RP43JYzG7jSI504lOwjE3aw2zUL+DLbzz/DNhjdTi0ndHpt\nywqlbqer12bBlMFgMBQ4oXHdGAwGg8EfjKE3GAyGAscYeoPBYChwjKE3GAyGAscYeoPBYChwjKE3\nGAyGAscYeoPBYChwjKE3GAyGAuf/AX9xwERO4Jk0AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "serial_time = [time]\n", "serial_x1 = [0]\n", "serial_y1 = [0]\n", "serial_x2 = [0]\n", "serial_y2 = [4]\n", "serial_x3 = [3]\n", "serial_y3 = [0]\n", "serial_energy = [E]\n", "\n", "dt = 0.1\n", "for i in range(6):\n", " for j in range(5000):\n", " minDist, maxVelo = evolve(dt)\n", " dt = min(minDist / maxVelo / 100, 0.1)\n", "\n", " x1 = x[0]\n", " y1 = x[1]\n", " x2 = x[2]\n", " y2 = x[3]\n", " x3 = x[4]\n", " y3 = x[5]\n", "\n", " serial_time.append(time)\n", " serial_x1.append(x1)\n", " serial_y1.append(y1)\n", " serial_x2.append(x2)\n", " serial_y2.append(y2)\n", " serial_x3.append(x3)\n", " serial_y3.append(y3)\n", " serial_energy.append(E)\n", "\n", " plt.subplot(321+i)\n", " plt.scatter(serial_x1, serial_y1, facecolor='r', edgecolor='none', marker='.')\n", " plt.scatter(serial_x2, serial_y2, facecolor='g', edgecolor='none', marker='.')\n", " plt.scatter(serial_x3, serial_y3, facecolor='b', edgecolor='none', marker='.')\n", "\n", " serial_time = [serial_time[-1]]\n", " serial_x1 = [serial_x1[-1]]\n", " serial_y1 = [serial_y1[-1]]\n", " serial_x2 = [serial_x2[-1]]\n", " serial_y2 = [serial_y2[-1]]\n", " serial_x3 = [serial_x3[-1]]\n", " serial_y3 = [serial_y3[-1]]\n", " serial_energy = [serial_energy[-1]]\n", " \n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 尝试动画的效果" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import base64\n", "\n", "from tempfile import NamedTemporaryFile\n", "\n", "VIDEO_TAG = \"\"\"\"\"\"\n", "\n", "def anim_to_html(anim):\n", " if not hasattr(anim, '_encoded_video'):\n", " with NamedTemporaryFile(suffix='.mp4') as f:\n", " anim.save(f.name, fps=200, extra_args=['-vcodec', 'libx264'])\n", " video = open(f.name, \"rb\").read()\n", " anim._encoded_video = base64.b64encode(video)\n", "\n", " return VIDEO_TAG.format(anim._encoded_video)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from IPython.display import HTML\n", "\n", "def display_animation(anim):\n", " plt.close(anim._fig)\n", " return HTML(anim_to_html(anim))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "# First set up the figure, the axis, and the plot element we want to animate\n", "fig = plt.figure()\n", "ax = plt.axes(xlim=(-3, 5), ylim=(-2, 8))\n", "\n", "liner, = ax.plot([], [], marker='.', color='r', lw=2)\n", "lineg, = ax.plot([], [], marker='.', color='g', lw=2)\n", "lineb, = ax.plot([], [], marker='.', color='b', lw=2)\n", "\n", "# initialization function: plot the background of each frame\n", "def init():\n", " liner.set_data([0], [0])\n", " lineg.set_data([0], [4])\n", " lineb.set_data([3], [0])\n", " return [liner, lineg, lineb]\n", "\n", "# animation function. This is called sequentially\n", "dt = 0.1\n", "def animate(i):\n", " global dt\n", " minDist, maxVelo = evolve(dt)\n", " dt = min(minDist / maxVelo / 100, 0.1)\n", "\n", " x1 = x[0]\n", " y1 = x[1]\n", " x2 = x[2]\n", " y2 = x[3]\n", " x3 = x[4]\n", " y3 = x[5]\n", "\n", " liner.set_data([x1], [y1])\n", " lineg.set_data([x2], [y2])\n", " lineb.set_data([x3], [y3])\n", "\n", " return [liner, lineg, lineb]\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%matplotlib inline\n", "\n", "from matplotlib import animation\n", "\n", "# call the animator. blit=True means only re-draw the parts that have changed.\n", "anim = animation.FuncAnimation(fig, animate, init_func=init,\n", " frames=30000, interval=5, blit=True)\n", "\n", "# call our new function to display the animation\n", "display_animation(anim)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.4.3" } }, "nbformat": 4, "nbformat_minor": 0 }