Background computation using ofilter

bg_ofilter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Background Computation with \"ofilter\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook illustrate background calculations using `ofilter` algorithm adapted from IRAF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "import ofiltsky\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set the seed for reproducibility\n",
    "np.random.seed(0)\n",
    "\n",
    "# Random Poisson data\n",
    "data = np.random.poisson(lam=1, size=(50, 50))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f348254cd68>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAEACAYAAAAN5psFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+0XVV17z8zgjYVchPuNUEI2HovAuUlqE8hPPuKKFUC\nDGF01PDqGxZ0iB0dtjoevFZ0lOeo0lHtECu2thq0ito+QWwl4g+QovVVCUUFkwEBvZcfAiGBQH6I\nig258/1xTpJz5ppn73nPOffmHO78jHFG7t537bXXXmfflb2++7vmFFUlSZJkmFhwoBuQJEkyU3Lg\nSpJk6MiBK0mSoSMHriRJho4cuJIkGTpy4EqSZOjIgStJkllFREZE5AsisklE7hSRk50yHxGRH4vI\nHSLy4ro6D5qdpiZJkuzjCuCrqvp6ETkI+NXWX4rIamBcVY9pDmofA1ZVVZhPXEmSzBoisgj476r6\nKQBVfVpVd5li5wCfaf7+VmBERJZV1ZsDV5Iks8mvA9tE5FMi8gMRWSsiC02ZI4EHW7Yfbu7rSA5c\nSZLMJgcBLwU+qqovBX4OXNKPSmcVEcnFkElygFBV6eX4XzvqYH3goaejxbeq6uFm30PAg6r6veb2\ntcA7TZmHgaNatpc393VkTsT50+V3mdI7GZcT0FUri98/vqJNq2Ns48+dWtrHv7M+8e2ixFUfWd22\nfdA524oyT39prP1cG35WnurWjQDVbV75XNO89vaNXbm+rNfU443oob6wC+NNe12ca4hg21h33QCj\n5tofv7DUWff2T2WbT17RfqrmdbYy9aH2F1QTF93aVT22jVX9PvXgNxk/6rR9/V5F0X8VfdGpvUDo\nXJab9NoZH2N54KGn2f3IeKjswc+fKnQpVd0qIg+KyItU9UfAq4G7TLF1wNuAq0VkFbBDVbdWnSvf\nKiZJUskene61ircD/ygiBwP3Am8SkT8AVFXXqupXReRMEZkEfga8qa7CHLiSJKlk2p0bxFHVHwIv\nN7s/bsr80UzqnLOBawnPA0DWbyh+N+rMqur4yn9ZXNaDqWitc+Cq9qnYzkvLKcHI2Y1/97YZp82s\naH/kt9Mjj2LaYKaFAKN26upMETpNPxZsHmPbEROAMwX24q4Fph9WIBnFTDmdeotjnOn43qOW8DyU\nDlMoc5wn1oxf84u27W2BqZg4U7FROzWs6K8lOg2bN6KmnmIa7Zy7OE+HNhftY+bn6he7dc+s1NsL\nczZwHSZL5+pUfWPY2nxoc9AaJpYMWR/D8N0XvdLrE9dskFPFJEkq2ZMDV5Ikw0Y+cSVJMnTsGcC8\nFDLbyTJERE+X360u5PlWLFZIDhwT8euMbiyF44hlVqVdKhbTj65Hywiqnmgtps1ePbuub/fVLLqs\nFGoLITsgzluxGco27zi+vR7XN9UF3nWK9b0511D0qSN+2+/GfSkR8bnV/K1491vxIiDwsiWC9a9B\n+V3cpNf2bEAVEX34oeeHyh65/JGezxcln7iSJKkkNa4kSYaO3YM3buXAlSRJNXtcB92BZU4GrlaD\nndVIoDQResbMMbpYb4ZjNLRmRKe9kfVlRfsiBtQT2+uZOq+8zgnbZuc6H/tR+3rLRetLnWnSruG7\nutR+atch4mlGZZsLulkb6PWf0ZSs/gelMdMrE8KYjCPa4sjZU+0FHI1wm+lT754MDQvmuuzfDDhG\n1rW9r1UEmM4nriRJho15+8SVJMnwkgNXkiRDx/TcOBxmxJz4uF59ynv3bbsLi42O4mktRbwrb+Gu\n1HdwEe+qIh7Xvnq9eqxGY+rxjtl2Yvt17TyuLLN4U/u2e52mfa5Ot8poP+tnHjsKSh+SPXdI/wtc\nQ2SR9UwWnM+0nrp6AR5/66rKQt1+V13FaQuc61/75OP6/gNH1RcE/usLHkwfV5Ikg8FufVZPx4vI\n/cBOYBrYraondSj3cuC7wHmq+s9VdebAlSRJJX3QuKaBV6rq9k4FRGQB8H7ghkiFOXAlSVLJHu05\np45Qn5jnj2nEo7cBB10yy0+SJJVMsyD0qUCBb4jIbSJyof2liBwBnKuqf0/Q1jYnT1ytwqIXKdQu\npvUWH1txvltT46jOfNGrJ/gW13FyvRHSmg/H1nZnlows3BX7/QcWpbt3jF1MburxFjUXL2Dc9rXj\nRkk1x01d7iwsvrjdfOu1Z9IYfQuTr9NGry9GN3hJXPZjzaYQNE6b/vKuocB5EVXcF+v7Y0DtNFW8\nY/3P+eH6QFvhFar6iIg8j8YAtklV/73l9x+mPfNP7eCVU8UkSSrpNFVccfIhrDj5kH3bn7niCbec\nqj7S/PcxEfkX4CSgdeB6GfB5ERFgDFgtIrtVdV2nNuXAlSRJJdM9iPMi8qvAAlV9UkSeC7wG+PPW\nMqr6wpbynwK+XDVoQQ5cSZLU8J/a0zCxDPiXZmLog4B/VNUbW9OTmfIhY+mcDFytc3Yvw4pt6dSa\nhUUZu0jYuzq7gHvMKVME6ltVtmdqTbvmMH61E/DP7ohkzLG6RCQYorOvCJ7naDae0dditRSv3y02\nWJ2rBdlMQIHv3M1aY9aOj19cLiYvMu041z1+TcCAGtINzbnNthdAYMwkiHL7K7Dwv8iE5QU+DBiw\nu6FGeK9EVe8DXuzs/7hTHFV9c6TefOJKkqSSPQO45CcHriRJKtkzgK6pHLiSJKlkuncDat/JgStJ\nkkrm7RNX2wzZM86ZbTe6oxVvV5am0PGLqk2hQGk0dNTvcZu5xhFCbbRJmy1o6vVOdFMrLkcyyzjs\nON7uCUQudaKA2BcV454x00SViKSLHzPn9qJ2REyXkRT3lrJvYMfxxoBqTchBihcK5hqK+w/K79i5\nJyMveooXAyZiKzhG6VvK5nRDr4usZ4N84kqSpJI+rFXsOzlwJUlSSS8G1NkiB64kSSqZv09crfpK\nYFGzl6W3MIFGsgZ7hryA0dASitq6ol1/WXx3WY+9zl2XlrrOyPvaz1VkYaaDlmIx2opr8LSalhMl\ntcz8XR+tdlvA/GqJ/J/uLdC37eumb1w8LcroSnbBvnsN3WTR7iazO8yaAXXeivNJkgwvgxhzPgeu\nJEkqGeonrmZo1e8BD6nq60RkCXA18ALgfmCNqu6clVYmSXLAGEQ7xEyG0ncAd7VsXwLcpKrHAjcD\n7+pnw5IkGQymdUHoM5eEnrhEZDlwJvAXwEXN3ecApzZ/vgr4Fo3BrKRVSAwY8AoDKH5UCaeh7dte\n6jUranrtMWVCUVsDhk+Lm8IsImw75ts6vGuw+FE/21+KFH0T0MIjRuAI3j3gRss11L1ggPIFjBfp\nwUZOtdE+IpFGFjvfXRGpN2Kc9vpibZ8cp4Zek2XY2Zr53SLgc8DRwLOAy1X103V1RqeKfw38CTDS\nsm+Zqm4FUNUtIrI0WFeSJENEH56m9s7WFjm/extwZ1N+GgPuEZHPqerTVRXWDlwichawVVXvEJFX\nVhTt+IgxpXfu+3nJzkM4bOTX606bJMkM2a6Psp3H+l5vLz6uDrO1VhQ4tPnzocDjdYMWxJ64XgG8\nTkTOBBYCh4rIZ4EtIrJMVbeKyOHAo50qGJcT9m/koJUks8ISWcoS9k987tNNFaXj9Oic92Zrrfwt\nsE5ENgOHAOdFKq0duFT13cC7AUTkVOBiVX2jiPwVcAHwAeB84LrICUORQt2d9dFD7cJTNx27NV0G\nTKpexh6rYRWRVcsjSq3lLE8gmnT2Gay51DPI2qiaXn+ZNrv9ZfWgSN9YIkbgCAHzppuRyVyDZ+q1\n/TXqnF5r+t3TERevWFVbpqBL4/SsZfnp8MR1/22P8cBt2zoe58zWvD/t1wK3q+qrRGScRhaglar6\nZFWbevFxvR+4RkTeDDwArOmhriRJBpROdogjX3Y4R77s8H3b/+9j99gi3mztM6r6+y1l3gT8JYCq\nTonIfcBxNMT8jsxo4FLVfwP+rfnzE8DpMzk+SZLho1vnfIfZ2u+bYg/QGEe+IyLLgBcB99bVnc75\nJEkq6SVZhofJ8HMZ8GkR2avz/GnzoaiSHLiSJKmkH8kyzGzt4y37H6Ghc82IgRi4IinvLYVpDyc9\nfCCtuytIBwyBNqqnrWfMOXeRDsyLgnGNTRnmRFI1qdpcg2xA/La3o2vmLNK51YvhRf850Toj33kR\nicIrY8+1oixjozq40SGMebSbfo+kHnP7y/a7018F3vc7W+nJcpF1kiTDRibLSJJk6Oh1yc9skANX\nkiSVPD09eNEhBmLgKrQWz6TXRVTIbqNqdmP2Cy1iNtqKl81ITKu9MlbLiOhZVtuDUmNbvKnssTJK\nqtFfHE3JLjbe9WdOpNezp9q2vYXPEQ2uKOPUE7l3Iv1uM02NmlgAXgYkq1957bP1jpUxBmKRe7vM\nGlVHxpxPkmTo6MdbxX6TA1eSJJWkOJ8kydCRdogkSYaO1Lgo06pDmeLKM29GIpcWZknn/IVR1BPV\nbWTLSJutiO2lYTPCtmtANdFfI4bFiPA+4Qi+hQAdiBjbzUuJ3evKa9h24fNq65m013C1I7wbvMgP\ndSZagKk1C9u2vSi89kVE0RcB4d3tr5WntB/TRZSO2SSfuJIkGTrSDpEkydCRU8UkSYaO+TtVbJmP\nRxaihvAMeDY6Z+Bcrg5mIlt6OthonaYVyB7kZjOy5+nWjBs5JqC5eVlpWim0PRzNzenkSIRWV2eq\nwdP7IrrmuM3g49RdGGDtuRztbMfxZttdWF+fqchmf3J14Mji7C7odeCqyfLzBuCdzc2fAn+oqrVh\ndfOJK0mSSvrwxFWV5ede4LdUdaeInAFcCdSGhxk8Z1mSJAPFtEro49GS5ecT3u9Vdb2q7mxurgeO\njLQpn7iSJKnk6d6c83VZflp5C/C1SKU5cCVJUkm3U8Vglp+9ZU+jkTjjNyN1z8nA1SqYRlKGuYJ5\nNyZVRywtRGBHzLXCsWsIrHuhEEl75h0XiXARiA4xfnG7sD15eb3wbqOvAqUptZs0bF5fmQgSO90I\nEqYpXor717eL1hMXl4K+G9nVYKNDRCKM2vvCu84R08fiDAKRCK3W3GqNrQDYFylr+5OerNPAtf2O\nB9l+x0NVh0ay/CAiK4G1wBmquj3SpnziSpKkkk4D18iJRzNy4tH7tu//jHljH8jyIyJHA18E3qiq\n7bGOKsiBK0mSSrTPPi6T5edS4DDg70REgN2qelJdHTlwJUlSST+c8xVZfi4ELpxpfXMycNVl7Sky\n5ngp5c0cf/K8ssy4kTciqemnHO3H6iQRjaTQRDwzYCSKZQCrK005fRHSrwyRBcr2GjyNy2pBkew8\n3kJsTm6/Lllf9tc47XpQZEF8oQUR0CwpI9hG7oudx7X3kGeqLTIeBe7bHWvqF+j3i/nrnE+SZGjZ\nMz14ds8cuJIkqaTfGlc/yIErSZJK5u1UsU0LCCyCdTMzFxlyAllhAhmex79Q1lNoDo5npvCemevy\nFixbxqU7X5L1/VidB2DXpe3XJRsWFmVsmz2/2qhpo83w7GUGGl17S7HPUnxXa8vvygb88/rUBhf0\nFllbTavbhes2EGRx3zoS08jd9XrfjuONDnZ1fTYjd4F+obP2x8flSZ8HmnziSpKkkozHlSTJ0JEa\nV5IkQ8e81biSJBlepqfn6cDVKm57ETOLaJheJUYhtNFEPbzFx/Y/Dz89vGmfYzS0InApupaKps2q\nI47oaRdHd7uweNFZda0pxe7Fm8oykdT0xTF9MtpaI+b4RfWiercRdq2h2Xvp4BlX6849Zl46eOK8\nfdHk3ttdLNDvFzlVTJJk6BjEqWKtJVZEniMit4rI7SKyUUTe09y/RERuFJF7ROQGEYkECkuSZMhQ\njX3mktqBS1V/CZymqi8BXgysFpGTgEuAm1T1WOBm4F2z2tIkSQ4IqhL6zCWhqaKq7p2EP6d5jALn\nAKc2918FfIvGYFbQarAbvTKQhSUSJNCjC20ltLg3oDnYReCjgeZGgsG5GWnMLs84ak2NXhbowsQY\nyA4e0a+scTSSjdvr44g2te3E9j48+NLxosyis9vDPLk6k6efWex12cXQXmYnsx0xV/crS3u/Agn2\nMiiJyHOAbwPPpjF2XKuqf+6UeyWNMM8HA4+p6mlV9YYGrmZ6oe8D48BHVfU2EVmmqlsBVHWLiCyd\nwfUkSTIk9DILVNVfishpqvpzEXkW8B0R+Zqq/sfeMk2Z6aPAa1T1YREZq6s3+sQ1DbxERBYB/yIi\nJ1BezwAuDEiSpFe0RztEhxlbK28AvqiqDzfLb6urc0ZvFVV1l4h8CzgD2Lr3qUtEDgce7XTc9q/f\nuP9n/SVL8uEsSfrOTzdP8uTmcPTjML3qV96MzRR5EXCwiHwTOAT4iKp+tqrO2oGr+di2u5mwcSHw\n28D7gXXABcAHgPOB6zrVseSM1+z/+cbZCXaWJPOdQ4+Y4NAjJvZtb/3BN/pSb69vDM2M7Usi8huq\neldLkYOAlwKvAp4L3CIit6jqZKc6I09czweuao6aC4CrVfWrIrIeuEZE3gw8AKzpVEGradEzhRaC\npZPlpIi86Yi5NhJoXfr4RkXOvogQWpN1KJIK3o2gGYjYYAV7T/AdtTsiGY8iRPrGieDZDbZ9buQH\n04BFlznnttfpZGCy382O48orKyLjauBlkOl3t8/NyBC6d7wXF15mqT7Q6YnrF3feyy/uum8G9eiu\n5lPVGTQyW+/lIWCbqj4FPCUi3wZOBLofuFR1I43R0O5/Ajg93OokSYaTDgPXwt8YZ+Fv7H+Lu+OL\n3yzKVMzYWrkO+JumeP8c4GTgQ1VNSud8kiSV9DhV7DRj25fpR1XvFpEbgA3AHmCtmUoW5MCVJEk1\nPQxcFTO2j5vtDwIfjNY7N4usazSPUbPweptT3uo6bvZmY6iMZJfxKDIUR7L8GNzFvgHTpV1UffC5\nzpvhDWbBbcBoG9LpIpFBA5m2i/4KaD9udp6AuXT0snZjrY3QCjBhFoaP/vDJsj2mjaOOJmjvizFj\nOnaNo162J0sXC8W9e/ugcx5r3/HxokhX9GqHmA3yiStJkkoyOkSSJMPHAFrLc+BKkqSGfOJKkmTY\nyCeuWMowr0woZZgx5XuievF/R8AQW6ze9+ruxrTqRIO1aeYXva8+UkAoqoPzTrvo50A9kQgc3fz/\n7ArSJnroNkpBempNIMW9jeoQSDfnieqjNVE53Kgmti3OvtB9YV8eeOe6sqNfszdy4EqSZNjIt4pJ\nkgwf+cSVJMnQkXYIXB3FRgK1WhWUxkwb4RNgzBr5nEWn1ugY0cq8NQ9Wcyg0OM/8ag2LAWyWHXCu\n4cTyGqxZc/d19bpOZHFvXR1RIpqlTWnv9cXI3aaemvZCzITsUavFBoIDRMyl7sJ6AqZj+7e1vj8R\nUL1sVAeafOJKkqSaHLiSJBk6cqqYJMnQkU9cSZIMHdMHugElczJwtYqPrjhpRNhI1M9xHCHUoI6J\nsIh0GRD5QwSOEdMea4wEmPrQybX1FCZLJ9KrTcllU64B9anHKF8oFJE+vKibXQRwinwPo2tvKQ+M\npKSLRCE1uC8qzAuEHefVf1eLN7X3j/eCJvJCwbYncp/Q3TsI5+Rzkp7sI8Bq4GfABap6R1W9+cSV\nJEklvbxVDKYnWw2Mq+oxInIy8DFwlkm0UJvJOkmSeY4GP50Or09Pdg7wmWbZW4EREVlW1aQcuJIk\nmVVEZIGI3A5sAb7hpCc7EniwZfvh5r6OzMlUsU1TCCzKjRgErV4EhLQMqzGMOpKINcR6FIZAqdcy\nIjpYRBMpjJCBMiECWWJ2HG/3lNpZJFpnqVkGFop3swic8r/3yfPKNtt+94y1tp9V6rMZWe1u7Mr6\nxdseIe3TRHq9v77aEJ2mir/48SRPTdbncQykJ5sxqXElSVJNpyw/E8ewcOKYfds7vl6dx7EiPdnD\nwFEt28ub+zqSU8UkSaqZDn4cRGRMREaaP+9NT2YWarEO+P1mmVXADlXdWtWkfOJKkqSSHtcqRtKT\nfVVEzhSRSRp2iDfVVZoDV5Ik1cxNerI/mkm9czJwtYqh407kByveupEkAxFGp9YsbNuecKJhhtK6\nrwikJ6sRhSMpw9zIqhttuq2FRRl7XV49aoR2L0VcEQXDeeFh1Y3xgJAceXlQRA/tIjWai/fyx9Tj\nRUktXiA4LypsFNmxtfXRIXYeZ74/5+VBxN4ZSb3nvrDqB7nkJ0mSYSPD2iRJMnxkdIgkSYaO+frE\nZc19dXi609Tl7Qa8iYtLncLqAJHF2t2mOy/0oYD+UkT9dHQna1j09Ji6er26QzpTIEtNXeRXcCLa\nOvrVqFkkb6O6ghMV1VlMXtTrnKvIouN0qe2fbc53U2Da4517xJQJmZu7NB3PWgTU+RodIkmS4SU1\nriRJho8cuJIkGTpy4EqSZNjIqSIdIjZYwdwpM35xvfAeiSphvwMr+kMp/LsG1Jp6XUHf1OPVO3pl\nu/jt1bPjuPazeS8qRrXeEGujErgivz1/NwJ5IFKFl56seFFhoh9A+RLCfdlihOyd148XZUYuM8dE\nTMcWz1xqOtWNgGquwYvca6NDTFztiPMRg+4zhHziSpKkmgF84qqNDiEiy0XkZhG5U0Q2isjbm/uX\niMiNInKPiNywdwV4kiTPLGQ69plLImFtngYuUtUTgFOAt4nIccAlwE2qeixwM/Cu2WtmkiQHjB5D\nN88GtVNFVd1CI+QqqvqkiGyiEejrHODUZrGrgG/RGMyqCczDPZOeNQ26elafon7O2gJgcy7XlmuN\nhgEdzNX7rKEyoLm5C85tvbY9gYi2kcw7HpGFxQWuztT+V1VkQHKO8zRBa5K1upxrCg0sdi803rKW\nkBF5tuhFnBeR5TTiyS+jEbXrSlX9SIeyLwe+C5ynqv9cVe+MNC4R+TXgxTQSHy3bG+xLVbeIyNKZ\n1JUkyZDQ29PU3hnbHSJyCPB9EblRVduCCTbjdb0fuCFSaXjgap70WuAdzScvezkdL++R7+1vywL9\nJYflGJckfecJfZTtPNb3entMT+bN2I6kjIL6xzTGl5dH6g0NXCJyULPSz6rqdc3dW0VkmapuFZHD\ngUc7Hf/8l71238+H3d6vLJVJkrRymCzlMPY/FNynm/pTcZ/0q5YZ261m/xHAuap6moicFKkrGnP+\nH4C7VPWKln3rgAuaP58PXGcPSpJk+OnHW0U7YzO//jDwztbidW2qfeISkVcA/xPY2MyNpsC7gQ8A\n14jIm4EHgDWd6mgVdCMmRzcdu21XIFXVqJQCcBkN04km8NZT2rYjUSa6wY3ialJneVFcI/8BWkE6\nlPItkM4tEuejTsQOE3jZUqYMK1tYRMpwUo8VES0cLbxIQWdeVLgmX61/iRPp01B8lVmKDtHphvvZ\nTyb52U8maw/vMGNr5WXA56XRwWPAahHZrarrOtUZeav4HeBZHX59em2rkyQZbjoMXM89aoLnHjWx\nb3vbd27sVIM3Y9tfveoL9/4sIp8Cvlw1aEE655MkqaFHO0SnGdsLaGb5MYeEzpYDV5Ik1fT2VrFq\nxuaVf3Ok3NwMXC3z+l3OAtfRs9v1l1Eva01AN/HSplusljHq+fqsPhSINmn1DU9TKoytTr3FQuJA\nVhjvXMVicieFe5F23qmn0L0i2YtqtKBwPdaM6+l05l7xDJ6FCdS5d2zdk05/2e+m0COd+6+baO0R\nHdjVWGdpkXVGh0iSZPjIgStJkmEjn7iSJBk+5u3AtWr/3P+xH5UaxAhm0as6PWU1pMCa026DFlpN\nK7IwtvCiHVdew/jFAS+V3eEtRjb94+ko9uxepqVCu3P6wvqiQgHtDF15kCgXfbvBGc01uHpkgMgC\nbnsfjF8TOFck+1NAvz343G3tx2ws+2tqjblP/1efsvzM24ErSZLhJQeuJEmGjhy4kiQZNnKqmCTJ\n8JEDF4wd+3i5s4vIpZHooV5Ez8LI16WYW2eOHPWOMduRtPPegvOJPplULZFIqoXI70UKtabQQETb\nfpkn7eJyt27nhYe9zoPOLU2zIxe1v0TqxpzrfleRF09GjJf15bnGzVd8f31rQsx1PPkI+cSVJEkl\nOVVMkmT4yIErSZKhYwAHLlFPE+jnCUT09FPeu2/bXUBqjH2ermPxjI/WLOkaR42WMbVmYW3dnlZW\np1e5JsJz2uOBjzjZZgrd6z/Kc++6fqJte+S95TV0oxl5vS5GD7L3ixc8L1Kv5fG3lvra+X/8tbbt\n6y88tShT1ON855Ggj6Uu5yzWrrmQ0Y2lLmZNoV427qnXt39/XmbybrhJr0VVu1nnvQ8R0Zde+KFQ\n2R9ceVFxPhH5JHA2sFVVC3FRRBYBnwOOphFF4nJV/XTduaKhm5Mkma/0llfxU8BrO/4W3gbcqaov\nBk4DLm9GTK0kp4pJklTivqkNoqr/LiIvqCoCHNr8+VDgcVV9uq7eHLiSJKlklu0QfwusE5HNwCHA\neZGDcqqYJEk1vU0V63gtcLuqHgG8BPhoMyNQJXPzxLV+w74fPaVw1GinnnmzNcIEEMqWYuv1yow7\nOujOL7dHaR0522mPFdHNiwHvBcP4ZUbwDZgRvfvhsXvae2j3yrJXbTRY96VIwOxqDafWyOpmQLLZ\nb5xzW0bXll/WVaxur3d9faTXHWvKyKVjAaNokRXJaY81LxeRRoojSlOox8jK9nq9PrXRRvol4Efo\n5OPatWWSn24pXzDNkDcBfwmgqlMich9wHPC9qoNyqpgkSTUdBq5FyyZYtGz/G+7NP+yY5Ufo/P/X\nAzSyhX1HRJYBLwLurWtSDlxJklTSY5affwJeCYyKyE+A9wDPZn+Gn8uAT4vI3mnZn6rqE3X15sCV\nJEk1vWX5eUPN7x+h2i7hMicDl5cxpZUJE0nya9d+tiizenn9eYrMNpeX5y20gZNKnWn3unbtR08u\nzYhWH7ImR0/biETZtNqUq5s42a0tRUYhx4RpM33vOL40shbRYO01OJmji6zQgew8biRaozNFFqWP\nX9NdRiFbxtVirXnZ3tfOH7i937wxYKfRr7zvV99ar4P1I7u6R65VTJJk6JDpwRu5cuBKkqSawRu3\ncuBKkqSajMeVJMnwMYBPXHMSHeLVLdEhdl5aRnV4+ktjbdtjzir7bSuM4BtJ6+4orGdt2N62/RUn\n4kAkAoJD+VadAAAOl0lEQVQVivsVJcFGY4ikagudKyBse3EE6oRZLx2YjdDgRU3YdWl7lITd5h6A\n+hRwjQaaRjv9VdwrXUYqtcd185fjXUPkOiNp16yJ9qb17+lLdIj/9rsfDJX97rX/u+fzRcknriRJ\nqpnlh5tuyIErSZJKUuNKkmTomLc+LmlZZO0tWH78Le36hkaW5brZZdp1ianXl8bRr6y0BtTyWyn0\nqsBi6CLDkHMJZTTMUrewGqAb3dQSyJLkRWTtNl19KxGjpnfjLzp7sm3b3gNA7UJ2oJjGeNmMiuim\nrobUroMV6eyBCesLDdxv1oA6+omyfbuub1/UP3pZ/ffSbaTXrsipYpIkw8a8feJKkmSIGcCBqzaQ\noIh8UkS2tqzeRkSWiMiNInKPiNwgIiOz28wkSQ4UorHPXBKJgOoFu78EuElVjwVuBt7V74YlSTIg\nTGvsM4fUThU7BLs/B9jr3LwK+BaNwcynRWT1jJrdGA0jkTfHnSiR9jhPzLXnd0XPOkHcSZFuI116\nAvDi64yILk4aNisKe0KyiTDgpXUvpG7nmmxqtiIKqCPw25Rvti1QRoMYjRiKVzqieiTaasvLIYBR\nrTeXLl5RH8nDPmV491vRp879v3vd80y9ZV/Y72HHGqd99m9i7bVlmS7oxQ4RSE/2BuCdzc2fAn+o\nqrUO4W5jzi9V1a0AqroFWNplPUmSDDqqsY9PXXqye4HfUtUTaQQVvDLSpH6J85XPiVMPfnPfz0t0\nmsMkx7kk6Tc/3TzJk5t7jgFf0It+VZeeTFVbpzPrgSMj9XY7cG0VkWWqulVEDgcerSo8ftRp+zc2\nzzzDcpIk9Rx6xASHHrE/BvzWH3yjPxXPnXz1FuBrtaWID1w22P064ALgA8D5wHVVB2uLPuXpV1ab\n8nSnCEXUT0dyiGRmsbpNJPKm1cG8iJ5jG6weNPM07+BoK17ETLPt9rvtZy/7jdWiIpltqDfEqtEs\nXUMl1efudP7iXOYaPDNucS5Pc+sCGw12x3llVN7xq9vPZfUsKKOtjtztnMwz6PaBTglhtz9xLzu2\n1+a1iJ1D5DQaGX9+M1K+duDqEOz+/cAXROTNNLJ0rOm2wUmSDDgdxPkli1/IksUv3Ld9/33/2lX1\nIrISWAucoarb68pD7K1ip2D3p8+gbUmSDCmdnrhmUgUdHo5F5Gjgi8AbVTUs0KVzPkmSanrwaAXS\nk10KHAb8nYgIsFtVT6qrNweuJEkq6fGtYl16sguBC2da79xEh2j52TMsFumsItnFvXRWXYjqHmOB\nqAm1LxAcodQ+co85ad4jhsVuoh3sON6JDBo4l8VelXdPP26j1XqGT1uPM5GYPK9dsJ9wRP86gyw4\n35VjZLX1FKnHvPPblGbePWnOveND5bntPemZhUfubu8gbzDxIs32hYwOkSTJsJGBBJMkGT7yiStJ\nkqFj8MatORq4Wha5uobPbvymnlnSbFstyCNkdg1EGJ00mogX3XTX/2nPbDPy3vrFvh5WK3ONmabN\nnkm1IJLZJhD1c/zi9vZE0sV7OtjY2tlZZbHjeGff5eb766K/rKkWSs3Nuy/sE433N+LqoYYiQ1Sf\n6IMdou/kE1eSJNXsyYErSZIhI5+4kiQZPnLgSpJk6Ji3A1eriBkQ1SNiuIut20S+BMBGtlxVnmvS\niShqsamqrJjrCayL3tcunm5bWZ7Hmg89I+TiTe21e4bFiHmz6C/vxUDxXbSfe/HdjtG2JqKEW29A\n2LZRRICiza7R1myPX+PUYyLWet9fXZQJz7hcpGpz6u2GULSPfpE+riRJho3UuJIkGT5y4EqSZOiY\nHry54pwMXK2LqN0F1NbI52gtEb0jEvUzFsEzYOQzbS7q8bS8gHHUYvUsKHUT7zqt0dEuZAfnu/C0\nRbGLe/uT8v7gc7eZg8r/1UfPbtcoPSOr1fe8/rKEoh04fWEzVNlze9Xuun68bXvkbCfklDWOeuc2\n2yHdcH1/svz0qnGJyBnAh2kk5/mkqn7AKfNK4K+Bg4HHVPU0W6aVfOJKkqSSXjQuEVkA/C3wamAz\ncJuIXKeqd7eUGQE+CrxGVR8WEfs+paDb9GRJkswXektPdhLwY1V9QFV3A5+nkZe1lTcAX1TVhxun\nU/M4XpIDV5Ik1fSWyfpI4MGW7YcoU5C9CDhMRL4pIreJyBvrmjQnU8WZZkxxF2IH9KCQDmCYvNwJ\nGGczEnuLV63vKPI4bfWigF8tsoA6ko3b1aLMtl0EDjDy3vas1JHgeZGJxaL3LawtU3yfXoDHwg9W\nFvnaQ3e0ba9eXnvqmKfNnMvNTL7OFDrZy+xkesw7t7kH3XtnlrL8dLq3H//5T3jiFz/pxxkOAl4K\nvAp4LnCLiNyiqpNVByRJknSmw8A1uvAoRhcetW976onvesUeBo5u2V7e3NfKQ8A2VX0KeEpEvg2c\nCHQcuHKqmCRJNXumYx+f24AJEXmBiDwb+B808rK2ch3wmyLyLBH5VeBkYFNVk/KJK0mSarR7P4Sq\n7hGRPwJuZL8dYpOI/AHNTD+qereI3ABsAPYAa1X1rqp6c+BKkqSaHp3zqvp14Fiz7+Nm+4PAB6N1\nis6ynV9E9PRT3rt/h3e+QNTPwqTqFCkWvQbEea+eQuLsIrKkOjVb8Xbiakdstn0REWGdPg2luDf9\nY6ObgmPWtBFRy9aVhlNPNLZGVkd499pssS8dHn/rKU6p9nMd/LrHihKFMTQQDdbeF+5fko1uGrjX\nI33qLSafMAv9b9JrUfXyJ8URET3jqHeEyn79wSt6Pl+UfOJKkqSaXKuYJMnQkQNXkiRDx549B7oF\nBXMzcLWO2BE9y6vC6CRTa0oD4/jVAaNrwLBoAwmOX1MaM+v0F09fG6/JhOy1z80cY87tBbArzJsz\nNAHvZWdNZqJuzcIWL2CiXSjufedjRpc7/+1fLcpcdcXqtu1FzkLnYoG+tyjdfn8mUGUsy0759HLW\nxh1t219ZsaQoY7+/0StnJwOSSz5xJUkydOTAlSTJ0NF5HeIBIweuJEkq0R4MqLNFDlxJklQzb5+4\nWgXmQDQBVxg1ArQrmFsDqpPZJiTym7o9Ib4QuwPG0VDkzRqTY6Oi9mvYFjFqOgL+pBHEJ5w+XfS+\n9ronz2vvL2t6BIrMSbv+rDy3NXx636cNYGczKXl85S2/VewbXW9MqoEMOd7LDBuZ16unqNfUo87r\njKv+pv3lwRhORJAuTNp9i4CaGleSJEPHANoheooOISJniMjdIvIjEXlnvxqVJMngoNPToc9c0vXA\n1RJL+rXACcDvichx/WpYkiQDQm+hm2eFXqaK+2JJA4jI3ljSd9uCrRqMZ5YscEyXVsfxsvwUxkfH\nHXn2hu1t29dfeGp9exwKzaEL4ygrHa3FXlcgG3c30WHByWbk3Hz2OsftqQJanmf4LI+pX3zvRv0s\nFn2XV2rNrRGjckRTivR7JOPR6AaTkcnT4CIG4i7N3bX0KM4Hs/x8BFgN/Ay4QFXvsGVa6WWqGIkl\nvY+fbu4YzHBg2b7zvgPdhBmxXR890E2YMU9kmwcfnY59HCIzMxFZDYyr6jHAHwAfq2vSnEVAfXJz\n/f+6g8b2XUM2cFGGaxl0ss2Dj05r6NOBSJafc4DPAKjqrcCIiCyralMvU8VILGkAHvneDTy5eYpH\nvncDC55azmEjv97DaZMk8XhCH52dQbU3A6o3MzuppszDzX1bO1Xay8C1L5Y08AiNWNK/5xV81e+8\nkY3/+s+sePXvsHjyqeL3dqxeerSTAcZoIIt/6dSz+8m2bS9+3dJfac81ecyxh3tNBmDHzw/hmGMP\nZ+nRv1L8bvFLzeB77PPb2+Kd26mnqNdcl72mRj1L27ZHmm3ZsXkzxxzh/6fgaVxqrt31mZnz2yJi\nrtvDu4a97dnXZq8ee1ygjL0mgIVL2/t9+bHlvSNOGzuxt80h5ce0x7sHRv7zqdoy3v1esLtdy/v7\n71eGbQ+jA2iHQFW7/gBnAPcAPwYu6VBG85Of/ByYTy9/382/3/tncL4tzvGrgK+3bF8CvNOU+Rhw\nXsv23cCyqnbNeujmJEnmLyLyLBoPN6+mMTP7D+D3VHVTS5kzgbep6lkisgr4sKpWLktI53ySJLNG\nMMvPV0XkTBGZpGGHeFNdvfnElSTJ0DEndohhWBokIp8Uka0isqFl3xIRuVFE7hGRG0Rk5EC2sRUR\nWS4iN4vInSKyUUTe3tw/kG0WkeeIyK0icnuzve9p7h/I9rYiIgtE5Acisq65PfBtfqYz6wPXEC0N\n+hSNNrZyCXCTqh4L3Ay8a85b1ZmngYtU9QTgFOBtzX4dyDar6i+B01T1JcCLgdUichID2l7DO4DW\nBKXD0OZnNHPxxBUxoB1wVPXfge1m9znAVc2frwLOndNGVaCqW/Yui1DVJ2mkLF/OYLd579qW59DQ\nV5UBbi80nmyBM4FPtOwe6DbPB+Zi4JrR0qABY6mqboXGQAEsrSl/QBCRX6PxFLOexmvkgWxzc8p1\nO7AF+Iaq3sYAt7fJXwN/QmOQ3cugt/kZz5wt+XmGMHBvMkTkEOBa4B3NJy/bxoFps6pON6eKy4GT\nROQEBri9InIWsLX5ZFuVoXlg2jxfmIuBK7w0aADZunfNlIgcDgzU6loROYjGoPVZVb2uuXug2wyg\nqruAb9EwMA9ye18BvE5E7gX+L/AqEfkssGWA2zwvmIuBa9/SIBF5No2lQevm4LzdILT/z7oOuKD5\n8/nAdfaAA8w/AHep6hUt+wayzSIytvftm4gsBH6bhi43kO0FUNV3q+rRqvpCGvftzar6RuDLDGib\n5wtz4uNqxuO5gv0GtPfP+klniIj8E/BKYJTG4s73AF8CvgAcBTwArFHVHZ3qmEtE5BXAt4GN7F9y\n8W4azuRrGLA2i8gKGkL2gubnalX9CxE5jAFsr0VETgUuVtXXDUubn8mkATVJkqEjxfkkSYaOHLiS\nJBk6cuBKkmToyIErSZKhIweuJEmGjhy4kiQZOnLgSpJk6MiBK0mSoeP/A0bem5IzMIeoAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3484c12a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(data, cmap='viridis', interpolation='none')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD3FJREFUeJzt3H/sXXV9x/Hnq9QWC6UaJmVpQZikrCxbFJPqwpZd4kSY\nGZD9QVhM0CnLEnCYLVlG/WNtl0XnH9vwj2GyiVCMypCosIVANexrosloVaBoK3QS6I/R6hKH1NJA\n53t/fE/hhrX9ftve7/fc3s/zkdz03M89534/h9LnPd9z77mpKiRJbVjQ9wQkSfPH6EtSQ4y+JDXE\n6EtSQ4y+JDXE6EtSQ2aMfpKVSR5J8oMkTya5pRt/c5JNSZ5K8nCSZUPbrE2yI8n2JFcMjV+aZGuS\np5PcNje7JEk6mtkc6R8C/ryqfg34TeDmJL8K3Ap8o6ouBh4B1gIkuQS4DlgNXAXcniTdc30G+EhV\nrQJWJXnfSPdGknRMM0a/qvZW1ePd8n5gO7ASuAbY2K22Ebi2W74auKeqDlXVs8AOYE2Sc4GlVbWl\nW+/uoW0kSfPguM7pJ7kAeDvwH8DyqtoH0y8MwDndaiuAXUOb7enGVgC7h8Z3d2OSpHky6+gnORO4\nD/hYd8T/+u9v8PscJGnMLZzNSkkWMh38z1fV/d3wviTLq2pfd+rmx934HuC8oc1XdmNHGz/Sz/MF\nRJJOQFXlWI/P9kj/c8C2qvr00NgDwIe65Q8C9w+NX59kUZILgYuAzd0poBeSrOne2L1haJsjTXxi\nb+vWret9Du6b++f+TdZt3bp1s4r5jEf6SS4DPgA8meQxpk/jfBz4FHBvkg8DzzH9iR2qaluSe4Ft\nwCvATVV1+Mj9ZuAu4HTgwap6aFazlCSNxIzRr6pvA6cd5eHfPco2nwQ+eYTx7wK/fjwTlCSNjlfk\n9mAwGPQ9hTkzyfsG7t+pbpL3b7b7ltfOvIyPJDWO85KkcZaEGtEbuZKkCWD0JakhRl+SGmL0Jakh\nRl+SGmL0JakhRl+SGmL0JakhRl+SGjKrr1buw6ZNm/qeQu/OOuss3v3ud/c9DUkTZGy/hmHZsvf2\nPY3eHTjwbbZv38rb3va2vqci6RQwm69hGNsj/Rde8Eh/6dLVvPzyy31PQ9IE8Zy+JDXE6EtSQ4y+\nJDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE\n6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtS\nQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDVkxugnuSPJviRbh8bWJdmd5Hvd7cqhx9Ym2ZFke5IrhsYv\nTbI1ydNJbhv9rkiSZjKbI/07gfcdYfzvq+rS7vYQQJLVwHXAauAq4PYk6db/DPCRqloFrEpypOeU\nJM2hGaNfVd8CfnqEh3KEsWuAe6rqUFU9C+wA1iQ5F1haVVu69e4Grj2xKUuSTtTJnNP/aJLHk3w2\nybJubAWwa2idPd3YCmD30PjubkySNI8WnuB2twN/XVWV5G+AvwNuHN20ANYPLQ+6myTpsKmpKaam\npo5rmxOKflX9ZOjuPwP/2i3vAc4bemxlN3a08WNYfyJTk6RmDAYDBoPBq/c3bNgw4zazPb0Ths7h\nd+foD/sD4Pvd8gPA9UkWJbkQuAjYXFV7gReSrOne2L0BuH+WP1uSNCIzHukn+SLT51bOTrITWAdc\nnuTtwC+AZ4E/AaiqbUnuBbYBrwA3VVV1T3UzcBdwOvDg4U/8SJLmT15r8vhIUjB+85pvS5eu5tFH\nv8Lq1av7noqkU0ASqupIn6x8lVfkSlJDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcTo\nS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JD\njL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4k\nNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcTo\nS1JDZox+kjuS7EuydWjszUk2JXkqycNJlg09tjbJjiTbk1wxNH5pkq1Jnk5y2+h3RZI0k9kc6d8J\nvO91Y7cC36iqi4FHgLUASS4BrgNWA1cBtydJt81ngI9U1SpgVZLXP6ckaY7NGP2q+hbw09cNXwNs\n7JY3Atd2y1cD91TVoap6FtgBrElyLrC0qrZ06909tI0kaZ6c6Dn9c6pqH0BV7QXO6cZXALuG1tvT\nja0Adg+N7+7GJEnzaOGInqdG9DxD1g8tD7qbJOmwqakppqamjmubE43+viTLq2pfd+rmx934HuC8\nofVWdmNHGz+G9Sc4NUlqw2AwYDAYvHp/w4YNM24z29M76W6HPQB8qFv+IHD/0Pj1SRYluRC4CNjc\nnQJ6Icma7o3dG4a2kSTNkxmP9JN8kelzK2cn2QmsA/4W+HKSDwPPMf2JHapqW5J7gW3AK8BNVXX4\n1M/NwF3A6cCDVfXQaHdFkjSTvNbk8ZGk5uRtglPM0qWrefTRr7B69eq+pyLpFJCEqsqx1vGKXElq\niNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGX\npIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYY\nfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlq\niNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqyElFP8mzSZ5I8liSzd3Ym5NsSvJU\nkoeTLBtaf22SHUm2J7niZCcvSTo+J3uk/wtgUFXvqKo13ditwDeq6mLgEWAtQJJLgOuA1cBVwO1J\ncpI/X5J0HE42+jnCc1wDbOyWNwLXdstXA/dU1aGqehbYAaxBkjRvTjb6BXw9yZYkN3Zjy6tqH0BV\n7QXO6cZXALuGtt3TjUmS5snCk9z+sqp6PslbgE1JnmL6hWDY6+/P0vqh5UF3kyQdNjU1xdTU1HFt\nc1LRr6rnuz9/kuRrTJ+u2ZdkeVXtS3Iu8ONu9T3AeUObr+zGjmL9yUxNkibeYDBgMBi8en/Dhg0z\nbnPCp3eSLElyZrd8BnAF8CTwAPChbrUPAvd3yw8A1ydZlORC4CJg84n+fEnS8TuZI/3lwFeTVPc8\nX6iqTUm+A9yb5MPAc0x/Yoeq2pbkXmAb8ApwU1Wd4Kmfdhw8eJCXXnqp72n07rTTTmPRokV9T0M6\n5WUcuzv9QjJ+85pvZ555LQcPPtz3NMbCokWL2bnzR5x99tl9T0UaW0moqmN+FP5k38jVHNq//2t9\nT2FsLF78Vvbv32/0pZPk1zBIUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOM\nviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1\nxOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkMW9j0BaXbe\nygUXXECSvifSuze96Rx27foRZ5xxRt9T0SnI6OuU8POffxMoqvqeSf8OHHgLL730ktHXCTH6OkWk\nu8nfdnQyPKcvSQ0x+pLUEKMvSQ0x+pLUEKMvSQ0x+pLUEKMvSQ0x+pLUEKMvSQ1JzfN17UmuBG5j\n+gXnjqr61BHWKfB6e+lIFi9+C+9//5UsXnx631Pp3fnnr+ATn/grFizw+BWmr9auqmNesj2v0U+y\nAHgaeA/wX8AW4Pqq+uHr1pvw6E8Bg57nMFemmNx9g/HYv28y/c9oLjwFXDxHzz16CxbczIsv/g9L\nliyZ1fpTU1MMBoO5nVRPpqamuPzyy2eM/nx/984aYEdVPQeQ5B7gGuCHx9xq4kzRfzjmyhSTu28w\nHvv3O91tLqwH/niOnnv0Fiz42HGtP+nRn435jv4KYNfQ/d1MvxBI0gkIN974pyxc+IZZrf3EE9/h\nmWf2zvGc+rFz53/Oar2x/ZbNs876/b6nMGcOHnyK00//bt/TmBOTvG/g/o2bF188xJe+9Lnj2mbr\n1lNn/47P7L59db6jvwc4f+j+ym7s//nZz/5tXibUl5df3tH3FObMJO8buH8aV7N7H3S+38g9jel3\nit4DPA9sBv6wqrbP2yQkqWHzeqRfVf+b5KPAJl77yKbBl6R5Mu+f05ck9WesrmhIcmWSHyZ5Oslf\n9j2fUUpyR5J9Sbb2PZe5kGRlkkeS/CDJk0lu6XtOo5RkcZJHkzzW7d+6vuc0akkWJPlekgf6nsuo\nJXk2yRPd39/mvuczakmWJflyku3dv8F3HXXdcTnSn+2FW6eqJL8F7Afurqrf6Hs+o5bkXODcqno8\nyZnAd4FrJuXvDyDJkqo60L039W3glqqamIAk+TPgncBZVXV13/MZpSTPAO+sqp/2PZe5kOQu4JtV\ndWeShcCSqvrZkdYdpyP9Vy/cqqpXgMMXbk2EqvoWMJH/wwFU1d6qerxb3g9sZ/q6jIlRVQe6xcVM\nvx82HkdMI5BkJfB7wGf7nsscCePVu5FJchbw21V1J0BVHTpa8GG8/iMc6cKtiYpGK5JcALwdeLTf\nmYxWd/rjMWAv8PWq2tL3nEboH4C/YIJeyF6ngK8n2ZLk1LnkeHYuBP47yZ3d6bl/SvLGo608TtHX\nBOhO7dwHfKw74p8YVfWLqnoH09eXvCvJJX3PaRSSvB/Y1/2mFmZ7lc+p5bKqupTp32Zu7k63ToqF\nwKXAP3b7eAC49Wgrj1P0Z33hlsZTdy7xPuDzVXV/3/OZK92vzv8OXNn3XEbkMuDq7rz3l4DLk9zd\n85xGqqqe7/78CfBVJuvrX3YDu6rqO939+5h+ETiicYr+FuCiJG9Nsgi4Hpi0TxFM6lHUYZ8DtlXV\np/ueyKgl+aUky7rlNwLvZUK+KLCqPl5V51fVrzD97+6Rqrqh73mNSpIl3W+gJDkDuAL4fr+zGp2q\n2gfsSrKqG3oPsO1o64/Nd+9M+oVbSb7I9Ncznp1kJ7Du8BsvkyDJZcAHgCe7894FfLyqHup3ZiPz\ny8DG7lNmC4B/qaoHe56TZmc58NXpr2xnIfCFqtrU85xG7RbgC0neADwD/NHRVhybj2xKkubeOJ3e\nkSTNMaMvSQ0x+pLUEKMvSQ0x+pLUEKMvSQ0x+pLUEKMvSQ35P9g/pL8cLJA6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f34825152b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h = plt.hist(data.flatten(), bins=5, histtype='stepfilled')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quick IRAF Detour"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To write this data out for IRAF, uncomment and run the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# from astropy.io import fits\n",
    "#\n",
    "# hdu0 = fits.PrimaryHDU()\n",
    "# hdu1 = fits.ImageHDU(data.astype(np.float32))\n",
    "# hdu = fits.HDUList([hdu0, hdu1])\n",
    "# hdu.writeto('im_poisson_0.fits', clobber=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the FITS file, you also have to create an one-liner `im_poisson_0.coo` file with the following contents:\n",
    "\n",
    "    25 25\n",
    "    \n",
    "Inside your IRAF session, set the following parameters (more or less) using `epar`:\n",
    "\n",
    "    datapars.sigma = 1.\n",
    "    fitskypars.salgo = \"ofilter\"\n",
    "    fitskypars.annulus = 1.\n",
    "    fitskypars.dannulus = 20.\n",
    "    fitskypars.smaxiter = 10\n",
    "    fitskypars.snreject = 10\n",
    "    \n",
    "Then, run the following command:\n",
    "\n",
    "    fitsky im_poisson_0.fits im_poisson_0.coo\n",
    "    \n",
    "The result would be in `im_poisson_0.fits1.sky.1`:\n",
    "\n",
    "    MSKY = 0.7501966\n",
    "    STDEV = 0.7515768\n",
    "    SSKEW = 0.5060266"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Back to Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this dataset, Python version of the `ofilter` algorithm gives slightly lower sky and skew values, but comparable sigma. The Python version uses third-party libraries like Numpy, Scipy, and Astropy. Thus, it is not shocking that we are not getting complete agreement here.\n",
    "\n",
    "Some questions that could be pursued:\n",
    "\n",
    "1. Is it good enough? (Also see the next sub-section.)\n",
    "2. Do we even care about the skew? Maybe not? In Python, it is calculated using `scipy.stats.skew()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSKY = 0.696395704502\n",
      "STDEV = 0.747370581312\n",
      "SSKEW = 0.3443684652164396\n"
     ]
    }
   ],
   "source": [
    "# NOTE: Sigma clipping does not matter much for this dataset.\n",
    "ofil_results = ofiltsky.fitsky_ofilter(data, sigclip_sigma=None)\n",
    "\n",
    "print('MSKY =', ofil_results[0])\n",
    "print('STDEV =', ofil_results[1])\n",
    "print('SSKEW =', ofil_results[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare with some other available statistics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MEAN = 0.9868\n",
      "MEDIAN = 1.0\n",
      "MODE = 1\n"
     ]
    }
   ],
   "source": [
    "sky_mean = data.mean()\n",
    "sky_med = np.median(data)\n",
    "sky_mode = scipy.stats.mode(data, axis=None).mode[0]\n",
    "\n",
    "print('MEAN =', sky_mean)\n",
    "print('MEDIAN =', sky_med)\n",
    "print('MODE =', sky_mode)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This sub-section attempts to generate a plot not unlike what was published in WFPC2 ISR 1996-03 (Ferguson 1996). Perhaps the plot here can answer, \"Is it good enough?\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "i=0, val=0.0, errmsg=Unable to construct histogram\n",
      "i=1, val=0.1, errmsg=Unable to construct histogram\n",
      "i=2, val=0.2, errmsg=Unable to construct histogram\n",
      "i=3, val=0.3, errmsg=Unable to construct histogram\n",
      "i=4, val=0.4, errmsg=Unable to construct histogram\n",
      "i=5, val=0.5, errmsg=Unable to construct histogram\n",
      "i=6, val=0.6, errmsg=Unable to construct histogram\n",
      "\n",
      "Number of data points for plotting: 38\n"
     ]
    }
   ],
   "source": [
    "# Populate desired background values\n",
    "vals = np.arange(0, 4.5, 0.1)\n",
    "\n",
    "# Initialize arrays to store results\n",
    "sky_vals = []\n",
    "sky_ofil = []\n",
    "sky_med = []\n",
    "sky_mean = []\n",
    "\n",
    "# Generate results\n",
    "for i, val in enumerate(vals):\n",
    "    np.random.seed(i)\n",
    "    data = np.random.poisson(lam=val, size=(50, 50))\n",
    "    \n",
    "    try:\n",
    "        msky = ofiltsky.fitsky_ofilter(data, sigclip_sigma=None)[0]\n",
    "    except ValueError as e:\n",
    "        print('i={0}, val={1:.1f}, errmsg={2}'.format(i, val, str(e)))\n",
    "        continue\n",
    "        \n",
    "    sky_vals.append(val)\n",
    "    sky_ofil.append(msky)\n",
    "    sky_med.append(np.median(data))\n",
    "    sky_mean.append(data.mean())\n",
    "    \n",
    "# Convert result to Numpy arrays\n",
    "sky_vals = np.asarray(sky_vals)\n",
    "sky_ofil = np.asarray(sky_ofil)\n",
    "sky_med = np.asarray(sky_med)\n",
    "sky_mean = np.asarray(sky_mean)\n",
    "    \n",
    "print()\n",
    "print('Number of data points for plotting:', sky_ofil.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f3481f2ea90>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAEPCAYAAADxtOYjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8FPW9N/DPNyGBUBJIIESuCSpaWsUroo+ii60KVMSD\nRgFFVMT2Jd6PFRHbLGobRZ5WqR6tSBUsBBsfzxERKqgsLZ5DDRWURhCPmoAEIRgggERy+T5/zO5m\nd7O72U1mdvbyeb9e85rM7OzMbyazO9/9XUVVQURERGSWNLsTQERERMmFwQURERGZisEFERERmYrB\nBREREZmKwQURERGZisEFERERmcr24EJExojIdhHZISKzwmw3QkQaRWRiLNNHRERE0bE1uBCRNADP\nArgCwI8BTBaRH4bY7gkA78Q2hURERBQtu3MuzgPwuapWq2ojgOUAJgTZ7i4ArwPYF8vEERERUfTs\nDi4GANjls/y1e52XiPQHcLWqPg9AYpg2IiIi6gC7g4tIPA3Aty4GAwwiIqI41sXm4+8GMNhneaB7\nna9zASwXEQHQB8BYEWlU1RWBOxMRDpRCRBQlVeWPNjKV3TkXFQBOFpFCEckEMAmAX9Cgqie6pyEw\n6l3cESyw8NmekypKSkpsT0M8TLwOvBa8FuEnIivYmnOhqs0icieANTACnUWquk1Efm68rC8GviXm\niSQiIqKo2F0sAlX9K4BTA9b9McS2t8YkUURERNRhdheLkEUcDofdSYgLvA6teC1a8VoQWUuSqcxN\nRDSZzoeIyGoiAmWFTjKZ7cUiREQU/7Kysr5paGgosDsdFD+6deu299ixYycEe405F0REKSzSnAt+\nv1KgcPcO61wQERGRqRhcEBERkakYXBAREZGpGFwQEREFMXfuXEydOhUAsGvXLuTk5LBX0wixtQgR\nEVEIxrBWwKBBg1BfX29zahIHcy6IiIjIVAwuiIjIMjt3hl82w5AhQzB//nycccYZyM7OxowZM7Bv\n3z6MGzcOOTk5uPzyy3Ho0CEAwMaNG3HhhRciNzcXZ511FtavX+/dT1VVFRwOB3r27IkrrrgC+/fv\n975WXV2NtLQ0tLS0AABeeeUV/OhHP0JOTg5OPvlkvPhi61BY69evx6BBg/C73/0OBQUFGDBgAF55\n5RXzTzyOMbiwUXNz+GUiokTW0gJMnAg8/rix/PTTwNixQFOT+cd644038N5772HHjh1YsWIFxo0b\nhyeeeAL79+9Hc3MzFixYgJqaGlx55ZX49a9/jQMHDmD+/Pm45ppr8O233wIApkyZghEjRmD//v14\n5JFHsHjxYr9jeIpIAKCgoACrVq1CfX09Xn75Zdx3333YsmWL9/VvvvkGhw8fRk1NDV566SXMnDnT\nG+CkAgYXNmluBi66CHC5jGWXy1hmgEFEySItDXjrLWDpUiA/H/jDH4DVq4EuFtT2u+uuu9CnTx/0\n69cPo0aNwsiRIzF8+HBkZmbi3/7t3/DRRx/hz3/+M372s5/hiiuuAAD85Cc/wbnnnotVq1Zh165d\n2LRpEx599FFkZGRg1KhRGD9+fMjjjR07FkVFRQCAUaNG4fLLL8ff//537+uZmZn41a9+hfT0dIwd\nOxY9evTAZ599Zv6JxykGFzZJTwdKS4HiYsDpNOalpcZ6IuZqUbLo1w+YMAHYvx/46U+BwYOtOU5B\nQWvP5FlZWW2Wjxw5gurqavzlL39BXl4e8vLykJubiw8++AB79uxBTU0NcnNzkZWV5X1fYWFhyOOt\nXr0aF1xwAXr37o3c3FysXr3arxild+/eSEtrfcR2794dR44cMet04x6DCxs5HMDMmcDcucacAzUS\nwFwtSi5PPw2UlwP/+Afwt7+1FpHEmohg8ODBuOmmm1BXV4e6ujocOHAAhw8fxoMPPoh+/frhwIED\nOHbsmPc9O0NUEDl+/DiuvfZaPPjgg6itrcWBAwcwduxYNlP1weDCRi4X8NxzQEmJMfc8TCi1MVeL\nkkVLC/Dll8C6dcB55wHvvw/s2WNNnYtI3HjjjVixYgXWrFmDlpYWNDQ0YP369aipqcHgwYNx7rnn\noqSkBI2NjdiwYQPeeustv/d7gofjx4/j+PHj6NOnD9LS0rB69WqsWbPGjlOKWwwubNLcDMyebUT0\nTqcxnz2bv07JwFwtSgZpacCCBa1FIf36GT+kzK5z4VvRMtiyx4ABA7BixQr89re/RX5+PgoLCzF/\n/nxvC5ClS5di48aN6N27Nx577DFMmzYt6H579OiBBQsWoLi4GHl5eVi+fDkmTJgQVRqTHUdFtVFz\ns/+v0cBlSl0ul5FjMXOm8WVcXp7aAQY/K9bhqKjUURwVNU4Ffjnyy5IA5moFYh0UosTDnAuiOMRf\n6v6Yk2Md5lxQRzHngijBMFfLH+ugECUWBhdEFPfYsooosbBYhIjimqfORWmpkWPhchl1UDZsYI6O\nGVgsQh0V7t5hcEFxi/UOyIP3gnUYXFBHsc4FJRy2ECBfrINClFgsGD6GOoO/0Ay+vVT6thBIxWtB\nRJRomHMRR/hr3R9bCPjjYGZE0WtoaMD48eORm5uL66+/HsuWLcOYMWO8r6elpeHLL7+0NA3vvPMO\nJk6caMm+rUr/s88+i4ceeqjjO1DVpJmM00ls69ap9umjWlJizNetszlBNuK1aNXUpHr++a3XYN06\nY7mpyc5UUTJwf28m7ffrq6++qiNHjtSWlpagr6elpekXX3yhqqo333yz/upXvzI9Deeee65++OGH\n3mUR0R49emh2drYOHDhQ77///pDp8+VwOHTRokV+63zTb6aGhgYdOHCg1tbWhtwm3L1je86FiIwR\nke0iskNEZgV5fYqIfOyeNojI6XakM1b4a93AXir9cTAzoo6prq7GKaecEnJsDzWxkmpzkC+oTZs2\nob6+HiNGjPCuExF88sknqK+vx3vvvYdly5Zh4cKFHTqmmen31bVrV4wbNw5Llizp2A5CRR2xmGAU\ny/wvgEIAGQC2APhhwDbnA+jp/nsMgI1h9tfhKC1e8Nd6q8Bf5fyVbtwXgDEnMgOszrloaVF9+23V\nefNU33zTWDbZtm3b1OFwaK9evfS0007TFStWqKpqSUmJZmZmakZGhmZnZ+uf/vQnfeWVV/Siiy7y\nvldE9IsvvtAXX3xRMzIytGvXrpqdna1XXXWVqqrW1NToNddco/n5+XriiSfqggULvO91Op167bXX\n6o033qg9e/Zsk6ugqvroo4/qjBkz/NZ5julRXFysd911lz711FN6zTXX+G17991367333qtz5szR\n9PR0zcrK0uzsbL3rrru8+3rhhRd06NChmpubqzNnzvS+t6WlRR977DEtLCzUgoICnTZtmh46dEhV\nVauqqlREdPHixTp48GDNz8/X3/zmN37HXrp0qV566aUhr3u4e8fu4OJ8AKt9lh8CMCvM9r0A7Arz\nesiLkAiY9U3hMPAkK1geXNx9t+oPfqCakWHMb73VjGR7NTY26sknn6xPPPGENjY26vvvv6/Z2dm6\nY8cOVTUCgKlTp3q3f+WVV3TUqFHeZd8HfWCxSEtLi55zzjn6+OOPa1NTk3711Vd60kkn6Zo1a7z7\nzszM9AYzDQ0NbdJXXFys8+fP91vne8zKyko94YQT9OWXX9Y9e/Zojx49vAFAU1OT9u3bVzdv3qyq\nwYtFRETHjx+v9fX1unPnTs3Pz9d33nlHVVUXLVqkQ4cO1aqqKj169KhOnDjRey08wcXtt9+u33//\nvX788cfatWtX3b59u3ffH330kfbu3TvktQ9379hdLDIAwC6f5a/d60K5DcBqS1Nko/R0o2MgT1GI\nw8GOgsjAYqK2WME1AezeDbz4InD0KNDYaMzLyoAdO0w7xMaNG3H06FHMmjULXbp0wejRo3HllVei\nrKys0/uuqKjA/v37MWfOHKSnp6OoqAi33XYbli9f7t3mggsuwPjx4wEYRQmBDh48iOzs7Dbrzz77\nbPTu3RsTJkzA7bffjptvvhknnHACLr74YpSXlwMAVq9ejfz8fJx55plh0zl79mxkZ2dj0KBBGD16\nNLZs2QIAWLZsGe6//34UFhaie/fuKC0txfLly71DzIsInE4nMjMzMXz4cJxxxhn4+OOPvfvNzs7G\noUOHorxqhoRpiioiowHcAuAiu9NiJbbnp2A8gafnfkj1wJO9diaIujogIwNoaGhdl5EBfPutaYeo\nqanBoEGD/NYVFhZi9+7dnd53dXU1du/ejby8PABGTn9LSwsuvvhi7zaBxw6Um5uLw4cPt1m/efNm\nDBkypM36m266CS+88AKmT5+OpUuXYurUqe2ms6CgwPt39+7dceTIEQDGtSksLPS+VlhYiKamJuzd\nu7fd9wLA4cOH0bNnz3aPH4zdwcVuAIN9lge61/kRkeEAXgQwRlUPhNuh0+n0/u1wOOBI1RqRlHQY\neLZiPygd53K54IrV4CxDhwJZWcCRI4Cn4mGXLsBpp5l2iP79+2PXrl1+63bu3IlTTz016n0FVvoc\nNGgQTjzxRHz22WcRvyfQ8OHDsSNITo2GqIh59dVX44477kBlZSVWrlyJp556KuJjBerfvz+qq6u9\ny9XV1cjIyEBBQUGbaxbMtm3bcMYZZ0R1TA+7g4sKACeLSCGAPQAmAZjsu4GIDAbw/wBMVdUv2tuh\nb3BBRMnLt2VVSUnqtqyKVuCPrrlz51p3sG7djGyla68FPv8cKCoyosAgxQQdNXLkSHTv3h3z5s3D\n/fffjw0bNmDlypUdehYUFBT49Rlx3nnnITs7G/PmzcPdd9+NjIwMbN++HceOHcO5554b0T7HjRuH\nSZMmRZyGrl274pprrsGUKVMwcuRIDBw4MGT62jN58mTMmzcPY8aMQZ8+fTBnzhxMmjQJaWlGjYhQ\nAY7H+vXrMXbs2IiP58vWOheq2gzgTgBrAFQCWK6q20Tk5yJyu3uzXwHIA/AfIrJZRD60KblEFEc4\nUmqCGDYMqKwEjh836lp08JdwKBkZGXjrrbewatUq9OnTB3feeSdeffVVDB06NKL3++YGTJ8+HZWV\nlcjLy8PEiRORlpaGlStXYsuWLRgyZAj69u2LGTNmoL6+PuL0nXXWWejVqxcqKiqCHjOYadOmYevW\nrbjpppv81t9zzz0oLy9H7969ce+99wbdl+/yrbfeiqlTp+Liiy/GSSedhO7du2PBggUh0+G73NDQ\ngFWrVmHatGkRnqk/DlxGRAmHdS7Mw4HLrLd27Vo8//zzeOONNyLafteuXRg2bBi++eYb9OjRw+LU\nBffss8/i66+/xhNPPBFyG46KSkRJh+PwmIPBRXxpaWnB/fffjyNHjuCll16yOzlhhbt37K5zQUTU\nIazgSsnmu+++Q0FBAYYMGYLVqxO71wXmXBARpTDmXFBHhbt37O5Ei4iIiJIMgwsiIiIyFYMLIiIi\nMhWDCyIiIjIVgwsiIiIyFYOLBMJRIImIotPQ0IDx48cjNzcX119/PZYtW4YxY8Z4X09LS4uqS+2O\neOeddzBx4kRLjxHM1q1bceGFF8b8uACDi4Th6ZHQ08Wxy2UsM8AgIgrt9ddfR21tLerq6vDaa69h\nypQp+Otf/+p93bfL61tuuQW//vWvTU/DI488gtmzZ3uXfQOauXPnIjMzEzk5OcjLy8NFF12EjRs3\nttlHVVUV0tPTMXPmzDavpaWlITs7Gzk5OcjOzvaO4nr66acjNzcXb7/9tunn1B4GFwnCdxRIp9OY\nl5ay4yAionCqq6txyimnhBzPw8y+O5qD/NrbtGkT6uvrMWLECO+6wLRMmjQJ9fX12L9/PxwOB4qL\ni9vsZ8mSJcjLy8Nrr72GxsZGv9dEBJ988gnq6+tx+PBh1NXVeV+bMmUKXnjhhc6eWtQYXCQQ31Eg\nZ87kKJAAi4qI4p2qYtWqVXjqqaewYsUKUx/mHtu3b8fo0aORm5uL008/HW+99RYAY5TsRx99FMuX\nL0dOTg5efvllLF68GKNGjWqzj4ULF2Lp0qWYN28ecnJyMGHCBADAnj17cO2116Jv37446aST8Ic/\n/MH7nrlz56K4uBhTp05Fr169sHjx4jb7Xb16NS655BK/daGuQVpaGm644QbU1NTg22+/9XttyZIl\nePzxx70DtQXuL9Q+HQ4H3nvvvTYBidUYXCQQjgLpj0VFRPHv3nvvxXXXXYc5c+ZgypQpuO2220zd\nf1NTE8aPH48xY8agtrYWCxYswA033IDPP/8cTqcTDz/8sDdn4JZbbgEQfFTSGTNm4IYbbsCDDz6I\n+vp6vPnmm1BVjB8/HmeddRb27NmD9957D8888wzWrl3rfd+KFStw3XXX4eDBg7jhhhva7Hfr1q04\n9dRTIzqX48ePY/Hixejduzdyc3O96//+979j9+7dmDRpEoqLi4MGMaH0798fGRkZ+OyzzyJ+jxkY\nXCSI5mZj1MfycqNYpLzcWE7lBymLioji2+7du/Hiiy/i6NGjaGxsxNGjR1FWVoYdO3aYdoyNGzfi\n6NGjmDVrFrp06YLRo0fjyiuvRFlZWaf3XVFRgf3792POnDlIT09HUVERbrvtNixfvty7zQUXXIDx\n48cDALp27dpmHwcPHkR2dnbY47z22mvIy8tD9+7dsWjRIrz++utIS2t9PC9ZsgTjxo1Dz549vXVG\n9u/f77ePs88+G7m5ucjLy/MOx+6RnZ2NgwcPRn3+ncHgIkGkpxvDSXuKQhwODi8NsKiIKJ7V1dUh\nIyPDb11GRkabLP/OqKmpwaBBg/zWFRYWYvfu3Z3ed3V1NXbv3o28vDzk5eUhNzcXpaWl2Ldvn3eb\nwGMHys3NxeHDh8Nuc/3116Ourg779u3Daaedhk2bNnlfa2hoQHl5OaZMmQIAOP/88zFo0CAsW7bM\nbx+bN2/GgQMHUFdXh6efftrvtcOHD6NXr14RnbNZGFwkEI4C2RaLioji19ChQ5GVleVXDNGlSxec\ndtppph2jf//+2LVrl9+6nTt3YsCAAVHvK7C4ZNCgQTjxxBNRV1eHuro6HDhwAIcOHfKr8xCqoqjH\n8OHDI86pycvLwx//+Ec4nU7s3bsXAPDGG2+gvr4ed9xxB/r164d+/fqhpqamTdFIqDoXNTU1aGxs\njLhoxiwMLihhsaioLVZwpXjSrVs3uFwuDBs2DBkZGRg6dCjef//9dosJojFy5Eh0794d8+bNQ1NT\nE1wuF1auXInJkydHva+CggK/Pi/OO+88ZGdnY968eWhoaEBzczMqKyv9chbaM27cOLii+NVzyimn\nYMyYMXjyyScBAIsXL8b06dOxdetWfPzxx/j444+xYcMGbNmyBZWVle3ub/369bj00kvb5CBZjcEF\nJSwWFfljBVeKR8OGDUNlZSWOHz+OHTt24IwzzjB1/57WE6tWrUKfPn1w55134tVXX8XQoUMjer9v\nzsP06dNRWVmJvLw8TJw4EWlpaVi5ciW2bNmCIUOGoG/fvpgxYwbq6+sjTt9ZZ52FXr16oaKiIugx\ng3nggQewcOFC7Ny5E+vWrcN9992Hvn37eqezzz4bY8eO9eZehNvf0qVL8Ytf/CLi9JpFrGgWZBcR\n0WQ6H6JouVxGxdaZM41iovJy1kOh8EQEqhr+aQd+v3bG2rVr8fzzz+ONN96I6XG3bt2KX/ziF/jg\ngw8s2X+4e4fBBVGScTqNCq4lJcbfROEwuKCOCnfvsFiEKImwgisRxQPmXBAlCU+di9JSoyjE5TIq\nuKZyPRRqH3MuqKNYLJLkmpv9Hx6By5Q6eC+0xWsSHoML6igWiyQxthAgX+wLxR8/H0T26GJ3Aqhz\nfLvA9m0hkOoPFSKAnw8iuzC4SAK+XWCXlLDpIZEvfj7M0a1bt70iUmB3Oih+dOvWbW+o11jnIgmw\nbwOi0Pj5CC/SOhdE0WBwkeDYQoAoNH4+2sfggqzA4CIJsDY8UWj8fITH4IKsYHtrEREZIyLbRWSH\niMwKsc0CEflcRLaIyJmxTmO8YwsBotD4+SCKPVuDCxFJA/AsgCsA/BjAZBH5YcA2YwGcpKpDAfwc\nwAsxTygRERFFzO6ci/MAfK6q1araCGA5gAkB20wAsAQAVPUfAHqyxjIREVH8sju4GABgl8/y1+51\n4bbZHWQbIiIiihN2BxemE5E2kzPE0JBOp5Pbc3tuz+1Tavubb74ZTqfTOxFZwdbWIiJyPgCnqo5x\nLz8EQFX1SZ9tXgCwTlVfcy9vB3CJqrbpvENStLUIhcaWAq14LSgYEbYWIfPZnXNRAeBkESkUkUwA\nkwCsCNhmBYCbAG8wcjBYYEEUiONKtOK18Bd43ql6HYisYmtwoarNAO4EsAZAJYDlqrpNRH4uIre7\nt1kF4CsR+V8AfwRwh20JpoTiO66E02nMS0tT89c6r0UrBlpE1mMnWpT0nM7WcSVSvYiZ18LALsFb\nsViErGB3sQiRpVwu4+FRUmLMPb9WUxGvRSvfwcxmzkzdwILIKsy5oKTFcSVa8Vr4Y85FK+ZckBUY\nXFBSYwuJVrwWBgZa/hhckBUYXBBRymGg1YrBBVmBdS6IKOVwMDMiazG4ICIiIlMxuCAiIiJTMbgg\nIiIiUzG4IKKUx+7AiczF4IKIUhq7AycyXxe7E0CxwaZ3rXgtyJfvuCu+nWrxniDquIhyLkQkXUT6\ni8hgz2R1wsg8/GXWiteCgmF34ETmarcTLRG5C0AJgL0AWtyrVVWHW5y2qLETrdDY3XErXgt/zMlJ\n7XuCnWiRFSLJubgHwKmq+mNVPd09xV1gQeHxl1krXotWzMkxznX2bCOgcDqN+ezZqXUNiMwWSc7F\nOgCXqWpTbJLUccy5CC2Vf5kF4rXwx+uR2rk3zLkgK0QSXCwCcCqAtwF871mvqr+zNmnRY3ARHAdq\nasVrEZzTaeTklJQYf1PqYHBBVogkuCgJtl5V51qSok5gcBFaKv8yC8Rr4Y85F6mNwQVZgaOiEqUw\n5uQQgwuyQiQ5F/kAHgTwYwDdPOtV9VJrkxY9BhdE0WNOTlupdE0YXJAVImktshTAdgBDAMwFUAWg\nwsI0EVEMcfhxf2xBQ9R5keRc/FNVzxGRTzxNUEWkQlVHxCSFUWDOBRGZIZXqoTDngqwQSc5Fo3u+\nR0R+JiJnAcizME0UAxyoiSg09oVC1DmRBBePi0hPAP8O4AEALwG4z9JUkaWY7euPgRYFcrmMHIuS\nEmPu+awQUWTYWiRFpVK2bzhsLUGBUu2eYLEIWSGSOhenAHgeQIGqniYiwwFcpaqPxyKB0WBwER12\nnGRgoEWB2FqEqHMiKRZZCGA23HUvVPUTAJOsTBRZj9m+rVi+ToHYgoaocyIJLrqr6ocB6+J+nBEK\njQM1+WOgRURkri4RbLNfRE4CoAAgItcC2GNpqshS6en+5ccOR/KWJ7fHN9ByOIwpmcvXiYhiIZI6\nFycCeBHA/wFwAMBXAG5U1apOHVgkF8BrAAphdMx1naoeCthmIIAlAAoAtABYqKoLwuyTdS4oaqlU\nvk4UiHUuyAoRtxYRkR8ASFPVw6YcWORJAN+q6jwRmQUgV1UfCtjmBAAnqOoWEekB4J8AJqjq9hD7\nZHBBRBQFBhdkhUhyLnoBuAlAEXyKUVT17k4dWGQ7gEtUda87iHCp6g/bec9/AfiDqr4X4nUGF0Rk\njtpaoKoKKCoC8vPtTo1lGFyQFSKp0LkKRmCxFUbOgWfqrL6quhcAVPUbAH3DbSwiRQDOBPAPE45N\nRBRaWRlQWAhcdpkxLyuzO0VECSWSCp3dVPX+juxcRNbCqC/hXQWjYugjQTYPmeXgLhJ5HcA9qnok\n3DGdPh02OBwOONiukIiiUVsLTJ8OHDtmTICx/NOfJkUOhsvlgotNoshikRSL3AfgCICVAL73rFfV\nuk4dWGQbAIdPscg6VR0WZLsu7mOvVtVn2tkni0WIqHMqKowci0M+9ctzcoB33wVGxN14jZ3GYhGy\nQiTFIscBPAXgf9BaJLLJhGOvAHCz++9pAN4Msd2fAHzaXmBBRGSKoiLg+HH/dY2NxnoiikgkORdf\nAjhPVfebemCRPAB/ATAIQDWMpqgHRaQfjCanV4rIhQD+BqO+h7qnh1X1ryH2yZwLIuq8sjKjKCQj\nwwgsFi0CJk+2O1WWYM4FWSGS4GINgKtV9bvYJKnjGFwAtbW1qKqqQlFREfKToHyYyDZsLULUYZEE\nF/8J4McA1sG/zkWnmqJaIdWDi7KyMkyfPh2ZmZk4fvw4Fi1ahMlJ+muLiMzB4IKsEElwMS3YelVd\nbEmKOiGVg4va2loUFhbimKd2O4CsrCxUV1czB4PILgmQ+8HggqzQblPUeAwiqK2qqipkZmb6BRcZ\nGRmoqqpicEGdxi7SO8BTbyMz06ggmsT1NogCRdJahBJAUVERjgfUcG9sbERRKtRwr601mg/W1tqd\nkqTU3AxcdFHraLEul7GcqqPoRsS3r4xDh4z59Om8RyllMLhIEvn5+Vi0aBGysrKQk5ODrKwsLFq0\nKPlzLczqSZEBSkjp6UBpKVBcDDidxry0NE5yLuL1/1ZVZeRY+MrIMNYTpYCIBy4DjIHE3F11x6VU\nrnPhEVVrkQQoDw6rttYIKHyKgpCVBVRXt38+vuf+7rvMvo6A0wnMnQuUlBh/x1zg/drZYgcr7//O\n3JsxxjoXZAlVjXgC8FE028d6Mk6HIrJsmWpWlmrPnsZ82TK7UxS9Dz800g+0Tjk5xvpwfM+9WzfV\nzEz/fWRlqe7b5928qcn/7YHLSWHfPuO6+Zy3r3XrVPv0US0pMebr1sUycdr2fn3hBWMe5v8W1f6s\nuP89x8jJievPmPt70/bvb07JNUW3MbDZ7gS3kz6lCOzb17kv5njRkfMI9p7AySdAaWpSPf/81ofp\nunXGclIFGO08aG2/BsH+Z127qmZnh/y/Rb0/q+7/doK2eMDggpMVU7R1LhaamGlCdkmW8uD8fCMr\nPCvLGPshK8tYDpftHOzcA/l09RzX9Q3MEEHFw/R0YMMGwDMGoMNhLMfsGoS6XzvaRXcs7//8fGM8\nkjgrCiGyWlTBhar+h1UJoRhKprETJk82yrHffdeYt1fmHuzcMzLCBigOB/DLm2vx9twK/PLmWiTV\nQLsRPmhS/U1xAAARuUlEQVTDNUO1vNVIsP9ZczPwzDPRBZbh9peo9z9RvLI768TMCSwWiVyClAdb\nkq0c7NzDHOdfc5bpUWTpsa499Siy9F+PBFyreMv6jiY9URYR2FZEEup+jfRcA7dLlPs/BsBiEU4W\nTLYnwNSTSYXgwswHWbw9FAOZVeku2HlGeO5Ne/bpsTT/h++xtCxt2hPwkIqXirEdSU+UD1rbKnd2\n9H4NdU3i/f6PEQYXnKyYQr8ADA7z2ii7Ex4iXZrU4u1BZiWzKt119pp9+KG25Pi3SGnxVByMt4qx\nnUlPlA/akhJj9yUlnUqx9eLtfxSHGFxwsmIKV+fCJSIPioi3dFVECkTkzwB+b0kZDYXWXsW7SDoT\n6kCHQ4Hl6THrldGMSndm9JJYVARp9C+fF0/5vB0VY8P9DzuTnigqHrpcwHPPGf1dPPdca8+dcSlZ\nKi8TJZhwwcU5AE4CsEVELhWRewB8COB/AJwXi8SRj3BfkpH0UtmBnixt7fbZjEp3ZjxYwrVI6Uga\nO9OjZHv/Q6srKtbWonljBZ58oBbl5UbrmfJyYPbsKO6JWPeoycqbRPZoL2sDwD0AWgB8DWCg3Vkt\n7aRVE040FdKCZe9++mn72b6dyBq2tfOkzla6MzNLPNT/KZo0dqaIJtJzaS89JtRbaPF0YuXeT8SV\nOe0q1mPlzbDAYhFOFkyhXwB6AfgjgC0ALgfwNICtAC61O9Fh0qwJJdov22BfkpH0UtnRnizdbC1f\n72ylu1g8WCJJY2cDnWj+h+0FQtE+3EN1PJadHfl+Ij1/qypZsvJmSAwuOFkxhX4B+BLAAwC6+Kw7\nE8B/AyizO+Eh0qwJI9yXbbgvwsDXIvnSTtScC7PEw4OlkwFep4OTzrw/WNqD5aCFu8aRnH+yVViO\nh/suAgwuOFkxhX4hTBEIgBl2JzxEujRhhPqyfewxa5oSduAXvO3dPieTjnZVblbfDJ0JbtrrMr1b\nN6M77nD3bHvnHw+tOswMBhIoUGJwwcmKyfYEmHoyiRRcBPsy7dbN2qaEHfjyTIlBu2LFjPoZHX0A\ndvbh7UlPjx4aMshob7/hzr+zOTudZWYwEA+BUhQYXHCyYopqyPV4l3BDrnuGjM7IMGqwP/wwMH++\n0WzSIyfH6Np6xAj70pnoQ7PHk0iupVXDdQfebx0dovyjj4D77jP28/33QFqaf1rD3bOhzt/OIcrN\nPnZFhdGiJ94+xyFwyHWyAoMLu/l+2QL2fcGG4nkgZWYaTfqifSBR9Kx8OJkVKHr206MHcM455tyz\nnQ1+Osrs621noNQBDC7ICgwu4o1dX7DBJNiXZNJIsOvesrQMaTNa79mWhYuQdkMH71k7csmsuN7x\n9DluB4MLsgKDi3gUL8UQCZa9m1QS5OHk6Wht/qxaXDigCh/sLsIDT+bHdkh2M1hxvePlc9wOBhdk\nBQYXFFqC/YJOOgnycHK5gOJiYOZMozvw8nIk5rD0CXK9zcbggqzA4ILCS5Bf0GQvpxOYO9cYb8Tp\ntDs19mtu9s+5CVyOJwwuyArhxhYhMgKJ6mqjKKS6moFFANsGdosjCTWQWQzYOiYPUZxgcEHti2LE\nzFTS6YdIrAfxskBzszFwWYcHMktC6elAaalRVOR0GvPS0vjNuSCygm3FIiKSC+A1AIUAqgBcp6qH\nQmybBmATgK9V9aow+2SxCMVUh+sbJFET30QqAoilRCkqYrEIWcHOnIuHALyrqqcCeB/A7DDb3gPg\n05ikiigKDocRWMyda8wjCixqa43A4tgxoyXOsWPGcoLmYAQGEgwsjKDz2Wf9i4pSOTeHUo+dwcUE\nAIvdfy8GcHWwjURkIIBxAF6KUbooAqxrYOhQfYOqKiPHwldGhrGeEl5zM/DQQ0BBgRFslpcDd9wB\nXHhh6n5OKPXYGVz0VdW9AKCq3wDoG2K73wP4JQCWd8QJVlgzdLi+QVGRURTiq7GxtZdWSmjp6cAH\nHxjBZnGx8fmorQWeeIK5OpQ6uli5cxFZC6DAdxWMIOGRIJu3CR5E5GcA9qrqFhFxuN8fltOncNPh\ncMCRkA3u45tvhTXfugap9sWZng6/zqIcDkTWeVR+vlHHIrCJLyvMJo30dP8is5KS+On7w+VywZXq\nTXrIcnZW6NwGwKGqe0XkBADrVHVYwDa/BXAjgCYAWQCyAbyhqjeF2CcrdMbQvF/Wonx+FYofKMKD\nT/HBGLUU7bQpVSRK52Ks0ElWsDO4eBJAnao+KSKzAOSq6kNhtr8EwL+ztUh8qHykDEN+Mx1pXTPR\n8v1xfPXIIvz4scRs7UBkNk/RYWmpEVC4XEaRWTx2i87ggqxgZ3CRB+AvAAYBqIbRFPWgiPQDsFBV\nrwzYnsFFnGj+phaNAwrRraW1W/CGtCxk7K5G+gn8BU4EJE4TXQYXZAV2/03Rq6iA/vQySH1rtySa\nkwPhgGZECYfBBVmBPXRS9IqKII3+rR2ErR2IiMiNwQVFz9PaISvLGII9K4utHYiIyIvFItRx7bV2\nYGsIorjHYhGyAoMLskYSjZ1BlMwYXJAVGFyQ+WprgcJCY8wMj6wsY8h25mAQxRUGF2QF1rkg83Hs\nDCKilMbggszHsTOIiFIagwsyH1uTEBGlNNa5IOuwtQhR3GOdC7ICgwsiohTG4IKswGIRIoqJ5ubw\ny0SUPBhcEJHlPKOEulzGsstlLDPAIEpOXexOABElv/R0Y/jx4mJg5kzgueeA8vL4HCWUiDqPORdk\niXazwGtrgYoKY04pweEwAou5c425w2F3iojIKgwuyHTtZoGXlRk9eF52mTEvK7MppRRLLpeRY1FS\nYsw99wcRJR+2FiFLuFxts8AdDrBr8BTlCThLS437wOUCZs8GNmxg0Yjd2FqErMA6F2QJ3yzwkhKf\nLHBP1+C+wYWna3AGF0krPd0/kHA4GFgQJTMWi5AlQmaBs2vwlBUYSDCwIEpeDC7IdM3NRpZ3eTng\ndBrz2bPddS7YNTgRUdJjnQuyRHOz/y/TwOVU7Bq83WtCZAPWuSArMOeCLNFuFnh+PjBiREoFFuxE\niohSBSt0EsUAO5EiolTCnAuiGGEnUkSUKhhcUPxI8l472YkUEaUKBhcUH5K8186wLWiIiJIMW4uQ\n/VKk1062FqF4xNYiZAXmXJD9PL12+vL02plE2IkUEaUKBhdkP/baSUSUVGwLLkQkV0TWiMhnIvKO\niPQMsV1PESkXkW0iUikiI2OdVrJYCvba2e6Q9ERECczOnIuHALyrqqcCeB/A7BDbPQNglaoOA3AG\ngG0xSh/F0uTJRh2Ld9815pMn250iy7BDLSJKdrZV6BSR7QAuUdW9InICAJeq/jBgmxwAm1X1pAj3\nyQqdlBBCDklPFGOs0ElWsDPnoq+q7gUAVf0GQN8g2wwBsF9EXhaRj0TkRRHJimkqiSzADrWIKJlZ\nGlyIyFoR+cRn2uqeXxVk82BZDl0AnA3gOVU9G8B3MIpTiBIaO9QiomRm6dgiqnpZqNdEZK+IFPgU\ni+wLstnXAHap6ib38usAZoU7ptPp9P7tcDjg4E/CuJPq/T34dqjlcBjT7NnAhg2pdR3IHi6XCy5G\ns2QxO+tcPAmgTlWfFJFZAHJVtU2uhIisBzBDVXeISAmA7qoaNMBgnYv456nMWFpqPFRdrtR8sKZ6\ngBUMr4k9WOeCrGBnnYsnAVwmIp8B+AmAJwBARPqJyEqf7e4GsFREtsBoLfLbmKeUTOM7OqjTacxL\nS1PvIcIOtfyxBQ1RcmH332QLp9OozFhSYvxNxBY09mDOBVmBPXRSzLEyIwXDFjREyYPBBcUURwel\nUBh0EiUPFotQzLHiHgViRV/7sFiErMDggojiAoNOezC4ICuwWISI4gJb0BAlDwYXREREZCoGF0RE\nRGQqBhdERERkKgYXREREZCoGFxQXAvu5YL8XRESJi8EF2Y7jShARJRdLh1wnioTvYGa+40qwKSIR\nUWJizgXFBY4rQUSUPBhcUFzguBJERMmD3X+T7TiuBJF92P03WYHBBcUFjitBZA8GF2QFFotQXOC4\nEkREyYPBBREREZmKwQURERGZisEFERERmYrBBREREZmKwQURERGZisEFERERmYrBBREREZmKwQUR\nERGZisEFERERmYrBBREREZmKwQURERGZyrbgQkRyRWSNiHwmIu+ISM8Q290nIv8SkU9EZKmIZMY6\nrURERBQ5O3MuHgLwrqqeCuB9ALMDNxCR/gDuAnC2qg4H0AXApJimMkG5XC67kxAXeB1a8Vq04rUg\nspadwcUEAIvdfy8GcHWI7dIB/EBEugDoDqAmBmlLePzyNPA6tOK1aMVrQWQtO4OLvqq6FwBU9RsA\nfQM3UNUaAP8XwE4AuwEcVNV3Y5pKIiIiikoXK3cuImsBFPiuAqAAHgmyuQZ5fy8YORyFAA4BeF1E\npqjqMguSS0RERCYQ1TbP9NgcWGQbAIeq7hWREwCsU9VhAdtcC+AKVZ3hXp4KYKSq3hlin/acDBFR\nAlNVsTsNlFwszbloxwoANwN4EsA0AG8G2WYngPNFpBuA7wH8BEBFqB3yA0JERGQ/O3Mu8gD8BcAg\nANUArlPVgyLSD8BCVb3SvV0JjBYijQA2A7hNVRttSTQRERG1y7bggoiIiJJTUvTQKSJjRGS7iOwQ\nkVl2p8cuIrJIRPaKyCd2p8VuIjJQRN4XkUoR2Soid9udJruISFcR+YeIbHZfixK702Q3EUkTkY9E\nZIXdabGTiFSJyMfue+NDu9NDySPhcy5EJA3ADhj1MWpg1MmYpKrbbU2YDUTkIgBHACxxdzqWstyV\nhE9Q1S0i0gPAPwFMSMX7AgBEpLuqfici6QA+AHC3qqbsw0RE7gNwDoAcVb3K7vTYRUS+BHCOqh6w\nOy2UXJIh5+I8AJ+rarW7LsZyGM1XU46qbgDALwkYfaeo6hb330cAbAMwwN5U2UdVv3P/2RVGRe7E\n/lXRCSIyEMA4AC/ZnZY4IEiO5wDFmWS4qQYA2OWz/DVS+CFCbYlIEYAzAfzD3pTYx10MsBnANwDW\nqmrIVlcp4PcAfokUDrB8KIC1IlIhIjPsTgwlj2QILohCcheJvA7gHncORkpS1RZVPQvAQAAjReRH\ndqfJDiLyMwB73bla4p5S2YWqejaMnJyZ7qJVok5LhuBiN4DBPssD3esoxbnHo3kdwKuqGqwflZSj\nqvUA1gEYY3dabHIhgKvcdQ3KAIwWkSU2p8k2qrrHPa8F8J8wipmJOi0ZgosKACeLSKF7OPZJMDro\nSlX8NdbqTwA+VdVn7E6InUSkj4j0dP+dBeAyAClZsVVVH1bVwap6IozvivdV9Sa702UHEenuztmD\niPwAwOUA/mVvqihZJHxwoarNAO4EsAZAJYDlqrrN3lTZQ0SWAfhvAKeIyE4RucXuNNlFRC4EcAOA\nS93N7D4SkVT9td4PwDoR2QKj3sk7qrrK5jSR/QoAbHDXxdkI4C1VXWNzmihJJHxTVCIiIoovCZ9z\nQURERPGFwQURERGZisEFERERmYrBBREREZmKwQURERGZisEFERERmYrBBREREZmKwQURERGZisEF\nUZTcXc1vE5GXReQzEfmziPxERDa4l891d628SEQ2isg/RWS8z3v/JiKb3NP57vWXiMg6ESl37/tV\ne8+SiKjj2EMnUZREpBDA5wDOVNVPRWQTgC2qeps7iLgVwKcAKlV1mXtcjw9hDPuuAFpU9biInAyg\nTFVHiMglAP4LwI9gDIv+AYAHVPW/Y3+GRESd08XuBBAlqK9U9VP335UA3nP//S8ARTBG5x0vIr90\nr8+EMXrvHgDPisiZAJoBDPXZ54eeUSrd44AUwRgrhogooTC4IOqY733+bvFZboHxuWoCcI2qfu77\nJhEpAfCNqg4XkXQAx0Lssxn8fBJRgmKdC6KOaW9Y+3cA3O3d2MipAICeMHIvAOAmAOnmJ42IyF4M\nLog6RkP87Vl+DECGiHwiIlsBPOp+7T8A3Owe5voUAEcj2D8RUUJhhU4iIiIyFXMuiIiIyFQMLoiI\niMhUDC6IiIjIVAwuiIiIyFQMLoiIiMhUDC6IiIjIVAwuiIiIyFQMLoiIiMhU/x/PcJbwihZaowAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f34b068b128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(sky_mean, sky_med - sky_mean, color='b', marker='x', label='median')\n",
    "plt.scatter(sky_mean, sky_ofil - sky_mean, color='r', marker='o', label='ofilter (Python)')\n",
    "plt.scatter([0.9868], [0.7501699 - 0.9868], color='k', marker='o', label='ofilter (IRAF)')\n",
    "plt.xlabel('mean')\n",
    "plt.ylabel('X - mean')\n",
    "plt.axhline(0, color='k', linestyle='--')\n",
    "plt.legend(bbox_to_anchor=(1.01, 1), loc='upper left', scatterpoints=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f3481eb5eb8>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAEPCAYAAADxtOYjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X981XX9///bYz+AEUw3GKCIGyIWlZo/0QydmoUUb8qk\nRPJXSPV5S1pW/kh7M8pSyfe7Mvtm4t4KpsOw3kUGiT84VL6jMEFoglq5iYC+R8hPmWxnj+8fr7Nx\nNrYxttc553W2+/VyOZdzXq/z2uv13JPD2eP1eP4yd0dEREQkLDmZLoCIiIj0LgouREREJFQKLkRE\nRCRUCi5EREQkVAouREREJFQKLkRERCRUGQ8uzGyimW0ws5fN7MZ23i80s8VmtsbM1pnZlRkopoiI\niHSRZXKeCzPLAV4Gzgc2A6uAS9x9Q9IxNwOF7n6zmQ0FXgKGu3tjJsosIiIinct05uJ04BV3r3X3\nBmAhMKXNMQ4MTrweDPxLgYWIiEh0ZTq4GAlsTNp+PbEv2T3Ae81sM/ACcF2ayiYiIiLdkOngois+\nCqx29yOBk4Afm9mgDJdJREREOpCX4etvAo5O2j4qsS/ZVcDtAO7+DzN7FXgP8Fzbk5mZFkoRETlE\n7m6ZLoP0LpnOXKwCjjWzUjPrB1wCLG5zTC3wYQAzGw4cB/yzoxO6e6Qfs2fPzngZVE6VU+VUOZsf\nIqmQ0cyFu8fNbBawjCDQqXT39Wb2heBtvw+4DXjQzNYmfuwGd9+WoSKLiIjIQWS6WQR3/x3w7jb7\nfpr0egtBvwsRERHJApluFulzysvLM12ELlE5w6VyhkvlFIm2jE6iFTYz8970+4iIpJqZ4erQKSHL\neLOIiIhEX0FBwRv19fXDM10OiY4BAwa8uXfv3hHtvafMhYhIH9bVzIW+X6Wtzj476nMhIiIioVJw\nISIiIqFScCEiIiKhUnAhIiIioVJwkSbxeOfbIiIivYWCizSIx+FDH4JYLNiOxYJtBRgiItIbKbhI\ng9xcuP12mDoVKiqC59tvD/aLiPRmr73W+XYYRo8ezV133cWJJ57I4MGDmTlzJv/3f//HpEmTKCws\n5CMf+Qg7duwAYOXKlZx11lkUFRVx0kknsWLFipbzPPjgg7z3ve+lsLCQY489lvvuu6/lvRUrVjBq\n1Cj+67/+i+HDhzNy5EgefPDB8H+Z3iLTK/KFvLqfR9ns2e4QPIuIREHiezMl36/xuPspp7h/+9vB\n9ve/7/7e97o3NIRS9BZlZWV+5plnel1dnW/evNmHDRvmp5xyir/wwgv+zjvv+Hnnneff+ta3fNOm\nTT5kyBD/3e9+5+7uTz31lA8ZMsS3bt3q7u5LlizxV1991d3df//73/vAgQN99erV7u4ei8U8Ly/P\nKyoqvLGx0ZcsWeIDBw707du3h/vLZJHOPjvKXKRJLAY//jHMnh08NzeRiIj0Vjk58JvfwMMPQ0kJ\n/OhHsHQp5KVgbugvfelLDB06lCOOOIIJEyYwfvx4TjjhBPr168cnP/lJnn/+eX72s5/xsY99jI9+\nNFgL8/zzz+fUU09lyZIlAFx44YWUlZUBMGHCBD7ykY/whz/8oeUa/fr145vf/Ca5ublceOGFDBo0\niJdeein8X6YXUHCRBvE43HwzLFoUNIssWhRsq8+FiPR2RxwBU6bA1q3w4Q/D0Uen5jrDh++fmbyg\noOCA7d27d1NbW8vPf/5ziouLKS4upqioiGeffZYtW7YAsHTpUs4880yGDBlCUVERS5cuZevWrS3n\nGTJkCDk5+/9sDhw4kN27d6fmF8pyCi7SIDcX/vhHaF4gsbw82FafCxHp7X7wg+CG6s9/ht//Hm67\nLTPlMDOOPvpoLr/8crZt28a2bdt466232LVrFzfccAP79u3j4osv5oYbbqCuro633nqLCy+8sLlJ\nSA6Rgos0aRtIKLAQkd6uqQn++U9YvhxOPx2eeQa2bIHGxsyU57Of/SyLFy9m2bJlNDU1UV9fz4oV\nK9i8eTP79u1j3759DB06lJycHJYuXcqyZcsyU9BeQMGFiIikRE4O3H33/qaQI44I+pyF3efCzDrd\nbjZy5EgWL17Md7/7XUpKSigtLeWuu+6iqamJQYMGcffddzN16lSKi4tZuHAhU6ZMOaTryn5aFVVE\niMdbZ9PabkvvpVVRpbu0KmoEaIZOiSpN8iYiYct4cGFmE81sg5m9bGY3dnBMuZmtNrO/mdnydJex\np/TlLVGmSd5EJGwZbRYxsxzgZeB8YDOwCrjE3TckHXMY8L/AR9x9k5kNdfetHZwvsmm7WCz40r7m\nmqDNcdGi/aNHRKKgogLmzAnmYqmoyHRpJF3ULCLdFeVmkdOBV9y91t0bgIVA2x40lwK/cPdNAB0F\nFlFXXh4EFnPmBM8KLCRKNMmbiIQp08HFSGBj0vbriX3JjgOKzWy5ma0ys8vSVroQ6cs7XOrDEh5N\n8iYiYUvBJKyhywNOBs4D3gX8ycz+5O5/b+/giqR8bnl5OeURSBEkf3mXlwePm2/WRFrd1dyH5fbb\ng7qMxVSfPdE8yVtz3WmSt94tFosR092NpFim+1ycAVS4+8TE9k0EC6HcmXTMjcAAd5+T2L4fWOru\nv2jnfJFtE9RQv3CpD4tIONTnQroryn0uVgHHmlmpmfUDLgEWtznm18CHzCzXzAYC44H1aS5nj2mG\nznCpD4uISHRlNLhw9zgwC1gGVAML3X29mX3BzD6fOGYD8ASwFlgJ3OfuL2aqzBIN6sMiIl1RX1/P\n5MmTKSoq4jOf+QyPPPIIEydObHk/JyeHf/7znyktwxNPPMFFF12UknOnqvz33HMPN910U/dP0NFa\n7Nn4CH4d6e0aG93POMN9+fJge/nyYLuxMZOlEslOie/NXvv9+tBDD/n48eO9qamp3fdzcnL8H//4\nh7u7X3nllf7Nb34z9DKceuqp/pe//KVl28x80KBBPnjwYD/qqKP8+uuv77B8ycrLy72ysrLVvuTy\nh6m+vt6POuoor6ur6/CYzj47mW4WETlkubmwYkXrVWZXrFBTk4gcqLa2luOOO67DdUA8xH4k8XaG\nWD333HPs3LmT0047rWWfmbF27Vp27tzJ008/zSOPPMK8efO6dc0wy5+sf//+TJo0iQULFnTr5xVc\nSNaJx+Gcc1rPeHrOORo6KRJJ7rBkCXzve7B4cbAdsg0bNnDuuedSVFTE8ccfz29+8xsgGD34rW99\ni4ULF1JYWMgDDzzA/PnzmTBhwgHnmDdvHg8//DBz586lsLCwZdGyLVu2cPHFFzNs2DDGjBnDj370\no5afmTNnDlOnTuWyyy7j8MMPZ/78+Qecd+nSpZxzzjmt9vn+bBDHHXccEyZM4G9/+xt33XUXF198\ncatjr7vuOr7yla9w66238oc//IFZs2ZRWFjItdde23LMk08+yXHHHUdxcTGzZs1qdZ3bbruNsrIy\nRowYwZVXXsnOnTuBIOjKyclhwYIFlJaWMmzYML773e+2uvY555zDb3/724P/A7Sno5RGNj7I0rSd\nHLrly92HDnWfPTt4bm4iEZFDQ6qbRa691v1d73LPzw+eP/e5MIrdoqGhwY899li/4447vKGhwZ95\n5hkfPHiwv/zyy+7uXlFR4ZdddlnL8Q8++KBPmDChZdvMOmwWaWpq8lNOOcVvu+02b2xs9FdffdXH\njBnjy5Ytazl3v379fPHixe4eNCW0NXXqVL/rrrta7Uu+ZnV1tY8YMcIfeOAB37Jliw8aNMh37Njh\n7u6NjY0+bNgwX716tbu33yxiZj558mTfuXOnv/baa15SUuJPPPGEu7tXVlb62LFjvaamxvfs2eMX\nXXRRS13U1NS4mfnnP/95f+edd/yFF17w/v37+4YNG1rO/fzzz/uQIUM6rPvOPjvKXEhW0mgRkSyw\naRPcdx/s2QMNDcFzVRW8/HJol1i5ciV79uzhxhtvJC8vj3PPPZePf/zjVFVV9fjcq1atYuvWrdxy\nyy3k5uZSVlbG1VdfzcKFC1uOOfPMM5k8eTIQNCW0tX37dgYPHnzA/pNPPpkhQ4YwZcoUPv/5z3Pl\nlVcyYsQIzj77bBYtWgQEWY+SkhI+8IEPdFrOm2++mcGDBzNq1CjOPfdc1qxZA8AjjzzC9ddfT2lp\nKQMHDuT2229n4cKFNDU1AUHzTEVFBf369eOEE07gxBNP5IUXXmg57+DBg9mxY8ch1logGybREjlA\n29EizZOTiUiEbNsG+flQX79/X34+/OtfoV1i8+bNjBo1qtW+0tJSNm3a1ONz19bWsmnTJoqLi4Eg\n09/U1MTZZ5/dckzba7dVVFTErl27Dti/evVqRo8efcD+yy+/nHvvvZcZM2bw8MMPc9llB5+Uevjw\n4S2vBw4cyO7du4GgbkpLS1veKy0tpbGxkTfffPOgPwuwa9cuDjvssINevz3KXKSJpqsOj6arFskS\nY8dCQQEkd6bMy4P3vz+0Sxx55JFs3Lix1b7XXnuNkSPbriRxcG07fY4aNYpjjjmGbdu2sW3bNt56\n6y127NjR0qejvZ9p64QTTuDldjI13kHfk0984hOsXbuW6upqHn/8caZPn97la7V15JFHUltb27Jd\nW1tLfn5+q4CiM+vXr+fEE088pGs2U3CRBlpyPVzN01UnjxbRdNUiETRgQPCFN25ckLEYOxaeeQba\naSborvHjxzNw4EDmzp1LY2MjsViMxx9/nGnTph3yuYYPH95qzojTTz+dwYMHM3fuXOrr64nH41RX\nV/Pcc891+ZyTJk06pOnW+/fvz6c+9SkuvfRSxo8fz1FHHdVh+Q5m2rRpfP/736empobdu3dzyy23\ncMkll5CTE/zp7yjAabZixQouvPDCLl8vmYKLNMjNDdbBmDo1uNOeOjXY1h/D7tOMp+FSZk1SZtw4\nqK6GffuCvhbdvBPuSH5+Pr/5zW9YsmQJQ4cOZdasWTz00EOMHTu2Sz+fnA2YMWMG1dXVFBcXc9FF\nF5GTk8Pjjz/OmjVrGD16NMOGDWPmzJktIy664qSTTuLwww9n1apV7V6zPVdccQXr1q3j8ssvb7X/\nuuuuY9GiRQwZMoQvf/nL7Z4reftzn/scl112GWeffTZjxoxh4MCB3H333R2WI3m7vr6eJUuWcMUV\nV3TxN20to2uLhC3qc99XVAQdEGfPDl6LRIEWguvbtLZI6j355JP85Cc/4Ze//GWXjt+4cSPjxo3j\njTfeYNCgQSkuXfvuueceXn/9de64444Oj+nss6PgIk200JZEmT6ffZeCi2hpamri+uuvZ/fu3dx/\n//2ZLk6nOvvsaLRIGmjJdYm65KG9s2crsBDJhLfffpvhw4czevRoli5dmuni9IgyF2miJdclypS5\n6LuUuZDuUrOIiHRIfS76NgUX0l0KLkSkU8qs9V0KLqS7OvvsaCiqiGhor4iESsGFiIiIhErBhYiI\niIRKwYWIiPRa9fX1TJ48maKiIj7zmc/wyCOPMHHixJb3c3JyDmlK7e544oknuOiii1J6jfasW7eO\ns846K+3XBQUXIiLSiz322GPU1dWxbds2Hn30US699FJ+97vftbyfPOX1VVddxX/8x3+EXoZbb72V\nm2++uWU7OaCZM2cO/fr1o7CwkOLiYj70oQ+xcuXKA85RU1NDbm4u11xzzQHv5eTkMHjwYAoLCxk8\neHDLKq7HH388RUVF/Pa3vw39dzqYjAcXZjbRzDaY2ctmdmMnx51mZg1mlv7wT0REslJtbS3HHXdc\nh+t5hDkCJt7OojzPPfccO3fu5LTTTmvZ17Ysl1xyCTt37mTr1q2Ul5czderUA86zYMECiouLefTR\nR2loaGj1npmxdu1adu7cya5du9i2bVvLe5deein33ntvT3+1Q5bR4MLMcoB7gI8C7wOmmdl7Ojju\nDuCJ9JYwPFoYSqTv0P/3/dydJUuW8L3vfY/FixeH+se82YYNGzj33HMpKiri+OOPb1kSvaKigm99\n61ssXLiQwsJCHnjgAebPn8+ECRMOOMe8efN4+OGHmTt3LoWFhUyZMgWALVu2cPHFFzNs2DDGjBnD\nj370o5afmTNnDlOnTuWyyy7j8MMPZ/78+Qecd+nSpZxzzjmt9nVUBzk5OUyfPp3Nmzfzr3/9q9V7\nCxYs4LbbbmtZqK3t+To6Z3l5OU8//fQBAUmqZTpzcTrwirvXunsDsBCY0s5xXwIeA/4vnYULi5Zc\nF+k79P+9tS9/+ct8+tOf5pZbbuHSSy/l6quvDvX8jY2NTJ48mYkTJ1JXV8fdd9/N9OnTeeWVV6io\nqOAb3/hGS2bgqquuAtpflXTmzJlMnz6dG264gZ07d/LrX/8ad2fy5MmcdNJJbNmyhaeffpof/vCH\nPPnkky0/t3jxYj796U+zfft2pk+ffsB5161bx7vf/e4u/S779u1j/vz5DBkyhKKiopb9f/jDH9i0\naROXXHIJU6dObTeI6ciRRx5Jfn4+L730Upd/JgyZDi5GAhuTtl9P7GthZkcCn3D3nwAHneglirTk\nukjfof/v+23atIn77ruPPXv20NDQwJ49e6iqquLll18O7RorV65kz5493HjjjeTl5XHuuefy8Y9/\nnKqqqh6fe9WqVWzdupVbbrmF3NxcysrKuPrqq1m4cGHLMWeeeSaTJ08GoH///gecY/v27QwePLjT\n6zz66KMUFxczcOBAKisreeyxx8jJ2f/necGCBUyaNInDDjuspc/I1q1bW53j5JNPpqioiOLi4pbl\n2JsNHjyY7du3H/Lv3xPZsHDZD4DkvhidBhgVSWuZl5eXUx6RBRK0MFS4NKOkRFmU/7/HYjFizWmV\nFNu2bRv5+fnU19e37MvPzz8g5d8TmzdvZtSoUa32lZaWsmnTph6fu7a2lk2bNrV0kHR3mpqaOPvs\ns1uOaXvttoqKiti1a1enx3zmM59hwYIFbNu2jU996lM899xzLdeor69n0aJFVFZWAnDGGWcwatQo\nHnnkEa699tqWc6xevZrRo0e3e/5du3Zx+OGHH/wXDlGmg4tNwNFJ20cl9iU7FVhoQR5rKHChmTW4\n++L2TpgcXERJLBYsCDV7dvDcvDqqHDqthSFR9/TTrf+/T5gA55+f6VIF2t50zZkzJ2XXGjt2LAUF\nBezevbulT0BeXh7vf//7Q7vGkUceycaNG1vte+2117rcFJGsbXPJqFGjOOaYYzptUuioo2izE044\nocuZmuLiYn76059y6qmnMn36dIYPH84vf/lLdu7cyb//+78za9YsAHbs2MH8+fNbBRcd9bnYvHkz\nDQ0N3aqPnsh0s8gq4FgzKzWzfsAlQKugwd2PSTxGE/S7+PeOAouoSl5yvaIieL755r7bBttTSjtL\nlO3bB5/4BNxyS/D5vOWWYHvfvkyXLP0GDBhALBZj3Lhx5OfnM3bsWJ555pmDNhMcivHjxzNw4EDm\nzp1LY2MjsViMxx9/nGnTph3yuYYPH95qzovTTz+dwYMHM3fuXOrr64nH41RXV/Pcc891+ZyTJk06\npEzRcccdx8SJE7nzzjsBmD9/PjNmzGDdunW88MILvPDCC/zxj39kzZo1VFdXH/R8K1as4LzzziM/\nP7/LZQhDRoMLd48Ds4BlQDWw0N3Xm9kXzOzz7f1IWgsYktzc4K66+WahvFx32T2VnHa+5hplgXpK\noxvC068f/OpX8J3vBMHFd74TbPfrl+mSZca4ceOorq5m3759vPzyy5x44omhnr959MSSJUsYOnQo\ns2bN4qGHHmLs2LFd+vnkzMOMGTOorq6muLiYiy66iJycHB5//HHWrFnD6NGjGTZsGDNnzmTnzp1d\nLt9JJ53E4YcfzqpVq9q9Znu+9rWvMW/ePF577TWWL1/OV77yFYYNG9byOPnkk7nwwgtbOnZ2dr6H\nH36YL37xi10ub2iah7D0hkfw60hfsHy5+9Ch7rNnB8/Ll2e4QFmssdH9jDP21+Hy5cF2Y2MmS5X9\nZs92h+A5yhLfm/p+TaFly5b5Jz/5ybRfd+3atf7BD34wZefv7LOjJdcl66jPRfhisaB56Zprgj4C\nixYpG9QT2VSfWnJduquzz46CC8lKGi0SvoqK/aMbItovOitkW/Cr4EK6S8GFiHQqm+60s0E2Bb8K\nLqS7FFyISIey7U5bwqXgQrpLwYWIdCqb7rQlXAoupLs6++xkep6LPkND/STK2gYSCixEpCcUXKSB\nFjISkWw3YMCAN80MPfRofgwYMODNjj4vahZJE3WYk8irq4OaGigrg5KSTJdG0sSsa80iIodCmYs0\n0YySEmVND1dBaSlccAGUlgbbIiLdpMxFmihzIVEVf6OOhpGlDGja27KvPqeA/E215I5QBqO3U+ZC\nUkGZizTQwmXhUwfZ8ORurCGvoPXCF3kF+eRurMlMgUQk6ym4SAMtXBYudZANWVkZeU2tl+zMa2oI\n+l6IiHSDgos00VC/8GjJ9ZCVlFD91UrepoD6/oW8TQHVX61Up04R6TYFF5KV1EE2PPE4zHhyGqv/\np5YBf3iK1f9Ty4wnpykTJCLdpuBCslIsFnSMnT07eG5uIpHuazi8BE47LXgWEekBjRaRrKO1MMKn\n0Ux9l0aLSCoouJCspLUwwqcl1/smBReSCmoWkayUm0swo+SqVVBXp8Cih9TMJCJhUnAh2amq9YyS\nVGlGye7SPCwiEraMN4uY2UTgBwSBTqW739nm/UuBGxObu4D/5+7rOjiXmkX6grq6IKDYu39GSQoK\noLZWwye7Sc1MfZeaRSQVMpq5MLMc4B7go8D7gGlm9p42h/0TONvdTwRuA+alt5QSOTU10K/1jJLk\n5wf7pVs0D4uIhCnTzSKnA6+4e627NwALgSnJB7j7SnffkdhcCYxMcxklasrK4O23W+/bu1czSoqI\nRESmg4uRwMak7dfpPHi4Glia0hJJdjDrfFtERDImL9MF6CozOxe4CvhQZ8dVJI2hKy8vp1yD9Xuf\nmpqgj8W+pPUwBgwI9qvPhUinYrEYMQ0HkhTLaIdOMzsDqHD3iYntmwBvp1PnCcAvgInu/o9OzqcO\nnX2BOnSmRl1dEKCVlake+xB16JRUyHSzyCrgWDMrNbN+wCXA4uQDzOxogsDiss4CC+lDSkqgsjII\nKAoLg+dKLbTVIxraKyIh6lLmwsw+CJSR1Izi7gtCKUAwFPWH7B+KeoeZfSG4hN9nZvOAi4BawIAG\ndz+9g3NFNnOhoX7hischd9v+O+14cYnqs7uUCerTlLmQVDho5sLMHgLuIujrcFricWpYBXD337n7\nu919rLvfkdj3U3e/L/F6prsPcfeT3f2kjgKLKGteC6O5mTMWC7Y1SVH3tNRndbDQVqy6RPXZExra\nmxpJM8iK9DVd6dB5KvDeyKYEskBubrDIVtuFoXSn3T2qz5CVlbXuHAvQ0KChvT1RVQUzZgRB2759\nQbPdtGmZLpVI2nSlz8XfgBGpLkhvV14e/CGcMyd41iCWnlF9hkh9WMJVV4fPmBE0M+3YAXv3BtvK\nYEgf0pXgYijwopk9YWaLmx+pLlhvo4WhwqX6DNm0aUEfi6eeCp51l91t8X/UsGdf62amPfvyif+j\nJjMFEsmAg3boNLNz2tvv7itSUqIeiGqHzuY+AnfdWMdZI2t4dlMZX7uzhD/+Uan87miuz9tvDzIW\nsViw0JbqUyKhro74qFJy39nfQTbev4DcjdHsIKsOnZIKGV+4LExRDS4Amh6uImfm/jbYpnmV5EzX\n3WF3afSNRFpVFQ2Xz+DtxnwG5jWQvyC6fS4UXEgqdGW0yBlmtsrMdpvZPjOLm9nOdBSu16irCwKL\npDbYnJlqg+0JLbQVrrYjbTTypmdiR0zj+MJaHr36KY4vrCV2RDQDC5FU6Uqfi3uAacArQAHB+h4/\nTmWhep2OhvRpqJ9EgIZKhyseD5rp7v1FCZ+fdxr3/qKEm29WfUrf0pU+F8+5+6lmttbdT0jsW+3u\nJ6WlhIcgss0i69fj730vyXlHB+zFF2HcuEyVKqupWSRcsdiBQ3s1Aqf7sunzqWYRSYWuZC7eTkzN\nvcbM5prZV7r4c5IQ37Gbd6yg1b53bADxHbszVKLs1nyn/eyvgkmKnv1Vne60e0hDe0UkTF0JEi4D\ncoFZwB5gFPCpVBaqt8kdU0Z+2wkQ+xm5Y8oyUp5sl5sL959fxUmfLKV+wgWc9MlS7v9wVWTvDLOB\nhvaGR81MIhotkj5Z1Hs88rQWRqg0tDd82dTMpGYRSYUOMxdm9vPE8zozW9v2kb4i9g7qPR6imhoa\nc1qnghpztBZGd+XmBoFE8x+/8nIFFj2lZibp6zrMXJjZEe6+xcxK23vf3WtTWrJuiGrmQpNohSv+\nRh0NI0sZ0LQ/c1GfU0D+plpyRyhzIZmnzIX0dWoWSRNNohWulvrMz4eGBtWnREa2NTMpuJBU6Cxz\nsYtgxGS73L0wVYXqrsgGF+ojkBp1dUFTSFmZ6lEiJR6H3G37P5/x4pJIBhag4EJSo8Ml1919MICZ\nfRvYAjwEGDAdOCItpestamqCjEVycJGf6COgP4rdV1Ki+pNIyv156yXXc7XkuvQxXZlE6wV3P/Fg\n+6JAmQsRybi6Ory0FEv6/+4FBVhE/78rcyGp0JV5LvaY2XQzyzWzHDObTjDfhXRVSQlUVgYBRWFh\n8FxZGckvmmyhtTDCFY8TBMGrVgWreqo+u01Lrot0Lbi4FPg08GbiMTWxTw7FtGlBpuKpp4JnpUi7\nTZMUhSsehznvqSI+qhQuuID4qFIq3lOl+uym3DFlFOTta7WvIK9Bk+ZJn5Lx0SJmNhH4AUGgU+nu\nd7ZzzN3AhQQZkyvdfU0H54pms4iELpuG+kVeXR3xUaXkvrM/jR/vX0Duxmim8bNCFk2ap2YRSYWu\nLLl+nJk9bWZ/S2yfYGa3hnFxM8shWHX1o8D7gGlm9p42x1wIjHH3scAXgHvDuLZkN01SFKKaGnIH\ntE7j5/bXpGQ9oUnzpK/rSrPIPOBmoAHA3dcCl4R0/dOBV9y91t0bgIXAlDbHTAEWJK79Z+AwMxse\n0vUlS2ktjBCVlRGvb53Gj7/TEAzxlUOmJddFuhZcDHT3v7TZ1xjS9UcCG5O2X0/s6+yYTe0cI31I\n85f3okUtN4RJAAAW7ElEQVRQURE868u7++LFJVSMqiTeP+hwHO9fEGwXq0mkOzSdukjXgoutZjaG\nxIRaZnYxwbwXkWRmBzwqKiraPbaioiIzx3/969EqT5Ydn5dnrFxpnHtucHzst19v98s7quWP2vF5\necZtf7+UvHf2Yjt38u0vXUPFhmmqzx4cn5fX+vj2AotMlf/KK6+koqKi5SGSCl2Z5+IY4D7gg8Bb\nwKvAZ929pscXNzsDqHD3iYntmwBP7tRpZvcCy9390cT2BuAcd3+znfNFt0NnVetJddCkOj2j+gyf\nZjztk8zUoVPC1+XRImb2LiDH3XeFdnGzXOAl4HyCbMhfgGnuvj7pmEnANe7+sUQw8gN3P6OD80Uz\nuNAkWuFSfYZPwVqfpeBCUqHD6b+bmdnhwOVAGZBnFnwG3f3anl7c3eNmNgtYxv6hqOvN7AvB236f\nuy8xs0lm9neCoahX9fS6aVdTg+f3az1jX34+pum/u0fTqYerri4ILPbu3V+nM2bAhz+s+uymeLx1\nM13bbZHerit9LpYQBBbrgL8mPULh7r9z93e7+1h3vyOx76fufl/SMbPc/Vh3P9Hdnw/r2ukSH1VG\nw669rfY17KonPqosMwXKdmVlwd11sgaNbui25mAtWb6GonaXJnkT6ULmAhjg7tenvCS9WG4uWI5D\n0pdLXo6TozuZ7mmeTn3G/iXXNZ16DyhYC1VubrDcettJ3pS5kL6kK5mLh8xsppkdYWbFzY+Ul6w3\nqakhZ9DAVrty3lWgO8Oe0HTq4dHaN6HTJG/S13VltMg1wHeA7SSGoxL0hzgmxWU7ZFHu0KnplSXy\n1q+Hv/wFTj8dxo3LdGmyWiwGX/xUHddfVMN//bKMe39REtkAQx06JRW6krn4KnCsu5e5++jEI3KB\nRZRpkiKJuqaHq+CUU+C66+CUU4Jt6ZZ4HJ6ZWUX1nlI+v+gCqveU8vRMLQQnfUtXMhfLgE+4+9vp\nKVL3RTZzQaK3+Lb98wjEi0vUBiuREH+jjoaRpQxo2p9Zq88pIH9TLbkjFAAfsro6vLS09eiwggIs\nokOllbmQVOhKh849wBozWw6807wzjKGofUluLsEXS+LLRXGFREXuxhq8oB/s2f/HMK8gn9yNNaDg\n4tDV1GBthkqbhkpLH9OV4OJXiYdIZGgegRCVlZHX1Hq0SF6TRot0m0bfiBy8z4W7z2/vkY7CibSn\neR6BZ39VB6tW8eyv6jSPQE+UlFD91UrepoD6/oW8TQHVX9VokW7T6BuRrk//nQ2i3OcC0NoNIaq+\ntYrR35lBTv9+NL2zj1dvreR939Zw1O5oDtbuurGOs0bW8OymMr52Z4lW8uypLPn/rj4XkgoKLtKk\n6eEqcmbuX7uhaV4lOdP1x7BbtLZI6NTM1HcpuJBU6MpQ1BZmNiJVBenN4m/Use/yxNoNO3bA3r3s\nu3wG8TfqMl207FRTQ2NO6+mqG3M0XXVPtA0kFFiISE8cUnBBsM6IHKLcjTXkFbT+Y9jSG18OWXxU\nGY17W3eYa9zboLVaREQi4lCDC6XOukO98UOVO6KEfgtad5jrt6BSczKIiETEoQYX81JSit5OvfFD\nlzO99doi6r8iIhId6tCZBuqNLyJRpQ6dkgoKLtJEvfFFJIoUXEgqHGqziHSTeuOnQF0wiRZ1GnUT\nCtWniISkw+DCzI7u5L0JqSmOSBdVVQVzXVxwQfBcpVU8e0T1KSIh6rBZxMz+CdwL/Ke7xxP7hgP/\nCbzH3U9NWym7KMrNIhIiTaIVLtVnn6ZmEUmFzppFTgHGEKyIep6ZXQf8BfgTcHpPL2xmRWa2zMxe\nMrMnzOywdo45ysyeMbNqM1tnZlqJVYLJsvq1njeEfE2i1W2qTxEJ2UE7dCaCiu8Dm4Ez3P31UC5s\ndifwL3efa2Y3AkXuflObY0YAI9x9jZkNAv4KTHH3DR2cU5mLvkB32uFSfYYumzpwK3MhqdBZn4vD\nzeynwFXAROAxYKmZnRfStacAzaurzgc+0fYAd3/D3dckXu8G1gMjQ7q+ZCutOhku1Weomoeex2LB\ndiyGVu2VPudgfS7+P+AH7t6Y2PeBxL5ad+/RrEVmts3dizvabuf4MiAGvD8RaLR3TLQzF1mySmLW\nUH2GS/UZmlgMpk6Fa66BH/8YFi2C8vJMl6p9ylxIKuR18t7ZbZtAElmED5rZzK6c3MyeBIYn7wIc\nuLWdwzuMChJNIo8B13UUWDSrqKhoeV1eXk55VP5HV1XBjP2rolJZCdM0q2SPlJToj2BI4nHITarP\nKKfxs0F5eRBYzJkDs2dHK7CIxWLEmtMqIimSsUm0zGw9UO7ubyb6Vix393HtHJcHPA4sdfcfHuSc\n0cxc1NXhpaVYUpu2FxRgatOWCGhO499+e/BHMBaDm29GM8j2gDIX0tdlchKtxcCViddXAL/u4Lj/\nBl48WGARZfF/1LBnX+ve+Hv25RP/R01mCtQLtG2/Vnt29+XmBoHF1KlQURE83367AovuiseD4GzR\noqA+Fy0KtvUZlb4kk8HFncAFZvYScD5wB4CZHWFmjydenwVMB84zs9Vm9ryZTcxYibspd0wZBXmt\nV0UtyGsgd0xZZgqU5dRhLnzJafxrronuXXY2yM0Nsj7l7wtmPC1/X52yQNLnaG2RdKmqYvNln+Ol\neC7vzo1z5EP/rT4XPZBNaedsoPoMWVUVfO5zQUQRj8N/R/f/u5pFJBW0tkia3FoNR8WNSfnB8zdf\nzHSJspvutMOjNH7I6urgiiugvh727Amer7hCa7ZIn6LMRRq88UYdI0eW0tS0v0NnTk4BmzbVMmKE\nOnR2h+60w5VNkz5F3rJl8NGPHrj/iSfgIx9Jf3kOQpkLSQVlLtJg48YaBg1q3aFz0KB8Nm6syUyB\nspzutMOnVXvD1fYWJ3q3PCKppcxFGtTV1VFaWsrepKGoBQUF1NbWUqKhqN0Sj0Putv2TPsWLS/QH\nUSIh/kYd8SOPop/v78S9z/qRu/l1ciOYqVTmQlJBmYs0KCkpobKykoKCAgoLCykoKKCyslKBRQ/k\n/rz1EuG5P9cS4RINuSNKeOWWB9lLAe/kv4u9FPDKLQ9GMrAQSRVlLtKorq6OmpoaysrKFFj0hBba\nCp36XIRv7tfrWHRXDVO/VsYN34vu51KZC0kFZS7SJB4PMhinnXYaJSUl6h/QE1oiPFTN84Y8+6tg\nXoZnf1WneUN6KBaD7z1Ywsdmn8b3HixBs21LX6PgIg006VPIysqC9VmSNTQE++WQ5ebC/edXcdIn\nS6mfcAEnfbKU+z9cpcxFN6nDsYiaRdJGQydD1rwQXH5+EFhoIbjuUzNT6LKpmUnNIpIKna2KKiGK\n8iqJWWnaNPjwh7VEeBhqamjM6Uce+4OLxpx88mpqVK/dpKG90tepWSRNYrEgYzF7dvCsNtgQlJTA\naafpD2APxUeV0bi3dTNT494G4qPKMlMgEcl6Ci7SQG2wEmW5I0rot6AyaAopLISCAvotqNTQSRHp\nNvW5SJNsaoOVPqquTs1MYcqS+lSfC0kFBRfplCVfNiLSQ80djvv1C0Y2RbjDsYILSQUFF+mSRV82\nWUPBmkRRlo2+UXAhqaA+F+lQVxcEFnv3wo4dwfOMGVqCuSeqWk//TZWm/5aI0CRvIgou0kJfNuFS\nsJYadcEMnarHHtIkbyIKLtJCXzbhUrAWPmWCwlNSEjR7Jo2+obIykk0iIqmiPhfpohklw5NlbdqR\np/pMjSzpE6Q+F5IKGctcmFmRmS0zs5fM7AkzO6yTY3PM7HkzW5zOMoZq2rTgy/qpp4JnBRbdpzvD\ncCkTlBqa5E36sIxlLszsTuBf7j7XzG4Eitz9pg6O/QpwClDo7v/WyTmjm7mQ8GXJnWHkKXPRpylz\nIamQyT4XU4D5idfzgU+0d5CZHQVMAu5PU7kkW+jOMBzKBIlIyDKZudjm7sUdbSftXwR8BzgM+Koy\nFyIpokxQn6TMhaRCSldFNbMngeHJuwAHbm3n8AOiAjP7GPCmu68xs/LEz3eqoqKi5XV5eTnlUVp+\nVF/eEmUlJfpc9gGxWIyYVk6UFMtk5mI9UO7ub5rZCGC5u49rc8x3gc8CjUABMBj4pbtf3sE5o5u5\n0AydIhJBylxIKmS6Q+c2d7/zYB06E8efQ7Y2i6jDnIhElIILSYVMdui8E7jAzF4CzgfuADCzI8zs\n8QyWK3wa6iciIn2IJtFKB2UuRCSilLmQVND03+mgoX4ifY/WapE+TJmLdNJoEZG+IYs6cCtzIamg\n4EJEJExZ1gyq4EJSQc0iIiJhUgduEQUXIiKhKisLmkKSNTQE+0X6CAUXIiJhUgduEfW5EBFJiSzp\nwK0+F5IKCi5ERPowBReSCmoWERERkVApuBAREZFQKbiQ7KUZEMOl+hSRkCi4kOxUVRVMVHTBBcFz\nVVWmS5TdVJ8iEiJ16JTsk2UzIEae6jM1NFpE+jBlLiT7aAbEcKk+w6dMkPRxylykU5bcyUSe7rTD\npfoMV5bVpzIXkgrKXKSL7mTCoxkQw6X6DJcyQSLKXKRFlt3JZA1lgsKl+gxHlv1/V+ZCUkGZi3TQ\nnUxqlJTAaadF8gs7K6k+w6FMkIgyF2mRZXcyIhKCLMkEKXMhqZCxzIWZFZnZMjN7ycyeMLPDOjju\nMDNbZGbrzazazManu6w9pjsZkb5HmSDpwzKWuTCzO4F/uftcM7sRKHL3m9o57kFghbs/YGZ5wEB3\n39nBOaOZuWiWJXcyItJ3KHMhqZDJ4GIDcI67v2lmI4CYu7+nzTGFwGp3H9PFc0Y7uBARiRgFF5IK\nmezQOczd3wRw9zeAYe0cMxrYamYPmNnzZnafmRWktZQiIiJySPJSeXIzexIYnrwLcODWdg5vL+WQ\nB5wMXOPuz5nZD4CbgNkdXbOioqLldXl5OeXl5Ydc7pRRs4iIZFgsFiMWi2W6GNLLZbJZZD1QntQs\nstzdx7U5ZjjwJ3c/JrH9IeBGd5/cwTmj2yxSVQUzZgRDUvftCzp0TpuW6VKJ7Kfgt09Ss4ikQiab\nRRYDVyZeXwH8uu0BiWaTjWZ2XGLX+cCLaSldmOrqgsBi717YsSN4njFDS1tLdGgGWREJUSYzF8XA\nz4FRQC3waXffbmZHAPPc/eOJ404E7gfygX8CV7n7jg7OGc3MxapVwZf2jqRiFxbCU08FQ9VEMknz\nsPRpylxIKqS0z0Vn3H0b8OF29m8BPp60/QKQ3X+By8qCppBkDQ3BfpFMa55BNjm4aJ5BVsGFiHSD\npv9OB02iJVGm4FdEQqbpv9NJHeYkqpo7HOfnB4GFOhz3GWoWkVRQcCEiAQW/fZKCC0kFBRciIn2Y\nggtJBfW5EBERkVApuBAREZFQKbgQERGRUCm4EBERkVApuBAREZFQKbgQERGRUCm4EBERkVApuBAR\nEZFQKbgQERGRUCm4EBERkVApuBAREZFQKbgQERGRUCm4EBERkVApuBAREZFQZSy4MLMiM1tmZi+Z\n2RNmdlgHx33FzP5mZmvN7GEz65fusoqIiEjXZTJzcRPwlLu/G3gGuLntAWZ2JPAl4GR3PwHIAy5J\naylDFovFMl2ELlE5w6VyhkvlFIm2TAYXU4D5idfzgU90cFwu8C4zywMGApvTULaUyZYvG5UzXCpn\nuFROkWjLZHAxzN3fBHD3N4BhbQ9w983AfwKvAZuA7e7+VFpLKSIiIockL5UnN7MngeHJuwAHbm3n\ncG/n5w8nyHCUAjuAx8zsUnd/JAXFFRERkRCY+wF/09NzYbP1QLm7v2lmI4Dl7j6uzTEXAx9195mJ\n7cuA8e4+q4NzZuaXERHJYu5umS6D9C4pzVwcxGLgSuBO4Arg1+0c8xpwhpkNAN4BzgdWdXRC/QcR\nERHJvExmLoqBnwOjgFrg0+6+3cyOAOa5+8cTx80mGCHSAKwGrnb3howUWkRERA4qY8GFiIiI9E5Z\nN0OnmU00sw1m9rKZ3djO++eY2XYzez7xaK/zaDrKWWlmb5rZ2k6OudvMXjGzNWb2gXSWL3H9TssY\nobo8ysyeMbNqM1tnZtd2cFym6/Og5YxCnZpZfzP7s5mtTpRzdgfHZaw+u1LGKNRlUllyEmVY3MH7\nGf1sJpWjw3JGqT6lF3D3rHkQBEN/Jxg9kg+sAd7T5phzgMURKOuHgA8Aazt4/0Lgt4nX44GVESxj\nVOpyBPCBxOtBwEvt/LtHoT67Us6o1OnAxHMusBI4PYL1ebAyRqIuE2X5CvCz9soThbrsYjkjU596\nZP8j2zIXpwOvuHutB/0uFhIMVW0r4x073f2PwFudHDIFWJA49s/AYWY2vJPjQ9eFMkI06vINd1+T\neL0bWA+MbHNYFOqzK+WEaNTp24mX/Qk6drdtH41CfR6sjBCBujSzo4BJwP0dHJLxuoQulRMiUJ/S\nO2RbcDES2Ji0/Trtf3mfmUg//tbM3pueoh2ytr/LJtr/XTItUnVpZmUE2ZY/t3krUvXZSTkhAnWa\nSI+vBt4AnnT3tqOwMl6fXSgjRKAuge8DX6f94AciUJcJBysnRKM+pRfItuCiK/4KHO3uHwDuAX6V\n4fJks0jVpZkNAh4DrktkBiLpIOWMRJ26e5O7nwQcBYyP4h+SLpQx43VpZh8D3kxkrIyI3vl3sZwZ\nr0/pPbItuNgEHJ20fVRiXwt3392cTnX3pUB+Ythr1GwiGIbb7IDfJdOiVJeJtWUeAx5y9/bmRIlE\nfR6snFGq00QZdgLLgYlt3opEfULHZYxIXZ4F/JuZ/ROoAs41swVtjolCXR60nBGpT+klsi24WAUc\na2alFiy9fgnBZFwtktsyzex0guG229JbzP3FoeM7mcXA5QBmdgbBuilvpqtgSTosY8Tq8r+BF939\nhx28H5X67LScUahTMxtqZoclXhcAFwAb2hyW0frsShmjUJfu/g13P9rdjyH4PnrG3S9vc1jGP5td\nKWcU6lN6j0zO0HnI3D1uZrOAZQSBUaW7rzezLwRv+33AxWb2/wgm3doLfCYTZTWzR4ByYIiZvQbM\nBvo1l9Pdl5jZJDP7O7AHuCpqZSQ6dXkWMB1Yl2iDd+AbBKOGolSfBy0n0ajTI4D5ZpZD8P/o0UT9\ntfw/ikB9HrSMRKMu2xWxuuxQttSnZB9NoiUiIiKhyrZmEREREYk4BRciIiISKgUXIiIiEioFFyIi\nIhIqBRciIiISKgUXIiIiEioFFyJpYmbLzezkxOvHzaww02USEUmFrJpES6S3cPePZ7oMIiKposyF\nSCcSU82vN7MHzOwlM/uZmZ1vZn9MbJ9qZgPNrNLMVprZX83s3xI/O8DMqsys2sx+CQxIOu+rzes2\nmNn/mNkqM1tnZlcnHbPLzG5LrFL5v2ZWkvYKEBHpBmUuRA5uDPApd3/RzJ4Dprn7h8xsMnAL8CLw\ntLvPSKyH8RczexL4IrDH3d9nZscDzyedM3lq3KvcfbuZDQBWmdkv3P0t4F3A/7r7rWZ2JzAT+G7q\nf10RkZ5R5kLk4F519xcTr6uBpxOv/waUAR8BbkqsJxIjWJ/laOBs4GcA7r4OeCHpnMmLxX3ZzNYA\nKwlWzByb2P+Ouy9JvP5r4loiIpGnzIXIwb2T9LopabuJ4P9QI0Fm45XkHzI7YLHZA3eYnQOcB4x3\n93fMbDn7m08akg6No/+vIpIllLkQObh2l6RP8gRwbcvBZh9IvPw9wSqpmNn7gRPa+dnDgLcSgcV7\ngDMO4boiIpGk4ELk4LyD183b3wbyzWytma0DvpV47yfAIDOrBiqA59o5z+8SP1tN0J/iT51cS0Qk\nK2jJdREREQmVMhciIiISKgUXIiIiEioFFyIiIhIqBRciIiISKgUXIiIiEioFFyIiIhIqBRciIiIS\nKgUXIiIiEqr/H7KUT22MH00nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3481f65668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(sky_med, sky_mean - sky_med, color='b', marker='x', label='mean')\n",
    "plt.scatter(sky_med, sky_ofil - sky_med, color='r', marker='o', label='ofilter (Python)')\n",
    "plt.scatter([1], [0.7501699 - 1], color='k', marker='o', label='ofilter (IRAF)')\n",
    "plt.xlabel('median')\n",
    "plt.ylabel('X - median')\n",
    "plt.axhline(0, color='k', linestyle='--')\n",
    "plt.legend(bbox_to_anchor=(1.01, 1), loc='upper left', scatterpoints=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f3481e39d68>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAEPCAYAAADxtOYjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4FPXZN/DvnZBIMAc2EJLIISgRqrTxCBSVEh6rcrJU\nAtZAVKyH53qKh1d9q/JUTbS+gj5QldrWqghYDhGV9zJaUVQI1VYUikhNAvhWCQmEkEBMiCbmdL9/\nbLKbw+4mu5nZ2cP3c1177e7s7My9dyYz9/7mt78RVQURERGRUSKsDoCIiIhCC4sLIiIiMhSLCyIi\nIjIUiwsiIiIyFIsLIiIiMhSLCyIiIjKU5cWFiEwXkf0iclBE7vcw3wQRaRaRuf6Mj4iIiLxjaXEh\nIhEAngVwFYDxALJF5Adu5lsG4F3/RkhERETesrrlYiKAL1W1VFWbAeQDmONivjsAvAbguD+DIyIi\nIu9ZXVwMB1DW6Xl5+zQHETkDwM9V9U8AxI+xERERkQ+sLi764mkAnftisMAgIiIKYAMsXv8RAKM6\nPR/RPq2ziwHki4gAGApghog0q2pB94WJCC+UQkTkJVXllzYylNUtF7sApItImohEA7gOQJeiQVXP\nar+dCXu/i1+5Kiw6zc+bKnJzcy2PIRBuzANzwVx4vhGZwdKWC1VtFZHbAWyFvdBZpaolIvKf9pf1\n+e5v8XuQRERE5BWrT4tAVd8BMK7btD+7mfeXfgmKiIiIfGb1aREySWZmptUhBATmwYm5cGIuiMwl\noXTOTUQ0lD4PEZHZRATKDp1kMMtPixARUeCLiYk51tjYmGx1HBQ4Bg4cWNnQ0JDi6jW2XBARhbG+\ntlxw/0rdedp22OeCiIiIDMXigoiIiAzF4oKIiIgMxeKCiIiIDMXigoiIiAzF4oKIiIgMxeKCiIhM\nc/iw5+dGOPPMM7F8+XKcd955iIuLw6233orjx49j5syZiI+Px5VXXona2loAwM6dO3HppZfCZrPh\nggsuwI4dOxzLWbNmDc4991zEx8cjPT0dzz/vvLzVjh07MHLkSPzud79DcnIyhg8fjjVr1hj/YUIE\ni4sAUP1dtdUhUADiduHEXASntjZg7lzgscfsz59+GpgxA2hpMX5dmzdvxgcffICDBw+ioKAAM2fO\nxLJly1BdXY3W1lasXLkSR48exezZs/Hwww+jpqYGy5cvR1ZWFk6cOAEASE5Oxttvv426ujqsXr0a\nd999N/bu3etYx7Fjx3Dq1CkcPXoUL774IhYvXuwoWqgrFhcWK6stQ+qKVJTXlVsdCgUQbhdOzEXw\niogA3nwTWL8eSEoCfv97YMsWYIAJY0PfcccdGDp0KFJTUzFlyhRMmjQJGRkZiI6OxjXXXIM9e/Zg\n3bp1mDVrFq666ioAwOWXX46LL74Yb7/9NgBgxowZGD16NABgypQpuPLKK/Hhhx861hEdHY2HHnoI\nkZGRmDFjBmJjY3HgwAHjP0wIYHFhkabWJmRtysL09dPR0taC6eumI2tTFppbm60OjSzE7cKJuQgN\nqanAnDlAdTXw058Co0aZs57kZOfI5DExMT2e19fXo7S0FJs2bUJiYiISExNhs9nw97//HRUVFQCA\nLVu2YPLkyRgyZAhsNhu2bNmC6mpnq9mQIUMQEeE8bA4aNAj19fXmfKAgx+LCItGR0Ui3paO4qhgA\nUFRVhHRbOqIioyyOjKzE7cKJuQgNTz8NvPoq8MknwN/+5jxF4m8iglGjRuGGG27AyZMncfLkSdTU\n1ODUqVO477770NTUhHnz5uG+++5DVVUVampqMGPGDHDIc9+wuLDQPZPvcTwWSJfnFL64XTgxF8Gt\nrQ346itg+3Zg4kRg2zagosKcPhd9kZOTg4KCAmzduhVtbW1obGzEjh07cPToUTQ1NaGpqQlDhw5F\nREQEtmzZgq1bt1oTaAhgcWGhmsYa5GTkoGRxCRZmLERNY43VIVEA4HbhxFwEt4gIYOVK56mQ1FTg\nD38wvs+FiHh83mH48OEoKCjA448/jqSkJKSlpWH58uVoa2tDbGwsVq5cifnz5yMxMRH5+fmYM2eO\nV+slJ14VlYgojPGqqOQrXhWViIiI/IbFBRERERmKxQUREREZisUFERERGYrFBRERERmKxQUREREZ\nisUFERERGYrFBRERERmKxQUREZELjzzyCK6//noAQFlZGeLj43mtkT4y4cK3REREoaFjiO+RI0ei\nrq7O4miCh+UtFyIyXUT2i8hBEbnfxesLROTz9ttHIvIjK+IkIiKivrG0uBCRCADPArgKwHgA2SLy\ng26zfQXgJ6p6HoDHALzg3yiJiKi/Cg4UmHZK4cwzz8Ty5ctx3nnnIS4uDrfeeiuOHz+OmTNnIj4+\nHldeeSVqa2sBADt37sSll14Km82GCy64ADt27HAs59ChQ8jMzERCQgKuuuoqVFdXO14rLS1FREQE\n2traAABr1qzBueeei/j4eKSnp+P55593zLtjxw6MHDkSv/vd75CcnIzhw4djzZo1pnz2QGV1y8VE\nAF+qaqmqNgPIB9DlMnSqulNVa9uf7gQw3M8xWqr6u+reZwoTzIUTc+HEXAS+4qpizMmfg91Hd5u2\njs2bN+ODDz7AwYMHUVBQgJkzZ2LZsmWorq5Ga2srVq5ciaNHj2L27Nl4+OGHUVNTg+XLlyMrKwsn\nTpwAACxYsAATJkxAdXU1HnzwQaxdu7bLOjpfBTU5ORlvv/026urqsHr1atx9993Yu3ev4/Vjx47h\n1KlTOHr0KF588UUsXrzYUeCEA6uLi+EAyjo9L4fn4uEWAFtMjSiAlNWWIXVFKsrryq0OxXLMhRNz\n4cRcBLaWthaMemoUznvuPERKJC556RKkLE9BY0uj4eu64447MHToUKSmpmLKlCmYNGkSMjIyEB0d\njWuuuQZ79uzBunXrMGvWLFx11VUAgMsvvxwXX3wx3n77bZSVlWH37t149NFHERUVhSlTpuDqq692\nu74ZM2Zg9OjRAIApU6bgyiuvxIcffuh4PTo6Gg899BAiIyMxY8YMxMbG4sCBA4Z/7kBldXHRZyIy\nDcBNAHr0ywg1Ta1NyNqUhenrp6OlrQXT101H1qYsNLc2Wx2a3zEXTsyFE3MRHAZEDMDqOauhqmjV\nVrS2teKlOS9h4ICBhq8rOTnZ8TgmJqbH8/r6epSWlmLTpk1ITExEYmIibDYb/v73v6OiogJHjx6F\nzWZDTEyM431paWlu17dlyxZMnjwZQ4YMgc1mw5YtW7qcRhkyZAgiIpyH2EGDBqG+vt6ojxvwrP61\nyBEAozo9H9E+rQsRyQDwPIDpqlrjaYF5eXmOx5mZmcjMzDQiTr+KjoxGui0dm0s2AwCKqoow6+xZ\niIqMsjgy/2MunJgLJ+bCd4WFhSgsLPTb+kYPHo1WbcWQmCE40XACI+JH+G3dnYkIRo0ahRtuuAF/\n/vOfe7x++PBh1NTUoKGhwVFgHD58uEuB0KGpqQnz5s3DunXrMGfOHEREROCaa67hz1Q7sbrlYheA\ndBFJE5FoANcBKOg8g4iMAvA6gOtV9d+9LTAvL89xC8bCosM9k+9xPBZIl+fhhrlwYi6cmAvfZGZm\ndtlPmi1hYAJWz1mNyv9dib9c8xckDUoyfZ3u5OTkoKCgAFu3bkVbWxsaGxuxY8cOHD16FKNGjcLF\nF1+M3NxcNDc346OPPsKbb77Z5f0dxUNTUxOampowdOhQREREYMuWLdi6dasVHylgWVpcqGorgNsB\nbAVQBCBfVUtE5D9F5Lb22R4CkAjgjyLymYh8alG4flXTWIOcjByULC7BwoyFqGn02GAT0pgLJ+bC\nibkIDkMHDcWi8xchMiISORk5SI1LNXwdnTtaunreYfjw4SgoKMDjjz+OpKQkpKWlYfny5Y5fgKxf\nvx47d+7EkCFD8Nvf/hY33nijy+XGxsZi5cqVmD9/PhITE5Gfn485c+b0WF9fYgpVEkrNOCKiofR5\niIjMJiJQ1V6PfNy/Uneeth2rT4sQERFRiGFxQURERIZicUFERESGYnFBREREhmJxQURERIZicUFE\nRESGYnFBREREhmJxQURERIZicUFERCGrsbERV199NWw2G37xi19gw4YNmD59uuP1iIgIfPXVV6bG\n8O6772Lu3LmmLNus+J999lk88MADPr+fxQUREYWs1157DVVVVTh58iReeeUVLFiwAO+8847j9c7D\nct900014+OGHDY/hwQcfxJIlSxzPIyIiEBcXh/j4eIwcORL33ntvny56Nm3aNLz00ktdppk1rPit\nt96K9evXd7nSqzdYXBARUcgqLS3F2LFj3R6EjRzSvLW1tce03bt3o66uDhMmTHBMExHs27cPdXV1\n+OCDD7Bhwwa88MILPq3TrCHZTzvtNMycORMvv/yyT+9ncRHEqr/zraIMRcyFE3PhxFwEAFXg7beB\n//kfoKDA/txg+/fvx7Rp02Cz2fCjH/3IcTXTvLw8PProo8jPz0d8fDxWr16NtWvXYsqUKT2W8cIL\nL2D9+vV48sknER8f77gQWUVFBebNm4dhw4ZhzJgx+P3vf+94zyOPPIL58+fj+uuvx+DBg7F27doe\ny92yZQumTp3aZZqqOoqCsWPHYsqUKfjiiy+wfPlyzJs3r8u8d911F+6++248+OCD+PDDD3H77bcj\nPj4ed955p2Oe9957D2PHjkViYiJuv/32Lut57LHHMHr0aKSkpGDRokWoq6sDYC+6IiIi8PLLLyMt\nLQ3Dhg3D448/3mXdU6dOxV//+tfe/wCudHzIULjZP054OPzNYR3w6AAtqy2zOhTLMRdOzIUTc9E3\n7ftN8/avd96pevrpqlFR9vtf/tKIsB2am5s1PT1dly1bps3Nzbpt2zaNi4vTgwcPqqpqXl6eXn/9\n9Y7516xZo1OmTHE8FxH997//raqqixYt0oceesjxWltbm1500UX62GOPaUtLi3799dc6ZswY3bp1\nq2PZ0dHRWlBQoKqqjY2NPeKbP3++Ll++vMu0zussKirSlJQUXb16tVZUVGhsbKzW1taqqmpLS4sO\nGzZMP/vsM1VVzczM1FWrVvVY1tVXX611dXV6+PBhTUpK0nfffVdVVVetWqVnn322Hjp0SL/99lud\nO3euIxeHDh1SEdHbbrtNv//+e/3888/1tNNO0/379zuWvWfPHh0yZIjb3HvadthyEWSaWpuQtSkL\n09dPR0tbC6avm46sTVlobm22OjS/Yy6cmAsn5iKAHDkCPP888O23QHOz/X7jRuDgQcNWsXPnTnz7\n7be4//77MWDAAEybNg2zZ8/Gxo0b+73sXbt2obq6Gr/5zW8QGRmJ0aNH45ZbbkF+fr5jnsmTJ+Pq\nq68GYD+V0N0333yDuLi4HtMvvPBCDBkyBHPmzMFtt92GRYsWISUlBT/5yU/w6quvArC3eiQlJeH8\n88/3GOeSJUsQFxeHkSNHYtq0adi7dy8AYMOGDbjnnnuQlpaGQYMGYenSpcjPz3dcYl5EkJeXh+jo\naGRkZOC8887D559/7lhuXFwcamtrvcyaHYuLIBMdGY10WzqKq4oBAEVVRUi3pSMqMsriyPyPuXBi\nLpyYiwBy8iQQ1S3vUVHAiROGreLo0aMYOXJkl2lpaWk4cuRIv5ddWlqKI0eOIDExEYmJibDZbFi6\ndCmOHz/umKf7uruz2Ww4depUj+mfffYZTpw4gS+//BKPPPKIY/oNN9yAdevWAQDWr1+P66+/vtc4\nk5OTHY8HDRqE+vp6APbcpKWlOV5LS0tDS0sLKisre30vAJw6dQoJCQm9rt8VFhdB6J7J9zgeC6TL\n83DDXDgxF07MRYA4+2wgJgbo3JlywADghz80bBVnnHEGysrKukw7fPgwhg8f7vWyunf6HDlyJM46\n6yycPHkSJ0+eRE1NDWprax19Oly9p7uMjAwcdNFSo276nvz85z/Hvn37UFRUhLfeegsLFy7s87q6\nO+OMM1BaWup4XlpaiqioqC4FhSclJSU477zzvFpnBxYXQaimsQY5GTkoWVyChRkLUdNYY3VIlmEu\nnJgLJ3/movsPBFz8YCB8DRwIFBYC55xjb7E4+2xg2zbAxWkCX02aNAmDBg3Ck08+iZaWFhQWFuKt\nt95Cdna218tKTk7uMmbExIkTERcXhyeffBKNjY1obW1FUVERdu/e3edlzpw5E4WFhX2e/7TTTkNW\nVhYWLFiASZMmYcSIEW7j6012djaeeuopHDp0CPX19fjNb36D6667DhER9kO/uwKnw44dOzBjxow+\nr68Ld50xgvGGMOrQSUTWa2lR/fGPVbdvtz/fvt3+vKXFyqi8A7M7dPpBcXGxTp06VRMSEnT8+PH6\nxhtvOF7rrUNnRESEo3Pll19+qeeff77abDa95pprVFW1oqJCs7OzNSUlRRMTE3Xy5Mn6wQcfuFy2\nOxMnTtRPP/3U5Tpd+eijj1REdO3atV2mf/zxxzp27FhNTEzUu+66y+WybrrpJken1La2Nv3tb3+r\nI0eO1GHDhukNN9yg33zzjaraO3RGRERoa2ur473Tpk1zdBhtaGjQESNG6PHjx93G6WnbEe2lcgkm\nIqKh9HmIKPAVFgLz5wOLFwN/+APw6qtAZqbVUfWdiEBVe21v5/7Vd++99x7+9Kc/YfPmzX2av6ys\nDOeccw6OHTuG2NhYk6Nz7dlnn0V5eTmWLVvmdh5P2w6LCyKifsrLAx55BMjNtT8OJiwuAktbWxvu\nuece1NfX48UXX7Q6HI88bTsD/B0MEQWH1lYgMtL980Dmz9gLC+0tFrm59vvMzOBquaDA8d133yE5\nORlnnnkmtmzZYnU4/cIOnUTUQ2srcNll9gMnYL+/7LLg6Kzoz9hbW4ElS+ynQvLy7PdLlgRHnijw\nDBo0CKdOncK+fft8+rVLIOFpESITBfO3/2DuS+DP2IP5bwzwtAj5ztO2w5YLIpME87d/wH4wXrzY\n3pdg8eLgKSwA/8bevZAIpsKCyCwsLohMEhkJLF1q/wadl2e/X7o0eA4+3fsSePFTfcP4OoZEIMRO\nFM5YXBCZyJ/foI0czCkQ+hL42vLjKXYOeEXkHywuyK/8tXMPlIOIv75BezoQ+5KLyEjgo4+cxVBm\npv25Wa0urmL0teXHXexAcJ+mIgoq7kbXCsYbAngEOfLfaIaBMmqiv+PYvl116FDV3Fz7/fbtgZML\nT3qLMTdXFbDf95erHIU7hMAInZ40NDTo7NmzdfDgwXrttdfq+vXr9aqrrnK83vny52Z55513HCN+\n+tO+ffv0kksuMW35nrYdywsCI2/BuvEHsu4Hob4clDy9x18790A5iPiSv/5wdSAOlFx44i5GM2I3\nslgJBaFeXPzlL3/RSZMmaVtbm8vXOw+fvWjRIsfQ2Ua6+OKLuwz/3bmgycvL06ioKI2Li1ObzaaX\nXnqpfvzxxz2W8fXXX2tERIT+6le/6vGaiGhsbKzGxcVpbGys2mw2x2uzZs3St956y/DPpOp527H8\ntIiITBeR/SJyUETudzPPShH5UkT2iojnC9uTYXw5593be/zVByFQfungyy8JjO7E6CkXgXL6yFWM\nvfX78CV2dvQMP6WlpRg7dqzbK4raj5HGaHWxEe7evRt1dXWYMGGCY1r3WK677jrU1dWhuroamZmZ\nmD9/fo/lvPzyy0hMTMQrr7yC5ubmLq+JCPbt24e6ujqcOnUKJ0+edLy2YMECPPfcc/39aN5zV3X4\n4wZ7n4//ByANQBSAvQB+0G2eGQD+2v54EoCdHpZneGVmJjO/1VZ9W2XIctx9c/S1dcLob6Lu4uhY\nz69zq3qsx9+tCd7w9TSGp/e5y0UgnTLxdjvzJfbO76n6tsr0zxvI21lnMLnloq2tTf/617/qk08+\nqW+88YbbFoT+KCkp0czMTB08eLD+8Ic/1IKCAlVVzc3N1ejoaEfLwEsvvaRr1qzRyy67zPHejlaE\n559/XqOiovS0007TuLg4/dnPfqaqqkePHtWsrCxNSkrSs846S1euXOl4b15ens6bN09zcnI0ISHB\ncdGvzh599FG99dZbu0zr3nLR+eJnxcXFGhERodXV1V3eM2bMGH3uuec0JSVFX3/9dbfL6+7IkSMa\nExOjTU1NvebRW562HauLix8D2NLp+QMA7u82z3MAftHpeQmAZDfLMzp3pjFzx374m8M64NEBWlZb\n1mV93dffV92bkfsSu6umZ6M/s7vlff+9/f6VLfY8bHqnzLGeQDqguuNrAebqb9zxeV3loj/r8oWR\nhULHfN7G3tLS9f/DzMIi0LezDmYXF3feeaeefvrpGhUVpaeffrr+8pe/NCTuDs3NzZqenq7Lli3T\n5uZm3bZtm8bFxenBgwdVtferonY+MHc/LdLW1qYXXXSRPvbYY9rS0qJff/21jhkzRrdu3epYdnR0\ntKOYaWxs7BHf/Pnzdfny5V2muSsuvv/+e73//vs1KSmpy9VK//a3v+nAgQP1m2++0TvuuMNR+Lha\nnivx8fH6r3/9y+3rvgrk4iILwPOdnucAWNltnjcBXNLp+fsALnSzPKNzZyqjd+zft3yvc1+Zq+f+\n4VxFHnT8H8br3FfmasP3TT7v6Hw5F+7pNU8HGFfTfYnv+5bv9Zr8rnm4Jn+uNrU09RqfP3n6zEb1\nC+gtF76uy9u/V28HW1///t7E7u7/o3MujBQo21lvzCwuysvLdeDAgQrAcYuJidEDBw4YFv+HH36o\nqampXaZlZ2frI488oqr9Ky4++eQTTUtL67LspUuXOgqkvLw8nTp1qsf4rrjiCv3zn//cZVr34iI6\nOlptNptGRkbq0KFDdceOHV3mv+WWW3Tu3Lmqar/senR0tFZVVXVZXkJCgg4ePFhtNpvjcuwdhg8f\nrh9++KHHOH0RVsWFq1uumz1Pbm6u5fN33jkatXxMhSLPfrtv632qat+xxcR4t/yHHnI9/0MP5faI\n3VM8w4fnujxYuJv/xhtzXRY/7uafOjW3xwHmvq332fPg5vO6OigZlf+O/PR1/tzcXJenMeyfy3M8\nXQsS1/NfknOJY3vovE0Y9Xl9+Xu5Oth6G8+NNxr7/9g9fnfbv9f/j8h1WfxYtf+58cYb2/8Hch3z\nqEnFxb59+zQuLq7L+uPj4/Uf//iH18ty55VXXtGJEyd2mfbAAw/obbfdpqr9Ky42bdqkAwYMUJvN\npjabTQcPHqzx8fE6e/Zsx7JzcnI8xnfttdf2ueXixIkTmpmZqStWrHDM29DQoAkJCfraa685po0Z\nM0afeeaZLsv76quv3MYQji0XPwbwTqfnfTktsh8hcFpE1djm7w7HTh1zHEQkT/TYqWOO14z8hupL\n64QnvubC3fs85cHIb5T9af729qejvq7LXS6Mjr0vjGqRCaRcuMOWC/uBcdiwYSoijuIiMTFR6+rq\nDIvfVcvFggULfGq5uOmmm7oUFx9//LGOHTvW7bq7L9uVxx57zFHouFpn92UcOHBA4+Li9Ngx+7a5\nfv16FREdNmyYpqSkaEpKisbExOiFF17ocnndhWufi0g4O3RGw96h85xu88yEs0PnjxEiHTp9bSbu\n7X0lVSWaszmny33HfIFwQPXE2wOPpzjc5aE/sftSaPXGXd8Ud+v1ZV3uctHbunyJ3ROjD7a+xO4p\nF8HwP2IGM4sLVXsHxXPPPVejoqL07LPP1r179xoSd4empiYdM2aMPvHEE9rc3Kzbt2/X+Ph4n/pc\nPPDAA7pw4ULHa62trXrRRRfpE088oQ0NDdrS0qJffPGF7tq1y+WyXdmzZ0+PAsVTcaFq76dx9913\nq6rqlVdeqbfccotWVlY6bv/85z81IiJCv/jiix7L627Dhg06a9YsjzH6KmCLC3tsmA7gAIAvATzQ\nPu0/AdzWaZ5n24uQz+HmlIgGWXGh6nsB4e1O0IwdndE94X1tCTF6HA5P7zF6oCerv/33RyBsg2Yw\nMrf8tYj/FBcX69SpUzUhIUHHjx+vb7zxhuO13oqLzuNcfPnll3r++eerzWZzDHpVUVGh2dnZmpKS\noomJiTp58mT94IMPXC7bnYkTJ3YZ56LzOl0t45NPPtHY2FgtLS3VqKgoLSoq6rHMWbNm6a9//ese\ny3M135tvvtlrjL4I6OLCyFsgb/ze6m3n7cu3fE/PrWTGqQAz+NK51R1//kLCaP35uayn51YLhNxa\nIRSKi0C3detWjtAZzLdQ2/jdFRBG93cIhJ2+0acCzIhD1bef5fq6LlfzB0qhZfSvfqwWSLn1NxYX\n5CsWF0HI3QHV6G/5wbJT9cepAF9PR/nzgBrIB+9g2ZbcCeTcmonFBfmKxUWQ6c+YAL58yw/05mB/\nxudLUUdOgb4tUU8sLshXLC6CkD978fv6Hn+w4qDuLhfh+s3WW4G6LZFrLC7IV562HbG/HhpEREPp\n8/iisBCYP99+8ac//MF+safeLtrly3v8qbW16wW/uj83UqDnItCFav78uQ36m4hAVV1f1avrfGG/\nf6WuPG477qqOYLwhzCvrUO5z4Q/MRf+Eav5C9XN1AFsuyEeeth22XIQYX75hhfK3Mm8xF/0TqvkL\n1RYZoO8tFzExMccaGxuT/RETBYeBAwdWNjQ0pLh6jcUFEVEf5OUBjzwC5ObaH4eKvhYXRN6IsDoA\nIqJAV1hob7HIzbXfFxZaHRFRYGPLBRGRB62twGWXAUuX2k+FFBYCS5YAH30UGqd82HJBZmBxQUTU\ni1DtSwKwuCBz8LRImKn+rtrqEAIGc+HEXDi5ykX3QiJUCgsis7C4CCNltWVIXZGK8rpyq0OxHHPh\nxFw4MRdExuBpkTDQ1NqE7Nezsb96P4qrijE+aTzGDR2H/Kx8REVGWR2eXzEXTsyFUzjngqdFyAxs\nuQgD0ZHRSLelo7iqGABQVFWEdFt6yO80XWEunJgLJ+aCyFi9tlyIyOkAGlS1TUTGAvgBgC2q2uyP\nAL3Blgv3KusrkbLCPtaJQFBxbwWSY8NzPBzmwom5cArXXLDlgszQl5aLvwEYKCLDAWwFcD2ANWYG\nRcaraaxBTkYOShaXYGHGQtQ01lgdkmWYCyfmwom5IDJOX1ou9qjqhSJyB4AYVX1SRPaq6vn+CbHv\n2HJBROR/0zDmAAAXXElEQVQdtlyQGfrSciEiMhnAQgB/bZ/GH2IRERGRS30pLv4XgCUA/q+qFonI\nWQC2mxsWERERBSv+FJWIKIzxtAiZYYC7F0TkTQBuj9Sq+jNTIiIiIqKg5ra4ALDcb1EQERFRyOBp\nESKiMMbTImQGTy0XAAARORvAUgDnAhjYMV1VzzIxLiIiIgpSffm1yGoAfwLQAmAagJcBrDMzKCIi\nIgpefSkuYlT1A9hPoZSqah6AWeaGRURERMGq19MiAL4XkQgAX4rI7QCOAIg1NywiIiIKVn1pubgL\nwCAAdwK4CEAOgBv7u2IRsYnIVhE5ICLvikiCi3lGiMg2ESkSkX+JyJ39XS8RERGZqy/FRauq1qtq\nuarepKpZqrrTgHU/AOB9VR0HYBvso4B21wLgHlUdD2AygMUi8gMD1k0uVH9XbXUIAYO5cGIunJgL\nor7pS3GxQkRKROS3IvJDA9c9B8Da9sdrAfy8+wyqekxV97Y/rgdQAmC4gTFQu7LaMqSuSEV5XbnV\noViOuXBiLpyYC6K+69M4FyKSAuBaAL8AEA/gFVV9rF8rFjmpqonunruYfzSAQgA/bC80XM3DcS68\n1NTahOzXs7G/ej+Kq4oxPmk8xg0dh/ysfERFRlkdnl8xF07MhVOo54LjXJAZ+tKhE6p6DMBKEdkO\n4D4ADwPotbgQkfcAJHeeBPuQ4g+6Wo2H5cQCeA3AXe4Kiw55eXmOx5mZmcjMzOwtzLAWHRmNdFs6\nNpdsBgAUVRVh1tmzQmKn6S3mwom5cAq1XBQWFqKwsNDqMCjE9dpyISLnwN5ikQXgBIBXALyuqsf7\ntWKREgCZqlrZ3jKyXVXPcTHfAABvAdiiqs/0sky2XPigsr4SKStSAAACQcW9FUiOTe7lXaGJuXBi\nLpxCORdsuSAz9KXPxUsAagBcpaqZqvqn/hYW7QoALGp/fCOANzysv7i3woJ8V9NYg5yMHJQsLsHC\njIWoaayxOiTLMBdOzIUTc0HkHcuuLSIiiQA2ARgJoBTAtar6jYikAnhBVWeLyKUA/gbgX7CfNlEA\n/62q77hZJlsuiIi8wJYLMgMvXEZEFMZYXJAZ+nJahIjIOlVVwK5d9nsj3ufr8oioz7wqLto7XlKA\nq6qqwq5du1DFnScFu40bgbQ04Ior7PcbN/bvfb4uz1csZChMeXVaRET2qOqFJsbTLzwtAmzcuBE3\n33wzoqOj0dTUhFWrViE7O9vqsIi8V1VlLwAaGpzTYmKA0lIgKcn79/3zn8BFF3m/PF9t3AjcfDMQ\nHQ00NQGrVgEB+L/I0yJkBm9Pi3ADDGBVVVW4+eab0dDQgNraWjQ0NODmm28O/RYMfjsMTYcO2Q/M\nnUVF2ad3cPW3d/e+Tz/tfXm+6h5HVZW9sGhoAGpr7fc338xtlMKGt8XFC6ZEQYY4dOgQorvtPKOi\nonDIiJ1noPJ3Mzf5z+jR9m/8nTU326cD7v/27t43caLn5fnKVRx9KYyIQphXxYWq/tGsQKj/Ro8e\njaZuO8/m5maM7u/OM1Dx22FPvrTieHqPla1CSUn2UwkxMUB8vP1+1Sr7dE9/e3fvO+cc98vrC1e5\ncBdHbKw5hQxRsFDVkLnZP05427Bhg8bExGh8fLzGxMTohg0brA7JPJ9+qpqQoAo4b/Hx9umqqseP\n2x8fP97zvZ5eM5q/1rVhg2pMjD0nMTH25/15jy/LM4Or/PX2t3f3Pk/TPXGXC09xdLwnPt7a/PWi\nfb9p+f6bt9C6cZyLEFRVVYVDhw5h9OjRSDKjo1qg8NTh7/333Xem82dHO3+ty5fOj57eA/jWmdJf\nfO3safS6AM9xVFXZT4WMHh0YeXOBHTrJDG5Pi4jIKA+vTTEnHDJCUlISJkyYENqFBeC++Rtw32Tu\nz1Mp/lyXL+f4Pb0n0PsMeDplYjRPuegtjqQkYMKEgC0siMzi6aqohSLyHIAVqtoKACKSDGAFgB8A\nuNgP8RF5lp0N/PSnXb8d7tplPxh0/jbZ+cDo7jWjDwAdByV/rKu3zo++vCfQ+wy4+tubobc8+SsO\noiDiqUPnRQDGANgrIv8hIncB+BTAxwAm+iM4oj7p/u3Q08HAl4Owr/y5Ll++yXt6jz9bBvrDHy0D\nfckFWyiIuujLJdfvAvAUgKMAfqyq5f4IzBfsc2Ge6u+qMXTQUKvD6LuOvg5RUfYDuqs+F65eA3o9\nT+5VLnpbl9Hcxe7pM/n6GoJwu+iPEM0F+1yQGdwWFyIyGMATACYBuA/ATACXA7hLVbf5LUIvsLgw\nR1ltGc5aeRa+vutrjIgfYXU4fefLQbOXDpg+5cLqTn0mdSo1fLuwOk/9ELT/I2BxQebwVFx8BeCP\nAJ5W1Zb2aee3TytV1YAbx5bFRT+42LE3tTYh+/Vs7K/ej+KqYoxPGo9xQ8chPysfUZFRloZrCg+/\nCmhKTAjOXJjwqwpTtosgGSq7u1D4H2FxQWbw1OfiJ6q6vKOwAABV3auqlwAIyJaLoOTPQYrcrcvN\nSIfRkdFIt6WjuKoYAFBUVYR0W3rQ7DS95uFXAR25OH6oGBcfASoPdctFoA5BbsKvPgzfLoJ4MLSw\n+x8h6iO3xYWnvhWqymHAjdCfoau9PZi5W1cvO/Z7Jt/jWIRAujwPOb10wFxSPhqlTwHvvQyUPgUs\nKT/TPk8gD0FuUqdSn7cLb64FEig/e+1FWP2PEPWV1aN4GXlDMI3Qefy4fdS+ziP7xcT0bdRAb0dO\n9LSuXkY6LKkq0ZzNOV3uLeHvUS67j6p4/Li2xgzskqfWmIGqxcW+/x39xYSRIkuqSvS/Vs/Tr7du\n0v9aPa9v24W77bY//wsBIGD+R3wEjtDJmwk3ywMw9MMEU3HRl+GLXfFlR+xpXcGwY/f3MNTeDDe9\nZo1vf0dfYujLa74szxdGFridlxfgQ2WHIhYXvJlxszwAQz9MMBUXvh7UfSlK+rtj9+d1OLyN3eo4\nzGi5CPTrfRhd4HZerlXbWRhjccGbGTfLAzD0wwRTcaHq27c1Xw+2vhYQVh/MfG3hMYO7HBr5rdvT\n3zdQCi0zClyyDIsL3sy48cJlVvPlt/2+DgLl7br8eXGoQI6hezxG5NadXbvsHUNra53T4uPtF2ID\n3L82YYLv6/SWr38Tfw8oRn3Cn6KSGVhcBCsfB4HyiqcDnT8PZuF0UAqWK5X6+jcJ9IGyAj0+E7C4\nIDOwuAhkVrc0BFKrQTjt9PszdLk/hdrfJEgH8uovFhdkBhYXgcqXHZ0ZLQ2BdDALJ/243gf5IJAK\naT9jcUFmYHERiHzd0Zm1g+TBjEJdoJwCtACLCzKDp+G/ySq+jlho1mWyeTlp8kagDoXuiUkjmRKF\nKxYXgag/O7rsbHtLxfvv2+95CoP8KZCHQvfErMKcKEzxtEigYl8HCjah0G8hDE8B8rQImcGy4kJE\nbABeAZAG4BCAa1W11s28EQB2AyhX1Z95WGboFBdAWO7oKIiFcb+FYMbigsxg5WmRBwC8r6rjYL+E\n+xIP894FoNgvUQWSIOjrUP1dtdUhBAyfchGM/RPcYb8Fl/g/QuHIyuJiDoC17Y/XAvi5q5lEZASA\nmQBe9FNc1EdltWVIXZGK8rpyq0OxnE+5CNb+Ce6091toixmI2tOAtpiBYd9vgf8jFK6sPC1yUlUT\n3T3vNP1VAP8HQAKAe8PqtEiAamptQvbr2dhfvR/FVcUYnzQe44aOQ35WPqIio6wOz698zkUo9E/o\npiMXlV9/ge//fRADx4zDsDPHc7sI8P8RnhYhMwwwc+Ei8h6A5M6TACiAB13M3qMqEJFZACpVda+I\nZLa/36O8vDzH48zMTGRmZnoVM/UuOjIa6bZ0bC7ZDAAoqirCrLNnBdxO0x98zkXHz407FxcdPzcO\n0uKiSy6GA2g8gPtsc7hdILD+RwoLC1FYWGh1GBTirGy5KAGQqaqVIpICYLuqntNtnscB5ABoARAD\nIA7AZlW9wc0y2XLhJ5X1lfjhoykY/Q1QOhj418PHkByb3PsbQ1BlfSVSVqQAAASCinsres9FCLZc\nAD7mIkQFSy7YckFmsLLPRQGARe2PbwTwRvcZVPW/VXWUqp4F4DoA29wVFuRfLev/gvJnIrEzPxZl\nz0SiZcM6q0OyTE1jDXIyclCyuAQLMxaiprGm9zeF6LgKPuUiRDEXFM6sbLlIBLAJwEgApbD/FPUb\nEUkF8IKqzu42/1Swz0VgCNFv3Zbgz43JYmy5IDNwEC3yHsczIAoZLC7IDBz+m7zXl/EMQmn8BiIi\n8gqLC/Jeb/0FQm38BiIi8gpPi5DvXPUXYH8MoqDC0yJkBlPHuaAQl5TUs2AIwfEbiIjIOzwtQsbi\n9SWIiMIeiwsyVoiO30BERH3HPhdkDo7fQBQU2OeCzMDigogojLG4IDPwtAgREREZisUF+R8H2CIi\nCmksLsi/OMAWEVHIY58LMk31d9UYOmioc0IYD7DVIxdhjLkILOxzQWZgywWZoqy2DKkrUlFeV+6c\n2DHAVmcdA2yFMJe5CFPMBVF44AidZKim1iZkv56N/dX70dLWgunrpmPc0HHIz8pHVJgNsOUxF5FR\nVofnV8wFUXhhywUZKjoyGum2dBRXFQMAiqqKkG5Ltx9AwmyALY+5CDPMBVF4YZ8LMlxlfSVSVqQA\nAASCinsrkByb7JwhjAbY6jUXYYS5CEzsc0FmYMsFGa6msQY5GTkoWVyChRkLUdNY03WGpCRgwoSQ\nLyyAPuQijDAXROGDLRdERGGMLRdkBrZcEBERkaFYXBAREZGhWFwQERGRoVhcEBERkaFYXBAREZGh\nWFwQERGRoVhcEBERkaFYXFDgqKoCdu2y3xMRUdBicUGBYeNG++XYr7jCfr9xo9URERGRjzhCJ1mv\nqspeUDQ0OKfFxAClpWExRDiRlThCJ5nBspYLEbGJyFYROSAi74pIgpv5EkTkVREpEZEiEZnk71jJ\nZIcOAdHRXadFRdmnExFR0LHytMgDAN5X1XEAtgFY4ma+ZwC8rarnADgPQImf4iN/GT0aaGrqOq25\n2T6diIiCjpXFxRwAa9sfrwXw8+4ziEg8gCmquhoAVLVFVev8FyKZpfq7aueTpCRg1Sr7qZD4ePv9\nqlVhc0qkSy7CHHNBFBqsLC6GqWolAKjqMQDDXMxzJoBqEVktIntE5HkRifFrlGS4stoypK5IRXld\nuXNidra9j8X779vvs7OtC9CPXOYiTDEXRKHD1A6dIvIegOTOkwAogAcBrFHVxE7znlDVId3efxGA\nnQAmq+puEXkaQK2q5rpZHzt0BrCm1iZkv56N/dX7UVxVjPFJ4zFu6DjkZ+UjKjLK6vD8irlwYi6s\nxQ6dZIYBZi5cVa9w95qIVIpIsqpWikgKgOMuZisHUKaqu9ufvwbgfk/rzMvLczzOzMxEZmamt2GT\nSaIjo5FuS8fmks0AgKKqIsw6e1ZYHkCYCyfmwr8KCwtRWFhodRgU4iz7KaqIPAHgpKo+ISL3A7Cp\n6gMu5tsB4FZVPSgiuQAGqarLAoMtF4Gvsr4SKStSAAACQcW9FUiOTe7lXaGJuXBiLqzDlgsyg5V9\nLp4AcIWIHABwOYBlACAiqSLyVqf57gSwXkT2wv5rkcf9HikZpqaxBjkZOShZXIKFGQtR01hjdUiW\nYS6cmAui0MJBtIiIwhhbLsgMHP6biIiIDMXigoiIiAzF4oKIiIgMxeKCiIiIDMXigoiIiAzF4oKI\niIgMxeKCiIiIDMXigoiIiAzF4oKIiIgMxeKCiIiIDMXigoiIiAzF4oICSvV31VaHQERE/cTiggJG\nWW0ZUlekoryu3OpQiIioH3hVVLJcU2sTsl/Pxv7q/SiuKsb4pPEYN3Qc8rPyERUZZXV4RCGNV0Ul\nM7DlgiwXHRmNdFs6iquKAQBFVUVIt6WzsCAiClJsuaCAUFlfiZQVKQAAgaDi3gokxyZbHBVR6GPL\nBZmBLRcUEGoaa5CTkYOSxSVYmLEQNY01VodEREQ+YssFEVEYY8sFmYEtF0RERGQoFhdERERkKBYX\nREREZCgWF0RERGQoFhdERERkKBYXREREZCgWF0RERGQoFhdERERkKBYXREREZCgWF0RERGQoy4oL\nEbGJyFYROSAi74pIgpv57haRL0Rkn4isF5Fof8dKREREfWdly8UDAN5X1XEAtgFY0n0GETkDwB0A\nLlTVDAADAFzn1yiDVGFhodUhBATmwYm5cGIuiMxlZXExB8Da9sdrAfzczXyRAE4XkQEABgE46ofY\ngh53nnbMgxNz4cRcEJnLyuJimKpWAoCqHgMwrPsMqnoUwAoAhwEcAfCNqr7v1yiJiIjIKwPMXLiI\nvAcgufMkAArgQRez97hWuogMhr2FIw1ALYDXRGSBqm4wIVwiIiIygKj2OKb7Z8UiJQAyVbVSRFIA\nbFfVc7rNMw/AVap6a/vz6wFMUtXb3SzTmg9DRBTEVFWsjoFCi6ktF70oALAIwBMAbgTwhot5DgP4\nsYgMBPA9gMsB7HK3QP6DEBERWc/KlotEAJsAjARQCuBaVf1GRFIBvKCqs9vny4X9FyLNAD4DcIuq\nNlsSNBEREfXKsuKCiIiIQlNIjNApItNFZL+IHBSR+62OxyoiskpEKkVkn9WxWE1ERojINhEpEpF/\nicidVsdkFRE5TUQ+EZHP2nORa3VMVhORCBHZIyIFVsdiJRE5JCKft28bn1odD4WOoG+5EJEIAAdh\n749xFPY+Gdep6n5LA7OAiFwGoB7Ay+2DjoWt9k7CKaq6V0RiAfwTwJxw3C4AQEQGqep3IhIJ4O8A\n7lTVsD2YiMjdAC4CEK+qP7M6HquIyFcALlLVGqtjodASCi0XEwF8qaql7X0x8mH/+WrYUdWPAHAn\nAfvYKaq6t/1xPYASAMOtjco6qvpd+8PTYO/IHdzfKvpBREYAmAngRatjCQCC0DgOUIAJhY1qOICy\nTs/LEcYHEepJREYDOB/AJ9ZGYp320wCfATgG4D1VdfurqzDwFIBfI4wLrE4UwHsisktEbrU6GAod\noVBcELnVfkrkNQB3tbdghCVVbVPVCwCMADBJRM61OiYriMgsAJXtrVrSfgtnl6rqhbC35CxuP7VK\n1G+hUFwcATCq0/MR7dMozLVfj+Y1AH9RVVfjqIQdVa0DsB3AdKtjscilAH7W3tdgI4BpIvKyxTFZ\nRlUr2u+rAPxf2E8zE/VbKBQXuwCki0ha++XYr4N9gK5wxW9jTi8BKFbVZ6wOxEoiMlREEtofxwC4\nAkBYdmxV1f9W1VGqehbs+4ptqnqD1XFZQUQGtbfsQUROB3AlgC+sjYpCRdAXF6raCuB2AFsBFAHI\nV9USa6OyhohsAPAPAGNF5LCI3GR1TFYRkUsBLATwH+0/s9sjIuH6bT0VwHYR2Qt7v5N3VfVti2Mi\n6yUD+Ki9L85OAG+q6laLY6IQEfQ/RSUiIqLAEvQtF0RERBRYWFwQERGRoVhcEBERkaFYXBAREZGh\nWFwQERGRoVhcEBERkaFYXBD5kYicsjoGIiKzsbgg8i8OLENEIW+A1QEQBTMRWQqgTFX/2P48F0AL\ngGkABgOIAvCQqhZ0e18KgFcAxMH+f/hfqvp3f8ZORGQWjtBJ1A8icj6Ap1U1s/15EezXaKhV1XoR\nGQJgp6qe3f56narGi8g9AE5T1aUiIgAGqeq3Fn0MIiJDseWCqB9Uda+IJLW3RAwDcBLAMQDPiMgU\nAG0AzhCRYap6vNNbdwFYJSJRAN5Q1c/9HjwRkUnY54Ko/14FMB/AL2A/1ZEDYAiAC1T1AgDHAQzs\n/AZV/RDATwAcAbBGRHL8GjERkYlYXBD13ybYL9+dBXuhkQDguKq2icg0AGmd5hUAEJFR7fOsAvAi\ngAv9GzIRkXl4WoSon1S1WETiAJSraqWIrAfwpoh8DmA3gJLOs7ffZwL4tYg0AzgF4AZ/xkxEZCZ2\n6CQiIiJD8bQIERERGYrFBRERERmKxQUREREZisUFERERGYrFBRERERmKxQUREREZisUFERERGYrF\nBRERERnq/wO4kjupiGFBcQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3481e51d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(sky_vals, sky_mean - sky_vals, color='b', marker='x', label='mean')\n",
    "plt.scatter(sky_vals, sky_med - sky_vals, color='g', marker='*', label='median')\n",
    "plt.scatter(sky_vals, sky_ofil - sky_vals, color='r', marker='o', label='ofilter (Python)')\n",
    "plt.scatter([1], [0.7501699 - 1], color='k', marker='o', label='ofilter (IRAF)')\n",
    "plt.xlabel('vals')\n",
    "plt.ylabel('X - vals')\n",
    "plt.axhline(0, color='k', linestyle='--')\n",
    "plt.legend(bbox_to_anchor=(1.01, 1), loc='upper left', scatterpoints=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Try with Other Skewed Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Number of data points for plotting: 45\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import gumbel_r\n",
    "\n",
    "# Populate desired background values\n",
    "vals = np.arange(0.1, 4.6, 0.1)\n",
    "\n",
    "# Initialize arrays to store results\n",
    "sky_vals = []\n",
    "sky_ofil = []\n",
    "sky_med = []\n",
    "sky_mean = []\n",
    "\n",
    "# Generate results\n",
    "for i, val in enumerate(vals):\n",
    "    np.random.seed(i)  # Does this control Scipy?\n",
    "    data = gumbel_r.rvs(loc=val, size=(50, 50))\n",
    "    \n",
    "    try:\n",
    "        msky = ofiltsky.fitsky_ofilter(data, sigclip_sigma=None)[0]\n",
    "    except ValueError as e:\n",
    "        print('i={0}, val={1:.1f}, errmsg={2}'.format(i, val, str(e)))\n",
    "        continue\n",
    "        \n",
    "    sky_vals.append(val)\n",
    "    sky_ofil.append(msky)\n",
    "    sky_med.append(np.median(data))\n",
    "    sky_mean.append(data.mean())\n",
    "    \n",
    "# Convert result to Numpy arrays\n",
    "sky_vals = np.asarray(sky_vals)\n",
    "sky_ofil = np.asarray(sky_ofil)\n",
    "sky_med = np.asarray(sky_med)\n",
    "sky_mean = np.asarray(sky_mean)\n",
    "    \n",
    "print()\n",
    "print('Number of data points for plotting:', sky_ofil.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f3481dad9b0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg4AAAEPCAYAAAAqD0wBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8FPW5P/DPEyBAmgQChKBcAgq2tqd4F1DR9YJcvIBG\nWlHAnnp7VaR6OBX08LNErVXU3tC2akWsVqClx3NEC17Lqm1F8QhoEcRWxXARQhO5CYYkz++P2V03\ny+5mdne+O5f9vF+vfSWzmZ19Znay88zz/c53RFVBREREZEeR2wEQERGRfzBxICIiItuYOBAREZFt\nTByIiIjINiYOREREZBsTByIiIrLNeOIgImNEZIOIbBSRWUn+/gMRWS0ib4vIuyLSLCLdTcdFRERE\nmROT4ziISBGAjQDOBrAVwCoAl6rqhhTznw/gRlU9x1hQRERElDXTFYeTAXygqptU9SCAxQDGp5l/\nEoBFhmMiIiKiLJlOHPoCqIub3hx57hAi0hXAGAD/bTgmIiIiypKXOkdeAOAvqvqZ24EQERFRch0N\nL38LgAFx0/0izyVzKdI0U4gIb6pBRJQFVRW3Y6DgMF1xWAVgsIhUi0gxrORgaeJMItINwBkAnk63\nMFUN7GPOnDmux8D14/oV2roVwvoROc1oxUFVW0TkegAvwEpS5qvqehG51vqzPhyZdQKA51V1v8l4\niIiIKDemmyqgqs8B+GrCcw8lTP8WwG9Nx0JERES58VLnyIIWCoXcDsEorp9/BXndgOCvH5HTjA4A\n5SQRUb/ESkTkFSICZedIcpDxpgoiIvK2rl27fnrgwIEqt+Mgb+nSpcv2/fv390l8nhUHIqIAs1Nx\n4PcrJZNq32EfByIiIrKNiQMRERHZxsSBiIiIbGPiQERERLYxcSAiIiLbmDgQERGRbUwciIgoK6pA\nXd2X03v2AA0Nzr/PoEGDcN999+GYY45BWVkZrr76auzYsQPjxo1DeXk5zj33XOzatQsAsHLlSpx6\n6qmoqKjAcccdh1deeSW2nMceewxf//rXUV5ejsGDB+Phhx+O/e2VV15B//798dOf/hRVVVXo27cv\nHnvsMedXJgCYOBARUVZWrwZOOglYu9ZKGsaNAx56qP3XZeOpp57Cyy+/jI0bN2Lp0qUYN24c7r77\nbuzcuRMtLS2YN28etm7divPPPx8//OEP0djYiPvuuw81NTX417/+BQCoqqrCsmXLsHv3bixYsAD/\n8R//gTVr1sTe49NPP8WePXuwdetWPPLII5g2bVosIaEvMXEgIqKsHH888MADwKmnAn37AkcfDcya\nZea9pk+fjl69euGwww7DyJEjMWzYMAwdOhTFxcW46KKL8Pbbb+N3v/sdzjvvPIwePRoAcPbZZ+PE\nE0/EsmXLAABjx47FwIEDAQAjR47Eueeei9deey32HsXFxbj11lvRoUMHjB07FqWlpXj//ffNrJCP\nMXEgInJYS0v66SAZPRrYt8+qOHzve0CRoaNKVdWXI2J37dr1kOm9e/di06ZN+MMf/oAePXqgR48e\nqKiowF//+lds27YNALB8+XKMGDECPXv2REVFBZYvX46dO3fGltOzZ08Uxa1ASUkJ9u7da2aFfIyJ\nAxGRg1pagNNOA8JhazoctqaDmDxEmyeuvhr4/e+BsWOtZgs3iAgGDBiAqVOnoqGhAQ0NDWhsbMSe\nPXswc+ZMNDU14ZJLLsHMmTNRX1+PxsZGjB07FhxqO3O8yRURkYM6dADuuguYOBGYNg345S+BJUus\n54Nm2zbglFOs9S0qsh7vvAMcc4w78UyePBknnngiampqcM4556CpqQlvvPEGhgwZgvLycjQ1NaFX\nr14oKirC8uXL8cILL+Cb3/ymO8H6GCsOREQOC4WspOG226yfoZDbEZlx1FHA3LlfNk9ccgkwZYrz\n7yMiaaej+vbti6VLl+LHP/4xKisrUV1djfvuuw+tra0oLS3FvHnzMHHiRPTo0QOLFy/G+PHjM3pf\nsvDumEREDguHD604uJU88O6YlK1U+w4TByIiB0X7ONx1l5UshMPALbcAf/mLO80VTBwoW0wciMgR\nLS1tD4CJ0+StbcTEgbKVat9hHwcisq0Qrhhw4lLKxCSBiRUFCRMHIp/L55gB8VcM1NZaP++6KzgH\nxkJIjIhyxcShwBTSwDSFwI0DXZCvGAh6YkTkBOOJg4iMEZENIrJRRJIORioiIRFZLSJ/F5EVpmPy\nAxMHeK+fTTGpyZwbB7pw2LpSYM4c62d0fwqKICdGRI5QVWMPWInJPwBUA+gEYA2AryXM0w3AOgB9\nI9O9UixLC0Vzs+rw4aorVljTK1ZY083NuS97xQrVXr1U58yxfkbfw20m17kQzJmjClg/TSqEz8mr\n/yPZinx3tvdd7Vp85F2p9h3TicNwAMvjpm8GMCthnu8BuN3GskxtG08y+eVl9yCTeDAwfXBwYp3z\nHbMX5PtAF+RtHMTEiIkDZSvVvmO6qaIvgLi7tWNz5Ll4RwHoISIrRGSViBgYd8x/TJVL7ZaZ/dh2\n7vWmGBNaWqwxApYssZoqliyxplOtM68YSK9DB2u8hei+Fwo5N/4Cm+K857bbbsOUyFCXdXV1KC8v\n570rbPDCvSo6AjgewFkAvgLgdRF5XVX/kThjbW1t7PdQKIRQgBsfEw/woVDuyUP8QSa6vFQD07gx\n3n6u61wI9whINh5A/OcXf6BLnLepCTjjjOQDEwHeGXfAbSYSo3wOChUOhxEOWscTg6LDSvfv3x+7\nd+92ORqfSFaGcOoBq6niubjpZE0VswDMiZt+BEBNkmUZKsZ4j8lyaaZlZj+2necr5nzLZBulmvel\nlw5t1ghied6L3Oo7ATZVpFRbW6tTpkxxOwzPSrXvmE4cOuDLzpHFsDpHHp0wz9cAvBiZtwTAuwC+\nnmRZZreQx+SzHTnVe/mx7TxZzEFqk8/kM0k1b7LEyssdAjP5/Lz+WbuR1OYjcfh0z6f6et3ruQWa\nxsCBA/Xee+/VoUOHamlpqV511VW6fft2HTt2rJaVlemoUaP0s88+U1XV119/XU855RTt3r27Hnvs\nsRoOh2PL+eijj/SMM87Q8vJyPffcc/X666+PJQ4ff/yxioi2tLSoquqCBQv06KOP1rKyMj3yyCP1\noYceii0nHA5rv3799Cc/+Yn27t1bDz/8cF2wYIGx9XeLK4mD9b4YA+B9AB8AuDny3LUAromb5wew\nrqx4B8D0FMsxuX0KVqqzzS++8N9ZaLJ1GTbMevhpPdqTycEncd50CYIXKzVOVFlSfdZ+7PybjXwk\nDjcsv0Grf1ad0zLSGThwoI4YMULr6+t169at2rt3bz3hhBN07dq1+sUXX+hZZ52lt99+u27ZskV7\n9uypzz33nKqqvvTSS9qzZ0/duXOnqqqOGDFCf/CDH2hTU5O++uqrWlZW1iZxKCoqiiUOy5Yt048+\n+khVVV999VUtKSnR1atXq6qVOHTs2FFra2u1ublZly1bpiUlJbHkJShcSxycejBxMCfVF5rXz96S\nSRZzkNYvl4rDSy+lPrCaOqiZqiLlOm++m2fcbA4ymTj8aeOftOTOEu14e0ctvqNYu/yoi0770zRn\nAo8zcOBAXbhwYWy6pqZGr7vuutj0/fffrxMmTNC5c+fq1KlT27x29OjR+vjjj+snn3yinTp10s8/\n/zz2t8suuyxl4pBowoQJOm/ePFW1EoeSkpI28/bu3VvfeOON3FfWQ5g4BIDJA50XzzadlLh+fmzX\nd+Ls+4svDp3P1LZwq9+K3XnzffmvW4mqycShqblJv/v0d7XLj7ooaqFfvf+runX3VmcCjzNw4EB9\n+eWXY9OTJ0/W2267LTb9yCOP6DnnnKPXXXeddunSRSsqKrSiokK7d++upaWlOnfuXF25cqX27t27\nzXJvueWWtBWH4cOHa48ePbR79+7auXNn/eEPf6iqVuLQv3//tDEGQap9h0NO+4TJSw2DPhJgsvXz\n49DCmVwqmGre4uJD5zN1CaJT2ziT/TOTefN9+W8QL2Pt1KETepf0RlNLE0qLS6FQHFZ2mCuxiAgG\nDBiAqVOnoqGhAQ0NDWhsbMSePXswc+ZMHHbYYWhsbMT+/ftjr/nkk0+SLqupqQmXXHIJZs6cifr6\nejQ2NmLs2LHRJIuSZRNefIAVh4zOkOye3fjhzDuXM7X21i/olRYvSLaNc90/naicOFFx8HKn0igY\n7uOw4qMV+sbmN3TH3h366NuP5h5wEnYqDqNGjdLNmzdrnz599Pnnn9eWlhbdv3+/hsNh3bJli6pa\nfRxuuukmbWpq0tdee03Ly8uTdo7cs2ePduzYUV999VVV1VgfhltvvVVVWXFwPSGw+2DiYLFzoPN6\nJ7FMOJHYZHLViBPbItUyvLydTUm1jXPZP9N13HUjYfZ68mk6cciHQYMGtTkoT5kyJWnioKr65ptv\n6hlnnKE9evTQ3r176/nnn691dXWqqvrhhx/qyJEjtaysTM8991ydPn16yqaKX/3qV1pVVaUVFRU6\ndepUnTRpUtrEITHGIGDiEAAmOon5gYl1MXUFRtCuUkk33d786dY5188009cnW5d8d9x0SxASB3IH\nEwefy+YMyetnQplI1rkxnlOVAZPlaz8cZKKyqQrYbVKIynX/zOSeK17v/GkSEwfKFhOHAMjkYJnv\ng5TJ8nwmlxU6wYmEK9Uy/JTMZboP5bMiZjK2TPih+YmJA2WLiUMB8co16k6U51MtO9nQyU5gxaGt\nTBMdE31wnHq9n5I2JzFxoGwxcSgwXhkVz4kOiKnmd/pAkM0ByW67vh/7OJg8q891/8z09X5M2pzC\nxIGyxcQhwLxSLrVTnneqGuKF0nOm7fpe+ZzscKqPgxfW0cux5QMTB8oWE4eA8sqXYiYVh1wP+m6s\ncyaXdJp8v3xyqjJk6v3yuWwvfB7ZYuJA2WLiEGBul2Gz6eOQazNDPr/I8z2IlFeSwXzy8jp7OTY7\nmDhQtpg4BJzbHb8yuarC7UQnG/nu8OjHbZQrL6+zl2NrDxMHyhYThwDz05ean8/e8n2jLLeTQTd4\neZ29HFs6QU8c9u/fr+eff752795dv/Wtb+mTTz6po0ePjv1dRPSf//yn0Riee+45veiii4ws21T8\n999/v86aNSvtPEwcAspL7f35er0bUiVnptbFT8mgU7y8zl6OrT1BTxyeeOIJHTZsmLa2tib9e1FR\nUezA+53vfCc2bLSTTjzxRH3zzTdj0yKipaWlWlZWpv369dMZM2akjC9eKBTS+fPnt3kuPn4nHThw\nQPv166f19fUp50m17/DumD5n6s6GqThxl06/3SWwpQW45RZgyRLrLo9LlljTLS1m1iXd+wWVl9fZ\ny7ERsGnTJhx11FEQkaR/t45/zmhJ8qG/9dZb2L17N0466aTYcyKCd955B7t378bLL7+MhQsX4je/\n+U1W7+lk/PE6d+6McePG4fHHH8/8xcmyCS8+4OOMOGj8fPaVrXxXSfxYlcmVl9fZy7G1B6YrDp99\npvrgg6o//anq+vU5x5vM+vXrNRQKaffu3fXf/u3fdOnSpaqqOmfOHC0uLtZOnTppWVmZPvroo/rY\nY4/paaedFntttNT/8MMPa6dOnbRz585aVlamF154oaqqbt26VWtqarSyslKPOOIInTdvXuy1tbW1\neskll+jkyZO1W7duh1QDVFVvv/12vfrqq9s8l9i8MHHiRJ0+fbree++9WlNT02be73//+3rjjTfq\n7NmztUOHDtq1a1ctKyvT6dOnx5b14IMP6pAhQ7SiokKnTZsWe21ra6vecccdWl1drVVVVXrFFVfo\nrl27VPXLu33+9re/1QEDBmhlZaXeeeedbd77ySef1LPOOivldk+177ieENh9MHHwFr+29xIVGqOJ\nQ0ODav/+qiUlqsXF1s9XXnEk7qiDBw/q4MGD9e6779aDBw/qn//8Zy0rK9ONGzeqqnVwj97hUlX1\nscce05EjR8am4w/iiU0Vra2tesIJJ+iPfvQjbW5u1o8++kiPPPJIfeGFF2LLLi4ujiUqBw4cOCS+\niRMn6n333dfmufj3XLdunfbp00cXLFig27Zt09LS0tjBvbm5WXv37q2rV69W1eRNFSKiF1xwge7e\nvVs/+eQTrays1Oeff15VVefPn69DhgzRjz/+WPft26cXX3zxIbcJv+aaa/SLL77QtWvXaufOnXXD\nhg2xZb/99tvas2fPlNs+1b7DpgrKWDgM/PKXwJw51s9oswVlLrHyyfI3+coDDwDbtwOffw40NVk/\nr7vO0bdYuXIl9u3bh1mzZqFjx44488wzcf7552PRokU5L3vVqlXYuXMnZs+ejQ4dOmDgwIG46qqr\nsHjx4tg8I0aMwAUXXADAKu8n+uyzz1BWVnbI88cffzx69uyJ8ePH45prrsF3vvMd9OnTB6effjqW\nLFkCAFi+fDkqKytx7LHHpo3zlltuQVlZGfr3748zzzwTa9asAQAsXLgQM2bMQHV1NUpKSnDXXXdh\n8eLFaG1tBWA1mdTW1qK4uBhDhw7FMcccg7Vr18aWW1ZWhl27dmW41YCOGb+CClp8e28oZD1uucVs\nv4qgivYXuesuazuGw9yW5DM7dlgJQ7yGBkffYuvWrejfv3+b56qrq7Fly5acl71p0yZs2bIFPXr0\nAGBV4FtbW3H66afH5kl870QVFRXYs2fPIc+vXr0agwYNOuT5qVOn4sEHH8SVV16JJ598ElOmTGk3\nzqqqqtjvJSUl2Lt3LwBr21RXV8f+Vl1djebmZmzfvr3d1wLAnj170K1bt3bfPxErDpSRfHfGDLIO\nHaykYeJEq9PdxInWNLcl+cZ55wElJV9Od+kCjBvn6FscfvjhqKura/PcJ598gr59+2a8rMQOlP37\n98cRRxyBhoYGNDQ0oLGxEbt27cIzzzyT8jWJhg4dio0bNx7yvKbo1DhhwgS88847WLduHZ599llc\nfvnltt8r0eGHH45NmzbFpjdt2oROnTq1SRbSWb9+PY455piM3hNg4kBZ8NtVEV4WCgHTpgG33Wb9\njCZkRL4wZgxw771A9+5W0jB+PHD//Y6+xbBhw1BSUoJ77rkHzc3NCIfDePbZZzFp0qSMl1VVVYUP\nP/wwNn3yySejrKwM99xzDw4cOICWlhasW7cOb731lu1ljhs3DuEM2ms7d+6MmpoaXHbZZRg2bBj6\n9euXMr72TJo0CT/72c/w8ccfY+/evZg9ezYuvfRSFBVZh/ZUyUvUK6+8grFjx9p+vygmDkQuYn8R\n8r3rrgMaG4H9+4HFi4GuXR1dfKdOnfDMM89g2bJl6NWrF66//no88cQTGDJkiK3Xx5/FX3nllVi3\nbh169OiBiy++GEVFRXj22WexZs0aDBo0CL1798bVV1+N3bt3247vuOOOQ/fu3bFq1aqk75nMFVdc\ngXfffRdTp05t8/wNN9yAJUuWoGfPnrjxxhuTLit++rvf/S6mTJmC008/HUceeSRKSkowb968lHHE\nTx84cADLli3DFVdcYXNN45bTXkaSKxEZA+DnsJKU+ao6N+HvZwB4GkA0zXpKVX+UZDlqOlaifGIf\nB8oHEYGqpj2S8fs1Ny+++CJ+/etf46mnnrI1f11dHY4++mh8+umnKC0tNRxdcg888AA2b96Mu+++\nO+U8qfYdo4mDiBQB2AjgbABbAawCcKmqboib5wwA/6mqF7azLO7YFDiJg0glG1SKKBdMHLyltbUV\nM2bMwN69e/HII4+4HU5aqfYd01dVnAzgA1XdFAliMYDxADYkzJdZjxCigGB/EaLC8fnnn6OqqgqD\nBg3C8uXL3Q4na6YTh74A4rvDboaVTCQaISJrAGwBcJOqvmc4LiIiorwqKSlJeumm33hhHIf/AzBA\nVT8XkbEA/hfAUclmrK2tjf0eCoUQYhd0IqI2wuFwRr38iTJluo/DcAC1qjomMn0zrCEs56Z5zUcA\nTlDVhoTn2QZHRJQh9nGgbKXad0xfjrkKwGARqRaRYgCXAliaEFhV3O8nw0pmnB16jIiIiBxhtKlC\nVVtE5HoAL+DLyzHXi8i11p/1YQCXiMj3ABwEsB/At03GREREbXXp0mV7/EkcEWDtF8meNz6Og1NY\nSiMiypydpgqiTHDkSCIiIrKNiQMRERHZxsSBiIiIbGPiQERERLYxcSAiIiLbmDgQERGRbUwciIiI\nyDYmDkRERGQbEwciIiKyjYkDERER2cbEgYiIiGxj4kBERES2MXEgIiIi25g4EBEZsvPznW6HQOQ4\nJg5ERAbU7arDYT85DJt3b3Y7FCJHiaq6HYMtIqJ+iZWICldTSxMm/fckbNi5Ae/Vv4dvVH4DX+31\nVSyuWYxOHTrlPR4RgapK3t+YAosVByIiBxV3KMbgisF4r/49AMC6+nUYXDHYlaSByARWHIiIHLZ9\n73b0+UkfAIBAsO0/t6GqtMqVWFhxIKex4kBEWWHHv9QaDzRi8tDJWD9tPS4fejkaDzS6HRKRY1hx\nIKKM1e2qwxHzjsBHN3yEfuX93A6H0mDFgZzGigMR2dbU0oSaP9RgzJNj0NzajDG/G4OaP9TgYMtB\nt0NzHCsqRMkxcSAi2wql4x8vpSRKjU0VBWrn5zvRq6SX22GQD3mp45/TvHYppRPYVEFOM15xEJEx\nIrJBRDaKyKw0850kIgdF5GLTMRU6nk1RLoLc8a9QKipEuTBacRCRIgAbAZwNYCuAVQAuVdUNSeZ7\nEcB+AI+q6lNJlsWKQ46CeDZF5LSgVVRYcSCnma44nAzgA1XdpKoHASwGMD7JfNMB/BHADsPxFDSe\nTdnnx45xfozZi4JcUSFyQruJg4h8JVIRgIgcJSIXiojdI01fAHVx05sjz8Uv/3AAE1T11wCYFduQ\nywFixogZsd8F0maaLH5syvFjzF71tV5fwxMXPdHmJxF9qaONeV4FMFJEKgC8AKu54dsALncohp8D\niO/7kDJ5qK2tjf0eCoUQCoUcCsE/cr1+Pno2NXvkbNz52p1oPNAYK8Om6jBZKB0p45tyopcaer0p\nx48xk1nhcBjhcNjtMCjA2u3jICJvq+rxIjIdQFdVvUdE1qjqse0uXGQ4gFpVHROZvhmAqurcuHk+\njP4KoBeAfQCuUdWlCcsq6D4OpvsnpEpIMk1U/J5kzHpxFu752z2x6ZmnzMTcUXPTvMJ92cTs98+J\n7GMfB3KanT4OIiIjYFUY/hR5roPN5a8CMFhEqkWkGMClANokBKp6ROQxCFY/h+sSkwYy1z8h1YA+\n+5r2ZTzQTxDK5X5sysk05iB8TkTkHjuJw40AbgHwP6q6TkSOALDCzsJVtQXA9bCaONYBWKyq60Xk\nWhG5JtlLbMZdkEwc1FIlJF8p/ortRCVIown6sWOc3ZiD9Dl5DTumUiHhAFA+smHnBtz52p2x/gmz\nR852pONWqsvPMrkszckSP8vo5nitWSPXZXthX/H6fTvYVEFOS1lxEJFnRGRpqkc+gySLqd7eqc5Y\nMzn7dqoawjK6WV5q1sh12dm83snKACs4VKjSXVVxX96iIFdFExEAsZ/pnk8m3dUadvDqgPyw+zmZ\n/DxyXXa2r3e6MhBt5ntqvTVe3br6dThvyHncXynwUlYcVPWVdI98BlnI/NJ2mms1hINT5Yfdz8nk\n55HrsjN9vcnKgB870xLlys4AUENE5I8i8p6IfBh95CO4QldoZXsnv4RNJVx+SeScYPKgmOuyM3m9\nk0lQ4uefTWfaQtqHKJjsXFWxAMCvATQDOBPA4wB+ZzKoQleobadOXdGQacJl94u8EBK5+G3h5BUm\nThxwc3m9E0lQss8/00pbIexDFHx2BoD6P1U9QUTeVdVvxj+Xlwi/jKOgrqrw40BEbstmkCw77d6F\ncnMwp/oAJF7pYGq5mcjliiQnPn839yFeVUFOs1Nx+CJyr4oPROR6EbkIQKnhuAoe204zl0lJOpOq\njt/7X7RXUXGywhV/Rm1qudnIpQ+OE5+/3/chonh2EocbAJQA+D6AEwBMBnCFyaDIuTKxV9pT8xWH\n3YQr0y/y9pabbP28sO3tHHCdOKglSxIm/fckDOo+yPHlutFs50Qiz5MBCgo7iUOLqu5V1c2q+u+q\nWqOqK41HVuCcGLPBK+2p+YzD1NgT6ZabbP3c3vaZHnBzPailSj5uOuUmI8vN95m6E4m8H0clJUrG\nTh+HFQD6wLqPxO9V9e/5CCxJHAXVxyEXXmmTz0ccbrV7A8nXb0jPIRAI3v/X+0a3vZ31zqSfjBOj\nkiYbabTxQKOR5bY3RogXRpT0CvZxIKe1W3FQ1TNhXU1RD+AhEXlXRP6f8cgoa145S2svjlxL+W62\newPJ1++oHkdhSI8hRre93fXOpIrgRIUr2Rm1qeWm43a1hyjoMrpXhYh8E8BMAN9W1WJjUSV/b1Yc\nMpDNWVo+48imp330LNIrFRUg+foBMLLtM11vU/c28Sov7RdewooDOc3OAFBHi0itiLwL4H4AfwPg\nvTu5UBteaU9NjGP7vu1ZdXaLP4v0SkUFSL6dTW37TNfb1L1NvCrb/cILnViJ/MROH4fXASwGsERV\nt+YlquRxsOIQEJm0vac6i/zF6F+g/8/7Ayisdm+vVJK8KtPt4/U7WzqBFQdymp0+DiNU9RduJg0U\nLE4MF7z34N6CbPf2SiXJq+xuH69c5knkRxn1cXATKw7BkWnbey5n2V5s9w5C5SMICmV0VlYcyGl2\nxnEgclSmbe+5nGV7qT8EEJzKRxBwQCai7GR6VUUfVf3UYDzp3psVB8qKF/oFeLHyUegK5aoTVhzI\naZlWHJYZiYLIIC/0C/Ba5YMK76oTIqdkWnFYrarHGYwn3Xuz4kC+5oXKBxUeVhzIaZlWHH5jJAqi\nApBN5YNjDBCR1/CqCsoarw5wTrJtWQhjDJB5rDiQ03hVRYDk8+yUVwc4J3FbcowBIvIy44mDiIwR\nkQ0islFEZiX5+4UislZEVovImyJyqumYgihfB3Ie1JyTalsKpKA7UrJ5hsjbUiYOIjIgzd9G2lm4\niBQBeADAaADfADBJRBK7Lr+kqsdEOl1eCeARO8smS74P5Lw6wDnptmWhjjHAShaR96WrOIRFZKaI\ndIg+ISJzsP7HAAAQfElEQVRVIvI7AD+zufyTAXygqptU9SCse16Mj59BVT+PmywF0Gpz2QR3DuSF\nelAzIdW29MIlpPnEShaRf6RLHE4AcCSANSJylojcAOBNAK/DSgjs6AugLm56c+S5NkRkgoisB/AM\ngO/aXDZF5PtAXmgHNZNSbctCG2OAlSwi/7Bzd8wbYFUYtgIYrqq2a4giUgNgtKpeE5meDOBkVf1+\nivlPAzBHVUcl+ZvOmTMnNh0KhRAKheyGEmiFMgIeBRvHuXBGOBxGOByOTd922228qoIclTJxEJHu\nAOYCGAZgJoBxAM4GcIOq/tnWwkWGA6hV1TGR6ZsBqKqmvJOMiPwTwEmq2pDwPC/HJAowJsBm8HJM\nclq6xOFDAL8C8HNVbY48d2zkuU2qOqndhVv9I96HlXBsg9XUMUlV18fNc6Sq/jPy+/EAnlbV/kmW\nxcSBiChDTBzIaR3T/O30xGYJVV0D4BQRudrOwlW1RUSuB/ACrP4U81V1vYhca/1ZHwZQIyJTATQB\n2A/gW9msCBEREZnHkSOJiAKMFQdyGkeOJCIiItuYOBAREZFtTByIiIjINiYORESJ6uuBVausn0TU\nBhMHIqJ4ixYB1dXAqFHWz0WL3I6IyFN4VQURUVR9vZUs7N//5XNduwKbNgGVle7FlQNeVUFOY8WB\niCjq44+B4uK2z3XqZD1PRACYOBAdiu3bhWvgQKCpqe1zBw9az2eK+xEFFBMHonhs3y5slZXA/PlW\n80R5ufVz/vzMmym4H1GAsY8DUVQA27cpS/X1VvPEwIGZf/Ye24/Yx4GcxooDURTbtymqshI46aTs\nDvTcjyjgmDgQRTnZvk2Fi/sRBRwTB6Iop9q3qbBxP6KAYx8HokS5tG8TRXlkP2IfB3IaEwciogBj\n4kBOY1MFERER2cbEwW84qAx5AfdDooLFxMFP3BhUhgcIf8jn58TBjYgKGvs4+IXpQWWSdeRatAi4\n8krrmvSmJqtn+KRJub8XOSufn5OXBjfySOdDr2MfB3IaKw5+YXJQmWRnkPX11sFo/35g1y7r55VX\nsvLgNfn+nLwyuBGrHkSuYeLgF6YGlUl14Fm92hsHCEov3wdyLwxuxKSWyFVMHPzC1KAyqQ48gDsH\nCPapyEy+D+ReGNzIK1UPogLFPg5+43S7bro265dess7kOnWyDkam+ziwT0V2otstX58T4G7/Ai/1\ns/AB9nEgpxlPHERkDICfw6puzFfVuQl/vwzArMjkHgDfU9V3kyyHiYMp6Q48+TpA8GCQm0LrKOhG\nsuQEFz4nJg7kNKOJg4gUAdgI4GwAWwGsAnCpqm6Im2c4gPWquiuSZNSq6vAky2Li4JRkX15uH3hW\nrbI6uu3a9eVz5eVW1eOkk/IfT67c3p5B44V9Ntf3c6mixsSBnGa6j8PJAD5Q1U2qehDAYgDj42dQ\n1ZWqGj1arATQ13BMhS1Vb/RcbiMclUv/BC90unNKkHr8e6HPicl9NtcY7GKHTgoQ04lDXwB1cdOb\nkT4xuArAcqMRFTKTX165frE62enOzYNdkA4Q2XymTm97L2xPJ2Jgh04KkI5uBxAlImcC+HcAp6Wa\np7a2NvZ7KBRCKBQyHlegRL+84vsRRL+8cq00RL9Yo8u+8krgnHOs5dot8U6aZL0mn+Vgp8vdprZx\nPsRvCyD9Z5qMiVK8F7anEzHksaIWDocRDocdXy5RjKoaewAYDuC5uOmbAcxKMt9QAB8AODLNspRy\ntGOHateuqsCXj65dredz8eabqt26tV1uebn1/MKF1nt062b9XLjQmXVJJtP1yzS2HTusdUq3vUxt\nY9MSt8Udd6T+TJPJZr3zvT3tvJ/JGKLbuLzc/P9CnMh3p9Hvej4K62F24UAHAP8AUA2gGMAaAEcn\nzDMgkjQMb2dZSg4w8eWV6ov1vffyexBNl8DYjdmJJMOlA4RtiQfQZNuiS5fMtk8m2141/9sz1wQ2\nXQyZJCTZJi85YOLAh9MP828AjAHwfiQ5uDny3LUAron8/hsA/wLwNoDVAN5MsRwtOKa+ZEwsN9kX\na6YHk1xlkgyYTDKir8nzAcKWZAfQVNvijjvsH7DTbSM7iYrJ7elUxSBZDPmsqGWJiQMfTj9cD8B2\noIWWOPjgC+kQThwgcmX37NTJJMOrSYJq29iyqQxlsm7Jtn0miYqphNLU+/mkWYqJAx9OP1wPwHag\nhZQ4+OQLyRY3yvZ2D3ZOJBluJHiZrp+dfgtOfU65JiomuNG3x0OYOPDh9MP1AGwHWkiJg0++kGxL\ndqDL91l6qvfLJclws6LSXqKSTb8Fpz8TO51m85VQ5rNvj8cSfCYOfDj9cD0A24EWUuJgsk3WC5w6\nS8/2zNup98smwctXW70T/RZy1V68Xkgec43B6x1hVZk48OH4w/UAbAdaSImDau5fSF7tI+H0pW3Z\nnHk7dVZo+vLPRE516MznAdvLB9Z8J7AuYeLAh9MP1wOwHWihJQ6q7l93bkK2zTB22s5zvXoiGyY6\nY6aSbaLi9kHbiwdWL/+POIyJAx9OPzwzciQlUVmZ3eh4XhhtL5VsRtBLHJHwv/7L/vqZHrHP7miX\nTnwm0WG5E+8Kmer1TozE6YRs9+Ns2RkN1Mv/I0QeZ/y22k7h3TEz4PVbVGdyS+Rk69KlCyBif/28\ncAtmJz8T3nkzNbvDXnv9f8RBvDsmOY2JQ1B54WCZjt2DX6rbbd90E/DjH9tfPy8cbL3+mfhdpslA\ngXweTBzIaUwcgswLB8tcpTsYAP5bvyB8Jl6VKsl86SXr9tvJFMDnwcSBnMbEgbyvQM4MCbkdyAuo\n+SETTBzIaUwcyB8K4Myw4DlxW24mmYdg4kBOY+JA2eGB3LxC2sbsPGoMEwdyWpHbAZAPLVpkfcmP\nGmX9XLTI7YiCp9C2cfTyyHjRyyMzVVlp9Wlg0kBkBCsOlBm2I5tXiNu4ENc5T1hxIKex4kCZcfLM\nkJIrxG0cHdyqa1frSoiuXdMPbkVErmHFgTLDM0PzCnkbs3+C41hxIKex4kCZ4Zmh8+rrrTEI6uut\n6ULexuyfQOR5rDhQdnhm6Ix0lyByG5MDWHEgpzFxIHJLITdJUN4wcSCnsamCyC2F2AmSiHyPiUNQ\nJLaTk/eZvuU3EZEBTByCoNAGCwqKQu4ESUS+xT4Ofsd2cv9jJ0gyiH0cyGnGKw4iMkZENojIRhGZ\nleTvXxWRv4nIARGZYTqewGE7uf/xEkQi8pGOJhcuIkUAHgBwNoCtAFaJyNOquiFutn8BmA5ggslY\nAovt5ERElEemKw4nA/hAVTep6kEAiwGMj59BVXeq6v8BaDYcSzCxnZyIiPLIaMUBQF8AdXHTm2El\nE+SkSZOAc85hOzkRERlnOnFwVG1tbez3UCiEUCjkWiyeU1nJhIGIEA6HEQ6H3Q6DAszoVRUiMhxA\nraqOiUzfDEBVdW6SeecA2KOqP02xLF5VQUSUIV5VQU4z3cdhFYDBIlItIsUALgWwNM383LmJiIg8\nzPg4DiIyBsAvYCUp81X1bhG5Flbl4WERqQLwFoAyAK0A9gL4uqruTVgOKw5ERBlixYGcxgGgiIgC\njIkDOY1DThMREZFtTByIiIjINiYOREREZBsTByIiIrKNiQMRERHZxsSBiIiIbGPiQERERLYxcSAi\nIiLbmDgQERGRbUwciIiIyDYmDkRERGQbEwciIiKyjYkDERER2cbEgYiIiGxj4kBERES2MXEgIiIi\n25g4EBERkW1MHIiIiMg2Jg5ERERkGxMHIiIiso2JAxEREdnGxIGIiIhsM544iMgYEdkgIhtFZFaK\neeaJyAciskZEjjUdExEREWXHaOIgIkUAHgAwGsA3AEwSka8lzDMWwJGqOgTAtQAeNBmTV4XDYbdD\nMIrr519BXjcg+OtH5DTTFYeTAXygqptU9SCAxQDGJ8wzHsDjAKCqbwDoJiJVhuPynKB/eXH9/CvI\n6wYEf/2InGY6cegLoC5uenPkuXTzbEkyDxEREXkAO0cSERGRbaKq5hYuMhxAraqOiUzfDEBVdW7c\nPA8CWKGqv49MbwBwhqpuT1iWuUCJiAJMVcXtGCg4Ohpe/ioAg0WkGsA2AJcCmJQwz1IA0wD8PpJo\nfJaYNADc8YmIiLzAaOKgqi0icj2AF2A1i8xX1fUicq31Z31YVZeJyDgR+QeAfQD+3WRMRERElD2j\nTRVEREQULL7qHCkil4jI30WkRUSOdzseJ9gZIMvPRGS+iGwXkXfcjsVpItJPRP4sIutE5F0R+b7b\nMTlJRDqLyBsisjqyfnPcjslpIlIkIm+LyFK3YzFBRD4WkbWRz/BNt+OhYPBV4gDgXQAXAXjF7UCc\nYGeArABYAGv9gqgZwAxV/QaAEQCmBenzU9UvAJypqscBOBbAWBE52eWwnHYDgPfcDsKgVgAhVT1O\nVYP22ZFLfJU4qOr7qvoBgKB0lLQzQJavqepfADS6HYcJqvqpqq6J/L4XwHoEbAwSVf088mtnWH2i\nAtO2KSL9AIwD8IjbsRgk8Nn3PHkfdyh32Rkgi3xARAbCOit/w91InBUp5a8G8CmAF1V1ldsxOehn\nAG5CgJKhJBTAiyKySkSudjsYCgbTl2NmTEReBBA/5LTA2vlnq+oz7kRFlJqIlAL4I4AbIpWHwFDV\nVgDHiUg5gP8Vka+rqu9L+yJyHoDtqrpGREIIThUz0amquk1EKmElEOsjVUCirHkucVDVUW7HkEdb\nAAyIm+4XeY58QkQ6wkoanlDVp92OxxRV3S0iKwCMQTD6BJwK4EIRGQegK4AyEXlcVae6HJejVHVb\n5Ge9iPwPrOZRJg6UEz83VQThDCE2QJaIFMMaICuIvbsFwfi8knkUwHuq+gu3A3GaiPQSkW6R37sC\nGAVgg7tROUNV/0tVB6jqEbD+7/4ctKRBREoi1TCIyFcAnAvg7+5GRUHgq8RBRCaISB2A4QCeFZHl\nbseUC1VtARAdIGsdgMWqut7dqJwlIgsB/A3AUSLyiYgEZoAvETkVwOUAzopc7va2iIxxOy4HHQZg\nhYisgdV343lVXeZyTGRfFYC/RPqorATwjKq+4HJMFAAcAIqIiIhs81XFgYiIiNzFxIGIiIhsY+JA\nREREtjFxICIiItuYOBAREZFtTByIiIjINiYORA4SkT1ux0BEZBITByJncWAUIgo0z92rgshLROQu\nAHWq+qvI9BwAzQDOBNAdQCcAt6rq0oTX9QHwewBlsP7Pvqeqf81n7EREJnDkSKI0RORYAD9X1VBk\neh2sMf93qepeEekJYKWqDon8fbeqlovIDACdVfUuEREAJaq6z6XVICJyDCsORGlEbrtcGakg9AbQ\nAOBTAL8QkZEAWgEcLiK9VXVH3EtXAZgvIp0APK2qa/MePBGRAezjQNS+JQAmAvg2rOaHyQB6AjhO\nVY8DsANAl/gXqOprAE6HdZv0x0Rkcl4jJiIyhIkDUfv+AOvWyzWwkohuAHaoaquInAmgOm5eAQAR\nGRCZZz6ARwAcn9+QiYjMYFMFUTtU9T0RKQOwWVW3i8iTAJ4RkbUA3gIQfyv0aKehEICbROQggD0A\npuYzZiIiU9g5koiIiGxjUwURERHZxsSBiIiIbGPiQERERLYxcSAiIiLbmDgQERGRbUwciIiIyDYm\nDkRERGQbEwciIiKy7f8DRcCtxlsLYn0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3481e79278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(sky_vals, sky_mean - sky_vals, color='b', marker='x', label='mean')\n",
    "plt.scatter(sky_vals, sky_med - sky_vals, color='g', marker='*', label='median')\n",
    "plt.scatter(sky_vals, sky_ofil - sky_vals, color='r', marker='o', label='ofilter (Python)')\n",
    "plt.xlabel('vals')\n",
    "plt.ylabel('X - vals')\n",
    "plt.axhline(0, color='k', linestyle='--')\n",
    "plt.legend(bbox_to_anchor=(1.01, 1), loc='upper left', scatterpoints=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display histogram of the data from the last iteration above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7f3481d2c400>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE5tJREFUeJzt3X+MXWWdx/H3ty1toaVAlBaloBhBwahoNvUHMXtdWBR0\nW1wDomxX7KqJiIJE0hYTGTauBbMbdd1lE3+sKQJbii5Sf+xSSnM1FbSKBVnLlipbxGoHFSw2tkNn\n+t0/7kGmpdO5c+dOz/SZ9yuZcOa5z3POlyfMh3Ofc+65kZlIkso1qe4CJEljy6CXpMIZ9JJUOINe\nkgpn0EtS4Qx6SSpcW0EfEVsi4v6I2BAR66u2YyJidURsiog7IuKoQf2XRsTmiHgwIs4eq+IlScNr\n94x+D9DIzFdl5ryqbQmwJjNfAqwFlgJExGnABcCpwDnA9RER3S1bktSudoM+9tN3AbC82l4OnFdt\nzwdWZGZ/Zm4BNgPzkCTVot2gT+DOiPhhRLy3apuTmb0AmbkNmF21Hw88Omjs1qpNklSDKW32OyMz\nfx0RxwKrI2ITrfAfzGcpSNI41FbQZ+avq3/+JiK+Tmsppjci5mRmb0QcBzxWdd8KnDBo+NyqbS8R\n4f8YJKkDmTmi657DLt1ExBERMbPangGcDTwArAIurrq9G7i92l4FXBgRUyPiJODFwPohivUnk6uv\nvnrU+zjgfO7nNXoOPP9w6M5FKT/OhXOxv59OtHNGPwe4rToDnwLclJmrI+JHwMqIWAQ8QutOGzJz\nY0SsBDYCu4FLstPqJEmjNmzQZ+b/Aafvp/1x4KwhxiwDlo26OknSqPnJ2HGg0WjUXcK44Vw8w7l4\nhnMxOlHXqkpEuKLTRREx9PpdBOzzWlwT5NVDz/9+hkgaB6q/9e5ejJUkHdoMekkqnEEvSYUz6CWp\ncAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMK1++XgqsGNN97M\nBz5wWdv9jzzy2L1+P+KIGdx77zrmdrswSYcUg34c27hxIzt2vBe4oo3es9mxY+NeLRFnsW3bNoNe\nmuAM+nFvBnDssL1a9u43adJhXa9G0qHHNXpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn\n0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqXNtBHxGTIuLHEbGq\n+v2YiFgdEZsi4o6IOGpQ36URsTkiHoyIs8eicElSe0ZyRn8ZMPhLSZcAazLzJcBaYClARJwGXACc\nCpwDXB8R0Z1yJUkj1VbQR8Rc4Fzgi4OaFwDLq+3lwHnV9nxgRWb2Z+YWYDMwryvVSpJGrN0vB/80\ncCVw1KC2OZnZC5CZ2yJidtV+PHDPoH5bqzbVYP369fwZcNdddz3rtf21PeNMHnroIU455ZQxq03S\nwTFs0EfEW4DezLwvIhoH6JojPXhPT8+fthuNBo3GgXavkXrqqfO56qqvcQnw9rd/cu8XP7Kftr2c\nyete90Z+97utY1mipGE0m02azeao9tHOGf0ZwPyIOBc4HDgyIr4CbIuIOZnZGxHHAY9V/bcCJwwa\nP7dqe5bBQa/u27lzMTt3LgaC7dv3PXvfX9ve+vr6xqw2Se3Z9yT4mmuuGfE+hl2jz8yrMvPEzHwR\ncCGwNjMXAt8ALq66vRu4vdpeBVwYEVMj4iTgxcD6EVcmSeqKdtfo9+daYGVELAIeoXWnDZm5MSJW\n0rpDZzdwSWaOeFlHktQdIwr6zPwO8J1q+3HgrCH6LQOWjbo6SdKo+clYSSqcQS9JhTPoJalwBr0k\nFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1Lh\nRvMNUxrGnj17WLFiBTt37uxo/P333we8trtFSZpwDPoxdM8997Bo0UeYPPmtHe5hDrCgmyVJmoAM\n+jGUmUyffjLbt3+p7lIkTWCu0UtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BL\nUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCjds0EfEtIj4QURsiIgHIuLqqv2YiFgd\nEZsi4o6IOGrQmKURsTkiHoyIs8fyX0CSdGDDBn1m9gFvzMxXAacD50TEPGAJsCYzXwKsBZYCRMRp\nwAXAqcA5wPUREWNUvyRpGG0t3WTmH6vNabS+fjBpfZnp8qp9OXBetT0fWJGZ/Zm5BdgMzOtWwZKk\nkWkr6CNiUkRsALYBd2bmD4E5mdkLkJnbgNlV9+OBRwcN31q1SZJq0O4Z/Z5q6WYuMC8iXkbrrH6v\nbt0uTpI0elNG0jkzn4yIJvBmoDci5mRmb0QcBzxWddsKnDBo2Nyq7Vl6enr+tN1oNGg0GiMpR5KK\n12w2aTabo9pHZB74RDwingvszsztEXE4cAdwLfDnwOOZeV1ELAaOycwl1cXYm4DX0FqyuRM4Ofc5\nUETs21ScdevW8da3LmH79nUH4WjBUG+qkiD2fa0noOfA8z9jxnPZseO3XapPUjdEBJk5ohtc2jmj\nfx6wPCIm0VrquSUzvx0R3wdWRsQi4BFad9qQmRsjYiWwEdgNXFJ8okvSODZs0GfmA8Cr99P+OHDW\nEGOWActGXZ0kadRGtEaviWXXrh28/vXndjz+gx98DxdddH4XK5LUCYNeQxoY+C733NPpGv1ajj32\ndoNeGgcMeh3AaD7n9jjw390qRNIo+FAzSSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BL\nUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwvn0So2RSaxf/z0WLPibjkZPnjyJz3zmE5x44oldrkua\neAx6jZG3s20brFq1p6PR06b9G29723dYuHBhl+uSJh6DXmNkGvCujkcfdtgd3StFmuBco5ekwhn0\nklQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9J\nhTPoJalwBr0kFc4vHhnGqlWr+NnPftbR2IcffrjL1UjSyA0b9BExF7gBmAPsAb6Qmf8cEccAtwAv\nALYAF2Tm9mrMUmAR0A9clpmrx6b8sXf++e8kcxFwWAejp9Lfv7jbJUnSiLRzRt8PXJGZ90XETODe\niFgNvAdYk5mfiojFwFJgSUScBlwAnArMBdZExMmZmWP07zCmMmH37muBGXWXIkkdGXaNPjO3ZeZ9\n1fYO4EFaAb4AWF51Ww6cV23PB1ZkZn9mbgE2A/O6XLckqU0juhgbES8ETge+D8zJzF5o/c8AmF11\nOx54dNCwrVWbJKkGbQd9tWzzVVpr7juAfZdiDsmlGUkqXVt33UTEFFoh/5XMvL1q7o2IOZnZGxHH\nAY9V7VuBEwYNn1u1PUtPT8+fthuNBo1GY0TFS1Lpms0mzWZzVPuIdq6RRsQNwG8z84pBbdcBj2fm\nddXF2GMy8+mLsTcBr6G1ZHMn8KyLsRFxSFyfnTp1Brt3P8b4vxgbDPWmKgli39d6AnrG7/zPnLmQ\n668/m4ULF9ZdijSuRASZGSMZ087tlWcAFwEPRMQGWmlyFXAdsDIiFgGP0LrThszcGBErgY3AbuCS\nQyLRJalQwwZ9Zn4PmDzEy2cNMWYZsGwUdUmSusRHIEhS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TC\nGfSSVDi/eETj1lNPPcXOnTs7Gjt58mSmTp3a5YqkQ1Nbj0AYkwP7CIQuK+sRCNOm/T0DA51/5m76\n9CP41a+2cOSRR3axKql+Y/IIBKkOfX0fBz7e8fj+/mPZtWuXQS/hGr0kFc+gl6TCGfSSVDiDXpIK\nZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAG\nvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFW5K3QVIY+XWW29l1qxZHY19\n5Stfyctf/vIuVyTVY9igj4gvAW8FejPzFVXbMcAtwAuALcAFmbm9em0psAjoBy7LzNVjU7o0tMyP\nsnjx3R2NHRj4HSec8Fk2bfphl6uS6tHOGf2Xgc8BNwxqWwKsycxPRcRiYCmwJCJOAy4ATgXmAmsi\n4uTMzC7XLR1QX99i+vo6HX0vAwPv72Y5Uq2GXaPPzHXAE/s0LwCWV9vLgfOq7fnAiszsz8wtwGZg\nXndKlSR1otOLsbMzsxcgM7cBs6v244FHB/XbWrVJkmrSrYuxHS3N9PT0/Gm70WjQaDS6VI4klaHZ\nbNJsNke1j06Dvjci5mRmb0QcBzxWtW8FThjUb27Vtl+Dg16S9Gz7ngRfc801I95Hu0s3Uf08bRVw\ncbX9buD2Qe0XRsTUiDgJeDGwfsRVSZK6pp3bK28GGsBzIuIXwNXAtcCtEbEIeITWnTZk5saIWAls\nBHYDl3jHjSTVa9igz8x3DfHSWUP0XwYsG01RkqTu8REIklQ4g16SCmfQS1LhDHpJKpxBL0mFM+gl\nqXAGvSQVzi8ekfZj+/bf8IUvfKHj8QsWLGD27NnDd5QOAoNeepaX8oc//DWXX97Z0zt2797Aww8/\nwrJln+hyXVJnDHrpWWbQ1/eZUYz/BzL/2LVqpNFyjV6SCmfQS1LhDHpJKlzxa/Tf+tYdrFz59Y7H\n9/d3/A3TkjQuFB/0n/7057nrrucAr+pwD18DZnSxIkk6uIoP+pY3AW+vuwhJqoVr9JJUOINekgpn\n0EtS4Qx6SSqcQS9JhZsgd91IB9fq1XeyY8fOjsYefvh0eno+xowZ3tar7jDopa77OzZsOJwNGzob\nPX3655g//xze8IY3dLcsTVgGvdR1xwFXdDx62rTbuleKhGv0klQ8g16SCmfQS1LhXKOXxpmBgelc\nfPGlzJx5dC3Hf9/7FnLppe+t5dgaGwa9NM7s2HEDO3Y8VNPRv8ttt33boC+MQS+NO8+rfurwO6DD\n+0I1brlGL0mFM+glqXAGvSQVzjV6SYPModn8JtOnH9nR6Ijgm9+8jTPPPLPLdWk0DHpJg5zBnj1P\n0Ne3p6PRhx/+IX7+858b9OPMmAV9RLwZ+Ayt5aEvZeZ1neznpz/9Ke9//0cZGMiO6njwwQ3AxR2N\nlSam0Tw1c2rXqlD3jEnQR8Qk4F+AM4FfAT+MiNsz839Huq+7776be+9N+vou77CayVUZ41kTaNRc\nw3jRxLl4WhPnoqXZbNJoNOou45A1Vmf084DNmfkIQESsABYAIw56gMmTTwDe3L3qxp0m/kE/rYlz\n8bQmh9pcDAwcx6WXXs7lly/uaPzOnb/nyiuvZNasWXu1r127lnXr1nWjxAN6xzvewcknn9zR2IGB\nAfr6+jo+9pQpU5g6dWzeEY1V0B8PPDro91/SCn9JBXvqqauBj7B7d6d7WMU//uNDwK69WjP7aTZ3\n7X9Il2T+gBtvXMm8ea/vaPzXv34Lf/jD75ky5fBOjs7RRz+H3/zmlx0dezjj/mLsYYcdxp49dzBr\n1l/VXcqY2bVrE9On3zuqfTz5JEPP0X5ee5ID9AeefPIbtcx5N+aiFM7FM3bt2sb06feP6TEGBnax\nadMDbNr0QMf7mDXrTcBhHYxMnnhidcfHHU5kdnaR84A7jXgt0JOZb65+XwLk4AuyEdH9A0vSBJCZ\nMZL+YxX0k4FNtK6C/hpYD7wzMx/s+sEkSQc0Jks3mTkQEZcCq3nm9kpDXpJqMCZn9JKk8eOgP+sm\nIuZGxNqI+GlEPBARHz7YNYw3ETEpIn4cEavqrqVOEXFURNwaEQ9W/328pu6a6hIRH4mI/4mIn0TE\nTRExYT6JFBFfiojeiPjJoLZjImJ1RGyKiDsi4qg6azxYhpiLT1V/I/dFxNciYtaB9gH1PNSsH7gi\nM18GvA74YES8tIY6xpPLgI11FzEOfBb4dmaeCrwSmJDLfRHxfOBDwKsz8xW0llgvrLeqg+rLwJv2\naVsCrMnMlwBrgaUHvap67G8uVgMvy8zTgc20MRcHPegzc1tm3ldt76D1x3z8wa5jvIiIucC5wBfr\nrqVO1VnJGzLzywCZ2Z+ZT9ZcVp0mAzMiYgpwBK1PmE8ImbkOeGKf5gXA8mp7OXDeQS2qJvubi8xc\nk5lPP4zo+8Dc4fZT62OKI+KFwOnAD+qso2afBq4EJvrFkpOA30bEl6tlrM9HRCefPDnkZeavgH8C\nfgFsBX6fmWvqrap2szOzF1oni8DsmusZLxYB/zVcp9qCPiJmAl8FLqvO7CeciHgL0Fu9w4nqZ6Ka\nArwa+NfMfDXwR1pv1yeciDia1hnsC4DnAzMj4l31VjXuTPQTIyLiY8DuzLx5uL61BH31dvSrwFcy\n8/Y6ahgnzgDmR8TDwH8Ab4yIG2quqS6/BB7NzB9Vv3+VVvBPRGcBD2fm45k5APwn0Nnn8svRGxFz\nACLiOOCxmuupVURcTGvJt60TgLrO6P8d2JiZn63p+ONCZl6VmSdm5otoXWxbm5l/W3dddajelj8a\nEadUTWcycS9Q/wJ4bURMj4igNRcT7cL0vu9wV/HM88bfDUykE8S95qJ6BPyVwPzMbOspanXcXnkG\ncBHwFxGxoVqPLfnRlGrfh4GbIuI+WnfdfLLmemqRmetpvaPZANxP64/887UWdRBFxM3A3cApEfGL\niHgPcC3wlxHx9Cfur62zxoNliLn4HDATuLPKz+uH3Y8fmJKksvnl4JJUOINekgpn0EtS4Qx6SSqc\nQS9JhTPoJalwBr0kFc6gl6TC/T/5Kx9j15F30gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3481d49d68>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h = plt.hist(data.flatten(), bins=20, histtype='stepfilled')\n",
    "plt.axvline(msky, color='r')\n",
    "plt.axvline(np.median(data), color='g')\n",
    "plt.axvline(data.mean(), color='b')\n",
    "plt.axvline(val, color='k')"
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

Hmm....The offset is worrisome, although agreement with IRAF seems good for the one iraf datapoint.
Can you plot X-vals vs. vals?

@hcferguson , I added the plot as you requested. I also added a new section at the bottom using a different skewed distribution, and for that perhaps ofilter is better?