Match_LOFAR_Q0_gaus.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": false
   },
   "outputs": [],
   "source": [
    "from mltier1 import Field, Q_0, parallel_process"
   ]
  },
  {
   "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": false
   },
   "outputs": [],
   "source": [
    "combined = Table.read(\"pw.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lofar_all = Table.read(\"data/LOFAR_HBA_T1_DR1_catalog_v0.1.gaus.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', 'E_RA_tot', 'DEC', 'E_DEC',\n",
       "       'E_DEC_tot', 'Peak_flux', 'E_Peak_flux', 'E_Peak_flux_tot',\n",
       "       'Total_flux', 'E_Total_flux', 'E_Total_flux_tot', 'Maj', 'E_Maj',\n",
       "       'Min', 'E_Min', 'PA', 'E_PA', 'Isl_rms', 'S_Code', 'Mosaic_ID',\n",
       "       'Gaus_id', 'Isl_id', 'Source_id'], \n",
       "      dtype='<U16')"
      ]
     },
     "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": false
   },
   "outputs": [],
   "source": [
    "lofar = field.filter_catalogue(lofar_all, colnames=(\"RA\", \"DEC\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sky coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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": 15,
   "metadata": {
    "collapsed": false
   },
   "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": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1339600"
      ]
     },
     "execution_count": 16,
     "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": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3071319"
      ]
     },
     "execution_count": 17,
     "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": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "821044"
      ]
     },
     "execution_count": 18,
     "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": 19,
   "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": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_iter = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rads = list(range(1,26))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### i-band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "q_0_comp_i = Q_0(coords_lofar, coords_combined[combined_i], field)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1 0.32293 +/- 0.00011 [0.32274 0.32311]\n",
      " 2 0.40706 +/- 0.00016 [0.40681 0.40736]\n",
      " 3 0.43463 +/- 0.00039 [0.43407 0.43531]\n",
      " 4 0.45163 +/- 0.00071 [0.45053 0.45277]\n",
      " 5 0.46427 +/- 0.00061 [0.46291 0.46513]\n",
      " 6 0.47457 +/- 0.00067 [0.47386 0.47586]\n",
      " 7 0.48394 +/- 0.00060 [0.48312 0.48503]\n",
      " 8 0.49283 +/- 0.00082 [0.49112 0.49394]\n",
      " 9 0.50143 +/- 0.00109 [0.50009 0.50333]\n",
      "10 0.50901 +/- 0.00136 [0.50735 0.51253]\n",
      "11 0.51483 +/- 0.00094 [0.51348 0.51624]\n",
      "12 0.52188 +/- 0.00130 [0.51980 0.52365]\n",
      "13 0.52403 +/- 0.00208 [0.51971 0.52710]\n",
      "14 0.53285 +/- 0.00126 [0.53122 0.53519]\n",
      "15 0.53359 +/- 0.00265 [0.52952 0.53965]\n",
      "16 0.53520 +/- 0.00201 [0.53171 0.53942]\n",
      "17 0.53920 +/- 0.00303 [0.53580 0.54704]\n",
      "18 0.54199 +/- 0.00291 [0.53741 0.54699]\n",
      "19 0.54320 +/- 0.00321 [0.53766 0.54784]\n",
      "20 0.54872 +/- 0.00190 [0.54557 0.55240]\n",
      "21 0.55626 +/- 0.00458 [0.54802 0.56240]\n",
      "22 0.55035 +/- 0.00712 [0.53887 0.56393]\n",
      "23 0.56126 +/- 0.00697 [0.54816 0.57107]\n",
      "24 0.56854 +/- 0.00476 [0.56130 0.57421]\n",
      "25 0.56970 +/- 0.00821 [0.55572 0.57858]\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": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVEAAAFBCAYAAADOoJzZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHqBJREFUeJzt3Xl0XOWZ5/HvI8myZFmLbcm2LO/GYBsHjK2AEwhLEggw\nkyYrAdJJIEwMmZCEnEkOdGaaztLp6Z7uzjB9mjRxJx5IJx0gB0IgDSHLJCZsieUFvACO8G7LtrxJ\ntnapnvnjXtuFsGyVXlWVSvp9zqlTd3PVcynz83vvfe97zd0REZGByct2ASIiuUwhKiISQCEqIhJA\nISoiEkAhKiISQCEqIhIgYyFqZivMbL+ZbehjvZnZP5lZvZm9YmaLM1WbiMhAZbIl+gBw9WnWXwPM\njV/LgH/JQE0iIkEyFqLu/ixw6DSbXAf8wCMvARVmVp2Z6kREBmYonROtAXYmze+Kl4mIDFkF2S5g\nIMxsGdEhPyUlJUvmzZuX5YpEZLhZvXr1AXevOtN2QylEdwPTkuanxsvewt2XA8sBamtrva6uLv3V\niciIYmbb+7PdUDqcfwL4ZHyVfinQ5O4N2S5KROR0MtYSNbMfA5cDlWa2C/grYBSAu98PPAVcC9QD\nrcAtmapNRGSgMhai7n7jGdY78LkMlSMiMiiG0uG8iEjOUYiKiARQiIqIBFCIiogEUIiKiARQiIqI\nBFCIiogEUIiKiARQiIqIBFCIiogEUIiKiARQiIqIBFCIiogEUIiKiARQiIqIBFCIiogEUIiKiARQ\niIqIBFCIiogEUIiKiARQiIqIBFCIiogEUIiKiARQiIqIBFCIiogEUIiKiARQiIqIBFCIiogEUIiK\niARQiIqIBFCIiogEUIiKiARQiIqIBFCIiogEUIiKiARQiIqIBFCIiogEUIiKiARQiIqIBFCIiogE\nUIiKiARQiIqIBFCIiogEUIiKiARQiIqIBFCIiogEyGiImtnVZva6mdWb2d2nWF9uZk+a2ctmttHM\nbslkfSIiqcpYiJpZPnAfcA2wALjRzBb02uxzwCZ3Px+4HPhHMyvMVI0iIqnKZEv0QqDe3be4eyfw\nEHBdr20cKDUzA8YCh4DuDNYoIpKSTIZoDbAzaX5XvCzZPwPzgT3AeuCL7p7o/UFmtszM6sysrrGx\nMV31ioic0VC7sPQ+YB0wBVgE/LOZlfXeyN2Xu3utu9dWVVVlukYRkRMyGaK7gWlJ81PjZcluAR7z\nSD2wFZiXofpERFKWyRBdBcw1s1nxxaIbgCd6bbMDeA+AmU0CzgG2ZLBGEZGUFGTqi9y928zuAJ4B\n8oEV7r7RzG6P198PfBN4wMzWAwbc5e4HMlWjiEiqMhaiAO7+FPBUr2X3J03vAa7KZE0iIiGG2oUl\nEZGcohAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBUR\nCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmg\nEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAV\nEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCZDREDWz\nq83sdTOrN7O7+9jmcjNbZ2YbzWxlJusTEUlVQaa+yMzygfuAK4FdwCoze8LdNyVtUwF8B7ja3XeY\n2cRM1SciMhCZbIleCNS7+xZ37wQeAq7rtc1NwGPuvgPA3fdnsD4RkZRlMkRrgJ1J87viZcnOBsaZ\n2e/MbLWZfTJj1YmIDEDGDuf7qQBYArwHKAZeNLOX3H1z8kZmtgxYBjB9+vSMFykiclwmW6K7gWlJ\n81PjZcl2Ac+4e4u7HwCeBc7v/UHuvtzda929tqqqKm0Fi4icyRlD1MyuNbMpg/Bdq4C5ZjbLzAqB\nG4Anem3zM+ASMyswszHARcCrg/DdIiJp0Z/D+Q8C3zCzScBrwMvAuvh9k7v39OeL3L3bzO4AngHy\ngRXuvtHMbo/X3+/ur5rZL4BXgATwPXffkPJeiYhkiLl7/zY0+yrRhaCNROctbwQOufvU9JV3ZrW1\ntV5XV5fNEkRkGDKz1e5ee6btUrmw9DF3P3F+0sy+A3xlIMWJiAwXqVxYajazJcdn3H01UZckEZER\nK5WW6K3AY2a2ClgNvA3oSktVIiJ9ONzSyasNzby69yhNrZ109CTo6ErQGb93dPfQ2Z2gozsRv/fQ\n0Z3gkdveQcnowe/V2e9PdPfNZrYY+ABRgL4KfHXQKxIRARIJZ+fhVjbtaWZTQ/OJ94am9jdtV1iQ\nx+gTr3xGF+QlLctnTGEB48bkkejn9Z9U9TtEzezdwMeBI8AGoivox9JSlYiMOIdbOnmu/gB12w6x\nqaGZVxuOcqyjG4D8PGNOVQkXzRrPgillLKguZ151KRNKCjGzrNadStt2BXAnMAo4j6hFei5wVhrq\nEpEhrrm9i67uBOMHGGQ9CeflXUdY+XojKzc38squIyQcSgrzmV9dxocW17CguowFU8o4e1IpRaPy\n07AX4VIJ0e3u/ng8/ZN0FCMiQ9OBYx1s3NPMht1NbNrTzIY9TWw/2ApAWVEBs6vGMruqhDlVY5ld\nWcLsqrHMmDDmLcG3r7mdlZsbeXZzI7//0wGa2rowg0XTKvj8u+dy2TlVnD+1gvy87LYuU5FKiD5r\nZl8C7vX+di4VkZySSDgNze1s3N3Exj3NbNzTxIbdzextPnkecvr4MSysKeP62mkUj8pny4FjbGls\n4YX6gzy25uSd3GYwdVwxsyvHUl1exLqdR3ht71EAqkpHc+WCSVx2dhWXnFXJuJLCjO/rYEklRBcQ\nXVC6y8xWE921tM7d1SoVySGJhLO3uZ1tB1vYfrCVbQdaTk4fbKG9KwFAnsGcqrEsnT2ehTXlnDul\nnAVTyigvHtXnZ7d0dLP1QAtvNEbBuuVAC1saj7Fu5xEWVJdx19XzuOzsKuZXl2b9XOZgSeXq/IcB\nzKyYk4G6FB3ai2RMIuEc7ejmSGsnR1q7aG7virv1nOzK09EVvx9fFq/f29zO9jgsO7oTJz6zMD+P\n6RPGMHPCGC4+q5KZlSUsqC5jfnUpYwpT6xJUMrqAhTXlLKwpH+xdH7JSuTo/AbgeaCe69fNhd38g\nTXWJjAjuTlNbF/uaO9jb3M6+5nb2N7dzqKUrCsq2Lg63dtLU2sWRtmhZIoWTaXkGRaPyKSzIo3Ls\naGZOKOGys6uYMaGEmRNKmFk5hury4pw6BznUpPLPzE+BXwOfBTYD7zCzN9x9floqE8kRXT0JWjt7\naO/qobWzh7bOHtq6umnrTNDa2U1bV7Ssuf1kWO5vbmdfcwf7mtvf1Co8buzoAirGjKJizCjGjSmk\npqL4xHR5cfReMWYUZcWjKCrIZ/SoN/eTHD0qj8L8PAry9SzKdEslREvd/Rtm9iF3v8zMPswpxvoU\nyUUd3T1s3NPMmu2HWb39MLuPtNHV43T3JOhOON2JBN09Hi07MR2t60mhaTimMJ/JZUVMLBvNBdMr\n4ukiJpWNZnJZEZPKiqgqHT1ku/PIW6USoscvz3WYWbG7P2pmXwHuSUNdIml14FjHicBcvf0wr+xu\nojNuEU4bH11RHpWfR0GeUZBvSdN5jMo3CvLyKMg3CvKM4lH5FBdGrzGF+fF8QdJ09F5aVMDY0QXD\n5oKKRFIJ0X8ws/HAw8AKM3sBqEhPWSIDk0g4LZ3dHOvo5lh7N0fj92Md3Rw81sHanUdYs/0w2+I+\njoX5eSysKeNT75jBkhnjWDx9HBPLirK8F5JLUrk6/2g8+W0z+wTR1fkPpaUqkT50dPfwp33HeG3v\nUV5taOa1vc3sb+44EZrHOrs5XS/myrGFLJ4+jhsvnM6SGeNYWFOuQ2cJMqAhTdz93wa7EJFk7s7+\nox1samjmtYaTgflGY8uJc5CjC/I4Z3Ipc6rGUlpUQGnRKMYWFVA6uoCx8aHz8fnSolGUFRcwuaxI\nh9MyqFLp4nQp8NfABKJHg/xvd1+VrsJkZOjo7mHHwda4U3YLWw8cY+uBFur3H+Nw68mRFmsqipk3\nuZSrFkxmXnUp8yaXMauyRF1zJOtSHYDks0R3Ki0B/o+Z3evuj6SlMhlWmlq72Linifr4TpatB6LX\nrsOtb+r3WFU6mlmVJbzv3MnMm1zK/Ooy5k0uo3xM33fJiGRTKiF6wN1/FU//wsyeA14CFKLyJo1H\nO9iwJx6oYncTG/Y0sfNQ24n1JYX5zKoq4fxpFXzwghpmV5Uwq7KEmZUllBUpLCW3nDFEzewHwBrg\nOTO7B/gbd+8GOjjZ7UlGqP3N7azbeYQNe5rZGAfmvuaOE+tnVZZw3tQKbrpwBudOKeOcyaVMLB2t\n85IybPSnJfp9ok7144ErgE+bWT0wC3gsjbXJENPZnWDjnibW7DjC2h2HWbvjCLuPRC3MPIOzJo7l\n4jmVnFtTzsIpZcyfUqaWpQx7ZwxRd18JrDw+b2YFwDxgEbpjaVhraGpj7Y6oX+WaHYfZsKf5RIf0\nKeVFXDBjHLdcPJMLplewoLqc4kJ1FZKRJ+UuTvGh/Ib49cNBr0iyoqsnwWsNR6nbfujEXTzHn2VT\nWJDHeTXlfOodM1g8fRwXTB/H5HJ1SBeBAfYTldzX1NrFmp2HWb0tCsx1O4/Q1tUDRK3M43fvLJ4x\njgXVZRQWaCALkVNJOUTN7P3u/mQ6ipH0cHd2HW7jj1sPnWhpbt4XPWMwP89YUF3Gx94+jSUzxrFk\nxjimVBRnuWKR3DGQlui3AIXoEObuvNF4jD9sPcQf49fxQ/PSogKWzBjH+8+bwpKZ4zh/akVansUt\nMlIM5P8e9U0ZYnoSzqsNzXFoHmTVtsMcaukEos7rF84az0WzxnPhrPGcPbGUPN3lIzJoBhKiekhd\nlrk79fuP8Vz9AZ6vP8AfthziaPx87unjx3DFORNPhOaMCWPUJ1MkjXQclyMamtp4vv4gz8fBuf9o\n1KF95oQxvH/RlBOhWV2u85kimaQQHaKa27t48Y0oNJ+rP8CWxhYAJpQU8s6zKrnkrAm8c04l08aP\nyXKlIiPbQEJ036BXIQAcbunkV5v28dSGBp6vP0BXj1M8Kp+LZo/npgunc/FZlZwzSec0RYaSgXS2\nvzIdhYxU+4+288uN+3h6QwMvbTlET8KpqSjm5nfO5D3zJ7F4+jj10RQZwnQ4nwUNTW38YsNenl6/\nl1XbD+EOsytLuO3S2VyzsJqFNWW6GCSSIxSiGdLe1cPP1u3moVU7WbvjCADnTCrlC++ey7Vvq+bs\nSWMVnCI5qD9D4X3N3b9mZhcDr7j70QzUNWzsa27nhy9t50d/2MGhlk7OmVTKV953DtcsnMzsqrHZ\nLk9EAvWnJfpM/P5FYKGZjQI2Aa8QhepP0lVcLlu/q4kVz2/l56/soTvhvHf+JD598SyWzh6vFqfI\nMNKfofBejN+vBzCz0cC5RE/7vAhQiMa6exL8ctM+/u/zW1m17TAlhfn8+dIZ3PzOmcyYUJLt8kQk\nDQZydb6DaKT7NYNfTm5qauvi4VU7ePCF7ew+0sa08cX85X9ewEdrp2pQYpFhbiCjON0E/BnQDeQB\nT7r7jwe7sFzx9PoG/uKn6znS2sXS2eO55/0LeO/8SXoKpcgIMZCr85e5+w3HZ8zsPmDEhWhLRzdf\nf3Ijj9Tt4ryp5fzw1rexsKY822WJSIYNJERHm9l/AnYCU4ERd7P2up1HuPOhtWw/1MrnrpjDne89\nm1H56hAvMhL1K0TNbArwHmAM8DBQSXRhaRdwR9qqG2J6Es53flvPvb/5E5PLinjoM0u5aPaEbJcl\nIlnUn36iVwEPAr8jekzy7URheou7v5DW6oaQnYda+dLD66jbfpjrFk3hG9ctpLxYF41ERrr+tET/\nGniXu9cfX2Bm7wCWm9mtQKu7b0hXgUPB42t385ePR7t478cW8YELarJckYgMFf05kVeYHKBwou/o\nh4me9tnvi0pmdrWZvW5m9WZ292m2e7uZdZvZR/r72enQ1NbFF368ljsfXse86lKe+uK7FKAi8ib9\naYm2m1mVuzcmL3T3zWbWQ3Su9IzMLB+4D7iS6FzqKjN7wt03nWK7vwN+2Z/PTZfX9jZz6wN17G1u\n58tXnc1nLz9L3ZZE5C360xL9e+Dx+OLSCWZWCXS4+/5+fteFQL27b3H3TuAh4LpTbPd54FGgv5+b\nFt/8+Sbau3p49LPv5I53z1WAisgp9ee2z0fjWz1fNLPVwMtAIXA90fnS/qoh6hZ13C6i20ZPMLMa\n4IPAFcDbU/jsQbV+VxPP1x/kL66Zx6JpFdkqQ0RyQL86N7r7vwPzgZ8D5UAXcJO7PzjI9dwL3OXu\nidNtZGbLzKzOzOoaGxtPt+mAfPfZNygdXcBNF00f9M8WkeGl353t3b0VWBHwXbuBaUnzU+NlyWqB\nh+JRjiqBa82s290f71XLcmA5QG1t7aA+fXT7wRaeWt/AskvnUKr73kXkDDI5KPMqYK6ZzSIKzxuA\nm5I3cPdZx6fN7AHg570DNN2+9/utFOTlccvFMzP5tSKSozIWou7ebWZ3EI1Pmg+scPeNZnZ7vP7+\nTNXSl4PHOnikbicfvKCGSWVF2S5HRHJARh8P4u5PAU/1WnbK8HT3mzNRU7IHX9hGZ0+CZZfNzvRX\ni0iO0qgZsZaObh58cTtXzp/EHD22Q0T6SSEae6RuJ01tXdx22ZxslyIiOUQhCnT1JPje77dy4czx\nLJkxLtvliEgOUYgC//FKA7uPtHGbzoWKSIpGfIi6O/evfIO5E8dyxTkTs12OiOSYER+iKzc38tre\no9x22RzydH+8iKRoxIfod1duYXJZEX92/pQzbywi0suIDtGXdx7hxS0HufWSWRQWjOj/FCIyQCM6\nOb777BuUFhVwowYaEZEBGrEhuu1AC09v2Msnls5g7OiM3rglIsPIiA3Rf/39Fkbl53GzBhoRkQAj\nMkQbj3bwk9W7+PDiqUws1UAjIjJwIzJEH3xhG109CT7zrlln3lhE5DRGXIi2dHTzgxe38b4Fk5mt\ngUZEJNCIC9Ef/3EHze3dusVTRAbFiArRrp4E339uKxfNGs8F0zXQiIiEG1Eh2ni0g+ryIm7XcHci\nMkhGVAfJKRXFPPZfL8Z9UJ9tJyIj2IhqiR4XP01URCTYiAxREZHBohAVEQmgEBURCaAQFREJoBAV\nEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJ\noBAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCaAQFREJoBAVEQmgEBURCZDR\nEDWzq83sdTOrN7O7T7H+42b2ipmtN7MXzOz8TNYnIpKqjIWomeUD9wHXAAuAG81sQa/NtgKXufvb\ngG8CyzNVn4jIQGSyJXohUO/uW9y9E3gIuC55A3d/wd0Px7MvAVMzWJ+ISMoyGaI1wM6k+V3xsr7c\nCjx9qhVmtszM6sysrrGxcRBLFBFJzZC8sGRmVxCF6F2nWu/uy9291t1rq6qqMluciEiSggx+125g\nWtL81HjZm5jZecD3gGvc/WCGahMRGZBMtkRXAXPNbJaZFQI3AE8kb2Bm04HHgE+4++YM1iYiMiAZ\na4m6e7eZ3QE8A+QDK9x9o5ndHq+/H7gHmAB8x8wAut29NlM1ioikytw92zUEqa2t9bq6umyXISLD\njJmt7k8jbkheWBIRyRUKURGRAApREZEAClERkQAKURGRAApREZEAClERkQAKURGRAApREZEAClER\nkQAKURGRAApREZEAClERkQAKURGRAApREZEAClERkQAKURGRAApREZEAClERkQAKURGRAApREZEA\nClERkQAKURGRAApREZEAClERkQAKURGRAApREZEAClERkQAKURGRAApREZEAClERkQAKURGRAApR\nEZEAClERkQAKURGRAApREZEAClERkQAKURGRAApREZEAClERkQAKURGRAApREZEAClERkQAKURGR\nAApREZEAGQ1RM7vazF43s3ozu/sU683M/ile/4qZLc5kfSIiqcpYiJpZPnAfcA2wALjRzBb02uwa\nYG78Wgb8S6bqExEZiEy2RC8E6t19i7t3Ag8B1/Xa5jrgBx55Cagws+oM1igikpJMhmgNsDNpfle8\nLNVtRESGjIJsFzAQZraM6HAf4JiZvQ5UAgeyV1VaaJ9yw3Dbp+G2PzCwfZrRn40yGaK7gWlJ81Pj\nZalug7svB5YnLzOzOnevHZxShwbtU24Ybvs03PYH0rtPmTycXwXMNbNZZlYI3AA80WubJ4BPxlfp\nlwJN7t6QwRpFRFKSsZaou3eb2R3AM0A+sMLdN5rZ7fH6+4GngGuBeqAVuCVT9YmIDERGz4m6+1NE\nQZm87P6kaQc+N8CPX37mTXKO9ik3DLd9Gm77A2ncJ4tyS0REBkK3fYqIBBgWIXqm20lzkZltM7P1\nZrbOzOqyXc9AmNkKM9tvZhuSlo03s1+Z2Z/i93HZrDEVfezP18xsd/w7rTOza7NZY6rMbJqZ/dbM\nNpnZRjP7Yrw8l3+nvvYpLb9Vzh/Ox7eTbgauJOqcvwq40d03ZbWwQGa2Dah195ztr2dmlwLHiO5C\nWxgv+1/AIXf/2/gfvHHuflc26+yvPvbna8Axd/+HbNY2UPEdgdXuvsbMSoHVwAeAm8nd36mvfbqe\nNPxWw6El2p/bSSUL3P1Z4FCvxdcBD8bTDxL95c4JfexPTnP3BndfE08fBV4lukswl3+nvvYpLYZD\niA7XW0Ud+LWZrY7v0BouJiX1/d0LTMpmMYPk8/GoYyty6bC3NzObCVwA/IFh8jv12idIw281HEJ0\nuLrE3RcRjWz1ufhQcliJu7Tl9vmkaKSx2cAioAH4x+yWMzBmNhZ4FLjT3ZuT1+Xq73SKfUrLbzUc\nQrRft4rmGnffHb/vB35KdNpiONh3fGSu+H1/lusJ4u773L3H3RPAv5KDv5OZjSIKmx+5+2Px4pz+\nnU61T+n6rYZDiPbndtKcYmYl8QlxzKwEuArYcPo/lTOeAD4VT38K+FkWawnWa6jGD5Jjv5OZGfB9\n4FV3/3bSqpz9nfrap3T9Vjl/dR4g7qpwLydvJ/1WlksKYmaziVqfEN1V9u+5uE9m9mPgcqIRdPYB\nfwU8DjwCTAe2A9e7e05crOljfy4nOjx0YBtwWy6N92BmlwC/B9YDiXjxV4nOIebq79TXPt1IGn6r\nYRGiIiLZMhwO50VEskYhKiISQCEqIhJAISoiEkAhKiISQCEqIhJAISqDxsx64iHGNpjZk2ZWkeKf\n/5qZfTme/oaZvXcQaio2s5XxaF9Dkpm9zcweyHYdMjAKURlMbe6+KB4m7hADf9QL7n6Pu/96EGr6\nNPCYu/f0Z+P4IYkZ/f/C3dcDU81seia/VwaHQlTS5UXi0bTMbKyZ/cbM1sQDTZ8YqtDM/ruZbTaz\n54BzkpY/YGYfiae3mVllPF1rZr+Lpy9LGmB37fFbZXv5OPEti33VYWYzLRrU+wdEtwJOs2ig7zVm\n9rKZ/eZ032dmXzGzVfHoQF9P2odPxsteNrN/i5d9NG6pv2xmzybV+STRLcuSa9xdL70G5UU04C1E\nt9/+BLg6ni8AyuLpSqKnuRqwhOjWvDFAWbz8y/F2DwAfiae3AZXxdC3wu3j6SeDieHosUNCrnkJg\nb9J8X3XMJLo9cGm8ropoeMVZ8fz4vr6PaFyD5fHn5AE/By4FziUaLLyy12esB2ri6Yqk2i4Gnsz2\nb6hX6i+1RGUwFZvZOk6OP/mreLkBf2NmrwC/JmqhTgLeBfzU3Vs9Gqos1YFjnge+bWZfIAqk7l7r\nK4EjSfN91QGw3d1fiqeXAs+6+1YAP3nP+Km+76r4tRZYA8wD5gLvBn7i8ZMJen3GA2b2GaJ/bI7b\nD0xJcf9lCFCIymBq82gM1BlEgXX8nOjHiVp3S+L1+4CiFD63m5N/V0/8OXf/W+C/AMXA82Y2r3c9\nvb7ndHW0nKmIPr7PgP/p0bngRe5+lrt//zSfcTvwP4iGb1xtZhOS9qvtTDXI0KMQlUHn7q3AF4D/\nZmYFQDmw3927zOwKopAFeBb4QHwFvRR4fx8fuY3o0B/gw8cXmtkcd1/v7n9HNCTim0LU3Q8D+WZ2\nPCj7qqO3l4BLzWxW/D3jT/N9zwCfjgcAxsxqzGwi8P+Ajx4PyV6f8Qd3vwdo5ORYuGeTY8PoSaQg\n2wXI8OTua+PD5huBHwFPmtl6oA54Ld5mjZk9DLxMdDi7qo+P+zrwfTP7JvC7pOV3xmGYADYCT5/i\nz/4SuITo8P2UdZyi9kaLHsnyWHylfj/RgxDf8n3u3mFm84EXo2EsOQb8ubtvNLNvASvNrIfocP9m\n4O/NbC5RC/Y38b4DXAH8Rx/7L0OYhsKTYc3MFgNfcvdPZLuWvpjZaGAl0SNhep/XlSFOh/MyrHn0\n1MffDuXO9kQDH9+tAM1NaomKiARQS1REJIBCVEQkgEJURCSAQlREJIBCVEQkwP8HvbROvmy+cBIA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0aa8c71160>"
      ]
     },
     "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": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "q_0_comp_w1 = Q_0(coords_lofar, coords_combined[combined_w1], field)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1 0.39366 +/- 0.00010 [0.39348 0.39380]\n",
      " 2 0.52524 +/- 0.00014 [0.52505 0.52550]\n",
      " 3 0.57438 +/- 0.00013 [0.57413 0.57458]\n",
      " 4 0.60209 +/- 0.00042 [0.60125 0.60280]\n",
      " 5 0.62362 +/- 0.00027 [0.62320 0.62401]\n",
      " 6 0.64216 +/- 0.00047 [0.64136 0.64289]\n",
      " 7 0.66026 +/- 0.00037 [0.65985 0.66092]\n",
      " 8 0.67976 +/- 0.00060 [0.67875 0.68062]\n",
      " 9 0.69942 +/- 0.00045 [0.69856 0.70000]\n",
      "10 0.71500 +/- 0.00084 [0.71368 0.71647]\n",
      "11 0.72777 +/- 0.00034 [0.72733 0.72835]\n",
      "12 0.73687 +/- 0.00060 [0.73603 0.73802]\n",
      "13 0.74425 +/- 0.00076 [0.74310 0.74531]\n",
      "14 0.75317 +/- 0.00063 [0.75214 0.75437]\n",
      "15 0.75844 +/- 0.00090 [0.75645 0.75973]\n",
      "16 0.76225 +/- 0.00069 [0.76083 0.76328]\n",
      "17 0.76460 +/- 0.00071 [0.76379 0.76585]\n",
      "18 0.76753 +/- 0.00098 [0.76559 0.76880]\n",
      "19 0.76903 +/- 0.00110 [0.76740 0.77138]\n",
      "20 0.76909 +/- 0.00088 [0.76772 0.77085]\n",
      "21 0.76947 +/- 0.00141 [0.76670 0.77236]\n",
      "22 0.76735 +/- 0.00157 [0.76504 0.77048]\n",
      "23 0.76757 +/- 0.00113 [0.76508 0.76955]\n",
      "24 0.76743 +/- 0.00232 [0.76382 0.77018]\n",
      "25 0.76172 +/- 0.00249 [0.75766 0.76667]\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": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVEAAAFBCAYAAADOoJzZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VfW57/HPk3kgDCGAMqOCDCoKkWJrHatFW+tQ6xFt\nbT09Rc+tHc7toNfb28naY6dTO+jh0pardtDWVlutWFtE0KpUAjLPApFJAiQQIIQk5Ll/rIVsUwJ7\nZyV7yvf9euW117T3ehY7+bKm32+ZuyMiIh2Tk+oCREQymUJURCQChaiISAQKURGRCBSiIiIRKERF\nRCJIWoia2UwzqzGz5e3MNzP7sZmtN7OlZjYhWbWJiHRUMvdEHwKmHGf+FcDI8Gca8N9JqElEJJKk\nhai7vwjUHmeRq4FHPDAf6G1mJyenOhGRjkmnc6KDgM0x41vCaSIiaSsv1QV0hJlNIzjkp7S0dOLo\n0aNTXJGIZJuFCxfucvd+J1ounUJ0KzAkZnxwOO2fuPsMYAZAZWWlV1VVdX11ItKtmFl1PMul0+H8\nU8At4VX6ycBed9+e6qJERI4naXuiZvYocBFQYWZbgK8B+QDuPh2YBVwJrAcagFuTVZuISEclLUTd\nfeoJ5jvw6SSVIyLSKdLpcF5EJOMoREVEIlCIiohEoBAVEYlAISoiEoFCVEQkAoWoiEgEClERkQgU\noiIiEShERUQiUIiKiESgEBURiUAhKiISgUJURCQChaiISAQKURGRCBSiIiIRKERFRCJQiIqIRKAQ\nFRGJQCEqIhKBQlREJAKFqIhIBApREZEIFKIiIhEoREVEIlCIiohEoBAVEYlAISoiEoFCVEQkAoWo\niEgEClERkQgUoiIiEShERUQiUIiKiESgEBURiUAhKiISgUJURCQChaiISAQKURGRCBSiIiIRKERF\nRCJQiIqIRKAQFRGJQCEqIhKBQlREJIKkhqiZTTGzNWa23szuOsb8Xmb2tJktMbMVZnZrMusTEUlU\n0kLUzHKBB4ArgLHAVDMb22axTwMr3X08cBHwAzMrSFaNIiKJSuae6CRgvbtvcPcm4DHg6jbLOFBm\nZgb0AGqBliTWKCKSkGSG6CBgc8z4lnBarJ8CY4BtwDLgc+7e2vaDzGyamVWZWdXOnTu7ql4RkRNK\ntwtL7wcWAwOBs4GfmlnPtgu5+wx3r3T3yn79+iW7RhGRtyUzRLcCQ2LGB4fTYt0KPOGB9cBGYHSS\n6hMRSVgyQ3QBMNLMRoQXi24EnmqzzJvApQBmNgA4HdiQxBpFRBKSl6wVuXuLmd0BPAfkAjPdfYWZ\n3R7Onw7cAzxkZssAA+50913JqlFEJFFJC1EAd58FzGozbXrM8Dbg8mTWJCISRbpdWBIRyShJ3RMV\nkczg7rQ6tLpzuNVpDccPtzr+9jTIzzVKC/PIz+2++2MKUZE05u7sP9TCnoZm6hqaqD3Q9PZwXUMz\n9QebOdh0mIPNwU9j8+F3jscMH2718DNjPj9mPbHTYpeJR1F+Dj0K8ykryqNHYfhTlEdZ+FpSkEeO\nQauD47gfDWoPwxqOvkJwUSRodwNmYFj4GoznmNGvrJBT+pUyoqIHg/sUpyTMFaIiXczd2Xuwmd0H\nmqg/2My+xhbqG4PXfY3N1B8MXxuPvu4JQ3JPQxPNh4+daGZQVhgEVHFBLkX5uRTn51BckEvvkvxw\nPJhelJ/zjoAJsykYxv5pWo4ZOWbk5kBOTjhsQYjlhuM5OUZzSysHDrWw/1AL+w61sL8xGN7f2MKW\nuoPsP9T89rQj6zI7GoIWrosj42FIHglydw+C3o9M87fnHXanqeVoW5y8HGNIeQkjKkoZUVHK8IpS\nTgmHT+pZRE5OzAZ2IoWoSAe5O7UHmnirvpGafYfYWX+Imn3BcE04vKP+EDv3H3rHH3tbOQZlRcFe\nXM/wdURFKRNLC+hdUkCfknz6lBQEP6VHh3sW55PbRcGQCdyduoZmNu7az8ZdDeHrATbuauCVN3bR\n2Hz037wwL4eX7ryY/mVFnV6HQlTkOJoPt7Jtz0GqdzdQXdvA5toGqncfoHp3MHyg6fA/vadnUR79\nexbRv6yQSSPK6V9WSL+yQip6FNKz+EhQ5tOzOI+yonxKC3LfPmyV+JkZ5aUFlJeWM3FY+TvmtbY6\nO/Y1hqEafF8VpYVdUodCVCRUU99IVXUdC6vrWPPWPqprD7BtT+Pb5xIBCvJyGNKnmGF9S5l8Sl+G\nlJcwsFcR/XsW0r+siH5lhRTl56ZwKwSCUxAn9yrm5F7FvPvUii5dl0JUuqXDrc7aHfuoqq5jUXUd\nVdW1bK49CASHfqNPKuPsIX24enwJQ/uWMLS8hGF9SxhQ1nXn1iQzKUSlW9jb0MzybXup2lTHwjfr\neL26jn3hxY6KHoVUDuvDx88bzoRhfThjYC8K8rrvLTuSGIWoZBV3Z9veRlZuq2fFtr3haz1b9wR7\nmWZw+oAyrjp7IJXD+lA5rJwh5cU6JykdphCVjOXubNh1gGVb9gaBub2eldvqqWtoBoLAHFFRyoRh\nffjo5GGMG9iT8UN606s4P8WVSzZRiErGaGppZfm2vSzcVMeCTbVUVddRe6AJCC74jD6pjPePO4lx\nA3sydmAvRp9URmmhfsWla+k3TNJWfWNzcNEnDM3Fm/dwKLzfcnjfEi4Z3Z/KYX04e2hvTu3Xo1s3\nPZTUUYhK2nB3Vm6v59llb/H86hpWv1WPe9BCZtzAntz8rmGcO7wPE4f36ZKbpkU6QiEqKeXuLNu6\nl1nL3uLZ5dup3t1AjsGkEeV87tKRnDu8nLOH9NZhuaQt/WZK0rW2Oou37OHZZduZtewttu45SF6O\n8e7TKvj3C0/lsrED6Nuja1qXiHQ2hagkhbuz6M09PLN0O88u3872vY3k5xrvHdmPz79vJJeNHUDv\nkoJUlymSMIWodKm9Dc088foWHn3tTdbu2E9BXg4XjurHl6ecziWjB+h2I8l4ClHpdMFeZx2//seb\nPLN0O4daWhk/uBff+fCZXHnmyZQVKTgleyhEpdPsPdjMk4u28Ohrm1mzYx89CvO4fuJgpk4ayhmD\neqW6PJEuoRCVSI6c63z0tTf589JtNDa3ctbgXtx33ZlcNX6grqpL1tNvuHSIu/PSul38cPZaXn9z\nD6UFuVw3YTA3aa9TuhmFqCTE3Xnljd388G9rqaquY2CvIu65ehzXThhMD+11Sjek33qJ26tv7OaH\ns9fy2sZaTupZxD3XnMENlYMpzFMnxNJ9KUTlhP6xIQjP+RtqGdCzkG98aBz/cu4Q9eAugkJUjqNq\nUy0/nL2Wl9fvpl9ZIV/94FhuetdQhadIDIWo/JOF1XXcP3stL63bRUWPAr7ygTHc/K5hFBcoPEXa\nUojK2xa9Wcf9s9fx4tqdlJcWcPeVo/no5GGUFOjXRKQ9+usQXg/Dc14Ynv/riiA8dY+nyInpr6Qb\nW7x5D/fPXsvcNUF43nXFaD6m8BRJiP5auqHFm/fwo9lreWHNTvqU5HPnlNHccp7CU6Qj9FfTjSze\nvIcfP7+OOatr6FOSz5ennM7Hzxuu8BSJQH89Wc7dmbd2J9PnvcH8DbX0DsPzlvOGq4WRSCfQX1GW\najncyjPLtjN93gZWba/n5F5FfOUDY5g6aaj2PEU6kf6asszBpsP8rmozP3tpA1vqDjKyfw++/5Hx\nfGj8QAry9DRMkc6mEM0SdQeaeOTVah5+dRO1B5qYOKwPX7tqHJeO7k9OjqW6PJGspRDNcNv2HOTn\nL23ksQVv0tB0mEtH9+f2i07l3OHlqS5NpFtQiGaojbsOMH3uGzzx+hbc4UPjB3Lbhady+kllqS5N\npFtRiGaYldvqeXDuemYt205+bg43TRrKpy44hcF9SlJdmki3pBDNEAura3nghTeYs7qGHoV5TLvg\nVD55/gj6len57CKppBBNY0cewfHAC+v5x8Za+pTk84XLRnHLecPpVaInZoqkA4VoGnJ3Zq+q4Sdz\n1rF0y15O6lnE//ngWKZOGqIelUTSjP4i08ySzXu4d9YqXttYy7C+Jdx33ZlcO2GQHsEhkqYUomli\nc20D33tuDU8t2Ubf0gK+dc0Z3HjuEPJydYO8SDo7YYia2ZXAYnffloR6up29Dc08MHc9D728iZwc\nuOPi07jtwlMoK9I5T5FMEM+e6LXAN81sALAaWAIsDl9XuvvheFdmZlOAHwG5wM/d/b5jLHMRcD+Q\nD+xy9wvj/fxM0tTSyq/mV/PjOevYe7CZD08YzBcuH8XJvYpTXZqIJOCEIerunwIws7uBQcAG4GJg\nBlALDI5nRWaWCzwAXAZsARaY2VPuvjJmmd7Ag8AUd3/TzPontjnpz935y/K3uO8vq6ne3cB7TuvL\n3VeOYdzAXqkuTUQ6IJFzov/i7uOPjJjZg8CXEnj/JGC9u28I3/8YcDWwMmaZm4An3P1NAHevSeDz\n096at/Zx95PLWFhdx6gBPfh/t57LRaP6Yaa27SKZKpEQrTezie6+EMDdF5rZqATePwjYHDO+BXhX\nm2VGAflmNhcoA37k7o8ksI605O489Mom/vPZ1ZQV5nHfdWdy/cTBumgkkgUSCdFPAk+Y2QJgIXAm\n0NwF9UwELgWKgVfNbL67r41dyMymAdMAhg4d2skldK6d+w7xxceXMG/tTi4d3Z/vXH8WFT3Uykgk\nW8Qdou6+1swmANcQBOgq4O4E1rUVGBIzPjicFmsLsNvdDwAHzOxFYDzwjhB19xkE52SprKz0BGpI\nqjmrd/Clx5ey/1AL91w9jo9OHqZDd5EsE3eImtklwM3AHmA5sBTYn8C6FgAjzWwEQXjeSHAONNaf\ngJ+aWR5QQHC4/8ME1pEWGpsP8+1Zq3jk1WpGn1TGo9MmM2qAelcSyUaJHM7PBD5PcOvRWQR7pOOA\n0+J5s7u3mNkdwHMEtzjNdPcVZnZ7OH+6u68ys78QBHQrwW1QyxOoMeVWbqvnc4+9zrqa/fzb+SP4\n0pTT1dpIJIuZe3xHw2Y2Lx3v2aysrPSqqqpUl0FrqzPz5Y189y9r6FWSzw8+Mp4LRvVLdVki0kFm\nttDdK0+0XCJ7oi+a2X8A93u8ydtN1NQ38oXHl/DSul28b8wAvnv9WZSXFqS6LBFJgkRCdCzBBaU7\nzWwhQaulxe7+eJdUliFq6hv50E9fZs/BJu699gxumjRUF49EupFErs5/GMDMijkaqJOBbhuiTS2t\n/I9fL2LvwWZ+f/u7OWOQWh2JdDeJXJ3vC9wANAIrgN+6+0NdVFdG+NYzK6mqruOnN52jABXpphJp\nMvMk0A/4NvA9YK+ZreqSqjLA7xdu4ZFXq5l2wSl88KyBqS5HRFIkkRAtc/dvAjvCq/RT6aaH8su3\n7uXuJ5fx7lP78uX3n57qckQkhRIJ0cbw9ZCZFbv7H4DLu6CmtFZ7oInbfrmQitICfjL1HLV/F+nm\nErk6/30zKwd+B8w0s1eA3l1TVnpqOdzKZx5dxM79h/j97efRV23gRbq9uHej3P0P7l7r7j8AZhG0\ng7+uyypLQ9/76xpeXr+bb11zBmcN7lb/f4hIO+IOUTO7wMxeNLMVwBXA72M7VM52zyzdzv+dt4Gb\n3zWUGyqHnPgNItItJHJCbyZwD3AR8DBwv5nd0BVFpZu1O/bxpd8v4ZyhvfnaVeNSXY6IpJFEzonu\ncve/hcPPmdnLwHyCc6RZa+/BZm775UJKCvKY/tGJFOTpQpKIHHXCRDCzR8zs88DfzeyrYTd1AIc4\nesU+K7W2Ov/zt4vZXNvAgzdPYEDPolSXJCJpJp7dql8QdEtXTvBMpPVmNpvgyZ8vdGFtKfeTOet5\nfnUNX/nAGCaNKE91OSKShuJ52uc8YN6R8XBPdDRwNkGv81np5fW7uP/5tVx3ziA+/u7hqS5HRNJU\nIudEgaBzZYKe7ZcDv+r0itLEr+ZX069HIfdee6Z6ZRKRdukqyTE0tbTy0rpdXDqmP8UF6pVeRNqn\nED2GBZtq2X+ohUtGD0h1KSKS5hIOUTO7qisKSSdzVtdQkJfDe07rm+pSRCTNdWRP9N5OryLNzFld\nw3mn9KWkIOFTxiLSzXQkRLP6KsuGnfvZuOsAl47pn+pSRCQDdCREs/ohdXNW1wBw8ekKURE5MV1Y\namPO6hpGDejBkPKSVJciIhlAIRqjvrGZ1zbW6qq8iMStIyG6o9OrSBN/X7eLllbX+VARiVvCIeru\nl3VFIeng+VU19CrO55wh6nBZROKjw/lQa6szd00NF53eT89NEpG4KS1CS7bsYfeBJi4ZrUN5EYlf\npBA1s1s7q5BUm7O6hhyDC0f1S3UpIpJBou6JfqNTqkgDz6+qoXJYOb1LClJdiohkkBO2azSzpe3N\nArLiXqC39jaycns9d10xOtWliEiGiadx+ADg/UBdm+kGvNLpFaXAkVZKOh8qIomKJ0T/DPRw98Vt\nZ5jZ3E6vKAXmrN7B4D7FjOzfI9WliEiGiSdEP+vuB441w91v6uR6kq6x+TAvr9/NDZWD1YO9iCQs\nngtLa8zsQTM7q8urSYFXN+zmYPNhLtahvIh0QDwhejrwOvBzM3vZzG4xs8Iuritp5qyqoTg/l8mn\nqANmEUncCUPU3Q+4+8/cfRLwaeA8YJWZ/cDMRnZ5hV3I3ZmzuobzR1ZQlK9nKYlI4k4YomZ2mplN\nMLMLgSHA34EHgQ8QPHs+Y63dsZ+tew7qqryIdFg8F5bWAluBJwluc9oP7APuCV8z1vOrgw6p1AGz\niHRUPCE6AbgNuBh4DHjE3bOiO7wXVtdwxqCenNSrKNWliEiGiuec6GJ3/3dgMlAD/NHMfmdml3R5\ndV2o7kATC6vruER7oSISQSKPs2wF/gS8AFwCPGhm7u5juqSyLjZv7U5aHS4ZkxUtV0UkReJpO19H\n8HC6A0B9+LMPWB4OZ6Q5q2uo6FHAWYN6pboUEclg8eyJlrt7Vj3hs+VwK3PX1HD5uJPIyVErJRHp\nuHjOiWZVgAIsrK6jvrGFS3Vrk4hElHB/omZ2k5k9Zma/NrPfmNnUBN47xczWmNl6M7vrOMuda2Yt\nZnZ9ovXFY86aGvJzjfNHVnTFx4tIN5LIhaUjLnT3G4+MmNkDwKMnepOZ5QIPAJcBW4AFZvaUu688\nxnLfAf7agdriMmdVDZNGlFNWlN9VqxCRbqIjPdsXmtkHzOwsM7sSKI7zfZOA9e6+wd2bCO45vfoY\ny30G+APB7VSdbnNtA+tq9uvZ8iLSKeLaEzWzgcClQAnwW6ACOJNgj/KOONc1CNgcM74FeFeb9QwC\nriW4sf/cOD83IUc6YNb5UBHpDPHc4nQ58DAwFzgE3E4Qpre6e2f3bH8/cKe7tx6vb08zmwZMAxg6\ndGhCK3h+dQ2nVJQyvKI0Sp0iIkB8e6LfAt7r7uuPTDCz84AZZvZJoMHdl8fxOVsJOjA5YnA4LVYl\n8FgYoBXAlWbW4u5/jF3I3WcAMwAqKyvjvnvgwKEW5r+xm1vOGxbvW0REjiueEC2IDVAAd3/VzD5M\n8OiQQwSH9ieyABhpZiMIwvNG4B0947v7iCPDZvYQ8Oe2ARpFU0srn3zvCC4fq/OhItI54gnRRjPr\n5+47Yye6+1ozO0xwrvSE3L3FzO4AngNygZnuvsLMbg/nT0+w9oT1KS3gzil6oqeIdJ54QvR7BJ2O\nfMTdtx2ZaGYVwCF3j/squrvPAma1mXbM8HT3T8T7uSIiqXLCEHX3P4SPA3nVzBYCS4AC4AaC86Ui\nIt1WXPeJuvtvgDEE50B7Ac3ATe7+cBfWJiKS9uJuseTuDcDMLqxFRCTjdKTFkoiIhBSiIiIRKERF\nRCJQiIqIRKAQFRGJQCEqIhKBQlREJAKFqIhIBApREZEIFKIiIhEoREVEIlCIiohEoBAVEYlAISoi\nEoFCVEQkAoWoiEgEClERkQgUoiIiEShERUQiUIiKiESgEBURiUAhKiISgUJURCQChaiISAQKURGR\nCBSiIiIRKERFRCJQiIqIRKAQFRGJQCEqIhKBQlREJAKFqIhIBApREZEIFKIiIhEoREVEIlCIiohE\noBAVEYlAISoiEoFCVEQkAoWoiEgEClERkQgUoiIiEShERUQiSGqImtkUM1tjZuvN7K5jzL/ZzJaa\n2TIze8XMxiezPhGRRCUtRM0sF3gAuAIYC0w1s7FtFtsIXOjuZwL3ADOSVZ+ISEckc090ErDe3Te4\nexPwGHB17ALu/oq714Wj84HBSaxPRCRhyQzRQcDmmPEt4bT2fBJ49lgzzGyamVWZWdXOnTs7sUQR\nkcSk5YUlM7uYIETvPNZ8d5/h7pXuXtmvX7/kFiciEiMvievaCgyJGR8cTnsHMzsL+DlwhbvvTlJt\nIiIdksw90QXASDMbYWYFwI3AU7ELmNlQ4AngY+6+Nom1iYh0SNL2RN29xczuAJ4DcoGZ7r7CzG4P\n508Hvgr0BR40M4AWd69MVo0iIokyd091DZFUVlZ6VVVVqssQkSxjZgvj2YlLywtLIiKZQiEqIhKB\nQlREJAKFqIhIBApREZEIFKIiIhEoREVEIlCIiohEoBAVEYlAISoiEoFCVEQkAoWoiEgEClERkQgU\noiIiEShERUQiUIiKiESgEBURiUAhKiISgUJURCQChaiISAQKURGRCBSiIiIRKERFRCJQiIqIRKAQ\nFRGJQCEqIhKBQlREJAKFqIhIBApREZEIFKIiIhEoREVEIlCIiohEoBAVEYlAISoiEoFCVEQkAoWo\niEgEClERkQgUoiIiEShERUQiUIiKiESgEBURiUAhKiISgUJURCQChaiISAQKURGRCJIaomY2xczW\nmNl6M7vrGPPNzH4czl9qZhOSWZ+ISKKSFqJmlgs8AFwBjAWmmtnYNotdAYwMf6YB/52s+kREOiKZ\ne6KTgPXuvsHdm4DHgKvbLHM18IgH5gO9zezkJNYoIpKQZIboIGBzzPiWcFqiy4iIpI28VBfQEWY2\njeBwH2C/ma0BKoBdqauqS2ibMkO2bVO2bQ90bJuGxbNQMkN0KzAkZnxwOC3RZXD3GcCM2GlmVuXu\nlZ1TanrQNmWGbNumbNse6NptSubh/AJgpJmNMLMC4EbgqTbLPAXcEl6lnwzsdfftSaxRRCQhSdsT\ndfcWM7sDeA7IBWa6+wozuz2cPx2YBVwJrAcagFuTVZ+ISEck9Zyou88iCMrYadNjhh34dAc/fsaJ\nF8k42qbMkG3blG3bA124TRbkloiIdISafYqIRJAVIXqi5qSZyMw2mdkyM1tsZlWprqcjzGymmdWY\n2fKYaeVm9jczWxe+9klljYloZ3u+bmZbw+9psZldmcoaE2VmQ8zsBTNbaWYrzOxz4fRM/p7a26Yu\n+a4y/nA+bE66FriM4Ob8BcBUd1+Z0sIiMrNNQKW7Z+z9emZ2AbCfoBXaGeG07wK17n5f+B9eH3e/\nM5V1xqud7fk6sN/dv5/K2joqbBF4srsvMrMyYCFwDfAJMvd7am+bbqALvqts2BONpzmppIC7vwjU\ntpl8NfBwOPwwwS93RmhnezKau29390Xh8D5gFUErwUz+ntrbpi6RDSGarU1FHZhtZgvDFlrZYkDM\nvb9vAQNSWUwn+UzY69jMTDrsbcvMhgPnAP8gS76nNtsEXfBdZUOIZqvz3f1sgp6tPh0eSmaV8Ja2\nzD6fFPQ0dgpwNrAd+EFqy+kYM+sB/AH4vLvXx87L1O/pGNvUJd9VNoRoXE1FM427bw1fa4AnCU5b\nZIMdR3rmCl9rUlxPJO6+w90Pu3sr8DMy8Hsys3yCsPm1uz8RTs7o7+lY29RV31U2hGg8zUkzipmV\nhifEMbNS4HJg+fHflTGeAj4eDn8c+FMKa4msTVeN15Jh35OZGfALYJW7/1fMrIz9ntrbpq76rjL+\n6jxAeKvC/RxtTnpvikuKxMxOIdj7hKBV2W8ycZvM7FHgIoIedHYAXwP+CPwOGApUAze4e0ZcrGln\ney4iODx0YBNwWyb192Bm5wMvAcuA1nDy3QTnEDP1e2pvm6bSBd9VVoSoiEiqZMPhvIhIyihERUQi\nUIiKiESgEBURiUAhKiISgUJURCQChah0GjM7HHYxttzMnjaz3gm+/+tm9sVw+Jtm9r5OqKnYzOaF\nvX2lJTM708weSnUd0jEKUelMB9397LCbuFo6/qgX3P2r7j67E2r6V+AJdz8cz8LhQxKT+nfh7suA\nwWY2NJnrlc6hEJWu8iphb1pm1sPMnjezRWFH0293VWhm/9vM1prZ34HTY6Y/ZGbXh8ObzKwiHK40\ns7nh8IUxHey+fqSpbBs3EzZZbK8OMxtuQafejxA0BRxiQUffi8xsiZk9f7z1mdmXzGxB2DvQN2K2\n4ZZw2hIz+2U47SPhnvoSM3sxps6nCZosS6Zxd/3op1N+CDq8haD57ePAlHA8D+gZDlcQPM3VgIkE\nTfNKgJ7h9C+Gyz0EXB8ObwIqwuFKYG44/DTwnnC4B5DXpp4C4K2Y8fbqGE7QPHByOK8fQfeKI8Lx\n8vbWR9CvwYzwc3KAPwMXAOMIOguvaPMZy4BB4XDvmNreAzyd6u9QP4n/aE9UOlOxmS3maP+Tfwun\nG/BtM1sKzCbYQx0AvBd40t0bPOiqLNGOY14G/svMPksQSC1t5lcAe2LG26sDoNrd54fDk4EX3X0j\ngB9tM36s9V0e/rwOLAJGAyOBS4DHPXwyQZvPeMjMPkXwn80RNcDABLdf0oBCVDrTQQ/6QB1GEFhH\nzoneTLB3NzGcvwMoSuBzWzj6u/r2+9z9PuDfgGLgZTMb3baeNus5Xh0HTlREO+sz4D89OBd8truf\n5u6/OM5n3A58haD7xoVm1jdmuw6eqAZJPwpR6XTu3gB8FviCmeUBvYAad282s4sJQhbgReCa8Ap6\nGXBVOx+5ieDQH+DDRyaa2anuvszdv0PQJeI7QtTd64BcMzsSlO3V0dZ84AIzGxGup/w463sO+New\nA2DMbJCZ9QfmAB85EpJtPuMf7v5VYCdH+8IdRYZ1oyeBvFQXINnJ3V8PD5unAr8GnjazZUAVsDpc\nZpGZ/RZYQnA4u6Cdj/sG8AszuweYGzP982EYtgIrgGeP8d6/AucTHL4fs45j1L7TgkeyPBFeqa8h\neBDiP63fWG/BAAAAiklEQVTP3Q+Z2Rjg1aAbS/YDH3X3FWZ2LzDPzA4THO5/AviemY0k2IN9Ptx2\ngIuBZ9rZfklj6gpPspqZTQD+w90/lupa2mNmhcA8gkfCtD2vK2lOh/OS1Tx46uML6XyzPUHHx3cp\nQDOT9kRFRCLQnqiISAQKURGRCBSiIiIRKERFRCJQiIqIRPD/AUD42mv2XIPwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0a9c6f5b00>"
      ]
     },
     "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 (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": 0
}