{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate  optimal extraction error due to cross-dispersion shift"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,u'Flux per pixel (1D only)')"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl81fd95/vXRxsCIVYJbQhjFrNL\ngCWb2km8J9ggbCe2IZlM0+mSZm6cNM20k7a5TdLMo51Ok2mnt83cqR9t773tJAWc2BgBNk5sx7Hj\n2BY2SOyYzWhFYheL1vO5f5wjRRBJHITO+R3pvJ+Phx46y0+/8+Y8hD7n913N3REREQFICTqAiIgk\nDhUFERHppaIgIiK9VBRERKSXioKIiPRSURARkV4qCiIi0ktFQUREeqkoiIhIr7SgA9yonJwcnzlz\nZtAxRERGlPfee++Uu+de77gRVxRmzpzJjh07go4hIjKimNmH0Ryn5iMREemloiAiIr1UFEREpJeK\ngoiI9FJREBGRXioKIiLSS0VBRER6qSiIJIgLbZ1sqDpBKKQtciU4KgoiCeKf3zzG1360m7eOnA46\niiQxFQWRBODuVFY3ALClpiHgNJLMVBREEsD+xlaOtFxiQmYaL+5poqMrFHQkSVIqCiIJoLKmgdQU\n4xsVizh/pZM3D7cEHUmSlIqCSMB6mo4+MieHNaWFTMhMo7K6MehYkqRUFEQCtqv2HHVnr1BRWkhG\nWgoPLy7g5b1NtHV2Bx1NkpCKgkjAKqsbyUhN4eOL8gCoKC3kUkc3rx1oDjiZJCMVBZEAhULO1t0N\n3DsvlwmZ6QCsmDWFnPEZbKlRE5LEn4qCSICqjp/h5IV2KkoLex9LS03hkSUFvHLgJBfbuwJMJ8lI\nRUEkQJU1DYxNT+WBBdOueryitJC2zhCv7D8ZUDJJVioKIgHp6g6xbXcTDyyYxriMq3fGvX3GZPIn\nZPZOaBOJFxUFkYC8deQ0Zy51XNV01CMlxVhdUsDrh1o4f7kzgHSSrFQURAJSWd1A9pg07rktt9/n\nK0oL6ex2tu9tinMySWYqCiIBaO/qZvveJj6+KJ/M9NR+jymZPpEZU8ZRqbWQJI5UFEQC8MahU1xo\n66KitGDAY8yMitIC3jpymlMX2+OYTpKZioJIACprGpg8Lp275+QMelxFaSHdIefFPWpCkvhQURCJ\nsysd3fx430lWLi4gPXXw/4Lz8rKZM228RiFJ3KgoiMTZqweaudzRPWjTUQ8zo6KkkKrjZ2g8fyUO\n6STZqSiIxFlldQO52WO489apUR2/urQAd9iqZS8kDlQUROKota2T1w42s2pJAakpFtXPzM4dz6LC\nCVoLSeJCRUEkjn6y/yTtXaF+J6wNpqK0kF2156g9czlGyUTCVBRE4qiyupGiSWNZPmPSDf3cqiXh\n/gfNWZBYU1EQiZNzlzv42aEWVpcUYBZd01GP4injWDZjknZkk5hTURCJk5f2NNEV8htuOupRUVLI\n/sYLHG5uHeZkIr+koiASJ5U1Ddyak8WiwglD+vlVJQWYoasFiSkVBZE4aGlt5xdHTlMxhKajHnkT\nMrnz1ilsqWnA3Yc5oUhYTIuCma00s4NmdtjM/miQ454wMzezsljmEQnKi3saCTlDbjrqUVFayJGW\nS+xvVBOSxEbMioKZpQLfAx4GFgKfNrOF/RyXDXwZeCdWWUSCVlndwPz8bObmZd/UeR5eHJ7foFFI\nEiuxvFK4Azjs7kfdvQNYDzzaz3H/BfgroC2GWUQC03DuClXHz7K65PrLWlzPlKwM7p6TQ2W1mpAk\nNmJZFIqA2j736yKP9TKzZUCxu28Z7ERm9nkz22FmO1paWoY/qUgM9SxPsbrk5pqOelSUFFB39gq7\nas8Ny/lE+oplUeivN633o42ZpQB/A/yn653I3Z9x9zJ3L8vN7X+XKpFEVVnTQMn0iczMyRqW8318\nUT4ZqSkahSQxEcuiUAcU97k/HejbEJoNLAZ+ambHgRXAZnU2y2hy/NQlaurOUzFMVwkAE8emc8+8\nXLbubiAUUhOSDK9YFoUqYK6Z3WpmGcA6YHPPk+5+3t1z3H2mu88E3gbWuPuOGGYSiautu8Of5lcN\nQ39CXxWlhZy80E7V8TPDel6RmBUFd+8Cnga2A/uBje6+18y+bWZrYvW6IomksrqB8pmTKZw0dljP\n++CCaYxNT9UoJBl2MZ2n4O7b3P02d5/t7n8eeewb7r65n2Pv1VWCjCaHTrZyoKn1pucm9GdcRhoP\nLJjGtt1NdHWHhv38krw0o1kkRrZUN5Bi4bkFsbC6pJAzlzp468jpmJxfkpOKgkgMuDuVNY382uyp\n5GaPiclr3Dsvl/Fj0rR/swwrFQWRGNjbcIFjpy4N66ija2Wmp/LxRXm8tLeJ9q7umL2OJBcVBZEY\nqKxpIC3FWLk4P6avU1FaSGtbF28cOhXT15HkoaIgMszcnS3VjXzstlwmjcuI6Wt9ZE4Ok8alaxSS\nDBsVBZFh9v6Jc9Sfu0JFaWw6mPtKT03h4cUF/HjfSa50qAlJbp6Kgsgwq6xuICMthQcX5MXl9SpK\nCrjc0c2rB5rj8noyukVdFMwsK7IctogMoDvkbN3dyP3zppGdmR6X17xz1lRyxo/RKCQZFgMWBTNL\nMbPPmNlWM2sGDgCNZrbXzL5jZnPjF1NkZHjn2GlaWttjMmFtIKkpxuqSAl492ExrW2fcXldGp8Gu\nFF4DZgN/DOS7e7G7TwM+Snidor80s8/GIaPIiLGlppFxGancP39aXF+3orSAjq4QP9l/Mq6vK6NP\n2iDPPejuv/Kxw93PAD8CfmRm8bk+FhkBOrtDvLi7kYcW5jE2I74trcuKJ1M0aSyV1Y08vmx6XF9b\nRpcBrxR6CoKZfdfMFg12jIjAzw+f4uzlzphOWBtISqQJ6WeHWjh3uSPury+jRzQdzQeAZ8zsHTP7\ngplNjHUokZGosrqR7Mw0PnpbTiCvv7qkkK6Q89KepkBeX0aH6xYFd/9Hd78b+HVgJlBjZj8ws/ti\nHU5kpGjr7OblvU2sXJTPmLRgBuktLprAzKnjNJFNbkpUQ1IjQ1HnR75OAdXAV81sfQyziYwYrx9q\nobW9K66jjq5lZlSUFvKLI+ERUCJDcd2iYGZ/DRwEHgH+wt1vd/f/5u4VwLJYBxQZCbbUNDIlK4O7\nZk8NNEdFaSEhhxf3aP9mGZporhT2ACXu/rvu/u41z90Rg0wiI8rlji5+su8kjyzJJy012EUCbsvL\nZl5etiayyZANOCTVzJZHbu4C5pvZVc+7+/vufj6G2URGhFf2N3OlszuQUUf9qSgt4LsvH6Lh3JVh\n3wZURr/B5in890Gec+D+Yc4iMiJVVjeQN2EM5TOnBB0FCI9C+u7Lh9ha08jvfGxW0HFkhBmwKLi7\nRheJXMeFtk5+erCFz664hZQUu/4PxMHMnCyWFE2ksqZBRUFu2GBXCr3M7C7Cw1F7j3f3f4lRJpER\n4+W9J+noDsVlmewbUVFawF9sO8DxU5eYmZMVdBwZQaIZffSvwHeBjwDlka+yGOcSGREqqxuYPnks\nS4snBR3lKqsi/RtbNGdBblA0VwplwEJ391iHERlJzlzq4OeHT/E7H5vFtQMxglY0aSxlt0xmS00j\nT9+vBY0letEOSY3tRrMiI9BLe5roCnnCjDq6VkVpIQeaWjl0sjXoKDKCRFMUcoB9ZrbdzDb3fMU6\nmEiiq6xuYHZuFgsKsoOO0q+Hl+STYrBFcxbkBkTTfPStWIcQGWmaL7Tx9rHTfPn+uQnXdNRjWnYm\nK2ZNpbKmkd9/6LaEzSmJJZoF8V4nvFJqduRrf+QxkaS1dXcj7iTcqKNrVZQWcuzUJfY2XAg6iowQ\n0Yw+egp4F3gSeAp4x8yeiHUwkURWWd3AgoIJzJmWmE1HPVYuyictxbTshUQtmj6FrwPl7v45d/91\nwusd/WlsY4kkrrqzl3n/xLmEv0oAmJyVwUfn5rClphENIJRoRFMUUty9uc/901H+nMiotLUmvAJp\noo46ulZFaSH1567w/olzQUeRESCaP+4vRUYe/YaZ/QawFdgW21giiauypoGlxZMonjIu6ChReWhh\nHhlpKWpCkqhE09H8h8A/ACVAKfCMu38t1sFEEtHRlovsqb/A6pLEbzrqkZ2Zzn3zctm6u5HukJqQ\nZHBRrX3k7s8Bz8U4i0jC21LTiFl4JdKRpKK0kO17T/LOsdPcNTuYPaRlZIhp34CZrTSzg2Z22Mz+\nqJ/nv2Bmu81sl5m9aWYLY5lH5GZVVjdQPnMK+RMzg45yQ+6fP41xGalUVmtHNhlczIpCZF/n7wEP\nAwuBT/fzR/8H7r7E3ZcCfwX8dazyiNysg02tfNB8MdB9mIdqXEYaDy7I46U9jXR2h4KOIwksqqJg\nZrlmlnuD574DOOzuR929A1gPPNr3AHfvO6Mmi/DmPSIJqbK6gdQU4+HFI3MpsIrSQs5e7uTnh08F\nHUUS2IBFwcK+ZWanCM9oPmRmLWb2jSjPXQTU9rlfF3ns2tf5opkdIXyl8OUBsnzezHaY2Y6WlpYo\nX15k+Lg7lTUN3DV7KjnjxwQdZ0g+dlsO2ZlpakKSQQ12pfAV4G7CE9emuvtk4E7gbjP7/SjO3d9C\nK79yJeDu33P32cDXgP+zvxO5+zPuXubuZbm5N3rBInLzdtef58PTl0fM3IT+jElL5ROL8nl5bxNt\nnd1Bx5EENVhR+HXg0+5+rOcBdz8KfDby3PXUAcV97k8HBhsovR54LIrzisRdZXUD6anGJxaNzKaj\nHhWlhbS2d/H6IV1xS/8GKwrp7v4rjY/u3gKkR3HuKmCumd1qZhnAOuCqJbfNrO/uH6uAD6I4r0hc\nhULOlppG7rktl4njovnVT1x3zZ7KlKwMTWSTAQ1WFDqG+BwA7t4FPA1sB/YDG919r5l928zWRA57\n2sz2mtku4KvA56LMLRI37584S+P5thE56uha6akpPLw4n1f2N3O5oyvoOJKABpu8Vmpm/a23a0BU\ng7TdfRvXLInh7t/oc/v3ojmPSJAqqxvITE/hwQV5QUcZFhWlhXz/nRO8sr95VBQ6GV4DXim4e6q7\nT+jnK9vdR/Y1tEiUurpDbN3dyAPz88gaE9UCAAmvfOYU8iaMUROS9GtIk9fMbPxwBxFJRO8cO8Op\nix0jYpnsaKWmGKuWFPLTgy1caOsMOo4kmKHOaN43rClEElRldQNZGancO29a0FGG1erSAjq6Q7y8\n92TQUSTBDHg9bGZfHegpQFcKMup1dIV4cU8TH1+UT2Z6atBxhtWy4kkUTRpLZXUDT9w+Peg4kkAG\nu1L4C2Ayv9ybuedr/HV+TmRUePNwC+evdI6qpqMeZkZFaSFvHj7FmUvXHUwoSWSwnrP3gU3u/t61\nT5jZb8cukkhi2FLdyMSx6XxkzuicRV9RWsD/ev0IL+1p4jN3zgg6jiSIwT7x/wfgwwGeK4tBFpGE\n0dbZzcv7TvLw4nwy0kbnhfHCggnMys3SKCS5ymBDUg/2N6M58px6p2RU++nBZi62d43qcfxmRkVJ\nIW8fO03zhbag40iCGGyV1GfMbMkAz2WZ2W+a2b+LXTSR4FRWN5IzPoM7b50SdJSYqigtwB227tbK\nqRI22HXx/wT+1Mz2m9mzZvY/zeyfzewN4C3Cnc4/jEtKkTi62N7FKwdO8siSAtJSR2fTUY8507KZ\nn5+tJiTpNWBHs7vvAp6KTFQrAwqAK8B+dz8Yp3wicffK/pO0dYZGddNRXxWlhXxn+0Hqzl5m+uRx\nQceRgF33Y5C7X3T3n7r7v7n7JhUEGe0qqxspmJjJ7TMmBx0lLnr2iNhaoyYk0XwDkaucv9zJ64ea\nWV1SQEpKf/tEjT4zpo6jtHgSlTVqQhIVBZGrbN/XRGe3J03TUY+KkgL21F/gaMvFoKNIwFQURPqo\nrG5gxpRxLCmaGHSUuFpVEp61vUVNSElvsLWPKulnT+Ue7r5moOdERqJTF9t568hpvnDPLMySo+mo\nR8HEsdwxcwqbqxv40v1zku7fL7802DIX341bCpEE8OKeJrpDydd01KOitIA/fWEvB0+2Mj9/QtBx\nJCCDDUl9vee2mY0FZmjkkYxmW6obmDttPPPysoOOEoiHlxTwzc172VLdqKKQxK7bp2BmFcAu4KXI\n/aVmtjnWwUTiqel8G+8eP0NFaWHSNp3kjB/D3XNyqKxpwH3AlmMZ5aLpaP4WcAdwDnontc2MXSSR\n+Nu6uxF3WF0y+pbJvhEVJYV8ePoyu+vPBx1FAhJNUehyd/2GyKhWWd3A4qIJzMpN7v2jPrEon/RU\n07IXSSyaorDHzD4DpJrZXDP7O8JrH4mMCidOX2ZX7TlWlyRnB3NfE8el87G5uVRWN9IdUhNSMoqm\nKHwJWAS0Az8ALgBfiWUokXh69r1aUgweXaqiAPDJ5dNputDGm4f7XTlfRrloikLI3b/u7uWRr6+j\nPZpllOjqDvHsjjrunTeNgoljg46TEB5cOI0pWRlsqDoRdBQJQDRFocrMVvTcMbNPoeYjGSV+9kEL\nTRfaeKqsOOgoCWNMWiqPLyvix/tOcupie9BxJM6iKQqfAf7OzL5jZt8Hfge4P7axROJj/bu15IzP\n4IEF04KOklDWlhfT2e08/3590FEkzqJZOns38OfAF4D7gKfdvS7WwURirbm1jVcONPOp26eTPso3\n07lRt+Vls3zGJNZXndCchSQTzeS1fyLcsVwC/Aeg0sy+GOtgIrH2o/fq6Q45a9V01K915TM40nKJ\n9z48G3QUiaOohqQC97n7MXffDqwAlsc2lkhsuTsbqk5wx61Tkn5uwkBWlRSQlZHK+qraoKNIHEXT\nfPQ33uf60d3Pu/tvxTaWSGy9c+wMx09fZl25rhIGkjUmjTVLC9la08iFts6g40icDFgUzGxj5Ptu\nM6vp87XbzGriF1Fk+G2oqiU7M42HFyf3shbXs7Z8Blc6uzXDOYkMtnT270W+r45HEJF4OX+5k227\nG3mqrJixGalBx0lopdMnMj8/mw1Vtfy7O28JOo7EwYBXCu7eGPn+IeHZzKWEO5vbI4+JjEgvVNfT\n3hVirZqOrsvMWFteTE3defY2aAm0ZBDN6KPfBt4FPgk8AbxtZr8ZzcnNbKWZHTSzw2b2R/08/1Uz\n2xdplnrFzPRRRGLK3fm3d2tZXDSBxUm25eZQPb6siIy0FDaqwzkpRDP66A+BZe7+G+7+OeB24GvX\n+yEzSwW+BzwMLAQ+bWYLrzlsJ1Dm7iXAD4G/upHwIjdqT/0F9jdeYG35jKCjjBiTxmWwclE+z++s\np62zO+g4EmPRFIU6oLXP/VYgmo8MdwCH3f2ou3cA64FH+x7g7q+5++XI3beB6VGcV2TI1ledIDM9\nhTVJuuXmUK0rL+ZCWxfb9zYFHUViLJqiUA+8Y2bfMrNvEv7jfTjS9PPVQX6uiKuLR13ksYH8FvBi\nFHlEhuRyRxebdzXwyJICJo5NDzrOiLJi1lRmTBnH+nfVhDTaRVMUjgCbgJ65Ci8AjUB25Gsg/e1p\n2O98eTP7LFAGfGeA5z9vZjvMbEdLS0sUkUV+1bbdTbS2d7FOTUc3LCUl3OH8i6OnOX7qUtBxJIYG\nG5IKgLv/2RDPXQf0Hd4xHfiVwc5m9iDwdeAed+93SUZ3fwZ4BqCsrEwLsciQbKg6waycLMpnTg46\nyoj0xO3T+e8vH2Tjjlr+88r5QceRGInlKmBVwFwzu9XMMoB1wOa+B5jZMuAfgDXu3hzDLJLkjrRc\npOr4WdaWF2PW30WsXE/ehEzunz+NZ9+ro6s7FHQciZGYFQV37wKeBrYD+4GN7r7XzL5tZmsih32H\n8IY9z5rZLjPbPMDpRG7Kxqpa0lKMTy7XWIabsbZ8Bi2t7bx2UM24o9WgzUeRYaVfdve/GcrJ3X0b\nsO2ax77R5/aDQzmvyI3o6Arxo/freHBBHrnZY4KOM6LdNy+Xadlj2FB1gocW5gUdR2Jg0CsFd+/m\nmmGkIiPNqwdOcupiB2vv0Azmm5WWmsITt0/n1QPNNJ1vCzqOxEA0zUc/N7O/N7OPmtnynq+YJxMZ\nJuuraimYmMnH5uYGHWVUeKqsmJDDj97XXluj0XVHHwF3Rb5/u89jjrbklBGg4dwVXj/Uwpfum0Nq\nijqYh8PMnCxWzJrChqpa/uM9s0nR+zqqRDMk9b54BBGJhWd3hD/NPqnd1YbVuvIZfGXDLt4+epq7\n5uQEHUeGUTQL4uWZ2T+Z2YuR+wvNTJvsSMLrDjkbd9TykTk5FE8ZF3ScUWXl4nwmZKZpV7ZRKJo+\nhf+X8LDSnsViDhHes1kkof388Cnqz13REtkxkJmeyuPLinhpTxNnL3UEHUeGUTRFIcfdNwIh6J1/\noKUSJeFtqKpl8rh0DZ2MkbXlM+joDrFpV33QUWQYRVMULpnZVCLrFpnZCkC7bUhCO32xnZf3NfHJ\n5dMZk6bd1WJhYeEESqZPZP27tfTZxl1GuGiKwlcJL08x28x+DvwL8KWYphK5Sc/vrKez29V0FGNr\ny4s5eLKV6jp9ThwtrlsU3P194B7CQ1N/F1jk7jWxDiYyVO7O+qpals+YxG15gy3kKzdrTWkhY9NT\n2VB1IugoMkyiGX2UCXwZ+C/AnwFfjDwmkpDeP3GWw80XtUR2HGRnprOqpIDNuxq41N4VdBwZBtE0\nH/0LsAj4O+DvCW+t+a+xDCVyM9a/W0tWRiqrSgqCjpIU1pUXc6mjm601jUFHkWEQzYzmee5e2uf+\na2ZWHatAIjejta2TLTWNPLaskKwx0fx6y826/ZbJzM7NYn3VCZ5SH86IF82Vws7IiCMAzOxO4Oex\niyQydJXVjVzp7Gatmo7ixsxYVz6D90+c44OTrdf/AUlo0RSFO4G3zOy4mR0HfgHcY2a7zUwdzpJQ\nNlSdYH5+NqXTJwYdJak8vryI9FRjg2Y4j3jRXF+vjHkKkWGwr+EC1XXn+WbFQu2uFmc548fw0MI8\nnttZzx+unKe5ISNYNENSPxzsKx4hRaKxcUctGWkpPL6sKOgoSWlt+QzOXOrgJ/u0s+5IFss9mkXi\npq2zm+fer2PlonwmjcsIOk5S+sicHIomjWW95iyMaCoKMips39vEhbYu1mn0S2BSU4wny6bz5uFT\n1J65HHQcGaJoJq8t7Oexe2OSRmSI1r9bS/GUsayYNTXoKEmtZ9+KZ9/TrmwjVTRXChvN7GsWNtbM\n/g74r7EOJhKt46cu8Yujp1lbVqxdwAJWNGksH52by7M7aukOaZG8kSjaIanFwFtAFdAA3B3LUCI3\nYuOOWlIMnrhdTUeJYF15MY3n2/jZBy1BR5EhiKYodAJXgLFAJnDM3UMxTSUSpa7uED98r4775k0j\nf6KW5EoEDy7IY0pWBhve1ZyFkSiaolBFuCiUAx8BPm1mP4xpKpEo/fRgC82t7VoiO4FkpKXwqeVF\n/GT/SVpa24OOIzcomqLwW+7+DXfvdPcmd38UeCHWwUSisb6qltzsMdw3f1rQUaSPteXFdIWc595X\nh/NIE01RaDazGX2/gNdjHUzkek5eaOO1g808cft00lM1ujqRzJmWTdktk9lQpV3ZRppo/idtBbZE\nvr8CHAVejGUokWj88L06ukPOU2VqOkpEa8uLOXrqElXHzwYdRW5ANMtcLHH3ksj3ucAdwJuxjyYy\nsFDI2bijlhWzpnBrTlbQcaQfq0oKGD8mTTOcR5gbvuaObM9ZHoMsIlF7+9hpPjx9WburJbBxGWms\nWVrItt2NnL/SGXQcidJ1V0k1s6/2uZsCLAc0AFkCtaGqlgmZaaxcnB90FBnEuvJifvDOCTZXN/Dv\nV9wSdByJQjRXCtl9vsYQ7lt4NJahRAZz7nIHL+5p4vFlRWSma4nmRLakaCILCiawQU1II8Z1rxTc\n/c/iEUQkWpt21tPRFdLuaiNAeFe2Yr65eS976s+zuEibHyW6AYuCmVUCA44lc/c1MUkkMgh3Z31V\nLSXTJ7KwcELQcSQKjy0t4s+37WfjjloVhRFgsCuF797syc1sJfC3QCrwj+7+l9c8/zHgfwAlwDp3\n10xpGVRN3XkONLXy548vDjqKRGniuHQeWZzP8zvr+ZNHFqjJL8ENVhSOufuQGwLNLBX4HvAQUAdU\nmdlmd9/X57ATwG8AfzDU15Hksr6qlrHpqawpLQw6ityAteUz2LSrgRf3NPL4sulBx5FBDNbRvKnn\nhpn9aAjnvgM47O5H3b0DWM81HdTuftzdawAtsCfXdam9i8276llVUkB2ZnrQceQGrJg1hZlTx7Fe\ni+QlvMGKQt+F6WcN4dxFQN/fgLrIYyJDsnV3I5c6urW72ghkZjxVXsw7x85wtOVi0HFkEIMVBR/g\ndrT62+1kSIugmNnnzWyHme1oadEUiWS1oaqW2blZ3H7L5KCjyBA8sXw6qSnGxh1aJC+RDVYUSs3s\ngpm1AiWR2xfMrNXMLkRx7jrCm/P0mE54g54b5u7PuHuZu5fl5uYO5RQywn1wspX3PjzL2vJizLS7\n2kg0bUIm982bxg/fq6OzWy3GiWrAouDuqe4+wd2z3T0tcrvnfjRjAauAuWZ2q5llAOuAzcMVXJLL\nhqpa0lKMTy5XJ+VItq68mFMX23n1QHPQUWQAMVtv2N27gKeB7cB+YKO77zWzb5vZGgAzKzezOuBJ\n4B/MbG+s8sjI1d7VzXM763loYR4548cEHUduwr3zcpmWPYYNVepwTlTXndF8M9x9G7Dtmse+0ed2\nFeFmJZEB/WRfM2cudWh3tVEgLTWFJ8um83//9AiN569QMHFs0JHkGtqZRBLe+qoTFE7M5KNz1Z80\nGjxVVkzI4YfqcE5IKgqS0GrPXObNw6d4sqyY1BR1MI8Gt0zN4q7ZU9mwo5ZQSLuyJRoVBUloz74X\n/jT5ZJlaGUeTteXF1J29wltHTgcdRa6hoiAJqzvkPLujlo/OzWX65HFBx5Fh9IlF+Uwcm65d2RKQ\nioIkrDc+aKHxfJtmMI9CmempPL6siJf3nuTMpY6g40gfKgqSsDZU1TIlK4MHF+QFHUViYG15MR3d\nIZ7fWR90FOlDRUES0qmL7fx430k+tbyIjDT9mo5GCwomUFo8iQ1VJ3BXh3Oi0P82SUjPvV9HV8g1\nN2GUW1dezKGTF9lZey7oKBJ9zlAVAAAMFUlEQVShoiAJp2d3tbJbJjNnWnbQcSSGKkoLGZeRygYt\nqZ0wVBQk4ez48CxHWy7pKiEJjB+TxuqSAiprGrjY3hV0HEFFQRLQ+ndrGT8mjVUlBUFHkThYWz6D\nyx3dbK0Z0iLKMsxUFCShXGjrZOvuBtYsLWRcRkyX5pIEsXzGJOZOG896LZKXEFQUJKFs3tVAW2dI\ncxOSiJmxtryYnSfOcbCpNeg4SU9FQRLKhqpaFhRMYEnRxKCjSBx9cvl00lNNS2onABUFSRh76s+z\nu/48a8uma3e1JDMlK4OPL8znuZ11tHd1Bx0nqakoSMLYuKOWjLQUHltWFHQUCcDa8mLOXe7k5b0n\ng46S1FQUJCG0dXbz/M56Hl6cz6RxGUHHkQB8ZE4ORZPGqgkpYCoKkhBe3NNIa1uX5iYksZQU46my\nYt48fIraM5eDjpO0VBQkIax/t5Zbpo5jxa1Tg44iAXqybDpm4aZECYaKggTuaMtF3jl2hqfKiknR\n7mpJrXDSWO65LZdnd9TR1R0KOk5SUlGQwG3cUUdqivHE7dpdTcKL5DVdaONnH7QEHSUpqShIoDq7\nQ/zwvTrumzeNvAmZQceRBHD//DxyxmewXovkBUJFQQL16oFmTl1s1wxm6ZWRlsKnlk/nlQPNNLe2\nBR0n6agoSGA6u0P877c/ZFr2GO6dlxt0HEkgT5UX0x1yvv+2NuCJN604JnHl7uysPcemnfVsqWnk\nzKUOfv/B20hL1ecT+aXZueP5yJwc/vaVD6isbuDRpUU8tqyQW6ZmBR1t1LORVoXLysp8x44dQceQ\nG3S05SKbdjXwwq56Pjx9mYy0FB5akMejSwt5YEEeqRp1JNe42N7FluoGNu2q5+2jZwBYNmMSjy8r\nYtWSAqaOHxNwwpHFzN5z97LrHqeiILHS0trOlpoGNu2sp7ruPGbwa7Om8tiyIlYuzmdCZnrQEWWE\nqD93hc27wr9LB0+2kpZifOy2XB5bVsRDC/IYm5EadMSEp6Iggbjc0cXLe0/y/M563jx8iu6Qs7Bg\nAo8tK2RNaRH5EzXCSG7O/sYLbNpVzws7G2i60EZWRiqfWJzP48uKuGt2jq46B6CiIHHT1R3izcOn\n2LSznpf3neRyRzdFk8ayZmkhjy0tYl6+9lmW4RcKOe8cO8OmnfVs291Ia3sXudljWFNayOPLilhU\nOEGr7fahoiAx5e7U1J3n+Z31bKlp4NTFDiZkprGqpJDHlhZSPnOKZidL3LR1dvPagWae31nPaweb\n6ex2Zudm8fiyIh5dWkTxlHFBRwycioLExIenL7FpZ7jz79ipS2SkpvDAgmk8urSI++bnMiZNbbsS\nrHOXO9i2u4lNO+t593i4g7rslsk8uqyI1UsKmJyVnKvwqijIsDl9sZ2tuxt5fmc9O0+cA2DFrCk8\nvqyIlYsLmDhWHcaSmOrOXuaFSAf1B80XSU817rltGo8tK+TBBXlkpifPhxgVBbkpVzq6+fH+k2za\nWc/PDrXQFXLm52fz2LIi1pQWUjhpbNARRaLm7uxrvMCmnfVsrm7g5IV2xo9J4+HF+Ty2rIgVs6aO\n+g7qhCgKZrYS+FsgFfhHd//La54fA/wLcDtwGljr7scHO6eKQux0h5y3jpzi+Z31bN/TxKWObgom\nZvZ2GC8omBB0RJGb1h1y3j56mk0763lxTxMX27vImxDuoH5sWRELC0ZnB3XgRcHMUoFDwENAHVAF\nfNrd9/U55v8AStz9C2a2Dnjc3dcOdl4VheHl7uxtuMDzkU9QLa3tZGemsWpJAY8uLeLOW9VhLKNX\nW2c3r+wPd1D/9GAzXSHntrzxPLq0iEeXFjJ98ujpoE6EovBrwLfc/ROR+38M4O7/tc8x2yPH/MLM\n0oAmINcHCZXMRSEUcrrd6Q45oZ7vIej28P2rno883vfYvre7Qs47R0/z/M56jrSEO4zvm5/LY0uL\nuG/+tKRqaxUBOHupg627G9m0s54dH54F4I5bp/QOq05NMVLNMCN8O8VIMet9PCWFPrf7fO/7fOT4\nIK5Eoi0KsVz7qAjou/ZtHXDnQMe4e5eZnQemAqeGO8zGqlqeeePocJ92SPr+AQ+Fwpez3b/yR73P\n85HHYuGOW6fw2x+dxSOLC5g4Th3GkrwmZ2Xw2RW38NkVt1B75jIv7Krn+Z31/Mnzu4f9tcy4qniE\nC8zVxaan4PQUk5QU4/cfvI2K0sJhz9NXLItCf6Xw2r9s0RyDmX0e+DzAjBkzhhRmclYG8/ISZBKV\n0ecXwUiNfMJIsWs+ffR9fsBPH0Zq5Jcp5ZpPKX0fS02hzy9Z+LHZ08ZTpA5jkV9RPGUcT98/ly/e\nN4cDTa00t7aHP6j19wHOne5Q/1fyVz3fzxV7z7lCztUtAL3HctWxk+LwwS2WRaEO6LtI/nSgYYBj\n6iLNRxOBM9eeyN2fAZ6BcPPRUMI8tDCPhxbmDeVHRSRJmRkLCiawoCDoJPETy/WKq4C5ZnarmWUA\n64DN1xyzGfhc5PYTwKuD9SeIiEhsxexKIdJH8DSwnfCQ1H92971m9m1gh7tvBv4J+FczO0z4CmFd\nrPKIiMj1xXSTHXffBmy75rFv9LndBjwZywwiIhI9bXclIiK9VBRERKSXioKIiPRSURARkV4qCiIi\n0mvELZ1tZi3Ah0HnuEk5xGApjxFM78cv6b24mt6Pq93M+3GLu+de76ARVxRGAzPbEc3CVMlC78cv\n6b24mt6Pq8Xj/VDzkYiI9FJREBGRXioKwXgm6AAJRu/HL+m9uJrej6vF/P1Qn4KIiPTSlYKIiPRS\nUQiYmf2BmbmZ5QSdJShm9h0zO2BmNWb2vJlNCjpTEMxspZkdNLPDZvZHQecJkpkVm9lrZrbfzPaa\n2e8FnSloZpZqZjvNbEssX0dFIUBmVgw8BJwIOkvAfgwsdvcS4BDwxwHniTszSwW+BzwMLAQ+bWYL\ng00VqC7gP7n7AmAF8MUkfz8Afg/YH+sXUVEI1t8A/5l+tiBNJu7+srt3Re6+TXiXvmRzB3DY3Y+6\newewHng04EyBcfdGd38/cruV8B/DomBTBcfMpgOrgH+M9WupKATEzNYA9e5eHXSWBPObwItBhwhA\nEVDb534dSfxHsC8zmwksA94JNkmg/gfhD5ChWL9QTDfZSXZm9hMgv5+nvg78CfDx+CYKzmDvhbu/\nEDnm64SbDb4fz2wJwvp5LKmvIAHMbDzwI+Ar7n4h6DxBMLPVQLO7v2dm98b69VQUYsjdH+zvcTNb\nAtwKVJsZhJtL3jezO9y9KY4R42ag96KHmX0OWA08kKT7dNcBxX3uTwcaAsqSEMwsnXBB+L67Pxd0\nngDdDawxs0eATGCCmf1vd/9sLF5M8xQSgJkdB8rcPSkX/jKzlcBfA/e4e0vQeYJgZmmEO9kfAOqB\nKuAz7r430GABsfCnpf8POOPuXwk6T6KIXCn8gbuvjtVrqE9BEsHfA9nAj81sl5n9r6ADxVuko/1p\nYDvhTtWNyVoQIu4G/j1wf+R3Ylfkk7LEmK4URESkl64URESkl4qCiIj0UlEQEZFeKgoiItJLRUFE\nRHqpKIj0w8y6I8Mg95jZs2Y2LvL4W0M830wz2zO8KUWGn4qCSP+uuPtSd18MdABfAHD3u4KNJRJb\nKgoi1/cGMAfAzC5Gvj9uZj+xsAIzO2Rm+ZE1779jZlWR/SF+99qTmdkiM3s3ciVSY2Zz4/zvERmQ\n1j4SGURk+YmHgZf6Pu7uz5vZp4AvAiuBb7p7k5l9Hjjv7uVmNgb4uZm9zNWL230B+Ft3/76ZZQCp\ncfnHiERBRUGkf2PNbFfk9hvAP/VzzJeAPcDb7v5vkcc+DpSY2ROR+xOBuYTXNerxC+DrkTXyn3P3\nD4Y9vcgQqSiI9O+Kuy+9zjFFhNe3zzOzFHcPEV4C+0vuvr3vgZE9AQBw9x+Y2TuEN03Zbma/7e6v\nDmt6kSFSn4LIEESalf4f4DOEF7D7auSp7cB/jCz7jJndZmZZ1/zsLOCou/9fwGagJG7BRa5DVwoi\nQ/MnwBvu/kakmanKzLYS3i5xJuH9MQxoAR675mfXAp81s06gCfh2/GKLDE6rpIqISC81H4mISC8V\nBRER6aWiICIivVQURESkl4qCiIj0UlEQEZFeKgoiItJLRUFERHr9/1cBOGknYfLAAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109c53810>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%pylab inline --no-import-all\n",
    "matplotlib.rcParams['image.origin'] = 'lower'\n",
    "matplotlib.rcParams['image.interpolation'] = 'nearest'\n",
    "\n",
    "FWHM = 2\n",
    "sig=FWHM/2.35\n",
    "flux = 1 \n",
    "x = np.arange(-5,5,1)\n",
    "gaussian_model = flux*np.exp(-x**2/2./sig**2)/sig/np.sqrt(2.*np.pi)\n",
    "plt.plot(x,gaussian_model)\n",
    "plt.xlabel('Pixels')\n",
    "plt.ylabel('Flux per pixel (1D only)')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct gaussian matched filter, offset by some amount"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('Testing matched filter: measured:', 0.99677660015688285, ' input: ', 1)\n"
     ]
    }
   ],
   "source": [
    "offset = 0.1\n",
    "matched_filter = np.exp(-(x-offset)**2/2./sig**2)/sig/np.sqrt(2.*np.pi)\n",
    "matched_filter /= np.sum(matched_filter**2)\n",
    "print(\"Testing matched filter: measured:\",np.sum(gaussian_model*matched_filter),\" input: \",flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.331882884126\n"
     ]
    }
   ],
   "source": [
    "matched_filter = np.exp(-(x-offset)**2/2./sig**2)/sig/np.sqrt(2.*np.pi)\n",
    "print np.sum(matched_filter**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do this for a range of offsets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10ca1cc90>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd4FWX6//H3nYRUSgIJivQmiA0w\nSBUQFYIoqFgAe0NAwIKoa93V3d+u/YsIuIAFGwi4UkRBFFgQRAlIRwQpEnqvIZDk/v1xJtljSBkg\n50zK/bquuXJm5jlzPhmP3Jn2PKKqGGOMMQAhXgcwxhhTdFhRMMYYk82KgjHGmGxWFIwxxmSzomCM\nMSabFQVjjDHZrCgYY4zJZkXBGGNMNisKxhhjsoV5HeB0xcfHa61atbyOYYwxxcrixYv3qGpCQe2K\nXVGoVasWycnJXscwxphiRUQ2u2lnp4+MMcZks6JgjDEmmxUFY4wx2awoGGOMyWZFwRhjTLaAFQUR\neV9EdonIyjzWi4i8LSLrRWS5iDQNVBZjjDHuBPJI4UMgKZ/1nYH6ztQbGBHALMYYY1wI2HMKqjpX\nRGrl06Qb8JH6xgNdKCKxIlJFVbcHIs+an2ZwYOVMJLI8IVEVCI2KpUxMLOFlY4ksG0dUuYrEVKhI\ndGQUISESiAjGGFPkefnwWlVgi998irPslKIgIr3xHU1Qo0aNM/qwg7/Np+WWUQW2S9Vwjkg0xySG\n1JAYUkPLcjKsLCfDypEeWRGp3JDyNZtQtf4lxFcoi4gVEGNMyeFlUcjtX1PNraGqjgRGAiQmJuba\npiAt7nyJjPQXOHpkP6kH95J6+ABpR/dx4shB0o/tJyP1IKQegrRDhJw4SOiJI5Q5eZhyGUeITNtF\nVOpRyh86RJldGbASTmgo66QqOyPrciyuAaFVLiauVmNq1a5HpXKRZxLRGGM852VRSAGq+81XA7YF\n8gNDw8IoH5tA+dgCu//IlaafYO8fq9iz/heOb1tBxN7VNDy6koTts33HN0vggMawWGqyJ6YeaRUb\nEl71YuLrNKFetXOIjQ4v3F/IGGMKmZdFYQrQX0TGAc2Bg4G6nlBYJCycSnWaUKlOkz8t12P72bdx\nGfs2LiF9+0oq7fuVRsdmEnV0iu8E2ULYnFmZJWF12HvOFcQ17kKziy+iQnQZb34RY4zJg/iu8wZg\nwyJjgfZAPLATeBEoA6Cq74rvZPw7+O5QOgbcq6oF9nSXmJioxaJDvMxM9MBm9m1cysFNS9Gdq4jb\nt4yK6bsAWJ1ZkzXlWpBRryMNml7JxdUr2gVuY0zAiMhiVU0ssF2gikKgFJuikBtV0nesYkfyFFj3\nLVUOLSOUTPZpWRaGNGF3lfZUujSJFhfWJ75shNdpjTEliBWF4iB1P4dXzeDAsmnEbfsvZTMOkqHC\nEq3P6rIt0XodubBJS5rUiCMs1B4+N8acOSsKxU1mBpkpi9mzZCqsn0HlI2sB2KqVmC9N2FOlPedc\n2pFOTepSNqLYDYNhjPGYFYXi7tB2jq2ezqHl04jb8QMRmamkajjf0Irt9XrQpl0Sl1SPteckjDGu\nWFEoSdLT0M0L2PvT55RbP4mIzFTWZNbg+5guVGp5O9c2a0iFKLuTyRiTNysKJVXaYY4vGcexH0dT\n8dCvHNMIvtaWpNTtQZu2HbmsVkU7ejDGnMKKQmmwdQn75o0k5jff0cPqzJrMjO5MbPNedL38AuJi\n7GE5Y4yPFYXS5Pgh0n75nGM/jibOOXqYpq3YXOs2WrW9mpZ14+3owZhSzopCaaQK25ZwYN4oon+b\nRHhmKqsyazIjsjPlL+9Ft+YNSShnzz8YUxpZUSjtjh/i5NLPObpgFLGH1nJUI/hKW7PtwofokdSO\nKhWivE5ojAkiKwrGRxW2LuHg/FFE//oFkpnOZL2ClIsfpkendpxT3np0NaY0sKJgTnV4J4e/f53I\nZR86xaEt2y55mNs6taWydfdtTIlmRcHk7fAODn3/OlHLxiCZ6UzStuy4tD+3dWxr1xyMKaGsKJiC\nHd7Boe9eI2r5R5CZwWRty64mA7jtmjZUsg75jClRrCgY9w5t5+B3rxG94iPIzGQybdnTZCC3XdPG\nnnUwpoSwomBO36HtHJj5GjErs4pDO/Y29R052KhxxhRvVhTMmTu0jf0zX6Xsyk8gM5NJtOfgZQO4\n5eo2NlqcMcWUFQVz9g5tY/+3r1B25aegmUySK8ls9ww3t2tKqI0SZ0yx4rYo2MgtJm/lzyPu5iGU\neWwZhy+6kxuZQ+c5XXjv9cEs2bTL63TGmAAIaFEQkSQRWSsi60Xk6VzW1xSR70VkuYjMEZFqgcxj\nzlCFqlS8ZQih/Rdy/NxEeh8bRcz7VzJqzAfsO3rC63TGmEIUsKIgIqHAMKAz0AjoKSKNcjR7HfhI\nVS8BXgL+Gag85uxJfH0q95lK6s2fUDkqkwc3Pkrya12ZNOdnMjOL12lIY0zuAnmkcDmwXlU3qOoJ\nYBzQLUebRsD3zuvZuaw3RY0IURddT9wTS9jT7Anas5iOs6/j09cfYeVmO6VkTHEXyKJQFdjiN5/i\nLPO3DOjuvL4RKCcilQKYyRSWMlHEd3meMo8ks/+8ttx5bAxl32/Dx2NGcvDYSa/TGWPOUCCLQm63\np+Q8x/AE0E5EfgHaAVuB9FM2JNJbRJJFJHn37t2Fn9ScMYmrSdWHJnL01onEREZy58bBLH+1I1//\nd76dUjKmGApkUUgBqvvNVwO2+TdQ1W2qepOqNgGedZYdzLkhVR2pqomqmpiQkBDAyOZMxTS6hoTB\nyWxv/hyXyRqumtWVL994iF//2OF1NGPMaQhkUVgE1BeR2iISDvQApvg3EJF4EcnK8Bfg/QDmMYEW\nFk6VzoOJfPQXtlfrTPejn1PhvZZM+Ggoh1PtLiVjioOAFQVVTQf6AzOANcB4VV0lIi+JSFenWXtg\nrYj8BpwD/CNQeUzwhFSoQq0HP+Fwr6/QqErcsuE5fn31SmbNnUNxe1jSmNLGnmg2gZWZQcp3w6nw\n47+IzExlelwvWt77T+IrlPM6mTGlij3RbIqGkFCqdRxA9KDlbKqSxPUHPmbPW21YMH+O18mMMbmw\nomCCIrRsJer3+YyUpPepLAdo9u1NzBj2OIePHvM6mjHGjxUFE1TVWnSn7GOLWZdwNZ12v8fW11uz\nbPECr2MZYxwFFgURqZzLsgaBiWNKg/Dy8TTqP57frxzBObqXC6Zcx+xRT3E8Lc3raMaUem6OFOaJ\nyK1ZMyIyCPgycJFMaVG3XS/CH/mZX2PbceXWd9n0amvWrVzkdSxjSjU3RaE9cKeITBCRucD5+Po1\nMuasxcSdyyWPfcnK1m9zTsZOakzozPwxz5F+0rrKMMYLBRYFVd0OTAdaArXw9Wp6JMC5TClz0TV3\nI/0WsrpcC1pvHMrvr7QhZd0yr2MZU+q4uaYwE2gOXARcC7wlIq8HOpgpfWIrV6XJoKkkJ77OOekp\nxH9yFcljX0IzTukOyxgTIG5OHw1T1btU9YCqrgRaAaf0T2RMoRAh8boHOdH7R1ZFJZK49g3WvXIF\nezav9jqZMaWCm9NHk3LMp6vqy4GLZAxUPq8GTZ/8mnkX/4Nz0jYT80E7Vn7xL8jM9DqaMSWam9NH\nh0XkkDMdF5EMEbEjBRNwEhLCFd37c+DeuawMv5SLVvyTtW8mcfyA9bxqTKC4OVIop6rlnSkS36A4\nwwIfzRifmrXq0eTJGXxb+ylqHV7C0SEt2fbLdK9jGVMinfYTzc7ppA4ByGJMnsLCQul49zMs6/wF\nhzSScyf14LexT4NdhDamUIUV1EBEbvKbDQESOXUENWOC4vIW7dhedx5z3nuIDmtHsOnNBZx73ydE\nVqrhdTRjSgQ3RwrX+02dgMNAt0CGMiY/VRLiuWLwBCbXfpGEI7+S9k5rdiVP9jqWMSWCjadgirX5\nPy2k4jcPcQGb2FT/Hmrd9hqEhXsdy5gix+14CnmePhKRoeRzmkhVB55hNmMKTevmLdhSezZfvT+A\n69Z9yNY3fyLhnk8Jr1zX62jGFEv5nT5KBhbnMxlTJFSvXJFrnviIcXX+Sdmjm0kf0Zq9P431OpYx\nxVJ+F5qvVtU7ReQRVR0StETGnIGIsFB63NWPWQsvp+I3fWn8TR+2rp1F1R5DIDza63jGFBv5HSlc\nJiI1gftEJE5EKvpPbjYuIkkislZE1ovI07msryEis0XkFxFZLiLXnukvYgxAhxaJxPabyfjIm6m6\nYTy732rNye2rvI5lTLGRX1F4F1/vqA059dRRgVd6RSQU30NunYFGQE8RaZSj2XPAeFVtAvQAhp/u\nL2BMTrXOiaXroJG8X+dNOLaHzH9fycEfRkMxu6nCGC/kWRRU9W1VvQB4X1XrqGptv6mOi21fDqxX\n1Q2qegIYx6m3sipQ3nldAdh2Br+DMaeILBPKfXfdz+LOU1mi9anw3SB2fXgnpB32OpoxRZqbbi76\nnuG2qwJb/OZTnGX+/grcISIpwNfAgDP8LGNyldSiMQl9v+aDiDuptOkr9g5pS+budV7HMqbIOu1u\nLk6D5LIs5/F7T+BDVa2Gb6yGj0XklEwi0ltEkkUkeffu3QGIakqyeudWoMegIbxb4w04upu0Ee1I\nXTnV61jGFEmBLAopQHW/+WqcenrofmA8gKr+CEQC8Tk3pKojVTVRVRMTEhICFNeUZFHhofS77z6+\nbfM569MrEzXxDg5M+6t1xW1MDoEsCouA+iJSW0TC8V1InpKjzR/AVQAicgG+omCHAiYgRISe17Tm\nUK9pTOJKYhe9xd7RN0Dqfq+jGVNk5FkUcoyjcMpU0IZVNR3oD8wA1uC7y2iViLwkIl2dZoOAB0Vk\nGTAWuEeLW78bpthp3bAqlw34lKFRfSm39QcODmmD7ljpdSxjioQC+z4SkZeAHcDH+K4T3A6UU9VX\nAx/vVNb3kSksx06kM+Ljz7jjj+eJDTkOXYcS0eRWr2MZExBu+z5yc/qok6oOV9XDqnpIVUfgG2jH\nmGItOjyMx++7kxmtP2dFRg0iJj/IoclP2RgNplRzUxQyROR2EQkVkRARuR3ICHQwY4JBRLirYwuO\n9ZzEZyRR/pd3OTCyCxzd43U0Yzzhpij0Am4FdjrTLc4yY0qMthdUpfWA93kt6lEidyzmyNut0RTr\n99GUPm4eXtukqt1UNV5VE1T1BlXdFIRsxgRVzUox9Hv0ed6sPpQDx9NJf68TJxaN8TqWMUFVYFEQ\nkfNF5HsRWenMXyIizwU+mjHBFxMRxl/u78GM1uP4KaMB4dMGcmRif0hP8zqaMUHh5vTRKOAvwEkA\nVV2O75kDY0okEeH+js042WMi79GNsis/5vC7neCQdc1lSj43RSFaVX/OscxuzzAl3pWNqtCh/wj+\nFvUUIbtXc+ydNugfP3kdy5iAclMU9ohIXZx+i0TkZmB7QFMZU0TUjo9h0KNP8mr14ew6HkbGB11I\nX/KZ17GMCRg3ReFh4N9AQxHZCjwK9AloKmOKkLIRYbx4X3e+bvEJP6fXJ2xKX1K/fg4y7c5sU/Lk\nWxScHksTVfVqIAFoqKptVHVzUNIZU0SEhAj9rr2cPTeO5bPMq4n6eShHxtxm4zOYEiffoqCqmfj6\nL0JVj6qq/R9gSrWuTWvR6IHRvBpyP5GbvufI8Ktgv/2NZEoON6ePZorIEyJS/XTHaDamJGpcI447\nH/kHL5b/G5kHtpA6vB26+UevYxlTKNwUhfvwXVeYy2mM0WxMSValQhTPDujHmzWHsz0tgowPryN9\n8SdexzLmrLl5orl2LpObMZqNKdGiw8N44Z4b+LrFJ/yY3pCwqQ+TOu0ZuwBtijU3TzRHi8hzIjLS\nma8vItcFPpoxRV9IiND/2mbsv/FTPsnsSNSiYRwZcyscL3DIEWOKJDenjz4ATgCtnPkU4O8BS2RM\nMdS1aS0ufGAk/wp5kMhNszgyogPs2+h1LGNOm5uiUNcZUCerm4tUfIPtGGP8NKkRx10DX+b5ci+T\nfmAbx0e0Rzf94HUsY06Lm6JwQkSi+N8TzXUB6x3MmFycFxvF8wP78GaNEWxNiyJjTDfSradVU4y4\nKQp/BaYD1UXkU+B74MlAhjKmOIsOD+Ov93bl6+afsCD9AsKmDeT41KfsArQpFtzcffQtcBNwDzAW\n3xPOc9xsXESSRGStiKwXkadzWf+WiCx1pt9E5MDpxTemaAoJEQZ0SeTAjZ8yJjOJyMXvcvTD7nD8\noNfRjMmXm7uPpgAdgTmq+pWquhqnUERCgWFAZ6AR0FNEGvm3UdXHVLWxqjYGhgL/Od1fwJiirGvT\nmlz8wLv8I+QhwjfP5ciIq+HAFq9jGZMnN6eP3gCuAFaLyAQRuVlEIl2873JgvapuUNUTwDigWz7t\ne+I7EjGmRGlaI457Bv6NF8r9jcwDKRwb3h62LvE6ljG5cnP66L+q2g+oA4zEN17zLhfbrgr4/0mU\n4iw7hYjUBGoDs/JY31tEkkUkeffu3S4+2piipWpsFM8O6MurVYey97hwcnQSGauneh3LmFO4OVLA\nufuoO74us5sBbm6nyO22Vc2jbQ9goqrmeiVOVUeqaqKqJiYkJLiJbEyRUzYijL890J3xjT9kZUZ1\nZPydpM17GzSv/y2MCT431xQ+B9YAHfBdI6irqgNcbDsFqO43Xw3IazzDHtipI1MKhIYIg25qw5qO\nnzI983Iivn+eo/8ZCBk2mKEpGtw+0VxXVfuo6iynO203FgH1RaS2iITj+4d/Ss5GItIAiAOsm0lT\navRq05DoXh8xWrsRs+IjDn1wk3WNYYoEN0Xhe+BhEZnoTANEpExBb1LVdHxjMczAd6QxXlVXichL\nItLVr2lPYJyqHUOb0qV9w3O5ou8w/hXWj+gt8zg84io4mOJ1LFPKSUH/FovIaKAM/7uOcCeQoaoP\nBDhbrhITEzU52XruNiXH7sNpDB09iicO/J3QiBii756IVG3idSxTwojIYlVNLKidmyOFZqp6t3Pq\naJaq3ovvYrMxphAklIvgmf59ebvWMPanwcnRSaSv/srrWKaUclMUMpz+jgAQkTqAPa9vTCGKLBPK\nM3ffxNRmH7M6oyoh4+8g9b92Z5IJPjdFYTAwW0TmiMh/8T1LMCiwsYwpfUJChL7XtWJTl8/5NrMZ\nUbOf5/B/HrE7k0xQuXl47XugPjDQmRqo6uxABzOmtLqheX3i7v6MD+hGuRVjOPB+d7szyQSNm+cU\nHgaiVHW5qi4DokWkX+CjGVN6Na+bQPv+I3gj4mHKpszl4HC7M8kEh5vTRw+qanbvpaq6H3gwcJGM\nMQC142O4/5G/8kqlvyMHt3DknXbotl+8jmVKODdFIUREsruscHo/DQ9cJGNMltjocAb368vIeiM4\ncAJOjErixJrpXscyJZibojADGC8iV4lIB3zdUdi30pggCQ8LYdAdNzCz5Sf8lnEuoZ/35Oj8kV7H\nMiWUm6LwFL47jvoCD2MjrxkTdCLCvUktSblhInMzLyVm5mAOTP4LZLrtdcYYd9zcfZQJfAg8q6rd\nVfXfefVmaowJrM5N61P+3glMkE7E/jKcPR/dCSePex3LlCBu7j7qCizFOWUkIo2d0diMMR64rHYC\nzR/+gH9H3E38pq/YM6IzHNvndSxTQrg5ffQivlHUDgCo6lKgVgAzGWMKUCM+hh6PvMGQuGcot3cF\n+95uh+7d4HUsUwK4KQrpqmqjjRtTxFSILkPfhwczus7/Qeo+jgy/khObFnodyxRzborCShHpBYSK\nSH0RGQosCHAuY4wL4WEh9Lvrdr6+/GP2noyAD6/n8C9feh3LFGNuisIA4EIgDfgMOAg8GshQxhj3\nRIQ7unRg7XX/YZXWJGbyvez57i2vY5liKt+i4Dyo9jdVfVZVmznTc6pqtzsYU8R0uvwi9K4pzKI5\n8T/8lR3jHoFMu1HQnJ58i4Jz6+llQcpijDlLTeuex/n9v2B8mRs499cP2TbyZjhx1OtYphhxc/ro\nFxGZIiJ3ishNWVPAkxljzkiN+LJ0emw071d4mHO2z2bH21ejh3d6HcsUE26KQkVgL9ABuN6Zrgtk\nKGPM2akQXYY7BvydMTX/H+UP/86+t9tyYsdqr2OZYqDAMZrPauMiScAQIBQYrar/yqXNrcBfAQWW\nqWqv/LZpYzQb456qMmHqVK5c3J+okAz0tk8o1/BKr2MZDxTmGM1nGiAUGAZ0BhoBPUWkUY429YG/\nAK1V9ULsriZjCpWIcGvXrizr9AU7MisQMe5mds3/2OtYpggLWFHA9xT0elXdoKongHFAtxxtHgSG\nOWM0oKq7ApjHmFLr6lbNOHL7NJZzPpVn9mfL5Jdt/GeTqzyLgog84vxsfYbbrgps8ZtPcZb5Ox84\nX0Tmi8hC53RTbll6i0iyiCTv3r37DOMYU7o1Pr82Cf2m8V1YW6r/8jobP3zAxn82p8jvSOFe5+fQ\nM9y25LIs558mYfjGf24P9ARGi0jsKW9SHamqiaqamJCQcIZxjDE1K1ek2WNf8GXZntTePJGNQ69D\nbfxn4ye/orBGRDYBDURkud+0QkSWu9h2ClDdb74asC2XNpNV9aSqbgTW4isSxpgAqRATTpdHhzP+\nvCepvv8ntr51JWn7thT8RlMq5FkUVLUn0AJYz/9uRc26HfV6F9teBNQXkdoiEg70AHJ2uT0JuBJA\nROLxnU6yrh6NCbDwsBBuefAZvr5kCLHHUzj8TnsOblrqdSxTBBT0RPMOVb0U2A6Uc6Ztqrq5oA2r\najrQH99wnmuA8aq6SkRecsZowFm3V0RWA7OBwaq698x/HWOMWyJC1+53sfiqsaRnZBL6YWe2L/na\n61jGYwU+pyAi7YCPgE34rhNUB+5W1bkBT5cLe07BmMK3fPUqIsf3pDYp/NHqn9Tt+JDXkUwhK8zn\nFN4EOqpqO1VtC3QCrAtGY0qQSxpdSORD37Is9GLqLniSNWOftltWSyk3RaGMqq7NmlHV34AygYtk\njPFCjSrnUv/Rb5gT1ZEL1o5g9fDb0fQ0r2OZIHNTFJJF5D0Rae9Mo4DFgQ5mjAm+CuWiafX4OL6p\nfD+Ndk9j3ZtJpB2x8Z9LEzdFoS+wChgIPAKsBvoEMpQxxjvhZUJJ6vsGMxu+TK2jy9j5f+05sO13\nr2OZIAloh3iBYBeajQmeBd99yUXz+nEyJILUW8ZSrVFLryOZM+R5h3jGmOKv1dU3suWmSaRpGBXH\nd+PXuRO8jmQCzIqCMSZfF17anIz7ZrIlpBr1v3+QJRNe8TqSCSArCsaYAlWvUZtzB87il6gWNF31\n/0h+tzeZ6daZXklUYFEQkfNFZJSIfCsis7KmYIQzxhQdFWJjuXTQVH6Iv5XEHZ+z4q3rST1inemV\nNGEu2kwA3gVGARmBjWOMKcrKlClD64dHMn9sHVqsfYWNb7Un9oH/EF+lltfRTCFxc/ooXVVHqOrP\nqro4awp4MmNMkSQitO71F5Zd8S5V0lNI//dVbFj5k9exTCFxUxSmikg/EakiIhWzpoAnM8YUaU2v\n7sH27pMIIZPKE7qxbPYXXkcyhcBNUbgbGAwswPck82LAHhQwxlDvklbwwPfsCjuXC+c8wILPX/M6\nkjlLBRYFVa2dy1QnGOGMMUVf5Wp1OPfR2ayObkarNX/nh+F9SLc7k4otN3cflRGRgSIy0Zn6i4h1\niGeMyRZdLo4LB01jceXutNk1ll/euIHDhw96HcucATenj0YAlwHDnekyZ5kxxmQLDSvDZX3fY8kF\ng7ns2A9s/b+r2ZpS4HhcpohxUxSaqerdqjrLme4FmgU6mDGmGBKh6W3Psbb9cGpmbEJGX8XqZT97\nncqcBjdFIUNE6mbNiEgd7HkFY0w+LriyF3tunkQE6VT7Tzd+nGl3JhUXborCYGC2iMwRkf8Cs4BB\nbjYuIkkislZE1ovI07msv0dEdovIUmd64PTiG2OKquoXtSa09/ccKJNA4g8P8t2nr1PcemUujVx1\nnS0iEUADfGM0/6qqBQ7HJCKhwG/ANUAKsAjoqaqr/drcAySqan+3ga3rbGOKl7Qj+9k04mYaHE3m\nu4o9adNnKJHhdq9KsJ1119ki0sH5eRPQBagH1AW6OMsKcjmwXlU3qOoJYBzQzU14Y0zJEVE2jvMf\n+4aVVW7i6n1jWfr69ezeu9frWCYP+Z0+auf8vD6X6ToX264KbPGbT3GW5dRdRJY7t7tWd7FdY0wx\nI2HhXNT7fdZc+gzN0hay/50OrP1tjdexTC7y7BBPVV90Xr6kqhv914lIbRfbltw2m2N+KjBWVdNE\npA8wBuhwyoZEegO9AWrUqOHio40xRY4IF9z4FJuqNKDq9D4c+zSJH68aRcu2Hb1OZvy4udCc220D\nE128LwXw/8u/GrDNv4Gq7vW7PjEK3zMQp1DVkaqaqKqJCQkJLj7aGFNU1WpxA8fvnk5GaCRNvu/F\njM+H2QXoIiS/awoNRaQ7UEFEbvKb7gEiXWx7EVBfRGqLSDjQA5iS4zOq+M12Bex40phSoFLtxsQO\nnMvW6IZ0WvMMM955hOMnrGuMoiC/I4UG+K4dxPLn6wlNgQcL2rCqpgP9gRn4/rEfr6qrROQlEenq\nNBsoIqtEZBkwELjnTH8RY0zxEhl7DnUe/44151xH0t4xLHrjRnbt2+91rFKvwFtSRaSlqv4YpDwF\nsltSjSlhVFn7xcvUX/Emv4bUJfT2sTSod77XqUqcs74l1U8fEYn123CciLx/VumMMSaLCA1ufoGU\nTiOprVso/3ES8+fZiL9ecVMULlHVA1kzqrofaBK4SMaY0qhGq1tJvWMaZUKFJt/14OsJo+wCtAfc\nFIUQEYnLmnFGXXMztrMxxpyWivWaUbb/XHZH1SFp5WC+Gv6kXYAOMjdF4Q1ggYi8LCIv4xuB7dXA\nxjLGlFaRFatS4/FZ/F75Gq7fPZIFb9zGrv02NkOwuBl57SPgZmAnsAu4SVU/DnQwY0zpJeHR1O83\nnvUX9KdD2ndsf7sjv/6+wetYpYKbIwVUdRUwHpgMHBERe6zYGBNYItS77R9s6fAODXUDZT/qyLwf\n5nidqsRzMxxnVxFZB2wE/gtsAr4JcC5jjAGgets7OdprEjGh6TSdeSuTP3uXzEy7AB0obo4UXgZa\nAL+pam3gKmB+QFMZY4yfiudLcL8zAAAQZklEQVS3Jrr/PPZG16Hbb0/x1ZD+HEotsAd/cwbcFIWT\nqroX311IIao6G2gc4FzGGPMnERWrU/3x2aw/rytdD37CijeuY8PW7V7HKnHcFIUDIlIWmAt8KiJD\nALtHzBgTdFIminoPfsSmZi/QPD0ZHXkVC3/+yetYJYqbotANOAY8BkwHfsfXB5IxxgSfCLW6DGJ/\n9/HEhxym0bQbmDpxjD3oVkjyLQrOkJqTVTVTVdNVdYyqvu2cTjLGGM8kXHwNEX3ncjiyCl1WPMLk\nYU9yLO2k17GKvXyLgqpmAMdEpEKQ8hhjjGuRlWtz3uNz2XjO1dywZyQ/v34TKTv3eB2rWHNz+ug4\nsEJE3hORt7OmQAczxhg3JKIsdftOYMOlg2h7ch5HRlzF4qXLvI5VbLnpw2iaMxljTNEkQp0bX2BH\n9Uup/lUfjn/ZhW+2vEnSdbcgktvIwCYveRYFEamhqn+o6phgBjLGmDN1bmI3jp53Pic+uJVrkh9i\n0taVdL7vBSLDrQ9Pt/I7fTQp64WI5DZOszHGFDkx513AOY/N44+KrbhxxxDmvdmTHXutQz238isK\n/sdcdQIdxBhjCktIdCx1Bkzl94Z9uOb4t+x652qWr7Eh4N3IryhoHq+NMaboCwmhbo9X2HrNu9TT\nPzhnXBLTv5lszzMUIL+icKmIHBKRw8AlzutDInJYRA652biIJInIWhFZLyJP59PuZhFRESlw/FBj\njDkdVVv3JP3eGUiZKK5aeC+T3n2e1DTrlCEveRYFVQ1V1fKqWk5Vw5zXWfPlC9qw8+DbMKAz0Ajo\nKSKNcmlXDhgI2LPqxpiAKF+zMfGP/8gflVpz486h/Px6NzZt2+l1rCLJ1XgKZ+hyYL2qblDVE8A4\nfF1m5PQyvpHcjgcwizGmlAuJjqNu/8n8fukTtDk5H/13e35YMM/rWEVOIItCVWCL33yKsyybiDQB\nqqvqVwHMYYwxPiEh1L3xefbdNIG4kGM0ndGdyR//H+kZmV4nKzICWRRye2Ik+wqPiIQAbwGDCtyQ\nSG8RSRaR5N27dxdiRGNMaZRwyTVEDZzPrrIN6Pb7i3z/xl3s3u/qUmmJF8iikAJU95uvBmzzmy8H\nXATMEZFN+AbymZLbxWZVHamqiaqamJCQEMDIxpjSIiKuGrUen8Vvde6h07Gp7Hy7A8tXrvA6lucC\nWRQWAfVFpLaIhAM9gClZK1X1oKrGq2otVa0FLAS6qmpyADMZY8z/hJbh/LuGsOXqf1NbU6g2IYnp\nkz8t1betBqwoqGo60B+YAawBxqvqKhF5SUS6BupzjTHmdFVv04PMB2dzLCKBjkseZto7j3Lk+Amv\nY3lCiltFTExM1ORkO5gwxhS+zLSjrHu/Nw12fsXPoU2Iv2sMdWrW9DpWoRCRxapa4LNggTx9ZIwx\nxUpIRAwN+nzC783/TuOMFUS+34F5s6d7HSuorCgYY4w/Eep2HsChntMICw2h+ZxeTHvvJU6czPA6\nWVBYUTDGmFzEN2hB3OML2VyhOV22vMGPb3Rnx+6SP6qbFQVjjMlDmbKVqP/oNH5t9AhtUueQOuwK\nFs6f7XWsgLKiYIwx+QkJoeGtL7HzxvGUk+M0+fZmvhn9AsdPlMxO9awoGGOMC+c17ki5x35iU2xL\nOqcMYdlrSWz8Y7PXsQqdFQVjjHEponxlGjw6lbVNX6DxyaWUfa8tc6ZPKFEPu1lRMMaY0yFCg66D\nOHzHDNLKlKftjw/y3dB+HD56zOtkhcKKgjHGnIH4epdRZfBCVle5gWv2fcaWN9ry6+rlXsc6a1YU\njDHmDIVGxHBRnw9Z334Y1TK3Ue3zjswaP4zMzOJ7OsmKgjHGnKV67e9AH5rHjsi6dFj9DPPfuJU9\n+/Z6HeuMWFEwxphCUKFKXeo+OYfldXrT+shMjr7dmqU/zfE61mmzomCMMYVEQstwyV2vsaXr50TL\nCRp9fROzPniRk+nFp4sMKwrGGFPIal7WiZhHFrKufEs6bP4/lr3aia0pf3gdyxUrCsYYEwDRsZW5\n8PGvWHHpc1yctpTw0Vfw44yxRf6ZBisKxhgTKCJcfONg9vWaTmpoBVr+2If5b93Ovv37vE6WJysK\nxhgTYFUaJHLekwtZUv1uWh38mtQhzVk8Z0rBb/SAFQVjjAmCsIhomt7/Nn/c8B8ICeWyOXcy7+37\nOHT4oNfR/sSKgjHGBFGtJh1IeGIRS869lSv2fcGBN5uz/MdvvY6VLaBFQUSSRGStiKwXkadzWd9H\nRFaIyFIR+UFEGgUyjzHGFAXh0eVo2mcU65I+I5wMLpx+K/NGPEzqMe/7TwpYURCRUGAY0BloBPTM\n5R/9z1T1YlVtDLwKvBmoPMYYU9TUb9GFCo//zLKE67hi5ydsf70Fq5fM8zRTII8ULgfWq+oGVT0B\njAO6+TdQ1UN+szFA0b5XyxhjCllUuTia9v+E1R1GUz7zEPUnd2Pe6CdISzvuSZ5AFoWqwBa/+RRn\n2Z+IyMMi8ju+I4WBuW1IRHqLSLKIJO/evTsgYY0xxkuN2t5CxCM/szK2A1ekjGLzq61Zv2pR0HME\nsihILstOORJQ1WGqWhd4Cngutw2p6khVTVTVxISEhEKOaYwxRUO5uMo0eWwiK1oPJSFjF9XHd2b+\nmBdIP3kyaBkCWRRSgOp+89WAbfm0HwfcEMA8xhhTLFx8zV1Iv4X8Wq45rTcO4bdX2rJ53YqgfHYg\ni8IioL6I1BaRcKAH8KenNUSkvt9sF2BdAPMYY0yxEVu5KpcO+opfEl+hWvpmEj65ikXTRgf8cwNW\nFFQ1HegPzADWAONVdZWIvCQiXZ1m/UVklYgsBR4H7g5UHmOMKXZEaHJdH070ns+6mKZUrH5h4D+y\nqHfOlFNiYqImJyd7HcMYY4oVEVmsqokFtbMnmo0xxmSzomCMMSabFQVjjDHZrCgYY4zJZkXBGGNM\nNisKxhhjsllRMMYYk82KgjHGmGzF7uE1EdkNbD7Dt8cDewoxTmGzfGfH8p29op7R8p25mqpaYI+i\nxa4onA0RSXbzRJ9XLN/ZsXxnr6hntHyBZ6ePjDHGZLOiYIwxJltpKwojvQ5QAMt3dizf2SvqGS1f\ngJWqawrGGGPyV9qOFIwxxuSjxBQFEUkSkbUisl5Ens5lfYSIfO6s/0lEavmt+4uzfK2IdPIo3+Mi\nslpElovI9yJS029dhogsdaYpOd8bpHz3iMhuvxwP+K27W0TWOVNABkpyke8tv2y/icgBv3XB2H/v\ni8guEVmZx3oRkbed/MtFpKnfuoDuPxfZbncyLReRBSJyqd+6TSKywtl3ARvIxEXG9iJy0O+/4wt+\n6/L9bgQp32C/bCud71xFZ11Q9mGhUdViPwGhwO9AHSAcWAY0ytGmH/Cu87oH8LnzupHTPgKo7Wwn\n1IN8VwLRzuu+Wfmc+SNFYP/dA7yTy3srAhucn3HO67hg58vRfgDwfrD2n/MZbYGmwMo81l8LfAMI\n0AL4KYj7r6BsrbI+E+iclc2Z3wTEF4H91x746my/G4HKl6Pt9cCsYO/DwppKypHC5cB6Vd2gqieA\ncUC3HG26AWOc1xOBq0REnOXjVDVNVTcC653tBTWfqs5W1WPO7EKgWiFnOKt8+egEzFTVfaq6H5gJ\nJHmcrycwtpAz5EtV5wL78mnSDfhIfRYCsSJShSDsv4KyqeoC57Mh+N+9rAwF7b+8nM1317XTzBf0\n719hKilFoSqwxW8+xVmWaxv1jR99EKjk8r3ByOfvfnx/VWaJFJFkEVkoIjcUcrbTydfdOcUwUUSq\nn+Z7g5EP57RbbWCW3+JA7z838vodgrH/TkfO754C34rIYhHp7VGmLC1FZJmIfCMiWYMVF6n9JyLR\n+Ir6F36Li9I+LFCY1wEKieSyLOdtVXm1cfPes+X6M0TkDiARaOe3uIaqbhOROsAsEVmhqr8HOd9U\nYKyqpolIH3xHXR1cvjcY+bL0ACaqaobfskDvPze8/P65IiJX4isKbfwWt3b2XWVgpoj86vzVHGxL\n8HXTcERErgUmAfUpQvvPcT0wX1X9jyqKyj50paQcKaQA1f3mqwHb8mojImFABXyHg27eG4x8iMjV\nwLNAV1VNy1quqtucnxuAOUCTYOdT1b1+mUYBl7l9bzDy+elBjkP3IOw/N/L6HYKx/wokIpcAo4Fu\nqro3a7nfvtsFfEnhn1p1RVUPqeoR5/XXQBkRiaeI7D8/+X3/PN2Hrnl9UaMwJnxHPBvwnTbIuth0\nYY42D/PnC83jndcX8ucLzRso/AvNbvI1wXfBrH6O5XFAhPM6HlhHIV9Ic5mvit/rG4GFzuuKwEYn\nZ5zzumKw8zntGuC7qCfB3H9+n1WLvC+UduHPF5p/Dtb+c5GtBr5raa1yLI8Byvm9XgAkBWLfuch4\nbtZ/V3z/qP7h7EtX341A53PWZ/2hGePVPiyU39PrAIX4H+xa4DfnH9ZnnWUv4furGyASmOB8+X8G\n6vi991nnfWuBzh7l+w7YCSx1pinO8lbACufLvgK436N8/wRWOTlmAw393nufs1/XA/d6kc+Z/yvw\nrxzvC9b+GwtsB07i++v1fqAP0MdZL8AwJ/8KIDFY+89FttHAfr/vXrKzvI6z35Y5/+2fDcS+c5mx\nv9/3byF+BSy370aw8zlt7sF304r/+4K2DwtrsieajTHGZCsp1xSMMcYUAisKxhhjsllRMMYYk82K\ngjHGmGxWFIwxxmSzomBKPL9eUleKyASnKwJEZMEZbq9WXr1l5vOe10RklfMzQXw99f4iIlec5nZu\nEJFGp5fYGPesKJjSIFVVG6vqRcAJfPeXo6qtgpjhIaCpqg4GrgJ+VdUmqjrvNLdzA76efY0JCCsK\nprSZB9QDEJEjzs8bReQ7Z8yDKs54DOeKSKjzl/0ipyPAh/LbsPP+15wjkhUicpuzfAq+p1l/EpGn\ngFeBa52jlxgR+dDvPY8576krItOdTtTmiUhDEWkFdAVec95bN2B7yZRaJaVDPGMK5PR51RmY7r9c\nVb8Uke74ukJJAl5U1R1Oj5YHVbWZiEQA80XkW/LucO0moDFwKb4uNRaJyFxV7SoiR1S1sZNjJ74n\nmvuLyGVAVecoBhGJdbY1Et/TsutEpDkwXFU7OAXmK1WdWGg7xhg/VhRMaRAlIkud1/OA93JpMwBY\nia9Pp6wOzToCl4jIzc58BXw9c/6Wx+e0wdeTbAawU0T+CzQD8hvtbQNQR0SGAtPwdbFcFl/3HBN8\nQ34Avr65jAk4KwqmNEjN+is9H1WBTOAcEQlR1Ux8/RUNUNUZ/g3FbyjXHHLrxjlfqrpffMNfdsJ3\npHIr8ChwwEVmYwqdXVMwpZ5zWukDoBewBnjcWTUD6CsiZZx254tITD6bmgvc5lyLSMA3hOPPBXx2\nPBCiql8Az+O7GH0I2CgitzhtRP43bvJhoNyZ/J7GuGFHCsbAM8A8VZ3nnGZaJCLT8PUeWgtYIr7z\nOLvx3f2Tly+Blvh6xFTgSVXdUcBnVwU+EJGsP9D+4vy8HRghIs8BZfANM7nM+TlKRAYCN2vwBwsy\nJZz1kmqMMSabnT4yxhiTzYqCMcaYbFYUjDHGZLOiYIwxJpsVBWOMMdmsKBhjjMlmRcEYY0w2KwrG\nGGOy/X+Bbnf5RhKaVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109c45510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "offsets = np.arange(0.0,2.,0.1)\n",
    "measured_values = np.zeros(len(offsets))\n",
    "for i in range(len(offsets)):\n",
    "    offset = offsets[i]\n",
    "    matched_filter = np.exp(-(x-offset)**2/2./sig**2)/sig/np.sqrt(2.*np.pi)\n",
    "    matched_filter /= np.sum(matched_filter**2)\n",
    "    measured_values[i] = np.sum(gaussian_model*matched_filter)\n",
    "\n",
    "plt.plot(offsets,measured_values/flux)\n",
    "plt.xlabel(\"Pixel offset\")\n",
    "plt.ylabel(\"Fraction of recovered flux\")\n",
    "mf = np.exp(-(offsets)**2/2./sig**2/2.)\n",
    "plt.plot(offsets,mf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do this in 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.09560103  0.15640978  0.10811484  0.17688322]\n",
      "[ 0.17802505  0.29126107  0.20132785  0.32938603]\n",
      "0.265553988831\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD8CAYAAABXXhlaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAACxtJREFUeJzt3W/InXd9x/H3p0mamLSdLVWnSVnr\n6LqJMCr3ZrVMRqPQTTE+2KBCpRMhY2xaRZC6J33qAxF9MITQVssslRELFifOEhUZjLD0D+ufqJXW\npdHYdOvWdt3a/Ol3D+4jdLfp7nCf3/lTv+8XlPuck6u/60uS932d69wn10lVIamXcxY9gKT5M3yp\nIcOXGjJ8qSHDlxoyfKkhw5caMnypIcOXGto8z52dm621jR3z3KXUygs8z4l6MettN9fwt7GDt2f3\nPHcptXKwDpzVdj7VlxoyfKkhw5caMnypIcOXGjJ8qSHDlxoyfKkhw5caMnypIcOXGjJ8qSHDlxoy\nfKkhw5caMnypIcOXGjJ8qSHDlxoyfKkhw5caMnypoXXDT3JbkuNJHnrZYxcluSfJo5OvF852TEkj\nnc0R/8vAtWseuwk4UFWXAwcm9yW9SqwbflV9H3h6zcN7gNsnt28HPjB4LkkztNFz/DdU1TGAydfX\njxtJ0qzN/CO0kuwF9gJsY/usdyfpLGz0iP9kkjcCTL4ef6UNq2pfVa1U1coWtm5wd5JG2mj4dwM3\nTG7fAHx9zDiS5uFsfpx3J/BPwBVJjib5CPAZ4D1JHgXeM7kv6VVi3XP8qvrgK/ySn3ctvUr5zj2p\nIcOXGjJ8qSHDlxoyfKkhw5caMnypIcOXGjJ8qSHDlxoyfKkhw5caMnypIcOXGjJ8qSHDlxoyfKkh\nw5caMnypIcOXGjJ8qSHDlxqa+UdoaTGy5dxB64z5K1InToxZ59SpIet05xFfasjwpYYMX2rI8KWG\npgo/ySeSPJzkoSR3Jtk2ajBJs7Ph8JPsBD4GrFTVW4FNwHWjBpM0O9M+1d8MvCbJZmA78LPpR5I0\naxsOv6p+CnwWOAIcA56pqm+PGkzS7EzzVP9CYA9wGfAmYEeS68+w3d4kh5IcOsmLG59U0jDTPNV/\nN/B4VT1VVSeBu4B3rt2oqvZV1UpVrWxh6xS7kzTKNOEfAa5Ksj1JgN3A4TFjSZqlac7xDwL7gfuA\nBydr7Rs0l6QZmupfYFTVzcDNg2aRNCe+c09qyPClhgxfasjwpYa8As+SOef884esc+L3fmvIOs/8\n5pgr+Vz8wH8NWSf/8qMh69SLvd9M5hFfasjwpYYMX2rI8KWGDF9qyPClhgxfasjwpYYMX2rI8KWG\nDF9qyPClhgxfasjwpYYMX2rI8KWGDF9qyCvwjJCMW+rXXzdknaN/fnLIOj/8g1uHrHP53/7FmHWO\n/NqQdU4/eXzIOq9WHvGlhgxfasjwpYYMX2poqvCTvDbJ/iQ/SHI4yTtGDSZpdqZ9Vf8LwLeq6k+S\nnAtsHzCTpBnbcPhJLgDeBfwZQFWdAE6MGUvSLE3zVP/NwFPAl5Lcn+SWJDsGzSVphqYJfzPwNuCL\nVXUl8Dxw09qNkuxNcijJoZP0/tgiaVlME/5R4GhVHZzc38/qN4L/o6r2VdVKVa1sYesUu5M0yobD\nr6qfA08kuWLy0G7gkSFTSZqpaV/V/yhwx+QV/ceAD08/kqRZmyr8qnoAWBk0i6Q58Z17UkOGLzVk\n+FJDhi815BV4lkxOnR6yzon/2DZknb//7zHrbHlu0FWKTo/5/enOI77UkOFLDRm+1JDhSw0ZvtSQ\n4UsNGb7UkOFLDRm+1JDhSw0ZvtSQ4UsNGb7UkOFLDRm+1JDhSw0ZvtRQqmpuO7sgF9Xbs3tu+3s1\nOmfbmCvecMVlQ5b5n53nDVln+6P/PmSdlx4/MmSdOnVqyDrL5mAd4Nl6et3LHXnElxoyfKkhw5ca\nMnypIcOXGpo6/CSbktyf5BsjBpI0eyOO+DcChwesI2lOpgo/yS7gvcAtY8aRNA/THvE/D3wKeOmV\nNkiyN8mhJIdO8uKUu5M0wobDT/I+4HhV3fv/bVdV+6pqpapWtrB1o7uTNNA0R/yrgfcn+QnwVeCa\nJF8ZMpWkmdpw+FX16araVVWXAtcB36mq64dNJmlm/Dm+1NDmEYtU1feA741YS9LsecSXGjJ8qSHD\nlxoaco6vcV564YUxCz34oyHLbHtk05B1Tp8+PWQdXhq0TnMe8aWGDF9qyPClhgxfasjwpYYMX2rI\n8KWGDF9qyPClhgxfasjwpYYMX2rI8KWGDF9qyPClhgxfasjwpYYMX2rIS2/9qhp0iaryUle/kjzi\nSw0ZvtSQ4UsNGb7U0IbDT3JJku8mOZzk4SQ3jhxM0uxM86r+KeCTVXVfkvOBe5PcU1WPDJpN0oxs\n+IhfVceq6r7J7eeAw8DOUYNJmp0h5/hJLgWuBA6OWE/SbE39Bp4k5wFfAz5eVc+e4df3AnsBtrF9\n2t1JGmCqI36SLaxGf0dV3XWmbapqX1WtVNXKFrZOsztJg0zzqn6AW4HDVfW5cSNJmrVpjvhXAx8C\nrknywOS/Px40l6QZ2vA5flX9I5CBs0iaE9+5JzVk+FJDhi81ZPhSQ4YvNWT4UkOGLzVk+FJDhi81\nZPhSQ4YvNWT4UkOGLzVk+FJDhi81ZPhSQ4YvNWT4UkOGLzVk+FJDhi81ZPhSQ4YvNWT4UkOGLzVk\n+FJDhi81NO3HZF+b5IdJfpzkplFDSZqtaT4mexPwN8AfAW8BPpjkLaMGkzQ70xzxfx/4cVU9VlUn\ngK8Ce8aMJWmWpgl/J/DEy+4fnTwmacltnuL/zRkeq1/aKNkL7AXYxvYpdidplGmO+EeBS152fxfw\ns7UbVdW+qlqpqpUtbJ1id5JGmSb8fwYuT3JZknOB64C7x4wlaZY2/FS/qk4l+SvgH4BNwG1V9fCw\nySTNzDTn+FTVN4FvDppF0pz4zj2pIcOXGjJ8qSHDlxpK1S+952Z2O0ueAv51nc0uBv5tDuOcLedZ\n37LN1Hme36iq16230VzDPxtJDlXVyqLn+AXnWd+yzeQ86/OpvtSQ4UsNLWP4+xY9wBrOs75lm8l5\n1rF05/iSZm8Zj/iSZmxpwl+26/cluSTJd5McTvJwkhsXPROsXvIsyf1JvrEEs7w2yf4kP5j8Pr1j\nwfN8YvJn9VCSO5NsW8AMtyU5nuShlz12UZJ7kjw6+XrhvOdaaynCX9Lr950CPllVvwNcBfzlEswE\ncCNweNFDTHwB+FZV/TbwuyxwriQ7gY8BK1X1Vlb/xeh1Cxjly8C1ax67CThQVZcDByb3F2opwmcJ\nr99XVceq6r7J7edY/Uu90EuLJdkFvBe4ZZFzTGa5AHgXcCtAVZ2oqv9c7FRsBl6TZDOwnTNcGGbW\nqur7wNNrHt4D3D65fTvwgbkOdQbLEv5SX78vyaXAlcDBxU7C54FPAS8teA6ANwNPAV+anHrckmTH\nooapqp8CnwWOAMeAZ6rq24uaZ403VNUxWD2gAK9f8DxLE/5ZXb9vEZKcB3wN+HhVPbvAOd4HHK+q\nexc1wxqbgbcBX6yqK4HnWeBT2Ml58x7gMuBNwI4k1y9qnmW3LOGf1fX75i3JFlajv6Oq7lrwOFcD\n70/yE1ZPha5J8pUFznMUOFpVv3gWtJ/VbwSL8m7g8ap6qqpOAncB71zgPC/3ZJI3Aky+Hl/wPEsT\n/tJdvy9JWD1/PVxVn1vkLABV9emq2lVVl7L6+/OdqlrYEa2qfg48keSKyUO7gUcWNQ+rT/GvSrJ9\n8me3m+V5EfRu4IbJ7RuAry9wFmDKS2+NsqTX77sa+BDwYJIHJo/99eRyY1r1UeCOyTfrx4APL2qQ\nqjqYZD9wH6s/kbmfBbxjLsmdwB8CFyc5CtwMfAb4uyQfYfUb1J/Oe661fOee1NCyPNWXNEeGLzVk\n+FJDhi81ZPhSQ4YvNWT4UkOGLzX0v5oxl2flAnCGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ce30110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.special import erf\n",
    "size = 12\n",
    "# add some static offset to reflect arbitrary sampling\n",
    "_x = np.arange(size)-size//2+0.1\n",
    "_y = np.arange(size)-size//2+0.4\n",
    "_x, _y = np.meshgrid(_x, _y)\n",
    "sigma = 2./2.35\n",
    "psflet = (erf((_x + 0.5) / (np.sqrt(2) * sigma)) - \\\n",
    "    erf((_x - 0.5) / (np.sqrt(2) * sigma))) * \\\n",
    "    (erf((_y + 0.5) / (np.sqrt(2) * sigma)) - \\\n",
    "    erf((_y - 0.5) / (np.sqrt(2) * sigma)))\n",
    "\n",
    "psflet /= np.sum(psflet)\n",
    "psflet*=flux\n",
    "\n",
    "plt.imshow(psflet)\n",
    "print psflet[psflet/np.amax(psflet)>0.5]\n",
    "f = psflet[psflet/np.amax(psflet)>0.5]\n",
    "f /= np.sum(f)\n",
    "print f\n",
    "print (np.sum(f**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1157cfbd0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xd4FOX2wPHvSSGFhJpQpCX0ngAB\nAgGkCYj0Ik2aSlHatXAvevXqDxsqNqpUQaSjFOlIr0LoHUIPTYr0kvb+/tglBkhZQjabhPN5nnnY\nmXln9mSNOTtvFWMMSimlFICTowNQSimVdmhSUEopFUuTglJKqViaFJRSSsXSpKCUUiqWJgWllFKx\nNCkopZSKpUlBKaVULE0KSimlYrk4OoAn5ePjY/z8/BwdhlJKpSvbt2+/bIzxTapcuksKfn5+hIaG\nOjoMpZRKV0TklC3ltPpIKaVULE0KSimlYmlSUEopFSvdtSkopewrMjKS8PBw7t275+hQVDK4u7uT\nP39+XF1dk3W9JgWl1EPCw8Px9vbGz88PEXF0OOoJGGO4cuUK4eHh+Pv7J+sedqs+EpGJIvKXiOxL\n4LyIyDARCRORPSJS0V6xKKVsd+/ePXLmzKkJIR0SEXLmzPlUT3n2bFOYBDRK5PyLQDHr1hMYbcdY\nlFJPQBNC+vW0/+3sVn1kjFknIn6JFGkO/Gws64FuEZFsIpLXGHPeHvEc/HMZ1/atQNyz4OSeBWeP\nLDh7ZsXNMxtuXllx985OZq/sZPbOgouL1qoppZ5Njvzrlw84E2c/3HrssaQgIj2xPE1QsGDBZL3Z\n9SMbqXZmXJLlYoxwAw/uiCd3nTy55+TJfWcvIl28iMyUjeicxXDPVwbfwoHky1cQVxfnZMWjlIrf\nmTNn6NKlCxcuXMDJyYmePXsyYMAAALp168batWvJkiULd+/eJTg4mC+++IJ8+fI9dp/atWtz/vx5\n3N3d8fLyYuLEiZQoUYKFCxfy4YcfEhMTQ2RkJAMGDKBXr158/PHHjBs3Dl9fy6DfRo0aMWTIEGrX\nrs3QoUMJCgpKMvaTJ0+yadMmOnbsGO/5gQMHsnjxYho3bkyRIkXw9PSkS5cuTJo0iQYNGvDcc889\nxSeXMhyZFOJ7xjHxFTTGjAXGAgQFBcVbJinBnQcTHfU/bt+6zp2bf3Pv5t/cv32diNvXiLpznei7\nN4i5dwNz/yZy/wZOETdxibyFa9QtMkffxO3uBbLd/pssf8+DMGAtXDVenHEpxN+ZixCRozhuecuQ\nw688BQsWIqtH8lr+lXrWubi48M0331CxYkVu3rxJpUqVeOGFFyhdujQAX3/9NW3atMEYw/fff0+d\nOnXYt28fmTJleuxeU6dOJSgoiLFjxzJw4EB+/fVXevbsydatW8mfPz/379/n5MmTseXfeust3n33\n3WTHfvLkSaZNm5ZgUhgzZgyXLl3Czc3toeOTJk2ibNmyz3xSCAcKxNnPD5yz5xs6u7iQJVtOsmTL\nmbwbGMOtK+FcDNvNrfC98NdBvG6EUeTGSrxuLICTwGa4ZLKwzakgVzwLcz9bcVzylCabX3kKF8xP\n3qweKfkjKZXh5M2bl7x58wLg7e1NqVKlOHv2bGxSeEBEeOutt5g7dy5LliyhefPmCd6zVq1afP/9\n99y8eZOoqChy5rT8DXBzc6NEiRJPHGN0dDSDBg1izZo13L9/nz59+tCrVy8GDRrEwYMHCQwMpGvX\nrrz11lux1zRr1ozbt29TtWpV3nvvPQ4ePIiXl1fs1D2dOnXCw8ODzZs34+HhuL8TjkwKC4C+IjID\nqApct1d7QooRwcunAF4+BYAm/xw3hqhrZ7l8fDfXT+8h5q+D5Lp2lDJ3VuB5ewGcBbbDOZODRS4V\nuPZcLbKXa0BQqcLk8nZ31E+jVJL+7/f9HDh3I0XvWfq5LHzUtIxNZU+ePMnOnTupWrVqgmUqVqzI\noUOHEk0Kv//+O+XKlSNHjhw0a9aMQoUKUa9ePZo0aUKHDh1wcrL0ufnuu+/45ZdfAPjyyy9p2LBh\nvPebMGECWbNmZdu2bdy/f5+QkBAaNGjAkCFDGDp0KAsXLnzsmgULFuDl5cWuXbsA+PjjjwFo06YN\nI0aMsLmKyt7slhREZDpQG/ARkXDgI8AVwBjzI7AYaIylMuYO0N1esdidCC7Z85OnUn7yVHrpn+PG\nYK6f4cbpfVw7tYfoM6HUvryZzGdWEn36Y3YuLMZijyDuF6pLoXIhBBfxIZvn44/ASj2Lbt26RevW\nrfn+++/JkiVLguUsfVXi9+Dbt5+fH8OHDwdg/Pjx7N27lz/++IOhQ4eyYsUKJk2aBNhefbR8+XL2\n7NnDnDlzALh+/TpHjx6NtworvbFn76MOSZw3QB97vX+aIIJkK0jWbAXJWr6x5Vh0FNHh27m8axEF\nwlZS8eZ0nI5M48phb9bElCMsazWci9YlsGRxKvvnwMtNe0Ipx7H1G31Ki4yMpHXr1nTq1IlWrVol\nWnbnzp3Uq1cv3nMP2hQeVa5cOcqVK0fnzp3x9/ePTQq2MsYwfPjwx54k1qxZ80T3SYv0L05qc3bB\nuVBVcheqCgyG25eJOroSs3cJDU6vocWtTbDrG/bs8GeyCeBszhBylgyhWtHcVCyUHXdX7e2kMjZj\nDK+99hqlSpXi7bffTrTc8OHDOX/+PI0aJTYk6h+3bt0iNDSU2rVrA7Br1y4KFSr0xDE2bNiQ0aNH\nU7duXVxdXTly5Aj58uXD29ubmzdvPvH9knudPWhScLTMPrgEtsMnsB3ExMCF3UQeXkGhA8soe+l3\nnK7N48ZmTzZsLMtgqUhEiWa8WLEotYr74uqs8xmqjGfjxo1MmTKFcuXKERgYCMDnn39O48aWp+2B\nAwfyySefcOfOHYKDg1m9erXN1TbGGL766it69eqFh4cHmTNntukp4aWXXoqdS6hatWrMnDmTkydP\nUrFiRYwx+Pr6Mm/ePMqXL4+LiwsBAQF069btoYbmxHTr1o3evXuniYZmSaw+Li0KCgoyz8wiO3ev\nwYm1RB5eTvSRFbjfvcgd3JgfVZ2Frg3wL1+D5hXyU6lgdpycdASqShkHDx6kVKlSjg5DPYX4/huK\nyHZjTJIt2fqkkJZ5ZIPSzXEt3RxXYyA8FPftk3h57690iF7NwV2FmBZahw+86lOvQjFaVMhH8dze\njo5aKZWOaVJIL0SgQGWcClSGRl/A3tkU3z6ZTy5M4n7EdBZurMp7a+twJ3cQLSrko1ngczomQin1\nxDQppEfuWaDyazhXfg3O7cJtx2Ra7plF64h1nL5ZkEnLn6fx0hqU8C9Ei8B8vFgur46wVkrZRNsU\nMor7t2D/XNgxGcK3ESWZWOMczPg7tdghZahTMhctK+SjfqncuGgDtUqEtimkf9qmoMDNCyp2tmwX\n9+OyfTL19sygfqZ1XHErwNQTtfnv/up8ki0P3UP8aF+loI6BUEo9Rr8yZkS5y0Djr5B3DkPLseTM\nXYD+MVPY5tmfD50mMm7RRqp9vpLPFx/k3LW7jo5WKZWGaFLIyFw9IKAdvLoE+mzDqUInGt1bwubM\n7zAs+wwWbNhBra9WM2DGTvadve7oaJWK9eqrr5IrVy7Kli370PFu3brh7+9PQEAAxYsXp0uXLpw9\nezbee0RFRfH+++9TrFgxAgMDCQwM5LPPPkuN8KlevXqK3OfSpUtUrVqVChUqsH79eho3bsy1a9e4\ndu0ao0aNSpH3eJQmhWeFb3Fo+gP024FTQDvqXJ/PJo+3+SX/PHYdOEyT4RtoP3YzKw9eJCYmfbUz\nqYynW7duLF26NN5zX3/9Nbt37+bw4cNUqFCBOnXqEBER8Vi5Dz74gHPnzrF371527drF+vXriYyM\ntHfoAGzatClF7rNy5UpKlizJzp07qVmzJosXLyZbtmyaFFQKyl4Img2HfttxKteG4EtzWJPpX8wv\ntoQbl8/x2uRQXvhuLdP+PM29yGhHR6ueUbVq1SJHjhyJlnkwdXaePHlYsmTJQ+fu3LnDuHHjGD58\nOO7ulpmIvb29Y2cmBWjRogWVKlWiTJkyjB07Nva4l5dX7Os5c+bQrVs3AGbPnk3ZsmUJCAigVq1a\nAOzfv58qVaoQGBhI+fLlOXr06EP3uHXrFvXq1aNixYqUK1eO+fPnA5bZX0uVKkWPHj0oU6YMDRo0\n4O7dh6tyd+3axb///W8WL15MYGAgd+/exc/Pj8uXLzNo0CCOHTtGYGAgAwcOtPVjtYm2ND6rcvhD\ni5FQ821k3dcE7JnKIpdfCSvbno+u1OP9uXv5ZvlhXgkuROdqhfDxckv6nirjWTIILuxN2XvmKQcv\nDkmx28U3dXZYWBgFCxbE2zvhwZwTJ04kR44c3L17l8qVK9O6devYdRbiM3jwYJYtW0a+fPm4du0a\nAD/++CMDBgygU6dOREREEB398Bcpd3d35s6dS5YsWbh8+TLBwcE0a9YMgKNHjzJ9+nTGjRvHyy+/\nzK+//sorr7wSe21gYCCDBw8mNDSUESNGPHTfIUOGsG/fvthpuFOSPik863IWgZY/Qp9tSMkmFAv7\niam3erKx0jpCnhN+WHmU6kNW8d5vewj765ajo1XqMbZ0q//pp58IDAykQIECnDljWQV42LBhBAQE\nEBwczJkzZ2K/5SckJCSEbt26MW7cuNg//tWqVePzzz/nyy+/5NSpU4/NWWSM4f3336d8+fLUr1+f\ns2fPcvHiRQD8/f1j53aqVKnSQyvAOZI+KSgLn6LQehzUGois/ZJ8+8YwLNMvfBzyKiPuNmLqjrPM\n2HaGloH5eOuF4hTI4enoiFVqSMFv9PYS39TZRYsW5fTp09y8eRNvb2+6d+9O9+7dKVu2LNHR0axZ\ns4Y//viDzZs34+npSe3atbl37x5gqZZ64MExsDwV/PnnnyxatIjAwEB27dpFx44dqVq1KosWLaJh\nw4aMHz+eunXrxl4zdepULl26xPbt23F1dcXPzy/2nnGX5HR2dn6s+shR9ElBPcy3OLSZAG9uhmIv\nkGP7MP53rD07a2ylf3UfFu09T71v1vJ/v+/nyq37jo5WPcOMMQwbNizeqbM9PT157bXX6Nu3b+wf\n4ejo6NgG6evXr5M9e3Y8PT05dOgQW7Zsib02d+7cHDx4kJiYGObOnRt7/NixY1StWpXBgwfj4+PD\nmTNnOH78OIULF6Z///40a9aMPXv2PBTH9evXyZUrF66urqxevZpTp06lyM9uz6m2NSmo+OUqBW0n\nwRuboHBtPLd8y1v727L1hWO0rpCHyZtO8vzXaxi28ii370c5OlqVwXTo0IFq1apx+PBh8ufPz4QJ\nE2LPDRw4MLZL6rZt2xKcOvuzzz4jb968lC1blgoVKlCzZk26du3Kc889R6NGjYiKiqJ8+fJ8+OGH\nBAcHx143ZMgQmjRpQt26dWPXin7wvuXKlaNs2bLUqlWLgIAAZs6cSdmyZQkMDOTQoUN06dLloRg6\ndepEaGgoQUFBTJ06lZIlS6bI55MzZ05CQkIoW7Zsijc06zQXyjYX9sLyD+D4GshTnjPVP+HT3V4s\n238RHy83BtQrSvsqBXWNhwxAp7lI/55mmgu7/h8sIo1E5LCIhInIoHjOFxKRlSKyR0TWiEh+e8aj\nnkKectB5HrT5CW5fosBvzRmTdTLzXy1JYd/MfDh/P/W/XcuC3ed0nINS6ZjdkoKIOAMjgReB0kAH\nESn9SLGhwM/GmPLAYOALe8WjUoAIlG0FfbdBtb6wcyoBc+sxs9JhfupaCQ9XZ/pP30mzkRtYf/SS\no6NVSiWDPZ8UqgBhxpjjxpgIYAbQ/JEypYGV1ter4zmv0iI3b2j4GfTeALlKIwsHUGdDRxa18ebb\nlwP4+3YknSdspdP4LewJv+boaFUypLdqZfWPp/1vZ8+kkA84E2c/3Hosrt1Aa+vrloC3iCQ8ekSl\nLblLQ7dF0HIMXDuN8/g6tDr/Hav6BPC/JqU5cO4GzUZspM+0HZy4fNvR0Sobubu7c+XKFU0M6ZAx\nhitXrsSO4k4OuzU0i0hboKEx5nXrfmegijGmX5wyzwEjAH9gHZYEUcYYc/2Re/UEegIULFiwUkp1\n61Ip6O41WP05bBsHHjmgwSfcLNGGcetPMH7DCSKiYmhfpQDvvFCC7JltW2RdOUZkZCTh4eEP9dFX\n6Ye7uzv58+fH1fXhhbVsbWi2Z1KoBnxsjGlo3X8PwBgTb7uBiHgBh4wxiTY2a++jNO78blj0DoRv\ng4LVoPFQLmUuxvBVR5n652myuLvwXuNStKmYHycnSfp+SqkUkRZ6H20DiomIv4hkAtoDC+IWEBEf\nEXkQw3vARDvGo1JD3gB4dbll0r1Lh2FMLXw3fszghgVZ1L8GRXy9+PecPbQbu5nDF+wz+EYplXx2\nSwrGmCigL7AMOAjMMsbsF5HBItLMWqw2cFhEjgC5gdSZ7FzZl5MTVOwC/bZbVoLbMhpGVKbkpeXM\n6hnMV63LE/bXLV4atp4vlhzkToQOflMqrdDBa8r+wrfDorcsVUslm0DTH7hKFoYsOcis0HDyZfPg\n42ZleKF0bkdHqlSGlRaqj5SyyF8JeqyGFz6Bo8thVDA5zvzBV20CmN27GpndnOnxcyivTw4l/O87\njo5WqWeaJgWVOpycIaQ/9FwDXnlgRgeY34fKeVxY1L8m771Yko1hl3nh23X8uPYYkdExjo5YqWeS\nJgWVunKXgR4rocbbsGsa/BiC65nN9Hq+CH+88zw1ivkwZMkhXhq2nq0nrjo6WqWeOUkmBRHJFc+x\nEvYJRz0TXNyg/kfQfQmIM0x6CZb9l3yZhXFdghjXJYjb96N5ecxmBs7ezdXbj6+/q5SyD1ueFNaL\nyMsPdkTkHWBuIuWVsk3BYMtUGZW6weYRMK4OnN/DC6Vzs+LtWvR+vghzd56l7jdrmLH1tE60p1Qq\nSLL3kYjkBcYC97B0Gz0IvGOMccjajNr7KIM6ugLm94E7V6H2IAj5Fzi7cOTiTT6Yu4+tJ69SqVB2\nvmpTniK+XknfTyn1kBTrfWSMOQ8sBaoBflhmNdXFelXKKvYCvLkFSr4Eqz6Bn16EK8contubmb2C\n+bqNZWxD4x/WM3HDCX1qUMpObGlTWAFUBcoCjYHvRGSovQNTzyDPHJbV3lqNh8uH4ccasG0CArQN\nKsCKt2oRUtSHwQsP0HH8Fs5c1e6rSqU0W9oURhpjuhhjrhlj9gHVgetJXaRUsohA+bbwxmYoUAUW\nvQ1T28CN8+TK4s6ErkF81bo8+87eoNH365ix9bTO5qlUCtIRzSrtiomBbeNhxf/A1R1e+tayyA9w\n5uodBs7ZzZbjV6lTwpchrcuTO0vypwtWKqNLsTYFEbkpIjes2z0RiRYRfVJQ9ufkBFV7Qu/1kKMw\nzOkOC/pB5F0K5PBk2uvBfNS0NJuOXaHBd+tYsPucPjUo9ZRsaWj2NsZksW7uWNY8GGn/0JSy8ilm\nmXm1xtuw42cYXx8uh+HkJHQP8WfxgJr4+2Sm//Sd9J22U8c1KPUUnnhEszFmHlDXDrEolTBnF8uA\nt46z4cZZGFsb9v0GQBFfL+b0rsbAhiVYfuACDb5by4oDFx0br1LplC3jFFrF2XUCgoDnjTHV7BlY\nQrRNQXHtDMx5FcK3QuUelvWiXdwAOHj+Bm/N3MWhCzdpUyk//2tamizurkncUKmML8VWXhORn+Ls\nRgEngXHGmL+eKsJk0qSgAIiOhD8+toyEzhto6cqawx+AiKgYhq08yqg1YeTJ4s7XbQMIKerj0HCV\ncjSHL8dpL5oU1EMOLYJ5b4ABWoyCUk1iT+08/TfvzNrN8cu36VKtEINeLIlnJhfHxaqUAz11UhCR\n4Vj+V4uXMaZ/8sNLPk0K6jF/n4TZ3eDcTgjuA/U/BpdMANyNiObrZYeZuPEEfjk9+aF9BQIKZHNg\nsEo5Rkokha6JXWiMmZzM2J6KJgUVr6j7sPxD2DoG8gVZqpOyFYg9vfnYFd6ZtYtLt+7zn0YleTXE\nHycncVy8SqWylEgKU4wxnUVkgDHmhxSPMJk0KahE7Z8L8/tZeiu1HAPFG8aeunYngn/P2cPyAxep\nXcKXb9oGkNPLzYHBKpV6UmLwWiURKQS8KiLZRSRH3M3GIBqJyGERCRORQfGcLygiq0Vkp4jsEZHG\nttxXqQSVaQm91kLW/DDtZVjxEURHAZDNMxNjOldicPMybDp2hRd/WM+mY5cdHLBSaUtiSeFHLLOj\nlgS2P7Il+VVdRJyxDHJ7ESgNdBCR0o8U+wCYZYypALQHRj3pD6DUY3IWgddWWNZp2Pg9TG4KN84B\nICJ0qebHvDdD8HJ3odP4P/l2+WGidPlPpYBEkoIxZpgxphQw0RhT2BjjH2crbMO9qwBhxpjjxpgI\nYAbQ/NG3AbJYX2cFziXjZ1Dqca4e0PQHaDUOzu+GH2tC2MrY06Wfy8LCfjVoXTE/w1aF0WHcFs5e\nu+vAgJVKG2yZ5uKNZN47H3Amzn649VhcHwOviEg4sBjol8z3Uip+5V+Gnmsgsy/80hrWDQVrO5pn\nJheGtg3g+3aBHDh3g8Y/rGfZ/gsODVcpR3viaS6eQHxdOx5t1e4ATDLG5MeyVsMUEXksJhHpKSKh\nIhJ66dIlO4SqMjTf4tBjJZRrY1nAZ1YXuP/POlEtKuRjYf+aFMjhQa8p2/lo/j7uRUY7MGClHMee\nSSEcKBBnPz+PVw+9BswCMMZsBtyBx4aeGmPGGmOCjDFBvr6+dgpXZWiZMluqkhp8CocWwoQX4Orx\n2NP+Ppn59Y3qvFbDn8mbT9Fy1CaOXdIFBtWzx55JYRtQTET8RSQTlobkBY+UOQ3UAxCRUliSgj4K\nKPsQger94JVf4eZ5GFsHwv6IPe3m4syHTUozoWsQF67fpenwDczZHu7AgJVKfQkmhUfWUXhsS+rG\nxpgooC+wDDiIpZfRfhEZLCLNrMXeAXqIyG5gOtDNpLd5N1T6U6Qu9FgNWfLB1Law4fvYdgaAeqVy\ns2RALcrly8q7s3fz1sxd3Lof5cCAlUo9tkyINxi4AEzB0k7QCfA2xnxl//Aep4PXVIqJuA3z+1gG\nvJVpBc1HWKqZrKJjDCNWhfHDyiMUzOHJiI4VKZsvqwMDVir5UmzlNaChMWaUMeamMeaGMWY0loV2\nlErfMmWGNj9Z5kraPxcmNLTMo2Tl7CQMqF+M6T2CuRcZQ6tRm5iy5ZSu7qYyNFuSQrSIdBIRZxFx\nEpFOgHbNUBmDCNR4CzrNgeunLYv3HF/zUJGqhXOyZEBNqhfNyYfz9vH2rN3cidDqJJUx2ZIUOgIv\nAxetW1vrMaUyjmL1Le0MXnlgSkvYPPKhdobsmTMxsWtl3n6hOPN2naXlyE0c195JKgPS9RSUiuv+\nTcv6DAd/h3IvQ7NhltHRcaw7cokBM3YSFW34um0AjcrmcVCwStkuxdoURKS4iKwUkX3W/fIi8kFK\nBKlUmuPmDW1/hrofwN7ZMLGhZfnPOGoV9+X3fjUo7JuZ3r9s54vFB3XuJJVh2FJ9NA54D4gEMMbs\nwTLmQKmMyckJag2EjjPh6gkY+zycWP9QkfzZPZnVuxqvBBdkzLrjdBz/J3/dvOeggJVKObYkBU9j\nzNZHjmkrm8r4ijeEHqvAMyf83Bz+HPNQO4ObizOftijHty8HsCf8Gi8N28DWE1cdGLBST8+WpHBZ\nRIpgnbdIRNoA5+0alVJphU8xeH2lJUEs+Tf83h+iIh4q0qpifub1CSFzJmc6jNvC+PXHtduqSrds\nSQp9gDFASRE5C/wL6G3XqJRKS9yzQLupliqlHT9beifdvvJQkZJ5srCgXw3ql8rFp4sO0mfaDm7e\ni3RQwEolX6JJwTpjaZAxpj7gC5Q0xtQwxpxKleiUSiucnCyNz63GQ/g2GF8X/jr0UJEs7q78+Eol\n3m9ckmX7L9J85EaOXLzpoICVSp5Ek4IxJgbL/EUYY24bY/Q3XD3byreF7osh4o5lptWjKx46LSL0\nrFWEqa9X5cbdKJqP2Mj8XWcdFKxST86W6qMVIvKuiBR40jWalcqQ8gdBz9WQvZBlHejNox5qgAYI\nLpyTRf1rUDZfFgbM2MVH8/cREaXdVlXaZ8uEeCfiOWxsXJIzxengNZVmRNyG33pa1meo2BUaDwWX\nTA8ViYyO4cslhxi/4QSBBbIx+pWK5M3qkcANlbIfWwev6YhmpZ5GTAys/gzWD4VCNaDdFPB8/EF6\n8d7zDJy9G49MLox+pSKV/fRhW6WulBzR7CkiH4jIWOt+MRFpkhJBKpXuOTlBvQ8tq7qFb4NxdeHS\n4ceKNS6Xl3l9QvB2d6HD2C1M2XxSu62qNMmWNoWfgAigunU/HPjUbhEplR6Vfxm6LbJUKY2vD0f/\neKxIsdzezOsTQs1iPnw4fz//+XWPrgWt0hxbkkIR64I6D6a5uItlsR2lVFwFKltGQGcrBNPawpbR\njzVAZ/VwZULXyvSrW5RZoeG0G7uFC9d1egyVdtiSFCJExIN/RjQXAe7bNSql0qtsBeDVpVCiMSwd\nBAv/BdEPD2JzchLeaVCCH1+pyNGLN2kyfAPbTur0GCptsCUpfAwsBQqIyFRgJfBvewalVLrm5gUv\nT4Ga78D2SZYR0Hce/6PfqKylncHLzdnSzqCruqk0IMmkYIxZDrQCugHTsYxwXmPLzUWkkYgcFpEw\nERkUz/nvRGSXdTsiIteeLHyl0ignJ6j3P2g5Fs5sTbABunhub+b3rWFpZ5i3j0G/7uV+lLYzKMex\npffRAqABsMYYs9AYc9mWG4uIMzASeBEoDXQQkdJxyxhj3jLGBBpjAoHhwG9P+gMolaYFtINuCyHi\nlqUBOmzlY0Wyergyvmtl+tYpyszQM7Qbo+0MynFsqT76BqgJHBCR2SLSRkTcbbiuChBmjDlujIkA\nZgDNEynfAcuTiFIZS4EqlqU+sxWEqW1h24THijg7Ce82tLQzHNF2BuVAtlQfrTXGvAkUBsZiWa/5\nLxvunQ+Iu2RVuPXYY0SkEOAPrErgfE8RCRWR0EuXLtnw1kqlMQ8aoIvWh0Vvw9L3IebxaiJtZ1CO\nZsuTAtbeR62xTJldGZhsy2XxHEvot7s9MMcYE29lqjFmrDEmyBgT5Ovra0vISqU9bt7QYTpUfQO2\njIQZneD+rceKFc/tzfw+NagIYHyJAAAgAElEQVSh7QzKAWxpU5gJHATqYmkjKGKM6WfDvcOBAnH2\n8wPnEijbHq06Us8CJ2d4cYhlnqSjy+CnRnD98VlUs3paxjP0qVNE2xlUqrJ1RHMRY0xvY8wq63Ta\nttgGFBMRfxHJhOUP/4JHC4lICSA7sNnWoJVK96r0gI6z4epJGF8Pzu18rIizkzCwYUlGd/qnnWH7\nKW1nUPZlS1JYCfQRkTnWrZ+IuCZ1kTEmCstaDMuwPGnMMsbsF5HBItIsTtEOwAyjFafqWVOsPry2\nDJxc4KfGcHBhvMVetM6blNnNmfZjtzBr25l4yymVEmyZOns84Mo/7QidgWhjzOt2ji1eOkuqynBu\nXoQZHeDsDmjwCVTrC/J4k9y1OxH0nbaTDWGX6Vbdjw9eKoWLs03Ngkql3CypQGVjTFdr1dEqY0x3\nLI3NSqmU4J3bMple6eaw/AP4fcBjU2MAZPPMxKTulXk1xJ9Jm07SZeJW/r4d4YCAVUZmS1KIts53\nBICIFAa0K4RSKcnVA9r8ZJkaY8dk+KU13H18gL+LsxP/a1qar9uUJ/Tk3zQfuZHDF3SVXJVybEkK\nA4HVIrJGRNZiGUvwjn3DUuoZ9GBqjOaj4NQmyxrQV+Nb+BDaBhVgRq9g7kZG02rURpbvv5DKwaqM\nyqaV10TEDSiBZezBIWOMw2ZJ1TYF9Uw4ucEyjsHJGdpPg4LB8Ra7cP0evaaEsjv8Om+/UJx+dYsi\n8bRHKJWSK6/1ATyMMXuMMbsBTxF5MyWCVEolwK8GvL4S3LPB5KawZ3a8xfJkdWdmr2q0rJCPb1cc\noc+0HdyJiErlYFVGYkv1UQ9jTGzlpjHmb6CH/UJSSgHgUxRe/wPyV4bfXofVXzy2aA+Au6sz374c\nwPuNS7J03wVajdrEmat3HBCwyghsSQpOEud51Dr7aSb7haSUiuWZAzrPg4COsHYI/NYToh6vvRUR\netYqwsRulTl77S7NR25ky/ErDghYpXe2JIVlwCwRqScidbFMR7HUvmEppWK5ZIIWo6DuB7B3Fvzc\nIt5FewBql8jF/D4hZPN05ZXxf/LLllOpHKxK72xJCv/B0uPoDaAPuvKaUqlPBGoNhNYT4Ox2y9QY\nV47FW7Swrxfz+oRQs5gPH8zbx/tz9xIRZevsNOpZZ2vvIw+goDHm8aWjUpn2PlLPvNN/WkZAmxho\nNxX8QuItFh1jGLr8MKPXHKOKXw5GvVIRHy+3VA5WpRUp2fuoGbALa5WRiARaV2NTSjlCwaqWBmhP\nH/i5OeyeGW8xZyfhP41K8kP7QHaHX6P5iI3sP3c9lYNV6Y0t1UcfYVlF7RqAMWYX4GfHmJRSSclR\nGF5fYRm/MLcnrBkSb88kgOaB+ZjTuzoxxtBm9GaW7D2fysGq9MSWpBBljNGvF0qlNR7Z4ZXfLD2T\n1nwBc3vF2zMJoFz+rMzvG0KpvN68MXUH3604QkyMTkysHmdLUtgnIh0BZxEpJiLDgU12jkspZYsH\nPZPqfAB7ZibaMymXtzvTewbTplJ+flh5lDen7uD2fR3oph5mS1LoB5QB7gPTgOvAv+wZlFLqCYjA\n8w96JoXC+PoJ9kxyc3Hm6zbl+eClUiw/cIHWo3Wgm3pYor2PrAPVhhhjBqZeSInT3kdKJeL0Fpje\nATCWOZMKVU+w6Nojl+g7bQeuzk78+EolqvjnSL04VapLkd5HxphooFKKRaWUsq+CwdaeSTktPZMS\nmDMJ4PnivrED3TqN38L0radTMVCVVtlSfbRTRBaISGcRafVgs3tkSqnkyVkEXlsB+atY5kxa82WC\nPZMK+3ox980Qqhfx4b3f9vLR/H1ERutAt2eZLUkhB3AFqAs0tW5N7BmUUuopeeaAznMhoAOs+Rzm\n9k6wZ1JWD1cmdqtMj5r+TN58iq66otszzaYRzcm+uUgj4AfAGRhvjBkST5mXgY8BA+w2xnRM7J7a\npqDUEzAG1n0Nqz+DQiHQ7hdLwkjAnO3hvP/bXvJmc2dclyCK5/ZOxWCVPaXkGs3JDcAZGAm8CJQG\nOohI6UfKFAPeA0KMMWXQXk1KpSwReP7f0Go8hG9LtGcSQJtK+ZnRK5g7EdG0GrWJPw5cTMVgVVpg\nt6SAZRR0mDHmuDEmApgBNH+kTA9gpHWNBowxf9kxHqWeXeXbQpcFcPdvS2I4tTnBohULZmdB3xD8\nfTLTY0ooo9aEYc8aBZW2JJgURGSA9d/4Z9tKWj7gTJz9cOuxuIoDxUVko4hssVY3xRdLTxEJFZHQ\nS5cuJTMcpZ5xhapZeiZ5ZIefmyXaMylvVg9m965Gk/LP8dXSw/xr5i7uRUanYrDKURJ7Uuhu/Xd4\nMu8d30Kxj37dcAGKAbWBDsB4Ecn22EXGjDXGBBljgnx9fZMZjlKKnEUeXs0tkZ5J7q7ODGsfyMCG\nJViw+xwvj9nMhev3UjlgldoSSwoHReQkUEJE9sTZ9orIHhvuHQ4UiLOfHzgXT5n5xphIY8wJ4DCW\nJKGUspcHPZPKt0+yZ5KI0KdOUcZ2DuLYX7doOmIDO07/ncoBq9SUYFIwxnQAgoEw/umK+qA7alMb\n7r0NKCYi/iKSCWgPPDrl9jygDoCI+GCpTjr+hD+DUupJubhByx+hzn9hzwyY0jLBOZMAXiidm7l9\nQvBwdab9mC3MDj2TYFmVviU1ovmCMSYAOA94W7dzxpgk1/gzxkQBfbEs53kQmGWM2S8ig61rNGA9\nd0VEDgCrgYHGGF1YVqnU8IQ9k4rn9mZ+nxAq+2dn4Jw9fLLwAFE60C3DSXKcgog8D/wMnMTSTlAA\n6GqMWWf36OKh4xSUsoNTm2GGdYhQ+2mWRukEREXH8Omig0zadJKaxXwY0aEiWT1dUylQlVwpOU7h\nW6CBMeZ5Y0wtoCHw3dMGqJRKQx70TPLMkWTPJBdnJz5uVoavWpdny/ErNB+5gbC/bqZisMqebEkK\nrnHXZjbGHAH0a4FSGc0TzJkE8HLlAkzvEcyt+1G0GLmJVYd0oFtGYEtSCBWRCSJS27qNA7bbOzCl\nlAM8wZxJAEF+OVjQtwZ+Pp68NjmU0WuO6UC3dM6WpPAGsB/oDwwADgC97RmUUsqBXDJBi9E290x6\nLpsHs3tVp0n55/hy6SEd6JbO2XVCPHvQhmalUtHeOTDvDchaADrNtlQxJcAYw6g1xxi6/DBln8vK\n2C6VyJvVIxWDVYlx+IR4SqkMoFwb6Pp7nDmTEl6e/cFAt3Gdgzhx+TZNh29k+6mEnzBU2qRJQSmV\nuNjV3HJYVnPbPSPR4vVL52bum9XJ7OZMh7F/MmubDnRLTzQpKKWS9mDOpAJVYW4vWPUpxCQ8cK2Y\ndaBbFf8c/PvXPfzf7/t1oFs6kWRSEJHiIjJORJaLyKoHW2oEp5RKQzyywyu/QYXOloV7fn0VIu8m\nWDybZyYmda9M9xA/ftp4km4/bePaHV3RLa2zZUTzbuBHLN1QY7sUGGMc0i1VG5qVcjBjYNMwWPER\n5KsEHaaDV65EL5kVeoYP5u4jT1bLim4l8uiKbqktJRuao4wxo40xW40x2x9sKRCjUio9EoGQAdBu\nClzcD+PqwcUDiV7yclABpvcM5l5kNC1HbWTpvvOpFKx6UrYkhd9F5E0RySsiOR5sdo9MKZW2lWoK\nry6B6AiY0ACO/pFo8UqFsvN7vxoUz+1N71928M3yw8TEpK8u8c8CW6qPTsRz2BhjCtsnpMRp9ZFS\nacz1szC9neWp4cWvoEqPRIvfj4rmw3n7mBUaTv1SufiuXSDe7jpzjr2lWPWRMcY/ns0hCUEplQZl\nzQfdl0KxhrD4XVjyH4hJeESzm4szX7Yuz+DmZVhz+BItRm7k2KVbqRiwSowtvY9cRaS/iMyxbn1F\nRNO6Uuofbl7QfioE94E/f4Tp7eF+wjOnighdqvnxy+tVuXYnkhYjNuqEemmELW0Ko4FKwCjrVsl6\nTCml/uHkDI0+hybfQdhKmNAQriU+cC24cE4W9KtBwZyWCfVGrg7TCfUczJakUNkY09UYs8q6dQcq\n2zswpVQ6FfQqvDIHrofDuLoQnnhnxXzZPJjTuzrNAp7j62WH6TNtB7fvR6VSsOpRtiSFaBGJnQVL\nRAoTZ7yCUko9pkhdeG05uHrApMawf16ixT0yOfN9u0D+27gUS/ddoPXoTZy+cieVglVx2ZIUBgKr\nRWSNiKwFVgHv2HJzEWkkIodFJExEBsVzvpuIXBKRXdbt9ScLXymVZuUqCT1WQd4AmN0V1g1NdNEe\nEaFHrcJM6l6F89fv0WzkBjYcvZyKASuwcepsEXEDSmBZo/mQMSbhVTf+ucYZOAK8AIQD24AOxpgD\nccp0A4KMMX1tDVi7pCqVzkTegwV9Ye9sy+I9TX8AF7dELzl15TY9f97O0b9u8n7jUrxWwx8RSaWA\nM6an7pIqInWt/7YCXgKKAkWAl6zHklIFCDPGHDfGRAAzgOa2BK+UykBc3aHVOMuiPbunw+SmcOtS\nopcUypmZ396sToPSefh00UHenrVbF+5JJYlVHz1v/bdpPFsTG+6dD4jb9SDceuxRrUVkj7W7awEb\n7quUSm9E4Pl/Q9vJcH4PjKsDF/YleklmNxdGdarIOy8UZ+7Os7T9cTNnryU8AZ9KGQkmBWPMR9aX\ng40x3eNuwCc23Du+Z71H66p+B/yMMeWBP4DJ8d5IpKeIhIpI6KVLiX/DUEqlYWVaWKbGiIm2TI1x\naFGixZ2chH71ijG+SxAnL9+m6fANbArTdgZ7sqWh+dd4js2x4bpwIO43//zAubgFjDFX4rRPjMMy\nBuIxxpixxpggY0yQr6+vDW+tlEqznqtgaYD2LQEzOsGG7xJtgAbLwj3z+oaQM3MmXpnwJ2PXHdPx\nDHaSWJtCSRFpDWQVkVZxtm6Auw333gYUExF/EckEtAcWPPIeeePsNgMOPvFPoJRKf7Lkhe6LoWwr\n+ONjyzrQUYn3Xyni68XcPiE0KpuHzxcfou+0ndzS8QwpziWRcyWwtB1kw9KO8MBNIPEZrwBjTJSI\n9AWWAc7ARGPMfhEZDIQaYxYA/UWkGRAFXAW6JeunUEqlP64e0HoC+JaE1Z/B1ePQbip4JVwb4OXm\nwsiOFRm3/jhDlhziyMWb/Ni5EkV8vVIx8IzNlllSqxljNqdSPEnSLqlKZUD758Hc3pDZBzrMgDxl\nk7xkU9hl+k7fSURUDN+8HEDDMnlSIdD0KyUX2ektItni3Di7iEx8quiUUiquJ2yABqhe1Iff+9Wg\niG9mek3ZztBlh4nW9Rmemi1Jobwx5tqDHWPM30AF+4WklHomJaMBOl82D2b2qkb7ygUYsTqM7pO2\n8fdtXQf6adiSFJxEJPuDHeuqa4m1RSilVPIkowHa3dWZIa3L80Wrcmw5doWmIzaw7+z11Ik3A7Il\nKXwDbBKRT0TkE2AT8JV9w1JKPbMeNEA/wQhogA5VCjKrdzWiYwytR2/i1+3hqRBsxmPLyms/A22A\ni8BfQCtjzBR7B6aUeoYlYwQ0QGCBbPzerwYVCmbjndm7+d/8fURExaRCwBmHLU8KGGP2A7OA+cAt\nESlo16iUUgoeb4A++HuSl/h4ufHLa1XpUdOfnzefosO4LVy8cS8Vgs0YbFmOs5mIHAVOAGuBk8AS\nO8ellFIWDxqgc5WEma/A6i8gJvFv/y7OTvz3pdIM71CBg+dv0GT4BraeuJpKAadvtjwpfAIEA0eM\nMf5APWCjXaNSSqm4suSFboshsBOsHWJJDomsAf1A04DnmPtmCF5uLnQct4UJG07o9BhJsCUpRBpj\nrmDpheRkjFkNBNo5LqWUepirOzQfCY2+hCNLYXx9uHIsyctK5PFmft8Q6pTMxScLD/Dm1B3cuBeZ\nCgGnT7YkhWsi4gWsA6aKyA9YpqVQSqnUJQLBvaHzXLj1l6UB+ugfSV6Wxd2VsZ0r8d/GpVh+4CLN\nhm9g/zntthofW5JCc+AO8BawFDjGw3MhKaVU6ir8PPRcDVkLwLS2sPGHJAe6PVjuc2bPYO5FxtBy\n1CZmbD2t1UmPSDQpWJfUnG+MiTHGRBljJhtjhlmrk5RSynGy+8Fry6FUM1jxP/j1dYi4k+RlQX45\nWNS/BlX9czDot728M3s3dyK08uOBRJOCMSYauCMiWVMpHqWUsl2mzNB2EtT7H+z7FSY2hGtnkrws\np5cbk7pX4V/1izF351lajNxI2F+37B9vOmBL9dE9YK+ITBCRYQ82ewemlFI2EYGa70DHmfD3SRhb\nG04m3UHS2Un4V/3i/PxqFS7fiqDZiA3M33XW7uGmdbYkhUXAh1gamrfH2ZRSKu0o3tAynsEjO/zc\nDLaOS7KdAaBmMV8W969J6bxZGDBjFx/M28v9qOhUCDhtSnA9BREpaIw5ncrxJEnXU1BKJeredfi1\nBxxdBhU6w0vfgItbkpdFRsfw9bLDjF13nPL5szKyY0UK5PBMhYBTR0qspzAvzs3iW6dZKaXSHves\n0GE61HwXdk6BSU3g5oUkL3N1duL9xqUY07kSJy7f5qVh6/njwMVUCDhtSSwpSJzXhe0diFJKpRgn\nZ6j3oaUR+uI+SztDuG01DA3L5GFRv5oUyOHJ6z+H8sWSg0RFPzuT6iWWFEwCr5VSKn0o09LSbdXZ\nFX56EXb8bNNlBXN68usb1elYtSBj1h6n47g/n5lJ9RJLCgEickNEbgLlra9viMhNEblhy81FpJGI\nHBaRMBEZlEi5NiJiRCTJ+i6llHoiecpBz7VQqDos6Afz+0Dk3SQvc3d15vOW5fi+XSB7z17npWHr\nWX806XUd0rsEk4IxxtkYk8UY422McbG+frCfJakbWwe+jQReBEoDHUSkdDzlvIH+wJ/J/zGUUioR\nnjngld+g1kDY+QtMeAGunrDp0hYV8rGgbwjZPTPRecJWvlhyMEOv0WDTegrJVAUIM8YcN8ZEADOw\nTJnxqE+wrOT2bDybKaUcw8kZ6n4AHWfBtdMw9nk4vNSmS4vl9mZB3xqx1Ultx2zm9JWkR0+nR/ZM\nCvmAuEMLw63HYolIBaCAMWahHeNQSql/FG8IvdZBtkIwvR2s/MSyiE8SPDJZqpNGdarIiUu3aDxs\nfYYc7GbPpCDxHIttsBYRJ+A74J0kbyTSU0RCRST00qWMX6enlLKzB/MmVegM64fClJZw+7JNlzYu\nl5fFA2pSMo83A2bsYuDs3dy+n3HmTrJnUggHCsTZzw+ci7PvDZQF1ojISSwL+SyIr7HZGDPWGBNk\njAny9fW1Y8hKqWeGqwc0HwHNhsPpLTCmFpzZZtOl+bN7MqNnMP3rFmXOjnCaDt/AvrMZYypueyaF\nbUAxEfEXkUxAe2DBg5PGmOvGGB9jjJ8xxg/YAjQzxuhwZaVU6qnYBV5fAU4ulm6rNk6P4eLsxNsN\nSjDt9WBuR0TRatQmJmaAld3slhSMMVFAX2AZcBCYZYzZLyKDRaSZvd5XKaWeWN4A6LUWitSFxe/C\nbz0g4rZNl1YrkpMlA2pRq7gPgxce4LXJoVy5dd/OAdtPgnMfpVU695FSym5iYmDDt7D6M/ApAe2m\ngE8xmy41xjB500k+X3yIbJ6ufN8ukOpFfewcsO1SYu4jpZR6tjg5Qa13LWMabv8FY+vAgfk2XSoi\ndAvxZ26f6ni5u9Bpwp98vewQkelsigxNCkop9agidSzdVn1LwKwusOy/EB1p06VlnsvKwn41aFsp\nPyNXH6PdmM2cuZp+xjRoUlBKqfhkzQ/dl0CVnrB5BExuBjfO23SpZyYXvmoTwA/tAzl60TKmYdEe\n2651NE0KSimVEJdM0PhraDUezu+CH0PgyDKbL28emI9F/WtSxNeLPtN28O7s3dy8Z9sTh6NoUlBK\nqaSUb2uZVM87L0x7GZa+D1G29TAqmNOT2b2r0bdOUX7bEU6j79ez+dgVOwecfJoUlFLKFr7F4fWV\nluqkLSMtk+pdOWbTpa7OTrzbsASze1fH1VnoOH4Lny48wL3ItLfspyYFpZSylau7pTqp/TT4+5Rl\nFPTumTZfXqlQdhYPqEmnqgUZv+FEmhwJrUlBKaWeVMmX4I2NkKc8zO0Jc3vD/Vs2XeqZyYVPW5Rj\nUvfKXL8bSYuRGxmx6miaWd1Nk4JSSiVH1vzQ9Xd4/j+wZ6blqeHcLpsvr10iF8v+VYuGZfMwdPkR\n2o7ZzInLto2itidNCkoplVzOLlDnfUtyiLxraWfYMtqmuZMAsmfOxMiOFfmhfSDH/rpF4x/WM2XL\nKYfOn6RJQSmlnpZfDei9AYrUg6WDYHp7uG17D6PmgflY/tbzBPll58N5++j60zaHrQmtSUEppVJC\n5pzQYTo0+hKOrbKMaTix3ubL82R15+dXqzC4eRm2nrhCg+/W8fvuc0lfmMI0KSilVEoRgeDe8Pof\nkCkzTG4Kqz6DaNsW4RERulTzY3H/mvj5ZKbf9J30n76Ta3ci7Bz4PzQpKKVUSssbYBnsFtAB1n0F\nk5vA9XCbLy/s68Wvvavx9gvFWbz3PA2/X8e6I6mz6qQmBaWUsgc3L2g5GlqOhQt7YXQI7J9r8+Uu\nzk70r1eMuW+G4O3uSpeJW5my+aTdwn1Ak4JSStlTQDvLjKs5CsPsbjDnVbhz1ebLy+W3zLra6/nC\n1CuV235xWmlSUEope8tZBF5bAXX+a1mfYVQ1OLLc5svdXZ1578VSPJfNw45BWmhSUEqp1ODsAs//\nG3qsAs8cMK0tLOgH9286OrKHaFJQSqnUlDcAeq6BkH/Bzl9gdPUn6rpqb3ZNCiLSSEQOi0iYiAyK\n53xvEdkrIrtEZIOIlLZnPEoplSa4uMEL/wfdl4KTi6V30tL3LKOiHcxuSUFEnIGRwItAaaBDPH/0\npxljyhljAoGvgG/tFY9SSqU5BataRkJX7gFbRsGPNSF8u0NDsueTQhUgzBhz3BgTAcwAmsctYIy5\nEWc3M+C4CT+UUsoRMmWGl4ZC53n/zJ+06lOISr0Ba3HZMynkA87E2Q+3HnuIiPQRkWNYnhT62zEe\npZRKu4rUgTc3Qfl2sO5rGF8XLu5P9TDsmRQknmOPPQkYY0YaY4oA/wE+iPdGIj1FJFREQi9dSp1R\nfUoplercs1oGvLWfBjcvwNjasOE7iEm9FdrsmRTCgQJx9vMDic3uNANoEd8JY8xYY0yQMSbI19c3\nBUNUSqk0qORL8OYWKN4I/vgYJjayeenPp2XPpLANKCYi/iKSCWgPLIhbQESKxdl9CThqx3iUUir9\nyOwDL/8MrcbD5cPwYw3Y95vd39bFXjc2xkSJSF9gGeAMTDTG7BeRwUCoMWYB0FdE6gORwN9AV3vF\no5RS6Y4IlG8LfiGw6B3IWdT+b+nIFX6SIygoyISGhjo6DKWUSldEZLsxJiipcjqiWSmlVCxNCkop\npWJpUlBKKRVLk4JSSqlYmhSUUkrF0qSglFIqliYFpZRSsTQpKKWUipXuBq+JyCXgVDIv9wEup2A4\nKU3jezoa39NL6zFqfMlXyBiT5ORx6S4pPA0RCbVlRJ+jaHxPR+N7emk9Ro3P/rT6SCmlVCxNCkop\npWI9a0lhrKMDSILG93Q0vqeX1mPU+OzsmWpTUEoplbhn7UlBKaVUIjJMUhCRRiJyWETCRGRQPOfd\nRGSm9fyfIuIX59x71uOHRaShg+J7W0QOiMgeEVkpIoXinIsWkV3WbcGj16ZSfN1E5FKcOF6Pc66r\niBy1bnZZKMmG+L6LE9sREbkW51xqfH4TReQvEdmXwHkRkWHW+PeISMU45+z6+dkQWydrTHtEZJOI\nBMQ5d1JE9lo/O7stZGJDjLVF5Hqc/47/i3Mu0d+NVIpvYJzY9ll/53JYz6XKZ5hijDHpfsOystsx\noDCQCdgNlH6kzJvAj9bX7YGZ1telreXdAH/rfZwdEF8dwNP6+o0H8Vn3b6WBz68bMCKea3MAx63/\nZre+zp7a8T1Svh+Wlf5S5fOzvkctoCKwL4HzjYElgADBwJ+p+PklFVv1B+8JvPggNuv+ScAnDXx+\ntYGFT/u7Ya/4HinbFFiV2p9hSm0Z5UmhChBmjDlujIkAZgDNHynTHJhsfT0HqCciYj0+wxhz3xhz\nAgiz3i9V4zPGrDbG3LHubgHyp3AMTxVfIhoCK4wxV40xfwMrgEYOjq8DMD2FY0iUMWYdcDWRIs2B\nn43FFiCbiOQlFT6/pGIzxmyyvjek/u/egxiS+vwS8jS/uzZ7wvhS/fcvJWWUpJAPOBNnP9x6LN4y\nxpgo4DqQ08ZrUyO+uF7D8q3yAXcRCRWRLSLSIoVje5L4WlurGOaISIEnvDY14sNa7eYPrIpz2N6f\nny0S+hlS4/N7Eo/+7hlguYhsF5GeDorpgWoisltElohIGeuxNPX5iYgnlqT+a5zDaekzTJKLowNI\nIRLPsUe7VSVUxpZrn5bN7yEirwBBwPNxDhc0xpwTkcLAKhHZa4w5lsrx/Q5MN8bcF5HeWJ666tp4\nbWrE90B7YI4xJjrOMXt/frZw5O+fTUSkDpakUCPO4RDrZ5cLWCEih6zfmlPbDizTNNwSkcbAPOD/\n27ujEKmqOI7j358lFpuYoqRYIZolGaaZREu9VJQamZZl1Evlg0YaFVTUFkUvRT5GBWnYW5CFIBka\nQuWimEYpq1QaK0SYW4hlglTkv4dzZrouzu5VZu/I7u8Dw9y5c8+d/5y5zLnnnpn/mco5VH/ZXcC2\niCj2Ks6VOixlsPQUfgYuKzy+FDjUaBtJ5wOjSN3BMmWriA9JtwEdwIKI+Ku2PiIO5ftu4AtgVtXx\nRcSRQkyrgdlly1YRX8ED9Oq6V1B/ZTR6D1XUX78kzQDWAHdHxJHa+kLd/Qqsp/mXVkuJiGMRcTwv\nfwoMlzSWc6T+Cvo6/lpah6W1elCjGTdSj6ebdNmgNtg0vdc2j3PqQPOHeXk6pw40d9P8geYy8c0i\nDZhN7bV+NDAiL48FDtDkgbSS8U0oLC8CduTlMcDBHOfovDym6vjydleRBvVUZf0VXmsSjQdK7+TU\ngeadVdVfidguJ42ltTXKugcAAAN1SURBVPda3waMLCxvB+YORN2ViHF87XMlfan+lOuy1LEx0PHl\n52snmm2tqsOmvM9WB9DED2w+sD9/sXbkda+SzroBLgDW5YN/JzC5ULYjl/sBmNei+LYAPcDufNuQ\n17cDXflg7wKWtii+14B9OY7PgWmFso/mev0ReKQV8eXHrwCv9ypXVf19APwC/EM6e10KLAeW5+cF\nvJXj7wKur6r+SsS2BjhaOPa+zusn53rbkz/7joGou5IxrigcfzsoNGCnOzaqji9v8zDpRyvFcpXV\nYbNu/kezmZnVDZYxBTMzawI3CmZmVudGwczM6twomJlZnRsFMzOrc6Ngg14hS+peSetyKgIkbT/L\n/U1qlC2zjzKrJO3L9+OUMvV+K+nmM9zPQklXn1nEZuW5UbCh4EREzIyIa4C/Sb8vJyLaK4xhGXBd\nRDwD3Ap8HxGzIqLzDPezkJTZ12xAuFGwoaYTuAJA0vF8v0jSljznwYQ8H8N4SeflM/tdORHgsr52\nnMuvyj2SLklL8voNpH+zfiXpOeANYH7uvbRJer9Q5qlcZoqkTTmJWqekaZLagQXAqlx2yoDVkg1Z\ngyUhnlm/cs6recCm4vqIWC/pXlIqlLnAyxFxOGe0/CMi5kgaAWyT9BmNE67dA8wEriWl1NglaWtE\nLJB0PCJm5jh6SP9oXiFpNjAx92KQdHHe17ukf8sekHQD8HZE3JIbmE8i4qOmVYxZgRsFGwoulLQ7\nL3cC751mm5XAXlJOp1pCs9uBGZIW58ejSJk59zd4nZtImWT/BXokfQnMAfqa7a0bmCzpTWAjKcXy\nRaT0HOvSlB9Ays1lNuDcKNhQcKJ2lt6HicBJ4BJJwyLiJClf0cqI2FzcUIWpXHs5XRrnPkXEUaXp\nL+8g9VTuB54Efi8Rs1nTeUzBhrx8WWkt8CDwHfB0fmoz8Jik4Xm7KyW19bGrrcCSPBYxjjSF485+\nXnssMCwiPgZeIg1GHwMOSrovbyP9P2/yn8DIs3mfZmW4p2AGLwCdEdGZLzPtkrSRlD10EvCN0nWc\n30i//mlkPXAjKSNmAM9GxOF+XnsisFZS7QTt+Xz/EPCOpBeB4aRpJvfk+9WSngAWR/WTBdkg5yyp\nZmZW58tHZmZW50bBzMzq3CiYmVmdGwUzM6tzo2BmZnVuFMzMrM6NgpmZ1blRMDOzuv8AVY1mVfAt\ncnUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1157cf150>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "offsets = np.arange(0.0,2.,0.1)\n",
    "measured_values_psflet = np.zeros(len(offsets))\n",
    "for i in range(len(offsets)):\n",
    "    offset = offsets[i]\n",
    "    matched_filter = (erf((_x + 0.5-offset) / (np.sqrt(2) * sigma)) - \\\n",
    "                            erf((_x - 0.5-offset) / (np.sqrt(2) * sigma))) * \\\n",
    "                            (erf((_y + 0.5) / (np.sqrt(2) * sigma)) - \\\n",
    "                            erf((_y - 0.5) / (np.sqrt(2) * sigma)))\n",
    "    matched_filter /= np.sum(matched_filter)\n",
    "    matched_filter /= np.sum(matched_filter**2)\n",
    "    measured_values_psflet[i] = np.sum(psflet*matched_filter)\n",
    "    \n",
    "plt.plot(offsets,measured_values_psflet/flux)\n",
    "plt.xlabel(\"Pixel offset\")\n",
    "plt.ylabel(\"Fraction of recovered flux\")\n",
    "plt.plot(offsets,measured_values)\n",
    "plt.legend([\"2D PSFLet fit\",\"1D Gaussian fit\"])\n",
    "# mf = np.exp(-(offsets)**2/2./sig**2/2.)\n",
    "#plt.plot(offsets,mf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Study the number of effective detector pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline --no-import-all\n",
    "matplotlib.rcParams['image.origin'] = 'lower'\n",
    "matplotlib.rcParams['image.interpolation'] = 'nearest'\n",
    "from scipy.special import erf\n",
    "def psf(fwhm=2.,dx=0.,dy=0.,sizex=10,sizey=10):\n",
    "    _x = np.arange(sizex)-sizex//2+dx\n",
    "    _y = np.arange(sizey)-sizey//2+dy\n",
    "    _x, _y = np.meshgrid(_x, _y)\n",
    "    sigma = fwhm/2.35\n",
    "    psflet = (erf((_x + 0.5) / (np.sqrt(2) * sigma)) - \\\n",
    "    erf((_x - 0.5) / (np.sqrt(2) * sigma))) * \\\n",
    "    (erf((_y + 0.5) / (np.sqrt(2) * sigma)) - \\\n",
    "    erf((_y - 0.5) / (np.sqrt(2) * sigma)))\n",
    "\n",
    "    return psflet/np.sum(psflet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1158f1c10>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAADKCAYAAACFWKrDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAADFpJREFUeJzt3X2sZPVdx/H3x91l16VbHhUpEAFj\nGtFEITeVFkua0iDQBtQ0hsYqtiabRlEwGsU0aRv/qw+ND2lqVopWJS2RUiUttZA+hPoHq8t2eepi\n2VJst2wXFLu71rgsy9c/5my9XGbuPduZM/f+yvuVTO7cOb8588lvz37umTMPJ1WFJKkd37PaASRJ\nx8filqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDVm/RArPSEbaxMnDrFqSfqu9L98\ni2frcPqMHaS4N3EiP5nLhli1JH1X2l6f7j3WQyWS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3\nJDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtS\nYyxuSWpMr+JO8ptJHknycJIPJ9k0dDBJ0ngrFneSs4DfABaq6seAdcC1QweTJI3X91DJeuB7k6wH\nNgNPDhdJkrScFYu7qr4O/BHwVWAfcKCq7h46mCRpvD6HSk4BrgHOA14BnJjkrWPGbU2yI8mOIxye\nfVJJEtDvUMkbgK9U1dNVdQS4A3jN0kFVta2qFqpqYQMbZ51TktTpU9xfBS5OsjlJgMuA3cPGkiRN\n0ucY93bgdmAn8FB3n20D55IkTbC+z6Cqejfw7oGzSJJ68JOTktQYi1uSGmNxS1Jjeh3jlmYimXoV\n6056+XQRtmyZOkMdOjT1Oo4eODj1Oqiafh1qknvcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEW\ntyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbv49bcTPtd2gDPvPFHprr/\n/tc+P3WGMz4//f7OqZ/YPfU6jn7zwNTrUJvc45akxljcktQYi1uSGmNxS1JjehV3kpOT3J7k0SS7\nk7x66GCSpPH6vqvkT4F/qqo3JzkB2DxgJknSMlYs7iQvBy4Ffhmgqp4Fnh02liRpkj6HSs4Hngb+\nKskXktyc5MSBc0mSJuhT3OuBi4APVNWFwLeAm5YOSrI1yY4kO45weMYxJUnH9CnuvcDeqtre/X47\noyJ/garaVlULVbWwgY2zzChJWmTF4q6qbwBfS/LK7qbLgC8OmkqSNFHfd5X8OnBr946Sx4G3DRdJ\nkrScXsVdVbuAhYGzSJJ68JOTktQYi1uSGmNxS1JjPJGC5iZbtky9jmlPhPCVq7dNneE8tk69jtPu\nnX4u8EQKL1nucUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUt\nSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGuOJFDQ3dejQ1Os44/PT7WvM4iQI02aA2cyFXrrc45ak\nxljcktQYi1uSGmNxS1Jjehd3knVJvpDk40MGkiQt73j2uG8Adg8VRJLUT6/iTnI28Ebg5mHjSJJW\n0neP+0+A3wGenzQgydYkO5LsOMLhmYSTJL3YisWd5E3AU1V1/3LjqmpbVS1U1cIGNs4soCTphfrs\ncV8CXJ3kCeAjwOuT/N2gqSRJE61Y3FX1e1V1dlWdC1wLfKaq3jp4MknSWL6PW5Iac1xfMlVVnwM+\nN0gSSVIv7nFLUmMsbklqjMUtSY3xRAqam6MHDk69jlM/Md23Lpx275apM8ziJAizmAu9dLnHLUmN\nsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiL\nW5IaY3FLUmMsbklqjCdS0PxUTb2Ko988MN0Kpr2/tAa4xy1JjbG4JakxFrckNWbF4k5yTpLPJtmd\n5JEkN8wjmCRpvD4vTj4H/FZV7UyyBbg/yT1V9cWBs0mSxlhxj7uq9lXVzu76IWA3cNbQwSRJ4x3X\nMe4k5wIXAtuHCCNJWlnv93EneRnwUeDGqjo4ZvlWYCvAJjbPLKAk6YV67XEn2cCotG+tqjvGjamq\nbVW1UFULG9g4y4ySpEX6vKskwAeB3VX1vuEjSZKW02eP+xLgF4HXJ9nVXa4aOJckaYIVj3FX1T8D\nmUMWSVIPfnJSkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklq\njMUtSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY\n3JLUGItbkhrTq7iTXJHk35LsSXLT0KEkSZOtWNxJ1gHvB64ELgDekuSCoYNJksbrs8f9KmBPVT1e\nVc8CHwGuGTaWJGmSPsV9FvC1Rb/v7W57gSRbk+xIsuMIh2eVT5K0RJ/izpjb6kU3VG2rqoWqWtjA\nxumTSZLG6lPce4FzFv1+NvDkMHEkSSvpU9z/CvxwkvOSnABcC9w5bCxJ0iTrVxpQVc8luR74FLAO\nuKWqHhk8mSRprBWLG6Cq7gLuGjiLJKkHPzkpSY2xuCWpMRa3JDUmVS96S/b0K02eBv59mSGnA/8x\n8weePXPOTgsZwZyzZs7+frCqvq/PwEGKe8UHTXZU1cLcH/g4mXN2WsgI5pw1cw7DQyWS1BiLW5Ia\ns1rFvW2VHvd4mXN2WsgI5pw1cw5gVY5xS5K+cx4qkaTGDFrcK53yLMnGJLd1y7cnOXfIPBMynpPk\ns0l2J3kkyQ1jxrwuyYEku7rLu+ads8vxRJKHugw7xixPkj/r5vPBJBfNOd8rF83RriQHk9y4ZMyq\nzGWSW5I8leThRbedmuSeJI91P0+ZcN/rujGPJbluFXL+YZJHu3/TjyU5ecJ9l90+5pDzPUm+vujf\n9qoJ953bqRAn5LxtUcYnkuyacN+5zedxq6pBLoy+kOrLwPnACcADwAVLxvwq8Bfd9WuB24bKs0zO\nM4GLuutbgC+Nyfk64OPzzjYm6xPA6cssvwr4JKPvUL8Y2L6KWdcB32D03tRVn0vgUuAi4OFFt/0B\ncFN3/SbgvWPudyrwePfzlO76KXPOeTmwvrv+3nE5+2wfc8j5HuC3e2wXy/bC0DmXLP9j4F2rPZ/H\nexlyj7vPKc+uAT7UXb8duCzJuBM3DKaq9lXVzu76IWA3Y87w04hrgL+pkfuAk5OcuUpZLgO+XFXL\nfRBrbqrqXuCZJTcv3v4+BPzMmLv+NHBPVT1TVf8F3ANcMc+cVXV3VT3X/Xofo+/EX1UT5rOPuZ4K\ncbmcXdf8PPDhoR5/KEMWd59Tnn17TLdhHgBOGzDTsrpDNRcC28csfnWSB5J8MsmPzjXY/yvg7iT3\nJ9k6Znmv08zNybVM/g+xFuYS4Iyq2gejP+DA948Zs5bmFODtjJ5VjbPS9jEP13eHdG6ZcOhpLc3n\na4H9VfXYhOVrYT7HGrK4+5zyrNdp0eYhycuAjwI3VtXBJYt3MnrK/+PAnwP/MO98nUuq6iLgSuDX\nkly6ZPmamM/uhBtXA38/ZvFamcu+1sScAiR5J/AccOuEISttH0P7APBDwE8A+xgdhlhqzcwn8BaW\n39te7fmcaMji7nPKs2+PSbIeOInv7OnXVJJsYFTat1bVHUuXV9XBqvrv7vpdwIYkp885JlX1ZPfz\nKeBjjJ52LrZWTjN3JbCzqvYvXbBW5rKz/9ihpO7nU2PGrIk57V4UfRPwC9UdgF2qx/YxqKraX1VH\nq+p54C8nPP5amc/1wM8Bt00as9rzuZwhi7vPKc/uBI69Sv9m4DOTNsqhdMe5Pgjsrqr3TRjzA8eO\nvSd5FaN5+8/5pYQkJybZcuw6oxesHl4y7E7gl7p3l1wMHDh2KGDOJu7JrIW5XGTx9ncd8I9jxnwK\nuDzJKd1T/8u72+YmyRXA7wJXV9X/TBjTZ/sY1JLXU352wuOvlVMhvgF4tKr2jlu4FuZzWUO+8sno\nXQ5fYvQq8ju7236f0QYIsInR0+k9wL8A58/71Vngpxg9VXsQ2NVdrgLeAbyjG3M98AijV8DvA16z\nCjnP7x7/gS7LsflcnDPA+7v5fghYWIWcmxkV8UmLblv1uWT0h2QfcITRXt+vMHo95dPAY93PU7ux\nC8DNi+779m4b3QO8bRVy7mF0XPjY9nnsnVivAO5abvuYc86/7ba7BxmV8ZlLc3a/v6gX5pmzu/2v\nj22Ti8au2nwe78VPTkpSY/zkpCQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4Jakx/weL\nVswuSTAzfQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cddd390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(psf(sizey=10,sizex=20))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construct microspectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x115a08d50>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADxCAYAAAD1LG0eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGPNJREFUeJzt3X+0HOV93/H3BwmEkQ3GKHWxJAel\nqDmVE4JtWUBdU2JiLNwUpa1IhJNYOOQoJKa1e2In0KbYkU96SpOYtDV1rBpcAk6AyHZ7Ty1bxhDH\naWoU/TABC5lwrVC4QEyEVDAmAl3dT/+Yucqy7I/Z1e7dHfF5+czxzswzM1/N5X73uc88zzOyTURE\nzL3jRh1ARMTLVRJwRMSIJAFHRIxIEnBExIgkAUdEjEgScETEiCQBR0SMSBJwRMSIJAFHRIzI/FEH\nEBFxNN75owv91P7DlcruvO/5rbZXDzmkypKAI6LWntp/mD/b+vpKZeed/tCiIYfTkyTgiKg1AzPM\njDqMvqQNOCJqzZhDPlxpqULSakkPSpqUdHWL/edL2iVpWtLaFvtPlvSYpI93u1YScETU3kzF/3Uj\naR5wA3AxsAK4TNKKpmKPAJcDv9/mNB8F/rhK3GmCiIhaM+bw4KbVXQVM2t4LIOk2YA3wwJHr2Q+X\n+16S0SW9GXgt8CVgZbeLpQYcEbU3gystFSwGHm1Ynyq3dSXpOOC3gQ9VjTs14IioNQOHqyVXgEWS\ndjSsb7K9qWFdbS5RxS8BW2w/KrU6zUslAUdE7VWs3QLss92paWAKWNqwvgR4vOK5zwPeJumXgFcC\nJ0h61vZLHuTNSgKOiFozcGhwbcDbgeWSlgGPAeuAd1eKw/7p2c+SLgdWdkq+kDbgiKg5Yw5XXLqe\ny54GrgK2AnuAO2zvlrRR0iUAkt4iaQq4FPikpN39xq68lDMi6uyss473xJZqA9yWLf2rnV2aIOZU\nmiAiotaKkXD1lAQcETUnDrfsvDD+koAjotaKh3BJwBERc67oB5wEHBExEjOpAUdEzL3UgCMiRsSI\nwzUd0pAEHBG1lyaIiIgRMOIFzxt1GH1JAo6IWisGYqQJIiJiJPIQLiJiBGxx2KkBR0SMxExqwBER\nc694CFfPVFbPqCMiSnkIFxExQofTDzgiYu5lJFxExAjNpBdERMTcKybjSQKOiJhzRhzKUOSIiLln\nU9uBGPWMOiLiCDFTcal0Nmm1pAclTUq6usX+8yXtkjQtaW3D9rMlfV3Sbkn3SfqpbtdKDTgias0M\nrgYsaR5wA/AOYArYLmnC9gMNxR4BLgc+2HT4c8B7bD8k6XXATklbbf+/dtdLAo6I2hvgQ7hVwKTt\nvQCSbgPWAEcSsO2Hy30zjQfa/ouGz49LehL4PiAJOCKOTUaDnJB9MfBow/oUcE6vJ5G0CjgB+Han\ncknAEVFrxWvpK6eyRZJ2NKxvsr2pYb1VJncv8Ug6HbgFWG97plPZJOCIqDn1Mh/wPtsrO+yfApY2\nrC8BHq8ciXQy8AXg12zf0618EnBE1JoZ6Ei47cByScuAx4B1wLurHCjpBODzwO/Z/sMqx6QbWkTU\n3uGyFtxt6cb2NHAVsBXYA9xhe7ekjZIuAZD0FklTwKXAJyXtLg//SeB84HJJ95bL2Z2ulxpwRNSa\nrYHOBWF7C7Cladu1DZ+3UzRNNB93K3BrL9caSgI+QQt8IguHceqIOIYc5Hu84OePqgtD8RAuQ5GP\nOJGFnKMLh3HqiDiGbPNdAzhL3gkXETESxUO4TMgeETESmY4yImIEBjwSbk4lAUdE7eWlnBERI2DD\noZkk4IiIOVc0QSQBR0SMRA9zQYyVJOCIqLV0Q4uIGJk0QUREjEzV972NmyTgiKi1ohdE5oKIiJhz\nGYgRETFCaYKIiBiB9IKIiBihY7oXhKR/Dfw8xZfN/cB7bR8cZmAREVXYYrqmCbhr1JIWA/8KWGn7\nh4B5FC+qi4gYCzNWpWXcVG2CmA+8QtIh4CR6eE1zRMQw1bkNuGsN2PZjwG8BjwBPAE/b/nJzOUkb\nJO2QtOMQzw8+0oiINupaA67SBHEqsAZYBrwOWCjpZ5rL2d5ke6XtlcezYPCRRkS0MNsP+JhMwMCP\nAX9p+69tHwI+B/zD4YYVEVHdDKq0VCFptaQHJU1KurrF/vMl7ZI0LWlt0771kh4ql/XdrlWlDfgR\n4FxJJwF/A1wI7Kj0L4mIGDIbpgc0IbukecANwDuAKWC7pAnbDzQUewS4HPhg07GvAT4MrKRomt5Z\nHnug3fWqtAFvAzYDuyi6oB0HbOrh3xQRMVQDbIJYBUza3mv7BeA2iibYI2w/bPs+YKbp2HcCd9re\nXybdO4HVnS5WqReE7Q9TZPaIiLHS41wQiyQ1/gW/yXZjhXIx8GjD+hRwTsVztzp2cacDMhIuImrP\n1RPwPtsrO+xvdSJXPHfPxyYBv9ypvyfDxy2Ym54uM8/32aXRVX9n4lgwwMl4poClDetLqD7uYQq4\noOnYr3Y6oJ7j9yIiSvZA24C3A8slLZN0AsWo34mKoWwFLpJ0atl996JyW1tJwBFRc+LwzHGVlm5s\nTwNXUSTOPcAdtndL2ijpEgBJb5E0BVwKfFLS7vLY/cBHKZL4dmBjua2tNEFERO310AZc4VzeAmxp\n2nZtw+ftFM0LrY69Cbip6rWSgCOi1uo8F0QScETUm+v7zDUJOCJqL68kiogYAZcP4eooCTgiai9N\nEBERIzLIXhBzKQk4ImrNTgKOiBiZdEOLiBiRtAHH6PUxsc78JR1ny2vrqfNbDgQauNO+NtXXcdNT\nj/V+UF1/i1/mjJhJL4iIiNGo61dnEnBE1FuNH8JVqrdLerWkzZK+JWmPpPOGHVhERGWuuIyZqjXg\n/wR8yfbaco7Mk4YYU0RET+paA+6agCWdDJxP8RZQyhfVvTDcsCIiqjEwM1PPBFylCeIHgL8GPi3p\nG5I+JWnhkOOKiKjGgFVtGTNVEvB84E3AJ2y/EfgecHVzIUkbJO2QtOMQfb7HKyKiD3a1ZdxUScBT\nwJTtbeX6ZoqE/CK2N9leaXvl8czNCxsjIoDaPoTrmoBt/xXwqKQfLDddCDww1KgiIioTdrVl3FTt\nBfEvgc+UPSD2Au8dXkgRET0aw9ptFZUSsO17gZVDjiUioncGH8O9ICIixpwqLhXOJK2W9KCkSUmt\nOhwskHR7uX+bpDPK7cdLulnS/eWAtWu6XSsJOCLqb0AP4STNA24ALgZWAJdJWtFU7ArggO0zgeuB\n68rtlwILbP8w8GbgF2aTczuZC+IYctyC3nuf9Dur2T2/+bt9Hdercz90ZV/HnfrZfT0fM3PwYF/X\nijEwuDbgVcCk7b0Akm4D1vDijgdrgI+UnzcDH5ekMoqFkuYDr6AYsPZMp4ulBhwR9dbbQIxFs+MV\nymVD09kWA482rE+V21qWsT0NPA2cRpGMvwc8ATwC/Jbt/Z1CTw04Imqvh0EW+2x36lDQqqG4+ezt\nyqwCDgOvA04F/kTSV2Zr062kBhwR9Tejakt3U8DShvUlwOPtypTNDacA+4F3U0xadsj2k8Cf0qX3\nWBJwRNSeXG2pYDuwXNKyctzDOmCiqcwEsL78vBa427Ypmh3ersJC4FzgW50uliaIiKi3AQ4ztj0t\n6SpgKzAPuMn2bkkbgR22J4AbgVskTVLUfNeVh98AfBr4JkUzxadt39fpeknAEVFzg53pzPYWYEvT\ntmsbPh+k6HLWfNyzrbZ3kgQcEfV3LA9FjogYazOjDqA/ScARUW+z/YBrKAk4ImqvYg+HsZMEHBH1\nV9MEnH7AEREjkhpwRNRemiAiIkbBVB1mPHaSgCOi/mpaA67cBixpnqRvSPpfwwwoIqJXA5wLYk71\n8hDu/cCeYQUSEdG3Y/W19ACSlgD/BPjUcMOJiOhDTRNw1Tbg3wF+BXjVEGOJiOjZuDYvVNG1Bizp\nx4Enbe/sUm7D7Gs+DvH8wAKMiOhqcBOyz6kqTRBvBS6R9DBwG8WEw7c2F7K9yfZK2yuPp/eXQ0ZE\n9OuYfQhn+xrbS2yfQTHx8N22f2bokUVEVHWMtwFHRIynMa3dVtFTArb9VeCrQ4kkIqJfL4cEHBEx\njpQJ2WPUZp7vvffJaV+b6uta537oyr6O61W/8U33cS8i5loScETUX02bIDIfcETUW8UuaFUf1Ela\nLelBSZOSrm6xf4Gk28v92ySd0bDvLElfl7Rb0v2STux0rSTgiKi/AXVDkzQPuAG4GFgBXCZpRVOx\nK4ADts8ErgeuK4+dD9wKXGn7DcAFwKFO10sCjoj6G1w/4FXApO29tl+gGHy2pqnMGuDm8vNm4EJJ\nAi4C7rP95wC2n7J9uNPFkoAjotZE0QuiygIsmp0yoVw2NJ1uMfBow/pUua1lGdvTwNPAacDfByxp\nq6Rdkn6lW+x5CBcR9dbbQIx9tld22N9qwojms7crMx/4R8BbgOeAuyTttH1Xu4ulBhwR9Te4Jogp\nYGnD+hLg8XZlynbfU4D95fY/tr3P9nPAFuBNnS6WBBwR9Te4BLwdWC5pmaQTKOa/mWgqMwGsLz+v\npZgfx8BW4CxJJ5WJ+R8DD3S6WJogIqL2BjUXhO1pSVdRJNN5wE22d0vaCOywPQHcCNwiaZKi5ruu\nPPaApI9RJHEDW2x/odP1koAjov4GOBDD9haK5oPGbdc2fD4IXNrm2FspuqJVkgQcEfXmzAURETE6\nNR2KnAQcEbX3spgPOMace/+vcHrqsb4udepn9/V1XK/6ntWsj3sRNVbTH3cScETU25i+bqiKJOCI\nqDVR3yaIKq+lXyrpjyTtKadYe/9cBBYRUVVd34pcpQY8Dfyy7V2SXgXslHSn7Y4jPCIi5swYJtcq\nqryW/gnbu8rP3wX28NLZgSIiRufl8Fr6cub3NwLbhhFMRETPxrR5oYrKCVjSK4HPAh+w/UyL/RuA\nDQAnctLAAoyI6OpYTsCSjqdIvp+x/blWZWxvAjYBnKzX1PR2REQdHbNDkctXbdwI7LH9seGHFBHR\nm7o2QVSZD/itwM8Cb5d0b7m8a8hxRURUU/UB3Bgm6a41YNv/m9av4IiIGA9jmFyryEi4iKi1Oo+E\nSwJ+uetz0pqZgwcHHEhE/zRTzwycBBwR9Tam7btVJAFHRO2lCSIiYlSSgCMiRqOuNeAq/YAjIsbb\nAPsBS1ot6UFJk5KubrF/gaTby/3byjlyGve/XtKzkj7Y7VpJwBFRb+Vbkass3UiaB9wAXAysAC6T\ntKKp2BXAAdtnAtcD1zXtvx74YpXQk4AjotZm+wEPaEL2VcCk7b22XwBuA9Y0lVkD3Fx+3gxcWE7Z\ngKSfAPYCu6tcLAk4IurPrrZ0txh4tGF9ipfOf36kjO1p4GngNEkLgV8Ffr1q2HkIFxG118NDuEWS\ndjSsbypncjxyqhbHNJ+9XZlfB663/WxZIe4qCTgi6q23gRj7bK/ssH8KWNqwvgR4vE2ZKUnzgVOA\n/cA5wFpJ/xF4NTAj6aDtj7e7WBJwRNTeAOcD3g4sl7QMeAxYB7y7qcwEsB74OrAWuNu2gbcdiUf6\nCPBsp+QLScARcQwYVAK2PS3pKmArMA+4yfZuSRuBHbYnKOZHv0XSJEXNd12/10sCjoh6M31PKtXy\ndPYWYEvTtmsbPh8ELu1yjo9UuVYScETUXl1HwiUBR0T9JQFHRMy9Ok/IXmkgRrex0RERI2OjmWrL\nuOmagCuOjY6IGJ2avpSzSg24ytjoiIiRGeBcEHOqShtwq7HR5zQXkrQB2ABwIicNJLiIiK4MjGHz\nQhVVasBVxkZje5PtlbZXHs+Co48sIqKqmjZBVKkBVxkbHRExMuPYvFBFlQRcZWx0RMTIjGMPhyq6\nJuB2Y6OHHllERBVj2rxQRaWBGK3GRkdEjINiIEY9M3BGwkVE/Q1uOso5NZQE/F0O7PuKN//fFrsW\nAfuGcc0ejUMc4xADjEcc4xADjEcc4xADzF0c3z+Ik6QG3MD297XaLmlHl9no58Q4xDEOMYxLHOMQ\nw7jEMQ4xjFMclRzrbcAREeNrPOd5qCIJOCLqL00QlWzqXmROjEMc4xADjEcc4xADjEcc4xADjE8c\n3Xmg74SbU3JNvzkiIgBOfuVin/Mjv1ip7Ff+z7/bOU5t22mCiIj6q2k9stKE7L3qNoG7pAWSbi/3\nb5N0xhBiWCrpjyTtkbRb0vtblLlA0tOS7i2Xa1ud6yjjeFjS/eX5d7TYL0n/ubwX90l60xBi+MGG\nf+O9kp6R9IGmMgO/F5JukvSkpG82bHuNpDslPVT+/6ltjl1flnlI0vohxPGbkr5V3vPPS3p1m2M7\n/vyOMoaPSHqs4Z6/q82xA3shQps4bm+I4WFJ97Y5diD3Yhg0M1NpGTcDT8AVJ3C/Ajhg+0zgeuC6\nQccBTAO/bPsfAOcC72szkfyf2D67XDYOIQ6AHy3P3+pPn4uB5eWyAfjEoC9u+8HZfyPwZuA54PMt\nig76Xvx3YHXTtquBu2wvB+4q119E0muAD1NMe7oK+HC7RH0UcdwJ/JDts4C/AK7pcHynn9/RxABw\nfcM9f8lo0yG8EOElcdj+qYb/Pj4LfK7D8YO4F4NlioEYVZYK+q1ASnqHpJ3ll9ROSW/vdq1h1ICr\nTOC+Bri5/LwZuFBSq2kv+2b7Cdu7ys/fBfZQzG08btYAv+fCPcCrJZ0+xOtdCHzbdquBMgNl+2vA\n/qbNjT/7m4GfaHHoO4E7be+3fYAiWbZKXn3HYfvLtqfL1XsoZvkbmjb3ooqBvhChUxzl7+BPAn/Q\n7/lHQRi52tL1XEdXgdwH/FPbPwysB27pdr1hJOBWE7g3J74jZcpfgqeB04YQCwDlN9QbgW0tdp8n\n6c8lfVHSG4ZweQNfLr8RN7TYX+V+DdI62v+CDfteALzW9hNQfEkCf6dFmbm+Jz8HfLHNvm4/v6N1\nVdkMclObWv5c3ou3Ad+x/VCb/cO+F/2zqy3d9V2BtP0N27NT9e4GTpTUcXL0YSTgKhO4V5rkfRAk\nvZLiz6oP2H6mafcu4Ptt/wjwX4D/MYQQ3mr7TRTfqO+TdH5ziC2OGda9OAG4BPjDFrvn4l5UNZf3\n5N9SNFd9pk2Rbj+/o/EJ4O8BZwNPAL/dKsQW24b1yOkyOtd+h3kvjs7gEvCgKpD/AviG7ec7XWwY\nCbjKBO5HykiaD5xCf3+edSTpeIrk+xnbL2nXsv2M7WfLz1uA4yUtGmQMs9+Itp+kaHdd1VRkLie8\nvxjYZfs7LeIc+r0ofWe2iaX8/ydblJmTe1I+3Ptx4Kfdpj9mhZ9f32x/x/Zh2zPAf2tz7rm6F/OB\nfw7c3q7MMO/FUemtDXiRpB0NS3NN/qgrkOVfj9cBv9At9GEk4CMTuJc1rnXARFOZCYo2EoC1wN3t\nfgH6VbZn3Qjssf2xNmX+7mzbs6RVFPfjqQHGsFDSq2Y/AxcB32wqNgG8R4Vzgadn/0QfgrY1nGHf\niwaNP/v1wP9sUWYrcJGkU8s/yy8qtw2MpNXArwKX2H6uTZkqP7+jiaGxrf+ftTl3ld+nQfgx4Fu2\np1rtHPa9OFo99ILYN/vqtHJpHnByVBVISUsovpzeY/vb3eIeeD/gdhO4S9oI7LA9QZEYb5E0WQa+\nbtBxAG8Ffha4v6Fbzb8BXl/G+bsUyf8XJU0DfwOsG/AXwWuBz5d5bT7w+7a/JOnKhhi2AO8CJil6\nJ7x3gNc/QtJJwDto+FZuimPg90LSHwAXUNQ6pih6NvwH4A5JVwCPAJeWZVcCV9r+edv7JX2UIvkA\nbLTd919IbeK4BlgA3Fn+fO6xfaWk1wGfsv0u2vz8BhjDBZLOpqg9PUz5s2mMYdAvRGgVh+0bafFs\nYFj3YvAqNy9UUeUNQLOViK/TUIFU0ZXxC8A1tv+0ysUyEi4iau2UV5zu8878uUplt37z33cdCaei\nP/bv8LdfeL/RWIGUdCJFD4c3UlYgbe+V9GsUX+yNDzEvKptsWl8rCTgi6uyUV5zu836gYgJ+oHsC\nnksZihwRtZcJ2SMiRiUJOCJiBGw4PH7zPFSRBBwR9ZcacETEiCQBR0SMgIG8Ey4iYhQMThtwRMTc\nM3kIFxExMmkDjogYkSTgiIhRGOhkPHMqCTgi6s3AGL5ws4ok4Iiov9SAIyJGIUORIyJGw+D0A46I\nGJGMhIuIGJG0AUdEjICdXhARESOTGnBExCgYHz486iD6kgQcEfVW4+kojxt1ABERR80z1ZYKJK2W\n9KCkSUlXt9i/QNLt5f5tks5o2HdNuf1BSe/sdq0k4IioNQOecaWlG0nzgBuAi4EVwGWSVjQVuwI4\nYPtM4HrguvLYFcA64A3AauC/ludrKwk4IurNHmQNeBUwaXuv7ReA24A1TWXWADeXnzcDF0pSuf02\n28/b/ktgsjxfW2kDjojaG+BDuMXAow3rU8A57crYnpb0NHBauf2epmMXd7pYEnBE1Np3ObD1K968\nqGLxEyXtaFjfZHtTw7paHNPcdtGuTJVjXyQJOCJqzfbqAZ5uCljasL4EeLxNmSlJ84FTgP0Vj32R\ntAFHRPyt7cByScsknUDxUG2iqcwEsL78vBa427bL7evKXhLLgOXAn3W6WGrAERGlsk33KmArMA+4\nyfZuSRuBHbYngBuBWyRNUtR815XH7pZ0B/AAMA28z3bHxmm5pkP4IiLqLk0QEREjkgQcETEiScAR\nESOSBBwRMSJJwBERI5IEHBExIknAEREjkgQcETEi/x/w2xgSECLxeQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cd37810>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dy=0.5\n",
    "dxmin = 0.5\n",
    "microspec = psf(sizey=10,sizex=20,dy=dy,dx=dxmin)\n",
    "dxlist = np.arange(dxmin)\n",
    "plt.imshow(microspec)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Equivalent noise pixel computation for matched filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.90257591762\n",
      "2.90257591762\n"
     ]
    }
   ],
   "source": [
    "dy=0.2\n",
    "dxmin = 0.2\n",
    "p = psf(sizex=100,sizey=100,dy=dy,dx=dxmin)\n",
    "p[p<np.amax(p)/2.]=0.0\n",
    "print(np.sum(p)**2/np.sum(p**2))\n",
    "V = np.sum(p)\n",
    "sharpness = np.sum(p**2/V**2)\n",
    "beta = 1./sharpness\n",
    "print beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10000\n",
    "vals=np.zeros(N)\n",
    "dxlist = np.zeros(N)\n",
    "dylist = np.zeros(N)\n",
    "\n",
    "for i in range(N):\n",
    "    dylist[i] = np.random.uniform(-0.5,0.5)\n",
    "    dxlist[i] = np.random.uniform(-0.5,0.5)\n",
    "    p = psf(sizex=20,sizey=20,dy=dylist[i],dx=dxlist[i])\n",
    "    #p[p<np.amax(p)/2.]=0.0\n",
    "    vals[i] = 1./np.sum(p**2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x115e5add0>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAFIlJREFUeJzt3X2wXHV9x/H3N+FBVBRiHpoq6bUa\nYaqlON0qbceWSjNkrBV8AKHKxDaY2sepdixxcAbbypg2U7BP005kkHSqNGrbkEIbht6WAR3AuZgU\nsK0EKtLETBKLFNoqCnz7x/2hy3Jvdu/e3XP27L5fM5ndc+45ez93k9zPnvM7D5GZSJK0pO4AkqTR\nYCFIkgALQZJUWAiSJMBCkCQVFoIkCbAQJEmFhSBJAiwESVJxTN0BFmL58uU5NTVVdwxJapS77rrr\na5m5ottyjSqEqakpZmZm6o4hSY0SEV/pZTl3GUmSAAtBklRYCJIkwEKQJBUWgiQJsBAkSYWFIEkC\nLARJUmEhSJIAC0GSVFgIkiTAQpAkFRaCJAnooRAi4pqIOBwR97bNWxYRN0fEvvJ48jzr7o6IRyLi\nho75r4+IL0TEvRGxPSIaddVVSRpHvWwhXAus75i3GZjOzLXAdJmey1bg4vYZEbEE2A5cmJmvAr4C\nbFhAZknSEHQthMy8FXi4Y/a5zP5SpzyeN8+608BjHbNfBDyemfeV6ZuBt/YaWJI0HP2OIazKzIMA\n5XHlAtb9GnBsRLTK9NuAU/rMIUkakMoHlTMzgQuBqyLi88xuQTwx3/IRsSkiZiJi5siRI1XFlKSJ\n028hHIqI1QDl8fBCVs7M2zPzdZn5GuBWYN9Rlt2Wma3MbK1Y0fWWoJKkPvVbCLv47kDwBuD6hawc\nESvL4/HApcCf95lDkjQgvRx2eh1wO3BqROyPiI3AFmBdROwD1pVpIqIVEVe3rXsb8Gng7LLuOeVL\n74+IfwPuBv4uM/9poD+VJGnBYnaXfjO0Wq2cmZmpO4YkNUpE3JWZrW7LeaayJAmwECRJhYUgSQIs\nBElSYSFIkgALQZJUWAiSJMBCkCQVFoIkCbAQJEmFhSBJAiwESVJhIUiSAAtBklRYCJIkwEKQJBUW\ngiQJsBAkSYWFIEkCLARJUmEhSJIAC0GSVFgIkiTAQpAkFRaCJAmwECRJhYUgSQIsBElSYSFIkoAe\nCiEiromIwxFxb9u8ZRFxc0TsK48nz7Pu7oh4JCJu6Jh/dkR8ISL2RsRnI+Lli/9RJEmL0csWwrXA\n+o55m4HpzFwLTJfpuWwFLp5j/p8B78jMM4BPAh/sKa0kaWi6FkJm3go83DH7XGB7eb4dOG+edaeB\nx+b6EvCC8vyFwFd7CStJGp5j+lxvVWYeBMjMgxGxcoHrXwL8fUR8A3gUOLPPHJKkAalrUPm9wBsy\n8yXAx4Er51swIjZFxExEzBw5cqSygJI0afothEMRsRqgPB7udcWIWAH8UGbeWWbtAH5svuUzc1tm\ntjKztWLFij7jSpK66XeX0S5gA7ClPF6/gHW/DrwwIl6RmfcB64B/6zOHJDXKzj0H+NCuL/LIN749\n59ePWRI88VQCEMA7zlzDh8/7wUqydS2EiLgOOAtYHhH7gcuZLYJPRcRG4CHg/LJsC3hPZl5Spm8D\nTgOeX9bdmJk3RcS7gb+OiKeYLYhfGPhPJklD9sGd9/DJOx+i/P4eiCfaXiyBv7zjIYBKSiEyB/iT\nDFmr1cqZmZm6Y0iaIK+94mYOPfatWjMsjeCBj7yh7/Uj4q7MbHVbrt9dRpI0NnbuOcBv7Nhbd4x5\nPVnRB3cLQdJE6bYPfxQtjajk+1gIksba6Zfv5tHHn6w7xqJc9NpTKvk+FoKksfHBnfd8ZxC2qUb6\nKCNJGmWjMOi7GM87bilXvPkHOe/VL647ioUgqVnWXXkL+w7/b90xerI04A8uOGMkftn3wkKQNNKa\nVAAvPukE3n/OqY0pgE4WgqSR05SxgJNOOJYPvemVjS2AThaCpJHQlKOBmr4VcDQWgqTavONjt/O5\nBzpvtzJa1q58Hje/76y6Y1TCQpBUqSaUAIz3lsB8LARJQ3faZX/PN58c/eumjduYwEJZCJKGoikD\nw89ZGvz7Ff1fOG6cWAiSBqpJh4lO0vhALywESYvWpBKA2UtCXPX25pwwVhULQVJfmjI43O6dFV4X\nqIksBEk927nnAO/dsZfRHx5+plG6XtAosxAkddWUk8Y6uWtoYSwESXMa9buIdePuoYWzECQ9Q9MG\niDtVfQ+BcWIhSGLnngO8b8denqo7yCIsjeCi155iESyChSBNsKacPNbNj79sGZ9494/WHaPxLARp\nAjX9LmPtVp14nGUwIBaCNEGaPj7Qzl1Eg2chSBOgqYeNzuejHko6FBaCNKbGYaC40xLgSstgaCwE\nacyMy0BxJ7cKhs9CkMbEOA0Ud3rnmWssgwpYCFLDNeXmM/3yjOPqdC2EiLgGeCNwODNfVeYtA3YA\nU8CDwAWZ+fU51t0NnAl8NjPf2Db/NuDEMrkS+Hxmnreon0SaMONeBKtOPI47L1tXd4yJsqSHZa4F\n1nfM2wxMZ+ZaYLpMz2UrcHHnzMx8XWaekZlnALcDf9NzYmmC7dxzgKnNNzK1+caxLoOPvv0My6AG\nXbcQMvPWiJjqmH0ucFZ5vh24Bbh0jnWnI+KszvlPi4gTgdcDP99LWGlSjdP5A938+MuWOV5Qk37H\nEFZl5kGAzDwYESv7fJ03M7ul8Wif60tjbdzOH+jG8YJ61T2ofBFw9dEWiIhNwCaANWvWVJFJGglT\nm2+sO0JlHC8YDf0WwqGIWF22DlYDhxf6AhHxIuA1zG4lzCsztwHbAFqt1vjuNJUY33MIjsYb3Y+O\nfgthF7AB2FIer+/jNc4HbsjMb/aZQRobTb8ZTb+WBpbBCOnlsNPrmB1AXh4R+4HLmS2CT0XERuAh\nZn+5ExEt4D2ZeUmZvg04DXh+WXdjZt5UXvrC8jrSxBr3Q0e7+YMLzqg7gtpEZnP+MbZarZyZmak7\nhrRo7/jY7XzugYfrjlEb73VcrYi4KzNb3Zare1BZmiiTOEbQyTGD0WUhSBWY1DGCTt7ZbLRZCNKQ\n7NxzgPfu2EtzdsoOl2Uw+iwEaQheuvlGi6DNCccusQwawEKQBmjSjxqaz0fecnrdEdQDC0EagEm6\n1tBCeS+D5rAQpEUY55vSLNbJzz2Wy3/2lZZBg1gIUh/cNXR0XqSumSwEqUceOtqbY5dgGTSUhSB1\nYREszNbzvRxFU1kI0jwsgoXz5jbNZiFIHSyC/nmuQbNZCFLhdYb6F8CXt/xM3TG0SBaCxGTdnayb\nE45dwkfecrq7fiaQhaCJNUm7hrzCqHphIWjijOO9CLy/gAbBQtBEafquoWOXzB7W6S9+DYOFoLHX\n5C0Cd/WoShaCxtbpl+/m0cefrDtGz7zcg+pmIWjsNOXKo94wRqPGQtDYGPUrj7oFoFFnIajRRvVk\nMo/6URNZCGqkUTuHYAlwpQWghrMQ1CijdMSQYwAaNxaCGmFUimDVicdx52Xr6o4hDYWFoJE2Cncm\n+6i7gjQhLASNpDrPKPZkME0qC0Ejo+4jhtwS0KSzEFS7unYLuSUgPVPXQoiIa4A3Aocz81Vl3jJg\nBzAFPAhckJlfn2Pd3cCZwGcz841t8wP4MHA+8CTwZ5n5R4v9YdQcdR02evwxS/i9t3qtf2kuvWwh\nXAv8CfAXbfM2A9OZuSUiNpfpS+dYdyvwXOAXO+a/CzgFOC0zn4qIlQvMrYaqY2vgBccv5e7fXl/p\n95SaqGshZOatETHVMftc4KzyfDtwC3MUQmZOR8RZnfOBXwJ+LjOfKssd7jWwmumlm2+kyhp4ztLg\n3694Q4XfUWq+fscQVmXmQYDMPNjHJ/yXAW+PiDcDR4Bfz8x9cy0YEZuATQBr1qzpM67qUMduIQeG\npf7VNah8PPDNzGxFxFuAa4DXzbVgZm4DtgG0Wq16D0hXT6reLeTgsDQY/RbCoYhYXbYOVgML3eWz\nH/jr8vxvgY/3mUMjospLTnvhOGk4+i2EXcAGYEt5vH6B6+8EXs/slsFPAvf1mUM1e/kHbuSJijYG\nvHy0NFy9HHZ6HbMDyMsjYj9wObNF8KmI2Ag8xOzho0REC3hPZl5Spm8DTgOeX9bdmJk3lfU/ERHv\nBf4HuGTQP5iGxxKQxlNkNme3fKvVypmZmbpjTKQqdwl5FVFpsCLirsxsdVvOM5V1VFVdU8gSkOpn\nIehZqrw5/YNbfqaS7yOpOwtBQLUnjlkC0miyECZYlZeY9oQxafRZCBOk6stLWwJSs1gIY67K8QDw\nMFGpySyEMVT13cYcE5DGg4UwBqougAC+bAlIY8dCaKAqzxRu55aANN4shAao64bzbglIk8VCGDFV\n30imk4PC0uSyEGr02itu5tBj36o7hruCJAEWQiXquI/w0XjdIElzsRAGpO5dPUez6sTjuPOydXXH\nkDTiLIR51HE/4EFxMFhSPyauEHbuOcDWm77EgUe+UXeUgXIcQNJiTVQh7NxzgA/8zT1849vVXcph\nGNwFJGkYJqoQtt70pUaWgZ/+JVVhogrhqyO+m+gFxy/l7t9eX3cMSRNqogrhe086YSTGDtzlI2kU\nTVQhvP+cUysbQ/BeAJKaZqIK4elf0P0cZeTJXJLG3UQVAsyWgp/cJenZltQdQJI0GiwESRJgIUiS\nCgtBkgRYCJKkomshRMQ1EXE4Iu5tm7csIm6OiH3l8eR51t0dEY9ExA0d86+NiC9HxN7y54zF/yiS\npMXoZQvhWqDzegqbgenMXAtMl+m5bAUunudr78/MM8qfZl5nWpLGSNdCyMxbgYc7Zp8LbC/PtwPn\nzbPuNPDYYgJKkqrR7xjCqsw8CFAeV/bxGldExN0RcVVEHN9nDknSgNQ1qPwB4DTgR4BlwKXzLRgR\nmyJiJiJmjhw5UlU+SZo4/RbCoYhYDVAeDy9k5cw8mLMeBz4OvOYoy27LzFZmtlasWNFnXElSN/0W\nwi5gQ3m+Abh+ISu3lUkwO/5w79HXkCQNWy+HnV4H3A6cGhH7I2IjsAVYFxH7gHVlmohoRcTVbeve\nBnwaOLuse0750ici4h7gHmA58OFB/lCSpIWLzKw7Q89arVbOzMzUHUOSGiUi7srMVrflPFNZkgRY\nCJKkwkKQJAEWgiSpsBAkSYCFIEkqLARJEmAhSJIKC0GSBFgIkqTCQpAkARaCJKmwECRJgIUgSSos\nBEkSYCFIkgoLQZIEWAiSpMJCkCQBFoIkqbAQJEmAhSBJKiwESRJgIUiSCgtBkgRYCJKkwkKQJAEW\ngiSpsBAkSUAPhRAR10TE4Yi4t23esoi4OSL2lceT51l3d0Q8EhE3zPP1P46I/+k/viRpUHrZQrgW\nWN8xbzMwnZlrgekyPZetwMVzfSEiWsBJvcWUJA1b10LIzFuBhztmnwtsL8+3A+fNs+408Fjn/IhY\nymxZ/NZCwkqShqffMYRVmXkQoDyuXOD6vwrsevo1JEn1O6bqbxgR3wucD5zV4/KbgE0Aa9asGV4w\nSZpw/W4hHIqI1QDl8fAC1n018HLg/oh4EHhuRNw/38KZuS0zW5nZWrFiRZ9xJUnd9FsIu4AN5fkG\n4PpeV8zMGzPzezJzKjOngP/LzJf3mUOSNCC9HHZ6HXA7cGpE7I+IjcAWYF1E7APWlWkiohURV7et\nexvwaeDssu45w/ghJEmL13UMITMvmudLZ8+x7AxwSdv063p4/ed3W0aSNHyeqSxJAiwESVJhIUiS\nAAtBklRYCJIkACIz687Qs4g4AnxlAC+1HPjaAF5n2Mw5WE3JCc3Jas7BGlbO78vMrmf2NqoQBiUi\nZjKzVXeObsw5WE3JCc3Jas7Bqjunu4wkSYCFIEkqJrUQttUdoEfmHKym5ITmZDXnYNWacyLHECRJ\nzzapWwiSpA5jWwgRsT4ivhQR90fEs+75HBHHR8SO8vU7I2Kq+pTfydIt609ExBci4omIeFsdGUuO\nbjnfFxH/GhF3R8R0RHzfiOZ8T0TcExF7I+KzEfEDo5izbbm3RUSW+5DXoof39F0RcaS8p3sj4pK5\nXqfunGWZC8q/0y9GxCerzlgydHs/r2p7L++LiEcqCZaZY/cHWAo8AHw/cBzwL8APdCzzy8Cfl+cX\nAjtGOOsUcDrwF8DbRjjnTwHPLc9/qY73tMecL2h7/iZg9yjmLMudCNwK3AG0Rvjv/l3An9SRb4E5\n1wJ7gJPL9MpRzNmx/K8B11SRbVy3EF4D3J+Z/5GZ3wL+Cji3Y5lzge3l+WeYvWdDVJjxaV2zZuaD\nmXk38FQN+Z7WS85/zsz/K5N3AC+pOCP0lvPRtsnnAXUMpPXybxTgd4HfB75ZZbgOvWatWy853w38\naWZ+HSAzF3K3x0FZ6Pt5EXBdFcHGtRBeDPxn2/T+Mm/OZTLzCeC/gRdVkm6eHMVcWUfBQnNuBP5h\nqInm1lPOiPiViHiA2V+2v15RtnZdc0bEq4FTMvOGKoPNode/+7eW3YWfiYhTqon2DL3kfAXwioj4\nXETcERHrK0v3XT3/Xyq7XV8K/FMFuca2EOb6pN/5KbCXZaowKjm66TlnRLwTaAFbh5pobj3lzMw/\nzcyXAZcCHxx6qmc7as6IWAJcBfxmZYnm18t7+nfAVGaeDvwj3936rlIvOY9hdrfRWcx+8r46Ik4a\ncq5OC/k/fyHwmcx8coh5vmNcC2E/0P4J5SXAV+dbJiKOAV4IPFxJunlyFHNlHQU95YyInwYuA96U\nmY9XlK3dQt/PvwLOG2qiuXXLeSLwKuCWiHgQOBPYVdPActf3NDP/q+3v+2PAD1eUrV2v/++vz8xv\nZ+aXgS8xWxBVWsi/0QupaHcRMLaDyscA/8HsptbTgzav7FjmV3jmoPKnRjVr27LXUt+gci/v6auZ\nHSxbO+J/92vbnv8sMDOKOTuWv4X6BpV7eU9Xtz1/M3DHiOZcD2wvz5czu+vmRaOWsyx3KvAg5Xyx\nSrLV8Q+sojf9DcB95RfUZWXe7zD7yRXgOcCngfuBzwPfP8JZf4TZTxX/C/wX8MURzfmPwCFgb/mz\na0Rz/iHwxZLxn4/2i7jOnB3L1lYIPb6nHynv6b+U9/S0Ec0ZwJXAvwL3ABeOYs4y/SFgS5W5PFNZ\nkgSM7xiCJGmBLARJEmAhSJIKC0GSBFgIkqTCQpAkARaCJKmwECRJAPw/8NOdfVSLF/kAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x115b29c90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.sqrt(dylist**2+dxlist**2),vals)\n",
    "# do colorgrid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aperture photometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### signal within 2x2 pixel aperture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}