diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/.idea/.gitignore b/examples/symba_swifter_comparison/1pl_1tp_encounter/.idea/.gitignore new file mode 100644 index 000000000..26d33521a --- /dev/null +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/cb.swiftest.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/cb.swiftest.in new file mode 100644 index 000000000..96c7f920c Binary files /dev/null and b/examples/symba_swifter_comparison/1pl_1tp_encounter/cb.swiftest.in differ diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py b/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py new file mode 100755 index 000000000..b292ed42f --- /dev/null +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py @@ -0,0 +1,178 @@ +#!/usr/bin/env python3 +""" +For testing RMVS, the code generates clones of test particles based on one that is fated to impact Mercury. +To use the script, modify the variables just after the "if __name__ == '__main__':" line +""" +import numpy as np +import swiftest +from scipy.io import FortranFile +import sys + +swifter_input = "param.swifter.in" +swifter_pl = "pl.swifter.in" +swifter_tp = "tp.swifter.in" +swifter_bin = "bin.swifter.dat" +swifter_enc = "enc.swifter.dat" + +swiftest_input = "param.swiftest.in" +swiftest_pl = "pl.swiftest.in" +swiftest_tp = "tp.swiftest.in" +swiftest_cb = "cb.swiftest.in" +swiftest_bin = "bin.swiftest.dat" +swiftest_enc = "enc.swiftest.dat" + +MU2KG = swiftest.MSun +TU2S = swiftest.YR2S +DU2M = swiftest.AU2M + +GMSun = swiftest.GMSunSI * TU2S**2 / DU2M**3 + +# Simple initial conditions of a circular planet with one test particle in a close encounter state +# Simulation start, stop, and output cadence times +t_0 = 0 # simulation start time +deltaT = 0.25 * swiftest.JD2S / TU2S # simulation step size +end_sim = 0.15 +t_print = deltaT #output interval to print results + +iout = int(np.ceil(t_print / deltaT)) +rmin = swiftest.RSun / swiftest.AU2M +rmax = 1000.0 + +npl = 1 +plid = 2 +tpid = 100 + +radius = np.double(4.25875607065041e-05) +mass = np.double(0.00012002693582795244940133) +apl = np.longdouble(1.0) +atp = np.longdouble(1.01) +vpl = np.longdouble(2 * np.pi) +vtp = np.longdouble(2 * np.pi / np.sqrt(atp)) + +p_pl = np.array([apl, 0.0, 0.0], dtype=np.double) +v_pl = np.array([0.0, vpl, 0.0], dtype=np.double) + +p_tp = np.array([atp, 0.0, 0.0], dtype=np.double) +v_tp = np.array([0.0, vtp, 0.0], dtype=np.double) + +Rhill = apl * 0.0100447248332378922085 + +#Make Swifter files +plfile = open(swifter_pl, 'w') +print(npl+1, f'! Planet input file generated using init_cond.py',file=plfile) +print(1,GMSun,file=plfile) +print('0.0 0.0 0.0',file=plfile) +print('0.0 0.0 0.0',file=plfile) +print(plid,"{:.23g}".format(mass),Rhill, file=plfile) +print(radius, file=plfile) +print(*p_pl, file=plfile) +print(*v_pl, file=plfile) +plfile.close() + +tpfile = open(swifter_tp, 'w') +print(1,file=tpfile) +print(tpid, file=tpfile) +print(*p_tp, file=tpfile) +print(*v_tp, file=tpfile) +tpfile.close() + +sys.stdout = open(swifter_input, "w") +print(f'! Swifter input file generated using init_cond.py') +print(f'T0 {t_0} ') +print(f'TSTOP {end_sim}') +print(f'DT {deltaT}') +print(f'PL_IN {swifter_pl}') +print(f'TP_IN {swifter_tp}') +print(f'IN_TYPE ASCII') +print(f'ISTEP_OUT {iout:d}') +print(f'ISTEP_DUMP {iout:d}') +print(f'BIN_OUT {swifter_bin}') +print(f'OUT_TYPE REAL8') +print(f'OUT_FORM XV') +print(f'OUT_STAT UNKNOWN') +print(f'J2 {swiftest.J2Sun}') +print(f'J4 {swiftest.J4Sun}') +print(f'CHK_CLOSE yes') +print(f'CHK_RMIN {rmin}') +print(f'CHK_RMAX {rmax}') +print(f'CHK_EJECT {rmax}') +print(f'CHK_QMIN {rmin}') +print(f'CHK_QMIN_COORD HELIO') +print(f'CHK_QMIN_RANGE {rmin} {rmax}') +print(f'ENC_OUT {swifter_enc}') +print(f'EXTRA_FORCE no') +print(f'BIG_DISCARD no') +print(f'RHILL_PRESENT yes') +sys.stdout = sys.__stdout__ + +#Now make Swiftest files +cbfile = FortranFile(swiftest_cb, 'w') +Msun = np.double(1.0) +cbfile.write_record(0) +cbfile.write_record(np.double(GMSun)) +cbfile.write_record(np.double(rmin)) +cbfile.write_record(np.double(swiftest.J2Sun)) +cbfile.write_record(np.double(swiftest.J4Sun)) +cbfile.close() + +plfile = FortranFile(swiftest_pl, 'w') +plfile.write_record(npl) + +plfile.write_record(plid) +plfile.write_record(p_pl[0]) +plfile.write_record(p_pl[1]) +plfile.write_record(p_pl[2]) +plfile.write_record(v_pl[0]) +plfile.write_record(v_pl[1]) +plfile.write_record(v_pl[2]) +plfile.write_record(mass) +plfile.write_record(radius) +plfile.close() +tpfile = FortranFile(swiftest_tp, 'w') +ntp = 1 +tpfile.write_record(ntp) +tpfile.write_record(tpid) +tpfile.write_record(p_tp[0]) +tpfile.write_record(p_tp[1]) +tpfile.write_record(p_tp[2]) +tpfile.write_record(v_tp[0]) +tpfile.write_record(v_tp[1]) +tpfile.write_record(v_tp[2]) + +tpfile.close() + +sys.stdout = open(swiftest_input, "w") +print(f'! Swiftest input file generated using init_cond.py') +print(f'T0 {t_0} ') +print(f'TSTOP {end_sim}') +print(f'DT {deltaT}') +print(f'CB_IN {swiftest_cb}') +print(f'PL_IN {swiftest_pl}') +print(f'TP_IN {swiftest_tp}') +print(f'IN_TYPE REAL8') +print(f'ISTEP_OUT {iout:d}') +print(f'ISTEP_DUMP {iout:d}') +print(f'BIN_OUT {swiftest_bin}') +print(f'OUT_TYPE REAL8') +print(f'OUT_FORM XV') +print(f'OUT_STAT REPLACE') +print(f'CHK_CLOSE yes') +print(f'CHK_RMIN {rmin}') +print(f'CHK_RMAX {rmax}') +print(f'CHK_EJECT {rmax}') +print(f'CHK_QMIN {rmin}') +print(f'CHK_QMIN_COORD HELIO') +print(f'CHK_QMIN_RANGE {rmin} {rmax}') +print(f'ENC_OUT {swiftest_enc}') +print(f'EXTRA_FORCE no') +print(f'BIG_DISCARD no') +print(f'ROTATION no') +print(f'GR no') +print(f'MU2KG {MU2KG}') +print(f'DU2M {DU2M}') +print(f'TU2S {TU2S}') + + + + + diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swifter.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swifter.in new file mode 100644 index 000000000..d1a0c9f27 --- /dev/null +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swifter.in @@ -0,0 +1,26 @@ +! Swifter input file generated using init_cond.py +T0 0 +TSTOP 0.15 +DT 0.0006844626967830253 +PL_IN pl.swifter.in +TP_IN tp.swifter.in +IN_TYPE ASCII +ISTEP_OUT 1 +ISTEP_DUMP 1 +BIN_OUT bin.swifter.dat +OUT_TYPE REAL8 +OUT_FORM XV +OUT_STAT UNKNOWN +J2 2.198e-07 +J4 -4.805e-09 +CHK_CLOSE yes +CHK_RMIN 0.004650467260962157 +CHK_RMAX 1000.0 +CHK_EJECT 1000.0 +CHK_QMIN 0.004650467260962157 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 1000.0 +ENC_OUT enc.swifter.dat +EXTRA_FORCE no +BIG_DISCARD no +RHILL_PRESENT yes diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in new file mode 100644 index 000000000..36937896f --- /dev/null +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in @@ -0,0 +1,29 @@ +! Swiftest input file generated using init_cond.py +T0 0 +TSTOP 0.15 +DT 0.0006844626967830253 +CB_IN cb.swiftest.in +PL_IN pl.swiftest.in +TP_IN tp.swiftest.in +IN_TYPE REAL8 +ISTEP_OUT 1 +ISTEP_DUMP 1 +BIN_OUT bin.swiftest.dat +OUT_TYPE REAL8 +OUT_FORM XV +OUT_STAT REPLACE +CHK_CLOSE yes +CHK_RMIN 0.004650467260962157 +CHK_RMAX 1000.0 +CHK_EJECT 1000.0 +CHK_QMIN 0.004650467260962157 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 1000.0 +ENC_OUT enc.swiftest.dat +EXTRA_FORCE no +BIG_DISCARD no +ROTATION no +GR no +MU2KG 1.988409870698051e+30 +DU2M 149597870700.0 +TU2S 31557600.0 diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swifter.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swifter.in new file mode 100644 index 000000000..95513c9fd --- /dev/null +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swifter.in @@ -0,0 +1,8 @@ +2 ! Planet input file generated using init_cond.py +1 39.476926408897625196 +0.0 0.0 0.0 +0.0 0.0 0.0 +2 0.00012002693582795244940133 0.010044724833237891545 +4.25875607065041e-05 +1.0 0.0 0.0 +0.0 6.283185307179586 0.0 diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swiftest.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swiftest.in new file mode 100644 index 000000000..6f4bc1337 Binary files /dev/null and b/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swiftest.in differ diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb b/examples/symba_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb new file mode 100644 index 000000000..2c1c7d294 --- /dev/null +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb @@ -0,0 +1,138 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import swiftest\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swifter file param.swifter.in\n", + "Reading in time 1.355e-01\n", + "Creating Dataset\n", + "Successfully converted 199 output frames.\n", + "Swifter simulation data stored as xarray DataSet .ds\n" + ] + } + ], + "source": [ + "swiftersim = swiftest.Simulation(param_file=\"param.swifter.in\", codename=\"Swifter\")\n", + "swiftersim.bin2xr()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file param.swiftest.in\n", + "Reading in time 1.506e-01\n", + "Creating Dataset\n", + "Successfully converted 221 output frames.\n", + "Swiftest simulation data stored as xarray DataSet .ds\n" + ] + } + ], + "source": [ + "swiftestsim = swiftest.Simulation(param_file=\"param.swiftest.in\")\n", + "swiftestsim.bin2xr()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "swiftdiff = swiftestsim.ds - swiftersim.ds" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "swiftdiff = swiftdiff.rename({'time' : 'time (y)'})\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ]" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvVUlEQVR4nO3dd3hUZdrH8e+dAqGEFhIIhCZF6RBCQEAQEaWoCChNkSqLiorIuuzq+trWuq66KxZkUVQwIgqCIFhAERElIIo0pUogEDqEEjKZ+/1jRjZigEzIyZTcn+vKlTnnPM/ML2Vy57TnEVXFGGOMya8wfwcwxhgTXKxwGGOM8YkVDmOMMT6xwmGMMcYnVjiMMcb4JMLfAQpT5cqVtXbt2v6OYYwxQWPlypX7VDXWlz4hVThq165Namqqv2MYY0zQEJHtvvaxQ1XGGGN8YoXDGGOMT6xwGGOM8Ymj5zhEpBvwAhAOTFbVJ8/Y3gt4FHADLmCsqi71btsGHAVyAJeqJhUkQ3Z2NmlpaZw8ebLAX0cgioqKIiEhgcjISH9HMcYUM44VDhEJByYCXYE0YIWIzFHVdbmafQ7MUVUVkWbADOCSXNs7q+q+C8mRlpZGdHQ0tWvXRkQu5KkChqqyf/9+0tLSqFOnjr/jGGOKGScPVSUDm1R1i6qeAlKAXrkbqGqm/m+UxTJAoY+4ePLkSWJiYkKmaACICDExMSG3F2WMCQ5OFo7qwI5cy2nedb8jIr1FZAMwDxiea5MCn4jIShEZdbYXEZFRIpIqIql79+49W5uC5A9oofg1GWOCg5OFI6+/bH/Yo1DVWap6CXA9nvMdv2mvqolAd+AOEemY14uo6iRVTVLVpNhYn+5hMcaY4JXjggV/hcNpRf7SThaONKBGruUEYNfZGqvqEqCuiFT2Lu/yfs4AZuE59BVw2rVrl+f6oUOHMnPmzCJOY4wpNr54HJa/BNuXFflLO1k4VgD1RaSOiJQABgBzcjcQkXriPeYiIolACWC/iJQRkWjv+jLAVcBPDmYtsGXLiv6HZowp5n75FL56FhJvgWb9ivzlHbuqSlVdIjIGWIjnctwpqrpWREZ7t78C9AVuEZFs4ATQ33uFVRVglremRADTVXWBU1kvRNmyZcnMzERVufPOO1m0aBF16tTBZlY0xjjicBp8MAqqNIHuT+PKcXPS5aZsyaIbQcrRV1LV+cD8M9a9kuvxU8BTefTbAjR3MlthmzVrFhs3bmTNmjXs2bOHRo0aMXz48PN3NMaY/MrJhveGeT7fOBUiS/H43HUs3bSXWbe3p0wRFQ+7c7yQLFmyhIEDBxIeHk61atW44oor/B3JGBNqFt4Pad/Bdf+GyvWY/u2vTPl6Kx3qxRZZ0QArHIXKLpE1xjjmhxT47lW4dAw06cOyTft48MOfuPziWP7W45Lz9y9EVjgKSceOHUlJSSEnJ4f09HQWL17s70jGmFCR/gPMvRtqXwZXPszWfce4bdoq6lQuw78HtiQivGj/lIfUfBz+1Lt3bxYtWkTTpk1p0KABnTp18nckY0woOH4A3r0ZSsfADa9zKMvNiDdWECbw3yGtKRdV9OPVWeG4QJmZmYDnMNWLL77o5zTGmJDizoGZw+Hobhi2gFNRMYye8i1pB0/w9sg21Iwp7ZdYdqjKGGMC1aJHYcti6PksWj2RB2avYfmWAzx1Q1OS61TyWywrHMYYE4jWfQhLn4NWQyHxFl75cgszUtO4q0t9erdM8Gs0KxzGGBNoMjbA7NshoTV0f5r5a9J5asEGrm1ejXuurO/vdFY4jDEmoBw/AO8MgMjS0O9NVqef4J53V5NYswLP3NAsIC77t5PjxhgTKHKy4b0hcGQnDJ3HTndFRk79mtjokky6JYmoyHB/JwRsj8MYYwLHwvth6xK49gWOxrZkxBsryMrO4fWhralctqS/051mhcNhO3bsoHPnzjRs2JDGjRvzwgsv+DuSMSYQrXzj9J3hrqYDuPOd7/klI5OXbk6kfpVof6f7HTtU5bCIiAieffZZEhMTOXr0KK1ataJr1640atTI39GMMYFi+zKYNx7qdoErH+axeev5YuNe/tG7CZfVD7wJ6myPw2Hx8fEkJiYCEB0dTcOGDdm5c6efUxljAsahXz13hlesBTdMYeq3abyxbBsjO9Thpja1/J0uT8Vqj+PhuWtZt+tIoT5no2rl+L9rG+er7bZt2/j+++9p06ZNoWYwxgSprEx4Z6BnGtiBKSzefoqH567lyoZV+GuPhv5Od1a2x1FEMjMz6du3L88//zzlypXzdxxjjL+53TD7NshYBzdOYX12FcZMX0XD+HK8MKAF4WH+v+z2bIrVHkd+9wwKW3Z2Nn379uWmm26iT58+fslgjAkwXz4F6+fAVf8gPbY9wyYuo2xUBP8d0rpI59YoCNvjcJiqMmLECBo2bMi4ceP8HccYEwjWzoYvn4QWN3G05SiGvb6CzCwXrw9Npmr5KH+nOy8rHA77+uuveeutt1i0aBEtWrSgRYsWzJ8///wdjTGhKf1HzyGqhGSyuz/L7dO/Z1NGJi/dlEijasFxGNvR/SER6Qa8AIQDk1X1yTO29wIeBdyACxirqkvz0zdYdOjQAVX1dwxjTCA4ku4ZTqRURbT/W/x1zs989cs+nr6hGR0bBN5lt2fj2B6HiIQDE4HuQCNgoIicefPC50BzVW0BDAcm+9DXGGOCx6ljnqJx4hAMTOGFb48wc2Uad3epT7+kGv5O5xMnD1UlA5tUdYuqngJSgF65G6hqpv7v3/EygOa3rzHGBA23Gz4YBbt/hBum8N7Oijz/2S/0TUxgbACMdusrJwtHdWBHruU077rfEZHeIrIBmIdnryPffY0xJih8/hBs+Aiufpyvwlrx1w/W0KFeZZ7o0zQgRrv1lZOFI6/vxh8O9qvqLFW9BLgez/mOfPcFEJFRIpIqIql79+4taFZjjHHGyqnw9QvQeiTraw7itrdXUS+uLC/dnEiJiOC8PsnJ1GlA7gN3CcCuszVW1SVAXRGp7EtfVZ2kqkmqmhQbGzwnl4wxxcCWL2DeOKjbhfR2DzHsjVTKlozg9WGtKRcV6e90BeZk4VgB1BeROiJSAhgAzMndQETqiXc/TUQSgRLA/vz0NcaYgLZ3I7x7C8TU58i1rzFs6veeezWGtSa+fCl/p7sgjhUOVXUBY4CFwHpghqquFZHRIjLa26wv8JOIrMZzFVV/9cizr1NZnTZ8+HDi4uJo0qTJ6XUHDhyga9eu1K9fn65du3Lw4MHT25544gnq1avHxRdfzMKFC/0R2RhzIY7tg+n9IKIE2QNSuH3mJjZlZPLyzYk0jA+OezXOxdEDbKo6X1UbqGpdVf2Hd90rqvqK9/FTqtpYVVuo6qW/3cNxtr7BaujQoSxYsOB365588km6dOnCL7/8QpcuXXjySc9tKuvWrSMlJYW1a9eyYMECbr/9dnJycvwR2xhTEK4sSLkJju5GB0znr4sOs3TTPp7o0zQgh0gviOA8MxNkOnbsSKVKlX637sMPP2TIkCEADBkyhNmzZ59eP2DAAEqWLEmdOnWoV68e3333XVFHNsYUhCp8OAZ2LIfrX+bZdeVP36txY5Ddq3EugT2SVmH7eALsXlO4z1m1KXT3/ab2PXv2EB8fD3jm7MjIyABg586dtG3b9nS7hIQEm7/DmGDx5VOwZgZc8XemHknkxcVrGZhcIyjv1TgX2+MIMHkNTxKM13kbU+x8Pw2+eAKaD2Je+UE85J1X49FeTULuPVy89jgKsGfglCpVqpCenk58fDzp6enExcUBnj2MHTv+d+9jWloa1apV81dMY0x+bPoc5t4FF13O8sYPcs/U1bSqWZEXB7UkIjz0/j8Pva8oSFx33XVMnToVgKlTp9KrV6/T61NSUsjKymLr1q388ssvJCcn+zOqMeZcdq+BGUOg8sVs6DiRW6f9SK2Y0kwekkRUZLi/0zmieO1x+MnAgQP54osv2LdvHwkJCTz88MNMmDCBfv368d///peaNWvy3nvvAdC4cWP69etHo0aNiIiIYOLEiYSHh+YvnzFB73AaTLsRSkazq+ebDH57A2WjIpg6PJkKpUv4O51jJJSG/E5KStLU1NTfrVu/fj0NGwbu3L0XIpS/NmMC3olDMKUbHNnJoQFz6PP+YfZlZjHztnY0qBLt73T5JiIrVTXJlz62x2GMMb5yZcG7N8P+TZzs/y5D5h9n56ETTBvZJqiKRkHZOQ5jjPHFb/dqbPsK17X/ZvSyaNakHeLFQYkk1a50/v4hoFgUjlA6HPebUPyajAkKnz8Ca2agnR/gvl8a8cXGvTzeuyldG1Xxd7IiE/KFIyoqiv3794fUH1pVZf/+/URFBf6k9saElNQpsPRf0GooTx7ryQerdjKuawMGJNf0d7IiFfLnOBISEkhLSyPU5uqIiooiISHB3zGMKT42LoB590L9q5lSfgyvzv+ZwW1rcecV9fydrMiFfOGIjIykTp06/o5hjAlmaakwcxhUbcaH9R/lkQ9+pnuTqjx0XeOQuys8P0L+UJUxxlyQvT977tUoG8eXSRMZN3sz7evF8Fz/FoSHFb+iAVY4jDHm7A7vhLf7QFg4qzq9zq2zdtCkenleHRy6d4XnR8gfqjLGmAI5cRDe7gsnDvFzjxRumbWXWpVK88bQ1pQtWbz/dBbvr94YY/KSfQKmD4ADm9nZ8y36zzlOhdKRvDWiDRXLhO5QIvllhcMYY3LLccF7w2DHt+zv8Sp9F0QSHqa8PaINVcvbJfBg5ziMMeZ/VOGju+Hnjzna5QluWFKF46dcvDUimdqVy/g7XcCwwmGMMb/5/BH4/m1OtruX/quasPvwSV4f1pqG8eX8nSygOFo4RKSbiGwUkU0iMiGP7TeJyI/ej2Ui0jzXtm0iskZEVotI6pl9jTGmUC1/GZb+C1eLIQze3IVfMo7yyuBWtKpVPMaf8oVj5zhEJByYCHQF0oAVIjJHVdflarYV6KSqB0WkOzAJaJNre2dV3edURmOMAWDNTFgwAffF1zDqwEBSfz3AiwMT6dQg1t/JApKTexzJwCZV3aKqp4AUoFfuBqq6TFUPeheXAzaGhjGmaG36DGaNRmu2Y6xrDIt+PsA/rm9Kz2bx/k4WsJwsHNWBHbmW07zrzmYE8HGuZQU+EZGVIjLqbJ1EZJSIpIpIaqiNR2WMcdj2byDlZjTuEu6P+htz1h7ggZ4NGdSmeA1a6CsnL8fN6178PIeoFZHOeApHh1yr26vqLhGJAz4VkQ2quuQPT6g6Cc8hLpKSkkJnCFxjjLPSf4Dp/dDy1Xmi0uNMX3WEcV0bMPKyi/ydLOA5uceRBtTItZwA7DqzkYg0AyYDvVR1/2/rVXWX93MGMAvPoS9jjLlw+36Bt/qgJcvxn2rPMGnVUUZ3qlssR7otCCcLxwqgvojUEZESwABgTu4GIlIT+AAYrKo/51pfRkSif3sMXAX85GBWY0xxcWgHvHk9iPB63ef414oTDG1Xm790u7hYjnRbEI4dqlJVl4iMARYC4cAUVV0rIqO9218BHgRigJe8PzCXd9L0KsAs77oIYLqqLnAqqzGmmMjMgDd7QdZRZjR9lUeWZtMvKYEHr2lkRcMHEkoz4yUlJWlqqt3yYYzJw4mD8Ma1cGAz81q8xB1fleC65tWK9fDoACKy0vsPe77ZWFXGmNB36hhM6wf7NvJFq/9wx5ISdG1UhWf7NS/WRaOgrHAYY0KbKwtSboKdqXzX+jmGfRVNxwaxvDioJZHhNupSQVjhMMaErhwXvD8CtixmTdITDFwaR+vaFXn15laUjCi+EzFdKCu3xpjQ5M6BD++A9XPZ0OJv9F5Wm2YJ5ZkytDWlSljRuBBWOIwxocftho/Gwo8pbGpyD9euaErj6uWZOjy52M/eVxjsO2iMCS2q8PGfYdWbbGt0Oz2+b0PD+GjeHJ5MuahIf6cLCbbHYYwJHaqw8H5YMZkdDUdy1Q+X0aBqWd4a3obypaxoFBYrHMaY0KDqmYhp+UR2XnwLXdZ0oV5cNG+PaEP50lY0CpMVDmNMaPjyaVj6L3bXG8AVa7tzUeWyTBvZhgqlS/g7WcixwmGMCX5Ln4MvHiejbl86b7yO2jGeolGxjBUNJ9jJcWNMcPvmJfjsIfbVvpbOP99AQsUyTLu1DTFlS/o7WciyPQ5jTPBaMRkW/pUDtbrRectAqlYozfRb21LZioajrHAYY4LTqrdg3r0cSOhC562DiS1XhndubUtstBUNp1nhMMYEn1VvwZw7ORh/GZdvH0JshWhS/tSWuHJR/k5WLFjhMMYEl1VvwpwxHIjvQMcdtxJfqSIpo9oSF21Fo6jYyXFjTPBY+QbMvZv98R3ptGMkNWIrMW1kGyrZ1VNFyvY4jDHBIfV1mHs3+6p2pNOvI6kVV4npVjT8wgqHMSbwpU6Bj8aSUbUTHX8dQd34GKaPbGv3afiJHaoyxgS2Ff+FeePIqNKJTr+O4JKEyky1AQv9ytE9DhHpJiIbRWSTiEzIY/tNIvKj92OZiDTPb19jTDGwYjLMG8fuKp3otGMEjWvE2ii3AcCxwiEi4cBEoDvQCBgoIo3OaLYV6KSqzYBHgUk+9DXGhLLvXoN595JepROX/zqcZrXimDo8mWgrGn7n5B5HMrBJVbeo6ikgBeiVu4GqLlPVg97F5UBCfvsaY0LYd6/B/PGkxXqKRuJFVXl9WGvK2CRMAcHJwlEd2JFrOc277mxGAB8XsK8xJlQsexHmj2db5U503jGCDhdXZ8rQ1pQuYUUjUDj5k5A81mmeDUU64ykcHQrQdxQwCqBmzZq+pzTGBAZVWPJPWPwYG2O6cE3aEK5uVoPn+rcgMtwuAA0kTv400oAauZYTgF1nNhKRZsBkoJeq7velL4CqTlLVJFVNio2NLZTgxpgipgqfPwyLH+OHSt3osXMovZNq88KAllY0ApCTP5EVQH0RqSMiJYABwJzcDUSkJvABMFhVf/alrzEmRKjCggmw9Dm+rXQd1++6mVva1+XJPs0ID8vr4IPxN8cOVamqS0TGAAuBcGCKqq4VkdHe7a8ADwIxwEsiAuDy7j3k2deprMYYP3HnwEf3wKqpLK54A8N29ebOK+ozrmsDvH8TTAAS1TxPHQSlpKQkTU1N9XcMY0x+5Lhg9m2wZgYflR/EmD09mdC9IaM71fV3smJFRFaqapIvfewyBWNM0XOdgvdHwPo5zCg3hPv2XM2j1zdhcNta/k5m8sEKhzGmaGWfhBm3wC8LmVJ2FI/tvZx/9WtOn8SE8/c1AcEKhzGm6GRlQsogdOsS/h11OxMPduSlm1rQrUm8v5MZH1jhMMYUjeMHYNqN6K5VPBZ5JynH2/PGsCTa1avs72TGR1Y4jDHOO7IL3uqN+8BWxst4vnQnkzIqmaYJ5f2dzBSAFQ5jjLP2b4Y3r8d1/AAjs//CpjItmTmiDXUql/F3MlNAVjiMMc5J/wHe7ktWtosBx//G8cpNeX9EMlXK2fzgwSxfd46LyIgzlsNF5P+ciWSMCQnbvoY3riEzJ4IemfcTntCSGX+61IpGCMjvkCNdRGS+iMSLSBM8Q6BHO5jLGBPMNn6Mvt2H/WExdD30N2o3aMFbI9pQvrTNpREK8nWoSlUHiUh/YA1wHBioql87mswYE5xWv4N+eAc7SzXg2gN3c0ViI57s29QGKwwh+T1UVR+4G3gf2AYMFpHSDuYyxgSjb16C2aPZGNWcqw+M58aOLfjnjc2saISY/J4cnwvcoaqfi2fksXvwjGDb2LFkxpjgoQqfPwJL/8U3Jdsz5OAo7uvZjJGXXeTvZMYB+S0cycBIEbkDz4RKS/EMdW6MKe5cp2DOnfBjCnMirua+Y0N4blAiPZvZ3eChKr+FYzJwFPiPd3kgcCnQz4lQxpggcfIIzBgMW75gYthAJrl689bI1rSuXcnfyYyD8ls4LlbV5rmWF4vID04EMsYEiSPpMO1G3Bnr+Zv7dpaW6sr7w5KpF1fW38mMw/JbOL4XkbaquhxARNoAdlWVMcXV3o3wdl+yM/czIms8B+I78MHQ1sRF2z0axUF+C0cb4BYR+dW7XBNYLyJrAFXVZo6kM8YEnu3foO8M4FhOOP2P309sg2TeHZRImZI2EEVxkd+fdDdHUxhjgsO6D9H3byUjvAp9M++lQ1Iij13fhAi73LZYye8NgNudDmKMCXDfvop+/Bd+jmxI/yN3M6JrK8ZcUc/mBi+GbN/SGHNu7hz49EH45kWWhrfh9uN38Ej/VvRuaTP2FVeO7l+KSDcR2Sgim0RkQh7bLxGRb0QkS0TGn7Ftm4isEZHVIpLqZE5jzFmcOgbvDoZvXuQdunGPjuP1Wy+zolHMObbHISLhwESgK5AGrBCROaq6LlezA8BdwPVneZrOqrrPqYzGmHM4sgveGYA7fQ2PuobwVaW+fDCkNTVjbLSh4s7JQ1XJwCZV3QIgIilAL+B04VDVDCBDRHo6mMMY46v0H9DpAzh17BB/OnUvOXW78v6gRMqXstFtjbOHqqoDO3Itp3nX5ZcCn4jIShEZdbZGIjJKRFJFJHXv3r0FjGqMOW3jx+iU7hw87qLXiQdJSO7F60NbW9Ewpzm5x5HXpRbqQ//2qrpLROKAT0Vkg6ou+cMTqk4CJgEkJSX58vzGmNxUYflL6ML7+Tm8HrccH8tt17RnSLvaduWU+R0nC0caUCPXcgKwK7+dVXWX93OGiMzCc+jrD4XDGFMIcrLh4/sgdQqLpS33ue7gmSFt6XxJnL+TmQDkZOFYAdQXkTrATjyj6Q7KT0cRKQOEqepR7+OrgEccS2pMcXbyMLw3FDYvYpL7Ot4sNYS3hyVzSdVy/k5mApRjhUNVXSIyBlgIhANTVHWtiIz2bn9FRKoCqUA5wC0iY4FGQGVglnf3OAKYrqoLnMpqTLG1fzP6zkDc+zcxIXsU22v0YfbNiVQuW9LfyUwAc/QGQFWdD8w/Y90ruR7vxnMI60xHgOZ5rDfGFJbNi9AZQ8nMVkZlTaBucnfevqYxJSJs+BBzbnbnuDHFzW8nwT95gK0kMOLUvYy6/goGJtf0dzITJKxwGFOcuLLgo3tg9TQ+12QejriL525tT5JNvGR8YIXDmOLi6G405WZk5wqed/VhcZVhzLilNfHlS/k7mQkyVjiMKQ52rsT9ziCyjx3i7lNjKdW8N+/2aUpUZLi/k5kgZIXDmFC3+h3cc+8mw12e4Vn/R5/u3RjRoY7d1GcKzAqHMaHKdQoW/hVWTGalNuK+sHE8MuxyLqsf6+9kJshZ4TAmFB3Zhc64BUlbwWuuHsyvMpq3bm5NQkUb2dZcOCscxoSabV+TM2MI2SeOMv7UnVRI7k/KNY0oGWHnM0zhsMJhTKhQheUvo588wE6N4/acRxl+Qw/6JNqkS6ZwWeEwJhScOobOuQv5aSafuZN4vuw9PHtLRxtvyjjCCocxwW7fJnLeHYzsXc8z2f3YcvGtvNOvJeWibP4M4wwrHMYEszUzyZlzF5muMO7KnkD7q2/k5csuskttjaOscBgTjLJPogv+iqycwmq9mAcj7+WB4Vdyad0YfyczxYAVDmOCzf7N5MwYQvieNbziupbltW/jjf5JxEbbUOimaFjhMCaYrJ1FzuwxHMuGcdl/JrHrAKZ0rEtYmB2aMkXHCocxwcCVhS68H1nxGj+66/NI1HjuH3qVjWpr/MIKhzGBbv9mXDOGEbHnB15z9SC13l283q8VFUqX8HcyU0xZ4TAmUKnC6unkzBvPMVcYf3HdS3K3wbzSvrZdNWX8ygqHMYHo5GFy5owlfN0HrHA35Nmy4/n7oCtpllDB38mMscJhTMDZsYLsGcMIO7qTZ7L7cTDxdqZe25TSJeztagKDo7PSi0g3EdkoIptEZEIe2y8RkW9EJEtExvvS15iQ485Bv3wG95Sr2XPkJMPlEZoOfJTH+7a0omECimO/jSISDkwEugJpwAoRmaOq63I1OwDcBVxfgL7GhI7DO8meeSuRO75mTk475tX8M0/3b0eVclH+TmbMHzj5b0wysElVtwCISArQCzj9x19VM4AMEenpa19jQoIqrHkP19xxZGdn84DrNupfdSuvdLjI7s0wAcvJwlEd2JFrOQ1oU9h9RWQUMAqgZs2avqc0xl+OHyB7zlgiN3zIancD/lNuHH8Z1JNG1WxEWxPYnCwcef27pIXdV1UnAZMAkpKS8vv8xvjXL59y6oPbkRMHeNrVH/eld/HqVQ2JirTJlkzgc7JwpAE1ci0nALuKoK8xgSsrE9fCB4hY9Tpb3DX4Z+mnGT2gt90BboKKk4VjBVBfROoAO4EBwKAi6GtMYNrxHVkzRhJ59FdedfUko9W9/Ltnc7tiygQdx35jVdUlImOAhUA4MEVV14rIaO/2V0SkKpAKlAPcIjIWaKSqR/Lq61RWYxyVfQLXon8Q9s1EMjSGJ0o+yk03D6R9vcr+TmZMgYhq6JwWSEpK0tTUVH/HMOZ/fl3OyZmjiTqylemuzmxo+hf+3CuJaJudzwQIEVmpqkm+9LF9ZGOccOoY2Z88TETqJPZqZZ4p8TB9Bg5i0MVx/k5mzAWzwmFMYdu6hBPv306pzB284bqanYnjebxnImVL2tvNhAb7TTamsGQdJevjByi5+g12u6vwfJnHubn/QIbaFVMmxFjhMKYQ6Ib5nPxwHCVP7GZyTg+OXvoXnura1O7LMCHJCocxF+LILo5/eC+lN89nu7sGkyo8w4gB/Whcrby/kxnjGCscxhSEOwfXt6/h/uwRxJXNszqICl3G8nSH+kSEOzrotDF+Z4XDGF+l/0DmzDsou38NX+Y049M693F7nyupVqGUv5MZUySscBiTX1mZnPjkUUqunMQJjeaZEuPo1P9PPNawqr+TGVOkrHAYcz6quH/6gJMfTaB0VgbTc7qwt80EJlzVilIl7OS3KX6scBhzLnvWcWTWPZTbvZwt7tqkxP6LW/r1o0GVaH8nM8ZvrHAYk5eTh8lc+Bilvp9Mjpbm6cg/0fCaO3m0eQIiNsGSKd6scBiTm9tN9qppZC98kNKnDvKuduFgm/sY07WVjWJrjJe9E4zx0rRUDn9wLxUOrGaNux4Laj7O4D69qFGptL+jGRNQrHAYc+hXDs39OxU2zyZby/N0qbtp12cMf2tgAxIakxcrHKb4OnmEI589RamVrxLlhslhfSl1xb3c074hkXYTnzFnZYXDFD85Lk58OwX3oscp5zrIh+4O7Gp1Hzdd3Y5yNk+GMedlhcMUH6pkb/yEzLkTqHhsC9+5L2Fp3ScZ0KuX3fVtjA+scJhiwb19OQfm3E/l/akcdlfhjdgH6dp7BOMSKvg7mjFBxwqHCWm6ew17P/w7cemLUS3PxNKjaXzNnYxtVN3uxzCmgBwtHCLSDXgBCAcmq+qTZ2wX7/YewHFgqKqu8m7bBhwFcgCXr3PimuJN929h75z/o/L2uZTUUrxWcjDVrx7LbS3rEhZmBcOYC+FY4RCRcGAi0BVIA1aIyBxVXZerWXegvvejDfCy9/NvOqvqPqcymhB0ZBcZ8x6j0sYUojWcaRG9KdvlXoa1aWTDnRtTSJzc40gGNqnqFgARSQF6AbkLRy/gTVVVYLmIVBCReFVNdzCXCUWH09jz8ZNU3JBCBXXzQXhXpON4+nVIpGSEDURoTGFysnBUB3bkWk7j93sTZ2tTHUgHFPhERBR4VVUn5fUiIjIKGAVQs2bNwklugoYe3M6e+U8S88t7VFI3c8M6c6Lt3fS5or2NXGuMQ5wsHHkdSFYf2rRX1V0iEgd8KiIbVHXJHxp7CsokgKSkpDOf34QoPbCV9HmPE7f5fSopzAnvguvSsVzXqa0VDGMc5mThSANq5FpOAHblt42q/vY5Q0Rm4Tn09YfCYYoXd8bPpM9/gqrbZhOjYXwYcRV0uIeeHZKIirSCYUxRcLJwrADqi0gdYCcwABh0Rps5wBjv+Y82wGFVTReRMkCYqh71Pr4KeMTBrCbAZW1bzt6Pn6bankXEaASzIntQotNYrr20FSUi7KS3MUXJscKhqi4RGQMsxHM57hRVXSsio73bXwHm47kUdxOey3GHebtXAWZ5r7OPAKar6gKnspoA5XaT+dM8Dn/2LNWPfE9ZLcO7pfpRodMdXJ/c1K6SMsZPxHNBU2hISkrS1NRUf8cwF8qVxYHl08hZ+gKxJ7eRppX5ouKN1L36NtpeUtNu3DOmEInISl/vk7M7x03A0KN72Pn5y5Rd8yaVcvazXmuxMOHvJPUYxs3VY/wdzxjjZYXD+N3Jbd+R/skLVN+1gARcfE1zdlz8dy7vPoCbbfBBYwKOFQ7jH65T7P32XbKWvUzCsbXEahQLoroR3mYUV3ToQHu7pNaYgGWFwxQp1/5t7Pj8VSpteIdY90G2aDzvxd3JRV1v5dr6dv7CmGBghcM4z3WKvakfcOyb/1Lz8ApqKSwLa0lGo4dof3U/bixvc3obE0yscBjHZKWv59fPXqHK1lnEug+TrTHMqTiYCu2G0b5VC5ue1ZggZYXDFCo9cZC0r9/F/f00ah37kdoazrKI1hxuMpDkLjdyfcUy/o5ojLlAVjjMhXOdYs+quRz+dhq19y+hBtls1ap8GPcnqnUcxmWNL7E5MIwJIVY4TMGocnDjV+z+6k2q71pAFT1KuJZjUdkehLfoT3L7rvQqXcLfKY0xDrDCYfLP7ebAL8vYvXwGsb8uIDZnD1FaguUlL+Vkwxto2ak33SpF+zulMcZhVjjMublzyFi7mL3fvkfVXZ8S495PWQ1nVURzVtS7jYsvH0DnhHh/pzTGFCErHOYPNPsEu1Z/yqFVs6i2exFxeohyGsnKEq04Wu8uLmrflzY1E+yeC2OKKSscBoDjGdvYvnw2smkhtY+kUp1TVNCSfB+VzIl619CgfR/aV4vzd0xjTACwwlFMqSuLXT8tYd/3HxGzczEJru00BHZoHMvK90AaXMUll3anQ0wlf0c1xgQYKxzFRY6LPT9/x54fPqHEjq+pdewHqpNFrIbzU0QT1tfoRUyLa2jcrBU1Iu3XwhhzdvYXIlS5c8jYtIr0Hz4l4tevqHV0NVU4ThVgMwl8W74HclFH6ib3JDG+ir/TGmOCiBWOEHHycAa//riEzE3LKJ2xihon1hPHSeKAX6nKqnKd0VqXUb3lVdStcxF17cS2MaaArHAEoVPHj7BrYyqHtqTCzpXEHvqR6u5dNABcGsaWsNqsqtANd0Jr4ptdSb16F1PT7tw2xhQSKxyBTJXDGTvYufFbjm1fTeTetcRmbiQ+J53a4pnyd7+WY2upRmyKvZ4yddtRq1l7GlSqRAM/RzfGhC4rHAEgO+s4u7eu5eD2dWTt2Uj4gU1EH9tG1ewdlOc45b3t0qjCrlL12RxzDSWqN6dKg9bUqFWPpAib9MgYU3QcLRwi0g14AQgHJqvqk2dsF+/2HsBxYKiqrspP32By8vhR9u/czKH0rZzYt42cgzuIOLqTUifSqXgqnSruvdQQpYa3/W5iyChRgzUVryKncgOiayWScElrEipXJsGvX4kxxjhYOEQkHJgIdAXSgBUiMkdV1+Vq1h2o7/1oA7wMtMln3yKlbjcnj2dy/NhhTmYe5uSxo5w6cZjs40fJPnYI19EM9Nhewo7vJzLrAFGnDlI25xDl3IcpTybVgere58pRYa/EcCAijrQyTdhS4SIi4hpQLqER8Rc1pmqFSlT11xdqjDHn4eQeRzKwSVW3AIhICtALyP3HvxfwpqoqsFxEKohIPFA7H30LzaZHEynpPkEYOYSpm3ByCMNNOG7CyCFCcyhFFqVEKXWO58lR4ZCU40hYeY5FVGRP6frsLFkJd3Q8ERVrUjq2FhXjL6Jy9dpULVHSioMxJig5WTiqAztyLafh2as4X5vqZ1l/Zl8ARGQUMAqgZs2aBQp6qMxFiLpQCUclApUwNCwCJAzCIjyPI8tAiTJIyWjCosoSHhVNZKloIktHE1W2IuVi4ilfKY6Y8HBiCpTCGGOCg5OFI6/rPzWfbfLT17NSdRIwCSApKSnPNueTNG5mQboZY0yx5GThSIPT53sBEoBd+WxTIh99jTHG+EGYg8+9AqgvInVEpAQwAJhzRps5wC3i0RY4rKrp+exrjDHGDxzb41BVl4iMARbiuaR2iqquFZHR3u2vAPPxXIq7Cc/luMPO1deprMYYY/JPPBc0hYakpCRNTU31dwxjjAkaIrJSVZN86ePkoSpjjDEhyAqHMcYYn1jhMMYY4xMrHMYYY3wSUifHRWQvsL2A3SsD+woxTlGwzEUj2DIHW16wzEUlr8y1VDXWlycJqcJxIUQk1dcrC/zNMheNYMscbHnBMheVwspsh6qMMcb4xAqHMcYYn1jh+J9J/g5QAJa5aARb5mDLC5a5qBRKZjvHYYwxxie2x2GMMcYnVjiMMcb4JOQLh4h0E5GNIrJJRCbksV1E5N/e7T+KSGJ++wZaZhGpISKLRWS9iKwVkbsDPXOu7eEi8r2IfBQMmb3THM8UkQ3e7/elQZD5Hu/vxU8i8o6IRAVI5ktE5BsRyRKR8b70DbTMAf4ePOv32bs9/+9BVQ3ZDzxDsm8GLsIzOdQPQKMz2vQAPsYz62Bb4Nv89g3AzPFAovdxNPBzoGfOtX0cMB34KNB/N7zbpgIjvY9LABUCOTOe6Zi3AqW8yzOAoQGSOQ5oDfwDGO9L3wDMHMjvwTwz59qe7/dgqO9xJAObVHWLqp4CUoBeZ7TpBbypHsuBCiISn8++AZVZVdNVdRWAqh4F1uP5gxGwmQFEJAHoCUwugqwXnFlEygEdgf8CqOopVT0UyJm92yKAUiISAZSmaGbVPG9mVc1Q1RVAtq99Ay1zIL8Hz/F99vk9GOqFozqwI9dyGn/8IZ6tTX76OuFCMp8mIrWBlsC3hR/xDy408/PAfYDboXx5uZDMFwF7gde9u/aTRaSMk2HPk+e8bVR1J/BP4FcgHc9sm584mPWceYqg74UolNcNwPfguTyPD+/BUC8ckse6M68/Plub/PR1woVk9mwUKQu8D4xV1SOFmO1sCpxZRK4BMlR1ZeHHOqcL+T5HAInAy6raEjgGFMXx9wv5PlfE8x9oHaAaUEZEbi7kfHm5kPdRIL8Hz/0EgfkezLtjAd6DoV440oAauZYT+OPu+dna5KevEy4kMyISiecXdpqqfuBgznzlyUeb9sB1IrINz+71FSLytnNRz5snP23SgDRV/e0/yZl4ConTLiTzlcBWVd2rqtnAB0A7B7OeL4/TfS/EBb1uAL8Hz8b396DTJ238+YHnP8MteP7L+u2EUeMz2vTk9ycTv8tv3wDMLMCbwPPB8n0+o83lFN3J8QvKDHwFXOx9/BDwTCBnBtoAa/Gc2xA8J/fvDITMudo+xO9PNAfse/AcmQP2PXi2zGdsy9d7sMi+MH994LnK5Gc8Vxzc7103Ghid6wc90bt9DZB0rr6BnBnogGf39EdgtfejRyBnLsgvbSBkBloAqd7v9WygYhBkfhjYAPwEvAWUDJDMVfH8x3wEOOR9XO5sfQM5c4C/B8/6fc71HPl6D9qQI8YYY3wS6uc4jDHGFDIrHMYYY3xihcMYY4xPrHAYY4zxiRUOY4wxPrHCYUwevKPf3p5ruZqIzHTota4XkQfP0+afInKFE69vjK/sclxj8uAdZ+gjVW1SBK+1DLhOVfedo00t4DVVvcrpPMacj+1xGJO3J4G6IrJaRJ4Rkdoi8hOAiAwVkdkiMldEtorIGBEZ5x3wcLmIVPK2qysiC0RkpYh8JSKXnPkiItIAyFLVfSIS7X2+SO+2ciKyTUQiVXU7ECMiVYvwe2BMnqxwGJO3CcBmVW2hqn/OY3sTYBCe4az/ARxXz4CH3wC3eNtMwjOsRytgPPBSHs/THsg9DPcXeIYNARgAvK+esaXwtmt/gV+XMRcswt8BjAlSi71/6I+KyGFgrnf9GqCZd3TUdsB7IqcHLi2Zx/PE4xmi/TeT8QxvPRsYBtyaa1sGnpFtjfErKxzGFExWrsfuXMtuPO+rMOCQqrY4z/OcAMr/tqCqX3sPi3UCwlX1p1xto7ztjfErO1RlTN6O4pn6s0DUMwfDVhG5EU7PBd48j6brgXpnrHsTeAd4/Yz1DfAMUGiMX1nhMCYPqrof+FpEfhKRZwr4NDcBI0TkBzxDmuc17ekSoKXkOp4FTAMq4ikewOk5HurhGZHXGL+yy3GN8TMReQGYq6qfeZdvAHqp6uBcbXoDiar6dz/FNOY0O8dhjP89jmeiJUTkP0B3PHMr5BYBPFvEuYzJk+1xGGOM8Ymd4zDGGOMTKxzGGGN8YoXDGGOMT6xwGGOM8YkVDmOMMT75fxGpsOxzhiZMAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "swiftdiff['px'].plot.line(x=\"time (y)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "swiftestOOF", + "language": "python", + "name": "swiftestoof" + }, + "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.7.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swifter.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swifter.in new file mode 100644 index 000000000..9c026369e --- /dev/null +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swifter.in @@ -0,0 +1,4 @@ +1 +100 +1.01 0.0 0.0 +0.0 6.252003053624663 0.0 diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swiftest.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swiftest.in new file mode 100644 index 000000000..e1506974a Binary files /dev/null and b/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swiftest.in differ diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/.idea/.gitignore b/examples/symba_swifter_comparison/9pl_18tp_encounters/.idea/.gitignore new file mode 100644 index 000000000..26d33521a --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/cb.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/cb.in new file mode 100644 index 000000000..81c636655 --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/cb.in @@ -0,0 +1,5 @@ +0 +0.00029591220819207774 +0.004650467260962157 +4.7535806948127355e-12 +-2.2473967953572827e-18 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/cb.swiftest.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/cb.swiftest.in new file mode 100644 index 000000000..81c636655 --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/cb.swiftest.in @@ -0,0 +1,5 @@ +0 +0.00029591220819207774 +0.004650467260962157 +4.7535806948127355e-12 +-2.2473967953572827e-18 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/init_cond.py b/examples/symba_swifter_comparison/9pl_18tp_encounters/init_cond.py new file mode 100755 index 000000000..321c79932 --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/init_cond.py @@ -0,0 +1,132 @@ +#!/usr/bin/env python3 +import numpy as np +import swiftest +import swiftest.io as swio +import astropy.constants as const +import sys +import xarray as xr + +# Both codes use the same tp input file +tpin = "tp.in" + +swifter_input = "param.swifter.in" +swifter_pl = "pl.swifter.in" +swifter_bin = "bin.swifter.dat" +swifter_enc = "enc.swifter.dat" + +swiftest_input = "param.swiftest.in" +swiftest_pl = "pl.swiftest.in" +swiftest_cb = "cb.swiftest.in" +swiftest_bin = "bin.swiftest.dat" +swiftest_enc = "enc.swiftest.dat" + +sim = swiftest.Simulation() + +sim.param['T0'] = 0.0 +sim.param['DT'] = 1.0 +sim.param['TSTOP'] = 365.25e1 +sim.param['ISTEP_OUT'] = 11 +sim.param['ISTEP_DUMP'] = 1 +sim.param['CHK_QMIN_COORD'] = "HELIO" +sim.param['CHK_QMIN'] = swiftest.RSun / swiftest.AU2M +sim.param['CHK_QMIN_RANGE'] = f"{swiftest.RSun / swiftest.AU2M} 1000.0" +sim.param['CHK_RMIN'] = swiftest.RSun / swiftest.AU2M +sim.param['CHK_RMAX'] = 1000.0 +sim.param['CHK_EJECT'] = 1000.0 +sim.param['OUT_FORM'] = "XV" +sim.param['OUT_STAT'] = "UNKNOWN" +sim.param['GR'] = 'NO' +sim.param['CHK_CLOSE'] = 'YES' + +sim.param['MU2KG'] = swiftest.MSun +sim.param['TU2S'] = swiftest.JD2S +sim.param['DU2M'] = swiftest.AU2M + +bodyid = { + "Sun": 0, + "Mercury": 1, + "Venus": 2, + "Earth": 3, + "Mars": 4, + "Jupiter": 5, + "Saturn": 6, + "Uranus": 7, + "Neptune": 8, +} + +for name, id in bodyid.items(): + sim.add(name, idval=id) + +ds = sim.ds +cb = ds.sel(id=0) +pl = ds.where(ds.id > 0, drop=True) +npl = pl.id.size + +ntp = 16 +dims = ['time', 'id', 'vec'] +tp = [] +t = np.array([0.0]) +clab, plab, tlab = swio.make_swiftest_labels(sim.param) + +# For each planet, we will initialize a pair of test particles. One on its way in, and one on its way out. We will also initialize two additional particles that don't encounter anything +tpnames = np.arange(101, 101 + ntp) +tpxv1 = np.empty((6)) +tpxv2 = np.empty((6)) + +p1 = [] +p2 = [] +p3 = [] +p4 = [] +p5 = [] +p6 = [] + +for i in pl.id: + pli = pl.sel(id=i) + rstart = 2 * np.double(pli['Radius']) # Start the test particles at a multiple of the planet radius away + vstart = 1.5 * np.sqrt(2 * np.double(pli['Mass']) / rstart) # Start the test particle velocities at a multiple of the escape speed + xvstart = np.array([rstart / np.sqrt(2.0), rstart / np.sqrt(2.0), 0.0, vstart, 0.0, 0.0]) + # The positions and velocities of each pair of test particles will be in reference to a planet + plvec = np.array([np.double(pli['px']), + np.double(pli['py']), + np.double(pli['pz']), + np.double(pli['vx']), + np.double(pli['vy']), + np.double(pli['vz'])]) + tpxv1 = plvec + xvstart + tpxv2 = plvec - xvstart + p1.append(tpxv1[0]) + p1.append(tpxv2[0]) + p2.append(tpxv1[1]) + p2.append(tpxv2[1]) + p3.append(tpxv1[2]) + p3.append(tpxv2[2]) + p4.append(tpxv1[3]) + p4.append(tpxv2[3]) + p5.append(tpxv1[4]) + p5.append(tpxv2[4]) + p6.append(tpxv1[5]) + p6.append(tpxv2[5]) + +tvec = np.vstack([p1, p2, p3, p4, p5, p6]) +tpframe = np.expand_dims(tvec.T, axis=0) +tpxr = xr.DataArray(tpframe, dims = dims, coords = {'time' : t, 'id' : tpnames, 'vec' : tlab}) + +tp = [tpxr] +tpda = xr.concat(tp,dim='time') +tpds = tpda.to_dataset(dim = 'vec') + +sim.ds = xr.combine_by_coords([sim.ds, tpds]) +swio.swiftest_xr2infile(sim.ds, sim.param) + +sim.param['PL_IN'] = swiftest_pl +sim.param['TP_IN'] = tpin +sim.param['CB_IN'] = swiftest_cb +sim.param['BIN_OUT'] = swiftest_bin +sim.param['ENC_OUT'] = swiftest_enc +sim.save(swiftest_input) + +sim.param['PL_IN'] = swifter_pl +sim.param['TP_IN'] = tpin +sim.param['BIN_OUT'] = swifter_bin +sim.param['ENC_OUT'] = swifter_enc +sim.save(swifter_input, codename="Swifter") diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/param.swifter.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/param.swifter.in new file mode 100644 index 000000000..aa33eeaa4 --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/param.swifter.in @@ -0,0 +1,26 @@ +! VERSION Swifter parameter file converted from Swiftest +T0 0.0 +TSTOP 3652.5 +DT 1.0 +ISTEP_OUT 11 +ISTEP_DUMP 1 +OUT_FORM XV +OUT_TYPE REAL8 +OUT_STAT UNKNOWN +IN_TYPE ASCII +PL_IN pl.swifter.in +TP_IN tp.in +BIN_OUT bin.swifter.dat +ENC_OUT enc.swifter.dat +CHK_QMIN 0.004650467260962157 +CHK_RMIN 0.004650467260962157 +CHK_RMAX 1000.0 +CHK_EJECT 1000.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 1000.0 +EXTRA_FORCE NO +BIG_DISCARD NO +CHK_CLOSE YES +J2 4.7535806948127355e-12 +J4 -2.2473967953572827e-18 +RHILL_PRESENT YES diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/param.swiftest.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/param.swiftest.in new file mode 100644 index 000000000..6504c9637 --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/param.swiftest.in @@ -0,0 +1,35 @@ +! VERSION Swiftest parameter input +T0 0.0 +TSTOP 3652.5 +DT 1.0 +ISTEP_OUT 11 +ISTEP_DUMP 1 +OUT_FORM XV +OUT_TYPE REAL8 +OUT_STAT UNKNOWN +IN_TYPE ASCII +PL_IN pl.swiftest.in +TP_IN tp.in +CB_IN cb.swiftest.in +BIN_OUT bin.swiftest.dat +ENC_OUT enc.swiftest.dat +CHK_QMIN 0.004650467260962157 +CHK_RMIN 0.004650467260962157 +CHK_RMAX 1000.0 +CHK_EJECT 1000.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 1000.0 +MU2KG 1.988409870698051e+30 +TU2S 86400 +DU2M 149597870700.0 +EXTRA_FORCE NO +BIG_DISCARD NO +CHK_CLOSE YES +FRAGMENTATION NO +ROTATION NO +TIDES NO +ENERGY NO +GR NO +YARKOVSKY NO +YORP NO +MTINY 0.0 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.in new file mode 100644 index 000000000..bd980fc4b --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.in @@ -0,0 +1,33 @@ +8 +1 4.9125474498983623693e-11 +1.6306381826061645943e-05 +0.33206272695596028566 0.07436707001147663254 -0.02438290851908785084 +-0.0115920916602103591525 0.028710618792657981169 0.0034094833969203438596 +2 7.243452483873646905e-10 +4.0453784346544178454e-05 +-0.7188115337296047125 -0.0118554711069603201795 0.041316403191083782287 +0.00021427347881133320621 -0.020313576971905909774 -0.00029114855617710840843 +3 8.9970113821660187435e-10 +4.25875607065040958e-05 +0.35677088372527121507 -0.95189300879814897627 4.4027442504036787155e-05 +0.015830039028334789986 0.0059737936889703449964 -3.3484113013969089573e-07 +4 9.549535102761465607e-11 +2.265740805092889601e-05 +-1.5233712071242269115 0.6723825347339112968 0.051459143378398922164 +-0.0051275613251079554117 -0.011607719813367209372 -0.000117479966462153095864 +5 2.825345908631354893e-07 +0.00046732617030490929307 +4.049944927347420176 -2.9910878677758190314 -0.078187280837353656526 +0.0043972077687938898594 0.006432188574295680597 -0.00012509257442073270106 +6 8.459715183006415395e-08 +0.00038925687730393611812 +6.298929503477405767 -7.706413024510769816 -0.11669919842191249504 +0.0040140666547768266703 0.0035242303011843410798 -0.00022097170940726839814 +7 1.2920249163736673626e-08 +0.00016953449859497231466 +14.856082147529010129 13.007589275314199284 -0.14417795763685259391 +-0.0026158276515510360365 0.0027821364817078499815 4.40781085949555924e-05 +8 1.5243589003230834323e-08 +0.000164587904124493665 +29.55744967800954015 -4.629377558152945049 -0.58590957207831262377 +0.00046987400245862169295 0.0031274056019462009859 -7.51415892482447254e-05 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swifter.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swifter.in new file mode 100644 index 000000000..701e9a14f --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swifter.in @@ -0,0 +1,36 @@ +9 +0 0.00029591220819207775568 +0.0 0.0 0.0 +0.0 0.0 0.0 +1 4.9125474498983623693e-11 0.0014751243077781048702 +1.6306381826061645943e-05 +0.33206272695596028566 0.07436707001147663254 -0.02438290851908785084 +-0.0115920916602103591525 0.028710618792657981169 0.0034094833969203438596 +2 7.243452483873646905e-10 0.006759104275397271956 +4.0453784346544178454e-05 +-0.7188115337296047125 -0.0118554711069603201795 0.041316403191083782287 +0.00021427347881133320621 -0.020313576971905909774 -0.00029114855617710840843 +3 8.9970113821660187435e-10 0.010044787321379672528 +4.25875607065040958e-05 +0.35677088372527121507 -0.95189300879814897627 4.4027442504036787155e-05 +0.015830039028334789986 0.0059737936889703449964 -3.3484113013969089573e-07 +4 9.549535102761465607e-11 0.007246743835971885302 +2.265740805092889601e-05 +-1.5233712071242269115 0.6723825347339112968 0.051459143378398922164 +-0.0051275613251079554117 -0.011607719813367209372 -0.000117479966462153095864 +5 2.825345908631354893e-07 0.35527126534549128905 +0.00046732617030490929307 +4.049944927347420176 -2.9910878677758190314 -0.078187280837353656526 +0.0043972077687938898594 0.006432188574295680597 -0.00012509257442073270106 +6 8.459715183006415395e-08 0.4376527512949726007 +0.00038925687730393611812 +6.298929503477405767 -7.706413024510769816 -0.11669919842191249504 +0.0040140666547768266703 0.0035242303011843410798 -0.00022097170940726839814 +7 1.2920249163736673626e-08 0.4695362423191493196 +0.00016953449859497231466 +14.856082147529010129 13.007589275314199284 -0.14417795763685259391 +-0.0026158276515510360365 0.0027821364817078499815 4.40781085949555924e-05 +8 1.5243589003230834323e-08 0.7812870996943599397 +0.000164587904124493665 +29.55744967800954015 -4.629377558152945049 -0.58590957207831262377 +0.00046987400245862169295 0.0031274056019462009859 -7.51415892482447254e-05 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in new file mode 100644 index 000000000..bd980fc4b --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in @@ -0,0 +1,33 @@ +8 +1 4.9125474498983623693e-11 +1.6306381826061645943e-05 +0.33206272695596028566 0.07436707001147663254 -0.02438290851908785084 +-0.0115920916602103591525 0.028710618792657981169 0.0034094833969203438596 +2 7.243452483873646905e-10 +4.0453784346544178454e-05 +-0.7188115337296047125 -0.0118554711069603201795 0.041316403191083782287 +0.00021427347881133320621 -0.020313576971905909774 -0.00029114855617710840843 +3 8.9970113821660187435e-10 +4.25875607065040958e-05 +0.35677088372527121507 -0.95189300879814897627 4.4027442504036787155e-05 +0.015830039028334789986 0.0059737936889703449964 -3.3484113013969089573e-07 +4 9.549535102761465607e-11 +2.265740805092889601e-05 +-1.5233712071242269115 0.6723825347339112968 0.051459143378398922164 +-0.0051275613251079554117 -0.011607719813367209372 -0.000117479966462153095864 +5 2.825345908631354893e-07 +0.00046732617030490929307 +4.049944927347420176 -2.9910878677758190314 -0.078187280837353656526 +0.0043972077687938898594 0.006432188574295680597 -0.00012509257442073270106 +6 8.459715183006415395e-08 +0.00038925687730393611812 +6.298929503477405767 -7.706413024510769816 -0.11669919842191249504 +0.0040140666547768266703 0.0035242303011843410798 -0.00022097170940726839814 +7 1.2920249163736673626e-08 +0.00016953449859497231466 +14.856082147529010129 13.007589275314199284 -0.14417795763685259391 +-0.0026158276515510360365 0.0027821364817078499815 4.40781085949555924e-05 +8 1.5243589003230834323e-08 +0.000164587904124493665 +29.55744967800954015 -4.629377558152945049 -0.58590957207831262377 +0.00046987400245862169295 0.0031274056019462009859 -7.51415892482447254e-05 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb b/examples/symba_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb new file mode 100644 index 000000000..d0d223ce7 --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb @@ -0,0 +1,753 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import swiftest\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swifter file param.swifter.in\n", + "Reading in time 3.652e+03\n", + "Creating Dataset\n", + "Successfully converted 333 output frames.\n", + "Swifter simulation data stored as xarray DataSet .ds\n" + ] + } + ], + "source": [ + "inparfile = 'param.swifter.in'\n", + "swiftersim = swiftest.Simulation(param_file=inparfile, codename=\"Swifter\")\n", + "swiftersim.bin2xr()\n", + "swifterdat = swiftersim.ds" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file param.swiftest.in\n", + "Reading in time 3.652e+03\n", + "Creating Dataset\n", + "Successfully converted 333 output frames.\n", + "Swiftest simulation data stored as xarray DataSet .ds\n" + ] + } + ], + "source": [ + "inparfile = 'param.swiftest.in'\n", + "swiftestsim = swiftest.Simulation(param_file=inparfile)\n", + "swiftestsim.bin2xr()\n", + "swiftestdat = swiftestsim.ds" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "swiftdiff = swiftestdat - swifterdat" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "swiftdiff = swiftdiff.rename({'time' : 'time (d)'})" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "swiftdiff['rmag'] = np.sqrt(swiftdiff['px']**2 + swiftdiff['py']**2 + swiftdiff['pz']**2)\n", + "swiftdiff['vmag'] = np.sqrt(swiftdiff['vx']**2 + swiftdiff['vy']**2 + swiftdiff['vz']**2)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "plidx = swiftdiff.id.values[swiftdiff.id.values < 10]\n", + "tpidx = swiftdiff.id.values[swiftdiff.id.values > 10]" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAElCAYAAADnZln1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAy2UlEQVR4nO3deXxcdb3/8dcne5ulpW1KKW1pi2ALhZa27IssgoB4lU1lUxQter0uP+Ui6lULXgT1J4IX9YoCZfHXqijKJrILlrWFAq2l0FKgaQtd06yTZJLP749zJp2mkzSZzpp5Px+PeeTMWeZ85iT5zGe+53u+x9wdEREZ/IqyHYCIiGSGEr6ISIFQwhcRKRBK+CIiBUIJX0SkQCjhi4gUCCX8AmNmc83sznB6gpk1mVlxtuPqi5kda2Yrsh0H7DqWTB5TM3vCzD4XTl9gZg/FLTvazN4IY/mYme1pZk+aWaOZ/TTdsUluUsLPM2b2lpl9sMe8i83snwN9LXd/x92r3L0zdREOjJm5mb2vr3Xc/Sl3f3+mYupLz1h6/j6ydUzd/XfufkrcrKuAG8NY/gLMATYBNe7+jUzGJrlDCV9ympmVZDuGPLUPsKzH8395Elda6ncweCjhD0JmNtbM/mRmG81stZl9pZf1JoYVdkncdveY2RYzW2lmn49bt9jMvm1mq8JmgcVmNj5cNsXMHg63W2FmH4/bbp6Z/cLM7g+3e87M9g2XPRmu9nLY9PAJMzvezOrM7Jtm9i5wa2xe3GuON7M/h+9vs5nd2Mv7m2tmd5nZ78N9v2hm0+OWTw2bRerNbJmZ/VvcstPN7F/hdmvN7LJwfncsZnYHMAG4N4z/8gEe07lm9gczuz3czzIzm93H7/VkM3vNzLaF79nilnV/yzOzVcDkuLjmA58GLg+ff9DMiszsivD3uTmMY0SPv4tLzOwd4LFw/mfNbLmZbTWzv5vZPnH7dzP7QtiMtDX8ncfH9/lw28bwuM6MOz4J/1bN7DAzW2RmDWb2npld19uxkX5ydz3y6AG8BXywx7yLgX+G00XAYuB7QBnBP/6bwIfC5XOBO8PpiYADJeHzfwC/BCqAGcBG4KRw2X8CrwLvJ0g004GRQCWwBvgMUALMJGg6ODDcbh6wBTgsXP47YEFc7A68L+758UAU+BFQDgwJ59WFy4uBl4GfhfuuAI7p5VjNBTqAc4BS4DJgdThdCqwEvh0epxOBRuD94bbrgWPD6T2AmXHx1fX2+xjgMZ0LRIDTw/d1DfBsL+9lFNAQ917+T3icPtfzb6CXuOYB/x33/GvAs8C48Dj/Gpjf4z3cHh7jIcDHwuM1Nfw9/hfwdI/f433AcIIPwY3AqeGyc4G1wKEEfzvvI/jGsau/1WeAi8LpKuCIbP//5fsj6wHoMcBfWPCP3ATUxz1a2J7wDwfe6bHNt4Bbw+m5JEj4wHigE6iO2+4aYF44vQL4aIJ4PgE81WPer4Hvh9PzgN/GLTsdeC3ueaKE3w5U9JgXS/hHhsmkpB/Hai5xCTRMMOuBY8PHu0BR3PL5wNxw+h3gUoI2bxLFEvf7SJjw+3FM5wKPxC07AGjt5b18qsd7MaCO5BP+csIPnvD5XgQfjiVx72Fy3PK/AZf0OJYtwD5xv8dj4pb/AbginP478NUE72lXf6tPAlcCo7L9fzdYHmrSyU8fc/fhsQfw73HL9gHGhs0U9WZWT1DF7rmL1xwLbHH3xrh5bwN7h9PjgVUJttsHOLzH/i4AxsSt827cdAtBtdaXje4e6WXZeOBtd4/u4jVi1sQm3L2LIEmODR9rwnkx8e/3bIIPp7fN7B9mdmQ/9xdvV8cUdj42FZa4zXxsj/fi8c+TsA9wd9zvbDnBh1P838maHuvfELf+FoIPnb7eS+z33NffTl9/q5cA+wOvmdkLZnbGgN+l7EAnYwafNcBqd99vgNutA0aYWXVcgppA8FU89rr7AksT7O8f7n5ysgEn0NeJxTXABDMr6WfSHx+bMLMigiaMdbFlZlYUl/QnAK8DuPsLwEfNrBT4D4KKtfu1+hnrro7pQKzv8V6sl3j6aw3wWXdf2HOBmU0MJ73H+le7+++S3Ne+vczv9W/V3d8Azgt/b2cBd5nZSHdvTiIGQSdtB6PngYbwpOcQC062TjOzQ/vayN3XAE8D15hZhZkdTFBhxf7Bfwv8wMz2s8DBZjaSoN12fzO7yMxKw8ehZja1n/G+R9B2O5D3tx641swqw1iP7mP9WWZ2Vlg1fw1oI2i7fg5oJjiRWWpmxwMfARaYWZkF/dqHuXsHQdt5b90se42/H8d0IO4HDox7L19hx29RA/W/wNWxE69mVmtmH93F+t8yswPD9YeZ2bn93NdvgcvMbFb4t/O+cL99/q2a2YVmVht+INeHr5W1LsSDgRL+IONB/++PEJwgXE1wAvW3wLB+bH4eQfvtOuBugnb4h8Nl1xFUuQ8RJMCbgSFh5XoK8Mlwu3fZfsK1P+YCt4Vf6T++q5Xj3t/7CNrZ6wjOI/Tmr+HyrcBFwFnu3uHu7cC/AacRHKNfAp9y99fC7S4C3jKzBuALwIW9vP41wH+F8V+WYHlfx7Tf3H0TwcnPa4HNwH7ATtX5ANwA3AM8ZGaNBB+Ch/ex/7sJfq8LwmOylODY9Sf2PwJXA/+P4MT4X4AR/fhbPRVYZmZNYbyf7KOpT/rBwpMjIoOOmc0lOCHcW7IWKSiq8EVECoQSvohIgVCTjohIgVCFLyJSIJTwZdCwBCOJDhbWY4wekWQo4UteCZNeswWDgK01s+ssw+P5Wz+GdBbJRUr4ko+mu3sVcBJwPvD5XawvIijhSx4LL5J6CpjWc1k4tO4z4QVR683sRjMri1u+q+F8Ew4FbImHdB5lZveF+9piZk+FwwHsxMyOCseF2Rb+PCpu2RNm9gMzW2jBMMIPmdmoBK9xrpkt7jHvG2b2l4EdQSk0SviSt8zsAIJRL19KsLiTYAjhUQQjbJ7EjoPMAZxBMGTvdODjwIfC1/0YwSBeZwG1BB8q8wHc/bhw2+ke3E3q98A3CK74rSUY+OvbJBhjx4Lx5u8Hfk4wtPR1wP3hEBUx5xMMNT2aYMjgRFfv3gNM6jF8xYXAHQnWFemW8wnfzG4xsw1m1nPQrmRf78GwEruvx/xJFtyc4w0LbphR1ttrSNa9aGZbgXsJLsW/tecK7r7Y3Z9196i7v0UwZPMHeqx2rbvXu/s7wOMEl/hDMCzyNe6+PByg7YfADIu74UcPHQTDC+8TDtvwlCfu7/xh4A13vyOMaz7wGsHwAjG3uvvr7t5KMJTFjJ4v4u5twO8Jh3sIx7eZSDCukUivcj7hE4zjfWoKX+8nBOOk9PQj4GfhyH1bCQa5ktw00933cPd93f2/egxxDICZ7R82s7wbjv3yQ4JqP15vw/n2ZyjgeD8huDnIQ2b2ppld0ct6YwmGR463q+GSextK+jbg/LAZ6iLgD+EHgUivcj7hu/uTBP9w3cxs37BSXxy2l04ZwOs9SjCAU/zrGcEdj+4KZ91GcIcfyV+/Iqie93P3GoJmFut7k25rgEvj7zng7kPc/elEK7t7o7t/w90nE1TrXzezkxKsuo7gwyReUsMlu/uzBDeKOZagGUjNObJLOZ/we3ET8GV3n0XQxvnL3Xy9kUB93PjqdfRezUl+qCYY1bMpLAi+OIBtdzUU8A5DIpvZGeGQv8b2oZQTDeP7AMFQ0uebWYmZfYLgLlfJNsXcDtwIRN39n0m+hhSQvLuIw8yqgKOAP8Z1qigPl50FXJVgs7Xu/qG+XjbBPI05kd8uIygMLic4qft7gm9xu+Tud4d/ZwvCdvttwMPAH8NV5hIM6TwEmENQHNxIcNJ2K/BLd38iwetutuCuTTcQfANZCZwRDn2cjDuAH4QPkV3Ki7F0LLgDz33uPs3MaoAV7r7Xbrze8cBl7n5G+NwI7pM6xt2jFtzObu4uPiREsir8wNlAcE7jjWzHI7kv75p03L0BWB37im2B6bv5mk7QS+OccNanCW6cIZLLvgi8oGQv/ZXzFb6ZzQeOJ+hh8R7wfeAxgq/EewGlwAJ3T9SUk+j1ngKmEPR+2Axc4u5/N7PJwAJgBEETwIXq9SC5yszeImiK/Ji7J7oOQWQnOZ/wRUQkNfKuSUdERJKT0710Ro0a5RMnTsx2GCIieWPx4sWb3L020bKcTvgTJ05k0aJF2Q5DRCRvmFnPq7m7qUlHRKRAKOGLiBQIJXwRkQKR0234iXR0dFBXV0ckEsl2KL2qqKhg3LhxlJaWZjsUEZFueZfw6+rqqK6uZuLEicSNpZMz3J3NmzdTV1fHpEmTsh2OiEi3vGvSiUQijBw5MieTPYCZMXLkyJz+BiIihSnvEj6Qs8k+JtfjE5HClJcJX0RksOjo7GLB8+/Q0bnTjdtSriAT/lFHHZVw/sUXX8xdd92VcJmISDrc+/I6rvjzq/zqiVVp31dBJvynn054pzoRkYwrLgqagP/x+sa07yvveumkQlVVFU1NTbg7X/7yl3nssceYNGkSGjlURDKtuS24G+bLa+px97SeAyzICj/m7rvvZsWKFbz66qv85je/UeUvIhnXEOkAINrlrK1vTeu+CjrhP/nkk5x33nkUFxczduxYTjyxX7c8FRFJmcYw4QPUt3T0sebuK+iED+pCKSLZ1dAa7Z5uae9M674KOuEfd9xxLFiwgM7OTtavX8/jjz+e7ZBEpMDEV/gt7dE+1tx9BXnSNubMM8/kscce46CDDmL//ffnAx/4QLZDEpEC0xCJUlxkdHY5kY70VvgFmfCbmpqAoDnnxhtvzHI0IlLIGiMdjKmpYG19q5p0REQGs8ZIlNE15YDa8EVEBrWG1qDCB2hNc8LPaJOOmb0FNAKdQNTdZ2dy/yIiuaYxEmV0dWYq/Gy04Z/g7puysF8RkZzS2eU0tkUZPrSMspIiWjrS20tHTToiIlnS1BYk+OqKEoaWFae9SSfTCd+Bh8xssZnNSbSCmc0xs0VmtmjjxvQPJiQiki0NrUEf/JohpQwtLR50J22PdveZwGnAl8zsuJ4ruPtN7j7b3WfX1tZmOLz++exnP8vo0aOZNm1atkMRkTzWGAkq/JqKEoYMtgrf3deFPzcAdwOHZXL/qXLxxRfz4IMPZjsMEclzsStrK8tLGFpWkvYrbTOW8M2s0syqY9PAKcDSTO0/lY477jhGjBiR7TBEJM81hxX90LKgwh9MvXT2BO4OBysrAf6fu+9WmXzlvcv417qGVMTW7YCxNXz/Iwem9DVFRBJpaYtV+MUMLStmS3N7WveXsYTv7m8C0zO1PxGRXBer8CvLgl46dVsHT4WfcqrERSSfxdrsh5YVM6S0ZHCdtBURke2a2uJP2hYPnpO2g8l5553HkUceyYoVKxg3bhw333xztkMSkTzU0tZJcZFRXlIUJnw16eSc+fPnZzsEERkEmtujDC0rxsyoKC2mLdpFZ5dTXJSeO/GpwhcRyZKWtk4qy4K6e2hZMQCtabwJihK+iEiWNLdHGVoeJPqaIaXA9uEW0kEJX0QkS1rat1f4IyrLANLaF18JX0QkS5rbot1NOSPDhL9ZCV9EZPBpbo9SWR5U+COrgpugbG5qS9v+lPBFRLKkpa2zO+GrSScHrVmzhhNOOIGpU6dy4IEHcsMNN2Q7JBHJU83tUSrDJp2aihJKi42rH1jOp295HndP+f7UD3+ASkpK+OlPf8rMmTNpbGxk1qxZnHzyyRxwwAHZDk1E8kxLWydDw5O2ZsaIyjLea2jjvYYI4UCTKaUKf4D22msvZs6cCUB1dTVTp05l7dq1WY5KRPKNu4dt+MXd80ZUBu34U8ZUp2Wf+V3h/+0KePfV1L7mmIPgtGv7tepbb73FSy+9xOGHH57aGERk0GuLdtHldFf4ACXhFbZT9qpJyz5V4SepqamJs88+m+uvv56amvT8ckRk8GqOGws/ZmNj0ENHFX4i/azEU62jo4Ozzz6bCy64gLPOOisrMYhIfmtu2363q5ht4VW2U9NU4ed3ws8Cd+eSSy5h6tSpfP3rX892OCKSp5rDoZCr4ir8mz89m78uWcfo6vK07FNNOgO0cOFC7rjjDh577DFmzJjBjBkzeOCBB7Idlojkme03P9ledx/1vlH86JyD09JDB1ThD9gxxxyTlv6xIlJYYk068W346aYKX0QkCxJV+OmmhC8ikgXdFb4SvojI4BY7aTtUTToiIoNbrMKvKleFLyIyqLW0RykyKC/JXBpWwhcRyYLm8H626eqCmYgS/gBFIhEOO+wwpk+fzoEHHsj3v//9bIckInmoJe5+tpmifvgDVF5ezmOPPUZVVRUdHR0cc8wxnHbaaRxxxBHZDk1E8khz3P1sM0UV/gCZGVVVVUAwpk5HR0dGv5KJyODQ0lYAFb6ZFQOLgLXufsbuvNaPnv8Rr215LTWBhaaMmMI3D/tmn+t0dnYya9YsVq5cyZe+9CUNjywiA9bUFs3oRVeQnQr/q8DyLOw3ZYqLi1myZAl1dXU8//zzLF26NNshiUieaWnvzGiXTMhwhW9m44APA1cDuz3U5K4q8XQbPnw4xx9/PA8++CDTpk3Laiwikl+a26PsUzY0o/vMdIV/PXA50NXbCmY2x8wWmdmijRs3Ziyw/tq4cSP19fUAtLa28sgjjzBlypTsBiUieaelbRCftDWzM4AN7r64r/Xc/SZ3n+3us2trazMUXf+tX7+eE044gYMPPphDDz2Uk08+mTPO2K1TESJSgJoHebfMo4F/M7PTgQqgxszudPcLMxjDbjv44IN56aWXsh2GiOQxd6dlMHfLdPdvufs4d58IfBJ4LN+SvYhIMm7+52ouuvm57udt0S46uzzjFb764YuIpNnTKzfx9KrNRDuD05cNkeDetdUVpRmNIysJ392f2N0++CIi+WLN1hY6u5z12yIANEWCoZGrM9wtUxW+iEgauTt1W1sBun82tYUJv0IJX0Rk0Nja0kFLezD2/dr6MOGHFX6mL7xSwhcRSaM1W1q6p+u2BtONYYVfpQo/P3R2dnLIIYeoD76I9CnWjAOwduuOFX51eWZP2u7y48XMJvTzterdvWE348kbN9xwA1OnTqWhoWDesogkYU1Y1e+/Z9VObfiZrvD7s7fbAAf6GgPYgXnA7SmIKefV1dVx//33853vfIfrrrsu2+GISA7b1NjG0LJi3je6itfebQSgMeyWWZlrV9q6+wk955nZGHd/Nz0h9d+7P/whbctTOzxy+dQpjPn2t/tc52tf+xo//vGPaWxsTOm+RWTw2dbawbAhpdRUlNIYNuU0tkUpKymivCQ/Lrz6VEqjyCP33Xcfo0ePZtasWdkORUTyQCzhV1eUdFf2TZFoxvvgQ/Jj6XzUzFqAh919RSoDGohdVeLpsHDhQu655x4eeOABIpEIDQ0NXHjhhdx5550Zj0VEct+21g5qhpRSXVFKpKOLjs4umtqiGW+/h+Qr/LOAlcCZZvbbFMaT86655hrq6up46623WLBgASeeeKKSvYj0Kr7CB2iMRIMKPwsJP6k9uvt7wIPhQ0REetHQnfCDLpiNkQ4a26IZv+gKkqzwzewXZjYvnD4lpRHlkeOPP5777rsv22GISA7rrcKvynAffEi+SacdeDOcPjFFsYiIDCodnV00t3d299KBYKTMxraOrDTpJJvwW4BhZlYK9PfCLBGRgtLQGvTK6VnhN0ay06ST7B63AK3AL4CFqQtHRGTw2BaX8GMV/tbmdra1djCisizj8Qyowjez4WZ2K3B2OOt2YHbKoxIRGQS2Jajw39rcgjuMqi7PeDwDqvDdvd7MrgUmApuAg4E/pyEuEZG8F0v4NUNKu/vdr97UBEBtVY4n/NAlwGp3/zuwOMXxiIgMGvEVfmlxEUNKi1m9qRmA2urMN+kkk/C3Al8ws/cDLwNL3P2l1IaV2yZOnEh1dTXFxcWUlJSwaNGibIckIjmoobvCD1JtdUUJb24MEv6ofKjw3f0aM3sUeB2YARwHFFTCB3j88ccZNWpUtsMQkRwWu9FJ7IRtdUUJGxrbgDxJ+GZ2FVAMLCGo7p9IcUwiIoNCUyRKSZFRXhL0j6kZEiT+oWXFVOZDt0x3/56Z7QkcApxtZvu6++dTH9quPfWH19m0pimlrzlqfBXHfnz/PtcxM0455RTMjEsvvZQ5c+akNAYRGRxig6SZBbcTmT5uOC+9U5+VZA/J98O/FPi1uxfkWDoLFy5k7NixbNiwgZNPPpkpU6Zw3HHHZTssEckxTT3GzDn5gD2Z9/RbbAybdTIt2YR/C/BFM6sEfufuS1IXUv/tqhJPl7FjxwIwevRozjzzTJ5//nklfBHZSVOPK2oPmzQCoLuJJ9OS3etXCD4sSoCfpy6c3Nfc3Nx9p6vm5mYeeughpk2bluWoRCQX9azwS4uL+P2cI7j/K8dkJZ5kK/xVwH7AX939/6Qwnpz33nvvceaZZwIQjUY5//zzOfXUU7MclYjkoua2KHv0GELh8MkjsxRN8gl/GbAGuMTMfuLuh6Ywppw2efJkXn755WyHISJ5oLEtyrgRQ7MdRrdkE/7+wEbgJoILsXbJzCqAJ4HycL93ufv3k9y/iEjOy9a9a3uTbBv+FIKLrS4D+tsnsQ040d2nE1ywdaqZHZHk/kVEcl7PNvxsSzbhDwe+CVwORPqzgQdineZLw4cnuX8RkZzW2eW0tHdmrc99Iskm/KsITtiuALr6u5GZFZvZEmAD8LC7P5fk/kVEclpzezCsQjbubNWbfiX8MFGvN7PPAbh7nbs/Ek5f0d+duXunu88AxgGHmdlO/RnNbI6ZLTKzRRs3buzvS4uI5JSmSJDw865Jx907gaXAvqnYqbvXA08AO/VndPeb3H22u8+ura1Nxe5ERDKuORw4LV+bdIYCl4fV9z3h46/93djMas1seDg9BPgg8NqAos0R9fX1nHPOOUyZMoWpU6fyzDPPZDskEckxsZEyq3KoSWcgkRwZ/pwZPmBgJ133Am4zs2KCD5o/uPt9A9g+Z3z1q1/l1FNP5a677qK9vZ2WlpZshyQiOSbWpJNL3TIHEsmk3dmRu79CMMJmXmtoaODJJ59k3rx5AJSVlVFWlvk714hIbmvK5wrf3d9OZyDJeHzeTWx4+82UvubofSZzwsW9X1rw5ptvUltby2c+8xlefvllZs2axQ033EBlZWVK4xCR/BZL+JVluZPwszNkWx6LRqO8+OKLfPGLX+Sll16isrKSa6+9NtthiUiO6W7SyccKPxf1VYmny7hx4xg3bhyHH344AOecc44SvojspCnPe+kAYGYfSUcg+WLMmDGMHz+eFStWAPDoo49ywAEHZDkqEck1zW1RykuKKC3OnYaUZD56rgbuTXUg+eR//ud/uOCCC2hvb2fy5Mnceuut2Q5JRHJMY1s0p5pzILmEbymPIs/MmDGDRYsWZTsMEclhPe92lQuS+a6hAc9ERHahuS2aU+33oF46IiJp0ZhjQyODEr6ISFo0RXKvDT+ZhP9eyqMQERlkmgZDk467n5yOQEREBpNmNemIiBSGxrZoTo2jA0r4A7ZixQpmzJjR/aipqeH666/PdlgikkPao120R7uoyqFxdCDJoRXM7Ovufl04/f7wVocF4f3vfz9LliwBoLOzk7333pszzzwzu0GJSE5pzsGRMmGACT+8gcnPgClmFgFeAS4BPpP60HLfo48+yr777ss+++yT7VBEJIvqtrZQVlLE6OoKIG5o5Bxrwx9QNOGtCT9jZh8G3gVOAf6chrj6pf7eVbSva07pa5aNrWT4R/p3J8cFCxZw3nnnpXT/IpJ/PnnTs9RtbeUPlx7JYZNG0JiD97OF5NvwP0DQPfMIoCB77bS3t3PPPfdw7rnnZjsUEckid6duaysAf1pcB0Bz+yBo0okzHPgmcDlBk05W9LcST4e//e1vzJw5kz333DNrMYhI9m1ubu+e3toSTDdGOoDcq/CTjeYqYIq7rzCzrlQGlC/mz5+v5hwR4d1tke7pLWHyb2gNKvyaIaVZiak3STXpuHuduz8STl+R2pByX0tLCw8//DBnnXVWtkMRkSxbVx8050wYMZQtPSr8mopBkPDN7BdmNi+cPiWlEeWBoUOHsnnzZoYNG5btUEQky9aHFf6BY2vYGqvwc/D2hpD8Sdt2IHb38BNTFIuISN5Zt62VsuIi3je6ivrWDjq7nIbWDspKiqgoLc52eDtINuG3AMPMrBSYkMJ4RETyyrvbIowZVsHIyjLcYVtrBw2RaM4150DyJ223AK3AL4CFqQtHRCS/xBL+HpVlQHDitiHSQU2ONefAACt8MxtuZrcCZ4ezbgdmpzwqEZE8Ud/SwR5DSxkRJvytLe00tHZQnWM9dCCJK23N7FpgIrAJOJgsXmkrIpJt21o7GDZke8Lf3NROYySakxV+MhFdAqx2978Di1Mcj4hIXumZ8Le2BE06ew8fkuXIdpbMSdutwBfM7Hoz+4yZHZLqoHLdz372Mw488ECmTZvGeeedRyQS2fVGIjLotEe7aO3oZNiQUvYYGteG3xqlZkjuVfjJ3PHqGuDzwFxgNXBcf7Yzs/Fm9riZLTezZWb21YHuOxesXbuWn//85yxatIilS5fS2dnJggULsh2WiGRBQ+wCqyGlVJQWU11ewqamNhojHYOjl46ZXQUUA0uAJe7+RD83jQLfcPcXzawaWGxmD7v7vwYaQ7ZFo1FaW1spLS2lpaWFsWPHZjskEcmCba1Bwh8WnqCtrS6nbmsrbdGunLvoCpJI+O7+PTP7HsG3g7PNbF93/3w/tlsPrA+nG81sObA3kHTC/9vf/sa7776b7OYJjRkzhtNOO63X5XvvvTeXXXYZEyZMYMiQIZxyyimcckrBXWwsImxP+LExc0ZVl7NqY9MO83JJshde3QJMBUYCvxzoxmY2ETgEeC7BsjlmtsjMFm3cuDHJ8NJn69at/PWvf2X16tWsW7eO5uZm7rzzzmyHJSJZkKjCf3NjcI+OQdGkE/oKwfAKJcAN9LMdH8DMqoA/AV9z94aey939JuAmgNmzZ3tfr9VXJZ4ujzzyCJMmTaK2thaAs846i6effpoLL7ww47GISHY19Ez4VeXdy0ZXlyfcJpuSrfBXARXAX919IMm+lCDZ/87d87L//oQJE3j22WdpaWnB3Xn00UeZOnVqtsMSkSzobtKp2F7hx+y3Z3VWYupLsgl/GfAYcImZvdCfDczMgJuB5bEboOejww8/nHPOOYeZM2dy0EEH0dXVxZw5c7IdlohkQc8KP76qH1VVlpWY+pJsk86+BP3xbwp/9sfRwEXAq2a2JJz3bXd/IMkYsubKK6/kyiuvzHYYIpJl21o7GFJaTFlJUDvHKvyy4iKCGje3JJvw17j7Y2a2F7ChPxu4+z+B3DsCIiJJil1lGzMqbMPfozL3TthC8k06p5rZOOB/gZ+lMB4RkbzR84raPWsqAPjE7PHZCqlPqbiJ+edSFk0/uXtOfl2Kce+zc5GIDBJNbdEdblReW13Oc98+aYfeOrkk2Qr/KoIeOiuAzhTGs0sVFRVs3rw5Z5Oqu7N582YqKiqyHYqIpFlze5TK8h3r5j1rKigqys2CtF8VvpkVA3XAd939t+5eFz7P+E3Mx40bR11dHbl4UVZMRUUF48aNy3YYIpJmzW1R9qzOn+KuXwnf3TvNbClB75ysKi0tZdKkSdkOQ0SE5rbOnSr8XDaQSIcCl5vZycC6cJ67+0dTH5aISO4L2vBz60blfRlIwj8y/DkzfADkZkO6iEgGtLRHGTpIK3y1o4iIhNqinXR0+g69dHLdLiM1swnhZMJqPm55faLB0EREBqPmtqCDYmXZ4GrSuY0g2ffVz8iBecDtKYhJRCTnNbdFAQbXSVt3PyETgYiI5JPm9vxL+MleeCUiUtDyscJXwhcRSUJT2IafT90ylfBFRJKgCl9EpEB0J/wyJXwRkUFNFb6ISIFobg/74asNX0RkcGtqi1JSZJQV508azZ9IRURySHNbMBZ+Lt+MqSclfBGRJDS3debVODqghC8ikpSgws+f9ntQwhcRSUpze5ShedQlE5TwRUSS0vMG5vlACV9EJAktbZ1q0hERGWxef6+RiVfcz4vvbO2e1xT20sknSvgiIrtwyz9XA/DY8g3d85rbo3k1rAIo4YuI7GD1pmZeqavvfu7uLFy1CYBo1/Yb/zWrwhcRyW8X/OZZ/u3GhWxpbgdg3bYIa7a0ArB+W/CzPdoV3s9WbfgJmdktZrbBzJZmap8iIgPVFu0C4NdPrgJgQ0Oke9n6+mA6HwdOg8xW+POAUzO4PxGRAdtreAUAi98KTtBubGwDYN/aStbWBxV+Ux4OjQwZTPju/iSwJVP7ExFJRlMkSObbWjsA2NgUJPzp44fzXkOEzi7Py/vZQg624ZvZHDNbZGaLNm7cmO1wRKTAxKr3WMLf1Bi05R+09zCiXc7Gxjaa2/JvaGTIwYTv7je5+2x3n11bW5vtcESkwDTuVOFH2GNoKRNGDAWCE7exNvx8u9I2v6IVEUmjjs4u2qJdlJcU0RbtItLRycbGNmqryxk+tAyA+tYOIuHNTzSWjohInopV7nvvMQSAhkhHd8KvqQiSe1Mk2t3sk28Vfia7Zc4HngHeb2Z1ZnZJpvYtItIfseacvYeHCb+1g01N7dRWlVNdUdq9zvZumfnVhp+xjyd3Py9T+xIRSUasch8XVvjbWoMKf1RVOVVhhd8Y6ei+4jbfeunkV7QiImkUS/hjhwUJf219hNaOTmqry6ksK6bIgnWiXU5JkVFekl+t4kr4IiKhWB/8sWGTzqoNTQDUVpdjZlSVl9AYiRLt6qK6Ir/uZwtK+CIi3Zp6nLRdtXF7wgeoriilIdJBZ5dTM6Q0O0HuBiV8EZFQd8IPK/yVG3om/LDC7wwq/HyTfxGLiKRJrEln+NBSKsuKuyv8UVXbE35TJEpHZxfV5flX4efXGQcRkTRqjBsUrWZIKR2dTnGRsUd40VV1RSmNbR00RqLUDMm/elkJX0Qk1BSJBr1xiozxewRDKYysLKO4KDg5Gztp2xDp6O6Xn0+U8EVEQs1t0e7+9jP32QNgh7b6WJNOYyRKTR4m/Pz7TiIikiZNbdHu4RJmhwl//bbtN0CpriilvjXopZOPJ21V4YuIhBrbolSFlXuswm8JB0qDoMLvDK+yVbdMEZE81hTpoDqs8EdUlnHpByZz1L6jupf3bN7JN/kXsYhImjS1Rbv73AN867SpOywfFlfV52Mbvpp0RERCzW2dVPXWv76ri5m+nIMtuLl5jSp8EZH81RjpSNxU0xGBP36a8a8/yD3lcHv0ZGoqjt69nbVsgYphUJS5IZZV4YuIAO6+Qy+dHfzjWnj9QTj5B9wcPY1PlTzMmNV/Sm5H7c3wu3Phx5Pgf2bB2sW7F/gAKOGLiACtHZ10Od398LttXgULfw4zLoSjv8LLB/wnz3VNYcTCqyDSMPAd/XkOrHwEjv4qdHXC/POh8b3UvIldUMIXEWH7ODo73dRk4fVQXAof/D4A//fjhzD6nP9LUaQeFt86sJ288TC8dh+c9D04+So4fwG0bIYnfrj7b6AflPBFRNg+UmZ1fMJv3gRL5sOMC6BqNABlJUVMOvhYmPQBePZX0Bnt3w7c4dGrYMRkOOJLwbw9D4TZn4UX74Atb6by7SSkhC8iAolvTP7yAujqgMM+v/MGh18KjeuDtv3+WPM8vPsKHPVlKCnbPv/Yr4MZvHDzbkTfP0r4IiJsb9LpbsN3hxdvh3GHwuipO2+w34egei9YPK9/O3jhN1BeAwd9fMf51WNgyodhye+gozX5N9APSvgiImwfGrm7wq97ATatgEMuSrxBcUmwbOUjUP9O3y/etAGW/QVmnA/lVTsvn30JtG4N1kkjJXwREeIq/FjCf/F2KK2EaWf1vtHM8MPgxTv6fvEXbwuahg79XOLlk46DkfvBovQ26+jCKxERoLk9rkknsg2W/hmmnQnl1QnXj0QjPLjpRf4xaSpvr55P11+eZ2zV3hw77lhOn3Q6w8qHBSt2dcKieTD5eBi1X+KdmwUnb//+LVj/Cux1cOrfIKrwRUQAaIyv8F/5A3Q0B0m4h86uTm5bdhsn/fEkvrvwu/yrrIS921qZVDSUuqY6fvjcDznxDyfys8U/o7mjGV7/OzTU9V7dx0z/JJRU9P+cQBJU4YuIAFub2ykvKaK82GDRrbDXdBg7c4d13ml4h/9a+F+8tOEljtn7GD477bPMHjUdu/4gaOyEC+5h+ebl3Ln8Tm5Zegv3rbqPuW2lHFs9FvY/re8Aho6AAz4WfNh88PvBsAsppgpfRAR4e0sL+4wcitW9ABuWBdW9Bbc27PIuFry2gHPuPYeVW1dy9TFX88uTfsmhYw7FSsrgkAth5cNQv4apI6dy9TFXc+fpd1JTMoR/71rLf0+cSktX+66DOOIL0N4Ii25Jy3tUhS8iAry9uZkJIyrh+V9BWTVMOweA9U3r+d7T3+PZ9c9y1NijuPKoKxlTOWbHjWdeBE/9NEjU4RW502uns6BkEj9vWMbtvMFz932CHx5+FVN8NNHNm4lu2kx08yY6t2ylK9KKt7XjkQj+2lSKlv0vex7+RSitSOl7VMIXkYLX1eW8s6WFj41vgWV/hiP+HS+r5C9v3M2PX/gxnd7Jd4/4Lufufy4WVv072GMiHHgmPPsrfPpFtDca7a8spP33d3PR0Gmc2lZF/arXKKo/n5WeIAAzrKKCorIyKIKSYeOgpDzBirsnownfzE4FbgCKgd+6+7WZ3L+ISCIbGtuIdHRx+qZ5UFLB8gM/zLUPXsyLG15k1p6z+MHRP2B89fju9bsiETrq6mh/Zw0da96h/e23aV/ZTvu/auj43YehO6lXU1S9maqJNVQfcRLPFr/N011vsMdeE7ngqC8xadIMikeMwMrLMTOaO5qZt2we/9r8L24EEny07JaMJXwzKwZ+AZwM1AEvmNk97v6vTMUgIpLI25ubObFoERvr/8GPJs9g2V/mML69mp/ueT6zt02m87Z7WFdXR/uaNXSsWUN0w4Ydti+qrqZsn30Ycshshm19irLKDsqGF1F68a8pnnFG97eCicDwtx7kB8/8gLtXXsGHoh/i4mkXM7prNPNfm8+CFQtobG/k1ImnEumMMKRkSErfp7kn+n6RemZ2JDDX3T8UPv8WgLtf09s2s2fP9kWLFg14X9d/90d0kpn3JSKSamVexJf/+/KktjWzxe4+O9GyTDbp7A2siXteBxzecyUzmwPMAZgwYUJSOxrSUYSn+ruQiAx+ZmEzivVoT0lfQnGcoA3IiO292NLTgTKTCT/REdupDHf3m4CbIKjwk9nRpdf+ZzKbiYgMapnsh18HjI97Pg5Yl8H9i4gUtEwm/BeA/cxskpmVAZ8E7sng/kVEClrGmnTcPWpm/wH8naBb5i3uvixT+xcRKXQZ7Yfv7g8AD2RynyIiEtBYOiIiBUIJX0SkQCjhi4gUCCV8EZECkbGhFZJhZhuBt5PcfBSwKYXhpEM+xAiKM5XyIUZQnKmU6Rj3cffaRAtyOuHvDjNb1Nt4ErkiH2IExZlK+RAjKM5UyqUY1aQjIlIglPBFRArEYE74N2U7gH7IhxhBcaZSPsQIijOVcibGQduGLyIiOxrMFb6IiMRRwhcRKRCDLuGb2almtsLMVprZFTkQz1tm9qqZLTGzReG8EWb2sJm9Ef7cI279b4WxrzCzD6UxrlvMbIOZLY2bN+C4zGxW+P5WmtnPLXbzzvTFONfM1obHc4mZnZ7lGMeb2eNmttzMlpnZV8P5uXYse4sz145nhZk9b2Yvh3FeGc7PmePZR4w5dSwTcvdB8yAYdnkVMBkoA14GDshyTG8Bo3rM+zFwRTh9BfCjcPqAMOZyYFL4XorTFNdxwExg6e7EBTwPHElwR7O/AaelOca5wGUJ1s1WjHsBM8PpauD1MJZcO5a9xZlrx9OAqnC6FHgOOCKXjmcfMebUsUz0GGwV/mHASnd/093bgQXAR7McUyIfBW4Lp28DPhY3f4G7t7n7amAlwXtKOXd/EtiyO3GZ2V5Ajbs/48Ff7+1x26Qrxt5kK8b17v5iON0ILCe4f3OuHcve4uxNtuJ0d28Kn5aGDyeHjmcfMfYmK8cykcGW8BPdKL2vP+pMcOAhM1tswQ3aAfZ09/UQ/CMCo8P52Y5/oHHtHU73nJ9u/2Fmr4RNPrGv9lmP0cwmAocQVHw5eyx7xAk5djzNrNjMlgAbgIfdPeeOZy8xQo4dy54GW8Lv143SM+xod58JnAZ8ycyO62PdXIwfeo8rG/H+CtgXmAGsB34azs9qjGZWBfwJ+Jq7N/S1ai/xZCvOnDue7t7p7jMI7nt9mJlN62P1rMTZS4w5dyx7GmwJP+dulO7u68KfG4C7CZpo3gu/zhH+3BCunu34BxpXXTjdc37auPt74T9bF/Abtjd5ZS1GMyslSKK/c/c/h7Nz7lgmijMXj2eMu9cDTwCnkoPHs2eMuXwsYwZbws+pG6WbWaWZVcemgVOApWFMnw5X+zTw13D6HuCTZlZuZpOA/QhO6mTKgOIKv1o3mtkRYe+CT8Vtkxaxf/rQmQTHM2sxhq95M7Dc3a+LW5RTx7K3OHPweNaa2fBwegjwQeA1cuh49hZjrh3LhNJ5RjgbD+B0gh4Iq4DvZDmWyQRn518GlsXiAUYCjwJvhD9HxG3znTD2FaTxjD0wn+BrZwdBpXFJMnEBswn+sFcBNxJevZ3GGO8AXgVeIfhH2ivLMR5D8DX8FWBJ+Dg9B49lb3Hm2vE8GHgpjGcp8L1k/2fSFWcfMebUsUz00NAKIiIFYrA16YiISC+U8EVECoQSvohIgVDCFxEpEEr4IiIFQglfCoKZDTezf497PtbM7krTvj5mZt/rZVlT+LPWzB5Mx/5FeqOEL4ViONCd8N19nbufk6Z9XQ78sq8V3H0jsN7Mjk5TDCI7UcKXQnEtsG84TvlPzGyihePsm9nFZvYXM7vXzFab2X+Y2dfN7CUze9bMRoTr7WtmD4YD4T1lZlN67sTM9gfa3H1T+HySmT1jZi+Y2Q96rP4X4IK0vmuROEr4UiiuAFa5+wx3/88Ey6cB5xOMf3I10OLuhwDPEFzyDsHNqL/s7rOAy0hcxR8NvBj3/AbgV+5+KPBuj3UXAccm+X5EBqwk2wGI5IjHPRgnvtHMtgH3hvNfBQ4OR5k8Cvhj3E2JyhO8zl7AxrjnRwNnh9N3AD+KW7YBGJua8EV2TQlfJNAWN90V97yL4P+kCKj3YEjcvrQCw3rM6238kopwfZGMUJOOFIpGglv7JcWDseNXm9m5EIw+aWbTE6y6HHhf3POFBKO2ws7t9fuzfURFkbRTwpeC4O6bgYVmttTMfpLky1wAXGJmsdFPE90+80ngENve7vNVghvfvMDOlf8JwP1JxiIyYBotUyTFzOwG4F53f2QX6z0JfNTdt2YmMil0qvBFUu+HwNC+VjCzWuA6JXvJJFX4IiIFQhW+iEiBUMIXESkQSvgiIgVCCV9EpEAo4YuIFIj/D1muhiVRzFDnAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "swiftdiff['rmag'].sel(id=plidx).plot.line(ax=ax, x=\"time (d)\")\n", + "ax.set_ylabel(\"$|\\mathbf{r}_{swiftest} - \\mathbf{r}_{swifter}|$\")\n", + "ax.set_title(\"Heliocentric position differences \\n Planets only\")\n", + "fig.savefig(\"rmvs_swifter_comparison-mars_ejecta-planets-rmag.png\", facecolor='white', transparent=False, dpi=300)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAElCAYAAADnZln1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxL0lEQVR4nO3deZxcdZnv8c/TeyfdSUjSITtZWBISICYRUDECCoKDF1lEA4yACI7jOHoVt5m5KnoV1OuCg8sgm4BDZmQGRTZZArIPBEgggJElgXTWztL03tVV9dw/zqnu6k6vldr7+3696tWnTp3lqdPdT/3qd37nOebuiIhI8SvJdQAiIpIdSvgiIqOEEr6IyCihhC8iMkoo4YuIjBJK+CIio4QS/ihkZt8ys1vC6dlm1mJmpbmOazBm9l4z25DlfbqZHbyf23jJzI5PT0T7bHvA36OZHWhmj5hZs5n9yAI3mNleM3s6E/FI/lPCL0BmtsnMPtBn3oVm9thIt+Xub7l7jbvH0hfhyAwnsbr7o+5+WLZiShd3X+TuD0PvBJ2B/fT9PV4K7ALGufuXgOOAk4CZ7n50JmKQ/KeEL3nPzMpyHUMBOgh42XuurDwI2OTurSPdkI5/8VDCL1JmNt3M/svMGsxso5n94wDLzQlb2GVJ691hZnvM7DUzuyRp2VIz+yczez3sKnjWzGaFry0ws/vD9TaY2TlJ691oZj83s7vC9f7HzOaHrz0SLrYu7JL4mJkdb2b1ZvZVM9sO3JCYl7TNWWb23+H7221mVw9wDNrNbGLSvHeY2S4zKw+ff9LMXgm7Ov5kZgcNcJzGm9lN4f7eNLN/MbOSpNcvCbfTbGYvm9nScP4mM/uAmZ0C/BPwsfB9rjOzj5rZs3328yUz+/0AMcw1sz+H+7gfmNzf79HMbgQuAL4S7uvTwLXAu8Lnl4frnGZma82s0cyeMLMjk7a3KTz+LwCt4XaPDZdrDOM/Pmn5h83sO2b2eBjffWaWHN9xSetuNrMLw/mVZvb/zOwtM9thZr8ys+rwtclmdme4zh4zezT5mEsK3F2PAnsAm4AP9Jl3IfBYOF0CPAt8A6gA5gFvAB8MX/8WcEs4PQdwoCx8/mfgF0AVsARoAN4fvvZl4EXgMMCAo4BJwFhgM3ARUAYsJehOWBSudyOwBzg6fP23wKqk2B04OOn58UAU+D5QCVSH8+rD10uBdcBPwn1XAccNcKxWA5ckPf8h8Ktw+iPAa8DCMK5/AZ7oLy7gJuAPQG14zP4KXBy+9lFgC/DO8LgcDBzU93eVfNzD55XhcVmYNO954KwB3suTwI/D9VYAzYP8Hm8E/m9/fx/h86XATuCY8HheEMZamRT3WmBWePxnALuBDxH8fZ0UPq8Ll38YeB04NFz+YeDK8LXZYawrgXKCv5kl4Ws/Be4AJobH9o/AFeFrVwC/CtcpB94LWK7//wr5kfMA9Ejhlxb8M7YAjUmPNnoS/jHAW33W+TpwQzjdnXiSE0X4zx0DapPWuwK4MZzeAJzeTzwfAx7tM+/fgG+G0zcC1ya99iHgL0nP+0v4EaCqz7xEwn8XwQdR2TCO1aeA1eG0EXwwrQif30OYtMPnJeFxPCg5LoKE2AkcnrTsp4GHw+k/AZ8f5HfVb8IP5/0S+G44vQjYS5h0+yw3m+BDcGzSvH/v7/eYdMwHS/i/BL7TZx8bgPclxf3JpNe+CtzcZ/k/AReE0w8D/5L02t8D9yb97d3ez3syoBWYnzTvXcDGcPrbBB+yB/ddV4/UHvp6VLg+4u4TEg+Cf7CEg4Dp4VfhRjNrJOhOOHCIbU4H9rh7c9K8NwladxB8ILzez3oHAcf02d95wNSkZbYnTbcBNUPE0uDuHQO8Ngt4092jQ2wD4DaCrozpBK1iBx5NivuqpJj3ECShGX22MZngm9KbSfOGc1yG4zfAuWZmwN8C/+nunf0sNx3Y67374N/sZ7nhOgj4Up/f2axwPwmb+yz/0T7LHwdMS1pmoN/xQMenDhgDPJu0zXvD+RB8G3sNuM/M3jCzr438bUoynYwpTpsJWkmHjHC9rcBEM6tNSvqzCborEtudD6zvZ39/dveTUg24H4OVcd0MzDazsqGSvrs3mtl9wDkEXTe3eth8DLfzXXf/7RCx7AK6CE+EhvP6Oy5D2ec9uftTZhYh6K44N3z0ZxtwgJmNTUr6s/vb5jAl3vt3hxnvZoIW/iUDLTzEvvobGbQLaCfo+tvS98Xwb/BLBB9Mi4CHzOwZd38whRgEnbQtVk8DTeFJt2oLTrYuNrN3DraSu28GngCuMLOq8CTexQR97hCc+PuOmR1igSPNbBJwJ3Comf2tmZWHj3ea2cJhxruD4DzDSN7fNuBKMxsbxvqeQZb/d+ATwFnhdMKvgK+HySRxYvajfVf2YKjjfwLfNbNaC07sfhFIDLG8FrjMzJaFx+Vg6//k7w5gTj8nHm8Crgai7t7v0Fp3fxNYA1xuZhVmdhzw4UHe81B+DfydmR0TxjzWzP7GzGoHWP4W4MNm9sHw76nKghPpM4exr98CHzCzc8KTv5PMbIm7x8M4fmJmUwDMbIaZfTCcPi08lgY0EXQ35mz4cDFQwi9CYYL6MMFJ140ELalrgfHDWH0lQX/wVuB2gn74+8PXfkyQ+O4j+Ae8DqgOW2InAx8P19tOzwnX4fgW8Jvwa/05Qy2c9P4OBt4C6gnOIwzkDuAQYIe7r0vazu1hnKvMrIngm8upA2zjcwT9zW8AjxF8cFwfbud3wHfDec3A7wlOQvb1u/DnbjN7Lmn+zcDi8OdgziU4P7MH+CbBB0VK3H0NcAnBB81egq6TCwdZfjNwOkHXYANBq/3LDCOHuPtbBOdtvhTGvpbghD8E5wZeA54KfwcPEAwKgOB39gDB+aongV94eE2DpMZ6vt2KSC6EwxB3Akvd/dVcxyPFSy18kdz7DPCMkr1kmk7aiuSQmW0iGBn0kdxGIqOBunREREYJdemIiIwSSvhSNKyfKqLFwvrUPBJJhRK+FJQw6bVaUARsi5n92LJcy9/SUCdfJBeU8KUQHeXuNcD7Ccamp3L1p8ioo4QvBcvd/0JQF2dx39fM7GgzezK8mGubmV1tZhVJr7uZ/Z2ZvWpBaeSfh1d0Jl7vt2yy9V/OedhlfM3s3Wb2jJm9Hf58d9Jrg5YYTlpuRGWVRRKU8KVgmdnhBDVonu/n5RjwvwkKn72L4NvA3/dZ5jSCksZHEdTaSVzS/xGCK0rPJCjk9ShwK4C7rwjXPcqDO0z9B8EVpPXhsgeG6+4z/M2Cuvx3AT8jKBH8Y+CusDxFwrkEZaanEBRsu6yf93YHMLdP6YrzGfpKXRnl8j7hm9n1ZrbTzPoW7Ep1e/eGLbE7+8z/Bwtu+OH9taokrzxnZnsJaqdfC9zQdwF3f9bdn3L3qLtvIijX/L4+i13p7o3hpf8PEZSigKD08RXu/kpYnO17wJIB6uNAUFhtGkFZ5S4PbsfY33jnvwFedfebw7huBf5C75o4N7j7X929naCMxZK+Gwmraf4HQZInrAU0h6CmkciA8j7hE9T1PiWN2/shQRnavh4HPsD+lZyV7Fjq7ge4+3x3/5ewCFcvZnZo2M2yPazR8j2S7hAVGqic73DLJicMt4zvdPb9+0ouszxYTH0Nt6yySLe8T/ju/gjBP1w3M5sfttSfDftLF4xgew8SFLjqO//5sCUoxeGXBK3nQ9x9HEE3iw2+SrfNwKeT7zfg7tXu/kR/C7t7s7t/yd3nEbTWv2hm7+9n0a0EHybJksssD5u7P0Vwk5hEWWV158iQ8j7hD+Aa4HPuvoygj/MXOY5H8k8tQUXPlrBB8JkRrDtU2eRe5ZxHUMb3boIy0ueGZYI/BhxO6l0xQ5ZVFklWcBdxmFkN8G7gd0mDKirD184kuC1aX1vc/YPZiVDyxGUEDYOvEJzU/Q/gxOGs6O63h39nq8J++7eB++kpb/wtgnLO1cClBF0yVxOctN3LAGV83X23mZ0GXEXwDeQ14DR335Xie7wZ+E74EBlSQdTSMbM5wJ3uvtjMxgEb3H3aEKsNtr3jgcvc/bR+XtsELN+Pf0KRrFBZZRmpguvScfcmYGPiK7YFjhpiNZFipLLKMiJ5n/DN7FaCu90cZmb1ZnYxwQ2yLzazdcBLBHfiGe72HiX4av7+cHuJsdf/aGb1wEzgBTO7Nt3vRSRdwm+inye4BkBkWAqiS0dERPZf3rfwRUQkPfJ6lM7kyZN9zpw5uQ5DRKRgPPvss7vcva6/1/I64c+ZM4c1a9bkOgwRkYJhZgNWC1CXjojIKKGELyIySijhi4iMEnndh9+frq4u6uvr6ejoyHUoA6qqqmLmzJmUl5fnOhQRkW4Fl/Dr6+upra1lzpw5JNXSyRvuzu7du6mvr2fu3Lm5DkdEpFvBdel0dHQwadKkvEz2AGbGpEmT8vobiIiMTgWX8IG8TfYJ+R6fiIxOWU34ZrbJzF40s7VmpgH2IiJAeyTGbc/Wk+lSN7lo4Z/g7kvcfXkO9g3Au9/97n7nX3jhhdx2221ZjkZERrvv3PUyl/1uHU+9sWfohfdDQXbp7K8nnuj3TnUiIjmxtbEdgPauaEb3k+1ROk5wo2cH/s3dr+m7gJldSnAXIWbPnp2RIGpqamhpacHd+dznPsfq1auZO3duxr9OiYj0JxYPck9Jhs//ZbuF/x53XwqcCnzWzFb0XcDdr3H35e6+vK6u3/o/aXP77bezYcMGXnzxRX7961+r5S8iOZFI+GUlmU3JWU347r41/LkTuB04Opv77+uRRx5h5cqVlJaWMn36dE48cVi3PBURSatomPAzPcAvawnfzMaaWW1iGjgZWJ+t/Q9EQyhFJNcSLfyuWDyj+8lmC/9A4LHwtoRPA3e5+71Z3P8+VqxYwapVq4jFYmzbto2HHnool+GIyCiVaOFHoplN+Fk7aevubwB5dbPxM844g9WrV3PEEUdw6KGH8r73vS/XIYnIKBTvbuFnduBIwdXSSYeWlhYg6M65+uqrcxyNiIx2xdilIyIi/Ugk/IgSvohIcYvGg0SvFr6ISJGLZemkrRK+iEiORdWHLyIyOmRrlI4SvohIjmVrHL4Sfgo++clPMmXKFBYvXpzrUESkCHR0xQB16eSlCy+8kHvvzelFwiJSRDq6gkSvFn4eWrFiBRMnTsx1GCJSBKKxePf4+0y38Av6StvL//gSL29tSus2D58+jm9+eFFatykiMpCOpFZ9RCdtRUSKV3sk1j2tFv4g1BIXkUKXOGELOmkrIlLUkhO+TtrmoZUrV/Kud72LDRs2MHPmTK677rpchyQiBao9iy38gu7SyZVbb7011yGISJFI7sPXSVsRkSLWGXbjlBh0qUtHRKR4JUojj6ko00lbEZFiliiYNqaiVDdAEREpZtHkhK8uHRGR4pXo0qlWl46ISHFL7tJRPfw8s3nzZk444QQWLlzIokWLuOqqq3IdkogUsGgscdK2VOPw801ZWRk/+tGPWLp0Kc3NzSxbtoyTTjqJww8/PNehiUgB6oqrDz9vTZs2jaVLlwJQW1vLwoUL2bJlS46jEpFC1dPCL8v4KJ3CbuHf8zXY/mJ6tzn1CDj1ymEtumnTJp5//nmOOeaY9MYgIqNGYpROdRa6dNTCT1FLSwtnnXUWP/3pTxk3blyuwxGRAtWVuPCqPPMnbQu7hT/Mlni6dXV1cdZZZ3Heeedx5pln5iQGESkOyePwY3EnFndKSywj+1ILf4TcnYsvvpiFCxfyxS9+MdfhiEiBS/ThV1WUApmtmKmEP0KPP/44N998M6tXr2bJkiUsWbKEu+++O9dhiUiB6oo7ZSVGRWmQjjN54jbrXTpmVgqsAba4+2nZ3v/+Ou6443DPbD+biIwe0VicslKjoixI+JmsmJmLFv7ngVdysF8RkbwTjTvlJSWUhy38TJ64zWrCN7OZwN8A12ZzvyIi+Soac8pKLSnhF08L/6fAV4AB35GZXWpma8xsTUNDQ9YCExHJhWg8TllpSXeXTmcxdOmY2WnATnd/drDl3P0ad1/u7svr6uqyFJ2ISG50xZzyEqOi1MLnRZDwgfcA/8vMNgGrgBPN7JYs7l9EJO8EJ21LiqtLx92/7u4z3X0O8HFgtbufn639i4jko6548fbhF7yOjg6OPvpojjrqKBYtWsQ3v/nNXIckIgUsGotTXpKdPvyclFZw94eBh3Ox7/1VWVnJ6tWrqampoauri+OOO45TTz2VY489NtehiUgB2neUTpEMyywGZkZNTQ0Q1NTp6urCLDN1L0Sk+AVdOiXdV9pm8sKrgi6e9v2nv89f9vwlrdtcMHEBXz36q4MuE4vFWLZsGa+99hqf/exnVR5ZRFIWdOkY5WXFNUqnaJSWlrJ27Vrq6+t5+umnWb9+fa5DEpEC1bdLp6hq6aTTUC3xTJswYQLHH3889957L4sXL85pLCJSmLricWrKy3qKpxXDhVfFoqGhgcbGRgDa29t54IEHWLBgQW6DEpGCFY2F1TLLMn/StqBb+Lmwbds2LrjgAmKxGPF4nHPOOYfTTiu4op8ikie6snjhlRL+CB155JE8//zzuQ5DRIpENO6UlxrlRVZaQURE+ojG4pQlXXiVyZO2SvgiIjnUlRilU6KTtiIiRWHznjb+sr1pn/nReJyyEqOkxCgrMfXhi4gUuvf+4CEANl35N73mB+Pwg7Z3eWmJSiuIiBSrrvBKW4DyUstol86QLXwzmz3MbTW6+77fV0REZECxeE8Lv6KsNOddOr8BHBisQpgDNwI3pSGmghCLxVi+fDkzZszgzjvvzHU4IlKgEvXwASpy3cJ39xP6zjOzqe6+PTMhFYarrrqKhQsX0tSkLzUikrpEPXyA8rKSvByH/4m0RlFg6uvrueuuu/jUpz6V61BEpADE4z0nYt291/y4093Cz/RJ21RH6ZxuZm3A/e6+IZ0BjcT2732PzlfSWx65cuECpv7TPw26zBe+8AV+8IMf0NzcnNZ9i0hxao1Eu6e7Yk5FohRyPGjNlyeN0snHC6/OBF4DzjCza9MYT9678847mTJlCsuWLct1KCJSIJo6ehJ+RzTWPR0NW/Nl4Sidigx36aTUwnf3HcC94SNnhmqJZ8Ljjz/OHXfcwd13301HRwdNTU2cf/753HLLLVmPRUQKQ3NHV/d0Z1ccqoLp7oSfGKWT4ZO2KbXwzeznZnZjOH1yWiPKc1dccQX19fVs2rSJVatWceKJJyrZi8igmtp7WvidSS38ni6d5D78PEv4QAR4I5w+MU2xiIgUpeQWfkdXT0Lv6dJJ7sPPv5O2bcB4MysHhnthVtE5/vjjOf7443MdhojkuabkLp3kFn7Ymu8eh19Wkpc3Md8DtAM/Bx5PXzgiIsWnOfmkbXILPxyuWd594VUedemY2QQzuwE4K5x1E7A87VGJiBSRpvb+W/jRRAu/u0vH8ucm5u7eaGZXAnOAXcCRwH9nIC4RkaKR3MLvTGrhJy6y6nXSNs+6dC4GNrr7n4Bn0xyPiEjRSR6Hn9zCT0xXlpcCQWmFfDtpuxf4OzM7DFgHrHV33eRVRGQAyUk+uQ+/M2zNV5YlxuHn2YVX7n6FmT0I/BVYAqwAlPBFRAYQicapLCuhMxrvk/zDFn5Z0MLPuyttzezbQCmwlqB1/3CaY8p7c+bMoba2ltLSUsrKylizZk2uQxKRPBaJxqmtKqezpXPQFn7Ob4DSl7t/w8y+QTDC5ywzm+/ul6Q/tPz20EMPMXny5FyHISIFIBKLM666jF0tnX368IPkXlUeJPyxlWVE405HV4yqsF8/nVK90vZ6YCEwCfjFcFYwsyoze9rM1pnZS2Z2eYr7FhEpKIkWPvQepdPZp0tnXLhM8jDOdEr1wqt/JCivUAZcRdCPP5RO4ER3bwmv0H3MzO5x96dSjIFH//Ov7Nrckurq/Zo8q4b3nnPooMuYGSeffDJmxqc//WkuvfTStMYgIsUlEo0zpryU0hLrVS2zu0snbOGPrw4S/tvtXUwZV5X2OFJN+K8DhwB/cPf/PZwVPKj6n8jO5eEjc+OPMujxxx9n+vTp7Ny5k5NOOokFCxawYsVwPvNEZDSKxOKMrSyjqqykVwu/70nb5ISfCakm/JeAzcDFZvZDd3/ncFYys1KCsfsHAz939/9Jcf8AQ7bEM2X69OkATJkyhTPOOIOnn35aCV9EBhSJxqkoK6GyvLT/Fn7Zvi38TEi1D38+wYfFNcBFw13J3WPuvgSYCRxtZov7LmNml5rZGjNb09DQkGJ4mdPa2tp9p6vW1lbuu+8+Fi/e522IiHRLJPy+LfxsJ/xUW/ib3X21mU0Ddo505bBEw8PAKcD6Pq9dQ/BBwvLly/Ouy2fHjh2cccYZAESjUc4991xOOeWUHEclIvmsMxqnsjTRwk9O+DEqykowC0or5GvCP8XM/kpQLfNNgpO4gzKzOqArTPbVwAeA76e4/5yZN28e69aty3UYIlJAIrGwS6espHtkDgQjdqrKejpaaquClJxvXToTgK8CXyEYfTMc04CHzOwF4BmCG6DfmeL+RUQKRu8+/N5dOpVJ4+3LSkuoqSzLuxb+t4EF7r7BzGJDLg24+wvAO1Lcn4hIwYpE41SUBn34Hb1a+LHu/vuE8dXluW/hm9lRiWl3r3f3B8Lpr2UiMBGRYpHo0qmuKO2d8MMaO8nGVZf3ugduOo2kS+d5M3vBzL5iZrMyEo2ISJGJxZ1Y3KksK6W6vJT2SO/yyIkx+Anjq8sydqXtSBL+j4CxwJXARjN7yMw+mZGoRESKRKIYWqKF396nhZ+oo5OQF1067v5ld59PcEvDawnKKVyTkahERIpEr4Tft4XfFe+nhZ8HCd/MJpnZp4DvEVxsZQRX2446jY2NnH322SxYsICFCxfy5JNP5jokEclTnbEgwXcn/KQWfkc01l1HJ+GUxVP51HvnZiSWkYzS2U7wAbEXuAG4xd0fy0hUee7zn/88p5xyCrfddhuRSIS2trZchyQieSrRwq8sLWFM2KXj7phZ2MLvnfBPXHAgJy7ITCwjSfi3A7cA97h7Zr5vFICmpiYeeeQRbrzxRgAqKiqoqKjIbVAikreSu3SqKkpxT/Tdl/Z70jaThp3w3f2cTAaSioduvIadb76R1m1OOWgeJ1w4cLnjN954g7q6Oi666CLWrVvHsmXLuOqqqxg7dmxa4xCR4hCJ9e7DB2iPxMKEv28LP5Oyt6ciEY1Gee655/jMZz7D888/z9ixY7nyyitzHZaI5KnuFn5pUsIP+/ETLf1sSeWeth929z9mIpiRGqwlnikzZ85k5syZHHPMMQCcffbZSvgiMqC+wzIhKeH3c6VtJqWyp++mPYoCMnXqVGbNmsWGDRsAePDBBzn88MNzHJWI5Ku+wzKB7qGZHdH4PqN0MimVWjqW9igKzL/+679y3nnnEYlEmDdvHjfccEOuQxKRPNUZ67+FH43Fu6/AzZZUEn7e1ajPtiVLlrBmzZpchyEiBSC5Dz+e1MLve/OTbEi1WqaIiAxDJCmxh/c5ob1LCV9EpOgk9+GXlgQZvz0So7UzqIg5pjJ7aTiVPe1IexQiIkUqeRx+QntXjMa24PrVCeFtDbNhxAnf3U/KRCAiIsUouQ8/uYWfKJA2YUz2rtRXl46ISAYld+kkWvntXTEa2yNAz43Ls0FX2oqIZFByl05FaQkl1reFn+cJ38y+mDR9WPrCyX8bNmxgyZIl3Y9x48bx05/+NNdhiUie6kzq0jGz7hLJiT78bLbwR9SlY2YTgJ8AC8ysA3gBuJigPv6ocNhhh7F27VoAYrEYM2bM4IwzzshtUCKStzrC8gkWjsmsriijvSto4VeWleRvLR13bwQuMrMPAruAI4H/zkBcBeHBBx9k/vz5HHTQQbkORUTyVFskypiKnqReXVFCeyRGLOZZ7c6B1E/adrn7s2a2FdiZzoBGovGPrxPZ2prWbVZMH8uED88f1rKrVq1i5cqVad2/iBSXtkiMMRU9qfaAMRXsbo1QXV7ChOrs3ksj1ZO2p5jZTOBXBF08o04kEuGOO+7gox/9aK5DEZE81h6JddfQAZg2voptje00tnVltf8eUm/hTwC+CnwF+FTaohlpEMNsiWfCPffcw9KlSznwwANzFoOI5L+ghd+T8KdPqOaxV3dRWmLMmjgmq7Gk2sL/NvB7d98AxIZauBjdeuut6s4RkSG1d8W6yyIDzJhQTWskxuY9bVm9yhZST/hfB/42nH4oTbEUjLa2Nu6//37OPPPMXIciInmuvU8Lf9r4agBaI7Gsn7RNNeFHgMTNZE9IUywFY8yYMezevZvx48fnOhQRySMv1DfS3NHVa14wSqen93z6hKru6QPHVZFNqSb8NmC8mZUDs9MYj4hIQWrtjPK/rn6cL6xa22t+35O20ydUd09/cNHUbIUHpJ7wvwm8Dvwc+G36whERKUybdgdDxF/e1tRrfltX7y6duprK7ulsn7RNdZTOP7r7j2H0lVYQEenPxl1Bwp9U03tsfVufFn5JifHPH1rIkTOz3yWcSmmFXwIHhaUV1hEMyxyytIKZzQJuAqYCceAad79qpAGLiOSjTWHCnzi2pwUfizuRaJwx5b1T7SUr5mU1toQRl1Yws3rgEeB/gKMYfmmFKPAld3/OzGqBZ83sfnd/eUQRi4jkoTfChB+Lx7vntUXCu1pVZK9ezmBS6cPfDfwd8Inwef1wVnL3be7+XDjdDLwCzEhh/yIieSfRpdPcEe2e1x4JLlOqLtSE7+5XApcA3wI2Au8d6TbMbA7wDoJvCX1fu9TM1pjZmoaGhpFuOit+8pOfsGjRIhYvXszKlSvp6OjIdUgikmOb97QBvRN+W5jwC7aFb2bfBk4HTgK2uPvPRrh+DfBfwBfcvanv6+5+jbsvd/fldXV1Iw0v47Zs2cLPfvYz1qxZw/r164nFYqxatSrXYYlIjjWFiT55HH6+JfxU7mn7DTM7kKCFfpaZzXf3S4azbjhu/7+A37p7wZZVjkajtLe3U15eTltbG9OnT891SCKSQ12xePetDJuSu3S6gunqivy4m2yqUXwa+Dd3v3e4K1hQ/f864JXEkM79dc8997B9+/Z0bKrb1KlTOfXUUwd8fcaMGVx22WXMnj2b6upqTj75ZE4++eS0xiAihSXRkp84toI9rRE6ozEqy0rzroWf6oVX1wOfMbMfmtmSYa7zHoL6Oyea2drw8aEU958ze/fu5Q9/+AMbN25k69attLa2csstt+Q6LBHJocRonESphEQ/fiLhV2fxrlaDSfnCK4J6OmXAz4AVQ63g7o8BluL++jVYSzxTHnjgAebOnUvi/MKZZ57JE088wfnnn5/1WEQkP7R2Bol96rhKXtkWJPzJNZV0dBX4KJ3Q60AV8Ad3HzLZF5PZs2fz1FNP0dbWhrvz4IMPsnDhwlyHJSI5lGjhTw0rYSZO3BZLl85LwGrgYjN7Jo3x5L1jjjmGs88+m6VLl3LEEUcQj8e59NJLcx2WiORQTwu/d5dOS/hzbGVhn7SdD+wFrgl/jiqXX345l19+ea7DEJE80dPCD8oqJFr4u1sjVJSWUFvgCX+zu682s2nk8CbmIiL5oDXsukmctG1qDz4Adrd0MnFsBcEgxdzTTcxFRPZTW2eQ4GeEte73tEWCn62Rfapn5lKqCX8CPTcx70xbNMPk7tne5Yjke3wikl6JFn5dbSXV5aXsag7S4q7WCJOS6t/n2rATvpkdlfT02wQjdLJ+E/Oqqip2796dt0nV3dm9ezdVVdm9dZmI5E57d1XMMibXVtDQEiT83S2dTB6bPy38kfThP29m64FbgFvd/QEAd/9aRiIbwMyZM6mvrydfC6tB8KE0c+bMXIchIlnSGolRUVpCRVkJk2sq2dWd8CNMLNCE/yPgTOBK4Htm9ihws7tfn5HIBlBeXs7cuXOzuUsRkUG1dUYZUxmMta+rqeTN3W20RaK0d8UKs0vH3b/s7vOB5cC1BFfXXpOpwERECkVrJMbYsEDa5NpKGlo62d0SnLjNp5O2w27hm9kk4AzgbOAEgjIJb2UoLhGRgtEWiXZfTTu5ppK9bRF2NneEzwsw4QPbCb4R7AVuAG4J6+OIiIxqrZ0xxoQXV9XVVuIOr+5oAXrf4zbXRpLwbyc4YXuPu3cNtbCIyGjRFokytiLRhx+06J98YzfQMzY/HwyZ8M1sdjh5Wfhz2gBXjTX2dwcrEZFi19oZY/qEINFPDk/S3vnCNuZOHktdbWG18H8DJAa9D3R9sAM3AjelISYRkYLSGokyNhylc+jUWqrKS+joivPOOQfkOLLehkz47n5CNgIRESlUTe1djKsqB2BcVTkrDqnjvpd3cMTMCbkNrI/8KOEmIlKg3J2mjijjq8u75/3fjyxmbGUZpx0xLYeR7UsJX0RkP7RGYsTi3ivhTxlXxU8+tiR3QQ0g1eJpIiICvN0eDFocV53/7WclfBGR/dAUJvzkFn6+UsIXEdkP3S38KiV8EZGi1tOlo4QvIlLU1KUjIjJKqIUvIjJKNHVEMYPaSo3SEREpak3tXdRWllFSMlDlmfyhhC8ish+a2rsKojsHlPBFRPbL2+1dBXHCFpTwRUT2S6MSvojI6LD97Q6mjqvKdRjDkrWEb2bXm9lOM1ufrX2KiGRSLO7saOpg2gQl/L5uBE7J4v5ERDJqV0sn0bgzdXz+3MZwMFlL+O7+CLAnW/sTEcm0rY3tAEwfrxZ+SszsUjNbY2ZrGhoach2OiMiAtr3dAcA0tfBT4+7XuPtyd19eV1eX63BERAbU3cJXH76ISHHb/nYHVeUlGpYpIlLstjS2M318NWb5X1YBsjss81bgSeAwM6s3s4uztW8RkUzYsL2Zg6fU5DqMYctaeTd3X5mtfYmIZFprZ5SNu1s5fcmMXIcybOrSERFJcse6rVz32MYhl3tlWxPusGj6uCxElR5K+CIiSf7jmbe4fhgJ/6WtTQAsmqGELyJSkBqaO2lo6cTde83f/nYHH7/mSba9HQzFfKH+bSaNrSiYOjqghC8i0ktDcyeRaJym9miv+Xe+sJWn3tjDT+7/KwBPb9rN8jkHFMwIHVDCFxHpFonG2dsW3KN2Z3NHr9faIjEANuxoYUtjO5v3tHPsvElZj3F/KOGLiIR2t3Z2T+9s7uz12lt72gBYv+Vt/rR+OwDHzFXCFxEpSDubkhN+7xb+W7vbKC0xYnHnFw+/zoQx5SyYWpvtEPeLEr6ISKghqVWfnPwhaOF/+MhpjK0oZVdLJysOqSuIG5cnU8IXEQk1tPTfpdPRFWN7Uwfz6mo4Juy3f9+hhVfcUQlfRCSUaOFPH1/VK+HX7w3672dPHMP7F06hqryEFQWY8LNWWkFEJN81NHcyYUw5Mw6opiGpDz9xwnb2pDEsmTmBUxZNZVJNZa7CTJla+CIioYbmTupqKplS27uF/+bunhZ+SYkVZLIHJXwRkW47mzuoq62krraShqSTtm/taWNMRSmTxlbkMLr9p4QvIhJqaOnsTvjNnVHaw4ut3trdxuyJYwrqqtr+KOGLiADuTkNzJ1NqK5lSG3TZJMbiv7UnSPiFTglfRARo6YzS0RWnrraSKWFBtJ3NncTjzlt72jhokhK+iEhRSAzJrEtu4Td1svXtdjqjcWZPGpvL8NJCwzJFREhK+DVVvbp01m0OXj9yxvhchZY2SvgiIvRcZVtXW8kBYyooKzEamjvZsredyrISFk4rnBudDEQJX0SEnto5dbWVlJQYk2sq2fZ2B2/ubuWIGeOpKCv8HvDCfwciImmwq6WT0hJjQnU5AO+cO5EHXtnB+i1NvGP2hNwGlyZK+CIiwN62Lg4YU95dAXPl0bNo7oiCwcqjZ+c4uvRQl46ICNDYFmHCmJ4rad81bxLvXzCF4w6ZzLy6mhxGlj5K+CIiwN62CAeMKe9+bmZcd+E7cxhR+qlLR0QEaGzrYnx1YdfKGYoSvogIQcJPbuEXIyV8ERHCLp0Cr4Y5FCV8ERn1OrpidEbjTCjyFr5O2orIqLe3LQLAAWP6aeHv2Qgv/x52vQrN26B8DNQcCLOOgXnvg9qpw99RPA4Nf4HGt6CjEarGw/hZULcASjOfjpXwRWTU29vaBdDTh9/ZDC//Adb+O7z5eDCvdlrwaN4OGx+BNdeBlcD898PST8BhH+o/absH23juJnj1Pmjfu+8yVeNh3gnwjvOD7ZVkpvMlqwnfzE4BrgJKgWvd/cps7l9EpD+NbRHKiDJ315/h9oeDZN/VBhPnw4n/B476OIyf2bNCPAY71sMrfww+FP7zb2HcDFh2ESy7AGqmQEsDrPv3INHvfg0qx8GC02Due2HyYVA9AdobYc8bsPHPsOGe4JvEhINg+UVw7GehLL3nFMzd07rBAXdkVgr8FTgJqAeeAVa6+8sDrbN8+XJfs2ZNVuITkdHr7he30fK7i1hQ/SwtlTW0zz6W1tlH01wzheauZpojvR9NkSZaulpo62oj7nFi0Q482oHHolTGnWovYUw0xhiPM6a0htrxs6gdP4+ayvHUlI2lpnwsNWVjGVs6hpqyMVSXVhHpaKHjrSdp3/QoFmnhpK+uhRTusGVmz7r78v5ey2YL/2jgNXd/IwxqFXA6MGDCT9Wvv/IjOipi6d6siBS1E9nReiK0AnuAtV3AFgCqgCpqqKMGmDbyTdf3TEZw9tDCHlr6WXAicDrVkZKUkv1QsjlKZwawOel5fTivFzO71MzWmNmahoaGFHcVT3E9EZHcKyktzch2s9nC7+/jap/+JHe/BrgGgi6dVHZ0yQ++nMpqIiJFLZst/HpgVtLzmcDWLO5fRGRUy2bCfwY4xMzmmlkF8HHgjizuX0RkVMtal467R83sH4A/EQzLvN7dX8rW/kVERrusjsN397uBu7O5TxERCaiWjojIKKGELyIySijhi4iMEkr4IiKjRNZq6aTCzBqAN1NcfTKwK43hZEIhxAiKM50KIUZQnOmU7RgPcve6/l7I64S/P8xszUAFhPJFIcQIijOdCiFGUJzplE8xqktHRGSUUMIXERklijnhX5PrAIahEGIExZlOhRAjKM50ypsYi7YPX0REeivmFr6IiCRRwhcRGSWKLuGb2SlmtsHMXjOzr+VBPJvM7EUzW2tma8J5E83sfjN7Nfx5QNLyXw9j32BmH8xgXNeb2U4zW580b8Rxmdmy8P29ZmY/M0vffdkGiPFbZrYlPJ5rzexDOY5xlpk9ZGavmNlLZvb5cH6+HcuB4sy341llZk+b2bowzsvD+XlzPAeJMa+OZb/cvWgeBGWXXwfmARXAOuDwHMe0CZjcZ94PgK+F018Dvh9OHx7GXAnMDd9LaYbiWgEsBdbvT1zA08C7CO5odg9waoZj/BZwWT/L5irGacDScLoW+GsYS74dy4HizLfjaUBNOF0O/A9wbD4dz0FizKtj2d+j2Fr43TdKd/cIkLhRer45HfhNOP0b4CNJ81e5e6e7bwReI3hPaefujxDcqjnluMxsGjDO3Z/04K/3pqR1MhXjQHIV4zZ3fy6cbgZeIbhXc74dy4HiHEiu4nR3T9zduzx8OHl0PAeJcSA5OZb9KbaEP6wbpWeZA/eZ2bNmdmk470B33wbBPyIwJZyf6/hHGteMcLrv/Ez7BzN7IezySXy1z3mMZjYHeAdBiy9vj2WfOCHPjqeZlZrZWmAncL+7593xHCBGyLNj2VexJfxh3Sg9y97j7kuBU4HPmtmKQZbNx/hh4LhyEe8vgfnAEmAb8KNwfk5jNLMa4L+AL7h702CLDhBPruLMu+Pp7jF3X0Jw3+ujzWzxIIvnJM4BYsy7Y9lXsSX8vLtRurtvDX/uBG4n6KLZEX6dI/y5M1w81/GPNK76cLrv/Ixx9x3hP1sc+DU9XV45i9HMygmS6G/d/b/D2Xl3LPuLMx+PZ4K7NwIPA6eQh8ezb4z5fCwTii3h59WN0s1srJnVJqaBk4H1YUwXhItdAPwhnL4D+LiZVZrZXOAQgpM62TKiuMKv1s1mdmw4uuATSetkROKfPnQGwfHMWYzhNq8DXnH3Hye9lFfHcqA48/B41pnZhHC6GvgA8Bfy6HgOFGO+Hct+ZfKMcC4ewIcIRiC8DvxzjmOZR3B2fh3wUiIeYBLwIPBq+HNi0jr/HMa+gQyesQduJfja2UXQ0rg4lbiA5QR/2K8DVxNevZ3BGG8GXgReIPhHmpbjGI8j+Br+ArA2fHwoD4/lQHHm2/E8Eng+jGc98I1U/2cyFecgMebVsezvodIKIiKjRLF16YiIyACU8EVERgklfBGRUUIJX0RklFDCFxEZJZTwZVQwswlm9vdJz6eb2W0Z2tdHzOwbA7zWEv6sM7N7M7F/kYEo4ctoMQHoTvjuvtXdz87Qvr4C/GKwBdy9AdhmZu/JUAwi+1DCl9HiSmB+WKf8h2Y2x8I6+2Z2oZn93sz+aGYbzewfzOyLZva8mT1lZhPD5eab2b1hIbxHzWxB352Y2aFAp7vvCp/PNbMnzewZM/tOn8V/D5yX0XctkkQJX0aLrwGvu/sSd/9yP68vBs4lqH/yXaDN3d8BPElwyTsEN6P+nLsvAy6j/1b8e4Dnkp5fBfzS3d8JbO+z7BrgvSm+H5ERK8t1ACJ54iEP6sQ3m9nbwB/D+S8CR4ZVJt8N/C7ppkSV/WxnGtCQ9Pw9wFnh9M3A95Ne2wlMT0/4IkNTwhcJdCZNx5Oexwn+T0qARg9K4g6mHRjfZ95A9UuqwuVFskJdOjJaNBPc2i8lHtSO32hmH4Wg+qSZHdXPoq8AByc9f5ygaivs219/KD0VFUUyTglfRgV33w08bmbrzeyHKW7mPOBiM0tUP+3v9pmPAO+wnn6fzxPc+OYZ9m35nwDclWIsIiOmapkiaWZmVwF/dPcHhljuEeB0d9+bnchktFMLXyT9vgeMGWwBM6sDfqxkL9mkFr6IyCihFr6IyCihhC8iMkoo4YuIjBJK+CIio4QSvojIKPH/AR0puYBOqeiSAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "swiftdiff['vmag'].sel(id=plidx).plot.line(ax=ax, x=\"time (d)\")\n", + "ax.set_ylabel(\"$|\\mathbf{v}_{swiftest} - \\mathbf{v}_{swifter}|$\")\n", + "ax.set_title(\"Heliocentric velocity differences \\n Planets only\")\n", + "fig.savefig(\"rmvs_swifter_comparison-mars_ejecta-planets-vmag.png\", facecolor='white', transparent=False, dpi=300)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "No handles with labels found to put in legend.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAElCAYAAADgCEWlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqT0lEQVR4nO3de5xdVX338c83MxOSQCCEhJB7uIqAgBgiPF6KKApUi/XSYi1Ui1K8VK1SpLYvRKpV61OtCEqpIqJWH1sR0XIRLwiiIAFDIEBCIEAmCXPJPYQkc2Z+zx97DTlzcmYy53Du+b5fr/OafVl779/ZM3N+Z62199qKCMzMzAaNqXcAZmbWWJwYzMxsCCcGMzMbwonBzMyGcGIwM7MhnBjMzGwIJwYrStKlkr6TpudI2iKprd5xjUTSqyQtrXccsPtYanlOJd0u6T1p+p2Sfpa37hWSHkuxvFnSNEl3SNos6d+qHZs1JieGFiXpSUmvK1j2Lkm/KXVfEfF0ROwTEf2Vi7A0kkLSYSOViYg7I+JFtYppJIWxFP4+6nVOI+K7EfH6vEWXAVekWG4Azgd6gX0j4mO1jM0ahxODtQRJ7fWOoUnNBZYUzD8cZdz56t9B63Bi2INJmiHph5J6JK2Q9KFhys1L39jb87a7UdI6ScslvTevbJukT0h6PDVH3Cdpdlp3pKTb0nZLJf1Z3nbXSrpS0v+m7e6RdGhad0cq9kBq8vhzSadI6pT0cUnPAN8cXJa3z9mSrk/vb62kK4Z5f5dK+h9J/y8d+35Jx+Wtf3FqjtkgaYmkP8lbd6akh9N2qyRdmJY/H4ukbwNzgJ+k+C8q8ZxeKukHkq5Lx1kiaf4Iv9fTJD0qaWN6z8pb93ytUdLjwCF5cX0P+CvgojT/OkljJF2cfp9rUxyTC/4uzpP0NPDLtPyvJT0iab2kWyXNzTt+SLogNV+tT7/z/Pjem7bdnM7rCXnnp+jfqqQFkhZK2iSpS9IXhzs3NkoR4VcLvoAngdcVLHsX8Js0PQa4D7gEGEv2AfEE8Ia0/lLgO2l6HhBAe5r/NfBVYBxwPNADvDat+3vgQeBFZB9IxwEHAHsDK4F3A+3ACWRNFken7a4F1gEL0vrvAt/Piz2Aw/LmTwFywOeBvYDxaVlnWt8GPAB8KR17HPDKYc7VpUAf8DagA7gQWJGmO4DlwCfSeToV2Ay8KG27BnhVmt4fOCEvvs7hfh8lntNLgW3Amel9fRa4e5j3MgXYlPde/i6dp/cU/g0ME9e1wKfz5j8C3A3MSuf5P4DvFbyH69I5Hg+8OZ2vF6ff4z8Bvy34Pf4UmESWLHuA09O6twOrgBPJ/nYOI6vB7O5v9XfAOWl6H+Ckev//Nfur7gH4VaVfbPYPvwXYkPfays7E8HLg6YJt/gH4Zpq+lCKJAZgN9AMT87b7LHBtml4KnFUknj8H7ixY9h/AJ9P0tcDX89adCTyaN18sMewAxhUsG0wMJ6cPnfZRnKtLyfugTR9Ea4BXpdczwJi89d8DLk3TTwN/Q9YmT7FY8n4fRRPDKM7ppcDP89YdBTw3zHs5t+C9COik/MTwCClBpfnpZEm0Pe89HJK3/mbgvIJzuRWYm/d7fGXe+h8AF6fpW4EPF3lPu/tbvQP4FDCl3v93rfJyU1Jre3NETBp8Ae/PWzcXmJGaRzZI2kD2rXjabvY5A1gXEZvzlj0FzEzTs4HHi2w3F3h5wfHeCRyUV+aZvOmtZN/+RtITEduGWTcbeCoicrvZx6CVgxMRMUD2YTojvVamZYPy3+9byZLYU5J+LenkUR4v3+7OKex6bsapeJv+jIL3EvnzZZgL/Cjvd/YIWRLL/ztZWVD+y3nl15Elp5Hey+DveaS/nZH+Vs8DjgAelXSvpDeW/C5tCHcW7blWAisi4vASt1sNTJY0Me+DbA5ZE8Dgfg8FHipyvF9HxGnlBlzESB2kK4E5ktpHmRxmD05IGkPWdLJ6cJ2kMXnJYQ6wDCAi7gXOktQBfJDsG/Dz+xplrLs7p6VYU/BeNEw8o7US+OuIuKtwhaR5aTIKyn8mIr5b5rEOHWb5sH+rEfEY8I70e3sL8D+SDoiIZ8uIwXDn857s98Cm1Hk7Xlmn8TGSThxpo4hYCfwW+KykcZKOJfvGNvhB8HXgnyUdrsyxkg4ga1c+QtI5kjrS60RJLx5lvF1kbculvL81wOck7Z1ifcUI5V8m6S3pW/hHgO1kbev3AM+Sdch2SDoFeBPwfUljld0XsF9E9JG17Q93+emw8Y/inJbif4Gj897LhxhaKyvVVcBnBjuQJU2VdNZuyv+DpKNT+f0kvX2Ux/o6cKGkl6W/ncPScUf8W5X0l5KmpsS9Ie2rbpdWtwInhj1UZNfPv4mso3MFWUfw14H9RrH5O8jal1cDPyLrJ7gtrfsi2bfmn5F9UH4DGJ++Cb8eODtt9ww7O45H41LgW6kp4c92Vzjv/R1G1g/QSdbPMZwfp/XrgXOAt0REX0TsAP4EOIPsHH0VODciHk3bnQM8KWkTcAHwl8Ps/7PAP6X4LyyyfqRzOmoR0UvWifs5YC1wOLDLt/0SfBm4EfiZpM1kyfLlIxz/R2S/1++nc/IQ2bkbTez/DXwG+C+yDv4bgMmj+Fs9HVgiaUuK9+wRmhhtFJQ6b8z2WJIuJevYHu5D3WyP4hqDmZkN4cRgZmZDuCnJzMyGcI3BzMyGcGIwqyIVDHM9QrnnhzlvBMrGrvp0veOw+nBisIahnc8oGHyFpGfz5l9Vxj53GX68YP0pkgbS/jcrG9zv3WXGP2RgPCg6zLVZw/Odz9YwIuJp8obBkBTAcRGxvMqHXh0Rs9JdwmeR3Tl7T0Q8PNodDDM8hVlTco3BmoKkvST9X0lPKxta+SpJ49O6KZJ+mm4eWyfpTmXDRe8y3PVIx4jMDWQ3uR0l6Y8l/UHZcM4r0/0Og/EUG3J6cHjwDel4J6vg4UiSjtbOoce7JH1imPd7kqTfpvf0QLrjenDduyQ9kWo4KyS9c4Rz9u+SVqfXv0vaK60bHLb8Y5K6Ja0ZrqYk6SFJb8qb75DUK+n4kc6nNS8nBmsWnycbKO14sruZZ5INwwzwMbI7m6eSDaz2CbLP+XPI7np+U2RPKPvXkQ6Qksmfkg0J/SDZUBjnpvk/Bt4n6c0Fm/0R2RDTbwBenZZNSsf7XcH+JwI/B24hG+zuMOAXReKYSTa0xaeByWTDgP8wDUexN3A5cEZETAT+D7BomLf0j8BJZOfsOLIhzf8pb/1BZHcPzyQbguNKSfsX2c91DL2j+0xgTUQMd1xrci2RGCRdk771FA7cVu7+5kj6mbIHhjysnYOFWR2kJp73An8XEYOjkP4L2fAakA0DPZ1saOe+yB6rWcp12DOUjdjZC3ySbGz/pRFxe0Q8GBEDEbGYbLjtPyrY9tKIeDYinhvFcd4IPBMR/xYR2yJic0TcU6TcXwI3RcRN6di3AQvJPpABBoBjJI2PiDURsaTIPiAbvfayiOiOiB6yoanPyVvfl9b3RcRNZMO0F3s06neAMyXtm+bPAb49ivdrTaolEgPZGPKnV3B/1wFfiIgXk33L6q7gvq10U4EJwH3aOezyLWk5wBfIHg7zs9TEcnGJ+1+dhiafHBHHR8T3ASS9XNKvlD01bCPZWEhTCrYtZUjr4YaVLjQXeLuGDjP9SmB6GjH0z1Msa5Q98e7IYfYzg2z47kFPpWWD1haMPFt0qPOIWE023tJbJU0iG/uonAH+rEm0RGKIiDvIxn1/nqRDJd2i7NGSd47wzzOEpKPIHu5yW9r3lojYWvmorQS9wHNkT3sbfL7EfhGxD0D65v2xiDiEbLC1j0p6bdr2hdzB+V9kA8jNjoj9yEYOVUGZGGa6mOGGlS5W7tv5z9KIiL0j4nMAEXFrGr58OvAo8J/D7Gc1WZIZNIedQ4mX6ltkNZm3A7+LiHKGBLcm0RKJYRhXA38bES8ja6P96ii3O4Ks8/D61PH4BUltVYvSdisNp/yfwJckHQhZO7ykN6TpNyobolnsHPp6cNjlUofrzjeR7AE62yQtAP5iN+V7yJp5hjveT4GDJH0kdQxPlFRspNLvAG+S9AZlQ0yPS53FsyRNk/Qnqa9hO1nzz3BDTH+PbETXqZKmkPXJlHuvxA1kj2P9MFmN2lpYSyYGSfuQdcr9t6RFZI+QnJ7WvSVdZVH4ujVt3k72OMcLyZ49ewjZ4xCtvj5O1lx0t7LhnH/Ozvbww9P8FrLn/341Im5P63Y33PVI3g9cpmy46UvIhhMfVqpZfga4Kx3vpIL1m4HTyGo1zwCPAa8psp+VZJfNfoIs2awke5b2mPT6GNk3/3VkfR7vL9xH8mmyvonFZJ3p96dlJUt9KD8EDgauL2cf1jxaZqyk1EH804g4JnWSLY2I6WXs5yTgcxFxSpo/h+zh4h+oZLxmzUbSJcARHp689bVkjSEiNgErlJ4cpcxxo9z8XmB/SYMdm6cCo77RyawVSZpMdknr1fWOxaqvJRKDpO+RNSG8KN20cx7ZpXrnSXoAWEJWNd+t9LSoC4FfSHqQrLNxuM49s5Yn6b1kzVk3pws9rMW1TFOSmZlVRkvUGMzMrHKafuCvKVOmxLx58+odhplZU7nvvvt6I2JqsXVNnxjmzZvHwoUL6x2GmVlTkfTUcOvclGRmZkM4MZiZ2RBODGZmNoQTg5mZDeHEYGZmQzgxmJnZEE4MZmY2hBODmVkT6rnySrbcdVdV9u3EYGbWZGJggN4rv8rWKt3c68RgZtZk+jduhIEB2vffvyr7d2IwM2sy/es3ANDmxGBmZgD9G9YD0Lb/5Krs34nBzKzJ9K9bB0Db/pOqsn8nBjOzJpNbn9UY3MdgZmaA+xjMzKxA/7p1aPx4xowfX5X9OzGYmTWZ/vXrq9a/ADVMDJLGSfq9pAckLZH0qSJlTpG0UdKi9LqkVvGZmTWL3Ib1tFfpiiSo7aM9twOnRsQWSR3AbyTdHBF3F5S7MyLeWMO4zMyaSv+69VXrX4Aa1hgisyXNdqRX1Or4ZmatImtKaoHEACCpTdIioBu4LSLuKVLs5NTcdLOko4fZz/mSFkpa2NPTU82QzcwaTm7tWtonV68pqaaJISL6I+J4YBawQNIxBUXuB+ZGxHHAV4AbhtnP1RExPyLmT506tZohm5k1lIGtW4nnnqNtygFVO0ZdrkqKiA3A7cDpBcs3DTY3RcRNQIekKTUP0MysQeXWrgWgfXILJAZJUyVNStPjgdcBjxaUOUiS0vSCFN/aWsVoZtbocr29ALRXscZQy6uSpgPfktRG9oH/g4j4qaQLACLiKuBtwPsk5YDngLMjwh3UZmbJ8+MkHVC9xpSaJYaIWAy8tMjyq/KmrwCuqFVMZmbNJtebmpJarY/BzMzKk1ubmpJa5aokMzN7Yfp71zJm333R2LFVO4YTg5lZE8mtW0f7AdVrRgInBjOzptLf2+vEYGZmO+XWrqVtSnVv73JiMDNrItUeDgOcGMzMmkbs2MHApk1VHQ4DnBjMzJpGLt3c1l7Fm9vAicHMrGnU4uY2cGIwM2sa/euyxNDmPgYzM4P8GoObkszMjLzhMHwfg5mZQTYchiZMYMyECVU9jhODmVmTyK1bV/V7GMCJwcysaeR6e6rejARODGZmTSPX1U37tGlVP44Tg5lZk8h1dTkxmJlZpn/Lsww8+ywd0w6s+rFqlhgkjZP0e0kPSFoi6VNFykjS5ZKWS1os6YRaxWdm1shy3d0ANakx1OyZz8B24NSI2CKpA/iNpJsj4u68MmcAh6fXy4GvpZ9mZnu0XHcXAO1TW6jGEJktabYjvaKg2FnAdans3cAkSdNrFaOZWaPKdaXE0EpNSQCS2iQtArqB2yLinoIiM4GVefOdaVnhfs6XtFDSwp6enqrFa2bWKPq6sqakjlbrfI6I/og4HpgFLJB0TEERFdusyH6ujoj5ETF/6tSpVYjUzKyx5Lq6GDNxYtXveoY6XZUUERuA24HTC1Z1ArPz5mcBq2sTlZlZ48p1d9ekGQlqe1XSVEmT0vR44HXAowXFbgTOTVcnnQRsjIg1tYrRzKxR9XV30XFg9ZuRoLZXJU0HviWpjSwh/SAifirpAoCIuAq4CTgTWA5sBd5dw/jMzBpWrqubvU4+tCbHqlliiIjFwEuLLL8qbzqAD9QqJjOzZhD9/eR6emg/sMWakszMrDy5tWuhv7/1+hjMzKw8uRpeqgpODGZmDe/5u55r1PnsxGBm1uB2jpPkpiQzMwP6urqgra0mD+kBJwYzs4bXt3o1HdOmoba2mhzPicHMrMH1rVpNx8xdho2rGicGM7MG17dqlRODmZllYscOcl1dTgxmZpbpe+YZiKBj1qyaHdOJwcysgfV1dgLQMXNGzY7pxGBm1sB2rFoFwFg3JZmZGUDf0yuho4P2Gg2HAU4MZmYNbfuKJxg7Zw5qr91TEpwYzMzqbOu999Jz+eVse/jhXdbtWPEkYw+eV9N4avmgHjMzKxADA6z66MfI9fSw9d6FzP32dTvX5XLsePppJp56ak1j2m1ikDRnlPvaEBGbXmA8ZmZ7lG1LHibX00PHrFlsXbgweyDP1KlAuiKpr4+xBx9c05hGU2P4FhCARigTwLXAdcMVkDQ7rT8IGACujogvF5Q5BfgxsCItuj4iLhtFjGZmTWnLr34FY8Yw47P/wlPnnMumn/2Mye98JwDbn8g+ChuuKSkiXlO4TNJBEfFMicfKAR+LiPslTQTuk3RbRBQ2qt0ZEW8scd9mZk3p2XvuYdxLjmHCiSfSMWcOz/72dzsTw9JHAdjrsMNqGlO5nc/nlrpBRKyJiPvT9GbgEaB2F+aamTWYiGD70qWMP/poACacOJ+tCxcSAwMAPPfAYsYecghtEyfWNK5yE8NZkj4o6UXlbCxpHvBS4J4iq0+W9ICkmyUdPcz250taKGlhT09POSGYmdVd36rVDGzZwl5HZB+ley9YwMDGjWxftoyI4LnFixl/7LE1j6vcxPAWYDnwp5K+XsqGkvYBfgh8pEhn9f3A3Ig4DvgKcEOxfUTE1RExPyLmT02dNGZmzWb7sqUA7PWiIwCYcOKJADx712/pW7WK/nXrGH9c7RNDWZerRkQXcEt6jZqkDrKk8N2IuL7IfjflTd8k6auSpkREbzlxmpk1su1Ls8Qw7ogsMXTMmMG4o45i0003ofHjABh/wgk1j6usGoOkKyVdm6ZfP8ptBHwDeCQivjhMmYNSOSQtSPGtLSdGM7NGt/2xx+iYNYsxe+/9/LJ93/Qmti1ZQtdl/8y4l7yEvVLSqKVym5J2AE+k6dHeefEK4BzgVEmL0utMSRdIuiCVeRvwkKQHgMuBsyMiyozRzKyhbX/yyV3uUdjvjX9M2377ATD53HNI35Vrqtw7n7cC+6WmoVHdABcRv2HkeyGIiCuAK8qMycysaUQEfU8+xYQTXjZkefvUqRz269vZ/thjjDvmmLrEVm5iWAc8B1wJ3FW5cMzM9gz9vb0MbN3K2Llzd1k3Ztw4xr/kJXWIKh2/lMKSJkn6JvDWtOg6YH7FozIza3E7nnoKgLHz5tU3kCJKqjFExAZJnwPmAb3AscAuVxeZmdnIdiaGXWsM9VZOU9J5wIqIuBW4r8LxmJntEXY8+RR0dNAxfXq9Q9lFOYlhPXBBuuv5AWBRRPyhsmGZmbW2HU8+ydhZs2r6AJ7RKjmiiPispF8Ay4DjgVcDTgxmZiXY8dRTRTueG0HJiUHSZUAbsIistnB7hWMyM2tpMTDAjqefZu+TT653KEWVU2O4RNI0skHw3irp0Ih4b+VDMzNrTbnubmLbtobseIby72P4G+A/IqKksZLMzCx1PEPrNCUl1wDvk7Q32YB4iyoXkplZa2vkexig/LGSPkSWVNrJxjQyM7NR6lv5NOrooP2gg+odSlHlJobHgXHAjyPi1RWMx8ys5e3oXEXHzJloTLkfwdVVblRLgF8C50m6t4LxmJm1vL6VK+mYNaveYQyr3D6GI4Ae4GqyG97MzGyU+jo7GXds/QbJ251yawxHkt3UdiFwfuXCMTNrbf2bN9O/cSNjG7jGUG5imAR8HLgI2FaxaMzMWlxfZycAHTMbNzGU25R0GXBkRCyVNFDJgMzMWtmOwcQwu3ETw6hqDJLaJK2R9B6AiOiMiJ+n6YurGaCZWSvpW5klhqZvSoqIfuAh4NByDyRptqRfSXpE0hJJHy5SRpIul7Rc0mJJJ5R7PDOzRtTX2cmYiROff65zIyqlKWkCcJGk04DVaVlExFmj3D4HfCwi7pc0EbhP0m0R8XBemTOAw9Pr5cDX0k8zs5awY1VnQzcjQWmJYXAYwBPSCyBGu3FErAHWpOnNkh4BZgL5ieEs4LqICODu9CjR6WlbM7Om17eyk70OLbvxpSZKSQwHV+qgkuaRjc56T8GqmcDKvPnOtGxIYpB0Puky2Tlz5lQqLDOzqoqBAfpWrWKfU06pdygjGnViiIinKnFASfsAPwQ+EhGbClcXO3SRWK4mu7mO+fPnj7rWYmZWT7meXmL7djpmzax3KCOq6UAdkjrIksJ3I+L6IkU6gdl587PY2Z9hZtbU+lY1/hVJUMPEIEnAN4BHIuKLwxS7ETg3XZ10ErDR/Qtm1ir6Vq0CoGNmY9cYynm055si4idlHOsVwDnAg5IWpWWfAOYARMRVwE3AmcByYCvw7jKOY2bWkPpWZ99zO6ZPr3MkIyvnzufPACUnhoj4DcX7EPLLBPCBMmIyM2t4fWtW07b//oyZMKHeoYyonKakET/czcysuL7Vqxu+tgDlJQZfBWRmVobcmjW0z2jNxGBmZiWKCPpWraZj+ox6h7JbTgxmZjUwsHkzA1u30jGjNRNDV8WjMDNrcX2rs1uyWrKPISJOq0YgZmat7PlLVd3HYGZmkF2qCi1aYzAzs9L1rV6Nxo6l7YAD6h3KbpWVGCR9NG/6RZULx8ysNeXWrKF9+kFoTON/Hy/pzmdJk4AvAUdK2gYsBs7DQ1eYmY2ob/WaprhUFUqsMUTEhoh4N/BpsmcpvAooNkqqmZnl6Vuzpin6F6D8PoY/Irts9STAVymZmY0gduwg193dFPcwQPmJYRLwceAiYFvFojEza0F93d0QQcf0g+odyqiUM7oqwGXAkRGxVNJAJQMyM2s1uWeeAaB9WgsnhojoJHvaGhFxcUUjMjNrMX1d2YAR7dMOrHMko1Pu5apXSro2Tb++ohGZmbWYXFc3AB3TptU5ktEpt49hB/BEmj61QrGYmbWkXFcXGjeOMfvuW+9QRqXcxLAV2E9SB+nRnLsj6RpJ3ZIeGmb9KZI2SlqUXpeUGZuZWUPp6+6iY9o0pOZ4zlm5nc/rgOeAK4G7RrnNtcAVwHUjlLkzIt5YZkxmZg0p19VNe5M0I0GJNQZJkyR9E3hrWnQdMH8020bEHWQJxcxsj5Lr6mqqxFBSjSEiNkj6HDAP6AWOpbJ3Pp8s6QFgNXBhRCwpVkjS+cD5AHPmjKoly8ysLiIiu7mtSa5IgvKaks4DVkTErcB9FYzlfmBuRGyRdCZwA3B4sYIRcTVwNcD8+fP9DGoza1j969cTfX20H9g8NYZyOp/XAxdI+ndJ75b00koEEhGbImJLmr4J6JA0pRL7NjOrl9zz9zA0T2IoucYQEZ+V9AtgGXA88GrgDy80EEkHAV0REZIWkCWttS90v2Zm9TR4c1tLNyVJugxoAxYBiyLi9lFu9z3gFGCKpE7gk0AHQERcBbwNeJ+kHNkVT2dHhJuJzKypDd7c1uo1hkvSPQZjgLdKOjQi3juK7d6xm/VXkF3OambWMnJdXSDRPqV5WsbLvcHtGuDFwAHAVysXjplZa+nr7qJtygGoo6PeoYxauYnhQ2S1jXbgy5ULx8ysteS6uulooiuSoPzE8DgwDvhxRLy6gvGYmbWUZru5DcpPDEuAXwLnSbq3gvGYmbWULDE0zxVJUP5YSYeS3c9wdfppZmYFBnbsoH/jRjoO3DMSw8qI+KWk6UB3JQMyM2sV/b29ALRPnVrnSEpTblPS6ZJmAVcBX6pgPGZmLSPX0wNAWxNdqgrlJ4ZJwMeBi4DtFYvGzKyF5AZrDFOaq8ZQblPSZcCREbFUUn8lAzIzaxW5nsGmpBasMUhqk7RG0nsAIqIzIn6epi+uZoBmZs3q+RrD5Ml1jqQ0o0oMEdEPPER2NZKZmY1CrreHtv33b6q7nqG0pqQJwEWSTiN7kA5ARMRZlQ/LzKz55Xp7m2qMpEGlJIaT088T0gvAo5+amQ2jv6e36foXoLTEcHDVojAza0G53l7Gzz1h9wUbzG4Tg6TBhyoXrR3krd8QEZsqFZiZWTOLiNSU1FyXqsLoagzfIksKGqFMANcC11UgJjOzpjeweTOxfXvT3fUMo0gMEfGaWgRiZtZKdt7c1nx9DOXe+WxmZiNo1pvboIaJQdI1krolPTTMekm6XNJySYslNV+PjZlZkuvNxklyjWFk1wKnj7D+DODw9Dof+FoNYjIzq4p+NyXtXkTcAawbochZwHWRuRuYlIb1NjNrOrneXujoYMx++9U7lJI1Uh/DTGBl3nxnWrYLSedLWihpYU8a1tbMrJHkerK7nqWRLuhsTI2UGIqdvaL3TkTE1RExPyLmT23CS8HMrPU163AY0FiJoROYnTc/i51jMpmZNRUnhsq4ETg3XZ10ErAxItbUOygzs3LkenqaNjGU+6Cekkn6HnAKMEVSJ/BJoAMgIq4CbgLOBJYDW4F31yo2M7NKilyO/nXrmvIeBqhhYoiId+xmfQAfqFE4ZmZVk1u3DiKa7lnPgxqpKcnMrCU8fw9Dk14c48RgZlZhzTxOEjgxmJlV3M5xklxjMDMz8moMBxxQ50jK48RgZlZhud5exuyzD2PGj693KGVxYjAzq7Bcb/PewwBODGZmFdfMN7eBE4OZWcXlenpoP7A5O57BicHMrOJyPb1Ne0USODGYmVVU/5Znia1bnRjMzCyT6+kGaNrhMMCJwcysopp9OAxwYjAzq6hceqqkE4OZmQFODGZmViDX04M6OmibNKneoZTNicHMrIJyPT20TZ2CVOwx9s3BicHMrIJyPT1N3YwENU4Mkk6XtFTSckkXF1l/iqSNkhal1yW1jM/M7IVqhcRQy2c+twFXAqcBncC9km6MiIcLit4ZEW+sVVxmZpWU6+ll/Pz59Q7jBalljWEBsDwinoiIHcD3gbNqeHwzs6qKHTvo37ChqQfQg9omhpnAyrz5zrSs0MmSHpB0s6Sji+1I0vmSFkpa2JMuDTMzq7dcC9zcBrVNDMW66KNg/n5gbkQcB3wFuKHYjiLi6oiYHxHzpzb5L8DMWkcr3MMAtU0MncDsvPlZwOr8AhGxKSK2pOmbgA5JzV0nM7M9xs7EcGCdI3lhapkY7gUOl3SwpLHA2cCN+QUkHaR08a+kBSm+tTWM0cysbK1SY6jZVUkRkZP0QeBWoA24JiKWSLogrb8KeBvwPkk54Dng7IgobG4yM2tIuZ5ekGg/YHK9Q3lBapYY4PnmoZsKll2VN30FcEUtYzIzq5QdnStpnzYNtdf0o7XifOezmVmZIpdj0y23MPDccwDsePwJ9jrkkDpH9cI5MZiZlWnDD69n1Uf+jtUXfZwYGGD7ihWMPfTQeof1gjkxmJmVIQYGWHfttWjCBDbfdhsbf/QjYutW9jrUNQYzsz3StkceYceKFUy76CLG7LMPPV/JukfHHuzEYGa2R9r24IMA7P3KVzLxta8l98wzAC1RY2jurnMzszp5bvGDtO2/Px0zZ7D/OefQt2oVe7/ylU0/ThI4MZiZlWXbg4sZd+xLkMT4Y45m7ne+Xe+QKsZNSWZmJRp49lm2L3+c8ce8pN6hVIUTg5lZibY/9hhEMO7FR9Y7lKpwYjAzK9G2ZcsA2OuII+ocSXU4MZiZlWj7ssfQhAl0zJpV71CqwonBzKxE25ctY6/DDkNjWvMjtDXflZlZlURElhiOOLzeoVSNE4OZWQlyPT30b9jAuBbtXwAnBjOzkmxf9hjQuh3P4MRgZlaS7S1+RRI4MZiZlWT7smW0TZ1C++TmfkrbSJwYzMxKsH3ZMsYd3rq1BahxYpB0uqSlkpZLurjIekm6PK1fLOmEWsZnZjaS3Pr1bFu6lHHHHVvvUKqqZolBUhtwJXAGcBTwDklHFRQ7Azg8vc4Hvlar+MzMdmfLL34B/f3se9pp9Q6lqmo5uuoCYHlEPAEg6fvAWcDDeWXOAq6LiADuljRJ0vSIWFPpYL7yF39N/0B/pXdrZq3u+AVw2efrHQUAbWPa+Nv/uqbi+61lU9JMYGXefGdaVmoZJJ0vaaGkhT09PRUP1MxsT1bLGoOKLIsyyhARVwNXA8yfP3+X9aNRjSxrZtYKallj6ARm583PAlaXUcbMzKqolonhXuBwSQdLGgucDdxYUOZG4Nx0ddJJwMZq9C+YmdnwataUFBE5SR8EbgXagGsiYomkC9L6q4CbgDOB5cBW4N21is/MzDI1feZzRNxE9uGfv+yqvOkAPlDLmMzMbCjf+WxmZkM4MZiZ2RBODGZmNoQTg5mZDaGsv7d5SeoBnipz8ylAbwXDqRbHWTnNECM4zkpqhhih9nHOjYipxVY0fWJ4ISQtjIj59Y5jdxxn5TRDjOA4K6kZYoTGitNNSWZmNoQTg5mZDbGnJ4ar6x3AKDnOymmGGMFxVlIzxAgNFOce3cdgZma72tNrDGZmVsCJwczMhthjE4Ok0yUtlbRc0sV1juVJSQ9KWiRpYVo2WdJtkh5LP/fPK/8PKe6lkt5QxbiukdQt6aG8ZSXHJell6f0tl3S5pGIPZKp0nJdKWpXO6SJJZ9YzTkmzJf1K0iOSlkj6cFreUOdzhDgb5nxKGifp95IeSDF+Ki1vtHM5XJwNcy6HFRF73Its2O/HgUOAscADwFF1jOdJYErBsn8FLk7TFwOfT9NHpXj3Ag5O76OtSnG9GjgBeOiFxAX8HjiZ7Al9NwNn1CDOS4ELi5StS5zAdOCEND0RWJZiaajzOUKcDXM+0/72SdMdwD3ASQ14LoeLs2HO5XCvPbXGsABYHhFPRMQO4PvAWXWOqdBZwLfS9LeAN+ct/35EbI+IFWTPrlhQjQAi4g5g3QuJS9J0YN+I+F1kf+HX5W1TzTiHU5c4I2JNRNyfpjcDj5A9z7yhzucIcQ6n5nFGZkua7UivoPHO5XBxDqdu/0OF9tTEMBNYmTffych//NUWwM8k3Sfp/LRsWqSn16WfB6bl9Y691LhmpunC5bXwQUmLU1PTYLNC3eOUNA94Kdk3yIY9nwVxQgOdT0ltkhYB3cBtEdGQ53KYOKGBzmUxe2piKNY+V8/rdl8REScAZwAfkPTqEco2WuyDhourXvF+DTgUOB5YA/xbWl7XOCXtA/wQ+EhEbBqp6DDx1CvOhjqfEdEfEceTPRd+gaRjRihet3M5TJwNdS6L2VMTQycwO29+FrC6TrEQEavTz27gR2RNQ12pCkn62Z2K1zv2UuPqTNOFy6sqIrrSP+UA8J/sbG6rW5ySOsg+bL8bEdenxQ13PovF2YjnM8W1AbgdOJ0GPJfF4mzUc5lvT00M9wKHSzpY0ljgbODGegQiaW9JEwengdcDD6V4/ioV+yvgx2n6RuBsSXtJOhg4nKxjqlZKiitV6TdLOildSXFu3jZVM/gBkfwp2TmtW5xpn98AHomIL+ataqjzOVycjXQ+JU2VNClNjwdeBzxK453LonE20rkcVjV7thv5BZxJdsXF48A/1jGOQ8iuRHgAWDIYC3AA8AvgsfRzct42/5jiXkoVr04AvkdW1e0j+9ZyXjlxAfPJ/vgfB64g3XFf5Ti/DTwILCb7h5tezziBV5JV/xcDi9LrzEY7nyPE2TDnEzgW+EOK5SHgknL/Z6p8LoeLs2HO5XAvD4lhZmZD7KlNSWZmNgwnBjMzG8KJwczMhnBiMDOzIZwYzMxsCCcGszySJkl6f978DEn/U6VjvVnSJcOs25J+TpV0SzWObzYcJwazoSYBzyeGiFgdEW+r0rEuAr46UoGI6AHWSHpFlWIw24UTg9lQnwMOTePkf0HSPKXnPEh6l6QbJP1E0gpJH5T0UUl/kHS3pMmp3KGSbkmDIt4p6cjCg0g6AtgeEb1p/mBJv5N0r6R/Lih+A/DOqr5rszxODGZDXQw8HhHHR8TfF1l/DPAXZOPbfAbYGhEvBX5HNlQBZA91/9uIeBlwIcVrBa8A7s+b/zLwtYg4EXimoOxC4FVlvh+zkrXXOwCzJvOryJ5TsFnSRuAnafmDwLFpVNL/A/x33kO29iqyn+lAT978K4C3pulvA5/PW9cNzKhM+Ga758RgVprtedMDefMDZP9PY4ANkQ21PJLngP0Klg03Ps24VN6sJtyUZDbUZrJHWpYlsmcXrJD0dshGK5V0XJGijwCH5c3fRTbKL+zan3AEO0fgNKs6JwazPBGxFrhL0kOSvlDmbt4JnCdpcMTcYo+NvQN4qXa2N32Y7CFN97JrTeI1wP+WGYtZyTy6qlmdSPoy8JOI+Pluyt0BnBUR62sTme3pXGMwq59/ASaMVEDSVOCLTgpWS64xmJnZEK4xmJnZEE4MZmY2hBODmZkN4cRgZmZDODGYmdkQ/x+1M5mJ0FlzqAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "swiftdiff['rmag'].sel(id=tpidx).plot.line(ax=ax, x=\"time (d)\")\n", + "ax.set_ylabel(\"$|\\mathbf{r}_{swiftest} - \\mathbf{r}_{swifter}|$\")\n", + "ax.set_title(\"Heliocentric position differences \\n Test Particles only\")\n", + "legend = ax.legend()\n", + "legend.remove()\n", + "fig.savefig(\"rmvs_swifter_comparison-mars_ejecta-testparticles-rmag.png\", facecolor='white', transparent=False, dpi=300)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "No handles with labels found to put in legend.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAElCAYAAADgCEWlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAriklEQVR4nO3deZxcdZ3v/9e7qzvpTrqzhy0JJGwicoGBgHhxAUdGYMYfOi6DOoqKIuPoeH9uMDP+HJfxqtfRmfHnwmW4iKgjDx03HBFcccGNBMISMpEshDQB0mRfupPuqs/945xOqitVne6iqk519fv5eNSjz/I9pz59urs+/f1+z/l+FRGYmZkNa8s6ADMzay5ODGZmNoITg5mZjeDEYGZmIzgxmJnZCE4MZmY2ghODlSXpg5K+ki4fK2m3pFzWcY1G0vMkrW7we4akE5/mOVZKuqA2ER1y7oo/R0lHSvqFpF2SPqXEFyVtk/T7esRjE4MTQ4uS9IikF5Vse4OkX433XBHxaER0R0S+dhGOz1g+gCPilxHxjEbFVCsR8ayIuBNGfpDX4X1Kf45XAU8BMyLi3cBzgYuAhRFxbj1isInBicFagqT2rGOYgI4DHoqDT7keBzwSEXvGeyJf/9bixDCJSTpG0jcl9UlaL+lvKpRbnP7H3l503K2StkpaI+ktRWVzkv5O0tq0iWK5pEXpvlMk/Sg9brWkVxUdd5Okz0n6fnrc7ySdkO77RVrsvrQp5C8kXSCpV9I1kp4Avji8reiciyR9K/3+tkj6bIVr0C9pTtG2P5L0lKSOdP1NklalTSx3SDquwnWaKenm9P02SHq/pLai/W9Jz7NL0kOSzkq3PyLpRZIuBv4O+Iv0+7xP0islLS95n3dL+k6FGJZI+nn6Hj8C5pX7OUq6CbgCeF/6Xm8FbgCek65/KD3mzyStkLRd0q8lnV50vkfS638/sCc973lpue1p/BcUlb9T0kck3ZXG90NJxfE9t+jYjZLekG6fKumfJD0q6UlJ10nqSvfNk/Sf6TFbJf2y+JpblSLCrxZ8AY8ALyrZ9gbgV+lyG7Ac+AAwBTgeWAe8ON3/QeAr6fJiIID2dP3nwOeBTuBMoA/443Tfe4EHgGcAAs4A5gLTgY3AG4F24CySZoxnpcfdBGwFzk33fxW4pSj2AE4sWr8AGAI+AUwFutJtven+HHAf8M/pe3cCz61wrX4KvKVo/ZPAdenyS4E1wDPTuN4P/LpcXMDNwHeBnvSa/QG4Mt33SuAx4Jz0upwIHFf6syq+7un61PS6PLNo273Ayyt8L78BPp0e93xg1yg/x5uAfyz3+5GunwVsBp6dXs8r0linFsW9AliUXv8FwBbgUpLfr4vS9flp+TuBtcDJafk7gY+n+45NY3010EHyO3Nmuu9fgFuBOem1/R7wsXTfx4Dr0mM6gOcByvrvb6K/Mg/Arzr9YJM/2t3A9qLXXg4mhmcDj5Yc87fAF9PlAx9QxR8o6YdAHugpOu5jwE3p8mrgsjLx/AXwy5Jt/xv4h3T5JuCGon2XAv9VtF4uMewHOku2DSeG55AkrPYxXKs3Az9Nl0WSwJ6frv+A9MM9XW9Lr+NxxXGRfHDuA04tKvtW4M50+Q7gnaP8rMomhnTbF4CPpsvPAraRfjiXlDuWJFlOL9r27+V+jkXXfLTE8AXgIyXvsRp4QVHcbyradw3w5ZLydwBXpMt3Au8v2vc24Pai371vl/meBOwBTija9hxgfbr8YZJkfGLpsX5V/3KVq7W9NCJmDb9I/hCHHQcck1bBt0vaTtKMceRhznkMsDUidhVt20Dy3yIkiWNtmeOOA55d8n6vBY4qKvNE0fJeoPswsfRFxECFfYuADRExdJhzAPwHSRPKMST/ZQfwy6K4/7Uo5q0kH1YLSs4xj6TmtaFo21iuy1h8CXiNJAGvA74eEfvKlDsG2BYj+wg2lCk3VscB7y75mS1K32fYxpLyrywp/1zg6KIylX7Gla7PfGAasLzonLen2yGp3a0BfihpnaRrx/9tWil3GE1eG0n+6zppnMdtAuZI6ilKDseSNJMMn/cE4MEy7/fziLio2oDLGG1o4I3AsZLaD5ccImK7pB8CryJpMvpapP+Opuf5aER89TCxPAUMknboptvKXZfDOeR7iojfStpP0kzymvRVzuPAbEnTi5LDseXOOUbD3/tHxxjvRpIaw1sqFT7Me5W7E+opoJ+kyfGx0p3p7+C7SRLYs4CfSbo7In5SRQyWco1h8vo9sDPtPOxS0ml8mqRzRjsoIjYCvwY+Jqkz7Yy8kqRPAJIOzI9IOkmJ0yXNBf4TOFnS6yR1pK9zJD1zjPE+SdIPMp7v73Hg45Kmp7GeP0r5fwdeD7w8XR52HfC36YfOcAfzK0sPjuQW0K8DH5XUo6SD+l3A8K2nNwDvkXR2el1OVPlO7CeBxWU6UG8GPgsMRUTZW44jYgOwDPiQpCmSngu8ZJTv+XD+Dbha0rPTmKdL+lNJPRXKfwV4iaQXp79PnUpuCFg4hvf6KvAiSa9KO7HnSjozIgppHP8s6QgASQskvThd/rP0WgrYSdLMmdlt1a3CiWGSSj/IXkLSebye5D+zG4CZYzj81STt1ZuAb5P0E/wo3fdpkg/IH5L8of4foCv9z+5PgMvT457gYMfxWHwQ+FLanPCqwxUu+v5OBB4Fekn6OSq5FTgJeDIi7is6z7fTOG+RtJOkJnRJhXO8g6Q9fB3wK5IEc2N6nm8AH0237QK+Q9KZWuob6dctku4p2v5l4LT062heQ9J/tBX4B5KEUpWIWAa8hSQhbSNpsnnDKOU3ApeRNEn2kdQC3ssYPmci4lGSfqV3p7GvILlxAZK+izXAb9OfwY9Jbm6A5Gf2Y5L+tN8An4/0mRCrng7WmM2sWaW3Z24GzoqIh7OOx1qbawxmE8NfAXc7KVgjuPPZrMlJeoTkTqiXZhuJTRZuSjIzsxHclGRmZiM4MZjVkaTXps9IHK5c3UZVrYaSsav+Mes4LBtODNY0dHC+gOFXSNpTtP68Ks55yPDjJfsvkFRIz79LyeB+b6wy/hGDDQJExFcj4k+qOZ9ZVtz5bE0jvZf9wDAYkgI4IyLW1PmtN0XEwvQhqcuA/5D0u4h46HAHDpOHnbYW4hqDTQiqYuhlSV8mGRLie2mN4H2jvUckvkPyMNep6VO+90raqWQY6A8WxTNcO7hS0qMkI7QODw++PX2/56hkciRJz9LBoceflPR3Fb7f0YavfoOScYF2KRku/bWjXLN/kbQpff2LpKnpvuFhy98tabOkxyvVlCQ9KOklResdSoYlP3O062kTlxODTRSfIBmu+UySp5kXkAwZDsnTsr0kA6sdSfLkbUTE60ieen5JJDOX/a/R3iBNJi8DZpEMHb6HZJiMWcCfAn8l6aUlh72AZHylF5MMwAcwK32/35Scv4fkKd3bSQaiOxE4ZEwfSQuA7wP/SPJ09HuAb0qaL2k68BngkojoAf47yVPC5fw9cB7JNTuDZCyi9xftP4rkSfcFJMOafE7S7DLnuRn4y6L1S4HHI6LS+9oE1xKJQdKN6X89pQO3VXOuC5VMTDL8GijzYWANlDbxvAX4fyNieGTX/0kyvAYkg9cdTTIU9mAkU3yO5z7sY5SM2vkUyTASr4uI1RFxZ0Q8EBGFiLgf+BpJIij2wYjYExH9Y3ifPwOeiIhPRcRAROyKiN+VKfeXwG0RcVv63j8iGQPp0nR/AThNUldEPB4RKyu832uBD0fE5ojoAz5EMjrrsMF0/2BE3EYyrES5qVG/AlwqaUa6/joOPzSHTWAtkRhIxpW/uBYnioifRcSZEXEm8EKSoYEPe1eJ1VW9h17elA5NPif92d8CoGTwuJ8pmZFtB3A1RTOipTYecrbKxjr0dsXhq9NRU/8ijeVxJTPenVLhPMdw6DDgxUNmbykZebbsUOcRsQm4C3i5pFkkY0UdbrRZm8BaIjFExC9IBt46QNIJkm5XMrXkL0f54xnNK4AfRMTemgRq1Soeenl4fomZEdENydDLEfHuiDieZOC8d0n64/TYp/ME57+TDK63KCJmkoy0qpIyUWG5nLEOvT08fPWsotf0iPg4QETckQ5ffjTwXySjj5aziSTJDDs23VaNL5HUZF4J/KbcENjWOloiMVRwPfCOiDibpI3281Wc43KS5gPL0NMcenm8w3UX6yGZlGhA0rlUngdhWB9JM0+l9/tP4ChJ/yPtGO6R9Owy5SoOXy3pSEn/T9rXsI+k+afSMNNfA96f9k3MI+mTqfZZie+QTPX5Tp7GiK02MbRkYpDUTdIp9w1JK0imkDw63ffn6V0Wpa87Ss5xNPDfSKYmtOxVO/Tyx0g+HLdLes843/NtwIcl7SL5UP36aIXTmuVHgbvS9zuvZP8uknmQX0Iy7PjDwIVlzjPa8NVtJJ3tm0hqyS9g5Mx8xf6RpG/ifpLO9HvSbeOW9qF8E1gCfKuac9jE0TJjJUlaDPxnRJyWdpKtjoijD3PYaOd7J0nTxVW1itFsIpP0AeDkiPjLwxa2Ca0lawwRsRNYr3SmLSXOOMxhpV6Nm5HMAJA0h+SW1uuzjsXqryUSg6SvkTQhPCN9aOdKklv1rpR0H7CSpGo+1vMtJrmD5Od1CNdsQpH0FpLmrB+kN3pYi2uZpiQzM6uNlqgxmJlZ7Uz4gb/mzZsXixcvzjoMM7MJZfny5U9FxPxy+yZ8Yli8eDHLli3LOgwzswlF0oZK+9yUZGZmIzgxmJnZCE4MZmY2ghODmZmN4MRgZmYjODGYmdkITgxmZjaCE4OZWRPYfddd7Fu3PuswgBZ4wM3MrBVsvPLNADzzv1ZlHIlrDGZmmYtCIesQRnBiMDPLWGHXrgPLQ1u2ZBhJwonBzCxj+W3bDiwPrFyZYSQJJwYzs4wNbXViMDOzIvntBxPDUN9TGUaSaFhikHSjpM2SHqywX5I+I2mNpPslndWo2MzMsnSgKUkiv3vX6IUboJE1hpuAi0fZfwlwUvq6CvhCA2IyM8vccGLoOHYRhV27M46mgYkhnUR86yhFLgNujsRvgVmSjm5MdGZm2Rnatg1NnUrH/CNG3KGUlWbqY1gAbCxa7023HULSVZKWSVrW19fXkODMzOolv207udmzaevpIb97EtUYxkBltkW5ghFxfUQsjYil8+eXnbLUzGzCyG/dSm7ObNp6ul1jKNELLCpaXwhsyigWM7OGyW/bRvusWeS6XWModSvw+vTupPOAHRHxeNZBmZnVW2HvHtq6e2jr6aGwezcRZRtLGqZhg+hJ+hpwATBPUi/wD0AHQERcB9wGXAqsAfYCb2xUbGZmWSr0D9DW1UmupxvyeWLvXjR9embxNCwxRMSrD7M/gL9uUDhmZk2jMDCAOrto6+4BIL97N20ZJoZmakoyM5uUor+fts5O2nq6ATLvgHZiMDPLUEQkNYauTnIzZgCQd2IwM5vEBgchn09qDN1pjSHjO5OcGMzMMlTYtw8AdXaS60n6GNyUZGY2iRX6+wFo6+yiLU0M+YzHS3JiMDPLUAwMACS3qx5oSnKNwcxs0ir0J4lBnV1o2jSQKOzZk2lMTgxmZhmKgbQpqasTSbRNm+bEYGY2mR2sMXQC0DZ9OnknBjOzyatwoMbQlXzt7qaw24nBzGzSOtD5XFRjcFOSmdkkVq4pyYnBzGwSi32uMZiZWZEDNYa0jyHXPd1DYpiZTWYHbledOjX56hqDmdnkVugfgI4O1NEBODGYmU16hYH+A/0LkCSGGByksH9/ZjE5MZiZZSj6B0oSQzpeUoa1BicGM7MMJZP0dB1YH57S04nBzGySioH+Ax3P4MRgZjbpFfpdYzAzsyKlnc+57jQxZPgsgxODmVmGYmAf6hp5VxK4xmBmNmklNYaipqR0Frd8hvM+OzGYmWXokNtVe2YAUMhw3mcnBjOzDCW3qxY3JU2Dtjbyh5n3Ob9zJ5HP1yUmJwYzswxF/8imJEm09fRQ2Dl6Ynj4BRew+Z8+VZeYnBjMzDISEUmNoXPqiO25nh7yu3ZWPK6wfz/R309u5sy6xOXEYGaWkRgchEJhRI0BoG3G6DWGwo4dAORmzqhLXA1NDJIulrRa0hpJ15bZP1PS9yTdJ2mlpDc2Mj4zs0aK/uH5njtHbM/1zBi1jyF/IDFM8BqDpBzwOeAS4FTg1ZJOLSn218BDEXEGcAHwKUlTGhWjmVkjFQaGp/UsqTH0dI9aYxhODG0zJnhiAM4F1kTEuojYD9wCXFZSJoAeSQK6ga3AUANjNDNrmEgTQ9kawyjPMeR3JP0PE77GACwANhat96bbin0WeCawCXgAeGdEFEpPJOkqScskLevr66tXvGZmdXWwxlCSGGb0UBg1MaRNSbMmfmJQmW1Rsv5iYAVwDHAm8FlJh/SuRMT1EbE0IpbOnz+/1nGamTXEwT6Gkqak7h4Ku3dXfE4hv2M7ALkZE7/zuRdYVLS+kKRmUOyNwLcisQZYD5zSoPjMzBpquMbQVqbGAJUH0ivs3Anp8w710MjEcDdwkqQlaYfy5cCtJWUeBf4YQNKRwDOAdQ2M0cysYQppjaG0KWl4WIxK/Qz57TvIzZiB2urzEd5el7OWERFDkt4O3AHkgBsjYqWkq9P91wEfAW6S9ABJ09M1EfFUo2I0M2ukqFBjaOtJp/eslBh27KCtTh3P0MDEABARtwG3lWy7rmh5E/AnjYzJzCwrhf6087mkjyE3XGPYUf7p5/zOnXW7Iwn85LOZWWYKA2nnc2kfw+xZwMG7j0rld+yoW8czODGYmWUm+ss/4JabPRuA/NYtZY/L79juGoOZWSsq7BvuYxg5iF57mhiGtmwte1x+67YDyaMenBjMzDIS/QOoowO1j+zuVUcHuZkzy9YYCgMDFHbtor2Oz3A5MZiZZSSZpKer7L7c3LkMbd12yPahp5IbNZ0YzMxaUAz00zZ1atl9uTmzyW85tMYwtDkZBqj9iPolhsPerirp2DGea3tEVJ5ZwszMRij0V64xtM+Zy761aw/ZPpSOD9c+b17d4hrLcwxfIhnTqNxYR8MCuAm4uQYxmZlNCoX+/kPGSRqWmzuH/O9/f8j2oafSxFDHpqTDJoaIuLB0m6SjIuKJ+oRkZjY5RP/eQ55hGNY+Zy75HTuIoaERndNDfX2QyzXlXUmvr2kUZmaTUKF/AE2rUGOYMxsiyG/fPmL7UF8f7XPmoFyubnFVmxguk/R2Sc+oaTRmZpNI0pQ0rey+9rlzARgq6YAeeuqpujYjQfWJ4c+BNcDLJN1Qw3jMzCaN6O+v2JTUcdRRAAxuGjk7wdDmvronhqoG0YuIJ4Hb05eZmVWh0N9fsSmpY+FCAAZ7HzuwLSIY3LCBaUuX1jWuqmoMkj4n6aZ02aOhmplVYbSmpNzcuairi8He3gPb8lu2UNi7lynHjvUpgupU25S0n4MT6LywRrGYmU0qhYGBik1JkuhYcAz7HzuYGPZv2ADAlMXH1TWuahPDXmCmpA6gvqnLzKwFxeAgDA7SVqEpCWDKgoUjmpL2b3g02d6kNYatwFrgc8BdtQvHzGxyODCtZ4UH3CDpZxjs7SUigLTGkMvRccwxdY1tXIlB0ixJXwRenm66GahvL4iZWQsanr2trXP0xFDYvfvAswz7N2ygY+EC1NFR19jGlRgiYjvwceBDwO+Ak4Bv1T4sM7PWFv17AUZtSup8xskADKx8KPn6wAN0nnxy3WOr5nbVK4H1EXEHsLzG8ZiZTQqFgfLzPRfrPP0MyOXov2c5U49fwuBjjzHnivoPPFFNYtgGXJ0+9XwfsCIi7q1tWGZmra2wd3i+58qJIdc9nc5TTmHv8nuYsmQJANPOOafusY07MUTExyT9BPgDcCbwfMCJwcxsHApjaEoC6DrrLLZ/4xuovZ22nh6mNqApadx3JUn6MHAZcBHwWET8a82jMjNrcZE2JVUadnvY7Fe9EiLYc9ddzHnDFXUdPG/YuBNDRHwA2Jce+3JJ/1bzqMzMWtxwU5JGaUoCmHrSSRzzT59k9mtezby3vrURoVU3VhJwI/BmYDrw+dqFY2Y2OYy1KQlgxkUXMeOii+od0gHVPuD2NyRJpR1wU5KZ2TiNtSkpC9UmhrVAJ/DdiHh+DeMxM5sUDjQltVBiWAn8FLhS0t01jMfMbFIoDPRDLlf3p5irUW0fwwkkzzNcn341M7NxKOzdS1tXF5KyDuUQ1dYYNkbErSSzuK0a60GSLpa0WtIaSddWKHOBpBWSVkr6eZXxmZk1tcKePbRNn551GGVVW2O4WNIfSEZX3UDSGT0qSbm0/EVAL3C3pFsj4qGiMrNI7nK6OCIelXRElfGZmTW1wp69TZsYqq0xzAKuAd5H8kzDWJwLrImIdRGxH7iF5EG5Yq8BvhURjwJExOYq4zMza2rNXGOoNjF8mOSOpNVAfozHLAA2Fq33ptuKnQzMlnSnpOWSyo4WJekqScskLevr6xtv7GZmmWuJxCDpjOHliOiNiB+ny2X7Csqdosy2KFlvB84G/hR4MfD/STpkYJCIuD4ilkbE0vnz54/x7c3Mmkdhzx7auid4YgDulXS/pPdJWlTFe/UCxcctBDaVKXN7ROyJiKeAXwBnYGbWYgq7d5Ob6DUG4FMkQ2B8HFgv6WeS3jSO4+8GTpK0RNIU4HLg1pIy3wWeJ6ld0jTg2Yzjriczs4miJZqSIuK9EXECyVSeN5AMt339OI4fAt4O3EHyYf/1iFgp6WpJV6dlVgG3A/cDvwduiIgHx/oeZmYTRTMnhjHfrippLvAy4BXAhSR9Bo+O580i4jbgtpJt15WsfxL45HjOa2Y2kcT+/cTg4MRPDMATJDWMbcAXga9ExK/qEpWZWQvL79kDQNu0iZ8Yvg18BfhBRAzWKR4zs5ZX2JMOud3dnXEk5Y05MUTEq+oZiJnZZFHYsxugaZuSqn3AzczMqlQYbkpqlcQg6SX1CMTMbLI4mBimZRxJedXUGD5a8yjMzCaRlqsxUH5oCzMzG6PhxNAKTz4PKx3fyMzMxuFAjaFJ70py57OZWYPld+0CWqspyczMnobCzl20TZuG2qudK62+qkkMT9Y8CjOzSSS/cydtM2dmHUZF404MEXFRPQIxM5ss8rt2kuvpyTqMityUZGbWYIUdO8nNmJF1GBU5MZiZNVh+1y7aWi0xSHpX0fIzaheOmVnry+/c0dQ1hnF1iUuaBfwzcIqkAZIJda4E3lj70MzMWlNhx07aZjRvH8O4EkNEbAfeKOnFwFPA6cC36hCXmVlLiqEhCnv2kJvRvHclVXsT7WBELJe0Cdhcy4DMzFrZ8MNtzdyUVG3n88WSFgLXkTQtmZnZGBSGn3pu4qakahPDLOAa4H3AvppFY2bW4vI7dgK0ZFPSh4FnRMRqSflaBmRm1sryO3cAkGviGkO1ieFvgenAT4Cf1S4cM7PWdrApqfX6GPYD69LlC2sUi5lZy8tvH64xtF5i2AvMlNQBHFvDeMzMWlp+21YAcnPmZBxJZdUmhn8A1gKfA75au3DMzFrb0NZttHV30zZlStahVFRtH8PfRMSnwUNimJmNR37r1qauLUB1Q2J8ATguHRLjPuDNeEgMM7MxyW/bSvvs2VmHMapxD4khqRf4BfA74Aw8JIaZ2ZgNbd1Gx9FHZx3GqKrpY9gCXA28Pl3vHeuBki6WtFrSGknXjlLuHEl5Sa+oIj4zs6aVNCW1UI0BICI+LumnwB+AM4HnAfce7jhJOZLO6otIksndkm6NiIfKlPsEcMd4YzMza2YRwdC2bbS3Uh8DgKQPAzlgBbAiIu4c46HnAmsiYl16nluAy4CHSsq9A/gmcM54YzMza2aFXbtgcJDc7OZODNXM+fwB4DPALuDlkv5tjIcuADYWrfem2w6QtAB4GcngfBVJukrSMknL+vr6xhy7mVmW8luHn2Fosaak1FuB/x0Rt4/jGJXZFiXr/wJcExF5qVzx9KCI64HrAZYuXVp6DjOzpjS0dRtA6zUlpW4E/krSdOCrEbFiDMf0AouK1hcCm0rKLAVuSZPCPOBSSUMR8Z0q4zQzaxr5rVsAWq8pKfU3JEmlnaRZaSzuBk6StETSFOBy4NbiAhGxJCIWR8Ri4D+AtzkpmFmrGEqbvtuPOCLjSEZXbWJYC3QC342I54/lgIgYAt5OcrfRKuDrEbFS0tWSrq4yDjOzCWPwySchl6N93tysQxlVtU1JK0k6kq+U9MmIGNMdRBFxG3BbybayHc0R8YYqYzMza0pDm/tonzcP5XJZhzKqahPDCcA2kg7gbbULx8ysdQ09+WTTNyNB9YlhY0T8VNLRwOZaBmRm1qqGNm+m47jmn6mg2j6GiyUtJHne4J9rGI+ZWcsa3LyZjglQY6g2McwCrgHeB+yrWTRmZi2qMDBAYccO2o84MutQDmvMiUHSGUWrHya5I2k1kK95VGZmLWZoc9Lq3n5kCyUG4F5J90t6H6CI+DFARFQcJdXMzBKDTzwBQMeRrdWU9ClgOvBxYL2kn0l6U33CMjNrLYO9jwHQsXBhxpEc3pgTQ0S8NyJOIBm24gbg+aTjFZmZ2egGe3uhra3pJ+mBcdyuKmkuycinrwAuJBkU79E6xWVm1lIGH+ul/agjUUdH1qEc1nieY3iCpIaxDfgi8JWI+FVdojIzazH7N/YyZUHzNyPB+BLDt4GvAD+IiME6xWNm1pIGe3uZfv75WYcxJodNDJKGH9N7T/r16ApzJWyPiJ21CszMrFUU9u1LnnpeuODwhZvAWGoMX+LghDqVZs8J4Cbg5hrEZGbWUgYfS+5ImjIB7kiCMSSGiLiwEYGYmbWqfWvXAjDl+OMzjmRsqh0Sw8zMxmj/2nUATFnixGBmZsD+9etoP+ooct3Tsw5lTJwYzMzqbN/adUydIM1I4MRgZlZXEcH+desmTP8CODGYmdXV4GObKOzdy9QTT8w6lDFzYjAzq6OBVQ8B0HnqMzOOZOycGMzM6mjfqlWQyzH15JOzDmXMnBjMzOpo4KFVTD1+CW2dnVmHMmZODGZmdTSwahVTnzlxmpHAicHMrG4GH3+coSefpOu0/5Z1KOPixGBmVid7l98DwLSlZ2ccyfg4MZiZ1Un/Pctpmz59QnU8gxODmVnd7F22nK4zz0Tt45n6JntODGZmdZDfsYN9Dz9M19lnZR3KuDkxmJnVQf+KFRDBtLOXZh3KuDU0MUi6WNJqSWskXVtm/2sl3Z++fi3pjEbGZ2ZWK3uXLYf2drpOn1h3JEEDE4OkHPA54BLgVODVkk4tKbYeeEFEnA58BLi+UfGZmdXS3nvuofPUU2nr6so6lHFrZI3hXGBNRKyLiP3ALcBlxQUi4tcRsS1d/S0wMebBMzMrUti3j4H772fa0onXjASNTQwLgI1F673ptkquBH5QboekqyQtk7Ssr6+vhiGamT19A/ffTwwOOjGMgcpsi7IFpQtJEsM15fZHxPURsTQils6fP7+GIZqZPX17ly0DiWkT8I4kgEbeXNsLLCpaXwhsKi0k6XTgBuCSiNjSoNjMzGpm793LmHryyeRmzsw6lKo0ssZwN3CSpCWSpgCXA7cWF5B0LPAt4HUR8YcGxmZmVhMxOMjeFSsmbDMSNLDGEBFDkt4O3AHkgBsjYqWkq9P91wEfAOYCn5cEMBQRE/fqmtmkM7BqFbF374QbH6lYQ5/TjojbgNtKtl1XtPxm4M2NjMnMrJb2LlsOQNfZEzcx+MlnM7Ma6r/3HjqOPZaOI47IOpSqOTGYmdXQwMqH6DrttKzDeFqcGMzMamRo2zYGN22i81mlgzpMLE4MZmY1sm/VKgA6T3ViMDMzYOChhwDonGBzPJdyYjAzq5F9Dz9M+5FHkps1K+tQnhYnBjOzGtm3bj1TTzg+6zCeNicGM7MaiAj2r1/PlMVLsg7laXNiMDOrgaG+Pgq7dzPleNcYzMwM2L/+EQCmLFmcaRy14MRgZlYD+9evA2CqawxmZgawf/161NVF+5FHZh3K0+bEYGZWA/vWrWfKksWobeJ/rE7878DMrMFicJCBhx4i9u8/sG3/+vVMbYE7ksCJwcxsXPb87veseeEfs/7PX07vu95F5PMUBgYYfOyxlrgjCZwYzMzG5fH3vx91dTHniivY/eOfsOO7t7J/w6MQ0RJ3JIETg5nZuAxt2ULPC1/IEddeQ8eiRey87Tb2r1sLtMYdSdDgGdzMzCayyOeJvXtp6+5GEjMuvpgtN95IbuZM1NXFlBNOyDrEmnCNwcxsjAp79gCQ6+kGYMafXgr5PDu//32mnXsObVOmZBlezTgxmJmNUWH3bgDaupPE0HnKKUw77zwAus8/P7O4as2JwcxsjPK70sQwvfvAtnlv+ytys2fT/cIXZhVWzbmPwcxsjAp7RtYYAKafey4n/+bXWYVUF64xmJmN0XBTUq57esaR1JcTg5nZGJX2MbQqJwYzszE60MfQ05NxJPXlxGBmNkYHagzTXWMwMzPSzmeJtmldWYdSV04MZmZjlN+9m7bp01tiaO3RtPZ3Z2ZWQ4Vdu1u+4xkanBgkXSxptaQ1kq4ts1+SPpPuv1/SWY2Mz8xsNIXduw8Mh9HKGpYYJOWAzwGXAKcCr5Z0akmxS4CT0tdVwBcaFZ+Z2eEU9uxu+Y5naOyTz+cCayJiHYCkW4DLgIeKylwG3BwRAfxW0ixJR0fE47UO5v9/zZvIF/K1Pq2ZtbwpcPkVWQcBQK4txzv+/caan7eRTUkLgI1F673ptvGWQdJVkpZJWtbX11fzQM3MJrNG1hhUZltUUYaIuB64HmDp0qWH7B+LemRZM7NW0MgaQy+wqGh9IbCpijJmZlZHjUwMdwMnSVoiaQpwOXBrSZlbgdendyedB+yoR/+CmZlV1rCmpIgYkvR24A4gB9wYESslXZ3uvw64DbgUWAPsBd7YqPjMzCzR0PkYIuI2kg//4m3XFS0H8NeNjMnMzEbyk89mZjaCE4OZmY3gxGBmZiM4MZiZ2QhK+nsnLkl9wIYqD58HPFXDcOrFcdbORIgRHGctTYQYofFxHhcR88vtmPCJ4emQtCwilmYdx+E4ztqZCDGC46yliRAjNFecbkoyM7MRnBjMzGyEyZ4Yrs86gDFynLUzEWIEx1lLEyFGaKI4J3Ufg5mZHWqy1xjMzKyEE4OZmY0waRODpIslrZa0RtK1GcfyiKQHJK2QtCzdNkfSjyQ9nH6dXVT+b9O4V0t6cR3julHSZkkPFm0bd1ySzk6/vzWSPiOp3IRMtY7zg5IeS6/pCkmXZhmnpEWSfiZplaSVkt6Zbm+q6zlKnE1zPSV1Svq9pPvSGD+Ubm+2a1kpzqa5lhVFxKR7kQz7vRY4HpgC3AecmmE8jwDzSrb9L+DadPla4BPp8qlpvFOBJen3katTXM8HzgIefDpxAb8HnkMyQ98PgEsaEOcHgfeUKZtJnMDRwFnpcg/whzSWprqeo8TZNNczPV93utwB/A44rwmvZaU4m+ZaVnpN1hrDucCaiFgXEfuBW4DLMo6p1GXAl9LlLwEvLdp+S0Tsi4j1JHNXnFuPACLiF8DWpxOXpKOBGRHxm0h+w28uOqaecVaSSZwR8XhE3JMu7wJWkcxn3lTXc5Q4K2l4nJHYna52pK+g+a5lpTgryexvqNRkTQwLgI1F672M/stfbwH8UNJySVel246MdPa69OsR6fasYx9vXAvS5dLtjfB2SfenTU3DzQqZxylpMfBHJP9BNu31LIkTmuh6SspJWgFsBn4UEU15LSvECU10LcuZrImhXPtclvftnh8RZwGXAH8t6fmjlG222IdViiureL8AnACcCTwOfCrdnmmckrqBbwL/IyJ2jla0QjxZxdlU1zMi8hFxJsm88OdKOm2U4pldywpxNtW1LGeyJoZeYFHR+kJgU0axEBGb0q+bgW+TNA09mVYhSb9uTotnHft44+pNl0u311VEPJn+URaAf+Ngc1tmcUrqIPmw/WpEfCvd3HTXs1yczXg907i2A3cCF9OE17JcnM16LYtN1sRwN3CSpCWSpgCXA7dmEYik6ZJ6hpeBPwEeTOO5Ii12BfDddPlW4HJJUyUtAU4i6ZhqlHHFlVbpd0k6L72T4vVFx9TN8AdE6mUk1zSzONNz/h9gVUR8umhXU13PSnE20/WUNF/SrHS5C3gR8F8037UsG2czXcuK6tmz3cwv4FKSOy7WAn+fYRzHk9yJcB+wcjgWYC7wE+Dh9OucomP+Po17NXW8OwH4GklVd5Dkv5Yrq4kLWEryy78W+CzpE/d1jvPLwAPA/SR/cEdnGSfwXJLq//3AivR1abNdz1HibJrrCZwO3JvG8iDwgWr/Zup8LSvF2TTXstLLQ2KYmdkIk7UpyczMKnBiMDOzEZwYzMxsBCcGMzMbwYnBzMxGcGIwKyJplqS3Fa0fI+k/6vReL5X0gQr7dqdf50u6vR7vb1aJE4PZSLOAA4khIjZFxCvq9F7vAz4/WoGI6AMel3R+nWIwO4QTg9lIHwdOSMfJ/6SkxUrneZD0BknfkfQ9SeslvV3SuyTdK+m3kuak5U6QdHs6KOIvJZ1S+iaSTgb2RcRT6foSSb+RdLekj5QU/w7w2rp+12ZFnBjMRroWWBsRZ0bEe8vsPw14Dcn4Nh8F9kbEHwG/IRmqAJJJ3d8REWcD76F8reB84J6i9X8FvhAR5wBPlJRdBjyvyu/HbNzasw7AbIL5WSTzFOyStAP4Xrr9AeD0dFTS/w58o2iSrallznM00Fe0fj7w8nT5y8AnivZtBo6pTfhmh+fEYDY++4qWC0XrBZK/pzZgeyRDLY+mH5hZsq3S+DSdaXmzhnBTktlIu0imtKxKJHMXrJf0SkhGK5V0Rpmiq4ATi9bvIhnlFw7tTziZgyNwmtWdE4NZkYjYAtwl6UFJn6zyNK8FrpQ0PGJuuWljfwH8kQ62N72TZJKmuzm0JnEh8P0qYzEbN4+uapYRSf8KfC8ifnyYcr8ALouIbY2JzCY71xjMsvM/gWmjFZA0H/i0k4I1kmsMZmY2gmsMZmY2ghODmZmN4MRgZmYjODGYmdkITgxmZjbC/wWUDs+mQLA8xQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "swiftdiff['vmag'].sel(id=tpidx).plot.line(ax=ax, x=\"time (d)\")\n", + "ax.set_ylabel(\"$|\\mathbf{v}_{swiftest} - \\mathbf{v}_{swifter}|$\")\n", + "ax.set_title(\"Heliocentric velocity differences \\n Test Particles only\")\n", + "legend = ax.legend()\n", + "legend.remove()\n", + "fig.savefig(\"rmvs_swifter_comparison-mars_ejecta-testparticles-vmag.png\", facecolor='white', transparent=False, dpi=300)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'rmag' (time (d): 333)>\n",
+       "array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "...\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.13180114e-12,\n",
+       "       6.30252092e-12, 1.12657932e-11, 1.70947866e-11, 2.35410127e-11,\n",
+       "       3.01486367e-11, 3.63634702e-11, 4.16224366e-11, 4.54289913e-11,\n",
+       "       4.74142910e-11, 4.73824194e-11, 4.53327404e-11, 4.14594589e-11,\n",
+       "       3.61300773e-11, 2.98446324e-11, 2.31845539e-11, 1.67548923e-11,\n",
+       "       1.11262399e-11, 6.78147816e-12, 4.07218435e-12, 3.25977426e-12,\n",
+       "       4.52137637e-12, 7.66342713e-12, 1.23344633e-11, 1.81013732e-11,\n",
+       "       2.44264806e-11, 3.07065663e-11, 3.63320360e-11, 4.07478190e-11,\n",
+       "       4.35128453e-11, 4.43475549e-11, 4.31649567e-11, 4.00801554e-11,\n",
+       "       3.53984592e-11, 2.95862328e-11, 2.32329074e-11, 1.70175537e-11,\n",
+       "       1.17040422e-11])\n",
+       "Coordinates:\n",
+       "    id        int64 2\n",
+       "  * time (d)  (time (d)) float64 0.0 11.0 22.0 ... 3.63e+03 3.641e+03 3.652e+03
" + ], + "text/plain": [ + "\n", + "array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + "...\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.13180114e-12,\n", + " 6.30252092e-12, 1.12657932e-11, 1.70947866e-11, 2.35410127e-11,\n", + " 3.01486367e-11, 3.63634702e-11, 4.16224366e-11, 4.54289913e-11,\n", + " 4.74142910e-11, 4.73824194e-11, 4.53327404e-11, 4.14594589e-11,\n", + " 3.61300773e-11, 2.98446324e-11, 2.31845539e-11, 1.67548923e-11,\n", + " 1.11262399e-11, 6.78147816e-12, 4.07218435e-12, 3.25977426e-12,\n", + " 4.52137637e-12, 7.66342713e-12, 1.23344633e-11, 1.81013732e-11,\n", + " 2.44264806e-11, 3.07065663e-11, 3.63320360e-11, 4.07478190e-11,\n", + " 4.35128453e-11, 4.43475549e-11, 4.31649567e-11, 4.00801554e-11,\n", + " 3.53984592e-11, 2.95862328e-11, 2.32329074e-11, 1.70175537e-11,\n", + " 1.17040422e-11])\n", + "Coordinates:\n", + " id int64 2\n", + " * time (d) (time (d)) float64 0.0 11.0 22.0 ... 3.63e+03 3.641e+03 3.652e+03" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "swiftdiff['rmag'].sel(id=2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "swiftestOOF", + "language": "python", + "name": "swiftestoof" + }, + "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.7.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/tp.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/tp.in new file mode 100644 index 000000000..c7cf002d6 --- /dev/null +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/tp.in @@ -0,0 +1,49 @@ +16 +101 +0.33208578766229190915 0.07439013071780828379 -0.02438290851908785084 +-0.008988542188201206762 0.028710618792657981169 0.0034094833969203438596 +102 +0.33203966624962866216 0.07434400930514498129 -0.02438290851908785084 +-0.014195641132219511543 0.028710618792657981169 0.0034094833969203438596 +103 +-0.7187543234391324809 -0.011798260816488121555 0.041316403191083782287 +0.0065615071841567274707 -0.020313576971905909774 -0.00029114855617710840843 +104 +-0.71886874402007694407 -0.011912681397432518804 0.041316403191083782287 +-0.006132960226534060408 -0.020313576971905909774 -0.00029114855617710840843 +105 +0.35683111163121072895 -0.9518327808922094624 4.4027442504036787155e-05 +0.022724479262608666269 0.0059737936889703449964 -3.3484113013969089573e-07 +106 +0.3567106558193317012 -0.95195323670408849015 4.4027442504036787155e-05 +0.008935598794060913702 0.0059737936889703449964 -3.3484113013969089573e-07 +107 +-1.5233391647104730371 0.6724145771476651712 0.051459143378398922164 +-0.0020480822268840624331 -0.011607719813367209372 -0.000117479966462153095864 +108 +-1.5234032495379807859 0.6723504923201574224 0.051459143378398922164 +-0.008207040423331847523 -0.011607719813367209372 -0.000117479966462153095864 +109 +4.050605826355517358 -2.9904269687677218492 -0.078187280837353656526 +0.041279424970441319642 0.006432188574295680597 -0.00012509257442073270106 +110 +4.049284028339322994 -2.9917487667839162135 -0.078187280837353656526 +-0.032485009432853539924 0.006432188574295680597 -0.00012509257442073270106 +111 +6.299479995832536261 -7.7058625321556393217 -0.11669919842191249504 +0.02612723553831041573 0.0035242303011843410798 -0.00022097170940726839814 +112 +6.2983790111222752728 -7.70696351686590031 -0.11669919842191249504 +-0.01809910222875676239 0.0035242303011843410798 -0.00022097170940726839814 +113 +14.856321905516212567 13.007829033301401722 -0.14417795763685259391 +0.010478935887110856981 0.0027821364817078499815 4.40781085949555924e-05 +114 +14.855842389541807691 13.007349517326996846 -0.14417795763685259391 +-0.015710591190212928187 0.0027821364817078499815 4.40781085949555924e-05 +115 +29.55768244045575699 -4.6291447957067299868 -0.58590957207831262377 +0.014905509815736753265 0.0031274056019462009859 -7.51415892482447254e-05 +116 +29.557216915563323312 -4.6296103205991601115 -0.58590957207831262377 +-0.0139657618108195089035 0.0031274056019462009859 -7.51415892482447254e-05