Match_LOFAR_Q0.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ML match for LOFAR and the combined PanSTARRS WISE catalogue: Compute the $Q_0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Configuration\n",
    "\n",
    "### Load libraries and setup"
   ]
  },
  {
   "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\n",
    "from IPython.display import clear_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from mltier1 import Field, Q_0, parallel_process, describe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import clear_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Area limits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Busy week Hatfield 2017\n",
    "ra_down = 170.\n",
    "ra_up = 190.\n",
    "dec_down = 46.8\n",
    "dec_up = 55.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "field = Field(170.0, 190.0, 46.8, 55.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "combined = Table.read(\"pw.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#lofar_all = Table.read(\"data/LOFAR_HBA_T1_DR1_catalog_v0.9.srl.fits\")\n",
    "lofar_all = Table.read(\"data/LOFAR_HBA_T1_DR1_merge_ID_optical_v0.8.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['AllWISE', 'objID', 'ra', 'dec', 'raErr', 'decErr', 'W1mag',\n",
       "       'W1magErr', 'i', 'iErr'],\n",
       "      dtype='<U8')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(combined.colnames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Source_Name', 'RA', 'E_RA', 'DEC', 'E_DEC', 'Peak_flux',\n",
       "       'E_Peak_flux', 'Total_flux', 'E_Total_flux', 'Maj', 'E_Maj', 'Min',\n",
       "       'E_Min', 'DC_Maj', 'E_DC_Maj', 'DC_Min', 'E_DC_Min', 'PA', 'E_PA',\n",
       "       'DC_PA', 'E_DC_PA', 'Isl_rms', 'S_Code', 'Mosaic_ID',\n",
       "       'Masked_Fraction', 'ID_flag', 'ID_name', 'ID_ra', 'ID_dec', 'ML_LR',\n",
       "       'LGZ_Size', 'LGZ_Width', 'LGZ_PA', 'LGZ_Assoc', 'LGZ_Assoc_Qual',\n",
       "       'LGZ_ID_Qual', 'AllWISE', 'objID', 'gFApFlux', 'gFApFluxErr',\n",
       "       'gFApMag', 'gFApMagErr', 'rFApFlux', 'rFApFluxErr', 'rFApMag',\n",
       "       'rFApMagErr', 'iFApFlux', 'iFApFluxErr', 'iFApMag', 'iFApMagErr',\n",
       "       'zFApFlux', 'zFApFluxErr', 'zFApMag', 'zFApMagErr', 'yFApFlux',\n",
       "       'yFApFluxErr', 'yFApMag', 'yFApMagErr', 'w1Flux', 'w1FluxErr',\n",
       "       'w1Mag', 'w1MagErr', 'w2Flux', 'w2FluxErr', 'w2Mag', 'w2MagErr',\n",
       "       'w3Flux', 'w3FluxErr', 'w3Mag', 'w3MagErr', 'w4Flux', 'w4FluxErr',\n",
       "       'w4Mag', 'w4MagErr', 'XrayClass', '2RXS_ID', 'XMMSL2_ID', 'IRClass',\n",
       "       'z_spec', 'z_source', 'specAGN', 'mqcAGN', 'EBV', 'objName',\n",
       "       'z_best', 'z_best_source', 'z1_median', 'z1_min', 'z1_max',\n",
       "       'z1_area', 'z2_median', 'z2_min', 'z2_max', 'z2_area', 'chi_r_eazy',\n",
       "       'chi_r_atlas', 'chi_r_cosmos', 'chi_r_stellar', 'stellar_type',\n",
       "       'z_gpz', 'z_gpz_err'],\n",
       "      dtype='<U15')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(lofar_all.colnames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filter catalogues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lofar_aux = lofar_all[~np.isnan(lofar_all['Maj'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lofar = field.filter_catalogue(lofar_aux[(lofar_aux['Maj'] < 30.) &\n",
    "                                         (lofar_aux['ID_flag'] == 1)], \n",
    "                               colnames=(\"RA\", \"DEC\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sky coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "coords_combined = SkyCoord(combined['ra'], \n",
    "                           combined['dec'], \n",
    "                           unit=(u.deg, u.deg), \n",
    "                           frame='icrs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "coords_lofar = SkyCoord(lofar['RA'], \n",
    "                       lofar['DEC'], \n",
    "                       unit=(u.deg, u.deg), \n",
    "                       frame='icrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary of galaxy types in the combined catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8196213"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "combined_matched = (~np.isnan(combined[\"i\"]) & \n",
    "                    ~np.isnan(combined[\"W1mag\"]))\n",
    "np.sum(combined_matched) # Matches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13454849"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "combined_panstarrs = (~np.isnan(combined[\"i\"]) & \n",
    "                      np.isnan(combined[\"W1mag\"]))\n",
    "np.sum(combined_panstarrs) # Only PanSTARSS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5023475"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "combined_wise =(np.isnan(combined[\"i\"]) & \n",
    "                ~np.isnan(combined[\"W1mag\"]))\n",
    "np.sum(combined_wise) # Only WISE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "combined_i = combined_matched | combined_panstarrs\n",
    "combined_w1 = combined_matched | combined_wise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "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": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_iter = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rads = list(range(1,26))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### i-band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "q_0_comp_i = Q_0(coords_lofar, coords_combined[combined_i], field)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1 0.39416 +/- 0.00018 [0.39395 0.39453]\n",
      " 2 0.48341 +/- 0.00015 [0.48313 0.48368]\n",
      " 3 0.50341 +/- 0.00029 [0.50283 0.50389]\n",
      " 4 0.51244 +/- 0.00057 [0.51171 0.51382]\n",
      " 5 0.51767 +/- 0.00047 [0.51698 0.51835]\n",
      " 6 0.52033 +/- 0.00074 [0.51903 0.52178]\n",
      " 7 0.52212 +/- 0.00072 [0.52085 0.52325]\n",
      " 8 0.52557 +/- 0.00093 [0.52453 0.52728]\n",
      " 9 0.52914 +/- 0.00102 [0.52748 0.53145]\n",
      "10 0.53035 +/- 0.00071 [0.52862 0.53130]\n",
      "11 0.53297 +/- 0.00227 [0.52914 0.53675]\n",
      "12 0.53529 +/- 0.00171 [0.53281 0.53775]\n",
      "13 0.53556 +/- 0.00153 [0.53397 0.53931]\n",
      "14 0.53923 +/- 0.00239 [0.53543 0.54248]\n",
      "15 0.53991 +/- 0.00301 [0.53618 0.54549]\n",
      "16 0.54070 +/- 0.00308 [0.53607 0.54550]\n",
      "17 0.54471 +/- 0.00267 [0.53900 0.54840]\n",
      "18 0.54894 +/- 0.00429 [0.54378 0.55702]\n",
      "19 0.54513 +/- 0.00398 [0.54033 0.55521]\n",
      "20 0.54732 +/- 0.00452 [0.54058 0.55561]\n",
      "21 0.55536 +/- 0.00431 [0.54458 0.55987]\n",
      "22 0.55809 +/- 0.00564 [0.54508 0.56579]\n",
      "23 0.55459 +/- 0.00627 [0.54407 0.56626]\n",
      "24 0.56449 +/- 0.00719 [0.54647 0.57282]\n",
      "25 0.56499 +/- 0.00871 [0.54897 0.57933]\n"
     ]
    }
   ],
   "source": [
    "q_0_rad_i = []\n",
    "for radius in rads:\n",
    "    q_0_rad_aux = []\n",
    "    for i in range(n_iter):\n",
    "        out = q_0_comp_i(radius=radius)\n",
    "        q_0_rad_aux.append(out)\n",
    "    q_0_rad_i.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": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVEAAAFBCAYAAADOoJzZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHIZJREFUeJzt3XuUXnV97/H3Z+aZW5KZXMiFSxISJJabgDBcPKCCt4Ie\ni7dSwGoFjxx6xKrr6JLl6VJbq6c9ba3tEsvKqSnSU0WtYMGD4uUIVARNghASEIiQQAIkISFzYe7z\nfM8fe094MmSSeeY38zzzDJ/XWs/a1+z93ezkw77+tiICMzObmLpqF2BmVsscomZmCRyiZmYJHKJm\nZgkcomZmCRyiZmYJKhaiktZI2ilp4xjTJekfJG2WtEHSaZWqzcxsoip5JHo9cMFBpl8IrMp/VwL/\nWIGazMySVCxEI+IuYM9BZrkIuCEy9wLzJB1RmerMzCZmOl0TPQp4qmR4Wz7OzGzaKlS7gImQdCXZ\nKT+zZ88+/bjjjqtyRWY206xfv/65iFh0qPmmU4huB5aVDC/Nx71ERKwGVgO0t7fHunXrpr46M3tZ\nkbR1PPNNp9P5W4D353fpzwY6IuKZahdlZnYwFTsSlfRN4DxgoaRtwGeBBoCIuA64DXgrsBnoAS6v\nVG1mZhNVsRCNiEsPMT2AD1eoHDOzSTGdTufNzGqOQ9TMLIFD1MwsgUPUzCyBQ9TMLIFD1MwsgUPU\nzCyBQ9TMLIFD1MwsgUPUzCyBQ9TMLIFD1MwsgUPUzCyBQ9TMLIFD1MwsgUPUzCyBQ9TMLIFD1Mws\ngUPUzCyBQ9TMLIFD1MwsgUPUzCyBQ9TMLIFD1MwsgUPUzCyBQ9TMLIFD1MwsgUPUzCyBQ9TMLIFD\n1MwsgUPUzCyBQ9TMLIFD1MwsgUPUzCyBQ9TMLIFD1MwsgUPUzCyBQ9TMLIFD1MwsgUPUzCyBQ9TM\nLIFD1MwsgUPUzCyBQ9TMLIFD1MwsgUPUzCxBRUNU0gWSHpG0WdI1B5g+V9Ktkh6QtEnS5ZWsz8ys\nXBULUUn1wLXAhcAJwKWSThg124eBhyLiFOA84G8lNVaqRjOzclXySPRMYHNEPB4RA8CNwEWj5gmg\nVZKAOcAeYKiCNZqZlaWSIXoU8FTJ8LZ8XKmvAMcDTwMPAh+NiOLoBUm6UtI6Set27do1VfWamR3S\ndLux9LvA/cCRwKnAVyS1jZ4pIlZHRHtEtC9atKjSNZqZ7VPJEN0OLCsZXpqPK3U5cFNkNgNPAMdV\nqD4zs7JVMkTXAqskrcxvFl0C3DJqnieBNwJIWgL8DvB4BWs0MytLoVIrioghSVcDtwP1wJqI2CTp\nqnz6dcDngeslPQgI+FREPFepGs3MylWxEAWIiNuA20aNu66k/2ngLZWsycwsxXS7sWRmVlMcomZm\nCRyiZmYJHKJmZgkcomZmCRyiZmYJHKJmZgkcomZmCRyiZmYJHKJmZgkcomZmCRyiZmYJHKJmZgkc\nomZmCRyiZmYJHKJmZgkcomZmCRyiZmYJHKJmZgkcomZmCRyiZmYJHKJmZgkcomZmCRyiZmYJHKJm\nZgkcomZmCRyiZmYJHKJmZgkcomZmCRyiZmYJHKJmZgkcomZmCRyiZmYJHKJmZgkcomZmCRyiZmYJ\nHKJmZgkcomZmCRyiZmYJHKJmZgkcomZmCRyiZmYJHKJmZgkcomZmCRyiZmYJKhqiki6Q9IikzZKu\nGWOe8yTdL2mTpDsrWZ+ZWbkKlVqRpHrgWuDNwDZgraRbIuKhknnmAV8FLoiIJyUtrlR9ZmYTUckj\n0TOBzRHxeEQMADcCF42a5zLgpoh4EiAidlawPjOzslUyRI8CnioZ3paPK/VKYL6kOyStl/T+ilVn\nZjYBFTudH6cCcDrwRqAFuEfSvRHxaOlMkq4ErgRYvnx5xYs0MxtRySPR7cCykuGl+bhS24DbI+KF\niHgOuAs4ZfSCImJ1RLRHRPuiRYumrGAzs0M5ZIhKequkIydhXWuBVZJWSmoELgFuGTXPvwPnSipI\nmgWcBTw8Ces2M5sS4zmdfyfw55KWAL8BHgDuz7sPRcTweFYUEUOSrgZuB+qBNRGxSdJV+fTrIuJh\nST8ENgBF4J8iYmPZW2VmViGKiPHNKH2a7EbQJrLrlpcCeyJi6dSVd2jt7e2xbt26apZgZjOQpPUR\n0X6o+cq5sfQHEbHv+qSkrwKfnEhxZmYzRTk3ljolnT4yEBHryR5JMjN72SrnSPSDwE2S1gLrgVcB\ng1NSlZlZjRh3iEbEo5JOA95BFqAPA5+eqsLMzMayq6uftVv20Nk7yFAxKEYwXCz5RVAsBsNFGC4W\nGY7gI29YRXND/aTXMu4QlfQG4L3AXmAj2R307kmvyMxslOdfGODex3dzz+O7uee3u3lsZ3nRU18n\nPvTaY6obosAa4GNAA3Ay2RHpicCxk16VmdWM4WLQ2TvI3t5B9vYMsLd3kI6eQTp6B9nbM8je3gE6\ne4dobS6wqLWJhXMaWTinKfvlw02F/cOto3eQXz2xh3t+mwXnw890AjCrsZ72FQt412lLOfuYBSxp\na6ZQJ+rqRL2ybqFO1NeJOmnftKlUTohujYjv5f3fmYpizGz6Ghou8uiObjZs28sD2zrYsG0vT+3p\nobNv6KB/rrW5QGtTgc6+Ibr7DzxvW3MhD9QmegeG2fR0B8WApkId7Svm84m3vJLXvOIwTl46j4b6\n6dUMcjkhepekjwNfjvE+XGpmUyIi6OwdorNvkAWzG5ndNLnNYBSLwZbdL7BhWwcPbNvLhm0dbHq6\ng77BIpCF3slL59F+9HzmzWpk3qyG7NfSyNxZDcxraWDerEbamgsUSkKvb3CYXV39PNfdz3PdA1m3\nZHhXdz9zmgp85A2reM0rDuPVy+e95Ch1uinnv/wJZDeUPiVpPdlbS/dHhI9KzSZRsRjs6Rng2Y4+\nnuno49mO3rybDe/ozLq9gy++LDi7sZ7Fbc0sam1icWsTi1ubWdy2f3+dRHf/EF19g3T3DdHVP0RX\n3xDdfUN09w/SlY/b3d3Ppqc76cqPMJsb6jjpyLlcdubRnLJsLicvncfRC2ZN6DS5uaGeZQtmsWzB\nrEn771Vt5dydfzeApBZeDNSz8am92bh19w+xo7OPHR197Ojq49mOfnZ09rGzKwvJHZ397OzqY3B4\n/5O9Qp1Y0tbM4XObOf7INt5w3GIOn9tMW0sDu7sH2NnVx86ufnZ19rNxewc7u3bSMzCuN7IBaGmo\nZ05zgdbmAnNbGnj7KUdyytIsMFctnrPf0aTtr5y784cBFwN9ZK9+fisirp+iuswqYrgYPN8zwK6u\n/iyE8l9X3yDDxWBwOBguFhksBkPDRYaKwdBw5NOKWbcYDA4VGSoWGRwOhopFhoaz6SPzDwwX6ekf\n4oUDBFtrU4HFbU0cPreZs1YuYMncZg7PA/OIuVl34eymso/8uvuH2NmZheuOzr5sXc0F5jQ1MKcp\nC8xsuOCQTFDO6fzNwE+APwYeBV4j6bcRcfyUVGaWGy4G3X3Z9b+uvux0tKtkeGCoSBAUA4oRRGSn\nxEE2XIzsGuLgcLC7u59d3S+G5e4XBhguvvQSf6FOFOpFoa4u72b99XWioV55Nxsu1NfRkA83N2Td\nQj48soyGetHSWJ8dTbZlp9eHtzWzpK150q9njpjTVGDOojkcs2jOlCzfMuXsvdaI+HNJ74qI10t6\nNwdo69OsVETQP1TMf8N09Q3R0Zs9/jLyGMzIozD7xuePxIyE5Fh3dMtRp+xZwcNmN7GotYklbc2c\ndORcFrU27fdbnN8hnqpgs5mnnL8pfXm3X1JLRHxX0ieBz0xBXTZNjNwF3tXdx66u7O7pc139+3X3\n9gzuC8n+wRcDs3+oyMBQcVzrmd1Yz9yWBtpaGpjb0sDRh82iraUhP+VsoK25QFtzyXBL1m1tLtBU\nqKNO2XOBEvt16wTS1D4naC9v5YTo30haAHwLWCPpF8C8qSnLplIxvw74XPfAvtPbkcdNdr/k0ZMB\nBoZfGoQN9dr3wPS8WQ0saqinqVBHU6Gepoa6F/sLdflwPY2FuiwM86As/U23Z//Mxqucu/PfzXu/\nJOl9ZHfn3zUlVdmY+gaH9536jrwd0tk7yAv9Q/QMDtPTP8wLA0P09A/nw0PZ8MAwPQPDdPQOsucg\n1wEPK3mbZNXi1n2nuQvnNGb9c7LhuS0NPsIzY4IfqouIf5nsQl7uOnoHeWpPD1t39/Dknuy3q6tv\nv+uFe3sHD3l6XCeY3VRgVmM9sxsLzGqqZ1ZDgQWzG1k2P7sbm4VkIwtbm/JrhFlwOhjNylfOI06v\nA/4COIzs0yB/FxFrp6qwmaizb5BN2zvZuvsFntzTw9Y9PTyVB+benv1bFTxsdiOL25qZ19LAsYvn\nZKe9s7JT33ktjVl31ounwyPB2VSocxCaVVC5DZD8MdmbSqcDfy/pyxHx7SmprMZFBNv39rJ+6/Os\n3bKHdVue55EdXYy8MFuoE0vnt7BswSz+88lHsHzBLJYvmJ11D5vFHN8dNqsJ5fxLfS4ifpz3/1DS\nz4F7AYco2bOMDz/TuS801299nmc6sgcaZjfWc9rR87ngpMN59fL5HLNwNkfMbfYDzmYzwCFDVNIN\nwH3AzyV9BvhiRAwB/bz42NPL1oPbOvjKzx7j7s279z3PeHhbM+0r5nPGigWcfvR8jju81YFpNkON\n50j0a2QP1S8AzgeukLQZWAncNIW1TWu/ebaTv/vxo9y+aQfzZjVw0alHcubKLDSPmtfi65JmLxOH\nDNGIuBO4c2RYUgE4DjiVl+EbS7/d1c2Xf/IY39/wNHMaC3z8Ta/kinNX0NrcUO3SzKwKyr57kZ/K\nb8x//2fSK5qmntzdw9//9DFu/vU2mhvq+W/nvYIPvfYY5s1qrHZpZlZFvgV8CE/v7eUrP9vMt9c+\nRX2duOKclVx13itYOKep2qWZ2TRQdohKentE3DoVxUwnnX2DfOlHj/KNXz5JEFx21nI+fP6xLGlr\nrnZpZjaNTORI9AvAjA/RL3z/Yf7tvm2857SlfOSNx7J0/sxpidvMJs9EQnTG33be3d3Pzfdv5w/O\nWMYX3/mqapdjZtPYRB5enPEfqfvmr55kYKjI5f9pRbVLMbNpzk+AjzIwVOSGe7by2lULWbWktdrl\nmNk05xAd5Qcbn2FnVz9XnLuy2qWYWQ2YSIjumPQqpomIYM3Pn+CYhbN5/apF1S7HzGpA2SEaEW+e\nikKmg/ue3MsD2zq4/JwVE/qmtpm9/Ph0vsSau5+gtbnAu05bWu1SzKxGOERzT+/t5Ycbn+XSM5f7\nS49mNm6HDFFJn8u750iasberb7hnKxHB+19zdLVLMbMaMp5Drtvz7keBkyQ1AA8BG4ANEfGdqSqu\nUnoGhvjmr57kd0883G8mmVlZxtMU3j1592IASU3AiWRf+zwLqPkQvfnX2+noHeTyc/xYk5mVZyJN\n4fWTtXR/3+SXU3kRwT/fvYWTjmrjjBXzq12OmdWYibTidBnwe8AQ2TXVWyPim5NdWKX8x2PPsXln\nN1+6+BS3Rm9mZZvIbejXR8QlIwOSrgVqNkTX3P0EC+c08baTj6h2KWZWgyYSok2S3gY8BSwFWia3\npMr57a5u7nhkFx9/0ytpKtRXuxwzq0HjClFJRwJvBGYB3wIWkt1Y2gZcPWXVTbHr795CY30dl521\nvNqlmFmNGs8nk98CfB24g+wzyVeRhenlEfGLKa1uCnX0DPJv67fxe6ceyaJWf+rDzCZmPEeifwG8\nNiI2j4yQ9BpgtaQPAj0RsXGqCpwq31r3JL2Dw1x+zopql2JmNWw8r302lgYo7Ht29N1kX/sc900l\nSRdIekTSZknXHGS+MyQNSXrPeJddjqHhIl//xVbOWrmAE4+cOxWrMLOXifGEaJ+kl7QLFxGPAsNk\n10oPSVI9cC1wIXACcKmkE8aY76+AH41nuRPx44d2sH1vrx+uN7Nk4wnRvwa+l99c2kfSQqA/InaO\nc11nApsj4vGIGABuBC46wHwfAb4LjHe5Zfvnu7ewdH4Lbz5hyVStwsxeJsbz2ud381c975G0HngA\naAQuJrteOl5HkT0WNWIb2Wuj+0g6CngncD5wRhnLHreN2zv41ZY9/OnbjqfebYaaWaJxNYUXEd8A\njge+D8wFBoHLIuLrk1zPl4FPRUTxYDNJulLSOknrdu3aVdYK1tz9BLMb67n4jGUpdZqZAWU8bB8R\nPcCahHVtB0qTa2k+rlQ7cGP++uVC4K2ShiLie6NqWQ2sBmhvbx/310d3dvVx6wNPc9mZy2lrbpjA\nJpiZ7a+SrQ+vBVZJWkkWnpcAl5XOEBH77vRIuh74/ugATdHa1MBn334i5x67cLIWaWYvcxUL0YgY\nknQ1Wfuk9cCaiNgk6ap8+nVTXUNLYz1/eLYbXTazyVPR72BExG3AbaPGHTA8I+IDlajJzCyFv7Fk\nZpbAIWpmlsAhamaWwCFqZpbAIWpmlsAhamaWwCFqZpbAIWpmlsAhamaWwCFqZpbAIWpmlsAhamaW\nwCFqZpbAIWpmlsAhamaWwCFqZpbAIWpmlsAhamaWwCFqZpbAIWpmlsAhamaWwCFqZpbAIWpmlsAh\namaWwCFqZpbAIWpmlsAhamaWwCFqZpbAIWpmlsAhamaWwCFqZpbAIWpmlsAhamaWwCFqZpbAIWpm\nlsAhamaWwCFqZpbAIWpmlsAhamaWwCFqZpbAIWpmlsAhamaWwCFqZpbAIWpmlsAhamaWwCFqZpag\noiEq6QJJj0jaLOmaA0x/r6QNkh6U9AtJp1SyPjOzclUsRCXVA9cCFwInAJdKOmHUbE8Ar4+IVwGf\nB1ZXqj4zs4mo5JHomcDmiHg8IgaAG4GLSmeIiF9ExPP54L3A0grWZ2ZWtkqG6FHAUyXD2/JxY/kg\n8IMDTZB0paR1ktbt2rVrEks0MyvPtLyxJOl8shD91IGmR8TqiGiPiPZFixZVtjgzsxKFCq5rO7Cs\nZHhpPm4/kk4G/gm4MCJ2V6g2M7MJqeSR6FpglaSVkhqBS4BbSmeQtBy4CXhfRDxawdrMzCakYkei\nETEk6WrgdqAeWBMRmyRdlU+/DvgMcBjwVUkAQxHRXqkazczKpYiodg1J2tvbY926ddUuw8xmGEnr\nx3MQNy1vLJmZ1QqHqJlZAoeomVkCh6iZWQKHqJlZAoeomVkCh6iZWQKHqJlZAoeomVkCh6iZWQKH\nqJlZAoeomVkCh6iZWQKHqJlZAoeomVkCh6iZWQKHqJlZAoeomVkCh6iZWQKHqJlZAoeomVkCh6iZ\nWQKHqJlZAoeomVkCh6iZWQKHqJlZAoeomVkCh6iZWQKHqJlZAoeomVkCh6iZWQKHqJlZAoeomVkC\nh6iZWQKHqJlZAoeomVkCh6iZWQKHqJlZAoeomVkCh6iZWQKHqJlZAoeomVkCh6iZWQKHqJlZAoeo\nmVmCioaopAskPSJps6RrDjBdkv4hn75B0mmVrM/MrFwVC1FJ9cC1wIXACcClkk4YNduFwKr8dyXw\nj5Wqz8xsIip5JHomsDkiHo+IAeBG4KJR81wE3BCZe4F5ko6oYI1mZmWpZIgeBTxVMrwtH1fuPGZm\n00ah2gVMhKQryU73AbolPQIsBJ6rXlVTwttUG2baNs207YGJbdPR45mpkiG6HVhWMrw0H1fuPETE\namB16ThJ6yKifXJKnR68TbVhpm3TTNsemNptquTp/FpglaSVkhqBS4BbRs1zC/D+/C792UBHRDxT\nwRrNzMpSsSPRiBiSdDVwO1APrImITZKuyqdfB9wGvBXYDPQAl1eqPjOziajoNdGIuI0sKEvHXVfS\nH8CHJ7j41YeepeZ4m2rDTNummbY9MIXbpCy3zMxsIvzap5lZghkRood6nbQWSdoi6UFJ90taV+16\nJkLSGkk7JW0sGbdA0o8lPZZ351ezxnKMsT2fk7Q930/3S3prNWssl6Rlkn4m6SFJmyR9NB9fy/tp\nrG2akn1V86fz+eukjwJvJns4fy1waUQ8VNXCEknaArRHRM0+ryfpdUA32VtoJ+Xj/hewJyL+Mv8f\n3vyI+FQ16xyvMbbnc0B3RPxNNWubqPyNwCMi4j5JrcB64B3AB6jd/TTWNl3MFOyrmXAkOp7XSa0K\nIuIuYM+o0RcBX8/7v072l7smjLE9NS0inomI+/L+LuBhsrcEa3k/jbVNU2ImhOhMfVU0gJ9IWp+/\noTVTLCl59vdZYEk1i5kkH8lbHVtTS6e9o0laAbwa+CUzZD+N2iaYgn01E0J0pjo3Ik4la9nqw/mp\n5IySP9JW29eTspbGjgFOBZ4B/ra65UyMpDnAd4GPRURn6bRa3U8H2KYp2VczIUTH9aporYmI7Xl3\nJ3Az2WWLmWDHSMtceXdnletJEhE7ImI4IorA/6YG95OkBrKw+deIuCkfXdP76UDbNFX7aiaE6Hhe\nJ60pkmbnF8SRNBt4C7Dx4H+qZtwC/FHe/0fAv1exlmSjmmp8JzW2nyQJ+BrwcER8qWRSze6nsbZp\nqvZVzd+dB8gfVfgyL75O+oUql5RE0jFkR5+QvVX2jVrcJknfBM4ja0FnB/BZ4HvAt4HlwFbg4oio\niZs1Y2zPeWSnhwFsAf5rLbX3IOlc4D+AB4FiPvrTZNcQa3U/jbVNlzIF+2pGhKiZWbXMhNN5M7Oq\ncYiamSVwiJqZJXCImpklcIiamSVwiJqZJXCI2qSRNJw3MbZR0q2S5pX55z8n6RN5/59LetMk1NQi\n6c68ta9pSdKrJF1f7TpsYhyiNpl6I+LUvJm4PUz8Uy9ExGci4ieTUNMVwE0RMTyemfOPJFb030VE\nPAgslbS8kuu1yeEQtalyD3lrWpLmSPqppPvyhqb3NVUo6X9IelTSz4HfKRl/vaT35P1bJC3M+9sl\n3ZH3v76kgd1fj7wqO8p7yV9ZHKsOSSuUNep9A9mrgMuUNfR9n6QHJP30YOuT9ElJa/PWgf6sZBve\nn497QNK/5ON+Pz9Sf0DSXSV13kr2yrLVmojwz79J+ZE1eAvZ67ffAS7IhwtAW96/kOxrrgJOJ3s1\nbxbQlo//RD7f9cB78v4twMK8vx24I++/FTgn758DFEbV0wg8WzI8Vh0ryF4PPDuftoisecWV+fCC\nsdZH1q7B6nw5dcD3gdcBJ5I1Fr5w1DIeBI7K++eV1HYOcGu196F/5f98JGqTqUXS/bzY/uSP8/EC\nvihpA/ATsiPUJcBrgZsjoieypsrKbTjmbuBLkv6ELJCGRk1fCOwtGR6rDoCtEXFv3n82cFdEPAEQ\nL74zfqD1vSX//Rq4DzgOWAW8AfhO5F8mGLWM6yV9iOx/NiN2AkeWuf02DThEbTL1RtYG6tFkgTVy\nTfS9ZEd3p+fTdwDNZSx3iBf/ru77cxHxl8B/AVqAuyUdN7qeUes5WB0vHKqIMdYn4H9Gdi341Ig4\nNiK+dpBlXAX8KVnzjeslHVayXb2HqsGmH4eoTbqI6AH+BPjvkgrAXGBnRAxKOp8sZAHuAt6R30Fv\nBd4+xiK3kJ36A7x7ZKSkV0TEgxHxV2RNIu4XohHxPFAvaSQox6pjtHuB10lama9nwUHWdztwRd4A\nMJKOkrQY+H/A74+E5Khl/DIiPgPs4sW2cF9JjTWjZ5lCtQuwmSkifp2fNl8K/Ctwq6QHgXXAb/J5\n7pP0LeABstPZtWMs7s+Ar0n6PHBHyfiP5WFYBDYBPzjAn/0RcC7Z6fsB6zhA7buUfZLlpvxO/U6y\nDyG+ZH0R0S/peOCerBlLuoE/jIhNkr4A3ClpmOx0/wPAX0taRXYE+9N82wHOB/7vGNtv05ibwrMZ\nTdJpwMcj4n3VrmUskpqAO8k+CTP6uq5Ncz6dtxktsq8+/mw6P2xP1vDxNQ7Q2uQjUTOzBD4SNTNL\n4BA1M0vgEDUzS+AQNTNL4BA1M0vw/wEGtsYWKQ+zfgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f32894d45f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (5,5)\n",
    "plot(rads, q_0_rad_i)\n",
    "xlabel(\"Radius (arcsecs)\")\n",
    "ylabel(\"$Q_0 i-band$\")\n",
    "ylim([0, 1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### W1-band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "q_0_comp_w1 = Q_0(coords_lofar, coords_combined[combined_w1], field)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1 0.48528 +/- 0.00013 [0.48502 0.48546]\n",
      " 2 0.63373 +/- 0.00010 [0.63358 0.63391]\n",
      " 3 0.67774 +/- 0.00018 [0.67742 0.67798]\n",
      " 4 0.69997 +/- 0.00025 [0.69936 0.70034]\n",
      " 5 0.71345 +/- 0.00042 [0.71300 0.71422]\n",
      " 6 0.72368 +/- 0.00034 [0.72299 0.72415]\n",
      " 7 0.73504 +/- 0.00033 [0.73461 0.73573]\n",
      " 8 0.74619 +/- 0.00047 [0.74556 0.74696]\n",
      " 9 0.75823 +/- 0.00030 [0.75797 0.75892]\n",
      "10 0.76719 +/- 0.00046 [0.76662 0.76809]\n",
      "11 0.77414 +/- 0.00034 [0.77366 0.77459]\n",
      "12 0.77706 +/- 0.00069 [0.77573 0.77812]\n",
      "13 0.78045 +/- 0.00086 [0.77838 0.78151]\n",
      "14 0.78377 +/- 0.00064 [0.78290 0.78479]\n",
      "15 0.78532 +/- 0.00083 [0.78416 0.78703]\n",
      "16 0.78580 +/- 0.00095 [0.78438 0.78740]\n",
      "17 0.78524 +/- 0.00114 [0.78276 0.78676]\n",
      "18 0.78352 +/- 0.00110 [0.78125 0.78569]\n",
      "19 0.78376 +/- 0.00123 [0.78201 0.78612]\n",
      "20 0.78243 +/- 0.00187 [0.77898 0.78465]\n",
      "21 0.77945 +/- 0.00138 [0.77685 0.78243]\n",
      "22 0.77808 +/- 0.00188 [0.77507 0.78093]\n",
      "23 0.77228 +/- 0.00195 [0.76921 0.77522]\n",
      "24 0.77022 +/- 0.00298 [0.76421 0.77438]\n",
      "25 0.76521 +/- 0.00238 [0.76185 0.76875]\n"
     ]
    }
   ],
   "source": [
    "q_0_rad_w1 = []\n",
    "for radius in rads:\n",
    "    q_0_rad_aux = []\n",
    "    for i in range(n_iter):\n",
    "        out = q_0_comp_w1(radius=radius)\n",
    "        q_0_rad_aux.append(out)\n",
    "    q_0_rad_w1.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": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVEAAAFBCAYAAADOoJzZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH2ZJREFUeJzt3XuUXWWd5vHvU/dUJSGpVBFyhQRiAogEKG6CAtog0NPS\ntkID3lC7kR61255pl66eWXa3Tl9tHaenRSatCPZFtFtUcHAQlYsG0CTcAwQrISEJJJWqSqAuqftv\n/ti7krJMJefUrjqnTuX5rHXW2fvdu875bU94fPft3YoIzMxsfMqKXYCZWSlziJqZZeAQNTPLwCFq\nZpaBQ9TMLAOHqJlZBgULUUm3SmqR9MwYyyXpHyQ1S3pK0pmFqs3MbLwK2RO9Dbj8MMuvAFakrxuB\nLxegJjOzTAoWohHxENB+mFWuAr4eiUeBOZIWFKY6M7PxmUrHRBcB20fM70jbzMymrIpiFzAekm4k\n2eWnrq7urFWrVhW5IjObbjZs2NAaEY1HWm8qhehOYMmI+cVp26+JiDXAGoCmpqZYv3795FdnZkcV\nSdtyWW8q7c7fBbwvPUt/HvBqRLxS7KLMzA6nYD1RSd8ALgYaJO0A/gyoBIiIW4B7gCuBZqAb+ECh\najMzG6+ChWhEXHeE5QF8pEDlmJlNiKm0O29mVnIcomZmGThEzcwycIiamWXgEDUzy8AhamaWgUPU\nzCwDh6iZWQYOUTOzDByiZmYZOETNzDJwiJqZZeAQNTPLwCFqZpaBQ9TMLAOHqJlZBg5RM7MMHKJm\nZhk4RM3MMnCImpll4BA1M8vAIWpmloFD1MwsA4eomVkGDlEzswwcomZmGThEzcwycIiamWXgEDUz\ny8AhamaWgUPUzCwDh6iZWQYOUTOzDByiZmYZOETNzDJwiJqZZeAQNTPLwCFqZpaBQ9TMLAOHqJlZ\nBg5RM7MMHKJmZhk4RM3MMnCImpll4BA1M8vAIWpmlkFBQ1TS5ZI2SWqW9KlDLD9G0t2SnpS0UdIH\nClmfmVm+ChaiksqBLwFXAKcA10k6ZdRqHwGejYjTgYuBz0uqKlSNZmb5KmRP9BygOSK2REQfcAdw\n1ah1ApglScBMoB0YKGCNZmZ5KWSILgK2j5jfkbaN9I/AycDLwNPAH0XE0OgPknSjpPWS1u/Zs2ey\n6jUzO6KpdmLpbcATwEJgNfCPkmaPXiki1kREU0Q0NTY2FrpGM7MDChmiO4ElI+YXp20jfQC4MxLN\nwIvAqgLVZ2aWt4oCftc6YIWkZSTheS1w/ah1XgLeCvxU0nxgJbClgDWaATA0FPQMDLK/b5D9/YP0\n9A+yv2+I/f3J/P6+tK1/kIHBIQaHgsFI/m5gKBiKSNrS6YGhIAIqykT5iFfFr02XUVEmqivLqK2q\noLaqnBlV5dSNmq6pLCM5dWDFVrAQjYgBSR8F7gXKgVsjYqOkm9LltwCfBW6T9DQg4JMR0VqoGm36\nGhwK9nX30dbVR2tnL22dfbR19tLWlbS1Dbelyzt6Ju58ZpmgvEwIMZiGa1YSzKgsp6aynIggSAI8\ngAgYijj4DkQE5WWiprKcmopyairLqB5+Tz+nuqKMmspyZlSWMbe2ijm1VdTXVabvVcytrTzQXl7m\nAB9WyJ4oEXEPcM+otltGTL8MXFbImqw0RQSdvQPs6+7/lRDccyAMew+EZWtnH+1dvRwqu8oEc2ur\nmDezinl11Zy6cDYNM6uZXVPBjOHeX2U5Nen7jMpyZlSVMaOyghlVSQhVlJUlvUmJ8vLkvayMZL5M\nv9ZjjDRIB0b0UgcH40DA9g8O0TswRHfvIN19A3T3Dx6Y3t8/SFfvIPv7BujuG6RnYBAhygSSkKBM\nQkBZWfI+3D44FPSkveqe/iF6+gfpHUjeX93fT2+6rLtvkH37++kb+LVzukAS4LNrKqmvq2JWTQU1\nFeVUVyYBnIT08HQZMyrLqa4sp7aqnEVzZrC8sY4l9bVUV5RP/D+KIiloiJodSVfvAJt2d7BlTxft\nXb3s7e5nX3cfe7v6ae/uS6bTtv7BQ/foZlZXpKFYxZL6Ws5YOoeGmdXU11XRMLOaeTPT97ri9Kok\nUVEupnKORAT7+wdp7+pjX3c/7V197O3uY29XH+3p//7tXX109g7Q0z9IR88Aezp6Dwb0wMGwHq1M\nsHhuLcsa6ljWUMfyxjqWN8xkWWMdC2bXUFZivVyHqBXF4FCwra2L53d1JK9XXuP5XR281N79K+tV\nlos5tcmu5JzaKpY3zGRuuos53NaQ9iKHe5MzqqZwOpUISekx2QoWzx3/50RE0qvuG2R7ezcvtnax\nZU8nW1q7eLG1i3Vb2+nuGzywfnVFGYvmzGBmTQUzq0e8aiqoS6dn1VRQV5W0HTe7hiX1tcytrSza\nMWKHqE26gcEhnt/VwYZte9n48qts2tXBpt0dB3opZYJlDXWctvgYrmlazKrjZnPSsTOZN7OKmdUV\nPoFSwiQd2M2vr6vi9CVzfmV5RNDS0cuWPV1sae3kxT1dvPJaD129A3T2DPBSVzedvQN09Q7Q0TPA\nwBjHk+uqyllSX8viubUsnjuDJfW1LBl+r69lZvXkRZ1D1CZcR08/j720jw1b21m/bS9PbN93oLcx\nr66KkxfM5t3nHs+q42Zx8oIkMGsq3Xs8Gkli/uwa5s+u4fwT5x123eFebVfvAJ1pqL7yag/b27vZ\nvreb7e372bG3m0c2t9I1oncLMLe2kh/+8UU0zqqe8G1wiFomEcGOvfvZsG0v67e1s37rXjbt7iAi\n6WGevGA2V5+1mLNOqKfp+LksnDOj2CVbiRrZq503MwnD1y865tfWiwj2dvcfCNeX2rvZsXc/9XWT\nMwyHQ9TyNjgUPP7SXn747G7ue3Y3L7Z2AckJnTOWzuHy1x9H0/H1rF46Z1J3o8wORRL1dVWHPHww\nGfwv3HLS0z/IT3/Zyn3P7uLHz7XQ1tVHZbk4b/k83n/+8Zy9rJ5Vx8329YN21HGI2pjau/r48XNJ\nb/Onv2xlf/8gs6oruHjVsVx2ynwuWtnI7JrKYpdpVlQOUTsgInh+VwcPbNrD/c+3sH5bO0MBC46p\n4eqmxVx6ynzOXTaPqoqpNm6NWfE4RI9yHT39rG1u5YFNe3hg0x52vdYDJCeEPnrJSVx6ynG8ftFs\nX2ZkNgaH6FEmInhhdycPbGrh/k0trN+6l4GhYFZ1BReuaOCSlcdy0cpG5s+uKXapZiXBIXoU6O4b\nYG1zG/dvauGB51t4+dWkt7nquFn83puWc/HKRs46fi6V5d5NN8uXQ3Sa2tbWxf3Pt/CTTXt4dEsb\nfQND1FWVc8FJDXzsrSu4eGUjC47xNZtmWTlEp4m+gSHWb23nJ8+38JNNLWzZk1y7ubyhjveedzxv\nWXUsTSfMnVaj55hNBQ7REhURbG3rZm1zK2ubW/npL1vp7B2gqryMc5fX897zjueSlcdyQkNdsUs1\nm9YcoiWkpaOHh5vbWNvcysOb29i5bz+QXIL0W6cv4JKVx3LBSQ3U+S4hs4Lxf21T2Gs9/fx8S3sa\nmq28sLsTgGNmVHL+8nncdPGJXHDiPJY11PkSJLMicYhOMW2dvdzzzC7ufvJlNmzby+BQUFNZxtkn\n1PM7Zy7mghMbOGWhb680myocolNAR08/927czV1Pvsza5lYGh4KTjp3JTRct58KTGjnz+Dk+IWQ2\nRTlEi6Snf5CfPN/CXU+8zE82tdA3MMTiuTO48c3LefvpC1l13CzvopuVAIdoAQ0NBQ/+cg93PfEy\nP9y4i66+QRpmVnP9OUt5++qFnLFkjoPTrMQ4RAvk51va+Mz3n2Xjy68xu6aC//SGhbx99ULOWz7P\nxzfNSphDdJJtb+/mr3/wHPc8vYuFx9TwP3/3dK48bYGPcZpNEw7RSdLZO8DN9zfzlZ+9SLnEf7n0\ndfz+m5b7SZRm04xDdIINDQX/sWEHf3fvJlo7e/mdMxbxictX+j51s2nKITqBRh73PHPpHL7y/iZW\nF+AZL2ZWPA7RCTD6uOf/unY1bz99oc+0mx0FHKIZbdrVwdW3PEz/YPi4p9lRyCGawa5Xe7jha7+g\nprKc73/sjSydV1vsksyswByi4/RaTz83fO0XdPQM8M0Pn+cANTtK+XkQ49A3MMQf/MsGmls6+fJ7\nzuTUhccUuyQzKxL3RPMUEXzy20+xtrmNz199Om9a0VjsksysiNwTzdPn7t3Edx7fySfetpJ3nrW4\n2OWYWZE5RPPwz49u4+YHNnP9uUv5zxefWOxyzGwKcIjm6L5nd/Nn33uGt646ls+8/VRfA2pmgEM0\nJ4+/tJePfeMxTlt0DP/7+jOo8PPZzSzlNDiCF1u7+NDt6zl2Vg1fveFsaqt8Ls7MDnKIHkZrZy83\nfO0XANz+wXNomFld5IrMbKpxiI6hu2+AD922jt2v9fDV9zexzM9vN7ND8L7pGNY8tIWndr7K/3nP\nWZyxdG6xyzGzKco90TE8+MIeVi+Zw2WnHlfsUsxsCnOIHkJHTz9P7XiVC05sKHYpZjbFOUQP4edb\n2hkcCt540rxil2JmU5xD9BDWbm6luqKMM30s1MyO4IghKulKSQsLUcxU8XBzG2efUE9NpQdXNrPD\ny6Un+g7gLknbJd0n6e8lvUfSaZLyShlJl0vaJKlZ0qfGWOdiSU9I2ijpwXw+fyLs6ehl0+4O78qb\nWU6OeIlTRPw+gKQ/BRYBW4BLgDVAO5DTUEZp4H4JuBTYAayTdFdEPDtinTnAzcDlEfGSpGPz25zs\nHt7cCsCFJ/mkkpkdWT7Xif5uRJw+PCPpZuATefz9OUBzRGxJ//4O4Crg2RHrXA/cGREvAURESx6f\nPyHWNrcyu6bCAy2bWU7yObH0mqSzhmciYgPwujz+fhGwfcT8jrRtpNcBcyU9IGmDpPfl8fmZRQRr\nm9s4/8R5lJd5lCYzO7J8eqIfAu6UtA7YAJwG9E9CPWcBbwVmAI9IejQiXhi5kqQbgRsBli5dOmFf\n/lJ7Nzv37efDFy2fsM80s+kt555oGmRnAj8A5gPPAVfm8V07gSUj5henbSPtAO6NiK6IaAUeAk4f\ntQ4RsSYimiKiqbFx4h7Psba5DYA3+iJ7M8tRziEq6S3Al4FzSU4u/RTozOO71gErJC2TVAVcC9w1\nap3vARdKqpBUm37Xc3l8RyZrN7cyf3Y1JzZ6sBEzy00+u/O3Ah8HKoE3AL8NnAqclMsfR8SApI8C\n9wLlwK0RsVHSTenyWyLiOUn/D3gKGAK+EhHP5FHjuA0NBY9sbuPi1zV61Hozy1k+IbotIr6bTv/7\neL4sIu4B7hnVdsuo+c8BnxvP52fx/K4O2rv6eKMvbTKzPORzdv4hSX+sadpNG74+9AJfZG9mecin\nJ3oKyRn5T0raADwBPBER4+qVTjVrm1tZ3lDHgmNmFLsUMyshOYdoRLwTQNIMDgbqeYxz134q6R8c\n4hcvtvOOM0dftmpmdng5h6ikecA1QA+wEfhmRNw2SXUV1JPb99HVN+jxQ80sb/kcE/0O0Aj8FcmJ\nn1clFezyo8m0trkNCc4/0cdDzSw/+YTorIj4DLA7Ii4CrmMa7MpDcn3oqQtnM6e2qtilmFmJySdE\ne9L3XkkzIuLbwGWTUFNBdfcN8PhLe70rb2bjks/Z+b+XVA98C7hV0sPAnMkpq3DWbd1L/2D4+lAz\nG5d87p3/dkS0R8TnSS6YXwL8zqRVViAPN7dSWS7OPsGPAjGz/OVz7/ybJT0kaSNwBfAfIwdULlVr\nN7dyxtK51Fbl0yk3M0vkc0z0VuCzwMXA7cAXJV0zGUUVyr7uPja+/JqPh5rZuOXT/WqNiPvS6Xsl\nrQUeJTlGWpIe2dxGhG/1NLPxy+Vpn1+X9HHgZ5I+LWk4eHs5eMa+JK3d3EpdVTmnLyn582NmViS5\n7M5/lWRYunqSZyI1S/oR8Dxw/yTWNunWNrdxzrJ6KsvzOaphZnZQLk/7fBA48OjitCe6CljNIUad\nLxUv79vPi61dvPvciXu8iJkdffI+JR0RA8Az6etfJryiAlnbnAx950eBmFkWR+1+7MOb26ivq2LV\ncbOKXYqZlbCjMkSTRyO3cv6J8yjzo5HNLIO8Q1TSb01GIYW0eU8nLR29XOhbPc0so/H0RP9ywqso\nsOFHI/siezPLajwhWvL7v2ubW1k8dwZL59UWuxQzK3HjCdGY8CoKaHAoeHRLm3uhZjYhjroTS8/s\nfJXXegZ4o2/1NLMJcNSF6NrNvj7UzCbOeEJ094RXUUAPN7excv4sGmdVF7sUM5sG8g7RiLh0Mgop\nhJ7+QdZtbfeuvJlNmKNqd76zd4DffMMCfuPk+cUuxcymiaNqOPeGmdV84ZrVxS7DzKaRTD1RSR+Y\nqELMzEpR1t35v5iQKszMStQRd+clPTXWIsAHF83sqJbLMdH5wNuAvaPaBTw84RWZmZWQXEL0+8DM\niHhi9AJJD0x4RWZmJSSXEP3DiOg61IKIuH6C6zEzKym5nFjaJOlmSW+Y9GrMzEpMLiG6Engc+Iqk\ntZLeJ8n3TJqZkUOIRkRXRPxTRJwDfAQ4H3hO0uclrZj0Cs3MprBcLnE6CZgNzErffwZsBn4P+DhQ\nPpkFmplNZbmcWHoB2Al8h+Qyp06gA/hs+m5mdtTKJUTPBD4MXALcAXw9Ikp6ODwzs4mSyzHRJyLi\nD4DzgBbgu5K+Jektk16dmdkUl88oTkPA94D7gbcAN0uKiDh5UiozMysBuZxY2kvycLou4LX01QE8\nk06bmR21cumJ1kdEST/h08xssuRyTNQBamY2hrzHE5V0vaQ7JP2rpH+TdF0ef3u5pE2SmiV96jDr\nnS1pQNK78q3PzKyQxvN4kIsi4trhGUlfAr5xpD+SVA58CbgU2AGsk3RXRDx7iPX+FvjhOGozMyuo\n8YRotaTfBLYDi4EZOf7dOUBzRGwBkHQHcBXw7Kj1PgZ8Gzh7HLWZmRVUTiEqaSHwVqAW+CbQAJxG\n0qP8aI7ftYgkeIftAM4d9T2LgHeQXNjvEDWzKS+XS5wuA24HHgB6gZtIwvQDETHRI9t/EfhkRAxJ\nOlxNNwI3AixdunSCSzAzy10uPdH/AbwpIpqHGySdD6yR9CGgOyKeyeFzdgJLRswvTttGagLuSAO0\nAbhS0kBEfHfkShGxBlgD0NTU5KsHzKxocgnRqpEBChARj0h6J8mjQ3pJdu2PZB2wQtIykvC8FviV\nkfEjYtnwtKTbgO+PDlAzs6kkl0uceiQ1jm6MiBeAQZJjpUcUEQMkx0/vBZ4DvhURGyXdJOmmPGo2\nM5sycumJfo5k0JGrI+Ll4UZJDUBvRLTk+mURcQ9wz6i2W8ZY94ZcP9fMrFiOGKIR8e30cSCPSNoA\nPAlUAdeQHC81Mztq5XTHUkT8G3AyyTHQY4B+4PqIuH0SazMzm/Jyvtg+IrqBWyexFjOzkpP3vfNm\nZnaQQ9TMLAOHqJlZBg5RM7MMHKJmZhk4RM3MMnCImpll4BA1M8vAIWpmloFD1MwsA4eomVkGDlEz\nswwcomZmGThEzcwycIiamWXgEDUzy8AhamaWgUPUzCwDh6iZWQYOUTOzDByiZmYZOETNzDJwiJqZ\nZeAQNTPLwCFqZpaBQ9TMLAOHqJlZBg5RM7MMHKJmZhk4RM3MMnCImpll4BA1M8vAIWpmloFD1Mws\nA4eomVkGDlEzswwcomZmGThEzcwycIiamWXgEDUzy8AhamaWgUPUzCwDh6iZWQYFDVFJl0vaJKlZ\n0qcOsfzdkp6S9LSkhyWdXsj6zMzyVbAQlVQOfAm4AjgFuE7SKaNWexG4KCJOAz4LrClUfWZm41HI\nnug5QHNEbImIPuAO4KqRK0TEwxGxN519FFhcwPrMzPJWyBBdBGwfMb8jbRvLh4AfHGqBpBslrZe0\nfs+ePRNYoplZfqbkiSVJl5CE6CcPtTwi1kREU0Q0NTY2FrY4M7MRKgr4XTuBJSPmF6dtv0LSG4Cv\nAFdERFuBajMzG5dC9kTXASskLZNUBVwL3DVyBUlLgTuB90bECwWszcxsXArWE42IAUkfBe4FyoFb\nI2KjpJvS5bcAnwbmATdLAhiIiKZC1Whmli9FRLFryKSpqSnWr19f7DLMbJqRtCGXTtyUPLFkZlYq\nHKJmZhk4RM3MMnCImpll4BA1M8vAIWpmloFD1MwsA4eomVkGDlEzswwcomZmGThEzcwycIiamWXg\nEDUzy8AhamaWgUPUzCwDh6iZWQYOUTOzDByiZmYZOETNzDJwiJqZZeAQNTPLwCFqZpaBQ9TMLAOH\nqJlZBg5RM7MMHKJmZhk4RM3MMnCImpll4BA1M8vAIWpmloFD1MwsA4eomVkGDlEzswwcomZmGThE\nzcwycIiamWXgEDUzy8AhamaWgUPUzCwDh6iZWQYOUTOzDByiZmYZOETNzDJwiJqZZeAQNTPLoKAh\nKulySZskNUv61CGWS9I/pMufknRmIeszM8tXwUJUUjnwJeAK4BTgOkmnjFrtCmBF+roR+HKh6jMz\nG49C9kTPAZojYktE9AF3AFeNWucq4OuReBSYI2lBAWs0M8tLIUN0EbB9xPyOtC3fdczMpoyKYhcw\nHpJuJNndB+iUtAloAFqLV9Wk8DaVhum2TdNte2B823R8LisVMkR3AktGzC9O2/Jdh4hYA6wZ2SZp\nfUQ0TUypU4O3qTRMt22abtsDk7tNhdydXweskLRMUhVwLXDXqHXuAt6XnqU/D3g1Il4pYI1mZnkp\nWE80IgYkfRS4FygHbo2IjZJuSpffAtwDXAk0A93ABwpVn5nZeBT0mGhE3EMSlCPbbhkxHcBHxvnx\na468SsnxNpWG6bZN0217YBK3SUlumZnZePi2TzOzDKZFiB7pdtJSJGmrpKclPSFpfbHrGQ9Jt0pq\nkfTMiLZ6SfdJ+mX6PreYNeZjjO35c0k709/pCUlXFrPGfElaIul+Sc9K2ijpj9L2Uv6dxtqmSfmt\nSn53Pr2d9AXgUpKL89cB10XEs0UtLCNJW4GmiCjZ6/UkvRnoJLkL7fVp298B7RHxN+n/4c2NiE8W\ns85cjbE9fw50RsTfF7O28UrvCFwQEY9JmgVsAH4buIHS/Z3G2qZrmITfajr0RHO5ndSKICIeAtpH\nNV8F3J5O307yj7skjLE9JS0iXomIx9LpDuA5krsES/l3GmubJsV0CNHpeqtoAD+StCG9Q2u6mD/i\n2t9dwPxiFjNBPpaOOnZrKe32jibpBOAM4OdMk99p1DbBJPxW0yFEp6sLI2I1ychWH0l3JaeV9JK2\n0j6elIw0thxYDbwCfL645YyPpJnAt4GPR8RrI5eV6u90iG2alN9qOoRoTreKlpqI2Jm+twDfITls\nMR3sHh6ZK31vKXI9mUTE7ogYjIgh4J8owd9JUiVJ2PxrRNyZNpf073SobZqs32o6hGgut5OWFEl1\n6QFxJNUBlwHPHP6vSsZdwPvT6fcD3ytiLZmNGqrxHZTY7yRJwFeB5yLiCyMWlezvNNY2TdZvVfJn\n5wHSSxW+yMHbSf+yyCVlImk5Se8TkrvK/q0Ut0nSN4CLSUbQ2Q38GfBd4FvAUmAbcE1ElMTJmjG2\n52KS3cMAtgIfLqXxHiRdCPwUeBoYSpv/lOQYYqn+TmNt03VMwm81LULUzKxYpsPuvJlZ0ThEzcwy\ncIiamWXgEDUzy8AhamaWgUPUzCwDh6hNGEmD6RBjz0i6W9KcPP/+zyX9STr9GUm/MQE1zZD0YDra\n15Qk6TRJtxW7Dhsfh6hNpP0RsTodJq6d8T/qhYj4dET8aAJq+iBwZ0QM5rJy+pDEgv53ERFPA4sl\nLS3k99rEcIjaZHmEdDQtSTMl/VjSY+lA0weGKpT03yS9IOlnwMoR7bdJelc6vVVSQzrdJOmBdPqi\nEQPsPj58q+wo7ya9ZXGsOiSdoGRQ76+T3Aq4RMlA349JelLSjw/3fZI+IWldOjrQX4zYhvelbU9K\n+ue07eq0p/6kpIdG1Hk3yS3LVmoiwi+/JuRFMuAtJLff/jtweTpfAcxOpxtInuYq4CySW/Nqgdlp\n+5+k690GvCud3go0pNNNwAPp9N3ABen0TKBiVD1VwK4R82PVcQLJ7YHnpcsaSYZXXJbO14/1fSTj\nGqxJP6cM+D7wZuBUksHCG0Z9xtPAonR6zojaLgDuLvZv6Ff+L/dEbSLNkPQEB8efvC9tF/BXkp4C\nfkTSQ50PvAn4TkR0RzJUWb4Dx6wFviDpD0kCaWDU8gZg34j5seoA2BYRj6bT5wEPRcSLAHHwnvFD\nfd9l6etx4DFgFbACeAvw75E+mWDUZ9wm6fdJ/s9mWAuwMM/ttynAIWoTaX8kY6AeTxJYw8dE303S\nuzsrXb4bqMnjcwc4+G/1wN9FxN8AvwfMANZKWjW6nlHfc7g6uo5UxBjfJ+CvIzkWvDoiToqIrx7m\nM24C/jvJ8I0bJM0bsV37j1SDTT0OUZtwEdEN/CHwXyVVAMcALRHRL+kSkpAFeAj47fQM+izgt8b4\nyK0ku/4A7xxulHRiRDwdEX9LMiTir4RoROwFyiUNB+VYdYz2KPBmScvS76k/zPfdC3wwHQAYSYsk\nHQv8BLh6OCRHfcbPI+LTwB4OjoX7OkpsGD1LVBS7AJueIuLxdLf5OuBfgbslPQ2sB55P13lM0jeB\nJ0l2Z9eN8XF/AXxV0meBB0a0fzwNwyFgI/CDQ/ztD4ELSXbfD1nHIWrfo+SRLHemZ+pbSB6E+Gvf\nFxG9kk4GHkmGsaQTeE9EbJT0l8CDkgZJdvdvAD4naQVJD/bH6bYDXAL83zG236YwD4Vn05qkM4E/\njoj3FruWsUiqBh4keSTM6OO6NsV5d96mtUie+nj/VL7YnmTg4085QEuTe6JmZhm4J2pmloFD1Mws\nA4eomVkGDlEzswwcomZmGfx/fASb+KV6+S4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3319f91898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (5,5)\n",
    "plot(rads, q_0_rad_w1)\n",
    "xlabel(\"Radius (arcsecs)\")\n",
    "ylabel(\"$Q_0 W1-band$\")\n",
    "ylim([0, 1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}