ChiSquare_Model_I.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",
    "%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": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHqhJREFUeJzt3XmYXHWd7/H3t/c16fSSrTvpDiELIUJImgBhFVHDInHuvSBclE2HGQdnxOsdx2109BlHnxlHZcYRjYLAmAsi4ogMKKAMOEgCSSQhISF0lk56X9NLdVdVV9Xv/lEVaDqddKfTnVN16vN6nnrq1K9OdT5Z+pPTv/qdU+acQ0RE/CvD6wAiIjK1VPQiIj6nohcR8TkVvYiIz6noRUR8TkUvIuJzKnoREZ9T0YuI+JyKXkTE57K8DgBQXl7uampqvI4hIpJStmzZ0uGcqxhrv6Qo+pqaGjZv3ux1DBGRlGJm9ePZT1M3IiI+p6IXEfE5Fb2IiM+p6EVEfE5FLyLicyp6ERGfU9GLiPhcUqyjFxFJZs45hqKOUCRKOBIjHI3F74dtD0UdQ9H446Fhj4eiMSIxRyQaH4vEjjx2RGKO2uoZXLJ4zHOeToqKXkRSmnOOwaEo/aEIA6EogXCEgXCUQCjCYDjKQDjKwFCUwXCEwXCMwSPbQ1EGh2IMhqOEIlGCQ1GCQ7H4fSRKaChGKBIjFIkSisSYqo/X/vhlC1X0IuJf0Zijd3CInsEheoND9A5GEvdD9AUj9AWH6AtF6AtG6A9G6A9F6AtF6A8O0R+KEEgU+4mUcG5WBnnZmRTkZJKXnUluVgb5OZnkZWVSXpT11lhe9tvP52ZlkJOVQW5WJjmJ7ZzM+H125tvPZ2XYW2Pxm5GdmUFWppGVEX+clRnfLyvDyMwwzGzq/oATVPQictKcc/SFInQHwnQGwnQHwnQPDNEdCHN4ML7dMzDE4cEwhweGODyQKPNQ5Lhf1wyKcrMozs2iOC+bwtxMpudnU1WST2FuJoW5WRTlZlF45JYTL/DC3CwKcjLJz86iMDeT/JxMCnKyyM/OJDNj6os12ajoRWRUQ9EYnf1h2vtCtPcH6egL094forM/TGcgft/RH6IrEKZ7IMxQdPTD6swMoyQ/m5KCbKbnZzN7Wh5LZhczPT/7rdu0vMR9fjbFeVlMy89mWl4WhTlZZKRhMU82Fb1ImolEY7T3h2jpCdLaG6K1N5i4hWjrC9LeF6KtL17goynKzaKsKIfSwhyqZuRzdlUJpUU5lBXmMKMgPl5SkE1pYQ4zCnMozs06JdMTcmwqehEfcc7R0R+m8fAgTW/dgjT3DNLcE79v7wsRG3HwnZVhVBTnMnNaHvNKC1hVPYOK4lwqinMpL4rfVxTFt/NzMr35zcmEqehFUohzjq5AmINdAxzqHuRQ1wAN3YM0dA/Q2D1Iw+FBwpHYO16Tn53JnJI85k7P55JFFcyenhe/TctjVuJWVpijKRIfU9GLJJlYzNHSG+RAR4ADnQPUdwWo7xigvmuAg50BAuHoO/YvLcxh3ox8zpgzjSuWzaKyJJ+5JfnMLcmjsiSf6fnZmjpJcyp6EY/0BofY1x5gb1s/+zr62d8RYF97gAOdAYJDbx+V52RmMK80n+qyQs5bUMr80gLmlxYwr7SAqhn5FObq21iOT/9CRKZYdyDMntY+9rT1U9faR117P2+29tPWF3prn8wMY35pAQvKC7no9HJqygtZUF5ITXkhs6flpeWSQJk8KnqRSRIcirKntY/dLX28ceTW2kf7sEIvzMnk9FnFXLyogtNnFrGwopDTKoqoLisgO1OXnpKpoaIXmYDO/hA7m3rZ2dTL68297GruZV97/1urWfKyM1g8q5hLF1ewZFYxi2YVsWhWMXOn52m+XE45Fb3IGDr6Q7zW0MP2hh5ea+xhZ1MPzT3Bt56vLIm/EXrV8tksnTONpbOLqS4r1HSLJA0VvcgwgVCE1xp72HboMNsaDrPtUA+NhweB+On4p5UXsnpBKcvnTufMymksmzONkoIcj1OLHJ+KXtKWc46DXQNsqe9mS303fzx4mN0tvW9Nv8wvLeCc+SXcuqaGs6qmc2bldIq0wkVSkP7VStqIRGPsbOrllQNdvHKgiy313XT0x0/zL87NYsX8Ej5x+SLOmVfC2fNKKC3Ukbr4g4pefCscibG94TAb93WyaX+82AcSJxvNLy3gkkUVrKqZQW11KYtmFunMUPEtFb34RjTm2NHYw4t7O3hpbyebD3QzOBQv9iWzivlfq6pYvaCUc2tKmTUtz+O0IqeOil5S2sHOAX5f187v93Twh70d9Abj1zdfPKuI62uruGBhGasXlGkaRtKail5SykA4wkt7O3l+TzvP72mnvnMAgLnT81i7fDYXnl7OmoXlVBTnepxUJHmo6CXpHegI8LvdbTz3Rhub9nURjsbIz87kgoVl3LamhosXV3BaeaFORBI5BhW9JJ1ozLH1YDfP7mrl2ddb2dseAGBhRSE3X1DNu5fOpLZmBrlZui66yHiMWfRmNg94EJgFOGC9c+5uMysFfgrUAAeA651z3RY/rLobuAoYAG51zm2dmvjiF8GhKC/WdfD0zlae3dVKZyBMdqZx3oIyPnx+Ne9ZOov5ZQVexxRJSeM5oo8An3bObTWzYmCLmT0D3Ar81jn3DTP7LPBZ4G+AK4FFidt5wD2Je5F3GAhHeG53O0/taOa53W0EwlGKc7N499KZvHfZLC5dUsG0vGyvY4qkvDGL3jnXDDQntvvMbBdQCawDLkvs9gDwX8SLfh3woHPOARvNrMTM5iS+jqS5gXCE3+1u48nXmvnd7jaCQzHKCnO4dsVc1i6fwwWnlZGTpas4ikymE5qjN7Ma4BxgEzBrWHm3EJ/agfh/AoeGvawhMfaOojezO4A7AObPn3+CsSWVhCJRXtjTwePbmnj29VYGh6KUF+Vy3ap5XPWuOaxeUKoLgIlMoXEXvZkVAT8H7nLO9Q5f4eCcc2bmjvniUTjn1gPrAWpra0/otZL8YjHHywe6+OWrjTz5Wgs9g0PMKMjmT1ZW8oGz5qrcRU6hcRW9mWUTL/kNzrnHEsOtR6ZkzGwO0JYYbwTmDXt5VWJM0sDe9n5+sbWRX/yxkcbDgxTkZPL+M2dz7Yq5XHR6uT5cQ8QD41l1Y8C9wC7n3LeGPfU4cAvwjcT9L4eNf8LMHib+JmyP5uf9rTc4xK+2NfGzzQ28eugwGQYXL6rgr9+/hPedOYuCHK3iFfHSeL4DLwQ+ArxmZq8mxj5PvOAfMbOPAvXA9YnnniS+tLKO+PLK2yY1sSQF5xwb93Xx01cO8tSOFkKRGItnFfH5q5bywRWVzNS1ZESSxnhW3fw3cKzJ1PeMsr8D7jzJXJKk2vqCPLqlgUdeOcSBzgGK87K4rraK61bN46yq6To7VSQJ6WdqGVMs5vjD3k7+38v1PL2zlUjMcd6CUj55xSKuXD6HvGydoSqSzFT0ckyHB8I8uqWBDZsOsr8jQElBNreuqeHG8+azsKLI63giMk4qejnKzqYeHvjDAX75ahOhSIxV1TP4q/ecrqN3kRSlohcg/jF7v97Zwv0vHmBzfTf52Zn8j5VV3HxBNWfMmeZ1PBE5CSr6NHd4IMxDLx/i3186QFNPkPmlBXzx6jO4btU8phfoOjMifqCiT1P7OwL8+MX9/GxzA4NDUdYsLOMr65Zz+dKZOmNVxGdU9GlmS30XP3h+H8/saiU7I4N1K+Zy+0ULND0j4mMq+jQQizme3dXKD17Yx5b6bqbnZ3PnZadz85pqZhbrxCYRv1PR+9hQNMbjrzbx/ef38mZbP5Ul+fzdB5Zx/bnzdFkCkTSi73YfCg5F+dnmQ3z/+X00Hh5k6exi7r5hBVe/aw5ZuqiYSNpR0fvIYDjKhk31rH9hH219IVbOL+Gr687k8qUzdWkCkTSmoveBgXCEDRsP8oMX9tLRH+aC08r4zg0ruOC0MhW8iKjoU1lwKMqGTQe557/q6OgPc9Hp5XzyikWcW1PqdTQRSSIq+hQ0FI3xyOZD/Otv62jpDbJmYRn3fHixCl5ERqWiTyHOOZ58rYVvPv0G+zsCrKqewbc+dDZrFpZ7HU1EkpiKPkVs2tfJPzy5i20NPSyZVcy9t9TqTVYRGRcVfZKr7wzw9Sd38+udLcydnsc/X3c2HzynUpcpEJFxU9EnqUAowr/+ro57/3sf2ZkZ/N/3LeZjF5+mywSLyAlT0SeZI/Pwf/+fr9PcE+R/rqzib9Yu0WewisiEqeiTyKGuAb7wHzt4YU87y+ZM47v/+xxWVWsljYicHBV9EojGHPf/4QDf/M0bZBh8+QPLuPmCGs3Di8ikUNF7rL4zwP95ZBtb6rt595IK/v5P3kVlSb7XsUTER1T0HnHO8dNXDvHVJ14nM8P49ofO5oMrKrVcUkQmnYreAz2DQ3zm0W38ZmcraxaW8c3rzmaujuJFZIqo6E+xHY09/MWGrTQdHuSLV5/B7RcuIENz8SIyhVT0p9Ajrxzii7/cQVlhDj/9swtYVT3D60gikgZU9KdALOb4+lO7+OHv93PR6eXcfcMKyopyvY4lImlCRT/FBsNR7vrpH/nNzlZuvqCaL12zTJ/yJCKnlIp+CnUHwtz645fZ3tjDl65Zxm0X1mhVjYiccmMeWprZfWbWZmY7ho2tMLONZvaqmW02s9WJcTOzfzGzOjPbbmYrpzJ8MmvvC3HjDzeyq6WP9R+p5faLFqjkRcQT45lDuB9YO2LsH4GvOOdWAF9KPAa4EliUuN0B3DM5MVNLS0+QG9a/RH3nAD++9Vzeu2yW15FEJI2NWfTOuReArpHDwLTE9nSgKbG9DnjQxW0ESsxszmSFTQVtvUE+tP4lWntDPHD7ai48XR8KIiLemugc/V3Ab8zsm8T/s1iTGK8EDg3bryEx1jzhhCmkNzjEzfe9THtfiJ987DxWztfySRHx3kSXf3wc+JRzbh7wKeDeE/0CZnZHYn5/c3t7+wRjJI/gUJQ/fWAzdW39fP/Dq1TyIpI0Jlr0twCPJbZ/BqxObDcC84btV5UYO4pzbr1zrtY5V1tRUTHBGMkhGnPc9fCrbNrfxT9ffzaXLE7t34+I+MtEi74JuDSxfTnwZmL7ceDmxOqb84Ee55zvp22+/cwefr2zhb+9ZhnrVlR6HUdE5B3GnKM3s4eAy4ByM2sAvgz8KXC3mWUBQeIrbACeBK4C6oAB4LYpyJxUnt7Zwnefq+P62ipuv7DG6zgiIkcZs+idczce46lVo+zrgDtPNlSq2Nfez6cf2cZZVdP56rrlWicvIklJ5+JP0EA4wp//ZAtZmcb3blqpD+0WkaSlSyBM0Nef3M2bbf08ePtqqmYUeB1HROSYdEQ/AS/WdfDvG+u5/cIFXLxIK2xEJLmp6E9QX3CIzzy6ndPKC/nr9y/xOo6IyJg0dXOCvvafu2juGeTRj6/RvLyIpAQd0Z+A37/ZzsOvHOKOSxbqzFcRSRkq+nEaisb4yq9ep6asgLuuWOR1HBGRcVPRj9OGjfXUtfXzhauXacpGRFKKin4cugNhvv3sm1x0ejlXnDHT6zgiIidERT8O33l2D33BIf72mmU6+1VEUo6Kfgxvtvbxk00Huem8apbMLvY6jojICVPRj+Fbz+yhICeTT713sddRREQmREV/HHta+3hqRwu3ramhtDDH6zgiIhOioj+Of3uujoKcTG67cIHXUUREJkxFfwz7OwL8alsTHzm/mhk6mheRFKaiP4bvPVdHdmYGH7v4NK+jiIicFBX9KA51DfCLPzZy4+r5VBTneh1HROSkqOhH8aPf78MM/uxSHc2LSOpT0Y8QCEX4+dZGrjlrLnOm53sdR0TkpKnoR/jVtib6QxFuOm++11FERCaFin4Y5xw/2VTP0tnFrKrWZYhFxB9U9MNsb+hhR2MvN503X9e0ERHfUNEPs2FTPQU5mXzwnEqvo4iITBoVfULP4BCPb2ti3YpKivOyvY4jIjJpVPQJj21tIDgU05uwIuI7KvqER7c0cHbVdJZXTvc6iojIpFLRE7+uzc6mXq5dobl5EfEfFT3wxLYmzODqd83xOoqIyKRT0QNPbG/m3OpSZk/P8zqKiMikG7Pozew+M2szsx0jxv/SzHab2U4z+8dh458zszoze8PM3j8VoSfTm619vNHaxzVn62heRPwpaxz73A98F3jwyICZvRtYB5ztnAuZ2czE+DLgBuBMYC7wrJktds5FJzv4ZPnV9mYyDNYun+11FBGRKTHmEb1z7gWga8Twx4FvOOdCiX3aEuPrgIedcyHn3H6gDlg9iXknlXOOJ7Y3cd6CMmYWa9pGRPxponP0i4GLzWyTmT1vZucmxiuBQ8P2a0iMJaVdzX3saw9o2kZEfG08UzfHel0pcD5wLvCImZ3QxdvN7A7gDoD58705SemJ7U1kZhhXLlfRi4h/TfSIvgF4zMW9DMSAcqARmDdsv6rE2FGcc+udc7XOudqKiooJxjg5v97RwpqFZZTqM2FFxMcmWvT/AbwbwMwWAzlAB/A4cIOZ5ZrZAmAR8PJkBJ1sBzsH2NcR4PKlM72OIiIypcacujGzh4DLgHIzawC+DNwH3JdYchkGbnHOOWCnmT0CvA5EgDuTdcXN82+2A3DpYm9+mhAROVXGLHrn3I3HeOrDx9j/a8DXTibUqfD8G+3MK81nQXmh11FERKZUWp4ZG47EeGlvB5curtAHjIiI76Vl0W+p7yYQjnLpYs3Pi4j/pWXRP7+nnawM44KFZV5HERGZcmlb9LU1MyjKnehpBCIiqSPtir61N8iu5l5N24hI2ki7on9hj5ZVikh6Sbuif35POxXFuZwxp9jrKCIip0RaFb1zjhfrOrhkkZZVikj6SKui39seoHtgiPMWlHodRUTklEmrot9SH7+s/qqaGR4nERE5ddKs6LuZUZDNabrsgYikkbQq+s313ayqnqH5eRFJK2lT9F2BMPvaA6ys1rSNiKSXtCn6rfXdANRW641YEUkvaVP0m+u7yc40zqqa7nUUEZFTKm2Kfmt9N2fOnU5edqbXUURETqm0KPpwJMa2hsPUan5eRNJQWhT9jqYeQpEYq1T0IpKG0qLoj7wRqxOlRCQdpUXRbz7QzfzSAmYW53kdRUTklPN90Tvn2HKwW9M2IpK2fF/0Lb1B2vtCrJhX4nUUERFP+L7odzf3AXDGnGkeJxER8Ybvi/715l4AluqDRkQkTfm+6Hc191I1I59pedleRxER8URaFL2mbUQknfm66INDUfZ3BFT0IpLWfF30b7T0EXOwTPPzIpLGfF30uxJvxOqIXkTSme+LvjAnk3kzCryOIiLimTGL3szuM7M2M9sxynOfNjNnZuWJx2Zm/2JmdWa23cxWTkXo8drV0seS2cVkZOijA0UkfY3niP5+YO3IQTObB7wPODhs+EpgUeJ2B3DPyUecGOecVtyIiDCOonfOvQB0jfLUt4HPAG7Y2DrgQRe3ESgxszmTkvQENR4epC8YUdGLSNqb0By9ma0DGp1z20Y8VQkcGva4ITE22te4w8w2m9nm9vb2icQ4rl269IGICDCBojezAuDzwJdO5hd2zq13ztU652orKipO5kuNaldzL2awdLaWVopIesuawGsWAguAbWYGUAVsNbPVQCMwb9i+VYmxU25Xcy/VpQUU5k7ktygi4h8nfETvnHvNOTfTOVfjnKshPj2z0jnXAjwO3JxYfXM+0OOca57cyOOjN2JFROLGs7zyIeAlYImZNZjZR4+z+5PAPqAO+CHwF5OS8gQFQhHquwZU9CIijGPqxjl34xjP1wzbdsCdJx/r5OzvCOAcLJpZ5HUUERHP+fLM2AOdAQBqygs9TiIi4j1/Fn1HvOiry3TpAxERfxZ95wCzpuVSkKMVNyIiviz6+s4A1WWathERAZ8W/f6OAWo0bSMiAviw6PtDETr6Q3ojVkQkwXdFX39kxY2mbkREAF8W/QCgFTciIkf4ruj3v7W0Ukf0IiLgw6Kv7wxQUZxLkS5mJiIC+LDoD3RqxY2IyHD+K/oOraEXERnOV0U/EI7Q1hfSEb2IyDC+KvojK260hl5E5G2+KvojFzPTGnoRkbf5q+i1hl5E5Ci+Kvr6zgDlRTkU52V7HUVEJGn4qugP6KqVIiJH8VfRdwxo2kZEZATfFP1gOEpLb1BvxIqIjOCboj/YpaWVIiKj8U3RNx0eBKCyJM/jJCIiycU3Rd/aGwRg1jQVvYjIcL4p+ra+EAAVxbkeJxERSS6+KfrW3iAzCrLJzcr0OoqISFLxUdGHNG0jIjIK3xR9e19Q0zYiIqPwTdHriF5EZHS+KPpozNHeH2LWNB3Ri4iMNGbRm9l9ZtZmZjuGjf2Tme02s+1m9gszKxn23OfMrM7M3jCz909V8OG6AmGiMcfMYh3Ri4iMNJ4j+vuBtSPGngGWO+fOAvYAnwMws2XADcCZidd8z8ymfBnM22vodUQvIjLSmEXvnHsB6Box9rRzLpJ4uBGoSmyvAx52zoWcc/uBOmD1JOYdVVtfvOhnao5eROQokzFHfzvwVGK7Ejg07LmGxNiUauuNnyylN2NFRI52UkVvZl8AIsCGCbz2DjPbbGab29vbTyYGrYmiryjS1I2IyEgTLnozuxW4BrjJOecSw43AvGG7VSXGjuKcW++cq3XO1VZUVEw0BgCtfUFKC3PIyfLFIiIRkUk1oWY0s7XAZ4BrnXMDw556HLjBzHLNbAGwCHj55GMeX1tviJk6WUpEZFRZY+1gZg8BlwHlZtYAfJn4Kptc4BkzA9jonPtz59xOM3sEeJ34lM6dzrnoVIU/oq0vqDdiRUSOYcyid87dOMrwvcfZ/2vA104m1Ilq7Q2yZFbxqfwlRURSRspPakdjjo7+sFbciIgcQ8oXfWcgFD8rVidLiYiMKuWL/sgael3+QERkdKlf9H26/IGIyPGkfNEfOVlKq25EREaX8kXfprNiRUSOK+WLvrUvSJnOihUROaaUb8e2Xn2EoIjI8aR+0ffpIwRFRI4n5Yu+tTeoFTciIseR0kUfjTna+0JaQy8ichwpXfSdgRAxpzX0IiLHk9JF36Y19CIiY0rpon/7Q8FV9CIix5LSRV9SkM3aM2czt0RFLyJyLGNejz6ZraouZdVHSr2OISKS1FL6iF5ERMamohcR8TkVvYiIz6noRUR8TkUvIuJzKnoREZ9T0YuI+JyKXkTE58w553UGzKwdqD+Bl5QDHVMUZ7Io4+RQxsmhjJMj2TJWO+cqxtopKYr+RJnZZudcrdc5jkcZJ4cyTg5lnBypkHE0mroREfE5Fb2IiM+latGv9zrAOCjj5FDGyaGMkyMVMh4lJefoRURk/FL1iF5ERMYp5YrezNaa2RtmVmdmn/U6z0hmNs/MnjOz181sp5l90utMx2JmmWb2RzN7wussozGzEjN71Mx2m9kuM7vA60wjmdmnEn/PO8zsITPz/FNwzOw+M2szsx3DxkrN7BkzezNxPyMJM/5T4u96u5n9wsxKki3jsOc+bWbOzMq9yHaiUqrozSwT+DfgSmAZcKOZLfM21VEiwKedc8uA84E7kzDjEZ8Ednkd4jjuBn7tnFsKnE2SZTWzSuCvgFrn3HIgE7jB21QA3A+sHTH2WeC3zrlFwG8Tj710P0dnfAZY7pw7C9gDfO5Uhxrhfo7OiJnNA94HHDzVgSYqpYoeWA3UOef2OefCwMPAOo8zvYNzrtk5tzWx3Ue8nCq9TXU0M6sCrgZ+5HWW0ZjZdOAS4F4A51zYOXfY21SjygLyzSwLKACaPM6Dc+4FoGvE8DrggcT2A8AHT2moEUbL6Jx72jkXSTzcCFSd8mDvzDPanyPAt4HPACnzBmeqFX0lcGjY4waSsESPMLMa4Bxgk7dJRvUd4v9YY14HOYYFQDvw48T00o/MrNDrUMM55xqBbxI/smsGepxzT3ub6phmOeeaE9stwCwvw4zD7cBTXocYyczWAY3OuW1eZzkRqVb0KcPMioCfA3c553q9zjOcmV0DtDnntnid5TiygJXAPc65c4AA3k83vENinnsd8f+U5gKFZvZhb1ONzcWX2iXt0aiZfYH4FOgGr7MMZ2YFwOeBL3md5USlWtE3AvOGPa5KjCUVM8smXvIbnHOPeZ1nFBcC15rZAeLTX5eb2U+8jXSUBqDBOXfkp6FHiRd/MrkC2O+ca3fODQGPAWs8znQsrWY2ByBx3+ZxnlGZ2a3ANcBNLvnWfi8k/p/6tsT3ThWw1cxme5pqHFKt6F8BFpnZAjPLIf7G1+MeZ3oHMzPi88q7nHPf8jrPaJxzn3POVTnnaoj/Gf7OOZdUR6LOuRbgkJktSQy9B3jdw0ijOQicb2YFib/395BkbxgP8zhwS2L7FuCXHmYZlZmtJT6deK1zbsDrPCM5515zzs10ztUkvncagJWJf6tJLaWKPvFGzSeA3xD/hnrEObfT21RHuRD4CPGj5FcTt6u8DpWi/hLYYGbbgRXAP3ic5x0SP208CmwFXiP+/eT5mZNm9hDwErDEzBrM7KPAN4D3mtmbxH8S+UYSZvwuUAw8k/i++X4SZkxJOjNWRMTnUuqIXkRETpyKXkTE51T0IiI+p6IXEfE5Fb2IiM+p6EVEfE5FLyLicyp6ERGf+/896iyGM/zGIwAAAABJRU5ErkJggg==\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.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import galkin\n",
    "import galkin.plots\n",
    "# import auxiliary modules\n",
    "import galkin.processdata   # routines to process kinematic data\n",
    "import galkin.readparsFile  # routines to read and check input parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Galactic parameters\n",
    "R0=8.0# Galactocentric distance (kpc)\n",
    "V0=220.# local circular velocity (km/s)\n",
    "UsunINUSE=11.10# solar motion in the U-direction (km/s), e.g. from Schoenrich+ '10, MNRAS 403, 1829 (2010)\n",
    "VsunINUSE=12.24# solar motion in the V-direction (km/s), e.g. from Schoenrich+ '10, MNRAS 403, 1829 (2010)\n",
    "WsunINUSE=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": 8,
   "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": 9,
   "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": 10,
   "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) = ( 11.1 , 12.24 , 7.25 ) 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": 11,
   "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": 12,
   "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": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100.35627126410797 396.87928478818867 3.4805916598031352 17.55446098248422\n"
     ]
    }
   ],
   "source": [
    "print min(vecvRp),max(vecvRp),min(vecRp),max(vecRp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "R_new = np.linspace(min(vecRp),max(vecRp),len(vecRp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD8CAYAAACCRVh7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xt8HPV57/HPs7u6WLaM5Au28Z2YkDhAChbYJCklgRJwKdDmwu20kJAaXiU9adPzShwgadqkbUhOw3FPEgIN9JAGgksaLi13AqnTBhtsQmwwBozxTZHx3ZYty9LuPuePnV2NVjvyyrvSrqzv+/WSvTvz29lHo9nfM7/LzJq7IyIiUkis0gGIiEj1UpIQEZFIShIiIhJJSUJERCIpSYiISCQlCRERiVSWJGFmf2Fmr5rZK2b2YzOrN7PZZrbCzNab2VIzqw3K1gXP1wfrZ5UjBhERKb+Sk4SZTQX+J9Di7qcAceAK4FbgNnefA+wBrgtech2wJ1h+W1BORESqULm6mxLAKDNLAA1AG/AR4CfB+nuAy4LHlwbPCdafZ2ZWpjhERKSMEqVuwN1bzex/A5uBQ8BTwCpgr7sng2JbganB46nAluC1STPbB4wHdka9x4QJE3zWrFmlhioiMqKsWrVqp7tPLGUbJScJM2sm0zqYDewFHgAuLMN2FwGLAGbMmMHKlStL3aSIyIhiZptK3UY5upvOB9529x3u3g38FPgg0BR0PwFMA1qDx63AdIBg/XHArvyNuvud7t7i7i0TJ5aUCEVE5CiVI0lsBhaYWUMwtnAesBZ4Dvh4UOYa4OHg8SPBc4L1z7ruMigiUpVKThLuvoLMAPRLwJpgm3cCXwQ+b2bryYw53BW85C5gfLD888DiUmMQEZHBYcPhJL6lpcU1JiEiMjBmtsrdW0rZhq64FhGRSEoSIiISSUlCREQiKUmIiEgkJQkpm8vveJ7L73i+0mGISBkpSYiISCQlCRERiaQkISIikZQkREQkkpKEiIhEUpIQEZFIShIiIhJJSUJERCIpSYiISCQlCRERiaQkMcLpVhoi0h8lCRERiaQkISIikZQkREQkkpKEiIhEUpIQEZFIShIiIhJJSUJERCIpSYiISKSyJAkzazKzn5jZOjN7zczONrNxZva0mb0Z/N8clDUz+0czW29mq83sjHLEICIi5VeulsQS4Al3fw/wfuA1YDHwM3c/CfhZ8BzgIuCk4GcRcHuZYhARkTIrOUmY2XHAOcBdAO7e5e57gUuBe4Ji9wCXBY8vBX7oGcuBJjObUmocIiJSfuVoScwGdgD/bGa/MrMfmNloYJK7twVltgGTgsdTgS2h128NlvViZovMbKWZrdyxY0cZwhQRkYEqR5JIAGcAt7v76cBBerqWAHB3B3wgG3X3O929xd1bJk6cWIYwRURkoMqRJLYCW919RfD8J2SSxjvZbqTg/+3B+lZgeuj104JlUmG6I6yI5Cs5Sbj7NmCLmZ0cLDoPWAs8AlwTLLsGeDh4/Ajwx8EspwXAvlC3lIiIVJFEmbbzZ8C9ZlYLbAA+RSYB/auZXQdsAj4ZlH0MWAisBzqCsiIiUoXKkiTc/WWgpcCq8wqUdeDGcryviIgMLl1xLSIikZQkREQkkpKEiIhEUpIQEZFIShIiIhJJSUJERCIpSYiISCQlCRERiaQkISIikZQkREQkkpKEiIhEUpIQEZFIShIiIhJJSUJERCIpSYiISCQlCRERiaQkISIikZQkREQkkpKEiIhEUpIQEZFIShIiIhJJSUJERCIpSYiISKSyJQkzi5vZr8zsP4Lns81shZmtN7OlZlYbLK8Lnq8P1s8qVwwiIlJe5WxJfA54LfT8VuA2d58D7AGuC5ZfB+wJlt8WlBMRkSpUliRhZtOA3wN+EDw34CPAT4Ii9wCXBY8vDZ4TrD8vKC8V1t7ZTeveQ6zatKfSoYhIlUiUaTv/B/gC0Bg8Hw/sdfdk8HwrMDV4PBXYAuDuSTPbF5TfWaZYRpTL73gegKXXn31Ur8tauGQZa9vaAfjE93/JeyY38tjnzilPkCIybJXckjCzi4Ht7r6qDPGEt7vIzFaa2codO3aUc9NSwP7OZO5x2ns/F5GRqxzdTR8ELjGzjcD9ZLqZlgBNZpZtqUwDWoPHrcB0gGD9ccCu/I26+53u3uLuLRMnTixDmNKfJVecTizo9KuvibHkitMrG5CIVIWSk4S7f8ndp7n7LOAK4Fl3vxp4Dvh4UOwa4OHg8SPBc4L1z7q7lxqHlGbezGbeM7mRac2juPczC5g3s7nSIYlIFRjM6yS+CHzezNaTGXO4K1h+FzA+WP55YPEgxiAD0Fhfw9SmUUoQIpJTroFrANz958DPg8cbgLMKlOkEPlHO9xURkcGhK65FRCSSkoSIiERSkhjmdAGciAymso5JyNDIXgjX3tnd5wK4xvqaAV9YN9D3Hazti0j1UUtiGNMFcJV1+R3P97lyXeRYo5bEMJQ9k1+1aQ+f+P4vSXvPBXCavioi5aSWxDDW3wVwC5cs40O3PquxChEpiVoSw1xjfQ2N9TXMm9k8JGMVa9v2c/kdz2tcQmSEUEuiCpTat52tuLNKHasoNh71yYsc+9SSOAak0k7r3kO5MYnBGqto7+zmcDJNe2d3GaIe/to7u9nfmWTVpj0aC5JjlpJEFVjx9u6iy+afua9p3UdHV4qOrkO9upbeM7mR/Z3JsiWI8PdNrNvW3musYyS0JvJ/x6guPdAU4SiaQj08qbtpmEume26gG+5aKvfN+vK7sJZv6HN39xFF049lpFBLYpjJPwsLn+GXq2upUDfKkitOz3VhATQ31Jb0HsNN/n4/2i49nU3LcKOWxDDXWF9DQ228bN8DsaZ1H2vb2tm6J9N9le1WmjezmRnjGnLlbnloDe/s76R176ERMUaRP0iv798YON1CZnhSS6KKHO0A6KlTjyvpzHRN6z4+dOuzjK1P9Om+Wr5hVy6m/HUbd3UAEDN69clXm8E6ew9PPy7WUAx2F/v7DuVYUn9jOINJLbbSKUmUydF+4MJn4R+7/ZfMnTLwD0+xH4RwjGvb9tNxOIk7pIGOrkMA1CV6GpcGLDhxfO752PqewyURs1zSSDts3XOIac1UbaIYSlHHQjEV5UAqteHUdVVoDKfQsbKmdR/JtDNn4mgdS1VCSaJMVm4sfoZSWCrvi1vXtrUTtzIEdKT3TTtph/zvjT2cTOce56/LfmgN+JtLT+GWh9bkxij2dyZZt62d90xuHLSYh7v1Ow7mHqcdXmtr56zZ4/qUK6byX7VpD617D/VK3AM1GMklKvb8MZxEzHLTtr/5xDogk0Q7ulIAuWOpGmaMDadkPBiUJMokv7Kvlm2FLb3+bC6/4/leH8awmEE8ZnSHAgh3N2U5cNX8GXz90bUcTqZJhVoURzvLZ+GSZSVN2e3vg1yOLp6BbCOqMglPMgAYVRuPLJt/gWR+LOHtLFyyrNdZd/ZsfOGSZTz2uXP6jXWofPOJddTXxEmmnclj63NdlR+7/ZcANNTGC87UU2ui8pQkjkGFKszssrVt+zkQUZFnWwHh2VLh7qal15/NrMWPApkzw3jMqEvEeiWcbfs7B1QZZ5PW0d5GJKoi7e8WJcVUnOHtlrs/vSZoKh5OprlvxWaumj+jYFfg2rb9BV8fbu0BvLH9AKNq4kCmhZj9e6xta+fcbz1HMu2MrU+UFHPU32IgZ9nxmHG4O8Vv9nX2WdfRlerV1VlNN6wc6RdNKkmUyfwC3QbFWL11L4e6ez70o2pinDatqVeZbGUxd8rYXs9TaSeZ9j4Hb3tnNzsPHM5V6LPGN9BQG2fjrg6Sae/TjRTWWF/D3CmZC/HG1if45hPrchXAybc8nisXPgM0erqmkinnY7f/kslj69h5sItxDTU01GYGxKM+9MX2VxeS/QBHVYKlXM+Q7UIMneDmuolioS7BU7/6ZMHXdxzuea/s2A/Qq6V204NruPnBNQX/Ju1FxppMOe2pwmWzZ+yQ+VvFY4X7MrPH1kDlV6BRFWq2Fbu2bT/Tm0exblt7r/0KMLGxjrH1CX7vtBNYcOL4gsfKUF+4WakB97BKd3OZ+yD1bZRRS0uLr1y5stJhDIpVm/bkKtyYwQM3fKDPhyPcFXPtP79QsPLIVgDFVixR6hKx3JnquIYaGuoSJGLG3kPd7O0ofarrrPENfKJlOg+s3EIy7fzpuXN6jW0ATB5bx7jRtbnf+csPreGN7QcYU5dg4alT+NgZ05g3s7lP940B7w0N/K9t209XMt3rzHvW+AYmja2PjC9bmWVfD73PziGzrw939+2u82CMp1o/UUbmGGuo631u2HE42WdZMfL3S/jYgf6TUirtdKfSuYSZ/7eLEtW6KsbRJMLWvYfYuudQ7vm05lFMbRp11DEcjVKShJmtcveWUt5fSaIfA2lK55/h5J/9R8nvX87/YOV/EPsTt8Ebz6g2c6c0sn7HQbryul4gU+HUJGIYmWm7qVAGOr6xltkTxgCF/0ZLrz87l5SzM7j+9Nw53PbM6xhw+/9oYd7M5j6th/y/UwywUP1YVxMv+u+Y/7uMqU9woDPZK/nkV8gxeloq/W0rnCTyW6YDdTiZ7rX/4zHrta9rE7FeXUgQ3bqC3sf+0bZs+nM0lW3+gPtwuyamHElC3U1lkv9By34YCn0As91EiZjRler90T6aiiQrv/l+LAsn1nwOBZMHwPb2Lra376apIZEbm3nx7d1YkGCzXXRhNz24Jvf4j+5aQTxmdBxOkvJMYm6o6319CcAJBc44X96yt894QrbiLjTTLLs+VaCL0MkkBiczAJ6NCXpOFMLdgJBJnKdP7+nKLLUbI9wKrq+J8ZWL35drFUZVqOGTqfyz9HGja3P7rNJdLFnZiybLeR+04UZJoh/h/tVr//kFANZ89aN9yt23YjPptHOoO0UsZrz7+DG5Smx686heTej2zu5cf2zS4LhRCfZ0lOe+PyMoR5Rsb2ife+6fIzvUlcqdjR/sTBKPx5jenKnYNuw8SHcyTU2i76Dr5Xc8T23o7D//rD6/WwwySWB0fYLjRtXkrmPJvnbOxNFsCSrYqLPu8LFG8JpyVr7zZjZTl4iRTDtfufh9XDV/Bj9avrHfCjX8/sPlmxWP5qLJY0nJScLMpgM/BCaR+ajd6e5LzGwcsBSYBWwEPunue8zMgCXAQqADuNbdXyo1jnJbtWlPrqK/8p+Wk4hZrimc7Y7YfbCLVNp7fbjTKe91lru2rb3XBXL7O5O5D23aYeGpJ3Dvis0lROr0nDOGzx0HcrFFuIYMbyu8bggu3hgGnMyZ+sHOJGkylfvatvbc/W0ciJnlTirCerpiMi2Dxpo0Tow0scKtiGCXj61P5FobiZjxN5ee0mtGVJ+KP9Wd6cuJ17Bq3Vss33yAR1/bT2NdIrPcev8tS7kQNHvs3/LQGn60fGMucWWvfcg3tXszJ6S2UOPdbEnM4nfHJ2nv7OIv/6DvWBwAu9+Gw+1w/HuhuwN2vAFNM6C2AV6+D951HkyYc1Tx9yvZBYmRdX+yKCWPSZjZFGCKu79kZo3AKuAy4Fpgt7t/w8wWA83u/kUzWwj8GZkkMR9Y4u7z+3uPwRqT6O/D8fbOA2xv78o9z575nTy5sd+ujkKmN9WxcNRa7tvzbto7U/SuhMMyw56NdHCABhzDcBwjW3k3cjC3riey/Io8/29aqOLPTyxpaknTlTtvyN9ufsyFtlFKIimUmKJiP5L8jpb+lhV6v/410EGSBF3U5Mpnjw+Ahrp4r21fwn/Szmi2Jo/jpeRswIiT4i8SP+GziYc57Al+4Jfw7a7LSBEPtuecaes4ra6N0RzinLkzWL7NWDDFmDcpBhv/GxrGwdQW8DSkk5nksO3XsP5ZSHbCmEnQ/hsAfpk+hRfT7+ZDtW8xb8ZYOPXjMO5dsL+Vbz+3mePSezgQG8umxGx2xifx3q41zO1azYc7nuTVuvfzq7ozOWQNHIyNoc4PY6TZsM94smMOTowYaRY2rGVOeiM1dDOnIZMsaryb8akdxElxUtc66jhceKcm6uHCv4dJp0DjFDhuGrz6IPzkU/3/MeqOgw/fBIk6OPF3oHMfWBwO789sY/TxsOM1OOEM2PE6NM/KrHt7GezbAnu3QMfOzL4bd2Lm/11vwqbn4cRzM8fGnPNh7FSYsSCzz4eRqhy4NrOHge8EP+e6e1uQSH7u7ieb2R3B4x8H5V/PlovaZiWSxJ53tvBGRwPZyiVOiu/UfIcX7DTu6fod0sSD5UlSvRpkfSujD9lqflR3K//Q/XHmxjbRRDtPp1uIk+KfUheHXgf5FVaMVJA6LHiexrFQ4gi/pue9T+UtdnIcbUwoEFOh5ymuan6Tnd01PHXgxNyaGCnSJPr8Tn231V8CCpeln3JR5QutK5S0wufz2fX5ZcIxh98nHSyN3q/XxJ/kt+wtvpT8DF0kSBPDyLQGmjnIjYmHmGOtnFi3jx3xSbxR+16mJjdzzqFnAViVPomru26imwQ1JPlU7TO8y7dwAjs5NfY230heyb2pj5AdbYjhfD1xF1clniONsY3xTGZXRLujR1v8BDptFPviTayuPYNNySZ+uOdUksRIkOb2+u9wPn1bOllJ4iTIjI29XvNeZiY3UO99K/dV6ZO4suvm4PdJ8S+1f8/8WKYF0Wl1dFNLnBRtiamkiXNS9zpW157Osw0Xsic+nrmHV7Mzfjxj0vu5uv1uEhTodo0l4AN/Bh274O1fZFoR01pg5xsw6xz4xT/AgW397o9+jRqXaaHUjILDBzItlGQXJA/B6IlwcEdP2Wsfg1kfPPr3qoCqSxJmNgtYBpwCbHb3pmC5AXvcvcnM/gP4hrv/V7DuZ8AX3X1l3rYWAYsAZsyYMW/Tpk1li7MYq1a/wifve4s0MWpI8tXaH3FV7BlWTfxDrt7y+7kP+r21f8fr6WncnbqQt/wEyFXgkK1cZtk2fl73l7ltb0xPwnC+l7yEpekPE1159shWRsWerRcz26VX+dD02/tWbObxV9q46JQpnDy5keUbdvEvz29k2/6Is8Aitl3aoHoxrZisIz3PX5ZNDb1bbD3rMu+XIMXS2q/xfHou305+nDRxjDSn2Vt8peZHACxPv5cFsdeYF3uzz2/wwJirebnuTHbvaKWtu5FtNVPZ0T0qN8hcG1SQXUFLIhxHA110WQ1OjDF00EWCRF0DCZKcPqmGNDFefaeDJHHmTh5Dijhu8dx75w8Qz2qq4fPx+5mQ2slTDRdzKDaKgzs3cxpv0kkdx7Ob52jhdWaxiyZG0Ukz+2mkg7EcoIYkXdTwUvok9nX1zF5qqI0xJtZFihjJIImGTWU722mmm54xuux4Sl26k3HpnUxOtrJ4fh08sRiaZsK1j0LT9D77M6f7EKx7FOrGwpoHMi2o/a1wysfhtX/P7MPJp8H6Z2D0hEzFv/83cNaiTMtg3OyeLrgCXXGkU5lEsb8VJpwMdWOiY6lCVZUkzGwM8J/A37r7T81sbzZJBOv3uHtzsUkirFJTYBcuWcbBQ4f59vGPMX7z46xnOuf/1VOs+tVKlr+2iQXN+5g3LuiSqmlg1StrWf7mb3gzfQIPpX87t50bxr3M4o5vwZjJ/P6BxWzkBKC46a35Ux2LVZuIkUqlSXlxHTVzpxz5SuR33/J4wVlDMaBlVjMvbCx8C+gbzjmRO5ZtyL1/TdxorEtw4HCSrmEyZ9dwrow/x3+mT6PVJwRLnRgpLrCVPO1nksKI43w9cTdXJp7jF+lTuS91Hoe8jl/4+/tci5Ad18hqqI33GeOCwlNJs/KnsRYaxA5Ps86/Y292PONoxiXKcQ1B5ED6rregYTyMaiq8XopSNVNgzawG+DfgXnf/abD4HTObEupu2h4sbwXCpwbTgmVVJzer4dO3cfItF5BMO197YQtXzT+TeWec2af8vHkw7+1fQMM4Vt/237QygU+dczKLL/o94GYAGu5cztzQa/K/ujRmkIjHSCbTmGWuQg1/ECFTmUD0dNnauPHjP1mQGzxcev3ZvWaS5J8rj6qNF3UV6Vd//329poNmXTl/Bic0jYpMEncs28DM4CK6R1f/hsb6mj7Xh4TlJ7W6RIyJjXUkU+mjbs0cSWN9os+FiL3bKsZ9qY/0iTRNnCd8QW5JCrg5+Wl+4e/nV+mT2E70jJhEItYr6WangIanysaMou6Imp1W3d7Z3ads9jtH+ru7ajmuISjr7KTx7yrPdqRk5ZjdZMBdwGvu/u3QqkeAa4BvBP8/HFr+WTO7n8zA9b7+xiPKIeq2CUeSnXd+8i2P5z60Nz24hr/+91epjTizC15JO1MBuPeFrdz7wtZea8PTaMNzzWsTMX78J73nlud/EE+ckPmQt+491Gta5PGNmZkYDnw/uOArLH++9x/dtSJXaWzJS0JR9nR09VlWXxPjD8+YBvR0KxUaPUimnRs/PIdlb2T6ePu7Pcbvzp3EhMY6DPjv9Tt73dLjxC89WvbrQRpq48ydMjb3feFZiXjmZoceKpe9YKynXO/uicyZehPf+9xfH/F9oyrZy+94njWt+xg3urbXXVKj9Hf31Kx4MDuv2HtiFUPXEIwM5WhJfBD4I2CNmb0cLLuJTHL4VzO7DtgEfDJY9xiZmU3ryUyBPcL0hcrJdg/kn7En086RJsc1FnkL53kzm3NnefkJIru+0Acxv4K5vUBiyL+PTrg7KVxpQHFJInyzP8gkpvD7huN8fVt77p5E2QoQes5YV23aw9U/WE53Mk08HqM7mc6dtT+19p0+7331D5Zz72cW8PXLTi3YmjkaxzfWMmFMHVv2HGJt234O5f2du/O6wpJpz+23WeMbet0XyYAJjbVMHFNX9L19+qtk4zEr+jvKS7n3ValG+jUEI0HJSSIYW4gaTT2vQHkHbiz1fQei0AVwA3Hfis29KqavBfPUyyVb8UR90AqNFfRXwRRzZ9Wjue3BvJnNuQHxhto4L9z8u73WhyuMeTObc4PehW7WNm9mM/d+ZkFuPcDn7v8Vh5MpdrT3bbF0J9Ms37CLGz+cmRP/5Ydf6XULiJgFc5QcJgatqvAU5qaGGmqClkH4/k/Z/QV9b5ESVqjbZ9Oujlwrw4G6RHHddmGFKtnw/aOyz/tTTLfPSP9OBDl6uuK6CFfNn8Ff//urmTGJMicIOPok1t9Z3GCdXY4O+u6jbtwWlk0Wxa7/ry9+pHcLI9bT3VOT6LlteXb/hxP3ot8+kcZRNbmElF9x3nXNmZGxhCvOc7/1XK6FUF8TY/LY+sg72M79yhO5VubR9slHVdoDqcyL6fYZrOSgpHPsU5IoUm0iRi2UPUGUor8KZiCDiqm0576g/kiV3NwpY3nx7d0cTqYH5f76+S2MLz+0pmDll/07ZKfq5v9djtTSgsL7b9LYera3H+4zHlDo95w1voE3th9gXENNwe6+oaRuHxksShLHqGIHFbODtR1dh4r6Up41rftyt6PIjhMMRqLIbrO/yu+q+TP6TdpRr+3vS2TyE0fUFNFw19T29i6+/NCaXi21oT7D1hm9DBYliWNYMWeXhb4ysj/h8tlxgsE8ex1o5Zf/DW/5y47mS2QKxZC/n/RVm3KsUpIY4eZMHN3r60qzs5CKKR+Px/rMeKomhQbnj3asJj9RDJc7mIqUSklihMteaJXtgz9SRRee9kkVfmFVOWYCFUPXCMhIoSRRpFKn0VazgczJD3c3pdI+6N1N5TZvZjMP3PCByKm5A6HBYhkJlCSOYYMxmPnjP1mQu0I8PC11ODnS1NxiabBYRgIlCRmQ8BXiw+37fkVk4JQkRqjsWfDCJcs4nEzT3tld9GtPnXocUPjaARE5tvR3lzo5hl1+x/MsXLKMtW3tua/gXLhkWaXDEpEqoyQxghWa6y8iEqYkMUItvf5sllxxeu57mYu5RkJERh6NSYxgmusvIkeiJDHCaa6/iPRHSUIGTNcHiIwcGpMQEZFIShIiIhJJSUJERCIpSYiISCQlCRERiaQkISIikZQkREQkUsWShJldaGavm9l6M1tcqThERCRaRZKEmcWB7wIXAXOBK81sbiViERGRaJVqSZwFrHf3De7eBdwPXFqhWEREJEKlksRUYEvo+dZgmYiIVJGqHbg2s0VmttLMVu7YsaPS4YiIjEiVShKtwPTQ82nBshx3v9PdW9y9ZeLEiUManIiIZFQqSbwInGRms82sFrgCeKRCsYiISISK3Crc3ZNm9lngSSAO3O3ur1YiFhERiVax75Nw98eAxyr1/iIicmRVO3AtIiKVpyQhIiKRlCRERCSSkoSIiERSkhARkUhKEiIiEklJQkREIilJiIhIJCUJERGJpCQhIiKRlCRERCSSkoSIiERSkhARkUhKEiIiEklJQkREIilJiIhIJCUJERGJpCQhIiKRlCRERCSSkoSIiERSkhARkUhKEiIiEklJQkREIilJiIhIpJKShJl9y8zWmdlqM3vQzJpC675kZuvN7HUz+2ho+YXBsvVmtriU9xcRkcFVakviaeAUdz8NeAP4EoCZzQWuAN4HXAh8z8ziZhYHvgtcBMwFrgzKiohIFSopSbj7U+6eDJ4uB6YFjy8F7nf3w+7+NrAeOCv4We/uG9y9C7g/KCsiIlWonGMSnwYeDx5PBbaE1m0NlkUt78PMFpnZSjNbuWPHjjKGKSIixUocqYCZPQNMLrDqZnd/OChzM5AE7i1XYO5+J3AnQEtLi5druyIiUrwjJgl3P7+/9WZ2LXAxcJ67ZyvzVmB6qNi0YBn9LBcRkSpT6uymC4EvAJe4e0do1SPAFWZWZ2azgZOAF4AXgZPMbLaZ1ZIZ3H6klBhERGTwHLElcQTfAeqAp80MYLm73+Dur5rZvwJryXRD3ejuKQAz+yzwJBAH7nb3V0uMQUREBklJScLd5/Sz7m+Bvy2w/DHgsVLeV0REhoauuBYRkUhKEiIiEklJQkREIilJiIhIJCUJERGJpCQhIiKRlCRERCSSkoSIiERSkhARkUhKEiIiEklJQkREIilJiIhIJCUJERGJpCQhIiKRlCRERCSSkoSIiERSkhARkUhKEiIiEklJQkREIilJiIhIJCUJERGJpCQhIiKREpUOQI4dS68/u9IhiEiZlaVUYd94AAAGnklEQVQlYWZ/aWZuZhOC52Zm/2hm681stZmdESp7jZm9GfxcU473FxGRwVFyS8LMpgMXAJtDiy8CTgp+5gO3A/PNbBzwV0AL4MAqM3vE3feUGoeIiJRfOVoStwFfIFPpZ10K/NAzlgNNZjYF+CjwtLvvDhLD08CFZYhBREQGQUlJwswuBVrd/dd5q6YCW0LPtwbLopaLiEgVOmJ3k5k9A0wusOpm4CYyXU1lZ2aLgEUAM2bMGIy3EBGRIzhiknD38wstN7NTgdnAr80MYBrwkpmdBbQC00PFpwXLWoFz85b/POJ97wTuBGhpafFCZUREZHAddXeTu69x9+PdfZa7zyLTdXSGu28DHgH+OJjltADY5+5twJPABWbWbGbNZFohT5b+a4iIyGAYrOskHgMWAuuBDuBTAO6+28y+BrwYlPsbd989SDGIiEiJypYkgtZE9rEDN0aUuxu4u1zvKyIig0e35RARkUhKEiIiEklJQkREIilJiIhIJCUJERGJpCQhIiKRlCRERCSSvnRohNMXBYlIf9SSEBGRSEoSIiISSUlCREQiKUmIiEgkJQkREYmkJCEiIpGUJEREJJKShIiIRFKSEBGRSJb5ErnqZmY7gE0DeMkEYOcghTOYFPfQUtxDS3EPrQnAaHefWMpGhkWSGCgzW+nuLZWOY6AU99BS3ENLcQ+tcsWt7iYREYmkJCEiIpGO1SRxZ6UDOEqKe2gp7qGluIdWWeI+JsckRESkPI7VloSIiJTBsE0SZjbdzJ4zs7Vm9qqZfa5AmXPNbJ+ZvRz8fKUSseYzs41mtiaIaWWB9WZm/2hm681stZmdUYk482I6ObQfXzaz/Wb253llqmZ/m9ndZrbdzF4JLRtnZk+b2ZvB/80Rr70mKPOmmV1T4Zi/ZWbrguPgQTNrinhtv8fUYIqI+6tm1ho6FhZGvPZCM3s9ONYXD13UkXEvDcW80cxejnhtJfd3wbpv0I5vdx+WP8AU4IzgcSPwBjA3r8y5wH9UOtYCsW8EJvSzfiHwOGDAAmBFpWPOiy8ObANmVuv+Bs4BzgBeCS37JrA4eLwYuLXA68YBG4L/m4PHzRWM+QIgETy+tVDMxRxTFYj7q8D/KuI4egs4EagFfp3/GR7quPPW/wPwlSrc3wXrvsE6vodtS8Ld29z9peBxO/AaMLWyUZXNpcAPPWM50GRmUyodVMh5wFvuPpALHIeUuy8DductvhS4J3h8D3BZgZd+FHja3Xe7+x7gaeDCQQs0pFDM7v6UuyeDp8uBaUMRy0BE7OtinAWsd/cN7t4F3E/mbzQk+ovbzAz4JPDjoYqnWP3UfYNyfA/bJBFmZrOA04EVBVafbWa/NrPHzex9QxpYNAeeMrNVZraowPqpwJbQ861UVwK8gugPTzXu76xJ7t4WPN4GTCpQppr3/afJtDALOdIxVQmfDbrJ7o7o+qjmff3bwDvu/mbE+qrY33l136Ac38M+SZjZGODfgD939/15q18i0yXyfuD/Ag8NdXwRPuTuZwAXATea2TmVDqhYZlYLXAI8UGB1te7vPjzT9h42U/vM7GYgCdwbUaTajqnbgXcBvwW0kem6GU6upP9WRMX3d391XzmP72GdJMyshsxOutfdf5q/3t33u/uB4PFjQI2ZTRjiMPtw99bg/+3Ag2Sa3WGtwPTQ82nBsmpwEfCSu7+Tv6Ja93fIO9luu+D/7QXKVN2+N7NrgYuBq4MPfx9FHFNDyt3fcfeUu6eBf4qIp+r2NYCZJYA/BJZGlan0/o6o+wbl+B62SSLoM7wLeM3dvx1RZnJQDjM7i8zvu2vooiwY02gza8w+JjMw+UpesUeAPw5mOS0A9oWakZUWeYZVjfs7zyNAdjbHNcDDBco8CVxgZs1BF8kFwbKKMLMLgS8Al7h7R0SZYo6pIZU3hvYHFI7nReAkM5sdtFCvIPM3qrTzgXXuvrXQykrv737qvsE5visxOl+mEf4PkWlOrQZeDn4WAjcANwRlPgu8SmbWxHLgA1UQ94lBPL8OYrs5WB6O24Dvkpn5sQZoqXTcQVyjyVT6x4WWVeX+JpPI2oBuMv2u1wHjgZ8BbwLPAOOCsi3AD0Kv/TSwPvj5VIVjXk+mDzl7jH8/KHsC8Fh/x1SF4/6X4NhdTabympIfd/B8IZnZOW9VQ9zB8v+XPaZDZatpf0fVfYNyfOuKaxERiTRsu5tERGTwKUmIiEgkJQkREYmkJCEiIpGUJEREJJKShIiIRFKSEBGRSEoSIiIS6f8DSficmruK3PgAAAAASUVORK5CYII=\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(R_new,Vc_Model_I(vecRp,params))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXucXGWZ579P3brTTefehjYkhAiKUZSYHkiU8QO4oxhnjDcMwigyaOKIu7pzUcRZjbrOR3HVzey6kAgqzKAEQZRBBBnEjX5IwASRQAMSQzDp7ZAm6SRF+lpV7/5xzql+69TtdN375Pl+Pv3pc3nPqbdOVf3Oc573eZ9HjDEoiqIo4SXS7A4oiqIo9UWFXlEUJeSo0CuKooQcFXpFUZSQo0KvKIoSclToFUVRQo4KvaIoSshRoVcURQk5KvSKoighJ9bsDgDMnz/fLFmypNndUKpkz+BxAJZ2d+YsV3MeRVGKs3PnzheNMd3l2rWE0C9ZsoQdO3Y0uxtKlazdtA2ALetX5SxXcx5FUYojIs8HaaeuG0VRlJCjQq8oihJyVOgVRVFCjgq9oihKyFGhVxRFCTkq9IqiKCFHhV5RFCXkqNAriqKEnMBCLyJREfmdiNztrp8mIg+LyG4R2SIiCXd7m7u+292/pD5dVxRFUYIwFYv+k8BT1vrXgG8ZY04HhoAr3e1XAkPu9m+57ZRpwtpN27IzU+vRXlGUxhNI6EXkFOAdwA3uugAXAre7TW4C3uUur3HXcfe/xW2vKIqiNIGgFv3/BD4NZNz1ecARY0zKXd8PLHSXFwL7ANz9R932yglEcnSC/iMj7Hx+qNldUZQTnrJCLyJ/CRw0xuys5QuLyDoR2SEiOwYHB2t5aqXJrN64lb6BJPuHRrj4+odU7BWlyQSx6N8EvFNE9gK34rhsNgKzRcTLfnkK0O8u9wOLANz9s4BD/pMaYzYbY3qNMb3d3WWzbCrTiGOjqexyxsD2PXkfv6IoDaSs0BtjPmuMOcUYswS4BPilMeYy4EHgfW6zy4Gfust3ueu4+39pjDE17bXS0my8ZDkRd1SmPR5h5VL13ClKM6kmjv4zwN+JyG4cH/yN7vYbgXnu9r8Drq6ui8p0Y8Wpczjz5C5OmTODWz6ykhWnzml2lxTlhGZKhUeMMb8CfuUu7wHOKdBmFLi4Bn1TpjFd7XG62uMq8orSAujMWEVRlJCjQq8oihJyVOgVRVFCjgq9oihKyFGhVxRFCTkq9IqiKCFHhV5RFCXkqNAriqKEHBV6RVGUkKNCryiKEnJU6BVFUUKOCr2iKErIUaFXFEUJOSr0iqIoIUeFXlEUJeSo0CuKooQcFXpFUZSQU1boRaRdRB4Rkd+LyJMi8kV3+/dF5DkRecz9O9vdLiLyLyKyW0QeF5E31PtNKIqiKMUJUkpwDLjQGPOSiMSB34jIz919/2iMud3X/u3AGe7fucB17n9FURSlCZS16I3DS+5q3P0zJQ5ZA9zsHrcdmC0iPdV3VVEURamEQD56EYmKyGPAQeB+Y8zD7q6vuO6Zb4lIm7ttIbDPOny/u81/znUiskNEdgwODlbxFhRFUZRSBBJ6Y0zaGHM2cApwjoi8FvgscCbwZ8Bc4DNTeWFjzGZjTK8xpre7u3uK3VamM2s3bWPtpm3N7oainDBMKerGGHMEeBC4yBgz4LpnxoDvAee4zfqBRdZhp7jbFEVRlCYQJOqmW0Rmu8szgL8Anvb87iIiwLuAJ9xD7gI+5EbfrASOGmMG6tJ7peYkRyfoPzLCzueHmt2V0KJPNEqjCRJ10wPcJCJRnBvDbcaYu0XklyLSDQjwGPAxt/09wGpgNzAMXFH7biu1xBOd5OgEfQNJAC6+/iHOPLmLrvY4W9avamb3FEWpkrJCb4x5HFheYPuFRdob4Krqu6Y0mmOjqexyxjjrXe3xJvZIUZRaEMSiV0KOZ7HvfH6Ii69/iIyB9niEjZcsZ8Wpc5rcO0VRqkWFXsmy4tQ5nHlyF8dGUyryihIiVOiVHLra43S1x1XkFSVEaFIzpSwaJaIo0xsVekVRlJCjQq8oihJyVOgVRVFCjgq9oihKyFGhV8pSj7QIfQPHdIBXURqEhlcqRVm7aZumRVCUEKAWvVKSQmkRFEWZXqhFH0I8l0ilVrfnVtmyflXVaRGC9KXa/iqKUhoVeqUkmhZBUaY/KvRKWTQtgqJMb9RHryiKEnLUog8hO/YebnYXFEVpIcoKvYi0A1uBNrf97caYL4jIacCtwDxgJ/BBY8y4iLQBNwMrgEPAWmPM3jr1X2kAjRok1bh6RakPQVw3Y8CFxpjXA2cDF7m1YL8GfMsYczowBFzptr8SGHK3f8ttpyhZkqMTjKUyJEcnmt2VpqB1eYujmVLrQ5BSggZ4yV2Nu38GuBC41N1+E7ABuA5Y4y4D3A78bxER9zxKA0hXeaXTGZMVoloPwK7euDU7AevpA8nQi51ftIpNQAMNL1XqRyAfvVsYfCdwOvBt4I/AEWOMN3tmP7DQXV4I7AMwxqRE5CiOe+fFGvZbqRF+IdrVf5Th8TTD4yN5QgTVi5F/Atb2PYeqOt90Q+vyKs0gkNAbY9LA2SIyG7gTOLPaFxaRdcA6gMWLF1d7OqVGpDKTjwP1EKKNlyzPTsCKCKxcOo+tfxis2flbDf+NUevyKs1gSlE3xpgjIvIgsAqYLSIx16o/Beh3m/UDi4D9IhIDZuEMyvrPtRnYDNDb26tunSbhFyLbtVIPIfImYO0ePM7p3Z0nnMhVMgFNZw4r1VJ2MFZEul1LHhGZAfwF8BTwIPA+t9nlwE/d5bvcddz9v1T//PShqz1ORyLKKXNmcMtHVtZFiLva46TTGfYNjdT83NOBrvY4C2fPOOFuckrzCGLR9wA3uX76CHCbMeZuEekDbhWR/w78DrjRbX8j8K8ishs4DFxSh34rdSQaERUiRQkRQaJuHgeWF9i+BzinwPZR4OKa9E5RFEWpGk2BoNQVjRlXlOajQq/Ujd8+d5i+gST7h5xQTRV7RWkOKvRK3bBH4E/EmHkPne2pNBtNaqbUjL6BYznrMxJRhsfTgBOquXLpvOw+Y2AslWl5K19DG5UwoBa9kseynpk1EbZoRAqGaiZHJ8gA46kMF1//0Amb80ZRGoUKvVJXCoVqah3a1qKVXEs6eF8f1HWj5FCpJb920zaGx1J0tJX/Ss1sn2zTHo/krCv1o1o3VL1vBqUSvtWLE8Ulpxa9UlOGx1JlBaGrPU4ESMQi3PKRlQD0HxlRF84Jjj7p1Q81pUJMPdIMl8IYJ9ImOTpR1hITgbaYY2d4VlxEaIgVp1RGva3fqSR800HyqaFC38JU8qhsW8WVPv4G+fF4ffMibdIZQ8bd1zeQpCMRZWwindOmEHbIpabtLU9YBK7Q+6gk4ZsSDBX6kNGsfOd2euNC68VYuXQeESGbtlj99fkkRyc4NpoqOUAZlhtAV3ucrvZ4VuSnw/uaDn3UX1ULU8kXZ+fzQ7z3uoeA+uY79/rmfcntgbSIwOndndnslMt6ZuYdnzaOP9+24ma2x05oa77QE5xXCAbgvdc9RATobI/ltG1lgakX9s1PLf/yqNCHDPtL/+FVSxr2I7AF2nMXpTPDpDKm6CCrV/LQs+KUfPxPRpXm+167aRt9A8cK3nRbnb6BYzk3tmLROSfiDS8oKvQtzlkb7ptS+/FUJrt8/dY9fO+hvSRiUwuu2rXhbYHbev53W0C62uM5lmjfQDKQ5eVZac2w7Fvh8bvQa9uFYCIC7fFozSa0TQfsQvLed0LLMU4dFfoWZ3hsaiFm/sLgY6kMqXSmcOMakM6Ygla73xLdvudQUaHvGzhGOmOyNwYv+kaZLAQztzPBzPZY6Iu1eIbD2k3bciz3voEky3ocy90uR6nlGIOhcfQK4Nwg/DcJj2IzJ1dv3MrweJrxVCb7g/Q4vbszu+zVhi1FoVq1U6WVZnjWkmhEODoyEXqR9+P/Dnjr3rhOPaughQ216FucIDNNbY6PpbCN6YgEO4f35DCV1yslxp4lmsqYQLVhU6ncp45Wib6px6BfM9wutgukGQR1jXkuwC3rV+XE1Qvw8fNP59JzFwP50TlKacr+mkRkEXAzsABnLGizMWajiGwAPgoMuk2vMcbc4x7zWeBKIA38F2PM1BzNSsWIkDNiJ1K/19p4yfJshI8XImkTjQjRiATyn4pAR9y5McQCHlMLij0B+F0H3qDfPZ98c0P6VQn+QUsPvwtk9catOdd3V/9RUhnD6o1bW+b9ee8jHo0wlspggH/6yS7+bfteutrj03ZguVkEMZtSwN8bYx4VkS5gp4jc7+77ljHmf9iNRWQZTp3Y1wAvB/5DRF5pjEnXsuNKMXxKT/2UfsWpc7KvFo9GAo8FbFm/iiVX/wyYzJEDkzcGG79wlbIIa219t/qUfPva9A0cY3gsVXBy2pjvaWn34PHsrGR7bKRvIMmyz9/LWQtnVdynWj+t+GsaeAOvJ9KAdC0IUjN2ABhwl5Mi8hSwsMQha4BbjTFjwHNukfBzgPA5TxvAVCJgAH7w8J+45s5d2fUvr3lt9nHXo9Bj9FSiTuzYee+H6BcTj0b8IOthfftdB96gXyU0Kua7oy1W0Mq1rw044yeeRd9/ZITh8Unff9CJbvXG+87YUUc68Fo5U3KEisgSnELhDwNvAj4hIh8CduBY/UM4N4Ht1mH7KXBjEJF1wDqAxYsX+3crFXLpuYuzQv/P7z4rT+QrpW/gGGdtuC9HSIpZuLv6j+ZYhf6biD8mGiYnUPnHCAqFbxaiHta3N+i3e/A4S+eXH2fwCBLzDdVbv1O5Ua/euJXdg8c5vbsz5wZo38y8iW6tZCnbUUcq8pUTOOpGRE4C7gA+ZYw5BlwHvAI4G8fi/8ZUXtgYs9kY02uM6e3u7p7KoUpAgoq8N5mmFMnRFElXPD0hKDZgWs4q3NV/NLvcN5DE+JoPj6XYsfdw4AiaLetXsfGS5Xhen2qsbz9d7XHaYpGKxwyqvQF5g6hp95qmMyYnX3upSCN7X1d7nFhE8lIp2BEsrZhQbsv6VZy1cFZeTQNlagSy6EUkjiPytxhjfgxgjHnB2v8d4G53tR9YZB1+irtNCRnFwv3s0MpClJrt6d1MIkz6nYNEApVKiGXPCm2EGwlyk73ZHD4+ztER5ymm3GQ4238+bm0fHh/hvdc9lJM4zuMVn3XGPuxQ2Vd89mfZdftYbzzEGyM57A6Ga1qBqTEd0jEEiboR4EbgKWPMN63tPa7/HuDdwBPu8l3AD0TkmziDsWcAj9S010rNsS3FUl/WYtV/IkCG0kO/ngjGIpIjXJ4m2ZZ9Bjjuiv5YKsMLx0ZLitDaTdvYNzTCsp6ZLfVji0aECM57nGGJayk84S3nLh8ZT+elRIhQONLKP0Q/Mp7OPgFN3gScm0a5rKf1vFkWO7cXUTSVMaRG0IxiKTD1zyCIRf8m4IPALhF5zN12DfABETkb5/uzF1gPYIx5UkRuA/pwInau0oib1mXtpm389rnDZHCsPX/uEL9bZ/9QvkUIZFMUGxx3TFvMicIpFAY3Ol746+AfzvXWx1MZ9h4aBnIHWsv5pYP84Hf1H2VkPM2MRDQv2iSIS6voea1B9FJjFDb+1/KLs5+C+wVedXIX+4ZGsk9Dy3pm5g3IzvBZ9Bkzea6MyY3MqUUYYzVWr/c9bEWmSzqGIFE3v6GwoXZPiWO+Anylin4pDaRYCNtZG+5jeCxVcMbs8Hi6ZA6dsVSm4ACQVxi8UjwRsn/4niAXCi8sJdTezF5w3o83MOwXlaBPO5Xi76Ptqjo+mip4vdpiERKxSEGRLjQO4J/AZouRPwUF5Ebm+Ad9pyq6hazevYeGs7H7QYRxKnHzjRxMnkqxlGbSGtMPTzAqTaBVq8RbO/Yezlm3LUMvJ3zfwLGsv7wY40VCKj0MjkjuefElBJh/UltNomFiESn4w2+LR5k1I54VZG8g03t//lmh/r54NxBbeD2h9T/tFIoG8j6Xcm4G//ZSn+vaTdtyEsQBLJnXwa/+8QLAuVkNvjTGi8lxDJNic+29T+eNS5y14b6iE9iiEaEjEc0+3dht/PH6U8UOvc0Y+MPBl0i51oNXpKacS8u7idfTsq/kdzVdiqWo0IeYoBaoN5s2EYtYllz1eVUMjqXsidTB5DizO6p7rI1AXuUq2xr1BhsLuTX6BpKc//UHWTCzHciPGjImPzlbsaedWlLOrXHWwlkkRyfYPXicWEQ4dHw8b+4AQDwqLJ3fybX3Pg1M/UkkGhEEsknqCr3PStw4/loFczviHExOjtLM7UywcPaMkufY1X80L4tlqzAd0jGo0Lc4haypUlZNkFKCnivG/9jshRHaPtFZM+L8ou8FasWR4epyrWQgq77HR1MYCvuxi/m29x4a5mByjLMWzsoWKffszYm04ekDSdrjk/763z53OHsu+9E8yNNVoc/Jvy1onL0X5ukN1PYNHMubpJZKGzcSaiR787OfRPx43yfvnMY418JLUtdhjVtU+xS5euPWrNXrvddS7g7/XATvRv70gaTmn68AFfqQ4R8csgfVIDfcz/sxZ6MurMdjTwRGJ1p3HL1SX/9EOpPjzkj6rtnweDrr9vBeIxGLFM2U6Bd97zoWs35tt08lg3nFBliFSYv70T9NRkYFjd+PxSI57rhazpLtao+zb2iEa+99mi3rV03J3TFdBjxbGRX6FmeqaQrsUoLC5KCaJy5zOxM5U97ndiYYHR8hg+Pj9oQiaIhfUBJRYbxYHuQ640+4NpGeTODlz/cfEcey9Wd6bItFsi4Rvy/eJjk6Qdp4E8zy3Qz201KxNAu28Nk+cX9fvScSg/NZe2GVydEJJqxrLThuKi9+38P7rL3XWDRnRl6qhHpZzeXcHfbrTpcBz1ZGhb7FsYXk4ecO523zUygNrReRAvkTeAaTY1mrdXg8nZe+oBJsd4hHs0QeCt+snhpIctaG+/Iiiry2/hz7Y6kMu/qP5gwaejN4PZZ9/t6cQVPv+KibKtqeTOWJbLWDeSKTIXFp45x71ozcm8v8rkSgcRc7MqeRGUTLMV0GPFsZFfomE8RKn0qkg/2Ya8h/ZLeFyvtR29jrlWizAO2JaI7gtSJOFE7wCKDxVCY7yatQSGDaOJOQCr5WmetYyrr1W92e28Z2PXkMj6VY1jOTT190Zk7+muv/ureqcYVWYDoMeLYyKvRNppyI24/6nkUf1HUTEQoOHnopgvu+dFFOe3CEuv/ICG3xaMGZl+XwIm3CjDe24eFdq0QsUjCLZ4b8J6lahgr6XTp2MrYgRV9syk3qUqYnKvQhw/5Rn3lyVzae2sP+AZ//9QdzBAscF8X+E6xkXSXYNzNvuViqZnCs/WIDtPXwg1ebjE0JFyr0LczO54fYvudQoNhhf/gaOJZnVHJdMI88N+lT9tIK1BqR8u6KsFIsbcGMRDRH5GuVYG3RnBnsHjxOKpWpazWxZqOhlNWhQt8E7AkyxfDnCS8UB223ffpAkoyBy27Ynrc/aguASJ4bodbMnhFjaLi1qjE1gohAe3xyfEKAuDUJzXPDTWUiUyFXiuf28YdXYnJdQq1UhWnL+lWs3ri1aFI8pb6o0DeYnc8PZX+cl92wnYjkl88D2L7nUDYCJGPIcaf4BcJuO+FzH/QumZuzvubshTkVqJbM6+BgcqymfvWeWTM4Mpycsn+/2XRUOYjsxeB7GJxB3GcOOJ/3WRvuy0mp4E8OV2rWqVfX9fTuzmycvj+6plFU4r+30zhcfP1DtMed9z2Vc7XKTWs6okLfYLbvOZRdnkhliEYj2R+6HX7nzyNjR8+897qHWNYzOYMyOTqRjRUvJ652Bap4VNh3eLhsdE1EnGidiSaGSNaak2e2cXh4InudozKZAqDW79K+bP5z2+mCiw3M2ykevJm74E/hYEgwzmIGWTnyGK+ceIq56RfhR/Nh8Bk42McWr+mG8n0u1XZL/qayfDvzTr7BxWSIIibNVZlbuSpyl1ukNCAF+tIKbIl3wn99onzDJqJC32BWLp2XXY5GI8SKJHNKl3Fy27MDu9rj2ThjJyHZ5OO8P/Z79cat2eWgwp0xkJmCyL/40hgzWjzEcjA5lhVgbxwjeLhlbvH1KGnWRn/JQjnE/02/jkfMq/G89REy/CjxJVZEngVgZ+IMLhu/hglixElxS+Kfs/uKUVQkh5zz3ZE+D0F4T/TXrOBZOGId/GTAt1RnVkaeIkGKCSBOipWRp5rdpdoxcRzGj0PH3PJtm4QKfR0IOqFpIpUhhTN46WUp9B7PX3xpLCfxk41/duDaTduyccZ+hsfTOcUodg8er+g9xaPBLfrB5DgzrGieCHBmTxeDL40xWOQ91R6vDEqujR4nza2JL+eI686MI75jxDHuMRHSvIwhDjAve463Rnbwy8xyUkSt82YQMtyePp8UEV8+byGDcEf6vOzrrYg8yy2Jf2Z75tWsjDxVVuShvEj+OP1mxolxR/rPA904Gk4kzsujx/kcP+TxyDIueW0nT/zpQn4TeQ+ffMc50D4L2mdCohPiHRCJQTQB0ThEouXPr5RFhb7B+Cc0gSMZ/sRNi+d2EItOZNO52vhzrtgzXwthJwLzf+CFZrEW4oyXncQfB4+XDCH0MDipFcYmRkgb6GyPZW9Ex8dyU+5GI/7BYf/7dQQ1SobTZIDdZiGTc0FNXrvc89gJApz/KSI5wguwPfNqxolhsu0zZIjwAnOtczlkCtw8DBFS4FjcZBCMtVe4PX0+743+JkfsA4tx58tYMXsm17xwN4+mT+eDr+/iF4+dx+28jTlz5nHf4BxGSQDCKBE+krmGFT0d3HDFGx2hjCYqEstapcT2+JRVznHF+1ZxrXv+T56hfvdGEKSU4CLgZmABzrd7szFmo4jMxXHXLcGpMPV+Y8yQW3pwI7AaGAY+bIx5tD7db03KTWiy83Z4g7HveN3L6Rt4BnBk5OLeRaxcOo87Ht3Pf/QdYP5JbVmXzDMHkjnW/FTwC3UQke9IRNk3NEIqncnOpi2Viz4iMLMtxqsiv+PG+NedjZ4vNgK0T7b9QeoCrsl8JLsuZLJWtS3QBliaFfpS+SrFfVe2yE+2MQWEd2XkKWKkGXdLpYjb0vhuFA9kljNZTkWy/9M5r5BfbmWcOJ/svJaFs2dkq3l1JKL0femiAu+jMD9zP+cV71rFtS9MivCvNm6FrKtOmDVrDjf87QVlz1drIVdamyAWfQr4e2PMoyLSBewUkfuBDwMPGGO+KiJXA1cDnwHejlMn9gzgXOA6979Cft6OD3/PKadr++7TBr794G6++Ytnsn5k241zzZ27ciJn6oknSGs3beO0iWf56ov/GSLwOfM33JJ+C7bF/eXYdxmiy3FJHHkWAgSFvCqynxgpUkSJuO4W27IGQciQIMX62N2cn3mMLekL+H9mLoPMwb4ZTFKottXk/nHivHf8iyyZ10EqY5jZHmPcGtcwOeeaFPRYLEEmlalosPbw8XEOHx/PySu09Oqf0WkNqNrlBwvhzyzq1cm1ny/2HhouWrVJRf3EJUgpwQFce8wYkxSRp4CFwBrgfLfZTcCvcIR+DXCzMcYA20Vktq+Q+AlPobwdK06dk5NQKkgR6dphuCz6gGvvCq+R53jSnIbB8F5+AxvW5kVavCf6a25NX4A7AoABhujiqthdgV7xutRfcWv6AgaiPaRcmcotPuhZ9BCPRbnlo+fx334yf7KABXZWytxrdfLMNg4cG8uu57uHHKY6YayY26qtSOoDj4jbB3+bWkT3LOuZye/2Hcl5wqpFGt9qarwGQW86jWVKPnoRWQIsBx4GFljifQDHtQPOTWCfddh+d5sKfRmKZY3sGzjG6EQ6ZzD05JltnDqvM68sYL5L3/DWyG85P/L7SWs78mx2ANKL/nhP9NeTkSGZM9gwfjkTxNiSvoAvx77HpbEHc866IvIsX459l8+nriCNIBj6E69gXeZzjC7o5ea/vTBrffrL7a3dtC2bt2dZ90k5UUK2aHoTvbx8+na7M3ucCWSFBq1j0UiOlZup8wSxTIkIqQjwZ6fNzaYktvMK+WfLlnLDlcqJFIsI4771QgQtYlOqGEql2DV9a5VHR28WwQks9CJyEnAH8CljzDGx5lsbY4yITOnXJCLrgHUAixcvnsqhJxzLembSf2QkZ9JULOqI3xmyn3sTn85u/+D4Z/h15nV4Lg3BsD72s7zBv1LRH3ekz2OcOI7/WbgmdSXfT72NP7GAUdoA6GqPkUxNDgCngR8M/5kjbBGnLFyQH6KXUz0Wi/DDj67My83jvX97/gHAO173crb+YZCu9jgvJp0qUKfO62DvoeG8XD31jv6f3REvGE3UFotw9qLZ2evgPbWNjKcRHCu/UO3ZUtgD796xY77iMF4VrUJPhUFeRwt9hI9AQi8icRyRv8UY82N38wueS0ZEeoCD7vZ+YJF1+CnuthyMMZuBzQC9vb3hmYlTBQWF8fgh+PpSxwJn0gLfePwLrBh5FhK5zd8eecQVes9vnWF75tUFozz2RE5le6SX7/F+xkhk48j9IYIgPMtiZ+KUu7VYGoUgH6SdQTM5OpG1eEuxcum8rKsmIs761j8MAk5a5FTGMDxe/7QL8agwrzOR4xrqPqktW5zbJpXO5N20ohEhIpPphj2CWqe2j95OYfzSaCrn9YvVYfWnVCj0uvUo9LF641Z2Dx5n0ZwZaok3gSBRNwLcCDxljPmmtesu4HLgq+7/n1rbPyEit+IMwh5V/3xp4kzAo/8Kd32iaJug8deXxh7kT+ZlbEr/lSP1EuVn3R9la/vHAeeHvezz92bj9bva4+zw5Tb//hXn5LgYElHh9JedlK1UZc/UtPFyvZTDnrT11ECS7q5Eni+4kAh6qXdjEcla/nZfSs07OHlmO88fGiZSxF8flIm0yRF5AbdOazDSGUPGFL9RlsNfocpj9catWXdLteJcrtDHVF0vtivo6QPJooPFpdCbQ3UEsejfBHwQ2CUij7nbrsER+NtE5ErgeeD97r57cEIrd+OEV15R0x5Pd4zhOy+sZWbmKGyAbOxMgDHMbPz1qefBX94C81+JP2Wh7f8GR3i9gtHgCIInjF5ede+H67kErr336WxPh/4mAAASuElEQVR8fXdXguv/ujdbRs9zI9nlCAFOmTODme2xQKLnn0twMDnOZTds55aPrMwRMvs9gePmWb5odnZ7sb4kYhFSqUx2ytQtH3FcQoeOj+eVy6sV/pnA9rCy1397roRXzQucYi9TEb9CA6Vedai5nYmaWOC1LPShrqDmEyTq5jf4n+YneUuB9ga4qsp+hZfhw47Il+JNn4I3/wO05Was9Nwde694R8Uv7684ZVeU6mib/DqIO5q5dP5JOT92f51TcERt0k1QXug3XrI854kBnFnC2/ccynktT9BmuhOu7D54+AuttMej2UFtb6DaK0jtCW7QSWJBMFDwCcebpGaXG/Qb8f6yg/YMZhv7c7Ffyx4o7Rs4xthEmqMjE9mbcilKDcYW2+9d96la11rztfnozNhG0zmPyxf8mHFp44cfOy/rLikXQ12KUo/Sw+PpnAyJM9tjOb5uuwC0l9P+0xedyYe/90jR2bZ2BaNK8qCvOHVOntjGY5HsXILk6ATPHnwpJ8rITuLmf99eWOrp3Z08cyCZjfTpXTKXvoFjeZEe3utGI5CugeLbs4Bhskbs8FiqqEhD/lzeWCySjTCysd1Y9hPMdLGOteZr81GhbwKjkY6KjotOUVA95nYmODri5NjxJ0DbNzSSl9v84usfAlPa6vUqGKUrVErvicHvbrDTONscG01xzyffnLPNs/iXzOugqz1OcnSirP/bLrxSqcjbeX88C/Xae5/OWu52ami/9bvs8/cyMp5mRiLKknkd2ToC7fFIXmqLQhSzju1Q1iAWd7mZsbWeOas1X5uLCn0IKCQmnuVoC5HX1kuCZrsz/H7UIPeUZT0zy+bZKYYder5w9oysAPjDKMFNqWDNIC1UdMO2jofH03z6ojMLFsR+xTX3FL0RBE1R7Im8AEvnT60m61kLZ2WT1wFTtnTVOlYqQYU+hHhumkSsvJXoCaDfUlw6f7IqUik8sZ8KO58fykkFYGf0tMMoBZjflaAzEcsbfPSPNfgF2u/v9/jymtfmpI84eWYbsWiEw8fHs/0JSkTIuk22rF/FWRvuy970vKcN/yCrXYDDGwyPulFEQa3nWljH9Z75qrQWKvRNZjyVIZUx/ODhP3HpubWbOBYVWL5oduAfcSMtRb/Vbou23Y/Dx8c5PpbOTka6+PqH+NHH3pg3IBwRWDy3gz8dHs6Jsy+Ed41//sQAb39tDz99zJni8ZvPXAhMxntPlMhp4yVn6GiL5cSle9WjHtt3JCfdgT1GkvI9TaQypi7pLsqlyi4383WqN2+ltVGhrxJvMHUqeFbfq/7p51lBuObOXXzx358kUWAwzsMb7Cv0muUGc4POUm2EH9UvwrmVkib7cXRkIkcwM2bSUrdvCF5UTkciGuhGdem5i7OC7wm9NxA9sz3G8kWz8+uxuszpiNEza0Y2jNQT1F39R7NPKf6cNvbkJfu8EYEffrS8X77WaLjjiYcKfRPwIjH8boJUxvgnuuZQ6WBsq2FH3USgqMgs65lJcnQiZ8DSvkn4B2e9sYepCqctvt6g6LX3Pp0dKB5MjjGeyjAjEeWVC7wImNww0rmdibx4fu98/iIxdpRQvUR+KqmyazEpSmltVOirpJqwyB88/Kccf/GX17y2pu6bVsaLuikVmukP+wzqUrLbl2u3e/A4s2dM/gy8eP7k6EQ2hfHC2TPyIlr8g7y2eAqwwPX9F+qz56pplo9cB3RPPFTom8il5y7mi//+JKmMqbnI2/7jeuAXvEoJ0s+puJTswc6Lr3+IeDRCd1dbnqDaYZyDyfFsxI0BfrRjXzaFsT1gWgj7/cejTuZNAxw4NkZHIppN1WDHwnvx9MPjI3k+8ko/s0Izikuh4Y4nFsUdwkpDSMQidCSi09qSDxq73Qjswc6Mcfzl+4dGuOyG7ex8fii7zx4QFuAD5y7mH9/2Ks48uStvwHRuZyJQ1se8pGZFwjhHLJed5yNXlHqiFv0JSCWiPB3C8fzROLGIZMXWn2JhTsfkaEgGeO3LZ3HpuYu56oLT8yJ67HkI9mv5sROL+Wcd2xSa59Cq11QJByr0IaSacQM/hWbNeq6GRlnxU3kd2//88fNP559+souMyU2xADA0PJnpMiK563aKB2/ANGjUkjeA68/PYxONCBHg5XNm1FTkW+WpqhCt3LcTARV6pSzTLRzPjsb5t+17Cw462hOzEr6bAExmypyqQEUjUjAPvB+R3BnBilJPVOiVkvjdIbVyNXS0xbKFTupJsUHHFafO4UcfeyPb9xxi5dJ5Krg+1AIPFyr0SpZiP+56hOMt65mZV++20XgTr6YbKsLKVFGhVwKh4XiKMn3R8EpFmSJrN22rav7Asp6ZOXnqFaXelBV6EfmuiBwUkSesbRtEpF9EHnP/Vlv7Pisiu0XkGRGpXfiHoiiKUhFBLPrvAxcV2P4tY8zZ7t89ACKyDLgEeI17zP8RkfLVohVFUZS6EaRm7FYRWRLwfGuAW40xY8BzIrIbOAfQDElK6Ak6SJocdbJy2nn4G40O6J5YVOMo/ISIfAjYAfy9MWYIWAhst9rsd7flISLrgHUAixdP3+n/Si7pjKH/yEigGbTGTVHQyrNta4Htz7cnn/UNJPMKk6gAK/Wg0sHY64BXAGcDA8A3pnoCY8xmY0yvMaa3u7u7wm4orYSXUGz/kJOsy84tU6htBqfwij8PTS3Zsn5VS4mnP6+N5rlRGkFFFr0x5gVvWUS+A9ztrvYDi6ymp7jblBMAf0KxYuX8/G39eWjChn2jKZRHJ6zvW2kdKrLoRaTHWn034EXk3AVcIiJtInIacAbwSHVdVKYLp3d3Zpf9RUJKtfXnoWl1kqMTWffUVPEmnyViEc48uUtFXmkIZS16EfkhcD4wX0T2A18AzheRs3Eys+4F1gMYY54UkduAPiAFXGWMCV5tWZnW2Em9ylmqXe1xIkAsQAHzRmH70r2aqf54+SD1Vm0KuY262uO0xSItnS9ICRdBom4+UGDzjSXafwX4SjWdUqYvXlKvIMItAm2xSEuIvJ9i+edrleCtlXL4K+FHp+cpiksQ4a1HgjdFqTcq9IoyBbTeqjIdUaFXlCmiCd6U6YYmNVOahj1hSlGU+qFC32SW9cwMVHg6bDRqwpSiKOq6aTonauTFiTRhSlGajVr0SlOYzhOmFGW6oRa9EohaP3m04oSpoFR7LU7UpzileahFrzSFLetX0dkea9kJU4oSJlToFUVRQo4KvaIoSshRoVcURQk5KvSKoighR4VeURQl5Gh4pdI0TsQZwYrSDFTolarQmHBFaX3Kum5E5LsiclBEnrC2zRWR+0XkWff/HHe7iMi/iMhuEXlcRN5Qz84riqIo5Qnio/8+cJFv29XAA8aYM4AH3HWAt+PUiT0DWAdcV5tuKoqiKJVSVuiNMVuBw77Na4Cb3OWbgHdZ2282DtuB2b5C4kqISY5OMJbKkBydaHZXFEWxqNRHv8AYM+AuHwAWuMsLgX1Wu/3utgGU0LJ207acotl9A0l2Pj+kqQ0UpUWoOrzSGGMAU7ahDxFZJyI7RGTH4OBgtd1QmoxdNBtg+55DTeqJoih+KhX6FzyXjPv/oLu9H1hktTvF3ZaHMWazMabXGNPb3d1dYTeUVmDL+lVsvGQ5EXHWI4KmHVaUFqJSob8LuNxdvhz4qbX9Q270zUrgqOXiUUKMVzQ7EYtw5sld6rZRlBairI9eRH4InA/MF5H9wBeArwK3iciVwPPA+93m9wCrgd3AMHBFHfqstChd7XGWL5qtsfWK0mKUFXpjzAeK7HpLgbYGuKraTimKoii1Q3PdKIqihBwVekVRlJCjQq8oihJyVOgVRVFCjgq9oihKyFGhVxRFCTkq9IqiKCFHhV5RFCXkqNAriqKEHBV6RVGUkKM1Y5WmoTlxFKUxqEWvKIoSclToFUVRQo4KvaIoSshRoVcURQk5KvSKoighR4VeURQl5FQVXikie4EkkAZSxpheEZkLbAGWAHuB9xtjhqrrpqIoilIptbDoLzDGnG2M6XXXrwYeMMacATzgriuKoihNoh6umzXATe7yTcC76vAaiqIoSkCqFXoD/EJEdorIOnfbAmPMgLt8AFhQ6EARWSciO0Rkx+DgYJXdUBRFUYpRbQqE84wx/SLyMuB+EXna3mmMMSJiCh1ojNkMbAbo7e0t2EZRFEWpnqosemNMv/v/IHAncA7wgoj0ALj/D1bbSUVRFKVyKhZ6EekUkS5vGXgr8ARwF3C52+xy4KfVdlJRFEWpnGpcNwuAO0XEO88PjDH3ishvgdtE5ErgeeD91XdTURRFqZSKhd4Yswd4fYHth4C3VNMpRVEUpXbozFhFUZSQo0KvKIoSclToFUVRQo4KvaIoSshRoVcURQk5KvSKoighR4VeURQl5KjQK4qihBwVekVRlJCjQq8oihJyVOgVRVFCjgq9oihKyFGhVxRFCTkq9IqiKCFHhV5RFCXkqNAriqKEnLoJvYhcJCLPiMhuEbm6Xq+jKIqilKYuQi8iUeDbwNuBZcAHRGRZPV5LURRFKU29LPpzgN3GmD3GmHHgVmBNnV5LURRFKUG9hH4hsM9a3+9uUxRFURpMxcXBq0VE1gHr3NUxEXmiWX2pgPnAi83uREAa3tfbPlbxoTl9reI8jUK/B/VhOvUVmtvfU4M0qpfQ9wOLrPVT3G1ZjDGbgc0AIrLDGNNbp77UnOnUX+1r/ZhO/dW+1o/p0N96uW5+C5whIqeJSAK4BLirTq+lKIqilKAuFr0xJiUinwDuA6LAd40xT9bjtRRFUZTS1M1Hb4y5B7gnYPPN9epHnZhO/dW+1o/p1F/ta/1o+f6KMabZfVAURVHqiKZAUBRFCTkNFfpyaRFEpE1Etrj7HxaRJY3sn9WPRSLyoIj0iciTIvLJAm3OF5GjIvKY+/f5ZvTV6s9eEdnl9mVHgf0iIv/iXtvHReQNTernq6xr9piIHBORT/naNPXaish3ReSgHfIrInNF5H4Redb9P6fIsZe7bZ4Vkcub1Nevi8jT7ud8p4jMLnJsye9Mg/q6QUT6rc96dZFjG55SpUh/t1h93SsijxU5tqHXtizGmIb84QzK/hFYCiSA3wPLfG0+DlzvLl8CbGlU/3z96AHe4C53AX8o0Nfzgbub0b8ifd4LzC+xfzXwc0CAlcDDLdDnKHAAOLWVri3wZuANwBPWtmuBq93lq4GvFThuLrDH/T/HXZ7ThL6+FYi5y18r1Ncg35kG9XUD8A8BvicltaNR/fXt/wbw+Va4tuX+GmnRB0mLsAa4yV2+HXiLiEgD+wiAMWbAGPOou5wEnmL6z+xdA9xsHLYDs0Wkp8l9egvwR2PM803uRw7GmK3AYd9m+7t5E/CuAoe+DbjfGHPYGDME3A9cVLeOUrivxphfGGNS7up2nHksTafIdQ1CU1KqlOqvq0vvB35Y737UgkYKfZC0CNk27hf1KDCvIb0rgus+Wg48XGD3KhH5vYj8XERe09CO5WOAX4jITnfWsZ9WTEtxCcV/KK10bQEWGGMG3OUDwIICbVrxGv8NzpNcIcp9ZxrFJ1w303eLuMRa8br+OfCCMebZIvtb5doCOhhbEhE5CbgD+JQx5phv96M4LofXA/8L+Emj++fjPGPMG3Ayhl4lIm9ucn9K4k6keyfwowK7W+3a5mCcZ/OWD1cTkc8BKeCWIk1a4TtzHfAK4GxgAMcdMh34AKWt+Va4tlkaKfRl0yLYbUQkBswCDjWkdz5EJI4j8rcYY37s32+MOWaMecldvgeIi8j8BnfT7k+/+/8gcCfO465NkOvfSN4OPGqMecG/o9WurcsLnqvL/X+wQJuWucYi8mHgL4HL3BtTHgG+M3XHGPOCMSZtjMkA3ynSh5a5rpDVpvcAW4q1aYVra9NIoQ+SFuEuwItUeB/wy2Jf0nri+t9uBJ4yxnyzSJuTvfEDETkH51o266bUKSJd3jLOYJw/SdxdwIfc6JuVwFHLFdEMilpErXRtLezv5uXATwu0uQ94q4jMcV0Qb3W3NRQRuQj4NPBOY8xwkTZBvjN1xzdO9O4ifWi1lCr/CXjaGLO/0M5WubY5NHLkFyfy4w84I+ifc7d9CecLCdCO8yi/G3gEWNqMEWrgPJxH88eBx9y/1cDHgI+5bT4BPIkTAbAdeGMz+ur2Zanbj9+7ffKurd1fwSkG80dgF9DbxP524gj3LGtby1xbnBvQADCB4w++Emes6AHgWeA/gLlu217gBuvYv3G/v7uBK5rU1904Pm3vu+tFsr0cuKfUd6YJff1X9/v4OI549/j76q7naUcz+utu/773XbXaNvXalvvTmbGKoighRwdjFUVRQo4KvaIoSshRoVcURQk5KvSKoighR4VeURQl5KjQK4qihBwVekVRlJCjQq8oihJy/j88sh66e9XieQAAAABJRU5ErkJggg==\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",
    "plt.plot(vecRp,Vc_Model_I(vecRp,params))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "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,params)\n",
    "    xi = np.sum((vecvRp-model)**2./(vecerrvRp)**2.)\n",
    "    return xi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5842.855534734712\n",
      "403203.80345974874\n"
     ]
    }
   ],
   "source": [
    "xi= np.sum((vecvRp-Vc_Model_I(vecRp,params))**2./(vecerrvRp)**2.)\n",
    "print xi\n",
    "xi= np.sum((vecvRp-Vc_Model_I(vecRp,params))**2.)\n",
    "print xi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "403203.80345974874\n"
     ]
    }
   ],
   "source": [
    "print xi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "LS = minimize(Xi2,params, method='BFGS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: 5842.855534734712\n",
      " hess_inv: array([[1, 0, 0, 0, 0, 0, 0, 0],\n",
      "       [0, 1, 0, 0, 0, 0, 0, 0],\n",
      "       [0, 0, 1, 0, 0, 0, 0, 0],\n",
      "       [0, 0, 0, 1, 0, 0, 0, 0],\n",
      "       [0, 0, 0, 0, 1, 0, 0, 0],\n",
      "       [0, 0, 0, 0, 0, 1, 0, 0],\n",
      "       [0, 0, 0, 0, 0, 0, 1, 0],\n",
      "       [0, 0, 0, 0, 0, 0, 0, 1]])\n",
      "      jac: array([0., 0., 0., 0., 0., 0., 0., 0.])\n",
      "  message: 'Optimization terminated successfully.'\n",
      "     nfev: 10\n",
      "      nit: 0\n",
      "     njev: 1\n",
      "   status: 0\n",
      "  success: True\n",
      "        x: array([ 4.27000e+08,  7.11000e+08, -2.00000e+00,  3.83000e+00,\n",
      "        2.96000e+00,  2.00000e+00,  1.42875e+09,  4.76250e+08])\n"
     ]
    }
   ],
   "source": [
    "print LS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "L2 = minimize(Xi2,params, method='TNC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: 5842.855534734712\n",
      "     jac: array([0., 0., 0., 0., 0., 0., 0., 0.])\n",
      " message: 'Local minimum reached (|pg| ~= 0)'\n",
      "    nfev: 1\n",
      "     nit: 0\n",
      "  status: 0\n",
      " success: True\n",
      "       x: array([ 4.27000e+08,  7.11000e+08, -2.00000e+00,  3.83000e+00,\n",
      "        2.96000e+00,  2.00000e+00,  1.42875e+09,  4.76250e+08])\n"
     ]
    }
   ],
   "source": [
    "print L2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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",
      "402000000.0 134000000.0\n"
     ]
    }
   ],
   "source": [
    "print params"
   ]
  },
  {
   "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
}