jmeyers314/GalSim vs Analytic Gaussian Moments.ipynb

## GalSim vs Analytic Gaussian Moments.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:02:58.627017Z",
     "start_time": "2019-01-24T21:02:57.936980Z"
    }
   },
   "outputs": [],
   "source": [
    "import galsim\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:02:58.631909Z",
     "start_time": "2019-01-24T21:02:58.629197Z"
    }
   },
   "outputs": [],
   "source": [
    "galMom = dict(Ixx=3.0, Iyy=5.0, Ixy=-2.0)\n",
    "psfMom = dict(Ixx=1.0, Iyy=1.0, Ixy=0.0)\n",
    "convMom = {k:galMom[k]+psfMom[k] for k in galMom.keys()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:02:58.635901Z",
     "start_time": "2019-01-24T21:02:58.633402Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'Ixx': 4.0, 'Iyy': 6.0, 'Ixy': -2.0}\n"
     ]
    }
   ],
   "source": [
    "print(convMom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:02:58.641455Z",
     "start_time": "2019-01-24T21:02:58.638299Z"
    }
   },
   "outputs": [],
   "source": [
    "# Super conservative GSParams\n",
    "gsp = galsim.GSParams(\n",
    "    folding_threshold=1e-15,\n",
    "    stepk_minimum_hlr=200,\n",
    "    maxk_threshold=1e-15,\n",
    "    kvalue_accuracy=1e-15,\n",
    "    xvalue_accuracy=1e-15,\n",
    "    realspace_abserr=1e-15,\n",
    "    realspace_relerr=1e-15,\n",
    "    integration_abserr=1e-15,\n",
    "    integration_relerr=1e-15,\n",
    "    shoot_accuracy=1e-15\n",
    ")\n",
    "# or uncomment for default GSParams\n",
    "# gsp = galsim.GSParams()\n",
    "\n",
    "# put draw kwargs here so can try different variations quickly\n",
    "drawKwargs = dict(\n",
    "    method='fft',\n",
    "    nx=64, ny=64,\n",
    "    scale=0.2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:02:58.646832Z",
     "start_time": "2019-01-24T21:02:58.643779Z"
    }
   },
   "outputs": [],
   "source": [
    "def momToSizeE(mom):\n",
    "    T = mom['Ixx'] + mom['Iyy']\n",
    "    e1 = (mom['Ixx'] - mom['Iyy'])/T\n",
    "    e2 = 2*mom['Ixy']/T\n",
    "    # This size is invariant under shear, T above is *not*\n",
    "    size = np.sqrt(np.sqrt(mom['Ixx']*mom['Iyy']-mom['Ixy']**2))\n",
    "    return size, e1, e2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:02:58.652484Z",
     "start_time": "2019-01-24T21:02:58.648058Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8211602868378718 -0.25 -0.5\n",
      "1.0 0.0 0.0\n",
      "2.114742526881128 -0.2 -0.4\n"
     ]
    }
   ],
   "source": [
    "galSize, galE1, galE2 = momToSizeE(galMom)\n",
    "psfSize, psfE1, psfE2 = momToSizeE(psfMom)\n",
    "convSize, convE1, convE2 = momToSizeE(convMom)\n",
    "\n",
    "print(galSize, galE1, galE2)\n",
    "print(psfSize, psfE1, psfE2)\n",
    "print(convSize, convE1, convE2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:02:58.657180Z",
     "start_time": "2019-01-24T21:02:58.653839Z"
    }
   },
   "outputs": [],
   "source": [
    "convGalSim = galsim.Convolve(\n",
    "    galsim.Gaussian(sigma=galSize, gsparams=gsp).shear(e1=galE1, e2=galE2),\n",
    "    galsim.Gaussian(sigma=psfSize, gsparams=gsp).shear(e1=psfE1, e2=psfE2),\n",
    "    gsparams=gsp\n",
    ")\n",
    "\n",
    "convAnalytic = galsim.Gaussian(\n",
    "    sigma=convSize,\n",
    "    gsparams=gsp\n",
    ").shear(e1=convE1, e2=convE2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:03:03.333511Z",
     "start_time": "2019-01-24T21:02:58.658571Z"
    }
   },
   "outputs": [],
   "source": [
    "imgGalSim = convGalSim.drawImage(**drawKwargs)\n",
    "imgAnalytic = convAnalytic.drawImage(**drawKwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:03:03.602245Z",
     "start_time": "2019-01-24T21:03:03.334924Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0cAAADSCAYAAACB+bvuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXu8bFV15/v9Ve19zgERjoAPnoGEYxJMt6ZzLiYxDyOxecR40mlsQdtgGtvbHUnSeQoxH7Xp0Bc7t/WSoDG0EMGIQIjE0wbFV7zERHlo8AFIPAHUI0TeCMI5Z++q0X/MuWrNqppVtfauqr1X7Rrfz6c+tWquueaaVbX3qLXG/I0xZGY4juM4juM4juPMO431noDjOI7jOI7jOE4d8Jsjx3Ecx3Ecx3Ec/ObIcRzHcRzHcRwH8Jsjx3Ecx3Ecx3EcwG+OHMdxHMdxHMdxAL85chzHcRzHcRzHAfzmqAtJx0gySQsV+v6kpDvXYl5rjaS3SvrzMY7/iKQzJzknx3G6Gef/dKPYL0mflvS6MY5/QtL3TnJOjuOUzOt1laT3SvqDVR77akkfm/ScnOpsuJsjSadLulHSdyXdH7d/RZJWMdbzJH1M0iOSHpX0eUmnApjZ35rZ90/+HcwWuQs0MzvFzC5brzk5Tt2IF/GPSNq8Tuc3SccVr+fRfuVupMzsADO7a73m5DizwEa+rlLgLkm3r+V5k/P33Tya2fvN7F+vx3ycwIa6OZL0W8CFwB8CzwGeDfwn4EXAplUM+b+Bj8dxngX8GvCdiUzWcZy5QNIxwE8CBrx8XSfjOI6zAubguuqn4jy+V9L/tY7zcGrEhrk5knQQcB7wK2Z2jZk9boF/MLNXm9ne2O/nJP2DpO9I+qaktw4Y71DgWOB/mdm++Pg7M/tM3P9iSbuT/vdI+h1JX4relUskPTtKzB6X9AlJzxgy/x2Sbo3z+idJJ8f2wyXtlPSwpF2S/mNyzFslXS3p8niO2yRtj/vOkXRNzzkulPRHo8btOabrfSbv9WfjHH8PeGWUp3wx7u94aCU1JP2+pK9Hj9Pl8btKPSZnSvqGpAclvWnQZ+Q4M8ovAZ8D3gt0yU2j9OKdkv46/g/fKOn7kv0XRjv1nehh/cncCeLxv9rT9iVJvyDphtj0xfh/+sqM/TpK0gclPSDpIUkXDThPU9LvRRv1eJzTUXHfj0u6WdJj8fnHk+M+Lem/Sfq7eNzHoo1F0kclnd1zni9K+sVR4/Yc07WKnXpkJZ1PuEG9KH4GF8U+nRU1SQdF+/RAtFe/L6kR971W0mck/b8KHu+7JZ2Sm4fjbBQ049dVFTkT+BBwHf32eaDdivv/QtI/R9t0g6TnDXjfX5H088nrxXi98wKgsM+PRtv0Y4W9Sfo/T9LHFa7Xvi3p98Z8z84INszNEfBjwGbCH/kwvku4WNkK/BzwnyX9QqbfQ8Au4M/jBcazK8zh3wIvBZ4L/DzwEcLNw6GEz/rXcgdJOgG4HPidOK+fAu6Juz8A7AYOB04D/rukE5PDXw5cGY/bCVyUHHeqpAPjOZrAvwOuqDjuSMzso8B/B66K8pTnZ7q9Nj5+Bvhe4IBkjgU/AXw/cCLwZkk/uJJ5OE7N+SXg/fFxUsaWnAH8V+AZBJtzfrLvZuAFwMGE/92/kLQlc47LgH9fvJD0fOAI4Doz+6nY/Pz4f3pVemC0DR8Gvg4cE4+7csB7+c0431OBA4H/ADwp6WDgr4E/Ag4B3g78taRDkmNfBfwywUu7Cfjt2H5FHLOYz/HA98Tjq4w7EjN7E/C3wNnxMzg70+2PgYMIduqnCd/bLyf7XwjcSbDn/wO4RFq5rMhxZoiZva6qgqT9Cdc/hX0+XVLvatggu0Wcy7a47wtxjByXk9hngv28z8xuJVzvAWyNtumzPXN8OvAJ4KOE67XjgE+u4G06q2Aj3RwdCjxoZstFg6S/V9C0PiXppwDM7NNm9mUza5vZlwg3CT/dO5iZGeGC/h7gfwL3Rc/AtiFz+GMz+7aZfYvwQ3xj9LDsBa4FfnjAcWcBl5rZx+O8vmVmX40e2Z8A3mhme+I/0nuA1yTHfsbMrjOzFvA+4Plx/l8n/LMWBuolwJNm9rmK406KVwNvN7O7zOwJ4FyCAUqDM/+rmT1lZl8Evli8B8eZdST9BOFC/2oz+zzwT4Qf25QPmtlN0Xa9n3AzBICZ/bmZPWRmy2b2PwkXKjlN/oeAbYl9eg3BabGvwjRPIPzo/o6ZfTfahM8M6Ps64PfN7M7oQf6imT1EuCD6mpm9L871A8BXCRczBX9mZv9oZk8BVyfv81rgBZK+J75+dfxM9lYcd2ziDeIrgXOjd/wegt1PbeLXzex/RVt7GXAYQRrkOBuVWb6uqsIvAnuBjxEcRAsEm5MyyG5hZpdGe7EXeCvw/Lja1sufkzirCXblfRXn+DLgn83sf0bb/LiZ3Vjx2JFIulRB1fOVCY330fj38eGe9mMVlBFfk3RV5ia0Vmykm6OHgEPVHdT242a2Ne4r5BEvlPQ3UTrxGEE7e2huQDPbbWZnm9n3ES5wvkvwAAzi28n2U5nXBww47ijCRVMvhwMPm9njSdvXCZ7dgn9Otp8EtiSfQeqRfRXlqlGVcSfF4XHs9DwLdF9U9L6HQZ+T48waZwIfM7MH4+sr6JFuMOTvX9JvSbojyjYeJaxs9Nmr+ON8NfDvoxTsDKr/+B5FuPBfHtlzuK36ek/bKFt1AEC0Q38NnB73nU7pga0y7iQ4lOAV7rVV2fmb2ZNx022Vs5GZ2euqKL17Ij5ePWDsMwmOq+VoQz9IRfusIDG+QEFi/B1KtU/OPt8L/B3wbyVtBU5h8CpTL4Ns7qR4L3DyBMf7Q/KO9rcB7zCzbcAjhEWB2rKRbo4+S/AA7BjR7wqC/OwoMzsIeDcwUhphZt8E3gn80JjzzPFN4Psy7fcCB8dl1YKjgW9VHPcvgBdLOhL4N5Q3RysZ97vA/sWL6GF9ZrLfRszhXoIBTM+zTLeBc5wNh6T9CFLWn4669H8GfoPgXRy5OqoQX/TGOMYz4gXJYwy2V5cRVl1OJKwSf3ZAv16+CRytCql2GW6rvqenbSW26gPAGZJ+DNgP+JtVjNtlqwjB4ynDbNWDwBL9tqrq/B1nIzKz11UWsuYeEB99NyLxuuglBIdSYZ9PI6zwZG/sengV4XP5WYLT6phi6AH9C+nzK4DPxpUwGH0NNcjmTgQzuwF4OG2T9H1xBejzkv5W0g+sYLxPAqnjnSg/fglQxMFfRqlqqiUb5ubIzB4l6PbfJek0SQcoJAN4AfC0pOvTCasme2KsT6/EBQBJz5D0XyUdF8c5lKCv/9wUpn8J8MuSToznOkLSD0TD8ffA/yNpi6R/SbjbruRxMLMHgE8DfwbcbWZ3xPaVjPuPhNWon5O0CPw+QdpT8G3gmOitzvEB4DfikuoBlDFKVbzUjjPL/ALQAo4nSDFeAPwgQRrySxWOfzrBkfAAsCDpzYQ4nyzxZqhNkKv0rhp9mxBLk+Mm4D7gAklPizbhRQP6vgf4b5K2KfAvY/zPdcBzJb1KIQHCKwnv+8MDxunlOsKNyXkE+9BO2quOeyvwU5KOjtKWc3v2D/wMolTuauB8SU+PEr/fJMhhHGcumfHrqlG8hnB98/2U9vm5hFjsM4YcV/B0wo3jQwSnzH8f0f+vgH8F/DrdK2UPEOz2IPv8YeA5kv6LpM3RPr2wwvzG4WLgV83sRwgxVu8ac7xDgEeT677dTEepNDE2zM0RgJn9D8IP2u8C9xN+DP+U4H39+9jtV4DzJD0OvJnwg5hjH8ET8AlCmsmvEP4RXjuFed9ECPh7B8Ez/P9TejDPiPO4l6CvfYuZfXwFw19B8Gxc0dNeaVwze4zwmb2H4EX9LuEPu+Av4vNDkr6QOf+lhAu1G4C7gT3Ar2b6Oc5G40yCXv0bZvbPxYOQkOTVFVZqricE/P4jQeK1h+BFHMblwL+g/6L+rcBlUQv+79Id8cbg5wmBvt8g/H+/csD4byfYzI8R7OIlwH4x7uhlwG8RLhZ+F3hZIiccSiJp6bJVKxk32q+rgC8Bn6f/BupC4DSFbHN/lJnGrxLs213AZ+I8Lq0yf8fZqMzqdVUFzgTeldrmaJ/fTb+0LsflBLv8LeB2RtzgxZilvyRk6/tg0v4kIQnP30X7/KM9xz1OSEjx8wSJ39cIcVtTITqxf5yQ/OdWwnd9WNz3iwqZ93of148aNtM2asVsXVGIj3Mcx3FmHUm/BLzezH5ivefiOI7jlMTV/+ea2b8f2XkNUajF92Ez+yGFpBF3mtlhY4z3YuC3zexl8bUIK2TPMbPlKJ9+q5mdNPbkp8SGWjlyHMeZVxTS0v4KQRLhOI7j1ASFsgRnUXP7bGbfAe6W9AoINzZVYmRHjGmEONLTYlNRW6q2+M2R4zjOjCPpJIJn7tv0S2gdx3GcdULSfyRIoj8SEyDUBkkfICTe+H5JuyWdRUjsc5akLwK3MTohRzre3xLCLU6M4xWrQ28EflPSLkIM0iUVxztKIRPiHZJuk/TrmT6S9EeSdikUDP5XVec78Lwuq3Mcx3Ecx3Ecp05IOgw4zMy+oJBh+fPAL5jZ7UmfUwkxo6cSinVfaGZjJa3wlSPHcRzHcRzHcWqFmd1nZl+I248Dd9Cf6W4HcLkFPgdsjTdVq8ZvjhzHcRzHcRzHqS0xccQPAzf27DqC7kyuY6cKr1L0b2Js0mbb0pUaf0wGlhgbXHssu0cja5WtDRmJ42jRY6aHKyVHsofvss/2Vv7iT/qZp9lDD7cG7v/8l/Zeb2aTrDLtrCFum0bgtmnNWKltArdPG51DDz3UjjnmmPWeRj/tdrndGOFrL/qO6jfOPNKxp3m+Oeaee+7hwQcfrGyfjpPsySH77wsxTXuSpovNrC9pRUwx/pfAf4lJI7p2Z4Ye69dmTW+OtvA0XqgTuxur/vhnaoyqMeDYXD3S2Fe58w3655nmhUku1is1NJ1usV97wPdsmWNyfTP9Ks9rA3KjfXJF/R98eJm//+hgR8SWw++uUlHbqSlumxLcNq0rK7VN4PZpo3PMMcdwy003rfc0+tmTXNNu2VKt76h+48wjHXua55tjtp9wwor6PwW8Ycj+34c9ZrZ92BiSFgk3Ru83sw9muuwGjkpeH0mo4blq1vTmaFVUvfBI+2X2dy48chcb6YVGZn/2omVMsokwivMk+1RclDTLOXQd2859PrmLksz7zl2UpO91Ti5GqmBA293eTorbprDhtmndcfvkVGJ5OTwvTOjSbyU3HtO8ScmNXbSt5AZu0p+PgxgvfifWSLoEuMPM3j6g207gbElXEhIyPGZm941x2hm4OXKcGmAYSzZYtuI4jrNeuH1yHKeujHmj8SLgNcCXJd0a234POBrAzN4NXEfIVLcLeBL45fFOuZ43R8M8npP2yELpdc20dfcbcZ5h/XIMkJxoiGe0y/ua89h2HdPuP6bwxCbnKD6/LlmL+vvl5uNe2oB7ZucEt02ZRrdNdcftkzOS1ayI1EGe9uij5fbWrdWOGTXv3P7i8ylWkNI2Z1WMu3JkZp9hWLAunSKzw9R7K6bSnCVtlXSNpK/GQkw/JulgSR+X9LX4/IxJTsxx6oQBS7QHPpz1wW2T47h9qitun5x5R4RVmEGPulL1hu5C4KNm9gPA8wl5xs8BPmlm24BPxteOsyExoGU28OGsG26bnLnH7VNtcfvkzD2NIY+6MnJukg4EfooQEIWZ7TOzRwlFly6L3S4DfmFak+yfVCM8Giof3XMOcpRGo3xI4ZG0qRkeNJsjHo3ysbAQHo1m+SjaVvNIz7O4AIsLaCF5LIZHV7/ivSTvp3jP4X3HR9cx4TNTQ51H3+eZk+iED7R8zCmGsTTkUQVJJ0u6U9IuSX0/iJI2S7oq7r8x5vQv9p0b2++UdFLSfqmk+yV9ZcA5f1uSSdpw2arcNrltctsUGNc+uW2aPLW0T6thy5bpS+oefbR85Ni6tXxUZdS8h+1P7aAzFht55eh7gQeAP5P0D5LeI+lpwLOLbBDx+VlTnKfjrCtmsDTkMQpJTeCdwCnA8cAZko7v6XYW8IiZHQe8A3hbPPZ44HTgecDJwLvieADvjW25cx4FvBT4xore7OzgtslxGM8+uW2aGm6fnLlnVm+OqsxtAfhXwK+a2Y2SLmQFy8CSXg+8HmAL+4/ovIJ6If3nKV8UAc7NZv/+Xq8kBK9rbg65+iMrqUXSS1etkDi3RPLQCVxuJVmHrJhDco4kSNla7b55qVN/JOlXHJ6kzR2ZVndYOt25k2qI1vCYwFGcAOwys7sAYsrJHcDtSZ8dwFvj9jXARTGN5Q7gSjPbC9wtaVcc77NmdkPqxe3hHcDvAh8aZ+I1xm2T26aSubVNMKZ9cts0HSZmn44++mja0ZfdmIUYsmIFaNRKT9V+dSBN0lDgK0uVqLN8bhBV5rwb2G1mN8bX1xD+4b8t6TCA+Hx/7mAzu9jMtpvZ9kU2T2LOjrPmGLBkGviowBHAN5PXu2Nbto+ZLQOPAYdUPLYLSS8HvmVmX6wyuRnFbZPjMLZ9cts0HSZmn575zGeuyYQdZ9LM6srRyJsjM/tn4JuSvj82nUjwKO0EzoxtZ7KxPUDOnGNAK3pncw/gUEm3JI/X9wyRu0LpdXEP6lPl2HIQaX/gTcCbB76hDYDbJscJjGmf3DZNAbdPjhOYxYQMVW/cfhV4v6RNwF2EAksN4GpJZxF0w6+oNNJqA2dHSUpS+Yhy++N2ImfpzCUncUnbc9KVxgg5S46ueh/W19aRnCTvxQr5SVoDpF1KWzqSlkTu0pHApHKWuN/Sv8YhNUdC85BaI3NWrT54Zof+Kz9oZtuH7N8NHJW8PhK4d0Cf3ZIWgIOAhysem/J9wLHAF+Pf85HAFySdEH+wNxJum3qPcdtUbs+BbYKx7ZPbpukxMfs0E3K6gqoyubWW042qlzSsNpJL6FZFsXI0a1Sas5ndCuQM64mTnY7j1BNDtMbzc9wMbJN0LPAtQhDzq3r6FB7FzwKnAZ8yM5O0E7hC0tuBw4FtwE0D52r2ZZIgX0n3ANvN7MFx3kAdcdvkOGPbJ7dNU8LtkzPvjFsEdr2o5Q1d1UDnnEcWQIVXtZHxtDYST2wR7JzzyEIZDJ3xAttqvLMp0aOpLq9ru2sfgIqg5q5A6MR7W+xv9Hti03kVIyr1pBYV7JOh85XrPRC6gmd2+PFmy5LOBq4nRLxfama3SToPuMXMdhJSvr4vBjU/TLhIIfa7miDJWAbeYBa+NUkfAF5MkM3sBt5iZpeseqLOUNw2uW2qI+PYJ7dNc06xWgLTT9k9iF27yu3jjhvcNi6jVqqqvv/iM0tXk3xlaSCz+MnM4pwdZx0QrTFujgDM7Drgup62NyfbexggsTCz84HzM+1nVDjvMSudq+M4s8R49sltk+M408BXjhxnA2PAEs2R/RzHcdYat0+O49SRDR1zNFUGVT7v3Z9KMwqpxKCg51yAc5SsKK0b0pGzpG3pOFGmkmnLBj2vRMISZR6Wyj1aQ6Qr7eS9LKdBzzGYuZ2ZTyuRoXRO2y8vSWdtadrXKKtJa47MayC0mVgyv/iYK9w2lW1um2qN2ydnRdRBSpeSk81NSko3KVbzmaW1keZYdjfOO5d0KfAy4H4z+6HM/hcTMj7eHZs+aGbnjXFKoA43R44zA4RUubO4OOw4zkbH7ZPjOHVkArK69wIXAZcP6fO3Zvay8U7Tjd8cOU4FDLFk/u/iOE79cPvkOE4dGVdWZ2Y3SDpmMrOpzmxa00b/fWhXDZBC7tKVJSq2pXKWhX7pii2kcpfimGR/p61fKmJdUpoR7yEqP7oyNMWMUGqVbZ1aIqlcJX3/rUY8JpGSdOaTtBX1V9KaI5lpqZ3IVIrTJBmhChlL5SxRsGFkLK3RleadecdtU3mM26Y1xe3THFPIt6pKt+ogpZs1cp/ZIKld7vtY6Xe0gRjxk3OopFuS1xeb2cUrPMWPSfoiocbab5vZbSs8vo/5+5YcZxW4Z9ZxnLri9slxnDpSYeVoWIHqKnwB+B4ze0LSqcBfEeqtjUWtrGmnhsioivOdfZlA53Q7V3G+65gY1LyQ8dgCVhy/kHhnm/3jFG1d3tlRTryiUHxXgHMMhE48pB1Pbfpe0roiy9E7m3hvi/cqlcGAXUHRPVMc5DstPLW5yvWVA6GT+cyyl9Y1/fON2ya3TXXG7dOcU9fViGnUKhqHSc9n0Apc7vuo63e0BkxzTdvMvpNsXyfpXZIOHbew9Px+W46zAoJn1rNBOY5TP9w+OY5TRwQsTnN86TnAt83MJJ1AUPE9NO64fnPkOBXwVLmO49QVt0+O49SRcRMySPoA8GJCbNJu4C3E+y0zezdwGvCfFeQITwGnW64uxAqZzZsjZeQsXXKXZv/+jgwlkaYsZNoWh0tX2gtFfZF07ChdSeQzHRlLV6GOzFtpJ9KV+H1qub+WiJYTWchyrh5Kkk8/BmmndVO0FPerf44ikbh0TS4TKD2ngdAGY1Wgd+YEt03978tt09Rx++TUkjpI6VLWej65hA1zWPtoHMtkZmeM2H8RIdX3RJmPb8ZxxsRlK47j1BW3T47j1JFxV47Wi/Wb86jq871kgpWzgc7pdhoo3MxUlI9e2S6PbGY7Paa92O+dbS/0e2c7XtlRkWhpltsiEDrxxDaKdLipxzbdvxQ9p8l8tBTnmKYLLraXkrb4nPpJlaTYtcLDkQZm5wKh4xSUeqdHeWpnEA94nhPcNsWDk0PcNtUet0+OUzNyCRvmZLUoZRaLDMzft+Q4q8A9s47j1BW3T47j1JFpJ2SYFu5qcpwKGNC2xsBHFSSdLOlOSbsknZPZv1nSVXH/jWlVaEnnxvY7JZ2UtF8q6X5JX+kZ6w8lfVXSlyRdK2nrat+74zj1Zlz75LbJcZxpIMKNxqBHXan9ylG2hkhnZyPt2H9MVu6S1gBpdj1Dt3SlE+C8WB7T2pSRrixmgp6LQ9I431TZUqhBEt1IoRpRUvejsRQDoZPK9I2l/jomjSQwu6hS31WtXku90ymDmpMmSwKgFf88RlarL76HgbVEMvtnrL7IuNmgJDWBdwIvBXYDN0vaaWa3J93OAh4xs+MknQ68DXilpOOB04HnAYcDn5D0XDNrAe8lBCNe3nPKjwPnmtmypLcB5wJvXPUbcPpw2+S2qS6MY5/cNjmOM01qf6ORodKNm6R7JH1Z0q2SboltB0v6uKSvxednTHeqjrN+GLBkzYGPCpwA7DKzu8xsH3AlsKOnzw7gsrh9DXCiwtX0DuBKM9trZncDu+J4mNkNwMN98zX7mJkVV5OfA45c0RueEdw2Oc7Y9slt05Rw++TMO0VChkGPurKSuf1MT8XZc4BPmtkFcRn+HCbl/WlkPLI572vX/oynNvWWFoHLacX5wsOaSYcL0N7UjM9p0LO6nsMxsS0Zuvg9GlmZvss7G140Wsl73RTb0my4+xLPaBFwnbQ14vtupCl9O97i4Z9jzlOrRJaR9dQW+1pJW1eV+iKaOxMIPTMpdDUqVe6hxY9f5GIzuzh5fQTwzeT1buCFPWN0+kSv6mPAIbH9cz3HHrGCyf8H4KoV9J813Dbhtml+bROMaZ/cNk2XtbNP4/Loo+F5qysd15QNnt573hIy7CAUZoLgUfo0dfkHd5wJU3hmh/CgmW0fsn/E5efQPlWOzZ9UehOwDLy/Sv8NgtsmZ64Y0z65bVpb3D45c8NGT8hgwMckfV7S62Pbs83sPoD4/KzcgZJeL+kWSbcssXf8GTvOOmCItg1+VGA3cFTy+kjg3kF9JC0ABxFkKVWO7UPSmcDLgFdPomJ0TXHb5Mw9Y9ont03TYyL26YEHHlij6TrOZNnoCRleZGb3SnoW8HFJX616grh0fzHAgTq4zwiqq/7GSuuLDKgl0qnM3h/0bGnQc5SptBf720J7v0yltTlst3LSleTT7EhXUjlLTnHTJV2JQcipBCTKWJpLZcdUXlO0pzKVZqM/CLtZSETU/5kN+unM1hqJcpbsr1n6G5etUj+7v4FmIz2zo7gZ2CbpWOBbhCDmV/X02QmcCXwWOA34lJmZpJ3AFZLeTgh63gbcNOxkkk4meCN/2syeHGfiNcdtU8Rt03zaJhjbPrltmh4TsU/bt29fmz9Ql9OtDxtQSldQxBzNGpV+8c3s3vh8P3AtIeDy25IOA4jP909rko5TB8ZZOYoByGcD1wN3AFeb2W2SzpP08tjtEuAQSbuA3yRo0TGz24CrgduBjwJviNmgkPQBwgXL90vaLemsONZFwNMJP8i3Snr3ZD6FeuG2yXECq7VPbpumh9snx9mgK0eSngY0zOzxuP2vgfMoPUkXxOcPTXOiXTT6PY1dAbw5T2QR9Jx4MS1Tmb7wtELptW1tSryzmzLe2U30tWW9s+lfQsb12UmX2+WdjedYKtsaSSX59r4wQLNZDlSk0G12nS8GQmfioFNyP6PdgdBFWxlAaDF1bvoddL3XdpG+Nw2Enq0UupMosmhm1wHX9bS9OdneA7xiwLHnA+dn2s8Y0P+4sSY7A7htctvktikwrn1y2zR5ammfHGeNGXflSNKlBAnu/Wb2Q5n9Ai4ETgWeBF5rZl8Y45RAtTk/G7g2/rgsAFeY2Ucl3QxcHb1B32CA4XScjYAhltvj3Rw5E8dtk+Pg9qmmuH1y5p4JyOreS75eWsEpBDnvNkKWzT+hP9vmihk5ZzO7C3h+pv0h4MRxJ+A4s0J7JhNSblzcNjlOiduneuH2yXEC48jnzOwGSccM6bIDuDwmdvmcpK2SDiuSnqyW+sdJZarHjz6mP+i3kEV0BT0X0pWMXCVsZ4Ke43Zrczl0u5CzbCrbLH6y7VHSlZSi5Eai5ihqiLQTuUpjX3LuznlS2YzF52TouH8h8zEO+mQ7I6ZSkmI7aStqoHQFR7dK/U3nfSeSnLLjbNQXMYMl98w6KW6bwjhum9Ydt0+OMycUNZFmJIlDhZWjUTUiR5Gr03YEsMFvjhynBhSpch3HceqG2yfHcepKtjh6gdmoGpEjh8+NOsZ4gN8cOU4lDFgeXoHecRxnXXD75DhOLZFgy5bB+596atwzrKrW2ig2zs1R7s60K0PBNVtcAAAgAElEQVRRLiNUlKZk2qCUhWSzP21O64vEfUkZ4KJfeyHN1JTMLfc7VkhXEolHYznW+0gyQjUTiUx7b1HHpD8jVCrTyd1b5778rmlFSYpSmYr1y1Q6/ZIm65K2xFoiiUxlFuuLtP3iw1kNbpv63oPbpsnj9slx5oAZkdN1kKY9553A2ZKuJCRieGzceCPYSDdHjjNFzOSeWcdxaonbJ8dxakmjMXzl6PHHhx4e66W9mBCbtBt4C7AIYGbvJpQgOBXYRUjl/csTmPWM3hwVXr5BVeuzFdfpbyu8mKlHNlddPvG6FrVE2qmHtOOxTdo2Wxwj8ZqmVeobhSs2mXcn6DmpkbIcq8wvJ/PaW243F4qxM3VM0nji+FlZxoud/hEodZbmApzjduFxDefO/BnlgpUzVepnpb6IActtv/hwRuC2qXwLbpvWDLdPjjMlPvKR8HzKKes7j1lllKxuBIPqpSX7DXjDqk8wgNm8OXKcdcADnh3HqStunxzHqR3Tl9VNhdmbseOsA4bLVhzHqSdunxzHqSWjZHU1ZePcHOVqjaQyDfXLVKxTXyQNji4P6UhXkk+pkLGkdUMKyUprcymvKKQrtimVriSSjFjvoyvKuFCKJNIVWnEOiXSl0RWEXcyxX7rSFfSdiwnP/Zh2TdH62goJSTYQOhkvHduizEXKfM6zUl/E3DPrrBK3TeUwbpumg9snx5kOLqcbD185cpyNi2v6HcepK26fHMepJWPGHK0XG/PmaEjBqa6g32IzdQY2+z2Iqeez3QkyLo9pRw9sO/XO7he9iptK12ZjsXRFNqJ3Ng367cT3Jt7ZdvzBay+Vk2wtltvtvTGYOfFOl+lyk/eV/UyKIOOyRa00zW14kw3rb6OVuGwLr0CuWj1Jlfp0f1GlfkZS6HqRRWciuG3qmj+4bZoEbp8cZ4Pw6KPl9tat6zePSTGjK0fuanKcirSsMfBRBUknS7pT0i5J52T2b5Z0Vdx/o6Rjkn3nxvY7JZ2UtF8q6X5JX+kZ62BJH5f0tfj8jFW/ccdxas849sltk+M4U6GIORr0qCl+c+Q4FTCDVrsx8DEKSU3gncApwPHAGZKO7+l2FvCImR0HvAN4Wzz2eOB04HnAycC74ngA741tvZwDfNLMtgGfjK8dx9mAjGOf3DY5jjM1ClndjN0crftaVypTKGUczXzn8qCVnygj3cgGPaextc3uZyglK+3FRKayWEhXknlFyUpzy3LZtKncXlwI0o2FZiJniedOlRvLrXDypeVyEvs2lV9bazFsL3dVoS9qraTvWT3PZd0QJT+eSst4FNvpd1RIWxLpijKB0F3yk0Z8j83yPRSzsDTqeVgAdK6+CKxhAPTYspUTgF1mdhdArOa8A7g96bMDeGvcvga4SCFSfAdwpZntBe6WtCuO91kzuyH14vaM9eK4fRnwaeCN47yBecNtk9um8s3QT21sE4xpn9w2OU6d2bOn3K7xDUWWGZXVzd6MHWcdMBjlgT1U0i3J64vN7OLk9RHAN5PXu4EX9ozR6WNmy5IeAw6J7Z/rOfaIEVN+tpndF8e6T9KzRvR3HGdGGdM+uW1yHGc6eCrvKVFUO2+O8NimVPXYZRxtXYHCjcFtuRS6JB7b5ubgaty8eanTtn+yvXlhuesZYEH9Xud97fC+97XK9//k3jJX71ML4eT7FhY7ba1GnJwyP5aJd7Godt/tkU32t2Kl+HbiVY1e2UbXBxDT4SaV6dUqXa0Wj1fisS08410pdJWrTF+TAGgb+Wf1oJltH7I/59btHXFQnyrHOmuN26bw7LZp/RnPPrlt2igsl/+zs+itd8gnYZjBm4sOM7pyVDnmSFJT0j9I+nB8fWwMzPxaDNTcNGoMx5lVjLETMuwGjkpeHwncO6iPpAXgIODhisf28m1Jh8WxDgPurzLJWcRtkzPvjGmf3DZNEbdPzlwzozFHK0nI8OvAHcnrtwHviEGVjxACNh1ngxI0/YMeFbgZ2BZ/GDcRgph39vTZCZwZt08DPmUhx/BO4PSYMepYYBtw04jzpWOdCXyoyiRnFLdNzpwzln1y2zRd3D4580uxcjToUWmIkdk0XyvpAUm3xsfrxp12pZlJOhL4OeB84DdjIOZLgFfFLpcRgjX/pPKZO9XFc5XQk1X5orZHIougmQmEba8iELqYyqBbxFytkYychRhwrMVyDguLYXl7y6ZSrnLA5r2d7act7gNg/4V9nbZNjf6o3+V4oj3LpTTlu4vlmI8vbg7PzfLcexphPi3KY8qJp9KV7mco5SpQxip31ReJwddKg56Xh9QXgVLvkX6vRVuu/kiXdr6/5kg2AHoNgp/b7Uo3QVmiTv9s4HpCVP+lZnabpPOAW8xsJ3AJ8L4Y1Pww4SKF2O9qQoD0MvAGM2sBSPoAIbj5UEm7gbeY2SXABcDVks4CvgG8YtWTrzFum/r7um2aP9sEq7dPbpumx1Ts0zBmUL60ofnIR8rtU06Z7NizkqRhzJijJJvmSwkr1TdL2mlmt/d0vcrMzl79RLup+p/0/wG/Czw9vj4EeNTMCoFrlSBMx5lZilS5441h1wHX9bS9Odnew4ALBTM7n/AD29t+xoD+DwEnjjPfGcFtkzP3jGuf3DZNDbdPznxTyOpWT5VsmhNnpDWV9DLgfjP7fNqc6Zp1j0l6vaRbJN2yxN5cF8eZCcwGP5y1x22T45S4faoXk7RPDzzwwFTm6DhTZ3xZXS6bZs6h8G8lfUnSNZKOyuxfEVVm9iLg5ZJOBbYABxK8IVslLUQPyMAgzJgu9GKAA3Xw2prpnCxiFLksUWnpiuJ2MlVSNMPYjYVSUrG4GCQX+y2W2WMKuQrA1k1PAXDgYrk0ul8j7G+onGuhF3+qXcZsfmepvAvf3AzjLyRZlB6L0pUnk/dQyFi6Mz6p6xmgsVyee7mzP5GzdKQrqZwlnjupd5LLDtWpcQIofpCW9itquyS/FYoFVtY7M5Qh2mOuHDkTx22T26Zk/3zaJnD7VFMmZp+2b9++/n9kK+HsqG666KL1ncd6M2kpXUqdpXQpo2V1o8qgVHEo/G/gA2a2V9J/IshVX7Kq+UZGWlMzO9fMjjSzYwg640+Z2auBvyEEZsJ8BFU6c44NeThrj9smxylx+1Qv3D45TmT4ytGDZrY9eVzcc/TIjJhm9lAsRA3wv4AfGXvKYxz7RuBKSX8A/AMhYHNtKDx1g7yvOU9splJ68avRVcJj1C9JEWPbSDpGb6qStsJbmtYKOWCxlO5sXXwqPpc+1P2jd3axkdQqiOxtlwHMT18o78L3a4YA6E2ZYyxxKz8ZvbKtxBNbVJ9vLJdty8m24pCNpcR7u7k4pryvtqXosV0u/5y6vK6F1zYNTC8CpBupmzvzvXa8oUlAeC4AetqV6Q1sjIQMzpritqnT5rYp9N/AtgncPs0W62efJkWRaGDQysi8rxg5JeMXge1k0wS+RXA0vCrtIOmworA08HK6s0OuihXdHJnZp4FPx+27CIFSjjMXWLWU3c464LbJmXfcPtUXt0/O3DJmEdiK2TR/TdLLCRkzHwZeO+60Pe+j41TAGC+Vt+M4zrRw++Q4Ti0ZP1tdlWya5wLnjnWSHmp/c2Q5ycmQfgDKSVsK1UQSPNsZs1+F0tc+tC2SBisvxNoeRVAylDITgKctBBnLQQuldOXpjRAAvaVR9mvEiS9Z+VUd0Cz/0Iq+i5k6JOlUix/OpxLJSSFTUXm6LhlLMY1WKl2Jxzf2JQHMmwrpSkauArBQBEqXcywCoJUUZSnkLkpkKEUAdBH8HPqtg4reZStOD26b3DZBDWxTmIzbJ2ftmGaiAWdjMb6sbl2o/c2R49QGj2x2HKeuuH1yHKdujCmrWy9qNePC66ZGGoXc7O9XeGwHeekywbPqeHnLbsXxXeNkvLNdVdqL/SP03YodU49tGphcBDgXzwBbm0/GtjI4uhlP3kq8mE9rl/u3RNdqI5l4K85tX7v8epda4XNcXio/z6XodW3tK/u1yunQ2Bf2L6dt0WPb3JQGPcfUt4tJutzkPNYqvLfJd1kEQKeV64sA6FwAe7/zOZ6oOGbalenlntk5xm2T26ZwQB1tE7h9clbEo4+W21u3rt88nI2Prxw5zgbGPODZcZya4vbJcZy64itHjrOB8YsPx3Hqitsnx3Hqhq8cVaSQE3TkBankIFOTtp3Tj2SObffLHZTIIjpB0Wk9i6IMRTp0GhQdZQpdtUYyxxQ/SqnnrthOpSvNZHtRQYvxtESmUkhWDmzsSfr11whJjynGSVmyIBFJpSt7W2F7z1LZVshY2kvl595KgpmLGiKFhCXsj8+bk37xGC0mgdALyXYxflqFPla2tySYuVOFXrm25NhE2rSmAdCu6d/YuG0C3DbNpG0Ct0/zzHL8X6zqoZ9hKd39zWAbntUapGetKaNqQ21UZjTmKPOL7zhOHwa0NfhRAUknS7pT0i5J52T2b5Z0Vdx/o6Rjkn3nxvY7JZ00akxJJ0r6gqRbJX1G0nHjvH3HcWrMmPbJbZPjOFOhSOU96FFTZu92LqXLI5t6XeN2GvQaPbVqJd7Xoi0ZJ806W2ynDtDCK9uVVrdzuvJHqBiyPULq0KCc96Z4oi1J/tr9k9S5BVus9Ng2o7uwRXmePRYq1j/VKivXP75pMwBPbinb9u4N2/s2l8HIrc1Jaty9Ybu5qTx3e7F4Tt5rDIC2xMubDYBOg56LD7eRaUsCvCm8Q4kXd2AA9JRJ/8RWiqQm8E7gpcBu4GZJO83s9qTbWcAjZnacpNOBtwGvlHQ8oSr084DDgU9Iem48ZtCYfwLsMLM7JP0K8PtMoDCaUxG3TW6b1pjV2ie3TRuAYZ75PeVqb50vRqsycytGBfO2YlQwo7I6XzlynKqYBj9GcwKwy8zuMrN9wJXAjp4+O4DL4vY1wIkK2p0dwJVmttfM7gZ2xfGGjWnAgXH7IODeVb1nx3Fmg9XbJ7dNjuNMh0JWN+hRU+o7M8epE9YTy7FyjgC+mbzeDbxwUB8zW5b0GHBIbP9cz7FHxO1BY74OuE7SU8B3gB8da/aO49SX8eyT2ybHcaZDIaubMWbm5sgydUF6OvRvtzJyliTouVNLZDmVs/Rvq0vOEgOhk1hkxbZ2q/TQLcf6GUUND4CldrndjlKTdrJ4V8hY0kDmxdjWFTCdyF2IAdL7rBz7yWaQqTyxWP5BfmfTfgA8vq9se2JzkMUsbS7/DNpJMHN7k8XnJOg5ylgSVQzNKGNpJ0HNjSTo2RYyQc+xboiSNoufTyolsqK+SLqUngmAtnayCFroS5Lg6fHrioz0wB4q6Zbk9cVmdnH3AH30TmpQn0HtuZXfYszfAE41sxsl/Q7wdsJFiTNh3Da5bUoOTjbXyjbBmPbJbdNGYTn55y+88jN4YepMgEJOud7f/4wmZJi9GTvOejHcM/ugmW0fsn83cFTy+kj65SRFn92SFgiSk4dHHNvXLumZwPPN7MbYfhXw0aGzdxxntlm9fXLb5DjOdJjRmKP1uznqTZsLHa9a6mkrU6P2VzDv8timHr0imHmhP8CZnPc1DXpeTre7n6H0yjaWk5SuMa1smnZ2KQb4FmlqoTsIeU+MHt5nSaX4+B5bGWfcYuLISz21bfrT7j49emwPaJaBmAcuPBX6Le7fadtvU/iDfWqxjGpeWkw8tXG6WU9s0tZeKNqSdMGp1zWm3tRC8h0W+5PPsfO30OhvU/J3YuuVs3a8094MbJN0LPAtQhDzq3r67ATOBD4LnAZ8ysxM0k7gCklvJwQ9bwNuInhtc2M+Ahwk6blm9o+EoOg7xpr9POG2CXDbNFO2KZx8tbht2iikHvph6b1XkqRh167wfJwnFZwp6nJDMoGVI0knAxcSfmzfY2YX9OzfDFwO/AjwEPBKM7tnnHP6ypHjVMHK2jKrOjzo9M8Grif8g19qZrdJOg+4xcx2ApcA75O0i+CVPT0ee5ukq4HbgWXgDWbWAsiNGdv/I/CXktqEC5L/sOrJO45Tb8awT26bHMeZGmPGHI2TTXOMWfvNkeNUZkynsJldB1zX0/bmZHsP8IoBx54PnF9lzNh+LXDteDN2HGdmGMM+uW1yHGcqjC+r62S+BJBUZL5Mb452AG+N29cAF0mS2eoDOmfn5iitG9KMHrK0orxlttNA2SLgOD1mOWw3lpJA6OUkcDfKWFKZSlHaIyn30anWnkpXigrve5bLj/jJ5VLvUchYnmyXspE9cXspqa/RjlqZVE7e6NoOc0wDpbc0Qqn4/eMzwH7NMOH9F8qJb1kIYy8slscuLZZnaseq8pbWDSlqiSR/OYV0xRbUdywkAdCNZObFdjORsywXxVsy/dLvcp3qimgdFTNOjXHbVJ6va9tt01ri9snpYpiUaSUXqy6nc8ZhtKxuVDKrcbJpPrjaaY+8OZK0BbgB2Bz7X2Nmb4la4iuBg4EvAK+J9QwcZ+NRVKB3aoPbJseJuH2qHW6fHCeE8O5bbgzrMiqZ1TjZNFdNlSKwe4GXmNnzgRcAJ0v6UYKm7x1mto2gGz5rnIk4Tu2xIQ9nPXDb5DgFbp/qhtsnZ+4xC7lBBj0qsJJsmvRk01w1I1eOombvifhyMT4MeAllRpvLCHq/PxlnMlUpZIRdNUVSaUvRnu4vMkItl1oHRTmEknofjaX+7FBdbVH50VxK6mt05CxJLZEoXXlqXylX+e7S5s72d5ZDbY8nWuXy9uPNsP30dtlvS5SuNAZkhOq0JfsXo55jU1LwZEuceCFhAdjcDPsXF8rPZE8iXbHFWEskyaxVyFO6pSvFc5K1KZGXWJQaFc9hwrGWSJLpqagrYokHtNhvqewlU1ekzByWZBRL5UzFecaoKTJmEVhnwrhtctvktikZxu1TrVh3+zQsW50zn6RZCou/iyn/fbTb3addBavOpjnOSausHCGpKelW4H7g48A/AY+aWfELl1bF7j329ZJukXTLEntzXRxnNnDPbO1w2+Q4EbdPtWNS9umBBx5Ymwk7zhQYZ+Uo/q8UmS/vAK4usmlKennsdglwSMym+ZvAOePOudItY0zN+QJJWwlZZn4w123AsRcDFwMcqIP7+6Q3dx0PWr+nLfW+deqKJAHMXXVFovfOkuBhNWN9kVZ/0LOWy7ZGsr8Rg5gThybNqAxuJQrhxr4w78bexGO7NwY97y29s09sKQOcH9kXvLMHLpS1PYraH/sndUE2FcHMjfLWezFT7W8pqUJfVLZvJP06wdGNpMJ9M2wvNMu2RrP8HFtx21JPbDxN6mm1XFvOU5sNek7aiuDytL5Mgfo9trB2dUU0ZipvZzq4bUrO7LYptM2ZbQK3T3VlUvZp+/btK/9jyq0I+GpS/fnIR8rtU06Z7NhpIo6KmrZxmcDK0VjZNFfLiv5DzOxRSZ8GfhTYKmkh3tXlNICOs6Fw2Up9cdvkzDtun+qL2ydnXpnEzdF6MFJWJ+mZ0euBpP2AnyUsbf0NQdsHQev3oWlN0nFqgctWaoXbJsdJcPtUK9w+OU5gzIQM60KVlaPDgMtildoGQe/3YUm3A1dK+gPgHwiav7WhCHBOg2jbmQDXdqaWSBowGwOgtZQEQi8lAdD7ilojSV2MGNiclOnoyFnaiXSlvTccs7S3/Igf31MGM++/GKQrDy8+rdNWBCRvSQqVFMHMrSRTYbq/YCnRl+yLWpL2iHvfRgyeTktzKA2oLmKHk/1WtCUlQDplWpppv1TaEoOZk++rCHDukqkUcpauWiJx/xrXDOnD3DNbQ9w2uW3qMLe2Cdw+1ZP62adJyel27QrPXgNp8kxaSjeINZJWzurKUZVsdV8CfjjTfhehcq3jzAfuga0VbpscJ8HtU61w++Q4ZSrvWWNmovIsSYdbBEBb4jZMU+d20ukmx1g7E/TczgQ9JxXpm3G7vTdp2xTO2UySW7X3xBSySbX21qbgYWwtli7LPYtl0PNjiyEwbnPzgE7bpkZMX6t+V2Qa1JwGRTeju7BlpUfzyZhu97tJ2t3i+FbymbXjduqR7fLOxm1LnbyFdzZp63hsU49s5piuoOciDW4S9NwJZu5yF2cq2KeB8sV2znubenmLQPpOYH2m/wi8Ar2Tw22T26Zy8PWxTeD2yZkgjz5abm/d2r9/zleM7m8Gm/WsVh2WjeuN2QZdOXIcJ+IXH47j1BW3T47j1AxfOXKcjYxr+h3HqStunxzHqSEbNuZoTSlkCGkgbEdykAngTSvPKxP0nFYhL5Y/EwmEYtCzNdOg53Q71uTYl0hX9oW2hVS6sql4Lucdy4LQ3lRKTpYWk7oii0FWspicu5CupFXm2zHYeY+Vxz49rSsSK82nAc572qHvk+1SKlNs70vKxxfSlVQClCWZj3Wqwqf7e57plbHEY9L9KmQzXRHX/edu9PfLOkiTv4+OtKk9OVeqGP/iQ9LJwIWEYjjvMbMLevZvBi4HfgR4CHilmd0T950LnEUQ6fyamV0/bEyFD+wPCLn/W8CfmNkfjfcO5hi3TeFYt00lNbFNML59cts0x6RXrkUdnJyUbpTUbo5wOV11XFbnOBudMa5nYsaidwIvJVRFv1nSTjO7Pel2FvCImR0n6XTgbcArJR0PnA48Dzgc+ISk58ZjBo35WuAo4AfMrC3pWaufveM4tWe1sUpumxzHmRIuq1tDCq+b0vSsaRX6uG1JgLOiW82aSVv0zrKQ9842oie2uVh6/ooA6PZC2dZcLIKey/ksxCrsiVOV5YUkAHoheEsfayTzpp+lmIN2z6bEs9ssb8M3NzKpc6MHNvXOfnc5eIOfapXj7G2Ffsut5P2N9NQObus6NOOJ7QpmLk7ZlS63qie2vyL91KvRjy9bOQHYFTMVIelKYAeQXoDsAN4at68BLope1h3AlWa2F7hb0i7KbEeDxvzPwKvMwhKFmd0/1uydSrhtcttUbq6RbYqTGcM+uW2aZ4rVolHM+WpRbZixNOqzKqsbWQTWcZzIeEUWjwC+mbzeHduyfWL19MeAQ4YcO2zM7yN4dm+R9BFJ2yrN0nGc2WT19sltk+M4U6FYOdqIRWAdx2GkZ/ZQSbckry82s4vTwzPH9F62DOozqD3n3CjG3AzsMbPtkn4RuBT4yezMHceZecawT26bHMeZCtOMOZJ0MHAVcAxwD/DvzOyRTL8W8OX48htm9vJRY9fz5iitFZELYC4CXNO2dlqnoqduRNo3DaSLAdAdCQtgjfJWVlGeon2l5KS5EKUrSd2Qhb1RKpN8mu0oXSmew9jlHJebofNTqXQlU6iiCFJOJScHLpZ/afs19oW3kgZKRw3JU4l05TtLYen8iaWyvsje5ShdST47a3dpTsK8utr6pliSqVafttuoAOeiLdcvlb2sRyzkaA/sg2a2fcj+3QSdfcGRwL0D+uyWtAAcBDw84thB7buBv4zb1wJ/NnT2TjXcNnVw21QT2wTj2ie3TY4zK8yInK5gyjFH5wCfNLMLJJ0TX78x0+8pM3vBSgZ2WZ3jVETtwY8K3Axsk3SspE2EIOadPX12AmfG7dOAT1kIWNkJnC5ps6RjgW3ATSPG/CvgJXH7p4F/XM17dhxnNhjDPrltchxnKhQxR4MeY7IDuCxuXwb8wtgjRuq5cuQ4NWScCvRmtizpbOB6QmrbS83sNknnAbeY2U7gEuB9Maj5YcIFBbHf1YRg5mXgDWbWAsiNGU95AfB+Sb8BPAG8bvWzdxyn7qzWPrltchxnWkx55ejZZnZfOI/dNyTz5ZYoK14GLjCzvxo18EzfHKW1InLZoZRIYIrsUErcaEV2qFS6okReQswO1di3nBwTJBTNVJIS29rNRM6yYF37+rbjeVqUkpQn43M7kYoUspIiexN0y0+2LISMUAsZ9+C+dvmhPLkcZCzfXSrlLE8thTGXlsp+tpzU5GjFeaSqoUIVlPshXsmPc06mUmeMrs9hVUOYXQdc19P25mR7D6H2R+7Y84Hzq4wZ2x8Ffm68GTurxW2T26Y1ZUz75LZpjkmvXBdm+pLQqSEVYo6GxmtL+gTwnMxxb1rBNI42s3slfS/wKUlfNrN/GnaA/yc4TgXEeCtHjuM408Ltk+M4daRCKu+h8dpm9rOD9kn6tqTD4qrRYUC2LICZ3Ruf75L0aeCHgRm/OapamT4XAJ2rTJ8GPS9nPITN1DsbPCpK6oY04v5m0lZ4XRdSD3Hslwb/dgcCF+cuGwtP7VOJd7YV63zsXSq/qic3lx7dTbGK/aZGfyTwcnLCIsC5eAbYsy+Ms5yMbUuJd3a5CHpOpl1sp17KzI/yWDWBUg+5MmFx6f7i+xwVFJ0LlF8hfvHhdOG2CXDbNHD/GtomcPvkrJJRq0XFlW3VekiOkzBlWV0RC3lBfP5QbwdJzwCeNLO9kg4FXgT8j1ED1//myHHqwpiyOsdxnKnh9slxnJoxzVTehJuiqyWdBXyDKP2VtB34T2b2OuAHgT9V0K03CDFHtw8asMBvjhynCuaeWcdxaorbJ8dxasg0V47M7CHgxEz7LcREL2b298C/WOnYI2+OJB0FXE4IiGoTgqUurFp8aa3oCoBuBBeaWSlnKAKgLQmEVpQ9WKPUnGQDoPeW32wjtlkicWlG2UR3UHM8X5d0ZUQtjdi51S6/lr0xCHl5uZxjKmNZXAjzXWiW8y5Ok3wkLLfC8UvJOPv2hXFaexPNzVI5r8ZSfE7+sDtBz7lA6EGey6o/2qsJgC6+o9b0C4yMJcdxJo7bJtw2uW3q4PapXsyKffKEDM40qRBzVEuq1DlaBn7LzH4Q+FHgDZKOpyy+tA34ZHztOBsXG/Jw1gO3TY5T4Papbrh9chzC/fegR10Z6SaIOcSLPOKPS7oDOIJQfOnFsdtlwKfJV6adDFUr06ck7kmLkbBd/r9inHZShX458ejG/coERTcST2zhDs0FOKeV56375D3PZYdOmlqgFefTSrymezYnHtbF8Bk0GuVnocIznHxkFgPB28nY7b4Ncw8AABpsSURBVCJN7r4kqHtvsh3P2UjOrfjHrMQZWmynsg6l3nLL/Drn2obR9f2OcJHGvko+k9R7vyrMPbN1w22T2ya3TcUgbp/qRm3s0yhGrRZ5Iob5Y4KribO6crSidy3pGEIKvBupXnzJcWYe4RcfdcZtkzPPuH2qN26fnHllw98cSToA+Evgv5jZd1RRgy3p9cDrAbaw/2rm6Dj1wOUptcRtk+Pg9qmmTMI+HX300dOboONMkSmn8p4alW6OJC0S/rnfb2YfjM1Viy9dDFwMcKAOXhPzXcgUUukCRJlGO5EzxEBZDQpGLrZz9UUSSUoj9msm43QCoLvsYL+8RmlgdhE83ErbouRkKZWcJPVHYrX7VjP5aBtF/ZXkRMXbTseOAdXalwQ67022Y7uSP+wiALqxnEhTOvNO5pCpNZKVsKSsVM6ylli3HMepB26b3DZ12ubVNoHbp5oyKfu0ffv2tf1yVyOrevTR/ratWyczH2dtmWBijimn8p4aIxMyKLg5LgHuMLO3J7uK4kswoPiS42wkZIMfztrjtslxStw+1Qu3T45TrhxtuIQMhGqyrwG+LOnW2PZ7DCi+tCYMq0wPQwOg08DkTgrdVhIwnAi3LW5r31KyP1O5Po7ZzLR1TTuznN71w9UJek6GKYKeN6fe2fKgdvTOWvJNWs47WzS1U+9s9zmg2zvb3Buf9yXzyaTQLTy13Sl0+4Oeu7y3xXeY88i26ymed01/7XDb5Lap3D/HtgncPtWQ+tmnqqxm5cBXiabHrl3l9nHHrcsU2pWSW2eO26gxR2b2GXoEGAl9xZccZ8PiHtha4bbJcRLcPtUKt0+OE6jzCtEgVncr6DjzRkyVO+hRBUknS7pT0i5JfbUtJG2WdFXcf2PMcFTsOze23ynppBWM+ceSnljNW3YcZ0YY0z65bXIcZxoUK0eDHnVltssh5+qLQEfGUtTPgCQAup1qVzK/Ghl5SVo9XkXQcy44OmkrDhn0AStqaJRoacrg4aRflJWkUpFWEqRsi2G7kLDAoIDrSCoviefpkq7kZCpJW3Ofde0LxxdtiVwllam0MjKVYh6WkbOkFJIkW1/NSEiVu3rXrKQm8E7gpcBu4GZJO83s9qTbWcAjZnacpNOBtwGvjIUDTweeBxwOfELSc+MxA8eUtB1wrcN64LbJbdMaMo59ctvkODVnnaR0KY1RNdwGMKuyOl85cpyKjBnwfAKwy8zuMrN9wJWEYoApOwhFAQGuAU6MQb07gCvNbK+Z3Q3siuMNHDNe8Pwh8LvjvGfHcWaDMeyT2ybHcabGRk3IsCHIpdC1Isg4k0IXEudm4r0sfmdS72xnO/Hi5u46sx92V+X26LFNUtoWns+0Cn1jU3lMezE+LyRzLE6e8852nS+O1+UNLreLYOfGvvKgwlPbTDyxhec49cg2WhlPbdpWfOaTqFC/Fli31zzDoZJuSV5fHFOxFhwBfDN5vRt4Yc8YnT5mtizpMeCQ2P65nmOPiNuDxjwb2BnTxQ6duLO+uG3KnS+O57apGuPZJ7dNjuOUdysTTOU9zZUjSa8A3gr8IHCCmd0yoN/JwIWEuhnvMbMLRo09NzdHjjM2w6+LHjSz7UP2j7gcHNpnUHvuOtckHU7IgPTiIfNxHGcjsXr75LbJcZypMOUisF8BfhH400EdKsqG+/CbI8epgo0Xc0T4pzwqeX0kcO+APrslLQAHAQ+PODbX/sPAccCu6JndX9IuM1t/4bLjOJNnPPvktslxnKkwzSKwZnYH9MTZ9tOR+Ma+hcR3Tm6OcgHQq6kvkkhXOjKVXJX6TNBz7usZFNS1kKncXsgiuqQrrf6g53YScNyKQc/WLNuK7fR9dSaXka50BVkn24U8pSvouWhLpSudtqQmS1qlPtZqSSVCxffV9YNe9cc9HWcN646MWUzxZmCbpGOBbxGCmF/V06coDvhZ4DTgU2ZmknYCV0h6OyHoeRtwE+Fb7RvTzG4DntOZt/SEX3ysI26bynm7bZoaY9gnt03O+KRXwFu2rOzYXB2fGtT2mTsmKKcrqLByNCokYVyqyIb72Dg3R44zRUI2qNUfH3X6ZwPXE3Svl5rZbZLOA24xs52Eaurvk7SL4JU9PR57m6SrCZ6OZeANZtYCyI25+lk6jjOLjGOf3DY5jjM9DLOlYR2GhiRI+gSJQyXhTWb2oQoTqCIb7sNvjhynCmbjyuows+uA63ra3pxs72FAtXQzOx84v8qYmT4HrGa+juPMCGPaJ7dNjuNMBwP2jew18Giznx1zAlVkw31szJujQsayivoiXTKWOI61EkkGYX0w/RnK3ZaOkrF0JBCJ5qSUkpQ9O9KVJCNUezGRtiz2DZOVrlhmQsUcBklXCrlMI5GhFNKVZpolKkpWGhm5Std2V32RTEao4vPO1RcZ88ZkItRgCs6M47YpPLttmjw1mYYzp6xUSpeSk82NktJ95CPl9imnrP7czhqwrrXgqsiG+/A6R45TBQsxGIMejuM464bbJ8dxakmbsHI06LF6JP0bSbuBHwP+WtL1sf1wSddBkA0TygdcD9wBXF1F4rsxV44qkqsvMqpKfScQmuVMWz9dbe1+r2Mz8UQWsgi1S1drO3pqG4lHtrUp8d5Gb2mn8jzQ7nhn08DszOQK72w6h9Q726lSnwtw7m/TcjtpSwOgM0HPxXbaVgScryYQOqXjiZ/wRYFfYzhrhNumuMtt0wrGnvyQjrNujErI4KtFa0p7rLWU4UXYVouZXQtcm2m/Fzg1eT1S4tvLXN8cOc5KGDfmyHEcZ1q4fXIcp34YMDQhQy3xmyPHqciYqbwdx3Gmhtsnx3HqhzGtlaNpsrFvjirWF0klDqmMxYraHrmhk21FiYslcpasjCUbzFs2NVr99TW0HOaoTUlwdCIlKSQraSB0IVmxRvr+475UzdIJvE7a0o+n3S9TKfTrzX2JNKUT9JwcnG4Xmvfl8h+kEwidkfN0fUedXZl+KVP2msrcM+tMELdNycnjPrdNq8btk7PupMVshtXLGVUP6dFHw7PXNqoVjVUnVRgvW916sbFvjhxngnhgs+M4dcXtk+M49cNldfVmSArdbIX6ZL8lLk1lqp6XQc8LSdtybMvMIW3vCnpeiKdNPJ+LwSub/vA1lpL0vgvRO7sv8c4uRK9zevJhsXSW3+54Z5f7vbNdQc1Flfml1Pua2R4R9NzxwHZ5Yod7KyznqZ0Ghgc8O9PBbdNg3DZVw+2Ts96kq0XFKlJuBWlUyu+tWyc3J6cGzKasbmT6CUmXSrpf0leStoMlfVzS1+LzM6Y7TcdZb0KRxUEPZ31w++Q44PapfrhtchwoZXWTT+U9Tark5nsvcHJP2znAJ81sG/DJ+NpxNjZmgx/OevFe3D45jtun+vFe3DY5c0+xcjToUU9GyurM7AZJx/Q07wBeHLcvAz4NvHGC81pbEnlEd5X6qP1IPG9FZfcuCUuUw1gSkKjYsSs4uuucMcA5+eEqZBhpBfciYLiZBBEX0pR0u9FM26LkRv2B0F0UTalcJVPbJK0e35GuZAKcu+QqiYylCHbufl9FZfqkX3G+XC2RnOxlJYyQwIw+3jX9dWTD2ye3Tcm83TYNHsPtU93Y8LbJcSoxXzFHzzaz+wDM7D5JzxrUUdLrgdcDbGH/VZ7OcWqAX3vMCpXsk9smZ0Ph9mkWWNW109FHH71G03OcaVDfFaJBTD0hg5ldDFwMcKAOXn/znUuh27W/31Obq1JviSCx45VMPa3FvrRafSbla+rlVSt+HQsZT2SrTJfLcpI6twhwTqrQ0+hv63hqc3l8UzKe2i6PZLvd35bxznalxi22k7bCK9v1mbSLz3FEIPQ6kQt4d2YXt01umzaKbQK3TxuN1D5t3769Hn9kVRmWytuZM2YzlXeVmKMc35Z0GEB8vn9yU3Kc+iEz1Br8qDSGdLKkOyXtktSnNZe0WdJVcf+NqSRD0rmx/U5JJ40aU9L7Y/tXYmDw4lgfwGzh9smZK8a1T26b1gy3Tc6cUcjqBj3qyWpvjnYCZ8btM4EPTWY6jlNjxgh4ltQE3gmcAhwPnCHp+J5uZwGPmNlxwDuAt8VjjwdOB55HCPB9l6TmiDHfD/wA8C+A/YDXjfPWZwy3T878sUr75LZpTXHb5MwZ00vIIOkVkm6T1Ja0fUi/eyR9WdKtkm6pMvbItU9JHyAEEB4qaTfwFuAC4GpJZwHfAF5R5WS1I1dfpGt/DB7uCoQuamAkQca5QOhiX7KdrWbfSGQoGZlGZ8w0YLiZSkBirZFGMp8oWelqK97jqNvhNJY7Jxtp9UtuOnNLZTipTCVTN6Qz5sig5xEV6Yu+XZKjKSgQjK7g71VwArDLzO4CkHQlITj39qTPDuCtcfsa4CJJiu1Xmtle4G5Ju+J4DBrTzK4rBpV0E3DkOJOvKxvWPrlt6sdt02DGs09um6bAhrVNjrMippqQ4SvALwJ/WqHvz5jZg1UHrpKt7owBu06sehLH2QhouAf20B6PxMVRM15wBPDN5PVu4IU9Y3T6mNmypMeAQ2L753qOPSJuDx0zSlZeA/z6sMnPKm6fHCcwhn1y2zQF3DY5DjDFIrBmdgeABjkRx8Cj5hynEjYqTe+DZjZwWZcBzvmKfQa153ztvWO+C7jBzP52yNwcx5lpxrJPbpscx5kSIxMyjHIsT2oSH5NkwJ9WGd9vjmA6WaJGyVgKeUUzrSXS6JuPWrGtmchCFpLsUFE2Ykktkc5ddK5tJXfYOelKpraHMtKVLqlNR35Teg+s2N+VJSu+l5wHdL0zMRnjZqXaDRyVvD4SuHdAn92SFoCDgIdHHDtwTElvAZ4J/N/jTNxZR9w25XHb1M149sltk7N2JDXXKme127On3N6yZbLzmQa7doXn445b33nUhqErR0Mdy5I+ATwns+tNZlY1Zu9FZnZvTJ3/cUlfNbMbhh3gN0eOU5ExiyzeDGyTdCzwLUIQ86t6+hTBup8FTgM+ZWYmaSdwhaS3A4cD24CbCF7b7JiSXgecBJxoNokqk47j1Jkx7JPbJsdxpkSbcWKOzOxnx52Bmd0bn++XdC0hLtJvjlZERU9tNhA6uTk2RS/uoNPEsbsr0/d7LIug6C49eRoo3Iz7m6lruAhwTgOzV+Gd7RycBmH3e2fLAObUO9vq229dHttWf78hY2cDnYGh1eUn+btrdL+/lR4edPpnA9cDTeBSM7tN0nnALWa2E7gEeF8Man6YcEFB7Hc1IUB6GXiDmbUAcmPGU74b+Drw2eiZ/6CZnbfqN+CsP26bMpN12xTGY9X2yW2Ts6bkVotGrSbNwmpRSm7FaG5Xk0bK6qaKpKcBDTN7PG7/a2CkvfGbI8epRLWU3UNHCFmarutpe3OyvYcB2YvM7Hzg/Cpjxnb/33acuWE8++S2yXGc6TGdhAyS/g3wxwSJ7l9LutXMTpJ0OPAeMzsVeDZwbXTELABXmNlHR43tRspxqlKH2ALHcZwcbp8cx6kd00vlbWbXAtdm2u8FTo3bdwHPX+nYfnM0jGG1RrKB0Orbb+kNcyoBiWN2+fpi0LNSaUYMirZ2MnZSf6T4QewER0MpWUnm3Ql6TudYVcaSrd1h/fvToOac1CbXlrxXy/Wr+IM/tfohnRPQLZlxnPXEbVPfvOfWNoHbJ2e2qZqYIWXWkjTUQU63mmQYYzO9VN7TxG+OHKcSVsYiOI7j1Aq3T47j1JH1jTlaLX5zVIVVBUKrax+Qr1yfCw5OzxHbuoKaUw9hPE/XMUUgtEZ4YpOg6KHkPKRpYHaxPSgYuZMaN/0sMsHMnXESz/eosVfKanX57pl16ojbpv62ebNN4PbJ2XhstCQNdWDNVotSfOXIcTY2rul3HKeuuH1yHKd2TC/maJr4zZHjVMGsO7Wv4zhOXXD75DhOLXFZ3XywwkDolFGV64sfN2WkK5Z6BRPJSUfSokb/MV0nV3+/YqmzMSL4OSfXSOUj2arxmeryucr0mQDnbBX6AWSDnadVV3DMVN6OM1XcNsVzzKFtArdPjjOKWUvisGGYvVVtvzlynErYWEVgHcdxpofbJ8dx6kgbXzmaJyoGQqfe0NSTmKtcX/RNPbYqvJhpgHKaljazP5sat+ccXYzyzqZU9YamaXAzwcxDA5xz40zT41oFA1vvOThOFdw2lcyDbQK3T87GY1LJA9LVonlnXVJ5gydkcJyNjHtmHcepK26fHMepHZ6QwXE2LmaeDcpxnHri9slxnFoyh6m8JZ0MXAg0gfeY2QUTmdWsUTEQOidjyVWu76o5UoQuZyrYZ+dAUlfEMoHXyshDxvy7zUtORgVFDwlwHlE3ZFWBzhMIVjbPBjVTuH3CbdOc2CZw+zRLuG2aADmJWK7NEy+UrFudo9mLOapYaa8fSU3gncApwPHAGZKOn9TEHKdWWAx4HvSogKSTJd0paZekczL7N0u6Ku6/UdIxyb5zY/udkk4aNaakY+MYX4tjbhrr/c8Ybp+cuWJM++S2ae1w2+TMF8XN0aDH6pH0h5K+KulLkq6VtHVAv//T3v2GyFHfcRx/f7hTCxabGv9Uc1EjHm2j2FqKtLQUSXwQbTA+UEj/QGgFERQULFXrs9I+yJOmQrUQ1BJEiOEUcohtTDQP+sTU0zwoaao90mrOvycmVhSTXu7bB/O73Ho3e7e7czc7s/t5wXA7c7/5zfx2dj8zv9ndmQXzLU/HnSPgWmA8Io5ExElgJ7CpQH31FzE75P5/et4Q0zFvyCvHdJwe4tT07BCRneGcnp4dTp1qOsQyDLnLbXwuGtctDbPrnTM0tvv0U9fw/OQ9p2XI2y4tLr/FHeJtwLGIuALYBmxN864FNgNXAhuARyQNLFLnVmBbRAwDx1Ld/cT51MjZ1NvZ1GQbtrIOzqbSOZuWwuDg7DA1lQ1505r57LP5QxnKXl7XBdkV65oNhewFroqIq4HXgQfmFuj0ZESRztEq4GjD+ESaZtZzImLhg7DFtbJD3ATsSI9HgPXKvqe0CdgZESci4t/AeKovt840z7pUB6nOmztufD05n6xvFMwnZ1O5nE3WR5bvk6OIeD4iZnrALwFDOcU6OhlRpHOUd43VeaclJd0uaUzS2P84UWBxZt2VdyY994xxvlZ2iKfLpDf8R8DKBeZtNn0lcLwhNPpx57toPjmbrJcUyCdnU7naPnaanJwsYbXMlsPMBRmaDUvm58Cfc6Z3dDKiyK+zJoDVDeNDwNtzC0XEdmA7gKTJfTHyCfBBgeVWyXk0a0uR39h253e1zdtSP6205dJ2KvyYY3v2Te86b4EiX5A01jC+Pb32Z7SyQ2xWptn0vJMbC5XvJ4vmU042vUHvvA+cTdW05NkEhfPJ2VSujo6dNDDQH8dO9dNvbWkznz7aA7s7PnaStA/4Ss58D0bE7lTmQWAKeDKnXEeZU6Rz9DIwLGkN8BbZ945/vNAMEXG+pLGI+HaB5VaG21JNy9GWiNhQsIpWdogzZSYkDQJfAj5cZN686R8AKyQNpjO0uTvfHtdWPkXE+dA774NeaQe4La0omE/OpnL52MltqaQqHjtFxPUL/V/SFmAjsD5yLzHa2smIuTr+Wl0KtruAPcBhYFdEHOq0PrMed3qHmK7OtBkYnVNmFNiSHt8CvJje7KPA5nTFqDXAMPC3ZnWmefanOkh17l7GtlWO88msZc6mEjmbzJaGskvi3wfcFBGfNinWSr7NU+ii5xHxHPBckTrM+kFETEma2SEOAI9HxCFJvwbGImIUeAx4QtI42VnZzWneQ5J2Af8g++j4zog4BZBXZ1rkfcBOSb8BDqa6+4rzyWxxzqbyOZvMlsQfgLOAvdm1XngpIu6QdDHZ/cNubJZvi1Ws/E+hlo+k2+f8FqO23JZq6qW2WLl65bXTK+0At8UMeuu147ZUUy+1pajSO0dmZmZmZmZVVORS3mZmZmZmZj2j1M6RpA2SXpM0Lun+MpddhKTVkvZLOizpkKS70/RzJe2V9K/098vdXtdWpbuYH5T0bBpfI+lAastT6YdrlSdphaQRSf9M2+e7dd4u1h11zSZwPlWZ88mKcjZVi7OpP5TWOZI0ADwM3ACsBX4kaW1Zyy9oCrg3Ir4OfAe4M637/cALETEMvJDG6+JusivlzNgKbEttOQbc1pW1at9DwF8i4mvAN8jaVOftYiWreTaB86nKnE/WMWdTJTmb+kCZnxxdC4xHxJGIOAnsBDaVuPyORcQ7EfFqevwx2YtoFdn670jFdgA3d2cN2yNpCPgh8GgaF7AOGElFatEWSecAPyBd7SgiTkbEcWq6XaxraptN4HyqKueTLQFnU4U4m/pHmZ2jVcDRhvGJNK1WJF0GXAMcAC6MiHcgCwHggu6tWVt+D/wSmE7jK4Hj6f4LUJ9tczkwCfwpfcz9qKSzqe92se7oiWwC51PFOJ+sKGdTtTib+kSZnSPlTKvVpfIkfRF4GrgnIv7b7fXphKSNwPsR8Urj5Jyiddg2g8C3gD9GxDXAJ/Txx8DWsbq+/j/H+VQ5zicrqq6v/c9xNlWOs2kRZXaOJoDVDeNDwNslLr8QSWeQvbmfjIhn0uT3JF2U/n8R8H631q8N3wNukvQfso/o15GdDVkhaeamwHXZNhPAREQcSOMjZG/4Om4X655aZxM4nyrK+WRFOZuqw9nUR8rsHL0MDKcre5xJdoft0RKX37H0vdLHgMMR8buGf40CW9LjLcDustetXRHxQEQMRcRlZNvgxYj4CbAfuCUVq0tb3gWOSvpqmrSe7E7ttdsu1lW1zSZwPlWV88mWgLOpIpxN/aXUm8BKupGspz0APB4Rvy1t4QVI+j7wV+DvzH7X9Fdk353dBVwCvAncGhEfdmUlOyDpOuAXEbFR0uVkZ0POBQ4CP42IE91cv1ZI+ibZjyPPBI4APyPr9Nd2u1j56ppN4HyqMueTFeVsqh5nU+8rtXNkZmZmZmZWVaXeBNbMzMzMzKyq3DkyMzMzMzPDnSMzMzMzMzPAnSMzMzMzMzPAnSMzMzMzMzPAnSMzMzMzMzPAnSMzMzMzMzPAnSMzMzMzMzMA/g9PJOeHT9LWzwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x216 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(15, 3))\n",
    "\n",
    "a0 = axes[0].imshow(imgGalSim.array, vmin=0, vmax=1.5e-3)\n",
    "fig.colorbar(a0, ax=axes[0])\n",
    "axes[0].set_title(\"GalSim convolution\")\n",
    "\n",
    "a1 = axes[1].imshow(imgAnalytic.array, vmin=0, vmax=1.5e-3)\n",
    "fig.colorbar(a1, ax=axes[1])\n",
    "axes[1].set_title(\"Analytic convolution\")\n",
    "\n",
    "a2 = axes[2].imshow(\n",
    "    imgGalSim.array - imgAnalytic.array, \n",
    "    vmin=-2e-10, \n",
    "    vmax=2e-10, \n",
    "    cmap='seismic'\n",
    ")\n",
    "fig.colorbar(a2, ax=axes[2])\n",
    "axes[2].set_title(\"GalSim - Analytic\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:03:03.607277Z",
     "start_time": "2019-01-24T21:03:03.603609Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.197057e-08\n",
      "1.1641532182693481e-10\n"
     ]
    }
   ],
   "source": [
    "print(np.max(np.abs(imgGalSim.array - imgAnalytic.array))/np.max(imgAnalytic.array))\n",
    "print(np.max(np.abs(imgGalSim.array - imgAnalytic.array))/convGalSim.flux)\n",
    "# 1e-7 residuals compared to peak, or\n",
    "# 1e-10 compared to integrated flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:03:03.611928Z",
     "start_time": "2019-01-24T21:03:03.608622Z"
    }
   },
   "outputs": [],
   "source": [
    "# Any different if we use \"g\" shape instead of \"e\" shape?\n",
    "\n",
    "def momToSizeG(mom):\n",
    "    T = mom['Ixx'] + mom['Iyy']\n",
    "    size = np.sqrt(np.sqrt(mom['Ixx']*mom['Iyy']-mom['Ixy']**2))\n",
    "    D = T + 2*size**2\n",
    "    g1 = (mom['Ixx'] - mom['Iyy'])/D\n",
    "    g2 = 2*mom['Ixy']/D\n",
    "    return size, g1, g2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:03:03.616349Z",
     "start_time": "2019-01-24T21:03:03.613251Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8211602868378718 -0.13667504192892002 -0.27335008385784004\n",
      "1.0 0.0 0.0\n",
      "2.114742526881128 -0.10557280900008412 -0.21114561800016823\n"
     ]
    }
   ],
   "source": [
    "galSize, galG1, galG2 = momToSizeG(galMom)\n",
    "psfSize, psfG1, psfG2 = momToSizeG(psfMom)\n",
    "convSize, convG1, convG2 = momToSizeG(convMom)\n",
    "\n",
    "print(galSize, galG1, galG2)\n",
    "print(psfSize, psfG1, psfG2)\n",
    "print(convSize, convG1, convG2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:03:03.620879Z",
     "start_time": "2019-01-24T21:03:03.617579Z"
    }
   },
   "outputs": [],
   "source": [
    "convGalSim2 = galsim.Convolve(\n",
    "    galsim.Gaussian(sigma=galSize, gsparams=gsp).shear(g1=galG1, g2=galG2),\n",
    "    galsim.Gaussian(sigma=psfSize, gsparams=gsp).shear(g1=psfG1, g2=psfG2),\n",
    "    gsparams=gsp\n",
    ")\n",
    "\n",
    "convAnalytic2 = galsim.Gaussian(\n",
    "    sigma=convSize,\n",
    "    gsparams=gsp\n",
    ").shear(g1=convG1, g2=convG2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:03:08.295867Z",
     "start_time": "2019-01-24T21:03:03.622149Z"
    }
   },
   "outputs": [],
   "source": [
    "imgGalSim2 = convGalSim2.drawImage(**drawKwargs)\n",
    "imgAnalytic2 = convAnalytic2.drawImage(**drawKwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:03:08.522691Z",
     "start_time": "2019-01-24T21:03:08.297350Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0cAAADSCAYAAACB+bvuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXu8bFV15/v9Ve19zgERjoAPnoGEYxJMt6ZzLiYxDyOxecR40mlsQdtgGtvbHUnSeQoxH7Xp0Bc7t/WSoDG0EMGIQIjE0wbFV7zERHlo8AFIPAHUI0TeCMI5Z++q0X/MuWrNqppVtfauqr1X7Rrfz6c+tWquueaaVbX3qLXG/I0xZGY4juM4juM4juPMO431noDjOI7jOI7jOE4d8Jsjx3Ecx3Ecx3Ec/ObIcRzHcRzHcRwH8Jsjx3Ecx3Ecx3EcwG+OHMdxHMdxHMdxAL85chzHcRzHcRzHAfzmqAtJx0gySQsV+v6kpDvXYl5rjaS3SvrzMY7/iKQzJzknx3G6Gef/dKPYL0mflvS6MY5/QtL3TnJOjuOUzOt1laT3SvqDVR77akkfm/ScnOpsuJsjSadLulHSdyXdH7d/RZJWMdbzJH1M0iOSHpX0eUmnApjZ35rZ90/+HcwWuQs0MzvFzC5brzk5Tt2IF/GPSNq8Tuc3SccVr+fRfuVupMzsADO7a73m5DizwEa+rlLgLkm3r+V5k/P33Tya2fvN7F+vx3ycwIa6OZL0W8CFwB8CzwGeDfwn4EXAplUM+b+Bj8dxngX8GvCdiUzWcZy5QNIxwE8CBrx8XSfjOI6zAubguuqn4jy+V9L/tY7zcGrEhrk5knQQcB7wK2Z2jZk9boF/MLNXm9ne2O/nJP2DpO9I+qaktw4Y71DgWOB/mdm++Pg7M/tM3P9iSbuT/vdI+h1JX4relUskPTtKzB6X9AlJzxgy/x2Sbo3z+idJJ8f2wyXtlPSwpF2S/mNyzFslXS3p8niO2yRtj/vOkXRNzzkulPRHo8btOabrfSbv9WfjHH8PeGWUp3wx7u94aCU1JP2+pK9Hj9Pl8btKPSZnSvqGpAclvWnQZ+Q4M8ovAZ8D3gt0yU2j9OKdkv46/g/fKOn7kv0XRjv1nehh/cncCeLxv9rT9iVJvyDphtj0xfh/+sqM/TpK0gclPSDpIUkXDThPU9LvRRv1eJzTUXHfj0u6WdJj8fnHk+M+Lem/Sfq7eNzHoo1F0kclnd1zni9K+sVR4/Yc07WKnXpkJZ1PuEG9KH4GF8U+nRU1SQdF+/RAtFe/L6kR971W0mck/b8KHu+7JZ2Sm4fjbBQ049dVFTkT+BBwHf32eaDdivv/QtI/R9t0g6TnDXjfX5H088nrxXi98wKgsM+PRtv0Y4W9Sfo/T9LHFa7Xvi3p98Z8z84INszNEfBjwGbCH/kwvku4WNkK/BzwnyX9QqbfQ8Au4M/jBcazK8zh3wIvBZ4L/DzwEcLNw6GEz/rXcgdJOgG4HPidOK+fAu6Juz8A7AYOB04D/rukE5PDXw5cGY/bCVyUHHeqpAPjOZrAvwOuqDjuSMzso8B/B66K8pTnZ7q9Nj5+Bvhe4IBkjgU/AXw/cCLwZkk/uJJ5OE7N+SXg/fFxUsaWnAH8V+AZBJtzfrLvZuAFwMGE/92/kLQlc47LgH9fvJD0fOAI4Doz+6nY/Pz4f3pVemC0DR8Gvg4cE4+7csB7+c0431OBA4H/ADwp6WDgr4E/Ag4B3g78taRDkmNfBfwywUu7Cfjt2H5FHLOYz/HA98Tjq4w7EjN7E/C3wNnxMzg70+2PgYMIduqnCd/bLyf7XwjcSbDn/wO4RFq5rMhxZoiZva6qgqT9Cdc/hX0+XVLvatggu0Wcy7a47wtxjByXk9hngv28z8xuJVzvAWyNtumzPXN8OvAJ4KOE67XjgE+u4G06q2Aj3RwdCjxoZstFg6S/V9C0PiXppwDM7NNm9mUza5vZlwg3CT/dO5iZGeGC/h7gfwL3Rc/AtiFz+GMz+7aZfYvwQ3xj9LDsBa4FfnjAcWcBl5rZx+O8vmVmX40e2Z8A3mhme+I/0nuA1yTHfsbMrjOzFvA+4Plx/l8n/LMWBuolwJNm9rmK406KVwNvN7O7zOwJ4FyCAUqDM/+rmT1lZl8Evli8B8eZdST9BOFC/2oz+zzwT4Qf25QPmtlN0Xa9n3AzBICZ/bmZPWRmy2b2PwkXKjlN/oeAbYl9eg3BabGvwjRPIPzo/o6ZfTfahM8M6Ps64PfN7M7oQf6imT1EuCD6mpm9L871A8BXCRczBX9mZv9oZk8BVyfv81rgBZK+J75+dfxM9lYcd2ziDeIrgXOjd/wegt1PbeLXzex/RVt7GXAYQRrkOBuVWb6uqsIvAnuBjxEcRAsEm5MyyG5hZpdGe7EXeCvw/Lja1sufkzirCXblfRXn+DLgn83sf0bb/LiZ3Vjx2JFIulRB1fOVCY330fj38eGe9mMVlBFfk3RV5ia0Vmykm6OHgEPVHdT242a2Ne4r5BEvlPQ3UTrxGEE7e2huQDPbbWZnm9n3ES5wvkvwAAzi28n2U5nXBww47ijCRVMvhwMPm9njSdvXCZ7dgn9Otp8EtiSfQeqRfRXlqlGVcSfF4XHs9DwLdF9U9L6HQZ+T48waZwIfM7MH4+sr6JFuMOTvX9JvSbojyjYeJaxs9Nmr+ON8NfDvoxTsDKr/+B5FuPBfHtlzuK36ek/bKFt1AEC0Q38NnB73nU7pga0y7iQ4lOAV7rVV2fmb2ZNx022Vs5GZ2euqKL17Ij5ePWDsMwmOq+VoQz9IRfusIDG+QEFi/B1KtU/OPt8L/B3wbyVtBU5h8CpTL4Ns7qR4L3DyBMf7Q/KO9rcB7zCzbcAjhEWB2rKRbo4+S/AA7BjR7wqC/OwoMzsIeDcwUhphZt8E3gn80JjzzPFN4Psy7fcCB8dl1YKjgW9VHPcvgBdLOhL4N5Q3RysZ97vA/sWL6GF9ZrLfRszhXoIBTM+zTLeBc5wNh6T9CFLWn4669H8GfoPgXRy5OqoQX/TGOMYz4gXJYwy2V5cRVl1OJKwSf3ZAv16+CRytCql2GW6rvqenbSW26gPAGZJ+DNgP+JtVjNtlqwjB4ynDbNWDwBL9tqrq/B1nIzKz11UWsuYeEB99NyLxuuglBIdSYZ9PI6zwZG/sengV4XP5WYLT6phi6AH9C+nzK4DPxpUwGH0NNcjmTgQzuwF4OG2T9H1xBejzkv5W0g+sYLxPAqnjnSg/fglQxMFfRqlqqiUb5ubIzB4l6PbfJek0SQcoJAN4AfC0pOvTCasme2KsT6/EBQBJz5D0XyUdF8c5lKCv/9wUpn8J8MuSToznOkLSD0TD8ffA/yNpi6R/SbjbruRxMLMHgE8DfwbcbWZ3xPaVjPuPhNWon5O0CPw+QdpT8G3gmOitzvEB4DfikuoBlDFKVbzUjjPL/ALQAo4nSDFeAPwgQRrySxWOfzrBkfAAsCDpzYQ4nyzxZqhNkKv0rhp9mxBLk+Mm4D7gAklPizbhRQP6vgf4b5K2KfAvY/zPdcBzJb1KIQHCKwnv+8MDxunlOsKNyXkE+9BO2quOeyvwU5KOjtKWc3v2D/wMolTuauB8SU+PEr/fJMhhHGcumfHrqlG8hnB98/2U9vm5hFjsM4YcV/B0wo3jQwSnzH8f0f+vgH8F/DrdK2UPEOz2IPv8YeA5kv6LpM3RPr2wwvzG4WLgV83sRwgxVu8ac7xDgEeT677dTEepNDE2zM0RgJn9D8IP2u8C9xN+DP+U4H39+9jtV4DzJD0OvJnwg5hjH8ET8AlCmsmvEP4RXjuFed9ECPh7B8Ez/P9TejDPiPO4l6CvfYuZfXwFw19B8Gxc0dNeaVwze4zwmb2H4EX9LuEPu+Av4vNDkr6QOf+lhAu1G4C7gT3Ar2b6Oc5G40yCXv0bZvbPxYOQkOTVFVZqricE/P4jQeK1h+BFHMblwL+g/6L+rcBlUQv+79Id8cbg5wmBvt8g/H+/csD4byfYzI8R7OIlwH4x7uhlwG8RLhZ+F3hZIiccSiJp6bJVKxk32q+rgC8Bn6f/BupC4DSFbHN/lJnGrxLs213AZ+I8Lq0yf8fZqMzqdVUFzgTeldrmaJ/fTb+0LsflBLv8LeB2RtzgxZilvyRk6/tg0v4kIQnP30X7/KM9xz1OSEjx8wSJ39cIcVtTITqxf5yQ/OdWwnd9WNz3iwqZ93of148aNtM2asVsXVGIj3Mcx3FmHUm/BLzezH5ivefiOI7jlMTV/+ea2b8f2XkNUajF92Ez+yGFpBF3mtlhY4z3YuC3zexl8bUIK2TPMbPlKJ9+q5mdNPbkp8SGWjlyHMeZVxTS0v4KQRLhOI7j1ASFsgRnUXP7bGbfAe6W9AoINzZVYmRHjGmEONLTYlNRW6q2+M2R4zjOjCPpJIJn7tv0S2gdx3GcdULSfyRIoj8SEyDUBkkfICTe+H5JuyWdRUjsc5akLwK3MTohRzre3xLCLU6M4xWrQ28EflPSLkIM0iUVxztKIRPiHZJuk/TrmT6S9EeSdikUDP5XVec78Lwuq3Mcx3Ecx3Ecp05IOgw4zMy+oJBh+fPAL5jZ7UmfUwkxo6cSinVfaGZjJa3wlSPHcRzHcRzHcWqFmd1nZl+I248Dd9Cf6W4HcLkFPgdsjTdVq8ZvjhzHcRzHcRzHqS0xccQPAzf27DqC7kyuY6cKr1L0b2Js0mbb0pUaf0wGlhgbXHssu0cja5WtDRmJ42jRY6aHKyVHsofvss/2Vv7iT/qZp9lDD7cG7v/8l/Zeb2aTrDLtrCFum0bgtmnNWKltArdPG51DDz3UjjnmmPWeRj/tdrndGOFrL/qO6jfOPNKxp3m+Oeaee+7hwQcfrGyfjpPsySH77wsxTXuSpovNrC9pRUwx/pfAf4lJI7p2Z4Ye69dmTW+OtvA0XqgTuxur/vhnaoyqMeDYXD3S2Fe58w3655nmhUku1is1NJ1usV97wPdsmWNyfTP9Ks9rA3KjfXJF/R98eJm//+hgR8SWw++uUlHbqSlumxLcNq0rK7VN4PZpo3PMMcdwy003rfc0+tmTXNNu2VKt76h+48wjHXua55tjtp9wwor6PwW8Ycj+34c9ZrZ92BiSFgk3Ru83sw9muuwGjkpeH0mo4blq1vTmaFVUvfBI+2X2dy48chcb6YVGZn/2omVMsokwivMk+1RclDTLOXQd2859PrmLksz7zl2UpO91Ti5GqmBA293eTorbprDhtmndcfvkVGJ5OTwvTOjSbyU3HtO8ScmNXbSt5AZu0p+PgxgvfifWSLoEuMPM3j6g207gbElXEhIyPGZm941x2hm4OXKcGmAYSzZYtuI4jrNeuH1yHKeujHmj8SLgNcCXJd0a234POBrAzN4NXEfIVLcLeBL45fFOuZ43R8M8npP2yELpdc20dfcbcZ5h/XIMkJxoiGe0y/ua89h2HdPuP6bwxCbnKD6/LlmL+vvl5uNe2oB7ZucEt02ZRrdNdcftkzOS1ayI1EGe9uij5fbWrdWOGTXv3P7i8ylWkNI2Z1WMu3JkZp9hWLAunSKzw9R7K6bSnCVtlXSNpK/GQkw/JulgSR+X9LX4/IxJTsxx6oQBS7QHPpz1wW2T47h9qitun5x5R4RVmEGPulL1hu5C4KNm9gPA8wl5xs8BPmlm24BPxteOsyExoGU28OGsG26bnLnH7VNtcfvkzD2NIY+6MnJukg4EfooQEIWZ7TOzRwlFly6L3S4DfmFak+yfVCM8Giof3XMOcpRGo3xI4ZG0qRkeNJsjHo3ysbAQHo1m+SjaVvNIz7O4AIsLaCF5LIZHV7/ivSTvp3jP4X3HR9cx4TNTQ51H3+eZk+iED7R8zCmGsTTkUQVJJ0u6U9IuSX0/iJI2S7oq7r8x5vQv9p0b2++UdFLSfqmk+yV9ZcA5f1uSSdpw2arcNrltctsUGNc+uW2aPLW0T6thy5bpS+oefbR85Ni6tXxUZdS8h+1P7aAzFht55eh7gQeAP5P0D5LeI+lpwLOLbBDx+VlTnKfjrCtmsDTkMQpJTeCdwCnA8cAZko7v6XYW8IiZHQe8A3hbPPZ44HTgecDJwLvieADvjW25cx4FvBT4xore7OzgtslxGM8+uW2aGm6fnLlnVm+OqsxtAfhXwK+a2Y2SLmQFy8CSXg+8HmAL+4/ovIJ6If3nKV8UAc7NZv/+Xq8kBK9rbg65+iMrqUXSS1etkDi3RPLQCVxuJVmHrJhDco4kSNla7b55qVN/JOlXHJ6kzR2ZVndYOt25k2qI1vCYwFGcAOwys7sAYsrJHcDtSZ8dwFvj9jXARTGN5Q7gSjPbC9wtaVcc77NmdkPqxe3hHcDvAh8aZ+I1xm2T26aSubVNMKZ9cts0HSZmn44++mja0ZfdmIUYsmIFaNRKT9V+dSBN0lDgK0uVqLN8bhBV5rwb2G1mN8bX1xD+4b8t6TCA+Hx/7mAzu9jMtpvZ9kU2T2LOjrPmGLBkGviowBHAN5PXu2Nbto+ZLQOPAYdUPLYLSS8HvmVmX6wyuRnFbZPjMLZ9cts0HSZmn575zGeuyYQdZ9LM6srRyJsjM/tn4JuSvj82nUjwKO0EzoxtZ7KxPUDOnGNAK3pncw/gUEm3JI/X9wyRu0LpdXEP6lPl2HIQaX/gTcCbB76hDYDbJscJjGmf3DZNAbdPjhOYxYQMVW/cfhV4v6RNwF2EAksN4GpJZxF0w6+oNNJqA2dHSUpS+Yhy++N2ImfpzCUncUnbc9KVxgg5S46ueh/W19aRnCTvxQr5SVoDpF1KWzqSlkTu0pHApHKWuN/Sv8YhNUdC85BaI3NWrT54Zof+Kz9oZtuH7N8NHJW8PhK4d0Cf3ZIWgIOAhysem/J9wLHAF+Pf85HAFySdEH+wNxJum3qPcdtUbs+BbYKx7ZPbpukxMfs0E3K6gqoyubWW042qlzSsNpJL6FZFsXI0a1Sas5ndCuQM64mTnY7j1BNDtMbzc9wMbJN0LPAtQhDzq3r6FB7FzwKnAZ8yM5O0E7hC0tuBw4FtwE0D52r2ZZIgX0n3ANvN7MFx3kAdcdvkOGPbJ7dNU8LtkzPvjFsEdr2o5Q1d1UDnnEcWQIVXtZHxtDYST2wR7JzzyEIZDJ3xAttqvLMp0aOpLq9ru2sfgIqg5q5A6MR7W+xv9Hti03kVIyr1pBYV7JOh85XrPRC6gmd2+PFmy5LOBq4nRLxfama3SToPuMXMdhJSvr4vBjU/TLhIIfa7miDJWAbeYBa+NUkfAF5MkM3sBt5iZpeseqLOUNw2uW2qI+PYJ7dNc06xWgLTT9k9iF27yu3jjhvcNi6jVqqqvv/iM0tXk3xlaSCz+MnM4pwdZx0QrTFujgDM7Drgup62NyfbexggsTCz84HzM+1nVDjvMSudq+M4s8R49sltk+M408BXjhxnA2PAEs2R/RzHcdYat0+O49SRDR1zNFUGVT7v3Z9KMwqpxKCg51yAc5SsKK0b0pGzpG3pOFGmkmnLBj2vRMISZR6Wyj1aQ6Qr7eS9LKdBzzGYuZ2ZTyuRoXRO2y8vSWdtadrXKKtJa47MayC0mVgyv/iYK9w2lW1um2qN2ydnRdRBSpeSk81NSko3KVbzmaW1keZYdjfOO5d0KfAy4H4z+6HM/hcTMj7eHZs+aGbnjXFKoA43R44zA4RUubO4OOw4zkbH7ZPjOHVkArK69wIXAZcP6fO3Zvay8U7Tjd8cOU4FDLFk/u/iOE79cPvkOE4dGVdWZ2Y3SDpmMrOpzmxa00b/fWhXDZBC7tKVJSq2pXKWhX7pii2kcpfimGR/p61fKmJdUpoR7yEqP7oyNMWMUGqVbZ1aIqlcJX3/rUY8JpGSdOaTtBX1V9KaI5lpqZ3IVIrTJBmhChlL5SxRsGFkLK3RleadecdtU3mM26Y1xe3THFPIt6pKt+ogpZs1cp/ZIKld7vtY6Xe0gRjxk3OopFuS1xeb2cUrPMWPSfoiocbab5vZbSs8vo/5+5YcZxW4Z9ZxnLri9slxnDpSYeVoWIHqKnwB+B4ze0LSqcBfEeqtjUWtrGmnhsioivOdfZlA53Q7V3G+65gY1LyQ8dgCVhy/kHhnm/3jFG1d3tlRTryiUHxXgHMMhE48pB1Pbfpe0roiy9E7m3hvi/cqlcGAXUHRPVMc5DstPLW5yvWVA6GT+cyyl9Y1/fON2ya3TXXG7dOcU9fViGnUKhqHSc9n0Apc7vuo63e0BkxzTdvMvpNsXyfpXZIOHbew9Px+W46zAoJn1rNBOY5TP9w+OY5TRwQsTnN86TnAt83MJJ1AUPE9NO64fnPkOBXwVLmO49QVt0+O49SRcRMySPoA8GJCbNJu4C3E+y0zezdwGvCfFeQITwGnW64uxAqZzZsjZeQsXXKXZv/+jgwlkaYsZNoWh0tX2gtFfZF07ChdSeQzHRlLV6GOzFtpJ9KV+H1qub+WiJYTWchyrh5Kkk8/BmmndVO0FPerf44ikbh0TS4TKD2ngdAGY1Wgd+YEt03978tt09Rx++TUkjpI6VLWej65hA1zWPtoHMtkZmeM2H8RIdX3RJmPb8ZxxsRlK47j1BW3T47j1JFxV47Wi/Wb86jq871kgpWzgc7pdhoo3MxUlI9e2S6PbGY7Paa92O+dbS/0e2c7XtlRkWhpltsiEDrxxDaKdLipxzbdvxQ9p8l8tBTnmKYLLraXkrb4nPpJlaTYtcLDkQZm5wKh4xSUeqdHeWpnEA94nhPcNsWDk0PcNtUet0+OUzNyCRvmZLUoZRaLDMzft+Q4q8A9s47j1BW3T47j1JFpJ2SYFu5qcpwKGNC2xsBHFSSdLOlOSbsknZPZv1nSVXH/jWlVaEnnxvY7JZ2UtF8q6X5JX+kZ6w8lfVXSlyRdK2nrat+74zj1Zlz75LbJcZxpIMKNxqBHXan9ylG2hkhnZyPt2H9MVu6S1gBpdj1Dt3SlE+C8WB7T2pSRrixmgp6LQ9I431TZUqhBEt1IoRpRUvejsRQDoZPK9I2l/jomjSQwu6hS31WtXku90ymDmpMmSwKgFf88RlarL76HgbVEMvtnrL7IuNmgJDWBdwIvBXYDN0vaaWa3J93OAh4xs+MknQ68DXilpOOB04HnAYcDn5D0XDNrAe8lBCNe3nPKjwPnmtmypLcB5wJvXPUbcPpw2+S2qS6MY5/cNjmOM01qf6ORodKNm6R7JH1Z0q2SboltB0v6uKSvxednTHeqjrN+GLBkzYGPCpwA7DKzu8xsH3AlsKOnzw7gsrh9DXCiwtX0DuBKM9trZncDu+J4mNkNwMN98zX7mJkVV5OfA45c0RueEdw2Oc7Y9slt05Rw++TMO0VChkGPurKSuf1MT8XZc4BPmtkFcRn+HCbl/WlkPLI572vX/oynNvWWFoHLacX5wsOaSYcL0N7UjM9p0LO6nsMxsS0Zuvg9GlmZvss7G140Wsl73RTb0my4+xLPaBFwnbQ14vtupCl9O97i4Z9jzlOrRJaR9dQW+1pJW1eV+iKaOxMIPTMpdDUqVe6hxY9f5GIzuzh5fQTwzeT1buCFPWN0+kSv6mPAIbH9cz3HHrGCyf8H4KoV9J813Dbhtml+bROMaZ/cNk2XtbNP4/Loo+F5qysd15QNnt573hIy7CAUZoLgUfo0dfkHd5wJU3hmh/CgmW0fsn/E5efQPlWOzZ9UehOwDLy/Sv8NgtsmZ64Y0z65bVpb3D45c8NGT8hgwMckfV7S62Pbs83sPoD4/KzcgZJeL+kWSbcssXf8GTvOOmCItg1+VGA3cFTy+kjg3kF9JC0ABxFkKVWO7UPSmcDLgFdPomJ0TXHb5Mw9Y9ont03TYyL26YEHHlij6TrOZNnoCRleZGb3SnoW8HFJX616grh0fzHAgTq4zwiqq/7GSuuLDKgl0qnM3h/0bGnQc5SptBf720J7v0yltTlst3LSleTT7EhXUjlLTnHTJV2JQcipBCTKWJpLZcdUXlO0pzKVZqM/CLtZSETU/5kN+unM1hqJcpbsr1n6G5etUj+7v4FmIz2zo7gZ2CbpWOBbhCDmV/X02QmcCXwWOA34lJmZpJ3AFZLeTgh63gbcNOxkkk4meCN/2syeHGfiNcdtU8Rt03zaJhjbPrltmh4TsU/bt29fmz9Ql9OtDxtQSldQxBzNGpV+8c3s3vh8P3AtIeDy25IOA4jP909rko5TB8ZZOYoByGcD1wN3AFeb2W2SzpP08tjtEuAQSbuA3yRo0TGz24CrgduBjwJviNmgkPQBwgXL90vaLemsONZFwNMJP8i3Snr3ZD6FeuG2yXECq7VPbpumh9snx9mgK0eSngY0zOzxuP2vgfMoPUkXxOcPTXOiXTT6PY1dAbw5T2QR9Jx4MS1Tmb7wtELptW1tSryzmzLe2U30tWW9s+lfQsb12UmX2+WdjedYKtsaSSX59r4wQLNZDlSk0G12nS8GQmfioFNyP6PdgdBFWxlAaDF1bvoddL3XdpG+Nw2Enq0UupMosmhm1wHX9bS9OdneA7xiwLHnA+dn2s8Y0P+4sSY7A7htctvktikwrn1y2zR5ammfHGeNGXflSNKlBAnu/Wb2Q5n9Ai4ETgWeBF5rZl8Y45RAtTk/G7g2/rgsAFeY2Ucl3QxcHb1B32CA4XScjYAhltvj3Rw5E8dtk+Pg9qmmuH1y5p4JyOreS75eWsEpBDnvNkKWzT+hP9vmihk5ZzO7C3h+pv0h4MRxJ+A4s0J7JhNSblzcNjlOiduneuH2yXEC48jnzOwGSccM6bIDuDwmdvmcpK2SDiuSnqyW+sdJZarHjz6mP+i3kEV0BT0X0pWMXCVsZ4Ke43Zrczl0u5CzbCrbLH6y7VHSlZSi5Eai5ihqiLQTuUpjX3LuznlS2YzF52TouH8h8zEO+mQ7I6ZSkmI7aStqoHQFR7dK/U3nfSeSnLLjbNQXMYMl98w6KW6bwjhum9Ydt0+OMycUNZFmJIlDhZWjUTUiR5Gr03YEsMFvjhynBhSpch3HceqG2yfHcepKtjh6gdmoGpEjh8+NOsZ4gN8cOU4lDFgeXoHecRxnXXD75DhOLZFgy5bB+596atwzrKrW2ig2zs1R7s60K0PBNVtcAAAgAElEQVRRLiNUlKZk2qCUhWSzP21O64vEfUkZ4KJfeyHN1JTMLfc7VkhXEolHYznW+0gyQjUTiUx7b1HHpD8jVCrTyd1b5778rmlFSYpSmYr1y1Q6/ZIm65K2xFoiiUxlFuuLtP3iw1kNbpv63oPbpsnj9slx5oAZkdN1kKY9553A2ZKuJCRieGzceCPYSDdHjjNFzOSeWcdxaonbJ8dxakmjMXzl6PHHhx4e66W9mBCbtBt4C7AIYGbvJpQgOBXYRUjl/csTmPWM3hwVXr5BVeuzFdfpbyu8mKlHNlddPvG6FrVE2qmHtOOxTdo2Wxwj8ZqmVeobhSs2mXcn6DmpkbIcq8wvJ/PaW243F4qxM3VM0nji+FlZxoud/hEodZbmApzjduFxDefO/BnlgpUzVepnpb6IActtv/hwRuC2qXwLbpvWDLdPjjMlPvKR8HzKKes7j1lllKxuBIPqpSX7DXjDqk8wgNm8OXKcdcADnh3HqStunxzHqR3Tl9VNhdmbseOsA4bLVhzHqSdunxzHqSWjZHU1ZePcHOVqjaQyDfXLVKxTXyQNji4P6UhXkk+pkLGkdUMKyUprcymvKKQrtimVriSSjFjvoyvKuFCKJNIVWnEOiXSl0RWEXcyxX7rSFfSdiwnP/Zh2TdH62goJSTYQOhkvHduizEXKfM6zUl/E3DPrrBK3TeUwbpumg9snx5kOLqcbD185cpyNi2v6HcepK26fHMepJWPGHK0XG/PmaEjBqa6g32IzdQY2+z2Iqeez3QkyLo9pRw9sO/XO7he9iptK12ZjsXRFNqJ3Ng367cT3Jt7ZdvzBay+Vk2wtltvtvTGYOfFOl+lyk/eV/UyKIOOyRa00zW14kw3rb6OVuGwLr0CuWj1Jlfp0f1GlfkZS6HqRRWciuG3qmj+4bZoEbp8cZ4Pw6KPl9tat6zePSTGjK0fuanKcirSsMfBRBUknS7pT0i5J52T2b5Z0Vdx/o6Rjkn3nxvY7JZ2UtF8q6X5JX+kZ62BJH5f0tfj8jFW/ccdxas849sltk+M4U6GIORr0qCl+c+Q4FTCDVrsx8DEKSU3gncApwPHAGZKO7+l2FvCImR0HvAN4Wzz2eOB04HnAycC74ngA741tvZwDfNLMtgGfjK8dx9mAjGOf3DY5jjM1ClndjN0crftaVypTKGUczXzn8qCVnygj3cgGPaextc3uZyglK+3FRKayWEhXknlFyUpzy3LZtKncXlwI0o2FZiJniedOlRvLrXDypeVyEvs2lV9bazFsL3dVoS9qraTvWT3PZd0QJT+eSst4FNvpd1RIWxLpijKB0F3yk0Z8j83yPRSzsDTqeVgAdK6+CKxhAPTYspUTgF1mdhdArOa8A7g96bMDeGvcvga4SCFSfAdwpZntBe6WtCuO91kzuyH14vaM9eK4fRnwaeCN47yBecNtk9um8s3QT21sE4xpn9w2OU6d2bOn3K7xDUWWGZXVzd6MHWcdMBjlgT1U0i3J64vN7OLk9RHAN5PXu4EX9ozR6WNmy5IeAw6J7Z/rOfaIEVN+tpndF8e6T9KzRvR3HGdGGdM+uW1yHGc6eCrvKVFUO2+O8NimVPXYZRxtXYHCjcFtuRS6JB7b5ubgaty8eanTtn+yvXlhuesZYEH9Xud97fC+97XK9//k3jJX71ML4eT7FhY7ba1GnJwyP5aJd7Godt/tkU32t2Kl+HbiVY1e2UbXBxDT4SaV6dUqXa0Wj1fisS08410pdJWrTF+TAGgb+Wf1oJltH7I/59btHXFQnyrHOmuN26bw7LZp/RnPPrlt2igsl/+zs+itd8gnYZjBm4sOM7pyVDnmSFJT0j9I+nB8fWwMzPxaDNTcNGoMx5lVjLETMuwGjkpeHwncO6iPpAXgIODhisf28m1Jh8WxDgPurzLJWcRtkzPvjGmf3DZNEbdPzlwzozFHK0nI8OvAHcnrtwHviEGVjxACNh1ngxI0/YMeFbgZ2BZ/GDcRgph39vTZCZwZt08DPmUhx/BO4PSYMepYYBtw04jzpWOdCXyoyiRnFLdNzpwzln1y2zRd3D4580uxcjToUWmIkdk0XyvpAUm3xsfrxp12pZlJOhL4OeB84DdjIOZLgFfFLpcRgjX/pPKZO9XFc5XQk1X5orZHIougmQmEba8iELqYyqBbxFytkYychRhwrMVyDguLYXl7y6ZSrnLA5r2d7act7gNg/4V9nbZNjf6o3+V4oj3LpTTlu4vlmI8vbg7PzfLcexphPi3KY8qJp9KV7mco5SpQxip31ReJwddKg56Xh9QXgVLvkX6vRVuu/kiXdr6/5kg2AHoNgp/b7Uo3QVmiTv9s4HpCVP+lZnabpPOAW8xsJ3AJ8L4Y1Pww4SKF2O9qQoD0MvAGM2sBSPoAIbj5UEm7gbeY2SXABcDVks4CvgG8YtWTrzFum/r7um2aP9sEq7dPbpumx1Ts0zBmUL60ofnIR8rtU06Z7NizkqRhzJijJJvmSwkr1TdL2mlmt/d0vcrMzl79RLup+p/0/wG/Czw9vj4EeNTMCoFrlSBMx5lZilS5441h1wHX9bS9Odnew4ALBTM7n/AD29t+xoD+DwEnjjPfGcFtkzP3jGuf3DZNDbdPznxTyOpWT5VsmhNnpDWV9DLgfjP7fNqc6Zp1j0l6vaRbJN2yxN5cF8eZCcwGP5y1x22T45S4faoXk7RPDzzwwFTm6DhTZ3xZXS6bZs6h8G8lfUnSNZKOyuxfEVVm9iLg5ZJOBbYABxK8IVslLUQPyMAgzJgu9GKAA3Xw2prpnCxiFLksUWnpiuJ2MlVSNMPYjYVSUrG4GCQX+y2W2WMKuQrA1k1PAXDgYrk0ul8j7G+onGuhF3+qXcZsfmepvAvf3AzjLyRZlB6L0pUnk/dQyFi6Mz6p6xmgsVyee7mzP5GzdKQrqZwlnjupd5LLDtWpcQIofpCW9itquyS/FYoFVtY7M5Qh2mOuHDkTx22T26Zk/3zaJnD7VFMmZp+2b9++/n9kK+HsqG666KL1ncd6M2kpXUqdpXQpo2V1o8qgVHEo/G/gA2a2V9J/IshVX7Kq+UZGWlMzO9fMjjSzYwg640+Z2auBvyEEZsJ8BFU6c44NeThrj9smxylx+1Qv3D45TmT4ytGDZrY9eVzcc/TIjJhm9lAsRA3wv4AfGXvKYxz7RuBKSX8A/AMhYHNtKDx1g7yvOU9splJ68avRVcJj1C9JEWPbSDpGb6qStsJbmtYKOWCxlO5sXXwqPpc+1P2jd3axkdQqiOxtlwHMT18o78L3a4YA6E2ZYyxxKz8ZvbKtxBNbVJ9vLJdty8m24pCNpcR7u7k4pryvtqXosV0u/5y6vK6F1zYNTC8CpBupmzvzvXa8oUlAeC4AetqV6Q1sjIQMzpritqnT5rYp9N/AtgncPs0W62efJkWRaGDQysi8rxg5JeMXge1k0wS+RXA0vCrtIOmworA08HK6s0OuihXdHJnZp4FPx+27CIFSjjMXWLWU3c464LbJmXfcPtUXt0/O3DJmEdiK2TR/TdLLCRkzHwZeO+60Pe+j41TAGC+Vt+M4zrRw++Q4Ti0ZP1tdlWya5wLnjnWSHmp/c2Q5ycmQfgDKSVsK1UQSPNsZs1+F0tc+tC2SBisvxNoeRVAylDITgKctBBnLQQuldOXpjRAAvaVR9mvEiS9Z+VUd0Cz/0Iq+i5k6JOlUix/OpxLJSSFTUXm6LhlLMY1WKl2Jxzf2JQHMmwrpSkauArBQBEqXcywCoJUUZSnkLkpkKEUAdBH8HPqtg4reZStOD26b3DZBDWxTmIzbJ2ftmGaiAWdjMb6sbl2o/c2R49QGj2x2HKeuuH1yHKdujCmrWy9qNePC66ZGGoXc7O9XeGwHeekywbPqeHnLbsXxXeNkvLNdVdqL/SP03YodU49tGphcBDgXzwBbm0/GtjI4uhlP3kq8mE9rl/u3RNdqI5l4K85tX7v8epda4XNcXio/z6XodW3tK/u1yunQ2Bf2L6dt0WPb3JQGPcfUt4tJutzkPNYqvLfJd1kEQKeV64sA6FwAe7/zOZ6oOGbalenlntk5xm2T26ZwQB1tE7h9clbEo4+W21u3rt88nI2Prxw5zgbGPODZcZya4vbJcZy64itHjrOB8YsPx3Hqitsnx3Hqhq8cVaSQE3TkBankIFOTtp3Tj2SObffLHZTIIjpB0Wk9i6IMRTp0GhQdZQpdtUYyxxQ/SqnnrthOpSvNZHtRQYvxtESmUkhWDmzsSfr11whJjynGSVmyIBFJpSt7W2F7z1LZVshY2kvl595KgpmLGiKFhCXsj8+bk37xGC0mgdALyXYxflqFPla2tySYuVOFXrm25NhE2rSmAdCu6d/YuG0C3DbNpG0Ct0/zzHL8X6zqoZ9hKd39zWAbntUapGetKaNqQ21UZjTmKPOL7zhOHwa0NfhRAUknS7pT0i5J52T2b5Z0Vdx/o6Rjkn3nxvY7JZ00akxJJ0r6gqRbJX1G0nHjvH3HcWrMmPbJbZPjOFOhSOU96FFTZu92LqXLI5t6XeN2GvQaPbVqJd7Xoi0ZJ806W2ynDtDCK9uVVrdzuvJHqBiyPULq0KCc96Z4oi1J/tr9k9S5BVus9Ng2o7uwRXmePRYq1j/VKivXP75pMwBPbinb9u4N2/s2l8HIrc1Jaty9Ybu5qTx3e7F4Tt5rDIC2xMubDYBOg56LD7eRaUsCvCm8Q4kXd2AA9JRJ/8RWiqQm8E7gpcBu4GZJO83s9qTbWcAjZnacpNOBtwGvlHQ8oSr084DDgU9Iem48ZtCYfwLsMLM7JP0K8PtMoDCaUxG3TW6b1pjV2ie3TRuAYZ75PeVqb50vRqsycytGBfO2YlQwo7I6XzlynKqYBj9GcwKwy8zuMrN9wJXAjp4+O4DL4vY1wIkK2p0dwJVmttfM7gZ2xfGGjWnAgXH7IODeVb1nx3Fmg9XbJ7dNjuNMh0JWN+hRU+o7M8epE9YTy7FyjgC+mbzeDbxwUB8zW5b0GHBIbP9cz7FHxO1BY74OuE7SU8B3gB8da/aO49SX8eyT2ybHcaZDIaubMWbm5sgydUF6OvRvtzJyliTouVNLZDmVs/Rvq0vOEgOhk1hkxbZ2q/TQLcf6GUUND4CldrndjlKTdrJ4V8hY0kDmxdjWFTCdyF2IAdL7rBz7yWaQqTyxWP5BfmfTfgA8vq9se2JzkMUsbS7/DNpJMHN7k8XnJOg5ylgSVQzNKGNpJ0HNjSTo2RYyQc+xboiSNoufTyolsqK+SLqUngmAtnayCFroS5Lg6fHrioz0wB4q6Zbk9cVmdnH3AH30TmpQn0HtuZXfYszfAE41sxsl/Q7wdsJFiTNh3Da5bUoOTjbXyjbBmPbJbdNGYTn55y+88jN4YepMgEJOud7f/4wmZJi9GTvOejHcM/ugmW0fsn83cFTy+kj65SRFn92SFgiSk4dHHNvXLumZwPPN7MbYfhXw0aGzdxxntlm9fXLb5DjOdJjRmKP1uznqTZsLHa9a6mkrU6P2VzDv8timHr0imHmhP8CZnPc1DXpeTre7n6H0yjaWk5SuMa1smnZ2KQb4FmlqoTsIeU+MHt5nSaX4+B5bGWfcYuLISz21bfrT7j49emwPaJaBmAcuPBX6Le7fadtvU/iDfWqxjGpeWkw8tXG6WU9s0tZeKNqSdMGp1zWm3tRC8h0W+5PPsfO30OhvU/J3YuuVs3a8094MbJN0LPAtQhDzq3r67ATOBD4LnAZ8ysxM0k7gCklvJwQ9bwNuInhtc2M+Ahwk6blm9o+EoOg7xpr9POG2CXDbNFO2KZx8tbht2iikHvph6b1XkqRh167wfJwnFZwp6nJDMoGVI0knAxcSfmzfY2YX9OzfDFwO/AjwEPBKM7tnnHP6ypHjVMHK2jKrOjzo9M8Grif8g19qZrdJOg+4xcx2ApcA75O0i+CVPT0ee5ukq4HbgWXgDWbWAsiNGdv/I/CXktqEC5L/sOrJO45Tb8awT26bHMeZGmPGHI2TTXOMWfvNkeNUZkynsJldB1zX0/bmZHsP8IoBx54PnF9lzNh+LXDteDN2HGdmGMM+uW1yHGcqjC+r62S+BJBUZL5Mb452AG+N29cAF0mS2eoDOmfn5iitG9KMHrK0orxlttNA2SLgOD1mOWw3lpJA6OUkcDfKWFKZSlHaIyn30anWnkpXigrve5bLj/jJ5VLvUchYnmyXspE9cXspqa/RjlqZVE7e6NoOc0wDpbc0Qqn4/eMzwH7NMOH9F8qJb1kIYy8slscuLZZnaseq8pbWDSlqiSR/OYV0xRbUdywkAdCNZObFdjORsywXxVsy/dLvcp3qimgdFTNOjXHbVJ6va9tt01ri9snpYpiUaSUXqy6nc8ZhtKxuVDKrcbJpPrjaaY+8OZK0BbgB2Bz7X2Nmb4la4iuBg4EvAK+J9QwcZ+NRVKB3aoPbJseJuH2qHW6fHCeE8O5bbgzrMiqZ1TjZNFdNlSKwe4GXmNnzgRcAJ0v6UYKm7x1mto2gGz5rnIk4Tu2xIQ9nPXDb5DgFbp/qhtsnZ+4xC7lBBj0qsJJsmvRk01w1I1eOombvifhyMT4MeAllRpvLCHq/PxlnMlUpZIRdNUVSaUvRnu4vMkItl1oHRTmEknofjaX+7FBdbVH50VxK6mt05CxJLZEoXXlqXylX+e7S5s72d5ZDbY8nWuXy9uPNsP30dtlvS5SuNAZkhOq0JfsXo55jU1LwZEuceCFhAdjcDPsXF8rPZE8iXbHFWEskyaxVyFO6pSvFc5K1KZGXWJQaFc9hwrGWSJLpqagrYokHtNhvqewlU1ekzByWZBRL5UzFecaoKTJmEVhnwrhtctvktikZxu1TrVh3+zQsW50zn6RZCou/iyn/fbTb3addBavOpjnOSausHCGpKelW4H7g48A/AY+aWfELl1bF7j329ZJukXTLEntzXRxnNnDPbO1w2+Q4EbdPtWNS9umBBx5Ymwk7zhQYZ+Uo/q8UmS/vAK4usmlKennsdglwSMym+ZvAOePOudItY0zN+QJJWwlZZn4w123AsRcDFwMcqIP7+6Q3dx0PWr+nLfW+deqKJAHMXXVFovfOkuBhNWN9kVZ/0LOWy7ZGsr8Rg5gThybNqAxuJQrhxr4w78bexGO7NwY97y29s09sKQOcH9kXvLMHLpS1PYraH/sndUE2FcHMjfLWezFT7W8pqUJfVLZvJP06wdGNpMJ9M2wvNMu2RrP8HFtx21JPbDxN6mm1XFvOU5sNek7aiuDytL5Mgfo9trB2dUU0ZipvZzq4bUrO7LYptM2ZbQK3T3VlUvZp+/btK/9jyq0I+GpS/fnIR8rtU06Z7NhpIo6KmrZxmcDK0VjZNFfLiv5DzOxRSZ8GfhTYKmkh3tXlNICOs6Fw2Up9cdvkzDtun+qL2ydnXpnEzdF6MFJWJ+mZ0euBpP2AnyUsbf0NQdsHQev3oWlN0nFqgctWaoXbJsdJcPtUK9w+OU5gzIQM60KVlaPDgMtildoGQe/3YUm3A1dK+gPgHwiav7WhCHBOg2jbmQDXdqaWSBowGwOgtZQEQi8lAdD7ilojSV2MGNiclOnoyFnaiXSlvTccs7S3/Igf31MGM++/GKQrDy8+rdNWBCRvSQqVFMHMrSRTYbq/YCnRl+yLWpL2iHvfRgyeTktzKA2oLmKHk/1WtCUlQDplWpppv1TaEoOZk++rCHDukqkUcpauWiJx/xrXDOnD3DNbQ9w2uW3qMLe2Cdw+1ZP62adJyel27QrPXgNp8kxaSjeINZJWzurKUZVsdV8CfjjTfhehcq3jzAfuga0VbpscJ8HtU61w++Q4ZSrvWWNmovIsSYdbBEBb4jZMU+d20ukmx1g7E/TczgQ9JxXpm3G7vTdp2xTO2UySW7X3xBSySbX21qbgYWwtli7LPYtl0PNjiyEwbnPzgE7bpkZMX6t+V2Qa1JwGRTeju7BlpUfzyZhu97tJ2t3i+FbymbXjduqR7fLOxm1LnbyFdzZp63hsU49s5piuoOciDW4S9NwJZu5yF2cq2KeB8sV2znubenmLQPpOYH2m/wi8Ar2Tw22T26Zy8PWxTeD2yZkgjz5abm/d2r9/zleM7m8Gm/WsVh2WjeuN2QZdOXIcJ+IXH47j1BW3T47j1AxfOXKcjYxr+h3HqStunxzHqSEbNuZoTSlkCGkgbEdykAngTSvPKxP0nFYhL5Y/EwmEYtCzNdOg53Q71uTYl0hX9oW2hVS6sql4Lucdy4LQ3lRKTpYWk7oii0FWspicu5CupFXm2zHYeY+Vxz49rSsSK82nAc572qHvk+1SKlNs70vKxxfSlVQClCWZj3Wqwqf7e57plbHEY9L9KmQzXRHX/edu9PfLOkiTv4+OtKk9OVeqGP/iQ9LJwIWEYjjvMbMLevZvBi4HfgR4CHilmd0T950LnEUQ6fyamV0/bEyFD+wPCLn/W8CfmNkfjfcO5hi3TeFYt00lNbFNML59cts0x6RXrkUdnJyUbpTUbo5wOV11XFbnOBudMa5nYsaidwIvJVRFv1nSTjO7Pel2FvCImR0n6XTgbcArJR0PnA48Dzgc+ISk58ZjBo35WuAo4AfMrC3pWaufveM4tWe1sUpumxzHmRIuq1tDCq+b0vSsaRX6uG1JgLOiW82aSVv0zrKQ9842oie2uVh6/ooA6PZC2dZcLIKey/ksxCrsiVOV5YUkAHoheEsfayTzpp+lmIN2z6bEs9ssb8M3NzKpc6MHNvXOfnc5eIOfapXj7G2Ffsut5P2N9NQObus6NOOJ7QpmLk7ZlS63qie2vyL91KvRjy9bOQHYFTMVIelKYAeQXoDsAN4at68BLope1h3AlWa2F7hb0i7KbEeDxvzPwKvMwhKFmd0/1uydSrhtcttUbq6RbYqTGcM+uW2aZ4rVolHM+WpRbZixNOqzKqsbWQTWcZzIeEUWjwC+mbzeHduyfWL19MeAQ4YcO2zM7yN4dm+R9BFJ2yrN0nGc2WT19sltk+M4U6FYOdqIRWAdx2GkZ/ZQSbckry82s4vTwzPH9F62DOozqD3n3CjG3AzsMbPtkn4RuBT4yezMHceZecawT26bHMeZCtOMOZJ0MHAVcAxwD/DvzOyRTL8W8OX48htm9vJRY9fz5iitFZELYC4CXNO2dlqnoqduRNo3DaSLAdAdCQtgjfJWVlGeon2l5KS5EKUrSd2Qhb1RKpN8mu0oXSmew9jlHJebofNTqXQlU6iiCFJOJScHLpZ/afs19oW3kgZKRw3JU4l05TtLYen8iaWyvsje5ShdST47a3dpTsK8utr6pliSqVafttuoAOeiLdcvlb2sRyzkaA/sg2a2fcj+3QSdfcGRwL0D+uyWtAAcBDw84thB7buBv4zb1wJ/NnT2TjXcNnVw21QT2wTj2ie3TY4zK8yInK5gyjFH5wCfNLMLJJ0TX78x0+8pM3vBSgZ2WZ3jVETtwY8K3Axsk3SspE2EIOadPX12AmfG7dOAT1kIWNkJnC5ps6RjgW3ATSPG/CvgJXH7p4F/XM17dhxnNhjDPrltchxnKhQxR4MeY7IDuCxuXwb8wtgjRuq5cuQ4NWScCvRmtizpbOB6QmrbS83sNknnAbeY2U7gEuB9Maj5YcIFBbHf1YRg5mXgDWbWAsiNGU95AfB+Sb8BPAG8bvWzdxyn7qzWPrltchxnWkx55ejZZnZfOI/dNyTz5ZYoK14GLjCzvxo18EzfHKW1InLZoZRIYIrsUErcaEV2qFS6okReQswO1di3nBwTJBTNVJIS29rNRM6yYF37+rbjeVqUkpQn43M7kYoUspIiexN0y0+2LISMUAsZ9+C+dvmhPLkcZCzfXSrlLE8thTGXlsp+tpzU5GjFeaSqoUIVlPshXsmPc06mUmeMrs9hVUOYXQdc19P25mR7D6H2R+7Y84Hzq4wZ2x8Ffm68GTurxW2T26Y1ZUz75LZpjkmvXBdm+pLQqSEVYo6GxmtL+gTwnMxxb1rBNI42s3slfS/wKUlfNrN/GnaA/yc4TgXEeCtHjuM408Ltk+M4daRCKu+h8dpm9rOD9kn6tqTD4qrRYUC2LICZ3Ruf75L0aeCHgRm/OapamT4XAJ2rTJ8GPS9nPITN1DsbPCpK6oY04v5m0lZ4XRdSD3Hslwb/dgcCF+cuGwtP7VOJd7YV63zsXSq/qic3lx7dTbGK/aZGfyTwcnLCIsC5eAbYsy+Ms5yMbUuJd3a5CHpOpl1sp17KzI/yWDWBUg+5MmFx6f7i+xwVFJ0LlF8hfvHhdOG2CXDbNHD/GtomcPvkrJJRq0XFlW3VekiOkzBlWV0RC3lBfP5QbwdJzwCeNLO9kg4FXgT8j1ED1//myHHqwpiyOsdxnKnh9slxnJoxzVTehJuiqyWdBXyDKP2VtB34T2b2OuAHgT9V0K03CDFHtw8asMBvjhynCuaeWcdxaorbJ8dxasg0V47M7CHgxEz7LcREL2b298C/WOnYI2+OJB0FXE4IiGoTgqUurFp8aa3oCoBuBBeaWSlnKAKgLQmEVpQ9WKPUnGQDoPeW32wjtlkicWlG2UR3UHM8X5d0ZUQtjdi51S6/lr0xCHl5uZxjKmNZXAjzXWiW8y5Ok3wkLLfC8UvJOPv2hXFaexPNzVI5r8ZSfE7+sDtBz7lA6EGey6o/2qsJgC6+o9b0C4yMJcdxJo7bJtw2uW3q4PapXsyKffKEDM40qRBzVEuq1DlaBn7LzH4Q+FHgDZKOpyy+tA34ZHztOBsXG/Jw1gO3TY5T4Papbrh9chzC/fegR10Z6SaIOcSLPOKPS7oDOIJQfOnFsdtlwKfJV6adDFUr06ck7kmLkbBd/r9inHZShX458ejG/coERTcST2zhDs0FOKeV56375D3PZYdOmlqgFefTSrymezYnHtbF8Bk0GuVnocIznHxkFgPB28nY7b4Ncw8AABpsSURBVCJN7r4kqHtvsh3P2UjOrfjHrMQZWmynsg6l3nLL/Drn2obR9f2OcJHGvko+k9R7vyrMPbN1w22T2ya3TcUgbp/qRm3s0yhGrRZ5Iob5Y4KribO6crSidy3pGEIKvBupXnzJcWYe4RcfdcZtkzPPuH2qN26fnHllw98cSToA+Evgv5jZd1RRgy3p9cDrAbaw/2rm6Dj1wOUptcRtk+Pg9qmmTMI+HX300dOboONMkSmn8p4alW6OJC0S/rnfb2YfjM1Viy9dDFwMcKAOXhPzXcgUUukCRJlGO5EzxEBZDQpGLrZz9UUSSUoj9msm43QCoLvsYL+8RmlgdhE83ErbouRkKZWcJPVHYrX7VjP5aBtF/ZXkRMXbTseOAdXalwQ67022Y7uSP+wiALqxnEhTOvNO5pCpNZKVsKSsVM6ylli3HMepB26b3DZ12ubVNoHbp5oyKfu0ffv2tf1yVyOrevTR/ratWyczH2dtmWBijimn8p4aIxMyKLg5LgHuMLO3J7uK4kswoPiS42wkZIMfztrjtslxStw+1Qu3T45TrhxtuIQMhGqyrwG+LOnW2PZ7DCi+tCYMq0wPQwOg08DkTgrdVhIwnAi3LW5r31KyP1O5Po7ZzLR1TTuznN71w9UJek6GKYKeN6fe2fKgdvTOWvJNWs47WzS1U+9s9zmg2zvb3Buf9yXzyaTQLTy13Sl0+4Oeu7y3xXeY88i26ymed01/7XDb5Lap3D/HtgncPtWQ+tmnqqxm5cBXiabHrl3l9nHHrcsU2pWSW2eO26gxR2b2GXoEGAl9xZccZ8PiHtha4bbJcRLcPtUKt0+OE6jzCtEgVncr6DjzRkyVO+hRBUknS7pT0i5JfbUtJG2WdFXcf2PMcFTsOze23ynppBWM+ceSnljNW3YcZ0YY0z65bXIcZxoUK0eDHnVltssh5+qLQEfGUtTPgCQAup1qVzK/Ghl5SVo9XkXQcy44OmkrDhn0AStqaJRoacrg4aRflJWkUpFWEqRsi2G7kLDAoIDrSCoviefpkq7kZCpJW3Ofde0LxxdtiVwllam0MjKVYh6WkbOkFJIkW1/NSEiVu3rXrKQm8E7gpcBu4GZJO83s9qTbWcAjZnacpNOBtwGvjIUDTweeBxwOfELSc+MxA8eUtB1wrcN64LbJbdMaMo59ctvkODVnnaR0KY1RNdwGMKuyOl85cpyKjBnwfAKwy8zuMrN9wJWEYoApOwhFAQGuAU6MQb07gCvNbK+Z3Q3siuMNHDNe8Pwh8LvjvGfHcWaDMeyT2ybHcabGRk3IsCHIpdC1Isg4k0IXEudm4r0sfmdS72xnO/Hi5u46sx92V+X26LFNUtoWns+0Cn1jU3lMezE+LyRzLE6e8852nS+O1+UNLreLYOfGvvKgwlPbTDyxhec49cg2WhlPbdpWfOaTqFC/Fli31zzDoZJuSV5fHFOxFhwBfDN5vRt4Yc8YnT5mtizpMeCQ2P65nmOPiNuDxjwb2BnTxQ6duLO+uG3KnS+O57apGuPZJ7dNjuOUdysTTOU9zZUjSa8A3gr8IHCCmd0yoN/JwIWEuhnvMbMLRo09NzdHjjM2w6+LHjSz7UP2j7gcHNpnUHvuOtckHU7IgPTiIfNxHGcjsXr75LbJcZypMOUisF8BfhH400EdKsqG+/CbI8epgo0Xc0T4pzwqeX0kcO+APrslLQAHAQ+PODbX/sPAccCu6JndX9IuM1t/4bLjOJNnPPvktslxnKkwzSKwZnYH9MTZ9tOR+Ma+hcR3Tm6OcgHQq6kvkkhXOjKVXJX6TNBz7usZFNS1kKncXsgiuqQrrf6g53YScNyKQc/WLNuK7fR9dSaXka50BVkn24U8pSvouWhLpSudtqQmS1qlPtZqSSVCxffV9YNe9cc9HWcN646MWUzxZmCbpGOBbxGCmF/V06coDvhZ4DTgU2ZmknYCV0h6OyHoeRtwE+Fb7RvTzG4DntOZt/SEX3ysI26bynm7bZoaY9gnt03O+KRXwFu2rOzYXB2fGtT2mTsmKKcrqLByNCokYVyqyIb72Dg3R44zRUI2qNUfH3X6ZwPXE3Svl5rZbZLOA24xs52Eaurvk7SL4JU9PR57m6SrCZ6OZeANZtYCyI25+lk6jjOLjGOf3DY5jjM9DLOlYR2GhiRI+gSJQyXhTWb2oQoTqCIb7sNvjhynCmbjyuows+uA63ra3pxs72FAtXQzOx84v8qYmT4HrGa+juPMCGPaJ7dNjuNMBwP2jew18Giznx1zAlVkw31szJujQsayivoiXTKWOI61EkkGYX0w/RnK3ZaOkrF0JBCJ5qSUkpQ9O9KVJCNUezGRtiz2DZOVrlhmQsUcBklXCrlMI5GhFNKVZpolKkpWGhm5Std2V32RTEao4vPO1RcZ88ZkItRgCs6M47YpPLttmjw1mYYzp6xUSpeSk82NktJ95CPl9imnrP7czhqwrrXgqsiG+/A6R45TBQsxGIMejuM464bbJ8dxakmbsHI06LF6JP0bSbuBHwP+WtL1sf1wSddBkA0TygdcD9wBXF1F4rsxV44qkqsvMqpKfScQmuVMWz9dbe1+r2Mz8UQWsgi1S1drO3pqG4lHtrUp8d5Gb2mn8jzQ7nhn08DszOQK72w6h9Q726lSnwtw7m/TcjtpSwOgM0HPxXbaVgScryYQOqXjiZ/wRYFfYzhrhNumuMtt0wrGnvyQjrNujErI4KtFa0p7rLWU4UXYVouZXQtcm2m/Fzg1eT1S4tvLXN8cOc5KGDfmyHEcZ1q4fXIcp34YMDQhQy3xmyPHqciYqbwdx3Gmhtsnx3HqhzGtlaNpsrFvjirWF0klDqmMxYraHrmhk21FiYslcpasjCUbzFs2NVr99TW0HOaoTUlwdCIlKSQraSB0IVmxRvr+475UzdIJvE7a0o+n3S9TKfTrzX2JNKUT9JwcnG4Xmvfl8h+kEwidkfN0fUedXZl+KVP2msrcM+tMELdNycnjPrdNq8btk7PupMVshtXLGVUP6dFHw7PXNqoVjVUnVRgvW916sbFvjhxngnhgs+M4dcXtk+M49cNldfVmSArdbIX6ZL8lLk1lqp6XQc8LSdtybMvMIW3vCnpeiKdNPJ+LwSub/vA1lpL0vgvRO7sv8c4uRK9zevJhsXSW3+54Z5f7vbNdQc1Flfml1Pua2R4R9NzxwHZ5Yod7KyznqZ0Ghgc8O9PBbdNg3DZVw+2Ts96kq0XFKlJuBWlUyu+tWyc3J6cGzKasbmT6CUmXSrpf0leStoMlfVzS1+LzM6Y7TcdZb0KRxUEPZ31w++Q44PapfrhtchwoZXWTT+U9Tark5nsvcHJP2znAJ81sG/DJ+NpxNjZmgx/OevFe3D45jtun+vFe3DY5c0+xcjToUU9GyurM7AZJx/Q07wBeHLcvAz4NvHGC81pbEnlEd5X6qP1IPG9FZfcuCUuUw1gSkKjYsSs4uuucMcA5+eEqZBhpBfciYLiZBBEX0pR0u9FM26LkRv2B0F0UTalcJVPbJK0e35GuZAKcu+QqiYylCHbufl9FZfqkX3G+XC2RnOxlJYyQwIw+3jX9dWTD2ye3Tcm83TYNHsPtU93Y8LbJcSoxXzFHzzaz+wDM7D5JzxrUUdLrgdcDbGH/VZ7OcWqAX3vMCpXsk9smZ0Ph9mkWWNW109FHH71G03OcaVDfFaJBTD0hg5ldDFwMcKAOXn/znUuh27W/31Obq1JviSCx45VMPa3FvrRafSbla+rlVSt+HQsZT2SrTJfLcpI6twhwTqrQ0+hv63hqc3l8UzKe2i6PZLvd35bxznalxi22k7bCK9v1mbSLz3FEIPQ6kQt4d2YXt01umzaKbQK3TxuN1D5t3769Hn9kVRmWytuZM2YzlXeVmKMc35Z0GEB8vn9yU3Kc+iEz1Br8qDSGdLKkOyXtktSnNZe0WdJVcf+NqSRD0rmx/U5JJ40aU9L7Y/tXYmDw4lgfwGzh9smZK8a1T26b1gy3Tc6cUcjqBj3qyWpvjnYCZ8btM4EPTWY6jlNjxgh4ltQE3gmcAhwPnCHp+J5uZwGPmNlxwDuAt8VjjwdOB55HCPB9l6TmiDHfD/wA8C+A/YDXjfPWZwy3T878sUr75LZpTXHb5MwZ00vIIOkVkm6T1Ja0fUi/eyR9WdKtkm6pMvbItU9JHyAEEB4qaTfwFuAC4GpJZwHfAF5R5WS1I1dfpGt/DB7uCoQuamAkQca5QOhiX7KdrWbfSGQoGZlGZ8w0YLiZSkBirZFGMp8oWelqK97jqNvhNJY7Jxtp9UtuOnNLZTipTCVTN6Qz5sig5xEV6Yu+XZKjKSgQjK7g71VwArDLzO4CkHQlITj39qTPDuCtcfsa4CJJiu1Xmtle4G5Ju+J4DBrTzK4rBpV0E3DkOJOvKxvWPrlt6sdt02DGs09um6bAhrVNjrMippqQ4SvALwJ/WqHvz5jZg1UHrpKt7owBu06sehLH2QhouAf20B6PxMVRM15wBPDN5PVu4IU9Y3T6mNmypMeAQ2L753qOPSJuDx0zSlZeA/z6sMnPKm6fHCcwhn1y2zQF3DY5DjDFIrBmdgeABjkRx8Cj5hynEjYqTe+DZjZwWZcBzvmKfQa153ztvWO+C7jBzP52yNwcx5lpxrJPbpscx5kSIxMyjHIsT2oSH5NkwJ9WGd9vjmA6WaJGyVgKeUUzrSXS6JuPWrGtmchCFpLsUFE2Ykktkc5ddK5tJXfYOelKpraHMtKVLqlNR35Teg+s2N+VJSu+l5wHdL0zMRnjZqXaDRyVvD4SuHdAn92SFoCDgIdHHDtwTElvAZ4J/N/jTNxZR9w25XHb1M149sltk7N2JDXXKme127On3N6yZbLzmQa7doXn445b33nUhqErR0Mdy5I+ATwns+tNZlY1Zu9FZnZvTJ3/cUlfNbMbhh3gN0eOU5ExiyzeDGyTdCzwLUIQ86t6+hTBup8FTgM+ZWYmaSdwhaS3A4cD24CbCF7b7JiSXgecBJxoNokqk47j1Jkx7JPbJsdxpkSbcWKOzOxnx52Bmd0bn++XdC0hLtJvjlZERU9tNhA6uTk2RS/uoNPEsbsr0/d7LIug6C49eRoo3Iz7m6lruAhwTgOzV+Gd7RycBmH3e2fLAObUO9vq229dHttWf78hY2cDnYGh1eUn+btrdL+/lR4edPpnA9cDTeBSM7tN0nnALWa2E7gEeF8Man6YcEFB7Hc1IUB6GXiDmbUAcmPGU74b+Drw2eiZ/6CZnbfqN+CsP26bMpN12xTGY9X2yW2Ts6bkVotGrSbNwmpRSm7FaG5Xk0bK6qaKpKcBDTN7PG7/a2CkvfGbI8epRLWU3UNHCFmarutpe3OyvYcB2YvM7Hzg/Cpjxnb/33acuWE8++S2yXGc6TGdhAyS/g3wxwSJ7l9LutXMTpJ0OPAeMzsVeDZwbXTELABXmNlHR43tRspxqlKH2ALHcZwcbp8cx6kd00vlbWbXAtdm2u8FTo3bdwHPX+nYfnM0jGG1RrKB0Orbb+kNcyoBiWN2+fpi0LNSaUYMirZ2MnZSf6T4QewER0MpWUnm3Ql6TudYVcaSrd1h/fvToOac1CbXlrxXy/Wr+IM/tfohnRPQLZlxnPXEbVPfvOfWNoHbJ2e2qZqYIWXWkjTUQU63mmQYYzO9VN7TxG+OHKcSVsYiOI7j1Aq3T47j1JH1jTlaLX5zVIVVBUKrax+Qr1yfCw5OzxHbuoKaUw9hPE/XMUUgtEZ4YpOg6KHkPKRpYHaxPSgYuZMaN/0sMsHMnXESz/eosVfKanX57pl16ojbpv62ebNN4PbJ2XhstCQNdWDNVotSfOXIcTY2rul3HKeuuH1yHKd2TC/maJr4zZHjVMGsO7Wv4zhOXXD75DhOLXFZ3XywwkDolFGV64sfN2WkK5Z6BRPJSUfSokb/MV0nV3+/YqmzMSL4OSfXSOUj2arxmeryucr0mQDnbBX6AWSDnadVV3DMVN6OM1XcNsVzzKFtArdPjjOKWUvisGGYvVVtvzlynErYWEVgHcdxpofbJ8dx6kgbXzmaJyoGQqfe0NSTmKtcX/RNPbYqvJhpgHKaljazP5sat+ccXYzyzqZU9YamaXAzwcxDA5xz40zT41oFA1vvOThOFdw2lcyDbQK3T87GY1LJA9LVonlnXVJ5gydkcJyNjHtmHcepK26fHMepHZ6QwXE2LmaeDcpxnHri9slxnFoyh6m8JZ0MXAg0gfeY2QUTmdWsUTEQOidjyVWu76o5UoQuZyrYZ+dAUlfEMoHXyshDxvy7zUtORgVFDwlwHlE3ZFWBzhMIVjbPBjVTuH3CbdOc2CZw+zRLuG2aADmJWK7NEy+UrFudo9mLOapYaa8fSU3gncApwPHAGZKOn9TEHKdWWAx4HvSogKSTJd0paZekczL7N0u6Ku6/UdIxyb5zY/udkk4aNaakY+MYX4tjbhrr/c8Ybp+cuWJM++S2ae1w2+TMF8XN0aDH6pH0h5K+KulLkq6VtHVAv//T3v2GyFHfcRx/f7hTCxabGv9Uc1EjHm2j2FqKtLQUSXwQbTA+UEj/QGgFERQULFXrs9I+yJOmQrUQ1BJEiOEUcohtTDQP+sTU0zwoaao90mrOvycmVhSTXu7bB/O73Ho3e7e7czc7s/t5wXA7c7/5zfx2dj8zv9ndmQXzLU/HnSPgWmA8Io5ExElgJ7CpQH31FzE75P5/et4Q0zFvyCvHdJwe4tT07BCRneGcnp4dTp1qOsQyDLnLbXwuGtctDbPrnTM0tvv0U9fw/OQ9p2XI2y4tLr/FHeJtwLGIuALYBmxN864FNgNXAhuARyQNLFLnVmBbRAwDx1Ld/cT51MjZ1NvZ1GQbtrIOzqbSOZuWwuDg7DA1lQ1505r57LP5QxnKXl7XBdkV65oNhewFroqIq4HXgQfmFuj0ZESRztEq4GjD+ESaZtZzImLhg7DFtbJD3ATsSI9HgPXKvqe0CdgZESci4t/AeKovt840z7pUB6nOmztufD05n6xvFMwnZ1O5nE3WR5bvk6OIeD4iZnrALwFDOcU6OhlRpHOUd43VeaclJd0uaUzS2P84UWBxZt2VdyY994xxvlZ2iKfLpDf8R8DKBeZtNn0lcLwhNPpx57toPjmbrJcUyCdnU7naPnaanJwsYbXMlsPMBRmaDUvm58Cfc6Z3dDKiyK+zJoDVDeNDwNtzC0XEdmA7gKTJfTHyCfBBgeVWyXk0a0uR39h253e1zdtSP6205dJ2KvyYY3v2Te86b4EiX5A01jC+Pb32Z7SyQ2xWptn0vJMbC5XvJ4vmU042vUHvvA+cTdW05NkEhfPJ2VSujo6dNDDQH8dO9dNvbWkznz7aA7s7PnaStA/4Ss58D0bE7lTmQWAKeDKnXEeZU6Rz9DIwLGkN8BbZ945/vNAMEXG+pLGI+HaB5VaG21JNy9GWiNhQsIpWdogzZSYkDQJfAj5cZN686R8AKyQNpjO0uTvfHtdWPkXE+dA774NeaQe4La0omE/OpnL52MltqaQqHjtFxPUL/V/SFmAjsD5yLzHa2smIuTr+Wl0KtruAPcBhYFdEHOq0PrMed3qHmK7OtBkYnVNmFNiSHt8CvJje7KPA5nTFqDXAMPC3ZnWmefanOkh17l7GtlWO88msZc6mEjmbzJaGskvi3wfcFBGfNinWSr7NU+ii5xHxHPBckTrM+kFETEma2SEOAI9HxCFJvwbGImIUeAx4QtI42VnZzWneQ5J2Af8g++j4zog4BZBXZ1rkfcBOSb8BDqa6+4rzyWxxzqbyOZvMlsQfgLOAvdm1XngpIu6QdDHZ/cNubJZvi1Ws/E+hlo+k2+f8FqO23JZq6qW2WLl65bXTK+0At8UMeuu147ZUUy+1pajSO0dmZmZmZmZVVORS3mZmZmZmZj2j1M6RpA2SXpM0Lun+MpddhKTVkvZLOizpkKS70/RzJe2V9K/098vdXtdWpbuYH5T0bBpfI+lAastT6YdrlSdphaQRSf9M2+e7dd4u1h11zSZwPlWZ88mKcjZVi7OpP5TWOZI0ADwM3ACsBX4kaW1Zyy9oCrg3Ir4OfAe4M637/cALETEMvJDG6+JusivlzNgKbEttOQbc1pW1at9DwF8i4mvAN8jaVOftYiWreTaB86nKnE/WMWdTJTmb+kCZnxxdC4xHxJGIOAnsBDaVuPyORcQ7EfFqevwx2YtoFdn670jFdgA3d2cN2yNpCPgh8GgaF7AOGElFatEWSecAPyBd7SgiTkbEcWq6XaxraptN4HyqKueTLQFnU4U4m/pHmZ2jVcDRhvGJNK1WJF0GXAMcAC6MiHcgCwHggu6tWVt+D/wSmE7jK4Hj6f4LUJ9tczkwCfwpfcz9qKSzqe92se7oiWwC51PFOJ+sKGdTtTib+kSZnSPlTKvVpfIkfRF4GrgnIv7b7fXphKSNwPsR8Urj5Jyiddg2g8C3gD9GxDXAJ/Txx8DWsbq+/j/H+VQ5zicrqq6v/c9xNlWOs2kRZXaOJoDVDeNDwNslLr8QSWeQvbmfjIhn0uT3JF2U/n8R8H631q8N3wNukvQfso/o15GdDVkhaeamwHXZNhPAREQcSOMjZG/4Om4X655aZxM4nyrK+WRFOZuqw9nUR8rsHL0MDKcre5xJdoft0RKX37H0vdLHgMMR8buGf40CW9LjLcDustetXRHxQEQMRcRlZNvgxYj4CbAfuCUVq0tb3gWOSvpqmrSe7E7ttdsu1lW1zSZwPlWV88mWgLOpIpxN/aXUm8BKupGspz0APB4Rvy1t4QVI+j7wV+DvzH7X9Fdk353dBVwCvAncGhEfdmUlOyDpOuAXEbFR0uVkZ0POBQ4CP42IE91cv1ZI+ibZjyPPBI4APyPr9Nd2u1j56ppN4HyqMueTFeVsqh5nU+8rtXNkZmZmZmZWVaXeBNbMzMzMzKyq3DkyMzMzMzPDnSMzMzMzMzPAnSMzMzMzMzPAnSMzMzMzMzPAnSMzMzMzMzPAnSMzMzMzMzPAnSMzMzMzMzMA/g9PJOeHT9LWzwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x216 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(15, 3))\n",
    "\n",
    "a0 = axes[0].imshow(imgGalSim2.array, vmin=0, vmax=1.5e-3)\n",
    "fig.colorbar(a0, ax=axes[0])\n",
    "axes[0].set_title(\"GalSim convolution\")\n",
    "\n",
    "a1 = axes[1].imshow(imgAnalytic2.array, vmin=0, vmax=1.5e-3)\n",
    "fig.colorbar(a1, ax=axes[1])\n",
    "axes[1].set_title(\"Analytic convolution\")\n",
    "\n",
    "a2 = axes[2].imshow(\n",
    "    imgGalSim2.array - imgAnalytic2.array, \n",
    "    vmin=-2e-10, \n",
    "    vmax=2e-10, \n",
    "    cmap='seismic'\n",
    ")\n",
    "fig.colorbar(a2, ax=axes[2])\n",
    "axes[2].set_title(\"GalSim - Analytic\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-24T21:03:08.527534Z",
     "start_time": "2019-01-24T21:03:08.524143Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.197057e-08\n",
      "1.1641532182693481e-10\n"
     ]
    }
   ],
   "source": [
    "print(np.max(np.abs(imgGalSim2.array - imgAnalytic2.array))/np.max(imgAnalytic2.array))\n",
    "print(np.max(np.abs(imgGalSim2.array - imgAnalytic2.array))/convGalSim2.flux)\n",
    "# Same residuals"
   ]
  },
  {
   "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
}