jmeyers314/Higher Order PSF Moments Notebook.ipynb

## Higher Order PSF Moments Notebook.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.625979Z",
     "start_time": "2019-01-22T22:52:26.955223Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import galsim\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.630297Z",
     "start_time": "2019-01-22T22:52:27.627940Z"
    }
   },
   "outputs": [],
   "source": [
    "psfKwargs = {'nx':32, 'ny':32, 'scale':0.2}\n",
    "galKwargs = {'nx':128, 'ny':128, 'scale':0.2}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.633713Z",
     "start_time": "2019-01-22T22:52:27.631621Z"
    }
   },
   "outputs": [],
   "source": [
    "def getImg(obj, **kwargs):\n",
    "    return obj.drawImage(**kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.639158Z",
     "start_time": "2019-01-22T22:52:27.636041Z"
    }
   },
   "outputs": [],
   "source": [
    "def getHSM(obj, **kwargs):\n",
    "    img = getImg(obj, **kwargs)\n",
    "    mom = galsim.hsm.FindAdaptiveMom(img)\n",
    "    sigma = mom.moments_sigma\n",
    "    shape = mom.observed_shape\n",
    "    flux = mom.moments_amp\n",
    "    center = mom.moments_centroid\n",
    "    scale = img.scale\n",
    "    return flux, center*scale, sigma*scale, shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.644845Z",
     "start_time": "2019-01-22T22:52:27.641774Z"
    }
   },
   "outputs": [],
   "source": [
    "def getMatchingTransform(star, target):\n",
    "    targetFlux, targetCenter, targetSigma, targetShape = getHSM(target, **psfKwargs)\n",
    "    starFlux, starCenter, starSigma, starShape = getHSM(star, **psfKwargs)\n",
    "\n",
    "    dFlux = targetFlux/starFlux\n",
    "    dCenter = targetCenter - starCenter\n",
    "    dScale = targetSigma / starSigma\n",
    "    dShear = targetShape - starShape\n",
    "    return dFlux, dCenter, dScale, dShear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.649025Z",
     "start_time": "2019-01-22T22:52:27.646168Z"
    }
   },
   "outputs": [],
   "source": [
    "def matchObjToTarget(obj, targetObj, niter=3):\n",
    "    for _ in range(niter):\n",
    "        dFlux, dCenter, dScale, dShear = getMatchingTransform(obj, targetObj)\n",
    "        obj = obj.shear(dShear).dilate(dScale).shift(dCenter)*dFlux\n",
    "    return obj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.682595Z",
     "start_time": "2019-01-22T22:52:27.650385Z"
    }
   },
   "outputs": [],
   "source": [
    "psf1 = galsim.VonKarman(lam=750, r0=0.15, L0=25)\n",
    "psf2 = galsim.Kolmogorov(lam=750, r0=0.15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.775426Z",
     "start_time": "2019-01-22T22:52:27.683936Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.8743342757225037, galsim.PositionD(x=3.3000000241894707, y=3.2999999946789575), 0.38043971061706544, galsim.Shear((-2.281066224796291e-09-7.065105445125397e-11j)))\n",
      "\n",
      "(0.9052954912185669, galsim.PositionD(x=3.300000027531583, y=3.299999997889095), 0.4547154903411865, galsim.Shear((-1.432827190228636e-09+3.0518898519460436e-09j)))\n",
      "\n",
      "(0.9052950739860535, galsim.PositionD(x=3.300000029210454, y=3.2999999991961455), 0.4547145366668701, galsim.Shear((2.1831827012874793e-08+3.5881382309810322e-09j)))\n"
     ]
    }
   ],
   "source": [
    "matchedPSF = matchObjToTarget(psf1, psf2, niter=3)\n",
    "print(getHSM(psf1, **psfKwargs))\n",
    "print()\n",
    "print(getHSM(psf2, **psfKwargs))\n",
    "print()\n",
    "print(getHSM(matchedPSF, **psfKwargs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.779570Z",
     "start_time": "2019-01-22T22:52:27.776804Z"
    }
   },
   "outputs": [],
   "source": [
    "def getGalShape(conv, psf):\n",
    "    convImg = getImg(conv, **galKwargs)\n",
    "    psfImg = getImg(psf, **galKwargs)\n",
    "    plt.imshow(convImg.array)\n",
    "    plt.show()\n",
    "    return galsim.hsm.EstimateShear(convImg, psfImg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.960203Z",
     "start_time": "2019-01-22T22:52:27.781131Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD8CAYAAAB+fLH0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFkpJREFUeJzt3X2MXFd5x/Hvb2Z2vbHB2A4kGDslQXIpIYKCLBqgrRCmJaGIpBJUoahYEMmqREt4kSCBP1D/A4F4k4DWAkpaRSE0QBNFvDSEINpKBBxAIcGEuAESE4MDcV7Ixrs7O0//OGdm7pmd9a53XvaF30dazdxz79x7fD3z3Oece+69igjMzNpqq10BM1tbHBTMrOCgYGYFBwUzKzgomFnBQcHMCg4KZlZwUDCzgoOCmRUaq10BgEltiim2rHY1zDa0xzjxm4h42lLLrYmgMMUW/kT7VrsaZhvaN+L6XyxnOTcfzKzgoGBmBQcFMys4KJhZwUHBzAoOCmZWcFAws8KSQUHSZyUdl3RnpeyDkn4i6Q5JX5a0rTLvKklHJN0t6ZWjqriZjcZyMoXPARf1lN0MXBARzwN+ClwFIOl84DLgufkzn5RUH1ptzWzklgwKEfFt4KGesv+KiGae/A6wO7+/BPh8RMxExM+AI8CLhlhfMxuxYfQpvBn4an6/C7i/Mu9oLjOzdWKgax8kvRdoAte0i/os1vce8pIOAAcAptg8SDXMbIhWHBQk7QdeDeyL7sMjjgLnVBbbDTzQ7/MRcRA4CLBVO/zwCbM1YkXNB0kXAe8GXhMR05VZNwKXSdok6TxgD/DdwatpZuOyZKYg6VrgZcBTJR0F3kc627AJuFkSwHci4u8j4i5JXwB+TGpWvCUi5kdVeTMbPq2Fx8Zt1Y7w/RTMRusbcf3tEbF3qeU8otHMCg4KZlZwUDCzgoOCmRUcFMys4KBgZgUHBTMrOCiYWcFBwcwKDgpmVnBQMLOCg4KZFRwUzKzgoGBmBQcFMys4KJhZwUHBzAoOCmZWcFAws4KDgpkVHBTMrOCgYGYFBwUzKzgomFnBQcHMCksGBUmflXRc0p2Vsh2SbpZ0T37dnssl6eOSjki6Q9ILR1l5Mxu+5WQKnwMu6im7ErglIvYAt+RpgItJD5XdQ3rM/KeGU00zG5clg0JEfBt4qKf4EuDq/P5q4NJK+b9F8h1gm6Sdw6qsmY3eSvsUzo6IYwD59axcvgu4v7Lc0VxmZuvEko+iP03qU9b3sdaSDpCaGEyxecjVMLOVWmmm8Ot2syC/Hs/lR4FzKsvtBh7ot4KIOBgReyNi7wSbVlgNMxu2lQaFG4H9+f1+4IZK+RvzWYgLgUfazQwzWx+WbD5IuhZ4GfBUSUeB9wHvB74g6XLgPuB1efGvAK8CjgDTwJtGUGczG6Elg0JEvH6RWfv6LBvAWwatlJmtHo9oNLOCg4KZFRwUzKzgoGBmBQcFMys4KJhZwUHBzAoOCmZWcFAws4KDgpkVHBTMrOCgYGYFBwUzKzgomFnBQcHMCg4KZlZwUDCzgoOCmRUcFMys4KBgZgUHBTMrOCiYWcFBwcwKDgpmVnBQMLPCQEFB0tsl3SXpTknXSpqSdJ6k2yTdI+k6SZPDqqytIdLgf7YmrTgoSNoFvBXYGxEXAHXgMuADwEciYg9wArh8GBW1VTDqH7QDxZo0aPOhAZwhqQFsBo4BLweuz/OvBi4dcBtmNkYrDgoR8UvgQ6SnTh8DHgFuBx6OiGZe7Ciwa9BK2hislaP2WqnH77FBmg/bgUuA84BnAFuAi/ssGot8/oCkQ5IOzTGz0mqY2ZAt+Sj6U3gF8LOIeBBA0peAlwDbJDVytrAbeKDfhyPiIHAQYKt29A0cNmKnexTWiE5WRWuJ7eZ6hr8m4zDI//J9wIWSNksSsA/4MXAr8Nq8zH7ghsGqaEO3VFquWv8/QDUN/HfK7Q1SbxuKQfoUbiN1KH4f+FFe10Hg3cA7JB0BzgQ+M4R6mtmYDNJ8ICLeB7yvp/he4EWDrNdGYKkjbJ+j9KJH9UFEa8F6o1VpFvSuv1/TQnJTYoQ8otHMCgNlCrYOLNV3sKBIi8+vzNMy2/ax4Iheh1ZZplo3G4ieeai2eLaQNrCsetjyOVMws4IzhY1qmRnCgn4D1YqMIK0qT9dq1cLlVaP3SN5qEb2HolY3G2jXp28/gzOGsXBQ2GgW+7EuFgja5blMUvfH3xMM0rw+QaRX9cebf9ydZoTUDRSttFzUSIGhWG2e19ucsJFz88HMCs4UNpJ+WcIymwqdJkK93s0M6rXucgD1Wv/MIlvQqdgKaM2n5dpH/Pn5ImsA0Px8t0mxoBNS3WxhqWaEmxBD4UzBzArOFDaCJTKENKmFZfVat/+gXu+W5ffUclmj3l2m3tPfUNl2p6+g/TrfgvmUKURzPq9TqF1WrUu7rJMxdLOCBZ2Pi52mtKFwUFjPlhsM2uW9KX+tGwDUyF+FRqMbBNplE+k1JhqQ50V94bY1X3YgMtdEORionn/0c3PdZgPpCvuATiBRtfMRijMTNh5uPphZwZnCerXEOIEFnYo1lZ2JgOr1ThbQzgo0MQGTEwDEVLq9Zkymea2pBq1GOo5Eo92M6HYiqplea7M5O5htwmzKBjQz116822zoZAdB9I6FyM0JGz9nCmZWcKawkfQbSNS+D0JlUJLaHYkTjW6GMJmzgqlJYvMmAOY3T+bXtExzc53mVO70q/QptPsSGifTa/2JdJRvPFGnNp222a6ZIjodkp2Ow3pUOifnu/UFov+Nu2yEHBTWm+WORegZT1DtVOy8NhrdYHBGCgStLVPMPymVzT4lNSNmnpLWNbtVNDen9bXyN0cB9ZN5dY+nH/Cmx9LyrcdqTLSrUzkj0e5M7JxxaHWbNnGKZlHfIdA2dG4+mFnBmcIGUoxF6D39KHVHKFY6FTudie0mw9ZJZranDGH6zJRRnHxaWsfMjqC5Nd+oe6J92rFG/bG03KYTOYuYaNeoTm2+7Hysz9Zhrs9Yh3ZG09u/WNPCMhspZwpmVnCmsF4sY6BSX9XOxXrPCMWJBrEpdyZuyf0IT57giR1p/vTT0zaf2J2yg607H+OZ208AsLkxC8BvT27h/t9uS8tNbknrz6MRGyehkTsaGxN50FOj1v9SbFsz/L9iZgVnChtAMVCpT19Cp7xWDl+OiQYxlU83npHmzW6tcXJH+szJp6fG/NP+IGUHr3jG3bz0yT8FYItSpnDXzC5ubpwPwB0n08PAmo/k/olN0MqnH6LWrU/7/bJu01I50+CzDuPhoLAeLLfpUC1rp+aVC546HY3t14kGrZzWz0+l17kzRPNJeXXb0w//uWf+CoC/3HonfzaVmhL1vK3NtXs5vPkZABzedDYArUYKClGtduViqc5l1K2FN2Ox1efmg5kVnClsML3NhmK6p4Mv6uqMTGyn+a0JmJ/MIxQnUvNhSz1lDJtrM9SVzjfO5ysXp1ubeHw+dVI251K2MTmb1lWbC2pz+fqG+cprO0OoZA/Rc9n1wrtA27g4UzCzwkCZgqRtwKeBC0gXv70ZuBu4DjgX+DnwNxFxYqBa2soVnZA5Q6j1yR4yBdTy+KTZmfT1ODqdTjl+74lnMd36JQCPxxQA//vYH3LXb5+e1nsiZQyN36XPTzwR1E+mbKM2l0cgzTXTzVeAyK/Mz3ezh1bPvRP63Uuhb5kzi2EZtPnwMeBrEfFaSZPAZuA9wC0R8X5JVwJXkp4vaUPW7/LoZX2ucofl9vt2ml+fgcZ0Ws9s/pH/ZOosAB5vTvLfU3sAmG6meb84sZ1Hjz0ZgDN+lZoPUw+ldU0+2qKRL47SyXwJ9VwTmjnq5Ps30oqFzQafdVg1K24+SNoK/Dn5AbIRMRsRDwOXAFfnxa4GLh20kmY2PoNkCs8CHgT+VdLzgduBK4CzI+IYQEQck3TW4NW0oejp4FOzheba1ySkeRPTNSYfzovl5sbMbDpHeeShMzhSueYBoPFonS0Ppcxi6sG03jMeSuucfGyO+uOpk1Izs/lzze79GttXSVabD6fiZsNYDNLR2ABeCHwqIl4APE5qKiyLpAOSDkk6NMfMANUws2EaJFM4ChyNiNvy9PWkoPBrSTtzlrATON7vwxFxEDgIsFU7HO6HoRXQe0PVVuVUX8+TmWjOd69enE7t/Im6gNQ3UEt3UGPid2md81ONzoCkdmdkYzqYfDStb9Mj6XXykfTB+u9m0XQK+DqZMoWYm+v0KXQ6GqunJDv19ROiVsuKM4WI+BVwv6Rn56J9wI+BG4H9uWw/cMNANTSzsRr07MM/AtfkMw/3Am8iBZovSLocuA943YDbsEW0j6ILzkLAgt58taJzZO7cdr3WhJl0XKi3BzsF1ObSco0nUsaw1C3YGtNltlGfzv0I0zPdDGE29yk0m92bsrZfW61TZwjuSxirgYJCRPwQ2Ntn1r5B1msr1wkC7YJWZSxA+3Zmze5ox/aIx849FFstarP57s05KEw2Khc1tZ/H0kzrrTVb1NqnG9t3bs6BgNm51FyA7mnIuWbqWKzULSIWf7aDn/kwdh7RaGYFX/uwUbUzhPYoxojOjVLp3Cm5m1Go+lSnnCnUKtdI9Oo8Dao5nwYk5c+mstyR2JyvvM/zKqcfVzxQyU2HkXKmYGYFZwobSbQ6D2ZtP4uxkwFInScotJ/hSES3s3I+36qtWe8e8XuvtMyfKV77PUS21e1ILK5vAGi1umWVei/IENy5uGocFNaraHVuqlKchWj/mHqDQ+UxbJ3g0Ir0IBYg2j/kWn3BNRTVoNB3PEFnvEHPhU6tSsCofq7nB+8zDWuLmw9mVnCmsB5E9L8lW/uoWskYFoxZaF8RWatkC+0xDFL3OojODVgqD1lor7ffNivrj96RkpXpxUYqtuu76Ho75c4Sxs2ZgpkVnClsMN3+hd4jeq3Tv1A8xLWdNbSXW+IR990NVY7gvacYO+Ur6D/ot34bKweF9aL9I1lGMyIV9QyBrpyZ6CxTo3tWgPbqlxcUigBwijMHfccdOBisaW4+mFnBmcJ6s1TG0PM8iOqRekGTYp6Fy3OaR+o+R/3Tzg7AGcIa4kzBzArOFNarpU5TtvXpZ+jO0tJH8KWqcaprFZwdrEvOFMys4ExhPTtV/0JnmcrR+hT9DYPV4zSzDWcIa5qDwkawnOAAq3vDEgeCdcPNBzMrOFPYSPodjZc7QnHU9bB1w5mCmRWcKWx0o8oenA1sWA4Kv4/8g7ZTcPPBzAoOCmZWcFAws8LAQUFSXdIPJN2Up8+TdJukeyRdlx8pZ2brxDAyhSuAw5XpDwAfiYg9wAng8iFsw8zGZKCgIGk38FfAp/O0gJeTHksPcDVw6SDbMLPxGjRT+CjwLqA9qP5M4OGIyE8T4Siwa8BtmNkYrTgoSHo1cDwibq8W91m070lxSQckHZJ0aI6ZlVbDzIZskMFLLwVeI+lVwBSwlZQ5bJPUyNnCbuCBfh+OiIPAQYCt2uHRNGZrxIozhYi4KiJ2R8S5wGXANyPiDcCtwGvzYvuBGwaupZmNzSjGKbwbeIekI6Q+hs+MYBtmNiJDufYhIr4FfCu/vxd40TDWa2bj5xGNZlZwUDCzgoOCmRUcFMys4KBgZgUHBTMrOCiYWcFBwcwKDgpmVnBQMLOCg4KZFRwUzKzgoGBmBQcFMys4KJhZwUHBzAoOCmZWcFAws4KDgpkVHBTMrOCgYGYFBwUzKzgomFnBQcHMCg4KZlYY5KnT50i6VdJhSXdJuiKX75B0s6R78uv24VXXzEZtkEyhCbwzIp4DXAi8RdL5wJXALRGxB7glT5vZOjHIU6ePRcT38/vHgMPALuAS4Oq82NXApYNW0szGZyh9CpLOBV4A3AacHRHHIAUO4KxhbMPMxmPgoCDpScAXgbdFxKOn8bkDkg5JOjTHzKDVMLMhGSgoSJogBYRrIuJLufjXknbm+TuB4/0+GxEHI2JvROydYNMg1TCzIRrk7IOAzwCHI+LDlVk3Avvz+/3ADSuvnpmNW2OAz74U+DvgR5J+mMveA7wf+IKky4H7gNcNVkUzG6cVB4WI+B9Ai8zet9L1mtnq8ohGMys4KJhZwUHBzAoOCmZWcFAws4KDgpkVHBTMrOCgYGYFBwUzKzgomFnBQcHMCg4KZlZwUDCzgoOCmRUcFMys4KBgZgUHBTMrOCiYWcFBwcwKDgpmVnBQMLOCg4KZFRwUzKzgoGBmBQcFMyuMLChIukjS3ZKOSLpyVNsxs+EaSVCQVAc+AVwMnA+8XtL5o9iWmQ3XqDKFFwFHIuLeiJgFPg9cMqJtmdkQjSoo7ALur0wfzWVmtsYN8ij6U+n3NOooFpAOAAcAptg8omqY2ekaVVA4CpxTmd4NPFBdICIOAgcBJD34jbj+ceA3I6rP6XgqrkeV61Faz/V45nIWUkQsvdRpktQAfgrsA34JfA/424i46xSfORQRe4demdPkergev+/1GEmmEBFNSf8AfB2oA589VUAws7VjVM0HIuIrwFdGtX4zG421NKLx4GpXIHM9Sq5HacPXYyR9Cma2fq2lTMHM1oA1ERRW4zoJSedIulXSYUl3Sboil++QdLOke/Lr9jHVpy7pB5JuytPnSbot1+M6SZNjqMM2SddL+kneLy9ejf0h6e35/+ROSddKmhrX/pD0WUnHJd1ZKeu7D5R8PH9v75D0whHX44P5/+YOSV+WtK0y76pcj7slvXKQba96UFjF6ySawDsj4jnAhcBb8navBG6JiD3ALXl6HK4ADlemPwB8JNfjBHD5GOrwMeBrEfFHwPNzfca6PyTtAt4K7I2IC0hnry5jfPvjc8BFPWWL7YOLgT357wDwqRHX42bggoh4HumU/1UA+Xt7GfDc/JlP5t/VykTEqv4BLwa+Xpm+CrhqFepxA/AXwN3Azly2E7h7DNveTfqyvRy4iTQi9DdAo98+GlEdtgI/I/czVcrHuj/oDpHfQTo7dhPwynHuD+Bc4M6l9gHwL8Dr+y03inr0zPtr4Jr8vvjNkIYCvHil2131TIE1cJ2EpHOBFwC3AWdHxDGA/HrWGKrwUeBdQCtPnwk8HBHNPD2OffIs4EHgX3Mz5tOStjDm/RERvwQ+BNwHHAMeAW5n/PujarF9sJrf3TcDXx1FPdZCUFjyOomRblx6EvBF4G0R8ei4tlvZ/quB4xFxe7W4z6Kj3icN4IXApyLiBcDjjK/p1JHb65cA5wHPALaQ0vRea+G02ap8dyW9l9T8vWYU9VgLQWHJ6yRGRdIEKSBcExFfysW/lrQzz98JHB9xNV4KvEbSz0mXmL+clDlsy8PFYTz75ChwNCJuy9PXk4LEuPfHK4CfRcSDETEHfAl4CePfH1WL7YOxf3cl7QdeDbwhclth2PVYC0Hhe8Ce3Ls8SeowuXHUG5Uk4DPA4Yj4cGXWjcD+/H4/qa9hZCLiqojYHRHnkv7t34yINwC3Aq8dYz1+Bdwv6dm5aB/wY8a8P0jNhgslbc7/R+16jHV/9FhsH9wIvDGfhbgQeKTdzBgFSRcB7wZeExHTPfW7TNImSeeROj6/u+INjbLT6DQ6VF5F6k39P+C9Y9rmn5JSrDuAH+a/V5Ha87cA9+TXHWPcDy8Dbsrvn5X/Y48A/wFsGsP2/xg4lPfJfwLbV2N/AP8E/AS4E/h3YNO49gdwLakvY450BL58sX1ASts/kb+3PyKdMRllPY6Q+g7a39d/riz/3lyPu4GLB9m2RzSaWWEtNB/MbA1xUDCzgoOCmRUcFMys4KBgZgUHBTMrOCiYWcFBwcwK/w+7BRQCow9PeAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD8CAYAAAB+fLH0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFkpJREFUeJzt3X2MXFd5x/Hvb2Z2vbHB2A4kGDslQXIpIYKCLBqgrRCmJaGIpBJUoahYEMmqREt4kSCBP1D/A4F4k4DWAkpaRSE0QBNFvDSEINpKBBxAIcGEuAESE4MDcV7Ixrs7O0//OGdm7pmd9a53XvaF30dazdxz79x7fD3z3Oece+69igjMzNpqq10BM1tbHBTMrOCgYGYFBwUzKzgomFnBQcHMCg4KZlZwUDCzgoOCmRUaq10BgEltiim2rHY1zDa0xzjxm4h42lLLrYmgMMUW/kT7VrsaZhvaN+L6XyxnOTcfzKzgoGBmBQcFMys4KJhZwUHBzAoOCmZWcFAws8KSQUHSZyUdl3RnpeyDkn4i6Q5JX5a0rTLvKklHJN0t6ZWjqriZjcZyMoXPARf1lN0MXBARzwN+ClwFIOl84DLgufkzn5RUH1ptzWzklgwKEfFt4KGesv+KiGae/A6wO7+/BPh8RMxExM+AI8CLhlhfMxuxYfQpvBn4an6/C7i/Mu9oLjOzdWKgax8kvRdoAte0i/os1vce8pIOAAcAptg8SDXMbIhWHBQk7QdeDeyL7sMjjgLnVBbbDTzQ7/MRcRA4CLBVO/zwCbM1YkXNB0kXAe8GXhMR05VZNwKXSdok6TxgD/DdwatpZuOyZKYg6VrgZcBTJR0F3kc627AJuFkSwHci4u8j4i5JXwB+TGpWvCUi5kdVeTMbPq2Fx8Zt1Y7w/RTMRusbcf3tEbF3qeU8otHMCg4KZlZwUDCzgoOCmRUcFMys4KBgZgUHBTMrOCiYWcFBwcwKDgpmVnBQMLOCg4KZFRwUzKzgoGBmBQcFMys4KJhZwUHBzAoOCmZWcFAws4KDgpkVHBTMrOCgYGYFBwUzKzgomFnBQcHMCksGBUmflXRc0p2Vsh2SbpZ0T37dnssl6eOSjki6Q9ILR1l5Mxu+5WQKnwMu6im7ErglIvYAt+RpgItJD5XdQ3rM/KeGU00zG5clg0JEfBt4qKf4EuDq/P5q4NJK+b9F8h1gm6Sdw6qsmY3eSvsUzo6IYwD59axcvgu4v7Lc0VxmZuvEko+iP03qU9b3sdaSDpCaGEyxecjVMLOVWmmm8Ot2syC/Hs/lR4FzKsvtBh7ot4KIOBgReyNi7wSbVlgNMxu2lQaFG4H9+f1+4IZK+RvzWYgLgUfazQwzWx+WbD5IuhZ4GfBUSUeB9wHvB74g6XLgPuB1efGvAK8CjgDTwJtGUGczG6Elg0JEvH6RWfv6LBvAWwatlJmtHo9oNLOCg4KZFRwUzKzgoGBmBQcFMys4KJhZwUHBzAoOCmZWcFAws4KDgpkVHBTMrOCgYGYFBwUzKzgomFnBQcHMCg4KZlZwUDCzgoOCmRUcFMys4KBgZgUHBTMrOCiYWcFBwcwKDgpmVnBQMLPCQEFB0tsl3SXpTknXSpqSdJ6k2yTdI+k6SZPDqqytIdLgf7YmrTgoSNoFvBXYGxEXAHXgMuADwEciYg9wArh8GBW1VTDqH7QDxZo0aPOhAZwhqQFsBo4BLweuz/OvBi4dcBtmNkYrDgoR8UvgQ6SnTh8DHgFuBx6OiGZe7Ciwa9BK2hislaP2WqnH77FBmg/bgUuA84BnAFuAi/ssGot8/oCkQ5IOzTGz0mqY2ZAt+Sj6U3gF8LOIeBBA0peAlwDbJDVytrAbeKDfhyPiIHAQYKt29A0cNmKnexTWiE5WRWuJ7eZ6hr8m4zDI//J9wIWSNksSsA/4MXAr8Nq8zH7ghsGqaEO3VFquWv8/QDUN/HfK7Q1SbxuKQfoUbiN1KH4f+FFe10Hg3cA7JB0BzgQ+M4R6mtmYDNJ8ICLeB7yvp/he4EWDrNdGYKkjbJ+j9KJH9UFEa8F6o1VpFvSuv1/TQnJTYoQ8otHMCgNlCrYOLNV3sKBIi8+vzNMy2/ax4Iheh1ZZplo3G4ieeai2eLaQNrCsetjyOVMws4IzhY1qmRnCgn4D1YqMIK0qT9dq1cLlVaP3SN5qEb2HolY3G2jXp28/gzOGsXBQ2GgW+7EuFgja5blMUvfH3xMM0rw+QaRX9cebf9ydZoTUDRSttFzUSIGhWG2e19ucsJFz88HMCs4UNpJ+WcIymwqdJkK93s0M6rXucgD1Wv/MIlvQqdgKaM2n5dpH/Pn5ImsA0Px8t0mxoBNS3WxhqWaEmxBD4UzBzArOFDaCJTKENKmFZfVat/+gXu+W5ffUclmj3l2m3tPfUNl2p6+g/TrfgvmUKURzPq9TqF1WrUu7rJMxdLOCBZ2Pi52mtKFwUFjPlhsM2uW9KX+tGwDUyF+FRqMbBNplE+k1JhqQ50V94bY1X3YgMtdEORionn/0c3PdZgPpCvuATiBRtfMRijMTNh5uPphZwZnCerXEOIEFnYo1lZ2JgOr1ThbQzgo0MQGTEwDEVLq9Zkymea2pBq1GOo5Eo92M6HYiqplea7M5O5htwmzKBjQz116822zoZAdB9I6FyM0JGz9nCmZWcKawkfQbSNS+D0JlUJLaHYkTjW6GMJmzgqlJYvMmAOY3T+bXtExzc53mVO70q/QptPsSGifTa/2JdJRvPFGnNp222a6ZIjodkp2Ow3pUOifnu/UFov+Nu2yEHBTWm+WORegZT1DtVOy8NhrdYHBGCgStLVPMPymVzT4lNSNmnpLWNbtVNDen9bXyN0cB9ZN5dY+nH/Cmx9LyrcdqTLSrUzkj0e5M7JxxaHWbNnGKZlHfIdA2dG4+mFnBmcIGUoxF6D39KHVHKFY6FTudie0mw9ZJZranDGH6zJRRnHxaWsfMjqC5Nd+oe6J92rFG/bG03KYTOYuYaNeoTm2+7Hysz9Zhrs9Yh3ZG09u/WNPCMhspZwpmVnCmsF4sY6BSX9XOxXrPCMWJBrEpdyZuyf0IT57giR1p/vTT0zaf2J2yg607H+OZ208AsLkxC8BvT27h/t9uS8tNbknrz6MRGyehkTsaGxN50FOj1v9SbFsz/L9iZgVnChtAMVCpT19Cp7xWDl+OiQYxlU83npHmzW6tcXJH+szJp6fG/NP+IGUHr3jG3bz0yT8FYItSpnDXzC5ubpwPwB0n08PAmo/k/olN0MqnH6LWrU/7/bJu01I50+CzDuPhoLAeLLfpUC1rp+aVC546HY3t14kGrZzWz0+l17kzRPNJeXXb0w//uWf+CoC/3HonfzaVmhL1vK3NtXs5vPkZABzedDYArUYKClGtduViqc5l1K2FN2Ox1efmg5kVnClsML3NhmK6p4Mv6uqMTGyn+a0JmJ/MIxQnUvNhSz1lDJtrM9SVzjfO5ysXp1ubeHw+dVI251K2MTmb1lWbC2pz+fqG+cprO0OoZA/Rc9n1wrtA27g4UzCzwkCZgqRtwKeBC0gXv70ZuBu4DjgX+DnwNxFxYqBa2soVnZA5Q6j1yR4yBdTy+KTZmfT1ODqdTjl+74lnMd36JQCPxxQA//vYH3LXb5+e1nsiZQyN36XPTzwR1E+mbKM2l0cgzTXTzVeAyK/Mz3ezh1bPvRP63Uuhb5kzi2EZtPnwMeBrEfFaSZPAZuA9wC0R8X5JVwJXkp4vaUPW7/LoZX2ucofl9vt2ml+fgcZ0Ws9s/pH/ZOosAB5vTvLfU3sAmG6meb84sZ1Hjz0ZgDN+lZoPUw+ldU0+2qKRL47SyXwJ9VwTmjnq5Ps30oqFzQafdVg1K24+SNoK/Dn5AbIRMRsRDwOXAFfnxa4GLh20kmY2PoNkCs8CHgT+VdLzgduBK4CzI+IYQEQck3TW4NW0oejp4FOzheba1ySkeRPTNSYfzovl5sbMbDpHeeShMzhSueYBoPFonS0Ppcxi6sG03jMeSuucfGyO+uOpk1Izs/lzze79GttXSVabD6fiZsNYDNLR2ABeCHwqIl4APE5qKiyLpAOSDkk6NMfMANUws2EaJFM4ChyNiNvy9PWkoPBrSTtzlrATON7vwxFxEDgIsFU7HO6HoRXQe0PVVuVUX8+TmWjOd69enE7t/Im6gNQ3UEt3UGPid2md81ONzoCkdmdkYzqYfDStb9Mj6XXykfTB+u9m0XQK+DqZMoWYm+v0KXQ6GqunJDv19ROiVsuKM4WI+BVwv6Rn56J9wI+BG4H9uWw/cMNANTSzsRr07MM/AtfkMw/3Am8iBZovSLocuA943YDbsEW0j6ILzkLAgt58taJzZO7cdr3WhJl0XKi3BzsF1ObSco0nUsaw1C3YGtNltlGfzv0I0zPdDGE29yk0m92bsrZfW61TZwjuSxirgYJCRPwQ2Ntn1r5B1msr1wkC7YJWZSxA+3Zmze5ox/aIx849FFstarP57s05KEw2Khc1tZ/H0kzrrTVb1NqnG9t3bs6BgNm51FyA7mnIuWbqWKzULSIWf7aDn/kwdh7RaGYFX/uwUbUzhPYoxojOjVLp3Cm5m1Go+lSnnCnUKtdI9Oo8Dao5nwYk5c+mstyR2JyvvM/zKqcfVzxQyU2HkXKmYGYFZwobSbQ6D2ZtP4uxkwFInScotJ/hSES3s3I+36qtWe8e8XuvtMyfKV77PUS21e1ILK5vAGi1umWVei/IENy5uGocFNaraHVuqlKchWj/mHqDQ+UxbJ3g0Ir0IBYg2j/kWn3BNRTVoNB3PEFnvEHPhU6tSsCofq7nB+8zDWuLmw9mVnCmsB5E9L8lW/uoWskYFoxZaF8RWatkC+0xDFL3OojODVgqD1lor7ffNivrj96RkpXpxUYqtuu76Ho75c4Sxs2ZgpkVnClsMN3+hd4jeq3Tv1A8xLWdNbSXW+IR990NVY7gvacYO+Ur6D/ot34bKweF9aL9I1lGMyIV9QyBrpyZ6CxTo3tWgPbqlxcUigBwijMHfccdOBisaW4+mFnBmcJ6s1TG0PM8iOqRekGTYp6Fy3OaR+o+R/3Tzg7AGcIa4kzBzArOFNarpU5TtvXpZ+jO0tJH8KWqcaprFZwdrEvOFMys4ExhPTtV/0JnmcrR+hT9DYPV4zSzDWcIa5qDwkawnOAAq3vDEgeCdcPNBzMrOFPYSPodjZc7QnHU9bB1w5mCmRWcKWx0o8oenA1sWA4Kv4/8g7ZTcPPBzAoOCmZWcFAws8LAQUFSXdIPJN2Up8+TdJukeyRdlx8pZ2brxDAyhSuAw5XpDwAfiYg9wAng8iFsw8zGZKCgIGk38FfAp/O0gJeTHksPcDVw6SDbMLPxGjRT+CjwLqA9qP5M4OGIyE8T4Siwa8BtmNkYrTgoSHo1cDwibq8W91m070lxSQckHZJ0aI6ZlVbDzIZskMFLLwVeI+lVwBSwlZQ5bJPUyNnCbuCBfh+OiIPAQYCt2uHRNGZrxIozhYi4KiJ2R8S5wGXANyPiDcCtwGvzYvuBGwaupZmNzSjGKbwbeIekI6Q+hs+MYBtmNiJDufYhIr4FfCu/vxd40TDWa2bj5xGNZlZwUDCzgoOCmRUcFMys4KBgZgUHBTMrOCiYWcFBwcwKDgpmVnBQMLOCg4KZFRwUzKzgoGBmBQcFMys4KJhZwUHBzAoOCmZWcFAws4KDgpkVHBTMrOCgYGYFBwUzKzgomFnBQcHMCg4KZlYY5KnT50i6VdJhSXdJuiKX75B0s6R78uv24VXXzEZtkEyhCbwzIp4DXAi8RdL5wJXALRGxB7glT5vZOjHIU6ePRcT38/vHgMPALuAS4Oq82NXApYNW0szGZyh9CpLOBV4A3AacHRHHIAUO4KxhbMPMxmPgoCDpScAXgbdFxKOn8bkDkg5JOjTHzKDVMLMhGSgoSJogBYRrIuJLufjXknbm+TuB4/0+GxEHI2JvROydYNMg1TCzIRrk7IOAzwCHI+LDlVk3Avvz+/3ADSuvnpmNW2OAz74U+DvgR5J+mMveA7wf+IKky4H7gNcNVkUzG6cVB4WI+B9Ai8zet9L1mtnq8ohGMys4KJhZwUHBzAoOCmZWcFAws4KDgpkVHBTMrOCgYGYFBwUzKzgomFnBQcHMCg4KZlZwUDCzgoOCmRUcFMys4KBgZgUHBTMrOCiYWcFBwcwKDgpmVnBQMLOCg4KZFRwUzKzgoGBmBQcFMyuMLChIukjS3ZKOSLpyVNsxs+EaSVCQVAc+AVwMnA+8XtL5o9iWmQ3XqDKFFwFHIuLeiJgFPg9cMqJtmdkQjSoo7ALur0wfzWVmtsYN8ij6U+n3NOooFpAOAAcAptg8omqY2ekaVVA4CpxTmd4NPFBdICIOAgcBJD34jbj+ceA3I6rP6XgqrkeV61Faz/V45nIWUkQsvdRpktQAfgrsA34JfA/424i46xSfORQRe4demdPkergev+/1GEmmEBFNSf8AfB2oA589VUAws7VjVM0HIuIrwFdGtX4zG421NKLx4GpXIHM9Sq5HacPXYyR9Cma2fq2lTMHM1oA1ERRW4zoJSedIulXSYUl3Sboil++QdLOke/Lr9jHVpy7pB5JuytPnSbot1+M6SZNjqMM2SddL+kneLy9ejf0h6e35/+ROSddKmhrX/pD0WUnHJd1ZKeu7D5R8PH9v75D0whHX44P5/+YOSV+WtK0y76pcj7slvXKQba96UFjF6ySawDsj4jnAhcBb8navBG6JiD3ALXl6HK4ADlemPwB8JNfjBHD5GOrwMeBrEfFHwPNzfca6PyTtAt4K7I2IC0hnry5jfPvjc8BFPWWL7YOLgT357wDwqRHX42bggoh4HumU/1UA+Xt7GfDc/JlP5t/VykTEqv4BLwa+Xpm+CrhqFepxA/AXwN3Azly2E7h7DNveTfqyvRy4iTQi9DdAo98+GlEdtgI/I/czVcrHuj/oDpHfQTo7dhPwynHuD+Bc4M6l9gHwL8Dr+y03inr0zPtr4Jr8vvjNkIYCvHil2131TIE1cJ2EpHOBFwC3AWdHxDGA/HrWGKrwUeBdQCtPnwk8HBHNPD2OffIs4EHgX3Mz5tOStjDm/RERvwQ+BNwHHAMeAW5n/PujarF9sJrf3TcDXx1FPdZCUFjyOomRblx6EvBF4G0R8ei4tlvZ/quB4xFxe7W4z6Kj3icN4IXApyLiBcDjjK/p1JHb65cA5wHPALaQ0vRea+G02ap8dyW9l9T8vWYU9VgLQWHJ6yRGRdIEKSBcExFfysW/lrQzz98JHB9xNV4KvEbSz0mXmL+clDlsy8PFYTz75ChwNCJuy9PXk4LEuPfHK4CfRcSDETEHfAl4CePfH1WL7YOxf3cl7QdeDbwhclth2PVYC0Hhe8Ce3Ls8SeowuXHUG5Uk4DPA4Yj4cGXWjcD+/H4/qa9hZCLiqojYHRHnkv7t34yINwC3Aq8dYz1+Bdwv6dm5aB/wY8a8P0jNhgslbc7/R+16jHV/9FhsH9wIvDGfhbgQeKTdzBgFSRcB7wZeExHTPfW7TNImSeeROj6/u+INjbLT6DQ6VF5F6k39P+C9Y9rmn5JSrDuAH+a/V5Ha87cA9+TXHWPcDy8Dbsrvn5X/Y48A/wFsGsP2/xg4lPfJfwLbV2N/AP8E/AS4E/h3YNO49gdwLakvY450BL58sX1ASts/kb+3PyKdMRllPY6Q+g7a39d/riz/3lyPu4GLB9m2RzSaWWEtNB/MbA1xUDCzgoOCmRUcFMys4KBgZgUHBTMrOCiYWcFBwcwK/w+7BRQCow9PeAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using true PSF\n",
      "0.1974998414516449 0.2962857484817505\n",
      "\n",
      "Using bad PSF\n",
      "0.1992892026901245 0.2989780604839325\n",
      "\n",
      "(delta e)/e\n",
      "0.009060064176900692 0.009086876489936328\n"
     ]
    }
   ],
   "source": [
    "# Now let's draw with one and fit with the other.\n",
    "gal = galsim.Sersic(half_light_radius=0.8, n=1.5).shear(e1=0.2, e2=0.3)\n",
    "conv = galsim.Convolve(psf2, gal)\n",
    "shapeTruePSF = getGalShape(conv, psf2)\n",
    "shapeBadPSF = getGalShape(conv, matchedPSF)\n",
    "print(\"Using true PSF\")\n",
    "print(shapeTruePSF.corrected_e1, shapeTruePSF.corrected_e2)\n",
    "print()\n",
    "print(\"Using bad PSF\")\n",
    "print(shapeBadPSF.corrected_e1, shapeBadPSF.corrected_e2)\n",
    "print()\n",
    "print(\"(delta e)/e\")\n",
    "print((shapeBadPSF.corrected_e1-shapeTruePSF.corrected_e1)/shapeTruePSF.corrected_e1, \n",
    "      (shapeBadPSF.corrected_e2-shapeTruePSF.corrected_e2)/shapeTruePSF.corrected_e2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:27.965350Z",
     "start_time": "2019-01-22T22:52:27.961609Z"
    }
   },
   "outputs": [],
   "source": [
    "def plotPSFProfile(psf, ax, **kwargs):\n",
    "    xs = np.linspace(-2, 2, 400)\n",
    "    ax.plot(xs, [psf.xValue(0, x) for x in xs], **kwargs)\n",
    "def plotPSFDiff(psf1, psf2, ax, **kwargs):\n",
    "    xs = np.linspace(-2, 2, 400)\n",
    "    ax.plot(xs, [psf1.xValue(0, x)-psf2.xValue(0, x) for x in xs], **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:28.088039Z",
     "start_time": "2019-01-22T22:52:27.966793Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x813463080>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOW9x/HPbyb7QkhIgLAm7ARIICSA7JssLuBaQFtX5Fqr9uqtVrt4rd1ua1u72VZqW63VgruAKArKpixhCxDWgAHCGiBkX2Z57h8T0hASMglJTjL5vXVeM+fMmXO+c5L8eObMc54jxhiUUkr5FpvVAZRSSjU+Le5KKeWDtLgrpZQP0uKulFI+SIu7Ukr5IC3uSinlg7S4K6WUD9LirpRSPkiLu1JK+SA/qzYcHR1t4uLirNq8Ukq1Slu3bj1rjImpaznLintcXBxbtmyxavNKKdUqicgRb5bTwzJKKeWDtLgrpZQP0uKulFI+yLJj7kqp1sHhcJCdnU1paanVUdqUoKAgunXrhr+/f4Ner8VdKXVF2dnZhIeHExcXh4hYHadNMMZw7tw5srOziY+Pb9A69LCMUuqKSktL6dChgxb2ZiQidOjQ4ao+LWlxV0rVSQt787vafa6HZZS6KDcLDnwCRTkQ1hF6T4YOva1OpVSDaHFXKu84LP8O7P8IqHZN4SFfg6nPQkRXC4Kpi8LCwigsLARg+fLlfPvb32bVqlX06NGjxuWfffZZwsLC+M53vtOcMVsULe6qbctcBW/fB65ymPBdGDoP2veEC0dg6yuw8c9w6DOYtwi6p1qdts1btWoVjzzyCJ988kmthb0lcTqd+PlZU2b1mLtquzJXwb/nQUQ3eHA9THoaIuNAxHM/9VnP/MBwePUGOPKltXnbuHXr1vHAAw/w4Ycf0ru353DZkSNHmDJlComJiUyZMoWjR49e9rqJEyfy2GOPMX78eAYOHEhaWhq33HILffv25Qc/+EHlcr/5zW8YPHgwgwcP5re//W3l/B//+McMGDCAa6+9lnnz5vGrX/0KgB07djBq1CgSExO5+eabyc3Nrdze9773PSZMmMDvfve7GjPm5eURFxeH2+0GoLi4mO7du+NwOBptf2nLXbVNpzNg0Z0Q3Q/uXgIhUZQ73Rw8U0BesYN2wf4MjG2HPbovzF8Jf58Bi+6A+1dCdB+r01vmR0sz2HMiv1HXmdClHf9746ArLlNWVsbs2bNZvXo1AwYMqJz/8MMPc9ddd3H33Xfz97//nUcffZT333//stcHBASwdu1afve73zF79my2bt1KVFQUvXv35rHHHiMrK4t//OMfbNq0CWMMI0eOZMKECbhcLt555x22b9+O0+kkOTmZ4cOHA3DXXXfxhz/8gQkTJvDMM8/wox/9qPIfhQsXLrBmzRoAbrzxxhozJiUlsWbNGiZNmsTSpUuZPn16g/u010Rb7qrtKc2Hxd+AoAj4+jtcIIznV+xj+I8/5frfr+eOlzdxwx/WM+rnq3h+xT6K/dvDnW+C2GDx18GhJ/M0N39/f0aPHs3f/va3S+Zv2LCBO+64A4BvfOMbrF+/vsbXz5o1C4AhQ4YwaNAgYmNjCQwMpFevXhw7doz169dz8803ExoaSlhYGLfccgvr1q1j/fr1zJ49m+DgYMLDw7nxxhsByMvL48KFC0yYMAGAu+++m7Vr11Zub86cOXVmnDNnDosXLwZg0aJFl7ymMWjLXbU9Hz/t6Rlz91LSzvnzrdfXklNYxnWDY5k5pDMdQgM5U1DKsp0nefHzQ7y//QR/uGMYyTcvhNdvhc9+DNN/avW7sERdLeymYrPZePPNN5k6dSo/+9nP+N73vlfjcrV1HwwMDKxcz8XHF6edTifGmBpfV9v8uoSGhtb63MWMs2bN4umnn+b8+fNs3bqVyZMnN2hbtdGWu2pbvloLO/4FYx7lo4JezFu4kZAAO0sfHsuLdyZzQ2IXrundgdlDu/LXu1J487+uwW4T7vjrRj53J0LK/bDhRTi22ep30uaEhISwbNkyXn/99coW/OjRo1m0aBEAr7/+OmPHjm3QusePH8/7779PcXExRUVFvPfee4wbN46xY8eydOlSSktLKSws5MMPPwQgIiKCyMhI1q1bB8Brr71W2YqvrraMYWFhjBgxgm9/+9vccMMN2O32BmWvjbbcVdvhLIdlj0FkHGti7+PRN7aT1L09/7g3lXZBNR/rHBEfxTvfHM09/9jMf/1zK/++5zGG7/8IPnoS5n8GNm0fNaeoqCg+/vhjxo8fT3R0NL///e+57777eP7554mJieEf//hHg9abnJzMPffcw4gRIwCYP38+w4YNAzwt7KSkJHr27ElKSgoREREAvPrqqzz44IMUFxfTq1evWrd9pYxz5szh9ttvZ/Xq1Q3KfSXS0I8dVyslJcXoxTpUs9r0Enz0JCev/yfXLguiZ4cQ/r1gVK2FvaoLxeXc/KcvyStx8OmUU3T45GGY/SIM+3ozBLfW3r17GThwoNUxLFNYWEhYWBjFxcWMHz+ehQsXkpyc3Czbrmnfi8hWY0xKXa/VZodqG8oKYM0vcfUcx73r2+NvF/56V4pXhR2gfUgAf78nFafLzfzt8ZiuqfDZT/TL1TZgwYIFDB06lOTkZG699dZmK+xXSw/LqLZhw5+g+Cyvhd7Dvv2F/OPeVLq0D67XKuKjQ/npzUN45N/beT/1Pm4+/k3Y9k8YuaCJQquW4I033rA6QoN41XIXkRkisl9EMkXkqRqef0FEdlTcDojIhcaPqlQDlRXCxj+R3/NantsezNzU7kzq37FBq7oxqQuzkrrwxNYISmJHwvrfaOtdtUh1FncRsQMvAjOBBGCeiCRUXcYY85gxZqgxZijwB+DdpgirVINsfQVKL/DjCzPoEBbI0zOv7vjx/96YQEiAH79x3AoFJyG9dbbslG/zpuU+Asg0xhw2xpQDi4DZV1h+HvDvxgin1FVzlsOGP5ITPZK3Tsfy9MwBRIRc3VmAHcICefzafvw1uyt5kYM9h3wqTiNXqqXwprh3BY5Vmc6umHcZEekJxAOf1fL8AhHZIiJbcnJy6ptVqfrbuwQKTvJ/edcyuGs7bhraOKM7fn1UT/p2DOcPJdPh3EE4+EmjrFepxuJNca/plK/a+k/OBd42xrhqetIYs9AYk2KMSYmJifE2o1INl/Yy+cHdebdgAN+bORCbrXEuOuFnt/E/0/rxyoUkioM6w8YXG2W96nITJ05kxYoVl8z77W9/y0MPPdSgdV3sgp2VlUXfvn0vW7ev8Ka4ZwPdq0x3A07Usuxc9JCMailO7YajG/h72WRG9opmdJ/oRl399EGdGdi1A685p3jOfD13qFHXrzzmzZtXeYbnRYsWLWLevHkNXmd2djbTp0/n17/+NdOnT/fqNS5XjW3WFsub4p4G9BWReBEJwFPAl1RfSET6A5HAhsaNqFQDbfkbTlsgrxSP5tEpfRt99SLCf0/ty8uFo3GLHba92ujbUHDbbbexbNkyysrKAE+L+8SJE4wdO5YnnniCwYMHM2TIkMpBuFavXs3EiRO57bbbGDBgAHfeeeclY8ScOnWKadOm8ZOf/KRyQLGsrCzGjRtHcnIyycnJfPnll5XrmjRpEnfccQdDhgwhKyuLAQMGMH/+fAYPHsydd97JypUrGTNmDH379mXzZs+wFJs3b2b06NEMGzaM0aNHs3//fgBeeeUVbrnlFmbMmEHfvn158sknm2y/1dnP3RjjFJGHgRWAHfi7MSZDRJ4DthhjLhb6ecAiY9Upr0pVVZqPSV/MChlDn549uKZXhybZzKT+HYno2J2NJalcs/11ZNIPwC+gSbbVInz0FJza1bjr7DwEZv5frU936NCBESNG8PHHHzN79uzKERTfffddduzYQXp6OmfPniU1NZXx48cDsH37djIyMujSpQtjxozhiy++qBzT5a677uInP/kJt99+e+U2OnbsyKeffkpQUBAHDx5k3rx5lYdvNm/ezO7du4mPjycrK4vMzEzeeustFi5cSGpqKm+88Qbr169nyZIl/OxnP+P9999nwIABrF27Fj8/P1auXMn3vvc93nnnHcAzDvz27dsJDAykf//+PPLII3Tv3p3G5lU/d2PMcmNMP2NMb2PMTyvmPVOlsGOMedYYc1kfeKUssetNxFHEwuKJLBjfq8ku8GyzCQvG9WJh0Xik+CzsX94k22nrqh6auXhIZv369cybNw+73U6nTp2YMGECaWlpAIwYMYJu3bphs9kYOnQoWVlZleuaOnUqr732GsXFxZXzHA4HDzzwAEOGDOH2229nz549lc+NGDGC+Pj4yun4+HiGDBmCzWZj0KBBTJkyBRGpbNmDZ0jg22+/ncGDB/PYY4+RkZFR+fopU6YQERFBUFAQCQkJHDlypCl2mZ6hqnxU+iKO+sVxPngwUwZ2atJNzR7WhV99nMJZOhK99RUYdFOTbs9SV2hhN6WbbrqJxx9/nG3btlFSUkJycjL//Oc/a12+6rC+drsdp9NZOf3kk0/yr3/9i9tvv50PPvgAPz8/XnjhBTp16kR6ejput5ugoKDK5asP31t9yOCqwwlf3M4Pf/hDJk2axHvvvUdWVhYTJ070Kltj0rFllO85mwnZafyr5BruHh2PvZF6yNQm0M/O3WN788/S8XD4c7hw+aXe1NUJCwtj4sSJ3HfffZVfpI4fP57FixfjcrnIyclh7dq1laM61uWFF16gXbt23H///RhjyMvLIzY2FpvNxmuvvXbVX57m5eXRtaun2+0rr7xyVetqKC3uyvfsXIQbG5/Yx/O11MY/llmTO0f24EPxHO9l19vNss22Zt68eaSnpzN37lwAbr75ZhITE0lKSmLy5Mn88pe/pHPnzl6tS0R49dVXOXnyJE8++SQPPfQQr776KqNGjeLAgQNXvNiGN5588kmefvppxowZY1kvGx3yV/kWtxvXbxP54kIUq4b/iR/NHtxsm3588Q6+sWcBSR1t2L61sdm229Ta+pC/VtIhf5W66OgG7PnHeNs5lm9c07NZNz1vZA/ecY7GlrPX08deKQtpcVc+xaT/m2KCyek2lT4dw5t12yk9I9kbOQUndtj1VrNuW6nqtLgr3+Fy4MpYwseu4dyU2qfZNy8iXD9qMGtciTh2vKmDiSlLaXFXvuPwGvzK81gpo7k+sYslEW5J7sqHjMG/6AQc1ZO1lXW0uCuf4dj1DgUmmHaDphEWaM0pHO1DApAB11NCAK7d71mSQSnQ4q58hbMcs3cZn7iHc9vI3pZGuX54b9a4knBkLNFDM8oyWtyVb/hqDQGOfNJCJjC8Z6SlUcb1jWG93zUElZyG41stzeKLnn32WX71q1/xzDPPsHLlSgDWrVvHoEGDGDp0KCUlJTzxxBMMGjSIJ554wuK01tHhB5RPKNnxNk4TTOywmU02joy3/O02woZcT3n6i5hd7xPYPdXSPL7queeeq3z8+uuv853vfId7770XgJdeeomcnJxLTvW/EqfTiZ+fb5VDbbmr1s/lwLb/Qz5xD+f65Dir0wAwM3UAX7oH49j9PuhAqVftpz/9Kf3792fq1KmVw+fec889vP3227z88su8+eabPPfcc9x5553MmjWLoqIiRo4cyeLFi8nJyeHWW28lNTWV1NRUvvjiC8DzCWDBggVMmzaNu+66C5fLxRNPPEFqaiqJiYm89NJLwJWHEE5LS2P06NEkJSUxYsQICgoKal1Pc/Otf6pU23TkCwKdBeyJmMCtzdy3vTaJ3SL4TchYJha/CKd3e4a19QG/2PwL9p3f16jrHBA1gO+O+G6tz2/dupVFixaxfft2nE4nycnJDB8+vPL5+fPns379em644QZuu+02wDMWzY4dOwC44447eOyxxxg7dixHjx5l+vTp7N27t3Ld69evJzg4mIULFxIREUFaWhplZWWMGTOGadOmATUPITxixAjmzJnD4sWLSU1NJT8/n+DgYP72t7/VuJ6qI0s2By3uqtXL2/EBgcafrsOvszpKJRGh3dDZuDb+iZId7xE2wzeKuxXWrVvHzTffTEhICEDlBTa8tXLlykuG8M3Pz6egoKByXcHBwQB88skn7Ny5k7ff9owNlJeXx8GDBwkICKgcQhioHEI4IiKC2NhYUlM9h93atWt3xfVocVeqPoxB9i9nvXsIM5Ot7SVT3aThg0jbMID+u96HGc9YHadRXKmF3ZSu5nsUt9vNhg0bKot4VVUHCDPG8Ic//OGyy+6tXr26xmF6jTE15qptPc3Nq2PuIjJDRPaLSKaI1HhBDhH5mojsEZEMEXmjcWMqVYtTu2hXdoqDUeOJjbj8j9dKfTqGsS1kLJFFh+DsQavjtFrjx4/nvffeo6SkhIKCApYuXVqv10+bNo0//vGPldMXD9dUN336dP785z/jcDgAOHDgAEVFRbWud8CAAZw4caLyAiEFBQU4nc56r6ep1NlyFxE78CJwLZ6LZaeJyBJjzJ4qy/QFngbGGGNyRaRjUwVWqqpz294j0ghRQ+v3Ub25BAyZBWl/pXDHe4RNbbrrZfqy5ORk5syZw9ChQ+nZsyfjxo2r1+t///vf861vfYvExEScTifjx4/nL3/5y2XLzZ8/n6ysLJKTkzHGEBMTw/vvv1/regMCAli8eDGPPPIIJSUlBAcHs3Llynqvp6nUOeSviFwDPGuMmV4x/TSAMebnVZb5JXDAGPOytxvWIX9VYzjz/AiOFhi6/c86OkcE1f2CZnbwdAFlL46lY1R7Ov73GqvjNIgO+Wudph7ytytwrMp0dsW8qvoB/UTkCxHZKCIzvFivUlfnwjE6Fu1nT/iYFlnYAfp2Cmd70CiiL+yEonNWx1FtiDfFvaZvMqo39/2AvsBEYB7wsoi0v2xFIgtEZIuIbMnJyalvVqUucX7bBwAEDW6Zh2Qusg+YgQ03+bv14tmq+XhT3LOBqtcq6wacqGGZD4wxDmPMV8B+PMX+EsaYhcaYFGNMSkxMTEMzKwVA8a4lHHLHcs3IUVZHuaKU0ZM5Y9pzfvsSq6M0mFVXbGvLrnafe1Pc04C+IhIvIgHAXKD6b+n7wCQAEYnGc5jm8FUlU+pKSvPpnLuF7SGj6R4VYnWaK+rXOYItAal0PL0eXA6r49RbUFAQ586d0wLfjIwxnDt3jqCghh9urLO3jDHGKSIPAysAO/B3Y0yGiDwHbDHGLKl4bpqI7AFcwBPGGD3AqJrM+d2fEoULe39r+xJ7y9HrWkL2f0re/jVEJEy1Ok69dOvWjezsbPRQavMKCgqqPHGqIbw6ickYsxxYXm3eM1UeG+DxiptSTe7M9g/xN8EkXTPN6ihe6TPqRsr2fZ/TaR+0uuLu7+/f7GdXqqunA4ep1scYok+uJd1/KL06Wzu8r7cS4mLZZhtMRPYqq6OoNkKLu2p18o/tJtqdQ3HPSVZH8ZqIcLbLJDo5jlN2qnEH3lKqJlrcVatzdLPn+/zY4ddbnKR+YlJmA3B0w7sWJ1FtgRZ31er4H17FIbqRMGCQ1VHqZdiQIRw03ZFDK62OotoALe6qVXGVFhJfnM7RqDHYbdZecam+Av3sHIm6hh6F6bhLC6yOo3ycFnfVqhxO+4gAnAQntI4ukNUFJ0wnACdfbfnY6ijKx2lxV61K/u6PKTaBJIxqncV98MgZFJlA8nZ9ZHUU5eO0uKtWJTbnC/YGDaVdWJjVURokol0Y+4KHEXtmvV5bVTUpLe6q1Tj1VQZd3Ccp6TnZ6ihXpSxuMrHmNNmZu6yOonyYFnfValzsAtl9xI0WJ7k6cSM9XSKPbW69A4mplk+Lu2o1gr5axTHpQo/eCVZHuSpd4gdwzNaVoKOfWx1F+TAt7qpVKCkqpG9JOidixlzVxZJbijMdxzGwNJ28vHyroygfpcVdtQr7Nq8gWMoJTfCNi3y1S5xJkDjYt0l7zaimocVdtQpFGR9TavzpO9I3inuv4dMoxZ+yfSusjqJ8lBZ31eIZY+h27gsOhQwlMLh1doGszh4YwuHQZHqc34DbrV0iVePT4q5avMwDe4gzxymPb91dIKtz9ZpMHCfYu2en1VGUD9Lirlq841uWAtBjZMu+EHZ99azoEnlq6zKLkyhf5FVxF5EZIrJfRDJF5Kkanr9HRHJEZEfFbX7jR1VtVfCRzzlt60SHHq1rFMi6tOs6gFP2zoRnr7Y6ivJBdRZ3EbEDLwIzgQRgnojU1NF4sTFmaMXt5UbOqdqo8/mFDCrbwamOY8EHukBeQoQzncYzuDydM+fzrE6jfIw3LfcRQKYx5rAxphxYBMxu2lhKeezeuIIwKSViyEyrozSJyMSZhEgZezZ9YnUU5WO8Ke5dgWNVprMr5lV3q4jsFJG3RaR7o6RTbV7Z3hU48KNHcuscBbIu3YZNoxw/HNolUjUyb4p7TZ+Fq/fdWgrEGWMSgZXAqzWuSGSBiGwRkS05OTn1S6raHIfLTc/cDWSFJmILbmd1nCYhgWEcCRtK3IWNlDvdVsdRPsSb4p4NVG2JdwNOVF3AGHPOGFNWMflXYHhNKzLGLDTGpBhjUmJiYhqSV7Uhu/bsoR9HcfWaYnWUJmV6T6GvHGNnxm6roygf4k1xTwP6iki8iAQAc4FLhrMTkdgqk7OAvY0XUbVVJ7d9CECPkb79FU/3EZ4unqe2Lbc4ifIldRZ3Y4wTeBhYgadov2mMyRCR50TkYsfjR0UkQ0TSgUeBe5oqsGo7wo6t5rw9mpCug62O0qSCuwzinD2GdsdXWx1F+RA/bxYyxiwHlleb90yVx08DTzduNNWWHc3JZ5hjBye7TifK17pAVifC2c7jGJr9EUfOXKBnx/ZWJ1I+QM9QVS3S7k0raSfFRCZeZ3WUZhGZeB3tpISMzausjqJ8hBZ31SI59n+CCxsdk6ZZHaVZdEyahhM7zgPa3101Di3uqsUpKnPSO38jx8MSIbiNHKIIiuB42BB6522kuNxpdRrlA7S4qxZn8+59DJavMH18uwtkdabPVAZJFmm79lkdRfkALe6qxcnZ7vnuvktK674Qdn11rXi/Z3Zol0h19bS4qxbFGEPE8TXk2yPx75JkdZxm5d81iTx7FBHH12CMXsBDXR0t7qpFycjOZYR7B7mx48DWxn49RTgfO45U1w72Hr9gdRrVyrWxvx7V0u3esoZIKSRq6PVWR7FEVNJMIqWQPVtXWx1FtXJa3FWL4j74KW6E8IS20QWyuohB03Fhw3XgU6ujqFZOi7tqMXIKyhhYuIkz4YMgJMrqONYIieJ0+CD6FWzifFG51WlUK6bFXbUYX+7cR5IcwtavbbbaL7L1mUqSHGLD7gNWR1GtmBZ31WKc3/kRNjHEJLetLpDVdRx2PTYx5Oz42OooqhXT4q5ahHKnm5hT6yjwi0Rih1odx1K2bskU2SOIOrkGp0sv4KEaRou7ahHSDucwhh0UdJ3Q9rpAVmezk9dlLNeYdLZmnbM6jWql2vhfkWop9m9dTaQU0mHYDVZHaREiE68jRvLI2L7e6iiqldLirloEv8MrcWMjsF/bGk+mNsEDK75UPrjS2iCq1dLirix3KKeQoWVpnG2f2Ha7QFYX1pGz4QMZVLKFY+eLrU6jWiGviruIzBCR/SKSKSJPXWG520TEiEhK40VUvu7LHRkk2r4icOAMq6O0KPZ+UxkuB1i/O9PqKKoVqrO4i4gdeBGYCSQA80QkoYblwvFcP3VTY4dUvq1gt6fLX0QbueqStyITr8NP3JzdqWerqvrzpuU+Asg0xhw2xpQDi4CaLkf/Y+CXQGkj5lM+Lq/YQc/zX1LoHw2dE62O07J0G0GpPYxOZ9brBTxUvXlT3LsCx6pMZ1fMqyQiw4DuxphlV1qRiCwQkS0isiUnJ6feYZXvWbPvBONsOymNmwy+fiHs+rL7UdhlLONkB18ePGt1GtXKeFPca/qLqxxsWkRswAvA/9S1ImPMQmNMijEmJSYmxvuUymcd3r6adlJMVFLbHAWyLu0TZxIr59mdrkc7Vf14U9yzge5VprsBJ6pMhwODgdUikgWMApbol6qqLg6Xm9Bjn+HCjq3PJKvjtEh+/a713B/6VC/goerFm+KeBvQVkXgRCQDmAksuPmmMyTPGRBtj4owxccBGYJYxZkuTJFY+Iy3rPGPc28mLToagCKvjtEwRXclt159URxp7TuZbnUa1InUWd2OME3gYWAHsBd40xmSIyHMiMqupAyrflZa+mwTbEcIGz7Q6SosWkHAdKbKfDbt0lEjlPa/6uRtjlhtj+hljehtjflox7xljzJIalp2orXblDef+FQAEaP/2KwodfCN2MRRm6CiRynt6hqqyxKGcQoYWb6AguAt0vOy0CVVVl2EU+XegT+56cgrKrE6jWgkt7soSa3ZnMda2G/pfp10g62KzUd77WibY0lm9J9vqNKqV0OKuLHFu5woCxUF4on5t4432Q2cRLiUc3a5nqyrvaHFXze5CcTlxZ9dQag+DnqOtjtMqSK9JOCSAmBOf69mqyita3FWzW7P/FJNs2ynuORns/lbHaR0CQijoMoaJbGXtfj27W9VNi7tqdpnbPida8mk/7Caro7QqEUNn0cOWQ/r2jVZHUa2AFnfVrEodLtofXek5K7XvVKvjtCr2/p7zAYK++kSvrarqpMVdNasvMs8ywWwhr9NIPSu1vtrFktd+EGNcW0jLyrU6jWrhtLirZrVlaxp9bCdoN7SmUaNVXYIH30CyHOSL9H1WR1EtnBZ31WycLjf+hzxnpfoN1AtzNETAoOuxicGx90MdSExdkRZ31Ww2Z51njGsz+REDoH0Pq+O0Tp0TKQzuyojSL9l7ssDqNKoF0+Kums369P2kyH6CB+vY7Q0mgi3hRsbadrFm5yGr06gWTIu7ahbGGJx7lmEXg/9gPd5+NUKSbiFQnOTv+tDqKKoF0+KumsXO7DxGl31BYUh3vVbq1eqWSnFANIPz15B1tsjqNKqF0uKumsXn6QcYY9uN3+CbdKCwq2WzYQZczyRbOit2HLY6jWqhtLirJmeMoXTXMvzFRVDSzVbH8QmhQ28hRMo4s+Mjq6OoFsqr4i4iM0Rkv4hkishTNTz/oIjsEpEdIrJeRHSAblUp80whKcVrKQzats7aAAAeHklEQVSKhS7JVsfxDT3HUOoXwaA8PTSjalZncRcRO/AiMBNIAObVULzfMMYMMcYMBX4J/KbRk6pWa+X2TMbZdiEJs/SQTGOx++PqN4Optm18lH7U6jSqBfKm5T4CyDTGHDbGlAOLgEu6Oxhjql65NxTQsysU4Dkkk5u+lEBxEjr0Vqvj+JTQpFtoJ8Wc3K6X31OX86a4dwWOVZnOrph3CRH5logcwtNyf7Rx4qnWbv/pApIL11IcGAPdUq2O41t6T6LML4zEvM84ck4PzahLeVPca/ocfVnL3BjzojGmN/Bd4Ac1rkhkgYhsEZEtOTk6JnVbsGLbISbadiAJs8Gm3983Kr9AnP1uYJotjY93ZFmdRrUw3vy1ZQPdq0x3A05cYflFQI0DdRtjFhpjUowxKTExMd6nVK2SMYb89A8IEgfBSbdYHccnhSZ/jXZSQs72ZVZHUS2MN8U9DegrIvEiEgDMBZZUXUBE+laZvB442HgRVWu1+3g+Y4o/pyioM/S4xuo4vil+AiX+kQzNW6WHZtQl6izuxhgn8DCwAtgLvGmMyRCR50Tk4tWNHxaRDBHZATwO3N1kiVWrsWprBuNtO7En3a6HZJqK3Q/XwNlMtW1jxbZMq9OoFsTPm4WMMcuB5dXmPVPl8bcbOZdq5YwxOHa9j5+48Rs21+o4Pi1s+FzY+Qq52z+Aa5OsjqNaCG1OqSax7egFJpZ/Tl54H+g0yOo4vq37SIqCOpFS8Bl7T+bXvbxqE7S4qyaxZtMWUm0HCEqeqycuNTWbDduQ2xhv28lHaXusTqNaCC3uqtGVO93Y974LQODQ2y1O0zYED7sdf3FRnv4uLreeQ6i0uKsm8Pn+M0xzrSMvOhki46yO0zbEDqUgvDfTHKvYePic1WlUC6DFXTW6TRvWMNB2jLAU/SK12YgQlPoNkm2ZrN/wpdVpVAugxV01qtyicnoceRen+GNP1EMyzcl/2Dxc2InKfJuScpfVcZTFtLirRrV8Rxazbespip8OIVFWx2lbwjuT13U8N7KWT/dc6SRy1RZocVeN6vimd4mUQiJG32t1lDap/eh76Cy5HPhiSd0LK5+mxV01mkM5hYzM/ZDCwM7Qa5LVcdokW/+ZlPi1o/+pJZy4UGJ1HGUhLe6q0Xzy5VbPRTmG3gE2u9Vx2ia/QJwJtzHNtpUPNmZYnUZZSIu7ahQOlxvS38AmhtBROrSQlcJH3U2gOCjeskj7vLdhWtxVo1iZcZLrnas433GU9m23WmwSF9oncH35x6w/qNdNaKu0uKtGkb5uCT1sObQfc7/VUZQIoWP+iwG2Y2xZ96HVaZRFtLirq3b0XDHJp96i2D8SW8Ksul+gmpx/0tcotYfR7+ibnC0sszqOsoAWd3XVPly3kSm2bbiH3QX+QVbHUQABIZQMmst02cRHG9OtTqMsoMVdXZVyp5uQ9FdAhLAxC6yOo6qIHP8gAeKidNMruPWL1TZHi7u6Kp/tPsIs9yrOdbsWIrpZHUdVFd2XnJhRXFf+MWv2nbI6jWpmXhV3EZkhIvtFJFNEnqrh+cdFZI+I7BSRVSLSs/Gjqpboq89fIVIK6TD5EaujqBpETniIrnKO9M8WWx1FNbM6i7uI2IEXgZlAAjBPRBKqLbYdSDHGJAJvA79s7KCq5dmdfYHxue9xPrQP9vixVsdRNfAbeD0FgZ0YdWYRmWcKrY6jmpE3LfcRQKYx5rAxphxYBMyuuoAx5nNjTHHF5EZAP5+3AZ9/+gGDbEcIHvdNvdpSS2X3Q655iFG2vaxaubzu5ZXP8Ka4dwWOVZnOrphXm/uBj2p6QkQWiMgWEdmSk6MnV7Rmp/NLGfTVKxTbIwhOnmd1HHUFYdfcT4ktjLj9L5Nf6rA6jmom3hT3mppkNX71LiJfB1KA52t63hiz0BiTYoxJiYmJ8T6lanGWr1rFZNs2yofPh4BQq+OoKwkMp2DIXUxlMx+v3WB1GtVMvCnu2UD3KtPdgMsGixaRqcD3gVnGGD1rwoeVlLuITn+JMgmi/cSHrY6jvNBx6rdxix37phc94wApn+dNcU8D+opIvIgEAHOBSwaLFpFhwEt4CvuZxo+pWpIVX25hhllPbv+5ekGO1iK8M2fib+I652es2Lzb6jSqGdRZ3I0xTuBhYAWwF3jTGJMhIs+JyMVzzZ8HwoC3RGSHiOiVAnyU0+XGsf732MTQafrjVsdR9dBl5hMESzm5n/9RT2pqA/y8WcgYsxxYXm3eM1UeT23kXKqFWrkpnVmOFZyKv4mukXo6Q2siHQdwoss0Zh9fwuc7HmdK8gCrI6kmpGeoKq+53YaS1b/CT1zE3vhM3S9QLU6nWc8SJqWc/fQ3GKOtd1+mxV15bc2WHVxXtoLsnjdj6xBvdRzVAPbOg8iOncb1xR+wMeOg1XFUE9LirrzidhuKVj2PXdx0m62t9tas86z/JUTKOPXR89p692Fa3JVX1m7axLTSjzna81bsUXFWx1FXISB2EFmxM5lR+D6b0rXnjK/S4q7q5HS5kVU/wiX+9LzlOavjqEbQ7ZafYRc3+R/9SFvvPkqLu6rTms8+ZILzS7ITHsAeEWt1HNUIAmLiyep1J1NLV7L+izVWx1FNQIu7uqIyh5OYL3/CeVskfWZ/1+o4qhH1vvVZiiSEwM9/hFPPWvU5WtzVFa199yUSzT7Op34HCQy3Oo5qRPbQKI4nPcII1zbWfvi61XFUI9PirmqVczaHpD3P81VAf/pM/6bVcVQT6H/j42T79aDftufIy8u3Oo5qRFrcVa32vvE00VwgYPYLYLNbHUc1AfELxDXjebpxhh2L/tfqOKoRaXFXNTqQvoHR595hR8eb6DpojNVxVBPqmTKD9PZTGXXiNY4cSLc6jmokWtzVZZzlZciSh8mXMPrcoVdMbAu6z32BMvGn6O1vYdwuq+OoRqDFXV1mx+Ln6OvK5PDIH9MusqPVcVQziOrcg32JT5FQvovt7/zK6jiqEWhxV5c4dXAbSZl/ZnPIBIbPuNvqOKoZpd70CNsDUxiY8WvOHdtndRx1lbS4q0rGUUrpmw+QTyjd7vwjohe9blPEZqP9nD/jNDYuvH4/xqXXW23NtLirSvv+9ThxjkzSh/2YLl17WB1HWSC+Vz82DPw+vUt3c2DxD6yOo66CV8VdRGaIyH4RyRSRp2p4fryIbBMRp4jc1vgxVVM7tfldBh55nRXhNzNpth6OacumfO1hVgVdS98DL5Gz61Or46gGqrO4i4gdeBGYCSQA80QkodpiR4F7gDcaO6BqeuU5XxHy0aPsJZ6h9/5OD8e0cXab0O+ev5BlYvF77wFceSesjqQawJuW+wgg0xhz2BhTDiwCZlddwBiTZYzZCegAFa1NWSG5f78Vt9vNmel/oVNUhNWJVAvQvXM0Bye8SICrmNN/vR0cpVZHUvXkTXHvChyrMp1dMU+1dm43J1+9l+jiwyzt91MmXDPK6kSqBZk+eTJv9/ghXQp3c/L1B0GHBm5VvCnuNX1Gb9BPWUQWiMgWEdmSk5PTkFWoRnR+6Q+IPfEJr4XPZ84cPc6uLjfnrod4LegOYrPeI/eTX1gdR9WDN8U9G+heZbob0KCDcMaYhcaYFGNMSkxMTENWoRpJ4Zo/ErX9Rd6xTWP6/OcI8NOOU+pyQf52xs//JcsYR+SGn1O86VWrIykvefMXnQb0FZF4EQkA5gJLmjaWakrlO94k5PMf8IlJpe89fyG2fYjVkVQL1jM6nI53/pV17kQCP/pvHBnLrI6kvFBncTfGOIGHgRXAXuBNY0yGiDwnIrMARCRVRLKB24GXRCSjKUOrhnOmv4Xf+//FZvcAuOVlEnt0sDqSagVG9I0l94aX2eWOQ966G/fe5VZHUnUQq66fmJKSYrZs2WLJttsq98634N0FpLn7kTn1H9w5vnqPVqWu7OVPt5Oy7j4G249in/NPZMD1Vkdqc0RkqzEmpa7l9EBrG+He+BK8+wCb3f3ZOeFvWthVg8y/dhhrRy5kt6snZtE3MFv1GHxLpcXd17nduD59FtvHT7LKlczma17igSmDrU6lWrFHrkvh4+SXWOcahCx9FPdnP9Vuki2QFndfVpqHc9Ed2L94gTeckzk06c88OjPJ6lSqlRMRvjs7hW1j/sJi50Rsa3+J693/AkeJ1dFUFVrcfdWZfThfmoQcWMGPnN/Aff0LPDi5v9WplI8QER6bMYj8a3/Drx23Yd+1GMdLk+HcIaujqQpa3H2NMbDtn7gWTiIv9yz3uX/IuK//kK9fE2d1MuWDHpjQm4S5P+EB19MUnT2G6y/jYfe7VsdSaHH3LQWnMG/MgSWPkFYex/ygX/O9hx5g8oBOVidTPmzmkFge++ZD3Bv4G9LLY+HtezFv3QOFeha6lfysDqAagdsF21/D/emzOEuL+JnjLs4m3M0rtyQREexvdTrVBiR0acdfH5nNd9/sRv9D/+CxjPewH1qD7bpfwJDbQUcabXbacm/tjm7C/HUyLP0220o6M9v1fwy86Qn+cMdwLeyqWUWHBfLyvdcQe8P3me38ORmlHeDdBzB/mw7ZW62O1+boSUyt1andsPrnsG8ZZ20d+FHpPC70upEf3zSEuOhQq9OpNi7zTAFPvbWdXieW8HTg20S6c2HwbTDhuxDTz+p4rZq3JzFpcW9tjm+F9b+FvUsotoXyUtkM3gm6mSduTGZWUhe90IZqMdxuwzvbsvn98u3MLX+bB/w/xt+UI4NugnHfgc56vkVDaHH3Jc5y2LcMNv4ZsjdTagthYfl0/m27gTnjE5k/rhdhgfr1iWqZ8ood/PHzgyzbsJO75EPu8/+UQHcJ9BwLqffBgBvBL8DqmK2GFvfWzhg4sQ3SF2N2v40Un+O0Xxf+UjKFZfbJzBoxgG9O7E10WKDVSZXyypn8Uv685hDLNu3hVrOS+4JX09F5ChPaERl6Bwy5DToN1i9f66DFvTVyuz0Fff9y2LsUzh7AKQF8RgpvlI1hX0gqd4/twx0jehARol+WqtYpp6CMNzYd5fWNX5FQnMb9QZ8xxr0dGy7o0BcG3wIDZ0GnQVroa6DFvbXIPwFZ6+GrtZiDnyCFp3FjZ6c9gUWlI1lhRjFiYDy3De/OxP4x+Nu1g5PyDeVONx9nnOK9bdnsPniYa2Uzc4I3k+jcjWAw4bFInynQZyrEjYPQaKsjtwha3FsiZznk7IOT6XBsE+bIF8j5wwAU20JZb5L4sGwo6xnGgF49uHZgJ25M6kIHPfSifFxOQRnLdp5gafoJso99xQRbOtMCdjFGdhHiLgTAdOiL9BgJ3UdB95HQoTfY7BYnb35a3K1kDBSchLMH4ewBOLULczIdzuxBXOUAFNnCSHMPYJ1jABvdCeSG9SW1dwyTB3RkYr+OethFtVnni8pZeyCHz/adYcPB03Qv2csI2z5G+x8k2XaAMHcBAG6/YKRTAtJpsOdYfefBEN0PQjr49OEcLe5NrTQP8rIvvV04ivvsQcy5TOyOospFCySMDBPPDmdPMtxxZBCPLaoXST06MLJXFCPjo+gRFaLdGJWqxhjDV2eL2JKVS1rWebZmnUPOZzLMlslAOcoQv6MMlKOEm4LK17gCwpGoXtg69ILIeIjqBRHdoF0XCI+FoHYWvqOr16jFXURmAL8D7MDLxpj/q/Z8IPBPYDhwDphjjMm60jpbVHE3xjNcaVk+lFyA4rNQdBaKcqD4HO6iszjyz+AqOANFZ/EvOom/s/CSVTixc5ooDrk6c8h04bCJ5bCJJS8kjpDoHvTuGEZClwgSYtsxMDackADtuqhUQxSUOth3qoC9J/PZezKf/SfzKT6XTefSQ8TJKXrKaeLkFL3sZ+jKGfxwXfJ6p18ozpBOmPDO2CK6EhAeg4RGeVr8wVEQUu2xX8s6LOptca+zwoiIHXgRuBbIBtJEZIkxZk+Vxe4Hco0xfURkLvALYE7DotfBUeJpNTtKPDdnxb2jFOMowlVWgrO8GFeZ5+YuL8ZdXoIpL/YU77J8pCwfKS/Ez1GAv6MQf1chduOqdZP5JpRzph3naMd5056TphcnTAdOmGiKgmMx4V0IjIylc/tQOkcE0bV9MKkxYcRFh7aZ/ucOl4PcslxyS3M5X3qe3NJc8svzKXOVUe4qp9xdTpmrDIfLgV3s+Nv9CbAF4G/3J8geRPug9kQGRhIZFElMcAzRwdFt8pOM27jJLc3ldPFpcopzyC/Pp6C8gEJHISXOEtzGjTEGg8Ft3Pjb/An2C/bc/D334f7htA9qT1RgFO2D2hPmH+ZT+zI8yJ/UuChS46IumZ9f6uDI2WKyzhWRca6IT/NKyblQhPvCMeyFxwkqOUNnOU9nZy4dy3LpfOEsnbMPEEkBoVJW6/YctkAc9jCc/qG4/UNxB4RDQBgmIAwCw5HAMGxB4dgCw/ALCMYeGIx/YCi2gGDwDwa/YPAPuvQ+KMLzuAl5U3lGAJnGmMMAIrIImA1ULe6zgWcrHr8N/FFExDTBMZ+d7/wfift+W+NzgucNVX9TJSaAEgIoMCEUEkwBIRSYUAqIpsAEUyQhlNlDcfiF4/APpzQgCkdgJBLWkcB20USEhhAVGkD7EH+iQgNIDQkgMjSA6LAAAv187wsdYwxFjiJyy3LJK8urLNa5pbmVBTy3NJfzZf+ZX+gorHO9AbYAAuwBuIwLh9uB0+2sddlAeyBdw7rSNawr8RHx9I/qT7/IfvSK6EWAvfWf8FLiLOHQhUOeW94hDl84zOG8w5wqOoXD7ajxNX7ih4hgExuCICKUu8pxXaFhAuBn86v8h7PyPiiSqKAoooKiKudfnG4X2A6btL5eWe2C/BnSLYIh3SKqPTMK8PTOOVNQyrnCcnKLyzlWXM7OIge5xeUUFBZQXnAOd9E5/Epz8S+/QLAzjyBnHsHlBYRSQmhpKWGUECb5hHKGMEoIlRLCKCVYyuuVdffQZxh80/800juvmTfFvStwrMp0NjCytmWMMU4RyQM6AGdrW+n+/fuZOHFivcICZJ8/yvnyi60QgYpf8uqPq95L5X0JIqXY5ELFH4lnvg3B87/nP6CypXOl6crlpcr8ym1zyXPVX1d9nbUtW+M6q6zHYCpbcrXdu40bl3HhMi7PY7frknku93+KrdM4cbqd1Pbvsojgb/PHz+aHn83vP4/F77L5dpsdm9iwYau15WiMwWVcldt2uB2eTwGuXE65TvGl60tKnaW4jbty+8F+wYT6hxLmH0aofyjBfsFe/e5YqcRZQpGjiCJHUWUr/OI+tomNIL8gguxBBNoDKz/VBNgDsNvs2MVeWdirCyGkch+6jbvy5+p0Oyt/piXuEnJN7n9+xlVuNRER/OTSn+/FHHbx/EyrT3v+jmyX/L1V/UeoPi7ul0t+lzF4/r/0d7v6vLpeV9vyF7db/fHF17grlndX3Iypsg0TgMGvch785/7Sx1TeR7/3F7r8dmm99kt9eVPca/rJVP/L92YZRGQBsAAgMLBhx7FCQsIoltJLflCejVWdNhgufnyt8tyVfphXXJ9vsonNUzgqCnKQX1Blob44r3oRt0vjflKpWkhqYzCUOcsodhZXFsnc0lxyij3jhfvZ/C4p9qH+oVdcX1Nzup0UOgorC3mRowiX23VJ1vah7QnxCyHYL5hAv8DKf7Qb4uI+rC9jTOU/qFX/Maj+uNhZfEmDwNdVb4DV1nC77LGArbKhdunrKh5UPg4NCWvy9+HNb0Q20L3KdDfgRC3LZIuIHxABnK++ImPMQmAheL5QXb16dQMiN6+LraKqraOL906385LpS+7d7staU1daptZ7d83zLz6uLM41FOSLRTnYL9hTSPyDCfULJcQ/hBC/EOytuI+wMYas/CzSc9JJz0lnZ85OMi9k4jRO8sgjrl0cg6MHMyR6CIkxifSP7I+/vfG7l5a5ysjMzfRkOLuTnTk7OVbg+aAbKqEMixxWmWFIzBDi2sW1ykMeF7ncLoqdxZWfQoocRRQ7iz3/ILgcnk9eVW8V8+rDLvb/fGKpaFBU/aRgt3k+yVSfV+vyFdMXGzM2seFn87tkuuoyLZ23n4Tq7C1TUawPAFOA40AacIcxJqPKMt8ChhhjHqz4QvUWY8zXrrTeFtVbRvmEIkcRu8/uZtfZXezM2cmus7s4W+I5Muhv86dnu570ad+HXu17Ed8unk6hnYgJjqFjSMcrHsd3uB3kFOdwsugkJ4tOcrzgOAcvHORg7kGO5B+pbM3GBMeQFJNEYkwiiTGJJHRIaBWHjFTr0thdIa8DfounK+TfjTE/FZHngC3GmCUiEgS8BgzD02Kfe/EL2NpocVdNzRjDqaJT7Dy7k4xzGZVfYB4vPH7ZsqH+oQTaAwmwBxBoD/QcjnAUU+wspsx1eU+KrmFd6RvZl77t+9I/qj9JMUl0CunkU71SVMukJzEpVYtiRzHZhdnkFOdwpvgMp4tPe7ptOssoc3luF4+NXzyc1TG4I7FhscSGxtI5tLO2yJVlGq2fu1K+JsQ/hH6R/egXqVcEUr6r5X97oJRSqt60uCullA/S4q6UUj5Ii7tSSvkgLe5KKeWDtLgrpZQP0uKulFI+SIu7Ukr5IMvOUBWRHOBIA18ezRWGE7aQ5qofzVV/LTWb5qqfq8nV0xgTU9dClhX3qyEiW7w5/ba5aa760Vz111Kzaa76aY5celhGKaV8kBZ3pZTyQa21uC+0OkAtNFf9aK76a6nZNFf9NHmuVnnMXSml1JW11pa7UkqpK2gVxV1EnheRfSKyU0TeE5H2tSw3Q0T2i0imiDzVDLluF5EMEXGLSK3ffItIlojsEpEdItLkVyipR67m3l9RIvKpiBysuI+sZTlXxb7aISJLmjDPFd+/iASKyOKK5zeJSFxTZalnrntEJKfKPprfTLn+LiJnRGR3Lc+LiPy+IvdOEUluIbkmikhelf31TDNk6i4in4vI3oq/xW/XsEzT7i9jTIu/AdMAv4rHvwB+UcMyduAQ0AsIANKBhCbONRDoD6wGUq6wXBYQ3Yz7q85cFu2vXwJPVTx+qqafY8Vzhc2wj+p8/8BDwF8qHs8FFreQXPcAf2yu36cq2x0PJAO7a3n+OuAjQIBRwKYWkmsisKyZ91UskFzxOBzPdair/xybdH+1ipa7MeYTY4yzYnIj0K2GxUYAmcaYw8aYcmARMLuJc+01xuxvym00hJe5mn1/Vaz/1YrHrwI3NfH2rsSb918179vAFGn6i6Ra8XPxijFmLZ5rJNdmNvBP47ERaC8isS0gV7Mzxpw0xmyreFwA7AW6VlusSfdXqyju1dyH51+76roCx6pMZ3P5zrSKAT4Rka0issDqMBWs2F+djDEnwfPLD3SsZbkgEdkiIhtFpKn+AfDm/VcuU9G4yAM6NFGe+uQCuLXio/zbItK9iTN5qyX/DV4jIuki8pGIDGrODVcczhsGbKr2VJPurxZzDVURWQl0ruGp7xtjPqhY5vuAE3i9plXUMO+quwJ5k8sLY4wxJ0SkI/CpiOyraG1YmavZ91c9VtOjYn/1Aj4TkV3GmENXm60ab95/k+yjOnizzaXAv40xZSLyIJ5PF5ObOJc3rNhf3tiG55T9QhG5Dngf6NscGxaRMOAd4L+NMfnVn67hJY22v1pMcTfGTL3S8yJyN3ADMMVUHLCqJhuo2oLpBpxo6lxeruNExf0ZEXkPz0fvqyrujZCr2feXiJwWkVhjzMmKj59nalnHxf11WERW42n1NHZx9+b9X1wmW0T8gAia/uN/nbmMMeeqTP4Vz/dQLUGT/E5drapF1RizXET+JCLRxpgmHXNGRPzxFPbXjTHv1rBIk+6vVnFYRkRmAN8FZhljimtZLA3oKyLxIhKA5wuwJutp4S0RCRWR8IuP8Xw5XOO3+s3Miv21BLi74vHdwGWfMEQkUkQCKx5HA2OAPU2QxZv3XzXvbcBntTQsmjVXteOys/Acz20JlgB3VfQCGQXkXTwMZyUR6XzxuxIRGYGn7p278quuepsC/A3Ya4z5TS2LNe3+as5vkBt6AzLxHJvaUXG72IOhC7C82rfPB/C08r7fDLluxvOvbxlwGlhRPReeXg/pFbeMlpLLov3VAVgFHKy4j6qYnwK8XPF4NLCrYn/tAu5vwjyXvX/gOTyNCIAg4K2K37/NQK+m3kde5vp5xe9SOvA5MKCZcv0bOAk4Kn6/7gceBB6seF6AFyty7+IKPciaOdfDVfbXRmB0M2Qai+cQy84qdeu65txfeoaqUkr5oFZxWEYppVT9aHFXSikfpMVdKaV8kBZ3pZTyQVrclVLKB2lxV0opH6TFXSmlfJAWd6WU8kH/D9cSRlN8GC+yAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6, 4))\n",
    "plotPSFProfile(psf2, ax, label='Kolmogorov')\n",
    "plotPSFProfile(matchedPSF, ax, label='VonKarman')\n",
    "plotPSFDiff(psf2, matchedPSF, ax, label='difference')\n",
    "ax.legend()\n",
    "ax.axhline(0, c='k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:28.092817Z",
     "start_time": "2019-01-22T22:52:28.089620Z"
    }
   },
   "outputs": [],
   "source": [
    "def imgDiff(psf1, psf2, ax, **kwargs):\n",
    "    img1 = psf1.drawImage(nx=320, ny=320, scale=0.02)\n",
    "    img2 = psf2.drawImage(nx=320, ny=320, scale=0.02)    \n",
    "    im = ax.imshow(img1.array - img2.array, **kwargs)\n",
    "    plt.colorbar(im, ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:52:28.316318Z",
     "start_time": "2019-01-22T22:52:28.094441Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWIAAAD8CAYAAABNR679AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvX+wHdV17/lZ9557dQVXIIEwwgIe+KFUBb83wfYdoCozLo/tgHAlkVPBiXgug/NwEXugkkxlaoyfJ7GDcZWZScwzFYc8OSY2rpgfYydjJbEfD2yoJK8MQtj4B8YEGTOxQICEJCKBfl1pzR+99z379t3dvfvXOX3O3Z+qrtNn9/7Vffp8zzprr71bVJVIJBKJDI+JYXcgEolEljtRiCORSGTIRCGORCKRIROFOBKJRIZMFOJIJBIZMlGII5FIZMhEIY5EIkNBRDaKyFMiskNEbvQcXyEi95jjj4jIec6xj5j0p0Tk8qI6ReR8U8fTps5pk/5WEfmOiMyLyJWp9q8x+Z8WkWuc9LeIyA9MG7eJiNS9FlGII5HIwBGRSeCzwBXAhcBVInJhKtu1wD5VvQC4FbjFlL0Q2Ay8EdgI/JmITBbUeQtwq6puAPaZugH+BXg/8OVU/04DPgZcAlwMfExE1pjDtwPXARvMtrHWxSAKcSQSGQ4XAztU9RlVPQrcDWxK5dkEfNHsfwV4h7E+NwF3q+oRVf0psMPU563TlHm7qQNT57sBVPVZVf0+cCLV9uXA/aq6V1X3AfcDG0XkLOAUVf22JrPh7rR11aFXt4ImWLt2rZ533nnD7kYkMtY89thje1T1jKrlLxDR1wLz7oIngMNO0hZV3eK8Xw/8zHm/k8T6xJdHVedF5BXgdJP+cKrserPvq/N0YL+qznvyZ+Hr33qz7cxouzKdEOLzzjuPbdu2D7sbkchYMzkp/1+d8q8Bvx2Y9+NwWFXncrL4/Krp9Ray8mSl+/7h5+XPo2zbtYiuiUgkEoSQCEbIFsBO4Bzn/dnA81l5RKQHnArszSmblb4HWG3qyGortH87zX5ev0sThTgCwPx8e1tkPBCSv9AhWwCPAhtMNMM0yeDb1lSerYCNVrgS+Jbxy24FNpuoivNJBsy2ZdVpyjxo6sDU+bWC/t0HXCYia8wg3WXAfaq6CzggIpca3/PVAXUV0gnXRGTwDFIgfW314p03kjRluRmf7w0kgjcJ3KGqT4jITcB2Vd0KfB74kojsILGEN5uyT4jIvcCPgHngelU9DuCr0zT5YeBuEbkZ+K6pGxH5H4G/AdYAvyIif6Sqb1TVvSLyCRJxB7hJVfea/Q8BXwBWAt8wWy2kC8tgzs3NafQRN8+oWqNRpNthclIeK/Db5nK2iP5uYN7/A2q1tdyIt/wYMarCmyZ9HlGYu4H1EUeaJ97iI8y4CG8RUZi7QxTidoi39AiyXAQ4C3v+UZAHS7SI2yPeyiPEchfgNFGQB8/ksDswphT+wInIjIhsE5HvicgTIvJHJj1rEY3MhToi5YhhYGHE6zQYhESIQ7ZIOUL+aRwB3q6qvwBcRDLf+lKyF9HwLtQRCScKSj3i9WuPBid0RBwKr5kmHDRvp8ymZCyiQfZCHZEcokXXPPGaNkvDM+siDkHXzCwx9zjwEskqRD8hexGNRQt1AHahjkgGUSjaJ17jZohC3A5Bwxxm1spFIrKaZBbKz/uymdegRTFE5DqSNT0599xzgzo7TkRhGDzuNY8DfOWxU5wjzVPqx0tV9wMPAZeSvYhG1kId6bq2qOqcqs6dcUbllflGkijCwyd+BtWIFnE7hERNnGEsYURkJfBO4EmyF9HIWqhjWRP9ld0jfiblkBJbpBwh/zTOAr5oHkMyAdyrqn8nIj/Cs4gGGQt1LGfiF737zM9Hd0UIMTStHQpvPfMYkTd50p8heTRJOv0w8J5GejfiRAEeLeIEkXzizLr2iLdcC0QBHm2iIGcThbgd4q3WIOMgwE2fwyiLWRTkxcSoifaI13WZ0/aPR1w5bbyIFnE7xK9FA3TdEu5S//L60lWRjpZxQvQRt8cyv7Xq0SWBS9PlvmXR9QkXUZBjaFpbxB+4MWNc4mLH5TzGjSZXXxORjSLylFmp8UbP8cyVHEXkIyb9KRG5vKjOsqtFish7ReRxZzshIheZYw+ZNuyx14VfQT9RiCvQJZEY9yUgu3h+XerLIGnyKc5mXsJngSuAC4GrROTCVDbvSo4m32bgjcBG4M/Mejh5dZZaLVJV/0pVL1LVi4D3Ac+q6uNO395rj6vqSwGnnEsU4pJ05Qu4XMWgS+fdlX4MioZXX7sY2KGqz6jqUeBukpUbXbJWctwE3K2qR1T1p8AOU5+3TlOmzmqRVwF3hZ1WNaIQl2DYX7yuWYbDpCvXYtjtD5oGhXhhlUaDu4LjkjyplRyzymaln0691SJ/k6VC/JfGLfEHTSzzG4U4kGF+4bogOF1m2NdnOX02JYR4rYhsd7brUlWFrNKYlaep9MJ+iMglwGuq+kPn+HtV9d8D/7PZ3uepoxTLePw3jGF9ycal3UFGGAwz6mI5RFSUDF/bo6pzOccXVmk0uCs4pvPsTK3kmFfWl74Hs1qksXp9q0Wm27BsJmUNq+pz5vWAiHyZxCVyZ865FhIt4hzGRQzz2ml7MGwQbWS1OwzG3TpucPW1R4ENJpphmkTwtqbyZK3kuBXYbCIezgc2ANuy6jRlSq8WKSITJOvm3L1w/iI9EVlr9qeAXwZca7kSY/z7XY9Bf6EGPcNtmPj60oYlOSwLeVxXchOS56Q1garOi8gNwH0kEW93qOoTInITsF1Vt5KxkqPJdy/wI2AeuN48vAJfnabJD1N+tci3AjvNAmeWFcB9RoQngQeAz9W9HtKFpYLn5uZ027btw+7GAuMiwl0S3zK0JWKDFseuifHkpDxW4C7I5edFNPT/98VQq63lRsduleEzSPFquq1RFd40ba1PMWg/7rhZxnGKc3uM0W1Sn1EV4XER4CyaFtBBCuS4iXEU4nYYo1ukHqM2iDSKQl5XkJr0+Q7SOh4XMY4WcXuMwe1Rn1ES4S4KeZU2uyKkgxLJcRHjKMTtMAa3Rj0GIUh126hTvqtui6x+lRWrJsR9UNbxqItxXBi+PZb1dR1nEe6qABdRRxTrCt0ghHLUxThaxO0wwrdEPcZVhJs8r2PHqpWbaiDYtKogj4LQjUIffUQfcXuM4O3QfQbtSqhSpqrINlV3qFhXcT3UcVd0fXH6YROFuB2W5a02qn/b03RNgMtg+1HGeh6n9RxG1SqOQtwOI3gr1KNNER6UP7dM/jrCe/x4tXKToY9oYGn/QoS5jCA3YR23JZijJsZxsK49ltV17aIlXKZPoXmriG9V0Q2tK1Sc3b4XiXJZoeyi8HWxT3nEZ9a1Q+E/DRE5R0QeFJEnReQJEfldk/5xEXnOeW7Tu5wy3udJDZO2LeEqVm1omZC8x471tzyOH/dvbVOl3dBzavpaNlGmbP2jgNDsM+sifUJ+i+eB31fV74jIKuAxEbnfHLtVVf/YzZx6ntTrgQdE5Ofs6kjjRhdWTQuxgKuIbdVzK2Ph2X4VWczHjoVZyG1al6NmvbZB9BG3Q+Ftpaq7gF1m/4CIPMnSR5q4LDxPCvipWV7uYuDbDfS3NF2zZELK1BXfUNEd1qpvPjELcWeEuC1CfcJdC48bhYHIGL7WHqWuq3nU9JuAR0zSDSLyfRG5Q0TWmLSQZ1GNPKMqwlUXaPct8N5EXXnknUvIv4CmrncTZcaFBp9ZF3EIvmYiMgt8Ffg9Vf1X4Hbg3wIXkVjMf2KzeoovWfRYRK6zz7PavXt36Y6H0KV1fovKFAlTnr80z+daJJhFAhsqtFXrKTqed24hPuSQ/ndJjLss8jZqImSLlCNIiM1q9F8F/kpV/xpAVV9U1eOqeoJkhfqLTfaQZ1GhqltUdU5V584444w659B56n65soQmb8CrSIDKWrJNWcR1+pV3vnXjo7ssgF3BuiaiRdw8hT9e5lHRnweeVNVPO+lnGf8xwK/Rf27TVuDLIvJpksE6+zypgdEVa6XuX+M8AS5bV9axpq5VaD3WB5rOn/aN5vl6swb48iaJhPhgq4TDlckfSpf9xVFk2yHko/5FksdF/0BEHjdp/wm4SkQuInE7PAv8NuQ/T2o5MUgRblOAywp1kXhkiUye+GQdO37cH22RF2ERMtgWoyOyiULcDiFRE/+E3+/79ZwynwQ+WaNflWnDGh6kJVxXgMsKbx3/aJ5oZlFkEedZyr52q1rHTYpxm5Zxl34Qmo6aEJGNwGdIQo//QlU/lTq+guQx9W8BXgZ+U1WfNcc+AlwLHAd+R1Xvy6vTPO35buA04DvA+1T1aFYbJjDhSeAp052HVfWDpq63AF8AVpLo4O9qzYd/xh+4hqkzKNeGCBcNkuX1scgfXMVvXLY/oWllfcdtDeKNO035iEVkEvgscAVwIck/7AtT2a4F9qnqBcCtwC2mrDtXYSPwZyIyWVDnLSTzHjYA+0zdmW0YfqKqF5ntg0767cB1JG7XDaYPtRgrIR72F6dO+3VE2CcqVQSv7MBbESHREiHpWWlp2hjIG+V7qmkajpq4GNihqs+o6lESa3VTKs8m4Itm/yvAO8yY1cJcBVX9KbDD1Oet05R5u6kDU+e7C9rwXwORs4BTVPXbxgq+06mrMh3649M9mhYlH3UFuKiNOmlt0+uFuSjSaVnvobyrokm3QpcH2ZoiR6MWo7pWRLY7KVtUdYvz3jff4JJULQt5VHVeRF4BTjfpD6fK2rkKvjpPB/ar6rwnf1YbAOeLyHeBfwX+T1X9R5N/Z0bblRmbW6ZpIWnDLxxKEyLclACXtSaLZrwV+Xyz0sq8h+yBvCqU9dU27dvtjK9YJLwjx47tUdW5vNo8aWk/a1aerHTfP/y8/Hlt7ALOVdWXjU/4/xWRNwb2uzRd+HhrM+y/b1WiFrImZhSVT/toi9oqyhOySFCbpK3iLIs4631Z67iqZTxsMRx2+wuEC3FRjpD5BjbPThHpAacCewvK+tL3AKtFpGesYje/tw3jdjgCoKqPichPgJ8z+c8u6HdpxspH3BRlhL2KCPtIC16oXzTrWNbAmiVkll4VES4qn263qJ82LYsQf3JTPyZN3Rcji7WIQ7ZiHgU2iMj5IjJNMvi2NZVnK3CN2b8S+JYRyK3AZhFZYaIh7FwFb52mzIOmDkydX8trQ0TOMIN/iMgbTBvPmLkTB0TkUuNLvtqpqzJd+I3tFG2LcIglXMbtEGohp9utMilkUBRZwCHWcNqCTLsqBmEZd8aKbYqJCZiZCct74EDuYeOPvQG4jyTU7A4zB+EmYLuqbiWZSPYls3DYXhJhzZ2r4KvTNPlh4G4RuRn4rqmbrDaAtwI3icg8SYjcB1V1rzn2Ifrha98wWy2kZvhbI8zNzem2bduLM3oYpm+4CSEuI8JVBdlts6kJIVnkCY/vmCuOaVF082ftlzmWbs/XZl5fQ47Vydt2fZOT8liB3zaXuelp3b52bVBe2bWrVlvLjZH+vR6WCFeJPy2yhJuwgu1+HbFP06aF3OuFuQ1CLOLQyIoylnGeVezWmcfYDdyNlYnfHeJVNTQlOKH1lPFbFglvXptlxL6orhCywtHcY6Ei6SuXFl3ffhGhERVNid7QxbMpykRNREoRr2pJylqPdSzhMgKc5X6oG2XRFaam/CKfZSkPyjJeVroUhbg1RvaqNikYTbkk0oSGqPnqD/X/uu3kuSCK6qjibnEp8ze+yLXgiiUkgulbyCfEOs4TyyZijUPFuEnRHtoPQBTi1hjJqzoMEW6ijjqWcJ4vOMQCLuODrkKVQb2QMnl+ZJ8VHCLMNj/4LeO2Zt+NvBiLhEdNREoxkkI8CoSGi0G48Lr7Pmu7qGzI+6L+5ZEnDEWugjwXQ7peK5Y+sQ0R5jQhYhwhWsQtEq9qAE34hbPKhIqwzx9sxT1UvLP6m5dehrw6fDPmfHl8gpzGdVUUCW/Vwb0iV0g63T3HsSUKcWuM3My6LrklmnBJpNPzRNidldaUCM/PD2aALqudKoOT9tzd6xH6I5bV9qCefN2l+7c0zc6sizgs2yvW1k2cN3mijiWcJcC+vGVcE4Mm3X6Ir9i1bNN509ZxkWXsc4mU8Rc3wdAG2+oSLeLWGCmLeBgiUqbNKuve5lmJRSLsWpmjIMI+8vroO7c869hXJq+tEMp8pl2/PxshWsStEK9YDmV9w6F1hYhvmfpHVYQtaQuxyNcbWl9e2Sx/cUhIW15fRtbaDaHMWhORUozrLZNLiCjVHaALWVQnS3xDLGHfq69vdQS4qqBUadPnMrDvy9bnc1MUDd65hIa01RHjkRTs6JpojXhVG6Ds2g5uWpEIl6HI6g4h73tWNGnDl6/qP4cqU5ctWT7jUAsZYkiblyjErTEyV7Wpv9aDHvUuM4pvKWsJ1xHg9PfKfR/6nfOJWhah/atjEWfNxEvXXybWOK+OqjRlFQ/Muo5C3BrxqnqoI9ah6zy4aT5BLXJt1B2MyhLgECEOiQnO6lvZsnnWa4hroYo17KunLCPpeghhLE9q+IzEVe3KQFORbzgkFrWMS6INSzgtuFmv7n6WGK1Ysfh92dC60H67Imr7lU5zsf1u2kVR1lc8SAbSjzhY1xoduIUGR91BujR5ftwqIum6JLLKVLWEfQLrRhqlX634uMJT5COemkrOYXKyeEq3zy2QRxlL1tYV4qLw1ZUnamV8xWM3aBddE61ReFVF5BzgTmAdcILksdifEZHTgHuA84Bngd9Q1X3mOU6fAd4FvAa8X1W/U7WDXbaGXYpcEllpaUvYV0+IRZxFWoDTIuwK79RUvkin90POa34+OTdX5OxmxThUkH2DgnkWsbVm0z+YWRM8ssS4yEXRFTFtvR9RiFsjZELHPPD7qvrzwKXA9SJyIXAj8E1V3QB807wHuILkQXsbgOuA2xvvdQWGIeghfuHQcnXJirufnEz+ba5cmbyfmcneTj45cUe428kn55fp9ZK6Z2aStrL60TQh13dUfuQ7RdYHWOEDFZGNIvKUiOwQkRs9x1eIyD3m+CMicp5z7CMm/SkRubyoTvNA0UdE5GlT53ReGyLySyLymIj8wLy+3anrIdPG42Z7XenrmKLwipmnlu4y+wdE5ElgPbAJeJvJ9kXgIZIH9G0C7jRPTn1YRFaLyFmmnk5Txh0QshB7aHvHjuVHSJT1CfssWOvac1993x/r982yjt060/1JW8G2HdcKttvhw8nxw4eT+ux73/XxpedZw1nH0yu25dWd3q/qK+6KtdwIDVrE5gnJnwV+ieQR9Y+KyFZV/ZGT7Vpgn6peICKbgVuA3zSG4GbgjcDrgQdE5OdMmaw6bwFuVdW7ReTPTd23Z7UB7AF+RVWfF5F/R/JA0vVO396rqtUetOmh1BRn82vxJuAR4EwrrubV/iqsB37mFNvJ4hMYecrG+Oa5FEKXs/QdLyLPArbvXet1xYrEerUW7MwMzM7CSTMnmO6dYJqjTHOUiflks++neyc4aeYEs7P9craeFSuWWsi23SwLucw1LXPtQgZCQ6kS5z3yNLvoz8XADlV9RlWPAneTGHEum0iMPICvAO8wrs9NwN2qekRVfwrsMPV56zRl3m7qwNT57rw2VPW7qvq8SX8CmBGR1PB0cwT/vInILPBV4PdU9V+Tc/Nn9aQteVS0iFxH4rrg3HPP9VY0yL9sTbQVEi3g++K71nAdv7Bv4M21dK2V2+slAmsF0X2d4IRjts77O5XR8ESvx7TZmEkamp2dWLCA3deDB5NidoAvr/oQyzjr1XYvJKQtaxCxrhE4SKu41bbKRU2sFRHXYtyiqluc9z6D7ZJUHQt5VHVeRF4BTjfpD6fKWmPPV+fpwH5Vnffkz2pjj1PPrwPfVdUjTtpfishxEk282XgAKhP0kYnIlGnwr1T1r03yi9blICJnAS+Z9J3AOU7xs4HnSWE+lC0Ac3NztU6iiDoiW2aQrk57If7kMvWmDRQrwq5V6lrH070Ti5US+vuBQrzI3DXvJ2ZmjDBPeF0W1iURMhvRR56IpvOUwVemzUG7kXFhhHdyj6rO5RwPMdiy8mSl+/7h5+Uv7IeIvJHEXXGZc/y9qvqciKwi0cX3kQQ0VKbwqhqz/vPAk6r6aefQVuAa4FPm9WtO+g0icjfJr9Ero+AfDqVM3HCeJVcUJVFWdN1XK67Qt4RdC9jdJuaPJo3tP7hYiH1buqNpx7FvM52Znp1lutdj5oxpDh9mYXMtV3tNbFrab1xGnN16XD+v77qVochXPNY0GzURYrDZPDtFpAecCuwtKOtL3wOsFpGesYrd/FltICJnA38DXK2qP7GVqupz5vWAiHyZxCXSrhADv0ii+D8QkcdN2n8iEeB7ReRa4F+A95hjXycJXdtBEr72W3U6WJc2reHQvGUs2nSestaw66Lr9frRED6f8MT8URapYlqI61jE6Y4AzMwwMQMzM9NLzs8WsTHUtnjoNfNZwyFWZtMW9Vhbxc0K8aPABhE5H3iOZPDtP6TyWGPv28CVwLdUVUVkK/BlEfk0yWDdBmAbiXW7pE5T5kFTx90sNRx9bawG/h74iKr+9/4lkB6wWlX3GE/BLwMP1L0YhVdVVf8Jv/kO8A5PfgWur9mvkfMP12krywKu6oqw793BMDs4tkSEDx7sO22zLGJfWIfbMCwORE4rqetwhUVibGOMbbVpP27Za5ElpG5X8/K2xVi01aAQG3/sDSTRCJPAHar6hIjcBGxX1a0k/8S/JCI7SKzUzabsEyJyL/AjYB64XlWPJ11cWqdp8sPA3SJyM/BdUzdZbQA3ABcAfyAif2DSLgNeBe4zIjxJIsKfq3s9uvz7O1R8X/yskfIyg3Qu6Rl0Wa6JUEPUHZizlrB1SdiIhmmOwsGUALuCfOhQv2NFIpzuiCvG9pfAduTw4b5juNdjYnaWk2Zm6K2ZXnBBuOd/6JC/qaxuePQ+63dgiZ+3zqBdl6c9N07DT3FW1a+T/IN20/7Q2T9M/592uuwngU+G1GnSnyFxIaTTvW2o6s3AzRldf0tGemXG8XZZYJCWbhlcbUtTp8/pgbm0O2K6dyIRYdcd4VrCx47BkSOLBdj9Vcj6JbJK5IYnuIHIrqUMi5zY07P9gTzrM7ZFivS/iCxBtPV20b/baRGPM+taY9lf1apf8rxoiSwr2B2kS6fnDez5cK1hq2u+gbkFS3j+tSQcbf/+vgC7r1lCnO6EzzVx5MjiTrlCbDvlWsTuKzDd6zE7e9JCcVeQ3dA26zXJ8h1niViWa8IdHMyqy1dn1VXZOi2yIUQhbo14VQdAnsjm+YaLfiRcnzAstYStT3i6dyIR4SxL2B2ss66J0I5kHXPVyv3VSqubEejpmRMwM8GhQ0vH+NI6X9SdvEE3tx5fuUgOUYhbI17VQOpOaXbz5bklQrTPF6RgB+bS4WkLPmFrCbuv7iBd2jdchXTnXVdF2qSFJfOup2dmWLVqesGAtkVWrOi7mW31eREV6eNZ7gnfwGCeILvHlmUYWxTi1ohXtWVCXA5VdC/9fXCtYVeIl/iE09av65Io6kyeQuWl29epqf4AnrvYhPM6M9sPbbMuCtdfbLtR9pqFiG20iguIQtwanbyqwx5kS7cfsq5ASJREE4N06RA110ecnrQxMX90sS/44MG+JWxf07HC6cbsq/X7+nAH9tIn4/srYEXfirDt+Pw8EzMzC6Ft9rxct7KtMsRN4XahzKBdXvREFmmreNii3kr7DUdNRPp0UogHRV3BLyofcrxqH9Iha2m/8EKcsOsHtg0eOtQfkMsTYTcEzaalR6mOH+9/433+Zfdkof+rduhQUq/bP+N/sHHGtqm0z7iKRex2oY5A1RW4YQt0LaJF3Bpje1UHaVVXjZYoqieNL27Y5xteNFnD3dJREj4Rtg7nqalkseGpqaVLp7nmqN3cyItXX+1bye48ZduWFWsr4q46mdeJGVi1anoh2OLwYZYIc0gERag169brKz9IAe2sWEchbo14VVMMQsDrWMJp0pESvR6LXQ4+qzjLEnb9HL0erFqVqN+qVY7TmaV+AVv3gQOJuLoCm27HtYzTI3BuX3s9er3pReeV9hVXpQnLOLSdsdKtKMStEa9qScoKQagPMyRaAvzWsOsrnu6dSBbwcS1id5Du0KF+rLCLtX5Xr+6L7tq1yatNM77cE2aRq4UlM20bNiLDmrBWZK2V7Lso9ottozccH+T02hmYmfCesy0aGsrmNpWXPySfpc5TnkeWKMSt0LmrOuyBulHCndC2IFCuFWq3Y8f6VrDvl8Sam+5q8LOzfQFeu5YTvWn271884Dg5OcHU1DSrV5+WuENg8YLDsNiKTn+4doDP9i9tRc/P0+tNLxLfqan+WF8kn8Yt8mgRt8ZYXtVQK3TQ+NaWKBstkd7cJ2tMcGJpmJq1gI8c6Q/SuY1a8bWiu25dsj87y4l1r+fgQXjh2X7AhbtcsW0/yT7NunWvZ3YdTMzM9AvA4hC19MISdlEJN5LCftkPH2ZitsfKlRMLXXUPp326RaQtXvt+0FZtiEB20q1RbmH4SAm69lGPDUOx7H3WsDtJI/1LAItH+6wVbF5374ZXXoEXXkh0dc+evgsX+u7kw4eTIgCnngpnrl7dr99axr6IC9d3kN5f5F+eZpB0UgS7QLSIWyNe1RI0MVDUFK5/eGF9HdfR7HNNpHELp1wSR3sn8fLLiVG7a1fi4n3hBb8Qz88nLma7v2bNSUxbLZ6d7StbevEfS5ZrwpxLrze9ZFG36MIaElGIWyFe1QJClr7Moimx8IWtpd/3evTXk3AjEKwlnLaGbaGTT06iItauhdWrObHu9ezeDS+/DD/+cSLEO3bAvn3w4ot9Twf0XSJnnglr1iRibY3h008/iTPWncSE+5jmAweWuibckDY3FM6W6fXoeSLn0lF0VWOLq+BazMtqqnO0iFsjXtWG8YlBekZdE+Fr6S987kBdVmM2XthZMf7gwcQdsX9/sr3wQrLZfet2hv7TniERZ+ue2L8/6c/KlXCK9UG7bWUF/rphdYsG7Jaee51/J2lfcdEMu4ghCnFrxKvaIULF2Z3YYLdFT192Rdi1NNMWsVXqjQnSAAAgAElEQVTSVasSU3btWl54NhHcXbsSS/iFFxLLeN8+eO65xOWrmvxNEJlidjY5tmZNUu3Bg4mRbZs65by1yc6qVUlfVq70T4e2qudxUUxwgl5vYpEVHDrAFsW0QeJgXWvEW7QmofG/IXVUZYnQuH/3s7BK5kzJO9GbXgg9fvXVRGD3709eX34ZDhw4BhwCEiFWneLAgZW8/HJiSu7fn7grXn21H8J8ojedRFG4E0LyVDRjinRdMa0jyK7rY9mL+rK/AO3ge/x0pIOk7/8llqFnkGvhNe3odgfOnIE6647YsyexhF98MXl97jnYvfsY8CLwsnnt7+/efYznnltcZs+evntjYXX6JaOLDq4LxTfoyNIpzlnXJtIS1jURsgVVJxtF5CkR2SEiN3qOrxCRe8zxR0TkPOfYR0z6UyJyeVGdInK+qeNpU+d0023UIQrxCFN4v6eDltOkpqkdO7Z0bSA7Qa9vCR8CDputn2Yn17mhynbpiUUjbD5ruKifZc450h4NCrGITAKfBa4ALgSuEpELU9muBfap6gXArcAtpuyFJA/5fCOwEfgzEZksqPMW4FZV3QDsM3U33UZlohC3QJuj9757PNM1EVphL5m2bF23dt6FnQOS+IRdET6U2j+G6rGFiAp3WYvjx0mmRJewlHznEHTeDRLD4zw0axFfDOxQ1WdU9SjJY+43pfJsAr5o9r8CvENExKTfrapHVPWnwA5Tn7dOU+btpg5Mne9uso2g65dDFOJIJBJGOSFeKyLbne26VG3rgZ8573eaNG8eVZ0HXgFOzymblX46sN/UkW6rqTZqEf/otUCbMa3z88mAWDpt0dyzspaniUyYnJxYcBm7U6dFplCdAlamCs+YtClEpli5Minjrpo5OZmK6AgldQ5ZEW9tEV0gS1GFo/PBttseVZ3LOS6+JgLzZKX7OpeXv8k2ahFvtxHGJ8qLcNXE9+uQCiObmppeNKZmo9tmZ+HAgURwF2PTVi6Mx61cuXgMcGqKxYNvWYsO+fYzzjkyHFQbvf47gXOc92cDz2fk2SkiPeBUYG9BWV/6HmC1iPSM1evmb6qNWhT+vInIHSLykoj80En7uIg8JyKPm+1dzjHvSGOkHukvwJIINd/UM/vqm6ngTp4wI3KrVy+EE7NuXTJjbt06WL8ezjhjCjiT5F/bmYv2zzhjivXrF5cxE/WSmXZ2JC9vkokbSZExlbDo4R+RdrFCHLIF8CiwwUQzTJMMjG1N5dkKXGP2rwS+papq0jebiIfzgQ3Atqw6TZkHTR2YOr/WZBtBZ5xDiEX8BeBPgTtT6beq6h+7CamRxtcDD4jIz6lqh1ZpaBY3vrSqINR1ZczPw7T7SYbMerCq5qzUNjF/lNnZaWZnk0kZa9Ykg292ssbhw1McPDi1ZELH6acneVavTl5PPrm/bMXC00LclddC4psDXBNlqONqSP8+LFeatIhVdV5EbgDuAyaBO1T1CRG5CdiuqluBzwNfEpEdJFbqZlP2CRG5F/gRMA9cbzXGV6dp8sPA3SJyM/BdUzcNt1GZwltLVf/Bja0rYGGkEfipObmLgW9X7uEyIlSQjx9PDEjXAjnBBBPpARM7pTi9MIN9tfFpBw4sLFm5bt3rF/ry6quLpy2vWWND1BIL27ou7MqZF1yQ7J9zTt+qXggotutM2Pg296TtZqdBp87jBBNLrK3QKc7LXTybpsl/H6r6deDrqbQ/dPYPA+/JKPtJ4JMhdZr0Z0i0KJ3eWBt1qHOb3iAiVwPbgd9X1X0ko4cPO3kaGVEcJXxiaj0DVjzSmliF9ByNxCr2jF7nLcxw/PjiZ8sdPszsumQpy/n5/gI+Bw/2fdF5i/5YQV69OqljdhbYYyxhty0f6Ukmzpa+TiHP/8sjHWGV9YDqKOKLadhHHHGoeqvdDnyCZLTwE8CfAP+RsJFQAEw4y3UA5557bsVutM/UVPYKkkU3ZVPRE2nXx/x8f+0cd5u2I2T21RacnOz/5U/PXnv11UQcTf6JmRnOXL2aNWtOAhJj9uST85fBXLcuyXPOOX3LeHr+NXhh/+Ipdr4Hlrr9c0f5nFef/zFrIuEgSIu4ZaxXXgNOnFj8LNhIc1QSYlV90e6LyOeAvzNvQ0ZCbR1bgC0Ac3NztcM/BsHkZHfWJHbX8rHr5DDjGeSylqbPMnYLu8+4A6ZXJ0tZ2h+Agwf7PwZpIV67NrF+165NLOHp+dcS4XWfmZc1SAf9/vkmBfR6zB9ePLbYlc9guREt4vaoJMQicpaq7jJvfw2wERVbgS+LyKdJBuvsSOOyo81Y4txG05tradp9t2PudDroP94IOGPdSaw0ocMHD/YfVWSLW6PVPtRj3Trjjnhhf/9RSa4QpxU0q3/lZ2k1SnRJZBOFuB0KbzkRuQt4G8lMmZ3Ax4C3ichFJG6HZ4HfhvyRxkEySLdBGdLegTK+Yt9yDFY/7aJms7MT/ZXObAarpNbJm44Bc5XVljt4EGZnmTh8mFNmZjjlvKyHh/Yf/DwxfzRxQ+w5nLxaa9i6JOxTmi1usPKKFX2Hs7tS28wMJ5hYNHW6RrjUEl232j9oQoS+iz8G0SJuj5Coias8yZ/3pNn83pHG4A4Nw5IcUdyle+1r5oCddU34Fma3bgNY/AQNa/KSPBD0tNnZhSgGcGbM7TdWrxVd+2q3LJeEtYCzXBPGP+zOO6k7ULecaFrMoxC3Rwd/d7tNWT9xqHUOYfHIrgVo87rzMpiZYNrGnMHixw7ZCtzKLK++unQixeHDyQmvWrVoPeGJtClvBffAgeTi7N/f71TaErb1py1gOzXPeT06P7HECs5a5z7vuoZ6OMp6QoZhTQ+TOFjXHlGIUwzCIm8ifM3ie7rQQvTE/PziV0i+Se7AXVqYrfhOTSXCavO6c5/dE3B/Baz1a199T412rXTbz1SURFa0RFPW8KBcz110L9QlWsTtMIa3SsIgXRxp90C6bfuFLBKSrPIWXxibfXV1cOaMaSZmWJrJFnafdunmgb7J41rE1s/sc6q6Ynvo0GKT1XeCK1YkInzyyUm9riVsN/O0kAP7+iHOWdZw1mfsWsNF+NapT5cfpKh2VcCja6I9OvqRD4a6Yl1UPuQ4VOuDT4xcT8HMjBFj1+q0YrxyZX9WRtECDlYBs2Y9wFLrN0uEXUs4vTqQYw2f6E0vzAGpMziX1YU6DLv8MIlC3B6dvC2GPWCXbj9rUoevTJ5Vm55hl9dmFmn/MPQjJ9wQMzBibAbYFjplC0D/cczuyWWFZ9hG8zrlw5bJsoRXr17kG7Yi7G7W4k//+JSJlvBRNKOurFUd0uagaKP9KMTt0UkhHieyBNZNr/LD4+oqLI6gcAdUTrK+4fTr4cP9JyrbirKsY7fRMrjxwdYCtq+pMDW7pQXY5xuuIgbpsLWiPJGlRCFuj3jrBeJaxW7kRKiIuvmyIi9CoyfSx21gxNRUf98KWG/NNNOzKUs4HaLmFjp0aKnPpKzpaV9dAXajI3q9viXsWMRH5yc4cCDxmrhWsH0f4htOdyPPGk7nT+9n1ZsuP+5Tmy2qMWqiLaIQD4A8qxj8gQVFYuyWSw8GuukLIW2uGEIiiO5kjnQFWb4YX2fT++4qatYlkWEBW9+wDVVL+4XTXpMyIWu+Y9Eark60iNtj2d9+Vf3RPqu2avRE3qBdSASFNWih3yc7Nmfzzs6exPTMiSTBjS22lRw8mHTMmp/uAJyt6Ngx/yNBrEnouiLSQmwt4Yx4YTsT2oYdF/mHfeSJrO+YL1oiXZfveNX44VEX+yjE7THit0Y+wx70y6KJQTsfrl664r4g1DMTTLvTn93BO/e9Va4F/0bBqJXPEk77htOuCfPeWsLpiRtufHRVyg7SdYEui3UU4vbo8Mc+XLKiHoqWxMyyin31+daecPWwaDAv70thI85seUjcv6tWTTMzO52sSeFaxNaE9i3o4C4wkeWayFjUfSEsLWUR2+gIOxnPWsL2fdmwtbSHxPfbYX8bssrmRUtkCeRyWsc4CnF7dPaWGaQ1O+y2inzFIaRdIOm/8e6AHjhxxvagzWwLu6Kc5ztxTyAtvu6+u5CPE6KWHohLhyOXtYrL+IYHKZjj0lYU4nYIfjb2qFLnpixTNmS0PaS+dJ4yZWGpL9U+DcnVVrud6E0v9temt9Rst0X76ffpvL66PHHCbr/cpyiVEeCsa1TmehcN8oUyqPttGNi1JkK2OojIaSJyv4g8bV7XZOS7xuR5WkSucdLfIiI/MA8xvk1EJK9eSbjN5P++iLw5rw0ROUlE/l5EfiwiT4jIp5z87xeR3c6DlT8Qcs4d/+i7R1YYm4+iaAlLOhSuyABNk85nrd/0AJ47uDczY9fzmabXm2Z6bco1YSsqM0rmWsH2vRMVMT/PkmnLvgG5MqFqWd3xCavrG64qessxbM0yQNfEjcA3VfVTInKjef9hN4OInEayJO8cyXK8j4nIVvPItttJnv7zMMmz5TYC38ip9wqStdM3AJeY8pdktQEcAf5YVR80T3L+pohcoarfMN27R1VvKHPCy0KI67geisrmRU/UbSfEz5xFOp+7yqUrVO5AXq83zcSs45pw3RS+St2O2tfUZh/8mV4vIj1hI9QXnNe0++rLU4Yq0RLjbA3DQIV4E8ka6ABfBB4iJcTA5cD9qroXQETuBzaKyEPAKar6bZN+J/BuEiHOqncTcKeqKvCwiKwWkbNM3iVtqOpdwIMAqnpURL5D8jSiynT64x+277ZOHUWDdm5eK+auJVxm0M9SlH7kSD+azNblunOTIIcJILGSe+bYBCfyGzAdX/S0ZWcWtTv4lhZgu+SF66YOtYSLBtncHxzfBI46g3RlGBf/MJT6jqwVke3O+y3m8WghnGmfAKSqu0TkdZ4864GfOe93mrT1Zj+dnldvXl2+9AVEZDXwK8BnnORfF5G3Av8M/G+q6tbhpdNC3FVC1p5wyRNTX11Fohv6o+HmcX8g7KqWsNhCdr/Eboxtrzdh9qcX6lhUv8eVYM8pPRjnCrJbR1k3RJbg5FnGrishxJLOY7m5JaC0RbxHVeeyDorIA8A6z6GPBtaf9aDi4AcY161LRHrAXcBtqvqMSf5b4C5VPSIiHySxvN9e0H4UYpcywld3yrObN/2lTot2miqWsW99eGsFu6LrunfTIpX11zwd2eaKstueb45IFX9wqEVrzw/CfcNZ1nCRb7joh2EcaHJheFV9Z9YxEXnRPhfTuAhe8mTbSd/NAIlr4CGTfnYq3T7AOKvenfgfepzVhmUL8LSq/mfnvF52jn8OuCXrPF3GPmrCMowvRJ6/Mq8/bfTVZ5Van60bsZA3Ev7qq4mF626vvppfxo3caHpZyzxCrm9XRLIr/SjCWsQhW022AjYK4hrga5489wGXicgaE/1wGXCfcT0cEJFLTbTE1U75rHq3Aleb6IlLgVdMPd42AETkZuBU4PfcThmBt/wq8GTICXf+Fhikn7hOP7KsYl8kRJa1e+zY0nrS+Mrl9Svdj3T7rigdObLU+vX9nc9qw54D+K1k++qzfqtYwnY/S1jdcyjrD86yhvP6NWwG0Y8BfRc/BdwrItcC/wK8B0BE5oAPquoHVHWviHwCeNSUuckOqgEfAr4ArCQZpPtGXr0kkRXvAnYArwG/BZDVhoicTeJC+THwHRMd96eq+hfA74jIrwLzwF7g/SEnLMlA4XCZm5vTbdu2Zx5v8sMv61stSnf9u+noCZ/QFL26QhZaJq/PaaoMVIWurZD19KWQ19B+5/U5xCVRZYCuTbdEk+JZVNfkpDyW57ct4nWvm9Mrr8z+nrrcfnu9tpYbHfktz6fLVnGZuGK3jqxXt74i/7DPyi4ince1mNNtFglMyA9WluCW+Tzz+pZnEYeKsK8tS1URHjSD6Eec4tweHbmNukUd4S9yUfjaSb+m68krU7XfaQEuc75FefMs9bLXtYqfPSRUraiNdD1l6YpAN0kU4vYYmdulKau4bj1ly2f5i/PqsFZxWSFIW8mhZUKPhfiIyxzLo6wrAsrNeivzDyCkj1VoSqwHJfpxYfj2GBkh7jLpWOC0NVvGGva5KEJJiz6UF0I3f5lBwSr5fKTbrCIyTbsk3DqXM9Eibo/C21xE7gB+GXhJVf+dSTsNuAc4D3gW+A1V3WfCRT5DMgL5GvB+Vf1OO12vTohVm5enyFcMxa4Fdz9PjMvgs4irCnLVMlXJE+BQSzhkwkaIBRwiwnk/EIMepBsUUYjbIySO+Aski2a42MUzNgDfNO9h8eIZ15EsnjGy1P2yZeUP8XuWqT9EZLr8xS8jwmXqC7WA3f0Qd1CT98UoMcA44mVHoRCr6j+QxMO5bCKZuod5fbeTfqcmPAysTgU412IYN3mZNqv8fQ0RZVtveqF3X8TAqIlxWUu4yBquMiiXR5nPtOv3Z12iELdH1Y+x7OIZu6p3sR2aGvxLkw5ng+woCsiOmsjyGactNl84Wxa+PIP+0pSxTvNE2Tcwl/djlCXQedZwW37hLv4ghhAH69qj6SnOwQtuiMh1IrJdRLbv3r07uIEmb+JBjXpnPbrdV1eeiExN5VvHvjZCXCKDEIasduqIsHs9Qq5fXtuhESqDumcGXVcI0SJuj6ofZdnFM5ZglsTbAsnMuor9GAh5A14hA3d5ZfKs4HSa7UNReJtvskcWaQs9L08RVQapyrgm8qzgorJ5+z6W68I+ecTBuvaoahGXXTxj2ZH+IudZXFkCkicmeUKRLpsun/c+q38hW0gdIe+zzt2SJ8JF181HDFULI1rE7VH4Wy4id5EsBbdWRHaSPDqk1OIZTdOkf7eJuvIsZhffgj55lrG7H2LhuovL5xHiT/aVySNEiPPS8kQzbynLOvtu3ZY2pzGPslsCokXcJoUfp6pelXHoHZ68Clxft1MhDEOMy7bpc1HkrUdRJMZZfUi7KoqmWPt+OIp+TKp+8YsEOP0+y29bJJBVXBB1pjCH1F8l36DrKksU4nZYJt6t5sgT4yyRhOzJHj7L2L4v2rdUmfwxKMoKcZ74lvH7FvmiQ90ReaK3XHzDliYXho8sZpndStk0ZWGH1hOyUpuvziK3RV47eaLvtuVS9pqUEa4icfTlqzr4FtJOXrt1GBfBjq6J9hjpW6RJ90SZ+kL+xjdpGafrDLWUfYT6katS1oIMsYDz9sscS7fnazOvryHH6uQdRn1liELcHiMtxNC8GA+aogWC8ny7oT7gvEkmbn63jO9YFbLK5wliqNAWiW6oxT0qdMGyHuXvWpfpwEfbLcoIe1l/MRRbxr6yZf3IvvfDIkv88qzfovdFeX1pTVjCIcer5h0FokXcHsvm4aFlaOrLVsdv2euFCU7WsXT59Ht3Vlq6H+5WlqLy6XaL+mnTsig6b9unJljOIgyDiyMWkdNE5H4Redq8rsnId43J87SIXOOkv0VEfiAiO0TkNrMqZGa9Zt7DbSb/90XkzQFtPCQiT4nI42Z7nUlfISL3mLoeEZHzQs55LIR42Dd9kVD4jvuE0CdeeUJl97PeF+Xp9fr9yNrSwlq0FdWXbr9Kv4uuS971zPoByvqc3OPDZNjtQz9qImSrSdbqjguYpXg/BlwCXAx8zBHs20lWf7QrQdrVI0utGlnQBsB7VfUis9nZxdcC+1T1AuBW4JaQEx4LIYbhD4o02b7PgivrD82yLEPSXIqENb1lUac/Zd9Ds/7gYd4LbdRXhwHNrMta3dHlcuB+Vd2rqvuA+4GNZsmFU1T122Zew50sXh2yzKqR3jZK9P0rwDusRZ5Hhz7i7tHkQGCezzj0SRx5A3n2PdRLc6l67lX8rFXTyohUVb9wGbokmk1T0ke8VkTcRz5vMevLhJC1uqNL1kqP681+Oj2v3ry6fOmWvxSR48BXgZuN8C+UUdV5EXkFOB3Yk3fCY3XbNCmcVeosylt2AM/ii3DwLdYTmmbTs0QjZCGgEKr81a8j0qEDgyH9CD1eNe8w66yD6onQrHtUdS7roIg8AKzzHPpoYP1ZKz0GrwBZs673qupzIrKKRIjfR2J9V2l/vIS4C4SIMZSzjn2TP3ztZKX52svrR9GXv45QlxHgrPQmRDik710TweGjQOAspKKaVN+ZdUxEslZ3dNlJsgaO5WzgIZN+dirdrgBZdtXIrDZQ1efM6wER+TKJD/lOp66dItIDTmXpgzWWMDY+YktbVknTllFWnqJohqJ++QbDyqRn5Stqt2x9ZdOz2g2NzkjXU0TZz3s5WMOJEB8N3GqRtbqjy33AZSKyxgygXQbcZ1wPB0TkUuObvZrFq0OWWTXS24aI9ERkLYCITJE80/OHnjauBL5lXBa5dO6jHhdCXBp5ecpax5Bv3ZaxiH35mqKsRZx3rA1XRGie5Uuwa6IO3tUdRWQO+KCqfkBV94rIJ4BHTZmbVNVanh8iedbmSuAbZsusl4xVI7PaEJGTSQR5CpgEHgA+Z/J8HviSiOwgsYQ3h5ywBIh168zNzem2bduLM5akreDzsvWG5M/LkzWYl7dWRVGbTZ9D05EFece7IsBtCXZb9U5OymN5ftsiRN6k8GBg7jW12lpuxN/+AVB3EDHPOga/IBdZvKEWcTp/XdoQYKi/+ly0gkNozkccWcxY335tRFFUrbfOIB4sFhrfGseWLLeFpepaEm25MIrK1RXftlwRo2YNN0cU4jbo/MfeVdoQ49A8ebHHRctr+kLbQhhGrO1yE+HuEy3ithj7W6rsX/CydVddrzfEOs7LF2ohW/LcF1m0NaHDR8hMuKbEt0y+psoNq95mUSBw9lGkFCPx8XeZttwfZeoPmZ2X50/Oa7ttQqciNynCVRkNsWyTaBG3xbK5tdoUzCpWd5kyIXnzLGSXLOErI9BVqLL2Q+gAXNMx3k2U6VL9zRKFuA1G6haoS9vWaxXK9ClUvENF2SXUnVG1rlDKRD90JdysDl3sUzbRIm6LkboNmqBrlrFbLrRsmXbSwhYqzDCYJ1pUCTtr2wJuomwX6m+HgUzoWHaM5K1Qly5axlWoIvy+BYaGQdsC3HVG81yiRdwWI3k7dJ06kRpVwsuqlAkRwqpiXXdyhcugfbqjKZCDwq41EWmaZXvbDcIqrttGW+FxoTQpqGXpWnhZk4xCH/1Ei7gtat0SIvIscIDk05lX1TnzeJF7gPOAZ4HfMKvbd45xFWNbDkbPBTNMa3YQAjm6ImyJPuI2aGIZzP/FPLPJLvBR+LypLjGoL19dgalaR3pZybaWbBxWv5o4p0Fdky5c93pYizhki5ShjfWIQ5431SkG9QVpop2mRGMYotxkm03WMwhGX4QtUYjboO7tocB/ExEF/ot5JlXI86YQketInpjKueeeW7Mb9RlUJEVTLoOqa0YU1dVlhrHexai11S5xsK4t6t4iv6iqzxuxvV9Efhxa0Ij2FkjWI67Zj0YYZFhbk22Nqj84lKaFLIpwVZToI26HWreJqj5vXl8Skb8heW5TyPOmOsugxRiaF2TLqArzuCysM14ibIluhzao7CMWkZPNE0wxjw65jOS5TSHPm+o04/KF7dLgXBFt93VcPtPhEgfr2qLOYN2ZwD+JyPeAbcDfq+p/JXku1C+JyNPAL5n3I8cwvrhtC1FXIicG1Zdhnet4ijAMSohF5DQRuV9EnjavazLyXWPyPC0i1zjpbxGRH4jIDhG5zTxENLNe89DQ20z+74vIm/PaEJFVIvK4s+0Rkf9sjr1fRHY7xz4Qcs6VbxlVfQb4BU/6y8A7qtbbJYY1FXrQA4cuTbc7LFFabu0OjoH4iG0I7KdE5Ebz/sNuBjNf4WPAHMkvxGMistXMWbidJBDgYZIHg24keYBoVr1XABvMdokpf0lBGxc5fXkM+Gune/eo6g1lTriN8LWxYljW4zCtuSa3YfV90Az7X8ZgsFETIVstQkJgLwfuV9W9RhjvBzaacalTVPXb5jH2dzrls+rdBNypCQ8Dq0093jbcTojIBuB1wD/WOeEoxIF04a98xM+wr8/y+WwG5iNeFAJLInRp1gM/c97vNGnrzX46Pa/evLp86S5XkVjAbuTXrxsXx1dE5Jy8E7Usm1uoCYblqnDbt4xqRERTdEX8utKPwREssmtFZLvzfosJWQVARB4A1nnKfTSwfvGkaU56W3VtBt7nvP9b4C5VPSIiHySxvN9e0H4U4rIMW4zdfkA3+jJIuiR8XerLYCgVR7zHWfZgaU2q78w6JiIhIbA7gbc5788GHjLpZ6fSnzf7WfXuBM7xlMlqw/bzF4Ceqj7mnNfLTv7PAbdknadLdE1UYNh/hV26EgnRFl08vy71ZbAMzDUREgJ7H3CZiKwx0Q+XAfcZl8MBEbnUREtc7ZTPqncrcLWJnrgUeMXU423D6cNVwF1up4zAW34VeDLkhJfl7TTOjIulvDyFbhQYSIzwp4B7ReRa4F+A9wCIyBzwQVX9gKruFZFPAI+aMjep6l6z/yHgC8BKkmiJb+TVSxJZ8S5gB/Aa8FsABW0A/IYp5/I7IvKrwDywF3h/yAnLYh/zcJibm9Nt27YXZ+woXRe9rvfP0nXx7Xr/ipiclMfy3AVFiPwbTUWR5XB9rbaWGyN+a3WDrluheQIy6D6PopiNYp/bIa410RbxFlvmtL0+RRSxcSNOX26D+DVpkK5bxiFE4ewTr0Wa+Kiktoi3WguMgyAvZ6IAZxGFuC3iLdciUZBHiyjARcSF4dsi3noDoCuTQCLZRBEOJQ7WtUG8/QZEnJ7cPaL4liW6Jtoi3opDIFrIwyeKcFWiELdBvB2HRLSQB08U37pEi7gt4q3ZAaKF3D5RhJsi+ojbIN6eHSFayM0TxbdpThCjJtoh3qodJIa91SMKcJtE10QbxFu2w7Q9/XhciMI7KKKPuC3iLTxCREt5MVGAh0H0EbdBvJVHkOUuyFGAh0W0iNsi3tIjzHJxXUTh7RJRiNsg3uJjxLgIcxTerhKjJtoi3vJjTJcWhE8TxXZUiRZxGxr+ggoAAAZcSURBVLT28FAR2SgiT4nIDhG5sa12ItVIP5SzbWHs4kNAI2WxT+gI2SJlaEWIRWQS+CxwBXAhcJWIXNhGW5Hm8IlzU1tkXGj/Kc4icpqI3C8iT5vXNRn5rjF5nhaRa5z0t4jID4wReJt5mnNmvebpzbeZ/N8XkTc7df1XEdkvIn+Xavt8EXnE1HWPiEyb9BXm/Q5z/LyQc27LIr4Y2KGqz6jqUeBuYFNLbUUikYFgoybaFWLgRuCbqroB+KZ5vwgROQ34GHAJid58zBHs24HrgA1m21hQ7xVO3utMecv/DbzP08dbgFtNXfuAa036tcA+Vb0AuNXkK6QtIV4P/Mx5v9OkLSAi14nIdhHZvnv37pa6EYlEmkOBY4FbLTYBXzT7XwTe7clzOXC/qu5V1X3A/cBGETkLOEVVv63JI+rvdMpn1bsJuFMTHgZWm3pQ1W8CB9yGjYX9duArGXXZNr4CvMNa5Hm09afR17AueqO6BdgCICK7JyflVWBPS/0pw1qG348u9AFiP9KMej/+Tb1mX7kP/nZtYOYZEdnuvN9ivvMhnKmquwBUdZeIvM6TJ8vYW2/20+l59WbVtSujf6cD+1V1PpV/UV2qOi8ir5j8uZ9XW0K8EzjHeX828HxWZlU9Q0S2q+pcS/0Jpgv96EIfYj9iP9Ko6sbiXGGIyAPAOs+hj4ZW4UnTnPQqdVXJX6X91oT4UWCDiJwPPAdsBv5DS21FIpERQ1XfmXVMRF4UkbOM1XoW8JIn207gbc77s4GHTPrZqXRrBGbVW8pwJLFuV4tIz1jFbn5b104R6QGnAntz6gJa8hGbzt0A3Ac8Cdyrqk+00VYkEhk7tgI2CuIa4GuePPcBl4nIGjNIdxlwn3E9HBCRS41v9mqnfFa9W4GrTfTEpcAr1oXhw/ieHwSuzKjLtnEl8C2TPx9V7cQGXDfsPnSlH13oQ+xH7McQz+90kqiGp83raSZ9DvgLJ99/BHaY7bec9Dngh8BPgD8FpKBeIQm3/QnwA2DOqesfgd3AIRJr93KT/gZgm2n7/wFWmPQZ836HOf6GkHO2HYxEIpHIkGhtZl0kEolEwohCHIlEIkNm6EI8zDUpRORZMxXycRvzGDq9sma7d4jISyLyQyet9PTLlvrxcRF5zlyTx0XkXc6xj5h+PCUilzfYj3NE5EEReVJEnhCR3zXpA70mOf0Y6DURkRkR2SYi3zP9+COTfr40OK020iGG7JSfJHGQvwGYBr4HXDjA9p8F1qbS/i/gRrN/I3BLC+2+FXgz8MOidoF3Ad8gGVC4FHik5X58HPjfPXkvNJ/PCuB887lNNtSPs4A3m/1VwD+b9gZ6TXL6MdBrYs5r1uxPAY+Y87wX2GzS/xz4kNn/X4E/N/ubgXuavmfj1u42bIu4i2tShEyvrIWq/gNLYwtLT79sqR9ZbALuVtUjqvpTklHhixvqxy5V/Y7ZP0AS8rieAV+TnH5k0co1Med10LydMpvS8LTaSHcYthAXrknRMgr8NxF5TESuM2mLpkECvumVbZDV7jCu0Q3mL/8djmtmIP0wf6vfRGIFDu2apPoBA74mIjIpIo+TTDq4n8TaDppWC9hptZERYdhCXGk6YIP8oqq+mWT1petF5K0DbDuUQV+j24F/C1xEMtf+TwbVDxGZBb4K/J6q/mte1jb74unHwK+Jqh5X1YtIZm1dDPx8TlvD/h5FajJsIS47tbBRVPV58/oS8DckN/yL9m9uzvTKNshqd6DXSFVfNCJwAvgc/b/arfZDRKZIxO+vVPWvTfLAr4mvH8O6Jqbt/SRTdy/FTKv1tLXQjzLTaiPdYdhCvLAmhRkB3kwyRbB1RORkEVll90mmSP6QsOmVbdDI9Mu6pHytv0ZyTWw/NpsR+vNJ1m7d1lCbAnweeFJVP+0cGug1yerHoK+JiJwhIqvN/krgnST+6man1Ua6w7BHC0lGwP+ZxAf20QG2+waSEe/vAU/YtsmYBtlw23eR/MU9RmLNXJvVLjnTL1vqx5dMO98n+YKf5eT/qOnHU8AVDfbjfyL5K/194HGzvWvQ1ySnHwO9JsD/AHzXtPdD4A+de7axabVx684WpzhHIpHIkBm2ayISiUSWPVGII5FIZMhEIY5EIpEhE4U4EolEhkwU4kgkEhkyUYgjkUhkyEQhjkQikSHz/wM/l8LMomQHdwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6, 4))\n",
    "imgDiff(psf2, matchedPSF, ax, vmin=-1e-5, vmax=1e-5, cmap='seismic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-22T22:54:04.528921Z",
     "start_time": "2019-01-22T22:54:04.515893Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.0009574]\n",
      "[1.]\n"
     ]
    }
   ],
   "source": [
    "# Quick and dirty PSF photometry\n",
    "result, _, _, _ = np.linalg.lstsq(\n",
    "    np.atleast_2d(getImg(psf2, **psfKwargs).array.ravel()).T, \n",
    "    getImg(matchedPSF, **psfKwargs).array.ravel(),\n",
    "    rcond=-1\n",
    ")\n",
    "print(result)\n",
    "# Check that using true PSF gives flux=1\n",
    "result2, _, _, _ = np.linalg.lstsq(\n",
    "    np.atleast_2d(getImg(psf2, **psfKwargs).array.ravel()).T, \n",
    "    getImg(psf2, **psfKwargs).array.ravel(),\n",
    "    rcond=-1\n",
    ")\n",
    "print(result2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}