DMHalo_v2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La densidad para el halo de materia oscura está dado por \n",
    "$$\\begin{equation}\n",
    "\\frac{\\rho_{h} (R,z)}{\\rho_{h,0}} =  (\\frac{\\sqrt{R^{2} + z^{2}/q_{h}^{2}}}{a_{b}})^{-\\alpha_b} (1 + \\frac{\\sqrt{R^{2} + z^{2}/q_{h}^{2}}}{a_{h}})^{\\alpha_{h} - \\beta_{h}}\n",
    "\\end{equation}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def densidad_h(R,z,q_h=0.8,alpha_h=-2.,a_h=3.83,beta_h=2.96):\n",
    "    m = np.sqrt(R**2. + z**2./q_h**2.) \n",
    "    t1 = (m/a_h)**(-alpha_h)\n",
    "    d = (1. + m/a_h)**(alpha_h-beta_h)\n",
    "    rho_h = t1*d\n",
    "    return rho_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = np.linspace(-15.,15.,100)\n",
    "z = np.linspace(-15.,15.,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEMCAYAAADqG+D0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd8leX5+PHPlb0HCYQMSMLee4mIAweIirOV2opVq9b67fy2tUtbf9bWtt9aW20ddeDeAxXBLajsJRsCJGQBScgge5z798dzDkQMkJzkec45yfV+vfJKcvKcXPfJONe513WLMQallFKqo4J83QCllFKBSROIUkopr2gCUUop5RVNIEoppbyiCUQppZRXNIEopZTyiuMJRERmi8hOEckRkdvb+Hq4iLzo/voqEcly354lInUistH99pDTbVdKKXVMiJPBRCQYeBA4DygA1ojIImPMtlaX3QCUG2MGicjVwL3AN91f22OMGedkm5VSSrXN6R7IFCDHGLPXGNMIvADMO+6aecBC98evALNERBxso1JKqXZwtAcCpAP5rT4vAKae6BpjTLOIVAJJ7q9li8gGoAr4rTFm+fEBROQm4CaA6OjoicOGDevaR6D8S3MDHNoGiVkQmejr1nRvh/dASzP0HurrliibrVu3rtQY0/tU1zmdQDqjGOhvjCkTkYnAGyIy0hhT1foiY8wjwCMAkyZNMmvXrvVBU5Vj9n4KT10CC56A7Jm+bk33tuh/YNdS+F/9n+ruRCSvPdc5PYRVCPRr9XmG+7Y2rxGRECAeKDPGNBhjygCMMeuAPcAQ21us/NuRYut9XLpv29ETxKZB9SFoafJ1S5SfcDqBrAEGi0i2iIQBVwOLjrtmEbDA/fGVwEfGGCMivd2T8IjIAGAwsNehdit/VeV+/RGb6tt29ARxaYCBIwd83RLlJxwdwnLPadwGLAWCgceNMVtF5C5grTFmEfAY8LSI5ACHsZIMwEzgLhFpAlzALcaYw062X/mhqmKISICwKF+3pPuLS7PeHymGhH4nv1b1CI7PgRhjFgOLj7vtjlYf1wNXtXG/V4FXbW+gCixVRcee2JS9PD/nquNHnf1LU1MTBQUF1NfX+7opfi8iIoKMjAxCQ0O9un8gTaIr9XVHNIE4xjNMWFXs23acQkFBAbGxsWRlZaE7AE7MGENZWRkFBQVkZ2d79T20lIkKbFVFOv/hlMhECIn0+x5IfX09SUlJmjxOQURISkrqVE9NE4gKXC1N1qogXYHlDBGISz228s2PafJon87+nDSBqMB15ABgrCc15Yy4dKvXpxSaQFQg0z0gzotN1QTisKysLEpLS726pj337QxNICpw6R4Q58WlWYnbGF+3RPkBTSAqcHlWA+kqLOfEpUFLI9SW+bolfis3N5dhw4Zx3XXXMWTIEK655ho++OADTj/9dAYPHszq1asBOHz4MJdeeiljxoxh2rRpfPnllwCUlZVx/vnnM3LkSG688UZMq2T9zDPPMGXKFMaNG8fNN99MS0uLTx6jhy7jVYGrqhBCIrSIopNa7wWJTvZtW9rhD29tZVtR1akv7IARaXHcefHIk16Tk5PDyy+/zOOPP87kyZN57rnn+Oyzz1i0aBH33HMPb7zxBnfeeSfjx4/njTfe4KOPPuLaa69l48aN/OEPf2DGjBnccccdvPPOOzz22GMAbN++nRdffJHPP/+c0NBQbr31Vp599lmuvfbaLn18HaEJRAWuI8XW8JWuuHFOrCeBFEPqWN+2xY9lZ2czevRoAEaOHMmsWbMQEUaPHk1ubi4An332Ga++au2NPueccygrK6Oqqoply5bx2muvATB37lwSE60XSB9++CHr1q1j8uTJANTV1dGnTx+HH9lXaQJRgauqSCfQnRYgu9E9TtVTsEt4ePjRj4OCgo5+HhQURHNzs1ff0xjDggUL+NOf/tQlbewKOgeiAldVkS7hdVpMH5DggNgL4u/OOOMMnn32WQA++eQTkpOTiYuLY+bMmTz33HMAvPvuu5SXlwMwa9YsXnnlFQ4dOgRYcyh5ee2qum4b7YGowGSM9SSmE+jOCgqG2L66lLcL/P73v+f6669nzJgxREVFsXChdRDrnXfeyfz58xk5ciTTp0+nf//+AIwYMYK7776b888/H5fLRWhoKA8++CCZmZk+ewxiuvFyPD1QqhurKYW/DoTZ98K0W3zdmp7l0VkQHgvXvuHrlrRp+/btDB8+3NfNCBht/bxEZJ0xZtKp7qtDWCowecbgtQfivLg07YEoQBOIClS6B8R3PJsJVY+nCUQFJu2B+E5cGjRUQcMRX7dE+ZgmEBWYjhSDBEG0b9fB90it94KoHk0TiApMlQXWJsJgXUjouHj33pvK/b5th/I5TSAqMJXnQoLvli/2aJ6fe7lv9yAo39MEogJTeR4kZvm6FT1TbCoEh0GFJhAndFU595iYmC5vmyYQFXia6q2z0BO1B+ITQUGQ0N/qBaoeTROICjyV+dZ77YH4TkKmDmGdgJZzV8qfeZ64dA7EdxKzoGi9r1txau/eDgc2d+337Dsa5vz5pJdoOXel/FX5Puu9DmH5TmIm1JVDfSVExPu6NX5Hy7kr5a8q8iA4HGL6+rolPZdn+LA8D1LH+LQpJ3WKnoJdtJy7Uv6qPM+axA3SP1+fObqUN9enzQhkWs5dKV8oz9UJdF/zDB/qUl6vaTl3P6fl3LupP/eH0VfB3P/zdUt6tj/3h9HfgLl/83VLvkLLuXeMlnNXPYdn4lZ7IL6XkKlDWD2cJhAVWHQJr/9IzNQhrB5OE4gKLJ4nLF3C63uJWVCxH1wuX7fka7rz0HxX6uzPSROICiyeIRMdwvK9hExorofqg75uyVdERERQVlamSeQUjDGUlZURERHh9ffQVVgqsJTnQUSCbl7zB4nZ1vuKPIhL9W1bWsnIyKCgoICSkhJfN8XvRUREkJGR4fX9NYGowKJLeP1HYqu9IP2n+bQprYWGhpKdne3rZvQIOoSlAktFns5/+Iv4ftZ7LarYYzmeQERktojsFJEcEbm9ja+Hi8iL7q+vEpGs477eX0SqReR/nWqz8hMulzVpqyuw/ENohHW8ra7E6rEcTSAiEgw8CMwBRgDzRWTEcZfdAJQbYwYB9wH3Hvf1vwPv2t1W5YeOFENLow5h+ZNE3QvSkzndA5kC5Bhj9hpjGoEXgHnHXTMPWOj++BVglogIgIhcCuwDtjrUXuVPdAmv/9FzQXo0pxNIOpDf6vMC921tXmOMaQYqgSQRiQF+CfzhZAFE5CYRWSsia3UVRjdzdAmvTpD6jcQsqCqE5gZft0T5QCBNov8euM8YU32yi4wxjxhjJhljJvXu3duZlilnlOcBAvHeLztUXSwxEzBQWeDrligfcHoZbyHQr9XnGe7b2rqmQERCgHigDJgKXCkifwESAJeI1BtjHrC/2covlOdCXBqEhJ/yUuWQo+eC7IOkgT5tinKe0z2QNcBgEckWkTDgamDRcdcsAha4P74S+MhYzjDGZBljsoB/APdo8uhhynIgaZCvW6Fa8/w+yvb4th3KJxxNIO45jduApcB24CVjzFYRuUtELnFf9hjWnEcO8FPga0t9VQ9kDJTuhuQhvm6Jai26t1UVoGSnr1uifMDxnejGmMXA4uNuu6PVx/XAVaf4Hr+3pXHKf1UfhIZKTSD+RsT6nZTu8nVLlA8E0iS66sk8T1DJg33bDvV1yUOs3qHqcTSBqMBwNIFoD8TvJA+G6gPWQV+qR9EEogJD6W4Ii7FWYSn/4knqpTm+bYdynCYQFRhKd1mvdK2iBMqfJA+13pfqRHpPowlEBQZdgeW/EjMhKFQn0nsgTSDK/zXWQGW+TqD7q+BQ6DVAJ9J7IE0gyv95npi0B+K/kgdrD6QH0gSi/J8mEP+XPAQO74WWJl+3RDlIE4jyf6W7QIKsYRLln5KHgKtZzwbpYTSBKP9Xussq2qdFFP2Xp3eoJU16FE0gyv+V7j62VFT5J88CB50H6VE0gSj/5mqxqvDqCiz/FhEHsam6EquH0QSi/FvFfmhp0An0QKArsXocTSDKv/nhCqyiijoO1zT6uhnkHDpCfVOLr5txjKeoojG+bolyiCYQ5d885TH8ZAhrx4EqLrhvGeff9ylrcw/7pA0ul+EfH+zi3L8v4/on19DY7PJJO74meYhVcr/6oK9bohyiCUT5t9JdEJUMUb183RKKK+u47vE1RIUHExMewvxHV/Limv2OtqGmoZkfPLeef3ywm6nZvfhiTxm/eGUTxh9e9etEeo+jCUT5t5Kd0Nv3K7Cq6pu47vE1VDc08+R3p/DmD2YwbUASv3x1Mw9+7EwV2sZmF998ZAVLtx7gt3OH88JN0/j5BUN5Y2MRf13qB8tnew+z3h/a4dt2KMdoAlH+y+WCA1sgZZSvW8IPn9/AnpJqHvr2RIanxhEfFcoT103mjMHJPL0iz5EewIb95WwprOJPl4/mxjMGICLcetZAvjW1P//+ZA+vbyiwvQ0nFZsKkYlwcLNv26EcowlE+a/yfdBUA31H+7QZ6/IO88nOEn45exgzBicfvT0kOIgLRvblQFU9eWW1trdj5d7DiMDskalHbxMR7rpkJGMz4rnv/d20uHw4lCVi/a4OaALpKTSBKP/leSLycQJ5dNk+4iNDuWZa/699bdqAJABW7i2zvR0r95Yxwt37aS0kOIhbzhzI/sO1vL/tgO3tOKm+Y+DQdmhp9m07lCM0gSj/dWAzBIUcG1v3gbyyGpZuO8A1U/sTFRbyta8P7B1Nckw4q/bZuyKrobmF9fvLjyas450/si/9ekXy3+X7bG3HKfUdDc311uZP1e1pAlH+68Bma2loaITPmvDE57mEBAkLpme1+XURYeqAXqzcW2brPMim/Eoaml1MzW57NVpwkHD96dmszStnw/5y29pxSp75Kh3G6hE0gSj/dWCzT4evKmubeGltPhePTSMl7sRJbNqAJIor69l/2L55kJV7yxCBKSdIIADfmNSP2IgQ3/ZCkodAcBgc+NJ3bVCO0QSi/FNNGRwp8mkCeW71fmobW7hxxsnLyE9zP6nbOQ+ycm8Zw/vGkRAVdsJrosND+NbU/ry7pZh8G5PZSYWEWUOOB7f4Jr5ylCYQ5Z8O+nYCvbnFxcIvcjl9UBIj0uJOeu2gPjEkRYexaq898yCnmv9o7brpWQSJsPCLXFva0i59x0Dxl1rSpAfQBKL8k2cMPcU3CeSLPWUcqKrn21MzT3mtiDBtQJJt8yBfFlRS3+Ri2oBT78ZPjY/k3OEpvLGxkOYWH5U46TsKaku1pEkPoAlE+acDmyE2DaJP/arbDm9sKCQuIoSzh/Vp1/XTBvSiqLKe/MN1Xd6WlXtOPf/R2mUT0imtbuTzPfYvLW6Tp9eoE+ndniYQ5Z8ObPHZ8FVtYzNLth7gwtGpRIQGt+s+U23cD7JyXxnDTjH/0dpZQ3sTFxHCGxsKu7wt7aIrsXoMTSDK/zTVW1V4+/qmhMn72w5S29jCpePT232fwX1i6BUdxsp9XZtAGptdrMsrb9fwlUd4SDBzx6SxZMsBahp8sKEvMgHi+2sC6QE0gSj/U7IDXM0+64G8vqGQtPgIpmS1/0lbRJiYmciXBZVd2pbdh45Q3+RiYmZih+532fh06ppaeH+bj+YhtKRJj6AJRPmfoyVMxjgeurS6geW7S5k3Pp2gIOnQfTMSIymuqOvSifSiinoA+iVGdeh+kzITSU+I5HVfDWP1HW3tRm+s8U185QhNIMr/HNwCodGQmO146Lc3FdHiMlzWgeErj7T4SGoaW6iq67pho+JKa1I+NaFju/GDgoRLx6exfHcJJUcauqw97dZ3NGCsuliq29IEovxP8ZeQMhKCnP/zfH1jESNS4xiSEtvh+3qe5Isqu24lVmFFHaHBQnJ0eIfve+m4dFwG3tpU1GXtaTfP/FXxJudjK8doAlH+paUZijdC+gTHQ+eV1bApv4J549K8un9qfCRwrNfQFYor6ukbH9Hh4TSAwSmxjEyL401fJJCETIjsBYXrnY+tHON4AhGR2SKyU0RyROT2Nr4eLiIvur++SkSy3LdPEZGN7rdNInKZ021XDji0DZpqIWOy46Hf2VwMwNwxqae4sm1pnh6Ie96iKxRX1h1NTN6YOyaVTfkVFJQ7XNpExPodFqxxNq5ylKMJRESCgQeBOcAIYL6IjDjushuAcmPMIOA+4F737VuAScaYccBs4GER+Xp9bRXYPE84GZMcD714czFj+yWQ0cEJa48+sREEB0mX9kCKKupJT+hEAhltJcN3N/vgnJCMydZy7LoK52MrRzjdA5kC5Bhj9hpjGoEXgHnHXTMPWOj++BVgloiIMabWGOOZnYwAtNBOd1SwBqJ7W0MgDsorq2FLYRUXjfau9wFWSfW+cREUd1EPpMVlOFhVT2q89+XsM5OiGZUed7R35ah+7l5k4TrnYytHOJ1A0oH8Vp8XuG9r8xp3wqgEkgBEZKqIbAU2A7e0SihHichNIrJWRNaWlJTY8BCUrQrWQMYUawjEQZ4n2Dmj+3bq+6TGR3TZJHppdQPNLkNqJ3ogABeOTmVjfgWFFV1fZuWk0iYAosNY3VhATaIbY1YZY0YCk4FficjXXpoZYx4xxkwyxkzq3bu3841U3qs9bO0dCMDhK4/UhEiKK7umB1LkfsJP60QPBFoPYzncC4mIgz7DNYF0Y04nkEKgX6vPM9y3tXmNe44jHvhKfQhjzHagGvBNrQtlj4K11nuHJ9D3l9WypbCKuZ3sfYD1ZF9cWY/L1fkRVk8i6swkOljDWCPT4nj7Sx8MY2VMtn6vLh9VBla2cjqBrAEGi0i2iIQBVwOLjrtmEbDA/fGVwEfGGOO+TwiAiGQCw4BcZ5qtHFGwBiQI0sY7Gvbo8NUo7+c/PFLjI2hsdlFW09jp73W0B9LBTYRt8dkwVsZkqK+Aw3ucjasc4WgCcc9Z3AYsBbYDLxljtorIXSJyifuyx4AkEckBfgp4lvrOADaJyEbgdeBWY0ypk+1XNitYY20gDI9xNOw7m4sY2y+Bfr06N3wFHJ2v6IqVWMWV9USGBhMfGdrp7+WzYSxPbzJ/tbNxlSMcnwMxxiw2xgwxxgw0xvzRfdsdxphF7o/rjTFXGWMGGWOmGGP2um9/2hgz0hgzzhgzwRjzhtNtVzZyuazVOgE8fAVWORPomr0gxZV1pCZEIF2woCAr2RrGcnw1VvIQCI/XeZBuKqAm0VU3VroTGqocTyBdOXwFx8qZdEUPpKii/mhC6goXjk5lw/6Ko0NjjggKgoyJx+a3VLeiCUT5h6MbCJ1NIO9uKWZsRnyXDF8BJEWHERYS1CUrsaxd6J2f//C40DOMtcXhTYUZk+HQVmiodjausl2HEoiIfEtEXhCRZ0XkORGZb1fDVA9TsAYiEiBpkGMh8w/X8mVBJXM6sXnweCJi7QXp5Kv8xmYXh440dHoPSGvZydEMT41jsePzIFPAuKBog7Nxle062gM50xhztTHmGmPMt7AmtpXqvPw11v4PBzcQep5I53ZhAgFrHqSzPZCDVfUYA+ldsAKrtbmj+7Iur7xLy62ckqcwZv4q52IqR3Q0gYSLyFwRGSMiFwJd9/JI9VzVJVCyHTKnOxp28eZiRqd33fCVR2pCBMWd7IF01R6Q413oi9pYUb2gzwjI+9y5mMoRHU0gtwKJwIXu97d1eYtUz5O73HqffaZjIfMP17KpoPLoE2pXSouP5OCRBppbvN885+khdMUekNYG9I5hWN9Y54exsmdC3gpo7vz+GOU/OpRA3AUNnzHG/NkY86wxxuEa0apb2rcMwmIhdZxjIZe4J5K7evgKrB5Ii8twqBMnAXqWAXd1DwSsXsjavHIOdFHJlXbJngnNdVCoq7G6k06twhKR60XkrC5qi+qpcpdD1ukQ7Fx1/nc2FzMqPY7+SV07fAXH9oJ0Zp6huLKOuIgQosO7/mfi6XUt2eJgLyTzdKvKwL5lzsVUtuvsMt63gY6ftamUR2WhVUAx6wzHQhZW1LExv6LL9n4cL7ULDpYqqqgnrQtXYLU2qE8MQ1NiWezkPEhkAvQdowmkm/EqgYhIHIAx5pAxZmnXNkn1KEfnP2Y6FvKdL60jXi/y8uTBU+mKo227eg/I8eaOSWVN3mFnV2Nlz7SWazfqyHd30eEEIiLPA38RkYUi4tx/veqe9i2DyERIca6w8lubihmTEU9mUrQt3z8uIoTosOBO9UCKK+u7dA/I8S4ak4ox8I6TFXqzz4SWRl3O24140wPZboy5xRizAKuarlLeMcZKIFlnWCUvHJBbWsPmwkouHpNmWwwRcZ8L4t2r+/qmFg7XNHb6HJCTGdA7xvkS7/2nQVCIDmN1Ix3dif4ocL6I/NA9ea6HHSvvledCZb6jw1dvu4ev5to0fOWR6j4XxBt27QE53sVj09iYX0H+YYeGlMJjIH2SJpBupN0JRESCgO8D52Kd6zEW91GzSnnF80TiYAJ5a1MxkzITbZug9kiNj/B6mezBKut+fW3sgcCxJcxvuZOqI7LPgKL1UF/pXExlm3YlEBG5DTgI5AGrgZHGmPuNMTfb2TjVze1bBjEpVslvB+w6eISdB49w8Vj7hq88EqLCqKxr8uq+FbVN7u/R+XNATqZfryjG90/g7U1OzoPMtOpi5a1wLqayTXt7ID8DRhtj0oELgNNF5E77mqW6PZfr2PyHQ/Wv3t5URJDAnC46++Nk4iNDaWh2Ud/U0uH7VtZZu7UTosK6ullfc/GYNLYVV5FzyKFKuRlTIDgc9n3qTDxlq/YmkGrgEIAxphi4AbjCrkapHuDAJqg5BIPOdSScMYa3vixm2oAk+sTaOzQEHD1F0JteiOc+XXES4anMHZOKyLG5IduFRkDWDNj9njPxlK3am0D+A7wsIp5a2/0BXcytvLdzCSAw+HxHwm0tqmJfaY0jw1dwbPjJMxzVERW1TYQECdFhwV3drK9JiYtgSlYv3tpUhDHG9ngADJ1jbR4tzXEmnrJNuxKIMebfwHPAf0WkHMgBdorIVSIy2M4Gqm5q1xLoNwWinVmH8fqGQkKDhTmj7B++gs73QOIjQ7vkKNv2mDcunT0lNWwprHIk3tEXDbuWOBNP2eakCUREThP3X7Ex5lVjzFlAb2Ai8BEwHXjY7kaqbqaqGIo3wpALHAnX3OLizY1FnD20jyPzCgAJkVYcbxJIRV0T8TZPoLc2d3QqYcFBvLahwJmAiZlWeXdNIAHvVD2Qa4F17lMIrxORvsaYZmPMJmPMQmPMT4wx5zjRUNWN7HZXvxkyx5Fwn+WUUlrdwOUT0h2JB8d6IBW1HS9fXuXugTglPiqUc4b14a1NRZ0qQd8hQ2bD/hVQp1vJAtlJE4gx5vvGmAnA77HO/3hSRFaIyD0iMlNE7B+kVd3PrqUQ3x/6DHck3BsbComPDOXsYX0ciQddM4TlpMsmpFNa3cjynFJnAg6ZDa5m2POhM/GULdo7B7LDGHOfMWY2cA7wGXAVoEVtVMc01cGej63hKwfG+Gsamlm69SBzx6QSHuLc653YiBBEvBzCqm0iweEEYg3vhfL6+kJnAmZMgqgk68WEClgdPmzAGFMnIp8ZYxbb0SDVze1bbh0sNHS2I+GWbDlAXVMLl493bvgKIChIiIsIDZgeSFhIEBeNSeWVdQVUNzQTY8M5JF8RFGxNpu9aAq4W63MVcLQar3LWriUQGg2ZMxwJ9/qGQvr1imRiZqIj8VpLiArt8DJel8tQVd9EvEOT/a1dNj6D+iYX7zp13O2QC6CuHPJXOxNPdTmtxqucY4w1ZDHwbGtDmc0OVNbz+Z5SLhuX7tiS2NbiIzveAzlS34wxzmwiPN6E/glkJUXx+gaHhrEGzrKq8+5615l4qstpNV7lnMJ1UFUAQy90JNzrGwoxBi6bkOFIvON5k0AqPGVMfJBARITLxmewYm8ZBeUO7BOOiLNqY21bZL24UAFHq/Eq52x5FYLDYNhc20MZY3hxzX6mZPUiO9meg6NOxZsE4mQZk7ZcMdGaK3p5rUN7QkZeDuX7oGiDM/FUl9JqvMoZrhbY8hoMOs86H9tmq/YdJreslm9O7md7rBPpVAJxcCNhaxmJUZwxuDcvr82nxeVAr2D4RRAUar24UAFHq/EqZ+xfAdUHYNTljoR7cU0+seEhXDja3oOjTiYhykogHakxdbSUu496IADfnNSPosp6lu8usT9YZKJVUHPr61aFZhVQtBqvcsaWVyE0yiqkZ7PKuiYWby5m3vg0Ih0oSHgi8ZGhtLgM1Q3N7b6Pr4ewAM4d0Yde0WG8tDbfmYCjroCqQj0rPQBpNV5lv5Zm2PamlTzC7J+PWLSxkIZmF1dP7m97rJPxph6W59o4HyaQ8JBgLh+fzvvbDlJa3WB/wKFzICRSh7ECUEeq8T6LVuNV3tj3KdSWWa80HfDCmnxGpMYxKj3ekXgnEhfZ8ZLulXVNRIQGERHq241135zcj6YW48zO9PAYa0/ItjesFxsqYLR7FZYx5rVW1XgnoNV4VXtteQ3C4x05PGpLYSVbi6q4eorvJs89PGeCVHWgB1JR23i05+JLg1NimdA/gRfW7HfmnJBRV0BNCeQutz+W6jId3kjorsb7pVbjVe3SVA/b37JW24SE2x7u2VV5RIQGMW+ss6VL2nK0Im8Hh7B8Of/R2tVT+rOnpIZV+w7bH2zweRAWC1tesT+W6jLe7ETvFBGZLSI7RSRHRG5v4+vhIvKi++urRCTLfft5IrJORDa732viCgQ73oaGShh9pe2hKmobeX1DIZeOS/fZMtjWvKnIW1HrPwnkkrFpJESFsvCLXPuDhUbCiEtg6xvQ4ND57KrTHE0g7vLvDwJzgBHAfBEZcdxlNwDlxphBwH3Ave7bS4GLjTGjgQXA0860WnXK+oWQkAnZZ9ke6uW1BdQ3ufjOaZm2x2oPzxBWRyfR/SH5AUSEBvPNSf14b9tBiirq7A844VporIatr9kfS3UJp3sgU4AcY8xeY0wj8AIw77hr5gEL3R+/AswSETHGbDDGFLlv3wpEioj9YyLKe4f3wr5lMOE7EGTvn1qLy/D0yjwmZyUyMs23k+cekaHBhAZftQc5AAAgAElEQVRLhybRnT5M6lS+PS0TlzE8t2q//cH6TYXkobBu4amvVX7B6QSSDrReXF7gvq3Na4wxzUAlXy+ZcgWw3hjztTWGInKTiKwVkbUlJQ5shFIntv4pkCAY923bQ32y8xD7D9eyYHqW7bHaS0SIjwzr2BBWnfNngZxMv15RzBqWwvOr91Pf1GJvMBGYuAAK18LBrfbGUl3C8TmQzhKRkVjDWm2WUTHGPGKMmWSMmdS7d29nG6eOaWmCDc/C4Asgzv7d4AtX5NEnNpwLRva1PVZHxEeGUFnXvmNtG5td1Da2+FUPBOC66VmU1TSy2Iky72OutuqlaS8kIDidQAqB1usrM9y3tXmNiIQA8UCZ+/MM4HXgWmPMHttbq7y3awnUHLJeUdpsb0k1y3aVcM3UTEKD/es1UUJU+3sgnusS/GQOxOP0QUkM7B3NwhV59geLToLhF8OXL1inVyq/5vR/2xpgsIhki0gY1nkii467ZhHWJDnAlcBHxhgjIgnAO8DtxpjPHWux8s66hRCbahVPtNmTX+QSGizMn+r7vR/Hi49s/6FS/rALvS0iwoLpWWzKr2BdXrn9AScsgPpKq8y78muOJhD3nMZtwFJgO/CSMWariNwlIpe4L3sMSBKRHOCngGep723AIOAOEdnofuvjZPtVO1Xsh5wPYPy3Idjeo1HLqht4aW0+l45Lp0+s/YdUdVRHKvJ6hrr8bQgL4MqJGSREhfLwpw50/LPOgMRsWPek/bFUp9h88PHXuc9SX3zcbXe0+rgeuKqN+90N3G17A1XnrXrYmjyfYP/w1VMr8qhvcnHTzAG2x/JGxxKIZwjL9zvRjxcVFsK10zL518c55ByqZlCfGPuCBQXB5Bvgvd9a54SkjbcvluoU/xowVoGvvtIavhp5GSTYO6RU29jMUytyOXd4HwanxNoay1vxkaEcqW9u19ka/lCJ92SunZ5FWHAQ/12+1/5gE661dqaveND+WMprmkBU11r/FDQegem32R7q5bUFlNc2ccuZA22P5a2O1MPyh7NATiY5JpyrJmXw2vpCDlXV2xssIt5agLHlNah06HRE1WGaQFTXaWmClQ9B5gzbhx2aW1w8unwvEzMTmZTVy9ZYndGRelj+Oone2o0zBtDscvGEE+VNprpX6q96yP5YyiuaQFTX2fYmVBU40vt4Z3MxBeV13Oyncx8eHSlnUlHbRGxECMFBYnezvJaVHM2cUak8szKPI/UdO663wxL6w4h51pBofZW9sZRXNIGormEMfPEvSBpsbR60UYvL8MBHOQzqE8O5w1NsjdVZHSmo6G9lTE7kljMHcqS+2Zkii9Nvg4Yq2KCl7/yRJhDVNXKXQ/FGOO1W2+tevbO5mN2HqvnRrMEE+fGrdWg1hFV76t3oFQGSQEZnxHPu8D48smwvVXb3QtInQv/psPI/0Ny+Hf3KOZpAVOcZAx//ydo4OPZbtoZqcRnu/2AXQ1JimDva/hIpnRXvPhyqPZPolXVNfrcL/UR+fO4QquqbefyzffYHm/kzqMzXXogf0gSiOm/vJ7D/CzjjZxBq72a+tzYVsaekhh+fO8Tvex/QugfSvgQSCD0QgFHp8VwwMoXHPttHZQeqDXtl4CyrUu/y/7MOKFN+QxOI6hxj4OM/QlyGtXbfRs0tLu7/cDfD+sYy28+KJp5IWEgQUWHB7Z5Ej/eD42zb68fnDuFIfTOPfWbzvhAROPvXUFVoLRNXfkMTiOqcnA+gYA3M/F/bj6x9Y2MR+0oDp/fh0Z7d6MaYgJlE9xieGseFo/vy+Oe5lNfYPD+RfSZknu7uhWiRRX+hCUR5z9P7SOgP466xNVR9Uwv3vb+LkWlxXDDSv1deHS8+MvSU+0DqmlpobHEFzByIx4/PHUJtYzMPfpxjbyBPL6T6AKx93N5Yqt00gSjv7XjHqlU08xcQYu/Qy5Nf5FJYUcdvLhyOSOD0PqB9PRB/L2NyIkNSYrlqYj8Wrshlf1mtvcGyZlg9keV/130hfkITiPJOc4NV7K73MBg739ZQZdUNPPhRDucM68P0Qcm2xrJDfGToKSeaPZPsgZZAAH56/hBCgoK4d8kO+4Od+3uoLbWGspTPaQJR3ln1EJTvgwvusb1k+z8/3E1tUwu/vnCYrXHskhDV/h6Iv9bBOpmUuAhuPnMA72wutv+8kPQJ1lLxlf+Gww4UdVQnpQlEdVz1Ifj0rzBkNgyaZWuoPSXVPLtqP1dP7segPv5ZcfdUrDmQk08yB0IdrJO5aeYA+sSG88d3tmHMqSsPd8qsOyAoFN77nb1x1ClpAlEd99H/g+Y6ON/+41n+tHg74SFB/PjcIbbHsktsRCj1TS6aWlwnvOZIfTMAcRGBmUCiwkL42flDWL+/gre+tPns9LhUOOMnsONt2LfM3ljqpDSBqI4p3gTrn4YpN0PyYFtDvb/tIB9sP8T/zBpM71h7lwjbKTbCGuKrdieJtngKE3quDURXTuzHqPQ47n57m/2FFk+7DeL7w5JfQcuJf67KXppAVPu1NMOiH0J0Mpz5C1tD1TY28/tFWxmSEsMNM7JtjWW3WHev4shJE4j1tZgATiDBQcIfLx1NSXUDf39/l73BQiPhgj/CwS2wUg+d8hVNIKr9Vv7bKph44V8hMsHWUP/8MIfCijruvnQ0ocGB/WcaE24lhZMVHjxS30REaFDAP9ax/RK4Zmp/Fn6Ry5bCSnuDDb8Yhs6Fj++BMgfOaldfE9h/rco5ZXusTYND58KIS20NtfvgEf67fC9XTMhgSrb/HhbVXnGeIayGE/dAqhuaj/ZUAt3Pzx9Gr+gwfvvGFlztOMrXayIw928QHAZv/cja2KocpQlEnZox1j9ocJj1D2vjRj6Xy/Cb17cQHR4SsMt2j9eeIayq+uaAnv9oLT4qlF9fOJyN+RU8u3q/vcHi0uC8u6zjBNYvtDeW+hpNIOrU1j1h/YOed5f1D2ujJ7/IZXXuYX5z4XCSYgJ34rw1T2I42cTykfru0wMBuGx8OjMGJfOnxdvJP2zzDvUJC6xjlN+7Q89Pd5gmEHVyJbtgya9hwFnWP6qN9pZU85elOzh7aG+umpRhaywnxRxNICdfhRUb3j16IAAiwr1XjiFIhJ+/ssneoaygILjkn+BqhtduAleLfbHUV2gCUSfW3ACvXg9hUXDZw7aeNNjiMvz8lS8JCw7iz1eMCbh6VycT2545kG40hOWRnhDJ7y4azsq9h3l6ZZ69wZIGWsOreZ9btbKUIzSBqBP74A9wYDPMexBi7T1/47HP9rIur5w/zBtJSpy9h1I5LTwkmLCQoFOswup+CQTgG5P6cdbQ3vz53R3kltbYG2zsfBh1BXzyJ8hfbW8sBWgCUSey+31rff2Um2DoHFtDbSms5G9Ld3HeiBQuHZduayxfiYsIOfUQVjeaA/EQEf58+RhCg4UfvbCBxuYT78bvgmBw0X0Qnw6v3gB1FfbFUoAmENWWw3utf8CUUXDe/7M1VHVDM7c9t55e0WHc282GrlqLCT9xAmlxGWoaW47uF+lu+sZH8Jcrx7CpoJK/LrW5Ym9EPFzxGFQVwWvf0/kQm2kCUV/VUA0vXAMIfPMZW884N8bw69c2s/9wLf+cP55e0YFznGtHxUaEUn2CISzP3Eh3HMLymD0qlWtPy+TR5fv4cPtBe4P1mwKz/wy737M2GSrbaAJRxxgDb94KJTvgqiegl70lRF5am8+iTUX89Lwh3WLD4MnEnmQIy7O8N1ALKbbXry8czojUOH728iaKK20+lnbyjTD+O7D8b7DtTXtj9WCaQNQxy9z/bOf+AQaeY2uozQWV3PHmVk4flMT3zxpkayx/cPIE0v17IAARocE88K3xNDW7+P4z66lvsnF4SQTm/h+kT4LXv28tBlFdThOIsmx4Fj6+G8Z8E6b/j62hSo40cNPTa0mOCef+q8cTHNQ95z1aiwkPPeFGwu5QSLG9BvSO4f++MZaN+RX89o0t9p4dEhJuDcNGJsAzV0KFzbvieyBNIAp2fwCL/sfaLHjJA7aWKmlobuH7z6yjvLaRR66dSHI32W1+KrERIRw5wT6Q6gZPKffuPYTlMXtUKj+aNZhX1hXwxOe59gaLS4VrXoGmOiuJ1B62N14PowmkpyvaAC9dC31GwDeehhD7JrKNMdz55lbW5pXzt6vGMjIt3rZY/iYuIoTqhuY2d2T3lCGs1n40azAXjEzhj4u3s3x3ib3BUkbA/OesI5ifn28lE9UlNIH0ZAe2wNOXQ1QSXPMyRMTZGu7fn+zhhTX5/ODsgVw0xt6aWv4mNiIUY6Cm8eu9kKoemECCgoS/f2Mcg/vEcOsz69lWVGVvwKwZcPkjkL8KXvw2NNXbG6+HcDyBiMhsEdkpIjkicnsbXw8XkRfdX18lIlnu25NE5GMRqRaRB5xud7dzcBs8dQmERMCCN62uvo1eXpvPX5fu5NJxafzsvKG2xvJHJ6uHdfQ0wvCeMYTlER0ewuPXTSYmIoTrnlhNQbnNRRdHXgaX/AtyPrB63c0N9sbrARxNICISDDwIzAFGAPNFZMRxl90AlBtjBgH3Afe6b68Hfgf8r0PN7b5KdlrJIygUrnsbeg2wNdzHOw9x+2ubOWNwMn+5cixBPWDS/Hgnq4dVXd9MSJAQEdrzBgTSEiJZeP0U6ptauPbx1ZTXNNobcMJ34KJ/wO6l8NICaLY5Xjfn9F/sFCDHGLPXGNMIvADMO+6aeYCnsP8rwCwREWNMjTHmM6xEorxVtBGeuBAkyEoeSQNtDbcm9zC3PrOeYX1j+c+3JxIW0vOeJKH1mSBfX4nlqYPVXXfhn8qQlFj+u2AyBeV1XPfkGvvPU5/0XWuJ76534YX50Ghzja5uzOn/5nQgv9XnBe7b2rzGGNMMVAJJ7Q0gIjeJyFoRWVtSYvPkXKDZtwyevAhCo+C6xZA82NZw6/LKue7x1aQmRPDEdyd321Id7eHpgVSdYAirp6zAOpEp2b148FsT2FpYyXVPrDlp5eIuMflGazhrz0fw1KW6OstL3e7loDHmEWPMJGPMpN69e/u6Of5j+1vwzBVWobkblkKyvZv3NuwvZ8Hjq+kTF8Hz35tGn9juVWG3ozxnfbQ9B9Lco5Orx3kjUvjX/PFszK/g+ifWUNvGgoMuNeFauGohFLt75ZWF9sbrhpxOIIVAv1afZ7hva/MaEQkB4oEyR1rXHRkDn/0DXvwOpI6F775r+6mC6/IOc+3jq+kVHcZz35va7cqze8PTw6huK4E0dM9S7t6YMzqVf3xzHGvzDnPdEw4MZ424xNonUlkAj54DhevtjdfNOJ1A1gCDRSRbRMKAq4FFx12zCPAcfXcl8JGxdbtqN9bcAG/+AD64E0ZeCtcugih7a059uquEb/93NUnu5JEaH2lrvEBxsmNtu9txtp118dg0/nH1eNbnlfOtR1dRVm3zaqkBZ1q98uAwqyey9XV743UjjiYQ95zGbcBSYDvwkjFmq4jcJSKXuC97DEgSkRzgp8DRpb4ikgv8HbhORAraWMGlPKqKYOElsPFZOPN2uPIJ62RBG721qYgbF64hOzmal2+ZTkaivfECSVRYMEFy4mW8cdoD+YpLxqbxyLUT2XXwCFc9vILCCps3/6WMhO99BKlj4OXr4KO7tRR8Ozg+B2KMWWyMGWKMGWiM+aP7tjuMMYvcH9cbY64yxgwyxkwxxuxtdd8sY0wvY0yMMSbDGLPN6fYHhL2fwENnWAXkrnwczv6VreVJjDH8d/lefvjCBsb1S+D5m6bRO7ZnlChpLxFxnwnSdg+kJ9TB6qhzhqXw9A1TKalq4Ip/f8HWokp7A8b0hgVvwbhvw7K/wtOXQbUuxDmZbjeJ3qO5WuDTv1irSqKT4aaPrSM+bdTU4uI3b2zh7ne2c/6IFJ66firxkToc05bYiNCv1cMyxlCtcyAnNCW7Fy/efBoicNVDK3h/m81niYSEw6UPWjXh8lfBw2dA7uf2xgxgmkC6i8P74Ik58PEfYcw3rO54b3t3fFfWNnH9k2t4btV+bjlzIP+5ZiKRYcG2xgxkbZV0r2tqocVldA7kJEakxfHmD05nUJ8Ybnp6LY8s22NvFV+wNhze+IG15P3JufD+nbpzvQ2aQAKdMbD+KXhoBhzaAZf/Fy57GMKibQ27pbCSix5Yzsq9ZfzlyjHcPmdYj9xh3hFxEV8v6d4TCyl6o09cBC/edBqzR/blnsU7+J/nN1Bj916RvqPh5mUwcQF8/g94dBYc3GpvzACjCSSQHd4HT19qlWJPGw/f/xzGXGXrfAfAi2v2c/l/vqCp2fDCTafxjUn9Tn0nRUwbPRBPQtF9IKcWGRbMg9+awC9nD2Px5mLmPfg5OYeO2Bs0PAYuvh/mvwDVB+DhM+GjP2pvxE0TSCBqaYLP/wn/Pg0K1lllGa5dBAn2PpFX1Tfxkxc38stXNzM5K5F3fjiDiZmJtsbsTmLdJd1b8ySU7n6cbVcJChK+f9ZAnrlhKuU1jVzywOe8uGa//UNaQ+fArausOcVlf7F6/Do3ogkk4Oz5GP5zOrz/Oxh4NvxglVWWIcjeX+Xa3MNceP9yFm0q4sfnDuap66eS1EMOg+oqbc2B6BCWd6YPSuadH57B2IwEfvnqZr7/zHr7CzFGJ8HlD8O3X7XKwT95Ibx6o7VkvofSBBIoyvbAC9dYQ1YtDVaX+urnrNIkNqpvauHeJTv4xsMrCBLhpZtP48fnDukRx9B2tVj3HEjrV8vHEoj2QDqqb3wEz944lV/NGcaHOw4y+/5lfGD3Ki2AQedaL9xm/gK2LYJ/TbKW/TbaXI7eD2kC8XfVh+Cdn8GDU6zexzm/s7rSQ+fYPtexep/V6/jPJ3u4amI/Fv/oDB2y6oSY8BCaWgwNza6jtx2dA9EeiFeCgoSbzxzI67eeTmJUGDc+tZYfPr/B/t3rYVFwzm/gttXWSMBHd8O/JsC6J6HF5sl9P6J/tf6qpgxW/AtWPWL1OCZeZ73iiU2xPXR5TSN/WbqT51fvp1+vSJ69cSqnD0q2PW53F9fqUKmIUGu5s2dORIewOmdUejyLbpvBfz7ZwwMf72b57hJunzOMqyb2s3d1YGIWXP0s5H1hLfV960fwxQNw5i9h1OUQ1L2Xtetfrb+pLoEVD8DqR6Gp1jpF7Zzf2n5uB0CLy/D86v387b2dHKlv5oYZ2fzs/CFEhemfSVdofSaIZ6d+VX0zIhCjP+NOCwsJ4kfnDmbO6L78+rXN/PLVzTy/Op+75o1kTEaCvcEzp8MN78GOt+Hje+C1G+HTe2Hmz62J9+Du+fvtno8qEJXtsRLHhmehpdH6o5v5c+gzzJHwy3eXcM/iHWwvrmJqdi/umjeKoX1jHYndU8S2caztkfomYsJCdA9NFxqSEsvLt5zG6xsKuWfxDuY9+DmXj8/gZ+cPIS3BxuKeIjD8Yhg6F3a8BZ/cC6/fZA1vTb8Nxn/b9v1ZTtME4kvGWIc8rX4Edi6GoBAYOx+m/9D28zo8thRWcu+SHSzfXUq/XpH8a/54LhqT2mNPx7NTTBtngmgdLHuICJdPyOC8ESk88HEOT3yey9tfFvHd07P5/pkDiY+ycdFCUBCMmAfDLraOzv3sH/DuL+CTP1lnkEy+ERL62xffQfqX6wt1FbD5ZVjzGJRsh6gkmPETmHKzI3McANuLq/jHB7tYuvUg8ZGh/HbucL5zWibhId17zNaXjp4J0nBsN3p1vdbBslNsRCi/mjOc70zL5O/v7eKhT/fw7Mo8rp+RzfUzsu2t2xYUZC12GToH9q+EFQ/CF/+y3oZeaB2tO+Ac25fg20n/cp3icsH+FbDhGeu8geY6SB0H8/5tDVeFOnPo0qb8Cv7zyR6WbD1AbHgIP5o12P5/JAW0faztkQY9ztYJGYlR/P2b4/jezAHc/8Fu7v9wN098vo8F07O4bnqW/Xua+k+z3iryYe1jVvmhHW9DfH+r7tbY+bZvBLaDJhC7le2BL1+ETS9ARR6ExcK4+TBhAaSNc6QJLpdh2e4SHv50Lyv2lhEXEcIPzxnEDTMG2NuVV18Rd3QS/atDWL2iw3zVpB5neGocD31nIluLKvnnh7t54OMcHl2+l29M6scNM7LJTLJ5jiKhH5z7ezjrV7DjHVi/0CqA+vEfIesMK5EMvwgi4u1tRxfRBGKH8jzY9gZseRWKNwECA86Cs39j/XE4NJFW3dDMa+sLePKLXPaW1NA3LoLfzh3O1VP6a+0lH4gOt4YHWxdUPFLfTP9eevCW00amxfPwdyaxp6SaRz7dy/Or9/P0yjxmDevDgulZzBiUbO88YEi4tcx31OXW88WXL8Km5+HNW+Htn8Dg86yRicHnW/W4/JQ+i3QFY6Bkh9Ul3f6WO2kA6RPhgnuspbg2n0Pe2uaCSp5bvZ9FGwupaWxhXL8E7r96HHNGpRIWErjjrYEuJDiIqLDgr5yLrsfZ+tbA3jHce+UYfnb+EJ5dtZ9nV+XxwWOrGZAczTcn9+OKiRkk2z28lZgJZ/7CWnVZsAa2vGYNc+94G0IiYOAs64XnkNm2H0ndUZpAvNVUB3mfw673YNcSa3gKIGMKnPf/rOV8vbIda05pdQOLNhbx6voCthZVEREaxEVj0rhman/G99fd4/7i+HpYepytf+gTF8FPzhvCrWcPZPHmYp5flc+f3t3B397byaxhKVwxMYOzhvYmNNjGF2Ai0G+K9XbBH6050+1vWW873wEJsp5fhs62eiZ9RthejeJU9C+3vVwuOLQN9n0Kez6C3M+gud56hTDgLGsV1ZDZEJfqWJOqG5p5f9sB3tpUzKe7SmhxGUanx3PXvJHMG5euE+N+yDqV0BrCamx20dDs0lVYfiQ8JJjLxmdw2fgMcg4d4fnV+byxoZAlWw/QKzqMi8akctGYNCZlJtq7dycoGLJmWG+z/wxF62HXUtj5Lnzwe+stNtXqnQw8G7JnQkwf+9pzAvqXeyLGQMlOyPvMKtu8bxnUllpfSxoME78Lg2ZB5ulWXRyHVNY18fGOQyzZcoCPdx6iodlFanwE3ztjAJdPSGdIim7+82fWuehWD0TPAvFvg/rE8ruLRnD7nGEs313Cq+sKeXFNPk+tyKNvXARzRvflgpF9mZzVy97ioiLWcHj6RDj711b135wPYc+H1obFjc9Y1/UZYU3EZ51uPS9F219+SP9y21K2Bx47/1jCiE21kkX2mVamd3i5XV5ZDR/tOMRHOw6xYk8ZzS5D79hw5k/pz0VjUpnQ3+ZXQ6rLtB7COlYHS3uK/iw0OIhzhqVwzrAUqhua+XD7Qd7aVMyzq/bzxOe59IoOY9awPpwzrA8zBifb//uMS7OW/k74DrhaoHij9QJ376ew4WlY/bB13eAL4JqXbG2KJpC2JPS3xhn7TbOyeWK2o2ONtY3NrNp3mOW7Svl01yH2lNQAMCA5mhtmZHP+yL6M75egSSMAxUWEUlRRB+hZIIEoJjyEeePSmTcuneqGZj7dWcLSrQdYsvUAL68rIDRYmJTZi5lDenPG4GRGpMbZP9Tl6Z3M+Ak0N1oJJfczCLWxbIub/uW2JTgU5j3oWLj6phY25lewYk8ZK/aWsXF/BY0tLsJDgpiS3YtrpmZyzrA+ZCV3rzo6PVHrHkiVewhLeyCBKSY8hLljUpk7JpWmFhfr8sr5eMchPt1Vwr1LdnDvEkiKDmPagCSmDUzitAFJDOwdbfPy4LBjE/EO0ATiA4drGlmfV866/eWs2XeYLwsqaWxxIQKj0uK57vQszhiczOSsXkfLfqvu4atzINoD6S5Cg4OsRDEgiV9dOJyDVfV8truUz3NKWbG3jHc2FwNWQpmUlcjkrF6M75/IqPS4gC4fpH+5NqtrbGFbcSWb8iv5sqCCjfkV5JZZJ5eFBAmjM6yEMTmrF1Oye+nKqW4uNiKUuqYW5ty/nKo6Tw9E/w27m5S4CK6YmMEVEzMwxpBXVsvKvWWsyS1nTe5hlm61Tk4MCw5iZHocYzMSGJMRz5iMBAYkRwfM8LTYfhi9D02aNMmsXbvWsXglRxrYcaCKHcVH2FpUydaiKvaUVONy/4hT4sIZm5HAhMxEJvRPZHR6PJFhgfvqQ3XcjgNV3P/BbprdfxTJMeHcfekoPSK4hzlUVc/6/RVs2F/Ohv0VbC6spK6pBYDosGCGp8YxKj2e4amxDOsbx5CUWEefK0RknTFm0imv0wTSMcYYDtc0knOompySanYfrGbXwSPsOniE0urGo9elxIUzKi2ekenxjEqLY2y/BFLinCmYqJQKLC0uQ86hajYVVLC10Hrxua24itpGK6kECWQmRTMkJYahKbEM7BPDoD4xDOwdY8swd3sTiPadT+BIfRN5ZbXkltWQW1rD3tIa9rnfKmqP1TKKCgtmcEosZw/tw9C+sYxIjWNo31j7q3sqpbqN4CBhaN9Y6xC3SdY2AZfLsP9wLduLq9h+4Ai7Dhxh16EjvL/t4NFRDRFIT4gkOzmaAcnRZCVHk5UUTWZSFP16Rdm7cx7tgbQpr6yGM//6yVdu6xsXwYDe1i9oYG9P9o8mLT4yYMYrlVKBr76phdyyGvYcqiHnUDV7S6vZV1rD3pKao3uLAMb2S+DNH5zuVQztgXRCWkIkt88ZRlZSNFnJUfTvFaXngiul/EJEaDDD+sYxrG/cV243xlBW00heWQ25pbWOrODUZ8U2hAYHccuZA33dDKWUajcRITkmnOSYcCZmOlO1V2t7K6WU8oomEKWUUl7RBKKUUsormkCUUkp5xfEEIiKzRWSniOSIyO1tfD1cRF50f32ViGS1+tqv3LfvFJELnGy3Ukqpr3I0gYhIMPAgMAcYAcwXkRHHXXYDUG6MGQTcB9zrvu8I4GpgJDAb+Lf7+ymllPIBp3sgU4AcY8xeY0wj8AIw77hr5gEL3R+/AswSq/7xPOAFY0yDMWYfkOP+fkoppXzA6X0g6UB+q88LgG/SQCYAAASvSURBVKknusYY0ywilUCS+/aVx903/fgAInITcJP702oR2dmJ9iYDpZ24v7/oLo8D9LH4q+7yWLrL44DOPZbM9lzU7TYSGmMeAR7piu8lImvbs53f33WXxwH6WPxVd3ks3eVxgDOPxekhrEKg9YHiGe7b2rxGREKAeKCsnfdVSinlEKcTyBpgsIhki0gY1qT4ouOuWQQscH98JfCRsSo+LgKudq/SygYGA6sdardSSqnjODqE5Z7TuA1YCgQDjxtjtorIXcBaY8wi4DHgaRHJAQ5jJRnc170EbAOagR8YY1psbnKXDIX5ge7yOEAfi7/qLo+luzwOcOCxdOty7koppeyjO9GVUkp5RROIUkopr2gCOY6IXCUiW0XEJSKTWt2eJSJ1IrLR/faQL9vZHid6LO6vBWxZGBH5vYgUtvpdXOjrNnXEqcr5BBIRyRWRze7fQ8eP//QhEXlcRA6JyJZWt/USkfdFZLf7faIv29heJ3gstv+faAL5ui3A5cCyNr62xxgzzv12i8Pt8kabj6WblIW5r9XvYrGvG9Ne7SznE2jOdv8eAm3/xJNYf/+t3Q58aIwZDHzo/jwQPMnXHwvY/H+iCeQ4xpjtxpjO7F73Gyd5LFoWxnfaU85HOcAYswxrpWdrrUspLQQudbRRXjrBY7GdJpCOyRaRDSLyqYic4evGdEJbJWW+VhbGz90mIl+6u+4BMczg1h1+9q0Z4D0RWecuIxToUowxxe6PDwApvmxMF7D1/6RHJhAR+UBEtrTxdrJXgsVAf2PMeOCnwHMiEneS6x3h5WPxe6d4XP8BBgLjsH4v/+fTxvZsM4wxE7CG5H4gIjN93aCu4t7AHMj7HGz/P+l2tbDawxhzrhf3aQAa3B+vE5E9wBDApxOH3jwWAqAsTHsfl4g8Crxtc3O6kt//7DvCGFPofn9IRF7HGqJra/4wUBwUkVRjTLGIpAKHfN0gbxljDno+tuv/pEf2QLwhIr09E80iMgCrlMpe37bKawFdFsb9j+1xGdZigUDRnnI+AUFEokUk1vMxcD6B9btoS+tSSguAN33Ylk5x4v+kR/ZATkZELgP+BfQG3hGRjcaYC4CZwF0i0gS4gFuMMY5PWnXEiR6Lj8rCdKW/iMg4rOGFXOBm3zan/U5UzsfHzfJWCvC6iID1XPKcMWaJb5vUfiLyPHAWkCwiBcCdwJ+Bl0TkBiAP+IbvWth+J3gsZ9n9f6KlTJRSSnlFh7CUUkp5RROIUkopr2gCUUop5RVNIEoppbyiCUQppZRXNIEopZTyiiYQpZRSXtEEopTDRORmETkgIptEZI+IXOvrNinlDd1IqJTDROQBYIsx5iERmQIsNsYk+7pdSnWU9kCUct4YwHNOyz6g0YdtUcprmkCUct5oYKdYRaRuA37j4/Yo5RUdwlLKQSLSD6vXsQXrIKkvgXON/iOqAKQ9EKWcNRpYZowZh3WezDDgNN82SSnvaAJRylljgA0Axphy4Dlgrk9bpJSXNIEo5azRuBOI21vAhT5qi1KdonMgSimlvKI9EKWUUl7RBKKUUsormkCUUkp5RROIUkopr2gCUUop5RVNIEoppbyiCUQppZRX/j+S0OImOzJIBgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(R, densidad_h(R,z),label='model I')\n",
    "plt.plot(R, densidad_h(R,z,alpha_h=1.63,a_h=1.90,beta_h=2.17),label='model II')\n",
    "plt.xlabel(r'$R$')\n",
    "plt.ylabel(r'$\\rho_{b}/\\rho_{b,0}$')\n",
    "plt.legend(loc='upper right', prop={'size':10})\n",
    "ylim(0.,0.05)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para calcular la velocidad circular volvemos a utiliar la expresión (2.132) del libro Galactic Dynamics, Binney & Tremaine\n",
    "$$\\begin{equation}\n",
    "v_{c}^{2} (R) = 4\\pi G\\sqrt{1-e^{2}} \\int_{0}^{R} dm \\frac{m^{2} \\rho (m^{2})}{\\sqrt{R^{2} - m^{2}e^{2}}}\n",
    "\\end{equation}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{equation}\n",
    "\\frac{v_{c}^{2}}{4\\pi U\\sqrt{1-e^{2}}} = y^{2-\\alpha_{h}}\\int_{0}^{1} \\frac{x^{2-\\alpha_{h}}}{\\sqrt{1-x^{2}e^{2}}} (1-xy)^{\\alpha_{h}-\\beta_{h}}\n",
    "\\end{equation}$$\n",
    "donde $U = G\\rho_{h,0}a_{h}^{2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Constants\n",
    "G = 4.302e-6\n",
    "q_h = 0.8\n",
    "a_h_I = 3.83\n",
    "a_h_II = 1.9\n",
    "rho_h_I = 0.711e9\n",
    "rho_h_II = 0.266e9\n",
    "U_1 = G*rho_h_I*a_h_I**2.\n",
    "U_2 = G*rho_h_II*a_h_II**2.\n",
    "C = 4.*np.pi*q_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def I_h_1(x,y,q_h=0.8,alpha_h=-2.,a_h=a_h_I,beta_h=2.96):\n",
    "    #e = np.sqrt(1.-q_h**2.)\n",
    "    t1 = x**(2.-alpha_h) \n",
    "    t3 = x*y\n",
    "    t2 = (1.+x*y)**(alpha_h-beta_h)\n",
    "    d = np.sqrt(1.-x**2. * (1.-q_h**2.))\n",
    "    n = t1*t2\n",
    "    return (n*y**(2.-alpha_h))/d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def I_h_2(x,y,q_h=0.8,alpha_h=1.63,a_h=a_h_II,beta_h=2.17):\n",
    "    e = np.sqrt(1.-q_h**2.)\n",
    "    n = x**(2.-alpha_h)*(1+x*y)**(alpha_h-beta_h)\n",
    "    d = np.sqrt(1.-x**2. * e**2.)\n",
    "    return y**(2.-alpha_h)*n/d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.linspace(0.,15.,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Vc_ad_b_1(l):\n",
    "    Integral_h_1 = []\n",
    "    for i in range (0,len(l)):\n",
    "        result = integrate.quad(I_h_1,0,1,args=(l[i]))[0]\n",
    "        Integral_h_1.append(result)\n",
    "    return np.array(Integral_h_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Vc_ad_b_2(l):\n",
    "    Integral_h_2 = []\n",
    "    for i in range (0,len(l)):\n",
    "        result = integrate.quad(I_h_2,0,1,args=(l[i]))[0]\n",
    "        Integral_h_2.append(result)\n",
    "    return np.array(Integral_h_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = np.linspace(0.,15.,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEMCAYAAAAF2YvKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmYXGWZ9/HvXdVburN3ZyM72SAJoUNCgCCLoIjK5qgMiAKjDDAO7+jrMiOjI+q4jDMu4zs6jiAIjuyrILLJKpAICYTshNBk6Wyd7qy9Vddyv3+c6qQIWbqqq7qqq3+f6zrXOec5dU7d3UnXXc9ynmPujoiISCZC+Q5ARER6LyURERHJmJKIiIhkTElEREQypiQiIiIZUxIREZGMKYmIiEjGlERERCRjSiIiIpKxknwHkAs1NTU+YcKEfIchItKrLF68uNHdh6VzTlEmkQkTJrBo0aJ8hyEi0quY2fp0z1FzloiIZExJREREMqYkIiIiGVMSERGRjCmJiIhIxpREREQkY0oiIiKSsaK8T0SkL3B3onGnLRonEo0TiSVoT64jsQQdsQTR+P6lI+7E4gliCScWd+KJYDveubjjHlw34ZBIrjncI7TNCBmEzDAgFDLMIGxGyIxQKDheEjLCoRDhEIRDoeS+URoOykvCRmkoRGnYKC0JBdslRmk4RFk4RGk4OFZWEgqWcAgz66lftRyGkohID3MPPvh3tUbZ1Rpld1uU3W0d7GmLsac9yt72GHvbYzRHojRHYjRH4rREYrR2xGntiNESidMejdMWjRNPHOYDPosO9nl9uNzSE8rCoX1JpXzfEqa8NGW7JERFabAuLw1TURrsV5QE2/3KkttlYfolj/crDdMvud+vLExlaQkVZUpch6IkIpIF7s7O1ijb9rTTsDfC9r0RGva209TcQVNzhKaWDhqbO9jZ0sHO1g4iscRhr9e/vCRYKkqoKi+hqizM0KoyqsrC9CsroTLlQy71g7JznfrhWpryTb4k+a2/JBQiHLKgRhA2whbUDCxZq0itWRxJIuE4Qc3F99VggtpNwoPj8eR+LOEkEk40niCeCGpS8YQTTSSIxX1fralzuyOeIBp3OmIJOmLxYDu+v6bVEUvQEY8TiQavjUQTRGJx2qMJWjti7GzdXztrjyaIROO0J6+TrnDI3pVgKsuSSaYsTL/SEqrK373dryxMVfLfqqr8gHVZCZXl+4/35uSkJCLSBW0dcTbtamXjjjbqd7WxeVcbW3a1sXlXO1v2tLFtT4SOgySGyrIw1f3LqK4q56hBFcw8aiBDqsoYUlnGkMpSBleWMrBfKYOSy4CKUvqXlxDuwod3oehMNGF6T8zxhO+rzbV1xInE4rR1JGiL7q/ltUfjtHYExztf17avLKgZdu43Nbfu226NxGiNxrtcUzODytIgwaQmm/7J7f7J8qBs/7H++8revV1R2rM1JiURkaT2aJx3Glt4e3sz6xpbWNfUyrrGFtbvaGX73si7XlsaNkYOqmDUoH7MGTeEEQMr9i3DB5YzfEA5Nf3LqSrXn1ghCods3wdzLrg7kVhiXzNkSzLptEaC7ZZIjJZkwuncTl03R2I07G2nJRKnOfma1o5413+2svC+LyRV5WH6V5QyIKV227+8hAEVKUmoooQBGf4u9D9c+pz2aJy1Dc28uXUva7YFy1sNzWza1faub48jBpYzvrqKM6cOY9zQSsYOrWTs0H6MGVJJTf/yXlVbkJ5lZkHfS2mY6ixdM57wfX1inYmlM+G0dMRobg/6z5ojUVoicfa27z++u7WDTTtb953bHIllKSolESlye9ujLN+0h2WbdrFi8x5WbdnD29tb9nVIl4VDHD2sitnjhvDJOWM5elgVRw+rYmJNFZVl+vOQwhEOGQMqgibP7kokPEg8kSD57EkmnDN+mP619FciRSMWT7BmWzOvbdjJaxt2smTjLuq2t+w7PmpQBdNHDeSc6SM5dtRApo0cwITqSkrCul1K+pZQakIa1L1r9XgSMbNbgPOABnefmSy7G5iWfMlgYJe715rZBGAV8Gby2EJ3v7ZnI5ZCFYnFeWPjbl55p4m/vLOD19bvpCXZblxdVcbscYO5qHY0x40ZxHGjB1HTvzzPEYsUn3zURG4Ffg78trPA3f+6c9vMfgzsTnn92+5e22PRScFKJJwVm/fw4tpGXn67kVfX7aA9GoyImjZiAB87YTQnThjK7LFDGDu0X68eNinSW/R4EnH3F5I1jPew4K/+YuCsnoxJCteu1g5eeKuR51Y38Pya7TS1dAAwdUR/LjlxHPMnVXPihKEMqSrLc6QifVOh9YmcBmxz97dSyiaa2evAHuAb7v7n/IQmPWXjjlaeWrmNJ1du5dV1O4knnCGVpZwxdRhnTBvGqZNrGD6gIt9higiFl0QuBe5M2d8CjHP3JjObAzxkZjPcfc+BJ5rZ1cDVAOPGjeuRYPucWAd0NEO0FaJt0NECsfZgu3Md7wi2Y5FgiXekLNH960QU4rFgnYjRFomwbVcrjXtaaW2PMI0Ec8qM6uowg/uVUFUawvYm4JU4vOLgCZITPQEp+xDsd4mlzOeR3LZQcruzLLS/fN+x5HYo/O5jFgILH3AsdT+cXB+4n1peklKW3A4ltztf17n/rrKSlOVQZSUQKj1gvwTCqWWlKedowIEcWcEkETMrAf4KmNNZ5u4RIJLcXmxmbwNTgUUHnu/uNwI3AsydOzfPs/oUqFgEWhqhtRHadkLrjmDdvgvadkH77mCJ7H330tEcJIxENPP3DpVCuCz4wAqXQqiUeKiElqixO5KgOQpxwgwsK2PskHIGVVVSUVaW8oGb8iHeub/vg9/2b3feNX2k/pDO5JO63ZmE9iWkZJkn2J+oEikJLAHxeHAsEU+Wx/cfTy1LxJPrA/eTr0ndT8ToeiLMJUtJMCUHJJ2UhHPga1KPpx5L/rsTCh9iuwTCJe+9bufxIx4rOfR7hg84pv6yrCmYJAJ8AFjt7vWdBWY2DNjh7nEzOxqYAtTlK8CCFeuAPfWwZ3Ow7K6HvVuheSs0NwTbLY3QsffQ1yipgIpBwVI+EMoHwICRwbqsP5RVQVllsF1amVz67V9KKoKltALC5VCSXMLlwR9u8o82nnCeX9PAvYvqeXpVAx3xBNNGDOCC2qM4b9YoxldX9dAvrcB1JqFEbP/SmWwS0fcmnX3raJCoErHgeDyaXMdSjsdTrpMsj6ceix3i+P6a4/79lCW1hhmLvKe2GRxPfW3yPTrP60kHJpxDJZt95Wm+7qCJ7nDXyPS9S/NeY8zHEN87gTOBGjOrB25w95uBS3h3UxbA6cB3zCwKJIBr3X1HT8ZbMDpaoGltcqkL1rvWw64NQeI48Jtr+UDoPyJIBEfNhv7DoaoGKmuCdb+h0G8IVA6FisHBh38Obd3dzt2vbuTuVzeweXc71VVlfPrk8Xx8zmimjxqokVQHMgs+LMKF9D0vh/bV7KIHJKro/uSzLwlFDzgWfXeSS/vYge9xsP1Y0Ex7sPJDXa+napMW6noN7EjlGcjH6KxLD1F+5UHK7gfuz3VMBSWRgJ3vwJY3YOtSaFgVLLs28K7/lAOOgqETYeIZMHhcsAwaDQNHw4BRUN4/bz9CJ3fntQ07ueXFdTy+YivxhHPalBq+ef50zj52BKW6yU86dfYfhcJAkQyaSMQPSErxriesTBNdl5JnLGi9SLS89xoZ6CNfcwpY206oXwQb/xIsm5dAJDluIFQKNVNg9ByY/WmomQrVk2Ho0UHTUoGKxRM8umwLt7z4Dm/U72ZQv1KuOm0il80bz7jqwo1bJKtCvTApfjH9FgElkZ4WbYMNC6DuuWDZshTw4FvYyOPguE/CUbUw6ngYdiyU9J77HyKxOPctrudXz9exYUcrRw+r4l8vmsnHTxiteahEipT+sntC6w5Y8zisfhTWPg2xtqCWMXYenHk9jD8FjjqhIJqgMtEejXPHXzbwP8+/TcPeCMePGcQ3PjqHDxw7oksPNRKR3ktJJFeibUHSWHJHUOPweNBfMfvTMPVDMO6UXps0OkXjCe5dVM9/PfMWW3a3c9LEofzk4lpOnVytjnKRPkJJJNsaVsErN8Ky+yGyGwaNhVO/ANMvgFG1RTE+3d15fPlWfvj4atY1tTJ73GB+/MnjmT+5Jt+hiUgPUxLJBnd45wV4+b9g7VNQ0g+mXwi1n4IJp+V9HHc2Lavfzb/+YSWvrNvBtBEDuPmKuZx1zHDVPET6KCWR7nrnz/DUN2Hza1A1DN7/DZj7WajK1vPMCkNTc4QfPLaa+xbXU9O/jO9/7Dj++sSxerqfSB+nJJKp7W8GyWPN4zBwDJz/M5h1Sc5v2utpiYRzz6KN/OCx1bR2xLjmjKO57v2Ts/J0NRHp/ZRE0hVth2e/Bwt+EUwFcvYNcPLfBVN/FJm3tu3l+geWsWj9TuZNHMr3LprJlBED8h2WiBQQJZF0bFkKD14DDSvhhMuDBFJVfJ3J8YRz05/r+MmTa6gqD/OjTx7Px08YrX4PEXkPJZGuSCTgpZ/Csz+Aymr41L0w9Zx8R5UTddub+fK9b/D6hl2cO2Mk3/3YTD1WVkQOSUnkSKLt8NC1sOJBmPEx+OhPgkkLi4y7c9erG/n2IysoLwnzs0tqueD4o1T7EJHDUhI5nNYdcNengmlKPvgdmP8PRXGfx4H2tEf55weW8YelW3jf5Bp+cvHxDB9YXAMERCQ3lEQOZec6+N0ngunWP3ELzPx4viPKiaX1u7jujtfZtKuNr35oGn93xiRNVSIiXaYkcjB7t8Gt5wez6V7+exg/P98R5cS9izby9YeWU1NVxt1Xn8zcCcXXTCciuaUkcqCOFrjj4uARslc+CqNPyHdEWReNJ/jeo6u49eV1zJ9Uzc8/dQJDq3rPbMEiUjiURFIl4nDf54KHQV1yZ1EmkJ0tHXz+9tdYUNfE5943kes/fAwlejiUiGRISaSTOzz2T7DmMfjIj2DaufmOKOs2NLVy5W9eoX5XGz/+5PF8fM6YfIckIr1cj38FNbNbzKzBzJanlH3LzDaZ2ZLk8pGUY9eb2Voze9PMPpSzwN64C169CU65Dub9bc7eJl+W1u/ir375EjtaO7jjqpOUQEQkK/LRjnErcLCv+T9199rk8kcAM5sOXALMSJ7z32YWznpELY3wxPUw9mT44L9m/fL59uzqBv76VwupKA1z37Xz1YEuIlnT40nE3V8AdnTx5RcCd7l7xN3fAdYC87Ie1BNfh0hzMIliEU3bDvDIG5u56reLmDS8igc+P5/Jw3v3g7BEpLAU0ifmdWa2NNncNSRZNhrYmPKa+mRZ9rz9DCy9C973RRh+TFYvnW/3L67nC3e9zpzxQ7jr6lMYPkA3EIpIdhVKEvklMAmoBbYAP073AmZ2tZktMrNF27dv79pJ0Tb4w5dg6CQ47SvpvmVBu/OVDXzlvjeYP6mG2/5mHv3LNYZCRLKvIJKIu29z97i7J4Cb2N9ktQkYm/LSMcmyg13jRnef6+5zhw0b1rU3fv7fYec7cP5/FtVzQP534Xquf2AZZ04dxq+vmEu/sux3I4mIQIEkETMblbL7MaBz5NbDwCVmVm5mE4EpwCtZedPmBljwczj+UzDx9KxcshDcv7ief3loOR84dgS/+sxcKkqVQEQkd3q8jcPM7gTOBGrMrB64ATjTzGoBB9YB1wC4+wozuwdYCcSAv3f3eFYCWXwbxDvgtC9l5XKF4PHlW/nqfW9w6uRqfv6p2ZSVFMR3BBEpYj2eRNz90oMU33yY138P+F5Wg4hHYdEtMOksqJmS1UvnywtrtvMPd75O7djB3KgaiIj0kL75VXX1o7B3M8y7Ot+RZMWSjbu45n8XM2l4f35z5Tyq1IkuIj2kbyaRV38Ng8fBlN7/dMKNO1q56rZXqRlQxm8/O49BlaX5DklE+pC+l0S2rYR1f4YTr4JQ727y2dMe5bO3vkokluA3V57IsAF6jK2I9Ky+l0RevQlKKmD2Z/IdSbdE4wn+/vbXeKexhV99eg6Thw/Id0gi0gf1rcbztl3BRIszP9Hrn5N+w8Mr+PNbjfz7x2cxf3JNvsMRkT6qb9VE3rgLoq29fpbeu1/dwB1/2cC1Z0zi4hPHHvkEEZEc6VtJ5K0nYNixcFRtviPJ2LL63fzL71fwvsk1fPVD0/Idjoj0cX0nicRjsPGVXv289J0tHVz7u8XUVJXxs0tqCYcs3yGJSB/Xd/pEti2DjuZem0TiCeeLdy9h+94I91x7CtX9NRJLRPKv79RE1i8I1r00ifzyubU8v2Y73zx/OrVjB+c7HBERoE8lkZdg8HgYeFS+I0nbko27+Omf3uL844/ispPG5TscEZF9+kYScYcNC3tlLaQlEuOLd73OyIEVfPeimZipH0RECkff6BNpfAtaG2HcKfmOJG3ffXQl63e0cuffnsygfprSREQKS9+oiWx4OVj3sprIEyu2cucrG7nm9EmcfHR1vsMREXmPvpFE1i+AqmFQPTnfkXRZY3OE6x9YxszRA/nSB6fmOxwRkYPqG0lkw8sw7mToRf0J33lkJc3tMX56ca0eLiUiBav4P512b4JdG2Bc72nKemb1Nh5+YzPXnTWZKSM0saKIFK7iTyIbOu8P6R2d6s2RGN94cDlTR/Tn2jMm5TscEZHD6vEkYma3mFmDmS1PKfsPM1ttZkvN7EEzG5wsn2BmbWa2JLn8T9pvuP5lKBsAI47L4k+ROz964k227GnnB381S81YIlLw8vEpdStw7gFlTwEz3X0WsAa4PuXY2+5em1yuTfvdNiyAsfMgXPijmRev38ltC9ZxxSkTmDN+SL7DERE5oh5PIu7+ArDjgLIn3T2W3F0IjMnKm7XugIaVvaIpKxZP8PUHlzFqYAVf0ey8ItJLFGJ7yWeBx1L2J5rZ62b2vJmddqiTzOxqM1tkZou2b98eFO6oC9YjZuYu2iy569WNrN66l2+cN53+5YVfaxIRgQJLImb2dSAG3J4s2gKMc/fZwJeAO8xs4MHOdfcb3X2uu88dNmxYUBiLBOvSfrkNvJt2t0b58ZNvctLEoXx45sh8hyMi0mUFk0TM7ErgPOAyd3cAd4+4e1NyezHwNtD1O+/iySQSLstusFn2n0+vYXdblG+eP11zY4lIr1IQScTMzgX+EbjA3VtTyoeZWTi5fTQwBajr8oVjHcE6XLjP3nhr215+u2A9l8wbx4yjBuU7HBGRtPR447uZ3QmcCdSYWT1wA8ForHLgqeQ38YXJkVinA98xsyiQAK519x0HvfDBdNZESgqzJuLufOcPK6ksC/NlTW0iIr1QjycRd7/0IMU3H+K19wP3Z/xmBV4TeW7Ndv78ViP/ct50PalQRHqlgmjOypkCrokkEs5/PP4m44ZWcvkp4/MdjohIRoo7iXSOzirAmsgfl29h5ZY9fOmDUykNF/c/g4gUr+L+9Ionm7NKCiuJxOIJfvLkGqaNGMD5x/e+x/WKiHQq7iQSK8whvg+8tom6xha+fM5UwiEN6RWR3qu4k0gB1kQisTj/+ac1HD92MB+cPiLf4YiIdEtxJ5FYBDAIFc40Inf8ZQObd7fzjx+aphsLRaTXK+4kEo8EtZAC+bBu64jzi2fXMn9SNadOrsl3OCIi3VbkSSRaUCOz7lm0kcbmDr5w9pR8hyIikhVpJxEz+6dcBJITsUjB3CMSjSe48YU65owfwryJQ/MdjohIVhyxs8DM7kndBWqBH+YsomyKdxRMTeSRNzazaVcb37lwhvpCRKRodKXHeY+7X9W5Y2a/zGE82VUgNZFEwvnlc28zbcQA3j9teL7DERHJmq40Z33vgP2v5yKQnIhHCuIekadXN/BWQzN/d+YkQrovRESKyBGTiLu/A2BmNcn9rs+im2+xjrwnEXfnv59by5gh/Thv1qi8xiIikm3pdKzfkrMocqVziG8eLazbwesbdnHN6UdTojmyRKTIpPOp1vvaYWL571i/6c91VFeV8cm5Y/Mah4hILqSTRDxnUeRKPL8d6+ubWnj2zQYuO2kcFaXhvMUhIpIrqonk0P8uWE/YjMtO1vNCRKQ4pZNErs/GG5rZLWbWYGbLU8qGmtlTZvZWcj0kWW5m9v/MbK2ZLTWzE9J6szzWRFo7YtyzaCPnzhzJiIEVeYlBRCTXupxE3H35kV/VJbcC5x5Q9jXgaXefAjyd3Af4MDAluVwNpHePSiySt5rIg69vYk97jCvnT8jL+4uI9IS0hguZ2QtmNjC5fa2ZfdHM0vqq7+4vAAcOE74QuC25fRtwUUr5bz2wEBhsZl0fJxvvyEtNxN257eV1TB81kDnjh/T4+4uI9JR0x5wOcvc9ZjYH+FtgCHBTFuIY4e5bkttbgc4HbYwGNqa8rj5Z1jV5qoksqGtizbZmrpw/QVOciEhRS/dBG1EzKwEuB37o7veY2aJsBuTubmZpjwQzs6sJmrwYN25cUBiP5uU+kd++vJ4hlaVcUKtH34pIcUu3JvJfwBvAecAjybL+WYhjW2czVXLdkCzfBKTeYDEmWfYe7n6ju89197nDhg0LCvMw7cnmXW08uXIrf32ihvWKSPFLK4m4+23AKcBMd28zs8nAgizE8TBwRXL7CuD3KeWXJ0dpnQzsTmn2OlKwyT6Rnq2J3Le4noTDZSeN69H3FRHJh7Sas8zsJuCjQNzMNgNLgdfSvMadwJlAjZnVAzcA/wbcY2afA9YDFydf/kfgI8BaoBX4my6/Uefz1XuwJpJIOPcs2sipk6sZO7Syx95XRCRf0u0TOR0Y6+5xMxsNHA/MSucC7n7pIQ6dfZDXOvD3acYYiEWCdQ/WRBbWNVG/s42vfmhaj72niEg+pZtE/gJUAw3uvomgf+KPWY8qG/JQE7ln0UYGVJTwoRkje+w9RUTyKd2O9V8Bz5vZV8zsNDMblIugsqKzJtJDSWR3W5THlm/lotrR6lAXkT4j3STyO+C3BDWYzwMvm9nbWY8qG+I925z1yBubicQSXKzZekWkD0m3Oave3X+QWmBmhfEQ8wPFerY5695FGzlm5ABmjh7YI+8nIlII0q2JLDGzL6QWuHski/FkTw/WRFZv3cMb9bu5eO5Y3aEuIn1KujWREcAHzOyfCIb2vgEscfd7sx5Zd+2rieQ+idy7qJ7SsHHR7K7PyCIiUgzSSiLufjHsa8KaARwHnAQUXhLZVxPJbXNWLJ7g90s28YFjRzC0Kr/PcxcR6Wnp3mz4XXf/RrIJ6zXgNTMrzKFI+0Zn5bYmsqCuicbmDi6sVS1ERPqedPtERpvZpzp3zGwY8KfshpQlnfeJ5Lgm8vCSzQwoL+HMacNy+j4iIoUo3T6Ra4AnzGwtwTPXfwP8U9ajyoYeqIlEYnEeX7GVc2aM1L0hItIndSmJmNlvCZqvXieYhuQOIAZc5O5rcxdeN8SjwTqHo7Oef3M7e9tjnH9815+TJSJSTLranHUrYAQTIP4OmADsBD5tZp/ISWTdFc/9HeuPLN3C0KoyTp1ck7P3EBEpZF2qibj7M8AznfvJB1MdSzAB40nAfTmJrjtyPAFja0eMP63cxsfnjKY0nG7XkohIcehqc9YpwMLkrLq4ewxYllx+l7vwuiHHEzA+tXIbbdE458/S0wtFpO/q6lfoywmG895lZleaWeFPU5vjmsgjb2xh5MAKTpwwNCfXFxHpDbranPV3AGZ2DPBh4NbkDL7PAo8DL7l7PGdRZiKHNZHdrVGeX9PAFadMIBTSNCci0nel+3jc1e7+U3c/FzgLeBH4JMFzRgpLLAIYhNIdxXxkT6zYSjTuXFCrpiwR6dsy/oRNPmP9aHf/PwBmNtjdd2UvtG6KR4KmrBxMiPjEiq2MGdKP40YX7uNURER6QneHFY1P2f7n7lzIzKaZ2ZKUZY+ZfdHMvmVmm1LKP9KlC8Y6cnKjYUskxp/XNnLO9JGasVdE+rwuJREzO/ZQ5yefcBgCutXD7O5vunutu9cCc4BW4MHk4Z92HnP3rj2ONx7JyZQnz6/ZTkcswYdmjMj6tUVEepuu1kQeNbPfmNm4A8q/CswCbgJ+n8W4zgbedvf1GV8hRzWRJ1dsZWhVGXPGD8n6tUVEepuuJpFjCKY9ed7MfpaceBF3T7j7L9z9c+7+SBbjugS4M2X/OjNbama3mFnXPr1zUBOJxhM8vbqBs48ZToluMBQR6VoScfcOd/8vgrvUNwKvmNm/mlnWnwVrZmXABex/RskvgUlALbAF+PEhzrvazBaZ2aLt27cHo7OyXBNZWNfE3vYY58wo/NtkRER6QrpDfNvd/UfATKANWGxmX8lyTB8GXnP3bcn33ObucXdPEDSbzTtEbDe6+1x3nzts2LDgPpEs10SeXLGNfqVhTpuiubJERCDNJGJmE8zsXOAqYBywF/h+lmO6lJSmLDNLnSL3Y8DyLl0lnt0+kUTCeWrlNk6fWqNp30VEkro6d9ZSYDSwAVgNrAKeBn4OrMlWMGZWBXyQ4Lklnf7dzGoJnl+y7oBjhxbryOqUJ8s27Wbrnnb+cca0rF1TRKS36+rNhhcB73ROwJgr7t4CVB9Q9pmMLhaPQOngbIQFBDcYhkPGWccMz9o1RUR6u67OnVWX60CyLss1kSdXbuOkiUMZXJnbx+2KiPQmxTtONR7J2uSL6xpbWNvQzAen6wZDEZFUxZtEYpGs1USee7MBQE1ZIiIHKN4kEu+AcGlWLvXcmu1MrKlifHVVVq4nIlIsijeJZOlmw/ZonAVvN3HG1GFZCEpEpLgUbxKJZ6djfUFdE5FYgjOnKYmIiByoeJNILDsd68+/uZ3ykhAnH1195BeLiPQxxZtEEtGs1ESee7OBUyZV6y51EZGDKM4k0nlPZDdrIusaW1jX1MqZ6g8RETmo4kwiJIJVN2sinUN7z5ymob0iIgdTnElkX02km0kkObR3Qo2G9oqIHEyRJpHOmkjmzVka2isicmRFmkS6XxNZqKG9IiJHVJxJhGQS6UZN5DkN7RUROaLiTCJZqIm8uLaRk47W0F4RkcMp0iTSvdFZW3e3s7ahmdMm6zG4IiKHU5xJhO7dJ/LS2kYA5k9WU5aIyOEUZxLpbM7KsCby0tuNDK0q49iRA7MYlIhI8enq43F7hJnwYmgRAAAO8UlEQVStA/YCcSDm7nPNbChwNzCB4BnrF7v7zsNeqLM5K4OaiLvz0tpG5k+qJhSytM8XEelLCrEm8n53r3X3ucn9rwFPu/sU4Onk/uF1Y9qTt7c3s21PhFPVHyIickSFmEQOdCFwW3L7NuCiI5+SeXPWS2ubAHifkoiIyBEVWhJx4EkzW2xmVyfLRrj7luT2VuDIDzrvRnPWi2sbGTu0H2OHVqZ9rohIX1NQfSLA+9x9k5kNB54ys9WpB93dzcwPdmIy6VwNMG1MclRVmjWRWDzBwromzps1KoPQRUT6noKqibj7puS6AXgQmAdsM7NRAMl1wyHOvdHd57r73P79kxMmplkTWbZpN3vbY+oPERHpooJJImZWZWYDOreBc4DlwMPAFcmXXQH8/ogXy/Bmw877Q07RVCciIl1SSM1ZI4AHzQyCuO5w98fN7FXgHjP7HLAeuPiIV8pw2pOX1jYxfdRAqvt3/4mIIiJ9QcEkEXevA44/SHkTcHaaVwMLQbjrP15bR5zF63dyxfzx6b2ViEgfVjDNWVnlibRrIa+u20FHPKH+EBGRNBRpEvG0p4FfWNdEScg4ccLQHAUlIlJ8ijOJ4GnXRBbWNTFrzCCqygumhU9EpOAVZxJxT2tkVkskxtL63XoAlYhImoo0iSTSukdk0fqdxBLOKZOURERE0lGkSSS9mkhnf8ic8UNyGJSISPEpziSCp1UTWVjXxPFjB1NZpv4QEZF0FGcSSaM5q3lff4hGZYmIpKtIk0jXm7MWrdtBPOHqVBcRyUCRJpGu10QW1u2gNKz+EBGRTBRnEqHrNZGFdU0cP0b9ISIimSjOJOJd61hvjsRYtkn3h4iIZKpIk0iiSzWRV9UfIiLSLcWZRLo4xHdhXROlYeOE8YN7ICYRkeJTnEmki6OzFtbtUH+IiEg3FG8SOcIEjHvboyxXf4iISLcUZxIhccSp4Bet20lc82WJiHRLcSaRLtREFtQ1URYOccI43R8iIpKpgkkiZjbWzJ41s5VmtsLMvpAs/5aZbTKzJcnlI1264BFqIgvrmqgdO5h+ZeHuBy8i0kcVUo9yDPiyu79mZgOAxWb2VPLYT939R2ld7TA1kT3J/pDrzpqScbAiIlJAScTdtwBbktt7zWwVMDrjCx5mdNar7+wg4XCKOtVFikY0GqW+vp729vZ8h1LwKioqGDNmDKWlpd2+VsEkkVRmNgGYDfwFOBW4zswuBxYR1FZ2HuScq4GrAeaMCh32PpEFbzdRVhJi9jjdHyJSLOrr6xkwYAATJkzAzPIdTsFyd5qamqivr2fixIndvl7B9Il0MrP+wP3AF919D/BLYBJQS1BT+fHBznP3G919rrvPBQ6bRBa+08QJ4wZTUar+EJFi0d7eTnV1tRLIEZgZ1dXVWauxFVQSMbNSggRyu7s/AODu29w97u4J4CZgXpcudoiO9d2tUVZs3qP7Q0SKkBJI12Tz91QwScSCn+pmYJW7/ySlfFTKyz4GLO/SBQ/Rsf6Xd5pw9YeISIGbMGECjY2NGb2mK+dmSyH1iZwKfAZYZmZLkmX/DFxqZrWAA+uAa7p0tUN0rC+s20F5SYha9YeIiHRbwdRE3P1Fdzd3n+Xutcnlj+7+GXc/Lll+QXIU15Edok9kQV0Tc8YPobxE/SEikj3r1q3jmGOO4corr2Tq1Klcdtll/OlPf+LUU09lypQpvPLKKwDs2LGDiy66iFmzZnHyySezdOlSAJqamjjnnHOYMWMGV111Fe6+79q/+93vmDdvHrW1tVxzzTXE4/G8/IwHU0g1kew6SE1kV2sHq7fu4UsfmJqHgESkp3z7kRWs3Lwnq9ecftRAbjh/xmFfs3btWu69915uueUWTjzxRO644w5efPFFHn74Yb7//e/z0EMPccMNNzB79mweeughnnnmGS6//HKWLFnCt7/9bd73vvfxzW9+k0cffZSbb74ZgFWrVnH33Xfz0ksvUVpayuc//3luv/12Lr/88qz+fJkq3iRykJrIwrodQX+I5ssSkRyYOHEixx13HAAzZszg7LPPxsw47rjjWLduHQAvvvgi999/PwBnnXUWTU1N7NmzhxdeeIEHHngAgI9+9KMMGRJMyfT000+zePFiTjzxRADa2toYPnx4D/9kh1a8SeQgNZEFbzfSrzTMrDHqDxEpZkeqMeRKefn+z51QKLRvPxQKEYvFMrqmu3PFFVfwgx/8ICsxZlvB9Ilk3QGjs9ydZ9/czvxJ1ZSVFO+PLSKF7bTTTuP2228H4LnnnqOmpoaBAwdy+umnc8cddwDw2GOPsXNncE/12WefzX333UdDQwMQ9KmsX78+P8EfRBHXRN7dnPX29hY27Gjl6tOPzlNAIiLwrW99i89+9rPMmjWLyspKbrvtNgBuuOEGLr30UmbMmMH8+fMZN24cANOnT+e73/0u55xzDolEgtLSUn7xi18wfvz4fP4Y+1jqCIBiMfeosC96czMMGLGv7KYX6vjeH1fx0tfOYvTgfnmMTkRyYdWqVRx77LH5DqPXONjvy8wW75v1o4uKt13ngJrI06u3cczIAUogIiJZVLxJJKVPZE97lEXrdvL+YwpnRIOISDEo3iSSMjrrz2saiSWcs5RERESyqkiTiEFo/x3pz6xuYFC/UmaP1dBeEZFsKs4kkjJDZSLhPL+mgTOmDqMkXJw/rohIvhTnp2pKElm6aTeNzR1qyhIRyYHiTCIpP9YzqxsIGZwxdVge4xERSU+2poLv379/TuLrVJxJJKUm8uzqBmaPG8KQqkM/6VBERDJT1Elk8642lm3araYsEck5TQVfTCzIjT98fDVl4RDnzzoqzwGJSI967GuwdVl2rznyOPjwvx32JZoKvmgYL69t5PdLNvMPZ01mXHVlvgMSkT5AU8EXKDM7F/gZEAZ+7e6H/zpgxjd+v5xxQyv5/Psn90SIIlJIjlBjyBVNBV+AzCwM/AL4MDCd4Jnr0w93Tkcc6ra38O0LZ1BRqsfgikjh0FTwPW8esNbd6wDM7C7gQmDloU7oiDsfnjmS908rnCqfiAhoKvgeZ2afAM5196uS+58BTnL36w51zrSjBvpzq7YxapBm7BXpKzQVfHo0FfwBzOxqM1tkZotiVqYEIiLSA3pDEtkEjE3ZH5Msexd3v9Hd57r73CGjJvRUbCIifVpvSCKvAlPMbKKZlQGXAA/nOSYREaEXdKy7e8zMrgOeIBjie4u7r8hzWCJSgNwdS5n2SA4um33hBZ9EANz9j8Af8x2HiBSuiooKmpqaqK6uViI5DHenqamJioqKrFyvVyQREZEjGTNmDPX19Wzfvj3foRS8iooKxowZk5VrKYmISFEoLS1l4sSJ+Q6jz+kNHesiIlKglERERCRjSiIiIpKxgp/2JBNmthd4M99xdEMNcPjnYha23hx/b44dFH++9fb4p7n7gHROKNaO9TfTnf+lkJjZIsWfH705dlD8+VYM8ad7jpqzREQkY0oiIiKSsWJNIjfmO4BuUvz505tjB8Wfb30u/qLsWBcRkZ5RrDURERHpAUWVRMzsXDN708zWmtnX8h1POsxsrJk9a2YrzWyFmX0h3zFlwszCZva6mf0h37Gky8wGm9l9ZrbazFaZ2Sn5jikdZvZ/k/93lpvZnWaWnRn2csTMbjGzBjNbnlI21MyeMrO3kush+YzxcA4R/38k//8sNbMHzWxwPmM8nIPFn3Lsy2bmZlZzpOsUTRIxszDwC+DDwHTgUjObnt+o0hIDvuzu04GTgb/vZfF3+gKwKt9BZOhnwOPufgxwPL3o5zCz0cA/AHPdfSbBYxMuyW9UR3QrcO4BZV8Dnnb3KcDTyf1CdSvvjf8pYKa7zwLWANf3dFBpuJX3xo+ZjQXOATZ05SJFk0SAecBad69z9w7gLuDCPMfUZe6+xd1fS27vJfgAG53fqNJjZmOAjwK/zncs6TKzQcDpwM0A7t7h7rvyG1XaSoB+ZlYCVAKb8xzPYbn7C8COA4ovBG5Lbt8GXNSjQaXhYPG7+5PuHkvuLiR4EmtBOsTvH+CnwD8CXeowL6YkMhrYmLJfTy/7EO5kZhOA2cBf8htJ2v6T4D9fIt+BZGAisB34TbI57tdmVpXvoLrK3TcBPyL49rgF2O3uT+Y3qoyMcPctye2twIh8BtNNnwUey3cQ6TCzC4FN7v5GV88ppiRSFMysP3A/8EV335PveLrKzM4DGtx9cb5jyVAJcALwS3efDbRQ2E0p75LsO7iQIBkeBVSZ2afzG1X3eDB0tFcOHzWzrxM0Ud+e71i6yswqgX8GvpnOecWURDYBY1P2xyTLeg0zKyVIILe7+wP5jidNpwIXmNk6gqbEs8zsd/kNKS31QL27d9b+7iNIKr3FB4B33H27u0eBB4D5eY4pE9vMbBRAct2Q53jSZmZXAucBl3nvuodiEsGXkDeSf8djgNfMbOThTiqmJPIqMMXMJppZGUGn4sN5jqnLLHie583AKnf/Sb7jSZe7X+/uY9x9AsHv/hl37zXfhN19K7DRzKYli84GVuYxpHRtAE42s8rk/6Wz6UUDA1I8DFyR3L4C+H0eY0mbmZ1L0KR7gbu35juedLj7Mncf7u4Tkn/H9cAJyb+NQyqaJJLszLoOeILgj+ced1+R36jScirwGYJv8EuSy0fyHVQf83+A281sKVALfD/P8XRZsgZ1H/AasIzgb7ug7542szuBBcA0M6s3s88B/wZ80MzeIqhd/Vs+YzycQ8T/c2AA8FTyb/h/8hrkYRwi/vSv07tqWyIiUkiKpiYiIiI9T0lEREQypiQiIiIZUxIREZGMKYmIiEjGlERERCRjSiIiIpIxJRGRHDOzmWb2csr+CWb2dD5jEskW3WwokmNmFiKYln20u8fN7DngS51T/4v0ZiX5DkCk2Ll7wsxWADPMbAqwXglEioWSiEjPWEgwP9rnOcjT5ER6KyURkZ6xkOBxpL9IPkBKpCioT0SkBySbsZ4Hprh7S77jEckWjc4S6RlfAK5XApFioyQikkNmNsnMVgP93P22fMcjkm1qzhIRkYypJiIiIhlTEhERkYwpiYiISMaUREREJGNKIiIikjElERERyZiSiIiIZExJREREMvb/AXoBAWQPATJYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Vc_DM_I = np.sqrt(C*U_1*Vc_ad_b_1(y))\n",
    "Vc_DM_II = np.sqrt(C*U_2*Vc_ad_b_2(y))\n",
    "plt.plot(r,Vc_DM_I,label='model I')\n",
    "plt.plot(r,Vc_DM_II,label='model II')\n",
    "xlim(0.,14.)\n",
    "plt.xlabel(r'$y$')\n",
    "plt.ylabel(r'$V_{c}/km s^{-1}$')\n",
    "plt.legend(loc='lower right', prop={'size':10})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integral_R(x,R,q_h=0.8,alpha_h=-2.,a_h=a_h_I,beta_h=2.96):\n",
    "    e = np.sqrt(1.- q_h**2.)\n",
    "    n = x**(2.-alpha_h)*(1.+ (x*R)/(a_h))**(alpha_h-beta_h)\n",
    "    d = np.sqrt(1-x**2.*e**2)\n",
    "    t = R**2.*(R/a_h)**(-alpha_h)\n",
    "    result = (t*n)/d\n",
    "    return result\n",
    "def integral_R_II(x,R,q_h=0.8,alpha_h=1.63,a_h=a_h_II,beta_h=2.17):\n",
    "    e = np.sqrt(1.- q_h**2.)\n",
    "    n = x**(2.-alpha_h)*(1.+ (x*R)/(a_h))**(alpha_h-beta_h)\n",
    "    d = np.sqrt(1-x**2.*e**2)\n",
    "    t = R**2.*(R/a_h)**(-alpha_h)\n",
    "    result = (t*n)/d\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Vc_R(l):\n",
    "    Integral_R = []\n",
    "    for i in range (0,len(l)):\n",
    "        result = integrate.quad(integral_R,0,0.9,args=(l[i]))[0]\n",
    "        Integral_R.append(result)\n",
    "    return np.array(Integral_R)\n",
    "def Vc_R_II(l):\n",
    "    Integral_R = []\n",
    "    for i in range (0,len(l)):\n",
    "        result = integrate.quad(integral_R_II,0,0.9,args=(l[i]))[0]\n",
    "        Integral_R.append(result)\n",
    "    return np.array(Integral_R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "R= np.linspace(0.1,15.,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  if sys.path[0] == '':\n",
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:12: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  if sys.path[0] == '':\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEMCAYAAAAF2YvKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VfWd//HXJySQQDYgECALAWTfNSqKWhUXXFrs5tRacWmrM9Wp/Y3T1tqOtmMXO7XbzFgrrVYdodaqVapoVVxwR0BkB4EEyEJWspE9+f7+OBeICJh7ucm59+b9fDzu4yz33JtPAsn7fs/3nO/XnHOIiIiEIs7vAkREJHopREREJGQKERERCZlCREREQqYQERGRkClEREQkZAoREREJmUJERERCphAREZGQxftdQE/IyMhweXl5fpchIhJVVq9eXemcGxbMa2IyRPLy8li1apXfZYiIRBUz2xXsa3Q6S0REQqYQERGRkClEREQkZAoREREJmUJERERCphAREZGQKURERCRkMXmfiIh0k3PgOqGzvcujI/BoB9d1vdNbP7DPdQbWOw+tHzzmwLbrsn2kh/v4Ptyh/bgux7mPbuO6v+TQ4tC2O2z9Iwcd+WcVDLPDd3zycwf3W5fD7bD9doTX2GHLLsda3GH7j/IaC61NoRAR6QnOQXtz4NECbU3esr0ZOloP7W9vgY4WaG89bNkCHW3esR2tXda77Ots77KvDTq7Lg8EQpv3R7yj7bCQ6LItQTj8j//RBBk4UUwhIn3TgT/yLQ3QWg+t+wOPhi7r+6Gtscuy0Vu2NXqh0NYUWG/2lu0HloGwCIe4BOjXH/odWPaHfvGH1uPiDz0XnwgDUgKvifeWB56P6xfY7tfl+a6PfmD9Prr9kX0H1uMCy8C2xQXW4w5b7/qcHdr+yHrcoU/AXdexjz/3sU/NR/r0/UlLOPIn/sP39bDDWzSHt4I+1kL6hJbTEVtdXfcRaNUdvv8Iyx/lBP3tKEQkOnW0QVMNNNccWjbXdlnWQnMdtNRBS31gvT7wqPPCIphP4QkDISEJEgZ5y/4DvX2J6ZCS5O2LTzzKcoC3fviy3wCI73/YckAgEALL3vrDJr3n8H/TKP83VoiI/zo7oWkf7K849GisgsZqaKwMLKugqRoa93nHttYf+z379YfENBiQComp3if0QWO85YAU6J/80fX+g2BAsreeMNDb7p/shUV8kvcJXEQ+RiEiPaejDer3Qn0p1JVAQ5m33VDW5VEO+yu9TtgjSUyDgUMhaQgkZ8KwyTBwCCQN9h6J6ZCU7h13YH1AKiQk9u73KtJHKUQkNM55LYKaXVCzB2p2Q20R1BVBbTHUFXsBcXgHY1y8FwbJwyE1C0bNhkHDve1BGTAw49By4BDvfL6IRCyFiBxdZyfUl0DVDqjeCdU7YF9h4LHL61voqn+yFwxpWZA5BVKzIXUkpIwKLEd6LQqdGhKJGQoR8U47VW2H8s1QsRUqt0HVh1C5HdqbDh3XbwAMzvMeuafD4NGQnnvokZge9Z2EIhIchUhf01ABez+AvRugbAOUbfRC4+CVSuaFw9DxkHcWDB3nPYaM81oZakWISBcKkVjWVAPFq6F4DZS8D6Vrvb6KA1KzIXMqTLjQ67AePgkyJniXpYqIdEOvh4iZ5QAPA5l4va6LnHO/NbMfAl8HKgKH3uacWxZ4zfeArwIdwDedc//o7bojnnOwrwB2vQ2734I9K70WBgAGQ0+A0afDyFkwcqYXHgOH+FqyiEQ/P1oi7cAtzrk1ZpYCrDazFwPP/do5d3fXg81sCvAlYCowCnjJzCY4d7RrQvuQfbug4DUoWAEFr0PDXm9/YjrknAozLoesfMg60bsEVkQkzHo9RJxzpUBpYL3ezDYDWcd4yQLgUedcC1BgZtuBU4C3e7zYSNPWBIVvwIcvwIcvei0P8C6RHXOW19IYfTpkTFTfhYj0Cl/7RMwsD5gNvAvMBW4ys4XAKrzWyj68gHmny8uKOHboxJbGatj6HGz+O+x8xRuTKT7JC41T/xnGfgqGTdJVUSLiC99CxMySgSeAbznn6szsXuBOvH6SO4FfAtcF8X7XA9cD5Obmhr/g3tS0DzYthQ2PQ+Gb3t3caTlw4tUw4QIYfYbuyBaRiOBLiJhZAl6ALHbOPQngnCvr8vwfgGcCm8VA16ElswP7PsI5twhYBJCfnx994zC3t8K252DdY97pqo5WrzP8jG/B5E97HeJqbYhIhPHj6iwD7gc2O+d+1WX/yEB/CcBngQ2B9aXAEjP7FV7H+nhgZS+W3LOqdsCah2HtYm/gweQRcPLXYcYXFRwiEvH8aInMBa4C1pvZ2sC+24ArzGwW3umsQuAGAOfcRjN7DNiEd2XXjVF/ZZZzXv/G2/fA9pe8+RUmzIeTroYTzvPmYBARiQJ+XJ31BkeeHmzZMV7zE+AnPVZUb+lo805XvX0PlG/0rqo6+zY4caE3tpSISJTRHeu9oaMd1j0KK37hDV44fCos+B1M/4I3+ZCISJRSiPQk52D94/DqT71RcEfOhCv+4g0zor4OEYkBCpGeUrwanrsVilZC5nT40p9h4kUKDxGJKQqRcNtfCS/e7l1tNWg4LLgHZn5Zd5CLSExSiITTxqfg2VuguRbm3gxn/rs3v7eISIxSiITD/kpY9u+w8W/edK8L/u7N7CciEuMUIser4HV4/DporoF5t8PpN0M//VhFpG/QX7tQOQdv/Q+89ENv5r+FT3lzdIiI9CEKkVC01MNT34DNS2HKAq/zfECK31WJiPQ6hUiwGsrhkc9B2Sa44Mdw2k26bFdE+iyFSDBqdsPDC6B+L3z5MRh/nt8ViYj4SiHSXRVb4eHLoG0/XPUU5J7qd0UiIr5TiHRH+Wb408UQFw/XLIMR0/yuSEQkIihEPkltMTzyeejXH65d5l2JJSIigELk2JprYfEXoblOASIicgQKkaNpb4FHr4TKrXDl4zByht8ViYhEHIXIkTgHS78Jha/DZxfBuHP8rkhEJCJpaNkj+eBRbxKps78HM//J72pERCKWQuRw1QWw7NuQezqc9W2/qxERiWgKka462uHJ68Hi4HP3QVw/vysSEYlo6hPp6vVfejMRfv5+SM/1uxoRkYinlsgBRavhtZ/D9Mth+hf8rkZEJCooRMC7Guv578KgYXDJ3X5XIyISNRQiAFuegaL34JzbIDHN72pERKKGQqSjHV76EWRMhFlX+l2NiEhUUcf6+/8HVR/Cl5ZoWlsRkSD17ZZI63549WeQMwcmXux3NSIiUadvf/R+53fQUAaX/59mJxQRCUGvt0TMLMfMXjGzTWa20cxuDuwfYmYvmtmHgeXgwH4zs/82s+1mts7MTgxLIU018MZvYdKlmmBKRCREfpzOagducc5NAeYAN5rZFOBWYLlzbjywPLANcBEwPvC4Hrg3LFVseAJa6+HMW8LydiIifVGvh4hzrtQ5tyawXg9sBrKABcBDgcMeAi4LrC8AHnaed4B0Mxt53IW8/whkToNRs4/7rURE+ipfO9bNLA+YDbwLZDrnSgNP7QUyA+tZwJ4uLysK7Dv8va43s1VmtqqiouLYX7hsI5SsgdlfUV+IiMhx8C1EzCwZeAL4lnOurutzzjkHuGDezzm3yDmX75zLHzZs2LEPfn8xxCV4Q5yIiEjIfAkRM0vAC5DFzrknA7vLDpymCizLA/uLgZwuL88O7AtNe6s3V8iki2HQ0JDfRkRE/Lk6y4D7gc3OuV91eWopcHVg/Wrg6S77Fwau0poD1HY57RW8bc9BYxXMXhjyW4iIiMeP+0TmAlcB681sbWDfbcBdwGNm9lVgF3DgXNMy4GJgO9AIXHtcX/39RyBllKa8FREJg14PEefcG8DRerPnHeF4B9wYli9eVwLbX4Iz/k0TTomIhEHfGvZk7RJwnTBbAy2KiIRD3wqRrcsg+xQYMtbvSkREYkLfCZHWRij9APLO8LsSEZGY0XdCpHg1dLZD7ml+VyIiEjP6Tojsfsdb5pzsbx0iIjGk74TInndg2GRIGux3JSIiMaNvhEhnB+xZCblz/K5ERCSm9I0QKd8MLXUKERGRMOsbIbIn0B+iEBERCau+ESK734HkEZA+2u9KRERiSh8JkXe9KXA1d4iISFjFfojUFkPtbt0fIiLSA2I/RA70h+Sc6m8dIiIxKPZDZPe7kDAIRszwuxIRkZjTB0Lkbcg+Cfr5MXWKiEhsi+0QaamHsg2Qo0t7RUR6QmyHSNF73vwhuj9ERKRHxHaIJKbB9MshW4Muioj0hNjuKMg6CT7/B7+rEBGJWbHdEhERkR6lEBERkZApREREJGRBh4iZfbcnChERkejziR3rZvZY101gFvDzHqtIRESiRneuzqpzzn3twIaZ3duD9YiISBTpzumsnxy2/f2eKERERKLPJ4aIc64AwMwyAtvVPV2UiIhEh2A61h/osSpERCQqBRMiYZkW0MweMLNyM9vQZd8PzazYzNYGHhd3ee57ZrbdzLaa2YXhqEFERMIjmBBxYfqaDwLzj7D/1865WYHHMgAzmwJ8CZgaeM3vzKxfmOoQEZHj1OstEefcCqC7/SoLgEedcy2BvpntwCnhqENERI5fMAMwfq/HqvDcZGYLgVXALc65fUAW8E6XY4oC+z7GzK4HrgfIzc3t4VJFRKJXfXMbJTXNlNQ2UVrTzN7aJkprm0N6r26HiHNuwycfFbJ7gTvxTpndCfwSuC6YN3DOLQIWAeTn54fr1JuISFRp7+iktLaZkpomSmqbKKlpprimiZIaLzBKapqob2n/yGviDIanJIb09YIaCt7MVgCXOufqzOyfgUTgd8651pC+eoBzrqzL1/gD8ExgsxjI6XJodmCfiEif1NLeQWlNM0X7mija10jRviaKa7z14n1N7K1rpvOwj9FDBvVnVHoio4cO5LRxQxmZlsio9CRGpScyMi2J4SkDiO8Xh4VwF2Cw84mkBQLkJODreH/s/wBcHfyXPsTMRjrnSgObnwUOtHqWAkvM7FfAKGA8sPJ4vpaISCTr6HTsrWtmd1Uje/Y1UlTdyJ59Teyp9gKjrL4Z1yUk+sUZI1ITyRqcxJyxQ8kanERWehKj0pPIGpzEqLQkkvr33PVIwYZIm5nFAwuBnzvnHjOzVcG8gZn9GTgbyDCzIuAO4Gwzm4V3OqsQuAHAObcxMHbXJqAduNE51xFkzSIiEWV/Szu7qxvZVdXI7ur97K5uZHf1gaBopK3jUErEGYxMSyJ7cBJzT8ggZ0gS2YMHkj3Y2zciNZH4fv4NyG7Odb/7wMyuBr6DdxprmnOuycy2OOcm9VSBocjPz3erVgWVbSIiYVXX3EZh5X4KKvezq6ox8NhPYVUjlQ0tHzk2LSmB3CEDyR0ykJzA0lv3WhQJvRQSZrbaOZcfzGuCaok45x4ys78BbYEAOQF4O5j3EBGJFY2t7RRWNlJQuZ+CygYKKhsprPKCo3r/R7uKR6R6fRLzJg0nd+hARg/1gmL0kEGkDUzw6Ts4fsF2rP8BuAToMLMSYB2wpicKExGJBJ2djpLaJnZU7GdHeQMFlfvZWdnAzor9H7ssNjN1AGMyBnHh1Ezyhg5i9NBBjMkYRO6QgT3aL+GnYPtEzgJynHMdZpYFzARmhL8sEZHe1dLeQUHlfraXN7C9vOFgaOysbKC5rfPgcSmJ8YwdlsycsUMZmzGIMcO8oMgbOohBA4L9kxr9gv2O3wWGAuXOuWK8y22Xhb0qEZEe0tTawY6KBj4sr2dbWcPB0NhVtf/gpbFmkD04iXHDkjlt3FDGDUtm7LBBjBuWTEZyf8zCMoBHTAg2RO4DXjOz+/ECZZ1zrjb8ZYmIHJ/W9k52VjawdW8928q8wNhWVs/u6saDl8jGxxl5GYOYPDKFT88YybjhyYwfnsLYYYNITIjN00/hFmyIPIJ3V3g88A1ghpklOufGhb0yEZFucM5RXNPEltJ6tpbVs2VvPVv31rGzYj/tgaZFfJwxJmMQ00al8dnZWUzITGH88GTyMgb12pVPsSrYEClyzv2s6w4zGxDGekREjqqptYOtZfVsLq07+NhSWv+RYTyyBycxaUQK50/JZEJmChNHpDA2I5n+8QqLnhBsiKw1s5udc789sMM513KsF4iIhKKqoYWNJXVsKq3zliW1FFQe6rdIHhDPpBEpXDY7i4kjUpg8MoUJmSmkJEbv5bLRKNgQyQTOM7Pv4l3a+wGw1jn317BXJiJ9gnPeMB/ri2rZEAiLDcV17K07dPlsVnoSU0alcumMUUwemcqUkankDElSB3cECPZmw8vh4CmsqcB04FRAISIin8g5R1ldC+uKalhfXMu6olo2ltRS2eDdmBdnMG5YMnPGDmHqqDSmZnmBkT6wv8+Vy9EEe7Phj51zPwicwloDrNFMgyJyNPv2t/JBUQ3rimpZV1TDB0W1VNR7Z8DjDMYPT+HsicOZNiqV6dlpTB6ZysD+fe9ei2gW7L9Wlpl92Tm3BMDMhgGPAeeEvTIRiSot7R1sLKlj7e4aPiiqYe2eGnZVNR58ftywQZxxQgYzstOYkZ3GlJFpMXsXd18SbIjcAPzDzLbjjbj7J+C7Ya9KRCKac46S2mbW7NrHmt37eH93DZtK6mjt8O7sHpGayMycNL50ci4zs9OYnp2mDu8Y1a0QMbOH8U5fvQ/cCCzBG5r9Mufc9p4rT0QiQWt7JxtLalkdCI3Vu/ZRVuedlkpMiGNGdjrXzs1jdm46s3IGMyIttFnyJPp0tyXyIN44WdfijZWVB7wHfMXMNjjnHu+R6kTEF7VNbazZtY/3CqtZtWsfH+ypoaXda2XkDPEmPzoxdzAn5g5m0sgU3bDXh3UrRJxzLwMvH9gOTEw1GS9YTgUUIiJRrLyumXcLqnmvsJqVBdVsLavHOe9O76lZaXxlzmjyRw/mpNGDGZ6qVoYc0t3TWacB77jADFbOuXZgfeDxSM+VJyI9obimiXd2VLGyoJqVhdUUVO4HYGD/fpw0ejAXTx9Jft5gZucMVue3HFN3T2ctBH5nZluB54HnnXN7e64sEQmn0tom3t5Rxds7qninoIo91U2AN6PeyXlD+PIpuZwyZghTR6X6OtWqRJ/uns76FwAzmwRcBDxoZmnAK3ih8qbmPheJHFUNLby9s4q3AsFxoKWRPjCBU8cM4bq5Yzh1zFAmjUghLk53fUvogr1jfQuwBfi1mSXh3R/yReBXQFDz8opI+DS1dvBuQRVvbq/kje1VbC6tAyBlQDynjh3Clafmcvq4DIWGhF3It4YG5lgf65z7VwAzS3fO1YSvNBE5ms5Ox8aSOlZ8WMHrH1awZlcNrR2d9O8Xx0mjB/PtCydy+rihTM9K0+kp6VHHO77A6C7rtwHfOc73E5GjKK9vZsW2SlZsq+CN7ZVU7/fGm5oyMpVr5+Yx94QMTs4boo5w6VXdvTprsnNu8xGeijOzM4E3gSFhrUykj2vv6GTN7hpe3VrOq1sr2BQ4RZWRPICzJwzjrAnDmHtCBsNSNKWP+Ke7LZFnzew14A7n3O4u+78N/AtwDfBUmGsT6XMqG1p4dWsFr2wt5/VtFdQ1txMfZ5w4ejDfmT+RT00YxuQRqerXkIjR3RCZhDdu1mtmthT4sXOuwjnXCdzTY9WJxDjnHJtL61m+uYzlW8r5oKgG52BYygDmTxvBOROHM3d8Bqkad0oiVHcv8W0F/sfM/gDcBKw0s0eAXzjn6nqyQJFY09reyds7q3hpUxnLN5dRUutNvjQzJ51vzZvAvMnDmTJSrQ2JDsFe4tsM3G1m9wI3A6vN7D7n3N09Up1IjKhrbuPVrRW8sHEvr22toL6lnaSEfpwxPoObzxvPOZOGMzxFw4lI9Al2Uqo8vFNbE4FcoB74KaAQETlMRX0LL24q4x8b9/LWjkraOhwZyf25ZMZIzp+SydwTMkhM0JVUEt26e3XWOiAL2I13s+FmYDnwv8C2YL6gmT0AXAqUO+emBfYNAf6CNzpwIXC5c26feRMo/xa4GGgErnHOrQnm64n0ptLaJp7fsJfn1u/lvV3VOAejhw7k2rljuHBqJrNyBtNPp6kkhnS3JXIZUHBgAMbj9CBe+DzcZd+twHLn3F1mdmtg+7t4Q6yMDzxOBe4NLEUiRnFNE8+tL+XZ9aW8v9u733ZiZgrfPHc886eNYNKIFLzPQyKxp7sd6zvD9QWdcysCp8W6WgCcHVh/CHgVL0QWAA8HwusdM0s3s5HOudJw1SMSitLaJp5dV8oz60pZu8cLjqmjUvn2hRO5aNoIxg5L9rlCkd5xvHesh0tml2DYC2QG1rOAPV2OKwrs+1iImNn1wPUAubm5PVep9FmVDS0sW1/KMx+UsrKwGvCC4zvzJ3LxtJHkZQzyuUKR3hcpIXKQc86ZWdCnzZxzi4BFAPn5+eE47SZCfXMb/9hYxtIPSnhzeyUdnY4Jmcnccv4ELpkxUi0O6fMiJUTKDpymMrORQHlgfzGQ0+W47MA+kR7T1tHJim0V/O39Yl7cVEZLeyfZg5O44ayxfGbWKCaNSPW7RJGIESkhshS4GrgrsHy6y/6bzOxRvA71WvWHSE9wzrGuqJYn1xSx9IMS9jW2MXhgAv90cg4LZmVxYm66OsdFjqDXQ8TM/ozXiZ5hZkXAHXjh8ZiZfRXYBVweOHwZ3uW92/Eu8b22t+uV2FZa28STa4p5ck0ROyr20z8+jvOnZPK52VmcNWEYCRpGXeSYej1EnHNXHOWpeUc41gE39mxF0tc0t3XwwqYyHl9dxBsfVtDp4JS8IXz9zLFcNH0kaUkap0qkuyLldJZIj9tQXMtjq/bw1PvF1DW3k5WexE3njufzJ2YxeqiurBIJhUJEYlptUxtPry3m0ZV72FRaR//4OC6aNoLL83M4bexQDXIocpwUIhJznHOs3rWPJSt3s2x9Kc1tnUwdlcqdC6bymZlZpA3U6SqRcFGISMyobWrjqfeLWfLubraW1ZM8IJ7PnZjNFSfnMj07ze/yRGKSQkSi3obiWh55ZxdPry2hqa2DGdlp3PW56Xx65igGDdB/cZGepN8wiUot7R0sW1/Kw2/v4v3dNSQmxLFgZhZfmTNarQ6RXqQQkahSWtvE4nd28+eVu6na38rYjEHcfukUPn9Sti7NFfGBQkQi3oGO8j+9WcjzG/fS6RzzJmVy9emjOeOEDN1JLuIjhYhErNb2Tp5ZV8Kf3ixkfXEtqYnxfPWMMVw1ZzQ5Qwb6XZ6IoBCRCFTT2Mrid3fz0FuFlNe3MG7YIO68bBqfPzGLgf31X1Ykkug3UiLGrqr9/PH1Av66eg/NbZ2cOT6D//rCDM4aP0w3BYpEKIWI+O793ftYtGInz2/cS0JcHAtmjeKrZ47RkOsiUUAhIr5wzvHq1grufW0HKwuqSU2M518+NY5rTs9jeGqi3+WJSDcpRKRXtXd08sy6Un7/2g627K1nZFoiP7hkMl86JZdk3RgoEnX0Wyu9ormtgyfWFPH713awp7qJ8cOTufuLM/nMzFH0j9ecHSLRSiEiPaqxtZ0l7+5m0YqdlNe3MDMnnf+4ZArnTc5UZ7lIDFCISI+ob27j4bd3cf8bBVTvb+W0sUP5zT/N4rRxQ3VzoEgMUYhIWNU1t/Hgm4Xc/0YBtU1tnD1xGP967gmcNHqI36WJSA9QiEhYHAiPP76+k7rmds6bnMk3553AjOx0v0sTkR6kEJHj0tDSzoNvFrBohRce50/J5OZ545mWpZF0RfoChYiEpLG1nYfe2sWiFTvY19jGeZMz+dZ5Cg+RvkYhIkFpae9gybu7ueeVHVQ2tHD2xGH8v/MmMDNHp61E+iKFiHRLe0cnj68u4r+Xf0hJbTNzxg7hvqtOVIe5SB+nEJFjcs7x3Ia93P3CVnZW7GdmTjr/9YWZzD1Bl+qKiEJEjuHN7ZX8/PktrCuqZfzwZBZddRLnT8lUeIjIQQoR+ZhNJXXc9fwWVmyrICs9ibu/OJPPzs6in+4wF5HDKETkoOKaJn75j638bW0xqYkJfP/iyVx12mgSE/r5XZqIRKiIChEzKwTqgQ6g3TmXb2ZDgL8AeUAhcLlzbp9fNcaiuuY2fvfKDh54swCA688ayzc+dQJpAxN8rkxEIl1EhUjAOc65yi7btwLLnXN3mdmtge3v+lNabGnr6GTxO7v47fIP2dfYxudmZ3HLhRPJSk/yuzQRiRKRGCKHWwCcHVh/CHgVhchxcc7x8pZyfrJsMzsr9nP6uKHcdvFk3SgoIkGLtBBxwAtm5oD7nHOLgEznXGng+b1Apm/VxYDNpXX8+NlNvLm9irHDBvHANfmcM3G4rrgSkZBEWoic4ZwrNrPhwItmtqXrk845FwiYjzGz64HrAXJzc3u+0ihT2dDCL1/Yxl/e201qUgI/+sxUvnxqLgn9NCGUiIQuokLEOVccWJab2d+AU4AyMxvpnCs1s5FA+VFeuwhYBJCfn3/EoOmLWts7efCtAv5n+Xaa2jq45vQx3DxvvDrNRSQsIiZEzGwQEOecqw+sXwD8J7AUuBq4K7B82r8qo8vLW8q485nNFFTu55yJw/j+JVM4YXiy32WJSAyJmBDB6+v4W+DcfDywxDn3vJm9BzxmZl8FdgGX+1hjVNhR0cCdz2zi1a0VjB02iAevPZmzJw73uywRiUEREyLOuZ3AzCPsrwLm9X5F0aehpZ3/Wf4hD7xZQGJ8P35wyWQWnpZH/3j1e4hIz4iYEJHQOed4em0JP122mfL6Fr54UjbfmT+JYSkD/C5NRGKcQiTKbSqp446lG3ivcB8zs9O476qTmJ072O+yRKSPUIhEqdqmNn794jYefruQtKQE7vrcdC7PzyFOgySKSC9SiEQZ5xxPrCnmruc2U7W/lStPzeXfL5hI+sD+fpcmIn2QQiSKbC6t4/anvVNXs3PTefDaUzRUiUhAW1sbRUVFNDc3+11KxEtMTCQ7O5uEhOO/X0whEgXqm9v4zUsf8uBbhaQmxvPzz0/niyfp1JVIV0VFRaSkpJCXl6dhfI7BOUdVVRVFRUWMGTPmuN9PIRLBnHM8s66UO5/ZREVDC1eckst3LtSpK5EjaW5uVoB0g5kxdOhQKioqwvJ+CpEItbOigduf3sj26veEAAALvElEQVQb2yuZlpXKooX5zMpJ97sskYimAOmecP6cdBdahGlu6+CXL2xl/m9e54OiGv5zwVSevvEMBYhIH5OXl0dlZWVIx3TnteGilkgEeWVLObcv3cCe6iY+OzuL7108ieEpiX6XJSJyVGqJRICSmiZu+L9VXPvge/TvF8efvz6HX//TLAWISBQpLCxk0qRJXHPNNUyYMIErr7ySl156iblz5zJ+/HhWrlwJQHV1NZdddhkzZsxgzpw5rFu3DoCqqiouuOACpk6dyte+9jWcOzQY+SOPPMIpp5zCrFmzuOGGG+jo6PDlezwStUR81NbRyQNvFPDb5R/S6RzfmT+Rr50xVmNdiRynH/19I5tK6sL6nlNGpXLHp6ce85jt27fz17/+lQceeICTTz6ZJUuW8MYbb7B06VJ++tOf8tRTT3HHHXcwe/ZsnnrqKV5++WUWLlzI2rVr+dGPfsQZZ5zB7bffzrPPPsv9998PwObNm/nLX/7Cm2++SUJCAt/4xjdYvHgxCxcuDOv3FyqFiE9WFlTzg6fWs62sgfMmZ3LHp6eQM2Sg32WJyHEYM2YM06dPB2Dq1KnMmzcPM2P69OkUFhYC8MYbb/DEE08AcO6551JVVUVdXR0rVqzgySefBOCSSy5h8GBv+KLly5ezevVqTj75ZACampoYPjxyRuVWiPSyyoYWfrpsM0+uKSYrPYk/LMzn/Cma8VcknD6pxdBTBgw4NOhpXFzcwe24uDja29tDek/nHFdffTU/+9nPwlJjuOm8SS/p6HT839uFnHv3q/z9gxK+cfY4Xvq3TylARPqYM888k8WLFwPw6quvkpGRQWpqKmeddRZLliwB4LnnnmPfvn0AzJs3j8cff5zycm9S1+rqanbt2uVP8UeglkgveH/3Pv7j6Q1sKK7j9HFD+c8FUzlheIrfZYmID374wx9y3XXXMWPGDAYOHMhDDz0EwB133MEVV1zB1KlTOf3008nNzQVgypQp/PjHP+aCCy6gs7OThIQE7rnnHkaPHu3nt3GQdb0CIFbk5+e7VatW+V0GVQ0t/OIfW3n0vT1kpg7gB5dM4dIZI3VDlEgP2Lx5M5MnT/a7jKhxpJ+Xma12zuUH8z5qifSA9o5OFr+7m1++sJXG1g6uP2ss35w3nuQB+nGLSGzRX7Uwe2dnFT9cupEte+s544QMfviZKTp1JSIxSyESJnuqG7nruS08u76UrPQk7r3yROZPG6FTVyIS0xQix6mhpZ37XtvBfSt2Emfwb+dP4PqzxpKY0M/v0kREepxCJETtHZ38dXURv3xhG5UNLSyYNYrvzp/EqPQkv0sTEek1CpEgOed4ZWs5dz23hW1lDeSPHswfFp7E7NzBfpcmItLrdLNhEFYWVHP5fW9z3YOraG3v5PdfOZG//vNpChARCbtwDQWfnJzcI/UdoJZIN6zdU8NvXtrGq1srGJ4ygJ98dhqX5+eQ0E8ZLCJ9m/4KHsOqwmoWPrCSy+55k7V7arj1okm89u1zuPLU0QoQEfkIDQUvgDfG1ctbyvnj6zt5t6CaoYP6c+tFk/jKnNG6WVAkWjx3K+xdH973HDEdLrrrmIdoKPg+rK65jSdXF/GntwrZVdXIqLREfnDJZL58ai4D++vHJCKfTEPBRygzmw/8FugH/NE5d+yPA93U2el4a0cVj6/ew/Mb99Lc1sns3HS+feFE5k8dQbxOWYlEp09oMfQUDQUfgcysH3APcBEwBbjCzKaE+n7NbR28urWc/3hqA3N//jJfuf9dlm8p5/MnZvP0jXP52zfmcumMUQoQEekRGgq+950CbHfO7QQws0eBBcCmT3phc1sHZXXNbCypY0NxLeuLa1lVuI+mtg6SEvpxxvgMvnfxZC6Ykqk7zEWkV2go+F5mZl8A5jvnvhbYvgo41Tl309Fek5w10WVf9xua2zoP7ouPMyZkpnDS6MGcO3k4p40dquAQiSEaCj44Ggr+MGZ2PXA9QNqosSw8LY+0pAQykvszaUQqE0ekKDRERMIsGkKkGMjpsp0d2PcRzrlFwCLwJqW67WJ9IhER6WnR0Hv8HjDezMaYWX/gS8BSn2sSERGioCXinGs3s5uAf+Bd4vuAc26jz2WJSARyzmkOn24IZ194xIcIgHNuGbDM7zpEJHIlJiZSVVXF0KFDFSTH4JyjqqqKxMTEsLxfVISIiMgnyc7OpqioiIqKCr9LiXiJiYlkZ2eH5b0UIiISExISEhgzZozfZfQ50dCxLiIiEUohIiIiIVOIiIhIyCJ+2JNQmFk9sNXvOo5DBnDseTEjVzTXDqrfb6rfXxOdcynBvCBWO9a3Bjv+SyQxs1XRWn801w6q32+q319mtirY1+h0loiIhEwhIiIiIYvVEFnkdwHHKZrrj+baQfX7TfX7K+j6Y7JjXUREekestkRERKQXxFSImNl8M9tqZtvN7Fa/6wmGmeWY2StmtsnMNprZzX7XFAoz62dm75vZM37XEiwzSzezx81si5ltNrPT/K4pGGb2/wL/dzaY2Z/NLDwj7PUQM3vAzMrNbEOXfUPM7EUz+zCwHOxnjcdylPp/Efj/s87M/mZm6X7WeDRHqr3Lc7eYmTOzjO68V8yEiJn1A+4BLgKmAFeY2RR/qwpKO3CLc24KMAe4McrqP+BmYLPfRYTot8DzzrlJwEyi6Pswsyzgm0C+c24a3rQJX/K3qk/0IDD/sH23Asudc+OB5YHtSPUgH6//RWCac24GsA34Xm8X1U0P8vHaMbMc4AJgd3ffKGZCBDgF2O6c2+mcawUeBRb4XFO3OedKnXNrAuv1eH/AsvytKjhmlg1cAvzR71qCZWZpwFnA/QDOuVbnXI2/VQUtHkgys3hgIFDicz3H5JxbAVQftnsB8FBg/SHgsl4tKghHqt8594Jzrj2w+Q7eTKwR5yg/e4BfA98But1ZHkshkgXs6bJdRJT9ET7AzPKA2cC7/lYStN/g/Qfs9LuQEIwBKoA/BU7H/dHMBvldVHc554qBu/E+QZYCtc65F/ytKiSZzrnSwPpeINPPYo7TdcBzfhfRXWa2ACh2zn0QzOtiKURigpklA08A33LO1fldT3eZ2aVAuXNutd+1hCgeOBG41zk3G9hPZJ9K+YhA38ECvDAcBQwys6/4W9Xxcd6lo1F5+aiZfR/vFPViv2vpDjMbCNwG3B7sa2MpRIqBnC7b2YF9UcPMEvACZLFz7km/6wnSXOAzZlaIdyrxXDN7xN+SglIEFDnnDrT+HscLlWhxHlDgnKtwzrUBTwKn+1xTKMrMbCRAYFnucz1BM7NrgEuBK1303EMxDu8DyAeB3+FsYI2ZjfikF8ZSiLwHjDezMWbWH69TcanPNXWbefN53g9sds79yu96guWc+55zLts5l4f3s3/ZORc1n4Sdc3uBPWY2MbBrHrDJx5KCtRuYY2YDA/+X5hFFFwZ0sRS4OrB+NfC0j7UEzczm453S/YxzrtHverrLObfeOTfcOZcX+B0uAk4M/F4cU8yESKAz6ybgH3i/PI855zb6W1VQ5gJX4X2CXxt4XOx3UX3MvwKLzWwdMAv4qc/1dFugBfU4sAZYj/e7HdF3T5vZn4G3gYlmVmRmXwXuAs43sw/xWld3+VnjsRyl/v8FUoAXA7/Dv/e1yKM4Su2hvVf0tLZERCTSxExLREREep9CREREQqYQERGRkClEREQkZAoREREJmUJERERCphAREZGQKUREeoGZ3WBme83sAzPbYWYL/a5JJBx0s6FILzCz/wU2OOd+b2anAMucc92a9EckkqklItI7ZgBbA+sFQKuPtYiEjUJEpHdMB7YGBke8Cfi+z/WIhIVOZ4n0sMCUowXABryJ0tYB50XRMOEiR6WWiEjPmw6scM7NAiYAk4DT/C1JJDwUIiI9bwbwPoBzbh+wBG8uepGopxAR6XnTCYRIwN8BzRUjMUF9IiIiEjK1REREJGQKERERCZlCREREQqYQERGRkClEREQkZAoREREJmUJERERCphAREZGQ/X8LQdNLu09cfgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "V_de_R = np.sqrt(C*G*rho_h_I*Vc_R(R))\n",
    "V_de_R_II = np.sqrt(C*G*rho_h_I*Vc_R_II(R))\n",
    "plt.plot(R,V_de_R,label='model I')\n",
    "plt.plot(R,V_de_R_II,label='model II')\n",
    "xlim(0.,14.)\n",
    "plt.xlabel(r'$R$')\n",
    "plt.ylabel(r'$V_{c}/km s^{-1}$')\n",
    "plt.legend(loc='lower right', prop={'size':10})\n",
    "plt.savefig('DMHalo_Models_I_II.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}