{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gravitation\n",
    "\n",
    "The gravitational force is the dominant force that shapes the structure of galaxies. So we will go through some of the basics of gravitation, and see how the concepts we learn can be applied to galaxies. Einstein's general theory of relativity represents our current best understanding of gravity. Einstein's equations:\n",
    "\n",
    "\\begin{equation*}\n",
    "G_{\\mu\\nu} = 8 \\pi G T_{\\mu\\nu}\\,.\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tensor $G_{\\mu\\nu}$ contains the second derivatives of the metric components, and in the weak field, small velocity limit, the deviations of the metric from the Minkowski metric come in the form of $\\Phi/c^2$. The field equations of GR were derived such that they reduce to Newtonian gravity in such a limit, and reduce to the familiar Poisson equation\n",
    "\n",
    "\\begin{equation*}\n",
    "\\nabla^2 \\Phi = 4 \\pi G \\rho\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mass distribution within an object can be specified by the density, $\\rho$. The solution to the above equation will then determine the potential $\\Phi$. The potential then specifies the forces that will act on test particles, and thus determine their orbits.\n",
    "\n",
    "\\begin{equation*}\n",
    "\\vec{F} = -m \\nabla \\Phi\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The orbits of test bodies in the galactic potential are observable. Thus our goal would be to compare observations to the orbits predicted from the solutions of the Poisson equation above. Agreement would imply that the Poisson equation which is so well tested in the solar system works also on galactic scales. Disagreement could be a result of anamolous terms on either side of the Poisson equation. A change to the left hand side would imply a change in our understanding of gravity, while a change to the right hand side would imply a missing link in our understanding of matter which contributes the gravitational force."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time scales important for dynamics \n",
    "\n",
    "Let us look at order of magnitudes of some important time scales for galactic dynamics.\n",
    "\n",
    "#### Hubble time scale: \n",
    "\n",
    "The Hubble parameter is of the order $H_0$~100 km/s/Mpc. This gives a time scale for the Universe of order 1/H0, which is close to $10^{10}$ years."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Orbital time scale:\n",
    "\n",
    "The orbital time scale can be obtained simply as,\n",
    "\n",
    "\\begin{equation}\n",
    "t_{\\rm orb} \\sim 2\\pi R/v_{\\rm circ}\n",
    "\\end{equation}\n",
    "\n",
    "The circular velocity for Milky Way flattens and reaches a value close to 200 ${\\rm kms}^{-1}$. The orbital time scale for the rotation at the position of the Sun (10 kpc) is about $10^8$ years. Thus the Milky way would make many rotations within the age of the Universe, and we can assume that it should be in dynamical equilibrium. If you compute the orbital time scales close to the edge of the dark matter halos, these come out to be comparable to the Hubble time scale, so the assumption of dynamical equilibrium may not be appropriate on those scales."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Direct collision time scale:\n",
    "\n",
    "The mean free path of particles within the galaxy can be computed as $\\lambda \\sim 1/(n\\sigma)$, where $n$ is the number density of stars and $\\sigma$ the cross section for direct collisions of stars. The collision time scales will then be given by\n",
    "\n",
    "\\begin{equation}\n",
    "t_{\\rm direct} \\sim \\lambda/v \\sim R^3/R_{\\odot}^2/N/v \\sim t_{\\rm orb} (R/R_{\\odot})^2/N\n",
    "\\end{equation}\n",
    "\n",
    "For the Milky Way, the average number density of stars in the galaxy is about 1 star per parsec cubed ($10^{11}$ stars within a disk of 10 kpc with 1 kpc thickness), and the cross section is $4 \\pi R_\\odot^2$, assuming all stars in the galaxy to be like the Sun. This mean free path then turns out to be $\\lambda \\sim 1.5\\times10^{11}$ kpc, which is orders of magnitude larger than the path length within a single orbit.  With velocities of the order $100 kms^{-1}$, the collision time scale is quite large even compared to the Hubble time scale and stellar collisions within the galaxy should be exceedingly rare."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Relaxation time scale:\n",
    "\n",
    "Direct collisions are rare, but particles can change their velocities due to long range interactions. The time scale over which a particle changes its velocity by factor of order unity is called the relaxation time scale. The effect of an encounter of a particle with another particle of mass m with an impact parameter b and moving with velocity v is to change its velocity by\n",
    "\n",
    "\\begin{equation}\n",
    "\\Delta v_\\perp \\sim 2Gm/(bv)\n",
    "\\end{equation}\n",
    "\n",
    "The particle will encounter all sort of impact parameters during its motion through the galaxy. The number of impacts between $b\\pm db$ will be $2 \\pi b db N/(\\pi R^2)$. Thus the cumulative kicks it will gather within an orbital time will be\n",
    "\n",
    "\\begin{eqnarray}\n",
    "|\\Delta v_\\perp^2| &=& \\int_{R_{\\rm min}}^{R} \\left(\\frac{2Gm}{bv}\\right)^2 \\frac{2 \\pi b db N}{\\pi R^2} \\,\n",
    "&=& \\left(\\frac{2Gm}{v}\\right)^2 2 \\ln \\left(\\frac{R}{R_{\\rm min}}\\right) \\frac{N}{R^2} \\,\n",
    "\\end{eqnarray}\n",
    "\n",
    "The value of $R_{\\rm min}$ \n",
    "\n",
    "\\begin{equation}\n",
    "R_{\\rm min} \\sim \\frac{2Gm}{v^2} \\sim \\frac{2GM}{v^2N} \\sim \\frac{R}{N}\n",
    "\\end{equation}\n",
    "\n",
    "would determine the scale where the velocity kicks become comparable directly to the typical velocity. Thus, within one orbital time through the galaxy, the velocity kicks will accumulate energy proportional to\n",
    "\n",
    "\\begin{equation}\n",
    "|\\Delta v_\\perp^2| = \\left(\\frac{2GM}{v}\\right)^2 2 \\frac{\\ln N}{N} \\frac{1}{R^2} = \\left(\\frac{GM}{v}\\right)^2 8 \\frac{\\ln N}{N} \\frac{1}{R^2}\n",
    "\\end{equation}\n",
    "\n",
    "Equating these to the energy of the particle per mass, we get a value for the relaxation time\n",
    "\n",
    "\\begin{equation}\n",
    "t_{\\rm relax} \\sim \\frac{R}{v} \\frac{v^2}{|\\Delta v_\\perp^2|} \\sim \\frac{R}{v} \\frac{v^2}{\\left(\\frac{GM}{v}\\right)^2 8 \\frac{\\ln N}{N} \\frac{1}{R^2} } \\sim t_{\\rm orb} \\left[\\frac{N}{8\\ln N}\\right]\n",
    "\\end{equation}\n",
    "\n",
    "For a system with size of the Milky way, $N \\sim 10^{11}$, which corresponds to about $10^{17}$ years!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Close encounter time scale\n",
    "\n",
    "The close encounter time scale can be obtained by putting $R_{\\rm min}$ instead of $R_{\\odot}$ in the direct collision time scale.\n",
    "\n",
    "\\begin{equation}\n",
    "t_{\\rm close} \\sim t_{\\rm orb} (R/R_{\\rm min})^2/N \\sim N t_{\\rm orb}\\,.\n",
    "\\end{equation}\n",
    "\n",
    "For a system of size of the Milky way, this time scale is close to $10^{19}$ years, even larger than the relaxation time scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Velocity: 207.3499457439042 km/s\n",
      "Orbital timescale: 0.09431398657877353 Gyr\n",
      "Hubble timescale: 13.968571428571428 Gyr\n",
      "Relaxation timescale: 46545504.47438488 Gyr\n",
      "Close interaction timescale: 9431398657.877354 Gyr\n",
      "Direct collision timescale: 185508302266.48767 Gyr\n",
      "==========================\n",
      "Velocity: 2.0734994574390413 km/s\n",
      "Orbital timescale: 0.009431398657877356 Gyr\n",
      "Hubble timescale: 13.968571428571428 Gyr\n",
      "Relaxation timescale: 1.2800013730455844 Gyr\n",
      "Close interaction timescale: 94.31398657877355 Gyr\n",
      "Direct collision timescale: 185508302266.4878 Gyr\n",
      "==========================\n",
      "Velocity: 655.6981012630737 km/s\n",
      "Orbital timescale: 2.9824701279947594 Gyr\n",
      "Hubble timescale: 13.968571428571428 Gyr\n",
      "Relaxation timescale: 53.96959662622539 Gyr\n",
      "Close interaction timescale: 2982.4701279947594 Gyr\n",
      "Direct collision timescale: 29.824701279947593 Gyr\n",
      "==========================\n",
      "Velocity: 119.71354699169737 km/s\n",
      "Orbital timescale: 4.900698498564149 Gyr\n",
      "Hubble timescale: 13.968571428571428 Gyr\n",
      "Relaxation timescale: 1.209287679203387e+20 Gyr\n",
      "Close interaction timescale: 4.900698498564149e+22 Gyr\n",
      "Direct collision timescale: 4.2015590694137074e+36 Gyr\n",
      "==========================\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "class timescales:\n",
    "    def __init__(self, N, mass, radius, Rmin, H0):\n",
    "        self.gee = 4.2994e-9\n",
    "        \n",
    "        # Unitless\n",
    "        self.N = N\n",
    "        \n",
    "        # In Msun\n",
    "        self.mass = mass\n",
    "        \n",
    "        # In Mpc\n",
    "        self.radius = radius\n",
    "        \n",
    "        # In km/s/Mpc\n",
    "        self.H0 = H0\n",
    "        \n",
    "        # Mpc s/km in Gyr\n",
    "        self.Mpc_s_by_km_to_year = 9.778e2\n",
    "        \n",
    "        # Rmin\n",
    "        self.Rmin = Rmin\n",
    "\n",
    "        # velocity\n",
    "        self.v = np.sqrt(self.gee*self.mass/self.radius)\n",
    "        print(\"Velocity:\", self.v, \"km/s\")\n",
    "        \n",
    "        print(\"Orbital timescale:\", self.torb(), \"Gyr\")\n",
    "        print(\"Hubble timescale:\", self.tHubble(), \"Gyr\")\n",
    "        print(\"Relaxation timescale:\", self.trelax(), \"Gyr\")\n",
    "        print(\"Close interaction timescale:\", self.tclose(), \"Gyr\")\n",
    "        print(\"Direct collision timescale:\", self.tdirect(), \"Gyr\")\n",
    "        print (\"==========================\")\n",
    "\n",
    "    def tHubble(self):\n",
    "        # Return in Gyr \n",
    "        return 1/self.H0 * self.Mpc_s_by_km_to_year\n",
    "    \n",
    "    def torb(self):\n",
    "        # Return in Gyr \n",
    "        return 2*self.radius/self.v * self.Mpc_s_by_km_to_year\n",
    "    \n",
    "    def tdirect(self):\n",
    "        # Return in Gyr, assumes the radius to be Rsun\n",
    "        return self.torb()*(self.radius/self.Rmin)**2/self.N\n",
    "    \n",
    "    def tclose(self):\n",
    "        # Return in Gyr \n",
    "        return self.torb()*self.N\n",
    "    \n",
    "    def trelax(self):\n",
    "        # Return in Gyr \n",
    "        return self.torb()*self.N/8/np.log(self.N)\n",
    "\n",
    "Rsun_to_Mpc = 1./4.435e13\n",
    "mwbaryons = timescales(1e11, 1e11, 0.01, Rsun_to_Mpc, 70.0)\n",
    "globular = timescales(1e4, 1e4, 1e-5, Rsun_to_Mpc, 70.0)\n",
    "cluster = timescales(1e3, 1e14, 1, 0.010, 70.0)\n",
    "\n",
    "# Mass of dark matter particle assuming primordial black holes\n",
    "mdm = 1e-10\n",
    "    \n",
    "# Interaction cross section of about 1 cm^2/gm\n",
    "rdm = (1*mdm)**0.5*3.24e-25\n",
    "    \n",
    "mwdmhalo = timescales(1e22, 1e12, 0.3, rdm, 70.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the long range nature of gravity, one has to compute the potentials by summing over all the other stars. Since the number of stars in a galaxy are quite big, one can consider smooth potentials instead of discrete summation over all the stars.\n",
    "\n",
    "Our goal in galactic dynamics is to use kinematic observables to understand the density distribution in galaxies. The observables we have is the light coming from galaxies projected on the plane of the sky, and the line of sight velocity at pixelated locations of the galaxy. From these we need to construct the density distribution of the galaxy. The only way forward is to forward model the system. Make an ansatz about the density distribution in the galaxy, obtain the projected light distribution assuming a mass to light ratio. Obtain the potential corresponding to the density distribution. Construct orbits in this potential. Give appropriate weights to orbits and see if they can reconstruct the density distribution.\n",
    "\n",
    "This goal is quite hefty, so let us start with something very simple to begin with. Galaxies show a variety of symmetries, we will solve the Poisson equation for simpler symmetric density fields. We will start of with simple spherically symmetric mass distributions, and slowly work our way towards more complicated distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Density potential pairs:\n",
    "\n",
    "### Point mass distribution\n",
    "\n",
    "Let us consider the simple case of the point mass distribution. Let us choose our coordinate system such that the mass is located at the origin ($\\rho(\\vec{r}) = M\\delta(\\vec{r})$). We seek solutions of the Poisson equation for this density distribution.\n",
    "\n",
    "One can show that the function $\\nabla^2 |\\vec{r}|^{-1}$ is zero for all $|\\vec{r}|\\ne 0$. To compute its value at $\\vec{r}=0$, one can consider a small sphere centered at the origin with radius $R$. The integral\n",
    "\n",
    "\\begin{equation*}\n",
    "\\int_{V} \\nabla.\\nabla \\frac{1}{|\\vec{r}|} dV = \\int_{S} \\nabla \\frac{1}{|\\vec{r}|} . d\\vec{A} = -4 \\pi\n",
    "\\end{equation*}\n",
    "\n",
    "This implies we can write $\\nabla^2 |\\vec{r}|^{-1}= -4\\pi \\delta(\\vec{r})$, and therefore\n",
    "\n",
    "\\begin{equation*}\n",
    "\\Phi = -\\frac{GM}{|\\vec{r}|}\n",
    "\\end{equation*}\n",
    "\n",
    "is indeed a solution of the Poisson's equation for a delta function density. These two functions $\\left[M\\delta(\\vec{r}), -GM|\\vec{r}|^{-1}\\right]$ form a potential and density pair which satisfies the Poisson equation.\n",
    "\n",
    "The force is just the gradient of the potential, and thus is given by\n",
    "\n",
    "\\begin{equation*}\n",
    "F(\\vec{r}) = -\\nabla \\Phi = -\\frac{GM}{|\\vec{r}|^3}\\vec{r}\n",
    "\\end{equation*}\n",
    "\n",
    "It is also easy to see that if there are two potential density pairs then any linear combination of the two will also form potential density pair. This makes some of the analytical potential density pairs quite useful to study."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a generic mass distribution $\\rho(\\vec{r})$ we should be able to show that\n",
    "\n",
    "\\begin{equation*}\n",
    "\\Phi(\\vec{r}) = -\\int \\frac{G\\, \\rho(\\vec{r'})\\, d^3\\vec{r'}}{\\left|\\vec{r}-\\vec{r'}\\right|}\n",
    "\\end{equation*}\n",
    "\n",
    "A lot of things get simplified for spherically symmetric mass distribution which help us to quickly compute the potentials. The Newton's theorems make it easier to compute these quantities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Newton's theorems\n",
    "\n",
    "#### Newton's first theorem\n",
    "\n",
    "A spherical shell of matter does not exert any gravitational force on bodies that lie anywhere inside (the gravitational potential inside the shell of a given radius is constant).\n",
    "\n",
    "#### Newton's second theorem\n",
    "\n",
    "The gravitational force on a body lying outside a spherical shell of matter with mass $M$ is the same as the force exerted due to mass $M$ lying at the center of the shell.\n",
    "\n",
    "These theorems are easiest to prove using the Poisson equation integrated over suitable volumes. They also provide a useful way to compute the potentials when combined with the result for the point mass potential above. For a spherical shell of matter with radius $R$ centered at the origin, we can therefore write down,\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(\\vec{r}) = -\\frac{GM}{|\\vec{r}|} \\,\\,\\,\\,\\, \\forall \\,\\, r \\geq R\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(\\vec{r}) = -\\frac{GM}{R} \\,\\,\\,\\,\\, \\forall \\,\\, r \\leq R\n",
    "\\end{equation}\n",
    "\n",
    "We can plot the gravitational potential due to such a shell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import pylab as pl\n",
    "%pylab inline\n",
    "conf = %config InlineBackend.rc\n",
    "conf['figure.figsize'] = (6, 6)\n",
    "conf['savefig.dpi']=100\n",
    "conf['font.serif'] = \"Computer Modern\"\n",
    "conf['font.sans-serif'] = \"Computer Modern\"\n",
    "conf['text.usetex']=True\n",
    "\n",
    "%config InlineBackend.rc\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f6c0f9c2490>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAFnCAYAAACmbT7/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU9b3/8fc3G2EPCWEnhLDvEAYoYrUX41a1t62pVevS2hZv7722VW+xtbXt77b3+pP+apfbX1u5Wm/V1qIIatWqRKsVi0IS9jUQkpCwZZtAAlnne/9gqDEFMsksZ87M6/l48GDmO8PwOcmZec/3+z3ne4y1VgAAJDhdAAAgOhAIAABJBAIAwI9AAABIIhAAAH5JThcQjKFDh9rs7GynywAAVykqKqqx1mZ2bXd1IGRnZ6uwsNDpMgDAVYwx5edqZ8gIACCJQAAA+IU1EIwxyzvdzjfG5PWmDQAQfmELBGNMnqQF/tu5kmStLZDkNcbkBtoWrvoAAB8WqSGjz0ry+m+XSsrrQRsAIALCEgjGmFz/t/yz0iTVdbqf0YO2rq+9zBhTaIwprK6uDmHVABDfwtVDSA/T68pau9Ja67HWejIz/+4wWgBAL/XqPARjzLJzNJdaawvO0TuQzgwDnQ2JNEm1/tuBtgEAwqxXgWCtXXmBh3OMMTmdbudKWiXJc7ZN0tnACLQNABBmIR8ystauttau1plv+mn+tmLpb0ceea21xYG2hbo+AMC5GTdfMc3j8ViWrgCAnjHGFFlrPV3bOVMZACCJQAAA1zjd2qEXtlTpS7/dpJPNbSF/fVevdgoAsc7ns3rvYK2eK6rSqzuOqKm1Q6PT+qq89pRmjh4c0v+LQACAKFRW06Tniiu1prhKVd7TGtAnSdfMHqlPzRujRePTlZBgQv5/EggAECWaWtr18vYjWl1YqY1ldUow0pKJQ7X8qim6YvoI9U1JDOv/TyAAgIOstSqu8OqZTYf00rbDamrt0Pih/fWNK6fo+twxGjE4NWK1EAgA4ID6plat2VylVZsqtO9Yo/omJ+ra2SN1w4Kx8owbImNCPyTUHQIBACLEWqtNZfX6/fvlemXHUbW2+zRnbJoe/PQsXTdnlAb0cfYjmUAAgDA70dymNUWV+t37FSo53qiBfZL0Wc9Y3bQwS9NHDXK6vL8hEAAgTHYebtBT75Xr+c2HdbqtQ3PGDNaK62fr2jkj1S8l+j5+o68iAHCx1nafXt15VE/8tUyF5fVKTU7QP84ZrVs+Mk6zxoT2vIFQIxAAIASqT7bo6Y0Veuq9ch0/2aJxGf30nWum6TPzx2pwv2SnywsIgQAAQdh5uEG/WV+mP249rNYOny6dnKmHrs/WpZMzw3LyWDgRCADQQz6f1Rt7juux9aV6r7RO/VISdePCsbr9omxNyBzgdHm9RiAAQICa2zq0uqhSv1l/UKU1TRo1OFXfunqqblyYpcF93TEsdCEEAgB0o76pVU9sKNdvN5SprqlVs8cM1s9vmqePzxyhpMTYWTSaQACA86isP6VH3zmoVZsO6XRbhy6bOkzLLsnRwvHpjpxJHG4EAgB0UXLspH711gG9sPWwjKRPzhutZZfkaPLwgU6XFlYEAgD4ba9s0C/+XKLXdh5T3+RE3b44W1/66HiNSuvrdGkRQSAAiHtF5fX6rzdL9Nbeag1KTdJXl07UF5aM15D+KU6XFlEEAoC4tamsTj8rKNH6/TVK75+ib1w5RbcuHqdBqe4/Yqg3CAQAcaewrE4/Kdind/fXauiAFH3749P0uY9kReX6QpEU31sPIK5sOeTVj1/fq3dKajR0QIq+c800fW7RuLBficwtCAQAMW/3kRP68ev7VLD7mNL7p+j+j0/VLR8ZF/c9gq74aQCIWWU1TXp43T79cdthDeiTpH+7YrI+v2S84xeiiVb8VADEnOMnmvXzN0v0h42HlJyYoK9cOkF3XjLBNauOOoVAABAzGlvatfLtA/rvdw6qrcOnmxZm6a6lEzVsUOQuVO9mBAIA12vv8OnpjRX6aUGJaptadc3skfrGFVOUPbS/06W5CoEAwLWstXpzz3H95yu7daC6SQvHp+uxj0/T3LFpTpfmSgQCAFfac/SEfvjSbq3fX6Ocof218tb5unz68JhcdC5SCAQArlLb2KKH1+3T0xsrNKhvsr533XTd8pFxSo6hZaidQiAAcIX2Dp+efK9cD6/bp1OtHbptcba+njdJaf3ia72hcCIQAES9DQdq9f0Xd2rvsZP66KSh+t510zVxWGwvRe0EAgFA1Dp2oln/8fJuvbj1sMYM6atHbp2vK5gnCBsCAUDUae/w6X/+WqafFpSotcOnr+dN0j9dOkGpyaw5FE4EAoCosrmiXvev3aHdR07oH6Zk6vufmKFxGZxPEAkEAoCo0HC6TSte3aPfb6zQ8IGp+vUtubpyxgiGhyKIQADguFd3HNF3X9ipmsYWff6ibN17xRQWoHMAP3EAjjl2olnffWGHXtt5TNNHDtKjt3s0ewxnGTuFQAAQcdZaPVtYqR+8vEut7T598+qp+uLF4zm5zGEEAoCIqqw/pW+t2a53Smq0cHy6Hrp+tsazCF1UIBAARIS1Vqs2HdIPX94ta61+8MmZ+tzCLCUkMGkcLQgEAGF3tKFZ9z23TW/vq9binAytyJ+tsen9nC4LXRAIAMLqxa2H9cDzO9Ta7tO//+MM3bJoHL2CKBWWQDDG5ErKkSRr7Wp/W74kr6Rca+2KnrQBcB/vqVZ95/kdemnbEc3LStPDN8xlriDKhauHcKe19k5jzHJjTI6kXEmy1hYYY3KMMXmS0gJps9YWdH5hY8wyScskKSsrK0zlAwjGXw/U6J5VW1XT2KJvXDlFd16SoySOIIp6If8N+T+wi4wxOdbaFdbaUkkLJJX6n1KqMwERaNuHWGtXWms91lpPZmZmqMsHEITWdp8efGW3Pvfo++qXkqi1/7xE//IPEwkDlwhHD2GC/+9njDGPSLpP/m/+nWT0oA2ACxysadJdTxdrR9UJ3bwoS9+5Zpr6pTBN6Sa9+m35ewFdlXYa3jlgrfUaY4p0ZnjHKym9y/MDbQMQ5dZurtR31u5QUmKCfn3LfF01c4TTJaEXehUI1tqVF3h4kz74UE/TmQ/5Un3w7T9H0jr//UDaAESpU63t+u4LO7W6qFILsofopzfO0+i0vk6XhV4KeX/OWrvaP5mc57+/UpK6tBX0pA1A9Ck5dlL//Lti7a9u1F1LJ+prl01irsDljLXW6Rp6zePx2MLCQqfLAOLO2s2Vun/NDvVLSdRPPjtXl0zmAA83McYUWWs9XduZ8QEQsJb2Dv3gpV166r0KLcxO13/dPE/DB6U6XRZChEAAEJDD3tP6yu+KtfWQV8suydHyK6cwRBRjCAQA3dpwoFb/8vtitbb79OtbcnXVzJFOl4QwIBAAnJe1Vo+/W6b/eGW3sjP6aeVtHk3IHOB0WQgTAgHAOTW3dej+tdu1prhKl08frodvmKOBqclOl4UwIhAA/J3jJ5r15SeLtPWQV3fnTdZdSyeyQmkcIBAAfMi2Sq++/EShTpxuZ74gzhAIAP7mle1HdM8zW5TRv4+e+8pFmj5qkNMlIYIIBACy1uqXbx3Qj17bq9ysNK28zaOhA/o4XRYijEAA4lxru0/fXrtdzxZV6hNzRmlF/mylJic6XRYcQCAAcexkc5u+8lSx1u+v0Vcvm6S78ybJGCaP4xWBAMSpIw2n9YXHN2n/8Ub9KH+2PuMZ63RJcBiBAMShfcdO6rbHNqqxpV2Pf2GBPjqJxelAIABxZ1NZnb74P5uUmpyoZ+5czJFE+BsCAYgj63Yd07/+vlij0/rqt3cs1Nj0fk6XhChCIABxYnVRpZav3qpZY9L0m9s9yuCwUnRBIABx4DfrD+rfX9qliycO1SO3zlf/Prz18ffYK4AYZq3VTwtK9LM3SnTVjBH62U1z1SeJcwxwbgQCEKOstfqPl3fr0fUHlT9/jP7vp2dxQRtcEIEAxCCfz+p7L+7Uk++V67bF4/T962awWim6RSAAMabDZ3X/mu1aVXhIyy7J0beunsrZxwgIgQDEkA6f1fLV2/RccaW+unSi7r58MmGAgBEIQIzo8Fl9Y/VWrSmu0t15k/W1vElOlwSXIRCAGNA5DO65fLK+ehlhgJ7jkAPA5Xw+q28+t40wQNAIBMDFrLV64IUderaoUl+7bBJhgKAQCIBLWWv17y/t0u/er9BXPjZBX2fOAEEiEACX+vHr+/T4u2W6Y8l4Lb9yCkcTIWgEAuBCj7x9QL/4837dtHCsHrh2GmGAkCAQAJd5emOFHvzTHl07e6R++MlZhAFChkAAXOTlbUd0/9rt+tiUTD18w1wlshwFQohAAFzirwdqdPeqLZqfNUS/+tx8pSTx9kVosUcBLrCjqkHLnihS9tB+evR2j/qmsIQ1Qo9AAKJcRe0pff7xTRqUmqTf3rFQaf1SnC4JMYpAAKJYfVOrPv/4RrV1+PTEFxdq5OC+TpeEGMZaRkCUam7r0LInC1XpPa3ffWmRJg4b6HRJiHH0EIAo5PNZ3fvsVm0qq9ePPzNHC7LTnS4JcYBAAKLQj17fq5e3HdE3r56q6+aMcrocxAkCAYgyzxYe0q/eOqCbF2XpzktynC4HcYRAAKLIxoN1un/tdi2ZmKH/84kZnIWMiCIQgChRUXtKdz5ZqLFD+umXN89XciJvT0QWexwQBRpb2vWlJzbJZ6XHPr9Ag/slO10S4hCHnQIO8/ms7lm1RQeqm/TEHQs1fmh/p0tCnApLIBhj8iV5JeVYa1d2acu11q7oSRsQy37+Zole33VM3712upZMHOp0OYhjIQ8EY0yepFJrbbExJs8YkyspR5KstQXGmBz/c9ICabPWFnR5/WWSlklSVlZWqMsHIuq1nUf104IS5c8foy8syXa6HMS5cMwhFEp69mwQWGuLJS2QVOp/vFRSbg/aPsRau9Ja67HWejIzM8NQPhAZ+4836t5ntmrO2DT98JMzOaIIjgt5IFhrvZIekfSspAn+5rQuT8voQRsQc5pa2vVPTxWpT1KCfn1LrlKTWb0UzuvVkJF/2KarUv9QT76kAmvtCmPMQ53mBLqeex9oGxBTrLW677ltKq1u1FNfXMSCdYgavQqEsxPF55FjrV3tv/2gpBskbdIH3/5zJK3z3w+kDYgpv3m3TC/5l6W4iElkRJFwHGW00t+DKNWHjzJa7p841tmJ4kDbgFhRXFGvB1/ZrStnDGdZCkQdY611uoZe83g8trCw0OkygIB4T7Xq4z97R0mJCXrpqxdrUConn8EZxpgia62nazsnpgER4PNZ3fvMVtU0tuq5r1xEGCAqsXQFEAGPri/VG3uO69vXTNOsMYOdLgc4JwIBCLMth7xa8epeXTVjhG5bPM7pcoDzIhCAMGpsadfX/rBZwwel6qHrZ3PyGaIacwhAGH33+R06VHdKq+5czAqmiHr0EIAweX5zldZsrtJdSydxTWS4AoEAhEFl/Sk98PwOecYN0V1LJzpdDhAQAgEIsQ6f1T3PbJWV9JPPzlUSVz6DS7CnAiH22PpSbTxYp+9dN11j0/s5XQ4QMAIBCKE9R0/o/722T1dMH678+WOcLgfoEQIBCJHWdp++/octGtQ3SQ9+ehaHmMJ1OOwUCJFfvFmiPUdP6tHbPMoY0MfpcoAeo4cAhMCOqgb9/7cO6NPzRitv+nCnywF6hUAAgtTa7tO/PbtVGf1T9L3rZjhdDtBrDBkBQeo8VMTZyHAzeghAEHYdPqFfMlSEGEEgAL3U4bP65pptSuuXrAeune50OUDQCASglx5/96C2VTboe9fN0JD+KU6XAwSNQAB64VDdKf349X1aOnWYrp090ulygJAgEIAestbq28/vUIKRfvDJmZyAhphBIAA99OLWw/rLvmp948opGp3W1+lygJAhEIAeONHcph++vFuzxwzWrYuznS4HCCnOQwB64OHX96mmsUWP3e5RYgJDRYgt9BCAAO2oatATG8p060fGafaYNKfLAUKOQAAC0OGz+vba7Urv30f3XjHF6XKAsCAQgAA8U3hIWysb9MC10zS4L8tTIDYRCEA3Gk616Uev7dXC8en6xJxRTpcDhA2BAHTjJwX75D3Vqu9fN4NzDhDTCATgAvYePakn3yvX5xaN0/RRg5wuBwgrAgE4D2utvv/iTg1MTdI9l092uhwg7AgE4Dxe23lUG0prde8VU1i8DnGBQADOoaW9Qw/+aY+mDB+omxaMdbocICIIBOAcntxQrvLaU7r/mmlKSuRtgvjAng50Ud/Uqp+/UaJLJ2fq0smZTpcDRAyBAHTxszdK1NjSrm9fM83pUoCIIhCATkqrG/XUe+W6cWGWJg8f6HQ5QEQRCEAnP3ptr/okJejuPA4zRfwhEAC/LYe8+tOOo/ryJTnKHNjH6XKAiCMQAJ05Ce2hP+1RRv8UfemjOU6XAziCQAAk/aWkRhtKa3XX0oka0IfrRiE+EQiIez7fmd7B2PS+unnROKfLARwTkkAwxuR2uZ9vjMkzxiwPRRsQTi9tP6JdR07o3sunKCWJ70iIX0Hv/caYPEnPdrqfK0nW2gJJXmNMbjBtwdYHXEh7h08/WbdPU0cM5FoHiHtBB4L/w7u0U9NnJXn9t0sl5QXZBoTN2s1VOljTpLsvn6yEBK51gPgWjv5xmqS6Tvczgmz7EGPMMmNMoTGmsLq6OmRFI/60dfj08zdLNGv0YF0xfbjT5QCOc92AqbV2pbXWY631ZGayzgx679nCSh2qO617Lp/MldAASd0eX2eMWXaO5lL/UNG5eCWl+2+nSar13w6mDQiplvYO/eLNEs3LStPHpvDFApACCARr7coevuYqSR7/7RxJZ4MjmDYgpFZtOqTDDc1akT+H3gHgF4qjjPIlefx/y1pb7G/Pk+S11hYH0xZsfUBXLe0d+tVbB7Qge4iWTPy7aSogbgV9Sqa1drWk1V3a/q5XEUwbEErPFVXpSEOzHrp+Nr0DoBPXTSoDwWjr8OmXb+3XnLFp+uikoU6XA0QVAgFx5fnNVaqsP62vLp1I7wDogkBA3OjwWf3yrQOaPnKQlk4d5nQ5QNQhEBA3Xtp2WAdrmnQXvQPgnAgExAVrrX711gFNGjZAV84Y4XQ5QFQiEBAX3tpbrT1HT+qfLp3AmkXAeRAIiAu/evuARg1O1SfmsqIpcD4EAmJecUW9Nh6s0x0Xj1dyIrs8cD68OxDzHnn7gAb3TdZNC7OcLgWIagQCYtqB6ka9vuuYbls8Tv25VjJwQQQCYtp//6VUKYkJuv2ibKdLAaIegYCYVdPYojWbq5Q/f4yGDujjdDlA1CMQELOeeq9cre0+3XHxeKdLAVyBQEBMam7r0FPvlesfpmRqQuYAp8sBXIFAQEx6ceth1TS26osX5zhdCuAaBAJijrVWv1l/UFNHDOQCOEAPEAiIOX89UKs9R0/qjiXjWcQO6AECATHnsfUHNXRACstUAD1EICCmVNSe0p/3HtfNC7OUmpzodDmAqxAIiClPvV+uBGN086JxTpcCuA6BgJhxurVDqzYd0lUzRmjE4FSnywFch0BAzPjj1sNqON2mWxfTOwB6g0BATLDW6rcbyjR5+AAtGp/udDmAKxEIiAnFFV7tPHxCty3O5lBToJcIBMSEJzeUaWCfJH1q3minSwFci0CA69U1teqV7Ud1/fwxXPMACAKBANdbU1yp1g4fV0QDgkQgwNWstfr9xgrNHzdEU0YMdLocwNUIBLjaxoN1Kq1uoncAhACBAFd7emOFBqYm6ZpZI50uBXA9AgGuVd/Uqld2HNWn541W3xTWLQKCRSDAtdZsrlJru083MlwEhASBAFey1urpjRWaOzZN00YOcrocICYQCHClLYe82n+8UTcuGOt0KUDMIBDgSs8WVSo1OUHXzGYyGQgVAgGu09zWoT9uPayrZ47UwNRkp8sBYgaBANd5bedRnWxu12fmj3G6FCCmEAhwndVFlRqd1lcfyclwuhQgphAIcJXD3tNav79G188fo4QElrkGQolAgKus3Vwla6X8XIaLgFAjEOAa1lqtLqrUovHpysro53Q5QMwhEOAamw95dbCmSdfTOwDCgkCAa7ywuUopSQm6atYIp0sBYhKBAFdo6/DppW1HlDdtmAZx7gEQFkEHgjEmzxizrtP9NGNMvv/PQ53a8/3PXd7TNmB9SY1qm1r1yblcMxkIl6ADwVpb0KXpBknp1trVkmSMWWaMye/0XK//Az+gtq7/n//1Co0xhdXV1cGWD5d4fkuVBvdN1semDHO6FCBmhXzIyFq70lq70n83R1KBpAWSSv1tpZJye9B2rtf3WGs9mZmZoS4fUaippV2v7zymj88aqZQkRjmBcAnbu8sYkyOpzlpbKimty8MZPWhDnFu365hOt3XoU/MYLgLCKam7Jxhjlp2jufQcQ0Vd5Vtr7/Tf9kpK7/J4oG2Ic2s3V2l0Wl95xg1xuhQgpnUbCJ2GfwJmjMm31q7w386TtEkffPvPkbTOfz+QNsSx2sYWrd9fo2WX5LBUBRBmoTjKKF+S5+yEsD8AHjLGFBljiiTJP8Gcc3aS2FpbEGhbsPXB3f6046g6fFafmDPK6VKAmNdtD6E7/g/x1Z3uF0iacI7nrehtG+LXy9uOKCezv6aOGOh0KUDM45ANRK3jJ5v1/sFaXTtrpIxhuAgINwIBUeu1HUfls9K1DBcBEUEgIGr9cdsRTRo2QJOHM1wERAKBgKh07ESzNpXV6ZrZI50uBYgbBAKi0p+2H5G10rUEAhAxBAKi0svbj2jqiIGaOIzhIiBSCAREnaMNzdpUVq9rZtE7ACKJQEDUeX3XUUnS1VwIB4goAgFR59UdRzUhsz/DRUCEEQiIKvVNrXr/YJ2unEHvAIg0AgFR5Y09x9Xhs7pqJoEARBqBgKjy6o6jGjU4VbNGD3a6FCDuEAiIGqda2/VOSbWumDGCtYsABxAIiBpv761WS7uP+QPAIQQCosarO49qSL9kLcjmymiAEwgERIXWdp/e3HNcedOGKymR3RJwAu88RIX3Smt1srmd4SLAQQQCosIbu48pNTlBF08a6nQpQNwiEOA4a63e2HNcSyYMVWpyotPlAHGLQIDj9h1rVGX9aV02bbjTpQBxjUCA497Yc0yStHTqMIcrAeIbgQDHvbH7uGaOHqQRg1OdLgWIawQCHFXX1KriinpdNpXhIsBpBAIc9ec9x2WtdNk0hosApxEIcNSbe45r2MA+mjmKxewApxEIcExru09v76vW0qnDlJDAYnaA0wgEOGZTWZ0aW9o53BSIEgQCHPP2vmqlJCboogkZTpcCQAQCHPT23motGD9E/fskOV0KABEIcMiRhtPae+ykLp2c6XQpAPwIBDjiL/uqJUmXTuZwUyBaEAhwxNv7qjViUKomDx/gdCkA/AgERFx7h0/vlNTo0smZXDsZiCIEAiJuyyGvTja369IpzB8A0YRAQMS9va9aiQlGSyZyMRwgmhAIiLi391Vr3tg0De6b7HQpADohEBBRtY0t2l7VwOGmQBQiEBBR6/fXyFrpEgIBiDoEAiJqfUmN0vola+ZoVjcFog2BgIix1urd/TVanJOhRFY3BaIOgYCIKas9pcMNzRxdBEQpAgER8+7+GkkiEIAoRSAgYt7dX6NRg1OVndHP6VIAnEPQgWCMyTPGrDvPYw91up3vf+7ynrbB/Tp8VhtKa7Vk4lCWqwCiVNCBYK0tOFe7MSZPUo7/dn6n53r9H/gBtZ3jdZcZYwqNMYXV1dXBlo8I2XX4hLyn2hguAqJYWIaMjDE5kko7NS3odL9UUm4P2j7EWrvSWuux1noyMzmW3S3ePXBm/uCiiVwdDYhW4ZpDyLHWdg6EtC6PZ/SgDTHg3f01mjx8gIYNTHW6FADn0e21C40xy87RXHqhoaJzPOaVlN7LNrhcS3uHNpXV6cYFWU6XAuACug0Ea+3KHr5mnX/sP01SjjEmV9ImffDtP0fSOv/9QNrgcsXlXjW3+XQx8wdAVAvFUUb5kjydJoSL/T2EdPk/3K21q3UmHPL89wsCbQu2Pjhvw4EaJRhpYQ6dPyCaGWut0zX0msfjsYWFhU6XgW7c8MgGNbd16MV/vdjpUgBIMsYUWWs9Xds5MQ1h1dzWoS0VXi0aT+8AiHYEAsJqc4VXrR0+LRrPAWNAtCMQEFbvH6yVMdICeghA1CMQEFbvl9Zp2ohBXC4TcAECAWHT0t6h4op6LeLoIsAVCASEzbbKBrW0M38AuAWBgLB5v7RWkrSQ+QPAFQgEhM37B+s0ZfhApfdPcboUAAEgEBAWbR0+FZUzfwC4CYGAsNhe1aBTrR3MHwAuQiAgLN4vrZPE/AHgJgQCwqKwrE45Q/src2Afp0sBECACASHn81kVVdRr/rghTpcCoAcIBIRcaU2TvKfa5MkmEAA3IRAQckXlZ+YP5o9j/gBwEwIBIVdYVq8h/ZI1IbO/06UA6AECASFXVH5m/sAY43QpAHqAQEBI1Ta2qLSmieEiwIUIBIRUcYVXkphQBlyIQEBIFZbXKSUxQbNGD3a6FAA9RCAgpIrK6jVz9CClJic6XQqAHiIQEDIt7R3aVtUgTzbzB4AbEQgImR1VDWpt93GGMuBSBAJCpqi8XpKUm0UgAG5EICBkisu9GpfRjwXtAJciEBAyWw55NW9smtNlAOglAgEhcbShWUdPNGsOgQC4FoGAkNhy6Mz8wVwCAXAtAgEhsfmQVymJCZo+apDTpQDoJQIBIbGlwqtpowapTxInpAFuRSAgaB0+q+1VDZo7huUqADcjEBC0kuMndaq1Q3OzmD8A3IxAQNC2+Fc4nTuWE9IANyMQELQth7wa3DdZ2Rn9nC4FQBAIBARtyyGv5oxN4wppgMsRCAhKU0u79h07yfkHQAwgEBCUbZUN8lmxZAUQAwgEBGVr5ZkJZZasANyPQECvtbR36MUthzV+aH+l909xuhwAQSIQ0Gv/+fJu7TpyQt+6eqrTpQAIAQIBvfLStsP67YZyfeni8bpixginywEQAklOF+CE6pMtOt3a4XQZrlXd2KJvPrdduVlpuo/eARAz4jIQHnh+h17defwcdPAAAARYSURBVNTpMlxtSL9k/eLmXCUn0skEYkXQgWCMyZN0n7X28k5tuZJyJMlau9rfli/JKynXWruiJ22hdvtF2bp8+vBwvHTc8GQP0ai0vk6XASCEgg4Ea22BMea+Ls13WmvvNMYsN8bkSMrt9Nwcf4ikBdJmrS3o/MLGmGWSlklSVlZWr2pePCGjV/8OAGJZyPv7/g/sImNMjrV2hbW2VNICSaX+p5TqTEAE2vYh1tqV1lqPtdaTmZkZ6vIBIG6FYwB4gv9PnTHmEWNMmvzf/DvJ6EEbACACuh0y8n/j76q061BOFwestV5jTJHODO94JaV3eU6gbQCACOg2EKy1K3v4mpv0wYd6ms58yJfqg2//OZLW+e8H0gYAiICgh4z8RwV5/H+fPaoozT9JfHbMf7WknE5tBYG2BVsfACAwxlrrdA295vF4bGFhodNlAICrGGOKrLWeru2cVQQAkEQgAAD8CAQAgCQCAQDgRyAAACQRCAAAP1cfdmqMqZZU3st/PlRSTQjLcQO2OT6wzbEv2O0dZ639u8XgXB0IwTDGFJ7rONxYxjbHB7Y59oVrexkyAgBIIhAAAH4EAgBAEoEAAPAjEAAAkuI7EHp6nYdYwDbHB7Y59oVle+P2sFMAwIfFcw8BANAJgYCYYIzJN8bkGWOWd/O8Cz4ORCNjTO4FHgto3w9EXARCdz+wUP5Ao0UA27zM/+ehSNcWamffLP5LrnrP9+bxX5r18kjWFk4B/I5z/c/Jj3Rt4dKD9/KySNcWLv799tnzPBbQvh+omA+E7n5gof6BRoMAtjlPUoG1dqU6XcPaxT4ryeu/XSrJ7dvTrQD32291uk55POzXuZJK/Y+XxsI2S3/b3tLzPBzSfT/mA0Hd/8Bi8cOku23K6dRW6r/vZmmS6jrdz+j6BGNMrv+NFSsu+Dv29wo2SZK1doW1tjiy5YVFIO/Vsz3enBjZ5u50u+/3RDwEQnc/sJD+QKPEBbfJWrvS3zuQpFxJhZEqzEHpThcQYt3ttwskZfiHjWJlKLS7/bpYZ3oG9V2ehwDFQyDgPPxd6uIY+Cbl1Qcf+GmSajs/GIO9g0DVnv3dxtI8wvkYY9J0Zl94UNJ/G2Pc3vMNxAX3/Z6Kh0Do7gcW0h9olAh0m/KstfdFpqSwWqUPhr1yJBVIf/uAkM6Moef7JxrTY2Rsubvfca0+GHf26kyPwe262+Zlkh601q6Q9GVJMRuCnfbtc+77vRUPgdDdh0VIf6BRorttljFmmf+NI7dPKnf6Fpwnydupx/OG//HV/slV6cwHSSzo7ne8utPjafLPJ7hct/v1Wf7ft7druxv5e3eeLr28s/v2+fb93v1f8XCmsv+bYanOTDSt9LcVWWvnn+9xt7vQNnc6jK1OZ75xfSZOh1RcLcD9uk7SghjpCQayzcv9j6fHyns5kuIiEAAA3YuHISMAQAAIBACAJAIBAOBHIAAAJBEIAAA/AgEAIIlAAAD4/S+DzGY29aUjMAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "class shell:\n",
    "    def __init__(self, M, R):\n",
    "        self.Rshell = R\n",
    "        self.Mass = M\n",
    "        # Put gravitational constant * mass of the sun /(km/s)^2 /Mpc in google search!\n",
    "        # What are the units of this constant?\n",
    "        self.gee = 4.3022e-9\n",
    "\n",
    "    def potential(self, r):\n",
    "        if np.isscalar(r):\n",
    "            r = np.array([r])\n",
    "        pot = r * 0.0\n",
    "        idx = r<self.Rshell\n",
    "\n",
    "        pot[idx] = -1/self.Rshell\n",
    "\n",
    "        pot[~idx] = -1/r[~idx]\n",
    "\n",
    "        return pot*self.gee*self.Mass\n",
    "\n",
    "# Mass in hinv Msun, shell radius in hinv Mpc \n",
    "Mass = 1e12 \n",
    "sh = shell(1e12, 0.270)\n",
    "R = np.linspace(0.0, 1.0, 100)\n",
    "\n",
    "plot(R, sh.potential(R))\n",
    "\n",
    "# Label the axis of these plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A spherically symmetric mass distribution can be decomposed in to smaller shells, and the potential can then be added up. Thus for a density distribution which depends only on the radius, $\\rho(r)$, we obtain\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(\\vec{r}) = -4\\pi G \\left[\\int_r^{\\infty} \\rho(r') r' dr' + \\frac{1}{r}\\int_0^{r} \\rho(r') r'^2 dr'\\right]\n",
    "\\end{equation}\n",
    "\n",
    "This can be shown to be equal to\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(\\vec{r}) = -G \\left[\\int_{r}^{\\infty} \\frac{1}{r'^2} M(\\lt r') dr'\\right]\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Circular velocity\n",
    "\n",
    "Given the potential, one can calculate the velocity of circular orbits\n",
    "\n",
    "\\begin{equation}\n",
    "v_{\\rm c}(r) = \\sqrt{\\frac{GM(\\lt r)}{r}}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "One can also compute a dynamical time estimate, this is the time required to cross the system. It can be calculated as\n",
    "\\begin{equation}\n",
    "t_{\\rm dyn} = \\frac{2\\pi r}{v_c} \\propto \\frac{1}{\\sqrt{G\\rho}}\n",
    "\\end{equation}\n",
    "There are a number of ways to define dynamical time but all of these will come out to be proportional to\n",
    "$1/\\sqrt{G\\rho}$.\n",
    "\n",
    "\n",
    "1. Compute the circular velocity of gas particles in a Milky way halo. Assume mass of the Milky way is $10^{12} h^{-1}M_{\\odot}$ halo and compute the velocity at $0.1$ Mpc from the center which is about 10 times further than the galactocentric distance of the Sun.\n",
    "\n",
    "2. Compute the dynamical time at the same distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spherically symmetric distributions\n",
    "\n",
    "#### Uniform sphere\n",
    "\n",
    "For a uniform homogeneous sphere with constant density, $\\rho_0$ and radius $R$,\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\rho(r) = \\rho_0 \\,\\,\\,\\,\\forall\\,\\,r\\lt R \\\n",
    "\\rho(r) = 0 \\,\\,\\,\\,\\forall\\,\\,r\\gt R\\,. \\,\n",
    "\\end{eqnarray}\n",
    "\n",
    "\n",
    "The potential is given by\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(r) = - \\frac{4\\pi G\\rho}{3}\\frac{ R^3}{r} \\,\\,\\,\\,\\forall\\,\\,r\\gt R \\\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "\\Phi(r) = - \\frac{4\\pi G\\rho}{3}\\frac{ 3R^2-r^2}{2}  \\,\\,\\,\\,\\forall\\,\\,r \\leq R\\,. \\,\n",
    "\\end{equation}\n",
    "\n",
    "The circular velocity will increase linearly with the radius for $r \\leq R$,\n",
    "\\begin{equation}\n",
    "v_c(r) = \\sqrt{\\frac{4\\pi G \\rho}{3}} r \\,\\,\\,\\,\\forall\\,\\, r\\leq R\\,.\n",
    "\\end{equation}\n",
    "and it falls off as $r^{-1/2}$ for $r\\gt R$. If there is a disk of test particles moving in a uniform density sphere, then this disk will rotate with a constant angular speed $\\omega=v_{\\rm c}/r=$constant like a solid body.\n",
    "\n",
    "The dynamical time on the other hand for the system will be constant in these inner regions where the density is uniform,\n",
    "\\begin{equation}\n",
    "t_{\\rm dyn}(r) = \\sqrt{\\frac{3\\pi}{G \\rho}} \\,\\,\\,\\,\\forall\\,\\, r\\leq R\\,.\n",
    "\\end{equation}\n",
    "Thus the orbital time period of particles are also independent of the radius in regions of uniform density.\n",
    "\n",
    "Next let us consider a couple of potentials which smoothly interpolate between the point mass Kepler potential and the homogeneous sphere - Plummer sphere and the isochrone."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plummer sphere\n",
    "\n",
    "Let us consider the Plummer sphere where the potential is smoothed in the inner regions compared to the point mass potential.\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(r) = - \\frac{GM}{\\sqrt{r^2+b^2}} \\\n",
    "\\end{equation}\n",
    "\n",
    "The Plummer model is often used to compute the potentials of particles in numerical simulations. The smoothing at small scales avoids large forces during close chance interactions. These happen in numerical simulations due to finite mass resolution, but such collisions are expected to play a very insignificant role for large astrophysical systems like galaxies (globular clusters are different).\n",
    "\n",
    "The density corresponding to this potential can be found by using the Poisson equation and it corresponds to\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho(r) = \\frac{3M b^2}{4 \\pi}\\frac{1}{[r^2+b^2]^{5/2}}\n",
    "\\end{equation}\n",
    "\n",
    "The circular velocity will be given by\n",
    "\\begin{equation}\n",
    "v_{c}(r) = \\sqrt{\\frac{GM}{(b^2+r^2)^{3/2}}}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Isochrone potential\n",
    "\n",
    "Yet another useful potential is the isochrone potential where,\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(r) = - \\frac{GM}{b+\\sqrt{r^2+b^2}} \\\n",
    "\\end{equation}\n",
    "\n",
    "By taking the limits $b\\rightarrow 0$, or $r\\rightarrow 0$, one can smoothly interpolate between a point mass potential and the potential of a homogeneous sphere.\n",
    "\n",
    "- Compute the behaviour of the density distribution corresponding to the isochrone at small radii $r\\lt b$.\n",
    "- Compute the behaviour of the density distribution corresponding to the isochrone potential at large radii $r\\gt b$.\n",
    "\n",
    "The density associated with the potential can be shown to be\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho(r) = M \\frac{3(b+a)a^2-r^2(b+3a)}{4\\pi (b+a)^3a^3}\n",
    "\\end{equation}\n",
    "\n",
    "where $a=\\sqrt{b^2+r^2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Power law density profiles\n",
    "\n",
    "Many galaxies have density profiles which fall off as $r^{-\\alpha}$. The spherically symmetric version for such a system would be,\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho(r) = \\rho_0 \\left(\\frac{r_0}{r}\\right)^\\alpha\n",
    "\\end{equation}\n",
    "\n",
    "The mass within the radius $r$ is given by\n",
    "\n",
    "\\begin{equation}\n",
    "M(\\lt r) = \\int_0^r \\rho(r') 4 \\pi r'^2 dr' = \\frac{4\\pi \\rho_0 r_0^{3}}{(3-\\alpha)} \\left(\\frac{r_0}{r}\\right)^{\\alpha-3}\\,.\n",
    "\\end{equation}\n",
    "\n",
    "The mass outside the radius $r$ is given by\n",
    "\n",
    "\\begin{equation}\n",
    "M(\\gt r) = \\int_r^\\infty \\rho(r') 4 \\pi r'^2 dr' = \\frac{4\\pi \\rho_0 r_0^{3}}{(\\alpha-3)} \\left(\\frac{r_0}{r}\\right)^{\\alpha-3}\\,.\n",
    "\\end{equation}\n",
    "\n",
    "You will see that depending upon the magnitude of $\\alpha$ compared to 3, you either have a infinite enclosed mass or a infinite mass beyond radius $r$. Clearly such cases are not possible. One needs to have $\\alpha \\lt 3$ in the inner regions and $\\alpha \\gt 3$ in the outskirts.\n",
    "\n",
    "Assuming $\\alpha\\lt 3$ to avoid divergence of mass at the center, the circular velocity curve is\n",
    "\n",
    "\\begin{equation}\n",
    "v_{\\rm circ}^2(r) = \\frac{G 4\\pi \\rho_0 r_0^{3}}{(3-\\alpha)} \\left(\\frac{r_0}{r}\\right)^{\\alpha-3}\\frac{1}{r}\n",
    "\\end{equation}\n",
    "\n",
    "For $\\alpha=2$, we get a flat rotation curve. For $\\alpha=0$, we get the homogeneous sphere, which gives an increase in the circular velocity with radius.\n",
    "\n",
    "The potential of the system can be computed to be equal to\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(r) =\n",
    "\\begin{cases} \\frac{4\\pi G\\rho_0 r_0^\\alpha}{(\\alpha-2)(\\alpha-3)} r^{2-\\alpha}, \\,\\,\\,&& {\\rm if}\\ 2\\lt\\alpha\\lt 3 \\\\\n",
    "\\infty \\,\\,\\,&&{\\rm otherwise}\\,.\n",
    "\\end{cases}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Double power law models\n",
    "\n",
    "Given that single power law density profiles cannot be realistic, we will now consider double power law models. There is a whole family of such models\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho(r) = \\frac{\\rho_0}{\\left(\\frac{r}{r_{\\rm s}}\\right)^{\\gamma}\\left(1 +\\left[\\frac{r}{r_{\\rm s}}\\right]^{\\alpha}\\right)^{(\\beta-\\gamma)/\\alpha}}\n",
    "\\end{equation}\n",
    "\n",
    "We have already encountered one such model previously, $\\alpha=2, \\beta=5, \\gamma=0$, the Plummer sphere. But there are many others which obey the conditions that the logarithmic slope $\\alpha \\lt 3$ in the inner regions and $\\alpha \\gt 3$ in the outskirts.\n",
    "\n",
    "\n",
    "| $\\alpha$\t| $\\beta$  \t|  $\\gamma$ \t|   Profile\t|\n",
    "|---\t|---\t|---\t|---\t|\n",
    "|   2\t|   5\t|   0\t|  Plummer sphere \t|\n",
    "|   1\t|   3\t|   1\t|  NFW profile \t|\n",
    "|   1\t|   4\t|   2\t|  Jaffe profile \t|\n",
    "|   1\t|   4\t|   1.5\t|  Moore profile \t|\n",
    "|   2\t|   2\t|   0\t|  Cored isothermal profile \t|\n",
    "|   2\t|   3\t|   0\t|  Modified Hubble profile \t|\n",
    "\n",
    "The logarithmic slope of the density profile at small radii is $-\\gamma$ while that at large radii is $-\\beta$. The third parameter defines how sharp the transition between these two ranges is.\n",
    "\n",
    "The Navarro Frenk and White profile is a good fit to the density profiles of dark matter halos in numerical simulations (at least in the inner regions).\n",
    "\n",
    "The mass enclosed within the radius r, for the NFW profile is\n",
    "\n",
    "\\begin{equation}\n",
    "M(\\lt r) = 4\\pi r_{\\rm s}^3 \\rho_0 \\left[\\ln\\left( 1+ \\frac{r}{r_{\\rm s}} \\right) - \\frac{r}{r+r_{\\rm s}} \\right]\n",
    "\\end{equation}\n",
    "\n",
    "The potential can be found by integrating this mass profile divided by $r^2$ from $r$ to $\\infty$. This results in the expression\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi = - 4 \\pi G \\rho_0 r_{\\rm s}^3 \\frac{\\ln(1+r/r_{\\rm s})}{r} \n",
    "\\end{equation}\n",
    "\n",
    "- Code up these different potentials.\n",
    "- Work out the circular velocity versus radius for each of these profiles and compare them to each other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-spherical mass distributions\n",
    "\n",
    "Galaxies do not show spherical distributions, but do show some symmetries. Given the observed symmetries, we can consider mass distributions which have triaxiality or are flattened. Ellipsoidal distributions are more general than spherical systems, but hopefully are closer to reality.\n",
    "\n",
    "Ellipsoidal distributions can be described as a density which is a function of a dimensionless parameter $m^2$, where\n",
    "\n",
    "\\begin{equation}\n",
    "m^2 = a_1^2 \\sum_{i=1}^{3} \\left(\\frac{x_i}{a_i}\\right)^2\\,,\n",
    "\\end{equation}\n",
    "\n",
    "where we have assumed that the axis are aligned with the principal axes of the coordinate system we use.\n",
    "\n",
    "The axes $a_1$, $a_2$ and $a_3$ can be chosen to be ordered such that $a_1\\gt a_2 \\gt a_3$. The triaxiality parameter can be defined as the ratio of eccentricities in the 1-2 and 1-3 plane.\n",
    "\n",
    "\\begin{equation}\n",
    "T = \\frac{1-a_2^2/a_1^2}{1-a_3^2/a_1^2} \\,.\n",
    "\\end{equation}\n",
    "\n",
    "The special cases of spheroids are those which have two of the three axes are equal. The prolate spheroids have the smaller axes equal, while the oblate spheroids have the larger axes equal. For these spheroids, we can define the ellipticity $\\epsilon=1-q$, and $e=\\sqrt{1-q^2}$, with q the axis ratio."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generalization of Newton's theorems\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two theorems which turn out to be very useful for ellipsoidal distributions. A shell of similar **concentric** ellipsoids is called a Homoeoid. The isopotential surfaces outside of a uniform density homoeiod shell are confocal with the shell, and the potential interior to the shell is constant. This implies that a mass inside the homoeoid stays will not experience a net gravitational force.\n",
    "\n",
    "One can see that the isopotential surfaces will become spherical because these surfaces are confocal with the homoeoid shell. As we go further and further away from the system, the distance between the two focii is much much smaller than the radius at which we want to compute the potential (and thus it behaves more and more like a point mass)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The potential for the ellipsoidal distribution can be obtained by summing the potentials of homoeoids that make up the distribution and can be shown to be equal to\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi({\\bf x})  = -\\pi G \\left(\\frac{a_2a_3}{a_1}\\right)\\int_{0}^{\\infty} \\frac{\\left[ \\psi(\\infty) - \\psi(m) \\right] d\\tau}{\\sqrt{\\prod(\\tau+a_i^2)}}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\psi(m)$ is given by\n",
    "\n",
    "\\begin{equation}\n",
    "\\psi(m) = \\int_0^{m^2} \\rho(m^2) dm^2\n",
    "\\end{equation}\n",
    "\n",
    "and the relation between $m$ and $\\tau$ is given by\n",
    "\n",
    "\\begin{equation}\n",
    "m^2 = a_1^2 \\sum_{i=1}^3 \\frac{x_i^2}{\\tau + a_i^2}\\,. \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inhomogenous distribution\n",
    "\n",
    "For a very inhomogeneous distribution the potential can be computed by using a multipole expansion. For this purpose one can start with a shell of matter with inhomogenous distribution and computing the potential for such a shell. The Poisson equation that needs to be solved then is $\\nabla^2\\Phi=0$ with appropriate boundary conditions. The solution can be found by a separation of variables where $\\Phi(r, \\theta, \\phi) = R(r)P(\\theta)Q(\\phi)$. \n",
    "\n",
    "The potential can be computed as\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(r, \\theta, \\phi) = -4\\pi G \\sum_{l=0}^{\\infty} \\sum_{m=-l}^{l} \\frac{Y_{l}^{m}(\\theta,\\phi)}{2l+1} \\left[  \\frac{1}{r^{l+1}} \\int_0^r \\rho_{lm}(r')r'^{(l+2)} dr' + r^l \\int_r^{\\infty} \\frac{\\rho_{lm}(r')}{r'^{(l-1)}} dr' \\right]\n",
    "\\end{equation}\n",
    "\n",
    "where \n",
    "\n",
    "\\begin{equation}\n",
    "\\rho_{lm}(r) = \\int d \\cos\\theta\\, d\\phi\\, Y_l^m(\\theta, \\phi) \\rho(r, \\theta, \\phi)\n",
    "\\end{equation}\n",
    "\n",
    "For every $l$ there are $(2l+1)$ terms. The monopole is the $l=0$ terms, the dipole the $l=1$ terms, the quadrupole term is the $l=2$ terms and so on. On large scales the monopole term dominates in gravity, the quadrupole corresponds to the flattening. The higher order terms fall off quicker with radius in contrast to the monopole."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Azimuthally symmetric distributions\n",
    "\n",
    "#### Kuzmin disk model\n",
    "\n",
    "Given that disks are azimuthally symmetric, we next consider azimuthally symmetric distributions. Before considering general expressions for azimuthally symmetric distributions, it is worthwhile studying a simple example of a simple potential density pair. We will work in the cylindrical coordinate system for this purpose.\n",
    "\n",
    "Let us consider the Kuzmin potential\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi_{\\rm K}(R, z) = - \\frac{GM}{\\sqrt{R^2+(a+|z|)^2}}\\,.\n",
    "\\end{equation}\n",
    "\n",
    "Note that if we consider $z\\gt 0$, this potential is same as that from a point mass situated at $z=-a$, while for $z\\lt 0$, it corresponds to point mass situated at $z=a$. Given that the potential from a point mass obeys the Poisson equation with $\\nabla^2\\Phi=0$, it should be the same for the Kuzmin potential everywhere except near z=0. Near z=0, we can use the Poisson equation to obtain the surface density which gives rise to the Kuzmin potential.\n",
    "\n",
    "Consider a infinitesimally small cylindrical volume at radius $R$. The Poisson equation can be integrated over this volume. \n",
    "\n",
    "\\begin{eqnarray}\n",
    "4 \\pi G \\Sigma(R) dA &=& \\int \\nabla^2 \\Phi_{\\rm K} dV \\,,\\\\\n",
    "&=& \\int \\nabla \\Phi_{\\rm K}.{\\hat{n}} dA \\,,\\\\\n",
    "&=& 2 \\left.\\frac{\\partial \\Phi_{\\rm K}}{\\partial z}\\right|_{z=0} dA\n",
    "\\end{eqnarray}\n",
    "\n",
    "Computing this derivative is straightforward and leads us to\n",
    "\n",
    "\\begin{equation}\n",
    "\\Sigma(R) = \\frac{Ma}{2\\pi(R^2+a^2)^{3/2}}\\,.\n",
    "\\end{equation}\n",
    "\n",
    "The above infinitesimal disk density potential pair, was introduced by Kuzmin (1956).\n",
    "\n",
    "One can compute the circular velocity of the Kuzmin disk in the mid-plane, by equating the centrifugal force with the gravitational force at a given distance R from the disk.\n",
    "\n",
    "\\begin{eqnarray}\n",
    "v_{\\rm c}(R) &=& \\sqrt{R \\frac{\\partial \\Phi_{\\rm K}}{\\partial R}}\\,,\\\\\n",
    "&=& \\sqrt\\frac{GMR^2}{(R^2+a^2)^{3/2}} \\\\\n",
    "&=& \\frac{1}{(1+a^2/R^2)^{3/4}} \\sqrt\\frac{GM}{R}\n",
    "\\end{eqnarray}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$v_{\\\\rm c}$ [km$s^{-1}$]')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAF3CAYAAABHZycdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3yV9d3/8dc3G0JIyGJDCHtDCLJBJCpuxD2rVXHPKq33r+3d3rZV1LraWoXaarW1Ko4qokIEWQ4IG9mEPTLJXiS5fn9cJxBIUAJJrnPO9X4+HudxkutkfI7E8z7fbSzLQkRE3C3A6QJERMR5CgMREVEYiIiIwkBERFAYiIgIEOR0AacrNjbWSkhIcLoMERGfsXLlymzLsuLqe8xnwyAhIYG0tDSnyxAR8RnGmN0ne0zdRCIiojAQERGFgYiIoDAQEREUBiIigsJARERQGIiICAoDERFBYSAiIigMREQEhYGIiKAwEBERfHijOmm46mqLHVlFrN+fT05RBbklFRSVVdIyJJCIsCAiWwSTGNeKPu0iiGkV6nS5ItKMFAZ+rryyitSNmXy4eh/Ld+ZSUFZ59LGgAENEWBAlFVWUV1Yf932xrUIZ0S2aCb3jmNArjratw5q7dBFpRgoDP1VYdoRXF6Xz7+V7yC2uoH1kGBcNak9SlzYM6RxFu8gwWoUGYYwBoKKymsMlFWzLKGLzoQI2Hihg2Y5sPl1/EIBBnSK5IqkTlw7uQJvwECefmog0AWNZltM1nJbk5GRL5xnUVV1t8W7aXp6dt5XsonLO69eW60d0YVzPOAIDTIN+lmVZbD5UyFdbspiz7gDfHyggONCQ0rctt4xO4Kxu0UfDRES8nzFmpWVZyfU+pjDwHzlF5Tz0zhqWbMtmWNc2/OrifgzpHNVoP3/jgQLeX7WPD1bt43DJEQZ3imTa+O6c378tQYGaiyDi7RQGLrBydy73/ms1uSUV/OaS/lx3Vucme9deWlHF+6v28bcl6ezKKaFrTEseTunFJYM7NLj1ISLNR2Hg5xZuzuTOt1bSrnUYL9+QxICOkc3ye6uqLeZvPMSLX25n08ECereN4JHzenFev7bqPhLxQj8UBmrb+7h53x9i2ptp9Grbiv/eO6bZggAgMMAweUB7Pr1/LH++fihHqqq5882VXPHXr1m7N6/Z6hCRM6cw8GFfbcnknn+ton+HSP51+0jHZvkEBBguHtSBeQ+PZ8YVA9mTW8plf1nGo++tJbOwzJGaRKRhFAY+amtGIff9ezW92kbw5m1nEdki2OmSCAoM4JrhXVj46ATumtCdj9ccYOIzXzFz8Q4qq6p//AeIiGMUBj4ot7iC295YQVhwIH/7STIRYc4HQW0RYcH84oI+zHt4PCMTY/jD3M1c8udlrN5z2OnSROQkFAY+xrIsHnpnDRkF5cy6eRgdolo4XdJJJcSG87efJPPKjcM4XFzB1L9+za//u4GCsiNOlyYiJ1AY+Jg3v93N4q1Z/Oqivgzt0sbpcn6UMYbJA9ox/5Hx/GRUAm9+u5vznlvMws2ZTpcmIrUoDHzI9swifv/pJs7uHceNI7s6XU6DRIQF85tL+/PRPWOIbBHMra+v4GfvriW/RK0EEW+gMPAR1dUWj763lpYhgTx9xSCfncc/uHMUH98/hvvP6cFHa/Zz7vOLmL8xw+myRFxPYeAj3lu5lzV78/j1Jf2I9/EdREODAvnZeb35771jiA4P4Y5/pvHYe2sp1FiCiGMUBj4gv+QIMz7fwvCENkwZ0tHpchrNgI6RfHzfWO6b2IP3V+1j8gtL+DY9x+myRFxJYeADnpu/hbySCn576QCf7R46mZCgAB49vzez7x5NcKDhulnf8rs5Gyk7UuV0aSKuojDwcruyi3nruz3cMKIr/Tq0drqcJpPUpQ1zHxzHjSO68relO7nsz8vYfKjA6bJEXENh4OVe/HIbwYGG+yf1cLqUJtcyJIgnpgzg9VuHk1NcwaV/Xsbfl+6kuto3N1MU8SXNFgbGmKQTPp/muc2ode1KY0yKMWZ6c9XlzbZlFPLRmv38ZHQC8RG+PWjcEGf3jueLh8Yxvmcs/zdnI7e8voLMAu1xJNKUmiUMjDEpwKwTPk+1LGsmkOgJgCQAy7JSgbwTw8ONXkjdRsvgQO4c393pUppdTKtQZt2czO+mDGD5zhwmv7iEVE1BFWkyzRIGnhf43FqXEoEUz8fpns+vAfJqXUvBxdKzipi74SC3julGtEvPHDbGcOPIrsy5fyztWodx+z/T+N//btDgskgTcGTMwLKsmZ5WAUASkAZEcXxgxDR7YV7ktaU7CQ4M4JYxCU6X4rge8RF8eO9ofjqmG298s5spf1nGtoxCp8sS8SuODiB7uoLmW5a1ysk6vE1OUTmzV+5j6tCOxLYKdbocrxAaFMivL+nHP24ZTlZhORf/aSlvfbsbXz2pT8TbOD2bKMWyrKc9H+cB0Z6Po4A6q488A85pxpi0rKys5qqx2b317R7KK6u5fVw3p0vxOhP7xPPZQ+M4q1s0v/xoA/f8a5X2NxJpBI6FgTFmWk0QeAaU38EeO8Bzn3ri93i6l5Ity0qOi4trvmKbUUVlNW9+u4tz+sTTIz7C6XK8UnxEGG/cehaPX9CH+RszuPClJaTtyv3xbxSRk2qu2URXAsme+5oX/xnGmB3GmMMANV1Fnsfy3Np1NH9jBtlFFdw0yrd2JW1uAQGGOyd0Z/bdowkMMFwz81v+9OU2qrQmQeS0GF/tc01OTrbS0tKcLqPR3fi379iZXczi6RMJDPCvrSeaSmHZEf7fhxv4eO0BRiZG88I1Q2kX6Z51GSKnyhiz0rKs5Poec3rMQGrZlV3M0u3ZXDu8s4KgASLCgnnx2iE8c+Ug1u7N54IXF/PlJq1JEGkIhYEX+c+KvQQGGK4e3tnpUnyOMYarkjvzyf1jaRfZgtveSOO3n3xPeaXWJIicCoWBl6iqtnh/1T4m9o6nrY+fV+CkHvGt+PCe0dwyOoF/LNvF1Je/Zmd2sdNliXg9hYGX+HpHNlmF5UxN8p/zCpwSFhzIby7tz8ybhrE/r5SLX1rC+yv3OV2WiFdTGHiJj1YfICI0iHP6xDtdit84r387PntwHP07RvKz99by8DtrKCqvdLosEa+kMPACpRVVfL7hIBcObE9YcKDT5fiV9pEtePuOkTyU0pP/rtnPxS8tYf2+fKfLEvE6CgMvMH9TBsUVVVw2tIPTpfilwADDQym9+M+0UVRUVjP1r8uYtThd5ySI1KIw8AKfrD1A29ahjOzm6r35mtxZ3aKZ++A4zukTz+/nbuLW11eQVVjudFkiXkFh4LDi8koWb83iggHtCdDagiYX1TKEV24cxhNTBvBteg4XvLiYr7ZkOl2WiOMUBg77aksW5ZXVTB7QzulSXMMYw00ju/LxfWOJCQ/lln+s4Ik5G7UmQVxNYeCwzzYcJLZVCMMTon/8i6VR9W4XwX/vG8PNo7ry2tKdXP6Xr9meWeR0WSKOUBg4qOxIFQs3Z3Juv3bafsIhYcGB/N9lA/jbzckczC/lkj8t5e3le3ROgrhOkNMFuNmy7dkUV1RxQWN2EVVWQN5uKMqE4iyorHWQfFgktIyB8DiI7AyB+uevkdKvLZ8/NJ5H3l3D4x+s56stmTw1dRBtXHrkqLiPXg0c9OXmTFqFBjEy8QxmEZXlQ/pXsD0V9q+GrM1QfQqHvQQEQZtuENcbOiZBp+HQYSiEuvcMhbatw3jzpyN4belOnv5iM5NfXMwfrxrC2J6xTpcm0uQUBg6xLIuFmzMZ2yOWkKAG9tZVV0P6Alj5Omz5DKorITQSOg+HHpMgrg9EtLNbACEta34hlOVBcQ4UZUDuDsjeZofH5jn215gA6Jhs/4zu50DHYRDgrkVwAQGGO8YnMqp7DA/+ZzU3vvYdd4zrxqPn9yY0yF3/LcRdFAYO2XSwkIP5ZTx8bgO2n6iuhs2fwMI/2C/iLWNgxF3Q52L7nf3pdvuU5ML+VbDnG9ixAL56Cr56EsLjoe8l0O8y6DrGVd1KAzpGMuf+cfxh7iZmLdnJ0u05vHTtEHq2dW/LSfybDrdxyJ8XbOPZeVtZ/v8mER9xCruUHloPnzwI+1dCbC8YPx36XQpBoY1fXEmuHQqbPoat86Cy1A6GQVfDkOuhbf/G/51e7MtNGUyfvY7C8kp+MbkPt4xO0JoQ8Uk/dLiNwsAhU19eRmW1xcf3jf3hL6yqhEVPwdLnoUUbOPf/YNA1zdd9U1EM2+bD+vdg6+d2l1T7wTDsFhh4NYS2ap46HJZVWM7P31/Hgs2ZjOsZyzNXDtZpauJzFAZeJq+kgqFPzOf+c3ryyLm9Tv6FhRnw/m2wawkMvg7O/wO0dHA9QnEObJgNq96EjPUQEgGDr4Hht0N8X+fqaiaWZfH28r08MWcjIUEB/G7KAC4ZrP2kxHfo2Esv8/WOHCwLxv/QLJXMTTDzbNiXBlNegctfcTYIAMJjYMSdcNcSuC0V+lxkB8PLI+HNy+0ZTT765uJUGGO4fkQX5j44jm6x4dz/9moe+s9q8ktPYfaWiJdTGDhgybZsWoUGMbhzVP1fsOc7+PtksKrhtnkw5LrmLfDHGGPPXJr6KjyyCSb9GjI2wltXwMujYPVb9noHP9UtNpzZd43ioZSefLLuIJNfWMzSbdlOlyVyRhQGzcyyLJZsy2JU9xiCA+v5z793hf0uu2WMHQTtBzV/kQ0RHgPjfgYPrYfLX7XXL/z3XnhpCHzzsj3m4IeCAgN4KKUXH9w9mpYhgdz42nf86qMNlFTo8BzxTQqDZrY7p4R9h0sZV18X0aH18K8roFU83DoX2nRt/gJPV1AIDL7W7kK64X1okwBfPA7P94dFT9uL4/zQ4M5RfPrAOG4b2423vtvNBS8uIW1XrtNliTSYwqCZLdludyeM7XFCGBQcgLeuhJBWcPN/7UVjvsgY6Jlih9lP50HnEbDw9/D8QFj4JJTmOV1howsLDuRXF/fj7TtGUlVtcdWr3/DkZ5soO6JdUMV3KAya2bc7cmgfGUa32PBjF4+Uwn9ugPJCuP5d32oR/JAuI+D6d+DOxdBtnD1F9oWB9qI5PwyFkYkxfP7QeK4d3plXF6Vz6Z+Xsm6f/z1P8U8Kg2ZkWRbfpucwKjEGY2otWpr7KBxYDVfMgnYDnCuwqbQfDNf+C+5cAokTYNEMeHEQLH7GDkA/0io0iCenDuIftw4nv/QIl7/8Nc98sVlnJYjXUxg0o+2ZReQUVxy/Md362fbsm/GP2lM1/Vn7QXDNW3ZLocsoWPA7eHEwLHvJbh35kYm945n38ASmDu3IXxbu4JI/qZUg3k1h0Iy+Tc8BOBYGh3fDnIftfvUJv3CwsmbWfrDdfXT7l/bH838FLw2FFa9Blf/M2Y9sEcwzVw3mH7cOp6C0kstf/poZn2/WWIJ4JYVBM/o2PZcOkWF0jm5hL8765AH7/oq/uWoTuKM6JcNNH8ItcyGqC3z6CPw5Gda9a2/K5ycm9o5n3iPjuTKpE3/9agcX/2kpq/ccdroskeMoDJpJzXjByJrxgtVv2ecQnPtb+4XQzRLGwE+/sAfPQyLggzvg1XH2Jnl+sqK5dVgwM64cxBs/PYvi8kqu+OvXPDFno9YliNdQGDSTHVnF5BRXcFa3aCjOhnn/z94WetitTpfmHYyBXufb4wlXvGYvVvv3VfD6RfZCPD8xoVcc8x4ez3VndeG1pTs5X6uXxUsoDJrJyt32QqTkhGhY8IT9Ynfx8xCgf4LjBATAwCvh3uVw4bP2ATyvpdhTb7O2OF1do4gIC+b3lw/knWkjCQ4I4MbXvuOx99aSV+K/W3iI99MrUTNZseswbVoG071yG6x8A8660z5yUuoXFAJn3QEPrIaJv4T0RfaGeB/fby/Q8wMjEmOY++A47jm7Ox+s3k/Kc4v5dN1BfHUnYfFtCoNmsnL3YYZ1aYOZ9yt799EJ050uyTeEtoIJj8GDa+0AXfM2vJQEqb/1iy0uwoIDmT65Dx/fN4Z2kaHc++9VTHtzJQfy/GuqrXg/hUEzyCosZ2d2MVMit9pnE4yfDi1OsmOp1C88Bi54Cu5Ps4/iXPqcvUbhm5ehstzp6s5Y/w6RfHTPGB6/oA9LtmVx7nOLeG3pTiqr/GdWlXg3hUEzWLn7MGAxcd9fIbILJGvQ+LS1SbBXat+5GDoMtTfD+3MyrHvP56ejBgUGcOeE7sx/eALDu0XzxJyNTHl5Gev3+X4LSLyfwqAZrNpzmPOD1hCesx7O/kXTnFvsNu0H22sUbvoQwiLhg9th5gTYsdDpys5Y5+iW/OOW4fzl+iQyCsq57C9L+e0n31NUrmmo0nQUBs1g9e5cfhb2X4jqah8qL42n+zkwbTFcPtPe/O7NKfZ5EAfXOV3ZGTHGcNGg9nz5swncMKIrr3+9i5Q/LuLzDYc0wCxNQmHQxI5UVRNxYCm9KrfC2IchMNjpkvxPQIB9FvN9K+C838P+VfDqePhgGuTtcbq6M9I6LJgnpgzgg7tHE9UymLveWslPX1/B7hz/PDRInGN89V1GcnKylZaW5nQZP2rD/nxyX7mI4eEZtHh0g7qImkNpHix9Hr57xT469Kxp9mlsTp8hfYaOVFXzxte7eH7+Vo5UW9w9oTt3n92dsOBAp0sTH2GMWWlZVnJ9j6ll0MR2blzO+MD1lCfdriBoLi2i7G0+7l8JA6+Cb/4CLw6xA8KHd0cNDgzg9nGJLHj0bCb3b8eLX27j3OcX8eWmDKdLEz+gMGhibb//O6WEEjn2DqdLcZ/ITjDlZbh7GXQdBam/sdcorPonVPnuYGzb1mG8dN1Q/n3HCEKDArntjTRufyONvbklTpcmPkxh0JSKcxiSN49vIs7H+HgXhU9r29/eMvuWT6F1B3sV8ytjYPOnPr0R3ujuscx9YBz/c2Efvt6RTcpzi3gxdZu2yJbTojBoQuUrXieESg70vMHpUgQgYSzcngpX/xOqK+E/18Pfz4fd3zhd2WkLCQpg2vjufPmzCaT0a8vzqVuZ9MdFzF2vbS2kYRQGTaW6Givtdb6r7kOn3klOVyM1jIF+l8E939obBR7eDf+YDP++BjK+d7q609Y+sgV/uT6Jt+8YSURYEPf8axXXzfqWTQcLnC5NfITCoKmkLySsaA9vVaYwsGOk09XIiQKDIfmn9kZ4k/7Xbh38dQx8cKcdED5qVPcY5tw/liemDGDzoUIuemkJv/xoPYeLtSOq/LBmCwNjTNIJn19pjEkxxkz/oWs+a/VbFAe0Zn3EeGJaaRaR1wppCeMegQfXwOj74fsP4U/DYO50KMp0urrTEhQYwE0ju/LVo2dz86gE3l6+l7Of/YrXl2mvIzm5ZgkDY0wKMKvW50kAlmWlAnnGmKT6rjVHbU2iJBc2z+GzgPH07RzrdDVyKlpGw3lP2C2FIdfDir/Z01EX/M5nd0eNahnCby7tz9wHxjGgY2t+88lGLnhxCQs3Z2o8QepoljDwvMDn1rp0DZDn+TgdSDnJNd+0fjZUVfBa8RgGdlIXkU+J7AiXvmQfrtPrPFj8jL076rIXfXaNQu92Ebx12whevWkYR6qqufX1Fdz02nI2HtB4ghzj1JhBFMeHQ8xJrvmmde9Q1KYvm6yuDOqorap9UmwPuOp1mLYIOiTB/F/DS0PtFkOl7/W/G2M4v3875j08gV9f3I/1+/O56E9LmD57LRkFZU6XJ15AA8iNLWcH7E9jQ/R5APTv0NrhguSMdBgCN31gr1GI6gqf/szeMnvNv6Ha9+bzhwQF8NOx3Vj82ERuG9OND1fv5+xnvuKF1K2UVPjuQjw5c06FQR5QsworCsg5ybXjGGOmGWPSjDFpWVlZzVJog62fDRjmWmPoEBlGm/AQpyuSxpAwFn76Odww297u4qO77WM4v//QJ89RiGwZzC8v7kfqIxOY2CeOF1K3cfYzX/Huir1UVWs8wY2cCoN3gETPx4lA6kmuHceyrJmWZSVblpUcFxfXLIU2iGXBhtnQdQzLskLp10HjBX7FGOh5rt11dPU/wQTAe7fAzPGwea5PrmbuGhPOyzcMY/Zdo+gQ1YLp769j8guLmfe9tsp2m+aaTXQlkOy5x7KsVZ7rKUCeZVmr6rvWHLU1qsyNkL2Vij6XkZ5drC4if1WzcO3ur+1zFCqK4T/XwayJsG2+T4ZCckI0H94zmpdvSKKq2mLamyu58pVvWLEr98e/WfyCtrBuTAt+D0ueZe0133HZ69uYedMwzuvfzumqpKlVVcLat2Hx0/b5CZ3Ogon/A4ln28HhY45UVfNe2j5eSN1KZmE5k/rE89jk3vRppzc3vk5bWDeXjR9B1zGsO2yPE/TXymN3CAyCpJvgvpX2FhcF++0T1/5xIexc7HMtheDAAK4f0YVFj01k+uTeLN+VywUvLuGRd9ew77B2RvVXCoPGkrUVsrdCv8vYeLCAyBbBdIgMc7oqaU5BIce2uLjwWchNhzcusUMhfZHPhUKLkEDuObsHS6ZPZNq4ROasO8g5zy7it598T1ZhudPlSSNTGDSWzXPs+94XsulgIX3bR2B8sItAGkFQKJx1Bzy4Fi54Bg7vgn9eCv+4ANK/8rlQiGoZwuMX9mXRY2dz+dCO/POb3Yx/eiFPfraJXO155DcUBo1l86fQIYmqiA5sOVRI3/bqX3W94DAYMe1YS+HwbvjnZfD3ybBjgc+FQvvIFsy4chCpj0xg8oB2zFyczrgZC/jjvC3klx5xujw5QwqDxlB4CPanQZ8L2ZNbQumRKvpqsE1qBId5Wgpr7FDI3wtvXg5/S4Etn/tcKHSLDef5a4Yw76HxnN07nj8t2M64GQv405fbKCrXwjVfpTBoDNvm2fe9L2SzZ/94tQykjpruowdWw8UvQHEmvH0NvDIOvv/I5xav9WwbwV9uSOLTB8ZyVrcY/jh/K+NmLOCVRTu0mtkHKQwaw9YvoHUniO/HpoMFBBjo2baV01WJtwoKheRb4f5VMOWvUFkK7/3EXtG87l2fO5+5f4dI/vaTZP577xgGdYriqc82M/7phby6aAfFain4DIXBmaoshx0L7R0ujWHToUK6xYYTFhzodGXi7QKD7e2y710OV7wGAYHwwR323kdpf4cjvrWB3ODOUbzx07OYfdco+rZvzZOfbWbsjAX8ZeF2Css0puDtFAZnavfXcKQYetob0205VKjFOdIwAYEw8Eq4axlc85a999Gch49tnV1e6HSFDZKcEM2bt43gw3tGM7RLG575YgtjnlrAi6nbNNDsxRQGZ2rHlxAQDAnjKKmoZE9uCb3bRThdlfiigADoewncsRBu+gjietlbZz/f3z5kpzjb6QobZGiXNvz9luF8ct9YRiTG8HzqVsY+Zc8+0jGc3kdhcKa2L4AuIyG0FdsyigDo1VZhIGfAGOg+EX7yCdy+ABLG2YfsPD8APvu5veWFDxnYKZJZNycz94FxjOsVy58XbmfsjAU8+dkmMgt9qyvMnykMzkTBQcj8HnpMAmBLht2cV8tAGk2nYXDtv+Ce76D/lGPHcb5/Oxxc53R1DdKvQ2tevmEYXzw0nkl92zJzcTpjZyzk/324nt05xU6X53oKgzOR/pV93/0cALYeKiQ0KIAu0S2dq0n8U3wfuPwVeGANjLgLtnwGr46Df06xJzD40FqFXm0jeOm6oSz42dlckdSJ99L2MfHZr7j/7dV8f8A3z5v2BwqDM7FzEbSMgbYDAbtl0LNtKwIDtA2FNJGozjD5D/DwBpj0v/a26W9OsYNh/WyfmpbaLTacJ6cOZOnPJ3LH+EQWbs7kopeWcvPfl/PNjhydp9DMFAany7LslkG38fbAH7A1o1DjBdI8WrSBcY/AQ+vh0j/ZU5zfvw1eGgJf/wnKfOcddnzrMB6/oC/LfnEOj53fm40H8rlu1rdc/vLXfPH9Iap18lqzUBicruxtUHgQuk0AIL/kCBkF5fSMVxhIMwoKhaSb7TGFa9+2z2me90t4rh98/ri9SZ6PiGwRzL0Te7D05+fwxJQB5BSXc+ebK0l5fhH//m4PZUd878xpX6IwOF27ltj33cYDsD3LHjzupZXH4oSAAOhzIdz6qX0sZ5+LYPlMeGkovHMT7PnOZ8YVwoIDuWlkVxb+7GxevHYILYID+Z8P1zP6qQU8N3+rts9uIgqD07VrKbTuCNH2sc0100rVMhDHdRgCU2faXUhjHrQP2Pn7efC3SZ5xBd9Y+BUUGMBlQzoy5/6x/GfaSJK6RPHSl9sY89QCps9ey9YM31qM5+2CnC7AJ1mWHQbdJx491nBbZhFhwQF0atPC4eJEPFp3gJTfwPjHYM2/4du/2uMK834Jw2+DYbdCeKzTVf4oYwwjE2MYmRhDelYRf1+2k9kr9/Fu2j7G94rjjnHdGNsjVueHnCG1DE5H9jZ7x8mEsUcvbcssokd8KwI0k0i8TUi4vVvqfWlw/XsQ389e0fxcP/joHjiwxukKT1liXCt+N2UgX/9iEo+e14tNBwu46bXlTH5hCe+m7dW4whlQGJyO3cvs+65jjl7anlGoLiLxbgEB9oaKN30A966wB56//whmToDXzocNH/hMF1J0eAj3ndOTpT+fyLNXDcYYmD57HaOfWsDTn2/mQF6p0yX6HHUTnY4930B4/NHxgqLySg7kl9EjXoPH4iPiesFFz8I5v7S7kJa/CrNvhYgO9vbaSTdDRDunq/xRoUGBXDmsE1ckdeSbHTm8/vUuXlm0g1cXp3Nev7b8ZHQCI7pFqwvpFCgMTsfub+z9iDx/YDsy7cHj7nEKA/ExLaJg1D0w4k7YNt8OhYW/h0Uz7E3zht9ut4C9/MXUGMPoHrGM7hHLvsMlvPntbt5ZsZfPNhyiT7sIfjI6gSlDOtIiRFvLn4y6iRoqfz/k74Guo49e2pFlh4FaBuKzAgKh92S46UP70J0Rd9nbXLx+kX3ozvJZUFbgdJWnpFObljx+QV++fXwSM64YiDGGxz9Yz8gnv+QPczexN7fE6RK9ksKgofZ+Z993HnH00o6sIoICDF1jtCeR+IGY7nD+7+GRTXDZXyAoDOY+Cs/1hTmPwKENTld4SsKCA7lmeBfmPjCW9+4axdiesby2dCfjn3Pe2M0AAB4sSURBVFnI7W+sYOHmTKq0uvkodRM11N7lENwS2g08emlHZjFdYloSHKhsFT8S0hKG3mjf9q+EFa/B6rcg7TXoNByG3QL9L7dnK3kxYwzDE6IZnhDNofwy/vXdbt5evpfUTSvoGNWC687qzNXJnYlvHeZ0qY4yvroZVHJyspWWltb8v3jmRPuP/5Y5Ry+d+9wiEmLDmXVzcvPXI9KcSnJh7duw8nXI3gqhrWHQNXYwtBvgdHWn7EhVNfM3ZvDv7/awdHs2QQGGc/u15YYRXRndPcZvp4gbY1ZallXvC5VaBg1RUQKH1sHoB45eqqyqZldOMZP6tnWwMJFm0jIaRt0LI++xZ9WtfB1W/RNWzIKOyXYoDJjq9a2F4MAALhzYngsHtmdndjFvL9/De2n2gHNCTEuuO6sLVw7rREyrUKdLbTbq12iIg2uhuhI6n3X00p7cEo5UWXSP8+4/fpFGZYw9iWLqTPjZZjj/Sfus5o/vgz/2sc9w3r/KJ/ZD6hYbzv9c2JdvHp/Ei9cOIT4ijCc/28yoJxfwwNur+TbdHdtpq2XQEPtW2Pcdj7Wy0rPsE5oSNa1U3KpltD09deTdsOdbWPkPe+1C2t+h7QB7zGHg1RAe43SlPygsOJDLhnTksiEd2ZZRyL++28P7q/bx8doDJMS05KrkzlyR1Il2kf45tqAxg4Z492Z76f5Dx44bnLU4nd/P3cSaX59LVMuQ5q1HxFuV5sGG92H1m3BgNQSGQO8LYehN9p5eAb4x37+0oorPNhzk3bS9fJueS4CBCb3iuGZ4Z87p05aQIN/qXNGYQWPZl2YvNqslPbuI6PAQBYFIbS2i7M3wht9mT0Vd/Rasewc2fmTv9jvkehhyA0R3c7rSH9QiJJCpSZ2YmtSJXdnFvLdyL7NX7uOut1YREx7C5UM7cvXwzn5xqJVaBqeqMAP+2AvO/4M9gOZxzavfUFVtMfvu0T/wzSJCZTlsmWsHw/YvActe3Tz4Wug3BcJaO13hKamqtli8NYt30/aSuimDI1UWQzpHcXVyZy4Z3J6IsGCnSzwptQwaw4FV9n2HpOMup2cXc3avOAcKEvExQaH2uoT+l0P+PnuK6pq34eP7Ye506HuxHQyJ3t2NFBhgmNgnnol94skpKufD1ft5N20v//Phev5vzvdM7t+OqUmdGNMj1qfOQ1cYnKr9K8EEQPtBRy8Vlh0hq7Bcg8ciDRXZyT5nYdyjdvfr2rftMYb170FEexh0NQy+DuL7Ol3pD4ppFcrt4xK5bWw31u7L5920vcxZe4CP1hwgPiKUKUM7MjWpI33aeX+rR91Ep+qtK6DgINzz9dFL6/blcemfl/HKjcOYPMD7d3gU8WqV5bD1c7u1sH2+PY27/RA7FAZMhVbxTld4Ssorq1iwKZMPVu9n4eZMKqst+rZvzRVJHbl0SAfiI5ybjaRuojNlWfaMiF4XHHd5Z3bNtFKtMRA5Y0Gh0O8y+1aUBRtm21NUP/85fPE4JJ5tT1HtezGEeu+AbWhQIBcMbM8FA9uTW1zBnHUHeH/Vfn736Sb+MHcT43rGMTWpI+f1a+dVu6gqDE5F/j4oybHPlq1lV3YJxkCXaG1QJ9KoWsXZ6xZG3g2Zm+3uo/Xvwkd3wZwW0PsCGHgV9EiBIO+dyRcdHsLNoxK4eVQCO7KK+HDVfj5cvZ8H/7OGVqFBXDCgHZcN6cio7jGOjy8oDE7FQc+xgO2PD4Od2UV0iGxBWLD3pLuI34nvA5N+ZR/Es3e5HQrffwjffwAt2tgzkQZeBV1G2ae5eanuca149PzePHJuL5bvyuWDVfuYu/4Q763cR2yrUC4e1J5LBncgqUuUI4fxKAxOxcG19uDxCRtx7cwpISFWrQKRZmEMdBlh3yY/ZZ+3sP5de/3Cyn9A607Qf4o9vtAhyWsP5AkIMIxMjGFkYgz/d9kAFm7O5OO1B/j38j28/vUuOka14JLBHbh0cAf6to9otmBQGJyKg+sgtjcEtzju8q7sYi4e1N6hokRcLDDYPs+513lQXmSvX1g/G757Fb75M7RJ8ExjnWpvN++lwRAWfGx8obDsCPM3ZvDx2gPMWpLOK4t20D0unEsH2wPP3WKbdmzyB8PAGPNXYCVgAKvWfZ0vBZIsy7q70Sv0BofWQbfxx106XFxBfumRJv8HEpEfEdrKnoo66GooPQyb5thdSMtegqXPQ0xPOxgGTPXqqaoRYcFHVzvnFlcwd/1BPll7gBe+3MrzqVsZ2DGSSwa35+JBHegQ1eLHf2AD/VjLINWyrPdP5QcZY65ohHq8T1EWFB6EdoOOu7wzx55JlBCjMBDxGi3aQNJN9q04GzZ9DBs+gCXPwuKnIa6vHQr9L4fYnk5Xe1LR4SHcOLIrN47sysH8Uj5dd5CP1x7gD3M38/z8baz+9bmNPlb5g2FwYhAYYxIsy9p1Kl/rNw55NqVrf3wY7K4JA7UMRLxTeCwk/9S+FWbAxv/aLYaFv7dv8f2g76X2VNb4vl7bldQ+sgW3j0vk9nGJ7MwuZtPBgiaZtHLSMDDGPMbxXUIGSAauafQqvNmh9fZ9rWMu4di00s7Rjd9cE5FGFtEWRkyzb/n7YfMcOxwWzYBFT0FMD08wXGrPGvTSYOgWG95kXdM/1DJIB1adcC2vSarwZhkb7FkKLdocd3l3TjEdIlsQGqRppSI+JbIjjLjTvhVlHguGZS/C0ucgqsuxFkPHZK+ertqYfigMUi3Lyj/h2qymLMYrHdpQ79muu3M1rVTE57WKP9aVVJJrz0ra+PGxWUkRHewVz30usndYDfTeHUnP1EnDoJ4gaFTGmCuxWxqJlmXNPOFakmVZTzfl7z8lleX2od99Lqrz0O6cEu1HJOJPWkbbp7INvRHK8mHrF3aLYdWbsHwmhEVCz/Pt14Mek7x6S4zTccrrDIwx51iWtaAxfqkxJglItyxrlTEmxfM5AJZlpRpjEo0xSZZlndhN1byytoBVBW37H3c5v/QIucUVdNU2FCL+KSzy2HTVihJI/wo2f+pZz/AuBIZC4gQ7GHpdYI9J+LiGLDqLauTfPQM4F7tlkGqMmQHM9zyWDqRQd8yieWVutO9PCIM9OSUAdI1RGIj4vZCW0OdC+1ZVCXu/s0Nh8xz4ZB7wEHQabgdDn4shtofTFZ8WR1Yge1oE6caYHcDPPZejgNxaX+b86dkZG+x3ANHdj7u8O9eeVtpVawxE3CUwCBLG2Lfzfme/Ydz8qX1L/V/7FtMDek2GXufb+yX5yDiDI2FgjInCHht4FZhljDmlFoAxZhowDaBLly5NV2CNjI0Q18v+A6hlt6dloN1KRVzMGLvXoG1/mDAd8vbCls/sMxmWz7QHoEMjocc5dldSjxQId/497sk0JAwac+LtNOBJy7LyPEFQM3Ac7Xk8Csg58Zs8A80zwT7cphHrqV/WZkgYV+fynpwSYluFEB6qrZ1ExCOq87G1DOVF9jjD1s9h2zx7l1UTAJ3OslsMvSZ73UK3U341a6oVxjUDxkAq9qI2gJrPnVOaBwX7693LZE9uiVoFInJyoa3sKal9L4bqansb/K1fwNbP4Mvf2rfILnYw9DwPEsbaYxMOcmrM4GljzHRjTDoQXWtqabIxJgXIc34m0Wb7Pr5fnYf25JZwVrfoOtdFROoICICOSfZt4uNQcMBuLWz9Ala/BStm2WOTCWOgx7l2d1Jsz2ZvNTRkaukkYCh2N83jQI5lWc+e7i+ubx1BTSh4hZqZRHG9j7tcXlnFgfxSOqtlICKno3UHGHaLfTtSBruXwfZU+/bF4/YtqsuxYOg23m5pNLGGtAxyLMt61rOtNcBTxpiplmV90BSFOS5rCwSHQ2Tn4y7vP1yKZWnwWEQaQXCYvYCtxyTgSTi8yxMMX8La/0DaaxAYYs9K6pECPc+FuD5N0mpoSBgMN8bEYG9Ul2RZVr4xZmejV+QtsjbbM4lO2JdkT67WGIhIE2mTAMNvt2+V5bDnGzsctqXC/F/Zt6gucF8aBIU26q9uyADyLE9X0STLsnYZY24H0hq1Gm+StQUSJ9a5vNcTBp3bKAxEpAkFhULi2fbtvN9B/j47GA7vavQggAYOIFuW9WWtT3fir7uYluXbB9rE9arz0N7DpYQEBRAf0fj/GCIiJxXZyR5naCKnvTfrCcHgX7K22vexves8tCenhM5tWhAQ4D3zg0VEzlRDZhMlAHdx7MAbA0wChjd6VU7L3mLfx9UNg72HSzSTSET8TkO6ie7C3kiu9v5Bjb15nXfI2mKP4Ed1rfPQntwShnVtU883iYj4roaEwfwTu4aMMdMbuR7vkLPd3mzqhD2J8kuOUFhWqcFjEfE7DQmDdGPMo9jbS9cMHF8F3N3oVTkte2u9K4/3HvbMJNK5xyLiZxoSBj/HDoHaI6fdT/K1vquyAnJ3Qr8pdR6qmVbaSS0DEfEzDQmDVy3LWl37gjFmdiPX47zDO+3TzWJ71nnoWMtAYSAi/qUhU0stY8wQY0zrWteuaOyCHJez3b6PqRsG+w6XEhEWRGQL3zisQkTkVDWkZbAAWAEYY4yF3V3UDTjtzeq80tEwqNsDtje3RIPHIuKXGhIGV9Uzm2hoI9fjvJztEB4HLerOmt17uJTucTrqUkT8T4O6ieq5NqmxCvEa2dvrnHkMYFkW+w6XaPBYRPxSQ8JgZs14gTGmtTHmFWBG05TloNwd9hqDE+QUV1B2pJpObTStVET8T0PC4CrgXGPMY8CXwHscO7PYP5QXQlFGveMF+w+XAppWKiL+6ZTDwLKs1Z5zkA9jb2P9Jf7WTZSzw76vJwz2HQ0DtQxExP/84ACyMSYX2MHx+xEZ4CpjjAGGAf5z0lluun1fz5jBPs8ag44KAxHxQz82m6jODKLaPIfd+I+jYdCtzkP7DpcS2SKY1mFaYyAi/ufHuol+cOpo7aDwnHzm23J3QkR7CKk7fXR/Xikdo9QqEBH/9GMtg1jPi/zJTnKpfbaB7+9TlLsD2tRtFYDdTdQ1RmsMRMQ//WAYWJb1i+YqxCvk7oQeKXUuW5bF/sOljO4e60BRIiJN77SPvfQ7FcVQdAiiE+o8lF96hOKKKs0kEhG/pTCocXiXfV9PN5GmlYqIv1MY1KgJg3pmEu3Ps8OgY5QWnImIf1IY1PiBlkHN6mOtMRARf3XaYeA52+Acz8cJjVWQYw7vgtBIaFH3sPv9eaW0CA6kTUutMRAR/3QmLYNk7POQARIboRZnHd4FbbqCqTuL9kBeKR2iwjD1PCYi4g/OJAzaAFG1PvZtNWFQj/15pXTUBnUi4sdOOwwsy3oGexfTd7D3L/JdlgV5eyCq/jA4kFdKx6iwZi5KRKT5nPJJZ8aYqZ4PUy3LKoCjgeD7ijKgsgzaJNR5qOxIFdlFFdqKQkT8WkNaBjHAaiDGGDPJGDOkiWpqfod32/f1hMEBz7TSDgoDEfFjp9wysCxrVq1Pd4LdWrAsy/e3sM7bY9/X0010IK8MUBiIiH9rSDfRY9gziFZgdxWt4eQb2PmWvF32fVTnOg8dOLrgTGEgIv6rId1Eq4Bp2K2Ca40xX3Bs11LflrcHwuMhuO4L/v68UoyBtq01gCwi/uuUWwZAGjDMc/Tl+01UjzPy9tTbKgA4mF9KfEQoIUFarC0i/qshZyDnW5a1oCmLcUzeXojqUu9DB/LKNF4gIn5Pb3erqyF/L0TW3zKwVx8rDETEvykMijOhqqLeloFlWRzIL6VDpMYLRMS/KQzy9tr39bQMDpccoexINe0j1TIQEf+mMMj3hMEPTCvtoK0oRMTPKQxqwiCyU52HtPpYRNxCYZC/zz7HICyyzkMH8+3Vx+omEhF/pzDI31dvqwDgQH4pIYEBxISHNHNRIiLNS2GQvxciO9b70IG8MtpFhhEQ4B+7boiInExDViA3KmNMEp4T0izLmu25diWQByRZlvV0sxSSvx86Jtf70KH8UtprWqmIuICTLYM7PSGQaIxJ8oQDlmWlAnk1nzepihIozT15N5FWH4uISzgSBp4WwA4Ay7KetixrFXANdqsA7LOVU5q8kIID9n09YVBVbZFRYHcTiYj4O6daBsOxD8lJMsZM91yLAnJrfU3Mid9kjJlmjEkzxqRlZWWdeRUF++z71nXHDHKKyqmstrT6WERcwcluohxPi6CmpfCjLMuaaVlWsmVZyXFxcWdewdGWQd0wOOCZVtpO00pFxAWcGkDewbFWQDp2SyEPiPZciwJymryK/P32fUSHOg8dyrcXnGkAWUTcwKmWQSqemUSe+xXAOydcS23yKgr2Q8sYCK77gq/jLkXETRwJA8uy0rFnDF3p+Xx2rS6jFCCv5vMmVXAAWtdtFQAcKigjJCiANi2Dm7wMERGnObbOwLKsmadyrUkVHDjpgrOD+WW0jwzDGC04ExH/5+4VyIUHIKJ9vQ8dzNOCMxFxD/eGwZEyKMk5aTeR3TLQeIGIuIN7w6DokH1fT8ugWgvORMRl3BsGBQft+9Z1wyC72F5wpm4iEXEL94ZBoWfBWb1rDDwLzlorDETEHVwcBp5uonpaBod0qI2IuIyLw+AgBIVBWFSdhw4V1GxFoZaBiLiDi8PgELRqC/WsIziYX0ZwoNEJZyLiGu4Og5NMK83ILyM+QieciYh7uDgMDkJEu3ofqll9LCLiFi4Og0PQqv4wyCgoo63CQERcxJ1hUF4EFUUQ0bbOQ5ZlcTC/TNNKRcRV3BkGRRn2fT0tg4KySkqPVKmbSERcxd1hUE/LoGaNQVu1DETERdwZBjULzuppGWiNgYi4kTvD4Gg3Ud2WQYa2ohARF3JvGAQEQYs2dR6qaRnEtw5t7qpERBzj0jDIslsFAXWf/qGCMqLDQwgNCnSgMBERZ7g0DDIgPK7ehzLyyzR4LCKu494wqGe8AOyWQTt1EYmIy7gzDIqzoNVJWgY64UxEXMh9YVBdbYdBeHydhyoqq8kuqiA+QmEgIu7ivjAoPQzVlfV2E2UVlQNaYyAi7uO+MCjOsu/r6SbScZci4lYuDINM+76e2UQZWmMgIi7lvjAoqgmDumMGNWGgloGIuI37wqA4276vp2VwqMA+7jJax12KiMu4MAyywATWuxVFZkE58RFhmHrORRYR8WfuDIPw2Pq3osjXGgMRcScXhkE2tIyt96GMwjLaavBYRFzIhWHgaRnUI7OgXPsSiYgruS8MSrLrDYOi8kqKyisVBiLiSu4Lg+KceruJMgtqjrtUN5GIuE+Q0wU0u7uXQWDdqaMZBfZWFNqXSETcyH1hENW53suZhTUtA4WBiLiP+7qJTiJD3UQi4mIKA4+MgnJahgTSKtR9jSUREYWBR0aBfdylVh+LiBspDDwyC8qJi1AXkYi4k8LAI6OwTLuViohrKQwAy7LIKCgjXi0DEXEphQFQWF5J2ZFqTSsVEddSGHBs9bFOOBMRt1IYYA8eg1Yfi4h7KQywB49BC85ExL0cDwNjzPRaH19pjEmpfa05HN2XSGMGIuJSjoaBMSYFGO75OAnAsqxUIK/m8+aQWVBOuFYfi4iLOd4yqOUaIM/zcTqQ0ly/OKOwTK0CEXE1x8LAGJPkaQXUiAJya30e01y1ZBWUa42BiLiaky2D6IZ+gzFmmjEmzRiTlpWV1WiFqGUgIm7nSBjU0yoAu4uoJiCigJwTv8+yrJmWZSVblpUcFxfXKLVYlkWmWgYi4nJOjZgmGmMSa32cBLwDJNdcA04MiyZRVF5J6ZEqhYGIuJojLQPLsmZbljUbuyUQ5bm2Co7OMMqr+byp1Uwr1VYUIuJmjs6ltCxrJjDzhM+bVc1xl2oZiIibedPUUkdkFdYsOFMYiIh7uT4MavYlitO+RCLiYgqDwjJCgwJoHabVxyLiXgqDwnLiW4fq7GMRcTWFQUG5tq4WEddTGBTquEsREYVBoVYfi4i4OgzKjlRRWFapfYlExPVcHQY1awzi1DIQEZdzdRjUrD5WGIiI27k7DGqOu1QYiIjLuToMsopqwkBjBiLibq4Og8yCcgIDDDHhIU6XIiLiKHeHQWEZMeEhBARo9bGIuJvLw6Bcu5WKiODyMMgq1FYUIiLg8jDQ6mMREZtrw6Cq2iKnqFxrDEREcHEY5BSXU21pjYGICLg4DLQVhYjIMQoDhYGIiHvDILNQq49FRGq4NgzUMhAROcbVYRARGkRYcKDTpYiIOM7VYRCn1cciIoDbw6CVwkBEBFwcBpmFZTruUkTEw7VhkFVYTmwrbV0tIgIuDYPi8kqKK6o0rVRExMOVYZBdpGmlIiK1uTIMtMZAROR4rg4DbVInImJzZRhkqmUgInIcV4ZBVmE5gQGGNi01m0hEBFwcBjHhIQQGGKdLERHxCq4Mg2ydcCYichxXhkFWUTmx2opCROQod4ZBYblmEomI1OK6MKiutsguKidWYSAicpTrwiC/9AhHqiy1DEREanFdGGRpKwoRkTrcFwaeBWcaQBYROca1YaBuIhGRY1wbBhpAFhE5xn1hUFROWHAAEaFBTpciIuI13BcGhfaCM2O0FYWISA3XhoGIiBzjWF+JMWaa58PulmX93HPtSiAPSLIs6+mm+L1ZheV0jWnZFD9aRMRnOdIyMMakAKmWZc0EEo0xKcaYJADLslKBvJrPG5tWH4uI1OVUN1EikOL5ON3z+TXYrYKaayn1fN8Z+/kFfbhscIem+NEiIj7LkW4iT4ugRhLwDjAMyK11PebE7/N0LU0D6NKly2n97quTO5/W94mI+DNHB5A9XUHzLctadSpfb1nWTMuyki3LSo6Li2vi6kRE3MPp2UQptQaK84Boz8dRQI4zJYmIuI9jYWCMmVYTBJ4B5Xewxw7w3Kc6VZuIiNs4OZtohjFmhzHmMEBNV5HnsbxT7ToSEZEz59QAcirQpp7rM+v5chERaWJOjxmIiIgXUBiIiIjCQEREFAYiIoLCQEREUBiIiAgKAxERAYxlWU7XcFqMMVnA7lP40lggu4nLaS7+9FzAv56PPz0X0PPxZmfyXLpallXvxm4+GwanyhiTZllWstN1NAZ/ei7gX8/Hn54L6Pl4s6Z6LuomEhERhYGIiCgMREQEhYGIiKAwEBER3BEG/rQttj89F/Cv5+NPzwX0fLxZkzwXv59aKiIiP84NLQMREfkRCgMfYoxJcroGcQ9jzHSna5Dm47dhYIy50hiT4i9/0J6zoWc5XUdjMMZM89xmOF1LY/D8naX4y/OBo39vw52u40zV/JsYY6Y5XUtjMMYkeV7brmzsn+2XYVDzDtpz1nKeP7yj9jyXXKfrOFOeF5lUz3nXiZ7PfZbnb+tcz79Pkj/8rfmZacaYHUC604U0kjsty5qN/f9Oo/6t+WUYANcAeZ6P0wGffsHxM4kc+/dI93zusyzLWmVZ1s89nyZalrXK0YIagTEmyRNu/uAqy7K6+8Pz8bQGdgBYlvV0Y/+tBTXmD/MiURz/LjrGqULkeJ4WQY0k4B2namlMnu7IO52uo5FEO11AI0oyxgAkWZb1tNPFnKHhcLQ1mtLYz8dfWwbi5Tx/0PP94Z002O/UgDuNMVFO13Im/KxVUPMOOhWI8fUuSY+cmv9nGnvcwF/DII9j726igBwHa5H6Nfo7Gyd4BvRq+m7TAV8fqEysNUDZ6P3SzemEgdYcfLxLEruLqGbsI51GHuD31zB4h2P/8ImAz7/T8fxRJzfFLILmZoyZVhMEfvBuLYXj33j49EClZVmzPQOU0djPx5elc+z//e5AmoO1NIZUjn9dW9GYP9xvVyB7ppKlYw/q+dNSdJ/mefF/D3tMJxp7gM9nw9rTLXQ19vM517Isfxk38Aue14Fc7NcBf2iJNtnz8dswEBGRU+ev3UQiItIACgMREVEYiIiIwkBERFAYiIgI/rsdhUijM8YkAjOw53i/6rncHXjVsqx61xd41oWkWpaVV9/jIt5CYSByiizLSjfGvIO9nmAmHA2IlUCbk3xbooJAfIG6iUTOgKdFUO9KXc+CNAWB+AS1DEROk+fFfhpw1Um+5Grg3VpfX7P1xrnAk2oxiDdRy0Ck4aI9L+xXA909e/nUJ+qEF/w7PVtvvAM83tRFijSEwkCk4XIty6o5rS2vviMVPWMJxw0qW5ZV04Lwh03gxM8oDETOzA7sbp8TpZzYYjDGzPDMLvL540vF/ygMRBrmxFPAcvFsK+xpDdQ47p2/p/WQUzsgfPmsAPE/CgORU+R58R7GsTEDPC/uuZ53/FG1vu7EE9xSge4nnN/gT8dLio/TFtYijcwYM90f9s4Xd1HLQERE1DIQaUw14wA1h5aL+AqFgYiIqJtIREQUBiIigsJARERQGIiICAoDEREB/j/tcGp7Ot6bkwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "class Kuzmindisk:\n",
    "    def __init__(self, mass, a):\n",
    "        self.mass = mass\n",
    "        self.a = a\n",
    "        self.gee = 4.2994e-9\n",
    "        \n",
    "    def vcirc(self, R):\n",
    "        return np.sqrt(self.gee*self.mass*R**2/(R**2+self.a**2)**(3./2.))\n",
    "    \n",
    "    def mass_within_R(self, R):\n",
    "        return -self.mass*self.a*(1/(R**2+self.a**2)**(1./2.)-1/self.a)\n",
    "    \n",
    "    def vcirc_mass_within_R(self, R):\n",
    "        return sqrt(self.gee*self.mass_within_R(R)/R)\n",
    "\n",
    "# 10^11 Msun disk with radius 10 kpc\n",
    "aa = Kuzmindisk(1e11, 0.010)\n",
    "RR = np.logspace(-3, np.log10(0.060), 100)\n",
    "\n",
    "ax = pl.subplot(111)\n",
    "plot(RR/0.010, aa.vcirc(RR))\n",
    "plot(RR/0.010, aa.vcirc_mass_within_R(RR))\n",
    "ax.set_xlabel(\"R/a\")\n",
    "ax.set_ylabel(r\"$v_{\\rm c}$ [km$s^{-1}$]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Miyamoto-Nagai disk model\n",
    "\n",
    "The density potential pair of Kuzmin can be generalized to thicker disks by considering the potential (Miyamoto and Nagai 1975)\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi_{\\rm M}(R, z) = -\\frac{GM}{\\sqrt{R^2+(a+\\sqrt{z^2+b^2})^2}} \\,.\n",
    "\\end{equation}\n",
    "\n",
    "Depending upon the magnitude of $a$ compared to $b$, one can get a spherical or a highly flattened system. One can obtain the density distribution corresponding to this potential via the Poisson equation and it turns out to be\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho_{\\rm M}(R, z) = \\left(\\frac{b^2M}{4\\pi}\\right) \\frac{aR^2+\\left(a+3\\sqrt{z^2+b^2}\\right)\\left(a+\\sqrt{z^2+b^2}\\right)^2}{\\left[ R^2+\\left(a+\\sqrt{z^2+b^2}\\right)^2\\right]^{5/2}\\left(z^2+b^2\\right)^{3/2}}\n",
    "\\end{equation}\n",
    "\n",
    "### General azimuthally symmetric distributions\n",
    "\n",
    "There are different ways to approach the problem of general azimuthally symmetric distributions. One approach for obtaining potential density pairs for azimuthally symmetric distributions is by considering oblate distributions, and taking the limit where the non-equal axis becomes very thin. \n",
    "\n",
    "Another approach is by directly performing the integral over the density distribution to get to the potential. This involves elliptical integrals.\n",
    "\n",
    "One can also approach the problem by obtaining solutions to the Poisson equation for a thin disk away from the disk, this approach involves Bessel functions and Hankel transforms. Let us take a brief look at this approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bessel function approach\n",
    "\n",
    "For an axially symmetric system, the Poisson equation far away from the disk reduces to\n",
    "\n",
    "\\begin{eqnarray}\n",
    "0 = \\nabla^2\\phi &=& \\frac{1}{R}\\frac{\\partial }{\\partial R} \\left(R\\frac{\\partial \\Phi}{\\partial R}\\right) + \\frac{1}{R^2}\\frac{\\partial^2}{\\partial \\phi^2} \\Phi + \\frac{\\partial^2 \\Phi}{\\partial z^2}\\\\\n",
    "&=& \\frac{1}{R}\\frac{\\partial }{\\partial R} \\left(R\\frac{\\partial \\Phi}{\\partial R}\\right) + \\frac{\\partial^2 \\Phi}{\\partial z^2}\n",
    "\\end{eqnarray}\n",
    "\n",
    "where the second term in the first equation above for the Laplacian is zero due to axial symmetry. Given that the two terms in the following equation are independent, we can solve this equation using separation of variables. Defining $\\Phi(R, z) = P(R)Q(z)$, we obtain\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1}{P(R)R}\\frac{d}{dR}\\left( R\\frac{dP}{dR}\\right) = -\\frac{1}{Q(z)} \\frac{d^2Q}{dz^2}\n",
    "\\end{equation}\n",
    "\n",
    "Given that the LHS is independent of $z$ while the RHS is independent of $R$, the above expressions can only be equal if they were equal to a constant. We can set both equal to -$k^2$. For the z component, we have\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1}{Q(z)} \\frac{d^2Q}{dz^2} = k^2 \\,,\n",
    "\\end{equation}\n",
    "\n",
    "which has the solution\n",
    "\n",
    "\\begin{equation}\n",
    "Q(z) = Q_0 \\exp(\\pm kz)\n",
    "\\end{equation}\n",
    "\n",
    "If $k$ is complex, this can result in oscillatory behaviour as $z\\rightarrow \\infty$. Thus physical solutions which would require $Q\\rightarrow 0$ as $z\\rightarrow \\infty$, demand $k \\gt 0$, with a minus sign in the above expression for $z\\gt 0$. For $z\\lt 0$, we would require the positive sign with $k\\gt 0$.\n",
    "\n",
    "The other equation,\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1}{R}\\frac{d}{dR}\\left( R\\frac{dP}{dR}\\right) = -k^2 P(R)\n",
    "\\end{equation}\n",
    "\n",
    "The solution to this equation is the Bessel's function of the first kind $J_{\\rm \\alpha}(kR)$, with $\\alpha=0$.\n",
    "\n",
    "Thus,\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(R, z) = \\exp\\left(-k|z|\\right) J_0(kR)\n",
    "\\end{equation}\n",
    "\n",
    "for positive $k$ is a solution to the Poisson equation $\\nabla^2\\Phi=0$. Now these solutions solve the Poisson equation for all $z\\neq 0$. The discontinuity at $z=0$ helps to obtain the surface density of the disk, by performing a volume integral over a small cylinder straddling $z=0$. \n",
    "\n",
    "\\begin{equation}\n",
    "\\Sigma_k(R) = -\\frac{k}{2\\pi G} J_0(kR)\n",
    "\\end{equation}\n",
    "\n",
    "Given that solutions to the Poisson equation can be added up to obtain new solutions, it is clear that a givne surface density $\\Sigma(R)$ could be decomposed in to functions made of the above basis.\n",
    "\n",
    "\\begin{equation}\n",
    "\\Sigma(R) = -\\int \\tilde\\Sigma(k) \\frac{k}{2\\pi G} J_0(kR) dk\\,.\n",
    "\\end{equation}\n",
    "\n",
    "and then use an integral over k of the Bessel functions to obtain the potential $\\Phi(R, z)$,\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(R, z) = \\int \\tilde\\Sigma(k) J_0(kR)\\exp(-k|z|) dk\n",
    "\\end{equation}\n",
    "\n",
    "Inverting the equation for $\\Sigma(R)$ gives\n",
    "\n",
    "\\begin{equation}\n",
    "\\Sigma(k) = -2\\pi G \\int \\Sigma(R) R J_0(kR) dR\\,.\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "Substituting in the equation for the potential we finally obtain\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi(R, z) = -2 \\pi G \\int \\int \\Sigma(R') R' J_0(kR') dR'  J_0(kR)\\exp(-k|z|) dk\\,.\n",
    "\\end{equation}\n",
    "\n",
    "The circular velocity in the mid-plane can be obtained by using\n",
    "\n",
    "\\begin{equation}\n",
    "v_{\\rm c}^2(R) = -R \\left.\\frac{d\\Phi(R,z)}{dR}\\right|_{z=0}\\,.\n",
    "\\end{equation}\n",
    "\n",
    "The derivative  $dJ_0(x)/dx$ equals $-J_1(x)$, which implies\n",
    "\n",
    "\\begin{equation}\n",
    "v_{\\rm c}^2(R) = -R \\int_0^\\infty \\tilde\\Sigma(k) J_1(kR) k dk\\,.\n",
    "\\end{equation}\n",
    "\n",
    "This equation can in principle be inverted in order to obtain $\\tilde \\Sigma(k)$ and thus $\\Sigma(R)$. This route involves usage of the derivatives of $v_{\\rm c}^2(R)$ with respect to $R$, which could potentially be obtained observationally. However the propagation of errors in the measurement of $v_{\\rm c}^2(R)$ to the derivatives of $v_c^2(R)$ is quite cumbersome, and hence it is better to just use a model for $\\Sigma(R)$, compute $\\Sigma(k)$ and fit to $v_{c}^2(R)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Vertical structure of the potential of infinitesimally thin disks\n",
    "\n",
    "The Poisson equation in cylindrical coordinate system is \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2\\Phi}{\\partial z^2} + \\frac{1}{R}\\frac{\\partial}{\\partial R}\\left[ R \\frac{\\partial \\Phi}{\\partial R}\\right] = 4\\pi G\\rho(R, z)\\,.\n",
    "\\end{equation}\n",
    "\n",
    "The second term on the left hand size is finite irrespective of the thickness of the disk, while the first term can grow up as we consider infinitesimally thin systems. This equation helps to approximately compute the potential of disks given their density distribution. First approximate the disk as thin, and compute the radial variation of the potential. Then we can use the equation to compute the variation of the potential in the z direction.\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2\\Phi}{\\partial z^2}  = 4\\pi G\\rho(R, z)\\,.\n",
    "\\end{equation}"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}