analyseScript/Joschka/Function testing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "end_time": "2023-07-25T12:34:53.272379700Z",
     "start_time": "2023-07-25T12:34:52.798126900Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "outputs": [],
   "source": [
    "def find_fwhm(f_x, x):\n",
    "    ind = 0\n",
    "    for i in range(0, len(f_x)):\n",
    "        if f_x[i] > np.max(f_x)/2 and ind ==0:\n",
    "            a = i\n",
    "            ind = 1\n",
    "\n",
    "        if f_x[i] < np.max(f_x)/2 and ind ==1:\n",
    "            b = i\n",
    "            break\n",
    "    return x[b] - x[a]\n",
    "\n",
    "\n",
    "def gaussian(x, x0, sigma, A):\n",
    "    return A * np.exp(-0.5 * (x-x0)**2 / sigma**2)\n",
    "\n",
    "def thermal(x, x0, amp, sigma):\n",
    "    order = 15\n",
    "    res = np.exp(-0.5 * (x-x0)**2 / sigma**2)\n",
    "    return amp/1.6 * polylog(2, res)\n",
    "\n",
    "def Thomas_Fermi_1d(x, x0, amp, sigma):\n",
    "    res = (1-(( x - x0 ) / sigma) **2) **(3/2)\n",
    "    return amp * np.where(res > 0, res, 0)\n",
    "\n",
    "def density_1d(x, x0_bec, x0_th, amp_bec, amp_th, sigma_bec, sigma_th):\n",
    "    return thermal(x, x0_th, amp_th, sigma_th) + Thomas_Fermi_1d(x, x0_bec, amp_bec, sigma_bec)\n",
    "\n",
    "\n",
    "def polylog(pow, x):\n",
    "    order = 500\n",
    "    sum = 0\n",
    "    for k in range(1,order):\n",
    "        sum  += x ** k /k **pow\n",
    "    return sum\n",
    "\n",
    "def ThomasFermi_2d(x, y=0.0, centerx=0.0, centery=0.0, amplitude=1.0, sigmax=1.0, sigmay=1.0):\n",
    "    res = (1- ((x-centerx)/sigmax)**2 - ((y-centery)/sigmay)**2)**(3 / 2)\n",
    "    return amplitude * np.where(res > 0, res, 0)\n",
    "    # return amplitude * 5 / 2 / np.pi / max(tiny, sigmax * sigmay) * np.where(res > 0, res, 0)\n",
    "\n",
    "def polylog2_2d(x, y=0.0, centerx=0.0, centery=0.0, amplitude=1.0, sigmax=1.0, sigmay=1.0):\n",
    "    ## Approximation of the polylog function with 2D gaussian as argument. -> discribes the thermal part of the cloud\n",
    "    return amplitude/1.643  * polylog(2, np.exp( -((x-centerx)**2/(2 * (sigmax)**2))-((y-centery)**2/( 2 * (sigmay)**2)) ))\n",
    "    # return amplitude / 2 / np.pi / 1.20206 / max(tiny, sigmax * sigmay) * polylog(2, np.exp( -((x-centerx)**2/(2 * (sigmax)**2))-((y-centery)**2/( 2 * (sigmay)**2)) ))\n",
    "\n",
    "\n",
    "def density_profile_BEC_2d(x, y=0.0, amp_bec=1.0, amp_th=1.0, x0_bec=0.0, y0_bec=0.0, x0_th=0.0, y0_th=0.0,\n",
    "                           sigmax_bec=1.0, sigmay_bec=1.0, sigmax_th=1.0, sigmay_th=1.0):\n",
    "    return ThomasFermi_2d(x=x, y=y, centerx=x0_bec, centery=y0_bec,\n",
    "                          amplitude=amp_bec, sigmax=sigmax_bec, sigmay=sigmay_bec\n",
    "                          ) + polylog2_2d(x=x, y=y, centerx=x0_th, centery=y0_th,\n",
    "                                          amplitude=amp_th, sigmax=sigmax_th, sigmay=sigmay_th)\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:13:41.713691200Z",
     "start_time": "2023-07-25T13:13:41.678934700Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Jianshun Gao\\AppData\\Local\\Temp\\ipykernel_13296\\3599946220.py:38: RuntimeWarning: invalid value encountered in power\n",
      "  res = (1- ((x-centerx)/sigmax)**2 - ((y-centery)/sigmay)**2)**(3 / 2)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FWHM TF = 1.2161216121612153\n",
      "FWHM Thermal = 1.8241824182418238\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABK+UlEQVR4nO3de3hU5bk3/u+aSWZyPpCQEwYCqCAiB7WkwXraplJrbW13LaW+xc1WWy10W2NbpVVSdVesVXS3mxa1Ir4/a0XdrT1oYVMq+iJRlEOrCCjHcMpAgGQmk2QmM/P8/liz5kBmklkza81aM/l+rivXDJN1eMaR5OZ+7ud+JCGEABEREZFBLEYPgIiIiEY2BiNERERkKAYjREREZCgGI0RERGQoBiNERERkKAYjREREZCgGI0RERGQoBiNERERkqByjB5CIQCCAo0ePori4GJIkGT0cIiIiSoAQAi6XC3V1dbBY4uc/MiIYOXr0KOrr640eBhERESXh0KFDOOuss+J+PyOCkeLiYgDymykpKTF4NERERJQIp9OJ+vr60O/xeDIiGFGmZkpKShiMEBERZZjhSixYwEpERESGYjBCREREhmIwQkRERIZiMEJERESGYjBCREREhmIwQkRERIZiMEJERESGYjBCREREhmIwQkRERIZSHYy89dZbuO6661BXVwdJkvDqq68Oe86GDRtw4YUXwm634+yzz8aqVauSGCoRERFlI9XBiNvtxvTp07F8+fKEjt+/fz+uvfZaXHnlldi+fTu+973v4ZZbbsHatWtVD5aIiIiyj+q9aa655hpcc801CR+/YsUKjB8/Ho899hgA4LzzzsPGjRvx+OOPY86cOWpvT0RERFlG95qRtrY2NDc3R702Z84ctLW1xT3H4/HA6XRGfRFR9jnuOIaNz92HrRv+YPRQiMhAuu/a29HRgerq6qjXqqur4XQ60dfXh/z8/EHnLF26FPfff7/eQyMiA/V5fHA89WV8xr8TgX0S3ha5uOTKLxg9LCIygClX0yxevBjd3d2hr0OHDhk9JCLS2Ma//xEX+HcCACySgLTxcQQCwuBREZERdA9Gampq4HA4ol5zOBwoKSmJmRUBALvdjpKSkqgvIsou/n+8BAA4WjIDAPAp3zZs++SggSMiIqPoHow0NTVh/fr1Ua+tW7cOTU1Net+aiEyqs8eD83q3AgAKr/ohHLaxyJX8OPT+awaPjIiMoDoY6enpwfbt27F9+3YA8tLd7du3o729HYA8xTJ//vzQ8bfddhv27duHH/7wh9i1axd+9atf4aWXXsKdd96pzTsgooyz7R/bMM5yHD5YUTr5crjOugwAYG1/2+CREZERVAcj77//PmbOnImZM2cCAFpaWjBz5kwsWbIEAHDs2LFQYAIA48ePx2uvvYZ169Zh+vTpeOyxx/Cb3/yGy3qJRrCuXW8BAI4VTQHsRaiYfCkAoKHvI/R5/UYOjYgMoHo1zRVXXAEh4heZxequesUVV2Dbtm1qb0VEWSqv8wMAwEC1/I+asnNnA68Dk6V2fHDIgQsn1hk5PCJKM1OupiGi7OXzB1DT+zEAoGj8RQAAqbQeXdYK5Ep+HN35jpHDIyIDMBghorTad8KF83AAAFB59iz5RUnC6ZLJAID+wx8YNDIiMgqDESJKq/Y9H6JI6ocHdlhGnxt63T96CgCgsGu3UUMjIoMwGCGitOpvl+vHHPkTAWu4bK3grAsAAFX9e4esSyOi7MNghIjSynpSrhdxl02Ker1yolzMeo5oxwlnf9rHRUTGYTBCRGlV6NoPAJAipmgAwFY9GT5YUSL1ov3gJ0YMjYgMwmCEiNJGCIEqj9zyvfisKdHfzLHBkXsWAMB58J/pHhoRGYjBCBGlzcmefjTgKACgsmHqoO87CxoAAP7OvekcFhEZjMEIEaXN0QOfIE8agBc5sFeOH/R9b5n8mq1rX7qHRkQGYjBCRGnjOrIDAODIGQNYrIO+b6k8BwBQ0ts+6HtElL0YjBBR2gROyIWpp/IbYn6/sE5eYVM1cDhdQyIiE2AwQkRpY+/aAwDoK50Q8/uVY+Wi1lpxAn297rSNi4iMxWCEiNKm2C1Pv4jyibG/X1ELl8iHRRI43r4rnUMjIgMxGCGitCnzyitp7NWxgxHJYsGxnDEAAOcRBiNEIwWDESJKD/8ARgdOAABKas+Je9hJ+1gA4foSIsp+DEaIKC36Og8iBwF4RC6q6sbGPa63qB4AYOk6mK6hEZHBGIwQUVqcOizvSXMEVSjJt8c9zlciByr2Hq6oIRopGIwQUVr0Hpe7qp7IrR3yOOsoufFZSf8R3cdERObAYISI0sLXKW+Q58qrG/K4/KoGAECFzwEEAnoPi4hMgMEIEaVFjlNe1ttXWD/kcaXVDfAJC3LhA3o60jE0IjIYgxEiSosCt1wDMlA6bsjjqsoKcVRUAAB8J/frPi4iMh6DESJKi9JgDYilfOhgpLLQjiMYDQDo6eCGeUQjAYMRItJffzeKAk4AgH107FbwCotFQmeOXOTaf4LBCNFIwGCEiPR3Wu4ZclIUo6y8YtjDu4NFrr5TB/QcFRGZBIMRItJfsIHZITEao4vj9xhR9BfKLeFznYd0HRYRmQODESLS3cDJAwCAw2I0RhcNH4x4g11Y891sfEY0EjAYISLdeU7KGY4OjEZJfs7wJwSLXIs8DsA/oOfQiMgEGIwQke58p+UeIy5bNSRJGvZ4e3kdPCIHFgQAJzuxEmU7BiNEpDspGFD0FgzdCl5RWZyHY8FeI+hmMEKU7RiMEJHubO6jAAB/0dCt4BUVhfZwMMLMCFHWYzBCRPryeZDv6ZSfl56V0CmjCm04ilHyH7pZxEqU7RiMEJG+nHJWpF/koqCsOqFTKotsocxIgNM0RFmPwQgR6Ss4zXJUVGBUAst6AaC80IajohIA4DvVrtvQiMgcGIwQkb6C0yzHRAVGFdoSOiXXaoHTVgUAEJymIcp6DEaISF/BYOKoqEBZQW7Cp/Xn1wAALK6jugyLiMyDwQgR6UsJRpB4ZgQABorlYtdcbxfgdesxMiIyCQYjRKSvUM1IJcoLEg9GCorL4RL58h9YxEqU1RiMEJGuAl1yK/hjYpSqaZqKIhuOieDyXifrRoiyGYMRItJXMKtxXBqNInsC+9IElRfY2IWVaIRgMEJE+ul3wuJ1AgB682sS2pdGUZqfiyPswko0IjAYISL9BIOILlGI/MJSVadGZ0YOaT0yIjIRBiNEpJ+IHiNq6kUAoKwgF8fAaRqikYDBCBHpJxiMHBEVqlbSAEBZgQ1Hgl1YOU1DlN0YjBCRfkINzypRrqLHCBDMjCirabqPAEJoPToiMgkGI0Skn2BG45ioQLnKaZqompEBN9B3WuvREZFJMBghIv2EgpFRqrqvAkBJXg48sOGkKA5ei23hibIVgxEi0o/zGADAgXKUqawZybFaUJKXA4cyVeM6pvXoiMgkGIwQkX6CAUSHGKV6mgaQi1g7RLn8B2ZGiLIWgxEi0ofHBXh7AAAOUa66gBUAygtyw8EIMyNEWYvBCBHpIzhF4xL56EWe6qW9AFBaYAtP0zAzQpS1GIwQkT4ipmgAJDVNU16Qiw6wZoQo2zEYISJ9BIMHhygDABTnJVEzkp8LR6hmhMEIUbZiMEJE+ghOqzgwCsX2HFgtiW+Sp5ALWJXMCKdpiLIVgxEi0oerAwDQIcpRkq8+KwLIXVhDBay9JwGfR6vREZGJMBghIn0EMxmOFIKR8gIbulAEL4Lns26EKCsxGCEifSgNz8QolObnJHWJ0oJcABJOWiqirklE2SWpYGT58uVoaGhAXl4eGhsbsXnz5iGPf+KJJzBp0iTk5+ejvr4ed955J/r7+5MaMBFliMhpmiSKVwGgNJhRcYB1I0TZTHUwsnr1arS0tKC1tRVbt27F9OnTMWfOHBw/fjzm8S+88ALuuecetLa2YufOnXjmmWewevVq/OhHP0p58ERkUoEA0CMHIw5RHgoq1FKCmGMBrqghymaqg5Fly5bh1ltvxYIFCzBlyhSsWLECBQUFWLlyZczjN23ahEsuuQTf+MY30NDQgKuvvhrz5s0bNptCRBmstxMI+CAg4QTKkq4ZKQlO7xz2lckvsGaEKCupCka8Xi+2bNmC5ubm8AUsFjQ3N6OtrS3mObNnz8aWLVtCwce+ffvw+uuv4/Of/3wKwyYiUwsu63XllMMPa9LTNMp53J+GKLupqirr7OyE3+9HdXV11OvV1dXYtWtXzHO+8Y1voLOzE5/5zGcghIDP58Ntt9025DSNx+OBxxNewud0OtUMk4iMFsxgnLZWAkDSBax5uVbYciw4FmAXVqJspvtqmg0bNuChhx7Cr371K2zduhW///3v8dprr+HBBx+Me87SpUtRWloa+qqvr9d7mESkpWDQ0CnJq2CSnaYB5OxIqPGZ80jKQyMi81H1z5XKykpYrVY4HI6o1x0OB2pqamKec9999+Gb3/wmbrnlFgDABRdcALfbjW9961v48Y9/DItlcDy0ePFitLS0hP7sdDoZkBBlkmCh6XHI0yvJFrACct2Io0fZubcDEAKQ1HdzJSLzUpUZsdlsuOiii7B+/frQa4FAAOvXr0dTU1PMc3p7ewcFHFarFQAghIh5jt1uR0lJSdQXEWWQ4BLcI/4yAKlnRpSgBn6v3ImViLKK6onclpYW3HTTTbj44osxa9YsPPHEE3C73ViwYAEAYP78+RgzZgyWLl0KALjuuuuwbNkyzJw5E42NjdizZw/uu+8+XHfddaGghIiyTLDHyKHgKpjUMiO5GEAO+m0VyPOelItYCyu1GCURmYTqYGTu3Lk4ceIElixZgo6ODsyYMQNr1qwJFbW2t7dHZULuvfdeSJKEe++9F0eOHMHo0aNx3XXX4ac//al274KIzCU4TXPQK2c1k11NAwDFefKPKbd9tByMuI4BtdNSHyMRmUZSJe6LFi3CokWLYn5vw4YN0TfIyUFraytaW1uTuRURZaJgAevRgAY1I8FApjtnNCqwi8t7ibIQ96YhIm0N9AN9pwDI/UFyrRLycpP/UaM0PjudE9yfhst7ibIOgxEi0lYwWAhY7XCiEKX5uZBSWP2iZEY6oWyWx8wIUbZhMEJE2goWr3oLagBIKdWLAOGVOB3KihpmRoiyDoMRItJWcFlvn300AKA4hXoRACgJFrAe9XOzPKJsxWCEiLQVDBZcNjkYSaV4FQhnRg75SuUXXJymIco2DEaISFvBaZTuHDkYUTIbyVKmefYPBIORvtPAQF9K1yQic2EwQkTaCgYjpyzyfjKpZkaUTfaO9tmBnPyoexBRdmAwQkTaChawHoccjKTSCh4IZ0ZcHh9ESW3UPYgoOzAYISJtBZfeOpRgRKPVNAEB+Atrou5BRNmBwQgRaUeI0BTKsYBc41GUYs2IPccCm1X+UeUtkLed4DQNUXZhMEJE2unvAnz9AMI79hbbUwtGJEkK7U/Tl1clv8jlvURZhcEIEWlHCRLyy3HaI/94KUoxGAHCUzU9tmAwwswIUVZhMEJE2lF6gBTXwdXvk5+mOE0DhJcHd3N/GqKsxGCEiLSjrHIpqYXLIwcjqdaMAEBxsAi2y1opv8ACVqKswmCEiLSjTNMU16BHyYzYU1tNA4SnejotSmakQy6WJaKswGCEiLSj7NhbVIu+AT8AbTIjyjVOKJvl+T1yJ1YiygoMRohIO8FgpD+/OvSSFgWsyjWcAxJQEMyOcKqGKGswGCEi7QQDBHdwkzx7jgW2nNR/zChFsD39PqC4Tn6RXViJsgaDESLSTjBA6AkGI1qspAHCmRGXxwcUB7uwcvdeoqzBYISItOH3Ae7jAICuXHnVixZTNEC4ZqSn3wco+9Ow8RlR1mAwQkTa6HEAIgBYcnAaJQDCS3JTpQQ1PZ7IaRoGI0TZgsEIEWlDqeEoqoHLE5CfapQZCdWMRE3TMBghyhYMRohIG6HuqzVy0ABtlvUCQFGwV4k8TRPMjHA1DVHWYDBCRNqI6L4abnimRwFrbfT9iCjjMRghIm049dmXJvI68tLeYDDiPgH4BzS5PhEZi8EIEWnDFdEKXvNpGvk6fQN++PLKAUsuAMHsCFGWYDBCRNpQgpGScGakSIN9aQCgMGK6x+0VnKohyjIMRohIG6FN8mrR45GnT7TKjNhyLLAHO7m6PANsfEaUZRiMEJE2XOFgRMmMlGgUjABnLO9l4zOirMJghIhS53UDHqf8PLJmRKPVNJHXit6fhsEIUTZgMEJEqVMyFLYiIK8ktLRX02AkmBlxRbaEZzBClBUYjBBR6iKmaIBgPxBoVzMCxOk1wsZnRFmBwQgRpS5iWS8AuPrlAtYSjfamAc7owsrVNERZhcEIEaUuYlnvgD+A/gFt96YBIgtYB8It4TlNQ5QVGIwQUeqc4cyIOzhFA0T3B0lVdAFrcGmvtwfod2p2DyIyBoMRIkpdaJom3PDMnmOBLUe7HzGhAlaPD7AVAvbS4L05VUOU6RiMEFHqlIAgYllvsYb1IsAZmREgYkUNi1iJMh2DESJKnRIQlGi/SZ4iqukZEJ6qYeMzoozHYISIUiPEGZmRYCt4DetFIq8XDkaUIlZmRogyHYMRIkpN7ynA75WfF9VEbJKnTzDiGjRNw5oRokzHYISIUqMUrxZUAjm2cCt4jadpigZN07DxGVG2YDBCRKkJ9RiRg4Nejx8AUGizanqb4simZ0BE4zPWjBBlOgYjRJSaM1rBu71ysFCg9TTNmZkRTtMQZQ0GI0SUGmd0MNLr1SczElnAGgiI6JbwAb+m9yKi9GIwQkSpOTMzEsxcFNj0KWAFgN4BP1BYBUgWQPgB9wlN70VE6cVghIhSc2bNiJIZsWubGcnLtcAiyc97PT7AmgMUVUePgYgyEoMRIkqN68xpGn0yI5IkodDGxmdE2YjBCBGlJl7NiMaZEfmacjDiDq7YYeMzouzAYISIkucfCNdr6FwzAgAFwQBHWbHDFTVE2YHBCBElr8cBQACWXKCgAkDkahrtg5GiUGaE0zRE2YTBCBElL2JPGljkHyfhPiPaT9MU2JTMCKdpiLIJgxEiSp7Sil3JUCCyA2saMiOcpiHKCgxGiCh5ocxIbeilUGZE46Zn8jXPnKbh/jRE2YDBCBElT5keCQYF/oBA/0AAQHjli5YGr6YJBiP9XcBAn+b3I6L0YDBCRMlTMiMl0T1GAH0yI0XBOpTQffJKgdyC4FhYxEqUqRiMEFHynNGZEWUljUUC7Dna/3gpOLPpmSRxRQ1RFmAwQkTJO6NmRKnlKLTlQJIkzW83qIAViFhRw2CEKFMlFYwsX74cDQ0NyMvLQ2NjIzZv3jzk8V1dXVi4cCFqa2tht9tx7rnn4vXXX09qwERkIoNawcuZET2W9UZeN7S0F4hYUcNghChTqa4wW716NVpaWrBixQo0NjbiiSeewJw5c7B7925UVVUNOt7r9eKzn/0sqqqq8Morr2DMmDE4ePAgysrKtBg/ERnF0wN4nPLzksGZET3EzoxwmoYo06n+ibFs2TLceuutWLBgAQBgxYoVeO2117By5Urcc889g45fuXIlTp06hU2bNiE3NxcA0NDQkNqoich4yhSNrQiwFwNIQ2ZEWdobmRlh4zOijKdqmsbr9WLLli1obm4OX8BiQXNzM9ra2mKe86c//QlNTU1YuHAhqqurMXXqVDz00EPw+/0xjwcAj8cDp9MZ9UVEJnPGFA0Q2WNEn8yIsvleVGaEjc+IMp6qYKSzsxN+vx/V1dVRr1dXV6OjI/YPgn379uGVV16B3+/H66+/jvvuuw+PPfYY/vM//zPufZYuXYrS0tLQV319vZphElE6hIKRiO6roX1p9MmMKNM/vbEKWNn4jChj6b6aJhAIoKqqCk899RQuuugizJ07Fz/+8Y+xYsWKuOcsXrwY3d3doa9Dhw7pPUwiUksJRkrqQi8pQUKBDg3PgHDTs55YNSOuDkAIXe5LRPpS9ROjsrISVqsVDocj6nWHw4GampqY59TW1iI3NxdWa/hfSueddx46Ojrg9Xphs9kGnWO322G329UMjYjSzTk4M+LWOTOiFLD2ev0QQsjLh5VpIr8H6D0FFFbocm8i0o+qzIjNZsNFF12E9evXh14LBAJYv349mpqaYp5zySWXYM+ePQgEAqHXPv74Y9TW1sYMRIgoQ4SmaSIyIzrXjCiFsb6AgMcX/JmSYwMKKqPHREQZRfU0TUtLC55++mk899xz2LlzJ26//Xa43e7Q6pr58+dj8eLFoeNvv/12nDp1CnfccQc+/vhjvPbaa3jooYewcOFC7d4FEaVfqOFZRGZE2bFXp9U0kUuGe6NW1LDXCFEmU/3Pl7lz5+LEiRNYsmQJOjo6MGPGDKxZsyZU1Nre3g6LJRzj1NfXY+3atbjzzjsxbdo0jBkzBnfccQfuvvtu7d4FEaWfspS2JH2ZEatFQl6uBf0DAbg9PowqDGZXS2oBxwcsYiXKUEn9xFi0aBEWLVoU83sbNmwY9FpTUxPeeeedZG5FRGYkROzMiM41I4BcN9I/4D2jiJXLe4kyGfemISL1ek8Bfq/8vChiaa/Oq2mAcNYlcofgcDDCzAhRJmIwQkTqKbUZBZVyAWlQODOiXzASXt4bY38atoQnykgMRohIvRjdV4HImhH9pmmUKaCYjc9YwEqUkRiMEJF6ziPyY+mYqJd7g9kKXYORIRufMRghykQMRohIPefglTRAeG+aQh1rRpRlw1FLe5VxuE8APq9u9yYifTAYISL1uoOZkTOCkbRkRmwxMiMFFYBF3hUcPY4YZxGRmTEYISL1lGmakvA0jRAiTZmRGKtpItvCc6qGKOMwGCEi9WJM03h8AQSC+9TpWzMiX9sduZoGiFhRw+W9RJmGwQgRqRcKRsKZkcgaDr06sALhzIg7cpoGYOMzogzGYISI1Ol3Al6X/DwiM6IEB3m5Flgtkm63V2pG3N54wQgzI0SZhsEIEamjZEXyygBbYejl3jQ0PAPiND0D2PiMKIMxGCEidWIUrwLhTEWBTjv2KmI2PQPY+IwogzEYISJ14vQYUZb1pi8zcmYwwsZnRJmKwQgRqTNMwzM9V9IAcZqeRY7HeUzeVZiIMgaDESJSJ840TW8aeoxEXn/wappgZmTADXhcuo6BiLTFYISI1ImXGUlD91VgiNU0tkLAXio/51QNUUZhMEJE6sSrGQlN06QnM9I/EIDPH4j+JhufEWUkBiNEpE681TTpyoxErNbpHTijboSNz4gyEoMRIkqc1w30d8nP42RG9K4ZsVktyAk2VYvfhZWZEaJMwmCEiBKnNBSzFQN5JVHfcnvTkxmRJCl0j8Eratj4jCgTMRghosSFpmjqBn1LaUKmd58RIFyX0ntmF1bu3EuUkRiMEFHi4hSvAuEshd4dWCPv0RtvfxoWsBJlFAYjRJQ452H58YziVSB9e9MAiD9NUxocF4MRoozCYISIEjdEZiRdHVjle8TpNVJylvzY4wB8Xt3HQUTaYDBCRIkbappG2ZtG59U0QMRmeWdmRgorAasdgOCKGqIMwmCEiBIXp8cIYExmZNDOvZIUnqrpPqz7OIhIGwxGiChxCRSwpiMzEqoZObPpGQCUBqdquo/oPg4i0gaDESJKzEA/0HtSfh6rZsSTzsxIMBg5c2kvEK4b6T6k+ziISBsMRogoMUoNRk4+kF8e9S2fPwCPT94nJi2raexxCliBcGbEycwIUaZgMEJEiYmcopGkqG9FTpeko8+IUsDad2YBK8CaEaIMxGCEiBKTwEqaHIsEm1X/Hyv5oaW9rBkhygYMRogoMUoNRmn9oG8p0yX5NiukM7ImeghnRmJM04RqRpgZIcoUDEaIKDHKL3cl8xAh1GMkDfUigBz0AIA7VgGrMk3j6Qb6nWkZDxGlhsEIESVmiGAk1GMkDfUiQDjoGbQ3DQDYi4G8Uvk5i1iJMgKDESJKzFCZEW/6duwFIjfKi5EZAcJTSawbIcoIDEaIaHhCAF3BmpGysYO+HdqxNw09RuT7KJmROMGI0iGWvUaIMgKDESIaXn834HXJz2Pt2JvGfWmAyL1pYkzTABG79zIzQpQJGIwQ0fCUKZqCCsBWMOjb6dyXBogoYI07TcPlvUSZhMEIEQ1viHoRIGJfmjTVjCj38foC8PkDgw9gS3iijMJghIiGN0SPESBiX5o0raaJvM+Qm+VxmoYoIzAYIaLhmSwzYrNaYLXIzdVibpYXagl/RC6+JSJTYzBCRMMbJhhJd2ZEkqTwzr2xiliL6wBIgN8DuDvTMiYiSh6DESIanskyIwAigpEYmZEcG1BULT9n3QiR6TEYIaLhhWpGBvcYAdK/mgYIBz5KVmYQ1o0QZQwGI0Q0NP8A4DomP4+XGUlznxEgogtrrAJWIKJuhBvmEZkdgxEiGprrGCACgNUGFI6OeUjkrr3pUpAb7MIaq4AViGgJz2CEyOwYjBDR0JRf5iVjAEvsHxmG1IzYh+nCWsLMCFGmYDBCREMbpngViFhNk87MyFAFrEBEF1YGI0Rmx2CEiIY2TMMzICIzks6aEaWANe7+NAxGiDIFgxEiGtowmREhRCggKEzrahr5Xn3xMiNl4+THng5goD9NoyKiZDAYIaKhKcFIWezMiMcXCDU5LUhjZiQ/tLQ3TjBSMArILZSfMztCZGoMRohoaAl2XwWA/FwDMiMDcaZpJAkoC/ZF6TqYplERUTIYjBBRfEIAXUPXjCj1Ivm51tB+MemgLCOOmxkBgPLgVE1XexpGRETJYjBCRPH1dwNel/xcWSp7hlC9SJr2pVEoxbJxl/YCEZkRBiNEZsZghIjiU6Y3CkcDtoKYhyiZiYI09hiR7zfM0l6AwQhRhmAwQkTxnQ4GI8rKlBh6DdiXRr6fsrSXwQhRpksqGFm+fDkaGhqQl5eHxsZGbN68OaHzXnzxRUiShOuvvz6Z2xJRuim/xMtib5AHhDMj6ewxAkQu7eU0DVGmUx2MrF69Gi0tLWhtbcXWrVsxffp0zJkzB8ePHx/yvAMHDuD73/8+Lr300qQHS0RppkzTlJsvM5JQASt7jRBlBNXByLJly3DrrbdiwYIFmDJlClasWIGCggKsXLky7jl+vx833ngj7r//fkyYMCGlARNRGiUwTeM2YF8aIMEC1vxywFYkP2evESLTUhWMeL1ebNmyBc3NzeELWCxobm5GW1tb3PMeeOABVFVV4eabb07oPh6PB06nM+qLiAyQwDRNr7IvTZpX0yRUwMpeI0QZQVUw0tnZCb/fj+rq6qjXq6ur0dHREfOcjRs34plnnsHTTz+d8H2WLl2K0tLS0Fd9ffw9MYhIJ0JETNM0xD1MyYwYVcDq8QXgD4j4B7JuhMj0dF1N43K58M1vfhNPP/00KisrEz5v8eLF6O7uDn0dOnRIx1ESUUzuTmCgF4A05I69SmYk3dM0kcEPe40QZTZVPz0qKythtVrhcDiiXnc4HKipqRl0/N69e3HgwAFcd911odcCgYB845wc7N69GxMnThx0nt1uh91uVzM0ItKakhUprgVy4v99DGdG0huM2HMssEhAQMhTNcV5ubEPZDBCZHqqMiM2mw0XXXQR1q9fH3otEAhg/fr1aGpqGnT85MmT8cEHH2D79u2hry9+8Yu48sorsX37dk6/EJlZAitpgHBWIt0dWCVJCmVjIvfHGYQ1I0Smp/qfMi0tLbjppptw8cUXY9asWXjiiSfgdruxYMECAMD8+fMxZswYLF26FHl5eZg6dWrU+WVlZQAw6HUiMpkEVtIA4QLSdGdGALlo1uXxDdOFlfvTEJmd6p8ec+fOxYkTJ7BkyRJ0dHRgxowZWLNmTaiotb29HRYLG7sSZTwlkzDEShrAuMwIoARAnsRawvc4gIE+IDc/LWMjosQl9U+ZRYsWYdGiRTG/t2HDhiHPXbVqVTK3JKJ0O53YNI1Re9PI91SW9w7Xa6RY3vCv+zBQeU6aRkdEiWIKg4hiC/UYSbBmJM1LewH2GiHKFgxGiGiwQADoDi6pH2aaJpQZSfPeNEDEZnlDFbACXFFDZHIMRohoMNcxwO8FJCtQMmbIQ43MjCh1Kn0DQ2RGgHAwcpqZESIzYjBCRIMpGYTSswDr0BmPUJ8RAzIj+blKZmSYYGTUePnx9H6dR0REyWAwQkSDJbiSZsAfgNcnNzI0NDMyVAErAJQHg5FT+3QeERElg8EIEQ2W4EqayMJRI1bT5AcDIPdQBaxAODNy6oC85w4RmQqDESIa7PQB+bGsYcjDlHqRXKsEW076f5woHViHXNoLBFcESfLy3t6T+g+MiFRhMEJEgynTGRUThjxMqdXIz03/FA2Q4NJeAMjNA0rq5OenWDdCZDYMRohoMCUYGTV0MBLuvpr+KRogcmnvMMEIEH4vLGIlMh0GI0QUzdMDuI/Lz5XCzzjC3VeNyYyEl/YOM00DAOUN8iMzI0Smw2CEiKIpmYP8UUB+2ZCHGp0ZUaaHEsuMcEUNkVkxGCGiaAlO0QCRO/YalRlJsIAVCGd5OE1DZDoMRogomjKNkVAwonRfNapmJMECViAiM8JghMhsGIwQUbRQZmToehHA2H1pgHABa0LBiJIZcR+X62KIyDQYjBBRNFXTNMbtSwNEZkYSmKbJLwPyy+XnSh8VIjIFBiNEFE35RZ1AMBLal8bgaZr+gQD8gQQ6q7ItPJEpMRghorCBfqD7sPw8kcyIR1lNY2wBK5BgdoS9RohMicEIEYV1HQQgAFsxUFAx7OFGZ0bsORZYJPl5H4tYiTIWgxEiCgutpBkPSNKwh4f7jBiTGZEkKdyFVU0RKzMjRKbCYISIwlQUrwKRHViNyYzI91ZRxMrMCJEpMRghojCVwYjRq2kAtb1Ggu+r+xDg8+g4KiJSg8EIEYWp6DECGN9nBIjcLC+BzEhRNWArAkSAy3uJTITBCBGFKbUUw2yQp1AyI0a1gwfC9SoJZUYkCag4W37e+YmOoyIiNRiMEJHM5wVOH5SfK7+wh+E2eG8a+d4qurACQOU58uNJBiNEZsFghIhkp/cDwi8v6y2uSeiUUJ8RAwtYw5mRBKZpgIjMyB6dRkREajEYISJZ58fyY+XZCS3rFUKgd0CpGTE+M6LUrwxLCUaYGSEyDQYjRCRTaigqzkno8P6BAESwA7uhmRE1S3uBiGkaZkaIzILBCBHJlF/OlecmdLg74pd/fq6BmRF7kpmR3pNA7ymdRkVEajAYISJZ5DRNAno94eJVi2X4aR29qM6M2AqBkjHyc2ZHiEyBwQgRAUKonqZxh5b1GjdFE3n/hNrBKyomyo9c3ktkCgxGiEiesujvAiCFf1EPd4rB+9IoQqtpEml6pqjg8l4iM2EwQkThDEFZPZCbn9ApZtiXBgDyQ5kRFcGIUsTKzAiRKTAYIaJwvUiCUzSAOfalibx/wk3PgIjMyF4dRkREajEYIaLwdEVl4sGIGfalAZLowAqEi3RP7QMCKs4jIl0wGCGicDdSFcGIaTIjydSMlNYDVjvg9wBd7TqNjIgSxWCEiJKapgnvS2OOzIiq1TQWKzBqgvycy3uJDMdghGik83mB0wfk52oyIx7jd+wFktibRjE62NztxC6NR0REajEYIRrpTu0LbpBXBBTXJnxaKDNi8NJeJTMy4Bfw+gKJnzj6PPnxOIMRIqMxGCEa6U7slB9HT05ogzxFuGbE6GmacDCkKjtSNVl+VN4/ERmGwQjRSHc8+Mu46jxVp7kj2sEbKddqgS1H/lGmqm5EyYyc2I3Qjn9EZAgGI0QjXZLBiLKUttDgpb1ARK8RVV1YJwKWXMDbA3Qf0mlkRJQIBiNEI93xiGkaFXq95ihglceQxIoaa254B1/WjRAZisEI0Ujm88gFrABQNUXVqcovfqNrRoAke40ArBshMgkGI0QjWecn8kqavFKguEbVqaGlvQavpgGSzIwAXFFDZBIMRohGslC9yBRVK2mAiJoRM2VG1PYaYWaEyBQYjBCNZCeSqxcBwrvkmqpmxJNkZuTEbiCgokcJEWmKwQjRSJbkShoA6PWYcDWN2szIqAmA1QYM9ALd3KOGyCgMRohGsiSDEa8vAK9fziSYYZpG2TlYdWbEmhPej4d1I0SGYTBCNFJ5e8N70oxW22MknIEwQwFr0pkRgHUjRCbAYIRopDqxC4AACiqAotGqTlVWrdhyLMi1Gv9jJLyaJplgJLik2bFDwxERkRrG/xQhImM4PpQfq6eqPlVZ1ltoguJVILLPiMppGgComSY/dnyg4YiISA0GI0QjlfLLt+YC1af2KD1GTFAvAqSYGVHef+fHwECfhqMiokQxGCEaqZRgpHa66lPD+9KYLDOitukZIDd7K6gERAA4/pHGIyOiRDAYIRqJAgGgIzhNk0RmxK1M05hgWS8Q2WckicyIJIX/G3CqhsgQDEaIRqLT+wGvC7Daw0tbVTBT91UgPI6kMiMAgxEigzEYIRqJlF+61VPkXhsqhWtGzDFNoywvTqpmBGARK5HBkgpGli9fjoaGBuTl5aGxsRGbN2+Oe+zTTz+NSy+9FOXl5SgvL0dzc/OQxxNRGqRQvAqE+3mYZZpGyYz0pZwZ+ZBt4YkMoDoYWb16NVpaWtDa2oqtW7di+vTpmDNnDo4fPx7z+A0bNmDevHl444030NbWhvr6elx99dU4cuRIyoMnoiSFgpFpSZ3u9pirgFXJ0KjuwKqoOBvIyQMG3PIUFhGllepgZNmyZbj11luxYMECTJkyBStWrEBBQQFWrlwZ8/jf/va3+M53voMZM2Zg8uTJ+M1vfoNAIID169enPHgiSlLHP+XHJIORUGbELDUjwQxN34Af/oBQfwFrTrj5mfLfhojSRlUw4vV6sWXLFjQ3N4cvYLGgubkZbW1tCV2jt7cXAwMDGDVqVNxjPB4PnE5n1BcRaaTnBOA6BkCSa0aSuUQwA2GePiPhDE3fAItYiTKNqmCks7MTfr8f1dXVUa9XV1ejo6MjoWvcfffdqKuriwpozrR06VKUlpaGvurr69UMk4iG0vEP+XHUBMBenNQlwjUj5pimsedYYLVIAMLdYVWrDWaJjm7XZlBElLC0rqZ5+OGH8eKLL+IPf/gD8vLy4h63ePFidHd3h74OHTqUxlESZbkj2+THMRcmfYlwzYg5MiOSJIXrRpItYq0L/vc4uhUQSUz1EFHSVP0kqayshNVqhcPhiHrd4XCgpqZmyHMfffRRPPzww/jb3/6GadOGnqe22+2w2+1qhkZEiTqyRX6sSz4YUTIjZlnaC8j1K65+X3KNzwB5jx6rDeg7LRexjpqg7QCJKC5VmRGbzYaLLrooqvhUKUZtamqKe94jjzyCBx98EGvWrMHFF1+c/GiJKDVChIORMRclfZlQB1aT1IwA4V4jSTc+y7GF60aObNVoVESUCNXTNC0tLXj66afx3HPPYefOnbj99tvhdruxYMECAMD8+fOxePHi0PE/+9nPcN9992HlypVoaGhAR0cHOjo60NPTo927IKLEOI8A7uOAZA3XSCRBmQopMEnNCBAOjJJufAaEAzQGI0RppfqfNXPnzsWJEyewZMkSdHR0YMaMGVizZk2oqLW9vR0WSzjG+fWvfw2v14uvfvWrUddpbW3FT37yk9RGT0TqKFmR6ilAbn7Sl1GKRItMUjMChKeMepPtNQKEp66U/05ElBZJ/SRZtGgRFi1aFPN7GzZsiPrzgQMHkrkFEelB+Rd/ClM0QERmxEzTNLYUW8ID4f8ux/4B+H1JtconIvW4Nw3RSKJBvYgQImLXXvNM0xQEszRJL+0F5E6s9hLA1wec2KnRyIhoOAxGiEaKQCDcQyOFlTRefwC+YJdTM2VGClNd2gsAFgtQN0N+zqkaorRhMEI0UnR+DHhdQG4BMHpy0peJrMkoNNHSXiUw6k1lmgaIKGJlMEKULgxGiEaKQ+/Kj3UzU6qFUGoy7DkW5FjN8yNEmTJKerM8xZhg+4FD76U4IiJKlHl+khCRvtrfkR/Hxu8JlAizdV9VaJYZGftp+fHETqD3VIqjIqJEMBghGinag5tZphqMmLD7KqBRzQgAFFYCFefIz5VsEhHpisEI0Ujg6pBbnEMC6j+V0qWUmhEzdV8FNFpNoxgXDNjaE9uNnIhSw2CEaCRQpmiqpwJ5pSldym2yHXsV4Q6sKWZGgHD2SPnvRkS6YjBCNBKE6kU+nfKlekPBiNkyI8reNBpkRpT/Tke2AgN9qV+PiIbEYIRoJAjVi6QejPR4lO6r5syMpNQOXlE+HiiqAQID3KeGKA0YjBBlO08P0PGB/DzF4lUgXJNhupoRLdrBKyQpHLixboRIdwxGiLJdexsg/EDZWKB0TMqXU2oyzDZNo4ynV4uaESAcuB18W5vrEVFcDEaIst2+DfLj+Ms1uZySGSkwXQGrUjPihxAi9QuOv0x+PNgG+DypX4+I4mIwQpTt9r0pP064QpPLhVbTmG2aJpgZ8QcEPL5A6hesOg8orJI3zTu0OfXrEVFcDEaIspm7E3AE60U0yoy4TVrAWpBrhSTJz3u06DUiSeEATskuEZEuGIwQZbP9b8mP1VOBotGaXNLtMefSXotFQlEwW9PTr0EwAjAYIUoTBiNE2UzjehEgnHUozjNXMAIARcExaZIZAYAJwf9uR7cCfV3aXJOIBmEwQpTN9iv1ItoHI0Umy4wA4TG5tMqMlJ4l71MjAsCBjdpck4gGYTBClK069wCnDwCWXGDcbM0uO6IyI0B4qmbv37W7JhFFYTBClK0+XiM/NnwGsBdrdlmlHqPInqvZNbWiZEZ6PAPaXfTsq+THT/4X0GLJMBENwmCEKFspwci5n9P0si6POTfKA8LZGs0KWAG53iYnD+g+BDh2aHddIgphMEKUjfq6wm3Mz71as8t6fH54gz08ik2cGXFpOU1jKwAmXCk///iv2l2XiEIYjBBlo71/BwI+oPJcYNQEzS7rjtiEzoyZEWXqSNPMCABMCmaXdjMYIdIDgxGibPTxWvnx3DmaXlb5JZ+fa0WO1Xw/PnQpYAXCU11HtgAuh7bXJiIGI0RZx+cNTydoXi8iF4YWmXAlDQAU23WoGQGA4hqgbqb8XKnFISLNMBghyjb73wT6u+V9VZSdZzWi/JIvNmGPESAcJGlaM6KYdK38uPNP2l+baIRjMEKUbXb8QX6c8iXAom1dR6jhmUkzI+GmZxou7VWc/2X5ce8bgPuk9tcnGsEYjBBlE58X2PUX+bnyy1NDZu6+CuhYMwIAlWcDtdMB4Wd2hEhjDEaIssm+DfIUTVENMPbTml/e1W/uYES3mhHF+V+RHz/8H32uTzRCMRghyiYfvCw/Tvmi5lM0QAZM0+iZGQHC2aYDGwFXhz73IBqBGIwQZYv+7vD0wbSv63IL0xewar1R3pnKxwFnzQIggA9e0eceRCMQgxGibPHh7wFfPzB6MjDmQl1uYfbMiNIV1uMLhDrFam7GPPlx2//HvWqINMJghChbbP+t/DjjRkCSdLlFuIDVfK3ggeiusG69pmqmfhXILQBO7AIObdbnHkQjDIMRomxwYjdw+D1AsgLT5up2m9COvSbNjORYLcjPlQMS3epG8krChaxb/68+9yAaYRiMEGWDzU/Lj+d+Diiu1u02yi94s9aMABGNz/SqGwGAC+fLjzt+L9fqEFFKGIwQZbq+LmD7C/Lzxm/peiuls2mhiYOR0PJevTIjAFA/Cxh9HjDQC2x7Xr/7EI0QDEaIMt2254EBN1A1BRh/ua636gl2NjVrnxEgcnmvDl1YFZIEfPo2+fk7KwC/joEP0QjAYIQok/l9wOan5OeNt+lWuKpwKkt7TVozAqRhea9i2teBgkqgux3Y+Ud970WU5RiMEGWyD/8H6DoI5I8Cpn1N99t198nZhtJ8c66mAcLBiK7TNACQmwfMulV+vumXXOZLlAIGI0SZKuAH3vq5/Hz2d4HcfF1v1z/gD/XuKC0wbzBSnCePzdmXhqmTT90C5OQDR7cBn6zT/35EWYrBCFGm+vD3wMlPgPzy8L/QdaRkRSwSUGQz7zRNWTBQUsarq8JKYNYt8vO/PwgEdGq0RpTlGIwQZSKfF9jwkPy8aRFgL9b9lsov95L8XFgs+tampEKZQuru86bnhpfcCdiKgY5/cjdfoiQxGCHKRJufBE7tA4qqgcZvp+WWzgyoFwHCmZGu3jRkRgCgsAJo+o78fP0DgM+TnvsSZREGI0SZxt0JvPmI/PyqJWnJigCZUbwKRGZG0hSMAHJ2qrAKOLUXaPvv9N2XKEswGCHKNP97L+BxArXTgenfSNttQ9M0eeYORsoKbADSmBkB5BbxV/+n/PzNnwNdh9J3b6IswGCEKJN8vBb4x+8ASMDnHwUs6fsrzMzIMKZ9DRg7G/D1AX+5k0t9iVRgMEKUKfpOA3/+nvy8aaHckjyNIgtYzawsX6kZSVMBq0KSgC88DljtwJ51wHu/Se/9iTIYgxGiTBAIAH+4DXAdBUZNBP7l3rQPIVMyI0oBq9vrx4A/zUttqyYDn31Afv6/9wKOj9J7f6IMxWCEKBNsfAz4eA2QkwfcsEr3BmexKE3EzB6MFEfUtKR9qgaQVzed3Qz4+oEX5wG9p9I/BqIMw2CEyOw+eAX4e7A48vOPArXTDBlGpmRGrBYJJcG9c9JaxKqQJODLTwFlY4HTB4CX5nO5L9EwGIwQmdmev8nTMwAw69vAhd80bCingzUYZSZuBa9QVtQYkhkB5N4j81YDtiLgwP8DXl4gN6ojopgYjBCZ1c4/A7+bBwQGgPO/AnzuYUOHc7JH/td9ZZHd0HEkItwS3sAAoHoKMPd5uaB192vAKwuAgX7jxkNkYgxGiMxGCODdJ+X0vt8LnHcd8OUVaV3GG8vJHvkXe0WRzdBxJEKZSjrtNigzoph4JTDvBTkg2fUX4LkvAD3HjR0TkQkxGCEyk34n8D+3AH/9ISACwMz/A3x1FZBjbDbC4/PD5ZELWCsLzZ8ZUbI3nT0mqNU4uxn4P68AeWXA4feAJy8H9m0welREpsJghMgMhAB2/gVY3gh8+ApgyQHmLAW++N+A1fgdck+55axIjkVCSb7x4xlOVYkcjBx3mSAYAYDxlwG3rAcqzpaXZ//fLwGvfZ8rbYiCGIwQGe3ARmDVtcDqG+VfVOXjgX97Td58TTLH7riRUzSSScY0lOriPACAw2miGo3Ks4FvvwVc/O/yn997GvjFTODt/wL6u40dG5HBkgpGli9fjoaGBuTl5aGxsRGbN28e8viXX34ZkydPRl5eHi644AK8/vrrSQ2WKGv0dQHvrwSevEwORA6+LdcVfKYF+E4bMPbTRo8wyongdMeoDJiiASIyI06TZEYUtkK5S+v8PwJV5wP9XcC6JcCy84G/3g0c2cI28jQiqc63rl69Gi0tLVixYgUaGxvxxBNPYM6cOdi9ezeqqqoGHb9p0ybMmzcPS5cuxRe+8AW88MILuP7667F161ZMnTpVkzdBZHoD/UDHB8Chd+T9ZdrbgIBcgwGrDbhwvhyIlI4xdpxxHOuSMwx1pXkGjyQx1SXBzIjLRJmRSBOuAG77f/I+Q2//AujcDby7Qv4qbwDO/iww/lJg3CVAYaXRoyXSnSSEujC8sbERn/rUp/Df/y1vkx0IBFBfX4/vfve7uOeeewYdP3fuXLjdbvzlL38JvfbpT38aM2bMwIoVKxK6p9PpRGlpKbq7u1FSUqJmuETp4fMCfacAdyfQexLoPgyc3g+c2gd0fgIc/ygcfCiqpsgFqtO+LvelMLFH1uzCrzbsxU1N43D/l8z/j4j2k7247OdvwJ5jwa4HP2fuqaVAANi7Xg5Mdv8VGOiN/n5RjbxMePR5QFk9UDJGDlqLqoG8UrmXiZnfH41oif7+VpUZ8Xq92LJlCxYvXhx6zWKxoLm5GW1tbTHPaWtrQ0tLS9Rrc+bMwauvvhr3Ph6PBx5POL3qdDrVDDNh77zwIHD6YNRrEoaPzWIeEyOmS/paMV6TYsaM0a/F/nEkYlwu0bGeef1Y73vwSxLEoJcTu36sQ5Iba8zrpXAtiwggR3iQG/AgJ+BBrvAGHz2w+fuQF3DHGn2UnpwyHC6Ygr3FF2N3ySU4ZR8DOACsOwrg6LDnD/NWdPXbd9sBAGeVF6T3xkmqLrUjxyLB4wvgtue3oK4s/e3z1akD7Hch9/zvYKJzMxpcWzDetQVV/fuBng75a+/fY57phxX9OcXotxbBZ7HDJ+XCL9ngs0Q+5kJAgpAs8iMsgCRFPAcELPJX8PV4P1FSpcv/ujoFY0Kn/wZmNfbz30ddwyRD7q0qGOns7ITf70d1dXXU69XV1di1a1fMczo6OmIe39HREfc+S5cuxf33369maEkp2/caJvt26n4fGhn8QsJpFOO0KIZDlOGgqMEBUY2DohofiQYc7q8EeiTgOAD4AbQbPGL1pteXGT2EhNhzrDh/TCn+cagLa3c4jB6OSuOCX19BIfpwrnQYkyyHMFE6ilrpJMZIJ1EndWIUXMiV/LDCj0JfFwp9XQaPmzLdrs6vZ0Ywki6LFy+OyqY4nU7U19drfp/uSV9FW9fhGN+JFQ0Pfk0kGI2fGV1LUqyIO8b1Y10sxj1jR++JvjZY7PeVxH2lBI8L3TfWaNSNJfI7yX4+se4hJIv8r87gl99ig8+SJ//Zmoe+3DJ4coohJGvUeYUApgS/tCal+V9tZ5Xn41MN5Wm9Zyoe/eo0/PkfR+HP+IJQeVqsK/gV+ueTEMgN9MPucyHP74Ld50JOwAur8MqPgQFYhfI4AAkBQAASApBEIJgVFJCEkF+LeA6R5t2Oz2COfESm/3+j3viqsYbdW1UwUllZCavVCocj+l8aDocDNTU1Mc+pqalRdTwA2O122O36V+033vB93e9BRMY4p7oYLVcb8688IlJH1dJem82Giy66COvXrw+9FggEsH79ejQ1NcU8p6mpKep4AFi3bl3c44mIiGhkUT1N09LSgptuugkXX3wxZs2ahSeeeAJutxsLFiwAAMyfPx9jxozB0qVLAQB33HEHLr/8cjz22GO49tpr8eKLL+L999/HU089pe07ISIiooykOhiZO3cuTpw4gSVLlqCjowMzZszAmjVrQkWq7e3tsERs6DV79my88MILuPfee/GjH/0I55xzDl599VX2GCEiIiIASfQZMQL7jBAREWWeRH9/c28aIiIiMhSDESIiIjIUgxEiIiIyFIMRIiIiMhSDESIiIjIUgxEiIiIyFIMRIiIiMhSDESIiIjIUgxEiIiIylOp28EZQmsQ6nU6DR0JERESJUn5vD9fsPSOCEZfLBQCor683eCRERESklsvlQmlpadzvZ8TeNIFAAEePHkVxcTEkSdLsuk6nE/X19Th06FDW7nmT7e+R7y/zZft75PvLfNn+HvV8f0IIuFwu1NXVRW2ie6aMyIxYLBacddZZul2/pKQkK/8Hi5Tt75HvL/Nl+3vk+8t82f4e9Xp/Q2VEFCxgJSIiIkMxGCEiIiJDjehgxG63o7W1FXa73eih6Cbb3yPfX+bL9vfI95f5sv09muH9ZUQBKxEREWWvEZ0ZISIiIuMxGCEiIiJDMRghIiIiQzEYISIiIkNlfTDy05/+FLNnz0ZBQQHKyspiHtPe3o5rr70WBQUFqKqqwg9+8AP4fL4hr3vq1CnceOONKCkpQVlZGW6++Wb09PTo8A4St2HDBkiSFPPrvffei3veFVdcMej42267LY0jV6ehoWHQeB9++OEhz+nv78fChQtRUVGBoqIi/Ou//iscDkeaRpy4AwcO4Oabb8b48eORn5+PiRMnorW1FV6vd8jzzP4ZLl++HA0NDcjLy0NjYyM2b9485PEvv/wyJk+ejLy8PFxwwQV4/fXX0zRSdZYuXYpPfepTKC4uRlVVFa6//nrs3r17yHNWrVo16LPKy8tL04jV+8lPfjJovJMnTx7ynEz5/IDYP08kScLChQtjHm/2z++tt97Cddddh7q6OkiShFdffTXq+0IILFmyBLW1tcjPz0dzczM++eSTYa+r9u+wWlkfjHi9Xtxwww24/fbbY37f7/fj2muvhdfrxaZNm/Dcc89h1apVWLJkyZDXvfHGG7Fjxw6sW7cOf/nLX/DWW2/hW9/6lh5vIWGzZ8/GsWPHor5uueUWjB8/HhdffPGQ5956661R5z3yyCNpGnVyHnjggajxfve73x3y+DvvvBN//vOf8fLLL+PNN9/E0aNH8ZWvfCVNo03crl27EAgE8OSTT2LHjh14/PHHsWLFCvzoRz8a9lyzfoarV69GS0sLWltbsXXrVkyfPh1z5szB8ePHYx6/adMmzJs3DzfffDO2bduG66+/Htdffz0+/PDDNI98eG+++SYWLlyId955B+vWrcPAwACuvvpquN3uIc8rKSmJ+qwOHjyYphEn5/zzz48a78aNG+Mem0mfHwC89957Ue9t3bp1AIAbbrgh7jlm/vzcbjemT5+O5cuXx/z+I488gl/84hdYsWIF3n33XRQWFmLOnDno7++Pe021f4eTIkaIZ599VpSWlg56/fXXXxcWi0V0dHSEXvv1r38tSkpKhMfjiXmtjz76SAAQ7733Xui1v/71r0KSJHHkyBHNx54sr9crRo8eLR544IEhj7v88svFHXfckZ5BaWDcuHHi8ccfT/j4rq4ukZubK15++eXQazt37hQARFtbmw4j1NYjjzwixo8fP+QxZv4MZ82aJRYuXBj6s9/vF3V1dWLp0qUxj//a174mrr322qjXGhsbxbe//W1dx6mF48ePCwDizTffjHtMvJ9FZtXa2iqmT5+e8PGZ/PkJIcQdd9whJk6cKAKBQMzvZ9LnB0D84Q9/CP05EAiImpoa8fOf/zz0WldXl7Db7eJ3v/td3Ouo/TucjKzPjAynra0NF1xwAaqrq0OvzZkzB06nEzt27Ih7TllZWVS2obm5GRaLBe+++67uY07Un/70J5w8eRILFiwY9tjf/va3qKysxNSpU7F48WL09vamYYTJe/jhh1FRUYGZM2fi5z//+ZDTalu2bMHAwACam5tDr02ePBljx45FW1tbOoabku7ubowaNWrY48z4GXq9XmzZsiXqv73FYkFzc3Pc//ZtbW1RxwPy38lM+awADPt59fT0YNy4caivr8eXvvSluD9rzOKTTz5BXV0dJkyYgBtvvBHt7e1xj83kz8/r9eL555/Hv//7vw+5KWumfX6K/fv3o6OjI+rzKS0tRWNjY9zPJ5m/w8nIiI3y9NTR0REViAAI/bmjoyPuOVVVVVGv5eTkYNSoUXHPMcIzzzyDOXPmDLvJ4De+8Q2MGzcOdXV1+Oc//4m7774bu3fvxu9///s0jVSd//iP/8CFF16IUaNGYdOmTVi8eDGOHTuGZcuWxTy+o6MDNpttUM1QdXW1qT6vWPbs2YNf/vKXePTRR4c8zqyfYWdnJ/x+f8y/Y7t27Yp5Try/k2b/rAKBAL73ve/hkksuwdSpU+MeN2nSJKxcuRLTpk1Dd3c3Hn30UcyePRs7duzQdUPQZDU2NmLVqlWYNGkSjh07hvvvvx+XXnopPvzwQxQXFw86PlM/PwB49dVX0dXVhX/7t3+Le0ymfX6RlM9AzeeTzN/hZGRkMHLPPffgZz/72ZDH7Ny5c9giq0yRzPs9fPgw1q5di5deemnY60fWulxwwQWora3FVVddhb1792LixInJD1wFNe+xpaUl9Nq0adNgs9nw7W9/G0uXLjVtu+ZkPsMjR47gc5/7HG644QbceuutQ55rhs9wpFu4cCE+/PDDIespAKCpqQlNTU2hP8+ePRvnnXcennzySTz44IN6D1O1a665JvR82rRpaGxsxLhx4/DSSy/h5ptvNnBk2nvmmWdwzTXXoK6uLu4xmfb5ZYqMDEbuuuuuISNXAJgwYUJC16qpqRlUFayssqipqYl7zpmFOz6fD6dOnYp7TiqSeb/PPvssKioq8MUvflH1/RobGwHI/ypP1y+yVD7TxsZG+Hw+HDhwAJMmTRr0/ZqaGni9XnR1dUVlRxwOhy6fVyxq39/Ro0dx5ZVXYvbs2XjqqadU38+IzzCWyspKWK3WQSuXhvpvX1NTo+p4M1i0aFGokF3tv45zc3Mxc+ZM7NmzR6fRaausrAznnntu3PFm4ucHAAcPHsTf/vY31dnETPr8lM/A4XCgtrY29LrD4cCMGTNinpPM3+GkaFZ9YnLDFbA6HI7Qa08++aQoKSkR/f39Ma+lFLC+//77odfWrl1rmgLWQCAgxo8fL+66666kzt+4caMAIP7xj39oPDJ9PP/888JisYhTp07F/L5SwPrKK6+EXtu1a5dpC1gPHz4szjnnHPH1r39d+Hy+pK5hps9w1qxZYtGiRaE/+/1+MWbMmCELWL/whS9EvdbU1GTKAshAICAWLlwo6urqxMcff5zUNXw+n5g0aZK48847NR6dPlwulygvLxf/9V//FfP7mfT5RWptbRU1NTViYGBA1Xlm/vwQp4D10UcfDb3W3d2dUAGrmr/DSY1VsyuZ1MGDB8W2bdvE/fffL4qKisS2bdvEtm3bhMvlEkLI/yNNnTpVXH311WL79u1izZo1YvTo0WLx4sWha7z77rti0qRJ4vDhw6HXPve5z4mZM2eKd999V2zcuFGcc845Yt68eWl/f7H87W9/EwDEzp07B33v8OHDYtKkSeLdd98VQgixZ88e8cADD4j3339f7N+/X/zxj38UEyZMEJdddlm6h52QTZs2iccff1xs375d7N27Vzz//PNi9OjRYv78+aFjznyPQghx2223ibFjx4q///3v4v333xdNTU2iqanJiLcwpMOHD4uzzz5bXHXVVeLw4cPi2LFjoa/IYzLpM3zxxReF3W4Xq1atEh999JH41re+JcrKykIr2L75zW+Ke+65J3T822+/LXJycsSjjz4qdu7cKVpbW0Vubq744IMPjHoLcd1+++2itLRUbNiwIeqz6u3tDR1z5vu7//77xdq1a8XevXvFli1bxNe//nWRl5cnduzYYcRbGNZdd90lNmzYIPbv3y/efvtt0dzcLCorK8Xx48eFEJn9+Sn8fr8YO3asuPvuuwd9L9M+P5fLFfo9B0AsW7ZMbNu2TRw8eFAIIcTDDz8sysrKxB//+Efxz3/+U3zpS18S48ePF319faFr/Mu//Iv45S9/GfrzcH+HtZD1wchNN90kAAz6euONN0LHHDhwQFxzzTUiPz9fVFZWirvuuisqOn7jjTcEALF///7QaydPnhTz5s0TRUVFoqSkRCxYsCAU4Bht3rx5Yvbs2TG/t3///qj3397eLi677DIxatQoYbfbxdlnny1+8IMfiO7u7jSOOHFbtmwRjY2NorS0VOTl5YnzzjtPPPTQQ1FZrDPfoxBC9PX1ie985zuivLxcFBQUiC9/+ctRv+DN4tlnn435/2tkEjMTP8Nf/vKXYuzYscJms4lZs2aJd955J/S9yy+/XNx0001Rx7/00kvi3HPPFTabTZx//vnitddeS/OIExPvs3r22WdDx5z5/r73ve+F/ltUV1eLz3/+82Lr1q3pH3yC5s6dK2pra4XNZhNjxowRc+fOFXv27Al9P5M/P8XatWsFALF79+5B38u0z0/5fXXml/IeAoGAuO+++0R1dbWw2+3iqquuGvS+x40bJ1pbW6NeG+rvsBYkIYTQbtKHiIiISJ0R32eEiIiIjMVghIiIiAzFYISIiIgMxWCEiIiIDMVghIiIiAzFYISIiIgMxWCEiIiIDMVghIiIiAzFYISIiIgMxWCEiIiIDMVghIiIiAzFYISIiIgM9f8DZ9HJBtOOo04AAAAASUVORK5CYII="
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(-10, 10, 10000)\n",
    "\n",
    "TF = ThomasFermi_2d(x, y=0, centerx=0, centery=0, amplitude=1, sigmax=1,sigmay=1)\n",
    "th = polylog2_2d(x, sigmax=1)\n",
    "\n",
    "print(f'FWHM TF = {find_fwhm(TF, x)}')\n",
    "print(f'FWHM Thermal = {find_fwhm(th, x)}')\n",
    "\n",
    "plt.plot(x, TF)\n",
    "plt.plot(x, th)\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:13:46.864959500Z",
     "start_time": "2023-07-25T13:13:46.226243500Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Jianshun Gao\\AppData\\Local\\Temp\\ipykernel_13296\\3599946220.py:23: RuntimeWarning: invalid value encountered in power\n",
      "  res = (1-(( x - x0 ) / sigma) **2) **(3/2)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FWHM TF = 1.2161216121612153\n",
      "FWHM Thermal = 1.8241824182418238\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABSVUlEQVR4nO3deXxU9b3/8ddMkpnsCRCyYdgEwRVQJA1eq15T0Vqr7W2l6K9aarW20Kpor9IquFTxWkRsLxW7qO1tbV2u2l61WkqLK24groCCrEICCGQlmWTm/P44c2Yhk2RmMnvez8cjj5nMnHPmOwwkHz7fz/fztRmGYSAiIiKSJPZkD0BEREQGNwUjIiIiklQKRkRERCSpFIyIiIhIUikYERERkaRSMCIiIiJJpWBEREREkkrBiIiIiCRVdrIHEA6Px8OuXbsoKirCZrMlezgiIiISBsMwaGlpobq6Gru99/xHWgQju3btoqamJtnDEBERkSjs2LGDI444otfn0yIYKSoqAsw3U1xcnOTRiIiISDiam5upqanx/R7vTVoEI9bUTHFxsYIRERGRNNNfiYUKWEVERCSpFIyIiIhIUikYERERkaRKi5oRERGJPcMw6O7uxu12J3sokqaysrLIzs4ecNsNBSMiIoOQy+Vi9+7dtLe3J3sokuby8/OpqqrC4XBEfQ0FIyIig4zH42HLli1kZWVRXV2Nw+FQQ0mJmGEYuFwu9u7dy5YtWxg/fnyfjc36omBERGSQcblceDweampqyM/PT/ZwJI3l5eWRk5PDtm3bcLlc5ObmRnUdFbCKiAxS0f4vViRQLP4e6W+iiIiIJJWCEREREUkqBSMiIpIRVq1ahc1m4+DBg8keSkb51re+xQUXXBDX11AwIiIiKc9ms/X5dfPNNyd7iDFlBVaHf914440JH8u9997LQw89FNfX0GoaEUme9v2w5iGomgTjzkz2aCSF7d6923f/kUceYcGCBWzcuNH3WGFhIW+99VYyhhZXGzduDNogtrCwMKrruFyuqPuAlJSURHVeJJQZEZHkMAx4eCasvAX+8B+wbXWyRzRoGYZBu6s7KV+GYYQ1xsrKSt9XSUkJNpst6LHAX9Jr1qxh6tSp5OfnM3369KCgBeC+++7jyCOPxOFwMGHCBP7nf/4n6Hmbzcb999/Pl770JfLz8zn66KNZvXo1mzZt4vTTT6egoIDp06ezefNm3zmbN2/m/PPPp6KigsLCQk4++WT+8Y9/BF33l7/8JePHjyc3N5eKigq+9rWv9fu+y8vLQ77PHTt2cOGFF1JaWsrQoUM5//zz2bp1q+88a2rl9ttvp7q6mgkTJrB161ZsNhuPPvoop556Knl5eZx88sl89NFHvPnmm0ydOpXCwkLOOecc9u7d2+Na8aTMiIgkx7ZXYOcb3m8MePkeGFWX1CENVoe63Byz4PmkvPaHt84g3xHbX0U/+clPuPvuuxk+fDhXXnkl3/72t3nllVcAePLJJ7nqqqtYunQp9fX1PP3008yePZsjjjiCM844w3eN2267jSVLlrBkyRKuv/56LrroIsaOHcv8+fMZOXIk3/72t5k7dy5/+9vfAGhtbeWLX/wit99+O06nk9///vecd955bNy4kZEjR/LWW2/xwx/+kP/5n/9h+vTp7N+/n5deeimq99fV1cWMGTOoq6vjpZdeIjs7m5/+9KecffbZvPvuu74MyMqVKykuLmbFihVB5y9cuJClS5f63sdFF11EUVER9957L/n5+Vx44YUsWLCA++67L6rxRUPBiIgkx7uPmrejTjEDk03/gI4myI1/Slgy2+23385pp50GwA033MC5555LR0cHubm5LF68mG9961t8//vfB2DevHm89tprLF68OCgYmT17NhdeeCEA119/PXV1ddx0003MmDEDgKuuuorZs2f7jp80aRKTJk3yfX/bbbfx5JNP8te//pW5c+eyfft2CgoK+NKXvkRRURGjRo1iypQp/b6XI444Iuj7bdu28be//Q2Px8NvfvMbX+fcBx98kNLSUlatWsVZZ50FQEFBAb/5zW98wYmVObnuuuuC3sesWbNYuXIlp5xyCgCXXXZZ3GtEDqdgRESS45NV5u0pV0PbXtj3EWz+Jxz7lWSOalDKy8niw1tnJO21Y+2EE07w3a+qqgJgz549jBw5kvXr13PFFVcEHX/KKadw77339nqNiooKAI4//vigxzo6Omhubqa4uJjW1lZuvvlmnnnmGXbv3k13dzeHDh1i+/btAHzhC19g1KhRjB07lrPPPpuzzz6br3zlK/12wH3ppZcoKiryfT9kyBDeeecdNm3aFPQ4QEdHR9DU0fHHHx+yTiSc97Znz54+xxVrCkZEJPH2b4GD28CeDaOmw7h6MxjZ+rKCkSSw2WwxnypJppycHN99K3Pg8XgGfI2+rnvdddexYsUKFi9ezLhx48jLy+NrX/saLpcLgKKiItauXcuqVav4+9//zoIFC7j55pt58803KS0t7XUcY8aM6fF8a2srJ510En/84x97HD98+HDf/YKCgqjfW6R/XgMVcQHriy++yHnnnUd1dTU2m42nnnqq33NWrVrFiSeeiNPpZNy4cQlP/4hIitnuLVYdMRWchVAzzfx+Z+athpDUcvTRR/vqRyyvvPIKxxxzzICu+8orr/Ctb32Lr3zlKxx//PFUVlYGFZQCZGdnU19fz1133cW7777L1q1b+ec//xnxa5144ol8/PHHlJeXM27cuKCvRKx8iYeIg5G2tjYmTZrEsmXLwjp+y5YtnHvuuZxxxhmsW7eOq6++mu985zs8/3xyiqVEJAXsfse8HXGSeXvEyeZt4/vg0pb2Ej8/+tGPeOihh7jvvvv4+OOPWbJkCU888QTXXXfdgK47fvx4nnjiCdatW8c777zDRRddFJRdePrpp/n5z3/OunXr2LZtG7///e/xeDxMmDAh4te6+OKLKSsr4/zzz+ell15iy5YtrFq1ih/+8Ifs3LlzQO8jWSLOy51zzjmcc845YR+/fPlyxowZw9133w2YUenLL7/MPffc4yugEZFBxgpGqrwFf8UjoLASWhvM57SqRuLkggsu4N5772Xx4sVcddVVjBkzhgcffJDTTz99QNddsmQJ3/72t5k+fTplZWVcf/31NDc3+54vLS3liSee4Oabb6ajo4Px48fzpz/9iWOPPTbi18rPz+fFF1/k+uuv56tf/SotLS2MGDGCM888M6gnSTqxGeEu8g51ss3Gk08+2ef6489//vOceOKJLF261PfYgw8+yNVXX01TU1PIczo7O+ns7PR939zcTE1NDU1NTWn7By0iXh4PLDoCutp474K/s3SdnYlVRVy77ybsH/8dvrgYpl2e7FFmtI6ODrZs2cKYMWOi3vJdxNLX36fm5mZKSkr6/f0d96ZnDQ0NvmpdS0VFBc3NzRw6dCjkOYsWLaKkpMT3VVNTE+9hikii7N8MXW0Y2Xn8v6f2s3LDHpb9azNrD1Wbz+/5MLnjE5GES8kOrPPnz6epqcn3tWPHjmQPSURixTtFszv3SJo6/HPqj+/0/q+pUcGIyGAT92CksrKSxsbGoMcaGxspLi4mLy8v5DlOp5Pi4uKgLxHJEPs+AuDN9koAfjFrCtUluazrtDIj681W8SIyaMQ9GKmrq2PlypVBj61YsYK6OhWoiQxK3mDkvc4KHNl2vnBMBeeeUMVmoxo3WdDZBE3puSJARKITcTDS2trKunXrWLduHWAu3V23bp2vy9z8+fO55JJLfMdfeeWVfPLJJ/znf/4nGzZs4Je//CWPPvoo11xzTWzegYikl30fA7DZqKZ2zFByc7I4ZVwZXWSzzebNjuzdkMQBDh4DWL8g4hOLv0cRByNvvfUWU6ZM8fXUnzdvHlOmTGHBggWAuc2zFZiA2T3umWeeYcWKFUyaNIm7776b3/zmN1rWKzIYedzw2SbADEbqjhwGwLQxQ8my29jYbU7d8Nnm3q4gMWB122xvV08XGTjr71FgF9dIRdxn5PTTT+8zCgrVXfX000/n7bffjvSlRCTTNO2A7g5cZLPTGM4JI0oByHdkM768kC37zH1ErIBF4iMrK4vS0lLf/iP5+fm+tuAi4TIMg/b2dvbs2UNpaSlZWdHvM5Q5mxGISOrzTtF84qnCg51jq/3F6ceNKGHLXiszomAk3iorzT/rRG+IJpmntLTU9/cpWgpGRCRxvMWrm40qRpTmMaTAv6Po8SNK+L+1mqZJFJvNRlVVFeXl5XR1dSV7OJKmcnJyBpQRsSgYEZHE8QUj1RxTHbxkf0JlET83vNM0TTugqwNy1B003rKysmLyy0RkIFKy6ZmIZKj9nwCw1VPJ+PLCoKeOHF7IZxTTbOQDBhzYkoQBikgyKBgRkcQ5sBWA7UY5o8sKgp4qK3RQlJvDFkN1IyKDjYIREUkMd5evmdl2o4IxhwUjNpuNscMLFYyIDEIKRkQkMZp2guGhw8hhLyWMHlbQ45AjywrYZgUj3iyKiGQ+BSMikhje4GKnMZwCp4OyQkePQ8YOL2CnUeY9flsCByciyaRgREQSI6BeZExZQcgmWzVD89lhlJvfHNze43kRyUwKRkQkMQ6amY7tRjkjh+WHPGREaR47PMPNb5p2gMeTqNGJSBIpGBGRxPBmRnYY5YwozQt5yIgheTQwlG7DDm4XtOxO4ABFJFkUjIhIYviCkeFUl4RuZlZelIvNns0uw9xAz8qmiEhmUzAiIonhqxmpoKqXzEiW3UZVaa7qRkQGGQUjIhJ/HU1w6ABgZUZCByNg1o3sNLx1I1pRIzIoKBgRkfjzBhWfGUW0kUd1ae97zowozWeHFYxomkZkUFAwIiLxF1C86sy2M7SgZ48Ry4ghef5gRJkRkUFBwYiIxF/TDsBseFZVkhuyx4ilotjpn6ZRzYjIoKBgRETiz7snzafGMKp7KV61lBcFFLA27zT3tBGRjKZgRETizxuM7DaGUdVH8SqYmZG9lNBJDhgeaP40ESMUkSRSMCIi8ecNRnYZw6godvZ5aHlRLgZ2Goyh3nMVjIhkOgUjIhJ/3uzGLmMYZYV9ByNlhQ5sNtjlGRZ0rohkLgUjIhJf3Z3Q2gjALqOM4UV9ByPZWXbKCp3swsqM7Ij3CEUkyRSMiEh8Ne8CoBMH+ynqNxgBKC9ysssoM7/RNI1IxlMwIiLxZRWvMgywhRWMVBTnstvQNI3IYKFgRETiyxtM7HSb0y791YyAlRlRAavIYKFgRETiy1vzsdsYhiPbTnFudr+nlBfn+qdpmnfGc3QikgIUjIhIfHkzG7soY3ihs8/uq5bhRU7/NM2hA+Bqi+cIRSTJFIyISHwF9BgJp14EYFiBgxbyabd5G6RpqkYkoykYEZH4iqDHiGWYdyO9PTZN1YgMBgpGRCS+osmMFJrByE6r8ZkyIyIZTcGIiMRPRxN0NgNmAWu4wcjQAvO47d1DzAealBkRyWQKRkQkfrwZjTZ7Ee3kMtyb8ehPaV4OdhsBvUYUjIhkMgUjIhI/3nqRvfbhAAwLs2bEbrcxtMChLqwig4SCERGJH6vHCGaGY0h+eJkRwAxGUBdWkcFAwYiIxI83o7HD2311SEFO2KcOK3Cy29eFdScYRsyHJyKpQcGIiMSPt/B0a5dZiDo0ksxIocNfM9LVbjY/E5GMpGBEROLHO73yqcfMcJRGEIwMK3DQiYP27FLvtXbFenQikiIUjIhI/LQ0ANDIUIqc2Tiyw/+RM9Tb+OxgtreItWV3zIcnIqlBwYiIxI83gGg0hlAaQb0I+FfefGaziliVGRHJVApGRCQ+OlvA1QqYwUgk9SLgbwnfgLfxmTIjIhlLwYiIxEezGTy4sgtpJ5chBZEFI9Y0zU63NxhRZkQkYykYEZH4aDGDhzaH2fAs0syIFYzs6CrxXk+ZEZFMpWBEROLDW7x6MNsMRiJZSQNmS3iALa5i84FmBSMimUrBiIjEh3daZb/dXNY7NMIC1mJvMLLbuyzYyrSISOZRMCIi8WGtpPEWoEZaM5Kbk0VeThYNhrdmpP0z6O6M6RBFJDUoGBGR+PAGI7vckXdftZTm53CQQjxZzqBrikhmUTAiIvHhrfHY3mXWfERaM+I/x0ZnXkXQNUUksygYEZH48BawbnGZq2GGRjhNA/4i1jZnufmAdu8VyUgKRkQk9jweaDWDkc2HCoHIduy1lOab57TkmCtyNE0jkpkUjIhI7LXtBU83hs1Oo1EKQGledDUjAPuzvPvTaJpGJCMpGBGR2PNmMNx5w3GTRWGEm+RZSrwBzF6blveKZDIFIyISe95gxJVv1noU52ZHdRkrM9JoLe9VZkQkIykYEZHY8zY8a/cWnloNzCI1xBuM+PanUWZEJCMpGBGR2POupGn17ktTEmUwYk3TWMuDaWkAwxj4+EQkpSgYEZHY82YwDmabhafRBiPWNM2WziLzAbfL7MQqIhklqmBk2bJljB49mtzcXGpra3njjTf6PH7p0qVMmDCBvLw8ampquOaaa+jo6IhqwCKSBryZkc/sw4Dop2msYGTfIaDAu7y3WVM1Ipkm4mDkkUceYd68eSxcuJC1a9cyadIkZsyYwZ49e0Ie//DDD3PDDTewcOFC1q9fz29/+1seeeQRfvzjHw948CKSoryFpnsxV8FEnRnxTtMcbHdhFFWZD6rXiEjGiTgYWbJkCZdffjmzZ8/mmGOOYfny5eTn5/PAAw+EPP7VV1/llFNO4aKLLmL06NGcddZZzJo1q99sioikMe80TYO3x8hAp2k8BnQXVJoPKjMiknEiCkZcLhdr1qyhvr7efwG7nfr6elavXh3ynOnTp7NmzRpf8PHJJ5/w7LPP8sUvfrHX1+ns7KS5uTnoS0TSRFcHHDoAwI7uUiD6pb25OVk4vf1JfPvTKDMiknEi+gmxb98+3G43FRUVQY9XVFSwYcOGkOdcdNFF7Nu3j3/7t3/DMAy6u7u58sor+5ymWbRoEbfcckskQxORVGEFC9m5NHbmAm2U5EeXGQEYku+gobmDVudwCkGZEZEMFPfVNKtWreKOO+7gl7/8JWvXruWJJ57gmWee4bbbbuv1nPnz59PU1OT72rFjR7yHKSKxYgUjRVU0dXQD0U/TABTnmf9naskpD76+iGSMiDIjZWVlZGVl0djYGPR4Y2MjlZWVIc+56aab+OY3v8l3vvMdAI4//nja2tq44oor+MlPfoLd3jMecjqdOJ3OSIYmIqnCChaKq2na22XezR1AMOI992CWuTJHXVhFMk9EmRGHw8FJJ53EypUrfY95PB5WrlxJXV1dyHPa29t7BBxZWVkAGGpeJJJ5rGChqJKmQ2YwMrDMiHd5r7VZnrqwimSciKvK5s2bx6WXXsrUqVOZNm0aS5cupa2tjdmzZwNwySWXMGLECBYtWgTAeeedx5IlS5gyZQq1tbVs2rSJm266ifPOO88XlIhIBvFmRozCKlo6YhCMeItf9xjezMihA9B1CHLyBjZOEUkZEQcjM2fOZO/evSxYsICGhgYmT57Mc8895ytq3b59e1Am5MYbb8Rms3HjjTfy6aefMnz4cM477zxuv/322L0LEUkd3mCkM78Cjzf5GW3TM4Ai7zTNvm4n5ORDV7v5GkPHDnioIpIaolpvN3fuXObOnRvyuVWrVgW/QHY2CxcuZOHChdG8lIikG+80jbUvjSPbTm5O9FlQXwFrpxuKKmH/J+ZrKBgRyRjam0ZEYsubGWnOGdi+NBargLX5UBcUVQe9hohkBgUjIhI7huELFA54V78MOBjxnt/c0QXFagkvkokUjIhI7Bw6AN3mJpif2cx9aaLtvmrxZ0a6zWka0PJekQyjYEREYse7Wy95QzngMn+8DDQzUuQNZpo7NE0jkqkUjIhI7Fg9QIqqzEwGsZumaeno1jSNSIZSMCIisWNlRoqrfA3PBrKsF/zTPGYBqzcY0f40IhlFwYiIxE5zwL40Mei+Cv5gptXVjafQWzPS0mAWy4pIRlAwIiKxEzBNE6tgxKoZMQxo8S4Xxt1pFsuKSEZQMCIisRMwTWO1gi8a4GoaZ3YWzmzzR1Vzlx3yrQ3zNFUjkikUjIhI7FgBQlG1WXCKv537QAT1GvGtqGkY8HVFJDUoGBGR2Gnx79jb2mkFIwPLjIC/iDV4RY0yIyKZQsGIiMSGuxta95j3i/2ZkUJnDIKRvMCW8Gp8JpJpFIyISGy0NgIG2LMhvyygZmTg0zTWNZo7ugOmaZQZEckUCkZEJDasGo7CSgybLS7TNM2HAvenUc2ISKZQMCIisWFlKoqraHe58XjbgMQkGAnswqrGZyIZR8GIiMRGs7941aoXybLbyMvJGvClfZvldQR0YVVLeJGMoWBERGLDt5KmmtZOs16k0JmNzWYb8KWL8wKnabw1I217wd014GuLSPIpGBGR2LCCkeIqs9CU2EzRmNcJyIzkDQW7tyhWdSMiGUHBiIjERrO/FXxrDJf1QmABazfY7ZqqEckwCkZEJDasLEVRlW8lTXEMlvUGXqfN1e19QMGISCZRMCIisdHi37HX6jFSGKNpGus6VsZFjc9EMouCEREZOFcbdDab94urAvaliVEw4p3uaem0ghE1PhPJJApGRGTgrAyFowicRTFtBR94HV9mRI3PRDKKghERGbiADfKAgO6rsakZsYKRQ11uut0eNT4TyTAKRkRk4AKW9QIB+9LEJjNSEJBhaet0azWNSIZRMCIiAxdQvArEdF8aAEe2HWe2+eOqpTOg8ZmmaUQygoIRERm45uBgJNY1I+APbFo7u/2raVyt0NEcs9cQkeRQMCIiA9cSOhiJVc0IHFbE6igAZ0nwa4tI2lIwIiIDd1gBq6/PSAwzI1avEd/yXjU+E8kYCkZEZOB8BaxmLUesa0YgxPJeNT4TyRgKRkRkYAwjoBW8lRmJRzBiTvm0qvGZSMZRMCIiA9O+H9wu835hJW6PQbvLDcS2ZqTo8JbwanwmkjEUjIjIwFiZifwyyHb4gwViXDPSoyW8Gp+JZAoFIyIyMFZmwmp41mkWrzqz7TiyY/cjpudmeSpgFckUCkZEZGCszESPZb2xy4pAQAGrN9jRNI1I5lAwIiID4ytePbz7auzqRczrBTQ9g4AC1gbwuGP6WiKSWApGRGRgWg7PjMR2XxqLr2bEmqYpGA42OxhuaNsb09cSkcRSMCIiA3N4zUgcWsGDf7M8X2YkKxsKK7xjUN2ISDpTMCIiA5OgmpGiw5ueBbymGp+JpDcFIyIyML3UjFhNymKl8PCakYDXVOMzkfSmYEREoufu8tdrJKhmJCgzUqzMiEgmUDAiItFrbQQMsOdA/jAA2jrNlS2xrhnxZUZc3RiGYT5YpOW9IplAwYiIRC9wTxq7+eOk3WVmLvKdWTF9qSLvtI9h4Gs3r2kakcygYEREondY8Sr4MyMFjthmRnJz7GTZbUBA3YimaUQygoIREYneYbv1ArRZmRFHbDMjNputZ68RX+MzBSMi6UzBiIhEz5oeKa72PdQep5qRwGv6u7B6g6COg9B1KOavJyKJoWBERKLXV2YkDsFI0eGb5eWWQE6+eV+794qkLQUjIhI9X82IPzPS5s1aFMR4mgZCbJZns2lFjUgGUDAiItELmRkxp2nyY1zACv7lvS2hurCqbkQkbSkYEZHoWQFAUM2INzMS46W9EKJmBAJW1GiaRiRdKRgRkeh0tkJns3nfmxnxeAzau7xLe+NYMxI6M6JpGpF0pWBERKJj/fJ3FIGzCICObjdWc9RY9xkJvGab9qcRySgKRkQkOtYUTUC9iDV9YrOZTcpizVqhY63YAdT4TCQDKBgRkej46kX83VfbA7qv2my2mL9kobcOxeryCqjxmUgGiCoYWbZsGaNHjyY3N5fa2lreeOONPo8/ePAgc+bMoaqqCqfTyVFHHcWzzz4b1YBFJEX4MiMBreBd8SteBf8KneBpGm9mpqUB3xyRiKSViCd1H3nkEebNm8fy5cupra1l6dKlzJgxg40bN1JeXt7jeJfLxRe+8AXKy8t5/PHHGTFiBNu2baO0tDQW4xeRZGnuOU1jbWAXj3oR8K+mCZqmsYIhdye074eCYXF5bRGJn4h/YixZsoTLL7+c2bNnA7B8+XKeeeYZHnjgAW644YYexz/wwAPs37+fV199lZwcc9fN0aNHD2zUIpJ8vsyIf1mvVTMS6x17LdYKnaBpmmwH5JdB+z5zTApGRNJORNM0LpeLNWvWUF9f77+A3U59fT2rV68Oec5f//pX6urqmDNnDhUVFRx33HHccccduN3ukMcDdHZ20tzcHPQlIikmRAGrVTMSj4Zn4O/qGjRNA2p8JpLmIgpG9u3bh9vtpqKiIujxiooKGhpCr/H/5JNPePzxx3G73Tz77LPcdNNN3H333fz0pz/t9XUWLVpESUmJ76umpiaSYYpIIoRoeOarGYlDK3gIzIwcFoyo8ZlIWov7ahqPx0N5eTm/+tWvOOmkk5g5cyY/+clPWL58ea/nzJ8/n6amJt/Xjh074j1MEYmEYYRsBe/vvhqnzIi1msZ1WGZVjc9E0lpEPzHKysrIysqisbEx6PHGxkYqKytDnlNVVUVOTg5ZWf7/KR199NE0NDTgcrlwOBw9znE6nTidzkiGJiKJ1L4f3C7zfmHPfWniVcAamBkxDMO/fFiNz0TSWkSZEYfDwUknncTKlSt9j3k8HlauXEldXV3Ic0455RQ2bdqEx+PxPfbRRx9RVVUVMhARkTRg/dLPLzMLSL3aElTA2u0xcLn9P1PU+EwkvUU8TTNv3jx+/etf87vf/Y7169fzve99j7a2Nt/qmksuuYT58+f7jv/e977H/v37ueqqq/joo4945plnuOOOO5gzZ07s3oWIJFZzz4ZnEP+lvfk5/iBHjc9EMkfEPzFmzpzJ3r17WbBgAQ0NDUyePJnnnnvOV9S6fft27HZ/jFNTU8Pzzz/PNddcwwknnMCIESO46qqruP7662P3LkQksZo/NW+LRwQ93BbnmpHsLDu5OXY6ujy0dXYztMCblfE1PlMwIpKOovqJMXfuXObOnRvyuVWrVvV4rK6ujtdeey2alxKRVGStWglYSQMBmZE4TdOAmXXp6HIdtj+Ndxxte6HbFTR1JCKpT3vTiEjkfMFIcGbE1/QsTtM00Evjs/xhYM/xDkIrakTSjYIREYlcL9M07XHuMwKQH6rxmc2m5b0iaUzBiIhErpdpGitbkR+nmhEI2J9Gjc9EMoaCERGJXC/TNFZmpDCeNSO+zfJ6a3ymIlaRdKNgREQi09EMrhbz/mFLe60AIb41I9qfRiTTKBgRkchYWZHcUnAUBD3lW9obz2DEYWVGepumUTAikm4UjIhIZHopXvV4DN/S3nh1YIU+NstT4zORtKVgREQi00vx6qEufw1HYRwLWP3TNIfVjBRrmkYkXSkYEZHI9LaSxjttYreBMzt+P1p6z4wETNMYRtxeX0RiT8GIiESm11bw/n1pfLvpxkGvNSNWMNLVBp3NcXt9EYk9BSMiEplee4zEd8deS8gOrACOfMgtMe+r8ZlIWlEwIiKR6W9fmjiupDGv38vSXgiYqlHjM5F0omBERCLT2zSNK7479lp6bXoG6jUikqYUjIhI+Fxt0HHQvH94ZsRqBR/HfWmgjwLWwDEpGBFJKwpGRCR8VkMxRxHkFgc95Wt4FvfMiBnstB9ewApQVGneqvGZSFpRMCIi4Wvead4elhUB/zRN3DMj3pqU1r5qRpQZEUkrCkZEJHy9FK+Cv4A1ng3PwJ956ejy4PYc1k/EGpcKWEXSioIREQlfL8WrELC0N96raQKWDvfcn0bBiEg6UjAiIuHrIzPirxmJ7zSNI8tOtt0W9Jo+JTXmbWsjdLviOg4RiR0FIyISvr6CEWuTvDhnRmw2W++Nz/KHQXYuYECLsiMi6ULBiIiEr49pmnZXYjIj4K9L6ZEZsdn8gVLTzriPQ0RiQ8GIiISvz2maxHRgBf+KnR41IwAlR5i3TZ/GfRwiEhsKRkQkPF0d0P6ZeT/kaprEZUZ6naYBKLaCkR1xH4eIxIaCEREJj1WDkZ0HeUN6PN3amZiaEein8ZmVGWlWZkQkXSgYEZHwBE7R2Gw9nk5oZqSvxmcl3noW1YyIpA0FIyISHqsGI8QUDQTUjMS56Rn0UcAKqhkRSUMKRkQkPFYreKuXx2F8mZFEFLB6sy9914woMyKSLhSMiEh4rF/uVuYhgMdj+NrBx3tvGuhn515rmqazCTqa4z4WERk4BSMiEp6D3tUpIYKR9i5/hiIR0zRW9sVqtBbEWQS5JeZ9FbGKpAUFIyISHiszUtpzmqbdm6HIsttwZsf/x0qfmRHwTyWpbkQkLSgYEZH+GYa/b0eImpG2gCkaW4iVNrFW6KsZ6SUYsTrEqteISFpQMCIi/etoAlereb+PHXsTUbwK/l4mITuwgnqNiKQZBSMi0j9riiZ/GDjyezxtBSP5CegxAoFNz0LUjIB6jYikGQUjItK/PlbSgD8oSFRmxFfA2m/NiIIRkXSgYERE+tdHvQj4p0sS0X3VfJ0+9qaBgJoRBSMi6UDBiIj0r7/MSAJ37IV+du0F/zRN8y7weBIyJhGJnoIREelfP5mRVl/NSGKCEasdfLvLjWEYPQ8oqgZs4O6E9n0JGZOIRE/BiIj0r9+aEWs1TWKmaaygx+0x6OwOkfnIdkBhhXlfUzUiKU/BiIj0r6nvfWmsPiOJ6L4KkJ/jD3p6L2LVHjUi6ULBiIj0zd0FLbvN+73WjCQ2M2K323x1I/0u71WvEZGUp2BERPrWshsMD2Q5oGB4yEN8HVgTlBmBcBqfaXmvSLpQMCIifbN+mRePAHvoHxltCc6MgH8Zcf8t4RWMiKQ6BSMi0rd+ilchcG+aJGRGeus1opoRkbShYERE+tbPsl4IqBlJ4DRNoa8lvApYRdKdghER6Zv1y7y092DEv5omcdM0VmaktbfMSOko87a1Abo6EjQqEYmGghER6dtBKzPSxzSN1fQsgdM0Bf1lRvKHgqPQvK/siEhKUzAiIn0Lo2akPcF700DgZnm9ZEZsNigdad4/uC1BoxKRaCgYEZHeGUZYNSNtCd6bBgI3y+slMwIBwcj2BIxIRKKlYEREetfRBK5W8761VPYwbo/Boa7EdmCFMDbLAwUjImlCwYiI9M6a3igYDo78kIdYgQj4A4REsAKf9t6maUDTNCJpQsGIiPTugPeXuLUyJQRrmiTLbsOZnbgfKQXKjIhkDAUjItI7K6Ng/VIPwb+SJgubzZaIUZmvp5oRkYyhYEREemf9Eh/Se2bE2qiuMIH1IhCwmqa3jfIgoNdII3QdSsCoRCQaCkZEpHcRTNMksl4EwugzApA3RL1GRNJAVMHIsmXLGD16NLm5udTW1vLGG2+Edd6f//xnbDYbF1xwQTQvKyKJFs40jSvxreADX6/XPiOgXiMiaSLiYOSRRx5h3rx5LFy4kLVr1zJp0iRmzJjBnj17+jxv69atXHfddZx66qlRD1ZEEsgwAqZpRvd6mBUMJDoz4lva21fNCKhuRCQNRByMLFmyhMsvv5zZs2dzzDHHsHz5cvLz83nggQd6PcftdnPxxRdzyy23MHbs2AENWEQSpG0fdLUDtrC6rya6ZsR6vfa+akZAwYhIGogoGHG5XKxZs4b6+nr/Bex26uvrWb16da/n3XrrrZSXl3PZZZdFP1IRSSxrWqOoCrKdvR7mz4wkNhjJ9xWwdmMYRu8HWvUuCkZEUlZEPz327duH2+2moqIi6PGKigo2bNgQ8pyXX36Z3/72t6xbty7s1+ns7KSzs9P3fXNzcyTDFJFYOLDVvO1jJQ34p0kSuS9N4OsZhtl4rddgSJkRkZQX19U0LS0tfPOb3+TXv/41ZWVlYZ+3aNEiSkpKfF81Nb3viSEicWL98u5jJQ34l9YmOjOSl5OF1dakzyJWBSMiKS+inx5lZWVkZWXR2NgY9HhjYyOVlZU9jt+8eTNbt27lvPPO8z3m8XjMF87OZuPGjRx55JE9zps/fz7z5s3zfd/c3KyARCTRrGmafjIjvh17E1zAarPZKHBk09rZ7R1DL1NJVjBi9RrJyUvYGEUkPBFlRhwOByeddBIrV670PebxeFi5ciV1dXU9jp84cSLvvfce69at8319+ctf5owzzmDdunW9BhhOp5Pi4uKgLxFJsAP9L+uFgB17E1zACv4VNa19rajJGwKOIvO+eo2IpKSIf3rMmzePSy+9lKlTpzJt2jSWLl1KW1sbs2fPBuCSSy5hxIgRLFq0iNzcXI477rig80tLSwF6PC4iKSbMaRorM5KfhGCkwJkNLZ19r6ixeo3s+cDM9pSNT9wARSQsEf/0mDlzJnv37mXBggU0NDQwefJknnvuOV9R6/bt27Hb1dhVJK15PNC0w7zfzzSNlZVI9DQN+ItYw+o1sucD1Y2IpKio/iszd+5c5s6dG/K5VatW9XnuQw89FM1LikgitewGtwtsWVBU3eeh7UkqYA18zT4LWME/1WStEBKRlKIUhoj0ZBWvlhwBWX0HGVZWItFNz8CfjWnra38agKFjzNv9W+I8IhGJhoIREekpjN16Lb7MSIL7jIC/aLa9v2maId5g5ICCEZFUpGBERHoKcyUNBDQ9S8I0TYGvC2s/0zS+zMhWs0uaiKQUBSMi0pOv++rofg+1pkgSvVEe+LMx/RewjgJs4GqB9s/iPzARiYiCERHpaf8n5u3Qvje2dHsMOrrMRobJqBkJe7O8nFwoHmHeV92ISMpRMCIiPYUZjLQHFI4mo2bEv5qmn8wI+KdqVDciknIUjIhIsM5WaNtj3rcKP3thZSSy7TYcWYn/ceLrM9LfahrwTzlZgZaIpAwFIyISzMoc5A2FvNI+D7UanuU7srBZu9YlUNh9RkDLe0VSmIIREQkW5hQNQHsS96UBKPRmRtrDyoxomkYkVSkYEZFgVuYgjGDEmh5JVjBiZUZalRkRSWsKRkQkmC8z0ne9CPgzEsnYlwb8NSMRZUba9ph1MSKSMhSMiEiwCKZprIxEMvalAX9GJqyakbxSsw4GNFUjkmIUjIhIsAimaaw27AVJWNYLAR1Yw1naC5qqEUlRCkZExK+rA5o/Ne+HVTOS3AJWq+vroS43bk8Ybd5VxCqSkhSMiIjfwW2AAY4iyB/W7+HtvqW9yZ2mATMg6ZcyIyIpScGIiPgFFq+G0TfElxlJUgGrM9tOlt0cZ1hTNcqMiKQkBSMi4hdBvQj4A4D8JE3T2Gw231RNRC3hlRkRSSkKRkTEL4KVNBDQZyRJmRHztcPcLA/876tpB3R3xnFUIhIJBSMi4hdBjxFIfgdW87UjyIwUVoCjEAyPsiMiKUTBiIj4WcFIPxvkWfwdWJOYGbF6jYTT+Mxmg2HjzPufbYrjqEQkEgpGRMTU3QkHt5v3y8aHdUpbklfTmK9tZUbCmKYB/3v77OM4jUhEIqVgRERM+7eA4TaX9RZWhHVKu281TfKCkUKnVTMSZuOzYd5gZJ8yIyKpQsGIiJisTEHZuLCW9UJqTNNEtFkewLAjzVtlRkRShoIRETHt8/5yHhbeFA2kVgFre7gt4a1pmn0KRkRShYIRETFZBZ1lR4V9ipUZyU+Bpb1t4SztBX8B66H90L4/TqMSkUgoGBER076PzNuycWEd3u320NHlAZJbM5LvjHCzPEcBFI8w72tFjUhKUDAiImAYEU/TtAfsBZOfzKW91mqacAtYwZ8d0VSNSEpQMCIi0P4ZdBwEbP4Cz/5O8daL5GTZcGYnv89Ie7gFrBDQa0TBiEgqUDAiIv4pmtIayMkL6xR/vUjypmggoANrJJkRFbGKpBQFIyIS1Uoaq0YjmfvSgD8YCrtmBPzv87PNcRiRiERKwYiIBPQYiSQYMadFkrVjryWijfIsVpHu/k/AE8F5IhIXCkZExN+NNIJgpN3X8CwNp2lKaiDLCe6AFvgikjQKRkTEXzMSyTSNrxV8cqdpfBvlRVLAas+CoWPN+1reK5J0CkZEBrtuFxzYat6PJDOSApvkma9vbZQXQWYEYLi3udveDTEekYhESsGIyGC3/xP/BnlFVWGf1tqZ/H1pwL9RXme3h263J/wThx9t3u5RMCKSbApGRAa7PR+at8MnhL1BHvgLRpOfGfG/ftgt4QHKJ5q3e9fHeEQiEikFIyKDnTVNUX50RKdZBaOFSc6MOLLt5GSZQVR7JEWsVmZk70azA62IJI2CEZHBbo83MxBhMGJ1PE12ZgSiLGIddiTYc8DVCk074jQyEQmHghGRwS7KYKQtRWpGILDXSASZkawcf1t41Y2IJJWCEZHBrKvDLGAF/7RFmFKlHbw5BjMgao10RY3qRkRSgoIRkcHss4/NlTS5JVBUGdGpVgFrYZKbnoG/C2xEm+WBVtSIpAgFIyKDmfVLuPyYiFbSgH+aJj/JTc/AX0QbURdWUGZEJEUoGBEZzKxfwsMnRnxqWwoVsPo3y4syM7J3I3gi6FEiIjGlYERkMIuyeBX8WYjUKGA1xxBRASuYLeGzHNDVDk3ao0YkWRSMiAxmAwlGOq0+I8nPjES1tBcgK9u/H4/qRkSSRsGIyGDlavfvSRPhShrwdzvNT6VgJNLMCKhuRCQFKBgRGaz2bgAMyB8GhcMjOrXL7cHVbdZYFKZEzUiUm+UBVBxr3jZ+EMMRiUgkFIyIDFaN75u3FcdFfGrgEtr8FKgZsaaK2iPZm8ZSeYJ52/BeDEckIpFQMCIyWFm/fCuPj/jUVu90iLkvTPJ/jFiraSJuegb+97/vI+g6FMNRiUi4kv9TRESSwwpGqiZFfGq71Qo+BXqMgH9FT8SraQAKK6BgOBge/w7GIpJQCkZEBiOPZ2CZEd++NMmvF4EB9BkBs9mb9WegqRqRpFAwIjIYHdhi7lab5fQvbY1AKrWChwFmRkDBiEiSKRgRGYysX7oVx5i9NiLUmkKt4MG/a29UmRGACgUjIsmkYERkMBrAFA34MxCpMk1TEO3eNBZfZuR9tYUXSQIFIyKDUcO75q21rDVCrd4MREEK9BgBf1AU8a69lmHjIDsXutrMKSwRSaiogpFly5YxevRocnNzqa2t5Y033uj12F//+teceuqpDBkyhCFDhlBfX9/n8SKSAL7MSHTBSFuKFrC6ApqxRSQr29y5GPyBmogkTMTByCOPPMK8efNYuHAha9euZdKkScyYMYM9e/aEPH7VqlXMmjWLf/3rX6xevZqamhrOOussPv300wEPXkSi0LoXWnYDNrNmJAq+pb0p0PAMgpcYq4hVJP1EHIwsWbKEyy+/nNmzZ3PMMcewfPly8vPzeeCBB0Ie/8c//pHvf//7TJ48mYkTJ/Kb3/wGj8fDypUrBzx4EYlCwzvm7dCx4CyK6hK+aZoUyYxkZ9lxZps/ztqi6cIKUOXNEu1aF5tBiUjYIgpGXC4Xa9asob6+3n8Bu536+npWr14d1jXa29vp6upi6NChvR7T2dlJc3Nz0JeIxMinb5u3I06M+hK+AtYUWU0DgTv3RpkZqfb+eexaC4YRo1GJSDgiCkb27duH2+2moqIi6PGKigoaGhrCusb1119PdXV1UEBzuEWLFlFSUuL7qqmpiWSYItKXT9eYtyNOivoSqdb0DAa4WR6Ye/RkOeDQARWxiiRYQlfT3Hnnnfz5z3/mySefJDc3t9fj5s+fT1NTk+9rx44dCRylSAYzDH8wUj2QzEhqraYB/1ii2iwPINvhrxv5dG2MRiUi4YgoGCkrKyMrK4vGxsagxxsbG6msrOzz3MWLF3PnnXfy97//nRNO6LuC3+l0UlxcHPQlIjHQ/Cm07QFblr9GIgqpmBnx9RqJNjMC/myRghGRhIooGHE4HJx00klBxadWMWpdXV2v5911113cdtttPPfcc0ydOjX60YrIwFhZkYpjICcv6stYNSP5KbKaBgJqRqJdTQMBwciaGIxIRMIV8X9r5s2bx6WXXsrUqVOZNm0aS5cupa2tjdmzZwNwySWXMGLECBYtWgTAf/3Xf7FgwQIefvhhRo8e7astKSwspLCwMIZvRUT6Zf2PfwD1IuBvu54qe9NAYM1IlNM04J+62v0OuLujapUvIpGL+F/azJkz2bt3LwsWLKChoYHJkyfz3HPP+Ypat2/fjt3uT7jcd999uFwuvva1rwVdZ+HChdx8880DG72IRCYGxasQ0PQslWpGrC6sA8mMDBsHzmLobIa966Nuly8ikYnqJ8ncuXOZO3duyOdWrVoV9P3WrVujeQkRiTWPx99DYwDFqxDYgTWFpmm8gVHrQDIjdjtUT4EtL5iBm4IRkYTQ3jQig8W+j8DVAjn5MHxi1JfxeAzau1Kr6Rn461faB1LACv7+KzvfGuCIRCRcCkZEBosdr5u31VMGVAtxqMvt6wmWStM0hQ6rgHUAmRGAEd4iewUjIgmjYERksNj+mnk7sveVb+GwpmjsNsjNSZ0fIfkD7cBqGfk583bvemjfP8BRiUg4UucniYjE13bvlg0DDUYCGp7ZbLaBjipmrNb0AypgBSgog2HjzftWNklE4krBiMhg0NLgbXFug5qTB3SpthRseAaBmZEBTtMAjPIGbNvD23NLRAZGwYjIYGBN0VQcB7klA7qUFYykUsMzgEJnjDIj4M8eWX9uIhJXCkZEBgNfvcjnBnwpq8NpKjU8A8j3Le2NYTDy6VroOjTw64lInxSMiAwGvnqRGAQj3mkQq+NpqhjwRnmBhoyGwkrwdGmfGpEEUDAikuk6W6HhPfP+AItXwT9Nk2qZkZhslGex2fyBm+pGROJOwYhIptu+Ggw3lI6EkhEDvpy1miY/hXqMQOBGeW4MqxHKQIyabt5ue2Xg1xKRPikYEcl0n6wyb8ecFpPLpexqGu+0kdtj0NntGfgFR59q3m5bDd2dA7+eiPRKwYhIpvvkBfN27OkxuZxVwFqQYjUjgZmamNSNlB8NBeXQfQh2vDHw64lIrxSMiGSytn3Q6K0XyfDMSJbdRl5OjOtGrADOyi6JSFwoGBHJZFteNG8rjoPC4TG5pLWaJpV27LX4ilhj0WsEFIyIJIiCEZFMFuN6EfD38Sh05sTsmrFiTdXEpAsrwFjvn9uutXDoYGyuKSI9KBgRyWRbrHqR2AUjLR1dABTlptY0DQSsqInFNA1AyRHmPjWGB7a+HJtrikgPCkZEMtW+TXBgK9hz/MtUY8CXGUnFYCRWm+UFsqZqNv8zdtcUkSAKRkQy1UfPmbej/w2cRTG7bGuH+Yu+KMUKWCHGm+VZxp1p3n78d4hF/xIR6UHBiEimsoKRo86O6WVbOlI3MxLTzfIsY06D7Fxo2gGNH8TuuiLio2BEJBMdOuhvY37UWTG9dIt3mqYoN3ULWFtjmRlx5MPYM8z7G/8Wu+uKiI+CEZFMtPmf4OmGsqNg6NiYXbaz243L29001famgTjVjABM8GaXPlIwIhIPCkZEMtFHz5u3R82I6WUDazFSMhiJR80I+Ke6Pl0DLY2xvbaIKBgRyTjdLv//4I86J6aXtpb15juyyLLbYnrtWIj50l5LUSVUTzHvW7U4IhIzCkZEMs2WF6CjCQorYOTnYnppq3g1FXuMgH+zvJh1YA004Vzz9sO/xP7aIoOcghGRTPPBk+bt0V8Ge2xbtvu7r6ZmMGJlRmKyUd7hjv2KefvJKmj7LPbXFxnEFIyIZJJuF2x42rxv/fKMIf+y3tRbSQNQ4FtNE4fMSNk4qJoEhhvWKzsiEksKRkQyyServFM0lTGfogFo7fS2gk/RzEh+PPqMBDr2q+bt+0/E5/oig5SCEZFM8t5j5u0xsZ+igYDuqylaM2JlRtpjvZrGYmWbtr4MLQ3xeQ2RQUjBiEim6GiC9X8175/wjbi8RHNHqteMxLGAFWDIKDhiGmDAe4/H5zVEBiEFIyKZ4v0noLsDhk+EESfG5SVSeZM8CKgZ6YhTMAIweZZ5u/b32qtGJEYUjIhkinV/NG8nXwy2+PQASeVN8gCK88zC2jaXm263Jz4vctzXICcf9m2EHW/E5zVEBhkFIyKZYO9G2Pkm2LJgUnymaMCfGUnFfWkAigMyNi3xyo7kFvsLWdf+Pj6vITLIKBgRyQRv/Nq8PepsKCyP28tYHVhTdZomO8vu25+m2TvWuDjxEvP2gyfMWh0RGRAFIyLp7tBBWPeweb/2iri+VEuKF7CCf6qm+VAc60ZqpsHwo6GrHd7+Q/xeR2SQUDAiku7e/gN0tUH5MTDmtLi+VKoXsAIUe6eQmg7FMTNis8HnrjTvv7Yc3HEMfEQGAQUjIunM3Q1v/Mq8X3tl3ApXLS0pXsAKUGJlRuI5TQPm8un8Mmjaro6sIgOkYEQknb3/v3BwG+QNhRMujPvLWdkG6xd+KirOMwOl5nhmRgBycmHa5eb9V3+hZb4iA6BgRCRdedzw4s/M+9N/ADl58X05j+HLNpTkp3AwkohpGsvJ34HsPNj1Nny8Iv6vJ5KhFIyIpKv3n4DPPoa8If7/ocdRS0e37z//qZ0ZSdA0DUBBGUz7jnn/n7eBJ069TUQynIIRkXTU7YJVd5j36+aCsyjuL2llGnJz7DizY7/vTawkZDVNoFOuAUcRNLzrb8cvIhFRMCKSjt64H/Z/AoUVUPvdhLykFYyU5jkS8nrRshqfJWSaBqBgGNR937y/8lbo7kzM64pkEAUjIummbR+8cJd5/8wFCcmKABw85AJSe4oGEjxNY6mbCwXlsH8zrP7vxL2uSIZQMCKSbv5+I3Q2Q9UkmHRRwl42HVbSQMDS3kRlRsBsEX/WT837L/wMDu5I3GuLZAAFIyLp5KPn4Z0/ATb44mKwJ+6fsBWMFKd4MJLQ1TSBTrgQRp0C3Yfg6Wu01FckAgpGRNLFoQPwf1eb9+vmmC3JE8hXM5LCy3ohoM9IvDbK643NBucugexc2LQC3vxNYl9fJI0pGBFJBx4PPPFdaNkFQ4+Ef78x4UNoak+vaZqEZ0YAyifCF2417//9Rmj8MPFjEElDCkZE0sHLd8PHz5v/6/76Q3FvcBZKutWMuLo9HHK5Ez+AaVfAuHro7oA/z4L2/Ykfg0iaUTAikureexz+6S2O/OJiqDohKcNIl2maQmc2jizzR9v+dlfiB2CzwVd+BaWj4MBWePQSLfcV6YeCEZFUtukf8KR3d9hp34UTv5m0oRxoT4+lvTabjWGFZi+Uz1qTFAQUDIOLHgFHIWx9CR6bbTaqE5GQFIyIpKr1/wd/mgWeLjj2q3D2nUkdzr5W85fp8EJnUscRDn8wksQAoPxomPkHyHLCxmfg8dnQdSh54xFJYQpGRFKNYcDr95vpfbcLjj4PvrI8oct4Q7GyDMPSIBgZWmCOcV+yMiOWI8+AbzwMWQ7Y8DT87jxo3ZPcMYmkIAUjIqmkoxn+9zvwt/8EwwNT/h987SHITm4A0OX2cMC7mqasMLXbwQOUFXgzI20pMDUyvh7+3/9CbinsfBPuPw0+eSHZoxJJKQpGRFKBYcD6p2FZLbz/ONizYcYi+PJ/Q1Z2skfHAe8vdbsNSvNTPxhJes3I4cZ8Hr6zEoaNM5dn//7L8Mx1Wmkj4qVgRCTZtr4MD50Lj1xs/qIaMga+9Yy5+ZrNluzRAbDX+0t9aIGTLHtqjKkv1lRSSmRGLGXj4LsvwtRvm9+/+Wv4+RR45V7oaEru2ESSLKpgZNmyZYwePZrc3Fxqa2t54403+jz+scceY+LEieTm5nL88cfz7LPPRjVYkYxx6CC89QDc/3kzENn2ilno+G/z4PurYeTnkj3CIFYhaDpM0QAM807T7EtmAWsojgL40j1wyV+h/FjoOAgrFsCSY+Fv18Ona9RGXgaliPO/jzzyCPPmzWP58uXU1taydOlSZsyYwcaNGykvL+9x/KuvvsqsWbNYtGgRX/rSl3j44Ye54IILWLt2Lccdd1xM3oRIyuvqgIb3YMdr5v4y21eDx9uuPMsBJ15iBiIlI5I7zl7sOmiuAqksyU3ySMJTVWI2hfv0QHuSR9KLsafBlS+Z+wy98nPYtxFeX25+DRljNk0bc6q5101BWbJHKxJ3NsOILAyvra3l5JNP5r//29wm2+PxUFNTww9+8ANuuOGGHsfPnDmTtrY2nn76ad9jn/vc55g8eTLLly8P6zWbm5spKSmhqamJ4uLiSIYrkhjdLji0H9r2Qftn0LQTDmyB/Z/Avo9hz4f+4MNSfoxZoHrCN8y+FCnsZ89vYNm/NvPNz43itgtS/z8RO/a3c+pd/8KRbWfDrWdjT+WpJY8HNq80A5ONf4OuwwKowkqoOBaGT4TSGigeASVHQMFwyCs1e5mkyHSeyOHC/f0dUWbE5XKxZs0a5s+f73vMbrdTX1/P6tWrQ56zevVq5s2bF/TYjBkzeOqpp3p9nc7OTjo7/YVnzc3NkQwzbK89fBsc2Bb0mI3+Y7OQx4SI6aK+VojHbCFjxuDHQv84MkJcLtyxHn79UO+750M2jB4Ph3f9UIdEN9aQ1xvAteyGh2yjkxxPJ9meTnIMl/e2E4f7ELmetlCjD9KaXcrO/GPYXDSVjcWnsN85AhqBFbuAXf2e389bias/vr4dgJqhiW9DH42qklxysmy4uj384E9vU16c6suRq8F5LTnHfp8jm99gdMsaxrSsobxjC7Q2mF+bV4Y8000WHdlFdGYV0GVz0m134LY56LY76Lbl0G134rFlYWDHsNkwsPnug9383mYLeN5Obz9NohXzv64xDr6MGL/fdDXyi9dRPXpCUl47omBk3759uN1uKioqgh6vqKhgw4YNIc9paGgIeXxDQ0Ovr7No0SJuueWWSIYWldJPnmFi9/q4v44MDm7DxgGKOGAU0WiUss2oZKtRwTajgg+N0ezsKINWG+wBcAPbkzziyE06ojTZQwhLdpadY6qKeWdnE8+8tzvZw4nQKO/XVyngEEfZdjLRvp0jbbuosn1GtW0/1bZ9DKEFh81NFm4Kug9S0H0wyeOWdLdh3zfSIxhJlPnz5wdlU5qbm6mpqYn56zRN+BqrD+4M8UyoKLnnY0aY0XnoqPuwx0JcK+T/JkIeF954Q76HEEeEfl/Rv64R8o+pn+vZej8u7LGE+2fax/WsZwybnW670/fltjvotufSbXfSZc/lUE4pndlFGLasoOvkA0d7vwYiFbLwNUPymTZmaLKHEbZ7Zk7mqXW7cHs8yR7KAJnTYge8X759gA2DHE8Hzu4Wct0tOLtbyPa4yDJcZHm6fPezPS7sRjdgYDMMbHiwee+Dx/uY93HDwGbEenPB2OZFYv9PQQXDljHlI5P22hEFI2VlZWRlZdHY2Bj0eGNjI5WVlSHPqaysjOh4AKfTidMZ/7Rq7devi/triEhyjB1eyLwvHJXsYYhIGCJa2utwODjppJNYudI/d+nxeFi5ciV1dXUhz6mrqws6HmDFihW9Hi8iIiKDS8TTNPPmzePSSy9l6tSpTJs2jaVLl9LW1sbs2bMBuOSSSxgxYgSLFi0C4KqrruK0007j7rvv5txzz+XPf/4zb731Fr/61a9i+05EREQkLUUcjMycOZO9e/eyYMECGhoamDx5Ms8995yvSHX79u3YAzb0mj59Og8//DA33ngjP/7xjxk/fjxPPfWUeoyIiIgIEEWfkWRQnxEREZH0E+7vb+1NIyIiIkmlYERERESSSsGIiIiIJJWCEREREUkqBSMiIiKSVApGREREJKkUjIiIiEhSKRgRERGRpFIwIiIiIkkVcTv4ZLCaxDY3Nyd5JCIiIhIu6/d2f83e0yIYaWlpAaCmpibJIxEREZFItbS0UFJS0uvzabE3jcfjYdeuXRQVFWGz2WJ23ebmZmpqatixY0fG7nmT6e9R7y/9Zfp71PtLf5n+HuP5/gzDoKWlherq6qBNdA+XFpkRu93OEUccEbfrFxcXZ+RfsECZ/h71/tJfpr9Hvb/0l+nvMV7vr6+MiEUFrCIiIpJUCkZEREQkqQZ1MOJ0Olm4cCFOpzPZQ4mbTH+Pen/pL9Pfo95f+sv095gK7y8tClhFREQkcw3qzIiIiIgkn4IRERERSSoFIyIiIpJUCkZEREQkqTI+GLn99tuZPn06+fn5lJaWhjxm+/btnHvuueTn51NeXs6PfvQjuru7+7zu/v37ufjiiykuLqa0tJTLLruM1tbWOLyD8K1atQqbzRby68033+z1vNNPP73H8VdeeWUCRx6Z0aNH9xjvnXfe2ec5HR0dzJkzh2HDhlFYWMh//Md/0NjYmKARh2/r1q1cdtlljBkzhry8PI488kgWLlyIy+Xq87xU/wyXLVvG6NGjyc3Npba2ljfeeKPP4x977DEmTpxIbm4uxx9/PM8++2yCRhqZRYsWcfLJJ1NUVER5eTkXXHABGzdu7POchx56qMdnlZubm6ARR+7mm2/uMd6JEyf2eU66fH4Q+ueJzWZjzpw5IY9P9c/vxRdf5LzzzqO6uhqbzcZTTz0V9LxhGCxYsICqqiry8vKor6/n448/7ve6kf4bjlTGByMul4uvf/3rfO973wv5vNvt5txzz8XlcvHqq6/yu9/9joceeogFCxb0ed2LL76YDz74gBUrVvD000/z4osvcsUVV8TjLYRt+vTp7N69O+jrO9/5DmPGjGHq1Kl9nnv55ZcHnXfXXXclaNTRufXWW4PG+4Mf/KDP46+55hr+7//+j8cee4wXXniBXbt28dWvfjVBow3fhg0b8Hg83H///XzwwQfcc889LF++nB//+Mf9npuqn+EjjzzCvHnzWLhwIWvXrmXSpEnMmDGDPXv2hDz+1VdfZdasWVx22WW8/fbbXHDBBVxwwQW8//77CR55/1544QXmzJnDa6+9xooVK+jq6uKss86ira2tz/OKi4uDPqtt27YlaMTROfbYY4PG+/LLL/d6bDp9fgBvvvlm0HtbsWIFAF//+td7PSeVP7+2tjYmTZrEsmXLQj5/11138fOf/5zly5fz+uuvU1BQwIwZM+jo6Oj1mpH+G46KMUg8+OCDRklJSY/Hn332WcNutxsNDQ2+x+677z6juLjY6OzsDHmtDz/80ACMN9980/fY3/72N8NmsxmffvppzMceLZfLZQwfPty49dZb+zzutNNOM6666qrEDCoGRo0aZdxzzz1hH3/w4EEjJyfHeOyxx3yPrV+/3gCM1atXx2GEsXXXXXcZY8aM6fOYVP4Mp02bZsyZM8f3vdvtNqqrq41FixaFPP7CCy80zj333KDHamtrje9+97txHWcs7NmzxwCMF154oddjevtZlKoWLlxoTJo0Kezj0/nzMwzDuOqqq4wjjzzS8Hg8IZ9Pp88PMJ588knf9x6Px6isrDR+9rOf+R47ePCg4XQ6jT/96U+9XifSf8PRyPjMSH9Wr17N8ccfT0VFhe+xGTNm0NzczAcffNDrOaWlpUHZhvr6eux2O6+//nrcxxyuv/71r3z22WfMnj2732P/+Mc/UlZWxnHHHcf8+fNpb29PwAijd+eddzJs2DCmTJnCz372sz6n1dasWUNXVxf19fW+xyZOnMjIkSNZvXp1IoY7IE1NTQwdOrTf41LxM3S5XKxZsyboz95ut1NfX9/rn/3q1auDjgfz32S6fFZAv59Xa2sro0aNoqamhvPPP7/XnzWp4uOPP6a6upqxY8dy8cUXs3379l6PTefPz+Vy8Yc//IFvf/vbfW7Kmm6fn2XLli00NDQEfT4lJSXU1tb2+vlE8284GmmxUV48NTQ0BAUigO/7hoaGXs8pLy8Peiw7O5uhQ4f2ek4y/Pa3v2XGjBn9bjJ40UUXMWrUKKqrq3n33Xe5/vrr2bhxI0888USCRhqZH/7wh5x44okMHTqUV199lfnz57N7926WLFkS8viGhgYcDkePmqGKioqU+rxC2bRpE7/4xS9YvHhxn8el6me4b98+3G53yH9jGzZsCHlOb/8mU/2z8ng8XH311Zxyyikcd9xxvR43YcIEHnjgAU444QSamppYvHgx06dP54MPPojrhqDRqq2t5aGHHmLChAns3r2bW265hVNPPZX333+foqKiHsen6+cH8NRTT3Hw4EG+9a1v9XpMun1+gazPIJLPJ5p/w9FIy2Dkhhtu4L/+67/6PGb9+vX9Flmli2je786dO3n++ed59NFH+71+YK3L8ccfT1VVFWeeeSabN2/myCOPjH7gEYjkPc6bN8/32AknnIDD4eC73/0uixYtStl2zdF8hp9++ilnn302X//617n88sv7PDcVPsPBbs6cObz//vt91lMA1NXVUVdX5/t++vTpHH300dx///3cdttt8R5mxM455xzf/RNOOIHa2lpGjRrFo48+ymWXXZbEkcXeb3/7W8455xyqq6t7PSbdPr90kZbByLXXXttn5AowduzYsK5VWVnZoyrYWmVRWVnZ6zmHF+50d3ezf//+Xs8ZiGje74MPPsiwYcP48pe/HPHr1dbWAub/yhP1i2wgn2ltbS3d3d1s3bqVCRMm9Hi+srISl8vFwYMHg7IjjY2Ncfm8Qon0/e3atYszzjiD6dOn86tf/Sri10vGZxhKWVkZWVlZPVYu9fVnX1lZGdHxqWDu3Lm+QvZI/3eck5PDlClT2LRpU5xGF1ulpaUcddRRvY43HT8/gG3btvGPf/wj4mxiOn1+1mfQ2NhIVVWV7/HGxkYmT54c8pxo/g1HJWbVJymuvwLWxsZG32P333+/UVxcbHR0dIS8llXA+tZbb/kee/7551OmgNXj8Rhjxowxrr322qjOf/nllw3AeOedd2I8svj4wx/+YNjtdmP//v0hn7cKWB9//HHfYxs2bEjZAtadO3ca48ePN77xjW8Y3d3dUV0jlT7DadOmGXPnzvV973a7jREjRvRZwPqlL30p6LG6urqULID0eDzGnDlzjOrqauOjjz6K6hrd3d3GhAkTjGuuuSbGo4uPlpYWY8iQIca9994b8vl0+vwCLVy40KisrDS6uroiOi+VPz96KWBdvHix77GmpqawClgj+Tcc1VhjdqUUtW3bNuPtt982brnlFqOwsNB4++23jbfffttoaWkxDMP8i3TccccZZ511lrFu3TrjueeeM4YPH27Mnz/fd43XX3/dmDBhgrFz507fY2effbYxZcoU4/XXXzdefvllY/z48casWbMS/v5C+cc//mEAxvr163s8t3PnTmPChAnG66+/bhiGYWzatMm49dZbjbfeesvYsmWL8Ze//MUYO3as8fnPfz7Rww7Lq6++atxzzz3GunXrjM2bNxt/+MMfjOHDhxuXXHKJ75jD36NhGMaVV15pjBw50vjnP/9pvPXWW0ZdXZ1RV1eXjLfQp507dxrjxo0zzjzzTGPnzp3G7t27fV+Bx6TTZ/jnP//ZcDqdxkMPPWR8+OGHxhVXXGGUlpb6VrB985vfNG644Qbf8a+88oqRnZ1tLF682Fi/fr2xcOFCIycnx3jvvfeS9RZ69b3vfc8oKSkxVq1aFfRZtbe3+445/P3dcsstxvPPP29s3rzZWLNmjfGNb3zDyM3NNT744INkvIV+XXvttcaqVauMLVu2GK+88opRX19vlJWVGXv27DEMI70/P4vb7TZGjhxpXH/99T2eS7fPr6Wlxfd7DjCWLFlivP3228a2bdsMwzCMO++80ygtLTX+8pe/GO+++65x/vnnG2PGjDEOHTrku8a///u/G7/4xS983/f3bzgWMj4YufTSSw2gx9e//vUv3zFbt241zjnnHCMvL88oKyszrr322qDo+F//+pcBGFu2bPE99tlnnxmzZs0yCgsLjeLiYmP27Nm+ACfZZs2aZUyfPj3kc1u2bAl6/9u3bzc+//nPG0OHDjWcTqcxbtw440c/+pHR1NSUwBGHb82aNUZtba1RUlJi5ObmGkcffbRxxx13BGWxDn+PhmEYhw4dMr7//e8bQ4YMMfLz842vfOUrQb/gU8WDDz4Y8u9rYBIzHT/DX/ziF8bIkSMNh8NhTJs2zXjttdd8z5122mnGpZdeGnT8o48+ahx11FGGw+Ewjj32WOOZZ55J8IjD09tn9eCDD/qOOfz9XX311b4/i4qKCuOLX/yisXbt2sQPPkwzZ840qqqqDIfDYYwYMcKYOXOmsWnTJt/z6fz5WZ5//nkDMDZu3NjjuXT7/KzfV4d/We/B4/EYN910k1FRUWE4nU7jzDPP7PG+R40aZSxcuDDosb7+DceCzTAMI3aTPiIiIiKRGfR9RkRERCS5FIyIiIhIUikYERERkaRSMCIiIiJJpWBEREREkkrBiIiIiCSVghERERFJKgUjIiIiklQKRkRERCSpFIyIiIhIUikYERERkaRSMCIiIiJJ9f8BCCikD/vygPcAAAAASUVORK5CYII="
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(-10, 10, 10000)\n",
    "\n",
    "TF = Thomas_Fermi_1d(x, x0=0, amp=1, sigma=1)\n",
    "th = thermal(x, x0=0, amp=1, sigma=1)\n",
    "\n",
    "print(f'FWHM TF = {find_fwhm(TF, x)}')\n",
    "print(f'FWHM Thermal = {find_fwhm(th, x)}')\n",
    "\n",
    "plt.plot(x, TF, label='Thomas Fermi')\n",
    "plt.plot(x, th)\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:13:56.638668600Z",
     "start_time": "2023-07-25T13:13:55.916676900Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "outputs": [
    {
     "data": {
      "text/plain": "1.6110390064904865"
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "thermal(0, x0=0, amp=1, sigma=1)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:02:21.123800100Z",
     "start_time": "2023-07-25T13:02:21.081333500Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5\n"
     ]
    }
   ],
   "source": [
    "print(1**3/2)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T12:58:17.116863300Z",
     "start_time": "2023-07-25T12:58:17.096107400Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "outputs": [
    {
     "data": {
      "text/plain": "1.0"
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = 0\n",
    "sigma = 1\n",
    "\n",
    "np.exp(-0.5 * (0-x0)**2 / sigma**2)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:03:00.095225200Z",
     "start_time": "2023-07-25T13:03:00.053138800Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.642932065514894\n"
     ]
    }
   ],
   "source": [
    "print(polylog(2,1))"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:06:46.781843300Z",
     "start_time": "2023-07-25T13:06:46.758628100Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.2166174009154453\n"
     ]
    }
   ],
   "source": [
    "print(np.sqrt(1-4**(-1/3))*2)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:23:29.360679300Z",
     "start_time": "2023-07-25T13:23:29.311184700Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}