{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "c5a9be36", "metadata": { "execution": { "iopub.execute_input": "2024-09-30T14:48:27.312422Z", "iopub.status.busy": "2024-09-30T14:48:27.312202Z", "iopub.status.idle": "2024-09-30T14:48:27.707950Z", "shell.execute_reply": "2024-09-30T14:48:27.707178Z" }, "nbsphinx": "hidden" }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "id": "69a01195", "metadata": { "execution": { "iopub.execute_input": "2024-09-30T14:48:27.710353Z", "iopub.status.busy": "2024-09-30T14:48:27.710046Z", "iopub.status.idle": "2024-09-30T14:48:27.718566Z", "shell.execute_reply": "2024-09-30T14:48:27.717985Z" }, "nbsphinx": "hidden" }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/runner/work/gala/gala/docs/tutorials/nb_setup:1: DeprecationWarning: `magic(...)` is deprecated since IPython 0.13 (warning added in 8.1), use run_line_magic(magic_name, parameter_s).\n", " get_ipython().magic('config InlineBackend.figure_format = \"retina\"') # noqa\n" ] } ], "source": [ "%run nb_setup" ] }, { "cell_type": "markdown", "id": "f1156cb4", "metadata": {}, "source": [ "# Defining a Milky Way potential model" ] }, { "cell_type": "code", "execution_count": 3, "id": "0a7d823b", "metadata": { "execution": { "iopub.execute_input": "2024-09-30T14:48:27.720558Z", "iopub.status.busy": "2024-09-30T14:48:27.720349Z", "iopub.status.idle": "2024-09-30T14:48:28.322989Z", "shell.execute_reply": "2024-09-30T14:48:28.322218Z" } }, "outputs": [], "source": [ "# Third-party dependencies\n", "import astropy.units as u\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from astropy.io import ascii\n", "from scipy.optimize import leastsq\n", "\n", "# Gala\n", "import gala.potential as gp\n", "from gala.units import galactic\n" ] }, { "cell_type": "markdown", "id": "ce029bfb", "metadata": {}, "source": [ "## Introduction\n", "\n", "`gala` provides a simple and easy way to access and integrate orbits in an\n", "approximate mass model for the Milky Way. The parameters of the mass model are\n", "determined by least-squares fitting the enclosed mass profile of a pre-defined\n", "potential form to recent measurements compiled from the literature. These\n", "measurements are provided with the documentation of `gala` and are shown\n", "below. The radius units are kpc, and mass units are solar masses:" ] }, { "cell_type": "code", "execution_count": 4, "id": "d4483974", "metadata": { "execution": { "iopub.execute_input": "2024-09-30T14:48:28.325721Z", "iopub.status.busy": "2024-09-30T14:48:28.325195Z", "iopub.status.idle": "2024-09-30T14:48:28.335526Z", "shell.execute_reply": "2024-09-30T14:48:28.334899Z" } }, "outputs": [], "source": [ "tbl = ascii.read(\"data/MW_mass_enclosed.csv\")" ] }, { "cell_type": "code", "execution_count": 5, "id": "9683e2c0", "metadata": { "execution": { "iopub.execute_input": "2024-09-30T14:48:28.337719Z", "iopub.status.busy": "2024-09-30T14:48:28.337474Z", "iopub.status.idle": "2024-09-30T14:48:28.344203Z", "shell.execute_reply": "2024-09-30T14:48:28.343662Z" } }, "outputs": [ { "data": { "text/html": [ "
Table length=16\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
rMencMenc_err_negMenc_err_posref
float64float64float64float64str27
0.0130000000.010000000.010000000.0Feldmeier et al. (2014)
0.12800000000.0200000000.0200000000.0Launhardt et al. (2002)
8.189502860861.524294994562473.7977144858963492.608627Bovy et al. (2012)
8.3110417867208.190554475949382.6968844387023236.020782McMillan (2011)
8.4102421035406.9035616733918715.62994415468328224.531876Koposov et al. (2010)
19.0208023299175.3043844317988008.3810134833267089.920685Kuepper et al. (2015)
50.0539884832748.4897519995734543.31433268490735257.4718Wilkinson & Evans (1999)
50.0529886965173.187269997867269.65930238536752776.21277Sakamoto et al. (2003)
50.0399914690706.92847109976539940.5371172696676468.2511Smith et al. (2007)
50.0419910425325.726839991469076.96850638172735113.64386Deason et al. (2012)
60.0399914690957.518869985070910.6335464344945146.92987Xue et al. (2008)
80.0689852841359.0248299936018002.3314110361048549.1029Gnedin et al. (2010)
100.01399701417307.7747899808054059.3811831336271726.1903Watkins et al. (2010)
120.0539884832260.29584199957345314.72906123854764645.32489Battaglia et al. (2005)
150.0750000000000.0250000000000.0250000000000.0Deason et al. (2012)
200.0679854974257.2006409912558030.05396313652012195.06256Bhattacherjee et al. (2014)
" ], "text/plain": [ "\n", " r Menc ... Menc_err_pos ref \n", "float64 float64 ... float64 str27 \n", "------- ------------------ ... ------------------ ---------------------------\n", " 0.01 30000000.0 ... 10000000.0 Feldmeier et al. (2014)\n", " 0.12 800000000.0 ... 200000000.0 Launhardt et al. (2002)\n", " 8.1 89502860861.52429 ... 4858963492.608627 Bovy et al. (2012)\n", " 8.3 110417867208.19055 ... 4387023236.020782 McMillan (2011)\n", " 8.4 102421035406.90356 ... 15468328224.531876 Koposov et al. (2010)\n", " 19.0 208023299175.30438 ... 34833267089.920685 Kuepper et al. (2015)\n", " 50.0 539884832748.48975 ... 268490735257.4718 Wilkinson & Evans (1999)\n", " 50.0 529886965173.18726 ... 38536752776.21277 Sakamoto et al. (2003)\n", " 50.0 399914690706.92847 ... 72696676468.2511 Smith et al. (2007)\n", " 50.0 419910425325.7268 ... 38172735113.64386 Deason et al. (2012)\n", " 60.0 399914690957.5188 ... 64344945146.92987 Xue et al. (2008)\n", " 80.0 689852841359.0248 ... 110361048549.1029 Gnedin et al. (2010)\n", " 100.0 1399701417307.7747 ... 831336271726.1903 Watkins et al. (2010)\n", " 120.0 539884832260.29584 ... 123854764645.32489 Battaglia et al. (2005)\n", " 150.0 750000000000.0 ... 250000000000.0 Deason et al. (2012)\n", " 200.0 679854974257.2006 ... 313652012195.06256 Bhattacherjee et al. (2014)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tbl" ] }, { "cell_type": "markdown", "id": "a049382b", "metadata": {}, "source": [ "Let's now plot the above data and uncertainties:" ] }, { "cell_type": "code", "execution_count": 6, "id": "c8a0247e", "metadata": { "execution": { "iopub.execute_input": "2024-09-30T14:48:28.346184Z", "iopub.status.busy": "2024-09-30T14:48:28.345974Z", "iopub.status.idle": "2024-09-30T14:48:29.036621Z", "shell.execute_reply": "2024-09-30T14:48:29.035864Z" }, "lines_to_end_of_cell_marker": 2 }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwcAAAMACAYAAABxXRv3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAB7CAAAewgFu0HU+AABev0lEQVR4nO3deXzU1b3/8fdkDwmrKIIEiCZhiSAoiaUJQkDjXnFBXIAsmou/an/8LGprahG9dUGEyq31WhKzgOKGikIuV1QWWWwhNgpGMJMKMREQUDAkmUBI5vcHkylD9jAz35nM6/l45NGZ7/fMdz7DMWXenPM9x2S1Wq0CAAAA4PP8jC4AAAAAgGcgHAAAAACQRDgAAAAAYEM4AAAAACCJcAAAAADAhnAAAAAAQBLhAAAAAIAN4QAAAACAJMIBAAAAABvCAQAAAABJhAMAAAAANoQDAAAAAJIIBwAAAABsCAcAAAAAJBEOAAAAANgQDgAAAABIkgKMLgDepba2Vjt37pQknXvuuQoI4D8hAAAAdzt58qQOHTokSRo5cqRCQkKccl2+2aFDdu7cqfj4eKPLAAAAgM22bdsUFxfnlGsxrQgAAACAJEYO0EHnnnuu/fG2bdvUv39/A6sBAAC+6uOPP1Ztba1CQkJ05ZVXGl2O2+3fv98+m+P072dni3CADjn9HoP+/ftr4MCBBlYDAAB8Vd++fWWxWBQaGurz30eceQ8o04oAAAAASCIcAAAAALAhHAAAAACQRDgAAAAAYEM4AAAAACCJcAAAAADAhqVMAQAA4FXMZrPy8/O1b98+DRgwQEOHDlV0dLTRZXUJhAMAAAB4jdzcXGVkZKi+vt5+7L333lNWVpbS0tIMrKxrYFoRAAAAvILZbG4SDCSpvr5eGRkZMpvNBlXWdRAOAAAA4BVycnKaBING9fX1ysnJcXNFXQ/hAAAAAF5hz549Z3UebSMcAAAAwCtERkae1Xm0jXAAAAAAr5Ceni5/f/9mz/n7+ys9Pd3NFXU9hAMP8+qrr2rWrFkaO3asgoODZTKZlJeX12zb77//Xi+88IKSk5M1aNAgBQUF6fzzz9ett96qf/zjH+4tHAAAwMWio6OVlZXVJCD4+/srKyuL5UydgKVMPcxjjz2msrIy9e3bV/3791dZWVmLbf/yl79o/vz5uuiii5ScnKxzzz1XZrNZK1eu1MqVK7V8+XJNmzbNjdUDAAC4VlpamhITE5WZmWnf5+Dpp58mGDgJ4cDDZGdnKzo6WoMHD9azzz6rRx99tMW28fHx2rBhgyZMmOBwfNOmTZo8ebL+z//5P5oyZYqCg4NdXTYAAIDbREdHKyUlRRaLRaGhoQQDJ2JakYe58sorNXjw4Ha1veWWW5oEA0kaP368kpKSdOTIEe3cudPZJQIAAKCLIhzYHDx4UKtXr9bcuXN17bXXqm/fvjKZTDKZTEpNTe3QtcrKyjRnzhwNGzZMYWFh6tOnj+Li4rRgwQLV1NS45gOcITAwUJIUEMDgEAAAANqHb442/fr1c8p1Vq1apenTp6uystJ+rKamRoWFhSosLFR2drYKCgoUFRXllPdrznfffaePP/5Y/fv318iRI132PgAAAOhaGDloxqBBg5ScnNzh1xUVFWnatGmqrKxUeHi4nnrqKW3dulWffPKJMjIyJEklJSW6/vrrdezYMWeXLUmqq6vTjBkzdPz4cc2fP7/F5b4AAACAMzFyYDN37lzFxcUpLi5O/fr10969ezu8kcbs2bNlsVgUEBCgtWvXaty4cfZzkyZNUnR0tB555BGVlJRo4cKFmjdvnlM/Q0NDg1JTU/Xpp58qIyNDM2bMcOr1AQAA0LUxcmDzxBNP6IYbbuj09KJt27Zp06ZNkqR77rnHIRg0mjNnjoYPHy5JWrx4serq6jpf8BkaGhqUnp6u5cuXa/r06Xr55Zeddm0AAAD4BsKBk6xcudL+OC0trdk2fn5+mjlzpiTp6NGjWr9+vVPeu6GhQWlpacrPz9edd96pvLw8+fnRtQAAAOgYvkE6yebNmyVJYWFhuuyyy1psd/rSo1u2bDnr920MBkuXLtW0adO0bNky7jMAAABApxAOnGTXrl2SpKioqFaXDx02bFiT13RW41SipUuXaurUqXr11VcJBgAAAOg0bkh2gtraWh0+fFiSNHDgwFbb9u7dW2FhYaqurlZ5eXmT89nZ2fZRiMYNzLKzs7VhwwZJUmJiou69915J0pNPPqn8/HyFh4crJiZGf/rTn5pcb8qUKRo9enS7P0tFRUWr5/fv39/uawEAAMC7EA6c4PRlScPDw9ts3xgOqqqqmpzbvHmz8vPzHY5t2bLFYQpSYzjYu3evJKmqqkpPPfVUs+81ZMiQDoWDiIiIdrcFAABA10I4cILa2lr746CgoDbbBwcHS5IsFkuTc3l5ecrLy2vX+3akLQAAAJpXUlKiuro6BQYGKiYmxmWv8QaEAycICQmxPz5x4kSb7Y8fPy5JCg0NdVlNndXcVKfT7d+/X/Hx8W6qBgAAwPVKSkpksVgUGhraoXDQ0dd4A8KBE3Tv3t3+uLmpQmeqrq6W1L4pSO7W1j0TAAAA6LpYrcgJQkJCdM4550hq+4beI0eO2MMB8/sBAADgSRg5cJIRI0Zo06ZNKi0t1cmTJ1tcznT37t32x427JXuy2NhYh+fO3NUZAAAAnoWRAydJTEyUdGrK0Oeff95iu40bN9ofJyQkuLwuAAAAoL0IB04yZcoU++Pc3Nxm2zQ0NGjp0qWSpF69eikpKckdpZ2V4uJih59169YZXRIAAABchHDgJPHx8Ro/frwk6ZVXXtFnn33WpM3ChQvtuyLPnj1bgYGBbq0RAAAAaA33HNhs3rxZpaWl9ueNOx5LUmlpaZP9BFJTU5tcY/HixUpISJDFYlFycrIyMzOVlJQki8WiN954Q0uWLJEkxcTEaM6cOS75HAAAAEBnEQ5ssrOzm+xM3OjMHYql5sPBmDFj9Oabb2r69OmqrKxUZmZmkzYxMTEqKChwWP4UAAAA8ASEAye78cYbtWPHDi1evFgFBQWqqKhQUFCQoqKiNHXqVD3wwAPq1q2b0WW2G6sVAQCArsxsNis/P1/79u3TgAEDNHToUEVHRxtdlmEIBzZ5eXlNpg511uDBg7Vo0SItWrTIKdcDAADwJCUlJaqrq1NgYKBTdwd21XVbkpubq4yMDNXX19uPvffee8rKylJaWprL398TEQ7QquLiYofnFRUVbN4GAICPKykpkcViUWhoqNPDgSuu2xyz2dwkGEhSfX29MjIylJiY6JMjCKxWBAAAAJ+Tk5PTJBg0qq+vV05Ojpsr8gyEAwAAAPicPXv2nNX5roppRQAAADCcu28MjoyM7PT5rnwTMyMHAAAAMFRubq6GDx+uFStWaOvWrVqxYoWGDx+u3Nxcl71nenq6/P39mz3n7++v9PR0j6nVnQgHaFVsbKzDz6RJk4wuCQAAdCFt3RhsNptd8r7R0dHKyspqEhD8/f2VlZXV7EiAUbW6E+EAAAAAhjHyxuC0tDTt2rVLt912m375y1/qtttu065du1pcxtQXbmLmngO0iqVMAQCAKxl9Y3B0dLRSUlLsS6i2du+A0bW6AyMHAAAAMMzZ3Bjsbt5Ua2cRDgAAAGCYzt4Y3JaSkhIVFxerpKTkbMpz4KpaPQnhAAAAAIbpzI3B7VFSUqKvv/7aqeHAVbV6Eu45AAAAgKHS0tKUmJiozMxM+94BTz/9tEd+2famWjuDcAAAAADDdeTGYKN5U60dRThAq2JjYx2e19XVGVQJAAAAXI17DgAAAABIYuQAbWCfAwAAcDqz2az8/Hz7fPuhQ4d63LQab6jRUxEOAAAA0C65ubm699571dDQYD/27rvvKjs7u8Vdhd0tNzdXGRkZDjsZv/fee8rKyvKYGj0Z04oAAADQJrPZrHvuucchGEhSQ0OD7rnnHpnNZoMq+zez2dwkGEhSfX29MjIyPKJGT0c4AAAAQJuef/55Wa3WZs9ZrVY9//zzbq6oqZycnCbBoFF9fb1ycnLcXJH3IRwAAACgTZs2bTqr8+6wZ8+eszoPwgEAAAC6iMjIyLM6D8IBAAAA2mH8+PFndd4d0tPT5e/v3+w5f39/paenu7ki70M4QKtiY2MdfiZNmmR0SQAAwAAPPfSQ/Pya/+ro5+enhx56yM0VNRUdHa2srKwmAcHf319ZWVksZ9oOLGUKAACANkVHRys7O7vJakBGffGOiYlRXV2dAgMDHY6npaUpMTFRmZmZ9n0Onn766Vbra+lavohwgFaxCRoAAGjU+MU7JydHe/bsUWRkpNLT0w35F/mYmJgWz0VHRyslJUUWi0WhoaFt1tfatXwN4QAAAADtFh0drWeeecboMuAihAMAAAB0SElJiX0aDv/q3rUQDgAAANAhJSUl9ik7zYUDwoP3IhwAAACg3cxms3Jzc3XgwAGdf/75Gjp0aJM5/W2FB3guwgEAAADaJTc3t8lqRe+//76ysrKUlpZmYGVwFvY5AAAAQJvMZnOTYCBJ9fX1ysjIkNlsNqgyOBPhAAAAAG3KyclpEgwa1dfXKycnR9KpEJGfn68XXnhB+fn5hAYvw7QiAAAAtGnPnj1tnm9u2tF7773HtCMvwsgBAAAA2tSzZ8822zDtyPsxcoBWxcbGOjyvq6szqBIAAODJvvnmmzanHbF5mudj5AAAAABt+vnnn1s9/+OPP7Z6vq1pSfAMjBygVcXFxQ7PKyoqFBERYVA1AADAKJGRka2ev+CCC1ReXt7p10tSTEyMffM0GIORAwAAALQpPT1d/v7+zZ7z9/fXk08+2er59PT0Nt8jJiZGsbGxbJxmIMIBAAAA2hQdHa2srKwmAcDf319ZWVm66qqrWj1/5i7K8ExMKwIAAEC7pKWlKTExUZmZmdq3b58GDBigp59+2v7Fv63z8HyEAwAAALRbdHS0UlJSZLFYFBoa2uSLf1vn4dmYVgQAAABAEuEAAAAAgA3TigAAANDlsCxq5xAOAAAA4FSe8MWc5VA7h3AAAAAAp+KLuffingMAAAAAkggHAAAAAGyYVoRWxcbGOjyvq6szqBIAAAC4GiMHAAAAACQxcoA2FBcXOzyvqKhQRESEQdUAAADAlQgHAAAA6BBPWKoUrkE4AAAAQIewVGnXxT0HAAAAACQRDgAAAADYMK0IAAAA6KCuet8F4QAAAADooK563wXTigAAAABIIhwAAAAAsCEcAAAAAJDEPQcAAACGKSkpsd/U2lXnsMO7EA4AAAAMUlJSIovFotDQUMIBPALTigAAAABIIhwAAAAAsGFaEQAAgAHMZrPy8/O1b98+DRgwQEOHDlV0dLTRZcHHEQ4AAADcLDc3VxkZGaqvr7cfe++995SVlaW0tDQDK4OvY1oRAACAG5nN5ibBQJLq6+uVkZEhs9lsUGUA4QAAAMCtcnJymgSDRvX19crJyXFzRcC/EQ4AAADcaM+ePWd1HnAl7jlAq2JjYx2e19XVGVQJAABdQ2Rk5FmdB1yJkQMAAAA3Sk9Pl7+/f7Pn/P39lZ6e7uaKgH9j5ACtKi4udnheUVGhiIgIg6oBAMD7RUdHKysrq8lNyf7+/srKymI5UxiKcAAAAOBmaWlpSkxMVGZmpn2fg6effppgAMMRDgAAAAwQHR2tlJQUWSwWhYaGEgzgEbjnAAAAAIAkwgEAAAAAG8IBAAAAAEmEAwAAAAA2hAMAAAAAkggHAAAAAGxYyhQAAMAgMTExqqurU2BgoNGlAJIIBwAAAIaJiYkxugTAAdOKAAAAAEgiHAAAAACwIRwAAAAAkEQ4AAAAAGBDOAAAAAAgiXAAAAAAwIZwAAAAAEAS4QAAAACADeEAAAAAgCTCAQAAAAAbwoGHefXVVzVr1iyNHTtWwcHBMplMysvLc1p7AAAAoCUBRhcAR4899pjKysrUt29f9e/fX2VlZU5tDwAAALSEkQMPk52drb179+rQoUO67777nN4eAAAAaAkjBx7myiuvdGl7AAAAoCWMHNgcPHhQq1ev1ty5c3Xttdeqb9++MplMMplMSk1N7dC1ysrKNGfOHA0bNkxhYWHq06eP4uLitGDBAtXU1LjmAwAAAAdms1mPPvqo7rjjDj366KMym81GlwR4PEYObPr16+eU66xatUrTp09XZWWl/VhNTY0KCwtVWFio7OxsFRQUKCoqyinvBwAAmsrNzVVGRobq6+vtxxYsWKCsrCylpaUZWBng2Rg5aMagQYOUnJzc4dcVFRVp2rRpqqysVHh4uJ566ilt3bpVn3zyiTIyMiRJJSUluv7663Xs2DFnlw0AAHRqxODMYCBJ9fX1ysjIYAQBaAXhwGbu3LlatWqVDhw4oLKyMv3tb3/r8DVmz54ti8WigIAArV27VpmZmRo3bpwmTZqkJUuW6LnnnpN0KiAsXLjQ2R8BAABIysnJaRIMGtXX1ysnJ8fNFQHeg3Bg88QTT+iGG27o9PSibdu2adOmTZKke+65R+PGjWvSZs6cORo+fLgkafHixaqrq+t8wQAAoFl79uw5q/OALyMcOMnKlSvtj1uay+jn56eZM2dKko4ePar169e7ozQAAHxKZGRkq+eZ2gu0jHDgJJs3b5YkhYWF6bLLLmux3YQJE+yPt2zZ4vK6AADwNenp6fLza/krzocffsh9B0ALWK3ISXbt2iVJioqKUkBAy3+sw4YNa/IaT1JRUdHq+f3797upEgAAOic6OlrXXHON/ud//qfZ8433HTzzzDNurgzwfIQDJ6itrdXhw4clSQMHDmy1be/evRUWFqbq6mqVl5c3OZ+dnW0fhdi5c6f92IYNGyRJiYmJuvfeezvdvi0RERHtbgsAgKfq3r17q+e57wBoHuHACU6fuxgeHt5m+8ZwUFVV1eTc5s2blZ+f73Bsy5YtDlOQTv+y39H2AAD4grbuO2jrPOCrCAdOUFtba38cFBTUZvvg4GBJksViaXIuLy9PeXl57X7vjrZvS3OjGafbv3+/4uPjnfZ+AAC4Qnp6uhYsWNDskqb+/v5KT083oCrA8xEOnCAkJMT++MSJE222P378uCQpNDTUZTV1VlvTogAA8BZXX3211qxZI6vVaj/m7++vrKwsRUdHG1gZ4LkIB05w+rzG5qYKnam6ulpS+6YgAQCAjsnNzW12h+SxY8dq+fLlBAOgFSxl6gQhISE655xzJLW92s+RI0fs4YCbfwEAcC6z2dxsMJCkf/7znwZUBHgXRg6cZMSIEdq0aZNKS0t18uTJFpcz3b17t/1x427Jniw2NtbhObs6AwA8WU5OTrPBQJIaGhpYwhRoAyMHTpKYmCjp1JShzz//vMV2GzdutD9OSEhweV0AAPiStpYoZQlToHWEAyeZMmWK/XFubm6zbRoaGrR06VJJUq9evZSUlOSO0s5KcXGxw8+6deuMLgkAgBaxhClwdggHThIfH6/x48dLkl555RV99tlnTdosXLjQvivy7NmzFRgY6NYaAQDo6tLT0+Xv79/sOZYwBdrGPQc2mzdvVmlpqf15447HklRaWtpkL4HU1NQm11i8eLESEhJksViUnJyszMxMJSUlyWKx6I033tCSJUskSTExMZozZ45LPgcAAL4sOjpaWVlZTW5KZglToH1M1tMX//VhqampTXYabk1Lf2yrVq3S9OnTVVlZ2ez5mJgYFRQUKCoqqlN1Gq2iosK+ylJ5eTn7IgAAPJLZbFZOTo727NmjyMhIpaenEwzQpbjqOxkjB0524403aseOHVq8eLEKCgpUUVGhoKAgRUVFaerUqXrggQfUrVs3o8tsN1YrAgB4o+joaFYlAjqBkQO0qrlwYDabJTFyAAAAYBRGDmCI4uJih+en/4cIAACAroXVigAAAABIIhwAAAAAsCEcAAAAAJBEOAAAAABgww3JaBVLmQIAAPgORg4AAAAASGLkAG1gKVMAAADfwcgBAAAAAEmEAwAAAAA2hAMAAAAAkggHAAAAAGwIBwAAAAAksVoR2sA+BwAAAL6DkQMAAAAAkhg5QBvY5wAAAMB3MHIAAAAAQBLhAAAAAIAN4QAAAACAJMIBAAAAABvCAQAAAABJhAMAAAAANixlilaxCRoAAIDvYOQAAAAAgCRGDtAGNkEDAADwHYwcAAAAAJBEOAAAAABgQzgAAAAAIIlwAAAAAMCGcAAAAABAEuEAAAAAgA3hAAAAAIAkwgEAAAAAGzZBQ6tiY2MdntfV1RlUCQAAAFyNkQMAAAAAkhg5QBuKi4sdnldUVCgiIsKgagAAAOBKjBwAAAAAkEQ4AAAAAGBDOAAAAAAgiXAAAAAAwIZwAAAAAEAS4QAAAACADeEAAAAAgCTCAQAAAAAbwgEAAAAASYQDAAAAADaEAwAAAACSCAcAAAAAbAKMLgCeLTY21uF5XV2dQZUAAADA1Rg5AAAAACCJkQO0obi42OF5RUWFIiIiDKoGAAAArsTIAQAAAABJhAMAAAAANoQDAAAAAJIIBwAAAABsCAcAAAAAJBEOAAAAANgQDgAAAABIIhwAAAAAsCEcAAAAAJBEOAAAAABgQzgAAAAAIIlwAAAAAMCGcAAAAABAEuEAAAAAgA3hAAAAAIAkwgEAAAAAG8IBAAAAAEmEAwAAAAA2AUYXAM8WGxvr8Lyurs6gSgAAAOBqjBwAAAAAkMTIAdpQXFzs8LyiokIREREGVQMAAABXYuQAAAAAgCTCAQAAAAAbwgEAAAAASYQDAAAAADYuuSH5ySefdMVlm5g7d65b3gcAAADwBSar1Wp19kX9/PxkMpmcfdkm6uvrXf4ecHT6akXl5eUaOHCgwRUBAAD4Hld9J3PpUqYuyB127ggfAAAAgC9x6T0HX331lRoaGpz6s2PHDleWDAAAAPgsr7shmREDAAAAwDW8LhwAAAAAcA2X3HOwfv16SVJkZKTTrx0ZGWm/PgAAAADncUk4mDBhgisuK0nq1q2bS68PAAAA+CqmFQEAAACQRDgAAAAAYEM4AAAAACCJcAAAAADAhnAAAAAAQBLhAAAAAIAN4QAAAACAJBftc+AMFotFa9eu1a5du9TQ0KCYmBhdc801Cg8PN7o0AAAAoEvyyJGDZcuWKTIyUjfffLP+9Kc/6dlnn9Xtt9+uQYMG6eWXXza6PJd69dVXNWvWLI0dO1bBwcEymUzKy8tr9TXbt2/Xddddp169eiksLEy/+MUv9NZbb7mnYAAAAHQZHhcOnnnmGaWkpGjChAnauXOnqqqqVFlZqd27d2vKlCm6//779eijjxpdpss89thjWrJkicrKytS/f/82269fv14JCQnavHmzbr/9dt133306cOCApk2bpoULF7qhYgAAAHQVHhUONm3apMcee0x/+MMf9Oabbyo2NtZ+LiYmRjk5OZo/f76ee+45rV692sBKXSc7O1t79+7VoUOHdN9997Xa9uTJk8rIyJCfn58+/fRTLVmyRAsXLtSXX36pmJgYZWZmqqyszE2VAwAAwNt5VDh48sknFR8fr//8z/9ssc1DDz2kyZMna+7cuW6szH2uvPJKDR48uF1t161bp3/961+66667NHr0aPvxnj17KjMzUydOnFB+fr6LKgUAAEBX4zHhoKqqShs2bNB//Md/2I+dPHlS3377rb799lvV1dXZj8+aNUtffvml9u3b57T3P3jwoFavXq25c+fq2muvVd++fWUymWQymZSamtqha5WVlWnOnDkaNmyYwsLC1KdPH8XFxWnBggWqqalxWs0bNmyQJCUnJzc5d/XVV0uSNm7c6LT3AwAAQNfmMasVlZeXq76+XiNGjLAfKysrU3R0tPz8/LR161bFx8dLki6++GJZrVaVlZVpwIABTnn/fv36OeU6q1at0vTp01VZWWk/VlNTo8LCQhUWFio7O1sFBQWKioo66/cym82SpOjo6Cbnzj//fIWHh9vbAAAAAG3xmJGDoKAgSdKJEyfsxwIDAzVo0CBFREQoODjYfvz48eMymUz21zjboEGDmv3X+LYUFRVp2rRpqqysVHh4uJ566ilt3bpVn3zyiTIyMiRJJSUluv7663Xs2LGzrvPnn3+WdGoaUXN69OhhbwMAAAC0xWNGDgYNGqTu3bvrs88+0/jx4+3H9u7d26Tt1q1bFRAQoJiYGKe9/9y5cxUXF6e4uDj169dPe/fuVWRkZIeuMXv2bFksFgUEBGjt2rUaN26c/dykSZMUHR2tRx55RCUlJVq4cKHmzZvntPoBAACAs+UxIweBgYG65ZZb9F//9V+qqqpqsV1tba0WLVqka665Rt27d3fa+z/xxBO64YYbOj29aNu2bdq0aZMk6Z577nEIBo3mzJmj4cOHS5IWL17scB9FZzSOGLQ0OlBZWdniqAIAAABwJo8JB9KpL+g1NTW65ZZbmg0Ix48f17Rp07R//37Nnz/fgApbtnLlSvvjtLS0Ztv4+flp5syZkqSjR49q/fr1Z/WejfcaNHdfwYEDB1RVVdXs/QgAAABAczwqHAwaNEjvvPOOtm3bposvvlh/+ctf9I9//EOFhYV6+eWXNWrUKH388cd64403NGzYMKPLdbB582ZJUlhYmC677LIW202YMMH+eMuWLWf1no3XWrt2bZNzH374YZP3AwAAAFrjMfccNEpKSlJhYaEeeeQRPfTQQzp58qSsVqv8/f111VVX6Z133tHFF19sdJlN7Nq1S5IUFRWlgICW/1hPDzWNr+msyZMn68ILL9Ty5cv1f//v/7XvdfDzzz/r6aefVlBQkH2kor0qKipaPb9///7OlgsAAAAP53HhQDr1Bfvdd9/VsWPHZDab1dDQoKioKPXq1cvo0ppVW1urw4cPS5IGDhzYatvevXsrLCxM1dXVKi8vb3I+OzvbPgqxc+dO+7HGPQ0SExN17733SpICAgKUnZ2tq6++WldccYXuuOMOde/eXe+8847Kysr0/PPPa8iQIR36LBERER1qDwAAgK7DI8NBo+7du+vSSy81uow2nb4saXh4eJvtG8NBc/dVbN68ucmuxlu2bHGYgtQYDqRTIy2bN2/W448/rjfffFN1dXUaOXKk5s+fr2nTpnXm4wAAAMBHeXQ48Ba1tbX2x+3Ze6FxzwaLxdLkXF5envLy8jr0/vHx8VqzZk2HXtOS5kYzTrd//377ZnQAAADoWggHThASEmJ/fPombi05fvy4JCk0NNRlNXVWW9OiAAAA0HV51GpF3ur0/RZa26OhUXV1taT2TUECAAAA3IVw4AQhISE655xzJLW92s+RI0fs4YCbfwEAAOBJmFbkJCNGjNCmTZtUWlqqkydPtric6e7du+2PG3dL9mSxsbEOz892V2cAAAB4LkYOnCQxMVHSqSlDn3/+eYvtNm7caH+ckJDg8roAAACA9nLZyMGnn37q9GteccUVTr+ms0yZMkXPPPOMJCk3N1eXX355kzYNDQ1aunSpJKlXr15KSkpya42dUVxc7PC8oqKC6VAAAABdlMvCwcSJE2UymZx2PZPJpJMnTzrtes4WHx+v8ePHa9OmTXrllVeUkpKicePGObRZuHChfVfk2bNnKzAw0IhSAQAAgGa5/J4Dq9Xq6rdwis2bN6u0tNT+vHHHY0kqLS1tsvdAampqk2ssXrxYCQkJslgsSk5OVmZmppKSkmSxWPTGG29oyZIlkqSYmBjNmTPHJZ8DAAAA6CyT1UXf3v38/GQymRQSEqKbbrpJV111lfz8zu4Wh5SUFCdV11RqamqTnYlb09If26pVqzR9+nRVVlY2ez4mJkYFBQWKiorqVJ1GO31aUXl5OfsiAAAAGMBV38lcNnLQvXt3HTt2TBaLRW+++aY2btyou+66SzNmzNCoUaNc9baGu/HGG7Vjxw4tXrxYBQUFqqioUFBQkKKiojR16lQ98MAD6tatm9FlthurFQEAAPgOl40c1NbW6v3339eyZcu0du1anTx50n4PwsiRIzVz5kzdeeed6t+/vyveHk7SXDgwm82SGDkAAAAwiqtGDlwWDk536NAhLV++XMuWLdM///nPU29sMsnPz0+TJ0/WzJkzdfPNNys0NNTVpeAsMa0IAADAeK76TuaWfQ7OPfdczZ49W4WFhSouLtbvfvc7DRw4UPX19Vq7dq1mzJihfv36KTU1VZ988ok7SgIAAABwBrdvgjZ8+HA988wzKisr07p165Samqrw8HBVVVVp6dKlSk5OVkREhP7whz+4uzQAAADApxm6Q/LEiROVk5OjH374QcuXL9e1114rf39/ff/99/rzn/9sZGkAAACAz3H5Pgft0Xj/gclkcurGaTh7rFYEAADgOwwNBxs3btSyZcv0zjvv2PcFsFqt6t+/v2bMmGFkaQAAAIDPcXs42LVrl5YtW6bly5ervLxc0qlA0K1bN918882aOXOmJk+efNYbpsE5iouLHZ6ffmc8AAAAuha3hIODBw/q9ddf17Jly1RUVCTpVCDw8/NTUlKSZs6cqVtuuUVhYWHuKAcAAABAM1wWDmpra7Vy5UotW7ZMH330kerr69W4pUJsbKxmzpypu+++WwMGDHBVCQAAAAA6wGXh4LzzzlN1dbWkU6ME559/vu68807NmDFDo0ePdtXbAgAAAOgkl4WDqqoqmUwmhYSE6Fe/+pWSk5Pl7++vHTt2aMeOHZ265syZM51cJQAAAIBGJmvjXB8na1ya1FlMJpNOnjzptOuhc1y1VTcAAADaz1XfyVx6Q7KLcgfciH0OAAAAfIfLwsH69etddWkAAAAALuCycDBhwgRXXRpuxD4HAAAAvoOdxgAAAABIIhwAAAAAsCEcAAAAAJDkonsOvvvuO0nSBRdcIH9/f6deu76+Xt9//70kadCgQU69NgAAAODLXDJyMGTIEF144YX65ptvnH7t3bt3268PAAAAwHlcNq3I1XscsIcCAAAA4FwuvefAmTskAwAAAHAtl+6QnJycrMDAQKdekx163YsdkgEAAHyHy8KB1Wq13zgMAAAAwPO5JBykpKS44rIwADskAwAA+A6XhIPc3FxXXBYAAACAC7EJGgAAAABJhAMAAAAANoQDAAAAAJIIBwAAAABsCAcAAAAAJBEOAAAAANgQDgAAAABIcuEOyegaYmNjHZ7X1dUZVAkAAABcjZEDAAAAAJIYOUAbiouLHZ5XVFQoIiLCoGoAAADgSowcAAAAAJBk4MjBoUOH9O233+rAgQOqrq5WYGCgevXqpUGDBikqKkr+/v5GlQYAAAD4JLeFg+rqar3//vtas2aNNm7cqO+//77FtsHBwRozZoySk5N18803a9SoUe4qEwAAAPBZJqvVanXlGxQVFekvf/mL3n77bdXU1EiS2vuWJpNJ0qkVc+6//37NmDFD3bp1c1mtaNvp9xyUl5dr4MCBBlcEAADge1z1ncxl4aCoqEh//OMftWbNGkn/DgTnn3++4uPjddlll+m8885Tnz591Lt3b1ksFv300086cuSISkpKtH37du3YscO+dKbJZNI555yjRx55RL/5zW8UHBzsirLRBsIBAACA8Vz1ncwl04rS0tK0bNkyNTQ0SJIuvfRS3X333br11ls1aNCgdl/nxIkT+vTTT/Xaa6/pvffe0+HDh/W73/1OL730kpYuXarExERXlA8AAAD4JJesVpSfn6+AgABlZGRo9+7dKiws1IMPPtihYCBJQUFBuvLKK5Wbm6sffvhBS5cu1dChQ7V3716tW7fOFaUDAAAAPsslIwe//vWv9bvf/c6p6+EHBwdr+vTpuvvuu/X222+rvr7eadcGAAAA4KJw8OKLL7rispJO3Xtw++23u+z6AAAAgK9y6yZokyZN0oMPPqi///3v7nxbAAAAAO3g1k3QNm7cqI0bNyopKcmdbwsAAACgHdw6cnDuuedKklPvRQAAAADgHG4dORg1apQ++eQT7du3T2PGjHHnW6OTYmNjHZ437jsBAACArsetIwcpKSmyWq1688033fm2AAAAANrBZTskt+SKK67Q1q1btWrVKl177bXufGs4ATskAwAAGM9V38ncOnLw5ptv6sknn9Tw4cN12223KT8/351vDwAAAKAVbg0Hd955pyZPnqyvv/5aFotF6enp+sUvfqH8/Hzt37/fnaUAAAAAOINbb0iWpNNnMVmtVm3fvl3bt2+XJF1wwQUaM2aMRo8erTFjxmjMmDEaPHiwu0sEAAAAfJJbw8GuXbv05Zdf2n+++OIL7du3z36+oqJC33//vVavXm0/1rNnT3tYWLhwoTvLBQAAAHyK229IPtOPP/7oEBa+/PJL7dq1q8mSmSaTSfX19QZViUbckAwAAGA8V30nc/u0ojOdc845mjRpkiZNmmQ/VldXp127dtnDQuMPAAAAANcxPBw0JzAwUKNGjdKoUaOMLgUAAADwGW5drQgAAACA5yIcAAAAAJBEOAAAAABgQzgAAAAAIIlwAAAAAMCGcAAAAABAEuEAAAAAgI3b9jkoKirSzp075efnp+nTp7frNW+//bYsFouio6M1btw4F1cIAAAA+Da3hYOffvpJqampMplMGjhwoCZOnNhq+927d2vatGkymUxaunQp4QAAAABwMbdNK5o0aZIiIiIkScuWLWuzfWOb7t2769Zbb3VpbQAAAADcGA5MJpOmT58uq9WqFStWqLa2ttX2r732mkwmk2677TaFhIS4qUqcKTY21uFn0qRJRpcEAAAAF3HrDcmpqamSpKqqKr333nstttu4caO+++47h9cAAAAAcC23hoPTbyzOz89vsd3SpUslSRdeeKESExPdUhuaV1xc7PCzbt06o0sCAACAi7h9KdOUlBRZrVZ98sknOnDgQJPztbW1WrFihUwmk2bOnOnu8gAAAACf5fZwcMcddygkJEQNDQ1avnx5k/Pvv/++jh07RjgAAAAA3Mzt4aBHjx666aabZLVa7dOHTte4StEVV1yhwYMHu7s8AAAAwGcZskNy403GO3fu1Jdffmk/fujQIa1du1Ymk0kpKSlGlAYAAAD4LEPCwVVXXaUBAwZIctzzYPny5Tp58qTCwsI0depUI0oDAAAAfJYh4cDPz0933323rFarli9froaGBkmngoLJZNKtt96qbt26GVEaAAAA4LMMCQeSlJaWJkn64Ycf9NFHH2nXrl365z//KUlMKQIAAAAMEGDUGw8bNkxxcXEqLCxUfn6+hgwZIkkaPHiwJk6caFRZAAAAgM8yLBxIp0YItm/frg8++EC9evWSyWTSjBkzjCwJAAAA8FmGTSuSpLvuukvBwcGyWCzat2+fpH+vZAQAAADAvQwNB7169dKNN94oq9Uqk8mkhIQERUZGGlkSAAAA4LMMDQfSv0cKrFYrowYAAACAgQy950CSrrvuOvtSpgAAAACMY/jIAQAAAADPQDgAAAAAIIlwAAAAAMCGcAAAAABAEuEAAAAAgA3hAAAAAIAkwgEAAAAAG8IBAAAAAEmEA6/X0NCgF198UZdeeqm6deumHj166IorrtAHH3xgdGkAAADwMoQDL2a1WnX77bfrN7/5jSorK3XPPffojjvu0DfffKObbrpJL774otElAgAAwIu4JBy8++67rris3b59+/T3v//dpe/hDd555x298847SkhI0M6dO/WXv/xFS5YsUXFxsQYPHqyHHnpIe/fuNbpMAAAAeAmXhIPbbrtNo0eP1ooVK5x63fLycv3617/WRRddpLVr1zr12t7o/ffflyRlZmYqNDTUfrxv37568MEHdfz4ceXm5hpVHgAAALyMS8LBRRddpB07dmjatGmKjIzUH/7wBxUXF3fqWtXV1Xr11Vd13XXX6aKLLtLLL7+s+vp6XXTRRU6uumMOHjyo1atXa+7cubr22mvVt29fmUwmmUwmpaamduhaZWVlmjNnjoYNG6awsDD16dNHcXFxWrBggWpqalp83YEDByRJkZGRTc41Hlu3bl2HagEAAIDvCnDFRb/++mu98MILeu6551RWVqZnn31Wzz77rKKjo/WLX/xCcXFxGjNmjM477zz17t1bvXv3lsVi0U8//aQjR46opKRE27dv17Zt27Rt2zbV1tbKarVKkm655RY9/fTTiomJcUXp7davXz+nXGfVqlWaPn26Kisr7cdqampUWFiowsJCZWdnq6CgQFFRUU1e27dvX0nSnj17NHz4cIdze/bskSSVlJQ4pU4AAAB0fSZr47duF6iqqtJLL72kv/71ryovLz/1hiZTu1/fWFpwcLBuueUWzZ49W/Hx8S6ptaNO/xyDBg3SsGHD7FOdUlJSlJeX1+Y1ioqKlJCQIIvFovDwcD366KNKSkqSxWLRG2+8oaysLElSTEyMCgsL1b17d4fXL126VCkpKRo/frzWrl2rkJAQSdKPP/6osWPHau/evQoKCtLx48ed9KmliooKRURESDo1zWvgwIFOuzYAAADax1XfyVwyctAoPDxcjzzyiB566CF99NFHeuutt7R+/fp23SQbEhKiyy+/XDfddJNmzpypPn36uLLUDps7d67i4uIUFxenfv36ae/evc1O72nN7NmzZbFYFBAQoLVr12rcuHH2c5MmTVJ0dLQeeeQRlZSUaOHChZo3b57D6++66y7l5eVp/fr1GjlypK655hrV1dVp5cqV9pENPz8WpAIAAED7uHTkoCXff/+9tm7dqoqKCh06dEg//fSTQkJCdO655+rcc8/VyJEjNXbsWAUGBrq7tE47PRy0Z+Rg27ZtuvzyyyVJs2bN0ssvv9ykTUNDgy6++GLt2rVLvXr10sGDB5v8mRw/flzPPvusli9frr1796pnz566+eab9dBDDykmJkaDBg1SWVmZcz6kGDkAAADwBF45ctCSCy64QFOnTjXirT3GypUr7Y/T0tKabePn56eZM2fq0Ucf1dGjR7V+/XolJyc7tAkODtbjjz+uxx9/3OH4hg0bJEljx451at0AAADoulw252T//v2uunSXsHnzZklSWFiYLrvsshbbTZgwwf54y5Yt7b7+a6+9Jkm64447OlkhAAAAfI3LRg4uuOACnXfeebrkkks0evRojR49WmPGjNHQoUM7dFNyV7Vr1y5JUlRUlAICWu6GYcOGNXnN6SorK9WjRw+HYytWrFBOTo7i4uJ0yy23dKiuioqKVs8T+gAAALoul04rOnTokD7++GN9/PHH9mOhoaG6+OKLHQLDqFGjHDbx6upqa2t1+PBhSWpzfljv3r0VFham6upq+4pPp7v88ssVERGh4cOHKyQkRNu2bdOGDRt04YUX6u2335a/v3+HamucuwYAAADf4/J7Ds6837mmpkbbt2/X9u3b7cf8/PwUFRXlEBhGjx6t8847z9XlGeLYsWP2x+Hh4W22bwwHVVVVTc5NmzZN7777rv7+97+rrq5OkZGReuyxx/Twww83GVEA0PWUlJSorq5OgYGBhu//AgDwfi4PB8HBwfrVr36l66+/Xnv27NEXX3yhL774wmEFnfr6en3zzTcqKSnRW2+9ZT/er18/e2AYPXq0br/9dleX6xa1tbX2x0FBQW22Dw4OliRZLJYm5+bNm9dkidOz0dzoxOn279/vMXtNADgVDiwWi0JDQwkHAICz5rJwsGrVKv3ud7/T119/rRUrVmjdunX64x//qBUrVsjf318///yzPSg0/nz99deqq6uzX+PAgQP68MMP9eGHH8pkMnWZcNC4WZkknThxos32jZuYuWPqFUuTAgAA+C6XrVZ0/fXXa8eOHXr55ZfVr18//fjjj3rwwQc1YsQIvffee+rZs6cmTJig2bNnKzc3V0VFRaqqqlJRUZFyc3M1e/ZsTZgwQT179pTVam0yPcmbnb7TcXNThc5UXV0tqX1TkAAAAIDOcum0Ij8/P/3Hf/yH7r77bj333HNatGiRzGazbrvtNv3yl7/U888/b98ITJICAwN1ySWX6JJLLlFKSor9eFlZmb744gtXlupWISEhOuecc/Tjjz+2uTrQkSNH7OHAiJuFY2NjHZ6fPrIDAACArsVlIwenCwsL0xNPPKGSkhKlp6fLz89PW7Zs0S9/+UtNmzZN3377bauvHzx4sG666SZ3lOo2I0aMkCSVlpbq5MmTLbbbvXu3/fHw4cNdXhcAAAB8l1vCQaP+/fsrOztbRUVFSk5OltVq1YoVKzRixAjNmTNHR44ccWc5hkpMTJR0asrQ559/3mK7jRs32h8nJCS4vK4zFRcXO/ysW7fO7TUAAADAPdwaDhpdfPHF+t///V99+OGHGjlypE6cOKEXXnhBUVFRWrRokU9MXZkyZYr9cW5ubrNtGhoatHTpUklSr169lJSU5I7SAHgJs9ms/Px8vfDCC8rPz5fZbDa6JACAlzMkHDS66qqrVFRUpJycHPXv319HjhzRww8/rKFDhzosadoVxcfHa/z48ZKkV155RZ999lmTNgsXLrTvijx79mwFBga6tUYAnis3N1fDhw/XihUrtHXrVq1YsULDhw9v8R8bAABoD5PVQ5YBqq2t1RNPPKEFCxaooaFBoaGh9htxPdHmzZtVWlpqf3748GE9/PDDkk5N/7n33nsd2qempja5RlFRkRISEmSxWBQeHq7MzEwlJSXJYrHojTfe0JIlSyRJMTExKiwsdFjlyCgVFRX2G6PLy8tZ+hQwgNls1vDhw1VfX9/knL+/v3bt2qXo6GgDKgMAuIurvpO5fBO05tTU1Ojrr79WcXGx/X+Li4tVXl5uX7LUQzJLi7Kzs5Wfn9/suS1btmjLli0Ox5oLB2PGjNGbb76p6dOnq7KyUpmZmU3axMTEqKCgwLBgwGpFgOfJyclpNhhIpzaVzMnJ0TPPPOPmqgAAXYFLw0F7QoDkGAQCAwM1dOhQXXrppa4szWPceOON2rFjhxYvXqyCggJVVFQoKChIUVFRmjp1qh544AF169bN6DIBeJA9e/ac1XkAAFrisnAQGRnZagiQTq1eNGrUKI0aNUojR47UqFGjNHz4cK+YW5+Xl6e8vDynXGvw4MFatGiRFi1a5JTrOVNxcbHD89OHsAAYIzIy8qzOAwDQEpeFg7KyMvvjbt26KTY21h4AGn/69OnjqrcHgC4rPT1dCxYsaPGeg/T0dAOqAgB0BS6dVmQymdStWzclJydrzJgxGj16tEaPHs1NrABwFqKjo5WVlaWMjAyHgODv76+srCxuRgYAdJrLVivy8/v3Kqkmk8nhXO/evXXJJZdo9OjR9v8dMWKEAgIMuT8aHcBqRYDnMJvNyszM1L59+zRgwAA9/fTTBAMA8BFet1rRf//3f+uLL77QF198oZ07d6qmpsZ+7qefftKGDRu0YcMG+7HAwEANHz7cPrrQGBp69erlqhIBwKtFR0crJSVFFotFoaGhBAMAwFlzWTiYNWuW/bHValVJSYk9LDT+/PDDD/Y2J06c0JdffqkdO3bYdwWWpIiICF1yySUaM2aM5s2b56py0QKWMgUAAPAdbpnHYzKZNHToUA0dOlTTpk2zH//hhx+aBAaz2ayGhgZ7m++++07fffedVq9eTTgAAAAAXMjQSf79+vXT1Vdfrauvvtp+zGKxaMeOHQ6BYefOnbJYLAZW6rtYyhQAAMB3eNwdwKGhobr88st1+eWX2481TksCAAAA4Dp+bTcxXuO0JAAAAACu4xXhAAAAAIDrEQ4AAAAASCIcAAAAALDxuBuS4VnY5wAAAMB3EA4AwIvFxMSorq5OgYGBRpcCAOgCCAdoFfscAJ4tJibG6BIAAF0I9xwAAAAAkEQ4AAAAAGBDOAAAAAAgiXAAAAAAwIZwAAAAAEAS4QAAAACADUuZolVsggYAAOA7GDkAAAAAIImRA7SBTdAAAAB8ByMHAAAAACQRDgAAAADYEA4AAAAASCIcAAAAALAhHAAAAACQRDgAAAAAYEM4AAAAACCJcAAAAADAhnAAAAAAQBI7JKMNsbGxDs/r6uoMqgQAAACuxsgBAAAAAEmMHKANxcXFDs8rKioUERFhUDUAAABwJUYOAAAAAEgiHAAAAACwIRwAAAAAkEQ4AAAAAGBDOAAAAAAgiXAAAAAAwIZwAAAAAEAS4QAAAACADeEAAAAAgCTCAQAAAAAbwgEAAAAASVKA0QXAs8XGxjo8r6urM6gSAAAAuBojBwAAAAAkMXKANhQXFzs8r6ioUEREhEHVAAAAwJUYOQAAAAAgiXAAAAAAwIZwAAAAAEAS4QAAAACADeEAAAAAgCTCAQAAAAAbwgEAAAAASYQDAAAAADaEAwAAAACSCAcAAAAAbAgHAAAAACQRDgAAAADYEA4AAAAASCIcAAAAALAhHAAAAACQRDgAAAAAYEM4AAAAACCJcAAAAADAJsDoAuDZYmNjHZ7X1dUZVAkAAABcjZEDAAAAAJIYOUAbiouLHZ5XVFQoIiLCoGoAAADgSowcAAAAAJBEOAAAAABgQzgAAAAAIIlwAAAAAMCGcAAAAABAEuEAAAAAgA3hAAAAAIAkwgEAAAAAG8IBAAAAAEmEAwAAAAA2hAMAAAAAkggHAAAAAGwIBwAAAAAkEQ4AAAAA2BAOAAAAAEgiHAAAAACwIRwAAAAAkEQ4AAAAAGBDOAAAAAAgiXAAAAAAwIZw4OWsVqveffddJSUlqX///urWrZuGDh2qWbNm6dtvvzW6PAAAAHgRwoGXe+ihh3Trrbfqm2++0ZQpU/Sb3/xGkZGRysrK0ujRo/XVV18ZXSIAAAC8RIDRBaDzDhw4oBdeeEGDBw/Wl19+qZ49e9rP/fnPf9Zvf/tbLVq0SDk5OQZWCQAAAG/ByIEX27t3rxoaGpSQkOAQDCTphhtukCQdOnTIiNIAAADghQgHnXTw4EGtXr1ac+fO1bXXXqu+ffvKZDLJZDIpNTW1Q9cqKyvTnDlzNGzYMIWFhalPnz6Ki4vTggULVFNT0+LroqOjFRQUpC1btqiystLh3OrVqyVJkydP7vBnAwAAgG9iWlEn9evXzynXWbVqlaZPn+7w5b6mpkaFhYUqLCxUdna2CgoKFBUV1eS155xzjp599ll7sLjpppvUo0cPffnll1q3bp1+/etf64EHHnBKnQAAAOj6CAdOMGjQIA0bNkxr167t0OuKioo0bdo0WSwWhYeH69FHH1VSUpIsFoveeOMNZWVlqaSkRNdff70KCwvVvXv3Jtd48MEHdcEFF+jee+/Vyy+/bD+emJiou+66SwEBdDEAAADah2lFnTR37lytWrVKBw4cUFlZmf72t791+BqzZ8+WxWJRQECA1q5dq8zMTI0bN06TJk3SkiVL9Nxzz0mSSkpKtHDhwmav8eSTT2r69OnKzMxUeXm5jh07pk2bNqm2tlYTJ07UBx98cFafEwAAAL7DZLVarUYX0RXs3btXkZGRkqSUlBTl5eW12n7btm26/PLLJUmzZs1y+Ff/Rg0NDbr44ou1a9cu9erVSwcPHlRgYKD9/Mcff6yrrrpKDz74oBYtWuTw2gMHDujCCy/UBRdcILPZfJaf7t8qKioUEREhSSovL9fAgQOddm0AAAC0j6u+kzFyYJCVK1faH6elpTXbxs/PTzNnzpQkHT16VOvXr3c4v2bNGklSUlJSk9eef/75GjZsmEpLS1VVVeWkqgEAANCVEQ4MsnnzZklSWFiYLrvsshbbTZgwwf54y5YtDudOnDghqeXlSg8dOiQ/Pz+H0QYAAACgJYQDg+zatUuSFBUV1epNw8OGDWvymkYJCQmSpEWLFunnn392OPfyyy+roqJC48aNU3BwsLPKBgAAQBfGUjYGqK2t1eHDhyWpzflhvXv3VlhYmKqrq1VeXu5wburUqfrv//5vffrpp4qJidGvfvUr9erVS//85z+1bt06hYaGNrkXoS0VFRWtnt+/f3+HrgcAAADvQTgwwLFjx+yPw8PD22zfGA7OvHfA399fa9eu1Z///Ge99dZbWr58uU6cOKF+/frZVzAaPnx4h2prvLEFAAAAvodwYIDa2lr746CgoDbbN04LslgszZ77/e9/r9///vfOKxAAAAA+iXBggJCQEPvjxpuKW3P8+HFJUmhoqMtqanTm1KUz7d+/X/Hx8S6vAwAAAO5HODDA6Tsdt2eZ0erqakntm4J0tti3AAAAwHexWpEBQkJCdM4550hq+wbgI0eO2MMB9wMAAADAlRg5MMiIESO0adMmlZaW6uTJky0uZ7p79277447eXOwMsbGxDs/r6urcXgMAAADcg5EDgyQmJko6NWXo888/b7Hdxo0b7Y8b9zUAAAAAXIFwYJApU6bYH+fm5jbbpqGhQUuXLpUk9erVS0lJSe4ozUFxcbHDz7p169xeAwAAANyDcGCQ+Ph4jR8/XpL0yiuv6LPPPmvSZuHChfZdkWfPnq3AwEC31ggAAADfwj0HnbR582aVlpbanzfueCxJpaWlysvLc2ifmpra5BqLFy9WQkKCLBaLkpOTlZmZqaSkJFksFr3xxhtasmSJJCkmJkZz5sxxyecAAAAAGpmsVqvV6CK8UWpqqvLz89vdvqU/5lWrVmn69OmqrKxs9nxMTIwKCgoUFRXVqTqdraKiwr5qUnl5OUufAgAAGMBV38kYOTDYjTfeqB07dmjx4sUqKChQRUWFgoKCFBUVpalTp+qBBx5Qt27dDKuP1YoAAAB8ByMHaFVz4cBsNkti5AAAAMAojBzAEMXFxQ7PT/8PEQAAAF0LqxUBAAAAkEQ4AAAAAGBDOAAAAAAgiXAAAAAAwIYbktEqljIFAADwHYwcAAAAAJDEyAHawFKmAAAAvoORAwAAAACSCAcAAAAAbAgHAAAAACQRDgAAAADYEA4AAAAASGK1IrSBfQ4AAAB8ByMHAAAAACQxcoA2sM8BAACA72DkAAAAAIAkwgEAAAAAG8IBAAAAAEmEAwAAAAA2hAMAAAAAkggHAAAAAGxYyhStYhM0AAAA38HIAQAAAABJjBygDWyCBgAA4DsYOQAAAAAgiXAAAAAAwIZwAAAAAEAS4QAAAACADeEAAAAAgCTCAQAAAAAbwgEAAAAASYQDAAAAADZsgoZWxcbGOjyvq6szqBIAAAC4GuEA8BElJSWqq6tTYGCgYmJijC4HAAB4IMIBWlVcXOzwvKKiQhEREQZVg7NRUlIii8Wi0NBQwgEAAGgW9xwAAAAAkEQ4AAAAAGBDOAAAAAAgiXAAAAAAwIZwAAAAAEAS4QAAAACADeEAAAAAgCTCAQAAAAAbwgEAAAAASeyQDPgEs9ms/Px87du3TwMGDNDQoUMVHR1tdFkAAMDDEA6ALi43N1cZGRmqr6+3H3vvvfeUlZWltLQ0AysDAACehmlFQBdmNpubBANJqq+vV0ZGhsxms0GVAQAAT8TIAVoVGxvr8Lyurs6gStAZOTk5TYJBo/r6euXk5OiZZ55xc1UAAMBTMXIAdGF79uw5q/MAAMC3MHKAVhUXFzs8r6ioUEREhEHVoKMiIyPP6jwAAPAtjBwAXVh6err8/f2bPefv76/09HQ3VwQAADwZ4QDowqKjo5WVldUkIPj7+ysrK4vlTAEAgAOmFQFdXFpamhITE5WZmWnf5+Dpp58mGAAAgCYIB4APiI6OVkpKiiwWi0JDQwkGAACgWUwrAgAAACCJcAAAAADAhnAAAAAAQBLhAAAAAIAN4QAAAACAJMIBAAAAABvCAQAAAABJhAMAAAAANoQDAAAAAJIIBwAAAABsCAcAAAAAJBEOAAAAANgQDgAAAABIIhwAAAAAsAkwugB4ttjYWIfndXV1BlUCAAAAVyMcAD4iJiZGdXV1CgwMNLoUAADgoQgHaFVxcbHD84qKCkVERBhUDc5GTEyM0SUAAAAPxz0HAAAAACQRDgAAAADYEA4AAAAASCIcAAAAALAhHAAAAACQRDgAAAAAYEM4AAAAACCJcAAAAADAhnAAAAAAQBLhAAAAAIAN4QAAAACAJMIBAAAAABvCAQAAAABJhAMAAAAANoQDAAAAAJIIBwAAAABsCAcAAAAAJBEOAAAAANgQDgAAAABIIhwAAAAAsCEceLG8vDyZTKZWfyZPnmx0mQAAAPASAUYXgM4bPXq0Hn/88WbPrVixQsXFxbr66qvdXBUAAAC8FeHAi40ePVqjR49ucvzEiRN68cUXFRAQoJSUFPcXBgAAAK/EtKIuaOXKlfrxxx91ww03qF+/fkaXAwAAAC9BOOikgwcPavXq1Zo7d66uvfZa9e3b1z7PPzU1tUPXKisr05w5czRs2DCFhYWpT58+iouL04IFC1RTU9Ph2rKzsyVJ9957b4dfCwAAAN/FtKJOcta/yK9atUrTp09XZWWl/VhNTY0KCwtVWFio7OxsFRQUKCoqql3XKysr0yeffKKBAwfqmmuucUqNAAAA8A2MHDjBoEGDlJyc3OHXFRUVadq0aaqsrFR4eLieeuopbd26VZ988okyMjIkSSUlJbr++ut17Nixdl0zNzdXDQ0NSk1Nlb+/f4drAgAAgO9i5KCT5s6dq7i4OMXFxalfv37au3evIiMjO3SN2bNny2KxKCAgQGvXrtW4cePs5yZNmqTo6Gg98sgjKikp0cKFCzVv3rxWr9fQ0KDc3FyZTCalp6d35mMBAADAhzFy0ElPPPHEWd3wu23bNm3atEmSdM899zgEg0Zz5szR8OHDJUmLFy9WXV1dq9f8+OOP9d1332nSpEkdDioAAAAA4cAgK1eutD9OS0trto2fn59mzpwpSTp69KjWr1/f6jW5ERkAAABng3BgkM2bN0uSwsLCdNlll7XYbsKECfbHW7ZsabHdjz/+qPfff199+vTRzTff7LxCAQAA4DO458Agu3btkiRFRUUpIKDlbhg2bFiT1zRn2bJlOnHihKZPn67g4OBO11VRUdHq+f3793f62gAAAPBshAMD1NbW6vDhw5KkgQMHttq2d+/eCgsLU3V1tcrLy1ts98orr0g6+ylFERERZ/V6AAAAeC+mFRng9GVJw8PD22wfFhYmSaqqqmr2/LZt2/TVV18pPj5eI0eOdE6RAAAA8DmMHBigtrbW/jgoKKjN9o3ThCwWS7Pn4+PjZbVanVJba6MT0qlpRfHx8U55LwAAAHgWwoEBQkJC7I9PnDjRZvvjx49LkkJDQ11WU6O2pjkBAACg62JakQG6d+9uf9zSVKHTVVdXS2rfFCQAAACgswgHBggJCdE555wjqe3VgY4cOWIPB9wsDAAAAFdiWpFBRowYoU2bNqm0tFQnT55scTnT3bt32x837pbsTrGxsQ7P29qlGQAAAN6LkQODJCYmSjo1Zejzzz9vsd3GjRvtjxMSElxeFwAAAHwX4cAgU6ZMsT/Ozc1ttk1DQ4OWLl0qSerVq5eSkpLcUZqD4uJih59169a5vQYAAAC4B+HAIPHx8Ro/frykUxuYffbZZ03aLFy40L4r8uzZsxUYGOjWGgEAAOBbuOegkzZv3qzS0lL788YdjyWptLRUeXl5Du1TU1ObXGPx4sVKSEiQxWJRcnKyMjMzlZSUJIvFojfeeENLliyRJMXExGjOnDku+RwAAABAI5PVWbtn+ZjU1FTl5+e3u31Lf8yrVq3S9OnTVVlZ2ez5mJgYFRQUKCoqqlN1OltFRYV91aTy8nL2RQAAADCAq76TMXJgsBtvvFE7duzQ4sWLVVBQoIqKCgUFBSkqKkpTp07VAw88oG7duhlWH6sVAQAA+A5GDtCq5sKB2WyWxMgBAACAURg5gCGKi4sdnp/+HyIAAAC6FlYrAgAAACCJcAAAAADAhnAAAAAAQBL3HKCDTp48aX+8f/9+AysBAADwXad/Dzv9+9nZIhygVWeuVlRTU2N/HB8f7+5yAAAAcIZDhw5pyJAhTrkW04rQIfX19UaXAAAAABdhnwN0SGlpqaKjoyVJW7duNWxZ00mTJkmS1q1bZ9i1OvK6ttp29nx7j+/fv98+0rNt2zb179+/zZpdwZv6rT3tWmtDnznvWvyudZzR/ebMPmurDX3mvGvxu9Yxzuyzzlzv5MmTOnTokCRp5MiRCgkJcUodTCtCh5z+H15ERIRhm6AFBgZKklPev7PX6sjr2mrb2fMdPS5J/fv3p9/a8br2tGutDX3mvGvxu9ZxRvebM/usrTb0mfOuxe9axzizzzp7PWdNJTod04oAAAAASCIcAAAAALAhHAAAAACQxA3J6KCKigr7Tcjl5eWGzc1Ex9Bv3oc+8070m/ehz7wT/eY6jBwAAAAAkEQ4AAAAAGBDOAAAAAAgiXsOAAAAANgwcgAAAABAEuEAAAAAgA3hAAAAAIAkwgEAAAAAG8IBAAAAAEmEAwAAAAA2hAMAAAAAkggH8CC1tbX67W9/qyuuuEIDBgxQSEiIzj//fCUkJCg3N1d1dXVGl4gzfP/993rhhReUnJysQYMGKSgoSOeff75uvfVW/eMf/zC6PLTi1Vdf1axZszR27FgFBwfLZDIpLy/P6LJ83vbt23XdddepV69eCgsL0y9+8Qu99dZbRpeFVvC75H34u6t1bIIGj3H48GFFREQoPj5eMTExOvfcc3XkyBGtWbNGZWVlSk5O1po1a+TnR6b1FL///e81f/58XXTRRZo4caLOPfdcmc1mrVy5UlarVcuXL9e0adOMLhPNGDJkiMrKytS3b1+FhYWprKxMubm5Sk1NNbo0n7V+/XpdffXVCgkJ0R133KHu3bvrnXfeUVlZmZ5//nnNmTPH6BLRDH6XvA9/d7XBCniI+vp66/Hjx5scr6urs06cONEqybp69WoDKkNL3nnnHeuGDRuaHP/000+tgYGB1t69e1tra2sNqAxt+eijj6x79+61Wq1W6zPPPGOVZM3NzTW2KB9WV1dnveiii6zBwcHWoqIi+/GjR49aY2JirEFBQfb+gmfhd8n78HdX6/gnWHgMPz8/BQUFNTkeEBCgm2++WZJUWlrq7rLQiltuuUUTJkxocnz8+PFKSkrSkSNHtHPnTgMqQ1uuvPJKDR482OgyYLNu3Tr961//0l133aXRo0fbj/fs2VOZmZk6ceKE8vPzjSsQLeJ3yfvwd1frCAddxMGDB7V69WrNnTtX1157rfr27SuTySSTydThoc2ysjLNmTNHw4YNU1hYmPr06aO4uDgtWLBANTU1rvkArWhoaND//u//SpIuvvhit7+/q3TlPpOkwMBASafCXVfS1fvNF3lCn27YsEGSlJyc3OTc1VdfLUnauHFjh2rp6jyh39Bxnt5vXfXvrg4xeugCziGpxZ+UlJR2X+eDDz6w9ujRo8VrxcTEWM1ms+s+iNVqPX78uPXxxx+3zp0713r//fdbhw0bZpVkTUtLc+n7ultX6rMzlZWVWYODg639+/e3njx50q3v7Wpdsd98fSqEJ/TpbbfdZpVkLSwsbPZ8eHi4NSIiojMfr8vyhH47k6//LrWHJ/Zbo678d1dHMHLQBQ0aNKjZf31qS1FRkaZNm6bKykqFh4frqaee0tatW/XJJ58oIyNDklRSUqLrr79ex44dc3bZdidOnNATTzyhJ598Un/961/1zTff6KGHHtKSJUtc9p5G8/Y+O11dXZ1mzJih48ePa/78+fL393fL+xqhK/UbTjGqT3/++WdJp6YRNadHjx72NmiK30Xv5En95kt/d7XFh8dMupa5c+cqLi5OcXFx6tevn/bu3avIyMgOXWP27NmyWCwKCAjQ2rVrNW7cOPu5SZMmKTo6Wo888ohKSkq0cOFCzZs3r8k15syZo+PHj3foPaOjox2OhYeHy2q1qqGhQfv27dOqVauUmZmpzz77TP/zP/+jHj16dOhzeaqu1GeNGhoalJqaqk8//VQZGRmaMWNGhz6PN+iK/ebrPKVP0TH0m3fyxH7zhb+7OsTooQu4xp49ezo0TPePf/zD3n7WrFnNtqmvr7cOHz7cKsnaq1cv64kTJ5q0CQsLa3XI8Myf9evXt+vzvPXWW1ZJ1kceeaRd7b2Rt/dZfX29NSUlxSrJOn36dGt9fX1HPr7X8vZ+s1qZCnEmI/qUaUVnz6jfxdPxu9RxRvebr/7d1RqmFUGStHLlSvvjtLS0Ztv4+flp5syZkqSjR49q/fr1TdpUVVXJarW2+2fixIntqq9x2LHxpj14Vp81NDQoLS1N+fn5uvPOO5WXl8d+FC3wpH6DczijTxtHdcxmc5PXHjhwQFVVVYz8OJmzfhfhXs7sN/7uah5/ApAkbd68WZIUFhamyy67rMV2py/9tWXLFpfX1Wjfvn2S/r2KADynzxr/z3Xp0qWaNm2ali1b5tNzNdviKf0G53FGnzaeW7t2bZPXffjhh01ej7PH76J3cla/8XdXywgHkCTt2rVLkhQVFdXq8l3Dhg1r8hpn+frrr5tdeqympka//e1vJUnXXXedU9/Tm3lCnzU0NCg9PV1Lly7V1KlT9eqrr/J/rm3whH6DczmjTydPnqwLL7xQy5cv1xdffGE//vPPP+vpp59WUFCQ/V9C4Rz8LnonZ/Qbf3e1jhuSodraWh0+fFiSNHDgwFbb9u7dW2FhYaqurlZ5eblT63jrrbe0aNEiJSYmasiQIerRo4e+//57rVmzRj/++KPGjx+vBx980Knv6a08pc+efPJJ5efnKzw8XDExMfrTn/7UpM2UKVMcNnXyZZ7Sb5KUnZ1t/xe4xs1+srOz7VP3EhMTde+99zr9fbsaZ/VpQECAsrOzdfXVV+uKK67QHXfcoe7du+udd95RWVmZnn/+eQ0ZMsRVH8PnOPN3kd8l93FWv/F3V+sIB3BY5is8PLzN9o2/bFVVVU6t44YbbtC+ffu0detWffbZZ6qqqlLPnj01atQo3XHHHUpPT/ftTUlO4yl9tnfvXkmn5r8/9dRTzbYZMmSIz/4f7Jk8pd+kU0PzZ+64u2XLFofhd77QtM2ZfZqUlKTNmzfr8ccf15tvvqm6ujqNHDlS8+fP17Rp05xat69zZr/xu+Q+zuo3/u5qHd+0oNraWvvjoKCgNtsHBwdLkiwWi1PrGDt2rMaOHevUa3ZVntJneXl5ysvLc+o1uzJP6TeJvnMWZ/dpfHy81qxZ45zi0CJn9hu/S+7jrH6jz1rHPQdQSEiI/fGJEyfabN+4tnpoaKjLakLr6DPvRL91PfSpd6LfvBP95h6EA6h79+72x+2ZvlBdXS2pfUN6cA36zDvRb10Pfeqd6DfvRL+5B+EACgkJ0TnnnCNJqqioaLXtkSNH7L9sERERLq8NzaPPvBP91vXQp96JfvNO9Jt7EA4gSRoxYoQkqbS0VCdPnmyx3e7du+2Phw8f7vK60DL6zDvRb10Pfeqd6DfvRL+5HuEAkk4ttSadGoL7/PPPW2y3ceNG++OEhASX14WW0WfeiX7reuhT70S/eSf6zfUIB5B0aj3fRrm5uc22aWho0NKlSyVJvXr1UlJSkjtKQwvoM+9Ev3U99Kl3ot+8E/3meoQDSDq1fN748eMlSa+88oo+++yzJm0WLlxo32Vw9uzZCgwMdGuNcESfeSf6reuhT70T/ead6DfXM1mtVqvRReDsbd68WaWlpfbnhw8f1sMPPyzp1HDamRuwpKamNrlGUVGREhISZLFYFB4erszMTCUlJcliseiNN97QkiVLJEkxMTEqLCx0WDUAHUefeSf6reuhT70T/ead6DcvYEWXkJKSYpXU7p+WfPDBB9YePXq0+LqYmBir2Wx24yfruugz70S/dT30qXei37wT/eb5mFYEBzfeeKN27NihBx98UDExMerWrZt69eqlsWPHav78+SoqKlJUVJTRZeI09Jl3ot+6HvrUO9Fv3ol+cx2mFQEAAACQxA3JAAAAAGwIBwAAAAAkEQ4AAAAA2BAOAAAAAEgiHAAAAACwIRwAAAAAkEQ4AAAAAGBDOAAAAAAgiXAAAAAAwIZwAAAAAEAS4QAAAACADeEAAAAAgCTCAQAAAAAbwgEAAAAASYQDAAAAADaEAwAAAACSCAcAAAAAbAgHAAAAACQRDgAAAADYEA4AAB4tNTVVJpOpyc/evXubtJ03b579vDfasGFDs5913rx5RpcGwEcQDgAAAABIkgKMLgAAgPYYMGCAPvzwQ/vzCy64wMBqXCMuLk47d+60Px85cqSB1QDwRYQDAIBXCAwM1MUXX2x0GS4VFhbW5T8jAM/GtCIAAAAAkggHAAAAAGwIBwAAp9u9e7d9pZ3XX39dVqtVr732mpKTk9WvXz/5+flp4sSJhtT2xRdfqF+/fjKZTOrfv7927NhhP3fmakdHjx7V448/rtjYWIWHh6tPnz5KSkrS66+/3q732rJli+69914NHTpUPXr0UFBQkAYOHKgbbrhBf/3rX3X06FFXfEQA6DTuOQAAON2XX35pf3z++edrwoQJ2rRpk0ObUaNGubssbdq0STfeeKN+/vlnDRkyRB9//LEuuuiiZtvu2bNHV111lf71r3/Zj1VXV2vDhg3asGGDVq5cqddee00BAU3/KrVYLLrnnnuaDRHff/+9vv/+exUUFOjQoUMsUwrAoxAOAABOd3o4+M1vfqOvv/5ad911l+644w4NGDBA5eXlOu+889xaU0FBgaZOnSqLxaLY2FitXbtWAwYMaLH9tGnTtGfPHt1333267bbb1LNnT+3YsUPz589XSUmJ3nrrLQ0YMEB//vOfHV7X0NCgm266SR999JEkKTo6Wr/+9a81duxYdevWTfv379fWrVv11ltvufTzAkBnEA4AAE53ejj417/+pQ8++EA33HCD/dhll13m1nqWL1+ulJQUnTx5UvHx8VqzZo369OnT6mu2b9+u5cuX684777QfGzt2rKZOnarx48fryy+/1H/913/pnnvucVhh6MUXX7QHg5tvvlmvv/66goODHa59/fXX6z//8z+1f/9+J35KADh73HMAAHC6L774wv74r3/9q0MwcLeXXnpJ06dP18mTJzV58mR98sknbQYDSbrhhhscgkGj7t27a8mSJZJOjRK8/PLL9nMNDQ1asGCBJGngwIFaunRpk2DQyM/Pr0vu1QDAuxEOAABOdfjwYe3bt0/SqU290tPTDavlT3/6k+6//35ZrVbdfPPNKigoUHh4eLtem5aW1uK5+Ph4xcbGSpI+/vhj+/EvvvhCFRUVkqSMjIx2vxcAeArCAQDAqU6fUnTfffcZVseDDz6oP/7xj5JOfdF/++23W/xX/ObExcW1ej4+Pl6SVFJSohMnTkiSioqK7OfHjx/f0ZIBwHCEAwCAU50eDq6++mrD6njhhRckSRdffLGys7Pl7+/fode3dcN0v379JElWq1VHjhyRdGrUpFH//v079H4A4AkIBwAAp2q83+CCCy4wdE79rbfeKkn66quvNHv27A6/vnGvAwDwJYQDAIBTNY4cjBkzxtA6Xn/9dU2ZMkXSqRWEHnzwwQ69/ocffmjXeZPJpN69e0uS+vbtaz/PSkQAvBHhAADgNCdOnNCuXbskSaNHjza0lsDAQL355pv2lZJeeOEFPfzww+1+/fbt29t1Pjo6WkFBQZKkSy+91H7+008/7WjJAGA4wgEAwGl27dqluro6ScaPHEhSUFCQ3nnnHV133XWSpOeff16///3v2/Xa/Pz8Fs9t375dX331lSTpyiuvtB+/5JJLFBERIUnKzs5WVVVVZ0sHAEMQDgAATnP6/gaeEA6kUwHh3Xfftd8cPX/+fD322GNtvu6DDz5odhfjqqoqzZo1S9KpvQoaHzc+bxydqKio0MyZM+0rGZ2poaHBvuQrAHgKwgEAwGka7zfo2bOnIiMjDa7m34KDg7Vy5UpdddVVkqSnnnpKjz/+eKuvGTt2rO666y7df//9Wr9+vT7//HPl5uZq7Nix9iVL77//fo0aNcrhdffff7/9fd577z2NHDlSixcv1pYtW1RUVKQ1a9bo8ccf17Bhw+ybqQGApwgwugAAQNfRGA6Mvt+gOSEhIXr//fd1ww03aN26dXryyScVGBjY4ijCW2+9pcmTJ+ull17SSy+91OT8rbfeqkWLFjU57ufnp5UrVyolJUUrVqxQSUmJ/t//+3/O/jgA4BKMHAAAnMaTw4EkhYaGatWqVZowYYIk6Y9//KOeeeaZZttGRkbq888/V2ZmpoYPH65u3bqpZ8+euuKKK/Tqq69qxYoVCgho/t/YunXrprffflvr1q3TjBkzFBkZqdDQUAUFBSkiIkI33nij/va3v2nOnDku+6wA0BmMHAAAnOb0TcCMMG/ePM2bN6/VNt26ddOGDRvadb3evXvrqaee0lNPPdWpepKSkpSUlNSp1wKAEQgHAACvUFdXZ18hSJKGDh2qwMBAAytyvurqau3Zs8foMgD4MMIBAMAr7Nu3TyNHjrQ/37Nnj4YMGWJcQS6wfft2RhoAGIp7DgAAAABIYuQAAODh8vLylJeXZ3QZbjFx4kRZrVajywDgwxg5AAAAACBJMln5JwoAAAAAYuQAAAAAgA3hAAAAAIAkwgEAAAAAG8IBAAAAAEmEAwAAAAA2hAMAAAAAkggHAAAAAGwIBwAAAAAkEQ4AAAAA2BAOAAAAAEgiHAAAAACwIRwAAAAAkEQ4AAAAAGBDOAAAAAAgiXAAAAAAwIZwAAAAAEAS4QAAAACADeEAAAAAgCTCAQAAAAAbwgEAAAAASdL/B2Fffnz+PUKcAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 384, "width": 387 } }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n", "\n", "ax.errorbar(\n", " tbl[\"r\"],\n", " tbl[\"Menc\"],\n", " yerr=(tbl[\"Menc_err_neg\"], tbl[\"Menc_err_pos\"]),\n", " marker=\"o\",\n", " markersize=2,\n", " color=\"k\",\n", " alpha=1.0,\n", " ecolor=\"#aaaaaa\",\n", " capthick=0,\n", " linestyle=\"none\",\n", " elinewidth=1.0,\n", ")\n", "\n", "ax.set_xlim(1e-3, 10**2.6)\n", "ax.set_ylim(7e6, 10**12.25)\n", "\n", "ax.set_xlabel(\"$r$ [kpc]\")\n", "ax.set_ylabel(\"$M(" ] }, "metadata": { "image/png": { "height": 384, "width": 387 } }, "output_type": "display_data" } ], "source": [ "xyz = np.zeros((3, 256))\n", "xyz[0] = np.logspace(-3, 3, 256)\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n", "\n", "ax.errorbar(\n", " tbl[\"r\"],\n", " tbl[\"Menc\"],\n", " yerr=(tbl[\"Menc_err_neg\"], tbl[\"Menc_err_pos\"]),\n", " marker=\"o\",\n", " markersize=2,\n", " color=\"k\",\n", " alpha=1.0,\n", " ecolor=\"#aaaaaa\",\n", " capthick=0,\n", " linestyle=\"none\",\n", " elinewidth=1.0,\n", ")\n", "\n", "fit_menc = init_potential.mass_enclosed(xyz * u.kpc)\n", "ax.loglog(xyz[0], fit_menc.value, marker=\"\", color=\"#3182bd\", linewidth=2, alpha=0.7)\n", "\n", "ax.set_xlim(1e-3, 10**2.6)\n", "ax.set_ylim(7e6, 10**12.25)\n", "\n", "ax.set_xlabel(\"$r$ [kpc]\")\n", "ax.set_ylabel(\"$M(" ] }, "metadata": { "image/png": { "height": 384, "width": 387 } }, "output_type": "display_data" } ], "source": [ "xyz = np.zeros((3, 256))\n", "xyz[0] = np.logspace(-3, 3, 256)\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n", "\n", "ax.errorbar(\n", " tbl[\"r\"],\n", " tbl[\"Menc\"],\n", " yerr=(tbl[\"Menc_err_neg\"], tbl[\"Menc_err_pos\"]),\n", " marker=\"o\",\n", " markersize=2,\n", " color=\"k\",\n", " alpha=1.0,\n", " ecolor=\"#aaaaaa\",\n", " capthick=0,\n", " linestyle=\"none\",\n", " elinewidth=1.0,\n", ")\n", "\n", "fit_menc = fit_potential.mass_enclosed(xyz * u.kpc)\n", "ax.loglog(xyz[0], fit_menc.value, marker=\"\", color=\"#3182bd\", linewidth=2, alpha=0.7)\n", "\n", "ax.set_xlim(1e-3, 10**2.6)\n", "ax.set_ylim(7e6, 10**12.25)\n", "\n", "ax.set_xlabel(\"$r$ [kpc]\")\n", "ax.set_ylabel(\"$M(" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "potential = MilkyWayPotential()\n", "potential" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.10.15" } }, "nbformat": 4, "nbformat_minor": 5 }