{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "4b1cde21", "metadata": { "execution": { "iopub.execute_input": "2025-12-11T03:38:06.427057Z", "iopub.status.busy": "2025-12-11T03:38:06.426877Z", "iopub.status.idle": "2025-12-11T03:38:06.821759Z", "shell.execute_reply": "2025-12-11T03:38:06.821011Z" }, "nbsphinx": "hidden" }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "id": "8c648b64", "metadata": { "execution": { "iopub.execute_input": "2025-12-11T03:38:06.823922Z", "iopub.status.busy": "2025-12-11T03:38:06.823699Z", "iopub.status.idle": "2025-12-11T03:38:06.830653Z", "shell.execute_reply": "2025-12-11T03:38:06.829935Z" }, "nbsphinx": "hidden" }, "outputs": [], "source": [ "%run ../tutorials/nb_setup" ] }, { "cell_type": "markdown", "id": "648eb16a", "metadata": {}, "source": [ "# Defining the MilkyWayPotential model\n", "\n", "## Introduction\n", "\n", "`gala` provides simplified mass models for the Milky Way to use in orbit integration or dynamical calculations. Some of these mass models come from other publications or packages (e.g., the Law and Majewski 2010 model `LM10Potential`). Some of the potential models are defined and provided by Gala. This document describes how we determined the parameters of the Gala Milky Way models.\n", "\n", "We determine parameters of the Gala Milky Way models using compilations of enclosed mass measurements of the Milky Way and measurements of the mass structure of the Galactic disk. We then fit for the parameters of a multi-component model (e.g., disk, bulge, halo, etc.) using these measurements." ] }, { "cell_type": "code", "execution_count": 3, "id": "0045f910", "metadata": { "execution": { "iopub.execute_input": "2025-12-11T03:38:06.832572Z", "iopub.status.busy": "2025-12-11T03:38:06.832403Z", "iopub.status.idle": "2025-12-11T03:38:07.724374Z", "shell.execute_reply": "2025-12-11T03:38:07.723510Z" } }, "outputs": [], "source": [ "import astropy.units as u\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from astropy.constants import G\n", "from astropy.io import ascii\n", "from scipy.optimize import leastsq\n", "\n", "import gala.potential as gp\n", "from gala.units import galactic" ] }, { "cell_type": "markdown", "id": "101d32c0", "metadata": {}, "source": [ "## `MilkyWayPotential` version 1 (circa 2017)\n", "\n", "This model was previously just known as `MilkyWayPotential` in Gala, now known as \"version 1,\" and represents an older model based on measurements that are now out of date. We still describe the process of fitting for this model, for completeness.\n", "\n", "The source data for this model was compiled from published values and is included with Gala:" ] }, { "cell_type": "code", "execution_count": 4, "id": "c70283af", "metadata": { "execution": { "iopub.execute_input": "2025-12-11T03:38:07.726223Z", "iopub.status.busy": "2025-12-11T03:38:07.725942Z", "iopub.status.idle": "2025-12-11T03:38:07.737437Z", "shell.execute_reply": "2025-12-11T03:38:07.736758Z" } }, "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": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mwdata1 = ascii.read(\"data/MW_mass_enclosed.csv\")\n", "mwdata1" ] }, { "cell_type": "markdown", "id": "5b45f1f0", "metadata": {}, "source": [ "We can now plot the above data and uncertainties:" ] }, { "cell_type": "code", "execution_count": 5, "id": "5db6f754", "metadata": { "execution": { "iopub.execute_input": "2025-12-11T03:38:07.739338Z", "iopub.status.busy": "2025-12-11T03:38:07.739141Z", "iopub.status.idle": "2025-12-11T03:38:08.466998Z", "shell.execute_reply": "2025-12-11T03:38:08.466098Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABJsAAAMMCAYAAAD0M6ZoAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAewgAAHsIBbtB1PgAAaA5JREFUeJzt3Xt4X1WdL/53er9CuVmoDVBtQksAi9A6TItQ0CICIypYdAq9QA+ecebhx1QZ7SiiZ7gNgnJGPR5aekMREQSFDkeGqxQYoU4RDMWkM21taBFQsLRNadrk9wc2Q+i93ck3SV+v58nD3t+9vmt/vmE9IXmz1tplTU1NTQEAAACAAnQpdQEAAAAAdB7CJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDdSl0AHcv69evz3HPPJUkOOuigdOtmCAEAAEBHtXHjxrzyyitJkqOPPjq9evXa4z4lBeyS5557LqNGjSp1GQAAAEDBnnrqqYwcOXKP+7GMDgAAAIDCmNnELjnooIOaj5966qkccsghJawGAACA9uaBBx7I+vXr06tXr3zoQx8qdTnswKpVq5pXML39b/49IWxil7x9j6ZDDjkkgwcPLmE1AAAAtDcHHnhg6uvr07t3b38zdjBF7ctsGR0AAAAAhRE2AQAAAFAYYRMAAAAAhRE2AQAAAFAYYRMAAAAAhfE0OgAAAKAQtbW1mTt3blauXJlBgwbliCOOSEVFRanLoo0JmwAAAIA9Nnv27EydOjWbNm1qfu2uu+7KjBkzMnny5BJWRluzjA4AAADYI7W1tVsETUmyadOmTJ06NbW1tSWqjFIQNgEAAAB7ZNasWVsETZtt2rQps2bNauOKKCVhEwAAALBHli5dukfX6VyETQAAAMAeGTJkyB5dp3MRNgEAAAB7ZMqUKenatetWr3Xt2jVTpkxp44ooJWETAAAAsEcqKioyY8aMLQKnrl27ZsaMGamoqChRZZSCsKmd+f73v5+LL744xx9/fHr27JmysrLMmTNnq21ffPHFfOtb38q4ceNy6KGHpkePHjn44IPzyU9+Mr/85S/btnAAAAD2apMnT87ixYtzzjnn5C//8i9zzjnnZPHixZk8eXKpS6ONdSt1AbT05S9/OcuXL8+BBx6YQw45JMuXL99m23/5l3/Jtddem/e+970ZN25cDjrooNTW1ubuu+/O3XffnVtvvTXjx49vw+oBAADYm1VUVGTixImpr69P7969zWjaS5nZ1M7MnDkzy5YtyyuvvJLPfvaz2207atSoPPLII1myZElmzpyZq6++OnfccUcefvjhdO3aNf/zf/7PvPnmm21UOQAAAICwqd350Ic+lMMOO2yn2n7iE5/ISSedtMXrJ554YsaOHZvXXnstzz33XNElAgAAAGyTsOnPXn755dx77725/PLLc/rpp+fAAw9MWVlZysrKMmnSpF3qa/ny5Zk2bVqGDRuWvn37Zv/998/IkSNz3XXXZd26da3zAd6he/fuSZJu3ayUBAAAANqOJOLPBg4cWEg/99xzTyZMmJDVq1c3v7Zu3bosXLgwCxcuzMyZMzN//vwMHTq0kPttze9+97s88MADOeSQQ3L00Ue32n0AAAAA3snMpq049NBDM27cuF1+36JFizJ+/PisXr06/fr1y5VXXpknnngiDz74YKZOnZokqampyRlnnJE33nij6LKTJA0NDTn//PPz5ptv5tprr93isZMAAAAArcnMpj+7/PLLM3LkyIwcOTIDBw7MsmXLMmTIkF3q45JLLkl9fX26deuW+++/PyeccELztVNOOSUVFRW57LLLUlNTk+uvvz5XXHFFoZ+hsbExkyZNyi9+8YtMnTo1559/fqH9AwAAAOyImU1/9rWvfS1nnnnmbi+ne+qpp/LYY48lSS688MIWQdNm06ZNy/Dhw5MkN954YxoaGna/4HdobGzMlClTcuutt2bChAn53ve+V1jfAAAAADtL2FSQu+++u/l48uTJW23TpUuXXHDBBUmS119/PQ8//HAh925sbMzkyZMzd+7cfPrTn86cOXPSpYt/tQAAAEDbk0gUZMGCBUmSvn375rjjjttmu5NOOqn5+PHHH9/j+24OmubNm5fx48fnlltusU8TAAAAUDLCpoIsXrw4STJ06NB067btrbCGDRu2xXt21+alc/Pmzcu5556b73//+4ImAAAAoKRsEF6A9evX59VXX02SDB48eLtt99tvv/Tt2zdr167NihUrtrg+c+bM5llSzz33XPNrjzzySJJkzJgxueiii5IkX//61zN37tz069cvlZWV+ad/+qct+jv77LMzYsSInf4sdXV1272+atWqne4LAAAA2PsImwrwxhtvNB/369dvh+03h01r1qzZ4tqCBQsyd+7cFq89/vjjLZbcbQ6bli1bliRZs2ZNrrzyyq3e6/DDD9+lsKm8vHyn2wIAAAC8k7CpAOvXr28+7tGjxw7b9+zZM0lSX1+/xbU5c+Zkzpw5O3XfXWkLAAAAe4Oampo0NDSke/fuqays7DB9dybCpgL06tWr+XjDhg07bP/mm28mSXr37t1qNe2urS3te7tVq1Zl1KhRbVQNAAAA7JqamprU19end+/erRI2tVbfnYmwqQD9+/dvPt7a0rh3Wrt2bZKdW3LX1na05xQAAADA9ngaXQF69eqVAw44IMmON9h+7bXXmsMm+yMBAAAAnY2ZTQU58sgj89hjj2XJkiXZuHFjunXb+rf2hRdeaD4ePnx4W5W326qqqlqcNzQ0lKgSAAAAoCMws6kgY8aMSfLWErlf/epX22z36KOPNh+PHj261esCAAAAaEvCpoKcffbZzcezZ8/eapvGxsbMmzcvSTJgwICMHTu2LUrbI9XV1S2+HnrooVKXBAAAALRjwqaCjBo1KieeeGKS5Oabb86TTz65RZvrr78+ixcvTpJccskl6d69e5vWCAAAANDa7Nn0ZwsWLMiSJUuaz1999dXm4yVLlmTOnDkt2k+aNGmLPm688caMHj069fX1GTduXKZPn56xY8emvr4+t912W2666aYkSWVlZaZNm9YqnwMAAACglIRNfzZz5szMnTt3q9cef/zxPP744y1e21rYdOyxx+ZHP/pRJkyYkNWrV2f69OlbtKmsrMz8+fPTv3//QuoGAAAAaE8soyvYWWedlWeffTaXXnppKisr06dPnwwYMCDHH398rr322ixatChDhw4tdZkAAAAArcLMpj+bM2fOFkvldtdhhx2WG264ITfccEMh/ZVSVVVVi/OGhoYSVQIAANB51NTUpKGhId27d09lZWWpy+k0amtrM3fu3KxcuTKDBg3KEUcckYqKilKXtdcRNgEAAEAbq6mpSX19fXr37t1uwqaOHoDNnj07U6dOzaZNm5pfu+uuuzJjxoxMnjy5hJXtfYRNbFd1dXWL87q6upSXl5eoGgAAAFpLewzAdlZtbe0WQVOSbNq0KVOnTs2YMWPMcGpD9mwCAAAAOrRZs2ZtETRttmnTpsyaNauNK9q7mdkEAAAAe7mOvtfR0qVL9+g6xRI2AQAAwF6sM+x1NGTIkD26vjM6eiDXliyjAwAAgL3UjvY6qq2tLVFlu2bKlCnp2rXrVq917do1U6ZM2aP+Z8+eneHDh+eOO+7IE088kTvuuCPDhw/P7Nmz96jfzkrYxHZVVVW1+DrllFNKXRIAAAAF6Sx7HVVUVGTGjBlbBE5du3bNjBkz9mgGUmcJ5NqSsAkAAAD2Up1pr6PJkydn8eLFOeecc/KXf/mXOeecc7J48eI9XgrYWQK5tmTPJrarurq6xXldXV3Ky8tLVA0AAABFaou9jtpSRUVFJk6cmPr6+vTu3buQPZU6UyDXVsxsAgAAgL1Ua+91tCtqampSXV2dmpqaNrvnzuhsgVxbEDYBAADAXqo19zraVTU1NXn++efbXdjUngK5jkLYBAAAAHux1trrqLNoT4FcR2HPJgAAANjLtcZeR53J5MmTM2bMmEyfPj0rV67MoEGDctVVV/k+bYOwCQAAAGAHBHI7T9jEdlVVVbU4b2hoKFElAAAAQEdgzyYAAAAACmNmE9tVXV3d4ryuri7l5eUlqgYAAKDjq62tzdy5c5v3/jniiCP2+iVZviedi7AJAAAA2sjs2bNz0UUXpbGxsfm1n/zkJ5k5c+Ze+/S32bNnZ+rUqdm0aVPza3fddVdmzJix135POjrL6AAAAKAN1NbW5sILL2wRNCVJY2NjLrzwwtTW1paostKpra3dImhKkk2bNmXq1Kl75fekMxA2AQAAQBv4xje+kaampq1ea2pqyje+8Y02rqj0Zs2atUXQtNmmTZsya9asNq6IIgibAAAAoA089thje3S9M1q6dOkeXad9EjYBAAAAJTFkyJA9uk77JGwCAACANnDiiSfu0fXOaMqUKenatetWr3Xt2jVTpkxp44oogrCJ7aqqqmrxdcopp5S6JAAAgA7p85//fLp02fqf4V26dMnnP//5Nq6o9CoqKjJjxowtAqeuXbtmxowZqaioKFFl7IlupS4AAAAA9gYVFRWZOXPmFk9f62zBSmVlZRoaGtK9e/edaj958uSMGTMm06dPz8qVKzNo0KBcddVVu/392NX7UzxhE9tVXV3d4ryuri7l5eUlqgYAAKBj2xyszJo1K0uXLs2QIUMyZcqUThM0JW+FPbuqoqIiEydOTH19fXr37r1H34/duT/FEjYBAABAG6qoqMjVV19d6jKg1QibAAAAoI3V1NQ0L/UyE4fORtgEAAAAbaympqZ5ydiuhk2CKto7YRMAAAC0odra2syePTsvvfRSDj744BxxxBG7tEfRngRV0BaETQAAANBGZs+evcXT6H76059mxowZmTx5cgkrg+J0KXUBAAAAsDeora3dImhKkk2bNmXq1Kmpra0tUWVQLGETAAAAtIFZs2ZtETRttmnTpsyaNWuHfdTW1mbu3Ln51re+lblz5wqoaJcsowMAAIA2sHTp0j26vrUleHfddZcleLQ7ZjYBAABAG9h33313+7oleHQkZjaxXVVVVS3OGxoaSlQJAADA3mtnluBdffXVbVwVbJ2ZTQAAANAG/vSnP+329T1dggdtycwmtqu6urrFeV1dXcrLy0tUDQAAQMc1ZMiQ3b6+J+/dWZWVlWloaEj37t33uC/2bmY2AQAAQBuYMmVKunbtutVrXbt2zZQpU1rlvTursrIyVVVVqays3OO+2LsJmwAAAKANVFRUZMaMGVuERl27ds2MGTNSUVHRKu+FtmYZHQAAALSRyZMnZ8yYMZk+fXpWrlyZQYMG5aqrrtqpsGhP3gttSdgEAAAAbaiioiITJ05MfX19evfuvUth0Z68F9qKZXQAAAAAFEbYBAAAAEBhLKMDAAAASq6ysjINDQ3p3r17qUthDwmbAAAAgJKrrKwsdQkURNgEAAAAHYgZQLR3wiYAAADoQMwAor2zQTgAAAAAhRE2AQAAAFAYy+jYrqqqqhbnDQ0NJaoEAAAA6AjMbAIAAACgMGY2sV3V1dUtzuvq6lJeXl6iagAAADoHT5SjMxM2AQAAQBvzRDk6M8voAAAAACiMsAkAAACAwlhGBwAAALAT7LW1c4RNAAAAADvBXls7xzI6AAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMDYIBwAAYI/V1NQ0P6XLJsqwdxM2AQAAsMdqampSX1+f3r17C5tgL2cZHQAAAACFETYBAAAAUBjL6AAAANgjtbW1mTt3blauXJlBgwbliCOOSEVFRanLAkpE2AQAAMBumz17dqZOnZpNmzY1v3bXXXdlxowZmTx5cgkrA0rFMjoAAAB2S21t7RZBU5Js2rQpU6dOTW1tbYkqA0pJ2AQAAMBumTVr1hZB02abNm3KrFmz2rgioD0QNgEAALBbli5dukfXgc7Jnk1sV1VVVYvzhoaGElUCAAC0N0OGDNmj60DnZGYTAAAAu2XKlCnp2rXrVq917do1U6ZMaeOKgPbAzCa2q7q6usV5XV1dysvLS1QNAADQnlRUVGTGjBlbbBLetWvXzJgxIxUVFSWsDigVYRMAAAC7bfLkyRkzZkymT5+elStXZtCgQbnqqqsETbAXEzYBAACwRyoqKjJx4sTU19end+/egibYy9mzCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKEy3UhcAAABAx1dZWZmGhoZ079691KUAJSZsAgAAYI9VVlaWugSgnbCMDgAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywqZ35/ve/n4svvjjHH398evbsmbKyssyZM6ew9gAAAACtqVupC6ClL3/5y1m+fHkOPPDAHHLIIVm+fHmh7QEAAABak5lN7czMmTOzbNmyvPLKK/nsZz9beHsAAACA1mRmUzvzoQ99qFXbAwAAALQmM5v+7OWXX869996byy+/PKeffnoOPPDAlJWVpaysLJMmTdqlvpYvX55p06Zl2LBh6du3b/bff/+MHDky1113XdatW9c6HwAAAACgHTCz6c8GDhxYSD/33HNPJkyYkNWrVze/tm7duixcuDALFy7MzJkzM3/+/AwdOrSQ+wEAAKVTW1ubWbNm5dlnn80hhxyST33qUxk3blypywIoKTObtuLQQw/drf9ALFq0KOPHj8/q1avTr1+/XHnllXniiSfy4IMPZurUqUmSmpqanHHGGXnjjTeKLhsAAGhDs2fPzvDhw3PNNdfkX//1X3PzzTfn9NNPz+zZs0tdGkBJmdn0Z5dffnlGjhyZkSNHZuDAgVm2bFmGDBmyS31ccsklqa+vT7du3XL//ffnhBNOaL52yimnpKKiIpdddllqampy/fXX54orrij4UwAAAG2htrY2U6dOzaZNm1q83tjYmKlTp2bMmDGpqKgoUXUApWVm05997Wtfy5lnnrnby+meeuqpPPbYY0mSCy+8sEXQtNm0adMyfPjwJMmNN96YhoaG3S8YAAAomVmzZm0RNG22adOmzJo1q40rAmg/hE0Fufvuu5uPJ0+evNU2Xbp0yQUXXJAkef311/Pwww+3RWkAAEDBli5dukfXATozYVNBFixYkCTp27dvjjvuuG22O+mkk5qPH3/88VavCwAAKN6OttywRyuwNxM2FWTx4sVJkqFDh6Zbt21vhTVs2LAt3gMAAHQsU6ZMSZcu2/5z6uc//3lqa2vbsCKA9sMG4QVYv359Xn311STJ4MGDt9t2v/32S9++fbN27dqsWLFii+szZ85sniX13HPPNb/2yCOPJEnGjBmTiy66aLfb70hdXd12r69atWqn+wIAgM6qoqIiH/nIR/Kv//qvW72+ed+mq6++uo0rAyg9YVMB3j5Ftl+/fjtsvzlsWrNmzRbXFixYkLlz57Z47fHHH2+x5O7t4dGutt+R8vLynW4LAAB7s/79+2/3un2bgL2VsKkA69evbz7u0aPHDtv37NkzSVJfX7/FtTlz5mTOnDk7fe9dbQ8AABRjR/s27eg6QGclbCpAr169mo83bNiww/ZvvvlmkqR3796tVtPu2trSvrdbtWpVRo0a1UbVAABA+zVlypRcd9112bRp0xbXunbtmilTppSgKoDSEzYV4O3TZ7e2NO6d1q5dm2Tnlty1tR3tOQUAAPy30047Lffdd1+ampqaX+vatWtmzJiRioqKElYGUDrCpgL06tUrBxxwQP7whz/scIPt1157rTlssj8SAAB0TLNnz87UqVO3mNV0/PHH59ZbbxU0AXs1YVNBjjzyyDz22GNZsmRJNm7cmG7dtv6tfeGFF5qPhw8f3lbl7baqqqoW5w0NDSWqBAAA2ofa2tqtBk1J8h//8R8lqAigfelS6gI6izFjxiR5a4ncr371q222e/TRR5uPR48e3ep1AQAAxZo1a9ZWg6YkaWxszKxZs9q4IoD2RdhUkLPPPrv5ePbs2Vtt09jYmHnz5iVJBgwYkLFjx7ZFaXukurq6xddDDz1U6pIAAKCkli5dukfXATo7YVNBRo0alRNPPDFJcvPNN+fJJ5/cos3111+fxYsXJ0kuueSSdO/evU1rBAAA9tyQIUP26DpAZ2fPpj9bsGBBlixZ0nz+6quvNh8vWbIkc+bMadF+0qRJW/Rx4403ZvTo0amvr8+4ceMyffr0jB07NvX19bntttty0003JUkqKyszbdq0VvkcAABA65oyZUquu+66rS6l69q1a6ZMmVKCqgDaj7Kmtz+jcy82adKkzJ07d6fbb+vbds8992TChAlZvXr1Vq9XVlZm/vz5GTp06G7VWWp1dXXNT9FbsWJFBg8eXOKKAACg7W3taXRdu3bNjBkzMnny5BJWBrBrWuPvfDObCnbWWWfl2WefzY033pj58+enrq4uPXr0yNChQ3Puuefmb//2b9OnT59Sl7nTPI0OAAC2NHny5IwZMyazZs3K0qVLM2TIkEyZMiUVFRWlLg2g5MxsYru2FjbV1tYmMbMJAAAAOjozm2hz1dXVLc7fPggBAAAA3snT6AAAAAAojLAJAAAAgMIImwAAAAAojLAJAAAAgMLYIJzt2trT6AAAAAC2xcwmAAAAAApjZhPbVV1d3eK8rq4u5eXlJaoGAAAAaO/MbAIAAACgMMImAAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMJ5Gx3ZVVVW1OG9oaChRJQAAAEBHYGYTAAAAAIUxs4ntqq6ubnFeV1eX8vLyElUDAAAAtHdmNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQmG6lLoD2raqqqsV5Q0NDiSoBAAAAOgIzmwAAAAAojJlNbFd1dXWL87q6upSXl5eoGgAAAKC9M7MJAAAAgMIImwAAAAAojLAJAAAAgMIImwAAAAAojLAJAAAAgMIImwAAAAAojLAJAAAAgMIImwAAAAAoTLdSF0D7VlVV1eK8oaGhRJUAAAAAHYGZTQAAAAAUxswmtqu6urrFeV1dXcrLy0tUDQAAANDemdkEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAUplupC6B9q6qqanHe0NBQokoAAACAjsDMJgAAAAAKY2YT21VdXd3ivK6uLuXl5SWqBgAAAGjvzGwCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDdSl0A7VtVVVWL84aGhhJVAgAAAHQEZjYBAAAAUBgzm9iu6urqFud1dXUpLy8vUTUAAABAe2dmEwAAAACFETYBAAAAUBhhEwAAAACFETYBAAAAUBhhEwAAAACFaZWn0X39619vjW63cPnll7fJfQAAAADYOWVNTU1NRXfapUuXlJWVFd3tFjZt2tTq96Clurq6lJeXJ0lWrFiRwYMHl7giAAAAYHe1xt/5rTKzabNWyLGatUWYBQAAAMCuadU9m37zm9+ksbGx0K9nn322NUsGAAAAYA90uA3CzWgCAAAAaL86XNgEAAAAQPvVKns2Pfzww0mSIUOGFN73kCFDmvsHAAAAoH1plbDppJNOao1ukyR9+vRp1f4BAAAA2H2W0QEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQmHYbNtXX1+enP/1prrnmmlx11VW54447smbNmlKX1eq+//3v5+KLL87xxx+fnj17pqysLHPmzNnue55++ul89KMfzYABA9K3b9/8xV/8RW6//fa2KRgAAADgbbqVuoCtueWWW/KFL3whL7/8cvr06ZMuXbpkzZo1GTBgQK666qp89rOfLXWJrebLX/5yli9fngMPPDCHHHJIli9fvt32Dz/8cE477bT06tUr5513Xvr3758777wz48ePz4oVKzJt2rQ2qhwAAACgHc5suvrqqzNx4sScdNJJee6557JmzZqsXr06L7zwQs4+++x87nOfy5e+9KVSl9lqZs6cmWXLluWVV17ZYai2cePGTJ06NV26dMkvfvGL3HTTTbn++uvz61//OpWVlZk+ffoOwyoAAACAIrWrsOmxxx7Ll7/85fzjP/5jfvSjH6Wqqqr5WmVlZWbNmpVrr702//zP/5x77723hJW2ng996EM57LDDdqrtQw89lP/8z//MZz7zmYwYMaL59X333TfTp0/Phg0bMnfu3FaqFAAAAGBL7Sps+vrXv55Ro0blf/2v/7XNNp///Odz6qmn5vLLLy/03i+//HLuvffeXH755Tn99NNz4IEHpqysLGVlZZk0adIu9bV8+fJMmzYtw4YNS9++fbP//vtn5MiRue6667Ju3brCan7kkUeSJOPGjdvi2mmnnZYkefTRRwu7HwAAAMCOtJs9m9asWZNHHnkkN910U/NrGzduzO9+97skSXl5ebp3754kufjii/OpT30qK1euzKBBgwq5/8CBAwvp55577smECROyevXq5tfWrVuXhQsXZuHChZk5c2bmz5+foUOH7vG9amtrkyQVFRVbXDv44IPTr1+/5jYAAAAAbaHdzGxasWJFNm3alCOPPLL5teXLl2fo0KGprKzMokWLml8/6qij0tTU1Gr7ER166KFbnS20I4sWLcr48eOzevXq9OvXL1deeWWeeOKJPPjgg5k6dWqSpKamJmeccUbeeOONPa7zT3/6U5K3ls1tzT777NPcBgAAAKAttJuZTT169EiSbNiwofm17t2759BDD01ZWVl69uzZ/Pqbb76ZsrKy5vcU4fLLL8/IkSMzcuTIDBw4MMuWLcuQIUN2qY9LLrkk9fX16datW+6///6ccMIJzddOOeWUVFRU5LLLLktNTU2uv/76XHHFFYXVDwAAANAetJuZTYceemj69++fJ598ssVry5Yty9KlS/O+972v+fUnnngi3bp1S2VlZWH3/9rXvpYzzzxzt5fTPfXUU3nssceSJBdeeGGLoGmzadOmZfjw4UmSG2+8MQ0NDbtfcP57RtO2Zi+tXr16m7OeAAAAAFpDuwmbunfvnk984hP53//7f2fNmjXbbLd+/frccMMN+chHPpL+/fu3YYXbd/fddzcfT548eattunTpkgsuuCBJ8vrrr+fhhx/eo3tu3qtpa/syvfTSS1mzZs1W93MCAAAAaC3tJmxK3ppdtG7dunziE5/YauD05ptvZvz48Vm1alWuvfbaElS4bQsWLEiS9O3bN8cdd9w225100knNx48//vge3XNzX/fff/8W137+859vcT8AAACA1tauwqZDDz00d955Z5566qkcddRR+Zd/+Zf88pe/zMKFC/O9730vxxxzTB544IHcdtttGTZsWKnLbWHx4sVJkqFDh6Zbt21vhfX2uje/Z3edeuqpec973pNbb701zzzzTPPrf/rTn3LVVVelR48ezTOpAAAAANpCu9kgfLOxY8dm4cKFueyyy/L5z38+GzduTFNTU7p27ZoPf/jDufPOO3PUUUeVuswW1q9fn1dffTVJMnjw4O223W+//dK3b9+sXbs2K1as2OL6zJkzm2dJPffcc82vPfLII0mSMWPG5KKLLkqSdOvWLTNnzsxpp52WD37wgznvvPPSv3//3HnnnVm+fHm+8Y1v5PDDD9+lz1JXV7fd66tWrdql/gAAAIC9S7sLm5K3Zgf95Cc/yRtvvJHa2to0NjZm6NChGTBgQKlL26o33nij+bhfv347bL85bNraUsEFCxZk7ty5LV57/PHHWyy52xw2JW+FcwsWLMhXv/rV/OhHP0pDQ0OOPvroXHvttRk/fvwuf5by8vJdfg8AAADAZu0ybNqsf//+ef/731/qMnZo/fr1zcc9evTYYfuePXsmSerr67e4NmfOnMyZM2eX7j9q1Kjcd999u/QeAAAAgNbQrsOmjqJXr17Nxxs2bNhh+zfffDNJ0rt371araXdtbWnf261atSqjRo1qo2oAAACAjkbYVID+/fs3H29tadw7rV27NsnOLblrazvacwoAAABge9rV0+g6ql69euWAAw5IsuMNtl977bXmsMn+SAAAAEBnY2ZTQY488sg89thjWbJkSTZu3Jhu3bb+rX3hhReaj4cPH95W5e22qqqqFucNDQ0lqgQAAADoCMxsKsiYMWOSvLVE7le/+tU22z366KPNx6NHj271ugAAAADakrCpIGeffXbz8ezZs7faprGxMfPmzUuSDBgwIGPHjm2L0vZIdXV1i6+HHnqo1CUBAAAA7VirLaP7xS9+UXifH/zgBwvvsyijRo3KiSeemMceeyw333xzJk6cmBNOOKFFm+uvvz6LFy9OklxyySXp3r17KUoFAAAAaDWtFjadfPLJKSsrK6y/srKybNy4sbD+3mnBggVZsmRJ8/mrr77afLxkyZLMmTOnRftJkyZt0ceNN96Y0aNHp76+PuPGjcv06dMzduzY1NfX57bbbstNN92UJKmsrMy0adNa5XMAAAAAlFJZU1NTU2t03KVLl5SVlaWo7svKyrJp06ZC+tqaSZMmZe7cuTvdfluf65577smECROyevXqrV6vrKzM/PnzM3To0N2qs9Tq6uqan6K3YsWKDB48uMQVAQAAALurNf7Ob/Wn0fXu3Tsf+9jH8uEPfzhdunT+LaLOOuusPPvss7nxxhszf/781NXVpUePHhk6dGjOPffc/O3f/m369OlT6jJ3mqfRAQAAALui1WY27bvvvnnjjTfeuklZWQ4++OB85jOfyfnnn59jjjmmNW5JK9ha2FRbW5vEzCYAAADo6FpjZlOrhU3r16/PT3/609xyyy25//77s3HjxuY9nI4++uhccMEF+fSnP51DDjmkNW5PK7GMDgAAADqP1vg7v9XWtfXq1Svjx4/PvffemxdffDHf/OY3c+yxx6apqSnPPvtsvvCFL+TQQw/NRz7ykdx6662pr69vrVIAAAAAaCNtsonSQQcdlEsuuSQLFy5MdXV1/uEf/iGDBw/Opk2bcv/99+f888/PwIEDM2nSpDz44INtURIAAAAAraDNd+wePnx4rr766ixfvjwPPfRQJk2alH79+mXNmjWZN29exo0bl/Ly8vzjP/5jW5cGAAAAwB4q6ePhTj755MyaNSu///3vc+utt+b0009P165dm5fdAQAAANCxdCt1AclbT6vr0qVLysrKmjcRp33Y2tPoAAAAALalpGHTo48+mltuuSV33nlnVq9enSRpamrKIYcckvPPP7+UpQEAAACwG9o8bFq8eHFuueWW3HrrrVmxYkWStwKmPn365OMf/3guuOCCnHrqqenSpaQr/Piz6urqFudvfyQiAAAAwDu1Sdj08ssv54c//GFuueWWLFq0KMlbAVOXLl0yduzYXHDBBfnEJz6Rvn37tkU5AAAAALSSVgub1q9fn7vvvju33HJL/u3f/i2bNm1KU1NTkrf2Abrgggvy13/91xk0aFBrlQAAAABAG2u1sOld73pX1q5dm+StWUwHH3xwPv3pT+f888/PiBEjWuu2AAAAAJRQq4VNa9asSVlZWXr16pW/+qu/yrhx49K1a9c8++yzefbZZ3erzwsuuKDgKgEAAAAoUlnT5rVtBevSpUvKysoK66+srCwbN24srD92z9s3CF+xYkUGDx5c4ooAAACA3dUaf+e36gbhrZRj0YaqqqpanDc0NJSoEgAAAKAjaLWw6eGHH26trgEAAABop1otbDrppJNaq2vaUHV1dYvzt0+vAwAAAHinLqUuAAAAAIDOQ9gEAAAAQGGETQAAAAAUplX2bPrd736XJHn3u9+drl27Ftr3pk2b8uKLLyZJDj300EL7BgAAAGDPtMrMpsMPPzzvec978tvf/rbwvl944YXm/gEAAABoX1ptGV1TU1Nrdd0m/QMAAACw61p1z6aysrLW7B4AAACAdqZV9mzabNy4cenevXuhfTY0NBTaH9tXVVXV4tz3HwAAANieVgubmpqamjfyBgAAAGDv0Cph08SJE1ujW0qgurq6xXldXV3Ky8tLVA0AAADQ3rVK2DR79uzW6BYAAACAdq5VNwgHAAAAYO8ibAIAAACgMMImAAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMMImAAAAAArTrdQF0L5VVVW1OG9oaChRJQAAAEBHYGYTAAAAAIUxs4ntqq6ubnFeV1eX8vLyElUDAAAAtHdmNgEAAABQmJLNbHrllVfyX//1X3nppZeydu3adO/ePQMGDMihhx6aoUOHpmvXrqUqDQAAAIDd1GZh09q1a/PTn/409913Xx599NG8+OKL22zbs2fPHHvssRk3blw+/vGP55hjjmmrMgEAAADYA2VNTU1NrXmDRYsW5V/+5V/y4x//OOvWrUuS7Owty8rKkrz1RLTPfe5zOf/889OnT59Wq5Ude/ueTStWrMjgwYNLXBEAAACwu1rj7/xWC5sWLVqUr3zlK7nvvvuS/HfAdPDBB2fUqFE57rjj8q53vSv7779/9ttvv9TX1+ePf/xjXnvttdTU1OTpp5/Os88+m4aGhrcKLSvLAQcckMsuuyx/93d/l549e7ZG2eyAsAkAAAA6j9b4O79VltFNnjw5t9xySxobG5Mk73//+/PXf/3X+eQnP5lDDz10p/vZsGFDfvGLX+QHP/hB7rrrrrz66qv5h3/4h3z3u9/NvHnzMmbMmNYoHwAAAIDd1CpPo5s7d266deuWqVOn5oUXXsjChQtz6aWX7lLQlCQ9evTIhz70ocyePTu///3vM2/evBxxxBFZtmxZHnroodYoHQAAAIA90Cozm/7mb/4m//AP/9A8DasIPXv2zIQJE/LXf/3X+fGPf5xNmzYV1jcAAAAAxWiVsOnb3/52a3Sb5K29mz71qU+1Wv8AAAAA7L5WWUa3LaecckouvfTS/Pu//3tb3hYAAACANtIqM5u25dFHH82jjz6asWPHtuVtAQAAAGgjbTqz6aCDDkqSQvdyAgAAAKD9aNOZTcccc0wefPDBrFy5Mscee2xb3prdVFVV1eK8oaGhRJUAAAAAHUGbzmyaOHFimpqa8qMf/agtbwsAAABAGylrampqassbfvCDH8wTTzyRe+65J6effnpb3poC1NXVNS+DXLFiRQYPHlziigAAAIDd1Rp/57fpzKYf/ehH+frXv57hw4fnnHPOydy5c9vy9gAAAAC0sjYNmz796U/n1FNPzfPPP5/6+vpMmTIlf/EXf5G5c+dm1apVbVkKAAAAAK2gTTcIT5K3r9pramrK008/naeffjpJ8u53vzvHHntsRowYkWOPPTbHHntsDjvssLYuEQAAAIDd1KZh0+LFi/PrX/+6+euZZ57JypUrm6/X1dXlxRdfzL333tv82r777tscPl1//fVtWS4AAAAAu6jNNwh/pz/84Q8twqdf//rXWbx4cRoaGlq0Kysry6ZNm0pUJZvZIBwAAAA6j9b4O7/Nl9G90wEHHJBTTjklp5xySvNrDQ0NWbx4cXP4tPkLAAAAgPat5GHT1nTv3j3HHHNMjjnmmFKXAgAAAMAuaNOn0QEAAADQuQmbAAAAACiMsAkAAACAwgibAAAAACiMsAkAAACAwgibAAAAACiMsAkAAACAwnRrqxstWrQozz33XLp06ZIJEybs1Ht+/OMfp76+PhUVFTnhhBNauUIAAAAA9lSbhU1//OMfM2nSpJSVlWXw4ME5+eSTt9v+hRdeyPjx41NWVpZ58+YJmwAAAAA6gDZbRnfKKaekvLw8SXLLLbfssP3mNv37988nP/nJVq0NAAAAgGK0WdhUVlaWCRMmpKmpKXfccUfWr1+/3fY/+MEPUlZWlnPOOSe9evVqoyoBAAAA2BNtukH4pEmTkiRr1qzJXXfdtc12jz76aH73u9+1eA+lUVVV1eLrlFNOKXVJAAAAQDvWpmHT2zf6njt37jbbzZs3L0nynve8J2PGjGmT2gAAAADYc222QfhmEydOzJNPPpkHH3wwL730Ug4++OAW19evX5877rgjZWVlueCCC9q6PN6hurq6xXldXV3z3lsAAAAA79SmM5uS5LzzzkuvXr3S2NiYW2+9dYvrP/3pT/PGG28ImwAAAAA6oDYPm/bZZ5987GMfS1NTU/Nyubfb/BS6D37wgznssMPaujwAAAAA9kCbh03Jf2/6/dxzz+XXv/518+uvvPJK7r///pSVlWXixImlKA0AAACAPVCSsOnDH/5wBg0alOS/ZzIlya233pqNGzemb9++Offcc0tRGgAAAAB7oCRhU5cuXfLXf/3XaWpqyq233prGxsYkbwVPZWVl+eQnP5k+ffqUojQAAAAA9kBJwqYkmTx5cpLk97//ff7t3/4tixcvzn/8x38kiSV0AAAAAB1Ut1LdeNiwYRk5cmQWLlyYuXPn5vDDD0+SHHbYYTn55JNLVRYAAAAAe6BkYVPy1gymp59+Oj/72c8yYMCAlJWV5fzzzy9lSQAAAADsgZIto0uSz3zmM+nZs2fq6+uzcuXKJP/9pDoAAAAAOp6Shk0DBgzIWWedlaamppSVlWX06NEZMmRIKUsCAAAAYA+UNGxK/nsmU1NTk1lNAAAAAB1cSfdsSpKPfvSjaWxsLHUZAAAAABSg5DObAAAAAOg8hE0AAAAAFEbYBAAAAEBhhE0AAAAAFEbYBAAAAEBhhE0AAAAAFEbYBAAAAEBhhE0AAAAAFEbY1ME1Njbm29/+dt7//venT58+2WefffLBD34wP/vZz0pdGgAAALAXEjZ1YE1NTfnUpz6Vv/u7v8vq1atz4YUX5rzzzstvf/vbfOxjH8u3v/3tUpcIAAAA7GVaJWz6yU9+0hrdNlu5cmX+/d//vVXv0RHceeedufPOOzN69Og899xz+Zd/+ZfcdNNNqa6uzmGHHZbPf/7zWbZsWanLBAAAAPYirRI2nXPOORkxYkTuuOOOQvtdsWJF/uZv/ibvfe97c//99xfad0f005/+NEkyffr09O7du/n1Aw88MJdeemnefPPNzJ49u1TlAQAAAHuhVgmb3vve9+bZZ5/N+PHjM2TIkPzjP/5jqqurd6uvtWvX5vvf/34++tGP5r3vfW++973vZdOmTXnve99bcNW75uWXX869996byy+/PKeffnoOPPDAlJWVpaysLJMmTdqlvpYvX55p06Zl2LBh6du3b/bff/+MHDky1113XdatW7fN97300ktJkiFDhmxxbfNrDz300C7VAgAAALAnurVGp88//3y+9a1v5Z//+Z+zfPnyXHPNNbnmmmtSUVGRv/iLv8jIkSNz7LHH5l3velf222+/7Lfffqmvr88f//jHvPbaa6mpqcnTTz+dp556Kk899VTWr1+fpqamJMknPvGJXHXVVamsrGyN0nfawIEDC+nnnnvuyYQJE7J69erm19atW5eFCxdm4cKFmTlzZubPn5+hQ4du8d4DDzwwSbJ06dIMHz68xbWlS5cmSWpqagqpEwAAAGBnlDVtTnFawZo1a/Ld73433/nOd7JixYq3blhWttPv31xaz54984lPfCKXXHJJRo0a1Sq17qq3f45DDz00w4YNa17aN3HixMyZM2eHfSxatCijR49OfX19+vXrly996UsZO3Zs6uvrc9ttt2XGjBlJksrKyixcuDD9+/dv8f558+Zl4sSJOfHEE3P//fenV69eSZI//OEPOf7447Ns2bL06NEjb775ZkGfOqmrq0t5eXmSt5Y1Dh48uLC+AQAAgLbVGn/nt8rMps369euXyy67LJ///Ofzb//2b7n99tvz8MMP79Sm1b169coHPvCBfOxjH8sFF1yQ/fffvzVL3WWXX355Ro4cmZEjR2bgwIFZtmzZVpezbc8ll1yS+vr6dOvWLffff39OOOGE5munnHJKKioqctlll6WmpibXX399rrjiihbv/8xnPpM5c+bk4YcfztFHH52PfOQjaWhoyN13390886pLFw8cBAAAANpOq85s2pYXX3wxTzzxROrq6vLKK6/kj3/8Y3r16pWDDjooBx10UI4++ugcf/zx6d69e1uXttveHjbtzMymp556Kh/4wAeSJBdffHG+973vbdGmsbExRx11VBYvXpwBAwbk5Zdf3uJ78uabb+aaa67JrbfemmXLlmXffffNxz/+8Xz+859PZWVlDj300CxfvryYDxkzmwAAAKAz6XAzm7bl3e9+d84999xS3LrduPvuu5uPJ0+evNU2Xbp0yQUXXJAvfelLef311/Pwww9n3LhxLdr07NkzX/3qV/PVr361xeuPPPJIkuT4448vtG4AAACA7Wm1NVarVq1qra47hQULFiRJ+vbtm+OOO26b7U466aTm48cff3yn+//BD36QJDnvvPN2s0IAAACAXddqM5ve/e53513velfe9773ZcSIERkxYkSOPfbYHHHEEbu0SXhntXjx4iTJ0KFD063btv81DBs2bIv3vN3q1auzzz77tHjtjjvuyKxZszJy5Mh84hOfKKhiAAAAgB1r1WV0r7zySh544IE88MADza/17t07Rx11VIsA6phjjknv3r1bs5R2Zf369Xn11VeTZIdrIffbb7/07ds3a9eubX6i39t94AMfSHl5eYYPH55evXrlqaeeyiOPPJL3vOc9+fGPf5yuXbvuUm11dXXbvW7GGgAAALA9rb5n0zv3H1+3bl2efvrpPP30082vdenSJUOHDm0RQI0YMSLvete7Wru8knjjjTeaj/v167fD9pvDpjVr1mxxbfz48fnJT36Sf//3f09DQ0OGDBmSL3/5y/nCF76wxYynnbF5UzAAAACA3dHqYVPPnj3zV3/1VznjjDOydOnSPPPMM3nmmWdaPCFt06ZN+e1vf5uamprcfvvtza8PHDiwOYAaMWJEPvWpT7V2uW1i/fr1zcc9evTYYfuePXsmSerr67e4dsUVV+SKK64orDYAqKmpSUNDQ7p3757KyspSlwMAQAfTamHTPffck3/4h3/I888/nzvuuCMPPfRQvvKVr+SOO+5I165d86c//ak5eNr89fzzz6ehoaG5j5deeik///nP8/Of/zxlZWWdJmzq1atX8/GGDRt22P7NN99MkjZZari1pXpvt2rVqowaNarV6wCgdGpqalJfX5/evXsLmwAA2GWtFjadccYZOf300zNz5sxcccUVeemll3LppZfmO9/5Tq655pp8/OMfz0knndTiaWsNDQ15/vnnWwRQv/71r/P666+3Vpkl0b9//+bjrS2Ne6e1a9cm2bkld3tqR3tIAQAAAGxPl1btvEuX/I//8T9SW1ubr3zlK+nTp09qa2tzzjnn5MQTT8wvf/nLFu27d++e973vfZk4cWK++c1v5uGHH84f//jHLF26ND/5yU9as9Q21atXrxxwwAFJdrwh92uvvdYcNtlPCQAAAGjvWjVs2qxv37752te+lpqamkyZMiVdunTJ448/nr/8y7/M+PHj81//9V/bff9hhx2Wj33sY21Raps58sgjkyRLlizJxo0bt9nuhRdeaD4ePnx4q9cFAAAAsCfaJGza7JBDDsnMmTOzaNGijBs3Lk1NTbnjjjty5JFHZtq0aXnttdfaspySGjNmTJK3lsj96le/2ma7Rx99tPl49OjRrV7XO1VVVbX4OuWUU9q8BgAAAKDjaNOwabOjjjoq/+///b/8/Oc/z9FHH50NGzbkW9/6VoYOHZobbrihxSbhndXZZ5/dfDx79uyttmlsbMy8efOSJAMGDMjYsWPbojQAAACA3VaSsGmzD3/4w1m0aFFmzZqVQw45JK+99lq+8IUv5Igjjsjtt99eytJa3ahRo3LiiScmSW6++eY8+eSTW7S5/vrrs3jx4iTJJZdcku7du7dpjUlSXV3d4uuhhx5q8xoAAACAjqOsqampqdRFJMn69evzta99Ldddd10aGxvTu3fv5o2x26MFCxZkyZIlzeevvvpqvvCFLyR5a7nbRRdd1KL9pEmTtuhj0aJFGT16dOrr69OvX79Mnz49Y8eOTX19fW677bbcdNNNSZLKysosXLiwxVPsSqWurq55o/IVK1Z4eh1AJ1NbW5vp06dn5cqVGTRoUK666qpUVFSUuiwAAFpJa/ydX5Kwad26dXn++edTXV3d/M/q6uqsWLEiTU1NaWpqSq9evbJu3bq2Lm2nTZo0KXPnzt3p9tv6Nt9zzz2ZMGFCVq9evdXrlZWVmT9/foYOHbpbdRZN2ATQec2ePTtTp07Npk2bml/r2rVrZsyYkcmTJ5ewMgAAWktr/J3fbY972I4dhUqbvf24e/fuOeKII/L+97+/NUtrN84666w8++yzufHGGzN//vzU1dWlR48eGTp0aM4999z87d/+bfr06VPqMgHo5Gpra7cImpJk06ZNmTp1asaMGWOGEwAAO6XVZjYNGTJku6FS8tbT6Y455pgcc8wxOfroo3PMMcdk+PDhJdmbiK2rqqpqcd7Q0JDa2tokZjYBdCZf+tKXcs0112zz+he/+MVcffXVbVgRAABtoUPNbFq+fHnzcZ8+fVJVVdUcKG3+2n///Vvr9gDALli6dOkeXQcAgM1adRldWVlZ+vTpk3HjxuXYY4/NiBEjMmLECLNhOpDq6uoW529PPAHoPIYMGbJH1wEAYLNWW0bXpUuX/75JWVmLa/vtt1/e9773ZcSIEc3/PPLII9OtW6tmXxTABuEAnVNtbW2GDx++xZ5NyVubhC9evNieTQAAnVCHWkb3f/7P/8kzzzyTZ555Js8991yLJ8v98Y9/zCOPPJJHHnmk+bXu3btn+PDhzbOfNodQAwYMaK0SAYA/q6ioyIwZM7b5NDpBEwAAO6vVZja9XVNTU2pqaprDp81fv//977cs6B2zoMrLy/O+970vxx57bK644orWLpUdMLMJoHOrra3N9OnTs3LlygwaNChXXXWVoAkAoBNrjb/z2yRs2pbf//73WwRQtbW1aWxs3KJtWVnZVqf207aETQCd37333pv6+vr07t07Z555ZqnLAQCgFXWoZXQ7Y+DAgTnttNNy2mmnNb9WX1+fZ599tkUA9dxzz6W+vr6Ele69qqqqWpw3NDSUqBIAAACgI2h3O3L37t07H/jAB/KBD3yg+bXNy/AAAAAAaN/aXdi0NWVlZTniiCNKXcZeqbq6usX526fXAQAAALxTl1IXAAAAAEDnIWwCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAK0yGeRkfpVFVVtThvaGgoUSUAAABARyBsAgBaqKysTENDQ7p3717qUgAA6ICETWxXdXV1i/O6urqUl5eXqBoA2kJlZWWpSwAAoAOzZxMAAAAAhRE2AQAAAFAYYRMAAAAAhRE2AQAAAFAYYRMAAAAAhRE2AQAAAFCYbqUugPatqqqqxXlDQ0OJKgEAAAA6AjObAAAAACiMmU1sV3V1dYvzurq6lJeXl6gaAAAAoL0zswkAAACAwgibAAAAACiMsAkAAACAwgibAAAAACiMsAkAAACAwgibAAAAACiMsAkAAACAwgibAAAAACiMsAkAAACAwnQrdQG0b1VVVS3OGxoaSlQJAAAA0BGY2QQAAABAYcxsYruqq6tbnNfV1aW8vLxE1QAAAADtnZlNAAAAABRG2AQAAABAYYRNAAAAABRG2AQAAABAYYRNAAAAABRG2AQAAABAYYRNAAAAABRG2AQAAABAYYRNAAAAABRG2AQAAABAYYRNAAAAABRG2AQAAABAYbqVugDat6qqqhbnDQ0NJaoEAAAA6AjMbAIAAACgMGY2sV3V1dUtzuvq6lJeXl6iagAAAID2zswmAAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMMImAAAAAAojbAIAAACgMN1KXQDtW1VVVYvzhoaGElUCAAAAdARmNgEAAABQGDOb2K7q6uoW53V1dSkvLy9RNQAAAEB7Z2YTAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGFTB9fU1JSf/OQnGTt2bA455JD06dMnRxxxRC6++OL813/9V6nLAwAAAPYywqYO7vOf/3w++clP5re//W3OPvvs/N3f/V2GDBmSGTNmZMSIEfnNb35T6hIBAACAvUi3UhfA7nvppZfyrW99K4cddlh+/etfZ999922+9s1vfjN///d/nxtuuCGzZs0qYZUAAADA3sTMpg5s2bJlaWxszOjRo1sETUly5plnJkleeeWVUpQGAAAA7KWETbvp5Zdfzr333pvLL788p59+eg488MCUlZWlrKwskyZN2qW+li9fnmnTpmXYsGHp27dv9t9//4wcOTLXXXdd1q1bt833VVRUpEePHnn88cezevXqFtfuvffeJMmpp566y58NAAAAYHdZRrebBg4cWEg/99xzTyZMmNAiLFq3bl0WLlyYhQsXZubMmZk/f36GDh26xXsPOOCAXHPNNc1B1cc+9rHss88++fWvf52HHnoof/M3f5O//du/LaROAAAAgJ0hbCrAoYcemmHDhuX+++/fpfctWrQo48ePT319ffr165cvfelLGTt2bOrr63PbbbdlxowZqampyRlnnJGFCxemf//+W/Rx6aWX5t3vfncuuuiifO9732t+fcyYMfnMZz6Tbt38KwYAAADajmV0u+nyyy/PPffck5deeinLly/P//2//3eX+7jkkktSX1+fbt265f7778/06dNzwgkn5JRTTslNN92Uf/7nf06S1NTU5Prrr99qH1//+tczYcKETJ8+PStWrMgbb7yRxx57LOvXr8/JJ5+cn/3sZ3v0OQEAAAB2RVlTU1NTqYvoDJYtW5YhQ4YkSSZOnJg5c+Zst/1TTz2VD3zgA0mSiy++uMWspM0aGxtz1FFHZfHixRkwYEBefvnldO/evfn6Aw88kA9/+MO59NJLc8MNN7R470svvZT3vOc9efe7353a2to9/HT/ra6uLuXl5UmSFStWZPDgwYX1DQAAALSt1vg738ymErn77rubjydPnrzVNl26dMkFF1yQJHn99dfz8MMPt7h+3333JUnGjh27xXsPPvjgDBs2LEuWLMmaNWsKqhoAAABg+4RNJbJgwYIkSd++fXPcccdts91JJ53UfPz444+3uLZhw4YkySuvvLLV977yyivp0qVLi9lQAAAAAK1J2FQiixcvTpIMHTp0u5t4Dxs2bIv3bDZ69OgkyQ033JA//elPLa5973vfS11dXU444YT07NmzqLIBAAAAtsujykpg/fr1efXVV5Nkh2sh99tvv/Tt2zdr167NihUrWlw799xz83/+z//JL37xi1RWVuav/uqvMmDAgPzHf/xHHnroofTu3XuLvZx2pK6ubrvXV61atUv9AQAAAHsXYVMJvPHGG83H/fr122H7zWHTO/de6tq1a+6///5885vfzO23355bb701GzZsyMCBA5ufUDd8+PBdqm3zpmAAAAAAu0PYVALr169vPu7Ro8cO229eBldfX7/Va1/84hfzxS9+sbgCAQAAAHaTsKkEevXq1Xy8eZPv7XnzzTeTJL179261mjZ751K9d1q1alVGjRrV6nUAAAAAHZOwqQT69+/ffPzOpXFbs3bt2iQ7t+RuT+1oDykAAACA7fE0uhLo1atXDjjggCQ73pD7tddeaw6b7KcEAAAAtHfCphI58sgjkyRLlizJxo0bt9nuhRdeaD7e1c2+AQAAANqaZXQlMmbMmDz22GNZu3ZtfvWrX+UDH/jAVts9+uijzcejR49uq/KaVVVVtThvaGho8xoAAACAjsPMphI5++yzm49nz5691TaNjY2ZN29ekmTAgAEZO3ZsW5QGAAAAsNuETSUyatSonHjiiUmSm2++OU8++eQWba6//vosXrw4SXLJJZeke/fubVpjklRXV7f4euihh9q8BgAAAKDjsIxuNy1YsCBLlixpPn/11Vebj5csWZI5c+a0aD9p0qQt+rjxxhszevTo1NfXZ9y4cZk+fXrGjh2b+vr63HbbbbnpppuSJJWVlZk2bVqrfA4AAACAIpU1NTU1lbqIjmjSpEmZO3fuTrff1rf5nnvuyYQJE7J69eqtXq+srMz8+fMzdOjQ3aqzaHV1dc1PxVuxYkUGDx5c4ooAAACA3dUaf+dbRldiZ511Vp599tlceumlqaysTJ8+fTJgwIAcf/zxufbaa7No0aJ2EzQBAAAA7IiZTWzX1p5GV1tbm8TMJgAAAOjozGwCAAAAoF2zQTjbVV1d3eL87YknAAAAwDuZ2QQAAABAYYRNAAAAABRG2AQAAABAYYRNAAAAABTGBuFsV1VVVYvzhoaGElUCAAAAdARmNgEAAABQGDOb2K7q6uoW53V1dSkvLy9RNQAAAEB7Z2YTAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIXxNDq2q6qqqsV5Q0NDiSoBAAAAOgIzmwAAAAAojJlNbFd1dXWL87q6upSXl5eoGgAAAKC9M7MJAAAAgMIImwAAAAAojLAJAAAAgMIImwAAAAAojLAJAAAAgMIImwAAAAAoTLdSF0D7VlVV1eK8oaGhRJUAAAAAHYGZTQAAAAAUxswmtqu6urrFeV1dXcrLy0tUDQAAANDemdkEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGG6lboA2reqqqoW5w0NDSWqBAAAAOgIhE0ASWpqatLQ0JDu3bunsrKy1OUAAAB0WMImtqu6urrFeV1dXcrLy0tUDbSempqa1NfXp3fv3sImAACAPWDPJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAKI2wCAAAAoDDCJgAAAAAK063UBQCUWm1tbebOnZuVK1dm0KBBOeKII1JRUVHqsgAAADokYROwV5s9e3amTp2aTZs2Nb921113ZcaMGZk8eXIJKwMAAOiYLKMD9lq1tbVbBE1JsmnTpkydOjW1tbUlqgwAAKDjMrOJ7aqqqmpx3tDQUKJKoHizZs3aImjabNOmTZk1a1auvvrqNq4KAACgYzOzCdhrLV26dI+uAwAAsCUzm9iu6urqFud1dXUpLy8vUTVQrCFDhuzRdQAAALZkZhOw15oyZUq6du261Wtdu3bNlClT2rgiAACAjk/YBOy1KioqMmPGjC0Cp65du2bGjBmpqKgoUWUAAAAdl2V0wF5t8uTJGTNmTKZPn56VK1dm0KBBueqqqwRNAAAAu0nYBOz1KioqMnHixNTX16d3796CJgAAgD1gGR0AAAAAhRE2AQAAAFAYYRMAAAAAhRE2AQAAAFAYYRMAAAAAhRE2AQAAAFAYYRMAAAAAhRE2AQAAAFAYYRMAAAAAhRE2AQAAAFAYYRMAAAAAhRE2AQAAAFAYYRMAAAAAhRE2AQAAAFCYbqUugPatqqqqxXlDQ0OJKgEAAAA6AmETQJLKyso0NDSke/fupS4FAACgQxM2sV3V1dUtzuvq6lJeXl6iaqD1VFZWlroEAACATsGeTQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETQAAAAAURtgEAAAAQGGETR3YnDlzUlZWtt2vU089tdRlAgAAAHuRbqUugN03YsSIfPWrX93qtTvuuCPV1dU57bTT2rgqAAAAYG8mbOrARowYkREjRmzx+oYNG/Ltb3873bp1y8SJE9u+MAAAAGCvZRldJ3T33XfnD3/4Q84888wMHDiw1OUAAAAAexFh0256+eWXc++99+byyy/P6aefngMPPLB5n6RJkybtUl/Lly/PtGnTMmzYsPTt2zf7779/Ro4cmeuuuy7r1q3b5dpmzpyZJLnooot2+b0AAAAAe8Iyut1U1Iyhe+65JxMmTMjq1aubX1u3bl0WLlyYhQsXZubMmZk/f36GDh26U/0tX748Dz74YAYPHpyPfOQjhdQIAAAAsLPMbCrAoYcemnHjxu3y+xYtWpTx48dn9erV6devX6688so88cQTefDBBzN16tQkSU1NTc4444y88cYbO9Xn7Nmz09jYmEmTJqVr1667XBMAAADAnjCzaTddfvnlGTlyZEaOHJmBAwdm2bJlGTJkyC71cckll6S+vj7dunXL/fffnxNOOKH52imnnJKKiopcdtllqampyfXXX58rrrhiu/01NjZm9uzZKSsry5QpU3bnYwEAAADsETObdtPXvva1PdqA+6mnnspjjz2WJLnwwgtbBE2bTZs2LcOHD0+S3HjjjWloaNhunw888EB+97vf5ZRTTtnl4AsAAACgCMKmErn77rubjydPnrzVNl26dMkFF1yQJHn99dfz8MMPb7dPG4MDAAAApSZsKpEFCxYkSfr27Zvjjjtum+1OOumk5uPHH398m+3+8Ic/5Kc//Wn233//fPzjHy+uUAAAAIBdIGwqkcWLFydJhg4dmm7dtr111rBhw7Z4z9bccsst2bBhQyZMmJCePXsWVygAAADALrBBeAmsX78+r776apJk8ODB22273377pW/fvlm7dm1WrFixzXY333xzkj1fQldXV7fd66tWrdqj/gEAAIDOTdhUAm+88Ubzcb9+/XbYfnPYtGbNmq1ef+qpp/Kb3/wmo0aNytFHH71HtZWXl+/R+wEAAIC9m7CpBNavX9983KNHjx2237wsrr6+fqvXR40alaampmKKAwAAANgDwqYS6NWrV/Pxhg0bdtj+zTffTJL07t271WrabHtL9ZK3ltGNGjWq1esAAAAAOiZhUwn079+/+XhbS+Pebu3atUl2bsndntrRHlIAAAAA2+NpdCXQq1evHHDAAUl2vCH3a6+91hw22U8JAAAAaO/MbCqRI488Mo899liWLFmSjRs3plu3rf+reOGFF5qPhw8f3lblNauqqmpx3tDQ0OY1AAAAAB2HmU0lMmbMmCRvLZH71a9+tc12jz76aPPx6NGjW70uAAAAgD0hbCqRs88+u/l49uzZW23T2NiYefPmJUkGDBiQsWPHtkVpLVRXV7f4euihh9q8BgAAAKDjEDaVyKhRo3LiiScmSW6++eY8+eSTW7S5/vrrs3jx4iTJJZdcku7du7dpjQAAAAC7yp5Nu2nBggVZsmRJ8/mrr77afLxkyZLMmTOnRftJkyZt0ceNN96Y0aNHp76+PuPGjcv06dMzduzY1NfX57bbbstNN92UJKmsrMy0adNa5XMAAAAAFKmsqampqdRFdESTJk3K3Llzd7r9tr7N99xzTyZMmJDVq1dv9XplZWXmz5+foUOH7ladRaurq2t+Kt6KFSsyePDgElcEAAAA7K7W+DvfMroSO+uss/Lss8/m0ksvTWVlZfr06ZMBAwbk+OOPz7XXXptFixa1m6AJAAAAYEfMbGK7qqqqWpw3NDSktrY2iZlNAAAA0NGZ2QQAAABAu2aDcLarurq6xfnbE08AAACAdzKzCQAAAIDCCJsAAAAAKIywCQAAAIDC2LOJXbJx48bm41WrVpWwEgAAAGBPvf1v+7f/zb8nhE1sV1VVVYvzdevWNR+PGjWqrcsBAAAAWskrr7ySww8/fI/7sYyOXbJp06ZSlwAAAAC0Y2VNTU1NpS6CjmPJkiWpqKhIkjzxxBMpLy8vcUXFOuWUU5IkDz30UKe7fxF970kfu/renW1fRLtVq1Y1z9R76qmncsghh+xUjR1FZx3XRfW7u/201pje2bY7atOZx3VnHdNF9t0Zx3VnHtOJcd2afRjXpWNct14ffrcunc44rjdu3JhXXnklSXL00UenV69ee9ynZXTskrcPuvLy8gwePLiE1RSve/fuSVKyz9Wa9y+i7z3pY1ffu7Pti253yCGHGNcd5P5F9bu7/bTWmN7ZtrvSX2cb1511TBfZd2cf151tTCfGdWv2YVyXjnHden343bp0Ouu4LmLp3NtZRgcAAABAYYRNAAAAABRG2AQAAABAYWwQzi6pq6tr3hR8xYoVnW79LXsn45rOyLimszGm6YyMazoj45rEzCYAAAAACiRsAgAAAKAwwiYAAAAACmPPJgAAAAAKY2YTAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGETAAAAAIURNgEAAABQGGET7cb69evz93//9/ngBz+YQYMGpVevXjn44IMzevTozJ49Ow0NDaUuEXbZiy++mG9961sZN25cDj300PTo0SMHH3xwPvnJT+aXv/xlqcuD3fL9738/F198cY4//vj07NkzZWVlmTNnTqnLgh16+umn89GPfjQDBgxI37598xd/8Re5/fbbS10W7DY/j+ls/O7ceZQ1NTU1lboISJJXX3015eXlGTVqVCorK3PQQQfltddey3333Zfly5dn3Lhxue+++9Kli4yUjuOLX/xirr322rz3ve/NySefnIMOOii1tbW5++6709TUlFtvvTXjx48vdZmwSw4//PAsX748Bx54YPr27Zvly5dn9uzZmTRpUqlLg216+OGHc9ppp6VXr14577zz0r9//9x5551Zvnx5vvGNb2TatGmlLhF2mZ/HdDZ+d+48hE20G42Njdm4cWN69OjR4vWNGzfmwx/+cB555JHce++9OeOMM0pUIey6n/zkJznggANy0kkntXj9sccey6mnnpp+/fpl1apV6dmzZ4kqhF33wAMPpKKiIocddliuueaafOlLX/LHDe3axo0bM2zYsNTV1eXf//3fM2LEiCTJn/70p4waNSrLli1LTU1NDjvssNIWCrvIz2M6G787dx6miNBudOnSZYugKUm6deuWj3/840mSJUuWtHVZsEc+8YlPbPEfyyQ58cQTM3bs2Lz22mt57rnnSlAZ7L4PfehD/iinQ3nooYfyn//5n/nMZz7THDQlyb777pvp06dnw4YNmTt3bukKhN3k5zGdjd+dOw9hUyfx8ssv5957783ll1+e008/PQceeGDKyspSVla2y/9nY/ny5Zk2bVqGDRuWvn37Zv/998/IkSNz3XXXZd26da3zAbajsbEx/+///b8kyVFHHdXm96d0OvO4TpLu3bsneStQZe/Q2cc0bE17GPePPPJIkmTcuHFbXDvttNOSJI8++ugu1cLerT2Mayhaex/XfnfuWPxb6iQGDhxYSD/33HNPJkyYkNWrVze/tm7duixcuDALFy7MzJkzM3/+/AwdOrSQ+23Nhg0bctVVV6WpqSl/+MMf8uCDD+aFF17I5MmTc+qpp7bafWl/OtO4fqff/e53eeCBB3LIIYfk6KOPbrP7UlqdeUzDtrSHcV9bW5skqaio2OLawQcfnH79+jW3gZ3RHsY1FK09j2u/O3c8ZjZ1QoceeuhW/8/djixatCjjx4/P6tWr069fv1x55ZV54okn8uCDD2bq1KlJkpqampxxxhl54403ii672YYNG/K1r30tX//61/Od73wnv/3tb/P5z38+N910U6vdk/avo4/rt2toaMj555+fN998M9dee226du3aJvelfelMYxp2VqnG/Z/+9Kckby2b25p99tmnuQ3sKj/P6Yza07j2u3PHZGZTJ3H55Zdn5MiRGTlyZAYOHJhly5ZlyJAhu9THJZdckvr6+nTr1i33339/TjjhhOZrp5xySioqKnLZZZelpqYm119/fa644oot+pg2bVrefPPNXbrnO/8vY79+/dLU1JTGxsasXLky99xzT6ZPn54nn3wy//qv/5p99tlnlz4XHVdnGtebNTY2ZtKkSfnFL36RqVOn5vzzz9+lz0PH1hnHNOxIexn3UCTjms6oPY5rvzt3YE10SkuXLm1K0pSkaeLEiTts/8tf/rK5/cUXX7zVNps2bWoaPnx4U5KmAQMGNG3YsGGLNn379m3uZ2e+Hn744Z36PLfffntTkqbLLrtsp9rTOXX0cb1p06amiRMnNiVpmjBhQtOmTZt25ePTCXX0Md3U1NR09dVXNyVpmj179k5+avZ2pRj355xzTlOSpoULF271/f369WsqLy/f5c8Cm5Xq5/nb+XlM0Uo9rv3u3LFZRkeS5O67724+njx58lbbdOnSJRdccEGS5PXXX8/DDz+8RZs1a9akqalpp79OPvnknapv8xTOzRt8ws5oT+O6sbExkydPzty5c/PpT386c+bMSZcufgSza9rTmIa2UsS43zwzb2v7Mr300ktZs2aN2Xu0qaJ+nkN7UuS49rtzx+ffFkmSBQsWJEn69u2b4447bpvt3v4Yyscff7zV69ps5cqVSf77CQSwM9rLuN78H8t58+Zl/PjxueWWW6w1Z7e0lzENbamIcb/52v3337/F+37+859v8X5obX6e0xkVNa797tw5CJtIkixevDhJMnTo0O0+SnLYsGFbvKcozz///FYfg7lu3br8/d//fZLkox/9aKH3pHNrD+O6sbExU6ZMybx583Luuefm+9//vv9Ystvaw5iGtlbEuD/11FPznve8J7feemueeeaZ5tf/9Kc/5aqrrkqPHj2a/087tAU/z+mMihjXfnfuPGwQTtavX59XX301STJ48ODttt1vv/3St2/frF27NitWrCi0jttvvz033HBDxowZk8MPPzz77LNPXnzxxdx33335wx/+kBNPPDGXXnppofek82ov4/rrX/965s6dm379+qWysjL/9E//tEWbs88+OyNGjCj0vnQ+7WVMJ8nMmTOb/+/lc8891/za5qXOY8aMyUUXXVT4fdn7FDXuu3XrlpkzZ+a0007LBz/4wZx33nnp379/7rzzzixfvjzf+MY3cvjhh7fWx4AWivx57ucx7UVR49rvzp2HsIkWj5zs16/fDttv/sGwZs2aQus488wzs3LlyjzxxBN58skns2bNmuy777455phjct5552XKlCnbTcjh7drLuF62bFmSt/bIufLKK7fa5vDDD/cfTHaovYzp5K1p8nPnzm3x2uOPP95iKrw/bihCkeN+7NixWbBgQb761a/mRz/6URoaGnL00Ufn2muvzfjx4wutG7anyHHt5zHtRVHj2u/OnYe/3Mn69eubj3v06LHD9j179kyS1NfXF1rH8ccfn+OPP77QPtl7tZdxPWfOnMyZM6fQPtk7tZcxnRjXtJ2ix/2oUaNy3333FVMc7KYix7Wfx7QXRY1rY7rzsGcT6dWrV/Pxhg0bdtj+zTffTJL07t271WqCPWVc09kY0+yNjHs6I+Oazsi45p2ETaR///7Nxzuz3GLt2rVJdm56JJSKcU1nY0yzNzLu6YyMazoj45p3EjaRXr165YADDkiS1NXVbbfta6+91vyDoby8vNVrg91lXNPZGNPsjYx7OiPjms7IuOadhE0kSY488sgkyZIlS7Jx48ZttnvhhReaj4cPH97qdcGeMK7pbIxp9kbGPZ2RcU1nZFzzdsImkrz1WNTkremMv/rVr7bZ7tFHH20+Hj16dKvXBXvCuKazMabZGxn3dEbGNZ2Rcc3bCZtIkpx99tnNx7Nnz95qm8bGxsybNy9JMmDAgIwdO7YtSoPdZlzT2RjT7I2Mezoj45rOyLjm7YRNJHnrUcAnnnhikuTmm2/Ok08+uUWb66+/PosXL06SXHLJJenevXub1gi7yrimszGm2RsZ93RGxjWdkXHN25U1NTU1lboI9tyCBQuyZMmS5vNXX301X/jCF5K8NTXxoosuatF+0qRJW/SxaNGijB49OvX19enXr1+mT5+esWPHpr6+PrfddltuuummJEllZWUWLlzY4okD0BqMazobY5q9kXFPZ2Rc0xkZ1xSqiU5h4sSJTUl2+mtbfvaznzXts88+23xfZWVlU21tbRt+MvZmxjWdjTHN3si4pzMyrumMjGuKZBkdLZx11ll59tlnc+mll6aysjJ9+vTJgAEDcvzxx+faa6/NokWLMnTo0FKXCbvEuKazMabZGxn3dEbGNZ2RcU1iGR0AAAAABTKzCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAAAAAKIywCQAAAIDCCJsAADqBSZMmpaysbIuvZcuWbdH2iiuuaL7eET3yyCNb/axXXHFFqUsDACJsAgAAAKBA3UpdAAAAxRk0aFB+/vOfN5+/+93vLmE1rWPkyJF57rnnms+PPvroElYDALyTsAkAoBPp3r17jjrqqFKX0ar69u3b6T8jAHRkltEBAAAAUBhhEwAAAACFETYBAJTICy+80PwktR/+8IdpamrKD37wg4wbNy4DBw5Mly5dcvLJJ5ektmeeeSYDBw5MWVlZDjnkkDz77LPN1975NLvXX389X/3qV1NVVZV+/fpl//33z9ixY/PDH/5wp+71+OOP56KLLsoRRxyRffbZJz169MjgwYNz5pln5jvf+U5ef/311viIAEArsWcTAECJ/PrXv24+Pvjgg3PSSSflsccea9HmmGOOaeuy8thjj+Wss87Kn/70pxx++OF54IEH8t73vnerbZcuXZoPf/jD+c///M/m19auXZtHHnkkjzzySO6+++784Ac/SLduW/7aWV9fnwsvvHCrodSLL76YF198MfPnz88rr7ySK664orDPBwC0LmETAECJvD1s+ru/+7s8//zz+cxnPpPzzjsvgwYNyooVK/Kud72rTWuaP39+zj333NTX16eqqir3339/Bg0atM3248ePz9KlS/PZz34255xzTvbdd988++yzufbaa1NTU5Pbb789gwYNyje/+c0W72tsbMzHPvax/Nu//VuSpKKiIn/zN3+T448/Pn369MmqVavyxBNP5Pbbb2/VzwsAFE/YBABQIm8Pm/7zP/8zP/vZz3LmmWc2v3bccce1aT233nprJk6cmI0bN2bUqFG57777sv/++2/3PU8//XRuvfXWfPrTn25+7fjjj8+5556bE088Mb/+9a/zv//3/86FF17Y4gly3/72t5uDpo9//OP54Q9/mJ49e7bo+4wzzsj/+l//K6tWrSrwUwIArc2eTQAAJfLMM880H3/nO99pETS1te9+97uZMGFCNm7cmFNPPTUPPvjgDoOmJDnzzDNbBE2b9e/fPzfddFOSt2Yxfe9732u+1tjYmOuuuy5JMnjw4MybN2+LoGmzLl265N3vfvfufCQAoESETQAAJfDqq69m5cqVSZKRI0dmypQpJavln/7pn/K5z30uTU1N+fjHP5758+enX79+O/XeyZMnb/PaqFGjUlVVlSR54IEHml9/5plnUldXlySZOnXqTt8LAOgYhE0AACXw9iV0n/3sZ0tWx6WXXpqvfOUrSd4Kjn784x9vc5bR1owcOXK710eNGpUkqampyYYNG5IkixYtar5+4okn7mrJAEA7J2wCACiBt4dNp512Wsnq+Na3vpUkOeqoozJz5sx07dp1l96/ow3MBw4cmCRpamrKa6+9luStWV2bHXLIIbt0PwCg/RM2AQCUwOb9mt797neXdE+iT37yk0mS3/zmN7nkkkt2+f1lZWVFlwQAdHDCJgCAEtg8s+nYY48taR0//OEPc/bZZyd56wlxl1566S69//e///1OXS8rK8t+++2XJDnwwAObr3vSHAB0PsImAIA2tmHDhixevDhJMmLEiJLW0r179/zoRz9qfhLet771rXzhC1/Y6fc//fTTO3W9oqIiPXr0SJK8//3vb77+i1/8YldLBgDaOWETAEAbW7x4cRoaGpKUfmZTkvTo0SN33nlnPvrRjyZJvvGNb+SLX/ziTr137ty527z29NNP5ze/+U2S5EMf+lDz6+973/tSXl6eJJk5c2bWrFmzu6UDAO2QsAkAoI1t3q8paR9hU/JW4PSTn/ykebPya6+9Nl/+8pd3+L6f/exnuf3227d4fc2aNbn44ouTJF26dGk+3ny+efZUXV1dLrjgguYn1b1TY2NjVq5cucufBwAoHWETAEAb27xf07777pshQ4aUuJr/1rNnz9x999358Ic/nCS58sor89WvfnW77zn++OPzmc98Jp/73Ofy8MMP51e/+lVmz56d448/PosWLUqSfO5zn8sxxxzT4n2f+9znmu9z11135eijj86NN96Yxx9/PIsWLcp9992Xr371qxk2bFhuuummVvi0AEBr6VbqAgAA9jabw6ZS79e0Nb169cpPf/rTnHnmmXnooYfy9a9/Pd27d9/mLKfbb789p556ar773e/mu9/97hbXP/nJT+aGG27Y4vUuXbrk7rvvzsSJE3PHHXekpqYm/9//9/8V/XEAgBIwswkAoI2157ApSXr37p177rknJ510UpLkK1/5Sq6++uqtth0yZEh+9atfZfr06Rk+fHj69OmTfffdNx/84Afz/e9/P3fccUe6ddv6/9/s06dPfvzjH+ehhx7K+eefnyFDhqR3797p0aNHysvLc9ZZZ+X//t//m2nTprXaZwUAimdmEwBAG3v11VdLev8rrrgiV1xxxXbb9OnTJ4888shO9bfffvvlyiuvzJVXXrlb9YwdOzZjx47drfcCAO2PsAkAoBNpaGhofgJckhxxxBHp3r17CSsq3tq1a7N06dJSlwEAbIOwCQCgE1m5cmWOPvro5vOlS5fm8MMPL11BreDpp582EwoA2jF7NgEAAABQGDObAAA6gTlz5mTOnDmlLqNNnHzyyWlqaip1GQDANpjZBAAAAEBhypr8byEAAAAACmJmEwAAAACFETYBAAAAUBhhEwAAAACFETYBAAAAUBhhEwAAAACFETYBAAAAUBhhEwAAAACFETYBAAAAUBhhEwAAAACFETYBAAAAUBhhEwAAAACFETYBAAAAUBhhEwAAAACFETYBAAAAUBhhEwAAAACFETYBAAAAUBhhEwAAAACFETYBAAAAUBhhEwAAAACFETYBAAAAUBhhEwAAAACF+f8BrZFr8CSRtQkAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 390, "width": 589 } }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(6, 4), layout=\"tight\")\n", "\n", "ax.errorbar(\n", " mwdata1[\"r\"],\n", " mwdata1[\"Menc\"],\n", " yerr=(mwdata1[\"Menc_err_neg\"], mwdata1[\"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(r\"$M(" ] }, "metadata": { "image/png": { "height": 390, "width": 589 } }, "output_type": "display_data" } ], "source": [ "r = np.logspace(-3, 3, 256) * u.kpc\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(6, 4), layout=\"tight\")\n", "\n", "ax.errorbar(\n", " mwdata1[\"r\"],\n", " mwdata1[\"Menc\"],\n", " yerr=(mwdata1[\"Menc_err_neg\"], mwdata1[\"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", "# Use symmetry coordinates for spherical mass_enclosed calculation\n", "fit_menc = init_potential.mass_enclosed(R=r)\n", "ax.loglog(r.value, 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(r\"$M(" ] }, "metadata": { "image/png": { "height": 390, "width": 589 } }, "output_type": "display_data" } ], "source": [ "r = np.logspace(-3, 3, 256) * u.kpc\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(6, 4), layout=\"tight\")\n", "\n", "ax.errorbar(\n", " mwdata1[\"r\"],\n", " mwdata1[\"Menc\"],\n", " yerr=(mwdata1[\"Menc_err_neg\"], mwdata1[\"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(R=r)\n", "ax.loglog(r.value, 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(r\"$M(" ] }, "metadata": { "image/png": { "height": 390, "width": 589 } }, "output_type": "display_data" } ], "source": [ "r = np.logspace(-3, 3, 256) * u.kpc\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(6, 4), layout=\"tight\")\n", "\n", "ax.errorbar(\n", " mwdata2[\"r\"],\n", " mwdata2[\"Menc\"],\n", " yerr=(mwdata2[\"Menc_err_neg\"], mwdata2[\"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(r\"$M(" ] }, "metadata": { "image/png": { "height": 390, "width": 589 } }, "output_type": "display_data" } ], "source": [ "r = np.logspace(-3, 3, 256) * u.kpc\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(6, 4), layout=\"tight\")\n", "\n", "ax.errorbar(\n", " mwdata2[\"r\"],\n", " mwdata2[\"Menc\"],\n", " yerr=(mwdata2[\"Menc_err_neg\"], mwdata2[\"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_v2.mass_enclosed(R=r)\n", "ax.loglog(r.value, 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(r\"$M(" ] }, "metadata": { "image/png": { "height": 390, "width": 589 } }, "output_type": "display_data" } ], "source": [ "mwpot_v1 = gp.MilkyWayPotential(version=\"v1\")\n", "\n", "r_grid = np.linspace(0.5, 30, 128)\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(6, 4), layout=\"tight\")\n", "ax.plot(\n", " r_grid,\n", " fit_potential_v2.circular_velocity(R=r_grid),\n", " marker=\"\",\n", " lw=2,\n", " label=\"v2\",\n", " color=\"tab:blue\",\n", ")\n", "ax.plot(\n", " r_grid,\n", " mwpot_v1.circular_velocity(R=r_grid),\n", " marker=\"\",\n", " label=\"v1\",\n", " linestyle=\"--\",\n", " color=\"tab:orange\",\n", ")\n", "ax.scatter(\n", " eilers[\"R\"][eilers[\"R\"] < 15],\n", " eilers[\"v_c\"][eilers[\"R\"] < 15],\n", " marker=\"o\",\n", " s=8,\n", " color=\"k\",\n", ")\n", "ax.legend(loc=\"best\", fontsize=16)\n", "ax.set_xlabel(\"$R$ [kpc]\")\n", "ax.set_ylabel(\"$v_c$\")" ] }, { "cell_type": "markdown", "id": "7cceab20", "metadata": {}, "source": [ "Note that the parameters of this fit were further tweaked (as described in [Hunt et al. 2022](https://ui.adsabs.harvard.edu/abs/2022MNRAS.516L...7H/abstract)) to provide a better match to the vertical phase spiral morphology." ] } ], "metadata": { "jupytext": { "custom_cell_magics": "kql" }, "kernelspec": { "display_name": "gala", "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.12.12" } }, "nbformat": 4, "nbformat_minor": 5 }