{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Primordial Earth\n", "There are a wide variety of problems in the conext of the Solar System requiring accurate integration of N-bodies undergoing close encounters and/or collisions. Standard integrators such as IAS15 and WHFast might be insufficient for these types of problems.\n", "\n", "In this example we investigate the primordial Earth embedded in a disk of planetesimals, integrating it for a short period of time using the MERCURIUS integrator. MERCURIUS is a hybrid integration scheme which combines the WHFAST and IAS15 algorithms in a similar way to the hybrid integrastor in the MERCURY package." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import rebound\n", "import numpy as np\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First let's choose the basic properties required for the MERCURIUS integrator. In particular, we are: \n", "* setting planetesimals to *semi-active* mode, which means they can influence active bodies but not other semi-active bodies.\n", "* merge any planetesimals that collide with a planets, conserving momentum and mass.\n", "* remove particles from the similation which leave our pre-defined box.\n", "* track the energy lost due to ejections or collisions. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "sim = rebound.Simulation()\n", "\n", "#integrator options\n", "sim.integrator = \"mercurius\"\n", "sim.dt = 0.025*2.*np.pi # we're working in units where 1 year = 2*pi\n", "sim.testparticle_type = 1\n", "sim.ri_ias15.min_dt = 1e-6 # ensure that close encounters do not stall the integration \n", "\n", "#collision and boundary options\n", "sim.collision = \"direct\"\n", "sim.collision_resolve = \"merge\"\n", "sim.collision_resolve_keep_sorted = 1\n", "sim.track_energy_offset = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that the setup is complete, it's time to add some particles! When using the MERCURIUS integrator it is important to add active bodies first and semi-active bodies later. The `N_active` variable separates massive bodies from semi-active/test bodies. Here, we add two active particles, the Sun and the Earth. Thus, `N_active` will be 2." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "sim.add(m=1.)\n", "sim.add(m=3e-6,r=5e-5,a=1,e=0.05,f=np.pi)\n", "sim.N_active = sim.N # sim.N= 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's create our planetesimal disk. First we define three different distribution functions - powerlaw, uniform and rayleigh." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def rand_powerlaw(slope, min_v, max_v):\n", " y = np.random.uniform()\n", " pow_max = pow(max_v, slope+1.)\n", " pow_min = pow(min_v, slope+1.)\n", " return pow((pow_max-pow_min)*y + pow_min, 1./(slope+1.))\n", "\n", "def rand_uniform(minimum, maximum):\n", " return np.random.uniform()*(maximum-minimum)+minimum\n", "\n", "def rand_rayleigh(sigma):\n", " return sigma*np.sqrt(-2*np.log(np.random.uniform()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let's set up the basic properties of our planetesimal disk. For this simple example we are assuming that all planetesimals have the same mass and radius." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "N_pl = 8500 # Number of planetesimals\n", "Mtot_disk = 10*sim.particles[1].m # Total mass of planetesimal disk\n", "m_pl = Mtot_disk / float(N_pl) # Mass of each planetesimal\n", "r_pl = 2e-5 # Radius of each planetesimal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's add our planetesimals to the simulation!" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "np.random.seed(42) # by setting a seed we will reproduce the same simulation every time\n", "while sim.N < (N_pl + sim.N_active):\n", " a = rand_powerlaw(0, 0.99, 1.01)\n", " e = rand_rayleigh(0.01)\n", " inc = rand_rayleigh(0.005)\n", " f = rand_uniform(-np.pi,np.pi) \n", " p = rebound.Particle(simulation=sim,primary=sim.particles[0],m=m_pl, r=r_pl, a=a, e=e, inc=inc, Omega=0, omega=0, f=f)\n", " # Only add planetesimal if it's far away from the planet\n", " d = np.linalg.norm(np.array(p.xyz)-np.array(sim.particles[1].xyz))\n", " if d>0.01: \n", " sim.add(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We move to the COM frame to avoid having the simulation drive away from the origin. In addition, it is always good practice to monitor the change in energy over the course of a simulation, which requires us to calculate it before and after the simulation." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "sim.move_to_com()\n", "E0 = sim.energy()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's quickly plot the location of the particles with matplotlib:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+UlNWd5/H3lx+NQEwAQaKtiCGsiQaGzvSIDjk5biLRyAodo6LBHZxVmZxdZ45xwgYCJ6gDirKrZmeyO0GTGbJyFHW10q5GxB+cnHGFsU2jLRoCqFFKBSKiBhCa5rt/1NNOdVPd9TxVT/14qj6vc/p01VP3qbpP/6hv3Xu/915zd0RERLoNqHQFRESkuigwiIhIDwoMIiLSgwKDiIj0oMAgIiI9KDCIiEgPCgwiItKDAoOIiPSgwCAiIj0MqnQFCjF69GgfP358pashIpIoL7zwwh/cfUy+cokMDOPHj6etra3S1RARSRQz+32YcupKEhGRHhQYRESkBwUGERHpQYFBRER6UGAQEZEeEpmVJFItUu1pVqzdwtt7DzBoAHQe6fn4yGGDWXLhGbQ0NVamgiIFsDh2cDOznwP/Adjl7l/K8bgBPwYuAPYDV7r7b4LH5gKLg6JL3X1Vvtdrbm52patKuSxOdXDPhjfL+pqDB8CKS6YooEiszOwFd2/OWy6mwPBV4I/AL/oIDBcAf00mMEwFfuzuU81sFNAGNAMOvAD8qbu/39/rKTBInFLtaRY+9BIHen/crzIKFlKssIEhlq4kd/+1mY3vp8gsMkHDgQ1mNsLMTgDOAda5+x4AM1sHnA/cG0e9RHKZc9dzPLt9T6WrEVnnEbhuzSauW7Ppk2PTJoxi9TVnV7BWUovKNcbQCLyVdX9HcKyv40cxs3nAPIBx48aVppZSc5LSGijUs9v3MH7Bo5/cN+CO2WpVSHESM/js7iuBlZDpSqpwdaSK1Xow6I+TaVXc8tgrbFw0vdLVkYQqV2BIAydn3T8pOJYm052UfXx9meokNWb67evZumtfpatRFXZ+dOiTlsTYYxsUJCSScgWGVuBaM7uPzODzB+7+jpmtBW42s5FBuW8AC8tUJ0mw7jTR9N4Dla5K1csOEkMGDeDWb09WV5P0K5bAYGb3kvnkP9rMdgBLgMEA7v6PwGNkMpK2kUlX/cvgsT1m9nfA88FT3dQ9EC3Sl0qkj9aKg4ePcN2aTXz/gRf5b5f8iQKE5BRLumq5KV21PqXa03xvzSaq6S+2ccRQ5p93WqQ32GrMilJ2U30o6zyGclNgqC+p9nSPFM1ymXj8cNZdf07ZXi/VnuaG1s3sPdBZttfMprGI2qfAIImWak/zt/dvoqtMf55XnDWOpS2TyvNiEaXa01y/ZhPlyLFScKhtCgySWOVqIVRzMOhLqj3NjY9s5v39pW9V3Kn5EDVHgUESqZT977XYj16O7KwkBlDJTYFBEucLix7j45j7jurpTS3VnmbRwx3sO9RVkudXCyL5FBik6pXy0+7QwQO45aL6ztcvRetLS24kmwKDVLVSvGnVU+sgilLM+9AgdTIpMEjVmrzkcT48GF93h7o4wilFV5N+9smiwCBVJ+5sI70pFS7O30UtDurXqrLuxyDSn7i7MtSNUbzugBpHC+LZ7Xv4/MJH2XbLjDiqJlVALQYpqTiDQiHLT0g4cY35qBVX3dSVJFUhexOZQunNpjziWotKXUvVK2xgGFCOykj9mX77+qKDwshhgxUUyqilqZE7Zk9h5LDBRT3Ps9v3MHnJ4zHVSipBLQaJXbET1ZR2Wh3i6AZUYK8uGnyWsoprZdA3lmsAs1p0B+digsN1azbR9vs9CvQJE0tXkpmdb2ZbzGybmS3I8fgdZrYp+Pqdme3Neqwr67HWOOoj5ZVqTzP/gReLCgpjj21QUKhCS1smcefsKTQMtIKf454NbzL99vXxVUpKruiuJDMbCPwOmA7sILMb2+Xu/kof5f8aaHL3/xTc/6O7fyrKa6orqbo03fREUat9qusoGYodnC73/hZytHIOPp8JbHP319z9EHAfMKuf8pcD98bwulIFFqc6Cg4KAyzTB62gkAwtTY28vnwGE48fXtD5W3ft44wfPU6qPR1zzSRucQSGRuCtrPs7gmNHMbNTgFOBp7MOH2NmbWa2wcxaYqiPlEkxg5MTjx/Oa7fM0MBkAq27/hyuOGtcQefuO9TFdWs2KThUuXKnq14GPOju2VMtTwmaNt8B7jSzCblONLN5QQBp2717dznqKv1ItacLDgpXnDVOXQoJt7RlEm8sn8ExBY49/O395d+qVcKLIzCkgZOz7p8UHMvlMnp1I7l7Ovj+GrAeaMp1oruvdPdmd28eM2ZMsXWWIsy567mC19lR11Ft+e2yCxhUQGzocjTXoYrFERieByaa2alm1kDmzf+o7CIz+wIwEngu69hIMxsS3B4NTANyDlpLdZi6bF3BSydccdY4dR3VoG23zGDssQ2Rz/vwYBdfWPRYCWokxSo6MLj7YeBaYC3wKnC/u282s5vMbGZW0cuA+7xnGtQXgTYzexF4BljeVzaTVN7iVAc7PzpU0LlqKdS2jYumFzTu8HGXK5W1Cmnms4RWyBIXWjenvhS6nLdWzC0PrZUksZq6bF3kc95YPkNBoc60NDXyRgEprTs/OqQxhyqiwCD9mnPXc4xf8GjkLqRC+pyldhSSdfbhwS7m3PVc/oJScgoM0qdC1+g3ULeAMG3CqMjnxL0PuBRGgUH6VMg/6aeHDOR1rXkkwOprzi5olvT4BY+yONVRghpJWAoMklMhTfqxxzbw0o3nl6A2klSFzpK+Z8ObCg4VpMAgR0m1pyO3FiYeP1zdR5LT0pZJBXUrxblPuESjwCA9FJJuqCUuJJ/V15xdUHDQmkqVocAgnyg0B10T1ySM1decHXnbUC24VxkKDAJkgsLCh6L36Ra6yqbUpyUXnhH5nELX5ZLCKTAIAD986CUOdHblL5hFG+xIVC1Njdw5ewpR190rZIKlFE6BQVic6mB/55FI57yxfIaCghSke8OfKHZ+dEhrKpWRAoNEzv64c/aUEtVE6smwwdHefrbu2qfxhjJRYKhzUT+FjT22QUtnSyxuvmhy5HM03lAeCgx1buuufZHKa66CxKWlqbGgFFatp1R6Cgx1LGprIWLLXySv1decHTmz7dntezQrusT0r16nUu3pyK2FFZdobEHi171/dBSaFV1asQQGMzvfzLaY2TYzW5Dj8SvNbLeZbQq+rs56bK6ZbQ2+5sZRH8kv6mbsd86eorEFKakhg6K9HWkgunSKDgxmNhD4CfBN4HTgcjM7PUfRNe4+Jfi6Ozh3FLAEmAqcCSwxs5HF1kn6l2pP0xVh4z4FBSmHW78dbTB6/gMaiC6VOFoMZwLb3P01dz8E3AfMCnnuecA6d9/j7u8D6wAtz1lCUZe9UFCQcmlpaoy0ZEbnEQ1El0ocgaEReCvr/o7gWG/fNrOXzOxBMzs54rmY2TwzazOztt27d8dQ7fpTyFpICgpSTlGXzHh2+x51KZVAuQafHwHGu/tkMq2CVVGfwN1XunuzuzePGTMm9grWulR7mvkPvBjpnEI2WREpRktTIyOGRltob9HDylCKWxyBIQ2cnHX/pODYJ9z9PXc/GNy9G/jTsOdKPG5o3UznkfADC4OssH17RYp1w8xorYZ9h6Kt8SX5xREYngcmmtmpZtYAXAa0ZhcwsxOy7s4EXg1urwW+YWYjg0HnbwTHJGZ7D3RGKr/tFm3PKZXRvdBeFJrXEK+iA4O7HwauJfOG/ipwv7tvNrObzGxmUOxvzGyzmb0I/A1wZXDuHuDvyASX54GbgmMSo6h9sFGb8iJxa2lqjDShUvMa4mXuEfIWq0Rzc7O3tbVVuhqJ0XTTE7y/P3yLQZlIUg2iJktoGfj8zOwFd2/OV04zn2tcqj0dKShMPH64goJUhahrKd2z4U1lKMVEgaHG3dC6OVJ5DThLNVl9zdmRsuNWrN1SwtrUDwWGGpZqT0cadI66q5ZIOUT5sJLee6B0FakjCgw1qpB5C3O0f7PUAGUoFU+BoUatWLsl0ryFgYYG7qRqRVmaWxlKxVNgqFFvR2xS//dLtaS2VC99aCkvBYYaNTTirjrKRJJaEnUTKulJgaEGpdrT7O88Erp844ihJayNSDyGRfiwE3UTKulJgaEGRU3Zm3/eaSWqiUh8br4o2n4NGoQunAJDDYqSsjdtwih1I0kiRN2v4d6Nb+UvJDkpMNS51decXekqiIQWZb+GrgQu91MtFBhEJDHUui0PBYYaE2WtGA06SxJF+bvV1p+FUWCoMTc+En5tJA06SxJF+bvV1p+FUWCoMWFXUtWgsyRV1L9bLawXXSyBwczON7MtZrbNzBbkePx6M3vFzF4ys6fM7JSsx7rMbFPw1dr7XAkvyicjDTpLkkXpTtLCetEVHRjMbCDwE+CbwOnA5WZ2eq9i7UCzu08GHgRuy3rsgLtPCb5mIgUL+8koyhr3ItVo/nmnaTXgEoqjxXAmsM3dX3P3Q8B9wKzsAu7+jLvvD+5uAE6K4XWll7DrI6m1IEnX0tTIHRH3hZbw4ggMjUD2TJIdwbG+XAX8Kuv+MWbWZmYbzKwlhvrUpVR7GmVtSz3RGFnplHXw2cyuAJqBFVmHTwn2IP0OcKeZTejj3HlBAGnbvXt3GWqbHIXsvSBST7Q8RjRxBIY0cHLW/ZOCYz2Y2bnAImCmux/sPu7u6eD7a8B6oCnXi7j7SndvdvfmMWPGxFDt2hF17wWRWjHQwo00rNYeDZHEERieByaa2alm1gBcBvTILjKzJuCnZILCrqzjI81sSHB7NDANeCWGOtWVqHsviNSKy6eenL8QqJs1oqIDg7sfBq4F1gKvAve7+2Yzu8nMurOMVgCfAh7olZb6RaDNzF4EngGWu7sCQ0TDGgaGLjs8QlmRaqcNfEpjUBxP4u6PAY/1OvajrNvn9nHe/wP0my3SvkNdocsu+5Z+3FJbGkcM1VyFmGnmc51RJofUmrBLZGhpjPAUGEQk0cJ+2LmhNfw6YvVOgaGOaDVVqWd7D4RbR0wUGOqKVlOVeqfupHAUGOqIxhekVg0dHO6tTCuthqPAICKJ93HnkVDllL0UjgJDwoVtGg8ZpF+11K4TNX4WK71bJFzYpvHBw+E+UYkkkcbP4qXAkHBhm8Zh15QRSSKNn8VLgaFOdLlWixGRcBQY6sTIYYMrXQURSQgFhjqhBoOIhKXAUCc+0KxPEQlJgaFOKJ1Pal3IOW6a/RyCAkOdUDqf1LpPHRNuHE2zn/OLJTCY2flmtsXMtpnZghyPDzGzNcHjG81sfNZjC4PjW8zsvDjqI0dTOp/Uur37w3WXasfD/IreqMfMBgI/AaYDO4Dnzay1105sVwHvu/vnzewy4FZgtpmdTmYr0DOAE4EnzezfuXv4nWfq2POtP+VfGm7jRPsDb/tobjt8Ka1HvlLpatW9VHuaFWu38PbeA5w4YijzzztNgbkMhjUMDLVp1Qhl6OUVR4vhTGCbu7/m7oeA+4BZvcrMAlYFtx8Evm5mFhy/z90PuvvrwLbg+SSP51t/ypdeWMxJA/7AAIOTBvyB5YPvZuaAf6l01epaqj3Nwoc6SO89gJOZgLjwoQ71a5fB/pA7GSpDL784AkMj8FbW/R3BsZxlgj2iPwCOC3mu5HDyb1Yw1A71ODbMDvFfB91foRoJZPqvD3T2fIM60Nmlfu0yCPt+rwy9/BIz+Gxm88yszczadu/eXenqVNzxnvtncKK9V+aaSLa++q/Vr109PjNUXUn5xBEY0sDJWfdPCo7lLGNmg4DPAO+FPBcAd1/p7s3u3jxmzJgYqp1suyz3z+BtP67MNZFsfaUFK124emjZsPziCAzPAxPN7FQzayAzmNzaq0wrMDe4fTHwtLt7cPyyIGvpVGAi8K8x1KnmvfXl+ez3hh7H9nsDtx2+tEI1EsikBQ8dPLDHsaGDBypduIqEzV6qZ0UHhmDM4FpgLfAqcL+7bzazm8xsZlDsZ8BxZrYNuB5YEJy7GbgfeAV4HPgvykgK589m/hULOq9mx5HRHHFjx5HRLOi8WllJFdbS1MgtF02iccRQjMw+27dcNElZSVVErbf8ik5XBXD3x4DHeh37Udbtj4FL+jh3GbAsjnrUm0f8K7QeUiCoNi1NjQoEVUytt/wSM/gsR1PanUh0Ctr5KTAkWKOaxCKf0JhyfBQYEixKk1gTrKTWqQEdHwWGBIvSJL6hdXMJayIitUSBIeHCdift1WxPEQlJgSHhlGEhoq7SuCkwJFzY7iTt+Sy1TGtRxUuBIeHCflJacuEZJa6JSOWEXYtKy2GEo8CQcAsfeqnSVRCpuLCzmTX3JxwFhoQ70HkkVLn5D2wqcU1EKmdYQ7i3Ms39CUeBoU6EjB8iibR1175Q5ZSsEY4CQ8INUJ+pSGhaDiMcBYaE+87UcaHLKqVPRMJQYEi4pS2TQpf9oQaqRSQEBYY6sl8DDVKDFqc6Kl2FmqPAICKJdu/GtypdhZpTVGAws1Fmts7MtgbfR+YoM8XMnjOzzWb2kpnNznrsn83sdTPbFHxNKaY+9SrKpB2NM0it6Qo5OUGpquEV22JYADzl7hOBp4L7ve0H/sLdzwDOB+40sxFZj8939ynBl5LtCzAnwgC0lg6QWhKlG0mpquEVGxhmAauC26uAlt4F3P137r41uP02sAsYU+TrSpalLZNC/yLTIZcOEEmC1RvfDF1WqarhFRsYxrr7O8Htd4Gx/RU2szOBBmB71uFlQRfTHWY2pMj61K3bZ6sXTuqPlrgojbyBwcyeNLOXc3zNyi7n7k4/myiZ2QnA/wb+0t2702MWAl8A/gwYBfygn/PnmVmbmbXt3r07/5XVmSifhubc9VwJayIiSTcoXwF3P7evx8xsp5md4O7vBG/8u/oo92ngUWCRu2/Ieu7u1sZBM/sn4Pv91GMlsBKgublZnxOK8Oz2PZWugkjRoiRSaNn5aIrtSmoF5ga35wK/7F3AzBqAh4FfuPuDvR47IfhuZMYnXi6yPiJSJ6IkUmjZ+WiKDQzLgelmthU4N7iPmTWb2d1BmUuBrwJX5khLXW1mHUAHMBpYWmR96po+FUk9iZJIoYHnaPJ2JfXH3d8Dvp7jeBtwdXD7HuCePs7/WjGvLz0tufAMrlujjF+pfZrtXFqa+VxDWpoamTZhVKiyGoCWpEq1p1m9IXyaqia2RafAUGNWX3N2qHIagJakWrF2S9/pjzloYlt0Cgx1TMtjSBJFGVsY3jBQ4wsFUGCoYxqPkFq37Fvhl6WXf6PAUIOGDdavVWQAykYqlN5BatDNF00OXXbqsnUlrIlI5WgWbOEUGGpQS1MjV5wVbsXVnR8d0liDJEaUNNUTlY1UMAWGGhVly88f/B9t+SnVL9We5p6QaaoDTNlIxVBgEA4ePqIJQ1LVUu3pSMkSt186ReMLRVBgEIBIE4ZEyinVnmbhQ+E/uEybMEpBoUgKDDUsSnaSBuqkWq1Yu4UDnV2hy4ed5Cl9U2CoYVGykwCm376+NBURKcLb2nWw7BQYalhLUyMjhoZfcXXrrn0lrI1IYT4T4W84bDae9E+BocbdMDPaOvSa1yDV5qOPO0OXjZKNJ31TYKhxLU2NmIUvr3kNUk1S7Wm6Qg6AqbUQn6ICg5mNMrN1ZrY1+D6yj3JdWZv0tGYdP9XMNprZNjNbE+z2JjGbMzXaP0yUnbFESilKiqpaC/EptsWwAHjK3ScCTwX3czng7lOCr5lZx28F7nD3zwPvA1cVWR/JYWnLpND7NEBm9Uq1GqTSJi95PHRZ7bkQr2IDwyxgVXB7FZl9m0MJ9nn+GtC9D3Sk8yWa1deczdhjwzfIrr9/k4KDVMziVAcfHgyfoqpZzvEqNjCMdfd3gtvvAmP7KHeMmbWZ2QYz637zPw7Y6+6Hg/s7AM1KKaGNi6aHLnvE4YbWzSWsjUhui1MdoZe+gMzYgia0xSvvns9m9iTw2RwPLcq+4+5uZn0NE53i7mkz+xzwtJl1AB9EqaiZzQPmAYwbp0GmQg0eAJ1HwpXdeyB8NohIHKIGBdDYQinkbTG4+7nu/qUcX78EdprZCQDB9119PEc6+P4asB5oAt4DRphZd3A6Ceiz78LdV7p7s7s3jxkzJsIlSrYVl0yJVF6T3qSc7t34VqWrIBTfldQKzA1uzwV+2buAmY00syHB7dHANOAVd3fgGeDi/s6XeLU0NRIhe5Wtu/ZpgT0pmy6PtjjLnbOjfdCRcIoNDMuB6Wa2FTg3uI+ZNZvZ3UGZLwJtZvYimUCw3N1fCR77AXC9mW0jM+bwsyLrIyHMiZjvfc+GNzUQLSUX9QOIxhZKxzxihK4Gzc3N3tbWVulqJNr029dHWgJj0ABj280XlLBGUs/m3PUcz27fE7q8Aa8vn1G6CtUoM3vB3ZvzldPM5zq17vpzGBChT+nwkeR9gJBkSLWnIwUFgDvUhVRSCgx17PZLo/1zjV/wqLqUJHZRdxAce2yDupBKTIGhjrU0NTI0wp4NoIlvEq9Ue5qDh0PmTweizMeRwigw1LlbIu7ZcMRh0cPKUpLiRd2uE7RQXrkoMNS5lqbGyP9s+w51KYVViva9iEFh2oRRmsxWJgoMwtKWSZHzwZXCKsWYc9dzkbeT1Zad5aPAIEBhLYfr1mi8QaKLmpoKRFodWIqnwCCfKKSZHrWPWOrb4lRH5KBwzEBTa6HMFBikh0IG99RqkLCiLpA38fjh/HaZJlaWW97VVaW+LG2ZFPmft7vVoNxyySXVnmbF2i2k9x6IdJ6RmYgp5acWgxylkFbDdWs2KVNJjpJqT7PwoY7IQQE0u7mSFBjkKEtbJhUUHJSpJL2tWLuFA53hd2LrpgXyKkuBQXIqJIUVlKkkPRXSUrjirHGar1BhCgzSp0JSWEHBQTJdSOMXPBr5PAWF6qDAIP0qtFvpujWbtPtbnSpkqQvQzOZqosAgeS1tmcTE44dHPm/rrn1MXvJ4CWok1ayQoDDINLO5mhQVGMxslJmtM7OtwfeROcr8ezPblPX1sZm1BI/9s5m9nvWY0hCq1Lrrz2HssQ2Rz/vwYJdaDnVk6rJ1kc8ZZLDtFm26U02KbTEsAJ5y94nAU8H9Htz9GXef4u5TgK8B+4EnsorM737c3TWNtoptXDS9oKUJtu7ap70catziVAfjFzzKzo8ORTrvmIGmoFCFig0Ms4BVwe1VQEue8hcDv3L3/UW+rlTI6mvOZuSwwQWdq0Hp2pNqT3Pa4l9FnhTZTbOaq1OxgWGsu78T3H4XGJun/GXAvb2OLTOzl8zsDjMb0teJZjbPzNrMrG337t1FVFmKteTCMxgYZV/QLFF365Lq1T3IHHWjnW7aW6F6mXv/i9+a2ZPAZ3M8tAhY5e4jssq+7+5HjTMEj50AvASc6O6dWcfeBRqAlcB2d78pX6Wbm5u9ra0tXzEpoVR7mkUPd7DvUPTJS6C0xFrw+R8+VvBe4NMmjNJgcwWY2Qvu3pyvXN61ktz93H5eZKeZneDu7wRv8rv6eapLgYe7g0Lw3N2tjYNm9k/A9/PVR6pDS1MjLU2NBacm3rPhTdZtflfbNCbQ4lRHwV1HAHfOnqJZzVWu2K6kVmBucHsu8Mt+yl5Or26kIJhgZkZmfOLlIusjZdbS1FjQDGmAnR8d0qB0gnRPWis0KAwyeGP5DAWFBMjbldTvyWbHAfcD44DfA5e6+x4zawa+6+5XB+XGA88CJ7v7kazznwbGkFlIcVNwzh/zva66kqrT1GXrImeldFPXQnUrtpXw6SEDeenG82OskRQibFdSUYGhUhQYqlfTTU/w/v7O/AVzGHtsg7qWqkyxAQE0nlRNwgYGzXyWWC258AwKTFj6pGtpyo1PqHupCkxdtq7ooHDn7CkKCgmkwCCxamlq5PZLi5vAvvdAp/Z3qLDJSx4vuFuw27QJozSekFAKDBK7lqZG3lg+o6AlNLJpf4fy657B/OHBwtKQu11x1jiNGSWYxhikpApNZ81FfdWlM3nJ40UHA9Agc7XTGINUhWLSWXu7Z8ObTFio9NY4daegxhEUrjhrnIJCjVCLQcoi1Z5m/gOb6Cxs9YSj6JNp4VLtaW5o3czeA4Vlj/WmbLLkULqqVKVUe5oVa7cUtOVjLiOGDuaGmWdokDOkOXc9x7Pb98T2fOreSxYFBql6xUyI602fWvuWak9z4yObC55fkkvDQOO2i/9EATlhFBgkEabfvp6tu/bF+pyaRf1v4hpUzqZWQnIpMEhixDG7NpfhDQNZ9q1JdfepNtWe5ntrNhH3f/bE44ez7vpzYn5WKScFBkmcUgWIbrXc3VTKn526jWqHAoMkVtwDpLkkvbupHD8j0BLZtUaBQRKv1C2IbEnoJinFeExfkvDzkOgUGKRmlPMNsdvgAbDiksp8Wo57nkFYwwYP4OaLJquFUMMUGKTmVCJA5GIGc6YWn5mzONXBvRvfoqvC/4NqHdSPsgQGM7sEuAH4InCmu+d8tzaz84EfAwOBu919eXD8VOA+4DjgBeA/unvexHYFhvpWLQEi6ZI+ziLRlWutpJeBi4Bf91ORgcBPgG8CpwOXm9npwcO3Ane4++eB94GriqyP1IF115/DG8tnMG3CqEpXJZEmHj+cN5bPUFCQPg0q5mR3fxUgs2Vzn84Etrn7a0HZ+4BZZvYq8DXgO0G5VWRaH/+rmDpJ/ch+YyvnQHVSaWKahFVUYAipEXgr6/4OYCqZ7qO97n4467hGvaQgS1smffKmF+dSG0lWrxP8pHh5A4OZPQl8NsdDi9z9l/FXqc96zAPmAYwbN65cLysJ1HsSW7ly/quBWgUSh7yBwd3PLfI10sDJWfdPCo69B4wws0FBq6H7eF/1WAmshMzgc5F1kjrS3eVUiwFCKaZSCuXoSnoemBhkIKWBy4DvuLub2TPAxWQyk+YCZWuBSP3JNdiatGCh5SmkHIoKDGb2LeDvgTHAo2a2yd3PM7MTyaSlXuDuh83sWmAtmXTVn7v75uApfgDcZ2ZLgXbgZ8XURySqXMEi7j0jCqV0UqkUTXATiaA7aLy99wAjhg3GHfYe6MQg9GqmBszRWIBUQNh5DOVGU52qAAAEa0lEQVToShKpGS1NjerGkZpX7AQ3ERGpMQoMIiLSgwKDiIj0oMAgIiI9KDCIiEgPiUxXNbPdwO9L+BKjgT+U8PnLQddQeUmvP+gaqkVc13CKu4/JVyiRgaHUzKwtTK5vNdM1VF7S6w+6hmpR7mtQV5KIiPSgwCAiIj0oMOS2stIViIGuofKSXn/QNVSLsl6DxhhERKQHtRhERKQHBQbAzC4xs81mdsTM+hz5N7M3zKzDzDaZWVUt7xrhGs43sy1mts3MFpSzjvmY2SgzW2dmW4PvI/so1xX8DjaZWWu565mjPv3+TM1siJmtCR7faGbjy1/L/oW4hivNbHfWz/3qStSzL2b2czPbZWYv9/G4mdn/CK7vJTP7crnrmE+IazjHzD7I+h38qGSVcfe6/wK+CJwGrAea+yn3BjC60vUt9BrI7IexHfgc0AC8CJxe6bpn1e82YEFwewFwax/l/ljpukb5mQL/GfjH4PZlwJpK17uAa7gS+IdK17Wfa/gq8GXg5T4evwD4FZlVz88CNla6zgVcwznA/y1HXdRiANz9VXffUul6FCPkNZwJbHP319z9EJmd82aVvnahzQJWBbdXAS0VrEtYYX6m2df1IPB1M7My1jGfav+7yMvdfw30txXfLOAXnrGBzLbCJ5SnduGEuIayUWCIxoEnzOwFM5tX6coUoBF4K+v+juBYtRjr7u8Et98FxvZR7hgzazOzDWZW6eAR5mf6SRnP7G/+AXBcWWoXTti/i28H3TAPmtnJOR6vZtX+tx/W2Wb2opn9yszOKNWL1M1GPWb2JPDZHA8tcvewe01/xd3TZnY8sM7MfhtE+bKI6Roqqr9ryL7j7m5mfaXMnRL8Hj4HPG1mHe6+Pe66Sg+PAPe6+0Ez+ysyLaCvVbhO9eY3ZP72/2hmFwApYGIpXqhuAoO7nxvDc6SD77vM7GEyTfCyBYYYriENZH/SOyk4Vjb9XYOZ7TSzE9z9naCZv6uP5+j+PbxmZuuBJjJ95JUQ5mfaXWaHmQ0CPgO8V57qhZL3Gtw9u753kxkPSpKK/+0Xy90/zLr9mJn9TzMb7e6xrwOlrqSQzGy4mR3bfRv4BpAze6CKPQ9MNLNTzayBzEBoxbN6srQCc4Pbc4GjWkFmNtLMhgS3RwPTgFfKVsOjhfmZZl/XxcDTHowmVom819CrP34m8GoZ6xeHVuAvguyks4APsrotE8HMPts9NmVmZ5J5/y7NB4xKj8RXwxfwLTJ9jgeBncDa4PiJwGPB7c+RydZ4EdhMpvum4nWPcg3B/QuA35H5hF1t13Ac8BSwFXgSGBUcbwbuDm7/OdAR/B46gKuqoN5H/UyBm4CZwe1jgAeAbcC/Ap+rdJ0LuIZbgr/7F4FngC9Uus696n8v8A7QGfwfXAV8F/hu8LgBPwmur4N+sg+r+BquzfodbAD+vFR10cxnERHpQV1JIiLSgwKDiIj0oMAgIiI9KDCIiEgPCgwiItKDAoOIiPSgwCAiIj0oMIiISA//H04T+bw4nzFeAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "coords = np.zeros((2,sim.N))\n", "for i in range(sim.N):\n", " coords[0][i], coords[1][i] = sim.particles[i].x, sim.particles[i].y\n", "fig, ax = plt.subplots()\n", "ax.axis('equal')\n", "ax.scatter(coords[0],coords[1])\n", "ax.scatter(sim.particles[1].x,sim.particles[1].y); # Planet" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, let us simulate our system for 100 years, and check that our final relative energy error is small. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "times = np.linspace(0.,1000.,1000)\n", "encounterN = np.zeros(len(times))\n", "totalN = np.zeros(len(times))\n", "errors = np.zeros(len(times))\n", "for i,t in enumerate(times):\n", " sim.integrate(t,exact_finish_time=0)\n", " totalN[i] = sim.N\n", " encounterN[i] = sim.ri_mercurius._encounterN\n", " errors[i] = abs((sim.energy() - E0)/E0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The default values in this notebook yields an error of $\\approx 10^{-8}$. The following plot shows the relative energy error, the number of particles within 3 Hill radii of the planet and the number of ejected/merged particles as a function of time. " ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAGtCAYAAABEC0OXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8nXX5//HXlT3aJN27TWlLS2kpo5RVtmwQFAVFFETh\ny/eLCCpqAQEHP0FREARURKaCDFlSdtmjlLa0dO+06UzSZs8zPr8/zsnpSZpxknNOTsb7+XjkkXOP\nc9/Xfe6Tc658pjnnEBEREZHuLSnRAYiIiIhI+5S0iYiIiPQAStpEREREegAlbSIiIiI9gJI2ERER\nkR5ASZuIiIhID6CkTURERKQHUNImIiIi0gMoaRMRERHpAVISHUA8DB482OXn5yc6DBEREZF2LVq0\nqMQ5N6S9/Xpl0pafn8/ChQsTHYaIiIhIu8xscyT7qXpUREREpAdQ0iYiIiLSA7SZtJlZspn9qKuC\niZaZnWNmD5SXlyc6FBEREYmTmgYvzjkA5m/cTU2DNybHrff68Pj8ANQ2+GJ23FixxotudQezBc65\nWV0UT0zMnDnTqU2biIhIz/fgBxu5de4qRuVlcsT4gWwtrWVBwR6mj8qlpKqeHeV1zBiTx1nTh7Ow\noJRvzBrDox9v5oj9BnLi5KHc8tIKRudlkp6azIQh2SQnGRuLq7nw8DE8sWALKUmG3znmrSpiR3kd\nB47MYezALD7ZuJvstBQ+mnNS3K/RzBY552a2u18ESdtdQCrwFFDduN45tzjaIONFSZuIiEjvkD9n\nbkLP/8oPj2XqyJy4niPSpC2S3qMHB3//OmydA+KfeoqIiEif1j8jhco6L/sNyWZjcTXjB2czYUg2\nq3ZUUlnn4cQpQ9lRVse0UbnsrKglLyuNJIPKOi+pyUlUBZ+7pLCMJDOG9k+nos6D1+/ITE1mZF4m\nH60vYdKw/uQPyuKFJdso3FPLMRMHUdPgw99O4VZXajdpc86d2BWBiIiIiDSXk5HKKVOH8fvzDwIg\nJTm+fSh/curkuB4/Gu1euZnlmtmdZrYw+PNHM8vtiuBERESkbyuv9ZCXmUZKclLcE7buLpKrfwio\nBC4I/lQAD8czKBGRvqy8xhN6vHxbOT5/oHrGG+zVFgsNXn+ol1xbCvfU0F7bZ+l6q3dWhO7f6p0V\n1Ht9CY4oPjw+P1X1XvKyUhMdSrcQSZu2Cc6588OWf2VmS+IVkIhIb+Lx+an1+FixrYINxVV4fH4+\nXFfCyLxMxg3K4vvH7sdry3dQ7/VT7/Ezf9Nunlu8jVnjB+Lx+fl8SxmH5w/g8PyB3P/uBm49bxpn\nTh/BppJqCkqqWV9chdfn57OCUhq8fjLTklm5vYLHvzcLM6iq9zF6QCZ/f38jJ0weSkqSMW/1LuZ+\nsYMpw3N4/PuzeHrhVvx+xxnTh9Pg9bNocykj8zJZtLmU219dzZXHT+CYiYP45Usr+NWXpzF70mC8\nPj8pyUmU1TTgHKwrquJv723goNF5eHx+fnzK/jhg8+5qxgzMYlNJNWt3VeL1Of67dDtDczIYnpNB\ng89HWY2HnMxUpo7IoazWQ029l28fNY6qei8fr9/NYeMGUF7r4a/vbeCakycxaVj/0Ou7obiKcQOz\nSE4yPisoZU91A68s28FBo3MpqWrg8U8KeOenJ1Be42H59nKOnjCY0poGdpTX4fc73l5dxLJt5YzM\nzeScGSPZXV2P1+e4bPb40DkWbd7D1tJaNhZXc9SEQTR4/azeWcFZB40kNcmo9/rJy0rlkY8KGJ6b\nQWFpLRW1gcT70017OG7/wQzPyWBdURU/P30KGalJvLO6iIlD+5GdnsLSwjKy01M4eEwej8/fjN/v\nuPjIcRRX1lMePM728joWby7lmImD2VBcxe2vruaSo8bx1UNHc+59H3Hm9OH83wkTKa6sZ9ygLHx+\nh9/Bu2uK+MassSws2MONzy/nlWuOZWB2GnuqGxiQlYrfwduri0hLSWLLnhrSk5NIS0ni3INHsru6\ngQavn7IaD1tLa3hvbTGjB2RRWefh+P2HMCI3kx3ltXy0voRTDxxOcpLxyYbdHLf/YNJTkvE7R5IZ\nJVX1DMxOo6iynvxB2fTPSCEjNZnS6gaeWLCFYyYOJiXJeG9tMclJRr/0FE6cMpQFm3azeXcNgJK2\noEh6j34C/NQ592Fw+RjgD865o7ogvg4xs3OAcyZOnHj5unXrEh2OiAhTb36NmobWS0Gu/dIk/vRW\nz/q8mji0H+uLqpgxOpelW7t+XMz9BmdTWFrDQaPzWLS5FID+6SlU1sduTK2vHDKKXRV1rNheEUqc\neoNB2WlYMJEalJ1GeW2gQX5zA7JSKa2J33U3Jo6RuvsbB3PuwaPiFk+ixXLIjxnAY0BjO7ZS4BLn\n3BdRRxknGvJDRLqLRA9XIJ134MgcVmyvSHQYArx6zbEcMCK+w24kUkyG/DCzJGCyc26GmeUAOOf0\nDhYRiYAvrATjrgtn0D89lQUFe5g4pB8Ox5ThOYwZmEXB7mr8fkdZjYfkpEApyModFWSlJVNcWU9B\nSQ3nzBhBSnIS8zfuZvqoXCrqvGwJVj1W1Ho4cGQu64oqyUhN5rQDh/P5llI+31LGIWPzyMlMZcbo\nPLx+x5srdzEiN4O0lCTqvT4yU1N4Y+VOJgzpx66KOraW1jJhSDbby+uo9/jIyUxlUHYah4wdQG2D\nj21ltRw6dgD1Xh8lVfVMH5UXrPqqwznwucAwCptKqhmQlUZWWjKfbyllWE4GGanJfLCuhOQkmDw8\nh3mrdrH/sP5MGNqPl5duJznJGJCdxqDsNOq9fpLMmDAkmy17apg+Kpf+GakU7qnhtGnD8fkdzyws\n5JiJg0lNTmJTSTWlNQ1kp6ewsbiKd1YX8e2j8qlp8PLMwq1ccPgYMlOTeWd1ETPG5DK4XzrriqqY\nMKQfQ/uns6O8lv2G9CPJjEHZaWSnpzCkfzrFlfW88Pk2TjpgKLmZqdz/zgZmTxqE3w8fbShhwpB+\nlNd62FPdwOyJg8lMSyY3M5WUJGN7eR1JBp9t2sOBo3IZkZvBlj01jB+cTV5WGlV1XkbkZZCalMSC\ngj3sNzibMQOzWL2zAq/Psa2sllF5mWwqCdznjNQkHvxgExcePoYBWak8/FEBJ04eitcfqEafMKQf\nRZV1LCks50sHDGVwv3SG5WTQLyOFgdmBYTBWbK9g6sgcfD7H3GU7OHXqMDLSkrn15ZUcOnYAB4zI\n4eMNu/H5/fgdJBlsK9v7Xli1o4Izpg1n9qTBzP1iJ1NG9GfcoCx+/NRSTjtwONnpyWSmJTMiN4Pt\nZXU458hKS2HGmFycgzW7Kjl6wuBAdfamPQzPDbwvXl+xk7EDs8jJTOWZhYXkZaWSPyib5dvKqar3\n9eqErSMiKWlbGEn2152opE1EuoOKOg8H/fINbjzzAC4/br9EhyMi3VSkJW2R9B59y8yuM7MxZjaw\n8ScGMYqI9GqVdYE2Vv0zIunzJSLStkg+SS4M/r4qbJ0D9G+jiEgbqoJJWz8lbSISA5G0abvYOfdR\nF8UjItJrVNYFet/1z9BwBSISvTarR51zfuDeLopFRKTHqfP4ePqzQmb/7m2q671sLa3h8fmbcc6F\nhqBQ9aiIxEIknyTzzOx84DmnYbFFRELqPD6m3PRaaPnAW14PPZ6VPzA0wGr/dCVtIhK9SD5J/gf4\nMeAzs1rAAOecU/9bEenTVmxvfWDZ0/70fuhxbqaqR0Ukeu0mbc65/u3t012EzYiQ6FBEpA9YUhjZ\nbABD+qfHORIR6QvaHfLDAi42s5uCy2PMbFb8Q+s459x/nXNX5Obmtr+ziEiUdlXUtbvP8JwMzKwL\nohGR3i6S6tH7AT9wEvAboAq4Dzg8jnGJiHR79Z7AnKKPfPdwhvbPIC8rlffXFjMsJ4OtZbWBEe4H\nZCU4ShHpLSJJ2o5wzh1qZp8DOOdKzSwtznGJiHR79V4/w3LSOWHy0NC6b8wam8CIRKQ3i2RGBI+Z\nJRMYUBczG0Kg5E1EpE+r9/pJT0lOdBgi0kdEkrTdAzwPDDWz/wd8CPw2rlGJiPQA9V4faSmRfIyK\niEQvkt6j/zKzRcDJBIb7OM85tyrukYmIdHMNXj/pStpEpItENOKjc241sDrOsYiI9Cj1StpEpAvp\n00ZEpJPqPWrTJiJdR0mbiPR4dR4fm3dXx+XYbc3eV+/1kZ6qj1ER6RqRDK57tZkN6IpgRKTv2F5W\nS1Erg9MuKSxjZ3kdG4qrWtzu9fl5d00R1z/3BS8u2cbFD37K8Xe8y/Jt5VTUefbZv8Hr5/21xby9\nehdLCsto8PrZUFxF4Z6aVuP7eEMJD36wkfHXv8JLS7fvs/2DdcUs3Vqu6lER6TKRtGkbBnxmZouB\nh4DXNXG8iETr6NvfBqDg9rNC6656YjHzVu2iztN0VKGTpgwlOcn4/fkHkZqSxLSwidmfXFAYenz2\nnz8EIC8rlU/mnExmWjJ3vrGG5z7fxtbS2hbjCD//sq3lrNlVyXXPLG2yz/3vrGdzSTUfri/hutMm\n8/W/fhLalpKkpE1EukYkvUd/EZzC6lTgu8C9ZvY08A/n3IZ4Bygivd/jnxRQVFnP3C92tLj97dVF\nAByy8k1G5Ga0e7yyGg+H3fomyWZU1nvb3PfpzwrJyUzhucXbeGPlrhb3Wb2zktU7KwGaJGwQqJoV\nEekKkfYedWa2E9gJeIEBwLNm9qZz7mfxDFBEeh+Pb29J2o+fWsJzn2+L+Lk7ytuf7xOgpiGyZOpn\n//ki4nO3pKqdpFBEJFbaTdrM7BrgO0AJ8CDwU+ecx8ySgHVAXJM2MxtLYIDfPcBa59zt8TyfiMRf\ncWV96HFHErbuKNLkUEQkWpE0xhgIfNU5d5pz7hnnnAfAOecHzm7riWb2kJkVmdnyZutPN7M1Zrbe\nzOa0c/7pwLPOucuAQyKIV0S6uZ2tdEBIhP0GZ3PT2VM57cBhnXp+dYNK2kSka0RSPXo3gJkNDFtX\n6ZzzRDAzwiPAvcBjjSuC85jeB5wCbCXQyeElIBm4rdnzLwPmE6iKvQx4PIJ4RaSbK6tpaLI8K38g\nCwr2dMm5j9pvEMu3lXP+YaP53uzxpCYnMTw3g8uOyae81sMbK3Zx1IRB3PTichZvLqWiru2krFrV\noyLSRSJJ2hYDY4BSAtNY5QE7zWwXcLlzblFrT3TOvW9m+c1WzwLWO+c2ApjZv4FznXO30ULJnZld\nB9wSPNazwMMRxCwi3Vhjk7ZjJw3m8PyB/PDkSTy9sJDPt5Tx5IIt/O786Vx4+FjW7qpkzIAsfv3y\nCs47eBQpycaf3lrH4fkDqff62FFeR1Wdl2E5GTz1WSENPj/TRuWQnpLMos2lAAzMTmNkXgbLt1UA\n8OQVR7YYk5mRl5XGBYePAeCR787C4/Pjd44Fm/Zw0wvLKdhdw01nT2VEbgZrdlZy97x1fOuIcfF/\nwUREAGtv9A4z+zuB6snXg8unAucTSJ7uds4d0c7z84GXnXPTgstfA053zn0/uPxt4Ajn3A9aef40\n4JcE2tRVOeeua2W/K4ArAMaOHXvY5s2b27wuEUmc15bv4Mp/LuaVHx7L1JE5TbbtrqpnUL/0Dh+z\nqt7Lmp0VHDZub6VAbYOPlGQjNTmJt1fvYmj/DKaNyu1UzFX1XirrPIzIzQytc85hZp06nohIIzNb\n5Jyb2d5+kZS0Hemcu7xxwTn3hpn9wTn3P2bW8U/WDnLOLQe+FsF+DwAPAMycOVPjyIl0Y40lbclJ\n+yY8nUnYAPqlpzRJ2AAy0/ZOMXXSlM61WQs/fr/0ph+ZSthEpCtF0hFhh5n93MzGBX9+BuwKtk3z\nt/fkFmwjUN3aaHRwXdTM7Bwze6C8vDwWhxOROPEFS/hbStpERKRlkSRtFxFIrF4AnieQcF1EoOPA\nBZ0452fAJDMbb2ZpwDeAlzpxnH045/7rnLsiN7dz1R8i0jV8/sD/e0raREQi12b1aLA0bY5z7upW\ndlnfzvOfBE4ABpvZVgIdCv5hZj8AXieQ+D3knFvR4chFpMcKVY+qelFEJGJtJm3OOZ+Zze7swZ1z\n32xl/SvAK509bmvM7BzgnIkTJ8b60CISQ35/sHo0WUmbiEikIqke/dzMXjKzb5vZVxt/4h5ZJySi\nevTmF5fz/triLjufSG/gbUzaVNImIhKxSJK2DGA3cBJwTvCnzZkQ+pLHPtnMdx5akOgwRHqUxo4I\nSZF8AomICBDBkB/Oue92RSA9UYO3M51nRaSxejRFWZuISMTa/cQ0s/3NbF7j/KFmdpCZ/SL+oXVc\nVw/5UauJokU6RdWjIiIdF8m/uX8HrgcaJ4r/gsAwHd1OV7dpq/EE5hxMS1ZpgUhHNJa0qaBNRCRy\nkcyIkOWcW9Bs5O8+P0PyhuIq/vbeBgDSU/XNI9IRjW3aVD0qIhK5SD4xS8xsAuAgNHfojrhG1QOc\nd+9HPL1wKwCVdV7y58ylss4T2l7n8eHx9a42b5tKqqn3qkpYoudTSZuISIdFUtJ2FYE5PaeY2TZg\nE3BxXKPqpK4cp62yft/Cxt1VDWwqqeamF1ewtLCMqSNyeOWaY+MeS1eoqPNw4h/e5auHjuLOCw5O\ndDhxsaSwjBmjczWfZBz9ed46RuRlhpI2tWkTEYlcu//nOuc2Oue+BAwBpjjnZjvnCuIeWSckehor\nB/zqvytZWlgGwModFV1y3jqPj7W7KgFYWLCHU+58L+adJGrqA8d7bvE2Zv/u7ZgeG6Cooo79f/Eq\nX2wti/mxI/HhuhLOu+8jHv6oICHn74yKOg8rtwfeY0UVdbzweUym8O2UNTsrufONNbhgtWdr/vjm\nWq57ZunejgiaxkpEJGKR9B5NN7OLgGuAH5nZzWZ2c/xD65na+9KKh+ufW8apd71PWU0Dv3l5JeuK\nqlgTTOJi4Y7XV3PkbfNCy1tLa2N27EYfri+hwevnoQ83tbi9zuMjf85cnlu8NebnBthaWgPA6p2B\nJGjBpj0Jqd6+9t+f8+V7P4xo30sfWsCZ93yAc47vPvIZ1z61hPIaT4v7+vyOxVtKmb9xNyu2R967\nendVPflz5vLyF9v32Vbn8bG9LPBe+NaD87nn7fVU1LXe3DX89bxn3joAlWqKiHRAJNWjLwLlwCKg\nPr7h9GwdSdgWbS5lzc5Kbnh+GU9dcSRH7Deo3WPfOncVXztsNAeMyGmy7eMNJQDUefw0RhCLr8Lq\nei9H/nZei1XB63ZVMmpAJllpkbyFWrZuVyWn3PU+GalJ1HkCX+j+Vl7CXRV1APz46aV8smE3d3x9\nRqfP25KkYPLgd/DF1jIu+NsnXHn8BOacMSWm52nPC0sCydGf563jO0fl43eOQ37zJn+68GDOO2RU\nk30XbwmUSnr9jl0VgT/Nep8PSN3nuPfMW8fdwUQJ4OqTJnLYuAFkpCZzZBvvvXVFVQA8+nEBZx80\nssm2/3l8Ee+tLabg9rNCJWctJbqPz9/M0sIynl0Un4RbRKSviOQbd7Rz7vS4R9ILPD5/c+iLtD3n\n/+Xj0ONXlu1oMWlzznHPvPUcP3kIe6rr+ceHm3j+820svumUJvs1fk+GF1oU7K5m1IBMBvdLb/H8\nu6sCX/KDWtkOsL6oqsWEDeCUu95nzMBMLj16PN+bPb7VY7Rl3uoigFDCBnt7FTbnC8vmnlm0NSZJ\nW3FlPT94YjF/vuiQUJbrXKBtIsCqGFZvbyiu4idPL+Wx780iJ6NpUvX+2mLqPD5W79xbOvrHN9fy\nxspdnDp1GAAPf7Rpn6StUYPXT0qwmrGuoeXSwcYSxEZ/fnt96HHB7Wc12Vbn8XHTC8u57rTJoeS/\npdvyXnD6Nq/PH+oFWufZt1r+pheWtxiTiIh0TCRJ28dmNt05tyzu0UQp0RPGR9IeqrbBxz1vr2t3\nPwiUctz11lruemttaJ2nhVkYfP7Auq//9RO27AlU813z7yVkpyWz4tenU1XvJSXJyEhNDj3nsFvf\nAuCd605g/OBsAJYWljF9VC5JSY2lTm2XHBbuqeU3L6/km7PGdKrEraXSwJZKK1ftqODVZU07LOfP\nmcucM6Zw2oHDQ/FDILnbVVHHz579gumjcxnaP51Lj85vsRrun/M38+mmPXz1/o9Dya3DhdpZ+Voo\n9vP7HX7nmL9xD1X1Xt5evYvff63tBPKo2+axozxQUvjmil2ce/BIUsLG9mttGrRl28pZti1Qldla\nCSQEei83xtw4diDAlt01/OW9DZwweUiH2le+vmInzyzayqqdFZw2dXjw/K0HcNurqykJ/hNQ5/FR\n5/Gxsbia/MFZfLppT8TnFRGRtkXyTTsbuNTMNhGoHjXAOecOimtkneCc+y/w35kzZ16e6Fiaq673\n8tLS7fz9/Y1sLKlusq21dj1e375flPU+P3UeX5MErDG5aEzYQucMdkaYdsvrjMzN4OPrT97neCf+\n4V3uunAGo/KyuOBvn/Dz06fwvydMANpOFMLVefxkpe273ud3rC+qYvLw/i0+L6mF624pUTrj7g9a\nfP7tr67m9ldXs+G3Z4aSlt+/tpq/vb8RCLSTg8BQJWdNH8ER+w1iV0UdKUnGoH7poXNtLa0NtdN7\nbvE21u2qajWW655dynOLmzb4v/W86aSl7Ns8tKbBy3trikMJG8BPnlnKws2l3PbV6S1eU2vCkyaf\n3/HV+z8KLYe3N6wJ64Dyo6eXsGhzKU8u2BLxeT7duJtXggny8m0VLN8WSPbaeiv8I6wd4sMfFbB6\nZyWLNpdy7Zcm8ae3IvsHRURE2hdJ0nZG3KPoxWbe+hZXHr8ft85d1eo+Ly3dzpiBWftUM7ZUutHg\n9TPlpteaVGm1lFw0t728jpoGL0nWtMQN4EdPLQ09Dm+k3l5JW6N5q3bx6vKd/P07M1m7q5J731nP\ngk17KK4MlL68fPVsDhiR06Sn4A3PL+OJT/dNJsIvxevzc+Pz7Vet1Xt9oZK+xirXcI99spnHPtlM\nwe1nccRvAwnOyVOGtrgvECrd2lhSxWcFezg8fyAQ6KzQPGED2F1dz9D+Gfv0hLzlxRU800I7ro3F\nVe1eU3ON9/jdNUW8tGQ7S7e23JmgLpi0bS2tobS6ocPn+ct7G3h3TfE+6yNtrvmvsHvamPCJiEhs\nRDLkx2ZgDHBS8HFNJM+TgJKq+jYTNoA91YFen821lTSFN/hurR1Yc1Nvfp3jfv8Ory1vfWzk5iU6\nkfjps1/w9uoiKus8XP3k58z9YkcoYQM4+88f8tWwNnxAiwkbBBKapYVl3DNvHW+vLuKphYXtnv/O\nN9aGhllp6zV7f+3eZKS1hC3crop6rv33EgCKKuuY/bt3WtzvqNve5qfPLN1n/eZmJZ+NKoM9LBu8\nfh7+qOXess01XtelD3/Gc20M7VHT4MM5x+zfvbNPiW5r7n93Pflz5nL5YwupaqX3Z/ir6vH5Ofb3\n7Q/7UtjK9YuISOe0W9JmZrcAM4HJwMMEuqb9EzgmvqH1PY/P38z5h44iKy2Fpz8rbFLt1Fytx0dq\nsF1UW8nVD5/8vMlyUWU9V/5zcav7+8OazLVUPdsWr9+1Oj5cY1LVng3F1Zx7X6Dqr6Uqx5Y8+OEm\nHvxwU6D0sY2QW2s71pbSmga2ltbw8frdbe733OfbOGRsHhcfOQ4zI3/O3Fb3rawPDMvx9w82csfr\nayKKw+8CyX17dlbURZysNfr9a4EY3ly5iymtVGXjHH6/42f/+QIItGdsT2GpkjYRkViKpHr0K8Ah\nwGIA59x2M2vlkz2xEt0RIVo3vbCcldvLue2rB4W+HFtT1+AjJyOVwj01eNpIrl5auu/4Wm1pLNHx\n+x01DR2bYva5xVupbaH3YGc1tNDpoi2//u9KPP7Yjq1W0+BrtYStuZteXMEJk4cyZmBWm/sV7qkl\nf85cDhyZ0+Z+4dYXVXHob95sd79fRNlTs7XE0O9gV2Vdh4btqInxAM8iIn1dJElbg3POmVnj3KPZ\n7T0hUbpzR4RIbd4dWenEW6uK+HTTbl5c0rGkrD0VdR4KSqq5LthgviN++8rqmMbSUQ9FWNUYTxuK\nqxiWkxHRviu2d782X0WVLQ/FuGxbOUfdFvuZMEREJHKRJG1Pm9nfgDwzuxy4DPh7fMPqu6rrvfwn\ngtKMG56Pzwgs8zfu4YQ/vBuXY0Ng4Nov3/tR+zv2UE+1U60tIiLSWe0mbc65P5jZKUAFgXZtNzvn\n2q+nkU5ZurWcn7TQqL032FRS3aHhJ3qiV5fvTHQIIiLSS0U0ImowSVOiJlE58Q/vMjzCqkMRERFp\nSkN3SJfaWVHX/k4iIiKyDyVtIiIiIj1AREmbmWWa2eR4ByPSU80aP5Ch/dMTHYaIiPRi7SZtwbHP\nlgCvBZcPNrOX4h1YZ5jZOWb2QHl5y1P8SOINyEpNdAhx8e/LjyQ9tfcWXL9w1TF8+PMT+du3D0t0\nKCIifVYk3zK/BGYBZQDOuSXA+LaekCjOuf86567Izc1NdCjSis9vPpXxg7vtUH8d9qUDhvLRnJNI\nSrLQ/Jx/+dahnHfwyMQGFkODstM4eEweowdkMXvi4ESHIyLSZ0WStHmcc82Lrjo2v5FImOz05PZ3\n6iEumz2eUXmZwN5J1ccMzKKkKjCzwNUnTWx9aqgeYkhYtW92ekQdzkVEJA4iSdpWmNlFQLKZTTKz\nPwMft/ckkdZkp/WeL/78QfuWGvZLT+Hnp0/hkLF5XHn8BF679jguOWpcXM5/1vQR+6z75PqTWlzf\nWd+b3S0L1kVE+pxIkrargQOBeuAJoBy4Np5BSXylJSe27VU8SmuevfIofnZ6+31lLjsmtglI+Lhz\njSVSaSlJTB+dy/P/d0zoWn917jSO2m9Qu8frF9x/YHYa/7hkJnN/OJv7LjqUrx02ep99/3rxYdx7\n0SEsuPHkJutTk5OorN87b2w0JX0bf3smX585psm6hy6dyZnTh3f6mCIi0jmRfHtOcc7dCNwY72Ck\na1x5wgRm5Q/k7nlrqfP4WbatnGu/NIl1RVXM/WJHXM550pShobZssUzapgzvz+qdlRw8Jo+Z+QO5\n8rgJvLeumHEDszjpj+8B8IMTJ7J6ZyUz8wdw5fETYjJHaXKSserXp5OUZKF1f/v2Yby1ahcjg9Wl\nzeVktn0eSz/7AAAgAElEQVTd5x86mh+ePJHj73iXX597ICcfMAyAA0fmctZBI7jw8DF8vH43d721\nlu8cNY7TpwUSp6H9M3jssllkp6ewtLCMwf3SqQ4mbQ9dOpOjJwxm8ZZSPt24h7vnrYvo+g4ancuf\nLjy4yfU1OmnKME7YfyjXpy/jqYWFER1PRESiF8m35x/NbDjwLPCUc255nGPq0846aMQ+idMPT5rI\nPW+vj9k5UpKM2ZMGM3tSoFG5x+cnNTmJldsrmPvFDm4+eypTRvTnor9/GpPznTB5CA9denhoOTst\ndm3a/nrxYeRlpZISLD1MSjJOnDwUgBvPPICJw/qFljsqMzWZaaNy+KygdJ9tYwZkkpbStMRyWE4G\n3zqi9WrQS47K5/UVu7jrwhk8/dlWPtm4u8n2lCRj3KBs1v+/M0LXE+7w/IEcnj+Qi44Yu08v3OP2\nHwLAYeMGAIFE9cp/LuLw/IFkpCZz9ITBHD1hMNecPIn9bnil3WsfmZvJfkP6tbo9Kcn43dcO4v11\nxewo79yAyfd/69BOPU9EpK9qt57MOXcicCJQDPzNzJaZ2S/iHlkfMWNMHu/99ITQ8n0XHcrNZ09l\n8rBAldb0Ubn8+NTOD5F3+1en8+yVRzVZl9ys9CQ1mCBMHZnDghtO5rvH5DMgK63T5ww3//qTeeS7\ns5qsy82M3bAfA/ulkddKrJcft1+nEzaAVb85nWeuPJoXrzqmyfpbz5vGvy4/ssPHO3riYApuP4uv\nHDKag0YHejj//PQp/ObcAwFISbbg77b/LIf0T293nxOnDGXNrWfQP6Ppa52UZLx27bFcc/KkNp/f\n/D3Smuf/7xj++PUZEe0b7v5vHcqZMWx3JyLSF0TUuMk5t9M5dw9wJYEx226Oa1R9yHWn7s+4Qdk8\n8t3DQ19+l80ez6vXHMv3Z48PlUY88f0jQqUpkbr82PF89dDRzMwf2GR9anLrX8hDczIwM/JiMJ5a\nwe1nMTx337lGrzppIt86Yiw/OWV/MqIc26x/J6paxw/OZsaYPNbeegYHjc4lf1BWaNuVx0/YZ//G\npDY12Xjo0plcfOS4UI/RzkpP2XtMrz/Q7TQlwkQpWlOG5zB1ZE6b+0SatA3PzeDLHRze5InLj1DC\nJiLSCe1+45nZAcCFwPnAbuAp4CdxjqtP2PjbM0Nthk5oViKUlGT84uypoeWjJw5maE4GX7rzvYiP\nf+NZU1tcn5zUfqI0IjeT+defzJG3zYv4fJHKyUjl/31lemj5j2+u7fAxvjlrLLecMxWzjic671x3\nQujxSz+YDUD+nLkAzDljCl85ZBR7qhtC+zSWgI0dmMVJU4Z1+Hwt+Z/jJ1Dn9XPxkeNYsT0wos4x\nXTgG2pTh/Rk/OJsj9xvE859vpc7jZ/Kw/hw4MofnPt/GqAGRJ6WpEXZs+dpho7nw8DEc3uyfCBER\niUwkxRQPEUjUTnPObY9zPFEJzt5wzsSJExMdSkRaauTdlvBSqV+cdQC3zl21zz7XnDypxcbmv//a\nQZTVNPDQhwURDwcxPDeD/hkpVNZ529+5k/bvRM/G9356AqMHZEVcGhSJ8w8dTVlNIFGb3CymxvHX\nOpMgtiY7PYUbzjwAgMPGDWTpzaeS24WzRYwblB1KXr9z1DjOuPsDbjt/OoeMyeP4yUM4Y1rHSsIK\nbj+LW15czqOfbG6y/tCxeUwdmUNptYdbz5tGRmrvGaNPRKSrmXO9b5zcmTNnuoULF8b1HI0lM9Eo\nuP2sDu1fUlXPzFvfIjczlfnXn8wBN7/W4jEbY+vo8VuyvayW7WW1XPTgpzR4/UBgFgCf33H2QSP5\n1X9XUNEsqfvWEWM5Z8ZIjoxgiAuAhQV7mDEmD4A3VuziqicWt7n/6t+c3qVf/qt3VnD6nz5g0tB+\nvPnj47vsvD2Nz++orPNw8K/fBGi1Q4WIiDRlZoucczPb26/VkjYze9o5d4GZLaPpDAgGOOfcQTGI\ns886+6COt+lJC7WDSmrSc/HW86bxixfi06l3ZF4mI/Myeee6E9hUXM3REwY1KSE8fvIQFhaUcuU/\nF4XWXX3SpBbbsrUmvM3dWQeN4Kon2t4/PaVrE4Hxg7OZNiqHm1qpbpaA5CQjLyuNn5yyP4ePH6iE\nTUQkxtqqHr0m+PvsrgikN2ocQwzg2i9NYndVA4/P38yVx09gzhlTOny8/ukp/OhL+3Pm9OGhqsFR\neZlcfOQ43ly5KzRu15OXH8mIDiRNkRiVl9li4/vB/dI5fdpwHrp0JmMHZlHn8XcoYWvP6AGZ3HDm\nASzfVs76oio+XF8S02rKSKSnJPPy1cd26Tl7sqvb6ZkqIiKd02rS5pxrHCzs/5xzPw/fZma/A36+\n77Mk3C/Omsoziwp5ccl2xg3KCvUO7GwBhJlxzZf2fiE+dtmsUPurRy/bO6zGURMiq5aMpVg10A93\n/RlT+PrMMQzMTuPM6SPYVFLNqh0VMT+PiIhITxBJ+nBKC+vOiHUgvVHzTpr+xgbtxKak6Lj9hzAs\nJ7Ylat3J/xw/gYHZe8dgGz84W0NFiIhIn9VWm7b/Bf4P2M/Mvgjb1B/4KN6B9QbJzarx/MFOH100\nHJeIiIj0Im21aXsCeBW4DZgTtr7SObcnrlH1Es2HpPDHYeiI3uj1a4+juLI+0WGIiIh0K221aSsH\nyoFvApjZUCAD6Gdm/ZxzW7omxJ4rKalZRWiopE1JW1smD++/z1hpIiIifV27bdrM7BwzWwdsAt4D\nCgiUwEkrGifz3rd6NPBbOZuIiIh0VCQdEW4FjgTWOufGAycD8+MaVQ9WcPtZDOmfDgSqR8MHuFOb\nNhEREemsSJI2j3NuN5BkZknOuXeAdkft7ctam9tTbdpERESksyKZe7TMzPoB7wP/MrMioDq+YfVM\n4wZlAZAanGDc599bzmYYTm3aREREpJMiKWk7F6gFfgS8BmwAzolnUD3RlOH9eekHs4G9vUa9/qbz\nuqp6VERERDqr3aTNOVftnPM557zOuUedc/cEq0u7hJlNNbOnzewvZva1rjpvR00Y0o/czEAHhMaZ\nD3x+xxnTAoPBThuVywmThwJEPJG6iIiISKO2BtetpIWJ4tk7YXxOewc3s4cIzF1a5JybFrb+dOBu\nIBl40Dl3exuHOQP4s3PuAzN7CXi2vfMmQniNZ3Z64GV1znH6tOFsuu1MzIyJQ/ux8bdnNplwXURE\nRCQSbY3TFouBsh4B7gUea1xhZsnAfQSmx9oKfBZMxpIJDOQb7jLgceAWM/sy0G2LqMLbqd3xtRk8\n+nEBh+cPBJp2PFDCJiIiIp0RSUcEzGw2MMk597CZDQb6O+c2tfc859z7ZpbfbPUsYL1zbmPw2P8G\nznXO3UagVK4lVwWTveciiTcRwkvahvRP57rTJicuGBEREel12k3azOwWAkN8TAYeBtKAfwLHdPKc\no4DCsOWtwBFtnD8fuAHIBu5oY78rgCsAxo4d28nQRERERLqnSEravgIcAiwGcM5tN7Mum2PIOVdA\nMBlrZ78HgAcAZs6c6drZXURERKRHiWTIjwYXGGDMAZhZdpTn3AaMCVseHVwnIiIiIq2IJGl72sz+\nBuSZ2eXAW8CDUZzzM2CSmY03szTgG8BLURwvJDhP6gPl5eWxOJyIiIhItxHJOG1/IDDMxn8ItGu7\n2Tl3TyQHN7MngU+AyWa21cy+55zzAj8AXgdWAU8751Z09gKaxfpf59wVubm5sThcm5p3Aq1t8MX9\nnCIiItJ3RdR71Dn3JvAmgJklmdm3nHP/iuB532xl/SvAKx0JtLsxM3B7m85tKtHMXiIiIhI/rZa0\nmVmOmV1vZvea2akW8ANgI3BB14UYuURWj64vruryc4qIiEjf0Vb16OMEqkOXAd8H3gG+DpznnDu3\nC2LrsK6sHm3uyzNGdvk5RUREpO9oq3p0P+fcdAAzexDYAYx1ztV1SWQ9xJ+/eQgnThlKRkokfTpE\nREREOqetTMPT+MA55wO2dveEraurR39w4kTOmTGSfukppCQraRMREZH4aSvTmGFmFcGfSuCgxsdm\nVtFVAXZEIqtHRUREROKprQnjk7syEBERERFpner0RERERHoAc673TNNpZucA5wAXAuvifLrBQEmc\nz9Gd6fr77vX35WsHXb+uv+9ef1++dojv9Y9zzg1pb6delbR1JTNb6Jybmeg4EkXX33evvy9fO+j6\ndf199/r78rVD97h+VY+KiIiI9ABK2kRERER6ACVtnfdAogNIMF1/39WXrx10/br+vqsvXzt0g+tX\nmzYRERGRHkAlbSIiIiI9gJI2ERERkR5ASVsnmNnpZrbGzNab2ZxExxNPZjbGzN4xs5VmtsLMrgmu\nH2hmb5rZuuDvAYmONZ7MLNnMPjezl4PL483s0+B74CkzS0t0jPFiZnlm9qyZrTazVWZ2VF+6/2b2\no+B7f7mZPWlmGb35/pvZQ2ZWZGbLw9a1eL8t4J7g6/CFmR2auMij18q13xF8739hZs+bWV7YtuuD\n177GzE5LTNSx09L1h237iZk5MxscXO5V9x5av34zuzr4HlhhZr8PW9/l919JWweZWTJwH3AGMBX4\npplNTWxUceUFfuKcmwocCVwVvN45wDzn3CRgXnC5N7sGWBW2/DvgLufcRKAU+F5CouoadwOvOeem\nADMIvA594v6b2Sjgh8BM59w0IBn4Br37/j8CnN5sXWv3+wxgUvDnCuAvXRRjvDzCvtf+JjDNOXcQ\nsBa4HiD4OfgN4MDgc+4Pfj/0ZI+w7/VjZmOAU4EtYat7272HFq7fzE4EzgVmOOcOBP4QXJ+Q+6+k\nreNmAeudcxudcw3Avwnc0F7JObfDObc4+LiSwBf2KALX/Ghwt0eB8xITYfyZ2WjgLODB4LIBJwHP\nBnfptddvZrnAccA/AJxzDc65MvrQ/ScwR3OmmaUAWcAOevH9d869D+xptrq1+30u8JgLmA/kmdmI\nrok09lq6dufcG845b3BxPjA6+Phc4N/OuXrn3CZgPYHvhx6rlXsPcBfwMyC852KvuvfQ6vX/L3C7\nc64+uE9RcH1C7r+Sto4bBRSGLW8Nruv1zCwfOAT4FBjmnNsR3LQTGJagsLrCnwh8YPmDy4OAsrAP\n8t78HhgPFAMPB6uHHzSzbPrI/XfObSPwn/UWAslaObCIvnP/G7V2v/va5+FlwKvBx33i2s3sXGCb\nc25ps0194vqB/YFjg80h3jOzw4PrE3L9StokImbWD/gPcK1zriJ8mwuMG9Mrx44xs7OBIufcokTH\nkiApwKHAX5xzhwDVNKsK7eX3fwCB/6jHAyOBbFqoPupLevP9bouZ3Uiguci/Eh1LVzGzLOAG4OZE\nx5JAKcBAAs2Dfgo8HaxtSQglbR23DRgTtjw6uK7XMrNUAgnbv5xzzwVX72osCg/+Lmrt+T3cMcCX\nzayAQFX4SQTaeOUFq8ugd78HtgJbnXOfBpefJZDE9ZX7/yVgk3Ou2DnnAZ4j8J7oK/e/UWv3u098\nHprZpcDZwLfc3sFN+8K1TyDwD8vS4GfgaGCxmQ2nb1w/BD4DnwtWAy8gUOMymARdv5K2jvsMmBTs\nPZZGoCHiSwmOKW6C/1H8A1jlnLszbNNLwCXBx5cAL3Z1bF3BOXe9c260cy6fwL1+2zn3LeAd4GvB\n3Xrz9e8ECs1scnDVycBK+sj9J1AteqSZZQX/Fhqvv0/c/zCt3e+XgO8EexIeCZSHVaP2CmZ2OoHm\nEV92ztWEbXoJ+IaZpZvZeAIN8hckIsZ4cc4tc84Ndc7lBz8DtwKHBj8Xev29D3oBOBHAzPYH0oAS\nEnX/nXP66eAPcCaBXkQbgBsTHU+cr3U2gaqQL4AlwZ8zCbTrmgesA94CBiY61i54LU4AXg4+3i/4\nB7oeeAZIT3R8cbzug4GFwffAC8CAvnT/gV8Bq4HlwONAem++/8CTBNrveQh8SX+vtfsNGIHe9BuA\nZQR62Sb8GmJ87esJtF1q/Pz7a9j+NwavfQ1wRqLjj8f1N9teAAzujfe+jfufBvwz+Pe/GDgpkfdf\n01iJiIiI9ACqHhURERHpAZS0iYiIiPQAStpEREREegAlbSIiIiI9gJI2ERERkR5ASZuIiIhID5DS\n/i49z+DBg11+fn6iwxARERFp16JFi0qcc0Pa269XJm35+fksXLgw0WGIiIiItMvMNkeyn6pHRURE\nRHoAJW0iIiIiPYCStijlz5nLLS8uT3QYIiIi0sspaYtCg9cPwKOfRFQVLSIiItJpStqiUFbbAECS\nJTgQERER6fWUtEWhrMYDQP+M1ARHIiIiIr2dkrYolFYHStr6pffKkVNERESkG+l00mZmx5hZdvDx\nxWZ2p5mNi11o3V9pqKRNSZuIiIjEVzQlbX8BasxsBvATYAPwWEyi6iHKagIlbUraREREJN6iSdq8\nzjkHnAvc65y7D+gfm7B6hrLaQElbtqpHRUREJM6iyTYqzex64GLgODNLAvpUi/zKukDSlpaspoEi\nIiISX9FkGxcC9cD3nHM7gdHAHTGJqgf4338u4r53NgDgdwkORkRERHq9TpW0mVky8KRz7sTGdc65\nLfShNm2vLt8Zeuz1+xMYiYiIiPQFnSppc875AL+Z5cY4nh7Jp6I2ERERibNo2rRVAcvM7E2gunGl\nc+6HUUfVzQX6X+zl9SlpExERkfiKJml7LvjTIWY2hkA16jDAAQ845+42s4HAU0A+UABc4JwrNTMD\n7gbOBGqAS51zi6OIO2rNcjaVtImIiEjcdTppc849amaZwFjn3JoOPNUL/MQ5t9jM+gOLgqV1lwLz\nnHO3m9kcYA7wc+AMYFLw5wgC48Md0dm4Y6F5iuZRmzYRERGJs2hmRDgHWAK8Flw+2Mxeau95zrkd\njSVlzrlKYBUwisB4b48Gd3sUOC/4+FzgMRcwH8gzsxGdjTsWmlePqqRNRERE4i2aIT9+CcwCygCc\nc0uA/TpyADPLBw4BPgWGOed2BDftJFB9CoGErjDsaVuD6xKmeY6mNm0iIiISb9EkbR7nXHmzdRHX\nE5pZP+A/wLXOuYrwbcGZFjqUCZnZFWa20MwWFhcXd+SpHeZQSZuIiIh0rWiSthVmdhGQbGaTzOzP\nwMeRPNHMUgkkbP9yzjV2ZtjVWO0Z/F0UXL8NGBP29NHBdU045x5wzs10zs0cMmRI564oQs07IqhN\nm4iIiMRbNEnb1cCBBGZFeAIoB65p70nB3qD/AFY55+4M2/QScEnw8SXAi2Hrv2MBRwLlYdWoCaHe\noyIiItLVohny4yzn3I3AjY0rzOzrwDPtPO8Y4NsExnhbElx3A3A78LSZfQ/YDFwQ3PYKgeE+1hMY\n8uO7UcQcE43Vo9+cNZaymga+2Nq8llhEREQktqJJ2q5n3wStpXVNOOc+BKyVzSe3sL8DrupMgPHS\nWLA2fnAWG4qcprESERGRuOtw0mZmZxAo+RplZveEbcohMAZbr9c45EeSGcnJpupRERERibvOlLRt\nBxYCXwYWha2vBH4Ui6C6u/AcLTXJ8CppExERkTjrcNLmnFsKLDWzJ5xznjjE1P0FczQzIzkpCZ/G\naRMREZE4i6ZN2ywz+yUwLngcI9AErUMD7PZEjR0RkgxSkk1DfoiIiEjcRZO0/YNAdegiwBebcHqG\nxtpQA5KT1KZNRERE4i+apK3cOfdqzCLpQUIdEZJMbdpERESkS0STtL1jZncAzxEYYBeAxsnge7Om\nJW1JOAd+vyMpqbWRTERERESiE03SdkTw98ywdQ44KYpj9ghub08EUpIDiZrH7yc9KTmBUYmIiEhv\n1umkzTl3YiwD6Ukap7FKMkgJlq6pXZuIiIjEU6eTNjO7uaX1zrlfdz6cnsGFqkeN5GDS5vE6SEtg\nUCIiItKrRVM9Wh32OAM4G1gVXTg9Q/iQH2kpSQAa9kNERETiKprq0T+GL5vZH4DXo46oB/DvbdJG\nanIgaXvqs0JOnDyUqSNzEhiZiIiI9FZJMTxWFjA6hsfrthqH/DCzUNJ2x+truOyRzxIZloiIiPRi\n0bRpW0ZoQieSgSFAr2/PBuFt2iA1ee8wH7mZqYkJSERERHq9aNq0nR322Avscs55o4ynR3Bhc4+m\nJe8trBw7KCtBEYmIiEhv1+nqUefcZiAPOAf4CjA1VkF1d+EdEVKTY1nDLCIiItKyTmccZnYN8C9g\naPDnX2Z2dawC686adERI2fsSNrZ1ExEREYm1aKpHvwcc4ZyrBjCz3wGfAH+ORWDdWWjuUbMmbdrq\nvRr2Q0REROIjmro9A3xhy77gul4vfPKD8DZtH6wr4eP1JQmISERERHq7aJK2h4FPzeyXZvZLYD7w\nj5hE1e3tO+RHo4se/DQRAYmIiEgvF83gunea2bvA7OCq7zrnPo9JVN1c+Nyj6oggIiIiXSGacdqO\nBFY45xYHl3PM7AjnXK8vavKHzT2altInaoRFREQkwaIpJvoLUBW2XBVc1+tpyA8RERHpalF1RHBh\nY1w45/xE1xu1x2icG95aSdry58xlT3VDF0clIiIivVk0SdtGM/uhmaUGf64BNsYqsO7MhWbv2rcj\nQqOlW8u6LiARERHp9aJJ2q4Ejga2AVuBI4ArYhFUd9e0I0LLbdoqaj1dGJGIiIj0dtH0Hi0CvhHD\nWHqM8LlHWytpU9ImIiIisRRN79EhwOVAfvhxnHOXRR9W9xZJR4Q91UraREREJHai6TjwIvAB8BZN\nZ0bo9ZrMPdpK9WhRZV0XRiQiIiK9XTRJW5Zz7ucxi6QHaew0axhmLSdtNQ19Ko8VERGROIumI8LL\nZnZmzCLpQUJ9R9sYV9fj0+TxIiIiEjvRJG3XEEjc6syswswqzawiVoF1Z6GStmDW9v3Z47n+jClN\n9lHSJiIiIrEUTe/R/rEMpCcJH/ID4BdnTwVgwpB+fP+xhQB4fa6lp4qIiIh0SqdL2izgYjO7Kbg8\nxsxmxS607it87tFweVmpoccev5I2ERERiZ1oqkfvB44CLgouVwH3tfckM3vIzIrMbHnYuoFm9qaZ\nrQv+HhBcb2Z2j5mtN7MvzOzQKOKNmcbq0aRmbdrCkzavqkdFREQkhqJJ2o5wzl0F1AE450qBtAie\n9whwerN1c4B5zrlJwLzgMsAZwKTgzxV0kwnp/XtnsWoiN3Pv5Te2aXt/bTHLtpZ3UWQiIiLSW0WT\ntHnMLJlgZ8rgYLvtFi85594H9jRbfS7waPDxo8B5YesfcwHzgTwzGxFFzDHROLhu8+rR3My9JW0N\nPsemkmq+89ACzrn3wy6NT0RERHqfaJK2e4DngaFm9v+AD4HfdvJYw5xzO4KPdwLDgo9HAYVh+20N\nrtuHmV1hZgvNbGFxcXEnw4hQs44IjdJS9r6cSwvLOPEP78Y3DhEREekzOp20Oef+BfwMuA3YAZzn\nnHsm2oBcoMFYh1vxO+cecM7NdM7NHDJkSLRhtMkfNvdopPLnzGVbWW2cIhIREZHeLpqSNpxzq51z\n9znn7nXOrYriULsaqz2Dv4uC67cBY8L2Gx1cl1Dhc482t+m2Mzlress1uAs27Y5nWCIiItKLRZW0\nxdBLwCXBx5cQmNe0cf13gr1IjwTKw6pREyZ87tHmzKzV+Ui3l9VRUFLNu2uKQj1QRURERCIRzdyj\nnWJmTwInAIPNbCtwC3A78LSZfQ/YDFwQ3P0V4ExgPVADfLer423J3oSr5eQsJbnlXPiO19dwx+tr\nALjrwhl85ZDR8QhPREREeqFOJ21m9rvmE8a3tK4559w3W9l0cgv7OuCqzsYYL40pW0vVo0CrJW3h\ndpTXxS4gERER6fWiqR49pYV1Z0RxvB6j+dyjzaUktf+yprSW8YmIiIi0oMMlbWb2v8D/AfuZ2Rdh\nm/oDH8UqsO6s+dyjzaVEUNKWHEFiJyIiItKoM9WjTwCvEhjqY07Y+krnXPNBc3ul1uYebZQcwVAg\nKmkTERGRjuhwcY9zrtw5VxBsm7YV8BBo5tXPzMbGOsDuaG/1aMvbvRFMFh+es20trWF9UWUsQhMR\nEZFeKpqOCD8AfgnsYu/0VQ44KPqwurdQ39FWkraGCCaLv+nFFXywroQHvjOT2b97B4CC28/qXDzO\n4XeQrNK7mPH5nV7PFviD/5Ak6bUREely0TSsuhaY7Jw70Dk3PfjT6xM2CCtpa6V61ONtP2kDeGPl\nrla31TR4Kamqb7JuW1kt3hYSwvvf3cCEG16hpsEb0XmlbVtLa5hwwys8//nWRIfS7Vz71BL2u+GV\nRIchItInRZO0FQLlsQqkJwl1RGjl1fNEUNLW6JKHFrS4/sv3fsTMW9/iiscWsnxbOXUeH8fc/jY/\ne/aLffZ98IONAFTVKWmLhc27awB46rPCdvbse15auj3RIYiI9FnRJG0bgXfN7Hoz+3HjT6wC687a\n64jg8UU+28F7a/dObl/n8XHPvHVc9Pf5rC+qAgKlcdc9s5SKWg8Az32+7yxetR5f4Lx+xzuri5h4\nwytU1nkijqG7qPf6Enbu7WW15M+Zy5LCMjJSkwGoaUhcPJ3hi6AtZaw45wLV8l14ToGnFxZSXtvz\n/rZFJDaiSdq2AG8CaQSG+2j86fUa5x6Npk1bSxZs2sOdb67l4w1N5yhNMqOqvuVStAc/2EidJ3C+\neo+PP81bh9fvQklfvOyqqKMhwmrgSGwvq2XyL17j1pdXxuyYHfHR+hIAHv9kcyj5qW7lNe+O/v7+\nRibc8EpcYw5/D3r9jl/9d2WoqrTB62enBoyOq4KSan727Bf84InFiQ5FRBKk00mbc+5XzrlfAXc0\nPg4u93rtjdPWkerRcN9ppap05Y4KrntmaWh5aWEZEChZuXXuqtD6eq8/VPYXz/IPv99xxG/nce1T\nn8fsmEWVgfZ7D364aZ9tJVX1rNxeEbNzAWzeXc3SwjIueWgBJVX1pAanHvP4/KESvw3F1dz84vJQ\nG8ZFm/fEpd1g4Z4aCkqqgUDyuqG44wn388ES2NU7Y/s6hQsv4an3+nnk4wIgUEJ8y0vLOfK2ee2+\nPkn0vEgAACAASURBVHe+uZZ/L9gStxjjbUd5LUuCf3/hPtmwm2v//Xlc5xT2B4+9bFvHW6U0eP18\nunF3+zuG6ex7UUTip9NJm5kdZWYrgdXB5Rlmdn/MIuvG/O3MPTogKy3m51y8Ze8Xxbn3fcTG4ioe\nDX5pNqrz+EKJ5Luri1i3K7bDiLy4ZBslVfXUB0vYXlm2M6rjbd5dTf6cuTyzsJAP1xXvs72izsP+\nv3iVC/76CRf/41Om3/I6D7y/gTdWRHdev99x/B3vcu59H/He2mIe/bigadLm2Zt0P/bJZj5YV8Ku\nijrO/8snzPnPshaP15kv6/VFVUy44RWO/f07nPCHdwE4+va3OfmP70X0/O1ltby6bAcAE4b2A2BF\nDJPbNTsr+Wh9SajzS21YQlbv2Vt1XFHn4f21gZLKooq9nWfKaz08u6hpZ4575q1jznP7voadsbW0\nhtfD3gvxTJgg8Pd11G1vc959+44h/u1/fMoLS7aHmiq0pbbBx5MLtrQYb2O1c0saS/DLajzt7tvc\nb15eyYUPzOfRjwtYs7P9zwXnXOi9+O8FW9TJqQ9YvbOC/Dlz2by7OtGhSBuimTD+T8BpwEsAzrml\nZnZcTKLqIVorafvVuQdyeP5Abng+Nl9OLTmphS/2eq8/NLXWPW+v556313d6GJFwSwrLeG9NMXe9\ntbbF7RuLq3hxyXau/dKkVqf2asmbwd6zP22hc8VH60t4emEhDV4/G0v2foj89pXVAPz89CmcMnXo\n/2/vvsPcqK6Hj3/PdnevK/a67NrGHRcwrhgXINhgaiDUEKpDQocU2o8QQhohAUJ9Te+9hg7GxhiM\nccG44G7ce2/bJJ33jxnJklbaqh1tOZ/n2WdXo9HqXo00OnPLuXRr04QnvlrF0C4t6ZvTrMznnLJ0\nK5c8PStiW4pIaBWL3QeL2bw3sptvzY4DHNYsC4CFYa0cSzbvZfLirfzrk6VcOiKPO07pHbpvx/5C\nGmak0SAjNW5Z3pizPmIc2t8+XBx331jOfmwGG3bns/JvJ9G8Qbpb1oMV+h/Rpi7dyo79RSzbuo//\n9+Wq0PY7JvRmUG526HZhWNf4vgIfLRplsGF3Pgs37uGwZllkpady46vzmLxkKwM6NqNbm6qNnJi1\neidLN+/jwqGdQ9tOe+hrdhwo4qe/n8R9ny3jv1+sYMVfx5OWWj2rjYQHpNHSUgVfQNlX4CMjNYXt\n+4tC75loD0xezmNfriS7YTrj+raLuG/c/V+x82ARs247nq+Wb2PL3kKWb93H5cd0KTEc4Za3FvDK\nrHXl+owHWwf/9N4iwEkv9O68DTRrkM7oHm0i9l20cQ8n/3d66PbNby1gwYY9/PWMI8p8ntrqg/mb\nyEpP4bhebcvcd9bqnSzbso8LhnQuc9/a5K25Tmv9hws285vRXZNcGhNPlc5uqho9va52jdyupEAZ\na482zUrn/CGVyzPctXWjSper0BeIGUiu2LqfjbvzeXjKCl6aWbGuqdmrd3L6w1/HDdhueu0HTv7v\ndB6YvLzCA6Qz0+MHNBc8MZN358WfqfjPj5dw9mMzUHW6iCc8OJ13vt/ANrebdenmfaG/w70+u+SM\n0NQUwedOHpmxage3v7Mw4v6DRX72uTNzw/OTXfn8HP71yVIAnvo6slv3qLs/5/wnvgXg/s+XhdKH\nrNt5kMufncXeguISX8KTpq2iNP6Act0r34e6xzfszgecWcMFbgtPVSdPXPz0LG56/YeIgA3grvd/\njPjfsYI2gKtf+j4UFP/kBtuJmKtw9mMzuP2dhaEgd+7aXew4UBQqy1NfrwZgbxkzqP87eXmJ1r/y\nOhDW2qSq7D5YFAri092p5PsKinlr7gbG3Ds1butU8HMS/v5UVaYv387SLYfet7988jt+5x6LW99e\nEPGaA7zizm5WVfbkFzN/fWS37R3vLuR/7mzfWONPr3tlHhdHXcAAfLuq5MI2tWW84vz1u9lbiUlY\nV700l8uenV2ufc9+bAa3vb2w2lt2vZYR1ttgaq6qtLStE5HhgIpIOnAdULGmglpKQ7NHE++NK4cz\n8C+fVeqx+UX+EjNa+935SYkvsvMGdyyzRezeT5byyqx1JXLFRXtz7qEvwA8WbGJ411bktSpf4Nkg\nTtC2y/0yLsveAl9Ed9T1r87j4uG53HlqH068fxpNs9KYf+eJoftnrNzB1KUlu2FTU4Qif/xg50Ch\nj90HnTKFL1EW3aLz7rwNnDYgJxRAfe92ad//+XIACooDvDV3PbNW76LfnZ/SLk5LDDhfxP/+dBnf\n/bST164cBjhj+96dt5GPF25myV/Ghb0OxaHXIb+MbqyZq3bQOCuNPu0PtUrOW7cbf0A5qnN2KY+M\nHtPmR8T5LOzNLw4FbeAEvqc8OD3UQnr6w19zxcgu3HBC94j/93t3nOY/f94PEeciaMH6PeQX+xmc\n1yJmGdbuPEiTrDTOfOSb0LaDRX4y01LYXwi7DhZFlCXoQKGPD+Zv4j+fORcfZx3VIXSfP6CkSPyL\nMHCCxPDnLPIHOOORb/hp+wFaN8lknztJ4/j/TKNlowzyi/1s21dI55YlT7ENM0rOTn54ygru/fTQ\nhVGP2z+KeMzG3fkR3fbh8ov9XPHsbL5bvZOc5g349aguXDQsl+dmrOG5GWtIT01haQWGSsTKBemv\nwQHKez9sZFiXlrRolMGpDzld160aZzDz1uNLJMj++aPfMKZHa64ee3iVn3fngSJaNs4s175OV3b1\nJKX2+QOc+tDXnNjnMK47vvL1Cg4RSeQEs0T5v3cWsnlvAY9fNCjZRUm6qrS0XQlcBeQAG4AB7u06\n79BEhNI/gH1zmpbr/y29+9AXcFO3m6syrnxhDt+tjrxKjtXykHfLh9z32TI+Xhh/bNhDU1aUGbBF\nu+3thZx4/zT+/enSUlNBrN91kD35xXHTkny+OH7S4XD+gJaYVbs3v5iNbgvU3gIf93y8hJdmrmXy\n4i2c9/i3MVuiiv2BUk9U//1iRegqPDVFKPYHyL35gxIzdK97ZR6rtu2P2ypxy1sLmLV6V+j2plJa\nL/78vx95aMqKiOMZLHuhL0DeLYcS3O4ro6Ut9+YPQvkAz5n0bUTXFzhB1c8f/abE46L9+vk5ob/v\n/WRZ6HNw0VPfhSZCBIUPlj9Y5OeByctDLYNBr89Zz+tz1nP8f77kdDcgOuWh6fzi/80A4JsV20Ot\nYk0yneBn2ZZ97IwK6g8U+shMc05l93++PNT69cRXq0Jj3u54dxF/eLNkNzzA6HuncM6kb2Pe5w8o\nf/9oMXe8G9n6WlAcCLUkRrfoBlsAt++PffHRIEZKmRe+jWwBj25V27SngE17Dr1+4YHVgUI/s9Y4\n75MNu/O5491FPPTF8tD9V75w6LiVR6yWlupOJ7MnvzjUQra3oJg97ri9Yn+Af3y0JO6F3Pb9hVz7\n8vf89sU5EeeC7fuL2JNfjM8fYIs73KHQ52fOml0RwXFVrNuVX2Lb1r0FHCzysXVvAdv2FYYmNXW9\n9UOueSVxE7fAmRQTCDitrD9u2hu3N6S80tOc77Mif4ADhT7++MZ83omRYioZnv92TWg4DTjjiP/9\n6dLQsd20J7/U9+j2/YWhc2Q8ewuK+fuHi5Oaeqo8KtXSJiKpwC9V9YIEl6dWONQ9Wvp+718zkj++\nMZ9XZ6+jQ3YDfjUslyuO7cKKrfsIKPzsvmkAZKYdanHyaumkByY7J/UTerflrtP60DQrnRtencfN\n43vG/bIpjyJfgAe/WMFpA3Lo5g6OD/fuvA1c98q8Uv9HrDFu8VzzUuSJ8K3vN0Tksntk6soy/8eB\nQl+JL8l4fty0l2tfjn/y/X7t7lJb0MrrmbBJJoU+P5lpqXHTeYS3tO06WMScNbtokJ5K6yaZtG7i\ntAR8uWxbidmDq7cfiHi/VaS7p7yBdbgR//gi9PcVzx3qigq2yAWTRAP8+vnZfLLIeY7Pf9wSasla\nvyufJlmRp63rXvmejW4A/L8fNvK/Hzay+h8nh2ZWXzCkU0SLcLg9+cWs25nPup357v8/iM+v5Lqt\nxUs37yvRVQxw1Ytlp93YEeeipzjgvNe27ivk9dnr+GTRlhLjKKPtPFAU8bk4GPYFdKDQR3pKSkSq\nobICk2/D3gsFxf7QsIlLj8mLmWfSV4Hck5XR/8+fkpoirPzbSZz0wFes35XP6n+czEcLN/PYlyvZ\nW1DMn07pzY2v/cCNJ3Sna2vn3BIM4DfuLijx+dh1sIjHvlzJpGmrmHP78WyLOh4PT1lB26ZZEa2u\nqoovoHy9YjvHdGsVak1/d94G1uw4yLXHHU7zhunsPljM2p0HGdCxOQDTlm3jo4WbeTlqZvSYHq15\n+pLBBNQZN/fw+RV/bWas3MHgvBYRn9Wt+woY9vcv+O3orpx79KGhOIGAVro1L3iMi3wBfli3m1dn\nr+ONues5fWBO3MdMXryFb1ft4LaTe8fdpzx+2n6AzLQU2jdvUOp+M1ft4KjO2SzYsIcHv1jBvHW7\neej8Ixn29y/41bDO/Pm0voBzHG97ZyGn9GvPsK4tGXT35wzv2pKXrhga8f/W7jjI3z9azH9+MYD/\nfLqMZ75ZTdfWjfnF0R2rVJ/qVKmgTVX9InI+cF+Cy1MrlLX2aLg/ju9Jg4xUbj2pFxlua0BwUPaD\n5w0sc/xAg/TUcs1Iq6zPftxCsT/AVWO68emPW0pdWqsiwluugl0DV74wJ2H/P2jmTyXH31TU41+V\nTDNSmo/KaKEMtsCkpkhCWigen7aKjLSU0CSMaPsKfOS7rTazVu+KaDULtlABJVqTgjNWg+79dGmV\ny1pen8V4H4SnrwkGbAAfh80Q3X2wqMRFRfjM6qCj//p56O8X44zjVNXQ+ECAkfd8EQreVv/jZFQ1\n7lX3dDevX2nCWwT3FRTz8cLNnHVUh1Bw8fJ3a0t8yZfXkL9ODv0dfRzL49yw98LWvYXc5eZHHNip\nOd/F+Ewt2byXr1dsZ0S3VqFtWsbY3qCpS7fSuWUj8lo14usV22nTJJPD25acmBL8rKx3W7A27ykI\ndcmrKt+v3c0H8zexeU8Bb/5mOADb3VbO9NSSuSzfnLOe52asdvbbX8TyLU7LeLD7PDgeNTxoO1jk\n56bXfuDjRZt54bIhDO/aktdmrwvNeH582qrQBcSO/YWoKm/O3RCRkinclKXbIi6GDhT6aJRZ8mt3\n1uqdZKWlckSHyMlUU5Zs5ZJnZnHHhN6cMTCHkfdM4baTe4Vmb787byOn9G8f2r/A52dvvo/v1+5i\n/BHtSjzHdS9/z4fXjaR5VIaDmat2MMtt1d9f6AvV0R9QfvPCHB698KjQvsu37GPL3kJGdGsZ6oH4\n/Yk9+WjhJh6ZspKPrhtZZuAYvbbzGPc9nJWewmc3jKJji4ah+8K/S86Z9C0pcqi1em/+oVbZF2eu\nDQVte/N9vDRzLa/OWsfyu8cDhPKfhr9vb393IdOWbaNTy2XMXev0ggQvflSVBRv2UOwve+iIl6oy\npm26iDwEvAqEpvepap3P/FjekxU4J4g7T+0T877wD9s5gzqGciJNPLYLXy7dxrtXj+DzxVu4+qXv\nOTo3mxtP6MF5j8fuxqmK79fuDn3pJ8qT03+icWYqz85Yg4gzFsxXB7Pnn3lkDqf2bx8a0P1T2ExX\nf0DpmoB1OstqNbnvs2UsjpOfbV+c1rlYg+QfnlJ2q2SyPfjFCk4N+9zEE2sSSrRTHprOwg2HXrdg\nwBZ07SvzQgP5K2PHgSJ8/gCPTl3Jm3PXs3rHQfJaNeJgYdU/a4m8kNt+4NBrdcYjsbvJdx0s5oIn\nZoZmqh4o9NHnT59w+8m9uHxkFxZv2stXy7fxtw+X8JvRXfnjuJ6Ac668+OlZZKSmsOyv47ngiZmA\nExT3u/MTeh7WNDRmM9rMn3awxW1BzUpPDXXHz1mzixtfncex3Vvz7IzVgDMea1/UUJDwVvYnvlrF\n68Gu9qiW2tybPwirZ1FoSMKWvQV8sGBTRIqa8M/Tn//3I/nFfu75uPSLnYKwsYg79hfFDNrOfswZ\nErD6HyeHAtW5a3aFujzvev/HUGB9S1h59uYXR3yW84v8nPXYN6zflc+MW8ayv8DHJ4s2c9WYbjw1\n/Sc27ing44WbOXdwJ3YdKCIjLYU9+cURF3RvzFkfMVnno4Wb2VdQTJMsZ+jOCW4PUbi1Ow9y+9sL\n2VfoY9HGvaHg85GpK7jn46Us+cu40Eozz81YzR3vLmLeHSeUCB4LigOMvGcKPdo24ZMbnGQU0cNN\nAgoH3O8rv2qoWz38+2XdLmcWfYP01IhWaXC6qsf3bcfDFxwZ6jINb02//Z2FJSajJSILQ6JUJWgb\n4P6+K2ybAmOr8D9rhbKS61bGP8/qF/r71pN6cetJvQDo5F5xpKYIrRonPv8bOF1EsRKGVkV4d5Qq\n+GrwQOaqyG6YQbtmpTfpV7cfN1U8N9sZD5c9hq2mSsT6p1OWbo0I2KKFd89W1r8+WRpqzQm6491F\nlTpe1akis76DLSTB8a4PfrGCZVv28drsQ5/3R6euZNu+Qu75eb/QF2qRP8BVYSs5XPXiXPYW+EqM\nwV2381DKmvBhFB8v3Mw3Kw516UYPg8hIS4m7agwQCtjAGd4SqzURnJQ/6W76nx0HCsucjV1WwAaR\nF0jhM1tXbN1Ho8y0iPPHmh0HOOXB6WXOgg7aV+jjQNhFQH6xP9RSOezvh4YjnD2oY2i4yg/r93DO\n0VqhCW8fzN9En/bNeGTqipj3n/nI105PUqGTuqhDdgPuOasfj7sz4rfvL6RV40wWbdzDHe86aWc2\n7M4vEbQFLd2yj29WbqdjdkPmrI3fm+IPUCJYh0OttY0yU0usyR1QZ9Lcw0Tmm6wtKjumLQV4VFVf\nS3B5aoWy1h5NpD7tm/HrUV049+hOpLlRogh8cdPoUJNyIgRn1dV1n984iuP/U77ktVN+V/Zr3LZp\nZszZiuXVrlkWfXOaxewurE4VmU1YF0Xn6otW1YAtnpoWsEH5Z2sDDPnbZC4ZkRuaGb0nvzgiYAt6\nY856pi3bFlrpBJwv/tDfCw79Hd79NfKeKTGf15mIEX/c36Y9BbwVZ9xitHU780OTXaLd++nSUHni\nDUeoqKPuPtRVP+HB6RzbvTU9D2vCpGmraN8si+l/PNTOMepfUyv8/8Nf12P+Gfv1G3nPlFC9dh4o\nZNmWiq10UVZC7PAgc4Y7XvL2Cb1Dk/We/3YNC9bviViicdeBYmat3hl3WbbzH59ZZrlWbtvPh1Hv\npYy0FNa7LW2NMtMigvnwCXLTl29nYTmTkatqhXKQVqdKzR5V1QDwhwSXpdYoa+3RREpNEW4Z34u8\nVo1CiVqz0lLJa9WIZjFmmrZqnMmfTokcFLr4rnEl9qtPbj+5F09ffDTvXDUi5uSIeKJTl2Q3TOe7\n246L2HbSEe1o3rDkcWjbNH4qgGO7t2ZU99aA0706vu9h5S6TVzpkN2DBnT8rsT3YCgEw9XejPSyR\nqS43vhZ7PFYs2/cX8q9PlpZrIsrWcnRRA3SPSm9SGTsPFJWa17G8pi7dxq6DFc/zVhHTlm0L5WTc\nuKegxASJino1Ru7JaOGB8a6DxZVaCq2i+v/509BM6v/35aoSa2rvPFjE71//gS2lJK0uS5EvwPPf\nrgndDqZmCra0oZFrSIfnWrzwyZnlHnN8IMHDh6qiKik/PheR34lIRxFpEfxJWMlqsFCeNo8D72Ba\ng2Dw9sOffsbZYYNo+7Rvyuzbj+eSEXkRX7ilZeX3St+cpvx6VBeO6daKfh3KXrkgluFdW3LdcRXP\nQ3T5yC6M6dkmNNPrsxuO5fUrhzEkr0XEIGSApy6OnwcoKz2VNk2yePJXg3jg3AG8OnEoHbIbkp6a\nEjHO6slfDeKZSwZHPPa5SwcztqeTef74Xm3o75YlVaTEzNWMasroXxG92jWlSVZ66MJgeNeWvHf1\nCL679fjQPrnlzMdXHoNzS546Rh7eKsaeVTe6R+tq+b+J8tIVQ+jVzkkXdM6gjlwyIje5BTLVasjf\nJpe9UwLtPljE3gomQq8O1778PauruIJLtJ0Hi/h21Q7ened0ne86WBQRtF1fRuaCuP+3ChkVEq0q\nY9rOcX+H52ZToEsV/metEJqI4EH3aLjGmWlcO7YbJ/c7FCDcfUZfrhl7OGmpEtHyFhw0GmzRSUsR\nDm/bhAfPG8jCDXvYsDs/NN7mgiGd4s6wS5THLjyKDtnO+Lx7P1nK/PUVv9Jr1iCdk/u1C6UrCerW\npnGJnGmlCc5ce/XXw9hXUMwnizaHxkWM7dmWpy4eFOqKGdOjNVPchLzBgbSxlrr573kD6dehGXd/\nsJiBnbLJjmp9a5iRGuoqaJCeGmqmT0kRTh+Qw6KNe7hiZBcaZKSSmZrKH9+cHzFrMtyU342mVeMM\n9uQXh7pD7j27Pw9MXlZiMH1lnNK/PXdMcFprv/rjGA4W+kNLMqkqFw/P5fhyLPdTXjnNG/DalcPY\neaAIf0A5Z9IMAgHl+cuGRAwSv/WknrRpksX1r1buxBt0y/heMZMsl+X2k3vRu11Tzn+i7G6bWIbk\ntaB/x+ZlrnxxVOfs0FCIvjlNad+8AU+7Kz547eLhucxbtzvhY14TZfJNo/h00Rb++XFiujLrg10H\ni2OuGpGWUvsni427/6uI23vyi9kdFqBOXrK1Uv93+4FCOrVsWPaOHqj0Jb2q5sX4qfMBGxxK+eFR\nSrUQEeHGn/Wgx2GHpstnpqXSqWVD2jdvUGJW0tK7x/HUxUcD8ONd4/jf1U734OkDc/ht2Npy1bWm\nYItGGfRo24Qzj8wJBWwAV4/txm3uRAtwvgzLo1mD9FBrY9DRudl8dkPJJW+z0lN4/KJBfHDtMaX+\nzyZZ6Tx43sCIbWN7tg2tK/jEr45m3h0nAHBjVEb/aJcdk8fiu8bRolEGIsKSv4wLTU3PSk8l2ICW\nlZ4ayjCfliI0yEjl7tOPoHPLRrRpkkWzhuk8dP7AmN3a7ZtlkdeqEU2y0umQ3TDUajmhXzt+97Me\ngNMCW5kWyaD7zxkQyu/WNCs9Yg1NEeHOU/twTIJawV64bAhTfz8acN4vrZtk8vF1x/LpDaMi9lvy\nl3FcMbILvduXL2F1aWJ1kS/5yzjOPqoDx3SLXa8f7zqRy47Ji0hF8MC5A2jVOIM/ndK7xPsyliF5\nLbh6bDduP7kXZ0Tlvrr95F68cNkQnrt0MJlpqeS4+apSUqTU3I2TfnlU3PvC/eW02DPYy3LTz7pz\n5ajEnNaPzs3msQvLV95Y3vzNMM50X7dR3Vuz5C/j6Nq6cagFHeDPp/YpMTykvJ64aBBZ6c5x/PfZ\n/bmzkv+nIqb/cUy1P0e0bfsKSySoPuuoDiy660Tev6b082XQnaf0jtk6XtME1OmOrqqKLtFYnSrd\n0iYiF8XarqrPVb44tUOwlaSmDEyMJzxpb0bUl4qI8PsTe4Ru/9+E3rw+ex1LNld9gPpdp/XhhN5t\n486qzEpP5Ypju/BXd4H0y0d2icjRFXR0bnbECgLNGqSHWrsAXpk4lO5tmyAivPmb4ezYX8hzM9Zw\n/fGH0yG7YdwFu6O1bBR//FlqitC8YUa5pnyLSERXdFZ6Kr8anstjX66kaVZ6xJdvb7f7K9gNFi0t\nNYW0qF7t6EXpAZ69ZDDLtuwjKz2V0wbkcNoA50st1njHWC47Jo+CYn9ES6tXCZ6PyGlGr3ZNQsvn\nBIW/V+8+vS8rtu4PHfeW7qSPcX0O4/KReZzlpkqIfq+UJlb9stJT+dfZ/QEnJUD0WJeGGc6pMnhh\ndN7gThGv9yUj8thzsJjFm/dG5EALl9e6EU2z0rl8ZJfQEl6XH5PHBUM7lxg/2b9jcz5etJktewsj\nLnjAST8QbIEc2CmbjLQUinwB7j69L1OWbC3RmtDzsCb8clguAzpm06JxBk9+9RNfr9herskoTbLS\nOebw1gzJa1GpnIhH52YzoGNzpi3bzrlHd6p0sH/JiFyO6tyC5g0zeOv7DQzs1Dz0nggG8uce3ZFf\nDc8FnHQc5TG8a8uIsVaju7fh40WbSUsVzhjYgU9/3MIROc04UOTjylFdWbvjII9+uZKvljt5+l64\nbAib9uTz3g8bQ9temTg07nsAnGA/ODM2eGyz0lNCqUHCj291eW7GoXFg14ztxk3uBV/fnGalPv/Z\nR3Vg7c6DDM5ryddRY9TAuai8/9yBPDdjNe+HTZCIdvIR7SImoyTacT3bcPagjlz78vcVmm1+5aiu\nDOqczeVu8u9VfzsJxbtzYnlUpXv06LC/s4DjgLlAnQ/aQsl1k1qKqrtqTLfQ35cdk8dxPdsw+t6p\nnD6gPY2z0kJL61wxMo9ju7dm0rRVoRNTaYZ3bVWuNBhzbj+e6Mb48Nmdr185POLk0TSqpW1ol5ah\nv4PJD3/Wp+KD+hPRehPPH07swVlH5dCpZUNS3UXF/QHl9IE59GrXJJRouTwaxhibmN0ogyFhr0PQ\nqB6t+dmqtozs3pr+HZoxb93u0FT7cH1zmrJ9X+XHa7RqnFnh5c6+u/U40lNTyC7HrNsLh3aOuN2y\ncSZTfjeaDtkNIoK9Fy8fyvgHprFy2wF+O7oro7q3pnWTTMb+O3KmcHAsYf8Ozfhh/R4eu/CoEuuc\nzr7teN6cuz7mhUSLRhmh54/WrGE6Q7u0ZPofx9AkK52dB4qYvHgLizbu5e3vNzCo86HnCda9f8fm\nMdfqvXh4Lht2H+SiYZ1DuaRO6d+eozo1j9ivcWYa5w/uxDPfrOb8wZ24cGhnHvpiOTN/2slXy7cz\n+aZRtG3qXLwEc2fdcUpvVJWV2/bz/drdPPHVT6UGcI0z03j118MiPotDu7Tg3rP78/TXqznpiHa0\nbpzJsf+awklHHMaHCw516z964VG0apzJbWVc85w/pFNoVYZoj1xwJCe5iWK7tm7M5JtG0TmsxbNZ\ng3Sm/m50xEVaeqrEXNkh2kPnH8nMVTv4zYtz6dCiAX8/8wiyG6UzukcbmjVIL5FBv0N2Q4Z318cm\n3QAAIABJREFUa8W+gmL2FfhCGfxPH5jD4bd9xOXH5DG0S0tymjdgw+583r/mGCY8eGjZuNE9WnNq\n//YR6Uxmup+H4//zZYkcctG6tG7Eqm0H6NamMT0Oa8IH8zfRtmkmD553ZNwZsb8a1pln3SAtPDhM\nTRH+d/Ux9GpX8hz03/MGEggofXOa8vs35rNg/R4+um4k3do0DjVWHOa+r246oTsn9GnLuPu/4oKh\nnRmc14LBeS1YsvlLVmzdz40ndGdPfjFPTneSmH9+4yi6tWnM2DnruSlGUuK3fzs8lC/w0xuODa0c\n9MuhnZm1emfchoUrRuZxSv/2NMlKD32mhnRpUeZ3VvDC0FkNoREtw1JrVcdasVVV6aBNVa8Jvy0i\nzYFXqlyiWqC8a4/WNrmtGrH4rnE0yEhlf6GP12av5z+/6M8Edwxdq8aZjH/gq7iPTxG467S+5Z6h\nGWux5ejHXjEyjzfmrKd1k0z6dWgWurruk8BAKzVFeHXiUAqqYaHklBQJBWbBiZfBVpyKBGwADTPL\nP6Gka+vGTApbXLlfh+a8/f0Gvl+7m9+O7soL365hb4GPrLRU2pQy07Usb/92OFOWbo0ZEAJce9zh\njO97GLktG9Hrjo8BaNO0ast8hQc5N4/vyX8nLycjLYXnLxvC9BXb+cWgyCVoMtJS+OGOn5GaIqHJ\nQ89dOoSXZ63luF5tSrT0ZTfK4JIReQAxA7dYQVa4YOtJswZOq5qq8tcz+oZa6wBuOL47bZpkhoKR\naMEu86DFd40jKz2lROt+VnoKd0zozc3je4a+YK4eezhXqVJQHIg7CUnEeV92a9OEU/q35/kZa0It\n34vvGsePm/bGXTbtrtP6cPIR7WjZOJP/m3Co5XfJX8aRkZpC36WfcKHbetgqxmd8yV/G8fcPF4cC\nCYC/nXEEG3fnh8YaXjWmK9cd1x1/QEvUIbiEVbjoSTHpqSkU+2PP+PvTKb1DLXEtGmUw/oh2/HjX\niaHj8/cz+8V8XLgmWemhccPB5wvWH+D1K4fx1fJt9M1pRt+cpizcsJfvbj2Olo0zSxzDYFD97S3H\nxZzc1qNtEzLTU3hl4lDem7eRm99awICOzZl4bBc+mL+JRplpHB513rxyVFfaNs3knKM70iA9lYuG\n5/LqrHWcNqA9+wt8nDPpW/wBjXvBGj6x6vVfD8OvGtFzA04+0ZzsBpw+IIfDmmWF3qNBwffP6B6t\n6dehOUd2yqZ728ahc/zPj+rA+CMO49JnZtG+eQN+2n6AozplM7DTodUHuoetnPGX0/sy8bnZrNp2\ngCJ/gP/3y6P4esV2npuxhhaNMmIupTW2Z5syg7Y+7ZtyxsAcOmQ34MwjO5QrMXcyVaWlLdoBIC+B\n/6/GCoSmjya3HNUheIJsnJnGMnf5j6Doxexf+/Uwzp00g4DC0xcfTZ+cprRpUrkv5C9/PzqU5T27\nYXpo2v1tJ/cu8WF86fIh9Mmp3AzUeGK1ViXatccdzopt+zmuV5tyP+brm8cyZ80urn35exqmV20W\ncFP3S+aoztms3+V06eQX+zljYA4NM9JYtmVfhcdudGzRkAuGdA4FbT/9/aTQYvaDOmdz1ZiuoZP9\nxGO7hK7OE+XKUV25cpQzPrN98wYlArZ3rxpBm6aZJb74mzVMDz0ultQU4fKRXUhLEVZsq1hOq2gi\nEhGwgfM5u3xk+ceKRZe/X4dmzF+/BxEnEM1Kibw/uqu+NMHhCmN6tg49V6xle4KtRxcNy437f8AZ\nP1vW890+oTd9c5pFrKf6zCWDeXTqSv758RJSRUoM6aiIn/VuyztRKUD+b0Jvxvc9jPbNG/Dn//1I\no7DXJ/r4VEb40I32zRtwjrsm6CsThzFr9c6Ii5ULh3aiR9RSXrHqGz0so7HbEneg0BdKNXTs4a1L\nHOurx3ajcdgY566tG4cStld09Zu01JSYgUKDjNSIz1B0GYJfj8Fel5P7lbxAaZiRxisTY6+IEfTH\ncT1D56Xje7Ulr3Ujbhnv1OXITtk8N2MNJ8e5+Bnbs02JrvKhXVrw7Sqnq/+piweFgsTg57G6ktgn\nSlXGtP2PsDH5QG+gXiXbrYEtp9UqOJ7o0hF5XDIil44tGvLRdcfy7rwNjO7Rukpj/Dq3PHSlPPX3\nY0o9sQyPM1i8puvSujHvXzOyQo/Jad6Aw5pmMW/tbsb0LH+wF8sZA3P4ctk2OrdsRMcWzolU1fmC\nP6F3W07oXbkZoakpwg3Hd2dMz8j3wBvu+pBBt55UvgknidS/Y/OydyrFxSNq5nXoy1cMZdfBxKYh\nKKvl96PrR8Ztfauo9NQUzh7UMSJoA7hoWGfW7TrI5cdWbfLDP8/qx/XHd+fDhZtCqxYM79oy1J35\n7S3HhSYJVbfGmWmM6RH52Q1vRY1lzu3HE2siZ3D8bQM3/dAXNznrdKanpjD1d6N5ZOoKXpu9PiIg\njRYMrqorpU7QkxcfzXs/bKxUEHT/OQNCf/8mbNJc9ELurZtk8vmNo8iNM7Ozc8tGXDO2G0PyWnLh\nk86s71cmOl39I7q1ZGzPkue8tBqQcqk0opVcXkhEwqd3+YA1qlq+lNTVbNCgQTp79uxq+/+Tpq3k\nbx8uYeGfT4y4mqkPCor9ZKSm1Mi+flO2g0U+GmakUVDs57kZq7l0RF7CT1Lvz99Iz8MqNl7P1E+n\nP/w1XVo14j9hX9KJdtVLc/lg/ia++sOYiNm/tZGq8sw3qzljYE7MJaACAaXIH4ho9YuloNhPWorU\n+AAlkYJjMlf/42QKfX5SJX79w/f1iojMUdX4iUJdFY44RKQb0FZVv4zaPkJEMlW15q86XUV1uHe0\nTGWdDEzNFuwGykpPZeKx8bsHq2JCv7IXdDcG4J2rRlT7c9zz8378YlDHWh+wgdMqfkkprb8pKVKi\nqzyW+ngen9CvXWhoSvT4vGif3nBsje1Jq0wz0f3ALTG273XvO6VKJaoFDuVpq6FH1RhjDOCkaQkm\nGTf110PnH1nufbu3rbm9BJVpG22rqiVWj3W35Va5RLVAcCKCxWzGGGOM8UplgrbSRvaWnZyrDkjW\n2qPGGGOMqb8qE7TNFpErojeKyOXAnKoXqSQRGSciS0VkhYjcXB3PURHJWnvUGGOMMfVXZca0XQ+8\nLSIXcChIGwRkAGckqmBBIpIKPAycAKwHZonIe6pavnVKqoG1tBljjDHGaxUO2lR1CzBcRMYAfd3N\nH6jqFwkt2SGDgRWqugpARF4BTgOSF7S5v20igjHGGGO8UpVlrKYAUxJYlnhygHVht9cDQ6J3EpGJ\nwESATp06VWuBQhMRqvVZjDHGGGMOqTOZYVV1EjAJnOS61flcFw3LZUK/dtY9aowxxhjP1IagbQMQ\nvnZFB3db0rRolEGLRjV7fTJjjDHG1C21YQ2LWcDhIpInIhnAucB7SS6TMcYYY4ynanxLm6r6RORq\n4BMgFXhKVRcluVjGGGOMMZ6q9ILxNZmIbAPWVPPTtAK2V/Nz1GRW//pb//pcd7D6W/3rb/3rc92h\neuvfWVXLXG+tTgZtXhCR2ao6KNnlSBarf/2tf32uO1j9rf71t/71ue5QM+pfG8a0GWOMMcbUexa0\nGWOMMcbUAha0Vd6kZBcgyaz+9Vd9rjtY/a3+9Vd9rjvUgPrbmDZjjDHGmFrAWtqMMcYYY2oBC9qM\nMcYYY2oBC9oqQUTGichSEVkhIjcnuzzVSUQ6isgUEflRRBaJyHXu9hYi8pmILHd/Zye7rNVJRFJF\n5HsRed+9nSciM933wKvuah11kog0F5E3RGSJiCwWkWH16fiLyA3ue3+hiLwsIll1+fiLyFMislVE\nFoZti3m8xfFf93WYLyJHJq/kVRen7v9y3/vzReRtEWkedt8tbt2XisiJySl14sSqf9h9N4mIikgr\n93adOvYQv/4ico37HlgkIveEbff8+FvQVkEikgo8DIwHegPniUjv5JaqWvmAm1S1NzAUuMqt783A\nZFU9HJjs3q7LrgMWh93+J3CfqnYDdgGXJaVU3ngA+FhVewL9cV6HenH8RSQHuBYYpKp9cVZlOZe6\nffyfAcZFbYt3vMcDh7s/E4FHPSpjdXmGknX/DOirqv2AZcAtAO558Fygj/uYR9zvh9rsGUrWHxHp\nCPwMWBu2ua4de4hRfxEZA5wG9FfVPsC97vakHH8L2ipuMLBCVVepahHwCs4BrZNUdZOqznX/3ofz\nhZ2DU+dn3d2eBU5PTgmrn4h0AE4GnnBvCzAWeMPdpc7WX0SaAccCTwKoapGq7qYeHX+c5f4aiEga\n0BDYRB0+/qo6DdgZtTne8T4NeE4d3wLNRaSdNyVNvFh1V9VPVdXn3vwW6OD+fRrwiqoWqupPwAqc\n74daK86xB7gP+AMQPnOxTh17iFv/3wD/UNVCd5+t7vakHH8L2iouB1gXdnu9u63OE5FcYCAwE2ir\nqpvcuzYDbZNULC/cj3PCCri3WwK7w07kdfk9kAdsA552u4efEJFG1JPjr6obcK6s1+IEa3uAOdSf\n4x8U73jXt/PhpcBH7t/1ou4ichqwQVV/iLqrXtQf6A6MdIdDfCkiR7vbk1J/C9pMuYhIY+BN4HpV\n3Rt+nzp5Y+pk7hgRmQBsVdU5yS5LkqQBRwKPqupA4ABRXaF1/Phn41xR5wHtgUbE6D6qT+ry8S6N\niNyGM1zkxWSXxSsi0hC4Fbgj2WVJojSgBc7woN8Dr7m9LUlhQVvFbQA6ht3u4G6rs0QkHSdge1FV\n33I3bwk2hbu/t8Z7fC03AjhVRFbjdIWPxRnj1dztLoO6/R5YD6xX1Znu7Tdwgrj6cvyPB35S1W2q\nWgy8hfOeqC/HPyje8a4X50MRuRiYAFygh5Kb1oe6d8W5YPnBPQd2AOaKyGHUj/qDcw58y+0G/g6n\nx6UVSaq/BW0VNws43J09loEzEPG9JJep2rhXFE8Ci1X1P2F3vQf8yv37V8C7XpfNC6p6i6p2UNVc\nnGP9hapeAEwBznJ3q8v13wysE5Ee7qbjgB+pJ8cfp1t0qIg0dD8LwfrXi+MfJt7xfg+4yJ1JOBTY\nE9aNWieIyDic4RGnqurBsLveA84VkUwRycMZkP9dMspYXVR1gaq2UdVc9xy4HjjSPS/U+WPvegcY\nAyAi3YEMYDvJOv6qaj8V/AFOwplFtBK4Ldnlqea6HoPTFTIfmOf+nIQzrmsysBz4HGiR7LJ68FqM\nBt53/+7ifkBXAK8DmckuXzXWewAw230PvANk16fjD/wZWAIsBJ4HMuvy8Qdexhm/V4zzJX1ZvOMN\nCM5s+pXAApxZtkmvQ4LrvgJn7FLw/PdY2P63uXVfCoxPdvmro/5R968GWtXFY1/K8c8AXnA//3OB\nsck8/raMlTHGGGNMLWDdo8YYY4wxtYAFbcYYY4wxtYAFbcYYY4wxtYAFbcYYY4wxtYAFbcYYY4wx\ntYAFbcYYY4wxtUBa2bvUPq1atdLc3NxkF8MYY4wxpkxz5szZrqqty9qv1gRtIpKKk+Bzg6pOKG3f\n3NxcZs+e7U3BjDHGGGOqQETWlGe/2tQ9eh2wONmFMMYYY4xJBk+DNhFpJCIp7t/dReRUdzHysh7X\nATgZeKK6y2iMMcYYUxN53T06DRgpItnApziLr58DXFDG4+7HWbC3SfUWzxhjjDE1yeuz1/HizLVJ\ne/7/m9CbozpnJ+35w3kdtImqHhSRy4BHVPUeEZlX6gNEJgBbVXWOiIwuZb+JwESATp06JbLMxhhj\njEmSjxduZsXW/RyZpMApPVWS8ryxeB60icgwnJa1y9xtqWU8ZgRwqoicBGQBTUXkBVW9MHwnVZ0E\nTAIYNGiQJrbYxhhjjEmG4oDStU1jnrt0cLKLknReT0S4HrgFeFtVF4lIF2BKaQ9Q1VtUtYOq5gLn\nAl9EB2zGGGOMqZt8/gDpKTWntSuZPG1pU9UvgS9FpKF7exVwrZdlMMYYY0zt4fMraTWoizKZvJ49\nOkxEfgSWuLf7i8gj5X28qk4tK0ebMcYYY+qO4kCA9NTalKGs+nj9KtwPnAjsAFDVH4BjPS6DMcYY\nY2oJn18taHN5/iqo6rqoTX6vy2CMMcaY2qHYHyDNxrQB3s8eXSciwwF1k+raKgfGGGOMiavYb92j\nQV6/ClcCVwE5wAZggHvbGGOMMaYEX8AmIgR5PXt0O2WvfmCMMcYYA7izR1OspQ08CtpE5EEgbsJb\nVbW0H8YYY4wpweketZY28K6lbbZHz2OMMcaYOsS6Rw/xJGhT1We9eB5jjDHG1C02EeEQr5PrfiYi\nzcNuZ4vIJ16WwRhjjDG1hwVth3j9KrRW1d3BG6q6C2jjcRmMMcYYU0s4ExGsexS8D9r8ItIpeENE\nOlPKBAVjjDHG1F+q6o5ps5Y28D657m3AdBH5EhBgJDDR4zIYY4wxphbwBZx2nXRraQO8z9P2sYgc\nCQx1N13v5m4zxhhjjIng8ztBm7W0OTx5FUSkp/v7SKATsNH96eRuM8YYY4yJUBwIAFieNpdXLW03\n4nSD/jvGfQqM9agcxhhjjKklgi1tNnvU4VWetuC4tfGqWhB+n4hkeVEGY4wxxtQuxX6npc2S6zq8\nDl2/Kec2Y4wxxtRzwaAt3dYeBbxbe/QwIAdoICIDcWaOAjQFGpbx2CxgGpCJU943VPVP1VhcY4wx\nxtQAhyYiWEsbeDem7UTgYqADzri24Ku/F7i1jMcWAmNVdb+IpOOkDPlIVb+trsIaY4wxJvl8gWD3\nqLW0gYdrj4rI88B5qvpiBR+rwH73Zrr7Ywl5jTG1wuJNe1m+dX/ZOxpjSti4Ox+wPG1BnuVpU9WA\niNwAVChoAxCRVGAO0A14WFVnxthnIm6i3k6dOkXfbYwxSTHx+dms25mf7GIYU6u1apKZ7CLUCF6v\niPC5iPwOeBU4ENyoqjtLe5Cq+oEB7mLzb4tIX1VdGLXPJGASwKBBg6wlzhhTIxwo9HNK//Zcd9zh\nyS6KMbVSg4xUcpo3SHYxagSvg7Zz3N9XhW1ToEt5Hqyqu0VkCjAOWFjW/sYYk2zF/gAtG2XQrU3j\nZBfFGFPLeb2MVV5FHyMirYFiN2BrAJwA/DPhhTPGmGrg86tlczfGJITXLW2ISF+gNxBKqquqz5Xy\nkHbAs+64thTgNVV9v3pLaYwxieELBGzmmzEmITwN2kTkT8BonKDtQ2A8MB2IG7Sp6nxgoBflM8aY\nRFJViv1qM9+MMQnh9eXfWcBxwGZVvQToDzTzuAzGGOMJfyCYGNRa2owxVef1mSRfVQOAT0SaAluB\njh6XwRhjPOELWDZ3Y0zieD2mbbabtuNxnLxr+4EZHpfBGGM8EVw3McNa2owxCeD17NHfun8+JiIf\nA03dMWvGGFPnhNZNtDFtxpgESMbs0TOBY3Dys00HLGgzxtRJwZY2G9NmjEkET88kIvIIcCWwACc5\n7q9F5GEvy2CMMV4pdse0WZ42Y0wieN3SNhbo5S4Cj4g8CyzyuAzGGOMJX7ClLcVa2owxVef1mWQF\nEL6ae0d3mzHG1DnFfps9aoxJHK9b2poAi0XkO5wxbYNxZpS+B6Cqp3pcHmOMqTa+gNPSlm5j2owx\nCeB10HaHx89njDFJY7NHjTGJ5HXKjy+9fD5jjEmm4OxRa2kzxiSCnUmMMaaa2IoIxphEsqDNGGOq\nibW0GWMSyc4kxhhTTYJj2ixPmzEmETwZ0yYiC3Bmi8akqv28KIcxxnip2PK0GWMSyKuJCBPc31e5\nv593f1/g0fMbY4znLE+bMSaRPAnaVHUNgIicoKoDw+66WUTmAjfHe6yIdASeA9ritNZNUtUHqrO8\nxhiTCJanzRiTSF6fSURERoTdGF6OMviAm1S1NzAUuEpEeldjGY0xJiEsT5sxJpG8Tq57GfCUiDRz\nb+8GLi3tAaq6Cdjk/r1PRBYDOcCP1VlQY4ypKps9aoxJJK+T684B+geDNlXdU5HHi0guMBCYmfDC\nGWNMBX20YBN3vLcI1djzrAqKLWgzxiSOp0GbiLQF/ga0V9XxbjfnMFV9shyPbQy8CVyvqntj3D8R\nmAjQqVOn6LuNMSbh5q3bzc4DRZx7dMe4+7RpkkXbppkelsoYU1d53T36DPA0cJt7exnwKlBq0CYi\n6TgB24uq+lasfVR1EjAJYNCgQXHTixhjTKIU+5WG6an89Ywjkl0UY0w94HWbfStVfQ0IAKiqD/CX\n9gAREZygbrGq/qf6i2iMMeVT7A9YOg9jjGe8DtoOiEhL3ES7IjIUKGtc2wjgl8BYEZnn/pxUzeU0\nxpgy+QIB0my8mjHGI153j94IvAd0FZGvgdbAWaU9QFWnA3Ypa4ypcYr9Srql8zDGeMTr2aNzRWQU\n0AMnEFuqqsVelsEYYxLF57eWNmOMd7yePXpm1KbuIrIHWKCqW70sizHGVFVxQG1MmzHGM8lIrjsM\nmOLeHg3MAfJE5C5VfT7eA40xpqbx+QOk22LwxhiPeB20pQG9VHULhPK2PQcMAaZxaCF5Y4yp8Xx+\na2kzxnjH60vEjsGAzbXV3bYTsLFtxphaxeketZY2Y4w3vG5pmyoi7wOvu7d/7m5rhLMOqTHG1Bo+\nf4AMa2kzxnjE66DtKuBM4Bj39nPAm+os3DfG47IYY0yV+PxKmo1pM8Z4xLOgTURSgc9VdQzOklTG\nGFOrFfkDNEn3+trXGFNfeXaJqKp+ICAizbx6TmOMqU6+QIB0G9NmjPGI15eI+4EFIvIZcCC4UVWv\n9bgcxhhTZU73qI1pM8Z4w+ug7S33xxhjar1iv7W0GWO84/UyVs+KSAOgk6ou9fK5jTEm0Xy2IoIx\nxkOeXiKKyCnAPOBj9/YAEXnPyzIYY0yi2OxRY4yXvD7b3AkMxs3JpqrzgC4el8EYYxLC6R61ljZj\njDe8DtqKVXVP1LaAx2UwxpiE8AXUxrQZYzzj9USERSJyPpAqIocD1wLfeFwGY4xJiGJ/wMa0GWM8\n4/Ul4jVAH6AQeBnYC1zvcRmMMSYhbPaoMcZLXs8ePQjc5v6Um4g8BUwAtqpq3+oomzHGVJTlaTPG\neMnr2aODROQtEZkrIvODP+V46DPAuGounjHGlJuquik/rKXNGOMNr8e0vQj8HlhABSYgqOo0Ecmt\npjIZYzxWUOzn5e/WcrDIn+yiVJqqApBuLW3GGI94HbRtU9VqycsmIhOBiQCdOnWqjqcwxiTIzJ92\n8uf//ZjsYlSZCOS1bpTsYhhj6gmvg7Y/icgTwGScyQgAqGqVl7ZS1UnAJIBBgwZpVf+fMab6FBQ7\nLWxv/3Y4vds3TXJpKk8QMtKse9QY4w2vg7ZLgJ5AOoe6RxVbj9SYesXnd66rGmWmkZmWmuTSGGNM\n7eB10Ha0qvbw+DmNMTWML+Bcs9nMS2OMKT+v2/W/EZHeFX2QiLwMzAB6iMh6Ebks8UUzxnilyOcE\nbZbjzBhjys/rlrahwDwR+QlnTJsAqqr9SnuQqp7nReGMMd7wBZzuUVtNwBhjys/roM1yrRlj8PmD\n3aPW0maMMeXlSdAmInOA6cBHwFRVLfDieY0xNVOxOxEh3VrajDGm3Ly6zB0CvA2MBr4UkQ9F5DoR\n6e7R8xtjapDQRAQb02aMMeXmSUubqvqAqe4PItIep6v0bhHpCsxU1d96URZjTPIFW9ps9qgxxpSf\n12uPng2gqhtV9SlV/QXwT5zlrYwx9YQv1D1qLW3GGFNeXp8xb4mx7WZV/drjchhjksgXCCACqdbS\nZowx5ebVRITxwElAjoj8N+yupoDPizIYY2qOYr9aK5sxxlSQVyk/NgKzgVOBOWHb9wE3eFQGY0wN\n4fMHSLdWNmOMqRCvJiL8APwgIi+pajGAiGQDHVV1lxdlMMbUHMX+gM0cNcaYCvL6rPmZiDQVkRbA\nXOBxEbnP4zIYY5KsOKCWo80YYyrI66CtmaruBc4EnlPVIcBxHpfBGJNkPn/AVkMwxpgK8vqsmSYi\n7YBfAO97/NzGmBrC51dbd9QYYyrI66DtLuATYKWqzhKRLsByj8tgjEkyp3vUWtqMMaYiPF0wXlVf\nB14Pu70K+LmXZTDGJJ/TPWotbcYYUxFer4jQQUTeFpGt7s+bItLByzIYY5Kv2K82e9QYYyrI67Pm\n08B7QHv353/uNmNMPeILBGz2qDHGVJDXQVtrVX1aVX3uzzNA67IeJCLjRGSpiKwQkZurv5jGmOrk\nsxURjDGmwrw+a+4QkQtFJNX9uRDYUdoDRCQVeBgYD/QGzhOR3h6U1RhTTYptTJsxxlSYpxMRgEuB\nB4H7AAW+AS4u4zGDgRXupAVE5BXgNODH6itm6Yr9AXx+TdbTG1PrFfoCNM70+vRjjDG1m9ezR9fg\nrD8aIiLXA/eX8rAcYF3Y7fXAkMSXrvwembKS+z5flswiGFPrHd+rTbKLYIwxtUpNuNS9kdKDtnIR\nkYnARIBOnTpV9d+VakS3lmSm96zW5zCmrhvVvczhrMYYY8LUhKCtrIEtG4COYbc7uNsiqOokYBLA\noEGDqrXvclBuCwbltqjOpzDGGGOMiVATpm+VFWDNAg4XkTwRyQDOxUkbYowxxhhTb3jS0iYi+4gd\nnAnQoLTHqqpPRK7GWf4qFXhKVRclvpTGGGOMMTWXJ0Gbqjap4uM/BD5MUHGMMcYYY2odUa17qStE\nZBuwppqfphWwvZqfoyaz+tff+tfnuoPV3+pff+tfn+sO1Vv/zqpa9mIDdTFo84KIzFbVQckuR7JY\n/etv/etz3cHqb/Wvv/Wvz3WHmlH/mjARwRhjjDHGlMGCNmOMMcaYWsCCtsqblOwCJJnVv/6qz3UH\nq7/Vv/6qz3WHGlB/G9NmjDHGGFMLWEubMcYYY0wtYEFbJYjIOBFZKiIrROTmZJenOolIRxGZIiI/\nisgiEbnO3d5CRD4TkeXu7+xkl7U6iUiqiHwvIu+7t/NEZKb7HnjVXa2jThKR5iLyhogsEZHFIjKs\nPh1/EbnBfe8vFJGXRSSrLh9/EXlKRLaKyMKwbTGPtzj+674O80XkyOSVvOri1P1f7nsh+wphAAAH\neklEQVR/voi8LSLNw+67xa37UhE5MTmlTpxY9Q+77yYRURFp5d6uU8ce4tdfRK5x3wOLROSesO2e\nH38L2ipIRFKBh4HxQG/gPBHpndxSVSsfcJOq9gaGAle59b0ZmKyqhwOT3dt12XXA4rDb/wTuU9Vu\nwC7gsqSUyhsPAB+rak+gP87rUC+Ov4jkANcCg1S1L86qLOdSt4//M8C4qG3xjvd44HD3ZyLwqEdl\nrC7PULLunwF9VbUfsAy4BcA9D54L9HEf84j7/VCbPUPJ+iMiHYGfAWvDNte1Yw8x6i8iY4DTgP6q\n2ge4192elONvQVvFDQZWqOoqVS0CXsE5oHWSqm5S1bnu3/twvrBzcOr8rLvbs8DpySlh9RORDsDJ\nwBPubQHGAm+4u9TZ+otIM+BY4EkAVS1S1d3Uo+OPs3JMAxFJAxoCm6jDx19VpwE7ozbHO96nAc+p\n41uguYi086akiRer7qr6qar63JvfAh3cv08DXlHVQlX9CViB8/1Qa8U59gD3AX8gcjnKOnXsIW79\nfwP8Q1UL3X22utuTcvwtaKu4HGBd2O317rY6T0RygYHATKCtqm5y79oMtE1SsbxwP84JK+Debgns\nDjuR1+X3QB6wDXja7R5+QkQaUU+Ov6puwLmyXosTrO0B5lB/jn9QvONd386HlwIfuX/Xi7qLyGnA\nBlX9IequelF/oDsw0h0O8aWIHO1uT0r9LWgz5SIijYE3getVdW/4fepMQa6T05BFZAKwVVXnJLss\nSZIGHAk8qqoDgQNEdYXW8eOfjXNFnQe0BxoRo/uoPqnLx7s0InIbznCRF5NdFq+ISEPgVuCOZJcl\nidKAFjjDg34PvOb2tiSFBW0VtwHoGHa7g7utzhKRdJyA7UVVfcvdvCXYFO7+3hrv8bXcCOBUEVmN\n0xU+FmeMV3O3uwzq9ntgPbBeVWe6t9/ACeLqy/E/HvhJVbepajHwFs57or4c/6B4x7tenA9F5GJg\nAnCBHsqTVR/q3hXnguUH9xzYAZgrIodRP+oPzjnwLbcb+DucHpdWJKn+FrRV3CzgcHf2WAbOQMT3\nklymauNeUTwJLFbV/4Td9R7wK/fvXwHvel02L6jqLaraQVVzcY71F6p6ATAFOMvdrS7XfzOwTkR6\nuJuOA36knhx/nG7RoSLS0P0sBOtfL45/mHjH+z3gIncm4VBgT1g3ap0gIuNwhkecqqoHw+56DzhX\nRDJFJA9nQP53yShjdVHVBaraRlVz3XPgeuBI97xQ54+96x1gDICIdAcycBaNT87xV1X7qeAPcBLO\nLKKVwG3JLk811/UYnK6Q+cA89+cknHFdk4HlwOdAi2SX1YPXYjTwvvt3F/cDugJ4HchMdvmqsd4D\ngNnue+AdILs+HX/gz8ASYCHwPJBZl48/8DLO+L1inC/py+Idb0BwZtOvBBbgzLJNeh0SXPcVOGOX\ngue/x8L2v82t+1JgfLLLXx31j7p/NdCqLh77Uo5/BvCC+/mfC4xN5vG3FRGMMcYYY2oB6x41xhhj\njKkFLGgzxhhjjKkFLGgzxhhjjKkFLGgzxhhjjKkFLGgzxhhjjKkFLGgzxhhjjKkFLGgzxtQaItJS\nROa5P5tFZEPY7W+q4fkuFpFtIvJEBR83WkTej3PfhyLS3P35bTn+1xQR2S8igypSBmNM3ZNW9i7G\nGFMzqOoOnGS/iMidwH5Vvbean/ZVVb26vDuHLW8Vk6qe5O6XC/wWeKSM/ceIyNTyPr8xpu6yljZj\nTJ0gIvvd36NF5EsReVdEVonIP0TkAhH5TkQWiEhXd7/WIvKmiMxyf0aU4zmyRORp9/98LyLB5W0u\nFpH3ROQLnJUDAJqKyAcislREHhORFHff1SLSCvgH0NVtJfyXiLQTkWnu7YUiMrI6XidjTO1lLW3G\nmLqoP9AL2AmsAp5Q1cEich1wDXA98ABwn6pOF5FOwCfuY0pzFaCqeoSI9AQ+ddcjBDgS6KeqO0Vk\nNDAY6A2sAT4GzgTeCPtfNwN9VTXYcngT8Imq/lVEUoGGVXsJjDF1jQVtxpi6aJa6i1eLyErgU3f7\nAtzFn4Hjgd7OOvCA0zLWWFX3l/J/jwEeBFDVJSKyBggGbZ+p6s6wfb9T1VVuGV52HxsetJUoM/CU\niKQD76jqvHLU0xhTj1j3qDGmLioM+zsQdjvAoYvVFGCoqg5wf3LKCNjKciDqdvTCzqUu9Kyq04Bj\ngQ3AMyJyURXKYoypgyxoM8bUV5/idJUCICIDyvGYr4AL3P27A52ApXH2HSwiee5YtnOA6VH37wOa\nhD1/Z2CLqj4OPIHT3WqMMSHWPWqMqa+uBR4Wkfk458JpwJVlPOYR4FERWQD4gItVtTCsizXcLOAh\noBswBXg7/E5V3SEiX4vIQuAjYCHwexEpBvYD1tJmjIkgqqW22BtjTL0lIhcDgyqS8qOayjEV+J2q\nzk5mOYwxyWXdo8YYE18+ML6iyXUTSUSmAF2A4mSVwRhTM1hLmzHGGGNMLWAtbcYYY4wxtYAFbcYY\nY4wxtYAFbcYYY4wxtYAFbcYYY4wxtYAFbcYYY4wxtcD/B9Oa8W5jkeJcAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(10,7))\n", "ax = plt.subplot(311)\n", "plt.yscale(\"log\")\n", "plt.plot(times/(2.*np.pi), errors);\n", "ax.set_ylabel(\"Relative energy error\")\n", "ax = plt.subplot(312)\n", "ax.set_ylabel(\"Current encounters\")\n", "plt.plot(times/(2.*np.pi), encounterN);\n", "ax = plt.subplot(313)\n", "ax.set_ylabel(\"Lost/merged particles\")\n", "ax.set_xlabel(\"Time [orbits]\")\n", "plt.plot(times/(2.*np.pi), -(totalN-N_pl-2));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us also plot the final positions of all particles. We can see that the planet stirred up the planetesimal disk." ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9wXNWVJ/DvUbuN2wwV2bGH4AZhIIwYvArWoMLOOLUV\nyA8TWIzGBAzBtclWsmxql5oikxIlLywYypQ1cW0Ns7XZzRA2tczCgI3j9IgxWUNiUlPxxmzktI3i\nBG/4aWiT8MPIM0QCtaSzf3Q/8dR6r/v9uK/7vX7fT5XK3a3X3Vft7tP33XvuuaKqICKidOlodQOI\niKj5GPyJiFKIwZ+IKIUY/ImIUojBn4gohRj8iYhSiMGfiCiFGPyJiFKIwZ+IKIUWtLoBbpYtW6Yr\nV65sdTOIiBLl0KFDb6vq8kbHxTb4r1y5EiMjI61uBhFRoojIq16O47APEVEKMfgTEaUQgz8RUQox\n+BMRpRCDPxFRCsU224eo3RSKJezYdwwnxiawojOHgfXd6O/Nt7pZlFIM/kRNUCiWsGXPKCbK0wCA\n0tgEtuwZBQB+AVBLcNiHqAl27Ds2G/gtE+Vp7Nh3rEUtorRj8CdqghNjE75uJ4oagz9RE6zozPm6\nnShqDP5ETTCwvhu5bGbObblsBgPru1vUIko7TvgSNYE1qctsH4oLBn9qe3FJsezvzdd93ri0k9KB\nwZ/aWhxTLJ2CPIDYtZPaG4M/tbV6KZZ+g6qJnrnbl9FpCzqMtZPICwZ/Srx6QdlEimWhWMI9TxzF\nu+Pl2dtqe+Zubai9fXxyyjHI194WpJ1Efoiqhn8Qke8B+FcA3lTVf+HwewHw1wCuAjAO4Cuq+ot6\nj9nX16fczIUaqe1JA5Usmu0be9Dfm8e6of0o1QmgHQLMKCAAwn8SzMt35nBg8IpWN4MSREQOqWpf\no+NMpXr+TwBX1vn9FwBcWP25BcB/N/S8lHKNVs46pVjazVQjfqsD/5LFWddU0EKxhHVD+3He4F6s\nG9qPQrHUolZSOzHS8wcAEVkJ4B9cev5/A+Anqvpo9foxAJ9W1TfcHo89f7KGTEpjE8iIYFoV+Zph\nnfMG97Y8cJskAqhi9u8cefUkHjl4fM7fmM0ITl+4AKcmyswKonm89vybNeafB/Ca7frr1dtcgz+1\nt0aTp7XDOdPVTkrtWPuKzlzdYZ2ksfpipbEJ3LbzsOMx5WnF2ER59jhmBVEQsZrwFZFbUBkWQldX\nV4tbQ1FplH5ZKJbwzV1HZgN+Lfuwztj4ZHMaHWPMCqIgmhX8SwDOsV0/u3rbHKr6AIAHgMqwT3Oa\nRs3WaJx+y55R18BvqdczTiNmBZFfzQr+wwBuFZHHAKwBcKreeD+1t3rpl05fDNSYVSCOq4TJKyPB\nX0QeBfBpAMtE5HUAdwPIAoCqfgfAk6ikeb6ASqrnvzHxvJRMbuP0Kzpz7MEGVBqbwMrBvfNu43wA\nuTGS6qmqN6nqWaqaVdWzVfV/qOp3qoEfWvEfVPUCVe1RVabxpFi9CpcscWwWN4whNyzpTE3X35vH\ndZfmkREBUElvFChu23m4rTJ34oKvKTlh8CdfTCw4KhRL+P6h0uykriowXp4x3dTQrC8nr3LZDDav\n7fJ9v2ZYfc9TXBxGczD4k2dWimZpbAKKD8eU/QQVK42z1ZO6GREIKoupNq/tchyGumnNOfNuz3YI\nspn5wb0zl8X2jT3Y1t+DGUMLJ00amyjjtp2HcfN3f9bqplBMxCrPn+ItbIVM68ujURqnCZvXduH7\nh0qOXzL22j+WvnOXOmbJON0O1N+UJc4Lzw68eBIrB/di89ouPPP8W8wKSjFj5R1MY3mH+HErpSAA\nXh66uuH9GxVZM8UK7gAaloeIQqFYwjd2HnZ9rW52+GLKZTO47tL8vFIOzeL0hUjJFLfyDtQG6qVo\nurHX52kW62zkwOAVLduxy6kmjxX4t/X3uJ5pAJh3P+uLYc+h1yObG+Eq4fRhz588a1Q+2TrGCmqd\ni7N47/0plGfMv8cEwF9tWl23h+3lbCRKQRdceblf950/xAdT5r8ImnFmRNFiz5+Ma7QJee2Xg33z\nExOsaVb787qdVcRhvUCjPXvD3G8ygsAPNF4YxhXE7YPBn3ypF5iiKs1g9fKdnndgfbfj2Yg1Mduu\nopxUdhsCiuN+yBQcUz3JGNOlGToEuH/Tarw8dLVrcOnvzWP7xh7kO3OzqZtpmLh0WyW97oKlRh7f\n6YulUUE+Shb2/MmYbEYwOW1mfN/P2HPQ4ZWks2/6vmRxFndfs8roJjdWSui2/krmlIn9kCk+2PMn\nI+4sjBoL/AK0LFMnCazhF2tDFwB4vyYLyNScx8MHj88uDHN7zDjMr5B/DP4USqFYwup7nsLDB48b\ne0wGk/q8DL8MrO9GtsNMmYkDL57EmvuerluQj5KHwZ8CKxRL+Itdh+f0QE24/KLl3LC8Di/DL/29\nefzBInOjur/750kM7j7iOr/CTeaTh2P+5Is91S+KFSK5bMec1a/MKJnP62K7McOptu9PKwZ3H8Hz\n910153ZmASUTe/7kWaFYwsDjR2YLu5mW7RAsymaYUdJAveEXew+8w6W6aJiqo+9PKz6+Ze6mMcwC\nSiYGf/Js6/BRI6t185053L9pNZYszs7e1pnLYsf1l7j2VplR8iG39FYAc6quOhXQc6tW6seUVjKB\nrKEdZgElE4d9qC7TwzxWD9UtPTPOK3bjxOn1Wze033GRXUYEM6qO1UrDLBS7bedh3PPEUdf3Bf/P\n4o3Bn1w51fIJw0vuflpX7Jrg1tOeUZ1X58j+5RGm+J5bCQ/+n8Ufgz+5MlGuwW+BtUb1g8hdkKqr\nwIdfBL33PmWkHlPtgjOKJwZ/chW2dsyijMzLDPEirSt2wwp71nT3NaswsPsIyiEX69UuOKN44oQv\nuQqzRmiBIFDgp+DC1jnq781jxxcvQT7kWD0zfZKBPX9yFSax54Xtra2ln1Zhz5rs9185uLfB0e6Y\n6RN/7PmTozAbfZspKkCtFmY9gAJc6RtzDP40R6FYwh/d8SQOvHgy8GMwxa893LTmnFD3t1b68gsg\nnhj8aVahWMJtOw+Hqs7JFL/2sa2/B5vXdoV6DI7/xxeDPwGoBP5v7Dwc6L6duWyqNlJJk239PaGH\n8UpjE+z9xxAnfGl2MVfQ/v7WDczpbmcmtoxkobf4Yc+fQi/m2jp81GBrKG6cCsn5xeGf+GHPP+UK\nxVLoXp3pev4UL/ZV12HeKyeqwz9cvR0P7PmnWKFYwl8EHOendOnvzePA4BW4f9PqwGcB2YzMqTrK\nbKDWYvBPsa3DR2FiIb69NDO1N2sVcZA1AJPTyrr/McJhnxSyTr2DDNd0yNyVv9mM4O5rVhlsHcWd\nNUwz8PgRI/s7cDVwazD4p4B9nLVzcRanxsuBevzrLliK6/u6OGZLs//ntxkYNuSiwNYQddjtJw76\n+vp0ZGSk1c1IPGvrxbA9tDPPWIhn7/icoVZRu1g3tD/UJLCgUgrCy14P5I2IHFLVvkbHccy/zZna\nepGBn5yEXc1tvTM5+dt8RoK/iFwpIsdE5AURGXT4/VdE5C0ROVz9+ZqJ56XGTKRhhl3iT+QFJ3+b\nK/SYv4hkAHwbwOcAvA7g5yIyrKq/qjl0p6reGvb5qLkWZQTb+nta3QyKKdPBmpO/zWNiwvcyAC+o\n6ksAICKPAbgWQG3wpwjZJ3U/kstCBBgbL8/LzvGLG7JQPaaDNSd/m8fEsE8ewGu2669Xb6t1nYg8\nJyK7RSRcrViaw6rNYy2eGZso493xMhThAv+6C5aaaiK1KZPBmhVhm6tZE75PAFipqp8A8DSAh5wO\nEpFbRGREREbeeuutJjUt+UxstF5rQYfg+j6O9VN9Jur+AKwI2womgn8JgL0nf3b1tlmq+o6qflC9\n+iCAS50eSFUfUNU+Ve1bvny5gaalQxTjpFMzyuwLashp3+DOnP8V35dftJyBv8lMBP+fA7hQRM4T\nkYUAbgQwbD9ARM6yXd0A4NcGnpeqohonZfYFeWHV/Xl56GoMrO9GkN0fH332tTnXC8US1g3tx3mD\ne7kdZERCB39VnQJwK4B9qAT1Xap6VETuFZEN1cP+XESOisgRAH8O4Cthn5c+ZOrU2wmzL8gra+7p\n3XH/6cXTtsWmtXNYXAMQDSNj/qr6pKr+kapeoKr3VW+7S1WHq5e3qOoqVb1EVS9X1edNPC9V1J56\nBzntdsPsC/Iq7NyTFdydHodnoeaxtk+b6O/NzxkzXX3PU6EXeDH7gvwIe5Zo1QlyexyehZrF8g5t\nanIqXPYPsy/ILxNniVuHj7o+Ds9CzWLwb0OFYgnj5WCV+nPZDO7ftBoHBq9g4CdfnOaectmMr/0e\nxibKro/Ds1CzGPzbUJCxUStNj719Csop7XP7xh7cfc0qZDu8pwCNvHrS8XH4vjSLY/5tKEiJ3ZeH\nro6gJZQ2tXNPlpFXT+Lhg8c9Pcajz76Gbf0M9lFj8G8zhWJptka6V/bTcm6wTVF45nnvK/anY7rH\nSLvhsE+b2bHvmK/AD2B2G0bmV1NU/Gbq8D0XPQb/NhMkHc7q2TO/mqLiN1Nn6/DRiFpCFgb/NtPp\nI7MCAO7ftHr2MvOrKSoD67vhp+qDiU2IqD4G/zZg1UFZObjX99J6+3g+86spKv29edzsc0c4Dv1E\ni8E/4ezj9H7Ze/2Ae54286vJhL5z/e0P8c1dhyNqCQHM9km8MPVUarN47GP/zPYh0/zOHU0z6SdS\nDP4JZE/HNP35cMvTJgoryNzRnYVR7iEdEQ77JExtOmZQm32OvxKFFWTuyOvCMPKPwT9hTG3ZuK2/\nhxtmUFNFue8E+cfgnzCm0i65oIuazV77x487C6MRtSjdGPwTxu3UOeNj77xctoMLuigxHnmWQz9R\nYPBPGLd0TD/1ULZv/AQXdFHTBU1LZqmfaDD4J4xT2dzrLvWXndPfm+eCLmq6MPNVnJ8yj6meCVSb\njtl771O+H2NgfTe27Bmd82Hkgi6KUpizSvt71ZqfAuavVSHvGPwTysr193sKve6CyipLLuiiZlvR\nmQu0Eh2A6/wU36/BMfgnkDV2GuQU+pF/+8nZy1zQRc3kdLYZBuenwuGYfwKZyvUnaian+aowOD8V\nDnv+CcQeDyVV7dnmysG9gR6H81PhseefQOzxUDvwm7HDDd3NYvBPoIH13chm/GyNQRQv1ryVH789\n9T5uXtuFA4NXMPAbwOCfQCOvnkQ5YL1b5kdTHASZt5pWxcMHj7PcgyEM/glTKJbwSIhKhyzfQHEQ\nZt7q71juwQgG/4TZse9YqFLOnCymOAgzbzWjPIM1gcE/YcIGb04WUxyELe98287DLPMQEoN/wnwk\nlw18X6bHUVxYOf9hsAx5OAz+CVIoljA2UQ58f6bHUZz09+ZDL/RiGfLgGPwTZOvw0cD3XbK4csbA\nyogUJwPru0MHIc5jBcMVvgkSptd/9SfOYmVEiqWZkPfnPFYw7PmnxDPPv8Wduyh2/L7/nDYy4jxW\nMOz5J8iSxVm8Ox6s98+duygurHLkJ6r7R/uxfWMPy5AbYiT4i8iVAP4aQAbAg6o6VPP70wD8LYBL\nAbwDYJOqvmLiudPk7mtW4ZuPH8H0jP9Mf7da6jxl9s8evBiA/AlTjhxgGXKTQg/7iEgGwLcBfAHA\nxQBuEpGLaw77KoB3VfXjAP4KwF+Gfd60Cvof5rb3L0+Z/bHvQ6tguqFfYcuR3/zdnxlsTbqZGPO/\nDMALqvqSqk4CeAzAtTXHXAvgoerl3QA+IyKsTObTjn3HUA7Q6weca6kz9dM/p+DFuRPvwg4zHnjx\npKGWkIlhnzyA12zXXwewxu0YVZ0SkVMAPgrgbftBInILgFsAoKury0DT2kvYDw5PmcPj3Ek4YbZy\ntBSKJb6PDYhVto+qPqCqfarat3z58lY3J3Y4Pt96bv8H/L/xxsQw45Y9o7izMMo1KyGZCP4lAOfY\nrp9dvc3xGBFZAOAjqEz8kg9h6qHwQ2IG507C6e/Nzy44DGqiPI1HDh7nvEtIJoL/zwFcKCLnichC\nADcCGK45ZhjAl6uXvwhgv6qGKU6ZSrXj9p25LDo8zpxYH5KB3Uf4IQmBcyfh3X3NqlBF3QDMSxHl\nvIt/YiIGi8hVAO5HJdXze6p6n4jcC2BEVYdFZBGA/wWgF8BJADeq6kv1HrOvr09HRkZCt63dXfyf\nfojxsr81kh0CqIJpitQyYXL93QiAl4euNvRoySUih1S1r+Fxce2AM/h/yPqglMYmINXADYRb9GXJ\nZTPsuVJL+d3EXTC/5w9UzsIODF5hpE1J5jX4x2rCl+az55UDHwZ+AKEDP8DTZWqtIEOQN6/t4ryL\nAQz+MRd2UYwXTFOkVvHb8ch35rCtv4fzLgawtk/MNSMwM02RWsXv+9vq3XPNSnjs+cdc1IGZp8vU\nSn7f3wz45jD4x1zYvU7tctkMNq/t4ukyxcblF3ExZ6tw2CfmrMDslO3jFwM9xUmhWML3D3HNSasw\n+CeA2/jmuqH9oeukELVKM5IZyB2HfRLM71g9UzopTvxO9p55xsKIWpJODP4J1t+bx4V/eLrn40tj\nE6zvQ7Hhd7J3QcbM3BdVMPgn3NvvTfo6nkWwKC4G1nfDz6YeHOI0i8E/wQrFUqBVvlzVS3HQ35v3\nVdcnw/2fjGLwT7AwAZyreikOsj4i0HRM65AlFYN/goUJ4FzVS61UKJbQe+9T8FOQNs/3rFEM/gnW\nGXBTjGyHYHxyihu8UEtYxQr9DllyJbpZzPNPqEKxhPfenwp03xl8WBHU2uAF4NJ5ao6g+f18f5rF\nnn9C7dh3DOWZYGOg0zX3K08r7nniqIlmETUUZLiSk73mMfgnlOkJWxN7AxB5EWS+ae35SyJoSbox\n+CcUJ2wpqYKM3b/yDrPTTGPwTyiT1T6BymbwRM3Q35vHYj85nmBqchQ44ZtQ1uTXN3YeDr0BdrZD\nsHXDqvCNIqojzKbtPNM1jz3/BPO7QtJJZy6LHddfwkwKipR9L2q/71luOBQN9vxT7oMpH6tsiAIo\nFEv45q4jgVfoch+KaLDnn3BLAi70srDOD0XJ6vEHDfyduSwDf0QY/BPu7mtWIZsJlwPNyTSKStgN\nWzgXFR0O+yRc7TaPQXAyjaLCjkV8Mfi3Afs2jysH9/q6bzYjnEyjyKzozIWqw79lzygAlnaIAod9\n2kChWMK6of04b3Cv72Xw5WnlB4siE3Y9CuekosOef8JZE2rWuCprnlOcWB2L/7jnOYz7qd9sw6Gj\naLDnn3BhJ9QA4Obv/sxQa4jm6+/NY8nppwW+P+ekosHgn3AmekUHXjxpoCVE7oKO+3OBV3QY/BPO\nrVfkt1bPx7f4mygm8irMZkFc4BUdBv+Ec5pQy2Uz8Fv+fErDfUgpnezJBm67wnHCNp4Y/BOuvzeP\n7Rt7kO/MQVDZ53T7xh6MBajPf9vOw/wCIM9q6/WUxiawZc/onPdQoVgKnerJ92Q0RGOaHdLX16cj\nIyOtbkZirRvaH+hDt2RxFsW7Ph9Bi6jduL3H8p05HBi8AncWRvHwweOhn8d6PPJGRA6pal+j49jz\nb1OXX7Q80P3eHS9zY3fyxC3Z4MTYBArFkpHAX+95KJxQwV9ElorI0yLym+q/jnutici0iByu/gyH\neU7yZu9zbwS+r9spPJGdW7LBis5coHF+tyKFnSGLF5KzsD3/QQA/VtULAfy4et3JhKqurv5sCPmc\n1EChWDKyJy9XV1I9bskGA+u7fQ85nnnGQriNQL/3/hQ7IREIG/yvBfBQ9fJDAPpDPh4ZsHX4qLHH\n4ik3uXFLNhh51f+6kX96fxpjE84dlvKMshMSgbDlHc5UVWt84bcAznQ5bpGIjACYAjCkqoWQz0su\nCsWS64coCK6upHrsRQUBBJrkzWU7Gq5SZyfEvIbBX0R+BOBjDr+6w35FVVVE3FKHzlXVkoicD2C/\niIyq6osOz3ULgFsAoKurq2HjaT7TPaTxycopNxfakBd/96z/Sd4JDzV/2Akxr2HwV9XPuv1ORH4n\nImep6hsichaAN10eo1T99yUR+QmAXgDzgr+qPgDgAaCS6unpL6A5TPeQ3h0vs6wueVIoljATwaeW\nJR6iEXbMfxjAl6uXvwzg72sPEJElInJa9fIyAOsA/Crk85KLKHpInPileqxVvrftPGzsMa0F6tY8\nAjse5oUd8x8CsEtEvgrgVQA3AICI9AH4uqp+DcAfA/gbEZlB5ctmSFUZ/CMysL57TolnoNJzuu7S\nPH7wixJ+PxmsAijHXMlJbUlxE/KdOQys72bAj1io4K+q7wD4jMPtIwC+Vr38fwD0hHke8s6+reOJ\nsQmssH2QtvX3oFAs4Z4njvpOBV28MPiGHNS+TJQUt+Nq3ubhZi5tqDYDw+l3q+76377OAn4/Oc2J\nXwJQ6e1bnYugQ/yb13bh+4dK885QObbfPAz+KTUeYPjn9t1HsHX46JxU0s5cFls3rOKXQkqYGObJ\nZTvQd+5S9J271PEMlZqDwT+lgmysPTmtmKxZQzA2UcbA40cAMBsoDUwM80yUZ7Blzyi2b+zhEE8L\nsbBbSg2s70a2w2fRfxflGcU3dx1hQbgUCFOe2Y4ZZK3H4J9S/b157Lj+EmOPN63KgnApkPG7S1Ad\nzCBrLQb/FOvvzWPzWvMrqdmra1/TBvf/4Krd1mLwT7lt/T2RfAGwV9ee3Mou+8XMntZj8Cds6+/x\nveF7I6zB3p5OhSgVzlW78cJsn5SzcrZNVgIFPqzBzg94sthz+O3pl4ViCQOPH0bjEmzOuGo3friH\nb4pFsTTfjqs1k8Xp/ZDLZnD2kkX4zZu/9/14GRHctOYcbOvnAv9m8rqHL3v+KWZ6aX4tjvsni9P7\nYaI8HSjwn3nGQjx7x+dMNY0iwOCfYlEH50XZDqwb2s8VnAlh6v0gAAN/AnDCN8WiTrWbKM+gVK3/\nwvz/+DP1fnh56Gojj0PRYvBPMacNuKPE/P94M/F+uH/TakOtoagx+KdY7QbcnbksshlzKzidcB4g\n3hZlg4eEzWu7OKyXIBzzT7na8s/2VD8RGN+Wj6s64ynIxut2AqDv3KXmGkSRY/CnOexfBqZTQbMZ\n4arOGCoUS6ECPwAoKtlC7PknB4d9yJV9WAj4cIVmUOVpxZY9z7H6Z8yY2nuXQ3rJwp4/1TX/TOA5\nTJSDrvPE7H2t7B/rOaj5CsUSvmFw03UO6SULgz95Zn0RrBzca+TxrOwfBv/mKhRLuH33EUxOm5vQ\nqS3U5lYmguKDwZ88CbrxeyNBNwdhcHHn9NoAwB0/GPW1b7NXSxZncfc1q1zniniWF08M/tRQoVjC\nwO4jKBvsKdrdWRj1XP+lUCzN20c4zcHFCvSlsQlkRDCtCgFmN1YvjU0YG9N3s3jhgjmvu1uZCJ7l\nxQsnfKmhHfuORRb4AeDhg8ex5r6nGx5n9SidKpCmcQGZ9XpYZ0/WRivNLtVYO9HrNvHLCeF4YfCn\nhprxof3dP09iZYMsoK3DR+umnaYtuERdmM+r2olet4lfTgjHC4d9qKEVnTnXsXlrqMEUa5ji8ZHj\neOWdidlx68svWt5wz4G0BZdmf9l15rL4YGpmXsnn2rUbA+u7HUtDc41HvDD4U0MD67sdx/yzHTK7\nCbzpfQEOvHhy9nJpbAKPNFiElJbgYp/M7TD8xVuPANi6YRUANJxot65zQj7eGPypIetDa8/26cxl\nsXXDqnkTfUGzdxqpF+Jqs03a1Z2FUTxy8Pjsa9HMwH+zrW6Pl9e5tmwIxQ938iLjzhvc27RJxyWL\nsyje9fkmPVt9UaafWguyon5dOwT40pouPPP8W+y1JxR38qKWqTdHYFKHAHdfsyry5/Ei6tz2HfuO\nNeULdUaBZ55/i9tvpgCzfci4gfXdyHZEWxoaqASqkVdPNj7Qp0KxhHVD+33VIKqX226iDc34MrWk\nLWsqrdjzJ+P6e/ORrAZ28uizr7kuEHNaAJVvMIwRtAdvMre9dlFdMwM/kL6sqbRi8KdIjNUJ/NkO\nwR8sWmDky2Fa1bWcgT2IW5OjtcG89r7jk1MNV6c6Pd9HclnHVNQVnTlfcwE3f/dnczKdoiAAFmUz\njtlZacmaIk74UkTchipEKplCY+NldC7O4tR4GcFrhFZkOwTlml1n7CUOnFhnAF5TVAWVvWlrM26A\nyj4F0zPquPHNuguW4hfHT83Led++sWfeJjpRl2EAKtssumVlZUTwn2+4hJO7CccJX2opp8CazQig\nmO3xvzteRjYjOK1DQpWJrg38QOMSB6WxCV8rZK0evNOmJ/VKXzj14ifK07ht5+GmBHu7jFTmYdyG\nomZUGfhThBO+FIna/YHznTmcvnDBvEBdnlZMTrXm7NPPWHozCqRFbVoVW/aMonNx1vH3HOtPF/b8\nKTK1C33Oc9kHwPRipQ0dP8XtC3ZhhbyNE7oM35q6AcMznzL6HEk1UZ7GaQs6kKsZ8+dYf/qE6vmL\nyPUiclREZkTEdYxJRK4UkWMi8oKIDIZ5Tkout56lNRxhwoaOn2Io+yDO7ngbHQKc3fE2hrIPYkPH\nT409R9KNTZTnnZXVzkFQ+ws77PNLABsB/KPbASKSAfBtAF8AcDGAm0Tk4pDPSwk0sL4buWxmzm25\nbAY3rTln3u1B3b5gFxbL5JzbFsskbl+wy8jjx11GBK8MXY37N62ue0x/bx4HBq/Ay0NX48DgFQz8\nKRQq+Kvqr1W10SqWywC8oKovqeokgMcAXBvmeSmZnOYBtm/swbb+Hmzf2IPOnPNYtB8r5G2X298J\n/dhJYA2h1QvmzaoJRPHWjDH/PIDXbNdfB7DG6UARuQXALQDQ1dUVfcuo6eoV/Ppgam7GTzYjOH3h\nApyaKM/mx4+8ehKPPvuaawA7octwtsMXwAn9aPjGJ4B9CC3vUmYjz4ldgoeev4j8SER+6fBjvPeu\nqg+oap+q9i1fvtz0w1OMOaVdlqcVp5+2YM7QxLb+Hry4/Sq4zRJ8a+oGjOvCObeN60J8a+qGiFoe\nL/YvRbdhNk7sEuCh56+qnw35HCUA59iun129jWiW3/IIbsXjhmc+BZRRzfZ5Byf0o22R7eN10xx7\nr5519an+3DWtAAAGDUlEQVSeZgz7/BzAhSJyHipB/0YAX2rC81KCuAVztwwht9W5HVL5AhieTHaw\nt/O6GtmpV8+6+uQmbKrnn4nI6wA+CWCviOyr3r5CRJ4EAFWdAnArgH0Afg1gl6oeDddsajd+hyic\nJo/v37QaX1rT5Tok5EQAvDJ09Zyf+zetDjQuvu6CpfP+Bqn5N4gTYxOOf+/mtV1M16TAWNuHYiPs\nZii1FTm9yHfmGtauLxRL2Dp81HUPYXutntpjT1+YwfjktKda/B0Cx/pAXtpIZGFtH0qcsEMUfmr1\nAP4mP2szkazCcU4lou3H/n7SW3ty2QyuuzSP7x8qceUtNQWDP7UNL7XzMyKYUfV1ZuH0pWIF/toe\nuZ8vIKe29J27lBO01BQM/tQ2Gm0f6VRK2Qs/mUheN29xawsnaKlZWNWT2obTpLE10RpmQtQt48jp\ndi+VMTMinJyllmPPn9pGVHntTmmWbmPxjVIyg559EJnG4E9tJYphEz9fKrXHdi7OQhVzSlQw8FMc\nMNWTiKiNeE315Jg/EVEKMfgTEaUQgz8RUQox+BMRpRCDPxFRCsU220dE3gLwaqvbYbMMgPMegfGU\ntPYCyWtz0toLJK/NSWsv0Po2n6uqDXfDim3wjxsRGfGSPhUXSWsvkLw2J629QPLanLT2AslpM4d9\niIhSiMGfiCiFGPy9e6DVDfApae0FktfmpLUXSF6bk9ZeICFt5pg/EVEKsedPRJRCDP4uROR6ETkq\nIjMi4jpzLyKviMioiBwWkZZVovPR3itF5JiIvCAig81so0NblorI0yLym+q/S1yOm66+vodFZLgF\n7az7monIaSKys/r7Z0VkZbPb6NCmRm3+ioi8ZXtdv9aKdtra8z0ReVNEfunyexGR/1L9e54TkT9p\ndhtr2tOovZ8WkVO21/euZrexIVXlj8MPgD8G0A3gJwD66hz3CoBlSWgvgAyAFwGcD2AhgCMALm5h\nm78FYLB6eRDAX7oc914L29jwNQPw7wF8p3r5RgA7W/xe8NLmrwD4r61sZ017/iWAPwHwS5ffXwXg\nh6jsz7MWwLMxb++nAfxDq1/Xej/s+btQ1V+r6rFWt8Mrj+29DMALqvqSqk4CeAzAtdG3ztW1AB6q\nXn4IQH8L2+LGy2tm/zt2A/iMiAhaJ27/zw2p6j8COFnnkGsB/K1WHATQKSJnNad183lob+wx+Ien\nAJ4SkUMickurG9NAHsBrtuuvV29rlTNV9Y3q5d8CONPluEUiMiIiB0Wk2V8QXl6z2WNUdQrAKQAf\nbUrrnHn9f76uOoSyW0TOaU7TAovbe9eLT4rIERH5oYisanVjaqV6Jy8R+RGAjzn86g5V/XuPD/Mp\nVS2JyB8CeFpEnq/2Cowz1N6mqtdm+xVVVRFxSz07t/oanw9gv4iMquqLptuaMk8AeFRVPxCRf4fK\nmcsVLW5TO/kFKu/b90TkKgAFABe2uE1zpDr4q+pnDTxGqfrvmyLyA1ROuSMJ/gbaWwJg7+GdXb0t\nMvXaLCK/E5GzVPWN6in8my6PYb3GL4nITwD0ojKm3QxeXjPrmNdFZAGAjwB4pznNc9Swzapqb9+D\nqMy/xFnT37thqOo/2S4/KSL/TUSWqWps6hRx2CcEETldRM6wLgP4PADH2f+Y+DmAC0XkPBFZiMrk\nZNOzZ2yGAXy5evnLAOadvYjIEhE5rXp5GYB1AH7VtBZ6e83sf8cXAezX6qxfizRsc814+QYAv25i\n+4IYBvCvq1k/awGcsg0Zxo6IfMya9xGRy1CJta3sEMzX6hnnuP4A+DNUxhU/APA7APuqt68A8GT1\n8vmoZFIcAXAUleGX2La3ev0qAP8PlZ5zy9pbbctHAfwYwG8A/AjA0urtfQAerF7+UwCj1dd4FMBX\nW9DOea8ZgHsBbKheXgTgcQAvAPi/AM5v5evqsc3bq+/ZIwCeAXBRi9v7KIA3AJSr7+OvAvg6gK9X\nfy8Avl39e0ZRJwMvJu291fb6HgTwp61+T9T+cIUvEVEKcdiHiCiFGPyJiFKIwZ+IKIUY/ImIUojB\nn4gohRj8iYhSiMGfiCiFGPyJiFLo/wMidWlG/YFlSAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "coords = np.zeros((2,sim.N))\n", "for i in range(sim.N):\n", " coords[0][i], coords[1][i] = sim.particles[i].x, sim.particles[i].y\n", "fig, ax = plt.subplots()\n", "ax.axis('equal')\n", "ax.scatter(coords[0],coords[1])\n", "ax.scatter(sim.particles[1].x,sim.particles[1].y); # Planet" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.6.5" } }, "nbformat": 4, "nbformat_minor": 1 }