psf_shape_simon.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to make a notebook that looks at the PSF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.daf.persistence import Butler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "butler = Butler('/home/shared/twinkles/output_data_v2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset = butler.subset('calexp', filter='i')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataId = subset.cache[6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_calexp = butler.get('calexp', **dataId)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "psf = my_calexp.getPsf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "shape = psf.computeShape()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.650111999169302"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shape.getDeterminantRadius()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.afw.geom import Point2D\n",
    "point = Point2D(50.1, 160.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "shape = psf.computeShape(point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.6500495193183045"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shape.getDeterminantRadius()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "shapes = []\n",
    "for x in range(100):\n",
    "    for y in range(100):\n",
    "        point = Point2D(x*40., y*40.)\n",
    "        shapes.append(psf.computeShape(point).getDeterminantRadius())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pylab as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([3879., 1679., 1000.,  831.,  718.,  644.,  585.,  399.,  197.,\n",
       "          68.]),\n",
       " array([2.64647811, 2.64928225, 2.65208638, 2.65489051, 2.65769464,\n",
       "        2.66049878, 2.66330291, 2.66610704, 2.66891118, 2.67171531,\n",
       "        2.67451944]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFLhJREFUeJzt3X+s3fV93/HnK4YQRkKB+cZybTMTyalkkOKUK5cpUZUVtbikmsn+QKZq8B8Mp4KlidRtgu6PUk2W6NQkK9pAIgnDdFmo1STDamCIeNmiSAN6yRzAJixeMcNXxnbTRCT7g87Oe3+cD8rJ5Zp77g/fw/Xn+ZCOzue8v5/v93w++tp63e/3fM/5pqqQJPXrHeMegCRpvAwCSeqcQSBJnTMIJKlzBoEkdc4gkKTOGQSS1DmDQJI6ZxBIUufOG/cA5rJ69erauHHjuIchSSvG6tWrefzxxx+vqm2j9B85CJKsAqaA6ar6rSSXAX8ObASOADdW1Q9b3zuBW4DTwO9V1eOtfjXwIHAh8CjwqZrjNy42btzI1NTUqMOUJAFJVo/adz6nhj4FvDD0+g5gf1VtAva31yTZDOwArgS2Afe2EAG4D7gV2NQeI6WVJOnsGSkIkqwHPgp8Yai8HdjT2nuAG4bqD1fV61X1EnAY2JpkLXBxVT3ZjgIeGlpHkjQmox4R/FvgXwI/Haqtqapjrf0qsKa11wGvDPU72mrrWntm/U2S7EoylWTq5MmTIw5RkrQQcwZBkt8CTlTVM2fq0/7CX7Lfs66q+6tqsqomJyYmlmqzkqRZjPJh8YeAf5zkeuBdwMVJ/iNwPMnaqjrWTvucaP2ngQ1D669vtenWnlmXJI3RnEcEVXVnVa2vqo0MPgT+r1X1O8A+YGfrthN4pLX3ATuSXJDkCgYfCj/dTiO9luSaJAFuHlpHkjQmi/kewd3A3iS3AC8DNwJU1cEke4FDwCng9qo63da5jZ9dPvpYe0iSxihv91tVTk5Olt8jkKT5SfJMVU2O0tefmJCkzr3tf2JiMTbe8fWxvO+Ruz86lveVpIXwiECSOmcQSFLnDAJJ6pxBIEmdMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjpnEEhS5wwCSeqcQSBJnTMIJKlzBoEkdc4gkKTOzRkESd6V5Okk301yMMkftfpdSaaTHGiP64fWuTPJ4SQvJrluqH51kufasnvavYslSWM0yo1pXgd+rap+kuR84NtJ3rjX8Oeq6k+GOyfZzOAm91cCvwh8I8n7232L7wNuBZ4CHgW24X2LJWms5jwiqIGftJfnt8db3eh4O/BwVb1eVS8Bh4GtSdYCF1fVkzW4UfJDwA2LG74kabFG+owgyaokB4ATwBNV9VRb9MkkzyZ5IMmlrbYOeGVo9aOttq61Z9YlSWM0UhBU1emq2gKsZ/DX/VUMTvO8D9gCHAM+s1SDSrIryVSSqZMnTy7VZiVJs5jXVUNV9SPgm8C2qjreAuKnwOeBra3bNLBhaLX1rTbd2jPrs73P/VU1WVWTExMT8xmiJGmeRrlqaCLJJa19IfDrwPfaOf83fAx4vrX3ATuSXJDkCmAT8HRVHQNeS3JNu1roZuCRJZyLJGkBRrlqaC2wJ8kqBsGxt6r+MsmfJdnC4IPjI8AnAKrqYJK9wCHgFHB7u2II4DbgQeBCBlcLecWQJI3ZnEFQVc8CH5yl/vG3WGc3sHuW+hRw1TzHKEk6i/xmsSR1ziCQpM4ZBJLUOYNAkjpnEEhS5wwCSeqcQSBJnTMIJKlzBoEkdc4gkKTOGQSS1DmDQJI6ZxBIUucMAknqnEEgSZ0zCCSpcwaBJHXOIJCkzo1y8/p3JXk6yXeTHEzyR61+WZInkny/PV86tM6dSQ4neTHJdUP1q5M815bd025iL0kao1GOCF4Hfq2qPgBsAbYluQa4A9hfVZuA/e01STYDO4ArgW3Ave3G9wD3AbcCm9pj2xLORZK0AHMGQQ38pL08vz0K2A7safU9wA2tvR14uKper6qXgMPA1iRrgYur6smqKuChoXUkSWMy0mcESVYlOQCcAJ6oqqeANVV1rHV5FVjT2uuAV4ZWP9pq61p7Zl2SNEYjBUFVna6qLcB6Bn/dXzVjeTE4SlgSSXYlmUoydfLkyaXarCRpFvO6aqiqfgR8k8G5/ePtdA/t+UTrNg1sGFptfatNt/bM+mzvc39VTVbV5MTExHyGKEmap1GuGppIcklrXwj8OvA9YB+ws3XbCTzS2vuAHUkuSHIFgw+Fn26nkV5Lck27WujmoXUkSWNy3gh91gJ72pU/7wD2VtVfJvkfwN4ktwAvAzcCVNXBJHuBQ8Ap4PaqOt22dRvwIHAh8Fh7SJLGaM4gqKpngQ/OUv8BcO0Z1tkN7J6lPgVc9eY1JEnj4jeLJalzBoEkdc4gkKTOGQSS1DmDQJI6ZxBIUucMAknqnEEgSZ0zCCSpcwaBJHXOIJCkzhkEktQ5g0CSOmcQSFLnDAJJ6pxBIEmdMwgkqXMGgSR1bpSb129I8s0kh5IcTPKpVr8ryXSSA+1x/dA6dyY5nOTFJNcN1a9O8lxbdk+7ib0kaYxGuXn9KeD3q+o7Sd4DPJPkibbsc1X1J8Odk2wGdgBXAr8IfCPJ+9sN7O8DbgWeAh4FtuEN7CVprOY8IqiqY1X1ndb+MfACsO4tVtkOPFxVr1fVS8BhYGuStcDFVfVkVRXwEHDDomcgSVqUeX1GkGQj8EEGf9EDfDLJs0keSHJpq60DXhla7WirrWvtmfXZ3mdXkqkkUydPnpzPECVJ8zRyECR5N/AV4NNV9RqD0zzvA7YAx4DPLNWgqur+qpqsqsmJiYml2qwkaRYjBUGS8xmEwJeq6qsAVXW8qk5X1U+BzwNbW/dpYMPQ6utbbbq1Z9YlSWM0ylVDAb4IvFBVnx2qrx3q9jHg+dbeB+xIckGSK4BNwNNVdQx4Lck1bZs3A48s0TwkSQs0ylVDHwI+DjyX5ECr/QFwU5ItQAFHgE8AVNXBJHuBQwyuOLq9XTEEcBvwIHAhg6uFvGJIksZsziCoqm8Ds13v/+hbrLMb2D1LfQq4aj4DlCSdXX6zWJI6ZxBIUucMAknqnEEgSZ0zCCSpcwaBJHXOIJCkzhkEktQ5g0CSOmcQSFLnDAJJ6pxBIEmdMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjo3yj2LNyT5ZpJDSQ4m+VSrX5bkiSTfb8+XDq1zZ5LDSV5Mct1Q/eokz7Vl97R7F0uSxmiUI4JTwO9X1WbgGuD2JJuBO4D9VbUJ2N9e05btAK4EtgH3JlnVtnUfcCuDG9pvasslSWM0ZxBU1bGq+k5r/xh4AVgHbAf2tG57gBtaezvwcFW9XlUvAYeBrUnWAhdX1ZNVVcBDQ+tIksZkXp8RJNkIfBB4ClhTVcfaoleBNa29DnhlaLWjrbautWfWJUljNHIQJHk38BXg01X12vCy9hd+LdWgkuxKMpVk6uTJk0u1WUnSLEYKgiTnMwiBL1XVV1v5eDvdQ3s+0erTwIah1de32nRrz6y/SVXdX1WTVTU5MTEx6lwkSQswylVDAb4IvFBVnx1atA/Y2do7gUeG6juSXJDkCgYfCj/dTiO9luSats2bh9aRJI3JeSP0+RDwceC5JAda7Q+Au4G9SW4BXgZuBKiqg0n2AocYXHF0e1WdbuvdBjwIXAg81h6SpDGaMwiq6tvAma73v/YM6+wGds9SnwKums8AJUlnl98slqTOGQSS1DmDQJI6ZxBIUucMAknqnEEgSZ0zCCSpcwaBJHXOIJCkzhkEktQ5g0CSOmcQSFLnRvn1Uc3Txju+Prb3PnL3R8f23pJWJo8IJKlzBoEkdc4gkKTOGQSS1DmDQJI6N8rN6x9IciLJ80O1u5JMJznQHtcPLbszyeEkLya5bqh+dZLn2rJ72g3sJUljNsoRwYPAtlnqn6uqLe3xKECSzcAO4Mq2zr1JVrX+9wG3ApvaY7ZtSpKW2ZxBUFXfAv52xO1tBx6uqter6iXgMLA1yVrg4qp6sqoKeAi4YaGDliQtncV8RvDJJM+2U0eXtto64JWhPkdbbV1rz6xLksZsoUFwH/A+YAtwDPjMko0ISLIryVSSqZMnTy7lpiVJMywoCKrqeFWdrqqfAp8HtrZF08CGoa7rW226tWfWz7T9+6tqsqomJyYmFjJESdKIFhQE7Zz/Gz4GvHFF0T5gR5ILklzB4EPhp6vqGPBakmva1UI3A48sYtySpCUy54/OJfky8BFgdZKjwB8CH0myBSjgCPAJgKo6mGQvcAg4BdxeVafbpm5jcAXShcBj7SFJGrM5g6Cqbpql/MW36L8b2D1LfQq4al6jkySddX6zWJI6ZxBIUucMAknqnEEgSZ0zCCSpcwaBJHXOIJCkzhkEktQ5g0CSOmcQSFLnDAJJ6pxBIEmdMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjpnEEhS5+YMgiQPJDmR5Pmh2mVJnkjy/fZ86dCyO5McTvJikuuG6lcnea4tu6fdxF6SNGajHBE8CGybUbsD2F9Vm4D97TVJNgM7gCvbOvcmWdXWuQ+4FdjUHjO3KUkagzmDoKq+BfztjPJ2YE9r7wFuGKo/XFWvV9VLwGFga5K1wMVV9WRVFfDQ0DqSpDE6b4HrramqY639KrCmtdcBTw71O9pq/6+1Z9ZnlWQXsAvg8ssvX+AQ+7Txjq+P5X2P3P3RsbyvpMVb9IfF7S/8WoKxDG/z/qqarKrJiYmJpdy0JGmGhQbB8Xa6h/Z8otWngQ1D/da32nRrz6xLksZsoUGwD9jZ2juBR4bqO5JckOQKBh8KP91OI72W5Jp2tdDNQ+tIksZozs8IknwZ+AiwOslR4A+Bu4G9SW4BXgZuBKiqg0n2AoeAU8DtVXW6beo2BlcgXQg81h6SpDGbMwiq6qYzLLr2DP13A7tnqU8BV81rdJKks85vFktS5wwCSeqcQSBJnVvoF8qknzOuL7KBX2aTFssjAknqnEEgSZ0zCCSpcwaBJHXOIJCkznnVkFY8f3pbWhyPCCSpcwaBJHXOIJCkzvkZgbRAfpta5wqPCCSpcwaBJHXOIJCkzi3qM4IkR4AfA6eBU1U1meQy4M+BjcAR4Maq+mHrfydwS+v/e1X1+GLeX+qV353QUlqKI4J/VFVbqmqyvb4D2F9Vm4D97TVJNgM7gCuBbcC9SVYtwftLkhbhbJwa2g7sae09wA1D9Yer6vWqegk4DGw9C+8vSZqHxQZBAd9I8kySXa22pqqOtfarwJrWXge8MrTu0VaTJI3RYr9H8OGqmk7yXuCJJN8bXlhVlaTmu9EWKrsALr/88kUOUZL0VhZ1RFBV0+35BPA1Bqd6jidZC9CeT7Tu08CGodXXt9ps272/qiaranJiYmIxQ5QkzWHBQZDkoiTveaMN/AbwPLAP2Nm67QQeae19wI4kFyS5AtgEPL3Q95ckLY3FnBpaA3wtyRvb+U9V9V+S/BWwN8ktwMvAjQBVdTDJXuAQcAq4vapOL2r0kqRFW3AQVNVfAx+Ypf4D4NozrLMb2L3Q95QkLT1/dE7SyPwi27nJn5iQpM4ZBJLUOYNAkjpnEEhS5wwCSeqcQSBJnTMIJKlzBoEkdc4gkKTOGQSS1Dl/YkLS2964ftoC+vh5C48IJKlzBoEkdc4gkKTOGQSS1DmDQJI6ZxBIUueW/fLRJNuAPwVWAV+oqruXewySNKoe7sq2rEcESVYB/x74TWAzcFOSzcs5BknSz1vuU0NbgcNV9ddV9XfAw8D2ZR6DJGnIcgfBOuCVoddHW02SNCZvy5+YSLIL2NVe/iTJi+Mcz1mwGvibcQ/iLDhX5wXn7tyc19tU/njW8qjzmtfclzsIpoENQ6/Xt9rPqar7gfuXa1DLLclUVU2OexxL7VydF5y7c3NeK8vZmtdynxr6K2BTkiuSvBPYAexb5jFIkoYs6xFBVZ1K8s+AxxlcPvpAVR1czjFIkn7esn9GUFWPAo8u9/u+zZyrp73O1XnBuTs357WynJV5parOxnYlSSuEPzEhSZ0zCBYpyYYk30xyKMnBJJ86Q7+PJDnQ+vz3ofqRJM+1ZVND9cuSPJHk++350uWYz9D7n6153ZVkutUPJLl+OeYz9P6LndclSf4iyfeSvJDkH7b6St9fZ5rXWPdXG8OC55bkl4bGfiDJa0k+3Zat2H02x7zmv8+qysciHsBa4Jdb+z3A/wI2z+hzCXAIuLy9fu/QsiPA6lm2+2+AO1r7DuCPz5F53QX88xW8v/YA/7S13wlcco7srzPNa6z7aynmNtRnFfAq8A/OhX32FvOa9z7ziGCRqupYVX2ntX8MvMCbvy3928BXq+r/tH4nRtj0dgb/OWnPNyzNiEdzFuc1VouZV5JfAH4V+GKr/11V/aits2L31xzzGrsl/Ld4LfC/q+rl9nrF7rMZZs5r3gyCJZRkI/BB4KkZi94PXJrkvyV5JsnNQ8sK+Ear7xqqr6mqY639KrDmLA17Tks8L4BPJnk2yQPLfTg+bAHzugI4CfyHJP8zyReSXNSWreT99VbzgrfJ/oIF/1t8ww7gy0OvV/I+GzZzXjDffbbch3nn6gN4N/AM8E9mWfbvgCeBixh8Rfz7wPvbsnXt+b3Ad4Ffba9/NGMbPzxH5rWGwaHsO4DdDL5LsiLmBUwCp4Bfaf3+FPjXK31/zTGvt8X+Wsy/xbb8nQx+dmHNUG3F7rM55jXvfeYRwRJIcj7wFeBLVfXVWbocBR6vqv9bVX8DfAv4AEBVTbfnE8DXGPxCK8DxJGvb9tcCy37a5WzMq6qOV9Xpqvop8Hl+Nt9ls4h5HQWOVtUbf7X9BfDLrb2S99cZ5/V22F+wuH+LzW8C36mq40O1lbzP3vCmeS1knxkEi5QkDM6tvlBVnz1Dt0eADyc5L8nfA34FeCHJRUne07ZzEfAbwPNtnX3Aztbe2baxbM7WvN74j9d8jJ/Nd1ksZl5V9SrwSpJfav2uZfBBHqzg/fVW8xr3/mpjWPDchpbfxJtPn6zYfTa0/E3zWtA+G8eh0Ln0AD7M4Hz4s8CB9rge+F3gd4f6/QsG/7meBz7dau9jcNrku8BB4F8N9f/7wH4Gh4LfAC47R+b1Z8Bzbbv7gLUrZV6tvgWYauv/Z+DSlb6/5pjXWPfXEs3tIuAHwC/M2O5K32dnmte895nfLJakznlqSJI6ZxBIUucMAknqnEEgSZ0zCCSpcwaBJHXOIJCkzhkEktS5/w+d/4H16eprbgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb79e3c14a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(shapes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "src = butler.get('src', **dataId)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "star_idx = numpy.where(src.get('base_ClassificationExtendedness_value') == 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "star_shapes = []\n",
    "psf_shapes = []\n",
    "for x, y in zip(src.getX()[star_idx[0]], src.getY()[star_idx[0]]):\n",
    "    point = Point2D(x, y)\n",
    "    psf_shapes.append(psf.computeShape(point).getDeterminantRadius())\n",
    "for xx, yy in zip(src.get('base_SdssShape_psf_xx')[star_idx[0]], src.get('base_SdssShape_psf_yy')[star_idx[0]]):\n",
    "    star_shapes.append(math.sqrt(xx + yy))\n",
    "star_shapes = numpy.array(star_shapes)\n",
    "psf_shapes = numpy.array(psf_shapes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "good_star_idx = numpy.where(numpy.isfinite(star_shapes))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([197., 116.,  48.,  40.,  42.,  31.,  34.,  23.,  11.,  11.]),\n",
       " array([2.64647888, 2.64905415, 2.65162942, 2.65420469, 2.65677996,\n",
       "        2.65935523, 2.6619305 , 2.66450577, 2.66708104, 2.66965631,\n",
       "        2.67223158]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAEPhJREFUeJzt3X+sX3V9x/HnSxCSoZNir03DDy8khQSWWd0NblEJG5siLgL+wdotrm5slYQRydyWoskkW0jQiWaJU1NGhS3IjwkIiWwOiZOYzB+3rEJLQQqW0aa0V9TBnMG1vPfHPZ1f6r29t98fvbcfno/km+85n/Pr/elpXzn303PuSVUhSWrXKxa6AEnSaBn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYdvdAFACxdurTGx8cXugxJOqJs3Ljx+1U1Ntd6iyLox8fHmZycXOgyJOmIkuSp+azn0I0kNc6gl6TGzRn0SU5O8tUkjyTZkuQDXfsJSe5L8nj3vaRnm6uSbEvyWJJ3jLIDkqSDm88V/V7gg1V1JvCrwOVJzgTWAfdX1Qrg/m6ebtkq4CzgfODTSY4aRfGSpLnNGfRVtauqHuymnwe2AicCFwI3davdBFzUTV8I3FpVL1TV94BtwNnDLlySND+HNEafZBx4I/BNYFlV7eoWPQMs66ZPBJ7u2WxH1yZJWgDzDvokrwLuAK6squd6l9X0a6oO6VVVSdYmmUwyOTU1dSibSpIOwbyCPskrmQ75m6vqzq55d5Ll3fLlwJ6ufSdwcs/mJ3VtL1FV66tqoqomxsbmvN9fktSn+dx1E+AGYGtVfaJn0T3Amm56DXB3T/uqJMcmORVYAXxreCVLkg7FfJ6MfQvwXuDhJJu6tg8B1wK3J7kUeAq4BKCqtiS5HXiE6Tt2Lq+qfUOvvMf4ui+Ncvez2n7tuxbkuJJ0KOYM+qr6OpBZFp83yzbXANcMUJckaUh8MlaSGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaN593xm5IsifJ5p6225Js6j7b979iMMl4kp/0LPvsKIuXJM1tPu+MvRH4FPAP+xuq6nf2Tye5DvivnvWfqKqVwypQkjSY+bwz9oEk4zMtSxKmXwr+G8MtS5I0LIOO0b8N2F1Vj/e0ndoN23wtydsG3L8kaUDzGbo5mNXALT3zu4BTqurZJL8CfDHJWVX13IEbJlkLrAU45ZRTBixDkjSbvq/okxwNvAe4bX9bVb1QVc920xuBJ4DTZ9q+qtZX1URVTYyNjfVbhiRpDoMM3fwm8GhV7djfkGQsyVHd9GnACuDJwUqUJA1iPrdX3gL8O3BGkh1JLu0WreKlwzYA5wAPdbdbfgG4rKp+MMyCJUmHZj533ayepf19M7TdAdwxeFmSpGHxyVhJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklq3HxeJbghyZ4km3vark6yM8mm7nNBz7KrkmxL8liSd4yqcEnS/Mzniv5G4PwZ2j9ZVSu7z70ASc5k+l2yZ3XbfHr/y8IlSQtjzqCvqgeA+b7g+0Lg1qp6oaq+B2wDzh6gPknSgAYZo78iyUPd0M6Sru1E4OmedXZ0bZKkBdJv0H8GOA1YCewCrjvUHSRZm2QyyeTU1FSfZUiS5tJX0FfV7qraV1UvAtfzs+GZncDJPaue1LXNtI/1VTVRVRNjY2P9lCFJmoe+gj7J8p7Zi4H9d+TcA6xKcmySU4EVwLcGK1GSNIij51ohyS3AucDSJDuAjwDnJlkJFLAdeD9AVW1JcjvwCLAXuLyq9o2mdEnSfMwZ9FW1eobmGw6y/jXANYMUJUkaHp+MlaTGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuDmDPsmGJHuSbO5p+5skjyZ5KMldSY7v2seT/CTJpu7z2VEWL0ma23yu6G8Ezj+g7T7gl6rql4HvAlf1LHuiqlZ2n8uGU6YkqV9zBn1VPQD84IC2f62qvd3sN4CTRlCbJGkIhjFG/4fAP/fMn9oN23wtydtm2yjJ2iSTSSanpqaGUIYkaSYDBX2SDwN7gZu7pl3AKVW1EvhT4PNJfnGmbatqfVVNVNXE2NjYIGVIkg6i76BP8j7gt4Hfq6oCqKoXqurZbnoj8ARw+hDqlCT1qa+gT3I+8BfAu6vqf3rax5Ic1U2fBqwAnhxGoZKk/hw91wpJbgHOBZYm2QF8hOm7bI4F7ksC8I3uDptzgL9K8r/Ai8BlVfWDGXcsSTos5gz6qlo9Q/MNs6x7B3DHoEVJkobHJ2MlqXEGvSQ1zqCXpMYZ9JLUuDn/M1azG1/3pQU57vZr37Ugx5V0ZPKKXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNmzPok2xIsifJ5p62E5Lcl+Tx7ntJz7KrkmxL8liSd4yqcEnS/Mzniv5G4PwD2tYB91fVCuD+bp4kZwKrgLO6bT69/2XhkqSFMWfQV9UDwIEv+L4QuKmbvgm4qKf91qp6oaq+B2wDzh5SrZKkPvQ7Rr+sqnZ1088Ay7rpE4Gne9bb0bVJkhbIwP8ZW1UF1KFul2Rtkskkk1NTU4OWIUmaRb9BvzvJcoDue0/XvhM4uWe9k7q2n1NV66tqoqomxsbG+ixDkjSXfoP+HmBNN70GuLunfVWSY5OcCqwAvjVYiZKkQcz5ztgktwDnAkuT7AA+AlwL3J7kUuAp4BKAqtqS5HbgEWAvcHlV7RtR7ZKkeZgz6Ktq9SyLzptl/WuAawYpSpI0PD4ZK0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY2b8w1Ts0lyBnBbT9NpwF8CxwN/DEx17R+qqnv7rlCSNJC+g76qHgNWAiQ5CtgJ3AX8AfDJqvr4UCqUJA1kWEM35wFPVNVTQ9qfJGlIhhX0q4BbeuavSPJQkg1JlgzpGJKkPgwc9EmOAd4N/FPX9Bmmx+tXAruA62bZbm2SySSTU1NTM60iSRqCYVzRvxN4sKp2A1TV7qraV1UvAtcDZ8+0UVWtr6qJqpoYGxsbQhmSpJkMI+hX0zNsk2R5z7KLgc1DOIYkqU9933UDkOQ44LeA9/c0fyzJSqCA7QcskyQdZgMFfVX9GHjtAW3vHagiSdJQ+WSsJDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGDfrO2O3A88A+YG9VTSQ5AbgNGGf6nbGXVNUPBytTktSvYVzR/3pVrayqiW5+HXB/Va0A7u/mJUkLZBRDNxcCN3XTNwEXjeAYkqR5GjToC/hKko1J1nZty6pqVzf9DLBspg2TrE0ymWRyampqwDIkSbMZaIweeGtV7UzyOuC+JI/2LqyqSlIzbVhV64H1ABMTEzOuI0ka3EBBX1U7u+89Se4CzgZ2J1leVbuSLAf2DKFO9Rhf96UFO/b2a9+1YMeW1J++h26SHJfk1fungbcDm4F7gDXdamuAuwctUpLUv0Gu6JcBdyXZv5/PV9W/JPk2cHuSS4GngEsGL1OS1K++g76qngTeMEP7s8B5gxQlSRoen4yVpMYZ9JLUOINekho36H300mHhLaVS/7yil6TGeUWvQ7KQV9aS+uMVvSQ1zqCXpMYZ9JLUOINekhpn0EtS47zrRlqkfHZAw+IVvSQ1zqCXpMY5dCPNwYfEdKTzil6SGmfQS1LjBnln7MlJvprkkSRbknyga786yc4km7rPBcMrV5J0qAYZo98LfLCqHuxeEr4xyX3dsk9W1ccHL0+SNKhB3hm7C9jVTT+fZCtw4rAKkyQNx1DG6JOMA28Evtk1XZHkoSQbkiwZxjEkSf0ZOOiTvAq4A7iyqp4DPgOcBqxk+or/ulm2W5tkMsnk1NTUoGVIkmYxUNAneSXTIX9zVd0JUFW7q2pfVb0IXA+cPdO2VbW+qiaqamJsbGyQMiRJB9H3GH2SADcAW6vqEz3ty7vxe4CLgc2DlSjpcFuoh8T8HTujMchdN28B3gs8nGRT1/YhYHWSlUAB24H3D1ShJGkgg9x183UgMyy6t/9yJEnD5pOxktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOF8OLmnReDm+iP1w/H4fr+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWrcyII+yflJHkuyLcm6UR1HknRwIwn6JEcBfwe8EziT6ffInjmKY0mSDm5UV/RnA9uq6smq+ilwK3DhiI4lSTqIUQX9icDTPfM7ujZJ0mG2YL/rJslaYG03+99JHluoWg6TpcD3F7qIw8S+tunl1Fc4TP3NRwfa/PXzWWlUQb8TOLln/qSu7f9V1Xpg/YiOv+gkmayqiYWu43Cwr216OfUV2urvqIZuvg2sSHJqkmOAVcA9IzqWJOkgRnJFX1V7k/wJ8GXgKGBDVW0ZxbEkSQc3sjH6qroXuHdU+z8CvWyGqbCvrXo59RUa6m+qaqFrkCSNkL8CQZIaZ9D3IcnJSb6a5JEkW5J8YJb1zk2yqVvnaz3t25M83C2b7Gk/Icl9SR7vvpccjv4czAj7enWSnV37piQXHI7+HMwQ+np8ki8keTTJ1iS/1rW3eF5n62tT5zXJGT192ZTkuSRXdssW3XmdVVX5OcQPsBx4Uzf9auC7wJkHrHM88AhwSjf/up5l24GlM+z3Y8C6bnod8NGG+3o18GcL3b8h9/Um4I+66WOA4xs+r7P1tbnz2rPOUcAzwOsX63md7eMVfR+qaldVPdhNPw9s5eef/P1d4M6q+s9uvT3z2PWFTP8Dovu+aDgV92+EfV10BulrktcA5wA3dO0/raofdds0dV7n6OuiM8S/w+cBT1TVU938ojuvszHoB5RkHHgj8M0DFp0OLEnyb0k2Jvn9nmUFfKVrX9vTvqyqdnXTzwDLRlR2X4bcV4ArkjyUZMNi+7G3j76eCkwBn0vyH0n+Pslx3bLWzuvB+gptnddeq4BbeuYX9Xl9iYX+keJI/gCvAjYC75lh2aeAbwDHMf0o9ePA6d2yE7vv1wHfAc7p5n90wD5+uNB9HGFflzH9o/ArgGuYftZiwfvZb1+BCWAv8OZuvb8F/rrF8zpHX5s6rz3Lj2H61yEs62lbtOf1wI9X9H1K8krgDuDmqrpzhlV2AF+uqh9X1feBB4A3AFTVzu57D3AX07/tE2B3kuXd/pcDi2IIZBR9rardVbWvql4ErudnfwYLaoC+7gB2VNX+K8UvAG/qpls7r7P2tcHzut87gQerandP26I8rzMx6PuQJEyPT26tqk/MstrdwFuTHJ3kF4A3A1uTHJfk1d1+jgPeDmzutrkHWNNNr+n2saBG1df9/0A6F/OzP4MFM0hfq+oZ4OkkZ3Trncf0f+5BY+f1YH1t7bz2LF/NS4dtYBGe11kt9I8UR+IHeCvTY88PAZu6zwXAZcBlPev9OdP/ADYDV3ZtpzE9hPEdYAvw4Z71Xwvcz/SPjV8BTmi4r/8IPNzt9x5g+ZHc1659JTDZbf9FYEmL53WOvrZ4Xo8DngVec8B+F915ne3jk7GS1DiHbiSpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN+z9DWNWbAbU94gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb794c79198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(psf_shapes[good_star_idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([192., 120.,  52.,  39.,  42.,  31.,  37.,  17.,  12.,  11.]),\n",
       " array([3.74783312, 3.75138452, 3.75493593, 3.75848733, 3.76203873,\n",
       "        3.76559014, 3.76914154, 3.77269294, 3.77624434, 3.77979575,\n",
       "        3.78334715]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAEUVJREFUeJzt3X+s3XV9x/Hny4KoSKTYa9MB9UJSScC4Om+Ym5OxoRNlE9kf2GYzZWOrJGp02Q+LJopbSNCJZsmmpowGtrACE/mxyH4gosREwVtWoeXHKFBGu9JeQYduBm1574/76Xbs7u29ved77j2U5yM5Od/v+/vrfT+56avfH+fcVBWSJL1ooRuQJA0HA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWpmDIQkJya5I8n9SbYm+WCrH5fktiQPt/fFPdtcnGRbkoeSvG2QP4AkqRuZ6YNpSZYBy6rqniTHAJuAdwEXAE9X1WVJ1gGLq+rDSU4FNgKnAz8DfAV4TVXtG+DPIUnq0xEzrVBVu4BdbfoHSR4AjgfOBc5sq10NfA34cKtfW1XPAo8l2cZkOHxzumMsWbKkRkdH5/xDSNIL0aZNm75bVSNd7W/GQOiVZBR4PXAXsLSFBcCTwNI2fTzwrZ7NdrTatEZHRxkfHz+UViTpBS/J413ub9Y3lZO8HLgB+FBVPdO7rCavOx3SlyIlWZtkPMn4xMTEoWwqSRqAWQVCkiOZDINrqupLrby73V/Yf59hT6vvBE7s2fyEVvspVbW+qsaqamxkpLMzHknSHM3mKaMAVwIPVNVnehbdAqxp02uAm3vqq5IcleQkYAVwd3ctS5IGYTb3EN4EvAe4L8nmVvsIcBlwfZILgceB8wGqamuS64H7gb3A+3zCSJKG32yeMvoGkGkWnzXNNpcCl/bRlyRpnvlJZUkSYCBIkhoDQZIEGAiSpOaQPqk8rEbXfXlBjrv9snMW5LiSNAieIUiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAmYRSAk2ZBkT5ItPbXrkmxur+37/9ZyktEkP+pZ9oVBNi9J6s5svv76KuAvgb/ZX6iqd++fTnI58J896z9SVSu7alCSND9mDISqujPJ6FTLkgQ4H/jVbtuSJM23fu8hvBnYXVUP99ROapeLvp7kzX3uX5I0T/r9i2mrgY0987uA5VX1VJI3ADclOa2qnjlwwyRrgbUAy5cv77MNSVK/5nyGkOQI4DeB6/bXqurZqnqqTW8CHgFeM9X2VbW+qsaqamxkZGSubUiSOtLPJaO3AA9W1Y79hSQjSRa16ZOBFcCj/bUoSZoPs3nsdCPwTeCUJDuSXNgWreKnLxcBnAHc2x5D/SJwUVU93WXDkqTBmM1TRqunqV8wRe0G4Ib+25IkzTc/qSxJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSM5u/qbwhyZ4kW3pqlyTZmWRze72jZ9nFSbYleSjJ2wbVuCSpW7M5Q7gKOHuK+meramV73QqQ5FRgFXBa2+ZzSRZ11awkaXBmDISquhN4epb7Oxe4tqqerarHgG3A6X30J0maJ/3cQ/hAknvbJaXFrXY88ETPOjtaTZI05OYaCJ8HTgZWAruAyw91B0nWJhlPMj4xMTHHNiRJXZlTIFTV7qraV1XPAVfwf5eFdgIn9qx6QqtNtY/1VTVWVWMjIyNzaUOS1KE5BUKSZT2z5wH7n0C6BViV5KgkJwErgLv7a1GSNB+OmGmFJBuBM4ElSXYAHwfOTLISKGA78F6Aqtqa5HrgfmAv8L6q2jeY1iVJXZoxEKpq9RTlKw+y/qXApf00JUmaf35SWZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqRmxkBIsiHJniRbemp/nuTBJPcmuTHJsa0+muRHSTa31xcG2bwkqTuzOUO4Cjj7gNptwGur6nXAvwEX9yx7pKpWttdF3bQpSRq0GQOhqu4Enj6g9i9VtbfNfgs4YQC9SZLmURf3EH4X+Mee+ZPa5aKvJ3lzB/uXJM2DI/rZOMlHgb3ANa20C1heVU8leQNwU5LTquqZKbZdC6wFWL58eT9tSJI6MOczhCQXAL8O/FZVFUBVPVtVT7XpTcAjwGum2r6q1lfVWFWNjYyMzLUNSVJH5hQISc4G/gR4Z1X9d099JMmiNn0ysAJ4tItGJUmDNeMloyQbgTOBJUl2AB9n8qmio4DbkgB8qz1RdAbwp0l+AjwHXFRVT0+5Y0nSUJkxEKpq9RTlK6dZ9wbghn6bkiTNv75uKr/Qja778oIcd/tl5yzIcSUd3vzqCkkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEjCLQEiyIcmeJFt6ascluS3Jw+19cc+yi5NsS/JQkrcNqnFJUrdmc4ZwFXD2AbV1wO1VtQK4vc2T5FRgFXBa2+ZzSRZ11q0kaWBmDISquhN4+oDyucDVbfpq4F099Wur6tmqegzYBpzeUa+SpAGa6z2EpVW1q00/CSxt08cDT/Sst6PVJElDru+bylVVQB3qdknWJhlPMj4xMdFvG5KkPs01EHYnWQbQ3ve0+k7gxJ71Tmi1/6eq1lfVWFWNjYyMzLENSVJX5hoItwBr2vQa4Oae+qokRyU5CVgB3N1fi5Kk+XDETCsk2QicCSxJsgP4OHAZcH2SC4HHgfMBqmprkuuB+4G9wPuqat+AepckdWjGQKiq1dMsOmua9S8FLu2nKUnS/POTypIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1M/4JzekkOQW4rqd0MvAx4Fjg94GJVv9IVd065w4lSfNizoFQVQ8BKwGSLAJ2AjcCvwN8tqo+3UmHkqR50dUlo7OAR6rq8Y72J0maZ10FwipgY8/8B5Lcm2RDksUdHUOSNEB9B0KSFwPvBP6+lT7P5P2ElcAu4PJptlubZDzJ+MTExFSrSJLmURdnCG8H7qmq3QBVtbuq9lXVc8AVwOlTbVRV66tqrKrGRkZGOmhDktSPLgJhNT2Xi5Is61l2HrClg2NIkgZszk8ZASQ5Gngr8N6e8qeSrAQK2H7AMknSkOorEKrqv4BXHlB7T18dSZIWhJ9UliQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKnp928qbwd+AOwD9lbVWJLjgOuAUSb/pvL5VfW9/tqUJA1aF2cIv1JVK6tqrM2vA26vqhXA7W1ekjTk+jpDmMa5wJlt+mrga8CHB3CcF6zRdV9esGNvv+ycBTu2pMHq9wyhgK8k2ZRkbastrapdbfpJYGmfx5AkzYN+zxB+qap2JnkVcFuSB3sXVlUlqak2bAGyFmD58uV9tiFJ6ldfZwhVtbO97wFuBE4HdidZBtDe90yz7fqqGquqsZGRkX7akCR1YM6BkOToJMfsnwZ+DdgC3AKsaautAW7ut0lJ0uD1c8loKXBjkv37+buq+qck3wauT3Ih8Dhwfv9tSpIGbc6BUFWPAj87Rf0p4Kx+mpIkzT8/qSxJAgwESVJjIEiSAANBktQM4qsrpM75dR3S4HmGIEkCDARJUuMlIx2Shbx0I2mwPEOQJAGeIUhDyxvpmm+eIUiSAANBktQYCJIkwECQJDUGgiQJ8CkjaUZ+9kIvFJ4hSJIAA0GS1Mw5EJKcmOSOJPcn2Zrkg61+SZKdSTa31zu6a1eSNCj93EPYC/xhVd2T5BhgU5Lb2rLPVtWn+29PkjRf5hwIVbUL2NWmf5DkAeD4rhqTJM2vTu4hJBkFXg/c1UofSHJvkg1JFndxDEnSYPUdCEleDtwAfKiqngE+D5wMrGTyDOLyabZbm2Q8yfjExES/bUiS+tRXICQ5kskwuKaqvgRQVbural9VPQdcAZw+1bZVtb6qxqpqbGRkpJ82JEkd6OcpowBXAg9U1Wd66st6VjsP2DL39iRJ86Wfp4zeBLwHuC/J5lb7CLA6yUqggO3Ae/vqUJI0L/p5yugbQKZYdOvc25EkLRQ/qSxJAgwESVJjIEiSAL/+WtIUFuorv7dfds6CHFeTPEOQJAEGgiSpMRAkSYCBIElqDARJEuBTRpKGiE83LSwDQdIL3kIFEQxXGHnJSJIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRIwwEBIcnaSh5JsS7JuUMeRJHVjIIGQZBHwV8DbgVOB1UlOHcSxJEndGNQZwunAtqp6tKp+DFwLnDugY0mSOjCoQDgeeKJnfkerSZKG1IJ9l1GStcDaNvvDJA8tVC89lgDfXegmZsleB8Neu/d86RMWoNd8cs6bLgFe3V0ngwuEncCJPfMntNr/qqr1wPoBHX9OkoxX1dhC9zEb9joY9tq950uf8LzsdbTLfQ7qktG3gRVJTkryYmAVcMuAjiVJ6sBAzhCqam+S9wP/DCwCNlTV1kEcS5LUjYHdQ6iqW4FbB7X/ARmqS1gzsNfBsNfuPV/6hBd4r6mqrvcpSXoe8qsrJEnAYRoISV6S5O4k30myNcknpljnj5Nsbq8tSfYlOa4t257kvrZsvGeb45LcluTh9r54iHu9JMnOnu3eMQS9Hpvki0keTPJAkl9o9WEc1+l6HapxTXJKT31zkmeSfKhtM1TjOkOvQzWubdkftO22JNmY5CWt3um4DrDPQx/TqjrsXkCAl7fpI4G7gDceZP3fAL7aM78dWDLFep8C1rXpdcAnh7jXS4A/GrJxvRr4vTb9YuDYIR7X6XodunHtqS8CngRePazjepBeh2pcmfwg7WPAS9v89cAFgxjXAfZ5yGN6WJ4h1KQfttkj2+tgN0tWAxtnsetzmfyHgvb+rjk32Qyw187102uSVwBnAFe2ff24qr7f1huqcZ2h1851+DtwFvBIVT3e5odqXGfotXMd9HoE8NIkRwAvA/6j1Tsd1wH2OadmDssXk/8D2Qz8kIMkeBvAp4HjemqPtW03AWt76t/vmU7v/BD2egnwOHAvsAFYvJC9AiuBu4GrgH8F/ho4ehjHdYZeh2pcD1i2AXj/MP++HqTXoRtX4INtuwngmkGO64D6POQx7XvAh/0FHAvcAbx2muXvBv7hgNrx7f1VwHeAMw78RWjz3xviXpe2X7IXAZcy+VmQBesVGAP2Aj/f5v8C+LNhHNcZeh2qce2pv5jJr1xY2lMbqnGdodehGldgMfBVYITJ/7HfBPz2oMe14z4PeUwPy0tGvWryVP8O4OxpVlnFAae0VbWzve8BbmTy21sBdidZBtDe9wxrr1W1u6r2VdVzwBU9P8NC9boD2FFVd7X5LwI/16aHbVyn7XUIx3W/twP3VNXuntqwjeu0vQ7huL4FeKyqJqrqJ8CXgF9sywY2rl32OZcxPSwDIclIkmPb9EuBtwIPTrHeK4BfBm7uqR2d5Jj908CvAVva4luANW16Te92w9br/l/Y5ryen2FBeq2qJ4EnkpzSSmcB97fpoRrXg/U6bOPaY6pr9UM1rgfrdQjH9d+BNyZ5WZIw+TvwQFvW6bgOqs85jWmXp2XD8gJex+S133vbIHys1S8CLupZ7wLg2gO2PZnJSy/fAbYCH+1Z9krgduBh4CtMcW10iHr9W+C+tt9bgGUL2WurrwTG2/Y30a5pDtu4ztDrMI7r0cBTwCsOqA/juE7X6zCO6yeY/Id5S+vvqEGM6wD7POQx9ZPKkiTgML1kJEk6dAaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJAD+Bxclr/CCdDlcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb794c97d68>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(star_shapes[good_star_idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.00075, 1.00175)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD8CAYAAABpcuN4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzsvX2UXGd95/l57kvdeq9+b3W31FLLerFkWY6NjBEBwwAGgYOTMGGAGQgnO2dN9izMHnZZltk9G8LsmRwm4cQneDJJvIETAmfJwBkYSAx2bANxDMJgW5bQi9WWWlZLre5Sv9b7rbovz/5xq0pV3VXV3Vb7RdLzOUeWuurWvc8tyc/v/t6+PyGlRKFQKBSKdmiv9QIUCoVC8fpGGQqFQqFQdEQZCoVCoVB0RBkKhUKhUHREGQqFQqFQdEQZCoVCoVB0RBkKhUKhUHREGQqFQqFQdEQZCoVCoVB0xHitF7AR9PX1yW3btr3Wy1AoFIprimeffXZOStm/2nHXhaHYtm0bzzzzzGu9DIVCobimEEKcX8txKvSkUCgUio4oQ6FQKBSKjihDoVAoFIqOKEOhUCgUio4oQ6FQKBSKjihDoVAoFIqOKEOhUCgUio4oQ6FQKBSKjihDoVAoFIqOKEOhUCgUio4oQ6FQKBSKjihDoVAoFIqOKEOhUCgUio4oQ6FQKBSKjihDoVAoFIqOKEOhUCgUio4oQ6FQKBSKjihDoVAoFIqOKEOhUCgUio4oQ6FQKBSKjihDoVAoFIqOKEOhUCgUio4oQ6FQKBSKjihDoVAoFIqOKEOhUCgUio4oQ6FQKBSKjihDoVAoFIqOrGoohBBfFUJcFkIcb/O+EEJ8WQhxRghxTAhxR8N7h4QQp6vvfa7h9Q8KIU4IIXwhxIFl59svhDhcff9XQojw1dygQqFQKK6OtXgUfwMc6vD+e4Gd1V/3A38BIITQgT+vvr8X+IgQYm/1M8eBDwBPNp5ICGEA3wB+X0p5C/B2wFnbrSgUCoXilcBY7QAp5ZNCiG0dDvlN4G+llBL4uRCiSwgxBGwDzkgpJwCEEH9XPfaklPJU9bXl53o3cExKebR67fl13Y1CoVAoNpyNyFGMABcafr5Yfa3d653YBUghxKNCiOeEEJ/dgPUpFAqF4ipY1aN4lTGAtwB3AkXgCSHEs1LKJ5YfKIS4nyDUxejo6Ku6SIVCobiR2AiPYgrY0vDz5upr7V7vxEXgSSnlnJSyCPwAuKPVgVLKh6SUB6SUB/r7+1/24hUKhULRmY0wFN8Hfrda/fQmICOlnAZ+CewUQowJIULAh6vHduJR4FYhRLSa2H4bcHID1qhQKBSKl8mqoSchxDcJqo/6hBAXgc8DJoCU8i8JnvrfB5whCBf9XvU9VwjxSYLNXwe+KqU8UT3nbwMPAv3Aw0KI56WU75FSLgoh/pTAyEjgB1LKhzfwfhUKhUKxTkRQrHRtc+DAAfnMM8+81stQKBSKa4pqDvjAasepzmyFQqFQdEQZCoVCoVB0RBkKhUKhUHREGQqFQqFQdEQZCoVCoVB0RBkKhUKhUHREGQqFQqFQdEQZCoVCoVB0RBkKhUKhUHREGQqFQqFQdEQZCoVCoVB0RBkKhUKhUHREGQqFQqFQdEQZCoVCoVB0RBkKhUKhUHREGQqFQqFQdEQZCoVCoVB0RBkKhUKhUHREGQqFQqFQdEQZCoVCoVB0RBkKhUKhUHREGQqFQqFQdMR4rReguLY4NZ3hkeNpnp6YY3KxRMX1EUBIFySjIW7f0sXHDm5lz1DqtV6qQqHYIJShUKyJU9MZvnH4PP98Zp5CucJCwUUuO6bkuDzt+ExnbD57aLcyFgrFdYIyFIq2NHoPxy9lKVY8pGSFgaixUPQIGQ4XF4t8/vsnGe2JMtIV4dC+QWU0FIprGGUoFE08fGyKv/zJBOfmC5QdD1MT2K7EX+PnZ7IVZrIVoMBz5xcY6Ypw+OwcOwbiVDzZZDhqhmhqqaQMikLxOkYZCkWdh49N8e+/c5x82cWvug2O385/WB3Xh8mFEtOZEovFCof2DZEpOTz05Dnetaefx0/NkoqYGBr85PRlvntkirfu6OWjKsehULyuEFK+/I3g9cKBAwfkM88881ov45qk9lT/xMkZTkzn2oaVrgZNgGVojPXFiYcNIobgxdkCUdMgbGrkbYdUNISUEiEE2/vj3H/3mDIWCsUrjBDiWSnlgdWOUx7FDczDx6Z48ImzpLMlFkvuK3YdX4KUELd0lgoVXsjYOJ7HyEiY8wtFKp4kHjaJhnRyZZdUxOTrh8/TnwirsJRC8TpA9VHcgJyazvB/fecYn/nWUcYv515RI1EjbGoIIchXXDRNoOsaFU/iy6C0drFYoez6JMMmtuPy1Jl5MiWHoVS4Hq46NZ15xdepUChWogzFDcap6Qx//Mhp/uHYNCVXchUpiHVRKHvMZkuUKh6akHRFTMqujy4ESEmx4lF2fXYMxDg1naM7apKKmGhCkIoEf37kePrVWaxCoWhChZ5uEGq5iL8/OsVCwSFrv/JehABMHcKGjuP7TC3ZJCIGiXCI/kSYnQMxjl7I8NJcgZhlcPtoClPXWSw6/PqOnqZzJcIGU0ulV3zNCoViJcpQXMfUmuR+8KtLLJY8BEEPRFhv3wvxctAEGJrA0AVbuqJouiBhGZRdH19KLi3ZhDUdx5X0xiyKFY9NCQvPl+iaIBEx6Y2FmM1VuGU4wlt39GIaetM1crbLSFekfl+qrFahePVQhuI65dR0ho//9c+5XLjiOdSMg+1t3HVMDQxdI2LqxMM60zmb3miImYzNcCpMPGwy3BXmcraMpkl0TeNz793BT8/M89SZebqjJm/f3YdlGGRKDof2DQLw0JPngMCTyNku5+cLDKfC/I9/+0suzJfYNRhna1+snr/ojuo8cuIymZJDKmLy8YOjfOJtOzfuRhWKGxhlKK5DTk1n+Jd//hTFVzi6JADPB02TREM62ZLHYMLiwLZufvirNBcWSmzpBUPT6I6F2D0YZ1tfnHv3jzCeLvCOmwdIRcymcz5yPM2n79nFu/b087XDk6SzNoaA+UKFn1Y8HE8igcmFAm/b1c+uTUl+dXGRfziaxTI0BJK5XJk/eXScdNbmD95/6yv7JSgUNwCrJrOFEF8VQlwWQhxv874QQnxZCHFGCHFMCHFHw3uHhBCnq+99ruH1DwohTgghfCHEihpeIcSoECIvhPjMy72xG5X/8Pe/4t4/e2WNhCaCfziC4D9hQ8PxJK4vkRKOXcwSC2mUPZ/J+SKWobEpaXE6nefEpQwPPDbOyekMiXDwnDKXt/n5xDw/n5jjsZNpHj42xeOnZtk7lOTg9h6mszZLJZdK1UgAFB2ffxqfZTZnczqdQ0qJJ4P3Q7pAE4Jv/uKiqpRSKDaAtVQ9/Q1wqMP77wV2Vn/dD/wFgBBCB/68+v5e4CNCiL3VzxwHPgA82eacfwr8cA1rUzTwv3zzGb7608k1y220QwCWIQJDsIzaPxhNCxLVvoSM7bFYqOD7PhNzBU5PZym7PgNxE8+XTC2VODK5hO9L+uIhMiWHC/MlJucLzOVtnj2/hO14WLqGqQse/NFZPM8nFTGZmCtSdvz6tTURrE8AZdfnzGwB25EIUTVgQiCEIKQLHM9XlVIKxQawauhJSvmkEGJbh0N+E/hbGbR4/1wI0SWEGAK2AWeklBMAQoi/qx57Ukp5qvraipMJIX4LOAcU1nUnNzgPH5vie0c3ZlOUQNltne72qwcIQGiCkC6RUuBJScWpfkZAvuxguzq6FmzYW3ujCCF4/kKGO0a72DUY53Q6j2VohPTg30HZk9wxmuLJF2d55vwiZ2YLXM7Z9RJeCegCGpe2kK9g6gLHo8mwOX7QDa4qpRSKq2cj+ihGgAsNP1+svtbu9bYIIeLA/wF8YbWLCiHuF0I8I4R4ZnZ2dt2Lvt74038c39Dzaa3ciSoC0DXwZKDnBFdCQoLAy6h44HgSU9couz5hUyds6liGxpnZAlv7YmzujlDxfBxPYpk6d4x2IQQUyx6Fikvc0tFFs2cjhKBqV9A1QU88xEfeuBmQQWhKBr97vs+eoUS9UkqhULx8Xm/J7D8EHpBS5lt5G41IKR8CHoJA6+mVX9rrm7NzxQ07lyAI4SBbq8ZqIjAGSUsnU3YRQqAhg9f9anWVhIip4UtJzDJYKjrkKy5lx0MIwaaExS3DKW4ZTtUrlQB+PjFPKmqSLbmUXZ++eIis7VByfCTgekGYydLh1s0pvnDf3npp7Dd/cRHb9bEMjVtHUmzpiderqBQKxctnIwzFFLCl4efN1dfMNq934i7gd4QQfwx0Ab4QwpZS/ucNWOd1y3/4+19t6PkG4iYLxWAwkSAI93iymh8QoAMVH5bsQGXWQxI2NEBSqcaJgvCVjyEE2/uiHLuYxdQ1dE3g+XDkwhLv2TfI9v54UynsXL6MpevcNZZgvuCQ9Xy298WYzZdxXEnR8QibGm/e3su/e9fOupH4g/ffygcPjKr+CoXiFWAjDMX3gU9WcxB3ARkp5bQQYhbYKYQYIzAQHwb+dacTSSnfWvuzEOIPgbwyEqvz/z19ccPONRAzqHigiSsbfi0n4MnghXpBVYMfJ/BxPJqMiy4EIVPncrbCpqRFpuRQcn1Gu6Pc1B9jPF3g3v0j3H/3WH2D741bbEpYjPXHqXVB1DyOT9+zq369WtPdV556qckoLDcMqjlPobh6VjUUQohvAm8H+oQQF4HPE3gLSCn/EvgB8D7gDFAEfq/6niuE+CTwKMFD6FellCeq5/xt4EGgH3hYCPG8lPI9G3trNw62e3V1TgKIh3TCIZ2QoVMqVZAIdCTusuNqtqHWjR2zDLJFB9eHeDjoxhbA9v44t21JkSu5/Gxinv64xda+OKaAcwtFLowXMc/Os2swxr37R+qb96npDA89eY5Myak322VKDndu6+KBx8aZWiph6YJLGZutvbEm0cCaNHnNOJyczrRszlMS5grF+lhL1dNHVnlfAv9zm/d+QGBIlr/+XeC7q5z3D1dbm2JjkEDZ8+kyTUxdYzAZZrFQIVd2iWgCz/Nxq1LhjR/yfInr+cTDOkXHJxUNMZiw6IuHmCs4HJnMELd0EpbBXdt7SWdK/PTsPCE9qHSSEr74w9MA3Lt/pL7B58sOU0slkmGDW4ZT3Lmtqz7kaCgV5snxWfK2y1AqjCbMen6jVgr70JPnSEVMMkUHBIxfzhMPG/QnwvXjlKFQKNbO6y2ZrXiN0AVkbY8t3Sb5ssumVBg9XyZTcqj4QXmcoQscL7AWtVoDX0psxydmGQwmw0wvlTg7W6AnZtIdDdXFB8/PFzh5KUtIF+iawPVhuCuM40m+dniynqtIRUxu3pSsexKH9g3yyPE0nufz1IuzTGVsyo5HSIefT0jef1tQSFcTDXzkeLquNpsruySrXs6Z2QL9ibASF1QoXgbKUCgA2NwTJaRrzGQC4+D7EtMAp6oL5QMhIZHVPgZPgqGB68vqnwVDSYsz6Ry+hIWCA4CuaezfnCJuGTw3uYQuJJYuGEhaREMGvu+TztpNGzzQ5CWcuJTh5KUMs7kyuibQCMpvz88VGZ/JsqtqWEa6IkwtlRhKBZ5DMNvCwzI08lWD1SgueLWo/IfiRkEZCgUCKFU8ZstlihUvkOIAvGUyIJ4PpqEhHR+v+klT1+iPGYRDJtv64vxqKoPt+ORsh5mMR3/CYjpjk4yY7N6UINtQCguQK3sMJsOcuJQhW3LIlz3iYYMd/TF64xZTSyWytsti0cHQNQxNoGkSu+IjBRy9mGEwFSFTcvjQnZt55Hi6nvzeMRDj2fNL1YFIgehg7bhWrGfjr+VSauEwlf9QXM8oQ6FAFzCbK5MIm0RDsNRi4l1NANDQwTQEyZDB7x7cBsA/HLuE73t8+9kLzOcr1T4MsEydiKlzdraA6/vcOpJksVAGIGHp5MoehbLLb98+xHePTAMEoSLH47nJJXYNBDO253I2JcfH0EAQTMoz9KD5LmMHRuFDd26ub9C1ctuemMXuagd4suqtNB7XyHo3/poHdDlb4rGTGfJll7ChU664/JePrTqCWKG4plCGQhGUv/qSkKFRaKMmKAnyEgMJi/l8hf64RTpb4tR0jkuLRRwfYiGdhKWzUHCqZbKSqYwNElJhg/F0Ht8PRAAzJZ3eWIhbh5N87/kZpJSUq2Eiywi6uY9NLXFTf4JE2CRsBPkR2/UJGxqJsAkIbhqIN5XN7hlKNZXbbuuL8/tvv2nVp/xOoa9Wn51aKpG3K/zs7AIhXSNq6lQ8n5+8OMfDx6a4d39HEQKF4ppCjUK9DtA7N7GvCU9Sf3JvpPHUGjCXLxMO6fhS8uPTgXSKZeoIgk08eNoPPlVyg3xEXzxEoSoRPtIVZqQryt6hJP0Ji4FUBF8Gm78QAs+X9SS0JgRbeqLcMpykNxaql+RCIA9iaPDxg6NXf/MEG39NzbZGp8T3SFeEI5MZQrqGZWpomsDQAoPxtcOTG7ImheL1gjIU1wG/sX9jZCoq/hVl1hqNFbGeBN+X3LN3gJhlEA8Z3DKcRBI8gRsaLBYruNXubF8GXkbJ8QBBJKTjVtVkn5tc5OSlLBXXIxUJgRAkIyZRy+DdezexdziFZeokqmWtb9s9wKZkGF9CyfEJGYL/4S3bVjy5N/ZhNIaQVpMbH+mKkFs2HrZT4vvQvkHyZRdNBNLqblVmfSAZIp211/aFKxTXCCr0dB3wZx85ADxzVeqxpibwq3pNlRYT8GrlsZ6EpYJDxfOJWzpnZgvELYOy42PqWpBLEFCpWpiXFkqYGoRNA19K5nJ2IOVRPea5ySW290U5MpknY7s4nuT8XIEtPRFu39JFznbrYaBY2GRH2CQZNtg7nOKFmQKnpjNNoaH1hpBqHNo3uGKqXi3xXRspe+RCBonk9i1dfOzgVsb6YsxkbCqeT8jQ6E9YuL5kMBl+2X8PCsXrEWUorhP+7CMHiFnHeHpigXzZIWe7VFyfNmrhK9A10IWGlGDpkrLX/MEgICWJmAZHpzJs6Y5SqrjMZm0ipsblbLnuSbgEXklQPgtlDzThIQUIBJauUXB8YprEMjROz+TI2l4wQ1sEjXwXF0vcs3eAF2YCtfkzl/P1tewcjJOKmCzky3z++ycZ7Yky0hVh12CMx06m8aVPKhJix0CMvvjaeieW5zZGuiL16qgvPTrOubkCcUtHIHh6YoHpjM37b9vEN35+gZhlkLB05vIV5vMVuiIGDzw2zq7BQKZElc8qrnWUobiO+NjBrUxnbC4tlZgvOJjVnXo1YyEATRMMpcJczlUou25dIbYRKSEa0imUPXYMxPjJC7MUHI+t3RGGUhaTi0HIRQNSUYNoKKhgmis4lD1ZNx5Fx8XzBQXHo1h2uJSxsYxAhny4K0w0FJSyPvHCHF+4by+PHE+TzpUZTFjsHIzTFw8zm7MZT+dxfJ+7xnp4aS7Pd567SMTUiZgatuPx7Pkl3rC1C1PX19Q70Uor6oHHxpnLl0mEDcKmXv3CBAuFCsWK5HPv3c3XDk9yYaFIxfW5YzTFrVu66+u5fUuXkg9RXPMoQ3EdsWcoxWcP7eaR42n+/ugUhbJHvuwQFoKy4+H6zTmHGroGtuOTzlYYiJucW3CRknqZa81oeNX8gKnB8aksC8UKUsJ4Ohf0ODQ04+Vsl7Lj4/t+/TxSBt6FLiBuafgSZgsOZcevzr/QmC8E50xYOums3bR5N8qRn5ktgIC+uIUmBDPZMjHLQNcEZU9iGcGUu+NTWbb3x9v2TqzG1FKJsuuRDF/p/bAMjZztMrVU4tP37OLe/SM88Nh40/pq65mYKzCTK5O3XUxd8PXD5/mjD+x/WWtRKF4rlKG4zqhtrLWY+8RsnlLFZTpjky97WHqwWQMYgqrst8SXkCu75MuBvLgBaLoIhgDJK4OMShWHeDhEqeIhBIR0jZLjoVU7tGv4PhT9IGClV4+TQNnx0asd3XHLxNAFpYqLqetEQzqeJ4MGvbDRFOs/tG+QP37kNAuFChXXZyZjk4wYvGFrFwDpbImK61Oo+GzujuD5Es+XIORVPcWPdEV4MZ2rD1+CQD49ZGhNXkqtI3w2Z3NmtsDpmSymLihVPHYMJIhbOmXH46kz8yvyKgrF6x1V9XSdUou57xtOYruS4a4IqbCOJJgYZwCGHgwWcmXgaWiAWf0X4QKVhp3flxDSBX1xi1+/qZfuWIiEZeL6El0LDErjVLxGzyXwRoKN26iOLS1WfHwZzLHojobwpaTiBp6F50vm85UVpa9aVWBKIjGNYK6FlDCbs8nZHsWKTzQUdG+7vmTnQJx37910VZvyoX2DuJ7PCzM5jkwucnxqienFIj2xUNNQpJGuCOfnCjw3uUTZ8YiFdHK2i+tD3g5EDicXS5Qdj68fPv+y16NQvBYoj+I6Zs9Qiv/4gf0rZLc9GYSDJhdK1Nom9OpQosoyxXINQEDU1PjN20f4p9PBPOtCxSWka5QdD9MQ+DIoE9VFUDVUcnwqno+UYOoQr06586shrd54iOGuCGcv59naF+WOnihHLwYdzrGQwWhPpKn09ZHjabb0RNk3Emz6c3mbw2cXOHEpi2VoJCMG8/kKCcusN+yNp/O8Z99gXZ785SSUJ2bzzBcqRAwN2w16QQr4vGl7d9N5Du0b5NN/NwMiCE0lwiaXcxXCBlzK2MRCBgJIRQzlVSiuOZShuAFojPOfms7w9cPneerMPKmIiVso48sr40uXIwEhIV/x+d7zl7AdD11Ab9yi4sp66KqmJhuqdla7vkQTGp7vo2vBU75R9TxMQ2MgYQFgGBoV12fnYJKdg0mgORdRW/M/npwBCYmIyY7+GP2JMHdt7+bohSxLpSDRvWdTdSqe7QQNe5rg8VOzeJ7PTNbm+cklHj0+w6feedOaO6e/dniSrmiIVO+V9dQS7Z942876a3uGUmzpjZApOuTKLl3RENt6JZeWSni+xDI1emIhNCGIhFYv11WCg4rXE8pQ3GDsGUrxR1Uv4+uHz/Pfnr2IrHoDy8eU15rvak6G63kIgj4LX0okgexH2fEZSISqqrNQdiVCBhPuNF2jL27heIHQoKHBSFcYicAydd401s3J6XzToKLz8wWGU2E+8+2j9SFFIV2ry3w8N7nEHaNdWIbBPXuD8E/NuDROxTs5ncXzfMYv57EMjZ6YSdZ2efBHZ9neH1/TxpvO2gzEQ02v1RLty9k71Dz/ey5v870jl+iKGox0RSi7PmXX5/bRVMdyXSU4qHi9oQzFDUrNYGztjfDlJ85QXB5zAkINiW9dQNmt9ltogoW8Q38iRMjQCJuS7piFEIJi2UMi0XTBQMJi50CcuYJD3nbxfEkyYvC+W4fr18iUHDQEJ6ezpLM2CcsgGtIxDZ2euFEfUrR7U5yJuSKWIbB0wYlLzdVMrZrlkmGDmayNZWj1RHQybDBfqKx5eNFgMlxVvL2Szqsp3i5nedOeqev0JywMXZAru+hCYOqCX55boiceaht+Wt40WHE9Jmbz/G/fOsY9eweVd6F41VHJ7BucT7xtJ3/8O/vpj4fQNIFlBPMewoaouxi6uCLrIauT7Sq+ZKHoMJy0ePfeQXYNJrhn7yY+enArv3X7Zg5s7SEVMZlassmVHExdsLU3yubuaOB5SEmm5HB+vkDZk+wdSnLfbcPommA+X6HiemhC1DvA5woOb9jaRdjUKXs+jnelmqmWuE9FTKYzNqmIyf13j3HLcIqFgoNlXPlnXnZ9emOhNQ8v+vjBUQplt+ot+WRKDoWy21JjqtU6PvOeXewb6WJnfxzXkzieRNdgKGm1lBaphdmenpjn8MQ84zNZnptcAinxpb9mSRKFYiNRHoWCe/ePNE2Ye+bcPDPZMhXfR8MjZOi4vk/F8XGqeYxAdlxydCpL2DQwdMHT5+YplD1ils5QIsRC0SUZMZHV5EcibPKBO4abupWHU2FMQ7/y9NwgDdKfCJMMm5QqLnnbpS8epi8e5txsnplcma889VJT/L7VU/ajx2fI2lcm3ZVdn2290TUPL6rlMr52eLLu8dw6nOS/P3+Jrx2erI9rbVzD8nVs74/z+e+fxPF9+uJWvWM8U3KaPJtayKkxzPaLl/L0xEwwDFKR0JokSVR+Q7HRKEOhAJolLIa6I7g+7BqM43ge//ziPBVP4hO4oD6BsxEksQVPn1uoq8bGQjoVx+dXl3J0RUzu+7UrSeNMyWE8XWiSBf/Mt4/SE7/yz7DRMADsGIhx+OwCiXCgFXV+rsCRC0GOoha//9Kj4wwmLSqeXGE4PvXOm3jwR2eZL1TojYXY1htF0zQO7Rvk4WNTdQMwmAzz8YOjLZPcNUP6jcPn+ecz85ydK2A7HpGQQbboEDV1Hnqy2DaHsGcoxWhPlLvGeuolvrBSnbYWcto3kuTZ80tYhsD3A89L1zT2jSRbfq6RU9OZpn6TF9M5jl1c4rOHditjoXjZqNCTos6eoRSfvmcX/+/v3skDH76Nsf44sXCIu3f10RUx8aoJ75AuiJg6pi7qfRhUBQUL1UY8KSVFp1ldsNUGt1y1dcdAjHzZq5/b1HXG+mLcMpxkOmMzkytzx2gX2/riaELgeB7n5gqcuJRtqRZ77/4RHvjQbfzG/mFGe2Ns64tz/91jTMzm+eIPT5MtOQzEQ2RLDl/84WkePja14nupPekfv5SlK2KwVKyQsV0MTWCZGjO5MqmIySPH24syrkWdtiZ13hcP18NsmhZIqr9haxd98XDLzzXy9cPnmZwv1r9vgMn5ourdUFwVyqNQtGR5COXUdIaP/vUvyJcddCEAWS2BDRrqXB+kKwFJznbQhFihFdVqg2uVAB7rizGYDEaojnRF+NB7dtXX8plvH63PxAY4czkQ63M8iSZEy9BMq3DQ579/kphlNKjMBs9Mf/mTiRVCfrUn/Yrnk7AMvGrzYTobKOEWKh5ISbKhpHc5ndRpa4x0RepVU7Uw26akxemZPKau1z2q8XSeLb0RHnhsfEVY6ciFJeKWXk/eh00dpOTIhaW2a1MoVkMZCsWaCEJT2/iTR18ETWIKUe+7CMxGEIqq/Z4Maxi63lT22mpedUvV1gbDsJzGzRQgazuENEG8YejQWtRiW5W96kJOls6kAAAgAElEQVQyMVfkDctmWeRshz1DSZJhE7s6ha9Ydik6wcAlT/qcTucwNK3tdLt26rTLm/aWGxNN0/jUO29iPF3gxKUMFxdL7N4UZ7S3tdCgQKxoh5HV11XuQvFyUYZCsWY+8badHJ3M8OzkErbrERaiHv4pOYF3IQg8jELF53cPbiIRNltujC9301q+mYZ0jbztcsvIlc92Cs3UaFX2ms5WiDd5GcHvU0slcrbLjoEYz55fIm4ZzOXLeL6k5Egipo4uBMmIUe/RAFre3/L7X56Qb2VMJmbzHJ6Y5/RMjnC1QbGd93T7lhSHJ4KcUa1DPV/2uHlTXPVmKF42Qso1Dix4HXPgwAH5zDPPvNbLuCFobAYruy5/f3SGihvkFCpuoOdkGYJUNMSH7hxtSly3Okejt7HWTavRyNQa8rb2xtZ1roePTfHFH56uz5LIlT0mF4q8dUcvNzd8zpeSF2ayxC2zfs8nL+UYT+cAiaFpJCMmPbEQEVNnvlDhTWM9FB2/7f2t5/4b13k5W0ITgpLjM5i0CJsGcUsnGTH564/fWf9uvvRoII1edj0sQ6cvbjGYtAg1VJfBlSbFVn9HihsDIcSzUsoDqx2nPArFumh+6nUZ643y0nyBkKHRFdWJmDo528HxfB47mW45vOflTqFrXMPy/EmnkE4rlpe9DibDvG1nH4OpZk8kZ7vsrarx1u757bsHuGU4yS/OLdITMxHVSibb8eiNhThyIcMbx3ra3t967v9rhyfruZSc7VCqeDieTzprc8twiqztkrXdevPenqEUn3nPrhXezFeeeone+NpngisUjShDoVg3y0Mon/6vR4FAnuPiog0C+mMmnu/zxR+ebhre8/nvneB0OgdAzDK4bSTFrk3J+qa11pLVdutZjeVG5Qv37V3RxwArE86tjNOPX0gztVTEl8G9J8MGbxzrZfxyvl5xVKNxU65Jkrd7v5HGXEpPLMSZQh5BdXKgG3TT7x6Mt0zeN4a3JheKVByPeNjgzGyhPh/jluHkmr43xY2NMhSKq2LPUIpPveMmHnziLFNLNiEjiJ3XhABjlmAmV2asP87lbCnQX/IlcSvot/jp2XkABlMRCrZTD7M0lqwCLY3FevMcjSEfQ4OfnL7Md49M8dYdvXz04NY1JZxrTMzmqXhB4h4knhRkSsEwo+WzvqE5b7I8Ib/8/cb1Op7P6XSeeNigJxrCqqrYCoKKpn0jSXpi1gojs1wvynE9fn5uAVMIeuMhDC24ZjpbVkq2ilVRhkJx1dQa0v7Xbx1tUng9cmGJhKXXm+eOXsxgGVp1UFKgJhtC8NzkEm/eYXBxsdSyZPVrhydXGIq1Cuc1GpPJhSKbEhYV1+P5C8FawobgsVNp/vnMPG/Z0cvHDm5titmfms60lCn/2uFJ+hMWw11hFgqV4OleBoOZPnZwa8dS2LWUyj58bIoHnziL5/nYjof0JXbFQ0owNI237uhj16YrarvLjczy8Na2vjgvTOdYKjlUfEkybHLr5hSmrq855Ke4cVENd4oNYc9Qinfv3cRd23s5uL23Lr+RK3v10tV82UUXELUMhrvCGNWmOtv1uf/uMbJll4SlUyi7TMzm+dVUhvPzBU5MLa3QNmrcCGsVQMub3mrGJFMtd53PlxlP5+sGy5eS+UIFx5N0RwxOXMo2Nest/3zNGD18bIrTM1lmszbzhQrd0RA39ccZ64uSLbtttacaQ0Pv2tPPyeks3z96iZPTWd61p78pBPbgj86CgK29MYZSYdzq92QZGjdvSjCYitT1sjIlp2mIElxp3muk7HpBN321hFZKladQrA3lUSg2jOVPypuSFpeWSuwaiFen2emUKh6jyRDRkEE0ZAQqrxGTPUMpBpNhLmdKLJZcbCcwKp4PvoDPf+8EOwbidZmOE5cy7Blqjq+3k8SoPVX3xS0yJYd01mZ7X6x6rCAS0rBMnXzZqxubPUMpvn74PBOzeRxPEg8b7OiP4fs+Dz5xlrChYzsepWKF+XyFZNggbl0Z39opb3JqOsPjp2bZO5TkrrEecrbL46dm69LnjxxP43iBeKEQgk2pCF3VsNNob4x/+5Ztq4bHloe3ZnM2S7aLIUR9LOtzk0vsGogzVi3nVSjaoQyFYsNYHuPf1hfn3bcM1que3jDaxbMXlnB9ie/75MoehbLLp95xExAotf777xyn7HromkBWpUF6oibj6Txl1+fuXf1kSg4XF0tETb1pk2slidGYNN4xEOOZlxbx/EBwr1gJFGp7oiHKrk88bNSNzanpDE+dmac7YjRtrLqgOmY1xtPnFtG0YJ541nYoVFw+/MbmhsJWrFb1NLVUojcWaprTbRkaCwWHgzdF1pS8bzTatuPy4xdmqbg+UhcsFgMvqDYF8H/6Fzet/pe7RlRT3/WJMhSKDaXVJnZvw5+XVzV96h3BtLlT0xnG0wUsQ6NQdvE8STikM9IVoVhxEUI2yXTsHoxzeiZPT9xakyQGQF88zFDSYqFQ4cxsAYmkNxpC1wRl1+eW4WTd2DxyPE3Y0JjJlvFkMKApHjKYyZUDsUSfen6i4gUNcG/c2k2xsnpf0mpVTyNdERzX43Q6DwRGIlvVlloeYur093D/3WN8/fB5fnZ2gYrnM9YbiD3O5Ss4ns9AIkwqam7YRq4GLl2/KEOheFW5d/9Ix8T0WF8MKcGVks1dEWKWwXy6QkhvlukY7Y1RdLzqzIvg6fXObV1N3c67BmM8fmoWCDbi83MFzs4VefNNvUQtnedeWuLCUolwSOcNW7sJGXp9ut6PT8+SKVXq4151DbLFCiDYlAzz4myewWSYTalIIOth6ty6pXtN8f6RrkhdKj1vu8TDBpsSVt07CryBIrsH40xngjyIqQdSHusdn9qfCPOOmwc4NZ3Fdjzi4UAHyjJ19g4lmyqvrpb19sco7+PaYdVkthDiq0KIy0KI423eF0KILwshzgghjgkh7mh475AQ4nT1vc81vP5BIcQJIYQvhDjQ8Po9QohnhRC/qv7+jqu9QcXrn8YNZudgnETYwPMkl3M2dsVFAqamsaM/Vv9MrRHu0/fs4ksfvI1D+wZ5/NRsU+L58VOzvGtPfz2p3Kg8O5CIcOjWIe7ZM0B3zML1qQ9LMg0dx/OpuD5CgBASz4eKJ+mOmui6RkjXsB0P2/Eouz47+mNrkg4B2DUYVIRlSw6xkEa25HDkwhK7BoP7q3kD2/oCTaff2D/MAx+6rWM/SbvE+8npDImwwY6BGGU3qKAK6YKFfKVlEvxqaJVAb5csb7deNZDp9claPIq/Af4z8Ldt3n8vsLP66y7gL4C7hBA68OfAPcBF4JdCiO9LKU8Cx4EPAH+17FxzwPullJeEEPuAR4HO3VaKa57GUExfPMybd/RydHKJlxaKTGVs4iEdQxfkbZfeuGwZZqorvLoeT5/L1hvKfmpq/NEH9gMrlWch8ExMQ+dLH7yNBx4br8tcVDwfTRNoIpj9HQsblF2Jrgvuv3usPpuiO2rya1tShAx9xZraPTGPpwvcMdrFTLZM1nbqobTxdKEepmvV4NeqTHf5/bfTqapJl5+5XGAuX6Y3bm14SGit/SGd1qtKdV+frGoopJRPCiG2dTjkN4G/lYFo1M+FEF1CiCFgG3BGSjkBIIT4u+qxJ6WUp6qvLb/WkYYfTwARIYQlpSyv+Y4U1xytcgk3DcRxfXjTTb0kwgaT8wVOp/MUHY9bhlMrBAb/8eQMtuORLbn0xMwgWet4PHVmvt5QttpG1miwTF1DAxzfx/HA0DWSYQ3TMNgzlOI/fmB/kyFIRcwVa2oXr59aKtVnY9TwpWz75N1olPYMJVrG/tvlPZLhoLIMoCdmsWdIX5eu1npYS39IjfV0pyteezYiRzECXGj4+WL1tVav37WO8/5L4DllJK5/Wm0wp9N5dg3GmxrGumPWChG7+kS3aihFCHA8H8vQ0TVBd/RKuetqG1mjIRlKhrmwWCJsGCTDGn1xi5zt8sYtVzbXTtVHnZ6Y19OZ/dCT55iYzdMVCf5Xff5ChjtGu5rKeJevvfGctTGt69XCejmsp7N9Pd6H4rXndZnMFkLcAvwn4N0djrkfuB9gdHTloHvFtUOrDWZzd4StfbGm4xpLVx984kV+8dIiS0UH15dYOjge6HrwdH5pschAKsLto6n6U+pqG1mjIdm3OclisUK+4hHXA893rC/GRw9uXdM9dXpi/rdv2bamJ+/lA5NqHviZ2QJ3jfU0PX13MoLr0cK6WtZ6rfV4H4rXno0wFFPAloafN1dfM9u83hEhxGbgu8DvSinPtjtOSvkQ8BAEMuPrX7bi9cTyDeaBx8ZbPnGGdMHnv3eCk5ey6FrQ0wBQrk5dlR7omqTk+rxhaxemrjOQMNteZ/kaGpVx94+kmFwskSu7WKbOB+4YXvOG2+mJea1P3jVjUxuYFDZ1LCOYv7H86Xs9T/OvB6619d7obISh+D7wyWoO4i4gI6WcFkLMAjuFEGMEBuLDwL/udCIhRBfwMPA5KeVPN2BtimuUdk+cUVPjpfkClqlTqjTPoBYEM70FgTKtqa9MMNdYnmiOhgRPvDBX7+945819vDBT4MC2nvr1G7unV2PXYIwHnziL60t6YiabkmF0XauvpZXC6/Ikdc3Y1AYmASAlpi7aTgu8ljbaa229NzKrGgohxDeBtwN9QoiLwOcJvAWklH8J/AB4H3AGKAK/V33PFUJ8kqBySQe+KqU8UT3nbwMPAv3Aw0KI56WU7wE+CewA/kAI8QfVJbxbSnl5Y25Xca3Q7onzK0+9FPQDWAaZkl8fw1rD0AIJ7kTYWJFgrrE80fzUi5c5djGLoYMmBHO5Ms9PLnLrcIK37ArKR9dTlVOT6Ni96UofRNYOOtCXVzJ1alCrGctUxOT20RTPvLRIOlumLx4iampN52lVXfVyJNsbudrPr4bqo7h2UBPuFNcUDzw2zjd/cR7Ph0ypgpTgNczutgyBoWt8+/cPrqoiO9YfZy5v853npnA9H19CyAg24IrrEzI0fvv2EfoTQa7Bl5LpjM2XPnhbxzX+n985xolL2SaNqFrZbWMivlV4rdaIN9oTrTcNNs3LHgx6K2oe1rv29PP4qdmmaXknL2W4tFRiOmMTMjSGkhYh06BQdvnce3e3bHhcvmFPzOZXTABs9/lWrGYErnbKoWJjUBPuFNclh/YNcvjsHCcvZYPyVUndqwgbAss0ePNNPU1hnafPzfNiOk/E1BnuCvPSfJ6zlwN12orn43iy7pYYWqBoK6qaTmdmC0CQQF7IV+iJhzrOb2inEfVrW1Lky82hsuUJ79mczXg6j+P73DXWU28avP/uMWZzNrO5Mi9eLjCdLbOjP0YqYvK1w5NNHda1mR+2E4ynBbi4ZDPaEyVmGSsk29t5NS+mc2uWfG/1HdTOaerB3I//+swkvTGLoVSYW4ZTzOZs1UdxDaEMheKaYs9Qii/85i08+MSLHD67QK7sYmiS7miIgWSEvrjFp965M5jn8KOzFMoOs7kyphZsdOmszVLRJaRrlA2NkuMBEldS31iFEAiCRrvppRLZogMCdA2GklZH/aJHjqfpjpr189RE/Z49v0h3zOIz3z5af8JenvA+M1sAEajc1jStFvJlPvvtY0zMFYhbwfzrRuOTztrcNdZTv35NQr1YCYQVtWql1EzWZs+mBOmsvWK9rTbsC4tFdg00q8omLH3F51tRO6fjeRyZzCBlMEtjxrXxPEnU1Hn+YoY339TT5E1tZB+FCmttLMpQKK459gyl+C8fDbzlVhsCwINPBPMcPD+YQuf6HmXPI1Pdh0qOh1udtGdqIqiekhKn6mFIoDdqUHJcQkbQR7FjIEZfPHjqbvfkO7VUYu9wgiOTgRSFZWgUyw5TmTJ7h5JNT+21sBEEm+RCvoKuBSq3cMXDSOdsYiEdz4eZbJmuiEG+7PHoiTSRkM75uQJj/XFmczZz+TICQDb0kwiB40lyZa8ug9643lZlvKaukSt7dU8CaPn5VtTO+YtzWSxDYy5fJqRr+IBlaszkynRHTU5N5xhMXqncWq2PYq2bvxIn3HiUoVBc07SqnHngsfF6tdFMplQf0uNVf+nVWJWuge36mLqGaYBd8Sm7QdipP24SDZkUKy53butu2tA6PfnWvIQ3bO3i6IUME3PBfOpYSCdmGXVPAWA8XWhK2PfEQwwlLfriwWZc8zB0TTCQtJjOlHE9n0sZm2hIx/V8dg2kOHJhiUypQjpbQRMC1/MJmxolxwcC/SpNiCZJ9+XrXV7G+8at3fzqUja434YcxfLPd/oOsrZDwjIouz6aAEvX6uW9t21J8tMzC2RKzpr6KNaz+St5kI1HTbhTXHdMLZXoiZmUXR/XD9RfITASgsBoCA3G+uLc1B9nqCvCzUMpRnoi9CdCbEpaGLoOSKQMwkaNdHryPbRvkEzJIVdycT1JX9zC1IM51c9NLjGbC0I3NWOzp0HY8Av37UXTNDIlB19KFvIVpJRBaa2mMdwVxvN9PF9iaBrbemPctqWbO0a7eHG2gOP7jHRFCBk6kZBBwtLwJTh+sHm3SkTX1lu7Zu3Pn3rXTj733t0kIyaX8xWSEXPNiezaOWvCibqoCSpemfthGQZv3dHbdgrgctYy0bDx73+t4oSKtaE8CsV1x0hXhIrjMX45H/RWEHgPvhd4E54MktYLhQolJxgP+u/euYMfn55jaqHA8xczeH4wgyIW0phasjk3m2drX6ztk29jWCRqajyfzuH4fmAotBiadqWruj8RbmlslpcE1zyMeNjg2fNLWIaGqetELUF3LMT+qpzIaG+M5y9meN+tQ2hC8GI6y/MXMzieoCsa4gv37W27wbcrQ4bA4xntiXJwe++6Yvy1c37j8Hl+dHoW15OUKy4zGZt42GB/T3TdFU7r0YZq5SUdu7DIi7MF3vYnP35FSn2vd5ShUFx31OY57BqIk62GQDxfYuoQt0yEkGRKLiXHw9QEA3GLx0/NcjlT4pnJJTQRJLY932eh6NMbNZnJlQmZessO4uVhkVy1c/rNN/UwmIwwl7d59vwSIV2Qqz6xtwuzNIbSHj42VW/aC5uCYsXFQxI3dExdcGQyQzxcYFPCYjAZXDeQak/SFQ3VS3TH0wW2d6jUaqVU+3//9+NcWChhux5hQ+enZ2b5f35r37qMxUcPbuVSxmahUCFrOxTLHiUnmKH+sYNb1xUGWo8+1lzObhJRnF4q8cvziwwmw/TFQmRLDl/84WkAZSzWiAo9Ka47ak+0Y/1xbhlJsrkryj17Brnv14aJWQa2I9ncHWG0J0pP3OLAWDepiMmLl/PIaqhKCIEmNASBBzLaE+VLH7yNT9+zq2NMvBYWqSVrgbrEt6i2jVdcj6ip8ZWnXuKBx8ZbzmBobNrriZmUHB9d0/idO4aBIDndOMvinTf31Q1QOlvi6YkFctVcwHpnPXz58Rd5YSaH5wcVSp4veWEmx5cff3Fdfw+PHE+ztTfG23cPcN9tI3z4jaO8d98Q/YnwunMF7UJkjfM0agbbNHR+fUdQCfazswucnM4xmAwzkAijaRqpiFkvFVasDeVRKK5LGp+SG8NCB7f3cHhiHl0LSlf3jSTpi4fxpaTgeMQtnULFw5dBEjlqathOEPtvV3XTKiyydzjRlKw1dZ3t/fGmBrmeuNE2KdtofGKWAdU+jkdOXGZXf4ySK5tmWRQrsh5CevrcPPGwUb+3xnM29pe0qx765flFwoaGVe3+trQg+//LZbma1dhIKfG1aEM1J7FNBnYHXsh3j1ykLxZqXscaS30VAcpQKK571io4GDV1ihUPTdMwq8mNQsUjGjLYNRhrW3XTKizSmKxt3NjWWpFT22RnczbPTQb5ie6oQXqmxEy2wh1bu5o6xmuJ8UbDpTXMe2lU3l2tesjxJCG9+TsU1YT0ethoKfHVtKHWUupbKLssFCv1EbSdmicVV1CGQnHD0U5wcPdAnOenMkRMrd5PAYJbR5KMpwstN/hvHD4P0BQTD5tG22TtV556CVOHn09kydoOybDJ9v4oU0vNXdu1udrPTC5SqnhEQzrRkE4iYoK4khSfzdmcuJQlZzv8q786TDJskLVdKo5Xn8ENVzbotRiq0e4IL80X0UQQhvN8ie1ItvVGN+R7fqWkxFcr9S07XjV0BUjJrv6Y6q9YIypHobjhqIUxlpdmbuqO8pYdvcQsA0PX6E9YvHNPP5u6oi1LLsuuyz+fmV8RE6+4XtvNx9IFT08sYjseCcvAdjyenljE0punPdbmaudtl7AhKDs+6WyZ7b1RZLV0tpaLmMuX8X1JtuRwcaFE1BAcubDEudn8inj+WkpHf//t24lZetVAuHi+JGbp/P7bt2/I9/xKbcqrlfqWq8Y/Zum8ZWcft432tC2xVTSjPArFDUmrMEbwRGqwe9OV1xufUJc/rZ68lKM7ujImnoqYbTfDwEe5QqnisVCs8LOz8zzw2PiKudq/PH/Fo+iNh6j4cPOmBNPZMkcvZoiHDaK+jlbNudiOR9GV3DHaxXS2TNHxyNouqYjBI8fTFGyHJ9O5FYKFjeGgWiXQRijHvtpDk9rlMfYMpfjx6bm2ITlFZ5ShUCiqrBYqWf7eYtHhzTf1NJ1jtY2n4knuHOtmYq7IXM4mU3Lpj5noutZyrnZjD0VIFyzkK4x0RfjCfXv5ylMvMZQK8/ipNAkjSCrUOp9Hx3ooVDwSYZPN3VESYYOX5vIcubCELoIGwHLF5emJBUZ7o3z20O6mdd67f+SaLB2tGaWasWj0FiYXijw/uURPPMSO/ljbfhbFSpShUCiqrFZZs/y9t+zoJWQ0Z31X23hqcfSD23s5PDFPIhyM5gubesu52rXS2jOXC8zly/TGrXr4pnZM4wS8WudzznbJ2i6bu6P1885ky3RFQ5i6IBIyyNoO8bDBcKp9ueprLa633uu3StZ/6dFxfCkZSlpkihWyJYfnzi+xazDeNExK0R5lKBSKBlYbldqq0Q5WeiDthv40ei25koOpCyqeZN9Isn6e5XO1e2IWe4b0FQny2rk2JS1emMlRdn2QsLXa+ZyKGE35iEB7SadQ8XnT9l7gyoyNVmyUuN56xPwaj9s1GKuXEq/1+q2S9UfyZQD2b+4iHjbqRncmV+YL9+1View1oJLZCsXLpF2ytjb0J1tyGIhf6QR++NhU02cQgRT5G7Z21fsdls/VXn5uCMp7P/PtozxyPM279vSzrS/OaG+UZMRkc0+Esf449989xt6hFDn7SjVVMmySKwed0QDjM1m+9csLPHYyzb/6q8M8fKx5pP169JXaUTM2mZLTtNkvb/5rddyDPzqL5/nrun7rogOPiusDQfPjm7b38r5bhxjtiSojsUaUR6FQXAWtPJDPf/9kx6E/tc/UPAJT1/Gl5PxcgfF0ni29kXpiu3EiXqsn/Npgo8bjGmn0eCKmxvGpApoQnJxaouQEKrPbeqMtZS02omHu64fPMzGbb0qe1zb79s1yVOdZ+Mxk7aYy307XPzWdaZmHsJaFB+Hq+jluRJRHoVBsMOmsTcJq3pxadQI3eg2nprOMX86ze1Ocmze1lt1Y7xN+4/lfmMlyfr5ILKQTMXWKjh+o6CIQorWsxUhXpMkjgfVtsLVpf0jZNO3PdtwVm30rT6A3FmKh4DCbszk8Mc9jJ9M8OT5LaFkpce1aDz15jqGkha5Rz0Ocm83TF7foiYU6yn8oOqMMhUKxwQwmw+TKXtNr7Yb+1GTGbxlO8abtvWzri7c0AqemM/zjyRmenpjn8MT8CrnydtTOv3coRTJiMtQVZaw/jqlrWKZAAgvFSnCuZcZsLfpKnahP+xOiPu3PMjROTedWGJtWRmkoFcb1PJ6eWKBccTG0wFCls+UVoauaEd3WF+fAtkC7y/F9ZnJlPvOeXXz20O5XrZ/jekSFnhSKDebjB0frYZzlQ3/Woxe1XHYjpGtIKetP5neMdq3ogWjH1FKJTLFCxZNUPB/pSzwhQZP1+P1yY7YWfaXVrrl82p+UkqWSu8LYNCb5bcfl1HSOxWIw00IIqPiSZNjk1s0pTF1vK3kCQR6ipt81nbHrxynD8PJRhkKh2GBaNax96h03sb0/vi69qOWyG/tGktWeCoGlC05cyrK9P76m8s6QLliyXQwhCBmCcEgjZ3tI6WNqgstZm/l8ha6I0dT4dzUNc43T/s5cLpC1g43/rTt6W4oTvmtPPz89M8/Pzi7QHTX59R09/OLcIobQuH3LSm2rVtfaKF0pRTMq9KRQvALcu3+Eb33iIP/0v/8LvvWJg9y7f6RjjqFTmKcWv6/1VHi+5OJSiZfmi0TNtf0vLIB4SMf1Ja4viZg6lhGMSHV8SdZ2uWM0xZ1jveuWJW9H7Z5MXeeNYz3cNdbL9v44Hz24tWWVU21++DtuHuDtuwcYSEToi1t1basarQzA1YbJFJ1RhkKheJXopLPUSRepMX4vJbi+pD8RZvdgHNPQ17Splz3JW3f1saU7gl+dHX5Tf5y7d/fzb+7ayvtvG+a20Z6XXQbbik731M5oHrmw1PQd7RiI1bWtOhmAV1tX6kZDhZ4UileJ1cIj7cI8jfH7M5fz9dd3DsbbypS3u/Y79lzZYGtruZoy2NWa6drdU7trCkR9Uh8E+YaattV0xu6YJ3k1daVuNJRHoVC8Srzc8Ejj03I6VyYZNpqa9NayqXe6dqcy2FPTmXqD3/JpfGttpmtFu2veviW1Yp2apvGF+/a2nTCoeOVRhkKheJW4mvBIrcz1t35thL3DqabJdWtJ2na6djsjUhvW1M4QXE3ndrtrfvTgVhVCeh2iQk8KxavI1YZHrmYYULtrtyuDXW3I0dWErFYrvVWG4fWFMhQKxTXE1fY2dDpvq2l8hgYnp7P10aHb+6Lky0HI6GpLUlVO4dpBGQqF4hrj1dpgQ7rg6YkF4mGjLsHxy3OL3LU9mMHxao86Vbx2qByFQqFoiSCYyNdI4zsxPScAAAtrSURBVIQ+VZJ646A8CoXiOmGjhwyVPcld27uZmC2StYMBSXuHE5S9K+ZDhY9uDJShUCiuAzZqyFAjtRxEbcgRBL0XAwmzw6cU1yPKUCgU1wGrVSitRitvROUgFDVUjkKhuA7oJA+yGu0a5wCVg1AAa/AohBD/f3t3HmNVecZx/PububPAICCLhEUFW2pValudIN1so7UF12pio03jEhtsonZNG03/sP2nsYlJV6M1bphWCNUajbXujSZW0UGpsrigIoIoiDJsM8xyn/5xzpCLDmdm7jLb/X2SmznnPeeeeZ4chuee877nvrcBZwJbI2JeL9sF/BE4HdgLXBIRL6TbFqbbaoFbIuK6tP184NfAMcD8iGgpON41wGVAN/CjiHi4lATNqkFfQ1V7rhjWbmmlta2L8Y05jpuRPGyXdTXiJ6EN+nfr6Q7gL8CdB9m+CJibvk4CbgROklQL3ACcBmwCnpd0f0SsBVYD5wF/LTyQpGOBC4DjgBnAY5I+ExEHzgJjZgfIuk3Uc8WQz+d57b1d7N7XTVd3no3b9/LMGx+wdVc7Y+pyHDKmbv/0oQOd8tRGtz4LRUQ8JWl2xi7nAHdGRADPSpooaTowG1gfEW8CSFqW7rs2Italbb0da1lE7APekrQemA88M5CkzKpN1oN4v3/0NSaMqaNlw4e0tnfRUCtqVcP2Pft4t7WN+hoxflId+zq7+e8b22lqyLGvM8+kcfWs29K6/4qi3KOqbOQoR2f2TOCdgvVNaVtv7Sf141jP9nKsT5C0GFgMcMQRRwwsYrNRqK9van1vZzv1tSIfQVtnnraObupzNXRE0NEd5PN5duztZFd7F1PG1TN9fMP+kVNA2UdV2cgxYkc9RcTNwM0Azc3NH38uyMxSPf0XAJ1dedrTqU8BIoLOPBw1ZSyrNrUCgSSmjqvn+bc/YldbFyve2s7nZ05g2oQxRY+qspGtHKOeNgOHF6zPStsO1l7MscysSD3f1Dq5qZ69nd3kI/lcVVdbQz6gqSHHB3uSjvAjJ41lQmOOlzbvpKMzz7iGGna3d/Hk6x+wZcfeA47rfozqUY5CcT9wkRILgNaI2AI8D8yVNEdSPUkn9f39ONYFkhokzSHpIH+uDDGaVa2e/osTjzyUGokaQX2uhqb6Wupqa5jcVM+utmQ+6937utnZ3kV9bQ0NdTXkQxzSmGNsXS2r3jlwngnPSV09+jM8dinwDWCKpE3AtUAdQETcBDxIMjR2Pcnw2EvTbV2SrgQeJhkee1tErEmPeS7wZ2Aq8C9JqyLi2xGxRtJyYC3QBVzhEU9mpTtm+gR+e97xCFj97k46uvOMb6xjclMdm3e009GdZ96M8bzb2k7Lho8Y11BDV3cyv/Zh4xvI5/Ns2N5Ga1unH76rQooY+bf3m5ubo6Wlpe8dzapc4Vd9FP6H39MpvW5LK5ctaWF3exeHNOaY1FTP2PocrW2d1NWIU46Z5lFPo4iklRHR3Nd+I7Yz28wGrj8TBv3q9M9y3b9fpakhR2Ouhta2Tvbs6+LqRUdzxvG9DkK0Uc6FwqzK9PWNrz3FYMkzG3l/ZzvTxjdy1SmfcpGoYi4UZvYJZxw/04XB9vOXApqZWSYXCjMzy+RCYWZmmVwozMwskwuFmZllcqEwM7NMLhRmZpbJhcLMzDK5UJiZWSYXCjMzy+RCYWZmmVwozMwskwuFmZllcqEwM7NMLhRmZpbJhcLMzDK5UJiZWSYXCjMzy+RCYWZmmVwozMwskwuFmZllcqEwM7NMLhRmZpbJhcLMzDK5UJiZWSYXCjMzy+RCYWZmmVwozMwskwuFmZllcqEwM7NMLhRmZpapz0Ih6TZJWyWtPsh2SfqTpPWSXpJ0QsG2hZJeTbddXdA+SdKjkl5Pfx6attdJWiLpZUnrJF1TjiTNzKx4/bmiuANYmLF9ETA3fS0GbgSQVAvckG4/FrhQ0rHpe64GHo+IucDj6TrA+UBDRHwOOBG4XNLs/qdjZmbl1mehiIingA8zdjkHuDMSzwITJU0H5gPrI+LNiOgAlqX79rxnSbq8BPhOz68DmiTlgDFAB7BzgDmZmVkZlaOPYibwTsH6prTtYO0A0yJiS7r8HjAtXb4b2ANsATYC10dEVpEyM7MKG/LO7IgIkisJSK5CuoEZwBzg55KO6u19khZLapHUsm3btsEJ1sysCpWjUGwGDi9Yn5W2Hawd4P309hTpz61p+/eAhyKiMyK2Ak8Dzb390oi4OSKaI6J56tSpZUjDzMx6U45CcT9wUTr6aQHQmt5Weh6YK2mOpHrggnTfnvdcnC5fDNyXLm8ETgGQ1AQsAF4pQ4xmZlakXF87SFoKfAOYImkTcC1QBxARNwEPAqcD64G9wKXpti5JVwIPA7XAbRGxJj3sdcBySZcBbwPfTdtvAG6XtAYQcHtEvFSGPM3MrEhKughGtubm5mhpaRnqMMzMRhRJKyOi19v7hYa8M9vMzIY3FwozM8vkQmFmZplcKMzMLJMLhZmZZXKhMDOzTKNieKykbSTPY5TbFOCDChx3OHGOo4NzHB0GO8cjI6LPr7YYFYWiUiS19GeM8UjmHEcH5zg6DNccfevJzMwyuVCYmVkmF4psNw91AIPAOY4OznF0GJY5uo/CzMwy+YrCzMwyVWWhkNQo6TlJ/5O0RtJvetnnF5JWpa/VkrolTSrYXivpRUkPDG70fSs1P0kbJL2cbhuWX8tbhhwnSrpb0iuS1kn60uBnka2UHCUdXdC+StJOST8ZijyylOE8/jR932pJSyU1Dn4W2cqQ44/TtjVDdg4joupeJHNdjEuX64AVwIKM/c8CnvhY28+Au4AHhjqfcucHbACmDHUeFc5xCfCDdLkemDjUOZU7x4L2WpK56Y8c6pzKmSMwE3gLGJOuLwcuGeqcypzjPGA1MJZk/qDHgE8Pdg5VeUURid3pal36yuqsuRBY2rMiaRZwBnBLxYIsQan5jQSl5ChpAnAycGt6rI6I2FHBcItSxvN4KvBGRFTiodSSlCHHHDBGUo7kP9N3KxJoCUrM8RhgRUTsjYgu4EngvIoFexBVWShg/62jVSTzdT8aESsOst9YYCFwT0HzH4BfAvmKB1qkEvML4DFJKyUtrny0xSkhxznANpLZFF+UdEs69e6wU+J57HEBw/iDQLE5RsRm4HqSKZS3kEzD/MjgRD0wJZzH1cDXJE1Ot50OHD4YMReq2kIREd0R8QVgFjBf0ryD7HoW8HREfAgg6Uxga0SsHKRQi1Jsfqmvpu9dBFwh6eQKh1uUEnLMAScAN0bEF4E9wNUVD7gIJZ5HlMxXfzbwj8pGWrwS/hYPBc4hKfwzgCZJ3x+MmAeq2BwjYh3wO+AR4CFgFdA9CCEfoGoLRY/0lsN/SKp4bz7+aewrwNmSNgDLgFMk/a2iQZagiPx6PqkREVuBe4H5lYyxVEXkuAnYVPCp7m6SwjFsFXMeU4uAFyLi/UrFVi5F5PhN4K2I2BYRncA/gS9XNsrSFPn3eGtEnBgRJwMfAa9VNspPqspCIWmqpInp8hjgNOCVXvabAHwduK+nLSKuiYhZETGb5KQ+ERHD6lNMKflJapJ0SM8y8C2Sy99hpcRz+B7wjqSj06ZTgbUVD3qASsmxwLDufyoxx43AAkljJYnkPK6rfNQDU+p5lHRY+vMIkv6Juyod88flBvsXDhPTgSWSakmK5fKIeEDSDwEi4qZ0v3OBRyJizxDFWaxS8psG3Jv83ZED7oqIhwYv9H4r9RxeBfw9vTXzJnDpIMU9ECXlmBb604DLBzHmgSo6x4hYIelu4AWgC3iR4flkc6n/Vu+RNBnoBK4YioEXfjLbzMwyVeWtJzMz6z8XCjMzy+RCYWZmmVwozMwskwuFmZllcqEwM7NMLhRmZpbJhcLMzDL9H314DGTPcQnFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb796681630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(star_shapes[good_star_idx], star_shapes[good_star_idx]/psf_shapes[good_star_idx]*(math.sqrt(2)/2), alpha=0.5)\n",
    "plt.ylim(1.00075, 1.00175)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST_Stack (Python 3)",
   "language": "python",
   "name": "lsst_stack"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}