diff --git a/20180708Fields/Mercury Interconnect Sept 23 2018.ipynb b/20180708Fields/Mercury Interconnect Sept 23 2018.ipynb index 2ffb5ba5d..a5d3f0031 100644 --- a/20180708Fields/Mercury Interconnect Sept 23 2018.ipynb +++ b/20180708Fields/Mercury Interconnect Sept 23 2018.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 87, "metadata": {}, "outputs": [ { @@ -20,7 +20,9 @@ "import ctypes\n", "from ctypes import POINTER, c_double, c_int,byref\n", "%matplotlib inline\n", - "import matplotlib.pyplot as plt" + "import matplotlib.pyplot as plt\n", + "from matplotlib.patches import Circle\n", + "from matplotlib.collections import PatchCollection" ] }, { @@ -45,6 +47,25 @@ "rebound.Simulation.move_to_hel = move_to_hel" ] }, + { + "cell_type": "code", + "execution_count": 134, + "metadata": {}, + "outputs": [], + "source": [ + "def sim_copy(self):\n", + " \"\"\"\n", + " Shallow copy of a simulation.\n", + " \"\"\"\n", + " sim2 = rebound.Simulation()\n", + " for i in range(self.N):\n", + " sim2.add(self.particles[i])\n", + " sim2.t = self.t\n", + " sim2.dt = self.dt\n", + " return sim2\n", + "rebound.Simulation.copy = sim_copy" + ] + }, { "cell_type": "code", "execution_count": 16, @@ -203,10 +224,11 @@ }, { "cell_type": "code", - "execution_count": 55, + "execution_count": 135, "metadata": {}, "outputs": [], "source": [ + "rcrit = [0.0, 0.02, 0.04]\n", "def setup_sim():\n", " sim = rebound.Simulation()\n", " sim.add(m=1)\n", @@ -218,13 +240,12 @@ }, { "cell_type": "code", - "execution_count": 56, + "execution_count": 136, "metadata": {}, "outputs": [], "source": [ "def step_manual(self,algo):\n", " if algo == \"mercuryhybrid\":\n", - " rcrit = [0.0, 0.2, 0.4]\n", " self.mercury_step(rcrit=rcrit, algor=10)\n", " elif algo == \"mercurymvs\":\n", " self.mercury_step(algor=1) \n", @@ -250,12 +271,291 @@ " self.integrate(self.t+self.dt,exact_finish_time=1) \n", " else:\n", " raise ValueError(\"Not defined\")\n", - "rebound.Simulation.step_manual = step_manual" + "\n", + " rebound.Simulation.step_manual = step_manual" ] }, { "cell_type": "code", - "execution_count": 57, + "execution_count": 124, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbsAAAGfCAYAAADGVHw+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXusbNld3/n9rb3rcc65j36624+23cYG26AkHS62edlgG2wkhJ0AwUCISZixIg3KH6MZ0SNmgmQUyWgihhmFBEyIMQyWeQyYnhGTlsFJFCEMvh0zhDbY/bBxP3D3fZ9z7jlVtfdav/ljrbX32rv2rlPnnjr3Uff7QdVV+1FVu+qa+pzfY60lqgpCCCFknTE3+gIIIYSQ44ayI4QQsvZQdoQQQtYeyo4QQsjaQ9kRQghZeyg7Qggha89KZCci7xaRz4vIkyLycMfxt4rIfxGRUkS+r3XMisifhdsjq7geQgghJEWOOs5ORDIAXwDwHQCeBfAZAD+oqp9Lznk1gFMA/gcAj6jqbyfHdlX1xJEughBCCFlAvoLXeBOAJ1X1aQAQkY8DeA+ASnaq+qVwzK3g/QghhJBDsQrZvRzAM8n2swDefIjnj0XkLIASwIdU9RNdJ4nIBwB8AAC2tra+/vWvf/01Xi4hhJBbhccee+y8qt571NdZheyOyqtU9TkReQ2AT4nIf1XVp9onqeqHAXwYAM6cOaNnz5693tdJCCHkOiMif72K11lFg8pzAB5Itl8R9i2Fqj4X7p8G8B8BPLSCayKEEEIqViG7zwB4nYg8KCJDAO8DsFRXpYjcKSKj8PgeAN+MpNZHCCGErIIjy05VSwA/DuBRAH8J4DdV9XER+aCIfA8AiMg3iMizAL4fwC+KyOPh6W8AcFZE/j8A/wG+ZkfZEUIIWSlHHnpwI2DNjhBCbg9E5DFVPXPU1+EMKoQQQtYeyo4QQsjaQ9kRQghZeyg7Qgghaw9lRwghZO2h7AghhKw9lB0hhJC1h7IjhBCy9lB2hBBC1h7KjhBCyNpD2RFCCFl7KDtCCCFrD2VHCCFk7aHsCCGErD2UHSGEkLWHsiOEELL2UHaEEELWHsqOEELI2kPZEUIIWXsoO0IIIWsPZUcIIWTtoewIIYSsPZQdIYSQtYeyI4QQsvZQdoQQQtYeyo4QQsjaQ9kRQghZeyg7Qgghaw9lRwghZO2h7AghhKw9lB0hhJC1h7IjhBCy9lB2hBBC1h7KjhBCyNpD2RFCCFl7KDtCCCFrD2VHCCFk7aHsCCGErD2UHSGEkLWHsiOEELL2UHaEEELWHsqOEELI2kPZEUIIWXsoO0IIIWsPZUcIIWTtoewIIYSsPZQdIYSQtYeyI4QQsvZQdoQQQtYeyo4QQsjaQ9kRQghZeyg7Qgghaw9lRwghZO2h7AghhKw9lB0hhJC1ZyWyE5F3i8jnReRJEXm44/hbReS/iEgpIt/XOvZ+EXki3N6/iushhBBCUo4sOxHJAPw8gO8C8EYAPygib2yd9mUAPwrgY63n3gXgpwC8GcCbAPyUiNx51GsihBBCUlYR2b0JwJOq+rSqzgB8HMB70hNU9Uuq+ucAXOu57wLwSVW9qKqXAHwSwLtXcE2EEEJIxSpk93IAzyTbz4Z9K32uiHxARM6KyNlz585d04USQgi5PbllGlRU9cOqekZVz9x77703+nIIIYTcQqxCds8BeCDZfkXYd9zPJYQQQpZiFbL7DIDXiciDIjIE8D4Ajyz53EcBfKeI3BkaU74z7COEEEJWxpFlp6olgB+Hl9RfAvhNVX1cRD4oIt8DACLyDSLyLIDvB/CLIvJ4eO5FAD8NL8zPAPhg2EcIIYSsDFHVG30Nh+bMmTN69uzZG30ZhBBCjhkReUxVzxz1dW6ZBhVCCCHkWqHsCCGErD2UHSGEkLWHsiOEELL2UHaEEELWHsqOEELI2kPZEUIIWXsoO0IIIWsPZUcIIWTtoewIIYSsPZQdIYSQtYeyI4QQsvZQdoQQQtYeyo4QQsjaQ9kRQghZeyg7Qgghaw9lRwghZO2h7AghhKw9lB0hhJC1h7IjhBCy9lB2hBBC1h7KjhBCyNpD2RFCCFl7KDtCCCFrD2VHCCFk7aHsCCGErD2UHSGEkLWHsiOEELL2UHaEEELWHsqOEELI2kPZEUIIWXsoO0IIIWsPZUcIIWTtoewIIYSsPZQdIYSQtYeyI4QQsvZQdoQQQtYeyo4QQsjaQ9kRQghZeyg7Qgghaw9lRwghZO2h7AghhKw9lB0hhJC1h7IjhBCy9lB2hBBC1h7KjhBCyNpD2RFCCFl7KDtCCCFrD2VHCCFk7aHsCCGErD2UHSGEkLWHsiOEELL2UHaEEELWHsqOEELI2rMS2YnIu0Xk8yLypIg83HF8JCK/EY7/iYi8Oux/tYjsi8ifhdsvrOJ6CCGEkJT8qC8gIhmAnwfwHQCeBfAZEXlEVT+XnPZjAC6p6mtF5H0AfgbAD4RjT6nq3znqdRBCCCF9rCKyexOAJ1X1aVWdAfg4gPe0znkPgI+Gx78N4B0iIit4b0IIIeRAViG7lwN4Jtl+NuzrPEdVSwBXANwdjj0oIp8Vkf8kIt/a9yYi8gEROSsiZ8+dO7eCyyaEEHK7cOQ05hH5GwCvVNULIvL1AD4hIl+rqtvtE1X1wwA+DABnzpzR63ydhNwUWGtRlg6lVVinKJyDg8I5//8SqkBMmuRGYESQZwZGgMwYGCPIjMAYAZMr5HZiFbJ7DsADyfYrwr6uc54VkRzAaQAXVFUBTAFAVR8TkacAfDWAsyu4LkJuefb2ptibWezOSkzLAmXhUFoH5xQOgP9/IY+IIIuiy8zcffs2CLc8M8iMIMvYnE3Wl1XI7jMAXiciD8JL7X0Afqh1ziMA3g/gjwF8H4BPqaqKyL0ALqqqFZHXAHgdgKdXcE2E3JLMZjNc3pvhwtUCu3sTTK1DUSqsKqw25ZYiIjAS7gEYAYz4CC4XeJmJF1pmguxyfz/MMwzC/SjsY+RH1o0jy05VSxH5cQCPAsgA/DtVfVxEPgjgrKo+AuCXAfyaiDwJ4CK8EAHgrQA+KCIFAAfgn6rqxaNeEyG3EmVZ4itXdvGVKyUu7e5hFlKU1ikUgHMKVaBLc8YIpONxlaoMjzPjJZgbad4yg0GQXZ5nGOQ+2hsPcoxyg+Ego/TIWiB9fynezJw5c0bPnmWmk9zavHDhMp46P8WFq/uYFq4SnAsRnOsRXEQASIjm0scmuY/RXXwcxTcnPSM+uhvklfCGQXjjQY5hbpjmJDcEEXlMVc8c9XVudIMKIbcV0+kUX3xxF198YRuXZ77+Zp1PUboQyan2R3IRE4KtGHWZDtGl91kQXhaaVmKjShTfwAgGxmFQ2irSG8xKjIY59mclxsO8ivYoPXIrQtkRch24evUq/uL5XXzpwg72pw6FrdOU1mlDcO1siwJIE4neb+IjObSjuy7R+aiuLb14K42giNJzBrlVDIxgmBnMrMMwMyisw7SwGA0ybDDSI7cglB0hx8jlK7v4ixd28MUXtzGZalWPi7U4pzGaA4DDR3RxXyo7kTTSCylM03wcI7xUepkRzKximHn5FU4xtA5FkN0wzzArM8xKh1FusDHMMWJNj9wiUHaEHAO7u7v4ixd28dRXrmB7zzWbTrROV4bhcQ3Z9SEN0fkTo9QATR77yE2gVb1OXJSbf9yUXRrpGVgnGGSC0inKzCB3itI5FNanXUvrUA5zf9w6jId5NcSBkJsVyo6QFbK/v4/HX9jGF7+yjfO7BfZnDtYhNJ6EyE3jAPB5uS0aWpCcBYQUpguik0R6mRFYpz51Ge8l7vcpTRvkZl1bdloJLzNeeF58BqMgt9I62Cg86yXIKI/c7FB2hKwAay2e/so2/ur5i3hhd4adfYtpmURyleiaQqvEd+A7+DMk/LdyiigkbosXm3Ma0pnaSGc61ZDKVBhXC8+pl6CXHpAZV8lP1R8vDaAq1Vg/qyFSHThYVZQhat0Y5jCGwiM3H5QdIUfk/OVdPP7MRXz5ygTbuzPszPx0XtX4OK1TlwAaj5HsA+alJ8mjOmjSENzFfVrV61xIXfpoT+AkNrAoMvWiq5pWVKsILkqvfYuyyyrZaSPKsxo+59BHrk698JjWJDcblB0h18hsNsPnvnIFX/rKNi7sFbi8M8VeqXC2Kbj6Hkg11yW9fqLgAKSSS+bCFA2Csz6Ci8IT8QJMI7tacvFapSG+OCRBs6bs2lEeUEBDitZq3V26OcwxHGSr+aIJWQGUHSHXwIXLu3j8yxfx3O4UV67OcHm/wGSm0NhpGSK1eL8osquaU3reay512Y7stE5lVsILYkO4j9sqCicCY/z7uuo4gtBq+Xm5GeQZ4NRVssszfz0uXJdq6UWn4fMP69QtZ2AhNwuUHSGHwDmHJ164hKee38aLV2fY2Stweb9AUdSRnGtEdbXG2sIDlqnVRdLIrm5GqSO7+nFbeBq2FT6lWYnOKIz6bk3fGRpre15+/ubg1HdnqprwmZx/nNzb8GFjObL6/AAbV8hNAWVHyJJMJhM8/vw2vnR+F5f3ZtjeK3Blv0BZ9omueQ90RXXLSa8pOYRoLmwnj6Po4r1W275uZyBQqxAjMCpQ8dJzKshCdBdevrpFnDr4SM7Ar7kgYehEjPIAoKw/W3L9FB650VB2hCzBxStX8blnL+H5kLa8sldge79AaZcT3bINKkBTEnMNKlp3ZR4kOukSXdhnXJiZJUR3aEjO1/jCK8MLzt8A05DeIHMAfHRXf5Yyuer6Nhry54bcOPi/PkIO4PkLu/jC8xfxN7szbF+dYXtSYme/QFG6agaUVHTtOl2f8AC0orz2nohvMomHJNnXlF4cblDX6qLwGqJDXOTVD0GA1JJzoSszSyK7mNIMVchKel5wBn7BkngfrrFIhBemKhMRNq2QGwZlR8gCnnrhEr74lW28uDvDzr6P5nYmJaaFrcTWJbq++4MjPKBrAEKaAUwbVirZIcx/qX6geRxwHoVXPwZUagHG6zDiIzmjqKK8LuoAtEt0/r66VLFh2EOsFYb19TgsgdwAKDtCOrDW4okXLuPLL+7i/N7MS26/xM6kxKSwnZJbLp2ZNqwA7S7MPtkByXRh4b+dskOcJ7OO8trCU8ToTqEh4lIAxmcnG/W69i3ia3UdohMDsV54WZzVRWxlPAGwNRpw4Dm57lB2hLSw1uILL17BM+d2cWHfd1xuT0rsTgrsz+o2+74UZvOx11lXV2a6D2iKbr5ut0h4zXtTyQ617ODH16lIJT0TJZe8v8Y0KPpu7ajMi04kFZ3xKVbrIFLCAH6VdBFMwvtujnI2rJDrCmVHSIK1Fk+8eAXPvrCDi/sFtq+W2N4vcHVSYjIrUTrXSF+25dYV1aXSA/plBzSjpyahMaWVzqwGlLdEl04T5relGlDux9vN32LUpwA0NLCkN0+M4jwmzM4Sj4WYESIOJkgPsBAj1c1UC8WyfkeuH5QdIQFVxdPnruDZczs4v19gZ6/AzqTA1anFpLAoXJgeK5Ga6xFdX4Tn3+co15g8buzQzujOSHqLkZ4JjxHWtquF1ojbHOZsp9UB0zrRVOITcShsKj0gK2wV3WUimGSG9TtyXaHsCAl88dxlPHtuFxdDJHd16m+TWYlZWaJMRNeWXC0/L7V2CrOL5aO6edoZQIHUQlXAhWirKbs6usvEJKnMsNxPEuXBJBJy7ffyO0VMeI80uqvvTZCen3jawVgHKS2yzGAa5WeE6UxyXaDsCAHw7PltPPviLi5OSuzsW+xNS+xOS0xnFjNrUbRE1751ia7NspHdUjNmalMQ8TmVN8JxqwoXVkaIcstE4EwtPEUUnamiPDgXxGfaowoqpBKqvwL/2EvQiENpTRjeoDBWkYlFJsAsRHVMZ5LrCWVHbnvOXdnDX7+4jYv7vgnFi85iMikxLUtMrW2kL+cl115x3LNIbvONKYejLcRkYhUf5SXy8/t8FCoCWBFkairhqQo0LPIaxQdj6pSlC5FeFF7VBZoKNkoOydJCXnSFVWTiUNjQqGJ8VJcbn87Mg/wIOU4oO3Jbs7s7wV+/uI0LkxI70wJX9y32ZiX2pyVmzmFqLUqnYY7Ivsiu+Zp9kltWfsnOflpuiJKLr9WO8Kp9YZkeJxZODZwonDFwEGRBdip+OrAsRnUA4FwtPAEs0tRlLb+27ErbTGdmziGzDrl1KEqHWWYxzQw2OLsKOWb4vzBy21IUBb54fgcXrk6xNy2xP7O4OiuxP7WYFRazspxrSokrjvtITuciuS7B9UkvbPSyMJ0ZDkliPa13dkZ4UKmk58fNOVgBclUvPqN+3ksTXsSgjuqAhvAEgE1kNxfRVfd1OjNzisJpJbsovCJ3GFjHZhVyrFB25LblmUu7OLe9j+3QiLIXui5nhcXUliGqqyO6Zidm/Tpd0gO6U5ptfy0rvj4az6+iuaboavmlaU1BmE4FRYhao+gcFIqwwoGEgl0QnqhWskvTmNZJj+wQGlT8auaZU5TOoUiENy0sBplhswo5Vig7clvywuU9vHBxHzszi/3CYn9qsT8rMZlaTK1FYeuhBraRxqxf46DIrktwfdJrHKt3LEciOb/ZFF0qP02OV0PFBbAKOGeRB+nF4QY+ykuEF6I7cfUbOvHdn9Z1RXiCzPjvMMquiM0qYehBbh1mpRdenlF25Hig7Mhtx/7+FM9d2MGlSVGlLfcLi8nMeckF0cVanW1Fc0tJ7jBR3bI1vD76TpXWsXQ7jfaSW5HUJpuDyWvhiXpRwsWGl2ZnZio7LzqE1dGD8MTfl9ZVtyLcmMokxwVlR247vnRpD1f2Zrg6s9ibWuyHel1R+mEG03BfR3VJ5yQWS29OckeU3oLdc9QB3uHyoYpkMrAkygPQWHoIcXaUKq3pIM4E4XnxORVYp8iMX2XBumYq0zpFGWRnncJG4TlFYR1GTtmZSY4Fyo7cVlzY3selK3vYnvqOy2lhsV8635BSRXW+A9Oqq9KWi6K5pSS3lPT6Ng7cPRfEpTsECwQYz5O6xhjFN3O2sVpDdUzq6E5UIc5HcS5IzwgqybVTmdYprInNPm4uuiutw9Bw3B1ZPZQduW0oigLPX9zF7swLblo6TAqH6cyGNJpPXxbWouxJWy4rubbQuut37Qcdm0sGaYeP5boJIwuqKLFQF2p9WVhNAfVqCzGlGWRnWynMdiozMwrn/C2Kr7Sull6I7ga5YaMKWTmUHblteP7yHq7sTXG1sJjMLKYzi0nh29/LI4juMJLrEtzCfb07FjAX5i0gFu5iLjO8mSZrLZQO8CPrAEFcrUD9PJgxwgvpR1/vq+9jKjMzUgtPQ80uRM+xfuesQ+kcrFM2qpCVQ9mR24LJZIaLO1Nsz3zX5TR0Yc4K24joClenLtNaXZXOOyia6436cGjhdWxe0/CE5agvTkPkFqM8AChcTGPWEvSTTSvE+HSmqxpWQipTmrU7H+FJlcq0qigdUFoHjQ0rIfIDM5lkxVB25Lbg+e097ExmmBQW08L58XRliOpS0YUf9VRuS0VznY/ReKCNnT3H5w72+O2QzZrS3o51OrQDO5070fm1FFA4B6BsrpMngFGTyE46JNcUnm9i8Y0wVrUSXWwIYiqTHAeUHVl79vamuLwzDUv1+DrdzDrMCofCJWPqnMJBO5tSuoR2LZLrFVyy/bnzh/t8b7hryRPb1mvvT6I6n9706cxGDc85oCzrBWPDKgoShiuYkAKel5x/HOt1jfqd1qKzQXqq8ys7EHIUKDuy9rx4dd+vMl5aTIs4FVis1blKdHHlAqBbdEtHc8s8bm0fVnApf3mxfvz6u9I3qDlUdJdIT5IXi/NAl85han2EVy8VFGaXCdGdmq6oTvzKEGECbasKq74m6JzCBdHFGh67MskqoezIWrO/P8Olbd+UMiudl5zTKqornfWiQ0ipoSW6jlpdt/xwTcI7iuS6+KtEfF9zZ/NYKrhqu3EQregu7BMNN4EL04bNrIURQW7qdfGi8EyrOSUOPXBO4UxsYEG4BfHZ+eiOkFVC2ZG15sLeFHuzsk5fFn4uxqK0oU7kz7POhZXF+0V3UDR3IyXXxecvAV99R3NfW3hRaHF/Kro6uqujPOMcXFgNYWptYxHYSnguik2SJZCaUV1jCEJSu9M48XaItFm3I6uCsiNrS1mWuLQ7wd4szIpS2KpWZ0MazafTfILOHVJ0XUI7SHifu3C8n7nNFy77+68+Xe+rrkua8quCu1Avq/YLQlOKnwNTjAOcgYXDpCz9HJfG+DXyEpmlacy29KobmmlNp4A6Zd2OrBxOREfWlu39Elf3Cz9DSukwKxWTwsFaP/O+at2Q0pgtZZHotCk6hZfDIvnF49dbdClfuNKxM163Aur8vXP1dhSSurAd9zmFC3XOwjlMyxJFGBzuXL2obUxjVq+TiFDTc6L0XLK6RGuaNkKOCmVH1pbzu/vYL0JUVypmcbYOZ/2PKeof3jjkoE90Pr0Zj2Gh8DQKJDl+I0UX6RReSvwMrpacS6TnggydS4XnF7idlGWyHFIizSA5l0qvfqu5hhWfyqwjO0JWBWVH1pLJpPBRnfURXZx3sSiSQeOhMaWe4b9fdHVkhznR9R4Lt7+8CUQXOVB4keT6K+mFm7aEV1ifzpxZizKujYc6OqsEl4itvoWJp5P0ZXUJtB1ZIZQdWUt2piX2p4VfCds6lKXDtKiXr3HBUC7U62JjyoGiQ39k15XSTIcF3Cx84Uoz+oy3TjqEZ5N7G4Q3sxaTsvDTfmkzlZmmLhtRYprajNfimMYkxwNlR9aSS1cnmIbJnkvrQgpTwzADrdKYql54S4uuJ7LrSmnejKKLPLEdHiS2awtwDm3V7xzgLGCtwjqfyiyc9bU71SqV6Vqvqa338vvm5UbVkVVC2ZG1YzabYW9aYOLUj6krHQqb1uX8eXUaE1XjiT+hea9d9wcI72YWXaQSXptUfpiXTmxkcTZGd0BZKqalxX5RhOafGC333xD/2EheO9btYnRHyKrg0AOyduxO/YoGs8KitIrShq5Ap9UPcExjxjpdrBXp3D0W34f/tO+vJ/mpz2L00t+CiGvst1e/CvvP/LcLn/vENvC6UwtOSD9nGAZQDdELkV5tRIfd2QxbwyGGmiFTmW9waURz/pWj8LpqdBxrR1YFZUfWjt1pgf2iROE0DB5XWFvLK63PRQEC6LhH9ZyF9+E/8f6vLh3rx6sY3fcJDO78NIDu8WjZ1lPYeOCXDhSeKubmzOzUS0t8cfC50/q7mBUWV2cFxnkONVnVhJLe4mtUKV8k4xtbTSqErArKjqwdu/slZlZRlKFOFxsmUHcHAvPia9brUN9jCeGF++MWXX7qsxjd/zsQUwBYPOhaxAvvIJ7cAV57srmv0k0rmktPSKWnCljr7y9Ppjg9HiWiSwqD6JBeeGy1+fqErBLKjqwVRVFgMp1hFqeisn7YAVDLDQhNKah/aHtF1yO8G5G+3HzwZ2FGL17fWUW0cTcvv0R6UF+/m04sdqb72MgH9R8T1a0WXyo92o0cN2xQIWvF/qTELK5mEOp19Y8s6h/vRHya/KL3SawtvK5zPn/5eD5Tfuqz2Pqah49VdE/uLHlicFXazxP3x+/BWuD8zqQxMLyO69CSXA8s05EVw8iOrBX71UwpYSHQpNEklZuIVNFcemxZ4fUdWzUbD/wSsq2nrklyqr5J5TDnVyxZw0tTmen2ZApc2N/Dy06crNKZ4ZTmC3RAz5HjgJEdWSuulgUKl3RfarLqeIfcBNJMXwKLhXcdU5ibD/7sNYkuRk7LdGOmPLWbvkjztnAMXjvSC49f3L6KSVlUX1C38CL+Q0pYPSE+ZicmWRWM7MhaMZmUYdaU5rg6TX5ipbWv6kZ0Swiv59jS03AtycYDv3SotGWdXs0x+ZvvRbn90GovqHqj5G5B/U4UsCVwcTLB5iDHQIfN62w+EwCQJUuiG3DFA7JaKDuyNlhrURS2iuys1sMK0g6LWnz+salSmjg4wrsOEd3ovk8sHdFVnYyHjOJ6X6+1vfAS+sQXvpuyBKazKS7PRhgOcmBw8MrjAkCMj+iYdiKrhLIja4O1tpqb0bpmt18zugutg+39ywivtW/V5Kc+i8Gdnz5QdBoMM3n+H6w2itOFmw37SceJaQ1PFdgrSpwoZrg6zbE5zueeZ8SnKzNppjCNYQqTrBbKjqwNs5mf9DmuuQYA0BAtiHQKz4T96crdC4V3jKIDgNF9//dSoltVJHdotOPhgnRmUQBFabFfWuxNLIYnfHTn63HzL+/FJxB44RGyKlaSKRCRd4vI50XkSRF5uOP4SER+Ixz/ExF5dXLsfwr7Py8i71rF9ZDbE6uKUusZ94FQ9wm5tbQDE6ijukYa05/aEF6X/I4LyfYWHlcFJs//wLGJ7umr1/CktIkFze/JOr8wbmkLXJ0VKEpbPS3+EZJJLTmIT2OaJMojZBUcWXYikgH4eQDfBeCNAH5QRN7YOu3HAFxS1dcC+N8A/Ex47hsBvA/A1wJ4N4B/HV6PkENTKuqVDQT1/JfJr28qvDTakwXCa0d2bVbVnDK67xOLT1DFh148jy9M/tfVvGH/28zfMH/rfjIa4rMKlM6vDF9Yi+1JCWcV7aAt1uqyIL3MdEd+hFwrq4js3gTgSVV9WlVnAD4O4D2tc94D4KPh8W8DeIf4hPx7AHxcVaeq+kUAT4bXI+TQ2DBJY4jpADSjuLhd14ZCVBcaVcKueeG1Hx8Tgzv/pP8HXhU/sL2D797bgwjw9PCHjvlq2u8/f+sSYfs56oCZK3wd1ZYoCoftWQnn4r9DiOqC6IwRLz1jWLMjK2UVsns5gGeS7WfDvs5zVLUEcAXA3Us+l5CliFGdi3W68F+EH9S4B/DNKyJBcojj7QId6cvjrNPVLH6D//min6JF5CZry28LEK1UpqpfjRzw694VDvszCwl/eGQShhoEyWVhHyGr5Jbp7hWRD4jIWRE5e+7cuRt9OeQmJDZFqEtiuyA+iMZhy0F8Go7F6T/qHs2u9OXxi24BIaq77m97wO2gJ/oB/TEC9OvUOefgSoeJdZiVrhKeiZGdwEd52S3z00RuEVbxv6jnADyQbL8i7Os8R0RyAKcBXFjWflexAAAgAElEQVTyuQAAVf2wqp5R1TP33nvvCi6brBsaIojKcvARkFbi0qolXqSWHKrnVE9bOn15vRwYo7rrQTVk44Dboppe9Vrh9Rx881BcVmlWOOyXDmWpMECVwswygzwzyNiJSVbMKmT3GQCvE5EHRWQI33DySOucRwC8Pzz+PgCfUt8u9wiA94VuzQcBvA7An67gmshtiMAPRPYRns+NxfqcINTpBFX3pcTTwvGY9WxEdUCv0a5ParPjfXsaZW4oXSJE/R3FNLFfgdzCOsWsVMzUD/4XAbLMIMsMBhnrdWT1HFl2oQb34wAeBfCXAH5TVR8XkQ+KyPeE034ZwN0i8iSA/x7Aw+G5jwP4TQCfA/DvAfx3qmrb70HIMih8zQ4I3lKENJnf0W5SqVKarXFiSzelhIMLV/peEf/P5mYjknrN7GPH8j6v3ljhi8XvW/3wAx9Z+w5YW1g4VZROMXN+lYo8iC43TGGS1bOSQeWq+vsAfr+1758njycAvr/nuf8CwL9YxXWQ2xsjvpOvigmkruO1IzzpiPDas/cvityua2Algp+4+5X48Yv//OBzV0DnZzvoA3cNEA//EcS/qsVHdlDAeAmW1qF0foiCsw65YQqTHA+cQYWsDVVXX5gtvxIb4MO8SGxaCQ+qrk0giNG3zPdxI9KXBw02XxXXnB7teJ4ifI/iOzGRdL46B7jSwtncd9AiLM1UOs6cQo4F5gvI2jDIBBrHbiHpxhSBgVQt+2kKM9bx4vlVyjM8t5NjEp26wcLj+anPHs8bHyP1D0z9pTkoxPi/P1T9avK+sUhQqsNkVt6AKyXrDmVH1gYjggyoWtnFhHugDtySFGb8v5jWFGnEeJ30eW4VdbvpV/5+b2Ql4ufNvFEs0ZzZeTNSNwHVU4AJ1IUZbtRHchrepVSHq7MS04Kle7JaKDuyNoyHBsZINRVVXBNNErHVItNG9BcjvJjG7BXeMdbwDlq94LhTma8a+/uGsJYchtB3S7PHAkDFCzAOMxAA6lz1b+AcMLNeeDMKj6wQyo6sDcYY5EFYxgAwMaJA1aRSR3JSiS2N8iQRXpuFWc0VpTZVF/+/5HGmMqPclhs5vhzZAOEPEKkiZxOGIQjq+TBFAGt9V2ZpLWalxc60oPDIyqDsyNowGAwwzDMf0SFOQYU6WmtFed3Sq2t3c8I7oIb32pNH/wzTv/n+G5LKfOUI1yS3g4K7PA8zooip0pgigjykl7MwJsQAYao3X8MrrMOstNhlSpOsCMqOrBWj8aCu1YXVrg3iYqBRZDJXs4t1u0p8pim769WrUm4/hOLSW/qFd526Mtt0pTYPSm8ajWP70/Rx8odEBuS5wSA3MPDDEGaFhYsrzTvFpLC4OmOER44OZUfWio1hhtxI1RiBJH1WLxiqjYaJdqQnjUgwvPCStltJdPfCexceP3ApoEPyytH8voU1uyUxAmQZkBsDYwwy4weMx1lScjG+MxMIq8u7KrorS4vCOlhrMS0sdltr4RFyWCg7slZsDQbIK2lJFdWlHZdG6iV90k7BtG6XRneHjeq+agXCU9s9lYkIMLjz08dSu0sFd5SaXZXCHNXNKHk10bMgywR5bjDMvfSyzKBUoLSKWZCcah3dFVYxKx22J4zwyLVD2ZG1YmtoMMjC+mghnVkJL9SJ0qguTV+mNaW28Nr0NaWsqlll+sL3XJfa3QOjaxNcX1ozvo5RX68zRpCLgRFTyS7WRrNMkGcShiD4MXfOKaz1EV5ZWljr4JxDUVrMSoddpjTJNULZkbViczzAODPIgFp2Ia1pJKnlybz4qiaVVirTSLfwOgmy+KoTR/sc5fZDvdEd4Gt3R01nPhCbUpYQ3GHTmgaAyYFBlmFgMuRikJkMeeYjutHAR3VGDKwDCqeYhUiutPX0NWXYLqxPczKlSa4Vyo6sFYPBABujgV/pGsmYrmQKsXrEXbNBJd6yUGMyxndzLtuZ2d51VOEdFN0N7vz0NQvvgSEOlFyn3JZkMPKR3TAILjcGuREMB1mQXYZBWLPOv4+P7mLdrigtSuvgnIaIzwXp+XXwdpjSJIeEsiNrx+nNEQYhRRbH2VXRXUcqM/2/uFJ2GuEZ00xn9v7mdxx4zRGEd2Bn5jHU7+bqdod9HpopTGPg05hGqkme88zX8JzWww3KsNSPg5dbfGvr4tg7V8lvVlpMSsuUJjkUlB1ZO+7aGiDP6ppdtTBoOoWYdEd0jbpS1cXZqt8tEdWlB16zde2fZfrCexenMw9Zv3tgGKK6Foep23XV69KdeQ7kA2CYZRhlOQZZhqHJMMgNxgOD8SBDlizjo86vaO4UsNbPoFKUoTPTpeLTID2/f1Jw4DlZHsqOrB0nN4bYyEKre1aLzgutWaNL63JRcJUAq5pfEt31/X/MAQI8mvD605nA8vW7AyW3gM5OzY7nZACGIz9zyijLMAy3QajTbQwzDPOs8z1c8iFVQ/pSFUXpo7rYrOKcH3Belj6tuUvhkSWg7MjaMRgMcOLkuOrKjB2ZWRx/Z+YnJ86M6YzyusbdtdcWXTatea3CWzaduUh4bdEdSnKHqNkNQlSXZwajLKsWY41NKXlmqiEF1vnUZXxcupDOdL4ZRcMHjvdO/TCEmM50qpgWdUqTM62QRXA9O7KW3D0e47zZw8QYZMbV0pOwwGtXJBduMAZwDjbU8DxhRTZF9Seii02Dh0hrPrjp7794yIlQ4kDzwZ2f7uwMjcKz+69qTCjdF80tEpfOPVhMPC1TYDAE8gwY5TnGeY5xlmM0zLExzDAaZHN/PLRxTqEisKGel15EHKMH+NlWVNX/UWIE+zO/8rlgiOGgO3IktzeM7Mhacu/JHOM8NEPEVKagOf6uQ3R5EuHNNask6cx4O0yzSrorSu8wLFO/G9//u9X2YaO56tA1pjVHIy+74TDH5mDgZTcYYGOYYWOUN+p0C18zbiehrKtSmhZFaUNkhyq6s6HOxxoe6YOyI2vJeDzGya0NDMPsHbFud5Ds0npdlkrPmLlU5sIB5137Wgce3ARevQm8ut9fcxxYvzMz3PWyT3SK7sBo7gDJddbtAgMAwyGQDwTjQY6NENltjf3NGFM1/MxfdN8FJZvaPBTFV8Y6niqmpR+Htz2ZcQFYMgdlR9aWe0+OMcrFj/PKWrJLuy1b4+3S6C4P0Z2v6ZmG5HoHnB8m/RfOffVGfVvEQfU7CFCc+jSu3vOJ6j0OjOaWkdyCup0BMN70Ud3GcIDN4QBbgwG2NgbYGg8aEV2v8DDvvL5Laogv1PhmpYWzDtY51vBIJ6zZkbXlnhNDnBgNsTebYBqiuTwTZKUXnjWCTH3DRBRbGQpxMaJT5/yPtXNw7fodEFZG9/W7sAbpHL1RVU+qMy6i2vsal9+LS6cfA7Ki+4NH4QHYPNc9qfRBkd6y5wiA8QgYjYHxMMfWcIATQy+5k+NBHdElJktXhE9vC69FtTorfewUyPwqQVV0N8wMZuKwrQVOKTAasoZHGNmRNWY4HOL0ySE2BmHS4Xz5VGaM7qooL6Q109lVum7LcuhaX2v/5vm/v1hEQXizE/MDzpdJaS47qHwAYDwGRgODreEQW6MhTo4GODEeIs+ySmJxvGJ1eSGyPvAa+o63DsahCgirJ8xKP5/mzpQpTeKh7Mhac//WBsZ5hsyIn2E/RneNeTPnhyAI5oWXpSlNY8KMLK3bEvWnvn2HEeDo6kMYbL/lQOFN7n50/qWOkNJMGQLYOgUMR4IToxFObYxwajTEyfEYg6wZ0aUD9BuXWA3vkIMFmI7DS/bVQxP8Phs+SJmkNPenFN7tDtOYZK05fXIDp06McXVWorCKLPPDEPzUVQqXyC6iyXZMbebGAOogIaUpxjR+aIH5tGaY7nGOw6Q6F9Wtts6/F1fhI7i+PKDmlxvP6WOZlGXKEMDmSd+Becd4jBPjEU6OhtgcDavIOdIlvD6nSeu+93rrTCbiwyq9qerTms6no6eFg6qX3caIP3m3K/yXJ2vPS09vYGdngr3CYWAMiswhsz7Csy5EdxrqQfCCS2fzyI2BOj+mS8J+qw5O43lhwJ1Uj/xrqZ/+ao7DRHoHCPAg4Ul5x4HR3GEkBwTRnQA2NgSnx2OcGo1wamOM0TD369YZCbW6+YgubVCZi/IOeF/tOachPlWISFXLs06RGaC0FrszharD5rhj8CFZeyg7svbcdWKMra0RNguLmXUYOgNrFdYJSuMHKjtX/4xaoNGsAgCZGIjxqbE4wE4VsJUUnTdd0qlpnZ8UuZo/EstHen0HuuTUKzyXQWWK7a96GFLegdH5d2G4+9DC1zqIIYCtk8DWpsHp8Rinx2NsjkfIc1OnisPA71pqYbtLeGhesjSed7hrA9KGFW/AeB8nnd4rHXR/hq0NCu92g7Ija48xBi+7YxO7VyfYLwyKTMOsKgZ5piitj+7SsExDrc7F2fiDreI+depTmyGlaQUAXJXCBACI/+GOKU3ncKT63aJIb/PcezHbfxUmdz8KzS9D7CbUTIB8358zuIzJfb8DABjuPnRNohsBOHESOLGV4c6NDZwYjzHM82oM4zDURE0VzdU1O9MnvGC3uY7MA0TX152Z4tR32TpVGPh/g9I5TASQaYlNpjRvK/ivTW4L7j45xqkTG9grHabWYTgwYT5GPw7PuWaBzbTCCg2Si9FejPQA/6MqVQTnOzxEAFhAMv8DH3+PnZ0fIN3msEMV4v7h7kNV5Lb9qg8BpjUnmSkwvedRDHYewmEZAzh9J3Dn1ggnx2OMhkPfuGMEWViQNQtRcltyBwkvld5S4xA6iOlLDfU6/9i/llPfOGRVkavAOsWk9P8Qm+PB4d+M3JJQduS2wBiD+09vYHdvgv3CorSKWeaQW9+oYoxAM/jozs3LLm1YiVFevDeqgAPUBPGJb3/P4KO5KtqT2EhRR3vHFemljSmNYz37+xjCN93cfy9wcmsLm+NxWIjVT+qc58ZPuB3mrYyL5jY6MYHFwmuNuwOWc14UXPzM8fy+hhUJac04zdjMOci0wHiYL+wCJesBZUduG+4+Oca5q5u4OrOYWsXQ1rU76wz8z2QYVK5Jdya02a2ZbFdpzXAszuEoAtiQ1vSRBiCulp2EiOO4Ij0p74AOusQm2Hmdr+ENz78LwwVR3hDA3XcAd5wcY3M8xiDLQhdrMvYw3rK6AzOVnN9uzppSCS+9qjjsINTbnPohBPH7AvwfFY1/B62HKhg1yEz9miarU5xpw4pTP2G1T08DM6eQwlaTVJP1hbIjtw3GGLz81Aa2dyeYhIVAC6sYOL90jP9RNRAJkYCq775UQKC1yFrbVn1niqCWHRwgma8XOectFxtXVH2Pi7U+vbaqml56YHT+Xb5GZ4rmMfEn6OAyprGG1xJeDmA8AB542SaGeV5JrhqLKII895HdIDMYhls7qvMrSvjXTB8DdXozig0uisqvSqEhtls0XKLxsUP0pgDUAfGPFqcGuQFEFJnJkoYVhVF/XziFlA7D3FB4awxlR24rTm4OcfepMWZFiaIwfmHQTDDMmj+szYYVaYQhCoE2fhP9cQlRXprmdKpwIrDqoKKV7FyM7mJzp/ruTU3Fh2sfkxdrd9N7Hg2pS6lEV3/IArN7Hq1kZwBsjYCX3reBQZZVkksH18c6XZx6rR5q0Izq6nRlc2A5pJYbTCLDjlubZTXUbl6pVgoqbRgS4ZuT4h8lRhSFdRAAAwpvbaHsyG2FMQYvPbWFvf0yNKv4iYR9a7ofFCcA1IQ6kvN7NPnBno/sAECrSFBCpFfJLvyoOqMhagwNLVItnVelNtXUQxXS+TYbEc6SQxKGuw9VzSg7r3u48/vQ/DI2xsAdJ4DN8QjDPK+klrVElwXR+Tpd3ZCSNqU0m1PivUChsM4vsRSj3Or7jF+6LJbdtdCeR1NUUVhA1UHUd5Bap4ABSlWIdRj0rKRObm0oO3LbcWJziHtOj7E3KzArwyrZVlHmBijD8IG0Zqf+51eqqK2Wn+mRX9xWAE5d1QJfneN8pBdTmk4B42rRxehPe/b560I9S8sBkV5/DQ+49IqfwWb53bhL3lxJrV7pIZ0uzXdexinX0vRl/C7Se2lJrtGhiWZNz7Sk1zlA/sCpxKSxXTWvxE7N+iDiXAAlnJcgBDZcg3HamAGGrAeUHbktuXtzjJ2TBSazEqUzKK3/C19gIK5ORYYwD1BNfoPFp+A0Cs3LLwpNFLBQ343Zkp+/ORjJqm0NEWEUaBrZqWnui7f2vjlilBSiKLP7Lmyf/h2oKebOK+QSvjz4OEaa46Xy5kpy6eK1eWaq1OUw1OuyrE5f1rW5+MI9kmt0YbZSl4n0VqmaSoNpejM8itHezAJOLSQ3KK1AYKrPRtYDyo7clozHA9x7YgPTaYmJ28fM+uhO4Ws3mvnuTBH1Y/CMSwaLByFWr+bH3MVBzHDo3DZVpOfHegFoRnomNlr482KkF4m1p7nuzWRboim0XnfPGOCU/F1sTgzOj/89Crk0ZxMnBb4kv4dXyDfOr9YeanKDKDkjGITZUrrSl/Fz9421awuvnb7smhC6s4a3KNJLnpMGfWl3ppG6scU/RzC1LkworRgKOzTXCcqO3LbcsTXEzmSMvVmJIqQzfe1MoYhdmH76MJeFWpMKpOqu9KnNaMEY2WXG1FFfyFOmTSvWASaLzRPaGMoQz0vlB4ThDq55/anoYgQXaa/kkIlgS96MlxffiLPDf9b5fUxxEX+oH8DY3o3XZ9+LV2XfVDWjDPJacsMwvi7W2+pZSvwFZFW0t5zw4naM9I5zKZY0yvP/1pJEef4zxHXxFILxgPW7dYGyI7ctWZbhrhNjTGYFitKidA4utKEbZD7KcwrjfKRlxU8TFnKWPmUJQEMkolUk50Kzoe/ClNDiDqB+DlBFcNY5ZFktOoUiS2pMGoSoLQukdal0sHa66rok0oviG+EuTHGx93uZ4AL+q/0V5Ebw6sG3VEsiRcnFKcHi+8XP5K+jfb9YeI1oDklk2uKwAVarhNe5vyG5KGwR2PC9ltZhCnBIwppA2ZHbmlObQ+xOhpiEJYDKkF4sRDFQP6eiAeBCd6YTAGIAq5Aw93MWl/sRgVPXiOxE6pXQ42BzgWkMQjdBlPU+1xraUIswUkWeUTrhP3H4diq4tvAexHvwBfw6HGa934vFDJ8rfxuv3fxWP54udmGGJpV4FU6bacyuMXZzkkNzfzt92Z69pvp86XY7zbmkjBpNrdF2Wh+M0Z2qT1NbVcxKV6Vvya0LZUdue+49OcbVqcWstJg6hYiFFH6Edyl+0LEKUIZhCn5uMC8+UQGsgwZhSVLPqyI746M0EyI0m/zkWnX18kJJTS9uAyGFqTGSah6rIyKpIr12s0dbeC+Xb0SuBk/idzFZEOHt6YWFogO8rNpdmKYlvnQf4Mcwtut2Iv6Pij6dHCWyanRqtppUIvUfHf5wPQ7SD0swqlCroeGG0rsVoezIbc9gMMBLTm36geZOQ1cmEJKNYbS3H2JgUEcg1oW0Jvx6d40x6In8nIYf3ER+AKrILwZtVf2uFda1ozp/XXVOM45jgzYjI9PxOG6/3HwjHpBvghHBp4r/Efu40Pnd/O7lf4Y3nXgf3jB8ayW6qgEFqeDQul+ubhe/z7Ru16a7OaXzcpeiSv+m4gsPmzOs+HqtgYR18QSlU2QCdmreglB2hAA4vTXE7myIaWlx0rqq5R6wML5EByM+zSkOgKqf6ssIULoqyhMXIrkkrRnrev6H1ITIoU5tptuKujsQQKf80lodkAx5EMzJLt53SS+e+8b8+/Bn5UdgO9KaV915/OedD2OQG3z11rf6cYFpfQ3+u+lOX9bXd6Dw0Lz2lGUFeC3UY/BS8bXm06yiO3+NVtXPvELh3VJQdoQE7ju5iav7FkVpofA/2k4VhYmTHMeGFf/YhShQM4GoQF1IY1ptyC1W49TnxQDEaAyNbQCV/CLtbS/DOlUJhBSgzgsu3jejrGZ60xjBq/NvRpYJ/mL2W9jT+Qiv1Bn++MrH8NWb39pshGmlK9vpy4V1O0FSp6vTmF10R3btml37eM+LoTUsIT5uDfHIBJUA6y5Nvx/wtTx1YB3vFoKyIySQ5znuP72JWWFRuClEckCA/ZnDLNZxrAvjswQ2SWc5VTgjsNb5SMH5weZpPU/j4DfUHZYStxFmWlkQs7Qjv1SQxjQbP+J9OjwgfRxFJ0Hkrx18K75m6634tYs/1PneO/Y8/tWz/wAns7vxTad/CG84+baG4NppS7OE8ExHGrNNKvVqX8eJcQ7MduNP5/eoyR8AKtUMMALn/z3Fb2mSfo4NKzG6i+8BCu+WgbIjJOHU1gh3zQpMixJ7hcXGIIfAIiv9j/Gs9C3pAGBjfU/8UkFV96VVqJFGPQ8icFCIq6M8RRiwHjBJ/Q7oTmE6bf64KhS5NKfrSptUukQnYeycCRMx++m/fLfhCXMPdt35nm9HsWPP4w8v/SLECL725Nsa75uOr4v3XcJLJSchRdwltUW4MBwkNuqkM7Wk13AQldDglxRSB8Ap1CiMMYACg7y1RFA78qbwbgmOc/wmIbck95/exOmTG35cWRDBxjDD1jDDeGAwGmQYDfzq3MPML146CJMjx8d56GDMjCDLDUxc4DSdezLdNqYxPm5uUHichLk9u4nUU3ulA8jb81v6xVUNstxP85VngkGY+msYOi4HRvCWUz+IXIYLv59Sp/ijSx8DkAoOjfv2wPIovDR1aZDWELvfy7TEYsNsM9WA/cAqOiTbM9OoKkp1KKzDrLSw1sE6N/c8P1FAu4mI3GwwsiOkhTEGrzi9gdm0wM6kwDgXTEtgNMgBWGTiQkTiYI1CrG/YdE6ROaAofWTnYs1OfZ3NmRDZqZ+ayjlFMj1xldKMzNfvUHVyCkJNSerEZzuia9TSQuQTF1v1YqxveYj03nDybTBG8MdXPoYdewGdM0zDpzVT0VUyk/4GFQmp0yp1aZq1uzYhm+jrY8lltKOozucukF97YHm1v+ucJJVZapjFprDVXKERp75xiU0rNy+UHSEdjMdDvOyuE/jr8zuYzEpY42BVMcjDem4WyI2gsApjgNIiNKzAT/oMhDSmg6pUM6U4kSpqiNryAxwU4gSNsXVojukyrUimPdYOCOk7ad5Lkq6MkmusMp7MiAIAbzjxVnzN1lsBAL/y/D/Fju1Kayp+4a8/gLfd9cP4utPf1tmY0lWvE6BKn9b1w/lXr2aTcS2xHfxPd2CNr9eDjTF4zfF4groOaNVP3G2tVtfPLs2bH8qOkB7uODHG1WmBc9v7UIR10ABAgVGew4hDJg6ZEczE+ZlSnEDgEtmF5hOnMGLgqsmeBWVcNVzrVbbT9vfKbRKjC1OdoaKAmnrmlBjNmfo+FU4UmqnkVu9rdHVKPQ2aCPDNd/ww/uDiL6DU6dz3s12ew/977t9AjOBvBeHNpS1lXmzp4xjlRVxoNAH8HxNtukR1LG5JU5rVGLz6UD08QatJpTMTRbiatCpZLZQdIQt46Z1b2C8dsAdASx8BqY/ynBGYQYbM+Q7N2LASOzedi2PovPhsrO2owMF3dnr8oGWXZDFT2UncjmeL/+VVJFFdh2CieLIkmotpRJPMBBLFFn+8gTo6ef2Jt8KI4D9f+nXs2HNz30+pU/yn8/8n/vYd39bbiTmXvmw1qABNyQH9Auseh7dk5NZF3ySa6SlJlFePj5xvVHG2HmyfZ5TdzQZlR8gCjDF41V2beMo6qAKz0voVEBwA8VOIwRiMBpJ0Z4bZVcJAbw3iq2Xn17szqIcdGKPNzsCkD0JEqq7NWg7SkFSj2zGKJdynokvrY+nzEa6kcTyI6w0n34o3nnor/uXT34uuGt6V8hx++q/ei9P5PXjnff8IfycRXyOyM7WE42MIqvlIG997Xw2vvW/J83ppPb9r0dfuKK9bkk4VM+v/LQc5+/9uJig7Qg5gMBjg5ac28BwATAGIAiWQqUDVQULkZkRgMmlEecaIXycviC6VXZwjU+GbMFxcUQGANfUvrF/xIBmBJwrnWoO4UQslRnPSEcUBXZJDkHI6Dq4+Hh+fyu/Bdjkf3cWrvFKewyPP/ysYAR6689u7hZc8BppSjyzTmZle9/y+xZHetcRc7ZlW0lRmF4VzEAvkGYV3s0DZEbIEJ0+McR8EL27v+bSiK2E1LMSqgEBQAtU0YTCCzACZc8iNT9FZ5/w0m+rH5JVSj4IuwtguhNfIksYMl6TRqpSm1JFb1YiCWnZpLS7SJ7lq2rJku+vxt939D/H7L/7rzvpdpNApPvnCr+Lr73579fws3ifX6+uC3VFZX1TXJZb2uUeO9BK6ZlpJu0Lbqcw2hfX/cpkIm1ZuAig7QpbkzhMjFLas6mmz0gIimJV+ALJTv+AnMvjlf8Jf/pkRZPDpxDgY2i8FVM/8EdfMizOrpMO5XHi/KDQbaoEiTdkBHVFN+E9bFpLsi9XAtKW/K7L7utNvgxHBfzj/a7hSnkffsITLxTl86HM/iu962Y/izN3f3ojsYqNH36DvvnrbMhFcH9e6HFAfVZR3QL3Pd+AqYMCpxW4CKDtCDsFLTm/BqTeRCGBKQI3zDSkQvwo5BEWc+suEmp3W9TB/X6clY4ozqiOuUo5qu166B/BiBPp/tGME13VOGq3FbSTpy0WRnQnC+1t3vA0igp974sdwpSetebk4h9/+8v8BI8Cb73l7NRTCz5jWfd1xkHnX5+na3/Uy7X0HDkPovBJPYzxea2xe+CgHEldPcADUKtOaN5AjffMicpeIfFJEngj3d/ac9/5wzhMi8v5k/38Ukc+LyJ+F20uOcj2EXA/uPbmJE5tDbI4HGA5yjPIMuREMMz8LSSZ1B2RmalHUY9piOtG/nkhzvbl0VpVq9pSYlpT6fKAWm5G6/T3W6JqNJmiMpxNpTrElreOxkSR9XD/Hf4Z33vePMJBR7/dU6BS///xH6+hxgeiA/mN9UV2XyI7UmbkkjfipMQ4AACAASURBVJlW9GDlxdlVFEBp3VLPIavnqJHdwwD+UFU/JCIPh+2fSE8QkbsA/BSAM/D/3o+JyCOqeimc8sOqevaI10HIdSPLMtx/cgMvIMyaYgQyK2EsgDDOCvBzLYqaMEYupCOjrFR8o0vyY1wPKG/S+MHWZoqxi7rJZP6Hv53OVO1PX0pLeHF/PP+hO78NRoBPvvCruFx0R3iXZuegoV1/URqv71BvDW9JibVPm2tWWfJ1FN0pyzC5ygHPTRtcfPdpbo6eTiWH46gx9XsAfDQ8/iiA93ac8y4An1TVi0FwnwTw7iO+LyE3lDzPcO/WGCc2BxgOcwwHOQZ5hkEuVZSXG3+rZyox1Tp5XVFeO8JLx63V263OSmlGbtXg8URYadTWrMOhI9KT1vnz++vnC/7uXd+Oh9/4K7hjcG/PN6X4qT9/P85e+FTvd9mXvozv073/4LRmZ6TXdYWhE9Y6Xz8trVt4K6xFaf0kAqXzzz0oWmvPnVku8RyyWo4qu/tU9W/C468AuK/jnJcDeCbZfjbsi3wkpDD/F1nwp46IfEBEzorI2XPn+tqfCbl+DIc57t3cwOZogM1RjtEgwzjPwsTKppLdIAxHiF15qeSMSdOXh7+lEgTm5dcWn0nkBTTPu5ZUZrz/rpf9aG9K89LsRXzsS/87/vT8vPAEC9KX6OnA7Di/S2x9RLGV1oXJnbWaWHoZASnikkwKGyYYKJ2iCDK0HQKM0V0KhXd9OTCNKSJ/AOD+jkM/mW6oqorIYf/lflhVnxORkwD+LwA/AuBXu05U1Q8D+DAAnDlzhv8LITcFo1GOewW4IH7IgTECRQnjJ8lE6bwoZkDV2eBcODeZDgxYPo2Zbkt1Tn+qrzOdKQDcfLrSn784ldneZ0TwDXd/O4wAv//8R3Fp9uLctczcFL/3zEfwpnve3ti/KH3Zl/ZcNq1pku+1b1GCxve2bL1vwa9PQ2oKCLT6g8CpzA1VYErz+nFgZKeq71TVr+u4/R6AF0TkpQAQ7uf/Vw48B+CBZPsVYR9UNd7vAPgYgDcd7eMQcv0ZDXPctTHCyY0BhnmG8cCnNeMSQQPjl9BppDXFL7MzF+WZVhqzZ9mfdlpzmXRmuj/WmtqR3EGpzHaEV42fE8Gb73k7fvpv/yr6ehwvzl7ET372R6oIL5P+H/ne9OUhBKghJZnOZpOyCsEcFJjFCLB0Wi0T1IbLA10fjprGfARA7K58P4Df6zjnUQDfKSJ3hm7N7wTwqIjkInIPAIjIAMB3A/iLI14PITeE8WiAO8ZDnNgYYGM8wHCQYTxMhJfVKc0skZ4JkssPSGealsgWpTP7xNeu39U1wcWS6xNelQJFIk0B7hz21e+88H79iz+Hz5z/VK9s+tKXsZ7ZJhVgOlONa5noILW1j6861lL4yQOKRMBxf9k1lQxZKUeV3YcAfIeIPAHgnWEbInJGRP4tAKjqRQA/DeAz4fbBsG8EL70/B/Bn8NHeLx3xegi5YYwGGU6Ph9ga5dgY+YaV8TDHIDcYZRnyTDBManl5j7hSuVXiMsn6b+n+DglWkZ9pyizuV2hzn2mtgoCex2k0lwivnuTZN5i+94F/jKHpH5Iwc1N84pmPdB5bnL7sfj0fqWojimuf21XTmz+n95IXcpi6WxRwjPbi9VYD0MmxIbdigfTMmTN69ixHK5CbD1XFZFZityhxdb9AaR2mhe/emxUlnPp5E53z04AVYdJgm9SV/Erc/nGcPcUvAdScWQXoaKVHf/0uHos/qum+uF2lPeNxkUYTSZXGTB4b8U9IZ0b50/Ofwu898xFc7KjfpVd71/Be/L1X/mO85d53AOhPbcaZaLo+FyAN4XSd27WvPcA7bdxpH28/Pz2Wm+WbY/rOj99jbgynFmshIo+p6pmjvg5nUCFkhYgIRoPMz5u4KdibloAIstJCRFBaC2MFhYRFXdXPn2jUN7Oo+lXOJZGaAl4mWCIVl/ZboCkOgV8HL/3RTqUWn98lvDihNJBILuwT1EvfRN50z9vxpnvejp/87I8sEJ7i4uxF/NrTPwcA+OaXvONQogszcaHdMdLliq5Ir+t9ko2eaz46XXNqxiWOVB2GYtiwcgxw7hpCVowxBuNBhkG43xhlGA5zbAwzDPMMw0FWNawMc5/WHLaGKnStJn6YoQhdNTkvpjoqq5b+kVYaM6ZHjW+kyWKqUnxU0t6XyrLNew5IaQJ1WnPZZpWYsuxqi7zmmVXaxxcfPhKLsmlWFdPCckjCMcDIjpBjIMsMNod5iMoERhwK42fQ8B15OaxzMFZRwlUrlbswj2JcEMEBVerSqS4M7TqjlfCfKspz6XYr8ktSlGmE1xXVRWFCupfpibzpnrfDAPjdZz6Ci7Nz6OvbvzDtHjvbXp6o7qzsaVbp+H66RNV+7qHqd0cUXzqjShcOwKx0GGRMaa4Syo6QYyLLDLaGOTAtAQTBQDArLUQsrDMojU9rWucgElKZEEhoWnAA/Cxk/T+OkXh0kfRMmrKs9jXTlp3C69nftfBq+p5GgLe85B14y0t8Te4nHvuHnWnNu0fzHZyp0GKXZaS7ftddO7uWZYAO2j4qB00zZlUhzsGocPLoFcFvkZBjJEZ4sSlhkBuMBhk2RoPwOEeeGQzCzc++4sfgDcyC7k2Zv9Xj9HqGKSSpx5jCzCVJTyYpyvg4TV2mqU4Ac639KXE4Qlssf++V82nNoRnhe1/5Txr70jpdW3T9Y+069vU2tTR3dKU++7a71uE7LMutmMDJo1cJIztCjpk8N9jUHCIW08J6WRhgNAipTBE/9spaFDa2ZFqoq6em0tilEn8mF/32SeOu+qGusqBpWhOYWxOvGkbQEcmldAV1MZrri0Jj1+UnnvkILkzP4e7RvfjeV/4TfGOI/CLRUa41Xq4vfdlXq+uK6ub2zb/Y/IdaMQethQfUywNBBKXT6o8Vcm1QdoRcB4aDDCKoVkCYFRZGAA3RHiS0sxelF5sqrGg1Q4drOy5dZiZ5n8ZPYSI9Eb+iOq5BeF10RXVpbW8R3/ySd+Bb7ntn7/EYEaYDr4FFww8WDUs4eN9cva79HLSFOf/+h2XZOM06RR5aN636gi6Fd21QdoRcJwZ5hk0FVEsgR5hJw9bRhtTRXmYERWn9cAOrYX5HnZNd34+mtB5UVbkouXqzjuCwWHApaVR3UDR3mPP6RAcsGlTe8V69DSxzpjr0YPOuy6hmQ1nij5D4esvUYdvNLBTetUPZEXIdGQ4yQARXp4VfvTw3KKwfo+bCzB8qApPV9SojSc0qCNAfbNzN0RYeXP2j2SW8ZYlR3aIUZ9e1LFrPDlgsuq76H9DflNId1XXJb/6Jc+nQ5LFfGUGqVRKuKcVZ/bs5ZOH9Fv2hkUZ3AIV3rVB2hFxnhrmBcxkAxdQCeZi+a1bGHv44w4n/MUub8YZ5BudcJbhl+xbib+hhxdaFnyllOckBB5+bLhfUtTZce/hB+ry+tGbX+3Wde3D9zn/eVG6SplaPUNBTBVTqqND/o2q1YG2UX9dQBQrv8FB2hNwAxsM8/LhZTNVCVTDITJWqdE4hcX1sEVQCRJ12rHYsgZ8P82g/jP69tT+f2HH+QdFcKqx2MwqAatB7F30v3Sm1Hin2vWa8FgnNRNVzDqjfrQJFkFlYIsiIwMn88kBW1S92fxwXsYZQdoTcIMajHJA6LViKgysVuRFYAEZMEu15Eu9dF9qpSrvEey+b3mwPLziM6PrSmn1SWybSi4P6NZmA9Frqd8uiS0wAF8VnS4tRns3VPLke3vJQdoTcQMbDHM75hV5RAsjCQO3QuFL9IGcGZai5WWcP/T5LpzvjfYesFLrwdQ5Tw2uLrD3j/8KIbkFac9mFXdPnx8VdRboWzl2c5ryekimtC1O+abVaBTBf0yPdUHaE3GA2x0PItMSelmGPm29cMYIsHNXMVMvCNLr/DjnwOBXbMrW8vpc/VKNKUp+LtNdy6xPdH73wB/itL/0yLkzP4Z7RS/ADr/lv8C33v7N63cPW76Lk4veWGTP33PYHbU4WvfizHsRhx4k7VRhIFe05Gyb1Dk09B6WMb3coO0JuAjZG/v8V91CGdJ7xKyI0UpnNxhWp/hOpU4JtxCH8TF57NNIeRH4YyQHdEnPJi3aJMPJHL/wB/t0TP4uZmwIAzk9fwC99/l8CAL7l/nceun7X7vjsSn/Oj79rR3nXVy7tf1WFzwIIvOgEnEtzEZwujJCbhI1Rjq1hjmGe+Ugu8+m6PDNhei9TCWERsR2/73YUJNzilGLLRnNd56Z1unhO3/X91pd+uRJdZOam+I2n/+2h6ncA5js+O6K/ru/qUJNFHxOuY9qaKL0ZpxVbCCM7Qv7/9s411raquuP/Mdfa55zro5U3VKiApRqaNDS9bUxbTcvDS5sqJiWpkehNsRroN5s2Skj7wbQR6wfSxARCKAppG6k2UWKaS/Cq1aapLVoUJNJ7RRtA4F5QWqCXc85ec/TDnHPvuddej7nWmmvt1/glh7sf63n2Yf73GHPM8V8iJnN4YIwzArOxe2EijBLjMKBIgZ1TwgCDG5FJ+ZUZq1btV7WIfGoiW5yCdChFpa4Iz++eCJ6/yzQXXk9hT825ybmi/YZXu6pPW1t7oO1RIgUrBYjYCcKSMZnDw9guKifwWE9E0IiDmvRXzM/fdRHASYo0N49XZeNTdIy6Tilac9B2LvV55vbZeG732bn3z9g5u/D8vni6RtKFkW3RnF5hVBexDLMDmhlJxcnFHqgcSWMKwhJyYDvFq0YpRokCKWPomiYKaaJm19lhGsU4dwO33cTQteKndJ/cgB8ioORdQ0gbrKrt8qnP37/4DwvdEt7zxg8U7ufQeip0xYvK589dVNg473+3OCGp+ywyZvNTYb+0iUhkJwhLyoHtFMwaIMIeAA2NFAqAnlmDVzb4lReyxMNES+GD/yQlWhF1FAmTq7q89/E78fzuCZyxczbe88YP4K3nXTWznb+bX4QS3D+zIKorXnze/Rfa9hB1XniAndtTANuKTUlritgJwlLzqp0t4JV9gM0AuzfOoDRBw7gkmF6NUxugcuFD68XoeecBIKzx8+z5p4U1VQvTq9bX/ca5V+JtP3NV0EJzfzlD2bq8IsEtjupmny9aOEI+xsxbpiALzw0idoKw5LxqZwQiwsv7Y2jNUCMjegAwznia1rT/iTWH5w+OScvBMj8vV5Zaq5u/q3u/TOhK3Q9K1uQVRXUhBS1taLt0IcQLD8DM2jsRPBE7QVgJ3Do8ANjdz5AqBU3O+mc++ppNXdoHViDrXBLyA2KT4hT/GvLixDzf5BmojubcsarSnlVC171RdH7n+IJRZA+Ux52yybnzxSybLngidoKwIkzm8MDYJ2B3rJEmys5NhaQyvTm8Hqibv5tblF4jYiHbNBU6lBWqBEZ6XebqZlPOpiuO5sDj5X2dMP1SQxUCrDXP3NsmC56InSCsEG5ZwktsHM0zAjLWM/50gPNbM/t0XYtXJ56h83fuOCFLDpqmNYOEDsVzcsGRXoNCHACFUXfucK1xgjl5UmANBExbjPlsquCJ2AnCiuFSmi9bP7xtmMHeH1id87nBzeOpmQ4cocsJ/EExv/4uFOb6dKV/zibVmlmu5Vip0DUoVCluHxYm6FUC59NFbIpczousgYiKt93E5tEidoKwghjBY2A/wx5sMYIX1elcl2gi44dGMytrZ+fT/G19ODTVVoBLbRp/vvqUZUjEl19H50eMZULXpCKzeElCtdgVefFV0rPO+MIHMEZJbv4U5gtSmmzOUmsRO0FYUQ5smyrNl9iMnfuZNnM0VjTcWOeEr2r5QcyUVn7uri7KaZvWLOqtWYRShK/96AHc89jteO6VEzhz52y870034PILDjVYklB87HxD6VC6NpEOq8c0jLUGkZqzBmJgo9wSROwEYYVxjucv7RtRU8STAdgVMKiZiK96IXpTytqL+ZRaA9UUtPjbFQ3ILn1Z5333tR89gE8+fAt2bSPpk688i08+cguUIlx+/qG57Yuvc/b1xpHc3DFb79oKV6jirIGc6GlmkC5vmr1ObE4MKwhrys52itdspdixDYBd+y9lvc4ciqYOCkXtxNygXvXjtxmrai/mkxcFv61YyJKDKqGrOoYTwXseu30idI7dbBd3f+/2ufMVpS/zaVMzP9rty8LQ7cY49zhjcx9sW4ttAhLZCcIasLNl/ldmMPbHplhFAXPpTEU86WAytxShZPztMmc3OXRgFBey7dRstVxg/feee+VE4TYnT3mNpak6femaSUdhAUFUUQNp111Fkfm7WPd0pkR2grAm7Gyl+KntLYzSZBoR0TRiSb2Ij4h6zaVNIkHC1I8vIF2pAiK+yibSNN+I+swCZwQAOOvAOXYfIFXzQ6G7B9dMOhaLaiJdlrrWzNgbZ8jadA9YIUTsBGGN2Bol+KntEVJlzF5TpWYGf5MatKlM5RvDegLYYDDOpzjz6c2QMCY0remuq3QbK+h5EXzfm27AdjLrmLCd7ODwm28AUFJ8YiO9GCnLgkMvhKrbYAB7mV5rpwRJYwrCmrE1SvBabOHlvX1k2jReThiYFqnwzMA3315s+sR0/Sgu2gihtDXZRLgCjmWXTQDlTaSrilR+6/xDUIpw9/dux8lTz+KsA+fg8JtvwOXnHyp3OcfsQvVohNxv2GEaUydjbjkCM63lkgQRO0FYQ7ZHZobm5f0xyFZnuk74W2kyiVjq3M6LKhGbMNezs6Jqs+Dkdj7JbFsYdbj1gyXHc1Hj5ecfmqu8LBM6wnxrs1gs3gev+vzub2Sc6bWzBhKxE4Q1ZWuUAABO7Y8xZoCYpqLnBjEb0vlOCWYdcsfR3i17UBQubt6+eQEraiJd15Gl6v2ygbzvJF4ssWsjQqH35tbejW2jgnVZliBiJwhrzNYogSLC/+2PodlUYhITmBh7Yz3TfWQ2lUnTbvxQlZ35Z3QT04GYmcEN5//K0pqaw7ZzF1QV7VUJXZ/D+jJESUWtw/L4lZsZM7Am6/BE7ARhzUlThVdTilP7GcjO2Wn2zV9tP0fwzNd/NygmBExmrxqMebXBobfovWoAnkR1NSLmrrm0hJ5MMUqZmPY9nkfzwRtAd3y3hHURPBE7QdgAkkThAIDdcYaxtt3xobGXmZ6Vrq/mbDrTNhUmxMnvkWdJgyZFLlRvHhs4d1fEEELXde5zBp7aBdmn5eednN9uyxQklnm3hHUQPBE7QdgQkkRhhwh7mcY40yClABAyPXVMmE1neulIuDZj5lhBA6wCFDcTtulBpsUppiqy3dxc0PsDjN9tF2znPfAAcz/ctFB08rlpJEyTz5lKviAw5r3wMmYQL0c6tg0idoKwQShF2Cbznd0UqjB2YQYxtzQhX6E500aryTinGwaEBWlNXVEWWWsZVJG2tG/3XpACNBcHrXlG3OaO1+FamAGmeS885RUUTa5jzbzwROwEYcMgImylxuF8TGzbRunp4OalNN1cHhGVDr5lBG3vik1QVpgyf4w6vzugRgjJrKMLWV5w9IkjuOvR2ybr866/9EZcccHV9Tt6hER1TRpL9yE07txZNjWBLUtZrqoXnoidIGwgRIQ0IRCZwQ3M0DQVOG1Tmm4uT4GRaT2ZzwPsnJ55UHcy84/3tM1yBH/NXdV91Tmcu84odRx94ghufehj2M1eAQCcOPUMbn3oYwAQLHhVc1xNjF4nDKAxrlF0lpn8qfk7odn3V9AaSMROEDaYRBEICkTA3lhDw4iRmjSPNmlNN7eTX57QJ05ww7qsVKcsgWm0F9oZ5a5Hb5sInWM3ewV3PXpbkNgVuiigpchZhl6UPtYaIGW/aHjp5RW0BhKxE4QNRynClkoAkElnalOOwpgOrkb04veJnIMwmT9yi+Drtq8TOT/tWTUHmGfGFSHg9Txl1kRdvASHlhZ3pW65iqLpfa1awYqInSAIAICtVIHGjDEABZrMIznRSxM1EbuiJQqNmNS7lLcQY10egQVFfFY4fcf0JmJ91oFzcOLUM4Wv15GPeLqKHIBofTWb4ldlas/8NVG0UvN369ftUxCE1qSJMq4FNHUwSG36zx/TXIrOuCeo6X6KKn/cdqnnyFBk/FooDGRExJ2nLmWZKjWT9mva0f/6S2/EdrIz89p2soPrL72xcj/KCawzSe0KLcgvoejKNU9NbJtEy4tEIjtBECZMilKUF93BzNklKgHGGTSmc3lF+zuUUq0HebdbE9NXc87iebI2A7Kbl2tUjemn+WJEcx6Lmh4rayDtzF81a7OcZcnTmSJ2giDM4BYtazbpKz+lmSgFeI2k8+nM1mN7Lq0JYG6NV9W+frqyiLZzjVdccHX4UgMyJrBRXc29Y0exBmpxiLo7Meav2qTBl1jwROwEQZgjUQToqUg40UsUY9c2kGbkKzS98nRm4wlXMlJOlyEUD466boglI4whlZpDpdkSu3SjDwPURVoDAfUNpDNm7Gcao2R5BU/EThCEQhJF4IxnZIeIMLKFKq6VlStiQW47oshFAX6lZoMBtfcKUtjfFfd0rgZp3IBDtcItP6lCM5baFkjEThCEUtJEmQIL7zVFZmCbCpqX0vT6OHZqIE3VlZqhDBHVJYpar5sLIWZU12fUpZmhGMgAsG7fD7QvROwEQagkL3iu7D8fxTjxs8/Mf71KxLDu/G5tXZyYsO+orm+hAxafwgTCv7NotpZQ1iVhmQRPxE4QhFrSxPTSnMzhUVhvSTOVN9subCj6juoGEbqIYjGEZvrGr8smeJ2+PhHR6UT0ABEds/+eVrLdESJ6gYi+mHv9IiL6BhEdJ6J7iWiry/UIgtAfiVfWX9YKK0+btFks7egtqqNhhC7mXJ05XPtjNblP/0uG7qlgpw1dcwUfAXCUmS8BcNQ+L+ITAN5b8PrHAdzKzD8H4CcA3t/xegRB6JFEERI7AId8Y1/Ud/reRMguL6gTuqNPHMF191+Dt3/+Lbju/mtw9IkjjU8VS+jYLvx2wjO2fob7FT9um0ybfXSDNYP5LxnLInhdxe4aAHfbx3cDeFfRRsx8FMCL/mtkvvJdDuBzdfsLgrA8mO4kJk5IagbkPtZ1hdDH2EpkurLUDfzOLeHEqWfA4IlbQiPB6xDVcU7U/PSzq5yt+/W4bbRtsZaBMdY8EcKqBfNFXnzL0Gmlq9idw8xP28fPAKhvGjflDAAvMPPYPn8SwOvLNiaiDxLRg0T04MmTJ9tdrSAIUSCatv6qGpIn6/AGJmpkZ9OWZnlBfX/NKreEUOq+ROTReipwhULU8UPwD+dE0Imfi/xmrqfgV5Rx3I4yTaktUCGiLwE4t+Ctm/0nzMxE1NudMPMdAO4AgIMHDy4+JhYEAYkibKfKLDQv2aaN8WsXYp7Ld0wIXTDe1S1B1fT9dDSxCuqzotM3flW2Ute8Uuz6sKjG0bVix8xXlr1HRM8S0XnM/DQRnQfgRINzPw/gdUSU2ujufABPNdhfEIQlQCmF7ZSwlxVbAIVWbsYiitbRvHVQ6LxTF7eEkPRlE1dz77CD4KyACMYUOElmk4cMYJxppMnwHgRdz3gfgMP28WEAXwjdkc3Xka8AuLbN/oIgLA9Kmc4qo0TNDdZDt4/qVIXpnBXUbNurJvNNbd0SXPFLGVpPnQYaEaGvZtPdXZPo/UzP/e4Yw7Vw8+kqdrcAuIqIjgG40j4HER0kojvdRkT0dQCfBXAFET1JRIfsWx8G8MdEdBxmDu9vOl6PIAgLInFFK7aAxRe9Jmm0hSS5PJHLX2tTH7wrLrgaH7rsJpx94FwQCGcfOBcfuuym2obSZfN0rUXOsqhF6ZP0prMD8gRuEfN3tMgJw7YcPHiQH3zwwUVfhiAIOdgWLvjPmYFMa+wHfpvXmpG1tgZqVuYeYgI7zspNZGNRZE0Uq6l0jJQh1USdZYxy53ZfhlzhUsi1EdE3mflg45PnEPNWQRCiYcrzaea5UoRRmiANXoje/vxBGknTwpM6E9gh0m1FQufWt3UlVgq5bU1tPphy6c1M8+DpTGkXJghCVIzgYSbCA2zLMdtOSrvBLnJmqfRoS+qYkO9EE9siaNGdusrcEjQzOGNAEZRKBrkWETtBEKLjHM/9dKRSBG0tg1zPxwQ055bQxS7B7L8ajgn+sgZ3vqjiGqEwZXKoloepuhsGsK8ZNNZI0/6TjCJ2giD0glIE4tkITxHNzccVuiVgGlW5refWSXsDsFcKE+HKDX1GdXmhq+pI0pamC9Or6DNA3NcapKn3htEidoIg9EY+pelHd3X7TWx/kH9QRRzB6DOqy8/R9VIAEzGqA9ofywh49b5u7R0z9br+TgpUBEHoFVe04oa8oCKVFueJGRn1FdX5Qse2JL8PokZ1A1kDTUWvp999L0cVBEHwcL00la3OrBs/F+lX2stgS7MWSbELUWZORWHtxoKP18UaqMF2bAWvj5QuIGlMQRAGJFEE0gArmqvW9KEORSpdia5BuTVqVUJ39IkjuOvR23Dy1LM468A5uP7SG2sXo+eJPfU11BcP53LuBC9RcVOxInaCIAyKUoQRKYA1xhXf4JvKXYxggCN39ihbQ1eEswZyjgnOGghAsOCFNpFuwpB9NRXbxuGI3zRa0piCIAwOEWFrlFTOLS2izVU0nculLR1Vc3SdrYEiO5sDJqrzxdN9GdB6/ifGF4W8lVDMOU2J7ARBWBhbqZp4ouWHyUVkMjsXpljBKRKdujm6rtZAsYpS/HWPBABcbt80u+PsA7fGMSGeE80yNDN8h8SYH79EdoIgLAzn/p0maqZi07234AYgjShrJA2g1t0cKLcACrEG6pq+zJu/ahulEbUXHNchJ+NmLud9IWInCMJCcfNarmIz8dfYDZjKbLu2TqlppWkRoa4Jba2BiuYF63CpyLEVIF2Ugoy8Vg8IdDnvqUpV0piCICycRBHYayWmbBsxArCXhQ1+XdOejXatSFfmCV1i4IpQmlRj5ntr1rEs7uaOOtkSTAAADPZJREFUQpfzns4rYicIwlKQJsosKrbPTaRnBr5Jaq3H89cKgO25WWcJ5NN0Ld0VF1wdvtTAFsGE0KbvZgyxa3II53KudQY1SqJHlSJ2giAsDYmiuWKVRJlSdFe4UOaY0MvavJZuCUD8ZQz56wrxl2vbXHpod3kfDdNejsBR+2WK2AmCsDSYBsmYETxTqDJ9XumY0JKYbgmOvjqkhAhd1w4tsTSmkw8eVTceaIqInSAIS0WR4BW5Jbht844JvvBVDZV+EYzWDIoYhS1S6Dq32+qhMKUpZT54XRCxEwRh6cgLXhO3BFMyP7UKCiG2NPXVX7NK6GL54S1DE+k+virI0gNBEJYSv3k00J9bAhBXnPqI6tx6xKpzxhC6+E2k2yGNoAVB2DgSRaZqQaE2umtboxJraO2jKKVqHV1s94T4TaSXpy2ARHaCICw9iSIkVO9mvejBNXZQV9Rf0xFb6KJHdR0O1UcaUyI7QRBWAueWwKwLi1Uczd0S4g2tsY5FNcIea35uesLwNXvhh1yeqA4QsRMEYYUgImyPEuztZ6X2QK5IZWhitbmqa/9VVG3Z1QcvZlGKo+sh/UKjGIjYCYKwcoxSBWTadNzIDfyKejBgDaBrpFUXzQHlQtfFBy92+tKxCIumKmTOThCElcMIg0KiCGkuElrEvF2XSNKJXBuhAzr64PWQvgSGczdvgkR2giCsJK55NIiQEKDYW5dHFHdOq4Y2kSRReJ/NKougLj54faQvgeWL6gCJ7ARBWGHSRM10QnG+eH0N4mU0iez8SC5E6Oosgtr64IWevw3LJ3UidoIgrDj5NBwRYZQqjHLeeH0RJHQ09b1rIjIhywva+OB1NXutwkV1zjPPGcOW/bhtem2cDUljCoKw4pjuIphrGqyIAFXvltCVssO5vp1dUnoh6+ia+uA19cCrI9+MO1EE3eb4PH1AABQlUfuVitgJgrDyFAleviqz1C0BQJdlzLoHxwSgWduxUB+8kIrPEEpNYCM1ke7jS4mInSAIa0Fe8PLWQPltfbcEt02IW4LZY3qcPqgqSGlLV6ELcTmPFTH28WsVsRMEYW3IC16iwjzRXNeVxm4JPc0xxY5qughdk24ty1iF6RCxEwRhrXCCl2mzLMGkM/sRpT4OG9s1oa4jSxlNW5KpiOv1+igrkmpMQRDWDrcMwRQ6hG3fhthaFzt9WdVIuuoaxplu9gWhYyHO3OEkjSkIghBOmijr9VY9eC+L71q0CJTMgvEmIt7FRSF2+rKPZKiInSAIa02iCGDCblYhdi198GISq5F0m/m5slZkYSfsQex6CO1E7ARBWHuSRGGHCHsl6TlqqXYx9bFzVDdwNOcYultNW0TsBEHYCJQipEzINAr98BYZ3HWN6gaP5ix9dGLpq6JTxE4QhI0hUWZNnQLZuTx/ETpVmsLmiTlf1zaqa9JM2mec6Vbnmz15P8LUV6AoYicIwsZA1iEhs22tEk/0FjVv1yqqa5GyBObTll1MX1fNMUHEThCEjUIpAutpNJUogmJAMyHT2eB61ySqaxvJAcVC19b0ta9G0n1O/8k6O0EQNo5EzS5bdnNeW2kyiFOCo4ljQhNboKLz5AtR2pq+xm4k7dNnBxaJ7ARB2EjSRGGc6ZlITllhUV6zaM3F/TVjUJrBtI2l20ZxPmUVl21MX2M1ki49fm9HlshOEIQNpsgLz2/y7DzonDeeihz1TSI7mp4vUYRUNfO9qzp+2dKCxqavNP/7iglFckwoQ8ROEISNxfTRnB1gy1JpEyGy4pd6qUXlhJDK553ce25bZcXDuKuryXFiDfh1a+gamb4SkKp+5aLvJtKSxhQEYaPJOyUoRcgquq04lFJQbEv4PbugULTutyCjbrF4sOnrAELX1eQ2BBE7QRA2njnBI+rNKcHR5/FD19HVmr4OIHTAMNZAInaCIAiYtQZKFEEHRHddlub1JXWxLIL6LkbxUTQtCAIwZ6YbAxE7QRAEi7EGMgvNk4YdVZrQl+mrEYzux27rgReKEzaGqYDdj9HRpQYRO0EQhBwmolHQ4+pF5mQjkqb0lcGMEdX1JXRamyUcM78vAoiGqZMUsRMEQSggUYStRGEvtxbPp7UPXtuLqqCz0LVsQVYFM0NzeSQ7xFydQ8ROEAShhCRRGMEUrhRbA6GVcsVOY3ZNX8aO5nTJ72uGASowfUTsBEEQKpg4JbC1xfHeWwYfPKBDVBc5mgsSOcvQPngidoIgCBVMnBIApAlBa54pXFm0yXlb1wS3uD3WNTRZStFXI+nKc3bZmYhOJ6IHiOiY/fe0ku2OENELRPTF3OufJqIfENFD9ueyLtcjCILQB8prHK0UYZSoiVC0sdmJSaP1erb3Z6pUFKFjZoxL3N8rr2EB7uZdy2A+AuAoM18C4Kh9XsQnALy35L0/ZebL7M9DHa9HEAShF/JrzhIreumALgl5gtOXkUXOnbtN+nTo9KWjq9hdA+Bu+/huAO8q2oiZjwJ4seO5BEEQFoZJZ84P1GmqkNpG0SHDeMzAri5KdAvDY4ocM+P+H/4T3nPknXj759+C6+6/BkefOBK0b4zm1m3pOmd3DjM/bR8/A6CkXXYlf0lEfw4bGTLzbtFGRPRBAB+0T3eJ6JEW51oHzgTw3KIvYoFs8v1v8r0Dm33/S3PvyWuT00enj94AMsHSMRzDv/K/6P0f7/939mL24x5O+aYYB6kVOyL6EoBzC9662X/CzExETb+z3AQjklsA7gDwYQAfLdqQme+w24CIHmTmgw3PtRZs8r0Dm33/m3zvwGbf/6bfe4zj1IodM19ZcRHPEtF5zPw0EZ0H4ESTk3tR4S4RfQrAnzTZXxAEQRBC6Dpndx+Aw/bxYQBfaLKzFUiQSeK+C8CmpiYFQRCEHukqdrcAuIqIjgG40j4HER0kojvdRkT0dQCfBXAFET1JRIfsW39HRA8DeBgmJ/0Xgee9o+N1rzKbfO/AZt//Jt87sNn3L/feEeqr+7YgCIIgLAvDtJsWBEEQhAUiYicIgiCsPUsrdpvciizCvV9ERN8gouNEdC8RbQ1z5XFocP+H7TbHiOiw9/pXiegx77M/e7irbwcRXW2v+TgRzXUiIqJt+1ket5/thd57N9nXH/Pmw1eGtvdORBcS0Snvc7596GvvSsC9v42IvkVEYyK6Nvde4d//KtHx/jPvs7+v9mTOGmLZfgD8Fcwic8C0Ift4yXZXAHgHgC/mXv80gGsXfR8Luvd/APBu+/h2ADcu+p5i3z+A0wE8bv89zT4+zb73VQAHF30fDe43AfB9ABfDrDn9NoBLc9v8EYDb7eN3A7jXPr7Ubr8N4CJ7nGTR9zTQvV8I4JFF30PP934hgF8EcI8/nlX9/a/KT5f7t++91OR8SxvZYbNbkbW+d7uM43IAn6vbf4kJuf9DAB5g5h8z808APADg6oGuLza/CuA4Mz/OzHsAPgPzO/Dxfyefg6lsJvv6Z5h5l5l/AOC4Pd6q0OXeV53ae2fmHzLzdwDo3L7r8Pff5f4bs8xiF6sV2XeI6FYi2o54bX3T5d7PAPACM4/t8ycBvD7mxQ1AyP2/HsAT3vP8fX7Kpjf+bAUGxrp7mdnGfrb/A/NZh+y7zHS5dwC4iIj+k4j+mYje2vfFRqbLZ7fqnzvQ/R52iOhBIvo3Iqr9Qr9QPztaklZki6Dne196er7/65j5KSJ6LYB/hHHcuKfdlQpLzNMAfpaZnyeiXwbweSL6BWb+30VfmDAIb7D/n18M4MtE9DAzf79s44WKHW9wK7Ie7/15AK8jotR+Cz4fwFMdLzc6Ee7/KQC/6T0/H2auDsz8lP33RSL6e5h0yTKL3VMALvCeF31mbpsniSgF8NMwn3XIvstM63tnM3GzCwDM/E0i+j6AnwcQpZfiAHT57Er//leITn+73v/njxPRVwH8EswcYCHLnMbc5FZkre/dDgBfAeAqlxr/7paAkPu/H8Dbieg0W635dgD3E1FKRGcCABGNAPwulv+z/w8Al9gq2i2YIox8dZn/O7kWwJftZ30fgHfbisWLAFwC4N8Huu4YtL53IjqLiBIAsN/uL4Ep1FgVQu69jMK//56usy9a37+97237+EwAvw7g0cqdFl2RU1GpcwaM7c8xAF8CcLp9/SCAO73tvg7gJIBTMDnfQ/b1L8O0IXsEwN8CeM2i72nAe78YZsA7DtOmbXvR99TT/V9v7/E4gD+wr70awDcBfAfAdwH8NVagOhHA7wD4L5hvpjfb1z4K4J328Y79LI/bz/Zib9+b7X6PAfjtRd/LUPcO4PfsZ/wQgG8BeMei76WHe/8V+//2yzCR/He9fef+/lftp+39A/g1O75/2/77/rpzSbswQRAEYe1Z5jSmIAiCIERBxE4QBEFYe0TsBEEQhLVHxE4QBEFYe0TsBEEQhLVHxE4QBEFYe0TsBEEQhLXn/wEiwuUYyUMLWAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "\"\"\"Plotting the test orbit\"\"\"\n", + "fig, ax = plt.subplots(1,1,figsize=(7,7)) \n", + "\n", + "sim = setup_sim()\n", + "patches = [[] for i in range(sim.N)]\n", + "steps = 200\n", + "xyz = np.zeros((steps,sim.N,3))\n", + "for k in range(steps):\n", + " sim.step_manual(\"ias15\")\n", + " for i in range(sim.N):\n", + " refid = 1\n", + " ref = sim.particles[refid]\n", + " xyz[k][i] = [sim.particles[i].x-ref.x, sim.particles[i].y-ref.y, sim.particles[i].z-ref.z]\n", + " if rcrit[i]>0.:\n", + " circle = Circle((xyz[k][i][0], xyz[k][i][1]), rcrit[i])\n", + " patches[i].append(circle)\n", + " \n", + "for i in range(sim.N): \n", + " p = PatchCollection(patches[i], alpha=0.02)\n", + " ax.add_collection(p)\n", + "\n", + "ax.set_aspect(\"equal\")\n", + "ax.set_xlim([-0.15,0.15])\n", + "ax.set_ylim([-0.15,0.15])\n", + "for i in range(sim.N):\n", + " ax.scatter(xyz[:,i,0],xyz[:,i,1]);" + ] + }, + { + "cell_type": "code", + "execution_count": 212, + "metadata": {}, + "outputs": [], + "source": [ + "def symplecticity(sim,algo):\n", + " Ndim = sim.N*6\n", + " Omega1 = np.concatenate((np.zeros((Ndim//2,Ndim//2)),np.identity(Ndim//2)),axis=1)\n", + " Omega2 = np.concatenate((-np.identity(Ndim//2),np.zeros((Ndim//2,Ndim//2))),axis=1)\n", + " Omega = np.concatenate((Omega1,Omega2))\n", + " J = np.zeros((Ndim,Ndim))\n", + " delta = 1e-9\n", + " sim_0 = sim.copy()\n", + " sim_0.move_to_hel()\n", + " sim_0.step_manual(algo)\n", + " sim_0.move_to_com()\n", + " for i in range(Ndim):\n", + " sim_p = sim.copy()\n", + " if i