{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# User Defined Rebound Collision Resolutions\n", "\n", "In the [CloseEncounter](https://rebound.readthedocs.io/en/latest/ipython/CloseEncounters.html) example, we discuss methods for resolving collisions in REBOUND through exceptions and the use of the `sim.collision_resolve = \"merge\"` method.\n", "\n", "Using the same 3-Body setup, let us explore how to define and implement the same collision resolution function in python and pass it to the `sim.collision_resolve` function pointer." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import rebound\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "def setupSimulation():\n", " ''' Setup the 3-Body scenario'''\n", " sim = rebound.Simulation()\n", " sim.integrator = \"ias15\" # IAS15 is the default integrator, so we don't need this line\n", " sim.add(m=1.)\n", " sim.add(m=1e-3, a=1., r=np.sqrt(1e-3/3.)) # we now set collision radii!\n", " sim.add(m=5e-3, a=1.25, r=1.25*np.sqrt(5e-3/3.))\n", " sim.move_to_com()\n", " return sim" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To reiterate the previous method, let's run the built-in `merge` collision resolution method" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Particles in the simulation at t= 0.0: 3\n", "System Mass: [1.0, 0.001, 0.005]\n", "Particles in the simulation at t= 100.0: 2\n", "System Mass: [1.0, 0.006]\n" ] } ], "source": [ "sim = setupSimulation()\n", "sim.collision = \"direct\"\n", "sim.collision_resolve = \"merge\" # Built in function\n", "\n", "print(\"Particles in the simulation at t=%6.1f: %d\"%(sim.t,sim.N))\n", "print(\"System Mass: {}\".format([p.m for p in sim.particles]))\n", "sim.integrate(100.)\n", "print(\"Particles in the simulation at t=%6.1f: %d\"%(sim.t,sim.N))\n", "print(\"System Mass: {}\".format([p.m for p in sim.particles]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see above that two particles merged into one with a combined mass of 0.006.\n", "\n", "Let's now try to implement this collision function ourselves!\n", "\n", "To do this, we need to write a function which we can pass to `sim.collision_resolve`. In this case let's define `my_merge`. \n", "\n", "Now, whenever a collision occurs, REBOUND will pass our function two parameters:\n", "\n", " - `sim_pointer`: a pointer to the simulation object which the collision occurred in.\n", " - Because it is a ctypes pointer, you will need to use the `.contents` attribute to access the simulation object\n", " - `collision`: this structure contains the attributes .p1 and .p2 which are the indices of the two particles involved in the collision\n", "\n", "Using these inputs, we can define the necessary logic to handle the collision. The return value of our function determines how REBOUND proceeds afterwards:\n", "\n", " - 0: Simulation continues without changes\n", " - 1: remove p1 from simulation\n", " - 2: remove p2 from simulation\n", "\n", "Let us look at how this information can be used to implement the logic of the `merge` method for colliding particles in a totally inelastic collision." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def my_merge(sim_pointer, collided_particles_index):\n", "\n", " sim = sim_pointer.contents # retreive the standard simulation object\n", " ps = sim.particles # easy access to list of particles\n", "\n", " i = collided_particles_index.p1 # Note that p1 < p2 is not guaranteed. \n", " j = collided_particles_index.p2 \n", "\n", " # This part is exciting! We can execute additional code during collisions now!\n", " op = rebound.OrbitPlot(sim, xlim = (-1.3, 1.3), ylim = (-1.3, 1.3), color=['blue', 'green'])\n", " op.ax.set_title(\"Merging particle {} into {}\".format(j, i))\n", " op.ax.text(ps[1].x, ps[1].y, \"1\"); \n", " op.ax.text(ps[2].x, ps[2].y, \"2\")\n", " # So we plot the scenario exactly at the timestep that the collision function is triggered\n", "\n", " # Merging Logic \n", " total_mass = ps[i].m + ps[j].m\n", " merged_planet = (ps[i] * ps[i].m + ps[j] * ps[j].m)/total_mass # conservation of momentum\n", "\n", " # merged radius assuming a uniform density\n", " merged_radius = (ps[i].r**3 + ps[j].r**3)**(1/3)\n", "\n", " ps[i] = merged_planet # update p1's state vector (mass and radius will need corrections)\n", " ps[i].m = total_mass # update to total mass\n", " ps[i].r = merged_radius # update to joined radius\n", "\n", " return 2 # remove particle with index j" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can set our new collision resolution function in the simulation object." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAHWCAYAAAAsM2MeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABheUlEQVR4nO3dd3iT5f4G8DtJm6S7pbultmUVoUAZUkFZUkERlMMRQVAQkKHoEcXF+XncRzwe11E4yhBwgCiCoKjIkKFsyypQEEqhQGnpTLqb8fz+eE8jhbakI3kz7s915SpN3zf5NtbceZ73GQohhAARERE1mVLuAoiIiJwdw5SIiKiZGKZERETNxDAlIiJqJoYpERFRMzFMiYiImolhSkRE1EwMUyIiomZimBIRETUTw5SoCR566CHExcXJXYasBg4ciIEDBzb6PIVCgZdffrnF62kqR6uHnBPDlBzKsmXLoFAooFAo8Ntvv13zcyEEYmJioFAoMHz4cBkqdC/Hjx/Hyy+/jLNnz8pdSp3Ky8sxf/58DBkyBJGRkfDz80P37t3x0UcfwWQy2fz533jjDaxdu7bFH/err77CAw88gPbt20OhUDTpQwvZF8OUHJJWq8WKFSuuuX/79u24cOECNBqNDFX9adGiRTh58qSsNdjD8ePH8corr9QZphs3bsTGjRvtX9QVzpw5g8cffxxCCDz11FN4++23ER8fj0cffRSTJ0+26jEqKirwwgsvNOn5bRWmH330EdatW4eYmBgEBQW1+ONTy/OQuwCiugwbNgyrVq3CBx98AA+PP/9MV6xYgZ49eyI/P7/FnksIgcrKSnh5eVl9jqenZ4s9vyOqrKyEWq1u8Jjr/dweIiIikJaWhs6dO1vumz59OiZPnoylS5fiH//4B9q1a9fgY2i1WluX2Wiff/45oqOjoVQqkZiYKHc5ZAW2TMkh3X///SgoKMCmTZss91VXV+Obb77BuHHj6jzHbDbj/fffR+fOnaHVahEeHo7p06ejqKio1nFxcXEYPnw4fv75Z/Tq1QteXl5YsGABAODcuXO4++674ePjg7CwMDz55JP4+eefoVAosG3bNstjXH3N9OzZs1AoFHj77bexcOFCtG3bFhqNBjfddBP2799/Ta2rVq1Cp06doNVqkZiYiG+//dbq67A19W/cuBFJSUnQarXo1KkT1qxZU+u4wsJCPP300+jSpQt8fX3h7++PO++8E4cPH6513LZt26BQKLBy5Uq88MILiI6Ohre3Nz744AOMHj0aADBo0CBL93vN61DXNdPKykq8/PLL6NChA7RaLSIjIzFq1ChkZGQ0+DtdvHgRkydPRnh4ODQaDTp37owlS5Zc97UICQmpFaQ1/vKXvwAA0tPTr/sYV18zffnll6FQKHD69Gk89NBDCAwMREBAACZNmoTy8vJa55WVleHTTz+1vDYPPfSQ5ecHDx7EnXfeCX9/f/j6+mLw4MHYs2fPdesBgJiYGCiVfHt2JmyZkkOKi4tDnz598OWXX+LOO+8EAPz000/Q6XQYO3YsPvjgg2vOmT59OpYtW4ZJkybhb3/7GzIzMzFv3jwcPHgQO3furNWaPHnyJO6//35Mnz4dU6dORUJCAsrKynDbbbfh0qVLeOKJJxAREYEVK1Zg69atVte9YsUKlJSUYPr06VAoFHjrrbcwatQonDlzxvL8P/zwA8aMGYMuXbpg7ty5KCoqwpQpUxAdHW3185w6dQpjxozBjBkzMHHiRCxduhSjR4/Ghg0bcPvttwOQukDXrl2L0aNHIz4+Hrm5uViwYAEGDBiA48ePIyoqqtZjvvbaa1Cr1Xj66adRVVWFIUOG4G9/+xs++OAD/P3vf8eNN94IAJavVzOZTBg+fDi2bNmCsWPH4oknnkBJSQk2bdqEo0ePom3btnWel5ubi5tvvhkKhQKPPfYYQkND8dNPP2HKlCnQ6/WYNWuW1a9LjZycHABS2DbVfffdh/j4eMydOxcHDhzA4sWLERYWhn/9618ApNbjww8/jN69e2PatGkAYPkdjx07hn79+sHf3x/PPvssPD09sWDBAgwcOBDbt29HcnJyk+siByWIHMjSpUsFALF//34xb9484efnJ8rLy4UQQowePVoMGjRICCFEbGysuOuuuyzn/frrrwKAWL58ea3H27BhwzX3x8bGCgBiw4YNtY595513BACxdu1ay30VFRWiY8eOAoDYunWr5f6JEyeK2NhYy/eZmZkCgAgODhaFhYWW+9etWycAiO+//95yX5cuXUTr1q1FSUmJ5b5t27YJALUesz419a9evdpyn06nE5GRkaJ79+6W+yorK4XJZKp1bmZmptBoNOLVV1+13Ld161YBQLRp08byWtdYtWrVNb97jQEDBogBAwZYvl+yZIkAIN59991rjjWbzZZ/AxAvvfSS5fspU6aIyMhIkZ+fX+ucsWPHioCAgGtqup6qqirRqVMnER8fLwwGw3WPv7qel156SQAQkydPrnXcX/7yFxEcHFzrPh8fHzFx4sRrHnPkyJFCrVaLjIwMy33Z2dnCz89P9O/fv1G/T+fOnWu9zuSY2I9ADuu+++5DRUUF1q9fj5KSEqxfv77eLt5Vq1YhICAAt99+O/Lz8y23nj17wtfX95rWZXx8PIYOHVrrvg0bNiA6Ohp333235T6tVoupU6daXfOYMWNqDRjp168fAKmVCADZ2dlIS0vDhAkT4OvrazluwIAB6NKli9XPExUVZenKBAB/f39MmDABBw8etLTKNBqNpavQZDKhoKAAvr6+SEhIwIEDB655zIkTJzbquvHVVq9ejZCQEDz++OPX/EyhUNR5jhACq1evxogRIyCEqPXfbujQodDpdHXW2pDHHnsMx48fx7x582pdb2+sGTNm1Pq+X79+KCgogF6vb/A8k8mEjRs3YuTIkWjTpo3l/sjISIwbNw6//fbbdR+DnA+7eclhhYaGIiUlBStWrEB5eTlMJhPuvffeOo89deoUdDodwsLC6vz55cuXa30fHx9/zTHnzp1D27Ztr3njv94AlivdcMMNtb6vCdaa67bnzp2r9zHbtWtndXC0a9fumjo7dOgAQLp+GxERAbPZjP/85z/473//i8zMzFpTRYKDg695zLpek8bIyMhAQkJCowIsLy8PxcXFWLhwIRYuXFjnMVf/t2vIv//9byxatAivvfYahg0bZvV5dWnov6W/v3+95+Xl5aG8vBwJCQnX/OzGG2+E2WzG+fPn67zWS86LYUoObdy4cZg6dSpycnJw5513IjAwsM7jzGYzwsLCsHz58jp/HhoaWuv75rTAGqJSqeq8Xwhhk+dryBtvvIF//OMfmDx5Ml577TW0atUKSqUSs2bNgtlsvuZ4W70mDamp44EHHsDEiRPrPKZr165WPdayZcvw3HPPYcaMGU2e6nIlR/pvSY6PYUoO7S9/+QumT5+OPXv24Kuvvqr3uLZt22Lz5s245ZZbmhwKsbGxOH78OIQQtVp9p0+fbtLj1fcc9T1mY57n9OnT19T5xx9/AIBlRPA333yDQYMG4ZNPPql1bnFxsdUDc+rrnq1L27ZtsXfvXhgMBqunDoWGhsLPzw8mkwkpKSlWP9fV1q1bh4cffhijRo3C/Pnzm/w4jVXX6xMaGgpvb+865yGfOHECSqUSMTEx9iiP7IjXTMmh+fr64qOPPsLLL7+MESNG1HvcfffdB5PJhNdee+2anxmNRhQXF1/3uYYOHYqLFy/iu+++s9xXWVmJRYsWNan2ukRFRSExMRGfffYZSktLLfdv374daWlpVj9OdnY2vv32W8v3er0en332GZKSkhAREQFAalld3YpatWoVLl68aPXz+Pj4AIBVr99f//pX5OfnY968edf8rL7WnEqlwl//+lesXr0aR48evebneXl5133eHTt2YOzYsejfvz+WL19u1yklPj4+17w2KpUKQ4YMwbp162otdpGbm4sVK1bg1ltvbbCbmJwTW6bk8Orr/rvSgAEDMH36dMydOxeHDh3CkCFD4OnpiVOnTmHVqlX4z3/+U+/11hrTp0/HvHnzcP/99+OJJ55AZGQkli9fbpnU35hWWkPeeOMN3HPPPbjlllswadIkFBUVYd68eUhMTKwVsA3p0KEDpkyZgv379yM8PBxLlixBbm4uli5dajlm+PDhePXVVzFp0iT07dsXaWlpWL58ea1BMdeTlJQElUqFf/3rX9DpdNBoNLjtttvqvDY9YcIEfPbZZ3jqqaewb98+9OvXD2VlZdi8eTMeffRR3HPPPXU+x5tvvomtW7ciOTkZU6dORadOnVBYWIgDBw5g8+bNKCwsrLe+mnnBCoUC9957L1atWlXr5127drW6m7gpevbsic2bN+Pdd99FVFQU4uPjkZycjNdffx2bNm3CrbfeikcffRQeHh5YsGABqqqq8NZbb133cXfs2IEdO3YAkD5QlJWV4fXXXwcA9O/fH/3797fZ70RNJN9AYqJrXTk1piFXT42psXDhQtGzZ0/h5eUl/Pz8RJcuXcSzzz4rsrOzr3uuEEKcOXNG3HXXXcLLy0uEhoaK2bNni9WrVwsAYs+ePZbj6psa8+9///uax8RVUy+EEGLlypWiY8eOQqPRiMTERPHdd9+Jv/71r6Jjx44N/t5X1v/zzz+Lrl27Co1GIzp27ChWrVpV67jKykoxe/ZsERkZKby8vMQtt9widu/efc2UlpqpMVefX2PRokWiTZs2QqVS1Zomc/XjCCFEeXm5+L//+z8RHx8vPD09RUREhLj33ntrTRGp6/XIzc0VM2fOFDExMZbzBg8eLBYuXNjga1FTe323q5+nLlcfVzM1Ji8vr9ZxNX+bmZmZlvtOnDgh+vfvL7y8vASAWtNkDhw4IIYOHSp8fX2Ft7e3GDRokNi1a9d167myhqb+TmR/CiF4NZ2oIe+//z6efPJJXLhwoVELKzRWUlISQkNDa636VJe4uDgkJiZi/fr1NquFiBqH10yJrlBRUVHr+8rKSixYsADt27dvsSA1GAwwGo217tu2bRsOHz7M3UGInBSvmRJdYdSoUbjhhhuQlJQEnU6HL774AidOnKh3yk1TXLx4ESkpKXjggQcQFRWFEydO4OOPP0ZERMQ1CwUQkXNgmBJdYejQoVi8eDGWL18Ok8mETp06YeXKlRgzZkyLPUdQUBB69uyJxYsXIy8vDz4+Prjrrrvw5ptv1rmYAhE5Pl4zJSIiaiZeMyUiImomhikREVEz8ZrpdZjNZmRnZ8PPz6/FJu0TEZFzEEKgpKQEUVFRDa6uxTC9juzsbK6jSUTk5s6fP4/WrVvX+3OG6XX4+fkBkF5IrqdJRORe9Ho9YmJiLFlQH4bpddR07fr7+zNMiYjc1PUu83EAEhERUTMxTImIiJqJYUrkBubOnYubbroJfn5+CAsLw8iRI+vcvJqImoZhSuQGtm/fjpkzZ2LPnj3YtGkTDAYDhgwZgrKyMrlLI3IJXE7wOvR6PQICAqDT6TgAiVxGXl4ewsLCsH37dm40TdQAazOALVMiN6TT6QAArVq1krkSItfAMCVyM2azGbNmzcItt9yCxMREucshcgkMUyIXtWnTJvTp0wexsbHo06cPNm3aBACYOXMmjh49ipUrV8pcIZHr4KINRC5o8uTJWLp0qeX7rKwsDBkyBDfeeCNKSkqwY8eOBpdGI6LGYZgSuZhNmzbVCtIrpaen44MFH6DUuxSbz2yGrlKHnNIcFFcWQ1elQ3FlMfRVeigVShRXFsNoNsIkTDCajTCajTCYDJZ/G81GGMzS94mhiThRcAJeHl7w8vSC1kMLrYcWXh7Sv6P9oiEgEOwVjGDvYAR7BSPEO8Ty72DvYGg9tHZ+pYhaDkfzXgdH85Kz6dOnD/bs2VP/AZGAYry0NFrPuJ5I16UjUBsIf40/ArWBCNAGIMw7DAazAR5KD3goPeCp9LT8W6VUWf5d8zM/jR90lTpUGCtQaaxEpbESFcYKVBgqUGWqQoAmACfyT6CgogAF5QUwmA3XlOXj6YM+rfugwliBuMA4xAbESl8DYxEbEIsovyiolCpbvWxEdbI2A9gyJXIx2dnZDR9wCRBvS5+hH170MKb/fbodqvqTEAJlhjLkl+ejoLzAErD55fmoMlYhPT8dGUUZ+CXzF+SW5VrO81R6orV/a0vQJkUkISYgBp1COyHaL5pbJJKsGKZETqCmA0ngz44kpaLu8YNRUVHIysqq97E6dOgg6+pHCoUCvmpf+Kp9ERcY1+CxFYYKZOmycLb4LM7pzuFc8Tmc1Z3FwZyDOHL5CA7lHAIA+Kn90Cm0k+XWJbwLuoR1ga/a1/a/EBEYpkQOp67grKGA1PpqqBX26quvYsiQIfX+fPz48c2s0H68PL2QEJKAhJCEa34mhMAF/QUczzuO43nHcSzvGPZd3IcvjnyBbuHdcCj3EDqHdkbv6N6WW2xALFuwZBO8ZnodvGZKtnTl/35Xh6c1wVmfhx9+GJ988kmdP/v2228xcuTIRj+mszCYDDhdeBq/Z/+OfRf3Ye/FvThVeAoAEOYTJgVrVG8kt05Gt/Bu0HhoZK6YHJm1GcAwvQ6GKbUUIQRMwgQllMBV+ai44o6Wajlt2bIFL7zwArKzsxEVFYXXX38dKSkpLh+mdSmsKMTv2b9j74W92Je9D6nZqagwVkCtUqNbeDckRyfjlphbkNw6GUFeQXKXSw6EYdpCGKbUHGZhhslsgkmYYBZmALCMjAVaLjitpVAo3DJMr2YwGSzdwjUBG+0Xjd+zf0f3yO4YFDcIg+IGoVdUL3iqPOUul2TEMG0hDFNqrJrwNJlNEBBQQAGVUgWVQgWlQmmXAN20aRNefPFFS6v01Vdfxe23384wbcB53XnsOLcDW89uxbaz21BUWYT+sf0R7hOOER1GYHCbwZwL64Y4NYbITmq6b2tCFJBG2tbMyaxv1K2t1Lf60ZQpU+xah7OJCYjB+K7jMb7reJjMJqRdTsOOczvwzfFvsDp9Nbw9vXF7m9txd8LdSGmTAm9Pb7lLJgfClul1sGVKdamr+1alUFlaoHKNGN20aVODI3kB1x+AZAsZhRn4/o/v8f0f3+NI7hFoPbRIaZOCER1G4PY2t8NP4yd3iWQj7OZtIQxTquEI3bfXc93VjwA888y7uOuuQYiLa4XY2BvsVJnrOFd8zhKsBy4dgFqlxm3xt2FEhxEY2nYoArQBcpdILYhh2kIYpu7NLMyW9WgFhKX16aH0sEv3rRB/3pRK6WsNZR1PHxsb2+CCDVfy9JyIkSOXQa0GOncGEhOlr3FxdT82Xeu87jx+OPUDvjv5HfZn74en0hMD4gbg7g534452d3BksAtgmLYQhql7MpgMqDZVwyRMUEABtUoNT5WnTQO0JjSv/jcAKBR/3mq+r8v1WqZJSTfjo4924+JF4Px5ICcH+P134OhRoKhIOsbbG7j3XqB1ayAlBejdG/DkgNbrulRyyRKsey7sQfeI7oj2j8b9ifdjUNwgeKg4RMUZMUxbCMPUfZiFGdWmahhMBggIyyLutpoaIQRgNv8ZnFcGZV3/tsb1rplu3rwZgwcPrrOWnBwpVI8eBf74A1i3TgpYPz9g4EBg8GDg9tuBNm0a8Uu6qctll7H2xFqsPLoSOaU58FP7YUK3CXig6wPsBnYyDNMWwjB1fUazEdWmahjNRiiggKfKE2qV2iatUCEAk6l2y1OplAKz5mtz1bf60ZQpU7B48WKrH8dkAg4dAjZvlm579gBGIxAfL7VYR4wA+vUDtJwt0qC03DQsObgEq9NXw0Ppgfs634eHezyMdq3ayV0aWYFh2kIYpq6p5lpotanaci20pivXFkwm6WaWBv5CpZJutro2WdfqR3W1SBujpATYvh3YskUKVy8vqTU7ejQwYQLQrVsLFe+i8sry8PmRz7Hs0DLkledhUNwgTOs5DQNiBzjE4DWqG8O0hTBMXUtdrVBPpadN9skUQmrJGY1/tjxrQtQVnDwJ/PADsGABkJ0N9OoFPPwwMGqUFLRUt2pTNb47+R0Wpi5E2uU0tGvVDlN7TMW9ne7l3FUHxDBtIQxT5yeEQLWpulYrtCZEbdEiEAIwGKQQBQAPD2kAj6s2PoxGYMMGYPFiqdUaGAg88AAweTLQvr3c1TkuIQT2Z+/HwtSF+On0T/BT++GBrg9gUtIkRPtHy10e/Q/DtIUwTJ2XEAJVpipUGasgIKBWqaFRaWzSCpWe788QVSikEPXwcN0QrcuZM8CSJcDnnwOFhcBDD0mh2r273JU5tvO681hycAlWHF2Bsuoy3NX+Ljzc42H0iurFLmCZMUxbCMPU+ZiFGVXGKlSbqgFAClEPjc2mtZjNtbtz3TFEr1ZZCXz7LfDuu9LI4LvvBv7xD6BDB7krc2xl1WX4+tjXWHxwMc4UnUFSRBKm9piKER1GcMF9mTBMWwjD1HkIIVBhrECVsQqAtLG0RqWx2Sd7IYDqaqk1qlJJXbkenEpYi8kErFwJvPGGdF113Djg+eeBmBi5K3NsZmHGL5m/YPGBxdh+bjtS4lMwrP0wjO482rLjENkHw7SFMEwdnxAClcZKVBorAQBaDy20Hlqbdo8ZDFKQAoBazUUNrqeqCli2DPj3vwGdTur6ffppIDRU7soc38n8k/jy6JdYcnAJ4oPi8fdb/44hbYew+9dOGKYthGHquGquiVYYKiAgLCFqy1WKzGagokJqcanVgEbj3t25jVVWBnz0EfCf/0jd4nPmAFOncvSvNY5dPoZ//vpP/Jr1K26Kugn/6P8PdI/kxWhbY5i2EIapY6oyVqHCWAGzMEOj0sDL08vma+VWVUk3pVJ683eVKS5yKCoC3n8fSEsDsrKADz8E+vSRuyrnsOPcDry+43Wk56djRIcReO6W5xAbGCt3WS6LYdpCGKaOpdpUjXJDOczCDLVKDS8PL5uNzq1xZWtUo5Fu1DLOnAEeewzYtw+YPh144QW2Uq1hMpuwJn0N3tr1FvLL8zGx20Q8kfwEF9a3AYZpC2GYOgaDyYByQzlMwmS3EAWk66IVFVIrlK1R2zCZpIUfXn8diI6WWqk33yx3Vc6hwlCBxQcW47+//xcKKPB478cxuftkaDz4ia+lWJsB3GiJHJrJbIK+Sg9dlQ4A4K/xh6/a1y5BWlEBlJdLg4t8fBiktqJSAY8+CuzYAQQHA8OHA//3f9LrTw3z8vTC48mP49dJv2LUjaPw1q63MGDZAKxJX2PZtJ7sw6nCdMeOHRgxYgSioqKgUCiwdu3a656zbds29OjRAxqNBu3atcOyZctsXic1nxAC5YZyFFcWw2g2wk/thwBtgF2mBQgBlJZK10e9vKQtyTjIyPbatZOWJ3z5ZWDpUqB/f2DvXrmrcg4h3iF4/bbXsWXCFnQN74onNjyB4SuGY2fWTrlLcxtOFaZlZWXo1q0b5s+fb9XxmZmZuOuuuzBo0CAcOnQIs2bNwsMPP4yff/7ZxpVSc1QZq1BUWYQKQwW8PL0QpA2yW7eVySQt6G4yAb6+vD5qbyqVdA11+3agVStpXup770n/Pej62gS1wcIRC7HmvjVQq9QYu3osJnw7ASfzT8pdmstz2mumCoUC3377LUaOHFnvMc899xx++OEHHD161HLf2LFjUVxcjA0bNlj1PLxmaj8mswml1aUwmA1Qq9Tw8fSxS3dujZrro0ql1K1rqx1dyDo111JfflnaS/Xjj4EAbgVqNSEEfjz1I+b+Nhfn9efxSM9HMLXnVAR7B8tdmlPhNVMAu3fvRkpKSq37hg4dit27d9d7TlVVFfR6fa0b2V5ZdRkKKwphEiYEaALgr/G3a5BWVv55fdTXl0HqCGqupa5cCezfL21MfuKE3FU5D4VCgbs63IWtE7fipQEvYUPGBsz8cSbWnVgHJ21DOTSXfsvIyclBeHh4rfvCw8Oh1+tRUc/ohrlz5yIgIMByi+G6ZzZlNBtRWFGIMkOZpUvX3muQlpZKLVJeH3VMt90m7Z+q1QJDhwLr18tdkXPxVHlicvfJWH3fakT7RWPWz7Mw88eZKKwolLs0l+LSYdoUc+bMgU6ns9zOnz8vd0kuSQiB0upSy//QrbxawVfta/cl0vR6aaCRjw+vjzqyuDjgp5+k1umkSdL0GTMHqzZKsHcw3hn6DuYPm489F/ZgyOdDsDFjo9xluQyXDtOIiAjk5ubWui83Nxf+/v7wqmdmuEajgb+/f60btaxqUzUKKwpRYaiAr9oXrbxaybJ4d0mJdJ3Uz09aGpAcm48PsGiRdA31/felLuCa9ZHJesPaD8PPD/yM7hHdMX39dMz+eTb0Vbyc1VwuHaZ9+vTBli1bat23adMm9OG6ZbIQQqCkqgTFlcVQKpRo5dUK3p7estRSUiK1SP39GaTORKEAZs6UwvSHH4CJE6XeBWqcUJ9QLByxEG/f/jY2ntmIoV8Mxa/nfpW7LKfmVGFaWlqKQ4cO4dChQwCkqS+HDh1CVlYWAKmLdsKECZbjZ8yYgTNnzuDZZ5/FiRMn8N///hdff/01nnzySTnKd2sGkwH55fmoMFTAT+2HIK8guw4wulLNHFK2SJ3XiBHAN99IH4rGjWOgNoVCocBfO/0VGx/YiLZBbTFh7QS88MsLKDeUy12acxJOZOvWrQLANbeJEycKIYSYOHGiGDBgwDXnJCUlCbVaLdq0aSOWLl3aqOfU6XQCgNDpdC3zS7ihkqoScankksgvyxdGk1HeWkqEyMsTorJS1jKohRw+LESHDkIMGyYE/xdtOrPZLD4//Lm4cd6Nov/S/mLvhb1yl+QwrM0Ap51nai+cZ9p0JrMJuiodqk3V8FX7wlftK2s9ZWXSqF1fX2lkKLmGI0eA0aOlFZS+/FLquqemOVd8Dk9vehqp2amY0n0KZvedDa2He//PwnmmJKtKYyUKKgpgMpssI3XlVFEhhamPD4PU1XTtCqxaBWRkAPffzy7f5ogNjMXKv67EnFvn4LMjn2H4iuE4kntE7rKcAsOUWpQQAvoqPYori6FWqRHsHQy1St4LkwaDdJ3Uy4vbe7mqKwN17FgGanOolCpM7TkV6+9fD29Pb4z6ahTe2fUODCaD3KU5NIYptRiDyYCCigJUGCoQoAlAoDbQ5ht2X4/ZDOh0gIeH1L1LrqtLFylQz5xhoLaE9sHtsWbMGjyR/AQWpi7EYz8+hixdltxlOSyGKbWIsuoyFJQXQAEFgr2D4eUpfxNQCClIAWlNV65s5PquDNQxYxiozeWh9MDjyY9j9ZjVOFV4CiO+HIGfT3OjkLowTKlZhBAoriyGrkoHradWtgUY6lJaChiNUpByrV33UROoGg3w8MPSNChqnsSwRKwbuw79buiHBakL8M6ud7hf6lX4FkNNZjKbLHNHA7WBCNQG2n05wPpUVEg3Pz9p8XpyL126AC+8APz+O/DEE1x6sCX4afzw4Z0f4q72d2FB6gLM/GEm56RegWFKTVJlrEJeeR7MwowQ7xDZVjKqi8EgTeb38uLIXXfWqxcwf760UpKVWyDTdSgUCkzpMQULRyzE7gu7MeabMbiovyh3WQ6BYUqNVlpVioKKAngqPRHqE2r3XV4aIgRQVCRt3+XnJ3c1JLc77wTeeAN46y1g9Wq5q3EdA+MGYtXoVSitLsWor0chNTtV7pJkxzClRimqKEJxVTF81b4I9g6WfbTu1fR6aVPpwEC5KyFH8cAD0qIOzz4rLfBALaN9cHusuW8N2ga1xQPfPoA16WvkLklWjvVOSA5LCIH88nyUGcoQpA2Cv8bxlpmprpYWZvD3l1qmRIA0invuXKBTJ2lAUn6+3BW5jiCvIHw68lOM6jgKz21+Dm/+9iZMZpPcZcmCYUrXZTKbcLnsMqqMVQjxDoGP2kfukq4hBFBcLC1c7+N45ZHMNBpp+zajEZg2TbquTi3DU+WJ1297HS/0ewFLDy3FjPUzUFpdKndZdscwpQYZTAZcLrsMszAjzCfMYdfpLC1l9y41LCICWLgQOHgQePFFuatxLQqFAhOTJmLxiMVIvZSK0atGu90CDwxTqlfNiF2lQokwnzCHGmh0pZrRu76+0kpHRPXp1Qv45z+BQ4eA77+XuxrX0y+2H7657xsYTAaM+moU9l7YK3dJdsMwpTqVG8qRX55vGbEr196j1uBygdQY998PtGkD/P3vwKVLclfjetoEtcHq+1ajU2gnTFw7ESuPrpS7JLtgmNI1SqpKUFhRCC9PL4R4hzjciN0rlZZKA48CA7lcIFlHoZBap97ewOzZXNDBFgK0AVhyzxKM6TwGL297GfP3zXf5FZMc912SZFFUUQRdlQ7+Gn+08mrlMCsa1cVkkqbC+PhIA4+IrBUYCLzzDrBzJ7B0qdzVuCYPpQdeGfQKXh74Mj7c9yGe3/w8jGaj3GXZDMOUAPxv6kuZY099uVrNIuZcnIGa4tZbgcmTgTffBP74Q+5qXNfYxLF4d+i7+OHUD5j5w0xUGCrkLskmGKYEszDjctllVBgrEOLlmFNfrmY0/jmnlIvYU1M99xwQGyut38vpMrYzrP0wLBi+AHsv7sVzm59zyakzfBtyczVBajAbpKkvno459eVqOp20MAPnlFJzaLXABx9ILdN33pG7Gtd26w234tORn+JwzmFM/X4qyqrL5C6pRTFM3VhNkBrNRoT5hEHjoZG7JKtUV0s7wvj7c9ARNV+nTsDTTwOffw4cOCB3Na6tW0Q3zBs2DxmFGZi2fppLBSrD1E2ZhRm5pbkwmo0I9wmHWuU8I3j0emlbNbZKqaVMmwb07Am8+qr0YY1sp0t4Fyy+ezH+KPgD09dPd5lAZZi6oauD1FEXY6hLVRVQWSm1SolaikoF/N//AWfPAp99Jnc1rq9reFd8cvcnOJF/AtPXT3eJfVEZpm5GCGG5Rhru61xBCvy5QIOXl9yVkKtJSABGjgQ+/BAoLJS7GtfXNbwrFt+92GUClWHqRmqCtMpY5XRdu4DU/cZWKdnSE09I1+Hfe0/uStxDUkQSFo1YhON5xzFj/QynnjbDMHUj+eX5qDRWOtVgoyvp9VJ3nLe33JWQqwoKAv72N2DlSuDkSbmrcQ/dI7tj0YhFOJZ3DDN+mIFKY6XcJTUJw9RN1OxFGuoTCi9P5+sjNZmkeaVcoIFs7cEHgRtukJYcFELuatxDj8geWDh8IdJy0zBjvXMGKsPUDRRWFKK0uhQh3iHw9nTOZl3p/+Z4czF7sjVPT2kR/J07ga1b5a7GffSM6olFIxbhSO4RPPLDI04XqAxTF1dcWQx9lR7BXsHwVTtvEpWWSlNhuNoR2cNttwG33AK88QZXRrKnnlE9sWD4AhzKOYRHf3jUqQKVb00uTFepQ3FlMYK0QfDTOG//aGWl9IbGVinZi0IhTZU5dw744gu5q3EvN0XfhAXDF6CgvABzNs+ByWySuySrMExdVElVCYoqixCgCUCANkDucpqlpETqetM6x0qH5CISEoAxY6TlBouK5K7GvfSO7o2n+z6NLZlb8Mr2VyCc4OI1w9QFlVWXoaCiAP4afwR5BcldTrOYTEB5OQcekTyefBLo0AH45BO5K3E//WL74bVBr2F1+mrM2zdP7nKui2HqYsoN5cgrz4Ov2hetvFrJXU6zlf1vpTEuHUhyCA4GkpOlVZG4kIP93dPxHjx181P4OPVjrDy6Uu5yGsQwdSFVxirkluZCq9Ii2CtY7nJaREmJNK9UpZK7EnJXkydL11DZOpXH5O6T8WDXB/Gvnf/Cr+d+lbucejFMXYTRbERuWS7UKjXCfcOhcIHtVKqqpIFH7OIlOQUGAlOmAKmpf25IT/ajUCjw7C3PYkSHEZi9cTZOFZySu6Q6MUxdgBACOaU5MJvNiPCNcIkgBaRWqVLJgUckv3HjgKNHga+/lrsS96RUKDHn1jm4IeAGPPbTYyiqcLwRYQxTF5Bfno8qYxUi/CKgUrpOf2hZGa+VkmMICQFGjJD2PDUa5a7GPXl5euGDOz9AhaECszfOhtHsWP8hGKZOrqSqBLoqHUK8Q6D1cJ0mXGWlNJKXc0vJUTz0EJCTA2zYIHcl7ivKLwrv3/E+Dlw6gLd2viV3ObUwTJ1Ytakal8suw0/t5/RzSa9WViYNOtI433r85KISEoC+fYGlS7lmr5x6RPbAC/1fwJdHv8Q3x7+RuxwLhqmTMgszLpVcglqlRphPmNzltLiyMu5ZSo5n0iTp2mlqqtyVuLd7O92LsZ3H4vUdr+PApQNylwOAYeq0cktzYRImlxpwVMNolPYu5VZr5GhuvRVo2xZYvlzuSui5W59Dj8gemLVhFrJLsuUuh2HqjIoqilBmKEO4Tzg8VZ5yl9PiysulrwxTcjRKpdQ63b4duHxZ7mrcm4fSA+8MeQdenl74209/k31jcYapkyk3lKOgogCtvFrBR+2aQ13Ly6XpMNwhhhzRHXcAZjPw3XdyV0JBXkGYd+c8ZOmy8MIvL8i6hi/frpyI0WxEbmkuvD29XWKpwLoIIYUpW6XkqPz8gJQUYM0aDkRyBO2D22Pu4LnYeGYjFqYulK0OhqmTEELgUsklKBQKhPuEy12OzVRVSZ/6GabkyEaNAjIzgSNH5K6EAGBwm8F4tNej2Hp2K34795ssNTBMnUReeR6qTdWI9I10qYUZrlZRIX3a55QYcmTJyUBEhNQ6Jccwrec0hHiH4P+2/h/yyvLs/vwMUyegr9JDX6VHqE8oNB6unTKVldL1UhcboEwuRqUC/vIX4Mcfpb9Zkp9KqcJrg16Dp9ITc7bMgVmY7fr8DFMHV2WsQl5ZHvw1/vDX+Mtdjs3VhCmRoxs5EigtBTZvlrsSqhHkFYQ3U95E6qVULEpdZNfnZpg6MJPZZFmYIdQ7VO5ybM5kknaJYZiSM7jhBmDIEKl1So6jV1QvTOsxDR+nfozUbPutrsEwdWC5pbkwmo2I9It0uYUZ6lLTXcbrpeQsbroJ2LlT2uGIHMf0XtPRI6IHnt/yPHSVOrs8J8PUQZVWl0JfLV0n9VB6yF2OXVRVSXNL1Wq5KyGyzu23Syt2/fKL3JXQlZQKJeamzEVrv9Z449c37POcdnkWahSzMCOnNAc+nj4ut4B9Q3i9lJxNeLi0iMPvv8tdCV0tzCcMYxPHYuOZjdh8xvYXthmmDiivLA9mIW307U4YpuSM2rUDtmyRrvmTYxnSdggGxw/GG7++YfMNxRmmDqbCUIGiyiKEeoe65Lq79TEapRuvl5KzueUW6Zrp0aNyV0JXUygUmHOrNE1m7m9zbfpcDFMHIoTApdJL8PLwQpBXkNzl2FXNykcMU3I2nTtLSwzu2iV3JVSXYO9gzLl1Djad2YRNGZts9jwMUweSX54Pg8ngdt27gBSmAAcfkfNRqYCbb2aYOrKa7t65v821WXcvw9RBVBmrUFhRiGDvYJdf5aguRiPg4R6DlskF9e0rrdPLKTKOyR7dvQxTB3GpVFqcIdgrWO5SZGEwMEzJefXtK12m2LtX7kqoPrbu7mWYOoDCikJUGisR4RvhFosz1MVgADzdZ7wVuZioKCA2ll29js4yuve3N1BYUdiij80wlZnBZEB+eT6CtEHw8vSSuxzZsJuXnF3fvtJqSNzj1HHVdPcKIfDmb2+26GMzTGWWU5oDlUKFUB/XX3u3IWyZkrPr2xfIzgbOn5e7EmrIld29O87uaLHHZZjKSFepQ5mhDOG+4VAq3Pc/hckkfZpnmJIz690bSEoC0tLkroSuZ0jbIbinwz1447c3UG4ob5HHdN93cJmZzCZcLrsMf40/fNW+cpcjK4NB+spuXnJm3t7SFK+DB+WuhK5HoVBgWq9p0FXpsDB1YYs8JsNUJrlluQCAcJ9wmSuRX02YsmVKzq5DB+DkSbmrIGtE+UVhSvcpWJG2AmeKzjT78ZwuTOfPn4+4uDhotVokJydj37599R67bNkyKBSKWjetAyz+WlpdCn2VHuG+4VApVXKXIzujEVAo2DIl55eQAJw6JU2TIcc3odsERPpF4q2db0E0c+SYU4XpV199haeeegovvfQSDhw4gG7dumHo0KG4fPlyvef4+/vj0qVLltu5c+fsWPG1hBDIKc2Bt6c3/DX+stbiKDjHlFxFQoLU1ZuVJXclZA21So1n+z6L/dn7m72zjFOF6bvvvoupU6di0qRJ6NSpEz7++GN4e3tjyZIl9Z6jUCgQERFhuYWHy9utWlBe4LZLBtanpmVK5OwSEqSv7Op1HrfccAsGxA7Au3vebdZgJKcJ0+rqaqSmpiIlJcVyn1KpREpKCnbv3l3veaWlpYiNjUVMTAzuueceHDt2rMHnqaqqgl6vr3VrKSazCXnleQjUBkKt4iK0NUwmaX1TImcXECDtccowdS6z+8xGuaEcq46tavJjOE2Y5ufnw2QyXdOyDA8PR05OTp3nJCQkYMmSJVi3bh2++OILmM1m9O3bFxcuXKj3eebOnYuAgADLLSYmpsV+h4KKApiF2e3nlF7NbAaUTvOXSNQwDkJyPtH+0RjbeSwWHViEvLK8Jj2GS7+F9enTBxMmTEBSUhIGDBiANWvWIDQ0FAsWLKj3nDlz5kCn01lu51toBrbJbEJ+eT6CvYPhoeQFwisxTMmVJCQAf/whdxXUWA92exAaDw0WHVjUpPOd5i0sJCQEKpUKubm5te7Pzc1FRIR11x89PT3RvXt3nD59ut5jNBoN/P39a91aQl659Gkn1Jut0qtx5CO5koQEID8fKCiQuxJqDF+1LyYlTcK6k+twrrjxA1WdJkzVajV69uyJLVu2WO4zm83YsmUL+vTpY9VjmEwmpKWlITIy0lZl1sloNqKwohAh3iGcClMHITgAiVxHhw7SV7ZOnc99ne9DqHcoPvr9o0af6zRhCgBPPfUUFi1ahE8//RTp6el45JFHUFZWhkmTJgEAJkyYgDlz5liOf/XVV7Fx40acOXMGBw4cwAMPPIBz587h4YcftmvdeWV5UEDhtturEbmT1q2BxETg4kW5K6HGUqvUmN5zOrZkbsHxvOONOtepLt6NGTMGeXl5ePHFF5GTk4OkpCRs2LDBMigpKysLyisuvhUVFWHq1KnIyclBUFAQevbsiV27dqFTp052q9lgMqCwohBhPmFsldaDLVNyJUqlNN3r1Cm5K6GmGNZ+GD4/8jnm7ZuH/971X6vPU4jmLvvg4vR6PQICAqDT6Zp0/fSi/iJKqkvQIbiDWy9m35BTp6R1TaOj5a6EqGW8+qq0eMM//yl3JdQU285uwzObnsG8O+fhRv8brcoAvrvbULWpGsWVxQj1DmWQNoAf58jVaDTA2bNyV0FNNSB2ALqEdcG8/fNgFtaNkOQ7vA1dLrsMD6UHWnm1krsUIrKj4GBpRC85J4VCgcd7P44T+SewLXObVecwTG2k0lgptUp9QqHgBcEGCcHpMeRagoOB4mL+XTuz7pHdcUvMLfjk4CdWHc8wtZHLZZehVqkRpA2SuxSHp1LxTYdcS0iI9DddVCR3JdQcM2+aiSpjlVXHMkxtoMJQAX2VHmE+YWyVWkGpZJiSawn+3yw4dvU6t/bB7dE5tLNVxzJMbeBy2WVoVBoEaALkLsUpMEzJ1dSEKVdBcn4PdHvAquMYpi2s3FCOkuoStkobQaWSdo4hchWt/jfmkGHq/NoHt7fqOIZpC7tcehlalRYBWrZKrcWWKbkaT09pOzaGqftgmLagSmMl9NV6BHtz2cDGYJiSK+rUCTAY5K6C7IVh2oLyy/LhqfREoDZQ7lKcikrFNx1yPaWlwOXLcldB9sIwbSFGsxEFFQVo5dWK10obqeaaKVdCIlfi7Q2Ul8tdBdmLUy1078iKKqQJZezibTyP//0VGo3StSYiVxAWxg0c3AnDtIXkl+cjUBsIDyVf0sZSq6U3nepqhim5jrIyjlJ3J+zmbQGl1aWoMlWxVdpEarX0plNl3UIjRE6DA+vcB8O0BRSUF0Cj0sBX7St3KU5JpZK6ehmm5EqUSo4DcCcM02Yymo0orixmq7SZNBqGKbkWhql7YZg2U2FFIQBwm7Vm0mqBykq5qyBqOQoFu3ndCcO0mQrKCzjwqAWwZUquRqFgy9SdMEybgQOPWo5GIy3cwNGP5CoYpu6FYdoMHHjUcrRa6Stbp+QqNBrA31/uKsheGKZNVDPwKMQ7RO5SXIJGI31lmJKr0OulhUjIPTBMm6hm4FGQV5DMlbgGDw9pigwHIZGrqKr680MiuT6GaRNx4FHL02rZMiXX4ecHREbKXQXZC8O0CTjwyDa0WmmnDSJXcOYMByC5E4ZpExSUF0DroeXAoxbm7S2tZ8q5eeQKioulDcLJPTBMG8lkNqGwohBBWl4rbWl+ftInebZOydmZTEBJCRAYKHclZC8M00bSVekgILjikQ14eUmDkEpK5K6EqHlKSqQPhgxT98EwbaTCikL4qn3hqeJeYS1NoQB8fRmm5PyKi6Wv7OZ1HwzTRjCajSipKmEXrw35+bGbl5yfTid9ZZi6D4ZpIxRXFAPg3FJb8vOTJrpXVMhdCVHT1Qw+Ypi6D4ZpI+iqdPDy9OLcUhvy/d8AaXb1kjPLy5PmTHM5QffBMLWSEAK6Kh0CNPyoaUsqFeDjwzAl55abC4SHS+MAyD0wTK1UUl0CszAjQMswtTU/P4YpOTedDrjhBrmrIHtimFpJX6mHWqWGt6e33KW4PF9f6ZqpwSB3JURN88cfQFiY3FWQPTFMrcQuXvsJCJDm6NWMiCRyJkYjcOECEB8vdyVkTwxTK1WZqtjFayeentJ106IiuSsharyLF6UVkGJj5a6E7IlhaiUFFPDXcGievQQEAAUFXCicnM/Zs9LXuDg5qyB7Y5hayU/jB6WCL5e9hIRI3WXs6iVnc+6c9GGQSwm6F6aDlXi91L58faUt2fLz5a6EqHEyM9kqdUcMUyuxi9f+goPZ1UvO59w5hqk7YphaSeOhkbsEtxMSIk2P0evlroTIOgaD9OEvIUHuSsjeGKbksPz8ALVaap0SOYOMDCAri9Ni3BHDlBxaTVcvkTM4eVJqlbZpI3clZG8MU3JoISHSguFcXpCcwYED0ib3arXclZC9MUzJofn7AxoNcPmy3JUQNcxsBo4cAZKS5K6E5MAwJYemUACtWklhylG95MgyMqSN7Rmm7olhSg4vIkIaJVlYKHclRPU7fFjq3u3YUe5KSA4MU3J4Pj7SyN6LF+WuhKh+hw4BiYnS2tLkfhim5BQiI6WWaXm53JUQXctkAtLS2MXrzhim5BTCwqRP/NnZcldCdK1Tp6QPet26yV0JyYVhSk5BqZRapzk5UiuAyJGkpUnXSjt0kLsSkgvDlJxGVJQUpDk5cldCVNuvv0oLNXh4yF0JyYVhSk5Dq5UWceBAJHIk2dlSN2+vXnJXQnJimJJTiY6Wrk0VFcldCZFkxw7pg17v3nJXQnJimJJTCQyUpsmcPy93JUSS7duB5GRppS5yXwxTcjrR0dKm4aWlcldC7u78eeDsWWDAALkrIbkxTMnphIdLi4mfPSt3JeTutm8HvL2Bnj3lroTkxjAlp6NQALGxbJ2SvISQrpf26cNdYohhSk4qIkJaxCEjQ+5KyF1lZkrdvOziJYBhSk5KqQTatgXy8gCdTu5qyB1t3ixdv+cSggQwTMmJhYcDvr7A6dNyV0LupqoK2LMHSEnhwvYkYZiS01IogHbtpDmnBQVyV0PuZOtW6Zr9wIFyV0KOgmFKTi0kRJp7evo0Nw8n+xAC+P57aZGGiAi5qyFH4XRhOn/+fMTFxUGr1SI5ORn79u1r8PhVq1ahY8eO0Gq16NKlC3788Uc7VUr20q4dUFLCNXttxWw248svv5S7DIdx5AiQlQXcfbfclZAjcaow/eqrr/DUU0/hpZdewoEDB9CtWzcMHToUly9frvP4Xbt24f7778eUKVNw8OBBjBw5EiNHjsTRo0ftXDnZUmCgtEXbqVPcUcYW3nvvPYwbNw7vvPOO3KU4hO++A+LigC5d5K6EHIlCiMZ1jk2cOBFTpkxB//79bVVTvZKTk3HTTTdh3rx5AKRPzDExMXj88cfx/PPPX3P8mDFjUFZWhvXr11vuu/nmm5GUlISPP/7YqufU6/UICAiATqeDv79/y/wi1OIqKoCdO4H4eGmUL7WcTp06IT09HZ06dcKxY8fkLkdW2dnAjBnAY48BQ4bIXQ3Zg7UZ0OgNg3Q6HVJSUhAbG4tJkyZh4sSJiI6Oblax1qiurkZqairmzJljuU+pVCIlJQW7d++u85zdu3fjqaeeqnXf0KFDsXbt2nqfp6qqClVVVZbv9Xp98wonu/DykloL585J+556e8tdkXObPHkyLl26BAA4ceIEACA9PR133nknACAyMhJLliyRrT65rF0r7VvKgUd0tUZ3865duxYXL17EI488gq+++gpxcXG488478c0338BgMNiiRgBAfn4+TCYTwsPDa90fHh6OnHouluXk5DTqeACYO3cuAgICLLeYmJjmF092ER8v7Sd59CgHIzVHaWkpli1bhg0bNmDDhg2o6bwSQljuW7ZsGUrdbPmpggJgyxZpUXuueERXa9I109DQUDz11FM4fPgw9u7di3bt2uHBBx9EVFQUnnzySZw6daql67SbOXPmQKfTWW7nuT2J01CpgM6dgcJC4MIFuatxXr6+vti3bx8CAgLq/HlAQAB+//13+Pr62rkyea1ZI221NmyY3JWQI2rWAKRLly5h06ZN2LRpE1QqFYYNG4a0tDR06tQJ7733XkvVCAAICQmBSqVCbm5urftzc3MRUc/49IiIiEYdDwAajQb+/v61buQ8goOB1q2BkyeBykq5q3FevXr1arDHp0ePHnauSF75+dIo3nvvlS4pEF2t0WFqMBiwevVqDB8+HLGxsVi1ahVmzZqF7OxsfPrpp9i8eTO+/vprvPrqqy1aqFqtRs+ePbFlyxbLfWazGVu2bEGfPn3qPKdPnz61jgeATZs21Xs8uYaEBKmVevy43JU4ty+++KLO+5cvX27nSuS3ciWg1wN33CF3JeSwRCMFBweLoKAg8eijj4qDBw/WeUxRUZGIi4tr7ENf18qVK4VGoxHLli0Tx48fF9OmTROBgYEiJydHCCHEgw8+KJ5//nnL8Tt37hQeHh7i7bffFunp6eKll14Snp6eIi0tzern1Ol0AoDQ6XQt/vuQ7eTkCPHTT0JcuiR3Jc6rX79+AoCIiYkRu3btEq1btxYARL9+/eQuza6ysoT4y1+EWL9e7kpIDtZmQKNH87733nsYPXo0tFptvccEBgYiMzOzyQFfnzFjxiAvLw8vvvgicnJykJSUhA0bNlgGGWVlZUGp/LOx3bdvX6xYsQIvvPAC/v73v6N9+/ZYu3YtEhMTW7w2cizh4dItPR1o1YoDRppi5syZuPnmm/Hmm29CqVTi3LlzeP7559HTzTbvXLtWGiE+dKjclZAja/Q8U3fDeabOq7pa2m8yKIibN1PTpKYCb7wBPP20tG8puR9rM8CpVkAiagy1GkhMlAaPnD0rdzXkbAwG4JNPgK5dgZtvlrsacnQMU3JpERHSnpPHjwPFxXJXQ85k7Vrpg9jDD0s7FBE1hGFKLq9TJ8DfX+qyq66WuxpyBrm5wOrV0mL2dljgjVwAw5RcnlIpXTM1mYCDB7k6EjVMCODzz4GAAGD0aLmrIWfBMCW34OUFdO8O5OVJu8sQ1WfXLiAjA5g6FdBo5K6GnAXDlNxGaKi0oMMff0ihSnS1vDxg8WLgxhuBXr3kroacCcOU3Eq7dtLepwcOAOXlcldDjsRsBubPl3YcmjJF7mrI2TBMya0oFFJ3r1YL/P67NP2BCJA2/T5xApg5E/DxkbsacjYMU3I7np5Ajx5ASQmwd680MInc25kzwNdfA/fcI43+Jmoshim5JT8/aSJ+UZE0ZYYjfN1XZSWwcCFwww0cvUtNxzAltxUcDNx0E3DpEnD4sNzVkByEkIJUpwOeeELaXJ6oKRim5NYiIoCkJGm5wRMn5K6G7G3TJmDPHmDiRGkxe6Km4ucwcnuxsUBVlbTkoEYDxMfLXRHZw+HDwGefAaNGce1daj6GKRGADh2kQD18WBqg1Lq13BWRLV28CHz4odQrMWqU3NWQK2CYEv1PYiJgNEoDksxmaUAKuR6dTlqYIThYmgaj5MUuagEMU6L/USiklorJBOzbJ31ll69rKSsD/vUv6es//iEtM0nUEhimRFdQKKRF8VUqqYVqMEhdwOT8qqqAt9+WpkO98AIQEiJ3ReRKGKZEV6lZJcnTEzhyRArUzp3lroqaw2gE/vMf4MIFYM4cbqtGLY9hSlSPxEQpUNPSpH1Qk5K4SbQzMpmka6Tp6cAzzwBt2shdEbkihilRAxISpEA9ehSoqACSk6UuYHIOBgPw3/9KywU+8QSXCiTb4Tg2outo00a6jnrpErB1q7T8HDm+ykrgvfek+cNTpkg9C0S2wjAlskJ0NHDbbVLrdOtWoLBQ7oqoIaWlwFtvSStbzZ4NdO0qd0Xk6himRFYKCgJSUqT1WzdvBjIz5a6I6lJYKLVI8/KA557jaGyyD4YpUSN4eQGDB0vzT/fulfZE5RZujiMzE3j9dWnRjTlzpKUiieyBA5CIGkmplHabadVKCtPcXOCWW4DAQLkrc2/79gFLlkgrVz32GODvL3dF5E4YpkRN1LatNPF/505gwwZpgEtCAqfP2JvJBHz7rfTf4aabgAkTpBHYRPbEMCVqhoAAYOhQaYH8AwekEb99+gBardyVuYeCAmkOqV4PDBsmXdPmhxmSA8OUqJlUKqBHD2k/zN27gZ9+kkaPtm0rd2Wubd8+YPlywMcHmDwZaNdO7orInTFMiVpIZCRw113Smr47dwKnT0v7ZAYEyF2Za6moAFaulDb1vukmYPx4LlhP8mOYErUgjQbo2xeIi5NG+373nbSub9eu0pQaap5jx4BVq6RdXyZNklakYrcuOQL+701kA1FRwN13S2/+aWnSlI1evThVo6l0OmmQ0ZEj0mv42GPc9YUcC8OUyEZUKqlFGh8vXd/bvVta2u6mmxgE1jKZgF9/Bb7/XmrZjxkD9O7N1ig5HoYpkY35+UkLPVy8KM1L3bxZCtPu3YHgYLmrc0xCSNee16+XFqvv1Utq6fv4yF0ZUd0YpkR2Eh0tdf+eOQMcPAisWyd1WXbpAoSFyV2dYzCbgRMnpJbohQtSy374cGlwF5EjY5gS2ZFCIU2ZiY8HMjKkbt/vvgPCw6X9U2NjpRWW3I0Q0muxebO0olRYGPDkk9x7lJwHw5RIBkol0L69NDcyK0vaL3XLFsDXV9pzMyFBGhns6qqrpa7vbdukhek7dwYefBDo2JHXRcm5MEyJZKRQSK3R2FhpNZ+jR6VwSUuTWmcdOgAxMa7VWhVCun68b5/UpZufL3XnjhsnTSkickYKIYSQuwhHptfrERAQAJ1OB3+unE12UFEB/PGHdCsokJYmbNtWasmGhztvi02vBw4dkkL00iVpIfpevaTlF1u1krs6orpZmwFsmRI5GC8voFs36VZQAJw6Jd0yMqSfx8VJLdnoaECtlrXUBgkhXf88fhxIT5e6cbVaaTDR8OFSq9uVWtzk3himRA4sOFi6JScDOTnS4g9nz0oBpVIBoaFSOEVESDc5l9UTQgr/zEzplpsrjchVq6VrwDffLF0L5fQWckXs5r0OdvOSI9LppIFLly4B585JXcOA1HUaESGFbHAwEBQk3dfSXcNGoxScubnA5cvSLStLqkuhkAK+XTtpNG7btlxKkZwXu3mJXFhAgDQ/tUsX6Xu9Xmq55uRI3akZGdJiByqVtIpQQIA0UtjLC/D2/vOmVkvHKJV/drmazdK51dXSzWiU1sItKpKeR6+Xjisu/rOWsDCp9RwdLXVBu8NIZKIrMUyJXIC/v3Tr0EH6XgigpAQoLJRCT6eTvi8rk0bPlpdLrVmTSQrUykrpPJVKCtIaKpX0uBqNtOF2cLA0R7ZVK6nVGxbGvVuJAIYpkUtSKP4M2PoI8WfL02yWvhdCaqF6ekohywFCRNZhmBK5KYVCanGyS5ao+fi5k4iIqJkYpkRERM3EMCUiImomhikREVEzMUyJiIiaiWFKRETUTAxTIiKiZmKYEhERNRPDlIiIqJkYpkRERM3EMCUiImomhikREVEzMUyJiIiaiWFKRETUTAxTIiKiZmKYEhERNRPDlIiIqJkYpkRERM3kNGFaWFiI8ePHw9/fH4GBgZgyZQpKS0sbPGfgwIFQKBS1bjNmzLBTxURE5C485C7AWuPHj8elS5ewadMmGAwGTJo0CdOmTcOKFSsaPG/q1Kl49dVXLd97e3vbulQiInIzThGm6enp2LBhA/bv349evXoBAD788EMMGzYMb7/9NqKiouo919vbGxEREfYqlYiI3JBTdPPu3r0bgYGBliAFgJSUFCiVSuzdu7fBc5cvX46QkBAkJiZizpw5KC8vb/D4qqoq6PX6WjciIqKGOEXLNCcnB2FhYbXu8/DwQKtWrZCTk1PveePGjUNsbCyioqJw5MgRPPfcczh58iTWrFlT7zlz587FK6+80mK1ExGR65M1TJ9//nn861//avCY9PT0Jj/+tGnTLP/u0qULIiMjMXjwYGRkZKBt27Z1njNnzhw89dRTlu/1ej1iYmKaXAMREbk+WcN09uzZeOihhxo8pk2bNoiIiMDly5dr3W80GlFYWNio66HJyckAgNOnT9cbphqNBhqNxurHJCIikjVMQ0NDERoaet3j+vTpg+LiYqSmpqJnz54AgF9++QVms9kSkNY4dOgQACAyMrJJ9RIREdXFKQYg3XjjjbjjjjswdepU7Nu3Dzt37sRjjz2GsWPHWkbyXrx4ER07dsS+ffsAABkZGXjttdeQmpqKs2fP4rvvvsOECRPQv39/dO3aVc5fh4iIXIxThCkgjcrt2LEjBg8ejGHDhuHWW2/FwoULLT83GAw4efKkZbSuWq3G5s2bMWTIEHTs2BGzZ8/GX//6V3z//fdy/QpEROSiFEIIIXcRjkyv1yMgIAA6nQ7+/v5yl0NERHZkbQY4TcuUiIjIUTFMiYiImolhSkRE1EwMUyIiomZimFpJX8k1eomIqG4MUyvllefJXQIRETkohqmVcktz5S6BiIgcFMPUSvpqPSoMFXKXQUREDohhaiUFFMgtY+uUiIiuxTC1UrBXMHJK6987lYiI3BfD1EphPmEoKC+AwWSQuxQiInIwDFMrhfmGQUDgctnl6x9MRERuhWFqJa2HFoHaQHb1EhHRNRimjRDhG4HLZZdhMpvkLoWIiBwIw7QRwn3CYRIm5Jfny10KERE5EIZpI/hp/ODj6cMpMkREVAvDtJGi/KJwuewyzMIsdylEROQgGKaNFOkXiXJDOUf1EhGRBcO0kfw1/vDX+OOC/oLcpRARkYNgmDZBjH8MckpzuIADEREBYJg2SWv/1hBC4GLJRblLISIiB8AwbQKNhwahPqHs6iUiIgAM0yaL8Y9BYUUhyqrL5C6FiIhkxjBtogjfCHgoPdg6JSIihmlTqZQqtPZvjSxdFoQQcpdDREQyYpg2Q2xALMoMZbhUeknuUoiISEYM02YI0AYg2CsYZ4rOyF0KERHJiGHaTG2C2iC/PB/6Kr3cpRARkUwYps0U6RcJrYeWrVMiIjfGMG0mpUKJ+MB4nNedR7WpWu5yiIhIBgzTFhAXGAcBgSxdltylEBGRDBimLUDjoUG0XzTOFJ3hNBkiIjfEMG0hbYLaoNxQzo3DiYjcEMO0hQR5BSFIG4SMwgy5SyEiIjtjmLagNkFtkFeeh5KqErlLISIiO2KYtqBo/2hoVBpOkyEicjMM0xakVCjRJqgNzus5TYaIyJ0wTFtYfFA8zMKMUwWn5C6FiIjshGHawtQqNdoEtUFGUQYMJoPc5RARkR0wTG2gXat2EELgdOFpuUshIiI7YJjagNZDi7jAOJwuPA2j2Sh3OUREZGMMUxvpENwBJmHiyF4iIidm7doBDFMb8fL0QmxALE4VnILJbJK7HCIiaoKNGRutOo5hakMdgjvw2ikRkZPKK8vD79m/W3Usw9SGfNQ+iPaPxon8E5x3SkTkZH46/RM6BHew6liGqY11Cu0EAYH0vHS5SyEiIitl6bKw+cxm9IrqZdXxDFMb03ho0DGkIzKKMlBaXSp3OUREdB1CCHxx5AtE+kZiYNxAq85hmNpB+1btoVFpcPTyUblLISKi6zhw6QDS89Mxrss4qJQqq85hmNqBSqlCYlgiLugvoKC8QO5yiIioHkazESvSVqBrWFd0i+hm9XkMUzu5IeAGBGgCcCT3iNylEBFRPX4+/TPyK/Ixrsu4Rp3HMLUThUKBruFdUVBRgIv6i3KXQ0REV9FV6rD2xFqkxKcg2j+6UecyTO0o3Dcc4T7hSLucBrMwy10OERFdYXX6aqiUKvzlxr80+lyGqZ11De+K0upSLjNIRORAzhWfw7az2zDqxlHwVfs2+nyGqZ0FaAMQFxiH9Lx0btFGROQAhBBYnrYckb6RGBw/uEmPwTCVQefQzjCajThZcFLuUoiI3F7qpdRGT4W5GsNUBl6eXugQ3AGZRZkoqy6TuxwiIrdlNBvxZdqX6BberVFTYa7GMJVJQkgCFAoFUi+lyl0KEZHbaupUmKsxTGXiofRAj8geyCnNwQX9BbnLISJyO1dOhYnyi2rWYzFMZRTlF4Vov2gcvHSQg5GIiOzsm+PfNHkqzNUYpjLrHtkdBrOB6/YSEdnRueJz2H5ue5OnwlyNYSozb09vdA7tjNOFp1FUUSR3OURELu/KXWGaOhXmagxTB9A+uD1aebfC/uz9XBmJiMjGUi+l4kTBCYzvOr7JU2GuxjB1AEqFEj0ieqC4spjdvURENlRlrMJPp35Cr8he6BretcUe12nC9J///Cf69u0Lb29vBAYGWnWOEAIvvvgiIiMj4eXlhZSUFJw6dcq2hTZRkFcQEsMSkZ6XjvzyfLnLISJySavTVyNLl4WxiWNb9HGdJkyrq6sxevRoPPLII1af89Zbb+GDDz7Axx9/jL1798LHxwdDhw5FZWWlDSttuhtDbkSwdzD2XNgDo9kodzlERC7leN5xbD+7Hfd3uR/hvuEt+thOE6avvPIKnnzySXTp0sWq44UQeP/99/HCCy/gnnvuQdeuXfHZZ58hOzsba9eutW2xTaRQKHBz65tRaazEwUsH5S6HiMhllFaXYlHqIrRr1Q6D4ga1+OM7TZg2VmZmJnJycpCSkmK5LyAgAMnJydi9e3e951VVVUGv19e62ZOv2hfdI7ojoygD2SXZdn1uIiJXJITAJwc+gdFsxNSeU6FQKFr8OVw2THNycgAA4eG1m/Lh4eGWn9Vl7ty5CAgIsNxiYmJsWmdd2rZqiyi/KOy7uA9Vxiq7Pz8RkSvZenYrDuUewpQeUxCoDbTJc8gaps8//zwUCkWDtxMnTti1pjlz5kCn01lu58+ft+vz1+gd3RtCCOzP3i/L8xMRuYKL+otYeXQlBscPRlJEks2ex8Nmj2yF2bNn46GHHmrwmDZt2jTpsSMiIgAAubm5iIyMtNyfm5uLpKSkes/TaDTQaDRNes6WpPXQoldUL+w8vxNni88iLjBO7pKIiJxKtakaH/3+EcJ8wnBf5/ts+lyyhmloaChCQ0Nt8tjx8fGIiIjAli1bLOGp1+uxd+/eRo0IllNMQAziSuKQmp2KUO9Q+Kh95C6JiMhpfH3sa+SW5uLlgS9DrVLb9Lmc5pppVlYWDh06hKysLJhMJhw6dAiHDh1CaWmp5ZiOHTvi22+/BSCNjJ01axZef/11fPfdd0hLS8OECRMQFRWFkSNHyvRbNF7PyJ7wVHli78W9EELIXQ4RkVM4nHMYWzK3YGziWET7R9v8+WRtmTbGiy++iE8//dTyfffu3QEAW7duxcCBAwEAJ0+ehE6nsxzz7LPPoqysDNOmTUNxcTFuvfVWbNiwAVqt1q61N4enyhM3t74Zv2T+gj8K/kBCSILcJRERObTiymIsPrAYSeFJuC3+Nrs8p0KwudMgvV6PgIAA6HQ6+Pv7y1bHwUsHkVGUgcHxgxHkFSRbHUREjkwIgbd3vY3skmy8dttrzd4RxtoMcJpuXnfXNbwr/NR++C3rN+59SkRUj58zfsbx/OOY2nNqi2ytZi2GqZNQKVXoG9MXlcZK7LmwR+5yiIgcztnis/jm+De4s92d6BTaya7PzTB1In4aP9zc+mac15/HiXz7zr8lInJkVcYqfPz7x2jt3xqjbhxl9+dnmDqZmIAY3BhyIw5eOoi8sjy5yyEicgjL05ajuLIYM3rNgIfS/mNrGaZOqFtEN4T7huPXrF9RVl0mdzlERLLad3Effs36FQ90fQARvhGy1MAwdUJKhRJ9Y/pCqVBi69mtHJBERG4rpzQHa9PX4taYW3FLzC2y1cEwdVJaDy0GxQ1CWXUZdpzbAbMwy10SEZFdVRgq8OHeD6FUKjGuyzib7AZjLYapEwvQBmBA3ADklOZg/0UuiE9E7sMszFiYuhC6Kh0e7/04vDy9ZK2HYerkInwjcHPrm3Gq8BTS89LlLoeIyC6+P/k9juUdw4xeMxDuG379E2zMaZYTpPq1bdUWJdUlSL2UCl+1L2IC7L8HKxGRvfyS+QvW/7Ee47uOR2JYotzlAGDL1GV0C++G2IBY/Jb1GwrKC+Quh4jIJg5cOoCVR1fi9ra3Y2DcQLnLsWCYugiFQoG+MX0R5BWEbWe3ccoMEbmcUwWnsCh1EXpF9cLoTqPlLqcWhqkLUSlVGBg3ECqlilNmiMilZJdkY96+eWjXqh0md58s68jdujBMXcyVU2Z+zfqVe6ASkdMrrCjE+3veRyuvVnj0pkdlWeHoehimLqhmysylkkvYn80pM0TkvMoN5fjPnv9AqVDiiZufkH0KTH0Ypi4qwjcCya2T8UfBH5wyQ0ROyWAyYN6+edBV6TDr5lkI1AbKXVK9HK+tTC2mXat2KKkqQWp2Krw8vRAXGCd3SUREVjELMxYfWIyzxWcxu89s2dbctRbD1MUlRSSh0liJ37J+g1qlRpRflNwlERE1SAiBlUdX4mDOQcy8aSbatmord0nXxW5eF6dQKJDcOhlRflHYmrkVl8suy10SEVGDfjr9E7ZmbsWDXR9Et4hucpdjFYapG1AqlBgQOwChPqH4JfMX7oNKRA5r1/ld+Db9W4xIGIF+sf3kLsdqDFM3oVKqMChuEEK8Q7AxYyNyS3PlLomIqJa03DR8dvgz9IvthxEdRshdTqMwTN2Ip8oTA2IHWAI1uyRb7pKIiAAARy8fxce/f4y+MX3xQNcHHG5RhuthmLoZT5UnUtqkIMovCpvPbEaWLkvukojIzR29fBQf7f8IncM6Y1yXcVAqnC+anK9iajaVUoVB8YMQGxCLrZlbkVGYIXdJROSmrgzSaT2nOeTqRtZwzqqp2ZQKJfrH9oeH0gO/Zv0Ko9mIhJAEucsiIjdSE6SJYYmY2nOq0wYpwDB1azU7zXiqPLH7wm4YzAaH2RuQiFzb4ZzDWHhgoUsEKcAwdXsKhQK9o3vDU+mJ37N/h8FkQPfI7nKXRUQubGfWTixPW45bYm7BmMQxTh+kAMOU/qd7ZHd4qqRANZqNuCn6JrlLIiIXI4TAhtMb8N3J79A/tj/GJI5xysFGdWGYkkViWCI8lB7Yc2EPDGYD+rTu43TD04nIMZmFGV8f+xrbz27HiIQRuLPdnS71/sIwpVo6hnSEp9ITv2X9BqPZiFtvuNVlPjkSkTwMJgOWHVqGgzkHMb7reNx6w61yl9TiGKZ0jbat2sJD6YHt57bDaDZiQOwAqJQqucsiIidUYajAx79/jMziTEzvOd1p1tptLDY5qE6xgbEYHD8YF/UXsSVzC4xmo9wlEZGT0VXq8O7ud3Fefx5/S/6bywYpwDClBkT7R+P2trcjvzwfGzM2otJYKXdJROQkcktz8e9d/0aZoQzP9H0G7Vq1k7skm2KYUoMifCMwtO1Q6Cp1+P7k9yiqKJK7JCJycGeLz+LtXW/DU+mJZ/o+g0i/SLlLsjmGKV1XsHcwRiSMgJenF9adXMflB4moXkcvH8V7u99DmE8Ynu77NIK8guQuyS44AIms4qv2xR3t7sBvWb9h69mtyC3LRXJ0MgcmEREAaQ7ppjObsO7EOtzc+maMSRwDtUotd1l2wzAlq3koPTAwbiDCfcKx+8Ju5JXlYXCbwfBV+8pdGhHJqNpUjc8Pf44Dlw7gjnZ3YHiH4S41h9QaCiGEkLsIR6bT6RAYGIjz58/D399f7nIcRn5ZPrad2wajyYh+sf0Q7R8td0lEJIOC8gIsObgEBeUFGNd1HLqGd5W7pBal1+sRExOD4uJiBAQE1Hscw/Q6Lly4gJiYGLnLICIiGZ0/fx6tW7eu9+cM0+swm83Izs6Gn5+fU3Rb1HyKYku6YXydrMPXyTp8nazjjK+TEAIlJSWIioqCUln/mF1eM70OpVLZ4KcRR+Xv7+80f6xy4utkHb5O1uHrZB1ne50a6t6twakxREREzcQwJSIiaiaGqYvRaDR46aWXoNFo5C7FofF1sg5fJ+vwdbKOK79OHIBERETUTGyZEhERNRPDlIiIqJkYpkRERM3EMCUiImomhqkL+Oc//4m+ffvC29sbgYGBVp0jhMCLL76IyMhIeHl5ISUlBadOnbJtoTIrLCzE+PHj4e/vj8DAQEyZMgWlpaUNnjNw4EAoFIpatxkzZtipYvuYP38+4uLioNVqkZycjH379jV4/KpVq9CxY0dotVp06dIFP/74o50qlVdjXqdly5Zd83ej1WrtWK08duzYgREjRiAqKgoKhQJr16697jnbtm1Djx49oNFo0K5dOyxbtszmddoCw9QFVFdXY/To0XjkkUesPuett97CBx98gI8//hh79+6Fj48Phg4disrKShtWKq/x48fj2LFj2LRpE9avX48dO3Zg2rRp1z1v6tSpuHTpkuX21ltv2aFa+/jqq6/w1FNP4aWXXsKBAwfQrVs3DB06FJcvX67z+F27duH+++/HlClTcPDgQYwcORIjR47E0aNH7Vy5fTX2dQKkVX6u/Ls5d+6cHSuWR1lZGbp164b58+dbdXxmZibuuusuDBo0CIcOHcKsWbPw8MMP4+eff7ZxpTYgyGUsXbpUBAQEXPc4s9ksIiIixL///W/LfcXFxUKj0Ygvv/zShhXK5/jx4wKA2L9/v+W+n376SSgUCnHx4sV6zxswYIB44okn7FChPHr37i1mzpxp+d5kMomoqCgxd+7cOo+/7777xF133VXrvuTkZDF9+nSb1im3xr5O1v6/6MoAiG+//bbBY5599lnRuXPnWveNGTNGDB061IaV2QZbpm4oMzMTOTk5SElJsdwXEBCA5ORk7N69W8bKbGf37t0IDAxEr169LPelpKRAqVRi7969DZ67fPlyhISEIDExEXPmzEF5ebmty7WL6upqpKam1vo7UCqVSElJqffvYPfu3bWOB4ChQ4e67N8N0LTXCQBKS0sRGxuLmJgY3HPPPTh27Jg9ynUqrvT3xIXu3VBOTg4AIDw8vNb94eHhlp+5mpycHISFhdW6z8PDA61atWrwdx43bhxiY2MRFRWFI0eO4LnnnsPJkyexZs0aW5dsc/n5+TCZTHX+HZw4caLOc3Jyctzq7wZo2uuUkJCAJUuWoGvXrtDpdHj77bfRt29fHDt2zCk3zrCV+v6e9Ho9Kioq4OXlJVNljceWqYN6/vnnrxnAcPWtvv+R3YmtX6dp06Zh6NCh6NKlC8aPH4/PPvsM3377LTIyMlrwtyBX06dPH0yYMAFJSUkYMGAA1qxZg9DQUCxYsEDu0shG2DJ1ULNnz8ZDDz3U4DFt2rRp0mNHREQAAHJzcxEZGWm5Pzc3F0lJSU16TLlY+zpFRERcM1jEaDSisLDQ8npYIzk5GQBw+vRptG3bttH1OpKQkBCoVCrk5ubWuj83N7fe1yQiIqJRx7uCprxOV/P09ET37t1x+vRpW5TotOr7e/L393eqVinAMHVYoaGhCA0Ntcljx8fHIyIiAlu2bLGEp16vx969exs1ItgRWPs69enTB8XFxUhNTUXPnj0BAL/88gvMZrMlIK1x6NAhAKj1IcRZqdVq9OzZE1u2bMHIkSMBAGazGVu2bMFjjz1W5zl9+vTBli1bMGvWLMt9mzZtQp8+fexQsTya8jpdzWQyIS0tDcOGDbNhpc6nT58+10ytctq/J7lHQFHznTt3Thw8eFC88sorwtfXVxw8eFAcPHhQlJSUWI5JSEgQa9assXz/5ptvisDAQLFu3Tpx5MgRcc8994j4+HhRUVEhx69gF3fccYfo3r272Lt3r/jtt99E+/btxf3332/5+YULF0RCQoLYu3evEEKI06dPi1dffVX8/vvvIjMzU6xbt060adNG9O/fX65focWtXLlSaDQasWzZMnH8+HExbdo0ERgYKHJycoQQQjz44IPi+eeftxy/c+dO4eHhId5++22Rnp4uXnrpJeHp6SnS0tLk+hXsorGv0yuvvCJ+/vlnkZGRIVJTU8XYsWOFVqsVx44dk+tXsIuSkhLL+w8A8e6774qDBw+Kc+fOCSGEeP7558WDDz5oOf7MmTPC29tbPPPMMyI9PV3Mnz9fqFQqsWHDBrl+hSZjmLqAiRMnCgDX3LZu3Wo5BoBYunSp5Xuz2Sz+8Y9/iPDwcKHRaMTgwYPFyZMn7V+8HRUUFIj7779f+Pr6Cn9/fzFp0qRaHzgyMzNrvW5ZWVmif//+olWrVkKj0Yh27dqJZ555Ruh0Opl+A9v48MMPxQ033CDUarXo3bu32LNnj+VnAwYMEBMnTqx1/Ndffy06dOgg1Gq16Ny5s/jhhx/sXLE8GvM6zZo1y3JseHi4GDZsmDhw4IAMVdvX1q1b63wvqnltJk6cKAYMGHDNOUlJSUKtVos2bdrUep9yJtyCjYiIqJk4mpeIiKiZGKZERETNxDAlIiJqJoYpERFRMzFMiYiImolhSkRE1EwMUyIiomZimBIRETUTw5SIiKiZGKZERETNxDAlIiJqJoYpEV0jLy8PEREReOONNyz37dq1C2q1Glu2bJGxMiLHxIXuiahOP/74I0aOHIldu3YhISEBSUlJuOeee/Duu+/KXRqRw2GYElG9Zs6cic2bN6NXr15IS0vD/v37odFo5C6LyOEwTImoXhUVFUhMTMT58+eRmpqKLl26yF0SkUPiNVMiqldGRgays7NhNptx9uxZucshclhsmRJRnaqrq9G7d28kJSUhISEB77//PtLS0hAWFiZ3aUQOh2FKRHV65pln8M033+Dw4cPw9fXFgAEDEBAQgPXr18tdGpHDYTcvEV1j27ZteP/99/H555/D398fSqUSn3/+OX799Vd89NFHcpdH5HDYMiUiImomtkyJiIiaiWFKRETUTAxTIiKiZmKYEhERNRPDlIiIqJkYpkRERM3EMCUiImomhikREVEzMUyJiIiaiWFKRETUTAxTIiKiZvp/s4PjK0SXh9IAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sim = setupSimulation()\n", "sim.collision = \"direct\"\n", "ps = sim.particles\n", "sim.collision_resolve = my_merge # user defined collision resolution function\n", "sim.integrate(100.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we were not only able to resolve the collision, but also to run additional code during the collision, in this case to make a plot, which can be very useful for debugging or logging. Now that you know the basics, you can expand the scenario here and resolve collisions according to the astrophysical problem you are working on." ] } ], "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": 1 }