ISM_v6.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": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{equation}\n",
    "v_{c}^{2} = \\frac{G\\pi\\Sigma_{g}R^{2}}{R_{d}}\\int_{0}^{\\pi/2}(\\cos2u + 1)du\\int_{0}^{\\infty} d\\beta exp(-\\frac{R\\cos u \\cosh\\beta}{2R_{d}} -\\frac{R_{m}}{R\\cos u \\cosh\\beta})(\\frac{2R_{m}R_{d}}{R^{2}\\cos^{2}u\\cosh^{2}\\beta} -1)\n",
    "\\end{equation}$$\n",
    "$$\\begin{equation}\n",
    "\\frac{v_{c}^{2}R_{d}}{G\\pi\\Sigma_{g}} = R^{2}\\int_{0}^{\\pi/2}(\\cos2u + 1)du\\int_{0}^{\\infty} d\\beta exp(-\\frac{R\\cos u \\cosh\\beta}{2R_{d}} -\\frac{R_{m}}{R\\cos u \\cosh\\beta})(\\frac{2R_{m}R_{d}}{R^{2}\\cos^{2}u\\cosh^{2}\\beta} -1)\n",
    "\\end{equation}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{equation}\n",
    "v_{c}^{2} = \\frac{4\\pi G\\Sigma_{g}R^{2}}{2R_{d}}\\int_{0}^{1}ds\\frac{s^{2}}{\\sqrt{1-s^{2}}}\\int_{1}^{\\infty}dh \\frac{e^{-\\frac{Rsh}{2R_{d}}-\\frac{R_{m}}{Rsh}}}{\\sqrt{h^{2}-1}}\\left( \\frac{2R_{d}R_{m}}{R^{2}s^{2}h^{2}} - 1 \\right)\n",
    "\\end{equation}$$"
   ]
  },
  {
   "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 model I\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",
    "R_d_II = 3.2#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",
    "#parameters model II\n",
    "rho_b_II = 0.3e9#sM/kpc^3\n",
    "a_h_II = 1.9#kpc\n",
    "rho_h_II = 0.266e9#SM/kpc^3\n",
    "alpha_h_II = 1.63\n",
    "beta_h_II = 2.17\n",
    "R_d_II = 3.2#kpc\n",
    "Sigma_d_SD_II = (536.0e6)*0.75#SM/kpc^2\n",
    "Sigma_g_ISM_II = (536.0e6)*0.25#SM/kpc^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = np.linspace(0.1,15.,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vc2_ISM(R,R_d,Sigma_g_ISM,R_m=R_m):\n",
    "    def Int_ISM(h,s,R):\n",
    "        t1 = 1./np.sqrt(1.-s**2.)\n",
    "        t2 = 1./np.sqrt(h**2. - 1.)\n",
    "        t3 = np.exp(-(R*s*h)/(2.*R_d) - (R_m)/(R*s*h))\n",
    "        t4 = (2.*R_m*R_d)/(R**2.*h**2.) - s**2.\n",
    "        res = R**2.*t1*t2*t3*t4\n",
    "        return res\n",
    "    Integral_ISM = []\n",
    "    for i in range(0,len(R)):\n",
    "        integrando = lambda h,s: Int_ISM(h,s,R[i])\n",
    "        result = integrate.dblquad(integrando,0.1,1.,lambda s: 1., lambda s: np.inf)[0]\n",
    "        Integral_ISM.append(result)\n",
    "    return -4.*np.pi*G*Sigma_g_ISM*np.array(Integral_ISM)/(2.*R_d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vc2_ISM_I(h,s,R,R_d=R_d_I,R_m=R_m):\n",
    "    t1 = 1./np.sqrt(1.-s**2.)\n",
    "    t2 = 1./np.sqrt(h**2. - 1.)\n",
    "    t3 = np.exp(-(R*s*h)/(2.*R_d) - (R_m)/(R*s*h))\n",
    "    t4 = (2.*R_m*R_d)/(R**2.*h**2.) - s**2.\n",
    "    res = R**2.*t1*t2*t3*t4\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The integral is probably divergent, or slowly convergent.\n",
      "  warnings.warn(msg, IntegrationWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0022032922812955306, 0.014082320851371932, 0.03759504295288965, 0.07607010383984217, 0.1322864178826217, 0.2045662548687354, 0.2871030572873858, 0.3725246893188333, 0.4539091357015693, 0.5257231390165441, 0.5840311653359817, 0.6263506403320434, 0.6513831532663833, 0.6587340993593556, 0.6486662415630994, 0.6218998553105853, 0.5794579135771909, 0.5225496892703804, 0.45248517242422914, 0.3706133594436808, 0.2782786433294339, 0.17679073222902164, 0.06740457134346749, -0.04869240645733982, -0.17038770462218866, -0.29664735909821843, -0.4265168140065935, -0.5591198954342682, -0.6936564886405715, -0.8293993565092674, -0.9656904271748801, -1.1019367747073936, -1.2376064832220892, -1.3722244792592369, -1.505368449369367, -1.6366648807158253, -1.7657852678881782, -1.892442507975462, -2.0163874974109692, -2.137405937047053, -2.255315336829234, -2.369962232219721, -2.4812195884166663, -2.588984396036149, -2.6931754431149404, -2.79373125811154, -2.890608207427051, -2.983778744270783, -3.073229795371481, -3.158961276465404, -3.2409847282480633, -3.319322063776937, -3.3940044193053116, -3.4650710996645646, -3.532568615702786, -3.596549798108589, -3.6570729930059764, -3.714201320817645, -3.768002005826919, -3.8185457718580333, -3.865906228886677, -3.9101594728920324, -3.9513835138669924, -3.9896579154682463, -4.02506339691225, -4.05768149451603, -4.0875942490667505, -4.114883926979926, -4.139632767822134, -4.1619227594375285, -4.181835435629327, -4.199451696222797, -4.214851647114813, -4.2281144583759165, -4.239318240496439, -4.248539933078399, -4.255855211591084, -4.261338404126418, -4.265062420819931, -4.2670986947224145, -4.267517131132376, -4.266386067587938, -4.263772240260704, -4.259740759140378, -4.25435508876004, -4.247677037408724, -4.239766747247288, -4.230682694833347, -4.220481692354501, -4.209218894381947, -4.196947807921514, -4.183720305775252, -4.169586642674786, -4.154595474073426, -4.138793877128035, -4.122227373639397, -4.104939954783383, -4.086974107263726, -4.068370840784363, -4.049169716650281]\n"
     ]
    }
   ],
   "source": [
    "Integ = []\n",
    "for i in range (0,len(R)):\n",
    "    integrando = lambda h,s: vc2_ISM_I(h,s,R[i])\n",
    "    result = integrate.dblquad(integrando,0.1,1.,lambda s: 1., lambda s: np.inf)[0]\n",
    "    Integ.append(result)\n",
    "print Integ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x117601d10>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD8CAYAAACCRVh7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd8lfX9/vHXOztA2GGEsFREhgyJgNa2bpCq6LdWrQuVSmu1rf3a+nW0arW1VuvopMVRsa1S6sSJaN0VIcheEhkmYSQQAiEh+/3749z6i0ggIeM+Sa7n43EeOfmc+5xzHSBcuefH3B0REZH9iQk7gIiIRC+VhIiI1EolISIitVJJiIhIrVQSIiJSK5WEiIjUSiUhIiK1UkmIiEitVBIiIlKruLADNFT37t19wIABYccQEWlRFi1atN3dUw+2XIsviQEDBpCZmRl2DBGRFsXMNtVlOW1uEhGRWqkkRESkVioJERGplUpCRERqpZIQEZFaqSRERKRWKgkREalViz9PQkRar9KKKnaWlFNYUkFhSQV7yiopKa+kuKyKssoqyiurKa+spsodd3B3zIy4GCM21kiIjSExPpbEuBjaJcTSITGOlKQ4UpLi6ZwcT8fkeJLiY8P+mFFNJSEioXB38orKWJ9fzKYdxeQW7iV3515yC/eSX1RGflEZRWWVTZ6jXUIs3Tok0K19IqkpifTsmEjPlCR6dUqiT+dk0jon06tTUpstE5WEiDS53aUVrMjdxarNu1m7tYiPtxWxLm8PJeVVny8TG2P06phEWuckhqR15GsdIv9pd2mXQJd28XRKjqdDUhztE+NonxBHUnwMCXExxMfGEGuGGZgZ7k5VtVNZ7ZRXVVNaUUVZRTUl5VXsKaukuKySXXsrPr8VFJezY08ZO4rLyS4oIXNjATtLKr70GXp1TKJf13b07dqOgd3bMbB7BwZ2b89hqe1bdYGoJESkUVVXOx/nFZG5cSeLNu1kSXYhG7YXf/54akoig3umcH5GXw5Pbc/A7h3o360dvTslERfb8N2kZkZcrBEXC0nxsXRMiq/3a5RWVLFtdymbC0vZXLiXnJ17+bSghOyCEt5dl8/TH5XVeD9I75LMoB4pDOrZgcE9UziyZ+R+YlzLLw+VhIg02Mbtxby7Lp/3s3Ywf8MOCoPfxFNTEjmmX2fOG5PO8D6dGJ7WkW4dEkNOe3BJ8bH079ae/t3a7/fx4rJKNu4oZn1+MZ/k7yErL3J7d10+FVUORNaMjkjtwJDeKQxL68SwPh0ZltaJTsn1L60wqSREpN4qqqpZsKGA11dv4621+Z+vKfTpnMxpQ3oy/rBuHDugK327JmNmIadtfO0T4yL/8ad1+sJ4RVU1m3YUs2ZrEau37Gb1liLmry/guSWbP1+mf7d2jEjvzMj0Tozs25nhaZ1ITojeNQ5z97AzNEhGRobrKrAiTa+iqpr31m3n5eVbmLd6G4UlFSTExXD84d04aXAPvn5kKv27tWuVpdBQ2/eUsXLzblbk7mJZTiHLcnaxZVcpEFnjGNI7hdF9uzCmfxeO6delWcrVzBa5e8ZBl1NJiEht3J2lObt49qMcXli2hYLiclIS4zhlSA8mDu/N147sTrsEbZA4FHlFpSzN3sWS7Mh+myWfFlIc7MhPTUnk2AFdyOjflWMHdGVoWkdiYxq3NOpaEvrbFZEv2VlczjOLc5m9MJu124pIjIvh1CE9mTwqja8PTm0VO2TD1iMlidOGJnHa0J4AVFU7a7cWsWhTAYs27WThxp28vHwrACmJcYwZ0IVxA7sx/rCuHN2nU6Ps5K8LrUmIyOeW5+zisf9u5IVlmymvrGZk385ceGxfvjGi9yEdJSQNs7lwLws3FrBgQ+S2Lm8PAO0TYskY0JUHLxhFl/YJh/TaWpMQkTqprnbmrd7GjHfWs2jTTtolxHJ+RjoXj+vPkN4dw47XpqV1TmbyqD5MHtUHiOzbWLChgPnrd7By8+5mOVJKJSHSRpVXVvPs4hz++s561ucX07drMj8/cyjnjUlvcYdpthXdOyQy6ejeTDq6d7O9p0pCpI0pr6zmqUU5/OnNLHIL9zIsrSN/+PZozhjeq9m2c0vLUed/EWbW18zeNLNVZrbSzH4UjN9uZrlmtiS4TarxnJvMLMvM1prZhBrjE4OxLDO7scb4QDP7MBj/l5kd2sY2EfmSqmrnqUU5nPTbt7j52eWkpiTytyuO5cUfnMBZI9NUELJf9VmTqASud/ePzCwFWGRm84LHHnD339Zc2MyGAhcCw4A04HUzOzJ4+E/AaUAOsNDM5rj7KuA3wWvNMrO/AFOB6Yf64UQkchjrm2vz+M0ra1m7rYgR6Z341bnD+fqRqTqnQQ6qziXh7luALcH9IjNbDfQ5wFMmA7PcvQzYYGZZwNjgsSx3Xw9gZrOAycHrnQxcFCwzE7gdlYTIIVu7tYg7XlzJ+1k7GNCtHX+66BgmHd1L5SB1dkj7JMxsADAa+BD4CnCtmV0GZBJZ29hJpEDm13haDv+/VLL3GR8HdAMK3b1yP8uLSD0UlpRz/7yP+cf8TaQkxXP7WUO5eHx/4rVJSeqp3iVhZh2Ap4Hr3H23mU0H7gQ8+HofcGWjpvxyhmnANIB+/fo15VuJtCjuzrOLc/nVS6sp3FvBpeP7c92pRx7ysfQi9SoJM4snUhD/dPdnANx9W43HHwJeDL7NBfrWeHp6MEYt4zuAzmYWF6xN1Fz+C9x9BjADIifT1ecziLRWG7YXc/Mzy/lg/Q5G9e3M4+cO/9IF6ETqq84lYZGNmI8Aq939/hrjvYP9FQDnAiuC+3OAJ8zsfiI7rgcBCwADBpnZQCIlcCFwkbu7mb0JnAfMAqYAzzfkw4m0BVXVziPvree+1z4mMS6GX507nG8f24+YRr7Wj7RN9VmT+ApwKbDczJYEYzcD3zazUUQ2N20Evgvg7ivNbDawisiRUde4exWAmV0LzAVigUfdfWXwev8HzDKzXwKLiZSSiNQiK28PP/n3UpZkF3La0J786pzh9OiYFHYsaUV07SaRFsjd+fv8Tdz18mqS42O5/exhnD0yTUctSZ3p2k0irVReUSk3PLWMt9bm8/UjU7n3vBFae5Amo5IQaUHeW7ed6/61mKLSSu6YPIxLx/fX2oM0KZWESAtQWVXN795Yxx/fzOKI1A48cdV4juyZEnYsaQNUEiJRbvueMn7wxGI+WL+D8zPS+cXZw6N6TmRpXVQSIlFsSXYhV/9jEQXF5dz3rZF8c0x62JGkjVFJiESp2Quz+dlzK+jRMZGnrz6e4X10Ypw0P5WESJSpqnbuenk1j7y3ga8O6s7vLxyty2pIaFQSIlGkqLSCHz65mDfX5nP58QP42TeGaJ4HCZVKQiRKbC7cyxV/W8gn+Xv41bnDuXhc/7AjiagkRKLBqs27ueKxBZSUVfHYFWM5YVD3sCOJACoJkdC9t2473/vHIjokxvHvq4/jqF4dw44k8jmVhEiIXly2mR//awmHde/AY1ceS+9OyWFHEvkClYRISJ748FNueW45Y/p14ZHLj6VTcnzYkUS+RCUhEoLpb33Cb15dw4mDU5l+8RidQS1RSyUh0ozcnQdeX8fv31jHWSPTuO9bI0mI0yGuEr1UEiLNxN25Z+5apr/1Cd8ak87d3xxBrGaPkyinkhBpBu7Or15azcPvbeDicf24c/JwTS8qLYJKQqSJuTu/fmUND7+3gcuPH8BtZw3VHBDSYmhjqEgTcnfunbuWGe+sZ8px/VUQ0uKoJESa0IOvr+PPb33CReP6cfvZw1QQ0uKoJESayIx3PuF3b6zjW2PS+eXk4SoIaZFUEiJN4MkFn3LXy2s4c0Rv7v7mCO2klharziVhZn3N7E0zW2VmK83sR8F4VzObZ2brgq9dgnEzs9+bWZaZLTOzY2q81pRg+XVmNqXG+BgzWx485/emX72kBXpx2WZufnY5Jw5O5f7zR+kwV2nR6rMmUQlc7+5DgfHANWY2FLgReMPdBwFvBN8DnAEMCm7TgOkQKRXgNmAcMBa47bNiCZa5qsbzJh76RxNpfu+uy+fH/1pCRv8uTL94jE6Ukxavzv+C3X2Lu38U3C8CVgN9gMnAzGCxmcA5wf3JwOMeMR/obGa9gQnAPHcvcPedwDxgYvBYR3ef7+4OPF7jtUSi3orcXXzv74s4PLUDD085VpfakFbhkH7NMbMBwGjgQ6Cnu28JHtoK9Azu9wGyazwtJxg70HjOfsZFol52QQmX/20hndsl8NgVY3WxPmk16l0SZtYBeBq4zt1313wsWAPwRsp2oAzTzCzTzDLz8/Ob+u1EDmhncTlTHl1ARVU1M688ll6dksKOJNJo6lUSZhZPpCD+6e7PBMPbgk1FBF/zgvFcoG+Np6cHYwcaT9/P+Je4+wx3z3D3jNTU1Pp8BJFGVVZZxbS/Z5JTuJeHp2RwRI+UsCOJNKr6HN1kwCPAane/v8ZDc4DPjlCaAjxfY/yy4Cin8cCuYLPUXOB0M+sS7LA+HZgbPLbbzMYH73VZjdcSiTruzg1PLWPhxp3c962RHDuga9iRRBpdfa7d9BXgUmC5mS0Jxm4G7gZmm9lUYBNwfvDYy8AkIAsoAa4AcPcCM7sTWBgsd4e7FwT3vw88BiQDrwQ3kaj04OvreH7JZn46YTBnjUwLO45Ik7DIboSWKyMjwzMzM8OOIW3M80ty+dGsJZw3Jp17zxuhs6mlxTGzRe6ecbDldBC3SD0tzS7kp08tY+zArtx17tEqCGnVVBIi9bBtdylXPZ5Jj5RE/nKJTpaT1k/zSYjUUWlFFdMez2RPWSXPTD2eru0Two4k0uRUEiJ14O7c8uwKlubs4q+XjuGoXh3DjiTSLLSuLFIHf5+/iac/yuG6UwcxYVivsOOINBuVhMhBZG4s4I4XVnHKUT344cmDwo4j0qxUEiIHkLe7lKv/+RHpXZK5/4JRmhdC2hztkxCpRWVVNdc+sZg9pZX8Y+o4XbRP2iSVhEgt7n1tLQs2FvDgBaMY3EvXZJK2SZubRPZj3qpt/PXt9Vw0rh/njNYV66XtUkmI7CO7oITrZy9hWFpHbj1zaNhxREKlkhCpobyymmufXIwDf774GJLiNbuctG3aJyFSw32vrWVpdiF/vvgY+ndrH3YckdBpTUIk8ObaPP76znouHtePSUf3DjuOSFRQSYgQuXDf9bOXclSvFH6u/RAin1NJSJtXXe387+wl7C2v4o8XjdZ+CJEaVBLS5j3y3gbez9rBrWcN1RzVIvtQSUibtnLzLu6Zu4YJw3py4bF9w44jEnVUEtJm7S2v4kezltC1fQJ3/4+mIBXZHx0CK23Wr19ZTVbeHv4xdRxdNIGQyH5pTULapLc/zufxDzYx9YSBnDCoe9hxRKKWSkLanF0lFdzw1FIG9ejATycMDjuOSFSrc0mY2aNmlmdmK2qM3W5muWa2JLhNqvHYTWaWZWZrzWxCjfGJwViWmd1YY3ygmX0YjP/LzLT+L03i58+vYMeecu4/f5QOdxU5iPqsSTwGTNzP+APuPiq4vQxgZkOBC4FhwXP+bGaxZhYL/Ak4AxgKfDtYFuA3wWsdAewEph7KBxI5kBeXbWbO0s388JRBHJ3eKew4IlGvziXh7u8ABXVcfDIwy93L3H0DkAWMDW5Z7r7e3cuBWcBkixxWcjLwVPD8mcA5dc0mUhd5RaX87LkVjEzvxPdPPDzsOCItQmPsk7jWzJYFm6O6BGN9gOway+QEY7WNdwMK3b1yn3GRRuHu/OzZFZSUV3Hf+SOJi9XuOJG6aOhPynTgcGAUsAW4r8GJ6sDMpplZppll5ufnN8dbSgs3Z+lmXlu1jetPO1JnVYvUQ4NKwt23uXuVu1cDDxHZnASQC9Q8fTU9GKttfAfQ2czi9hmv7X1nuHuGu2ekpqY25CNIG5BXVMptc1Yyul9nvvPVw8KOI9KiNKgkzKzm9ZTPBT478mkOcKGZJZrZQGAQsABYCAwKjmRKILJze467O/AmcF7w/CnA8w3JJgKRzUy3BJuZ7j1vJLExOqtapD7qfMa1mT0JnAh0N7Mc4DbgRDMbBTiwEfgugLuvNLPZwCqgErjG3auC17kWmAvEAo+6+8rgLf4PmGVmvwQWA480+NNJm/fisi3MW7WNm844iiN6dAg7jkiLY5Ff4luujIwMz8zMDDuGRKGC4nJOu/9t0rsk8/TVx2tntUgNZrbI3TMOtpyu3SSt1p0vrmLX3gr+edU4FYTIIdJPjrRKb67J49nFuXz/pCM4qlfHsOOItFgqCWl19pRVcsuzyzmyZweuOUknzYk0hDY3Savz27lr2bK7lKcvPp7EOF2bSaQhtCYhrcqS7EJmfrCRy8b355h+XQ66vIgcmEpCWo2KqmpufHoZPVOS+IkuAS7SKLS5SVqNh9/dwJqtRcy4dAwpSfFhxxFpFbQmIa3CpztKePD1j5kwrCenD+sVdhyRVkMlIS2eu/Pz51cQF2P84uzhYccRaVVUEtLivbJiK29/nM/1pw+mV6eksOOItCoqCWnRikor+MULKxmW1pHLjusfdhyRVkc7rqVFu++1j8krKmPGpRm69IZIE9BPlbRYK3J38fgHG7lkXH9G9u0cdhyRVkklIS1SdbXzs+dW0LV9os6JEGlCKglpkWZnZrMku5CbJx1Fp2SdEyHSVFQS0uLsLC7nN6+uYezArpw7uk/YcURaNZWEtDj3zF3D7tJK7pw8HDNNRyrSlFQS0qIsyS5k1sJsrjh+AIN7pYQdR6TVU0lIi1Fd7dz6/ApSOyTyo1MHhR1HpE1QSUiLMTszm2U5u7h50hBdwE+kmagkpEXYVVLBPXPXcuyALkwelRZ2HJE2QyUhLcL989ZSWFLOL87WzmqR5lTnkjCzR80sz8xW1BjrambzzGxd8LVLMG5m9nszyzKzZWZ2TI3nTAmWX2dmU2qMjzGz5cFzfm/6n0ACqzbv5u/zN3HJ+P4MTesYdhyRNqU+axKPARP3GbsReMPdBwFvBN8DnAEMCm7TgOkQKRXgNmAcMBa47bNiCZa5qsbz9n0vaYPcndtfWEmn5Hj+97Qjw44j0ubUuSTc/R2gYJ/hycDM4P5M4Jwa4497xHygs5n1BiYA89y9wN13AvOAicFjHd19vrs78HiN15I27KXlW1iwoYCfTBhM53YJYccRaXMauk+ip7tvCe5vBXoG9/sA2TWWywnGDjSes59xacP2lldx10urGdq7Ixce2y/sOCJtUqPtuA7WALyxXu9AzGyamWWaWWZ+fn5zvKWEYPrbn7B5Vym3nz2M2BjtohIJQ0NLYluwqYjga14wngv0rbFcejB2oPH0/Yzvl7vPcPcMd89ITU1t4EeQaJRdUMJf3/6Es0amMXZg17DjiLRZDS2JOcBnRyhNAZ6vMX5ZcJTTeGBXsFlqLnC6mXUJdlifDswNHtttZuODo5ouq/Fa0gbd/coazOCmM44KO4pIm1bnmenM7EngRKC7meUQOUrpbmC2mU0FNgHnB4u/DEwCsoAS4AoAdy8wszuBhcFyd7j7ZzvDv0/kCKpk4JXgJm3Qh+t38NLyLVx36iDSOieHHUekTbPIroSWKyMjwzMzM8OOIY2kqto56w/vUVhSzhvXn0hyQmzYkURaJTNb5O4ZB1tOZ1xLVPl3ZjartuzmxklDVBAiUUAlIVFjd2kFv31tLRn9u3DWiN5hxxERVBISRf70nyx2FJdz21nDdH0mkSihkpCosGlHMX97fyPfPCado9M7hR1HRAIqCYkKv355DXGxxk8nDA47iojUoJKQ0H3wyQ5eXbmV7594OD07JoUdR0RqUElIqKqqnTtfXEWfzsl856uHhR1HRPahkpBQPb0oJ3LI6xlHkRSvQ15Foo1KQkKzp6ySe19byzH9OnOmDnkViUoqCQnNX976hPyiMn5+5lAd8ioSpVQSEorcwr089O56Jo9KY3S/Lgd/goiEQiUhobjn1TUA3DBRV3kViWYqCWl2iz/dyfNLNvOdrw6kj67yKhLVVBLSrNydX760mu4dErn6xCPCjiMiB6GSkGb18vKtLNq0k5+cfiQdEus8nYmIhEQlIc2mtKKKu19dzVG9UvhWRt+DP0FEQqeSkGYz878byS7Yy8++MZTYGB3yKtISqCSkWezYU8Yf/5PFyUf14IRB3cOOIyJ1pJKQZvHg6+soqaji5kk65FWkJVFJSJPLyiviiQWfctHYfhzRIyXsOCJSDyoJaXJ3vbyGdgmxXHfqoLCjiEg9qSSkSb27Lp//rMnjBycfQbcOiWHHEZF6UklIk6mqdn710mr6dk1myvEDwo4jIoegUUrCzDaa2XIzW2JmmcFYVzObZ2brgq9dgnEzs9+bWZaZLTOzY2q8zpRg+XVmNqUxskl4/p2ZzZqtRdw4cQiJcZorQqQlasw1iZPcfZS7ZwTf3wi84e6DgDeC7wHOAAYFt2nAdIiUCnAbMA4YC9z2WbFIy7OnrJLfvvYxY/p3YdLRvcKOIyKHqCk3N00GZgb3ZwLn1Bh/3CPmA53NrDcwAZjn7gXuvhOYB0xswnzShKa/lcX2PZorQqSla6yScOA1M1tkZtOCsZ7uviW4vxXoGdzvA2TXeG5OMFbb+JeY2TQzyzSzzPz8/Eb6CNJYcnaW8NC7GzhnVBqj+nYOO46INEBjXWHtBHfPNbMewDwzW1PzQXd3M/NGei/cfQYwAyAjI6PRXlcaxz2vrsXQXBEirUGjrEm4e27wNQ94lsg+hW3BZiSCr3nB4rlAzau7pQdjtY1LC/LRpzuZs3Qz0752GGmaK0KkxWtwSZhZezNL+ew+cDqwApgDfHaE0hTg+eD+HOCy4Cin8cCuYLPUXOB0M+sS7LA+PRiTFsLdufPFVfRISeR7Xz887Dgi0ggaY3NTT+DZYOdkHPCEu79qZguB2WY2FdgEnB8s/zIwCcgCSoArANy9wMzuBBYGy93h7gWNkE+ayZylm1n8aSH3nDeC9porQqRVMPeWvUk/IyPDMzMzw47R5u0tr+KU+96iS/sEXrj2BGJ0KXCRqGZmi2qcslArnXEtjeKhd9ezeVcpt545VAUh0oqoJKTBtu0uZfpbn3DG8F6MO6xb2HFEpBGpJKTB7nl1LVXVzk1nDAk7iog0MpWENMiS7EKe/iiHK04YQL9u7cKOIyKNTCUhh6y62rl9zkpSUxL5wcmaK0KkNVJJyCF7fmkuS7ILuWHCYDrokFeRVkklIYekuKySu19Zw4j0TnzzmPSw44hIE1FJyCGZ/tYnbNtdxm1nDdMhryKtmEpC6u3THSXMeHc9k0elMaa/pvwQac1UElJvd7y4irgY0yGvIm2ASkLq5a21eby+ehs/OHkQvTolhR1HRJqYSkLqrLyymjteWMXA7u258oQBYccRkWagkpA6e/T9DazfXsytZw0lMS427Dgi0gxUElInW3eV8oc31nHqkB6cNLhH2HFEpJmoJKRO7nxpFZXVzq1nDgs7iog0I5WEHNQ7H+fz0rItXHPSEbo+k0gbo5KQAyqrrOK2OSsZ0K0d0752WNhxRKSZ6YI7ckAPvbOeDduLmXnlWJLitbNapK3RmoTU6tMdJfzxzSwmHd2Lrx+ZGnYcEQmBSkL2y9255bnlxMXE8PMzh4YdR0RCopKQ/ZqzdDPvrtvOTycMpnen5LDjiEhIVBLyJYUl5dzxwipG9e3MJeP7hx1HREIUdSVhZhPNbK2ZZZnZjWHnaYt+/fIaCvdW8Ov/OZpYXQZcpE2LqpIws1jgT8AZwFDg22amDeLN6L9Z2/lXZjbf+epAhvTuGHYcEQlZVJUEMBbIcvf17l4OzAImN8Ubrd1aRN7uUty9KV6+RSouq+SGp5cxsHt7fnzqkWHHEZEoEG3nSfQBsmt8nwOM23chM5sGTAPo16/fIb3RtU98xLq8PXRIjOPw1PaMSO/MqUN7Mv6wrm324nX3vLqG3MK9zP7ucTonQkSA6CuJOnH3GcAMgIyMjENaFfjF5GGs27aH9fl7yMrfw1OLcvj7/E10SIxj4vBeXHPSEQzs3r5Rc0ezD9fvYOYHm7j8+AEcO6Br2HFEJEpEW0nkAn1rfJ8ejDW64w/vzvGHd//8+9KKKt7P2s5rK7fx/NJcnl2cyzmj+vDDU46gf7fWXRZ7y6v4v6eX0a9rO26YODjsOCISRaJtn8RCYJCZDTSzBOBCYE5zvHFSfCynDOnJb84bwTs3nMQVxw/gxWWbOf2Bd/jnh5ta9b6Lu15ezcYdJdz9zaNplxBtvzeISJiiqiTcvRK4FpgLrAZmu/vK5s7RIyWJn505lHduOImxA7tyy7MruPaJxezaW9HcUZrcm2vy+Pv8TXznhIFfWLMSEQGwlv4bckZGhmdmZjbZ61dXO399Zz2/fW0tfTon84+p41rN5bJ37CljwoPv0r1DAs9d8xXtrBZpQ8xskbtnHGy5qFqTiEYxMcbVJx7O7O+OZ3dpBef95b98vK0o7FgN5u7c9Mxydu+t4IELRqkgRGS/VBJ1NKZ/V/417TgcuOCvH7A8Z1fYkRrkiQWf8tqqbfx0wmCdNCcitVJJ1MPgXin8+7vH0S4hjosems+K3JZZFCtyd/GLF1bxtSNTmXrCwLDjiEgUU0nU04Du7fn3946jY3I8Ux5dwPr8PWFHqpfdpRVc88RHdG2XwIMXjCJG12YSkQNQSRyCtM7J/H3qWAAufWQBmwv3hpyobtydG59eRs7OvfzxotF0bZ8QdiQRiXIqiUN0WGoHZl45lt17K7j0kQ8pKC4PO9JBPfzuBl5evpUbJgwmQ2dVi0gdqCQaYHifTjxy+bHk7NzL1JkL2VteFXakWr25Jo+7XlnNGcN7cdVXDws7joi0ECqJBho7sCu/u3A0S7MLufaJj6isqg470pes21bED55czNDeHbnv/JHaDyEidaaSaAQTh/fijsnDeWNNHrc8uyKqLuGxs7icqTMzSYqP5aHLMnTZDRGpF/2P0UguGd+fvN2l/P4/WXRpn8CNZxwVdiT2lFVy+WML2bq7lFnTxpPWWXNVi0j9qCQa0Y9PO5KCknL+8vYnpCTFcc1JR4SWpbSiiqtmZrIidxd/uWQMx/TrEloWEWm5VBKNyMy44+yK4NKoAAAH30lEQVTh7Cmt5N65a0lJiuOy4wY0e46KqmqufWIxH6zfwQMXjOS0oT2bPYOItA4qiUYWE2Pc+62RFJdXcevzK4kx45Lx/Zvt/Usrqrhu1hJeX72NOyYP49zR6c323iLS+mjHdROIj43hD98ezalDevCz51bw57eymuV9i0oruPxvC3h15VZuPXNoKGsxItK6qCSaSFJ8LNMvGcPkUWnc8+pa7n5lTZMe9ZRXVMqFM+aTuXEnD14wiit1TSYRaQTa3NSE4mNjeOD8UXRIjOMvb3/Cph3F/Oa8EXRMim/U9/ngkx38aNZiikoreXhKBicO7tGory8ibZfWJJpYTIzxy3OGc9MZR/Haqm2c/Yf3WLV5d6O8dlW18+DrH3Pxw/PpkBjH01cfr4IQkUalkmgGZsZ3v344T141nr0VVZz75/f505tZDbqMx4INBfzPn9/nwdfXcc6oPrzwgxMYmqZ5IUSkcWn60ma2fU8Ztzy7nLkrt9GrYxI/mTCYc0f3IbaOl8pYvWU3D77+8efPv2nSUUwe1aeJU4tIa1PX6UtVEiFZsKGAX720iqU5u+jWPoHTh/XijOG9GJHeiU7J8ZhFSmN3aQWf7ijh7Y/zmbNkM2u3FdE+IZarTzycqSccRnKCph0VkfpTSbQA1dXO66u38cKyLbyxehslweanhNgYUlMSKS6vpLCk4vPlM/p34exRaZw5Ik1zQYhIg9S1JHR0U4hiYozTh/Xi9GG9KK2o4v2s7WzYXkz+njLyd5eRnBBL367t6Ne1HSPSO5HepV3YkUWkjWlQSZjZ7cBVQH4wdLO7vxw8dhMwFagCfujuc4PxicDvgFjgYXe/OxgfCMwCugGLgEvdPfpn8mkkSfGxnDJEl88QkejSGEc3PeDuo4LbZwUxFLgQGAZMBP5sZrFmFgv8CTgDGAp8O1gW4DfBax0B7CRSMCIiEqKmOgR2MjDL3cvcfQOQBYwNblnuvj5YS5gFTLbIXtqTgaeC588EzmmibCIiUkeNURLXmtkyM3vUzD67HnUfILvGMjnBWG3j3YBCd6/cZ1xEREJ00JIws9fNbMV+bpOB6cDhwChgC3BfE+f9LNM0M8s0s8z8/PyDP0FERA7JQXdcu/updXkhM3sIeDH4NhfoW+Ph9GCMWsZ3AJ3NLC5Ym6i5/P4yzQBmQOQQ2LrkExGR+mvQ5iYz613j23OBFcH9OcCFZpYYHLU0CFgALAQGmdlAM0sgsnN7jkdO1ngTOC94/hTg+YZkExGRhmvoeRL3mNkowIGNwHcB3H2lmc0GVgGVwDXuXgVgZtcCc4kcAvuou68MXuv/gFlm9ktgMfBIA7OJiEgD6YxrEZE2qM1clsPM8oFN9Xxad2B7E8RpTMrYcNGeD5SxsShj/fV399SDLdTiS+JQmFlmXRo0TMrYcNGeD5SxsShj09F8EiIiUiuVhIiI1KqtlsSMsAPUgTI2XLTnA2VsLMrYRNrkPgkREambtromISIiddCmSsLMJprZWjPLMrMbw86zLzPra2ZvmtkqM1tpZj8KO1Ntgku/LzazFw++dPMzs85m9pSZrTGz1WZ2XNiZ9mVmPw7+nleY2ZNmlhQFmR41szwzW1FjrKuZzTOzdcHXLgd6jZAy3hv8XS8zs2fNrHO0Zazx2PVm5mbWPYxs9dVmSuIgc1lEi0rgencfCowHronCjJ/5EbA67BAH8DvgVXc/ChhJlGU1sz7AD4EMdx9O5AoEF4abCoDHiMwBU9ONwBvuPgh4I/g+TI/x5YzzgOHuPgL4GLipuUPt4zG+nBEz6wucDnza3IEOVZspCWqZyyLkTF/g7lvc/aPgfhGR/9ii7pLpZpYOfAN4OOws+2NmnYCvEVzaxd3L3b0w3FT7FQckm1kc0A7YHHIe3P0doGCf4clE5niBKJjrZX8Z3f21GlMNzCdykdDQ1PLnCPAAcAORSxm1CG2pJGqbyyIqmdkAYDTwYbhJ9utBIv/Qq8MOUouBRKbU/VuwSexhM2sfdqia3D0X+C2R3yi3ALvc/bVwU9Wqp7tvCe5vBaJ9nt0rgVfCDrGvYHqFXHdfGnaW+mhLJdFimFkH4GngOnffHXaemszsTCDP3ReFneUA4oBjgOnuPhooJvxNJF8QbNefTKTQ0oD2ZnZJuKkOLrhic9T+FmxmtxDZbPvPsLPUZGbtgJuBW8POUl9tqSQONMdF1DCzeCIF8U93fybsPPvxFeBsM9tIZJPdyWb2j3AjfUkOkOPun62FPUWkNKLJqcAGd8939wrgGeD4kDPVZttn0wIEX/NCzrNfZnY5cCZwsUffsf2HE/mFYGnws5MOfGRmvUJNVQdtqST2O5dFyJm+IJjr+xFgtbvfH3ae/XH3m9w93d0HEPkz/I+7R9VvwO6+Fcg2s8HB0ClELlsfTT4FxptZu+Dv/RSibOd6DXOIzPECUTrXi5lNJLIJ9Gx3Lwk7z77cfbm793D3AcHPTg5wTPBvNaq1mZIIdmp9NpfFamB2jbksosVXgEuJ/Ha+JLhNCjtUC/UD4J9mtozI9Lp3hZznC4K1nKeAj4DlRH4WQz8j18yeBD4ABptZjplNBe4GTjOzdUTWgO6Owox/BFKAecHPzV+iMGOLpDOuRUSkVm1mTUJEROpPJSEiIrVSSYiISK1UEiIiUiuVhIiI1EolISIitVJJiIhIrVQSIiJSq/8Hsf9F1XYdPU8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "integ_array = np.array(Integ)\n",
    "C_integ_I = -2.*np.pi*G*Sigma_g_ISM_I/R_d_I\n",
    "plt.plot(R,C_integ_I*integ_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in sqrt\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.0, 15.0)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8nXWd9vHPN/vSpmmadKFJm+4LZSmkpcIUERDLomXEQVAREKw6qIzDqKA+4jwuw2yijI+MDGCrQBFZK+DCUkAEulu609CkbULbJG3TZk9O8n3+OAcmdMtyktwn51zv1yuvc87v3Ce5KMmVX373fe7b3B0REYlfSUEHEBGR/qWiFxGJcyp6EZE4p6IXEYlzKnoRkTinohcRiXMqehGROKeiFxGJcyp6EZE4lxJ0AID8/HwvLi4OOoaIyKCyZs2aGncv6Gq7mCj64uJiVq9eHXQMEZFBxcx2dmc7Ld2IiMQ5Fb2ISJxT0YuIxDkVvYhInFPRi4jEORW9iEicU9GLiMS5mDiOXkQST2uog7rmNupbQtQ1h2hsbaexNURTazvNoXZaQx20hDpoDXUQ6nDaO5xQu9P+7uVPI7dmRnKSkWSQkpxESpKRmpxEWkoS6SlJpKckk5GaRGZqMplpyWSlpZCdnsyQ9BSy01NITY7/+a6KXkT6TH1LiD21Tew93My+wy3sO9xMTX0L++tb2d/QwoGGNg41tnKwsY2mtvag4wKQkZrE0IxUcjJSyMlMJTczldysNIZlppKXncbw7DTystLIy04jf0ga+UPSGZaZSlKSBR2921T0ItJtHR3OnsPN7Kiup3x/I7v2N7BzfyMVB5uorG3iUFPbUa8Zkp7CiCFpjMhOY2xuBieflMPwrFSGZaYyNCOVIekpDMlIITstJTLjTiYzNZm0lKT3PlKTkkhOsvc+OnN3OpzwjL+jg7Z2J9TeQWt7By1t4b8KmtvaaWxtp7mtnYbWEI0t7dS3hGhoCVHXEqKuuY3DTSEON7dRU99KaXU9tQ1t1LWEjvnvkJxk5A9JY+TQDEYOTWdkTgajctIZlZPB6JwMRg/LYMywDIZlpmIW/C8EFb2IHMXd2Xu4ma176ti6t45tew+zvaqeHdUN75uJp6ckMS4vi6K8LM4cP5yTcjM5KTdcdqNyMhiZk05WWv/WjJmRbOHyTevj3Y6toQ5qG1vZ39DKwYZWqiN/ndTUt1Bd10JVXQvvHGpmfUUtNfWtR70+MzWZMbkZjM3N5KRhmYwdnknh8EwKh2dRODyTUTkZR/3i6g8qehFh3+Fm1u2qZX1FLRsrD7HpncMcaPjf4hozLIOpo4Zy1oQRTBqZzcT8IUzIz2bk0PRBtYTRU2kpSYzMyWBkTkaX27aGOqipb2HPoWb2Hmpmz6Em3ql997aJLXsOH/XLIDXZ3iv9cXlZjB+Rxbi8bIrzsxifl01mWnKf/Hd0WfRmdj9wGVDl7rM6jX8FuAloB55x929Exm8DboiMf9Xd/9gnSUWkT7R3OFv2HGZ1+QFWlR9kzc6D7D3cDISLZ8rIoVw4YyQnnzSMGWNymDZqKMOyUgNOHfvSUpIif9FkHneb5rZ2KmubqDzYxO6Djew+8O5tI89s2ENt4/uXvkblpDN+RDYT87Mpzs9mQn74/vgR2aSldP+vl+7M6BcDPwN+9e6AmX0IWAic5u4tZjYyMj4TuAo4GTgJeN7Mprp7bOx1EUlA7s62fXW8Vrqf197ez4qy/dQ1h9eeTxqWwZwJecwuyuX0cbnMHJNDRmrfzCLlaBmpyUwqGMKkgiHHfP5QUxu79jey80AD5TUNlO9vpKymgec272N/p7+wkgzG5WV1++t2WfTu/oqZFR8x/CXgDndviWxTFRlfCDwcGS8zs1JgLvB6txOJSNQONbXx6vYaXtpWxctvVVNV1wLA+BFZXHbqGOZNHEFJcR5jTzD7lIE3LDOVUwqHcUrhsKOeO9TURnlNAztqwvtKdlQ38HI3P29v1+inAvPN7IdAM/BP7r4KGAu80Wm7isiYiPSz3QcaeX7LPp7bvI+VZQcIdTg5GSnMn1rAB6cUcPbkERQO7/4sUGLLsMxUTivK5bSi3PfGfv6Z7r22t0WfAuQB84A5wCNmNrEnn8DMFgGLAMaNG9fLGCKJbdf+8Nrusxv2sKHyEABTRg7h8+dO5MIZIzmtMJeUBHhDkJxYb4u+Anjc3R1YaWYdQD5QCRR12q4wMnYUd78HuAegpKTEe5lDJOHU1Lfw9Pp3eOKv77B+dy0ApxUO47aLp/ORk0dTnJ8dcEKJNb0t+ieBDwHLzWwqkAbUAMuAh8zsx4R3xk4BVvZFUJFE1hrq4MWt+3hkdQUvv1VNe4czY0wOt108nUtOGUNRD3bMSeLpzuGVS4HzgHwzqwBuB+4H7jezjUArcG1kdr/JzB4BNgMh4CYdcSPSe2U1DTy0YiePr61kf0Mro3LS+fz8ifzt7LFMGz006HgySJh78KsmJSUlrouDi4SF2jt4YWsVD7yxkz9vryElybhwxig+OaeI+VPyteYu7zGzNe5e0tV2emesSIyoa27jN6t2s/i1cioONjFmWAa3fHgqn5xT1K13Zoocj4peJGB7DzVz36s7WLpyN/UtIeYUD+c7l87gwhmjNHuXPqGiFwlIeU0Dv3jlbR5bU0m7O5eeMoYb50/g1MLcrl8s0gMqepEBVl7TwF0vbufJdZWkJCdx5ZxCvnDuJB05I/1GRS8yQHYfaOS/XtzOY2srSUkyPnfOBBZ9cCIjh2r9XfqXil6knx1saOVny0v59es7weCzHxjPlz44STtYZcCo6EX6SXNbO/f/pYy7l79NQ2uIvzuziH/48BTGDNOJxGRgqehF+pi78/uNe/nRs1uoONjEhTNG8o0F05k6Sm9wkmCo6EX60La9dXz3qY2sKDvA9NFDeejGszh7cn7QsSTBqehF+kBDS4ifvrCd+14tIycjhR9cPour5hTpOHiJCSp6kSj9adNebl+2iT2HmrlqThHfXDCd4dlpQccSeY+KXqSXquta+N7vNvHMm3uYPnooP/vUbM4cnxd0LJGjqOhFesjdeeqv7/C9322isaWdf7poKl/44CRStUwjMUpFL9ID++tb+PYTG/nDpr3MHpfLv11xKlN0NI3EOBW9SDc9v3kftz6+gcNNbdx28XRunD+R5CQLOpZIl1T0Il1obmvn+09v5sEVu5gxJocHbpzL9NE5QccS6TYVvcgJvLWvjq88tI5t++pYdO5EbrloKukpyUHHEukRFb3Icfxm1S5uX7aJIekpLPncXD44tSDoSCK90uVhAmZ2v5lVRa4Pe+Rzt5iZm1l+5LGZ2V1mVmpmb5rZGf0RWqQ/Nbe18/Xfruebj23gzPHDefbm+Sp5GdS6M6NfDPwM+FXnQTMrAi4CdnUavhiYEvk4C7g7cisyKOzc38CXHljL5j2H+er5k7n5wqna4SqDXpdF7+6vmFnxMZ66E/gG8FSnsYXArzx8xfE3zCzXzMa4+56+CCvSn155q5ovP7QWgPuvK+H86aMCTiTSN3q1Rm9mC4FKd19v9r7Zzlhgd6fHFZExFb3ELHfnvlfL+NGzW5g6aij3XFPCuBG62pPEjx4XvZllAd8ivGzTa2a2CFgEMG7cuGg+lUivtYTa+fYTG3l0TQULTh7Nf155GtnpOkZB4ktvvqMnAROAd2fzhcBaM5sLVAJFnbYtjIwdxd3vAe4BKCkp8V7kEInKwYZWvvDrNawsP8DNF0zh5gumkKT1eIlDPS56d98AjHz3sZmVAyXuXmNmy4Avm9nDhHfCHtL6vMSinfsbuO6Xq6isbeK/rp7NR087KehIIv2my6I3s6XAeUC+mVUAt7v7fcfZ/FngEqAUaASu76OcIn1m7a6D3LhkNR3uPHjjWcwp1hknJb5156ibq7t4vrjTfQduij6WSP94ces+/v7BtYzKyWDx9XOZkJ8ddCSRfqe9TpIwnlhXwT/99k1mjBnK4uvnkj8kPehIIgNCRS8J4b5Xy/j+05s5e9IIfnHNmQzNSA06ksiAUdFLXHN37nqhlDuff4uLZ43mzk+eTkaqTkomiUVFL3HL3fn3P27j5y+9zRVnFPJvnzhVpzOQhKSil7jk7nz/6S3c/5cyrp47jh9ePkvHyEvCUtFL3HF3vrdsE0te38l1Zxdz+0dncsSpOkQSiope4oq784NntrDk9Z3c+DcT+PalM1TykvB02XqJG+7OHX/Yyn2vlnHd2cUqeZEIFb3EjTufe4tfvLyDT581Tss1Ip2o6CUu3PvnHdz1YilXlhTy/YWzVPIinajoZdB7ZPVufvDMFi45ZTT/8vFTdXSNyBFU9DKo/XHTXm597E3mT8nnzk+eruPkRY5BRS+D1ood+/nK0nWcVpTLf3/mTNJT9I5XkWNR0cugVFpVx+d/tZqi4Zn88ro5uiqUyAmo6GXQqTrczLX3ryI9NZnF188lNyst6EgiMU1FL4NKQ0uIzy1ZxcHGVu6/dg5FebqIt0hX9PeuDBrtHc7ND69jy5467v1sCacUDgs6ksigoBm9DBr/9oetPL+liu99dCYfmj6y6xeICNCNojez+82sysw2dhr7dzPbamZvmtkTZpbb6bnbzKzUzLaZ2Uf6K7gklt+u3s0vXtnBNfPGc80HioOOIzKodGdGvxhYcMTYc8Asdz8VeAu4DcDMZgJXASdHXvNzM9MxbxKVVeUH+NYTGzhn8gi++9GZQccRGXS6LHp3fwU4cMTYn9w9FHn4BlAYub8QeNjdW9y9DCgF5vZhXkkw79Q28cVfr6FweBY//9SZpCZrtVGkp/rip+ZzwO8j98cCuzs9VxEZO4qZLTKz1Wa2urq6ug9iSLxpbmvnSw+sobmtnf/57JkMy9J1XkV6I6qiN7NvAyHgwZ6+1t3vcfcSdy8pKCiIJobEqX/+3SbWVxziP688nckjhwYdR2TQ6vXhlWZ2HXAZcIG7e2S4EijqtFlhZEykRx5euYulK3fz9+dNYsGs0UHHERnUejWjN7MFwDeAj7l7Y6enlgFXmVm6mU0ApgAro48pieTNilq++9Qm5k/J55aLpgUdR2TQ63JGb2ZLgfOAfDOrAG4nfJRNOvBc5Lzfb7j7F919k5k9AmwmvKRzk7u391d4iT+Hm9v48kPryB+Sxk+vmq2zUYr0gS6L3t2vPsbwfSfY/ofAD6MJJYnJ3bn1sTeprG3ikS/MIy9b57AR6Qs6Vk1ixq/f2MmzG/by9Y9M48zxeUHHEYkbKnqJCRsrD/GDp7dw3rQCFs2fGHQckbiiopfANbaG+OrSdeRlp/HjK0/XpQBF+pjOXimB+/7TWyjb38CDN56ldXmRfqAZvQTqT5v2snTlLhbNn8jZk/KDjiMSl1T0EpiqumZufXwDM8fk8I8XTQ06jkjcUtFLINydr//2TRpaQtx19em6sLdIP1LRSyCWrtzNy29V861LZug8NiL9TEUvA67iYCM/fGYzZ08awTXzxgcdRyTuqehlQIXf/boBgH+94lQdSikyAFT0MqAeWrmLV0truO2SGRTlZQUdRyQhqOhlwOw+0MiPntnC30zO59NnjQs6jkjCUNHLgHB3vvNk+Pryd1xxCpGznorIAFDRy4BYtv4dXn6rmq9/ZBqFw7VkIzKQVPTS7w42tPJ/f7eZ04pyueYDxUHHEUk4OteN9LsfPbuFQ01tPPDxU3QhEZEAaEYv/eq10hp+u6aCz587kRljcoKOI5KQuix6M7vfzKrMbGOnsTwze87Mtkduh0fGzczuMrNSM3vTzM7oz/AS21pC7XznyY2MH5HFzRdMCTqOSMLqzox+MbDgiLFbgRfcfQrwQuQxwMWELwg+BVgE3N03MWUwuu/VMnbUNPDPHzuZjFSdy0YkKF0Wvbu/Ahw4YnghsCRyfwlweafxX3nYG0CumY3pq7AyeFTWNvFfL5Ry0cxRnDdtZNBxRBJab9foR7n7nsj9vcCoyP2xwO5O21VExiTB/PCZzXS4838umxl0FJGEF/XOWHd3wHv6OjNbZGarzWx1dXV1tDEkhry6vYZnN+zlpg9N1mkORGJAb4t+37tLMpHbqsh4JVDUabvCyNhR3P0edy9x95KCgoJexpBY0xrq4PZl4R2wi87VRb5FYkFvi34ZcG3k/rXAU53GPxs5+mYecKjTEo8kgF+/sZO3qxv47mUztQNWJEZ0+YYpM1sKnAfkm1kFcDtwB/CImd0A7ASujGz+LHAJUAo0Atf3Q2aJUbWNrdz1wnbmT8nn/OnaASsSK7osene/+jhPXXCMbR24KdpQMjj95Pnt1DW38Z1LZ+qkZSIxRO+MlT7xdnU9D7yxk6vmjmPaaF0aUCSWqOilT/zLs1vISE3maxdODTqKiBxBRS9R+0tpDc9vqeKmD02mYGh60HFE5AgqeolKR4dzx++3MjY3k+vPKQ46jogcg4peovL7jXvZUHmIr314qg6nFIlRKnrptbb2Dv7jT9uYOmoIfztbZ7oQiVUqeum1R9dUUFbTwNc/Ml0XFBGJYSp66ZXmtnZ+8vxbnDEulwtn6M1RIrFMRS+9suS1cvYdbuGbC6brzVEiMU5FLz1W3xLi7pff5rxpBZw1cUTQcUSkCyp66bElr5VT29imN0eJDBIqeumR+pYQ//PnHZw/fSSnFeUGHUdEukFFLz3y7mxeF/sWGTxU9NJtms2LDE4qeuk2zeZFBicVvXSLZvMig5eKXrpl6Ypd1Da28ZXzJwcdRUR6SEUvXWoJtXPvqzs4e9IIZo8bHnQcEemhqIrezL5mZpvMbKOZLTWzDDObYGYrzKzUzH5jZml9FVaC8eS6SvYdbuGLH5wUdBQR6YVeF72ZjQW+CpS4+ywgGbgK+FfgTnefDBwEbuiLoBKM9g7nF6/s4OSTcpg/JT/oOCLSC9Eu3aQAmWaWAmQBe4DzgUcjzy8BLo/ya0iAntu8lx3VDXzpvEk6p43IINXronf3SuA/gF2EC/4QsAaodfdQZLMKQCcqH6TcnbtfepvxI7K4eNaYoOOISC9Fs3QzHFgITABOArKBBT14/SIzW21mq6urq3sbQ/rR62/vZ33FIRadO1HnmxcZxKJZurkQKHP3andvAx4HzgFyI0s5AIVA5bFe7O73uHuJu5cUFBREEUP6y//8eQf5Q9K44ozCoKOISBSiKfpdwDwzy7Lw4u0FwGZgOfCJyDbXAk9FF1GCsKO6nuXbqvnMvPG6FqzIIBfNGv0Kwjtd1wIbIp/rHuCbwD+aWSkwArivD3LKAFv8WjlpyUl8+qzxQUcRkSildL3J8bn77cDtRwzvAOZG83klWIea2nh0TQUfPe0kCoamBx1HRKKkd8bKUR5ZtZvG1nauP6c46Cgi0gdU9PI+ofYOFr9WztwJecwaOyzoOCLSB1T08j7Pb9lHZW0TnztnQtBRRKSPqOjlfe7/SzmFwzP58MxRQUcRkT6iopf3vLWvjpVlB7hm3ni9QUokjqjo5T0PrdhFWnISf1dSFHQUEelDKnoBoKm1ncfWVrBg1mjysnVmaZF4oqIXAJ5+8x3qmkN86qxxQUcRkT6mohcAHlq5i0kF2Zw1IS/oKCLSx1T0wpY9h1m3q5ar547TOedF4pCKXsI7YVOS+MSZOkulSDxS0Se4xtYQT66r5NJTxpCbpZ2wIvFIRZ/gnt2wl7qWEFfP1U5YkXilok9wj6+tYPyILOYUDw86ioj0ExV9AqusbeL1Hfv5+OxC7YQViWMq+gT25LpK3OHjZ+j67SLxTEWfoNydx9ZWMHdCHkV5WUHHEZF+pKJPUOsrDrGjuoErNJsXiXtRFb2Z5ZrZo2a21cy2mNkHzCzPzJ4zs+2RW+3li0GPrakgPSWJS04ZE3QUEeln0c7ofwr8wd2nA6cBW4BbgRfcfQrwQuSxxJCWUDu/e/MdPnLyaIZmpAYdR0T6Wa+L3syGAecC9wG4e6u71wILgSWRzZYAl0cbUvrW8q3V1Da2cYXeCSuSEKKZ0U8AqoFfmtk6M7vXzLKBUe6+J7LNXuCYlyoys0VmttrMVldXV0cRQ3rqyXWVFAxN55xJI4KOIiIDIJqiTwHOAO5299lAA0cs07i7A36sF7v7Pe5e4u4lBQUFUcSQnmhoCbF8WxWXzBpNSrL2xYskgmh+0iuACndfEXn8KOHi32dmYwAit1XRRZS+9OLWKlpCHdoJK5JAel307r4X2G1m0yJDFwCbgWXAtZGxa4GnokooferZDXsoGJpOSbHOOy+SKFKifP1XgAfNLA3YAVxP+JfHI2Z2A7ATuDLKryF9pLE1vGxzZUmRLv4tkkCiKnp3/ytQcoynLojm80r/eHFrFc1tWrYRSTTaG5dAnt2wh/wh6czRso1IQlHRJ4jG1hAvbq1iwaxRWrYRSTAq+gTx0rZqLduIJCgVfYJ4ZsMe8oekcdYEvUlKJNGo6BNAc1s7y7dWcdHJo7VsI5KAVPQJYGXZARpb2/nwjGOejUJE4pyKPgEs31ZFekoS8yZq2UYkEanoE8BL26r5wKQRZKYlBx1FRAKgoo9zZTUNlNU08KFpI4OOIiIBUdHHuZe2hc8pp6IXSVwq+ji3fFs1EwuyGTdCFwAXSVQq+jjW2BrijR37NZsXSXAq+jj2+tv7aQ11qOhFEpyKPo4t31ZFVloycyYMDzqKiARIRR+n3J3lW6s5Z3I+6Sk6rFIkkano41RpVT2VtU1athERFX28erW0BoBzp+YHnEREghZ10ZtZspmtM7OnI48nmNkKMys1s99ELjMoA2xV+QHG5mZSOFyHVYokur6Y0d8MbOn0+F+BO919MnAQuKEPvob0gLuzsuwgc4q1E1ZEoix6MysELgXujTw24Hzg0cgmS4DLo/ka0nPl+xupqW9hzgRdMlBEop/R/wT4BtAReTwCqHX3UORxBTA2yq8hPbSq7AAAc3VtWBEhiqI3s8uAKndf08vXLzKz1Wa2urq6urcx5BhWlR9geFYqk0cOCTqKiMSAaGb05wAfM7Ny4GHCSzY/BXLNLCWyTSFQeawXu/s97l7i7iUFBQVRxJAjrSo/QElxHuGVNBFJdL0uene/zd0L3b0YuAp40d0/DSwHPhHZ7FrgqahTSrdV1TVTvr9RyzYi8p7+OI7+m8A/mlkp4TX7+/rha8hxrCo7CKAdsSLynpSuN+mau78EvBS5vwOY2xefV3puVfkBMlOTOfmknKCjiEiM0Dtj48zKsgOcMT6X1GT9rxWRMLVBHDnc3MbWvYcpGa9lGxH5Xyr6OLJ250E6HOZqfV5EOlHRx5FV5QdISTJmj8sNOoqIxBAVfRxZv/sQM8bkkJXWJ/vYRSROqOjjyPaqOqaNHhp0DBGJMSr6OHGoqY19h1uYotMeiMgRVPRxorSqHkDntxGRo6jo40RpVR0AU0Zq6UZE3k9FHye276snIzWJscMzg44iIjFGRR8nSqvrmVQwhOQknbFSRN5PRR8ntu+r1/q8iByTij4ONLSEqKxt0hE3InJMKvo48Hb1u0fcaEesiBxNRR8Htu8LF/2UUZrRi8jRVPRxoLS6ntRkY3xeVtBRRCQGqejjwPZ99UzIzyZF56AXkWNQM8SB0qo6vVFKRI6r10VvZkVmttzMNpvZJjO7OTKeZ2bPmdn2yO3wvosrR2pua2fXgUYdWikixxXNjD4E3OLuM4F5wE1mNhO4FXjB3acAL0QeSz8pq2mgw3WOGxE5vl4Xvbvvcfe1kft1wBZgLLAQWBLZbAlwebQh5fi2V+mIGxE5sT5ZozezYmA2sAIY5e57Ik/tBUb1xdeQYyvdV0eSwYT87KCjiEiMirrozWwI8BjwD+5+uPNz7u6AH+d1i8xstZmtrq6ujjZGwtpeVU/xiGzSU5KDjiIiMSqqojezVMIl/6C7Px4Z3mdmYyLPjwGqjvVad7/H3UvcvaSgoCCaGAlte1U9k7Q+LyInEM1RNwbcB2xx9x93emoZcG3k/rXAU72PJyfS1t5BeU2DznEjIicUzVWkzwGuATaY2V8jY98C7gAeMbMbgJ3AldFFlONpaAlx2aljmFOcF3QUEYlhvS56d38VON7Jzy/o7eeV7svNSuMnV80OOoaIxDi9M1ZEJM6p6EVE4pyKXkQkzqnoRUTinIpeRCTOqehFROKcil5EJM6p6EVE4pyFzzsWcAizOmBb0Dl6KB+oCTpEDwy2vKDMA2Gw5QVl7my8u3d5srBoToHQl7a5e0nQIXrCzFYPpsyDLS8o80AYbHlBmXtDSzciInFORS8iEudipejvCTpALwy2zIMtLyjzQBhseUGZeywmdsaKiEj/iZUZvYiI9JPAi97MFpjZNjMrNbNbg85zImZWZGbLzWyzmW0ys5uDztRdZpZsZuvM7Omgs3SHmeWa2aNmttXMtpjZB4LOdCJm9rXI98RGM1tqZhlBZzqSmd1vZlVmtrHTWJ6ZPWdm2yO3w4PMeKTjZP73yPfFm2b2hJnlBpmxs2Pl7fTcLWbmZpY/0LkCLXozSwb+H3AxMBO42sxmBpmpCyHgFnefCcwDborxvJ3dDGwJOkQP/BT4g7tPB04jhrOb2Vjgq0CJu88CkoGrgk11TIuBBUeM3Qq84O5TgBcij2PJYo7O/Bwwy91PBd4CbhvoUCewmKPzYmZFwEXAroEOBMHP6OcCpe6+w91bgYeBhQFnOi533+PuayP36wiXz9hgU3XNzAqBS4F7g87SHWY2DDiX8DWJcfdWd68NNlWXUoBMM0sBsoB3As5zFHd/BThwxPBCYEnk/hLg8gEN1YVjZXb3P7l7KPLwDaBwwIMdx3H+jQHuBL4BBLJTNOiiHwvs7vS4gkFQnABmVgzMBlYEm6RbfkL4m6wj6CDdNAGoBn4ZWW6618yygw51PO5eCfwH4dnaHuCQu/8p2FTdNsrd90Tu7wVGBRmmFz4H/D7oECdiZguBSndfH1SGoIt+UDKzIcBjwD+4++Gg85yImV0GVLn7mqCz9EAKcAZwt7vPBhqIvSWF90TWtRcS/gV1EpBtZp8JNlXPefgQvEFzGJ6ZfZvwcuqDQWc5HjPLAr4FfDfIHEEXfSVQ1OlxYWQsZplZKuGSf9DdHw86TzecA3zMzMoJL42db2YPBBupSxV07Lr8AAABaElEQVRAhbu/+9fSo4SLP1ZdCJS5e7W7twGPA2cHnKm79pnZGIDIbVXAebrFzK4DLgM+7bF9jPgkwhOA9ZGfwUJgrZmNHsgQQRf9KmCKmU0wszTCO7CWBZzpuMzMCK8bb3H3Hwedpzvc/TZ3L3T3YsL/vi+6e0zPNt19L7DbzKZFhi4ANgcYqSu7gHlmlhX5HrmAGN55fIRlwLWR+9cCTwWYpVvMbAHhpciPuXtj0HlOxN03uPtIdy+O/AxWAGdEvscHTKBFH9mh8mXgj4R/MB5x901BZurCOcA1hGfFf418XBJ0qDj1FeBBM3sTOB34UcB5jivyl8ejwFpgA+Gfq5h796aZLQVeB6aZWYWZ3QDcAXzYzLYT/svkjiAzHuk4mX8GDAWei/wM/negITs5Tt7A6Z2xIiJxLuilGxER6WcqehGROKeiFxGJcyp6EZE4p6IXEYlzKnoRkTinohcRiXMqehGROPf/AVoScrxYvKYAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(R,np.sqrt(C_integ_I*integ_array))\n",
    "plt.xlim(0.,15.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x117720e50>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEOCAYAAADGy2O9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd8FVX6x/HPkw4JSUihhk4ooUNoCiICUhXsHURX7LKr/hR1V9S164q6gIoioqKAKIIiKCCIdIL0HnpoCUkINf38/pjBDRgggSQzN3ner9e8cu+5M/c+AcI3Z+bMOWKMQSmllHKal9MFKKWUUqCBpJRSyiU0kJRSSrmCBpJSSilX0EBSSinlChpISimlXEEDSSmllCtoICmllHIFDSSllFKu4ON0AZ4kIiLC1K5d2+kylFLKo6xcufKwMSbyQvtpIBVC7dq1iYuLc7oMpZTyKCKyuyD7eeQpOxHZJSLrRGS1iMTZbWEiMltEttlfK9rtIiLvi0i8iKwVkdZ53meQvf82ERnk1PejlFLKQwPJ1tUY09IYE2s/HwbMNcZEA3Pt5wC9gWh7GwJ8AFaAAcOB9kA7YPjpEFNKKVXyPDmQztYfGG8/Hg8MyNP+ubEsBUJFpCrQE5htjEkxxqQCs4FeJV20Ukopi6deQzLALyJigI+MMWOAysaYA/brB4HK9uPqwN48xybYbedqV0qpP2VlZZGQkEB6errTpbheQEAAUVFR+Pr6XtTxnhpInYwx+0SkEjBbRDbnfdEYY+ywumQiMgTrVB81a9YsirdUSnmQhIQEKlSoQO3atRERp8txLWMMycnJJCQkUKdOnYt6D488ZWeM2Wd/TQSmYl0DOmSfisP+mmjvvg+okefwKLvtXO1nf9YYY0ysMSY2MvKCoxaVUqVMeno64eHhGkYXICKEh4dfUk/S4wJJRAJFpMLpx8DVwHpgOnB6pNwgYJr9eDow0B5t1wFIs0/t/QxcLSIV7cEMV9ttSil1Bg2jgrnUPydPPGVXGZhqf+M+wFfGmFkisgKYLCL3AruBm+39fwL6APHASWAwgDEmRUT+Dayw93vJGJNSct+GsxKPprP10HESj6WTeCwDAWqEladmWHnqRgZS3s8T/2kopTyZx/2vY4zZAbTIpz0Z6JZPuwEePsd7fQp8WtQ1utWho+n8uPYAM9cdIG536jn38/PxolujSlzbohpdG1UiwNe7BKtUShWn0zf4R0REFHqfghx7KTwukFThpZ3KYvS8eMYt3kVmdi6NqwbzRI8GtK0TRqUK/lQKDiAn17A35SR7U06ybGcKP67dz8z1BwkL9GNot2hub18TX2+PO8OrlPIgGkilmDGGL5fu5j+zt5J2KovrW0XxcNd61I0Mynf/kOohNK0eQu9mVfln38Ys3p7MB/O3M3z6Bj5bvItn+zSmR0zlfI9VShWPXbt20atXLzp06MDixYtp27YtgwcPZvjw4SQmJjJhwgTatWtHSkoK99xzDzt27KB8+fKMGTOG5s2bk5yczG233ca+ffvo2LEj1kkjy5dffsn7779PZmYm7du3Z/To0Xh7O3dGRAOplEo7lcX/fbOGXzYe4vL64TzbpzFNqoUU+Hgfby+uaBBJ5+gI5m1J5NWfNnPf53Hc2rYGw69pQjk/PY2nyp4Xf9jAxv1Hi/Q9Y6oFM/yaJufdJz4+nm+++YZPP/2Utm3b8tVXX7Fw4UKmT5/Oq6++yvfff8/w4cNp1aoV33//Pb/++isDBw5k9erVvPjii3Tq1Innn3+eGTNmMHbsWAA2bdrEpEmTWLRoEb6+vjz00ENMmDCBgQMHFun3VxgaSKXQ+n1pPDThD/YfOcW/+sVwz+UXf/+EiHBVo8pcER3JiDlbGT1/O3/sSWXU7a2JrlyhiCtXSuWnTp06NGvWDIAmTZrQrVs3RIRmzZqxa9cuABYuXMi3334LwFVXXUVycjJHjx5lwYIFfPfddwD07duXihWtGdLmzp3LypUradu2LQCnTp2iUqVKJfydnUkDqZSZtzmRB75cSVigH5Pu70ibWkUzPZ+Ptxf/17MR7euE849Jq7l25CJG39marg2d/QesVEm6UE+muPj7+//52MvL68/nXl5eZGdnX9R7GmMYNGgQr732WpHUWBT0KnUp8sOa/dz3eRzRlYP44dFORRZGeV3RIJKfhnambmQg942P4/tVf7mXWCnlgM6dOzNhwgQA5s+fT0REBMHBwVxxxRV89dVXAMycOZPUVGuEbbdu3ZgyZQqJidYcAikpKezeXaBVIoqNBlIp8fXyPTw2cRWta1bkq/s6EBHkf+GDLlLl4AAmDulA29ph/H3SasYu3Flsn6WUKpgXXniBlStX0rx5c4YNG8b48dZc08OHD2fBggU0adKE77777s8p0GJiYnj55Ze5+uqrad68OT169ODAgQPn+4hiJ3lHXKjzi42NNW5coG/Sij08/e06ujaMZPQdbUpswEF6Vg5/n7iaWRsO8s++jflb57ol8rlKlaRNmzbRuHFjp8vwGPn9eYnIyjxLBZ2T9pA83I9r9zPsu3Vc0SCSD+8quTACCPD1ZtQdrenTrAovz9jE5BV7L3yQUkqdgw5q8GDzNify94mria1VkY/ubIO/T8kPxfb2Ekbc0pLjGSsZ9t1aggJ86NOsaonXoZTyfNpD8lBxu1J44MuVNKpagbF3t3X0viB/H28+vLM1rWtWZOjEVSzZnuxYLUopz6WB5IG2HTrGvePjqBZajvGD2xEccHGLYRWl8n4+jL27LbXCA3lwwkp2J59wuiSllIfRU3Ye5kDaKQZ9uhw/Hy8+v6cd4Zc6mi47ExI3QOJmSNkBKdvh6H44mQKnUiHzOJhcMAa8vME/GAKCoVxFCImC0JoQWgsqNyGkUmM+GRjLgNGLuHd8HFMfuowKLghLpZRn0EDyIGkns7j70xUcTc9m4pAO1AgrX/g3yUqHPYshfi7sXQYH1kJOhvWaeEFIDStkIhtC+TDwC7LaRSA3BzKOQvpROJkM+1bCxumQm/Xn29cOq8svtZvz3tYIXv08iZfvvR5vnZRVKVUAGkge4mRmNveMX8GOw8f5bHA7mlYv+Lx0pKfBlpmwYSrs+A2yT4G3H1RvA+3ug+qtoUpzq6fj41e4wnJz4Og+OLQBDq6HA6uptGchr/gchv2fcPSNZwlu1geir4Z6XcG3XOHeXylVpIpq+YmgoCCOHz9epLVpIHmAzOxcHvzyD1btSWXk7a25vH4B1iLJyYbtc2HVF7D1F6sXFFID2gyCet2g9uXgF3jpxXl526ftakLD3labMZAcz7fffk1Qwm90XzMJ75XjrN5Wg17Q5DqI7gE+xXfzrlLK82gguVxOruHxyav5bWsSb9zQ7MJDqo/uh7hxsOpLOLYfAiOh7b3Q5HqIirVOvRU3EYiIps/g5xgw6kqeP3acn64XwnfNhE0/wPop1jWoZjdBqzuh6l/WW1RK2XT5CeUKWTm5PDF5DT+uPcAzvRtxS9ua595573JY9iFsnGadRqvfHfq8afVIvJ0ZWFDOz7px9tqRC3lgcTBf3/cuPn3/Y502XD0BVo6H5WOgaktoNwSaXq+n9JS7zRwGB9cV7XtWaQa9Xz/vLrr8hHJUelYOj3z1B3M2JTKsdyPu71Lvrzvl5sDmH2HxSEhYDv4h0P4BaPs3CKtT8kXno36lIF67vhlDJ67m3TnbeLJnQ4jubm0nU2DdFFjxCUx7CH55DtoMhvb3Q4UqTpeulGvo8hPKMcfSs7j/i5Us2ZHMvwc05a4Otc7cIeO41cNYOhpSd1mDEXq/CS3vAP/8V4N1Uv+W1fl922FGz4+na6NI2tQKs14oHwbth1gDK3b9Dss+goUjYMlIaH4LXD4UIqKdLV6pvC7QkykuuvyEcsTahCP0++9Clu1M4Z2bW5wZRkf2wpwXYEQMzHzKuj508+fw2CqrV+HCMDpt+DUxVAstxz8mreF4xlk/QCJQ5wq4dQI8uhJa3QXrvoGRbeGbwXBoozNFK+VBSsPyE9pDconcXMOni3byxqzNRAb5/7m8A7m5sHM+LP8Ets60dm7UDy57FGq0c7TmwqgQ4MuIW1pyy0dLeOmHDbx54zkGMoTXg37vQNdnrZ7S8o9hw3fQ+FqrrZLOuqxUfl544QXuuecemjdvTvny5c9YfuK2226jSZMmXHbZZfkuP5Gbm4uvry+jRo2iVq1a5/uYYlWml58QkV7Ae4A38Ikx5rz98eJYfsIYw+yNh3h3zjY2HjjK1TGVefOGZoQe2wprJ8P6b637fMpHQOuBEDvYGmLtod76eTOj5m3no7va0LNJAa4TnUyBpR9YW+ZxaH4zXDkMwnSpC1UydPmJwrmU5SfKbA9JRLyBUUAPIAFYISLTjTElcn4o8Vg6czcl8uXS3WzYn0abiulM7JRG+9yfkTG/Qtpe8PKxRstd/W+rV1QK7tsZ2q0B87ck8dzU9bSrHUbFwAvciFs+DK56Djo8CIvehWVjrJCOvRe6PAWBBbgnSynlEcpsIAHtgHhjzA4AEZkI9AeKPJBSk5PYk5DAwcSDJB1O4tCBBLJS9lCNw7zgn0Tz4N34n0qFOKy54upcAZ2fsE5TBYYXdTmO8vPx4q0bW3DtyIW89ONGRtzSsmAHlg+DHi9Bh4fgtzeskXmrv4JOQ6HjIzpcXKlSoCwHUnUg74pyCUD74vigHZ/eQ5sTCzjjqokP5PiF4BVeG6nSF6q0gGotoVpr8C7dfy0x1YJ5qGt93p+7jX7Nq9KtceWCH1yhCvQbAe0fhLkvwq8vW/czdX8Bmt5QMjf+qjLHGIPov60LutRLQKX7f74iICJDgCHAnxcDCyvo8iFsOtqPiPBIwsIj8Q4Mh5AovAOCi7JUj/JI1/r8suEgz05dxy+1wwgpV8ibdyMbWKPydi2EWc/At/daw8b7vAnVWhVP0apMCggIIDk5mfDwcA2l8zDGkJycTEBAwEW/R5kd1CAiHYEXjDE97efPABhjzjkovzgGNZRl6xLSGDB6ETe2juKNG5tf/Bvl5lin7+a+BCeSrMEf3Z7X60uqSGRlZZGQkEB6errTpbheQEAAUVFR+Pqe+QtmQQc1lOVA8gG2At2AfcAK4HZjzIZzHaOBVPRe+2kTHy3YwaQhHWhf9xKvl6WnwW9vWlMo+QXCVf+C2HusCWCVUo4paCCV2RtjjTHZwCPAz8AmYPL5wkgVj6Hdo6kRVo5npq4jIzvn0t4sIAR6vgIPLrYmbP3pSfj4KkhYWTTFKqWKVZkNJABjzE/GmAbGmHrGmFecrqcsKu/nw8sDmrEj6QQfzN9eNG8a2RAGTocbxsKxg/BJN5jxhNWDUkq5VpkOJOUOXRpEcm2Laoyet534xCJa8EsEmt0Ij6ywplWK+9Saimj9t9Z6TUop19FAUq7wr34xlPPz5p/fr7vkoaNnCAiG3m/A3+ZaQ8an3ANf3QxH9hTdZyilioQGknKFyAr+PNWrIUt3pDBt9f6i/4DqreFvv0LP12DXIhjVAZaMtkboKaVcQQNJucatbWvSIiqEl2ds4mh6VtF/gLcPdHwIHl4KtS6Dn5+BsT3gkI5lUcoNNJCUa3h7CS8PaEbyiQze+WVr8X1QaE244xtr0EPqLvjoCvj1FcjOKL7PVEpdkAaScpVmUSHc2b4Wny/Zxfp9xTgq7vSgh4dXWFMOLXjTCqYEvc9MKadoICnXefLqhoQF+vH8tPXk5hbziLjAcLh+DNwxBTKOWafwfn4OMk8W7+cqpf5CA0m5Tkh5X57q1Yg/9hxh6qp9JfOh0T3goaXQepC1MOCHnWD3kpL5bKUUoIGkXOrG1lG0rBHKazM3c6w4BjjkJyAYrnnXuqk2NwvG9YaZwyDzRMl8vlJlnAaSciUvL+HFa5uQfCKD9+duK9kPr9sFHlwC7e6DZR/AB5fD7sUlW4NSZZAGknKtFjVCublNDcYt2kV84rGS/XD/IOjzFtw9AzAwrg/MfFp7S0oVIw0k5Wr/16sh5fy8efGHjUU7g0NB1e5kTdba7j5rFnG9tqRUsdFAUq4WEeTPP7o34Pdth5m7KdGZIvwC/9dbys2xri3NekZH4ilVxDSQlOvd1bEW9SsF8fKMjZe+RMWlON1bansvLB1t9Zb2LHOuHqVKGQ0k5Xq+3l78q18Mu5JP8tmiXc4W4x8Eff8DA6dBTiaM6wW//AuydDVRpS6VBpLyCF0aRNKtUSX++2s8icdc8J9/3Sut3lKru2Dx+9YsD/t0IUClLoUGkvIYz/VtTEZ2Dm/N2uJ0KZaAYLj2fbjzW2uWh096wNyXdE48pS6SBpLyGHUjg7j7stpM+SOheOe5K6z63eGhJdDiVvj9PzCmKxxY43RVSnkcDSTlUR7tFk1YeT9e+tGhYeDnUi4UBoyG2ybBycPw8VUw/3XIKaFZJpQqBTSQlEcJDvDl8asbsHxnCrPWH3S6nL9q2MuaE6/J9TD/NSuYDq53uiqlPIIGkvI4t8TWoGHlCrw6cxPpWS5c8bV8GNzwMdzyJRw7AGOuhAVvQU6205Up5WoaSMrj+Hh78c9+jdmbcorPFu9yupxza3wNPLTM+vrryzC2OyRucroqpVzLYwJJRF4QkX0istre+uR57RkRiReRLSLSM097L7stXkSG5WmvIyLL7PZJIuJX0t+PujSdoyPp3rgSI3+NJ+mYi0e1BYbDTePgps/gyB5rePjv72hvSal8eEwg2UYYY1ra208AIhID3Ao0AXoBo0XEW0S8gVFAbyAGuM3eF+AN+73qA6nAvSX9jahL92yfxqRn5TBiTjEud15Umlxn9ZYa9IK5L8KnV0PiZqerUspVPC2Q8tMfmGiMyTDG7ATigXb2Fm+M2WGMyQQmAv1FRICrgCn28eOBAQ7UrS5R3cgg7upYi4nL97DlYAnPBn4xgiLh5s/hxk8hZSd81Fl7S0rl4WmB9IiIrBWRT0Wkot1WHdibZ58Eu+1c7eHAEWNM9lntygMN7RZNhQBfXp7hsmHg5yICTW+Ah5dBg55Wb2lsD722pBQuCyQRmSMi6/PZ+gMfAPWAlsAB4D8lVNMQEYkTkbikpKSS+EhVCKHl/RjaLZrftx1m/lYP+vsJqgQ3fwE3joMju61rSwve0vuWVJnmqkAyxnQ3xjTNZ5tmjDlkjMkxxuQCH2OdkgPYB9TI8zZRdtu52pOBUBHxOav9XDWNMcbEGmNiIyMji+YbVUXqzg61qBMRyCszNpGVk+t0OQUnAk2vh4eXQ6O+1ki8j6+Cg+ucrkwpR7gqkM5HRKrmeXodcPpuw+nArSLiLyJ1gGhgObACiLZH1PlhDXyYbqzzOvOAG+3jBwHTSuJ7UMXDz8eLZ/s0Jj7xOBOX73G6nMILjLBG4d38xf/uW/r1FZ0TT5U5HhNIwJsisk5E1gJdgX8AGGM2AJOBjcAs4GG7J5UNPAL8DGwCJtv7AjwNPC4i8VjXlMaW7Leiilr3xpXoWDecEXO2cTTdQ097xVxr9Zaa3ggL3oSPukCCziCuyg7xiAvBLhEbG2vi4uKcLkOdw/p9aVwzciFDrqjLM70bO13Opdn6C/z4d6vH1OEh6Poc+JV3uiqlLoqIrDTGxF5oP0/qISl1Xk2rh3BD6yjGLdzF3hQPX168wdXWnHht7oYlI+GDjrBjvtNVKVWsNJBUqfLk1Q3x9hJen1UKbjoNCIZ+I+DuGSBe8Hl/mPYwnEp1ujKlioUGkipVqoQEcH+XusxYe4C4XSlOl1M0aneyVqe9/O+w+msY2Q42TAU93a5KGQ0kVeoMuaIulYP9+feMTeTmlpL/tH3LQY8XYcg8CK4K39wNX98KaQlOV6ZUkdFAUqVOeT8f/q9nI9bsPcIPa/c7XU7RqtoC/vYrXP0y7FwAo9rD0g8h14XLcChVSBpIqlS6vlV1mlYP5s1ZW9y5ZtKl8PaByx61lk2v0R5mPQ2fdNNl05XH00BSpZKXl/Bcnxj2HTnF2IU7nS6neFSsDXd+CzeMhbR9MKYr/PwcZHjARLNK5UMDSZVaHeuFc3VMZUbPc/maSZdCBJrdCI8sh9Z3WUPER7WHTT/ooAflcTSQVKn2TJ/GZGTn8s7sLU6XUrzKVYRr3oN7Z1uPJ90JX91iLXOhlIfQQFKlWp2IQAZ2rM2kFXvZdOCo0+UUvxrtYMhv1qCH3YtgdAf47S2dF095BA0kVeoN7RZNcDlfXpmxyTPWTLpUpwc9PLzcWqF23sswuiPEz3G6MqXOSwNJlXoh5X0Z2i2ahfGHmbcl0elySk5Idbh5vDXwAeDLG2DiHXDEA2dEV2WCBpIqE+7sUIu6EYG87GlrJhWF+t2tIeLdnof4udZMD7+9CVnpTlem1Bk0kFSZ4OttrZm0I+kEE5budrqckufjD52fgEdWWEunz3sFRrWDTT/qaDzlGhpIqszo1rgSnepHMGLONo6czHS6HGeE1rBO4w2cDr7lYdId8MUAOLTR6cqUKnwgicjTxVGIUsVNRPhnv8YcS8/i3TnbnC7HWXW7wAO/Q+83Yf8q+LATzHgSTpaSCWmVR7pgIInI5DzbN8DfSqAupYpFoyrB3NauJl8u3U184nGny3GWty+0vx8eXQWxgyFuLLzfEpZ+ADkeuuqu8mgF6SEdNcbcbG83ATp2VHm0x3s0oJyvN6/+tMnpUtwhMBz6/gceWATVWsOsYdb9S5tn6PUlVaIKEkivnPX8ueIoRKmSEh7kz2Pdovl1cyK/bU1yuhz3qBwDd02F2ydbCwJOvB0+62ed0lOqBFwwkIwxOwFEJMJ+rieZlccbeFktaoeX598/bix7w8DPR8QahffgYqvXlLQZxlwJ3/4NUsvg6ERVogozqOHTYqtCqRLm7+PNc31jiE88XjaHgV+Ity+0/Rs89oc1XHzTDzAy1ppNXAc+qGJSmECSYqtCKQd0zzMMPPVEGR0GfiEBIdYNtY/+Ac1vhiWj4L0WsOBtyDzhdHWqlClMIJXI1U0RuUlENohIrojEnvXaMyISLyJbRKRnnvZedlu8iAzL015HRJbZ7ZNExM9u97efx9uv1y6J7025i4jwr34x9jDwrU6X424h1aH/KOtUXq3L4dd/w/utYcUnkK1hroqGG3tI64HrgQVnfLhIDHAr0AToBYwWEW8R8QZGAb2BGOA2e1+AN4ARxpj6QCpwr91+L5Bqt4+w91NlUMMqFbijfS2+XLaHrYd0YbsLqhwDt0+EwbMgrA7MeMI6lbdmoi6jri5ZYQLpmWKrIg9jzCZjTH6L1/QHJhpjMuyBFvFAO3uLN8bsMMZkAhOB/iIiwFXAFPv48cCAPO813n48Behm76/KoMd7NCDI34cXf9hQNmYDLwq1OsLgmXDHFAgIhqn3WzOKb5gKuTpIRF2cAgeSMWa9iCwQkWAAEXlARP5++jRYCagO7M3zPMFuO1d7OHDEGJN9VvsZ72W/nmbvr8qgioF+PN6jAYvik/l5wyGny/EcIhDdA4YsgJs+s9q+uRs+6qxz5KmLUtipg0KMMUdFpA1wH1AR+LiwHyoic0RkfT5b/8K+V3ETkSEiEicicUlJes9KaXVH+5o0rFyBl2dsJD1LTz0VipcXNLnOmlH8+o8h65Q1R54GkyqkwgZSloj4AAOBN4wxw7Gu6RSKMaa7MaZpPtu08xy2D6iR53mU3Xau9mQg1K43b/sZ72W/HmLvn1+tY4wxscaY2MjIyMJ9o8pj+Hh7MfyaGBJST/Hxgh1Ol+OZvLytkXgPL4cBH1qj8CbdAR92ho3T9FSeuqDCBtL7wBqgH/CD3RZUpBWd23TgVnuEXB0gGlgOrACi7RF1flgDH6Yb62LAPOBG+/hBwLQ87zXIfnwj8KvRiwdl3mX1I+jdtAqj5sez/8gpp8vxXN4+0PI2eHiFFUzZp2DyQPjgMlg3BXKyL/weqkwqVCAZYz4H2gNNjTGnRKQ+sKQoCxKR60QkAegIzBCRn+3P3gBMBjYCs4CHjTE59jWgR4CfgU3AZHtfgKeBx0UkHusa0Vi7fSwQbrc/Dvw5VFyVbc/2aYwx6Dx3ReHPYFoON9g/et/ea43KW/kZZGc4Wp5yHylsx0BEfI0xZXIq4NjYWBMXF+d0GaqYvTtnK+/O2cZX97XnsnoRTpdTeuTmwpaf4Pe3rfnxKlSFDg9Bm7utkXqq1BKRlcaY2AvtV6gekoh8AuwRkb32DaUfi8ijF12lUi70QJd6RFUsxwvTN+g8d0XJywsa94P75lmTuEY0gNn/ghFNYc4LcOyg0xUqhxX2GlJnIMoYUwPr5tWpQGCRV6WUgwJ8vXm+XwxbDx3n8yU6z12RE4F6V8Gg6XDfr9ZigQvftYLp+4d19doyrLCBtAz7fh1jzD5jzE/GmNeLviylnNUjpjJdGkTy7uytJB5Ld7qc0qt6G7jlC2sS1zZ3w/pv4YOO8MV1sG2ODhkvYwobSB8Bv4nIkyLSWURCiqMopZwmIgy/Job07Bxen7nZ6XJKv7C60PdteHwjXPUvq5c04QYY1d6aLy+jjK/uW0YUNpC+BD4HfICHgMUisr3Iq1LKBepGBnFf57p898c+lu/UJRdKRPkwuOJJ+Ps6uO4j8C1nzZf3Toy19EWy/ndTmhVqlJ2I/G6M6XxWm78xpkyM39RRdmXPycxseryzgCB/H358rBO+3oX9HU5dEmMgYQUs/QA2TYfcbKjfHdreZ01b5OXtdIWqAIpllB2wWkSG5m0oK2Gkyqbyfj48f00MWw4dY/ziXU6XU/aIQI12cNM4+McGuPIZOLgevr7FXpfpLTim8w+WFoXtIU0GmmFNtfMH1qwNq40x3xRPee6iPaSyyRjDPZ+tYMWuVOY+0YXKwQFOl1S25WRZ9zOtGAs7fwMvH2jQC9oMhnpdtdfkQgXtIRX6xlj7zf2x5rBrBjQzxjxZ+BI9jwZS2bU7+QQ9Rizg6pjKjLy9tdPlqNOSt8PKcbD6KziZDCE1oOUd0PJ2qFjL6eqUrVgCSUReNsb886w2b2NMmZgeWQOpbHt/7jbemb2V8fe0o0sDnWjXVbIzYcsM+ONz2D4PMFCnixVOjfuBn94u6aTiuoZUXURuz/MhkcCcwhanlCe6v0td6kYE8vy09bpEhdv4+FlLYNyzOnTrAAAdaElEQVQ11Rqh1/U5SN0FU4fA2w2sG253/q4zjrtcYQPpfuA+EWknIm2xZtN+u+jLUsp9/H28efm6puxOPsmoefFOl6POJbQGdHkKHlsNd/8ETQbAxu9hfD94txnMHg6HNlz4fVSJK9ApOxH5HGsQwyogCfgKyAZuNcaUmZ9MPWWnAB6ftJof1u5n5tArqF+ppFZfUZck86Q1EGLtJIifCyYHIhtDsxugyfUQXs/pCku1Ir2GJCJXAS3srTlQF2sdokXAemPMlEsr1zNoICmAw8czuOrt+TSuGszEIR0QEadLUoVx4jBsmGpNU7THXj2nagsrmGL6Q1gdZ+srhYo6kDoCS08vYmevstoYO6SMMf93ifV6BA0kddrE5XsY9t063ryhOTe3rXHhA5Q7Hdlrnc5b/x3s/8Nqq9LcCqaY/hAR7Wx9pURRB9IHWAvzbcVaHG+WMabMzRWvgaROy8013PrxUrYcPMacx7sQWcHf6ZLUpUrdBRunWzNCJKyw2iIaQKN+1latlbWEhiq04hr23QjoDfTEujl2HlZALSoLQ781kFRe8YnH6fPe7/RqWoX3b2vldDmqKKXtg80zYPOPsGuhdc0pqLJ1A27D3taQcr/yTlfpMYr1xlj7A8oBXbECqmNBPszTaSCps703Zxsj5mxl3OC2dG1YyelyVHE4mQLxc6yAip8LmcfAJwBqd4boqyG6uzVbuTqnYg+kskgDSZ0tIzuHvu8v5FRmDj//4wqC/H2cLkkVp+wM2L0Ytv0CW2dByg6rPayuNelrvaugdifwr+BsnS5TZIEkIj2Am4FRxpjVIjLEGDOmiOr0KBpIKj8rd6dw44dLGNihFi/2b+p0OaokJW+3ek3xs61Te1knrbn1otpZK+HW6QJRseDt63SljipoIBXk17l7gAeBf4pIGNDyrA96xBgz0n4caow5cjEFK+Wp2tQKY1DH2oxfsot+LarRtnaY0yWpkhJez9raD7F6T3uXwfZfremL5r8O818D30Co2cHqOdW5whpiXsYD6lwK0kMaY4wZYj9+HehmjGmb5/W3Tg/7FpE3jTFPFWfBTtIekjqXExnZ9Hx3AX4+Xvz0WGcCfHXG6TLvZIrVa9q5AHb9Dkn2ysO+gdaSGrUuh5rtrWXcS/lce0XZQ5ohIo2NMZuMMcNE5NGzXvcSkc5YN8nqr4aqTAr09+G165tx19jlvDd3G0/3auR0Scpp5cMg5lprAzieCLsXWdegdi+Gea8AxjrFV6W5FVJRba2vITWstaDKmILeh7QD+A143hiz96zXvLBO6bUGvjfG/HBJBYncBLyAdeNtO2NMnN1eG9gEbLF3XWqMecB+rQ3wGVAO+AkYaowx9inGSUBtYBdwszEmVaxb698D+gAngbuNMX9cqDbtIakL+b9v1vDdqn1MfegymkeFOl2OcrNTqbB3BexdCnuWWTfmZp20XguqbPWcqreGaq2te6DKe+7v+0V9Y6wf1sSqjwPTgZeNMUmXXGX+n9UYyAU+Ap48K5B+NMb85aqxiCwHHgOWYQXS+8aYmSLyJpBijHldRIYBFY0xT4tIH+BRrEBqD7xnjGl/odo0kNSFpJ3KoueIBQSX8+GHRzvh76On7lQB5WRZk74mrICEONi3EpK3/e/1irWt609VW0CVFlClGVSo7Fi5hVGUp+wwxmQC/xWRj4FHgOUi8iXwljHm6KWV+pfP2gQUeH4wEakKBBtjltrPPwcGADOB/sCV9q7jgfnA03b75/ZUSEtFJFREqhpjDhTdd6LKopByvrx2QzMGj1vBe3O28ZSeulMF5e0L1VpaW7v7rLZTR+DAGti/yupBHVgLG6f975jASCuYKjeBSk2srxENwNczVzUu1E0Txph04G17KqGhwEoR+cgYU1JLUNQRkVXAUeCfxpjfgepAQp59Euw2gMp5QuYgcPrXierA3nyO+UsgicgQYAhAzZo1i+jbUKVZ14aVuDk2ig9/207PJlVoUUNP3amLVC7UGj5et8v/2tLTrGA6tB4OrrO2ZR9BTqb1unhZ90VFNoLIhhDRECIbQHg0+Lt7dvpCBZJ92qwR0BCoCRwDXqWQayKJyBygSj4vPWeMmZZPO1hhUdMYk2xfM/peRJoU9DPta0qFvgvYvudqDFin7Ap7vCqb/tkvht+3HeaJb9bw46OddNSdKjoBIVCns7WdlpMNKdutU36Jm6wRfUmbYctMa9qj0ypUhfD61hZW939bxVquGOlXoEASkbVYPYg9wGaswQVzgZFYE64WijGm+0UckwFk2I9Xish2oAGwD4jKs2uU3QZw6PSpOPvUXqLdvg+ocY5jlLpkwQG+vHFDcwZ+upy3f97CP/vFOF2SKs28fazeUGTDM9uzMyF1JyRtsa5HHY63vm783hpUkVdQZes6VWjN/20hNewtqkTm7itoD2kAsNMUZAREMbGXS08xxuSISF0gGthhjEkRkaMi0gFrUMNA4L/2YdOBQcDr9tdpedofEZGJWIMa0vT6kSpqVzSI5K4OtRi7aCfdGlemY71wp0tSZY2PX/5BBVYgpeyAlJ3WTOent73LrOU4zp4v+55frPumipHr5rITkeuwAiUSOAKsNsb0FJEbgJeALKxReMNPDzEXkVj+N+x7JvCofYouHJiMdXpxN9aw7xR72PdIoBfWsO/Bp0fznY+OslOFdTIzm77vLyQzO5dZf+9MhQC9Q195gJxsOLYf0hLsbS+0HgSBERf1djq5ajHQQFIX4489qdz4wWKubx3F2ze1cLocpUpcQQNJV5tSqpi1rlmRh66sz5SVCcxcp2eGlToXDSSlSsDQ7tE0jwph2HfrOJB2yulylHIlDSSlSoCvtxfv3dqKzOxcnpi8htxcPVWu1Nk0kJQqIXUiAnnh2hgWb0/mk4U7nC5HKdfRQFKqBN0cW4NeTarw1s9bWJeQ5nQ5SrmKBpJSJUhEeP2GZkQG+fPI139wLD3L6ZKUcg0NJKVKWGh5P967rRUJqad4bup69NYLpSwaSEo5oG3tMP7RPZrpa/bzzcqECx+gVBmggaSUQx68sj6X1Qtn+LQNbD10zOlylHKcBpJSDvH2Et69pSWB/j48NOEPTmRkO12SUo7SQFLKQZWCA3j/tpbsSDrOs1PX6fUkVaZpICnlsMvqRfB4jwZMW72fCcv2OF2OUo7RQFLKBR66sj5XNozkpR82smbvEafLUcoRGkhKuYCXlzDi5pZEVvDnwS9Xcvh4htMlKVXiNJCUcomKgX58dFcbkk9k8uhXq8jOyXW6JKVKlAaSUi7StHoIr1zXjCU7knlj1many1GqRBV0CXOlVAm5sU0UaxOO8PHvO2laPYT+Las7XZJSJUJ7SEq50D/7xtCuThhPTVmrgxxUmaGBpJQL+fl48cEdrYms4M+QL+I4dDTd6ZKUKnYaSEq5VHiQPx8PjOVYejZDvlhJelaO0yUpVaw0kJRyscZVgxlxS0vW7D3Ck9/oSrOqdNNAUsrlejapwtO9GvHj2gO8M3ur0+UoVWxcF0gi8paIbBaRtSIyVURC87z2jIjEi8gWEemZp72X3RYvIsPytNcRkWV2+yQR8bPb/e3n8fbrtUvye1SqsB7oUpdbYmswcl4838TtdbocpYqF6wIJmA00NcY0B7YCzwCISAxwK9AE6AWMFhFvEfEGRgG9gRjgNntfgDeAEcaY+kAqcK/dfi+QarePsPdTyrVEhJeva8rl9cN5duo6FscfdrokpYqc6wLJGPOLMeb0PPxLgSj7cX9gojEmwxizE4gH2tlbvDFmhzEmE5gI9BcRAa4CptjHjwcG5Hmv8fbjKUA3e3+lXMvX24vRd7ShTkQg93+xko37jzpdklJFynWBdJZ7gJn24+pA3nMVCXbbudrDgSN5wu10+xnvZb+eZu//FyIyRETiRCQuKSnpkr8hpS5FSDlfPhvcjkB/H+4et5yE1JNOl6RUkXEkkERkjoisz2frn2ef54BsYIITNZ5mjBljjIk1xsRGRkY6WYpSAFQLLcf4e9pxKiuHQZ8uJ/VEptMlKVUkHAkkY0x3Y0zTfLZpACJyN9APuMP8b8WyfUCNPG8TZbedqz0ZCBURn7Paz3gv+/UQe3+lPELDKhX4eGAse1NPMfizFbrarCoVXHfKTkR6AU8B1xpj8p6PmA7cao+QqwNEA8uBFUC0PaLOD2vgw3Q7yOYBN9rHDwKm5XmvQfbjG4FfjS7VqTxMh7rh/Pe2Vqzbl8aQL+L0xlnl8VwXSMBIoAIwW0RWi8iHAMaYDcBkYCMwC3jYGJNjXwN6BPgZ2ARMtvcFeBp4XETisa4RjbXbxwLhdvvjwJ9DxZXyJD2bVOHNG5qzKD6Zx77WJSuUZxPtGBRcbGysiYuLc7oMpf5i3KKdvPjDRq5rVZ23b2qBt5cOGlXuISIrjTGxF9pPl59QqhQYfHkdTmRk8/YvW/H1Fl6/vjleGkrKw2ggKVVKPHJVNJk5hvfnbsPH24tXBjRFb69TnkQDSalS5B/do8nKyeWD+dvxFuHFa5toT0l5DA0kpUoREeGpng3JzTV8tGAHWTm5vHpdMw0l5RE0kJQqZUSEYb0b4evtxch58WTlGN68sbkOdFCup4GkVCkkIjzZsyF+Pl68M3srGdk5vHNzS/x83Hinh1IWDSSlSrHHukXj7+PFazM3czwjmw/uaEM5P2+ny1IqX/rrklKl3P1d6vHa9c34bWsSd41dRtqpLKdLUipfGkhKlQG3tavJyNtasybhCLd8tIQDaaecLkmpv9BAUqqM6Nu8KuPubkdC6imuH72YzQd1PSXlLhpISpUhnaIjmHx/R3KN4aYPlujKs8pVNJCUKmNiqgXz3UOXUzU0gIGfLmfi8j1Ol6QUoIGkVJlUPbQcUx68jMvqRzDsu3W89MNGcnJ1omXlLA0kpcqo4ABfPh0Uy+DLa/Ppop3c89kK0k7qCDzlHA0kpcowH28vhl/ThFeva8bi7Ye5dtRCHeygHKOBpJTi9vY1mTikA6cyc7hu1GKmrd7ndEmqDNJAUkoB0KZWGD8+1omm1YMZOnE1z01dp8uiqxKlgaSU+lOlCgF8dV8H7u9SlwnL9nDd6MXsSDrudFmqjNBAUkqdwdfbi2d6N2bc3W05mHaKfv9dyKQVezBGR+Gp4qWBpJTKV9dGlfhpaGdaRIXy9LfreODLlaScyHS6LFWKaSAppc6pakg5JvytPc/1acy8zUn0fHcBszcecrosVUppICmlzsvLS7jvirp8//DlRAT5c9/ncQyduEp7S6rIuS6QROQtEdksImtFZKqIhNrttUXklIistrcP8xzTRkTWiUi8iLwvImK3h4nIbBHZZn+taLeLvV+8/TmtnflulfIcMdWCmfbw5fyjewN+WneAHu/8xver9um1JVVkXBdIwGygqTGmObAVeCbPa9uNMS3t7YE87R8A9wHR9tbLbh8GzDXGRANz7ecAvfPsO8Q+Xil1AX4+XgztHs0Pj3YiKqw8f5+0mjvHLtOReKpIuC6QjDG/GGOy7adLgajz7S8iVYFgY8xSY/2q9jkwwH65PzDefjz+rPbPjWUpEGq/j1KqABpVCea7By/j3wOasjYhjV7v/s6bszZzIiP7wgcrdQ6uC6Sz3APMzPO8joisEpHfRKSz3VYdSMizT4LdBlDZGHPAfnwQqJznmL3nOEYpVQDeXsJdHWox94ku9G1eldHzt9P17fl8uzKBXJ2oVV0ERwJJROaIyPp8tv559nkOyAYm2E0HgJrGmFbA48BXIhJc0M+0e0+F/ikRkSEiEicicUlJSYU9XKlSr1KFAEbc0pLvHrqMqqHleOKbNVwzciG/b9OfF1U4Pk58qDGm+/leF5G7gX5ANztIMMZkABn245Uish1oAOzjzNN6UXYbwCERqWqMOWCfkku02/cBNc5xzNm1jgHGAMTGxuqvfUqdQ+uaFZn64GVMX7Oft3/Zwl1jl9OpfgRP9mxIyxqhTpenPIDrTtmJSC/gKeBaY8zJPO2RIuJtP66LNSBhh31K7qiIdLBH1w0EptmHTQcG2Y8HndU+0B5t1wFIy3NqTyl1kby8hAGtqjP3iS483y+GDfvTGDBqEXePW86qPalOl6dcTtw2ZFNE4gF/INluWmqMeUBEbgBeArKAXGC4MeYH+5hY4DOgHNY1p0eNMUZEwoHJQE1gN3CzMSbFDq6RWKPxTgKDjTFxF6otNjbWxMVdcDellO14RjZfLNnNmAXbST2ZRaf6ETzQpR6X1w/HvjtDlQEistIYE3vB/dwWSG6mgaTUxTmRkc0XS3czduFOko5l0LR6MH/rVJc+zari5+O6EzWqiGkgFQMNJKUuTUZ2Dt+v2sdHC3awI+kEEUH+3N6+Jre3q0mVkACny1N5GGNISD3Fkh3JLN2RzL/6xlAx0O+i3ksDqRhoIClVNHJzDQu2JTF+8S7mb01CgK4NK3FL2xp0bVQJX2/tNZU0Yww7D59g2c4Ulu9MYdmOZPanpQMQHujHJ4NiaVWz4kW9d0EDyZFRdkqpss3LS7iyYSWubFiJ3cknmLRiL9+sTGDu5kQigvzo17waA1pVp0VUiF5rKiaZ2bls2J/Gyt2prNiVQtyuVJLt+QkjgvxoVyeMB+uG06FuOPUrBZXI34P2kApBe0hKFZ/snFzmbUli6qoE5mxKJDM7l5ph5endtAo9m1ahZVQoXl4aThfrQNopVu05wuq9R1i1J5W1CWlkZOcCUCOsHG1rhxFbK4z2dcOoGxFYpAGkp+yKgQaSUiUj7VQWP68/yI/rDrA4/jDZuYZKFfzp2rASXRtVolN0BEH+eoLnXA4fz2D9vjTWJpzejpB4LAMAP28vYqoFE1urIq1rVaR1zYrFfv1OA6kYaCApVfLSTmUxd9Mh5m5KZMHWJI5lZOPjJbSsEcpl9cLpWC+CljVCKefn7XSpJS4n17A7+QSbDx5j4/6jbDpwlA37j3LwaPqf+9SNDKRFVCjNo0JoVbMijatWwN+nZP+sNJCKgQaSUs7Kysll5e5UftuaxOLtyaxLOEKuAR8voUm1YNrUCqN5VAhNq4dQNyKw1Jziy8zOZW/qSeITjxOfeJzticfZmniMbYeO/3nazdtLqBcZSEzVYJpWD6FJtRCaVA8mOMDX4eo1kIqFBpJS7nI0PYs4+4J83O5U1iYcIT3L+g860M+bBlUq0LByBRpUrkDdyEDqRARSPbQcPi4bxWeM4cjJLPYdOcW+I6fYm3KSPSkn2Z18kl3JJ0hIPUVOnglrqwQHEF05iEZVKhBduQKNqwQTXTmIAF939hI1kIqBBpJS7padk0t80nHWJaSxYf9Rthw8xpZDx85Y3dbXW6gaUo7qoeWoFlqOysH+VKrgT2SFACoG+hJazo+Kgb4E+fsQ6OdzUb2snFzDicxsTmRkcyw9m7RTWRw5mUXqiUwOn8gg+Xgmh49ncOhoOoeOZnAwLZ1TWTlnvEcFfx9qhpendkQgdcIDqR0RSP1KQdSLDKSCC3o9haGBVAw0kJTyPMYYkk9ksvPwCXYePsGuw1aPY9+RU+xLPUXS8Ywzeh9nK+frTYCvF34+1ubj5YUInI6pnFxDdq4hKyeX9Kxc0rNy/jyNdi7l/byJCPKncrA/lYMDqBwcQLXQ0yEZQM2w8oSU8y01Q971PiSllAJEhIggfyKC/GlbO+wvr+fmGlJPZpJ0PIOUE5mkncwi9WQWxzOyOJGRw4mMbDKyc8nMziUzJ5fsXIMxhtO/y/t4C95egq+XFwG+XgT4ehPg602FAB8C/X0I8vehYnk/Qsr5Elrel/AgP8r76X+9+dE/FaVUmeblJYQH+RMe5O90KWWeu67sKaWUKrM0kJRSSrmCBpJSSilX0EBSSinlChpISimlXEEDSSmllCtoICmllHIFDSSllFKuoFMHFYKIJAG7C3lYBHC4GMopSlpj0XB7jW6vD7TGouK2GmsZYyIvtJMGUjETkbiCzOHkJK2xaLi9RrfXB1pjUfGEGvOjp+yUUkq5ggaSUkopV9BAKn5jnC6gALTGouH2Gt1eH2iNRcUTavwLvYaklFLKFbSHpJRSyhU0kIqRiPQSkS0iEi8iw5yuJy8RqSEi80Rko4hsEJGhTtd0LiLiLSKrRORHp2vJj4iEisgUEdksIptEpKPTNZ1NRP5h/z2vF5GvRSTABTV9KiKJIrI+T1uYiMwWkW3214ourPEt++96rYhMFZFQt9WY57UnRMSISIQTtRWWBlIxERFvYBTQG4gBbhORGGerOkM28IQxJgboADzssvryGgpscrqI83gPmGWMaQS0wGW1ikh14DEg1hjTFPAGbnW2KgA+A3qd1TYMmGuMiQbm2s+d9Bl/rXE20NQY0xzYCjxT0kWd5TP+WiMiUgO4GthT0gVdLA2k4tMOiDfG7DDGZAITgf4O1/QnY8wBY8wf9uNjWP+JVne2qr8SkSigL/CJ07XkR0RCgCuAsQDGmExjzBFnq8qXD1BORHyA8sB+h+vBGLMASDmruT8w3n48HhhQokWdJb8ajTG/GGOy7adLgagSL+zMevL7cwQYATwFeMxAAQ2k4lMd2JvneQIu/A8fQERqA62AZc5Wkq93sX6ocp0u5BzqAEnAOPu04iciEuh0UXkZY/YBb2P9pnwASDPG/OJsVedU2RhzwH58EKjsZDEFcA8w0+kiziYi/YF9xpg1TtdSGBpIZZyIBAHfAn83xhx1up68RKQfkGiMWel0LefhA7QGPjDGtAJO4PxppjPY12H6Y4VnNSBQRO50tqoLM9YQYNf+di8iz2Gd+p7gdC15iUh54FngeadrKSwNpOKzD6iR53mU3eYaIuKLFUYTjDHfOV1PPi4HrhWRXVinPK8SkS+dLekvEoAEY8zp3uUUrIByk+7ATmNMkjEmC/gOuMzhms7lkIhUBbC/JjpcT75E5G6gH3CHcd+9M/WwfvlYY//sRAF/iEgVR6sqAA2k4rMCiBaROiLih3URebrDNf1JRATruscmY8w7TteTH2PMM8aYKGNMbaw/v1+NMa76zd4YcxDYKyIN7aZuwEYHS8rPHqCDiJS3/9674bKBF3lMBwbZjwcB0xysJV8i0gvrNPK1xpiTTtdzNmPMOmNMJWNMbftnJwFobf9bdTUNpGJiX/R8BPgZ64d/sjFmg7NVneFy4C6sXsdqe+vjdFEe6lFggoisBVoCrzpczxns3tsU4A9gHdbPveN38ovI18ASoKGIJIjIvcDrQA8R2YbVs3vdhTWOBCoAs+2fmw9dWKNH0pkalFJKuYL2kJRSSrmCBpJSSilX0EBSSinlChpISimlXEEDSSmllCtoICmllHIFDSSllFKuoIGklMuIyP0iclBE1ojIdhEZmOe1D0XkcvtxNxH5wrlKlSpaGkhKuU8z4AVjTAvgNiDv1E4dsJY8AGvtpVUlXJtSxUYDSSn3aQ5ssR/vBDIBRKQxsNUYk2O/1gJYJSL+IvKZiLwqlq9FZJKILBeR3SLS1z6+moh8ay+TsVlE2pX0N6bU+WggKeU+zYAt9kSojwDP2e29gVl59muONRv2z8AcY8yz9szTLYAdxph2wB3AcHthvpnAOHuZjNa4d4JVVUbpXHZKuYi97PROYD3Wgo5rge7GGCMiPwODjTH77aVDDgO7gfuNMUvs4wOwFoasYYxJF5EwrIUXnwFuNMa4YelypfKlPSSl3KUZsMAY0xJoADQCOtqLroUaY04vPd4Ya4mTbCAnz/FNgW3GmHT7eWtgDdYs5EtRysU0kJRyl+bYAxWMManAV0BfoCswL89+LYDFWOtEjRORynnaa4pIgL2U+ovACKzlwJucPlhEIov5+1Cq0DSQlHKXZpw5cu4HoA9/vX7UAlhvjNkKPA1Mtk/jtcBaEXYZVg/qA2PMIuAzoLKIbBCR1UDH4v5GlCosvYaklAcQkT+A9vYS5Ofb7zdgiDFmy/n2U8qNNJCUKkVEJAGoaYzJdboWpQpLA0kppZQr6DUkpZRSrqCBpJRSyhU0kJRSSrmCBpJSSilX0EBSSinlChpISimlXEEDSSmllCtoICmllHKF/wc9oHtwIdzPEAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(R,vc2_ISM(R,R_d=R_d_I,Sigma_g_ISM=Sigma_g_ISM_I,R_m=R_m),label='model I')\n",
    "plt.plot(R,vc2_ISM(R,R_d=R_d_II,Sigma_g_ISM=Sigma_g_ISM_II,R_m=R_m),label='model II')\n",
    "plt.xlabel(r'$R/kpc$')\n",
    "plt.ylabel(r'$V_{c}^{2}/km s^{-1}$')\n",
    "#plt.title(r'$V_{c}^{2}$)\n",
    "plt.legend(loc='upper right', prop={'size':10})\n",
    "#plt.savefig('Model_II.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in sqrt\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n",
      "/Users/atalia/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in sqrt\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.0, 15.0)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEaCAYAAAD65pvjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl83XWd7/HXJ0uTtlnaNGm6pPtKV8CClE20bAoCMyOOK8XRAZdxdJxxBmfuFX1cB3XGO+o43plBQVHAEREVFBcom2Uv0A3SJU1bmjZpli5ZumT73D9+JyVkac5Jzjm/c07ez8ejj5Nz8jvn924p/eS7m7sjIiLSW1bYAUREJPWoOIiISD8qDiIi0o+Kg4iI9KPiICIi/ag4iIhIPyoOIiLSj4qDSAzMbI+ZXWpmY8zs/5pZjZm1Rl7/Vp/r2s2stM/7XzEzN7PZyc4uEgsVB5Hh+QKwCjgXKAQuAV7uc81u4P09T8xsOTAuSflERkTFQWR4zgF+4e4HPLDH3X/U55ofAzf0er4W6HuNSEpScRAZnueAz5nZJ81suZnZINcUmdkZZpYNvA+4O6kpRYZJxUFkeL4KfB34ILAB2G9mawe4rqf1cBlQCexPWkKREcgJO4BIOnL3LuC7wHfNbCzwF8CdZvaCu1f2uvTHwFPAHNSlJGlELQeREXL34+7+XeAwsKTP9/YSDEy/C3gghHgiw6KWg8gwmNlngY3A80AHQfdSIfDKAJd/FJjo7m1mpv/nJC3oL6rI8BwD/i8wH3BgB/Bn7l7d90J335XkbCIjZjrsR0RE+tKYg4iI9KPiICIi/ag4iIhIPyoOIiLST9rOViotLfXZs2eHHUNEJK289NJLje5eNtR1aVscZs+ezYYNG8KOISKSVsxsbzTXqVtJRET6UXEQEZF+VBxERKSftB1zEJHRp6Ojg5qaGk6cOBF2lJSXn59PRUUFubm5w3q/ioOIpI2amhoKCwuZPXs2A5+vJADuTlNTEzU1NcyZM2dYn6FuJRFJGydOnGDSpEkqDEMwMyZNmjSiFpaKg4ikFRWG6Iz0z0nFIUk6u7r5484GfvTsHppaT4YdR0TktDTmkGDb6pq5+7m9/HZLHU1t7QDc9nAl7ztnJh+7aA4VE8eFnFBEwtKzmLe0tDTma6J570ioOCTQjoMtXP+fz9LR3c2aM8p594ppzCgZyw+e3sPdz+3lx8/t5Wt/upzrV80IO6qIyJuoWylBGlpO8pEfvEj+mGwe+9tL+O4HzubKZVNYOq2Yb1y/kqf+/u28dU4J//SLrbzy+uGw44pIFPbs2cPixYu58cYbWbhwIR/84Ad59NFHueCCC1iwYAEvvPACAIcOHeK6665jxYoVnHfeeWzevBmApqYmLr/8cpYuXcrHPvYxeh+2dvfdd3Puuedy5plncvPNN9PV1RXK77GHWg4JcKKji5t+vIGmtpPcd/Nqpk0Y2++aaRPG8t0PnM27/2M9n7j7ZR769IWUFeaFkFYkPX35oVd57UBzXD9zybQibn330tNeU1VVxc9+9jPuvPNOzjnnHO69917Wr1/Pgw8+yG233cYvf/lLbr31Vs466yx++ctf8thjj3HDDTewceNGvvzlL3PhhRfyxS9+kd/85jfccccdAFRWVvLTn/6Up59+mtzcXD75yU9yzz33cMMNN8T19xcLtRzizN35/P2beeX1I3zrz89kRcWEQa+dOH4M//Wht3D4WDt/de/LdHR1JzGpiAzHnDlzWL58OVlZWSxdupQ1a9ZgZixfvpw9e/YAsH79ej784Q8D8I53vIOmpiaam5t56qmn+NCHPgTAVVddxcSJEwFYt24dL730Eueccw5nnnkm69ato7q633HkSaWWQ5w9WlnPQ5sO8PkrFnHlsqlDXr9sejFf+7Pl/M1PN/Gvv9/OP77rjCSkFEl/Q/2Enyh5eW+08LOysk49z8rKorOzc1if6e6sXbuWr371q3HJGA9qOcSRu/Odx3Yyo2QsN108N+r3/clZFfz5qhn84Ond7D9yPIEJRSQZLrroIu655x4AnnjiCUpLSykqKuLiiy/m3nvvBeC3v/0thw8H441r1qzh/vvvp76+HgjGLPbujWpn7YRRcYijJ3c0sLnmKJ+8ZD652bH90f71pQsA+O8ndyUimogk0Ze+9CVeeuklVqxYwS233MJdd90FwK233spTTz3F0qVLeeCBB5g5cyYAS5Ys4Stf+QqXX345K1as4LLLLqO2tjbM3wLWe7Q8naxatcpT6bAfd+c9//UstUeO88Tn386YnNjr7i0/38wDr+xn/d+/nclF+QlIKZLeKisrOeMMdb1Ga6A/LzN7yd1XDfVetRzi5NldTby09zAfv2TesAoDwCcumUdnVzff+2O4A1EiIioOcfLtdTuZXJjHe0ewoG3WpPFce+Z07n7udQ5FVlOLiIRBxSEOXtxziOd3H+Lmt80jPzd7RJ/1qbfP40RnF3eu3x2ndCIisVNxiIN7nttL8dhcPnDuzBF/1vzJhbxr2VTuemYPzSc64pBORCR2oRQHM5tgZveb2TYzqzSz1WZWYmaPmNnOyOPEMLLF6kRHF4+8dpArl05h7JiRtRp6fPSiObSc7OR3W+vi8nkiIrEKq+XwbeB37r4YWAlUArcA69x9AbAu8jzlPbG9nrb2Lq5eOfSCt2idNWMCM0vG8dCmA3H7TBGRWCS9OJhZMXAxcAeAu7e7+xHgWuCuyGV3AdclO9twPLS5lknjx7B67qS4faaZ8e6VU3m6qpFGnf0gkrFmz55NY2PjsK7p/XpBQUHcs4XRcpgDNAA/MLNXzOz7ZjYeKHf3nlUfdUB53zea2U1mtsHMNjQ0NCQx8sCOtXfyWGU9Vy6bQk6Mi96Gcs3K6XQ7PLwl3IUwIjI6hVEccoCzgf9097OANvp0IXmwMq/f6jx3v93dV7n7qrKysqSEPZ3HttVzvKOLq1dMi/tnL5pSyKLyQh7cqK4lkVShLbsTqwaocffnI8/vJygOB81sqrvXmtlUoD6EbDH59aZaygrzOHdOSUI+/90rp/KNP+xg/5HjTB9g22+RUe23t0Ddlvh+5pTl8M6vnfYSbdmdIO5eB+wzs0WRl9YArwEPAmsjr60FfpXsbLFoPdnJ49vruWr5VLKzEnPg+btXBi0SDUyLpA5t2Z1YnwbuMbMxQDXwEYJCdZ+ZfRTYC7w3pGxRefS1g5zs7ObqFfGbpdTXrEnjWTljAg9tOsDH3zYvYfcRSUtD/ISfKNqyO4HcfWNk7GCFu1/n7ofdvcnd17j7Ane/1N0PhZEtWr/efICpxfmcPTOxyzGuWTmNVw80s6uhNaH3EZH40Zbdo9TJzi7WVzVyxdIpZCWoS6nH1SumYga/2axZSyLpQlt2hyjMLbufq27ifbc/x/duWMVlS/rNuI27a/5jPXk5Wfzs4+cn/F4iqUxbdsdGW3Yn2TO7msgyEjZLqa8L5pfyyutHaDs5vP5MEZFYqTgMw7O7Glk+vZjisblJud8F80rp7HZe2J3SwzAikkFUHGJ0rL2TV14/wvnzS5N2z1WzJzImJ4v1VadfZi8yGqRrV3iyjfTPScUhRi/sPkRnt3P+vPjtpTSU/Nxszpk9kadVHGSUy8/Pp6mpSQViCO5OU1MT+fnDP244rHUOaevZXU3kZhurZiVnvKHHBfNL+Zffbaeh5SRlhXlDv0EkA1VUVFBTU0Mq7K2W6vLz86moqBj2+1UcYvTMribOmjkxbmc3ROuCeaXAdp7Z1ci1Z05P6r1FUkVubi5z5swJO8aooG6lGBw51s7WA0cj/1An17LpxRTl56hrSUSSQsUhBs9VH8Idzp+fvPGGHtlZxvnzSnm6Sv2tIpJ4Kg4xeHZXI2Nzs1lZMSGU+1+woJT9R46zt+lYKPcXkdFDxSEGz+xq4tw5JYzJCeeP7cLI9FlNaRWRRFNxiFJ9ywl21rcmdQprX7MnjWNacb7GHUQk4VQcovTSnmD3xGRtmTEQM2P1vFKe331I4w4iklAqDlHaVHOU3GxjybSiUHOcNXMCh9raqTl8PNQcIpLZVByitGnfEc6YWkReTnLXN/TVMxi+qeZIqDlEJLOpOEShu9vZsv9oaLOUels0pZAx2VlsrjkadhQRyWAqDlGobmyl9WQnK2eEXxzG5GSxZFoRG/ep5SAiiaPiEIWN+4Kf0ldWFIecJLCyopit+4/S1a1BaRFJDBWHKGzad4SCvBzmlhWEHQWAFRUTONbepXOlRSRhVByisLnmCMunF5Od4POio7VyRtCC2aSuJRFJEBWHIZzs7OK12mZWzEiNLiWAuaUFFOTlaMaSiCSMisMQKmtb6OhyzkyBmUo9srKM5dOLNWNJRBImlOJgZnvMbIuZbTSzDZHXSszsETPbGXmcGEa2vjZHfjpPhZlKva2YUUxlbTMnO7vCjiIiGSjMlsPb3f1Md18VeX4LsM7dFwDrIs9Dt3HfEUoL8phaPPzj9hJhZcUEOrqcytqWsKOISAZKpW6la4G7Il/fBVwXYpZTNu07wpkzijFLjcHoHj0tmc0adxCRBAirODjwBzN7ycxuirxW7u61ka/rgPK+bzKzm8xsg5ltSMYZss0nOqhubEuJldF9TSvOp7RgDJv2adxBROIvrDOkL3T3/WY2GXjEzLb1/qa7u5n1W+Hl7rcDtwOsWrUq4SvAttYcxR1WpNh4AwQ7tK6omKCWg4gkRCgtB3ffH3msB34BnAscNLOpAJHH+jCy9baxZzA6RVZG97WiopiqhmBrDxGReEp6cTCz8WZW2PM1cDmwFXgQWBu5bC3wq2Rn6+vVA83MKBnLhHFjwo4yoJUVE3CHV/era0lE4iuMbqVy4BeRAd4c4F53/52ZvQjcZ2YfBfYC7w0h25tsq21m8ZRwz284nYVTCgHYWd/KW+eGd0KdiGSepBcHd68GVg7wehOwJtl5BnOio4vdjW1ctXxq2FEGNa04n/Fjsqmq1x5LIhJfqTSVNaVU1bfS7bB4auq2HMyM+ZML2FmvtQ4iEl8qDoOorG0GYHGk6yZVzZ9cyM6DajmISHypOAxie10L+blZzJo0Puwop7WgvID6lpMcPd4RdhQRySAqDoPYVtfCwvLClNmmezALJgdnTGjcQUTiScVhENvqmlO+SwlgweQgY5XGHUQkjlQcBtDQcpLG1nYWpfA01h7TJ44lPzdL4w4iElcqDgPYXhf8FH5GGrQcsrOMeWUF7FS3kojEUVh7K6W0bXXBTKVFySoObY2w9edQXwldHdDVDlk5MPdtsOByGFdy2rfPn1zAhj2Hk5NVREYFFYcBVNa2MLkwj0kFeYm7iTvs/AO8dBfs/D10d8K4SZAzFrJz4WQLbLoXLBtmnQ/nfxoWXjHgRy2YXMCvNh6g9WQnBXn6TyoiI6d/SQaw/WBzYlsNJ5rh138DW++H8ZPhrR+HMz8A5UvfuKa7Gw68Att+Da8+APe+F5b+CVz5dSh8827m8yOD0rvqW1PuxDoRSU8qDn10dnWz42ArN54/OzE32P8y3P8XcGQvvP1/wYV/A9kD/GfIyoKKtwS/LvkCPP1teOpfYNdj8M5/gZXvO3XpgvI3prOqOIhIPGhAuo89TW20d3YnZhrr1p/DHZcH4wo3Pgxv+/zAhaGvnDHBtZ94BiYvhV/cDOu/derbs0rGkZttGpQWkbhRcehjW2SmUty7lXY9Dg/cDBXnwMf/CLNWx/4ZpQtg7YOw7M/g0Vvh8a+COznZWcwtLdBaBxGJG3Ur9bGttoXsrGBDu7ip3QQ//RCULoT3/wTGjqDrJzsX/vR7wcD1k1+Djja47P8wv7yArTrXQUTiRMWhj211zcwrG09eTnZ8PvDQbrj7PTB2Inzo/pEVhh5Z2XDNdyB3LDzzHSgoZ8HkK3h4Sy0nOrrIz41TdhEZtdSt1Me2upb4rYzuPAn/8wHo7oAP/RyKpsXncyEYsH7Xv8IZ18Ajt/LWrG24w64GjTuIyMipOPTSdrKTmsPHWVQepy6lp/4V6l+DP7kdyhbF5zN7M4Nrvwslcznnxc8xmcPagE9E4kLFoZfqhjaA+Iw31G6CP/4brPwALLx85J83mPwi+PMfk93Zxv8b8+/sqj2SuHuJyKih4tBLT5fMvLIRFofOdvjlp2B8KVx5WxySDWHyGdg132FV1naW7fhO4u8nIhlPxaGXqvpWsrNs5Af8rP8mHNwCV38zGIhOhuXv4YlxV7Dm8E/h4KvJuaeIZCwVh16q6luZNWkcY3JG8MfSWBWMNSx7Dyy+Kn7hovDkzE/TTEGwNUd3d1LvLSKZRcWhl6qG1pF3KT3x1WAtwpVfjU+oGJSUTeGfO94P+56HV36c9PuLSOYIpTiYWbaZvWJmv448n2Nmz5tZlZn91MzGJDtTR1c3e5vaRjYYffDVYIuMt94MBZPjFy5KFSVjub/rYo5POy9YQd3WmPQMIpIZwmo5fAao7PX868A33X0+cBj4aLIDvX7oGB1dzvyRtBwevw3yCuH8v45fsBhUTBwHGFvPujXY8vsP/zuUHCKS/pJeHMysArgK+H7kuQHvAO6PXHIXcF2yc/WsD5g33JbD/peD7bVX/9WQh/MkyvQJYwHY2V0R5Nh0L9RtDSWLiKS3mIuDmf3DCO/5LeDvgZ4R00nAEXfvjDyvAaYPcu+bzGyDmW1oaGgYYYw3O1UcyoY5U+mxr8DYEjjvE3FMFZvyonxysoyaw8fgws9CXlGw/5KISIyGLA5mdl+vXz8DPjbcm5nZ1UC9u780nPe7++3uvsrdV5WVlQ03xoB2NbQypSifwvzc2N+891nYtS74Bzk/TltvDEN2ljFtwlhqDh8PptCe9wmofAhqN4eWSUTSUzQth2Z3f2/k1/XAoyO43wXANWa2B/gfgu6kbwMTzKxnE8AKYP8I7jEsu+pbhz8Y/cy/w7hSOOcv4xtqGComjg1aDgDnfRLyiuEJtR5EJDbRFId/7vP8n4Z7M3f/grtXuPts4H3AY+7+QeBx4D2Ry9YCvxruPYaZi10NbcPrUjq8F7b/Ft5yI4wZF/dssZre03KAYAfY1Z+C7b+BAxvDDSYiaWXI4uDuuwHMrDTy/FACcvwD8DkzqyIYg7gjAfcYVF3zCVpPdg6v5bDhDrAsWPWR+AcbhoqJ46hvOcnJzq7ghfM+DvlqPYhIbGIZkL4znjd29yfc/erI19Xufq67z3f36939ZDzvNZRd9cGGezHPVOo4Di//KFgJXVyRgGSxq5gYzFg6cORE8EJ+Maz+NOz4LRx4JcRkIpJOYikOlrAUIes5XjPmNQ5b7ofjh+HcmxKQanh6isOpcQcIFuWNKYDn/zukVCKSbmIpDp6wFCGramilMD+HssK86N/kDi/8N0xeArMvTFy4GFWUBOMep8YdIJhBtfJ9weptrZoWkSio5UCwxmH+5AKC9XhR2vcC1G2Bc/8yOHQnRZQX5pGdZezvXRwgmEnV1R50g4mIDCGW4vCFhKUIWVV9W+xdSi/cHkwTXf7exIQappzsLKYW57+5Wwlg8mKYfRFsuBO6Ogd+s4hIRNTFwd0zch+Go8c6aGw9Gdtg9PEjUPkgrPxzyIvTkaJxFKx1ON7/G+feBEf3wY7fJT+UiKSVmLbPMLOnzKwo8vXHzeyzYeygGk9VkdPfYmo5VD4YdNGsfF+CUo1MxcRxAxeHRe+Counw4veSH0pE0kqseysVu3uzmb0F+EtgIpDW/9JUR4rD3FgWwG2+D0rmwbSzE5RqZComjuVgywnaO/sc+JOdE6zHqH4CGnaEkk1E0kOsxaEjss3FDcDX3f1WYGn8YyVPdWMbOVnGjJIoVzc3H4A962HFe1NqILq36RPG4g61RwdoPZx9I2SPgRe/n/RcIpI+Yi0O3wE2AVcDD0VeS71O9xjsbmhj5qRx5GZH+Uex5X7AYfn1Cc01EsG5DgzctVRQFnQvbf05dHUkOZmIpIuYioO73wWsBpa5+3Ezmw88m5BkSVLd2Mrc0hi6lLbcB9PfApPmJS7UCA24EK63Fe+FY41Q/WQSU4lIOol1QPp7wDZgh5k9T7An0suJCJYMXd3OnqZjzI12MLp+W7C2IcWmr/Y1tTif7CwbuOUAMP/SYFuNLfclN5iIpI2coS95k4uBGe7eZWbTgZXAivjHSo4DR47T3tnNnGhbDlvuA8uGZX+a2GAjlJOdxZSi/MGLQ04eLLkWtvwc2o+lxG6yIpJaYh1zeJ5g11Tcfb+7P+zuabvdZ3VjsOFeVN1K7rDlZzD3EiiYnNBc8TB94tj+q6R7W349dLQFG/KJiPQRa3H4b+BJM/s7M7vIzIoTESpZdkemsc6JZhprzYtw5PWgvz4NvOnQn4HMugAKp0UG2EVE3izW4nA38COC7qhPAs+Y2a64p0qS6sY2CvJyKCuIYsO9bb+GrBxYeGXig8XB1OJ8DracpLt7kP0SsyLdYzsfgWOJOKJDRNJZrMWhxt2/6u5fc/f3u/tSYEkigiXD7sY25paNj27DvW0PB7uvjp2Q+GBxUDI+j65up/nEaaarLr8eujuCFd8iIr3EWhw2mtlner+Q7IN54qm6oS26wejGndC0ExZdlfhQcTJpfLCrSVNb++AXTV0Jkxaoa0lE+om1OJQDHzezA2b2azP7ZzNL3dVgp3Gio4sDR48ztzSKaazbfhM8LnpnYkPFUUlPcWg9TXEwg+XvCVZ8t9YnKZmIpINYF8G9193PAOYAXwR2AG9NRLBE29PUhnuUg9HbH4YpK2DCjMQHi5NJBUFxONQ2RMNu0TsBD8YeREQiYl0E9xUIupLc/eXIiul/SEiyBKtuiHIaa2t9cLDP4vTpUgKYND4YZD9ttxIERa9wKuz8QxJSiUi6iLVbabqZfaDniZmVAY/GN1Jy7I6scRhyzGHH7wAP9iNKIxPH5wJw6HTdShB0LS24DHY9pr2WROSUWIvDzcBfmtm5ZnYO8DjwjfjHSrxdDa2UF+UxPm+IReLbHobimTBleXKCxUleTjaFeTlDtxwAFlwBJ5vh9ecSH0xE0kJU22eY2Y8I9lB6BfgUcC/QCVzn7lWJi5c4uxvbhh6Mbm+D6sfh7LUpuz336ZQUjOFQNMVh7iWQlQs7fw9zLkp0LBFJA9G2HH4IGPARgoVws4HDwIfM7D2x3NDM8s3sBTPbZGavmtmXI6/PMbPnzazKzH6a6BPmdje2DT0Yvetx6DwBi9OrS6lHyfgxNA01IA3BUaezL4AdGncQkUBUxcHdH3P3b7r7je5+NlAKfA6oIvbZSieBd7j7SuBM4EozOw/4OvBNd59PUHg+GuPnRu1QWztHjnUMPRi943eQVxxsNZGGJo0fc/qprL0tuAIat8PhPQnNJCLpIariYGarrdcyYnfvdPct7n63u38+lht6oDXyNDfyy4F3AD2rse4Crovlc2OxuzGKo0Hdg0HaeZdAdm6ioiTUpPF50XUrASy8InjUlFYRIfpupRuAl83sf8zsRjObMpKbmlm2mW0E6oFHgF3AEXfvjFxSA0wf4H03mdkGM9vQ0NAw7PvvOjWN9TRjDo07oHk/zHvHsO8TtpKCMRw+1o77IPsr9TZpXnAu9o7fJz6YiKS8aLuVPuHuZwFfAiYCPzSzZ83sNjO72MyyY7mpu3e5+5lABXAusDjK993u7qvcfVVZWVkst3yT3ZFzo3tOTBtQ1brgMY2Lw6TxY+jocppPdA59MQSthz1/DM54EJFRLdYV0tsiYw9XEnQDrQeuJzjnIWbufoRgOuxqYIKZ9cyeqgD2D+czo1Hd0MrMSePIOd250bseC/YdmjAzUTESrmcLjai7lhZcFgzA734qgalEJB3Eus7hFHc/Dsx190+7+yozi2q7UjMr67nWzMYClwGVBEWiZ+bTWuBXw802lCGnsXacCPYbSuNWA/QuDlHujTjrAsjOC1oPIjKqDbs4RMzq9fU/RvmeqcDjZrYZeBF4xN1/TbANx+fMrIrgtLk7RphtQD3nRs873WD0vueg8zjMX5OICEnTs4VGY7QzlnLyYPrZWgwnIlEvgjvD3SsH+FaWmV0EPA2URPNZ7r4ZOGuA16sJxh8Sav/hKM6NrloXLApL0ymsPUoKYuxWApi5Gp7592AB4Jgoz9YWkYwTbcvhN2b2AzPr2wH/eWAF8D0S2A0UT7tOTWM9TbfSrsdh5nnB4rA0NinWMQeAWedDdyfUbEhQKhFJB9EWh8UE22c8aWbfjmy4h7t3u/t33f2j7v5QwlLG0e6eaayDdSu1HISDW9J+vAEgPzeb8WOyo18IB1BxDmDqWhIZ5aKdytru7t8BzgD2AS+Y2f8xs6KEpkuA6sZWivJzTv1U3f+Cx4PHDCgO0LO/UgyH9Y2dAOXL4PVnEhdKRFJerFNZT7j7N4BlwHHgJTP7u4QkS5DqhjbmlBUMfm501ToYVxqcc5ABSsbnRbcza28zz4N9L0JXlOsjRCTjxHrYz2wzuxL4GDATaAFuS0SwRNnd2Ma8wQaj3YOWw9xLIGukE7lSQ0z7K/WYtRo62qBuc2JCiUjKi3Zvpc1m1gT8gmANQjGwjmBbjbQZtT3W3knt0RODjzfUV0JbQ1AcMkTJ+Ci37e5t5urg8fVn4x9IRNJCVFNZCTbB2+1RbdKTunqOBp0z2AK4nsVfGXSmwaTImQ7uPnhXWl9F02DCLNj7DKz+VGIDikhKinZAujrdCwNAdeMQM5V2PxWc+jZxdvJCJdik8WNo7+qm9WSM4wczVwczltL/P7uIDENmdKxHaXdDG2aDnBvd3Q17n86oVgMEA9IQ41oHCMYdjjVCU1oe9CciIzSqikN1YyvTiseSnzvAJrIHt8LxwzA7s4pDz5Td2GcsnR88atxBZFQaXcWhoW3wLqUMHG+AXpvvxTpjqXQBjJsEe1UcREajUVMc3D2yG+tg4w1/hJK5UFyR3GAJVnKq5RDDQjgAs2C19IFXEpBKRFLdqCkODS0naT3ZOfCeSt1dwcycDOtSgmC2EgyjWwmgfGlwIl7HiTinEpFUN2qKw65T01gHaDnUboKTR2HOxUlOlXjjxuQwNjc79m4lCIqDd0Hj9vgHE5GUNmqKQ/Wp3VgHKA494w2zL0xiouQZ1kI4gPLlwePBV+MLh/vIAAAQAklEQVQbSERS3qgpDrsb2sjLyWJa8QDnRu/+I5QuhMIpyQ+WBJMKxgyvW6lkLuTkqziIjEKjpjhUN7Yxp3Q8WVl9Vgl3dQTTNTNwvKHHsFsO2TlQtjiY5isio8roKQ4NrQN3KdVugvbWjJvC2lvJ+DE0tcY4W6lH+TK1HERGoVFRHNo7u9l3+DhzB9pTqWe8YVZmjjdAZGfWyP5KMStfGmxG2Fof/2AikrJGRXHY09RGV7czf/IAxWHvM1C6CArKkh8sSUrG53Gys5tj7V2xv7l8afCoriWRUWVUFIeq+mCmUr/i0N0VbC43+4IQUiVPz1qH4c1Y6ikO6loSGU1GRXHYebAVM5jXdwFc3WY42QyzMrw4DHd/JYDxpVAwRcVBZJRJenEwsxlm9riZvWZmr5rZZyKvl5jZI2a2M/I4MV73rGpopWLiWMaO6bPh3t7IOckZXhxO7a8U6xYaPcqXqltJZJQJo+XQCfytuy8BzgM+ZWZLgFuAde6+gOCUuVvidcOdB1uYP9C2GXueDubyF02N161S0syScXz20gXMLBlkX6mhlC+Fhu3BtF8RGRWSXhzcvdbdX4583QJUAtOBa4G7IpfdRXD63Ih1dTvVjW0sKC988ze6u+H1ZzK+1QAwqSCPz166cOAB+WiUL4Oudp3tIDKKhDrmYGazgbOA54Fyd6+NfKsOKB/g+pvMbIOZbWhoaIjqHvsOHaO9s7t/y6H+teD8hlFQHEZMg9Iio05oxcHMCoCfA5919+be34scSdpvUr673+7uq9x9VVlZdFNPT81UKu9THHrGGzJ8plJclC6ErByNO4iMIqEUBzPLJSgM97j7A5GXD5rZ1Mj3pwJxWXW1c7BprHvXB+dFT5gZj9tktpwxwVoQtRxERo0wZisZcAdQ6e7/1utbDwJrI1+vBX4Vj/tV1bdSXpRHUX7uGy+6By2HWefH4xajQ/lSFQeRUSSMlsMFwIeBd5jZxsivdwFfAy4zs53ApZHnI1ZV38KCyX0Goxt3BltCqEspeuVLoXl/ME4jIhkvJ9k3dPf1gA3y7TVxvhdV9a1cv2rGm7+xd33wqMHo6JXMCR6P7IOxcVuCIiIpKqNXSNcePUFbe1f/8YY9TwerfkvmhhMsHRVNDx6bD4SbQ0SSIqOLw4B7KrnD3qeDU99ssAaM9FMYWSjYvD/cHCKSFBldHHpmKi3oXRwOVUNLbcYeCZowBeVgWcGfnYhkvIwuDlX1LUwcl8ukgrw3Xszw86ITJjsn6IpTt5LIqJDhxaG1/0ylPU8HPwVPmh9OqHRWNE3dSiKjRMYWB3dnZ33rm1dGu8Oe9RpvGK6iaWo5iIwSGVscmtraOXKs4817Kh2qhpYDmsI6XEXToFljDiKjQcYWh50HI4PRvVsOe58OHmdfFEKiDFA0Ddpb4ETz0NeKSFrL2OJQVd8C9JnGumc9jJ8MpQtCSpXmtNZBZNTI2OJQWddCUX4OU4rygxdOjTdcoPGG4SqaFjxqUFok42Vscdhe18LiqUVYTyE4vDv4R01TWIfvVHFQy0Ek02VkcXB3dtS1sHhKr2msezTeMGI9q6S1EE4k42Vkcdh/5DgtJztZ9KbisB7GlwUH18jw5OTBuFJ1K4mMAhlZHLbVBoPRp1oOPeMNszTeMGJa6yAyKmRkcdh+MCgOC8sjxeFQNTTXaLwhHoqmqziIjAIZWRy21bVQMXEshT2nv1U/ETzOfXtomTKGWg4io0JGFoftdc1vHoyufgKKKmDSvNAyZYyiqXD8EHQcDzuJiCRQxhWHk51d7Gpoe2Mwursr2Il17iUab4gHLYQTGRUyrjjsqm+jq9tZNKUoeKFuc3Du8dxLwoyVObTWQWRUyLjisP1gsO/PGT0th1PjDW8LJ1CmUctBZFTIuOKwra6FMdlZzC4dH7xQ/QRMXgoFk0PNlTFOLYRTcRDJZJlXHGpbmDe5gNzsrGDQdO+z6lKKp7wCyCtWy0Ekw2Vccdjee9uMfc9D10kVh3jTdFaRjJf04mBmd5pZvZlt7fVaiZk9YmY7I48Th/PZR491UNd84o2ZStVPQFYOzDo/LtklQseFimS8MFoOPwSu7PPaLcA6d18ArIs8j9m2umAwenHv4lBxbtAVIvGjE+FEMl7Si4O7PwUc6vPytcBdka/vAq4bzmdvq+vZU6kIjh2CAxvVpZQIRdOg9SB0dYSdREQSJFXGHMrdvedH0TqgfKCLzOwmM9tgZhsaGhr6fX9bXQvFY3MpL8oLFr7hmsKaCEXTAIeWurCTiEiCpEpxOMXdHfBBvne7u69y91VlZWX9vr+trplFUwqDA352PhLMqpm+KtGRRx+tdRDJeKlSHA6a2VSAyGN9rB/Q1e1sr2thydSiYIvuqnUw7xLIzol3VtFxoSIZL1WKw4PA2sjXa4FfxfoBe5vaONbeFRSH+teCRVrzL4trSInoKQ46EU4kY4UxlfUnwLPAIjOrMbOPAl8DLjOzncClkecxqYwc8LNkWlHQpQQwf02cUsub5E+Ai/4Opr8l7CQikiBJ73Nx9/cP8q0R/Uv+Wu1RsrOM+ZML4JFHoXzZGz/hSnyZwZr/HXYKEUmgVOlWGrHK2hbmlxWQ39UGrz8L8y8NO5KISNrKoOLQzBlTC2H3k9DdCQs03iAiMlwZURwOt7VTe/QEZ0yNjDfkFcGMt4YdS0QkbWVEcaisDbbNWDK1EKoeDRa+ZeeGnEpEJH1lRHF4LVIcluUeCObeawqriMiIZExxmFyYx8QDTwYvaDBaRGREMqM4HGh+Y7xh8lIonh52JBGRtJb2xaG9s5tdDa2cNdmCKawLLw87kohI2kv74rCzvoWOLudC2xxMYV3Y96gIERGJVdoXh55tMxYefRrGToSKc0JOJCKS/jKgODQzLhcK9z0OCy6HrOywI4mIpL20Lw6vHWjm6pID2PFDsPCKsOOIiGSEtC4O7k5lXTPvzNsElg3ztAuriEg8pHVxqD16giPHOlh57HmYdT6MnRB2JBGRjJDWxaG0II8HPzyLktYdwXiDiIjERVoXhzE5Waw49lzwRFNYRUTiJq2LAwA7fg8T50DpgrCTiIhkjPQuDu3HYPdTQavBLOw0IiIZI72Lw4mjsPhqOOPdYScREckoST9DOq6KpsJ77gg7hYhIxknvloOIiCSEioOIiPSTUsXBzK40s+1mVmVmt4SdR0RktEqZ4mBm2cB3gXcCS4D3m9mScFOJiIxOKVMcgHOBKnevdvd24H+Aa0POJCIyKqVScZgO7Ov1vCby2ilmdpOZbTCzDQ0NDUkNJyIymqRScRiSu9/u7qvcfVVZWVnYcUREMlYqFYf9wIxezysir4mISJKZu4edAQAzywF2AGsIisKLwAfc/dVBrm8BticvYVyUAo1hh4hRumVOt7ygzMmQbnkhcZlnufuQXS8ps0La3TvN7K+A3wPZwJ2DFYaI7e6+Kjnp4sPMNihzYqVbXlDmZEi3vBB+5pQpDgDu/jDwcNg5RERGu1QacxARkRSRzsXh9rADDIMyJ1665QVlToZ0ywshZ06ZAWkREUkd6dxyEBGRBFFxEBGRftKyOKTb7q1mNsPMHjez18zsVTP7TNiZomFm2Wb2ipn9Ouws0TCzCWZ2v5ltM7NKM1sddqbTMbO/ifx92GpmPzGz/LAz9WVmd5pZvZlt7fVaiZk9YmY7I48Tw8zY1yCZ/zXy92Kzmf3CzCaEmbGvgTL3+t7fmpmbWWkyM6VdcUjT3Vs7gb919yXAecCn0iAzwGeAyrBDxODbwO/cfTGwkhTObmbTgb8GVrn7MoK1Pe8LN9WAfghc2ee1W4B17r4AWBd5nkp+SP/MjwDL3H0FwWLbLyQ71BB+SP/MmNkM4HLg9WQHSrviQBru3urute7+cuTrFoJ/tKaf/l3hMrMK4Crg+2FniYaZFQMXA3cAuHu7ux8JN9WQcoCxkd0BxgEHQs7Tj7s/BRzq8/K1wF2Rr+8CrktqqCEMlNnd/+DunZGnzxFsz5MyBvlzBvgm8PdA0mcOpWNxGHL31lRmZrOBs4Dnw00ypG8R/KXsDjtIlOYADcAPIl1h3zez8WGHGoy77we+QfATYS1w1N3/EG6qqJW7e23k6zqgPMwww/AXwG/DDjEUM7sW2O/um8K4fzoWh7RlZgXAz4HPuntz2HkGY2ZXA/Xu/lLYWWKQA5wN/Ke7nwW0kXrdHadE+umvJShq04DxZvahcFPFzoO58GkzH97M/omgm/eesLOcjpmNA/4R+GJYGdKxOKTl7q1mlktQGO5x9wfCzjOEC4BrzGwPQbfdO8zs7nAjDakGqHH3nhbZ/QTFIlVdCux29wZ37wAeAM4POVO0DprZVIDIY33IeaJiZjcCVwMf9NRf4DWP4AeHTZH/DyuAl81sSrICpGNxeBFYYGZzzGwMwSDegyFnOi0zM4K+8Ep3/7ew8wzF3b/g7hXuPpvgz/cxd0/pn2rdvQ7YZ2aLIi+tAV4LMdJQXgfOM7Nxkb8fa0jhAfQ+HgTWRr5eC/wqxCxRMbMrCbpJr3H3Y2HnGYq7b3H3ye4+O/L/YQ1wduTveVKkXXGIDCr17N5aCdw3xO6tqeAC4MMEP4FvjPx6V9ihMtCngXvMbDNwJnBbyHkGFWnh3A+8DGwh+H8x5bZ4MLOfAM8Ci8ysxsw+CnwNuMzMdhK0gL4WZsa+Bsn8H0Ah8Ejk/7//CjVkH4NkDjdT6reuREQk2dKu5SAiIomn4iAiIv2oOIiISD8qDiIi0o+Kg4iI9KPiICIi/ag4iIhIPyoOIkMws5vNrM7MNpnZLjO7odf3/svMLoh8vcbMfhxeUpH4UXEQGdpy4EvuvhJ4P9B7C5TzCLaAhuAMiVeSnE0kIVQcRIa2Atge+Xo30A5gZmcAO9y9K/K9lcArZpZnZj80s9ss8BMz+6mZvWBme83sqsj7p5nZzyNbjG8zs3OT/RsTGYyKg8jQlgPbIxvk/RXwT5HX3wn8rtd1Kwh2KP098Ki7/2Nk98+VQLW7nwt8ELg1csDPb4EfRLYYP5v02XhPRgHtrSRyGpFjGncDWwkOldoMXOrubma/Bz7i7gciW7I3AnuBm9392cj78wkOp5rh7ifMrITgoKcvAO9x91Q8GlRELQeRISwHnnL3M4GFwGJgdeQwlgnu3nO05xkE28l3Al293r8M2OnuJyLPzwY2Eewa+xwiKUrFQeT0VhAZZHb3w8C9BGdrvx14vNd1K4FnCM6/+IGZlfd6faaZ5UeOLf0ywbnAdcDSnjebWVmCfx8iMVFxEDm95bx5BtJDwLvoP96wEtjq7juAfwDui3Q1rSQ45e15gpbFf7r708APgXIze9XMNgKrE/0bEYmFxhxEhsHMXgbeGjni83TXPQnc5O7bT3edSKpRcRBJIDOrAWa6e3fYWURioeIgIiL9aMxBRET6UXEQEZF+VBxERKQfFQcREelHxUFERPpRcRARkX5UHEREpJ//D8uI7MuUFqa3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(R,np.sqrt(vc2_ISM(R,R_d=R_d_I,Sigma_g_ISM=Sigma_g_ISM_I,R_m=R_m)),label='model I')\n",
    "plt.plot(R,np.sqrt(vc2_ISM(R,R_d=R_d_II,Sigma_g_ISM=Sigma_g_ISM_II,R_m=R_m)),label='model II')\n",
    "plt.xlabel(r'$R/kpc$')\n",
    "plt.ylabel(r'$V_{c}/km s^{-1}$')\n",
    "plt.title('ISM')\n",
    "plt.legend(loc='upper right', prop={'size':10})\n",
    "plt.xlim(0.,15.)\n",
    "#plt.savefig('Model_II.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Y finalmente comparando con el modelo aproximado"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def M_g_I(u,R,R_d=2.,R_m=4.):\n",
    "    return u*np.exp(-((R*u)/(2.*R_d))-(R_m/(R*u)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def V_c_g_I(l):\n",
    "    Integral_M_g = []\n",
    "    for i in range (0,len(l)):\n",
    "        result = integrate.quad(M_g_I,0,1,args=(l[i]))[0]\n",
    "        Integral_M_g.append(result)\n",
    "    Integral = np.array(Integral_M_g)\n",
    "    V = l*Integral\n",
    "    return V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vc2_ISM_ap(R,R_d,Sigma_g_ISM,R_m=R_m):\n",
    "    def Int_ISM(h,s,R):\n",
    "        t1 = 1./np.sqrt(1.-s**2.)\n",
    "        t2 = 1./np.sqrt(h**2. - 1.)\n",
    "        t3 = np.exp(-(R*s*h)/(2.*R_d) - (R_m)/(R*s*h))#t4 = (2.*R_m*R_d)/(R**2.*h**2.) - s**2.\n",
    "        res = R**2.*t1*t2*t3#*t4\n",
    "        return res\n",
    "    Integral_ISM = []\n",
    "    for i in range(0,len(R)):\n",
    "        integrando = lambda h,s: Int_ISM(h,s,R[i])\n",
    "        result = integrate.dblquad(integrando,0.1,1.,lambda s: 1., lambda s: np.inf)[0]\n",
    "        Integral_ISM.append(result)\n",
    "    return 4.*np.pi*G*Sigma_g_ISM*np.array(Integral_ISM)/(2.*R_d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(R,vc2_ISM_ap(R,R_d=R_d_I,Sigma_g_ISM=Sigma_g_ISM_I))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vc_I = np.sqrt(G*2.*np.pi*Sigma_g_ISM_I*V_c_g_I(R))\n",
    "plt.plot(R,Vc_I,label='aprox')\n",
    "plt.plot(R,np.sqrt(vc2_ISM_ap(R,R_d=R_d_I,Sigma_g_ISM=Sigma_g_ISM_I)),label='exact without t4')\n",
    "plt.xlabel(r'$R/kpc$')\n",
    "plt.ylabel(r'$V_{c}/km s^{-1}$')\n",
    "plt.title('ISM model I')\n",
    "plt.legend(loc='upper right', prop={'size':10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vc2_v4_ISM_II(h,s,R,R_d=R_d_II,R_m=R_m):\n",
    "    t1 = s**2./np.sqrt(1.-s**2.)\n",
    "    t2 = 1./np.sqrt(h**2. - 1.)\n",
    "    t3 = np.exp(-(R*s*h)/(2.*R_d) - (R_m)/(R*s*h))\n",
    "    t4 = (2.*R_m*R_d)/(R**2.*s**2.*h**2.) - 1.\n",
    "    res = R**2.*t1*t2*t3*t4\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Integ = []\n",
    "for i in range (0,len(R)):\n",
    "    integrando = lambda h,s: vc2_v4_ISM_II(h,s,R[i])\n",
    "    result = integrate.dblquad(integrando,0.1,1.,lambda s: 1., lambda s: np.inf)[0]\n",
    "    Integ.append(result)\n",
    "print Integ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integ_array_II = np.array(Integ)\n",
    "C_integ_II = -2.*np.pi*G*Sigma_g_ISM_II/R_d_II\n",
    "plt.plot(R,C_integ_II*integ_array_II)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(R,np.sqrt(C_integ_II*integ_array_II))\n",
    "plt.xlim(0,15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def M_g_II(u,R,R_d=3.2,R_m=4.):\n",
    "    return u*np.exp(-((R*u)/(2.*R_d))-(R_m/(R*u)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def V_c_g_II(l):\n",
    "    Integral_M_g = []\n",
    "    for i in range (0,len(l)):\n",
    "        result = integrate.quad(M_g_II,0,1,args=(l[i]))[0]\n",
    "        Integral_M_g.append(result)\n",
    "    Integral = np.array(Integral_M_g)\n",
    "    V = l*Integral\n",
    "    return V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vc_II = np.sqrt(G*2.*np.pi*Sigma_g_ISM_II*V_c_g_II(R))\n",
    "plt.plot(R,Vc_II)\n",
    "plt.plot(R,np.sqrt(C_integ_II*integ_array_II))"
   ]
  },
  {
   "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
}