{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convergence of Green function calculation\n",
    "We check the convergence with $N_\\text{kpt}$ for the calculation of the vacancy Green function for FCC and HCP structures. In particular, we will look at:\n",
    "\n",
    "1. The $\\mathbf{R}=0$ value,\n",
    "2. The largest $\\mathbf{R}$ value in the calculation of a first neighbor thermodynamic interaction range,\n",
    "3. The difference of the Green function value for (1) and (2),\n",
    "\n",
    "with increasing k-point density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.extend(['../'])\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "%matplotlib inline\n",
    "import onsager.crystal as crystal\n",
    "import onsager.GFcalc as GFcalc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an FCC and HCP lattice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#Lattice:\n",
      "  a1 = [ 0.   0.5  0.5]\n",
      "  a2 = [ 0.5  0.   0.5]\n",
      "  a3 = [ 0.5  0.5  0. ]\n",
      "#Basis:\n",
      "  (fcc) 0.0 = [ 0.  0.  0.]\n",
      "#Lattice:\n",
      "  a1 = [ 0.5       -0.8660254  0.       ]\n",
      "  a2 = [ 0.5        0.8660254  0.       ]\n",
      "  a3 = [ 0.          0.          1.63299316]\n",
      "#Basis:\n",
      "  (hcp) 0.0 = [ 0.33333333  0.66666667  0.25      ]\n",
      "  (hcp) 0.1 = [ 0.66666667  0.33333333  0.75      ]\n"
     ]
    }
   ],
   "source": [
    "a0 = 1.\n",
    "FCC, HCP = crystal.Crystal.FCC(a0, \"fcc\"), crystal.Crystal.HCP(a0, chemistry=\"hcp\")\n",
    "print(FCC)\n",
    "print(HCP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will put together our vectors for consideration:\n",
    "\n",
    "* Maximum $\\mathbf{R}$ for FCC = (400), or $\\mathbf{x}=2\\hat j+2\\hat k$.\n",
    "* Maximum $\\mathbf{R}$ for HCP = (440), or $\\mathbf{x}=4\\hat i$, and (222), or $\\mathbf{x}=2\\hat i + 2\\sqrt{8/3}\\hat k$.\n",
    "\n",
    "and our sitelists and jumpnetworks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "FCCR = np.array([0,2.,2.])\n",
    "HCPR1, HCPR2 = np.array([4.,0.,0.]), np.array([2.,0.,2*np.sqrt(8/3)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "FCCsite, FCCjn = FCC.sitelist(0), FCC.jumpnetwork(0, 0.75)\n",
    "HCPsite, HCPjn = HCP.sitelist(0), HCP.jumpnetwork(0, 1.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use $N_\\text{max}$ parameter, which controls the automated generation of k-points to iterate through successively denser k-point meshes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kpt\tNkpt\tG(0)\tG(R)\tG diff\n",
      "6x6x6\t   216 (16)\t-1.344901582401\t-0.119888361621\t-1.225013220779\n",
      "10x10x10\t  1000 (48)\t-1.344674624975\t-0.084566077531\t-1.260108547444\n",
      "14x14x14\t  2744 (106)\t-1.344663672542\t-0.084541308263\t-1.260122364278\n",
      "18x18x18\t  5832 (200)\t-1.344661890661\t-0.084539383601\t-1.260122507060\n",
      "22x22x22\t 10648 (337)\t-1.344661442418\t-0.084538941204\t-1.260122501213\n",
      "26x26x26\t 17576 (528)\t-1.344661295591\t-0.084538798573\t-1.260122497018\n",
      "30x30x30\t 27000 (778)\t-1.344661238153\t-0.084538742761\t-1.260122495392\n",
      "34x34x34\t 39304 (1095)\t-1.344661212587\t-0.084538717850\t-1.260122494737\n",
      "38x38x38\t 54872 (1491)\t-1.344661200055\t-0.084538705591\t-1.260122494464\n",
      "42x42x42\t 74088 (1971)\t-1.344661193423\t-0.084538699082\t-1.260122494341\n",
      "46x46x46\t 97336 (2545)\t-1.344661189691\t-0.084538695410\t-1.260122494281\n",
      "50x50x50\t 125000 (3218)\t-1.344661187483\t-0.084538693232\t-1.260122494251\n"
     ]
    }
   ],
   "source": [
    "FCCdata = {pmaxerror:[] for pmaxerror in range(-16,0)}\n",
    "print('kpt\\tNkpt\\tG(0)\\tG(R)\\tG diff')\n",
    "for Nmax in range(1,13):\n",
    "    GFFCC = GFcalc.GFCrystalcalc(FCC, 0, FCCsite, FCCjn, Nmax=Nmax)\n",
    "    Nreduce, Nkpt, kpt = GFFCC.Nkpt, np.prod(GFFCC.kptgrid), GFFCC.kptgrid\n",
    "    for pmax in sorted(FCCdata.keys(), reverse=True):\n",
    "        GFFCC.SetRates(np.ones(1), np.zeros(1), np.ones(1)/12, np.zeros(1), 10**(pmax))\n",
    "        g0,gR = GFFCC(0,0,np.zeros(3)), GFFCC(0,0,FCCR)\n",
    "        FCCdata[pmax].append((Nkpt, g0, gR))\n",
    "    Nkpt,g0,gR = FCCdata[-8][-1]  # print the 10^-8 values\n",
    "    print(\"{k[0]}x{k[1]}x{k[2]}\\t\".format(k=kpt) + \n",
    "          \" {:5d} ({})\\t{:.12f}\\t{:.12f}\\t{:.12f}\".format(Nkpt, Nreduce, \n",
    "                                                               g0, gR,g0-gR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kpt\tNkpt\tG(0)\tG(R1)\tG(R2)\tG(R1)-G(0)\tG(R2)-G0\n",
      "6x6x4\t  144 (21)\t-1.367909503563\t-0.192892722514\t-0.131552967388\t-1.175016781049\t-1.236356536175\n",
      "10x10x6\t  600 (56)\t-1.345034474341\t-0.087913619020\t-0.089866654871\t-1.257120855321\t-1.255167819470\n",
      "16x16x8\t 2048 (150)\t-1.344668575390\t-0.084546609595\t-0.088212957806\t-1.260121965795\t-1.256455617584\n",
      "20x20x12\t 4800 (308)\t-1.344662392185\t-0.084539941251\t-0.088166498574\t-1.260122450934\t-1.256495893611\n",
      "26x26x14\t 9464 (560)\t-1.344661615456\t-0.084539088966\t-0.088165768509\t-1.260122526490\t-1.256495846946\n",
      "30x30x16\t14400 (819)\t-1.344661401027\t-0.084538892419\t-0.088165529659\t-1.260122508608\t-1.256495871368\n",
      "36x36x20\t25920 (1397)\t-1.344661260564\t-0.084538764009\t-0.088165374312\t-1.260122496555\t-1.256495886252\n",
      "40x40x22\t35200 (1848)\t-1.344661230214\t-0.084538734661\t-0.088165342770\t-1.260122495553\t-1.256495887444\n",
      "46x46x24\t50784 (2600)\t-1.344661210808\t-0.084538715598\t-0.088165322977\t-1.260122495211\t-1.256495887832\n",
      "50x50x28\t70000 (3510)\t-1.344661197817\t-0.084538703416\t-0.088165309065\t-1.260122494400\t-1.256495888752\n",
      "56x56x30\t94080 (4640)\t-1.344661192649\t-0.084538698279\t-0.088165303871\t-1.260122494370\t-1.256495888778\n",
      "60x60x32\t115200 (5627)\t-1.344661189980\t-0.084538695678\t-0.088165301128\t-1.260122494302\t-1.256495888852\n"
     ]
    }
   ],
   "source": [
    "HCPdata = []\n",
    "print('kpt\\tNkpt\\tG(0)\\tG(R1)\\tG(R2)\\tG(R1)-G(0)\\tG(R2)-G0')\n",
    "for Nmax in range(1,13):\n",
    "    GFHCP = GFcalc.GFCrystalcalc(HCP, 0, HCPsite, HCPjn, Nmax=Nmax)\n",
    "    GFHCP.SetRates(np.ones(1), np.zeros(1), np.ones(2)/12, np.zeros(2), 1e-8)\n",
    "    g0,gR1,gR2 = GFHCP(0,0,np.zeros(3)), GFHCP(0,0,HCPR1), GFHCP(0,0,HCPR2)\n",
    "    Nreduce, Nkpt, kpt = GFHCP.Nkpt, np.prod(GFHCP.kptgrid), GFHCP.kptgrid\n",
    "    HCPdata.append((Nkpt, g0, gR1, gR2))\n",
    "    print(\"{k[0]}x{k[1]}x{k[2]}\\t\".format(k=kpt) + \n",
    "          \"{:5d} ({})\\t{:.12f}\\t{:.12f}\\t{:.12f}\\t{:.12f}\\t{:.12f}\".format(Nkpt, Nreduce,\n",
    "                                                                                  g0, gR1, gR2, \n",
    "                                                                                  g0-gR1, g0-gR2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, look at the behavior of the error with $p_\\text{max}$(error) parameter. The k-point integration error scales as $N_\\text{kpt}^{5/3}$, and we see the $p_\\text{max}$ error is approximately $10^{-8}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pmax\tGinf\talpha (Nkpt^-5/3 prefactor)\n",
      "-1\t-1.3622362852792858\t203.75410596197204\n",
      "-2\t-1.345225052792947\t24.479334937158068\n",
      "-3\t-1.3446883432274557\t3.322356166765774\n",
      "-4\t-1.3446627820660566\t1.186618137589885\n",
      "-5\t-1.344661289561045\t1.0554418717631806\n",
      "-6\t-1.3446611908160067\t1.1378852023509276\n",
      "-7\t-1.3446611836353533\t1.2547347330078438\n",
      "-8\t-1.344661182870995\t1.403893171115624\n",
      "-9\t-1.344661182337601\t1.6081890775249583\n",
      "-10\t-1.3446611814380371\t1.8951451743398244\n",
      "-11\t-1.3446611800056545\t2.291136009713601\n",
      "-12\t-1.344661177927421\t2.8185182806123508\n",
      "-13\t-1.3446611751184294\t3.4945030859983652\n",
      "-14\t-1.3446611715189616\t4.331229854988949\n",
      "-15\t-1.3446611670905142\t5.336529960337354\n",
      "-16\t-1.3446611618106175\t6.514974709092111\n"
     ]
    }
   ],
   "source": [
    "print('pmax\\tGinf\\talpha (Nkpt^-5/3 prefactor)')\n",
    "Ginflist=[]\n",
    "for pmax in sorted(FCCdata.keys(), reverse=True):\n",
    "    data = FCCdata[pmax]\n",
    "    Nk53 = np.array([N**(5/3) for (N,g0,gR) in data])\n",
    "    gval = np.array([g0 for (N,g0,gR) in data])\n",
    "    N10,N5 = np.average(Nk53*Nk53),np.average(Nk53)\n",
    "    g10,g5 = np.average(gval*Nk53*Nk53),np.average(gval*Nk53)\n",
    "    denom = N10-N5**2\n",
    "    Ginf,alpha = (g10-g5*N5)/denom, (g10*N5-g5*N10)/denom\n",
    "    Ginflist.append(Ginf)\n",
    "    print('{}\\t{}\\t{}'.format(pmax, Ginf, alpha))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the error in the Green function for FCC (at 0, maximum R, and difference between those GF). We extract the infinite value by fitting the error to $N_{\\mathrm{kpt}}^{-5/3}$, which empirically matches the numerical error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAE2CAYAAACZXx5aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXdclWUbx78HcE9QRFFxoNyuzD1Qk9LcE+frADOtTFNz\na2VlWmnDXjIXuUf6gpq5M2eucpWzW9NcDPdABRE47x/P4XSAg3CAM5D7+/mcD+fczz1+z/Oc81zc\n67p0er0ehUKhUCgswcneAhQKhUKR/VDGQ6FQKBQWo4yHQqFQKCxGGQ+FQqFQWIwyHgqFQqGwGGU8\nFAqFQmExyngociRCiAQhRDsbtrdLCDEjnXldhBBDbNm+EOKwEGJyVrapeL5xsbcAhSKH0BV4ms68\nfYCPgTl2al+hSBNlPBQKGyClvGdB9iwfEbCwfYUiTZTxUOR4hBANgF+AD6SU/00lzz9AENADqAX8\nCQyTUh41HM8FTAICgVLAEWC0lPJ3w/FdwGEp5TghxIdATeAf4DUgBggFRgIvAQsNZeKBl6WUe83o\nGQaMBdyAdYAzIKWUU4QQi9B+25WBSkBP4IPE9g3lRxvaKwTMBnQZuniKHIua81DkaIQQVYFNwKep\nGQ4TpgBLgNrAX8DPQoiihmOzgIHAEDTjchrYLoTwSKWuDkABoBHwITAU6ATsR3uo3wZKAgfMaO4N\nTAcmAvXQhqN6J8v2HzRj9wpwMFn5AYY23zW0Xx6ok8a5KxRJUMZDkZMpB2wD5kkpP09H/h+klPOk\nlBJ4E+2h3UsIUQTNcIyUUm4zHB8CXAWGpVLXI+AdKeV5KeX3aD2Z+lLKOOA+oJdS3jR8Ts47wGwp\n5UoTLWHJ8kjD8RNSyuhkx94ylA+VUv4FvA7cScf5KxRGlPFQ5GS+BjyBK4kJQoimQogow+uBEGK2\nSf59iW+klLFoD/wagA/ab+mQyXE9Wq+heiptX0lmGB4AudKpuybasFhiW3Gmnw1cfEb5GsBxk/LR\naD0lhSLdqDkPRU5mBdrw0+dCiB+llDeAw8CLJnkemLxP3gtwBuLR5izMzRk4kfo/aLFm0tI77/D0\nGfUmkry3YYreTFvm9CgUqaJ6HoqczBrgGyDc8Bcp5RMp5UWT1y2T/HUT3wgh8gEvAH8Af6M90H2T\n1d8IOJsBXWnFSTiVTIsT2jxMejkJNDQpnxutN6JQpBvV81DkZHRSyjghxFBglxBisZTy52fkf0MI\ncQRtyOcDtP/WQ6SU0UKIb4GvhRCP0FZRjUCbiA7OgK6HQEHDZP5FKeWTZMdnAkuFEMfQhqtGAl6k\nbXQS+RpYIoQ4jjbUNg4ongGdihyM6nkocirGB62Ucg+wEpgthMjzjDILgFFoD+wSQAsp5SPDsYnA\namARcBSoCvhJKRPnHtJ6sJse34HWozkGpNgFL6VcB0wGvjC0lRttRdWzhp5MzzcUzbh9bCgfB+xM\nQ59CkQSdiiSoUKSNYZ/HF1LK2Wlmtr6W5sAlKeVlk7RTwOdSyuX2U6bISahhK4Ui+9EJeEUIMQi4\nC/RFWzW21a6qFDkKZTwUivThSF30D4AiwGa0jYbHgFbJJvcVCquihq0UCoVCYTFqwlyhUCgUFqOM\nh0KhUCgsRhkPhUKhUFiMmjC3MUKIGmjeUl+RUk6ytx5HITteF0fQbCsN1monK+vNiroyU0dGylpa\nxpL81v5uKONhZYQQrdHiKujRNpB5o20yS+5C+7lFCKEDvpJSjjJJ+wA4AVSXUn6Kg10Xc5oN6UWB\niVLK8VhZc3INBjck3YEngKeUck5mNRjikLwB5AVcpJTTzaVZ2k4qdWTo92Bwn9ILzRNxB7Q4Ko8t\n+Q6lVofh2DPvqbmyaP7M0rwXpnWn57om0/Ks6zXEcNwbuCalfGLr35QatrIiQgg3IEBKOQttR7KQ\nUq4HipHSC+pziRDClX+DHCWmtQAwXItcQoimjnRdzGk2oQ/gDkb9VtGcigY/4C9Du9eFELWyQEN3\nYKWU8iugnhCioZm0BhloJ0Udmfg91AdelVKuBQqj7XGx9DuUog6TY2ndU3Nl/UjfvTDWbe6amClj\nzJ/G9WothAgDOhkMh81/U8p4WJdewG+G91OllMeFEOOAC4C3EKKy/aTZBinlXSnlTJJ6p23Cvy7B\njwMthBBjcZDrkopmDLpMd3Vb7V6mouEp8LEQogBatMJ/skCDQPuegubGvTSai3nTtDIZuD/J6y1D\nBn8PUsr9aDFMQHuwHsbC71AqdaTrnqZSNs17kbxuc9ckmeaWyfI/63rtRnN/M9Nw3Oa/KWU8rEsN\ntC9IOzSfSKA5oquL1u29nFrB5xBTF+Al0IYAQHMC6IH2I3H061IdzaNtIja9l1LKX9F2lJ8GHkkp\n72eBhs/RoiOC5or+N7QohcnTLL0/5urNzO8hlxBiFLBISnmdjH2HktcB6b+nScqm814kr9vcNTHV\nXDRZ/mddr+JAVUM4YTJ4PTKFmvOwLs7AfSnlZiFENSFEWynlFsOxFHGpcxBOaHEwwBATQ/4bp9sh\nr4sQwhctRGx+DIbQ1pqFECUNGn4Fpgghfs6sBilljKHupsAOKaUxImGytMT0dLVjrl4hRIZ/D4bd\n818LIUKFEOfJwHfITB160nlPk5U9h9ZzSPVepPJ9MXetwwxpccB5Q/5EUr1eQojGUkq9EOJtIUSr\njFyPzKJ6HtYl3PACLcxnTo6ZYOrK4DqaWw3QxpBv2l6OxQigLdq4tbcQopEdNAwGFkspFwEBQJes\nqNQwSdtUSvnFs9KyoN6s+D1ItHsQSca/Q4l1+GD5PZVAD2AQz74XZr8vz7iupvkrGfKbvV5Ci0E/\n0JAegxZZMjPXI0Mo42FddqKN9QK4oa2EyKmYDlvtQ/vCAzTAJHyrg2HULKVcJKVcCoQCF6SUttKc\nPOJfosv4k8DtLGqjNzBDCOGSOPGaSlpm691BBn4PQoiJQogPDR8rov3Xb9F3yFwdUsrF6bmnqbSv\n4xn34hnfF7PXNZX8qT0/bgEbDO/Lo81x2Pw3pXxbWRnD8rlrQFGTya0cg2FCcTAwHm0sfT5aiNQv\n0L7g9aSUE+ynMCXmNBuWhuYFpgJtgLdNhgVsogEtbsdgDP+NSilXZEE7g4EZaDE9dEBztIiISdKk\nlBbFODdXr5TydEZ+D0KIimhRGfMAdYDhhkPp/g6Zq8Mw7JPmPU2l/SKkcS+S143Wu0j1uprTYu56\nCW0J9ztoiylySynnG9Js+ptSxkOhUCgUFpPtJ8yFEI2BN9HG1EdIKR+kUUShUCgUmeR5mPN4w/Ba\ngIPsTlYoFIrnHYfueRh2u34upXzZMKY3G219dAwwSGrxoZ2klLFCiEiS7hpVKBQKhZVw2J6HYXdk\nMP+uaOgC5JFS+gITga8N6Y+F5numFNpyNYVCoVBYGYc1HsDfQFeTz00xxGiWUibunATNwMxDG7pa\nbkuBCoVCkVNx2GErKeU6IUQ5k6TCwH2Tz/FCCCcp5THgtbTqO3r0qFpWplAoFBZSt27d5HuNNPR6\nvcO+fHx8yvn4+BwwvP/Kx8enu8mxK5bUdeTIEb0jEBQUZG8JDsO0adPsLcEiHOXeWVNHWFiY1dvJ\nynqzoq6M1mGrckFBQUnuS1bWnRaG56bZZ6ojD1slZz/QDsCwdf+kfeUoFApFzsVhh63MsA54VQix\n3/A5zaEqhUKhUFgHhzYeUsrLaK4SkFLqMUTPUigUCoV9yU7DVgqFQqFwEJTxsDENGjSwtwSHoVat\nWvaWYBGOcu9spcNa7WRlvVlRV0brsFU5S/Lb8juaYxwjHj16VF+3bt20MypsRnh4OJ6envaWkSlO\nnjxJbGysvWUoFBkmd+7cvPDCC2aPHT16NNWlug4956FQODqxsbGof0oU2ZmjR49mqJwatlIoFIoc\nzjfffMOKFZaFh1HGQ6FQKHI4JUqU4OZNyyLXKuOhUCgUCotRxkOhUCgUFqOMh0KhUCgsRhkPhUKh\nyMZIKe3SrjIeCoVC4cCcOXOG+Ph4Ll68mGJP0aVLl7hz5w5Pnz5l3bp1bN68mQkTJhAdHW11Xc+N\n8RBCvCyECLa3DoVCochKAgMDadKkCbt27SJ37txJjv322280aNCAEydOsH//ftq1a8ejR484ePCg\n1XU9F5sEhRDeQB3+DVmrUCgUzwWTJ0+mY8eOZo/FxcXh7OxM3bp1EUIAcPv27VR3jGclDms8hBAN\ngc+llC8LIXTAbOBFIAYYJKW8mJhXSnkB+EoIsdQ+ahUKhcI6nDp1isKFCyOl5I033jCm37hxAw8P\nD+Pn2NhYFixYQPfu3XF3d7e6LoccthJCjEWLTZ7Yk+gC5JFS+gITga8N+aYIIVYKIYoa8pkPl6hQ\nKLKU9AyL3L9/n4ULF9pAzfPNhAkTaN68Ofnz52ffvn3G9H379tG0aVPjZzc3N15//XV27drFsWPH\nrK7LUXsefwNdgWWGz02BrQBSyt+EEPUM7ycnK/dML4/h4eFZLFORGaKiotQ9yYbMnz+fDh06ABAZ\nGcns2bMRQnDmzBmGDBlCmTJlAChSpAhVqlRhw4YNqQ675HSePHnCli1bkqTp9Xpy585N+/btWbNm\nDXFxcfTq1Yu8efNy7tw5o8GIjo4mb968KeqsWLEi27Zto06dOunWERUVxePHjy36PTqk8ZBSrhNC\nlDNJKgzcN/kcJ4RwklImJCsX8Kx6s7sH1+eN58GrbkREhL0lWJ34+HicnZ0B+OOPPyhQoIDxvo0d\nO5aRI0dSt25d/vjjD0aPHs3q1auNZX19fRkxYgStW7dOMdmrgDx58tClS5dUj7u5uVGzZk0AwsLC\njC7XHzx4QJEiRYz55s2bx9OnTxk2bBhXr16ldu3aFukoVKgQ0dHRKX6Pz/p+O6TxMMMDoJDJ5xSG\nQ6FQZBy9Xs+WLVs4fPgwrq6u5MuXD4DWrVszZ84cPvvsMwAWLlxofH/p0iWklEavwrVq1eKff/5J\n8U9By5Yt2bRpE127drXxWdmH2NhYgoOD2bRpk/G6JrJx40ZmzZpF+fLl6d69Oy1btnxmXX5+fixb\ntoyCBQtSsmRJGjduDMDevXtp1qyZMV+7du34888/WbNmDa6urvTr1886J2eCQ855mGE/0A5ACNEI\nOGlfOQrF80NERAT9+vXjxo0bfPjhhwwfPpzBgwczcOBAxo4di6urKwB3794FoECBAoC2Oa1kyZJJ\n6vLw8OD48eNJ0urVq8eOHTtSbf/WrVuMGDGCOnXq0LRpU95//30ePnwIaP9tV6lShTlz5tCwYUPe\nfPNNs2mgrTIaP348TZo0oW7duowcOdLo7C+1MtYgd+7cDB48mL59+/LkyRN++eUX47EOHTrg7+/P\n3Llz0zQcADqdjoCAAPz9/enVq5cx/e7du0l6HmXLlqVDhw5069aNDz74AJ3O+tO/2aXnsQ54VQix\n3/D5NXuKUSieF+7cuUP//v0ZNGgQvXv3TnLM2dmZZs2aGYdKfv/9d6pVq2Y8fvv27RRj7vny5Uvh\nnbVUqVJcvXo1VQ3Dhg3D09OTkJAQYmJimDFjBu+++y7Bwf9u29q7dy8hISE8ffrUbFp8fDyBgYG4\nuroSHByMXq9n6tSpDBs2LMkwmrl6rMGRI0do2LAhMTExBAcHJzEUTk6Z+589JiaG/PnzZ1ZipnFY\n4yGlvAz4Gt7rgSH2VaRQPH9MmzYNd3f3FIYjEV9fX1588UUAzp8/T8WKFY3HHjx4QK5cuZLkz5Ur\nF48ePUpRT3x8PDExMSmMzaFDhzh37hzLly/HxUV7HM2YMYPmzZtz4cIFY/7AwEC8vLwArReRPG3X\nrl1cvnyZxYsXU7x4cUCLUdGiRQsOHjxozGda5llcuXKFiRMn8vDhQ15//XU6deqU5Pj169e5fft2\nEmNqyrlz5/D19aVkyZLMnj2bI0eOUK9ePS5evEj58uXTbP9ZHD58OMkqK3uRXYatFIpsRY0aNdDp\ndDZ71ahRw2KN9+7dY+vWrakaDoA6deoYJ8vv3btHoUL/Tj0mDl+ZEh0dTdGiRVOkFypUiKioqBTp\nFy5cIDo6mvr161O7dm1q165NmzZtcHJy4uJF41Yu4wouU0zTLly4gKenp9FwgDaEVrp0ac6fP//M\neswxa9YsZs6cSWhoKH///Tfff/99kuN//PGHcTjPHInhvQsWLEi3bt2M5Q8dOkSjRo3SpSE1mjVr\nlmR/h71w2J6HQpGdOXXqlL0lpMmVK1dISEigevXqSdL37NnDsWPHOHHiBJ6engwcOBBvb29iY2ON\nvQOA4sWL8/jx4yRlo6OjKVasWIq2YmNjzS4rjYuLo3Tp0ixatCjFsWLFihnnWcyVNU3Lk8e8c4mE\nhAQSEv5dW2OunuRERkbi5+dHiRIlABg1ahTBwcHMnDmToUOHcu/ePX755Rdat25ttvydO3dwc3Mz\nfg4ICKBNmzacP3+eqKgoChYsyLZt24iJiaFz585p6nFUVM/DhuzaBfHx9lahUGh4eHig0+lSjP83\nb96cLl26cOjQIYYOHYq3tzcArq6uPHjwwJivTp06REZGGj/Hx8cTHh5O5cqVU7T1+PHjJL2WRLy9\nvblx4wYFChSgbNmylC1bFicnJ6ZNm8bt27dT1Z58Qtjb25vw8PAk8y2RkZFEREQY9aeXu3fvpljq\nOnjwYMqVK0e3bt2My5NTY//+/Ul6F2XKlKFFixbMnTvXuFzZ19eX0NDQdGvaunWrRedgC5TxsCH9\n+kGdOhAXZ28lCoVmPPz9/Vm0aJFxmCWRFStW4OPjk2TJbcmSJZM80D08PChTpoxxt/m+ffuoVq1a\niod1XFxckpVBpjRp0oRKlSoxcuRITp8+zdmzZxkzZgwRERHPHGJKrtfX1xchBKNGjeL06dOcPHmS\n0aNHU6FCBePy1vRStWpVSpUqlSLd39+fDRs2sGTJEkqXLp1q+WvXrqUYVnrttdfYtGkTtWrVArRh\nvPSuiAoPD2f79u0WnIFtUMNWNmT/fqhRA154Af74A1LpaSsUNmPKlCksXryYMWPGUKFCBQoXLoxO\np6Nz584p/vuuX78+CxYsSJI2ffp05syZw7lz5zhz5gxffvllijbOnDljfGgmR6fTMWfOHKZNm0Zg\nYCBOTk40btyYSZMmGR+u5h6y5tJmz57NtGnTCAgIwNnZmZdeeomgoCDjUJu1l6+ePn2aVatWsXPn\nTuLj4xk2bJjxWK1atfDz8zMuPjCla9eutGnThjfffJNVq1axevVqBg0aRFhYGOXKlTPWHRoaSvv2\n7Y17cOyOXq/PEa8jR47oHYFLl/T6QoX0+kqV9PqoKHursS9hYWH2lpBpHOV7ZSsGDBhgcZmgoCD9\nqVOnrKAme9OvXz/9nj179H/99VeS9L59+xrf+/v76+/fv6/v37+/1XQcOXJEv2LFCv3MmTPNHtOn\n8kxVw1Y2plw5OH0abtyA6tXhzh17K1Io0k/btm3Zv39/2hkNxMbGEhkZmWJSXqEt912/fn2SeSRI\nug/Ey8uLf/75x/jZFkGe0osyHnagbFk4eRLc3eGllyAHuEdSPCf4+/uzY8eOFHMOqbFkyRKGDh1q\nZVXZEy8vL6ZOncrcuXOJM5kIjTdZVRMREWHcWxMXF+dQq/iU8bATXl5w+DD85z/QrBlcumRvRQpF\n2ri4uDBkyBDWrVuXZt779+/TtGnTbO/80hqsX7+eixcvcv78efLkycP48eONLllu3brFtm3bmD9/\nPkOGDKFQoUJUqVKFTZs2pbop0R6oCXM7otPBe+9BkSKaAdm2DRzou6FQmMXd3R1/f/808xUpUiTV\nVVY5nc6dOxv3eMyePTvJsRIlSqTYQzJp0iSbaUsv6TIeQog6gCsQLaU8YF1JliGEeAUIBPIBn0gp\ns53TxGHDNAPyyiuwcSPUq2dvRQqFwh5s2LCBq1evcuLECaMrdkclzWErIYQfcF1KuQP4Qwgx3Oqq\nLCOflDIQ+BRoZW8xGaV/f/jqK2jRAnbvtrcahUJhDzp27Mju3bsd3nBA+uY8ikgpwwCklI+Bkmnk\nzzRCiIZCiF2G9zohxBwhxAEhxE4hREXTvFLKTUKI/MA7wBJra7MmdeuCiwt07qz1QBQKhcJRSY/x\ncBJCeAMIIV4FylpTUAbil7sB/wUmSylvWVObtalSBXbu1AxIv36wcqW9FSkUCoV50pzzMISEfUsI\nMQmohvWHhiyKXy6EWAIUBz4TQvwopVybWsXZIV62uzusXJmL3r2LMXSonsuXowgMfJx2wWyIimGu\nUDgGVolhLoTIg/bfvwR+QZuYTulbOYuwNH65Yb4jXWSXJYOentrKq7Zt4bvvigJFmTBBW531PKFi\nmCsUjkFGYpinZ9hqJLBYSjkC+BDonxmRGSBHxi9v0AAOHICDB2HFChg/HtK5L0uhUCisTnqMRy4p\n5SoAw0Pb1uMMOTZ+uRDabvQ9e7TXm28ql+4KhcIxSI/x+E0I0QxACOFO0iEkW7AOeGKIX/4V8K6N\n27c7xYrBL7/A339D374QG2tvRQqFIqeTngnz7UKIckKIlsBVKeVma4tS8ctTUqgQbN4MvXpBly4Q\nGgr589tblUKhyKmka4e54WF+2cpaFGmQOzc8eQLOztC6tbYXRHl/UCgU9iBdjhGFEM5CiEVCiBlC\niOZCCGdDehXrylOY4uQEb78NR45A6dLw8suaa3eFQqGwNekyHlLKeGAq0AzwAxoaDo0RQnhZR5rC\nHJ06wezZmguTevU0l+5Xr9pblUKhcASklDZryxKX7N2AplLKj02cI5YD/hFC/CGE6J318hTm6NYN\n/vtfbdiqUyfNI++5c/ZWpVAorM2ZM2eIj4/n4sWLxCZbOXPp0iXu3LnD06dPWbduHZs3b2bChAlW\nCyBlUTwPQw/ElN+AqsBetJ3hChvRqxd88YU2lDV5Mvj5aXHRFQrF80tgYCBNmjRh165d5M6dO8mx\n3377jQYNGnDixAn2799Pu3btePToEQcPHrSKFkvieRQwk3ZASnkOGC6EeAc4kjWyFOmhb99/3xcu\nrE2ir10LTZrYT5NCobAekydPpmPHjmaPxcXF4ezsTN26dRFCAHD79m1eeOEFq2ixxHgUT56QbNnu\nc+Y8I3vRvbu2nLdLF1i+XDMkCoXi+eLUqVMULlwYKSVvvPGGMf3GjRt4eHgYP8fGxrJgwQK6d++O\nu7u7VbRYMmy1RwjxrA16RTMrRpE5WreGH3+EgABtH4hCYS0sHQq5f/8+CxcutJKanMOECRNo3rw5\n+fPnZ9++fcb0ffv20bRpU+NnNzc3Xn/9dXbt2sWxY8esosWSnkcImgF5LKWcZ3pACKEDRJYqU2QI\nHx9YuhReew3u34fXX7e3IsXzxvz58+nQoQMbN25k2bJl/Pnnn+TOnZtatWoRFxdHdHQ0HTp0YODA\ngegM3jyLFClClSpV2LBhQ6rDLjmdJ0+esGXLliRper2e3Llz0759e9asWUNcXBy9evUib968nDt3\nzmgwoqOjyZs3b4o6K1asyLZt26hTp06W60238ZBS6oUQfYHdQoj+wFzgGJAbGGV4r7AzGzbA55/D\n6tVadMJ792D0aHurUmRn4uPjcXZ2BuCPP/6gQIECeHp64unpSe3atWnRogXu7u4sXboUgHXr1jFx\n4kQuX77MlClTjPX4+voyYsQIWrdunWKyVwF58uShS5cuqR53c3MzRhgMCwujQYMGADx48CBJrPh5\n8+bx9OlThg0bxtWrV6ldu7ZV9FrS80BKedXgnPC/wCK0YS+94f3MrJeXNob46qOBWGCclPKmPXQ4\nCgMHQmQkvPGGNnnep49mQKZMef5cuiuyDr1ez5YtWzh8+DCurq7ky5cPgNatWzNnzhw+++wzABYu\nXGh8nxqJD6vQ0FDGjRtHwYIFjcdatmzJpk2b6Nq1q5XOxLGIjY0lODiYTZs2Ga9xIhs3bmTWrFmU\nL1+e7t2707Jly2fW5efnx7JlyyhYsCAlS5akcePGAOzdu5dmzZoZ87Vr144///yTNWvW4OrqSr9+\n/axybhYZDwDDw7mPEGIYUAm4JqW0Z0SfPGi+r1oDjYGf7KjFIZg0CZ4+1eY+1q7VohLeu6ftDXGy\naHG2IicQERHBmDFjePXVV/nwww+N6fHx8fTp04e6desCcPfuXQAKFDC38PJf4g2un3U6nXHYKpF6\n9eoxbdq0ZxqPW7du8cknn/Drr7+SP39+/Pz8mDBhAgULFiQsLIwWLVowYsQIFi9eTK1atZg8eXKK\ntE8//ZQZM2awb98+YmJiaNasGe+99x7u7u5m65g3b16qejJD7ty5GTx4MEWLFmXBggX88ssvRiPR\noUMHwsPDk0x8PwudTkdAQECK9Lt37ybpeZQtW5ayZbWAr926dcuCszCPxcYjESnlHeD3LNRiRAjR\nEPhcSvmyYT5lNvAiEAMMklJeNNFxUAjRGK330cMaerIjkydr3nf79NGGsvr0gcBAWLgQcuWytzqF\no3Dnzh369+/PoEGD6N076T5fZ2dnmjVrZhwe+f3336lWrVqadW7fvh2dTsfAgQNTGJpSpUpxNQ2X\nCMOGDcPT05OQkBBiYmKYMWMG7777LsHBwcY8e/fuJSQkhKdPnyZJCw0N5dGjRwwYMICiRYsSHByM\nXq9n6tSpDBs2jNWrVz+zDmtw5MgRGjZsSExMDMHBwUl6GE6Z/G8uJiaG/HbykJph42EtDDHM+wMP\nDUnGGOYGo/I10EUIMQWt5zMTbX9JO7QY52Ntr9rx0Olg6lQtqFSZMlpkwu7dtdfq1WBmbk2RA5k2\nbRru7u4pDEcivr6+vPjiiwCcP3+eihUrms1369YtAgICuHbtGuHh4UyYMIEBAwaYzRsfH09MTIzZ\nCd5Dhw5x7tw5li9fjouL9niaMWMGzZs358KFC8YygYGBeHlpnpHCwsKMaWXLlmXXrl1cunSJXbt2\nUby4tsPgm2++oUWLFhw8eNBYzrSOtLhy5QoTJ07k4cOHvP7663Tq1CnJ8evXr3P79m2zxvXcuXP4\n+vpSsmRJZs+ezZEjR6hXrx4XL16kfPny6Wo/NQ4fPpxklZUtccRBjMQY5okkiWEOGGOYSyn7oEUZ\nXIg2DxPPhiT9AAAgAElEQVRiW6mOjU4HnTtrf/Pn15bx5ssH7dpBlNUCCSsAatTQrrutXjVqWK7x\n3r17bN26NVXDAVCnTh3jZPm9e/coVKiQ2XzFixdn6dKl/PTTT7i7uxMcHExkZKTZvIUKFSIqlS/g\nhQsXiI6Opn79+tSuXZvatWvTpk0bnJycuHjROOBAmTJlUpRNTLtw4QKenp5GwwHg4eFB6dKlOX/+\n/DPrSI1Zs2Yxc+ZMQkND+fvvv/n++++THP/jjz9wdXU1W1ZvCAFasGBBunXrZix76NAhGjVqlG4N\n5mjWrFmS/R22xOF6HhmIYb4T2Jmeui0J7v688sUXMHFiEV56KRfLlt3Gzc1+sW2joqKe23ty6pS9\nFaTNlStXSEhIoHr16knS9+zZw7Fjxzhx4gSenp4MHDgQb29vYmNjjb2B1ChYsCC1a9dm+/btbN26\n1WzvIzY21myvA7Rd0qVLl2bRokUpjhUrVsw472KufGJanjx5zNadkJBAQkJCivxpERkZiZ+fHyVK\nlABg1KhRBAcHM3PmTIYOHcq9e/f45ZdfaG1mZ+6dO3dwc3Mzfg4ICKBNmzacP3+eBw8eGBcTbNu2\njZiYGDp37pwuTVlNVFQUjx8/tuj3mCnjIYT4TEo5MTN1pIMsi2GePLh7TmXZMpgwAXr3LsXPP4O9\nLkt4eHi2vycRERH2lpBhPDw80Ol0Kcb8mzdvjpeXF/Pnz2fHjh3Ge+Tq6sqDBw/SrDdXrlzo9Xou\nXzYfAujx48ep9mC8vb25ceMGBQoUMD50w8LC+OSTT3jvvfdSnSMwnZj39vYmPDycmzdvGndXR0ZG\nEhERgbe3d5r6k3P37t0Uy10HDx7M2rVr6datG25ubnz66admy+7fvz9J76JMmTK0aNGCuXPnUsOk\nu+jr68vbb7+dbuOxdetW2rRpY/G5pEahQoWIjo5O8Xt81vc7s8NWtthVnmNjmFuD27fhrbfgww+1\nVVjNmoHJaIAiB+Hh4YG/vz+LFi0yDq0ksmLFCnx8fJI8TEqWLMnt27eT5EteDv4dDjp79iwA06dP\nNx6Li4tLsjIoOU2aNKFSpUqMHDmS06dPc/bsWcaMGUNERMQzh5lMdfj6+iKEYNSoUZw+fZqTJ08y\nevRoKlSoYFzeaglVq1alVKlSKdL9/f3ZsGEDS5YsoXTp0mbLXrt2LcWw0muvvcamTZuoVauWMa1Q\noUIpVqalRnh4ONu3b7fgDKxDZoetbDHmsQ541RDDHOA1G7T53FK0KDx6BF27wvr1WiTCl16CrVsz\nNm6uyN5MmTKFxYsXM2bMGCpUqEDhwoXR6XR07tw5xX/b9evXZ8GCBcbPiTvMdTqdccL8ww8/5I03\n3uDOnTvs3buXnj170rZtW2OZM2fOJHloJken0zFnzhymTZtGYGAgTk5ONG7cmEmTJhkfruYessnT\nZs+ezbRp0wgICMDZ2ZmXXnqJoKAg47Bbeh/UGeX06dOsWrWKnTt3Eh8fz7Bhw4zHatWqhZ+fn3Eh\nQnK6du1KmzZtKFKkCKtXr2bQoEGEhYVRrlw5WrduzcmTJzl9+jShoaG0b9/euCfH5uj1+gy/fHx8\nZmemvC1fR44c0Ss0nj7V63v21OvbtdPrY2L0+hUr9PoSJfT6336zrY6wsDDbNmgFctr3asCAAZkq\nHxQUpD916lQWqXm+6Nevn37Pnj36v/76y5jWt29f43t/f3/9/fv39Xq9Xt+/f/8sa/fIkSP6FStW\n6GfOnGn2mD6VZ6ojrrZSWBkXF83zbt680LOntnx3wQLo0AF2pmvpgSKn0rZtW/bv3592RjPExsYS\nGRmZYoJeoXH9+nXWr1+fZF7JdI7Hy8uLS5cuAf8O01kr0FN6UMYjh5IrF/zwA+j12t8OHSAkBHr3\n1oazFApz+Pv7s2PHDrNzHWmxZMkShg4dagVVzwdeXl5MnTqVuXPnGnfpJ/4FbfK6QoUKgDbsFhcX\nxyk7LutTxiMHkzs3rFmjuTEBaN4cNm/WJtSXLbOvNoVj4uLiwpAhQ1i3bp1F5e7fv0/Tpk2z/eo6\na7F+/XouXrzI+fPnyZMnD+PGjSMqKopbt26xbds25s+fz5AhQ4yr1KpUqcKmTZvStePfWmR2wly5\n2svmJHdVUq8e7NihxQa5dw/eecc+uhSOi7u7O/7+/haVKVKkyDNXWeV0OnfubFymO3v2bGN6iRIl\nzO4fmTRpks20pUZmex5nskSFwqGoVg1+/VVzpDh1qja0pVAobMuGDRu4evUqJ06csLcUs2Sq5yGl\n/DarhCgch9u3oVw5zYC0bg1378KXXyqX7gqFLenYsaNDB85Scx425Oeff+bmTccPNzJ8uPYqWRJ2\n74YDB2DQIDCZu1MoFDkcZTxsyI8//kjlypXp2rUrP/30k9VdQWeU2bPht9+0CISurrB9O1y5Ar16\nwZMn9lanUCgcAWU8bMjs2bO5cuUK7du3Z8aMGZQtW5bRo0dz8qRjeVwpUkRz4b57t+YDq0AB2LgR\nEhKgUydth7pCocjZpMt4CCGchRCLhBAzhBDNhRDOhvQq1pX3/FG4cGEGDRrEvn372Lt3L3nz5qVd\nu3bUq1ePWbNmpfAdZC8Sexxbt2qBpfLkgf/9D0qVglattJVYCoUi55Iu4yGljAemAs0AP6Ch4dAY\nIUT6oqkoUuDj48O0adO4dOkSn376KQcOHMDb25sePXqwadMm4uLi7KqvWDH45Rdwd9dWXLm4aJEI\n69cHPz+4ft2u8hQKhR2xZNiqG9BUSvmxlPKAIa0c8I8Q4g8hROoRZayMEKKEEOKwvdrPLM7OzrRq\n1YqVK1dy6dIlWrZsydSpUylbtizjxo3jzBn7rYh2d9cmzxNXWjk5wcyZmmPFZs0gFa/bCoXiOcei\nOQ9DD8SU34CqwF60CID2YixwyY7tZxlFixblzTff5ODBg+zcuRMnJydatmxJw4YNmTNnjjEYjj3R\n6TSX7kOHagbkr7/srUihUNgaS4xHATNpB6SU56SUwwHLHeWnghCioRBil+G9TggxRwhxQAixUwhR\nMVnet4AVQExWte8oVK1alc8//5wrV67w0UcfsXv3bipUqEDv3r3ZunVrEr839mDECPjkE3j5ZTh2\nzK5SFAqFjbHEeBRPniCl3GzyMUu2kAkhxgLBQGIsyS5AHimlLzAR+NqQb4oQ4gegO/Am0EAI0S0r\nNDgaLi4utG3bltWrV3Px4kVeeuklJk+ejJeXFxMnTkRKaVM9ly796/sqMBC++w7atNE2FSoUipyB\nJTvM9wgh3pVSzkzleFZFFfwb6AokuuZrCmwFkFL+JoSoZ3g/2bSQEGKplHLNsyp+XuJld+nShS5d\nuvDXX38REhLCSy+9RNmyZenZsyedOnWicOHCVm0/IsKZCROKcePGQ/7zn8c0agRBQbnp0sWVb765\nR4sW6dsM8jzHMFcoshPWjmEegmZAHksp55keEELoAGFBXakipVwnhChnklQYuG/yOU4IkSKOuZQy\nIK26nzePnp6enrzyyit8++23bN26lcWLF/Ppp5/Svn17BgwYwCuvvIKzs7MV2tX2gLz8clHc3YsS\nEKC5ci9fHjp3LkZQkLahMC1UDHOFwjHISAzzdBsPKaVeCNEX2C2E6A/MBY4BuYFRhvfW4AFQyORz\nCsOR03FxcaFDhw506NCB27dvs3LlSiZMmMDNmzcJCAggMDCQypUrZ2mblStry3hfeUVbwtunDzRq\npO0NadsW7t+HN97I0iYVCoUDYelqq6tAI+AKsAg4CRwBngCpDWdllv1AOwAhRCNDm4pUKFasGO+8\n8w5Hjx5l48aNPH78mKZNm9KsWTMWLFiQJEpZZqlSBX7+WXNjcvq0llazJuzZA599BjNmZFlTCgfj\n4MGDaea5f/8+CxcutIGa9JMe3aY44jk4Cha7J5FS3pRS9gE80FZYeUkpB1uxN7AOeCKE2A98Bbxr\npXaeO2rWrMnXX3/NtWvXGDNmDBs3bsTLy4uAgAB27txJQkLmb1mNGvDnn2AaWbRSJdi3DxYvhokT\nlUv354358+dTrpw2shwZGcnkyZNZsWIF7733HteuXTPmK1KkCFWqVGHDhg32kpqERN0bN26kV69e\nVKlShZo1axIQEECfPn3o2rUrCxYsSBIl0dHOwaFILbj58/YyBHLP8Vy/fl0/c+ZMfc2aNfXlypXT\nT548WX/hwgWrtHXzpl5ft65e/9Zben18fMrjYWFhVmnXluSE71VcXJzx/fHjx/XLly83fu7Xr5/x\nGhw/flzfs2fPFOWHDx+uf/LkifWFmmCqWa9PqfvatWt6IYT+lVdeMaatXbtWL4TQf/DBBynqs8c5\n2IojR47oV6xYoZ85c6bZY/pUnqnKMWIOo0SJEowcOZI///yTdevWce/ePRo2bIifnx+LFy/m4cOH\nWdZW8eKwcyecPQv9+4ODOhFWoP0TuXnzZj7++GOCgoIIDg4mODiYK1eu8P777xvzLVy4kC5dugBw\n6dIlpJTUrVsXgFq1avHPP/+kWLHTsmVLNm3aZDfNyXWnRu3atQEIDQ1N8Tuw1jlkZ5TxyMHUrl2b\n//73v4SFhTF8+HDWrl1L2bJlee2119izZ0+mhrUS3XIVLgxbtmgT6P7+EB2dReIVWUZERAT9+vXj\nxo0bfPjhhwwfPpzBgwczcOBAxo4di6urK4DRu0GBAtp+YSklJUuWTFKXh4cHx48fT5JWr149duzY\nYRfN5nSnRuKmW51Ohy5Z5LP0nMOtW7cYMWIEderUoWnTprz//vtGIxQWFkaVKlWYM2cODRs25M03\n3zSbdvv2bcaPH0+TJk2oW7cuI0eONMYAMpffnijjoSB37tz4+/vz008/cfbsWWrUqMHQoUOpXLky\nU6ZM4dKlSxbVFxcHDRpoE+cA+fLBunWaIWnbFrJwzl6RSe7cuUP//v3p2LEjAwYMSHLM2dmZZs2a\n4efnB8Dvv/9OtWrVjMdv375N3rx5k5TJly9fioBnpUqV4urVq3bRbE53amzfvh2dTsfAgQNTGJr0\nnMOwYcNwdnYmJCSEefPmcfXqVd59N+kU7d69ewkJCWHcuHFJ0kJDQ3n33XcZMGAA4eHhBAcHs3Tp\nUq5fv86wYcPSrMMeWGQ8hBAbhRBVrSVGYX9KlixpjDGyevVqbty4Qb169WjRogXLli3jUTqCebi4\naGFru3eH/fu1tFy5tF3p1appy3tv3bLyidiZjz7SfIAlf330kXXyp5YvLaZNm4a7uzu9e5v3a+rr\n62scljp//jxeXv860X7w4AG5cuVKkj9XrlxmvyPx8fHExGSNByFLNJvTbcqtW7cICAjglVde4Ztv\nvmH8+PGMHj3abN5nncOhQ4c4d+4cM2bMwNvbm+rVqzNjxgx+/fVXLly4YMwXGBiIl5cX3t7eSdLK\nli1LREQEly5dYubMmVSrVo3q1avzzTffcPr06SSrxMzVYQ8s7Xk0BmKtIUThWOh0OmOMkbCwMN56\n6y1WrVpFmTJljPFI9M9YRvXKK7BiheZ997fftDQnJ82VSatWmkPFq1ezfgOjo/DRR9oqs+SvZxmD\nzOTPiPG4d+8eW7duTfUhDFCnTh3jRtN79+5RqNC/W67MDQNFR0dTtGhKZxOFChUiKirKcpGZ1JxY\nxlS3KcWLF2fp0qX89NNPuLu7ExwcTGRkpNm8zzqHCxcuEB0dTf369alduza1a9emTZs2ODk5cfHi\nRWO+MmXKpCibmHbhwgU8PT0pXvxfT1AeHh6ULl2a8+fPP7MOe2DJDnPQ/EotEUJ8A/wDJBnBllLa\nz3e4wmrkyZOHHj160KNHD8LDw1m+fDmDBw8mLi6OAQMGEBAQQNmyZVOUa9UKFi2Cjh1h82aoV0/7\nL/nTT6FECWjTxh1/fxg7VtszorAtV65cISEhgeqm66yBPXv2cOzYMU6cOIGnpycDBw7E29ub2NhY\nXFz+fWQUL16cx48fJykbHR1NsWLFUrQVGxubYogLYNmyZfz9998p5hhAmxAvXLhwkp6ApZoT2zbV\nbY6CBQtSu3Zttm/fztatW1MMhz3rHADi4uIoXbo0ixYtSnGsWLFixnkXc+UT0/LkyZPiGEBCQkKS\n+cfUNNgaS43HJ4a/vmaO6YHn919JBaC5RBk3bhxjx47l999/Z/HixdSqVYu6desyYMAAunbtSr58\n+Yz527eH4GA4cEAzHomMHAktW15n7dpSNG8Ovr4wfry2S11hGzw8PNDpdDxNtgyuefPmeHl5MX/+\nfHbs2GF0WeHq6ppkk2mdOnWS/JceHx9PeHi4WW8Gjx8/Nvvff//+/a2q2Zzu1MiVKxd6vZ7LqQSp\nSe0cALy9vblx4wYFChTAzc0N0Ca4P/nkE9577z2cnMwP8pgaTW9vb8LDw7l58ybu7u6Ato8mIiLC\n7kNU5rB02KrCM14Vn1FO8Zyh0+mMMUauXbvGwIEDWbp0KaVLlzbGI0kc1urcWQsolRw3Nz2TJ8M/\n/0CLFvCf/0Dz5lovRW0stD4eHh74+/uzaNGiFEOQK1aswMfHJ8lDuGTJkknCJHt4eFCmTBnjePy+\nffuoVq1aigddXFwcRYoUsYtmc7oBs0OuicNBZ8+eBWD69OnpPocmTZpQqVIlRo4cyenTpzl79ixj\nxowhIiLimcNMpjp8fX0RQjBq1ChOnz7NyZMnGT16NBUqVKBx4yyLeJFlWNTzkFJeBhBCtACqoxmf\ns8AOKaV9Y6Yq7Ea+fPno3bs3vXv35tq1ayxbtowBAwag0+kYMGAA/fv3p3Tp0qmWz58fhg2Dt97S\n4qRPnAgTJsC4cZqDxWRzsoosZMqUKSxevJgxY8ZQoUIFChcujE6no3PnzsZ9D4nUr1+fBQsWJEmb\nPn06c+bM4dy5c5w5c4Yvv/wyRRtnzpyhVq1adtFsTvfGjRtZtmwZOp3OOGH+4Ycf8sYbb3Dnzh32\n7t1Lz549adu2bbrPQafTMWfOHKZNm0ZgYCBOTk40btyYSZMmGXsX5obmkqfNnj2badOmERAQgLOz\nMy+99BJBQUHGYTdzddiN1HYPmnv5+PiU9PHxOejj4xPr4+NzzsfH52/D+xM+Pj4lLKnL1q+csBPY\nkUhISNDv379fP3jwYL2rq6u+TZs2+lWrVumjo6ONeVLbYZ6QoNdv3qzXN2+u15crp9cHBen1Dx/a\nRrel5LTv1YABAywuExQUpD916pQV1KSfjOg2xRHOwVrYaof5f4F4oIKU0kdKWQltuOoehiBNtkYI\n8aIQYo8QYpEQork9NChSotPp8PX1Zf78+Vy7do1+/frx/fffU7p0ad5++21CQk5w8WJq48DafpDd\nu2H1ati1CypUgI8/hmSjDwob07ZtW/Ynrr9OB7GxsURGRqaY4LY1luo2xVHOwdGw1Hi0AYZLKcMS\nE6SU19BcsrfLSmEW0ACIAOKA03bSoHgG+fPnp2/fvmzfvp3jx4/j6enJsGFL8fPT8+abszhx4kSq\ny34bNoS1a7UohVevaq7gR4yAVOY0FVbG39+fHTt2PHOZtilLlixh6NChVlaVNpbqNsVRzsHRsNR4\nRKOtqkpOAlm40sqSGObAr8BgYDowNqs0KKyDl5cX77//PpGRXzB6dATLlwfSqtWHVK1alQ8++ICT\nJ0+a/YELAd9/D6dOQZ48UKeO5i/rpHLQb1NcXFwYMmQI69atSzPv/fv3adq0qUME/LJEtymOdA6O\nhqXGYxvwtRDCIzFBCFESzVX6tqwQZGEM85VALTTDdQ+1VDjboNPpGDGiJGvWFCIhYS1vv72e6Oho\nOnToQNWqVZk8ebJZQ+LpqcUJuXBBcwPfqpW2HHjvXrVCy1a4u7vj7++fZr4iRYpQtarjOKRIr25T\nHO0cHAlL93mMBXYCl4UQiQMH5YATQJ8s0mRRDHMhRGPgW7Sd7x8/q2IVL9uxiIqKombNcIKDczNo\nUCU2bBjLu+++y/Hjx9m4cSNt2rQhf/78dOzYkQ4dOiCESLLaJCAAevaEkJD8BAYWxM0tgaFDH9Kq\nVQypLKtXKBRmyEgMc52lY4BCCGegLVAViAHOSil/saiStNsoB/wgpfQVQgQDoVLKbYZjl4CKlgaf\nOnr0qN7U343C/pjGML91S3Phboper+e3334jJCSEkJAQChYsSI8ePejZs2eKycv4eG1uZPp0ePRI\nW+bbty/kzm3dczh69Cjqe6XIzhw9ehQpJTdu3GDkyJEpjtWtW9fs+mCLHSMCPlLKjVLKL6SU32a1\n4TCDimGeA0huOEAb2mrUqBFfffUVly5dYuHChURFRdGmTRuqVavGRx99xGlD/FtnZ+jRAw4fhlmz\n4IcfoGJF+OoryAKXSgqFIhnZwTGiimGuwMnJiUaNGvH1119z+fJlYzz2Nm3aUL16dT7++GPOnDmD\nTqftVv/5Z/jpJ82YVKgA770H16/b+ywUiueH7OAYcR3wqiGGOcBrVmhD4YBcvAheXpqLd1MSd+82\nbtyYL7/8kkOHDhESEkKrVq0oWrSocWirTp2qrFqlTa5/9RVUrartWB8zBhzQVZBCka2waM5DCGFu\nuEgP6AC9lNJhVzupOQ/Hw3TOwxz9+sHjx7ByJaTHkWhCQgKHDh3if//7H6GhoRQtWpSePXvSo0cP\nqlatyvXrEBQE8+ZpvZPx47Ulv5nh5MmTxMaqKAWK7Et8fDx///23xXMelvY8tqCtuEo7IpBCkUkW\nLND2crRvDz/+CKk4NDXi5OSEr68vvr6+fP311xw8eJCQkBBatmyJm5sbPXv2pH//HkyYUIX586FT\nJ603MmGCFn8kI26DXnjhBQC++eYbSpQokYGzdCyioqJS9RyrsB+OeF8sNR6NgJhEB4kKhTXJk0eb\n+B4yROspbN5sfmLdHE5OTjRp0oQmTZrw9ddfc+DAAUJCQmjRogXFihWjR48ebNnSkyNHBMOGQYEC\n2gqtbt20yXdLcXd358aNG5YXdDAeP35MtAo073DY4r4kuoFPL5YOW72Htkw32wWDUsNWjkdaw1aJ\n6PUwaRJs2ABHj2pGJaMkJCRw4MAB49BW8eLF6d69J8WLv8by5aW5eVObEwkMTN9Q2fNGeu+JwrbY\n675k2VJdtGBQvsD/gMPAKZOXWgWlsAo6HXz2mRYDPTOGA7QeSdOmTQkKCuLatWt899133Lp1g6lT\n6/Pw4Yv4+S3lhx8eUqGC1ua9e1lzDgrF84YKBqXINpgJ1ZApnJycaNasGUFBQVy9epVZs74lb97D\nSFmZwoW7s3r1SSpUiGfsWAgLS7s+hSInYZHxkFJeNsx35AfqALcMdVxR8yCK7Exi4J1vv/2Wa9eu\nERw8nKZN55ErVyMWL15J5crRdO9+n717IUFtUVUoLN5hXlgIsRnN9XkI4IEW4+O4EEINlCpsjjV2\njycaklmzZhERcYjQ0NL06fMR27Z9S+vW53Bze8DgwXf480/ljFGRc7F02OorNG+3Zfh3svwd4CHa\nJLpCYTPi4qBBA1ixwnptODs707x5c77/fjr37k1ky5Zw2rQJ4ocfltKgQRglS95k1KjbXLxoPQ0K\nhSNiqfFoD4yVUhpdLxqGq4YBLbNSmEKRFi4uEBqq7dP49lvrt+fs7Iyfnx+rVr3P/fvvsHXreZo0\nWcLcuesR4g5eXmFMnnxTuUFR5AgsNR4FSbY816Qe5QRbYXOqV9eiDAYFwUcf2W4YydnZmZdf9mPt\n2jFERQWyZctJatZcxxdf7KZ06Sh8fC4yY8Z1HjywjR6FwtZY+sDfCkwWQiRuLtQLIdyBL4HtWaos\nnQghqgoh5gohvhNCVLOHBoV9KV8e9u2D9eth1Cjbt+/s7EzLls3ZuHEYDx/6s2HDMcqX38nkycdw\ndY3ihRf+Ys6cSGJibK9NobAWlhqPdwAv4A7aiqtfgCtAYWBE1kpLN0OAMLRzuWQnDQo74+EBu3ZB\nSzsPnjo7O9O2bXN+/nkQjx61Yu3aPyhW7CgjRlygYMH71K//J8uXhxMfb1+dCkVmsXSp7nUpZROg\nMzAcCDK8r286D5JZLIxhXg4tkmAoEJhVGhTZj6JFNT9YjoKzszOdOzdj9+6+REc3YuXK0+TKdZ7X\nXrtJvnw3adr0d3788ZpasaXIlljq2woAKeUuYFcWawGMMcz7o63gApMY5kKIhmhu4bsIIaYAlYGb\nwGO03lAGXNspFNbH2dmZnj196dkT4uLiWL78CLNm3aZ7d3dcXC7TtOlVJk4sR4sWZe0tVaFIFxky\nHlbG0hjmdYFgtCBV455VsYph7lhERUXZ5J4kJOBwMc1btfKiVSsvnj6NY8UKybJl8bz6am7y5j2N\nr+9lhg/3oF69UjbXZat7orAMR7wvDmc8pJTrDDHMEykM3Df5HCeEMIailVIeJZ3DVcrhm2NhC2dv\n4eGa6/W1a7WJdUdk0iQvJk2CJ0/iCAoKZ8GCInTpUpGCBf/m1Vdv8f77Valdu7xNtCjHiI6Jve5L\nREREqscc7P8xs6gY5ooM4+kJAwZAs2ZwxmF9PmvkyePC2LG1+euvJjx4UIhRo/Jw/LgHdeu6UbTo\nbvr1W8/Zs5fsLVOhADLY8zAs1c1FsjkGKeXjrBCVjP1AByBUxTBXZIThw8HVVQv49NNP2q50R6dg\nQRc++qgmH30Et2/HMX16CX74oQQrVxbF1XULHTtGM2FCbapUqWBvqYociqW+rRoJIY4DT9AmtKOS\nvazBOuCJIYb5V8C7VmpH8RzTvz8EB2ursXbssLcayyhWzIUZM6px9Wo1rlwpSN++Ffj55+pUq+ZK\n8eKbeO21Nfz11z/2lqnIYVja8/gGzUh0QRtOsgoGlye+hvd6tL0cCkWm6NhRc2dy5469lWScMmVc\nCAqqQlAQXL0ax/Tp3qxZ48ySJUVxc9tIp04xjB9fFyFUj0RhXSw1HjWAxlJKNXSkyJY0b25vBVlH\n2bIuzJpVhVmz4PLlOGbMqMyaNc4sXlwEN7ef6Nz5CePH18PHRxkSRdZj6YT5X4Dt1w8qFIpnUq6c\nC5TZvDAAABplSURBVN99J4iMrMTFi4Xp2VOwaVMdqlQpjLv7jwwa9D/OnVNDW4qsw9Kex7fAfCHE\nLLT9GLGmB6WUm7NKmEKhyBjly7swe7Zg9my4eDGOGTOqsW6dCwsXFqJYsXV06RLLuHENqFxZ9UgU\nGcfSnsciNN9WM4C1wEaT14aslaZQ2IZDh+Cdd3gu/U1VrOjC3Lk+XL9ekXPniuLvX5316xsgREFK\nlFjDG2+s4vx51SNRWI5FPQ8pZXbYF6JQWET16toekP/8B5Yvh9y57a3IOlSq5MK8eT7Mmwfnz8cx\nY0ZN1q3Lzfff56d48VC6dn1KQEAFtUlQkS4yZAyEEC2EEMOFECOFEK1NXLQrFNmOQoVg0yZ4+lTb\njf7okb0VWZ/KlV0IDq7MzZvl+OsvVzp1epG1a31p1qwCJUr8j7feWsnff6seiSJ1LN3nUVIIcQjY\nghY9cBjacNUxIUQJK+hTKGxC3rwQEgKlSsGrr8Ldu/ZWZDt8fFz4/nvNkOzeHU/HjnUIDW2KEPkM\nhuQHZUgUKbC05/FfIA6oIKX0kVJWAioC99C83SoU2RYXF1iwQHNl8uef9lZjHypVggULKnHrlhcn\nTxY3GJIm+PgUoESJEN58U82RKDQsNR5tgOFSyrDEBCnlNWAU0C4rhSkU9sDJCaZPBz8/eyuxP9Wq\nuRgNyZkzbnTqVJu1axsjREHc3dcwaNAqzp27aG+ZCjthqfGIBsyFrkkAnDMvR6FQOCJVqrjw/feV\njHMkXbvW5KefGlGlSiHc3dfy+uurkVIZkpyEpcZjG/C1EMIjMUEIURLN59S2rBSWXoQQI4QQi4QQ\n+4QQb9lDg0KRk/DxcWH+/MrcuFGec+dc6dbtBTZubEDVqkUoXvxHBg78nzIkOQBLjcdYwB24LISQ\nQgiJFje8AHaKYS6l/C/wBnBKSjnXHhoUzz+bN8PWrfZW4XhUquTC3LmVuX69AufPF6FHj+ps3lyP\nqlWLUrz4egYM+B9nz16wt0yFFbA0hvkN4EWgO/A9MAvoIKVsIKVMPWqIhVgYwxzgP2ibFhUKq1C0\nqBYXZNQoiLKW/+hsjre3C3PmVCYysiIXLhSmV6+qbN1al+rV3ShWbCOBgSGcOaMMyfOCxfs8pJTx\nUsqNUsovpJTfSil/yUpBhhjmwUAeQ5IxhjkwEcOqLiHEFCHESiGEK9BMSvlzVupQKEzx9YWTJ7Ul\nvNWqad559eZm/xQAVKjgwnff+RAZ6c0//xSmTx8ffv65NjVqFMPNbRP9+4dy6pQyJNmZNDf3CSHa\nAdullE8N71Mli3xbWRTD3KBRTdYrrI67OyxaBL/+CkOGwO7dMGuWvVU5PuXKOfPttz58+y1cvRrP\njBmVWbsWVqwoTpEiW2nb9hETJrxIzZqV7C1VYQHp2Rm+ESgJ3DC8Tw09WbDiytIY5oYyA9NTt6MF\nkM/pREVFZct74u2t7UiPiHAmPPz5cohl7Xvi7AwTJxZk4kSIjHzIvHnubNhQhlWrSlCw4Db8/G4z\ndGh5XnihvNU0ZEcc8beSpvEw9WdlJ99WWRbDXPnscSzCw8Oz9T0pVy7tPNkNW94TT0+YN6808+bB\n9evxfPVVeVavLkPbtl4UKnSYV199wIQJ1alXr7JN9Dgy9vqtRESkPpVtqXuSnUKIombS3YUQRzOg\nLT3sx7ABUcUwVzg6Dx9CZKS9VWQ/PDycmTFDcPlydW7cyM/QoZ4cPepNgwYlKVx4D127/sShQ+fs\nLVNhQnrmPPyAav9v796jo6rOPo5/JyEBEkAEAQWEKsgDomhJioIXBKVSQEVQihdkKSAX74KKWrVq\ny+ulgKgLL2hVrL6Ar0WhXla1oLZgUaOCNzaCUkHiDUQCCBTJ+8cedIgJ5ISZOXP5fdaaNWfOTCZP\nuxd53Gef/TzRl92BkWZW8X6TDkCb+Ib2o9lAr2gPc4DzE/R7RPbaq6/6u7JuvhlGjvSXaSSY/fbL\nZcKE9kyYAN9++wOTJzfjySf3oVu3/SksXEDPnt9yzTVGt26akYQpUr6HW0bM7HDgWSACtAZWA7EX\nesuBjcAU59zDCYpzr5WUlJQXFRWFHYbESPfLVlV5/30YMwa2bIH774fOncOOqPpSeUw2bNjBlCmO\n6dO/Z8WKthQUfMQJJ6xl3Li2nHBCu7DDS6iwxqWkpISioqJIZe9VZ83jPXzxQ6J7LwY457Ko5qhI\nMIcd5mcgjz0GffrAoEEwebJmIXurQYMcbrihAzfcAJs27WDKlPpMn55Lz55NqVPnHY477iuuvPIX\n/PrX7YhEKv17J3EUdJNgj8oSh5nlm1nX+IUlkt4iEX/56oMP4PDDlTjirbAwh+uuO5SlS4vZtKkB\nt95amzVr9qFPn/2oW/dDevZ8mTlzHHu6siI1F6iJk5l1AR4ADuPniac86PeJZLrGjWHEiLCjyGx1\n6+YwduyhjB0L27btYNq0L3nooboMGNCQ3NyP6dLlcy6++ADOPNPIydGMJF6C3np7N7AWGIyvsHsu\nftf3RnyJEBGpph01uuFcdic/P4eLLjqUd945hq1bmzJ16g9s2ZLDuecWUqfOZxx99Gs89thSduzQ\njGRvBU0eRwBXOueeBt4GvnbO3YEvinh5vIMTyVT/+Y+/nPXSS2FHkrlycyMMG9aBN9/sztatLXnk\nkU3k5GxnxIha5Od/QXHxAh54YCnbtyuR1ETQ5LEdv2kPYBnQKXo8H+gYr6BEMl3r1nDHHXDhhTB4\nMKTY5uGMk5MT4ZxzDmXhwp5s3dqGmTPXU1i4mUsv3UHt2mvp1GkRkyY5tmxRIqmuoMljETDGzHKA\nxcBvoucPA7bFMzCRTNe3r19Qb9MGjjgC7r4bfsisaicpKRKJMHBgB159tRdbtnRg7tyvadZsHb/7\nXRkFBRto1+4dbrnFsX69rivuTtDkMR6/Se9KfOHC9mb2CTALeDLOsYlkvIIC+OMf4bXX/CWsL78M\nO6LsEolE6NOnAy+99Bs2bSpi/vxS2rdfzZ13fkWjRpv4xS8+4Oqrl1FaqqxeUdDk8XvgRGC6c249\n8CvgT8AFwBXxDU0ke3ToAHPn+npPEo5IJEL37u2ZM+cUysqO4403SjnqqGVMm/YpLVpspHnzjxkz\n5mOWL98edqgpIeittV2BjdGmUDjnvgCmxj0qEZGQFRe3Y+ZMv3P9gw9WcOed7zJrVh7339+QRo2+\np2/fbVx6aWs6d84jG/ckBk0ek4DpZnYX8Cn+dt0fOec+jFdgIuIbTl1zDVxwAbRvH3Y02atjxzY8\n+qgv3/fxx58yadIinnmmnL/8pRYFBYX06FHGmDEtOfHEfPLyQg42SYImj1ujz91izpXj617FpZ9H\nUNHmUOfj+378yTm3ONkxiCTKjh3QogUceyyMGgXXXw9164YdVXY75JCDuO++g7jvPvjss1VMnfp3\nZsz4nn79fklurtGlyzpGjGhK//51aNAg7GgTJ+iax0GVPA6OeQ5DEb6qbwtgVUgxiCREbi5cdhks\nXgzLlkHHjjBuHCxcGHZkAtCq1YHcdts5rFw5nNWrW3DzzbP55pvZDBv2Go0bb6ZTp1ImTvyeVRn4\nlynQzMM5959EBRLLzI4CbnPO9TCzCH5d5QhgCzDcOfdJzMffBh4GegL9gOnJiFEkmVq0gFmz4PXX\nffvbtWsr/9zGjVBYSFZegw/b/vvvz/jxQxg/HtauXcvMmU/z0EOrGD/+YK69tg/Nm29j8OACBg0q\n4Je/TP8xClrbagf+8lRltgGfA/8L/N45V6N728zsKmAIvuQJQH+gtnOuWzSpTAL6m9ktwCHAvkBf\n4Bt+6jsikpG6dvWPqlx7LTzxBBQXQ1GRfy4uhlat0v+PVTpp3LgxY8YMYcwY+O6773j22b8xbdpH\nTJq0P/fccxp5efvSt2+EM88s4KSToF69sCMObo/9PGKZ2Qj8usdNwOvR012Am4FpwPv423mfdc5d\nX5OAzOx0YAnweDRhTAQWOedmRd9f7ZxrGfP5U/G1trYCVznnvqnse0tKSsoPOOCAmoQkCVJWVkb9\n+vX3/EEJ5KuvcliyJI8lS/JYvDifJUvymDx5PSecsHWPP6sxSazNmzfzyiuvMGvWu/zznw3Iy+vP\n1q1HcOSR33PKKRFOPHELrVv//L+7wxqX0tLSmvfzqGAcMMw591zMuSVmtgrfDKq9mZUCM4AaJQ/n\n3Gwzi+0O3QD4Lub1djP7sY+5c24OMKc6352qTW6yVSo3HkpnzZvDkUfueq68vHGlM4/LLoOGDX+a\npdSvrzFJtLZt2zJ8OGzbto158+YxY8ZVzJ69keXLT+X220+iSZM8BgzIp1+/CMccA3l5GdDDHL8o\n/Ukl51cDrWKOGwX83t3ZAMSm3B8Th4hUT1WXrHr18iVRpk6FTp2gc+dmnHIKlFVsNC1xl5+fT+/e\nvXn00btZt24as2c34oILbmDTpoE88si9nHXWaho33s4ZZ5TzxBMFfPZZ2BHvKujM41/AHWY21Dm3\nDsDMGgO3ATt7jA/AF02MlwX4hfD/M7Ojgffi+N0iWa1fP/8Av6fkzTe/4fPPm6XlNfh0lpuby/HH\nH8/xxx/PXXeVU1JSwl//OpVZs15j3rwi3n57MBMmFNGsWR69e0c4+WTo3t2XtwlL0OQxAngB+NzM\nVuJnLq2Aj4ABZtYbmACcEccYZwO9zGxncjo/jt8tIlGRCLRs+QNdulT+/rJl8MILcNZZ0LRpcmPL\nJpFIhOLiYoqLi5kwAZYuXcr06dP5xz+u5KOP6vDKK6P5+997sGrVfnTt6hPJySf727iTeVNEoAVz\nADPLxde3Ohz4L/C+c25e9L0mQHlVi9ZhKikpKS8qKgo7DImhNY/Us7sxcc4XcZwzB447DoYMgVNP\nhTp1khxkFto5LmvWrGHOnDk888wzLFjwHmajqF9/IJ980g6oxeLFfg0rXkpKSqpcMA+65kH0FtxV\n+LWPh4FPo3sxcM59nYqJQ0T2nhlMnw6rV8OgQTBtml+cnzs37MiyR/PmzRk1ahQvvvgiq1d/yLhx\nh9C06c2sX78fBxxwGpHIhj1/SZwE3efRAH8nVW9gB9AOuAtoZWZ9nHNqaSOS4erV87OOIUN8Iqld\nO+yIstM+++zD4MGDGTx4MFu3buWtt96iXr3CpP3+oDOPiUBtoCU/FUW8BL+h7644xiUiaaBlS2jS\n5Ofny8v9LGXduuTHlI1q167NMcccQ25u8soLBk0effEb8X6cYURLllwMnBTPwEQkfW3cCM8/Dwcf\nDGecAc8+C9vUazSjBE0e9ahQhj3mewKvn4hIZqpfH2bMgJUroXdvmDjR1+eaODHsyCRegv7BfxG4\n0cx2rpWUR++w+hPwUlwjE5G017AhDB/u2+wuWgTduu35ZyQ9BE0elwCtgXVAAfAy8BmwD3B5fEMT\nkUxy8MFVF3XU2kj6CVqS/Uugm5n1wFewrYXfIPiScy7YhhERkajevSE/H0aPhoEDtXckHQS9VXce\nMMA5Nx+YH3O+iZm96JzTLjwRCWzhQr9f5P774YorYOhQuPBCOOSQsCOTquwxeZjZCfzUJ6M7MNLM\nKpZN6wC0iW9oIpItatWC00/3j+XL/QbE666Dp54KOzKpSnVmHmvxpdgj0cdFQGzB+XL8Po+xcY+u\nGszst/jCid8Av3PObQojDhGJj7Zt4fbbw45C9mSPycM59x7R/uRmNh9/2erbRAcWwCnAecCRwFB8\ny1oRyVAPPAAHHuiLASZxT5xUEHTBvAdA9FbdPPxMJPb9zfEIKmAP83uBh/B3fdWo9a2IpI/69eGm\nm+Cii2DkSL820iieHYSkWgLdqmtmR5nZO/iWrxuBsuhj5/Fei/Ywn4YvgwIxPcyBa/E9zDGzW8zs\nSWB/YDjwT3wCEZEMdvbZ8Oabfj1k6VJo0wYuvdSXRJHkCdrPYwo+SfTHd/hLhOXA6cDj0dfH4jcn\n4pxbZGbF0eMbAaK3DT8KbAdG7u6L16xR3cZUUlZWpjFJMek0Js2bw4QJcPnlOfz73/mUlm4JO6SE\nScVxCZo8DgO6RtdBEqIGPcx3uW14d9Q7IrWon0fqSccxqaxne6bJhB7mS4ED9iqa4NTDXERqZPRo\nuPFG+OKLsCPJPEFnHvcAD5rZvfjLS7vUyXTOPR+vwGKoh7mI1Mhll8GUKdChA5x2Glx+eebPUpIl\naPJ4JPp8RyXvlQOJuHFOPcxFpEbat4f77vPtcx98EPr1g86dfYn4ZPb7zkRBb9VNStn1aI+QbtHj\ncmB0Mn6viGSmRo1g/HgYOxaWLFHiiIfqlCfpgy98+N/ocVXKnXMvxC80EZH4ysuDoioq8O3YATnq\nSlRt1Zl5/A2/l+Kr6HFVEnXZSkQk4U47zd+5NXYstGsXdjSprzrlSXIqOxYRySQPPwz33gvHHusf\nV11Vdf8RUetYEREAmjaFW26BTz+Fnj3hnHNg0KCwo0pdQe+2EhHJaIWFcPHFMGoUrFgRdjSpSzMP\nEZFK1KoFZpW/t0PblJU8RESCKC/3ayJjx0KKlZtKKiUPEZEAIhGYOdPPPg47zJeG/ywL63kreYiI\nBHTggTB5si8JX6+eL3nyhz+EHVVyKXmIiNRQ06a+Ze7HH/vSJ9kkLe+2ivbwONs5N8LMuuL7eJQD\nlznnEtVnRESkUo0b+0c2SbuZh5m1ATrzU6fBC6OPh4HBYcUlIlLR5s0wbBi8+27YkcRfSiSPaHvb\n+dHjiJndZ2YLzWyemR0c+1nn3Arn3MSYU7nOuW3AF/gyKiIiKSESgY4doU8fX/7krbfCjih+Qk8e\nQXuWm1nDCl+xyczy8U2q1PJFRFJG3bpw5ZV+s+FJJ0H//jBwICxfHnZkey/05MFPPct32qVnOfBj\nz3Ln3NnOufUVfn4a8AD+0tVfEh+uiEgwdevCJZf4hfUuXeDzz8OOaO+FvmAetGd5zM+dF31+m2o2\niEq1BvLZrqysTGOSYjQmiTdkiH8O8n9zKo5L6MmjEgnrWR5GA3mp2po1azQmKUZjEp7NmyE/35dF\nqSiscSktLa3yvVS4bFXRAqAPgHqWi0i2eOQR6NQJ5s71JVBSXSomj9nA1mjP8onAFSHHIyKScGPG\nwJ13+na5PXrAG2+EHdHupcRlK/UsF5FsF4lA375w8snw2GMwYIAvwPj442FHVrlUnHmIiGStWrX8\nxkLn4Le/9X3XU5GSh4hICioshNNP3/PnwqLkISKSZlJhQV3JQ0QkjTgHxcXhlzpR8hARSSPt2sEV\nV/jF9XHjYNOmcOJQ8hARSSORCJx7Lrz/PpSW+r0hL7+c/DiUPERE0lCTJvDEE3D33X4GsiHJnYyU\nPERE0ljfvvDOO9CgQXJ/r5KHiEiai0SS/zuVPEREJDAlDxERCSxtk4eZ9TCzaVW9FhGRxEnL5GFm\nbYDORFvXVnwtIiKJlRJVdQHM7CjgNudcDzOLAFOBI4AtwHDn3Cc7P+ucWwFMNLPplb0WEZHESomZ\nh5ldhe9FvnPm0B+o7ZzrBlwLTIp+7hYze9LMGkY/V/EegxDuORARyT6pMvNYDpwO7KxcfyzwIoBz\nbpGZFUePb6zwcxXLg+22XFhJScneRypxtbs2lxIOjUlqSrVxSYnk4ZybbWatY041AL6Leb3dzH7W\ny9w5d97uXscqKirSrEREJE5S4rJVJTYA9WNe/yxxiIhIeFI1eSwA+gCY2dHAe+GGIyIisVLislUl\nZgO9zGxB9PX5YQYjIiK7ipSnQksqERFJK6k680goM+sJDAXqArc653RZLGRm1hkYC2wDrnbOfR1y\nSBJlZk2B55xzvwo7FgEzOwK4G/gEeNQ592oYcaTqmkei1XXODQUmAL8OOxgB/B6f0cDzQNeQY5Fd\nXQWsDDsI+VEXoBTYDnwQVhAZN/Oozk5159xzZlYAXAJcE2K4WaGaY/K6mXXFzz7ODDHcrFGdcTGz\nUcAT+HGRBKtmpY1/ATOAZvjEHsrfsIyaeQTYqd4YmALc6Jz7JoxYs0WAMSkG3sLfZXdpCKFmleqO\nC9ALGAl0MbOBSQ80iwQYkyOBXGB99DkUGZU8+Gmn+k677FQHiqLnJwHNgf8xswFJjTD7VHdMGgB/\nxif1p5IZYJba07jsrOow0Dk3GljknHs66VFml+r+W1kJ3APcjl/7CEVGXbaqxk71H6I71YcmObSs\nFWBM5gHzkhtd9gpa1WF31RskPgL8W3kdeD250f1cps08KtJO9dSjMUlNGpfUk9JjkunJQzvVU4/G\nJDVpXFJPSo9JRl22qoR2qqcejUlq0riknpQeE+0wFxGRwDL9spWIiCSAkoeIiASm5CEiIoEpeYiI\nSGBKHiIiEpiSh4iIBKbkISIigSl5iIhIYEoeIiISWKaXJxEJhZnVBq4GzgIizrkOMe+dBdwELAP+\n7Jx7JpwoRWpO5UlEEsTM8oHh+CRyeWySMLPxzrnbQgtOZC9p5iGSOMcDrwAF+FahsTOMH8IISCRe\ntOYhkjiHO+c+BB4EDjWzYwHMzIBlZjbQzM4NNUKRGlLyEEmcCIBzbgO+xe7V0fM98V0TXwaGVffL\nzOyMeAcoUlNKHiIJYGZNgK9iTk0BTjazjsC+zrky59x3QLUWHc2sFbv2txYJldY8RBKjFzE92Z1z\nK83sWeB64K2KHzazt4GngHXASOAOoDWw3Dn3NPAroMjMLgBmOOc2J/5/gkjVlDxEEuMg59yTFc5N\nAhYC98Sci5hZb2Coc+49ADM72zk3I3r8ppn9wzn3tJld5Jz7c1KiF9kDXbYSiSMz62xmDwKXmNmN\nse855/4NPAcsijndEjgP2DfmXOydWCuAdtHjSPR3FMQ7bpGgNPMQiSPn3NvAhdFHZe+fUuHUcvxe\nkL+a2QLn3A/s+u+yFeCix+VmVgsoBl6La+AiASl5iIQkeptue6AjsAV43MxGA/ub2UCgLfCH6MI6\nwGJgMLvuFxEJhXaYi6QYM5vvnOsRdhwiu6M1D5EUYmZnA23MrEvYsYjsjmYeIiISmGYeIiISmJKH\niIgEpuQhIiKBKXmIiEhgSh4iIhKYkoeIiASm5CEiIoEpeYiISGD/Dxm9jFQOmDUTAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111cad208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the errors from pmax = 10^-8\n",
    "data = FCCdata[-8]\n",
    "Nk = np.array([N for (N,g0,gR) in data])\n",
    "g0val = np.array([g0 for (N,g0,gR) in data])\n",
    "gRval = np.array([gR for (N,g0,gR) in data])\n",
    "\n",
    "gplot = []\n",
    "Nk53 = np.array([N**(5/3) for (N,g0,gR) in data])\n",
    "for gdata, start in zip((g0val, gRval, g0val-gRval), (0,1,2)):\n",
    "    N10,N5 = np.average(Nk53[start:]*Nk53[start:]),np.average(Nk53[start:])\n",
    "    denom = N10-N5**2\n",
    "    g10 = np.average(gdata[start:]*Nk53[start:]*Nk53[start:])\n",
    "    g5 = np.average(gdata[start:]*Nk53[start:])\n",
    "    Ginf,alpha = (g10-g5*N5)/denom, (g10*N5-g5*N10)/denom\n",
    "    gplot.append(np.abs(gdata-Ginf))\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.plot(Nk, gplot[0], 'k', label='$G(\\mathbf{0})$ error $\\sim N_{\\mathrm{kpt}}^{-5/3}$')\n",
    "ax1.plot(Nk, gplot[1], 'b', label='$G(\\mathbf{R})$ error $\\sim N_{\\mathrm{kpt}}^{-5/3}$')\n",
    "ax1.plot(Nk, gplot[2], 'b--', label='$G(\\mathbf{0})-G(\\mathbf{R})$ error')\n",
    "ax1.set_xlim((1e2,2e5))\n",
    "ax1.set_ylim((1e-11,1))\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_yscale('log')\n",
    "ax1.set_xlabel('$N_{\\mathrm{kpt}}$', fontsize='x-large')\n",
    "ax1.set_ylabel('integration error $G-G^\\infty$', fontsize='x-large')\n",
    "ax1.legend(bbox_to_anchor=(0.6,0.6,0.4,0.4), ncol=1, \n",
    "           shadow=True, frameon=True, fontsize='x-large')\n",
    "ax2 = ax1.twiny()\n",
    "ax2.set_xscale('log')\n",
    "ax2.set_xlim(ax1.get_xlim())\n",
    "ax2.set_xticks([n for n in Nk])\n",
    "ax2.set_xticklabels([\"${:.0f}^3$\".format(n**(1/3)) for n in Nk])\n",
    "ax2.set_xlabel('k-point grid', fontsize='x-large')\n",
    "ax2.grid(False)\n",
    "ax2.tick_params(axis='x', top='on', direction='in', length=6)\n",
    "plt.show()\n",
    "# plt.savefig('FCC-GFerror.pdf', transparent=True, format='pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the error in Green function for HCP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEyCAYAAAACzlxMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8TFf7wL8TsYtUat+bVE5iSRGUUKqU1lJbrLF1QdVS\n9b6lqrWVtpbyI6WUqhZ9kaitgr62otaE2B3EkqBvkVpDaeT+/jiTMUkmyYxkkuB8P5/7yZ17z/Lc\nm5n73HPOs5gMw0Cj0Wg0GkdwyW4BNBqNRvP4oZWHRqPRaBxGKw+NRqPROIxWHhqNRqNxGK08NBqN\nRuMwWnloNBqNxmG08tA4DSFEghCiRRb2t0UIMcnOsq5CiP4Otr9PCDHKvD9aCLHXznp23QchxPdC\niGWOyJSZCCFqCSE2mvfHCCH2ZZMcvYQQV8z7BYQQB4UQRbJDFk3qaOWheZJoB4y1s2w3B8raYjLQ\n3M6yJYH/2lFuMPDOI0uUAYQQJmAO8LH5kGHesgsDQEp5B5gNTMlGWTQ2cM1uATSazEJKed2B4hl6\ncTI/1O7YWfayneVuZUSmDNIOeCCltGs0lcV8D3whhPhMSnkuu4XRKLTy0GQJQog6wEbgUynl9FTK\nnAVmAB2B6sBBYKCUMsJ8PjfqzbgXUAoIB/6V+MATQmwB9kkphwkhRgN+wFngTeBvIBQYAjQE5pvr\nPAAaSym32ZDnX+bybsAswGR1bgzQQkpZRwhxDpgipfza6vwvwDkp5UAhRALQSkoZJoRoAHwFVAOu\nAYuAj6SUhhDie6CglLKTuY0mwHhz2SvATCnlFPO5XsBAYBkwFPVbXg/0lVLeFUIUQo0kmgP5gO3A\nICnladv/IQYAPyc7lksI8X/AW8At4HMp5Uyra/wAeBeoCNwGfgHetad/IUQQMAwQwD1gC9BHSnkl\nuWBSyr+FEP81y/hhKvJrshg9baVxOkIIX2At6uFjU3FYMQ74AagBnAB+FUI8Yz73NepB1h+lXI4C\n/xVClEilrVZAQaAuMBr18HkD+B2lFGJRU0o7bcjc21znA3P9ikBNqyLWUzpLgE5Wdd2BpsDiZG26\nACtRSlQA3VHTVL1t9P8SsA5Ybb7Wj4FPk63T+AEBwCvmdtoDfc3nxpv7eNlcPx74Lnk/5r4KoRTq\n+mSnaqDuTx2UgpoohOhsrtMVGIW6j8+jFHpbe/oXQtRDKe+JQCWgjbmvxCkzW6wHsmz9TJM+euSh\ncTYVUHPWc6SUX9pR/j9SyjkAQoh+wOtAZyHEEpTi6CSl3GA+3x9ogHoD/9RGW3Got9144JS5fG0p\n5SohxA3AsPWma+ZdYJaUMtTc19tAs1TKLgb+LYQoJaX8A/UQ/0NKuStZOXegCHBZShkDxAghmgG2\nprUGA+ullF+YP58WQpRDPWC/MR9zRb2tXwWOCyHWA7XN5yqarz9aSnlTCNHHfMwWNVGjqmPJjv8F\n9JZS/g2cEELURingpcAf5nPrzGVjhBC/AVXt6P8e8I6U8ieruqus6triGOAjhChgnjLUZDN65KFx\nNlOB0kB04gEhRAMhxC3zdlMIMcuq/I7EHSnlfdTUVVXAG/V93W113kCNGqqk0ne0WXEkchPIbafc\nVYEDVn3dRY10UiClPIx6uHU0H+oE/MdGuWvABGCaEOIPIcR84FmzIklOZSC58tkBlBZCFDZ/vmVW\nHIlYX98EcxuXhRC/AoHAkVSutQQQZ77f1hw2K45EIjA/4KWUW4FoIcQ4IcQyIcRR1MggV3r9Syn3\nA7uEEJ8IIX4SQhxAjSZzkTqx5r/F0yijyUK08tA4m8XAR8CXQojEH/4+4AXzVh01/ZFIfNLq5AIe\noNYsTKTEhdS/x8kfhqTShi0MG2VttZfIT0Ans0lpE/PnFEgpR6GmaiYC5YC1QohPbBT928Yxl2R/\nU70+KeU+1Jt+D+ACMAb1wM5ro46tawV135P3fx8say6/Ax6o6bWuwJrEgmn1b17LOQx4Ab8B/VBr\nSmmRqFiSy6TJJrTy0Dib5cD/AZfMf5FS3pNSnrHarN+e/RN3hBD5UYvFkcBp4B/UHL81dYHjjyBX\nemaoh4EXrWTJQ9rTKj+Zy78DnJBSpnjLF0IUF0J8DfxPSvl/UspXUUqkq432jpPyWuujprzStSoT\nQrwPNJRShkgp30LdpyqodZLk/AEUtKFYqpjXaRIJ4OHU1gBgkpRyoJTyeynlIZRSNNnRfz9gqZTy\nTSnlHLPBg6VuKhQ1//1feteuyRr0mofG2ZiklPFCiAHAFiHEAinlr2mU7yuECEdNGX2KetMNMVvw\nBANThRBxKCuq91Fvt3MfQa7bQCHzYv4ZKeW9ZOenAj+Yp1R2oyyDipIKUspoIcRus8wTUin2F2pR\nuYAQ4gugEGodZY+NspOBfUKIkSiLqlrACMCedSOAMsAgIcSbqDf/t4AbgLRR9hBqxPcCYG2qWwyY\nb3a8rGtu4w3zuVigsRBiKeol9APA19xWev3HAo2EEDVQ/4fewGtYTUna4AXUNNo/dl6/xsnokYfG\nmVje7qWUv6HezmelMnWSyHcoy55w1Px2EyllnPncCNRi7feo+Xdf4GUp5Znk/aUnD7AJNaLZjw0r\nHvNC+fsoR8II1MN1czrtL0ZZdy2x1a95/aUFSuGFo6yujpj7Sd7/QdTCeyDqgTwOGG2n0QEoJbbO\nLMtRoDHKtPimjb5uodZTGiU79StqcXsvyvJsoJQy0dnxffN17TOXyw18wUOLtLT6Hw1EoaasdqDW\nRv4FVDaP8GzREGUKrMkhmHQmQU1OweznMVlKmd78tyaTEUIEonxwXshuWZJjNhC4AFSTUp7Pbnk0\nCj3y0Gg0oNamEELUz25BbPA2sFwrjpyFVh6anIQeBmcTZrPnPsBn2S2LNUKIAijloT3Lcxh62kqj\n0Wg0DqNHHhqNRqNxGG2qmw5CiI9Q5om5UeEqvs9mkZyOEOJF4EspZWMhhBewAEgAjkgpB2SrcBkg\n2XUVQ5n4PoNyQOsppTybxfK4omI8VQTyoEx8j+Hk+51Kv7tx8H6YfUDmomJYxaMCULpktvyp9HPb\nHnnNjqnhqFhjDx5FtmRtXLOz3/1Aoj/OWeBbYDrKV+m/UspxdtT5HBXaJzfK6q0LykfG8jwCtj3K\nNWUGeuSRBkKIRkA9KWUAKsBbueyVyPkIIT5E/TgSzWmnAh9LKRsBLkKINtkmXAawcV2TgEVSypdR\nZqU+2SBWd+CqlLIhKobX12TN/bbut4W530e5H61R8cEaoMxvpzlJflv9pCuvWUnO5mHofIdls9GG\nPf3mNcv7inl729xGFynlS8CLQojqdtT5Fhhp7ms20Jmkz6Pyj3JNmYVWHmnTHDgihFiJim76NNiZ\nn0bldkjEX0q53by/DvX29TiS/LrqA2XNob67AVuzQaZlPAzo6IJ6q66ZBffbul+Tud8AoJwj90NK\nuYqHUXQroLy/M13+VPqxR94pqCCSl1DX+SiyWbcB9n1vXkB57G8QQmw0R0jOIx/mItmACmGTVp26\nKD+nN4RKNVAXpSySP4+y4vtiE6080qYoKlxGICpwm814RU8SUsoVJI0vZR0y4hYqMuxjh43rqgj8\nZQ4REoOKv5XVMt2RUsYJIdyAEGAkWXC/U+n3OSDW0fshpUwQQixA5WFZjpPkt9FPmvIKFVL/stmp\nMVEm6+ddurLZaMOEfd+bOyh/peao58b3JE0cZqvv5HV+QjlP/iqlbAw8i1KY1s+jxY5eU2ailUfa\nxAIbpJTxUsqTwN9CiFRDVDyhJFjtu/FwTvZxJ5aHgfzWYBVTKysxh1nfDPwgpVxCFt1vG/1e5RHv\nh5SyNyrq8Twgv9WpTJU/WT/XSFveN4FXzW/tLwA/osKtOCKbrTbi0+kX4CTmXC5SylOosCwe6fSd\nvM5V835ikrJfUOFskjyPSKossvT3qZVH2uxAxdxBCFEaKMDD0NBPC/uFEA3N+6+jMsI9CWznYViS\nhqQSbt2ZCJXEagMwTEr5g/nwAWff71T63YGD90MI0d1sUALqQfYACDevFUImyZ9KP78BLVOTV0rZ\nSErZ2PzWHomK7rvOkXubShu/pNWvmbdQ2SKtnxtxQojnhMoV39xG38nruAERQmWeTOxrD0mfRwWB\nTZl9v+1FW1ulgZRyrRDiJSHEXtSQ9T2zM9XTxL+BuUKlgD2OSuX6JPBvYJ5QCaJuoOavs5oRKKud\nT4UQo1BOku8DwU6+37b67QV85+D9+Bn4XqgkUK6oBFYnUPc1M+W31c9Bs7zvOiBvZnyXE783afX7\nnVne7aiR5Jvmvz+hXth/lSpkfXp14lCx4HKhrK96AeOtnkf9gXNk/v22C+0kqNFoNBqHeexHHkLl\nQ+6H+a3NVtRQjUaj0WQuT8KaR1/z9h3KiUaj0Wg0TiZHjzySeQSbUB6VL6AWzd4x53FwkVLeF0L8\nD3glG8XVaDSap4YcO/Kw4RHcFshr9q4cgfKsBLhjTiBTCp2iUqPRaLKEHKs8SOkR3ABYDyCl3MND\n++q5wBzU1NWirBRQo9FonlZy7LSVlHKFEKKC1aHCKNO4RB4IIVyklPtRZm1pEhERoc3KNBqNxkH8\n/f1NNk8YhpFjN29v7wre3t47zftfeXt7B1qdi3akrfDwcCMjzJgxI0P1H2ecde0XL150Srv2khP/\np9klU2K/Gf2fZJX8j9pPRuXLyn6t66T2f3H2/TY/N20+U3PytFVyfsfsAWsOGnY4e8XRaDSap5cc\nO21lgxWoODO/mz+nO1Wl0Wg0GueQo5WHOeF9gHnfQLnjazQajSabydHKQ6PRZB2HDx/m/v37/PHH\nH4/cxqVLl4iIiMhEqTK3n4zKl5X9Jq9j6/+Smfc7T548VKtWze7yWnnYSZ06dbJbhGzjSb32nHhd\n2SVTnTp1uH//Pv7+GYtMHx8fn+E2nNlPRuXLyn7tqZOZ99tRJfTUBEaMiIgwsuJLrbGfS5cuUbp0\n6ewWQ2MmIiIiSx78mpxJREQE27dvp1ixYgQFBVmOpWaq+zhZW2k0Go3GiRQvXpwrV67YVVYrD41G\no9E4jFYeGo1Go3EYrTw0Go1G4zBaeWg0Gs1jRHh4OKNHj2bjxo2WY/v27eP+/fucPHmSnTt3MnXq\n1DRayBy0qa5Go9E8RhiGwdixY5Mcu3PnDnny5CE6Oppq1aqxdu1ap8uhRx4ajeaxICEhwbJ///79\nbJQkeylevDinTp0iJCQEgLi4OAoVKgRA06ZNuXbtGlWrVnW6HE+M8hBCNBZCzM1uOTQaTca5evUq\nM2bMYNGiRaxatYo1a9awdOlSALZs2UJcXBwAM2fOZOPGjcyePRsAKSWHDh3KNrkdZcGCBVSpUsVi\nHrt//37at2/PqlWr+PPPP9m+fTs7duxg+/btREZGAhAZGUnFihW5cOECd+/eZf/+/dSsWROAuXPn\nUr58eaKjozl37pxTZX8ipq2EEF5ATR5mHdRoNI8pFy5cYOzYsXz11VcULlwYgIkTJ9KkSROuXLlC\nXFwcRYoUYdeuXYB62z5x4gTh4eHUqlWL+fPnU7lyZVxdc/7jzdfXlyZNmrB27Vp69+5NzZo16dOn\nD6+//joAJUqUSFGnSpUqHD9+nDJlypA/f34SEhIwmZQfX/Xq1Tl69Ch58+Z1ugNujh15CCFeFEJs\nMe+bhBDfCCF2CiE2CyE8rctKKaOklF9lj6QajSYzGT58OP369bMoDgB/f3/8/PxYvnw5TZs2BZT3\nc+XKlQH1EN69ezcAAQEBrF+/PusFfwQuXrzIoEGDWL16NZB0Cio1nn/+efz8/OjUqRN//vknJUuW\ntJyrXbs2tWvXZsiQIeTJk8epsudI1WzOX94DuG0+ZMlfLoR4EZW/vK0QYhzwPPCelPI6YDvjVQaJ\njoZSpSB3bme0rtFoEomMjOTmzZvUqlUryfFEhREbG0u+fPkA+Ouvv8ifPz8ABQoU4OrVqwD4+PgQ\nGhpKq1atkrQRHx/P6NGjiY6OJiEhgffff5+LFy+yfPlyDMNg0KBB/PHHH5bP/fv3Z/Xq1URHR2MY\nBr169aJFixasWLEiSZ26des+8vW6uLhQqVIlTCYTZ86c4erVq1SpUsXu+qdPn6Z+/fqP3H9GyKkj\nj/Tyl9cy74+SUnYzKw4ApwTqGj8evL1hzhy4d88ZPWg0OY+qVatiMpkybbNnEffAgQNJIrtGREQw\nfvx4xo0bx40bN7hn9QNMSEggV65cKfZTIyQkBA8PDxYuXMjMmTMZN24cJpMJd3d3Fi9ebFECiZ/P\nnz+Ph4cHS5YsYf78+UyfPp3r168nKZOa4pg3bx5ffvkl0dHRlmMxMTGWtZpEEqeb2rVrx+rVq4mN\njcXDwyPd+5RIdikOyKEjDzvyl8eb85cnJKvXM612L1269EjyjBoFLVrkJjjYjTFjcvPee7fp1u0O\n+fM/HUElncWtW7ce+X+icT5HjhzJ8j4TH+aJ+Pv7s3HjRvz9/XF3dyc+Pt5yrmjRoty5cweA27dv\nU6RIEcu5xOPWnDx5koiICA4ePIhhGDx48IDr16/z3HPPJSmX+DkqKoqAgAAAChYsiJeXl0UZJK9j\nza5du2jevDnlypVj8eLF1KxZ0zKt1rFjR0u5K1euUKpUKQBatmxJ165d8fb2tu9GOYlbt25x584d\nu36XOVJ52OAm4Gb1OYXisIdHXUBasADmzoWPPoKSJeHzz92ZNcudoUOhf39IZ4pSkwo6qm7OIiN5\nPDKLhg0bMnLkSAzDwGQyYRgGu3btonfv3gBJFsH9/f05fPgwjRo14tChQ9SrV89yztYoxNPTk1Kl\nStG3b1/u3bvH7NmzcXNz46+//kpSzsVFTch4eXkRHh5O06ZNuX37NqdOnaJs2bJERUVZytjimWee\noVy5cgAEBQURFhbG6tWradKkSZJyR44c4cUXXwSgSJEieHl5WWRZtWoVlSpVsqzpZBVubm7cvXvX\n8rtM6zuRU6etkpOt+ct79ID334fRo+HNN6F9e1i7FiIiwNNTTWtdv55+OxqNJm08PT3p2bMnEyZM\nYOXKlaxdu5a+fftarI7y5n1oUFm3bl2uXbvG+vXrMZlMNGjQwHIucV3Ems6dOxMVFUWPHj3o2rUr\npUuXTnOqq1OnTly/fp1u3brRq1cvBg4caNeUkq+vb5LPLVq0YPjw4UnWcXbv3k1wcDDbtm2zHOvQ\noYNFWZQqVSrFFJc1UVFR7Nu3L11ZnEmOzedhnrb6j3mR3ATMAvzMp9+UUp50pL3MyOdhGPDf/8KX\nX0JUFOzYAXfuwOefK2XSvz8MGQLPPpuhbp4a9MgjZ/E45POYP38+gYGBSSyxkhMdHc3evXsJDAzM\nQskyl7179/LgwQP27NlDYGAgv/32G66urhQoUIAKFSogpeTu3bu0b98+Xesse4mIiEBKyeXLlxky\nZIjlWGr5PHLstFVOzF9uMkGzZmo7cADKllXHfvgBzpxRSsXbG95+G4YOVVNcGo0m8+jYsSNhYWF0\n7tw51TJbt25N8/zjwqZNm2jYsCFly5a1WGS98MILTJ06NcUUWHbwuExb5Thq1FCKIxFPT/j2W4iM\nhLt3oXJlNdV14UL2yajRPGm4ubnh5eWV6lx8TEwMPj4+Saa3HkfOnj2Lr68vYWFhnDhxAlBTVVu3\nbuWNN96gXLlynD592mKenB3k2JFHjiI2FooUgTQWyRJZvFhNZa1cCatXg58fdOoEw4dDGgYaGo3G\nTpL7gFhTokQJy2L140ziyKlDhw4ALF26lFy5ctGsWTNLma5du2aLbInokYc9vP8+VKkC8+bB33+n\nWbRPH6hQAQID4exZWLpUrYHUqqUW2086tFKj0Wgcwdle1dlF586dk5j55gS08rCHhQth1ixYsUIN\nH8aPV6MRGzz7rPILOXsWGjVS6x+7d8Phw6pq/frQrRtkgwm9RqPRZBpaediDyQSNGyuTqo0b1ep4\npUowaJDat0HBgjB4sLLKGjYMSpdWSuXMGXjhBWjaFDp0gP37s/haNBqNJhPQysNRqlSB+fPh6FFw\nc4M6daBjR9izx2bx3LmhefOHn93c1PrHmTPw0kvQujW0aqVGJxqNRvO4oJXHo1KqlHLwOHcOGjSA\nzp2hYUO1Sp6QvvN7gQLw118wYAC88oqq/uqr8Ntvzhddo9FoMopWHnYwefLkJPmCk1CokFpQP31a\naYJx45Sd7ty56S6ud+ig1j4+/xy6doWWLdUaScOG8OuvyilRo9FociJaedhBrbJlGf/WW7Rt25ao\nqCjbhVxd1fBh3z6YPRtWrYKKFeGzzyAVW+wXXoCffoK9e+HmTaV3WrSAvn2Vp3rdurBmjVYiGo0m\n56GVhx00dndny99/M/3wYb6qXp1Phw7l1q1btgubTPDyy/DLL7B5M5w/rxbXBwxQq+c28PRUxlzH\nj6sZsO7d1Yjkww/h00+VQ2JoqF2zYRqNRpMlaOVhDy1aYIqJocKkSUytU4fhM2eyumRJlk+eTEJa\nT/TKlZVvyPHj8MwzaigRGJjq6niJEsqhEJQ/YmCgCoPy2WcwaRJUraqcEK2iUms0Gk22YJfyEELU\nFEI0EUIEOFsgRxFCvCKE+EEIsUwIUS39Go9I3rzQoQP5Nm2i0NmzvNS7Nz8tWkS9evUs6S9TpWRJ\nmDDhofNH165qiLFqVbrDCZNJWWQ1bw5dusA334CvrzL4+uefTLw+jUajcYB0lYcQ4mXgTynlJiBS\nCDHY6VI5Rn4pZS/gc6BZeoUzhdKlKT9zJiEHDjBw4EA6dOhAz549VQIVw0j9qV6okPINOXVKOYGM\nH680wZw5KiBWGvj5qZAnN26opZXFi9Vs2DffpLsur9FonhLCw8MZPXp0CgOfffv2cfLkSXbu3MnU\nqVMzpS97Rh7uUsqLAFLKO4DTY8UKIV4UQmwx75uEEN8IIXYKITYLITyty0op1wohCgCDgB+cLZs1\nLi4u9OjRgxMnTlC2bFmqVavG7KFDMcqWVWF1D6eSdsTVVc1P7d2rrLLWrlXu5+PGpbq43rGjyh8y\nZQr8/rtaPunWTVX18oL/+z8VU0uj0Ty9GIbB2LFjLTnfE4mLiyM6OhovLy9iU4mO4Sj2KA8XIYQX\ngBDiVcCpUceEEB8Cc4HEsJhtgbxSygBgBDDVXG6cEOInIYQHMB0YJaXMlhCTbm5ufP755+zbt48N\n587xSp48nIiJwXj9dahdWw0PbGWLMpke+oZs2QIxMWo48d57anRio3jz5qrokiVq0PLLL8oia9s2\ntfA+cSKktpav0TzOWK8v3r9/PxslybkUL16cU6dOERISYjkWFxeHm5sbTZs25dq1a3blkreHdKPq\nmvOJvyuE+BiojPOnhk4D7YCF5s8NgPVmWfYIIWqZ90cBCCF+AIoCXwghVkopf06tYWfny86XLx8z\nZ85k27ZttB0zhuIVKjC9dWt81q3j/uXLxPXpk3pld3cYOxaXAQMouGABBerW5f6LL3L73Xf5x0YU\n0fLl1XbpklpS+fprOHHCleDgQkyalJe33orjrbficHfPuXa+Ooe5JjWuXr3KTz/9hIeHB25ubri4\nuHD79m26du3Kli1bqF69OqtXryY4OJgiRYpQq1YtpJQEBgbSpk0bm21KKbl37x5+fn42z+c0FixY\nwOTJk9m6dSvFihVj//79jB8/nl69etGmTRv+/PNPTp48icmcG6JgwYJER0fTokULfvnlF+7evUv+\n/PnZv38/9evXZ+7cuQQFBbFq1SrOnTtHxYoVU/TpSA5zDMNIc/P29s7r7e293tvbe7q3t/dn3t7e\nxdOrk9HN29u7gre3907z/lxvb+/mVufOeXt7uzjaZnh4uJGV/PPPP8bXX39tFCtWzBgwYIBx9epV\nxxq4fdswgoMN47nnDCMgwDB+/tkw4uPTrZaQYBgffGAY7doZhoeHYXz8sWFcufKIF+FkLl68mN0i\naKzI6t9IasTExBjvvPOOcePGDcuxL7/80ti3b59x+fJlY82aNZbj3bt3Nz777DPDMAzjl19+MXx9\nfY2TJ0+m2vb8+fONf/75x3nCZyK7d+82Bg0aZHz//feWY2FhYWnWOXXqlHHw4EFj6dKllmNbt241\nDMMw9u7da+zdu9eYNm2ace/evRR1w8PDjcWLFxvTpk1LcsxI5Zlqz7TVEGCBlPJ9YDTQwwHlmRnc\nBNysPrtIKXO8x4OrqysDBgzg+PHjgMpr/PXXXxNvbWebkKAWLhYuTLlgUbAgDByopq8++EClKfT1\nVQ6IaSyuP3igBjE7dqiwW8ePq+yG//43/O9/zrhSjSZzGT58OP369UuSatbf3x8/Pz+WL1+eYj4/\nEQ8PDxISEjh79myqbQcEBLB+/fpMl9kZXLx4kUGDBrF69WpATT+ll3L2+eefx8/Pj05mm/8///yT\nkuaUprVr16Z27doMGTIkU0LX26M8cksplwCYH9pZPc/wO9ACQAhRF0hlFTpn8uyzz/L111+zadMm\nVqxYQfXq1dm0aZM6aRjQvr1awChbFvr1Uz4g1i7luXI99A357jtYt055ro8ZA1eupOjP1RVGj1aB\nF5s3Vw7vVaqoEFyVKytjr5iYrLhyjcZxIiMjuXnzZoqET02bNiVPnjzExsaSL18+m3XPnDmDm5sb\nNWrUSLV9IQSRkZEpjsfHxzNy5Eh69OhBUFAQ+/btA2DFihV0796doKAgQkNDLfs7duxg2LBhdOnS\nhc6dOxMWFpaifLom/Ong4uJiST975swZjh49SpUqVRxq4/Tp0wghMiRHqvLZUWaPEOIlACFEMeCG\nUyRJnRXAPSHE78BXwAdZ3H+mUK1aNTZu3Mhnn31Gnz59aNeuHVHnzinFsHatssyqWBF69oSgoJQN\nmEwqDO+qVSp64h9/gBDQv7/NDFOFCqkQJ1FRKl5W27Zw7Bjky6fCovTtm2o0eY0GUE6pJlPmbfas\n0x44cIBq1R66a0VERDB+/HjGjRvHjRs3uHfvXoo6J06cYNasWUybNo3g4GCKFSvm8LWGhITg4eHB\nwoULmTlzJmPHjrWcc3d3Z/HixeTKlcuyf/78eTw8PFiyZAnz589n+vTpXDcbxSSWqVu3rs2+5s2b\nx5dffkmVnXbWAAAgAElEQVR0dLTlWExMDHFxcUnKJa5ltGvXjtWrVxMbG4uHh4dD11W/fn2HyjtC\nuspDSvlfIFoI0RTwkFKGOU2ah32eN1tXIaU0pJT9pZT1zdtjm4vPZDLRrl07jh07Rp06dXjxxRf5\n+OOPVaiTMmVgxAiQEqZPT7shHx/lG3L8OBQrphwO27WDnTtTFM2TB3r3ViFPSpaEyZOVrilRQk1r\n9eqlutRoknPkiBoEZ9ZmTwI0k8mEu7u75bO/vz+5c+cmICAAd3f3pNO+Znx8fHjvvfeoWrUqc+bM\nSXJu06ZNXL58OcmxOzZs2k+ePMlvv/1Gz549GTRoEA8ePODGDfWe/JxV/ujE/aioKMvoqGDBgnh5\neVmUwXNp5JvetWsXzZs356OPPmL79u2Wae3du3dTsGBBS7krV65QqlQpAFq2bMn69esxcliQO7s8\nzM0P841S6sdMZpAvXz5GjBjBwYMHuXDhAj4+Pvz444/KFNFkUgrBFqtWJf0FliihfEPOnlXZpXr0\ngIAA+PlntfiRCkWLqpAnUsK1ayq7YZcuqbulaDRZRcOGDYmMjLQ8KA3DYNeuXZbRiKtr6gaipUqV\n4sCBAxiGQWxsLFevXmXFihUpHrq5cuVKUdfT05NWrVrx448/Mm/ePF577TWLEnNxefiYTNz38vIi\nPDwcgNu3b3Pq1CnKli2bonxynnnmGUuO9aCgIM6ePcvEiRNTKJwjR45YTGqLFCmCl5cXf/31FwCr\nVq3i2LFjqfaRVdgbniSXEOJ7IcQkIUQjIUQu83Ef54r3ZFOmTBl+/PFHli9fzsyZMwkICGBPKkml\nADVkeO01NWSYPfuh70jBgirw4smT8K9/qUBYPj7KvyQNz0EXF5VXxDDU8kmTJmoAExGRyReq0diJ\np6cnPXv2ZMKECaxcuZK1a9fSt29fSpQoAUDevHktZUNDQ4mJiSEiIoKNGzfStm1bSpYsyaRJk9iz\nZw9FixbFxyflI8rWmknnzp2JioqiR48edO3aldKlS6cpZ6dOnbh+/TrdunWjV69eDBw40K4pJV9f\n3ySfW7RowfDhw5Os8ezevZvg4GC2bdtmOdahQwcqV64MKCWZfIrLmqioKMuajVNJzQwr+ebt7e3l\n7e29y9vbe7S3t3eA+dg8b2/v8va2kZ1bTjFDTI0HDx4YCxYsMEqXLm307NkzdTPW+HjDWLfOMDp2\nNAx3d8Po3t0wHjxIWiYhwTC2bzeMNm0Mo1gxwxg1yjD+/DPVvk+dMox+/QyjSBHDeOklwyhRwjBe\nf90wdu7MxAu0gTbVzVnk9N+IYRjGd999l8SENz2Cg4ON//3vf5bP58+fN0JCQpwhWpaxZ88eY+fO\nnca0adOMmJgYY9GiRcaSJUuM1atXGwcPHjSWLVtm/PDDD8atW7ccatcZprqJdAAaSCnHSikTJ9cr\nAGeFEJFCiC6Zr9qeHlxcXOjVqxcnTpygdOnS+Pn58cUXX/B38sBVuXKp0ceyZWo1vH17NYSwxmRS\n6yArV8L27fDnn2pxvV8/mwsczz+vBjJHj6pZr2nToE0bFb+xWbNUI6ZoNFlOx44dWbdunV1lY2Nj\nOXv2bBKrp61bt9K6dWtniZdlbNq0iZo1a1K2bFkqVaqEp6cnzZs3JywsjIoVK6YY4TgDh0KySymT\nT6TvAXyBbSjPcE0GcXNz44svvmDPnj3s2bOHKlWqsHLlStuLZc8+q+aZbPG//6kpKyGUZpBSrZi/\n9JIyvdqxI0WWqVKllDtJ165Kz5w6pXKJvPIKJFtz1GiyBTc3N7y8vPjjjz/SLfvss8/y1VdfWTzO\nY2Ji8PHxSTL19Thy9uxZfH19CQsL48SJE4Caqtq6dStvvPEG5cqV4/Tp01x19ltfakOS5Ju3t/dY\nG8daWO0Psret7NgehyG5LX799VejcuXKRpMmTYzDhw/bX3HKFDUP1bevYezeraayDMMw4uIMY9Ys\nw3j+ecOoW9cwQkPT9FxPSDCMTz81jMqVDcNq9J8p6GmrnMXj+huxF1te1U8CS5YsMZYtW5bhdpw5\nbVU0+YFkZrumjKsyTXJeffVVDh48SJs2bXjllVcYNGiQxeoiTf71Lzh0CCpUUHa6VaqokLzx8co3\n5MQJlapwyhTlgj5zJthYhDOZYOxYNX318svKvUSjeRzJDK/qnEjnzp3p2LFjlvfriPL4TQiRloPe\nMxkVRmMbV1dXBg0axLFjx0hISMDHx4eZM2fatHlPQtmy8PHHygprzhzlJZhYJ1cutV6ya5cKj7Jp\nkwoLP2qUWiOx4vBh5QTfpg00bqyCMWo0mqcbR5RHCNBOCNEv+QkhhAlwjg+8xkLRokWZOXMmmzZt\nYvny5dSoUYPNmzenXzHRO33+fLBlTlivnvIN2bFD2ez6+CgXdPN8qp8fvPmmUiK9eqkRyIULmXtt\nGo3m8cJu5SGlNIAgYJgQYocQorsQorIQojoqCdN+ZwmpSUq1atXYtGkTY8eO5e2336Z9+/acyUis\nkbVr4cUXYfNm+OILNVIpU0alzG3bFo4fZ8wYtWheuDD06aMUiI6RpdE8vThqbRUD1AWige9RQQrD\ngXvAtEyXzg7M+dUXm50YHQ9q85hiMplo3749x48fp1atWtSpU4eRI0dy+/Ztxxt77TUVTXHjRhVf\n64MPlKnvmTNKgTRsSJ7hH7D4m5uMGQMtW6p8VY0awfnzmX1lGo3mccAh5QEgpbwipewGlADqAeWl\nlH2yMUx6XqA/EGaW56kiX758fPzxxxw8eJDo6Gh8fHxYuHBhkqxr6eLqCi1aQGgonD6tsh8OHQph\nYUqRHDsGcXF4t/Lmi9e38eUXCQwdCu+/r0YgaUTA1mg0TygOK49EpJR/SSn3SikzffnUwRzmu4Aq\nwL+AA5kty+NCmTJlWLhwISEhIQQHB1O/fn327t3reENFiyqtEBkJHTqoY8WKwbffQlgYb0d9zHdH\n68Hvv/P++ypPSOPGOkKvRvO08cjKw1k8Qg7z2qipsxbA4GwQOUdRr149du/eTb9+/Wjbti29e/e2\ny6EqBSZTSs/1mjUxbd9G7g+HqEiKQUEMaHOBjz5SI5DT2k1Uo3lqyHHKg4c5zBNJksMcsOQwN0+f\nuQHzgekoi7CnHhcXF3r37s2JEycoWbIk1apVY+LEiTZzITjM7NmwdCn88gt4ekL16rwbO4ExH/1N\n48Y2U4toNJonkNTjG2cTUsoVQogKVocKkzQBVbwQwpKKVkq5GbDDXhX7kro/YQwePJjWrVszbtw4\nZs+ezahRo2jWrJkl0YzDvPYahS5domDjxtzt2JE7//kPbtOm0fPotxRtOp5GDbuydNlfPP98Oj4o\nwK1bt57K/4lGk1O5desWd+7cset3mSHlIYT4Qko5IiNt2EGm5TBPL8zyk0rp0qXZsGEDv/76K0OG\nDOE///kP06ZNczilpYXx42HQIAp9+ikFgrpzqutoxAcf8MYHH1Cr6I90C5zOrK2VMUeQTpVLly49\ntf+TnMgjTW9qchTh4eGsWbOGl156KUmu93379uHu7s7Vq1fZvXs3Q4cOtVnfzc2Nu3fvWn6XaX0n\nMjptlRVe5Y91DvOcRLNmzTh48CCtWrXi5ZdfZvDgwfaFOrFFiRLw7bec/WYDv849x6kKTSEyktL9\n3mD93y+zs/b7HPv9WuZegEajSRPDMBg7dmwSxQEQFxdHdHQ0Xl5exMbGZkpfGVUeWZEX8YnIYZ5T\nyJ07N4MHD+b48ePEx8fj6+vLN998k36ok1Tw6lAdJk4iqLuJfwxXGDSIfGeO0ajefYo19OHCp3PS\nzGqo0diLtfn5/fv3s1GSnEvx4sU5deoUISEPl3/j4uJwc3OjadOmXLt2zZKhMKPkxAXzJzaHeU6i\naNGizJo1i//+978sW7aMmjVrsmXLlkdqa+BAFR1+7FhL41Ta+A2R49cSM3Exd6rUUnlFNBo7uXr1\nKjNmzGDRokWsWrWKNWvWsHTpUgC2bNlCXFwcoaGhNGrUiLZt2zJ+/Hh69OjBqlWrUm1TSsmhQ4ey\n6hIyzIIFC6hSpQpXrlwBYP/+/bRv395yjX/++Sfbt29nx44d7NixgwMHDhAZGUnFihW5cOECd+/e\ntdSrUaMGc+fOpXz58kRHR3Pu3LkMy5cjlYcm6/Dz82Pz5s2MHj2at956iw4dOnDWQa8/kwm+/x6+\n+85KR9y6xatzO5G/R0f+9ce/ud8xSCUK0TFNNOlw4cIFRowYQe/evenevTtt2rTh2LFjVKpUiStX\nrhAXF0eRIkUIDAykfPny1KpVi08++YQuXbowYsQITp06ZbNdIQQRERGPPMrOanx9fWnSpAlr164F\noGbNmvTp08eSn6REiRK89NJLNGjQgAYNGlCjRg2qVKnC8ePHKVOmDPnz5wfUiM3FxYXq1atz9OhR\n8ubNmylrjVp5aDCZTHTo0IFjx45Ro0YNy4/RkVAnJUvC3Lnw22/mA25usGoV1aNXM7HwBPr8PYNL\nhbxVdqnPPgPzW5FGk5zhw4fTr18/ChcubDnm7++Pn58fy5cvTzGfn4iHhwcJCQlpvvwEBASwfv36\nTJfZGVy8eJFBgwaxevVqQE0/FSpUKM06zz//PH5+fnTq1AlQo5OSJUsCULt2bWrXrs2QIUMyJTy9\nVh4aC/nz5+eTTz7h4MGDnDt3Dh8fHxYvXmw7i6ENWrWCTz6xOlCtGvz6K4W/mcTMgsM5/uNejn0W\nqvKMVK5MvrVrU2Qz1DzdREZGcvPmTWrVqpXkeNOmTcmTJw+xsbHky5fPZt0zZ87g5uZGjRo1Um1f\nCEFkZGSK4/Hx8YwcOZIePXoQFBTEvn37AFixYgXdu3cnKCiI0NBQy/6OHTsYNmwYXbp0oXPnzoSF\nhaUob53+9lFwcXGhUqVKmEwmzpw5w9GjRx22kDx9+jRCOCfgeUaVh04A9QRStmxZFi1axLJly/jq\nq69o1qwZUVFRj9aYyQStWlHo7GFK9mrOOyNLsmtoCMyfj9vUqdCkiYr1rsl5VK2q/n+ZtdmxUHvg\nwAGqVatm+RwREcH48eMZN24cN27csOnoeuLECWbNmsW0adMIDg6mWDHH46OGhITg4eHBwoULmTlz\nJmMtC3jg7u7O4sWLyZUrl2X//PnzeHh4sGTJEubPn8/06dO5fv16kvJ169a12de8efP48ssviY6O\nthyLiYkhLlkytkRfrHbt2rF69WpiY2PxsJVSIQ3q16/vUHlHyKjyOJYpUmhyJAEBAezdu5fmzZvz\n4osvMmXKlEefL86ThyrfDuHTxT60aQO/52nMlQ0bVPysJk3Uqvujmg1rnMORI2pkmFnbkSPpdmky\nmXB3d7d89vf3J3fu3AQEBODu7m7z++fj48N7771H1apVmTNnjuX44cOH2blzJ3Pnzk1S/s6dOyna\nOHnyJL/99hs9e/Zk0KBBPHjwgBs3lG/yc889ZymXuB8VFWUZHRUsWBAvLy+LMrAun5xdu3bRvHlz\nPvroI7Zv387x48cB2L17NwULFrSUu3LlCqVKlQKgZcuWrF+/3u4ZgKwiQ8pDShmcWYJociaurq78\n+9//Zs+ePWzYsIEXX3yR/fsfPXXL66/DokXQrh3siSgAAwbA/v2QkAC+vvDNN9q09ymmYcOGREZG\nWh6UhmGwa9cuy2jE1TV1v+ZSpUpx4MABDMMgNjaWo0ePUqtWLa5du5ZEYeTKlStFXU9PT1q1asWP\nP/7IvHnzeO211yxKzMUqxlvivpeXF+Hh4QDcvn2bU6dOUbZs2RTlk/PMM89Qrlw5AIKCgjh79iwT\nJ05MoXCOHDliMaktUqQIXl5eFp+sVatWcexY9r+36zUPjV14eXnx66+/MnjwYF5//XWGDRtm8w3O\nmoMHVSis5DRrBj/9BH36FGHrVlQe9dOnYdYsWLYMata0WnnXPE14enrSs2dPJkyYwMqVK1m7di19\n+/alRIkSAOTNm9dSNjQ0lJiYGCIiIti4cSNt27alZMmSTJo0iT179tClSxdcXV0xDIMCBQpY6tla\nM+ncuTNRUVH06NGDrl27pmuN1KlTJ65fv063bt3o1asXAwcOtGtKydfXN8nnFi1aMHz48CRrPLt3\n7yY4OJht27ZZjnXo0IHK5pANpUqVSjHFZU1UVJRlzcaZmHLaUMhZREREGP7+/tktxhPB5cuXGTJk\nCHv27OHbb7+lSZMmNstduAD+/rB6tUpUmJyQkKsMGFCUpYv+obGcraywOnSAWrVg3DhVafJkqFAh\nZWVNphMREUFO/43Mnz+fwMDAJJZYaREWFka9evVwc3PD1dWV6Oho9u7dS2BgoJMldR579+7lwYMH\n7Nmzh8DAQH777TdcXV0pUKAAFSpUQErJ3bt3ad++fbrWWdZEREQgpbT8vhOP+fv721zb1iMPjcMU\nL16cn376iRkzZvDWW2/x5ptv2gx5ULYszJwJ3buDLavf+vXvExICnbvn5r8+g1TO9Dx5YPhwGDwY\nqlRR2mfsWEhnlKN5OujYsSPr1q2zq+yaNWv4/fffmTp1qmUqaevWrbRu3dqZImYJmzZtombNmpQt\nW5ZKlSrh6elJ8+bNCQsLo2LFiilGOM7ALuUhhMhlTvM6SQjRSAiRy3zcx7niaXIyLVu25MiRIxQu\nXJiqVauyZMmSFIt6gYHw0ksqv5QtGjWC5cshKAg27POA6dNhxw4oXlylxt2/X2UyrFwZQkK0ae9T\njpubG15eXnYFcWzdujUTJkzgs88+w8XFhZiYGHx8fJJMfT2OnD17Fl9fX8LCwjhx4gSgpqq2bt3K\nG2+8Qbly5Th9+jRXr151riCGYdi1eXt7e3l7e+/y9vYe7e3tHWA+Ns/b27u8vW04a/P29i7u7e29\nL60y4eHhhsZ57Nq1y6hatarRsmVL4/z580nO3bxpGF5ehhEamrTOxYsXLfu//24YxYoZxtq1qXSw\ndath+PkZRqNGhnHwYOYKrzEMwzCe9N/IvXv3slsEp7BkyRJj2bJlGW4nPDzcWLx4sTFt2rQkx4xU\nnqmOTFt1ABpIKcdKKXeaj1UAzgohIoUQXTJftdnNh8C5bOz/qadu3bpERERQt25datasSXBwMA/M\nVlNubmqBPHfu1OsHBKi1kd69VZ6pFDRqBBER0Lo1vPoqvPceZFJ0UM3TQWZ4VedEOnfuTMeOHbO8\nX4fWPKSUyW0o9wC+wDZUBsBMwZEc5kKId4HFwN+Z1b/m0ciTJw+ffPIJO3bsICQkhAYNGnDEbNtf\npw688Uba9evWVYrj7bfBZny7s2fhiy+U4gBl2vv11/CYxCrSaJ4kHFEeBW0c2ymlPCmlHAzUywyB\nHMxh/h8gEOgH1BFCdMgMGTQZw8fHh61bt9K7d28aN27MqFGj+Ptv+3R7nToQFgZ9+8KKFclOVqoE\nO3eqdZANG+Djj+Hnn1W8rM12JZPUaDSZhCPKo2jyA1LKMKuPmRWqxJEc5l2llE2llP2BPVLK5Zkk\ngyaDuLi40K9fPyIjIzly5AjVq1dnu51h2f39Yf166N8fQkOTnfT2VsOSOXNUGN+EBBg0SA1XAgMh\nE0JNazSa9HEkDe1vQogPpJTTUjmfKVkFHc1hblWvZ3pt63zZWY/JZOLrr78mLCyMTp068eqrr/Lx\nxx9TuHDhNHOYlygBCxe60r37s1y5coM2bZKNXCpXhrVryR8ayt9NmmA0aUKhOXMoWLMmd3r14vbA\ngRjmkNQajcY+nJXDPASlQO5IKedYnxBCmADnhG7UOcyfCN555x0CAwMZPnw4TZs2pXfvEC5frsa3\n36buxFS6NGzcCM2aeeDuDt262ShknYt50iQYNAi3YcNwa9xYfe7cWQXl06SLzmGucSSHud3KQ0pp\nCCGCgK1CiB7AbGA/kAcYat53Br8DrYBQncP88eaZZ55hzpw5bNu2jbfeGsSlS2upX/8evXo9m2qd\natUSFYgKedWjRzqdlCsH//mPWjAZO1aFPJkxA6pXz9yL0Wiechy1tooB6gLRwPeoB3k4cA9IbTor\no+gc5k8YDRs25MiRHbRosZQ334znq68WJclPnZwqVZQC+egjWLDAzk4WLYJChaBBA3jtNXj3XXC2\n05RG8xThyLQVAFLKK0A3IcRA4HnggpQyUxcTpJTnAUsOc6B/ZravyX7y5cvHjBmBFCli4rPPnmPV\nqsbMnfttqolrfH2VQVXTpmoE8vbb6XSwbJnKjfvpp/DKK8qc19dXfe7fP22nE02OJDGdKsD9+/ez\n1G/Duu/s6D8n4rDySERK+RewNxNl0TyFzJpVnMjIYpQoMYb69evzwQcf8OGHH9r8YQqhFEiTJkqB\n9O2bRsO5csE770CnTso35NtvVc6QNWuUpdb06UoTaXIkV69e5aeffsLDwwM3NzdcXFy4ffs2Xbt2\nZcuWLVSvXp08efIwc+ZMhBCcPn2ad999Fykl9+7dw8/Pzyl9A5b+V69eTXBwMEWKFKFWrVpIKQkM\nDLTkGE9OZsiWk9CBETXZSu7csHSpicmTGxMREcHvv/9OrVq12LvX9ntJpUqwZQtMmKBSf6RL4cJK\neYSHq+mrX39Vld9+G0aN0rGyciAXLlxgxIgR9O7dm+7du9OmTRuOHTtGpUqVuHLlCnFxcRQpUoRd\nu3YBKkVtfHw84eHhCCEIDw9/5KRlafUNJOk/MDCQ8uXLU6tWLT755BO6dOnCiBEjOHXqlM22hRBE\nREQ8ekK1HIZWHppsx9MTKlaEChUqsHbtWj766CPeeOMNhgwZwm0b4Xi9vJQCmThROZjbxXPPQb16\nyvKqbVvYt08pkh49wEZqU032MXz4cPr165ck7Lq/vz9+fn4sX76cpuYRY0REhCXHha+vryVneEBA\nAOvXr8/Uvl944QWAJP0nx8PDg4SEBM6ePZtq+xmRLaehlYcmR2EymejWrRtHjhzh2rVrVK1a1WYI\nbk9P2LoVpk5VM1AOU7y4mgO7e1eZcukUuDmCyMhIbt68mSQ5EqjRRZ48eYiNjbUkc/rrr7/Ib/bl\nKVCggCWKrI+PD5GRkZnad27zGpl1/8k5c+YMbm5u1KhRI9U+hBA2ZYuPj2fkyJH06NGDoKAgSzKn\nFStW0L17d4KCgggNDbXs79ixg2HDhtGlSxc6d+5MWFhYivKJytRZOKQ8hBC/CCGcHyhe89RTtGhR\nfvjhB7799lsGDBhAUFAQV65cSVKmYkWlQGbMUErEYfLlgyVLVEyUgACIisoM0Z8YxoxRA7Xk25gx\nmVPeFgcOHLCknAU1uhg/fjzjxo3jxo0b3LMaJSYkJFhSylrvPypp9X3z5k2AJP0ncuLECWbNmsW0\nadMIDg6mWLFilnObNm3i8uXL6fYdEhKCh4cHCxcuZObMmYwdO9Zyzt3dncWLF5MrVy7L/vnz5/Hw\n8GDJkiXMnz+f6dOnc/369STl69at+8j3wh4cXTCvB9x3hiAaTSKGoQYEBQpAs2bNOHz4MKNHj6Zq\n1apMmTKF7t27YzI7/pUvrxRIokHVsGEOdPTuu1C/vspW6OmpzHp//llNb2kYM8axB7+j5W1hMpks\nucNBTRlt3LgRf39/3N3dk6wXFC1a1JIK+fbt2xQpUsRyLnmK5Lt377Jhw4YkxwxzetrmzZun23fi\nNFZipGhrfHx8eO+999i7dy9z5syxPLSvXr3KihUrLLnIU5MN4OTJk0RERHDw4EEMw+DBgwfcuKEC\na1jnN0/cj4qKIiAgAICCBQvi5eVFdHR0ivLOxFHlMRX4QQjxf8BZ4K71SSll9mdl1zz2zJsH69ap\nJFEmk/pxTJkyha5du/LOO++wcOFC5syZY/mRlCv3UIE8eAAjRtjZUceOKktVjx7KfLdCBRX6d9Ys\ndU6T5TRs2JCRI0diGAYmkwnDMNi1axe9e/cGwNX14SPL39+fw4cP06hRIw4dOkQ9K6WffBSSP39+\n2rZt63Dfu3fvtvRtq11rSpUqxbp16zAMg7/++ouiRYvi45MyX56tNjw9PSlVqhR9+/bl3r17zJ49\n26LIrE2EE/e9vLwIDw+nadOm3L59m1OnTlG2bFmioqKSlHcmjvbyGcr/YhmwDzhitWnPb02m0LOn\nmkGaPz/pcX9/f/bu3UuTJk2oXbs2X331leVNtEwZpUB++EGlQreLpk0hf/6HCURatFCL6EOHqtAm\n2hIry/H09KRnz55MmDCBlStXsnbtWvr27UuJEiUAkmQBrFu3LteuXWP9+vWYTCYaNGhgOZfauoSj\nfffp08fSd/L+Q0NDiYmJISIigo0bN9K2bVtKlizJpEmT2LNnT6r92JKtc+fOREVF0aNHD7p27Zpu\nKKVOnTpx/fp1unXrRq9evRg4cCAeHh4OX3OGSC1LlK3N29u7QlqbI21l9fakZ0l7HLHOJJicI0cM\no2hRw5DS9vlTp04ZjRs3Nvz9/Y0DBw5Yjv/xh2H4+hrG6NGGkZBghxBLlxpGvXpJC8fEqKyF/foZ\nxj//2HcxTwCPw2/ku+++M27cuJFmmfPnzxshISHZ1n8iV69eNYYOHWqsXLkyS2TLKM7MJIiU8rzZ\n+/t5oA0qdLoPcNF8XKPJFKpUUfPn3bvDP/+kPP/888+zadMmBgwYQLNmzfjoo4+4e/cuJUuqEUho\nqJ1uHB06wOXLKm96ImXLwvbtcP48tGoF5sVSTfbTsWNHm9Z31mzdupXWrVtnW/+JPPvss3z11VdJ\nnAadKVtW46i1VUkhxC5gHTDQvK0B9gshijtBPntkekEI8ZsQ4nshRKPskEHjHN57D4oVUw7htjCZ\nTLz55pscOnSIs2fP4ufnx+bNmyleXPmBrFoFI0emo0By5YIvv4Tki5iFCytv9IoV4aWX4MKFzLos\nTQZwc3PDy8sr1WivMTEx+Pj4JJleysr+08LZsmU1ji6YTwceAM9JKS8CCCHKAj+hFtO7Z654dlEH\n+AOIB45mQ/8aJ2EyweLFyuoqLUqWLMnSpUtZs2YNvXv35tVXX2Xy5Mls3uzBq68qK6yJE9OIzB4Y\naMo9hTQAACAASURBVPu4q6tyY58yRVlgrV6tshZqspXkfhjWlChRgnLlymVb/2mRFbJlJY4umL8G\nDE5UHABSyguokOwtMksoR3KYA9uBPsBE4MPMkkGTM3jmGbA3/lzr1q05cuQIBQoUoGrVqmzevIyN\nGw02bYJ///sR179NJvjwQ5g2TTkThoWlX0eTbeTkYIU5WbZHwVHlcRew9RNMADLmoWPGwRzmPwHV\nzX1fzywZNI8vhQsXJjg4mNDQUMaOHcubb7ZhwYILbNsGH3yQAQOqwEA18nj7bWXKq9E85TiqPDYA\nU4UQFts1IURJVJ6NDanWcgxHcph3A84DwaiRx4xMkkHzmBMQEMD+/fupVasWjRtXp1OnuezaZTB4\ncAYUSL16amF9xgw1lEkjB4lG86Tj6JrHh8Bm4LwQItG6qgJwCLCVJNRhHM1hLqXcBeyyp22dwzxn\nkVYO89S4ds3E3bsmSpe278H9zjvv0LBhQz788EPi4xexefNqevd2ZcKEG9j0pbp3D1N8PEbBgrYb\nzJ8f0/LleLzzDgmtW3N9xgydK13zxOCsHOZIKS8LIV4AXgd8gb+B41LKjY8kqX3oHOZPKJcuXXL4\nf7Jihcoyu3WrWs+2h9KlS7Nnzx7mzJnDyJEvsGXLNkaNKkNwcC5S6Ih//UtlILSKLWSjQSXA22+T\nv1s3NZ1l5Uj2uKJzmGscyWHucGBEwFtK+YuUcrKUMtjJigNUDvMW5v51DvOnnP79VTzDL790rJ6L\niwv9+/fn0KEdVKv2IaGhv1Ky5N8MHRqPOSSQol8/ZWEVF5d2g3nzwsKF0Ly5ms46ftzha9FoHmcc\nXfPIjsCIOoe5xoKLiwpB8vXXkEYEiFQpW7Ysv/yyhKVLoXr1d5k9ex4+Pndo2fIO27eDUckbGjZU\nAbbSw2SCceOUN+LLLyvnEo3mKSFHBkbUOcw1aVGmjDJ4CgqCAwfAzS39OtaYTCZef/11Xn/9dU6c\nOMG0aaNYuNCFbduGUKqUG191GUarKYGY3nvPvlznvXur8L6dOyufkJ49H+m6NGmjc5jnLBxVHokh\n5wKsjhmAyfxXm8pqsoT27VXk3bAw9cx+VHx8fJgzZwqTJ9/k++9/YNKkQ3Sb9ia/3PXkcvuFvDT3\nLUqWtKOhV15R6yAtW8KZMzB6dBpeiZr00DnMcz6OTlutA6oCz1ltnlZ/NZosY86cjCkOawoXLsz7\n7w8iJmYOoaG3WOLnwbwNJ6hY8Q7t2t0mPNyORipXht27lUbr1Uunt31EdA7zxwNHlUdd4O/EAInJ\nN2cIqNGkhjPSFri4uNC8eXO+2f8zs4734+23P2fDhsk0bHiZKlWus2SJYTNQo4USJdQI5PZttZh+\n7VrmC/mEo3OYPx44+vObCvwohAgUQvgLISpbb84QUKPJLry8vJg5czxXrgxjypSV3Lo1hrfeiqB4\n8duMGXOP2NhUKhYoACEhULOmssQ6cyZL5X6c0TnMn9Ac5uhkUJqnkIIFC/Lee305f34aa9feombN\nMXzxxXLKlLlDp043OHLERqVcuVRi9UGDVKpbJ/+QnUI2JDF/knKYHz58mJ07dzJ37tx0+34acphn\nTXJcjeYR2LlTvfRXr+6c9k0mE40bN6Zx48ZER0czZcpUvvsuF6tX98XHx2DMGA9at3YhyTNswACV\n3rZ1a5g9W+UPeVzIhiTmT1IO86NHj9K+fXt27NjBnTt3KGAOD/2k5DB/1GRQBYCawFVzG9F6zUOT\n3Zw7B127pkzNkWGuXYNk4RrKly/PjBmfcPXqEIKDf+H69Wl06XKYEiWu88UXd7lhHVCnVSvYsEHl\nS58yRae3TYOGDRsSGRmJYb5HhjmHeeKIIHkOcyklAIcOHaK61VtDajnMrbd27dpZFEdqfe/evTvJ\naCS9HOYHDhzAMAxiY2Pp0qULrq7/396dh9lcvg8cf58ZY1C2GFshKXf2ZZBkTXbZkpAle4kW0vL1\n/Ulp+xYiVEqpKUq0iEoRlUQyqGxPIpIZhDAau/P74zk0pmHmM3OWz5y5X9d1rrPOOfdcn2vOPc/n\neZ77znUuSV3sPa666iratWtHXFwc06dPp1WrVhnqYQ6c18M89esDyekO8wIi8im2b8YcoDi2x8da\nEdHaHyqkevSA2Fhbs9Cvpk2Dhx9O86m8efMycGAffvvtCRYvPkKNGuN59NHPKF78b3r23M+5hTe1\nasGKFRAXZ7tchcmKG38Ltx7mCxcuZNCgQZxMscoiXHqYe7wO/gsSkVexS3J7AQaojm0ONRNIMMZ0\nDUSQ/hAfH++NjY0NdRgqhczUtkrPoUP2tNULL9gzRX5x8CCULw9r1thTUOlISEjg2WdnMX16Lk6c\n6EO1ascZOzaGVq0i8SQdhltvtXMis2c73+EYQPHx8bj9b+T111+nS5cu562GSu33339n1apVdLlQ\nk68Af/5Z8+fPZ+XKlURERPDYY48RERER0NiyKj4+HmMMe/fu5b777jv3WGxsbJoblpyOb9oCI40x\n58bwvtNVQ4G0168pFUQFC9qSUwMHwu7dfnrTQoVsH48JEzL08lKlSjFx4gPs338X06Z9zoEDr9Gh\nwy+UKLGPca9E8PfsBbZPesOGsGtX+m+ozslOPcxvvvlmnnzyScaOHXvuVFKO7WEOXEqqkiQp3ic4\nJ9qUSkeDBrY47q+/+vFN77vPZiXfctCMiI6Opm/fbmzbNopvvz1C9eqv8sgjX1OkxDHuOD6Wv1p1\nh3r14Mcf/RhoeNMe5u7hdLXVQmC0iPTy3feKSAwwDljk18gySEQqAvdiT59N9Vd9LZW9jfR3Q+JS\npexKqSlTMrWiqG7dOnzxRR327t3LM8+8ySuv5Gbm8Tt58Mr8PNr4JqJmxeFp09rPQYcn7WHuDk5H\nC8OAMsAB7IqrxcDv2IZN9/o3tAy7C9iF/V22hygGlROMGmVrV2VBsWLFmDBhKH/91ZfXXlvCbO8J\nmiW9yJ8392Hpbc9z7JifYs2h3Fys0M2xZYbTZlB7gBtEpClQyffzm4BFvuq3fiEi1wHPGGOaiogH\neBE7OX8MGGCMSblltywwGogF+gAv+SsOpc5z5ZX24gdRUVH07t2J3r1hzZp1PHTvQB55bzyfzPmG\njRUaUG14R1rfUZ4w+75RYcTpaSsAjDFLgYA0LxCRkdjVXEd8D3UEoo0x9X1JZQLQUUQeB64B/gSS\nsaMhLWOqsp1atWowY1kNDmztR4EHn+L6xTOpPXgMGweX46eStYnp0ZKbHu1AVH7ny0+VCpRMJY8A\n+xXoBLzlu98AO9eCMeZ7Eantuz0aQERigVexTaoevNgbaw9zd8lMD/PMeuutfJQte5pGjVxc6TZv\nXipPtl0P9h8+zOapczg2dw2XTXiCo+MHsOLSKiTfUJcqw9sTUflaLfmu/C5gPcyDwRjzoYikXExf\nAEi5X/eUiJzrY26MiceerkqX9jB3l0Ds87iQ2rVtlfR166Bo0aB8ZNaUKkX5yf8Hk23dps/f+ZYf\n/vcFVyzeTMXP7yAy13ESq9RBht3KJR1uhiJFsvyR2sNcOelh7rrkkYbDQMqdVOcSh1IZ1ayZLV0y\nYAB8+KEf/mn//Xe7eTAIjX0iIiJofXsjWt/eCIBvl/3BrDHLuOS7jTTuP5PGA4ZwoNjlFOzSmkK3\n3WqX/2akA2IquXPnJj4+3t/hq2wirbpdF5Op5CEiuYAoUs0xGGP8XVUIYDnQDpgrIvXQ6r0qk554\nwn6vTp9uNxFmyXffwdSpsGyZX2JzokHDK2jwpe1qt3ZtEkP/9wsHPtnCDVO/o/XLfbgmYjfH69Wh\n4G1diWjZ0u6Oz0C2rFq1KgkJCbz33nsUK1Ys0L+GciApKYn8LqpGAA6Th+/L+yXgQv9uBaIN7YdA\ncxFZ7rvfNwCfoXKA6GiYNctu7G7UCESy8GZdutilu8uX25LrIVKzZn7efDcWiGXjxi5MmDCcbz88\nQu3l39JqRRxtIh8hulA+om++mVxt2th2uSkqx6YlJiaGvXv3BucXUBmSnJzM0aNp7c/2v7Ml5dPj\ntLbVSuzE9HPY00nnMcZ8neE3CzKtbeU+wZzzSGnePKhZE8qUyeIbvfSSbaT+8cd+icufNm2Cl17a\nx3uzvZQ8sJHm3kl0KbCWGsmJUK0audu1gxYtoE4dUtaQD9UxURcXquNysdpWTpPHEeB6Y0y2O3Wk\nycN9sv0X1dGjUK4cfPklVK4c6mguaMMGeOONZGbOPMnxg4e47uQz3FZoGa0j93HZ0WQiW7TA07Yt\ndOhAQorJUuUebkweTneYbwZKZj0kpcJA3rxwzz3w7LOhjuSiKleG557Lx65dBfnq+zLUGDmZUdEr\nueboz1Q+8xD3L97CT088wanSpSncqxe8+aZdDKDURTgdefQBHgOmYPdjnEj5vDHmU79G50c68nAf\nN408jh2Dhx6yNbF8PXUy5tAhWyyxfPmAxRYIXi/89BPMnu3l7bdPcOTI3xTwvEuL469wd7EkKu/d\nQ2Tjxni6doUOHWxlYRUy4TDymIGtbfUs8AGwIMVlflaCVCqUvF4oUMD2Ann6aUijVXXaChbMdokD\n7OKr6tXhqac87NgRzZIll9H77iEsillOo/1ruSb6BQZ9/SdrH32MU5dfzpm2bW0jKx2RKB9HI4/s\nTEce7uOmkcdZ27bB/ffbeYKJE20H2ZwkISGBw4dL8cEHMGvWcbZvP0VM9DxuPDqZu4smUnX/n0Q2\naULEbbdB+/Y6IgkSN448MpU8RKQZUBk7ctkEfGmMcXVfTU0e7uPG5HHWwoW25fiMGVC/fqijCZ7U\nx2T7dvjgA3jnnWNs2ADF835K078nMKToTqr/tY+IJk2I7NbNJpJ0lgCrzMv2yUNESgAfAbWw5c8j\nsKexNgM3GWNcuzhck4f7uDl5AJw8Cbly5awSUhc7JomJdnf+O+8cY/VqDyXyfUmTpGe5K+Y3ah7c\nh6dJE3J17w4dO8KllwY58vDmxuThdM5jEnAKKGeMqWCMuRrb0/wgttqtUmEjKsph4jDGDlnCVMmS\nMGQILFuWh507oxk9rg1/NF5I4wNbqJDnYwYuOsGqBx7hRPHinL79dli8GByWvFDZh9Pk0Qq4xxhz\nrvGyMeYPYDjQxp+BKeVWcXG2wOK/HDoEgwfbIUuYK1oU+vaFRYvysGdPbp6a0ox9LRfQ+NBWquZ+\nn4fm/M6223pwrEQJzjz4oJ1EUmHFafI4CqR1nusMgSlNopTrnDkDLVva/8L370/xRN26duXVu++G\nLLZQKFDAFp2cPz+afftyM3ZaKzY0+4LKyTtpfHQck174kr+uu56/K1XCO3EiaOmTsOA0eXwOTBCR\n4mcf8M2DjPc9F3Qicq+IzBCRb0XkzlDEoHKWO+6AzZttVY9KlWyVknNnZx56yG4azCGrGFO75BLo\n2hU++ywve/ZEM+zlPsyr9xWXn0zktm29mPufaRwtU5YjN94Ic+agfXezL6fJYyQQA+wQESMiBjtx\nfgkh6mFujJkEDALWG2NeDkUMKucpXBgmT4ZFi2D2bBgxwvdEixZ2lv1T1+6XDZoCBaBnT/jqq/zs\n2n0JnaY+zPhqK7j8zA7uXVaF5X0e4GiRIiT36mULTObQhJtdOe1hvldEqgOtgYrYnuKbjDGL/RmU\nwx7mAN2xmxaVCqpq1WDpUjhytmmyx2NHH5MmQdu2IY3NTQoXhv79PfTvX4g//4Q5c55nyLQxJG3a\nR7dZYxgwpwNFCkQS3b8feQYOhKuuCnXIKh2O+3kYY07zz65yv3PYw/xq4G6goTEmqx0alMoUjwfO\na7XQpQs0bx6yeNwuJgaGDIlkyJDLSEy8jHfeeYM20w5SaNsmbn/2SXqOn8yp8mUoNGwoUT166EZE\nl0o3eYhIG2CRMeak7/YF+am2laMe5r4YdbJeuUeuXFCkCFu32tW7bXQd4gWVLAnDh+dm+PBi7NxZ\njDfeqEPd1/6iypav6TlsAi3vHcGRG+oRM/x+Ilu3zlSHRBUYGRl5LABKAHu5+GjDix9WXDntYe77\nmX4Zee+MNHVXwZOUlBTWx2Tz5ijuuacwzz9/ijFjDlGunPv3PITymERGQv/+9vTWjh3NmDmzEf95\nL4lmyz6h17LhVMx1OwdbNyPfnXdyqmrVHLV7041/K+kmD2NMRFq3g8hvPczdvJs5J3L7DvOsKlXK\nnr2aODEXHTrkYdAg+M9/3L352i3HpFQpuP56YAoYU5EpU/qwauZW2n08m97zexFVIIpc3TtRbMhd\ntuZ8mCeSUB2XxMTECz7nKBmIyBIR+dcJSBGJEZH4TMSWEcvxbUDUHuYqu8mdGx580JY/37kTatXK\nEXsI/UoEJk+OYeX+enRYN4HHeq2h36lXmPXycf6o3pjEQiXZPWAg3rVrdcVWEGVkzqMJUMl3tzEw\nWESSUr2sIhCoutTaw1xlXxs2wLJllLrzTt56CxIS9LR9Znk8UK2ah9ffKI13Rml++qkNY154lB1z\n19Pi9QV0mdGSfPm8nO7YhpLD7sZTp07Yj0hCKd3CiCJSFZgHeICywB9AypO3XuzKqEnGmNcCFGeW\naWFE93HLKZKA2rULqlaFLVugSJFQR5Ou7HpMNm70MmVyAmb2Bm46+BldeJdCeU9wvF0LSg27m4j6\n9SEiFGfd/cONhRGdVtVdCnQ2xvzlr+CCRZOH+2TXLyrH+veHsmVh9Og0n/Z6YcECuyorMsTrBsPh\nmGzZ4uXFqYn8OGsTTfZ9QRdmUTw6ieTWN1Fq8EAimzSB6OhQh+lItk8eFyIiuYFYY8yKLL9ZgGjy\ncJ9w+KLKkM2boXFj+O03yJfvX08nJdn9hEeOwH332cnimJh/rgMuKclOxFx2Wdgdk+3b4aWX9vBD\n3AYa7FlCG96jSsRO9letQtGe3bikSxeb2F0u2ycPEakLTAOq8O/Jdq8xxvGmw2DR5OE+4fZFdVGd\nOkGzZjB0aJpPe722zMlHH8Gff9ragbVqwZtv/vu169fb18bEQLFi9jomBi6/PBNnxubNg7vuguRk\naNeOfbfeStH27cNyruCPPyAu7gCfvbWDK39ZR0vvO7TyrOBUoUuJbNeCor164WnY0JWjknBIHiux\n8xsvAW8CA4HSwChs2ZA5WQ83MDR5uE+OSh4rV9phxcqVWX4rY2zy2LvXJpqzyaZhQ1ukMbW1a21S\nOptkziacy4udpHC/TracSuXK8MYbnJoyhVz9+sF//5vlON0sORk+/fQ4M15LJOmrbTQ78TFt+JhK\nEYkcrlmdIj17kLtDB9eMSsIheRwFrjPG/CQi3wCPG2MWi8gdwEBjzA1+iTgANHm4T45KHgAnTti1\nu0G2fj3MnftPkjl73bIlPP/8+a9N+OMPShUubMvj5hBeL6xZ42XGjD/55oO9VN/zPa28cbSOiOdM\nkQJEtW9LwW7doEGDkI1K3Jg8nJ5mOoXdtAfwC1ANWAwsBSZmOkKlcoIQJA6AKlXsJUMiIlj03SUU\nLQo1a6Z6buFCaNrUlad1ssLjgdhYD7GxxWBKMRITqzBnTndaxe0nct12Wr0eR7s3+lHRs5fkpo0o\nMmqUPb0Vhqf2nHC6du17YIiIRAA/Yqvrgp0DOeHPwJRSfrRwYYrSvxe3daudoqlZ05ad378fOHoU\nxo2DMmXg4Yft5H+YKlkS7rknH6tWl2bpkYbU/ngaT3dYReXo9Ty1uA5bm3Xhz5KlOTp+PBw+nP4b\nhimnyeNh7Ca94djChdeKyDbgPWCWn2NTSmXVX39Bnz627eHvv2foR+68E7Zts7li5UrbHLFrn7wk\nf7wYli2zp9/q1LFri7/8MsC/QGjlyQNt20Ywd25xth0uT4clj3NX47Xcvm8an42cS1KREuy5pYst\nIZDDOE0eY4BmQJwx5iBQBxgH9APu929oSqks+eQTu0Exf3775VapUvo/4xMRYReHzZxpl7t27Ohb\nZVyhAkyYYGutdO1qr3MIjwcaN45g0eLLeW9fW355aiFNinzLyx+WJaFmY3ZeKRydPh2OHw91qMHh\n9XozfKlQocL+ChUqXOXkZ9xyWb16tVe5y65du0IdQmhs2uT13n9/4N7/5Emvt08fr7dcOa93yRJH\nP+r0mOzf7/UeOeLoR8LO6tWnvR3b/uq9NfJt7yJPTe+B3Pm9u3rd4fVu2+a3zwjV34rvezPN71Sn\nI48JQJyIdBGRWBGplPISiOSWHhGpLSJTReQtX5dDpdytbFmYNQs2bgzM++fKBTfeaEcbTZsG5jN8\n5s2DK66AAQPgu+98dQlPn4ZXXskx/4HHxkbw4YLyvJl0O5smfcktMfOY/VY0+8tX4xeJ5ejcuSma\n3IcPp8ljLFAfO8fxA7AeW+X27HUoxGILM14O5JwxtMq+8ua1mwWfey5wn9G7d1Bqv/fta3NghQrQ\nr59d1fX53CTbw71aNVjs1w7VrpY3LwwbVpglfzSl9aap3NtxBeO23cGGW/+PfXmLsKVRS47NmhU2\nk+xOl+qWC0gUqTjsYb4GeA24EWgHxAUjRqWy5O677Uz0zp1QunSoo8mSkiVt2fmRI+2irhEjC1F5\n4UdcsXY+DBwI9erB+PG23koOce21kbz9QRVOnqxCXFx3Hh+/hquXr6XVshe4IaIfe8peTZGenSnY\no4etOZ8Nl/06Sh7GmB2BCuQshz3MrwEKA22BffxTOl4pdytc2P7b/vzzdgI6Mw4dst/YQ4fa//JD\nzOOB1q2hVSvfd+EVN9tZ96eesvF99ZWDDSfhISoK+vcvSv/+LThxogULFtxNt8m/kue79TQfu4C2\nT95A1CWReFvfSLE7+uBp2tQu8coGHCUPETmDLcGelhPALuAdYIwxJrMn+Rz1MBeR9r7XHgdGZvIz\nlQq++++H226DM2eclwv//HP7X32bNlAuKCcEMuy8f6Lz5YMnnrCJ0mVxBlvu3NC586V07lwDr7cG\n33/fjUem7mTHgq00mPMVbebeT03PDlsepVcPojp1cvWo1Gl5koHYeY9HgbMVdOsCjwGvYuc+xgDz\njDGjMhuUr4f5O77RxqvAXGPM577ntgNXOW1FGx8f7y1ZsmRmQ1IBkJSURP78+dN/YTjzeh2dsvAk\nJVHg8ceJ/vprDo0bx/FGjfwaTiCPyYsvXkq7dkcpUyb8Jo+zavv2CGbN+ptvPjpM1YR1tCGOthGr\nSL6iBLm7d+VQ8+bkvfbaoMeVmJjot/IkDwD9jTGfpHjsJxHZiW0Gda2IJALvYosl+oP2MA9TOa62\nVVadOWPPCdWrBxs3UqRAAb9/RKCOyZkztlxWu3YFGDbMzpHkzYstV1+hQrZu1OQPpUpB/frAFNi/\nvyazZ3eh+St7ifl5PV3+N5XOz73AqSsvp2C/O4jq3j1oozi/9TDHrmjalsbjfwBlUty+zOH7Xoz2\nMFcK7BfsF1/YZbABSByBFBEBjzwCa9bAzz/bIr4LFgAjRsANN8C6daEO0TWKFIEhQ/Lzw7ryzNzX\ngSMT36d+ubX0/G0cb/53CQcrVOXA1ddw6sknbYXLEHGaPL4FnhWRc8lBRIoAz2C/5AE6Y4sm+suH\nwHFfD/Px6E52lZOVKBHqCLKkTBlb4ffll23eeLPLfLtJpGVLW7I+TJax+kvhwnDPPflZ/PWlvLW7\nPQee+YDrr9xAt9/GM2P0EuaPeD79NwkQp3MepYHPgPLAdmzyKQNswiaNa7Ff9l1SndoKOS3J7j56\n2uoijhyx53mCvIQzmMfk+HF7OitvXmDfPltw8bPPYMYMaNEiKDFkF6mPS0ICvPbaISpX9tK5c6GA\nfe7FSrI7GnkYY3Zi91t0AKYDU4G2xphaxpjtQDxQ2m2JQynX277d/vft9dpig1Wq2KWtYSw62pc4\nAIoWhenTYc6cbD+6CoZSpeD//q9gQBNHehy3jTXGnPZNkF8CfAEUExGPMcZrjPnT7xEqlROULQu7\nd9ult+vX23mNAJcWcaX69Vm1Copug6uuCnUw6mIcjTxEpICIfApsAOYAxYFJwFoR0fMPSmWWxwPP\nPGO/MX/+2e62y6F+/hnq1oUxY2wbEcCWgVeu4nTCfDwQDVwBnD2sw7C7wbWToFJZ0bo1TJ0KhUJ3\nKsIN+ve3fdc3bLCrsubPB0aPtnXhdwS8yIXKIKfJoy0w0hiTcPYBX8mSocBN/gxMKZVzlS5tpz+m\nTbMVWIYdeMw2oIqNhaef1pGICzhNHpfyz4gj9fvk7F0+Sim/a97cVpYfOiIaRo2CH36wtd8rVoRn\nnw11eDma0y/8hcBoETk70e4VkRhsN8FFfo1MKaWwNaFEfHfKlbPnsWbNylFVet3IafIYBpQFDgD5\ngMXA70BB4D7/hqaUUhdw3XUkd+7J4MG2RuTJkymeW7fOrlhzsIdNOed0n8ceY0x97D6Pe4AXfLdr\nG2N2BSA+pZRK0+nTcO21dlVWyZJ2o/oXX8Cp1WvtkufKleHRRzWRBIjTpbpLRKSQMWapMWaqMWaS\nMeYLoKiIxAcoRqWU+pf8+W1V+xUrID7eToOMHg0Dl/e1my5ff93u1D+bSLZsCXXIYSXdTYIi0oR/\nmiw1BgaLSFKql1XElixRSqmgK1vW1soaMcJ3CisiwlYfrlcPxo2DVatsYa20PPAAXHaZ3WNz9lKk\nSLbs7hdMGdlhvh9bit3ju9wNpCzI78Xu8xjh9+gyQERuw7af3Qf81xjzdyjiUEq5Q1RUqgc8Hrju\nOgYNgsREuPpquOYae3311VC2QkUit22BDz6Abdvs5dQpW7E2m3T1C4V0k4cx5mfgKgARWQp0Nsb8\nFejAHLgZ6A3UAPpg+50rpdR5Ro2ymw9//RV+/BHef9/enj27P/UGpXrxwYMkHMhDTEyqZHTgAHzz\njd2wmMM57WHeFMC3VDcKOxJJ+XyyP4Ly9Sp/xhjTVEQ82IRQHTgGDDDGpOwpMgVbpPF3zh8RKaXU\nOWXL2kuGFCrE7Z3slpIrrvhntNK79m7qPvUgvPsuTJ4MMTEBjdnNnE6YXycia7H9wo8ASb7Le1q9\ngwAACQtJREFU2dtZJiIjsS1to30PdQSifau8HgEm+F73uIjMAkoAA4Bl2ASilFJZtnQpJCXZKvH3\n3munQgZNrMTf3/1oM0rVqnYbfA7ltKruJGyS6IhtDxsIvwKdgLd89xtgNydijPleRGr7bo8GEJGm\nwBvAKWDwxd44ISHhYk+rIEtKStJj4jJ6TP7t0kuhRg176doVDp2AQ8OHE9WoEYVGjOBUXBx/vfgi\nREYGLAY3HhenyaMKcL1vHiQgjDEfikjKwWUB4FCK+6dE5Fwfc2PMUmBpRt5bGw+5izaDch89Jg60\nbw8tWhD1+efkLV06oB8VquPizx7mm4GSWYrGucNA/hT3zyUOpZQKqTx5OHNzB07nwNlWpyOPycAr\nIjIFe3rpvNKWxphP/RVYCsuxS3Hnikg9IGCjHqWUcmrmTLsfcdYsu9M9p3CaPGb4rtMqZ+kFAnHS\n70OguYgs993vG4DPUEqpTOnRA377zVaLj4uDm27CVv99+mnbnyVMM4rTpbpBKbvu6xFS33fbC9wV\njM9VSimnIiNtWZQGDaBnTxg4EEY/VI3IypWhenW7w71Xr7DbsZ6R8iRtgEXGmJO+2xfiNcZ85r/Q\nlFIq+7jxRlizBm6/HR76O5px48ZCp07Qty+MHWtLoLz3XtqbTSZOtHW48uX751K0KDRpYm+7UEZG\nHguweyn2+m5fSKBOWymlVLZQooSt7Hv47EaGWrVg9WrYvBmSk6FYsbR/MFcu27B9/357nZwMu3fb\nn08reRw6BAULBuz3yIiMlCeJSOu2Ukqpf4uMhMKFUzwQFWU3FF7M0KEZ/4AzZ+x298hIW2flQgkp\nwJxOmCullMqEOXNgyhQoVMheCha01w0b2na7GRYRAXv2wM6dIS2PoslDKaWCoGFD+11/8KC9HDpk\nr48dS/v1EyfC229Do0ZQuXIe2rdPkSs8nguXmA8STR5KKRUEJUrYS0bdeSfUrg3LlsHMmfkYMcKu\n+h0/3va3CjVNHkop5UJ58tjlvw0aQJ8+ByhevBQ//eSeQr6aPJRSKhuIjISaNUMdxT909ZRSSinH\nNHkopZRyLFuetvL18OhhjBkoItdj+3h4gXuNMYHqM6KUUson2408RKQ8UIt/Og0O8l1eA7qFKi6l\nlMpJXJE8fO1tl/pue0TkJRH5TkSWiMhVKV9rjNlqjBmf4qFIY8wJYDe2jIpSSqkAC3nycNqzXEQK\npXqLv0UkN7ZJ1e4gha2UUjmaG+Y8HPUsT8OrwDTs73LRHuZKKaX8I+TJw2nP8hQ/19t3vYYMNohy\nWwP5nC4pKUmPicvoMXEnNx6XkCePNASsZ3koGsirC0tISNBj4jJ6TNwpVMclMTHxgs+FfM4jDcuB\nNgDas1wppdzJjSMP7VmulFIu54rkoT3LlVIqe3HjaSullFIup8lDKaWUY5o8lFJKOabJQymllGOa\nPJRSSjmmyUMppZRjmjyUUko5pslDKaWUY5o8lFJKOabJQymllGPZNnmISFMRefVC95VSSgVOtkwe\nqfuYp9HXXCmlVAC5ojAi2D7mwDPGmKYi4gFeBKoDx4ABxphtZ19rjNkKjBeRuLTuK6WUCixXjDyy\n0Mfck+qtUt9XSikVAG4ZeWS2j7k3nftKKaUCwBXJI6t9zC90P7X4+Pgsx6r862JtLlVo6DFxJ7cd\nF1ckjzT4vY95bGysntJSSik/ccWcRxq0j7lSSrmYW0ce2sdcKaVczOP16hyzUkopZ9x62koppZSL\nufW0VUCJyI1AHyAvMNYYo3MqISYitYARwAngQWPMnyEOSfmISDHgE2NMnVDHokBEqgMvANuAN4wx\nX4cijpw68shrjOkDPAW0CHUwCrAbRO8CPgWuD3Es6nwjge2hDkKdUxdIBE4BG0IVRNiNPDJS5sQY\n84mI5AOGAQ+FMNwcIYPHZIWIXI8dfdwawnBzjIwcFxG5E5iJPS4qwDJYpulb4F2gODaxh+Q7LKxG\nHg7KnBQBJgGjjTH7QhFrTuHgmNQGVmOXaN8TglBzlIweF6A5MBioKyK3BD3QHMTBMakBRAIHfdch\nEVbJg3/KnJx1XpkTINb3+ASgFPC0iHQOaoQ5T0aPSQHgdWxSnxPMAHOo9I7L2ZJAtxhj7gK+N8a8\nH/Qoc5aM/q1sByYD/8POfYREWJ22ykCZk9O+Mid9ghxajuXgmCwBlgQ3upzLaUmg9Er/qKxz8Ley\nAlgR3Oj+LdxGHqn5vcyJyjI9Ju6kx8V9XH1Mwj15aJkT99Fj4k56XNzH1cckrE5bpUHLnLiPHhN3\n0uPiPq4+JlqeRCmllGPhftpKKaVUAGjyUEop5ZgmD6WUUo5p8lBKKeWYJg+llFKOafJQSinlmCYP\npZRSjmnyUEop5Vi47zBXKiREJBp4EOgOeIwxFVM81x14FPgFeN0Y81FoolQq83SHuVIBIiK5gQHY\nJHJfyiQhIg8bY54JWXBKZZGOPJQKnEbAV0A+bLe3lCOM06EISCl/0TkPpQKnqjFmI/AKUElEGgCI\niAC/iMgtItIzpBEqlUmaPJQKHA+AMeYwtkvig77Hb8Q2vloM9M/om4lIF38HqFRmafJQKgBEJAbY\nm+KhSUBLEakMFDbGJBljDgEZmnQUkTKc36JUqZDSOQ+lAqM5KdrqGmO2i8g8YBSwOvWLRWQNtnf7\nAWAw8CxQFvjV1zu8DhArIv2Ad40xyYH/FZS6ME0eSgVGOWPMrFSPTQC+AyaneMwjIq2APsaYnwFE\npIcx5l3f7R9E5EtjzPsicrcx5vWgRK9UOvS0lVJ+JCK1ROQVYJiIjE75nDFmJfAJ8H2Kh68AegOF\nUzyWciXWVqCC77bH9xn5/B23Uk7pyEMpPzLGrAEG+S5pPX9zqod+xe4F+UBElhtjTnP+32UZwPhu\ne0UkF1Ab+MavgSvlkCYPpULEt0z3WqAycAx4S0TuAkqIyC3A1cATvol1gB+Bbpy/X0SpkNAd5kq5\njIgsNcY0DXUcSl2Mznko5SIi0gMoLyJ1Qx2LUhejIw+llFKO6chDKaWUY5o8lFJKOabJQymllGOa\nPJRSSjmmyUMppZRjmjyUUko5pslDKaWUY5o8lFJKOabJQymllGP/D5sdVMZax/23AAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11224b128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the errors from pmax = 10^-8\n",
    "data = HCPdata\n",
    "Nk = np.array([N for (N,g0,gR1,gR2) in data])\n",
    "g0val = np.array([g0 for (N,g0,gR1,gR2) in data])\n",
    "gR1val = np.array([gR1 for (N,g0,gR1,gR2) in data])\n",
    "gR2val = np.array([gR2 for (N,g0,gR1,gR2) in data])\n",
    "\n",
    "gplot = []\n",
    "Nk53 = np.array([N**(5/3) for (N,g0,gR1,gR2) in data])\n",
    "for gdata, start in zip((g0val, gR1val, gR2val, g0val-gR1val, g0val-gR2val), (3,3,3,3,3)):\n",
    "    N10,N5 = np.average(Nk53[start:]*Nk53[start:]),np.average(Nk53[start:])\n",
    "    denom = N10-N5**2\n",
    "    g10 = np.average(gdata[start:]*Nk53[start:]*Nk53[start:])\n",
    "    g5 = np.average(gdata[start:]*Nk53[start:])\n",
    "    Ginf,alpha = (g10-g5*N5)/denom, (g10*N5-g5*N10)/denom\n",
    "    gplot.append(np.abs(gdata-Ginf))\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.plot(Nk, gplot[0], 'k', label='$G(\\mathbf{0})$ error $\\sim N_{\\mathrm{kpt}}^{-5/3}$')\n",
    "ax1.plot(Nk, gplot[1], 'b', label='$G(\\mathbf{R}_1)$ error $\\sim N_{\\mathrm{kpt}}^{-5/3}$')\n",
    "ax1.plot(Nk, gplot[2], 'r', label='$G(\\mathbf{R}_2)$ error $\\sim N_{\\mathrm{kpt}}^{-5/3}$')\n",
    "ax1.plot(Nk, gplot[3], 'b--', label='$G(\\mathbf{0})-G(\\mathbf{R}_1)$ error')\n",
    "ax1.plot(Nk, gplot[4], 'r--', label='$G(\\mathbf{0})-G(\\mathbf{R}_2)$ error')\n",
    "ax1.set_xlim((1e2,2e5))\n",
    "ax1.set_ylim((1e-11,1))\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_yscale('log')\n",
    "ax1.set_xlabel('$N_{\\mathrm{kpt}}$', fontsize='x-large')\n",
    "ax1.set_ylabel('integration error $G-G^\\infty$', fontsize='x-large')\n",
    "ax1.legend(bbox_to_anchor=(0.6,0.6,0.4,0.4), ncol=1, \n",
    "           shadow=True, frameon=True, fontsize='medium')\n",
    "ax2 = ax1.twiny()\n",
    "ax2.set_xscale('log')\n",
    "ax2.set_xlim(ax1.get_xlim())\n",
    "ax2.set_xticks([n for n in Nk])\n",
    "# ax2.set_xticklabels([\"${:.0f}$\".format((n*1.875)**(1/3)) for n in Nk])\n",
    "ax2.set_xticklabels(['6','10','16','20','26','30','36','40','46','50','56','60'])\n",
    "ax2.set_xlabel('k-point divisions (basal)', fontsize='x-large')\n",
    "ax2.grid(False)\n",
    "ax2.tick_params(axis='x', top='on', direction='in', length=6)\n",
    "plt.show()\n",
    "# plt.savefig('HCP-GFerror.pdf', transparent=True, format='pdf')"
   ]
  }
 ],
 "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.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}