PanSTARRS_WISE_Q0.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  PanSTARRS - WISE cross-match: Compute the $Q_0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.table import Table\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord, search_around_sky"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mltier1 import generate_random_catalogue, Field, Q_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "field = Field(170.0, 190.0, 45.5, 56.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data\n",
    "\n",
    "### Load the catalogues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "panstarrs_full = Table.read(\"panstarrs_u2.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wise_full = Table.read(\"wise_u2.fits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Restrict the study to the well sampled area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "panstarrs = field.filter_catalogue(\n",
    "    panstarrs_full, \n",
    "    colnames=(\"raMean\", \"decMean\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Free memory\n",
    "del panstarrs_full"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wise = field.filter_catalogue(\n",
    "    wise_full, \n",
    "    colnames=(\"raWise\", \"decWise\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Free memory\n",
    "del wise_full"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coordinates\n",
    "\n",
    "As we will use the coordinates to make a cross-match we to load them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords_panstarrs = SkyCoord(panstarrs['raMean'], panstarrs['decMean'], unit=(u.deg, u.deg), frame='icrs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "coords_wise = SkyCoord(wise['raWise'], wise['decWise'], unit=(u.deg, u.deg), frame='icrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Cross-match and random cross-match\n",
    "\n",
    "We create an instance of Q_0 using as the input the two catalogues and the area. It will compute the $Q_0$ for different radius.\n",
    "\n",
    "The following function is not used but shows the internal code use for the computing of the $Q_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example function (not used, we use a class that contains this code)\n",
    "def q_0_r(coords_wise, coords_panstarrs, field, radius=5):\n",
    "    \"\"\"Compute the Q_0 for a given radius\"\"\"\n",
    "    random_wise = field.random_catalogue(len(coords_wise))\n",
    "    \n",
    "    idx_random_wise, idx_panstarrs, d2d, d3d = search_around_sky(\n",
    "        random_wise, coords_panstarrs, radius*u.arcsec)\n",
    "    \n",
    "    nomatch_random = len(coords_wise) - len(np.unique(idx_random_wise))\n",
    "    \n",
    "    idx_wise, idx_panstarrs, d2d, d3d = search_around_sky(\n",
    "        coords_wise, coords_panstarrs, radius*u.arcsec)\n",
    "    \n",
    "    nomatch_wise = len(coords_wise) - len(np.unique(idx_wise))\n",
    "                                          \n",
    "    return (1. - float(nomatch_wise)/float(nomatch_random))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_0_comp = Q_0(coords_wise, coords_panstarrs, field)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.629392983348576"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_0_comp(radius=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $Q_0$ dependence on the radius\n",
    "\n",
    "We will iterate 10 times for each radius. However, the error is so small that a direct computation can be accurate to 4 significant figures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_iter = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The radius tested ranges from 1 to 25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rads = list(range(1,26))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1 0.51223 +/- 0.00002 [0.51220 0.51226]\n",
      " 2 0.60789 +/- 0.00004 [0.60783 0.60795]\n",
      " 3 0.63058 +/- 0.00004 [0.63050 0.63063]\n",
      " 4 0.63359 +/- 0.00006 [0.63350 0.63368]\n",
      " 5 0.62944 +/- 0.00012 [0.62920 0.62962]\n",
      " 6 0.62105 +/- 0.00012 [0.62086 0.62123]\n",
      " 7 0.61010 +/- 0.00011 [0.60994 0.61036]\n",
      " 8 0.59739 +/- 0.00014 [0.59709 0.59756]\n",
      " 9 0.58558 +/- 0.00027 [0.58532 0.58608]\n",
      "10 0.57676 +/- 0.00011 [0.57657 0.57691]\n",
      "11 0.57078 +/- 0.00027 [0.57038 0.57132]\n",
      "12 0.56704 +/- 0.00017 [0.56661 0.56728]\n",
      "13 0.56498 +/- 0.00046 [0.56393 0.56564]\n",
      "14 0.56440 +/- 0.00053 [0.56357 0.56518]\n",
      "15 0.56420 +/- 0.00052 [0.56345 0.56501]\n",
      "16 0.56449 +/- 0.00052 [0.56374 0.56525]\n",
      "17 0.56498 +/- 0.00063 [0.56407 0.56621]\n",
      "18 0.56534 +/- 0.00058 [0.56443 0.56628]\n",
      "19 0.56610 +/- 0.00041 [0.56539 0.56684]\n",
      "20 0.56615 +/- 0.00055 [0.56536 0.56738]\n",
      "21 0.56695 +/- 0.00071 [0.56604 0.56807]\n",
      "22 0.56742 +/- 0.00069 [0.56615 0.56908]\n",
      "23 0.56797 +/- 0.00143 [0.56518 0.56984]\n",
      "24 0.56760 +/- 0.00140 [0.56549 0.57007]\n",
      "25 0.56979 +/- 0.00074 [0.56866 0.57078]\n"
     ]
    }
   ],
   "source": [
    "q_0_rad = []\n",
    "for radius in rads:\n",
    "    q_0_rad_aux = []\n",
    "    for i in range(n_iter):\n",
    "        out = q_0_comp(radius=radius)\n",
    "        q_0_rad_aux.append(out)\n",
    "    q_0_rad.append(np.mean(q_0_rad_aux))\n",
    "    print(\"{:2d} {:7.5f} +/- {:7.5f} [{:7.5f} {:7.5f}]\".format(radius, \n",
    "            np.mean(q_0_rad_aux), np.std(q_0_rad_aux), \n",
    "            np.min(q_0_rad_aux), np.max(q_0_rad_aux)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVEAAAFBCAYAAADOoJzZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGuhJREFUeJzt3X+UX3V95/Hne2YyCfkdkiEQkiA/IhBBEGaRrSj4A0v0\nKLZVj2hta7vLsiu23dPuym63ate2a+vaaldXTrZy0B7U2koVLJQWFGkrcEgQSAIkRiS/zI8JJCE/\nyI+Zee8f3zvhy5hfM5+Z78x35vk453vm3s+9c+/75mZe38+993vvNzITSdLgtIx0AZLUzAxRSSpg\niEpSAUNUkgoYopJUwBCVpAINC9GIuCUitkXEyqNMj4j4i4hYGxFPRMQljapNkgarkT3RW4FrjjF9\nCbCoel0PfLEBNUlSkYaFaGY+ADx/jFmuBb6SNQ8BMyPitMZUJ0mDM5rOiZ4ObKgb31i1SdKo1TbS\nBQxGRFxP7ZCfKVOmXHreeeeNcEWSxprly5dvz8yO4803mkJ0E7Cgbnx+1fYzMnMpsBSgs7Mzly1b\nNvzVSRpXImLdicw3mg7n7wB+pbpKfzmwKzM3j3RRknQsDeuJRsTXgKuAORGxEfg4MAEgM28G7gLe\nBqwF9gEfalRtkjRYDQvRzLzuONMT+HCDypGkITGaDuclqekYopJUwBCVpAKGqCQVMEQlqYAhKkkF\nDFFJKmCISlIBQ1SSChiiklTAEJWkAoaoJBUwRCWpgCEqSQUMUUkqYIhKUgFDVJIKGKKSVMAQlaQC\nhqgkFTBEJamAISpJBQxRSSpgiEpSAUNUkgoYopJUwBCVpAKGqCQVMEQlqYAhKkkFDFFJKmCISlIB\nQ1SSChiiklTAEJWkAoaoJBUwRCWpgCEqSQUMUUkqYIhKUgFDVJIKGKKSVMAQlaQChqgkFTBEJamA\nISpJBRoaohFxTUSsjoi1EXHTEabPiIg7I+LxiFgVER9qZH2SNFANC9GIaAW+ACwBFgPXRcTifrN9\nGHgyMy8CrgI+ExHtjapRkgaqkT3Ry4C1mflMZh4Evg5c22+eBKZFRABTgeeB7gbWKEkD0sgQPR3Y\nUDe+sWqr93ngfOCnwArgtzKzt/+CIuL6iFgWEcu6urqGq15JOq7RdmHp54HHgHnAxcDnI2J6/5ky\nc2lmdmZmZ0dHR6NrlKTDGhmim4AFdePzq7Z6HwJuz5q1wE+A8xpUnyQNWCND9BFgUUScWV0seh9w\nR7951gNvBoiIucC5wDMNrFGSBqStUSvKzO6IuBG4B2gFbsnMVRFxQzX9ZuCTwK0RsQII4KOZub1R\nNUrSQDUsRAEy8y7grn5tN9cN/xR4ayNrkqQSo+3CkiQ1FUNUkgoYopJUwBCVpAKGqCQVMEQlqYAh\nKkkFDFFJKmCISlIBQ1SSChiiklTAEJWkAoaoJBUwRCWpgCEqSQUMUUkqYIhKUgFDVJIKGKKSVMAQ\nlaQChqgkFTBEJamAISpJBQxRSSpgiEpSAUNUkgoYopJUwBCVpAKGqCQVMEQlqYAhKkkFDFFJKmCI\nSlIBQ1SSChiiklTAEJWkAoaoJBUwRCWpQNtIFzBeHOrpZce+g+w70MOcaROZOtF/emks8C95CCxf\n9zyrt+xhx76DPL/3IDv2HqwN7ztUG957kN0Hul/2O9MmtjF3xiROmzGJudMncer0SZw6o+7njEnM\nntJORIzQVkk6EYZogWe37+UP//5J7n1q2+G2ye2tzJrczslT2pk1pZ0zZ09mZt34SRNa6dp9gK0v\n7GfLrv1seWE/P9q6na49B+jpzZctf8ZJE7howUwuXjCT11Q/Z01pb/RmSjoGQ3QQ9h7o5vPfW8uX\n/vknTGgNblpyHu+8aB4nT2ln0oTWQS2zpzfZvucAW3btZ/Ou/WzZ9SKrt+7mh+t38vnv/oi+fH3F\n7MlcXAXqxQtnsfi06bS3eWpbGimG6ABkJt96bBOfuvtptr5wgF+85HRuuuY8Tpk+qXjZrS3B3Om1\nQ/uLFrx82t4D3TyxcRePbdjJYxt28IMfP8e3HvspAO1tLVwwbzpXnDOHK8/t4KL5M2lrNVSlRonM\nPP5co1hnZ2cuW7Zs2NfzxMadfOKOVTy6ficXzZ/Bx9/5Ki5ZOGvY13skmcnmXfurUN3JI88+z+Mb\ndtKbMH1SG69f1MGV53Zw5Ss7mDsEAS+NRxGxPDM7jzefPdHj2L7nAJ/+h9V8Y/kGZk+ZyJ+++9W8\n+5L5tLSM3AWfiGDezJOYN/Mk3nbhaQDs2neIf1m7nftXb+P7a7r4+xWbATjv1GmHA7XzjJM99JeG\nWEN7ohFxDfA5oBX4y8z81BHmuQr4LDAB2J6ZVx5rmcPVEz3U08uXf/Asn7v3R7x4qIdfv+JMPvKm\nc5g2acKQr2uoZSZPb9nN99d08f3VXSxb9zyHepIp7a1csWgOSy44jTedfwrTm2BbpJFyoj3RhoVo\nRLQCa4CrgY3AI8B1mflk3TwzgR8A12Tm+og4JTO3HXGBleEK0c/842r+z3fXcuUrO/jYOxZzdsfU\nIV9Ho+w50M2DP36O+1dv476ntrHlhf20t7bw+kVzuOaCU7l68VxmTvaqv1RvNB7OXwaszcxnACLi\n68C1wJN187wfuD0z1wMcL0CHy4HuHr768HquXjyXpR+8tOk/qzl1YhtXL57L1Yvn8slrkx9u2Mnd\nKzZz98ot3Pf0Ntpagp87Zw5LLjiVty6ey+ypE0e6ZKlpNDJETwc21I1vBF7bb55XAhMi4n5gGvC5\nzPxKY8p7yT2rtvLc3oN88PIzmj5A+2tpCS49YxaXnjGL33v7+azYtIu7Vmzh7pWb+W+3r+D3/m4F\nl581myUXnMrVi0/l1BlemJKOZbRdWGoDLgXeDJwEPBgRD2XmmvqZIuJ64HqAhQsXDnkRtz20joUn\nT+aKc+YM+bJHk4jg1fNn8ur5M/noNefy1Obd3L1yM3et2Mzvf3sVv//tVVx4+gzecv5c3rL4FBaf\nNn3MvalIpRoZopuA+k9Azq/a6m0EnsvMvcDeiHgAuIjaudTDMnMpsBRq50SHssi123bz8E+e56Yl\n543oFfhGiwgWz5vO4nnT+Z23nsuPtu7mn57ayn1PbeOz963hz+9dw7wZk3jz+XN5y+K5XH7WyUxs\nG9yNBdJY0sgQfQRYFBFnUgvP91E7B1rv28DnI6INaKd2uP/nDayR2x5ez4TW4D2Xzm/kakedRXOn\nsWjuNP7TVeewfc8Bvvv0Nu59cit/u3wjf/XQOqa0t/KGV3bwlvPnctW5HZ5H1bjVsBDNzO6IuBG4\nh9pHnG7JzFURcUM1/ebMfCoi/gF4Auil9jGolY2q8cWDPXxz+UaWXHCaoVBnztSJvLdzAe/tXMD+\nQz08+OPnql7qVu5euQWAszumHD7XeukZszhrztRx1ZPX+OUdS3W+sWwD//Vvn+Cvr7+c1541e0iW\nOZZlJis3vcA/r+3i0XU7WL5uBzv2HQJqD0+5ZOFMLj1jFpecMYuLF8xkcvtoOwUvHd1o/IjTqHfb\nw+s555SpXHbmySNdSlOICC6cP4ML588AaqH6k+17Wb5uB4+ur4Xq91Z3AbVnA5x/2jTOP3U6Z58y\nlbPmTOGsjqmcMXsyE7zXX03MEK2s3LSLxzfs5OPvWOwV6EGKCM7qmMpZHVN5T2ftGuKufYf44YYd\ntZ7q+h3cv6aLv1m+8fDvtLUEC0+ezFkdUzi7Y+rhn6+YM4WTJ7d7SkCjniFaue3h9Uya0MIvvmZ8\nX1AaajMmT+Cqc0/hqnNPOdz2wv5DPNO1l2e69vDjrj3V8F4eWLOdgz29h+drawk6pk2kY9pETql+\ndkybdHi8r23W5HYmt7f65tcgmUl3b9Ldk3T39tLTm7VXJr29VD/r22rz18/X0/uzr755MpOWlqAl\ngpbg8HBrv/GWgBcP9bBnfze793ez+0A3u/cfYs/+bvYc+Nm2b9/4umE5pWSIArv3H+Lbj23iHa+e\nx4zJ3k8+3KZPmnD4maj1enqTTTte5Mdde1j33F627T7Att0H6Np9gE07a0+tem7vQY50Gr8lYEp7\nG1MntTF1YhtTJrYxbVJbv7ZW2lpaaGsJWltrf5StLVGNV+1VW2tL0JfJLVEbDmp/uLX2vuEgM+lN\nDv/szaQ3k0xIasHS19ZThUxPTy89SS1s+oVK/byZfW0vLaO3bz114XOop/dwqB2q/9nTW02vhVOf\n+n/D5OX/oJnQ3Zsc7O6lu7e23IP9lj+aTW5vZVq1z6dNmsC0SW2cOn0S3b3DU7chCnzrsZ+y72AP\nH7j8jJEuZVxrbQkWzp7MwtmTjzpPd08vz+09SNfuA2zbvZ9tLxxg14uH2HOg1vvo64X0vba+sP9w\n296DPT/z7QGjUQRVrytoaaHqdb3UC2uNIKrx1pagrTWY0NJCW2vQ1tLChNagrXpTmNzedri9/6nn\nIF62znoTWltettwJrS8td0JrCxNaXlpH35tOy+E3oVrNrXXTWqM2vX7+vva21qC1peXwG1gEh9+E\n6t9UeuvHqzemk9pbq7BsY9rECbU3ygafYx/3IZqZ3PbQOl41bzoXVRdINHq1tbYcfng1DHx/9fXs\nuut6cb0vG++ltxe6e3tJ+npsfT3Nl3qWSdXTzFoAvTzwaj3U4KUAjLoAbGmhFmrVcH2g1ELE0xLN\nZNyH6KPrd/L0lt388S9c6H/ecSCqno83W2mojPvPltz28DqmTmzjnRfPG+lSJDWhcR2iO/cd5DtP\nbOZdr5nn98BLGpRxHaJ/u3wjB7t7ef9lXlCSNDjjNkQzk68+vJ5LFs5k8bzpI12OpCY1bkP0wWee\n45nte/nAa+2FShq8cRuitz28npmTJ/D2V5820qVIamLjMkS7dh/gnpVbePcl85k0wc+6SBq8cRmi\n31i2ge7e5LrXDv1Xi0gaX8ZdiPb01i4o/dzZs5v6a5AljQ7jLkQfWNPFpp0vekFJ0pAYdyF628Pr\nmDN1IlcvnjvSpUgaAwZ8m05EvB94J9ADBHBnZn5tqAsbDj/d+SLffXob//Gqs2lvG3fvH5KGwWDu\ndbwyM9/XNxIRXwCaIkRPmTaRL/7ypVx4uk9rkjQ0BhOiEyPi7cAGat8df9LQljR82lpb+PlXnTrS\nZUgaQ04oRCNiHvBmYDLw18Ac4EJgI3DjsFUnSaPccUM0It4KfBm4HzgA3EAtTD+UmT8Y1uokaZQ7\nkZ7oHwKvz8y1fQ0R8W+BpRHxG8C+zFw5XAVK0mh2IiHaXh+gAJn5YET8EvAdar3TC4ejOEka7U7k\ncz77I6Kjf2NmrqH2Mac3D3lVktQkTiREPw18q7q4dFhEzAEOZOa2YalMkprAcQ/nM/ObETEReDAi\nlgOPA+3Ae6mdL5WkceuEbtvJzK8C51M7BzoDOAS8PzO/PIy1SdKod8Ifts/MfcAtw1iLJDUdbyCX\npAKGqCQVMEQlqYAhKkkFDFFJKmCISlIBQ1SSChiiklTAEJWkAoaoJBUwRCWpgCEqSQUMUUkqYIhK\nUgFDVJIKNDREI+KaiFgdEWsj4qZjzPdvIqI7It7dyPokaaAaFqIR0Qp8AVgCLAaui4jFR5nvT4B/\nbFRtkjRYjeyJXgaszcxnMvMg8HXg2iPM9xHgm4BfgCdp1GtkiJ4ObKgb31i1HRYRpwO/AHyxgXVJ\n0qCNtgtLnwU+mpm9x5opIq6PiGURsayrq6tBpUnSzzrhL6obApuABXXj86u2ep3A1yMCYA7wtojo\nzsxv1c+UmUuBpQCdnZ05bBVL0nE0MkQfARZFxJnUwvN9wPvrZ8jMM/uGI+JW4Dv9A1SSRpOGhWhm\ndkfEjcA9QCtwS2auiogbquk3N6oWSRoqjeyJkpl3AXf1aztieGbmrzWiJkkqMdouLElSUzFEJamA\nISpJBQxRSSpgiEpSAUNUkgoYopJUwBCVpAKGqCQVMEQlqYAhKkkFDFFJKmCISlIBQ1SSChiiklTA\nEJWkAoaoJBUwRCWpgCEqSQUMUUkqYIhKUgFDVJIKGKKSVMAQlaQChqgkFTBEJamAISpJBQxRSSpg\niEpSAUNUkgoYopJUwBCVpAKGqCQVMEQlqYAhKkkFDFFJKmCISlIBQ1SSChiiklTAEJWkAoaoJBUw\nRCWpgCEqSQUMUUkqYIhKUgFDVJIKNDREI+KaiFgdEWsj4qYjTP9ARDwRESsi4gcRcVEj65OkgWpY\niEZEK/AFYAmwGLguIhb3m+0nwJWZeSHwSWBpo+qTpMFoZE/0MmBtZj6TmQeBrwPX1s+QmT/IzB3V\n6EPA/AbWJ0kD1sgQPR3YUDe+sWo7mt8A7j7ShIi4PiKWRcSyrq6uISxRkgZmVF5Yiog3UgvRjx5p\nemYuzczOzOzs6OhobHGSVKetgevaBCyoG59ftb1MRLwa+EtgSWY+16DaJGlQGtkTfQRYFBFnRkQ7\n8D7gjvoZImIhcDvwwcxc08DaJGlQGtYTzczuiLgRuAdoBW7JzFURcUM1/WbgY8Bs4P9GBEB3ZnY2\nqkZJGqjIzJGuoUhnZ2cuW7ZspMuQNMZExPIT6cSNygtLktQsDFFJKmCISlIBQ1SSChiiklTAEJWk\nAoaoJBUwRCWpgCEqSQUMUUkqYIhKUgFDVJIKGKKSVMAQlaQChqgkFTBEJamAISpJBQxRSSpgiEpS\nAUNUkgoYopJUwBCVpAKGqCQVMEQlqYAhKkkFDFFJKmCISlIBQ1SSChiiklTAEJWkAoaoJBUwRCWp\ngCEqSQUMUUkqYIhKUgFDVJIKGKKSVMAQlaQChqgkFTBEJamAISpJBQxRSSpgiEpSAUNUkgoYopJU\noKEhGhHXRMTqiFgbETcdYXpExF9U05+IiEsaWZ8kDVTDQjQiWoEvAEuAxcB1EbG432xLgEXV63rg\ni42qT5IGo5E90cuAtZn5TGYeBL4OXNtvnmuBr2TNQ8DMiDitgTVK0oA0MkRPBzbUjW+s2gY6jySN\nGm0jXcBgRMT11A73AfZExGpgDrB95KoaFm5Tcxhr2zTWtgcGt01nnMhMjQzRTcCCuvH5VdtA5yEz\nlwJL69siYllmdg5NqaOD29Qcxto2jbXtgeHdpkYezj8CLIqIMyOiHXgfcEe/ee4AfqW6Sn85sCsz\nNzewRkkakIb1RDOzOyJuBO4BWoFbMnNVRNxQTb8ZuAt4G7AW2Ad8qFH1SdJgNPScaGbeRS0o69tu\nrhtO4MODXPzS48/SdNym5jDWtmmsbQ8M4zZFLbckSYPhbZ+SVGBMhOjxbidtRhHxbESsiIjHImLZ\nSNczGBFxS0Rsi4iVdW0nR8Q/RcSPqp+zRrLGgTjK9nwiIjZV++mxiHjbSNY4UBGxICK+FxFPRsSq\niPitqr2Z99PRtmlY9lXTH85Xt5OuAa6m9uH8R4DrMvPJES2sUEQ8C3RmZtN+Xi8i3gDsoXYX2gVV\n258Cz2fmp6o3vFmZ+dGRrPNEHWV7PgHsycz/PZK1DVZ1R+BpmfloREwDlgPvAn6N5t1PR9um9zIM\n+2os9ERP5HZSjYDMfAB4vl/ztcCXq+EvU/vP3RSOsj1NLTM3Z+aj1fBu4Clqdwk283462jYNi7EQ\nomP1VtEE7o2I5dUdWmPF3LrP/m4B5o5kMUPkI9VTx25ppsPe/iLiFcBrgIcZI/up3zbBMOyrsRCi\nY9UVmXkxtSdbfbg6lBxTqo+0Nff5pNqTxs4CLgY2A58Z2XIGJyKmAt8EfjszX6if1qz76QjbNCz7\naiyE6AndKtpsMnNT9XMb8HfUTluMBVv7nsxV/dw2wvUUycytmdmTmb3A/6MJ91NETKAWNrdl5u1V\nc1PvpyNt03Dtq7EQoidyO2lTiYgp1QlxImIK8FZg5bF/q2ncAfxqNfyrwLdHsJZi/R7V+As02X6K\niAC+BDyVmX9WN6lp99PRtmm49lXTX50HqD6q8Fleup30j0a4pCIRcRa13ifU7ir7ajNuU0R8DbiK\n2hN0tgIfB74FfANYCKwD3puZTXGx5ijbcxW1w8MEngX+QzM97yEirgD+GVgB9FbN/53aOcRm3U9H\n26brGIZ9NSZCVJJGylg4nJekEWOISlIBQ1SSChiiklTAEJWkAoaoJBUwRDVkIqKnesTYyoi4MyJm\nDvD3PxERv1sN/8+IeMsQ1HRSRHy/etrXqBQRF0bErSNdhwbHENVQejEzL64eE/c8g/+qFzLzY5l5\n7xDU9OvA7ZnZcyIzV1+S2NC/i8xcAcyPiIWNXK+GhiGq4fIg1dO0ImJqRNwXEY9WD5o+/KjCiPi9\niFgTEf8CnFvXfmtEvLsafjYi5lTDnRFxfzV8Zd0Ddn/Yd6tsPx+gumXxaHVExCui9lDvr1C7FXBB\n1B70/WhEPB4R9x1rfRHxXyLikerpQH9Qtw2/UrU9HhF/VbW9p+qpPx4RD9TVeSe1W5bVbDLTl68h\neVF74C3Ubr/9G+CaarwNmF4Nz6H2ba4BXErt1rzJwPSq/Xer+W4F3l0NPwvMqYY7gfur4TuB11XD\nU4G2fvW0A1vqxo9Wxyuo3R54eTWtg9rjFc+sxk8+2vqoPddgabWcFuA7wBuAV1F7WPicfstYAZxe\nDc+sq+11wJ0jvQ99DfxlT1RD6aSIeIyXnj/5T1V7AH8cEU8A91Lroc4FXg/8XWbuy9qjygb64Jh/\nBf4sIn6TWiB195s+B9hZN360OgDWZeZD1fDlwAOZ+ROAfOme8SOt763V64fAo8B5wCLgTcDfZPXN\nBP2WcWtE/HtqbzZ9tgHzBrj9GgUMUQ2lF7P2DNQzqAVW3znRD1Dr3V1aTd8KTBrAcrt56f/q4d/L\nzE8B/w44CfjXiDivfz391nOsOvYer4ijrC+A/5W1c8EXZ+Y5mfmlYyzjBuB/UHt84/KImF23XS8e\nrwaNPoaohlxm7gN+E/idiGgDZgDbMvNQRLyRWsgCPAC8q7qCPg14x1EW+Sy1Q3+AX+prjIizM3NF\nZv4JtUcivixEM3MH0BoRfUF5tDr6ewh4Q0ScWa3n5GOs7x7g16sHABMRp0fEKcB3gff0hWS/ZTyc\nmR8DunjpWbivpMkeo6eatpEuQGNTZv6wOmy+DrgNuDMiVgDLgKereR6NiL8GHqd2OPvIURb3B8CX\nIuKTwP117b9dhWEvsAq4+wi/+4/AFdQO349YxxFq74raV7LcXl2p30btixB/Zn2ZeSAizgcerD3G\nkj3AL2fmqoj4I+D7EdFD7XD/14BPR8Qiaj3Y+6ptB3gj8PdH2X6NYj4KT2NaRFwC/OfM/OBI13I0\nETER+D61r4Tpf15Xo5yH8xrTsvatj98bzR+2p/bg45sM0OZkT1SSCtgTlaQChqgkFTBEJamAISpJ\nBQxRSSrw/wFJVSyanZnIPAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb53b48cbe0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (5,5)\n",
    "plot(rads, q_0_rad)\n",
    "xlabel(\"Radius (arcsecs)\")\n",
    "ylabel(\"$Q_0$\")\n",
    "ylim([0, 1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (py36)",
   "language": "python",
   "name": "py36"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}