{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ed71de7e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:28.462971Z",
     "iopub.status.busy": "2025-07-31T14:29:28.462722Z",
     "iopub.status.idle": "2025-07-31T14:29:28.877557Z",
     "shell.execute_reply": "2025-07-31T14:29:28.876987Z"
    },
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fa859498",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:28.879503Z",
     "iopub.status.busy": "2025-07-31T14:29:28.879284Z",
     "iopub.status.idle": "2025-07-31T14:29:29.044105Z",
     "shell.execute_reply": "2025-07-31T14:29:29.043575Z"
    },
    "nbsphinx": "hidden"
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'ZMQInteractiveShell' object has no attribute 'magic'",
     "output_type": "error",
     "traceback": [
      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
      "\u001b[31mAttributeError\u001b[39m                            Traceback (most recent call last)",
      "\u001b[36mFile \u001b[39m\u001b[32m~/work/gala/gala/docs/tutorials/nb_setup:1\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m1\u001b[39m \u001b[43mget_ipython\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m.\u001b[49m\u001b[43mmagic\u001b[49m(\u001b[33m'\u001b[39m\u001b[33mconfig InlineBackend.figure_format = \u001b[39m\u001b[33m\"\u001b[39m\u001b[33mretina\u001b[39m\u001b[33m\"\u001b[39m\u001b[33m'\u001b[39m)  \u001b[38;5;66;03m# noqa\u001b[39;00m\n\u001b[32m      3\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mmatplotlib\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mpyplot\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mplt\u001b[39;00m\n\u001b[32m      5\u001b[39m plt.style.use(\u001b[33m\"\u001b[39m\u001b[33mdefault\u001b[39m\u001b[33m\"\u001b[39m)\n",
      "\u001b[31mAttributeError\u001b[39m: 'ZMQInteractiveShell' object has no attribute 'magic'"
     ]
    }
   ],
   "source": [
    "%run nb_setup"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "394aaa0d",
   "metadata": {},
   "source": [
    "# Defining a Milky Way potential model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4f932b92",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:29.045813Z",
     "iopub.status.busy": "2025-07-31T14:29:29.045650Z",
     "iopub.status.idle": "2025-07-31T14:29:29.808604Z",
     "shell.execute_reply": "2025-07-31T14:29:29.808034Z"
    }
   },
   "outputs": [],
   "source": [
    "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",
    "import gala.potential as gp\n",
    "from gala.units import galactic\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19d56fe4",
   "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": "071ec9f9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:29.810631Z",
     "iopub.status.busy": "2025-07-31T14:29:29.810383Z",
     "iopub.status.idle": "2025-07-31T14:29:29.818600Z",
     "shell.execute_reply": "2025-07-31T14:29:29.818162Z"
    }
   },
   "outputs": [],
   "source": [
    "tbl = ascii.read(\"data/MW_mass_enclosed.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3ffc325f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:29.820169Z",
     "iopub.status.busy": "2025-07-31T14:29:29.819963Z",
     "iopub.status.idle": "2025-07-31T14:29:29.825826Z",
     "shell.execute_reply": "2025-07-31T14:29:29.825341Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><i>Table length=16</i>\n",
       "<table id=\"table139906022094704\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>r</th><th>Menc</th><th>Menc_err_neg</th><th>Menc_err_pos</th><th>ref</th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>str27</th></tr></thead>\n",
       "<tr><td>0.01</td><td>30000000.0</td><td>10000000.0</td><td>10000000.0</td><td>Feldmeier et al. (2014)</td></tr>\n",
       "<tr><td>0.12</td><td>800000000.0</td><td>200000000.0</td><td>200000000.0</td><td>Launhardt et al. (2002)</td></tr>\n",
       "<tr><td>8.1</td><td>89502860861.52429</td><td>4994562473.797714</td><td>4858963492.608627</td><td>Bovy et al. (2012)</td></tr>\n",
       "<tr><td>8.3</td><td>110417867208.19055</td><td>4475949382.696884</td><td>4387023236.020782</td><td>McMillan (2011)</td></tr>\n",
       "<tr><td>8.4</td><td>102421035406.90356</td><td>16733918715.629944</td><td>15468328224.531876</td><td>Koposov et al. (2010)</td></tr>\n",
       "<tr><td>19.0</td><td>208023299175.30438</td><td>44317988008.38101</td><td>34833267089.920685</td><td>Kuepper et al. (2015)</td></tr>\n",
       "<tr><td>50.0</td><td>539884832748.48975</td><td>19995734543.31433</td><td>268490735257.4718</td><td>Wilkinson &amp; Evans (1999)</td></tr>\n",
       "<tr><td>50.0</td><td>529886965173.18726</td><td>9997867269.659302</td><td>38536752776.21277</td><td>Sakamoto et al. (2003)</td></tr>\n",
       "<tr><td>50.0</td><td>399914690706.92847</td><td>109976539940.53711</td><td>72696676468.2511</td><td>Smith et al. (2007)</td></tr>\n",
       "<tr><td>50.0</td><td>419910425325.7268</td><td>39991469076.968506</td><td>38172735113.64386</td><td>Deason et al. (2012)</td></tr>\n",
       "<tr><td>60.0</td><td>399914690957.5188</td><td>69985070910.63354</td><td>64344945146.92987</td><td>Xue et al. (2008)</td></tr>\n",
       "<tr><td>80.0</td><td>689852841359.0248</td><td>299936018002.3314</td><td>110361048549.1029</td><td>Gnedin et al. (2010)</td></tr>\n",
       "<tr><td>100.0</td><td>1399701417307.7747</td><td>899808054059.3811</td><td>831336271726.1903</td><td>Watkins et al. (2010)</td></tr>\n",
       "<tr><td>120.0</td><td>539884832260.29584</td><td>199957345314.72906</td><td>123854764645.32489</td><td>Battaglia et al. (2005)</td></tr>\n",
       "<tr><td>150.0</td><td>750000000000.0</td><td>250000000000.0</td><td>250000000000.0</td><td>Deason et al. (2012)</td></tr>\n",
       "<tr><td>200.0</td><td>679854974257.2006</td><td>409912558030.05396</td><td>313652012195.06256</td><td>Bhattacherjee et al. (2014)</td></tr>\n",
       "</table></div>"
      ],
      "text/plain": [
       "<Table length=16>\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": "be9d818a",
   "metadata": {},
   "source": [
    "Let's now plot the above data and uncertainties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f73b651e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:29.827347Z",
     "iopub.status.busy": "2025-07-31T14:29:29.827185Z",
     "iopub.status.idle": "2025-07-31T14:29:30.396292Z",
     "shell.execute_reply": "2025-07-31T14:29:30.395721Z"
    },
    "lines_to_end_of_cell_marker": 2
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<>:21: SyntaxWarning: invalid escape sequence '\\o'\n",
      "<>:21: SyntaxWarning: invalid escape sequence '\\o'\n",
      "/tmp/ipykernel_3838/2928631988.py:21: SyntaxWarning: invalid escape sequence '\\o'\n",
      "  ax.set_ylabel(\"$M(<r)$ [M$_\\odot$]\")\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAGGCAYAAAB/gCblAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKpBJREFUeJzt3X10VPWdx/HPzGASHqdgaEIgIYrIaUASCEmk1WMioWlUInSlVKiEgNSHgHimu67srrB7jsoeWwFXZtdVFyMWdgNWY1SCHrLJco7EBpKGItEqZ2PhAAlQJCFpCDIz+4dl1gt5mqfMQ96vc+bIvXPzu98fQ+bjvb/fvdfkcrlcAgDgL8zBLgAAEFoIBgCAAcEAADAgGAAABgQDAMCAYAAAGBAMAAADggEAYDAk2AWEOqfTqZMnT2rkyJEymUzBLgcAvOJyuXThwgUlJCTIbO79mIBg6MPJkyeVmJgY7DIAwC+OHz+uCRMm9LoNwdCHkSNHSvrmL3PUqFFBrgZAOPrqq69UXV2t7OxsjR49Oig1tLW1KTEx0f2d1huCoQ9XTh+NGjWKYADgFYfDoWHDhoXE90h/Tokz+AwAMCAYACDAKioqVFJSooqKimCX0i+cSgKAACovL9eSJUtkNpu1e/dujRgxQgUFBcEuq1ccMQBAAFVVVcliscjpdMpisai6ujrYJfWJYACAAMrJyZHD4ZDZbJbD4VB2dnawS+oTwQAAAVRQUKDt27crPz9f27dvD/nTSFIEBsOCBQs0evRo3Xfffe51x48fV3Z2tlJSUjR9+nTt2rUriBUCGGzy8/NVWFio/Pz8YJfSLxEXDGvWrNG2bdsM64YMGaLNmzersbFRH374oR5//HF1dHQEqUIACG0RFwzZ2dnXXNk3btw4paWlSZLi4+MVGxurc+fOBaE6AAh9IRUM+/bt07x585SQkCCTyaSysrJrtrHb7UpOTlZMTIyysrJUW1vr0T7q6urkcDi4/xEA9CCkgqGjo0Opqamy2+3dvl9aWiqbzab169ervr5eqampysvL0+nTp/vV/rlz57R06VK9/PLL/iwbACJKSF3glp+f3+vgzMaNG7Vy5UoVFRVJkl566SW9//772rp1q5588sle2+7q6tL8+fP15JNP6vvf/36v23V1dbmX29raPOwFAIS3kDpi6M2lS5dUV1en3Nxc9zqz2azc3FzV1NT0+rMul0vLli3TnXfeqQceeKDXbTds2CCr1ep+ccoJwEDp7OzUkSNH1NnZ6dF7/hY2wXD27Fk5HA7FxcUZ1sfFxam5udm9nJubq4ULF2r37t2aMGGCampq9NFHH6m0tFRlZWVKS0tTWlqaDh8+3O1+1q5dq9bWVvfr+PHjAe0XAFxx8eJFNTY26uLFix69528hdSrJH/bu3dvteqfT2a+fj46OVnR0tOx2u+x2uxwOhz/LA4CQFzZHDLGxsbJYLGppaTGsb2lpUXx8vN/3V1xcrMbGRh04cMDvbQNAKAubYIiKilJ6eroqKyvd65xOpyorKzV79uwgVgYAkSWkgqG9vV0NDQ1qaGiQJDU1NamhoUHHjh2TJNlsNr3yyit6/fXX9emnn+qRRx5RR0eHe5aSP9ntdqWkpCgjI8PvbQMIL54M/PoySBwqz20IqTGGgwcPKicnx71ss9kkSYWFhSopKdGiRYt05swZrVu3Ts3NzUpLS9OePXuuGZD2h+LiYhUXF6utrU1Wq9Xv7QMIH1cGfhMSEjR06FC/bfttofTchpAKhuzsbLlcrl63WbVqlVatWjVAFQGAZyoqKrR9+3Z9/fXXWrx4cb9/7spzGxwOh/u5DcEKhpA6lQQA4ezK//Xv2bNHS5YsUXl5eb9/tq/nNgzkaSaCoQeMMQDwlC9Pa+vtuQ2+BI43CIYeMF0VgKf6+7S2ngaoe3puw0A/HpRgAAA/6e/T2jy9inmgHw8aUoPPABDu8vPzdd111xnu6+arK4GzY8cOLV68OOCD0gRDD7glBoArvJ1p5M/2AhE4PSEYesB1DACk/x/4leSX6wtC6XqFnjDGAAC9ePXVVw3L//Ef/+FTewM9kOwNggEABtBADyR7g2DoAdcxAJCkBx98UJJkMpkkSStWrPC4jZiYGKWkpCgmJqbXmUvf3i6YGGPoAWMMAKRvZgS98847qq6uVnZ2tlfjAUOHDtXUqVPdyz0NJF+9XbAQDADQh7lz52rSpEm68cYbJX1zgdr//u//6sYbb/ToRnnhglNJANCHnTt36oknntDOnTslDexjNoOBIwYA6EV5ebmWLVsmk8mk3bt3a/To0Wpvb/frdQ2hhmAAgF5cma565ZEAzzzzjGpra0P6OgRfcSqpB8xKAtCds2fPhvx1CL4iGHrA3VUBSP8/XfWKn/70p71ehxAqU059wakkAOjF1TewW7x4saZOndrjDe1CZcqpLwgGAOjD1dcdDOQN7YKBU0kAMMBC/XQTRwwA4CFfv9hD/XQTwQAAHgr1L3ZfcSqpB0xXBTBYEQw9YLoqgCtCfUzA3ziVBAB9iPRTR1fjiAEAwsBAHrVwxAAAYWAgj1o4YgAAGBAMACJCZ2enjhw5os7OzmCXEvYIBgARIdIfnjOQCAYAEaGiokIlJSWqqKgIdilhj8FnAGGvvLxcS5YsieiH5wwkjhgAhL2qqqqIf3jOQCIYesAtMYDwkZOT0+vDc+AZk+vKg0zRrba2NlmtVrW2tmrUqFHBLgdAD3bs2GF4mA6MPPkuY4wBQESI9IfnDCROJQEADAgGAIABwQAgIgy2W2MHEmMMACLCYLs1diBxxAAAMCAYAAAGBAMAwIBgAAAYRFwwLFiwQKNHj9Z9993Xr/UAAKOIC4Y1a9Zo27Zt/V4PADCKuGDIzs7WyJEj+70eQHDs2rVLS5cu1a5du4JdCq4SUsGwb98+zZs3TwkJCTKZTCorK7tmG7vdruTkZMXExCgrK0u1tbUDXygAn5SXl+snP/mJtm/frp/85CcqLy8Pdkn4lpAKho6ODqWmpsput3f7fmlpqWw2m9avX6/6+nqlpqYqLy9Pp0+fHuBKAfiiqqpKZrNZTqdTZrOZ5yeEmJAKhvz8fD399NNasGBBt+9v3LhRK1euVFFRkVJSUvTSSy9p2LBh2rp16wBXCsAXw4YNk9PplCQ5nU4NHTo0yBXh20IqGHpz6dIl1dXVGW6pazablZubq5qaGr/tp6urS21tbYYXAP/685//LLP5m68fs9mszs7OIFeEbwubYDh79qwcDofi4uIM6+Pi4tTc3Oxezs3N1cKFC7V7925NmDDBHRo9rb/ahg0bZLVa3a/ExMTAdQoYpK4cMZhMJjmdTp64FmIi7iZ6e/fu9Wj91dauXSubzeZebmtrIxwAPyovL9ezzz7rHmOw2WwqKCgIdln4lrAJhtjYWFksFrW0tBjWt7S0KD4+3m/7iY6OVnR0tOx2u+x2uxwOh9/aBvDNwLPFYnE/o5nTSKEnbE4lRUVFKT09XZWVle51TqdTlZWVmj17tt/3V1xcrMbGRh04cMDvbQODWU5OjjsUnE6nbrvttmCXhKuE1BFDe3u7jh496l5uampSQ0ODxowZo6SkJNlsNhUWFmrWrFnKzMzU5s2b1dHRoaKioiBWDcATBQUF2rlzp9555x3NmjWrx1mICB6Ty+VyBbuIK6qrq5WTk3PN+sLCQpWUlEiStmzZol/+8pdqbm5WWlqa/uVf/kVZWVl+r+Xbp5I+//xztba2atSoUX7fDwAMhLa2Nlmt1n59l4VUMIQiT/4yASBUefJdFjZjDACAgUEwAAAMCIYe2O12paSkKCMjI9ilAMCAYoyhD4wxAIgEjDEAALxGMAAADAiGHjDGAGCwYoyhD4wxAIgEjDEAALxGMAAADAiGHjDGAGCwYoyhD4wxAIgEjDEAALxGMAAADAgGAIABwQAAMCAYesCsJACDFbOS+sCsJACRgFlJAACvEQwAAAOCAQBgQDAAAAwIBgCAAcHQA6arAhismK7aB6arAogETFcFAHiNYAAAGBAMAAADggEAYEAwAAAMCAYAgAHBAAAwIBgAAAYEAwDAgGDoAbfEADBYcUuMPnBLDACRgFtiAAC8NsSTjcvLyz3ewdy5czV06FCPfw4AEBweBcP8+fM9atxkMumLL77QjTfe6NHPAQCCx+NTSc3NzXI6nf16DRs2LBA1AwACyKNgKCws9Oi00M9+9jMGbAEgzDArqQ/MSgIQCZiVBADwGsEAADDwazB0dHTod7/7nc6ePevPZgEAA8hvwfDss8/qxz/+sd566y09+uijWr58uTo7O/3VfL8sWLBAo0eP1n333WdY/95772nKlCmaPHmyXn311QGtCQDCjUfXMfTktdde07lz5/TBBx+41+3atUuPP/64/v3f/90fu+iXNWvWaPny5Xr99dfd6y5fviybzaaqqipZrValp6drwYIFuv766wesLgAIJ345YigtLdWTTz4pSVq+fLlaWlq0cOFCHThwwB/N91t2drZGjhxpWFdbW6upU6dq/PjxGjFihPLz8/Xhhx8OaF0AEE78EgzXXXedLl26JEn6wQ9+4L6wzel09ruNffv2ad68eUpISJDJZFJZWdk129jtdiUnJysmJkZZWVmqra3ts92TJ09q/Pjx7uXx48frxIkT/a4LAAYbvwTDgw8+qCeeeEIOh0MrVqzQyJEjtWnTJuXn5/e7jY6ODqWmpsput3f7fmlpqWw2m9avX6/6+nqlpqYqLy9Pp0+f9kcXAAB/4ZcxhnvvvVcdHR2aM2eOJk6cqFOnTun73/++nn766X63kZ+f32uQbNy4UStXrlRRUZEk6aWXXtL777+vrVu3uk9jdSchIcFwhHDixAllZmb2uH1XV5e6urrcy21tbf3uAwBEAr8EgyQtXrxYixcv1rlz5zR69GiZTCZ/Na1Lly6prq5Oa9euda8zm83Kzc1VTU1Nrz+bmZmpTz75RCdOnJDValVFRYWeeuqpHrffsGGD/umf/slvtQNAuPH7BW5jxozxayhI0tmzZ+VwOBQXF2dYHxcXp+bmZvdybm6uFi5cqN27d2vChAmqqanRkCFD9PzzzysnJ0dpaWn6xS9+0euMpLVr16q1tdX9On78uF/7AgChzm9HDKFg79693a4vKChQQUFBv9qIjo5WdHS0P8sCgLASFrfEiI2NlcViUUtLi2F9S0uL4uPjA7JPnvkMYLDy+Ihh+fLl/dpu69atHhfTk6ioKKWnp6uystL9sCCn06nKykqtWrXKb/v5tuLiYhUXF7vvSAgAg4XHwVBSUqKJEydqxowZ8ucdu9vb23X06FH3clNTkxoaGjRmzBglJSXJZrOpsLBQs2bNUmZmpjZv3qyOjg73LCUAgJ+4PPToo4+6Ro8e7UpLS3O98MILrj/96U+eNtGtqqoql6RrXoWFhe5tXnzxRVdSUpIrKirKlZmZ6fr444/9su/ubNmyxfW9733PdfPNN7skuVpbWwO2LwAItNbW1n5/l3n1oJ6uri699dZb2rp1q/bv36+7775bK1as0A9/+EO/z0gKNh7UAyASePJd5vMT3P74xz+qpKRE27Zt0+XLl3XkyBGNGDHClyZDCsEAIBIM6BPczGazTCaTXC6XHA6Hr80BAILMq2Do6urSf/7nf2ru3Lm6+eabdfjwYW3ZskXHjh2LmKMFpqsCGKw8PpX06KOP6r/+67+UmJio5cuXa8mSJYqNjQ1UfUHHqSQAkSCgYwxms1lJSUmaMWNGrwPNb731lifNhiyCAUAk8OS7zOPrGJYuXRpxM4+6Y7fbZbfbGTcBMOj4PCsp0nHEACASDOisJABAZPEoGH7/+9979LjOI0eO6PLlyx4XBQAIHo+CYcaMGfrTn/7U7+1nz56tY8eOeVwUACB4PBp8drlceuqppzRs2LB+bX/p0iWvigoFDD4DGKw8GnzOzs72eEbSjh07NG7cOI8LCxUMPgOIBAGbrlpdXe1LXQCAMMCsJACAAcEAADAgGAAABgRDD7i7KoDByqdbYnz99ddqbm7Wn//8Z40dO1ZjxozxZ20hgVlJACJBQG+JceHCBf3bv/2b7rjjDo0aNUrJycn63ve+p7Fjx2rixIlauXKlDhw44HXxAIDg8igYNm7cqOTkZL322mvKzc1VWVmZGhoa9Pnnn6umpkbr16/X5cuX9cMf/lA/+tGP9MUXXwSqbgBAgHh0Kun+++/XP/zDP2jq1Km9btfV1aXXXntNUVFRWr58uc9FBhOnkgBEgoA+qEeSHA6H3n33Xc2ZM0cjR470utBwQDAAiAQBv+22xWLR/fffrzNnznhVIAAgdHk9XTUjI0NNTU3+rCWkMF0VwGDl9XTVXbt26Ve/+pXefPNNJSYm+ruukMGpJACRIKDPfL5i0aJFkqSpU6eqoKBA2dnZmjFjhm655RZFRUV52ywAIMi8DoampiYdOnRIDQ0NOnTokDZs2KAvv/xSQ4YM0ZQpU/T73//en3UCAAaI18EwceJETZw4UQUFBe51Fy5cUENDA6EAAGHMp1tiDAaMMQCIBAGfrgoAiFwEAwDAwKtg6Ojo0OHDh7t978iRI2pvb/epKABA8HgVDF9//bWysrJUW1trWN/Y2KgZM2YQDAAQxrwKhu985zu65557tG3bNsP6N954Q3PmzFF8fLxfigMADDyvxxgKCwtVWlqqy5cvS5JcLpe2b9+uoqIivxUXTNwSA8Bg5XUw/OhHP9KQIUP0/vvvS5Kqq6vV3t6u+fPn+6u2oCouLlZjYyMPHQIw6HgdDBaLRUuWLHGfTnrjjTe0aNEibocBAGHO6yufpW9OJ2VmZurEiRP6zW9+ow8++MBfdQEAgsTnK5/T09M1cuRINTc367PPPvNXXSGDK58BRIIBvfJ56dKl2rdvn5YuXeprUwCAEODTqSRJeuCBB3T+/Pmwf7YzAOAb3ESvD5xKAhAJuIkeAMBrBAMAwIBgAAAYDIpg+NWvfqWpU6dq2rRp+vWvfx3scgAgpHkUDMeOHfOo8RMnTni0fSAcPnxYO3bsUF1dnQ4cOKAtW7bo/PnzwS4LAEKWR8GQkZGhhx56qNf7B7W2tuqVV17RtGnT9Jvf/MbnAn316aefavbs2YqJidHQoUOVmpqqPXv2BLssAAhZHgVDY2Ojhg8frrlz5yo+Pl533323Vq5cqdWrV+tnP/uZZs6cqe9+97vaunWrnnvuOT322GM+F7hv3z7NmzdPCQkJMplMKisru2Ybu92u5ORkxcTEXPOciGnTpqm6ulrnz5/XV199perq6pA4kgGAUOVRMFx//fXauHGjTp06pS1btmjy5Mk6e/asvvjiC0nSkiVLVFdXp5qaGt11111+KbCjo0Opqamy2+3dvl9aWiqbzab169ervr5eqampysvL0+nTpyVJKSkpeuyxx3TnnXfqxz/+sW699VZZLBa/1AYAkcjjC9yefvpppaena+bMmYqLiwtUXd0ymUx6++23Dbf2zsrKUkZGhrZs2SJJcjqdSkxM1OrVq/Xkk09e08aDDz6oBQsW6O677+52H11dXerq6nIvt7W1KTExkQvcAIS1gF7gtm7dOt1zzz1KSEjQ+PHjNW/ePP3jP/6j3n333QE/RXPp0iXV1dUpNzfXvc5sNis3N1c1NTXudVeOHv7whz+otrZWeXl5Pba5YcMGWa1W9ysxMTFwHcCA6ezs1JEjR9TZ2RnsUoCQ5/G9kjIyMnTq1CkVFRUpNjZW9fX1euutt/Tss8/K4XBo7Nixmjlzpnbv3h2Ieg3Onj0rh8NxzZFLXFyc4U6v9957r1pbWzV8+HC99tprGjKk526vXbtWNpvNvXzliAHh7eLFi2psbFRCQoKGDh0a7HKAkOZxMPz2t79VSUmJ/u7v/k4ZGRnauHGjJk2apK6uLjU0NKi+vl6/+93vAlGr17599NCX6OhoRUdHB7AaAAhtXl3gtmzZMn3++eeaMmWKZs6cqbVr18rhcCgrK0uPPPKIXn75ZX/X2a3Y2FhZLBa1tLQY1re0tCg+Pt6ntnnmM4DByusrn0eMGKHnnntOBw8e1CeffKKbbrrJ/ZjPgRIVFaX09HRVVla61zmdTlVWVmr27Nk+tc0znyNLRUWFSkpKVFFREexSgJDn0/MYLl++rK6uLt1///1qaWlRUVGR7rnnHo0ZM8Zf9am9vV1Hjx51Lzc1NamhoUFjxoxRUlKSbDabCgsLNWvWLGVmZmrz5s3q6OhQUVGR32pAeCsvL9eSJUtkNpu1e/dujRgxQgUFBcEuCwhZHgfDP//zP+vw4cM6fPiwPvvsM8XExGj69OnKzMzUQw89JKvV6tcCDx48qJycHPfylYHhwsJClZSUaNGiRTpz5ozWrVun5uZmpaWlac+ePT5PpbXb7bLb7XI4HD61g+CrqqqSxWKRw+GQxWJRdXU1wQD0wuPrGMxms5KTk1VYWKj7779fN998c6BqCwk8qCf8lZeX695775XZbJbT6dQ777xDMGDQ8eS7zONguOOOO9TQ0KALFy5o+PDhmj59umbOnOl+TZs2LaKuLCYYIsOOHTu0Y8cOLV68WIsXLw52OcCA8+S7zONTSf/zP/8jSfriiy9UV1en+vp61dfXa/v27Tp//ryio6N1yy23GO5XFI44lRRZ8vPzdd111xkuhgTQPa8HnydPnqzJkyfrpz/9qXtdU1OTDh48GHLXMXijuLhYxcXF7pQFgMHCp1lJV7vhhht0ww03aOHChf5sFgAwgAbFE9wAAP1HMPSAK58jS0xMjFJSUhQTExPsUoCQ5/GspMGGWUkAIkFAb7sNAIhsBAMAwIBgAAAYEAw9YPAZwGDF4HMfGHwGEAkYfAYAeI1gAAAYEAwAAAOCAQBgQDD0gFlJAAYrZiX1gVlJACIBs5IAAF4jGAAABgQDAMCAYAAAGBAMAAADgqEHTFcFMFgxXbUPTFcFEAmYrgoA8BrBAAAwIBgAAAYEAwDAgGAAABgQDAAAA4IBAGBAMAAADAiGHnDlM4DBiiuf+8CVzwAiAVc+AwC8RjAAAAwIBgCAAcEAADAgGAAABgQDAMCAYAAAGBAMAAADggEAYEAwAAAMBkUwbNq0SVOnTlVKSooee+wxcRcQAOhZxAfDmTNntGXLFtXV1enw4cOqq6vTxx9/HOyyACBkDQl2AQPh8uXLunjxoiTp66+/1ne/+90gVwQAoSvkjxj27dunefPmKSEhQSaTSWVlZddsY7fblZycrJiYGGVlZam2ttb93tixY/XXf/3XSkpKUkJCgnJzczVp0qQB7AEAhJeQD4aOjg6lpqbKbrd3+35paalsNpvWr1+v+vp6paamKi8vT6dPn5YkffXVV3rvvff05Zdf6sSJE9q/f7/27ds3kF0AgLAS8sGQn5+vp59+WgsWLOj2/Y0bN2rlypUqKipSSkqKXnrpJQ0bNkxbt26VJO3du1c33XSTxowZo6FDh+ruu+/udYyhq6tLbW1thhcADCYhHwy9uXTpkurq6pSbm+teZzablZubq5qaGklSYmKi9u/fr4sXL8rhcKi6ulpTpkzpsc0NGzbIarW6X4mJiQHvBwCEkrAOhrNnz8rhcCguLs6wPi4uTs3NzZKkW2+9VXfddZdmzJih6dOna9KkSSooKOixzbVr16q1tdX9On78eED7AAChZlDMSnrmmWf0zDPP9Gvb6OhoRUdHB7giAAhdYX3EEBsbK4vFopaWFsP6lpYWxcfH+9S23W5XSkqKMjIyfGoHAMJNWAdDVFSU0tPTVVlZ6V7ndDpVWVmp2bNn+9R2cXGxGhsbdeDAAV/LBICwEvKnktrb23X06FH3clNTkxoaGjRmzBglJSXJZrOpsLBQs2bNUmZmpjZv3qyOjg4VFRUFsWoACF8hHwwHDx5UTk6Oe9lms0mSCgsLVVJSokWLFunMmTNat26dmpublZaWpj179lwzIO0pu90uu90uh8PhUzsAEG5MLu4o16u2tjZZrVa1trZq1KhRwS4HALziyXdZWI8xAAD8j2DoAbOSAAxWnErqA6eSAEQCTiUBALxGMAAADAiGHjDGAGCwYoyhD4wxAIgEjDEAALxGMAAADAgGAIABwdADBp8BDFYMPveBwWcAkYDBZwCA1wgGAIABwQAAMCAYAAAGBEMPBsOspM7OTh05ckSdnZ3BLgVACCEYelBcXKzGxkYdOHAg2KUEzMWLF9XY2KiLFy8GuxQAIYRgAAAYEAwAAAOCAQBgQDAMYhUVFSopKVFFRUWwSwEQQoYEuwAER3l5uZYsWSKz2azdu3drxIgRKigoCHZZAEIARww9iPTpqlVVVbJYLHI6nbJYLKqurg52SQBCBMHQg0ifrpqTkyOHwyGz2SyHw6Hs7OxglwQgRBAMg1RBQYG2b9+u/Px8bd++ndNIANwYYxjE8vPzdd111yk3NzfYpQAIIRwxAAAMCAYAgAHBAAAwIBgAAAYEAwDAgGAAABgQDD2I9CufJSkmJkYpKSmKiYkJdikAQojJ5XK5gl1EKGtra5PValVra6tGjRoV7HIAwCuefJdxxAAAMCAYAAAGBAMAwIBgAAAYEAwAAAOCAQBgQDAAAAwIBgCAAcEAADAgGAAABoMiGP7whz8oLS3N/Ro6dKjKysqCXRYAhKRB8cznKVOmqKGhQZLU3t6u5ORkzZ07N7hFAUCIGhRHDN9WXl6uOXPmaPjw4cEuBQBCUlgEw759+zRv3jwlJCTIZDJ1exrIbrcrOTlZMTExysrKUm1tbbdt7dy5U4sWLQpwxQAQvsIiGDo6OpSamiq73d7t+6WlpbLZbFq/fr3q6+uVmpqqvLw8nT592rBdW1ub9u/fr7vuumsgygaAsBR2z2MwmUx6++23NX/+fPe6rKwsZWRkaMuWLZIkp9OpxMRErV69Wk8++aR7uzfeeEMffPCBfv3rX/fYfldXl7q6utzLbW1tSkxM5HkMAMLaoHoew6VLl1RXV6fc3Fz3OrPZrNzcXNXU1Bi27c9ppA0bNshqtbpfiYmJAakbAEJV2AfD2bNn5XA4FBcXZ1gfFxen5uZm93Jra6tqa2uVl5fXa3tr165Va2ur+3X8+PGA1A0AoWpQTFeVJKvVqpaWlj63i46OVnR09ABUBAChKeyPGGJjY2WxWK750m9paVF8fLzX7drtdqWkpCgjI8PXEgEgrIR9MERFRSk9PV2VlZXudU6nU5WVlZo9e7bX7RYXF6uxsVEHDhzwR5kAEDbC4lRSe3u7jh496l5uampSQ0ODxowZo6SkJNlsNhUWFmrWrFnKzMzU5s2b1dHRoaKioiBWDQDhKSyC4eDBg8rJyXEv22w2SVJhYaFKSkq0aNEinTlzRuvWrVNzc7PS0tK0Z8+eawakPWG322W32+VwOHyuHwDCSdhdxzDQPJn7CwChypPvsrA4YgimK7nZ1tYW5EoAwHtXvsP6cyxAMPTgyqmkK1dBc6EbgEhw4cIFWa3WXrfhVFIfzp8/r9GjR+vYsWN9/mX2JiMjo18znHrarr/re1u++s+VlZVKTEzU8ePHfTpN1p++9bZNd+/1ta6nfl7575VbmfjSN18/s/70o7t1ff05nPoWqv8ee9suUn/XXC6XLly4oISEBJnNvU9I5YihD1f+Aq1Wq08fqMVi6dfP97Rdf9f3ttzTn0eNGhXwvvW2TXfv9bWup35evd6Xvvn6mfX0nid9663P4dC3UP332Nt2kfy71t//uQ376xjCRXFxsU/b9Xd9b8s9/dlX/Wmrt226e6+vdT31c6D71dd2vvYtmJ9ZX9v1p2+h+u+xt+0i+XetvziV1IdInpVE38JTpPYtUvslhV/fOGLoQ3R0tNavXx+R90+ib+EpUvsWqf2Swq9vHDEAAAw4YgAAGBAMAAADggEAYEAwAAAMCAY/On/+vGbNmqW0tDRNmzZNr7zySrBL8pvjx48rOztbKSkpmj59unbt2hXskvxmwYIFGj16tO67775gl+Kz9957T1OmTNHkyZP16quvBrscv4qkz+nbQvF3i1lJfuRwONTV1aVhw4apo6ND06ZN08GDB3X99dcHuzSfnTp1Si0tLUpLS1Nzc7PS09P1+eefa/jw4cEuzWfV1dW6cOGCXn/9db355pvBLsdrly9fVkpKiqqqqmS1WpWenq79+/dHxL8/KXI+p6uF4u8WRwx+ZLFYNGzYMElSV1eXXC5Xv+5kGA7GjRuntLQ0SVJ8fLxiY2N17ty54BblJ9nZ2Ro5cmSwy/BZbW2tpk6dqvHjx2vEiBHKz8/Xhx9+GOyy/CZSPqerheLv1qAKhn379mnevHlKSEiQyWRSWVnZNdvY7XYlJycrJiZGWVlZqq2t9Wgf58+fV2pqqiZMmKC/+Zu/UWxsrJ+q791A9O2Kuro6ORyOAbnj7ED2K9h87evJkyc1fvx49/L48eN14sSJgSi9T5H8OfqzbwP5u9WbQRUMHR0dSk1Nld1u7/b90tJS2Ww2rV+/XvX19UpNTVVeXp5Onz7t3ubK+MHVr5MnT0qSvvOd7+jQoUNqamrSjh071NLSEjF9k6Rz585p6dKlevnllwPeJ2ng+hUK/NHXUEXf+u7bQP9u9co1SElyvf3224Z1mZmZruLiYveyw+FwJSQkuDZs2ODVPh555BHXrl27fCnTK4Hq28WLF1233367a9u2bf4q1SOB/Myqqqpcf/VXf+WPMv3Cm75+9NFHrvnz57vfX7NmjWv79u0DUq8nfPkcQ+1zupq3fQv279bVBtURQ28uXbqkuro65ebmuteZzWbl5uaqpqamX220tLTowoULkqTW1lbt27dPU6ZMCUi9nvBH31wul5YtW6Y777xTDzzwQKBK9Yg/+hUu+tPXzMxMffLJJzpx4oTa29tVUVGhvLy8YJXcb5H8Ofanb6H4u0Uw/MXZs2flcDgUFxdnWB8XF6fm5uZ+tfHHP/5Rt99+u1JTU3X77bdr9erVuuWWWwJRrkf80bePPvpIpaWlKisrU1pamtLS0nT48OFAlNtv/uiXJOXm5mrhwoXavXu3JkyYEJJfRv3p65AhQ/T8888rJydHaWlp+sUvfhEWM5L6+zmGw+d0tf70LRR/t3hQjx9lZmaqoaEh2GUExG233San0xnsMgJi7969wS7BbwoKClRQUBDsMgIikj6nbwvF3y2OGP4iNjZWFovlmsHilpYWxcfHB6kq/4jUvkVqv7oTyX2lb6GHYPiLqKgopaenq7Ky0r3O6XSqsrJSs2fPDmJlvovUvkVqv7oTyX2lb6FnUJ1Kam9v19GjR93LTU1Namho0JgxY5SUlCSbzabCwkLNmjVLmZmZ2rx5szo6OlRUVBTEqvsnUvsWqf3qTiT3lb6FWd+CPS1qIFVVVbkkXfMqLCx0b/Piiy+6kpKSXFFRUa7MzEzXxx9/HLyCPRCpfYvUfnUnkvtK38Krb9wrCQBgwBgDAMCAYAAAGBAMAAADggEAYEAwAAAMCAYAgAHBAAAwIBgAAAYEAwDAgGAAABgQDAAAA4IBCKDs7GyZTCaZTCb3Q5yys7P1+OOPB2yfy5Ytc++zrKwsYPtB5CIYgABbuXKlTp06pWnTpg3I/l544QWdOnVqQPaFyDSonscABNLly5c1ZMi1v1LDhg0b0Kd1Wa1WWa3WAdsfIg9HDIAXvvzyS5lMJu3cuVO33367oqOjVV5e7lVb77//vqxWq7Zv3y7pm1NNq1at0qpVq2S1WhUbG6unnnpK375DvtPp1HPPPaebbrpJ0dHRSkpK0jPPPOOXvgEcMQBeOHTokCTpl7/8pZ599lndcMMNGjt2rMft7NixQw8//LB27Nihe+65x73+9ddf14oVK1RbW6uDBw/q5z//uZKSkrRy5UpJ0tq1a/XKK69o06ZNuu2223Tq1Cl99tln/ukcBj2CAfBCQ0ODhg8frl27dik5OdmrNux2u/7+7/9e7777ru644w7De4mJidq0aZNMJpOmTJmiw4cPa9OmTVq5cqUuXLigF154QVu2bFFhYaEkadKkSbrtttt87RYgiWAAvHLo0CEVFBR4HQpvvvmmTp8+rY8++kgZGRnXvH/rrbfKZDK5l2fPnq3nn39eDodDn376qbq6ujRnzhxvywd6xRgD4IWGhgZlZ2d7/fMzZszQ2LFjtXXrVnn6dN2hQ4d6vV+gPwgGwENtbW368ssvNWPGDK/bmDRpkqqqqvTOO+9o9erV17z/29/+1rD88ccfa/LkybJYLJo8ebKGDh2qyspKr/cP9IZTSYCHDh06JIvFoltuucWndm6++WZVVVUpOztbQ4YM0ebNm93vHTt2TDabTQ899JDq6+v14osv6vnnn5ckxcTE6G//9m/1xBNPKCoqSj/4wQ905swZHTlyRCtWrPCpJkAiGACPHTp0SFOmTFFMTIzPbU2ZMkX//d//rezsbFksFveX/9KlS9XZ2anMzExZLBatWbNGP//5z90/99RTT2nIkCFat26dTp48qXHjxunhhx/2uR5AkkwuT09wAui37OxspaWlGY4GAvEz3TGZTHr77bc1f/58n9rB4MMYAxBg//qv/6oRI0bo8OHDA7K/hx9+WCNGjBiQfSEyccQABNCJEyfU2dkpSUpKSlJUVFSfP+PrEcPp06fV1tYmSRo3bpyGDx/uVTsYvAgGAIABp5IAAAYEAwDAgGAAABgQDAAAA4IBAGBAMAAADAgGAIABwQAAMCAYAAAGBAMAwIBgAAAY/B8Nfcnx1Lj2MAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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(<r)$ [M$_\\odot$]\")\n",
    "\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f0d00f7",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "We now need to assume some form for the potential. For simplicity and within\n",
    "reason, we'll use a four component potential model consisting of a Hernquist\n",
    "([1990](https://ui.adsabs.harvard.edu/#abs/1990ApJ...356..359H/abstract))\n",
    "bulge and nucleus, a Miyamoto-Nagai\n",
    "([1975](https://ui.adsabs.harvard.edu/#abs/1975PASJ...27..533M/abstract))\n",
    "disk, and an NFW\n",
    "([1997](https://ui.adsabs.harvard.edu/#abs/1997ApJ...490..493N/abstract))\n",
    "halo. We'll fix the parameters of the disk and bulge to be consistent with\n",
    "previous work ([Bovy\n",
    "2015](https://ui.adsabs.harvard.edu/#abs/2015ApJS..216...29B/abstract) -\n",
    "please cite that paper if you use this potential model) and vary the scale\n",
    "mass and scale radius of the nucleus and halo, respectively. We'll fit for\n",
    "these parameters in log-space, so we'll first define a function that returns a\n",
    "`gala.potential.CCompositePotential` object given these four parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "19a38ad6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:30.398164Z",
     "iopub.status.busy": "2025-07-31T14:29:30.397952Z",
     "iopub.status.idle": "2025-07-31T14:29:30.401772Z",
     "shell.execute_reply": "2025-07-31T14:29:30.401322Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_potential(log_M_h, log_r_s, log_M_n, log_a):\n",
    "    mw_potential = gp.CCompositePotential()\n",
    "    mw_potential[\"bulge\"] = gp.HernquistPotential(m=5e9, c=1.0, units=galactic)\n",
    "    mw_potential[\"disk\"] = gp.MiyamotoNagaiPotential(\n",
    "        m=6.8e10 * u.Msun, a=3 * u.kpc, b=280 * u.pc, units=galactic\n",
    "    )\n",
    "    mw_potential[\"nucl\"] = gp.HernquistPotential(\n",
    "        m=np.exp(log_M_n), c=np.exp(log_a) * u.pc, units=galactic\n",
    "    )\n",
    "    mw_potential[\"halo\"] = gp.NFWPotential(\n",
    "        m=np.exp(log_M_h), r_s=np.exp(log_r_s), units=galactic\n",
    "    )\n",
    "\n",
    "    return mw_potential"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6f0d706",
   "metadata": {},
   "source": [
    "We now need to specify an initial guess for the parameters - let's do that (by\n",
    "making them up), and then plot the initial guess potential over the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2077040e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:30.403374Z",
     "iopub.status.busy": "2025-07-31T14:29:30.403214Z",
     "iopub.status.idle": "2025-07-31T14:29:30.414866Z",
     "shell.execute_reply": "2025-07-31T14:29:30.414387Z"
    }
   },
   "outputs": [],
   "source": [
    "# Initial guess for the parameters- units are:\n",
    "#     [Msun, kpc, Msun, pc]\n",
    "x0 = [np.log(6e11), np.log(20.0), np.log(2e9), np.log(100.0)]\n",
    "init_potential = get_potential(*x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "317ef75c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:30.416424Z",
     "iopub.status.busy": "2025-07-31T14:29:30.416256Z",
     "iopub.status.idle": "2025-07-31T14:29:31.061460Z",
     "shell.execute_reply": "2025-07-31T14:29:31.060915Z"
    },
    "lines_to_end_of_cell_marker": 2
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<>:27: SyntaxWarning: invalid escape sequence '\\o'\n",
      "<>:27: SyntaxWarning: invalid escape sequence '\\o'\n",
      "/tmp/ipykernel_3838/1286690888.py:27: SyntaxWarning: invalid escape sequence '\\o'\n",
      "  ax.set_ylabel(\"$M(<r)$ [M$_\\odot$]\")\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAGGCAYAAAB/gCblAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQUpJREFUeJzt3Xl8U3W+P/5Xku5baGlp6c5abIEWugAC0kqRKUoBleEKSi2IGy7z7dw7I3O/wtffz5H5ekfEGTOXGR1E58IM4oIoi0ptLQLajUKhslTK1iVtoWto0zY53z8KBwPd0iY5WV7Px6MPPSenJ+8PafLOZ5cJgiCAiIjoBrnUARARkXVhYiAiIgNMDEREZICJgYiIDDAxEBGRASYGIiIywMRAREQGmBiIiMiAk9QBWDu9Xo+qqip4e3tDJpNJHQ4R0aAIgoCWlhYEBwdDLu+7TsDE0I+qqiqEhYVJHQYRkUlcvnwZoaGhfV7DxNAPb29vAN3/mD4+PhJHQ0S2qKGhAbm5uUhOToavr68kMTQ3NyMsLEz8TOsLE0M/bjYf+fj4MDEQ0aDodDp4eHhYxefIQJrE2flMREQGmBiIiMxs//792LZtG/bv3y91KAPCpiQiIjPas2cPVqxYAblcjn379sHLywvp6ekWj8OYHRZYYyAiMqOcnBwoFAro9XooFArk5uZa9Pk7dXrkV1zDn7J/GvDvMDEQEZlRSkoKdDod5HI5dDodkpOTLfK8bR06fHVKjVc+/xHbf7iMysa2Af8um5KIiMwoPT0d27dvx44dO7B8+XKzNyO1tHfh2zN1OFR+Fe2dukHdw+5qDEuWLIGvry8efvhh8dzly5eRnJyM6OhoTJ48Gbt27ZIwQiJyNGlpacjIyEBaWprZnqNB04GPiyvxyuc/4usfa8WkIIMMcWFKPJM8esD3srsaw4svvohVq1bh/fffF885OTlh8+bNiIuLQ01NDeLj47FgwQJ4enpKGCkR0dBdbe3Al6fUKLzQAN3POpid5DIkjfLDvRMCEODtiubm5gHf0+4SQ3Jy8h2dOyNHjsTIkSMBAEFBQfD398e1a9eYGIjIZjW1deKrU2p8f/4auvS3EoKLQo6ZY4cjOSoAwzycB3Vvq2pKysvLw8KFCxEcHAyZTIbdu3ffcY1KpUJkZCTc3Nwwbdo05OfnG/UcRUVF0Ol0XP+IiGxSq7YLu0uq8P9/cRrflV8Vk4K7swLzYwLxf9LvwuIpwYNOCoCV1Rg0Gg1iY2OxatUqPPjgg3c8vnPnTmRlZWHLli2YNm0aNm/ejPnz5+PMmTMYMWJEv/e/du0aVq5ciXfeeccc4RMRmU1bpw45p+uQe6YO2i69eN7VSY7kqACkRAXA3UVhkueyqsSQlpbWZ+fMpk2bsGbNGmRmZgIAtmzZgr1792Lr1q146aWX+ry3VqvF4sWL8dJLL+Huu+/u8zqtViseG9MuR0Rkajq9gO/PX8O+0hq0arvE805yOWaPH465E0bA2820H+VWlRj60tHRgaKiIqxbt048J5fLkZqaiqNHj/b5u4Ig4PHHH8e9996Lxx57rM9rN27ciFdeecUkMRMRGaOtrQ3nz5/H6NGj4e7ujtPVLfi0pAo1Te3Q6/TQXNfA28sTM8cF4L7owCE1F/XFZhJDfX09dDodAgMDDc4HBgbi9OnT4nFqaiqOHz8OjUaD0NBQ7Nq1CzqdDjt37sTkyZPFfot//OMfmDRp0h3Ps27dOmRlZYnHN5eqJSIyt/b2dpSVlcHZxx/fnK/Gj9Ut4mM6vQ7DdA14dvZ4jA3tv+l8KGwmMQzUwYMHezyv1+t7PH87V1dXuLq6QqVSQaVSQacb3AQRIiJjdXTpcaLRFd9+ewUKp1sfz+F+Hpg7xgvniisw3NM8tYSfs5nE4O/vD4VCAbVabXBerVYjKCjI5M+3du1arF27Fs3NzVAqlSa/PxHRTYIg4MSVZvzr+8v4qcUF/m4CFACGeTgjPXYkpoYPQ2NjI85ZKB6rGq7aFxcXF8THxyM7O1s8p9frkZ2djRkzZkgYGRHR4NW1aLHl2wpsPXwBTe3dncsKuQzzYwLxnwsmID7C1+L7zVtVjaG1tRXl5eXicUVFBUpKSuDn54fw8HBkZWUhIyMDCQkJSEpKwubNm6HRaMRRSqbEpiQiuun2TmFTXKvTC8g5XYf9J9XoutHUXVH8LWqO52Huow9hwaR4k5bBGFaVGAoLC5GSkiIe3+wEzsjIwLZt27Bs2TLU1dVh/fr1qKmpQVxcHA4cOHBHh7QpsCmJiG662SkcHBzcb2IYyLWVDW3YkX8ZVxpurXiqLj2EA5v/HXK5HE/n7MFI/2GS7NsAWFliSE5O7nczieeeew7PPfechSIiIjLO/v37sX37dnR2dmL58uUGj3XpBOw9UYPsH2vFdY1kkGFOlD9yjp6DQqGATqcT921gYiAisnF97dZ2rUMO1aEraNDe+vIb5OOGR5JCEenvCcXlufjzn/7U674NfSUcU2Ni6AX7GIjIWDd3a/v5t/77H1iI3HMN+KbWE8P9O+Ds7AyFTIbU6BGYFz0CzoruMUB97dtg6e1BbWZUkqWtXbsWZWVlKCgokDoUIrIRt+/WNmXaTPzpm5/w9ZlruNlKHurrjufmhCFCfhVdHVqD3+9t3wZLbw/KxEBEZCI3v/WnpaXh/7z5DkoV43GhXtP9oAxIHuuL/zVvLHxdgbKyMrS3tw/ovpbeHpRNSUREJpSceh8Km72h9hwJ5xuroPp5OGNSwHXMm+AHJ7nx38ctvT0oE0Mv2MdARDcNtOP30tXr+OuhSly57gz/G/uAzRgzHCmR7sjLPWv0/X4uLS0Nzs7OSE1NHVJZBoKJoRecx0BEwK2OXwC9dvwKgoBvz9bj8+PVaNN2AgDcnOXImBWJ2FAlGhoa7rifpTqSB4N9DEREfXj33XcNjv/+978bHGu0XXj30AV8eqxK3E3Nz1WH52aHIjb0zi+Vlu5IHgwmBiKiQapsaMMbX53DyapbG3rNHjMMKQHX4dvLXgmW7kgeDCaGXqhUKkRHRyMxMVHqUIhIQk888QQAiAvZrV69GgBQdLEBbx4sx1VNBwDAy9UJT80ZhV/cNRzy29a8c3NzQ3R0NNzc3AxGLm3fvt2gGenn10mJfQy9YB8DEQHdI4I+++wz5ObmIjk5GQseeACfHqtC7pk68ZpwPw+smhUBXw8Xg/6Em9zd3RETEyMe99aRfPt1UmFiICLqx7x58zBmzBiMCAnHf+dU4GxNMzTXNfD08MTd4/zxcHyIOIPZHjAxEBH148MPP8QH/9oF5cR7ERY/Fzq9DtdbW7EsMQz3xYZafL8Ec2NiICLqw549e/D4448DkAEH9uKBX2+Gs1xA04mDuBq4HLK4CKlDNDkmBiKiXgiCgNc2qW4eAQBOfLEVl86cgFwuxzdf7rXKeQhDZT+NYibGUUlEjq1Tp8f7Ry+htsVwoTtFR4vVz0MYKiaGXnB1VSLH1dzWiT9/8xOOXWpETMoSg8ceeeSRPuchWMuQ06FgUxIR0c9UNrThb4cq0Hi9e2mLu6bfi1lv/x2H9n+C5cuXY/ny5YiJiel1QTtrGXI6FEwMREQ3lFU1Y9uRi9DeWBV1mIcznpw9Ch4IxegR3uK8A0suaCcFJgYiIgDfn7+GnQVXoL+xo07EcA88MSsSPu7OaGgY2L4JA2XtzU1MDETk0ARBwIGTahw4pRbPxYUp8ej08F4nrQ31g93am5uYGIjIYXXp9dhVWInvz18Tz80Z748lU4L7nLRm7R/sQ8XE0Atu1ENk37RdOrx3+CJ+rG4Rzy2OC0bKhAAJo7IOHK7aCw5XJbJfLe2d+PM358Wk4CSXIePuiF6TgrX3CZgaawxE5FBqm7XY8u15cblsd2cFVs+KxLhAr15/x96bjm7HxEBEDqOiXoN38i5A09EFoHs46tNzRmOk0vprApastTAxEJFDKK1swvtHLqFT1z1HIXiYO566ZxSG9bLTmrWxZK2FiYGI7N535+rxUVEVhBsL4Y0b4YXVsyLh7qKQODLrxM5nIrILbW1tOHXqFNra2sRzgiDgi+PV2FVUKSaF+AhfPJ08ikmhD0wMRGQX2tvbUVZWhvb27lnKXXo9tv9wGV//WCtek3rXCDw2PQxOcn709YVNSURkF/bv34/t27ejs7MTDy5dhvcOX8Tpmu7hqDLI8ODUYNwz3l/iKG0DEwMR2bw9e/ZgxYoVkMvl2LdvHzLONkE54W4AgJNcjpUzwhEbppQ4StvB+hQR2bycnBxx8xyZXIHjBUcAAB4uCqxNGc2kYCQmhl5wBzci25GSkgKdTgeZXA5Br0PoXQnw9XDBi3PHYnSAp9Th2RyZINxYY5Z61NzcDKVSiaamJvj4+EgdDhH14MSVJvznpr/jTNF3GDN1JlLmL8RTc0ZB6W4bcxQswZjPMvYxEJFNO1x+FbsKKxE25R64h0Zj2l1heHbuGLg7czjqYDExEJFNEgQBe0tr8HXZreGo4R6dyEgayaQwREwMRGRzuvR67CyoRH7FrX0UZo8ZBsWVCijkve+jQAPDxEBENuV6Rxe2fncR52pbAXTPUVgyNRhJYZ44r9Q6zNLY5sTEQEQ2o65Fi7/lVaC2RQvgzjkKjrQ0tjkxMRCRTThfp8Hfv7uAVm33ktlerk54YnYkRvlzOKqpMTEQkdUrutiAHT9cRpe+e3R9kI8bnrxnFIZ7uUgcmX1iYiAiqyUIAr48pcb+k2rxXFSgNzJnRnB1VDNiYiAiq9Sp0+Nf+VdQeLFBPDdjzHA8HB/M1VHNzO7+dZcsWQJfX188/PDDAzpPRNanVduFv+ScF5OCDDIsigvGsoQQJgULsLt/4RdffBEffPDBgM8TkXWpbdbiza/Lcb5eAwBwVsiROSsC904IgEzGOQqWYHeJITk5Gd7e3gM+T0TS2LVrF1auXIldu3aJ586qW7Hp63Oob+0ejurj5owX7h2D2FCujmpJVpUY8vLysHDhQgQHB0Mmk2H37t13XKNSqRAZGQk3NzdMmzYN+fn5lg+UiIZkz549+OUvf4nt27fjl7/8Jfbs2YMfKq5hS+55tHXqAADBw9yRdd9YhA/3kDhax2NViUGj0SA2NhYqlarHx3fu3ImsrCxs2LABxcXFiI2Nxfz581FbW9vj9URknXJyciCXy6HX6yGXy7F1117s+OEydDcWe75rpDdenDsGvh4cjioFq0oMaWlpePXVV7FkyZIeH9+0aRPWrFmDzMxMREdHY8uWLfDw8MDWrVstHCkRDYWHhwf0ej0AQK/Xo7JFLz42e5w/1sweBTcuhCcZq0oMfeno6EBRURFSU1PFc3K5HKmpqTh69KjJnker1aK5udngh4hM6/r165DfGF0kk8nR1dl+Y1/mEDwcH8KF8CRmM4mhvr4eOp0OgYGBBucDAwNRU1MjHqempmLp0qXYt28fQkNDxaTR2/nbbdy4EUqlUvwJCwszX6GIHJSgcOmuMchkEAQ9Iicm4onZkZgz3l/q0Ah2OMHt4MGDRp2/3bp165CVlSUeNzc3MzkQmdCWD3birTdeB2QyQBAwfclq/Om3TyDMj53M1sJmEoO/vz8UCgXUarXBebVajaCgIJM9j6urK1xdXaFSqaBSqaDT6Ux2byJHV3SxAe/s/AIyuQKCvnuP5hh/JyYFK2MzTUkuLi6Ij49Hdna2eE6v1yM7OxszZsww+fOtXbsWZWVlKCgoMPm9iRzRt2fr8Y+jlxF8V4KYFAS9Hvcm3yN1aHQbq6oxtLa2ory8XDyuqKhASUkJ/Pz8EB4ejqysLGRkZCAhIQFJSUnYvHkzNBoNMjMzJYyaiPoiCAL2larxVVl3bX90QjJ+9Ye/oq40DwkJCb2OQiTpWFViKCwsREpKinh8s60/IyMD27Ztw7Jly1BXV4f169ejpqYGcXFxOHDgwB0d0qbApiSiodMLAnYVVuLIT1fFc/dFB2LBsjWQyZ6UMDLqi0wQbswooR41NzdDqVSiqakJPj4+UodDZDP0goB/5l8x2Jf5wakhHHkkEWM+y6yqxkBE9kEvCNjxw2UUXOheHVUhk2HF9DDER/hKHBkNBBMDEZmUTi9g+w+XUXTxVlLImBnBhfBsiM2MSrI0lUqF6OhoJCYmSh0Kkc3QC4ZJwUkuQ+YsJgVbwz6GfrCPgWhgBEHAx8VVOHSuHkB3Unh8ZgQmhTApWAP2MRCRxX11qlZMCnKZDJkzIzExhF+mbBGbkohoyA6XX8W+k7fWLHskKZRJwYYxMfSCfQxEA3P8ShN2FVaKx4vigpE0yk/CiGio2MfQD/YxEPXu0rXr+FP2T+jUde+nkDIhAIvjgiWOinpizGcZawxENChNbZ34+6ELYlJIjPTFotiREkdFpsDEQERG69TpsfW7C2hs6wQAjPL3xLLEUMhk3GDHHjAx9IJ9DES9+7i4CheuXgcADPNwxupZEXBW8OPEXrCPoR/sYyAyVHChAf/z/SUAgLNCjl+ljkWor7vEUVF/2MdARGZR3dSODwuuiMdLE0KYFOwQEwMRDYi2S4f3Dl9Ex43O5umj/TCNw1LtEhMDEQ3IZyXVUDe3AwBChrnj4fgQiSMic2FiIKJ+naxsxuHy7s12XBRyPH43O5vtGV/ZXnBUElG3lvZO/DP/sni8eEowRvi4ShgRmRtHJfWDo5LIkQmCgL9/dwGllc0AgEkhPlg9K5LzFWwQRyURkUmUXG4Sk4KXqxOWJYYxKTgAJgYi6pFG24WPim4tjrc0IQTeblyp3xEwMRBRjz49VoVWbRcAYHKokruwORAmBiK6Q1l1MwoudG/P6e6swMPxIWxCciBMDERkoL1Thw8LbjUhLZ4SDKW7s4QRkaUxMfSCw1XJUe0/qUbD9Q4AwLgRXpg2ylfiiMjSOFy1HxyuSo6kqrEN//XlOegFAc4KOX77i/EI8OacBXvA4apEZDRBEPBxcRX0N74rzrtrBJOCg2JiICIAwLHLTSivbQUA+Hu54N67AiSOiKTCxEBE0Hbp8NmxKvF4yZQQroXkwPjKExG+OlUrbtMZPdIHMcHeEkdEUmJiIHJwtc1a5J6pAwA4yWV4cGow5yw4OCYGIge353gVuvTdHc4pEwLY4UxMDESOrLy2VVwkT+nujHnRIySOiKwBEwORgxIEAZ+VVIvHCyYFwdVJIWFEZC2YGIgcVPGlRly6dh0AEDzMHUmc4Uw3MDH0gktikD3r1OnxxYka8Tg9diTk7HCmG7gkRj+4JAbZo29O1+Gzku55CxOCvPFM8miJIyJz45IYRNQrjbYLX51SAwBkkGFR3EiJIyJrY9R2THv27DH6CebNmwd3d3ejf4+IzOOrslq0deoAAEmjfBE8jO9PMmRUYli8eLFRN5fJZDh37hxGj2Y1lcga1LVo8d25egCAs0KOBZOCJI6IrJHRTUk1NTXQ6/UD+vHw8DBHzEQ0SHtP1NyazBYVgGEe3ICH7mRUYsjIyDCqWejRRx9lhy2Rlaio1+DY5UYAgJerE+Zy9VTqBUcl9YOjksgeCIKAt7J/QkW9BgDwy4RQzBw7XOKoyJI4KomIDJy40iwmhUAfN0wf7SdxRGTNmBiI7FyXXo89x28tfZEeOxIKOSezUe9Mmhg0Gg2OHTuG+vp6U96WiIbgcPk11LdqAQBjR3hxrwXql8kSw2uvvYYHH3wQn3zyCZ599lmsWrUKbW1tprr9gCxZsgS+vr54+OGHDc5/8cUXiIqKwrhx4/Duu+9aNCYiKbV16PDlSbV4vChuJPdaoH4ZNY+hN++99x6uXbuGL7/8Ujy3a9cu/OpXv8Jf//pXUzzFgLz44otYtWoV3n//ffFcV1cXsrKykJOTA6VSifj4eCxZsgTDh7PjjezfV2VqaDq6AAAJEb4I9+MQcuqfSWoMO3fuxEsvvQQAWLVqFdRqNZYuXYqCggJT3H7AkpOT4e1tWE3Oz89HTEwMQkJC4OXlhbS0NHz11VcWjYtICtc0Hcg7292s6ySX4/7JnMxGA2OSxODs7IyOjg4AwMyZM8WJbXq9fsD3yMvLw8KFCxEc3L2t4O7du++4RqVSITIyEm5ubpg2bRry8/P7vW9VVRVCQkLE45CQEFRWVg44LiJb9fnxanEyW3KUP/w8XSSOiGyFSRLDE088gd/85jfQ6XRYvXo1vL298eabbyItLW3A99BoNIiNjYVKperx8Z07dyIrKwsbNmxAcXExYmNjMX/+fNTW1pqiCER25eLV6yi+1AigezJb6l3cmY0GziR9DIsWLYJGo8HcuXMRERGB6upq3H333Xj11VcHfI+0tLQ+E8mmTZuwZs0aZGZmAgC2bNmCvXv3YuvWrWIzVk+Cg4MNagiVlZVISkrq9XqtVgutViseNzc3D7gMRNZAEATsPlYlHv9iYiDcXbgzGw2cSRIDACxfvhzLly/HtWvX4Ovra9KRDx0dHSgqKsK6devEc3K5HKmpqTh69Gifv5uUlISTJ0+isrISSqUS+/fvx8svv9zr9Rs3bsQrr7xistiJLO3ElWacvzGZbYS3K2aM4WQ2Mo7JJ7j5+fmZfDhcfX09dDodAgMDDc4HBgaipubWLlSpqalYunQp9u3bh9DQUBw9ehROTk544403kJKSgri4OPz617/uc0TSunXr0NTUJP5cvnzZpGUhMqc7JrPFjYSTnPNYyTgmqzFYg4MHD/Z4Pj09Henp6QO6h6urK1xdXU0ZFpHFGExmC/DCxGCu70XGs4mvEv7+/lAoFFCr1Qbn1Wo1goLMMwSPez6Trbne0cXJbGQSRtcYVq1aNaDrtm7danQwvXFxcUF8fDyys7PFzYL0ej2ys7Px3HPPmex5fm7t2rVYu3atuCIhkbX76lStOJktMdIX4cM5mY0Gx+jEsG3bNkRERGDKlCkw5Yrdra2tKC8vF48rKipQUlICPz8/hIeHIysrCxkZGUhISEBSUhI2b94MjUYjjlIicmQ1Te3iZDZnhRz3c2c2GgKjE8MzzzyDf/7zn6ioqEBmZiYeffRR+PkNfdRDYWEhUlJSxOOsrCwA3ZsDbdu2DcuWLUNdXR3Wr1+PmpoaxMXF4cCBA3d0SJuKSqWCSqWCTqczy/2JTEUQBHxSXAXdjS9qcycEwJeT2WgIBrVRj1arxSeffIKtW7fiyJEjuP/++7F69Wrcd999dtemyY16yNodv9KErd9dAAD4ebrgdwui4Kywie5DsiCzb9Tj6uqKRx55BF9//TXKysoQExODZ599FpGRkWhtbR1U0ERkvE6d3mAy2+K4YCYFGrIh/wXJ5XLIZDIIgsBmFyIL++bHOlzTdK9TFhXojcmhrNXS0A0qMWi1Wvzzn//EvHnzMH78eJSWluLtt9/GpUuX4OXlZeoYJcHhqmTtrmk68PWP3WuFKWQyPDg12O6ackkaRvcxPPvss/jXv/6FsLAwrFq1CitWrIC/v7+54pMc+xjIGgmCgHcOXcCpqu61vJKjArBkSrDEUZE1M+azzOjEIJfLER4ejilTpvT57eSTTz4x5rZWi4mBrNGxS43YduQiAEDp7ox1aVFcKI/6ZMxnmdHDVVeuXOkQ1VUOVyVrdb2jCx8X31ox+KGpIUwKZFKDGq7qSFhjIGvzr4IrOPrTVQDApBAfrJ4V6RBf1mhozD5clYikUV7bKiYFNycFHooPYVIgkzMqMZw4ccKo7TpPnTqFrq4uo4Miojt16vTYWXBFPH4gNgi+HpzhTKZnVGKYMmUKrl69OuDrZ8yYgUuXLhkdFBHdaV9pDWpbupfUjhzugZlje99XhGgojOp8FgQBL7/8Mjw8BrZqY0dHx6CCsgbsfCZrUl7bipzT3YvkOcllWJYYCjmbkMhMjOp8Tk5ONro9c8eOHRg5cqTRgVkLdj6T1No6dPi/B86i4Xr3F61FccG4d0KAxFGRrTHbcNXc3NyhxEVEg/DJsSoxKYwN8EJylP1OKCXrwFFJRFbs+OUm5FdcA9A9CmnF9DA2IZHZMTEQWamG6x0Go5Aeig+GH/dZIAtgYiCyQl16PbYduSRu1RkbqkRipK/EUZGjYGLoBVdXJSl9frwGF+o1AABfDxcsSwzlRDaymCEtidHZ2Ymamhpcv34dAQEBJtni09pwVBJZWsnlRrx3uHuBPCe5DC/OHYvw4QMbIk7UG7MuidHS0oL//u//xpw5c+Dj44PIyEjcddddCAgIQEREBNasWYOCgoJBB0/kyGqbtfjnD7f6FZZMCWZSIIszKjFs2rQJkZGReO+995Camordu3ejpKQEZ8+exdGjR7FhwwZ0dXXhvvvuwy9+8QucO3fOXHET2Z22Th22Hr6A9q7uSZUJEb6c3UySMGoeQ0FBAfLy8hATE9Pj40lJSVi1ahW2bNmC9957D4cOHcK4ceNMEiiRPdPpBbx/5CKqm9oBAEE+bvhlIhfII2kMqo9Bp9Ph888/x9y5c+Ht7W2OuKwG+xjIEj4urkTe2e4lLzxcFPhfqeMwwsdV4qjInph92W2FQoFHHnkEdXV1gwqQiG7JO1svJgWFTIZVMyOZFEhSgx6umpiYiIqKClPGYlU4XJUsoay6GZ8UV4nHyxJDMS7QS8KIiIYwXHXXrl344x//iI8++ghhYWGmjstqsCmJzOVCvQZ/yT0PbVf3Hiepd43AwljbXXCSrJtZ93y+admyZQCAmJgYpKenIzk5GVOmTMGkSZPg4sJp+0R9qWpsw1/zKsSkEBuqxAOTgySOiqjboBNDRUUFjh8/jpKSEhw/fhwbN27EhQsX4OTkhKioKJw4ccKUcRLZjfpWLf47twLXO7qHpY4b4YXHZoRzBBJZjUEnhoiICERERCA9PV0819LSgpKSEiYFol40tXXiLznn0dzeCQAI9/PAE7Mj4azg6jRkPYa0JIYjYB8DmUpTWydUOeehbr41V+GFuWPg6Tro72dEA2aRPgYiGrimtk68/c1P4p7Nfp4ueCZ5NJMCWSXWX4nMrKek8FzKGAzzcJY4MqKeDSoxaDQalJaW9vjYqVOn0NraOqSgiOxFT0nh+XvHYLgXR+6R9RpUYujs7MS0adOQn59vcL6srAxTpkxhYiACUNeixVsHy+9ICtyFjazdoBLDsGHD8MADD+CDDz4wOP+Pf/wDc+fORVAQx2OTY6tqbMOfsn/CVU0HAGA4kwLZkEH3MWRkZGDnzp3o6ureelAQBGzfvh2ZmZkmC05KXBKDBut8nQZ/yv5JHJIapHTDC3PHMimQzRj0cFWdTofQ0FBs2bIFixYtQk5ODh566CHU1NTY1cxnDlclY5RVN2PrdxfRqeue0Rw53ANP3jOKo49IcmZfXRXoXmF1xYoVYnPSP/7xDyxbtsyukgKRMYovNeLvhy6ISSEq0BvPpnBIKtmeIf3FZmRkICkpCZWVlfj444/x5ZdfmiouIptyuPwqdhVWQkB3BTw2VInHZoRzRjPZpCHPfI6Pj4e3tzdqampw+vRpU8VlNdiURH0RBAFfnlJj/0m1eG76aD8sSwyFnGsfkRWxSFPSTStXrkReXh5Wrlw51FsR2RS9IODDwkqDpDD3rhH4NyYFsnFDbvx87LHH0NjYiFWrVpkiHiKb0KnT4/0jl1Ba2SSeWxwXjJQJARJGRWQaXESvH2xKottptF1459AFVNRrAABOchmWTwtDfISvxJER9Y6L6BGZScP1DmzJrUDNjRVSXZ3kWD0rElFB3hJHRmQ6TAxEA1Td1I4t355H4/XuiWterk54as4ohPt5SBwZkWkxMRANwPk6Df6WV4G2zu5d1/y9XPH0nFEI8HaVODIi03OIQdZ//OMfERMTg4kTJ+J//ud/pA6HbMyJK034S+55MSmE+XngV6ljmRTIbhmVGC5dumTUzSsrK4263hxKS0uxY8cOFBUVoaCgAG+//TYaGxulDotsxOHyqwZLXEQFeuO5lNHwdmNlm+yXUYkhMTERTz31FAoKCnq9pqmpCe+88w4mTpyIjz/+eMgBDtWPP/6IGTNmwM3NDe7u7oiNjcWBAwekDousnCAI+OqUGh8WXhFnMydE+OLJOZFwc1ZIHB2ReRmVGMrKyuDp6Yl58+YhKCgI999/P9asWYPnn38ejz76KKZOnYoRI0Zg69ateP311/HCCy8MOcC8vDwsXLgQwcHBkMlk2L179x3XqFQqREZGws3N7Y59IiZOnIjc3Fw0NjaioaEBubm5VlGTIeslCAI+P1GDvaU14rmUCQF4dHoYnOQO0fpKDs6ov/Lhw4dj06ZNqK6uxttvv41x48ahvr4e586dAwCsWLECRUVFOHr0KBYsWGCSADUaDWJjY6FSqXp8fOfOncjKysKGDRtQXFyM2NhYzJ8/H7W1tQCA6OhovPDCC7j33nvx4IMPYvr06VAo+I2PeiYIAnYVVSL7x1rx3MLJI7E4rvuLCZEjMHqC26uvvor4+HhMnToVgYGB5oqrRzKZDJ9++ikWL14snps2bRoSExPx9ttvAwD0ej3CwsLw/PPP46WXXrrjHk888QSWLFmC+++/v8fn0Gq10Gq14nFzczPCwsI4wc0BdOn12PHDFRRdbAAAyCDD0oQQzBw7XOLIiIbOrBPc1q9fL35zCgoKwtSpUxEfHy8mi5CQkMFFPQgdHR0oKirCunXrxHNyuRypqak4evSoeK62thYjRozAmTNnkJ+fjy1btvR6z40bN+KVV14xa9xkeW1tbTh//jxGjx4Nd3f3Ox7vXuLiIkormwEAcpkMj07nbGZyTEYnhsTERFRXVyMzMxP+/v4oLi7GJ598gtdeew06nQ4BAQGYOnUq9u3bZ454DdTX10On091RcwkMDDRY6XXRokVoamqCp6cn3nvvPTg59V7sdevWISsrSzy+WWMg29be3o6ysjIEBwffkRg6dXq8e+gCTte0AACc5HJkzozAxBDWEMkxGZ0YfvjhB2zbtg2/+93vkJiYiE2bNmHMmDHQarUoKSlBcXExjh07Zo5YB+3ntYf+uLq6wtWV49Mdxe1JwdVJjidmj8L4QC+JIyOSzqCGWDz++OM4e/YsoqKiMHXqVKxbtw46nQ7Tpk3DM888g7/97W+mjrNH/v7+UCgUUKvVBufVajWCgoKGdG/u+Wz/ekoKT88ZzaRADm/QY++8vLzw+uuvo7CwECdPnsTYsWPFbT4txcXFBfHx8cjOzhbP6fV6ZGdnY8aMGUO699q1a1FWVtbnnA2yHfv378e2bduwf/9+AL0nhdEBnlKGSWQVhjR9s6urC1qtFo888gjUajUyMzPxwAMPwM/Pz1TxobW1FeXl5eJxRUUFSkpK4Ofnh/DwcGRlZSEjIwMJCQlISkrC5s2bodFokJmZabIYyLbt2bMHK1asgFwux759++Dm4QH1sMlMCkS9MDox/OEPf0BpaSlKS0tx+vRpuLm5YfLkyUhKSsJTTz0FpVJp0gALCwuRkpIiHt/sGM7IyMC2bduwbNky1NXVYf369aipqUFcXBwOHDgw5KG0KpUKKpUKOp1uSPch6eXk5EChUECn00GhUOAvOz5HzIOjADApEPXE6HkMcrkckZGRyMjIwCOPPILx48ebKzarwI16bN+ePXuwaNEiyOVy6PV6PPDrzRidkMykQA7FmM8yoxPDnDlzUFJSgpaWFnh6emLy5MmYOnWq+DNx4kS7mlnMxGAftm/fjrfe+xDDJ92D8dNS4SSX4+k5ozCOHc3kIMw6we3bb78FAJw7dw5FRUUoLi5GcXExtm/fjsbGRri6umLSpEkG6xXZIjYl2ReXMdMRuWA4/AP8IZfJ8PjMcCYFol6YdM/niooKFBYW4tixY3jttddMdVtJscZg+3LP1OHD/Iuor6uHf4A/MmePQWIkZzSTY5Fsz+dRo0Zh1KhRWLp0qSlvSzRo+RXX8OmxKvH4/hh/JgWifnANYbJbZ2pa8K/8K+JxtFKLu0eZdtQckT3iNlS9YB+DbatqbMPWwxehu9FSOnucPyZ5ucHNzU3iyIisn0n7GOwR+xhsT1NbJzZ9fQ6N1zsBAJNCfLBqViTk3E+BHJgxn2VsSiK7ou3S4W95FWJSCPfzwGMzwpkUiIzAxEB2QxAEfHDkEq40tAEA/Dxd8OQ9kXB1sp95NUSWwMRAdmNvaQ1OVnVvtOPurMBT94yCt5uzxFER2R4mhl5w2W3bUnypEV+Xde/TLIMMmTMjEKRkRzPRYLDzuR/sfLZ+l69dx1vZP6FTpwcAPDglGHOiAiSOisi6sPOZHEZLeyfe/e6CmBSmj/bDPeP9JY6KyLYxMZDN0ukFbDt8SRyBNMrfEw/Hh0DGEUhEQ8LEQDZrb2kNyutaAQDD3J2xamYEnBX8kyYaKr6LyCYdv9KE7B+7O5sVMhkenxkBH3eOQCIyBSaGXnBUkvWqa9Fixw+XxeNFU4Ixyp+b7RCZCkcl9YOjkqxLR5cebx4sR1Vj9yS2KeHDkDEjnP0KRP3gqCSyW58cqxKTQqCPGx5JCmVSIDIxJgayGSWXG3H0p6sAABeFHKtmRnC5CyIzYGIgm9Cg6TDYW+Gh+BDObCYyEyYGsno6vYAPjl5CW2f33hhTwoZh2ijuwkZkLkwMZPW+OqXG+XoNgO4VU5clsl+ByJyYGHrB4arW4afaVnx5qnu+glwmw2PTw+Huwn4FInPicNV+cLiqdDTaLrz+5VlxyYsFk4IwPyZQ4qiIbBOHq5Jd2FVUKSaFsQFemBc9QuKIiBwDEwNZpWOXGnHsUiMAwMNFgUdnhHF7TiILYWIgq9PS3oldhZXi8cPxIfD1cJEwIiLHwsRAVkUQBOwsuAJNRxcAIDZUianhw6QNisjBMDGQVSm40IDSyu59m71cnbA0gUNTiSyNiYGsRsP1DnxSXCUeL0sMhbebk4QRETkmJgayCoIg4F/5V8TZzYmRvpgcqpQ4KiLHxMRAVuHIT9dwuqYFAKB0d8aDU4MljojIcTEx9IIzny2n8Xon9pRUi8f/lhQKDxc2IRFJhTOf+8GZz+YlCAL+/t0FscN52ig/LJ8WJnFURPaHM5/JZpy40iwmBW83JyyKGylxRETExECSaevQ4aOiWxPZHpoaAk9XNiERSY2JgSTz2fFqNLd3r4UUE+yDuDCOQiKyBkwMJIny2lZxm05XJzmWJoRwIhuRlWBiIIvr1Omxs+DWNp0PTB7JtZCIrAgTA1ncl6fUqG3RAgAi/T0xa9xwiSMiop9jYiCLqmlqR87pOgCAk1yGf0sM5XLaRFaGiYEsRhAEfFRUiS5999SZlAkBGKl0kzgqIrodEwNZTPGlRpyrbQUA+Hm64L5obtNJZI2YGMgi2jp1+Oxny148NDUELk788yOyRg7xznzzzTcRExOD6OhovPDCC+AqIJZ34KQaTW3dcxYmBvtgYgiXFyGyVnafGOrq6vD222+jqKgIpaWlKCoqwvfffy91WA6lqrENeWfrAQDOCjkenBoicURE1BeHWH+gq6sL7e3tAIDOzk6MGDFC4ogchyAI+LCwEvobtbT7okdguBfnLBBZM6uvMeTl5WHhwoUIDg6GTCbD7t2777hGpVIhMjISbm5umDZtGvLz88XHAgIC8O///u8IDw9HcHAwUlNTMWbMGAuWwLHlX2hARb0GADDC2xUpEwIkjoiI+mP1iUGj0SA2NhYqlarHx3fu3ImsrCxs2LABxcXFiI2Nxfz581FbWwsAaGhowBdffIELFy6gsrISR44cQV5eniWL4LDaOnQG+yw8HB8CZ4XV/8kROTyrf5empaXh1VdfxZIlS3p8fNOmTVizZg0yMzMRHR2NLVu2wMPDA1u3bgUAHDx4EGPHjoWfnx/c3d1x//3399nHoNVq0dzcbPBDg3PglBqt2i4AwJSwYYgK8pY4IiIaCKtPDH3p6OhAUVERUlNTxXNyuRypqak4evQoACAsLAxHjhxBe3s7dDodcnNzERUV1es9N27cCKVSKf6EhXHTmMGobdbi0M86nLnPApHtsOnEUF9fD51Oh8BAw4lSgYGBqKmpAQBMnz4dCxYswJQpUzB58mSMGTMG6enpvd5z3bp1aGpqEn8uX75s1jLYq0+PVUF3o8N57oQA+Hqyw5nIVjjEqKTf//73+P3vfz+ga11dXeHq6mrmiOxbWXUzyqq7m+CGeThj7l0cBUZkS2y6xuDv7w+FQgG1Wm1wXq1WIygoaEj3VqlUiI6ORmJi4pDu42i69Hp8WlwlHqfHjuQMZyIbY9PvWBcXF8THxyM7O1s8p9frkZ2djRkzZgzp3mvXrkVZWRkKCgqGGqZD+e7cVXFJ7VH+npgaPkzagIjIaFbflNTa2ory8nLxuKKiAiUlJfDz80N4eDiysrKQkZGBhIQEJCUlYfPmzdBoNMjMzJQwasfU0t6FAydv1d4enBLMXdmIbJDVJ4bCwkKkpKSIx1lZWQCAjIwMbNu2DcuWLUNdXR3Wr1+PmpoaxMXF4cCBA3d0SBtLpVJBpVJBp9MN6T6OZP/JGrR1dv97TRvlh/DhHhJHRESDIRO4olyfmpuboVQq0dTUBB8fLvzWm6rGNrx+4BwECHB1kuN/3z8BPu7OUodFRDcY81lm030MZB0EQcAnxVUQcHM9pEAmBSIbxsTQC45KGrgTlc3iBjz+Xq6YE+UvcURENBRsSuoHm5L61qXX4w/7z6Luxkik1bMiMTlUKXFURHQ7NiWRxRwpvyYmhbEBXpjEDXiIbB4TAw1aW4cOX566NTx1UdxIDk8lsgNMDL1gH0P/sn+sFVdPnRo+jMNTiewE+xj6wT6GnjVc78Dv955Bp04PJ7kMv1swgTuzEVkx9jGQ2e0rVaNTpwcAzB7nz6RAZEeYGMholQ1tKKhoAAC4OyswL5qrpxLZEyYGMtpnJdW3JrPFBMLT1epXViEiIzAx9IKdzz37sboFZ9QtAAA/TxfMHjdc4oiIyNTY+dwPdj7fohcE/NeX51DV2AYAWDkjHPERvhJHRUQDwc5nMouCigYxKYT7eXCvBSI7xcRAA9Kp02NfaY14nB7LyWxE9oqJgQbk0LmraGzrBADEBPtgXKCXxBERkbkwMVC/2jp0OFhWCwCQQYYHJg9tP20ism5MDL1whFFJbW1tOHXqFNra2vq8LudMHTQd3UtfJEQOQ/Awd0uER0QSYWLoxdq1a1FWVoaCggKpQzGb9vZ2lJWVob29vddrWto7kXumDgCgkMnwi4lD2zKViKwfEwP16auyWmi7upe+uHvscPh7uUocERGZGxMD9eqapgNHyq8CAFwUctzHpS+IHAITA/Vq/0k1uvTd8x/njPfnPs5EDoKJwYHt378f27Ztw/79++94rKapXVwoz8NFgbl3sbZA5Ci4+pmD2rNnD1asWAG5XI59+/bBy8sL6enp4uN7S2vEhfLmThgBdxeFVKESkYWxxtALex+umpOTA4VCAb1eD4VCgdzcXPGxC/UanLjSBABQujvjnvH+EkVJRFJgYuiFvQ9XTUlJgU6ng1wuh06nQ3JyMgBAEAR8ceLW0hfzYwLh4sQ/EyJHwne8g0pPT8f27duRlpaG7du3i81IZ9WtOFfbCgDw93LFtNFcPZXI0bCPwYGlpaXB2dkZqampAO6sLSyYFAgnOb87EDkavutJdOJKMy5duw4ACBnmzmW1iRwUEwMB6N6E5+fLaj8wOYjLahM5KCYGAgAUX2xETXP3mkmj/D1x10hviSMiIqmwj4Gg0ws4cKpWPF4wibUFIkfGGgOhpLIFdS1aAMDYEV4Yz014iBwaE4OD0wvAN2cbxOP7J3ETHiJHx8TQC3uf+QwAbm5u0A8fgxZt99IXE4K8MTrAU+KoiEhqTAy9sPeZzwDg5OKK0xo3yBXdfwYLWFsgIjAxOLSjP11D4/VOAEBMsA8ihntIHBERWQMmBgfVqdPj4I+GI5GIiAAmBof13bmraGrrri3EhioR6usucUREZC2YGByQtksn1hZkkOEXEwMljoiIrAkTgwM6dPYqWrVdAIC4cCWCh7G2QES3MDE4mLZOHb45XQegu7aQxtoCEd2GicHB5J2ph6aju7aQEDkMgT5uEkdERNaGicGBXO/oQs6Z7tqCXCbD/BjWFojoTkwMDiT3TD3aOnUAgKRRvgjwdpU4IiKyRkwMDuJ6Rxe+PVsPAFCwtkBEfXCIxHDmzBnExcWJP+7u7ti9e7fUYVnUt2fq0f6z2oKfp4vEERGRtXKI/RiioqJQUlICAGhtbUVkZCTmzZsnbVAW1NahM6gtzItmbYGIeucQNYaf27NnD+bOnQtPT8dZRTTv7K2+hcRRvhjuxdoCEfXOJhJDXl4eFi5ciODgYMhksh6bgVQqFSIjI+Hm5oZp06YhPz+/x3t9+OGHWLZsmZkjth5tHTqDkUjzokdIHBERWTubSAwajQaxsbFQqVQ9Pr5z505kZWVhw4YNKC4uRmxsLObPn4/a2lqD65qbm3HkyBEsWLDAEmFbBYPaQqQv/L04EomI+mYTfQxpaWlIS0vr9fFNmzZhzZo1yMzMBABs2bIFe/fuxdatW/HSSy+J13322We477774ObW+6QurVYLrVYrHjc3N5ugBNJo69Qh9yxrC0RkHJuoMfSlo6MDRUVFSE1NFc/J5XKkpqbi6NGjBtcOpBlp48aNUCqV4k9YWJhZ4raEQ2frcb2ju7aQEMl5C0Q0MDafGOrr66HT6RAYaDjSJjAwEDU1NeJxU1MT8vPzMX/+/D7vt27dOjQ1NYk/ly9fNkvc5tbWadi3cB9rC0Q0QDbRlGQKSqUSarW63+tcXV3h6mr736x/XluIjxjG2gIRDZjN1xj8/f2hUCju+NBXq9UIChr8rmQqlQrR0dFITEwcaogW196pQ+6Z7nkLMnCWMxEZx+YTg4uLC+Lj45GdnS2e0+v1yM7OxowZMwZ937Vr16KsrAwFBQWmCNOivjt31WAFVdYWiMgYNtGU1NraivLycvG4oqICJSUl8PPzQ3h4OLKyspCRkYGEhAQkJSVh8+bN0Gg04iglR6LtMtxv4T7OciYiI9lEYigsLERKSop4nJWVBQDIyMjAtm3bsGzZMtTV1WH9+vWoqalBXFwcDhw4cEeHtDFUKhVUKhV0Ot2Q47ekQ2dv1RbiI4ZhhA9rC0RkHJkgCILUQViz5uZmKJVKNDU1wcfHR+pw+qTt0uH/+/w0WrVdkEGGdQvGcyMeIgJg3GeZTdQYpHQzb9rCRLfcM3W42tAIAIgLU8IdHWhu7pA2KCKyCjc/wwZSF2Bi6MXNpqSbs6BtcaLbM1IHQERWp6WlBUqlss9r2JTUj8bGRvj6+uLSpUv9/mP2JTExcUAjnHq7bqDn+zq+/f+zs7MRFhaGy5cvD6mZbCBl6+uanh7r71xv5bz53+bm5iGXbaiv2UDK0dO5/v7flspmrX+PfV1nr+81QRDQ0tKC4OBgyOV9D0hljaEfN/8BlUrlkF5QhUIxoN/v7bqBnu/ruLf/9/HxMXvZ+rqmp8f6O9dbOW8/P5SyDfU16+0xY8rWV5ltoWzW+vfY13X2/F4b6Jdbm5/HYCvWrl07pOsGer6v497+f6gGcq++runpsf7O9VZOS5erv+uGWjYpX7P+rhtI2az177Gv6+z5vTZQbErqhy2NSjIWy2ab7LVs9louwPbKxhpDP1xdXbFhwwa7WD/pdiybbbLXstlruQDbKxtrDEREZIA1BiIiMsDEQEREBpgYiIjIABMDEREZYGIwocbGRiQkJCAuLg4TJ07EO++8I3VIJnP58mUkJycjOjoakydPxq5du6QOyWSWLFkCX19fPPzww1KHMmRffPEFoqKiMG7cOLz77rtSh2NS9vQ6/Zw1vrc4KsmEdDodtFotPDw8oNFoMHHiRBQWFmL48OFShzZk1dXVUKvViIuLQ01NDeLj43H27Fl4enpKHdqQ5ebmoqWlBe+//z4++ugjqcMZtK6uLkRHRyMnJwdKpRLx8fE4cuSIXfz9AfbzOt3OGt9brDGYkEKhgIeHBwBAq9VCEIQBrWRoC0aOHIm4uDgAQFBQEPz9/XHt2jVpgzKR5ORkeHt7Sx3GkOXn5yMmJgYhISHw8vJCWloavvrqK6nDMhl7eZ1uZ43vLYdKDHl5eVi4cCGCg4Mhk8mwe/fuO65RqVSIjIyEm5sbpk2bhvz8fKOeo7GxEbGxsQgNDcV//Md/wN/f30TR980SZbupqKgIOp3OIivOWrJcUhtqWauqqhASEiIeh4SEoLKy0hKh98ueX0dTls2S762+OFRi0Gg0iI2NhUql6vHxnTt3IisrCxs2bEBxcTFiY2Mxf/581NbWitfc7D+4/aeqqgoAMGzYMBw/fhwVFRXYsWMH1Gq13ZQNAK5du4aVK1fib3/7m9nLBFiuXNbAFGW1Vixb/2Wz9HurT4KDAiB8+umnBueSkpKEtWvXisc6nU4IDg4WNm7cOKjneOaZZ4Rdu3YNJcxBMVfZ2tvbhdmzZwsffPCBqUI1ijlfs5ycHOGhhx4yRZgmMZiyHj58WFi8eLH4+Isvvihs377dIvEaYyivo7W9TrcbbNmkfm/dzqFqDH3p6OhAUVERUlNTxXNyuRypqak4evTogO6hVqvR0tICAGhqakJeXh6ioqLMEq8xTFE2QRDw+OOP495778Vjjz1mrlCNYopy2YqBlDUpKQknT55EZWUlWltbsX//fsyfP1+qkAfMnl/HgZTNGt9bTAw31NfXQ6fTITAw0OB8YGAgampqBnSPixcvYvbs2YiNjcXs2bPx/PPPY9KkSeYI1yimKNvhw4exc+dO7N69G3FxcYiLi0Npaak5wh0wU5QLAFJTU7F06VLs27cPoaGhVvlhNJCyOjk54Y033kBKSgri4uLw61//2iZGJA30dbSF1+l2AymbNb63uFGPCSUlJaGkpETqMMxi1qxZ0Ov1UodhFgcPHpQ6BJNJT09Henq61GGYhT29Tj9nje8t1hhu8Pf3h0KhuKOzWK1WIygoSKKoTMNey2av5eqJPZeVZbM+TAw3uLi4ID4+HtnZ2eI5vV6P7OxszJgxQ8LIhs5ey2av5eqJPZeVZbM+DtWU1NraivLycvG4oqICJSUl8PPzQ3h4OLKyspCRkYGEhAQkJSVh8+bN0Gg0yMzMlDDqgbHXstlruXpiz2Vl2WysbFIPi7KknJwcAcAdPxkZGeI1f/7zn4Xw8HDBxcVFSEpKEr7//nvpAjaCvZbNXsvVE3suK8tmW2XjWklERGSAfQxERGSAiYGIiAwwMRARkQEmBiIiMsDEQEREBpgYiIjIABMDEREZYGIgIiIDTAxERGSAiYGIiAwwMRARkQEmBiIzSk5Ohkwmg0wmEzdxSk5Oxq9+9SuzPefjjz8uPufu3bvN9jxkv5gYiMxszZo1qK6uxsSJEy3yfG+99Raqq6st8lxknxxqPwYic+rq6oKT051vKQ8PD4vu1qVUKqFUKi32fGR/WGMgGoQLFy5AJpPhww8/xOzZs+Hq6oo9e/YM6l579+6FUqnE9u3bAXQ3NT333HN47rnnoFQq4e/vj5dffhk/XyFfr9fj9ddfx9ixY+Hq6orw8HD8/ve/N0nZiFhjIBqE48ePAwD+67/+C6+99hpGjRqFgIAAo++zY8cOPP3009ixYwceeOAB8fz777+P1atXIz8/H4WFhXjyyScRHh6ONWvWAADWrVuHd955B2+++SZmzZqF6upqnD592jSFI4fHxEA0CCUlJfD09MSuXbsQGRk5qHuoVCr853/+Jz7//HPMmTPH4LGwsDC8+eabkMlkiIqKQmlpKd58802sWbMGLS0teOutt/D2228jIyMDADBmzBjMmjVrqMUiAsDEQDQox48fR3p6+qCTwkcffYTa2locPnwYiYmJdzw+ffp0yGQy8XjGjBl44403oNPp8OOPP0Kr1WLu3LmDDZ+oT+xjIBqEkpISJCcnD/r3p0yZgoCAAGzduhXG7q7r7u4+6OclGggmBiIjNTc348KFC5gyZcqg7zFmzBjk5OTgs88+w/PPP3/H4z/88IPB8ffff49x48ZBoVBg3LhxcHd3R3Z29qCfn6gvbEoiMtLx48ehUCgwadKkId1n/PjxyMnJQXJyMpycnLB582bxsUuXLiErKwtPPfUUiouL8ec//xlvvPEGAMDNzQ2//e1v8Zvf/AYuLi6YOXMm6urqcOrUKaxevXpIMREBTAxERjt+/DiioqLg5uY25HtFRUXhm2++QXJyMhQKhfjhv3LlSrS1tSEpKQkKhQIvvvginnzySfH3Xn75ZTg5OWH9+vWoqqrCyJEj8fTTTw85HiIAkAnGNnAS0YAlJycjLi7OoDZgjt/piUwmw6efforFixcP6T7keNjHQGRmf/nLX+Dl5YXS0lKLPN/TTz8NLy8vizwX2SfWGIjMqLKyEm1tbQCA8PBwuLi49Ps7Q60x1NbWorm5GQAwcuRIeHp6Duo+5LiYGIiIyACbkoiIyAATAxERGWBiICIiA0wMRERkgImBiIgMMDEQEZEBJgYiIjLAxEBERAaYGIiIyAATAxERGWBiICIiA/8P/b1ndegymygAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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(<r)$ [M$_\\odot$]\")\n",
    "\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c736e26",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "It looks pretty good already! But let's now use least-squares fitting to\n",
    "optimize our nucleus and halo parameters. We first need to define an error\n",
    "function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6858896e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:31.063275Z",
     "iopub.status.busy": "2025-07-31T14:29:31.062989Z",
     "iopub.status.idle": "2025-07-31T14:29:31.066321Z",
     "shell.execute_reply": "2025-07-31T14:29:31.065834Z"
    }
   },
   "outputs": [],
   "source": [
    "def err_func(p, r, Menc, Menc_err):\n",
    "    pot = get_potential(*p)\n",
    "    xyz = np.zeros((3, len(r)))\n",
    "    xyz[0] = r\n",
    "    model_menc = pot.mass_enclosed(xyz).to(u.Msun).value\n",
    "    return (model_menc - Menc) / Menc_err"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60972d2a",
   "metadata": {},
   "source": [
    "Because the uncertainties are all approximately but not exactly symmetric,\n",
    "we'll take the maximum of the upper and lower uncertainty values and assume\n",
    "that the uncertainties in the mass measurements are Gaussian (a bad but simple\n",
    "assumption):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7ce76299",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:31.067974Z",
     "iopub.status.busy": "2025-07-31T14:29:31.067810Z",
     "iopub.status.idle": "2025-07-31T14:29:31.305290Z",
     "shell.execute_reply": "2025-07-31T14:29:31.304690Z"
    }
   },
   "outputs": [],
   "source": [
    "err = np.max([tbl[\"Menc_err_pos\"], tbl[\"Menc_err_neg\"]], axis=0)\n",
    "p_opt, ier = leastsq(err_func, x0=x0, args=(tbl[\"r\"], tbl[\"Menc\"], err))\n",
    "assert ier in range(1, 4 + 1), \"least-squares fit failed!\"\n",
    "fit_potential = get_potential(*p_opt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "813108e9",
   "metadata": {},
   "source": [
    "Now we have a best-fit potential! Let's plot the enclosed mass of the fit potential over the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "7e0b7cfd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:31.306891Z",
     "iopub.status.busy": "2025-07-31T14:29:31.306725Z",
     "iopub.status.idle": "2025-07-31T14:29:31.851250Z",
     "shell.execute_reply": "2025-07-31T14:29:31.850641Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<>:27: SyntaxWarning: invalid escape sequence '\\o'\n",
      "<>:27: SyntaxWarning: invalid escape sequence '\\o'\n",
      "/tmp/ipykernel_3838/1667721107.py:27: SyntaxWarning: invalid escape sequence '\\o'\n",
      "  ax.set_ylabel(\"$M(<r)$ [M$_\\odot$]\")\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAGGCAYAAAB/gCblAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQYVJREFUeJzt3XtYVPW+P/D3zADDHUEQ5CZeiAIVFAEpTUjMKKVs53ZnJWG5q223Q2efk+ec9Ol5Ks9v753ZydnZrm12sb3Ndl4qtZIkyksIipF4I1GQywAKDIwwwMz6/TG6dJTbwAxrLu/X8/DUWjOs+Sxh5s1a35tMEAQBREREl8mlLoCIiGwLg4GIiEwwGIiIyASDgYiITDAYiIjIBIOBiIhMMBiIiMgEg4GIiEy4SF2ArTMYDKipqYGPjw9kMpnU5RARDYogCGhtbUVoaCjk8r6vCRgM/aipqUFERITUZRARWURVVRXCw8P7fA6DoR8+Pj4AjP+Yvr6+EldDRPaoqakJ+fn5SEtLg7+/vyQ1aDQaREREiJ9pfWEw9OPK7SNfX18GAxENil6vh6enp018jgzkljgbn4mIyASDgYjIynbt2oWNGzdi165dUpcyILyVRERkRTt27MBDDz0EuVyOnTt3wtvbG1lZWVKX1SdeMRARWdHevXuhUChgMBigUCiQn58/rK+v69bjwK8X8FberwP+Hl4xEBFZUXp6OtauXQu5XA69Xo+0tLRhed2GVh0KTjeisKIJHV166C5dGvD3MhiIiKwoKysLmzZtwieffILFixdb9TaSIAg4Xa9FwakG/FLdCgGDW6DT4YJhwYIFyM/Px+zZs/HZZ58BMI5BeOSRR1BfXw8XFxe89NJLWLhwocSVEpGzyMzMhKurKzIyMqxy/C69AYfPNSP/VCNqmttNHnNVyJE4ZgQmjwrGOwM8nsMFw3PPPYelS5figw8+EPe5uLhg7dq1SEhIQF1dHRITE3H33XfDy8tLwkqJiIamvUuPfacvIP9UA1o7uk0eG+HhipnRgUgdHwAvpQs0Gs2Aj+twwZCWlnZD487o0aMxevRoAEBISAgCAwNx8eJFBgMR2aU2XTe+P9mIH043or1Lb/LYmJGemHVTIOIj/ODSz5xIvbGpXkkFBQWYP38+QkNDIZPJsG3bthueo1KpEBUVBXd3d6SkpKCwsNCs1yguLoZer+f8R0Rkd5oudWLrkRq8vOM4vilTi6EggwwJEX74t4wJyJ0TjcQx/oMOBcDGrhi0Wi3i4+OxdOlS3H///Tc8vnnzZuTm5mL9+vVISUnB2rVrMXfuXJw8eRKjRo3q9/gXL17EkiVL8O6771qjfCIiq2hp78I3x9Q4eOYiug1XG5QVMhmSxvpj9s2jMMpXabHXs6lgyMzMRGZmZq+Pr1mzBsuWLUNOTg4AYP369fjqq6+wYcMGvPjii30eW6fT4b777sOLL76IW2+9tc/n6XQ6cduc+3JERJbUpuvGnuP1+PH0BXTpDeJ+V4UcqeMDcEdMEPy93Cz+ujYVDH3p7OxEcXExVqxYIe6Ty+XIyMjAgQMH+vxeQRDw6KOP4o477sAjjzzS53NXr16Nl19+2SI1ExGZo729HWfOnMHo8DE4WNmG7082oqPbeLvIoDegU3cJcydHYM7EUPi4W+/j226CobGxEXq9HsHBwSb7g4ODceLECXE7IyMDR48ehVarRXh4OLZs2QK9Xo/Nmzdj8uTJYrvFRx99hEmTJt3wOitWrEBubq64fWWqWiIia2vVXsLWn8rRXNqJrmuagF3kckyN8IKh6lekj4+1aigAdhQMA7Vnz54e9xsMhh73X0+pVEKpVEKlUkGlUkGv1/f/TUREQyAIAkqqWrClsArlzUoEuurh6iqHQibD9PEBuDM2GIKuDXtqBjdgzVx2EwyBgYFQKBRQq9Um+9VqNUJCQiz+esuXL8fy5cuh0Wjg5+dn8eMTEQHAuQuXsPVIDSoatejqMo5FkAFIivLH3LhgBPkYG5WbdH0cxMJsqrtqX9zc3JCYmIi8vDxxn8FgQF5eHlJTUyWsjIjIfC3tXfjoQCXWfHsaFY1acf8o924svz0cD0+PFENhuNnUFUNbWxvKy8vF7YqKCpSUlCAgIACRkZHIzc1FdnY2pk2bhuTkZKxduxZarVbspWRJvJVERFdcaRQeN24cPDw8hvRcvUHAD6cbsatULTYsA8AoHyVczxXj+6/+hZIRDyF2zGKLn8dA2VQwFBUVIT09Xdy+0gicnZ2NjRs3YtGiRWhoaMDKlStRV1eHhIQE7N69+4YGaUvgrSQiuqKjowNlZWUIDQ3tNxj6eu7Zi+345qcGk/mMvNxccNfEYFwo+xH3L3/MJtZtsKlgSEtLgyD03bjy9NNP4+mnnx6mioiIzLNr1y5s2rQJXV1dWLzY+Fe/tlOPwovu2Lu/Bq6urgCMo5VTxwdg3uQQeCld8G9vfw+FQgG9Xi+u28BgICKyc9ev1ubl5YXwhNvxz4NVOKd1RaCn8Xnh/h5YmBiGqMCr87X1t25DT4FjLQyGXrCNgYjMdWW1tit/9b+1aQfi2qPQdXlOI3dXORYkhuG2CSMhl8lMvrevdRuGe3lQu+mVNNyWL1+OsrIyHDp0SOpSiMhOpKenQ6/XQ3b5r36PyMniY2GeXXh+VgRmRgdC19GBY8eOob3ddO2EzMxMZGdn3zA10HAvD8pgICKykJkZdyFn5VuYMONe3PX8XzBuWhp83V2xODEYt47sEEcsX2mg7ujoGNBxrwTOcC0PyltJREQWUHyuCVuKquE+YToS/CYgMCgQ08cF4N6E0dBpW1F7bPDHHs7lQQEGQ6/YxkBEV/TV8KvVdWNLcTWOVDaL+zwUBjyaMhrJN4UBAHRak28ZVEOytZcHvRaDoRccx0BEwNWGXwA3NPyerGvFpp+q0NLeJT5/cqg3fNwqEB3k2efxbGG8Qm/YxkBE1If33nvPZPvvf/879AYBXx6txdv5FWIoeLopkH3rGCyaGgy3Pj5Zh7sheTAYDEREZujSG/DWd7/i2+P1EGAckBsT7IP/vCsGUyNH9Pv9w92QPBgMhl6oVCrExsYiKSlJ6lKISEKPP/44AEB2edyB56Q7xUnvFDIZ7k0IxVNpYzHC07XXY7i7uyM2Nhbu7u5iQ3JmZiY2bdpkchvp2udJiW0MvWAbAxEBxh5Bn2/dig2f7QRC4hAWfzsAIMDLDdmpkSajl3vj4eGBuLg4cbu3huTrnycVBgMRUR9a2rvwq/vNCJoVCC9PYwjEhXgh0b8DwV6OedPFMc+KiMgCzjRo8ZevT2PP7i/x48ev49dD32FhYhgemByAM6eOD3iAmr3hFQMR0XUEQUDB6QvYfqQGpw/txddv/hGADKcLtuHJO6Kxu61t2Ca0kwKDgYjoGp3dBvzz0HkUn2sCABzbu/XyI8YeSK+++ioKCwttehzCUPFWUi/YK4nI+Vxo68Qbe8rFUACAwOuW12xsbLT5cQhDxWDoBWdXJXIuv9a3Yc23p8XV1ZQuciy9LQr/k7vc5Hm/+93v+hyHYCtdToeCt5KIyOkdPHMRW4rOo9tgvF00ykeJx2ZEIcTPHfERphPYLV68GHFxcb1OaGcrXU6HgsFARE7LIAjYUVKLvScbxH03h/jg0VvHwMNNIe67ftzBcE5oJwUGAxE5pfYuPT7cX4myWo247/abAnFfQigUclkf3zl0tn67icFARE6nsU2HdwvOok5jHIegkMnwm8tLbg7EUD/Ybf12E4OBiJzKuQuX8LeCCrTpugEYZ0XNuS0KNwV7D/gYtv7BPlQMhl5woR4ix1Na3YIP9leiS28AAAT7umPZzCgEXdcl1dmxu2ov2F2VyLH8eLoRf//hnBgKE4K88XzG+AGFgq23CVgarxiIyKEJgoAvf67DnuP14r6pkSOwOCUCroqB/W3s6LeOrsdgICKH1aU34B+F501GMmfcMgrzJoeI6yvYi+G8amEwEJFDau/U470fz6K8vg0AIIMMv0kMxczoQIkrG5zhvGphMBCRw2m61In131egrsXYHdVVIUf2rZGYFMZFtwaCwUBEDqG9vR1nzpyBe0AoPiisRUt7FwDAW+mCZTOjBrTSGhkxGIjIIXR0dGDv4ZP41UUHg8w4nUWgtxJPzhrL7qhmYndVInII6z76HP/33kc4UbgXABA10hPPZ0xgKAwCrxiIyK4JgoCX132El597HDK5HKcLtmH8a3/D8gceg5sL//YdDAYDEdktvUHAlqLz2L7rG8jkCggGPeRyBVzryxgKQ8B/uV5wBTci26br1uPdHypw4MxFhMcmQTDoIZPLYTDokZ6eLnV5dk0mCIIgdRG2TKPRwM/PDy0tLfD19ZW6HCICoGnvwjsFFTjfZFxtzUUuw0j1Iez/equ4mA6ZMuezjLeSiMiuqDUdWP99BS5qOwEAHq4KPD4zCiNdIxAd4uuwi+cMJwYDEdmNMw1avPtDBS51Gmc9HuHpiidnjcNoP3c0NTX18900UAwGIrILRyqb8fHBKnQbjLOjhvt74Pe3j4Wfh6vElTkeBgMR2bz8kw3YdqQWAoxNojeH+CDntjFwd726LrOzTY1tTQwGIrJZBkHAtiM1+P5Uo7gvZWwAfpsUBhe5aadKZ5sa25oYDERkk7r0Bnx0oBJHz7eI++6KC8ZdE4Ptbspse8NgICKbo9V1490fzqKiUQsAkMtkWJQUjunjAiSuzDkwGIjIplxo68T678+gvlUHAFC6yPHobWMQO5rjiIYLg4GIbMb5pna8830FNB3GKbN93V3xxKyxCPf3kLgy5+JwU2IsWLAA/v7+eOCBBwa0n4hswyl1G97K+1UMhWBfd/zbnAkMBQk4XDA899xz+PDDDwe8n4ikd6SyGe98fwYd3caBa1GBXnhu9ngEeLlJXJlzcrhbSWlpacjPzx/wfiKSxpYtW/DFF19gQlI6LgYlimMU4kJ98eitYzg7qoRs6l++oKAA8+fPR2hoKGQyGbZt23bDc1QqFaKiouDu7o6UlBQUFhYOf6FENCQ7duzAb3/7W3y8aRNWPbsUvxYZF9eZPi4Aj82IYihIzKb+9bVaLeLj46FSqXp8fPPmzcjNzcWqVatw+PBhxMfHY+7cuaivrx/mSoloKPK++w4ymRyCwQDIZDh/vAh3xgbjd0nhUMg5RkFqNhUMmZmZeOWVV7BgwYIeH1+zZg2WLVuGnJwcxMbGYv369fD09MSGDRuGuVIiGqwuvQGnLnRBEIxzHkEQMCkyCPdMDuHANRthU8HQl87OThQXF5tMqSuXy5GRkYEDBw5Y7HV0Oh00Go3JFxFZhq5bj/XfV6DuQgtkMuPHj0wuRyA7HtkUuwmGxsZG6PV6BAcHm+wPDg5GXV2duJ2RkYGFCxdi586dCA8PF0Ojt/3XW716Nfz8/MSviIgI650UkRNp79Tjr/kVKK9vg4ubBwTBAJlMBsFgQFpamtTl0TUcrlfSnj17zNp/vRUrViA3N1fc1mg0DAeiIWrTdePt/DM439SOM0X5KNr+njEUBAG5ubnIysqSukS6ht0EQ2BgIBQKBdRqtcl+tVqNkJAQi72OUqmEUqmESqWCSqWCXq+32LGJnJGmvQuq/DOoa+kAAKhPFkOhUECv10Mul6O9vV3iCul6dnMryc3NDYmJicjLyxP3GQwG5OXlITU11eKvt3z5cpSVleHQoUMWPzaRs2hp78Jb3/0qhoKfhyuWL84SQ8FgMGDGjBkSV0nXs6krhra2NpSXl4vbFRUVKCkpQUBAACIjI5Gbm4vs7GxMmzYNycnJWLt2LbRaLXJyciSsmoh60trRjb/uvToZnr+nG5anj0OQTyw8XT7F9u3bMW3atF57IZJ0ZIIgCFIXcUV+fj7S09Nv2J+dnY2NGzcCANatW4c///nPqKurQ0JCAv7v//4PKSkpFq/l2ltJp06dQktLC3x9Obsj0UBc6uyGaq+xTQEAArzc8Owd4+HPKS4ko9Fo4OfnN6DPMpsKBltkzj8mEQHtXXq8nX8G5y5cAgCM8HDFM7PHI9BbKXFlzs2czzK7aWMgItvXpTfgvYKzYih4K12wPJ2hYG8YDERkEYIgYNPBKpQ3tAEAvNxcsDx9HEb5MhTsDYOhFyqVCrGxsUhKSpK6FCK78MXPdThS1QwAcFPI8eSssQgdwSHN9ohtDP1gGwNR//aVX8CnRecBADLIsOz2KMSF8v1iS9jGQETDpqxGgy1F1eL2A4mhDAU7x2AgokFTazqwcf85cZGdO24OwozoQImroqFiMPSCbQxEfevo0mPDj+eg6zZOnz0lYgSy4kdLXBVZAtsY+sE2BqIbCYKAjfvPoaSqBQAQOsIDz2eMh9JFIXFl1Bu2MRCRVe092SiGgoerAktvG8NQcCAMBiIyy2l1G744WituPzw9EkE+HKvgSBgMvWAbA9GN2nTd+PBAJQyX70DfGRuMiWG8xepo2MbQD7YxEBkJgoD3953D0fPGW0gxwT54Mm0s5Fyn2S6wjYGILK7wbJMYCl5uLnhoegRDwUExGIioX41tOvyr+OogtkXJ4fDzcJWwIrImBgMR9UlvEPDxwSpxvML0cQGID/eTuCqyJgYDEfXpuxMNqGjUAgACvd2wYEqoxBWRtTEYesFeSURAvUaH3b/UATBOjvfw9Ei4u3K8gqNjr6R+sFcSOStBELDuuzPi+gp33ByEexN4tWCv2CuJiIbs4JmLYigEershc2KIxBXRcGEwENENNO1d2F5ydXTzb6eFw82FHxfOgj9pIrrB54dr0N6lBwAkRfkjJsRH4opoODEYiMjEsRqNuESnt9IF97EXktNhMBCRqEtvMBnIdt+UUHgrXSSsiKTAYOgFu6uSM9p7ogEXtJ0AgOhR3pg2ZoS0BZEk2F21H+yuSs6i6VInXvvqJDr1BshlMvxxbjRCR3hIXRZZCLurEpHZtpfUolNvnPZiRvRIhoITYzAQEcrr23CkshmAscE5c2KwtAWRpBgMRE5ObxDwr8M14vY9k0Pg6cYGZ2fGYCBycgfPXERNczsAINzfA9PHBUhcEUmNwUDkxHTdeuy6PEkeAPxmahgX3yEGA5Ez23uiAa0d3QCAhAg/jAvykrgisgUMBiInpWnvwncnGgAACpkM8yaPlrgishUMBiIn9fUxtbgq260TRiLIRylxRWQrGAxETqheo8P+Xy8CAJQucsyNY/dUuorB0AtOiUGO7Mufa2G4POnB7JtHwced3VPpKk6J0Q9OiUGO5kyDFm/mlQMAfN1d8T/zYqB04XKdjo5TYhBRjwRBwI6jVxfguWtiMEOBbmDW9eOOHTvMfoE5c+bAw4NzrhDZgtJqDSoatQCAYF93DmajHpkVDPfdd59ZB5fJZDh9+jTGjRtn1vcRkeXpDQK+/PnqYLb58SFQyDmYjW5k9q2kuro6GAyGAX15enpao2YiGoSDZy5CrekAAIwN9MLEULaZUc/MCobs7Gyzbgs9/PDDbLAlsgG6bj12/6IWt+9NGA0Zp76gXph1K+n999836+Bvv/22Wc8nIuvIP9EITUcXACA+3A9jAzn1BfWOvZKIHFxrRzfyTtQDAOSc+oIGgMFA5OBMpr4YH4BRvpz6gvpm0WDQarU4cuQIGhsbLXlYIhqkhlYd9pdfAGCc+uIursxGA2CxYHjttddw//334/PPP8cf/vAHLF26FO3t7ZY6/IAsWLAA/v7+eOCBB0z2f/nll4iJiUF0dDTee++9Ya2JSEpf/lwL/eXJDe64OQg+7q4SV0T2wCITpLz//vu4ePEivv76a3Hfli1b8Pzzz+Odd96xxEsMyHPPPYelS5figw8+EPd1d3cjNzcXe/fuhZ+fHxITE7FgwQKMHDly2OoiksLZRi1KqloAAD7uLki/OUjiisheWOSKYfPmzXjxxRcBAEuXLoVarcbChQtx6NAhSxx+wNLS0uDj42Oyr7CwEHFxcQgLC4O3tzcyMzPxzTffDGtdRMNNEARsL7k69UXmxBBOfUEDZpFgcHV1RWdnJwDgtttuEwe2GQyGAR+joKAA8+fPR2hoKGQyGbZt23bDc1QqFaKiouDu7o6UlBQUFhb2e9yamhqEhYWJ22FhYaiurh5wXUT2qLRagzOXp74Y5aNEyjh/iSsie2KRYHj88cfxH//xH9Dr9Xjsscfg4+ODN954A5mZmQM+hlarRXx8PFQqVY+Pb968Gbm5uVi1ahUOHz6M+Ph4zJ07F/X19ZY4BSKH0W0wmEyUl5UwGi5ydkCkgbNIG8O9994LrVaL2bNnY8yYMaitrcWtt96KV155ZcDHyMzM7DNI1qxZg2XLliEnJwcAsH79enz11VfYsGGDeBurJ6GhoSZXCNXV1UhOTu71+TqdDjqdTtzWaDQDPgciW3Dg14toaDX+Do8P4tQXZD6Lrc6xePFiLF68GBcvXoS/v79Fh9t3dnaiuLgYK1asEPfJ5XJkZGTgwIEDfX5vcnIyfvnlF1RXV8PPzw+7du3CSy+91OvzV69ejZdfftlitRMNp/bO66e+COXUF2Q2i19fBgQEWPwXsbGxEXq9HsHBpn2wg4ODUVd3dbbIjIwMLFy4EDt37kR4eDgOHDgAFxcXvP7660hPT0dCQgJeeOGFPnskrVixAi0tLeJXVVWVRc+FyJryjtejTdcNAJgaOQJjRnIiSzKfQ63nt2fPnh73Z2VlISsra0DHUCqVUCo5MpTsT9OlTuSfMg4udZFz6gsaPLtokQoMDIRCoYBarTbZr1arERISYpXX5JrPZG+++rkOXXpjT8CZNwVipLebxBWRvTL7imHp0qUDet6GDRvMLqY3bm5uSExMRF5enrhYkMFgQF5eHp5++mmLvc61li9fjuXLl4vrpBLZsvNN7Sg62wwA8HRT4M7YUdIWRHbN7GDYuHEjxowZgylTpkC4PNTeEtra2lBeXi5uV1RUoKSkBAEBAYiMjERubi6ys7Mxbdo0JCcnY+3atdBqtWIvJSJnJQgCth2pgQDj+3FuXDA83RzqLjENM7N/e5566in84x//QEVFBXJycvDwww8jIGDo68YWFRUhPT1d3M7NzQVgXBxo48aNWLRoERoaGrBy5UrU1dUhISEBu3fvvqFB2lJUKhVUKhX0er1Vjk9kKT9Xa3C6vg0AMNLLDbdN4HQvNDQyYRB/9ut0Onz++efYsGED9u/fj3vuuQePPfYY7rzzTofrGnflVlJLSwtXoyOb06U34LWdJ3FRa5x5YOmMKMSH89Yn3cicz7JBNT4rlUo8+OCD+Pbbb1FWVoa4uDj84Q9/QFRUFNra2gZVNBGZ77vjDWIoxAT7YHIY/3ihoRtyryS5XA6ZTAZBEHjbhWgYNWk78e3xqyuz3T+Vg9nIMgYVDDqdDv/4xz8wZ84c3HTTTSgtLcW6detQWVkJb29vS9coCXZXJVu3vaRW7J56+02BCPFzl7gichRmtzH84Q9/wD//+U9ERERg6dKleOihhxAYGGit+iTHNgayRafVbVi391cAgLfSBf9zz83wcOO02tQ7cz7LzA4GuVyOyMhITJkypc/L1s8//9ycw9osBgPZmm6DAX/efRp1mg4AwIPJEZg+bug9A8mxmfNZZnZ31SVLljjFfUx2VyVblXe8QQyFyABPpIzlWgtkWYPqrupMeMVAtqReo8P/230S3QYBcpkML9wZjXB/D6nLIjtg9e6qRDT8BEHA5kPn0W0w/i2XfnMQQ4Gswqxg+Pnnn81arvPYsWPo7u42uygiutFPFU0obzCOEwr0dsNdcdYZ9U9kVjBMmTIFFy5cGPDzU1NTUVlZaXZRRGSqtaML247UiNu/nRYONxde8JN1mNX4LAgCXnrpJXh6Dmzxj87OzkEVZQvY+Ey2QhAEfFpUjfYu4+9iUpQ/YkJ8JK6KHJlZjc9paWlm90j65JNPMHq0/S4YwsZnktpPFRfxyU/GlQS9lS54MTMGPu6cPZXMY7Xuqvn5+UOpi4jMdKGtE58XX3MLKSmcoUBWx5uURDbKIAjY9FMlOrqNt5BSxgZw5lQaFgwGIhu190QDfm3QAgACvNywYGqoxBWRs2AwENmg803t2FlaBwCQQYZHpkfCw5VzIdHwYDD0grOrklTaO/V4f985cSDb7FuCMC7IS+KqyJkMaUqMrq4u1NXV4dKlSwgKCrLIEp+2hr2SaDgJgoCN+8+hpKoFgHEupGdnj4ergn/D0dBYdUqM1tZWvP3225g1axZ8fX0RFRWFW265BUFBQRgzZgyWLVuGQ4cODbp4ImdWcPqCGAoergo8eusYhgINO7N+49asWYOoqCi8//77yMjIwLZt21BSUoJTp07hwIEDWLVqFbq7u3HnnXfirrvuwunTp61VN5HDOa1uw/ZrRjc/ND0CI73dJKyInJVZHaIPHTqEgoICxMXF9fh4cnIyli5divXr1+P999/HDz/8gOjoaIsUSuTIGlp1eH/fOeiFqxPkTQpj11SSxqDaGPR6Pb744gvMnj0bPj6OPTSfbQxkbe1deqz9tlxcY+GW0T5YNnMsFHLHX/eEho/Vp91WKBR48MEH0dDQMKgCicjIIAj4cH+lGArBvu7ITh3DUCBJDbpVKykpCRUVFZasxaawuypZmyAI2F5Si7JaDQDAy80Fy2ZGce1mktygu6tu2bIFf/nLX/DZZ58hIiLC0nXZDN5KImvZU1aPL36uBQAoZDI8mTYONwV7S1wVOSqrrvl8xaJFiwAAcXFxyMrKQlpaGqZMmYJJkybBzY09KYj6sq/8ghgKALBwWhhDgWzGoIOhoqICR48eRUlJCY4ePYrVq1fj7NmzcHFxQUxMDH7++WdL1knkMI5UNmNLUbW4PX/yaKSOHylhRUSmBh0MY8aMwZgxY5CVlSXua21tRUlJCUOBqBdltRp8fLASAq52S519S5DEVRGZGtKUGM6AbQxkKWU1Gvz9x3Povrxu+vRxAfhdUrjZi18RDcawtDEQ0cCVVrdg4zUT48WH++G30xgKZJsYDERWVlLVjA/3V4qjmqdEjMDDqREcq0A2a1DjGLRaLUpLS3t87NixY2hraxtSUUSO4tDZJnxwTShMG+OPR1Ij4SLnxHhkuwb129nV1YWUlBQUFhaa7C8rK8OUKVMYDOT0BEFA3vF6fHywEobLoZA8NgAPTeeVAtm+QQXDiBEjMG/ePHz44Ycm+z/66CPMnj0bISEhFimOyB4ZBAGfH6nBjqNXxynMmDASi5PDIWebAtmBQV/PZmdnY/Pmzeju7gZg/Atp06ZNyMnJsVhxUuKUGDQYXXoDPtxfiYJTjeK+eyaF4IHEMDY0k90YdHdVvV6P8PBwrF+/Hvfeey/27t2L3/zmN6irq3Ookc/srkoD1XSpExt+PIfKi5cAAHKZDL9LCkfKOMdb2ZDsj9VnVwWMM6w+9NBD4u2kjz76CIsWLXKoUCAaqNPqNvzl69NiKLgp5Fg2M4qhQHZpSN1Vs7OzkZycjOrqavzrX//C119/bam6iOyCIAj4/lQjtpfUio3MI73c8NiMKIT5e0hcHdHgDHnkc2JiInx8fFBXV4cTJ05Yqi6bwVtJ1JvObgP+eeg8is81iftuDvHBktRIeCk5RIhsy7DcSrpiyZIlKCgowJIlS4Z6KCK7caGtE2/mlZuEwpxbRuGJWWMZCmT3hvwb/Mgjj6C5uRlLly61RD1ENu9kXSs+2F8JbaexR57SRY6HUiIRH8E1mskxcBK9fvBWEl0hCAK+O9GAL47WibOjBvko8diMKIz2c5e4OqK+cRI9IgvTdevxj5/O40hVs7gvLtQXj0yP5FKc5HAYDET9aGjV4b0fz6KupUPcd1dcMO6aGMxBa+SQGAxEfThWo8FHByrR3qUHALi7KvDw9AhMCmN7Ajkup5ji8S9/+Qvi4uIwceJEfPzxx1KXQ3ZAEAR8/Ysa7xacFUMhxNcdL8yJZiiQwzPriqGyshKRkZEDfn51dTXCwsLMLsqSSktL8cknn6C4uBiCICA9PR3z5s3DiBEjJK2LbFd7lx6bDlaitFoj7osP98PilAi4u7I9gRyfWVcMSUlJeOKJJ3Do0KFen9PS0oJ3330XEydOxL/+9a8hFzhUx48fR2pqKtzd3eHh4YH4+Hjs3r1b6rLIRqk1HVjzzWkxFGSQYd7k0ci5bQxDgZyGWcFQVlYGLy8vzJkzByEhIbjnnnuwbNkyPPPMM3j44YcxdepUjBo1Chs2bMCf/vQnPPvss0MusKCgAPPnz0doaChkMhm2bdt2w3NUKhWioqLg7u5+wzoREydORH5+Ppqbm9HU1IT8/HxUV1cPuS5yPD+fb8Gab8pR36oDAHi6KfD7WVGYEzuKjczkVMwKhpEjR2LNmjWora3FunXrEB0djcbGRpw+fRoA8NBDD6G4uBgHDhzA3XffbZECtVot4uPjoVKpenx88+bNyM3NxapVq3D48GHEx8dj7ty5qK+vBwDExsbi2WefxR133IH7778f06dPh0LBv/zoKkEQsLO0Dn//8Sw6uo3tCaEjPPDCndGIHc2xK+R8zB7g9sorryAxMRFTp05FcHCwterqkUwmw9atW3HfffeJ+1JSUpCUlIR169YBAAwGAyIiIvDMM8/gxRdfvOEYjz/+OBYsWIB77rmnx9fQ6XTQ6XTitkajQUREBAe4OShdtx4fH6zCz+dbxH1TI0fgd8nhULrwDwhyHFYd4LZy5UrxsjokJARTp05FYmKiGBbD2djc2dmJ4uJirFixQtwnl8uRkZGBAwcOiPvq6+sxatQonDx5EoWFhVi/fn2vx1y9ejVefvllq9ZNw6+9vR1nzpzBuHHj4OFhnPX0Qlsn3vvxLGqa2wEY2xOyEkYjPSaQt47IqZkdDElJSaitrUVOTg4CAwNx+PBhfP7553jttdeg1+sRFBSEqVOnYufOndao10RjYyP0ev0NVy7BwcEmM73ee++9aGlpgZeXF95//324uPR+2itWrEBubq64feWKgexbR0cHysrKEBoaCg8PD5xSt2HjvnPifEcergosuTWSt46IMIhg+Omnn7Bx40b813/9F5KSkrBmzRqMHz8eOp0OJSUlOHz4MI4cOWKNWgft2quH/iiVSiiVSitWQ1LbV34BnxVXi+snjPJRYtnMsRjly587ETDIAW6PPvooTp06hZiYGEydOhUrVqyAXq9HSkoKnnrqKfztb3+zdJ09CgwMhEKhgFqtNtmvVqsREhIypGNzzWfHIwjA18cv4NOi82Io3DLaB7lzohkKRNcY9Mhnb29v/OlPf0JRURF++eUXTJgwQVzmc7i4ubkhMTEReXl54j6DwYC8vDykpqYO6djLly9HWVlZn2M2yH58+dVO/L93PsYHW7aL+9JjgvD728dyEjyi6wxprqTu7m7odDo8+OCDUKvVyMnJwbx58xAQYLl1btva2lBeXi5uV1RUoKSkBAEBAYiMjERubi6ys7Mxbdo0JCcnY+3atdBqtcjJybFYDWTftvxrK5Y88jBkcjmE73bA9YU38ccnFmNmdKDUpRHZJLOD4X//939RWlqK0tJSnDhxAu7u7pg8eTKSk5PxxBNPwM/PsvPIFBUVIT09Xdy+0jCcnZ2NjRs3YtGiRWhoaMDKlStRV1eHhIQE7N69e8hdaVUqFVQqFfR6/ZCOQ9K6qO3Emx9th0yugGDQQyZXYITmJEOBqA9mj2OQy+WIiopCdnY2HnzwQdx0003Wqs0mcKEe+1XX0oG/5p/BkR++xZevP2+8YjAYsH37dmRlZUldHtGwMuezzOxgmDVrFkpKStDa2govLy9MnjwZU6dOFb8mTpzoUCOLGQz2qfLCJaz/vkLsjlp3NB/NR79F9sOLsXjxYomrIxp+Vg2GK06fPo3i4mIcPnxY/GpuboZSqcSkSZNM5iuyR9feSjp16hSDwY6U17fh3YKr01tEBHjid/H+OFCwFxkZGfD395e4QqLhNyzB0JOKigoUFRXhyJEjeO211yx1WEnxisG+lNVqsOHHc+jSGwAAE4K88fjtUeho02DPnj0MBnJakq35PHbsWIwdOxYLFy605GGJBqSsVoO//3AO3QZjKMSO9sXSGWPgqpCjo5/vJaKruLQnOYTjta0moZAQ4YdHUiPhIneKRQqJLIrvml5w5LP9OF7bivd+OCuGwpSIEViSOsYkFNzd3REbGwt3d3epyiSyGxZtY3BEbGOwbafqWvFOwVleKRD1w5zPMr57yG6dbdTivR+vhkJ8OEOByBL4DiK7VNPcjncKKqDrNobCpDBfLLmVoUBkCXwXkd1pbNNh/fcVuNRpHKcQPcob2beOYSgQWQjfSb1g47Nt0rR34e38M2hp7wIAjBnpicdnRsFVwV9lIkth43M/2PhsO3Tderz13RlUXbwEAAjxdcezs8fDS8le10T9YeMzORyDIODD/ZViKPh7uuGptHEMBSIrYDCQzRMEAZ8frsEvNRoAxvWZn5g1FiM8XSWujMgxMRjI5uWfbMQPpxsBAC5yGZbOiMJoPw5UI7IWBgPZtNLqFmwvqRW3FyWF46ZgbwkrInJ8DIZesFeS9GpbOvDRgUoIMPaPuCsuGMljLbdsLBH1jL2S+sFeSdLQ6rqx5ttyNLbpAABTIkcgOzUSMplM4sqI7BN7JZFd0xsEfLC/UgyFcH8PLE6OYCgQDRMGA9mcHUdrcVLdCgDwVrrg8RlRcHPhryrRcOG7jWxK8bkm5J9sAGDsgfTYjCj4e7lJXBWRc2EwkM2oa+nA5kPnxe37p4ZhXJCXhBUROScGA9kEXbce7+87J86Wmjw2ALeOZw8kIikwGHrB7qrDRxAE/LPwPOo0xpWZQ0d4YGFiGBubiSTC7qr9YHdV6/vxdCO2FFcDANxdFPj3udEI8lFKXBWRY2F3VbIb55vasfVIjbi9OCWCoUAkMQYDSaaz24APD1Si22C8aE2LCUJ8hJ/EVRERg4Eks62kBurL7Qrh/h6YNzlE4oqICGAwkER+Pt+CfeUXAACuCjmWpEZyFTYiG8F3Ig275ktd+GfhteMVQhHsy2m0iWwFg4GGlSAI+PhgJbSd3QCA+HA/pI7jeAUiW8JgoGH13YkGnK5vAwCM8HDFoqRwjlcgsjEMBho255va8dXPdQAAGWR4aHok12wmskEMhl5w5LNldRsM2PRTFfSXx1POviWIK7ER2SiOfO4HRz5bxs7SOnx9TA3AOOXFC3dOgIucf5cQDReOfCabUnXxEr4tqwcAKGQyPJQSwVAgsmF8d5JVdemNt5AMly9M58SOQri/h8RVEVFfGAxkVV8fU6O2xTi6OWyEB+bEjZK4IiLqD4OBrKbywiXkHTeuxsZbSET2g+9SsorrbyHNjQtGGG8hEdkFBgNZxbdl9eLCO+H+HpgdGyRxRUQ0UAwGsri6lg7kHWcvJCJ7xXcrWZQgCPjnofPiGgt33BKE0BG8hURkTxgMZFH7f72IikYtACDQW4m5ccESV0RE5mIwkMW0tHfhi6O14vaipHCusUBkh5ziXfvGG28gLi4OsbGxePbZZ8FZQKxj6+EatHfpAQDJYwM4FxKRnXL4YGhoaMC6detQXFyM0tJSFBcX4+DBg1KX5XCO1WhwpKoZAOCtdMG9CaOlLYiIBs0p5jzu7u5GR4ex62RXVxdGjeLoW0vSdevxadHVFdnumxIKb06nTWS3bP6KoaCgAPPnz0doaChkMhm2bdt2w3NUKhWioqLg7u6OlJQUFBYWio8FBQXh3//93xEZGYnQ0FBkZGRg/Pjxw3gGjm/XL2o0X+oCAMQE+2DamBHSFkREQ2LzwaDVahEfHw+VStXj45s3b0Zubi5WrVqFw4cPIz4+HnPnzkV9vbEffVNTE7788kucPXsW1dXV2L9/PwoKCobzFBxabUsHCk42AgBcFXIsnBbGFdmI7JzNB0NmZiZeeeUVLFiwoMfH16xZg2XLliEnJwexsbFYv349PD09sWHDBgDAnj17MGHCBAQEBMDDwwP33HNPn20MOp0OGo3G5It6JggCPiuuFhffybglCEE+SomrIqKhsvlg6EtnZyeKi4uRkZEh7pPL5cjIyMCBAwcAABEREdi/fz86Ojqg1+uRn5+PmJiYXo+5evVq+Pn5iV8RERFWPw97dbiyGeWX128O9HbD7FvYdkPkCOw6GBobG6HX6xEcbDqIKjg4GHV1xrWFp0+fjrvvvhtTpkzB5MmTMX78eGRlZfV6zBUrVqClpUX8qqqqsuo52KuOLj22l1wds3D/1DCOWSByEE7RdeTVV1/Fq6++OqDnKpVKKJW8HdKf3cfUaGk3NjhPDPVFXCiXPSVyFHb9J15gYCAUCgXUarXJfrVajZCQkCEdW6VSITY2FklJSUM6jiO6vsH5/qlhEldERJZk18Hg5uaGxMRE5OXlifsMBgPy8vKQmpo6pGMvX74cZWVlOHTo0FDLdCg9NTiP9HaTuCoisiSbv5XU1taG8vJycbuiogIlJSUICAhAZGQkcnNzkZ2djWnTpiE5ORlr166FVqtFTk6OhFU7LjY4Ezk+mw+GoqIipKeni9u5ubkAgOzsbGzcuBGLFi1CQ0MDVq5cibq6OiQkJGD37t03NEibS6VSQaVSQa/XD+k4jkTXrceOo2xwJnJ0MoEzyvVJo9HAz88PLS0t8PV17gbWXaV12H3M2J4TO9oXT8waK3FFRDRQ5nyW8c89GpCmS53IO9EAwLgq24IpoRJXRETWwmDoBXslmdpRUosuvQEAMPOmQIzyZZdeIkfFW0n94K0k4EyDFm/mGTsAeLm54H/mxcDTzeabp4joGryVRBYjCAK2HqkRt++eFMxQIHJwDAbq06GzTai8eAkAMNrPHanjR0pcERFZG4OhF2xjMHZP/fLnOnH7/imhUMg5pTaRo2MbQz+cuY3hy59r8W2ZcV2LSWG+eHwmu6cS2Su2MdCQNV3qRP7l+ZBc5DLcm8DuqUTOgsFAPdpZqr7aPTU6kAvwEDkRBgPdoLqpHYcqmgAAHq4KzInlfEhEzoTB0AtnbnzeXlILAcampzvjguGlZPdUImfCYOiFs067fby2FSfVrQCAAC83zIxm91QiZ8NgIJFBEExmT503OYSzpxI5Ib7rSXSoogk1ze0AgMgAT0yNHCFtQUQkCQYDAQA6uw34qvTqYLZ7E0ZDJuNgNiJnxGAgAED+yQa0tHcBACaG+mLCKG+JKyIiqTAYeuEMvZLa29tx7NgxNDS3Iu+4ca0FuUyG+fGjJa6MiKTEKTH64chTYjQ1NWHPnj24FByPwzXGtoVbx4/EoqRwiSsjIkvjlBg0YK1dchyqNHZPVbrIkTlxaGtlE5H9YzA4uV80bjBcvmi84+Yg+Hq4SlwREUmNweDEqpt1OH/JGATeShek3xwkcUVEZAsYDE5szQefoWT7O6g4/D3ujAuG0kUhdUlEZAM4CY6TevvDzXj/5Wcgk8txumAbls0aD9x0n9RlEZEN4BVDLxy5u6ogCPjntt2QyRUQDAbIFQr8WFAgdVlEZCMYDL1w5En0Sqs18B2XAMGgh0wuh0GvR1pamtRlEZGN4K0kJ2MQBHz1cx3GTUvDXc//Be3Hv8fvl/wOWVlZUpdGRDaCweBkDp1tQp2mAwAwa/ZcjEuNxJw5GRJXRUS2hLeSnEiX3oBdpWpx+86bA8B58ojoegwGJ7Kv/AKaLnUCAG4Z7YOxIz0kroiIbBGDwUl0dOnxbVm9uD1vMifKI6KeMRicRP7JBrTpugEAUyNHINyfVwtE1DMGgxPQ6rqx92QjAOO02ndPCpG4IiKyZQwGJ/DdiQZ0dOkBACnjAhDko5S4IiKyZQyGXjjKyOfWjm4UnDJeLbjIZZgbO0p8zN3dHbGxsXB3d5eqPCKyQVyopx/2vlDP1iM1yD9pXJ1tZnQgHkgMk7giIpICF+ohAEDzpS7sK78AAHBVyDHnmqsFIqLeMBgc2LdlanTpDQCAGdEj4cdFeIhoABgMDuqithMHz1wEYFyyc/bNvFogooFhMDior4+p0W0wNh/dflMgfNw5LRYRDQyDwQE1tOpQWNEEAPBwVeAOLtlJRGZgMDigr4+pYbjc2Sw9JgiebrxaIKKBYzA4mLqWDhSdbQYAeLm54PaYQGkLIiK7w2BwMLt/UUOA8Wph9i1B8HBVSFwREdkbBoMDqW5qx5GqZgCAj7sLZkSPlLYgIrJLDAYHsuuXOvH/59wyCkoXXi0QkfmcIhhOnjyJhIQE8cvDwwPbtm2TuiyLqrxwCaXVGgDACA9X3DqBVwtENDhO0V0lJiYGJSUlAIC2tjZERUVhzpw50hZlYbuPXbNkZ1wwXBVOkflEZAVO9+mxY8cOzJ49G15eXlKXYjGVFy7hWI3xasHf0w0p4/wlroiI7JldBENBQQHmz5+P0NBQyGSyHm8DqVQqREVFwd3dHSkpKSgsLOzxWJ9++ikWLVpk5YqH17VXC3NiR8FFbhc/ViKyUXbxCaLVahEfHw+VStXj45s3b0Zubi5WrVqFw4cPIz4+HnPnzkV9fb3J8zQaDfbv34+77757OMoeFtdeLYzwdOXVAhENmV20MWRmZiIzM7PXx9esWYNly5YhJycHALB+/Xp89dVX2LBhA1588UXxedu3b8edd97Z58I0Op0OOp1O3NZoNBY4A+sxaVuIDebVAhENmd1/inR2dqK4uBgZGRniPrlcjoyMDBw4cMDkuQO5jbR69Wr4+fmJXxEREVap2xJ4tUBE1mD3wdDY2Ai9Xo/g4GCT/cHBwairu9qvv6WlBYWFhZg7d26fx1uxYgVaWlrEr6qqKqvUbQlf82qBiKzALm4lWYKfnx/UanW/z1MqlVAqlcNQ0dBUXryEX665Wkgey6sFIrIMu/8TMzAwEAqF4oYPfbVajZCQkEEfV6VSITY2FklJSUMt0Sq+/uWanki3jOK4BSKyGLv/NHFzc0NiYiLy8vLEfQaDAXl5eUhNTR30cZcvX46ysjIcOnTIEmVaVNV1Vwsp4wIkroiIHIld3Epqa2tDeXm5uF1RUYGSkhIEBAQgMjISubm5yM7OxrRp05CcnIy1a9dCq9WKvZQczbVtC7xaICJLs4tgKCoqQnp6uridm5sLAMjOzsbGjRuxaNEiNDQ0YOXKlairq0NCQgJ27959Q4O0OVQqFVQqFfR6/ZDrt6Sqi9fMicSrBSKyApkgXF7qi3qk0Wjg5+eHlpYW+Pr6Sl0O3vuhQgyGhYlhmBHNhXiIqH/mfJbZxRWDlK7kpi0MdKtuakfR6RoAxhlUbwl0sYm6iMj2XfmsGMi1AK8YenHlVpJOp8OZM2ekLoeIyCKqqqoQHh7e53MYDP1obm6Gv78/Kisr4efnN+jjJCUlDaiHU2/PG+j+vrav//+8vDxERESgqqpqSLfJBnJufT2np8f629fbeV75r0ajGfK5DfVnNpDz6Glff/9vT+dmq7+PfT3PUd9rgiCgtbUVoaGhkPczGJa3kvpx5R/Qz89vSD9QhUIxoO/v7XkD3d/Xdm//7+vra/Vz6+s5PT3W377ezvP6/UM5t6H+zHp7zJxz6+uc7eHcbPX3sa/nOfJ7baB/3LKf4zBZvnz5kJ430P19bff2/0M1kGP19ZyeHutvX2/nOdzn1d/zhnpuUv7M+nveQM7NVn8f+3qeI7/XBoq3kvpha72SLInnZp8c9dwc9bwA+zs3XjH0Q6lUYtWqVXYxf5K5eG72yVHPzVHPC7C/c+MVAxERmeAVAxERmWAwEBGRCQYDERGZYDAQEZEJBoMFNTc3Y9q0aUhISMDEiRPx7rvvSl2SxVRVVSEtLQ2xsbGYPHkytmzZInVJFrNgwQL4+/vjgQcekLqUIfvyyy8RExOD6OhovPfee1KXY1GO9HO6li2+t9gryYL0ej10Oh08PT2h1WoxceJEFBUVYeTIkVKXNmS1tbVQq9VISEhAXV0dEhMTcerUKXh5eUld2pDl5+ejtbUVH3zwAT777DOpyxm07u5uxMbGYu/evfDz80NiYiL279/vEL9/gOP8nK5ni+8tXjFYkEKhgKenJwBAp9NBEIQBzWRoD0aPHo2EhAQAQEhICAIDA3Hx4kVpi7KQtLQ0+Pj4SF3GkBUWFiIuLg5hYWHw9vZGZmYmvvnmG6nLshhH+TldzxbfW04VDAUFBZg/fz5CQ0Mhk8mwbdu2G56jUqkQFRUFd3d3pKSkoLCw0KzXaG5uRnx8PMLDw/HHP/4RgYHDs17CcJzbFcXFxdDr9YiIiBhi1f0bzvOS2lDPtaamBmFhYeJ2WFgYqqurh6P0fjnyz9GS5zac762+OFUwaLVaxMfHQ6VS9fj45s2bkZubi1WrVuHw4cOIj4/H3LlzUV9fLz7nSvvB9V81NZfXSRgxAkePHkVFRQU++eQTqNXqHl/LHs8NAC5evIglS5bgb3/7m9XPCRi+87IFljhXW8Vz6//chvu91SfBSQEQtm7darIvOTlZWL58ubit1+uF0NBQYfXq1YN6jaeeekrYsmXLUMocFGudW0dHhzBz5kzhww8/tFSpZrHmz2zv3r3Cb37zG0uUaRGDOdd9+/YJ9913n/j4c889J2zatGlY6jXHUH6OtvZzut5gz03q99b1nOqKoS+dnZ0oLi5GRkaGuE8ulyMjIwMHDhwY0DHUajVaW1sBAC0tLSgoKEBMTIxV6jWHJc5NEAQ8+uijuOOOO/DII49Yq1SzWOK87MVAzjU5ORm//PILqqur0dbWhl27dmHu3LlSlTxgjvxzHMi52eJ7i8FwWWNjI/R6PYKDg032BwcHo66ubkDHOHfuHGbOnIn4+HjMnDkTzzzzDCZNmmSNcs1iiXPbt28fNm/ejG3btiEhIQEJCQkoLS21RrkDZonzAoCMjAwsXLgQO3fuRHh4uE1+GA3kXF1cXPD6668jPT0dCQkJeOGFF+yiR9JAf4728HO63kDOzRbfW1yox4KSk5NRUlIidRlWMWPGDBgMBqnLsIo9e/ZIXYLFZGVlISsrS+oyrMKRfk7XssX3Fq8YLgsMDIRCobihsVitViMkJESiqizDUc/NUc+rJ458rjw328NguMzNzQ2JiYnIy8sT9xkMBuTl5SE1NVXCyobOUc/NUc+rJ458rjw32+NUt5La2tpQXl4ubldUVKCkpAQBAQGIjIxEbm4usrOzMW3aNCQnJ2Pt2rXQarXIycmRsOqBcdRzc9Tz6okjnyvPzc7OTepuUcNp7969AoAbvrKzs8XnvPXWW0JkZKTg5uYmJCcnCwcPHpSuYDM46rk56nn1xJHPledmX+fGuZKIiMgE2xiIiMgEg4GIiEwwGIiIyASDgYiITDAYiIjIBIOBiIhMMBiIiMgEg4GIiEwwGIiIyASDgYiITDAYiIjIBIOByIrS0tIgk8kgk8nERZzS0tLw/PPPW+01H330UfE1t23bZrXXIcfFYCCysmXLlqG2thYTJ04cltd78803UVtbOyyvRY7JqdZjILKm7u5uuLjc+Jby9PQc1tW6/Pz84OfnN2yvR46HVwxEg3D27FnIZDJ8+umnmDlzJpRKJXbs2DGoY3311Vfw8/PDpk2bABhvNT399NN4+umn4efnh8DAQLz00ku4doZ8g8GAP/3pT5gwYQKUSiUiIyPx6quvWuTciHjFQDQIR48eBQD8+c9/xmuvvYaxY8ciKCjI7ON88sknePLJJ/HJJ59g3rx54v4PPvgAjz32GAoLC1FUVITf//73iIyMxLJlywAAK1aswLvvvos33ngDM2bMQG1tLU6cOGGZkyOnx2AgGoSSkhJ4eXlhy5YtiIqKGtQxVCoV/vu//xtffPEFZs2aZfJYREQE3njjDchkMsTExKC0tBRvvPEGli1bhtbWVrz55ptYt24dsrOzAQDjx4/HjBkzhnpaRAAYDESDcvToUWRlZQ06FD777DPU19dj3759SEpKuuHx6dOnQyaTidupqal4/fXXodfrcfz4ceh0OsyePXuw5RP1iW0MRINQUlKCtLS0QX//lClTEBQUhA0bNsDc1XU9PDwG/bpEA8FgIDKTRqPB2bNnMWXKlEEfY/z48di7dy+2b9+OZ5555obHf/rpJ5PtgwcPIjo6GgqFAtHR0fDw8EBeXt6gX5+oL7yVRGSmo0ePQqFQYNKkSUM6zk033YS9e/ciLS0NLi4uWLt2rfhYZWUlcnNz8cQTT+Dw4cN466238PrrrwMA3N3d8Z//+Z/4j//4D7i5ueG2225DQ0MDjh07hscee2xINREBDAYisx09ehQxMTFwd3cf8rFiYmLw3XffIS0tDQqFQvzwX7JkCdrb25GcnAyFQoHnnnsOv//978Xve+mll+Di4oKVK1eipqYGo0ePxpNPPjnkeogAQCaYe4OTiAYsLS0NCQkJJlcD1vienshkMmzduhX33XffkI5DzodtDERW9te//hXe3t4oLS0dltd78skn4e3tPSyvRY6JVwxEVlRdXY329nYAQGRkJNzc3Pr9nqFeMdTX10Oj0QAARo8eDS8vr0Edh5wXg4GIiEzwVhIREZlgMBARkQkGAxERmWAwEBGRCQYDERGZYDAQEZEJBgMREZlgMBARkQkGAxERmWAwEBGRCQYDERGZ+P/4xZO08KLA8gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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(<r)$ [M$_\\odot$]\")\n",
    "\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "001e19fe",
   "metadata": {},
   "source": [
    "This potential is already implemented in `gala` in `gala.potential.special`,\n",
    "and we can import it with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0dcdc3d5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:31.853107Z",
     "iopub.status.busy": "2025-07-31T14:29:31.852904Z",
     "iopub.status.idle": "2025-07-31T14:29:31.855420Z",
     "shell.execute_reply": "2025-07-31T14:29:31.854940Z"
    }
   },
   "outputs": [],
   "source": [
    "from gala.potential import MilkyWayPotential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e6008c2f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-31T14:29:31.856975Z",
     "iopub.status.busy": "2025-07-31T14:29:31.856812Z",
     "iopub.status.idle": "2025-07-31T14:29:31.870488Z",
     "shell.execute_reply": "2025-07-31T14:29:31.869899Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<CompositePotential disk,bulge,nucleus,halo>"
      ]
     },
     "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.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}