{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Assigning particles unique IDs and removing particles from the simulation\n", "\n", "For some applications, it is useful to keep track of which particle is which, and this can get jumbled up when particles are added or removed from the simulation. It can thefore be useful for particles to have unique IDs associated with them.\n", "\n", "Let's set up a simple simulation with 10 bodies, and give them IDs in the order we add the particles (if you don't set them explicitly, all particle IDs default to 0):" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]\n" ] } ], "source": [ "import rebound\n", "import numpy as np\n", "\n", "def setupSimulation(Nplanets):\n", " sim = rebound.Simulation()\n", " sim.integrator = \"ias15\" # IAS15 is the default integrator, so we don't need this line\n", " sim.add(m=1.,id=0)\n", " for i in range(1,Nbodies):\n", " sim.add(m=1e-5,x=i,vy=i**(-0.5),id=i)\n", " sim.move_to_com()\n", " return sim\n", "\n", "Nbodies=10\n", "sim = setupSimulation(Nbodies)\n", "print([sim.particles[i].id for i in range(sim.N)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's do a simple example where we do a short initial integration to isolate the particles that interest us for a longer simulation:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUIAAAE4CAYAAAAuFPo7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXe4LEW1xX/7Xi45pwuKJEEQBSVjQEoERVBEQTAHkGBO\n+CSIc0efARVzQjFjzqAimMoEBkDFjAEQVMAEDxRUYL0/qnumz9w5oXb3zEm9vu98Z2ZO1a46M9Or\nd1XtvbZJokWLFi0WM5bM9gRatGjRYrbREmGLFi0WPVoibNGixaJHS4QtWrRY9GiJsEWLFoseLRG2\naNFi0aM2EZrZe83sOjP7aeW1Dc3sK2Z2uZldYGbr1x2nRYsWLUaFJjzC9wEHDbx2EvAVSXcDvlY8\nb9GiRYs5CWsioNrMtgbOlbRz8fxXwH6SrjOzzYAoacfaA7Vo0aLFCDCqPcLlkq4rHl8HLB/ROC1a\ntGhRGyM/LFFyOds8vhYtWsxZrDIiu9eZ2WaSrjWzzYHrBxuYWUuOLVq0GAkkWU77URHhOcCTgdOL\n358b1ih3sosBZrZC0orZnsdcQ/u+DEf7vqwMj5PVRPjMR4ELgR3M7GozeyrwauBAM7sc2L943qJF\nixZzErU9QkmPneRPB9S13aJFixbjQJtZMvcQZ3sCcxRxticwRxFnewILAY3EEboGNlO7R9iiRYum\n4eGW1iNs0aLFokdLhC1atFj0aImwRYsWix4tEbZo0WLRoyXCFi1aLHq0RNiiRYtFj5YIW7Rosegx\nqlzjFi2Gw8yATYGdip97FL/vVry+NNPiv4E/A38CfgH8pPcj/V9Ds26xwNEGVLdoHmarAw8BDiXl\nmm+d0fsm4K/A34qff5FWLkuLn2XAWsCGwObAuhm2LwfOL36+jnRLRt8W8wQebmmJsEU9mG0JPAM4\nDthgklbXAF8vfiLwB0b9xTNbAtwV2IuU9/4QEnEOwxXAB4APIF050nm1GDlaImwxWpgtBY4CXgTc\ne0iLbwEfBT6N9JdxTi0bZpsAjwWeBOw+pMVFwGuBzyPdMc6ptaiHlghbNA+z/YFXkTyrKn4MvBn4\n2IJZYpqtChwJ/A+w88BfvwuciPS9sc+rRRZaImxRH2ZrAC8BThn4y7eBU5C+M/5JzSLM7kcSGL7f\nwF9eDbxswdwEFhBaImzhQ1omvgF4fOXVm4HjSR5fuzQEMFsGPA14PbB65S/vAE5qT6nnBloibDFz\nmK1G8nSeW3n158AxSN+fnUnNM5jtC3wI2Kry6muBU5H+OzuTatESYYvpYfZ44OzKK78GHoX0i1ma\n0cKA2b2AT5DiIUs8HOkLszSjRYuWCFsMh9kGwKeBB1ZebS/SUcHsMaTT8xJfAw5HunGWZrSo0Aqz\ntpgIswNJFb3+TiLBdwOrIVlLgiOE9DHShbgGcBbwIOAGzITZfWd3ci2GoSXChQizpxcEeEHxyoEF\n+R2H9J/ZnNqignQr0rEFKT68ePW7BSEeO5tTazERLREuFJgtweyNBQG+HbgS2LwgwK/O7uRaIH2h\nIMQ7AVcD7yoIsVPkX7eYRbREON9hZpi9C7iddAL8ZWAtpG2Qrp3dybVYCdKfkbYk5Uh/D1gB3IHZ\n21pCnD20hyXzFemieSVwUvHK+4CnzYeYP+vaRsA+pDS9exY/dwNWdZq8haQ8U/5cBvxAHf21/mxH\njJTN8hHg8OKVlyC9YhZnNO/RnhovFpi9ADijePZ54Aik22ZxRivBurY+cBgpN/mgGXb7K3B95efv\nJPWZktzLA4gNgI2BTYA7F79ngi+QQlzOUWeOneCmYO0vAA8uXnkK0gdmcUbzFi0RLnSkdK8yxe07\npEOQW2dxRgBY1zYlZVw8G9hsSJPbSUv2bxQ/P1FHt494TkuBXYBAUp85eJKmfwbeArx7TniQZusA\nF9OPR7xPm9+ch5YIFyrMNgaqai6boNm7aK1rWwCnAicM+fNXScv0c9TRzWOd2AxhXVsbeARwNEkv\ncRDXAAeqo1+NdWJVmG1GImmA/wPu0qbwzQwtES40pH3AT9LfP9p3NkQPrGurAM8E3jjwp1uALvDO\nObfUzIR1bT0Ssb96yJ/fAjx/1F7sUJgdAHylePZ2pGeOfQ7zDC0RLiSY7UcSMYWUu/rKsQ7ftXWB\n15CEF6p4BfA6dXTDOOczbljX7gycx8pyXD8D9h37/2/2JuA5xbN2uTwFWiJcCEhiCL8FtiDJ1m+G\n9K+xDN211UlE94LKy1cBx6qjrwzvtfBhXTPSe/K6gT9dBOyvzpj2adP+4Q2ksLcfAvvMhyiBcaMl\nwvkOs0eRcoIBDkH60liG7drDgXMqL/0BOEId/XAc4883WNceB3x44OVnqqO3j2cCE3KZD2wD5iei\nJcL5iiSB/2tSjY2fALuN+k5fhLd8lImhLYepo8+PctyFhMJTfBXw4oE/bTLyE+i0crgKWA78CrhH\n6x0mtEQ4H2G2F1Dq/43cC7Su7UFaVpX4CHCcOvrnKMdd6CiCxAfJ70B1RuytmT0MOLd4tjPSz0Y6\n3jxAS4TzDWYfBR4D/AdYb5QxgUOWc4eoM56l92KDda1DSp0rcYo6etXoBrQ1gfJG9mak507VfKGj\nJcL5ArO1SQchAM9CetvIhuraccCZxdMbgN3V0e9HNV6LPqw7IQAe4K3q6NmjG3ACAa+5WOuptEQ4\nHzAxO2RbpCtGMkzXTiDV0gD4DbC3OvrHKMZqMTWsa1uR1IBKnKGOThzNYLY9qZA9wK5IPx7JOHMY\nLRHOdZi9Hng+aS9p+Sg2t61rjwA+Vzz9FXCfhR7zN19gXbsT8MfKS8eqo7OaH8hWAcqaKS9Gek3j\nY8xhtEQ4V5EyRP4CbASchvS/jQ/RtXsCPy2e/gXYUR39velxmoAZ6wGPIlXNe1DD5r9IElb4nMSc\nTEmz7gSvDWAfdUZQMMvsraSMoN8AOzBbF/uY0RLhXEQKgi0vyMYzAqxrawK/oy92sIM6unyKLmOD\nGTuRKuU9LLPr34ArSN7TDcBt9JVn1iOFjGxNurHk4JPAyRK/y+w3EljX9ifVMymxduOn92YPIuV/\nwyLZN2yJcK7BbEfgl8Wz5UjXN2q+a13gpcXTQ9XRuVO1HyXMWI0UU/f8KZpdRcpX/qDESL3Vwut8\nEvBCJpbbHMQZwKkS/x7lfKaCde1lwGnF03PV0aHNDmB3JglJAGyJdHWj9ucYWiKcS7AJe3XLmtQL\ntK7dnSRACqmu7pPVGe8HaYYBLwFeNkmT04GXSYwlPXCmMGNd0s3jhZM0eR7wZonxvp8pOLu6Z7yn\nOrq4uQFsNaAMz3oA0rcbsz3H0BLhXIHZc0mez7eRHtCY2e5KajTL1WnWy5xyfGNtUkGo+wz580Ok\nXrGoeQUzDiGJog7ie8ABEmMLNreu3YMk7ABJlHbtxm5yaa/6CpKHfCwawUHNHEBbznMuwOwNJBJ8\ne8MkuDPJYzgceJ46snGQoBmrm/EDM0SKfSxJ8BXAEgkrfuYlCQJIfLH8P4BV6Isr7APcbIbMiMXy\nf7Rz6ejn6shIhz5rAncU8YgNGJeQtgY+CLwbs1MbsbsA0HqETcLsHFLZxhchDSqV+M127f3Ak4un\nG4wjHMaMU0hkV8UTJc4e9dhzCWY8FXjvwMvPl1bSZmx+7K5tDvypePpLdbRTc8btlcDJwFuQnjNd\n8/mEdmk8mzC7kOQtPRbpY42YnJi/eqI6OmOq9rXHM4bly54iMbr0sHmCYk+0WiyrxAYSI70xWde+\nDdy/eNrcdojZ84A3AJ9CenQjNucAWiKcLZh9D9ibBkUTrGtHAh8vnt5FHV0zVftaYxnVEAtIytOb\nzdU4vNmGGesD1zGx6t79JC4c2ZgTxTJeqI5e34xhewLpwO2bSKERm7OMdo9wNmD2QxIJHtwgCX6R\nRILfBJaMigTNeGax91eS4Ipir2zNlgQnh8QNEqsVe4qlcvh3i73EJ0/V1z1mOkEur9czrGvNeDDS\n2cChwH6Yzdt93rpoPcI6MLsY2B04COn82ua6E1REnqCOBsU/G4EZ5al2iQdKvbIALRww46FA9Ub4\nLImRiGlY195Bv3DWOo0UyeqLAp+D9Ija9mYR7dJ4nDD7Kik97CFIte+k1rVdSKKsANuq07wYgxnV\n5TbAzhKLXr+uSZixO6kcZ4nDJT7T+DjdCTqWe6ujH9Q32lO+/jjSY2rbmyW0S+NxwexMEgk+qiES\nPIo+CS5rmgTN2KVYApckuHuxBG5JsGFIXFIsmfcsXvp0sWTesdFxEvGtVTz9vnVtqoyeGRrVx4Cn\nAEdh9tra9uYRWo8wF2Ynk/aFGtERtK69BngRcJ46mqwIuc+2sQZMyOw4VGLW0vAWI8w4CqhGEawm\n8Z/G7E/MSInq6IH1jdpJpHTJ45HeVdvemOHilhRjOf6fNPTsjO3+gScIJDi9kfdgBd9lBWIFpzb/\n/ur9aaoS6IxZf+8mfvZG8qi/D6jhn0sK2zbb/+fA5/G+yudxZuP2V3BR8V1SIzbh/cVkHzzb753j\n+5X9HrQe4Uxhdl/gu8CnkY6oZSrdxf8FrE7DYglm7EDSISyxisT4C5P35mM7khRW7pTZ9Q8k+ajr\nSAo0/wE2KOxsSyp0lWtvf0mzpjxjRlUnEOCuEo2phQ+INyypnZpndgmwG7AL0k+naz5X0B6WjApm\ny4FrgeuQNpuu+ZSmuraUJCsFsJs6+lHd6fVsG78A7l483Vui/gZ69hxs8PR0EP8BngZ8RNJICNrM\nlgBHAmcDS6do+nBJw3KMRwozypsqwIUSzaTQsVL86VJ1aor/Wi9MZyM0N/UtB9ES4SgwUe13CTXe\nMOvaMujtD22nTjPeiRnbQk9j7xKJPZqwO7OxbVXgUuAeQ/78N2A3SX8Y13ymgpndhSSkMMw7/Tmw\nq6T/DvnbiObD70jeLcBdJBqJFx2olbKqOjX+p3RTKW9YS5kHJUPbU+PRoPwSrdMgCW7WIAm+hz4J\n7jQOEjSzVc3sRkvewr/pk+C/geWSrPjZeK6QIICkqyXduZwfsCn05LbuAfzHzGRmN5vZstHPh7tC\n7/O62oy3NGK3o+8COxdP/1N895zGdAd9Ady/1JzanEVLhFPBrAySvjvyB61a11ahT4Ibq6Pr6k+N\npUVIzNFQ7vf2RGBHAjN7a4X81i1efi+wpCCX1dWw+OwoIekvkpYUpLiEpMoCKSylJMWRiitIXEL/\nOnxWEWoz1XJ+ZnY7+hlwt+JpXTL8O0mJZ0PM3jFd8/mIlggng9mTgAcDj0H61XTNJzWT9gRLr3K5\nOvpb/amxG/19xqdIo/sczWy9ghBEqn8B8A365HeMZmt/pUEo4ckVUvxm8afnlv+/ma03mrHTjQw4\nrnjpNrOhWw15djv6DbB98fQ/xXfRaUzfJ6mPn4DZI+vOba6h3SMchrSX9Afgi0i59Tb6ZibGeG2u\njq6tPzVeR19deUOJkZToNLN9gW8NvLx8Pnl8TcAmytyX2E/S4HvT0HhsTH8JeppE7UJf1rWdSHug\nUPc0ua+ytAnSoFLRnEB7WNIEqpvDNedXSYzfUp36dSLM+DeF4knhQTQOMzuCpIJd4lJgj4Xg9dWF\n9cNJShwuqfH0uTRWb+/yOolakQoA1rX7QFLHKYRfvROr3txrHR6OCu1hSTP4bfG71jLIunZZ8XDn\nuiRohhUXxqrA20ZBgmb2sGL5W5LgqcXSd/eWBBOK98KAFcVLny6WzAc1PxZGCv9ZXuwb1vrM1dFF\nFCUeainXpO/C5sWzL9eZ01zCSInQzK40s8vM7EdmNvaYtmyYHQ1sAxyA5Jahsq59gHRqd0CxaV1j\nSqxL/w58qMSz6thb2b5tWxBgGdT9woIAXzlVv8UMSd2CEE8sXjqvIMStmx2HJwJl8P4dZr3cYp+9\njj5DSufEuvbrGhO7FngC8GDMDqszp7mCkS6NzewKYHcNCcScc0tjsw1JcW9fQjrEbaZrzwDeBhyv\nTr08zYH4wO3UYD1emxgfBnCmpBMma99icpjZWcAxlZeWqsF4u4FsoS0k/ljLXtc+AzwSeL06mqya\n30wmVsrQrYt0U505NYk5t0dYEOEe0sonpXOQCMs3wr3vYV3bmxSw+3519NR606Eqs7SuRGNfNLMJ\nNVBEunDb5W9NmE1Ycr5b0nGTNs62zfrQOxjbRaJWyltlefxQdeRb4lb3C+fQtTwX9wgFfNXMLjaz\nY0c8lh9mZS2QbWuQ4EYkEryhARJ8MH0SXKUpEjSzzYqLtSTBdYs4upYEG0Bx8ZV7y8cWy+Xlzdjm\nBuhV0bvMrF5aXuXA5Dzr2ibOSQnYFQCzZ9eZz2xj1B7h5pL+bGabAF8Bnq2isPSc8QjNtgKuBF6J\n5CpvOBAmUys8wYyDSaUcky01U2jczL4DvYvnhVJDNS+ciNE2J8lTNVbytMC3gCNC0KxmQZjZi4DX\nFE+/KunAZuxS3dLYX+Ibblvdiemj7u+t2adIBzFzIh95zi2NJwxk1gFullIltsIz6VaaRElxLJOZ\nOLH0BtQgZevaz0gpWhvXCZg24wDSDaOx8BgzWwuoZsWM3QOMcWgs3nS4kCRM8AvSknBVUkrc9qQ4\ntr0ybN0GbB7CeOPebGKoCcCakm6pb5eq3bpk2Nt/dIfVzPIS2cwCECovdeYMEZrZmqS9p5uKi/EC\noKtC0XlOeIRmZe3erXDmxFrXngm8FThYHZ3nnwr7ABcVTxvxBM3seOCdxdOTJb26rs2ZIEZbBvwR\nmGrJtVsIzSnvDIxf3V8dhl8DO4UwHgEB63/PAI5XA2KnA2S4m4T7vbRuT2z4berIF5VgdjfS+/oc\npEZypr2YUx6hmW0DfLZ4ugrwYUmvqvx9donQbF3gRuDDSE9wmejadiTNvLPV0RP9U2FH6OUJN0WC\nVRtrS/rnpI0bQIwreZ5VHBmCPjnJ38aCGG1YofYSq4WgxlSjh8HM1ob+Xm8T3/0BMtxG4kq3rf7h\nyS7qOLUHzc4DDgJWR/q3dy51MaeIcNqBZ58Ia50SV/cF60Tqm7EJUKatLZWo5aXYRNkwSRrpgViM\n9i1g34GXfwjsHcLcPISJ0ZYAfyYttav4Ygj+lMqZYOAGtYpqajIOiL2upxplWCtk6Nsv7IdkXYW0\ntXcedeHhllVGNZk5jSQeCrB3jRShMkB83SlbTTkNltEnwTUaIMFqvNnzJL2pjr3JEKPdCVaKZbsc\n2HGukl8VxZJ4OfRI8QZgHeCQGHtksF4I/qD6ySDJzOyFwOuA28xsO9VQzZa4rQi6/z/gRrNaiuRb\nkPZy7wDHHrV0B2YPB87FbB+k7znnMXYsTo8w3ZX/gbShq3u3p8J8hDr6tH8avSXw5hK1BBnM7HFA\nWQd5G0lX1rE3DDFagJU25u8cgv7U9FizgRh7EQRV7BqCftz0WMXWUSnTf5SkT9SzR89enYM269oZ\nwAuAp6ijDzgnU/sAsg7apfHMBi6rxq3nSaOzrq1K0uP7rTrafrr2k0+jR4IPkPi2106yZe8Eji+e\nLpN021TtcxGj7U+qO1LithA0cuHS2UKMK532Auwegi5tcpyBbYw3SXpePXs8hCL/tyYZlt/N1dRx\n7J2mcLnrgRcyC2FaLRFOP+jqwC3AB5GePF3zoSa6di1pWbWKOr79HTM+DDwOWCFNCCFy2OqlOUHD\noTExTvBaAL4VgvZryv58QIx2FbBl5aXNQqgvrFtFZd/wB5L2rmeLVwInA++SejfHPBtd65WBrRFS\n8zVgf2ZB3r8lwukH/TmwE/4DkgcD5wMHqqOv+qZAaeMv0kqb9Zm27F/AGtDMKWSJYt+sSvJXhTB7\nm99zATFaTwKtwJIm90MrZPg3SRvXs8WtpCwU92rDur2A8Keqo/c7JlF6u59FepRnDl60RDj1gJsD\nfwKehPSh7O7dHjlcq442n6798CnQC6GoGzBdPX1smAQ/C1QVRRq94OczhiyZ3xmCnt6U/cpn+h9J\nq03ZeFpbva2X1SVcoSyVJbKvGp7ZS0lJExsg3eCZgwctEU494C2k+CbXmNa1COwHrK6OL0aqiS9n\nsmP/pTjxb+o9jLGnvlPiTiHoz03YXmgYsmWwRgi6tQnbFTK8VdIafjv04jq9N13r2nqkE/U6S2QB\n1yDdxdXfNeTcE12YGzDbllRM3SWvZV3bhkSCT6tBgqUu4QNqkuCvaZ4Ev0ifBM8OQdaS4OQIQVeE\nICOJbADcEqO9c6o+M0XlM13dzNzZIhL/BB4OYMbXXTY6upFCp9K6drdpmk+GJwJbFAcocxaLwyOs\neZxfLhG8d0UzHkAqBvQdaaXg4ww79lHgMdAMCQ5Z6jXm2SwWxNjLUAKgIMjaqHiGH5D0FL8dfgns\nCOwt4RJHrvv9L66/a5FvSyl/uHZpPGygsnDNfjgK7ljXHg18Ari7OvnV7KppULVCGqyX09wUCe4O\nXFw8/V0I2q6uzew5WPSIMUyFnwG7B4WRpssNQyUQG2CHEHR5XZsVMjxO0rv9dnpbMq70zUrc7AvU\n0RscE3gs8BFgM9Tsifvw4VoiHDaQ2xuspNH9Th0fUVS+hJtKvgLZZrYnRSZLQyRYyiYB7BWCfljX\n5ozGtXgS8KppG/bxUdJN6A5Sjd7tgV1INXZnisODwkgKLA0iRuupBwGvCUEvrmuzQoa7ST6RCjPK\ng8I6+4V10+8E/BVp5EvklghXHmQL4Gpgf6RsqSLr2huB5wJrqpMvn2TGIcAXgLdKuIQrCxWfUjCh\ndpzggOcy0hPhaHEp/frLw3BUUKiVUTEwnpFCPk6cpMntwLKgMLr/OU4Ikq69VB6Q8lpdTjEDMz5E\nqjMSpF7N5pn379qmwHXAN9VRcEzgKJL+5Mg1C1siXHmQ64BNnd5g+YV2Kcs0uCQuP6BNVLOObJUE\nm9rLGjqOxTfDUOK/a1D4/ZDXR4Zo8UHAsJjP84NC49XneuM2+F4XKtfXQr0VQQNL5LKPL5kgfZcv\nRtozu2/WMK3oQh9m65DURbyy+WUOsSsDhX59CbdUe4UEH12HBAcORS4PQTt4bU06hsWq51rFsqDQ\naMpfDoLC1yDdiKLFqljEQ6LFHjE07SWGIIvR7gAsRlMdMpR0nSWP6uM1HYitgKuAy0hVFnOxKvAf\n+iIVuXg5cBpmq9BwGmhdLFyP0OybwAOc3uBqwK3Aq9TRKflD9xLgvy7xoNz+yYaVXtWlknafrv1k\nGMgSeUMIeoHX1lD7FjcABpc6hwSFLzU5TtOIFo8Bzhp4ed2g0Gg1thjtHIowFmpuRZj1lNBfK+l/\nfDa4nLTXeieJ7BAp69qVJELNz0PuL/PfiZoLRB8yTLs0Lo2T3vA3Iz03u3vXLgAOxLkxXC5B3BvT\nZncB/pBs+N+jARJ8Vgh6m9fWSrYt9vJRK1glKNTS1xs3osXBQHKApUGhsfzYGK3M/4X6ZFhRLJJL\nsajO97OSYeULpynFW0foBLUB1X2U+1PZNVuLvcEDgdOdJHh08fCRuX0rKMsGuLcuiuVwSUpPbpgE\nxUQSXBoUbL6RIEBQ+HtQMIqc7QK3V5bN9ccIOgUoPbi6BFt+J+oEvD8BwIzH53YsUu2ugN61kotU\nsN7Mu+U0EixUj1DAdUibZXft2tnA45k9b/CvwEbAwZK/Bkpls/6YEDSZRH2eTYulhFmJ1YPCrEmy\njwID+4gAjwgK5zRiO1oHWAH1DlDM7JHAZ8C/YmjIK7xOnfxrDLPbSKo0I7n+W48QIKk0g6NMZBE3\n+Hjgg04SLIVRt8rtm/rbvUgk+LOGSHBFgyQo+iT4iMIDXFAkCBAU/lR4iGWt68835R2GoC7wblgp\njCkLkspaQJiZd//4bql/CtLPGr8vwLC8uGZycX/S4LXUl5rEwvMIrdAL9B2SvAw4DVimTt6pVhPh\nMuX+T819wfIDPS8EHey107NnscysAaAgiUWDARLcNyh8p7bNaFcAW0Ntz7DW96WmV1iKOnxUHT3O\nMbiAbyM1Xde69QiLQ5LlgDei/zTgklwSLHBZ8Xstz8BmFouHO3n6A8Rol5SPGyJB0SfBBy02EoQe\n8Zf7Wd9uwjsMQduUj2O0j9cw9SAAM/MGKG+Y+idV6xyo06uK+Fjn2GexctGvWcPC8gjNjgPOxKGK\na10rA283UScvZq+uN2jWk8Byi3LGaI8APgfNBEtXL/jFSICDiBYniNU28Z5UvPd9QtBUdZgnReUU\neUNJ/5iy8dD+/iBr69o9gZ8Cj1ZHn8ocuFSLPwzp81l9pzW92MNn6uUVuxU2zLgE2A1YTSI74b/y\nZXal0MXYi3uEmuEZ0eL69IPBZ5UELUYjFYnfi/T/XQr8Q2F0KXLTYcAjrBVmMxDovtRTcL6ague5\nnipiwRdIPCS7fx1lmhEVeVrcmSXpDgMOzcFCgBJS2IwHuwE4SbCso/uMGnnEJQnepSYJ3g34dfH0\nBUEhX2kkExZXOqWdSZ/J/rSxQhiMCWwUQcGixYtIwg+3R4sbBAWX+nIIUoy2C2lb5XbIX01Ikpl9\nCHiimT1KUpbAhMTNxXHHg3PHLnAa8HLr2vrqZKtQHwe8C7Ml465rMoiFtEd4OgCSJ6PhAwCeOiRm\nvK94uN6UDSfHuQCS3uHpXFlevTcEuSWtosU96ZPg3UdJghajyh+mJsEVwEYKwYClwDbAVPP664Dt\nkSAo3AcoK879I1q8l9tW0E8pKgTGaNkybwCSnlQ89JaW3QjAjDdnj93R/xYPPd+99xS/R5ZlMlMs\nnKVxcrN/g5StpFu49+eqo0Pzh61x8mbWBV4K3FvST3L7VzUF6+wLDggTbB4UatVYHoYpiOlXCuHu\nDdg3Uh7tMEn4WxWCW/Z+MkSLu5KW6wDbBIUr3bb6N7QtQlCWhwxgZqXS0VmSjs3vX+sEeU4tjxfv\nHqH1yg/ug/I2na1r9wEuBDZUJ2+z2YxSWugeEr/I6Zv61wt/qFw87n3BwhMslYvXCwrZtZ4ng8U4\nWA2vxCYKoZaSzgzG7im2DGBJk3uM0eKOwC+LpxsH+Zbm1XRI702tzvfJjPKmepDE+Vl9u7Y1Kdvk\nYHUy41/NDgM+CyxrSohhMRNhiv/zHZLcDKzlPCSp4w2+nbQkuJOUXx+kQoJPDcFRbpGVBBPWDgrD\n1GOyUXiww34rAAAgAElEQVRng3s+z1YI2cG7Dc2n9LyraIwQo8V9gIuKp6t5FbJjtJcDLwGurIbY\nzBRmdm/gR8CnJD06v/8seIX9w56nIzVS92UxE6GAW3FU/So+wBero9fkDclGwF+BF0m8LnvcGnfv\nGG0T4Hrwew8D4SDbBYXfeewMYsgS+E4KYU4UgrIYtwKurL5W7D/WRrR4KPB5qHfSXrnBrR5Cvghr\nTa/wdFJO9NpF8aeZ9+3apcCuc2F5vDgDqq2X4pOtO2hdKwtPnzFlw+H4K4CTBF9RPLyzY1woSJB6\np/4lCR7ZBAlajM8dIMHtFYLNFRIEUAhXFcS3a/lacbByTF3bRT7y+2ClEJtclAH53iJauwGYWXZq\npdRLRLjKMe4DAKxrxzn6nubo0yjmv0do9iTSqe8SMv8Z69p/SOl0414W1/EGy6DxH4egXadrP9RG\n/0JtRKV5gAB/rBBc8xo3LMYrqeSFN+EdNhGIXvEKDwpBWft1UNsrnI3lcRlcnb3HP9zcYvQI4fUA\nuSRYYBl9nbgZw4wy9m/L/L5WinR65crPBKhBgr0tgLokaDFuWiXBwgOcFyQIoBC2pnINFN7h2nVs\nVskvWnyWy0Z/uyM79a3AYQBmdqSj7w6pb75oCfBdRx9Qr4SsK4SsCSwEItyIIg4wB9a1Mnjasywu\nYv+42tH3nNRXF0/XcBAx2rnFw6OnbDhZ/3Q48iKonzFiMb6DVMwHaG6vbdxQCBqY+00W46l1bFbe\n27cUBaU8eB30DlCyoH7KWnYes0RZhvSbuX2Bg4CyBG4uzqeyZTFuzG8iTHVJALqO3u8CUEf/na7h\nJMgulm3Wy2DxzBdInmgIet90DSdBeUK8obM/0FsKn1A83Wm+kmAVxf9QekH/20BAdultuzImQlAp\nefYS5/gRwMzWdfT9w/RNVoY6url46KlM2AHAXGKvtTG/iRCeD4B0haPv1tDTD5wxzCiXO8Ex5nUA\nklbkdozRvlg8zE4hBIgWf148vDAoZCfnlxggiCUK4ZeTNp5nUAjfpnIAVYcMg0Jvby9a/JzTzNEA\nMVpWREOB/YvfNzr63hvArFCTHg9Kx8KznK+N+X1Y4jx2t65tDPwFuIs6eWlps3hIIvCFy0SLy0jV\nx2otiQf3A7125gOa+l/Lw5O6Byeez302Dk0qCQobq6O84PI0398ibZ/VbyUzi/mwJA+nAuSSYAU/\ny+1gZuVSKTtQNkYrl0lPy+1boAzw3cTZf3ZI0GwpZi/HTMXPyZgtG8fQ1f+x5jL5/lArpOYN0JNZ\ny8W+AGaWrSpDchSyoY7KwPJPOrp/CdjOM25dzF+PsJ9WtyVS1qFFccx/hzpamjck+5H2XjaRyNQs\nnDVvsMxBdXslFuO3KC6qkZFgX4rMh1HVv+iT4EcUQnaxI5hAgq585NnwCs3YArga2FMi62CvRhjN\n3sD3gFWQo4B8z8zi8gjTlzKTBCvIPo2j2IDOJcEKsvfTYrQy6Nq7F1eXBHenryTc7Pel7+lNRoLX\nAWsiWUF0awHDc6Gnt+VF+T8/zmLc1mmjvOF69rJ7iNHWrNM/B1JPTeaHju6XTd9kKMp9wuDs78Z8\nJsLsesUA1rVyefiWBucy9ZhmRxUP7+Pofg1ACMqW8I8Wy/zaOkGqpTewSWNiBcMJ6xk9wuv/bIZ0\nS6+F9C+k9VZqB2+fgX0Xiv+53M5wZeBUxVujxQMcJsrr1JMLfhcAM3N5s048CsC6tnlWr/7y9IQp\n240A85kI7wl5LnuB5wLkbuSa9e7qnhSijwFI8pzg1UEXICjs4+lcWRae3YhazDCC6hOaP5hWemaF\nFKcez2M+hCspvBXvfmHFI/9Kdt8aYrtSbx/8bEf3k4GyFMXMx+yovGG82zHmrTDW02pgfhMh+CLR\nXzR9k6E4pfh9lrN/dgJ9jFbGkGWXbIwWSw/gg7l9ASzGXkEghfBEj42+MTt0gJA+MZS4msDkhLh1\nLbMh7N0zF+NFU7WdDtHivR3dHg4QoytbxIsybMebi+0J9WpEgSYX85MIzco0qCxZ8gKrAh919HsZ\nQHaBm/4F6ImafzlACLp0uoZDcDZAUHjydA0HUchobQANHI4kAuwX50lEddTkHRrCyoR4RV3vsPJe\nuDzsilf4o+y+QV8oHnqq3h0IYGZb5HSSesHgHs/OizOBsQdWz08iLO6OKLtGQokzm5vKtPgagKSx\nBR5Hi65KeBWUF4AnZKOPicSz8ahOdqdEGrPnzTWwVD4eaofUUCP1LhtSrwSF92DRg1OmbzIEUlmu\noPF6x1NhvhKhq5aqda2sLfEt57geuSDXSWOMVua7DpOenw5lDFj252sx9rwdhXCOY+zCUIVwknc2\n0qJKU0L6wQQSrkGGCuFdPTMxegJ/y9WMJ/UuxSTGnnDHOPAmZ78zAKxr2ds6BR41fZPmMF+J8OGQ\nt0Qt8FQAdTLluqx32utJdYIUKJqL/wWoU5ApyHXKW+5/+b8bgyQ4V9AQGdIPh7l8ylZDUEcFPASV\n6i7+G1Q+XgJgxi45ndRRGcTvETWBlghnjI85+jzOOdbbIL9cp5mVUfJPcY6bjWjxhcXDbIkti/Hw\n8rE7VGaukmCJBshQoR8OU/WgM3AuQLS473QNG8QuAGZ215xOEqWQgjfcbD9Hnx8DeaE3NTGfidDj\nZW2CL+TGKw/0HgBJWelKMfYUQwbrbMwEr4OJSf8Z+BTUOCCZ6yRYohnPsLx2sk+Qg0JZLdGzRfNl\n6BV7mjEk/bR46D3xHuee3WfHOBYwH4nQel+AbzgtfH76Jo3B++Up65F4sl9csBhrSXNh1s9AmMsk\nWGIiGR6f3b3iMVuMqzU0q5ng4OL315z93fnmDlzi7OcVpHVj/hEhpAwL5dd+LeAlQq8GoAeuCyta\nLAUxPQRcHmZk5V9XsEfxez59p0o9S2/s2urFb099kcsAosXlOZ0qwdXBMaYXX5y+yVB4BFGgTNEz\n28DZPxvz6Utb4oGeTta1sjZFlnKMWY+UXuUZF3i1s58HnwAICt/2Gqjuf80YE5fEjeX6WtfMunaN\ndU3Wtf9a1w6xrjXnbaonJOpaIiuE7CD5EkGhjGAYVnt5VPCmWr4KwCybLz4LYF3Li3zoS/d7y1lk\nY9EQIcWmbe6JMXAUgMRvcjqZ9fIss06aY+xd6DGnXx1Y7NUxyVcWNuuHkDSwJC5IT6VCEP1Kf6uQ\nBCTuqLR5Wd3xBpbIHm/4SoAmKuGNAU8AMLP1M/tdWPzOuvbU6eWJPyZzvBL3d/bLxnwkwvs5+93X\n2e9YZ7+U0yzlqkE/svg9zlixFwEoBE/GRxlCUqvoUYX8hmEysYPTpuk3U4Ti9225HRVCKcjgTb30\nYHuAGC0rjlHSb4uHr8zr1wtVe0ZOvwq8KZreaz0bs1IfoCY2pV97Iwce5Rfw35Vc6jjApwFCqCzb\nZoBo8enFw/Fthptt1HssueLjiqXuhOX4THTsrGvbUiFI65pcxcUBpG9SOuJm1uTyfhqcBTwtWlwz\nKPxrpp1C0G9TFiTfAbL2GAs8HR+peWP7dnb28zov2ZiPHiHALxx9dsFxx6+B1adv0ijeDhCUpxJj\nsVY6XjmWyxu0ru3IRBJcfaZkpo5+X7Tt1xhJ3qFXs68UCHAVW4L80+OgUK42vEHWmzr7zQdczRiv\noflKhN68Xe9x/kKGOx2vB4c3aF1bQuVzVEemjrIPH9TR7QPk+U/XgYrkiUstUe4tek6PWwxH1p58\nXcxXIvR4hHX6ZS1TK7jJ2W/syM4kMfPUva2iJ8XuXtJWMGDD7dUBYJaVLeI6aW8xHX47fZPmMF+J\n0OsRet/cr07fZCi8ZRy9MZLjRIpV9FRIqxxuNEGCw2y5DlD6/0stvcE5Dm/Iz4XTNxkKTyYXtB7h\njOCt/eCSWidtSnvgjf5/q7NftvrxuFEsiYFmSXCUNucovGT9IWe/707fZCg8NU+g9QhnBFepQeAq\nZz/vXS0rgDXGnhhllqx6tFjGhWWdVBcCrOD/snrgrk6WiwbCavLHjNk6g2+fvslQeAP8PWLG4BCT\nLfDT6ZsMxR+c/VyYX0TYr2nrrf1xnbPfr6Zv0oeZrVc8zJVp2gdc0lsnAgSF3C2DMrUsL0zBrAyQ\n9aglA6P13BqxnR9cXeZq58bMrQCIFnNPgL8EEKPl6l1+D8AsT7QBRy3vAt4l7ljr+4yMCM3sIDP7\nlZn9xsxe3JDZFLcmeTensxStzVireHh95jjbAyh/ng/KbF/i1OmbDMVxAAohN6wolTqQsjIGrGve\nWM5xogyByXpPFEIZOP+BnH5BoVzdPDWrX+jV/c2qC1IJ8M+S46LYjqoUMZspvHqa6YDSGkypnAIj\nIUJLd9O3kjTxdgIea2Z3b8D0RtM3GTKffjjF8Jq4k2M55NcpAbJqQ1TgKeozn1BuuO8/rgGzl8dS\nluZkg3iWs9/e0zcZih1yGld0CXNViq6HCdfgTFEe6iybslVDGJVHuBfwW0lXSvovSUS1Xv2LhDWc\n/daGFHOW2c+rfuGJ9gefLP+8gzrySqjlYL3pm8wpeG+ennIB4Bc+XWf6JhNQ5huvmtmvJMLcfi6M\nigjvzMRCMdfQT56vA+/m91rTNxkKL/F6sxvm28U7FphxX7O8C1cd5Xr/8xWbOft5b9a52R7l9lCu\nZ1d65mMhwlHlGs+IsMxsReVplBSn6eLdG/SeVHrfH+++xvw6vBofngP8Hm9ltIUNr3Pg/W7nfkfL\n+eVdE9LtRf73tERoZoGa+oyjIsI/MnGZdxeGbJpKWpFp10uE3iDS/zr7eXOaZ5x4v5gguWWcFgO8\nYS3eLKvcg0MA1JEny+pGZnDtFg5ULJ+bWSd3oFF5IBcD25vZ1ma2KknTr4nKW7WI0LFhe8v0TYYi\nV3qrhDe8Z16hUXHVycfwbmvMFrx5ylkiGxV4tw5yrwmv4jkkR83rjGRhJEQo6TbSKdj5pDvPxxsq\ncO5d4pYeWq4H7P2SeUMGxhpNP4sYR25u6V17P8Nxw5vC5vXsspSxzXrXTi4R1tGpXMaYFKNGticl\n6TxJO0jaTpI3Cn4QHh3C6mlx7onXnwDMsg9NvFHxlzr7ZQV8NwazrBvLbKS/qaM8fUazI4tH3tjX\nrDzxaLF8D1/nHC8r/dP6geK56aYbAkjZxFRH5m0Z89kjHCFS8Gk/wyQXWSdslQ89N/j0KgAzyyXe\nCyC/VCNFhki07PSu3QEsxt0y+5Xv/1i+pLmwru1eo3vKlpGySixYjGV1udw6G2XBraw88RjtTsXD\n3D3CLQGk7D07b7iNTzMxbandjrJD3lyYX0TYf1PGHad3r+mb9FFsDUBBNBkoc6FzVbHfVvw+cspW\nA1AIpQeap9PY///cGHEecJkb3tRKZCb4IoBC+HNmv7cABGVn9xwOEEL2Z7HH9E2GwnvteK/V9Rlj\nmt38IsI+vHcn74fpraaVRWiVUo1ZdXYrF9HHcvrNBgaksrxV1SbFgMTXfAi3cWVLAU9z9vOmcXoD\nt7d09luPzJTYOpivROgNIt1m+iZD4a1b8lBnv8c5+40faQmTi/Ii3qvRqXQnbEW8IN+A7dvcbEaO\nXZz9DnD2u4ezX1YqXwWtRzgDbDV9k6HwfpjePaexFZ+ZBZTZOh55/feUjxteIvdCQtTRGxz9v5U6\nL2hNw9z97hK5+8gldnL2W4+WCKeF9264a6OzWBj4LYDFmOctS7WCvwfVpK1r+3ltWdeeMyrV6xmN\n3z9setM4xx0zvNeOlwhbj3AG8Ki03IB/v2LBQiGUez+/dxsxcwWeDxBW9HiHRZ8eAblJ0Gp5ppcA\nKITn5XSKFstthfqF6scDTxbL+vhKTyxnjAkGi4kIvRXsvOl5Xmlzr0z/ugDR4snO/vnoLyHdZRcH\niass2G5dm1SNxbp25JDC7hc14gmOd1l8DUBQyEoJi7GXmeOV7vq8s5+3ZMX3HH3ujD8xIRvzkQh/\nhU+j7AfO8d7l7Pd6Z7//AYjRHpDTKSiUcWGvdI6Lxej/PtTwqIpSnoMEdHWFFCf8MKCMXfT378fW\nmHtuLeMB5AV791GKuGbJ/Ftf5PTVznHjGPvdmTEWMZuPRPhjZz9v+cm3AJhlywF9LvWzrAs0BJXL\nzLrlMnOwbvE7P3i1QQ+qQog/n0HzR0xCoDUm4LJV5giPSsBkGN4DE8KtZopDASRleWhmPVk5b3Gw\n6OizBWMkwnF+eE3hQnCpkXwLwLq2tjqacZ1iid8U99HHAB+ceT/dUdyAX4M//MaFaNGCZl6nWCHc\nZDHWH9hMTRCjOrpn/cnMEH1vsJYoiELIuolEi+V2QhNiJDPFa539DgCQ8up0W9fKeF9PPnTrEU6D\n84HsWgbq9DytBzrHPc3Z737Ofh6UCe7u3GOL8UXZnarkZzZbUvf5qC6JpWwFdYuxzFy5esqGw1F+\nHw9z9PXCGxT9KGe/AwHUyazdk67tLWj3CKdEWRXrbs7+D3P2287Zz4PDAWK0rFIBQeGfxcPs90Yh\nlGSWlWPbN9Ajw2WYNaFGPlqY9YPd/V7sSQAKwR2NkOO5A8Ro5ffQG9fnwROc/Q519tsc+CdSGz4z\nKdTbFznIacF7d/PgfzydQlBZe/bXDc5lxrAYvQWByv/3mnFVH3PBbC2Kcpg4C2ZZjA/2Dh8tlp67\nJyLhNwAhKCuUpXJQklUtr8BSij3vTBzu6AMpGyW3FG4tzD8i7MNDhN/BJwvkDZR9A4CZee+MnlPF\nLQGiRY8OX/l98IQ7gFTdg7pjTpJhUi4q94h/hPQTp6XzYYInnYObAIKCO+zIgTJ/PavcaAXvcPbz\n7IHuwJidgMVGhGc7x3opgFleil5FhcYTt+XK3AgK5X5VdjK/Qn+ZZjG6vNmBZeYdmHkLWTWPVFK2\nv4cpuZaXFuMZTU0pB5X4QU9g8zsApLyTZjN2Lh5+1TEm+IRAWiKcITxuOsAnYKXk/Gkh9XJYZ3xq\n3AA2B4jRPIc0vweIFrfN7VjxcE53jFsYmUCG/8TM64U0B7PrqZ5e1jvdfgH4vMFosfTM3jZlw+E4\nq/jt3brw4OkAUp6qeOXE2HOttkQ4Q7gISR2VtUS86i4eD2IFTNijmRFC6JWjzE6/CgplYn2uCnGJ\nfwFYjP60s0Q05bLoSTVT2OohjV1uMzyxDglW3pMLnCbeCRAUPFkhRwOEoCxB3Ioq9bGOMZ/u6ANw\nHEyI1sjBPfGXIHBhvhJh2ug28+oSPsfR56rpmwzFy4vfpzr710Iln3XGUAi9OtAWozdpvgxJ6Y9v\nJszyw3O8MHvvAAEvQfJuj2Ax9nQpFcJDcvtHiz61ZiYsi7NqjRQ4A0DSWdM1nASebKVnu0Yy24QU\n4O+9ibswP4lQKk/bHu/ofSE+RYxDAMzypMelXgzVy6dsOBzrA8ToWlqWdVa8udLl3t5Msjwmh/Tf\nAQ/sNQUh3r2W3alg9vCCAPsnpJKRuT82BGWaprdUxHUAQa4DlvJQxyMu/FxHH8woywF46qlsRKHa\nnYndgUsb+KyyMD+JsA+PSq+rSI7UI4RPObq79jtC6MVRPSm7r0KvPGS0/BxihXALxelmrSVyz6CM\nid+3XxSE2NwXvm+vf1KZCLD26XXlPThPIVtWn2jRqwNYYmdwSfOX+Iajz6kAUl552kq51jc7xtwd\nv0CKG/OZCC/Ep357DoB1zRuQ6lExvi+AmXmyCN4IEKN5gnbXK367CuAohDIHuSkyLFPwJqZ2lgSW\nfmZOWmarTOg7EUuayoOu/u8K4eCp2k6BslRr9jUXo5WxjtlbOmb2lOLhw3P7As/AF70Qit+e3OSW\nCDPh8+z6pT1f4uh+nmtMqSxD+tncviHo+cXD7D3KoNBTbI4WXZkx1ZNRi/HDHhsrG9XtU3hqdwwQ\n4+Q/K1fR+2PPbkNLK4ux95k5YwaJ1g++zs0kKfAjgBD0Fkff9wFI+ud0DSdBlsZigZcBqOP6DPag\nJcIspOWPuTy7K4FHOvo9Og3JUY6+tU/BHGU+of8Z/2bKVjOz8TiL8Sk17KyMPnEZfRWcHKxRsTGp\nhqEHFuPx9HOB61wr54NvbzDGXk2YfNHavnedHQNolvbEgffm9iWJjPwwu5fZliRty99O17RpzF8i\n7Jf2PNHR+0UwYS9jhkNS3lU9QaJ7ApjZSx19S8277CVu4YFcCRAtvtMxdhloXeY9v89iHE2aonTT\nBGKc2c+t0xvOh8V4DEWoC7BmNdg8B9F6y2oPoUD/sMujFFWumg6ZstVwFJ5k3nfOur2btWfF9QDg\nW+M+KIH5TIQJNwOPdfT7dPHbU2Xub44+qF/jo5vbNwTVUnQJCmU9kuM9BycACqFa6uDTFmOWqvJ8\ngsX4VvrBy5sUB0fZiBbLPVqCwjHZ/fshM4SQqeCSkAK/5fr+bEKxP52J8nr07A8+gPHqcPYw34nQ\nVeuhsnfhOdXaGcAsSQxl4tmpr3nynTcFiNH+4OgL/eI7roMTAIVwNf3A5BUWo+umMJdRHIw8s3i6\njkLw5GyXKOvyetMMS/Jba8pWQ2B9BaDsiAOznn6mx6tLqXy+/cHkEc4C5jsRpvoeZjtP024YPoCj\ntKHEn4uH2ZkFksp6JH/J7RuCyj53qXoKM+6v0FP2jha/NFXbqVAQQ5mpsGEjp8lzBAP/yxKFMGMB\n30FUlsS/CMr3KKv7wSG4KgZeAyDpQ46+n0198RywrEMRwJ0FszuRbvY/c4xZG/ObCNVL3/Fo6J0I\nYF3ziFVe6ehTIkvldwAbFr89y6TqZv1Do0W3ZqBCuGPgNFkW4ziVVBqFxbjOQIiMefcEAaLF3p5c\nUPDW0i4992xvsnJIki3Qa4aRFJqyvUHrWvm/elZqBwEXVPb+x4r5TYQJP8GhRKOOyiXPmY4xdwEw\nc+kNbpr6WvaBSwi9XGlitPUdY0O6YwNcEy3WirMbCCe5xWL8+6SN5ygKAuwXhneGyJQo3tMvFE9d\n11eMVmZ0VGvYZJkofntWSkcXvz0Fnt4LoE4vTz4HB+MMT2sCC4EIXwhAP7E8BxfgkO6v1G7IVmhR\n/5TTE4ID/c8sK9q/RFC4mX5GjsuzrKIgjrLGyAaFd+hJAxsrLMbtBpbC29YlwQLle7qXM2YQ+rU6\nvNfnA2CCDFwOzkp9XXvJe1EoPGUhaUQ+CPiyY8xGMP+JUPpa8ejoKdsNx3EA1jXPhZuW1tZbruZg\n29TXsoUYisplnwSI0bI3wgGCwnsoiLSyl+WGQvj5AIn8oSDEOSfMajEuKQiwF1dZLIWvqGu78l6e\nHxTy4+iAGHuiFK9zVKnDzMp96GyBB7NeiFT2QaB17V7FQ49azX2B3yKNraD7IGwWQnbSwGZSU6Ug\ny/Qqh72iTu631VFWHeE0bApylcgft5iz9z2IMfUPwf8eVknQKQSwEizGVVg542N1heAVf2gEFuNa\n9JWpSyxVCLW9YmjmvSwOwe4A/+da53tlxreAfV3f5679CtjBVV7V7HTgVqRGQrI83DL/PcKEpK7h\nk4b/CL78YejXtPVgOwAze4Wz/w7QJ0QPqhdstOipxLYSFMJthXe4VeXlWwsP0RNMXgsW4xsLD7BK\ngpsXXuCcIcEC5XzuNGWrSWBmpZx+dl56cUiyL76QMkjfx3yx2XTNHs54y5quPI0F4hEuBW4DjkT6\nZFbXrq0L3Ajsq46+kzcsa5NOgc+UOCGnb+pf2yu8lZR18vAQ9IXp2k9qp38hXxwU9pyycSaK0+Rh\nG/4nKATPQdVMxjyR4TV8V1UIWaKm06EpEozRnkgSHP5dCHLlhdf0Bk8CXgWs4sgmOZgkubVmthCr\n2e7Ax4HtG8sPd3DLwiDCZFDAbUjZWnHF8vi/6ihbxLTm8ng5SWjzHDnq6sIEj3CJZ0+pZ6d/Qf8l\nKLgFRKeCxfhJ4IhJ/vwl4GG5YSvFPuQFFEXIh+D9CsFTuW1aNEiCyyhqqdRYEn+PJOG/UUXkI6M/\nAn4tsWN2325BwP5l8e1Ip2T3ndTk4ibCU4BXOPcJjwbeA6ymTl46khm7kEJ4niSRHbxqffmoVaU8\nCXboyXNdBfX2C2E0e4aTwWJ8PvD6EZk/RiF4c3tnhCbfq8rNbKMQPCRma5Dksm6TxxFIAgtfADaV\n8oL9rWtrAv8EDlNHeUXK0rL498BhNaoJDjO7qImw3KR/HNJHs7unu9qZ6sixxK3lFfYOF2oskc8C\njgEuCkH39djo2RojGVZRLKGvpx/nOFNcB2ypEGrlY88URZxgb2+xARK8FlgOvDQEeVTMqzfTpRVF\n9Iz+iKQWmX1mYF37CPBYpze4F/AhYMcmhRYWNxEmo3VOjyOwn+cDNeOBwNeBwyU+M137lfvb+4En\nA/eW885Y8SoOCyHzzjxoa2JIzQZB4YZJGy8iRItb0c8quiwo3GuK5tPbi3YMRdxejSXxvqT83DdK\nPe3KjP697+62EtkhRIUDcY46jq0dszcD/2jqtLhvtiXCQ0k1hNfIlWeyrm1IUpY5WB1lR7jX8QpT\n/3oHJzCBDDet5Cb7bFl8PVBeWKcHhZPq2JvviBY/RTrdBDgiKHx6qvbT2ou2NSTiqbOlUfd7U2s1\n07UnkXL211YnU/jVbHVSPvQeSFfmjj216cUbPpMglUfw2Unf6vT2ZryCBPsBmFUKBuXh3qm/nevs\nD32Vkus9wgxVBIUX0FeaeXETgdfzEdGiFf97SYJrN0CCy6DnfWUf0JUws+8XD7MPOFJ/HlY89NZT\n+QBwUzYJJjwS+FHTJOjFwiLChB+Rai14cC8A61p2lTupJx/k2qQvlsQ3Aw8zc2W6lColpdxW7Ri5\noPDXgVhDRYtPrmt3viBaPJmB/cCg4JW8r6Lcz9wut0ZxCTPbnpTS9ntJ3mLo5wI3Svw+e/xur47K\n7s6xe9sCcwELkQjTXc5sn9yO6uiy4qG3hOXd09AuQUsklQcFXs1BQtCPKQp51wm2nmAzkeHjiqfv\nLwhxIX53AIgW1y+8wLKe732aOjiqfCaPDEF1avdeDiDJ5c2Z8eLi4TZTNpwcPwJQR/klIMy2JTkd\nnwcuT1MAAB5USURBVHOO3TgW3pdZ+lPx6CKnhYcBWNeW5w/dkz16bhGp78FeAGb2I2d/QtBZFNsD\nDZLhRwfI4PaFuFwu/qeeoEXhBX6vEdv9z+LFIchNAmZWamLec8qGk/ZnKUld5vzcUp0A1rVti4cP\n8oxPEv34cKU++axj4RFhQirMZJYbioE6KotSe+Oa1i5+u+TdJf2QROL3NrPJgoSnRQg6kbSH0xgZ\nQs873Kx8XniHtQ5m5gKK/6P6Pi1rMnyo8hm8JgR59DMBMLNHkN7/z0nyrlzK7/jDpmw1OX4HoI6+\nnt0zxTw+DU863gixsE6NJw4g4EtI2YVrrGuPAT4KbK6Ors0fmm+Qarve1bP/kmz0LpzVVePOGaN9\ngEKuvW7A9Uq2Ld4X+O7Ay0tqyE+NFYMxgQV2Cwpub3zoOH0SfG0I8mhYAmBmZUpnnVPiLYCrgedL\n+Vs41rUyhOgh6ihbpR2zY4FDkTx1lmc4xGIPn5k4wPuAp5AKfefLGaX4qJvUkafEZBPhND21lLrv\nU4xWDYWplYo31L7FR9EviFXivkHBuz0xUkSLT4CVsoDuGhRcN61Jx0ly+2Xe7ikh6FV17DVxc6z9\nvayXTmckKf7nVOTzGkcbPjMRxxW/T3P2fwSwjnXNu5n8IAAz3uTpXBTkPjDZsDoFhAhBL6CvE3dH\njOYi90ntK3ymWEauV3n5wnK5GS261FSaRLR438ryt0qCaxf7gE2T4Eb0SfAxDZLgfWuQYFmUylOe\noirF790bPID0nuQvqUeMhesRpkEuAvbxZJpAzbsf/bsvsI60khbeDG30lravlnSyx0aJGK26lL1f\nCLqwjr0px7L4DhiqyHM+8NBRL5+LZe+5DK/pe1pQ+N+RjR3tQfSLqt8rhF40ggtm9k7geODNkp7r\ns8E6pJIEX5JcdY5rXw+YfQ34ENL7Xf1nPEy7NB4cpNxTOR7pXdndu7YbcAlwgDr5rrwZvaWRdymS\n7PS8gYMkne+1AxCjbQCUweNnhaBj69ib0ZgWL2Bq1eO3Ac/2kmO0uA7wYWCqfafjg0L2dyB7LtHO\no19DZ50Q5K6EB2Bm5bbDLZK8ZUGrN+UlEp6tolJqayd19EvHBO5H8sR3wCEukjdUS4TDBroC2LqG\nV3gLsHoNr/BxpIv0dRIvmq795HZ6ZHg3yRG7VUFVCRmaP0SZcmyLq5FO1Ec95veB+wWFsVRFG8V7\namalslGtfWIzXkAKp9pV4sfTtR9qI3mDN6ojX9Ewsy8Dn0Z6t6t/1lAtEQ4baBOSqskRSNmpUZUc\n5NPUkWs5Vbkb31niT1M2ntTGhI33jSXVLq4eY0/DDhpYwrnnkYKzjwfe7jRx/6AweHo9NsRo9wHK\nbYZPhaBH17VpZptBr4b2EjkvVDNKO2dLPNFlo9s7bFtfHd3omMRewKeA7VCezJ0HLRFOPtjfgQ1q\neIWfBQ4jeYbZG9VFcPUdUHuJXGq/AawluQp/T0CMtjtwcfH0byFo47o2FxNitNvpHzruGII73a0H\nM1sPKBV/aoVPNXBKXCq4v10dPXO69pNM4lzgPCTvjS5zuJYIJxtsY+AvwDM9H4Z1e8ueP6ijraZr\nP3wK7AxcBtwmkS2e2bdjm5I0+ADWUKbKzmQYCLreJQT9tAm7CxUDXmBj2wvVWEFgA0luCTQzvkiq\nF+xfiXTtX8AaNQ5I7kMq8bl9riKUF234zGSQ/krKjXwbjgJP6kikjfgtrWv3802Bn5KSzFcxw3Xy\nl+zoeqAUZbilUCeujeJCLg8bLovRVFfBZiEiRlta3DRKEjywQRJclz4JLq9Jgo8hkeATapDgAcAa\nOGp/F5MwUu2Y08ZFgl4sDo8wDVjKmb8G6cXTNR9qomu3AUuBJQU5OqbR2y90CWH27fQl+qnpOQxi\nwDv8agjKrnO7EBGj/ZhCoQi4LgRtNlX7HFh/1QKwuZSf0dS31cse+YrEg102ur096b+ro42cEzkM\neBmwK9JYDq3SsK1HODmkW4DPAP9DqnrnQVkAu07pwfI9/30RXuOCpD/Qz/n9R0GMjaDwcMrC9QcU\n3uErp+qzkBGjvaO4OZQkuHbDJHg3+iS4UU0SXIVEgnhJsEB5+LSFcyLLgNOB/xknCXqxeDzCNGhZ\n9vObSMFlomsnAO8A9lBHl/im0TvJq3V4kmz1NrMB9pP0rana5yJG2xuoqq+8KQQ9r8kx5ipi7KVp\nltg5BP2syTHM7EBSFT5o4ACssuJYW8KlnWjd3v7n49XRR5wTOYFUsfDAJuuRzGzo9rBkJgM/nRSm\nsRXJq8o30e2X0KyxRC4rhzVBhiXBA5woKVuhezrEaA8Bvlx56Spgm6bzlmcbxb7o3+h7/wD7hpBX\n83omMLMOsKJ46g6R6dsjkpTS7yn5NDWt2/su/Vsdre6cyIbAL4CHIjUqYDGz4VsinOng7iJPANbt\nCSJcqo68Cr2YFUuHGns5E+31CPrHknadsrETMdquwKUDL+8Zgi4e1n6+IEYLwDcGXt4pBEcWxQxg\nZr8DtoX6ohrJHq8ATgEeKfkFT61rVwFbkk6KfQccKSXwv0jP9s6jDuYMEZrZCpLmWLnvcbKkLw+0\nmU0ivAdJBeNIpE+6THTtMOCzwFHq6BP+qfADYE/glRKneu307dlX6Bc7d5V3nAlinBDTWMXyEHT9\nKMZsGjHanYA/DvnTaiGMJvC3Wr4V+LSkyQreZ9ikrMvdkXiZ207XngK8DzhQHX11muaTTWYvkvL0\nTjR4gJc3hblDhB3gJkmTFvCeVSJME/gFSVp/Fe9mblkCFNhUHX/VuMq+zjGSr+bJRHtWpvUB7CKN\nNiYwRjsCGHZDGamwgwcx2kHAsCqFh4Qgb+GuGcHM9gB+WDx9pORXqe7bpLwhf0Ti8W47SZH9WuA8\ndXSwczJLgR8Ab0A62zuXuphrRHjzVHtVc4AIl5GK6FyMtKfLRD/Q2q/I0ZtOjwwPk6hVlzjZmxCO\n8WFJT6hrcyaI0Z4NvHmSP4+cbIbMZ5hWYomnh6B3jmMeZvYF+ko4jYQ7mVGq3HxdcktjNfc9Nnsm\n8GjggeM+IJk4jblFhE8lnWZeDLxw8IOfdSJMkzgK+BiwN9IPXCb6d9JvquM7ie5Pp0eGD5UmHEzU\nsJkyA4qn60j11FByEKNtB0wnEHEBqYaHSwygGMeAPUjBu/tN0fQ/wLYhaNhyeCQYSJe7RpKrQuHK\ndikl1X4msXMtW127nKRRuLE6zhx2s61I1/p+SL+oM5+6GCsRFntRw2KpTiWFW5TeyMtJAaLHDPSf\nfSJME+mdAHvvYtbtLQ1PUEdn1ptOjwwfIdWKV6zYtEeSYigBXiap04TdXMRoawFngn8Jl4l3Ac8N\nYXayGszsNdBTHDpEasYbrpDgNRK1iNW6dhLwKurtCy4BvgJcgHR6nfk0gTnjEU4YwGxr4FxJOw+8\nLqBbeSlKiiOdzDD0czvPRTrUbaZrZ5Mu8HtVyoI6p9QjwydLfLCOrb7NCZv0AJsopR7OOmK0pcCD\nSbVuD5+meRUi5bG+l5QBM5KDoVyYrXQI09ihlRllGNNVElvXstW1PUl7em9RR8+pMamnA08G7o90\n23TNm4aZBVKNoBKdOUGEZra5pD8Xj58P7CnpcQNt5oZHCGB2PPBO4L5I7joblfjC9dTR/9WbUo8M\nT5GoJfM+0a4dCXy8ePojSbs1ZXuxw1Ju7eXAdsVLh0o6tzn7lEXFfiilsq9uW13bCPgrcLM6yq72\nWJnUtiQyvT/Sr6ZrPg7MpRS7083sMjP7CWnP5vnTdZhVSGeSvMILi0MUL1Ypft9Y5GrWmFIvyPqV\nZisVGqphV58gfe5/AnY1M5nZrMR7LSSY2UmkA4ftgMslWcMkeDKJBM9pgARXIZEgTKwzkzuppaRw\nm1fNFRL0YnEGVA9D/xT5BqQNpms+qZlu/7S27klymhZXk/I9b5TwqQNPanuC+CfAA2dle2Iew6wn\nYV+iEdHciWPwOVIxsdMlTqptr79yWVcd3TRl46kndhqwP3DAXMonnkse4fxDqqNwf2D9Ik/SZ6aj\nvwL3ALCu1c62KDbD3w2sZ4YKkddGIOna4guzf/HSNwoP8T5NjbFQYWah2OcuSfB+hRfYGAmaYWbc\nSiLBxzdEguUp/g41SXBf4JnA4+cSCXrREmEV0ndJp6vvoIaaizpKeZawe3GIUnNaHAc8tnh6hxnL\n69qcaF/fKAixPM29sCDERzY5zkKAmR1VEGCZjndEQYCNBo6bsT5pqb0asJuET/ygarNrXyYt3R+o\nji6vMbmNSAH7RyO5tA7nGtql8TD0Q2rcWScA1rVjSGKsb1KnvmKLGdsAZf3dI6Wh2Ry1MXCgAvDe\nwfCnxYTiEOQTJDWVEo0ehEwcj2o50E0l3FlLPZtdew9wNHCEOvm1eyqTWwJ8Hvg10ol15zUKzMnw\nmUkHnttEWNaA/T8k/2YyYF0rK4h11dGK+lOjqjRzsYQrK2ZmY01ICSuxvaTfjmrMuQQz2550ClzF\nzlKzUlwTx+RM4DjSDW87T+nNlWx2rRRkeLo6NTNpko7AAcD+4yjE5EG7R9gUpJuAfYB1MasVIKqO\nXk9S6e1Y12rfQSVuL06UvwvsUewbrl3X7vCxdHHxhVqTflzcb4pl86VmttYoxp1NmNm6ZnZ5sfwt\nSfAXwGrFEngkJGjGWkXI1HHAiRJ3bYgEOyQSfGkDJPgIUqznEXOVBN2QNCs/aejZGXvGP7BCICXB\n03r/7wpezwrECl7c3HuoB6fpSaBnjOlzewgpkLn682tg01n/vPz/03KSBzb4f+0/nvH18MrnuGNj\ndlfwiuI79/La9mBHwfWCvWf785rB56nsPvNpsrPyA1cU39Atav/PK3hN8cV8ZXPvo1apXEQCLR3j\nZ3jCEPIQcBJJaHT2P7/h814KvHSSuT9xfPPQUtDPi8/t56DG3jNWcEbxXTuttj1YT/BrwdGz/dnN\n8PNVdp/5NNlZ+QGrsMzqtf/vFZxafEHPbPb91GkVMnzSLHye95uEWEQSVthptj5D4J5AnGJ+e45/\nThO8+UMbtb2CjxTfsRfVtvf/7Z1/sFxlecc/X/MDDDi0GgQMkQCCCRVFJ1FsIIZiMECAhAKGgmQU\nkLEGBKH80IyHA4opEIqUHy0/BIdCGAdEQzFAgmNBOpaGUskkpBpiaiKQdqrFglakPv3jOZu7u/dm\n7727Z/ec3X0+M2ey5z3nvO+b3b3Pvj+e5/vAeIPVBjcU9fk18XnbaJ+JzZKRIO0M/CY7a1qcYXt1\nqT4L3Ag8YYk1UksZXb2+VljtG7aLGS0ngR99PyTgdBg2Tvop4F7gfrPWxFwl7QXMx2NePzTM7aea\n2X2ttNcs2We0DV93fQmYYkZu621K9SxwCPApS+zO1iqTgLvw6JM/pUv8BWPXuJ1Ik4CtwDas9Qxm\nSnUCruT7e2Bss7lPhqxbXAhcm51+x4z5edXdLFn+5XNxoY3mcmGMnleBJcCt5lkMC0UiYSBHyRFm\nfD+3ugfSbwIcZYmtar1SXQnMwXeIO/6D2ixhCNuNVJE/Wok1qeJbXV2q9wEVHb4JluT3x5qlCq3+\nBV9oVuMbWBrk4Y0zgMOBPwYOBvYd5rFNwLPAM7hz8z9bSUcsEocDleyC9wCf8FlnTvWnNZkMW1Y/\n8kp1NnAJ8GGsefX1IghD2AkGnI1vxeyclqsbEHYF2N8S29To/lHXL/YDXqgqOtiMtvnBBQNITAF+\nWlW0uxm5Sp8p1XuASiqGllJGDFSqE4GbgFmYDSesWzrCj7ATuHrLYuDTmXNpa9Ultg0PowJ4QalO\nbrXOmvqNTeZ+h5V612a+h7koJQeDkXi7xC8ZMILTfT0+dyN4GgNGcGxORvAYPG/3Md1oBJslRoTN\n4usnS4BzMbsxlypTrQaOBO62xM7Io85BbYjLgWqF6nebDYqeCJpAYm/c+bqi75dL/pkh20r1AHAi\n8JAlzQsK11aqI3Gpr+Mx+2EudRZATI07jXQL7kt3PmZfy6XKdHsCeoBdLGnPIrXEV6FGzeRjZjzW\njrZ6HYn6XM8nme0wYVRrbaWq9gw40xJrOeuhV6yZeDa8kzB7Yrjby0xMjTuN2WfwtZTrkf5iuNtH\nVGVitwBTs9PXlOqwPOod1I5xWTZlXpwVPZpNmf+6He31GplE1gVZWFzFCB6VTYHbZQSPYMAI7puj\nETwMN4KndbsRbJYYEeaBdC1wIbAEs6/kUmWqMbif2e7AvZZYWxMeScwEflBX/C6zmo2WvkdiL1x9\npSJ2sQWYZcbmtrXp6TYfxqXdtgL7WJJTfhbpo7gv52lYDi43JSCmxkUiXQVchsuWfyG3alOdB1Sm\n3VMssX/Pq+4h2xO74K4pb68qfgY43IzCffGKQGI88GUGMtKBO8R/3qwmIVb+baeqll47yxK7I7/K\ndRxwB+4s/WRu9RZMGMKikZbg6Uu/hdlosrE1rrbWxeZqS+ySvOpu2K6YA4PWDR8HjjXjt53oQ1Fk\nPwhfAi6uKt6Ep1ntiPuRUl0DVBSL3maJ/SK/ynUqcD0wD7N6qbWuJgxhGZBOB+4GNmG2f65VpzV5\ncvexxH6WZ/0N2xYVXcVqHsKz7PWEX6LEVGApLo1f4RXcGf2RjvUj1YG4og/AMktyFED1sLmL8Cif\nYzFbO8wTXUcYwrIgzQL+ITsbQ045bQGU6p1AZXr8AHBynuF5I+qDOAzfJHpv3aXbgavN6Ar/M4lp\nePTEorpL/wqcZcYzHe2PrwWuAOZlRXtZYi83eGSUDWgMcAMewXMMZltzq7tEhCEsE6r5Vf9DzP47\n1+pTXQxURGOPtcS+m2f9I+6H2Bn4DD6SGl93eTOe7vGeojddstHeqbiw6KS6y6/g0+BbzfjfTvcN\nQKlOhO27zZ+zxG7ItwFNwDdFdsXXBF8Z5omuJQxh2dD2JNoAMzBrOatdTfWpJuDRC5WNjbZvpgyH\nxE64o+9iPG54KF7BY7b/CfghnrD8ly20KeBtuOrKB/BshEfADpW7vw/8LfBAuzc7hkOp9oHtO86b\ngWmWWL7G2AVDHgSeB86m19Sl6whDWEY82c2rwJuBxZjdlHsTqQ7BBQjAIxum5yng0CqZAMR03EDO\nBd7XoaafBv4eX8t8zozclihaRanegu/GH5AVvccSW5d/QzoMj42/EVjaqoRcNxCGsMxIXwc+CTyC\n2dFtaSLdvlEDcD+w0JJyKrJUk43oJuFGYXJ2TMR18MbjitJjAQGv4Xp+lWMrsBF4sUyGbkco1Thg\nJR5KCXCGJXZ3g0eabEjCo55SYBFmK3Nvo6SEISw70kI8lhNgN8x+lXsTvuB+De7gDS6O+sncHHCD\nplCqscDfAR/PipZYko/z/eDGtDM+AjwUmE+fZB2sEIawG1DNmtBczB5tSzMu1HkLnhUNXAR2oSXW\n0/5/ZUOp3oxHoszJim4GFrdtp1+aCtyHZ+D7FGavtqWdEhOGsFvwact6PKb4fsxyld6qacoN4o34\nzi74TvYsS1qTxg8ao1TvAJ4E9suKrgEuaaMBFPAobnDPAW7rh/XAoQhD2G1In8Z3LwH2xuznjW5v\nublUFwDXVRXNs8Qebmeb/YZSHQ1UuzKdb0k+ykQ7blSTgYpz/Z9htrzR7b1OGMJuxJMOvZidJZhd\n0fYmU/0JHipXYTlwtiX2Wrvb7kWUajfcX3JBVfFMS+wf29+4UtwHEmAiZv/V9jZLThjCbmZA2xDg\nrZg17Vc34iY918W9wLFVxecAt3U6WqXbyDalKjJsFR4CTrck/02wwR3YnkwM4G8ySbiAMITdj/Ru\nYEN2lmJ2eceaTnUMLvVU4ffAQuD+MIpOZvzq05S+ARxpSQd1/KRlwOezsylYsU70ZSMMYa8w4HMI\nMA2zDY1uz7Vp31w5n8ECC5cC1/fbrrNS7YILXSR1l84DbuqoW5I0Dd9kA1iG5SjG0EOEIewlatcO\nnwRm5yneMKIuuFE8E7i17tKPcOOwutdGi9mobw7wVTxcr5pzgNs77pPp6U7XMCByMQmzFxs80deE\nIexFpLOA27Kz8zArTEo/C+VbCnys7tI6fAS5PPc42TaT+fmdik81/6ju8iPARW0JfRsp0udw3UCA\n0zG7p7C+dAlhCHsVHxGsBmZlJbOKVhTORk7H4VPmDw9xyzo83O++ooUgYHt/p+Jag6cA7x/itieA\nrwCrCh/peka51dnZI7h2YEQHjYAwhL2OapSqAd6J2ZaiulOPUk3Gtf0+ARy4g9sMn+qvwUUH/gV4\nwRJrSQVGqXYG9senjzNw43xog0c241P+uyyxl1ppO1c8MuT57Ox3wGTMthXYo64jDGG/IH0Ql7AC\nFyE4sKxrRlkSqunAUfiUemaHu/AUHnHxHWBt4SO9HeFrwhuBCVnJwZj1hPJ3pwlD2G9IJ+AxxOAa\nf1OxHBWNO4RS7YSrz0wG9s6OiXii9J2y403Ab4FfVR2/wPOIvABs6QalnUFIU/C44HFZybFYMSK7\nvUIYwn5FNerGvwbei1mk4SwztT6jAMdj9lBR3eklIsF7v2L2LfyDPwmfWm1EsixnbVAmpOOQjAEj\n+FHMFEawWMIQ9hJmD2QGsZJ8fFVmEK/MlLKDIpDGIH0tM4ArstIZmQF8vNGjQWeIqXEvI+0JfA+Y\nlpW8DBzRyUiVvsYjQR7F1z7BMxvOzzuRV1BLTI2DWsxexuwg/HO+CNgTeD4bJS7LlIyDPJF2Qro5\nG/2tx43gF4E3YTY7jGA5iRFhvyHthzvoHlBVeiFwfTjsNok0Fg85vKqqdA1wYpn8PPuF2DUORoc0\nF1+zGldVeinwV72e8rFl3PhdAFxdVfo6cFJsfBRLGMKgOVzmfREuLlrNN4ELMds6+KE+xF1ergXm\nVZW+AZyM2beHfijoNLFGGDSHmWF2V7aLKeBwPFb4FGBLtqZoSFcg/UGxne0g0h5IV23//7vLyzxc\nfecj2fs1Loxg9xMjwqAx0kR8/eviIa5WEoc/1fWJgnxUPB04F4+VrmcpcBVm/9PRfgWjJqbGQftx\nifgLcIMxfog7HsKTy68o7Q6ptCse93wGcPwQd/wG1yO8sRMpE4J8CUMYdB4fSc0EPotL+++Ibbjq\nzA+y4zmsNcWZBn3aFRdVnZEdHwKm7ODu1/G10TuBp7t+ZBuEIQxKhDQBmAvMx3UU92myplcB4aPP\nccPc24gteE6W7wLfwyJjX68ShjDoHtyZewruz3gA8K7s2BPYreqo53VchebXeBa3LXhO359lr/8N\nWIfZG+39DwRlJQxhEAR9T7jPBEEQNEEYwiAI+p4whEEQ9D1hCIMg6HvCEAZB0PeEIQyCoO8JQxgE\nQd/TtCGUdLKkdZL+T9IH6q5dJuknkjZIOqr1bgZBELSPVkaEa4EFwBPVhZIOAj4OHISHWN2sSBw0\nYiTNLroPZSTel6GJ9yUfmjZQZrbBzH48xKUTgOVm9jsz2wxsBD7YbDt9yOyiO1BSZhfdgZIyu+gO\n9ALtGKm9A48BrbAVmNSGdoIgCHJhbKOLklbhQfD1fMFGl5chpI2CICgtDQ2hmc1pos6fM5DHFWDv\nrGwQcvnzoA5JSdF9KCPxvgxNvC+t09AQjoJqpYcVwL2SrsOnxAcAT9c/EMozQRCUhVbcZxZI2gIc\nCjwsaSWAma3Hs5+tB1YCf26h+hsEQYkpTI8wCIKgLHTcvy8csYdH0uWStkp6NjvmFt2nopA0N/s+\n/ETSJUX3pyxI2izpuez7MWjpqV+Q9HVJ2yStrSp7q6RVkn4s6TGNIAVtEY7O4Yg9PAZcZ2bvz45H\niu5QEUgag6cLnYt/L06VNK3YXpUGA2Zn349+9tO9E/9+VHMpsMrMDgQez84b0nFDE47YIyY2k/zz\n32hmm80z3t2Hf08Cp++/I2b2JFCfcvV44BvZ62/gCcQaUqYRVzhi13KupB9JumMkQ/seZRKekKlC\nv38nqjFgtaQ1ks4uujMlYw8z25a93gbsMdwDebnP1BCO2MPT4D36InALcEV2fiWwDDizQ10rEz37\n+efATDN7SdLuwCpJG7LRUVCFmdlI/JXbYgjb7YjdC4z0PZJ0OzCaH49eov47MZnaWUPfYmYvZf/+\np6QH8WWEMITONkl7mtnLkvYC/mO4B4qeGtc7Yi+UNF7SvuzAEbsfyD68CgvwDaZ+ZA1wgKQpksbj\nm2krCu5T4UiaIOkt2etdgKPo3+/IUKwAFmWvFwHfHu6BtowIGyFpAXADMBF3xH7WzI42s/WSKo7Y\nb9Dfjth/KekQfGr4U+CcgvtTCGb2hqTFwKPAGOAOM3u+4G6VgT2AByWB/w3fY2aPFdulYpC0HPgI\nMDEL8PgSsBT4pqQzgc3AKcPW07+2JgiCwCl6ahwEQVA4YQiDIOh7whAGQdD3hCEMgqDvCUMYBEHf\nE4YwCIK+JwxhEAR9TxjCIAj6nv8HtcdGszlcFVwAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Noutputs = 1000\n", "xs = np.zeros((Nbodies, Noutputs))\n", "ys = np.zeros((Nbodies, Noutputs))\n", "times = np.linspace(0.,50*2.*np.pi, Noutputs, endpoint=False)\n", "for i, time in enumerate(times):\n", " sim.integrate(time)\n", " xs[:,i] = [sim.particles[j].x for j in range(Nbodies)]\n", " ys[:,i] = [sim.particles[j].y for j in range(Nbodies)]\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "fig,ax = plt.subplots(figsize=(15,5))\n", "for i in range(Nbodies):\n", " plt.plot(xs[i,:], ys[i,:])\n", "ax.set_aspect('equal')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At this stage, we might be interested in particles that remained within some semimajor axis range, particles that were in resonance with a particular planet, etc. Let's imagine a simple (albeit arbitrary) case where we only want to keep particles that had $x > 0$ at the end of the preliminary integration. Let's first print out the particle ID and x position." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ID\tx\n", "0\t6.82422856846e-05\n", "1\t0.991001690909\n", "2\t-0.938416458567\n", "3\t-2.16909675239\n", "4\t-0.0073018056096\n", "5\t-4.9263375846\n", "6\t-4.91601116735\n", "7\t-2.1198298552\n", "8\t1.96480987036\n", "9\t5.29695349399\n" ] } ], "source": [ "print(\"ID\\tx\")\n", "for i in range(Nbodies):\n", " print(\"{0}\\t{1}\".format(i, xs[i,-1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let's use the `remove()` function to filter out particle. As an argument, we pass the corresponding index in the particles array." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of particles after cut = 4\n", "IDs of remaining particles = [0, 1, 8, 9]\n" ] } ], "source": [ "for i in reversed(range(1,Nbodies)):\n", " if xs[i,-1] < 0:\n", " sim.remove(i)\n", "print(\"Number of particles after cut = {0}\".format(sim.N))\n", "print(\"IDs of remaining particles = {0}\".format([p.id for p in sim.particles]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By default, the `remove()` function removes the `i`-th particle from the `particles` array, and shifts all particles with higher indices down by 1. This ensures that the original order in the `particles` array is preserved (e.g., to help with output).\n", "\n", "By running through the planets in reverse order above, we are guaranteed that when a particle with index `i` gets removed, the particle replacing it doesn't need to also be removed (we already checked it).\n", "\n", "If you have many particles and many removals (or you don't care about the ordering), you can save the reshuffling of all particles with higher indices with the flag `keepSorted=0`:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of particles after cut = 3\n", "IDs of remaining particles = [0, 1, 9]\n" ] } ], "source": [ "sim.remove(2, keepSorted=0)\n", "print(\"Number of particles after cut = {0}\".format(sim.N))\n", "print(\"IDs of remaining particles = {0}\".format([p.id for p in sim.particles]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the `particles` array is no longer sorted by ID. Note that the default `keepSorted=1` only *keeps* things sorted (i.e., if they were sorted by ID to start with). If you custom-assign IDs out of order as you add particles, the default will simply preserve the original order.\n", "\n", "You might also have been surprised that the above `sim.remove(2, keepSorted=0)` succeeded, since there was no `id=2` left in the simulation. That's because `remove()` takes the index in the particles array, so we removed the 3rd particle (with `id=4`). If you'd like to remove a particle by `id`, use the `id` keyword, e.g." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of particles after cut = 2\n", "IDs of remaining particles = [0, 1]\n" ] } ], "source": [ "sim.remove(id=9)\n", "print(\"Number of particles after cut = {0}\".format(sim.N))\n", "print(\"IDs of remaining particles = {0}\".format([p.id for p in sim.particles]))" ] } ], "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 }