diff --git a/examples/Basic_Simulation/initial_conditions.py b/examples/Basic_Simulation/initial_conditions.py index 95fdb608e..78900ce1f 100644 --- a/examples/Basic_Simulation/initial_conditions.py +++ b/examples/Basic_Simulation/initial_conditions.py @@ -15,14 +15,7 @@ Returns ------- -param.in : ASCII text file - Swiftest parameter input file. -pl.in : ASCII text file - Swiftest massive body input file. -tp.in : ASCII text file - Swiftest test particle input file. -cb.in : ASCII text file - Swiftest central body input file. +Updates sim.data with the simulation data """ import swiftest @@ -72,5 +65,5 @@ sim.add_body(name_tp, a_tp, e_tp, inc_tp, capom_tp, omega_tp, capm_tp) -# Save everything to a set of initial conditions files +# Run the simulation sim.run() diff --git a/examples/Basic_Simulation/param.in b/examples/Basic_Simulation/param.in index 014bf1fb8..0ee870562 100644 --- a/examples/Basic_Simulation/param.in +++ b/examples/Basic_Simulation/param.in @@ -22,7 +22,6 @@ MU2KG 1.988409870698051e+30 TU2S 31557600.0 DU2M 149597870700.0 GMTINY 9.869231602224408e-07 -MIN_GMFRAG 9.869231602224408e-10 RESTART NO CHK_CLOSE YES GR YES @@ -35,3 +34,4 @@ RHILL_PRESENT NO INTERACTION_LOOPS TRIANGULAR ENCOUNTER_CHECK TRIANGULAR TIDES NO +MIN_GMFRAG 9.869231602224408e-10 diff --git a/examples/Basic_Simulation/run_from_file.py b/examples/Basic_Simulation/run_from_file.py index cede6e2ea..9c477410a 100644 --- a/examples/Basic_Simulation/run_from_file.py +++ b/examples/Basic_Simulation/run_from_file.py @@ -1,3 +1,3 @@ import swiftest sim = swiftest.Simulation() -sim.run() \ No newline at end of file +sim.run(tstop=20.0) \ No newline at end of file diff --git a/examples/Basic_Simulation/run_simulation.ipynb b/examples/Basic_Simulation/run_simulation.ipynb new file mode 100644 index 000000000..8f2ab51b3 --- /dev/null +++ b/examples/Basic_Simulation/run_simulation.ipynb @@ -0,0 +1,159 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "86c845ce-1801-46ca-8a8a-1cabb266e6a6", + "metadata": {}, + "outputs": [], + "source": [ + "import swiftest\n", + "import xarray as xr\n", + "import numpy as np\n", + "import os" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d716c371-8eb4-4fc1-82af-8b5c444c831e", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/param.in\n" + ] + } + ], + "source": [ + "sim = swiftest.Simulation()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "ec7452d6-4c9b-4df3-acc0-b11c32264b91", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "tstop 10.0 y\n", + "Writing parameter inputs to file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/param.in\n", + "Running a Swiftest symba run from tstart=0.0 y to tstop=10.0 y\n", + "\u001b]2;cd /home/daminton/git_debug/swiftest/examples/Basic_Simulation\u0007\u001b]1;\u0007\u001b]2;/home/daminton/git_debug/swiftest/bin/swiftest_driver symba \u0007\u001b]1;\u0007 Parameter input file is /home/daminton/git_debug/swiftest/examples/Basic_Simulation/param.in\n", + " \n", + " Warning! NPLM variable not set in input file. Calculating.\n", + " *************** Main Loop *************** \n", + "Time = 1.00000E+00; fraction done = 0.100; Number of active plm, pl, tp = 13, 14, 10\n", + "Integration steps: Total wall time: 6.30684E-01; Interval wall time: 5.29890E-01;Interval wall time/step: 2.70141E-03\n", + "Time = 2.00000E+00; fraction done = 0.200; Number of active plm, pl, tp = 13, 14, 10\n", + "Integration steps: Total wall time: 1.66455E+00; Interval wall time: 5.27720E-01;Interval wall time/step: 2.99766E-03\n", + "Time = 3.00000E+00; fraction done = 0.300; Number of active plm, pl, tp = 13, 14, 10\n", + "Integration steps: Total wall time: 2.64051E+00; Interval wall time: 5.20805E-01;Interval wall time/step: 2.88832E-03\n", + "Time = 4.00000E+00; fraction done = 0.400; Number of active plm, pl, tp = 13, 14, 10\n", + "Integration steps: Total wall time: 3.60585E+00; Interval wall time: 5.24579E-01;Interval wall time/step: 2.82311E-03\n", + "Time = 5.00000E+00; fraction done = 0.500; Number of active plm, pl, tp = 13, 14, 10\n", + "Integration steps: Total wall time: 4.58823E+00; Interval wall time: 5.37595E-01;Interval wall time/step: 2.96439E-03\n", + "Time = 6.00000E+00; fraction done = 0.600; Number of active plm, pl, tp = 13, 14, 10\n", + "Integration steps: Total wall time: 5.55600E+00; Interval wall time: 5.27663E-01;Interval wall time/step: 2.83397E-03\n", + "Time = 7.00000E+00; fraction done = 0.700; Number of active plm, pl, tp = 13, 14, 10\n", + "Integration steps: Total wall time: 6.64194E+00; Interval wall time: 5.83431E-01;Interval wall time/step: 3.11589E-03\n", + "Time = 8.00000E+00; fraction done = 0.800; Number of active plm, pl, tp = 13, 14, 10\n", + "Integration steps: Total wall time: 7.63044E+00; Interval wall time: 5.52408E-01;Interval wall time/step: 2.95843E-03\n", + "Time = 9.00000E+00; fraction done = 0.900; Number of active plm, pl, tp = 13, 14, 10\n", + "Integration steps: Total wall time: 8.62339E+00; Interval wall time: 5.46944E-01;Interval wall time/step: 3.01578E-03\n", + "Time = 1.00000E+01; fraction done = 1.000; Number of active plm, pl, tp = 13, 14, 10\n", + "Integration steps: Total wall time: 9.67297E+00; Interval wall time: 6.00750E-01;Interval wall time/step: 3.20243E-03\n", + "\n", + "Normal termination of Swiftest (version 1.0)\n", + "------------------------------------------------\n", + "\n", + "Creating Dataset from NetCDF file\n", + "Successfully converted 11 output frames.\n", + "Swiftest simulation data stored as xarray DataSet .ds\n" + ] + } + ], + "source": [ + "sim.run(tstop=10.0)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "5b0a57a6-dbd5-4d34-8e6e-91fc8ad8789f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ,\n", + " ]" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAACEU0lEQVR4nOzdfXyN9f/A8de12c7uN5vZDWez2TDGRiKjkHszFaEIaygJCSWiKFnubyIhG2VuSlTfiuKbkbvITS18VWyNH2uRtsy2s51z/f6YnXZsY9M456z389H16Fyfz+f6XO9zOTvX+3yuO0VVVRUhhBBCCCtlY+4AhBBCCCH+CUlmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBBCCGHVJJkRQgghhFWrYe4A7jSDwcCFCxdwdXVFURRzhyOEEEKIClBVlb/++gt/f39sbG4+9lLtk5kLFy6g1WrNHYYQQgghbsO5c+eoW7fuTdtU+2TG1dUVKNoYbm5uZo5GCCGEEBWRnZ2NVqs17sdvptonM8WHltzc3CSZEUIIIaxMRU4RkROAhRBCCGHVJJkRQgghhFWTZEYIIYQQVq3anzMjhCXT6/UUFBSYOwzxL2dnZ4etra25wxDitkkyI4QZqKpKRkYGf/75p7lDEQIADw8PfH195X5cwipZTDITHx/PlClTeO6551i0aBFQ9IU/Y8YMVq5cyZUrV2jdujXLli2jSZMm5g1WiH+oOJGpXbs2Tk5OsgMRZqOqKteuXSMzMxMAPz8/M0ckROVZRDJz+PBhVq5cSbNmzUzK58yZw4IFC1izZg0NGjRg5syZdOnShdOnT1founMhLJFerzcmMl5eXuYORwgcHR0ByMzMpHbt2nLISVgds58AfPXqVQYNGsSqVauoWbOmsVxVVRYtWsTLL79Mnz59CA8PZ+3atVy7do3169ebMWIh/pnic2ScnJzMHIkQfyv+PMo5XMIamT2ZefbZZ4mOjqZz584m5ampqWRkZNC1a1djmUajoX379uzfv7/c/vLz88nOzjaZhLBEcmhJWBL5PAprZtbDTBs3buTo0aMcPny4VF1GRgYAPj4+JuU+Pj78+uuv5fYZHx/PjBkzqjZQIYQQQlgss43MnDt3jueee45169bh4OBQbrsbfy2oqnrTXxCTJ08mKyvLOJ07d67KYhZCCCGE5THbyMyRI0fIzMzknnvuMZbp9Xr27NnD0qVLOX36NFA0QlPy7PrMzMxSozUlaTQaNBrNnQtcCCGEEBbFbMlMp06dSElJMSl78sknadSoEZMmTSI4OBhfX1927NhB8+bNAdDpdOzevZvZs2ebI+S7TlVVoGhSVcP1/xeXGa6/NlxvazBthwrX64teG0r0WXLESwFFQUEpen3DfFG7vyeT5cqsL66zMamT4/FCCCHuFLMlM66uroSHh5uUOTs74+XlZSwfN24cs2bNIjQ0lNDQUGbNmoWTkxMDBw40R8gmzv/fetLTV/2dXKjXk4brSYRKyeTi72Tj7ySkYklK9XOrJKh43ub66+LEyxZFsUFRbP+esAXFBkWpUaLu+jw2ULJtmXUlluXvdrdfZ2MSm2JTAxvFDkWpgY2N3fVl7CgsrIHBYIten4deT6ntIUmgEEJUjkXcZ6Y8L774Irm5uYwaNcp407yvvvrKIu4xU1iQTW5uurnDuAXTxODvnaRaNGhDyaRKLTF/JxWv74bSO71aC2Jj44+H+6vk5iro9RVIVm4yclZUXV5d6XnTtpRfd5ME62Z9mvZx8ytkOnToQLNmzXBwcODdd9/F3t6ekSNHMn36dAAWLFhAYmIiZ8+exdPTk5iYGObMmYOLiwsAa9asYdy4caxbt44JEyZw7tw5evbsydq1a9m8eTOvvvoqWVlZPPHEEyxatMh47xSdTsfUqVNJSkrizz//JDw8nNmzZ9OhQ4db/1sIISySRSUzycnJJvOKojB9+nTjl5sl8fV9CI+a917/Ure54dDM9VEFxebvL31FKRoRKPGF//foQ4lDMhVtV2aSUlxv849/0RePDpkmObdOgsoadeL6mBUl64yHwUovZ7psyVj0f08YUNXConYly6/XoeoxqPrrdSWW5e92xjr0qIbiuuL+Cv+uU8tZF6Z9c0NcqlqIwVBwfb7A+Bq1ZolRI+WG7VHmP0bxFrw+f+dTzqpTfiKk1+eydm0io0fHsWvXhxw6dIynn36Re+4JoVOnB9Drs5g7dyqBgVrSfj3PuOemMGHCNRYvfhNFUSgoyObatWssWjSf999/h7/+usqAAbE8/HAMHh4efPLJJlJT0xgwYAht2txD//6PAgpDh8bx66/pJCWtxd/fn48//oTu3bvz/fdHCA1tUHbcyKXLwrxUVb3+XVLy+67w7+8cgx6uf78Z1MKiowVltS81X047KrY8qgGDWoiXZztq1XrQbNtHUdXq/Zs4Ozsbd3d3srKycHNzM3c4QpCXl0dqaipBQUGlruQzTQbLSxpvnlSWbltWEslN6m623K36rLjo6Dj0ej3bt681lnXsOJAHHmjFjBnjSrXfuvUrxo+fSWrqHgCSkj5h1KhpHDv2OcHBWgDGjXudTZv+w88/J+PiUnQTuD59RhIQUIdFi6Zx9uw5WrToxalTO/Dzq23su3fvEdxzTzivvvrcTSK+cTSrqOzGpOfGkbNSo2XGMkr0U8Zr47rKqTe+LNnq9tvk5RWQnn4ed/ffsLPXXz8ca4Ni8iPq+g+M4nTa+DksLisuL/EZK1V2wzIl2t6s//KWMUn0jespbmu4vlP/+4eJ8QcJaukfNer19iV+qFD8g8XYV/GPlus/Yq73c2Nb47yxraHED56S9fq/6411NyYZhSW2g2UKDHyGkPoTq7TPyuy/LWpkRoh/u793GmBtAwEVTqiu74hsbBwID2+Eo2Ogsc7fX8uVK3k4OPizK/kb5sxeyP/+9xPZ2X9RWKgnLy+PggJHnJwdsbV1wsnJkYYNw43Jna+vH4GBWtw9vI07PB8fHy5d+hMbGw0pKT+jqir33NObkjuH/PwCPL08KLnDLuMdXv/PtL66/BzU6VQKCq7wy5l4DIYL5g5HVIKi2JU4l6/keYI1riejNVBsbAFbbIrLjecA3jjd2McNfdnU+Pu8wBLtPTxamXUbSDIjhKgSlU3EFMUWjcYZO7u/f3HZ2toDdly8eJWHH3qckSNH8sYbs/H09GTv3r0MGzYMW9taODp4YG/viZ2dPY6OAcbl7e1rYm/vhLNTsLHMzs4NGxsDLi4NsLc/hq2tLUeOHC31/CEXFxfc3HxRbxgpUEv80i9rvmSCdvP5G/oprjZJjtQb/l88EnFjfRltbtrfrdep1xdiY6PB3a0FKIGYXiFZdMi56Dz9v0eliv+vlBhxMh35uT66pJTR1mSZstuW3X/JuoosU+LiAEqerF80X7q+eETq74sOjCf+F49WGXfm1+uxKdG2ZL1StCx/Xxzw92hX8bpuqEcpkVAUX9xQRpJS4sIGIcmMEMICfffddxQWFjJ//nxsbIq+rD/44IN/3G/z5s3R6/VkZmZy//33l9lGuWHHaG0jZLfLxiYPe/sCgoJm3/RGpkJYIknphBAWp379+hQWFvLWW29x9uxZ3n//fd55551/3G+DBg0YNGgQQ4YMYcuWLaSmpnL48GFmz57NF198UQWRCyHMQZIZIYTFiYyMZMGCBcyePZvw8HCSkpKIj4+vkr4TExMZMmQIEyZMoGHDhvTu3Ztvv/0WrVZbJf0LIe4+uZpJiLvsZlczCWEu8rkUlqYy+28ZmRFCCCGEVZNkRgghhBBWTZIZIYQQQlg1SWaEEEIIYdUkmRFCCCGEVZNkRgghhBBWTZIZIYQQQlg1SWaEEEIIYdUkmRFCCCGEVZNkRghRYbGxsSiKwsiRI0vVjRo1CkVRiI2NvfuBCSH+1SSZEUJUilarZePGjeTm5hrL8vLy2LBhAwEBAbfdr6qqFBYWVkWIJvR6PQaDocr7FUJYDklmhBCV0qJFCwICAtiyZYuxbMuWLWi1Wpo3b24sU1WVOXPmEBwcjKOjIxEREWzevNlYn5ycjKIofPnll7Rs2RKNRsM333yDwWBg9uzZhISEoNFoCAgI4I033jBZ5s8//zT2c/z4cRRFIS0tDYA1a9bg4eHBZ599RuPGjY392tnZkZGRYfJeJkyYwAMPPHAHtpIQ4m6qYe4AhBBFO/7cAr1Z1u1oZ4uiKJVa5sknnyQxMZFBgwYBkJCQQFxcHMnJycY2U6dOZcuWLSxfvpzQ0FD27NnDE088gbe3N+3btze2e/HFF5k3bx7BwcF4eHgwefJkVq1axcKFC2nXrh0XL17kf//7X6Xiu3btGvHx8bz77rt4eXlRt25dgoODef/993nhhRcAKCwsZN26dbz55puV6lsIYXkkmRHCAuQW6Gn8ypdmWffJ17rhZF+5r4LBgwczefJk0tLSUBSFffv2sXHjRmMyk5OTw4IFC/j6669p06YNAMHBwezdu5cVK1aYJDOvvfYaXbp0AeCvv/5i8eLFLF26lKFDhwJQv3592rVrV6n4CgoKePvtt4mIiDCWDRs2jMTERGMy8/nnn3Pt2jX69+9fqb6FEJZHkhkhRKXVqlWL6Oho1q5di6qqREdHU6tWLWP9yZMnycvLMyYpxXQ6ncmhKICWLVsaX586dYr8/Hw6der0j+Kzt7enWbNmJmWxsbFMnTqVgwcPct9995GQkED//v1xdnb+R+sSQpifJDNCWABHO1tOvtbNbOu+HXFxcYwePRqAZcuWmdQVn3D7+eefU6dOHZM6jUZjMl8ymXB0dLzpOm1sik7zU1XVWFZQUFCqnaOjY6lDZ7Vr1yYmJobExESCg4P54osvTA6LCSGslyQzQlgARVEqfajH3Lp3745OpwOgWzfTRKz4xNv09HSTQ0q3EhoaiqOjI//9738ZPnx4qXpvb28ALl68SM2aNYGiE4Aravjw4Tz22GPUrVuX+vXr07Zt2wovK4SwXNb17SmEsBi2tracOnXK+LokV1dXJk6cyPPPP4/BYKBdu3ZkZ2ezf/9+XFxcjOfD3MjBwYFJkybx4osvYm9vT9u2bfn99985ceIEw4YNIyQkBK1Wy/Tp05k5cyY///wz8+fPr3DM3bp1w93dnZkzZ/Laa6/d/psXQlgUSWaEELfNzc2t3LrXX3+d2rVrEx8fz9mzZ/Hw8KBFixZMmTLlpn1OmzaNGjVq8Morr3DhwgX8/PyMN+mzs7Njw4YNPPPMM0RERHDvvfcyc+ZM+vXrV6F4bWxsiI2NZdasWQwZMqTib1QIYdEUteTB52ooOzsbd3d3srKybvrFK8TdkpeXR2pqKkFBQTg4OJg7nH+dESNG8Ntvv/Hpp5+aOxSLIp9LYWkqs/+WkRkhxL9CVlYWhw8fJikpiU8++cTc4QghqpAkM0KIf4WHHnqIQ4cO8fTTT5e6ZFwIYd0kmRFC/CvIZdhCVF/ybCYhhBBCWDVJZoQQQghh1cyazCxfvpxmzZrh5uaGm5sbbdq0Ydu2bcb62NhYFEUxme677z4zRiyEEEIIS2PWc2bq1q3Lm2++SUhICABr167loYce4tixYzRp0gQoustoYmKicRl7e3uzxCqEEEIIy2TWZCYmJsZk/o033mD58uUcPHjQmMxoNBp8fX3NEZ4QQgghrIDFnDOj1+vZuHEjOTk5tGnTxlienJxM7dq1adCgASNGjCAzM/Om/eTn55OdnW0yCSGEEKL6Mnsyk5KSgouLCxqNhpEjR7J161YaN24MQI8ePUhKSuLrr79m/vz5HD58mAcffJD8/Pxy+4uPj8fd3d04abXau/VWhBBCCGEGZn+cgU6nIz09nT///JOPPvqId999l927dxsTmpIuXrxIYGAgGzdupE+fPmX2l5+fb5LsZGdno9Vq5XEGwmJY623jY2JiyM3NZefOnaXqDhw4QFRUFEeOHKFFixZmiE78U9b6uRTVl1U9zsDe3t54AnDLli05fPgwixcvZsWKFaXa+vn5ERgYyM8//1xufxqNBo1Gc8fiFeLfatiwYfTp04dff/2VwMBAk7qEhAQiIyMlkRFCmIXZDzPdSFXVcg8jXb58mXPnzuHn53eXoxJC9OrVi9q1a7NmzRqT8mvXrrFp0yaGDRvG/v37eeCBB3B0dESr1TJ27FhycnKMbevVq8esWbOIi4vD1dWVgIAAVq5caaxPTk5GURT+/PNPY9nx48dRFIW0tDQAfv31V2JiYqhZsybOzs40adKEL7744k6+dSGEhTNrMjNlyhS++eYb0tLSSElJ4eWXXyY5OZlBgwZx9epVJk6cyIEDB0hLSyM5OZmYmBhq1arFI488Ys6whah6qgq6HPNMFTzSXKNGDYYMGcKaNWsoeXT6ww8/RKfTERERQbdu3ejTpw8//PADmzZtYu/evYwePdqkn/nz59OyZUuOHTvGqFGjeOaZZ/jf//5X4U317LPPkp+fz549e0hJSWH27Nm4uLhUeHkhRPVj1sNMv/32G4MHD+bixYu4u7vTrFkztm/fTpcuXcjNzSUlJYX33nuPP//8Ez8/Pzp27MimTZtwdXU1Z9hCVL2CazDL3zzrnnIB7J0r1DQuLo65c+eSnJxMx44dgaJDTH369GHVqlUMHDiQcePGARAaGsqSJUto3749y5cvN56H0bNnT0aNGgXApEmTWLhwIcnJyTRq1KhCMaSnp9O3b1+aNm0KQHBwcGXerRCiGjJrMrN69epy6xwdHfnyyy/vYjRCiFtp1KgRUVFRJCQk0LFjR86cOcM333zDV199xXPPPccvv/xCUlKSsb2qqhgMBlJTUwkLCwOgWbNmxnpFUfD19b3lLRdKGjt2LM888wxfffUVnTt3pm/fviZ9CiH+fcx+ArAQArBzKhohMde6K2HYsGGMHj2aZcuWkZiYSGBgIJ06dcJgMPD0008zduzYUssEBAT8vTo7O5M6RVEwGAwA2NgUHfkueRiroKDApP3w4cPp1q0bn3/+OV999RXx8fHMnz+fMWPGVOp9CCGqD0lmhLAEilLhQz3m1r9/f5577jnWr1/P2rVrGTFiBIqi0KJFC06cOGG8OvF2eHt7A0W3YahZsyZQdALwjbRaLSNHjmTkyJFMnjyZVatWSTIjxL+YxV3NJISwbC4uLgwYMIApU6Zw4cIFYmNjgaLzXw4cOMCzzz7L8ePH+fnnn/n0008rlWSEhISg1WqZPn06P/30E59//jnz5883aTNu3Di+/PJLUlNTOXr0KF9//bXxEJYQ4t9JkhkhRKUNGzaMK1eu0LlzZ+MhpGbNmrF7925+/vln7r//fpo3b860adMqdSsFOzs7NmzYwP/+9z8iIiKYPXs2M2fONGmj1+t59tlnCQsLo3v37jRs2JC33367St+fEMK6mP0OwHdaZe4gKMTdIHdaFZZIPpfC0lRm/y0jM0IIIYSwapLMCCGEEMKqSTIjhBBCCKsmyYwQQgghrJokM0IIIYSwapLMCCGEEMKqSTIjhBBCCKsmyYwQQgghrJokM0IIIYSwapLMCCEsUlpaGoqilPmgSSGEKEmSGSFEhcXGxqIoSqmpe/fu/7jfhx9+uGqCFEL869QwdwBCCOvSvXt3EhMTTco0Gs1t9aXX61EUpSrCEkL8i8nIjBCiUjQaDb6+viZTzZo1AViwYAFNmzbF2dkZrVbLqFGjuHr1qnHZNWvW4OHhwWeffUbjxo3RaDQ8+eSTrF27lk8++cQ40pOcnGxc5uzZs3Ts2BEnJyciIiI4cODA3X7LQggLJyMzQlgAVVXJLcw1y7odazhW2eiIjY0NS5YsoV69eqSmpjJq1ChefPFF3n77bWOba9euER8fz7vvvouXlxe+vr7k5eWRnZ1tHPHx9PTkwoULALz88svMmzeP0NBQXn75ZR5//HF++eUXatSQry8hRBH5NhDCAuQW5tJ6fWuzrPvbgd/iZOdU4fafffYZLi4uJmWTJk1i2rRpjBs3zlgWFBTE66+/zjPPPGOSzBQUFPD2228TERFhLHN0dCQ/Px9fX99S65s4cSLR0dEAzJgxgyZNmvDLL7/QqFGjCscshKjeJJkRQlRKx44dWb58uUmZp6cnALt27WLWrFmcPHmS7OxsCgsLycvLIycnB2dnZwDs7e1p1qxZhddXsq2fnx8AmZmZkswIIYwkmRHCAjjWcOTbgd+abd2V4ezsTEhISKnyX3/9lZ49ezJy5Ehef/11PD092bt3L8OGDaOgoODv9TlW7rCWnZ2d8XXxcgaDoVIxCyGqN0lmhLAAiqJU6lCPJfruu+8oLCxk/vz52NgUXVvwwQcfVGhZe3t79Hr9nQxPCFGNSTIjhKiU/Px8MjIyTMpq1KhB/fr1KSws5K233iImJoZ9+/bxzjvvVKjPevXq8eWXX3L69Gm8vLxwd3e/E6ELIaopuTRbCFEp27dvx8/Pz2Rq164dkZGRLFiwgNmzZxMeHk5SUhLx8fEV6nPEiBE0bNiQli1b4u3tzb59++7wuxBCVCeKqqqquYO4k7Kzs3F3dycrKws3NzdzhyMEeXl5pKamEhQUhIODg7nDEQKQz6WwPJXZf8vIjBBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqZk1mli9fTrNmzXBzc8PNzY02bdqwbds2Y72qqkyfPh1/f38cHR3p0KEDJ06cMGPEQgghhLA0Zr0DcN26dXnzzTeNz3lZu3YtDz30EMeOHaNJkybMmTOHBQsWsGbNGho0aMDMmTPp0qULp0+fxtXV1ZyhWy1VVUtN5ZVbUtuS8Zf3viypvJiiKMbnCRW/Ll4mLy/PZPnynld0s+cYVXaZipbfbL4yz1USQoi7weJumufp6cncuXOJi4vD39+fcePGMWnSJKDoNuo+Pj7Mnj2bp59+ukL93amb5qWlpXH27FlUVcVgMJj835LLhPm5uLjQtm1b6tSpQ40a1v1EkZslOVVdVxWvzdXWGhJAuWmesDSV2X9bzDepXq/nww8/JCcnhzZt2pCamkpGRgZdu3Y1ttFoNLRv3579+/eXm8zk5+eTn59vnM/Ozr4j8Z47d449e/bckb4tTVmjC7eaqrJtcX3JeG4Wa2XK7+Qy5Y04aTQaatSogZ2dXbnJzM1+Y9ytUaWKqMiomTD1T5Kyf7r8zV7rdDry8vI4cuSI8W/PxsYGGxsb4+uK/k1b0nQ729Aaks877cZR8eKysl5D0TaztbW9K7GVxezJTEpKCm3atCEvLw8XFxe2bt1K48aN2b9/PwA+Pj4m7X18fPj111/L7S8+Pp4ZM2bc0ZgB/P39adWqlckf/Y1/8JZYduNruHXSIapW8S/gmjVrWuQv4JJJ2I3lcXFxvPfee4wYMYK3337bpG7MmDGsWLGCwYMH8+6775bbR1mvb1Z3O+3Kel3Rstvtq6IsNQEsLCwkLy+PAwcOcPXqVXOHY3FuJ6m83XbFbvbZv3G+qupuV7t27ejcuXOV9HU7zJ7MNGzYkOPHj/Pnn3/y0UcfMXToUHbv3m2sv/EfWFXVm+5gJ0+ezPjx443z2dnZaLXaKo+7fv361K9fv8r7FcLcbvbrVFEUtFotH3zwAYsXL8bR0REoStA2bdpEQEAANjY22Nvb39a6CwoKsLOzu/3gzayiSdKdSNb+6fI6nQ57e3saNmxIQUGByeHpkoetLWEqGc/dcqvEQpiX2ZMZe3t74wnALVu25PDhwyxevNh4nkxGRgZ+fn7G9pmZmaVGa0rSaDRoNJo7G7QQ/2ItWrTg7NmzbNmyhUGDBgGwZcsWtFotwcHBxnbbt29n5syZ/Pjjj9ja2tKmTRsWL15s/BGQlpZGUFAQmzZt4u233+bgwYMsX76cBx98kNGjR7N37150Oh316tVj7ty59OzZ0yzvtzKs+TBFXl4eTk5OdOnSxSJHDMtzq8SnZLsbl7tbdZXto6KjO7eav5ttzXmICSwgmbmRqqrk5+cTFBSEr68vO3bsoHnz5kDRL4fdu3cze/ZsM0cpRNVSVRU1N9cs61YcHSu9833yySdJTEw0JjMJCQnExcWRnJxsbJOTk8P48eNp2rQpOTk5vPLKKzzyyCMcP34cG5u/7woxadIk5s+fT2JiIhqNhqeeegqdTseePXtwdnbm5MmTuLi4VMl7FdWPHA4XYOZkZsqUKfTo0QOtVstff/3Fxo0bSU5OZvv27SiKwrhx45g1axahoaGEhoYya9YsnJycGDhwoDnDFqLKqbm5nG5xj1nW3fDoERQnp0otM3jwYCZPnkxaWhqKorBv3z7j32+xvn37miyzevVqateuzcmTJwkPDzeWjxs3jj59+hjn09PT6du3L02bNgUwGe0RQoiymDWZ+e233xg8eDAXL17E3d2dZs2asX37drp06QLAiy++SG5uLqNGjeLKlSu0bt2ar776Su4xI4SZ1apVi+joaNauXYuqqkRHR1OrVi2TNmfOnGHatGkcPHiQS5cuGW8NkJ6ebpLMtGzZ0mS5sWPH8swzz/DVV1/RuXNn+vbtS7Nmze78mxJCWC2zJjOrV6++ab2iKEyfPp3p06ffnYCEMBPF0ZGGR4+Ybd23Iy4ujtGjRwOwbNmyUvUxMTFotVpWrVqFv78/BoOB8PBwdDqdSTtnZ2eT+eHDh9OtWzc+//xzvvrqK+Lj45k/fz5jxoy5rTiFENWfxZ0zI8S/kaIolT7UY27du3c3JibdunUzqbt8+TKnTp1ixYoV3H///QDs3bu3wn1rtVpGjhzJyJEjmTx5MqtWrZJkRghRLklmhBC3xdbWllOnThlfl1SzZk28vLxYuXIlfn5+pKen89JLL1Wo33HjxtGjRw8aNGjAlStX+PrrrwkLC6vy+IUQ1YckM0KI21beLcZtbGzYuHEjY8eOJTw8nIYNG7JkyRI6dOhwyz71ej3PPvss58+fx83Nje7du7Nw4cIqjlwIUZ1Y3LOZqtqdejaTELdLnoEjLJF8LoWlqcz+2+amtUIIIYQQFk6SGSGEEEJYNUlmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGE2aSlpaEoCsePHzd3KEIIKybJjBCiwmJjY3n44YerrD+tVsvFixcJDw8HIDk5GUVR+PPPP6tsHUKI6k8eNCmEMBtbW1t8fX2rvF9VVdHr9dSoIV9xQvwbyMiMEBZAVVUK8vVmmW73WbP16tVj0aJFJmWRkZFMnz7dOK8oCsuXL6dHjx44OjoSFBTEhx9+aKwveZgpLS2Njh07AlCzZk0URSE2Nta4febMmUNwcDCOjo5ERESwefNmYz/FIzpffvklLVu2RKPR8M0339zW+xJCWB/52SKEBSjUGVj53G6zrPupxe2x09jesf6nTZvGm2++yeLFi3n//fd5/PHHCQ8PJywszKSdVqvlo48+om/fvpw+fRo3NzccHR0BmDp1Klu2bGH58uWEhoayZ88ennjiCby9vWnfvr2xjxdffJF58+YRHByMh4fHHXtPQgjLIsmMEOKO6tevH8OHDwfg9ddfZ8eOHbz11lu8/fbbJu1sbW3x9PQEoHbt2sZkJCcnhwULFvD111/Tpk0bAIKDg9m7dy8rVqwwSWZee+01unTpchfelRDCkkgyI4QFqGFvw1OL29+64R1a951UnICUnK/M1UsnT54kLy+vVJKi0+lo3ry5SVnLli1vO04hhPWSZEYIC6Aoyh091HMn2NjYlDrfpqCgoELLKopS4fUYDAYAPv/8c+rUqWNSp9FoTOadnZ0r3K8QovqQZEYIcVu8vb25ePGicT47O5vU1NRS7Q4ePMiQIUNM5m8cUSlmb28PgF6vN5Y1btwYjUZDenq6ySElIYQoJsmMEOK2PPjgg6xZs4aYmBhq1qzJtGnTsLUtPbr04Ycf0rJlS9q1a0dSUhKHDh1i9erVZfYZGBiIoih89tln9OzZE0dHR1xdXZk4cSLPP/88BoOBdu3akZ2dzf79+3FxcWHo0KF3+q0KISycJDNCiAozGAzGe7dMnjyZs2fP0qtXL9zd3Xn99dfLHJmZMWMGGzduZNSoUfj6+pKUlETjxo3L7L9OnTrMmDGDl156iSeffJIhQ4awZs0aXn/9dWrXrk18fDxnz57Fw8ODFi1aMGXKlDv6foUQ1kFRb/cmE1YiOzsbd3d3srKycHNzM3c4QpCXl0dqaipBQUE4ODiYO5xK6d69OyEhISxdurRC7RVFYevWrVV612BxZ1jz51JUT5XZf8tN84QQt3TlyhU+//xzkpOT6dy5s7nDEUIIE3KYSQhxS3FxcRw+fJgJEybw0EMPmTscIYQwIcmMEOKWtm7delvLVfOj2EIICyGHmYQQQghh1SSZEUIIIYRVk2RGCCGEEFZNkhkhhBBCWDWzJjPx8fHce++9uLq6Urt2bR5++GFOnz5t0iY2NhZFUUym++67z0wRCyGEEMLSmDWZ2b17N88++ywHDx5kx44dFBYW0rVrV3Jyckzade/enYsXLxqnL774wkwRCyGEEMLSmPXS7O3bt5vMJyYmUrt2bY4cOcIDDzxgLNdoNPj6+t7t8IQQQghhBSzqnJmsrCwAPD09TcqTk5OpXbs2DRo0YMSIEWRmZpbbR35+PtnZ2SaTEKLqZGZm8vTTTxMQEGD8odGtWzcOHDhQoeXXrFmDh4fHnQ1SCPGvYjE3zVNVlfHjx9OuXTvCw8ON5T169KBfv34EBgaSmprKtGnTePDBBzly5AgajaZUP/Hx8cyYMeNuhi7Ev0rfvn0pKChg7dq1BAcH89tvv/Hf//6XP/74467HUlBQgJ2d3V1frxDCwqgWYtSoUWpgYKB67ty5m7a7cOGCamdnp3700Udl1ufl5alZWVnG6dy5cyqgZmVl3Ymwhai03Nxc9eTJk2pubq6xzGAwqLrcXLNMBoOhwrFfuXJFBdTk5ORy28yfP18NDw9XnZyc1Lp166rPPPOM+tdff6mqqqq7du1SAZPp1VdfVVVVVQF169atJn25u7uriYmJqqqqampqqgqomzZtUtu3b69qNBo1ISFBHTp0qPrQQw+pc+fOVX19fVVPT0911KhRqk6nq/D7EmV/LoUwp6ysrArvvy1iZGbMmDF8+umn7Nmzh7p16960rZ+fH4GBgfz8889l1ms0mjJHbISwZIX5+SwZ+qhZ1j127WbsKviUZBcXF1xcXPj444+57777yvxbs7GxYcmSJdSrV4/U1FRGjRrFiy++yNtvv01UVBSLFi3ilVdeMV656OLiUql4J02axPz580lMTESj0bB792527dqFn58fu3bt4pdffmHAgAFERkYyYsSISvUthLBOZj1nRlVVRo8ezZYtW/j6668JCgq65TKXL1/m3Llz+Pn53YUIhRAl1ahRgzVr1rB27Vo8PDxo27YtU6ZM4YcffjC2GTduHB07diQoKIgHH3yQ119/nQ8++AAAe3t73N3dURQFX19ffH19K53MjBs3jj59+hAUFIS/vz8ANWvWZOnSpTRq1IhevXoRHR3Nf//736p740IIi2bWkZlnn32W9evX88knn+Dq6kpGRgYA7u7uODo6cvXqVaZPn07fvn3x8/MjLS2NKVOmUKtWLR555BFzhi5Elaqh0TB27Wazrbsy+vbtS3R0NN988w0HDhxg+/btzJkzh3fffZfY2Fh27drFrFmzOHnyJNnZ2RQWFpKXl0dOTg7Ozs7/ON6WLVuWKmvSpAm2trbGeT8/P1JSUv7xuoQQ1sGsyczy5csB6NChg0l5YmIisbGx2NrakpKSwnvvvceff/6Jn58fHTt2ZNOmTbi6upohYiHuDEVRKnyoxxI4ODjQpUsXunTpwiuvvMLw4cN59dVX6dixIz179mTkyJG8/vrreHp6snfvXoYNG0ZBQcFN+1QUpdRTtstapqyE6MaTgBVFwWAw3MY7E0JYI7MmMzd+cd3I0dGRL7/88i5FI4S4XY0bN+bjjz/mu+++o7CwkPnz52NjU3QUu/gQUzF7e3v0en2pPry9vbl48aJx/ueff+batWt3NnAhRLVgEScACyGsw+XLl+nXrx9xcXE0a9YMV1dXvvvuO+bMmcNDDz1E/fr1KSws5K233iImJoZ9+/bxzjvvmPRRr149rl69yn//+18iIiJwcnLCycmJBx98kKVLl3LfffdhMBiYNGmSXHYthKgQi7ppnhDCsrm4uNC6dWsWLlzIAw88QHh4ONOmTWPEiBEsXbqUyMhIFixYwOzZswkPDycpKYn4+HiTPqKiohg5ciQDBgzA29ubOXPmADB//ny0Wi0PPPAAAwcOZOLEiTg5OZnjbQohrIyi3upYj5XLzs7G3d2drKws3NzczB2OEOTl5ZGamkpQUBAOVnSejKje5HMpLE1l9t8yMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQogK69ChA+PGjStV/vHHH6Moyt0PSAghkGRGCHGHFRQUmDsEIUQ1J8mMEBZAVVUMOr1Zpqp+1uz06dOJjIwkISGB4OBgNBoNqqqyfft22rVrh4eHB15eXvTq1YszZ84Yl0tLS0NRFLZs2ULHjh1xcnIiIiKCAwcOlOq7pEWLFlGvXj3jfHJyMq1atcLZ2RkPDw/atm3Lr7/+WqXvUQhhWWqYOwAhBKgFBi68st8s6/Z/LQrF3rZK+/zll1/44IMP+Oijj7C1Leo7JyeH8ePH07RpU3JycnjllVd45JFHOH78ODY2f/+uevnll5k3bx6hoaG8/PLLPP744/zyyy/UqHHrr6vCwkIefvhhRowYwYYNG9DpdBw6dEgOgQlRzUkyI4Socjqdjvfffx9vb29jWd++fU3arF69mtq1a3Py5EnCw8ON5RMnTiQ6OhqAGTNm0KRJE3755RcaNWp0y/VmZ2eTlZVFr169qF+/PgBhYWFV8ZaEEBZMkhkhLIBiZ4P/a1FmW3dVCwwMNElkAM6cOcO0adM4ePAgly5dwmAwAJCenm6SzDRr1sz42s/PD4DMzMwKJTOenp7ExsbSrVs3unTpQufOnenfv7+xHyFE9STnzAhhARRFwcbe1ixTZQ7BuLm5kZWVVar8zz//xM3NzTjv7Oxcqk1MTAyXL19m1apVfPvtt3z77bdA0ShOSXZ2dibbBTAmPjY2NqXO8bnxBOPExEQOHDhAVFQUmzZtokGDBhw8eLDC71EIYX0kmRFCVFijRo347rvvSpUfPnyYhg0blrvc5cuXOXXqFFOnTqVTp06EhYVx5cqVSq/f29ubjIwMk4Tm+PHjpdo1b96cyZMns3//fsLDw1m/fn2l1yWEsB6SzAghKmzUqFGcOXOGZ599lu+//56ffvqJZcuWsXr1al544YVyl6tZsyZeXl6sXLmSX375ha+//prx48dXev0dOnTg999/Z86cOZw5c4Zly5axbds2Y31qaiqTJ0/mwIED/Prrr3z11Vf89NNPct6MENWcJDNCiAqrV68e33zzDWfOnKFr167ce++9rFmzhjVr1tCvX79yl7OxsWHjxo0cOXKE8PBwnn/+eebOnVvp9YeFhfH222+zbNkyIiIiOHToEBMnTjTWOzk58b///Y++ffvSoEEDnnrqKUaPHs3TTz99W+9XCGEdFLWqbzJhYbKzs3F3dycrK8vkmL4Q5pKXl0dqaipBQUE4ODiYOxwhAPlcCstTmf23jMwIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqZk1m4uPjuffee3F1daV27do8/PDDnD592qSNqqpMnz4df39/HB0d6dChAydOnDBTxEIIIYSwNGZNZnbv3s2zzz7LwYMH2bFjB4WFhXTt2pWcnBxjmzlz5rBgwQKWLl3K4cOH8fX1pUuXLvz1119mjFwIIYQQlsKsycz27duJjY2lSZMmREREkJiYSHp6OkeOHAGKRmUWLVrEyy+/TJ8+fQgPD2ft2rVcu3ZNnoIrhBnExsaiKApvvvmmSfnHH3+MoihVui5FUfj444+rtE8hRPVU458sfPLkSdLT09HpdCblvXv3vq3+srKyAPD09ASKnoCbkZFB165djW00Gg3t27dn//79ZT48Lj8/n/z8fON8dnb2bcUihCibg4MDs2fP5umnn6ZmzZrmDkcIIW5vZObs2bNEREQQHh5OdHQ0Dz/8MA8//DCPPPIIjzzyyG0Foqoq48ePp127doSHhwOQkZEBgI+Pj0lbHx8fY92N4uPjcXd3N05arfa24hHiblJVFZ1OZ5apss+a7dy5M76+vsTHx5fbZv/+/TzwwAM4Ojqi1WoZO3asyeHjevXq8frrrzNw4EBcXFzw9/fnrbfeMqkHeOSRR1AUxTgfGxvLww8/bLKucePG0aFDB+N8hw4dGDt2LC+++CKenp74+voyffp0k2WysrJ46qmnqF27Nm5ubjz44IN8//33ldoOQgjLcVsjM8899xxBQUHs3LmT4OBgDh06xOXLl5kwYQLz5s27rUBGjx7NDz/8wN69e0vV3Th8rapquUPakydPZvz48cb57OxsSWiExSsoKGDWrFlmWfeUKVOwt7evcHtbW1tmzZrFwIEDGTt2LHXr1jWpT0lJoVu3brz++uusXr2a33//ndGjRzN69GgSExON7ebOncuUKVOYPn06X375Jc8//zyNGjWiS5cuHD58mNq1a5OYmEj37t2xtbWt1Htau3Yt48eP59tvv+XAgQPExsbStm1bunTpgqqqREdH4+npyRdffIG7uzsrVqygU6dO/PTTT8aRYSGE9bitkZkDBw7w2muv4e3tjY2NDTY2NrRr1474+HjGjh1b6f7GjBnDp59+yq5du0y+GH19fQFKjcJkZmaWGq0pptFocHNzM5mEEFXrkUceITIykldffbVU3dy5cxk4cCDjxo0jNDSUqKgolixZwnvvvUdeXp6xXdu2bXnppZdo0KABY8aM4dFHH2XhwoUAeHt7A+Dh4YGvr69xvqKaNWvGq6++SmhoKEOGDKFly5b897//BWDXrl2kpKTw4Ycf0rJlS0JDQ5k3bx4eHh5s3rz5djeJEMKMbmtkRq/X4+LiAkCtWrW4cOECDRs2JDAwsNSl1Tejqipjxoxh69atJCcnExQUZFIfFBSEr68vO3bsoHnz5gDodDp2797N7Nmzbyd0ISySnZ0dU6ZMMdu6b8fs2bN58MEHmTBhgkn5kSNH+OWXX0hKSjKWqaqKwWAgNTWVsLAwANq0aWOyXJs2bVi0aNFtxXKjZs2amcz7+fmRmZlpjO/q1at4eXmZtMnNzeXMmTNVsn4hxN11W8lMeHg4P/zwA8HBwbRu3Zo5c+Zgb2/PypUrCQ4OrnA/zz77LOvXr+eTTz7B1dXVOALj7u6Oo6MjiqIwbtw4Zs2aRWhoKKGhocyaNQsnJycGDhx4O6ELYZEURanUoR5L8MADD9CtWzemTJlCbGyssdxgMPD000+XOUobEBBw0z5vdUWUjY2N8Ryf4v8XX4BQ8tyfGjVqmMwrioLBYEBVVfR6PX5+fuzatQtuOF3Iw8MD1VC5c4jKVMnzkMzqeqhqoQHVoKL/S0dhfsWvTCuzZUULK7qaCrcrv6FJVfHMjc2VG17cUF+02O0ti3Lrz7e4fbeVzEydOtV4Mt/MmTPp1asX999/P15eXmzatKnC/SxfvhzA5OQ9gMTEROOX44svvkhubi6jRo3iypUrtG7dmq+++gpXV9fbCb3K5J35k7zTf/z9ZVjWd9eNX2hltCl18uWNbcr6UizVppy+bre/Eu3UCrSpaF+oZTS5ne12Q32pMEoWlFd34ypKvtFytu+t4lNvUmdcToUCRxV9RA0KLl3D1k5PKZXYD5ZqWu6yt+i0gus0XCvAkFeI7sJVAF4f/wr3do2ivl89AHT/d5XIRs1IOfoDAY6+pTv4XYcOHehV9ifvRTfs7/tF7d/1DQ3qhaA7X1RmZ2dHfuZV4zyAp4M7Kb9+b1J2/PBR7GrYUfB/RTGp+XoMOQXGeQBDbiEGu6KyZgFhZGRkoF7Kp5420DQ+HRRcuMq/UUGhDn1WPr9//AM1/rKiZMya3SwhUm4sLyow5kMVaFuUQJVRd7P2t+q7ZIclqp1b++H6gOn5c3fTbSUz3bp1M74ODg7m5MmT/PHHH9SsWbNSmWdFrqJQFIXp06eXuhrB3HTn/uLqnv8zdxjCChW6Kqh6Z9RCFRWDucOpHPX6dH30IrxhYx5/pD9vJ664Xq8ycdQ47u/dibFTnmfYwFicnJz4388/8d9vvmbR639fIHDgu2+Z9/YieneL5r97dvHRZx/z8ZoPjfWBdQP4em8ybVreh8benpoeNekQ9QAL3lnMus3raX1PKzZs2cSJ06eIbGJ6WOlmOt3fkftatKLf8IG8MXkGDeqHcvG3DLZ//SW9u/XinogWVbGlrI9yfbJRoHLnW/+tItl1RfKkf0suVe4PpzJ/5ZVbU5lV3CmGa4V3aU1l+0f3mSnp33YFgH1dV1weqHN9Tinrfzdk3bcamiyeV242W/YwalmZehl9K2UFVlbuWcl2pm+znPdZumGF2ikVGdItr+56LKa/NMpY/sb6G+pK/xuUs1wZ6y4VoqKQr9dxreASNWpqqKFx4KZu+tugnMrbHcmuwHKKYw2UfFtq+DgZy15/8w02f7YVADsfZ1r4tib5611MfWUqDz7aHVVVqR9cn/79+2Pn61y0kK3C+PHj+f5kCm8sehNXV1fmzZlH9GMPGWOZv2A+E16YSMKGtdSpU4fUM2eJfvwhpp6eypQ3XyUvL48nY59k8JDB/JjyI3Z+RX0r9rbYONsZ5wFsHGyxcaqBnX/RuX5f7NjOy1Nf5ulJo/n999/x9fXlgfsfoE54kLFNVW+723E3D0sY8vKocc0B3+fDcHC4xefSAtxyVLvMupuNvBYfvryNZf/RaPENL9Qbi0oUlLU+k/ZlrKecdatqZdvfvK2tmwZzUtTK3mTCymRnZ+Pu7k5WVpZc2SQsQl5eHqmpqQQFBVnFTuNOqFevHuPGjWPcuHHmDkVcJ59LYWkqs/+Wp2YLIYQQwqpJMiOEEEIIq1Zl58wIIURFpaWlmTsEIUQ1IiMzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEqBLJyckoisKff/5p7lCEEP8ykswIISosNjYWRVFQFAU7OzuCg4OZOHEiOTk5t93fww8/XLVBCiH+deSmeUKISunevTuJiYkUFBTwzTffMHz4cHJychgwYIC5QxNC/EvJyIwQFkBVVfT6a2aZKvusWY1Gg6+vL1qtloEDBzJo0CA+/vjjUu2mT59OZGSkSdmiRYuoV6+esX7t2rV88sknxtGe5ORkAFJSUnjwwQdxdHTEy8uLp556iqtXr97GlhVC/BvIyIwQFsBgyCV5d1OzrLtD+xRsbZ1ue3lHR0cKCgoqvdzEiRM5deoU2dnZJCYmAuDp6cm1a9fo3r079913H4cPHyYzM5Phw4czevRo1qxZc9txCiGqL0lmhBC37dChQ6xfv55OnTpVelkXFxccHR3Jz8/H19fXWL527Vpyc3N57733cHZ2BmDp0qXExMQwe/ZsfHx8qix+IUT1IMmMEBbAxsaRDu1TzLbuyvjss89wcXGhsLCQgoICHnroId566y1OnjxZJfGcOnWKiIgIYyID0LZtWwwGA6dPn5ZkRghRiiQzQlgARVH+0aGeu6ljx44sX74cOzs7/P39sbOzAyiVzNjY2JQ6H6cih6NUVUVRlDLryisXQvy7yQnAQohKcXZ2JiQkhMDAQGMiUxZvb28yMjJMEprjx4+btLG3t0ev15uUNW7cmOPHj5tc7r1v3z5sbGxo0KBB1bwJIUS1IsmMEOKO6NChA7///jtz5szhzJkzLFu2jG3btpm0qVevHj/88AOnT5/m0qVLFBQUMGjQIBwcHBg6dCg//vgju3btYsyYMQwePFgOMQkhyiTJjBDijggLC+Ptt99m2bJlREREcOjQISZOnGjSZsSIETRs2JCWLVvi7e3Nvn37cHJy4ssvv+SPP/7g3nvv5dFHH6VTp04sXbrUTO9ECGHpFLWyN5mwMtnZ2bi7u5OVlYWbm5u5wxGCvLw8UlNTCQoKwsHBwdzhCAHI51JYnsrsv2VkRgghhBBWTZIZIYQQQlg1SWaEEEIIYdUkmRFCCCGEVZNkRgghhBBWTZIZIYQQQlg1SWaEEEIIYdUkmRFCCCGEVZNkRgghhBBWzazJzJ49e4iJicHf3x9FUfj4449N6mNjY1EUxWS67777zBOsEMJi1KtXj0WLFpk7jApJS0tDUZRSD9kUQlQdsyYzOTk5RERE3PSZK927d+fixYvG6YsvvriLEQohSir+gTFy5MhSdaNGjUJRFGJjY+94HIcPH+app56qsv6KE47iyd7enpCQEGbOnImlPfElPz+fMWPGUKtWLZydnenduzfnz583afPGG28QFRWFk5MTHh4e5glUiLuohjlX3qNHD3r06HHTNhqNBl9f37sUkRDiVrRaLRs3bmThwoU4OjoCRc/12bBhAwEBAXclBm9v7zvS786dO2nSpAn5+fns3buX4cOH4+fnx7Bhw+7I+m7HuHHj+M9//sPGjRvx8vJiwoQJ9OrViyNHjmBrawuATqejX79+tGnThtWrV5s5YiHuPIs/ZyY5OZnatWvToEEDRowYQWZm5k3b5+fnk52dbTIJYelUVSVHrzfLVNmRhxYtWhAQEMCWLVuMZVu2bEGr1dK8eXNj2fbt22nXrh0eHh54eXnRq1cvzpw5Y6zX6XSMHj0aPz8/HBwcqFevHvHx8cb66dOnExAQgEajwd/fn7FjxxrrSh5mevzxx3nsscdMYiwoKKBWrVokJiYat++cOXMIDg7G0dGRiIgINm/eXOq9eXl54evrS2BgIIMGDSIqKoqjR48a6w0GA6+99hp169ZFo9EQGRnJ9u3bTfo4dOgQzZs3x8HBgZYtW3Ls2DFjnaqqhISEMG/ePJNlfvzxR2xsbEy2T1mysrJYvXo18+fPp3PnzjRv3px169aRkpLCzp07je1mzJjB888/T9OmTW/anxDVhVlHZm6lR48e9OvXj8DAQFJTU5k2bRoPPvggR44cQaPRlLlMfHw8M2bMuMuRCvHPXDMYqL8nxSzrPvNAU5yv/6KvqCeffJLExEQGDRoEQEJCAnFxcSQnJxvb5OTkMH78eJo2bUpOTg6vvPIKjzzyCMePH8fGxoYlS5bw6aef8sEHHxAQEMC5c+c4d+4cAJs3b2bhwoVs3LiRJk2akJGRwffff19mLIMGDaJ///5cvXoVFxcXAL788ktycnLo27cvAFOnTmXLli0sX76c0NBQ9uzZwxNPPIG3tzft27cvs9/vvvuOo0ePMnToUGPZ4sWLmT9/PitWrKB58+YkJCTQu3dvTpw4QWhoKDk5OfTq1YsHH3yQdevWkZqaynPPPWdcXlEU4uLiSExMZOLEicbyhIQE7r//furXr3/T7X7kyBEKCgro2rWrsczf35/w8HD2799Pt27dbrq8ENWVRSczAwYMML4ODw+nZcuWBAYG8vnnn9OnT58yl5k8eTLjx483zmdnZ6PVau94rEL8mwwePJjJkycbzzXZt28fGzduNElmihOJYqtXr6Z27dqcPHmS8PBw0tPTCQ0NpV27diiKQmBgoLFteno6vr6+dO7cGTs7OwICAmjVqlWZsXTr1g1nZ2e2bt3K4MGDAVi/fj0xMTG4ubmRk5PDggUL+Prrr2nTpg0AwcHB7N27lxUrVpgkM1FRUdjY2KDT6SgoKOCpp55iyJAhxvp58+YxadIk40jQ7Nmz2bVrF4sWLWLZsmUkJSWh1+tJSEjAycmJJk2acP78eZ555hljH08++SSvvPIKhw4dolWrVhQUFLBu3Trmzp17y+2ekZGBvb09NWvWNCn38fEhIyPjlssLUV1ZdDJzIz8/PwIDA/n555/LbaPRaModtRHCUjnZ2HDmAfMcEnCyqfzR5lq1ahEdHc3atWtRVZXo6Ghq1apl0ubMmTNMmzaNgwcPcunSJQwGA1CUqISHhxMbG0uXLl1o2LAh3bt3p1evXsYRh379+rFo0SKCg4Pp3r07PXv2JCYmhho1Sn9l2dnZ0a9fP5KSkhg8eDA5OTl88sknrF+/HoCTJ0+Sl5dHly5dTJbT6XQmh8UANm3aRFhYGAUFBaSkpDB27Fhq1qzJm2++SXZ2NhcuXKBt27Ymy7Rt29Y4anTq1CkiIiJwcnIy1hcnUMX8/PyIjo4mISGBVq1a8dlnn5GXl0e/fv0qvP1vpKoqiqLc9vJCWDurSmYuX77MuXPn8PPzM3coQlQpRVEqfajH3OLi4hg9ejQAy5YtK1UfExODVqtl1apV+Pv7YzAYCA8PR6fTAUXn3qSmprJt2zZ27txJ//796dy5M5s3b0ar1XL69Gl27NjBzp07GTVqFHPnzmX37t3Y2dmVWtegQYNo3749mZmZ7NixAwcHB+PFBcVJ1Oeff06dOnVMlrvxh49WqyUkJASAsLAwzp49y7Rp05g+fbqxzY1JQ8lEoqLnHw0fPpzBgwezcOFCEhMTGTBggEkCVB5fX190Oh1XrlwxGZ3JzMwkKiqqQusWojoy6wnAV69e5fjx48b7L6SmpnL8+HHS09O5evUqEydO5MCBA6SlpZGcnExMTAy1atXikUceMWfYQgiKbpug0+nQ6XSlztW4fPkyp06dYurUqXTq1ImwsDCuXLlSqg83NzcGDBjAqlWr2LRpEx999BF//PEHAI6OjvTu3ZslS5aQnJzMgQMHSEkp+7yiqKgotFotmzZtIikpiX79+mFvbw9A48aN0Wg0pKenExISYjLd6hC0ra0thYWF6HQ63Nzc8Pf3Z+/evSZt9u/fT1hYmHFd33//Pbm5ucb6gwcPluq3Z8+eODs7s3z5crZt20ZcXNxN4yh2zz33YGdnx44dO4xlFy9e5Mcff5RkRvyrmXVk5rvvvqNjx47G+eJzXYYOHcry5ctJSUnhvffe488//8TPz4+OHTuyadMmXF1dzRWyEOI6W1tbTp06ZXxdUs2aNfHy8mLlypX4+fmRnp7OSy+9ZNJm4cKF+Pn5ERkZiY2NDR9++CG+vr54eHiwZs0a9Ho9rVu3xsnJiffffx9HR0eT82pKUhSFgQMH8s477/DTTz+xa9cuY52rqysTJ07k+eefx2Aw0K5dO7Kzs9m/fz8uLi4mJ/hevnyZjIwMCgsLSUlJYfHixXTs2BE3NzcAXnjhBV599VXq169PZGQkiYmJHD9+nKSkJAAGDhzIyy+/zLBhw5g6dSppaWmlrlwq3l6xsbFMnjyZkJCQUoeiyuPu7s6wYcOYMGECXl5eeHp6MnHiRJo2bUrnzp2N7dLT0/njjz9IT09Hr9cbfzCGhIQYT5IWolpRq7msrCwVULOysswdihCqqqpqbm6uevLkSTU3N9fcoVTa0KFD1Yceeqjc+oceekgdOnSoqqqqumPHDjUsLEzVaDRqs2bN1OTkZBVQt27dqqqqqq5cuVKNjIxUnZ2dVTc3N7VTp07q0aNHVVVV1a1bt6qtW7dW3dzcVGdnZ/W+++5Td+7caVxPYGCgunDhQpN1nzhxQgXUwMBA1WAwmNQZDAZ18eLFasOGDVU7OzvV29tb7datm7p7925VVVU1NTVVBYyTra2tWrduXXXEiBFqZmamsR+9Xq/OmDFDrVOnjmpnZ6dGRESo27ZtM1nXgQMH1IiICNXe3l6NjIxUP/roIxVQjx07ZtLuzJkzKqDOmTPnVpvdRG5urjp69GjV09NTdXR0VHv16qWmp6ebtBk6dKjJ+ymedu3addN+rfVzKaqnyuy/FVW1sNtbVrHs7Gzc3d3Jysoy/roSwpzy8vJITU0lKCgIBwcHc4cjzGTfvn106NCB8+fP4+PjY+5w5HMpLE5l9t9WdQKwEEJYu/z8fM6dO8e0adPo37+/RSQyQlg7i78DsBBCVCcbNmygYcOGZGVlMWfOHJO6pKQkXFxcypyaNGlipoiFsHwyMiOEEHdRbGxsuQ/j7N27N61bty6zrqxL0oUQRSSZEUIIC+Hq6ipXawpxG+QwkxBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgirU69ePRYtWmTuMCokLS0NRVGMz0cSQlQ9SWaEEBUWGxuLoiiMHDmyVN2oUaNQFKXce6hUpcOHD/PUU09VWX/FCUfxZG9vT0hICDNnzsTSnviSn5/PmDFjqFWrFs7OzvTu3Zvz588b69PS0hg2bBhBQUE4OjpSv359Xn31VXQ6nRmjFuLOkmRGCFEpWq2WjRs3kpubayzLy8tjw4YNBAQE3JUYvL29cXJyqvJ+d+7cycWLF/n555+ZMWMGb7zxBgkJCVW+nn9i3LhxbN26lY0bN7J3716uXr1Kr1690Ov1APzvf//DYDCwYsUKTpw4wcKFC3nnnXeYMmWKmSMX4s6RZEYIC6CqKtd0hWaZKjvy0KJFCwICAtiyZYuxbMuWLWi1Wpo3b24s2759O+3atcPDwwMvLy969erFmTNnjPU6nY7Ro0fj5+eHg4MD9erVIz4+3lg/ffp0AgIC0Gg0+Pv7M3bsWGNdycNMjz/+OI899phJjAUFBdSqVYvExETj9p0zZw7BwcE4OjoSERHB5s2bS703Ly8vfH19CQwMZNCgQURFRXH06FFjvcFg4LXXXqNu3bpoNBoiIyPZvn27SR+HDh2iefPmODg40LJlS44dO2asU1WVkJAQ5s2bZ7LMjz/+iI2Njcn2KUtWVharV69m/vz5dO7cmebNm7Nu3TpSUlLYuXMnAN27dycxMZGuXbsSHBxM7969mThxosm/lxDVjdwBWAgLkFugp/ErX5pl3Sdf64aTfeW+Cp588kkSExMZNGgQAAkJCcTFxZGcnGxsk5OTw/jx42natCk5OTm88sorPPLIIxw/fhwbGxuWLFnCp59+ygcffEBAQADnzp3j3LlzAGzevJmFCxeyceNGmjRpQkZGBt9//32ZsQwaNIj+/ftz9epVXFxcAPjyyy/Jycmhb9++AEydOpUtW7awfPlyQkND2bNnD0888QTe3t60b9++zH6/++47jh49ytChQ41lixcvZv78+axYsYLmzZuTkJBA7969OXHiBKGhoeTk5NCrVy8efPBB1q1bR2pqKs8995xxeUVRiIuLIzExkYkTJxrLExISuP/++6lfv/5Nt/uRI0coKCiga9euxjJ/f3/Cw8PZv38/3bp1K3O5rKwsPD09b9q3ENZMkhkhRKUNHjyYyZMnG8812bdvHxs3bjRJZooTiWKrV6+mdu3anDx5kvDwcNLT0wkNDaVdu3YoikJgYKCxbXp6Or6+vnTu3Bk7OzsCAgJo1apVmbF069YNZ2dntm7dyuDBgwFYv349MTExuLm5kZOTw4IFC/j6669p06YNAMHBwezdu5cVK1aYJDNRUVHY2Nig0+koKCjgqaeeYsiQIcb6efPmMWnSJONI0OzZs9m1axeLFi1i2bJlJCUlodfrSUhIwMnJiSZNmnD+/HmeeeYZYx9PPvkkr7zyCocOHaJVq1YUFBSwbt065s6de8vtnpGRgb29PTVr1jQp9/HxISMjo8xlzpw5w1tvvcX8+fNv2b8Q1kqSGSEsgKOdLSdfK/tX9d1Yd2XVqlWL6Oho1q5di6qqREdHU6tWLZM2Z86cYdq0aRw8eJBLly5hMBiAokQlPDyc2NhYunTpQsOGDenevTu9evUyjjj069ePRYsWERwcTPfu3enZsycxMTHUqFH6K8vOzo5+/fqRlJTE4MGDycnJ4ZNPPmH9+vUAnDx5kry8PLp06WKynE6nMzksBrBp0ybCwsIoKCggJSWFsWPHUrNmTd58802ys7O5cOECbdu2NVmmbdu2xlGjU6dOERERYXI+T3ECVczPz4/o6GgSEhJo1aoVn332GXl5efTr16/C2/9GqqqiKEqp8gsXLtC9e3f69evH8OHDb7t/ISydJDNCWABFUSp9qMfc4uLiGD16NADLli0rVR8TE4NWq2XVqlX4+/tjMBgIDw83XlXTokULUlNT2bZtGzt37qR///507tyZzZs3o9VqOX36NDt27GDnzp2MGjWKuXPnsnv37jKfHj1o0CDat29PZmYmO3bswMHBgR49egAYk6jPP/+cOnXqmCyn0WhM5rVaLSEhIQCEhYVx9uxZpk2bxvTp041tbkwaSiYSFT3/aPjw4QwePJiFCxeSmJjIgAEDKnRCs6+vLzqdjitXrpiMzmRmZhIVFWXS9sKFC3Ts2JE2bdqwcuXKCsUlhLWSE4CFELele/fu6HQ6dDpdqXM1Ll++zKlTp5g6dSqdOnUiLCyMK1eulOrDzc2NAQMGsGrVKjZt2sRHH33EH3/8AYCjoyO9e/dmyZIlJCcnc+DAAVJSUsqMJSoqCq1Wy6ZNm0hKSqJfv37Y29sD0LhxYzQaDenp6YSEhJhMWq32pu/R1taWwsJCdDodbm5u+Pv7s3fvXpM2+/fvJywszLiu77//3uRKr4MHD5bqt2fPnjg7O7N8+XK2bdtGXFzcTeMods8992BnZ8eOHTuMZRcvXuTHH380SWb+7//+jw4dOtCiRQsSExOxsZGvelG9WddPQSGExbC1teXUqVPG1yXVrFkTLy8vVq5ciZ+fH+np6bz00ksmbRYuXIifnx+RkZHY2Njw4Ycf4uvri4eHB2vWrEGv19O6dWucnJx4//33cXR0NDmvpiRFURg4cCDvvPMOP/30E7t27TLWubq6MnHiRJ5//nkMBgPt2rUjOzub/fv34+LiYnKC7+XLl8nIyKCwsJCUlBQWL15Mx44dcXNzA+CFF17g1VdfpX79+kRGRpKYmMjx48dJSkoCYODAgbz88ssMGzaMqVOnkpaWVurKpeLtFRsby+TJkwkJCSl1KKo87u7uDBs2jAkTJuDl5YWnpycTJ06kadOmdO7cGSgakenQoQMBAQHMmzeP33//3bi8r69vhdYjhNVRq7msrCwVULOysswdihCqqqpqbm6uevLkSTU3N9fcoVTa0KFD1Yceeqjc+oceekgdOnSoqqqqumPHDjUsLEzVaDRqs2bN1OTkZBVQt27dqqqqqq5cuVKNjIxUnZ2dVTc3N7VTp07q0aNHVVVV1a1bt6qtW7dW3dzcVGdnZ/W+++5Td+7caVxPYGCgunDhQpN1nzhxQgXUwMBA1WAwmNQZDAZ18eLFasOGDVU7OzvV29tb7datm7p7925VVVU1NTVVBYyTra2tWrduXXXEiBFqZmamsR+9Xq/OmDFDrVOnjmpnZ6dGRESo27ZtM1nXgQMH1IiICNXe3l6NjIxUP/roIxVQjx07ZtLuzJkzKqDOmTPnVpvdRG5urjp69GjV09NTdXR0VHv16qWmp6cb6xMTE03eS8npVv1a6+dSVE+V2X8rqmpht7esYtnZ2bi7u5OVlWX8dSWEOeXl5ZGamkpQUBAODg7mDkeYyb59++jQoQPnz5/Hx8fH3OHI51JYnMrsv+UwkxBC3EX5+fmcO3eOadOm0b9/f4tIZISwdnJWmBBC3EUbNmygYcOGZGVlMWfOHJO6pKQkXFxcypyaNGlipoiFsHwyMiOEEHdRbGxsuQ/j7N27N61bty6zrqxL0oUQRSSZEUIIC+Hq6oqrq6u5wxDC6shhJiGEEEJYNUlmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBBWp169eixatMjcYVRIWloaiqJw/Phxc4ciRLUlyYwQosJiY2NRFIWRI0eWqhs1ahSKopR7D5WqdPjwYZ566qkq66844Sie7O3tCQkJYebMmVjaE1/y8/MZM2YMtWrVwtnZmd69e3P+/HmTNr179yYgIAAHBwf8/PwYPHgwFy5cMFPEQtx5Zk1m9uzZQ0xMDP7+/iiKwscff2xSr6oq06dPx9/fH0dHRzp06MCJEyfME6wQAgCtVsvGjRvJzc01luXl5bFhwwYCAgLuSgze3t44OTlVeb87d+7k4sWL/Pzzz8yYMYM33niDhISEKl/PPzFu3Di2bt3Kxo0b2bt3L1evXqVXr17o9Xpjm44dO/LBBx9w+vRpPvroI86cOcOjjz5qxqiFuLPMmszk5OQQERHB0qVLy6yfM2cOCxYsYOnSpRw+fBhfX1+6dOnCX3/9dZcjFeIOU1XQ5ZhnquTIQ4sWLQgICGDLli3Gsi1btqDVamnevLmxbPv27bRr1w4PDw+8vLzo1asXZ86cMdbrdDpGjx6Nn58fDg4O1KtXj/j4eGP99OnTCQgIQKPR4O/vz9ixY411JQ8zPf744zz22GMmMRYUFFCrVi0SExOvb16VOXPmEBwcjKOjIxEREWzevLnUe/Py8sLX15fAwEAGDRpEVFQUR48eNdYbDAZee+016tati0ajITIyku3bt5v0cejQIZo3b46DgwMtW7bk2LFjxjpVVQkJCWHevHkmy/z444/Y2NiYbJ+yZGVlsXr1aubPn0/nzp1p3rw569atIyUlhZ07dxrbPf/889x3330EBgYSFRXFSy+9xMGDBykoKLhp/0JYK7PeAbhHjx706NGjzDpVVVm0aBEvv/wyffr0AWDt2rX4+Piwfv16nn766TKXy8/PJz8/3zifnZ1d9YELUdUKrsEsf/Ose8oFsHeu1CJPPvkkiYmJDBo0CICEhATi4uJITk42tsnJyWH8+PE0bdqUnJwcXnnlFR555BGOHz+OjY0NS5Ys4dNPP+WDDz4gICCAc+fOce7cOQA2b97MwoUL2bhxI02aNCEjI4Pvv/++zFgGDRpE//79uXr1Ki4uLgB8+eWX5OTk0LdvXwCmTp3Kli1bWL58OaGhoezZs4cnnngCb29v2rdvX2a/3333HUePHmXo0KHGssWLFzN//nxWrFhB8+bNSUhIoHfv3pw4cYLQ0FBycnLo1asXDz74IOvWrSM1NZXnnnvOuLyiKMTFxZGYmMjEiRON5QkJCdx///3Ur1//ptv9yJEjFBQU0LVrV2OZv78/4eHh7N+/n27dupVa5o8//iApKYmoqCh5JIKotiz2nJnU1FQyMjJM/mg1Gg3t27dn//795S4XHx+Pu7u7cdJqtXcjXCH+VQYPHszevXtJS0vj119/Zd++fTzxxBMmbfr27UufPn0IDQ0lMjKS1atXk5KSwsmTJwFIT08nNDSUdu3aERgYSLt27Xj88ceNdb6+vnTu3JmAgABatWrFiBEjyoylW7duODs7s3XrVmPZ+vXriYmJwc3NjZycHBYsWEBCQgLdunUjODiY2NhYnnjiCVasWGHSV1RUFC4uLtjb23PvvffSv39/hgwZYqyfN28ekyZN4rHHHqNhw4bMnj2byMhI4yhRUlISer2ehIQEmjRpQq9evXjhhRdM1vHkk09y+vRpDh06BBSNIq1bt464uLhbbveMjAzs7e2pWbOmSbmPjw8ZGRkmZZMmTcLZ2RkvLy/S09P55JNPbtm/ENbKYp/NVPyH6ePjY1Lu4+PDr7/+Wu5ykydPZvz48cb57OxsSWiE5bNzKhohMde6K6lWrVpER0ezdu1aVFUlOjqaWrVqmbQ5c+YM06ZN4+DBg1y6dAmDwQAUJSrh4eHExsbSpUsXGjZsSPfu3enVq5fxx0u/fv1YtGgRwcHBdO/enZ49exITE0ONGqW/suzs7OjXrx9JSUkMHjyYnJwcPvnkE9avXw/AyZMnycvLo0uXLibL6XQ6k8NiAJs2bSIsLIyCggJSUlIYO3YsNWvW5M033yQ7O5sLFy7Qtm1bk2Xatm1rHDU6deoUERERJufztGnTxqS9n58f0dHRJCQk0KpVKz777DPy8vLo169fhbf/jVRVRVEUk7IXXniBYcOG8euvvzJjxgyGDBnCZ599VqqdENWBxSYzxW78wyvrj7YkjUaDRqO502EJUbUUpdKHeswtLi6O0aNHA7Bs2bJS9TExMWi1WlatWoW/vz8Gg4Hw8HB0Oh1QdO5Namoq27ZtY+fOnfTv35/OnTuzefNmtFotp0+fZseOHezcuZNRo0Yxd+5cdu/eXeahkkGDBtG+fXsyMzPZsWMHDg4OxkPYxUnU559/Tp06dUyWu/G7QqvVEhISAkBYWBhnz55l2rRpTJ8+3djmZt9JFb3yafjw4QwePJiFCxeSmJjIgAEDKnRCs6+vLzqdjitXrpiMzmRmZhIVFWXStlatWtSqVYsGDRoQFhaGVqvl4MGDpZIrIaoDi01mfH19gaIRGj8/P2N5ZmZmqdEaczj25Wcc+uT6CYQlv8BKvDb5WrteXt6XnVqJPkq3KVleRn8l6m/5XXuTBiq3WPhm1bdY8S37/gcU/sEv0X/yK7aMRRUUHGt60uzRJ7hkZ4udrW3VdF2hMP/5L/K8q3+Rfy2HzLSztGjUgLzrVzQ1bxhKZtpZ8q/lkFfDlv8dO8KpU6d4c/qrNK0fBMC3330HwJ+ZGWT+etbYZ8f7WtHxvlZ0vr8djw2N5fT3x/D08ACgTURT2kQ05fGHHyKqU2f2fLWdZuHhGAoLufrHZX7/NRWA0Dp+1PHzY/Xyt/lvcjK9uncj6+L/AVDbxQmNvT0/Hj1C43o3XG1lKOT39DQu/1/Rpc1/XLzA754exurc7CwKCwu5cOZnXF1d8fXx4cvP/kNYiX72JO+ieUQEv6enofWtzXtrj5H+0/9wdHAAYMcXn5fq+97wxjg6ODBv1hts27aNTz/YxKX0NMr9N7peXM/HGzs7Oz7asJ6HY3oBkPFbJj/++CMvvzCBS+fKHrW+fP2y7Mxzv3KpbtnnZhUUFnL1ymW2fvAeur/KOdewnL+Hm36yyvsbKrevm/RWblc3i6B0XZnNyygsM5Z/tGwFY7lh2VLvr9TsrddVej03dnKLdZahZJuGbe4nvGOXm7S+syw2mQkKCsLX15cdO3YYh4J1Oh27d+9m9uzZZo4OdLm5XL18ydxhCCtkp9OBqqIaDKi3mSyZ684nqqqiqioGvR4F2LN9G1D0tWjQ6431bi4u1KzpwXvr1+Pt5cX/XbzAG3OLruBRDQYMhXpWJCTiU9ubJmFh2NjY8Olnn1Hb2xtXJyfWb/oAvV5Pi8gIHB0c2PTRR0X3TPHxQV9YiErRiIu+sNAY28MxvVizLomzaWlsXveesc7RwYGRw4cx7fXX0RcW0Oqelvx19SrfHT2Ks7MT/fv0wVBYdJXP5Uu/c/HCBQoL9Zz66TSrEhJpe999ODk4oC8o4Jnhw5i3eAkBdevQJCyMTZs/4seTp1g6fx76ggIe7tmTWXPmMW7iCzw3ahTn/+88b69YCYChsAB9iauJ+vd5hJmz51AvMIDmzZpSWIErjZwcHHi836O8+vpM3F1d8HD34LU33ySsYQPatmpFoU7Hse+/59j3P9Cq5T24u7uTnn6OOYsXUy8ggMjwphReHxm7UaFej6FQz58ZF7n2h3y3icqpHXTzk9fvNLMmM1evXuWXX34xzqempnL8+HE8PT0JCAhg3LhxzJo1i9DQUEJDQ5k1axZOTk4MHDjQjFEXCe/QmXrNShxvL7FTKi+jNZaXbGvaoJw+lLJfltPGWGyyHpMFy4yvYtW3Wrbyv6oq0eA2/IPd/j9ZtMyRqKIyXUEBmVeyqOlXB4frv94rtdI7lsmU7vjGEo2zM/YFhXjVLRqZ8Lqh3t7RCY2TM94B9Vi/LonnJ0ygY3QvGjRowML58+nctSuuXt541Q3Ap05d3lm5kl9++QVbW1ta3nMP//nPf/AOqEedoGDmzpvLjDdno9frCQ9vwsdbthDatBkAtra2OHt44FX37/Ph4p56iiXL3yEwIIAevR82+SzOmTefwOD6LFu5kokvT8PDw4PmkZFMevFFPOtoyS4sOhTVf0issX8/X196REfz2vTpeHp7A/DilJfR29jw+uy5ZP7+O2Fhjfho84e0bNsOAE/g461bGT12DF0fepiwsEa8+WY8Ax4fiHttHzzr1DVu2GdGj2HJ8neIixuGp3/dCv8rvbV0GS9NmcLIcePJzc2lY4cOrFn7Ht516wIqPn9cYceSt5i/dBk5OTn4+vrStUsXJk+ahF+dv0dlbvy3zc/P58/cfLo9M44aNqX/Fss9hHaTz2O5I67l9nWzkeHyKio5mlxW0U3+Xm+5qjIKy15vBWO5sfCG5Up3U1act+jjVrGV+Z5uvoyXNrDUMneToprx9pbJycl07NixVPnQoUNZs2YNqqoyY8YMVqxYwZUrV2jdujXLli0jPDy8wuvIzs7G3d2drKws3NzcqjJ8IW5LXl4eqampBAUFlZHMiH+Lffv20aFDB86fP28Rh87lcyksTWX232ZNZu4GSWaEpZGdxr9bfn4+586d46mnnsLPz4+kpCRzhwTI51JYnsrsvy32PjNCCFEdbdiwgYYNG5KVlcWcOXNM6pKSknBxcSlzatKkiZkiFsLyWewJwEIIUR3FxsaW+zDO3r1707p16zLr5O69QpRPkhkhhLAQrq6uuLq6mjsMIayOHGYSQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBBCCGHVJJkRQgghhFWTZEYIYXXq1avHokWLzB1GhaSlpaEoCsePHzd3KEJUW5LMCCEqLDY2FkVRGDlyZKm6UaNGoShKufdQqUqHDx/mqaeeqrL+ihOO4sne3p6QkBBmzpxZ/jOJzCQ/P58xY8ZQq1YtnJ2d6d27N+fPny+3bWRkpCRTotqTZEYIUSlarZaNGzeSm5trLMvLy2PDhg0EBATclRi8vb1xcnKq8n537tzJxYsX+fnnn5kxYwZvvPEGCQkJVb6ef2LcuHFs3bqVjRs3snfvXq5evUqvXr3Q6/Wl2r744ov4+/uX0YsQ1YskM0JYAFVVuVZwzSxTZUceWrRoQUBAAFu2bDGWbdmyBa1WS/Pmfz9Jfvv27bRr1w4PDw+8vLzo1asXZ86cMdbrdDpGjx6Nn58fDg4O1KtXj/j4eGP99OnTCQgIQKPR4O/vz9ixY411JQ8zPf744zz22GMmMRYUFFCrVi0SExON23fOnDkEBwfj6OhIREQEmzdvLvXevLy88PX1JTAwkEGDBhEVFcXRo0eN9QaDgddee426deui0WiIjIxk+/btJn0cOnSI5s2b4+DgQMuWLTl27JixTlVVQkJCmDdvnskyP/74IzY2NibbpyxZWVmsXr2a+fPn07lzZ5o3b866detISUlh586dJm23bdvGV199VWpdQlRHcgdgISxAbmEurdeXfRv7O+3bgd/iZFe5UY4nn3ySxMREBg0aBEBCQgJxcXEkJycb2+Tk5DB+/HiaNm1KTk4Or7zyCo888gjHjx/HxsaGJUuW8Omnn/LBBx8QEBDAuXPnOHfuHACbN29m4cKFbNy4kSZNmpCRkcH3339fZiyDBg2if//+XL16FRcXFwC+/PJLcnJy6Nu3LwBTp05ly5YtLF++nNDQUPbs2cMTTzyBt7c37du3L7Pf7777jqNHjzJ06FBj2eLFi5k/fz4rVqygefPmJCQk0Lt3b06cOEFoaCg5OTn06tWLBx98kHXr1pGamspzzz1nXF5RFOLi4khMTGTixInG8oSEBO6//37q169/0+1+5MgRCgoK6Nq1q7HM39+f8PBw9u/fT7du3QD47bffGDFiBB9//PEdGcESwtJIMiOEqLTBgwczefJk47km+/btY+PGjSbJTHEiUWz16tXUrl2bkydPEh4eTnp6OqGhobRr1w5FUQgMDDS2TU9Px9fXl86dO2NnZ0dAQACtWrUqM5Zu3brh7OzM1q1bGTx4MADr168nJiYGNzc3cnJyWLBgAV9//TVt2rQBIDg4mL1797JixQqTZCYqKgobGxt0Oh0FBQU89dRTDBkyxFg/b948Jk2aZBwJmj17Nrt27WLRokUsW7aMpKQk9Ho9CQkJODk50aRJE86fP88zzzxj7OPJJ5/klVde4dChQ7Rq1YqCggLWrVvH3Llzb7ndMzIysLe3p2bNmiblPj4+ZGRkAEWjP7GxsYwcOZKWLVuSlpZ2y36FsHaSzAhhARxrOPLtwG/Ntu7KqlWrFtHR0axduxZVVYmOjqZWrVombc6cOcO0adM4ePAgly5dwmAwAEWJSnh4OLGxsXTp0oWGDRvSvXt3evXqZRxx6NevH4sWLSI4OJju3bvTs2dPYmJiqFGj9FeWnZ0d/fr1IykpicGDB5OTk8Mnn3zC+vXrATh58iR5eXl06dLFZDmdTmdyWAxg06ZNhIWFUVBQQEpKCmPHjqVmzZq8+eabZGdnc+HCBdq2bWuyTNu2bY2jRqdOnSIiIsJkNKQ4gSrm5+dHdHQ0CQkJtGrVis8++4y8vDz69etX4e1/I1VVURQFgLfeeovs7GwmT5582/0JYW0kmRHCAiiKUulDPeYWFxfH6NGjAVi2bFmp+piYGLRaLatWrcLf3x+DwUB4eDg6nQ4oOvcmNTWVbdu2sXPnTvr370/nzp3ZvHkzWq2W06dPs2PHDnbu3MmoUaOYO3cuu3fvLvPp0YMGDaJ9+/ZkZmayY8cOHBwc6NGjB4Axifr888+pU6eOyXIajcZkXqvVEhISAkBYWBhnz55l2rRpTJ8+3dimOGkoVjKRqOj5R8OHD2fw4MEsXLiQxMREBgwYUKHDQb6+vuh0Oq5cuWIyOpOZmUlUVBQAX3/9NQcPHiz13lq2bMmgQYNYu3ZthWIUwprICcBCiNvSvXt3dDodOp3OeK5GscuXL3Pq1CmmTp1Kp06dCAsL48qVK6X6cHNzY8CAAaxatYpNmzbx0Ucf8ccffwDg6OhI7969WbJkCcnJyRw4cICUlJQyY4mKikKr1bJp0yaSkpLo168f9vb2ADRu3BiNRkN6ejohISEmk1arvel7tLW1pbCwEJ1Oh5ubG/7+/uzdu9ekzf79+wkLCzOu6/vvvze50uvgwYOl+u3ZsyfOzs4sX76cbdu2ERcXd9M4it1zzz3Y2dmxY8cOY9nFixf58ccfjcnMkiVL+P777zl+/DjHjx/niy++AIpGnd54440KrUcIayMjM0KI22Jra8upU6eMr0uqWbMmXl5erFy5Ej8/P9LT03nppZdM2ixcuBA/Pz8iIyOxsbHhww8/xNfXFw8PD9asWYNer6d169Y4OTnx/vvv4+joaHJeTUmKojBw4EDeeecdfvrpJ3bt2mWsc3V1ZeLEiTz//PMYDAbatWtHdnY2+/fvx8XFxeQE38uXL5ORkUFhYSEpKSksXryYjh074ubmBsALL7zAq6++Sv369YmMjCQxMZHjx4+TlJQEwMCBA3n55ZcZNmwYU6dOJS0trcyriWxtbYmNjWXy5MmEhISUOhRVHnd3d4YNG8aECRPw8vLC09OTiRMn0rRpUzp37gxQ6vL44pOi69evT926dSu0HiGsjSQzQojbVryTv5GNjQ0bN25k7NixhIeH07BhQ5YsWUKHDh2MbVxcXJg9ezY///wztra23HvvvXzxxRfY2Njg4eHBm2++yfjx49Hr9TRt2pT//Oc/eHl5lRvLoEGDmDVrFoGBgaXOa3n99depXbs28fHxnD17Fg8PD1q0aMGUKVNM2hUnBLa2tvj5+dGzZ0+T0YyxY8eSnZ3NhAkTyMzMpHHjxnz66aeEhoYa39N//vMfRo4cSfPmzWncuDGzZ88udTI0wLBhw5g1a1aFR2WKLVy4kBo1atC/f39yc3Pp1KkTa9asKZVQCvFvoqiWdnvLKpadnY27uztZWVnlfvEKcTfl5eWRmppKUFAQDg4O5g5HmMm+ffvo0KED58+fx8fHx9zhyOdSWJzK7L9lZEYIIe6i/Px8zp07x7Rp0+jfv79FJDJCWDs5AVgIIe6iDRs20LBhQ7KyspgzZ45JXVJSEi4uLmVOTZo0MVPEQlg+GZkRQoi7KDY2ttyHcfbu3ZvWrcu+E3RZl6QLIYpIMiOEEBbC1dUVV1dXc4chhNWRw0xCCCGEsGqSzAghhBDCqkkyI4QQQgirJsmMEEIIIayaJDNCCCGEsGqSzAghrE69evVYtGiRucOokLS0NBRF4fjx4+YORYhqS5IZIUSFxcbGoigKI0eOLFU3atQoFEUp9x4qVenw4cM89dRTVdZfccJRPNnb2xMSEsLMmTOxtCe+5OfnM2bMGGrVqoWzszO9e/fm/PnzJm3q1atn8n4URSn1oE8hqhNJZoQQlaLVatm4cSO5ubnGsry8PDZs2FDqic13ire3N05OTlXe786dO7l48SI///wzM2bM4I033iAhIaHK1/NPjBs3jq1bt7Jx40b27t3L1atX6dWrF3q93qTda6+9xsWLF43T1KlTzRSxEHeeRd80b/r06cyYMcOkzMfHh4yMDDNF9LdfT1zm7PHfzR2GsEI29gbcggq5eiUPnb1l/eq/lYJ8Pc3CI0j7NY2ktRvp3+8xAD7cvIk6/nWpF1iPgnw92Zdz2fnfr5i7YDanTp3ExtaWVi1b8easeQQHBQOg0+mYMm0Sn/7nY/7M+hOf2j7EDh3GhHEvABA/eybr1r9H5u+ZeNb05KHejzAnfj4ATZs34pmnRzNq5GjiRgxFVVUS333v7zgLCmjQJJjXp7/BEwOHoKoqi99aSMLad/nttwxC6ofywoSXeLj3IwD8dSUPAI2tM0527jjZQUz3PrS6910O7j9Ev4cHAmAwGJg7/03WvJfApcuXaNigIdOnvU7nTl2N6z5y9DDPTRjDTz+dJqxRYyaOfxGAq3/mkXXpGs1bNSVu6HCeGz3OuMzJUydo80Arjh1KMW6fsmRlZ7F69WpWvr2a1s2Lngz+zpJVhEU04NMtn9P5wS4AqAYVOxsHnO3cjcuq+fBXfm6Z/Rb9e+STl1PAof+cxaCT37nmYl3fCH8LbOxFcHNvs63fopMZgCZNmrBz507jvKU85v7Sub84+c0Fc4chrJCDuw1N/dzJv1Zo3Gmoqgr5eeYJSOOAoigVaqovMKAvVOnfZyDvrVtL7x59AVj7/hoG9B3E/oN70RcYyLtawJ9/ZPPUk6MIa9iEa9dymL1wFgMHD+DrL/ZiY2PD2yvf4ottn7Ny6Rrq+NflwsX/4/8u/B95Vwv4zxcfs+ydpaxYkkDDBo3I/D2TE6dSyLtaAIBqKEqs8q4W8HCvR3nq2Vgu/3YFZ2cXAL7673au5eTQrWM0eVcLmDX3Nb7Y/h9mvzafoKD6HPx2P0+NjMPNyYOo+9qRn1PUr+5aoXEdx384yvffH6Pfw48Zy955dxlvvb2EuW8spGmTCDZ88D6PPdGPPV99S3BQfXKu5dDv8b60a/MAS+evJP3cr7w8bZKx7/ycQh57dBDrkt7jqdhnjds1cc0a7rs3Cj9vLbnX11WWbw8epqCggKhWDxjbebh606hBY/bt3U/bVh2M22fh4vnMmRePv19dYqIf4tmnnsPe3r7cvgsKCynI0/Pzd3+Ql2Wo0OdBiGIOznaSzNxMjRo18PX1rXD7/Px88vPzjfPZ2dl3Iiz8QzxoFRN0R/oW1UO5+YGtnhpOuTi62aOx1wBgyM3lfPcOdy22kuomH8DG0aFCbWvY21DDzoahTw5h1twZXMq6iKIoHD7yLe+/l8ShI/upYW+Ds4eGAY/3N1l2RcOVBIVoSc84S5PGTfjt0kVCQkJ4sEsHFEWhUXiose3vf2Tg6+NDj17dsLOzoyEh3N8hyliv2IDGsQbOHhp69e6J04vO/Peb7Tz+2CAAPt32ET16ROOr9SYnJ4cVq5fx+afbad3qPgCaNGvE0R8OsX7zWrp074RTVtFOvtejXbGxsUGn01FQUMCTscN4cliscb3vrF7K+HETeGJI0Xoi7pnNgcP7SEhawcJ5i9n08fuoBgOrVr6Lk5MT97SK5HJWJuPGj8HR1R5nDw3DhsUxZ8EsTv3yAy3vuZeCggI++uQD3nhtFs7umptu/+yrf2Bvb0+dANPvRF8/H/7IumRc/tlRo4mIiKSme02+O3qYV197hYu/nWfZknfK7TtfB/ZZNWjWsS7oLeNHo7AefvXdb93oDrL4ZObnn3/G398fjUZD69atmTVrFsHB5Q/DxsfHlzo0dSf4hXjgF+Jxx9cjqp+8vDxSU1NxcrXHweF6MmOnv8VSd46zmz02TjffiRarYW+LrZ0NgcF1iI6O5sOtG1FVlejoaAKD62BrZ0MNe1uc3TWcOXOGadOmcfDgQS5duoTBUPRr/9KVDJzdWzDi6WF06dKFFq2a0b17d3r16kXXrkWHawYNeZy3VyylafMwunfvTs+ePYmJiaFGjaKvLMVGwd6xxvWdt4b+/fuxeesHDH86jpycHD7/4jPWr1+Ps7uGkz/9QF5eHr0fiTZ5LzqdjubNm+PsrsHJrej9b9q0ibCwMAoKCkhJSWHs2LHU9qnFm2++SXZ2NhcvXqBjp/YmScf9D7Tj+++/x9ldw9m0X4iIjMDbr6axvsOD9wMUJTPuGuq71yM6OpoNH66jfad2bN36Bfn5eTwxdCBOt/h30DgXPWzS2cO0nY2tgr1DDWP5pCkvGOtat2uJb53aPProo8xfOA8vL68y+7bNU9E41qBRh7o4OFQsuRXCUlh0MtO6dWvee+89GjRowG+//cbMmTOJiorixIkT5f5BTp48mfHjxxvns7Oz0Wq1dytkIW6L4uhIw6NHzLbu2xEXF8fo0aMBWLZsWan6mJgYtFotq1atwt/fH4PBQHh4ODqdDoAWLVqQmprKtm3b2LlzJ/3796dz585s3rwZrVbL6dOn2bFjBzt37mTUqFHMnTuX3bt3l/n06EGDBtG+fXsyMzPZsWMHDg4O9OjRA8CYRH3++efUqVPHZDmNxjQp0Gq1hISEABAWFsbZs2eZNm0a06dP/3t73TDkpqqqsayiVz4NHz6cwYMHs3DhQhITExkwYECFTmj29fVFp9Nx5coVatb8O2HKzMwkKiqq3OXuu69oROqXX34p97tTCGtm0clM8ZcRQNOmTWnTpg3169dn7dq1JglLSRqNptQXlBCWTlEUlDtwdc6d1L17d2Ni0q1bN5O6y5cvc+rUKVasWMH99xeNTOzdu7dUH25ubgwYMIABAwbw6KOP0r17d/744w88PT1xdHSkd+/e9O7dm2effZZGjRqRkpJCixYtSvUTFRWFVqtl06ZNbNu2jX79+hnPD2ncuDEajYb09HTat29fqfdoa2tLYWEhOp0ONzc3/P392bt3Lw888ICxzf79+2nVqpVxXe+//z65ubk4Xk8SDx48WKrfnj174uzszPLly9m2bRt79uypUDz33HMPdnZ27Nixg/79iw7jXbx4kR9//JE5c+aUu9yxY8cA8PPzq9gbF8LKWHQycyNnZ2eaNm3Kzz//bO5QhPjXs7W15dSpU8bXJdWsWRMvLy9WrlyJn58f6enppe5zsnDhQvz8/IiMjMTGxoYPP/wQX19fPDw8WLNmDXq9ntatW+Pk5MT777+Po6MjgYGBZcaiKAoDBw7knXfe4aeffmLXrl3GOldXVyZOnMjzzz+PwWCgXbt2ZGdns3//flxcXBg6dKix7eXLl8nIyKCwsJCUlBQWL15Mx44dcXNzA+CFF17g1VdfpX79+kRGRpKYmMjx48dJSkoCYODAgbz88ssMGzaMqVOnkpaWxrx588rcdrGxsUyePJmQkBDatGlToW3u7u7OsGHDmDBhAl5eXnh6ejJx4kSaNm1K586dAThw4AAHDx6kY8eOuLu7c/jwYZ5//nl69+591y6dF+KuU61IXl6eWqdOHXXGjBkVXiYrK0sF1KysrDsYmRAVl5ubq548eVLNzc01dyiVNnToUPWhhx4qt/6hhx5Shw4dqqqqqu7YsUMNCwtTNRqN2qxZMzU5OVkF1K1bt6qqqqorV65UIyMjVWdnZ9XNzU3t1KmTevToUVVVVXXr1q1q69atVTc3N9XZ2Vm977771J07dxrXExgYqC5cuNBk3SdOnFABNTAwUDUYDCZ1BoNBXbx4sdqwYUPVzs5O9fb2Vrt166bu3r1bVVVVTU1NVSm6KlYFVFtbW7Vu3brqiBEj1MzMTGM/er1enTFjhlqnTh3Vzs5OjYiIULdt22ayrgMHDqgRERGqvb29GhkZqX700UcqoB47dsyk3ZkzZ1RAnTNnzq02u4nc3Fx19OjRqqenp+ro6Kj26tVLTU9PN9YfOXJEbd26teru7q46ODioDRs2VF999VU1Jyfnlv1a6+dSVE+V2X8rqmpht7csYeLEicTExBAQEEBmZiYzZ85k9+7dpKSklPsL7UbZ2dm4u7uTlZVl/HUlhDkVnwAcFBQkJ1r+i+3bt48OHTpw/vx5fHx8zB2OfC6FxanM/tuiDzOdP3+exx9/nEuXLuHt7c19993HwYMHK5zICCGEpcnPz+fcuXNMmzaN/v37W0QiI4S1s+hkZuPGjeYOQQghqtSGDRsYNmwYkZGRvP/++yZ1SUlJPP3002UuFxgYyIkTJ+5GiEJYHYtOZoQQorqJjY0t92GcvXv3pnXr1mXWlXVJuhCiiCQzQghhIVxdXXF1dTV3GEJYHXmamBBmYsHn3ot/Ifk8CmsmyYwQd1nx4YJr166ZORIh/lb8eZTDWcIayWEmIe4yW1tbPDw8yMzMBMDJyanCT60Woqqpqsq1a9fIzMzEw8Oj1A0QhbAGkswIYQbFT4IvTmiEMDcPDw/j51IIayPJjBBmoCgKfn5+1K5dm4KCAnOHI/7l7OzsZERGWDVJZoQwI1tbW9mJCCHEPyQnAAshhBDCqkkyI4QQQgirJsmMEEIIIayaJDNCCCGEsGqSzAghhBDCqkkyI4QQQgirJsmMEEIIIayaJDNCCCGEsGqSzAghhBDCqkkyI4QQQgirJsmMEEIIIayaJDNCCCGEsGqSzAghhBDCqslTs2/TH++9z6UVK8wdhhBCCGF2noMHU2vk02ZbvyQzt8mQl4f+8mVzhyGEEEKYneHaNbOuX5KZ22Tr6YlT69bmDkMIAaCYO4B/AdXcAQhLVqN2bfOu36xrt2L6P/7g2rffmjsMIYQQwuwcIyLMun5JZm7T783q8n+jeps7DCGEEMLsCiPqYc6xGUlmbtMau7/YVM/b3GEIIYQQZtfLIYelZly/JDO3KdfWnzzXBuYOQwghhDC7nBrZZl2/JDO3yeWCG057fjV3GEIIIYTZubWuC/eYb/2SzNymuvY1MMg9B4UQQgjq2Js3nZBk5jbVr3WCIfftMHcYQgghhNmF1b4faGq29Usyc5v2ff8hWzll7jCEEEIIs8u7/H/E3DvYbOu3imTm7bffZu7cuVy8eJEmTZqwaNEi7r//frPGVMcugOjfJZkRQggh6ngEmnX9iqqqFn1fx02bNjF48GDefvtt2rZty4oVK3j33Xc5efIkAQEBt1w+Ozsbd3d3srKycHNzuwsRCyGEEOKfqsz+2+KTmdatW9OiRQuWL19uLAsLC+Phhx8mPj6+VPv8/Hzy8/ON89nZ2Wi1WklmhBBCCCtSmWTGoi/H0el0HDlyhK5du5qUd+3alf3795e5THx8PO7u7sZJq9XejVCFEEIIYSYWncxcunQJvV6Pj4+PSbmPjw8ZGRllLjN58mSysrKM07lz5+5GqEIIIYQwE6s4AVhRTB+Jq6pqqbJiGo0GjUZzN8ISQgghhAWw6JGZWrVqYWtrW2oUJjMzs9RojRBCCCH+nSw6mbG3t+eee+5hxw7Tm9Pt2LGDqKgoM0UlhBBCCEti8YeZxo8fz+DBg2nZsiVt2rRh5cqVpKenM3LkSHOHJoQQQggLYPHJzIABA7h8+TKvvfYaFy9eJDw8nC+++ILAQPPeoEcIIYQQlsHi7zPzT8lN84QQQgjrU23uMyOEEEIIcSuSzAghhBDCqkkyI4QQQgirJsmMEEIIIayaJDNCCCGEsGoWf2n2P1V8sVZ2draZIxFCCCFERRXvtyty0XW1T2b++usvAHl6thBCCGGF/vrrL9zd3W/aptrfZ8ZgMHDhwgVcXV3LfTjl7crOzkar1XLu3Dm5h80dJNv57pDtfHfIdr47ZDvfHXdyO6uqyl9//YW/vz82Njc/K6baj8zY2NhQt27dO7oONzc3+WO5C2Q73x2yne8O2c53h2znu+NObedbjcgUkxOAhRBCCGHVJJkRQgghhFWTZOYf0Gg0vPrqq2g0GnOHUq3Jdr47ZDvfHbKd7w7ZzneHpWznan8CsBBCCCGqNxmZEUIIIYRVk2RGCCGEEFZNkhkhhBBCWDVJZoQQQghh1SSZuU1vv/02QUFBODg4cM899/DNN9+YO6RqJT4+nnvvvRdXV1dq167Nww8/zOnTp80dVrUXHx+PoiiMGzfO3KFUS//3f//HE088gZeXF05OTkRGRnLkyBFzh1WtFBYWMnXqVIKCgnB0dCQ4OJjXXnsNg8Fg7tCs2p49e4iJicHf3x9FUfj4449N6lVVZfr06fj7++Po6EiHDh04ceLEXYtPkpnbsGnTJsaNG8fLL7/MsWPHuP/+++nRowfp6enmDq3a2L17N88++ywHDx5kx44dFBYW0rVrV3JycswdWrV1+PBhVq5cSbNmzcwdSrV05coV2rZti52dHdu2bePkyZPMnz8fDw8Pc4dWrcyePZt33nmHpUuXcurUKebMmcPcuXN56623zB2aVcvJySEiIoKlS5eWWT9nzhwWLFjA0qVLOXz4ML6+vnTp0sX4fMQ7ThWV1qpVK3XkyJEmZY0aNVJfeuklM0VU/WVmZqqAunv3bnOHUi399ddfamhoqLpjxw61ffv26nPPPWfukKqdSZMmqe3atTN3GNVedHS0GhcXZ1LWp08f9YknnjBTRNUPoG7dutU4bzAYVF9fX/XNN980luXl5anu7u7qO++8c1dikpGZStLpdBw5coSuXbualHft2pX9+/ebKarqLysrCwBPT08zR1I9Pfvss0RHR9O5c2dzh1Jtffrpp7Rs2ZJ+/fpRu3ZtmjdvzqpVq8wdVrXTrl07/vvf//LTTz8B8P3337N371569uxp5siqr9TUVDIyMkz2ixqNhvbt29+1/WK1f9BkVbt06RJ6vR4fHx+Tch8fHzIyMswUVfWmqirjx4+nXbt2hIeHmzucamfjxo0cPXqUw4cPmzuUau3s2bMsX76c8ePHM2XKFA4dOsTYsWPRaDQMGTLE3OFVG5MmTSIrK4tGjRpha2uLXq/njTfe4PHHHzd3aNVW8b6vrP3ir7/+eldikGTmNimKYjKvqmqpMlE1Ro8ezQ8//MDevXvNHUq1c+7cOZ577jm++uorHBwczB1OtWYwGGjZsiWzZs0CoHnz5pw4cYLly5dLMlOFNm3axLp161i/fj1NmjTh+PHjjBs3Dn9/f4YOHWru8Ko1c+4XJZmppFq1amFra1tqFCYzM7NUVir+uTFjxvDpp5+yZ88e6tata+5wqp0jR46QmZnJPffcYyzT6/Xs2bOHpUuXkp+fj62trRkjrD78/Pxo3LixSVlYWBgfffSRmSKqnl544QVeeuklHnvsMQCaNm3Kr7/+Snx8vCQzd4ivry9QNELj5+dnLL+b+0U5Z6aS7O3tueeee9ixY4dJ+Y4dO4iKijJTVNWPqqqMHj2aLVu28PXXXxMUFGTukKqlTp06kZKSwvHjx41Ty5YtGTRoEMePH5dEpgq1bdu21O0FfvrpJwIDA80UUfV07do1bGxMd222trZyafYdFBQUhK+vr8l+UafTsXv37ru2X5SRmdswfvx4Bg8eTMuWLWnTpg0rV64kPT2dkSNHmju0auPZZ59l/fr1fPLJJ7i6uhpHwtzd3XF0dDRzdNWHq6trqfOQnJ2d8fLykvOTqtjzzz9PVFQUs2bNon///hw6dIiVK1eycuVKc4dWrcTExPDGG28QEBBAkyZNOHbsGAsWLCAuLs7coVm1q1ev8ssvvxjnU1NTOX78OJ6engQEBDBu3DhmzZpFaGgooaGhzJo1CycnJwYOHHh3Arwr10xVQ8uWLVMDAwNVe3t7tUWLFnLJcBUDypwSExPNHVq1J5dm3zn/+c9/1PDwcFWj0aiNGjVSV65cae6Qqp3s7Gz1ueeeUwMCAlQHBwc1ODhYffnll9X8/Hxzh2bVdu3aVeZ38tChQ1VVLbo8+9VXX1V9fX1VjUajPvDAA2pKSspdi09RVVW9O2mTEEIIIUTVk3NmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBAWKTk5GUVR+PPPP80dihDCwskdgIUQFqFDhw5ERkayaNEioOhBdX/88Qc+Pj4oimLe4IQQFk0eNCmEsEj29vb4+vqaOwwhhBWQw0xCCLOLjY1l9+7dLF68GEVRUBSFNWvWmBxmWrNmDR4eHnz22Wc0bNgQJycnHn30UXJycli7di316tWjZs2ajBkzBr1eb+xbp9Px4osvUqdOHZydnWndujXJycnmeaNCiDtCRmaEEGa3ePFifvrpJ8LDw3nttdcAOHHiRKl2165dY8mSJWzcuJG//vqLPn360KdPHzw8PPjiiy84e/Ysffv2pV27dgwYMACAJ598krS0NDZu3Ii/vz9bt26le/fupKSkEBoaelffpxDizpBkRghhdu7u7tjb2+Pk5GQ8tPS///2vVLuCggKWL19O/fr1AXj00Ud5//33+e2333BxcaFx48Z07NiRXbt2MWDAAM6cOcOGDRs4f/48/v7+AEycOJHt27eTmJjIrFmz7t6bFELcMZLMCCGshpOTkzGRAfDx8aFevXq4uLiYlGVmZgJw9OhRVFWlQYMGJv3k5+fj5eV1d4IWQtxxkswIIayGnZ2dybyiKGWWGQwGAAwGA7a2thw5cgRbW1uTdiUTICGEdZNkRghhEezt7U1O3K0KzZs3R6/Xk5mZyf3331+lfQshLIdczSSEsAj16tXj22+/JS0tjUuXLhlHV/6JBg0aMGjQIIYMGcKWLVtITU3l8OHDzJ49my+++KIKohZCWAJJZoQQFmHixInY2trSuHFjvL29SU9Pr5J+ExMTGTJkCBMmTKBhw4b07t2bb7/9Fq1WWyX9CyHMT+4ALIQQQgirJiMzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqkswIIcT/t1sHJAAAAACC/r9uR6ArBNZkBgBYkxkAYE1mAIA1mQEA1mQGAFgLiLFkSit7X7AAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "sim.data.where(sim.data['particle_type'] == 'Massive Body',drop=True)['a'].plot(hue=\"name\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f46b22e-dbd0-4562-a4cf-ea55247a2126", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python (My debug_env Kernel)", + "language": "python", + "name": "debug_env" + }, + "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.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/Basic_Simulation/test_io.ipynb b/examples/Basic_Simulation/test_io.ipynb deleted file mode 100644 index ce92a8229..000000000 --- a/examples/Basic_Simulation/test_io.ipynb +++ /dev/null @@ -1,143 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "86c845ce-1801-46ca-8a8a-1cabb266e6a6", - "metadata": {}, - "outputs": [], - "source": [ - "import swiftest\n", - "import xarray as xr\n", - "import numpy as np\n", - "import os" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "d716c371-8eb4-4fc1-82af-8b5c444c831e", - "metadata": {}, - "outputs": [], - "source": [ - "sim = swiftest.Simulation()" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "83cebbc1-387b-4ef5-b96e-76856b6672e5", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "integrator symba\n", - "codename Swiftest\n", - "driver_executable /Users/daminton/git/swiftest/build/swiftest_driver\n", - "t0 0.0 y\n", - "tstart 0.0 y\n", - "tstop NOT SET\n", - "dt NOT SET\n", - "istep_out NOT SET\n", - "istep_dump NOT SET\n", - "init_cond_file_type NETCDF_DOUBLE\n", - "init_cond_format EL\n", - "init_cond_file_name init_cond.nc\n", - "output_file_type NETCDF_DOUBLE\n", - "output_file_name bin.nc\n", - "output_format XVEL\n", - "rmin 0.004650467260962157 AU\n", - "rmax 10000.0 AU\n", - "qmin_coord HELIO\n", - "MU: MSun 1.988409870698051e+30 kg / MSun\n", - "DU: AU 149597870700.0 m / AU\n", - "TU: y 31557600.0 s / y\n", - "close_encounter_check True\n", - "general_relativity True\n", - "fragmentation True\n", - "rotation True\n", - "compute_conservation_values False\n", - "extra_force False\n", - "big_discard False\n", - "rhill_present False\n", - "restart False\n", - "interaction_loops TRIANGULAR\n", - "encounter_check_loops TRIANGULAR\n", - "ephemeris_date 2027-04-30\n" - ] - }, - { - "data": { - "text/plain": [ - "{'! VERSION': 'Swiftest parameter input',\n", - " 'T0': 0.0,\n", - " 'TSTART': 0.0,\n", - " 'IN_TYPE': 'NETCDF_DOUBLE',\n", - " 'IN_FORM': 'EL',\n", - " 'NC_IN': 'init_cond.nc',\n", - " 'OUT_TYPE': 'NETCDF_DOUBLE',\n", - " 'BIN_OUT': 'bin.nc',\n", - " 'OUT_FORM': 'XVEL',\n", - " 'CHK_RMIN': 0.004650467260962157,\n", - " 'CHK_RMAX': 10000.0,\n", - " 'CHK_QMIN_COORD': 'HELIO',\n", - " 'CHK_QMIN': 0.004650467260962157,\n", - " 'CHK_QMIN_RANGE': '0.004650467260962157 10000.0',\n", - " 'MU2KG': 1.988409870698051e+30,\n", - " 'DU2M': 149597870700.0,\n", - " 'TU2S': 31557600.0,\n", - " 'CHK_CLOSE': True,\n", - " 'GR': True,\n", - " 'FRAGMENTATION': True,\n", - " 'ROTATION': True,\n", - " 'ENERGY': False,\n", - " 'EXTRA_FORCE': False,\n", - " 'BIG_DISCARD': False,\n", - " 'RHILL_PRESENT': False,\n", - " 'RESTART': False,\n", - " 'INTERACTION_LOOPS': 'TRIANGULAR',\n", - " 'ENCOUNTER_CHECK': 'TRIANGULAR'}" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "sim.get_parameter()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ec7452d6-4c9b-4df3-acc0-b11c32264b91", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.13" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 2ec86f108..2b8c47cb3 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -852,6 +852,10 @@ def swiftest2xr(param, verbose=True): ds = fix_types(ds,ftype=np.float64) elif param['OUT_TYPE'] == "NETCDF_FLOAT": ds = fix_types(ds,ftype=np.float32) + # Check if the name variable contains unique values. If so, make name the dimension instead of id + if len(np.unique(ds['name'])) == len(ds['name']): + ds = ds.swap_dims({"id" : "name"}) + ds = ds.reset_coords("id") else: print(f"Error encountered. OUT_TYPE {param['OUT_TYPE']} not recognized.") return None diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 4dbc05610..7c17c3d32 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -23,6 +23,7 @@ import shutil import subprocess import shlex +import warnings from typing import ( Literal, Dict, @@ -80,7 +81,7 @@ def __init__(self,read_param: bool = True, **kwargs: Any): The stopping time for a simulation. `tstop` must be greater than `tstart`. Parameter input file equivalent: `TSTOP` dt : float, optional - The step size of the simulation. `dt` must be less than or equal to `tstop-dstart`. + The step size of the simulation. `dt` must be less than or equal to `tstop-tstart`. Parameter input file equivalent: `DT` istep_out : int, optional The number of time steps between outputs to file. *Note*: only `istep_out` or `toutput` can be set. @@ -275,7 +276,7 @@ def __init__(self,read_param: bool = True, **kwargs: Any): self._getter_column_width = '32' self.param = {} - self.ds = xr.Dataset() + self.data = xr.Dataset() # Parameters are set in reverse priority order. First the defaults, then values from a pre-existing input file, # then using the arguments passed via **kwargs. @@ -316,7 +317,7 @@ def __init__(self,read_param: bool = True, **kwargs: Any): # Let the user know that there was a problem reading an old parameter file and we're going to create a new one if read_param and not param_file_found: - print(f"{self.param_file} not found. Creating a new file using default values for parameters not passed to Simulation().") + warnings.warn(f"{self.param_file} not found. Creating a new file using default values for parameters not passed to Simulation().") self.write_param() # Read in an old simulation file if requested @@ -326,7 +327,7 @@ def __init__(self,read_param: bool = True, **kwargs: Any): if os.path.exists(binpath): self.bin2xr() else: - print(f"BIN_OUT file {binpath} not found.") + warnings.warn(f"BIN_OUT file {binpath} not found.") return @@ -352,35 +353,46 @@ def run(self,**kwargs): self.write_param() if self.codename != "Swiftest": - print(f"Running an integration is not yet supported for {self.codename}") + warnings.warn(f"Running an integration is not yet supported for {self.codename}") return if self.driver_executable is None: - print("Path to swiftest_driver has not been set!") - print(f"Make sure swiftest_driver is compiled and the executable is in {self.binary_path}") + warnings.warn("Path to swiftest_driver has not been set!") + warnings.warn(f"Make sure swiftest_driver is compiled and the executable is in {self.binary_path}") return print(f"Running a {self.codename} {self.integrator} run from tstart={self.param['TSTART']} {self.TU_name} to tstop={self.param['TSTOP']} {self.TU_name}") # Get current environment variables env = os.environ.copy() + driver_script = os.path.join(self.binary_path,"swiftest_driver.sh") + shell = os.path.basename(env['SHELL']) + with open(driver_script,'w') as f: + f.write(f"#{env['SHELL']} -l {os.linesep}") + f.write(f"source ~/.{shell}rc {os.linesep}") + f.write(f"cd {self.sim_dir} {os.linesep}") + f.write(f"{self.driver_executable} {self.integrator} {self.param_file}") try: - cmd = f"{self.driver_executable} {self.integrator} {self.param_file}" - p = subprocess.Popen(shlex.split(cmd), + cmd = f"{env['SHELL']} -l {driver_script}" + with subprocess.Popen(shlex.split(cmd), stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=env, - universal_newlines=True) - for line in p.stdout: - print(line, end='') - res = p.communicate() - if p.returncode != 0: - for line in res[1]: - print(line, end='') - raise Exception ("Failure in swiftest_driver") + universal_newlines=True) as p: + for line in p.stdout: + print(line, end='') + res = p.communicate() + if p.returncode != 0: + for line in res[1]: + print(line, end='') + raise Exception ("Failure in swiftest_driver") except: - print(f"Error executing main swiftest_driver program") + warnings.warn(f"Error executing main swiftest_driver program") + return + + # Read in new data + self.bin2xr() return @@ -499,7 +511,7 @@ def set_simulation_time(self, if tstop is not None: if tstop <= tstart: - print("Error! tstop must be greater than tstart.") + warnings.warn("tstop must be greater than tstart.") return {} if tstop is not None: @@ -512,8 +524,8 @@ def set_simulation_time(self, if dt is not None and tstop is not None: if dt > (tstop - tstart): - print("Error! dt must be smaller than tstop-tstart") - print(f"Setting dt = {tstop - tstart} instead of {dt}") + warnings.warn("dt must be smaller than tstop-tstart") + warnings.warn(f"Setting dt = {tstop - tstart} instead of {dt}") dt = tstop - tstart if dt is not None: @@ -521,26 +533,27 @@ def set_simulation_time(self, if istep_out is None and tstep_out is None: istep_out = self.param.pop("ISTEP_OUT", None) - - if istep_out is not None and tstep_out is not None: - print("Error! istep_out and tstep_out cannot both be set") + elif istep_out is not None and tstep_out is not None: + warnings.warn("istep_out and tstep_out cannot both be set") return {} + else: + update_list.append("istep_out") if tstep_out is not None and dt is not None: istep_out = int(np.floor(tstep_out / dt)) if istep_out is not None: self.param['ISTEP_OUT'] = istep_out - update_list.append("istep_out") if istep_dump is None: istep_dump = self.param.pop("ISTEP_DUMP", None) if istep_dump is None: istep_dump = istep_out + else: + update_list.append("istep_dump") if istep_dump is not None: self.param['ISTEP_DUMP'] = istep_dump - update_list.append("istep_dump") time_dict = self.get_simulation_time(update_list, verbose=verbose) @@ -769,7 +782,7 @@ def set_integrator(self, if codename is not None: valid_codename = ["Swiftest", "Swifter", "Swift"] if codename.title() not in valid_codename: - print(f"{codename} is not a valid codename. Valid options are ",",".join(valid_codename)) + warnings.warn(f"{codename} is not a valid codename. Valid options are ",",".join(valid_codename)) try: self.codename except: @@ -783,7 +796,7 @@ def set_integrator(self, self.binary_path = os.path.realpath(os.path.join(os.path.dirname(os.path.realpath(_pyfile)),os.pardir,os.pardir,os.pardir,"bin")) self.driver_executable = os.path.join(self.binary_path,"swiftest_driver") if not os.path.exists(self.driver_executable): - print(f"Cannot find the Swiftest driver in {self.binary_path}") + warnings.warn(f"Cannot find the Swiftest driver in {self.binary_path}") self.driver_executable = None else: self.binary_path = "NOT IMPLEMENTED FOR THIS CODE" @@ -793,7 +806,7 @@ def set_integrator(self, if integrator is not None: valid_integrator = ["symba","rmvs","whm","helio"] if integrator.lower() not in valid_integrator: - print(f"{integrator} is not a valid integrator. Valid options are ",",".join(valid_integrator)) + warnings.warn(f"{integrator} is not a valid integrator. Valid options are ",",".join(valid_integrator)) try: self.integrator except: @@ -804,9 +817,9 @@ def set_integrator(self, if mtiny is not None or gmtiny is not None: if self.integrator != "symba": - print("mtiny and gmtiny are only used by SyMBA.") + warnings.warn("mtiny and gmtiny are only used by SyMBA.") if mtiny is not None and gmtiny is not None: - print("Only set mtiny or gmtiny, not both!") + warnings.warn("Only set mtiny or gmtiny, not both!") elif gmtiny is not None: self.param['GMTINY'] = gmtiny update_list.append("gmtiny") @@ -851,13 +864,13 @@ def get_integrator(self,arg_list: str | List[str] | None = None, verbose: bool | try: self.integrator except: - print(f"integrator is not set") + warnings.warn(f"integrator is not set") return {} try: self.codename except: - print(f"codename is not set") + warnings.warn(f"codename is not set") return {} if verbose is None: @@ -994,19 +1007,19 @@ def set_feature(self, if fragmentation is not None: if self.codename != "Swiftest" and self.integrator != "symba" and fragmentation: - print("Fragmentation is only available on Swiftest SyMBA.") + warnings.warn("Fragmentation is only available on Swiftest SyMBA.") self.param['FRAGMENTATION'] = False else: self.param['FRAGMENTATION'] = fragmentation update_list.append("fragmentation") if fragmentation: if "MIN_GMFRAG" not in self.param and minimum_fragment_mass is None and minimum_fragment_gmass is None: - print("Minimum fragment mass is not set. Set it using minimum_fragment_gmass or minimum_fragment_mass") + warnings.warn("Minimum fragment mass is not set. Set it using minimum_fragment_gmass or minimum_fragment_mass") else: update_list.append("minimum_fragment_gmass") if minimum_fragment_gmass is not None and minimum_fragment_mass is not None: - print("Warning! Only set either minimum_fragment_mass or minimum_fragment_gmass, but not both!") + warnings.warn("Only set either minimum_fragment_mass or minimum_fragment_gmass, but not both!") if minimum_fragment_gmass is not None: self.param["MIN_GMFRAG"] = minimum_fragment_gmass @@ -1051,8 +1064,8 @@ def set_feature(self, if interaction_loops is not None: valid_vals = ["TRIANGULAR", "FLAT", "ADAPTIVE"] if interaction_loops not in valid_vals: - print(f"{interaction_loops} is not a valid option for interaction loops.") - print(f"Must be one of {valid_vals}") + warnings.warn(f"{interaction_loops} is not a valid option for interaction loops.") + warnings.warn(f"Must be one of {valid_vals}") if "INTERACTION_LOOPS" not in self.param: self.param["INTERACTION_LOOPS"] = valid_vals[0] else: @@ -1062,8 +1075,8 @@ def set_feature(self, if encounter_check_loops is not None: valid_vals = ["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] if encounter_check_loops not in valid_vals: - print(f"{encounter_check_loops} is not a valid option for interaction loops.") - print(f"Must be one of {valid_vals}") + warnings.warn(f"{encounter_check_loops} is not a valid option for interaction loops.") + warnings.warn(f"Must be one of {valid_vals}") if "ENCOUNTER_CHECK" not in self.param: self.param["ENCOUNTER_CHECK"] = valid_vals[0] else: @@ -1194,13 +1207,13 @@ def set_init_cond_files(self, return {} def ascii_file_input_error_msg(codename): - print(f"Error in set_init_cond_files: init_cond_file_name must be a dictionary of the form: ") - print('{') + warnings.warn(f"in set_init_cond_files: init_cond_file_name must be a dictionary of the form: ") + warnings.warn('{') if codename == "Swiftest": - print('"CB" : *path to central body initial conditions file*,') - print('"PL" : *path to massive body initial conditions file*,') - print('"TP" : *path to test particle initial conditions file*') - print('}') + warnings.warn('"CB" : *path to central body initial conditions file*,') + warnings.warn('"PL" : *path to massive body initial conditions file*,') + warnings.warn('"TP" : *path to test particle initial conditions file*') + warnings.warn('}') return {} if init_cond_format is None: @@ -1220,21 +1233,21 @@ def ascii_file_input_error_msg(codename): else: init_cond_keys = ["PL", "TP"] if init_cond_file_type != "ASCII": - print(f"{init_cond_file_type} is not supported by {self.codename}. Using ASCII instead") + warnings.warn(f"{init_cond_file_type} is not supported by {self.codename}. Using ASCII instead") init_cond_file_type = "ASCII" if init_cond_format != "XV": - print(f"{init_cond_format} is not supported by {self.codename}. Using XV instead") + warnings.warn(f"{init_cond_format} is not supported by {self.codename}. Using XV instead") init_cond_format = "XV" valid_formats = {"EL", "XV"} if init_cond_format not in valid_formats: - print(f"{init_cond_format} is not a valid input format") + warnings.warn(f"{init_cond_format} is not a valid input format") else: self.param['IN_FORM'] = init_cond_format valid_types = {"NETCDF_DOUBLE", "NETCDF_FLOAT", "ASCII"} if init_cond_file_type not in valid_types: - print(f"{init_cond_file_type} is not a valid input type") + warnings.warn(f"{init_cond_file_type} is not a valid input type") else: self.param['IN_TYPE'] = init_cond_file_type @@ -1261,7 +1274,7 @@ def ascii_file_input_error_msg(codename): elif type(init_cond_file_name) is dict: # Oops, accidentally passed a dictionary instead of the expected single string or path-like for NetCDF # input type. - print(f"Only a single input file is used for NetCDF files") + warnings.warn(f"Only a single input file is used for NetCDF files") else: self.param["NC_IN"] = init_cond_file_name @@ -1402,7 +1415,7 @@ def set_output_files(self, if output_file_type is None: output_file_type = "NETCDF_DOUBLE" elif output_file_type not in ["NETCDF_DOUBLE", "NETCDF_FLOAT"]: - print(f"{output_file_type} is not compatible with Swiftest. Setting to NETCDF_DOUBLE") + warnings.warn(f"{output_file_type} is not compatible with Swiftest. Setting to NETCDF_DOUBLE") output_file_type = "NETCDF_DOUBLE" elif self.codename == "Swifter": if output_file_type is None: @@ -1410,7 +1423,7 @@ def set_output_files(self, if output_file_type is None: output_file_type = "REAL8" elif output_file_type not in ["REAL4", "REAL8", "XDR4", "XDR8"]: - print(f"{output_file_type} is not compatible with Swifter. Setting to REAL8") + warnings.warn(f"{output_file_type} is not compatible with Swifter. Setting to REAL8") output_file_type = "REAL8" elif self.codename == "Swift": if output_file_type is None: @@ -1418,7 +1431,7 @@ def set_output_files(self, if output_file_type is None: output_file_type = "REAL4" if output_file_type not in ["REAL4"]: - print(f"{output_file_type} is not compatible with Swift. Setting to REAL4") + warnings.warn(f"{output_file_type} is not compatible with Swift. Setting to REAL4") output_file_type = "REAL4" self.param['OUT_TYPE'] = output_file_type @@ -1431,7 +1444,7 @@ def set_output_files(self, self.param['BIN_OUT'] = output_file_name if output_format != "XV" and self.codename != "Swiftest": - print(f"{output_format} is not compatible with {self.codename}. Setting to XV") + warnings.warn(f"{output_format} is not compatible with {self.codename}. Setting to XV") output_format = "XV" self.param["OUT_FORM"] = output_format @@ -1606,7 +1619,7 @@ def set_unit_system(self, self.param['MU2KG'] = 1000.0 self.MU_name = "g" else: - print(f"{MU} not a recognized unit system. Using MSun as a default.") + warnings.warn(f"{MU} not a recognized unit system. Using MSun as a default.") self.param['MU2KG'] = constants.MSun self.MU_name = "MSun" @@ -1629,7 +1642,7 @@ def set_unit_system(self, self.param['DU2M'] = 100.0 self.DU_name = "cm" else: - print(f"{DU} not a recognized unit system. Using AU as a default.") + warnings.warn(f"{DU} not a recognized unit system. Using AU as a default.") self.param['DU2M'] = constants.AU2M self.DU_name = "AU" @@ -1649,7 +1662,7 @@ def set_unit_system(self, self.param['TU2S'] = 1.0 self.TU_name = "s" else: - print(f"{TU} not a recognized unit system. Using YR as a default.") + warnings.warn(f"{TU} not a recognized unit system. Using YR as a default.") self.param['TU2S'] = constants.YR2S self.TU_name = "y" @@ -1832,7 +1845,7 @@ def set_distance_range(self, if qmin_coord is not None: valid_qmin_coord = ["HELIO","BARY"] if qmin_coord.upper() not in valid_qmin_coord: - print(f"qmin_coord = {qmin_coord} is not a valid option. Must be one of",','.join(valid_qmin_coord)) + warnings.warn(f"qmin_coord = {qmin_coord} is not a valid option. Must be one of",','.join(valid_qmin_coord)) self.param['CHK_QMIN_COORD'] = valid_qmin_coord[0] else: self.param['CHK_QMIN_COORD'] = qmin_coord.upper() @@ -1946,7 +1959,7 @@ def add_solar_system_body(self, >*Note.* Currently only the JPL Horizons ephemeris is implemented, so this is ignored. Returns ------- - ds : Xarray dataset with body or bodies added. + data : Xarray dataset with body or bodies added. """ if type(name) is str: @@ -1955,7 +1968,7 @@ def add_solar_system_body(self, if type(ephemeris_id) is int: ephemeris_id = [ephemeris_id] if len(ephemeris_id) != len(name): - print(f"Error! The length of ephemeris_id ({len(ephemeris_id)}) does not match the length of name ({len(name)})") + warnings.warn(f"The length of ephemeris_id ({len(ephemeris_id)}) does not match the length of name ({len(name)})") return None else: ephemeris_id = [None] * len(name) @@ -1968,11 +1981,11 @@ def add_solar_system_body(self, try: datetime.datetime.fromisoformat(date) except: - print(f"{date} is not a valid date format. Must be 'YYYY-MM-DD'. Setting to {self.ephemeris_date}") + warnings.warn(f"{date} is not a valid date format. Must be 'YYYY-MM-DD'. Setting to {self.ephemeris_date}") date = self.ephemeris_date if source.upper() != "HORIZONS": - print("Currently only the JPL Horizons ephemeris service is supported") + warnings.warn("Currently only the JPL Horizons ephemeris service is supported") body_list = [] for i,n in enumerate(name): @@ -2032,6 +2045,7 @@ def add_solar_system_body(self, J2=J2, J4=J4, t=t) dsnew = self._combine_and_fix_dsnew(dsnew) + self.save() return dsnew @@ -2076,8 +2090,8 @@ def set_ephemeris_date(self, datetime.datetime.fromisoformat(ephemeris_date) except: valid_date_args = ['"MBCL"', '"TODAY"', '"YYYY-MM-DD"'] - print(f"{ephemeris_date} is not a valid format. Valid options include:", ', '.join(valid_date_args)) - print("Using MBCL for date.") + warnings.warn(f"{ephemeris_date} is not a valid format. Valid options include:", ', '.join(valid_date_args)) + warnings.warn("Using MBCL for date.") ephemeris_date = minton_bcl self.ephemeris_date = ephemeris_date @@ -2112,7 +2126,7 @@ def get_ephemeris_date(self, arg_list: str | List[str] | None = None, verbose: b try: self.ephemeris_date except: - print(f"ephemeris_date is not set") + warnings.warn(f"ephemeris_date is not set") return valid_arg = {"ephemeris_date": self.ephemeris_date} @@ -2181,7 +2195,7 @@ def add_body(self, Adds a body (test particle or massive body) to the internal DataSet given a set up 6 vectors (orbital elements or cartesian state vectors, depending on the value of self.param). Input all angles in degress. - This method will update self.ds with the new body or bodies added to the existing Dataset. + This method will update self.data with the new body or bodies added to the existing Dataset. Parameters ---------- @@ -2214,7 +2228,7 @@ def add_body(self, Returns ------- - ds : Xarray Dataset + data : Xarray Dataset Dasaset containing the body or bodies that were added """ @@ -2262,16 +2276,16 @@ def input_to_array(val,t,n=None): J2 = input_to_array(J2,"f",nbodies) J4 = input_to_array(J4,"f",nbodies) - if len(self.ds) == 0: + if len(self.data) == 0: maxid = -1 else: - maxid = self.ds.id.max().values[()] + maxid = self.data.id.max().values[()] if idvals is None: idvals = np.arange(start=maxid+1,stop=maxid+1+nbodies,dtype=int) - if len(self.ds) > 0: - dup_id = np.in1d(idvals,self.ds.id) + if len(self.data) > 0: + dup_id = np.in1d(idvals, self.data.id) if any(dup_id): raise ValueError(f"Duplicate ids detected: ", *idvals[dup_id]) @@ -2284,6 +2298,7 @@ def input_to_array(val,t,n=None): J2=J2, J4=J4,t=t) dsnew = self._combine_and_fix_dsnew(dsnew) + self.save() return dsnew @@ -2302,7 +2317,7 @@ def _combine_and_fix_dsnew(self,dsnew): """ - self.ds = xr.combine_by_coords([self.ds, dsnew]) + self.data = xr.combine_by_coords([self.data, dsnew]) def get_nvals(ds): if "Gmass" in ds: @@ -2314,14 +2329,14 @@ def get_nvals(ds): return ds dsnew = get_nvals(dsnew) - self.ds = get_nvals(self.ds) + self.data = get_nvals(self.data) if self.param['OUT_TYPE'] == "NETCDF_DOUBLE": dsnew = io.fix_types(dsnew, ftype=np.float64) - self.ds = io.fix_types(self.ds, ftype=np.float64) + self.data = io.fix_types(self.data, ftype=np.float64) elif self.param['OUT_TYPE'] == "NETCDF_FLOAT": dsnew = io.fix_types(dsnew, ftype=np.float32) - self.ds = io.fix_types(self.ds, ftype=np.float32) + self.data = io.fix_types(self.data, ftype=np.float32) return dsnew @@ -2350,7 +2365,7 @@ def read_param(self, param_file, codename="Swiftest", verbose=True): self.param = io.read_swift_param(param_file, verbose=verbose) self.codename = "Swift" else: - print(f'{codename} is not a recognized code name. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn(f'{codename} is not a recognized code name. Valid options are "Swiftest", "Swifter", or "Swift".') self.codename = "Unknown" return @@ -2400,7 +2415,7 @@ def write_param(self, elif codename == "Swift": io.write_swift_param(param, param_file) else: - print( 'Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') return def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", tpname="tp.swiftest.in", @@ -2430,10 +2445,10 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t """ oldparam = self.param if self.codename == newcodename: - print(f"This parameter configuration is already in {newcodename} format") + warnings.warn(f"This parameter configuration is already in {newcodename} format") return oldparam if newcodename != "Swift" and newcodename != "Swifter" and newcodename != "Swiftest": - print(f'{newcodename} is an invalid code type. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn(f'{newcodename} is an invalid code type. Valid options are "Swiftest", "Swifter", or "Swift".') return oldparam goodconversion = True if self.codename == "Swifter": @@ -2454,7 +2469,7 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t if goodconversion: self.write_param(param_file) else: - print(f"Conversion from {self.codename} to {newcodename} is not supported.") + warnings.warn(f"Conversion from {self.codename} to {newcodename} is not supported.") return oldparam def bin2xr(self): @@ -2466,25 +2481,24 @@ def bin2xr(self): Returns ------- - self.ds : xarray dataset + self.data : xarray dataset """ # Make a temporary copy of the parameter dictionary so we can supply the absolute path of the binary file # This is done to handle cases where the method is called from a different working directory than the simulation # results param_tmp = self.param.copy() - param_tmp['BIN_OUT'] = os.path.join(self.dir_path, self.param['BIN_OUT']) + param_tmp['BIN_OUT'] = os.path.join(self.sim_dir, self.param['BIN_OUT']) if self.codename == "Swiftest": - self.ds = io.swiftest2xr(param_tmp, verbose=self.verbose) - if self.verbose: print('Swiftest simulation data stored as xarray DataSet .ds') + self.data = io.swiftest2xr(param_tmp, verbose=self.verbose) + if self.verbose: print('Swiftest simulation data stored as xarray DataSet .data') elif self.codename == "Swifter": - self.ds = io.swifter2xr(param_tmp, verbose=self.verbose) - if self.verbose: print('Swifter simulation data stored as xarray DataSet .ds') + self.data = io.swifter2xr(param_tmp, verbose=self.verbose) + if self.verbose: print('Swifter simulation data stored as xarray DataSet .data') elif self.codename == "Swift": - print("Reading Swift simulation data is not implemented yet") + warnings.warn("Reading Swift simulation data is not implemented yet") else: - print( - 'Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') return def follow(self, codestyle="Swifter"): @@ -2500,7 +2514,7 @@ def follow(self, codestyle="Swifter"): ------- fol : xarray dataset """ - if self.ds is None: + if self.data is None: self.bin2xr() if codestyle == "Swift": try: @@ -2515,10 +2529,10 @@ def follow(self, codestyle="Swifter"): i_list = [i for i in line.split(" ") if i.strip()] nskp = int(i_list[0]) except IOError: - print('No follow.in file found') + warnings.warn('No follow.in file found') ifol = None nskp = None - fol = tool.follow_swift(self.ds, ifol=ifol, nskp=nskp) + fol = tool.follow_swift(self.data, ifol=ifol, nskp=nskp) else: fol = None @@ -2560,17 +2574,17 @@ def save(self, param = self.param if codename == "Swiftest": - io.swiftest_xr2infile(ds=self.ds, param=param, in_type=self.param['IN_TYPE'], framenum=framenum) + io.swiftest_xr2infile(ds=self.data, param=param, in_type=self.param['IN_TYPE'], framenum=framenum) self.write_param(param_file=param_file) elif codename == "Swifter": if codename == "Swiftest": swifter_param = io.swiftest2swifter_param(param) else: swifter_param = param - io.swifter_xr2infile(self.ds, swifter_param, framenum) + io.swifter_xr2infile(self.data, swifter_param, framenum) self.write_param(param_file, param=swifter_param) else: - print(f'Saving to {codename} not supported') + warnings.warn(f'Saving to {codename} not supported') return @@ -2596,7 +2610,8 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil Returns ------- - frame : NetCDF dataset + frame : NetCDF dataset + A dataset containing the extracted initial condition data. """ if codename != "Swiftest": @@ -2608,13 +2623,13 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil if codename == "Swiftest": if restart: - new_param['T0'] = self.ds.time.values[framenum] + new_param['T0'] = self.data.time.values[framenum] if self.param['OUT_TYPE'] == 'NETCDF_DOUBLE': new_param['IN_TYPE'] = 'NETCDF_DOUBLE' elif self.param['OUT_TYPE'] == 'NETCDF_FLOAT': new_param['IN_TYPE'] = 'NETCDF_FLOAT' else: - print(f"{self.param['OUT_TYPE']} is an invalid OUT_TYPE file") + warnings.warn(f"{self.param['OUT_TYPE']} is an invalid OUT_TYPE file") return if self.param['BIN_OUT'] != new_param['BIN_OUT'] and restart: @@ -2631,7 +2646,7 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil new_param.pop('TP_IN', None) new_param.pop('CB_IN', None) print(f"Extracting data from dataset at time frame number {framenum} and saving it to {new_param['NC_IN']}") - frame = io.swiftest_xr2infile(self.ds, self.param, infile_name=new_param['NC_IN'], framenum=framenum) + frame = io.swiftest_xr2infile(self.data, self.param, infile_name=new_param['NC_IN'], framenum=framenum) print(f"Saving parameter configuration file to {new_param_file}") self.write_param(new_param_file, param=new_param) diff --git a/src/io/io.f90 b/src/io/io.f90 index aa797ddbe..86cf55728 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -566,6 +566,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 ! Internals logical :: t0_set = .false. !! Is the initial time set in the input file? + logical :: tstart_set = .false. !! Is the final time set in the input file? logical :: tstop_set = .false. !! Is the final time set in the input file? logical :: dt_set = .false. !! Is the step size set in the input file? integer(I4B) :: ilength, ifirst, ilast, i !! Variables used to parse input file @@ -593,6 +594,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) case ("T0") read(param_value, *, err = 667, iomsg = iomsg) param%t0 t0_set = .true. + case ("TSTART") + read(param_value, *, err = 667, iomsg = iomsg) param%t0 + tstart_set = .true. case ("TSTOP") read(param_value, *, err = 667, iomsg = iomsg) param%tstop tstop_set = .true. @@ -743,6 +747,12 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) read(param_value, *, err = 667, iomsg = iomsg) param%maxid_collision case ("PARTICLE_OUT") param%particle_out = param_value + case ("RESTART") + if (param_value == "NO" .or. param_value == 'F') then + param%lrestart = .false. + else if (param_value == "YES" .or. param_value == 'T') then + param%lrestart = .true. + end if case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(*,*) "Ignoring unknown parameter -> ",param_name @@ -929,7 +939,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = 0 ! Print the contents of the parameter file to standard output - call param%writer(unit = OUTPUT_UNIT, iotype = "none", v_list = [0], iostat = iostat, iomsg = iomsg) + ! call param%writer(unit = OUTPUT_UNIT, iotype = "none", v_list = [0], iostat = iostat, iomsg = iomsg) end associate diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 2cca37386..4c037a291 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -351,13 +351,15 @@ module subroutine netcdf_open(self, param, readonly) ! Internals integer(I4B) :: mode, status character(len=NF90_MAX_NAME) :: str_dim_name + character(len=STRMAX) :: errmsg mode = NF90_WRITE if (present(readonly)) then if (readonly) mode = NF90_NOWRITE end if - call check( nf90_open(param%outfile, mode, self%ncid), "netcdf_open nf90_open" ) + write(errmsg,*) "netcdf_open nf90_open ",trim(adjustl(param%outfile)) + call check( nf90_open(param%outfile, mode, self%ncid), errmsg) call check( nf90_inq_dimid(self%ncid, TIME_DIMNAME, self%time_dimid), "netcdf_open nf90_inq_dimid time_dimid" ) call check( nf90_inq_dimid(self%ncid, ID_DIMNAME, self%id_dimid), "netcdf_open nf90_inq_dimid id_dimid" )