fjaviersanchez/Ivezic_depth_check.ipynb

## Ivezic_depth_check.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import astropy.table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "example_table = astropy.table.Table.read('LSST_r_512_512.fits',hdu=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def calculate_depth(band,make_plot=False, exptime=6900./30.):\n",
    "    nominal_depth0 = {'u':23.68, 'g':24.89, 'r':24.43, 'i':24.00, 'z':24.45, 'y':22.60}\n",
    "    instrumental_noise_corr = {'u':0.67 , 'g':0.21, 'r':0.14, 'i':0.08, 'z':0.05, 'y':0.04}\n",
    "    exptime_corr = {'u':0.24 , 'g':0.09, 'r':0.05, 'i':0.04, 'z':0.03, 'y':0.02} #This is the correction if exptime = 60\n",
    "    airmass_corr = {'u':0.21 , 'g':0.15, 'r':0.14, 'i':0.13, 'z':0.13, 'y':0.15}\n",
    "    exptime_corr[band] = instrumental_noise_corr[band]-1.25*np.log10(1+(10**(0.8*instrumental_noise_corr[band])-1)/exptime)\n",
    "    fid_depth = nominal_depth0[band]+instrumental_noise_corr[band]+exptime_corr[band]+airmass_corr[band]\n",
    "    source_table = astropy.table.Table.read('LSST_%s_512_512.fits' %(band))\n",
    "    if make_plot:\n",
    "        plt.hist(source_table['ab_mag'][source_table['snr_grpf']>5.],alpha=0.5,label='%s-band' %band)\n",
    "        plt.xlabel('mag$_{AB}$')\n",
    "        plt.ylabel('Counts')\n",
    "        plt.legend(loc='best')\n",
    "    maxdepth = np.nanmax(source_table['ab_mag'][source_table['snr_grpf']>5.])\n",
    "    mediandepth = np.nanmedian(source_table['ab_mag'][source_table['snr_grpf']>5.])\n",
    "    meandepth = np.nanmean(source_table['ab_mag'][source_table['snr_grpf']>5.])\n",
    "    mad = np.nanmedian(np.fabs(source_table['ab_mag'][source_table['snr_grpf']>5.]-mediandepth))\n",
    "    print 'Band: ',band, 'median 5-sigma depth = ', mediandepth, '+-', mad, ' Fiducial value = ', fid_depth\n",
    "    return maxdepth, mediandepth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Band:  u median 5-sigma depth =  25.4252 +- 0.661346  Fiducial value =  25.2242815359\n",
      "Band:  g median 5-sigma depth =  26.0608 +- 0.747641  Fiducial value =  25.4588863457\n",
      "Band:  r median 5-sigma depth =  25.4173 +- 0.716911  Fiducial value =  24.8493060544\n",
      "Band:  z median 5-sigma depth =  24.2508 +- 0.789758  Fiducial value =  24.6797723306\n",
      "Band:  y median 5-sigma depth =  23.3362 +- 0.722896  Fiducial value =  22.8298195494\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAERCAYAAACZystaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8lOWZ//HPNQkQEIFwCichKhYVq+KqRSsaW7UeFsUT\nKrZyUHetVephbaVWDLrrua2u667bn62mVqpYFWgVtS4Giy0iVgSxoaBCAiGBEEIIIZlkcv3+mCEk\nIQmTZCaThO/79eLFPKf7uW4Snmue032ZuyMiIge3QKIDEBGRxFMyEBERJQMREVEyEBERlAxERAQl\nAxERoR2SgZn90MxWR/7MjMxLNbN3zGytmb1tZn3jHYeIiDQtrsnAzMYC1wMnAycC/2xmRwJ3A++6\n+xhgMTArnnGIiEjz4n1mcAzwobtXunsIeB+4DLgYyIqskwVMinMcIiLSjHgng8+ACZHLQr2AC4HD\ngDR3LwRw9wJgcJzjEBGRZiTHs3F3zzGzR4A/AWXAJ0CosVXjGYeIiDQvrskAwN2fA54DMLP/APKA\nQjNLc/dCMxsCbG1sWzNTkhARaQV3t5as3x5PEw2K/D0SuBSYCywEpkVWmQosaGp7d++yf+67776E\nx6D+qW/qX9f70xpxPzMAXjWz/kAVcLO7l0YuHc0zsxnARmByO8QhIiJNaI/LRGc2Mq8YOCfe+xYR\nkejoDeQEysjISHQIcdWV+9eV+wbq38HIWnt9qT2YmXfk+EREOiIzw1t4A7k97hnEXHp6Ohs3bkx0\nGF3OqFGj2LBhQ6LDEJEE6JRnBpGsl4CIujb9u4p0Da05M9A9AxERUTIQERElAxERQckg4Q4//HAW\nL17c5fYlIp1Lp3yaqDGzZz9Bbm5J3NofObIf999/W9zaF5H4mP3QbHILc+PS9si0kdw/6/64tN3e\nukwyyM0tIT09M27tb9gQv7ZFJH5yC3NJn5Qel7Y3zN8Ql3YTQZeJYiwQCPDll1/WTk+fPp3Zs2c3\nu83y5csZO3YsAwYM4PrrrycYDAJQUlLCxIkTGTx4MAMGDGDixIls3ry5druzzz6b2bNnc8YZZ9Cn\nTx/OP/98iouLa5e/8MILpKenM2jQIB588MEY91REuhIlgxgza9GjvQDMnTuXP/3pT3zxxResXbuW\nf//3fwegpqaGGTNmkJeXR25uLr169eKWW26pt+3vfvc7srKy2LZtG5WVlTz++OMAfP7559x88828\n+OKL5Ofns3379nqJRESkLiWDGGvNS1u33norw4YNo1+/ftxzzz387ne/A6B///5ceuml9OjRg0MO\nOYRZs2bx/vvv19t2+vTpHHnkkfTo0YPJkyezcuVKAF599VUmTpzIN7/5Tbp168YDDzzQqkQlIgcH\nJYN2dOGFF3LooYfSp0+f2gM+wIgRI2o/jxo1ivz8fAD27NnDv/7rv5Kenk6/fv0466yzKCkpqZdw\nhgwZUvu5V69elJWVAZCfn89hhx1Wb9mAAQPi1jcR6dy6zA3kjqJXr16Ul5fXThcUFNQelN98881G\nt8nLy6v9vHHjRoYNGwbA448/zrp16/joo48YNGgQn376KSeddBLufsBv+UOHDiUnJ6d2ury8nO3b\nt7e6XyLStenMIMbGjRvH3Llzqamp4a233mLJkiUH3Obpp59m8+bNFBcX8+CDD3L11VcDUFZWRs+e\nPenTpw/FxcVkZmZGHccVV1zBH//4R/7yl79QVVXF7NmzNe6QiDQp7mcGZnY7cD1QA6wGpgOHAC8D\no4ANwGR339mW/Ywc2S+uj3+OHNkvqvWeeOIJpk6dytNPP82kSZO49NJLm13fzJgyZQrnnXceW7Zs\nYdKkSdxzzz0A3HbbbUyZMoWBAwcyfPhw7rzzThYuXFhv26Yce+yxPP3001xzzTWUl5dzxx131Lsc\nJSJSV1xHLTWzYcBS4Gh3D5rZy8CbwLHAdnd/1Mx+DKS6+92NbK9RS9uR/l2lK5p227S4vmfw/BPP\nx6Xttuio9QySgEPMrAboCWwGZgFnRZZnAdnAfslARKQj+2TlJ0y7bVpc2m7vt5vjmgzcPd/Mfgbk\nAuXAO+7+rpmluXthZJ0CMxsczzhEROJhd3B3l3m7Oa7JwMz6AZcQvjewE3jFzK4FGl6LaPLaRN2b\nphkZGapdKiLSQHZ2NtnZ2W1qI96Xic4BvnT3YgAzex04HSjce3ZgZkOArU010JInaEREDkYNvyjP\nmTOnxW3E+9HSXGC8maVY+NGXbwOfAwuBaZF1pgIL4hyHiIg0I973DJab2e+BT4CqyN+/BA4F5pnZ\nDGAjMDmecYiISPPi/jSRu88BGp6zFBO+hCQiIh2A3kAWERElg0RT2UsR6Qi6zEB18SxtB12rvJ2I\nSENdJhnEs7QddK3ydiIHs09//R6+tU1DodXqua6QlQ/PB8AG9+WEGWfHpN1E6DLJoCP529/+xg03\n3MAXX3zBd77zHQKBAF/72te4//7GzyyWL1/OrbfeSkFBAZMmTeJ//ud/6N69OyUlJXzve9/jww8/\nJBQKcfrpp/PMM88wfPhwIFz2csKECSxevJhVq1Zx+umnM3fuXPr37w+Ey17ee++97N69m9tvv73d\n+i/SkfnWndw2JLqBJw9kVUo3jo+09URBSUzaTBTdM4ixqqoqLrvsMmbMmEFxcTHXXHMNr7/+erPb\nqOyliCSakkGMLVu2jFAoxC233EJSUhKXXnopp556arPbqOyliCSakkGM5efn117G2WtvpTOVvRSR\njkr3DGJs6NCh+12OycvLY/To0Sp7KSIdls4MYuy0004jKSmJp59+mlAoxIIFC1i+fHmz26jspYgk\nWpc5MxiZNjKuj3+OTBsZ1XrdunXjtdde4/rrr2fWrFlccMEFTJw4kR49ejS6vspeikhHENeyl23V\nVcpejh8/nu9///tMnTo10aE0q7P9u4pEo2HZy5UPz4/do6X/t4rjv308EH609MS7J8WkXWhbSc3W\nlL3UZaI4eP/99yksLCQUCpGVlcXq1as5//zzEx2WiEiTusxloo5k7dq1TJ48mfLyco444gheffVV\n0tLSEh2WiEiTlAzi4MYbb+TGG29MdBgiIlFTMhCRDueJ2bMpyY3NwJObli+lJGdl7XTxui0Qo3sG\nXUlck4GZfQ14mXDBewOOAO4FXojMHwVsACa7e2xGjhKRTq8kN5fM9PSYtDU/ZyX96hz8Z6+O3+jG\nnVlcbyC7+z/cfZy7nwT8E7AbeB24G3jX3ccAi4FZ8YxDRESa155PE50DfOHuecAlQFZkfhYQu+ex\nRESkxdozGVwFzI18TnP3QgB3LwAGt2McIiLSQLvcQDazbsDFwI8jsxq+2dTkm051h2DIyMggIyMj\nxtG1rzlz5rB+/XpeeOGFLrUvEUmc7OxssrOz29RGez1NdAHwsbsXRaYLzSzN3QvNbAiwtakNox2P\nJ5ZPHzSm38iR3NZEcZqWas+hpDVstUjX1/CL8pw5c1rcRnslg2uA39WZXghMAx4BpgIL2rqDWD59\n0JjMDRtavE0oFCIpKSn2wYiIxFjc7xmYWS/CN49fqzP7EeBcM1sLfBt4ON5xtJfDDz+cRx99lBNO\nOIHevXtTU1Oz3zp79uzh6quvpk+fPpx88smsWrWqdtkjjzzC6NGj6dOnD8cddxzz58+vXZaVlcWE\nCRO466676N+/P0ceeSRvvfVW7fINGzaQkZFB3759+c53vkNRUREiItGIezJw93J3H+Tuu+rMK3b3\nc9x9jLuf5+6du3hoAy+99BKLFi2ipKSEQGD/f+KFCxdy1VVXsWPHDq655homTZpEKBQCYPTo0Xzw\nwQeUlpZy33338d3vfpfCwsLabZcvX84xxxzD9u3bueuuu7j++utrl02ZMoVTTjmFoqIifvrTn5KV\nlbXfvkVEGqOB6uLghz/8IcOGDWty2Op/+qd/4tJLLyUpKYk77riDiooKli1bBsDll19eO47RlVde\nyVFHHVWvHsKoUaOYMWMGZsbUqVPZsmULW7duJS8vjxUrVnD//ffTrVs3JkyYwMSJE+PfWRHpEpQM\n4mBv3YC5c+fWlrm86KKLapfXLUdpZowYMaK21OVvfvMbxo0bR2pqKqmpqaxZs6be5Z66ZS579uwJ\nhIvg5Ofnk5qaWjsPwolDRCQaGpsoDvY+wTNlyhSmTJmy3/K6ZS7dnU2bNjFs2DByc3P5l3/5F957\n7z1OO+00AMaNGxdVjYGhQ4eyY8cO9uzZU5sQcnNzG71MJSLSkI4UCfDxxx8zf/58QqEQv/jFL0hJ\nSWH8+PHs3r2bQCDAwIEDqamp4bnnnuOzzz6Lqs2RI0dy8sknc99991FVVcXSpUv5wx/+EOeeiEhX\n0WXODPqNHNmqxz9b0n40onmu/5JLLuHll1/muuuu46ijjuL1118nKSmJY445hjvvvJPx48eTlJTE\nddddxxlnnBH1/l588UWmTp3KgAEDOO2005g6dSolJV3q3ryIxInKXkot/btKosx+aDa5hfteGt30\nzlJu6dM7Jm1vKdzCMWcds29ff1rF/eceH5O2u1LZyy5zZiAinVduYW69OsUlDYadblPbmzVkdTR0\nz0BERJQMREREyUBERFAyEBERlAxERAQlAxERQclARERQMkiYJUuW1BuwrqvsS0Q6py7z0tnsxx4j\nN47FXEYOHMj9d90V0zZV/lJEOoq4JwMz6ws8CxwH1AAzgH8ALwOjgA3AZHff2Zb95BYVkX7TTW0L\nthkbnnkmbm2LiCRae1wmehJ4092PAU4AcoC7gXfdfQywGJjVDnG0i3nz5tXWMOjTpw8pKSl861vf\nanRdd+ehhx5i0KBBHHHEEcydO7d22ZtvvslJJ51E3759GTVqVL0C1xs3biQQCPCb3/yGUaNGMXjw\nYB588MHa5RUVFUybNo3+/ftz3HHH8dFHH8WvwyLSJcT1zMDM+gAT3H0agLtXAzvN7BLgrMhqWUA2\n4QTR6U2ePJnJkycDsGvXLsaPH99oTQOAgoICiouLyc/P569//SsXXnghp5xyCkcddRS9e/fmhRde\nYOzYsXz22Wece+65jBs3josvvrh2+w8++IB169aRk5PDqaeeyuWXX86YMWPIzMzkq6++4quvvqKs\nrIzzzz+/XfouIp1XvM8MDgeKzOw5M/ubmf3SzHoBae5eCODuBcDgOMfR7tyda665hrPPPpsbbrih\n0XXMjAceeIBu3bpx5plnctFFFzFv3jwAzjzzTMaOHQvAcccdx9VXX82SJUvqbZuZmUn37t05/vjj\nOeGEE/j0008BeOWVV/jpT39K3759GT58ODNnzoxzb0Wks4v3PYNk4CTgB+6+wsx+QfgMoOE4yU2O\nm5yZmVn7OSMjg4yMjNhHGQc/+clP2L17N08++SR5eXkce+yxQPggXlpaCkBqaiopKSm124waNaq2\n/OWHH37IrFmz+OyzzwgGgwSDQa688sp6+9hbKxmgV69elJWVAZCfn19benNvuyLSdWVnZ5Odnd2m\nNuKdDDYBee6+IjL9KuFkUGhmae5eaGZDgK1NNVA3GXQWL730Ei+//DIrVqwgKSmJww47jF27du23\nXmNlKr/+9a8DcO211zJz5kzefvttunXrxu2338727duj2v/QoUPJy8vjmGPCY7hv3LgxRj0TkY6o\n4RfluvcYoxXXy0SRS0F5Zva1yKxvA2uAhcC0yLypwIJ4xtGePvnkE2bOnMn8+fPp379/s+u6e22Z\nyj//+c+88cYbtfcbysrKSE1NpVu3bixfvrzezeW92zZl8uTJPPTQQ5SUlLBp0yb+67/+q+0dE5Eu\nrT3eM5gJvGhm3YAvgelAEjDPzGYAG4HJbd3JyIED4/r458iBA6Nab+HChZSUlHDGGWfg7pgZEyZM\n4I033thv3aFDh5KamsqwYcM45JBD+N///V+OOuooAP77v/+bO+64g1tuuYWzzjqLq666ql4Jy4bv\nDdSdvu+++7jppps4/PDDGT58ONOnT+fJJ59sTbdF5CChspdSS/+ukijTbptWr9LZyofnc1uMKp3V\nLU0JKnvZFA1HISIiSgYiIqJkICIiKBmIiAhKBiIigpKBiIigZCAiIigZiIgISgYJlZWVxYQJE7rc\nvkSk8+kyZS8fm/0YRbnxK3s5cORA7ro/tmUvQaUvRaRj6DLJoCi3iJvS41f28pkNKnspIl2XLhPF\n2OOPP84VV1xRb97MmTO5/fbbG12/pqaGW2+9lX79+nHssceyePHi2mXPP/88xx57LH369GH06NH8\n8pe/rF22ZMkSDjvsMH7+85+TlpbG8OHDef7552uXFxcXc/HFF9O3b1/Gjx/PF198EduOikiXomQQ\nY9/97nd5++23awvYhEIhXn75ZaZOndro+h9++CFHHXUU27dvJzMzk8suu6x2dNK0tDTefPNNSktL\nee6557j99ttZuXJl7bYFBQXs2rWL/Px8nn32WX7wgx+wc+dOAG6++WZ69epFYWEhv/rVr/j1r38d\n556LSGemZBBjQ4YM4cwzz+SVV14BYNGiRQwaNIgTTzyx0fXT0tKYOXMmSUlJTJ48mTFjxtQOd33B\nBReQnp4OwIQJEzjvvPP485//XLtt9+7duffee0lKSuKCCy6gd+/erF27lpqaGl577TUeeOABUlJS\nGDt2bJPJSEQElAzi4rrrruO3v/0tAC+++CLf+973WLp0KYceeih9+vSprWYGMHz48Hrb1i19uWjR\nIk477TQGDBhAamoqixYtoqho303yAQMGEAjs+xHuLX25bds2QqGQSl+KSNTingzMbIOZfWpmn5jZ\n8si8VDN7x8zWmtnbZtY33nG0p0mTJrFq1SrWrFnDH//4R6699lrOOOMMdu3aRWlpKatXr65dd/Pm\nzfW2zc3NZdiwYQSDQa644gp+9KMfsW3bNnbs2MEFF1wQVb2BQYMGkZSURF5eXr12RUSa0h5nBjVA\nhruPc/dTI/PuBt519zHAYmBWO8TRbnr06MHll1/OlClT+MY3vlHvG3pDhYWFPPXUU1RXV/PKK6+Q\nk5PDRRddRDAYJBgMMnDgQAKBAIsWLeKdd96Jav+BQIDLL7+czMxM9uzZw+eff05WVlasuiciXVCL\nHy01s1TgMHdfFe0m7J90LgHOinzOArIJJ4hWGzhyYFwf/xw4Mrqyl3tNnTqVZ599tt4TPo0ZP348\n69atY+DAgQwZMoRXX32Vfv3CFZ7+8z//kyuvvJJgMMjEiRO55JJLmm2r7nsETz31FNOnT2fo0KEc\nffTRzJgxg/fee69FfRCJl9mPPUZunUueS1etZ2VJQe108bYgD5eVNLZpkwZbiBlHDohZjAebqMpe\nmlk2cDHh5PExsBX4wN3viGLbL4ESIAT8r7s/a2Y73D21zjrF7r5f9fjOXPYyLy+PY445hoKCAnr3\n7p3ocKLSGf5dpWuY9uMfk37TvveC5r+9gH5H7ytzWbA0h/G9U1rUZsHv53P36P1LZarsZXSiPTPo\n6+6lZnYD8Bt3v8/Moj0z+Ka7bzGzQcA7ZrYWaHjE6VJHoJqaGn72s59x9dVXd5pEICIHt2iTQbKZ\nDQUmA/e0ZAfuviXy9zYzmw+cChSaWZq7F5rZEMJnGo3KzMys/ZyRkUFGRkZLdt/uysvLSUtL4/DD\nD2fRokWJDkdEDgLZ2dlkZ2e3qY1ok8Ec4G1gqbt/ZGZHAOsOtJGZ9QIC7l5mZocA50XaWghMAx4B\npgILmmqjbjLoDHr16sWuXbsSHYaIHEQaflGeM2dOi9uINhlscffai2zu/qWZ/TyK7dKA183MI/t6\n0d3fMbMVwDwzmwFsJHzGISIiCRJtMngKOCmKefW4+1fAfq/eunsxcE6U+xYRkThrNhmY2WnA6cAg\nM6v75FAfICmegYmISPs50JlBd6B3ZL1D68wvBa5odIt2MGrUKI3NHwcaskLk4NVsMnD3JcASM3ve\n3Te2U0wHtGHDhkSHICLSpUR7z6CHmf0SSK+7jbt/Kx5BiYhI+4o2GbwCPAM8S/hNYhER6UKiTQbV\n7v4/cY1EREQSJtpk8Aczuxl4HajcOzPyiKiIdGENB5WLhU/WrCE9pi1KW0WbDPaWybqrzjwHjoht\nOCLS0eQWFdUbVC4Wlt54Y0zbk7aLKhm4++HxDkRERBInqmRgZtc1Nt/dfxPbcEREJBGivUx0Sp3P\nKcC3gb8BSgYiIl1AtJeJbq07bWb9gJfiEpGIdArLsrOp2LmzVdsWb91K9oJ9gxWXbs6lomhfpbM9\n23dBC4vbSNu0uOxlxG5A9xFEDmIVO3eS0W//ymLR2NStW71tc4oKSKlz8H+vsHVJRlov2nsGf2Bf\nNbIk4BhgXryCEhGR9hXtmcHjdT5XAxvdfVMc4hERkQQIRLNSZMC6HMIjl6YCwXgGJSIi7SuqZGBm\nk4HlwJWEq5J9aGZRD2FtZgEz+5uZLYxMp5rZO2a21szeNrO+rQleRERiI6pkANwDnOLuU939OsJF\n7e9twX5+CHxeZ/pu4F13HwMsBma1oC0REYmxaJNBwN231pneHu22ZjYCuJDwiKd7XQJkRT5nAZOi\njENEROIg2hvIb5nZ28DvItNXAW9Gue0vCI9pVPdSUJq7FwK4e4GZDY6yLRERiYMD1UAeTfjAfZeZ\nXQacEVn0V+DFAzVuZhcBhe6+0swymlnVm1kmIiJxdqAzgyeIXM9399eA1wDM7OuRZRMPsP03gYvN\n7EKgJ3Comb0AFJhZmrsXmtkQYGtTDWRmZtZ+zsjIICMj4wC7FBE5uGRnZ5Odnd2mNg6UDNLcfXXD\nme6+2szSD9S4u/8E+AmAmZ0F3Onu3zOzR4FpwCOEh8de0FQbdZOBiIjsr+EX5Tlz5rS4jQMlg+be\nNe/Z4r3t8zAwz8xmABsJP64qIh3ce0uz2bk7PFRE6eZccuqMJ9QSZbvLyFmfUzu9q2wXKQOjG4to\nx45SQqGa/ebv3FlGTk7J/vsq20NOzob9pnv0SObww0e0PPgu6kDJYIWZ3eju/6/uTDO7Afi4JTuK\nvLi2JPK5GDinJduLSOLt3L2TfkeHvyNWNBhPqCWSuifVO/jvLI1+LKJQqIbk5P2/pyYl9SYlpbH5\npaSkpO83XVGxoUUxd3UHSga3Aa+b2bXsO/ifDHQHLo1nYCIi0n6aTQaRxz9PN7OzgeMis99w98Vx\nj0xERNpNtPUM3gPei3MsIiKSING+gSwiIl2YkoGIiCgZiIhI68teikgHM/uxx8gtKop5u5+sWUN6\nzFuVjkbJQKSLyC0qIv2mm2Le7tIbb4x5m9Lx6DKRiIgoGYiIiC4TichBaldZWb0xi1qjrGwPX321\nqUuMcaRkICIHpZoQ9cYsao2kpFIqK6tjE1CC6TKRiIgoGYiIiJKBiIigZCAiIsQ5GZhZDzP70Mw+\nMbPVZnZfZH6qmb1jZmvN7G0z6xvPOEREpHlxTQbuXgmc7e7jgBOBC8zsVOBu4F13HwMsBmbFMw4R\nEWle3C8TuXt55GMPwo+yOnAJkBWZnwVMinccIiLStLgnAzMLmNknQAHwJ3f/CEiLVFHD3QuAwfGO\nQ0REmhb3l87cvQYYZ2Z9CNdTHkv47KDeak1tn5mZWfs5IyODjIyMOEQpIg2t+PmzVOeXYGu2s+zf\nHgcgmJ9HSd/uAHjhTtZ3S4qqLUtN5sgLTmbHjlJCoRqCwSqKikpqlx9ouq7KYJBkvS5bT3Z2NtnZ\n2W1qo93+Sd291MyygfOBQjNLc/dCMxsCbG1qu7rJQETaT3V+CVePuIFVPdM4fsREAHIq15IyMAWA\nsmARI7pHdwj5ffFLAIRCNSQn98MC3UhO7le73AKVzU7XVVHR5OHioNXwi/KcOXNa3Ea8nyYauPdJ\nITPrCZwL/B1YCEyLrDYVWBDPOEREpHnxPjMYCmSZWYBw4nnZ3d80s2XAPDObAWwEJsc5DhERaUZc\nk4G7rwZOamR+MXBOPPctIiLR0xvIIiKiZCAiIkoGIiKCkoGIiKBkICIiKBmIiAhKBiIigpKBiIjQ\njmMTicj+Hpv9GEW5RTFpa/2Kj1m1+GYGjzwiJu0VrVkHI2LSlHQCSgYiCVSUW8RN6TfFpK2h697l\nhU0vcvXpN8SkvadW3B6TdqRzUDIQ6cI2bdhAdWVlq7bds3s3G9aurf0boLJ0J9XlZQBUlVdClKOW\nSsenn6RIF1ZdWUl6Skqrtu0ZCJCekkJp5G+AovIykiMJ4KuyipjFKYmnG8giIqJkICIiSgYiIoKS\ngYiIEP+ylyPMbLGZrTGz1WY2MzI/1czeMbO1Zvb23tKYIiKSGPE+M6gG7nD3scBpwA/M7GjgbuBd\ndx8DLAZmxTkOERFpRlyTgbsXuPvKyOcy4O+E32m8BMiKrJYFTIpnHCIi0rx2u2dgZunAicAyIM3d\nCyGcMIDB7RWHiIjsr11eOjOz3sDvgR+6e5mZeYNVGk7XyszMrP2ckZFBRkZGPEIUaVezH3uM3KIi\n1q/4mKHr3o1Jm1u2bYtJO9L5ZGdnk52d3aY24p4MzCyZcCJ4wd0XRGYXmlmauxea2RBga1Pb100G\nIl1FblER6TfdRMH6x+k34uTYtPmHP8SkHel8Gn5RnjNnTovbaI/LRL8GPnf3J+vMWwhMi3yeCixo\nuJGIiLSfuJ4ZmNk3gWuB1Wb2CeHLQT8BHgHmmdkMYCMwOZ5xiIhI8+KaDNz9AyCpicXnxHPfIiIS\nPb2BLCIiSgYiIqJkICIiqLiNiERhx45SQqEaKoJBAlQDEArVUF4eXYGbioogRUUlVAaDJOuo0yHp\nxyIiBxQK1ZCc3I9AdTWBQPiwYVZNIBBdFbVAoDvJyf2oqGjylSJJMF0mEhERJQMREVEyEBERlAxE\nRAQlAxERQU8TibTYY7Mfoyi3qE1trF/xMQXrH6dozbpwuSeRBFMyEGmhotwibkq/qU1tDF33Lv1G\nnMxTK26PUVQibaPLRCIiomQgIiJKBiIigu4ZiBzQ3nrFe8WibvGWbdvo19bARGIo3pXOfgX8M1Do\n7sdH5qUCLwOjgA3AZHffGc84RNpib73ivWJRt1j1iqWjifdloueA7zSYdzfwrruPARYDs+Icg4iI\nHEBck4HxkRvSAAAMvElEQVS7LwV2NJh9CZAV+ZwFTIpnDCIicmCJuIE82N0LAdy9ABicgBhERKSO\njnAD2ZtbmJmZWfs5IyODjIyMOIcj0r42bdhAdWVlm9vZs3s3e0K7Wf3xCmpqagAIVVXRu2d0NQca\nClYFKSouIlgVpCJYTqC6mupQJd07xGFD6srOziY7O7tNbSTip1poZmnuXmhmQ4Bmq13UTQYiXVF1\nZSXpKa07YNdVGgjQ0wMMDRjJvXsA8FVxFcm9Wvff3JICJPdKxpICBFIC4aI2u6OrbCbtq+EX5Tlz\n5rS4jfa4TGSRP3stBKZFPk8FFrRDDCIi0oy4JgMzmwv8BfiameWa2XTgYeBcM1sLfDsyLSIiCRTX\ny0TuPqWJRefEc78iItIyGo5CRESUDEREpGM8WioicVJTU0N5eeueAAqFQpSXVxAKhagOheiur45d\nmpKBSBfmDoFA6x5bNUsiEEjBLOkAbwNJV6BcLyIiOjOQg0Nr6havXLWS3RW7KSwuZfn7f6+dX5VX\nSE7lhDbFU75nT7PLd+woJRSqaVGbwWA1wVA1FRVBAlQD4PpKL1FSMpCDQmvqFi/IWUC/4/uRs3IT\nKWOOq50/L/chUga27Y3hGm/+QB8K1ZCc3LKKBxZIxjyZQKB7+G1hAG/7MBdycFAyEIlC2ZYdUBUC\noHpPkLIWnmU0FAoG+fSj5ZTuLKF4e9F+4wdVBIMEqqtb1mYoSCgU1PhB0ir6jRGJRlWIEd3D/11S\nAoHaz61VDozq04OeyUmYsd/4QQGq9327j5IlBTAC6G6vtIZuIIuIiM4MJHZac5O2KQNHDuSu++9q\n0TbfuegiCnc2XkE1+I9SLKVHi9or211G0sbNVFSHSElJhTaeDYh0ZPrtlphpzU3apjyz4ZkWb1O4\ncyeTbrih0WXLnv0jRw88qUXtFRUXkdwrmS8+WQM1uvQiXZuSgTTqidmzKcnNbdE2ny5dz5KVQxtd\ntqVwC1XVVQDUdO/Gnh49qGrm+Prqtld57a1FUe03OSmZAan9Kdm8nWXP/rHRdYo25sPAqJprkz17\nKnE/cOJwp/bt3lBo/7eE9cavtDclA2lUSW4umenpLdrmxysLOKtf449Dri0qICWtDwDLyiqo6N2d\nPice1+i6AIf86W0mnzsrqv1WFFVw9JFj+PCll/jGwNMaXeepfzwUVVtt5e5Rv/G79+1es8D+21Tv\njkN0Ik1TMoiz/Px8lixaUm9eqDrEx6s+pbKqZY8OAqQO7c8RR49udNmhvXpx2UUXkZzc8X+slZVB\nqnY1fcCrrg6xq5nleyUnB7B6tZNEpDUSdtQws/OBJwg/0fQrd38kUbHE07p169j+1nZGD9p3AK8o\nr6Dbmr70SktrUVt7grv5aNd6dpx2SqPLd//lL5x31ln0a+LbeUdSWVlFoDIZs8YP5DU1RjDY7YDt\n7NpVQk1pNR/u/JhgMMim/PxG16uorGxyGYAZ9OjRvX6MwcpWl4wU6WwS8ptuZgHgvwhXOssHPjKz\nBe6ek4h44m3AIQMYM3hM7XT57nI29tnFlh4VHD3i5KjbKS0v5quUHQweM4ZQKMSeivrXmUtCIbZu\n3UplK4urFxUV4e4kJSWxZ88eysrKWrR9dVUVwWCwdvrDnWv5Rt9wvxu7jp6U1I3wr8L+zAIkJXVv\ndNle7k5NDRBIoseAgQQDAZJ79Wq8vaSkJpcB1NRU7Hfgr6hserTP9cU5jO5/dLPxdWZflnTJ/4q1\nckpKSE10EB1Mor72nAqsc/eNAGb2EnAJ0LV/AxvI2bSiRckAYNOyj9nzb0XsLCmhsrSUQJ1v1pWb\nN3PVLxdS3cS37cZUp4ToO6o/VVVVFJcUQwAqgZpQDR/lrAagoGw3Q3ofcsC2ivND9ClaVzv9Xtlf\n+ap3AHfHa5xuxeGDa35ZBUEgqWh7k2cGHgxSU9T8Y6ruTrfqSqgmvG4UN26bUh0ZrrmuvTd2Q6EQ\nXucmbygU4u+FnzEsJb3RdrrCjd+DIRmc1jPRUXQsiUoGw4G8OtObCCcIOQDbVcXVI25gW7dCevYo\npXf3fd+ed+5aS3VlJQOOPz7q9l4q+i3jr/9nysvLyduWR4++PfhsdwU7q0KM7BdOAGt/PY9Tp08+\nYFvBucvo239s7XTK+nX0HT2W6spqQsEQPQ4NP+ffrbgMvvySYclNXyZKCQQYdoB7H+5OlQegxume\nnMy6NiQDGhnq2ay60Zu84enkxm8U68avdFK6IBpngUCA9RXrKc4rrp1XUVHB8u1rWF+0mreWPxt1\nW9U11QTLdrFz3TrKi4spLS2lKFDna2h5OaGaGqoroi9msiMYZPXOnVRUVlKwuxKrCFKKUxWqocIi\n34SrQ1TsPHCboWCI6op9N8VrqmuorqjGG3lG3wD3qia/zLvXUFMTbHyhiMScRfNMdMx3ajYeyHT3\n8yPTdwPe8CaymelNHxGRVnD3Fj1ml6hkkASsJXwDeQuwHLjG3f/e7IYiIhIXCblM5O4hM7sFeId9\nj5YqEYiIJEhCzgxERKRj6TAPwZnZr8ys0MxWNZh/q5n93cxWm9nDiYqvrRrrn5mdYGZ/NbNPzGy5\nmbXsOdMOwsxGmNliM1sT+TnNjMxPNbN3zGytmb1tZn0THWtrNNK/WyPzH438bq40s1fNrE+iY22N\npn5+dZbfaWY1ZtY/UTG2VnN96wrHlmb+77X82OLuHeIPcAZwIrCqzrwMwpeSkiPTAxMdZ4z79zZw\nXuTzBcB7iY6zlX0bApwY+dyb8P2go4FHgB9F5v8YeDjRsca4f+cAgcj8h4GHEh1rLPsXmR4BvAV8\nBfRPdKwx/Nl1iWNLI/3LAY5pzbGlw5wZuPtSYEeD2d8nfACpjqwTm8HyE6CJ/tUAe78t9wM2t2tQ\nMeLuBe6+MvK5DPg74YPIJUBWZLUsYFJiImybJvo33N3fda8tZryMcJ87nab6F1n8C6BlhSU6kGb6\n1iWOLY30LwcYRiuOLR0mGTTha8CZZrbMzN7rrJdRmnE78LiZ5QKPAtEN09mBmVk64TOgZUCauxdC\n+JcWGJy4yGKjTv8+bLBoBhDdmNsdWN3+mdnFQJ67r05oUDHS4GfX5Y4tDfrX4mNLR08GyUCqu48H\nfgTMS3A8sfZ94IfuPpLwD+/XCY6nTcysN/B7wn0qY/9ivJ36aYVG+rd3/j1AlbvPTVhwMVC3f0AI\n+AlwX91VEhFXLDTys+tSx5ZG+tfiY0tHTwZ5wGsA7v4RUGNmAxIbUkxNdff5AO7+ezrxkBxmlkz4\nl/EFd18QmV1oZmmR5UOArYmKr62a6B9mNg24EJiSoNBiopH+HQmkA5+a2VeEL4F9bGad7uyuiZ9d\nlzm2NNG/Fh9bOloyMOp/+5gPfAvAzL4GdHP37YkILEYa9m+zmZ0FYGbfBv6RkKhi49fA5+7+ZJ15\nC4Fpkc9TgQUNN+pE9utfZBj2u4CL3b11Q8V2HPX65+6fufsQdz/C3Q8nPH7YOHfvjAm9sd/NrnRs\naax/LT62dJj3DMxsLuE7/AOAQsKnpy8AzxG+DlYJ3OnuS5pqoyNron9rgf8EkoAK4GZ3/yRRMbaW\nmX0TeB9YTfhSkBO+xLCc8On3YcBGYLK7lyQqztZqon/3EP7ZdQf2HkSWufvNCQmyDZr6+bn7W3XW\n+RI42d2LG2+lY2rmd/P/CB9EO/WxpZn+ldLCY0uHSQYiIpI4He0ykYiIJICSgYiIKBmIiIiSgYiI\noGQgIiIoGYiICEoGIiKCkoGIiKBkIAcpMxsVKWzyXKT4zm/N7NtmtjQyfXJkvdfN7KNI4ZAb6mx/\nr5nlmNn7ZjbXzO5IXG9E2i4hNZBFOogjgcvd/XMzWwFc4+5nRIZuvge4FJju7iVmlgJ8ZGavRra7\nFPg60AP4G7AiMV0QiQ2dGcjB7Ct3/zzyeQ3h8WogPM7LqMjn28xsJfuK1xwFfBNY4O5VkeGC/9BY\n45GymIdEPh9pZgvM7Ewz+4WZHRanPom0ipKBHMzqjjRaU2e6BkiOjPr4LeAb7n4isBJIiabhyEiY\nwwiXJYTwGcRb7v4+sAc4pO3hi8SOkoEczJor1mJAH2CHu1ea2dHA+MiyD4CJZtYjUlTknxvZfhyw\nhH3J4Axgu5ldRLh6WE5MeiASI0oGcjDzJj7vnX4L6GZma4AHgb8CuPsKwrUaPgXeAFYBO/duaGbn\nAmWEa8/uTQZj3X2eu7/BvqQi0mFoCGuRVjCzQ9x9t5n1JDye/I3uvtLMxgNnuvujZjad8E3mZ4H/\nIFxcfiSw1d3fSVjwIo1QMhBpBTN7ETiW8L2A59390QSHJNImSgYiIqJ7BiIiomQgIiIoGYiICEoG\nIiKCkoGIiKBkICIiKBmIiAhKBiIiAvx/tr9aB/zmo0AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c95da50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for band in ['u','g','r','z','y']:\n",
    "    calculate_depth(band, make_plot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
	{
	"cells": [
	{
	"cell_type": "code",
	"execution_count": 1,
	"metadata": {
	"collapsed": false
	},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"Populating the interactive namespace from numpy and matplotlib\n"
	]
	}
	],
	"source": [
	"%pylab inline"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 2,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": [
	"import astropy.table"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 4,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": [
	"example_table = astropy.table.Table.read('LSST_r_512_512.fits',hdu=1)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 61,
	"metadata": {
	"collapsed": false
	},
	"outputs": [],
	"source": [
	"def calculate_depth(band,make_plot=False, exptime=6900./30.):\n",
	" nominal_depth0 = {'u':23.68, 'g':24.89, 'r':24.43, 'i':24.00, 'z':24.45, 'y':22.60}\n",
	" instrumental_noise_corr = {'u':0.67 , 'g':0.21, 'r':0.14, 'i':0.08, 'z':0.05, 'y':0.04}\n",
	" exptime_corr = {'u':0.24 , 'g':0.09, 'r':0.05, 'i':0.04, 'z':0.03, 'y':0.02} #This is the correction if exptime = 60\n",
	" airmass_corr = {'u':0.21 , 'g':0.15, 'r':0.14, 'i':0.13, 'z':0.13, 'y':0.15}\n",
	" exptime_corr[band] = instrumental_noise_corr[band]-1.25np.log10(1+(10(0.8instrumental_noise_corr[band])-1)/exptime)\n",
	" fid_depth = nominal_depth0[band]+instrumental_noise_corr[band]+exptime_corr[band]+airmass_corr[band]\n",
	" source_table = astropy.table.Table.read('LSST_%s_512_512.fits' %(band))\n",
	" if make_plot:\n",
	" plt.hist(source_table['ab_mag'][source_table['snr_grpf']>5.],alpha=0.5,label='%s-band' %band)\n",
	" plt.xlabel('mag$_{AB}$')\n",
	" plt.ylabel('Counts')\n",
	" plt.legend(loc='best')\n",
	" maxdepth = np.nanmax(source_table['ab_mag'][source_table['snr_grpf']>5.])\n",
	" mediandepth = np.nanmedian(source_table['ab_mag'][source_table['snr_grpf']>5.])\n",
	" meandepth = np.nanmean(source_table['ab_mag'][source_table['snr_grpf']>5.])\n",
	" mad = np.nanmedian(np.fabs(source_table['ab_mag'][source_table['snr_grpf']>5.]-mediandepth))\n",
	" print 'Band: ',band, 'median 5-sigma depth = ', mediandepth, '+-', mad, ' Fiducial value = ', fid_depth\n",
	" return maxdepth, mediandepth"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 62,
	"metadata": {
	"collapsed": false
	},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"Band: u median 5-sigma depth = 25.4252 +- 0.661346 Fiducial value = 25.2242815359\n",
	"Band: g median 5-sigma depth = 26.0608 +- 0.747641 Fiducial value = 25.4588863457\n",
	"Band: r median 5-sigma depth = 25.4173 +- 0.716911 Fiducial value = 24.8493060544\n",
	"Band: z median 5-sigma depth = 24.2508 +- 0.789758 Fiducial value = 24.6797723306\n",
	"Band: y median 5-sigma depth = 23.3362 +- 0.722896 Fiducial value = 22.8298195494\n"
	]
	},
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAERCAYAAACZystaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8lOWZ//HPNQkQEIFwCichKhYVq+KqRSsaW7UeFsUT\nKrZyUHetVephbaVWDLrrua2u667bn62mVqpYFWgVtS4Giy0iVgSxoaBCAiGBEEIIIZlkcv3+mCEk\nIQmTZCaThO/79eLFPKf7uW4Snmue032ZuyMiIge3QKIDEBGRxFMyEBERJQMREVEyEBERlAxERAQl\nAxERoR2SgZn90MxWR/7MjMxLNbN3zGytmb1tZn3jHYeIiDQtrsnAzMYC1wMnAycC/2xmRwJ3A++6\n+xhgMTArnnGIiEjz4n1mcAzwobtXunsIeB+4DLgYyIqskwVMinMcIiLSjHgng8+ACZHLQr2AC4HD\ngDR3LwRw9wJgcJzjEBGRZiTHs3F3zzGzR4A/AWXAJ0CosVXjGYeIiDQvrskAwN2fA54DMLP/APKA\nQjNLc/dCMxsCbG1sWzNTkhARaQV3t5as3x5PEw2K/D0SuBSYCywEpkVWmQosaGp7d++yf+67776E\nx6D+qW/qX9f70xpxPzMAXjWz/kAVcLO7l0YuHc0zsxnARmByO8QhIiJNaI/LRGc2Mq8YOCfe+xYR\nkejoDeQEysjISHQIcdWV+9eV+wbq38HIWnt9qT2YmXfk+EREOiIzw1t4A7k97hnEXHp6Ohs3bkx0\nGF3OqFGj2LBhQ6LDEJEE6JRnBpGsl4CIujb9u4p0Da05M9A9AxERUTIQERElAxERQckg4Q4//HAW\nL17c5fYlIp1Lp3yaqDGzZz9Bbm5J3NofObIf999/W9zaF5H4mP3QbHILc+PS9si0kdw/6/64tN3e\nukwyyM0tIT09M27tb9gQv7ZFJH5yC3NJn5Qel7Y3zN8Ql3YTQZeJYiwQCPDll1/WTk+fPp3Zs2c3\nu83y5csZO3YsAwYM4PrrrycYDAJQUlLCxIkTGTx4MAMGDGDixIls3ry5druzzz6b2bNnc8YZZ9Cn\nTx/OP/98iouLa5e/8MILpKenM2jQIB588MEY91REuhIlgxgza9GjvQDMnTuXP/3pT3zxxResXbuW\nf//3fwegpqaGGTNmkJeXR25uLr169eKWW26pt+3vfvc7srKy2LZtG5WVlTz++OMAfP7559x88828\n+OKL5Ofns3379nqJRESkLiWDGGvNS1u33norw4YNo1+/ftxzzz387ne/A6B///5ceuml9OjRg0MO\nOYRZs2bx/vvv19t2+vTpHHnkkfTo0YPJkyezcuVKAF599VUmTpzIN7/5Tbp168YDDzzQqkQlIgcH\nJYN2dOGFF3LooYfSp0+f2gM+wIgRI2o/jxo1ivz8fAD27NnDv/7rv5Kenk6/fv0466yzKCkpqZdw\nhgwZUvu5V69elJWVAZCfn89hhx1Wb9mAAQPi1jcR6dy6zA3kjqJXr16Ul5fXThcUFNQelN98881G\nt8nLy6v9vHHjRoYNGwbA448/zrp16/joo48YNGgQn376KSeddBLufsBv+UOHDiUnJ6d2ury8nO3b\nt7e6XyLStenMIMbGjRvH3Llzqamp4a233mLJkiUH3Obpp59m8+bNFBcX8+CDD3L11VcDUFZWRs+e\nPenTpw/FxcVkZmZGHccVV1zBH//4R/7yl79QVVXF7NmzNe6QiDQp7mcGZnY7cD1QA6wGpgOHAC8D\no4ANwGR339mW/Ywc2S+uj3+OHNkvqvWeeOIJpk6dytNPP82kSZO49NJLm13fzJgyZQrnnXceW7Zs\nYdKkSdxzzz0A3HbbbUyZMoWBAwcyfPhw7rzzThYuXFhv26Yce+yxPP3001xzzTWUl5dzxx131Lsc\nJSJSV1xHLTWzYcBS4Gh3D5rZy8CbwLHAdnd/1Mx+DKS6+92NbK9RS9uR/l2lK5p227S4vmfw/BPP\nx6Xttuio9QySgEPMrAboCWwGZgFnRZZnAdnAfslARKQj+2TlJ0y7bVpc2m7vt5vjmgzcPd/Mfgbk\nAuXAO+7+rpmluXthZJ0CMxsczzhEROJhd3B3l3m7Oa7JwMz6AZcQvjewE3jFzK4FGl6LaPLaRN2b\nphkZGapdKiLSQHZ2NtnZ2W1qI96Xic4BvnT3YgAzex04HSjce3ZgZkOArU010JInaEREDkYNvyjP\nmTOnxW3E+9HSXGC8maVY+NGXbwOfAwuBaZF1pgIL4hyHiIg0I973DJab2e+BT4CqyN+/BA4F5pnZ\nDGAjMDmecYiISPPi/jSRu88BGp6zFBO+hCQiIh2A3kAWERElg0RT2UsR6Qi6zEB18SxtB12rvJ2I\nSENdJhnEs7QddK3ydiIHs09//R6+tU1DodXqua6QlQ/PB8AG9+WEGWfHpN1E6DLJoCP529/+xg03\n3MAXX3zBd77zHQKBAF/72te4//7GzyyWL1/OrbfeSkFBAZMmTeJ//ud/6N69OyUlJXzve9/jww8/\nJBQKcfrpp/PMM88wfPhwIFz2csKECSxevJhVq1Zx+umnM3fuXPr37w+Ey17ee++97N69m9tvv73d\n+i/SkfnWndw2JLqBJw9kVUo3jo+09URBSUzaTBTdM4ixqqoqLrvsMmbMmEFxcTHXXHMNr7/+erPb\nqOyliCSakkGMLVu2jFAoxC233EJSUhKXXnopp556arPbqOyliCSakkGM5efn117G2WtvpTOVvRSR\njkr3DGJs6NCh+12OycvLY/To0Sp7KSIdls4MYuy0004jKSmJp59+mlAoxIIFC1i+fHmz26jspYgk\nWpc5MxiZNjKuj3+OTBsZ1XrdunXjtdde4/rrr2fWrFlccMEFTJw4kR49ejS6vspeikhHENeyl23V\nVcpejh8/nu9///tMnTo10aE0q7P9u4pEo2HZy5UPz4/do6X/t4rjv308EH609MS7J8WkXWhbSc3W\nlL3UZaI4eP/99yksLCQUCpGVlcXq1as5//zzEx2WiEiTusxloo5k7dq1TJ48mfLyco444gheffVV\n0tLSEh2WiEiTlAzi4MYbb+TGG29MdBgiIlFTMhCRDueJ2bMpyY3NwJObli+lJGdl7XTxui0Qo3sG\nXUlck4GZfQ14mXDBewOOAO4FXojMHwVsACa7e2xGjhKRTq8kN5fM9PSYtDU/ZyX96hz8Z6+O3+jG\nnVlcbyC7+z/cfZy7nwT8E7AbeB24G3jX3ccAi4FZ8YxDRESa155PE50DfOHuecAlQFZkfhYQu+ex\nRESkxdozGVwFzI18TnP3QgB3LwAGt2McIiLSQLvcQDazbsDFwI8jsxq+2dTkm051h2DIyMggIyMj\nxtG1rzlz5rB+/XpeeOGFLrUvEUmc7OxssrOz29RGez1NdAHwsbsXRaYLzSzN3QvNbAiwtakNox2P\nJ5ZPHzSm38iR3NZEcZqWas+hpDVstUjX1/CL8pw5c1rcRnslg2uA39WZXghMAx4BpgIL2rqDWD59\n0JjMDRtavE0oFCIpKSn2wYiIxFjc7xmYWS/CN49fqzP7EeBcM1sLfBt4ON5xtJfDDz+cRx99lBNO\nOIHevXtTU1Oz3zp79uzh6quvpk+fPpx88smsWrWqdtkjjzzC6NGj6dOnD8cddxzz58+vXZaVlcWE\nCRO466676N+/P0ceeSRvvfVW7fINGzaQkZFB3759+c53vkNRUREiItGIezJw93J3H+Tuu+rMK3b3\nc9x9jLuf5+6du3hoAy+99BKLFi2ipKSEQGD/f+KFCxdy1VVXsWPHDq655homTZpEKBQCYPTo0Xzw\nwQeUlpZy33338d3vfpfCwsLabZcvX84xxxzD9u3bueuuu7j++utrl02ZMoVTTjmFoqIifvrTn5KV\nlbXfvkVEGqOB6uLghz/8IcOGDWty2Op/+qd/4tJLLyUpKYk77riDiooKli1bBsDll19eO47RlVde\nyVFHHVWvHsKoUaOYMWMGZsbUqVPZsmULW7duJS8vjxUrVnD//ffTrVs3JkyYwMSJE+PfWRHpEpQM\n4mBv3YC5c+fWlrm86KKLapfXLUdpZowYMaK21OVvfvMbxo0bR2pqKqmpqaxZs6be5Z66ZS579uwJ\nhIvg5Ofnk5qaWjsPwolDRCQaGpsoDvY+wTNlyhSmTJmy3/K6ZS7dnU2bNjFs2DByc3P5l3/5F957\n7z1OO+00AMaNGxdVjYGhQ4eyY8cO9uzZU5sQcnNzG71MJSLSkI4UCfDxxx8zf/58QqEQv/jFL0hJ\nSWH8+PHs3r2bQCDAwIEDqamp4bnnnuOzzz6Lqs2RI0dy8sknc99991FVVcXSpUv5wx/+EOeeiEhX\n0WXODPqNHNmqxz9b0n40onmu/5JLLuHll1/muuuu46ijjuL1118nKSmJY445hjvvvJPx48eTlJTE\nddddxxlnnBH1/l588UWmTp3KgAEDOO2005g6dSolJV3q3ryIxInKXkot/btKosx+aDa5hfteGt30\nzlJu6dM7Jm1vKdzCMWcds29ff1rF/eceH5O2u1LZyy5zZiAinVduYW69OsUlDYadblPbmzVkdTR0\nz0BERJQMREREyUBERFAyEBERlAxERAQlAxERQclARERQMkiYJUuW1BuwrqvsS0Q6py7z0tnsxx4j\nN47FXEYOHMj9d90V0zZV/lJEOoq4JwMz6ws8CxwH1AAzgH8ALwOjgA3AZHff2Zb95BYVkX7TTW0L\nthkbnnkmbm2LiCRae1wmehJ4092PAU4AcoC7gXfdfQywGJjVDnG0i3nz5tXWMOjTpw8pKSl861vf\nanRdd+ehhx5i0KBBHHHEEcydO7d22ZtvvslJJ51E3759GTVqVL0C1xs3biQQCPCb3/yGUaNGMXjw\nYB588MHa5RUVFUybNo3+/ftz3HHH8dFHH8WvwyLSJcT1zMDM+gAT3H0agLtXAzvN7BLgrMhqWUA2\n4QTR6U2ePJnJkycDsGvXLsaPH99oTQOAgoICiouLyc/P569//SsXXnghp5xyCkcddRS9e/fmhRde\nYOzYsXz22Wece+65jBs3josvvrh2+w8++IB169aRk5PDqaeeyuWXX86YMWPIzMzkq6++4quvvqKs\nrIzzzz+/XfouIp1XvM8MDgeKzOw5M/ubmf3SzHoBae5eCODuBcDgOMfR7tyda665hrPPPpsbbrih\n0XXMjAceeIBu3bpx5plnctFFFzFv3jwAzjzzTMaOHQvAcccdx9VXX82SJUvqbZuZmUn37t05/vjj\nOeGEE/j0008BeOWVV/jpT39K3759GT58ODNnzoxzb0Wks4v3PYNk4CTgB+6+wsx+QfgMoOE4yU2O\nm5yZmVn7OSMjg4yMjNhHGQc/+clP2L17N08++SR5eXkce+yxQPggXlpaCkBqaiopKSm124waNaq2\n/OWHH37IrFmz+OyzzwgGgwSDQa688sp6+9hbKxmgV69elJWVAZCfn19benNvuyLSdWVnZ5Odnd2m\nNuKdDDYBee6+IjL9KuFkUGhmae5eaGZDgK1NNVA3GXQWL730Ei+//DIrVqwgKSmJww47jF27du23\nXmNlKr/+9a8DcO211zJz5kzefvttunXrxu2338727duj2v/QoUPJy8vjmGPCY7hv3LgxRj0TkY6o\n4RfluvcYoxXXy0SRS0F5Zva1yKxvA2uAhcC0yLypwIJ4xtGePvnkE2bOnMn8+fPp379/s+u6e22Z\nyj//+c+88cYbtfcbysrKSE1NpVu3bixfvrzezeW92zZl8uTJPPTQQ5SUlLBp0yb+67/+q+0dE5Eu\nrT3eM5gJvGhm3YAvgelAEjDPzGYAG4HJbd3JyIED4/r458iBA6Nab+HChZSUlHDGGWfg7pgZEyZM\n4I033thv3aFDh5KamsqwYcM45JBD+N///V+OOuooAP77v/+bO+64g1tuuYWzzjqLq666ql4Jy4bv\nDdSdvu+++7jppps4/PDDGT58ONOnT+fJJ59sTbdF5CChspdSS/+ukijTbptWr9LZyofnc1uMKp3V\nLU0JKnvZFA1HISIiSgYiIqJkICIiKBmIiAhKBiIigpKBiIigZCAiIigZiIgISgYJlZWVxYQJE7rc\nvkSk8+kyZS8fm/0YRbnxK3s5cORA7ro/tmUvQaUvRaRj6DLJoCi3iJvS41f28pkNKnspIl2XLhPF\n2OOPP84VV1xRb97MmTO5/fbbG12/pqaGW2+9lX79+nHssceyePHi2mXPP/88xx57LH369GH06NH8\n8pe/rF22ZMkSDjvsMH7+85+TlpbG8OHDef7552uXFxcXc/HFF9O3b1/Gjx/PF198EduOikiXomQQ\nY9/97nd5++23awvYhEIhXn75ZaZOndro+h9++CFHHXUU27dvJzMzk8suu6x2dNK0tDTefPNNSktL\nee6557j99ttZuXJl7bYFBQXs2rWL/Px8nn32WX7wgx+wc+dOAG6++WZ69epFYWEhv/rVr/j1r38d\n556LSGemZBBjQ4YM4cwzz+SVV14BYNGiRQwaNIgTTzyx0fXT0tKYOXMmSUlJTJ48mTFjxtQOd33B\nBReQnp4OwIQJEzjvvPP485//XLtt9+7duffee0lKSuKCCy6gd+/erF27lpqaGl577TUeeOABUlJS\nGDt2bJPJSEQElAzi4rrrruO3v/0tAC+++CLf+973WLp0KYceeih9+vSprWYGMHz48Hrb1i19uWjR\nIk477TQGDBhAamoqixYtoqho303yAQMGEAjs+xHuLX25bds2QqGQSl+KSNTingzMbIOZfWpmn5jZ\n8si8VDN7x8zWmtnbZtY33nG0p0mTJrFq1SrWrFnDH//4R6699lrOOOMMdu3aRWlpKatXr65dd/Pm\nzfW2zc3NZdiwYQSDQa644gp+9KMfsW3bNnbs2MEFF1wQVb2BQYMGkZSURF5eXr12RUSa0h5nBjVA\nhruPc/dTI/PuBt519zHAYmBWO8TRbnr06MHll1/OlClT+MY3vlHvG3pDhYWFPPXUU1RXV/PKK6+Q\nk5PDRRddRDAYJBgMMnDgQAKBAIsWLeKdd96Jav+BQIDLL7+czMxM9uzZw+eff05WVlasuiciXVCL\nHy01s1TgMHdfFe0m7J90LgHOinzOArIJJ4hWGzhyYFwf/xw4Mrqyl3tNnTqVZ599tt4TPo0ZP348\n69atY+DAgQwZMoRXX32Vfv3CFZ7+8z//kyuvvJJgMMjEiRO55JJLmm2r7nsETz31FNOnT2fo0KEc\nffTRzJgxg/fee69FfRCJl9mPPUZunUueS1etZ2VJQe108bYgD5eVNLZpkwZbiBlHDohZjAebqMpe\nmlk2cDHh5PExsBX4wN3viGLbL4ESIAT8r7s/a2Y73D21zjrF7r5f9fjOXPYyLy+PY445hoKCAnr3\n7p3ocKLSGf5dpWuY9uMfk37TvveC5r+9gH5H7ytzWbA0h/G9U1rUZsHv53P36P1LZarsZXSiPTPo\n6+6lZnYD8Bt3v8/Moj0z+Ka7bzGzQcA7ZrYWaHjE6VJHoJqaGn72s59x9dVXd5pEICIHt2iTQbKZ\nDQUmA/e0ZAfuviXy9zYzmw+cChSaWZq7F5rZEMJnGo3KzMys/ZyRkUFGRkZLdt/uysvLSUtL4/DD\nD2fRokWJDkdEDgLZ2dlkZ2e3qY1ok8Ec4G1gqbt/ZGZHAOsOtJGZ9QIC7l5mZocA50XaWghMAx4B\npgILmmqjbjLoDHr16sWuXbsSHYaIHEQaflGeM2dOi9uINhlscffai2zu/qWZ/TyK7dKA183MI/t6\n0d3fMbMVwDwzmwFsJHzGISIiCRJtMngKOCmKefW4+1fAfq/eunsxcE6U+xYRkThrNhmY2WnA6cAg\nM6v75FAfICmegYmISPs50JlBd6B3ZL1D68wvBa5odIt2MGrUKI3NHwcaskLk4NVsMnD3JcASM3ve\n3Te2U0wHtGHDhkSHICLSpUR7z6CHmf0SSK+7jbt/Kx5BiYhI+4o2GbwCPAM8S/hNYhER6UKiTQbV\n7v4/cY1EREQSJtpk8Aczuxl4HajcOzPyiKiIdGENB5WLhU/WrCE9pi1KW0WbDPaWybqrzjwHjoht\nOCLS0eQWFdUbVC4Wlt54Y0zbk7aLKhm4++HxDkRERBInqmRgZtc1Nt/dfxPbcEREJBGivUx0Sp3P\nKcC3gb8BSgYiIl1AtJeJbq07bWb9gJfiEpGIdArLsrOp2LmzVdsWb91K9oJ9gxWXbs6lomhfpbM9\n23dBC4vbSNu0uOxlxG5A9xFEDmIVO3eS0W//ymLR2NStW71tc4oKSKlz8H+vsHVJRlov2nsGf2Bf\nNbIk4BhgXryCEhGR9hXtmcHjdT5XAxvdfVMc4hERkQQIRLNSZMC6HMIjl6YCwXgGJSIi7SuqZGBm\nk4HlwJWEq5J9aGZRD2FtZgEz+5uZLYxMp5rZO2a21szeNrO+rQleRERiI6pkANwDnOLuU939OsJF\n7e9twX5+CHxeZ/pu4F13HwMsBma1oC0REYmxaJNBwN231pneHu22ZjYCuJDwiKd7XQJkRT5nAZOi\njENEROIg2hvIb5nZ28DvItNXAW9Gue0vCI9pVPdSUJq7FwK4e4GZDY6yLRERiYMD1UAeTfjAfZeZ\nXQacEVn0V+DFAzVuZhcBhe6+0swymlnVm1kmIiJxdqAzgyeIXM9399eA1wDM7OuRZRMPsP03gYvN\n7EKgJ3Comb0AFJhZmrsXmtkQYGtTDWRmZtZ+zsjIICMj4wC7FBE5uGRnZ5Odnd2mNg6UDNLcfXXD\nme6+2szSD9S4u/8E+AmAmZ0F3Onu3zOzR4FpwCOEh8de0FQbdZOBiIjsr+EX5Tlz5rS4jQMlg+be\nNe/Z4r3t8zAwz8xmABsJP64qIh3ce0uz2bk7PFRE6eZccuqMJ9QSZbvLyFmfUzu9q2wXKQOjG4to\nx45SQqGa/ebv3FlGTk7J/vsq20NOzob9pnv0SObww0e0PPgu6kDJYIWZ3eju/6/uTDO7Afi4JTuK\nvLi2JPK5GDinJduLSOLt3L2TfkeHvyNWNBhPqCWSuifVO/jvLI1+LKJQqIbk5P2/pyYl9SYlpbH5\npaSkpO83XVGxoUUxd3UHSga3Aa+b2bXsO/ifDHQHLo1nYCIi0n6aTQaRxz9PN7OzgeMis99w98Vx\nj0xERNpNtPUM3gPei3MsIiKSING+gSwiIl2YkoGIiCgZiIhI68teikgHM/uxx8gtKop5u5+sWUN6\nzFuVjkbJQKSLyC0qIv2mm2Le7tIbb4x5m9Lx6DKRiIgoGYiIiC4TichBaldZWb0xi1qjrGwPX321\nqUuMcaRkICIHpZoQ9cYsao2kpFIqK6tjE1CC6TKRiIgoGYiIiJKBiIigZCAiIsQ5GZhZDzP70Mw+\nMbPVZnZfZH6qmb1jZmvN7G0z6xvPOEREpHlxTQbuXgmc7e7jgBOBC8zsVOBu4F13HwMsBmbFMw4R\nEWle3C8TuXt55GMPwo+yOnAJkBWZnwVMinccIiLStLgnAzMLmNknQAHwJ3f/CEiLVFHD3QuAwfGO\nQ0REmhb3l87cvQYYZ2Z9CNdTHkv47KDeak1tn5mZWfs5IyODjIyMOEQpIg2t+PmzVOeXYGu2s+zf\nHgcgmJ9HSd/uAHjhTtZ3S4qqLUtN5sgLTmbHjlJCoRqCwSqKikpqlx9ouq7KYJBkvS5bT3Z2NtnZ\n2W1qo93+Sd291MyygfOBQjNLc/dCMxsCbG1qu7rJQETaT3V+CVePuIFVPdM4fsREAHIq15IyMAWA\nsmARI7pHdwj5ffFLAIRCNSQn98MC3UhO7le73AKVzU7XVVHR5OHioNXwi/KcOXNa3Ea8nyYauPdJ\nITPrCZwL/B1YCEyLrDYVWBDPOEREpHnxPjMYCmSZWYBw4nnZ3d80s2XAPDObAWwEJsc5DhERaUZc\nk4G7rwZOamR+MXBOPPctIiLR0xvIIiKiZCAiIkoGIiKCkoGIiKBkICIiKBmIiAhKBiIigpKBiIjQ\njmMTicj+Hpv9GEW5RTFpa/2Kj1m1+GYGjzwiJu0VrVkHI2LSlHQCSgYiCVSUW8RN6TfFpK2h697l\nhU0vcvXpN8SkvadW3B6TdqRzUDIQ6cI2bdhAdWVlq7bds3s3G9aurf0boLJ0J9XlZQBUlVdClKOW\nSsenn6RIF1ZdWUl6Skqrtu0ZCJCekkJp5G+AovIykiMJ4KuyipjFKYmnG8giIqJkICIiSgYiIoKS\ngYiIEP+ylyPMbLGZrTGz1WY2MzI/1czeMbO1Zvb23tKYIiKSGPE+M6gG7nD3scBpwA/M7GjgbuBd\ndx8DLAZmxTkOERFpRlyTgbsXuPvKyOcy4O+E32m8BMiKrJYFTIpnHCIi0rx2u2dgZunAicAyIM3d\nCyGcMIDB7RWHiIjsr11eOjOz3sDvgR+6e5mZeYNVGk7XyszMrP2ckZFBRkZGPEIUaVezH3uM3KIi\n1q/4mKHr3o1Jm1u2bYtJO9L5ZGdnk52d3aY24p4MzCyZcCJ4wd0XRGYXmlmauxea2RBga1Pb100G\nIl1FblER6TfdRMH6x+k34uTYtPmHP8SkHel8Gn5RnjNnTovbaI/LRL8GPnf3J+vMWwhMi3yeCixo\nuJGIiLSfuJ4ZmNk3gWuB1Wb2CeHLQT8BHgHmmdkMYCMwOZ5xiIhI8+KaDNz9AyCpicXnxHPfIiIS\nPb2BLCIiSgYiIqJkICIiqLiNiERhx45SQqEaKoJBAlQDEArVUF4eXYGbioogRUUlVAaDJOuo0yHp\nxyIiBxQK1ZCc3I9AdTWBQPiwYVZNIBBdFbVAoDvJyf2oqGjylSJJMF0mEhERJQMREVEyEBERlAxE\nRAQlAxERQU8TibTYY7Mfoyi3qE1trF/xMQXrH6dozbpwuSeRBFMyEGmhotwibkq/qU1tDF33Lv1G\nnMxTK26PUVQibaPLRCIiomQgIiJKBiIigu4ZiBzQ3nrFe8WibvGWbdvo19bARGIo3pXOfgX8M1Do\n7sdH5qUCLwOjgA3AZHffGc84RNpib73ivWJRt1j1iqWjifdloueA7zSYdzfwrruPARYDs+Icg4iI\nHEBck4HxkRvSAAAMvElEQVS7LwV2NJh9CZAV+ZwFTIpnDCIicmCJuIE82N0LAdy9ABicgBhERKSO\njnAD2ZtbmJmZWfs5IyODjIyMOIcj0r42bdhAdWVlm9vZs3s3e0K7Wf3xCmpqagAIVVXRu2d0NQca\nClYFKSouIlgVpCJYTqC6mupQJd07xGFD6srOziY7O7tNbSTip1poZmnuXmhmQ4Bmq13UTQYiXVF1\nZSXpKa07YNdVGgjQ0wMMDRjJvXsA8FVxFcm9Wvff3JICJPdKxpICBFIC4aI2u6OrbCbtq+EX5Tlz\n5rS4jfa4TGSRP3stBKZFPk8FFrRDDCIi0oy4JgMzmwv8BfiameWa2XTgYeBcM1sLfDsyLSIiCRTX\ny0TuPqWJRefEc78iItIyGo5CRESUDEREpGM8WioicVJTU0N5eeueAAqFQpSXVxAKhagOheiur45d\nmpKBSBfmDoFA6x5bNUsiEEjBLOkAbwNJV6BcLyIiOjOQg0Nr6havXLWS3RW7KSwuZfn7f6+dX5VX\nSE7lhDbFU75nT7PLd+woJRSqaVGbwWA1wVA1FRVBAlQD4PpKL1FSMpCDQmvqFi/IWUC/4/uRs3IT\nKWOOq50/L/chUga27Y3hGm/+QB8K1ZCc3LKKBxZIxjyZQKB7+G1hAG/7MBdycFAyEIlC2ZYdUBUC\noHpPkLIWnmU0FAoG+fSj5ZTuLKF4e9F+4wdVBIMEqqtb1mYoSCgU1PhB0ir6jRGJRlWIEd3D/11S\nAoHaz61VDozq04OeyUmYsd/4QQGq9327j5IlBTAC6G6vtIZuIIuIiM4MJHZac5O2KQNHDuSu++9q\n0TbfuegiCnc2XkE1+I9SLKVHi9or211G0sbNVFSHSElJhTaeDYh0ZPrtlphpzU3apjyz4ZkWb1O4\ncyeTbrih0WXLnv0jRw88qUXtFRUXkdwrmS8+WQM1uvQiXZuSgTTqidmzKcnNbdE2ny5dz5KVQxtd\ntqVwC1XVVQDUdO/Gnh49qGrm+Prqtld57a1FUe03OSmZAan9Kdm8nWXP/rHRdYo25sPAqJprkz17\nKnE/cOJwp/bt3lBo/7eE9cavtDclA2lUSW4umenpLdrmxysLOKtf449Dri0qICWtDwDLyiqo6N2d\nPice1+i6AIf86W0mnzsrqv1WFFVw9JFj+PCll/jGwNMaXeepfzwUVVtt5e5Rv/G79+1es8D+21Tv\njkN0Ik1TMoiz/Px8lixaUm9eqDrEx6s+pbKqZY8OAqQO7c8RR49udNmhvXpx2UUXkZzc8X+slZVB\nqnY1fcCrrg6xq5nleyUnB7B6tZNEpDUSdtQws/OBJwg/0fQrd38kUbHE07p169j+1nZGD9p3AK8o\nr6Dbmr70SktrUVt7grv5aNd6dpx2SqPLd//lL5x31ln0a+LbeUdSWVlFoDIZs8YP5DU1RjDY7YDt\n7NpVQk1pNR/u/JhgMMim/PxG16uorGxyGYAZ9OjRvX6MwcpWl4wU6WwS8ptuZgHgvwhXOssHPjKz\nBe6ek4h44m3AIQMYM3hM7XT57nI29tnFlh4VHD3i5KjbKS0v5quUHQweM4ZQKMSeivrXmUtCIbZu\n3UplK4urFxUV4e4kJSWxZ88eysrKWrR9dVUVwWCwdvrDnWv5Rt9wvxu7jp6U1I3wr8L+zAIkJXVv\ndNle7k5NDRBIoseAgQQDAZJ79Wq8vaSkJpcB1NRU7Hfgr6hserTP9cU5jO5/dLPxdWZflnTJ/4q1\nckpKSE10EB1Mor72nAqsc/eNAGb2EnAJ0LV/AxvI2bSiRckAYNOyj9nzb0XsLCmhsrSUQJ1v1pWb\nN3PVLxdS3cS37cZUp4ToO6o/VVVVFJcUQwAqgZpQDR/lrAagoGw3Q3ofcsC2ivND9ClaVzv9Xtlf\n+ap3AHfHa5xuxeGDa35ZBUEgqWh7k2cGHgxSU9T8Y6ruTrfqSqgmvG4UN26bUh0ZrrmuvTd2Q6EQ\nXucmbygU4u+FnzEsJb3RdrrCjd+DIRmc1jPRUXQsiUoGw4G8OtObCCcIOQDbVcXVI25gW7dCevYo\npXf3fd+ed+5aS3VlJQOOPz7q9l4q+i3jr/9nysvLyduWR4++PfhsdwU7q0KM7BdOAGt/PY9Tp08+\nYFvBucvo239s7XTK+nX0HT2W6spqQsEQPQ4NP+ffrbgMvvySYclNXyZKCQQYdoB7H+5OlQegxume\nnMy6NiQDGhnq2ay60Zu84enkxm8U68avdFK6IBpngUCA9RXrKc4rrp1XUVHB8u1rWF+0mreWPxt1\nW9U11QTLdrFz3TrKi4spLS2lKFDna2h5OaGaGqoroi9msiMYZPXOnVRUVlKwuxKrCFKKUxWqocIi\n34SrQ1TsPHCboWCI6op9N8VrqmuorqjGG3lG3wD3qia/zLvXUFMTbHyhiMScRfNMdMx3ajYeyHT3\n8yPTdwPe8CaymelNHxGRVnD3Fj1ml6hkkASsJXwDeQuwHLjG3f/e7IYiIhIXCblM5O4hM7sFeId9\nj5YqEYiIJEhCzgxERKRj6TAPwZnZr8ys0MxWNZh/q5n93cxWm9nDiYqvrRrrn5mdYGZ/NbNPzGy5\nmbXsOdMOwsxGmNliM1sT+TnNjMxPNbN3zGytmb1tZn0THWtrNNK/WyPzH438bq40s1fNrE+iY22N\npn5+dZbfaWY1ZtY/UTG2VnN96wrHlmb+77X82OLuHeIPcAZwIrCqzrwMwpeSkiPTAxMdZ4z79zZw\nXuTzBcB7iY6zlX0bApwY+dyb8P2go4FHgB9F5v8YeDjRsca4f+cAgcj8h4GHEh1rLPsXmR4BvAV8\nBfRPdKwx/Nl1iWNLI/3LAY5pzbGlw5wZuPtSYEeD2d8nfACpjqwTm8HyE6CJ/tUAe78t9wM2t2tQ\nMeLuBe6+MvK5DPg74YPIJUBWZLUsYFJiImybJvo33N3fda8tZryMcJ87nab6F1n8C6BlhSU6kGb6\n1iWOLY30LwcYRiuOLR0mGTTha8CZZrbMzN7rrJdRmnE78LiZ5QKPAtEN09mBmVk64TOgZUCauxdC\n+JcWGJy4yGKjTv8+bLBoBhDdmNsdWN3+mdnFQJ67r05oUDHS4GfX5Y4tDfrX4mNLR08GyUCqu48H\nfgTMS3A8sfZ94IfuPpLwD+/XCY6nTcysN/B7wn0qY/9ivJ36aYVG+rd3/j1AlbvPTVhwMVC3f0AI\n+AlwX91VEhFXLDTys+tSx5ZG+tfiY0tHTwZ5wGsA7v4RUGNmAxIbUkxNdff5AO7+ezrxkBxmlkz4\nl/EFd18QmV1oZmmR5UOArYmKr62a6B9mNg24EJiSoNBiopH+HQmkA5+a2VeEL4F9bGad7uyuiZ9d\nlzm2NNG/Fh9bOloyMOp/+5gPfAvAzL4GdHP37YkILEYa9m+zmZ0FYGbfBv6RkKhi49fA5+7+ZJ15\nC4Fpkc9TgQUNN+pE9utfZBj2u4CL3b11Q8V2HPX65+6fufsQdz/C3Q8nPH7YOHfvjAm9sd/NrnRs\naax/LT62dJj3DMxsLuE7/AOAQsKnpy8AzxG+DlYJ3OnuS5pqoyNron9rgf8EkoAK4GZ3/yRRMbaW\nmX0TeB9YTfhSkBO+xLCc8On3YcBGYLK7lyQqztZqon/3EP7ZdQf2HkSWufvNCQmyDZr6+bn7W3XW\n+RI42d2LG2+lY2rmd/P/CB9EO/WxpZn+ldLCY0uHSQYiIpI4He0ykYiIJICSgYiIKBmIiIiSgYiI\noGQgIiIoGYiICEoGIiKCkoGIiKBkIAcpMxsVKWzyXKT4zm/N7NtmtjQyfXJkvdfN7KNI4ZAb6mx/\nr5nlmNn7ZjbXzO5IXG9E2i4hNZBFOogjgcvd/XMzWwFc4+5nRIZuvge4FJju7iVmlgJ8ZGavRra7\nFPg60AP4G7AiMV0QiQ2dGcjB7Ct3/zzyeQ3h8WogPM7LqMjn28xsJfuK1xwFfBNY4O5VkeGC/9BY\n45GymIdEPh9pZgvM7Ewz+4WZHRanPom0ipKBHMzqjjRaU2e6BkiOjPr4LeAb7n4isBJIiabhyEiY\nwwiXJYTwGcRb7v4+sAc4pO3hi8SOkoEczJor1mJAH2CHu1ea2dHA+MiyD4CJZtYjUlTknxvZfhyw\nhH3J4Axgu5ldRLh6WE5MeiASI0oGcjDzJj7vnX4L6GZma4AHgb8CuPsKwrUaPgXeAFYBO/duaGbn\nAmWEa8/uTQZj3X2eu7/BvqQi0mFoCGuRVjCzQ9x9t5n1JDye/I3uvtLMxgNnuvujZjad8E3mZ4H/\nIFxcfiSw1d3fSVjwIo1QMhBpBTN7ETiW8L2A59390QSHJNImSgYiIqJ7BiIiomQgIiIoGYiICEoG\nIiKCkoGIiKBkICIiKBmIiAhKBiIiAvx/tr9aB/zmo0AAAAAASUVORK5CYII=\n",
	"text/plain": [
	"<matplotlib.figure.Figure at 0x10c95da50>"
	]
	},
	"metadata": {},
	"output_type": "display_data"
	}
	],
	"source": [
	"for band in ['u','g','r','z','y']:\n",
	" calculate_depth(band, make_plot=True)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": null,
	"metadata": {
	"collapsed": true
	},
	"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.12"
	}
	},
	"nbformat": 4,
	"nbformat_minor": 0
	}