diff --git a/Makefile.Defines b/Makefile.Defines index ed862668b..6a52a013f 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -51,7 +51,7 @@ ADVIXE_FLAGS = -g -O2 -qopt-report=5 -vec -vecabi=cmdtarget -simd -shared-intel VTUNE_FLAGS = -g -O2 -vec -simd -shared-intel -qopenmp -debug inline-debug-info -parallel-source-info=2 -parallel -DTBB_DEBUG -DTBB_USE_THREADING_TOOLS -qopenmp -fp-model no-except -mp1 -xhost -traceback #Be sure to set the environment variable KMP_FORKJOIN_FRAMES=1 for OpenMP debuging in vtune -IDEBUG = -O0 -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin +IDEBUG = -O0 -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin -init=snan,arrays STRICTREAL = -mp1 -fp-model strict -prec-div -prec-sqrt -assume protect-parens SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -prec-div -prec-sqrt -assume protect-parens @@ -70,7 +70,7 @@ GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporari #FFLAGS = -init=snan,arrays -traceback -no-wrap-margin -O3 -g -CB -nogen-interfaces -no-pie -fp-speculation=safe $(SIMDVEC) $(PAR) #$(HEAPARR) #FFLAGS = $(IDEBUG) $(HEAPARR) -FFLAGS = -init=snan,arrays -no-wrap-margin -fp-model strict -fp-model no-except -traceback -g $(OPTIMIZE) -O3 $(SIMDVEC) $(PAR) $(HEAPARR) +FFLAGS = -init=snan,arrays -no-wrap-margin -fp-model strict -fp-model no-except -traceback -g -O3 $(SIMDVEC) $(PAR) $(HEAPARR) FORTRAN = ifort #AR = xiar diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/.ipynb_checkpoints/swiftest_rmvs_vs_swifter_rmvs-checkpoint.ipynb b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/.ipynb_checkpoints/swiftest_rmvs_vs_swifter_rmvs-checkpoint.ipynb new file mode 100644 index 000000000..105fe9dc8 --- /dev/null +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/.ipynb_checkpoints/swiftest_rmvs_vs_swifter_rmvs-checkpoint.ipynb @@ -0,0 +1,249 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import xarray as xr\n", + "import swiftestio as swio\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" + ] + } + ], + "source": [ + "inparfile = 'param.swifter.in'\n", + "param = swio.read_swifter_param(inparfile)\n", + "swifterdat = swio.swifter2xr(param)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file config.swiftest.in\n" + ] + } + ], + "source": [ + "config_file_name = 'config.swiftest.in'\n", + "config = swio.read_swiftest_config(config_file_name)\n", + "swiftestdat = swio.swiftest2xr(config)" + ] + }, + { + "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": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAElCAYAAADgCEWlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABRb0lEQVR4nO3dd3hUZfbA8e9JD4TeCYSE3nsRpIllEWQt2LuLoq5ucXXV1V3b/lxdt4mKfe2u2BUVVKQIgvQOoYSaQCAhkIQUUs/vj3vREJKQhGlJzud55snM3Dv3nrlJ5sz73vueV1QVY4wx5rggfwdgjDEmsFhiMMYYcwJLDMYYY05gicEYY8wJLDEYY4w5gSUGY4wxJ7DEYMokIo+IyDvu/RgRyRKRYH/HVRERGSUiW/0dB5w6Fl8eUxFZICI3u/evEZFvSyw7U0S2u7FcJCKtRGShiBwVkX95OzYTmCwx1FIisltEzin13I0i8kNVt6Wqe1U1SlWLPBdh1YiIikjnitZR1UWq2s1XMVWkdCylfx/+Oqaq+q6qnlfiqceA59xYPgOmAoeAhqp6ty9jM4HDEoOpFUQkxN8x1FAdgE2lHm/Waox8td9B7WGJoQ4TkbYi8rGIpIrILhH5bTnrxbrf2ENKvG6miBwWkQQRuaXEusEi8oCI7HC7I1aJSHt3WXcRmeO+bquIXF7idW+IyHQR+cp93TIR6eQuW+iuts7t8rhCRMaKSJKI3CciB4DXjz9XYpvtReQT9/2lichz5by/R0TkIxF53933ahHpV2J5D7c7Jl1ENonIL0ssmyAim93X7RORe9znf4pFRN4GYoAv3PjvreIxfUREPhCRt9z9bBKRwRX8Xs8VkS0ikuG+Zymx7KdWo4jsADqWiOs94AbgXvfxOSISJCL3u7/PNDeOpqX+LqaIyF5gnvv8r0QkXkSOiMg3ItKhxP5VRG5zu6+OuL/zkvHd4r72qHtcB5Y4PmX+rYrIUBFZKSKZInJQRP5d3rExlaSqdquFN2A3cE6p524EfnDvBwGrgIeAMJwPiJ3AL9zljwDvuPdjAQVC3MffA88DEUB/IBU42132R2AD0A3nA6kf0AyoDyQCNwEhwECcLote7uveAA4DQ93l7wIzSsSuQOcSj8cChcDfgXAg0n0uyV0eDKwD/uPuOwIYWc6xegQoAC4FQoF7gF3u/VAgAXjAPU7jgKNAN/e1ycAo934TYGCJ+JLK+31U8Zg+AhwDJrjv6wlgaTnvpTmQWeK93OUep5tL/w2UE9cbwP+VePx7YCnQzj3OLwHvlXoPb7nHOBK4yD1ePdzf45+BJaV+j18CjXGSZSow3l12GbAPGILzt9MZpwVzqr/VH4Hr3PtRwBn+/v+r6Te/B2A3L/1inX/4LCC9xC2HnxPDMGBvqdf8CXjdvf8IZSQGoD1QBDQo8bongDfc+1uBC8uI5wpgUannXgIedu+/AbxaYtkEYEuJx2UlhnwgotRzxxPDcPdDJ6QSx+oRSnzQuh9EycAo93YACCqx/D3gEff+XuBWnD55yoqlxO+jzMRQiWP6CPBdiWU9gdxy3sv1pd6LAElUPzHE4yYo93EbnCQaUuI9dCyxfDYwpdSxzAE6lPg9jiyx/APgfvf+N8DvynhPp/pbXQg8CjT39/9dbblZV1LtdpGqNj5+A35dYlkHoK3bPZIuIuk434pbnWKbbYHDqnq0xHN7gGj3fntgRxmv6wAMK7W/a4DWJdY5UOJ+Ds63v4qkquqxcpa1B/aoauEptnFc4vE7qlqM82Ha1r0lus8dV/L9TsZJYntE5HsRGV7J/ZV0qmMKJx+bCCm7T79tqfeiJR9XQwfg0xK/s3icJFby7ySx1PrTSqx/GCc5VfRejv+eK/rbqehvdQrQFdgiIitE5IIqv0tzAjtZVHclArtUtUsVX7cfaCoiDUp8kMXgdAEc324nYGMZ+/teVc+tbsBlqOgEaSIQIyIhlUwO7Y/fEZEgnK6T/ceXiUhQieQQA2wDUNUVwIUiEgrcifMN+KdtVTLWUx3Tqkgu9V6knHgqKxH4laouLr1ARGLdu1pq/cdV9d1q7qtTOc+X+7eqqtuBq9zf2yXARyLSTFWzqxGDwU4+12XLgUz35G2kOCeNe4vIkIpepKqJwBLgCRGJEJG+ON/Yjn8QvAr8VUS6iKOviDTD6VfuKiLXiUioexsiIj0qGe9BnL7lqry/ZOBJEanvxnpmBesPEpFL3G/hvwfycPrWlwHZOCdkQ0VkLDAJmCEiYeKMC2ikqgU4ffvlXX5abvyVOKZV8RXQq8R7+S0ntsqq6kXg8eMnkEWkhYhceIr1/yQivdz1G4nIZZXc16vAPSIyyP3b6ezut8K/VRG5VkRauIk73d2W3y6trg0sMdRR6lw/PwnnROcunBPBrwKNKvHyq3D6l/cDn+KcJ5jjLvs3zrfmb3E+KP8LRLrfhM8DrnRfd4CfTxxXxiPAm25XwuWnWrnE++uMcx4gCec8R3k+d5cfAa4DLlHVAlXNB34JnI9zjJ4HrlfVLe7rrgN2i0gmcBtwbTnbfwL4sxv/PWUsr+iYVpqqHsI5ifskkAZ0AU76tl8F04CZwLcichQnWQ6rYP+f4vxeZ7jHZCPOsatM7B8CjwP/wznB/xnQtBJ/q+OBTSKS5cZ7ZQVdjKYSxD15Y0ydJSKP4JzYLu9D3Zg6xVoMxhhjTmCJwRhjzAmsK8kYY8wJrMVgjDHmBJYYTJ0jZVSerS2kVA0mY6rDEoOpldwPx2xxisHtE5F/i4/nk5BKlAo3JhBZYjC1WT9VjQLOBq4GbjnF+sYYLDGYOsAdjLYI6F16mVuy+Ud34FmyiDwnImEllp+qTHSZJaal7FLhzUXkS3dfh0VkkVvG4SQiMsKt+5Ph/hxRYtkCEfmriCwWpzz1tyLSvIxtXCYiq0o9d7eIfFa1I2jqGksMptYTkZ44VVLXlLG4CKc0dXOciqxnc2KxQYALcEpB9wMuB37hbvcinGJulwAtcJLPewCqOtp9bT91Zkd7H7gbZwR2C5wCcA9QRg0lceY7+Ap4Bqdk+b+Br9zSIsddjVPCvCVOKeqyRlPPBOJKlR25Fni7jHWN+UmtSAwi8pqIpIhI6cJt1d3e1+63ui/LWf6sO/zeBLbVInIE+AKnhMLrpVdQ1VWqulRVC1V1N04p8DGlVntSVdNVdS8wH6c0Azjltp9Q1Xi3UN/fgP5SYmKaUgpwylZ3cMttLNKyrxefCGxX1bfduN4DtuCUhTjudVXdpqq5OCVI+pfeiKrmAe/jlulw6xfF4tStMqZctSIx4NSQH+/B7f0DpwbOScSZOauxB/dlvGegqjZR1U6q+udSpbMBEJGubvfOAbe2z99wWg8llVcmujIlpkv6B84kNt+KyE4Rub+c9drilN0u6VRluMsrUf4mcLXb/XUd8IGbMIwpV61IDKq6EOef8ici0sn95r/K7cvtXoXtzcUp4nUC96qWfwD3nm7MJmC8gPNtvIuqNsTp3pGKX/KTRODWknNeqGqkqi4pa2VVPaqqd6tqR5xv/38QkbPLWHU/TtIpqVpluFV1Kc6ERqNwup+sG8mcUq1IDOV4GfiNqg7C6X993gPbvBOYqarJHtiWCQwNcKrAZrlfHm6vwmtPVWL6hFLbInKBW0pa+LlEd1nloWfhlCi/WkRCROQKnFnbqtsF9BbwHFCoqj9UcxumDqmVg2BEJAoYAXxY4gKScHfZJcBjZbxsn6r+ooJttsUpZzzWo8Eaf7sH50vEvTgnp9/Hmdf5lFT1U/dvbYZ7XiEDmAN86K7yCE6p8EhgKk5X0HM4J5+PAM+r6oIytpsmzixk03BaNAnABW5J7ep4G/irezPmlGpNrSRxZpP6UlV7i0hDYKuqtjmN7Y0F7lHVC9zHE3HmFjhe5z0G2KmqNoDJBDQ3MaXgnHPZ7u94TOCrlV1JqpoJ7DrerBdHv9Pc5leq2lpVY1U1FsixpGBqiNuBFZYUTGXViq4kEXkPp4unuYgkAQ/jTDT/goj8GQgFZgDrKrm9RUB3IMrd3hRV/cYbsRvjTSKyG+dk+kX+jcTUJLWmK8kYY4xn1MquJGOMMdVX47uSmjdvrrGxsf4OwxhjapRVq1YdUtUWZS2r8YkhNjaWlStX+jsMY4ypUUSk9Oj6n1hXkjHGmBNYYjDGGHMCSwzGGGNOUOPPMZSloKCApKQkjh07duqV/SgiIoJ27doRGhrq71CMMeYntTIxJCUl0aBBA2JjYylRKymgqCppaWkkJSURFxfn73CMMeYntbIr6dixYzRr1ixgkwKAiNCsWbOAb9UYY+oenycGEQkWkTVlzY7m1jR6RkQSRGS9iAw8jf2cXqA+UBNiNMbUPf5oMfwOiC9n2flAF/c2FafksDHGVNu6xHSem7ed91fsJTuv0N/heMzihENkHivwyrZ9mhhEpB3OfLavlrPKhcBb6lgKNBaRapfOPl0jRowo8/kbb7yRjz76yMfRGGOqIjuvkLveX8uF0xfzz2+3cd/HGxjzjwXM23LQ36GdtozcAn71xgr+8fVWr2zf1y2Gp3EmRDlp7l1XNM50icclUcb8uSIyVURWisjK1NRUjwd53JIlZc7QaIwJcElHcpj8whI+X7uP357dhXUPn8fHtw+nZYNwbnlrFV9vPHDqjQSwmev2k1dYzOWD23tl+z5LDO6MVCmquqqi1cp47qTyr6r6sqoOVtXBLVqUWerDI6Kioo7vjzvvvJOePXsyceJEUlJSvLZPY8zpOVZQxE2vr2Bfei5v3DSUP5zblUaRoQzq0JQPbhtOn+hG/OGDtWw/eNK07jXG+yv20rNNQ3pHN/TK9n3ZYjgT+KVbH34GME5E3im1ThJQMgW2w5kY3a8+/fRTtm7dyoYNG3jllVesJWFMAPvXt1vZnpLFc1cPZHTXE784RoWH8PJ1gwgPCeL+TzZQXFzzph3YuC+DjfsyuXJoe69dwOKzxKCqf1LVdu7sZ1cC81T12lKrzQSud69OOgPIUNVkX8VYnoULF3LVVVcRHBxM27ZtGTeuUlMCG2N8bNWeI7z6wy6uPSOGMV3L7k1o2TCCByf2ZNWeI7y3Yq+PIzx97y7bS3hIEBf2O6mX3WP8Po5BRG4Tkdvch7OAnTiTn78C/NpvgZVil5YaE/j++c1WmtUP50/n96hwvckDoxnRqRlPzt5Cek6+j6I7ffvTc/l4VRKTB7WjUT3vVUzwS2JQ1QWqeoF7/0VVfdG9r6p6h6p2UtU+qhoQ9bRHjx7NjBkzKCoqIjk5mfnz5/s7JGNMKT/uSOPHnWncPrYT9cMrLuogIjw8qRdZeYW8+P1OH0V4+l5YsANFueMs70437/cWQ01w8cUX06VLF/r06cPtt9/OmDFj/B2SMaaU/3y3jZYNwrlmWEyl1u/WugEX9Y/mjSW7SMkM/AoE+9NzeX9FIpcPbk9040iv7qtW1krylKysLMD5dvHcc8/5ORpjTHk27c9g+a7D/HliDyJCgyv9ut+f04Uv1u1n+vwEHr2wtxcjPH3HWwu/9nJrAazFYIypBWYsTyQsJIhLB7Wr0us6NKvP5IHtmLEikUNZeV6K7vQdby1c5oPWAlhiMMbUcLn5RXy2dh8Teremcb2wKr9+6piO5BcV8+aS3Z4PzkN+ai2M7eST/VliMMbUaLM2JHP0WCFXDq3cuYXSOrWI4hc9W/Pmkt1kBWAtpZKthXZN6vlkn5YYjDE12owVe4lrXp9hcU2rvY3bxnYi81ghM5YH3rgGX7cWwBKDMaYGS0g5yordR7hiyOmNAu7fvjFndGzKq4t2kV9YXik330vOcFoLlw7yXWsBLDEYY2qwGcsTCQkSJg+s2knnstw+tjMHMo/x2dp9HojMM15YsINiVe44y3etBbDE4DWJiYmcddZZ9OjRg169ejFt2jR/h2RMrZJXWMQna/Zxbs9WtGgQftrbG92lOT3bNOSl73cERA2l5IxcZiz37bmF4ywxeElISAj/+te/iI+PZ+nSpUyfPp3Nmzf7Oyxjao05mw9yODu/2iedSxMRbh3TkR2p2cyJ9/+cDU/P2e6OcvZtawEsMXhNmzZtGDjQmZm0QYMG9OjRg337AqeJakxNN2N5ItGNIxnVubnHtjmxTxvaN43k+QU7UPVfq2FDUgYfrErkpjPjfN5agDow8vnRLzaxeX+mR7fZs21DHp7Uq9Lr7969mzVr1jBs2DCPxmFMXbU3LYcfEg5x1zldCQryXIHLkOAg7hjbmfs/2cCczQc5r1drj227slSVR7/YRLP6Ydw5zvujnMtiLQYvy8rKYvLkyTz99NM0bOidSTWMqWs+WJlIkMDlQ07/pHNplw5qR1zz+jz93Xa/tBreX5HIyj1H+OMvutEwwnsVVCtS61sMVflm72kFBQVMnjyZa665hksuucRvcRhTmxQWFfPhqkTGdmtJm0aeLw8REhzE7WM7ce9H61mwLZWzurX0+D7Ksz89l//7Kp7hHZtx2SDvTNtZGdZi8BJVZcqUKfTo0YM//OEP/g7HmFpj/tZUDmbmceUQ731wXtQ/mjaNInhh/g6v7aM0VeVPn2ygqFj5++S+Hu0iqypfzvkcISLLRWSdiGwSkUfLWGesiGSIyFr39pCv4vO0xYsX8/bbbzNv3jz69+9P//79mTVrlr/DMqbGm7F8Ly0ahHNWd+99kw8LCWLq6I4s332YpTvTvLafkj5clcT321K5b3w3Ypr5/oRzSb7sSsoDxqlqloiEAj+IyGxVXVpqvUXHJ/GpyUaOHOnXqxqMqY0OZBxj/tYUbhvTidBg736vvXJIDK8s3MljX2zmi9+MJNiL3+C3HMjk4c83MSyuKdcPj/XafirLl3M+q6pmuQ9D3Zt9chpjKu3DlYkUK1zhxW6k4yLDgnlgYg82J2cyw4tzQ2fkFnDb26toEBHCs1cP8GsX0nE+PccgIsEishZIAeao6rIyVhvudjfNFhH/nTk2xgSU4mLl/ZWJnNm5GR2a1ffJPif2acMZHZvyj2+2emWWt2MFRdzy1kr2pecy/ZqBtGwQ4fF9VIdPE4OqFqlqf6AdMFRESk+ZtBrooKr9gGeBz8rajohMFZGVIrIyNTXVmyEbYwLEDwmHSDqSy5VDPDPSuTJEhP+7qA/HCoq4+8N1Hi2VUVSs/H7GWpbvOsy/Lu/PkNjqV4f1NL9claSq6cACYHyp5zOPdzep6iwgVEROGtaoqi+r6mBVHdyiRQsfRGyM8bcZK/bSpF4o5/Vq5dP9dm4ZxV8u6Mmi7Yf47w+7PLJNVeWhzzfy9aYD/OWCnvyyX1uPbNdTfHlVUgsRaezejwTOAbaUWqe1uLVzRWSoG59vLgkwxgSsQ1l5zNl8kEsGtiM8pPJzOnvK1UNjGN+rNU/Mjuf7baffS/HM3ATeXbaXW8d0ZMrIOA9E6Fm+bDG0AeaLyHpgBc45hi9F5DYRuc1d51Jgo4isA54BrlS7tMeYOu+T1UkUFClXDfXPoC8R4V+X96Nb64bc+e5qNu7LqNZ2VJXp8xP4z3fbuGRgNPeP7+7hSD3DZ5erqup6YEAZz79Y4v5zwHO+ismbjh07xujRo8nLy6OwsJBLL72URx89aeiGMeYUVJUZKxIZ3KEJnVs28Fsc9cNDePWGwVz2whKufHkpL147iJFdKl/ALz0nnz99soHZGw8wqV9b/j6572lNLuRNNvLZS8LDw5k3bx7r1q1j7dq1fP311yxdWnrIhjHmVFbvPcLO1GyfXKJ6KtGNI/nk12fSrkkkN76+nLd/3F2p8UrfbDrAuf9ZyLebD/LAhO48c2V/r4/DOB2BG1kNJyJERUUBTs2kgoKCgP12YEwg+2JdMuEhQYzv7ftKp2Vp3SiCD24bzuiuLfjL55u48701bNyXUWaCSM7I5c7/rebWt1fRPCqcz+84k6mjOwX8Z0GtL6LH7PvhwAbPbrN1Hzj/yVOuVlRUxKBBg0hISOCOO+6wstvGVFFRsfLl+mTGdW9JAz9VGi1Lw4hQXr1+MNPnJzB9QQJfrU9mbLcWnNOjFc2jwknOyGXpzjS+i08hSODuc7ty21jvj9b2lNqfGPwoODiYtWvXkp6ezsUXX8zGjRvp3bv00A1jTHmW7UzjUFYekwLsck6AoCDhN2d34dozOvD+ykRe/H4HC7b+fMVSq4bh/OrMWK4fHkv7pv6tfVRVtT8xVOKbvbc1btyYsWPH8vXXX1tiMKYKvli/n/phwT4tfV1VTeqHcduYTkwd1ZGDR4+RkplHk3phfi+EdzpqRrumBkpNTSU9PR2A3NxcvvvuO7p3D8xL04wJRMcKipi14QDn9mxFZJjvxy5UVVCQ0KZRJP3aN67RSQHqQovBT5KTk7nhhhsoKiqiuLiYyy+/nAsuqPFFY43xmbnxKWTkFnDJQM/P0mYqZonBS/r27cuaNWv8HYYxNdbHq5No3TCCMztXfqyA8QzrSjLGBJyUo8f4flsqlwyM9uo8CKZslhiMMQHn8zX7KSpWJg+ybiR/sMRgjAkoqspHq5IYENOYTi2i/B1OnWSJwRgTUDbuy2TrwaNcaq0Fv7HEYIwJKB+vTiIsJIgL+gbeoLa6whKDMSZg5BcW8/nafZzXsxWNIgOnBEZdY4nBy4qKihgwYICNYTCmEr6LP8iRnALrRvIzSwxeNm3aNHr06OHvMIypEd5bvpfoxpGM6mJT9vqTL6f2jBCR5SKyTkQ2ichJs9aI4xkRSRCR9SIy0FfxeUNSUhJfffUVN998s79DMSbg7U3LYdH2Q1w+uL2NXfAzX458zgPGqWqWiIQCP4jIbFUtOXvN+UAX9zYMeMH9WW1/X/53thzecuoVq6B70+7cN/S+U673+9//nqeeeoqjR496dP/G1EYzVuwlSODyIdaN5G8+azGoI8t9GOreSs9scSHwlrvuUqCxiLTxVYye9OWXX9KyZUsGDRrk71CMCXgFRcV8uCqJs7q1pE2jSH+HU+f5tFaSiAQDq4DOwHRVXVZqlWggscTjJPe55FLbmQpMBYiJialwn5X5Zu8NixcvZubMmcyaNYtjx46RmZnJtddeyzvvvOOXeIwJZHPjD5J6NI+rhlb8/2x8w6cnn1W1SFX7A+2AoSJSenKCsjoWT5ovT1VfVtXBqjq4RYvAPEn1xBNPkJSUxO7du5kxYwbjxo2zpGBMOd5dtpc2jSIY2y0w/5/rGr9claSq6cACYHypRUlAyRm/2wH7fROVMcYfjp90vmJIe0JqyNSXtZ0vr0pqISKN3fuRwDlA6bPCM4Hr3auTzgAyVDWZGm7s2LF8+eWX/g7DmIB0/KTzFUPan3pl4xO+PMfQBnjTPc8QBHygql+KyG0AqvoiMAuYACQAOcBNPozPGONjBUXFfLAyiXHd7aRzIPFZYlDV9cCAMp5/scR9Be7wVUzGGP+as/kgh7LyuHqYnXQOJNahZ4zxm/8t20vbRhGM6drS36GYEiwxGGP8Yk9aNj8kHOKKITE20jnAWGIwxvjFe8sTCQ4SO+kcgCwxGGN8Lr+wmI9WJTKue0taN4rwdzimFJ+OfK5rYmNjadCgAcHBwYSEhLBy5Up/h2RMQJi1IZlDWflce0YHf4diynDKxCAilb1cIF1VM08znlpn/vz5NG/e3N9hGBNQXl+ym47N6zOqs/1vBKLKtBjexClLUdHZIQXeAN7yQEzGmFpsbWI66xLTeWRST4LspHNAOmViUNWzSj8nIq1V9YB3QvKsA3/7G3nxni27Hd6jO60feOCU64kI5513HiLCrbfeytSpUz0ahzE10ZtLdhMVHsJkm6UtYFX3HMP1wFOeDKQ2Wrx4MW3btiUlJYVzzz2X7t27M3r0aH+HZYzfpB7N48v1+7lmWAcaRNiczoGquonhQhHJAeao6lZPBuRplflm7y1t27YFoGXLllx88cUsX77cEoOp095bvpeCIuX64XbSOZBV93LVS3DqGV0sIq96MJ5aIzs7+6eZ27Kzs/n222/p3bt0lXFj6o78wmLeWbqHMV1b0LFFlL/DMRWoVotBVQ8CX7s3U4aDBw9y8cUXA1BYWMjVV1/N+PGlq4wbU3d8vekAKUfz+PvkWH+HYk6hWolBRKYD9VX1RhE5T1W/9XBcNV7Hjh1Zt26dv8MwJmC8uWQ3sc3qMaarTcYT6KrblZQP7HTvj/NQLMaYWmpdYjqr9hzhuuGxdolqDVDdxJADNBKRUMDq5RpjKvTc/AQaRYZaXaQaorqJ4TCwA5gOLPZcOMaY2mbLgUzmbD7IjSNiiQq3Kjw1QZUSg4g0FpHXgcnuU28Bgyv52vYiMl9E4kVkk4j8rox1xopIhoisdW8PVSU+Y0zgeWHBDuqHBXPTmbH+DsVUUpXSt6qmi8iTQCxwCOgLfFLJlxcCd6vqahFpAKwSkTmqurnUeotU9YKqxGWMCUy7D2Xzxbr93DKqI43rhfk7HFNJ1WnXTQF2qeo3wKrKvkhVk4Fk9/5REYkHooHSicEYU4Gc/EI27sukoKiYFg3C6dqqgb9DKteL3+8gJDiIKSPj/B2KqYLqJIYjwG0i0g1YB6xV1TVV2YCIxOLM/7ysjMXDRWQdsB+4R1U3lfH6qcBUgJiYwD33nZ6ezs0338zGjRsREV577TWGDx/u77BMDVVUrLz4/Q6mz08gJ7/op+cHxjTmoUm96N++sf+CK8P+9Fw+Xp3ElUNiaNnQ5lyoSaqcGFT1CRGZC2wD+gOjgUonBhGJAj4Gfl9Gme7VQAdVzRKRCcBnQJcyYngZeBlg8ODBWtX34Cu/+93vGD9+PB999BH5+fnk5OT4OyRTQxUWFfPbGWuYteEA43u15vIh7agfFkJ8ciYvfL+DyS8s4clL+nDZ4MC56ueFBTtQhVvHdPR3KKaKqpwYROQxIBhYi9NaWFCF14biJIV3VfWkcxMlE4WqzhKR50Wkuaoeqmqc/paZmcnChQt54403AAgLCyMszPpYTdWpKg9+upFZGw7wwITu3DKqIyLOWIBhHZtx8cB23PHuau77eD0NIkIZ37u1nyN2zi28t3wvVwxpT7sm9fwdjqmi6rQYHhKRVjhdQZNFpJOq3nKq14nzl/xfIF5V/13OOq2Bg6qqIjIU56qptKrGWNKiD7ZxKDHrdDZxkubtoxh1edcK19m5cyctWrTgpptuYt26dQwaNIhp06ZRv359j8Ziar+/f72V91cm8ttxnZk6utNJyxtFhvLy9YO45tVl/Pa9NXxw23C/dyv9e842QoOD+N3ZJzX4TQ1Q3XEMtwJrVPXJyiQF15nAdcC4EpejThCR20TkNnedS4GN7jmGZ4ArVTVgu4oqUlhYyOrVq7n99ttZs2YN9evX58knn/R3WKaG+WzNPl78fgfXDIvhrnPL/zJSLyyE124YQosG4dzx7mqOZOf7MMoTbdyXwcx1+/nVyFg7t1BDVXe0yWvA7SJSH6dbaO2pXqCqP1DxLHCo6nPAc9WMqUyn+mbvLe3ataNdu3YMGzYMgEsvvdQSg6mSpCM5/OWzjQzu0ITHLuz9U/dReZrUD+OFawdy6Qs/cs+H63j1hsGnfI03PPXNVhrXC+XWMSe3bkzNUN0Ww29xkkoIzjd7U0rr1q1p3749W7c601XMnTuXnj17+jkqU1MUFSt3f7AOBf5zRX+CK1lfqG+7xtx3fnfmbknh87X7vRtkGZbsOMTCbancMbYzDW0inhqruolhBxABfK6qNvNMOZ599lmuueYa+vbty9q1a3nAj5MGmZrllUU7WbbrMA9P6kn7plU7eXvjiFj6t2/Mo19sIi0rz0sRnkxV+fvXW2nbKILrbCKeGq26iWETMA+YIiIrPBhPrdK/f39WrlzJ+vXr+eyzz2jSpIm/QzI1wOb9mfzr262c37s1l1ZjXuTgIOGpS/uSlVfIo1/4bvzo1xsPsC4xnd+f25WI0GCf7dd4XnUTQ1ecS1ZfBm7yXDjG1G2qyp8/20CjyDD+dnGfap8j6NqqAXee1YWZ6/Yzf2uKh6M8WU5+IY/PiqdLyygmD6x6MjOBpbqJoTvOoLZ7cEcgG2NO3+yNB1i9N517zutKk/qnN+7l9rGd6Ni8Po99sZn8wmIPRVi2f3+7jaQjuTx+cZ9Knw8xgau6iaExcB9wL3DMY9EYU4flFxbz96+30K1VA4+MYA4LCeIvk3qy61A2ry/e5YEIy7Y+KZ3XFu/immExDI1r6rX9GN+pbmJ4DOfE81bAu19FjKkj3l22hz1pOdw/obvHvnWf1a0lZ3dvybS529l1KNsj2yypoKiY+z7eQIsG4dx3fnePb9/4R6USg4gEi0iyiNwMoKpJqvqde/9+bwZoTF2QkVvAtLnbObNzM8Z6eE7kv17U2xmFPGONx7uU/vntVuKTM3nswt52eWotUqnEoKpFwEbARqwY4wXPL0ggI7eAByb08PigtLaNI/n75D6sT8rgX3O2emy787em8NL3O7l6WAy/6OX/+kzGc6rSlVQPuFdEVorITPf2ubcCq+m2bt1K//79f7o1bNiQp59+2t9hmQCUdCSH1xfv5uIB0fRq28gr+xjfuw1XD4vhpe938sP2069JeSDjGHd/sI7urRvw0AU2cLO2qUpJjOMTCQx0bwA1so6RL3Tr1o21a9cCUFRURHR0NBdffLF/gzIB6Z/fbEWAe87r5tX9/GViT1bsOsxdH6xl9u9G0TwqvFrbyc0v4vZ3V3GsoIjnrh5oYxZqoaq0GOLKuFmh9UqYO3cunTp1okMHGw1qTrQhKYPP1u5nysg42jaO9Oq+IsOCefbqAWTmFnD3B+soKq7697qMnAJufH05axPT+ffl/ejcMsoLkRp/q3SLQVX3eDMQb5n/xsuk7Nnp0W227NCRs26s/PCNGTNmcNVVV3k0BlM7PP3dNprUC+W2sb45fde9dUMemtSTBz/dyO3vrOKZqwZU+hv/usR0/vDBWhIP5/L0Ff0Z37uNl6M1/lLdy1VNJeXn5zNz5kwuu+wyf4diAsy2g0eZuyWFG0fE+fSKnmuGdeDRX/ZiTvxBrn11Gek5FZfoLipWnv5uGxc/v5jMY4W8NWUoF/aP9lG0xh+qW3a7xqjKN3tvmD17NgMHDqRVq1Z+jcMEnpcX7iQyNJjr/VBw7oYRsTSPCueu99cy/ulF/N9FvTmn58l/o3vTcvjjR+tYtuswFw+I5rELe9HALkut9aozteckVf2iGq9rD7wFtMYZFPeyqk4rtY4A04AJQA5wo6quruq+Asl7771n3UjmJMkZuXy+dh/XDOtw2qUvqmti3za0axLJvR+t5+a3VtKvfWPO6taCPtGNOJiZx6LtqXyz6QDhIcH887J+1SroZ2qm6rQYHgeqnBiAQuBuVV0tIg2AVSIyR1VLln88H+ji3oYBL7g/a6ScnBzmzJnDSy+95O9QTIB5ffFuihWmjIzzaxz92jfmi9+M5L3le/lwVSLT5m7n+JyJTeuHccvojvzqzDha2UxsdUp1EkO1Rt+oajKQ7N4/KiLxQDRQMjFcCLzlTue5VEQai0gb97U1Tr169UhLO60pq00tlHmsgP8t28vEPm2qPNeCN4SFBHHDiFhuGBFLRk4BOw9l0aReGB2a1fPLDHDG/6qTGE577IKIxAIDgGWlFkUDiSUeJ7nPnZAYRGQqblXXmJiY0w3HGJ/637K9ZOUVMnV04F3t3aheKANibN6Qus7nVyWJSBTwMfB7Vc0svbiMl5yUiFT1ZVUdrKqDW7TwbF0ZY7wpr7CI137YxcjOzekd7Z1RzsacLp8mBhEJxUkK76rqJ2WskgSUrDfcDvD9xLXGeMnna/aTcjSPW8cEXmvBmOOqkxgOVmdH7hVH/wXiVfXf5aw2E7heHGcAGTX1/IIxpRUXKy8t3EHPNg0Z2bm5v8MxplxVPsegqudWc19nAtcBG0RkrfvcA0CMu90XgVk4l6om4FyuatOGmlpj7pYUdqRmM+3K/nZS1wQ0nw1wU9UfOMUVTe7VSHf4JiJjfOul73cQ3TiSiX2slIQJbFYSw4v+85//0KtXL3r37s1VV13FsWM2C2pdtWrPYVbuOcIto+IICbZ/OxPYqvUXKiJ/KHHfu7WCa6h9+/bxzDPPsHLlSjZu3EhRUREzZszwd1jGT176fieN64Vy+ZDTn8vZGG+rUleSiDQG/gN0F5FjwHpgCnYuoEyFhYXk5uYSGhpKTk4Obdu29XdIxg8SUrKYE3+Q34zrQr2wWl+ezNQCVforVdV04CYRmQgcAM4DyrrsNGCkf7GD/P2enQQ9rG19Gk+quExydHQ099xzDzExMURGRnLeeedx3nnneTQOUzO8umgnYcFB3OCHYnnGVEd1OzvH4Fy2egZQ3auUarUjR47w+eefs2vXLvbv3092djbvvPOOv8MyPnYg4xifrN7HZYPb0ayaM6YZ42vVbdc2Bu4D7sXpSgpYp/pm7y3fffcdcXFxHB+Zfckll7BkyRKuvfZav8Rj/OP5BQkUq3LraP/8HRpTHdVtMTwGfK6qW3FKaJtSYmJiWLp0KTk5Oagqc+fOpUePHv4Oy/jQ/vRcZixP5LLB7QKiWJ4xlVWtFoOqJuGUr0BV7/doRLXEsGHDuPTSSxk4cCAhISEMGDCAqVP9O2mQ8a0XFuxAUe44q7O/QzGmSqqVGERkOlBfVW8UkfNU9VsPx1UrPProozz66KP+DsP4wf70XN5fkchlg9vTrom1FkzNUt2upHxgp3t/nIdiMabWmD4/wVoLpsaqbmLIARq51VJtQgRjSkg6ksMHKxO5fHB7ohtH+jscY6qsulclHQZygenAYs+FY0zNN33+DgSx1oKpsarUYnCn2nwdmOw+9RYw2ONRGVNDJR3J4cOViVwxpD1trbVgaqgqj3wWkSeBWOAQ0JcAH/lsjC9Nn59AkAi/PsvGLZiaqzpdSVOAXar6DbDKw/EYU2MlHs7hw5VJXD0shjaNrLVgaq7qnHw+AtwmIk+LyE0iMsDTQdUW06ZNo3fv3vTq1Yunn37a3+EYL3tuXgJBQcKvx9q5BVOzVTkxqOoTwC3AI8AuYHRlXicir4lIiohsLGf5WBHJEJG17u2hqsYWSDZu3Mgrr7zC8uXLWbduHV9++SXbt2/3d1jGS/am5fDR6iSuHhpD60YR/g7HmNNS5cQgIo8BF+IUz9unqtMq+dI3gPGnWGeRqvZ3b49VNbZAEh8fzxlnnEG9evUICQlhzJgxfPrpp/4Oy3jJM/O2Exwk3D7Wzi2Ymq86cz4/5H6bDwImi0gnVb2lEq9bKCKx1YjxtMyePZsDBw54dJutW7fm/PPPr3Cd3r178+CDD5KWlkZkZCSzZs1i8GC7gKs22rQ/g49XJzHlzDhaNbTWgqn5qjuO4TXgZqA+8LznwmG4iKwD9gP3qOqmslYSkanAVHCK1QWiHj16cN9993HuuecSFRVFv379CAmxSVpqG1Xl8a/iaRwZym/GdfF3OMZ4RHU/qX6LUxYjBJhGJc8znMJqoIOqZonIBOAzoMz/NFV9GXgZYPDgwVrRRk/1zd6bpkyZwpQpTlXyBx54gHbt2vktFuMdc+NTWLIjjUd/2YtG9UL9HY4xHlHdkhg7gAic0tueSAqoaqaqZrn3ZwGhItLcE9v2l5SUFAD27t3LJ598wlVXXeXniIwn5RcW87dZ8XRqUZ+rhwVmy9WY6qhui2ETkAhMEZF/qOqQ0w1ERFoDB1VVRWQoTtJKO93t+tPkyZNJS0sjNDSU6dOn06RJE3+HZDzo3WV72Hkom9dvHEJocHW/YxkTeKqbGDrhjGd42f15SiLyHjAWaC4iScDDQCiAqr4IXArcLiKFOHWYrlTVCruJAt2iRYv8HYLxkiPZ+Tz93XZGdWnO2G4t/B2OMR5V3cSQqKrzRKQNkFKZF6hqhf0oqvoc8Fw14zHGp576ZitZeYX85YKeiIi/wzHGo6rb/h0vIu2AF4H/eDAeYwLe2sR0ZqzYy6/OjKVrqwb+DscYj6tuYmgM3AfcC+R5LBoPqgm9UDUhRnOiwqJiHvx0Ay2iwvndOV39HY4xXlHdxPAYzhVJW4EiD8bjEREREaSlpQX0B6+qkpaWRkSEDYiqSf77wy427c/kkV/2IircxqWY2qlSf9kiEgwkAX9R1VdVNcl9jKre78X4qqVdu3YkJSWRmprq71AqFBERYWMbapDdh7L595xtnNezFef3bu3vcIzxmkolBlUtcovf1YhCMKGhocTFxfk7DFOLqCp/+mQDYcFB/PWi3nbC2dRqVWkL1wPuFZFzcUpWAKiqXuj5sIwJLB+sTOTHnWn87eI+Vg/J1HpVSQzD3Z8D3RtA4HbiG+Mhe9KyefyreIbFNeXKIe39HY4xXleVxGB9M6bOyTxWwJQ3VxIUJPzj0n4EBVkXkgkcquqVbs1TJgYROV4EpszWQYnl6aqa6anAjPE3VeWuGWvZfSibt6cMI6ZZPX+HZOq43MJcftj3A3N2z2H9ofVM7jKZW/qectaDKqtMi+FNnKRQUVpSnIl43vJATMYEhLd+3MPcLSk8PKknwzs183c4pg7KK8pjyb4lzE+cz4/JP5KWm0ZBcQFNI5oytPVQ4hp5pyPnlIlBVc/yyp6NCWAbkjJ4fFY8Z3VrwY0jYv0djqlj8ovy+WT7J7yy/hVSclOICo1iRNsRREdFc2b0mQxqNYiQIO+No7EROsaUsnFfBle9spTm9cN46tJ+dmmq8ZkD2Qd4fePrfLfnO1JyUxjYciAPj3iY4W2GExrsu/k+LDEYU8KO1CxufnMlDSNC+Oj2EbRoEO7vkEwdsDR5KZ9s/4S5e+aiKCOjR/LXbn9leNvhfvliYonBGNfetByueOlHVOGdm4fRtnGkv0MytZSqsj97P59s/4RVB1ex6uAqmkY0ZVKnSdzS9xaio6L9Gp8lBmOArQeO8qs3VlBYrHx02wg6t4zyd0imFiooKmB+4nxe2fAKWw5vIViCiW0Yyz2D7+HK7lcSHhwYLVRLDCZg5RcWsz4pnc3JmWTkFNAgIoQ+7RrRv30Tgj04nmDVnsPc+NoKIsOCeWfKMEsKxuPi0+L5ePvHzNo1i6P5R4ltGMt9Q+5jdLvRxDQMvGlhfZYYROQ14AIgRVV7l7FcgGnABCAHuFFVV/sqPhM4MnIK+O8PO3l98W6O5hWetDy2WT2mju7EJQOjiQgNPq197UjN4ta3V9EsKoz/3XKGdR8Zj1BVEtITWHlwJZ9u/5T4w/GEB4dzTodzGB87nlHRowgOOr2/XW/yZYvhDZwZ2sob63A+0MW9DQNecH+aOuR/y/byxOx4jh4r5PzerbmwfzQDYhrTtH4YaVn5LNuVxquLdvHApxt4dt52/jShB5P6tqnyCbriYmXGikSemB1PWHAQr94w2JKCOW2Hcg/x0baP+GrnV+zO3A1A1yZd+dPQPzGx40QahTfyb4CVJL6cs0BEYoEvy2kxvAQsUNX33MdbgbGqmlzRNgcPHqwrV670RrjGh4qLlSe/3sLLC3dyZudmPDihJz3bNixzXVVlcUIaT8yOZ9P+TIbGNeWxC3vRvXXZ65eWlpXHb2esYXFCGmd0bMo/Lu1H+6Y2qtlUX05BDov3L+bJ5U+SkpPC4FaDmdBxAoNaDiKuUVxAXvIsIqtUdXBZywLpHEM0kFjicZL73EmJQUSmAlMBYmICr3/OVM2xgiLuen8tszce4PrhHXh4Uq8KzyGICCO7NGdmp5F8sDKRp77ewsRnfuCG4bHcdGZsuR/yqUfzePvH3by8aCfFxfDEJX24ckj7gPynNYGvqLiILUe2MG/vPN7e/Da5hbm0qteKDyd9SPem3f0d3mkJpMRQ1n9nmc0ZVX0ZeBmcFoM3gzLelZaVx81vrWRtYjp/ntiDKSMr/+0qOEi4amgM43u15qlvtvD6kl28vmQXIzs3P6G1kZlbwNYDR1mbmE6xwqR+bfnNuM42X7Oplp0ZO5m+ZjpLk5eSme+Uhzs75myu7n41A1oNIDTIdwPRvCWQEkMSULKmcTt+nvfB1EKpR/O44uUf2XcklxeuGcj43m2qtZ0m9cN44pK+/GZcFz5YmchHq5JYvuvwT8vrh4fQqUV9fj22M5P6taVba0sIpuqKiot4e/PbvLDuBUKCQhgXM44z2pzBsDbDaB7Z3N/heVQgJYaZwJ0iMgPnpHPGqc4vmJorLSuPa15dSnL6Md6eMoyhcU1Pe5ttG0fy+3O68vtzunogQmMc2QXZvLXpLWbvns2ujF2Mbjeav5zxF1rXr73Tu/ryctX3gLFAcxFJAh4GQgFU9UVgFs6lqgk4l6ve5KvYjG+l5+Rz7X+Xsycth9dvGuKRpGCMp+3O2M33Sd/z2sbXOHzsMENbD+XmPjczqeOkWn9eymeJQVWvOsVyBe7wUTjGTzKPFXD9a8vZkZLFqzcMZkSn2tUEN7XDgsQF3DX/Lgq1kCGth/DcwOfo06KPv8PymUDqSjK1XFZeITe8tpz45Exeum4Qo7u28HdIprJUIS8Tsg9BzmHIOQQ5ac7zYfUhLMr9WR8im0D9FhBW8y4BVlW+2/sd9y28j+5Nu/P4qMeJaxiYl5t6kyUG4xPHCoqY8sYKNiRlMP2agYzr3srfIRmAowchNAIiGkFeFqCgxVBwDA5uhN0/QOoWSFoJ2SlV23ZYFDRsCw2joW1/aNkLtAiadoR2QyDAPmzXpa7jP6v+w6qDq+jRtAcvnvtijRmQ5mmWGIzXFRUrv5uxhmW7DjPtyv78olftPWkX8I4egF2LYM9i2LMEDm0FCYLGMXBk98nrB4VA827QYbjzYV6vOdRrBvWbQWRT57X52e4ty/mZewSyUyErBY7uh/REWPyMkxSOq9ccOoyANv2gw5nOtoP983G0O2M3c/bM4dk1z9I0oikPDnuQyV0m+3T+g0BjicF4lary8MyNfLPpIA9d0JML+/u3nHCdlJcF696D7d9CwndOiyC8IcScAf2vhoIcSF4P/a6G4x+GYVHQrJOzTlh9D8RwFDL3O4kmcbnTEkn4DuJnOssjm0DvydBxLES1guZdIbLx6e+3AtkF2UxfO51349+lWIsZFT2Kf475J/VCa14XmKdZYjBe9dLCnbyzdC+3ju7Ir0Z6Z35aU46sVJj3GGz+HI5lQNNOMPwO6H0ptO4DviziFt4AWnRz7jfrBP3da1GOZcCOebBlFqx6A1a8+vNrmnaE/tfAoBuhvmcuUsgvyichPYHVB1fzxqY3OJhzkMu6XsYV3a6gS5MuBEmQR/bjE3lZTtL2QpecT2sleYPVSgpcry/exaNfbGZSv7ZMu6I/QR4slW0qkHPYSQbfP+WcIO55IQy9BdoP9XdkFcvcD1kHne6u1C2QMBd2L3LOf8SNgVF/gLYDqrzZ/KJ8Vh5cyeJ9i/ks4bOfRiv3bd6XPw75I/1b9vfwG/GBjCR4+2Lod5VzXKqhptRKMrXI8wsSeOrrrfyiVyv+eVlfSwq+smUWzLzTSQit+8JV7zknfmuChm2dG0C382HkXZASDwuedLqets6Gwb+Cs//itEAqcDT/KN/u/pZZu2axOW0zWQVZhEgIZ3c4m3M7nEuPpj0Cch6ESknd5iSFvExo750C1JYYjEepKv+Zs41n5iVwYf+2/OuyfoQE16DmeU2Vcxhm3wcbPoBWfeDqDyF6YMBd+VNlLXvA5W9CdhrM+yuseAW2zYa+VzjJI3rQT6smHElg4b6FbEnbwrzEeeQV5RHXKI7z485nbPux9G3el8YRjf33Xjxh32p4Z7LTDXjjV9Cmr1d2Y11JxmNUlSdmO6Wzrxjcnr9d0sejM62ZMhQVwtp3YN7/OVcDjboHRt0NIWH+jsw79i5zWkSHtgFQHNWaj5q14r2wIhIKnS6iphFNObfDuVzU+SJ6NetVe8Yg7PweZlwN9ZrCdZ8552pOg3UlGa/LyC3gN++tYeG2VG4cEctDF/S07iNvOnoA1rwDGz+GlM3Q/gyY8A+vfYMMGDHDKLh9MauSfuC7NS8xLzOBVD1C36w8HsjK5txG3Wg+6iHn8ldPXE0VCFRh1etOi7BpJ7juU2hYvYKTlWWJwZy29Jx8rvvvcrYcyOSvF/Xm2mExtedbWqBJiYclz8KWr+BYunNZ5+VvQY9f1vxuowrkFeXx4/4f+W7PdyxIWkBGXgaRIZGMjBnDhLgJnB0Vh+yY57Sc3roQQus5l+IOux2ad/Z3+NWXnw1f3gXr34dOZ8PkV50Wg5dZYjCnxamSuoydh7J5+brBnNW9pb9Dqp0OboZZf4Q9PzhjDLqcC2P/9PMloLVQTkEOC/ctZO6euSxMWkhOYQ4NQhswtv1Yzu5wNiPajiAypMR0rM06OZe37lzgJM7VbzmXv7YdCINucK7gCQn32/upskPb4f3rnCu0znrQ6SYM8s35OjvHYKrtcHY+V728lN1p2bxy/WCrfeRpqrBzPix5zvmwi2jkXJrY7yqPXdcfaIqKi1h2YBlf7PiCuXvnkluYS9OIpoyLGcc5MecwtPXQyo9IzkqB1W86l+4e2ABRreGM250kEdnEu2/kdG2eCZ/92jlXNPlV6DTO47uo6ByDJQZTLUey87n61WXsTM3itRuHcGbn2vlB5Re5R5xvu6vfgrQEp9ZQn0vhzN/7pBvBH7Yf2c4XO77gq51fkZKbQoOwBvwi9hdMjJvIgJYDCD6dwXiqsOt7+OE/ToINCnVGWZ9xW7XGRXhVXhZ8+6Az2C96kNNN2KidV3ZlicF4VMnuo1etpeA5Rw86l2Ou+C/kHnZOKA+6AXpd4hS6q2USMxP5cteXfLv7WxLSEwiREEZGj2RSp0mMaT+G8GAvdPskr4M17zpJtzDXGesx8Hroe7nTIvOnvUvh01vhyB4Y8RsY92evdn0FTGIQkfHANCAYeFVVnyy1fCzwObDLfeoTVX2som1aYvCt1KPOzGt7D+fw6vVDGNnFWgqn7fBOWPqi0+1RlA+dz4VxDzoF5mqRguIC1hxcw6a0TSxIXMDqlNUIwoCWAzgv9jzOjzufphE+ahHlpsOGD51jfmADhERCr4udlkTcKN+ei8g5DAuegOWvOMUML37RKTDoZQGRGEQkGNgGnIszv/MK4CpV3VxinbHAPap6QWW3a4nBd9YnpXPH/1Zz6Gg+r904hOGdmvk7pJpL1Sn3sPQFZ0RvUIjzrXXU3ad9fXogKdZi1qasZfau2Xy751sOH3Pm4u7UqBOTOk1iYseJ/p0iUxX2r3ESxIaPnAqxYQ2gyznQ/QLnJL83WxK7FsEH1zk1owZPgXMePuWobk8JlHEMQ4EEVd3pBjUDuBDYXOGrTED4YGUif/5sIy2iwvnfLcMYEBPgJ+8CVXGRczJ00b/h4AanhPXoP8KQKdCgdpQjV1U2HtrI17u/5pvd33Aw5yDhweGMbjeaiR0nMrDlQJpEBMjfj4gzQjx6IIz/u3MuYstXsHUWbPrUOR8RNwq6T4RuE5wrwvKznJbd8XLl1bHze1j0L2d/zbvBjbOgVU/PvrfT4MsWw6XAeFW92X18HTBMVe8ssc5Y4GOcFsV+nNbDpjK2NRWYChATEzNoz549Xo+/riooKuavX27mrR/3cGbnZjx71UCa1q+lo2q9qajA6bpY9G9I2+6MPxjxG+hzGYRGnvr1NUB+UT5vbnqTj7d/zL6sfYQEOecMxseOZ2z7sdQPrUEDzoqLnMmJtn4F8V/C4R1lr9d2ADTr4nT9tOzpJPcGbcoeeX5kj1PYMHGZ8zfQMNqp/TTkZq+XGC9LoHQlXQb8olRiGKqqvymxTkOgWFWzRGQCME1Vu1S0XetK8p4DGcf4zXurWbH7CLeMiuO+8d2t7lFVFebB2nedK2LS9zp1jEbfAz0m+bbstRcUFReRXZjNiuQVzEucx8KkhaTnpTO8zXAmdJzAuJhxNAxr6O8wT5+qU4Jj62xnsqHIps68FZnJsHshHEpwJiQ6ToKcD/3iQueKovAGTjXUtB0QEuHMORF7ptN15MeLCgKlKykJaF/icTucVsFPVDWzxP1ZIvK8iDRX1UM+itG45sYf5I8fredYQRHTruxvE+xUVW66kwzWvO1UOo0eBOc/BV3H1/gRykfzjzJt9TS+3Pkl2QXZADQMa8jodqP5ZadfMrztcD9H6GEizkDCsgYTjvnjz4kjPRGOJkNGIhze5SSP9L3OvBgtujvnLAb/Chq3P3k7AcaXiWEF0EVE4oB9wJXA1SVXEJHWwEFVVREZCgQBaT6Msc7LyC3g/77czIerkujeugHPXT2Qzi2j/B1WzZF7xCnLsOoN5xtjzwudiWY6nlXjE0KxFjNzx0yeXvU0R/KOMKnjJDo27kivZr0Y2GogoUF1dCrMihJHDeWzxKCqhSJyJ/ANzuWqr6nqJhG5zV3+InApcLuIFAK5wJVa0wda1CDfb0vl/o/XczDzGHec1Ynfnt2F8JCa3d3hM5n7nZPKC//pjEEYdCMMvKHmzIVwCpvTNvP4ssdZn7qevi36Mv2c6fRq1svfYRkvsQFuhqPHCvjbrHjeW55I55ZR/OuyfvRr39jfYdUMeUdh2YvOSeWCHGfilAn/rPFVTvOK8pifOJ96IfVYlLSID7Z9QOPwxtw16C5+2emXNWsKTFOmQDnHYALQ4oRD3PvRepIzcrl1TEfuOqcrEaHWSqhQXpZT8nr1W5C+x7l8sdtEOPshpzuhBncZ7UzfyVMrnmJT2ibS89IBCJIgruh2BXcOuLN2nEw2p2SJoY7KzivkidnxvLN0Lx2b1+fD20YwqEOAXFseaPKznUSQl+mUc178DGQdgHZDnZnE+l8D7QadejsBIr8on9c3vs73Sd/TKLwRPZv1pHvT7mw9vJX/bvwv9UPrM6bdGM6JOYdmkc1oFtmM6Ci7+KAuscRQx6gq329L5S+fbyTpSC43j4zjnl90s1ZCWYoKYOMnMPcxyEz6+fnoQc50kzFn+C+2atibuZdNaZuYvnY6ezL3MKDlAFJyUvhx/48UaREAE+ImcO+Qe2kWaaPa6zJLDHXI4oRDPP5VPJuTM+nQrB7vTx3O0LjaWa2z2nIOO5OibJ3lTLqedQBa9XZKH7fpC1kHoUlcjekuKtZiNhzawP/i/8fXu7+mWIvp0LADL537EiPaOvV4cgtz2ZWxi4iQCDo26ujniE0gsMRQBxwrKOLJ2Vt4Y8luYpvV46nJffll/7bWSjjuWKZzvfnS553RyUX5TjJoPwQGXOcUtTs+QUrTwPvgVFUKiwv5YNsHfLnjS4q0iNb1W9Mkoglf7fyKvKI8GoY15KruVzE+djw9m/UkLPjnkbmRIZH0bBY45RiM/1liqMUycgr4ZvMBXliwg12HsrlxRCz3n9+9bieE4mJnjuS9Pzqjkjd+5BRRAwgOdy4xHXQjtO7t1zBLSs1JRUSICI5ARMgvyic+LZ5NaZv4YNsHZOZlEhocSkZeBn2a96FlZEu2HdnG4WOHOTvmbPq16MeFnS+sWSUpjF9ZYqiFioqVlxbu4Nm5CeQWFNG5ZRTvTBlWN0tkF+bBwY1OvZtD22DbN1Bc8PPy5t3grD87k6t3OQ+i/D81qapSqIWsOriKVze8yrLkZeWuO7zNcOIaxZFbmMs5Hc5hVPQom2/bnDZLDLXMou2pTPtuOyv3HOEXvVrx67Gd6RPdiKCgOvRhUVwM8TOdy0mT10HOIZBgpxTBwOugRQ/oPsFpIdRr5rN5dEsrKC5g25Ft5Bbkkl+Uz7GiYxzIPsBH2z9i+5HtADSPbM6d/e+kUXgj8oryfnptt6bdaFu/LTENq1nd05gKWGKoJbLyCpk+P4EXFuygeVQ4/7ysH5MHRtetb48HNsCPz8Ouhc5VRE3iIG409PylM/CsYVu/hKWq7Mvax5L9S8gtzKWguID0Y+nM3jWblNyUk9bv1qQbt/S5hXYN2jGx40TvzGRmTAUsMdRwOfmFvPXjHl76fgdHcgq4dFA7/nZxH8JC6sjI1OJi2DkPljwHO+dDaH1ncpXuDzuzcXmxgumezD2sPLCSo/lHOZJ35KTlWw9vJf5wPJl5mRRq4QnLBGFE2xHcPfhumkY2JSI4grDgMOqH1iemQUzdSugm4FhiqKFUlS/XJ/PYl5tJPZrH6K4t+MO5XelfF0pZFBVAWoJTBnn1m3BkN0S1hnMecU4cR1Z/oJ6qsu3INlYeXElh8c8f5nlFeWxO20xBcQGpOakkZyf/NDIYIDQoFOHED/Pmkc0ZFzOOxuGNaRrRlNHtRtM0oilhwWGEBoVaWQkTsCwx1DDFxcrsjQd4dt52thw4Su/ohrxwzUAGx9bC8QhFhZCyySlbnRIPBzY6J5JTtziXlAJ0GAlnPehUMa3EPL2Hjx1mQeICjhUeo6C4gMLiQnIKc0jMTCT+cDwHcw6SW5hb5mtjGsQQFRZFo/BG9Gneh46NOzIyeiRNIppYqQhTq1gRvRpCVflifTLPzt3O9pQsOraoz51ndeaX/drWmslzcgpyWBH/IXm7FvBRRjz7C446E6McJ8HOh39wmPMzJMKpeV8FablpZBVknfCcILSNakuXJl1oF9WOTo07MTJ6JA3CGpywTr3Qeqf1/owJJFZEr4Y7kp3PvR+vZ87mg3RpGcUzVw1gYp82BNfEK40SVzizXQWHQ0g4R4oLSNy3lIOJi/lnYTL7g5wvKo1VGB7Rwp0mMQIiGlaqRXAqEcERXNHtCtpEtSE0KPSnW3ANn03NGE+yxBCgVJVN+zOZtSGZ95bvJSuvkD9P7MGvzoyr3KWnxcXOB3BxEaCgxc5MU3r8fjFkpzolHkLCnXmHg8MgKMS55aY7M1FJUIlbif0W5TuzVIVEOBVFG7X/eTL71C1waDu07OHcD410zgFs+BA2fQLAosgIPo+qz+J6kWS5l4vGBoUyvf35tOg2ifYtehMVZhMEGeMPlhgCSEJKFt9sOsCczQfZl55L6tE8goOEMV1bcPd5XenVtpGzoqrT1565/8QNpG51TsTuW+ksy071bsDB4c4sZSW7e1wKFAAqxx8LBIeSNeou/l2cwhf7F9EirBFDotpzSdwFhDWNY1CrQXZppjEBwKeJQUTGA9NwZnB7VVWfLLVc3OUTgBzgRlVd7csYfamoWPl0zT6S03NZl5TBd/EHATgjOowHmy6gc2wxHVs2oF5YCGxfAAkAAglzYe+Ssjca0ciZcL7zOU4V0NBI59s+cuI3fxFncFdUa2ckcH6O0wrQIufDPiQSmnUq0cIo9eEvQU43T0GuO8ftTqeVARSEN+D23R+xLHXNyfElfUyIhHBbv9u4pc8tJ9TsMcYEBp8lBhEJBqYD5wJJwAoRmamqm0usdj7Qxb0NA15wfwasYwVFbEs9SHLmUdJzCsjIzedITgGZuQXUD25CVEQYxaoUFBXTtH44ocFCdk4uyRvmEXVkC3mFzgdu5+AgHmy7i9jMVciRAufD+SCwvYydRjZxJpaPHswJV0hGtYJG7Xzxtn/yQvzbvLvlXUpexFCkRWQXZHN9z+tpEuFcOioIIoIgjIweSZcmXXwap/EP/an78sSbOgt/flxUhObmosXFHt1/SNOmSGj156Iuysoma/48NL/g1CtXIKxjHPUGDDitbRyXv3s3Oauc78vhXbsS2cfzdb18dlWSiAwHHlHVX7iP/wSgqk+UWOclYIGqvuc+3gqMVdXk8rZb3auSXr3/n+SGePaP0Jg6pyZc1Hi612h46j166lqREvHULwrmV0/eXb1wKrgqyZfXOUYDiSUeJ7nPVXUdRGSqiKwUkZWpqdXsR7eRpcacPqkBt0B5j55SYpsS4p1OH1+eYyjr0JTOxZVZB1V9GXgZnBZDdYK5+YnqZVljjKntfNliSALal3jcDthfjXWMMcZ4kS8Twwqgi4jEiUgYcCUws9Q6M4HrxXEGkFHR+QVjjDGe57OuJFUtFJE7gW9wLld9TVU3icht7vIXgVk4l6om4FyuepOv4jPGGOPw6TgGVZ2F8+Ff8rkXS9xX4A5fxmSMMeZEtaP6mjHGGI+xxGCMMeYElhiMMcacwBKDMcaYE9T4iXpEJBXYU82XNwcOeTAcb7E4PacmxAgWp6fVhDh9HWMHVW1R1oIanxhOh4isLK9WSCCxOD2nJsQIFqen1YQ4AylG60oyxhhzAksMxhhjTlDXE8PL/g6gkixOz6kJMYLF6Wk1Ic6AibFOn2MwxhhzsrreYjDGGFOKJQZjjDEnqLOJQUTGi8hWEUkQkfv9HU9JIrJbRDaIyFoRWek+11RE5ojIdvdnEx/H9JqIpIjIxhLPlRuTiPzJPbZbReQXfo7zERHZ5x7PtSIywZ9xikh7EZkvIvEisklEfuc+H1DHs4I4A+14RojIchFZ58b5qPt8wBzPCmIMqGP5E1Wtczecst87gI5AGLAO6OnvuErEtxtoXuq5p4D73fv3A3/3cUyjgYHAxlPFBPR0j2k4EOce62A/xvkIcE8Z6/olTqANMNC93wDY5sYSUMezgjgD7XgKEOXeDwWWAWcE0vGsIMaAOpbHb3W1xTAUSFDVnaqaD8wALvRzTKdyIfCme/9N4CJf7lxVFwKHKxnThcAMVc1T1V0482sM9WOc5fFLnKqarKqr3ftHgXicuc0D6nhWEGd5/BWnqmqW+zDUvSkBdDwriLE8fvsfgrrblRQNJJZ4nETFf/C+psC3IrJKRKa6z7VSdzY792dLv0X3s/JiCsTje6eIrHe7mo53Kfg9ThGJBQbgfIMM2ONZKk4IsOMpIsEishZIAeaoasAdz3JihAA7llB3E4OU8VwgXbd7pqoOBM4H7hCR0f4OqIoC7fi+AHQC+gPJwL/c5/0ap4hEAR8Dv1fVzIpWLeM5f8YZcMdTVYtUtT/OPPFDRaR3Bav7Jc5yYgy4Ywl1NzEkAe1LPG4H7PdTLCdR1f3uzxTgU5wm5EERaQPg/kzxX4Q/KS+mgDq+qnrQ/acsBl7h5ya53+IUkVCcD9t3VfUT9+mAO55lxRmIx/M4VU0HFgDjCcDjWTrGQD2WdTUxrAC6iEiciIQBVwIz/RwTACJSX0QaHL8PnAdsxInvBne1G4DP/RPhCcqLaSZwpYiEi0gc0AVY7of4gJ8+FI67GOd4gp/iFBEB/gvEq+q/SywKqONZXpwBeDxbiEhj934kcA6whQA6nuXFGGjH8ie+OssdaDdgAs5VFjuAB/0dT4m4OuJcjbAO2HQ8NqAZMBfY7v5s6uO43sNp6hbgfJuZUlFMwIPusd0KnO/nON8GNgDrcf7h2vgzTmAkTrfAemCte5sQaMezgjgD7Xj2Bda48WwEHnKfD5jjWUGMAXUsj9+sJIYxxpgT1NWuJGOMMeWwxGCMMeYElhiMMcacwBKDMcaYE1hiMMYYcwJLDMaUICKNReTXJR63FZGPvLSvi0TkoXKWZbk/W4jI197YvzHlscRgzIkaAz8lBlXdr6qXemlf9wLPV7SCqqYCySJyppdiMOYklhiMOdGTQCe3Nv4/RCRW3LkdRORGEflMRL4QkV0icqeI/EFE1ojIUhFp6q7XSUS+dosgLhKR7qV3IiJdgTxVPeQ+jhORH0VkhYj8tdTqnwHXePVdG1OCJQZjTnQ/sENV+6vqH8tY3hu4GqemzeNAjqoOAH4ErnfXeRn4jaoOAu6h7FbBmcDqEo+nAS+o6hDgQKl1VwKjqvl+jKmyEH8HYEwNM1+duQmOikgG8IX7/Aagr1uJdATwoVNqCHAmWymtDZBa4vGZwGT3/tvA30ssSwHaeiZ8Y07NEoMxVZNX4n5xicfFOP9PQUC6OuWVK5ILNCr1XHn1aSLc9Y3xCetKMuZER3GmsawWdeYr2CUil4FToVRE+pWxajzQucTjxThVfuHk8wld+bnqpjFeZ4nBmBJUNQ1YLCIbReQf1dzMNcAUETleIbesaWMXAgPk5/6m3+FMyrSCk1sSZwFfVTMWY6rMqqsa4yciMg34QlW/O8V6C4ELVfWIbyIzdZ21GIzxn78B9SpaQURaAP+2pGB8yVoMxhhjTmAtBmOMMSewxGCMMeYElhiMMcacwBKDMcaYE1hiMMYYc4L/B96BbIx8+1IGAAAAAElFTkSuQmCC\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": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAElCAYAAADgCEWlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABM50lEQVR4nO3dd3gc1dX48e+RVr1asiTLlrtlLHds4QIGTI1pIbQEAkkgEAIhCXlTCckvhDQgeRNKICH0BBL8JgSIA6YbsDE2xr3JvUqW1Xvf3fv7Y0b2SlZdb9PqfJ5nn92dnZ05Gklz5pa5V4wxKKWUUu0igh2AUkqp0KKJQSmlVAeaGJRSSnWgiUEppVQHmhiUUkp1oIlBKaVUB5oYVJdE5Oci8oL9epSI1ItIZLDj6omInCkiOwO8TyMiE05yG9tEZKFvIjph293+HkUkS0SWi0idiPxeLM+KSJWIrPFHPGpg0MQQpkTkgIic32nZjSLyUX+3ZYw5ZIxJNMa4fBdh//TlBGyMWWGMOSVQMfmKMWaKMeYD6Hgi98N+Ov8ebwXKgWRjzPeABcAFQI4xZo4/YlADgyYGFRZExBHsGAag0cB2c/wu19HAAWNMQ383pMc/vGhiGMREZLiI/FtEykRkv4h8u5v1xthX7A6P7y0RkUoR2SMiX/NYN1JE7haRvXYVxToRGWl/NklE3rG/t1NEPu/xvedE5DERed3+3iciMt7+bLm92ia7KuQLIrJQRApF5EcichR4tn2ZxzZHisjL9s9XISKPdnMMmkQkzWPZqSJSLiJR9vuvikiBXcXyloiM7uY4pYjI3+z9HRSRn4pIhMfnX7O3Uyci20Vklr38gIicLyKLgLuBL9g/5yYRuUZE1nXaz/dE5NVuYhgrIh/a+3gHGNrV71FEngO+AvzQ3tfXgaeA+fb7e+3vXCoiG0WkWkQ+FpHpHts7YB//zUCDvd159nrVdvwLPdb/QER+KSIr7fjeFhHP+BZ4fPewiNxoL48Rkf8VkUMiUiIij4tInP3ZUBF5zf5OpYis8DzmykvGGH2E4QM4AJzfadmNwEf26whgHfAzIBoYB+wDPmN//nPgBfv1GMAADvv9h8CfgFhgJlAGnGd/9gNgC3AKIMAMIB1IAA4DNwEOYBZWNcYU+3vPAZXAHPvzvwOLPWI3wASP9wsBJ/AAEAPE2csK7c8jgU3Ag/a+Y4EF3RyrZcDXPN7/Dnjcfv05YA+QZ8f1U+DjruIC/gb8B0iyj9ku4Gb7s2uAIuA0+7hMAEZ3/l15Hnf7fYx9XPI8lm0ArurmZ1kF/MH+3llAXQ+/x+eAX3X192G/nwWUAnPt4/kVO9YYj7g3AiPt4z8CqAAuxvr7usB+n2Gv/wGwF5hor/8BcL/92Sg71uuAKKy/mZn2Zw8BS4A0+9j+F7jP/uw+4HH7O1HAmYAE+/9voD+CHoA+/PSLtf5p64Fqj0cjxxPDXOBQp+/8GHjWfn3sBOV5QrFPAi4gyeN79wHP2a93Apd3Ec8XgBWdlv0FuMd+/RzwlMdnFwM7PN53lRhagdhOy9oTw3yshOXow7G6BVhmvxasBHaW/f4N7JO7/T7CPo6jPePCOnG2AJM91v068IH9+i3gzh5+V10mBnvZn4Ff26+nAFXYJ+dO643CSpYJHsv+0dXv0eOY95QY/gz8stM+dgJne8T9VY/PfgQ832n9t4Cv2K8/AH7q8dk3gDc9/vZe6eJnEqABGO+xbD6w3379C6xkPKHzd/Xh/UOLXOHtc8aY1PYH1j9iu9HAcLsIXi0i1VjVGFm9bHM4UGmMqfNYdhDrahGsxLG3i++NBuZ22t/1wDCPdY56vG4EEnuJpcwY09zNZyOBg8YYZy/bAHgJqwplONZVtgFWeMT9sEfMlVgnqxGdtjEUq+R10GNZX45LX/wV+KKICPAl4J/GmJYu1hsOVJmObQQHu1ivr0YD3+v0Oxtp76fd4U7rX9Np/QVAtsc63f2Ouzs+GUA8sM5jm2/ay8Eq3e0B3haRfSJyV/9/TNWZNhgNXoexrrpy+/m9I0CaiCR5JIdRWNUk7dsdD2ztYn8fGmMu8DbgLvQ0NPBhYJSIOHpLDsaYahF5G/g8VpXRi8a+HLW382tjzN97iaUcaMNu0LWXdXVcenPCz2SMWS0irVjVJF+0H10pBoaISIJHchjV1Tb7qP1n/3Uf4z2MVWL4Wncr97KvrnpClQNNWFWORZ0/tP8Gv4eVwKYA74vIp8aY97yIQdm0xDB4rQFq7cbDOLEajaeKyGk9fckYcxj4GLhPRGLtxsibsdoEwGrA/KWI5IpluoikA68BE0XkSyISZT9OE5G8PsZbgtUO0p+frxi4X0QS7FjP6GH9fwBfBq6yX7d7HPixfdJpb2C+pvOXjdUF9J/Ar0UkSawG6u8C7V1PnwK+LyKz7eMyQbpuxC4BxnTRgPo34FHAaYzpssuxMeYgsBa4V0SiRWQBcFkPP3NvngRuE5G5dswJInKJiCR1s/4LwGUi8hn77ylWrA4BOX3Y19+B80Xk83YjdrqIzDTGuO04HhSRTAARGSEin7FfX2ofSwFqsao5g9atOlxoYhik7BPZZViNx/uxrsyeAlL68PXrsOqrjwCvYLUTvGN/9gesE+TbWP+oTwNx9pXdhcC19veOcrzhuC9+DvzVrk74fG8re/x8E4BDQCFWO0d3lgC5QIkxZpPHdl6x41wsIrVYJaGLutnGt7Dqw/cBH2ElmGfs7fwL+LW9rA54FasxtbN/2c8VIrLeY/nzwFT7uSdfxGo/qgTuwUooXjHGrAW+hpWQqrCqbG7sYf3DwOVYVZJlWKWAH9CH84wx5hBWu9L37Ng3YnVcAKvtYg+w2v4dvIvVuQGs39m7WO1pq4A/GfueEOU9OV5iVkqFKrt7ZikwyxizO9jxqPCmJQalBobbgU81KahA0MZnpUKciBzA6gn1ueBGogYLrUpSSinVgVYlKaWU6kATgxp0pIuRZ8OFdBrXSilvaGJQYck+OTaINSBckYj8QQI8n4T4YK4GpYJBE4MKZzOMMYnAeVj9+725I1epQUcTgwp7xpgdWGMfTe38mYjMEZFV9o1zxSLyqIhEe3xuROQ2Edkt1rDbj9l32bZ/3uWQ3NL1UOF9HiJaRE4XkU9FpMZ+Pt3jsx6Hr/ZYr19DdivVThODCnsiMhlrnKENXXzsAv4HaxC8+Vili290WudSrOGyZ2CNp9Q+HMPnsO7yvRJrULcVwIsAxpiz7O/OMNasaf+HdVdvob1ulv3dE7oFijU3xOvAI1jDT/8BeN0eWqTdF7GGMM/EGrzv+138bEuAsZ2GHbmB3u+eVoNcWCQGEXlGREpFpPPAbd5u7037qu61TsufE2tCm432Y6Yv9qf8Zr2IVGGN3/8U8GznFYwx64wxq40xTmPMAayhwM/utNr9xphqe9iG97GGEQFrWO37jDEF9kB9vwFmdjMGEliD7GVjDdndZqypSLvqL34JsNsY87wd14vADjqOe/SsMWaXMaYJawiSmZ03Yo/A+n9YyQB7vKcxWONWKdWtsEgMWOPKL/Lh9n6HNbxxV35gjJlpPzb6cJ/K92YZY4YYY8YbY35qD8jWgYhMtKt3jtrj8PwGj1nPbN0NFd3XIbnb9XWI6OGcOFy25xDePcXUWV+H7FbqmLBIDMaY5Vj/lMeIyHj7yn+dXZc7qR/bew9roDMV/v6MdTWea4xJxqrekZ6/csxh4Ouec14YY+KMMR93tbIxps4Y8z1jzDisq//vish5Xax6BCvpePIcwrvPjDGrsSY0ah+yW6uRVK/CIjF04wngW8aY2Vj1r3/y0XZ/LSKbReRBEenryKAqdCVhjQJbb1883N6P7/Y2JHeHocL7MUT0Uqwhyr9oD0H9BWAy3lcB9Tpkt1KewjIxiEgicDrwLxHZiFVvnG1/dqWIbO3i8VYfNv1jYBJWQ2Qa1nDAamD7PtaVdB3WuP//19cv9mFI7p/TcajwPg0RbYypwGrw/h7WnMk/BC41xpT394ez9XXIbqWAMBorSUTGAK8ZY6aKSDKw0xiT3cvXetreQuD7xphLvflcqVChQ3ar/grLEoMxphbY316sF8uMXr7WKxFpL3W0j3Tpk15QSvmZDtmt+iUsSgwi8iKwEKs3SQnWzFXLsBoWs4EoYLEx5hd93N4KrCqjRKyi/M3GmLdEZBlWH3TBmmHqNmNMvU9/GKV8SDyG7DbGdHUfh1InCIvEoJRSynfCsipJKaWU9wb80LxDhw41Y8aMCXYYSik1oKxbt67cGJPR1WcDPjGMGTOGtWvXBjsMpZQaUESk8931x2hVklJKqQ40MSillOpAE4NSSqkOBnwbQ1fa2tooLCykubk52KH0KDY2lpycHKKiooIdilJKHROWiaGwsJCkpCTGjBmDx2RbIcUYQ0VFBYWFhYwdOzbY4Sil1DFhWZXU3NxMenp6yCYFABEhPT095Es1SqnBJywTAxDSSaHdQIhRKTX4hG1iUEopT6v3VfDSukIaWpzBDiXkaWLowemnn97l8htvvJGXXnopwNEopbz19Ef7ufaJ1Xz/X5u49onVNLV2NT+SaqeJoQcff9zlDI1KqQFkb1k9D7yxg/PzsnjkulPZeqSGP7yzM9hhhbSw7JXkK4mJidTX12OM4Vvf+hbLli1j7Nix6Ii0Sg0cL6y2Rn64/6ppDE2M4cOdZfxt1UG+dtY4MpNigxxdaNISQx+88sor7Ny5ky1btvDkk09qSUKpAaLN5WbJxiOcPzmToYnWFO13nDOeFqebl9YVBjm60KWJoQ+WL1/OddddR2RkJMOHD+fcc88NdkhKqT74ZF8lFQ2tXD5zxLFl4zISmT16CP/ZcCSIkYU2TQx9pF1LlRp4lu0oJdoRwZm5Qzssv3zmcHaW1LGntC5IkYU2TQx9cNZZZ7F48WJcLhfFxcW8//77wQ5JKdUHy3aUcPr4dOKjOzanLpyYCcCqvRXBCCvkaWLogyuuuILc3FymTZvG7bffztlnnx3skJRSvdhXVs+BikbOnZR5wmcj0+IYnhLL6n2VQYgs9GmvpB7U19cDVjXSo48+GuRolFL9sWxHKQDnnHJiYhAR5o1L58NdZRhjtKq4Ey0xKKXC0rIdpUzMSmRkWnyXn88bn05FQyu7S+sDHFno08SglAo7Tpeb9YeqOH380G7XmT8uHbCGylAdaWJQSoWdXSX1NLe5OXVUarfr5AyJY0RqnCaGLgQsMYhIrIisEZFNIrJNRO7tYp2FIlIjIhvtx88CFZ9SKnxsKqwGYEZOarfriAhzxqbx6YGqwAQ1gASy8bkFONcYUy8iUcBHIvKGMWZ1p/VWGGMuDWBcSqkws+lwNSlxUYxO77p9od30nBRe2VBESW0zWck6PEa7gJUYjKW9lSfKfuigQ0opn9tUWMOMkam99jaanpMCwObCmkCENWAEtI1BRCJFZCNQCrxjjPmki9Xm29VNb4jIlG62c6uIrBWRtWVlZf4M2WuHDx/mnHPOIS8vjylTpvDwww8HOySlBoXGVie7SuqYaZ/0ezI5O4UIgc121ZOyBDQxGGNcxpiZQA4wR0SmdlplPTDaGDMD+CPwajfbecIYk2+Myc/IyPBnyF5zOBz8/ve/p6CggNWrV/PYY4+xffv2YIelVNjbdqQWl9swY2Rqr+vGRUcyMStJSwydBKVXkjGmGvgAWNRpeW17dZMxZikQJSLd9zcLYdnZ2cyaNQuApKQk8vLyKCoqCnJUSoW/TYerAZjeQ8Ozp+k5KWwpqtHh9D0ErPFZRDKANmNMtYjEAecDD3RaZxhQYowxIjIHK3GdVF+ye/+7je1Hak9mEyeYPDyZey7rsparSwcOHGDDhg3MnTvXp3EopU5UUFxHZlIMGUkxfVp/Wk4q/1xbSFF1EzlDem6sHiwC2SspG/iriERinfD/aYx5TURuAzDGPA5cDdwuIk6gCbjWDPA0Xl9fz1VXXcVDDz1EcnJysMNRKuztLq1jYlZSn9efPsJqi9hSWKOJwRawxGCM2Qyc2sXyxz1ePwr4dFCi/lzZ+1pbWxtXXXUV119/PVdeeWXQ4lBqsHC7DXtK6/nCaSP7/J1J2Uk4IoTNRTVcNC3bj9ENHHrns58YY7j55pvJy8vju9/9brDDUWpQKKpuorHV1a8SQ4wjkgmZiRQU+7bKeSDTxOAnK1eu5Pnnn2fZsmXMnDmTmTNnsnTp0mCHpVRY221PvJObmdiv703OTtbE4EGH3faTBQsWaC8HpQJsV4l1D21uP0oMYHUoeXlDERX1LaQn9q3ROpxpiUEpFTZ2l9STlRxDSlxUv76Xl211DCko1qk+QRODUiqM9LdHUrv2xLC9WG90A00MSqkw4XYbdpfUM6Gf7QsAaQnRDEuO1RKDTRODUiosFFU30dTWvx5JniYPT/b5zbADlSYGpVRYaO+RNDGr/yUGgLzsJPaW1dPc5vJlWAOSJgalVFho75E0IdPLEkN2Ck77BrnBThODnzQ3NzNnzhxmzJjBlClTuOeee4IdklJhbVdJnVc9ktrlZVsJRauT9D4Gv4mJiWHZsmUkJibS1tbGggULuOiii5g3b16wQ1MqLO0uqfe6fQFgdHoC8dGRbNcb3bTE4C8iQmKiVdfZ1tZGW1tbr7NJKaW80z5GUq6X1UgAkRHCKcOSNDEwGEoMb9wFR7f4dpvDpsFF9/e6msvlYvbs2ezZs4c77rhDh91Wyk/aeyTletnw3G5ydjJLNh3BGDOoL+S0xOBHkZGRbNy4kcLCQtasWcPWrVuDHZJSYWlXycn1SGo3eXgydc1OCquafBHWgBX+JYY+XNn7W2pqKgsXLuTNN99k6tTOs5kqpU7WyfZIajf52B3QtYxMG7xzM2iJwU/Kysqorq4GoKmpiXfffZdJkyYFNyilwtTu0jqGJcd63SOp3aRhyUSI9kwK/xJDkBQXF/OVr3wFl8uF2+3m85//PJdeemmww1IqLO0uqT/p9gWAuOhIxmUksk0TQ2CISCywHIix9/uSMeaeTusI8DBwMdAI3GiMWR+oGH1p+vTpbNiwIdhhKBX22nskXTdnlE+2Nzk7mXUHq3yyrYEqkFVJLcC5xpgZwExgkYh07tR/EZBrP24F/hzA+JRSA1BhVfsYSSdfYgCrAbqouonqxlafbG8gClhiMJb2e82j7EfnmWwuB/5mr7saSBURnYRVKdWtY7O2ncTNbZ6mDD/eAD1YBbTxWUQiRWQjUAq8Y4z5pNMqI4DDHu8L7WWdt3OriKwVkbVlZWV+i1cpFfqO90jyTYnh2NwMg7idIaCJwRjjMsbMBHKAOSLSue9mV3eUnDA/pjHmCWNMvjEmPyMjww+RKqUGit0lvumR1G5oYgxZyTGDugE6KN1VjTHVwAfAok4fFQIjPd7nAEcCE5VSaiDaVVrnkx5JnqaNSGFL0eCdzS1giUFEMkQk1X4dB5wP7Oi02hLgy2KZB9QYY4oDFaNSamBp75F0MoPndWV6Tip7y+qpb3H6dLsDRSBLDNnA+yKyGfgUq43hNRG5TURus9dZCuwD9gBPAt8IYHx+4XK5OPXUU/UeBqX8oLCqieY2N7k+al9oNy0nBWNg6yAtNQTsPgZjzGbg1C6WP+7x2gB3BCqmQHj44YfJy8ujtnbw1lcq5S/tYyT5qkdSu+kjUgDYUljDvHHpPt32QKBDYvhRYWEhr7/+OrfcckuwQ1EqLO061lXVtyWG9MQYRqTGsamw2qfbHSjCfkiMB9Y8wI7Kzk0ZJ2dS2iR+NOdHva73ne98h9/+9rfU1dX5dP9KKcueknqyU2JJjvVNjyRP03MGbwO0lhj85LXXXiMzM5PZs2cHOxSlwtau0jqf3b/Q2fScVA5WNA7KO6DDvsTQlyt7f1i5ciVLlixh6dKlNDc3U1tbyw033MALL7wQlHiUCjftPZKunzvaL9ufnmO1M2wurOGsiYPrfiktMfjJfffdR2FhIQcOHGDx4sWce+65mhSU8qHDVY00t7l9NkZSZ1PbG6AHYXWSJgal1IC020eT83QnJS6KsUMT2HS42i/bD2VhX5UUChYuXMjChQuDHYZSYWXHUasL+CnD/JMYwLoD+tMDlX7bfqjSEoNSakAqOFrHqLR4EmP8d307PSeF4ppmSuua/baPUKSJQSk1IO0ormWSH0sLAKeOGgLA+kE2cY8mBqXUgNPc5mJ/eQOT7CGy/WXaiBRioyJYvW9wVSdpYlBKDTi7S+pxG8jzc4kh2hHBrFFDWLNfE4NSSoW0Arvh2d8lBoA5Y9MoOFpLTVOb3/cVKnpttRGRvs6wXW2M0ZHilFJ+V1BcS1xUJKPS4v2+rzlj0zAG1h2s5NxJWX7fXyjoS3P+X7FmUetqdrV2BngO+JsPYgobY8aMISkpicjISBwOB2vXrg12SEqFhR3FdUwclkRkRE+nJd+YNWoIUZHCJ/s1MRxjjDmn8zIRGWaMOeqfkMLL+++/z9ChQ4MdhlJhwxjDjqO1fGbKsIDsLzYqkhk5qXwyiBqgvW1j+LJPo1BKqT4qrWuhqrHN711VPc0fn86WohpqmwdHO4O3d4ZcLiKNWLOw7fRlQL529De/oaXAt8Nux+RNYtjdd/e6nohw4YUXIiJ8/etf59Zbb/VpHEoNRgXFgWt4bnfGhKH8cdkeVu2tCFhJJZi8LTFciTX95hUi8lRfviAiI0XkfREpEJFtInJnF+ssFJEaEdloP37mZXwhYeXKlaxfv5433niDxx57jOXLlwc7JKUGvIJia36TvGGBSwyzRg0hPjqSj3aXB2yfweRVicEYUwK8aT/6ygl8zxizXkSSgHUi8o4xZnun9VYYY3w2QXJfruz9Zfjw4QBkZmZyxRVXsGbNGs4666ygxaNUONhxtJbhKbGkxPt+cp7uRDsimDs2jZV7QiMxGGO44elPuHhatl+GHfeqxCAij4nIc/brC/vyHWNMsTFmvf26DigARniz/4GgoaHh2MxtDQ0NvP3220ydOjXIUSk18BUU1/p14LzuLMjNYF95A0XVTQHfd2d7SutZuacCt/HP9r2tSmoF9tmvz+3vl0VkDHAq8EkXH88XkU0i8oaITOnm+7eKyFoRWVtWVtbf3QdESUkJCxYsYMaMGcyZM4dLLrmERYsWBTsspQa0plYXe0rrj82VEEhn5lq9Cz/aHfxzztvbSwC4IM8/3We9bXxuBFJEJAro6w1wAIhIIvBv4Dtd3BC3HhhtjKkXkYuBV4HcztswxjwBPAGQn5/vp5x5csaNG8emTZuCHYZSYWV7cQ1uY41hFGi5mYlkp8Ty/o4yvnBav057Pvf29hJm5KQwLCXWL9v3tsRQCewFHgNW9vVLdiL5N/B3Y8zLnT83xtQaY+rt10uBKBHRmwCUUgBsKbRmU5uWE/jEICKcl5fJ8t1lNLe5Ar7/dkdrmtl0uJoL/dg7ql+JQURSReRZ4Cp70d+A/D5+V4CngQJjzB+6WWeYvR4iMseOr6I/MSqlwteWolqGJkYzLNk/V8q9OS8vi8ZWF6v2Be+09E6BVY104WT/3YXdr6okY0y1iNwPjAHKgenACVf+3TgD+BKwRUQ22svuxq6KMsY8DlwN3C4iTqAJuNYYE5JVRUqpwNtSVM20ESnY148BN39cOvHRkbxXUMI5p2QGJYZ3tpcwdmgCEzL9M9c1eNfGcDOw3xjzFrCur18yxnxEz+MtYYx5FHjUi5iUUmGusdXJntJ6FgXxBrPYqEjOzB3Ku9tL+eXlJuAJqra5jVV7y/nqGWP9um9v2hiqgNtE5CERuUlETvV1UEop1VlBca3V8JyTGtQ4zs/L4mhtM5vt9o5A+mBnGW0uwwV+rEYCLxKDMeY+4GvAz4H9gN6xpZTyu/YTcTB6JHm6YHIWjgjh9S3FAd/3ko1HyEyKOTblqL/0OzGIyC+Ay4ELgCJjzMM+jypMVFdXc/XVVzNp0iTy8vJYtWpVsENSasDaeLiazKQYspJjghpHanw0Z+YO5fXNxQSyCbSyoZUPdpbyuVNH+H24cW9KDD8DWuzvXiUiT/o8qjBx5513smjRInbs2MGmTZvIy8sLdkgqzBytaeYbf19H3v97k/P/8CHrDobv0NDrD1Uxe/SQoDU8e7p0+nCKqpvYcLg6YPt8bfMRnG7DFaf6f8AIb+9jeAbIA9KBP/kunPBRW1vL8uXLufnmmwGIjo4mNTU1uEGpsNLqdPOlpz/h/R1lXDFrBG0uN196eg27S+qCHZrPldY1c7iyiVl+rkLpqwumZBEdGcFrmwJXnfTy+iImDUsiLwCjynp75/O3sYbFcAAPE8LtDCv+uYvyw/U+3ebQkYmc+fmJPa6zb98+MjIyuOmmm9i0aROzZ8/m4YcfJiEhwaexqMHr6Y/2s7u0nmduzOfcSVkcrWnmkkdW8K0XN/D6t88MyOxmgbL+YDUAs0anBjWOdsmxUSw8JYMlm45w10WTiHZ4e43dN3tK69l4uJq7L57k1/208/an2QvEAv8xxoRsUggmp9PJ+vXruf3229mwYQMJCQncf//9wQ5LhYmmVhdPrtjHwlMyjk03OSwllnsvn8KOo3X8d9ORIEfoWxsOVREdGcGU4cFtePZ03dxRlNe38OY2/09m+cLqg0RFClfOyvH7vsD7EsM24DBws4j8zhhzmg9j8qneruz9JScnh5ycHObOnQvA1VdfrYlB+cwrG4qobGjljnMmdFh+8dRsJg3bwyPLdnPZjOFhU2pYf6iKKSOSiY2KDHYox5ydm8Ho9HieX3WAz84Y7rf9NLY6+ff6Qi6els3QxMA0vHtbYhiPlVSeAG7yXTjhY9iwYYwcOZKdO60J7t577z0mT54c5KhUuHh1QxETMhPJH92xzj0iQrjzvFz2lTWETamh1elmU2FNyLQvtIuIEG6YO5pPD1Sx/Ujn8UB9Z8nGI9Q1O/nSPN/Pu9AdbxPDYWPMEqxZ3Ap8GE9Y+eMf/8j111/P9OnT2bhxI3cHcdIgFT4KqxpZc6CSz80c3mUPnc9MGUZuZiJPf7Q/CNH53tYjNbQ63cweHVqJAeCa/BxioyJ4asW+3lf2gttteO7jA0walhTQn9/bxLBIRHKAx4EHfRhPWJk5cyZr165l8+bNvPrqqwwZEnp/2Grg+c9GqyRw+cyuuy1GRAhfmj+aLUU1bC6sDmBk/rFqrzVg3ZyxaUGO5ESp8dHcMHc0r24sYn95g8+3//b2EnYcrePWs8YFtJuut4khFfgR8EOsexqUUgFgjOE/G4vIHz2EkWnx3a73uVNHEOOI4OX1RQGMzj9W76tgYlZiwOrX++vWs8cRFRnBo8v2+HS7xhgeeW83Y9Lj/dqG0RVvE8MvsHok7QSCNzC5UoPM3rIGdpXU89mZPZ8o2rtTvrG1GLe/5n8MgFanm7UHqpg/Lj3YoXQrMymW6+1Sw8EK35Ua3tpWwvbiWr55bi6OSP92h+2sz3sTkRntr40xhcaYd+3Xd/kjMKXUiZbtsMbiP68PUzpePC2bktoW1h2q8ndYfrO5sJqmNhfzx4duYgC47exxOCLEZ6WGxlYnv3xtOxMyE/lcLxcB/tCfNLRBRDaLyA9FZKTfIlJKdeu9glLyspMZkRrX67rn5WUR7Yjg9c2BH+zNV1btrUAE5o4N7cSQmRzLF+eO4uUNRWzywTAZD7+3m6LqJn5zxbSAlxagf4nh90ACcD+wX0TeF5Gv+icspVRnNY1trD1YxXmT+jZBTGKMg4UTB3Z10qp9FUwalsyQhOhgh9Kr75w/kYzEGL73r00nNfXnhkNVPLViP9eeNjJoDe59TgzGmB8YY8ZjTeX5FNYwGE/4KzClVEcf7CrF5Tacm9f3mcMumW5VJ60fgNVJDS1O1h6oYsGE0C4ttEuJi+K3V09nT2k9D7y5w6ttVDe28s1/bCA7JZYfXxy8QTf708aQLiK3AL/BuqlNsO5+7uv3R9qljAIR2SYid3axjojIIyKyx662mtXX7YeanTt3MnPmzGOP5ORkHnrooWCHpQawZTtKSU+IZmY/Jqo5Vp0UhLkDTtZHe8ppdbmPDfkxEJw1MYMbTx/DsysP8PL6wn59t8Xp4tuLN1Ja18xjX5xFSlyUn6LsXX+GxDiKlUiqgGeBF+zpOvvKCXzPGLNeRJKAdSLyjjFmu8c6FwG59mMu8Gf7ecA55ZRT2LhxIwAul4sRI0ZwxRVXBDcoNWC53YYVu8tZODGDiH4Mc5EY4+DsiRm8seUo/++Syf36brAtKyglKdZB/piBdf/PTy7JY+fROu56eQvDUmI5ffzQXr/T4nTxjRfWs3xXGQ9cNY0ZI1P9H2gP+tPG8ApwBZBtjLmtn0kBY0yxMWa9/boO647pznfoXA78zVhWA6kikt2f/YSi9957j/HjxzN6dOBuaVfhZXtxLZUNrZw5sfeTTGeXTMvmaG0zGw4PnOokt9vw/s5SzpqYQVQQGl9PRlRkBI9dP4vRafF85Zk1/O9bO2lsdXa7flF1E9c+sZr3dpTyq89N5QunjQpgtF3rc4nBGPN5X+1URMYApwKfdPpoBB2rpwrtZR3KwSJyK3ArwKhRPR/E9597gtKDvr1dPXP0OM658dY+r7948WKuu+46n8agBpflu8sAOGNC/xPDeXmZREdG8Prmo8weHXp3D3dl25FaSutaOPeUvrenhJK0hGj+ddt8fr5kG4++v4d/rDnENbNzuGzGcCbb8ykcrmrk3+sKee7jA7gN/On6WVw8LTSug70dXdVrIpII/Bv4jjGm88hTXZVzT+hOYYx5ArvhOz8/P6S7W7S2trJkyRLuu+++YIeiBrAVu8rJy04mMym2399Nio3irIlDeWNrMT+9JG9AVCct21GKCCw8JSPYoXgtNT6ah649lRvmjebJFft46qP9/GX5PpJiHLS53TS3uRGB8/Oy+MnFeYwZGjpztfQ7MYjIZcaY/3qzMxGJwkoKfzfGvNzFKoWA5z0SOcBJDRHZnyt7f3jjjTeYNWsWWVkDpwFNhZbGVidrD1Zy0xljvd7GxdOyebeglA2Hq0NyMLrOlu0oYebIVNJDdBiM/sgfk0b+mDTK61tYvquMDYeqiXFEMD4zkfnj0kMqIbTzpsTwa6DfiUGsEaCeBgqMMX/oZrUlwDdFZDFWo3ONMWbgdafw8OKLL2o1kjopn+yvpM1lODO3/9VI7c6fbE1FuXRLccgnhrK6FjYV1vD9C4Mzl4q/DE2M4cpZOQGbbOdkeNOq42059AzgS8C5IrLRflwsIreJyG32OkuBfVjDeT8JfMPLfYWExsZG3nnnHa688spgh6IGsI92lxPtiOC0Md63DyTHRrEgdyhvbj2KMSFd+8p7BdawH+f08UY+5XvelBi8+quyezH1mFSM9Rd7hzfbD0Xx8fFUVFQEOww1wK3cU07+6CEnPXvZoqnDWLajlC1FNUzvx70QgbZk0xHGDk041kirAm9g9QNTapApr29hx9E6r3ojdXZBXhaREcKbW/0/R7G3SmqbWbWvgs/O6HoSIhUYmhiUCmHtk9Sc7oPRRYckRDN/XHpIVye9trkYY+h1WHHlX94khhKfR6GU6tLHe8tJinEwbUSKT7b3manD2FfewM6SOp9sz9eWbCxi6ohkxmckBjuUQa3ficEYc4E/AlFKnWjlngrmjkv32dDLF00dhiNCQnJmtwPlDWwqrAn4bGXqRFqVpFSIOlzZyKHKRs7w4eiiQxNjOC8vk3+vK6TV6fbZdn1hySbrlqVLp2tiCDZNDEqFqI/3lgPeDYPRk2tPG0VFQ+ux2eBCgTGGVzcWMWdMGsP7MAmR8i+vEoOIfNfj9Sm+Cye8PPjgg0yZMoWpU6dy3XXX0dzcHOyQ1ACyck8FGUkx5Gb6tr79rIkZDEuOZfGnfR413+9W7atgX1kDXzhNJ4cMBf1KDCKSKiLPAteIyDdEZAGgcz53oaioiEceeYS1a9eydetWXC4XixcvDnZYaoBwuw0f7y3n9PHpPu+2GRkhXJOfw4e7yjhS3eTTbXvrhdUHSY2P4pLpoTGI3GDXr8RgjKk2xtwE/BxrZNRcoKsxjxTgdDppamrC6XTS2NjI8OFad6r6ZntxLeX1rZyV659B5K6ZPRJj4KV1/ZtMxh9Kapt5a1sJX8gfedI38Snf8HZ01TZjzDoROQKU+jIgX6v+715ajzT4dJvRwxNIvWx8j+uMGDGC73//+4waNYq4uDguvPBCLrzwQp/GocLXh7usYbbPmuifxDAqPZ4zJqTzj08OcdvZ44l2BK+58cU1h3AbwxfnBn8eAmXx9q9hkYjkAI8DD/ownrBRVVXFf/7zH/bv38+RI0doaGjghRdeCHZYaoD4YGcpU0ckk5Hkv9FFbzlzHEdrm4/1BgqGNpebF9cc4uyJGYxOD71RRgcrb0sMqcCPgB8Ct/gsGj/o7creX959913Gjh1LRoZ1xXfllVfy8ccfc8MNNwQlHjVw1DS1sf5QNbedPc6v+1k4MYNJw5L4y4d7ufLUEUGZp+GtbUcpqW3hN1fo7IahxNsSwy+AV40xOwGXD+MJG6NGjWL16tU0NjZijOG9994jLy8v2GGpAeDjPeW43IaFfp69TES47ezx7C6tZ9mOwNcIG2N47P29jMtI8PvPqvrH28TwY6whtAHe91EsYWXu3LlcffXVzJo1i2nTpuF2u7n11uBOGqQGhg92lpEU6+DUAEwIf+n0bEakxvH4h3v9vq/O3tp2lILiWr6xcAKRA2BWucHE28TQijVvAsA5Pool7Nx7773s2LGDrVu38vzzzxMTM/Bno1L+ZYzhw11lnJk71GfDYPTEERnB184cy9qDVazeF7gh4ttcbh54cycTsxL5nA6YF3K8/ctrBFLsqTq1K4FSPrKrpJ6jtc2c7afeSF35wmmjyE6J5Vevb8ftDsyoqy+uOcT+8gbuumhSQBKg6h9vfyP3AHuBx4C/+y4cpQa3D3Zadf3+6qbalbjoSO66aBJbi2p5ab3/72uoa27j4Xd3M29cGudo20JI8jYxfNsY8ydjzK1Y03D2SkSeEZFSEdnazecLRaTGY9rPn3kZm1ID1oe7ypg0LInslMCOF/TZGcM5dVQqv31zJ9WNrX7d1+Mf7qWioZW7L87TyXhClDdDYrwIXG0PiXEGfR8S4zlgUS/rrDDGzLQfv+hPbEoNdA0tTj49UBnQaqR2IsIvL59KdWMr9/53u9/2s7ukjieW7+OKU0eE9PSig12/h8QACoHngdXARPo4JIYxZjlQ2c/4lBo0Vuwuo81lOPuUwCcGgKkjUvjGORN4ZUMR72z3/cirbrfhrpe3kBDj4CeXaNftUOZNVVIFcBvwZfu9Lysl54vIJhF5Q0Sm+HC7SoW8pVuOkpYQzZwxaUGL4ZvnTCAvO5m7X9ni8yqlvyzfx7qDVfz0kskMTdQeeqHMmxnc7ge+hjWQ3n7gTB/Fsh4YbYyZAfwReLW7FUXkVhFZKyJry8rKfLR733v44YeZOnUqU6ZM4aGHHgp2OCqENbe5eK+ghM9MyQpqL51oRwT/e810qhpauevfW3w2N/Tb247y27d2cMn0bK6aNcIn21T+0++/QBH5BXA5cAFQZIx5xBeBGGNqjTH19uulQJSIdDlDiTHmCWNMvjEmv33IiVCzdetWnnzySdasWcOmTZt47bXX2L17d7DDUiFqxe5yGlpdXDQ1+MNOTxmewo8WTeLNbUe5/80dJ50c1h+q4lsvbmB6Tir/e/UMbXAeALwpMfwMeASoA64SkSd9EYiIDBP7L0ZE5tixBe6OGx8rKChg3rx5xMfH43A4OPvss3nllVeCHZYKUW9sLSYlLor54303jefJuOXMsdwwbxR/+XAfv3q9wOvksP5QFTc9+ynDUmJ55iv5xEXrsNoDgbeD6H0d+Isx5s2+fsHuzbQQGCoihVj3QkQBGGMeB64GbhcRJ9AEXGt8UI594403OHr06MlupoNhw4Zx0UUX9bjO1KlT+clPfkJFRQVxcXEsXbqU/Px8n8ahwkOr080720v4zJRhRIXIzV7tvZQcERE8/dF+Wp1u7v3slH4NtLdsRwnf/McGMpJieOHmuaRru8KA4W1ieAbrJJ4A/N0Ys7G3Lxhjruvl80eBR72MJ+Tk5eXxox/9iAsuuIDExERmzJiBw+Ht4VbhbOXecuqanVw8bViwQ+lARLjnssnEOCL4y/J9HKxs5AcXnsK0nJQev1da28wDb+7k3+sLmZydzHM3nUZmcmyAola+4O2Z6ttY4yU5sKqVzvJZRD7W25W9P918883cfPPNANx9993k5OQELRYVupZuLiYpxsEZE7psUgsqEeGuiyaRlRzLw+/t5rJHP+LM3KFcP3cUk4YlkzMkjoZWFwcrGth5tI4PdpbxTkEJGLh94XjuPC9XZ2UbgLxNDHuxpvX8jzHmf3wYT1gpLS0lMzOTQ4cO8fLLL7Nq1apgh6RCTFOrize2HuUzU4YR4wjNE6iI8NUFY7k6P4fnVx3kuY8PcNsL6wGIEPAcXik9IZovzhnFjaePYcxQnXhnoPI2MWwDDgM3i8jvjDGn+TCmsHHVVVdRUVFBVFQUjz32GEOGDAl2SCrEvLmtmPoWJ9fkh35pMjk2ijvOmcAtZ45la1EN+8oaOFjRSHKcg1Fp8eRmJTEmPUGH0A4D3iaG8UAV8IT9rLqwYsWKYIegQtw/Py1kVFo8c8cG76a2/opxRDJ7dBqzRw+cmFX/eNsF4rAxZgnWAHoFPoxHqUHjcGUjq/ZVcPXsHO3br0KKt4lhkYjkAI8DD/owHqUGjZfWFSICV80O/WokNbh4mxhSgR8BPwRafBaND/nqVn5/GggxKv9wuw0vrStkwYShjEgN7BDbSvWmz4lBRGZ4vP0FVo+knYDL51GdpNjYWCoqKkL6xGuMoaKigthY7d89GK3cW05RdRNXa2lBhaD+ND5vsCfZeQF40RjzLoAxpq/zMQRMTk4OhYWFhPIAe2AlML23YXB6duUBhibGsGhqaN3UphT0LzH8HrgSuB/4jYisAJ43xjzjl8hOQlRUFGPHjg12GEp1aV9ZPct2lPKd83ND9t4FNbj1uSrJGPMDY8x4IB94Cutu5yf8FZhS4eqvHx8gOjKC6+eODnYoSnWpzyUGEUkHrsAa7O4cQIBDfopLqbBU1dDKv9YVcumMbDKSdFA5FZr6U5V0FKuEUQU8C7xgjPnIL1EpFaae+/gAja0ubjt7fLBDUapb/UkMr2A1PL9hjGnzUzxKha2GFid/XXWA8/OymJiVFOxwlOpWr4lBREbZL79vP2d3c5dmtTGm1leBKRVu/vHJIaob2/jGOVpaUKGtLyWGvwLtNwR0d9++AZ4D/uaDmJQKO9WNrTz6/h7OzB3KrFE6mKIKbb0mBmPMOYEIRKlw9sCbO6hrbuMnl+QFOxSlehWweQRF5BkRKbVvkuvqcxGRR0Rkj4hsFpFZgYpNKX9ava+CF9cc5pYzxzFpWHKww1GqV4GcYPY5YFEPn1+ENflPLnAr8OcAxKSUXzW3ufjxy1sYlRbP/5w/MdjhKNUnAUsMxpjlQGUPq1wO/M1YVgOpIpIdmOiU8o8/LtvN/vIGfnPFNOKi9S5nNTAEssTQmxFYs8K1K7SXKTUgrT1QyeMf7uPq2TksyA29+ZyV6o63M7j5Q1c9nrocHlVEbsWqbmLUqFFdraLCSFldC+8WlFBc3cTYjAROHz+UrOTQHpW2prGNOxdvZERqHPdcNjnY4SjVL6GUGAqBkR7vc4AjXa1ojHkCe5ym/Pz80B1bW52U+hYn9y0t4B9rDuE5gnqEwL2XT+VL80JzrKE2l5tvL95ASW0zL91+OkmxUcEOSal+CaXEsAT4pogsBuYCNcaY4iDHpIJka1ENX39+HUdqmvjyvNFcN3cUuZlJFBTX8uA7u/h/r26ltqmNbywcH1LTYhpjuPvlLXy4q4z7r5zGzJGpwQ5JqX4LWGIQkReBhcBQESkE7gGiAIwxjwNLgYux5pFuBG4KVGwqtOwvb+C6J1aTFOvgpdtOZ/bo4zeETR2RwuNfms0PX9rM797aSU1TGz++aFLIJIcH39nFv9YVcud5uVw7R6s51cAUsMRgjLmul88NcEeAwlEhqqnVxe0vrCMyUvjnbfPJGRJ/wjpRkRH8/poZJMU6eGL5PgT48cXBv3HsH58c4pFle/hC/ki+c35usMNRymuhVJWkFL96fTs7S+p49sbTukwK7SIihHs/OwW3Mfxl+T5ys5KCOk3mPz89zE9e3cLCUzL41RVTQ6YEo5Q3NDGokLF6XwV//+QQtywYy8JTMntdX0T42aVT2F/ewA9f2kRsVASXTh8egEg7emH1QX766lbOmpjB4zfMJioylHqBK9V/+hesQkKr081PXtlCzpA4vnfhKX3+XrQjgie/nE/+6DTuXLyRt7cd9WOUJ3p25X5++upWzpuUyRNfmk1slN7EpgY+TQwqJDz90X72ljVw72en9PsO4fhoB0/fmM+0ESnc8Y/1vL+j1E9RHtfidHHPf7Zy73+3s2jKMP58gyYFFT40MaigK6xq5JH3dnPh5CzOy8vyahtJsVH89aY5TMxK4pa/reWx9/fQ6nT7OFLLgfIGrvrzx/x11UFuXjCWP37xVKId+q+kwof+Naug+8V/twPws5O8QzglPorFt87jwslZ/O6tnSx6eDkrdpf5IkQAXG7DUyv2sejh5RyubOLJL+fz/y6drG0KKuzoX7QKqvcKSnh7ewnfPi+3x15IfZUUG8Wfb5jNszeehttt+NLTa7jt+XUcrmz0epttLjfLdpRw9eMf86vXC5g/Lp037jyTCyZ7V7pRKtRpryQVNE2tLu5Zso0JmYncvGCsT7d9zqRMTp+QzlMr9vPosj28U1DCxdOyuXp2DvPHpfda9eNyG/aU1vPPtYd5eX0hVY1tZCTF8NAXZnL5zOHaHVWFNU0MKmj+9MEeCquaePFr8/xSRx/jiOSOcyZw1awcnlm5nxc/OcR/Nx0hKcbBvPHpTMxKZEJmImPSE4iKjKCu2cmesnre2nqU9YeqaGx1ERkhXDR1GJ+dMZyFp2RqW4IaFMSYgT0GXX5+vlm7dm2ww1D9tLesnoseWsHF04bx0LWnBmSfzW0uVu4p561tR1l7sIqDFY243Cf+/Y9Oj2fhxAxmjExlztg0n1RxKRVqRGSdMSa/q8+0xKACzhjDPf/ZRowjgrsDOAdybFQk5+Ud7/nU6nRzsKKBQ5VWgoiPdjA6PZ6cIXFaVaQGNU0MKuBe21zMR3vKufezU8hMCt68CtGOCHKzksjNSgpaDEqFIq0wVQFV29zGL1/bzpThydwQovMpKDXYaYlBBdR9Swsor2/hiS/nExmh1TVKhSItMaiA+Wh3OS+uOczXzhynE9goFcI0MaiAaGhx8qN/b2bc0AT+54KJwQ5HKdUDrUpSAfHAmzs4UtPEv74+XwebUyrEaYlB+d0n+yr426qD3Hj6GPLHpAU7HKVULwKaGERkkYjsFJE9InJXF58vFJEaEdloP34WyPiU7zW2WlVIo9Li+cFn+j7PglIqeAJWlSQikcBjwAVAIfCpiCwxxmzvtOoKY8ylgYpL+U+by83tL6znUGUjf79lHvHRWnOp1EAQyBLDHGCPMWafMaYVWAxcHsD9qwD737d28uGuMn5zxTTmj08PdjhKqT4KZGIYARz2eF9oL+tsvohsEpE3RGRKVxsSkVtFZK2IrC0r8914+8p3XlxziL8s38eX5o3m2jmjgh2OUqofApkYurqbqfMIZuuB0caYGcAfgVe72pAx5gljTL4xJj8jI8O3UaqT9u72En7yyhYWnpJx0pPvKKUCL5CJoRAY6fE+BzjiuYIxptYYU2+/XgpEicjQwIWoTtb6Q1V888X1TBuRwp+un6Wzmyk1AAXyv/ZTIFdExopINHAtsMRzBREZJvawliIyx46vIoAxqpOwt6yem5/7lGHJsTx942na2KzUABWw/1xjjFNEvgm8BUQCzxhjtonIbfbnjwNXA7eLiBNoAq41A33CiEGitLaZrzyzhsgI4a9fncPQxJhgh6SU8lJAL+ns6qGlnZY97vH6UeDRQMakTl5pXTNffmYNlQ2tLL51HqPTE4IdklLqJGhZX52UNfsr+faLG6hpauPJL+czPSc12CEppU6SJgblFZfb8Kf39/Dgu7sYlRbPU1+Zz9QRKcEOSynlA5oYVL+V1jXzP/+3kZV7Krh85nB+fcU0EmP0T0mpcKH/zarP2lxu/vHJIR58dxfNbS5+e9V0rsnP0fmRlQqG1kZwtUJcqs83rYlB9crtNry2pZgH39nF/vIGTh+fzr2fnaJzJSsVSMZAxR7Y/Q7seRcOfAQLvgPn3O3zXWliUN0yxvBeQSm/f2cXBcW1TBqWxNNfyefcSZlaSlDK39qa4OgWKFoHhWuhcA1UH7I+GzoRTrsFJlzgl11rYlAnaHG6eH1zMX/5cB87S+oYnR7Pw9fO5LLpw4nQeZqV8o/6MihaC4WfwoGV1mu30/osaTjk5MMZ34EJ58OQ0X4NRRODAqC5zcXyXWW8sfUo724voa7FySlZSfzh8zO4bMZwHdpCKV9wtUHlfijfCeW7oHy39Vy5D5qqrHUkEobPhPl3QM4cGDELkocHNExNDINYc5uLD3aWsXRLMe8VlNDQ6iI1PoqLpg3j0unDOTN3qFYZqY7cLmiuseq72/82REAiALGXCWDAuK31IiLtz7HeY+xnICbp+HYGMmOgpQ6aKqGx0n6uOv6+sQIOrYayguOlALBKAkNzYfLnIH0CjJgN2TMgOr6LXRianE3UtNRQ3VJNdUs1WQlZjEsZ5/MfRxPDIFNa18yHO8v4YGcZ7+8spbHVxZD4KC6bMZyLp2Uzf3y6lg7CidsNjeVQVwx1R8HZApFREBEFkQ5obbDWczZDU7V11dpUBc3VHu+roaXWWre52jrh+4oj1koakdGQmAXGZW3f7bIShiMOojwejliIioeo9uc46yq8ucbahiPG+vmi4iE+3fpcIjwSV8TxRHbsudPnGOt4OFush6vV+tkby6HZPg6tDdBmPzeUWcfW1dr9zxmbgjtrCu7Tv4U7fQJNQ0bTlJxNnRiqW6qpaK6gsqmSmuqN1JYsp7a1lpqWmhOe29xtHTZ709Sb+O7s7/ru99H+a/H5FlVIKa5pYuWeClbuKWfN/kqKqpsAyEiK4fKZI7hkWjbzxqXh0GQQuoyB1nrrJNVSa12JNlbYV6MV9sN+3VR1/H1rPbQ19v9E7oiDuCH2IxXSxkJMsnWSjU+H+DTrBNrh6t+jhIDxONly/EQPx0sUIta69SXW8tYG6+eRyOMlDGPA2QRtzVZDbFOV9dz+cNrPEgGxqeBuA2erdYJ2tfjiyHcUGQNxQ3BFx1ERHU91dBwmJg6TMo2quNMpdkRyFBfNkVG0RERSbVqpdDZR0VpLRXMF1S1HcBf9G4p63k1CVAIp0SkkxySTEp3C+NTxJEcnkxyTTGpMKinRKdZzTAo5STm+/znRxBBWjDEUVTexpbCGj/dWsHJvOfvKrCvC9IRo5o1L54Z5ozkzdyiTs5O1ITnQnK3W1WWHq/Kq41flnZc311on96aqjtUPnUkExKUdP2mnjbOqJNpP5olZkJwNSdnWFbWrzdqe22lfsYv1HDfEOsFGxQbmePiTy2klGmdzx8Rk2hOY9XAZF3WtddS21VHbWk9tWz21bY3UulupdTdT62ym1tlEjbOB0uYKjjaWUNFUgcvUAXXWjDJN9sMWHRFNTGQMyTHJpMelMyJpBNMzppMWm0ZMZAwiQrwjnjhHHInRiSRHJzM0bihpsWkkxyQTFREVhAPWkSaGAcrlNpTUNnOospFNh6tZd7CKdQerqGiwirPx0ZHMGZvGF+eM4owJQzklK0kTgb+017uXbofSAqtxMS4VIhxQtd96X7kPao9w4txUtogo6zvtV+qJw6wuidGJx6/cHXEQnWAnADsJxKdBTApEDO4SX7Oz2TrBt9Zajxbr+WjDUYrqi2h0NtLibKG2tZa61rpj69a31fe43eiIaJJjkkmOTiYzPpP52fPJjM8kKz6LIbFDiLBLRakxqWQlZJGdkI0jYuCfVgf+TzBINLQ4WbO/kqVbilm9v4Li6mac7uMnmdHp8Sw8JZOZo1KZMjyZqcNTiHYM7pPFSXG7rKv4Y1U15dZzQxnUFEFtkXWirym0rvA9OeKsag6AhAzrCn7MmVYXw+ThHtU09hV63BDrhB8OjbC2Vlcrlc2VNLQ1ECmRNLQ1UNVSRaREEuuIRRCaXc20OFtocbXQ6m6l1WU92txtON1O2txtNDubaXI20dDWYD2cDTS2NdLY1kh9W/2x5851757SY9NJiEogxhFDUlQS2YnZnBJ9CsnRySRFJ5EUnURydDIpMSlWlY1dbZMcnUysIwxKT17QxBDCqhtbeX7VQd4tKGHrkVpcbkNSjIOzJmbw2RnxDE+NY0RqHFNHpOj8B33hclon+PpS+1ECDR6v25c3llvVN93VzcelQcoISBkJo+ZZJ//oBMjIg6zJVpWNy67rjhm4d4fXtdZR3VxtNXy21hy7ym5obTh+Ine3Utti1aFXNVdR2VxJVXNVr1fifRUhEcQ74omPiichKoEER4JVB5+YQmJUorUsKoHEqMQOJ/T2k31mfCbxUSf28FE908QQgmqb2/jnp4f5y/J9lNW1kD96CN9YOJ7TxqQxd1waMY7IYIcYXMYc7yXSUm81NLY1WVfz9WUeJ/tS+7Vdr99a1/X2ohMhMRMSMiFjIiScAfFDO1bZJHi8j4rrPUZHjPUIMXWtdVQ0VVDdUt2h22NNSw2VzZVW75jmSsoayyhuKO51e+1VLWmxaaTFpjF16FTSY9MZEjuEtNg0EqMScRonCY4EhsQOwWVctDhbMBhiImOIdcQSHRl9rF4+KjKKqIjjD0eEQ7tMB4EmhhBSXNPEXz7cxz/XHqax1cWcMWk8e+Np4Tuctctp9bJp723iarEGBjtWfWM/Gso9qnMqrfc99TqRCOvEnphpXc2nT7BO6DHJkJhhJYDELOvzxEzraj8EGGNocbXQ6LSqShqdjTQ5m469bmhroKalhpqWGpqcTVa9elsdlc2VVDZX0uxsps3dRpurjciISJKjk0mISsAYQ5u7jfKmcsqayrrcd6REMiR2COmx6aTFpjEyayQTUieQGZ95rHql/So8MTqR6MhoHKIn7XAV0MQgIouAh7Gm9nzKGHN/p8/F/vxioBG40RizPpAxBprT5eatbSU8v/oA6w9W4zaGz84YzlcXjA2dhGCM1ZPlWN/uJvvZft/m+d5+tDZYDbLNNdbJv7nWWt5cY9XR93QF7ykmBRLsK/XkHBg2o+MVfEyS1aPGEWMngyzr8wjvSlXGGNzGjRs3xhhcxtXlMrdx0+pqpcnZ1OHR6Gw8Vi9+7NHWcZ0m1/FlnZOAuw9dSwUhzhFHrCOWhKgE0mLTyEnMIc4RR3RkNFERUTjdTupa66hvqydCInBEOMgdksu4lHFkxmeSGpN67JESa1XLtDekKhWwxCAikcBjwAVAIfCpiCwxxmz3WO0iINd+zAX+bD8HRavTzYGKBlqd1j/r0ZpmqhpbiRAhIgIiRBARIkWIEBD7OTJCiBDBESkkxjjIy04mNso6UdU0trFqXwWf7q9g/+FCqkoLiW2tYGpSM18eH8G8CRmkJZZD2XYot/t0Rzism3faGq0r5pZa66o4IrJjv2+JsE6Kjlirfru973dTlfXe7bKej53gPZ7bmjst7/Tey5uaTFQCrpgknLFJtDliaItJwpkzi7a4FNpiEnFGJ9AWFU9LpINGgWaBJkc0zZFRNEVE0GycuNyuYydml9uFy7hoaKujpbEcd4Mbt93t0POk3dXDZVzHGjadbifNruYOJ/FWVyumu15DJyE2MpY4R1yHR3xUPFnxWcdetz+316fHOeKOvT5Wx+5IIDnGunLXk7jyp0CWGOYAe4wx+wBEZDFwOeCZGC4H/maMMcBqEUkVkWxjTO+Vnf30wo9/Q0mAeu0s62JZPDCl/U0k0BhNwX4o2F8FVJ3E3iq6WOYgaLWGLqAZqOlppRb70a6h183G4SAurGpCDa000EoD1cEORQ0Yw1qjuf6BgX3n8wjgsMf7Qk4sDXS1zgigQ2IQkVuBWwFGjRrlVTBRMQ4SWrR+NDC6Oc59PPz6W1Kqa3Hx/mkfC2Ri6Or/u3O5vS/rYIx5AngCID8/36uy/xd+/kNvvqaUUmEvkBWVhcBIj/c5wBEv1lFKKeVHgUwMnwK5IjJWRKKBa4ElndZZAnxZLPOAGn+0LyillOpewKqSjDFOEfkm8BZWc+szxphtInKb/fnjwFKsrqp7sLqr3hSo+JRSSlkC2q3DGLMU6+Tvuexxj9cGuCOQMSmllOpIO0MrpZTqQBODUkqpDjQxKKWU6kATg1JKqQ7Eau8duESkDDjo5deHAuU+DMdfNE7fGQgxgsbpawMhzkDHONoYk9HVBwM+MZwMEVlrjMkPdhy90Th9ZyDECBqnrw2EOEMpRq1KUkop1YEmBqWUUh0M9sTwRLAD6CON03cGQoygcfraQIgzZGIc1G0MSimlTjTYSwxKKaU60cSglFKqg0GbGERkkYjsFJE9InJXsOPxJCIHRGSLiGwUkbX2sjQReUdEdtvPQwIc0zMiUioiWz2WdRuTiPzYPrY7ReQzQY7z5yJSZB/PjSJycTDjFJGRIvK+iBSIyDYRudNeHlLHs4c4Q+14xorIGhHZZMd5r708ZI5nDzGG1LE8xhgz6B5Yw37vBcYB0cAmYHKw4/KI7wAwtNOy3wJ32a/vAh4IcExnAbOArb3FBEy2j2kMMNY+1pFBjPPnwPe7WDcocQLZwCz7dRKwy44lpI5nD3GG2vEUINF+HQV8AswLpePZQ4whdSzbH4O1xDAH2GOM2WeMaQUWA5cHOabeXA781X79V+Bzgdy5MWY5UNnHmC4HFhtjWowx+7Hm15gTxDi7E5Q4jTHFxpj19us6oABrbvOQOp49xNmdYMVpjDH19tso+2EIoePZQ4zdCdr/EAzeqqQRwGGP94X0/AcfaAZ4W0TWicit9rIsY89mZz9nBi2647qLKRSP7zdFZLNd1dRepRD0OEVkDHAq1hVkyB7PTnFCiB1PEYkUkY1AKfCOMSbkjmc3MUKIHUsYvIlBulgWSv12zzDGzAIuAu4QkbOCHVA/hdrx/TMwHpgJFAO/t5cHNU4RSQT+DXzHGFPb06pdLAtmnCF3PI0xLmPMTKx54ueIyNQeVg9KnN3EGHLHEgZvYigERnq8zwGOBCmWExhjjtjPpcArWEXIEhHJBrCfS4MX4THdxRRSx9cYU2L/U7qBJzleJA9anCIShXWy/bsx5mV7ccgdz67iDMXj2c4YUw18ACwiBI9n5xhD9VgO1sTwKZArImNFJBq4FlgS5JgAEJEEEUlqfw1cCGzFiu8r9mpfAf4TnAg76C6mJcC1IhIjImOBXGBNEOIDjp0U2l2BdTwhSHGKiABPAwXGmD94fBRSx7O7OEPweGaISKr9Og44H9hBCB3P7mIMtWN5TKBauUPtAVyM1ctiL/CTYMfjEdc4rN4Im4Bt7bEB6cB7wG77OS3Acb2IVdRtw7qaubmnmICf2Md2J3BRkON8HtgCbMb6h8sOZpzAAqxqgc3ARvtxcagdzx7iDLXjOR3YYMezFfiZvTxkjmcPMYbUsWx/6JAYSimlOhisVUlKKaW6oYlBKaVUB5oYlFJKdaCJQSmlVAeaGJRSSnWgiUEpDyKSKiLf8Hg/XERe8tO+PiciP+vms3r7OUNE3vTH/pXqjiYGpTpKBY4lBmPMEWPM1X7a1w+BP/W0gjGmDCgWkTP8FINSJ9DEoFRH9wPj7bHxfyciY8Se20FEbhSRV0XkvyKyX0S+KSLfFZENIrJaRNLs9caLyJv2IIgrRGRS552IyESgxRhTbr8fKyKrRORTEfllp9VfBa7360+tlAdNDEp1dBew1xgz0xjzgy4+nwp8EWtMm18DjcaYU4FVwJftdZ4AvmWMmQ18n65LBWcA6z3ePwz82RhzGnC007prgTO9/HmU6jdHsANQaoB531hzE9SJSA3wX3v5FmC6PRLp6cC/rKGGAGuylc6ygTKP92cAV9mvnwce8PisFBjum/CV6p0mBqX6p8XjtdvjvRvr/ykCqDbW8Mo9aQJSOi3rbnyaWHt9pQJCq5KU6qgOaxpLrxhrvoL9InINWCOUisiMLlYtACZ4vF+JNcovnNieMJHjo24q5XeaGJTyYIypAFaKyFYR+Z2Xm7keuFlE2kfI7Wra2OXAqXK8vulOrEmZPuXEksQ5wOtexqJUv+noqkoFiYg8DPzXGPNuL+stBy43xlQFJjI12GmJQang+Q0Q39MKIpIB/EGTggokLTEopZTqQEsMSimlOtDEoJRSqgNNDEoppTrQxKCUUqoDTQxKKaU6+P/UHkAdbKu2qQAAAABJRU5ErkJggg==\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": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAElCAYAAADnZln1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA2mElEQVR4nO3deZxb5XXw8d+RRrN7Fs+M7bHH+4oB2xgHr2HHGEKwgUBwWULiJG3TNEmbNE3SvinN27xJ2rdp836atCUpIQTCTtgXG4MxYBt7DGO87/t4GXtWzz7Sef+4d4wYj2Y0skbSSOf7+ciWdK90j66kM4+e+9zziKpijDEm+XniHYAxxpjYsIRvjDEpwhK+McakCEv4xhiTIizhG2NMirCEb4wxKcISfooRkftF5BH3+igROSMi3njH1RMR+bSI7Ix3HNB7LLHcpyKySkS+7F6/S0SWBy2bLyK73ViWiMhQEVktIg0i8q/9HZtJTJbwBxgROSAi13a57z4Rebevz6Wqh1Q1V1X90Yuwb0RERWRCT+uo6juqOjlWMfWkayxd34947VNVfVRVFwbd9SPgP9xYngO+CpwC8lT127GMzSQOS/gmoYlIWrxjGKBGA1u73N6mEZxpae9B8rCEn4REZLiIPCMiVSKyX0S+EWK9MW4LOy3ocS+ISLWI7BGRrwSt6xWRH4jIXrdbYKOIjHSXTRGRFe7jdorIHUGPe0hEfikiL7uPe19ExrvLVrurbXK7Hj4vIleKyBER+VsROQ78tvO+oOccKSLPuq/vtIj8R4jXd7+IPC0iT7jb/kBEpgctv8DtFqkVka0icnPQshtFZJv7uKMi8h33/rOxiMjvgVHAi2783+3jPr1fRJ4UkYfd7WwVkVk9vK/XicgOEalzX7MELTv7K09E9gLjguJ6DPgC8F339rUi4hGR77nv52k3jsFdPhfLROQQ8KZ7/5dEZLuI1IjI6yIyOmj7KiJ/5nYj1bjveXB8X3Ef2+Du15lB+6fbz6qIXCYi5SJSLyInROTnofaNCZOq2mUAXYADwLVd7rsPeNe97gE2Aj8E0nG++PuA693l9wOPuNfHAAqkubffBn4FZAIzgCrgGnfZ3wCbgck4iWY6UATkAIeBLwJpwEycroML3cc9BFQDl7nLHwUeD4pdgQlBt68EOoCfARlAlnvfEXe5F9gE/Ju77UxgQYh9dT/QDnwO8AHfAfa7133AHuAH7n66GmgAJruPPQZ82r1eCMwMiu9IqPejj/v0fqAFuNF9XT8B1oV4LcVAfdBr+St3P32562cgRFwPAf8UdPtbwDqgzN3P/w081uU1POzu4yxgibu/LnDfx78H1nR5H18CCnD+CFYBi9xltwNHgU/hfHYm4Pzi6O2zuha4x72eC8yJ9/dvoF/iHoBd+viGOV/kM0Bt0KWJjxP+bOBQl8d8H/ite/1+ukn4wEjADwwKetxPgIfc6zuBxd3E83ngnS73/TfwD+71h4DfBC27EdgRdLu7hN8GZHa5rzPhz3WTSVoY++p+ghKom2COAZ92L8cBT9Dyx4D73euHgD/F6fOmu1iC3o9uE34Y+/R+4I2gZVOB5hCv5d4ur0WAI0Se8Lfj/uFxb5fi/HFMC3oN44KWvwos67Ivm4DRQe/jgqDlTwLfc6+/Dnyzm9fU22d1NfCPQHG8v3fJcrEunYFpiaoWdF6ArwUtGw0Md7spakWkFqcVO7SX5xwOVKtqQ9B9B4ER7vWRwN5uHjcamN1le3cBw4LWOR50vQmntdaTKlVtCbFsJHBQVTt6eY5OhzuvqGoAJ0kOdy+H3fs6Bb/e23D+OB0UkbdFZG6Y2wvW2z6Fc/dNpnTfZz68y2vR4NsRGA38Meg9247zxyn4c3K4y/q/CFq/GuePTk+vpfN97umz09NndRkwCdghIhtE5KY+v0rzCXYwJvkcBvar6sQ+Pq4SGCwig4IS1Cicn+Kdzzse2NLN9t5W1esiDbgbPR1YPAyMEpG0MJP+yM4rIuLB6cKo7FwmIp6gpD8K2AWgqhuAxSLiA76O02I9+1xhxtrbPu2LY11ei4SIJ1yHgS+p6ntdF4jIGPeqdln/x6r6aITbGh/i/pCfVVXdDSx137dbgadFpEhVGyOIwWAHbZPReqDePeiZJc7B1otE5FM9PUhVDwNrgJ+ISKaITMNpYXV+wX8D/G8RmSiOaSJShNNvO0lE7hERn3v5lIhcEGa8J3D6bvvy+o4BPxWRHDfW+T2sf6mI3Oq2mr8FtOL0Xb8PNOIcyPSJyJXAZ4HHRSRdnHHt+arajtN3HmqYZcj4w9inffEycGHQa/kGn/wV1Vf/Bfy488CriJSIyOJe1v++iFzorp8vIreHua3fAN8RkUvdz84Ed7s9flZF5G4RKXH/INe6zxW3IcTJwBJ+klFn/PdncQ4Q7sc5gPobID+Mhy/F6b+tBP6I0w+/wl32c5xW7nKcBPg/QJbbcl0I3Ok+7jgfH3ANx/3A79yf9Hf0tnLQ65uA089+BOc4QijPu8trgHuAW1W1XVXbgJuBG3D20a+Ae1V1h/u4e4ADIlIP/Blwd4jn/wnw92783+lmeU/7NGyqegrn4OdPgdPAROCc1nkf/AJ4AVguIg04fwRn97D9P+K8r4+7+2QLzr4LJ/angB8Df8A5MP4cMDiMz+oiYKuInHHjvbOHrj4TBnEPjhiTdETkfpwDwqGStTEpxVr4xhiTIizhG2NMirAuHWOMSRHWwjfGmBRhCd+YCEiXcsQ9rHe2HHUiEKe20T/FOw4TH5bwTb+Tj2vEd15URBqDbn86guc8p0x0l+VXikjAff4GcYq6fTHC+D9REA26LUdsTMKzM21Nv1PVQwSVUxARBaar6p5+3nSlqpa5Z6UuxjlT831V3RbuE4Qoc2DMgGQtfBNXIpIhIv9XRA6JUwL3v0Qky11WLCIvuSc1VYvIO+KU9T2nLHFP21DHczgnX00Vkc+IyIfilN097I7X74ynu9LAnWWca93tzZUuk86IyIXycYnoEyLygxCvd46IrHFf0yb3DN/OZfeJyD73F8l+Ebmrh3327yJS6V7+XUQy3GWd5aW/LSInReRYqF82IrJFRD4bdNsnIqdEZEZP+9MMXJbwTbz9DKdA1gycs2dH4JTLBfg2zpm0JTgFtX6Ak7/vwTnL9rPqzOj0zz1twP0jcQtO6d7NOCUV7nVvfwb4cxFZ0uVhV+CUAr4euNy9r8Dd3touzz8IeAN4DafI2QRgZTdxjMApkfBPwGCccs3PuGUNcoD/B9ygqoOAeUBFiJf0d8AcnH02Haf09N8HLR+Gc7bqCJxSDr8UkcJunudhPnkG8Y3AMVUNtV0zwCV8wheRB92WSteiXZE+32tu6+qlLvd/XZwJKlREiqOxLdMzt6vlK8BfqWpnVcn/g1OmAZxyvaU4JXjb1ZlesC/jiIeLU4HxFPAPOLXVd6rqKlXdrKoBVf0IpyzyFV0ee7+qNqpqcxjbuQk4rqr/qqotqtqgqu93s97dwCuq+oq77RVAOU6iBQgAF4lIlqoeU9Wt3TwHONVIf6SqJ1W1CqeE8D1By9vd5e2q+gpOOe3upoh8BLhRRPLc2/cAvw/j9ZoBKuETPk4d70VRfL5/4ZNfjk7vAdfilK81sVECZAMb5ePyuK+594PzXu3BqfeyT0S+18fnr3RLSA9W1Rmq+jiAiMwWkbfEmWWpDqdWTtc/8n0pPRyq/G9Xo4Hb5ZPlgBcApW4FyM+7sRwTZ4awKSGeZzif/JwedO/rdLpLJdFuS1KraiXO5/42ESnAqY0TSWE3M0AkfMJX1dU4tbfPEpHxbkt9o9uvG+qL0d3zrcQp4NT1/g9V9cB5B2z64hTQjDM7Vmd9/3xVzQVwW8rfVtVxOEW2/lpErnEfez5nDP4Bp3DYSFXNx6kEKV3W0RDXuxOq/G936/0+eC4DVc1R1Z8CqOrrbpnpUmAH8OsQz1OJ88ej0yg+LvncV7/D+eVxO7BWVSMp3WwGiIRP+CE8APylql6K0w/6qzjHYyLglr39NfBvIjIEnH5uEbnevX6TOKV0hY9LFHeWx+1rWeVgg3AmJmkRkcuAP+ll/Sqc7pZQ23sJGCYi33IPqA4Ske4qTz4CfFZErhenFHCme5C1TESGisjNbl9+K043TKhSwI/hVOgscbsff+g+dySew5mW8ps4ffomiQ24hC8iuTgHtJ4SkQqc6fRK3WW3uiMPul5ej2PIpmd/i9Nts06csrtv8HF/80T39hmc+U1/paqr3GW9lSXuydeAH4lTFviHOGWfQ1LVJpzyvu+525vTZXkDcB3Or5DjwG7gqm6e5zDO8NAf4PwROYwzV7DHvXwbp6VejXNM4Wtdn8P1Tzh9/x/hHIT+wL2vz9xjFM8AY4FnI3kOM3AMiFo64szA85KqXuQeYNqpqqXn8XxXAt9R1XOmTBORA8Ast/64MUlPRH4ITLIy0slvwLXwVbUe2C/ubDvimB7nsIwZkERkMM7QzQfiHYvpfwmf8EXkMZyf85PdE0qW4QxLWyYim4CtOD+Tw32+d4CngGvc5+vsL/6GiBzBmfP0IxH5TbRfizGJRES+gtOt9Ko7OMIkuQHRpWOMMeb8JXwL3xhjTHQkdGGo4uJiHTNmTLzDMMaYAWPjxo2nVLWku2UJnfDHjBlDeXl5vMMwxpgBQ0RCVguwLh1jjEkRlvCNMSZFWMI3xpgUYQnfGGNShCV8Y4xJEZbwjTEmRVjCN8aYFGEJ3xhjEsjhw4d57733+uW5LeEbY0yC2L9/Pw8//DAbN26ktbU16s+f0GfaGmNMqti1axdPPvkkhYWF3HvvvWRkZER9G5bwjTEmzrZt28bTTz/N0KFDufvuu8nJyemX7VjCN8aYONq0aRPPPfccZWVl3HXXXWRmZvbbtizhG2NMnGzYsIGXX36ZsWPHsnTpUtLT0/t1e5bwjTEmDtasWcPy5cuZNGkSt99+Oz6fr9+3aQnfGGNiSFV5++23WbVqFRdeeCG33norXq83Jtu2hG+MMTGiqqxYsYI1a9YwY8YMbr75Zjye2I2Ot4RvjDExEAgEeOWVVygvL+eyyy5j0aJFMU32YAnfGGP6XUdHB88++yzbtm1jwYIFXHPNNYhIzOOwhG+MMf2otbWVJ554gn379rFw4ULmzZsXt1himvBF5K+ALwMKbAa+qKotsYzBGGNipbGxkUcffZRjx46xZMkSZsyYEdbjNKCIJ/q/AGLWgSQiI4BvALNU9SLAC9wZq+0bY0ws1dbW8uCDD3Ly5EnuvPPOsJK9qlL/1mFO/XYL6teoxxTrLp00IEtE2oFsoDLG2zfGmH538uRJfv/739PW1sY999zD6NGje31MoM1PzTO7ad5URdb0EggoeKPbyo9ZwlfVoyLyf4FDQDOwXFWXx2r7xhgTC4cPH+YPf/gDXq+XL37xiwwbNqzXx3TUtnD64W20H2skb9EYBl1R1i8HdWPZpVMILAbGAsOBHBG5u5v1vioi5SJSXlVVFavwjDHmvO3evZuHH36YzMxMli1bFlaybz1Qx8n/qKDjdAtFX7iQvCtH9tsInlgOAr0W2K+qVaraDjwLnHO4WlUfUNVZqjqrpKQkhuEZY0zkNm/ezGOPPUZRURHLli2jsLCw18ecWX+Mql9vxpOZxpC/mEHWlMH9GmMs+/APAXNEJBunS+caoDyG2zfGmKhTVdasWcOKFSsYPXo0S5cu7bXipfoD1L60j8a1x8iYVEjR0il4svo/HceyD/99EXka+ADoAD4EHojV9o0xJtoCgQCvvvoqGzZs4MILL2TJkiW9FkHzn2mj+g87aN1XR+7lZeQvGtMvQzC7E9NROqr6D8A/xHKbxhjTH9ra2njmmWfYuXMn8+bN49prr+21VELroXqqH92Ov7GDws9PJueSITGK1mFn2hpjTB+dOXOGxx57jKNHj3LDDTcwe/bsHtdXVRrfP07ti3vx5mcw5M+nkz4iN0bRfswSvjHG9MHp06d55JFHaGho4M4772TKlCk9rq/tfmqe20vTxhNkTi5k8Ocn48nu/9r33bGEb4wxYeocYy8i3HfffZSVlfW4fkd1C6cf2UZ7ZSODrhlF3jWjYtZf3x1L+MYYE4bt27fzzDPPkJeXx1133UVRUVGP67fsrOb04ztBoegLU8m6oOf1Y8ESvjHG9EBVWbt2LcuXL6esrIylS5eSk5MTev2A0vDmIepXHsI3NIeiey4grSgrhhGHZgnfGGNC6Ojo4OWXX+bDDz9k6tSpLFmypMeJxv0NbVQ/voPWvXVkXzKEglsm4EmPzfSF4bCEb4wx3WhqauLJJ5/kwIEDXH755Vx55ZU9Drts2VND9eM70VY/hbdNJHvW0LhMctITS/jGGNNFVVUVjz32GHV1ddx6661MmzYt5LoaUOpXHqLhzUOklWRR9JWL8Q0N3eUTT5bwjTEmyN69e3nyySfxer184QtfYNSoUSHX9de7XTj76sieOYSCJYnVhdOVJXxjjHFt2LCBV155hZKSEpYuXdpjAbSW3TVUP+F24XxuEjmzhsYw0shYwjfGpDy/38/rr7/O+vXrmThxIrfddlvIAmjqV+pXHqThrcOklWQndBdOV5bwjTEp7cyZMzz11FMcPHiQOXPmsHDhwpAHZztON1P9xE7aDjWQfelQChaPT+gunK4s4RtjUtbRo0d54oknaGpq4pZbbmH69OndrqeqNG08Se0Le8EDg5dOIXv6wJuvwxK+MSYlVVRU8OKLL5Kbm8uXvvQlhg8f3u16gaZ2av64h+bNp0gfm8fgz08mraDneveJyhK+MSalBPfXjxkzhttvvz3kmbMte2upeXIn/oZ2Z67Zy8viWgvnfFnCN8akjK799ddddx1e77l98NoRoH7FQRpWHyGtKIshX5tKetmgOEQcXZbwjTEp4ciRIzz55JM0NTX1eDJVW+UZap7aRfuxRnIuG0b+TeMG1IHZnljCN8YkNVVl3bp1rFixgry8PJYtW0Zpaem56/kDNKw6Qv3KQ3hy0ii6dypZU+Nf4TKaLOEbY5JWc3Mzzz//PDt27GDy5MksWbKErKxzK1e2H2+k+qldtB89Q9aMEgo+Ox5vTnwmKelPlvCNMUmpsrKSp556irq6OhYuXMjcuXPPKWamfqVh9RHq3ziIJzONorsvIOui4jhF3P8s4SeBQJufgz9bQ9rwbIbffjHevNDlW41JdqpKeXk5r732Gjk5Odx3333d1sNpO9ZIzTO7aD9yhqyLiylYPB5vbnJ/dyzhJ4Ez1fX4GoHdTRz52ToGXTeKwk+PRrwDd/iYMZFoaWnhpZdeYsuWLUyYMIFbbrnlnCGXgTY/DSsP0fDOUTxZaQP2JKpIWMJPAhpQAPbmnSLjjOB9Tah//yild15Mxui8OEdnTGwcPHiQZ599lvr6eq655hrmz59/TomElt011PxxD/7qFrJnDSX/hrFJ2VcfiiX8JKDqJPziCaUM/vRoyp94m0nHizj5nxV4pxUybMkUPNmp86E2qcXv97Nq1SreffddCgoKWLZs2TmTi/vPtFH38n6aPjxJWnEWxV+5mMzxBfEJOI4s4ScDN+EDlJaW8plv3EHF+g84+OoepnykHNq+lsE3jWfQp4YP6LMEjenq9OnTPPPMM1RWVnLJJZewaNEiMjIyzi7XgNL0wQnqXtlPoNXPoKtHknfVKMQXeuaqZGYJPwl0tvA7RyB4PB5mzpnFmYum8O7zb1KyRUj74z5q1hxm+NKL8Q0bGKVcjQlFVfnwww959dVX8Xq93HHHHUydOvUT67QdbqD2hb20HW4gfXQehbdOGDBljPuLJfwk0NmHT5fGe25uLovuupn9+/bxwZNrmXpiGMf/fSMZs4dQfMMEPJn29puBp6GhgZdeeomdO3cyduxYlixZQn5+/tnl/jNt1L9+kMby43hyfBTeMYnsS4Yk3Pyy8WDf+GTQ2aMT4vM8dtw4Rn5nFOtWraF11TEmvy8crqhi8I3jyf1UqXXzmAFBVdm8eTOvvvoqbW1tLFy4kDlz5pw9MKt+pXFdJXUrDqFtfnIXjCDvmlHWsAlieyIJdHbp0EMLJi0tjQXXXk71zGreefoNyvZn4f3jXmrfOcSQ2y4gY2x+yMcaE28NDQ28/PLL7Nixg7KyMhYvXkxJiTOUUlVp2VlD3Sv76TjZRMaEAgpuHo9vSHaco048lvCTQGeXTjg/WQcPHszNX7mdnTt28P5z67nw1HCq/vsj0i7Ip/jmSaQVDsw63yY5qSofffQRr732Gm1tbVx33XXMnTv3bKu+7UgDda/sp3VfHWnFWRTdfQGZFxZZ900IlvCTgAaN0gmHiDDlggsYP2EC7616l4bVB7h4+0gqd25g0OUjyL96dNJUBzQD1+nTp3n55ZfZt2/fOa36juoW6pYfoLmiCk9OGgWLx5Nz2TDEm5qjb8JlCT8pdH/Qtjc+n48rr7uK6kureeuFFRTt8jBhFTSsP0bRTRPInjHE+vdNzHV0dLBmzRrefvtt0tLSuPHGG5k1axYejwd/fRsNbx/mzLpjIMKgq0Yy6Ioy66cPk+2lJNCXLp3uDB48mNvu+zy7du3i7RfeYerpYXie3EXt6kMU3zqZjFF2tq6JjUOHDvHiiy9SVVXF1KlTWbRoEXl5efgb2qh9+4iT6AMBsmcOJe+60aTlZ/T+pOYsS/hJINSwzL6aNGkS4741jvING1izYhszjo9Gf7UJ39QCij4zgbSic8vKGhMNDQ0NrFy5koqKCvLz81m6dCmTJ0/Gf6aN2lf207i2Eu0IkH3JEPKuGWWfxQhZwk8mUThQlZaWxpy5c5k2fTqr31hF2/pKpm0bxbHtteTMGUb+NaOTvqKgiZ2Ojg7WrVvH6tWr6ejoYN68eVxxxRV4m6H2pX00rj+GtgfInjGEQVePxFdiI2/OhyX8JHC+XTrdyc7OZtHNN3JyzklWvbySwj0wea3SWH6C/KtGkbtghB3YNRFTVXbs2MHy5cupqalh8uTJLFy4kLyOTBqeP0BTRRWgZE8rYdDVo2yIZZRYwk8GGp0une4MGTKEO764lN27d/PmS6uYUFXEmOVK3btHKFw0lpxLh9rICNMnlZWVrFixgv3791NSUsI9d9/NCCnmzEtHObG9GvF5yJ1TSu6CEaQNtmHC0RTThC8iBcBvgItwhpZ8SVXXxjKGpHT2TNv+G1EzceJExn1jHB9++CEr39jIRXXDkWf91L5xgMKF45xT163+vulBVVUVb731Ftu2bSMrK4vPLLyByZTR9OJxTp04hic7jUHXjCJ33vCUKlkcS7Fu4f8CeE1VPyci6YD9TouCj4un9e92vF4vs2bNYvr06ax//33eWrWVixrK4OkOalfup3DhOLKml9hQTvMJtbW1rFq1ik2bNuHz+bhu5hVM9o+g9fVT1Lfuwzcil8LPTSJ7eknKVrGMlZglfBHJAy4H7gNQ1TagLVbbT2baj1063fH5fMxfsIBLZ81i7Zq1bHlnK9NqR6FPtCNvHKBw4ViyLi62xJ/iqquree+996ioqCBd07ihbD4jGwvxr2mixXuS7IuLyZk3nPSRg+zM2BiJZQt/HFAF/FZEpgMbgW+qamPwSiLyVeCrQLfzUJpuaPQP2oYjMzOTq66+isbZl/HuO++yae1WLqkegz7WiueNTAqvH0vm1CJL/Cnm5MmTvPvuu2z9aAsjtZjFg+YxuCYDdiueYULuTePInlFio73iQHo7LV9Ews26tapa38PzzALWAfNV9X0R+QVQr6r/K9RjZs2apeXl5WFuPnUdrdiPPn6E5qvzmLhwetziqKur4+1Vb1NXfpSZ/nHkB7LxDsuiYOFYMi8YbK24JKaqHDp0iHXvraVx+ynG6zDG6BC8fsGTnUbW9BJyZg3DNzzHPgf9TEQ2quqs7paF08L/Hc5hwZ7eJQUeAh7uYZ0jwBFVfd+9/TTwvTC2b3oRrROvzld+fj43L76ZU/NPseqtVbR8tJ9LT47D/3Az3mHZ5F89iqyLrKsnmbS3t7N1/SaOrNlNXo2PTwWGks4IJNNL1kXFZF9cTMaEAhvJlSB6TfiqelXX+0RkmKoe78uGVPW4iBwWkcmquhO4BtjWl+cwIcSpSyeU4uJiPnf75zi24BhvvfkWur2BGSfH4P9DE57BGeRdOZKcmUORNEsCA1Ggzc+pzUeo3LAPOdxCsT+XYsrwZwo5U4vJmT6EzPEF9v4moEj78O8F/jmCx/0l8Kg7Qmcf8MUIt2+ChVEPPx5KS0v5k7v+hMrKSt59512atuzlktoxBJ5tpW75AfIuH0nO7GF4Mux0kESmAaVlVw1NO09Rv7MKb7UfD8JgPDTk+PBfmE/pnLH4SnMTptFhuhfpN22xiDQBK9zWelhUtQLotm/JRE5VERKnhd/V8OHDuePzd1B1dRXvvvMu6yv2ML1pDPpKB3VvHmLQvOHkzh9hY68TVN3m45x5bA9+ApySOhoGtZE/dRgTFkxjVElhvMMzfRBpwr8VuAS4RUQmqOqXoxiT6atepjhMFCUlJdxy6y3UXFXDe++9x8byjUxrG8XoN/00rD5Czsyh5M4fnvITTSeakwcqyQa2zahj2qcvZfjw4QnbuDA9iyjhq+oJ4DX3YuIs0Vv4XRUWFnLTTTfRcEUDa9eu5bl15UxpLWXihgCN64+TMbGA3AUjyJxYaAd4E0BHUzsA865aQP7QwXGOxpyPiBK+iPwSyFHV+0Rkoaouj3Jcpg9ifeJVtAwaNIiFCxeyYMECysvLeWHdB4ysL+CivaNo3V2LtziTQfNHkD1zKJ4MK9QWL4GWDgAyBllJ4oEu0i6dNuCEe/1qwBJ+PMWglk5/ys7O5vLLL2fevHls3bqVlWvWkV0ZYFrNGPzPt1D76n5yLhlCzuxS0ofnxjvclBNo7aADyMiyQmYDXaQJvwnIFxEfYKfDxpkm2LDMSKWlpTF9+nSmTZvGgQMHWLtmLbU7dzDVX8a49X4a3z+OryyX3NmlZE0vsfLMMaKtAdpFB/zny0Se8KuBZuCXwHvRC8dEJBCb4mmxIiKMHTuWsWPHcvr0acrLy3nmgw2UNRZw4fGRtD9zhpqX9pJzyVByZg3FN8KGA/arNj8dNqQ+KfQp4bvljf8NmAw8gnNm7bLoh2UikoRJr6ioiOuvv56rr76abdu2sXZDOW2H6rmgo4xx7/tpXHeMtCFZZM8cSvYlQ2yO034g7YrfEn5S6FPCV9VaEfkpMAY4BUwDnu2HuEwfDNSDtn3h8/mYPn0606dP58SJE1RUVPBsRTlDG3KZXD2CIa81U//aATImFJA9cwhZFxbbgd4okQ4IpPVcc8sMDJF06SwD9qvq6zgVL02c9ccUh4ls6NChXH/99Vx77bXs2bOHiooK3t6+lfHtQ5l8oIzWPbXU+PaQdWER2RcXkzmpEPFZ8o+Up0Noz7CEnwwiSfg1wJ+JyGRgE1Chqh9GNyzTJwN8lE6kvF4vkydPZvLkyTQ1NbF582ZWb/qIjiNnmOgvZdxHbTRXVIHPQ9YFg8m6qJjMyYOt5d9HXr/Q5kutz1ay6nPCV9WfiMhKYBcwA2dSE0v4CSDF8v0nZGdnM3v2bGbPnk1NTQ3btm3jza3bCBxuYmxgCOO2tNH80SlIEzInFpI5eTCZUwpJK7Chhr1JC3jAStcnhT4nfBH5EeAFKnBa96uiHJPpo7NdOnZWKuCcyTt//nzmz59PbW2tk/y3bMV/pJGx/iGM2dVKzvZqANKGZJE5pYjMyYVkjMmzMr7d8KkXsYSfFCJp4f9QRIbi1NK5TUTGq+pXoh+aCdfHk9hYwu+qoKCAefPmMW/ePGpra9m+fTvrdu6i9kAVIzoKGXW6hGHvNHFm9RFI95A5roCMcflkjMt3qj+m+MTs7S1tePEgNvgpKUQ6Dv9Pgf9WVaulk0Cshd+zgoIC5s6dy9y5c2ltbeXAgQPs3r2b9bsqyKnxUOYvpmxPE4N2OK1/0j1kjnWSf/rIPHxluSl3sldLQzMAnkwrYZ0MIn0XHwT+XERygEfdsscmXhJkxquBJCMj4+wBX1Xl1KlT7Nmzh41793Lq4HEGt+RQ6i9gxJ4i8nbWOA8SSBuaTcaoPNJHDsJXNghfSVbSTfTR0drOyR1Hqd51jLYD9Qwhm7Qs69NJBpEm/G/g1NNJA/4fzoFbEyeaoBOgDBQiQklJCSUlJcydO5dAIMCJEyc4dOgQWw4e5MSBbWQ3eCkJ5DO0Kp8hJ/PxrXdb+h7wFmWRXpqDb2gOvmHZpA3JJq0wM+H/EGhHgKYT9VTvO0HDoWr8J5pJqw2Q0+bDg4cCoJ0ManNbKJt1QbzDNVEQacLfC0wEnlfVv4piPOY8pMo4/P7m8XgoLS2ltLSU2bNno6pUV1dz9OhRKisr2Va5nzNHa8lvy2RwIIfC07kUVeeR+9HHHd0KePJ8+IqySCvKIq0oE29BJt5cH95B6XhyfXiyff3WDaeqaKsff30bbTXNnDlRS1NVPa2nm9CaNtIaIaPNiwdBgDygUQI0ZrXTPEzIHlVAyZQRDB87BI8dyE4akSb8rcBhYJmI/IuqfiqKMZm+slE6/UpEKCoqoqioiGnTpgEQCASorq7mxIkTnDp1iq2nTlF94jT+qmZyOzLIC2Qx6EwWeQ1Z5B3IJkvP7RJRAbI8SKYXyfDiyUjDk+n8783yOQeMPZ0XQAT1B/C3ddDR1kGg3f/xpbmDQEsH2uLH06542wWvfjJR+4AA7TR4WqjL9MPgNHxDshlUNpjiCcMYPqwIj8eSezKLNOFPAqqAB3BOxDJx9PEoHRMrHo+H4uJiiouLP3F/IBCgvr6empoaamtrqaur40hdHQ3VdbTVNKONHaS1CVmaTramk9WaTnqLDx9efJqGDy/p7v9O61vwIHj4OBH7CeAnQIAAfhS/BGijgzbpoI12Aj5gkAfN9kJOGumFWeQMySe/tJAhw4oZl2vF5lJVpAl/Ck4N/F8CB3H69E2c2Zc4/jweDwUFBRQUFIRcx+/309zcTHNzM01NTTQ3N9Pe3k57ezstHR2c6eigvd2ZZUpEnPdVnWPy4hHSM7JJT0/H5/ORnp5OVno6GRkZ5ObmkpWVZa10E1KkCb8A+Fvgu1i1zLg7e8zWunQGBK/XS25uLrm5NpmLia1IE/6PgCmqulNEAtEMyETARukYY8IQ1m8/EfGKyDER+TKAqh5R1Tfc69/rzwBNGCzfG2PCEFbCV1U/sAUY37/hmEjYOHxjTDj60qWTDXxXRK4DKt37VFUXRz8s0ydqwzKNMb3rS8Kf6/4/073Ax5XYTTzZu2CMCUNfEv7YfovCnBc928K34XjGmNDCTviqerA/AzHn4WwXvnXpGGNCsyZhMrBy+MaYMFjCTwJqB22NMWHoc8IXkc/2RyDmPHQmfOvSMcb0IJIW/o+jHoWJDsv3xpgeRJLwLa0kGjtoa4wJQyQJ30Z9JxgblmmMCYdliGRgo3SMMWGwhJ8M7KCtMSYMkST8E1GPwpwXtT58Y0wY+pzwVfW689mgW2r5QxF56XyexwQ5m/HjG4YxJrHFo0vnm8D2OGw36dlBW2NMT2KaIUSkDPgM8JtYbjfpWR++MSYMESV8EfnroOuT+/DQf8eZBzfktIgi8lURKReR8qqqqkjCSzlqo3SMMWHoU8IXkQIR+S1wu4h8TUQWAGFNcSgiNwEnVXVjT+up6gOqOktVZ5WUlPQlvJRntXSMMT3p0yTmqloLfFFEPgMcBxYCz4b58PnAzSJyI5AJ5InII6p6d19iMN0IdHbpWB++MSa0SDPEFTjDM+cAYY3aUdXvq2qZqo4B7gTetGQfXdbCN8b0JNKEXwD8LU5/fEvUojGRsVGZxpgw9KlLJ8iPgCmqulNEQh6ADUVVVwGrIty26cLq4RtjwhFRwlfVI8AR93pYB21NP+ps4ds4fGNMDyIdlvlLEXnIvb4wqhEZY4zpF5E2CduAfe71q6MUi4mUdekYY8IQacJvAvJFxAeMimI8JgLWh2+MCUekB22rgWbgl8B70QvHRELOTkljCd8YE1qkZ9re5t71MDAr6lGZPjlbWcFa+MaYHvT5TFsR+SkwBjgFTCP8M21Nfzk7SscSvjEmtEi6dJYB+1X1daDHujgmRtyE77FhmcaYHkSS8GuAP3OrZG4CKlT1w+iGZfpGUZtb3hjTiz4nfFX9iYisBHYBM4DLAUv48aRYujfG9KrPCV9EfgR4gQqc1v2qKMdk+krP/mOMMSFFMqftD4FW97G3icivox6V6SPr0DHG9C7So3wPAhcARcCvoheOiYRatjfGhCHShP8NnO6gNOAX0QvHRELU2vjGmN5FmvD34sxa9byqXh7FeEwELNkbY8IRacLfCrwJLBORDVGMx0TCRukYY8IQaS2d8Tjj8R9w/zfxZKN0jDFhiDThH1bVN0WkFDgZzYBMZCzdG2N6E2mXziIRKQP+C/i3KMZjImHDdIwxYYjGJOatUYvGREZBrW6aMaYX0ZjE3B/NgEykrJVvjOlZWC18EfGKyDER+TI4k5ir6hvudZvEPN5slI4xJgxhJXxV9QNbcEbnmERjffjGmDD0pUsnG/iuiFwHVLr3qaoujn5Yxhhjoq0vCX+u+/9M9wLWk5AYrLSCMSYMfUn4Y/stCmOMMf2u14QvIqPcq902IYOW16pqfbQCM32ggA3LNMb0IpwW/u/oPaUo8BDwcBRiMn1lo3SMMWHoNeGr6lWxCMScD2viG2N6F+mZtibB2EFbY0xvLOEnA8v1xpgwWMJPBtajY4wJgyX8JGGNfGNMbyzhJwObAMUYEwZL+EnBkr0xpneW8JOBiqV8Y0yvLOEnBbWDtsaYXsUs4YvISBF5S0S2i8hWEflmrLad9Kx5b4wJQ6QzXkWiA/i2qn4gIoOAjSKyQlW3xTCGpCRYzjfG9C5mLXxVPaaqH7jXG4DtwIhYbT+Z2fwnxphwxKUPX0TGAJcA73ez7KsiUi4i5VVVVTGPbcCyPnxjTC9invBFJBd4BvhWd+WUVfUBVZ2lqrNKSkpiHd6AJDYBijEmDDFN+CLiw0n2j6rqs7HctjHGpLpYjtIR4H+A7ar681htNyVYLR1jTBhi2cKfD9wDXC0iFe7lxhhu3xhjUlrMhmWq6rtYO7TfWA++MaY3dqZtMrAuHWNMGCzhG2NMirCEnwTE+nOMMWGwhJ8k1Lp0jDG9sISfDKy2gjEmDJbwk4I1740xvbOEnyws5xuTNLSffrVbwk8G1qNjTFJob2/n1Vd/zYMP/qhfkn4s6+GbfmL18I0Z2FSVioo3WbFiJU1N6QwdWktrawOZmXlR3Y4l/GRhXTrGDEhHjuzkxRf/wIkTQk5OCzfcOJbLPvW/EPFGfVuW8JOBnWlrzIBz6tRxXn31IfbubSYtrY1LL83guut+QGZmcb9t0xJ+klA7+8qYAaGhoY7XX/89W7eeRCTAxIktXH/9lyguntrv27aEnwREQayJb0xCq6+v5c03H2Xz5uP4/cLIUdUsvO5zjBx5ZcxisIRvjDH96NSpY6xc+SQ7d54mEBBKS2u58qqrmTRxMSKxHShpCT9JWGkFYxKH3+9n67b1rFv7BpWVHYgEKBtZw+WfXsiECTfFPNF3soRvjDFRoKocOXKAjRuXs2PHEVpavKSnNzJlSoB5829kZNnVOBP/xY8l/GRgo3SMiYvm5mb27dvBzp3vs3fvMRobBRE/RUWnmT17OJdd9hfk5IyJd5hnWcJPApbrjYlMTc1JHn/8N+TkZDNkyAhKS8dTUFBEdnY2mZmZiAgiQmtrK01NTdTVneb48f2cOHGYqqpqamoCgODxdFBQcIoLpuYx85KFDBt2JR5P4qXXxIvIGGNiZM+e1Zw40UZ6ei379lUDW8N6XEZGI9nZ9UycKIwdO5aJk66gaPCleDy+/g34PFnCTwbWpWNMRFrbmgG47ba55OQIJ08eor6+msbGetraAqh2AF7SM9LIysokP6+EIUNHUpA/idzcSXg8GfF9AX1kCT8JCGLTXhkTgfa2FgAKCqYwdOgFlJXFOaB+ZtUyk4Sle2P6rq2tDYD0jNw4RxIblvCTgWV7YyLS3uEk/Iz0nDhHEhuW8JOF9eEb02ftbgs/I2NQnCOJDUv4ScByvTGRae/owOPpwOtN7NE10WIJP1lY1jemz9rb2/F4AvEOI2Ys4ScDBeJ8yrYxA1FHux+v1xK+GUCcVG9Hbo3pq46OgCV8MxBZC9+Yvuro8OON/kyCCcsSflIQy/fGRKCjQy3hm4HFcr0xkfH7lbQUqjdgCT8ZWC0dYyLS0QFeb+p8eSzhJ4vU+cwaEzV+P6SlpU4aTJ1XmsRsAnNjIuP3C2m+1EmDqfNKk53lfGP6zO/34EtLnaO2lvCTgAT9a4wJnz/gJS2Fjtpawk8GdtDWmD4LBAIE/F58Pkv4/UJEFonIThHZIyLfi+W2k5nlemP6rqOjGRB8vtQonAYxTPgi4gV+CdwATAWWisjUWG0/udmJV8b0VUvLGQB8vvQ4RxI7sfwtcxmwR1X3AYjI48BiYFu0N/SLv/8Z/hSa8s/j89B2qo0Xvv9mvEMxZmDJSGfz87vY9bt74x3JJ/jUy9ee+G3UnzeWCX8EcDjo9hFgdteVROSrwFcBRo0aFdGGMjs8BFIo4QO0d7TgDbTHOwxjBhRvaxtpjS1oglWb9Xv6J55YJvzuXsE5WVlVHwAeAJg1a1ZEWftPf/o3kTzMGGOSWiwP2h4BRgbdLgMqY7h9Y4xJabFM+BuAiSIyVkTSgTuBF2K4fWOMSWkx69JR1Q4R+TrwOuAFHlTVrbHavjHGpLqYnnGgqq8Ar8Rym8YYYxx2pq0xxqQIS/jGGJMiLOEbY0yKsIRvjDEpQlQT94xUEakCDkb48GLgVBTD6Q8DIUawOKNtIMQ5EGIEi7M7o1W1pLsFCZ3wz4eIlKvqrHjH0ZOBECNYnNE2EOIcCDGCxdlX1qVjjDEpwhK+McakiGRO+A/EO4AwDIQYweKMtoEQ50CIESzOPknaPnxjjDGflMwtfGOMMUEs4RtjTIpIuoSfyBOli8gBEdksIhUiUu7eN1hEVojIbvf/wjjE9aCInBSRLUH3hYxLRL7v7t+dInJ9HGO8X0SOuvuzQkRujGeM7nZHishbIrJdRLaKyDfd+xNmf/YQY0LtTxHJFJH1IrLJjfMf3fsTZl/2EmdC7U8AVDVpLjhll/cC44B0YBMwNd5xBcV3ACjuct8/A99zr38P+Fkc4rocmAls6S0unAnoNwEZwFh3f3vjFOP9wHe6WTcuMbrbLgVmutcHAbvceBJmf/YQY0LtT5xZ8nLd6z7gfWBOIu3LXuJMqP2pqknXwj87UbqqtgGdE6UnssXA79zrvwOWxDoAVV0NVHe5O1Rci4HHVbVVVfcDe3D2ezxiDCUuMQKo6jFV/cC93gBsx5nPOWH2Zw8xhhKv91xV9Yx70+delATal73EGUrcPp/JlvC7myi9pw9yrCmwXEQ2upO1AwxV1WPgfBGBIXGL7pNCxZVo+/jrIvKR2+XT+dM+IWIUkTHAJTgtvoTcn11ihATbnyLiFZEK4CSwQlUTcl+GiBMSbH8mW8IPa6L0OJqvqjOBG4C/EJHL4x1QBBJpH/8nMB6YARwD/tW9P+4xikgu8AzwLVWt72nVbu6LSazdxJhw+1NV/ao6A2cO7MtE5KIeVk+0OBNufyZbwk/oidJVtdL9/yTwR5yfcSdEpBTA/f9k/CL8hFBxJcw+VtUT7hctAPyaj38WxzVGEfHhJNJHVfVZ9+6E2p/dxZio+9ONrRZYBSwiwfZlsOA4E3F/JlvCT9iJ0kUkR0QGdV4HFgJbcOL7grvaF4Dn4xPhOULF9QJwp4hkiMhYYCKwPg7xdX7ZO92Csz8hjjGKiAD/A2xX1Z8HLUqY/RkqxkTbnyJSIiIF7vUs4FpgBwm0L3uKM9H2J5Bco3TUOQJ+I86og73A38U7nqC4xuEcmd8EbO2MDSgCVgK73f8HxyG2x3B+crbjtD6W9RQX8Hfu/t0J3BDHGH8PbAY+wvkSlcYzRne7C3B+nn8EVLiXGxNpf/YQY0LtT2Aa8KEbzxbgh+79CbMve4kzofanqlppBWOMSRXJ1qVjjDEmBEv4xhiTIizhG2NMirCEb4wxKcISvjHGpAhL+CYliEiBiHwt6PZwEXm6n7a1RER+GGLZGff/EhF5rT+2b0wolvBNqigAziZ8Va1U1c/107a+C/yqpxVUtQo4JiLz+ykGY85hCd+kip8C49265P8iImPEra0vIveJyHMi8qKI7BeRr4vIX4vIhyKyTkQGu+uNF5HX3OJ374jIlK4bEZFJQKuqnnJvjxWRtSKyQUT+d5fVnwPu6tdXbUwQS/gmVXwP2KuqM1T1b7pZfhHwJzj1Tn4MNKnqJcBa4F53nQeAv1TVS4Hv0H0rfj7wQdDtXwD/qaqfAo53Wbcc+HSEr8eYPkuLdwDGJIi31KkN3yAidcCL7v2bgWluZcl5wFNOKRrAmcCiq1KgKuj2fOA29/rvgZ8FLTsJDI9O+Mb0zhK+MY7WoOuBoNsBnO+JB6hVpwRuT5qB/C73hapfkumub0xMWJeOSRUNONP5RUSdevH7ReR2cCpOisj0blbdDkwIuv0eTtVWOLe/fhIfV1A0pt9ZwjcpQVVPA++JyBYR+ZcIn+YuYJmIdFY87W76zNXAJfJxv883cSa72cC5Lf+rgJcjjMWYPrNqmcZEmYj8AnhRVd/oZb3VwGJVrYlNZCbVWQvfmOj7P0B2TyuISAnwc0v2JpashW+MMSnCWvjGGJMiLOEbY0yKsIRvjDEpwhK+McakCEv4xhiTIv4/+bJ4Eib1UgwAAAAASUVORK5CYII=\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": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAElCAYAAADnZln1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnd0lEQVR4nO3deZxkVXn/8c+3enp6RkBGYMBhHQQBlQCScQuIEheQaNAYl7jEhaDGGE1EcYkvg6gR488tv7iEoIJiNCouSBQUgaCGKIMMm4iyzzADDMs4AzPT3dP15I9zqrump6u6605113K/79erum7d9al7q54+de655yoiMDOz/lfpdABmZjY3nPDNzErCCd/MrCSc8M3MSsIJ38ysJJzwzcxKwgm/hCSdJuncPLyvpIckDXQ6rmYkPV3STXO8zZB04Hau4wZJz2xPRNusu+FxlLSHpMslbZD0cSVfkvSgpF/ORjzW/Zzwe5Ck2yU9e9K410r6Wavriog7I2LHiBhrX4StmUlijYifRsTBcxVTu0TEEyLiMtg6Qc/CdiYfxzcA9wGPjIhTgKOB5wB7R8STZyMG635O+Nb1JM3rdAw9aD/g1zFxZeV+wO0R8XCrK/L+7x9O+H1K0p6SzpO0VtJtkt7aYL6luYQ9r2658yU9IOlmSSfXzTsg6b2SbslVBVdJ2idPO0TSj/NyN0l6ad1yZ0v6jKT/ysv9QtIBedrlebZrcpXEyyQ9U9IqSe+SdDfwpdq4unXuI+nb+f3dL+lfG+yDTZJ2qRv3REn3SRrMr18v6cZc1XGRpP0a7KedJX05b+8OSe+TVKmbfnJezwZJv5Z0ZB5/u6RnSzoeeC/wsvw+r5H0EklXTdrOKZK+2yCG/SX9d97Gj4HdpjqOks4GXgOcmrf1RuAs4Gn59QfyMs+XtELSOkn/I+mwuvXdnvf/tcDDeb1PzfOty/E/s27+yyR9UNLPc3w/klQf39F1y66U9No8fkjS/5N0p6R7JH1e0sI8bTdJF+RlHpD00/p9bgVEhB899gBuB549adxrgZ/l4QpwFfB+YD7wGOBW4Lg8/TTg3Dy8FAhgXn7938BngQXAEcBa4Fl52juB64CDAQGHA7sCOwArgdcB84AjSdUJT8jLnQ08ADw5T/8q8PW62AM4sO71M4EtwEeBIWBhHrcqTx8ArgE+mbe9ADi6wb66BDi57vXHgM/n4RcCNwOPy3G9D/ifqeICvgx8D9gp77PfAiflaS8B7gKelPfLgcB+k49V/X7Pr4fyfnlc3birgRc3eC9XAJ/Iyx0DbGhyHM8GPjTV5yO/PhK4F3hK3p+vybEO1cW9Atgn7/+9gPuBE0ifr+fk14vz/JcBtwAH5fkvA87I0/bNsf4FMEj6zByRp30KOB/YJe/b7wMfydM+Anw+LzMIPB1Qp79/vfzoeAB+FDho6cv4ELCu7rGRiYT/FODOScu8B/hSHh5PPPWJIn+5x4Cd6pb7CHB2Hr4JOHGKeF4G/HTSuH8D/jEPnw2cVTftBOA3da+nSvgjwIJJ42oJ/2mkf0TzZrCv/gq4JA+L9I/pmPz6h+SknV9X8n7crz4uUkIcBh5fN+8bgcvy8EXA25ocqykTfh73OeDDefgJwIPkpDtpvn1J/wR3qBv3H1Mdx7p93izhfw744KRt3AQ8oy7u19dNexfwlUnzXwS8Jg9fBryvbtqbgQvrPnvfmeI9CXgYOKBu3NOA2/Lw6aR/sgdOXtaPYg//POpdL4yIRbUH6QtWsx+wZ/4pvE7SOlJ1wh7TrHNP4IGI2FA37g5S6Q7SP4RbplhuP+Apk7b3SuDRdfPcXTe8EdhxmljWRsTmBtP2Ae6IiC3TrAPgW6SqjD1JpeIAfloX96frYn6AlIT2mrSO3Ui/lO6oGzeT/TIT5wCvkCTg1cA3ImJ4ivn2BB6Mrevg75hivpnaDzhl0jHbJ2+nZuWk+V8yaf6jgSV18zQ6xo32z2LgEcBVdeu8MI+H9GvsZuBHkm6V9O7W36bV88mY/rSSVEp6bIvLrQZ2kbRTXdLfl1RdUVvvAcD1U2zvvyPiOUUDnkKzblxXAvtKmjdd0o+IdZJ+BLyUVHXztcjFx7yeD0fEV6eJ5T5glHwiNI+bar9MZ5v3FBH/K2mEVF3xivyYyhrgUZJ2qEv6+061zhmqvfcPzzDelaQS/smNZp5mW1O1DLoP2ESq+rtr8sT8GTyF9I/pCcClkq6MiJ8UiMHwSdt+9UtgfT7ptlDpZOuhkp7UbKGIWAn8D/ARSQvySbyTSHXukE78fVDSY5UcJmlX4ALgIEmvljSYH0+S9LgZxnsP6TxDK+9vDXCGpB1yrEc1mf8/gL8EXpyHaz4PvCcnk9qJ2ZdMXjhSU8dvAB+WtJPSid23A7UmlmcB75D0h3m/HKipT/7eAyyd4sTjl4F/BbZExJRNayPiDmA58AFJ8yUdDbygyXuezr8Db5L0lBzzDpL+RNJODeY/F3iBpOPy52mB0on0vWewra8Cz5b00nzyd1dJR0RENcfxSUm7A0jaS9Jxefj5eV8KWE+qbuxY8+F+4ITfh3KCegHppOttpJLUWcDOM1j8L0j1wauB75Dq4X+cp32ClPh+RPoCfgFYmEtizwVenpe7m4kTrjNxGnBO/ln/0ulmrnt/BwJ3AqtI5xEaOR94LHBPRFxTt57v5Di/Lmk96ZfL8xqs429J9c23Aj8j/eP4Yl7PN4EP53EbgO+STkJO9s38fL+kX9WN/wpwaH5u5hWk8zMPAP9I+kdRSEQsB04m/aN5kFR18tom868ETiRVDa4lldrfyQxySETcSTpvc0qOfQXphD+kcwM3A/+bj8HFpEYBkI7ZxaTzVVcAn418TYMVo4lft2bWCbkZ4r3AkRHxu07HY/3LJXyzzvtr4Eone5ttPmlr1kGSbie1DHphZyOxMnCVjplZSbhKx8ysJJzwzQqQ9Mrcvn+6+Wath8wilPo1+lCn47DOcMK3WaeJvtprj5D0cN3rpxdY5zZdRE+a/kxJ1bz+DUodur2uYPxbdTAHEBFfjYjnFlmfWaf4pK3NutwOe7wrBUkBHB4RN8/ypldHxN75wp0TgW9J+kVE/Hq6BWvkroGtj7iEbx2lAt3jSvoKqVuB7+cS/KnNthHJd0kXGD0+X1F6taT1Sl31nlYXT600f5KkO0m9bda6cF6Xt/c0TbrhjKQnaKJ76HskvbfB+23WxfBrlfqM2aDUpfUrm+yzT0lanR+fkjSUp9W6lj5F0r2S1jT6ZSPpekkvqHs9qNR19BHN9qf1Lid867SPkrrUPYJ05exepG6dIV2ZuYrUmdYepKs8IyJeTbrC9gWR7vL0z802kP9JvAhYROre+WFSVwuLgD8B/lrSCyct9gxS3zvHkTpdA1iUt3fFpPXvRLoi9EJS52MHAtv09yJpL+C/gA+RrsR9B3CepMWSdgD+BXheROwE/BHpitSp/APwVNI+O5zUT8376qY/mnRV9V6krjE+I+lRU6zny8Cr6l6fAKyJiEbbtR7X9Qlf0hdzSWVyh11F13dhLl1dMGn8F3KJ61pJ35I0XW+Otp1yVcvJwN9HRK2Xzn8iddEAqcOyJaTuikcj3eawlXbEeyr1wHgfqSuCV0fETRFxWURcFxHViLgW+Bopwdc7LSIejohNM9jO84G7I+LjEbE5IjZExC+mmO9VwA8i4gd52z8m9Y9zQp5eBQ6VtDAi1kTEDQ2290rg9Ii4NyLWAh8g9bRZM5qnj0bED0hdE0x1e8hzgRMkPTK/fjXTd+9gPazrEz6pX+/j27i+j7H1l6Pm7yPi8Ig4jFR6fEsbt2lTm+3ucVfn7qN3iYgjIuLrAEodhl2qdPeq3wNvou7uUdnKbdbW2Ey7R27YxXDuAfNlOZY1SncHO6TBevZk266a67s1vn9SL6JTdkcdEauBnwMvlrSI1I/QdD2HWg/r+oQfEZeTOlwaJ+mAXFK/KtfrNvpiTLW+n5A6uJo8fn1et0h37PEVabOvvnvcWt/+O0fEjpC6x42IUyLiMaTO0t4u6Vl52e05Pv9B6lBtn4jYmdRrpibNEw2GpzLT7pFrXQwvqnvsEBFnAETERbmL6SXAb0g9SU5lNemfR82+eVwR55B+ebwEuGKqboqtf3R9wm/gTOBvI+IPSfWgn23HSiV9idTT4yHA/2/HOq2x7ewet9UulevtRLrRy2ZJT6ZxH/Q1a0nVLY22dwHwaEl/l0+o7iTpKVPM17CLYUl7SPrTXJc/TKqGadQV8NeA9+W6/91I5zyKtvX/Lul2h29jO3rftN7Qcwk/163/EfBNSStIt9Jbkqf9WW55MPlx0UzWHRGvI/00vpHm3e1a+xTtHvcjpKS3TtI7Wtzmm4HTJW0gJctvNJs5IjaSuj/+ed7eUydN30C6x+sLSAWG3wHHTrGeZl0MV0gnqVeTftE+g63vYlbvQ6S6/2tJJ6F/lce1LJ+jOA/YH/h2kXVY7+iJvnQkLQUuiIhD8wmmmyJiyTSLNVvfM4F3RMTzG0x/BvDORtPN+omk9wMHRcSrpp3ZelrPlfBzXfttyncmUnL4NIs1lddxYG2YVFL7zXYHa9blJO1Carp5ZqdjsdnX9Qlf0tdIP+cPzheUnERqlnaSpGuAG0g/k2e6vp+S7jz0rLy+40gn7M6RdB3pJ/IS4PQ2vxWzriLpZFK10g9z4wjrcz1RpWNmZtuv60v4ZmbWHl3dMdRuu+0WS5cu7XQYZmY946qrrrovIhZPNa2rE/7SpUtZvnx5p8MwM+sZku5oNM1VOmZmJeGEb2ZWEk74ZmYl4YRvZlYSTvhmZiXhhG9mVhJO+GZmJdHV7fDNrPOGh4e58sorGR0d7XQoc65TXc/Mnz+fo48+uu3rdcI3s6ZuvfVWLr744k6HUSo77rijE76Zzb2xsXTjrTe/+c3svvvuHY7Gtofr8M2sqWq1CsDAwECHI7Ht5YRvZk3VEn6l4nTR63wEzawpJ/z+4SNoZk054fcPH0Eza8oJv3/4CJpZU074/cNH0MyacsLvHz6CZtaUE37/8BE0s6ac8PuHj6CZNVVL+JI6HIltLyd8M2uqWq0iySX8PuAjaGZNVatVJ/s+4c7TzKypysYqzxo+lLVfuK7TobSuM70bb7fKwnns+srHtX29Tvhm1tSCdWLfLbtRfWgUDfZoSb/Hzj/EvOqsrNcJ38yaG0vF5F1ecQiDix/R4WBse/Tov2szmzPVlPA14HTR63wEzaypyAmfgd6qFrFtOeGbWXO5SkcVJ/xe54RvZs3l84dyCb/nOeGbWXPjVTpOF71uzo+gpAFJV0u6YK63bWYFjJ+0dQm/13XiX/bbgBs7sF0zK6LWJNx1+D1vThO+pL2BPwHOmsvtmllxqgZB+KRtH5jrEv6ngFOZKDOYWberQlU92keBbWXOEr6k5wP3RsRV08z3BknLJS1fu3btHEVnZg1Ve7ZLGptkLkv4RwF/Kul24OvAH0s6d/JMEXFmRCyLiGWLFy+ew/DMbCqqQriE3xfmLOFHxHsiYu+IWAq8HLgkIl41V9s3s4LCVTr9wg1rzaypVMLvdBTWDh3pLTMiLgMu68S2zaw1Clfp9AuX8M2sOZfw+4YTvpk1pYBwpugLPoxm1pRCLuH3CSd8M2tKATjh9wUnfDNrquIqnb7hw2hmTSnkhN8nfBjNrCmFXKXTJ5zwzawpt9LpHz6MZtZUBTlT9AkfRjNrSiHf/KRPOOGbWVOVcAm/X/gwmllTwiX8fuGEb2ZNVVyl0zec8M2soWq16pO2fcSH0cwaioiU8Adcwu8HTvhm1lAq4VdcpdMnnPDNrCFX6fQXH0Yza2hsbIwBKq7S6RNO+GbWUHVLFQBVnCr6gY+imTVU3TKWBlzC7wtO+GbWUHU0JXw54fcFJ3wza2i8hO9WOn1hXqcDsPb41a9+xfr168dfS9t+Qds9rtn47eX1dsd6N9//MPszgAZcNuwHTvh94OGHH+b888/vdBjWh3aIIfbnaBbssLDToVgbOOH3gWq1yuO27MVhj3kCS5YsSSNzAS4itinNBdF8hePLwuSC4JRLqsGLrcbHNNMbra9+vLaZZZv30qjkOoN1biUa7KM2FLin2futb6LueLVbbNrC5svWsHiPxe1fuc05J/w+MPbQKEdtOYT47RgP3XJXGrnVlz+mHKTBLGaTDew81OkQrA2c8PtAVFNb6Q1HDvL4lz61feuNVv9RFJ+/Yek0thmY2bobrHC7SsGzUYTeZhuzv4mWtzEgBnYYnJVQbG454feD8dqS9p7g26oqaEarLr59twExm33TJnxJ+85wXesiYv30s1m7RdX1MWY2vZmU8M8hlSGbFcICOBv4chtishZF1C5/dznZzBqbNuFHxLGTx0l6dETcPTshWavGS/iz1GbbzPpD0asp/rKtUZiZ2awretL2REkbgR9HxE3tDMhaVyvhu4BvZs0ULeH/GXAz8CJJZ7UxHitgvPmk6/DNrIlCJfyIuAe4MD+s08ar8J3wzayxQiV8SZ+RdHYefu4Ml1kg6ZeSrpF0g6QPFNm2batWpePGmWbWTNEqnRHg1jz8xzNcZhj444g4HDgCOF5S+y4LLbNas0yX8M2siaInbTcCO0saBGZ0YVakiuaH8svB/HChtA3GT9q6B1sza6JoingAuAX4DPDzmS4kaUDSCuBeUgufX0wxzxskLZe0fO3atQXDK5lZ6lrBzPpLSyV8SYuATwIHA+eSrqw9aabLR8QYcERez3ckHRoR10+a50zgTIBly5b5F8AMjLfScb43a0m1uoUbbriA4eGNnQ5lK/PmDXHEES9u/3pbmTki1kk6A1gK3AccBny71Y3m9VwGHA9cP83sNo2JdvjO+GatuOHXF3HeeSs6HcY25s8f7nzCz04CbouIi4CrZrqQpMXAaE72C4FnAx8tsH2bZKKE74Rv1orhzZsAeMYz9mPPPQ/qcDQTBgZmpzvqIgn/QeBNkg4GrgFWRMTVM1huCXCOpAHSuYNvRMQFBbZvk7kdvlkh1dzCba+9DuCgg47qcDSzr+WEHxEfkfQT4Lek5pXHANMm/Ii4Fnhiq9uz6bkO36yY2s2DytLTbMsJX9LpwACwglS6v6zNMVmrXIdvVkjtnsiVkrRpbvldRsT7SRdRVYAXS/r3tkdlLXEJ36yYiWtYBjocydwo+m/ti8DjgF2Bz7YvHCuklu9L8rPUrF2CXKVTktJS0YT/VlJ10Dzg0+0Lx4qIubi5tlk/inJVhxZN+LcAC4DvRcQxbYzHCnCzTLNiat+diqt0mroBuAQ4SdKVbYzHCpioh3TCN2tFNbfSKcsJsKKdpx1Aao9/Zn62LlCWn6Vm7TJRwi9HK52iCX9lRFwiaQmpIzTrpKpb6ZgVUWuWKTfLbOp4SXsDnyd1pmYdFCU78WTWNuMNHsrx3Sma8BcB7wJOJbXJtw7yPW3Niql1rVCWk7ZFq3ROBw6JiJskjbUzICvArTLNihlv4OYqna1IOrw2HBGrIuLiPPzu2QjMWhBupWNWRIzfHtQJf7KrJV0r6VRJ+8xaRNYy94dvVkzZzn+1kvA/DuwAnAHcJulSSa+fnbCsFWX70Jq1iy+8aiAi3hkRBwDLgLNI3SKfOVuBWQvK1dDArG0mCkvlqNKZ8UlbSbsCLwL+HDiWlF7unKW4rBXuWsGsECf8xu4m/SJ4EPgScG5E/GxWorKWhE/amhUy3lumE/42vgOcC/wwIkZnKR7bDq7DN2tNrcGDu1aYJCJeOpuB2HZw1wpmhZStSqcc77LPRckuHjFrHyf8piS9YDYCseJch29WzHivJCWp0inyLj/c9ihs+4w3y3TCN2tFrS+dslR2FHmXzirdxjcxNyumZP3hF3mX7qqry7hKx6yYifNfvtLWeoX70jErxK10rGe5hG/WmhjvD78cqbDIu7yn7VHYdglXspkV4iqdaUTEc2YjENsOrsM3K2Tinrbl+O6U43dMn3N/+GbFTHx3XMK3HqOS1EOatY+bZU5L0tvrhg9uXzhWiCvxzQpJrXTK8/1p6SbmkhYBnwQOkbQZuBY4CXhd+0OzmRo/8eQ6fLOWOOE3ERHrgNdJOg64DzgM+PYsxGWt8Elbs0ICkJzwpzMaEVdJWg3c286ArIDxpmVO+GYtKVl1aNEzFcdL2hv4PKmKZ1qS9sk3Pr9R0g2S3lZw2zaZb3FoVkg1olQl/KIJfxHwLuBUYHiGy2wBTomIxwFPBf5G0uMLbt/quA7frKCSlfCLVumcDhwcETdJGpvJAhGxBliThzdIuhHYC/h1wRisJtwO36yIiHLV4Rct4b8HeHUevrTVhSUtBZ4I/GKKaW+QtFzS8rVr1xYMr5xcwjdrTZSohQ4UT/gjwK15+NhWFpS0I3Ae8HcRsX7y9Ig4MyKWRcSyxYsXFwyvXKLqOnyzIqIapfraFE34G4GdJQ0C+850oTz/ecBXI8LNOdvFdfhm26E8pfyiCf8fgVuAzwBfnckCShXMXwBujIhPFNyuTcV1+GaFlOycbeGTtm+tJe0WulY4ilTvf52kFXnceyPiBwVjsMmc781aEiVrllmka4XPAfvlrhWuAf6KGXStEBE/wylpVkw0wy9HB1Bm7VOeZA8FulaQtAq4nNTC5nDctULnuWsFs0LShVedjmLuFKnSuR94E3AwqYS/qq0RWetyIaXihG/Wmhj/UwotJ/yIOEPSJcBvgSOApwNXtzkuK6JMRRWzNkgXXnU6irnTcsKXdDowAKwAVkTEZW2OyVpVDUBupWPWorJdeFWkhP9+SXuQrpR9saQDIuLk9odmM5U+tHIdvlmr3EpnRt4I/FtEXNjOYKwgX3hlVojb4c/MF4G/lrQD6arZFe0LyVpW+9C6SsesJWWrwy/acPutpH8W84B/aV84VoivtDUrpGx1+EUT/i3AAuB7EXFMG+OxIlylY1aMu0eekRuAS4CTJF3ZxnhsO7iEb9aa8IVXM3IA8CBwZn62DoqAasl+mpq1Q9m+NUUT/sqIuETSEnwT886LoHwfXbPt55O2M9PyTcxtFoXTvVkhJWuXOZc3MTcz6you4Tcg6fC6l6eTWujcBMzoJuY2iyJK17zMrB3K9q1ppYR/taRrJZ0KKCIuBoiId89OaDZjZfvUmrWLS/gNfRzYATgDuE3SpZJePzthWUtch29WSLgOf2oR8c6IOABYBpwFHENqlmldoVwfXLN2KVMJf8bNMiXtCrwI+HPgWNLtCu+cpbisFRGE7x5p1rKynbRtpR3+3aRfBA8CXwLOzfeptU5z4d6skLJ9dVpJ+N8BzgV+GBGjsxSPFREQJeoPxKxt3LXC1iTtmwffkZ+XNOizZV1ErG9XYGZms61k52xnVMI/h7oe1xvME8DZwJfbEJO1ynX4ZoW5hF8nIo6di0CsuJIVUszaJtLdQUujaNcK1k2ifDdyMGuHiFLleyf8fqDxP2bWqjJV6Tjh94Nw+d6siLKdtHXC7wNl+9CatUu68Ko8RXwn/H4Q43/MrEUlyvdO+P1AON2bFVG2X8dO+P2gZE3LzNolcAnfzKwcStZ5mhN+P3ArHbNC0vemPBnfCb9vOOWbtcwl/Nkh6YuS7pV0/VxtszQCokQfWrN2csKfHWcDx8/h9srDhXuzQtxKZ5ZExOXAA3O1PTOz6fjCqw6T9AZJyyUtX7t2bafD6RklK6iYtU2J8n33JfyIODMilkXEssWLF3c6nN7gdvhmhaQqnfJ8ebou4VsBLt6bFSSX8K23pK4VnPXNWhVuljk7JH0NuAI4WNIqSSfN1bb7nqt0zAor00nbmdzTti0i4i/maltmZjPhZpnWe8r2qTVrozKV8J3w+4J8pa1ZARE+aWtmViLlyfhO+H1ArtExK8StdKznBOEqHbNC5Dp86y1ytjcrxCV8600l+tCatY9KlfGd8PuB6/DNCnFvmdaDXIdvVlSZvjpO+H1ApfrImrVPaodfnu+PE34/iPE/ZtYqJ3zrOSX60Jq1S5Ts97ETfj9w4d6sGJ+0tV4j8Elbs4Kc8M3MSsAnba33+AYoZgWpVN8dJ3wzK62IcjVrdsLvA+X5uJq1m1ClPN8gJ/x+ED5pa9aqarVKucr3Tvh9QbgZvlnrUntmqTxpsDzvtM+5Kb5ZayKqQLkKS074/SDK1dLArB2q1bE8VJ4vjxO+mZVSREr4FZ+0tV6i8T9mNlMxXg9ani+PE36/KFNFpFkb1Er4vtLWeorcPbJZy2p1+E741oPK86E1a4eIWrPM8nx3nPD7glvpmLVqollmeb48Tvh9wCdtzVrnKh3rTa6+NyvAVTrWq8rzmTVrizKW8Od1OgDbfuXq/smsPWrNMsvUpNkl/H5Rns+sWVvULrwqU4HJJfw+IN/xyqxl4xdeNelaISK4e2SU2zaOcNumYW7dNMztm4Z5eEuV+RWxx9Agew0NsveC+SxdOMTShUPsOjjQtdVEc5rwJR0PfBoYAM6KiDPmcvtmNjMPbRljzfAo94yMctfmUVYPj3DX5lHuys/3jIwyFkFF8IjKALvPn8ceQ4PsMX+Q3Yfm8ej5g+Ov9xiax26Dgwx2WZ81W3Id/vqxAa7ZsJG7Nqf3tmp4hFWbR7ht4zC3bRphU7U6vsygxH4L5/PIeQMMj1a5av1G7h/dstV6dxyosHThEPstTP8E9pg/j0cNpscu8wZYNDiPHQcqDFXEUCU9z9U/iDlL+JIGgM8AzwFWAVdKOj8ifj1XMfQrodrltnMmIghSO4eI/EwaVw2oRhpTjYnpKdb0CKCSBypKAyK9rn34az+5Iy898XoihtrriIltVQOq1G0/x5SeqxMxAqMRjFSD0QhGq1VGIhiNifHDUWVTtcrmarBpLNgUVTaPpXGbqpGmjVXZFMFwtcqgxHyJ+RWl4UqFQYmhgfR6UOlLPlipzVdhfmVimfRIy4ynAOX3rvr3PrG/tx6XHlsiGKlWJ95f7T1GGh7J0zdsqbJhbIz1WyYe945s4eGxKooqFapUIqhQZfd5FfYcGuAP5s/jOY+oMB8RjLFpy2buGx7h/nUjrBgeZf3oyPhyAzFGhfS8E6M8kjEWaQuPZAs7xRYewSgLqiMMxWh6VEcYGhthqDrM0Ngw86sjDI5tZt5Yeh4Ym3hdiSqjA0OMDgyxJT8PV4bYPLCAzZX5bB5YwMbKAjYNDLGxsoCNlfk8VBniQebzQAyycSQ4FPjuqip3X3oBY6owpgHmVQbYdWiIP1iwkON2XMCeC4fYe+FC9lm4gEfPn0el1n9+/pwOV4PVwyOs3LyFlcMjrBoeZeXwRu58YB1XjGxhNKL26Z54rj+4wDylxD9YqTBPYtfBQX74tMNm/oWcIcVED0KzStLTgNMi4rj8+j0AEfGRRsssW7Ysli9f3vK2PvW+M6h2V2FiVlVUYaQ6zOaxjZ0OxayHiOrQfBauXsPQ73/f6WC2UiF463/+Z6FlJV0VEcummjaXVTp7ASvrXq8CnjJ5JklvAN4AsO+++xba0MItFapd9vNxVgWMbhlloO6np5lNb2B4I4ObKlB55DbTamXyzhidlbXOZcKfKgNvszcj4kzgTEgl/CIbeuMZpxZZzMysr81ls8xVwD51r/cGVs/h9s3MSm0uE/6VwGMl7S9pPvBy4Pw53L6ZWanNWZVORGyR9BbgIlKzzC9GxA1ztX0zs7Kb03b4EfED4AdzuU0zM0vctYKZWUk44ZuZlYQTvplZSTjhm5mVxJx1rVCEpLXAHQUX3w24r43hzIZeiBEcZ7v1Qpy9ECM4zqnsFxGLp5rQ1Ql/e0ha3qg/iW7RCzGC42y3XoizF2IEx9kqV+mYmZWEE76ZWUn0c8I/s9MBzEAvxAiOs916Ic5eiBEcZ0v6tg7fzMy21s8lfDMzq+OEb2ZWEn2X8CUdL+kmSTdLenen46kn6XZJ10laIWl5HreLpB9L+l1+flQH4vqipHslXV83rmFckt6T9+9Nko7rYIynSbor788Vkk7oZIx5u/tIulTSjZJukPS2PL5r9meTGLtqf0paIOmXkq7JcX4gj++afTlNnF21P4F8M+o+eZC6Xb4FeAwwH7gGeHyn46qL73Zgt0nj/hl4dx5+N/DRDsR1DHAkcP10cQGPz/t1CNg/7++BDsV4GvCOKebtSIx520uAI/PwTsBvczxdsz+bxNhV+5N0l7wd8/Ag8Avgqd20L6eJs6v2Z0T0XQn/ycDNEXFrRIwAXwdO7HBM0zkROCcPnwO8cK4DiIjLgQcmjW4U14nA1yNiOCJuA24m7fdOxNhIR2IEiIg1EfGrPLwBuJF0P+eu2Z9NYmykU8c8IuKh/HIwP4Iu2pfTxNlIxz6f/Zbwp7pRerMP8lwL4EeSrso3awfYIyLWQPoiArt3LLqtNYqr2/bxWyRdm6t8aj/tuyJGSUuBJ5JKfF25PyfFCF22PyUNSFoB3Av8OCK6cl82iBO6bH/2W8Kf0Y3SO+ioiDgSeB7wN5KO6XRABXTTPv4ccABwBLAG+Hge3/EYJe0InAf8XUSsbzbrFOPmJNYpYuy6/RkRYxFxBOke2E+WdGiT2bstzq7bn/2W8Lv6RukRsTo/3wt8h/Qz7h5JSwDy872di3ArjeLqmn0cEffkL1oV+HcmfhZ3NEZJg6RE+tWI+HYe3VX7c6oYu3V/5tjWAZcBx9Nl+7JefZzduD/7LeF37Y3SJe0gaafaMPBc4HpSfK/Js70G+F5nItxGo7jOB14uaUjS/sBjgV92IL7al73mRaT9CR2MUZKALwA3RsQn6iZ1zf5sFGO37U9JiyUtysMLgWcDv6GL9mWzOLttfwL91Uon0hnwE0itDm4B/qHT8dTF9RjSmflrgBtqsQG7Aj8Bfpefd+lAbF8j/eQcJZU+TmoWF/APef/eBDyvgzF+BbgOuJb0JVrSyRjzdo8m/Ty/FliRHyd00/5sEmNX7U/gMODqHM/1wPvz+K7Zl9PE2VX7MyLctYKZWVn0W5WOmZk14IRvZlYSTvhmZiXhhG9mVhJO+GZmJeGEb6UgaZGkN9e93lPSt2ZpWy+U9P4G0x7Kz4slXTgb2zdrxAnfymIRMJ7wI2J1RPz5LG3rVOCzzWaIiLXAGklHzVIMZttwwreyOAM4IPdL/jFJS5X71pf0WknflfR9SbdJeoukt0u6WtL/Stolz3eApAtz53c/lXTI5I1IOggYjoj78uv9JV0h6UpJH5w0+3eBV87quzar44RvZfFu4JaIOCIi3jnF9EOBV5D6O/kwsDEinghcAfxlnudM4G8j4g+BdzB1Kf4o4Fd1rz8NfC4ingTcPWne5cDTC74fs5bN63QAZl3i0kh9w2+Q9Hvg+3n8dcBhuWfJPwK+mbqiAdINLCZbAqyte30U8OI8/BXgo3XT7gX2bE/4ZtNzwjdLhuuGq3Wvq6TvSQVYF6kL3GY2ATtPGteo/5IFeX6zOeEqHSuLDaTb+RUSqb/42yS9BFKPk5IOn2LWG4ED617/nNRrK2xbX38QEz0oms06J3wrhYi4H/i5pOslfazgal4JnCSp1uPpVLfPvBx4oibqfd5GutnNlWxb8j8W+K+CsZi1zL1lmrWZpE8D34+Ii6eZ73LgxIh4cG4is7JzCd+s/f4JeESzGSQtBj7hZG9zySV8M7OScAnfzKwknPDNzErCCd/MrCSc8M3MSsIJ38ysJP4P4BeURiIeUWkAAAAASUVORK5CYII=\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": 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_mars_disk/param.full.in b/examples/symba_mars_disk/param.full.in new file mode 100644 index 000000000..ca73467f0 --- /dev/null +++ b/examples/symba_mars_disk/param.full.in @@ -0,0 +1,35 @@ +!Parameter file for the SyMBA-RINGMOONS test +T0 0.0 +TSTOP 2400.0 +DT 600.0 +PL_IN mars.in +TP_IN tp.in +IN_TYPE ASCII +ISTEP_OUT 1 +ISTEP_DUMP 1 +BIN_OUT bin.dat +PARTICLE_FILE particle.dat +OUT_TYPE REAL8 +OUT_FORM EL +OUT_STAT REPLACE +J2 0.0 +J4 0.0 +CHK_CLOSE yes +CHK_RMIN 3389500.0 +CHK_RMAX 3389500000.0 +CHK_EJECT 3389500000.0 +CHK_QMIN 3389500.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 3389500.0 338950000000.0 +ENC_OUT /dev/null +EXTRA_FORCE no +BIG_DISCARD no +RHILL_PRESENT yes +MTINY 1000.0 +ENERGY yes +FRAGMENTATION yes +ROTATION yes +MU2KG 1.0 +DU2M 1.0 +TU2S 1.0 +SEED 2 3080983 2220830 diff --git a/examples/symba_mars_disk/param.in b/examples/symba_mars_disk/param.in index 119d43069..b96460880 100644 --- a/examples/symba_mars_disk/param.in +++ b/examples/symba_mars_disk/param.in @@ -1,6 +1,6 @@ !Parameter file for the SyMBA-RINGMOONS test T0 0.0 -TSTOP 6e8 +TSTOP 6000.0 DT 600.0 PL_IN mars.in TP_IN tp.in diff --git a/examples/symba_mars_disk/param.partial.in b/examples/symba_mars_disk/param.partial.in new file mode 100644 index 000000000..482074183 --- /dev/null +++ b/examples/symba_mars_disk/param.partial.in @@ -0,0 +1,35 @@ +!Parameter file for the SyMBA-RINGMOONS test +T0 0.0 +TSTOP 1200.0 +DT 600.0 +PL_IN mars.in +TP_IN tp.in +IN_TYPE ASCII +ISTEP_OUT 1 +ISTEP_DUMP 1 +BIN_OUT bin.dat +PARTICLE_FILE particle.dat +OUT_TYPE REAL8 +OUT_FORM EL +OUT_STAT REPLACE +J2 0.0 +J4 0.0 +CHK_CLOSE yes +CHK_RMIN 3389500.0 +CHK_RMAX 3389500000.0 +CHK_EJECT 3389500000.0 +CHK_QMIN 3389500.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 3389500.0 338950000000.0 +ENC_OUT /dev/null +EXTRA_FORCE no +BIG_DISCARD no +RHILL_PRESENT yes +MTINY 1000.0 +ENERGY yes +FRAGMENTATION yes +ROTATION yes +MU2KG 1.0 +DU2M 1.0 +TU2S 1.0 +SEED 2 3080983 2220830 diff --git a/examples/symba_mars_disk/param.restart.in b/examples/symba_mars_disk/param.restart.in new file mode 100644 index 000000000..2ae441a01 --- /dev/null +++ b/examples/symba_mars_disk/param.restart.in @@ -0,0 +1,37 @@ +TSTOP 2.40000000000000000E+03 +T0 1.20000000000000000E+03 +DT 6.00000000000000000E+02 +PL_IN dump_pl2.bin +TP_in dump_tp2.bin +IN_TYPE REAL8 +ISTEP_OUT 1 +BIN_OUT bin.dat +PARTICLE_FILE particle.dat +OUT_TYPE REAL8 +OUT_FORM EL +OUT_STAT APPEND +ENC_OUT /dev/null +ISTEP_DUMP 1 +CHK_RMIN 3.38950000000000000E+06 +CHK_RMAX 3.38950000000000000E+09 +CHK_EJECT 3.38950000000000000E+09 +CHK_QMIN 3.38950000000000000E+06 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 3.38950000000000000E+06 3.38950000000000000E+11 +MTINY 1.00000000000000000E+03 +MU2KG 1.00000000000000000E+00 +TU2S 1.00000000000000000E+00 +DU2M 1.00000000000000000E+00 +SEED 2 3080983 2220830 +EXTRA_FORCE F +BIG_DISCARD F +RHILL_PRESENT T +CHK_CLOSE T +FRAGMENTATION T +ENERGY T +ROTATION T +TIDES F +GR F +YARKOVSKY F +YORP F +RINGMOONS F diff --git a/src/io/io_conservation_report.f90 b/src/io/io_conservation_report.f90 index 820569bd3..6a666e0ff 100644 --- a/src/io/io_conservation_report.f90 +++ b/src/io/io_conservation_report.f90 @@ -12,19 +12,16 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, type(symba_pl), intent(inout) :: symba_plA !! Swiftest planet data structure integer(I4B), intent(in) :: npl !! Number of massive bodies real(DP), intent(in) :: j2rp2, j4rp4 !! Central body oblateness terms - type(user_input_parameters), intent(in) :: param !! Input colleciton of user-defined parameters + type(user_input_parameters), intent(inout) :: param !! Input colleciton of user-defined parameters logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen ! Internals real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now - real(DP), dimension(NDIM), save :: Ltot_orig, Lorbit_orig, Lspin_orig real(DP), dimension(NDIM), save :: Ltot_last, Lorbit_last, Lspin_last - real(DP), save :: Eorbit_orig, Mtot_orig, Lmag_orig real(DP), save :: ke_orbit_last, ke_spin_last, pe_last, Eorbit_last real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now real(DP) :: Eorbit_error, Etotal_error, Ecoll_error real(DP) :: Mtot_now, Merror real(DP) :: Lmag_now, Lerror - logical, save :: lfirst = .true. character(len=*), parameter :: egyfmt = '(ES23.16,10(",",ES23.16,:))' ! Format code for all simulation output character(len=*), parameter :: egyheader = '("t,Eorbit,Ecollisions,Lx,Ly,Lz,Mtot")' integer(I4B), parameter :: egyiu = 72 @@ -35,7 +32,9 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, associate(Ecollisions => symba_plA%helio%swiftest%Ecollisions, Lescape => symba_plA%helio%swiftest%Lescape, Mescape => symba_plA%helio%swiftest%Mescape, & Euntracked => symba_plA%helio%swiftest%Euntracked, mass => symba_plA%helio%swiftest%mass, dMcb => symba_plA%helio%swiftest%dMcb, & - Mcb_initial => symba_plA%helio%swiftest%Mcb_initial) + Mcb_initial => symba_plA%helio%swiftest%Mcb_initial, Eorbit_orig => param%Eorbit_orig, Mtot_orig => param%Mtot_orig , & + Ltot_orig => param%Ltot_orig(:), Lmag_orig => param%Lmag_orig, Lorbit_orig => param%Lorbit_orig(:), Lspin_orig => param%Lspin_orig(:), & + lfirst => param%lfirstenergy) if (lfirst) then if (param%out_stat == "OLD") then open(unit = egyiu, file = energy_file, form = "formatted", status = "old", action = "write", position = "append") diff --git a/src/modules/io.f90 b/src/modules/io.f90 index 469d37cda..79ab37ee1 100644 --- a/src/modules/io.f90 +++ b/src/modules/io.f90 @@ -37,7 +37,7 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, type(symba_pl), intent(inout) :: symba_plA !! Swiftest planet data structure integer(I4B), intent(in) :: npl !! Number of massive bodies real(DP), intent(in) :: j2rp2, j4rp4 !! Central body oblateness terms - type(user_input_parameters), intent(in) :: param !! Input colleciton of user-defined parameters + type(user_input_parameters), intent(inout) :: param !! Input colleciton of user-defined parameters logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen end subroutine io_conservation_report diff --git a/src/modules/user.f90 b/src/modules/user.f90 index 7cba21f51..2f7bf92e9 100644 --- a/src/modules/user.f90 +++ b/src/modules/user.f90 @@ -49,6 +49,14 @@ module user logical :: ltides = .false. !! Include tidal dissipation logical :: lringmoons = .false. !! Turn on the ringmoons code logical :: lenergy = .false. !! Track the total energy of the system + logical :: lfirstenergy = .true. + real(DP) :: Eorbit_orig = 0.0_DP + real(DP) :: Mtot_orig = 0.0_DP + real(DP) :: Lmag_orig = 0.0_DP + real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP + real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP + real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP + logical :: lfirstkick = .true. !! Initiate the first kick in a symplectic step ! Future features not implemented or in development logical :: lgr = .false. !! Turn on GR diff --git a/src/symba/symba_step_eucl.f90 b/src/symba/symba_step_eucl.f90 index d0971004e..0634820d1 100644 --- a/src/symba/symba_step_eucl.f90 +++ b/src/symba/symba_step_eucl.f90 @@ -87,73 +87,74 @@ SUBROUTINE symba_step_eucl(t,dt,param,npl, ntp,symba_plA, symba_tpA, & INTEGER(I4B) :: irec, nplm, ipl, ipl1, ipl2 INTEGER(I8B) :: i, k, counter INTEGER(I8B), DIMENSION(:), ALLOCATABLE :: pltp_encounters_indices - LOGICAL, SAVE :: lfirst = .true. LOGICAL(LGT), ALLOCATABLE, DIMENSION(:) :: pltp_lencounters LOGICAL(lgt), ALLOCATABLE, DIMENSION(:) :: plpl_lvdotr, pltp_lvdotr ! Executable code - call symba_step_reset(npl, symba_plA, symba_tpA, plplenc_list, pltpenc_list, mergeadd_list, mergesub_list) - nplplenc = 0 - npltpenc = 0 - irec = 0 - nplm = count(symba_plA%helio%swiftest%mass(1:npl) > param%mtiny) + associate(lfirst => param%lfirstkick) + call symba_step_reset(npl, symba_plA, symba_tpA, plplenc_list, pltpenc_list, mergeadd_list, mergesub_list) + nplplenc = 0 + npltpenc = 0 + irec = 0 + nplm = count(symba_plA%helio%swiftest%mass(1:npl) > param%mtiny) - CALL symba_chk_eucl(npl, irec, symba_plA, dt, plplenc_list, nplplenc) - - if(ntp>0)then - allocate(pltp_lencounters(symba_tpA%helio%swiftest%num_pltp_comparisons)) - allocate(pltp_lvdotr(symba_tpA%helio%swiftest%num_pltp_comparisons)) - pltp_lencounters = .false. - pltp_lvdotr = .false. + CALL symba_chk_eucl(npl, irec, symba_plA, dt, plplenc_list, nplplenc) + + if(ntp>0)then + allocate(pltp_lencounters(symba_tpA%helio%swiftest%num_pltp_comparisons)) + allocate(pltp_lvdotr(symba_tpA%helio%swiftest%num_pltp_comparisons)) + pltp_lencounters = .false. + pltp_lvdotr = .false. - CALL symba_chk_eucl_pltp(symba_plA, symba_tpA, dt, pltp_lencounters, pltp_lvdotr, npltpenc) - - ! npltpenc = count(pltp_encounters > 0) - ! print *,'step npltpenc: ',npltpenc - if(npltpenc>0)then + CALL symba_chk_eucl_pltp(symba_plA, symba_tpA, dt, pltp_lencounters, pltp_lvdotr, npltpenc) + + ! npltpenc = count(pltp_encounters > 0) + ! print *,'step npltpenc: ',npltpenc + if(npltpenc>0)then - allocate(pltp_encounters_indices(npltpenc)) + allocate(pltp_encounters_indices(npltpenc)) - counter = 1 - do k = 1,symba_tpA%helio%swiftest%num_pltp_comparisons - if(pltp_lencounters(k))then - pltp_encounters_indices(counter) = k - counter = counter + 1 - endif - enddo + counter = 1 + do k = 1,symba_tpA%helio%swiftest%num_pltp_comparisons + if(pltp_lencounters(k))then + pltp_encounters_indices(counter) = k + counter = counter + 1 + endif + enddo - symba_plA%ntpenc(symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))) = symba_plA%ntpenc(symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))) + 1 - symba_tpA%nplenc(symba_tpA%helio%swiftest%k_pltp(2,pltp_encounters_indices(:))) = symba_tpA%nplenc(symba_tpA%helio%swiftest%k_pltp(2,pltp_encounters_indices(:))) + 1 + symba_plA%ntpenc(symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))) = symba_plA%ntpenc(symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))) + 1 + symba_tpA%nplenc(symba_tpA%helio%swiftest%k_pltp(2,pltp_encounters_indices(:))) = symba_tpA%nplenc(symba_tpA%helio%swiftest%k_pltp(2,pltp_encounters_indices(:))) + 1 - pltpenc_list%status(1:npltpenc) = ACTIVE - pltpenc_list%lvdotr(1:npltpenc) = pltp_lvdotr(pltp_encounters_indices(:)) - pltpenc_list%level(1:npltpenc) = 0 - pltpenc_list%indexpl(1:npltpenc) = symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:)) - pltpenc_list%indextp(1:npltpenc) = symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:)) + pltpenc_list%status(1:npltpenc) = ACTIVE + pltpenc_list%lvdotr(1:npltpenc) = pltp_lvdotr(pltp_encounters_indices(:)) + pltpenc_list%level(1:npltpenc) = 0 + pltpenc_list%indexpl(1:npltpenc) = symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:)) + pltpenc_list%indextp(1:npltpenc) = symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:)) - deallocate(pltp_encounters_indices) - endif + deallocate(pltp_encounters_indices) + endif - deallocate(pltp_lencounters, pltp_lvdotr) - endif - -! END OF THINGS THAT NEED TO BE CHANGED IN THE TREE + deallocate(pltp_lencounters, pltp_lvdotr) + endif + + ! END OF THINGS THAT NEED TO BE CHANGED IN THE TREE - ! flag to see if there was an encounter - lencounter = ((nplplenc > 0) .OR. (npltpenc > 0)) + ! flag to see if there was an encounter + lencounter = ((nplplenc > 0) .OR. (npltpenc > 0)) - IF (lencounter) THEN ! if there was an encounter, we need to enter symba_step_interp to see if we need recursion - CALL symba_step_interp_eucl(t, npl, nplm, ntp, symba_plA, symba_tpA, dt, nplplenc, npltpenc, & - plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list, param) - lfirst = .TRUE. - ELSE ! otherwise we can just advance the particles - CALL symba_step_helio(lfirst, param%lextra_force, t, npl, nplm, ntp,& - symba_plA%helio, symba_tpA%helio, & - param%j2rp2, param%j4rp4, dt) - END IF + IF (lencounter) THEN ! if there was an encounter, we need to enter symba_step_interp to see if we need recursion + CALL symba_step_interp_eucl(t, npl, nplm, ntp, symba_plA, symba_tpA, dt, nplplenc, npltpenc, & + plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list, param) + lfirst = .TRUE. + ELSE ! otherwise we can just advance the particles + CALL symba_step_helio(lfirst, param%lextra_force, t, npl, nplm, ntp,& + symba_plA%helio, symba_tpA%helio, & + param%j2rp2, param%j4rp4, dt) + END IF - RETURN + RETURN + end associate END SUBROUTINE symba_step_eucl !********************************************************************************************************************************** diff --git a/src/user/user_udio_reader.f90 b/src/user/user_udio_reader.f90 index 027eadefb..ea10222e1 100644 --- a/src/user/user_udio_reader.f90 +++ b/src/user/user_udio_reader.f90 @@ -179,6 +179,37 @@ module subroutine user_udio_reader(param, unit, iotype, v_list, iostat, iomsg) param%seed(nseeds_from_file+1:nseeds) = [(param%seed(1) - param%seed(nseeds_from_file) + i, i=nseeds_from_file+1, nseeds)] end if seed_set = .true. + case ("FIRSTKICK") + call util_toupper(param_value) + if (param_value == "NO" .or. param_value == 'F') param%lfirstkick = .false. + case ("FIRSTENERGY") + call util_toupper(param_value) + if (param_value == "NO" .or. param_value == 'F') param%lfirstenergy = .false. + case("EORBIT_ORIG") + read(param_value, *) param%Eorbit_orig + case("MTOT_ORIG") + read(param_value, *) param%Mtot_orig + case("LTOT_ORIG") + read(param_value, *) param%Ltot_orig(1) + do i = 2, NDIM + ifirst = ilast + 1 + param_value = user_get_token(line, ifirst, ilast, iostat) + read(param_value, *) param%Ltot_orig(i) + end do + case("LORBIT_ORIG") + read(param_value, *) param%Lorbit_orig(1) + do i = 2, NDIM + ifirst = ilast + 1 + param_value = user_get_token(line, ifirst, ilast, iostat) + read(param_value, *) param%Lorbit_orig(i) + end do + case("LSPIN_ORIG") + read(param_value, *) param%Lspin_orig(1) + do i = 2, NDIM + ifirst = ilast + 1 + param_value = user_get_token(line, ifirst, ilast, iostat) + read(param_value, *) param%Lspin_orig(i) + end do case default write(iomsg,*) "Unknown parameter -> ",param_name iostat = -1 diff --git a/src/user/user_udio_writer.f90 b/src/user/user_udio_writer.f90 index 37b5bd0b4..d138ab584 100644 --- a/src/user/user_udio_writer.f90 +++ b/src/user/user_udio_writer.f90 @@ -21,9 +21,10 @@ module subroutine user_udio_writer(param, unit, iotype, v_list, iostat, iomsg) ! Internals character(*),parameter :: Ifmt = '(I0)' !! Format label for integer values character(*),parameter :: Rfmt = '(ES25.17)' !! Format label for real values + character(*),parameter :: Rarrfmt = '(3(ES25.17,1X))' !! Format label for real values character(*),parameter :: Lfmt = '(L1)' !! Format label for logical values character(len=*), parameter :: Afmt = '(A25,1X,64(:,A25,1X))' - character(25) :: param_name, param_value + character(256) :: param_name, param_value type character_array character(25) :: value end type character_array @@ -106,6 +107,15 @@ module subroutine user_udio_writer(param, unit, iotype, v_list, iostat, iomsg) if (param%lringmoons) then write(param_name, Afmt) "RING_OUTFILE"; write(param_value, Afmt) trim(adjustl(param%ring_outfile)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) end if + if (param%lenergy) then + write(param_name, Afmt) "FIRSTENERGY"; write(param_value, Lfmt) param%lfirstenergy; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "EORBIT_ORIG"; write(param_value, Rfmt) param%Eorbit_orig; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "MTOT_ORIG"; write(param_value, Rfmt) param%Mtot_orig; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(unit, '("LTOT_ORIG ",3(1X,ES25.17))') param%Ltot_orig(:) + write(unit, '("LORBIT_ORIG",3(1X,ES25.17))') param%Lorbit_orig(:) + write(unit, '("LSPIN_ORIG ",3(1X,ES25.17))') param%Lspin_orig(:) + end if + write(param_name, Afmt) "FIRSTKICK"; write(param_value, Lfmt) param%lfirstkick; write(unit, Afmt) adjustl(param_name), adjustl(param_value) iostat = 0