1-GaussQuadrature.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3df22115-eb02-41bf-a2ee-18d72da91080",
   "metadata": {},
   "source": [
    "# Introduction to gauss quadrature rule for integration\n",
    "\n",
    "- [x] 1d gauss quadrature rule\n",
    "- [x] 2d gauss quadrature rule\n",
    "- [x] integration in 1d/2d\n",
    "\n",
    "Author  : Yang Bai @ M3 Group\n",
    "\n",
    "Date    : 2021.12.22\n",
    "\n",
    "QQ group: 628204857"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c14356e4-cb9e-45d9-9ebb-e5946f483cff",
   "metadata": {},
   "source": [
    "Question: how can we calculate the integration looks like: $\\int_{-1}^{1}cos(x)dx$ ?\n",
    "\n",
    "We know the analytical solution is: $\\mathrm{sol}=sin(x)|_{-1}^{1}$, how can we get the answer numerically?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9f41f221-7196-485a-b051-e182dfb0796a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1-th gauss point: xi=  -9.06180e-01, weight=   2.36927e-01\n",
      "2-th gauss point: xi=  -5.38469e-01, weight=   4.78629e-01\n",
      "3-th gauss point: xi=   0.00000e+00, weight=   5.68889e-01\n",
      "4-th gauss point: xi=   5.38469e-01, weight=   4.78629e-01\n",
      "5-th gauss point: xi=   9.06180e-01, weight=   2.36927e-01\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAESCAYAAAA8BeghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWp0lEQVR4nO3df4xlZ33f8feHDMnsr8F1GWxjvFg2NpQlMoghSWPHNtAWbJXalakSduVAKnDWFjRpBIX+YSGtHUAkqpCQ42hT0ThWYkCRXWxBgqIQh+6CIbNEDt0k3qa72IC861lohp2Vx66Xb/+413B7PbN755n7Y1d+v6SjmXnOc5/7PWfPvZ99zrk/UlVIktTiBZMuQJJ05jJEJEnNDBFJUjNDRJLUzBCRJDWbmnQB4/biF7+4LrzwwkmXIUlnlH379h2tqtn+9uddiFx44YXMz89PugxJOqMkeXSldk9nSZKaGSKSpGaGiCSpmSEiSWpmiEiSmhki0mng8ccf56qrruLw4cOTLkVaE0NEOg3cdttt7Nmzh127dk26FGlN8nz7KPi5ubnyfSI6XWzYsIHl5eXntE9PT/Pkk09OoCJpZUn2VdVcf7szEWmCDh48yPbt29m4cSMAGzduZMeOHRw6dGjClUmDMUSkCTrvvPOYmZlheXmZ6elplpeXmZmZ4dxzz510adJADBFpwo4cOcLOnTt56KGH2LlzpxfXdUbxmogk6ZS8JiJJGjpDRJLUzBCRJDUzRCRJzcYeIknOTnJfkuNJHk2yfZV+70pyIslSz3J1z/oHkyz3rHtkXNsgSeqYxDcb3gE8DZwDvBb4fJKHq2r/Cn2/WlVXnGSs91bVfx1BjZKkAYx1JpJkE3ADcGtVLVXVHuB+4MZx1iFJGo5xn866FDhRVQd62h4Gtq3S/3VJjiY5kOTWJP0zp4921+/tPdXVL8lNSeaTzC8sLKxrAyRJPzbuENkMLPa1LQJbVuj7ZeA1wEvozF7eAXygZ/0HgYuA84HdwANJLl7pTqtqd1XNVdXc7Ozs+rZAkvQj4w6RJWCmr20GONbfsaoOVtWhqvphVX0T2AW8vWf916rqWFU9VVV3AXuBa0dYuySpz7hD5AAwleSSnrbLgJUuqvcrIOtYL0kasrGGSFUdB+4FdiXZlORy4Drg7v6+Sa5Jck7391cBtwKf6/59VpK3JJlOMpVkB3Al8MVxbYskaTJvNrwF2AA8AdwD3FxV+5Ns7b7fY2u335uBv0lyHPgCnfD5SHfdC4HbgQXgKPA+4Pqq8r0ikjRGfoqvJOmU/BRfSdLQGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqNvYQSXJ2kvuSHE/yaJLtq/R7V5ITSZZ6lqvXOo4kaXSmJnCfdwBPA+cArwU+n+Thqtq/Qt+vVtUVQxhHkjQCY52JJNkE3ADcWlVLVbUHuB+4cRLjSJLWZ9ynsy4FTlTVgZ62h4Ftq/R/XZKjSQ4kuTXJszOnNY2T5KYk80nmFxYW1rsNkqSucYfIZmCxr20R2LJC3y8DrwFeQmfW8Q7gAw3jUFW7q2ququZmZ2cbS5ck9Rt3iCwBM31tM8Cx/o5VdbCqDlXVD6vqm8Au4O1rHUeSNDrjDpEDwFSSS3raLgMGuRheQIYwjiRpSMYaIlV1HLgX2JVkU5LLgeuAu/v7JrkmyTnd318F3Ap8bq3jSJJGZxJvNrwF2AA8AdwD3FxV+5Ns7b4XZGu335uBv0lyHPgCndD4yKnGGddGSJIgVTXpGsZqbm6u5ufnJ12GJJ1Rkuyrqrn+dj/2RJLUzBCRJDUzRCRJzQwRndTjjz/OVVddxeHDhyddijQUHtPDZYjopG677Tb27NnDrl27Jl2KNBQe08Plq7O0og0bNrC8vPyc9unpaZ588skJVCStj8f0+vjqLK3JwYMH2b59Oxs3bgRg48aN7Nixg0OHDk24MqmNx/RoGCJa0XnnncfMzAzLy8tMT0+zvLzMzMwM55577qRLk5p4TI+GIaJVHTlyhJ07d/LQQw+xc+dOL0TqjOcxPXxeE5EknZLXRCRJQ2eISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpqNPUSSnJ3kviTHkzyaZPsAt/lSkkoy1dP2YJLlJEvd5ZHRVi5J6jeJmcgdwNPAOcAO4M4k21brnGQHMLXK6vdW1ebu8srhlypJOpmxhkiSTcANwK1VtVRVe4D7gRtX6f8i4MPAfxpflZKkQY17JnIpcKKqDvS0PQysNhP5CHAnsNrnNX80ydEke5NcvdqdJrkpyXyS+YWFhYayJUkrGXeIbAYW+9oWgS39HZPMAZcDn1xlrA8CFwHnA7uBB5JcvFLHqtpdVXNVNTc7O9tauySpz7hDZAmY6WubAY71NiR5AfA7wK9V1TMrDVRVX6uqY1X1VFXdBewFrh1BzZKkVYw7RA4AU0ku6Wm7DNjf128GmAM+k+Qw8Ffd9u8k+YVVxi4gwyxWknRyq73qaSSq6niSe4FdSd4NvBa4Dvj5vq6LwEt7/r4A+DrwemAhyVnAzwJ/CTwD/CJwJfDrIyxfktRnrCHSdQvwKeAJ4HvAzVW1P8lW4G+BV1fVY/RcTE8y3f31SFU9033V1u3Aq4ATwN8D11eV7xWRpDEae4hU1feB61dof4zOhfeVbvMtek5VVdUC8IbRVChJGpQfeyJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGYDhUgSv3ZWkvQcg85E/iHJf0iy4pdGSZKenwYNkfuADwPfTvLxJBeMsCZJ0hlioBCpqvcDLwP+M3AN8L+TfDqJX1ErSc9jA19Yr6onq+p3q+qn6QTJBuArSf5Hkn87sgolSaetQS+svzTJtiRXJHkbcD7wJeDTwM8BfzzoHSY5O8l9SY4neTTJ9gFu86UklWRqPeNIkoZr6tRdAPgO8ENgCTgGLAI/6P68t/tzUHcATwPnAK8FPp/k4arav1LnJDtWqXNN40iShm/QEPlr4FLg94FPVNW3Wu4sySbgBuA1VbUE7ElyP3Aj8KEV+r+IzgX9Xwa+2jqOJGk0Br2w/nrg3wCvAB5J8pkkP9Nwf5cCJ6rqQE/bw8C2Vfp/BLgTOLyecZLclGQ+yfzCwkJD2ZKklazlwvpfVNW/pnPq6Bjwl0n2JLl+DW9G3MxzT30tAlv6OyaZAy4HPrmecbq1766quaqam52dHbBUSdKprPljT6rq74D3AJcB/wB8Fjhw0hv92BIw09c2QyeUfiTJC4DfAX6tqp5pHUeSNFoDXRNJso/O//630Hmy3vDsqu7Piwa8vwPAVJJLqup/ddsuA/ovhs8Ac8BnupOcn+i2fyfJvwO+MeA4kqQRGvTC+jeA/3OK5ZSq6niSe4FdSd5N59TYdcDP93VdBF7a8/cFwNeB1wMLVfX0gONIkkZooBCpqvcM8T5vAT4FPAF8D7i5qvYn2Qr8LfDqqnqMnovpSaa7vx7pOb214jhDrFOSdAqDzkSGpqq+D1y/QvtjdE6ZrXSbb/HjU2cnHUeSND5+n4gkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKajT1Ekpyd5L4kx5M8mmT7Kv1+KckjSRaTPJHkriQzPesfTLKcZKm7PDK+rZAkwWRmIncATwPnADuAO5NsW6HfXuDyqnoRcBEwBdze1+e9VbW5u7xylEVLkp5rrCGSZBNwA3BrVS1V1R7gfuDG/r5V9e2qOtrTdAJ4xXgqlSQNYtwzkUuBE1V1oKftYWClmQhJrkiyCByjEz6f6Ovy0SRHk+xNcvXwy5Ukncy4Q2QzsNjXtghsWalzVe3pns56GfBbwLd6Vn+Qzmmu84HdwANJLl5pnCQ3JZlPMr+wsLC+LZAk/ci4Q2QJmOlrm6Ez01hVVX0X+FPg0z1tX6uqY1X1VFXdRecayrWr3H53Vc1V1dzs7Oy6NkCS9GPjDpEDwFSSS3raLgP2D3DbKWDFmUZXAVlHbZKkNRpriFTVceBeYFeSTUkuB64D7u7vm2RHkq3peDnwm8Cfd9edleQtSaaTTCXZAVwJfHF8WyNJmsRLfG8BNgBPAPcAN1fV/m5gLCXZ2u33auArdE6B7QUeAd7TXfdCOi/3XQCOAu8Drq8q3ysiSWOUqpp0DWM1NzdX8/Pzky5Dks4oSfZV1Vx/ux97IklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSp2dhDJMnZSe5LcjzJo0m2r9Lvl5I8kmQxyRNJ7koys9ZxJEmjM4mZyB3A08A5wA7gziTbVui3F7i8ql4EXARMAbc3jCNJGpGxhkiSTcANwK1VtVRVe4D7gRv7+1bVt6vqaE/TCeAVax1HkjQ6456JXAqcqKoDPW0PAyvOIJJckWQROEYnND7ROM5NSeaTzC8sLKxzEyRJzxp3iGwGFvvaFoEtK3Wuqj3d01kvA34L+FbjOLuraq6q5mZnZxtLlyT1G3eILAEzfW0zdGYaq6qq7wJ/Cnx6PeNIkoZr3CFyAJhKcklP22XA/gFuOwVcPIRxJElDMtYQqarjwL3AriSbklwOXAfc3d83yY4kW9PxcuA3gT9f6ziSpNGZxEt8bwE2AE8A9wA3V9X+bmAsJdna7fdq4Ct0Tl3tBR4B3nOqcca0DZIkIFU16RrGam5urubn5yddhiSdUZLsq6q5/nY/9kSS1MwQkSQ1M0QkSc0MEUlSM0NEktTMEJEkNTNEJEnNDBFJUjNDRJLUzBCRJDUzRCRJzQwRSVIzQ0SS1MwQkSQ1M0QkSc0MEUlSM0NEktTMEJEkNTNEJEnNDBFJUrOxh0iSs5Pcl+R4kkeTbF+l3zuT7EvygyTfSfLxJFM96x9Mspxkqbs8Mr6tkCTBZGYidwBPA+cAO4A7k2xbod9G4NeBFwM/C7wZeH9fn/dW1ebu8srRlSxJWsnUqbsMT5JNwA3Aa6pqCdiT5H7gRuBDvX2r6s6eP7+b5A+BN46tWEnSKY17JnIpcKKqDvS0PQysNBPpdyWwv6/to0mOJtmb5OrVbpjkpiTzSeYXFhbWWrMkaRXjDpHNwGJf2yKw5WQ3SvIrwBzw2z3NHwQuAs4HdgMPJLl4pdtX1e6qmququdnZ2dbaJUl9xh0iS8BMX9sMcGy1GyS5HvgYcE1VHX22vaq+VlXHquqpqroL2AtcO/ySJUmrGXeIHACmklzS03YZzz1NBUCStwK/B7ytqr55irELyFCqlCQNZKwhUlXHgXuBXUk2JbkcuA64u79vkjcBfwjcUFVf71t3VpK3JJlOMpVkB51rJl8c/VZIkp41iZf43gJsAJ4A7gFurqr9SbZ23++xtdvvVuBFwBd63gvyJ911LwRuBxaAo8D7gOuramTvFXn88ce56qqrOHz48KjuQpJGYpTPX2MPkar6flVdX1WbqmprVf1Rt/2x7vs9Huv+/caqmup5H8jmqrqmu26hqt5QVVuq6qyq+rmq+rNR1n3bbbexZ88edu3aNcq7kaShG+XzV6pq6IOezubm5mp+fn7g/hs2bGB5efk57dPT0zz55JPDLE2ShmqYz19J9lXVXH+7n511CgcPHmT79u1s3LgRgI0bN7Jjxw4OHTo04cok6eTG8fxliJzCeeedx8zMDMvLy0xPT7O8vMzMzAznnnvupEuTpJMax/OXITKAI0eOsHPnTh566CF27tzpxXVJZ4xRP395TUSSdEpeE5EkDZ0hIklqZohIkpoZIpKkZoaIJKmZISJJava8e4lvkgXg0UnXcQovpvPBkqc76xwu6xwu6xyul1fVc77V73kXImeCJPMrvR77dGOdw2Wdw2Wd4+HpLElSM0NEktTMEDk97Z50AQOyzuGyzuGyzjHwmogkqZkzEUlSM0NEktTMEJEkNTNEJiDJ2UnuS3I8yaNJtp+k7+8mWepZnkpyrGf9g0mWe9Y/MqE635XkRF+tV7eMNeI635lkX5IfJPlOko8nmepZP9T9ucba/mOSw0kWk3wqyU+1jDPKOse9/9ZR58SOxzXWObHH99BUlcuYF+Ae4DPAZuAKYBHYNuBtfx/4VM/fDwLvnnSdwLuAPaPY5iHXeTPwC8BPAucD+4APjWp/Dlob8BbgCLAN+CfdOj42jv23xjrHuv/WUefEjsf1jD/Ox/fQtnXSBTzfFmAT8DRwaU/b3b1PGKe47THgqp62kRxka63zZA/a9WzzKPdnt+9vAA+MYn+upTbgj4CP9Pz9ZuDwqPffEI7Jke2/de7PiRyP6xl/nI/vYS6ezhq/S4ETVXWgp+1hOv8DPZUbgAXgy33tH01yNMne3in7OrXU+bpuHQeS3NpzmmM92zyKOntdCezvaxvW/lxLbdu663r7nZPkn65xnFHX2W+U+6/fWuucxPG4nvHH+fgeGkNk/DbTmdr2WgS2DHDbdwJ/UN3/onR9ELiIzqmF3cADSS6eQJ1fBl4DvITOg+EdwAcaxxplnT+S5FeAOeC3e5qHuT/XUlt/32d/37LGcVo0jT+G/beeOid1PK5n/HE+vofGEBmy7oWwWmXZAywBM303m6EzjT3ZuBcAVwF/0NteVV+rqmNV9VRV3QXsBa4dd51VdbCqDlXVD6vqm8Au4O3d1U3bPIo6e8a9HvgYcE1V/egTVFv35yrWUlt/32d/P7bGcVqsefwx7b/mOkd1PA67zmcN+/E9TobIkFXV1VWVVZYrgAPAVJJLem52Gc89JdDvl4GvVNXBU5UAZIJ1rlRH81ijqDPJW4HfA97WfYIZdDvWai217e+u6+13pKq+t8ZxRl3nOPdfv/Xsh6EcjwNqGX+oj++xmvRFmefjAnyazqs3NgGXM8ArN4BHgH/f13YWnVf1TANTwA7gOPDKcdcJXAOc0/39VcD/BD68nm0eUZ1vAr4HXLnCuqHvz0FrA94KHAZeTefVWV/i/3911sj23xrrHOv+W0edEzseW8afxON7aNs66QKejwtwNvDfuwfEY8D2nnVb6UyHt/a0/fNu3y1948wCf0VnmvyPwEPAv5xEnXTOix/p9j1I5/TBCwcZa8x1/gXwTLft2eVPRrU/V6ttlX/n3+juwx8A/w34qXHsv7XUOe79t446J3Y8Nvy7T+TxPazFD2CUJDXzmogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiDRhSd60ygdMfnfStUmnMnXqLpJG7O+AzwL/Cvhj4KvAt+l8lpZ0WjNEpMl7P51Pef1nVWVw6Izi6Sxp8t4FfNIA0ZnImYg0ef8I/GqSbwPfAI5W1fJkS5IG40xEmrxrgf8LfI7OtZAnk5w72ZKkwTgTkSYoyU/T+QKjnwRuB/4aeMJTWzpT+H0i0oQk2Qj8PZ3g+EVPYelM5ExEmpxrgQvofFudAaIzktdEpMm5oPvzmYlWIa2DISJNzoN0AuSBJDuT/Iskl0y4JmlNvCYiTVCSN9J5s+EbgLOBnwB2V9WvTrQwaUCGiHSaSBLgvwDvrKqzJ12PNAhPZ0mnj4uBnwH+bNKFSIMyRKTTx2eBQ8DOSRciDcrTWZKkZs5EJEnNDBFJUjNDRJLUzBCRJDUzRCRJzQwRSVIzQ0SS1Oz/AZWR1r7KJu9yAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from FEToy.fe.gaussrule import gausspoint1d\n",
    "my1dgp=gausspoint1d(ngp=5) # up to 5\n",
    "my1dgp.creategausspoint() # generate gauss points\n",
    "my1dgp.print() # print out basic message\n",
    "my1dgp.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "69fe0670-0d6e-4173-a6b7-41124e7f5ce9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error=4.70247e-08 %\n",
      "1-th gauss point: xi=  -9.06180e-01, weight=   2.36927e-01\n",
      "2-th gauss point: xi=  -5.38469e-01, weight=   4.78629e-01\n",
      "3-th gauss point: xi=   0.00000e+00, weight=   5.68889e-01\n",
      "4-th gauss point: xi=   5.38469e-01, weight=   4.78629e-01\n",
      "5-th gauss point: xi=   9.06180e-01, weight=   2.36927e-01\n"
     ]
    }
   ],
   "source": [
    "analyticalsol=np.sin(1)-np.sin(-1) # this is our analytical solution\n",
    "\n",
    "ngp=5\n",
    "my1dgp.setnum(ngp) # set up the number of gauss points\n",
    "my1dgp.update()    # you must call update after the set function\n",
    "sum=0.0\n",
    "for i in range(ngp):\n",
    "    x,w=my1dgp.gpcoords[i,1],my1dgp.gpcoords[i,0] # 0->weight, 1-> coordinate\n",
    "    sum+=np.cos(x)*w\n",
    "\n",
    "error=100.0*np.abs(sum-analyticalsol)/analyticalsol # in percentage %\n",
    "\n",
    "print('error=%g %%'%(error))\n",
    "my1dgp.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58a31249-3038-4c0d-891a-bdc9e42477e8",
   "metadata": {},
   "source": [
    "Now we calculate the 2D integration, i.e., $\\int_{-1}^{1}\\int_{-1}^{1}x^{2}y^{4}dxdy=\\int_{-1}^{1}\\frac{x^{3}}{3}y^{4}|_{-1}^{1}dy=\\int_{-1}^{1}\\frac{2}{3}y^{4}dy=\\frac{2}{3}\\frac{y^{5}}{5}|_{-1}^{1}=\\frac{4}{15}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d336fc1a-484d-47e8-b68b-5cdb2dd3d685",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1-th gauss point: xi=  -7.74597e-01, eta=  -7.74597e-01, weight=   3.08642e-01\n",
      "2-th gauss point: xi=  -7.74597e-01, eta=   0.00000e+00, weight=   4.93827e-01\n",
      "3-th gauss point: xi=  -7.74597e-01, eta=   7.74597e-01, weight=   3.08642e-01\n",
      "4-th gauss point: xi=   0.00000e+00, eta=  -7.74597e-01, weight=   4.93827e-01\n",
      "5-th gauss point: xi=   0.00000e+00, eta=   0.00000e+00, weight=   7.90123e-01\n",
      "6-th gauss point: xi=   0.00000e+00, eta=   7.74597e-01, weight=   4.93827e-01\n",
      "7-th gauss point: xi=   7.74597e-01, eta=  -7.74597e-01, weight=   3.08642e-01\n",
      "8-th gauss point: xi=   7.74597e-01, eta=   0.00000e+00, weight=   4.93827e-01\n",
      "9-th gauss point: xi=   7.74597e-01, eta=   7.74597e-01, weight=   3.08642e-01\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAADzCAYAAAB+IVisAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABx7klEQVR4nO19d3hc5Zn9uVM06ppRt3qzbHVblkzLAqHjgG0CwZBdyC9ACKGEOJUUEqcBWRKybMgSNgkLIZTFJlQbAwmwAeK4q1qSJVm9zYxmRtP79/tD/j7fGU25M3NVbM95Hj2PPbr3zh3Nvee+3/ue97wcIQRxxBFHHKEgWe4TiCOOOFY+4kQRRxxxhEWcKOKII46wiBNFHHHEERZxoogjjjjCIk4UccQRR1jIwvw+XjuNI47FB7fcJxAO8YgijjjiCIs4UcQRRxxhESeKOOKIIyziRBFHHHGERZwo4ogjjrCIE0UcccQRFnGiiCOOOMIiThRxxBFHWMSJIo444giLOFHEEUccYREnijjiiCMs4kQRRxxxhEWcKOKII46wiBNFHHHEERZxoogjjjjCIk4UywBCCJxOJ9xuN+LjEuI4HRDOuCYOkeH1euF0OmG329lrUqkUcrkcMpkMUqkUHLfifUziOMvAhXmixR93IoEQArfbDbfbDY7j4HK52OuEEHi9XkYQceI467Div+A4USwB6FKDTwaUKAJtS7cfHh7G6tWr48Rx5mPFf6Hxpcciw263w+l0Qi6Xg+M4cBwHr9cLvV6PlJQUJCQk+GxPt5FKpbDZbJBIJPB6vbDZbIwgZDIZ+4kTRxxLgThRLBLoUmNqago2mw2VlZUAAKfTic7OTnAcB6fTCY/Hg4yMDKhUKiiVSsjlcgDzhEEIYcQhkUjYcT0eD9xuN3svmUzGIg6JRBInjjhER5woFgFerxcul8tnqQEAOp0OPT09qKqqglKpBMdx8Hg8mJubg16vx+joKAghUCqVSE9PD1gRocRB4U8cHMf5RBxx4ohDDMRzFCKC3rQ0/yCRSDAzMwOj0QgA0Ov1aGxshEKhgNPpBIAFN7Hb7YbBYIBOp8PU1BRSU1OhVCqhUqmQkZEBqVQa9hzoz+joKAoLC5GYmBgnjpWNFf+FxCMKkUAIgcvlgsfj8Xnqu1wuTE5OorCwEK2trex1urTwh0wmQ3Z2NlQqFUwmExobG2EwGKDVajE4OAipVOpDHHRJQsF/b4PBgMLCQrjdbkZe8YgjjmgQJwoRQLUR/JwCAKjVagwODkKpVKKqqiqqY8vlcuTk5CAnJwfAfI7DYDBArVZjYGAAMpkMKpUKKpUKaWlpAYmD/xolNH/ikMvlkEqlceKIIyDiRBED/LUR9Ib0er3o6+uD1WpFdXU1W3pEgmARR0JCAnJzc5GbmwsAcDgc0Ov1mJychMlkQkJCAiOOYDkO/vIlEHHQxKhMJluQE4nj7EScKKKEvzaC3kwWiwWdnZ3Iz8/H2rVrodPpFlWmrVAokJ+fj/z8fADz5Vi9Xo/x8XGYTCZ0d3cjKysLSqUSqampC276QMThdDrhcDjgcrlgMpmwatUqFnHEiePsRJwoooDH48HMzAwMBgMqKirYjTM5OYnh4WHU1dUhIyMDQPDIIByi3S8xMRGrVq3CqlWrcOTIEVRUVMBsNmN0dBRmsxnJycks4khOTg5JHA6HA2q1GllZWXA4HCxq8l+qxHHmI04UEYC/1KD/5jgObrcbPT098Hq92LhxI2SyU3/WaG94McBxHJKSkpCeno6CggIQQmC1WqHX63HixAlYrVakpKQw4khKSgpJHPRzOJ1OVrWJE8fZgThRCIS/NoIqJo1GI7q6ulBSUoLCwsKAN1q0EYXY4DgOKSkpSElJQVFREQghsFgs0Ov1GBgYgM1mQ1paGhN/0eSs/znFiePsQ5wowsBfG8FfoxuNRhgMBjQ2NiI1NTXg/ssZUYQDx3FITU1FamoqiouLQQiByWSCwWDA8ePHYbPZQAjB9PQ0VCoVFArFgv0BX+KgOQ5KHB6PB1KpFCkpKXHiOI0RJ4oQCKWNGBgYgMvlwnnnnRdSBLWSicIfHMchPT0d6enpKCkpwdzcHIaHh2G323Hs2DG4XC4fuXmwPhUKQgg0Gg1sNhtKSkoAzEcc/Aa3OHGcHogTRRAE00bo9XocO3YMBQUFMJvNYZWSpxNR+EMikSAhIQFlZWUoKyuD1+vF3NwcDAYDxsfHg/apUPAjDqlUyiIOh8MBh8MBQojPMoWWY+NYeYgThR/4Sw3/ZqwTJ05Aq9WiubkZbrcbJpMp7PFOZ6Lwz1FIJBKW+CwvLw/Yp8InDplM5nOMQBGH1+uNm/icBogTBQ80VJZIJD6aA4fDgY6ODmRkZKC1tRUSiQQWi0UQAQQjitOZQCikUikyMzORmZkJYL5PhRLH8PAwq5gkJSWxXAUfceI4fRAnipOgSw2tVoukpCSkpaUBADQaDY4fP441a9YgOzubbU99JcIhGCF4PB72+5UK/4giHGQyGbKyspCVlQVgPpdDy7BHjhzx6VNJT08XTBxTU1Mwm80oLi6OE8cy4awnCn8ZNl1Le71e9Pf3w2QyoaWlJWDGP9qIYmpqCv39/azqkJmZCZVKhcTERFE/23JDLpcjJSUFaWlpKCgoCNinwieOUA1uhJCAJj5x4lganNVE4a+NoDkJu92OgwcPIjc3Fxs2bAh48UkkkoiJwuv1ore3Fw6Hg3WSUh1Db28vnE4nW+OrVCrRP2+kiDSiCHeMYH0q09PTOH78uE+fSmpqqk/vDP1+/PNGcfevpcFZSRShtBEmkwkajQbr1q2DUqkMeoxIlx5WqxUdHR1YtWoVampq4Ha74fV6kZaWhrS0NJSUlDABl06nw9jYGCwWCwYGBlhyMFyFRWyITRT+CNenkpSUBKVSucAAiCKY+xdNRANx9y+xcNYRhf9Sg144Ho8Hvb29MBqNKC4uDkkSQGQRhd1ux9GjR1FfX896QIIdU6lUsvc+cOAAlEoldDodhoaGfKoOgUL1lQiv1yv4PPl9KoQQRhyTk5Ow2+0wGo1h+1SC2QZS4VhxcXGcOKLAWUUUgZYawHwU0dXVhaKiIqhUKp+sezAIyVHQPIfD4cCFF164QGcg5D2ys7NZEtXpdPqE6gqFAiqVCpmZmUhJSRH9ol/siCIUaJ9KUlIS01+oVCro9XoMDQ3BYrEI6lPh5zgmJydRUFAQtw2MAmcFUYTSRoyPj2NsbAwNDQ1IS0vDzMyM4CVFqO1sNhs6OjqQnZ2NlJSUiEkiEBISEpCXl4e8vDz2Hnq9HiMjIzCbzezGyczMRFJSUszvJwbEIhsqAw/Xp0KTo/6f3+v1LlCC0ugy7v4VHmc8UVAZ9tGjR9HU1OQjw+7u7oZMJsM555zD1v/RJCn9oVar0d/fj9raWmRkZECj0Yj3gXigT1zaGUpvnP7+fjYmQK1WIycnZ4HcWgiWM6Lgw+v1+nTkAoH7VMxmM/R6PY4fPw6Hw8Ea3FQqFSQSSciqCj3XuPtXYJzRRMGXYVssFvYlGwwGdHd3o6KiAqtWrfLZJ5IkZaD36+/vh9lsRmtrKxISEuD1epdEWOV/43i9Xhw+fBh2ux1dXV3weDzsaUtVk0sBsYgiXJ6D47gFiWGTycQk93TW68zMDFQqVUDi5LfU03P3Jw5+OfZsIo4zkigCWdTRCGB4eBgzMzNYv349kpOTF+xLa/WRwm63o729HTk5OWhubvaRLS+HApO2e5eUlEAulzO5tU6nY6pJ+rQNZNILrJyIgmooIoFEIkFGRgYyMjJQVlYGi8WC3t5eWK1WTExMhO1TAQITx/DwMGQyGfLy8nwa3M5028AzjiiCWdR5vV4cOXIEKSkp2LhxY9ALT+jSgw+q3qypqWFyZopgF060hBQt/OXWLpcLer2eiZ/kcrmPSa9YF7xYEYUY55OYmIjy8nLBfSr+oKRPyYFvG0ivNblcfkZ6cZxRROGfmKIX1+zsLDO6pWKfYBC69ADmb4Ljx4/DaDQGVG+G23c5IZfLfcRPfA2D2WxGUlIS5HI5I85ob9SlWnpEeoxwfSoAGGnwNSz8nhV+xOFv4nPvvffiBz/4AdauXRvTea8UnBFEEUwb4fV6MTg4CIPBgNTUVNaDEApCn/R2ux1WqxVSqTSoevN0gr+GwWq1YmxsDDqdDgcPHkRqaiqrqERKiCth+UKrHsEQqE+FLtVOnDjBNCxWqzWgxoYvKQfmp8KtlMqTGDjtiSKYNsJms6GzsxNZWVloaWnBwYMHw14sgDCimJ2dRW9vLxQKBZspeiaB4+Yt8zIzM5GYmIjS0lKWGOzp6YHL5UJ6ejoyMzODru8pVkpE4fF4IjqGXC5foGGhfSr9/f0YHR0N2adisViCup4FA8dxVwF4HIAUwB8IIY/4/T4DwJ8BlGD+3v0lIeR/InqTKHHaEgW/Bs7XRgDAzMwMBgYGUFtby3omhEYKoZKPhBAMDg5Cr9ejpaUFhw8fjvkznA6RCMedcr4qLS1lBjY6nY6t7/mJUf8E4HIkM/0RK9nQPhWNRoOysjLIZDIYDIaAfSopKSmwWq0REQXHcVIAvwVwOYBxAAc5jnuDEHKMt9k9AI4RQq7lOC4HQB/Hcc8TQpxRfzCBOC2Jgpat2tvbUVpaymTRHo8HfX19sNvtrDxJIZQogm1HPSmUSiVaWlpOixs8VgS7yflScuDUvFQ69pA/vUyMRKRYxxAjuUgjE4VC4SN+4+d4tm/fDrVajd/+9re45JJLfPQ7IbARwAAh5AQAcBz3EoAtAPhEQQCkcfMHSwWgA+D2P9Bi4LQjCr42graEA4DZbEZnZycKCgpQU1Oz4IuJhSjoFHJ/T4pY4PF4RLt4FwtCowE6L5Ufput0OkxOTkKn08HlcrF5qtFIzRcjmRnLcQItX/k5nnfeeQfnn38+MjIy8Nvf/hZPPfWUkM9cCGCM9/9xAOf4bfMEgDcATAJIA7CNELIkpbPThigCaSOkUik8Hg8mJiYwMjKC+vp6pKenB9w/EqKg5EPt72ZnZ7FhwwbR/CL0ej26u7sBzHdQUj+KxejXWA4kJCSwrlCXy4Xi4mJYrVYMDw+ztTu/RyMcVsLSgyKQU5c/qBfobbfdhttvv13ooQN98f5r4CsBtAG4BEAlgPc4jvuIEBL5zMoIcVoQRTBtBAAMDg4iMTFxweAdf0SSo6BRS2dnJ1JTU9HS0hL0IotkDU4IwcjICGZmZrBu3TrIZDI4HA4mgrJYLEx27K/HWA6IlUOh08kKCwsDSq35idFAikkxlh5CbnAhEEo4UfztxgEU8/5fhPnIgY8vAniEzD/JBjiOGwKwFsCBSN4oGqx4oqDNXP5u2EajEVNTU8jLy0NtbW3Y40QSUTidThw8eBDV1dVsiniwbYVeEG63G11dXZDL5WhtbYXX64Xb7UZSUhIKCwvZTWQymaDT6XDs2DFYLBb09/ezm2ip/SjEgP/fJ5DU2mg0svW91+tFRkYGMjMzkZGRwYRNYkQUYjTmCfm+3W53NN/VQQCrOY4rBzAB4CYAn/fbZhTApQA+4jguD8AaACcifaNosGKJItBSg74+MjKCqakpFBQUMG/LcBBCFFSi63A48KlPfSpsWEyjj3AXscfjwYEDB1BaWorCwkL2XoGOR6sLZWVlOHDgADIzM1ktny8SElM9GQxLoarke3BQxaTBYGDt5BKJBDabDUajEUqlMmrCEFIaFwu0BT4SEELcHMfdC+AdzJdHnyaEdHMcd9fJ3/8OwE8BPMNxXCfmlyrfIYRoxT37wFiRRBFMG+F0OtHV1YWkpCScc845bLaEEEil0pBEQY+dnJyM5ORkQWtnIX0c09PTsNlsOO+88wSTGv/4mZmZTAREk4TUASo5OZkRx0oV90RKNlKp1Ef45HQ6cejQISY1px4c1C5P6LGXMnEcjYYCAAghewDs8Xvtd7x/TwK4IuYTjAIriihCWdTRysPq1auZ7FgikQgmilDb0m7Sqqoq5OXl4R//+IfgYwYjCq/Xy0q11GA2VvCThFQ9qdPpFqz1VSrVkoXZi32MhIQEyOVyJoW22+1Mv8H34AhmXkOx1EQRaUSx0rFiiIJqIzo7O1FUVMRkslTkpNPp0Nzc7PPklEgkjFTCIdDSgy5jpqeng3aThkKwvhC73c5Ma9auXYt9+/ZFdcPQiCXQflQ9mZKSwtrK+X6bNAFsMBiQlZW1bGVYsUVliYmJKCgoWODBMTAwALvd7uNBwZeai0EUQj+L2WyOKqJYyVgRRMHXRvBvaLvdjs7OTiZy8v+iaXlUCPyJwuVyoaurCwqFImQ3abhj+kcUNPJZu3YtC5+XouTp77fpcrlw5MgRaLVaDA0NRWWbJ0YScTHVp4E8OPgeFG63m3WEulwuUWTgQvIcNNI5k7CsRBHIok4mk8Hj8bDWbf4N549IWrX5287NzTHjGuoAHQ34EQVNhKrV6gWai1CRQbjjR9tlKpfLkZCQgKqqKsjlcmabx9cy0PxGJE1e0WCptCH+HhT8VnKdTgebzcY+c6ABROEglCgilW+fDlg2ogimjZBIJBgdHQWABTJsf0QaUTidToyOjmJiYgJNTU0xsz69kd1uNzo7O5GYmMhGDgbabqlgc3kwqLFi1OhG7cn39bfNM5vNrAzrdrsDul+dLr0owcCvElksFlRUVMBms0Gj0bABRDTKElJFErp8iecoREKwSeHUfUipVKKxsTHsFxdJREFdmJVKJTZu3ChKqYzjOJjNZnR0dKC8vHyBrR5/u2iIIpr9po123PtyN4w2N+wOJz7RDeChLTWQSQJrGUpLSxeUJKVSKVQqFdxu9xnzZPR6vUhISEBqairTxlCxG7+KFGocgNCIItqqx0rGkhJFMG0EMD9mb2hoCPn5+QG/pEAQGlEYjUYMDg4iPT0ddXV1gs813DnY7Xb09/dj/fr1IS+MpYwoHv3rCcyanUhOkMIjITg4Moc9XWpsbswLuk+gkiRNiqrVami1WlZNCVVZWMkIFA0oFAofDw6bzcY0K1ar1ScxmpiYGBFRrARlrZhYMqIIpo1wu93o7e2F2+3Gxo0bMTMzE1Mlgw9qxz8+Ps7CzkiOG+yioKMBbTYb6uvrwz49ljKiGNHZIJee8uv0EIIRnTWiY9AyrMPhgEKhQFpaGnQ6HassiF2GXQqEWzZwHMc0NHQcAE2M0nGPNJfjcrlCfm6LxYLi4uKgvz8dsehEEUobYTKZ0NnZiZKSEhQWFjJrMSEDeIDQEYXb7caxY8fAcRw2btwIg8EAi8Ui6Lih9BF0XkdeXh5ycnIErVkjsdeLFWtzU/DxCT1kEnKyw5ZDdV70YbCQMixd5wcz6Y0FYkVikeZbAnlwjI6OYnZ2Fh0dHSCE+OR1+A+VRTKt+RaAfz35XxmAGgA5hBBdRG8UJRaVKOhSo729HfX19eyLIoRgbGwMExMTaGxs9PmjxlLypKAExJdMR5LPCHZjU2craojT29sbtbW/kIs2moji65dWYGKuB2N6G5xeDp+py8bla6NrjQ90c/mXYd1uNzPp7e/vR0JCAksgipHQE1MoFcuSSSKRICkpCdnZ2SgtLWUeHHx5PfXfiFRHIcS0hhDyKIBHT25/LYDtS0USwCISBX+pYTKZ2JdE9QsJCQkBk4qxEgVtOaeTv0JtK/S4/HZzvomu0Bs5EPFYLBYYjUZkZmaK2oOQmZKAP/xrI2ZMDhw/1oVPtZZDsog5BZlMhpycHJYgpMpJWoZ1OByYmpqKugy7kiov/BxFIA8OvV6P3//+93jrrbfQ3d2N66+/HrfeequQfIUQ0xo+bgbwoggfSTBEJ4pA2gh6oxiNRnR3d6OysjKofiFcT4b/tpRUPB4Pjh07BkJIwJbzaImCqkVTUlIWiL4imSrGB7Xqy8jIYHMi6FOY378QbW5DKuFQkJGIacXSi6X8lZP79++H0+n0EUDRblghQ4hWkrlPqHOh4x5/8IMfYGBgAHfeeSdGR0fZnNMwEGJaAwDgOC4ZwFUA7o3s7GODqERBZdgej2eBNoLKsMNJpSPp36DHp+5WRUVFKCoqCmrfFunSw2g0oqurC5WVlczyLNB2Qo9HCEF/fz9MJhNaWlqY8pGW6Wj/AnUMX8q5H4sBmnMqLS1lZVj+ECJqqUd1DMGGEK0UovB4PIKStxaLBbW1tbjyyiuFHlqIaQ3FtQA+WcplByAiURBC4HA4FmgjHA4HzGYz0tLSAoqR/BHJ0gOYf+J3dHQsWGoEOm4kEcX09DQ0Gk1IYVYkSw8qqU5PT0dzczMTnAELy3Rmsxmzs7Mwm81ob29n0UYsbdbRQOyw33+WBg3XJycnYTQakZSU5NMNSwl2JS49QoGWViOAENMaipuwxMsOQESioOTA/1K1Wi36+vqQmpqK0tJSQRe5UKLweDzMOv6CCy4IG8YKjVToU8/pdKK1tVUU1yxqWlNdXc0ik2AEwxdDGY1GVFZWwmazQavVsjZrejMJ1ZvEgsU8Pn86O+2G9W/wEiMhKlblZBGVmUJMa8DN2/VfBODfIjm4GBB16UFvHDqsl07Q6u/vj8g3Ity2FosFHR0dKCwsxNzcnKC1rpCb2mq1oqOjAwkJCaioqAh7XCERxdTUFPR6Perq6gIuX8IdXyKR+CTN+KIgm822qJqGWG+wSPbnl2GLiopYEnx6ehoGgwGHDh2CUqlkZdhIEsBL6ZdJt4tkCLRA0xoAuA7Au4QQYXV+ESF6MpPqDHJycpitvVQqFZrUCUsU09PTOHHiBOrq6pCRkYHx8XFBxw1HFLQJra6uDlNTU4KOGYooKFlarVbk5uaKZizDt86jN9Ps7CzGxuZzYSqVCllZWaL4X8S69Ihlf9rgRf/GVVVVPiMB5HI5I8hwBjaRDv8JdRwhxrrRIJxpzcn/PwPgmajeIEaIShTT09Po7+/3GbwDRJZ3CGU80tvbC4fDgdbWVp+np5ALMthNTQjBwMAADAYDa0KbmZmJaQaI0+lkuYV169ahp6cnqsRkuIiF3y0JnBo8PDk5CZPJxEqT2dnZy+KAJZbxDe0q5kdW/gY2/JGH/m7pYtngCSUc/yX4mQBRiUIqlWLjxo0LQuBIE5T+oEuC/Pz8BTM7wsmtKQJ9cU6nEx0dHcjIyPAZ6hOpYzcfc3NzLB9BtQXBLhqxLyb+4GFCCA4dOgSv18scsGjorlKpBN04yxlRUARLZvqXYWk3bG9vL1wul08ZdrFnevhjKTuFlwqiEkVubm5AQoiFKKjmoK6uLuBwWHrsSJ8Y9IbmW+tRCNVHSCQSnyXVxMQERkdHsW7dOp9kVrQS7liayeiSr7CwcEFpcmhoKKh2Q0yIRRThbvJA3bD8Miwh83L2ubm5oGVYIRBynTmdztOm/yUSLElTmEwmE5yjoKCekzabLaQvRST6CMBXPh5M0xGpPoKeK10W+SeyltqPIhD8S5OBtBvUyJf+rVdCRBHNMfw/6+zsLEZHR9mSLDExcUEZVgiELD3ORC8KQGSiCPYHl0qlTDMgBF6vFwcOHEBeXh7Wrl0b8ouM1JOis7OTNYoFezpEsvRwuVw4dOgQcnJygp5rLN2ji4Vg2o3u7m54PB4olUrY7faYRF9LFVGEg1QqRUpKCqqrq33ayQcGBhZUjkIZJQEQRBRnmhcFsEQRRSRLD7VaDavVipaWFkE9/UKPbbFYYLVaUVpaiqKiopDbCl162Gw21tgWaiZpqERquBtpKSIRfuheVlbGGp5mZ2fR19fnMxYgEu3GSiEKfiTg305OK0c6nQ4TExPwer1Rl2GBeEQRE4TczLScaDabkZGRIZiVhTz9aWdjYmJiWJIAhC09qM9FXl5e2MHFS5XMFAu0wjA7O4tVq1ZBLpf7aDfS0tKQlZUVVruxEuz+gdBJSH7lqLy8POBkdkqSQkg7ThQCEOwLDZejoPb2WVlZaG5uxpEjR6JqDPOH1+vFwMAATCYTWltbcejQIUEXXqgxAF6vFz09PfB4PKiurobBYAh7joGIZ7HazMUEfe9ItBv8p/9KiSgiOUagMqxer8fo6CisViu6uroYcQQaWh1fesSAUDczFTrV1NSwpYYYnhQOhwMdHR1QqVRobm5mKkchJa5gSw+73Y729nbk5+ejpKQEOp1OcC5juZOZ0SCYH0Uo7Qa/X0OsiCJWDUQsZJOYmIhVq1YhPz8fBw8eRFlZGSvDOp1OZl6jUqkgk8misuoPZ1pzcpuLAfwHADkALSHkoqg+UJRYNqLgC534Hg/Bto/k2HSug/+QYdoYJkRz4U8A9Jh8QoukzZy/HSEEQ0NDzOw3Kysr4BDi04Fg/LUb/Oll1HpQq9UKbiv3h9frjWo//2OIEZVIpVI2R6SkpMRnHMDIyAi6u7vxwQcfIDs7O6xdHsXJazekaQ3HcUoA/wXgKkLIKMdxuYGOFQgcx/0TwCsnjW/AcdwfAdwGYBUhZJrjuHTM95dcQwj5v2DHWRaioIN9VCqVj9Ap2PahwI8oCCEYHR3F1NTUgqli/ttGcsyxsTFMTU0FndcRDvztPB4Pm2q+fv16FsafOHECcrkcWVlZS9bwFQ7R2MfxbfMMBgNGRkZgMBhYWzktwQrVboh1k8dKNoFKo/5l2LKyMnR1daGjowMtLS148MEHccMNN4Q87oEDB4DwpjWfB/AXQsgoABBC1BGcugFA2slj5wK4AYAOgArANIAvnHz/oCQBLEOOgtrJhRrsEylReDweuN1udHd3QyaTBZ38JbSDlBIFNcPhOA4tLS0LnviRRhQ2mw3t7e0oLCxEUVERnE6nz4VGy3aDg4Ow2+0ghCAxMRHJyckxX+iLDZfHizG9HW4vQZEyEckJUnAch8TERFRVVQEIrt0I5X4lxvLF4/HEPORIiNgqOzsbZWVl2LhxI774xS8KutYmJiaA8KY11QDkHMd9iPmb/nFCyJ8EnroeAE2a3APgVQDrANCS4t0AFix1/LEkVx+98QYGBqDX6xc8mf0RqcuVzWbDgQMHfDwyQ52HkPN1Op04ePAgCgoKUFxcHFQfITRHYbFYMDIywvpgAhGMf9Kwu7sbFosFR48eZdn3rKwswSMBY4XQm9Tp9uJPB8YxPGuDBEBaogy3nV8Mqd/+gbQb/CFEfJNeelMudTIz1DGEWvXT70cIuQd50Pi/KAOwAcClAJIA7OM47p+EkONh3+BkRMFxXCKAuzA/Df0/Aag4jrsM85HFS+EOsiRE4XA4YLVaQQjBhg0bRDWvMZlMmJmZwYYNG5Cenh5yW6FEYTQaoVar0dzc7NPc5g+hSw86XGfjxo2Cm7MkEgkUCgVyc3OhVCrZ05h6UaanpwsqUS4F2ieMGNJaUaSa/2wakwPv9mhxVWVw1WMg2bVer2dlSWrS63Q6V0TlZLGG/5ws14czrRnHfALTAsDCcdzfATQBEEIUegDlmPew6CKEtHMcZ8Q8QVwP4ElCiCPcQRZ96UGH9ioUCqxevVrQcYQsEWiz09zcHAoLC8OSBD1uuDkgIyMjmJycRHZ2dkiSEHI82vFqsVhQUlISUwcn/2nML1GOjo6C4ziW2xAyGk8ohEYUBpsLctmpGzFFIYPO6gQhiYLPRSqVBvTdMBgMMBgMLOkbDTGKLdoKhUjnjra2tgLhTWteB/AEx3EyAAmYX5r8WuBbGDC/XPkagO+cfG0O88uPywHcKeQgixZR8J2rN2zYgCNHjgjeN1xEQXUX2dnZqKiogNFoFHTcUDc2P8lYV1eHkZGRsMcLFVHQVvPs7Gykp6dH3OsS6vj8EmVFRQWb7EVH46WlpSEzM3PJPDdLM5PxwfFZuDxeSCUcdBYnLq7Oiim/QJdhJpOJjWr0127QYcPhbuDlWHoIxcnlSUjTGkJID8dxewF0APBivoTaJfAt9Jh3xZrAKb8LI+YJ4mWhiVHRiYLjODgcDnR2diI1NdXHuVrohROqN4RGKDQZOjs7G7Pmwmq1or29HcXFxSgqKoLZbI6q7ElhNBrR2dnJyrOTk5OLWuakk73y8/PZhCudTgeLxYK2tjaW20hPT4/oxhX6fVXnpuDa+ly826uFx0vQWqrExauzMKfXieZHkZaWtkC7MT09jePHjy/w2vTHSl56AIJNa9hcjwhBk5n/QU5dhHP0NaEHEZ0o9Ho9uru7F7Rv0yhBSIInmO5ieHgYarXaJxkay7wO4JSvZ319PbsQoymjUlAHrqamJnbBxNIUFul+HHdqwpVOp0NNTQ1MJhMmJyfR29uLlJQUtkyJtRLAx3kVmTi3XAUvmR8XACyeMjOUdoPvu0G1G0tJFNRIeqWAELILfi7fhJDvAvhuJMcR3a5/eno6oIZBJpNFTRRutxudnZ1ITExc4OQdi+ZiaGhowVAf/+1CgX8jUwGZ0Whc4MC1nMIp/5vKYrFgdnaWVRrokzjQOMBodBRS3uZiOlyFek//kYd0ghfVbjidTthstpgGLAslmygcuE8LiJ7MrK2tDXiTRXtD0/GA5eXlbK0abFuhx+UTT6AqTKT6CLfbjY6ODqSkpDC5eKDtIoXYBMNxHFMW0rF4/HGANITPyspiEdty+1FEatdPRV18342jR49iamoKAwMDgrQbgeDxeMK2oAPzRBFqbs3piiVT8URjsDs5OYnh4eGQMzsinddBNRdlZWUoKCgIuJ1QfQRtHhNyvJUoxeaPA+SH8NROzuVywWAwIDs7O6p+i5XQFKZQKCCXy1FTUwOJRCJIuxEIQpceYvlzrjQsKVFEYrCr1+vh8XgCjgfkI5LJYhaLhekjQpVThUYpOp0OVqsV55xzDstvBMJKiSjCvRc/hPd4PDhy5Ah0Oh1GRkagUCh85OVCsBRLDyGgZBNKuzEwMICEhISgEnoh5dGV+DAQC4tS9QgEoURht9vR1dUFjuPQ1NQU9kITclxCCBtpWFxcvIAkDFYX3uiYgtrowNr8NFxRmxPyS6d6i5mZGSQnJ4ckCWDlRhShIJVKIZfLUVlZiYSEBNhsNszOzrLhPEKMepdj6REMgY4RTLvB992gyxShkQLHnXkO3MASRhQ0mRkKtA9k9erVGB4eFvQHD/f05w8ZplO3+LC7PPjVX/uhNTuRlCBF34wZGrMDa4KMfuT3f7S2tuKf//xn2HM8HSKKQODf6ElJSWy2Kz9hSI16s7KykJWV5fMkFqPaIEZEIRSBfDeoPsVqtbLrLJh2g5r4nolYETkKvjirpaUFcrkcg4ODgo4b6mYym83o6OhARUUF8vPzoVarF5DV8KwVGpMT+RnzybtUhQz/PKFHRWlwP4pVq1YF7f8Qeo42mw0nTpxgtvLLLcOOBP4JQzpjgz+9LCsrC263O+YSrBhRSTTwd76iyepQ2g2n0xnx5927dy+uvvrqPgTxojjpQ/E6gKGTL/2FEPKTGD9exFj2pQd94icnJzNxViTMHOz9qJ6BnwgNFH1IOA4EhF2QhADgAP/DGgwGdHd3h+x6DXWO/M9DvS1KSkpgtVrZtDP6VKYt2LFGFLE+3YTepPwZG3QK/OzsLKanpyGVSuFwOCJqLV+pyMnJQXFxcUDtRkpKCg4cOBCRKtPj8eCee+4BgKsRxIviJD4ihFwj3ieJHEsaUTgcvr0nRqMRXV1d7IlPEWtJrr+/n9nf8Z/UgYiiPDsZZVnJOKGxQiGXwOb0YFN9HmR2K9uGzusIZu8fDvwbnnptNjc3QyKRwOr0wCjPhNnuxJzLBsvJFuz09HR4vV5RRVFLAYlEAqVSycROMpkMEomEtZZTv83TLYriVz0CaTfGx8exf/9+dHZ24tOf/jRuu+023HLLLSGPeeDAAVRVVWFwcDCUF8WKwJLmKKzWhTdfY2OjaB6DLpcLHR0dSE9PD6hnCEQUcqkEX7ukCh/0aaAxO7E6NwXnlmfin/8cYY1ndLZItJ4QtNza29vLjiWRSGAw2/DXPi3cHgKZhMOoS4oLKspRq0qC0WjE4OAghoeHMT09HTAHIPS9o0WsYT8hhAm+aGs5bWajURTVbYjZzOZ/DmIgVHlUIpGgpKQEX/3qV+H1evH4449jeno67DEnJiZQXMxvHA3oRQEA53Ec1475rtJvEkK6I/8EsWFJlx5utxsej4eZ08Zy8/mDCrOqqqoWTP6iCJb4TEqQYlNDvs9rhBAcOXIEKpUKa9asieki9ng8mJ2dRWFhIdatW8cijEm9BWa7C6tO5kcSZBL0TJtRkjlfSaHCJ6VSyVywbDYb+53Q0YDLBX+i4cvLy8vL4XK5fJrZUlNTWbQhRNwUzTks5nFoQxhVwgo5ZqCX/f5/BEApIcTMcdwmAK8BENaGLSKWfOlx8OBBrFq1CiUlJaI9QWieI1x0IlQfYTabYbFYUFVVhby8vJjOzWq14tixY0hKSmJt9lTNKZXKwOGkuIvMEwo5eY607k8IgUKh8MkBzM3NYXZ2FkNDQ8w+j0YbYkKMiCLU/nK5HHl5ecjLy2NGNrOzs+jq6oLX64VKpYLb7Y6peiLW3FEgfHQWjRcF7YalL8HPi4IQYuT9ew/Hcf/FcVw2IUQr+I1EwJIRhdFohEajQUtLS8AZov6g4XqoL5kuDVwuF8477zxBU56EzAAZGBhAUlJSzCRBk5aVlZVQq9UghLARhBzHoTAzBWnTFszZPZBLOZjsXpxblgaPxwOPxwOv17vgqSORSJjrMzBfcfDXN1Cz3uVGJETDF0PRIUQ6nQ6Tk5M4ePAgkpOTGSFGkrcRkyjCIVIH7tbWVvT39yOUFwXHcfkAZgghhOO4jQAkAGZFPG1BWPSlBxU7abVaqFQqwRewhPdkDQTq90CH7Aq5GEIRBS3R6nQ6tLS0CJ4BEgw0ablhwwZ4PB7MzMz4kATHcUhOkOKymhwcV5vhdBNsVCWiUJkEr9cLs9kMtVqN1atXsxkjNMrgf9bExESf2j+d8EV9N6kJTzSmOYsdUYSCTCZDbm4uRkZG0NLSAovFskB6nZWVFbCZjY9oBlhHi0gjCplMhieeeAKf+cxngnpRYN4M9yscx7kB2ADcRMRKvESARY0onE4nSy6uX78eHR0dgvcN1ZZOJ5FTvwcq9w6X7whGFG63G11dXVAoFKxJjJZphVzo/O0IIT7DlaVSKSwWC5xOJ1wuF2Qymc8xUxUyNBcrfY5nNBrR09ODuro6pKamskiElpc9Hg8jDP5N4q9vOHDgADsfl8sl+ObiY7mTofQc/G3y+c1siYmJLLfhT4hLGVFE6m4FAJs2bQIhpJr/Gt+LghDyBIAnxDnD6LFoREFvZupLQcNpoQimuwhUqhTaGBaIKGw2G9ra2lBSUuJjzCtk6UO3ozeE2+1Ge3s70tPTfZKWcrkcaWlpOHz4MBQKBZMNB3rKT01NYWxsDOvXr/cxIKaSakoYhBDmPk5/7x9tSKVSFBYWBry5kpKSogrlI0GsRBHs7+8vvbZarZidnWWaBtropVQqRXO3EvI5zGZz0MbA0x2iEwU5OQdjfHzc52aOpB0cWEgUtLzodDoXVEsiteGnoG5ZdXV1C5ZEkZrXUNVmWVkZa4enSUuO41BdPf/QsFqt0Gq16OnpgdPpRFZWFrPLGx4ehslkQnNzc9DoiB9F0BwGzWcAvtEGH/ybiwqGqC+Fx+Nh0YYQazmhWKqlCx06TJvZqLz8xIkTbJlntVqj9qMQunyxWq1n5NxRYBGIwuv1wmq1YuPGjT5/3Ei/ID5ROBwOtLe3IycnBzU1NYL0EYHA95mgg4KCjQ6IxLxGr9fj+PHjPi5Z9Ob1bxJKTk5GSUkJe8rPzs5icnISbW1tSEhIQHl5eURt8wB8rO3p+9I2cZrf4JMHXzBUUlLCfCmmp6fR19fHXLBi9dyMtU8jmmhAKpWySAmYV+hSLwq73c5Ky5FMLhN6HnRWyZkI0YlCJpNh7dq1MQtdKFEIkU5H0sJOCEF3dzc8Hk/AoT4UQs1rHA4H+vv7GeH4VzZCEaRUKoVSqcTIyAiqqqqgVCqh1WrR3t4OYF7SnZOTI1j6TKMNmnMpKChgyxXgVKu0f8Th70thsVig1Wpht9tx+PDhqEVRsXZ+itEQJpfLkZ6ejsrKygWlZdrMlpmZGXJWSiQRRZwolhhSqRQzMzMwGo0BrfX4EPr0p/NFCgsLUVpaGvIiDmdeQ5OEbrcbzc3NEZMEMP8EonkcSoJpaWkoLy+H0+lkF7TFYkFGRgays7ORmZkZ8knIXwLxy7v8aIOfGKU3AD/aoIlDjUaDxsbGBQ7fQiXYYuQoxDS+8S8tOxwOzM7OYnh4mFnYBRoJEImxbnzpIRKEXDxerxcajQZSqZRVDkJBSERBk6sKhQJlZWVhzzMU+fCTlnRNT3MF9POF+4yzs7Po7+9HfX19wKdQQkKCzxyPubk5aLVaJrKi+Qa+yMpkMqGrqws1NTUBcy783IY/YQSLNvxFUXwJNsdxIWeJLlYyMxKEusn9hWz8WSm0gpSVlRURUZyJfpnAIhFFsK5HIdPE6RMxMTERubm5gr6gcBHF1NQUhoaGsG7dOhbWCzlmoM9Arf1p0rKtrY35cNL9wmF8fJwNUhYiVfZ/EtpsNuYeTrP8CoUCU1NTaGxsDPtUC5QQ9SePQN+hvwSbzhMZ5TWy0cFJMplsWXUYFELJJtislNHRUczNzUEqlWJ6ejqkvDyeoxAJ9Mkf7OanSsaamhpYLBbBybRg5VFCCI4fPw6LxRLWUs8fgciHnh8/aQnMP8mzsrIEWaX19/fDbrejubk5aiFQUlISiouLWZa/v78fIyMjkMvlGBgYYNGGkLJnsISoXq8HgJBiL/48EfpE1mq1GBkZYZJ9m82GxEThE8P4EKu0Gc0x+J9NrVZDp9PBZrOhs7MTAHyqRPSzRePAHc6PgoLjuFYA/wSwjcxb8C8plpwo3G73AkYmhGB0dBTT09MsKWi32xe0pQdDoPIov5N0/fr1EV+o/jkKvtKSn48oLCzE2NgYBgYGoFKpkJOTA5VKteDi9Hg86OzsRFpaGhoaGkRrVBoaGoLD4cCnPvUpSCQSlojs7OyE1+v1Kb8KTYjSJc66desgl8uDll/9xV70iQzMr//b2trY3yaaRjax/DJjbTz0er2sWkWb2fR6vc+slP7+frjd7ogUsEL9KDiOkwL4BeaniS0LliWi4MPj8aC7uxsSicRnZkek9v589yyLxYL29vYFPheAcCEVjSho0tJut7N8CV8pScU9Xq8Xer0eGo0Gx48fR3JyMnJycpCdnQ2v14uOjg4UFxcHHDkQDbxeL44dOwa5XI7GxkZGAjQRWVZWBpfLxcbwmUwmtjTIysoKevOMjY0xA2Ka0POPNmg+hv47kNhLoVBAoVCgpqYGMpksqkY2sZOZYh0j0KyUN998E2NjY7jgggtw5ZVXYseOHWHPPQI/ivsAvAKgNaYPEgMWLUcR8M38fDNtNhva29tRWFjo35cfkbs2fwQhvVEbGhoCOm2H6yHhb+d2u3HkyBGkp6czo19+0pKG43R7evHzS4xHjhyBxWJBfn4+k2PHevHTaCknJwclJSVBt5PL5T6jBmlClC4N/BOig4ODsFqtWL9+fcC/Tzixl9vtZtvwJfCBciw6nS5gIxs/2ljOpQcfodoDaJXogQcewNtvv4133nkHbW1tgr5jIX4UHMcVArgOwCU404giGPhRAjXSra2tDTg1nI6Ci+S4Q0ND0Gg0aG1tDZpwElpKdbvdOH78OFavXu2jtBRS2aAXDzXqaWlpgdVqZaVOpVLJSp2R5ilsNhs6OjpQXl4uyPOAf07UeQqYTxprtVr09/fDZrOBEILk5GQ0NDQITv4BgcVeNNrgl2P5x+Sb2Ho8HhZtDA4O+vRtLGUyMxQ8Ho+gfA8hBBkZGbjooosEHVegH8V/APgOIcSznDaCS04ULpeL3dD+o/z4iCSiAMBcoPhDkYMdNxxR6HQ6TE9Po6SkhJEEX2kpJGk5OjoKrVaLDRs2QC6XIyMjg5U6DQYDtFotBgcHoVAo2BIlkEKUD6PRiO7u7oDlz0iRmJiIoqIirFq1Cu3t7VAoFJBKpThw4ACSkpJYtBHunCj8o43x8XG2JKGkESghKpVKfRrZqLS8r68PVqsVCoUCer0+okY2PsToHl2smR5C/CgAtAB46SRJZAPYxHGcmxDyWsRvGAOWdOnBcRyGh4eRnp4e9oaOZA7I8ePHoVAoUFtbG3b7cERBk5ZFRUVISkqKWERFe1IIIQFDeP8OT7pEoWrRYMlHOqSmqalJNIMa2qpfWFjImploHwg9J7fbzc4pIyMj7OcnJ2eeGI1GnyRyMLGX/9+H37cxOTmJubm5mBrZxFq+hCMbu90ecSu/ED8KQkg5/TfHcc8AeGupSQJYwoiCuk2rVCpBN7QQoqDlytLSUszNzQk6j2BE4Z+0HB8f97m4hZAEddrKzMwMq/ykoD0XpaWlzBqOn3zMycmB3W5nCUaxLOJofqiqqop1YQK+fSB0Puns7CwmJibQ09ODtLQ0lhD1V2bScrTb7V6whPEXe/GXJ/T3/tEGXcLxna/9G9nCVXQWW7RFEY0qU6AfxYrAkhAFTTAWFhYKLlWFI4rx8XGMjY2hubmZNTUJQSCi4CstadKSto0LJYlocwd8+KsgDQYD+vv7YbFYkJ6ejpmZmahNaPigCs7a2tqwU85kMpnPORmNRmi1WqZe5LfM9/T0sMgu2N+LLlFoDiqU1wb/Jg/WyEZLlMH8NpeKKCJ1t6II50fh9/r/i/gNRMKiLj34rlGtra3Q6/UwmUyCjhFMRMVvN6cdqlarNaIKCf+4/kpLet7p6eno7e1lN2dOTk7QxiGDwcDa1UPNNI0EhBCMj49DqVSitbWVJR97enrgcrmQmZmJnJwcQcsBPugsimiWMBzHMa1EZWUl65UYGBiATqdDamoq8vLyBIXqgG9CNJDXht1uR0JCQsDWef9GNr7fJiHER369FJWTM7nPA1jEiMLtdqOzsxNJSUnMNSpSbYT/tnRNnZWV5dNuHonXBV9IRf0oArWHp6amorW1FU6nkyUebTYbE1UplUpIJJKgRjOxgJY/c3NzWfmMr8akfpL85UBOTk7YRq3p6Wlm+iOGWQ1NxE5NTaG6uhrJyck+SdpQBj2BwE+ITk1NwWAwoK6uLqzYi/Pz2+S7e9OoTMjfJxiELj3OVPk2sEhEYbVacfDgQZSVlfk4/kQ60ZwPk8mEjo4O5pjFh1CHK+AUqQRTWvovNRISEnwah3Q6HdRqNfr6+piWoqmpSTSSoLmDioqKoEsY6idJBT90OeCvj+A/4WgVJpQpTqRwOp1oa2tDWVkZO1d+9cI/AqIJ0XBP58nJSUxOTvqcayRiL/4SzmKxoKioCHq9XlAjWyDEiWKRiEKhUPiM8mNvJmBQcSDQ8YBNTU0Bv4xISqkSiYStr/2VluHyEXRNrlKp0N3dDY7jkJycjM7OTnaD5uTkRF2VmJubw7FjxwTlDij8lwN8fYTdbodKpYLD4QDHcVi3bp1o7lWU0Kqrqxk58ME36KER0NTUFJM85+TkICsra0FydmxsDBqNBuvXr/e5OSMVe/GRkZEBpVLp074/MjLCcj/8RrZAiC89FokoZDJZwOaYUIOKA4E2URmNxgXjAfmIREQ1PT2N1NTUgEpLIUlLahicl5fHlgX8G5R2dFLTGaE5BI1Gg8HBQaxbty6mZCXVRxQVFcHlcqG9vZ01dnV1dQW9QSOB2WxGZ2enYELzj4DMZjM0Gs0Cgx6NRgOTyRSQ0Ix2F372dj+OjBqgTJbjgStWo6VUCSCw2As41QHM//v7t+/7N7IFm8gmxFskHlGIhEiWHm63GzabDV6vN+B4QD6E3Ig0aUkNYKIhCYvFwqaR8UuKgO8Nyre44+cQgvVYjI2NYWZmhomzxADNEeXk5KC0tDTgDUqFXqHcnfxhMBjQ29srqJ09EPj5BNrOTZvYHA4HcnNzodVqkZWV5RNR/PDNPhwZm0OyXAKdxYUHXjuG/7l1PYpVSUG9Nux2OzweD1wuV8Dya6BGtkAT2YSIqeJLDxEhlCjoTS2Xy1FdXR2zjJeftJybm2PkEImHRDijGT6kUumCHIJGo8Hw8DDkcjnL1isUCp+2c7GWBdRjtKSkhDXFBbtBaX9HqM5X/7/BunXrRMvJyOVyGI1GqFQqVFdXs3wL36BHlZmFI6NzSFHM3+yJEg42pxcdE0YUq3yjL0oaDocDPT09WLNmDcthUd8QjuMCir34RjZTBhsmtQZYZ3SwWq1oa2sL2chmtVoDLsHOFCypMlMIUdAekPr6evT29goutQWDf9LSaDQym3shUQQ9xtTUVFTVAn4OoaqqCjabDRqNBt3d3TCZTEhOTo55tikfVqsVHR0dQXMHFP5JWv/OV5pvoUsUWjERU/RFCGEdsKtXrwbHcQENevqP90HidcHmkEAhl0HCcQAHpCkCX752ux1tbW0L/gZCxV7v9mjwXq8WEg6QSji0pCbinDVrfBrZqB8FrX5Fk6PYu3cv7r//fhw/fnwAAbwoOI7bAuCnALwA3AC+Rgj5OKI3EQmLFlEEc0gKBir9nZmZYT0gNEkZDVEEaw9PSEjAwMAAgPnQO9STUSyjGT6SkpKYGUp5eTkSExOZEjMjI4OV8aJ5L5oMra+vj8hAJVjna0dHB7xeL+RyOZxOp0/reazwer3o6upCamoqysvLA14b/JLwAwnTePidfphsLgBeVKrkKE+a9yzhkzcliTVr1ixoNhQi9po2OfHuMQ1yMxSQSThYHB68N+LB1k8n+jSy0X6dgYEBHDhwgBkuCQX1onjvvfdQWVlZi8BeFH8D8AYhhHAc1wjgZQBrBb+JiFgR5rper5dVEaL1pOCDKi0zMjJ8kpZerxeZmZlobGz06a+gT09+uczj8aCrqwspKSmiGc0Ap574lZWVyMnJAQAmUpqbm2NJzcTERJZDEBLF0Is21mQolU2npqaitLSUiamSkpJw6NChmDpfKTweDzo6OpjUXQiuqs9HSVYKuqdMSE+U4ZzCRBj0Oh+DnrS0NAwODgpqmgsm9jLZPeA4Agkh8Hq8SJQCTi/gcHuRJJeyffgjAdLT0/Hhhx/iiSeewOOPP45du3aF9WWlXhQVFRUghDgDeVEQQsy8XVKwsLN0ybDsREE9MvPz8xdMOI+UKAghrGznr7TkJy35pTuXy8XWxBaLBZmZmcjIyMDIyAiKi4tFnfxEn/iBFJz+ng0WiwUajQadnZ0ghIS07p+cnMTExIToywLat7Fx40YmVIu285WCknheXh6KiooiOqfaVWmoXXUqUlJmpDOB1eTkJLq7uyGXyzExMcEqT0I1IzTaKMxMmbfx8xAkyiSY1JmQnSyFlHjgdpOA5de1a9ciMzMTv/jFL7B69WpBxC7EiwIAOI67DsDDAHIBfEbQh1kELOnSg4LesNQZO9jMjkiIQiqVsvZkvtIyXNJSLpf7lMuo2lEul2N2dpY9PWIVKanVamYvJ+SJT/sa6I3AJzO+OnRsbAx6vV60pRFwyj3Lv2/Dv/PVarWyfEuozlcKl8uFtrY21t4uFpxOJyYnJ7FhwwakpaUxg57h4WHIZLKAArRgyExJwBfPLcafD05AM2eB1GXF/Z9pYsuVYLkN6sAtNE8h0IsChJBXAbzKcdyFmM9XXCboDUTGkkcUtLY9NTWFkZERn7GD/ohEcel0OqMexENBuyTPOeccJCUlwWQyQa1WL6hWRJrxHx0dhUajiXqN709mNPFIhV4VFRWi+C4Ap7w9lUpl2PA5OTkZpaWlITtf6RySQCpOMWCxWNDR0eGTl6EGPVVVVQEFaFRgFay6U1uQjq9d4MKxvn6c27IRibwIIZjYi/a6CIVALwoGQsjfOY6r5DgumxCiFfxGIoELUyOOek3kcrkC3uQHDhxAamoqHA4HGhoaQj6p+/r62FMq6AmeTFpOTU1h48aNSElJiZgk+EYzjY2NAW9mWq3QaDRB8xqBjnv8+HE4nU7U1dWJVv6kiUCaGNVqtdBqteA4zqeBLVJQgdaqVat8BjZHCr7tHo3KrFYrqqqqRF3KUeFXQ0ODoJuUT7IGgwGJiYkBDXpo45yQKpfX68W+fftw6623MntCIXC73aiursbf/vY3VFRUKAAcBPB5Qkg33YbjuCoAgyeTmc0A3gRQRKJxyYkRi0YUbrc7oDP2Rx99hIKCAkElwYGBAaSlpflMvPJ/D5q0NJlMrDEpUqOZvr4+eL1e1NTUCLqZ6VJAo9GwvAa/UQzwTYZWVlaKlgwN1DBG4XA42HnZ7XafDtNwn4tqL8R+4ttsNhw9ehTZ2dkwm80xdb7yQVvloxV+8Q16tFotPB4PMjMzoVAoMD4+jubmZkG5hsOHD+O+++7D66+/LjgxS7Fnzx587WtfQ39//wnMe1H8nO9FwXHcdwDcCsAFwAbgW8tVHl0yojCbzejo6IBUKkVdXZ2gJ8DQ0BATwfjDvz28o6MDZWVlSE5OFqy0jMZoxh/0KaVWq2EwGJCamorMzExMTEywcppYsNvt6OjoQGlpaVDypPB4PNDpdNBoNJibmwtpOBOubyNa0Cc+P3lLw3StVhv2vIKBWgJGSxKB4Ha7MTo6ipGRESgUCh/H8mDn1d7ejrvuugt/+ctfUFlZGcvbL58ZpkAsCVFQK7PGxkaMjo6iqKhIUI/A6OgoOI5b8OQM1B7e3d3NynZCogJqNOM/ozMWEEKg0Wjw4t+78fIgYPcAzYUpeHhrLbIzYpP3Uvl4IH2AkPMymUzQaDRsKUDzLTQnIaaXBnDqiR9qWcDvfKXnFS7xSCtHYloCAqek6evWrYNCofA5L75BD5W7d3V14Y477sDOnTuxZs2aWN/+7CYKt9vNjGuampqQkJCA3t5e1vcQDhMTE3C5XD5JNaq05MuIPR4PqwoAYDdBsIuNGs1E0qUpBAaDAXv3d+NXRz1wewkkICAA6rOluHedQlBeI9T5Cl2Lh4PdbodGo8HU1BRMJhNWrVqFgoKCmJYC/udL+0EiuZlp4pFOUvf3/qDHbWpqitnlK9D5BpOm0x4QjUaDjz/+GB9//DG6urqwc+dONDc3i3EKZy9R0ElRCoUCa9asYU/5/v5+ZGRkCFoHT09Pw2KxoLKy0kdp2dDQELQ93Ol0QqPRQK1Ww+FwsJuTluyoFLmxsVG0fgUAmJmZwfDwMLpcufjtR6PzMmPMPzWlEg6ffP08n7yGSqVCbm6uT14jEDQaDWuxF/N8qUCroaGBaTaMRiOrVvg3ZQkFTQTG2g9CxxrSxKNcLofNZsO6detEHQQ8NzeHnp4ewefb3d2Nr371q1i9ejW6u7vxs5/9DFdffXWsp3H2EsXY2BicTucCUc3Q0BASExMF1dHpRVJeXs6SljQxKKSyQY1hafsyvSHXr18vmjAJAEZGRjA7O4uGhgbsPqbFQ3v74fHO50k8XgJVshzvf+18tn2gvAYVLvGrQBMTE5iamkJTU5No0mlgnoDHxsZYlEdBqxUajQY6nQ4JCQkRCaooqdHwXSzQ/p+cnBxmohzOnlAIKEkIjVCGhoZw880345lnnmGRhBienDibiYI2XvmDmsYIUeXRVm2z2RxSaSmkstHd3Q2v18vmRKSmpiI3NzcmMRWNctxuN2prayGRSGB3efBvzxzBqM4Gt5dAJuHw8JYaXLo2cNmMnz/QarWsY9Jut8Nms6G+vl40IRVwyhimsbEx7OemVQFaEqbq0LS0tAV/c0o+dF6pWKCenHzyoZ2vGo1GcOerPyIlidHRUWzbtg2///3vsXHjxpg+UwCcvUTh9XqZYQofgfIOwTA6Oor+/n60tLT4KC35I/3CIZDRDL051Wo1Zmdn2ZOTtn4LAX/ocEVFhc+NY3d58M4xDYx2FzaUKH1kx+FgtVrR3d3NpoBHm9fwBzk50NhkMgmeBsYHnWOq0WhgNpuhVCqZoGpqagrT09NoamoSzWYPAJOKh4oAaXSm1Wqh0+lY52uoHhmj0cgSokJIYmJiAjfeeCOeeOIJXHDBBTF9piCIE4U/+HmHUBgfH8fIyAiSk5Oxfv16AL7TuoTcNKGMZvigUmSNRgNCCHJycpCbmxs0ERdocI4YoNqLtLQ0lJeXw+12R5XX8AcVfnk8Hh9T4mhBez40Gg1mZmZACEFlZSVyc3NFW9JpNBomeRd6TKqNoNGZ1+tdQLSRksT09DRuuOEGPPbYY7j44otj/FRBEScKf2g0Guj1elRXVwfYy7c9vLKykoWdkcqxaVJNiNEMH+GSoZR8Vq9eLahyIxThVJFC8xqB9qN9G1VVVaIJv4D5Nfvc3BwqKytZtAHEnj9Qq9UYGRmJeRlDhXFarRZmsxnJyckwGo2CE6JqtRrXX389fvGLX+Cyyxa1xeLsJQpCCJswzoder8f09HTA3n2+0rKyshJOpxOdnZ1Yv369z+zKcJiYmMDk5CQaGxtjSqr5J0Np/0djY6OoZVXaQVtRUSFIAhworxGoD4W2c6tUKkFLPaEghGBgYAAOh4PlZij4+YNA4w3CYbFyHUajkbW2G43GsJ2vs7Oz+OxnP4uf/OQnYlQ1wiFOFP4wGo0YGRlBQ0ODz+tUaVleXs7s21wuF/75z38y/YCQno2BgQHYbDbU1dWJmgSkTuDp6ekwmUyiJEOBU+rFWAYP+/eh0IlZAwMDKCgoEFUdSiM+QgjWrl0b8juh4w1o9Yq6b2dnZwckgampKUxMTGDdunWi5jqo+Isv0vJP1FJZeXp6OgwGA66//np873vfw+bNm0U7jxCIE4U/LBYL812kCKS0pEnL6elpzMzMwG63L1gG8LFYvRXUeUun07FKQaBkaG5uro91nBDo9Xr09fWhoaFBNCmyy+XC9PQ0BgcHfRSYkVQEgoFa1yUkJES8jOGb+2q12gXjDSYnJ1kpWEySoEQcSu5NI0e1Wo0vfelLsNvt2Lp1K374wx8ulWHu2UsUwLzoyh92ux3d3d3YsGEDgOBKS/98BFVfqtVqmM1mlthTqVQsuVhUVCRqclFow5jVaoVarWZr9HDJUGBeoDUyMiK68Ivft6FUKploiZaEheQ1AkGIdV0koCpMWkUBgNraWmRmZopG8kJIgg+TyYTPfe5zuPzyy2G32zE1NYWnn35alHMJgzhR+MPtduPIkSNobW0VpLQMBJrYm5mZgV6vh9PpRFlZGUpKSkRbbtD1fUZGRkQ3Bu3gVKvVcDqdyMrKQm5uro/2YGxsDGq1OmhLe7QI1IRF4d/vwZ/dGY6o6N8iKysLJSUlop0vcErXUVBQAJ1OxxrFQo03EIJIScJisWDbtm34f//v/+HWW2+N6j1jwNlNFE6nc4GTDyEE+/btg0Kh8FFauj1eaM12JMslSE6QCR6aMzAwgLKyMphMJszOziIlJQW5ublRPTUpaMt1cXFxTE5M/JCWag/ohPRotAyhQPsVhC5jAuU1cnNzF+SCaII5Pz9f1FwHMK+ToUs6+rfgjzeYnZ2NyjCImtkI7Y2x2WzYtm0bbr75Ztx+++0xfaYoEScK/+NbrVZ8/PHHaGxsZEnLSYMNP3qzBxqzExyA2y8owdV1wTs6CSE+CkP6VKbrYLVaDa1WG1XuwGw2o6urS/SWa7fbjY6ODhZlifHUpKDCpGj7Qfz9NWilIjU1FR0dHSguLmbflVgYHh7G3NxcWMKkhKbVauFyuZg6NJjdXqQkYbfb8a//+q/YvHkz7rrrroiXPdRy3+Px4I477sADDzzg8/tHH30Uzz//PID5a6CnpwcajQaZmZkoKytDWloaurq62gG4CSEtEb35EmJJiYImLb1eL/7lX/5l/g0Iwf3/24FRvQ2qZDncHi9MDg/+/bpaVOYsfDJGYjTDzx1wHMdII5jQJlrtRTjwy5TU3IQ+NbVaLRQKRVTJUCB430a0oMu66elpTE9PIz09HcXFxRH5RYQDVYjW19dHFFX5l6v9xxtEShJOpxO33HILLrvsMnz1q1+NmCQ8Hg+qq6vx3nvvoaioCK2trXjxxRdRW1sbcPs333wTv/71r/H+++8DAMrKynDo0CFkZ2ev+IhiUT0z+Qa7/EE8R44cAXBSlOV248SsFTkpcnAA5FIJCDwY09sWEAU1mqG6gHBfbHJyMsrKylBWVsZaq48dOwaPx8PCWXpBTU1NYWxsLKohP6HAT7TylzH8oUC0e7O9vR0cx7FzC9eizR/oK1alQCKRICUlBSaTiZGPRqPByMhITL6hFHQyWaQkAczPMKVTyv3HG8hkMlitVsGmSC6XC7fddhsuvPDCqEgC8LXcB4CbbroJr7/+elCiePHFF3HzzTdH/D4rAYturhtoEA9wqrIhlUiQl6aA0e5CmkIGj5cABMhO9X06xmo0k5iYyIbJuFwudoHZbDbmsLx+/XpRk4u0AhFOxcl33HY4HNBoNOjr6wuaDCWE4MSJEzCbzaJOKAdOzR1Zu3Yt03Wkp6ejsrLSZ8pZqLxGIBBCMDg4CLvdjvr6+pgrG/zxBlarFUePHkVeXh6GhoYwODgYskfG7XbjzjvvxIYNG/DNb34z6nPxt9wvKirC/v37A25rtVqxd+9ePPHEE+w1juNwxRVX4OjRo4cBPEUI+e+oTmQJsKhEQe3Z/QfxcBwHrVbLSmHfvKwSO3b3wWBzw+Ml+Ex9Lup4jVRiG83I5XIUFBQgPz8fx44dg8PhYANuou2n8AcV+UTqHKVQKNiwYxpqj4yMsJJwTk4O1Go1CCFobGwUVZIdqmoCzE/u4s9DmZ2dXTBCIJBegwrhXC4X6urqRD1nKtRrbGxksuxg4w1UKhUIIbj77ruxZs0afO9734u50c4fwY735ptv4oILLvDJe33yyScoKCgAx3FXA3iP47heQsjfoz6hRcSiEsXAwAAKCwtZIoyWP6urqzE5OeljYvPbbfWYmHMgLVGG0sxTIff09DTT/YvpakSnfSuVSnbx8suufX19SE9PR25ubsRTsegw31jt2vxDbZ1Oh97eXpbUU6vVyM7OFqUkTH0oha7v5XI58vPzkZ+fv2B2KdVr0ETt8ePHWU5pMUiirq7Op3cj2HiDu+++GyMjIygpKcF9990X87n4W+6Pj48H1fG89NJLC5YddFtCiJrjuFcBbASwIoliUZOZfMv+QB4ShBDW6KTX65GWlsZk0RKJhDUdhbP1jxTUpDZU+ZOauKjVata+LKTsSh20Iul6FAL/hKhYyVAgeuu6QPDvQ3E6nUhKSkJdXZ2oRE+XdbW1tYIiNq/Xi69//etwOp1YvXo13nvvPezZsyemz8u33C8sLERrayteeOEF1NXV+Ww3NzeH8vJyjI2NsdK1xWKB1+ulS8pUAO8B+AkhZG/UJ7SIWFSioAa74aZ1Aafq5zMzM5idnYXb7UZKSgrq6+tFveFoeL127VrBJrWByq5Ufck/N+p0JcQUJhLQzlI6fdwfNBlKqztCk6HAqegnVus6fxBC0NPTAwBITU0Nq9eIBJQkampqBC1FvV4vvvvd78Lr9eI3v/mNqDkdarnv8Xhw22234fvf/z5+97vfAQDuuusuAMAzzzyDvXv34qWXXmL7nThxAtdddx0AoKOj4xiAFwghPxftxETGohOF0+mMqD2cVgnoLAqtVovExET2xIwl2ShW+dO/7JqTk8OeEP7dlLEi0nkbNBmq0WjgdDpZUi+QK1U0ng9CEKwnhOY11Gp12LxGMERDEjt27MDc3ByeeuopUb8bEbHiy6OLShTPPvssKioqsG7dOkHr6GBGMxaLBTMzM9BqtZDJZMjNzY3YJIWWP5uamkQtf9JqjNPphEIx77Sdm5sbk5cjBa1ARCv+8leG8hO1arV6Udq5qf9FUlLSAucv/+34fShUURtKr2G329HW1iaYJAgh+PnPf46JiQk8/fTTUeVywgmqPvzwQ2zZsgXl5eUAgM9+9rP44Q9/KGhfHs5uonj11VfxwgsvoK+vD5dccgm2bNmC1tbWgKxOn/b+iSl/2Gw2zMzM+IiocnNzg4bN1ALOaDSynhKxQNWW2dnZPpPR1Wo1bDYbC7ODqQhDIdqqSTDwjW+ok1d1dTVyc3NFHW5MG8eotkAI+Es7/7kjNK9BSYJftg13zEcffRT9/f149tlno1oKChFUffjhh/jlL3+Jt956K+J9eVjxRLGoVY/rrrsO1113HWw2G/bu3Ys//OEPuO+++3DRRRdhy5YtOO+88yCVSnH8+HHMzc0JEjslJSX5iKjUajUzzs3JyUFeXh67uLxeL3p6eiCRSFh5VizQJUFJSQmr6vCz7R6PB7Ozs2xobyRlV35yUaz2c4lEgqysLFgsFthsNpSXl0Or1WJkZASJiYnsxox2CeL1etHZ2YmMjIyITXI4jkNaWhrS0tKYXkOr1aKnpwculwtKpRJarRa1tbWCSeLxxx9Hd3c3XnjhhajzRZEKqsTadyViSaaZJyUlMdJwOBx477338MILL2D79u3sBvrDH/4Q8ZIgMTGR1fWdTifUajV6enrgdruRmZkJnU6H3NzcqMcFBgNdIoVaEkilUhbtRFJ2pX0bi5FcHBoa8hFpqVQqrF69OqAyNDc3V3CVwuv1LpCox4KkpCQmjjObzTh69CiSk5PR29sbNq9BCMGTTz6JAwcO4OWXX45pWSVUULVv3z40NTWhoKAAv/zlL1FXVxeRGOt0wJIQBR8KhQLXXHMNLr/8ctx8881ISEhAWloa/uVf/gUbN27E1q1bcdFFF0X8ZEtISGBCJbPZjLa2NsjlckxPT8Ptdi9QN0YLOtKuvr5e8CAa+jTPysryKbsODg4iOTmZPc01Gg3Gx8dFnztCBU9OpzOgdDqQMpQ+zWnOJViVgpZts7OzF4x+jBUOhwNdXV2or6+HSqVaoNfwz2sQQvDHP/4RH3zwAf7yl7/E/DcUIqhqbm7GyMgIUlNTsWfPHmzduhX9/f0RibFOByw5UVDI5XLce++9uOSSSwDMr/c/+ugj7Ny5E9/73vewfv16bN26FZdccklET1b+2l6lUrGE3vDwMCwWC8sbRDM+j07XikX8xXEclEollEqlz9p837598Hg8Ea3thYBK6IF5Y5hwn9lfGarVajE8PAyz2bxgarvH40F7eztyc3MFzWmJBA6HA0ePHvWZtepPuPRvNzo6iueeew4OhwOjo6N45513RElYCxFU8fNHmzZtwt133w2tVhuRGOt0wKImM6OFx+PBJ598gldeeQXvv/8+amtrsXXrVlx++eUhtQFUExDMk4FO+J6ZmfHJG6hUqrA30OTkJCYmJkTr0qSgfRt0hIG/HiKSJYA/xHTgDuQCbjabUVRUtCiRRFtbG1avXi242vPb3/4W//u//wulUgmLxYK//vWvMed3hAiqpqenkZeXB47jcODAAdxwww0YGRlhycxwYqyTWPGhxookCj68Xi8OHjyInTt34r333kNVVRU2b96Mq666yif0j/RG5ucN5ubmmJTcfxo6IQTDw8MwGAxobGwUtWpCn/aB5M10CaBWq+F2uyMuu9LkYnp6OivdiQWXy4XDhw9DoVDA4XAwnUt2dnbMJOp0OnH06NGISGLnzp14+umnsXv3bkZgYtkEhBNUPfHEE3jyySchk8mQlJSExx57DOeff37QfYMgThRiwuv1oq2tDTt37sTevXtRXFyMa6+9Fm1tbbjoootw9dVXR3UjE0JgMBiYtR6VklMna6/Xi7Vr14oq1qFjDpOSksKaAUdadl1M6zq3242jR4/6mNkEUoZGEwlRkqiqqhI8M+W1117Dk08+ibfeekvUEQpLjDhRLBYIITh69Chuu+02OJ1OlJSUYOvWrfjMZz4T02AeKiWfnp7G5OQkFAoFKioqYrLW8we9kTMzMyOuEtCyq1qtDlh2XUzrOtoNXFJSErTVnx8JCUmGUkRDErt378avf/1r7N69W7Acf4UiThSLiZdffhkTExP42te+hr6+PuzatQtvvvkm0tPTsXnzZlx77bXIycmJeG3ucrnQ0dHBbkCqChVDSh6ubyMSBMobmEwmlJWViZ44c7lcOHr0qGApOQCWDKWSbf9kKIXT6URbWxsqKysFk8S7776Lhx9+GHv27In6wRBOOfn888/jF7/4BYD5fpUnn3wSTU1NAMBs7KRSKWQyGQ4dOhTVOZxEnCiWGtQg5ZVXXsHrr7+OhIQEbN68GVu2bEF+fn5Y0qCdpYFuCIvFwpSN0UjJaZKuvLxc8M0mFHa7HUeOHEFKSgpsNhuSkpJY3iBWiTa9kcvLywVNMgsE/2FA6enprAels7OTRW1C8MEHH2DHjh3YvXt31H9HIcrJf/zjH6ipqYFKpcLbb7+NHTt2MC0Ez8Yuqvf3Q5wolhOEEIyOjuKVV17Bq6++CkIIrr32WmzduhVFRUULSIMa6/JLcsEQqZQ81r6NUKDyZnpsQohPf4xcLmeRUKRlQ0pukSwJwoG/vBsfH0daWhqKiooEJUM/+ugjfO9738Pu3btjMvzdt28fduzYgXfeeQcA8PDDDwMAvvvd7wbcXq/Xo76+HhMTEwDiROGP05oo+CCEYGpqipGG1WrFNddcgy1btqCiogK9vb2Ym5uLqrOU+nGq1WomJecPABK7b4OPQNZ1/rDZbFCr1VCr1RElG6mWYTHIjS5lysvLkZyczFr4Q5kg79u3D9/85jfx1ltvxZx/2bVrF2srAIDnnnsO+/fv97Gq4+OXv/wlent72fbl5eWsrP7lL38Zd955Zyyns+KJYtkEV0sNjuNQUFCA++67D/fddx9mZmbw6quv4hvf+AbGx8fhcDjwzDPPRFV75/txUik5daJKS0uDXq/HunXrROvboAhnXUeRlJSE0tJSlJaW+igvQ5VdaZQiJLqKFHySoEuZ8vJylJeXswlifGWoyWSC0+nEN77xDbzxxhuiJGkjUU5+8MEH+OMf/4iPP/6YvUZt7NRqNS6//HKsXbsWF154YczntVJx1hCFP/Ly8nDXXXchKSkJv//977Ft2zY88sgjmJ6expVXXonrrrsu7DiAQOBLyaempjAwMICUlBR0dnaym1IMKXmk1nUUfOUlLbtSk2FadpXL5WGjlGhBKydlZWUB8x2JiYk+5zczM4P77rsPHR0duO666zA6OipKyVeocrKjowN33HEH3n77bZ+lF902NzcX1113HQ4cOBAnijMZGzZswI033oikpCTcd999MBgMePPNN/Hzn/8cw8PDuPzyy7F161Y0NTVFRBpTU1MYHx/HOeecg4SEBFGl5LS7NFZPzkDdrsPDw9BqtcjOzmZO6WLpR9xuN9ra2lBaWiooCSmXy2EwGKDT6fD3v/8d4+Pj2L9/Pz71qU/FfC6tra3o7+/H0NAQCgsL8dJLL+GFF17w2WZ0dBSf/exn8dxzz6G6upq9zrexs1gsePfdd5kHxZmKsyZHEQ1MJhN2796NV155BX19fbj00kuxZcsWtLS0hLx5RkdHodVqg1rixSIlXyzrOuCUWW1NTQ08Hg8zleF7mUarTKVCrVAaDH/09PTgi1/8Il566aVFac8Op7q844478MorrzCtCy2D8m3s3G43Pv/5z4dSXQrBis9RxIlCIOhchl27dqGzs5N5apx77rns5uH3bQgdcCNUSg4AarUaw8PDolvXAadG8fl3xdIKBTWViabsGg1JHD9+HLfeeiuef/55NDQ0RPWZTiPEieJMhN1ux3vvvYddu3bh8OHDOP/887F582bs3r0b27ZtQ2tra1Q5iGBScmrNPz4+Lrp1HXAqKRou30HLrrRCQbUkocqudLlRVFQkuJw5NDSEm2++Gc8++yzWr18f1WcKJ6YihOD+++9nTtzPPPMMmpubBe27CDh7iEKn02Hbtm0YHh5GWVkZXn755QXZ8r6+Pmzbto39/8SJE/jJT36Cr33ta9ixYwd+//vfswTXQw89hE2bNkX0YZYDTqcT7777LrZv346kpCRs2LAB1113HS688MKYnvr8J/nU1BQIIaiqqkJeXp6oDt+0dBtpUhQ4VXal1npUS0LLmtGQxOjoKLZt24Y//OEPaG1tjfjzAMLEVHv27MFvfvMb7NmzB/v378f999+P/fv3RzxPVCSseKIQrcvpkUcewaWXXor+/n5ceumleOSRRxZss2bNGrS1taGtrQ2HDx9GcnIyW+sBwPbt29nvTweSAOarHIODg7j77rtx5MgR3HLLLXj77bdxwQUX4Mtf/jLefvtt2O32iI/LcRwyMjKQkJCA1NRUNDU1wW634/Dhw2hra8Pk5CRcLldM5240GtHV1YXGxsaoui1p2bWlpYV11vb09GD//v0YGBjA4cOH2UQ2IZiYmMBNN92EJ598MmqSAHxt6BISEpgNHR+vv/46br31VnAch3PPPRcGgwFTU1OC9j0bIdqj6fXXX8eHH34IAPjCF76Aiy++mOnkA+Fvf/sbKisrRbFOW27wh9xecskluOSSS5inxq5du/CjH/0I9fX12Lp1Ky677DJBlYpA1nVKpRKVlZUs/D969GjUruTUqSvWygkFv+zqcDhw+PBhSKVSjI6Owmw2h63wTE9PY9u2bXj88cdZm3a0EGJDF2ibiYmJM87CTiyIRhQzMzNs6taqVaugVqtDbh9oxNoTTzyBP/3pT2hpacGvfvWr06YjMNDFL5VKceGFF+LCCy+E1+vFgQMHsHPnTjz00ENYvXo1tm7diiuvvDLgk5wQgv7+frhcLjQ0NCw4fkpKChMo0fCf+l2Gk5IDp8qrYo9pBObD/q6uLtaYRsuuExMT6OnpgVKpZBUemqxVq9X43Oc+h0cffRQXXXRRzOcgREwVbJszzcJOLEREFJdddhmmp6cXvP7zn0c24MjpdOKNN95g+noA+MpXvoIHH3wQHMfhwQcfxDe+8Q08/fTTER13pUIikeDcc8/FueeeC6/Xi6NHj2Lnzp147LHHUFJSgs2bN2PTpk3IyMiAx+NBb28vpFKpIOs6vuqSSsn5ruR8KTkw37PQ19e3KOVVj8eDtrY2rFq1igmS/E2GDQYD1Go1jh8/ziKmP/3pT3jooYdw6aWXinIeQsRUwbZxOp1nlIWdWIiIKP76178G/V1eXh6mpqawatUqTE1NhRTUvP3222hubvYplcnlclx11VUYHh5GXl4eZmdnA+4brL1XSDJ1JUAikWDDhg3YsGEDHnroIXR1dWHXrl245pprkJ2dDZvNhmuvvRZ33313xE8yfym5RqNhUnJamaBDf6Ilib/1afDI3gFYnB5cXpON719VjQTZKf/M/Pz8oDeWRCJBZmYma1xzOp344Q9/iLm5Ofz3f/838vLyoq5y8CFETLV582Y88cQTuOmmm7B//35kZGRg1apVyMnJCbvv2QjRkpmbN2/Gs88+C2B+QtiWLVuCbvviiy8uWHb84Ac/YMnQcLLhDz74AG1tbT4eAEKSqSsNEokEjY2N+MlPfoJPPvkEwHyEsGvXLlx33XV4+umnoVarA4bD4ZCQkIDCwkI0NzejubkZbrcbx48fByEE4+PjMBqNER+3fXwO33rlGCbm7DDYXHizYwY/e/s4I4m8vDzBfRhGoxE//elP8d3vfhcDAwP42c9+JlYnJmQyGZ544glceeWVqKmpwY033oi6ujr87ne/Y4KqTZs2oaKiAlVVVfjSl76E//qv/wq579kO0cqjs7OzuPHGG5kWf+fOncjMzMTk5CTuuOMO7NmzB8C8cKm4uBgnTpzwsS5LT09HQUEBEhISkJ+fj8HBQQwODi54n2DtvWvWrMGHH37IIpqLL76YuU+fDpicnMTevXtx2223MU+NXbt24Y033oBCocC1114r2FPDH7SfY/369ZBKpdBqtZiZmYlYSv4f7w/i9x+P+lwUGUky/NdlqcjJyRHsxG0ymXDDDTfg3nvv9SmXn04ghIiZu1jxSZAVI7hSKpUwGAzs/yqVCnq9fsF2wdp7he5/uoEQgpGREdYez3EcrrnmmqCeGv4INYiYSsnVajWMRmPYaWZ//GQE//nhEFyeU5dFViKHFz9fJZgkLBYLbrzxRtx222245ZZbBO2zUrBr1y4cO3YMW7ZsYU5XIiFOFHyESoZ+4QtfEHSjT05O+rT3/uY3v8GFF154xhIFH4QQTE5OMtKw2+3MU6O8vHwBaVDJ9/r168OqOYVIyQ1WF7b87gAMNhdcHoIEKfDtC/Px+X+pEXT+NpsN27Ztw+c//3ncdttt0f0R/CAkNzU2NoZbb70V09PTkEgkuPPOO3H//fcDgGCh3yOPPIIf//jHKCgowOjoKN566y1ceeWVYjXNxYlCKKJZOuzYsQOpqan45je/edovPSIFIQRqtRp/+ctf8Je//AVzc3PYtGkTtmzZgurqauzfvx9SqTQqyTeVkqvVauh0Oh8pucnhxcuHJ3BibBIXV2fjqpbq8AfEvOz985//PLZu3Yovf/nLooXt3/72t5GZmYkHHngAjzzyCPR6/QL9ztTUFKamptDc3AyTyYQNGzbgtddeQ21trc81FAydnZ148MEH8dWvfhUtLS346U9/it/97nc4dOgQ1qxZI8YyZMUThXj+8zFCSDLUYrHAZDKxf7/77ruor69n+z/55JO4/PLL0dDQAIfDETCiGBsbw6c//WnU1NSgrq4Ojz/+OPvdjh07UFhYiHXr1mHdunUsr7ISwXEc8vLy8JWvfAXvvfce9uzZg+LiYnz/+9/Hhg0bcM8990Aul0cl9+Y4DiqVCmvWrMG5556L4uJizM3N4cCBAxjpP4ZG+TTuOr9AMEk4nU7ceuut2LRpk6gkAcwL/b7whS8AmBf6vfbaawu2WbVqFevjSEtLQ01NDbO0C4c333wT//mf/4mCggKcf/75SE9Px49//GNcddVV+MpXvoKpqamzQ2dBCAn1s2TQarXkkksuIVVVVeSSSy4hs7OzhBBCJiYmyNVXX00IIWRwcJA0NjaSxsZGUltbS372s5/57F9SUkKysrLIJZdcQh588EHy7W9/e8H7TE5OksOHDxNCCDEajWT16tWku7ubEELIj370I/Loo48u9kddVPz5z38mF110EXnqqafI5s2bybp168i3vvUt8sknnxCTyUQsFkvUP0ajkfz9738nH330EXn//ffJJ598Qvr7+4nBYAi6j8FgIFu3biWPPvoo8Xq9on/ejIwMn/8rlcqQ2w8NDZHi4mIyNzdHCJn/zktLS0lDQwP54he/SHQ6nc/2jzzyCOE4jnz2s58lhBDi8XgIIYTs3buXXHrppeThhx8mVqs11o8R7j5c9p8VQxRioLq6mkxOThJC5gmhuro67D6bN28m7777LiHkzCCK9vZ2YrFY2P/n5ubICy+8QK6//nrS0NBAtm/fTj788MOIScNkMpF//OMf5NixY+w1tVpNurq6yAcffEA++ugj0tfXR/R6Pfv93NwcueGGG8jPf/7zmEji0ksvJXV1dQt+XnvttYiIwmQykebmZvLKK6+w16anp4nb7SYej4d873vfI1/84hcX7Ld9+3aiUqnI+Pi4z+vf//73SX19PTl06BAhhBC32x3tR1x2Igj3s2JyFGIg0oTm8PAwLrzwQnR1dSE9PR07duzAM888g/T09NNORi4EVqsVe/bswSuvvIKuri5cfPHF2LJlC84555yQhjRer5f9jcrKygJuwzfw1el0OHToEPr7+1FTU4MdO3YsWnguNDflcrlwzTXX4Morr8TXv/71gMcaHh7GNddcg66uLgBgicqZmRncdNNNSElJwXPPPedzTVx66aWYm5tjmp4ok5srfu2yYnIUQnHZZZehvr5+wU+kHX5msxnXX389/uM//oMZ037lK1/B4OAgkyF/4xvfWIyPsGxITk7GDTfcgBdffBEHDx7E5Zdfjj/96U8477zzsH37dvz973+H2+322UcISQCnpOStra1obGxEW1sbDh48iL/97W94/vnnF+0zCcltEUJw++23o6amZgFJTE1NsX+/+uqrLOcFgN3weXl5uPvuu2EymfDnP/8ZLpeLidVeeOEFlqjl73PGIUzIcVpB6NLD6XSSK664gvzqV78KeqyhoSFSV1e3KOe50uBwOMiePXvI7bffTmpra8ltt91G3njjDaJWq8mjjz5Kurq6IlqifOlLXyLbt28nHo+HqNVqcuTIkUU7dyG5rY8++ogAIA0NDaSpqYk0NTWR3bt3E0II+bd/+zdSX19PGhoayLXXXsuuH3+43W5y9913k40bN5Kenh6f3z311FPk05/+NJmZmYn2Yyz70iLczxlFFN/85jfJww8/TAgh5OGHHybf+ta3Fmzj9XrJLbfcQu6//36f199++21SUVFBKisrycMPP0wee+wxsm3bNp/97rvvPlJZWUkaGhpYQpTuW11dzfY9neF0Oslf//pX8qUvfYnk5+eTyy67jOzatYvodDpBJHHPPfeQe+65hyX9zgTw8yvFxcXk5ptv9vn97OwsmZiYWLBtBFh2Igj3c0YRRbRPlzfeeINUVFSQLVu2kNraWpKYmEguvvhin6fL7t27yVVXXUW8Xi/Zt28f2bhxIyFk/klTUVFBBgcHicPhII2NjayKcjrj7rvvJj/+8Y/Jhx9+SO69915SV1dHbrrpJvLiiy8SjUYTkCS2b99OvvSlL51RJEFBE5WTk5NEo9Gw10Wq5Cw7EYT7OaOIIlr84x//IFdccQX7/0MPPUQeeughn23uvPNO8sILL7D/02WOkH1PR9AnJIXH4yH/+Mc/yPbt20lDQwO5/vrryXPPPUdmZmaI2Wwm3/nOd8gXvvCFWDL/PpidnSWXXXYZqaqqIpdddtmCsiVFaWkpqa+vJ01NTWTDhg0R7x8txPqcJ7HsRBDu5wzNvESGYG5HQrYRsu/pCP9WcYlEgvPOOw+PPfYY2tra8MADD6CzsxNXXHEFNm7ciL6+Pvzxj3+M2s7fH5F0Ay9HN7FYn/N0QZwoEHdEihQSiQQtLS34xS9+gSNHjuDhhx/Gc889J+rNI0RxuZj7x+GLs35SGCCeI9LevXvx9a9/HV6vF4WFhQts3p9//nnWh5Camoonn3ySdSEGM+RZ6ZBIJNi8ebPoxxVqrchxHK644ooF3cSRWjPGEQZh1iZnBVwuFykvLycnTpxgCcmuri6fbd566y2fZGZra6vPvv39/aS8vJysWbOGHD16NGBS85NPPmFr5T179rCEKCHza21+kuxsgBiKS5pLmZmZIY2NjeT//u//CCGRS7uXGcuegwj3E48o4OtqRMfLUUckYH683KZNm7Bnzx5UVVUhOTkZ//M//+Oz76c//WnMzs7i+9//PtatW8ds3vnzIPju0ueeey7Gx8eX9oOuMIhhrRhsWHAk1oxxCEAYJolDIHbu3Eluv/129v8//elP5J577gm6/aOPPuqzfVlZGVm/fj1pbm4mTz311KKe6+kAIZoYs9lMjEYj+/d5551H3n77bcH7ryAse8QQ7idOFCLh5ZdfXkAU9957b8Bt33//fbJ27Vqi1WrZa8FC6LMVYnQTB9p/hWLZiSDcT3zpIRKEJEQBoKOjA3fccQfefvttZGVlsdeDhdBnK7KysvC3v/1twesFBQXMJ6SiogLt7e0R7R9HlAjDJGcVaBgbDYQkREdGRkhlZSX55JNPfF73D6HXrl1LCgsLg0rCP/jgA5Kens6UpT/+8Y/Z784kOflZhGWPGML9xIniJOx2O7nxxhtJc3MzeeSRR8jQ0FDEx9i9ezdZvXo1qaioYGHwk08+SZ588klCCCG33347USqV7AanSkJ+CF1TU0NUKlVISfgHH3xAPvOZzyx4/9NFTi5ENdnb28v+Tk1NTSQtLY38+te/JoTM+4YUFBQsaPA6jbHsRBDu56wnCirFHRkZIQ8++CDhOI4kJSURjuNISUkJefbZZ5f0fIRIwoMRxekiJ//Wt77lk2gM5ETGh9vtJnl5eWR4eJgQcmYYDPlh2Ykg3M9Zr8yUSqVwuVy48sor0dfXh9dffx2ffPIJDh8+jB07dvjkEZYCQiXh+/btQ1NTE66++mp0d3dHtO9yI1LV5Jk00Pp0RTyZCeCdd95BX18fOjs7fcxoxRhvFykICS8Jb25uxsjICFJTU7Fnzx5s3boV/f39gvZdCTibB1qfrjjrIwpgfkrZqlWrsG/fvuU+FUHVk/T0dDYFfdOmTXC5XNBqtYIrL0sBsZzI6EDrz33uc+y1M92JbEUizNrkjEdbWxv50Y9+ROrr60lKSgq56aab2Fp4OSCkejI1NcV8EPbv30+Ki4uJ1+tdsG95eTkpKysLWgH593//d5YQrKurIxKJhOkNgrVvi4FITJBfe+01cvnllwf9/RniRLbsOYhwP8t+Asv2weeXXTcBmATwBoBvA/gOgP8DcPUyn9smAMcBDAL4/snX7gJw18l/3wugG0A7gH8COD/IvrMAKgAknNy2NsR7Xgvgfd7/hwFkL9LnexTAAyf//QCAfw+x7UsAvuj32irev7cDeGm5r6cz/SecC/cZC47jWgE8BuCfhJBvLff5iA2O484DsIMQcuXJ/38XAAghDwfZ/gUAHxBCfn/y/8MAWggh2kU4tywALwMoATAK4HOEEB3HcQUA/kAI2XRyu2QAYwAqCCFzvP2fA7AO8y7xwwC+TAiZQhyLhrM5mekC0ADgTQDgOE4CAIQQ73KelIgoxPxNRjEO4JxAG568Ia/CfKRCQQC8y3EcAfAUIeS/xToxQsgsgEsDvD6J+YiI/t8KYEHZiRByek03PgNw1hIFIaSN47jnAWzjOM4M4GlCiH25z0tEBCp3BAsfrwXwCSFEx3vtAkLIJMdxuQDe4ziulxDyd9HPMo7TAmd71eOHAHYC+C6AMY7jrlzm8xET4wCKef8vwnw+JhBuAvAi/4WTT3cQQtQAXgWwcRHOMY7TBGc1URBCZgkhjxBCijF/E50HnFqGnOY4CGA1x3HlHMclYJ4M3vDfiOO4DAAXAXid91oKx3Fp9N8ArgDQtSRnHceKxFm79AAAbl6NtAHAZswn1naf/NVpn+ElhLg5jrsXwDsApJhfWnVzHHfXyd//7uSm1wF4lxBi4e2eB+DVk2ItGYAXCCF7l+7s41hpOGurHgDAcZwcwK0ArgHwK0LIx8t8SnHEsSJxVhNFHHHEIQxnwlo8Jpwh+Yg44lhUxCOKOOKIIyziT9M44ogjLOJEEUcccYRFnCjiiCOOsIgTRRxxxBEWcaKII444wiJOFHHEEUdYxIkijjjiCIv/D6GlE/+RcA2dAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from FEToy.fe.gaussrule import gausspoint2d\n",
    "\n",
    "my2dgp=gausspoint2d(ngp=3) # generated in the tensor-product way, ngp means the number of gauss point in each direction\n",
    "my2dgp.creategausspoint() # generate gauss points\n",
    "my2dgp.print()\n",
    "my2dgp.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "431daf1f-f9f2-47c7-89b4-9e58ddfbd7d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error is 0.00188868 %\n",
      "1-th gauss point: xi=  -9.06180e-01, eta=  -9.06180e-01, weight=   5.61343e-02\n",
      "2-th gauss point: xi=  -9.06180e-01, eta=  -5.38469e-01, weight=   1.13400e-01\n",
      "3-th gauss point: xi=  -9.06180e-01, eta=   0.00000e+00, weight=   1.34785e-01\n",
      "4-th gauss point: xi=  -9.06180e-01, eta=   5.38469e-01, weight=   1.13400e-01\n",
      "5-th gauss point: xi=  -9.06180e-01, eta=   9.06180e-01, weight=   5.61343e-02\n",
      "6-th gauss point: xi=  -5.38469e-01, eta=  -9.06180e-01, weight=   1.13400e-01\n",
      "7-th gauss point: xi=  -5.38469e-01, eta=  -5.38469e-01, weight=   2.29085e-01\n",
      "8-th gauss point: xi=  -5.38469e-01, eta=   0.00000e+00, weight=   2.72287e-01\n",
      "9-th gauss point: xi=  -5.38469e-01, eta=   5.38469e-01, weight=   2.29085e-01\n",
      "10-th gauss point: xi=  -5.38469e-01, eta=   9.06180e-01, weight=   1.13400e-01\n",
      "11-th gauss point: xi=   0.00000e+00, eta=  -9.06180e-01, weight=   1.34785e-01\n",
      "12-th gauss point: xi=   0.00000e+00, eta=  -5.38469e-01, weight=   2.72287e-01\n",
      "13-th gauss point: xi=   0.00000e+00, eta=   0.00000e+00, weight=   3.23635e-01\n",
      "14-th gauss point: xi=   0.00000e+00, eta=   5.38469e-01, weight=   2.72287e-01\n",
      "15-th gauss point: xi=   0.00000e+00, eta=   9.06180e-01, weight=   1.34785e-01\n",
      "16-th gauss point: xi=   5.38469e-01, eta=  -9.06180e-01, weight=   1.13400e-01\n",
      "17-th gauss point: xi=   5.38469e-01, eta=  -5.38469e-01, weight=   2.29085e-01\n",
      "18-th gauss point: xi=   5.38469e-01, eta=   0.00000e+00, weight=   2.72287e-01\n",
      "19-th gauss point: xi=   5.38469e-01, eta=   5.38469e-01, weight=   2.29085e-01\n",
      "20-th gauss point: xi=   5.38469e-01, eta=   9.06180e-01, weight=   1.13400e-01\n",
      "21-th gauss point: xi=   9.06180e-01, eta=  -9.06180e-01, weight=   5.61343e-02\n",
      "22-th gauss point: xi=   9.06180e-01, eta=  -5.38469e-01, weight=   1.13400e-01\n",
      "23-th gauss point: xi=   9.06180e-01, eta=   0.00000e+00, weight=   1.34785e-01\n",
      "24-th gauss point: xi=   9.06180e-01, eta=   5.38469e-01, weight=   1.13400e-01\n",
      "25-th gauss point: xi=   9.06180e-01, eta=   9.06180e-01, weight=   5.61343e-02\n"
     ]
    }
   ],
   "source": [
    "#analyticalsol=4.0/15.0\n",
    "analyticalsol=0.431111\n",
    "\n",
    "ngp=5\n",
    "my2dgp.setnum(ngp) # set up the number of gauss point in each direction\n",
    "my2dgp.update()\n",
    "sum=0.0\n",
    "for i in range(my2dgp.ngp2):\n",
    "    # ngp2=ngp*ngp\n",
    "    x,y,w=my2dgp.gpcoords[i,1],my2dgp.gpcoords[i,2],my2dgp.gpcoords[i,0]\n",
    "    # 0-> weight, 1-> xi, 2-> eta\n",
    "    #sum+=(x**2)*(y**4)*w\n",
    "    sum+=np.sin(x*x*y*y)*w\n",
    "\n",
    "error=100.0*np.abs(sum-analyticalsol)/analyticalsol # in percentage\n",
    "\n",
    "print('error is %g %%'%(error))\n",
    "my2dgp.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0f7c038-5500-49c6-ba35-45f4af5c4b21",
   "metadata": {},
   "source": [
    "$(x,y)$ -> $\\xi,\\eta$ ($\\in\\{-1,1\\}$)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}