notebooks/2_plot_psf.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.lines as mlines\n",
    "from utils_shear_ana import pltutil\n",
    "\n",
    "import astropy.io.fits as pyfits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAJlCAYAAAASFxJOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOU0lEQVR4nO3df3Cb92Hn+Q9ISqREhoZI2ZZNakODzXRWm7g1ocTXbLXx1GCbZHPnRAtI2+buqt6uhGwxd1K1HSHupY7S1JXBbFX7bGxDuu25003mbCJKO9Nt0wVziattcq4FOGcnzt0leqTUon+JIiFalEiKIO4PFTBBAiB+PM8DPMD7NeMRCTzP9/kS4GN8+P3pymQyGQEAADhEW70rAAAAUAnCCwAAcBTCCwAAcBTCCwAAcBTCCwAAcBTCCwAAcJSOelegWayurmpmZkaStH37drlcrjrXCAAAZ8hkMrp+/bokaefOnWprK922QngxyczMjO688856VwMAAEd76623dMcdd5Q8hm4jAADgKLS8mGT79u25r3/84x9r165dll9zfn5ed999t15//XX19vY2XLm1lFPNuZWcU+6xL7zwgh588EF985vf1P3331/Rz9CMrPqdM4PddWvk+6+R771Kjuf+e1cj33tS7fVbWFjI9V6s/TwthvBikrVjXLZv367u7m7Lr5lOpyVJ3d3dpl7PrHJrKaeacys5p9xjt23blvvXjve00Vn1O2cGu+vWyPdfI997lRzP/feuRr73JHPrV86YUbqNAACAoxBeLPDAAw9oz549ikajll6ns7NTn//859XZ2dmQ5dZSTjXnVnKOVa9ds2vk183uujXy/dfI916t9WtVjf6a1VK/aDQqr9db0TkudpU2x8LCgnp6eiRJ165da8hmPVTuu9/9rj784Q/rO9/5jn7u536u3tUBWgr3X+uo9DOUlheghF27dikUCtkyABtAPu4/FEPLi0loeWlOi4uLunjxooaGhtTV1VXv6gAthfuvddDyAgAAmhrhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhBQAAOArhxQJer1d79uxRNBqtd1UAAGho0WhUXq+3onNcmUwmY1F9WsrCwoJ6enokSdeuXVN3d3edawQzLC4u6uLFixoaGlJXV1e9qwO0FO6/1lHpZygtLwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFEILwAAwFE66l2BRpZKpTQxMSG3261EIqHh4WGdOHGi3tUCAKClEV5KOHz4sB5++GGNjIxIknbs2CG3260jR47UuWYAALQuuo1KMAxDU1NTue89Ho8SiUQdawQAAGh5KWF9UEkmk4pEIhWXc/r0ac3Pz6u3t1fHjx83q3oAALSkhg8vqVRKhw8f1sGDB+X3+4seNzY2pitXrqi/v1/nz5/X6OhoyeMrrUM4HNb4+Lh8Pl/F558+fVrT09MaGBggvAAAUKOGDS+BQEB9fX2SpFgspoMHDxY9NhgManh4OK9VZHR0VLOzszWPT0mlUpqampLb7db58+drKssqtOwAAFqJK5PJZOpdiVIMw9Dw8LAmJycLtqQkk0l5vV6t/zGKPV6L4eFheTwexePxDc8tLCyop6dHknTt2jV1d3fnnhscHMy1vFy6dMm0+thVfitbXFzUxYsXNTQ0pK6urnpXB2gp3H+to9RnaCEN2/JSrvHx8dxsoLWyj8VisVzoicViBYPHWl6vV0eOHFEqlVIgEND4+Lg8Ho8kyefzaWJiwuSfoLHRqgMAaDSODy9TU1MFw4skud1uxePxXHjx+/1lj4M5d+6czp07l/fY7OxsLsi0CsbrAAAajePDi2EYRQfR9vX1bQgg5fL5fDpw4EBeWJmamtLTTz+96blvv/12XpNXtusqk8locXGxqvqUYmX5Vte90S0tLWllZUVLS0v1rgrQcrj/mtvS0lLuvV1YWKjoXEeHl1QqVfJ5t9u96TGlhMNhhcPh3Aymp59+uqyWm/WtM9kgk06ndfHixarrU0w6nbasfCvLlqRnnnlG165dU09Pjw4dOmR6+bVaWVnJ/Q51dDj6dgEch/uvuT311FOKRqNVnctvQwkej6eqdV0Mw8hrebnvvvu0sLCg9vZ2DQ0NmVjDW9rb23P/ml2+lWVL0p/92Z/p9ddf1913362TJ0+aXn6tsn8V7N69W52dnXWuDdBauP+a26OPPqpHHnlE0q2Wl0qGZTg6vLjd7pLP19LqUos77rgjL7y4XK7cv1aMmLeyfCfX3azBxh0dHers7GS2A1AH3H/Na+172lLdRptpxQG2eBeDjQGgOTk+vLjdbs3OzhZ8LpVKae/evTbXKF86nc41fS4tLSmdTue6YuBsTCMHgPpw/MaMBw4ckGEYRZ8fHR21sTa3eL1e7dmzR4cPH9bQ0JBmZmYkSTMzMxoaGtKZM2dMu1ahcAR7nD59Wl/4whd0+vTpelcFABwrGo3K6/VWdI7jw0sgEFAymdwwviW7G3Q1exHVKpFI6Hd/93f1x3/8xxtWvJ2enpbf7zclwJw5c8bScEQwAgBYLRQKbdgIeTMNH16yoaRY15DP55Pf79epU6fyHo9EIpqcnNx0UK8V0um0jh49WnBrguxjx44dqykMnDlzRn6/37JwZHUwkghHpZw+fVonT56kVQcACmjYvY3C4bAMw1AymZRhGHK73fL5fOrr69P4+PiG463cVboca/dl+Ku/+it9/OMf3/Scb33rW3rggQcqvlY6ndbQ0FDRfYxcLpcGBwd14cKFqsbXZIPR+l+N7MygWCym/fv3V1zu+mscPXo072cYHBzUE088UXPZ0q3XaNeuXZqZmdHOnTv15ptvVvValNpbxco9pdivCmBvo1ZS6d5GysAU165dy0jKSMr8yZ/8Se7rUv999atfrepa3/rWt8oq/1vf+lbFZa+srGQGBweLlulyuTK7d+/OrKysVFX3TCaT+drXvpZxuVwFy3a5XJmvfe1rVZedLX/9zzA4OFhVuTdu3Mj88Ic/zNy4cWPDcwMDAxlJmYGBgZrqW4iVZQNOUer+Q3NZ+xl67dq1TY9v+G4jJ9q1a1dZx911111Vlf/GG2+YetxaZ8+eLfmXfiaT0WuvvaazZ89WXLZkfZea1d1pWXR5AUD9EF4scOzYsZJLWbtcLu3evVv79u2rqvxyQ0814cjKYCRZG47sGGsk2TMeyGqMqQHQKFpytlEjSiaTevbZZ+VyuXLjRLKy3z/++ONVr/eyb98+DQ4Obih77TWqDUdWBiPJ2a1Gkj0tO3a06jDNG0CjaMrZRk61f/9+xWIxDQwM5D0+ODhY84DX9vZ2PfHEE5JkejiyMhhJzm41smsWmdNbdQDAaoQXC+3fv18XL17Uzp07JUk7d+7UhQsXTJlNY1U4sjIYSc5uNbK6Zceu8ToA4HSEF4u1t7fndkPt7Ow0dWsAq8IRrUaFWdmyY9d4ney1rOyWYjwNAKsRXhzOqnBEq9FGVrbs2DFeR7KnW4rxNACsRnhBUbQa5bOyZcfq8ToS3VIAmgfhxQLZjRmj0Wi9q9LQnNZqZGXLjtXjdezslgKASjBVukEkEgm9+uqrCoVC9a5Ky7IyGFnRsmP1eB27uqVYvA9ApZgqDdjAipYdq8fr2NUtZcc0bwYEAyC8AFWwomXHyvE6VndL2TmehgHBAAgvQAOxaryOld1SjKcBYDfCC9BgrGjVsbJbyq7xNBJjagDcQngBWoRV3VJ2jKeR2DoBwLsILxZgqjQalRXdUlaPp5HsH1PDgGDAPtVMle6wqC4tLZFIqLu7u97VAAoyu1sqO55menq64LgXl8ulwcHBqqd5bzamxuVy6dixY3rooYdM6WI7ffq0pqenNTAwoOPHj9dcHoDSQqGQDh06pJ6enrLPoeUFQE2snuZt55gaAM5Aywvq6vjx45qfn1dvb6+jyrajfCfJjqc5evRoXtAYHBzU448/XlO3lF1jagA4B+EFJdkRAKxidZO/1XV3WjDav3+/HnroIe3atUszMzO58TS1duXYMaYmq9BsJjP39AJgDsJLE7C69QL2c+rrbsU0b6vH1GSdOXNGR48e3TCb6YknnjBlN3UA5iG8NAGnftAB5ciOqfH7/XK5XHkBxowxNdK7s5nWh6PsbKZaVzjOOn36dO4PDe5boHqEFxs4sQsAzc1pv5NWjqmxczYTM5kAcxBebMD/pNBonPg7adWYmkpmMz3wwAM1XQuAOZgqbQEWqQOsYcWYGmYzAfXFInUNgkXq0Mqc1iVl12wmZjIBhVWzSB3hBYCpnNYlZcdsJmYyAeai2whAS7N6hWD2ZQLMR3gB0PKs2nF7s5lMknTs2DGl0+mqyl/v9OnT+sIXvkB4QdOj2wiAo1g1psaK2UzMZAKsQXgB4ChWjqkxezYTM5kAa9BtBAAWqfe+TECzIrwAgEWyM5nWDwTOcrlc2r17tyn7Mg0NDW2YzcRgYDQrwgsAWMTqmUySfbOZGAyMRkJ4AYA1jh8/rs9//vOmja2xaiaTZP9sJqBRMGDXAl6vV21tbQqFQgqFQvWuDoAKWDEgmH2ZgOKi0aiefPLJis4hvFiA7QEArOfkfZnY2gBWqmZ7ALqNAMCh7JjNZMdgYKBShBcAcCirZzPZubUBUAnCCwA4lJWzmRgMjEZGeAEAB7NqNlMlg4EBuzFgFwBs4qR9meze2uD06dO518bKLSDQHAgvAGATJ+3LZOfWBtKt8DI9Pa2BgQHCCzZFtxEAYAO7tjYAqkF4AQBsYMfWBkC1CC8AgIKs3NoAqAVjXgAARVm1tcFarOCLShFeAKBJWDWbyYqtDbLOnDmjo0ePbljB94knntDHP/5x066D5kJ4sQAbMwKoB6fN0smu4Lt+IbzsCr5f/epX9bM/+7M1X4dp2I2tmo0ZXZlCyyeiYgsLC7lNpa5du8bGjE1icXFRFy9e1NDQkLq6uupdHaBuBgcHc1OZSy1eV650Oq2hoaGiZblcLg0MDOgb3/iGhoeHa7r/zK47zFfpZygDdgEAtitnBd9Lly4pkUjYWCs4BeEFAGC7clfmvXz5ssU1gRMRXgAAtit3Zd7bb7/d4prAiRiwCwDYlNkzmbIr+E5PTxfcuTo75sXr9ZpyPTQXwgsAYFNmz9LJruDr9/vlcrnyAkx2Bd8vfelLNU/LZg2Z5kS3EQCgLjZbwfeTn/xkTeWfOXNGQ0NDG9aQOXPmTE3lrnf69GmdPHlSp0+fNrVcFMdUaZMwVbo5MVUasF46nc5bwffNN99Ue3t7TfdfsTVksq06Zm5vwFTs2jFVGgDgKGav4JtOp3X06NGCY2myjx07dkzpdLqm66B+CC8AgKZSzhoyr732ms6ePWtjrWAmwgsAoKmUu4ZMuceh8RBeAABNpdw1ZMo9Do2H8AIAaCrZNWSyg3PXc7lc2r17t/bt22dzzWAWwgsAoKlk15CRtCHAZL9//PHHTVnvpdA6MrAe4QUA0HQ2W0PGjGnSdq0jg41YYRcAUHdmbz8g3QowDz30UN4aMhcuXDClxaXYOjLT09Py+/2mriODjQgvAIC6M3v7gSyz15CRNl9HxuVy6dixY3rooYfYisAihBcAACpQyToyDzzwQE3XOn36dK5FyqqA50SEFwAAKmDnOjKnT5/ObT1AeHkXA3Yt4PV6tWfPHkWj0XpXBQBgMtaRMVc0GpXX663oHFpeLJBIJNiYEQCaVHYdmenp6YLjXlwulwYHB1lHpkyhUEiHDh3KbcxYDlpeAACogJ3ryKAwwgsAABWyYx0ZFEe3EQAAVbByHRmURngBADQ1KxbAy7JiHRlsjvACAGhqTp1iXGjfJMLRLYx5AQCgwbBvUmmEFwAAGkh236T1q/hm900iwBBeAABoGJvtmyRJx44dUzqdtrtqDYXwAgBAg6hk36RWRngBAKBB2LlvkpMx2wgAgAZh975JTt21mvACAECDsHvfJKfuWk23EQAADYJ9k8pDywsAADUwewXf7L5JR48ezRu8Ozg4qMcff5x9k0R4AQCgJlZ0t7BvUmmEFxs883xS1xaX1NPVqUMfGal3dQAADsC+ScURXmzwzPNJvXX1mu68rYfwAgBAjQgvAADAdFZOwya8AAAA01k5DZup0gAAwFEILwAAwFEIL2VKpVIKBoP1rkZBzzyf1FN/810983yy3lUBAMByjHkpUzgclmEY9a5GQcxmAgC0ElpeypBM0qIBAECjILyU4dy5cxodHa13NQAAME06ndbS0pIkaWlpSel0us41Kl/Dh5dUKqVAIKBYLFbyuLGxMYXDYY2NjSkYDG56fLlisZgOHDhgSlkAADSCM2fOaGhoSDMzM5KkmZkZDQ0N6cyZM3WuWXkadsxLIBBQX1+fpFsB4uDBg0WPDQaDGh4eViQSyT02Ojqq2dlZHTlypOo6pFIpSZLb7a66DKdjawMAaC5nzpyR3+9XJpPJe3x6elp+v1+xWKzhN39s2PAyOTkpSTIMQxMTE0WPSyaTmpiY2PAmRCIReb3emsLLc889V9P5zYDBwADQPNLptI4ePbrhM1OSMpmMXC6Xjh07poceeqih91Jq2PBSrvHxcY2MbPxQzT4Wi8Xk9/tzX8fj8ZLlZQPP1NSUfD6f+RUGAKBOzp49q0uXLhV9PpPJ6LXXXtPZs2f1wAMP2FexCjk+vExNTRUML9Kt7p54PJ4LL36/P/d1OdaOm4nH4zIMQ2NjY/L7/fJ4PLVVHACATRw/fjy3P5AZ3njjDVOPqxfHhxfDMIq2kPT19encuXNVlevz+TaUG4/HdeLEiU3Pffvtt9Xd3Z37fnV1VdKtRLu4uFhVfUrJNv9ZUb6VZTvB0tKSVlZWciPyAdiH+0/69V//9dzXZvw/uL+/v+zjar3eZp8fS0tLufd2YWGhorIdHV6yA2qLcbvdmx5TjrGxMT377LO5lpcjR46UHMS7vlXmZ37jS5KrQ+n0ii5evFhzfdZLp1dy/5pdvpVlO8HKykrud6ijw9G3C+A43H/mGxgY0K5du/TWW28VHPficrl05513amBgoOb/52enXqfT6YJlPfXUU4pGo1WVzW9DGU6cOFFWi0uWYRh5LS/+//jnuvzOdbW3d2hoaMj0+rW3d0hasqR8K8uWpD/7r/+3FpaW1d25Vf/Dz/+M6eXXKvtXwe7du9XZ2Vnn2gCthfvPGn/wB3+gX/mVX5HL5coLMC6XK/f88PBwzdfJDvhtb28v+Pnx6KOP6pFHHpF0q+WlkuEYjg4vm01hNqPVpRp33HFHXnhpa7u1nI7L5VJXV5fp18v+wllRvpVlS9JXvvNKbjbTYd/9ppdvho6ODnV2dlry8wMojfvPfP/6X/9rbd26VUePHs0bvDs4OKjHH3/clGnS6XRay8vLkqTl5WVt2bJlw+ylte9ppd1GDb9IXS1mZ2dbeo0WAAAK2b9/vy5evKidO3dKknbu3KkLFy6YElzsWADP8eHF7XZrdna24HOpVEp79+61uUYAADS+9vb2XHdcZ2enKeu6ZBfAWz8dO7sAnlkBxvHh5cCBAyV3e67HnkRer1d79uypeiASAABOs9kCeJJ07NixDXsoRaNReb3eiq7l6DEv0q1tBCYmJpRKpfK6iKampiSpLgvNJRKJvDEvqA+2NgAA+1S7AF4oFNKhQ4fU09NT9rUavuUlO+i2WNeQz+eT3+/XqVOn8h6PRCKanJxkzEsLe+b5pKL/5QU983yy3lUBgKZn5wJ4DdvyEg6HZRiGkslk7vt4PK6+vj6Nj4/nHTs5OZnbVbq/v1/nz59XMBisaDVdAABQvbvuusvU40pp2PCydofoclSyDgsAADDXvn37NDg4qOnp6aIL4A0ODmrfvn01X6vhu40AAEDja29v1xNPPCHp3TXCsrLfP/7446bMaiK8WIDZRgCAVrR//37FYjENDAzkPT44OKhYLFZwHZmWnG3UiJht1BqYzQQAG+3fv18PPfSQdu3apZmZmdwCeMVaXKqZbUR4Aar0zPPJ3NYGhBcAeJcVC+CtRbcRAABwFMILAABwFMILAABwFMILAABwFMKLBZgqjVo983xST/3Nd9naAEDTY6p0g2CqNGrFTCYAraIpN2YEAABYi/ACAAAchfACAAAchTEvQAtiawMATkZ4scF9r13R8sINbZ1fqndVAEkMCAbgbHQbWWD9VOkP/eSyHnxrXh/6yeU61wwAgMbCVOkGsX6q9Jb0at6/AADgFqZKAwCApkfLCwAALer48eOan59Xb29vvatSEcILAFMxkwlwjuPHj9e7ClUhvDSBT7zyD9LSstS5td5VAZjJBMByhJcmMDR7TW2SVheW610VAAAsx4BdAADgKIQXAADgKIQXC6xfpA6AeZ55Pqmn/ua7eub5ZL2rAsAELFLXINYvUgfAPAwIBppLNYvUEV5QEjOZAACNhvCCkpjJBABoNIQXAIDlMpmMbt68qdXV8vd4W1pa0urqqpaWliysGczQ3t6ujo4OuVwuW65HeAGANVgh2FzXr1/X1atX9c477yidTld0biaT0crKiqanp237UET1Ojs75Xa7tWPHDsvfL8ILAKzBgGDzvPPOO7p06ZK2bNkit9ut7u5utbW1lf3Btrq6quXlZW3dulVtbUyObVTZkHn16lW99dZbWl5e1q5duyy9JuEFdXXfa1e0vHBDW+dpFgaayfXr13Xp0iX19vbq7rvvruov8dXVVblcLnV2dhJeHOA973mP5ubm9Oabb2rbtm2WXovwgrr60E8uqzO9qqX2a/WuCgATXb16VVu2bKk6uMCZduzYofn5ec3Pz1t6HaIs6mpLejXvXwDOl8lk9M4776i3t5fg0oJ6enp0/fp1dXRY1z5CeLEAK+wCKKRVVge+efOm0uk0i3W2qK6uLq2uruq2224r63hW2G0QrLALoJBWGQycnQ7NOJXWlH3fu7q6yjqeFXYBAA2jEbuMTp8+rfn5efX29ur48eP1rk5Tyr7vhw4d0sc+9jH19vaafg3CC5oWWxsAWO/06dOanp7WwMAA4cViv/Zrv1Z260ulCC9oWmxtAADNifACAE2C1YHRKggvANAkWmVAMMBQcAAA4CiEFwAA4Ch0GwFVYjYTWgnjadBICC9AlZjNhFbCeBo0ErqNAAAtIZ1Oa2np1g72S0tLSqfTda4RqkV4AQA0vTNnzmhoaEgzMzOSpJmZGQ0NDenMmTN1rhmqQXgBADS1M2fOyO/369KlS3mPT09Py+/3E2AciDEvFvB6vWpra1MoFFIoFKp3deBA9712RcsLN7R1fqneVQFsMfPOglZXM2prc2nne8zb2DadTuvo0aPKZDIbnstkMnK5XDp27Jgeeughtbe3m3ZdO6RSKQUCARmGIcMwNDc3J8Mw9Oyzz+ael6RwOCyPx1PzeVaJRqN68sknKzqHlhcLJBIJvfrqqwQXVO1DP7msB9+a14d+crneVQFsceWd63p7fkFX3rluarlnz57d0OKyViaT0WuvvaazZ8+ael07uN1uxeNx+f1+SdJzzz0nSYpEIopEIhofH9fo6Ki8Xq+mpqZqPs8qoVBIiUSionMIL0AD2pJezfsXQHXeeOMNU49rRKOjo7mvR0byZ4L5/X4dOHBAgUAg16JS63mNgPACAGhad911l6nHNTKfz1fw8XA4rFQqpXA4bOp59UR4AQA0rX379mnXXXfL5XIVfN7lcmn37t3at2+fzTWzj8fjkcfjqbgLqNrz7EB4AQA0rfb2dv3W7/yuJG0IMNnvH3/8cccN1q2Ux+ORYRi2nWc1wgvQgj7xyj/o3/zozVtbHABN7hc//gmd/vIf6c5d+V1Dg4ODisVi2r9/f51qhmoxVRpoQWxtgFbj+9i/1C9+7OP6yMgHNDMzo507d+rChQumtLhYNc3bTIZhbBiUa+V5VqPlBQDQEtrb29XZ2SlJ6uzsNK2ryKpp3mbJrudSbGCu2efZgfACAEATSCaTBR+PRCJyu92KRCKmnldPhBcAABrUzDsLevvqNc28s7DpsS+++OKGwbUTExN67rnn9M1vftP08+qJMS8ATMXWBoB5rrxzXTfTq9rS3rbpeJpgMCjDMBSLxW6de+WKJOnChQtyu92mn1dPhBcApvrQTy6rM72qpfZr9a4K0HJ8Pl/FY1Rm3lnQz97/cw05tqUYwgsAU7G1AeAsV965odveuW76TKmv//0PdOXGsnq6OnXoI+bOWCK8AAAA0339xVf1yutXdOdtPaaHFwbsAgAAR6HlBQDQMo4fP675+Xn19vbWuyqmye76XOnuz9nj35m/am6FbEB4AeAon3jlH6SlZalza72rAgc6fvx4vatgmlQqpUAgoHPnzkmSAoGAPB6P4vF4Ref9+393WP/kve/V3/3t85bX2SyEFwCOwtYGwC1ut3vToLLZef/v65dzU7GdxFm1BQAALY/wYgGv16s9e/YoGo3WuyoAADS0aDQqr9db0Tl0G1kgkUiou7sxdxYFAKCRhEIhHTp0SD09PWWfQ3gBgDXY3gBofIQXAFiD7Q2AxseYFwBYg+0NgMZHeAEAAI5CeAEAAI7CmBcAQMt45vmkri0uWbLTMexDeAEAm7C1Qf0983xSb129ZslOx7AP4QUAbMLWBoA5GPMCAAAchfACAAAcxZbwcvHiRTsuAwAAWoAtY17C4bCeffZZOy4FAC2LAcFoFbaEl0wmY8dlAKClMSAYrcKWbiOXy2XHZQAAQAuwJbzQ8gIAqLf06qqWV1YkScsrK0qvsn+VU9HyAgBoei/8+DUdefovNLewKEmaW1iU73f/RP/l5R/XuWaohuljXj772c/qwoULue8zmYxeeuklHTx4UJlMRi6XK/dvJBLR0NCQ2VUAAJjsvteuaHnhhrbOL9W7KhX77o/+Qf/hL//rhsffunpNx/70L/X4r35Cv3jvT9WhZqiW6eHlscce2/DYgQMHmG0EAA72oZ9cVmd6VUvt1+pdlYqkV1f1R/9nouBzGUkuSaf+4tt68P0etbc5a+mzVCqlQCAgwzBkGIbm5uZkGEbu8zaVSkm6NePX4/EULWf+6lX973/4pNpWluV2uyVJo6Oj8vl8isVicrvd8vl8Vv84FaHbCACwqS3p1bx/nSJhTOvKtetFn89IejN1TQlj2r5KmcTtdisej8vv90uSnnvuOUlSJBJRJBLR+Pi4RkdH5fV6NTU1VbCMH7zysj76zz8k944dGh8fz52bSqVkGIbi8bhtP08lnBUzAQCowOX5BVOPa0Sjo6O5r0dG8jeb9Pv9OnDggAKBQK4lZq1fO/iv9N/8/D7921//nzecF4vFcoGo0RBeAABN6/beblOPa2TFunbC4bBSqZTC4fCGx+evXtVvPPzbBc8bGRkpGHgagS2L1O3YscOOy5huampK8XhcH/zgByVJL774ooLBYMm+QwBA5W67sSytZqQ2c4cZeD0D6u/ZXrTryCXpTnePvJ4BU6/bSDwejzwez4auo6mpKfXedpt2v/e9Bc/bu3evHdWrii0tL1/+8pftuIzpDMPQxMSEAoGATp06pdHRUYILAFhg68qqOlcz2rpi7pia9rY2/dtf8BZ8LhuTHn7oAccN1q2Ux+ORYRh5jyWTSb3/3p+pU41qU/G7dfDgQX3ve9+zoCqNaW5uTplMRolEouFGWwMANvdz7/sn+s1P/Lz6e7bnPX6nu4dp0g5VcXi55557io5anp+fL6sMdpkGANjp/p/arYnDD2lHd5ckaUd3l6b+1/+pZYKLYRgbBvN6PB699g8/qVONalNxeCk2eOdrX/uaDh8+vOn5H/zgBzU8PKz3ve99ZYWY7Dz2WCxW8rixsTGFw2GNjY0pGAxueny5JiYmNDExobGxMU1MTJhSJgDAPrfdWNbOxZvqW1rR1o5bQz23dnQ0fVdRVnYdmPW9Bz6fT6/95Ceav3q14Hmzs7N2VK8qZQ/Y/d73vqdMJqPJyUklEhsX/DEMIzewtZinn35aiURCJ06c0MGDBxUOh4suXhcIBNTX1ydJisViOnjwYNFyg8GghoeHFYlEco+Njo5qdnZWR44cKefHK2jv3r3yeDy5RXu8Xq/6+vpyc+oBAI1v68qqXJIyq829z14ymSw4LjMSicjtdud9RmYf/z+efVZ/FP3fFP7tz284b3x83LK61qrs2JnJZDQ7O6tAIKBvfvObBY+5cuVKyTLGx8e1Y8cOPfbYY7rvvvtKjmSenJzU+Pj4hqld6yWTSU1MTOjEiRN5j0ciEQWDwZLnbmZkZCQXXKRbYaaR30wAQOt68cUXNwzKnZiY0HPPPVfwc9vtduvxL/+RYl/9T/rO3z6f91wsFlN/f7+l9a1F2S0v9913n6R3Wx/6+vr0qU99Kvf8yMiIPvOZz+jUqVMFz3/ppZeUTCY3DSOVGh8f39CPl62PdOsNyLaUxGKxTVcL9Hq9udaaiYmJvJab4eHhhl2wBwBQH9uXV7S6mlFbnVcfDgaDMgwjN2wi26Bw4cKFvD/E1/rwv/iIvvF3f68nHvtd/T8vvajh4WGlUin5fD6NjIyY/pltlorXeXG73Xrsscfk9/s1Ojoqv9+vVCqlU6dO5bp5ChkfH5fL5cq1hrz00ksFQ0elpqamipazfulkv99fdpePYRgKBoPy+Xy5ZrgrV64wVRoAkGf78oraMtJqA+yE4/P5Kp4Z23vbbfqdsd/XT999e97j61txGklVi9SdOHFCV65c0Ze+9CXF43FlMhn5fD6Nj4/r4MGDG8axXL16VRMTE/J6vbldpLMtObUqNAgpq6+vT+fOnauqXI/HoxMnTuSFlampKT388MObnvv222+ru7vwao2Li4tV1adcVpbfinVfWlrSysqKlpZK76TbiHVvhPKdXHery2/mui8tLSmTyWh1dVWrq5W3RmTPWf+vOTK5f80q15V5919z67q2tsXLXvs6VXr9UuVXW+7q6qoymYwymVulZzKZgr8zS0tLuf+3LixUtj1D1SvsRiIRfeYzn9HU1JQ8Ho8efPBBSbcG2v7SL/2SJiYm9N5/XLUv2+xUzgd/JTZbttjtdte0tHEwGFQ4HFZ/f7/Onz+vhx9+uKyWm/WtM5O/8Bu5r62eJm5l+a1Y95WVldzvUEdH8dulEeveCOU7ue5Wl9/MdV9dXdXKyoqWl5er2ph37R8LKysrFZ9fyj9+niqT0aZ/lFTD9DLXVLhY2Tdv3pQkLS8vV379EuUvLy/nyq+k3OXlZa2srCi9mpYkpdMrBX9nnnrqKUWj0crq+49q2h7gnnvu2TA92u/358KMy+XS7Oys5ubmFAwGtX///louZzuPx7NhdHY5DMPIa3n5219+LPd1tuXJTN9f87XZ5VtZttXlm1F29obdvXu3Ojs7TS+/mFZ/3Zu1/Fap+9LSkqanp7V169YN900xN9Z8nT0n+9f+li1b1FbDtOa1ZWezlMulsutWSflmlZnlXpzNbZtQrOwtW7ZIUkWvd47rmqSM5NpY/tatW3PlV1JuJpNRR0eH2tvaJUnt7R0Ff2ceffRRPfLII5JutbxUMizDkr2NRkZG9OMf/1gvvfSSzp07p71795rWTbRWsQFIWfXaUOqOO+4o2m3U1dVl6bWtLL9V697R0aHOzs6SZTRq3etdvpPrbnX5zV53l8ultra2qkLH+nOqLaeQQx/x6triknq6Ok0rcy2zy1w7zbtY2dkFYufn5yu+vmvNv+vPrbbctrY2uVyuXKuby+Uq+Duz9jHbuo3Kcd9991kSWso1OzvLAFsAQM6hj9Q+UaRRZBdxzY7tDAQC8ng8m86qXavYTKnR0dHcgN3Dhw9rfHxck5OTmzYa2MWWXaWt5Ha7i64CmEqlGnpXTAAAqpWdUVuLYjOlai3Xao5fG/nAgQMlp3ONjo7aWJtbvF6v9uzZU/VAJAAA7LB2ppTZHvj/3tC/+dGb+sQr/1DyuGg0Kq+38M7fxTi+5SUQCGhiYkKpVCqvOSu7eWQ9doJOJBJFx7wAANAKBq5eV3phWasLyyWPC4VCOnTokHp6esouu+FbXrKDbot1Dfl8Pvn9/g0r+0YikYbqnwMAAOZo2JaXcDgswzCUTCZz38fjcfX19W3YX2hycjK3q3R2TZZgMMgGigAANKGGDS+Vrq+yfmNGAADQnBq+2wgAAGAtwosFmG0EAEB5WnK2USNithEAAOWpZrYR4QUA0DJ+NPlt3by+qC3bu/S+wAP1rg6qRHgBALSMH8We1+LMVXXtvI3w4mCMeQEAAI5CeAEAAI5CeLEAs40AACgPs40aBLONAAAoT1PubQQAALAW4QUAADgK4QUAADgK4QUA0BIyq6tKL69IktLLK8qkV+tcI1SL8AIAaHqXX/qR/q/P/bFuzi9Ikm7OL+gbn/6ips++XOeaoRqEFwBAU7v80o/0/af/Ukupa3mP37h8VS+cfIYA40BMlbaA1+tVW1ubQqGQQqFQvasDAC0rk17Vjye/XfKYl6N/rrs//H652p3793wqldKpU6eUSqXkdrslSaOjo/L5fIrFYnK73fL5fBvOOfkffk/z8/Pq7e3Vtr7eTc+xQjQa1ZNPPlnROc59pxpYIpHQq6++SnABgDqbecXY0OKy3o3LKc28YthUI/Mlk0ndc8896u/v1/j4uCKRiCKRiFKplAzDUDweL3pOn3uH/uB3HtPJ3/ytTc+xSigUUiKRqOgcwgsAoGktzs6belwjevDBB+Xz+XTixIm8x/1+v2KxmJ577rmi5/wv//bflX1OIyG8AACaVldfr6nHNZpwOKxUKqVIJFLw+ZGREaVSqZrPaTSEFwBA09r5AY863aWXnd92u1s7P+CxqUbmmpqaktvtlsdTuP579+415ZxGQ3gBADQtV3ubfirwQMlj7g190rGDdZPJZMVho5pzGo0z3y0AAMp0+33v0/sPf2JDC8y22926/+QhDey7t041Q7WYKg0AaHq33/c+9f/MsL7zW3+km/ML2tLbrY9+5XOObXHJ8ng8MozKZkpVc06jcfa7BgBAmVxtbWrfeutv9vatHY4PLpLk8/lkGEbRAbazs7OmnNNonP/ONSCv16s9e/YoGo3WuyoAgCYWiUTkdrt16tSpgs+Pj4+bco6VotGovF5vRecQXizAInUAADu43W5NTk5qYmJCU1NTec/FYjH19/eXPOfb3zlb1jlWqmaROsa8AADgYD6fTxcuXFA4HFY8Htfw8LBSqZR8Pp9GRkYUDoeLnnMsGNLz3zmrod3v1c2utpLnNBLCCwAADud2uwt295QamOt2u3X6dx6TS1JGUt9P7970nEZBtxEAAHAUwgsAAHAUuo0AAC3jff6P6Ob1RW3Z3lXvqqAGhBcAQMt43yZbBTSb7FoulWy0WM05dqPbCACAJjQ6OqpAICBJOnz4sEZHRzcNJNWcUw+0vAAA0ITi8bgt59QDLS8WYIVdAADKU80Ku7S8WCCRSKi7u7ve1QAAoOGFQiEdOnRIPT09mx/8j2h5AQAAjkJ4AQBYIpPJ1LsKqIN333fr3n/CCwDAVG1ttz5aVldX61wT1EP2fV9dXrHsGoQXAICptmzZovb2di0sLNS7KqiDxcVFtbW1KX1tybJrEF4AAKZyuVx6z3veo/n5ebqOWtC1a9e0fft2KW1dyxvhBQBguttuu003b97U66+/ToBpIXNzc7p+/bp6e3stvQ5TpQEAptu+fbsGBwd16dIl3bhxQ729vdq+fbva29vlcrk2HL+8uiKXbg3xXFxclHRr7MTy8rIymUxuHE01CpVtJivLd0LdM5mMVlZWdPXqVb3zzjvasWOHbrvtNlPruR7hBQBgife85z1673vfq6tXryqVSunKlStFj73+1mzu6znd+hDNfih2dHQUDDzlKlS2maws30l17+zs1J133qkdO3bUWq1NEV4AAJbZvn27tm/frl27dunmzZtFZyD9zReezbUA/NIzn5UkLS0taXp6WgMDA+rs7Ky6DoXKNpOV5Tul7u3t7TWHzEoQXgAAlnO5XNq6dWvR529enlebpFVJXV1ducfb2trU2dmZ91ilipVtFivLd3LdrcSAXQAA4CiEFwAA4CiEFwuwqzQAAOVhV+kGwa7SAACUh12lAQBA0yO8AAAARyG8AAAARyG8AAAARyG8AAAARyG8AAAARyG8AAAARyG8AAAARyG8AAAAR2GFXQAAYIrX5+Y1t7B46+ttW+SSlJH0g0tvS5J2dHfp7h29NV+H8AIAQAuxKmC8Pjevjz32p1peSd964Kfvyj33h3/wVUnS1o52/fVnf7XmAEN4AQCgRVgZMOYWFt8tt4jllbTmFhZrDi+MeQEAoEVUEjAaGeEFAAA4Ct1GAAA0ELsGvToZ4QUAgAZh56BXJyO8AACajtWtF1aVb+egVycjvAAAmorVrRe0jtQfA3Yt4PV6tWfPHkWj0XpXBQBajtUzapplxo7ZdnR3aWtHe8ljtna0a0d3V95j0WhUXq+3omvR8mKBRCKh7u7uelcDABraZl0v2ztc9atck8oGjFLhq1DAKMfdO3r115/91dx7+s3P/H7uPX3wy/8+d/31rVGhUEiHDh1ST09P2dcivAAAbFdu18tTgZ/XUB3q16yqDRiVlJ8994c3bqpN0qqkfzZ4R+2VX4PwAgCwXbldL/OLyzbVqHXYFTCsxJgXAADgKLS8AAAKYrE0+1k5JqWZEF4AABswHbg+rB6T0iwILwCADZy8WJrVrRdWl98MY1KsRngBADQVO2bU0DpSX4QXAEDTsbr1gtaR+iK8AICDOXVQbbldL71dW22sFZyC8AIADuXkQbXldL1s73BpKXWlfpVEwyK8AIBDOXlQrbR518vi4qIuEl5QAIvUAQAARyG8AAA2qHaHYMAOdBsBADZgOjAaGeEFACzk1NlAEtOB0bgILwBgESfPBgIaGWNeAMAilcwGAlA+wgsAOBSDatGq6DYCAIdiUC1aFeEFAByMQbVoRXQbAQAAR6HlZROpVEqnTp3SBz/4QUnSyMiIPB5PnWsFAEDrIryUYBiGAoGAEomEJCkYDCoej2t8fLzONQMAoHURXkoIh8MKBoN53wNoPlYtJJedDVRqujSzgYDKEV5KiMViikQiSiaTkm51GQFoLlYuJMdsIMAaDR9eUqmUDh8+rIMHD8rv9xc9bmxsTFeuXFF/f7/Onz+v0dHRksdvJhtYkslkLrQEAgE9/fTTcrvdVZcLoLFUspBcNSGD2UCA+Ro2vAQCAfX19Um61QJy8ODBoscGg0ENDw8rEonkHhsdHdXs7KyOHDlS1fUNw8h9nR2g6/F4FA6HGfMCAEAdNWx4mZyclHQrRExMTBQ9LplMamJiQplMJu/xSCQir9dbdXjJtq6s7Sr64Ac/qMOHDxNeAACoo4YNL+UaHx8vOBYl+1gsFst1H8ViMcXj8ZLlZQNPtrUl2/qTlUqlTKg1AAColuPDy9TUVNGBtG63W/F4PBde/H5/2eNgPB6P3G63DMPIlT87O8saLwAA1Jnjw4thGPL5fAWf6+vr07lz56ou++GHH84LR/F4vKzp0m+//ba6u7sLPre4aO3usVaW34p1X1pa0srKipaWliwpvxyt+LrbWf7yJu/t2uPM/Hmc8NrUo+y15Zd7/1VTtlWa4XW3q+ylpaXce7uwsFBReY4OL5t14bjd7pq6eU6cOKGxsTGNjY1JujUIuJwxNOtbZyZ/4TdyX1+8eLHq+pTDyvJbse4rKyu536GOjuK3SyPWvRHKd0LdX5+5Wt5xb7yubTcr+x9sKU54bepR9tryy73/qinbKs3wuttV9lNPPaVoNFpVeY4OL3Y4ceJExecYhpHX8vK3v/xY7uuhoSEzqpXn+2u+Nrt8K8u2unwzys7+VbB79251dnaaXn4xrf6621l+p/sdbe34+00XktvzUx7d5X5PTddy2mtjV9nFyi91/9Vatpma7XW3q+xHH31UjzzyiKRbLS+VDMtwdHjZbL2Veg2uveOOO4p2G3V1WbuSppXlt2rdOzo61NnZWbKMRq17vct3Qt3v2dVVl4XknPDa1KPs9eWXc/9VW7YVmuV1t6PstY+1VLfRZhhgC6AcLCQHOEtbvStQK7fbrdnZ2YLPpVIp7d271+YaAQAAKzk+vBw4cCBvNdz1RkdHbazNLV6vV3v27Kl6IBIAAK0iGo3K6/VWdI7ju40CgYAmJiaUSqXyxsBMTU1JUtFp1FZKJBJFx7wAAIB3hUIhHTp0SD09PWWf0/AtL9lBt8W6hnw+n/x+v06dOpX3eCQS0eTkJJsoAk3g9bl5/eDS2/rBpbf1+rYtmt62Ra9v2/LuY3Pz9a4iABs1bMtLOByWYRi53Z3D4bDi8bj6+vo27C00OTmpsbExhcPh3K7SwWCwpl2lATSG1+fm9bHH/vTdqcw/fVfuuT/8g69KujWV+a8/+6umzwgC0JgaNrys3SG6HNWsxwKg8c0tLJZcg0WSllfSmltYJLwALaLhu40AAADWIrxYgNlGAACUpyVnGzUiZhsBAFCeppxtBAAAsBbhBQAAOArhBQAAOArhBUBD29Hdpa0d7SWP2drRrh3d1u64C6BxMGDXAl6vV21tbQqFQgqFQvWuDhzi9bl5zS0s3vp62xa5JGUk/eDS25JufYi34jomd+/o1V9/9ldzr803P/P7udfmwS//e0mt+9oAzSAajerJJ5+s6BzCiwWYbVQaH9IbsYpsaXfv6M393D+8cVNtklYl/bPBO+paLwC1q2a2EeEFtnL6h7RVwYtVZAGgfIQX2MrJH9JOD14A0CwILxaha6T5ODl4AUAzIbxY4I25ef2rJ2P8hV4HhMb64HUHYCfCiwXmrvMXej3QrVMfvO4A7EZ4cSgr/9J16l/RdOvUB687ALsRXizw33/602r7FwcsK9/Kv3T5KxoAYKdq1nlhhV0L/KevfMXS8iv5S7eRykZxrCILoFWFQiElEomKzqHlBbbKfkiXCkit+CFtxyqyTu0OBID1CC+wlZOXerc6eFm5iizdgQCaCeEFtnPqUu9ODl4MqgXQTAgvFtixna6RZuXU4AUAzYTwYoG7HPwXupMxnqY+eN0B2I3wYhH+Qrefk7t1nIzXHYDdCC8OZOVfuk7/K5rQWB+87gDsRHixgNfrVVtbm0KhkEKhkOnlW/mXLn9FAwDsVM0idYQXCyQSCXV3d1t6DSv/0uWvaACAXUKhkA4dOqSenp6yz2GFXaAFsIIvgGZCywvQAugOBNBMCC9Ai6A7EECzoNsIAAA4CuEFAAA4CuEFAAA4CuEFAAA4CuEFAAA4CuHFAl6vV3v27FE0Gq13VQAAaGjRaFRer7eic5gqbQE7VtgFAKAZsMJuA8qkV0t+DwAAKkN4sdD02Zf1jU9/Mfcit0n6xqe/qOmzL9ezWgAAOBrhxSLTZ1/WCyef0Y3LV/Mev3H5ql44+QwBBgCAKhFeLJBJr+rl6NdLHvNy9M/pQgIAoAqEFwvM/ODChhaX9W5cTmnmFcOmGgEA0DwILxZYmn2nrOMWZ+ctrgkAAM2H8GKBzr73lHVcV1+vxTUBAKD5EF4ssPOf3aNtt99W8phtt7u18wMem2oEAEDzILxYwNXepntDnyp5zL2hT8rVzssPAECl+PS0yMC+e3X/yUMbWmC23e7W/ScPaWDfvXWqGQAAzsb2ABYa2Hev7v7w+/W1X/xNtUlalfTRr3yOFhcAAGrAp6jF1gcVggsAALXhk9QC7CoNAEB52FW6QbCrNAAA5WFXaQAA0PQILwAAwFEILwAAwFEILw63fmdqdqoGADQ7wouDTZ99Wd/49Bdzb2KbpG98+ouaPvtyPasFAIClCC8ONX32Zb1w8hnduHw17/Ebl6/qhZPPEGAAAE2L8OJAmfSqXo5+veQxL0f/nC4kAEBTIrw40MwrxoYWl/VuXE5p5hXDphoBAGAfwosDLc7Om3ocAABOQnhxoK6+XlOPAwDASQgvDrTzAx5tu/22ksdsu92tnR/w2FQjAADsQ3hxIFd7m+4NfarkMfeGPlnzDtasIQMAaESEF4ca2Hev7j95aEMLzLbb3br/5CEN7Lu3pvJZQwYA0KgILw42sO9effQrv61se8iqpI9+5XOmBBfWkAEANCrCi8Ot7xoyo6uINWQAAI2M8II8rCEDAGh0hBfkYQ0ZAECjI7wgj51ryDCbCQBQDcKLBbxer/bs2aNoNFrvqlTMrjVkmM0EAJCkaDQqr9db0TmEFwskEgm9+uqrCoVC9a5KxexYQ4bZTACArFAopEQiUdE5hBdsYOUaMsxmAgDUivCCgqxaQ4bZTACAWhFeUJTZa8hIzGYCANSO8AJbsSM2AKBWhBfYih2xAQC1IrzAVnbtiA0AaF58QsB2Vu+IDQBobh31rgBa08C+e3X3h9+vr/3ib6pN785mosUFALAZPilQN1bMZgIAND8+LQAAgKMQXgAAgKMQXgAAgKMQXgAAgKMQXgAAgKMQXtCU1u9KzS7VANA8CC9oOtNnX9Y3Pv3F3C93m6RvfPqLmj77cj2rBQAwCeEFTWX67Mt64eQzunH5at7jNy5f1QsnnyHAAEATILygaWTSq3o5+vWSx7wc/XO6kADA4QgvaBozrxgbWlzWu3E5pZlXDJtqBACwAuEFTWNxdt7U4wAAjYnwgqbR1ddr6nEAgMbErtIlTExMyOPxqK+vT263O/f4+u/RGHZ+wKNtt99Wsuto2+1u7fyAx8ZaAQDMRngpIRKJyDA2jo9IJBIaGRmpQ41Qiqu9TfeGPqUXTj5T9Jh7Q580ZffqQuvIsCs2ANiD/9uWEAwGlclkcv/Nzc0pEokQXBrYwL57df/JQ9p2+215j2+73a37Tx7SwL57a74G68gAQH0RXkrw+/153586dUonTpyoU21QroF99+qjX/ltZdtGViV99CufMy24sI4MANRXw3cbpVIpHT58WAcPHtwQJtYaGxvTlStX1N/fr/Pnz2t0dLTk8eXweN4dGzExMaFgMFhTebDP+i4cs7qKyllH5u4Pv58uJACwUMOGl0AgoL6+PklSLBbTwYMHix4bDAY1PDysSCSSe2x0dFSzs7M6cuRIzXVJpVKKx+OmlAXnqmQdmdt/9qdsqhUAtJ6GDS+Tk5OSJMMwNDExUfS4ZDKpiYkJZTKZvMcjkYi8Xq8pgWNiYkKjo6M1lwNnYx0ZAGgMDRteyjU+Pl5wAG32sVgslus+isViisfjJcsrFHjGx8c1Pj5uUo3hVHatI8NMJgAozfHhZWpqqujsH7fbrXg8ngsvfr+/qnEwhmHkurDQuuxYR2b67Mt6Ofr1DTOZ7g19ypQBxwDQDBwfXgzDkM/nK/hcX1+fzp07V3P5kipalO7tt99Wd3d3wecWFxdrqs9mrCy/Feu+tLSklZUVLS0tSZL+6eF/qeTvfbXo8f/08Me1dHNZuln5td78u+8XLDs7k2nkt35Fu/75+ysveJ31LTs3Fq5b1rLj5N8Zs8u383WXGvN+qrT89fefmWVbpRled7vKXlpayr23CwsLFZXn6PCSSqVKPu92uzc9phzZVXYrOX6tyV/4jdzXFy9erLk+pVhZfivWfWVlJfc71NHRIQ30aPDwqN6c/I5WUu/ebB07urXL/2EtDvRUda3M6qp+9B//ouQxL//hX+jGXdvlaqv+A2/+exf05uR38lp24v/jKe0KfFi9P3tP1eVmZVbzP6AvGEZN9d2M2b8zVtXf6te9EDNfm3q9rxvuPxPLtkoj/n+sEcovVPZTTz2laDRaVXmODi928Hg8On/+fEXnGIaR1/Lyt7/8WO7roaEhs6qW8/01X5tdvpVlW12+GWVn/yrYvXu3Ojs7s4XpA//tR/Sf/7vPqU231pEZ/dOHa/or+srLRl4YKmRlbkG916T+e4equsabf/d9vfr0xjFfK6kFXXo6XnPLzpt/9329OvGXeR/QF77wnPYc+YQpLUZZVv3OWFV/q1/3tax4ber5vha8/0wq20yN/v+xepW/WdmPPvqoHnnkEUm3Wl7W/+FfiqPDy2ZdOWa0ulTjjjvuKNpt1NXVZem1rSy/Veve0dGhzs7OkmVs695edfmStHqtvOba1WuLVf0smfSqfvj0fy55zA+f/iu994GRqkLY9NmXC3Z5Lc7MK/l7XzVtdeP1XS+dW7aa0vViVf2tft1LMeN+sut9XW9t3cu5/6ot2wqN+v+xepdfqOy1j1XabdTUUxhmZ2fZQBGOYPVMpkrWqKlUuYv3rQ8elbJqWwYr62/l6241u95XoBqODy9ut1uzs7MFn0ulUtq7d6/NNQIql53JVEotM5msXKPGjg9oK7dlsLL+dq4NVGiKfS3sDF5m1x3Nz/Hh5cCBAwV3fs6qx+JyXq9Xe/bsqXogElpPdkfsUmrZEdvKlh2rP6CtbgGwsv52rQ1kRauUXcHL6o1OrQ5GBK/aRaNReb3eis5xfHgJBAJKJpMbxrdMTU1JUtFp1FZKJBJ69dVXFQqFbL82nMvKHbGtbNlxcpeXZG39rW5Rk6xrlbIjeG1W9zf/7vtFziy/fCuDkdODV6MIhUJKJBIVndPw4SUbSop1Dfl8Pvn9fp06dSrv8UgkosnJSca8wFGs2hHbypYdJ3d5SdbW3+oWNStbpax+X8up+6sT/3nDNO1yWb0DvB3lWxmMJGeHo4YNL+FwWIFAQIFAIO/7Qjs7T05Oqr+/X+FwWGNjYwoGgwoGgzXvKg3UgxU7YkvWtew4uctLsr7+VraoWdkqZfXrUk7dF2eu6vqP36y4bKu7Gq0u3+pglL2G1eHISg07VXrtDtHlOHHihEU1AZrHwL57dfeH36+v/eJv5tao+ehXPldzQMp+QL8c/Xre/3C33e7WvaFPmtLlZeW2DFbWP1u+Fa+71a1SVr4u5dZpZf56xWVbvQO8leWXG4zu/vD7q/79yYaj9bLhyIwp8Fbv0dawLS8ArGFly47TurzWsqr+WVa87naMS7HqdSm3Th29la+hZHWoc/LsPTumwNvRqkN4sQCzjdCqnNbltZ5V9beKHQOCJWtel3Lq3rXzNm3/qV0Vl211qHPy7D2rw1E1XV4tOduoETHbCDCf1S0jTmRXq5QVyqn7niP/sqo9lKwOdU6evWdlOKq2VacpZxsBQJbTWkbsYFerlBU2q3u1eydZHeqcPHvPynBk58KG3PkA4HBObpWyqu5Whzqnzt6zMhzZuaI04QUAmoCTW6WcNojc6vKtDF5WhiO7VpSWCC8AgCZmdahzYvCyKhzZNYBcIrxYgtlGAIBaWRm8rAhH1bbqVDPbqGEXqXOyRCKh7u7uelcDAICirAhH1SxsGAqFdOjQIfX09JR9HcILAAAwjVUrSq9FtxEAADCV1WONCC8AAMBRCC8AAMBRCC8AAMBRCC8AAMBRCC8WYJ0XAADKwzovDYJ1XgAAKE8167zQ8gIAAByF8AIAAByF8AIAAByF8AIAAByF8AIAAByF8AKU8MYbb+ipp57SG2+8Ue+qAC2H+w/FEF6AEt58801Fo1G9+eab9a4K0HK4/1AM4cUCIyMjtixSt7S0pJMnT1pW7tLSUt3KqebcSs4x62dsNY38utldN6uuZ1a5jXrvVXsNVP+e2qGW97SaReoILxb49re/rVdffVWhUMjS6ywtLekLX/iCZeWaEV6qLaeacys5x6yfsdU08utmd92sup5Z5TbqvVftNVD9e2qHWt7TUCikRCJR0TmssGuSTCaT+/r69etaWFjIfb+Yvqk2SatS3uO1ypZldvnZMhYWFmoqe2057e3tBY8pVn455252vVJ1L7f8Gzdu5P5dX4ZV76vVZddSfi3vqVnM/J2ppPz1qrleOWXX8nOsLb/aMopd34z7qdTxhcovdf9VolHvp0rLlmr//S5Vfi11L/c93ex8Kf/ztBhXppyjsKm3335bd955Z72rAQCAo7311lu64447Sh5DtxEAAHAUWl5Msrq6qpmZGUnS9u3b5XK56lwjAACcIZPJ6Pr165KknTt3qq2tdNsK4QWogWEYGh8flyRFIpE61wZoHalUShMTE5KkeDyuQCCgI0eO1LlWsAsDdoEqTU1NKZVKyTAMeTyeelcHaCnhcDj3h8ORI0d0zz33qK+vT36/v841gx0Y8wJUyefzye/3q6+vr95VAVpKKpXSuXPnlEqlJElut1sHDhzQqVOn6lsx2IbwAgBwHMMwZBhG7nu32533PZob3UZoKalUSocPH9bBgwdLNi+PjY3pypUr6u/v1/nz5zU6OkpzNFAjs+4/t9utubm5vHOmpqa0d+9ey+qOxkJ4QUsIBAK57p1YLKaDBw8WPTYYDGp4eDhvAO7o6KhmZ2cZEAhUwer7L9sKMzk5aW7F0bCYbYSWYhiGhoeHNTk5WfAvv2QyKa/Xu2GFx2KPS7f+Z+t2u5ltBGzCivtPkrxeryKRiHw+nyX1RuNhzAuwxvj4uEZGRjY8nn0sFovZXSWgZVRz/4XDYYJLCyK8AGtMTU0VnfbsdrsVj8dtrhHQOiq9/yYmJjQ6OpoLLlNTU5bXEY2B8AKsYRhG0anPfX19OnfunM01AlpHJfdfdp0lj8ejZDKpZDLJHxcthAG7wD/KrhlRjNvtzjsmmUxqamoq99def3+/fD5fwWZvAKVVcv+lUimNjo5KutVtlMWA+tZBeAGqNDIyopGREZ04caLeVQFaitvtLjp4F62BbiPgH7nd7pLPb/aXIYDqcf+hEoQXoEyzs7Ob/g8WgDW4/7AW4QVYw+12a3Z2tuBzqVSKFTwBC3H/oVyEF2CNAwcOlNwfJTtIEID5uP9QLsILsEYgEFAymdzQv56dUcRCWIB1uP9QLsILWkr2f4rFmqZ9Pp/8fr9OnTqV93gkEtHk5CR97kANuP9gFvY2QksIh8MyDEPJZFKGYcjtdsvn86mvr0/j4+MbjmdXacA83H8wG+EFAAA4Ct1GAADAUQgvAADAUQgvAADAUQgvAADAUQgvAADAUQgvAADAUQgvAADAUQgvAADAUQgvAADAUQgvAOoiFotpeHhYw8PD7BZcQDAYzL0+Y2Nj9a4O0FA66l0BAK0rGAzqxIkT9a5GQ8ru+ROLxfTiiy/WuTZAYyG8AKjaxMSEzp8/L0nq7+8niACwBeEFQFUCgYA8Ho8ikYgkaXR0VCMjI/L5fHWuWf0NDw/L5/MV3DEZQO0Y8wKgYsFgUIZh5IKLJI2MjPBhvYbb7a53FYCmRcsLgIpMTU1pYmJCiUQi7/H+/n5NTEzUqVaNJduVBsAatLwAqEgwGNTIyIhGRkbyHr9y5YpSqVR9KgWgpdDyAqBsU1NTMgxDk5OTG54zDKMONQLQiggvAMqWHdPi9/s3PJdMJk0Z55FKpTQxMZErK5FI5Fp7sgzDUCAQkGEYOnDggCKRSK7LKh6PKx6P55WZnW7c39+vK1euaHh4WEeOHNlQ1t69exWPx3NlnT9/XslkUpFIRCMjI0omkzp37lyuXsPDwxtmWI2OjsowDLnd7lzXWqlrpFIpvfjii4pEIvJ4PDW/fkBLyABAGebm5jKSMj6fr6LnipmcnMxEIpENj584cWJD2W63OxOPxzcc6/P5Mn6/P1dOPB7PSMqcP38+rzy/37+hzPXX8fv9GZ/PlxkfH8/Mzc3l1dPtdmcSicSGOng8nsz4+PiGeh05ciQzMjKy4fHNrlHI5OTkhroCrY4xLwDKkm0pMAxDo6Ojef89+OCDkrRhHEylksmkYrFYXheU2+3WkSNHFA6HNxyfbQ3JTs/2+Xyam5vLtWBMTU1pbGxMTz/9dN55hmFsuI7H49G5c+e0d+/evBYkn8+nVCql8fHxDdPAfT5fwS604eHhgj/fZtdIJpNFXhkAa9FtBKAs2VVe4/H4hu6NYDCoZDKpgwcP1nQNt9ut2dlZGYaRd43h4eGiM5lmZ2fzQtPaUBAOh+X3+zd0Z83Ozmp2djbvsf7+fqVSqQ0BLHtuoUCSrW+5NrtGJWUBrYzwAqAs2TEthcZlPPfccwVnIFXK4/Fobm4u971hGEqlUhumZa8/p1SdCy2al22hqaSsWn+2cq4BoDx0GwEoS3aw6XrJZFKpVEoPP/ywKddJpVIKh8MKBoOampqS2+2W1+stenxfX1/R+kq3WjvKZcfCcixeB9SOlhcAZSvUanDq1Cm53e6CM5AqZRiGvF6vIpFIbjaQVDyglJKt65UrV2quF4DGQssLgLK43e4NrQbZga/rB8RWK7tf0trgIm0cCzI1NVVWeR6Pp+T6MyyqBzgT4QVAWQp1GQWDQfn9flNaXaTiY1QSiURe0Ch3QbxIJKJYLFYwpBiGUXYIAtBYCC8AyhIOh/Om8o6NjUlSwanC1fL5fBsCRXZqtqTcdOJskEqlUiVn6Pj9fh05ckSBQGDDc7FYLC90bba9QaWtNIWON/saQKtizAuAsvh8PiWTydx6K/39/RtWsq1VPB5XMBhUIBDIBRaPxyO/369IJKJAIKBAICCfz6dAIKCpqSmlUqlcd9PaXa6zxsfHFYvFFAwGNTw8LI/Ho9nZ2dzKuKlUSocPH84r6+DBg/L7/RobG8v9jOFwWM8++6wmJydlGIbC4XDeOcFgUHv37s0ra3R0NLdybqXXAFCcK5PJZOpdCQCtJ7tI3Prl9ZEvu7VBoWAGtCq6jQAAgKMQXgAAgKMQXgAAgKMQXgAAgKMQXgDUzalTpzQ8PJybWYR3ZWdHHT58uN5VARoOs40AAICj0PICAAAchfACAAAchfACAAAchfACAAAchfACAAAchfACAAAchfACAAAchfACAAAc5f8H6mIR3NGpRREAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 600x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "corpp = fitsio.read(\"../psf_sys/star_gpgp_used.fits\")\n",
    "corpq = fitsio.read(\"../psf_sys/star_gpgq_used.fits\")\n",
    "corqq = fitsio.read(\"../psf_sys/star_gqgq_used.fits\")\n",
    "\n",
    "\n",
    "plt.close()\n",
    "fig=plt.figure(figsize=(6,7))\n",
    "ax1=fig.add_subplot(1,1,1)\n",
    "\n",
    "ax1.errorbar(corpp[\"r_nom\"],corpp[\"xip\"],np.sqrt(corpp[\"sigma_xip\"]),fmt='o',label='pp',color='black')\n",
    "msk=(corpq[\"xip\"]>0)\n",
    "ax1.errorbar(corpq[\"r_nom\"][msk],corpq[\"xip\"][msk],np.sqrt(corpq[\"sigma_xip\"])[msk],fmt='o',label='pq',color='#28728f')\n",
    "ax1.errorbar(corpq[\"r_nom\"][~msk],-corpq[\"xip\"][~msk],np.sqrt(corpq[\"sigma_xip\"])[~msk],fmt='s',color='#28728f')\n",
    "\n",
    "msk=(corqq[\"xip\"]>0)\n",
    "ax1.errorbar(corqq[\"r_nom\"][msk],corqq[\"xip\"][msk],corqq[\"sigma_xip\"][msk],fmt='o',label='qq',color='#ad4271')\n",
    "ax1.errorbar(corqq[\"r_nom\"][msk],-corqq[\"xip\"][msk],corqq[\"sigma_xip\"][msk],fmt='s',color='#ad4271')\n",
    "\n",
    "ax1.set_yscale('log')\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_ylabel(r'$\\xi_+$')\n",
    "ax1.set_xlabel(r'$\\theta$ [arcmin]')\n",
    "ax1.set_xlim(3,360)\n",
    "ax1.legend()\n",
    "ax1.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([(2.82214539e-01, 2.85500638e-01, -1.25587868, 9.30272218e-04, -2.32540829e-06,  6.39782553e-23, -3.35039636e-06, 1.09403318e-06, 1.09403318e-06, 4.06714000e+05, 4.06714000e+05),\n",
       "       (3.59631839e-01, 3.63851634e-01, -1.01334572, 9.26764557e-04,  1.22387671e-06, -8.37356249e-25,  1.14988747e-06, 8.67141653e-07, 8.67141653e-07, 6.47396000e+05, 6.47396000e+05),\n",
       "       (4.58286309e-01, 4.63669940e-01, -0.77089311, 9.15119753e-04, -4.97102290e-07, -5.36452617e-23, -1.66603792e-06, 6.85883920e-07, 6.85883920e-07, 1.03478200e+06, 1.03478200e+06),\n",
       "       (5.84003746e-01, 5.90913677e-01, -0.52836451, 9.08797432e-04,  6.68559935e-07, -1.05633098e-22,  1.26168165e-06, 5.41750410e-07, 5.41750410e-07, 1.65863800e+06, 1.65863800e+06),\n",
       "       (7.44208083e-01, 7.53007413e-01, -0.28592465, 8.99160621e-04, -2.25365610e-07,  2.35116851e-23, -4.79432647e-07, 4.28105250e-07, 4.28105250e-07, 2.65612800e+06, 2.65612800e+06),\n",
       "       (9.48359791e-01, 9.59768293e-01, -0.04327767, 8.90677689e-04,  5.97492749e-08, -1.10471970e-22,  2.31845434e-07, 3.37599289e-07, 3.37599289e-07, 4.27117200e+06, 4.27117200e+06),\n",
       "       (1.20851455e+00, 1.22328890e+00,  0.19935156, 8.79730508e-04, -4.47425843e-07, -2.71628657e-22,  1.92085742e-07, 2.66655365e-07, 2.66655365e-07, 6.84619800e+06, 6.84619800e+06),\n",
       "       (1.54003514e+00, 1.55921600e+00,  0.44200658, 8.71956512e-04, -1.02449725e-06, -9.35845556e-22, -3.49359572e-07, 2.10620754e-07, 2.10620754e-07, 1.09735660e+07, 1.09735660e+07),\n",
       "       (1.96249872e+00, 1.98710302e+00,  0.68450884, 8.63179543e-04, -9.21785010e-07,  2.77720116e-22, -4.79830454e-08, 1.66359371e-07, 1.66359371e-07, 1.75895840e+07, 1.75895840e+07),\n",
       "       (2.50085281e+00, 2.53254600e+00,  0.92706495, 8.52880634e-04, -1.17899466e-07, -7.09204363e-22,  2.02367324e-07, 1.31437529e-07, 1.31437529e-07, 2.81780760e+07, 2.81780760e+07),\n",
       "       (3.18688859e+00, 3.22771421e+00,  1.16962127, 8.40837724e-04, -8.90110274e-08, -1.73570582e-21,  6.41637091e-07, 1.03973262e-07, 1.03973262e-07, 4.50305000e+07, 4.50305000e+07),\n",
       "       (4.06111821e+00, 4.11384657e+00,  1.41220878, 8.25579106e-04,  5.54906664e-08, -1.48330241e-21,  2.46033819e-07, 8.23093504e-08, 8.23093504e-08, 7.18541340e+07, 7.18541340e+07),\n",
       "       (5.17516714e+00, 5.24250603e+00,  1.65464826, 8.07857366e-04,  4.44441913e-07, -1.81374827e-21,  2.07443986e-07, 6.51786742e-08, 6.51786742e-08, 1.14587980e+08, 1.14587980e+08),\n",
       "       (6.59482279e+00, 6.68100845e+00,  1.89711822, 7.89419274e-04,  4.54296846e-07,  1.92853731e-21, -1.29952756e-07, 5.16643982e-08, 5.16643982e-08, 1.82375862e+08, 1.82375862e+08),\n",
       "       (8.40391943e+00, 8.51596484e+00,  2.13980259, 7.64830133e-04,  1.59270626e-06, -4.36965828e-21,  5.15798048e-07, 4.08956555e-08, 4.08956555e-08, 2.91068870e+08, 2.91068870e+08),\n",
       "       (1.07092888e+01, 1.08535564e+01,  2.38234724, 7.38169865e-04,  2.69208208e-06, -1.66473642e-20,  4.35978298e-07, 3.23177690e-08, 3.23177690e-08, 4.66087400e+08, 4.66087400e+08),\n",
       "       (1.36470687e+01, 1.38306846e+01,  2.62474878, 7.07615201e-04,  4.47149238e-06, -4.36121527e-20,  6.66250533e-07, 2.55588652e-08, 2.55588652e-08, 7.45189920e+08, 7.45189920e+08),\n",
       "       (1.73907427e+01, 1.76264570e+01,  2.86726897, 6.64666360e-04,  6.10314598e-06, -1.46358037e-19,  1.59431552e-07, 2.02312597e-08, 2.02312597e-08, 1.18933516e+09, 1.18933516e+09),\n",
       "       (2.21613840e+01, 2.24632618e+01,  3.10975037, 6.14594680e-04,  7.32728489e-06,  1.01103361e-19, -9.21187132e-07, 1.60300555e-08, 1.60300555e-08, 1.89443652e+09, 1.89443652e+09),\n",
       "       (2.82407111e+01, 2.86247112e+01,  3.35211641, 5.58966495e-04,  9.18181858e-06,  4.27575215e-20, -1.61152303e-06, 1.27287865e-08, 1.27287865e-08, 3.00452697e+09, 3.00452697e+09),\n",
       "       (3.59877238e+01, 3.64665841e+01,  3.59426351, 5.02713204e-04,  1.06736242e-05, -2.42670969e-21, -2.05308328e-06, 1.01430411e-08, 1.01430411e-08, 4.73166224e+09, 4.73166224e+09),\n",
       "       (4.58599026e+01, 4.64688303e+01,  3.83664733, 4.40281395e-04,  8.26045833e-06, -4.58325658e-20, -2.35087535e-06, 8.11383288e-09, 8.11383288e-09, 7.39431575e+09, 7.39431575e+09),\n",
       "       (5.84402248e+01, 5.92067672e+01,  4.07887761, 3.88290097e-04,  9.45737214e-06,  8.09901342e-20, -2.18160468e-06, 6.54324995e-09, 6.54324995e-09, 1.13700686e+10, 1.13700686e+10),\n",
       "       (7.44715902e+01, 7.53944127e+01,  4.32058482, 3.48289263e-04,  4.51898319e-06, -1.52184496e-19,  1.62229170e-06, 5.32188802e-09, 5.32188802e-09, 1.71877348e+10, 1.71877348e+10),\n",
       "       (9.49006916e+01, 9.60228602e+01,  4.56244034, 2.83283477e-04, -3.66447163e-06, -8.62318233e-22,  6.07721209e-06, 4.38533859e-09, 4.38533859e-09, 2.53130132e+10, 2.53130132e+10),\n",
       "       (1.20933919e+02, 1.22309126e+02,  4.8043821 , 2.39898389e-04,  5.24022035e-06,  4.17037735e-20,  2.58880430e-06, 3.68927238e-09, 3.68927238e-09, 3.57658596e+10, 3.57658596e+10),\n",
       "       (1.54108601e+02, 1.55666368e+02,  5.04554451, 2.03460029e-04,  8.70759031e-06, -8.84783460e-20, -3.18793030e-06, 3.18105032e-09, 3.18105032e-09, 4.81070838e+10, 4.81070838e+10),\n",
       "       (1.96383785e+02, 1.98272171e+02,  5.28748398, 1.57617274e-04,  1.91248778e-05,  3.88029192e-20, -1.02643827e-05, 2.77142942e-09, 2.77142942e-09, 6.33785737e+10, 6.33785737e+10),\n",
       "       (2.50255929e+02, 2.52466720e+02,  5.52912448, 1.05929680e-04,  2.64874020e-05,  1.72091057e-20, -1.53292527e-05, 2.48135168e-09, 2.48135168e-09, 7.90630322e+10, 7.90630322e+10),\n",
       "       (3.18906320e+02, 3.21168492e+02,  5.76979058, 8.73475694e-05,  3.24354939e-05,  1.59150341e-20, -2.17947666e-05, 2.33080384e-09, 2.33080384e-09, 8.96063231e+10, 8.96063231e+10)],\n",
       "      dtype=[('r_nom', '>f8'), ('meanr', '>f8'), ('meanlogr', '>f8'), ('xip', '>f8'), ('xim', '>f8'), ('xip_im', '>f8'), ('xim_im', '>f8'), ('sigma_xip', '>f8'), ('sigma_xim', '>f8'), ('weight', '>f8'), ('npairs', '>f8')])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corpp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}