{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hybrid integrations with MERCURIUS\n",
    "REBOUND comes with several integrators, each of which has its own advantages and disadvantages. MERCURIUS is a hybrid integrator that is very similar to the hybrid integrator in John Chambers' Mercury code (J. E. Chambers 1999). It uses a symplectic Wisdom-Holman integrator when particles are far apart from each other and switches over to a high order integrator during close encounters. Specifically, MERCURIUS uses the efficient WHFast and IAS15 integrators internally."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start out by showcasing the problem with traditional fixed timestep integrators such as WHFast. We setup a simulation of the outer solar system and increase the masses of the planets by a factor of 50. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<rebound.plotting.OrbitPlot at 0x107a69be0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbIAAAHDCAYAAABBDZ94AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB5f0lEQVR4nO3dd3yT5fo/8E920pXupnvQAbLLatkbAfFwQMAtyEFUXMBRwQHiguNCRQ84AQWBAwoqIFDKKtIyWsruoHTvkSbpyH5+f/hrvlQKdCR5Mq7365UXNkmfXIaST+/7uZ775jAMw4AQQgixU1y2CyCEEEK6goKMEEKIXaMgI4QQYtcoyAghhNg1CjJCCCF2jYKMEEKIXaMgI4QQYtcoyAghhNg1CjJCCCF2jYKMEEKIXbOrICstLcWjjz4KHx8fSCQS9O7dG+fOnTM9zjAMVqxYgcDAQEgkEowfPx65ubksVkwIIcTS+GwX0F5yuRzDhg3DmDFj8Mcff8DPzw+5ubnw8vIyPeeDDz7A559/js2bNyMyMhJvvvkmJk2ahKtXr0IsFt/1NYxGI8rKyuDu7g4Oh2PJ/x1CCCF3wDAMVCoVgoKCwOXeZczF2IlXX32VGT58+G0fNxqNjEwmYz788EPTffX19YxIJGK2bdvWrtcoLi5mANCNbnSjG91s5FZcXHzXz267GZH99ttvmDRpEmbNmoXjx48jODgYzz77LBYsWAAAyM/PR0VFBcaPH2/6HqlUiiFDhiA1NRUPPvjgLcfUaDTQaDSmr5n/vxFAcXExPDw8LPx/RAgh5HaUSiVCQ0Ph7u5+1+faTZDduHED69evx5IlS/Daa6/h7NmzeOGFFyAUCvHEE0+goqICABAQENDq+wICAkyP/d3q1auxatWqW+738PCgICOEEBvQntM8dtPsYTQaER8fj/fffx/9+/fHU089hQULFmDDhg2dPuby5cuhUChMt+LiYjNWTAghxBrsJsgCAwNxzz33tLqvR48eKCoqAgDIZDIAQGVlZavnVFZWmh77O5FIZBp90SiMEELsk90E2bBhw5Cdnd3qvpycHISHhwMAIiMjIZPJkJycbHpcqVTi9OnTSExMtGqthBBCrMduzpEtXrwYQ4cOxfvvv4/Zs2fjzJkz+Prrr/H1118D+Gse9aWXXsK7776LmJgYU/t9UFAQpk+fzm7xhBBCLMZugmzQoEHYvXs3li9fjrfffhuRkZH49NNP8cgjj5ie88orr6CxsRFPPfUU6uvrMXz4cBw4cKBd15ARQgixTxympeecQKlUQiqVQqFQ0PkyQghhUUc+j+3mHBkhhBDSFgoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF3js10AIcT6GIaBVqsFAHA4HAAAl8sFn08fCcT+0E8tIXZOpVKhsrISCoWi3TeVSgU/Pz9UVVWZjtO3b1/k5OTAzc2tzZurqyuCgoIgk8kQEhKC0NBQSKVSUxASwhYKMkLshFarRV5eHrKyslrdAMDb2xs3btwwPVcikcDT0xMeHh6QSqWQSqUICAgw/beHhwdcXFxMIcQwDDgcDpRKJRoaGtDY2AiVSoXGxkY0NDSgqKgIDQ0NqK2tRX19vel13NzcEBoaagq2lj8jIyPh7+9PIUesgoKMEBvDMAyqqqpuCay8vDwYDAYAQGhoKOLi4vDggw8iNjYWgYGB8PT0NIWUQCCwWH0NDQ0oLi5GcXExSkpKTH8ePHgQ5eXlMBqNcHFxgVgsRkJCAiIjIzF48GD06tULQqHQYnUR58VhGIZhuwhboVQqIZVKoVAo4OHhwXY5xEkwDIPc3FykpaXh1KlTUKlUSEtLAwC4uLige/fu6NGjB7p3747u3bsjNjYWbm5uLFfdNr1ej7KyMuTn5yMrKwvl5eU4evQolEolRCIR+vfvj8GDB2Pw4MHo0aMHuFzqNyNt68jnMQXZTSjIiLWUl5cjNTUVR48eRWpqKmprayEQCNC/f3+MHj0aMTExiIuLQ3BwsN1/2BuNRuTk5ODs2bM4ffo0MjIy0NzcDDc3NwwYMACDBw/GoEGDEB0dTVORxISCrJMoyIiltIy6Dh48iEOHDuHixYuIj48Hh8PB0KFDkZiYiPj4eEgkErZLtTi9Xo8rV67gzJkzOHv2LDIzM6HT6eDl5YXBgwdj9OjRGDp0KP0bdHIUZJ1EQUbMiWEYXLlyBb/++isOHTqEgoICuLq6YvTo0Zg4cSJGjx4NT09PtstknUajwYULF3DmzBmcOXMG169fh16vx+DBgzFt2jSMGDECLi4ubJdJrIyCrJMoyIg51NXV4ZdffsG2bdtQVVUFkUiE8ePHY9KkSRg6dChEIhHbJdq02tpaHDlyBElJSaipqUFtbS0mTpyIf/7zn+jRowdNPzoJCrJOoiAjnWU0GvHnn39i27ZtOHDgABiGweTJk/Hggw9i2LBh4PF4bJdolyoqKvDbb7/ht99+Q2VlJWJiYjB9+nRMnjwZ7u7ubJdHLIiCrJMoyEhHlZWV4X//+x927NiBkpISxMTE4KGHHsLMmTPh7e3NdnkOw2g0Ii0tDbt370ZKSgp4PB7Gjx+P6dOno1+/fjRKc0AUZJ1EQUbaw2Aw4NixY9i8eTOOHTsGiUSC+++/Hw8++KCpgYNYTm1tLX7//Xf8+uuvKCkpQUJCAu6//36MGzfO7js8yf+hIOskCjJyJxqNBnv27MFnn30Gg8EAf39/PPzww5g2bZrNXtflyIxGI86ePYtt27YhNTUV4eHhmDt3LiZNmkRTuQ6AgqyTKMhIW3Q6Hf73v//h448/hru7O3r27ImFCxeib9++bJdG/r+rV6/i+++/R0pKCoKDg/HEE09gypQpFl3hhFgWBVknUZCRmxmNRvz666/48MMPkZ+fj+nTp+Pll19GVFQU26WR28jJycH333+Po0ePQiaT4bHHHsP9999PS2PZIQqyTqIgI8Bf138lJSXhP//5D65evYoJEybg1VdfRc+ePdkujbTTjRs3sHHjRiQlJSEwMBBPPPEEpk2bRlOOdoSCrJMoyMiff/6J1atXIz09HUOHDsXy5csxcOBAtssinVRYWIhdu3Zh165diIiIwPPPP4+hQ4eyXRZpBwqyTqIgc14ZGRlYs2YNUlJS0K9fPyxfvhwjRoygDkQHkZ2djc8//xwqlQohISFYunQpfHx82C6L3EFHPo+pV5U4tZycHMybNw9Tp05FVVUVNm7ciP3792PkyJEUYg4kLi4OX3zxBebOnYuLFy/i4Ycfxt69e0G/xzsGGpHdhEZkzkOlUuGjjz7Cnj174ObmhiVLlmD69Ol0DsUJKBQKrFu3Dn/88QcGDRqEV155BUFBQWyXRf6GphY7iYLMORw7dgyLFy9GeHg4JkyYgAULFlBXmxM6ffo0/vOf/0ChUOCpp57CrFmz6IJqG0JTi4S0QaVSYcmSJZg1axYiIyOxbt06LFq0iELMSQ0ZMgRbt27FtGnTsG7dOixcuBD5+flsl0U6gYKMOIWjR49ixIgR+OWXX/Dhhx/i559/RmhoKNtlEZZJJBK89NJLWL9+PRobGzF37lxs3LgROp2O7dJIB1CQEYfWMgqbPXs2oqKikJKSgrlz51IjB2mld+/e2LRpEx599FFs3LgRTz75JK5cucJ2WaSd6BzZTegcmWM5evQoFi9ejPr6eqxatQqPP/44BRi5q+vXr+P9998HAIwZMwaPPPIInTtjAZ0jI05NpVJh8eLFmD17Nrp164aUlBQ88cQTFGKkXaKjo/H1118jISEB33zzDV5++WXI5XK2yyJ3QEFGHMrZs2fx9NNPY/fu3fjoo4+wa9cuOhdGOozP5+Opp57Cxx9/jOrqajz//PO4fv0622WR26AgIw6BYRhs2LAB9913H+RyOY4dO0ajMNJlgwYNwkcffQSJRIJFixYhJSWF7ZJIGyjIiN1TKpWYO3cuXn/9dSxcuBC//vorIiIi2C6LOAh/f3+sW7cOCQkJeOONN7BlyxZaEcTG8NkugJCuyMrKwvLly5GZmYkff/wRU6ZMYbsk4oDEYjFWrlyJ8PBwfPPNNygoKMArr7xC1yDaCBqREbt1+PBhTJgwATqdDkePHqUQIxbF5XLx5JNPYuXKlTh+/DhefPFF1NXVsV0WAQUZsUMt58PmzJmD4cOHY/v27TSVSKxm7NixWLduHSorK7Fw4UJqArEBdhtka9asAYfDwUsvvWS6T61WY9GiRfDx8YGbmxtmzpyJyspK9ookZqfT6bB06VIsX74czz33HLZs2QI3Nze2yyJOpnv37vj666/h5eWFjz76CBkZGWyX5NTsMsjOnj2Lr776Cn369Gl1/+LFi/H7779j586dOH78OMrKyjBjxgyWqiTmVl9fj1mzZmHLli1Yt24dVq1aRavVE9b4+vris88+g4eHB5YvX44zZ86wXZLTsrsga2howCOPPIJvvvkGXl5epvsVCgW+++47fPLJJxg7diwGDBiAjRs34tSpU0hLS2OxYmIOBQUFmDBhAi5evIjdu3fj0UcfZbskQiCRSPDee+9h4MCBeP3113Hq1Cm2S3JKdhdkixYtwtSpUzF+/PhW96enp0On07W6v3v37ggLC0Nqamqbx9JoNFAqla1uxPZcuXIFEyZMQExMDJKTkzFs2DC2SyLERCAQYNWqVRg6dChWrFiBEydOsF2S07GrINu+fTsyMjKwevXqWx6rqKiAUCiEp6dnq/sDAgJQUVHR5vFWr14NqVRqutEKELbn3LlzmDJlCmQyGT7//HNERkayXRIht+Dz+VixYgVGjRqFVatW4ciRI2yX5FTsJsiKi4vx4osvYuvWrRCLxWY55vLly6FQKEy34uJisxyXmMfJkydx//33o3v37ti7dy98fX3ZLomQ2+LxeHj99dcxfvx4vPvuuzh06BDbJTkNu7kgOj09HVVVVYiPjzfdZzAYcOLECXzxxRc4ePAgtFot6uvrW43KKisrIZPJ2jymSCSCSCSydOmkEw4ePIjHH38ciYmJ+Omnn+Di4sJ2SYTcFZfLxbJlyyAQCLBmzRrodDpMnTqV7bIcnt0E2bhx43Dp0qVW982bNw/du3fHq6++itDQUAgEAiQnJ2PmzJkAgOzsbBQVFSExMZGNkkkn/fLLL1iwYAHuvfdefP/99/TLBrErHA4HS5cuhUAgwLZt2yASiW45p0/My26CzN3dHb169Wp1n6urK3x8fEz3z58/H0uWLIG3tzc8PDzw/PPPIzExEQkJCWyUTDph586dWLZsGWbOnIn//ve/4PPt5keUEBMOh4MXXngBH330Ef7zn//Azc2NPocsyG7OkbXH2rVrcd9992HmzJkYOXIkZDIZfvnlF7bLIu30888/41//+hceeOABbNiwgUKM2DUOh4MlS5YgMTERH330Ea5du8Z2SQ6Ldoi+Ce0QzZ59+/bh0UcfNYUYXehMHIVWq8WKFSuQn5+PL7/8kpqW2ol2iCZ25fDhw3j88ccxdepUrF+/nkKMOBShUIhXX30VPB4Pb731FrRaLdslORwKMsKqM2fOYPXq1Rg7diy+//57mk4kDsnLywurVq1Cfn4+Pv30U9rPzMzoU4OwJj8/H7NmzUJcXBx++OEHp97biWEY1NbWoqysDFVVVZDL5aYtQiorK6HRaKDVaqHRaG773xKJBFqtFtHR0airq4O7uzvc3Nzg7u5uurm5ucHDw6PV156envD29qaRsIXFxMRg6dKlWL16NaKjo2kdWDOiICOsqK+vx4wZM+Dp6Ynt27dDIpGwXZLFGY1GVFdXo6SkBKWlpa1uZWVlUKvVAAAXFxdwOBx4eXmhW7duaG5uhlAohLu7O3x9fSEUCiEWiyEUCk3XQgqFQvB4PBiNRgB/vb8qlQoNDQ1QqVQoLCw0fa1UKluNCAIDAyESiSAQCBAREYGIiAiEh4cjIiICvr6+4HA4rLxfjmjs2LHIy8vDhg0bEBER0eq6WNJ51OxxE2r2sA6tVov7778fly9fxrFjxxAdHc12SWan1+tx48YN5OTkIDs7Gzk5OdBqtbhx4waAvy6clclkCAkJQVBQEIKDgxESEoLg4GD4+/tbNNgZhkFTU5Mp2ORyOSorK5Gfn4/CwkIUFBSgubkZwF+h2hJqNwecu7u7xepzdEajEa+//jqysrLw5ZdfIigoiO2SbFJHPo8pyG5CQWZ5DMNg4cKF2LlzJ/bt24ehQ4eyXZJZ6HQ6ZGVlITMzExkZGbh48SKCg4NRWFiI8PBwxMbGIi4uDuHh4QgODkZAQIDNng9kGAbV1dWmUCsoKEBhYSGKi4uh1+sBADKZDP3790dUVBTi4+Nvu3oOaVtDQwMWLVoEoVCIzz77jFauaQMFWSdRkFnemjVr8M4772Djxo2YPXs22+V0mk6nw9WrV5GRkYHMzExcunQJGo0GLi4u6NOnD/r374++ffsiJibGbGuDsk2v16O0tNQUagUFBUhPT4fBYEBQUBDi4+MRHx+PXr160Wos7VBYWIjnn38e8fHxWLlyJU3h/g0FWSdRkFnWjh078OSTT2LFihV49dVX2S6nwxQKBU6dOoWTJ08iLS0N7u7uaG5uRt++fdGvXz/Ex8cjJibGqZommpqacOnSJWRkZCAjIwNVVVUQCATo2bOnKdhCQkLoQ/o2UlNTsWLFCjzxxBO0x97fUJB1EgWZ5fz555+47777MGvWLHz11Vd288FWV1eHpKQkpKSk4Pz58zAajejRowdGjBiBoUOHIiYmBlwuXcUC/DUlWVZWZgq1y5cvQ6vVwtfXF/Hx8UhISECfPn0gEAjYLtWmbN26FZs2bcJ7772HwYMHs12OzaAg6yQKMssoLCzEiy++CLVajd9++83m2+ybm5tx7NgxHDx4EKdPnwaXy8XEiRPRr18/DBs2jFZmaCetVosrV66Ygq2pqQkajQYTJ05EQkIC4uLi7OYXGktiGAZr167FyZMn8dVXX8HPz4/tkmwCBVknUZCZX3NzM8aOHYv6+nqkpKTYbAgYDAacPn0aBw8exLFjx6BWq9G3b19MnjwZY8eOhVQqZbtEu1daWoqUlBSkpaWhuLgYoaGhmDBhAkaNGgU3Nze2y2OVSqXC008/jaCgIPznP/+hUT4oyDqNgsz8XnzxRfzwww84duwY+vbty3Y5t6ipqcHu3btx5MgR3LhxAxEREZg8eTImTZqEwMBAtstzSAzD4OLFi0hKSsKZM2fA5XIxbNgwTJw4EbGxsU47Srtw4QL++9//YsqUKfjHP/7Bdjms68jnsW32/xKH8PPPP+OXX37BF198YVMhxjAMzp8/j507d+Lo0aMQCoWYPHky3nzzTfTo0cNpP0ithcPhoG/fvujbty/q6+tx9OhRJCUl4dixYwgLC8OECRMwcuRIpxul9e3bF71798amTZuQmJgIf39/tkuyGzQiuwmNyMynrKwM8fHxGDt2LLZu3WoT4aDX63HgwAEcOXIEKSkpCA8PxwMPPID77rvP6T40bU3LKO3QoUM4e/YsuFwupkyZgrFjxyIkJITt8qymsbERCxcuRFRUFFatWmUT/27YQlOLnURBZh5GoxH33Xef6Torb29vVuvRaDT47bffsHnzZpSXl+Pee+/F9OnTMXDgQKf+oLBV9fX1SE5ORnJyMqqrqzFo0CBMnz4dsbGxbJdmFampqXj77bexbNkyjBo1iu1yWENB1kkUZObxxRdfYOnSpdi/fz/GjRvHWh2NjY3YtWsXtmzZgvr6ekycOBFz585FTEwMazWR9tPr9UhJScGePXug1WoRGhqKhx56CJGRkWyXZnHvvfceLl26hK+//tppP4soyDqJgqzrrl69ioSEBPzrX//CJ598wkoNzc3N2LZtG3bs2IH6+npMmzYNTzzxBEJDQ1mph3QNwzA4ffo0duzYgbKyMgwdOhSzZ8926GYcuVyOBQsWYOjQoViyZAnb5bCCgqyTKMi6RqvVYvjw4dBqtUhNTbX6ivZ6vR579uzBhg0boFQqMXfuXMycORMBAQFWrYNYhsFgwPHjx7Fr1y7I5XKMHTsWDzzwALy8vNguzSIOHjyITz/9FO+//z769+/PdjlWR0HWSRRkXfP666/js88+w8mTJ9GvXz+rvS7DMEhOTsYXX3yBoqIiTJkyBc8++yytKu6gtFotDh06hN27d0Oj0Zja1V1dXdkuzawYhsGyZctQVVWFDRs2ON36lRRknURB1nkpKSmYMGEC3n33Xfz73/+22uteuXIFP/74Iw4ePIihQ4fixRdfdJqmAGfX1NSEvXv3Yu/evfD398f48eMxYcIEh1rrsrS0FM888wxmzJiBuXPnsl2OVVGQdRIFWcckJSVhxYoVKCkpQXV1Nbp374709HSrfJAolUqsW7cOO3fuRM+ePfHCCy9gyJAhFn9dYnsUCgX27NmDgwcPIiwsDPPnz3eohp5du3bhf//7Hz777DOHPi/4dxRknURB1n5PPvkkNm7ceMv98+fPx7fffmux1zUajfjtt9+wdu1a6HQ6LFq0CA8++KBD/RZOOicvLw/fffcd8vPzMWbMGDz00EMOsQGoVqvFv/71L/Ts2dMud43oLAqyTqIga5+kpCRMnDjxto8fPnzYIm33eXl5ePvtt5GZmYkpU6Zg6dKlNrt2I2GH0WhEcnIytm/fDplMhqlTpyIxMdHurxc8dOgQPvvsM3z66acONdq8EwqyTqIga5/ExESkpaXd9vGEhASkpqaa7fUMBgO+++47bNq0CVFRUXjhhRdouwtyR0qlElu3bsXJkycxePBgzJs3z67/TRsMBjz33HPw8vLCe++9Z/fB3B4d+TymJZZJh5WVlXXp8Y64ceMGHn30UXz55ZeYM2cOvv/+ewoxclceHh545pln8MILL+Dq1at49dVXce7cObbL6jQej4e5c+fiwoULOH/+PNvl2BwKMtJhd2trN0fbu9FoxKZNmzBr1iw0NDTgxx9/xIsvvmjze5kR2zJkyBB88MEHiImJwdq1a7F+/Xo0NjayXVanDB48GD179sT3338PmkhrjYKMdNjbb799x8fffffdLh2/vLwcy5YtwyeffIIHH3wQu3btQp8+fbp0TOK8pFIpFi9ejGeeeQYZGRn4+OOPkZOTw3ZZHcbhcDBv3jzk5+fj2LFjbJdjUyjISIdNmDAB8+fPb/Ox+fPnd6nR4/jx45gxYwZycnKwceNGvPzyy053ISgxPw6Hg+HDh2PNmjXgcDh47733sH//frsb2fTo0QOJiYn47bffoNfr2S7HZlCQkU759ttvcfjwYSQkJCAsLAwJCQk4fPhwp1vv9Xo9PvroIzz77LMYMGAANm/ejAEDBpi5auLsfHx88Nprr2Hy5MnYtm0bPv/8czQ1NbFdVoc89thjuHHjBo4ePcp2KTaDuhZvQl2L7CgvL8fSpUtx+fJl/Pvf/8Zjjz3mFF1ZhF3p6en4+uuv4ebmhhdeeAHh4eFsl9Ru77//PgoKCrBhwwZwuY45HqGuRWI3UlNT8eKLL6K6uhpbtmzB448/TiFGrGLAgAF45513IJFIsGrVKhw/fpztktpt1qxZKC8vx8mTJ9kuxSZQkBHW/PTTT1iwYAECAwOpoYOwwt/fHytWrMCwYcOwd+9ebN26FUajke2y7iomJgb9+/e3y/N8lkBBRqxOr9fj7bffxjvvvINHHnkEn3zyCaRSKdtlESclFAoxf/58TJo0CUlJSfj000+hVqvZLuuuZs6ciaysLFy7do3tUlhHQUasSqFQYMGCBfjf//6Hd955B8uXL6d1EolNGD9+PJYuXYrs7Gy8++67qKurY7ukO+rTpw9kMhkOHjzIdimsoyAjVlNYWIhFixYhKysLGzduxAMPPMB2SYS00rt3b7z55ptobm7G+vXrUVpaynZJt8XhcDBx4kT8+eefaGhoYLscVlGQkbtiGKbL8/BXr17FQw89BJ1Oh//9738YNGiQmaojxLxCQkKwYsUKqNVqrFmzBgUFBWyXdFvjxo2D0Wh0+lZ8CjJyV10NsTNnzuCxxx5DYGAg1q9fj9DQUDNVRohlSKVSLFu2DP7+/vjggw+Qm5vLdkltkkqlSEhIwIEDB5y66YOCjNxRyz+OzrbEJycnY/78+ejduzd++OEHeHt7m7M8QizG1dUVL7/8MsLCwvDRRx/h6tWrbJfUpkmTJqGkpMSpmz4oyMgdMQzT6RDbs2cPnnvuOYwZMwZff/01XF1dzVwdIZYlFouxZMkSxMbGYu3atbhw4QLbJd2ipenj0KFDbJfCGgoycltdGY3t2LEDn3/+OR588EGsXbuWVq0ndksoFOLFF19Enz598Pnnn+Ps2bNsl9RKS9PHyZMnnbbpg4KM3FZnR2O7du3CG2+8gVGjRmHFihXUXk/sHp/Px7PPPovBgwdj/fr1NreihrM3fVCQkTZ19sTxL7/8gtdeew0PP/wwVqxYQctNEYfB4/GwYMECDB8+HL///vsdd0m3Nk9PTyQkJODgwYNO2fRBQUbuqCNBtHv3bixbtgxz5szBypUrKcSIw+FyuZg3bx7i4uLw3Xff2dRuzZMmTUJxcbFTNn1QkJE2dXRa8cCBA/jkk08wa9YsrFq1ymFX5CaEw+Fg7ty5iI+Px+7du3H9+nW2SwLg3E0f9GlDbtHRJo+0tDS89NJLGDZsGN555x0KMeLwuFwuFixYAHd3d3z55Zeorq5muyRwOBzce++9KCwsRHNzM9vlWBV94pBbdGQ0lpWVhYULF2LIkCEUYsSp8Pl8PPPMM3BxccG6detsIjyGDRsGrVaLzMxMtkuxKvrUIa10ZDRWWlqKefPmITw8HP/9738hEAgsXR4hNsXNzQ3PP/88FAoFvvrqK9a3gPH394dYLMaZM2dYrcPaKMhIK+0djdXX12PevHkQiUT47rvv6GJn4rRkMhmefvppXLt2DTt27GC7HPTv3x+ZmZlO1b1IQUZMjEYjDAbDXYNMq9Xi7bffRn19PTZt2gQ/Pz8rVUiIberRowcefvhhHDlyBMeOHWO1lvj4eCiVSuTl5bFahzXx2S6A2I72TIswDIMVK1Zg79692LZtGyIiIixfGCF2YNSoUSgvL8e2bdvg5+eHnj17slJHXFwcXFxccP78eURHR7NSg7XRiIyYGI3GuzZrbNmyBTt27MB7772HAQMGWKkyQuzD7Nmz0bNnT3z11VcoLy9npQYej4e+ffsiIyODlddnAwUZAfB/o7E7BVlaWhpWrVqFuXPnYtasWdYqjRC70dKW7+3tjXXr1rG29mF8fDxyc3OhVCpZeX1royAjAP4KMg6Hc9vzY6WlpXj22WcxZMgQvP7661aujhD7IZFI8Pzzz8NgMGDr1q2sdDL269cPDMPY5Gr9lkBBRgDceVpRr9fjhRdeQJ8+ffDFF1+Az6dTq4TciY+PD5588kmUlJSwstKGt7c3IiIinGZ6kYKM3HVa8cMPP0RmZiZeeukleHl5WbM0QuxWXFwcBg4ciL1796KkpMTqr+9MbfgUZOSO04rHjx/Hhg0b8Morr6Bfv37WL44QOzZ16lTIZDJs3rwZer3eqq89YMAAp2nDpyAjt51WrKqqwuLFizFq1CgsWLCAhcoIsW98Ph9z585FWVkZ/vjjD6u+dmxsrKkN39FRkDm5200rGo1GLF68GFwuF2vXrqU1FAnppJCQEEyZMgUHDhxAYWGh1V6Xx+MhISEBRUVFVntNttCnk5O73bTili1bcPLkSXz66afw8fFhqTpCHMO9996LkJAQbNq0CTqdzmqvGx4ejvT0dBgMBqu9JhsoyJxcW9OK2dnZWLlyJRYvXozhw4ezVBkhjoPH42Hu3LmoqanB77//brXXjY6Ohk6ns+pIkA0UZE6srWlFvV6PJUuWICIiAs888wxbpRHicAIDAzFt2jQcPnwYN27csMprRkVFgcfj2czmn5ZCQebEdDqdaWqxxVdffYXLly/j448/hkgkYrE6QhzP+PHjERERgc2bN0Or1Vr89YRCIcLDw5Gbm2vx12ITBZkTMxqN4PF4pq9zc3Px8ccf46mnnkJ8fDyLlRHimLhcLp544gnI5XLs2bPHKq8ZHR1NQUYck9FoBMMwpmlFhmHw/vvvY9iwYVi6dCnL1RHiuAICAjB9+nQcPXrUKtd4xcTEQKvVoqmpyeKvxRa7CbLVq1dj0KBBcHd3h7+/P6ZPn47s7OxWz1Gr1Vi0aBF8fHzg5uaGmTNnorKykqWKbdvfz4/t3r0bhw4dwsKFCyEWi9ksjRCHN2bMGPTu3Rt79uyxeEdhZGQkamtrWVldxFrsJsiOHz+ORYsWIS0tDUlJSdDpdJg4cSIaGxtNz1m8eDF+//137Ny5E8ePH0dZWRlmzJjBYtW2q6VbkcPhQKVS4Z133sF9991HXYqEWAGHw8G0adNQWFiItLQ0i75WUFAQOBwOiouLLfo6bLKb1V8PHDjQ6utNmzbB398f6enpGDlyJBQKBb777jv89NNPGDt2LABg48aN6NGjB9LS0pCQkMBG2TbLYDCYzo99+umnUKlUWLFiBctVEeI8QkNDMWDAAPzxxx8YNGgQhEKhRV5HIBBAJpPRiMwWKRQKAH+t8gwA6enp0Ol0GD9+vOk53bt3R1hYGFJTU9s8hkajgVKpbHVzBgzDgGEY8Hg85Obm4ttvv8ULL7yA4OBgtksjxKlMnToVKpUKJ06csOjrBAcHU5DZGqPRiJdeegnDhg1Dr169AAAVFRUQCoXw9PRs9dyAgABUVFS0eZzVq1dDKpWabqGhoZYu3Sa0zMlzOBy88cYbCAkJwdNPP81yVYQ4H19fXwwfPhyHDh2yaDNGaGgoSktLLXZ8ttllkC1atAiXL1/G9u3bu3Sc5cuXQ6FQmG6OPId8s5bzY0eOHMHp06fxzjvvWGxagxByZ5MmTYLBYMDhw4ct9hohISGora112M5Fuwuy5557Dnv37sXRo0cREhJiul8mk0Gr1aK+vr7V8ysrKyGTydo8lkgkgoeHR6ubMzAajTAYDFi5ciXGjRtnOqdICLE+Dw8PjBkzBseOHTOdMjG34OBguLu7o7y83CLHZ5vdBBnDMHjuueewe/duHDlyBJGRka0eHzBgAAQCAZKTk033ZWdno6ioCImJidYu12YxDAOj0Yiff/4ZN27cwJIlS9guiRCnN27cOAiFQott9RIQEIDGxkZUV1db5Phss5uuxUWLFuGnn37Cr7/+Cnd3d9N5L6lUColEAqlUivnz52PJkiXw9vaGh4cHnn/+eSQmJlLH4k2MRiOam5vxySefYMaMGejZsyfbJRHi9CQSCSZOnIhff/0V48aNg5+fn9mPLxQKIZfLzXpcW2E3I7L169dDoVBg9OjRCAwMNN127Nhhes7atWtx3333YebMmRg5ciRkMhl++eUXFqu2PQaDAT///DNqa2vx8ssvs10OIeT/GzFiBKRSKfbu3Wv2Y3M4HHh5eTlskNnNiIxhmLs+RywW48svv8SXX35phYrsU0NDA9auXYu5c+ciIiKC7XIIIf+fQCDAlClTsHXrVowfP97sXdSenp4OG2R2MyIj5rF9+3bI5XLMnz+f7VLIHTAMA71eD41GA7VajebmZjQ1NaGxsbHVrbm5GWq1GhqNBlqtFnq93rSOJrE/gwcPRnR0tEU6GL28vG5phnMUdjMiI12n0+nwzTffYNq0aTQasyF6vd4UQgaDAQaDwbQWJvDXpow8Hs+03U7Ln1wu13RNYMtF7jfv9s3j8cDlck1//n0DVWJ7uFwuEhIS8NNPP6GyshIBAQFmO7aXl5fDXmJEQeZEdu/ejZKSEmzatIntUpya0WiERqMx3VrCSCQSQSAQQCgUmsKrJYTaq6UrtSXYjEYj9Hq96fGWQLs5GIltGTBgAPbu3Yvjx49j9uzZZjsunSMjdo9hGKxbtw6jR49Gnz592C7H6RiNRjQ1NaG5uRkajQYcDgdCodDUTSYSicwSLBwOp9Uecy1aQq0l2Fq28GkZwVGo2Q4+n48xY8bg4MGDuO++++Di4mKW43p5eaGpqQkajcbhNs2luQYncfjwYWRnZ+Ppp5+mDy0rUqvVqK2tRWlpKerq6sAwDLy8vBAYGAg/Pz94eHhALBZb/O+kJeAEAgHEYrFpRNZycbzBYGjzvFpSUhISExMRHh6OxMREJCUlWbRO8pdBgwZBr9fj3LlzZjtmy/J9jniejILMSXz++ecYMGAAXVNnBQzDoKGhAWVlZaioqIBarYZUKkVwcDD8/f3h6urK+vmqlilGPp9vGsH9PcyefPJJTJw4EWlpaSgqKkJaWhomTpyIf/3rX2yV7TTc3d3Rp08fpKammq1xx8vLCwAccnqRgswJnDt3DkajEQsXLmxz2omYj0qlQnFxMWpqasDj8RAQEIDg4GB4eHjY7HvfMlrj8/mmkWFSUhI2btzY5vO/++67VivoEMtITExERUUFCgoKzHI8CjJi1zZu3IjKykqMHj2a9ZGAo2psbERRURGqq6shEokQFBSEgIAASCQStkvrlLvtTffGG29YqRLnFRMTAx8fH5w6dcosx5NIJHBzc2u1GbGjoE81B1dbW4s9e/bgiSeeoBZsC9BqtSgrK0N5eTl4PB6Cg4MREBBg97sJlJWVdelx0nUcDgeJiYnIzMw0y6r1HA4HXC7XIVfAp081B/frr79CIBBgzpw5AEBBZiYMw6CmpgaFhYXQ6XQIDAxEcHCww3SDBQUFdelxYh6DBw+G0Wg0W9NHcHCwQzZ70aeaA2MYBt988w0mT55s6liiIOs6jUaDgoIC1NXVwdvbG2FhYXB1dWW7LLN6++237/j4Sy+9ZJ1CnJy5mz7kcjk0Go0ZKrMt9KnmwC5duoScnBzMnDkTRqPRIX8Tsza5XG46+R4eHg4fHx+HfF8nTJhw22XMXF1dsXTpUpw4ccLKVTmnlqaPwsJCtkuxWRRkDmznzp3w9fXF6NGjYTAYHPID11oMBgOKi4tRWVkJT09PREREOMw04u18++23OHz4MBISEhAWFoaEhAQcPnwY+fn5uOeee3Dvvfdi7dq1tK6jhcXExCA4ONgs04uO+ndFK3s4qJbtWv75z3+Cz+ejubmZphU7SaPRoKioCAaDASEhIXBzc2O7JKsZN24cxo0bd8v9+/fvx8qVK/H777+jpqYGq1atAp9PHyeWwOFwEBcXh3PnzmHGjBn077gN9I44qJMnT6KystK0VhtNLXZOU1MTbty4AQCIiopyqhC7Ez6fj/feew8LFy7Exo0bMX/+fOh0OrbLclj9+vVDQ0OD6WexKxzxc4B+hXJQ//vf/xAZGYn+/fsDgGlldNJ+CoUCJSUlcHFxQVhYmM1e0MymOXPmQCKR4PHHH4der8fGjRvt/tIDWxQSEgJvb29kZmYiOjqa7XJsDo3IHFBzczP27t2L2bNnm8KLgqxjamtrUVJSAqlUioiICAqxO7j//vvx008/Yf/+/XjsscccsiuObRwOB3379kVOTk6rLX46ylE/ByjIHNCBAwfQ0NCAWbNmme5z1B9gS2gJMU9PT4SEhND71g5TpkzBjh07cPjwYbz22ms0zWgB99xzD+rq6hx2T7GuoCBzQKdOncK0adMQGRkJ4P86legD+e7kcjmKi4vh4+OD4OBgtsuxKxMnTsT27duRmZmJN998k+1yHE5ERAQkEgmys7M7fQxH7VqkIHMwarUaW7ZsMZ0bAyjI2kupVKKoqAje3t4ICQlhuxy7NGHCBDz66KNYv349Nm/ezHY5DoXL5SI2NhY5OTlsl2JzKMgczIkTJ6BWqzFlyhTTfRRkd9fQ0ID8/Hx4eHggNDSU7XLs2rx58zB//nwsXboUaWlpbJfjUGJjY1FUVNSl9RId8XOAgszBHDp0CEOGDEFMTIzpvpaTw474A2wOWq0WeXl5kEgkiIiIoPfJDNasWYPBgwfjlVdeQUVFBdvlOIy4uDgwDIPc3NxOfT9NLRK7kJKSgtjY2Fb3tfzw0oWUt2IYBnl5eeBwOIiKiqIQMxOhUIhNmzahoaEBixYt6lKnHfk/np6eCAgI6NJ5MkdEn2wOpKqqCllZWRgxYkSr+xmGoa7F22iZpunWrRutTGFm/v7++PDDD3HkyBH88MMPbJfjMOLi4pCdnd2p0ZWrqytcXFwsUBW7KMgcyMmTJwEAw4cPb3V/S5CR1mpra1FdXe2Qq9fbinHjxuHll1/GqlWrUFRUxHY5DiE2NhY8Hg81NTUd/t6qqiqHvCaSgsyBpKSkIDo6GoGBgWyXYvM0Gg0KCwvh4+MDPz8/tstxaC+88ALc3Nzw8ssv0y9UZhAeHm7qsO0ItVoNhmHg7u5uocrYQ0HmQFJSUm6ZViRta1mzLiwsjOVKHJ+bmxv+85//QC6X49ChQ2yXY/fEYjH8/Pw6HGT19fUAYNqb0JFQkDmI8vJy5ObmYuTIkbc8RlOLrdXU1EChUCAqKsohp1ls0eTJk+Hq6op3330XBoOB7XLsXlhYWIeDTC6XAwC8vLwsURKrKMgcxO3Oj7WgIPuLXq9HQUEBfH19HfI3U1vF4XDw5ptvIisrCzt37mS7HLsXFhaGioqKDi0F1jIik0qlFqqKPRRkDiIlJQXdu3eHv7//LY/dvHCws2vZZTc8PJzlSpxPfHw8pk2bhjVr1tDCwl0UFhYGo9GIkpKSdn9PfX09JBKJQ24IS0HmIC5fvkznx+6ioaEB1dXVCAkJgUAgYLscp7R8+XLIZDLs2LGD7VLsWkBAAIRCYYemF+vr6x1yWhGgIHMINTU1OHv2LIYMGdLm4zQi+0tRUREEAgECAgLYLsVpxcTEoF+/fvj++++d/uexK7hcLrp164aqqqp2f099fb3DTqd3OMieeOIJnDhxwhK1kE66evUqEhISWi0UfDO6EBpQqVSoq6tDWFgYvR8smzZtGrKzs2kdxi7y9PREWVlZu58vl8spyFooFAqMHz8eMTExeP/991FaWmqJukgHXLlyBenp6aZtW27HmX8DLiwshIuLC3x9fdkuxekNHToUUVFR2LJlC9ul2LWAgABUV1e3e/kvGpHdZM+ePSgtLcUzzzyDHTt2ICIiApMnT8auXbtoMz2W5OTkoFu3brdtJXf2EYhKpYJSqaTRmI3gcDh49NFH8fvvv6Ouro7tcuyWv78/DAYDamtr7/pchmFoRPZ3fn5+WLJkCS5cuIDTp08jOjoajz32GIKCgrB48eJOr8xMOic7O/uWhYJvxuFwYDQanXbh1tLSUgiFQpsejRmNRuj1etPNYDA49Ah6zpw5AECt+F3Q0qHcnvNkTU1N0Ov11OzRlvLyciQlJSEpKQk8Hg9TpkzBpUuXcM8992Dt2rXmqpHcRW5u7l2DDHDOqUWdTofKykqbavDQaDRQKBSoqqpCaWkpiouLUVJSgoqKClRVVaGqqgrV1dWmm1wuR0NDA9RqNfR6Pdvlm4W3tzceeughnD59mu1S7JabmxtcXV3bNSKrq6sDh8Nx2BFZh5f71ul0+O2337Bx40YcOnQIffr0wUsvvYSHH34YHh4eAIDdu3fjySefxOLFi81eMGmtoaEBZWVldwyylu1bnHFEVlFRAQ6HA5lMxmodTU1NUKlUaGxshNFoBI/Hg1AohFgshkAgAJ/PB4/HA4fDMe1UwDAMDAaDaTTd1NQEHo8HLpdr+h57XrE/MTERzz//vM39omEvOBwOvL292xVkNTU14HA4DruuaIf/FQQGBsJoNOKhhx7CmTNn0K9fv1ueM2bMGIdNflvTMo3bnhGZMwZZeXk5fH19WblujGEY1NfXo66uDjqdDgKBAB4eHnBzc4NYLO7U+TqDwQCdTmeaftTr9XYbaKNGjQKHw8Hx48cxe/ZstsuxS1KpFEql8q7PKy8vh0QiMQ02HE2Hf/rXrl2LWbNmQSwW3/Y5np6eyM/P71JhpH1ycnIA3DnIgL9GZc42tahQKKBWqxEdHW3115bL5aiqqoJer4eHhwcCAwPNsg8Uj8czNfW0BJlOp4PRaASfz7erzVO9vLzQp08fHDt2jIKskzw8PNr1WVtWVobAwECHbXbq8E/9Y489dscQI9aVk5MDmUx2160ZuFyu0y3WWl1dDR6PZ9UT3Gq1Grm5uSgpKYFYLEZ0dDRCQ0Mtspkhj8eDSCSCWCwGl8s1NYq0/MKSlJSExMREhIeHIzExEUlJSWavoatGjx6NEydOON3Pprl4eHi0e0QWFBRkhYrYYT+/vpE25eTkICYmpl3PdbYPi8rKSvj4+Fjtt9DKykpkZ2dDr9cjMjIS4eHhVlnXjsPhmKYXWzpUn3zySUycOBFpaWkoKipCWloaJk6ciH/9618Wr6cjRo4cidDQUGRnZ7Ndil3y8PBAc3PzHS99YhgG5eXlDr1PIQWZncvJyUFcXNxdn8fj8ZwqyJqamtDc3GyVlnuDwYC8vDwUFxfDx8cHcXFxcHNzs/jr/h2XywWPx0NycjI2btzY5nO+++47JCcnW7my2+vduzeys7ORmZnJdil2qWUl+zuNympqauDm5oaQkBBrlWV1FGR2zGg0wt/fH927d7/rc50tyKqrq8Hlci0+rajT6ZCVlQWlUono6GiEhISwfp5q5cqVd3z8jTfesFIldyeRSDBs2DBaIaiTWpo37hRkhYWFkMvlDr3jAwWZHWtqasKhQ4fg4+Nz1+fyeDyn6lqsq6uDl5eXRTfO1Gq1uHLlCtRqNWJjY22mU/du6+91ZH0+a/D19cWZM2fYLsMutSfI8vPz4ePjc9fz6PaMgsyOKRQKAO3bupzH44FhGKcIM2ssx6PX63Ht2jUYjUb07NkTrq6uFnutjrrbSX1bO+nfrVs35OXlsV2GXRIKhfDy8kJjY+Ntn1NQUHDXdVjtHQWZHevIjq8tQeYoK0PcSWNjIzQajcV2wjUajbh27Rqam5vRo0cPm+viffvtt+/4+LvvvmulStonOjoaQUFBkMvlbJdil7hc7m2DzGg0orCwEBEREdYtysooyOxYS5C1Z+TB5/PBMIxTLOysVCrB4XAsdvHnjRs3oFQqcc8990AikVjkNbpiwoQJmD9/fpuPjR07FuPGjbNyRXcmk8lw7do1VFZWsl2KXRKLxVCr1W0+Vl5eDo1GQyMyYrs6MrXI4XDA4/GcYkSmVCohkUgsstpFTU0NSkpKEBkZ2amg1Ov1aG5uRkNDAxoaGtDY2Ijm5mZotVqzXrD+7bff4vDhw0hISEBYWBgSEhIwdepU+Pj42FxjRUtnaXV1NcuV2CeRSASNRtPmY/n5+eBwOAgLC7NyVdZlf+vaEJOOTC0CgEAgcJoRmSXOWel0OmRnZ8PX17fd55mMRiMaGhqgUqmgVqvBMAy4XK6pVZ7D4Zi+5nK5EAqFpltXg3jcuHGtRl/19fUYOHAgXn31VZvaC6wlyGpqaliuxD65ubndtlO2oKAAgYGBNjf9bW4UZHZMoVCYVnZoDz6ff9spCEfS1NRkkUWCr1+/DoZh7rocGPDXtWXV1dWmBV1dXFwglUohFoshFotNIQbAtDhwy/qJWq3WtIaiSCQyW+elp6cn1qxZg3nz5uHIkSMYO3asWY7bVWKxGPHx8U7xs2kJRqMRzc3NbT5WU1ODHj16WLki66OpRTsml8s71NAgEAhaLWHkiBiGgVarNfu5K5VKhYqKCoSFhUEoFN7xuTU1Nbhy5Qqqq6vh5eWFmJgYREZGwtfXF25ubqYVOFq0rMwhFovh6uoKT09PuLi4gMvlQqPRmKYdzbHk1MyZM/Hkk0/igw8+sKkO1rKyMhqRmZlKpcK1a9ccvtEDoCCzax3durxlqsqRz5O1fPCbO8jy8/PB5/MRGhp62+cYDAbk5uYiLy8P7u7u6N69O4KCgu4afG1pCTaRSAQOh2O2Jac4HA7mzJmD06dP4+DBgx2uy1Isfc2fM8rKygKAdi2YYO8oyOyYQqHoUJC1bGXiyEGmVqtNe32ZS2NjI8rKyhAZGXnbdRvVajUuXrwIhUKBuLg4REZGmmXrGB6Ph2PHjmHTpk1tPt6ZJaeGDh2Ke++9F7///nuX6zOXlvO9xHyuXbuGoKAgm7lQ35IoyOxYfX19h6YWWxoKtFqtBatil06ng06nM2uQFRUVgc/n37bBQ61W49KlSzAajejVq5fZPzhWrFhxx8c7s+TUo48+iitXrpi2AWKbTqezyz3VbEFbpwoYhsG1a9ec4vwYQEFm19zd3Ts8/83n8x06yFpGm+b6UDQajaioqEBwcHCbnWF6vR5Xr14FwzDo3bu3RbrDLLHk1Pjx41FSUoKff/65s2WZlV6vZ2XzU0dVVVWFuro6CjJi+1r2o+oIoVDo8EHW0jxhDrW1tVCr1QgODm7z8aysLKhUKvTq1cuso8CbWWLJKZFIhMcee8xm1jikEZl5Xbt2DVwut10dto6AgsxOJSUlYceOHdiwYUOHOtiEQiF0Op3DniczGAwwGAxmW4G+oqICfD6/zYufKyoqTO3Nltg4s4WllpwaNGgQtFot64sIGwwG9OrVq12LX5P2uXbtGrp162aV/fBsAQWZHWrpYFMoFGhqaupQB1vL1JcjX7PD4XDMdolBXV0d/P39b7lfr9fj+vXr8PLyavNxc7rTklP9+/fv9JJTgwcPRkZGBjIyMrpSXpdVVVXhwoULDr06uzXp9XpUV1ejd+/ebJdiNRRkdiYpKalLmya2rB7hqEHWsu+aOa6R0ul0qK2tbbN5o7i4GA0NDe3enbur2lpyauHChfD09ERxcXGnjunn54eYmBhcuHDBzNV2TMuI0NZW5bdXubm5KC8vR8+ePdkuxWooyOyMOTrY7rTIqL1ruRbJHJuItuzx9PfNORmGQUlJCQIDAy06pfh348aNQ2pqKgoLC5GamooPP/wQAG7bmt8effv2veNeVtZAQdY1fD6/1c9hZmYmfHx8bnte1xFRkNkZc3SwicVi6PV6hzxP1hJk5vh/UyqV4PF4cHNza3V/bW0t5HI56wuxuru7Y9y4cbh+/XqnR6ABAQGsLyJcWloKb29vm9xJwB40NjaaptIZhkFmZib69+9/22seHREFmZ0xRwebI58nEwgEZuvMbG5uNi3qe7Py8nK4urreMlJjw9ixY3HgwAGcP3++U98vlUpx9epVM1fVMeXl5QgMDGS1BnvG5XJN5xfz8/OhVCrRt29flquyLocMsi+//BIREREQi8UYMmSIzbQYm4M5Otgc+TyZWCxGU1PTbRdR7QitVttm11dVVRW8vb27fHyj0Yht27Z16RiDBg3CuHHjcOnSpU59v1QqRUNDA6vrb+bl5TlNm7glVFdXm345vXTpEqRSKaKioliuyrocLsh27NiBJUuWYOXKlcjIyEDfvn0xadIkVFVVsV2aWdypg23+/Pnt7mBz1PNkEokEHA4HTU1NXT5Wc3PzLaMxg8EApVJplk07165di4cffhgff/xxp4/RcoH7uXPnOvX9Li4u6N69O2sLCOv1ety4ccOpOuzMyWAwoKmpCR4eHjAajUhLS8PAgQPNdvmJvXC4/9tPPvkECxYswLx583DPPfdgw4YNcHFxwffff892aWbTVgfb4cOH8e2337b7GI56nozD4cDFxeW2W7935ng3a2hogNFoNEuQfffddwDQ5Z/NQYMGdXqfuebmZuTk5LC2YG9WVhYaGxudZgUKc1OpVAD+Ol969epVKJVKDBkyhOWqrM+hLqXXarVIT0/H8uXLTfdxuVyMHz8eqamptzxfo9G02lmV7e6tjvj7pokd1TJlplarb2lmsHeurq5oaGjo8nFu3jOshVqthl6v7/RSVE8++STKy8sB/N/q5NeuXcPkyZMBAIGBgR0ONldXV1RUVHSqnsbGRotsQtpe58+fh0gkoiDrpJbPLA8PDyQnJyM4OPiOOzQ4KocKspqaGhgMBgQEBLS6PyAgwPShcbPVq1dj1apV1irPpvB4PAgEAjQ2NjpckHl4eJhlbysOh3NLIOr1enC53E4tp9TQ0IBNmzbdcj6KYRgcOHDA9Jqff/55h/5OXF1dOz2y1uv1GDp0aKe+1xwyMjLQu3dvWmexk1qCjMvl4tKlS/jnP//JckXscLipxY5Yvnw5FAqF6dbZC0vtlUQiadW66yg8PDxQW1vb5WlTiURyy5Qdl8uFWq3u1DklNzc3nDlz5rY7FkilUpw7d67Dv1gwDNPpX0auXbtmlvOJnWE0GlFdXY2RI0ey8vqOQKVSQSQS4eLFi2AYBoMGDWK7JFY4VJD5+vqCx+OhsrKy1f2VlZWQyWS3PF8kEsHDw6PVzd505SS9m5vbHbdJt1c+Pj5gGAZ1dXVdOo6Liwu0Wm2rVn6BQACRSNTp92zgwIG3nQasqKhAfHx8h49ZVVXV6XOC+fn5iIyM7NT3dtXFixdx5coV9O/fn5XXdwRKpRLu7u5ITU1F7969HW52pb0cKsiEQiEGDBjQapkmo9GI5ORkJCYmsliZZXV2RCUUCiEQCMxyPsmWeHp6gs/no7a2tsvHUavVkMvlpvtazr915T3bsmVLm/dv3bq1U8fT6XSdCoPm5mZ4eHigX79+nXrdrkpNTcWgQYPQq1cvVl7fEahUKuj1epSUlDj0Z9zdOFSQAcCSJUvwzTffYPPmzbh27RqeeeYZNDY2Yt68eWyXZpPc3NwcbnqRw+EgICCgy5dceHp6QiKRtDqORCKBl5dXl87B/fDDDwCA0NBQnDp1CiEhIQCAzZs3d+p4J0+e7NRv4ufPn8elS5dwzz33dOp1u8JoNOL3339HVFSU07WKm5NSqURJSQnc3d2dam3Fv3OoZg8AmDNnDqqrq7FixQpUVFSgX79+OHDgwC0NIOQvbm5ukMvlaG5utuq6gZbm4+ODq1evwmg0dvqDksvlwtvbGxUVFa2uc3J1dUVpaSkGDhzYqeMuWrQICQkJWLNmDbhcLgoLC7Fs2TIMGDCgw8eqqKiAv78/EhISOvy9aWlp8PX1ZWVq8cKFC6isrDR1a5LOkcvlKCoqwpw5c5z6FwKHCzIAeO655/Dcc8+xXYZVdHU9NZFIZJpedKQgCwoKQnp6Oqqqqto8P9peMpkMV65cgV6vN3UqhoeHIy0tDc3NzZ1aH3DOnDmYM2eO6Wsul4sPPvigU/UdO3YMV69e7fC5NYZhkJubi+nTp7OyJt+JEyfQr18/uhC6CxiGQV5eHgQCQad+kXEkzhvhDqSr04KOOL3o7e0NNzc3FBUVdek44eHhUCgUrY4THh4OrVaL7OzsrpbZJQzD4ODBg5g1a1aH9/K6dOkSTp06xcp5lYaGBuzatQuJiYlOPYroKpVKhbKyMoSFhTn9WpX0U0Tg5uYGg8HgcEtWRUREoLi4uEsB7enpCZlMhtzcXNN9YrEYMTExOH/+vFm2i+ms9PR0yOVyjB49usPfm5ycjO7du2Pw4MHmL+wu9u3bB51Oh3/84x9Wf21HcvHiRdTU1GDChAlsl8I6CjJiml5sWe7GUYSHh5t+a+2KuLg45ObmQqFQmO7r06cPjEYja5tSGgwGfPnllwgMDOxwGNXW1iIjIwNTp061+ojIaDTi2LFjmD59Ovz8/Kz62o7m6NGjcHd3x4gRI9guhXUUZHaOw+GYZUrQzc0NKpXKoaYXfX19IZVKkZeX16XjxMbGQiqVIjMz03Sft7c3IiMjkZqaysrSZnv27IFQKMSCBQs6HEabN29GUVERK6tAHDlyBBkZGbj//vut/tqOpKamxrRAMFvrZNoSCjIHwDBMl1cvd3d3h8FgcLhrymJjY3Hjxo0uTZvy+Xz06tULFy5caHVN2bBhw+Dp6Yl9+/ZZdfX469evY/v27YiOju5wk0dZWRmOHDmCBQsWdPi8WlcZjUYcPHgQQ4cOZaXl35EkJSVBp9Nh1KhRbJdiEyjIHIA5us6EQiEkEkmr6TNHEB0dDZFI1OXNI1tWTThx4oTpPpFIhFGjRkGhUODw4cNdLbVdamtr8frrr0Mmk3W4M5dhGHzyyScIDw/HjBkzLFTh7SUlJeHkyZN0TWcX1dXVISUlBWFhYQgPD2e7HJtAQeYgzDEl6OHhgaamJrPsrmwrRCIRIiIicOXKlU5vdQL8NSobOXIksrOzce3aNdP9ISEhGDZsGM6dO4cjR45YdGq2rq4O//73vxEWFoY33nijw63/SUlJyMnJwezZs61+qYVOp8NXX32FkSNHok+fPlZ9bUdz6NAhGAwGhIWFtWtHeGdAQeYAzHWezMPDAxwOB/X19V0vyob06dMHOp2uzR0QOiI6OhqxsbE4dOhQq8aY3r17Y/LkyTh9+jT27NnTamsgc7l+/TreeOMNcLlcvPDCCx2+wL+kpATfffcdBgwYwMp01O7du1FZWYlnn33W6q/tSOrq6nD69GlER0fDx8eH1S14bAkFmYMwGo1dDjMOhwNPT0/U19eztmOwJbi6uiI6Otp0YXNXTJo0CVKpFL/88kur1vv+/ftj+vTpKC8vx/r161FQUNDFqv9iMBiwe/durF27FlwuF6tXr+7wflONjY14//33IRQK8e9//9ssdXVEbW0tvv32Wzz66KM0FdZFhw8fhkQigY+PD43GbkJB5gBautbMMSrz9PSE0Wh0uHNl/fr1g1Kp7HK7vEQiwYQJE1BWVmbaQ6xFjx49MGfOHPj5+eGHH37Azp07b9mJob2MRiNOnTqFZcuWYfv27RgyZAg+/vhj+Pv7d+g4Op0Ob731FsrLy/H222+z8hv8N998g5CQEDz88MNWf21HIpfLkZaWhtGjR6O2tpaC7CYOuUSVMzLX9KJAIIC7uzvkcjm8vLzMUJlt8PDwQJ8+fXD+/HnExcV1abuL4OBgTJ8+HX/88Qd4PB4mTZpkarjx8fHBQw89hIyMDKSnp+Pjjz9Gr169EBMTg+7du9/xPTUajSgoKMCFCxeQnJwMoVAIT09PrFmzplPrIer1evznP/9BeXk53njjDVbWVExOTsa+ffvw1ltv3XYfNtI+hw8fNu2mfe7cOQqym1CQOQgOhwOj0WiWa0q8vLxQVFSExsZGh5qDj4+PR3Z2Nk6fPo1x48Z16Vg9evRAU1MTkpKSoNfrMWXKFNPImMvlYuDAgejfvz8uX76M8+fP49ChQ/j111/h7e0NsViMgIAAGI1G03YzSqUSpaWlqK+vh0wmQ69evTB+/HhER0d3qr7m5masXLkSeXl5ePnllzu1IHFX1dbW4pNPPsHYsWO7/H47u/r6eqSmpmLy5Mmora0Fh8Pp0hqijoaCzEFwudwun/9p4eLiApFIhJqaGocKMqFQiCFDhiA1NRWlpaUIDg7u0vEGDBgADoeDI0eOoK6uDrNmzWrVScjj8dC3b1/07dsXjY2NuH79OiorK1FYWAiVSoWKigp4e3ujqakJISEhppFbVFRUl34hKS8vx4cffgitVos33niDlY0rGYbBxx9/DD6fjyVLllj99R1Ny2hsxIgROHLkCPz8/CAQCNguy2ZQkDmIlqmtrmxbcjMfHx+UlJR0eoV3W9Wy3FRKSgpmzJgBoVDYpePFx8fD3d0dBw4cwH//+1/Mnj27zWYMV1dX9O3bt0uv1R4pKSn4+OOPERYWhqVLl7K2+/Mvv/yCK1euYPny5TSl2EVyuRxZWVkYPXo0xGIxamtrTXvYkb9Qs4eDuDnIzMHDwwNCobBLG0jaIg6HY7qIOSUlxSzHjImJwSOPPAIvLy9s3boV+/fvt/q1eNXV1Xj//ffx6aefYty4cVi9ejVrIXbu3Dls2LABM2fOdPrtRczhwIED0Gg0GDVqFBobG1FfX4+wsDC2y7IpNCJzMOZsm/fz80NpaSnUajXEYrHZjss2Dw8PjBw5EidOnEBISAji4uK6fExvb2888cQTSElJQVJSEvLy8jBo0CAMGjTIomvhqVQq7N+/HydPnkRtbS2WLFmCoUOHsrLHGAAUFRXhrbfeQnx8PB555BFWanAkN27cQHp6OubMmQOxWIzMzEyo1WrWfkmxVRRkDoTH43Vp9Yq/8/DwQHV1NWpqahxuKqNHjx4oLi5GcnIyfH194ePj0+Vj8ng8jB49Gn369MGRI0ewa9cuJCUlYejQoRgyZEiXOiX/rqX2/fv3w8XFBWPHjsXMmTNZPaepUqnw2muvwdfXFytWrKDFbLvIYDBg165dCAsLM+1GnpOTg9DQUIea7jcHCjIH0nJuzFznyTgcDnx9fVFWVgaNRgORSNTlY9qSMWPGoLq6GklJSZg+fbrZRp3e3t544IEHMHz4cPz555/Yt28f9u3bh/j4eHTr1u2ubfhtYRgGhYWFuHjxIjIzM5GRkYGwsDDMmjULEydOZP08lFarxaefforGxkb897//dagmIbYcO3YM1dXVWLJkCTgcDtRqNYqKijB27Fi2S7M5HMaR9u3oIqVSCalUCoVCAQ8PD7bL6RSNRgMul2u2jqaW7dSFQqFDzssrFAps27YNMpkM9913H/h88/9u19jYiHPnziEnJwcXLlyAVCqFh4cHPDw8EBISAg8PD7i7u0MoFILL5cJgMKCxsRFyuRyVlZVQKpUoLi5GQUEBIiIiEBUVhb59+2LQoEE20bmm0WiwcuVKqFQqLFq0iFa2N4Pa2lp88MEHGDFiBO677z4AwJUrV7Bv3z48++yzZh3d26qOfB5TkN3EEYJMp9PBaDSadfQkl8tRWlqKbt26OeSURllZmWkK5/7777foZpMtbfj5+fkoKCiAXC5HbW0txGIx6uvrodPpEBAQgKamJjQ0NCAkJAQhISGQyWSIiYlBbGysTYRXC41GgzfffBNXr17Fu+++i379+rFdkt1jGAbffPMNqqqq8Morr5g6a3/55Rc0Nzc7zbnHjnwe09Sig+HxeNDr9WAYxmwn/D09PVFTU4OKigqHPMkcFBSE+++/H3v27MEff/yBKVOmWKxZoqUN/+ZWfKPRiKamJqjVajAMAy6XC4lEAolEwlrTRns0NzfjzTffRHZ2Nt5//31a1d5MMjMzkZ2djfnz55tCTKvVIj8/HyNHjmS5OttE7fcOhsvlgmEYs10cDcC0ioBKpWJlN2RriIiIwNSpU3H9+nXs37/fqosmc7lcuLm5wdfXF35+fvDx8YGLi4tNh5hKpcJHH32EnJwcrF69mkLMTJqbm/Hrr7+iT58+raZob9y4AYPBgNjYWBars10UZA6Ix+O1WpndHNzd3eHu7o6ysjKHWhn/ZjExMbj33ntx6dIl7Ny506wdoI6koKAAL774IuRyOT788EP06tWL7ZIcRss1iNOnT291f05ODgICAlhv6rFVFGQOiM/nm2Vbl78LCgqCTqdDdXW1WY9rS+Li4jBnzhwUFRXhp59+QnNzM9sl2ZRTp05h8eLFEIlEeOWVV8xyDR75S2FhIVJTUzFlypRWgaXX65GXl0ejsTugIHNALc0K5h6ViUQi+Pn5oaqqyqF2kf67yMhIPProo2hqasLGjRtRVVXFdkmsMxqN2LJlC9555x0MHDgQn3zyCS1aa0YGgwE7d+5ESEgIhg4d2uqxgoIC6HQ6CrI7oCBzQBwOx+wXR7fw9/cHn89HWVmZ2Y9tSwIDA/Hggw/CxcUF33//PS5evMh2SaypqqrC+++/j1OnTuHxxx/Ha6+95pDdq2w6ceIEKioq8MADD9zSNZudnQ0fHx+zXLTvqCjIHFTLeTJzn8/icrkIDAyEUqmESqUy67FtjZeXFx555BF0794de/bswb59+6DRaNguy2oYhsH+/fuxaNEi5OXl4emnn8ZDDz1k000o9qiiogJnz57FuHHjbllBR6fTITc3l85D3gW13zsoPp8PDocDnU5n9hU5PD09UVtbi7KyMsTGxjr0B5tAIMD06dMRFRWFvXv3IicnB9OmTev0PmH2ory8HF988QUuXryIe++9F/PmzYOLiwvbZTkcrVaLLVu2QCAQtLlnW1ZWFvR6PXr06MFCdfaDgsyB8fl86PV6iywtFRwcjNzcXFRVVSEgIMDsx7c1ffr0QVhYGH799Vfs3r0bsbGxGDNmjN1eOH87TU1N2L17N06fPg2NRoN3333XKtvPOCOGYfDzzz9DoVDghRdeaHNLoYsXL6J79+4O93NmbhRkDkwgEECn00Gv15t96SWxWAwfHx9UVFRAKpU61Or4t+Pp6YnHHnsMmZmZOHnyJD799FMMGzYMw4cPt/t1KHU6Hfbv349du3ZBq9XigQcewLRp05zi75Utp0+fRmZmJh555BH4+fnd8nhpaSkqKysxbNgwFqqzLxRkDozL5Zp2jrbEGoIymQz19fUoLCx0+CnGFlwuF/Hx8bjnnnuQkpKCkydP4uzZsxgxYgQGDRrU5Y06rU2j0eD48ePYs2cPKisrMWHCBMyePRve3t5sl+bQSkpK8Ntvv2Ho0KG3vZg8PT0d3t7eiIiIsG5xdoiCzMEJBAJoNBqzLlnVgsvlIjw8HDk5OaioqEBgYKBZj2/LxGIxJkyYgMGDB+PEiRPYt28fDh8+jISEBCQmJtr8VFBlZSUOHTqE5ORkeHt7o3///pg6dSqCgoLYLs3hNTc3Y8uWLQgMDMTUqVPbfE59fT2uX7+OCRMmWLk6+0RB5uD4fD7UajW0Wq1Fpr9cXV0hk8lQUVEBDw8Pp9u+QyqVYtq0aRg5ciT+/PNPpKam4vTp04iKisKAAQMQGxtrM/tyabVaZGRk4MyZMzh58iRcXFwwbtw4TJgwga4JsxKGYbBjxw6o1WosXLjwtjMl58+fh0QioSaPdqIgc3AcDgd8Pt8i3YstZDIZFAoFCgoK0KNHD4uuHm+rpFIppkyZgrFjx5r2C/v+++9NiwTfc8896Natm0WmeO+ksbER58+fx9mzZ3HhwgV4eHggODgYTz/9NIYNG2b35/bszYkTJ3Dt2jXMmzfvtnvSaTQaXL58GQMGDLD6z4u9onfJCQiFQjQ2NkKn01lkCxAOh4PIyEhcu3YNJSUlDrlvWXuJxWIkJCRgyJAhKCsrQ2ZmJgoKCnDy5Enw+XzExMQgKioKYWFhCA4ONvuFxXK5HHl5ecjLy8Ply5ehUChQXV2Nbt264Z///CcGDx7sVFPAtiQ/Px8HDhzAmDFj0L1799s+79KlSzAajdQt2gEUZE6Ax+OBz+dDo9FYbC8rkUiEkJAQFBUVQSqVOv3iphwOB8HBwQgODgbDMKisrERWVhaysrJw+vRp/PrrrwAAHx8fREVFwcXFBV5eXpBKpfD09DRttMnj8Uw3g8EAtVqNxsZGNDY2QqlUora2FgqFAvn5+SgpKUFDQwMMBgO8vb3Ro0cPjB07Fn379qXmDZapVCps3boVERERmDhx4m2fZzQacf78ecTFxTndNH1XUJA5CYFAgMbGRot1MAKAr68vFAoFCgsL0b17d7vr4LOUlm1wZDIZRo8eDaPRiKqqKpSUlKC4uBj19fW4cOECFAqFaVmxgIAAlJaWmo7h6uoKlUoFFxcXNDQ0APhruTC5XI5u3brBzc0N48aNQ1hYGKKiom47bUWsz2g0Ytu2bWAYBg8//PAdp95zc3OhUqkwYMAAK1Zo/yjInIRAIACPx4NarbboNunh4eG4evUq8vLy0L17d6doye8oLpdrCraBAwea7mcYBk1NTVAoFGhsbIRarYbBYIDBYDDtZCASieDq6gpXV1e4ubnZ/L5lzq5lma+6ujo88sgjcHd3v+PzL168iNDQ0DavKyO3R0HmRMRiscVHZXw+H926dUNWVhaKi4ud+nxZR3E4HFNIEcdw9OhRpKam4h//+AeioqLu+Nz8/HyUlpZixowZVqrOcThfe5kTu3lUZkmurq4ICwtDVVUVamtrLfpahNiqtLQ0JCcnm643vJvU1FQEBQXRL3+dQEHmZMRiMfR6PfR6vUVfx8/PD76+vigsLERjY6NFX4sQW3Px4kXs3bsXw4YNw6hRo+76/Ly8PFRVVSExMdEK1TkeCjInY61RGQCEhYVBJBIhNzfXoTfiJORm2dnZ2LlzJ/r374/Jkyff9RymwWDAn3/+iaioKISGhlqpSsdCQeaErDUq43K5pl1tc3Nzzb43GiG2prCwENu2bUNcXBz++c9/tqsRp6VjlRYH7jwKMidkzVGZQCBAbGws1Go1cnNzTd13hDiaiooK/PDDDwgJCcGDDz7YrhVumpubcebMGfTq1Yt2gO4CCjInJRKJoNVqTdctWZKLiwu6deuG+vp63Lhxw+KvR4i11dXVYePGjfD29sZjjz3W7q7g1NRUAEBCQoIly3N4FGROqmXViJaLay3N09MT3bp1Q3V1NQoLC63ymoRYg1KpxPfffw+xWIy5c+e2e/3KmpoaXLlyBUOGDDH7UmXOhoLMibm6usJgMKC5udkqr+fr64uIiAiUl5e3WrWCEHvV1NSETZs2wWg0Yt68eR26BvDEiROQSqW33Y+MtB8FmRPj8/kQi8VoamqyWiOGTCZDSEgIiouLUVlZaZXXJMQSGhoaTEtPzZs3D56enu3+3ry8PJSUlGDkyJE2s82PPaOVPZycq6srtFotGhoarLYZZEhICPR6PfLz88Hn8+kkN7E7tbW1+PHHH6HX6/HYY491aEmplnb7sLAw2v3ZTGhE5uRalkXSarVWvdYrIiICvr6+uH79Ourq6qz2uoR0VWlpKb799lvweDz861//QkBAQIe+/9y5c1Cr1Rg5cqSFKnQ+FGQEIpHItDq+Ndvju3XrBk9PT2RlZdFSVsQu5ObmYuPGjfDx8cH8+fM7NJ0I/DWSS09PR79+/WhrHTOiICMAADc3NxiNRqs1fgB/jQZjY2Ph7e2Na9eu0TkzYtPOnz+Pn376CVFRUXjiiSfg4uLSoe83Go1ITk6Gp6cnbdNiZnSOjAD4a/NNiUSCpqYmiEQiq52A5nA4iIuLA5/PR25uLgwGA4KCgqzy2oS0B8MwSElJQXJyMgYOHIipU6e262Lnv7tw4QJqamowY8YMavAwMwoyYiKRSKDRaNDQ0GDVHZ45HA6io6PB5/Nx48YN6PV6WgGc2ASj0Yg//vgDZ86cwZgxYzBq1KhO7f9WXV2N06dPY8CAAZDJZBao1LlRkBETDocDNzc31NfXo6mpqcNTJ10VEREBPp+PgoIC6PX6u+7fRIgl6fV6/Pzzz7h27Rruv//+Tk8H6nQ6JCUlwcfHB4MGDTJzlQSgICN/IxAIwOfzIZfLIRAIIBAIrPr6ISEh4PP5uH79OnQ6HWJjY2kHZGJ1TU1N2LlzJ4qKivDQQw8hLi6u08f6888/0dDQgNmzZ3dqSpLcHQUZuYWnpye0Wi3q6urg7+9v9SCRyWTg8XjIycmBRqPBPffcY7EdrQn5u5KSEuzbtw8ajQZPPPFEl6a58/LycPXqVYwZM6bDHY6k/ejXA3ILDocDb29v6PV61NfXs1KDn58fevXqhcbGRpw/f96q3ZTEOTEMg9TUVGzatAl8Pr/LIdbQ0IBjx46hW7du6NGjhxkrJX9HQUbaJBAI4OnpicbGRjQ1NbFSg1QqRb9+/cDhcJCRkUEXThOLaWpqwvbt25GUlISEhATMmzevSw1PDMPg8OHDEAgEGD16tPkKJW2i+RpyW66urtBoNJDL5RAKhaxM70kkEvTr1w9ZWVm4ePEiIiMjERYWRufNiNkUFhbil19+gcFgwMMPP4zo6OguH/PcuXOoqanBlClT2r0aPuk8CjJyRzefL/Pz82MlQPh8Pnr27ImCggLk5+dDpVKhe/fudN6MdEnL9WHHjx9HWFgYZsyYAXd39y4ft6CgAOnp6Rg8eDBdE2kl9ElA7ojL5cLHxwdVVVVQKBSsnbDmcDiIjIyEh4cHrl27hnPnzqFHjx5Wvd6NOI6Ghgbs3r0b+fn5GDVqFEaOHGmWX9LkcjmSk5MRGRmJ+Ph4M1RK2oPOkZG7EggEkEqlaGhoYL3pwsfHBwMHDoRAIEBGRgZu3Lhh1fUhif3Ly8vDV199herqajz++OOdvsj577RaLQ4cOAB3d3eMHTvWDJWS9rKLICsoKMD8+fMRGRkJiUSCbt26YeXKlbes1n7x4kWMGDECYrEYoaGh+OCDD1iq2PG4ublBIpGgtrYWer2e1VrEYjHi4+MRGRmJwsJCnDt3jrWGFGI/tFot/vjjD/z2228IDAzEwoULzbaNitFoxPHjx6FWqzFp0iSrX3/p7OxiajErKwtGoxFfffUVoqOjcfnyZSxYsACNjY346KOPAPy13fjEiRMxfvx4bNiwAZcuXcKTTz4JT09PPPXUUyz/HzgGLy8vaLVaVFVVQSaTsXpxJ4fDQUREBLy9vXH16lWcOXMGMTExCA4OZq0mYruys7Nx+PBhaDQajBgxAgMHDjTr+d6TJ0+isLAQkydPpuluFnAYO52X+fDDD7F+/XrcuHEDALB+/Xq8/vrrqKiogFAoBAAsW7YMe/bsQVZWVruOqVQqIZVKoVAorLbJpL3R6/UoLy+HQCBAQECATXQPGgwGXL9+HaWlpfDx8UGPHj1MPwPEucnlchw8eBA3btxAdHQ0Jk6caPbzvBkZGTh37hxGjx6N2NhYsx7bmXXk89guphbbolAoWu3nk5qaipEjR7b6AJs0aRKys7Mhl8vbPIZGo4FSqWx1I3fG5/Ph7+8PjUaD6upqtssB8NfK/XFxcejbty9UKhX+/PNPFBUV0bkzJ6bX63HixAl8/fXXqK2txQMPPIDZs2ebPcSys7Nx7tw5DBw4kEKMRXYZZNevX8e6deuwcOFC030VFRW37NTa8nVFRUWbx1m9ejWkUqnpFhoaarmiHYhIJIKfnx+amppsakNMHx8fDBkyBD4+Prh8+TJOnTrF2sokhD15eXn4+uuvcerUKQwZMgRPPfWURUKmuLgYJ06cQI8ePahDkWWsBtmyZcvA4XDuePv7tGBpaSnuvfdezJo1CwsWLOjS6y9fvhwKhcJ0Ky4u7tLxnImLiwt8fX2hUqmgUCjYLsdEIBCgT58+SExMBMMwOHXqFC5dunRLYxBxPAqFArt27cKOHTvg5eWFBQsWYPTo0RZpvKipqUFSUhLCwsIwbNgwsx+fdAyrzR5Lly7F3Llz7/icm7fyKCsrw5gxYzB06FB8/fXXrZ4nk8lu2WG45evb7f8jEonoqvsucHNzg16vh1wuB4/Hg5ubG9slmXh5eWHYsGEoLi5GdnY2KioqEBsbS6uCOCC1Wo20tDRcu3YNOp0O06dPxz333GOx16uvr0dycjJ8fX0xbtw4WtHeBrAaZH5+fvDz82vXc0tLSzFmzBgMGDAAGzduvOWHJzExEa+//jp0Op3pN7CkpCTExcXBy8vL7LWTv3h6esJgMKCmpgZcLtfqe5jdCYfDQVhYGGQyGbKzs3HlyhWUlJQgLi4Ovr6+bJdHukiv1yMjIwOnTp2CwWDAoEGDkJCQYNFGn/r6euzbtw8SiQSTJk2i1WVshF10LZaWlmL06NEIDw/H5s2bW20T3jLaUigUiIuLw8SJE/Hqq6/i8uXLePLJJ7F27dp2t99T12LnVVVVobm5GQEBARCLxWyX06b6+npkZWWhpqYGvr6+6N69O22tYYe0Wi0yMzNx5swZ8Hg8REZGYvjw4RafEVAoFNi3bx9EIhGmTp1qsz/njqIjn8d2EWSbNm3CvHnz2nzs5vIvXryIRYsW4ezZs/D19cXzzz+PV199td2vQ0HWeQzDoKKiAmq1GoGBgTb9j7yiogJZWVloaGhAYGAgoqOj6dofO6DRaJCeno5z585Bq9Wid+/eGDJkiFV+GVEqldi3bx8EAgGmTp0KiURi8dd0dg4XZNZCQdY1LWGm0Wggk8lsOswYhkFpaSlyc3OhUqkQEBCA2NhYmoa2Qc3NzTh37hzS09NhNBrRp08fDBkyxCwL/LaHQqFAcnIyDAYDpk6dalPT546MgqyTKMi6riXMmpqaEBgYaPP/6P8eaH5+foiIiIBMJqOmEJbV1NTg/PnzKCgoQENDA/r3749BgwbB1dXVqjUcPHgQEokEkydPppGYFVGQdRIFmXncHGYBAQE21c14Oy01X79+HbW1tZBIJIiMjER4eDh1tlqR0WhEXl4eMjIyUFxcDFdXV/Tv3x99+/a1eoiUl5fj8OHD8PLywoQJE+jnwMooyDqJgsx8GIZBVVWVaZRjT+eg6uvrkZ+fj5KSEjAMg+DgYERGRrZaSYaYV319Pa5evYrc3FxUV1cjODgY/fv3R0xMTKvmLmspKCjAsWPHEBgYiHHjxlF3IgsoyDqJgsz8ampqUF9fD29vb7sLAq1Wi6KiIuTn56OxsRFSqRQREREIDg6mtRzNQKvVIicnB1euXEFpaSmEQiF69uyJnj17wt/fn7W6cnJy8OeffyIiIgKjRo2i68RYQkHWSRRkllFXV4e6ujp4enra5fVbLaPLGzduoKqqCnq9Hv7+/ggJCaFQ6yCNRoP8/Hxcv34dlZWVUCgUCA8PR8+ePREdHc3qyIdhGJw/fx7Z2dmIiIhAQkICnSdlEQVZJ1GQWY5CoUBNTQ0kEgnrW8B0hVqtRmlpKUpKSlBTUwMAFGp30dTUhLy8PFy/fh3FxcUwGo0ICAhAXFwcYmJirNZ9eCc6nQ4pKSkoKirCwIED0atXL7ZLcnoUZJ1EQWZZDQ0NqKiogEAgQHBwsN2fd7g51Kqrq8HhcBAQEICAgADIZDKb+IBmA8MwqKmpQWFhIW7cuIGysjJwOByEhISgW7du6Natm029N42NjUhOToZKpcLIkSNp8XAbQUHWSRRklqfRaFBaWgoACAoKsulrzTqiJdQqKipQXl4Oo9EIFxcXyGQy+Pr6ws/Pz6pt49ak0+lQXV2NiooKVFRUoLa2FnK5HHw+H5GRkQgLC0NUVJRNtq5XVVXhyJEj4PP5GDduHF1HaEMoyDqJgsw6DAYDSktLoVarERAQYFcdje2h1+tNH+w1NTWmrW4kEgn8/Pzg4+MDT09PeHp62mVLd319vSm0KisrUVNTA4ZhTJutBgUFISgoCIGBgax0HLZXVlYWLly4AA8PD4wZM8ZhfqlyFBRknURBZj0Mw6CyshL19fWQSqUOfQGyVqtFdXU1ampqUF1dDblcDp1OB+Cv7XBaQs3T09O0Nx7b5xANBsMtm86qVCoUFxdDo9EAALy9vSGTyUxTqd7e3nbxd6jVapGamorCwkL06NEDAwYMYP39JreiIOskCjLrUygUqKiogFAodJpmCYZhoFKpIJfLUV9fb7o1NjaCy+XCYDBAKBTC1dUVEomk1c3FxcX032KxuMMfwAaDAXq9HjqdDmq1Gs3NzVCr1aagUigUUCqVaGxsNH0Pj8eDh4cHfH194eHhAZlMBn9/f7scTdbW1uLEiRPQaDQYOnQowsLC2C6J3AYFWSdRkLGj5byZXq+HTCZz2vdeq9Wivr4eSqUSDQ0NaG5uRnNzM5qamkyBc/M/V4PBAB6PB4lEAp1OZ9qMtiXcOBwO+Hw+mpubodPpoNfrTd8vFovR1NQEAODz+eByufDw8LjlJpVK4eLiYhcjrbvJysrCuXPn4OXlhZEjR9pUwwm5VUc+j+27bYw4BJFIhIiICJSXl6OsrAwqlQoymcymz69YglAohL+//20vBmYYBmq12hRsGo3GFFAGgwEMw8BoNJqeazQawePxTIF2800oFEIkEkEsFkMsFtt9B+mdNDc3Iy0tDdXV1YiNjcWAAQOc7mfL0TnuTy+xK1wuF8HBwaZGgqamJgQFBTlsp19ncDgc07QiaZ/8/HycPXsWPB4Pw4YNQ3BwMNslEQugICM2xdPTE66urigrK0NhYSG8vb3h7+9PJ+NJh6jVapw+fRrFxcWIiIjAoEGD7PKcHmkfCjJicwQCAcLDw1FXV4fKyko0NDQgKCjI5reEIbYhPz8f586dAwCMHDmSGjqcAAUZsVne3t5wc3NDaWkpbty4AW9vbwQEBND5DdImpVKJc+fOoaKiAhEREYiPj6drw5wEBRmxaUKhEBEREabRmVKpRGBgoMNdRE06T6/X48qVK7h27RpcXFwwevRoBAUFsV0WsSIKMmLzOBwOfHx84OHhgfLychQXF6OmpsYudqAmllVSUoL09HQ0NzejZ8+euOeee2jE7oQoyIjdEAgECAsLg0qlQnl5Oa5fvw5PT0/IZDKnuJCa/J/a2lpcuHDBtNfduHHj7GIncmIZFGTE7ri7u8PNzQ1yuRyVlZXIzs6Gr68v/P396bdxB6dSqXDx4kUUFxdDKpUiISGBphEJBRmxTxwOB97e3vD09ER1dTWqq6tRV1eHgIAA+Pj4OMRKFOT/NDc348qVK8jLy4NEIsGQIUMQERFBf88EAAUZsXNcLhcBAQHw9vZGZWUlysrKUFNTA39/f3h5edEHnZ1rampCdnY2CgoKAAB9+/ZFTEwMjbxJKxRkxCEIBAKEhITA19cXFRUVKCoqQmVlJXx8fODj40MffHamoaEB165dQ2FhIfh8PmJjYxEbGwuBQMB2acQGUZARhyIWixEREYHm5mZUVlaitLQU5eXl8PHxgb+/PzWF2Di5XI7s7GwUFxdDJBKhd+/e6Natm0OvBUm6jn46iEOSSCSIiIhAcHCw6RxaVVUVPD094e/vTx1uNsRoNKK0tBTXr19HTU0N3N3dER8fj4iICBpJk3ahICMOTSAQICgoCDKZDHV1daiqqkJOTg5cXFzoPBrLmpubcePGDeTn50OtVsPf3x9Dhw5FUFAQ/Z2QDqEgI06By+XC19cXvr6+UCqVqKqqQkFBAUpKSuDt7Q0fHx9aVd4KDAYDysvLUVRUhKqqKgBAeHg4unXr5rT70JGuoyAjTqdl00i1Wo2qqirU1NSgrKwMEokE3t7e8Pb2plAzI4ZhUFtbi6KiIpSWlkKn08Hb2xt9+/ZFSEgINXCQLqMgI05LLBYjLCwMISEhUCgUqKurQ0VFBUpLSynUuohhGMjlcpSVlaGkpARNTU1wdXVFt27dEBYWRucoiVlRkBGnx+Vy4eXlBS8vLxiNRigUCsjlclOoicViU6jR2o63ZzAYUF1djfLyclRUVECj0UAkEiEwMBChoaHw8fFhu0TioCjICLnJzaHGMIxppNZysbVYLIa7uzs8PDzg7u7u1O38Le9PTU0NampqUF9fD7VaDTc3N4SGhiIwMBDe3t7UuEEsjoKMkNvgcDjw9PSEp6cnGIaBUqlEXV0d5HI5ysvLAeCWYHPkERvDMFCpVKitrUV1dTVqa2uh0+nA4/Hg7e2NmJgY+Pv7w93dne1SiZOhICOkHTgcDqRSKaRSKSIjI6HVaqFSqaBUKqFSqVBTUwOGYcDn8+Hu7m4KN1dXV7u8FophGDQ1NUGhUKC+vh719fVQKBQA/ppC9PLyQlRUFPz8/ODp6Qkul8tyxcSZUZAR0glCodC0/BXw14d7Q0ODKdhKS0tRVFQEg8EAFxcXiEQiSCQSiMViiMVi03+zHQAGgwGNjY1oaGi45U8OhwOdTgcXFxdIpVLExsaaRqj2GM7EcVGQEWIGPB7PNGID/m9E09DQgKamJjQ1NUEul0OtVsNoNJq+TygUmoJNJBJBKBRCIBCAz+eDx+OBy+Xe8uedws9gMECv10On00Gv17e6abVaaDQaqNVqqNVqaDQaNDU1mb5XJBLB1dUVnp6eCA4ONo0qnfk8ILEPFGSEWACHw4GrqytcXV1veUyr1UKtVqO5udkUKs3NzVAqlWhqagLDMOBwODAYDLc9dkuw6fV6MAwDAK3+uy0ikQg8Hs90Xs/Pzw8SicRUJ13PRewVBRkhViYUCiEUCttcyYJhGNOIymg0wmg0wmAwtPnnzSO7lnBrGcm1jOpabjwej7oHicOiICPEhnA4HAgEAhodEdIB1GpECCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErtldkGk0GvTr1w8cDgeZmZmtHrt48SJGjBgBsViM0NBQfPDBB+wUSQghxGrsLsheeeUVBAUF3XK/UqnExIkTER4ejvT0dHz44Yd466238PXXX7NQJSGEEGvhs11AR/zxxx84dOgQfv75Z/zxxx+tHtu6dSu0Wi2+//57CIVC9OzZE5mZmfjkk0/w1FNPsVQxIYQQS7ObEVllZSUWLFiAH3/8ES4uLrc8npqaipEjR0IoFJrumzRpErKzsyGXy9s8pkajgVKpbHUjhBBiX+wiyBiGwdy5c/H0009j4MCBbT6noqICAQEBre5r+bqioqLN71m9ejWkUqnpFhoaat7CCSGEWByrQbZs2TJwOJw73rKysrBu3TqoVCosX77crK+/fPlyKBQK0624uNisxyeEEGJ5rJ4jW7p0KebOnXvH50RFReHIkSNITU2FSCRq9djAgQPxyCOPYPPmzZDJZKisrGz1eMvXMpmszWOLRKJbjkkIIcS+sBpkfn5+8PPzu+vzPv/8c7z77rumr8vKyjBp0iTs2LEDQ4YMAQAkJibi9ddfh06ng0AgAAAkJSUhLi4OXl5elvkfIIQQwjq76FoMCwtr9bWbmxsAoFu3bggJCQEAPPzww1i1ahXmz5+PV199FZcvX8Znn32GtWvXWr1eQggh1mMXQdYeUqkUhw4dwqJFizBgwAD4+vpixYoV1HpPCCEOjsMwDMN2EbZCqVRCKpVCoVDAw8OD7XIIIcRpdeTz2C7a7wkhhJDboSAjhBBi1yjICCGE2DWHafYwh5bThbRUFSGEsKvlc7g9bRwUZDdRqVQAQEtVEUKIjVCpVJBKpXd8DnUt3sRoNKKsrAzu7u7gcDi3PK5UKhEaGori4mLqarwLeq/ah96n9qH3qX0c6X1iGAYqlQpBQUHgcu98FoxGZDfhcrmmC6zvxMPDw+5/SKyF3qv2ofepfeh9ah9HeZ/uNhJrQc0ehBBC7BoFGSGEELtGQdYBIpEIK1eupBXz24Heq/ah96l96H1qH2d9n6jZgxBCiF2jERkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBVkHaTQa9OvXDxwOB5mZma0eu3jxIkaMGAGxWIzQ0FB88MEH7BTJkoKCAsyfPx+RkZGQSCTo1q0bVq5cCa1W2+p5zv4+tfjyyy8REREBsViMIUOG4MyZM2yXxKrVq1dj0KBBcHd3h7+/P6ZPn47s7OxWz1Gr1Vi0aBF8fHzg5uaGmTNnorKykqWKbcOaNWvA4XDw0ksvme5ztveJgqyDXnnlFQQFBd1yv1KpxMSJExEeHo709HR8+OGHeOutt/D111+zUCU7srKyYDQa8dVXX+HKlStYu3YtNmzYgNdee830HHqf/rJjxw4sWbIEK1euREZGBvr27YtJkyahqqqK7dJYc/z4cSxatAhpaWlISkqCTqfDxIkT0djYaHrO4sWL8fvvv2Pnzp04fvw4ysrKMGPGDBarZtfZs2fx1VdfoU+fPq3ud7r3iSHttn//fqZ79+7MlStXGADM+fPnTY/997//Zby8vBiNRmO679VXX2Xi4uJYqNR2fPDBB0xkZKTpa3qf/jJ48GBm0aJFpq8NBgMTFBTErF69msWqbEtVVRUDgDl+/DjDMAxTX1/PCAQCZufOnabnXLt2jQHApKamslUma1QqFRMTE8MkJSUxo0aNYl588UWGYZzzfaIRWTtVVlZiwYIF+PHHH+Hi4nLL46mpqRg5ciSEQqHpvkmTJiE7OxtyudyapdoUhUIBb29v09f0PgFarRbp6ekYP3686T4ul4vx48cjNTWVxcpsi0KhAADTz096ejp0Ol2r96179+4ICwtzyvdt0aJFmDp1aqv3A3DO94mCrB0YhsHcuXPx9NNPY+DAgW0+p6KiAgEBAa3ua/m6oqLC4jXaouvXr2PdunVYuHCh6T56n4CamhoYDIY23wdneQ/uxmg04qWXXsKwYcPQq1cvAH/9fAiFQnh6erZ6rjO+b9u3b0dGRgZWr159y2PO+D45dZAtW7YMHA7njresrCysW7cOKpUKy5cvZ7tkVrT3fbpZaWkp7r33XsyaNQsLFixgqXJirxYtWoTLly9j+/btbJdic4qLi/Hiiy9i69atEIvFbJdjE5x6G5elS5di7ty5d3xOVFQUjhw5gtTU1FvWLxs4cCAeeeQRbN68GTKZ7JauoJavZTKZWeu2tva+Ty3KysowZswYDB069JYmDkd+n9rL19cXPB6vzffBWd6DO3nuueewd+9enDhxotW2SjKZDFqtFvX19a1GG872vqWnp6Oqqgrx8fGm+wwGA06cOIEvvvgCBw8edL73ie2TdPagsLCQuXTpkul28OBBBgCza9cupri4mGGY/2ti0Gq1pu9bvny50zUxlJSUMDExMcyDDz7I6PX6Wx6n9+kvgwcPZp577jnT1waDgQkODnbqZg+j0cgsWrSICQoKYnJycm55vKWJYdeuXab7srKyHLqJoS1KpbLV59GlS5eYgQMHMo8++ihz6dIlp3yfKMg6IT8//5auxfr6eiYgIIB57LHHmMuXLzPbt29nXFxcmK+++oq9Qq2spKSEiY6OZsaNG8eUlJQw5eXlplsLep/+sn37dkYkEjGbNm1irl69yjz11FOMp6cnU1FRwXZprHnmmWcYqVTKHDt2rNXPTlNTk+k5Tz/9NBMWFsYcOXKEOXfuHJOYmMgkJiayWLVtuLlrkWGc732iIOuEtoKMYRjmwoULzPDhwxmRSMQEBwcza9asYadAlmzcuJEB0ObtZs7+PrVYt24dExYWxgiFQmbw4MFMWloa2yWx6nY/Oxs3bjQ9p7m5mXn22WcZLy8vxsXFhfnnP//Z6hclZ/X3IHO294m2cSGEEGLXnLprkRBCiP2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIcRDV1dWQyWR4//33TfedOnUKQqEQycnJLFZGiGXRWouEOJD9+/dj+vTpOHXqFOLi4tCvXz/84x//wCeffMJ2aYRYDAUZIQ5m0aJFOHz4MAYOHIhLly7h7Nmzt2wKS4gjoSAjxME0NzejV69eKC4uRnp6Onr37s12SYRYFJ0jI8TB5OXloaysDEajEQUFBWyXQ4jF0YiMEAei1WoxePBg9OvXD3Fxcfj0009x6dIl+Pv7s10aIRZDQUaIA3n55Zexa9cuXLhwAW5ubhg1ahSkUin27t3LdmmEWAxNLRLiII4dO4ZPP/0UP/74Izw8PMDlcvHjjz8iJSUF69evZ7s8QiyGRmSEEELsGo3ICCGE2DUKMkIIIXaNgowQQohdoyAjhBBi1yjICCGE2DUKMkIIIXaNgowQQohdoyAjhBBi1yjICCGE2DUKMkIIIXaNgowQQohdoyAjhBBi1/4ftxxTLMPJCl0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "import rebound, rebound.data\n",
    "%matplotlib inline\n",
    "sim = rebound.Simulation()\n",
    "rebound.data.add_outer_solar_system(sim) # add some particles for testing\n",
    "for i in range(1,sim.N):\n",
    "    sim.particles[i].m *= 50.\n",
    "sim.integrator = \"WHFast\" # This will end badly!\n",
    "sim.dt = sim.particles[1].P * 0.002 # Timestep a small fraction of innermost planet's period\n",
    "sim.move_to_com()\n",
    "E0 = sim.energy() # Calculate initial energy \n",
    "rebound.OrbitPlot(sim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us integrate this system for a few hundred years. An instability will occur. We can then measure the energy error, which is a good estimate as to how accurate the integration was."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative energy error with WHFast: 42.011050\n"
     ]
    }
   ],
   "source": [
    "sim.integrate(600*2.*math.pi)\n",
    "E1 = sim.energy()\n",
    "print(\"Relative energy error with WHFast: %f\"%((E0-E1)/E0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An energy error that large means we basically go it wrong completely. Let's try this again but use MERCURIUS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative energy error with MERCURIUS: -1.629869e-07\n"
     ]
    }
   ],
   "source": [
    "sim = rebound.Simulation()\n",
    "rebound.data.add_outer_solar_system(sim) # add some particles for testing\n",
    "for i in range(1,sim.N):\n",
    "    sim.particles[i].m *= 50.\n",
    "sim.integrator = \"mercurius\" \n",
    "sim.dt = sim.particles[1].P * 0.002 # Timestep a small fraction of innermost planet's period\n",
    "sim.move_to_com()\n",
    "E0 = sim.energy() # Calculate initial energy \n",
    "sim.integrate(600*2.*math.pi)\n",
    "E1 = sim.energy()\n",
    "print(\"Relative energy error with MERCURIUS: %e\"%((E1-E0)/E0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, MERCURIUS is able to integrate this system with much better accuracy. When a close encounter occurs, it automatically (and smoothly!) switches to the IAS15 integrator. When there is no close encounter, you still get all the benefits in terms of speed an accuracy from a symplectic integrator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a few options to adjust MERCURIUS. First of all, because it uses IAS15 internally, you may want to set a minimal timestep for IAS15. This ensures that IAS15 never stalls while it is tries to resolve one very close encounter and can be done with the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sets the minimal timestep to a fraction of the global timestep\n",
    "sim.ri_ias15.min_dt = 1e-4 * sim.dt  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You also may want to change the critical distance at which MERCURIUS switches over from pure WHFast to IAS15. This is expressed in units of Hill radii. The default is 3 Hill radii, in the following we change it to 5 Hill radii:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.ri_mercurius.hillfac = 5"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}