ChiSquare_Model_I_v4_8p.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",
    "from scipy import integrate\n",
    "import galkin\n",
    "import galkin.processdata   # routines to process kinematic data\n",
    "import galkin.readparsFile  # routines to read and check input parameters\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Constants\n",
    "G = 4.302e-6#kpc/SM(km/s)^2\n",
    "q_b = 0.6\n",
    "r_b = 1.9#kpc\n",
    "a_b = 1.#kpc\n",
    "q_h = 0.8\n",
    "alpha_b = 1.8\n",
    "R_m = 4.#kpc\n",
    "#parameters\n",
    "rho_b_I = 0.427e9#sM/kpc^3\n",
    "a_h_I = 3.83#kpc\n",
    "rho_h_I = 0.711e9#SM/kpc^3\n",
    "alpha_h_I = -2.\n",
    "beta_h_I = 2.96\n",
    "R_d_I = 2.0#kpc\n",
    "Sigma_d_SD_I = (1905.0e6)*0.75#SM/kpc^2\n",
    "Sigma_g_ISM_I = (1905.0e6)*0.25#SM/kpc^2\n",
    "params = np.array([rho_b_I,rho_h_I,alpha_h_I,a_h_I,beta_h_I,R_d_I,Sigma_d_SD_I,Sigma_g_ISM_I])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Vc_Model_I(R,params):\n",
    "    rho_b,rho_h,alpha_h,a_h,beta_h,R_d,Sigma_d_SD,Sigma_g_ISM = params\n",
    "    def Vc_b(R):\n",
    "        def I_b(x,R):\n",
    "            e = np.sqrt(1. - q_b**2.)\n",
    "            n = x**(2.-alpha_b)*np.exp(-(R**2.*x**2.)/r_b**2.)\n",
    "            d = np.sqrt(1.-x**2.*e**2.)\n",
    "            t = R**2.*(R/a_b)**(-alpha_b)\n",
    "            return (t*n)/d\n",
    "        Integral_b = []\n",
    "        for i in range (0,len(R)):\n",
    "            result = integrate.quad(I_b,0,1,args=(R[i]))[0]\n",
    "            Integral_b.append(result)\n",
    "        return np.array(Integral_b)\n",
    "    def Vc_DM(R):\n",
    "        def I_h(x,R):\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",
    "        Integral_DM = []\n",
    "        for i in range (0,len(R)):\n",
    "            result = integrate.quad(I_h,0,1,args=(R[i]))[0]\n",
    "            Integral_DM.append(result)\n",
    "        return np.array(Integral_DM)\n",
    "    def V2_c_SD(R):\n",
    "        result = (R_d/R) - np.exp(-R/R_d) - (np.exp(-R/R_d)*R_d/R)\n",
    "        return result\n",
    "    def Vc_ISM(R):\n",
    "        def M_g_ISM(u,R):\n",
    "            return R*u*np.exp(-((R*u)/(2.*R_d))-(R_m/(R*u)))\n",
    "        Integral_M_g = []\n",
    "        for i in range (0,len(R)):\n",
    "            result = integrate.quad(M_g_ISM,0,1,args=(R[i]))[0]\n",
    "            Integral_M_g.append(result)\n",
    "        return np.array(Integral_M_g)\n",
    "    suma = 4.*np.pi*G*q_b*rho_b*Vc_b(R) + 4.*np.pi*G*q_h*rho_h*Vc_DM(R) + 2.*np.pi*G*Sigma_d_SD*R_d*V2_c_SD(R) + 2.*np.pi*Sigma_g_ISM*G*Vc_ISM(R)\n",
    "    return np.sqrt(suma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = np.linspace(0.,15.,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  \n",
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  \n",
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The occurrence of roundoff error is detected, which prevents \n",
      "  the requested tolerance from being achieved.  The error may be \n",
      "  underestimated.\n",
      "  warnings.warn(msg, IntegrationWarning)\n",
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:29: RuntimeWarning: divide by zero encountered in divide\n",
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:29: RuntimeWarning: invalid value encountered in subtract\n",
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:33: RuntimeWarning: divide by zero encountered in double_scalars\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4XOV99vHvT7u1WJItWdiSbNnYeAHbLDZmCXvSsDWmSZOSkoQQiLOQraVNQugL7fu+aUmbNk3alJaENSFQkkAgLEmAEKAJBmzjDbxgbMuWJVuyrM3aNfPrH3MMQpZsjZHmjDT357rmmpkzM/JtjHTrPOc8zzF3R0REZLjSwg4gIiJji4pDRETiouIQEZG4qDhERCQuKg4REYmLikNEROKi4hARkbioOEREJC4qDhERiUtG2AFGQ0lJiVdVVYUdQ0RkTFm9evV+dy892vvGZXFUVVWxatWqsGOIiIwpZlY9nPdpqEpEROKi4hARkbioOEREJC4qDhERiYuKQ0RE4qLiEBGRuKg4REQkLuNyHoeISDKLRp2eSJTu3ijdkQg9fdHYLRKlty/2Wk9flN7IoZu/9bgvEns9Eo1t64v6W48jUee698yiMDdzVPOrOERE+umLRGnvidDe3UdHTx/t3RHae/ro7InQ3hOhs6ePjp4Inb0Runoibz3u7I3Q1RuhsydCV2+Urr7YfXdfJFYQwX1XX4TeiI9KdjP48GmVKg4RkeHq6o3Q0tlLa2cvrV29tHb2Bfe9tHb10dbVx8Hu3th9Vx9t3bH79p7Y/cHuPrr7osP+88wgNzOdnOA2ISudnMw0cjLSycvKYHJeGtmZ6eRkpJMdbM/OTCMrPY3szDSyM9LJykgjOz2NrIw0MoP72GMjKz22LXYzMtPTyDh0n2ZkBNsz0mLP09JsFP/rvk3FISJJpy8SpamjlwPtPRxo76GpI3Zr7uilqb2H5s5emjt6aemMbWvu7KWls5eeo/zQz0w3CnIyKcjJID87dptWlENedgZ52RkUBPd52RnkZcWKID87gwlZsSLIzUonNzuD3KAksjPSMEvMD+tkouIQkVHn7rT3RGho66ahrZv9Bw/deth/sJvGg900HoyVRGN7Dy2dvUN+rQmZ6RTnZlKYm0XRhExmT8mnKDeTiRMyKZyQycSc2H3hhNi2gpwMCnIymJiTSU5megL/1uOXikNE3pX27j72tnaxr6WLfW1d7GvtZl9rF/Wt3dS3dVHf1k19azedvZHDPmsGk3KzmJyfxaS8LOZPm0hJXhaT8rKZlJ9FcW4mxbmx14pzsyjK1Q//ZKDiEJEhdfdF2NvSxZ6mTmpbuqht7qS2uZO6li7qWmL3bV19h30uPzuDKQXZTJmYzeKKIkoLsplSkE1pQTYl+W/fT8rLIj1B4/IyclQcIimsNxKltrmT3Qc62d3Uwe4DHdQ0dVLTFLtvONiNDzgBqCQ/i6mFE6ianMeZsyZTVpjD1MIcyia+fcvP1o+W8Uz/uiLjXFdvhOrGDnY2tlPd2E51Y0fsdqCd2uYuItG3myEjzZhWNIGK4gmcd0IpFcW5TCvKobxoAuXFEyibmKOhIlFxiIwH0ahT19rFm/UHebPhIDv2t7O9oZ0d+9upbel8x15DcW4m0yfncUplMcsX5zJ9Ui6Vk3KpnDSBqYUTNHQkR6XiEBlDIlFn94EOtuxrY1v9QbbVH+SN+jberG9/x8HnguwMZpXmsbSqmJkllVSV5FI1OY+qyXmjPjlMxj8Vh0iSqm/rYnNdG1v2trF5bxtb9rXyxr6D75igNq0wh9llBSw9fRKzp+RzfGk+s0rzKM3PTsn5BZIYKg6RkEWizo79B3mttpXXa1t5va6VTXWt7D/Y89Z7phRkM/e4Aj5+xgxOKCtgTlk+s6fkU5CjvQdJPBWHSAJFos6bDQdZX9PChppmNgZlcWiYKSs9jTll+Vwwdwrzp05k3tQC5h83keK8rJCTi7xNxSEyStyd2pYu1u5qZl1NM2t3N7NxTwsdPbGSyM1K58RpE/mzpZWcVF7IidMmMntKPpnputqBJDcVh8gI6eqNsHFPC6urm1izq4lXdzVT39YNQFZGGidOm8iHT6tgUUURiysLmVmSrzOYZExScYgco+aOHlbtbOKV6gO8suMAG/e00hOJHbieMTmXs2eXcMr0Ik6uLGLecRPJytCehIwPKg6RYdp/sJuXth/gpR2NvLzjAJv3tgGxFVcXVRRxzdlVnDajmFNnFFOSnx1yWpHRo+IQGUJLZy8rtzfy4puN/OHN/WzddxCIHZs4bUYxly2cyukzJ7G4skizqSWlqDhEAj19UV7d1cT/bNvPC2/sZ31NM1GHnMw0llZNYvnJ5Zx5/GQWlhfqALakNBWHpLTdBzr43dYGntvSwItv7qe9J0J6mrG4opAvXDA7OE5RrOMTIv2oOCSl9EaivLLjAM9uqee3m+t5s6EdgIriCVxxSjnnzCnlzOMnUzhBE+tEhqLikHGvpaOXZ7fU8/SmfTy3pYG27j6y0tNYNmsSVy2bwflzS5lZkqclOkSGKeHFYWaVwL1AGeDA7e7+XTObBPw3UAXsBD7i7k0W+27+LnAp0AF80t3XJDq3jC21zZ089fo+fvP6XlZuP0Ak6pTkZ3PpwqlcOH8K75ldQp6uGSFyTML4zukDbnD3NWZWAKw2s6eATwLPuPutZvZ14OvA14BLgDnBbRlwW3Av8g7Vje08uXEvT27cy7rdzQDMnpLPZ86dxXsXlHFyRRFpmnAn8q4lvDjcvQ6oCx63mdkmoBxYDpwfvO0e4HfEimM5cK+7O7DSzIrMbGrwdSTFVTe28/iGOh5fX8drta0ALKoo5KsXz+X9Jx7H8aX5IScUGX9C3Vc3syrgFOAloKxfGewlNpQFsVLZ3e9jNcG2dxSHma0AVgBMnz591DJL+Pa2dPHY+lp+ua6WdTUtAJxcWcRNl87n4pOOo3JSbsgJRca30IrDzPKBnwNfcffW/gcm3d3NzIf88CDc/XbgdoAlS5bE9VlJfi2dvfxqYx2/eLWWlTsacYeTyidy4yXzuGzRVCqKVRYiiRJKcZhZJrHSuM/dHwo27zs0BGVmU4H6YPseoLLfxyuCbTLO9UaiPL+1gYfW7OGpTfvo6YsysySPL180hw8snsYsDUOJhCKMs6oMuAPY5O7/0u+lR4GrgVuD+0f6bf+CmT1A7KB4i45vjG9b9rbx01W7+cXaPew/2MOkvCz+/PTp/Mkp5SyqKNRpsyIhC2OP42zg48AGM1sbbPsGscJ40MyuBaqBjwSvPUHsVNxtxE7HvSaxcSURWrt6eXRtLQ+u2s36mhYy042L5pXxodMqOH9uqZb4EEkiYZxV9T/AUL8yXjTI+x24flRDSSjcndXVTdz/8m4e31BLV2+UeccVcPPlC7jilHIm6ap3IklJM6Ak4Vo6e3loTQ33v7yLrfsOkp+dwQdPreDKpZUsLNdQlEiyU3FIwmzc08KPXqzmkXV76OqNsriikG99aCGXL5qmWdwiY4i+W2VUdfdFeHLDXu7+w07W7m5mQmY6f3JKOVctm8FJ5YVhxxORY6DikFFR39rFj1dW85OXd7H/YA+zSvK45Y8X8KHTKpiYo5VnRcYyFYeMqA01Ldz5+x08tr6Wvqhz4dwpXH1WFe+ZXaJ1okTGCRWHvGuRqPP0pn3c8cIOXt55gPzsDD52xgyuPrOKqpK8sOOJyAhTccgx6+qN8NPVNdzxwnZ2NnZQXjSBv7lsPn+2tJICDUeJjFsqDolbU3sP975Yzb0v7qSxvYfFFYV8/89P5f0nlpGhiXoi456KQ4atrqWTH76wg/tf3kVHT4QL501hxbmzWDZzkuZeiKQQFYccVXVjO7f97k1+vqaGqMMHFk/js+cdz9zjCsKOJiIhUHHIkLbVH+T7z27jkbV7yEhP48ql01lx7ixd70Ikxak45DDb6tv43jPb+OX6WnIy0rn2PTP59DmzmDIxJ+xoIpIEVBzylh372/nXp7fy6LpaJmSm85lzj+fT58xkcn522NFEJImoOIQ9zZ187+k3+NmaGjLTjRXnzmLFObNUGCIyKBVHCmtq7+E/freNe16sBodPnDmDz51/PFMKNCQlIkNTcaSgrt4Id/5+B7c9+ybtPX186NQKvvK+EygvmhB2NBEZA1QcKcTdeWx9Hbc+uZk9zZ1cNG8KX7tkHieU6bRaERk+FUeKeL22lZsf2ciq6ibmHVfAT65bxlmzS8KOJSJjkIpjnGvp7OU7T23l3hd3UpSbxT98cCEfWVJJulaqFZFjpOIYp9ydJzbs5ZZHX+NAezcfO2MGN7xvLoW5WnxQRN4dFcc4VNfSyf/5xWs8vWkfC8sLufuapbranoiMGBXHOOLuPLRmD3/76Gv0RqPcdOl8rjm7SivWisiIUnGMEwfae7jp4Q08uXEvS6uK+faHFzNjsi6iJCIjT8UxDqzc3siX7n+Vpo4evnbxPFacO0sHv0Vk1Kg4xrBo1Pmv57fzT7/eTNXkPO66ZiknTtOxDBEZXSqOMaq1q5e//O+1PL2pnssWTeVbH1pEfrb+OUVk9OknzRi0+0AH197zCtsb2rnljxfwybOqdAU+EUkYFccY8+quJj597yq6+6Lc86nTOVuzv0UkwVQcY8gzm/bx+fvWMGViNg+sOIPZU7TGlIgkXsJP8DezO82s3sw29tt2spmtNLO1ZrbKzE4PtpuZfc/MtpnZejM7NdF5k8Xj6+v4zI9WM/e4An7x+bNVGiISmjBmht0NXDxg2z8Cf+fuJwM3B88BLgHmBLcVwG0JyphUHn61hi/ev4aTK4u477plusCSiIQq4cXh7s8DBwZuBiYGjwuB2uDxcuBej1kJFJnZ1MQkTQ4PranhLx9cx7KZk7nnU6dTkKO1pkQkXMlyjOMrwK/N7NvEyuysYHs5sLvf+2qCbXWJjReO327ex1//bD1nHT+ZO65eSk5metiRRERCGaoazOeAv3D3SuAvgDvi/QJmtiI4PrKqoaFhxAMm2ppdTXz+vjUsmDqR//r4EpWGiCSNZCmOq4GHgsc/BU4PHu8BKvu9ryLYdhh3v93dl7j7ktLS0lELmgjb6tv41N2vcNzEHO66Zqkm9olIUkmW4qgFzgseXwi8ETx+FPhEcHbVGUCLu4/rYaqWzl6uu2cVGWlp3PupZZToQLiIJJmE/yprZvcD5wMlZlYD3AJ8GviumWUAXcTOoAJ4ArgU2AZ0ANckOm8iRaPODQ+upaapkwdWnMH0yblhRxIROUzCi8PdPzrES6cN8l4Hrh/dRMnj+89u4+lN9fztHy9gSdWksOOIiAwqWYaqUt5zWxv4l6e3csXJ07j6rKqw44iIDEnFkQQOtPdww4NrmVtWwD98cJEWLBSRpKbTdZLAzY9spKWzlx9du4wJWTrtVkSSm/Y4Qvb4+joeW1/Hly+aw/ypE4/+ARGRkKk4QtTQ1s3f/GIDiysK+ex5x4cdR0RkWFQcIbrl0Y2090T4548sJiNd/xQiMjbop1VI/vDmfp7YsJcvXThbS6SLyJii4ghBJOr831++TkXxBK47Z1bYcURE4qLiCMEDr+xi8942vnHpfC1eKCJjTtzFYWZfG40gqaKls5d//s1WTp85iUtOOi7sOCIicTvqPA4ze7D/U+Bk4Fujlmic+7dn3qCpo4ebL1+giX4iMiYNZwJgq7tfd+iJmaXk5VtHwt6WLu59sZo/PbWCk8oLw44jInJMhjNU9c0Bz28ajSCp4PbntxNx50sXzQk7iojIMTtqcbj7DgAzKwmeD7xeuAzD/oPd/OTlaq44uZzKSVouXUTGrngOjt85ailSwA9f2EF3X5TPX6AZ4iIytsVTHDqSe4yaO3r40Ys7uXzRNI4vzQ87jojIuxJPcfiopRjn7vr9Ttp7IlyvvQ0RGQe0xzHKOnr6uOv3O3jfgjLmHafVb0Vk7IunOG4ctRTj2GPr6mjt6uPTWlpERMaJYReHu28czSDj1Y9fquaEsnyWVhWHHUVEZETEteSImT1vZhODx581s6+YWdboRBv71tc0s76mhauWzdAscREZN+Jdq6rQ3VvN7DTg00Ax8IORjzU+/OSlXUzITOdPTi0PO4qIyIiJ95rjvWaWAXwC+Ja7P2hmq0Yh15jX2tXLI2trWX7yNCbmZIYdR0RkxMRbHP8GrANygK8H2zQxYRAPr9lDZ2+Eq5bNCDuKiMiIiqs43P0eM3sY6HX3TjObDbw4OtHGLnfnxyurWVxRyMIKLWYoIuNLXMVhZj8ALgMiZlYLrAfWjEawsWx9TQtv1B/k1g8uDDuKiMiIi3eo6lyg0t0jZlYOLAYWjXysse2x9bVkphuXnDQ17CgiIiMu3uJ4CZgM1Lv7HmAP8MSIpxrDolHn8fV1nDunlMJcHRQXkfEn3tNx/wt4zsz+yszOMTMN4A/w6u4malu6uHyx9jZEZHyKtzh+DNxLbE/l88AfzOzNeL6Amd1pZvVmtnHA9i+a2WYze83M/rHf9hvNbJuZbTGz98eZN+F+ua6OrIw03ju/LOwoIiKjIt6hqhp3/4f+G8wsO86vcTfw78QK6NDXuABYDix2924zmxJsXwBcCZwITAOeNrMT3D0S55+ZEJGo88SGOi6YW0qB5m6IyDgV7x7HWjP7cv8N7t4dzxdw9+eBgVcR/Bxw66Gv5e71wfblwAPu3h1ciXAbcHqcmRPmlZ0HqG/r5vJF08KOIiIyauItjjLgs2ZWa2aPmdk3zezDI5DjBOAcM3vJzJ4zs6XB9nJgd7/31QTbktJj62uZkJnORfOnhB1FRGTUxDsB8CPw1vDUicBCYBnw0xHIMQk4A1gKPGhmca1DbmYrgBUA06dPf5dx4tcXifKrjXu5cP4UcrPiHQEUERk74l0d9/9DbHjK3de4+z3A10YgRw3wkMe8DESBEmKn+1b2e19FsO0w7n67uy9x9yWlpaUjECk+q6qb2H+wh8sX6mwqERnf4h2qKjezPz/0xMxKgadHIMcvgAuCr3kCkAXsBx4FrjSzbDObCcwBXh6BP2/EPbe1gYw045wTEl9aIiKJFO+YymeAX5vZNmLXIL+LOPc4zOx+4HygxMxqgFuAO4E7g1N0e4Cr3d2B18zsQeB1oA+4PlnPqHpuSwNLqorJz9YwlYiMb8P6KWdm9xJbk+pV4HrgJ8R+kF/h7tvi+QPd/aNDvPSxId7/TeCb8fwZiVbf1sXrda189eK5YUcRERl1wx2quhsw4BpikwCrgCbgY2b2p6OSbAx5Yet+AM7TMJWIpIBh7XG4+2+B3x56HlzMaT6xRQ6XAT8blXRjxHNbGygtyGbB1IlhRxERGXXDHao6E1gZHHfA3fuADcHtx6MXL/lFos4LbzRw4bwyXVdcRFLCcIeqPgGsMbMHzOyTZnbcaIYaSzbsaaGpo5fz5mqYSkRSw3CHqj4HYGbzgEuAu4OVcZ8FfgX8PlnPdhptz21pwAzOmV0SdhQRkYSIax6Hu2929++4+8XAhcD/AB8mdp2OlPTc1noWVRRRnJcVdhQRkYSIdwLgW9y9E5jl7l909yVmVjSCucaElo5e1u5u1tlUIpJSjrk4AjP6Pf7Gu/xaY86L2/cTdTjvBA1TiUjqGFZxmNn8oT4fXAkwjdgihSnllZ1NZGeksbA85Xa2RCSFDXeP43Ezu8vMBi47+9fAIuAHwCMjmmwMWFXdxOLKIrIy3u2Om4jI2DHcn3jziC058pyZfTdY3BB3j7r79939Wnf/5ailTEJdvRFe29PCaTOKw44iIpJQwyoOd+9x938jNlt8N/Cymf0/M0vZqdLrdjfTF3WWqDhEJMXEezpul7t/GzgJ6ARWm9lfjUqyJLd6VxMAp05XcYhIaon3Qk5VZnYxcB0wHWgD/n40giW71TubOL40T/M3RCTlDPesqvVm1gg8DFwNFALPEFuKJH/04iWnaNRZvauJJTNS7kQyEZFhX8jpCmDHoUUOU932/e00d/TqwLiIpKThrlW1fbSDjCWrqw8AcFqVikNEUo8mIByDVTubKM7NZFZJXthRREQSTsVxDFbvauK0GcW6/oaIpCQVR5wOtPewvaGd03RgXERSlIojTmuqY/M3dGBcRFKViiNO62uaSTNYWF4YdhQRkVCoOOL0el0bM0vymJCVHnYUEZFQqDjitHlvK/OnpuwSXSIiKo54tHb1UtPUqeIQkZSm4ojD5ro2ABaoOEQkhak44rCprhVAexwiktJUHHHYVNdKcW4mZROzw44iIhIaFUccNtXFDoxrxriIpDIVxzBFos6WfW0aphKRlJfw4jCzO82s3sw2DvLaDWbmZlYSPDcz+56ZbQuuCXJqovMesmN/O129URWHiKS8MPY47gYuHrjRzCqBPwJ29dt8CTAnuK0AbktAvkG9fWC8IKwIIiJJIeHF4e7PAwcGeek7wFeB/heLWg7c6zErgSIzm5qAmIfZVNdKRpoxe0rKXfBQROQdkuIYh5ktB/a4+7oBL5UDu/s9rwm2DfY1VpjZKjNb1dDQMOIZN9W1MntKPtkZWmpERFJb6MVhZrnAN4Cb383Xcffb3X2Juy8pLS0dmXD9bKrTgXERERj+NcdH0/HATGBdcJprBbDGzE4H9gCV/d5bEWxLqKb2Hva2dun4hogISbDH4e4b3H2Ku1e5exWx4ahT3X0v8CjwieDsqjOAFnevS3RGzRgXEXlbGKfj3g+8CMw1sxozu/YIb38C2A5sA34AfD4BEQ/zuopDROQtCR+qcvePHuX1qn6PHbh+tDMdzZsN7RTnZlKSr6VGRERCH6oaC6ob26kqyQs7hohIUlBxDMPO/e3MnKziEBEBFcdRdfVGqG3pYoaKQ0QEUHEc1e4DHQBUleSGnEREJDmoOI5iZ2OsOLTHISISo+I4ip372wF0jENEJKDiOIqdje0U5WZSmJsZdhQRkaSg4jiK6sYODVOJiPSj4jiKnY3tVE3WgXERkUNUHEfQ3RehtrmTKu1xiIi8RcVxBLsPdBJ1nYorItKfiuMIqhtjZ1TpGIeIyNtUHEdwaA6HhqpERN6m4jiCnfvbmZiTQbFOxRUReYuK4wh2BqviBlcmFBERVBxHpDkcIiKHU3EMoacvSk1TBzM1h0NE5B1UHEOoaeog6jqjSkRkIBXHEKobtZy6iMhgVBxD2Kk5HCIig1JxDKG6sYOC7Awm52WFHUVEJKmoOIawp7mTaUUTdCquiMgAKo4h1Ld2UVaYE3YMEZGko+IYQn1bN1MKssOOISKSdFQcg4hGnfq2bsomqjhERAZScQyisb2HSNQpm6ihKhGRgVQcg6hv6wLQUJWIyCBUHIOob+0GYIr2OEREDqPiGMS+1tgeh4aqREQOl/DiMLM7zazezDb22/ZPZrbZzNab2cNmVtTvtRvNbJuZbTGz9yciY31bbI+jNF9DVSIiA4Wxx3E3cPGAbU8BJ7n7ImArcCOAmS0ArgRODD7zH2aWPtoB97V2MSkvi6wM7ZCJiAyU8J+M7v48cGDAtt+4e1/wdCVQETxeDjzg7t3uvgPYBpw+2hn3tWoOh4jIUJLxV+pPAU8Gj8uB3f1eqwm2jar6ti4dGBcRGUJSFYeZ3QT0Afcdw2dXmNkqM1vV0NDwrnLUt3ZTpj0OEZFBJU1xmNkngcuBq9zdg817gMp+b6sIth3G3W939yXuvqS0tPSYc0SiTsPBbqZo1riIyKCSojjM7GLgq8AH3L2j30uPAleaWbaZzQTmAC+PZpbG9m7NGhcROYKMRP+BZnY/cD5QYmY1wC3EzqLKBp4KljFf6e6fdffXzOxB4HViQ1jXu3tkNPO9NfmvQMUhIjKYhBeHu390kM13HOH93wS+OXqJ3umt5UY0VCUiMqikGKpKJvuCPQ4NVYmIDE7FMcChoSrNGhcRGZyKY4B9bZo1LiJyJPrpOEB9a5dmjYuIHIGKY4DYlf90fENEZCgqjgH2aY9DROSIVBz9RKJOg/Y4RESOSMXRT+PBbqIOZZrDISIyJBVHP29dwEmzxkVEhqTi6OftS8Zqj0NEZCgqjn40a1xE5OhUHP0cWqeqRLPGRUSGpOLoZ19rN5M1a1xE5Ij0E7Kf+lZdMlZE5GgSvqx6MltcWcSCadGwY4iIJDUVRz9fumhO2BFERJKehqpERCQuKg4REYmLikNEROKi4hARkbioOEREJC4qDhERiYuKQ0RE4qLiEBGRuJi7h51hxJlZA1Adx0dKgP2jFGekKOPIUMaRoYwjI9kyznD30qO9aVwWR7zMbJW7Lwk7x5Eo48hQxpGhjCNjLGQcjIaqREQkLioOERGJi4oj5vawAwyDMo4MZRwZyjgyxkLGw+gYh4iIxEV7HCIiEpeULg4zu9jMtpjZNjP7eth5BjKzSjN71sxeN7PXzOzLYWcaipmlm9mrZvZY2FkGY2ZFZvYzM9tsZpvM7MywMw1kZn8R/DtvNLP7zSwpLkdpZneaWb2Zbey3bZKZPWVmbwT3xUmY8Z+Cf+/1ZvawmRUlW8Z+r91gZm5mJWFki1fKFoeZpQPfBy4BFgAfNbMF4aY6TB9wg7svAM4Ark/CjId8GdgUdogj+C7wK3efBywmybKaWTnwJWCJu58EpANXhpvqLXcDFw/Y9nXgGXefAzwTPA/T3Rye8SngJHdfBGwFbkx0qAHu5vCMmFkl8EfArkQHOlYpWxzA6cA2d9/u7j3AA8DykDO9g7vXufua4HEbsR925eGmOpyZVQCXAT8MO8tgzKwQOBe4A8Dde9y9OdxUg8oAJphZBpAL1IacBwB3fx44MGDzcuCe4PE9wBUJDTXAYBnd/Tfu3hc8XQlUJDzYO/MM9t8R4DvAV4Exc8A5lYujHNjd73kNSfhD+RAzqwJOAV4KN8mg/pXY//jJesH2mUADcFcwnPZDM8sLO1R/7r4H+Dax3zrrgBZ3/024qY6ozN3rgsd7gbIwwwzDp4Anww4xkJktB/a4+7qws8QjlYtjzDCzfODnwFfcvTXsPP2Z2eVAvbuvDjvLEWQApwK3ufspQDvhD628Q3CMYDmxkpsG5JnZx8JNNTweOzUzaX9bNrObiA373hd2lv7MLBf4BnBz2FnilcrFsQeo7Pe8ItiWVMzug+PfAAADjklEQVQsk1hp3OfuD4WdZxBnAx8ws53EhvsuNLMfhxvpMDVAjbsf2lv7GbEiSSbvBXa4e4O79wIPAWeFnOlI9pnZVIDgvj7kPIMys08ClwNXefLNPTie2C8K64LvnwpgjZkdF2qqYUjl4ngFmGNmM80si9iByEdDzvQOZmbExuU3ufu/hJ1nMO5+o7tXuHsVsf+Gv3X3pPpN2d33ArvNbG6w6SLg9RAjDWYXcIaZ5Qb/7heRZAfwB3gUuDp4fDXwSIhZBmVmFxMbQv2Au3eEnWcgd9/g7lPcvSr4/qkBTg3+f01qKVscwUGzLwC/JvYN+qC7vxZuqsOcDXyc2G/xa4PbpWGHGqO+CNxnZuuBk4G/DznPOwR7Qz8D1gAbiH1vJsWsYjO7H3gRmGtmNWZ2LXAr8D4ze4PY3tKtSZjx34EC4Knge+c/kzDjmKSZ4yIiEpeU3eMQEZFjo+IQEZG4qDhERCQuKg4REYmLikNEROKi4hARkbioOEREJC4qDpFjZGafMbO9ZrbOzN40s0/0e+0/zezs4PFFZvaj8JKKjCwVh8ixWwj8rbsvBj4K9F8W5gxiS3lD7PofryY4m8ioUXGIHLtFwJbg8Q6gB8DM5gNb3T0SvLYYeNXMss3sbjP7e4u538z+28xeNrNqM7ss+Pw0M/t5sAT8ZjM7PdF/MZEjUXGIHLuFwJZgUcIvADcF2y8BftXvfYuIrR77a+Bpd/9GsFLrYmC7u58OXAXcElzE6UngrmAJ+FNJ7sUOJQVprSqRYxBc7nMHsJHYBcDWA+91dzezXwPXuHttsCz+fqAa+Iy7vxh8PofYhcQq3b3LzCYRu0jXjcCfunuyXDZW5DDa4xA5NguB5939ZOAEYB5wZnBxniJ3P3TZ1/nElvDvAyL9Pn8S8Ia7dwXPTwXWEVu5dyUiSUzFIXJsFhEc8Hb3JuAnxK67fgHwbL/3LQb+QOxaJXeZWVm/7dPNLCe4jO3fEbv29F7gxEMfNrPSUf57iMRNxSFybBbyzjOlfglcyuHHNxYDG919K/A14MFg+Goxsav8vURsj+Q2d/89cDdQZmavmdla4MzR/ouIxEvHOERGkJmtAZYFl3890vueA1a4+5YjvU8kGak4REJgZjXAdHePhp1FJF4qDhERiYuOcYiISFxUHCIiEhcVh4iIxEXFISIicVFxiIhIXFQcIiISFxWHiIjERcUhIiJx+V8IUJv4nIIJoQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(R,Vc_Model_I(R,params))\n",
    "plt.xlabel(r'$R/kpc$')\n",
    "plt.ylabel(r'$V_{c}/km s^{-1}$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#10.3, 15.3, 7.7from the paper\n",
    "# Galactic parameters\n",
    "R0=8.0# Galactocentric distance (kpc)\n",
    "V0=220.# local circular velocity (km/s)\n",
    "UsunINUSE=10.3#11.10# solar motion in the U-direction (km/s), e.g. from Schoenrich+ '10, MNRAS 403, 1829 (2010)\n",
    "VsunINUSE=15.3#12.24# solar motion in the V-direction (km/s), e.g. from Schoenrich+ '10, MNRAS 403, 1829 (2010)\n",
    "WsunINUSE=07.7#07.25# solar motion in the W-direction (km/s), e.g. from Schoenrich+ '10, MNRAS 403, 1829 (2010)\n",
    "SYSTDISP=0.# systematic dispersion due to spiral arm streaming (km/s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flags\n",
    "flagPROPERMOTIONS=0\t \t\t\t# proper motions not supported in current version - please keep at 0\n",
    "flagHITERMINAL=0\t\t\t\t# whether to use HI terminal velocities\n",
    "flagFich89tab2=0\t\t\t\t# whether to use Fich+ '89, ApJ 342, 272 (1989) (Table 2)\n",
    "flagMalhotra95=0\t\t\t\t# whether to use Malhotra '95, ApJ 448, 138 (1995)\n",
    "flagMcClureGriffithsDickey07=0\t\t\t# whether to use McClure-Griffiths & Dickey '07, ApJ 671, 427 (2007)\n",
    "flagHITHICKNESS=0\t\t\t\t# whether to use the HI thickness method\n",
    "flagHonmaSofue97=0\t\t\t\t# whether to use Honma & Sofue '97, PASJ 49, 453 (1997)\n",
    "flagCOTERMINAL=0\t\t\t\t# whether to use CO terminal velocities\n",
    "flagBurtonGordon78=0\t\t\t\t# whether to use Burton & Gordon '78, A&A 63, 7 (1978)\n",
    "flagClemens85=0\t\t\t\t\t# whether to use Clemens '85, ApJ 295, 422 (1985)\n",
    "flagKnapp85=0\t\t\t\t\t# whether to use Knapp+ '85, AJ 90, 2 (1985)\n",
    "flagLuna06=0\t\t\t\t\t# whether to use Luna+ '06, ApJ 641, 938 (2006)\n",
    "flagHIIREGIONS=1\t\t\t\t# whether to use HII regions\n",
    "flagBlitz79=0\t\t\t\t\t# whether to use Blitz '79, ApJL 231, L115 (1979)\n",
    "flagFich89tab1=0\t\t\t\t# whether to use Fich+ '89, ApJ 342, 272 (1989) (Table 1)\t\n",
    "flagTurbideMoffat93=0\t\t\t\t# whether to use Turbide & Moffat '93, AJ 105, 5 (1993)\n",
    "flagBrandBlitz93=0\t\t\t\t# whether to use Brand & Blitz '93, A&A 275, 67 (1993)\n",
    "flagHou09tabA1=1\t\t\t\t# whether to use Hou+ '09, A&A 499, 473 (2009) (Table A1)\n",
    "flagGMC=0\t\t\t\t\t# whether to use giant molecular clouds\n",
    "flagHou09tabA2=0\t\t\t\t# whether to use Hou+ '09, A&A 499, 473 (2009) (Table A2)\n",
    "##\n",
    "flagOPENCLUSTERS=0\t\t\t\t# whether to use open clusters\n",
    "flagFrinchaboyMajewski08=0\t\t\t# whether to use Frinchaboy & Majewski '08, AJ 136, 118 (2008)\n",
    "flagPLANETARYNEBULAE=0\t\t\t\t# whether to use planetary nebulae\n",
    "flagDurand98=0\t\t\t\t\t# whether to use Durand+ '98, A&AS 132, 13 (1998)\n",
    "flagCEPHEIDS=1\t\t\t\t\t# whether to use classical cepheids\n",
    "flagPont94=1\t\t\t\t\t# whether to use Pont+ '94, A&A 285, 415 (1994)\n",
    "flagPont97=0\t\t\t\t\t# whether to use Pont+ '97, A&A 318, 416 (1997)\n",
    "flagCSTARS=0\t\t\t\t\t# whether to use carbon stars\n",
    "flagDemersBattinelli07=0\t\t\t# whether to use Demers & Battinelli '07, A&A 473, 143 (2007)\n",
    "flagBattinelli12=0\t\t\t\t# whether to use Battinelli+ '12, Ap 56, 68 (2013)\n",
    "###\n",
    "flagMASERS=0\t\t\t\t\t# whether to use masers\n",
    "flagReid14=0\t\t\t\t\t# whether to use Reid+ '14, ApJ 783, 130 (2014)\n",
    "flagHonma12=0\t\t\t\t\t# whether to use Honma+ '12, PASJ 64, 136 (2012)\n",
    "flagStepanishchevBobylev11=0\t\t\t# whether to use Stepanishchev & Bobylev '11, AstL 37, 4 (2011)\n",
    "flagXu13=0\t\t\t\t\t# whether to use Xu+ '13, ApJ 769, 15 (2013)\n",
    "flagBobylevBajkova13=0\t\t\t\t# whether to use Bobylev & Bajkova '13, AstL 39, 809 (2013)\n",
    "##\n",
    "flagastropy=0\t\t\t\t \t# whether to use astropy for equatorial-to-galactic conversions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputpars=(R0,V0,UsunINUSE,VsunINUSE,WsunINUSE,SYSTDISP, flagPROPERMOTIONS,flagHITERMINAL,flagFich89tab2,flagMalhotra95,flagMcClureGriffithsDickey07, flagHITHICKNESS,flagHonmaSofue97,flagCOTERMINAL,flagBurtonGordon78,flagClemens85,flagKnapp85,flagLuna06, flagHIIREGIONS,flagBlitz79,flagFich89tab1,flagTurbideMoffat93,flagBrandBlitz93,flagHou09tabA1, flagGMC,flagHou09tabA2,flagOPENCLUSTERS,flagFrinchaboyMajewski08,flagPLANETARYNEBULAE,flagDurand98,flagCEPHEIDS,flagPont94,flagPont97, flagCSTARS,flagDemersBattinelli07,flagBattinelli12, flagMASERS,flagReid14,flagHonma12,flagStepanishchevBobylev11,flagXu13,flagBobylevBajkova13,flagastropy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "checking validity of input parameters...\n",
      "printing input parameters...\n",
      " R0=  8.0  kpc\n",
      " V0=  220.0  km/s\n",
      " (Usun,Vsun,Wsun) = ( 10.3 , 15.3 , 7.7 ) km/s\n",
      " systematic dispersion =  0.0  km/s\n",
      " use HI terminal velocities?            0\n",
      " use HI thickness method?               0\n",
      " use CO terminal velocities?            0\n",
      " use HII regions?                       1\n",
      "  use Blitz 79?                           0\n",
      "  use Fich+ 89 (Table 1)?                 0\n",
      "  use Turbide & Moffat 93?                0\n",
      "  use Brand & Blitz 93?                   0\n",
      "  use Hou+ 09 (Table A1)?                 1\n",
      " use giant molecular clouds?            0\n",
      " use open clusters?                     0\n",
      " use planetary nebulae?                 0\n",
      " use cepheids?                          1\n",
      "  use Pont+ 94?                           1\n",
      "  use Pont+ 97?                           0\n",
      " use C stars?                           0\n",
      " use masers?                            0\n",
      " use astropy?                           0\n"
     ]
    }
   ],
   "source": [
    "galkin.readparsFile.CheckAndPrintParameters(inputpars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "processing HII regions...\n",
      " processing Hou+ 09 (Table A1)...\n",
      "  selected  316  out of the total sample of  815  HII regions\n",
      "processing classical cepheids...\n",
      " processing Pont+ 94...\n",
      "  selected  245  out of the total sample of  278  cepheids\n"
     ]
    }
   ],
   "source": [
    "vecout=galkin.processdata.ProcessData(inputpars)\n",
    "totallistvc=vecout[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "vecRp    = np.array([row[0] for row in totallistvc])# galactocentric distance [kpc]\n",
    "vecerrRp = np.array([row[1] for row in totallistvc])# error in galactocentric distance [kpc]\n",
    "vecvRp   = np.array([row[2] for row in totallistvc])# rotation velocity [km/s]\n",
    "vecerrvRp= np.array([row[3] for row in totallistvc])# error in rotation velocity [km/s]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "128.218523307317 386.01829276368164 3.4805916598031352 17.55446098248422\n"
     ]
    }
   ],
   "source": [
    "print min(vecvRp),max(vecvRp),min(vecRp),max(vecRp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEOCAYAAABSLcpPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XucFNWd9/HPr7tnGMABBkEgXESDUTFeIjwKuRiNbqLEqLlq4m7UXMSs2VfyyuaJRrM+PslmN5cn65rdrEqiRnc1orkou/ESb1l3VzGCNxRUCJEAQUEFGWWGme7+PX9UdVPd093Tw3RX98x836/XvOiuqu7+dVF9fnVOnXPK3B0REZE4JBodgIiIjBxKOiIiEhslHRERiY2SjoiIxEZJR0REYqOkIyIisVHSERGR2CjpiIhIbJR0REQkNqlGB9AIkyZN8tmzZzc6DBGRIWXlypWvuPvkwbzHiEw6s2fPZsWKFY0OQ0RkSDGzDYN9DzWviYhIbJR0REQkNko6IiISGyUdERGJjZKOiIjERklHRERio6QjIiKxUdIREZHYKOlIzZx5zSOcec0jjQ5DRJqYko6IiMRGSUealmpOIsOPko6IiMRGSUdERGKjpCMiIrFR0hERkdgo6YiISGyUdEREJDZKOiIiEhslHRERiY2SjoiIxEZJR0REYqOkIyIisVHSERGR2CjpiIhIbJR0REQkNk2ZdMxsgpn93MyeM7M1ZrbQzCaa2b1mtjb8tyPc1szsh2a2zsyeNrOjGx2/iIiU1pRJB7gSuNvdDwGOBNYAFwP3u/tBwP3hc4BTgIPCv/OBq+IPV0REqtF0ScfMxgPHAdcCuHuPu+8ATgduCDe7ATgjfHw6cKMHlgMTzGxazGGLiEgVmi7pAAcA24DrzewJM/uJmY0Fprj7lnCbl4Ap4ePpwMbI6zeFywqY2flmtsLMVmzbtq2O4YuISDnNmHRSwNHAVe7+DuBN9jSlAeDuDvhA3tTdl7j7fHefP3ny5JoFKyIi1WvGpLMJ2OTuj4bPf06QhF7ONZuF/24N128GZkZePyNcJiIiTabpko67vwRsNLODw0UnAquBZcA54bJzgDvCx8uAT4e92BYAr0ea4UREpImkGh1AGX8F3GRmrcB64DyCBHmrmX0W2AB8Itz2TmARsA7YFW4rTeDMax4BYOnihQ2ORESaRVMmHXd/EphfYtWJJbZ14MK6ByUiIoPWdM1rIiIyfCnpiIhIbJR0REQkNko6IiISGyUdERGJjZKOiIjERklHRERio6QjIiKxUdIREZHYKOmIiEhslHRERCQ2SjoiIhIbJR0REYmNko6IiMRGSUdERGKjpCMiIrFR0hERkdgo6YiISGyUdEREJDZKOiIiEhslHRERiY2SjoiIxEZJR0REYqOkIyIisVHSERGR2KQaHUA5ZpYEVgCb3f1UMzsAuAXYF1gJ/IW795jZKOBGYB7wKnCmu7/YoLBHpDOveaTRIYjIENHMNZ0vAWsiz78LXOHuc4DtwGfD5Z8FtofLrwi3ExGRJtSUScfMZgAfBH4SPjfgfcDPw01uAM4IH58ePidcf2K4vYiINJmmTDrAPwJfA7Lh832BHe6eDp9vAqaHj6cDGwHC9a+H28teOPOaR9RcJiJ103RJx8xOBba6+8oav+/5ZrbCzFZs27atlm8tZXR297J5RxcrN2xvdCgi0iSaLukA7wJOM7MXCToOvA+4EphgZrmODzOAzeHjzcBMgHD9eIIOBQXcfYm7z3f3+ZMnT67vNxjhzrzmERZd+RCrt3SyaXsXH7/6YSUeEQGaMOm4+9fdfYa7zwbOAh5w97OBB4GPhZudA9wRPl4WPidc/4C7e4whSwk7u9P5x1mH5ev7nAeIyAjUdEmngouAr5jZOoJrNteGy68F9g2XfwW4uEHxSWjp4oVcedY7SITdOdpaEiw4UJfZRKSJx+kAuPtvgd+Gj9cDx5TYphv4eKyBSb/m7d/BIVPb2dmd5sqz3sG8/TsaHZKINIGmTjoytLW3tdDe1qKEIyJ5Q6l5TUREhjglHRERiY2SjhTQ2JrG0KBcGSl0TWeEixZ0nd29rN7SCcDHr36YQ6a2097WwtLFC+v2ufV4bxFpXqrpSF7x2JrocxGRWlBNZ4SL1jRWbtjOx69+mKwHY2uiXZ0XXfmQuj+LyKAp6Uhe8dia7939HBBvs5uIDG9qXpMC7W0tTJ8wuqA2U89mN11AD6gDh4wUqukMM7W4QL96y07OvOaR/HtUanarZUwjqXNBNR04ckbC/hiokXSsDDdKOsPMihdfG/R7ZLKeP+uet38H8/bv4LYL3sny9a+y4MB9a3JNp7O7l53daca1pQoKWBh5t78uVZMs3iciw4WSzjCTGeD82sUF/KrNr7OrJ8Ounq6SZ90XnjBn0DGu3LA9f2YPMKY1Oej3HGqq7cAhMtzoms4wtbe1hXR2T9aq1fWb4usVxbc5iH7mSJTrwDGjYzQ3fW5B1QlH18NkKFJNZ4QrbhPP3XwNanPWnas5QXC94rYL3smCA/clYUFSA8i6j/gL6JocdWByzbO5JmAZOpR0pEB7WwtjWpNMHNu6VwmnuBNCcc1p+fpXufCEORwytZ21W9+gN+OkM87ZP1nO1HFtpLNe8jpPs6jVBezheCE8rlpXfx0v6mk4/X81ipLOMNXZ3bvXrz18+vgB/bhyZ52vvdlDd0+moCltzuSxBTWn3M3c2ttaMDMgSEo96SwvvroLgIQRa0EyVNX6bL9UImzG5DiQjherNr9OOuvMmTy2aY6nZtyncVLSaUJ7e8YYTTSrt3Sy6MqH9uqHVu2P4fDL7yGT9XzzWc6m7UEnhLaWJMmEkUstU8e1FRSO0YSUMCMb3mVcPbgKlToeipstcTCD+bMnFmwXZ8FW688qVziX6njxN7evKjmoObePnnupU4Oam4SSTpPam67PxT3XVm/pJGk1CqiMcp0Asg5dPRn21GXgxVd3FZyVb9zeRQJIpRJc/qHD+Mbtq8g6GNCTztDZ3bvXiade0/Y0y7WE4mZLI/grVu1Z9aIrH2LdtjeZM3lsDaMcuGrinbd/B20tSdJZ58BJY/mb21flT14+etXDjEolaE0l2J3O5l+jE5nmMeCkY2YXuft36xGMBJYuXsjsi39dk/caaBfqUsoVBKsu/wArN2zno1c93Oc1bS0JEmaks05P5Me/fP2r+cJ61+40WWBUKsEdT26mrSXJ2FFJtnX2sDX8mzutfcBxlmrzv/NLxw3qO595zSM1nQ6oOHlVeo9S647//oMFzZG5WmXxtrnvUan2HP1exTXkXPPUoisfGvA+rKdkwtjdm2Hj9q6C5AKwO50lmTBSCaMnXNZMXdGb5cSlUfpNOmZ2a/QpcBSgpDPCRQv4UqaOa2PKuLaCAi16TWfp4oXMvexudvVk6OxOs2LDdlIJ4/Wuwi7ar7yxG9hTc0klgkRWqQAZzLQ9lQrnwQziPPOaR1i9ZSdAQZPkR696mDGtQcIAmDttXL/v1dndm084AGaWf7/jv/8gr74ZFLVvdgdJPUHQ9FZOcWX1ha1vMLolWRDn6i2dHPyNu3AY9PWRSgk2WiBXMnfaOFZv2cncaeMKjrGciWNbmT5hNJ3dvXzwiLeUHdQcd5fzRnaCyGl082I1NZ2d7v653BMzu6qO8Ujo2AMm9r9RkbUvd/LarsIkMKNjNNMnjC5Yliv8cj9cCArC3AXXYp3dvbzyxu587euMo97Cqs2v09WTIVmm/e6VN3bne6LNndbOum1vcuCksfkf/soN2wuuBWWyTqZEU93Wzh4OuPjXFK8pVdOITtvzsasexoHWVHCGmxNtRip15l5upoSlixcWXEsw4C+Pn8Onjp1V8vtH5RLOrt1B0ir+ml09GcKck///yMm9BsrXWqP7LZqMcrJAnx1YQTrjdGb6JupcjWL1ls6CRDlQ5Qr64gK5JZnAoey1yVyCbm9rYfa+Y/LfvZlqNcU0+wSYe+Wj0cwOcPc/RJ5PdPfBz7XSQPPnz/cVK1Y0Ooy6iDbLtbUkSg42zNUaxrWl+pwhAvkCpbM7TYKw0KqB/dpbGdOaYldPmte70n2aRQZqcnsrb+7OMHFsK6mEsasnjQM7dvXSGymhW5JGSzLRp8PD1HGjSCUTjGtLsXF7V59OEbnXJROWT9BvdKfz5Xe0l100kUctXbywTyEbLVxz77Fxe1fBNm+Gn2PsqaXUoqm0HpIGY0YN/vLw7nS2oCk2qpokl8l6/kTobfvtU9fCfG9rC8WdIAYyGLgZmNlKd58/mPfo90jJJRwzm+Turwz1hNMoA+0mGS2oyhVoxYqbutzh3Ot/V7CsVG+zYl09mXzBWquEA0GthXwr++Bt6wzea1dPV8XtejNOb6bvd35p5+7841GpRJ9OEbnXpZJGZ3cvmawXVBiyDn96vZvs9i560llaU6Un+Ojs7mXNlk6SSWPGhNGks87sfcfkmwmD3lZ7vkMm6/n97sCYlsICd+60cTz2h9fK/t/kOhUUr6+0fG/yWa6SO2ZUqqpmwf5Ek3FxTLnmsmLFNUMnqKkV18ZqEV8tFN8+ZCglnFoZyOnJdcBp9QpECkV/TLkmluIfWLGu3sKCdXc6u1e1iWTSSDfraXWdVNpPuUKslB2R5szd6SyP/qHwnCzaNJjOeEHzV4LgpCD3/+seFJqpoibL4gJ36eKFHH75PcCeZtFcDWH2vmOYMq4NoE8so8NBv5sitSoDDp0W1NZKjWnJHXMzO0b3mS9v9TdPLrlPBuPgb9xFOut8/t0HsOS/1vc7H1305Gzzjq6CE5Dofmv0dYyokT77xECSTp073w5v0Quk37v7ufxF0FI/hpsf/SNGUIjlzp5zOaD4jC16DeKVN3aHtYnAfu2tvPpGT/61SQsKtv7S0EhLOPVUaU9mgc4SnRx6i/b/uLbCn2m0oO0Oa6VjWpPs7s3kE05xrXdUKmgmHNeWyk9BZMC3P3x4/rpUf7XxlRu2s3z9q/z66T9VbLoazDizXPJf8l/r89d0Dpw0Nj/2ppLofmrm6zoj3UCSTiwlkZnNBG4EpoSfucTdrzSzicBSYDbwIvAJd99uwbD2K4FFwC7gXHd/PI5Yq7Vyw3bWbOnEgU/+eHmfi/W5H/OCA/fl+Zc6ueRXq/LronOUzewYXXF24stOPazgtdEElKN00pxKNXHlmsOC6z19mxCjzXDdvRnaWpJ9ajc5uVrv8+Egyb1p3snd5uKhF7ZV/ZqquGNkebNrN4RXsoKkmGVMCiaOypLIdpEkQ5IMCc/Q4r200EvKe2nxXhJk2Z0cxdTJxrxDD2LBzDbmjX0Vtm2DsZP6/fwKKyu8psLrKrzn0k8d2H9Mw1i/HQnyG5r9u7t/qM7xYGbTgGnu/riZtQMrgTOAc4HX3P07ZnYx0OHuF5nZIuCvCJLOscCV7n5spc+oR0eCSmd3f3jljYIE0JI0sllnzKhUn2aL/oxKJThq5gQAnti4o+yF10bb2+sEtVeqKC+3rlq5b5e71B9dFvybIEOWZNFr9nxukiy3tn6TeYm1rMwexC8z7+G2zHFkSNJCmpta/455ibUVo/hR+jR+kP44WZIkyfCV1G1cmFoGwKW9n+GmzIn5eN6fWMGS1iv28vsWWpk9iOXZQ1mQWNNvjAN937N7LqGXVNX7YMg64Dg4598bHcWAxdKRIOLrg/mgarn7FmBL+LjTzNYA04HTgePDzW4AfgtcFC6/0YPsudzMJpjZtPB9mkJx0RZ9PtBbB5S6blA/QaRGhkPYyBpmUfpuGMUpxjAynJRYyQPZo8mEr7FwGy97R41owU3R48KY9mxPmW2Ln5dbVxi7kesskCja1vtsW+49APZlJ9voKFrn+e0dWJ49FICzey6hhxQpspyZfICPJv+7qsJ2QWINraTpBVpIsyCxJr/uI8n/YmnmeNJh4nsw+w4u6T2v6vcuJ5cYekjRWuPEMC+xlstSN3JX9hhOSfyutgmnbTx9jiUrd3wUKTvIySoPgKr0nu/5aoXXDW9V13QawcxmAw8Bbwf+6O4TwuUGbHf3CWb2H8B33P2/w3X3Axe5e9mqTNxdpgvGjSSNOfvtw8btXcydNo6vnXxIwYj+/rsoO+9KPMMPWq5mqm0vOOv8bM9fs4N2+hbawb9JMpyYWMl92Xlk8403hpHFK/xAkmS4tfVbAPwi825e8fFs9334vb+F14j+mHNJKsuosEAC8vHlHnfQyW+zR/FsdhZvMpo3aCNNC6UK71KxHM7veZK3VdiqmnpWdJvipJLvrFzifUolnr7P59gm1vmMovcmv20C5+DWrbya3Yet6bH57ZIEnUGyJDCMsW0tFXte5cZL5bpW5y4HJSLfpFipwaLRHmiVektu3tFV0BGh0jiwgSruWbk3Y4FKxdxMnQiGurhrOpjZQ8Cp7r7TzC4A2oB/cffa9YPd81n7AL8Avhx+Xn6du7uZDShbmtn5wPkAs2b1P6Cvlubt38ERU0cxc+dKzpv1CqP/9AhzWRXU566HX7T2ba64OX0CSzMn8Iq3s5n92FMoZflK6udMtWDE9rzE2vxrjk88ye3Z91C6kMuQJMtvs0fhQJJsmOCMFtKcm7yb+zJHg0GXt0Y+EzIkWJ49lAtTy/LNQWf3XMLu/OGzp/A2srw78QxfTv2SeYm1bEzN4sKT3wfTvwrt05g3eiIkUxxcZvocMI6Z3cFBU9q5f83LdHanybqTdXjv2yaz+L1v5Uu3tEOk4GtJBpOFZrKF79O/Ume6xbOYlXqf4mV9n6/zmf1+zgu9U4sGihqZyE/SCcbrlJuHb8yoFMmEsU9bil2704xqSeYL7dz0QqV65aVSCUaV6d7dn+jF+oQVPo9OGbQ3KvU+q5YSTPMb6Iiu8WECmAd8HvgP4MfAObUMysxaCBLOTe7+y3Dxy7lms/C6z9Zw+WYg+gufES4r4O5LgCUQ1HRqGW81zu29lQ9nl8KLQRPFpZnP4Hi+uaO4KeFTqQf5VOpBVmYP4pM9l9JDiiRZvpW6vmyzwz+2Xs0T3XPYwLT8soNsE1PtNTb6FHYkxrMj00ZQuIEljWzGSbS2ciMfYVcmU/LUOJVIsODzP+TMu/+S1Vt2Mn5sC91lxsaYJXh1yrv5XtsJLF28sOA/Jqr47qE5CYP3HrwfF54wh3Vb3wjO5HuDgvM3q1/mgee38s3T3p6fHDRhMH3C6JIj8fMxhf/mvlpujEw6ky0Yq1NL7W2pgoGkUfnabBVHYRao5mbeY0alGD+6paDQntw+itfe7KE3k833iEtY/9PY5Lphl5pwNXe/pUq3C6jVwEn1PhueBpp0es0sBXwa+K6732pmNW2nCpvOrgXWuPs/RFYtI0hu3wn/vSOy/ItmdgtBR4LX63k9JzqH1kCkdy/i8cQYtiSncn/PYfnmrJszJ4XNCGXOPBPQQziOgwR/n/oCf88XCjZZdfkH8o+3XXY3hGe7rakk3/n8mXzv7ueYAvzg5EP6jIbOdUUtPsucMKaF3nSWd86ZxAXvfWv+x59rEoy+T8KM7p4MqVSi6nm5cnOwFWtNJQrWFQ/YTGecb9y+ilkTx+Sn2envuthJc6dwwsH7ccV9z2PADz5xVP77lJpiZ7DGtCaZO21cwe0HorLsaTqqNFg3N1tBtRNtliq0o9P8l5rep1i52wFEJRNGMmE1v02ABk6ODANNOv8EPEXQrHZxuGyfmkYE7wL+AlhlZk+Gyy4hSDa3mtlngQ3AJ8J1dxL0XFtH0GX6vBrHUxM9tHBz9s9IWgIvumqTzkKyitaOitcsQ8mEkQDe0jE6/8ONFgzFP+rofGXRAuvac/5X2R99ceGQG3cEVD31yLz9OwrmyzKCaW2u+vP5Ze+5k5P1IBn990Xvy8d+9k+W05vO0pJK0N6Wys9WAHDv6pe5b/XL+eSSu232vP07+PaHDy/oZh6VACwR9DR0gkS8Y1ffCU5z/y2T21uZtM8oNm7vYvWWnXRVmPmhqyfDPm1B81hxM5gBk9pbmbzPKNrbWgruxFpJpUI7NyCxP42eG2ykD5wcCQaUdNz9BjP7FdDr7l1mNgeo6TStYYeAcsXriSW2d+DCWsZQyWDP7KKdCiDoWPCzz9d2/qWxbal8gVys3I+6mrPM6ADXWkxzP2VcWz7p/PwL7yz5me1tLfziC+/kl49v4uZH/4iz5yw+GvtNn1uQH+sE5BNoTvH0NblbLOQGRpZKPK0tCaaOa+PNnjTbOnvyCSeZCJrPWpMJvnzSwWzf1ZMfMBmda63U7Mew51YE0Yvej/9xe74JzIFRqeReFfbF/7/F11n6O36raeIa6Xe+lMEZaEeCHwMfBDJm9ifgaaCpBmI2u3n7d3DotHbWbn2DCWNauLrozH6wok1tpVQqKEolpGrvUTOYua2SRsV9kBuU+JGjZ+QTS6mkGV122wVBorptxUYyWSeZMHoznk9a0Sa8Tx07i3+87/mCsVRHzhjPZR86jO/d/RybdxRev5o2fnSfpH7hCXPyj6P7OHrfm3K1OQgmYY3eAqLWzUvVJohqTj7qmWyUyIa/gTavHQfMdPeMmU0HjgSOqH1Yw1t7WwttLUkOnLRPUzUjVPrBD+QeNQO5SVUCSCYTJbctfp/ixFJJqUQVvaVx8ftc9efz+eSP9zTRXfahw/LNj6XO/qs1ZVwbWzt3M3Fsa/56SrkanQGJhHHZqYft9XFRi0JbTVxSTwNNOo8C+wJb3X0zQS+xO2selTSNctd8ShW8q7fsLLgw3t9dOzu7e8kCPelswXWW3OflzvzP/snyvZ4CPpqoKhWm8/bv4GefX1CyJtXf2X9/zU2HTx/f5xYHpW534AQ9x75x+yr+bfmLfe7nIzIcDDTpXAP8p5ldS5CAnnb312sfljSbappdMlkveV/6coprT9FbWUe7VPemswXr9lZ/BXelmlSlhFWpZhf9zIHWJBt1cy8lOKmngSadfyMY65IC/hI4wsza3P2tNY9Mmk6lgrdU9+D+mqKKZwWOXmdZcOC++clOW4q6Ucel1D2NKt2QrdrbDy9dvLBPwb4yMlhWY1RkOBto0tnk7n8fXWBmo2oYjwxRxeNpkgnrt0ms0nWMZhuzUa6jRK1qKPP272DutOb5viL1MtCk86SZfcndr8wtcPf6DOmWIWXO5LE891JnfpaAg6f030mi+DrGwVPbS44taVQBXE0zUy1H0Tf6+4rEYaBJZwpwkpldRNBV+ingSXe/reaRDWPDsc28va0lXzPpb9R7TqVrOkNFLWtkw/G4ECk20MGhn4B8k9phwOEEU88o6YwA5QrF6ADEgTQtVbqmM5TUYqCsyEgxoKlmzexvIWhSc/fH3f0GgnvaiAzYnV86jrnT2pnRMXqvu0SLyNAy0PnNp5vZp3JPzGwycF9tQ5KRpL2thekTRivhiIwQA72msxi4x8zWEUwRdT2q6YiISJWqSjpmdiNBx4EnCCbXvBlIA2e4+7r6hSdDkS6Ii0g51Tav/ZRgvsLzCAaIzga2A39uZh+rS2QiIjLsVFXTcfcHgAdyz8MbuR1KMOHnscDP6xKdiIgMK9U2ry0Elof3rsHd08Cq8O/f6heeiIgMJ9U2r30aeNzMbjGzc81saj2DEhGR4ana5rUvAJjZIcApwE/NbDzwIHA38D/uXv7evCIiIgx8RoLngOeAK8xsNHAC8HHgH4D5tQ9Phjv1dBMZWQY6TifP3bvM7EB3/ysAM5vg7jtqF5qIiAw3A52RoNj+kceXDPK9RERkmKsq6ZjZoeVeb2bvMbMEMLF2YYmIyHBUbU3n12Z2vZnNKlr+v4EjgB8Dd9Q0MhERGXaqvaZzCMG8a/9pZsuAv3X3be6eBX5Ut+hERGRYqaqm4+497v5PBLMQbAR+Z2bfMrPS9/AVEREpYUAdCdy9293/H/B2oAtYaWZfrUtkIiIy7Az0Jm6zzexk4HPALKAT+Lt6BDYQZnaymT1vZuvM7OJGxyMiIqVVO/fa08B04I8Eg0PXAPcD/wy8ULfoqmBmSYLrSn8GbAIeM7Nl7r66kXGJiEhf1XYkOAP4Q27CzyZzDLDO3dcDmNktwOmAko6ISJOptiPB+iZNOBDUwDZGnm8KlxUws/PNbIWZrdi2bVtswYmIyB6DnZFgyHD3Je4+393nT548udHhiIiMSMMh6WwGZkaezwiXiYhIkxkOSecx4CAzO8DMWoGzgGUNjklERErY61mmm4W7p83si8A9QBK4zt2fbXBYIiJSwpBPOgDufidwZ6PjEBGRyoZD85qIiAwRSjoiIhIbJR0REYmNko6IiMRGSUdERGKjpCMiIrFR0hERkdgo6YiISGyUdEREJDZKOiIiEhslHRERiY2SjoiIxEZJR0REYqOkIyIisVHSERGR2CjpiIhIbJR0REQkNko6IiISGyUdERGJjZKOiIjERklHRERio6QjIiKxUdIREZHYKOmIiEhslHRERCQ2TZV0zOz7ZvacmT1tZr8yswmRdV83s3Vm9ryZfSCy/ORw2Tozu7gxkYuISDWaKukA9wJvd/cjgBeArwOY2VzgLOAw4GTgX8wsaWZJ4EfAKcBc4JPhtiIi0oSaKum4+2/cPR0+XQ7MCB+fDtzi7rvd/Q/AOuCY8G+du6939x7glnBbERFpQk2VdIp8BrgrfDwd2BhZtylcVm65iIg0oVTcH2hm9wFTS6y61N3vCLe5FEgDN9Xwc88HzgeYNWtWrd5WREQGIPak4+4nVVpvZucCpwInuruHizcDMyObzQiXUWF58ecuAZYAzJ8/30ttIyIi9dVUzWtmdjLwNeA0d98VWbUMOMvMRpnZAcBBwO+Ax4CDzOwAM2sl6GywLO64RUSkOrHXdPrxz8Ao4F4zA1ju7he4+7NmdiuwmqDZ7UJ3zwCY2ReBe4AkcJ27P9uY0EVEpD9NlXTcfU6Fdd8Gvl1i+Z3AnfWMS0REaqOpmtdERGR4U9IREZHYKOmIiEhslHRERCQ2SjoiIhIbJR0REYmNko6IiMRGSUdERGKjpCMiIrFR0hERkdgo6YiISGyUdEREJDZKOiIiEhslHRERiY2SjoiIxEZJR0REYqOkIyIisVHSERGR2CjpiIhIbJR0REQkNko6IiISGyUdERGJjZKOiIhOqObaAAAIoUlEQVTERklHRERio6QjIiKxacqkY2Z/bWZuZpPC52ZmPzSzdWb2tJkdHdn2HDNbG/6d07ioRUSkP6lGB1DMzGYC7wf+GFl8CnBQ+HcscBVwrJlNBP4PMB9wYKWZLXP37fFGLSIi1WjGms4VwNcIkkjO6cCNHlgOTDCzacAHgHvd/bUw0dwLnBx7xCIiUpWmSjpmdjqw2d2fKlo1HdgYeb4pXFZuuYiINKHYm9fM7D5gaolVlwKXEDSt1eNzzwfOB5g1a1Y9PkJqbOnihY0OQURqLPak4+4nlVpuZocDBwBPmRnADOBxMzsG2AzMjGw+I1y2GTi+aPlvy3zuEmAJwPz5873UNiIiUl9N07zm7qvcfT93n+3uswmayo5295eAZcCnw15sC4DX3X0LcA/wfjPrMLMOglrSPY36DiIiUlnT9V4r405gEbAO2AWcB+Dur5nZt4DHwu2+6e6vNSZEERHpT9MmnbC2k3vswIVltrsOuC6msEREZBCapnlNRESGv6at6cjQo95mItIf1XRERCQ2SjoiIhIbJR0REYmNko6IiMRGSUdERGKjpCMiIrFR0hERkdgo6YiISGyUdEREJDYWTGs2spjZNmDDAF4yCXilTuHUk+KOl+KOl+KO1yRgrLtPHsybjMikM1BmtsLd5zc6joFS3PFS3PFS3PGqVdxqXhMRkdgo6YiISGyUdKqzpNEB7CXFHS/FHS/FHa+axK1rOiIiEhvVdEREJDZKOiEzm2lmD5rZajN71sy+VGKb483sdTN7Mvy7rBGxFjOzF81sVRjTihLrzcx+aGbrzOxpMzu6EXEWxXRwZD8+aWY7zezLRds0xf42s+vMbKuZPRNZNtHM7jWzteG/HWVee064zVozOye+qMvG/X0zey48Dn5lZhPKvLbiMVVPZeK+3Mw2R46FRWVee7KZPR8e6xfHF3XZuJdGYn7RzJ4s89pG7u+SZV/djnF311/QxDgNODp83A68AMwt2uZ44D8aHWuJ2F8EJlVYvwi4CzBgAfBoo2Muii8JvATs34z7GzgOOBp4JrLse8DF4eOLge+WeN1EYH34b0f4uKPBcb8fSIWPv1sq7mqOqQbEfTnw1SqOo98DBwKtwFPFv+G44y5a/wPgsibc3yXLvnod46rphNx9i7s/Hj7uBNYA0xsbVc2cDtzogeXABDOb1uigIk4Efu/uAxmwGxt3fwh4rWjx6cAN4eMbgDNKvPQDwL3u/pq7bwfuBU6uW6BFSsXt7r9x93T4dDkwI654qlVmf1fjGGCdu6939x7gFoL/p1hUitvMDPgE8LO44qlWhbKvLse4kk4JZjYbeAfwaInVC83sKTO7y8wOizWw8hz4jZmtNLPzS6yfDmyMPN9EcyXUsyj/Y2zG/Q0wxd23hI9fAqaU2KbZ9/tnCGrApfR3TDXCF8NmwevKNPU08/5+D/Cyu68ts74p9ndR2VeXY1xJp4iZ7QP8Aviyu+8sWv04QRPQkcA/AbfHHV8Z73b3o4FTgAvN7LhGB1QtM2sFTgNuK7G6Wfd3AQ/aGYZUN1AzuxRIAzeV2aTZjqmrgLcCRwFbCJqqhpJPUrmW0/D9Xansq+UxrqQTYWYtBDv9Jnf/ZfF6d9/p7m+Ej+8EWsxsUsxh9uHum8N/twK/ImhmiNoMzIw8nxEuawanAI+7+8vFK5p1f4dezjVRhv9uLbFNU+53MzsXOBU4OyxM+qjimIqVu7/s7hl3zwI/LhNPs+7vFPARYGm5bRq9v8uUfXU5xpV0QmGb67XAGnf/hzLbTA23w8yOIdh/r8YXZcmYxppZe+4xwYXiZ4o2WwZ8OuzFtgB4PVJtbrSyZ4DNuL8jlgG5njrnAHeU2OYe4P1m1hE2B70/XNYwZnYy8DXgNHffVWabao6pWBVdg/wwpeN5DDjIzA4Ia9BnEfw/NdpJwHPuvqnUykbv7wplX32O8Ub0lmjGP+DdBNXHp4Enw79FwAXABeE2XwSeJegVsxx4ZxPEfWAYz1NhbJeGy6NxG/Ajgp49q4D5jY47jGssQRIZH1nWdPubICluAXoJ2qw/C+wL3A+sBe4DJobbzgd+EnntZ4B14d95TRD3OoI2+NwxfnW47VuAOysdUw2O+1/DY/dpgsJwWnHc4fNFBL2vft8McYfLf5o7piPbNtP+Llf21eUY14wEIiISGzWviYhIbJR0REQkNko6IiISGyUdERGJjZKOiIjERklHRERio6QjIiKxUdIRqSMzW2xmL4WTlv7ezD4dWXe1mb0rfHyimf1r4yIViYeSjkh9HQ5c7sGkpZ8EotOMLCCYaQHgSOCJmGMTiZ2Sjkh9HQE8Hz7+A9ADYGaHAi+4eyZcdyTwhJmNMrOfmtnfhXPl/Sy8++TvzGyDmX0wfP1bzOwXZvaEBXcCbeiEnCLVUtIRqa/DgefDSRW/CFwaLj8FuDuy3REEs/jeA9zn7pd4MEfVkcB6dz8GOBv4P+GsxXcB17v7OwjuVrkmlm8jMkiae02kTsxsJkHt5hmCG1s9DZzk7m5m9xBMjvincFr5V4ANwGJ3fyR8fRvB5Jwz3b3bzCYS3Fzr68DH3P2s+L+VyOCopiNSP4cDD7n7UcDbgEMI7oQ6Bpjg7n8KtzuUYFr+NJCJvP7twFp37w6fH00wE/FR7LkWJDKkKOmI1M8RhJ0DPLh//M3AB4ETgAcj2x0JPExw/5frzWxKZPksM2sL77Pyf4ErCG4dnL91t5lNrvP3EKkZJR2R+jmcwh5p/05wn5Li6zlHAs+4+wvARcCtYZPbkcAvCZrUHgOucvf/Ibg/yxQze9bMngQW1vuLiNSKrumIxMzMHgeOdffefrb7T+B8d3++0nYiQ4mSjkiTMrNNwCx3zzY6FpFaUdIREZHY6JqOiIjERklHRERio6QjIiKxUdIREZHYKOmIiEhslHRERCQ2SjoiIhIbJR0REYnN/wfu3bVvuy2V0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(vecRp, vecvRp, xerr=vecerrRp, yerr=vecerrvRp, fmt='.')\n",
    "plt.plot(vecRp,Vc_Model_I(vecRp,params))\n",
    "plt.xlabel(r'$R/kpc$')\n",
    "plt.ylabel(r'$V_{c}/km s^{-1}$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEOCAYAAACaQSCZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XucXGWZ4PHfU5fuTiede0xCyIUsKGQmI5AWg7dlZXUQZ8UZxYCMgDICM7irO7PjMOAl4+iuoqMbd1xIRkQYEeJlRFYZGAZxs0qCJtwCSZQYkgltQ0LSSTrp7nRX1bN/nPNWnzp9qvtUdV27n+/n00nVqVPVb1dXv895b88rqooxxhgzlkS9C2CMMaY5WMAwxhgTiwUMY4wxsVjAMMYYE4sFDGOMMbFYwDDGGBOLBQxjjDGxWMAwxhgTiwUMY4wxsaTqXYBKmjt3ri5btqzexTBl2HPwRMH95fOmFhx398t5zXKea8xksm3btldUdd5Y502ogLFs2TK2bt1a72KYMqxZv7ng/sbrzi847u6X85rlPNeYyURE9sU5z7qkjDHGxFLzgCEiSRF5UkR+5N8/TUQeF5HdIrJRRFr8463+/d3+48tqXVZjjDHD6tHC+CiwM3D/C8BXVPV0oAe4xj9+DdDjH/+Kf54xxpg6qWnAEJFTgXcCX/fvC/BW4Hv+KXcC7/ZvX+Lfx3/8Qv980+DWrN88YkzCGNP8at3C+J/Ax4Gcf38OcERVM/79F4FF/u1FwH4A//Gj/vnGGGPqoGYBQ0T+ADigqtsq/LrXishWEdl68ODBSr60McaYgFq2MN4IvEtE9gL34nVFrQNmioib3nsq0OXf7gIWA/iPzwAOhV9UVTeoaqeqds6bN+Y0YmOMMWWqWcBQ1b9W1VNVdRlwGfATVb0CeBR4r3/aVcAP/dv3+/fxH/+J2n6yxhhTN42wDuOvgD8Xkd14YxS3+8dvB+b4x/8cuLFO5TPGGEOdVnqr6k+Bn/q39wDnRZwzAFxa04IZY4wpqhFaGMYU1TswRNeRfrbt66l3UYyZ9CZULikzcaxZv5negSF2dPcCcOltj/Hd69/AqqWz6lwyYyYva2GYhnVsIJO/nVPYsmfEJDljTA1ZwDANaeN157PusnNI+Gv729IJVi+3dZvG1JN1SZmGtWrpLM5c0MGxgQzrLjvHuqOMqTMLGKahdbSl6WhLW7AwpgFYl5QxxphYLGAYY4yJxQKGMcaYWCxgGGOMicUChjHGmFgsYBhjjInFAoYxxphYLGAYY4yJxQKGMcaYWCxgGGOMicUChjHGmFhqFjBEpE1EfiEiT4vIcyLyN/7xb4rICyLylP91tn9cROSrIrJbRJ4RkXNrVVZjjDEj1TL54Engrap6XETSwM9E5J/9x/5SVb8XOv8dwBn+1+uBW/3/jTHG1EHNWhjqOe7fTftfOspTLgHu8p+3BZgpIgurXU5jjDHRajqGISJJEXkKOAA8rKqP+w99zu92+oqItPrHFgH7A09/0T9mjMHbxnbN+s31LoaZRGoaMFQ1q6pnA6cC54nI7wJ/DZwJvA6YDfxVKa8pIteKyFYR2Xrw4MGKl9kYY4ynLrOkVPUI8Chwkap2+91OJ4E7gPP807qAxYGnneofC7/WBlXtVNXOefPmVbvopoHYFbYxtVXLWVLzRGSmf3sK8DZglxuXEBEB3g086z/lfuBKf7bUauCoqnbXqrzGGGMK1XKW1ELgThFJ4gWq76jqj0TkJyIyDxDgKeB6//wHgIuB3UAf8MEaltUYY0xIzQKGqj4DnBNx/K1FzlfghmqXy1Re78AQxwYybNvXY3txGzOB1LKFYSYwN5bQOzDEju5eAC697THOXNBBR1sagI3XnV+38hljxs9Sg5iKOjaQyd/OaeF9Y0xzsxaGqQjXeti2r4dLb3uMnEJbOsG6y86xbiljJggLGKaiVi2dxZkLOjg2kLFgYcwEYwHDVFxHW5qOtvSIYOHGOWwsw5jmZGMYxhhjYrGAYYwxJhYLGMYYY2KxgGFqpndgiK4j/Wzb11PvohhjymCD3qbq1qzfXHRB33gHwHd0H2PN+s02kG5MDVgLw9SELegzpvlZC8MUValpsBuvO98W9BkzAVjAMFUTDDjjXdBnaziMqT8LGKZmii3oqxQLKsaxz0J12BiGMcaYWCxgGNOkbJqyqTXrkjJFbd17uN5FML7w3uW274ipBwsYpmZqVZGFK9eJKGqasgsYxlRLzQKGiLQBm4BW//t+T1U/LSKnAfcCc4BtwAdUdVBEWoG7gFXAIWCNqu6tVXmNaSThYGvTlE091HIM4yTwVlV9LXA2cJGIrAa+AHxFVU8HeoBr/POvAXr841/xzzM1lFXvq1zV7mPvHRjiZCZH78BQVV6/kblpyqfOmsLdf7LagoWpiZq1MFRVgeP+3bT/pcBbgff7x+8E1gK3Apf4twG+B/y9iIj/OqbBBLuBdnQfI5tT+gazQHX62C9etynfh7+ju3dSDvxWe5qyMWE1nSUlIkkReQo4ADwM/AY4oqquQ/ZFYJF/exGwH8B//Chet1X4Na8Vka0isvXgwYPV/hFMTJnccFyvRiqQ8Ott2XOooq8/Ea1Zv3lSjO+Y6qnpoLeqZoGzRWQm8APgzAq85gZgA0BnZ6e1Puok2GpwyQZ3vdRbtT72dZedk+/DB1i9fA6bfm0XDMZUU13WYajqEeBR4Hxgpoi4wHUq0OXf7gIWA/iPz8Ab/DZNoKMtXdU+dteHL0B7S9K6ZYypgZoFDBGZ57csEJEpwNuAnXiB473+aVcBP/Rv3+/fx3/8JzZ+0Vw62tIsmjmlapV5R1uahEAyIVV5fWNMoVp2SS0E7hSRJF6g+o6q/khEdgD3ishngSeB2/3zbwf+UUR2A4eBy2pYVmOMMSG1nCX1DHBOxPE9wHkRxweAS2tQNGOMMTFYLinT8HZ0H2N711HLm2RMnVnAMA3vxECGvsEsL/b0c+ltj03aoGHTYuOzxIzVYbmkTMMJ79MdnOmQU2/NRbPNirL9GYZVO+iNlpixWibL79VaGKbhTWlJ5m+3pROsXj68flMVTmZydiVp8mz/+OqxFoapmkpddSUTQntLktlTWwoWAPYODJEDBjM5rvj6FpbPnWoZW5tAta/GLTFj9VjAME0hmZARazqCV45DmZyl+K6RqO61RupyG+/+8fXQSO/faCxgmIayo/sYfSfjdSFMbxv++KZTiYL7ZnTNUkGVK05ixon+HlSDjWGYitt43fnj+iPMKjy5/8iY4xIdbWkSQEsqwd1/shqAriP9kzLd+VjWrN/Mju5j9S5GxdnMsdqySzIzpm37eqrWrHd/7Nu7jpLJKUOZHOCNS1x622O0pZOcHPLSpEdVeCLQmkr4j3szYxJCTWbGNJtsTuk60s/F6zaxv6efFQun17tINWOticqwgDEJlHMFFrxKL3dqYpw/zvDeGUE5LUyTPppgevPJvmVp1O97e9dR+gaz9A32A17Xgpu+7Ey2yrR3YIhjA5mqXhBNNBYwTKRa7Rm9YuF0uo705yuyoLZ0guVzp7K/pz9/bjGrl88hIV5ZE4KNZ4SEA2+xMDzWlbjr2mrW1on7+UZbqzHZAmcp7K9qEijnD2Dbvh7ec+tjQHWnJm687vyCaZCOG5e45cFdZHN9ZHI6YmxiR/ex/BaywZkx09tSk7Z1AdG/7+AOhQmBtnSSFQunT6rK0XVprlg4vWYXRBONBQwTKRgcrj5/WVWb7MHK/vCJQfoGs2SzOVYtnZXvSoHhrViLlcXNjBmv8fR3N2pfeUdbOr+WZXpbKt9qazblvr9u//fT500dcZFiazXis4AxCaxc+1DJzxn0B58Bbtu0hzse20tLqrRJddvX/n7sc10FFt7bItyVUiwtSLCrpHdgiN0HT5BKCCsXzSipzM2k1EozmRCO9g9xtL/8WWSu4m30mWiuNeF2fwzu/37xuk35Db6aaa1GI7CAYSKFK+pMTmmJ8Ty3hqK9tTIfrdPnTS3oSgmmBYlarxGsHAah4Su2RhQeDHeKVbzuOUBBDrBGEU4N4rqf4qzVMIUsYEwScRfDOdnQqGg2p7Fewz0v6txyuhNcV0omp5w+b+qYf9y7D54ouG95hCqnWMXbiNygfLj7SYA/u+B03v/6JbZ+owwWMEzDKBaQkgkhmZBYlVMqIQwG7pczW6qc6ZZjVT7NNIWz2GB4sOJNCAVdOcGLgUYZxwm2lNLJBCczORT4xH3b+daWvXS0petexmZTs4AhIouBu4D5eLP6NqjqOhFZC3wYOOifepOqPuA/56+Ba4As8F9UtfTOeMP2tb9f8tXU8y/3crhvuDtndnuaM+Z3jPm84EyU8SilS6vvZCb/fV2iwkxOSQWCTPjnj6oowv3dbrrlAx99S7k/RlVesxrifD7c5ITdB0+MaO0FA2IjCqfIb7TWUbNcUNSyhZEB/kJVnxCRDmCbiDzsP/YVVf1S8GQRWYG3j/fvAKcA/yoir1bVkSu8zJhKTQvRP1T4Nh8LVMqjca2EUr/fiYEMOUrLVdPemqI31E0y3jGUaqTGbvR021FpQ4oFkP09/WSzOfb39Bdd09CWTpJMyLi6fCpx5R9sKQWnFcedFVWrLqt67N/hlPo+13JP726g27/dKyI7gUWjPOUS4F5VPQm8ICK78fb+to7HMpR6xf/CK8c50DvcuTO7Pc1pc6cVnBPVmiilhRGscNycLPf/iVClGtVN4oJDVr3bwSAR7t6KU65i0y3HI/yariunHOO5CnWzm1L+LLShbI49rxyP1SIY670LB8ShbI5MTugdGKrbVXz4sxKcVtxos6KaaU1IXcYwRGQZcA7wOPBG4CMiciWwFa8V0oMXTLYEnvYiowcYM4pSrySCC/daUglu/ePOEX9k4b5qd3+0BWFb9x5m5dqH8lNud3Qf42RgCq+T88sw2vfTUbKGuMH37V1HS7pSrEZq7NG6ckazZv1mtu49DHg/q3uX3nPrYyTw8mhB8daUC5rB5wbHdw72DuZfK/wbcEGi72SGrEJShjMJZxWOD2TY3nWUZELIhmbUDWUVUHa91Dvq1XI5rYhyg2ZwbCXu82o1vtFMa0JqHjBEZBrwfeBjqnpMRG4F/havm/Fvgb8DPlTC610LXAuwZMmSyhd4kgp+YO/58OqKfYCzSr4baeN157Ny7UMjBqqdYH6osG37egoquaySr7iCgaRvMMv2rqOcHMrG7qYqNt0yuNaj1Mqkoy1NayoR+8ox3E0Ujo2KN+MnjrGycUU9Hmy9uf/7Tmby763ivbfBwBU22tVysam7o4nqunHrd+K8VqOmNGmm/TtqGjBEJI0XLO5W1X8CUNWXA4//A/Aj/24XsDjw9FP9YwVUdQOwAaCzszNepjpTkmp/gMOL9ZwfP/PbyBXJa9ZvpuvIyONRCQwB+gezCN5Wri8fGyCT06JXqOXmSqrUYH9QMDCFr0Lv/pP4QTwq9YrTkhTuufZ8Pnnf9nxalWAmW9eqaG9N5fN+vRj4nZwyawqLZk7J//yLZ02JNVZQzvjAeMeC+mKOw9VDs6wJGTNgiMjsGK+TU9UjY7yOALcDO1X1y4HjC/3xDYA/BJ71b98PfFtEvow36H0G8IsYZTF14Cpal0J7rC6D4KKvKOnkcBBxf+jByrjU6bJuK9e9h/qAwtlK4a6uYMUSd4poNqf5fFfBn6tSA6elXIVGfc+2dDIyoKaSCT5w++MFj7UH9lBfsXB6QauqWPdJ8H1a8akHyeR01KBWbn6z8Pd2ge7jF53ZMIPYE1mcv7rf+l+jtYCTwFj9QW8EPgBsF5Gn/GM3AZeLyNl4rdy9wHUAqvqciHwH2IE3w+oGmyHVmFywcDOd+gb7C2Z7uP7j7V1H88/Z0d2L4H2oUhEpR4ayyo7u3oLKq1xRzc6ceov8oiqR1nSSGVPS+TGU0bpP3LTZYL6rFQs7Is8PB9NS1yuMdhU61mslE4JQ+F6Ifzw8huRSf0QF9FVLZ9GW9qYtL587lVse3AUUrvR2+5e4xyop+L0/ed/2kmYXhS8CqmGir+uIEzB2quqo0zpE5MmxXkRVf0Z00HlglOd8DvjcmCU0DaHYXHcXUMJXuOp/DUYMejsnM7nIwe3wiu5SJcRb5Oeunnd0H2Pl2ofye3O4oLdkdntB7iQ32BsU7hrZffBEflMnJ/y6o/W/x6l0ws+JUxmGA8aUlmQ+95areF1Z3YD1/p5+TgxkeHL/kXyr0AUEr/zezxCcTu3GPX75wmFEhgflK9Fdl0wIJ4ey7O/pZ8aU4eDQ6LOLJoo4ASNOyJzYYbVBlbuitlIrcd0Mnt/8j3eycu1DBRWS25NiR/cxegcysQdow9xAdu9AhsdfOMyrOlpob0mho02RiiG4k5/7OdpbUwVBLafku7CCFWrfYJYLvvgo86e3ASO7x1IRYzL9odeN0//uWjZxfk/FusSCpralyOaU/sEsCpwcyhbt03ctsEwml+/OK9biczOnoHA6dA5IhH5N4/3MBceYPn7RmSXNLrp43abY3Vcm2pgBQ1UHAETkUuBBfw3FJ/GmxX5WVZ9w55iJKe4URhGvgkilEpw+b6pfcXlXoJWajeCtDYmaU1WaqP788GLF0ew91Me+Q32cd9rsggp62Zx2MjkdsSfHL184nH8PXOXmumzC05LzZfTHUuKkHXE/T9RU1uCUUhg5kA1e0AxOv02Il/hx98ETBS3A2VNb8tluwwPjAMlkouD8VCpRtZlJccZ1bMOkyipl5PCTqvpdEXkTcCHwReBW4PVVKZmpmHK6L+Js0equKi9etyl/zO2xHRy72Lr38Ihkho0oU2IhFfJXq8tu/DEw3CIBb/A4mRBWLJzuTT1Vb01LsO9/NG4qazHu91hO14wLFq6ydCnwXUvFTaW+4IuP5n+mYoEOhj9Lwat+F3SqWSHHnV3UTIvjGlkpAcNdfr0TLw/Uj0Xks1Uok2kA4T+wcL98cLHWju7egsVfJwYKr4yjUng0m3RSmDO1hZeOnSw4vmXPoaKVf99gdkQXTmsqke/7d2MnUd1OLmAX63kLthridM0Eu56C4w3uNcKB6ZYHd9E7MFQQABdMb+OWB3dFXnC4Y7c8uCs/MJ2KmTCymopNS2709Q6NqpSA0SUi64G3AV8QkVZKS/1j6iRcGcUZwwiu9HZXivt7+gu6HoJyodvZnBasDC6Vm1Ubfm544LZWhrI6IlgAfO3R3fmxkCh9g1l++cLhESupw4IL9YJ5j3LAL144zLS24gPHlVr4FXyvd3QfK2i5wMg9Uopx2YUbSTMtjmtkpQSM9wEXAV9S1SMishD4y+oUy5Sq0imlg39Qn333yvz+AVELtBLiVRJDgRpn9tQWFs2cwuMvHC7r+7e3pvIzi4IarWer2GLBoGCw8GZ9KUNZ5an9R8hkcyPScOwMDLCD9zMf9ycOBFsFwbUpY3XNuNYMFM5aCo6fBMciwoPKwXTmUZ+10Y41imZZHNfI4izcOx/Yoqp9wD+548FkgqaxBSv6FQun5yvxuGMYN/1geP8AJ3j7zAUd/ObgCYLVeZc/ZTTB8PTZUmRzyuypLfQNNufe08UEB4Td+gfXXTdac929h667D+IPisfhkhOqDqf6KDcHVjWFB/BNbcVpYVwJfE1Efg08iDdT6qXqFsuUqtIpD6J2Vwt2SQW/395DfSMWfykUpJAoVd9glqHs5Jp8597B0YKreyzYYoDyWpYu2ITXYaDDj8XNgVXLmUaNvvfGRBZnWu2fAojImcA7gG+KyAzgUbwA8nNbgd3YwruguRk9cccwwFtrcLR/qKA7w7VUqtUSGGqGqVU1NiUwiL541pT8uoKxuM9AsBvLBZx67M9RbiuhGntvNIJmmdobewxDVXcBu4CviMgU4D8AlwJfBjqrUzxTCdv29dB1pH/M/EtuYVN4ULC9JUlHW5oVCwvn9rvAk8kWti7qNTDdzMZ6z1waFTfAHhxIf8+tj+Wn8ELxFdXhVmhUfijwusbcY2vWby4rO2+1BFf459QbiG+0AfaJrKxstaraj5fSo2haD9MYwpv3FGvGb9vXw66XeskpXPH1Ldz9J6vzj/UNZvOroZ3XfOKf87fDs4cqFSyaMfAkpXD/ibimtCQZGMyO+rxgGvHw+1JqxemCgZvSe+aCDn594DjZGrXqyg1A4Z3zls+dantz11DsgCEincDNwNLg81T196pQLjOGuK2GLXsO5a8cc1p8j4ngeUOZ3Kh7UUD8KZZh7S3RWVOjLJ3Tzr5DfQ0dNJZFlNFtNVus8netheDjA2OsMm9JJWjx18G4/E8uwMdNdx5MjQGFCzj3HurLL1xUSt94qlbCO+dVI8GhKa6UFsbdeNNot1P6BZSpoLitBoDVy+cU3J/V3lLwOq6SWb18Dm3pBEOZHOlUouB5ArwmtNL7krMXcdMPtpdc9rjBAqhIptpqakklmD+9jUMnBvMzndwFulI8aBTLnjuak5lcfmKBG4dw6wqmt6XGrDijUmMEJzD0h34v/YPZht1wKJkQFs2c0hCztiabUgLGQVW9v2olMbGN1WoIztMPZ4INVvLvufUxViwcDgQJEZLJBAkRrr5jeOsRhXziOdftsfG68/OvtWB6Kwd7vW6pSvZovHL8JFNKaJHUmpeQ71hk+o6xphKP521y38+tK4gjKjVGcALD9q6jBe+zy2Rb0NWjCv09LB36DWkdgp2H4FgXSAIO7IShftAsHNwF0xbwl4ePcCC5AH66CY6/BMd+C0e7vOf0l7c+Z6O70Q2sDdxfW9bLeSQBkgTUX1rvon45E8LH4bXvh0u+BonGXQ9dSsD4tIh8HXgEyHdau53zTO2sXj6HhHh/+MLIVkOwq2qsrqNgTp2xVuj2D2ZxDwfzR0WtgA4T4KyFHbxy/KSfQHBsB3sHC2YFtbckWTanPd/X3ghdVZVKebJgeiuH+wb9AD9y54oEWf4i9V1uSPnXbCWugNqWO4MruIkhUqTJsO7Ep1mVeH74tRKwreUMtuTOYnVip/eYXymH3eJubBz52LCnh2fC/PSHpRW21jTnfdXb09+Gd/0vGjmBhsRNEy0i3wLOBJ4jMGVcVWPvv11tnZ2dunXr1noXo2JG60N++dhAPs9PIlC/BLONnrmggz9evaxo11Gw7zu8h3TUKuugVFJKStYnwLS2VMEKccEbyC0W05ICC2dOoftIP1mFjrZUfv1BMOVIbQfHh79TCxnuafksqxLP87XMu/hS5lKUZP6cJFkUIZc/JghZkiiKIP7/6r/Wp1J38ZnMlQySIldQaShtDHF3y39nVeJ5tuVCFXtM5T6vaUxbANNPgRmLYMZi7/a0BTBtHrTPhdZpkJoCyTQkWyCRJL9FT35GQYn3i21o3mREZJuqjjnbtZQWxutU9TXjKJOpoGDLwbU0gly3w/tfv4SvPvLrfCsgWLkGs6aWmlG21MyuTkdbmgRe2utUQkb0nQelUwmmt6Xoxgsewf70YELDeAHDq7ATZLk8+SgHdQY7cktQoJu5gQraezXBDQALwe+QIOufmyBLgi25swD4rc4lgZLNB4YcFyae5Ke5sxlCUP97t5Dh6uSDPKfLmMJJ/jW3ihxJhoAeOvhU6i42ZP+AvTqf4XXywtXJB/PB4orBmxgkRQuZfBCJY1Xi+RHnHknMYub8pTBtPptezPCzE6ewK7eYV3cM8Ym3nwZt06G1A9JToaUd0v6XX+leccc2sqS49/o3jvh+xaZpj0dwLUl4fZGpvlJaGHcAX1TVHdUtUvkmWgtjNOHsmwn/SmdgKDti5szKtQ9xYiCDQsGYQGtg5k25SQLjWrHQ2z8bhlNpg5cjKerbnjprSn4WjJvS27lsNqcNPc/nX/nPBeduy53B5YM3M0SKJFkSQIYEORL+v4VX8sFK9muZd/F3mUsLWgGgLOIA/ym5hduy78p/n/NkJ+9O/pzPZK70v1eOf594kv+bO4cMCRIoOT84pMmyNvVN1mauyp/7vuSj/K7sZW3maoZIBloY3u+uI5WlNxO46g1Jysj9JoKPAQV7XISF12G4sRA3luG22HWCY1ZO+LWDKWeCgoPsrrVbicy14b08in3/uCzQeKrRwlgNPCUiL+CNYQhel1SsabUishi4C5iP95e5QVXXichsvN7QZXh7er9PVXtERIB1wMVAH3C1qj5RQnkntHD2TTdI/d3r38Dl/7ClYOC67+RwRRC8oj+ZyTGYyZGQyg5WO8Er/4PHT7JtX0++ReP+wKcc389jB1sY9D+KLQxxT8vnWNX/PNzhd5O3+i9SpN9+S+4sMiRRPzS8N/kTFskhZtFLDx2sTuzMn+dufy3zLlYndrI6sZMWMgyE+o2PMo2uUy4itT9DhhQgPJ38XX6RWZH/2bIkeTh3Xv45wcp2kAQ3ZT6cv58hwbezbyv4HuG3fEDTEUc9rakEZy+eWZAQ0HFddeMV/s5D2RzJRHkz1Wq1/0QjzuKayEoJGB8ADoaO/acSnp8B/kJVnxCRDmCbiDwMXA08oqqfF5EbgRuBv8JLQ3KG//V6bLOmEaKyb65aOmvEftLtranhmTwiBXtZJBJCe0sy1uBt8Frze3/6hoLUIWEXJJ7i7Ylf8snMh8iS4GDvSS6/9adeMHADqr6bkh/knuyFKAmyJNmSO6ukPnZX6Z9EyJHgh9k3k2xr59iAV4W7q+/OZbP5cShv0vf/9A0sv287vz06wJE+l3BROM5UfrS/sBIN5suqRmNsZnuag0UmBJy9eGb+atglBEz5kxTirsQOdt+EE1ICI3JKZbLK4ldNiazow7sEhr9/tfafsC6o+iolYGwArlTVZwFE5DK8IPL3cZ4czG7rb/O6E1gEXAJc4J92J/BTvIBxCXCXen1mW0Rkpogs9F/HjCLqqstVDlNbkgWzmuZNa2HpnKkF54A3ZTScUHCK30VxNrtYdcf72ettaR05mLotdwZrBj/pdxB5NfYQqchg8J7kz/in7FsYJIWgzMKrtH6bXMSWKW9mW+tqPvdnH4Bk4cd15dqH6B3I0NGWItGiMJjlrMA04WCFWCxv0pY9h+hoS/OatjS/8LdRXTC9lZeOnazqQHraj2IuX1ZbOsFtf9zJB25/fMRkg/aWZEEF6RICOqUXehl0AAAX0ElEQVTs/T2ajrZ0wWQGpfyWQbX2nwgmHrR1GLVXSsB4L/A9EXk/8Ga8LLZvL+ebisgyvD3BHwfmB4LAS3hdVuAFk/2Bp73oHysIGCJyLXAtwJIlS8opzqTggkjXkcIkgalkeAqf8jqe4xupz/Dh3Md4ONcJJEiS5YbcvdyQGLkUJ2owdUvuLHKhAeM0GVYndvJCbgEPJN7Cw7lOjrSfxoG+HAP+ho454KbMh/nWvD8vrKiSo39UVy6aAYyd/nrdZefkW0ZtaW+B4qZfew1n14LKxRzXG49wYkVVuPqOXzAQMQlgtA2aShEMKMUGjMOpN8ZT2Y+1/0SpK8lH25c7LmuZjE8pyQf3+K2K+4B/A97u55QqiYhMA74PfExVj0lgWpqqqoiU9NeqqhvwWj90dnY2wtT8uhs7C+3PSZIbno/fH939c33qx/y/wdfm5+67/v9RrXg3/3nHq/ll4tXeYigFEGa2p7j9qjdyy4ML2Lr3MK1pf/B9cCjyZUq9sg1eeYa7S6Bws6D2Fm8LUTdLbEf3MbI5zY9BRHULnTprCqmE5LceDW5dWg6Rwu1XNfT/WNze26kK724XTr3RSFfxti93/cXZQGk7hZ/j2UASeFxESsolJSJpvGBxd2DB38uuq8nfxe+Af7wLWBx4+qn+MRNHdgh2PwL3rCk4vAr4fniBVhGrEs9zd8t/z5+bbpvGpzu+xK/SZ3Hv9W/Kz3ZyrRdXIR9Yv5mlwOyBoXxf+2vmT89XPu2tKWZMSRdNiZ4QxsyRVfCj5jR/5ekSJ0btae26bVxrxFmxcDpdR/rz5XEf9pbU8Iyko/1DBT/rgd6T9A9my+62CjdiTp/nJdELr7h2gsEveM4gXsvoyf1HuHjdprKDbPD9qlXqjVKv9m1f7vqL81f5B5X4Rv6sp9uBnar65cBD9wNXAZ/3//9h4PhHRORevMHuozZ+MdKZg8/CHZ9lOz/zDqwd+zlRXUgA/OEGWHEJpL3BCZe+3K3SDldG7iq3d2Co4DFXEU1vS40YgIfo7T9dzZtKJfKVJ4xeqbjB+6FAenWXODFYkbjd5IK7CIZfO7yHeVs6ycpFM0YMDAdbL7/09wNpSSUKxnvKWUi491AfyYREdj9llYJsweEZbYpLU9LrjRgVaXC46bNQuDAz2L3jpq3G3ckvmMAw7mPldgvZvtz1F2cDpX0V+l5vxBsk3y4iT/nHbsILFN8RkWuAfXh7h4OXOv1iYDfetNoPVqgcE8NLz7Kx+yLv9uiJZYddvhHOeBvLbnoQgL2ff2fBw2vWb4YtTwKFA+Aul1SwMrp43aZ8heMeW7n2ocgV4smICiz4xz+9LcWvXvJaCHF2dwOvgnffZyir+Uo6nDgxeDW+o7s3MlGjq9BWLCwsT9RmQ+7c4H4UUbsNlsrtgb6j+1h+xlpSCit5J7xeIiiVSkQGaSicDBFsUTVT947ty11fcbqknlDVc8d7jqr+jGIrkuDCiPMVuGGs8k1a015VcPffmM+XuJKvfvrmqqUrSKUSZP2r+fBsI1dJlpL2vJTkeWHhhIsz21NMbU2PuPIMl2fLnkMjrnCDLaJFM6fE+v6Vzjw0vS2VH4gOLlR0gmVes34z27uOksl5gdIFrLhpzqF49054JfVYRpvmalNgJ544XVJnicgzozwuwIxRHjfVMO1VrFnotRQ2Xnc+73Srp8cRLMJ/2Mtv/HG+YnSVkVt4F+xSClZU4dTrbgvNSgunbT/Sl2HhjJH97qfPm5rfNyKcqDFqP2u3utmtJC623uD0mx8omh6lNdRFNRa3o6F7/ZVrH4rMgOu4wOKCRmsqwbyO1pK6aax7x5QjTsAYe8NgaMz80yZSVBdRlKltXr6mltTIK9diFU64qynObmhxKsmwVUtnRa4ZCOtoS5NOehW4Ap+4bzvf2rI3n6Yk/BzXInFjM8W8+lXT8oGmJSmcMnMKmZxy+MRgPutv3EHxlYtmFH2PXOsnPKAdHhw/fGIwH8zjXtFXonvH1kVMLrUcwzBV5AagK/2HmxQ4Z/HMyNcsVuG446V0RQQX18UVrLSh+Myqgt3wdDhIhPezTojXCgmPeUT97A989C1s29fDlj2HWL18Tr6y/tlfvRUYTpLn0q9EaU0lRuyHvmb95nxer3BQCOZ2Cne1lbsD4mjGGvQea11EcNzHuqUmhrL29DaVVe5WmO4PMjgA/d5bH4uc0RTkelKivm8z/WEHf8Zgt07Y6fOmjliM5oRbRLsPnih4bnjGVdCqpbMKHtvedZQ3feEnrLvsnPz7GFwI56STwhmvmsb+nn5aUokR77mbcjsYCiZuYBxGJve758Pxxi4qydZFTD4WMBpEqVfXMJxt9PkDx/PH1L/flh47aVw537NRA8po4yRjLUZz3VMwcjA4PFYS5eJ1m/j1geNkskrfYH/BWpDg957elmJ/T/+I1NzBwN07MJQfNwqOkYTXHaxZvzm/APH0eVOrEizG+l2PtS7CBr0nHgsYDWCsdBZj2fPK8YLVyTPb0yyfO63o+eNNCd1sgu/vWBVrsMWRSggfvffJghZJWFQLIrwWxC2E80QvVnS/k6jBcsHbuyRcdhck6zWGYAPnk09ZAUNEPqKqf+/fnqmqRypbrMlnPFdhwUVnLUnhtj/uHPWPt5Qrv/bWVFNcIcadBhpHR1uavYf6Cha2LZndHjlGFDXQHlwLElw0+MBH31K0G9AdD8/acq/X0ZYe8Vy35qVvsH/EGEKtfme2LmJyKbeFsTRw+ybg4xUoiynTqqWz8t0TcfqyGyUANEo5ooR3NHS5o8KpR4LJDFMJWPO6JfzRuadyy4O7Cir/Hd29o6buCL4XF3zx0YLtd0+fNzWyFRrc26RSYwjjbe2aia3cgJEQkTcDP8fLLWXqzE3ltCu9ygiv33DhI9zd9Mn7hvdLTyQS/NG5p+YfC7c+whV6sYA5f3obB3pP5sc9irUaVnzqwXwryHIrmVqIHTBE5OfAJ1T1UeAvgT/F2/zovuoUzUw21RxTKbU109GWzvfP/9kFp/OJ+7aT05GpR4JBIZsdDiZRU3ZLqdALxz2Kn5MATvG3s50MwaKRW6WTQSktjOuAvxGRTwA3q+rXqlQm0yC2r/39UR8v54+3GRZ6Rf1c39qyN3Jwd91l5+SDQjiY1GJQWISKZ5atVKVslfvEU8p+GM8C7xGRc4HP+PtY3KyqT1ercGbiCKfhCA7Sjqdi6SghDfp4FBvcXbV0Ft+9/g35BXzFFjE2anA0phTl/LXtBv4W+BCwrczXMJPQRF3oFV7A1yysBWBKVcoYxqPAGXgTyXf4X1dXp1immUVVROE+fRukNab5lNI6+Au8zY9K3pbVGLCFXsY0u+idViKo6hMWLMx4dbSla7L952SwYuH0yA2WjKmW2AHDGFM5a9ZvtgVypunULGCIyDdE5ICIPBs4tlZEukTkKf/r4sBjfy0iu0XkVyIy+vxOM+H0DgzRdaQ/ckvVoGxOI/frLuU1GlElU50YUym1bM9+E/h74K7Q8a+o6peCB0RkBXAZ8DvAKcC/isirVdU2aprgdnQfK9gX/NLbHuO7178hsgurd2Co6N4V2/b15KfwhtN5NLpgoBgtaATzVNWLBbXJpWYBQ1U3iciymKdfAtyrqieBF0RkN3AeYG34SSCcx6nYnhTh1BvB84J7fofTeZSj3hVjuPtqrDxV9S6vmZgaYQzjIyLyjN9l5f6iFwH7A+e86B8zE9yKhdM5fd5U3PYWo+1JEdxhLyGF+3yvXj4n/xrhFdgTQVSeKmOqrd5TLG7FWwSo/v9/h7cgMDYRuRa4FmDJkiWVLp+pg2Aep9Gm37rNiaI2EZpoU3jDLYbx5Kkyplx1bWGo6suqmlXVHPAPeN1OAF3A4sCpp/rHol5jg6p2qmrnvHnzqltgUzNxp98mE0Krv19Eua/RjFxAbEklOHNBx4T8GU3jqWvAEJGFgbt/CLgZVPcDl4lIq4ichrfC/Be1Lp+pPTezqZ4DubVQiRlcHW3posHSmGqoWZeUiNwDXADMFZEXgU8DF4jI2XhdUnvxMuKiqs+JyHfw0o9kgBtshtTEt73raH7WE0B7y9j7kjei4AC123o1vG93VBLG0RQbxLbpt6aWajlL6vKIw7ePcv7ngM9Vr0Sm0QRnR0Xdb0ZRe3xM1CSMZuKr96C3MXnBXe7c1qTNaKwrfkvCaJqVBQzTMIKzo4Jbk040E20Gl5k8LGCYhuI2HJrobGMl04waYeGeMcaYJmAtjAlirP23J7JsTvNTVO2K3ZjqsRaGaWonBjL0DWZ5saefK76+pSkz0xrTLCxgmKYWnHjrkgwaY6rDuqRMU5vSkswv9mumJIOVWGxnC/ZMrVnAMDVV6UoumRDaW5LMntpiU1SNqTLrkjINYeN155cdTJIJmbBJBo1pJBYwjDHGxGIBwzStjdedH5mryRhTHRYwTFOzoGFM7VjAMMYYE4sFDGOMMbFYwDDGGBOLBQxjjDGx1CxgiMg3ROSAiDwbODZbRB4Wkef9/2f5x0VEvioiu0XkGRE5t1blNMYYE62WLYxvAheFjt0IPKKqZwCP+PcB3gGc4X9dC9xaozIaY4wpomYBQ1U3AYdDhy8B7vRv3wm8O3D8LvVsAWaKyMLalNQYY0yUeo9hzFfVbv/2S8B8//YiYH/gvBf9Y8YYY+qkYZIPqqqKiI59ZiERuRav24olS5ZUvFymPiwTqzGNp94tjJddV5P//wH/eBewOHDeqf6xEVR1g6p2qmrnvHnzqlpYY4yZzOodMO4HrvJvXwX8MHD8Sn+21GrgaKDryhhjTB3UrEtKRO4BLgDmisiLwKeBzwPfEZFrgH3A+/zTHwAuBnYDfcAHa1VOY4wx0WoWMFT18iIPXRhxrgI3VLdExhhjSlHvLiljjDFNwgKGMcaYWCxgGGOMicUChjHGmFgsYBhjjInFAoYxxphYLGAYY4yJxQKGMcaYWCxgGGOMicUChjHGmFgsYBhjjInFAoYxxphYGmYDJWPKZZstGVMb1sIwxhgTiwUMY4wxsVjAMMYYE4sFDGOMMbFYwDDGGBNLQ8ySEpG9QC+QBTKq2ikis4GNwDJgL/A+Ve2pVxmNMWaya6QWxn9Q1bNVtdO/fyPwiKqeATzi3zfGGFMnjRQwwi4B7vRv3wm8u45lMcaYSa9RAoYC/yIi20TkWv/YfFXt9m+/BMyPeqKIXCsiW0Vk68GDB2tRVmOMmZQaYgwDeJOqdonIq4CHRWRX8EFVVRHRqCeq6gZgA0BnZ2fkOcYYY8avIVoYqtrl/38A+AFwHvCyiCwE8P8/UL8SGmOMqXvAEJGpItLhbgNvB54F7geu8k+7CvhhfUpojDEGGqNLaj7wAxEBrzzfVtUHReSXwHdE5BpgH/C+OpbRGGMmvboHDFXdA7w24vgh4MLal8gYY0yUundJGWOMaQ4WMIwxxsRiAcMYY0wsFjCMMcbEYgHDGGNMLBYwjDHGxGIBwxhjTCwWMIwxxsRiAcMYY0wsFjCMMcbEYgHDGGNMLBYwjDHGxGIBwxhjTCwWMIwxxsRiAcMYY0wsFjCMMcbEYgHDGGNMLA0dMETkIhH5lYjsFpEb610eY4yZzBo2YIhIEvga8A5gBXC5iKyob6mMMWbyatiAAZwH7FbVPao6CNwLXFLnMhljzKTVyAFjEbA/cP9F/5gxxpg6SNW7AOMlItcC1/p3T4rIs/UsTwnmAq/UuxAlqGl5v3P9uJ5eUNZxvla12eegeqys8S2Nc1IjB4wuYHHg/qn+sQKqugHYACAiW1W1szbFG59mKis0V3mtrNXTTOW1slZeI3dJ/RI4Q0ROE5EW4DLg/jqXyRhjJq2GbWGoakZEPgI8BCSBb6jqc3UuljHGTFoNGzAAVPUB4IESnrKhWmWpgmYqKzRXea2s1dNM5bWyVpioar3LYIwxpgk08hiGMcaYBtKUAWOslCEi0ioiG/3HHxeRZbUvJYjIYhF5VER2iMhzIvLRiHMuEJGjIvKU//WpepQ1UJ69IrLdL8vWiMdFRL7qv7fPiMi5dSrnawLv2VMickxEPhY6p27vrYh8Q0QOBKd5i8hsEXlYRJ73/59V5LlX+ec8LyJX1bG8XxSRXf7v+QciMrPIc0f9zNSorGtFpCvwu764yHNrmm6oSFk3Bsq5V0SeKvLcmr6vsahqU33hDYD/BlgOtABPAytC5/wZcJt/+zJgY53KuhA417/dAfw6oqwXAD+q9/saKM9eYO4oj18M/DMgwGrg8QYocxJ4CVjaKO8t8BbgXODZwLFbgBv92zcCX4h43mxgj///LP/2rDqV9+1Ayr/9hajyxvnM1Kisa4H/FuNzMmrdUYuyhh7/O+BTjfC+xvlqxhZGnJQhlwB3+re/B1woIlLDMgKgqt2q+oR/uxfYSfOvVr8EuEs9W4CZIrKwzmW6EPiNqu6rcznyVHUTcDh0OPi5vBN4d8RTfx94WFUPq2oP8DBwUdUK6osqr6r+i6pm/Ltb8NZC1V2R9zaOmqcbGq2sfp30PuCeapahkpoxYMRJGZI/x//AHwXm1KR0RfjdYucAj0c8fL6IPC0i/ywiv1PTgo2kwL+IyDZ/FX1YI6ZsuYzif3SN9N7OV9Vu//ZLwPyIcxrx/QX4EF7LMspYn5la+YjfffaNIt19jfbevhl4WVWfL/J4o7yvec0YMJqOiEwDvg98TFWPhR5+Aq8r5bXA/wLuq3X5Qt6kqufiZQm+QUTeUufyjMpf1Pku4LsRDzfae5unXp9DU0xRFJGbgQxwd5FTGuEzcyvw74CzgW68rp5Gdzmjty4a4X0t0IwBI07KkPw5IpICZgCHalK6EBFJ4wWLu1X1n8KPq+oxVT3u334ASIvI3BoXM1ieLv//A8AP8JrxQbFSttTQO4AnVPXl8AON9t4CL7vuO///AxHnNNT7KyJXA38AXOEHuRFifGaqTlVfVtWsquaAfyhShoZ5b/166Y+AjcXOaYT3NawZA0aclCH3A252yXuBnxT7sFeT30d5O7BTVb9c5JwFbnxFRM7D+53UK7hNFZEOdxtv0DOczPF+4Ep/ttRq4Gigm6Ueil6lNdJ76wt+Lq8CfhhxzkPA20Vklt+t8nb/WM2JyEXAx4F3qWpfkXPifGaqLjSO9odFytBI6Yb+I7BLVV+MerBR3tcR6j3qXs4X3kydX+PNeLjZP/YZvA82QBteF8Vu4BfA8jqV80143Q7PAE/5XxcD1wPX++d8BHgOb8bGFuANdXxfl/vleNovk3tvg+UVvI2tfgNsBzrrWN6peAFgRuBYQ7y3eEGsGxjC6yu/Bm8c7RHgeeBfgdn+uZ3A1wPP/ZD/2d0NfLCO5d2N1+fvPrtu5uEpwAOjfWbqUNZ/9D+Pz+AFgYXhsvr3R9QdtS6rf/yb7nMaOLeu72ucL1vpbYwxJpZm7JIyxhhTBxYwjDHGxGIBwxhjTCwWMIwxxsRiAcMYY0wsFjCMMcbEYgHDGGNMLBYwjIlBRK4TkZf8RIa/EZErI85ZJiL9bn8D//64VueKyBR/P4TBOqc1McYChjExrQTWqpfI8HIgMtULXpr1syv1TVW133+931bqNY0plwUMY+L5PeBX/u0XgMFSniwiy0XkSRF5nb+L3d0islNEvici7YHzrvRTdD8tIv9YwfIbM24WMIyJZyXwKz+Z4UeAm+M+UUReg5ex+GrgIPAa4H+r6lnAMbwdIvH36/gE8Fa/JTNiS19j6skChjFjEJHFeFvsPoCXkvxNeMnj4piHl5X2ClV92j+2X1V/7t/+lv96AG8FvquqrwCoajm7yhlTNRYwjBnbSmCTP5bwauBM4PyYzz0K/BvDQQFGbpxkGUBNU7CAYczYfg94EkC9fba/Dbwz5nMH8fZnuFJE3u8fWyIiLuC8H/iZf/snwKUiMgdARGZXoOzGVIwFDGPGthI/YPj+D96+CrGo6gm8Xev+K952sr/C23JzJzALb3tRVPU54HPA/xWRpyk+E8uYurD9MIypEBFZBvxIVX93POcUed5evM2qXim/hMaMj7UwjKmcLDDDLdyrBLdwD0gDuUq9rjHlsBaGMcaYWKyFYYwxJhYLGMYYY2KxgGGMMSYWCxjGGGNisYBhjDEmFgsYxhhjYrGAYYwxJhYLGMYYY2L5/zRa8jQljpKbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "ax1.errorbar(vecRp, vecvRp, xerr=vecerrRp, yerr=vecerrvRp, fmt='.')\n",
    "ax1.axis([0., max(vecRp)+0.1*(max(vecRp)-0.), min(0.,min(vecvRp)-0.1*(max(vecvRp)-min(vecvRp))) , max(vecvRp)+0.1*(max(vecvRp)-min(vecvRp))])\n",
    "ax1.set_xlabel('$R$ [kpc]')\n",
    "ax1.set_ylabel('$v_c$ [km/s]')\n",
    "plt.plot(vecRp,Vc_Model_I(vecRp,params))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Xi2(parameters):\n",
    "    rho_b,rho_h,alpha_h,a_h,beta_h,R_d,Sigma_d_SD,Sigma_g_ISM = parameters\n",
    "    model = Vc_Model_I(vecRp,parameters)\n",
    "    xi = np.sum((vecvRp-model)**2./(vecerrvRp)**2.)\n",
    "    return xi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7142.2492346595245\n"
     ]
    }
   ],
   "source": [
    "print Xi2(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.27000e+08  7.11000e+08 -2.00000e+00  3.83000e+00  2.96000e+00\n",
      "  2.00000e+00  1.42875e+09  4.76250e+08]\n"
     ]
    }
   ],
   "source": [
    "print params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.array([np.log(1.0e8),np.log(1.0e8),-3.,1.,1.,1.,np.log(0.15e9),np.log(1.0e8)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
      "  If increasing the limit yields no improvement it is advised to analyze \n",
      "  the integrand in order to determine the difficulties.  If the position of a \n",
      "  local difficulty can be determined (singularity, discontinuity) one will \n",
      "  probably gain from splitting up the interval and calling the integrator \n",
      "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
      "  warnings.warn(msg, IntegrationWarning)\n"
     ]
    }
   ],
   "source": [
    "LS = minimize(Xi2,x0,method='L-BFGS-B',bounds=((np.log(1.0e8),np.log(3.0e9)),(np.log(1.0e8),np.log(2.0e9)),(-3.,3.),(1.,9.),(1.,3.),(1.,4.),(np.log(0.15e9),np.log(3.0e9)),(np.log(1.0e8),np.log(20.0e8)),))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: 751388.4414566102\n",
      " hess_inv: <8x8 LbfgsInvHessProduct with dtype=float64>\n",
      "      jac: array([-1.16415322e-02, -3.19781248e+02, -1.37382734e+05, -2.23935349e+03,\n",
      "        2.36101914e+02, -4.19095159e-01, -4.65661287e-02, -2.32830644e-02])\n",
      "  message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
      "     nfev: 63\n",
      "      nit: 3\n",
      "   status: 0\n",
      "  success: True\n",
      "        x: array([19.06096501, 21.41641302,  3.        ,  9.        ,  1.        ,\n",
      "        4.        , 21.82187813, 20.51096785])\n"
     ]
    }
   ],
   "source": [
    "print LS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
      "  If increasing the limit yields no improvement it is advised to analyze \n",
      "  the integrand in order to determine the difficulties.  If the position of a \n",
      "  local difficulty can be determined (singularity, discontinuity) one will \n",
      "  probably gain from splitting up the interval and calling the integrator \n",
      "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
      "  warnings.warn(msg, IntegrationWarning)\n"
     ]
    }
   ],
   "source": [
    "L2 = minimize(Xi2,x0, method='TNC',bounds=((-3.,3.),(1.,4.),))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: 2016.0642953358938\n",
      "     jac: array([0.0001819 , 0.00102318])\n",
      " message: 'Converged (|f_n-f_(n-1)| ~= 0)'\n",
      "    nfev: 47\n",
      "     nit: 11\n",
      "  status: 1\n",
      " success: True\n",
      "       x: array([-0.32976035,  3.7299834 ])\n"
     ]
    }
   ],
   "source": [
    "print L2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(vecRp,Vc_Model_I(vecRp,np.asscalar(LS.x)),'.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Grid\n",
    "R_d_test = np.linspace(1.,4.,10)\n",
    "Sigma_d_test = np.linspace(0.1e9,3.0e9,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "#g=np.meshgrid(R_d_test,Sigma_d_test)\n",
    "#np.shape(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix = []\n",
    "#for i in range(len(R_d_test)):\n",
    "#    row=[]\n",
    "#    for j in range(len(Sigma_d_test)):\n",
    "#        row.append(Sigma_d_test[j])\n",
    "#        matrix.append(row)\n",
    "#print np.shape(matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Vc_Model_I_grid(R,parameters,rho_b=rho_b_I,alpha_h=alpha_h_I,beta_h=beta_h_I,rho_h=rho_h_I,a_h=a_h_I,Sigma_g_ISM=Sigma_g_ISM_I):\n",
    "    R_d,Sigma_d_SD = parameters\n",
    "    def Vc_b(R):\n",
    "        def I_b(x,R):\n",
    "            e = np.sqrt(1. - q_b**2.)\n",
    "            n = x**(2.-alpha_b)*np.exp(-(R**2.*x**2.)/r_b**2.)\n",
    "            d = np.sqrt(1.-x**2.*e**2.)\n",
    "            t = R**2.*(R/a_b)**(-alpha_b)\n",
    "            return (t*n)/d\n",
    "        Integral_b = []\n",
    "        for i in range (0,len(R)):\n",
    "            result = integrate.quad(I_b,0,1,args=(R[i]))[0]\n",
    "            Integral_b.append(result)\n",
    "        return np.array(Integral_b)\n",
    "    def Vc_DM(R):\n",
    "        def I_h(x,R):\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",
    "        Integral_DM = []\n",
    "        for i in range (0,len(R)):\n",
    "            result = integrate.quad(I_h,0,1,args=(R[i]))[0]\n",
    "            Integral_DM.append(result)\n",
    "        return np.array(Integral_DM)\n",
    "    def V2_c_SD(R):\n",
    "        result = (R_d/R) - np.exp(-R/R_d) - (np.exp(-R/R_d)*R_d/R)\n",
    "        return result\n",
    "    def Vc_ISM(R):\n",
    "        def M_g_ISM(u,R):\n",
    "            return R*u*np.exp(-((R*u)/(2.*R_d))-(R_m/(R*u)))\n",
    "        Integral_M_g = []\n",
    "        for i in range (0,len(R)):\n",
    "            result = integrate.quad(M_g_ISM,0,1,args=(R[i]))[0]\n",
    "            Integral_M_g.append(result)\n",
    "        return np.array(Integral_M_g)\n",
    "    suma = 4.*np.pi*G*q_b*rho_b*Vc_b(R) + 4.*np.pi*G*q_h*rho_h*Vc_DM(R) + 2.*np.pi*G*Sigma_d_SD*R_d*V2_c_SD(R) + 2.*np.pi*Sigma_g_ISM*G*Vc_ISM(R)\n",
    "    return np.sqrt(suma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Xi2_grid(parameters):\n",
    "    R_d,Sigma_d_SD = parameters\n",
    "    model = Vc_Model_I_grid(vecRp,parameters)\n",
    "    xi = np.sum((vecvRp-model)**2./(vecerrvRp)**2.)\n",
    "    return xi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[101060.48867827347, 90840.46525295783, 81487.38721596487, 72917.945183866, 65062.85317741078, 57863.439172412625, 51269.307478658826, 45236.66813004723, 39727.10326638539, 34706.63183290136, 89065.88069781233, 73645.81667669886, 60367.812755251536, 48942.98765786897, 39146.3270506917, 30797.383080628577, 23748.22674455907, 17875.459422641405, 13074.672816742655, 9256.470774554136, 77955.76709218594, 58028.67632112704, 42094.45473625616, 29481.97329386064, 19696.87439465023, 12360.282674742026, 7173.188406215191, 3894.160133609738, 2324.629404586478, 2298.7669272521744, 68363.3207971151, 44970.324921145984, 27847.37033525709, 15765.758237775524, 7862.66000113316, 3501.282806970273, 2194.330622628417, 3558.370879618789, 7284.9038568062515, 13121.140296521722, 60312.40106506189, 34530.34510184152, 17554.797066984887, 7455.367353449926, 2930.709422011213, 3049.9712168468795, 7118.380227822683, 14600.308459670992, 25072.10567912343, 38191.631366889145, 53621.108518231864, 26410.35265079562, 10662.588066943428, 3642.951866106989, 3570.9341827307157, 9205.586303381873, 19639.248767039506, 34183.47449359466, 52300.99659361248, 73562.7719469557, 48065.96299178258, 20220.126704014205, 6513.0783943432, 3341.788293717722, 8429.484146456036, 20223.354863523717, 37606.30040244991, 59742.063995683704, 85984.95673312806, 115823.90279781561, 43439.80495283859, 15586.04607362592, 4499.554065223086, 5677.288993778083, 16345.823845432798, 34647.76207963754, 59265.064691776715, 89221.5681219166, 123770.0321672906, 162323.1865973865, 39567.67350898526, 12187.142997521694, 4114.086860912868, 9937.470801976286, 26397.759061831606, 51345.4603281264, 83272.63018867385, 121072.31092637958, 163903.1144244737, 211107.5920842766, 36307.33888097055, 9760.835128823986, 4950.6902312882585, 15566.77574228605, 37880.6569451511, 69465.89940840362, 108636.66890852794, 154164.5719339123, 205121.39477946048, 260785.3057443381]\n"
     ]
    }
   ],
   "source": [
    "Xi_grid = []\n",
    "for i in range(len(R_d_test)):\n",
    "    for j in range(len(Sigma_d_test)):\n",
    "        param=np.array([R_d_test[i],Sigma_d_test[j]])\n",
    "        Xi_grid.append(Xi2_grid(param))\n",
    "print Xi_grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2194.330622628417, 260785.3057443381)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(Xi_grid),max(Xi_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "#np.std(Xi_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "36"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xi_grid.index(min(Xi_grid))"
   ]
  },
  {
   "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
}