{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1e7fff1e",
   "metadata": {},
   "source": [
    "作业1:• 单位平方度内，有多少r<20.5 mag 的星系？ • SDSS给出的本地星系的Schechter光度函数的参数为Mr*=-20.83，α=-1.0，请问单位体积范围内Mr=-16 mag的星系数目是Mr=-22mag星系数目的多少倍？SDSS是星等极限样本（r<17.77）,请问在SDSS的观测样本中，Mr=-22 mag的星系数目是Mr=-16mag星系数目的多少倍？ • 根据以上光度函数，画出一个Mr=20 mag星系的红移分布概率不考虑星系的K改正和演化效应"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2107561e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number density of galaxies with r < 20.5 mag: 4.17e+00 per steradian\n",
      "Ratio of galaxy counts (M_r = -16 / M_r = -22): 18.65\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\lenovo\\AppData\\Local\\Temp\\ipykernel_21756\\3770505998.py:31: RuntimeWarning: divide by zero encountered in log10\n",
      "  return m - 5 * np.log10(d_l * 1e6) + 5\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAHWCAYAAACVPVriAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAgq9JREFUeJzt3Xd4k9XfBvA73bule9IWCpQO9l4FKZQ9FQWVJSgCoiig7LIFFw6WKCAqoiBbBCqyZSPaQUtLy+6g0D3T5Lx/9Ne8xLbQlKRp0vtzXVyak2d809OEmyfnnEcihBAgIiIiItJTBtougIiIiIhIkxh4iYiIiEivMfASERERkV5j4CUiIiIivcbAS0RERER6jYGXiIiIiPQaAy8RERER6TUGXiIiIiLSawy8RERERKTXGHiJqmnLli2QSCS4efNmlfe5efMmJBIJPv7446duGx4eDolEotRWUlKCWbNmwcvLCwYGBhgyZIiKVVfP2LFj4ePjUyPn8vHxwdixYxWPy37Oly5dqpHzd+/eHd27d6+Rc/1XfHw8evfuDVtbW0gkEuzZs0crdRA9K02/j8o+S7ds2aKxc5B+YeAlvVIWjsr+GBkZwcPDA2PHjsW9e/e0Xd4z27RpEz766CM8//zz+O677zB9+nTExMQgPDy8ysG7LEiX/bGwsED9+vUxcOBAbN68GUVFRWqpVdW6alJtrW3MmDGIjIzEsmXL8P3336NNmzYaO1dZYJBIJFi6dGmF27z88suQSCSwsrLSWB01JTY2FrNmzUKLFi1gbW0NNzc39O/fv9J/SN27dw8jRoyAnZ0dbGxsMHjwYCQmJtZw1TVj//79GDhwIFxcXGBiYgJ7e3t069YNn3zyCbKzs7VdHpFaGGm7ACJNWLx4MXx9fVFYWIhz585hy5YtOH36NKKiomBmZqbt8qpk3rx5+OCDD5Ta/vzzT3h4eOCzzz5TtO3cuROLFi1C9+7dVboKu27dOlhZWaGoqAj37t3D4cOHMX78eKxevRoHDhyAl5eXYtuNGzdCLperVH9MTEy16oqLi4OBgWb/Lf6k2o4cOaLRc1emoKAAZ8+exdy5czF16tQaO6+ZmRl++uknzJs3T6k9Ly8Pe/fu1Zn3y9N88803+PbbbzF8+HBMnjwZWVlZ2LBhAzp06IBDhw4hNDRUsW1ubi569OiBrKwszJkzB8bGxvjss88QEhKCq1evwsHBQYuvRH3kcjlee+01bNmyBcHBwZg8eTK8vLyQk5ODs2fPYt68eTh48CCOHj2q7VLL8fb2RkFBAYyNjbVdCukIBl7SS3379lVcHZswYQIcHR2xcuVK7Nu3DyNGjNBydVVjZGQEIyPlt2haWhrs7OzUcvznn38ejo6OiscLFizAjz/+iNGjR+OFF17AuXPnFM9p+i8VIQQKCwthbm4OU1NTjZ7raUxMTLRy3gcPHgCA2voXKA2tlpaWT9ymX79+2LVrF/755x80b95c0b53714UFxejT58++PPPP2ukFk0aOXIkwsPDla5Wjx8/Hk2bNkV4eLhS4F27di3i4+Nx4cIFtG3bFkDpZ0pQUBA++eQTLF++vMbr14RVq1Zhy5YtmD59Oj755BOlIVRvv/02kpOTsXXrVi1WWDmJRKI3/xijmsEhDVQndO3aFQBw48YNpfbY2Fg8//zzsLe3h5mZGdq0aYN9+/aV2z86OhrPPfcczM3N4enpiaVLl1Z4xfPSpUsICwuDo6MjzM3N4evri/Hjx1dY09dff42GDRvC1NQUbdu2xcWLF5Wef3wMb9nXz8eOHUN0dLTiq+gtW7bghRdeAAD06NFD0X78+HGVf0ZA6VfYEyZMwPnz5xEREaFor2gM7/bt29G6dWtYW1vDxsYGwcHB+PzzzwHgqXX5+PhgwIABOHz4MNq0aQNzc3Ns2LBB8dzjY3jL5Ofn44033oCDgwNsbGwwevRoZGRkKG0jkUgQHh5ebt/Hj/m02ioae5iWlobXXnsNLi4uMDMzQ/PmzfHdd98pbfP4+Oyn9e1/hYeHw9vbGwAwc+ZMSCQSpZ/333//jb59+8LGxgZWVlbo2bOn0j9Iyl6XRCLBiRMnMHnyZDg7O8PT0/OJ5wWAjh07wtfXF9u2bVNq//HHH9GnTx/Y29s/9RgVvR6JRIKYmBiMGjUK9erVQ5cuXVQ+jjq1bt263NAMBwcHdO3aFdeuXVNq37lzJ9q2basIuwDg7++Pnj174pdffnnquSQSCaZOnYodO3YgICAA5ubm6NixIyIjIwEAGzZsgJ+fH8zMzNC9e/dyQ2tOnTqFF154AfXr14epqSm8vLwwffp0FBQUlDtX2TnMzMwQFBSE3bt3V2nMfX5+PlauXInAwEB89NFH5eYLAICbmxvef/99pbbNmzfjueeeg7OzM0xNTREQEIB169Y99WdSXFyMBQsWoHXr1rC1tYWlpSW6du2KY8eOKW23cOFCGBgYlLuq/Prrr8PExAT//PMPgMrH8FblM10qlWLRokVo1KgRzMzM4ODggC5duih95pH+4RVeqhPK/kKpV6+eoi06OhqdO3eGh4cHPvjgA1haWuKXX37BkCFD8Ouvv2Lo0KEAgJSUFPTo0QMlJSWK7b7++muYm5srnSMtLQ29e/eGk5MTPvjgA9jZ2eHmzZvYtWtXuXq2bduGnJwcvPHGG5BIJFi1ahWGDRuGxMTECq+mOjk54fvvv8eyZcuQm5uLFStWAAAaNWqEadOm4YsvvsCcOXPQtGlTAFD8tzpeffVVfP311zhy5Ah69epV4TYREREYOXIkevbsiZUrVwIArl27hjNnzuDtt99Gt27dnlpXXFwcRo4ciTfeeAMTJ05EkyZNnljX1KlTYWdnh/DwcMTFxWHdunW4desWjh8/XuFf1pWpSm2PKygoQPfu3ZGQkICpU6fC19cXO3bswNixY5GZmYm3335baXtV+xYAhg0bBjs7O0yfPh0jR45Ev379FOEsOjoaXbt2hY2NDWbNmgVjY2Ns2LAB3bt3x4kTJ9C+fXulY02ePBlOTk5YsGAB8vLyqvQzGTlyJH744Qd8+OGHkEgkSE9Px5EjR/D999/j0KFDVTpGRV544QU0atQIy5cvhxCi0u3kcjkePXpUpWPa2tqq9RuHlJQUpW865HI5/v333wr/odquXTscOXIEOTk5sLa2fuJxT506hX379mHKlCkAgBUrVmDAgAGYNWsW1q5di8mTJyMjIwOrVq3C+PHjla6i79ixA/n5+XjzzTfh4OCACxcu4Msvv8Tdu3exY8cOxXa//fYbXnzxRQQHB2PFihXIyMjAa6+9Bg8Pj6e+7tOnTyMzMxMzZsyAoaHhU7cvs27dOgQGBmLQoEEwMjLC/v37MXnyZMjlcsVrrUh2dja++eYbjBw5EhMnTkROTg6+/fZbhIWF4cKFC2jRogWA0qFc+/fvx2uvvYbIyEhYW1vj8OHD2LhxI5YsWaL0LcR/VfUzPTw8HCtWrMCECRPQrl07ZGdn49KlS7hy5Uqln3mkBwSRHtm8ebMAIP744w/x4MEDcefOHbFz507h5OQkTE1NxZ07dxTb9uzZUwQHB4vCwkJFm1wuF506dRKNGjVStL3zzjsCgDh//ryiLS0tTdja2goAIikpSQghxO7duwUAcfHixUrrS0pKEgCEg4ODePTokaJ97969AoDYv3+/om3hwoXiv2/RkJAQERgYqNS2Y8cOAUAcO3asSj+jsuM+ePCgwuczMjIEADF06FBF25gxY4S3t7fi8dtvvy1sbGxESUlJped5Ul3e3t4CgDh06FCFz40ZM0bxuKxPW7duLYqLixXtq1atEgDE3r17FW0AxMKFC596zCfVFhISIkJCQhSPV69eLQCIH374QdFWXFwsOnbsKKysrER2drYQQrW+rUjZ/h999JFS+5AhQ4SJiYm4ceOGou3+/fvC2tpadOvWTdFW9nPq0qXLE/ulovNFRUUJAOLUqVNCCCHWrFkjrKysRF5enhgzZoywtLR86vEeV/Y7NnLkyCptX1ZLVf5U9fe8Kk6ePCkkEomYP3++ou3BgwcCgFi8eHG57desWSMAiNjY2CceF4AwNTVVfDYIIcSGDRsEAOHq6qr4nRFCiNmzZyt9jgghRH5+frljrlixQkgkEnHr1i1FW3BwsPD09BQ5OTmKtuPHjwsASu/Xinz++ecCgNizZ49Se0lJiXjw4IHSH7lc/sTawsLCRIMGDZTa/vs+KikpEUVFRUrbZGRkCBcXFzF+/Hil9sjISGFiYiImTJggMjIyhIeHh2jTpo2QSqWKbcp+ZzZv3qxoq+pnevPmzUX//v2f8NMhfcQhDaSXQkND4eTkBC8vLzz//POwtLTEvn37FF/xPnr0CH/++SdGjBiBnJwcpKenIz09HQ8fPkRYWBji4+MVqzocPHgQHTp0QLt27RTHd3Jywssvv6x0zrKxlwcOHIBUKn1ifS+++KLS1eayIRe1YRZ42ZXFnJycSrexs7NDXl7eM30F6Ovri7CwsCpv//rrrytd2XvzzTdhZGSEgwcPVruGqjh48CBcXV0xcuRIRZuxsTGmTZuG3NxcnDhxQml7dfatTCbDkSNHMGTIEDRo0EDR7ubmhlGjRuH06dPlZtFPnDhRpSt2ABAYGIhmzZrhp59+AlB6lXrw4MGwsLBQuebHTZo0qUrbubq6IiIiokp/nnSFTxVpaWkYNWoUfH19MWvWLEV72bCBisaSl40ZrWhowX/17NlTaVhB2ZX44cOHK10dLmt//Pfj8W+P8vLykJ6ejk6dOkEIgb///hsAcP/+fURGRmL06NFKQzVCQkIQHBz81PrKfm/+O8wjMjISTk5OSn8ePnxYYW1ZWVlIT09HSEgIEhMTkZWVVen5DA0NFePjy67ol5SUoE2bNrhy5YrStkFBQVi0aBG++eYbhIWFIT09Hd999125OQ2PU+Uz3c7ODtHR0YiPj3/qz4n0B4c0kF5as2YNGjdujKysLGzatAknT55U+gssISEBQgjMnz8f8+fPr/AYaWlp8PDwwK1bt8p9bQyg3FfwISEhGD58OBYtWoTPPvsM3bt3x5AhQzBq1Khyf3nWr19f6XFZQPrvmFRtyM3NBYAnfmU7efJk/PLLL+jbty88PDzQu3dvjBgxAn369KnyeXx9fVWqq1GjRkqPrays4ObmpvGlxW7duoVGjRqVWzmibAjErVu3lNrV2bcPHjxAfn5+hcM9mjZtCrlcjjt37iAwMFDRrurPtcyoUaPwySefYPr06fjrr78wZ86cah3ncVWtxczMTGnS2LNISUlRemxra1tu+FFeXh4GDBiAnJwcnD59Win0lW1b0fJ8hYWFSts8yX9/D2xtbQFAafWTx9sf//24ffs2FixYgH379pX7vSkLlWW/d35+fuXO7efnVy5E/lfZ+7vs/f74vmX/kN26dSu+//57pefPnDmDhQsX4uzZs8jPzy9XW9nrqch3332HTz75BLGxsUoXBSr6PZk5cya2b9+OCxcuYPny5QgICHji61HlM33x4sUYPHgwGjdujKCgIPTp0wevvvoqmjVr9sRzkG5j4CW91K5dO8UqDUOGDEGXLl0watQoxMXFwcrKSjHhbMaMGZVeZazoL5InkUgk2LlzJ86dO4f9+/crlvn65JNPcO7cOaW/VCu7AieeMM6xpkRFRQF48ut3dnbG1atXcfjwYfz+++/4/fffsXnzZowePbrcZK7KVCU0qItMJquxc2m7b6v7cx05ciRmz56NiRMnwsHBAb17966xWmQymWKViqext7d/4koabm5uSo83b96sNAmyuLgYw4YNw7///ovDhw8jKCio3PFNTU2RnJxc7thlbe7u7k+ts7Lfg6f9fshkMvTq1QuPHj3C+++/D39/f1haWuLevXsYO3asyssDVsbf3x9A6ft98ODBinYrKyvFPz5Onz6ttM+NGzfQs2dP+Pv749NPP4WXlxdMTExw8OBBfPbZZ0+s7YcffsDYsWMxZMgQzJw5E87OzjA0NMSKFSvKTSYGSq94l12BLZvs9ySqfKZ369YNN27cwN69e3HkyBF88803+Oyzz7B+/XpMmDDhqeci3cTAS3qv7EO1R48e+Oqrr/DBBx8ovh42NjZ+6pUlb2/vCr/6iouLq3D7Dh06oEOHDli2bBm2bduGl19+Gdu3b9fYB6kqE7aqouyKztOGG5iYmGDgwIEYOHAg5HI5Jk+ejA0bNmD+/Pnw8/NTe13x8fHo0aOH4nFubi6Sk5PRr18/RVu9evWQmZmptF9xcXG58KJKbd7e3vj3338hl8uVrvLGxsYqntcUJycnWFhYVPi7FhsbCwMDg3JXDKurfv366Ny5M44fP64YLlJT7ty5U+WrwceOHXviHbz+O8zm8avfcrkco0ePxtGjR/HLL78gJCSk3P4GBgYIDg6u8IYU58+fR4MGDZ46Ye1ZREZG4vr16/juu+8wevRoRft/X1fZ711CQkK5Y1TU9l9du3aFra0ttm/fjtmzZ1dp7ev9+/ejqKgI+/btU7qC/d+VFiqyc+dONGjQALt27VJ6/y1cuLDctnK5HGPHjoWNjQ3eeecdLF++HM8//zyGDRtW6fFV+UwHSv9hM27cOIwbNw65ubno1q0bwsPDGXj1GMfwUp3QvXt3tGvXDqtXr0ZhYSGcnZ3RvXt3bNiwocIrOY9fberXrx/OnTuHCxcuKD3/448/Ku2TkZFR7ipe2cxjdd29rCJla5v+N+hVx7Zt2/DNN9+gY8eO6NmzZ6XbPT6mDygNCWVfB5a9VnXWBZQu4/b416Dr1q1DSUkJ+vbtq2hr2LAhTp48WW6//17hVaW2fv36ISUlBT///LOiraSkBF9++SWsrKwqDE3qYmhoiN69e2Pv3r1KQzdSU1Oxbds2dOnSBTY2Nmo739KlS7Fw4UK89dZbajtmVahzDG9oaKjSn8ev+L711lv4+eefsXbt2ieGp+effx4XL15UCr1xcXH4888/FUvaaUrZFeDHP0uEEIol/8q4u7sjKCgIW7duVRqWcOLEiSpdEbWwsMCsWbMQFRWFDz74oMJvIP7bVlFtWVlZ2Lx5c7Ve1/nz53H27Nly23766af466+/8PXXX2PJkiXo1KkT3nzzTaSnp1d6fFU+0//7+WVlZQU/Pz+Nfk6T9vEKL9UZM2fOxAsvvIAtW7Zg0qRJWLNmDbp06YLg4GBMnDgRDRo0QGpqKs6ePYu7d+8q1nucNWsWvv/+e/Tp0wdvv/22Ylmysit/Zb777jusXbsWQ4cORcOGDZGTk4ONGzfCxsZG6SqkurVo0QKGhoZYuXIlsrKyYGpqqlgn80l27twJKysrFBcXK+60dubMGTRv3lxp6aOKTJgwAY8ePcJzzz0HT09P3Lp1C19++SVatGihGNta3boqU1xcjJ49e2LEiBGIi4vD2rVr0aVLFwwaNEiprkmTJmH48OHo1asX/vnnHxw+fFhp2SlVa3v99dexYcMGjB07FpcvX4aPjw927tyJM2fOYPXq1Rq92geUhtCIiAh06dIFkydPhpGRETZs2ICioiKsWrVKrecKCQnRaICvjDrH8FZm9erVWLt2LTp27AgLCwv88MMPSs8PHTpU8Q+hyZMnY+PGjejfvz9mzJgBY2NjfPrpp3BxccF7772n0Tr9/f3RsGFDzJgxA/fu3YONjQ1+/fXXCseAL1++HIMHD0bnzp0xbtw4ZGRk4KuvvkJQUFC5sbkV+eCDD3Dt2jV89NFHOHLkCIYPHw5PT09kZGTgypUr2LFjB5ydnRWT9Xr37q34ZueNN95Abm4uNm7cCGdn5wpD5uMGDBiAXbt2YejQoejfvz+SkpKwfv16BAQEKNV67do1zJ8/H2PHjsXAgQMBlK4x3aJFC8XcgcpU9TM9ICAA3bt3R+vWrWFvb49Lly5h586dNXqHQ9ICrawNQaQhZUszVbQ0mEwmEw0bNhQNGzZULNt048YNMXr0aOHq6iqMjY2Fh4eHGDBggNi5c6fSvv/++68ICQkRZmZmwsPDQyxZskR8++23SssJXblyRYwcOVLUr19fmJqaCmdnZzFgwABx6dIlxXEqW3pKiPJLalV1WTIhhNi4caNo0KCBMDQ0fOrSTWXHLftjZmYmPD09xYABA8SmTZuUlvQp899lyXbu3Cl69+4tnJ2dhYmJiahfv7544403RHJycpXq8vb2rnRZoMqWJTtx4oR4/fXXRb169YSVlZV4+eWXxcOHD5X2lclk4v333xeOjo7CwsJChIWFiYSEhHLHfFJt/11OSQghUlNTxbhx44Sjo6MwMTERwcHBSsshCaFa31bkSftfuXJFhIWFCSsrK2FhYSF69Ogh/vrrL6VtnvS7r+r5Hvcsy5JVtvSdNowZM+aJy509viyYEELcuXNHPP/888LGxkZYWVmJAQMGiPj4+CqdC4CYMmWKUltlP+9jx44JAGLHjh2KtpiYGBEaGiqsrKyEo6OjmDhxovjnn3/KLcMlhBDbt28X/v7+wtTUVAQFBYl9+/aJ4cOHC39//yr/bHbv3i369esnnJychJGRkbCzsxNdunQRH330kcjMzFTadt++faJZs2bCzMxM+Pj4iJUrV4pNmzaV+xn+930kl8vF8uXLhbe3tzA1NRUtW7YUBw4cUPpsKSkpEW3bthWenp7lzlu2jNrPP/+s9PP878+jKp/pS5cuFe3atRN2dnbC3Nxc+Pv7i2XLlikte0j6RyJELZglQ0RERGrRokULODk58c5hRI/hGF4iIiIdJJVKUVJSotR2/Phx/PPPP0+c2EdUF/EKLxERVUlubu5Tx4Y6OTmpfOMLqp6bN28iNDQUr7zyCtzd3REbG4v169fD1tYWUVFRcHBw0HaJRLUGJ60REVGVfPzxx1i0aNETt0lKSlK6wxhpTr169dC6dWt88803ePDgASwtLdG/f398+OGHDLtE/8ErvEREVCWJiYlPvUVyly5dFLP6iYhqCwZeIiIiItJrnLRGRERERHqNY3grIJfLcf/+fVhbW6v99qhERERE9OyEEMjJyYG7u/tTb4/NwFuB+/fvq+3+9ERERESkOXfu3IGnp+cTt2HgrUDZrULv3Lmj1vvUV0YqleLIkSPo3bs3jI2NNX4+Uj/2oe5jH+o+9qFuY//pvpruw+zsbHh5eVXpFu8MvBUoG8ZgY2NTY4HXwsICNjY2fJPrKPah7mMf6j72oW5j/+k+bfVhVYafctIaEREREek1Bl4iIiIi0msMvERERESk1ziGt5qEECgpKYFMJnvmY0mlUhgZGaGwsFAtxyPVGBoawsjIiEvQERER6SkG3mooLi5GcnIy8vPz1XI8IQRcXV1x584dhi4tsbCwgJubG0xMTLRdChEREakZA6+K5HI5kpKSYGhoCHd3d5iYmDxzSJXL5cjNzYWVldVTF04m9RJCoLi4GA8ePEBSUhIaNWrEPiAiItIzDLwqKi4uhlwuh5eXFywsLNRyTLlcjuLiYpiZmTFsaYG5uTmMjY1x69YtRT8QERGR/mC6qiYGU/3C/iQiItJf/FueiIiIiPQaAy8RERER6TUGXiIiIiLSawy8dcjYsWMhkUgwadKkcs9NmTIFEokEY8eOrfnCVLRixQq0bdsW1tbWcHZ2xpAhQxAXF6e0TWFhIaZMmQIHBwdYWVlh+PDhSE1N1VLFREREpE0MvHWMl5cXtm/fjoKCAkVbYWEhtm3bhvr161f7uMXFxeoor0pOnDiBKVOm4Ny5c4iIiIBUKkXv3r2Rl5en2Gb69OnYv38/duzYgRMnTuD+/fsYNmxYjdVIRERU1wghUFhL75/FZcnUQAiBAmn1e1gul6OgWAaj4hKVVgswNzZUeQ3gVq1a4caNG9i1axdefvllAMCuXbtQv359+Pr6Vvk43bt3R1BQEIyMjPDDDz8gODgYx44dU6mW6jp06JDS4y1btsDZ2RmXL19Gt27dkJWVhW+//Rbbtm3Dc889BwDYvHkzmjZtinPnzqFDhw41UicREVFdcethHubvicT9VAMMHSC0XU45DLxqUCCVIWDB4Ro/b8ziMFiYqN6F48ePx+bNmxWBd9OmTRg3bhyOHz+u0nG+++47vPnmmzhz5kyl2/z444944403nnic33//HV27dlXp3I/LysoCANjb2wMALl++DKlUitDQUMU2/v7+qF+/Ps6ePcvAS0REpCZFJTJ8fSIRXx1LQFGJHIYSCeJScxHsZa/t0pQw8NZBr7zyCmbPno1bt24BAM6cOYPt27erHHgbNWqEVatWPXGbQYMGoX379k/cxsPDQ6XzPk4ul+Odd95B586dERQUBABISUmBiYkJ7OzslLZ1cXFBSkpKtc9FRERE/++vG+mYtycKiQ9KhxR2amiPHtZp8He11nJl5THwqoG5sSFiFodVe3+5XI6c7BxY21irPKShOpycnNC/f39s2bIFQgj0798fjo6OKh+ndevWT93G2toa1tbV+8U/deoU+vbtq3i8YcMGxVXpMlOmTEFUVBROnz5drXMQERGRah7kFGH5wWvY/fc9AICjlSnmD2iKvgFO+P3337VcXcUYeNVAIpFUa2hBGblcjhITQ1iYGNXYHb/Gjx+PqVOnAgDWrFlTrWNYWlo+dZtnGdLQpk0bXL16VfHYxcVF6fmpU6fiwIEDOHnyJDw9PRXtrq6uKC4uRmZmptJV3tTUVLi6uj61ZiIiIipPLhf46eJtrPw9FtmFJZBIgFc7eOO93k1ga24MqVSq7RIrxcBbR/Xp0wfFxcWQSCQIC6v+1emneZYhDebm5vDz8yvXLoTAW2+9hd27d+P48ePlJtu1bt0axsbGOHr0KIYPHw4AiIuLw+3bt9GxY8dqvhIiIqK6K/p+FubujsLVO5kAgCAPGywbEozmXnZarauqGHjrKENDQ1y7dk3x/5ryLEMaKjNlyhRs27YNe/fuhbW1tWJcrq2tLczNzWFra4vXXnsN7777Luzt7WFjY4O33noLHTt25IQ1IiIiFeQWleDTI9ex5a8kyAVgZWqEGb0b49WOPjA0UG2lKG1i4K3DbGxstF1Ctaxbtw5A6dJoj9u8ebPixhmfffYZDAwMMHz4cBQVFSEsLAxr166t4UqJiIh0kxACh6JSsGh/DFKyCwEAA5q5Yf6AALjYmGm5OtUx8NYhW7ZseeLze/bsqfKxVF3RQZ2EePr6fmZmZlizZk21xycTERHVVbcf5mPBvigcj3sAAPB2sMDiwUEIaeyk5cqqj4GXiIiIiFBUIsPGk4n48s/SNXVNDA0wKaQBJvfwg1k1V4aqLRh4qZzbt28jICCg0udjYmKe6TbEREREVLucvfEQ8/ZE4oZiTV0HLBkShIZOVlquTD0YeKkcd3d3peXAKnqeiIiIdF96bhGW/3YNuxRr6ppg/oAADGruDolEdyalPQ0DL5VjZGRU4XJgREREpB/kcoHtF+9g5aFYZBVIIZEAr7T3xoyw0jV19Q0DbzVVZeIU6Q72JxER1RUx97Mxd08k/r6dCQAIdLfBsqHBaKEja+pWBwOvioyNS//Vk5+fD3Nzcy1XQ+qSn58P4P/7l4iISN/kFpXgs4jr2PLXTcjkAlamRni3V2OM7ugNI8OaudOrtjDwqsjQ0BB2dnZIS0sDAFhYWDzzGBe5XI7i4mIUFhbW2K2FqZQQAvn5+UhLS4OdnZ1Gb8JBRESkDUIIHI5OQfi+/19Tt39w6Zq6rra6t6ZudTDwVoOrqysAKELvsxJCoKCgAObm5no1QFyX2NnZKfqViIhIX9x5lI+F+6LxZ2xpZqlvb4HFgwPRvYmzliurWQy81SCRSODm5gZnZ2dIpdJnPp5UKsXJkyfRrVs3fqWuBcbGxryyS0REeqW4RI6NpxLx5Z/xKJTKYWwowaSQhpiiB2vqVgcD7zMwNDRUS1AyNDRESUkJzMzMGHiJiIjomZxLfIh5e6KQkJYLAOjYoHRNXT9n/VhTtzoYeImIiIj0QHpuEZYfvIZdV/5/Td15/QMwuIV+ralbHQy8RERERDpMLhf4+dIdfPj7/6+pO6pdfcwK84etBb85Bhh4iYiIiHTWteRszN0diSv/W1M3wM0Gy4YGoWX9etotrJZh4CUiIiLSMXlFJfj8aDy+PZ0EmVzA0sQQ03s1xthOPnq/pm51MPASERER6ZAj0SkI3xeN+1mla+r2DXLFgoEBcLPlDbEqw8BLREREpAPuZuQjfF8M/riWCgDwrGeOxYMD8Zy/i5Yrq/20es375MmTGDhwINzdS2cP7tmzR+l5IQQWLFgANzc3mJubIzQ0FPHx8U897r179/DKK6/AwcEB5ubmCA4OxqVLlzT0KoiIiIg0RyqTY8OJG+j16Un8cS0VRgYSTO7eEBHTQxh2q0irgTcvLw/NmzfHmjVrKnx+1apV+OKLL7B+/XqcP38elpaWCAsLQ2FhYaXHzMjIQOfOnWFsbIzff/8dMTEx+OSTT1CvHgdvExERkW65dPMRBnxxGit+j0WBVIZ2PvY4+HZXzOrjD3OTuncDierS6pCGvn37om/fvhU+J4TA6tWrMW/ePAwePBgAsHXrVri4uGDPnj146aWXKtxv5cqV8PLywubNmxVtvr6+6i+eiIiISEMy84vx4e+x2H7xDgCgnoUx5vRriudbe9b5NXWro9aO4U1KSkJKSgpCQ0MVbba2tmjfvj3Onj1baeDdt28fwsLC8MILL+DEiRPw8PDA5MmTMXHixErPVVRUhKKiIsXj7OxsAKW3/FXHrYOfpuwcNXEu0gz2oe5jH+o+9qFuY/+VEkJgz9VkrDgUh4z80p/FC609MLN3I9SzMEFJSYmWK6xcTfehKueptYE3JSUFAODiojw2xcXFRfFcRRITE7Fu3Tq8++67mDNnDi5evIhp06bBxMQEY8aMqXCfFStWYNGiReXajxw5AgsLi2d4FaqJiIiosXORZrAPdR/7UPexD3VbXe6/lHxgR5IhErJLr+C6mguMaCBDQ5NbOHv8lparq7qa6sP8/Pwqb1trA291yeVytGnTBsuXLwcAtGzZElFRUVi/fn2lgXf27Nl49913FY+zs7Ph5eWF3r17w8bGRuM1S6VSREREoFevXjA25h1RdBH7UPexD3Uf+1C31eX+K5TKsPZEIr65cBNSmYCZsQGmdm+IcZ28YWKkO2vq1nQfln0jXxW1NvC6uroCAFJTU+Hm5qZoT01NRYsWLSrdz83NDQEBAUptTZs2xa+//lrpPqampjA1NS3XbmxsXKNvupo+H6kf+1D3sQ91H/tQt9W1/jsel4YFe6Nx+1Hp1crn/J2xaFAgvOxr7htmdaupPlTlHLU28Pr6+sLV1RVHjx5VBNzs7GycP38eb775ZqX7de7cGXFxcUpt169fh7e3tybLJSIiIqqy1OxCLN4fg98ikwEAbrZmWDgwEGGBLpyUpgFaDby5ublISEhQPE5KSsLVq1dhb2+P+vXr45133sHSpUvRqFEj+Pr6Yv78+XB3d8eQIUMU+/Ts2RNDhw7F1KlTAQDTp09Hp06dsHz5cowYMQIXLlzA119/ja+//rqmXx4RERGREplcYOvZm/jkyHXkFpXA0ECCcZ188E6vxrAyrbXXIXWeVn+yly5dQo8ePRSPy8bRjhkzBlu2bMGsWbOQl5eH119/HZmZmejSpQsOHToEMzMzxT43btxAenq64nHbtm2xe/duzJ49G4sXL4avry9Wr16Nl19+ueZeGBEREdF//Hs3E3N2RyLqXunY0xZedlg2NAiB7rZarkz/aTXwdu/eHUKISp+XSCRYvHgxFi9eXOk2N2/eLNc2YMAADBgwQB0lEhERET2T7EIpPjkch63nbkEIwMbMCO/39cfItvVhYMDhCzWB186JiIiINEAIgQP/JmPxgRg8yCld739IC3fM7R8AJ+vyk+VJcxh4iYiIiNTsZnoe5u+Nwqn40mGXDRwtsWRIEDr7OWq5srqJgZeIiIhITYpKZNhwIhFfHUtAcYkcJkYGmNLdD5O6N4CpkaG2y6uzGHiJiIiI1OCvhHTM2xOFxPQ8AEDXRo5YPDgIvo6WWq6MGHiJiIiInsGDnCIsP3gNu/++BwBwtDLFgoEBGNjMjWvq1hIMvERERETVIJcL/HTxNlb+HovswhJIJMCrHbzxXu8msDWvO3eL0wUMvEREREQqirmfjbl7IvH37UwAQKC7DZYPDUZzLzut1kUVY+AlIiIiqqK8ohKs/uM6Np25CZlcwMrUCO/1boxXO3jDyNBA2+VRJRh4iYiIiJ5CCIEjMakI3xeN5KxCAED/YDfMHxAAV1uzp+xN2sbAS0RERPQEdzPyEb4vGn9cSwMAeNmbY/GgIPTwd9ZyZVRVDLxEREREFZDK5Pj2dBI+/yMeBVIZjA0leL1bA0zt0QjmJlxTV5cw8BIRERH9x8WbjzBvdxTiUnMAAO187bFsSBAauVhruTKqDgZeIiIiov/JyCvGh7/H4udLdwAA9pYmmNOvKYa38uCaujqMgZeIiIjqPCEEdl6+ixW/x+JRXjEA4KW2Xni/jz/qWZpouTp6Vgy8REREVKclpOVgzu4oXEh6BABo4mKNZUOD0MbHXsuVkbow8BIREVGdVFAsw1fH4vH1yURIZQLmxoZ4J7QRxnfxhTHX1NUrDLxERERU5xyLTcOCfVG486gAABDa1BnhgwLhWc9Cy5WRJjDwEhERUZ2RklWIxQeicTAyBQDgZmuG8EGB6B3gwklpeoyBl4iIiPReiUyOrWdv4ZMjccgrlsHQQILxnX3wTmhjWJoyDuk79jARERHptat3MjF3dySi72cDAFrWt8OyIcEIcLfRcmVUUxh4iYiISC9lFUjx8eE4/HD+FoQAbMyM8EHfpniprRcMDDh8oS5h4CUiIiK9IoTA/n+TseRADB7kFAEAhrX0wJz+TeFoZarl6kgbGHiJiIhIbySl52HB3iicik8HADRwssTSwUHo5Oeo5cpImxh4iYiISOcVlciw/ngi1hxPQHGJHCZGBpjaww9vhDSAqZGhtssjLWPgJSIiIp12JiEd8/dEITE9DwDQtZEjlgwOgo+jpZYro9qCgZeIiIh00oOcIiz7LQZ7rt4HADhZm2LBgAAMaObGNXVJCQMvERER6RS5XGDbhdtYeSgWOYUlkEiA0R288V5YE9iYGWu7PKqFGHiJiIhIZ0Tfz8Lc3VG4eicTABDkYYPlQ4PRzNNOq3VR7cbAS0RERLVeblEJPou4js1nkiAXgJWpEWb0boxXO/rAkGvq0lMw8BIREVGtJYTA4ehULNofjeSsQgBA/2ZuWDAgAC42ZlqujnQFAy8RERHVSnce5SN8XzSOxqYBAOrbW2Dx4EB0b+Ks5cpI1zDwEhERUa0ilcnxzakkfH70OgqlchgbSvBGt4aY+pwfzIy5pi6pjoGXiIiIao2LNzMQfuAarqfmAgA6NLDH0iFB8HO21nJlpMsYeImIiEjrHuUVY1uCAc6fvQgAsLc0wdx+TTGslQfX1KVnxsBLREREWiOEwI7Ld7Hi4DVk5BsAAEa288L7ffxhZ2Gi5epIXzDwEhERkVZcT83BvN1RuHDzEQDAzUJg9cvt0b6hk5YrI33DwEtEREQ1qqBYhi/+jMfGk4kokQuYGxti2nMN4ZIZg1b17bRdHukhBl4iIiKqMX/GpmLB3mjczSgAAIQ2dcGiwYFwtjTCwYMxWq6O9JWBNk9+8uRJDBw4EO7u7pBIJNizZ4/S80IILFiwAG5ubjA3N0doaCji4+OrfPwPP/wQEokE77zzjnoLJyIiIpUkZxVg0veXMX7LJdzNKIC7rRm+frU1vhnTBh525touj/ScVgNvXl4emjdvjjVr1lT4/KpVq/DFF19g/fr1OH/+PCwtLREWFobCwsKnHvvixYvYsGEDmjVrpu6yiYiIqIpKZHJ8ezoJoZ+cwKHoFBgaSPBGtwaIeDcEvQNdtV0e1RFaHdLQt29f9O3bt8LnhBBYvXo15s2bh8GDBwMAtm7dChcXF+zZswcvvfRSpcfNzc3Fyy+/jI0bN2Lp0qUaqZ2IiIie7OqdTMzZFYmY5GwAQKv6dlg2NBhN3Wy0XBnVNbV2DG9SUhJSUlIQGhqqaLO1tUX79u1x9uzZJwbeKVOmoH///ggNDa1S4C0qKkJRUZHicXZ26RtTKpVCKpU+w6uomrJz1MS5SDPYh7qPfaj72Ie1R3aBFJ/8EY+fLt6FEICtuRFm9m6MF1p5wMBAUmEfsf90X033oSrnqbWBNyUlBQDg4uKi1O7i4qJ4riLbt2/HlStXcPHixSqfa8WKFVi0aFG59iNHjsDCwqLKx3lWERERNXYu0gz2oe5jH+o+9qH2CAFcTpdgzy0D5EhLbxbR1kmOwd6FsE77F4cO/fvUY7D/dF9N9WF+fn6Vt621gbc67ty5g7fffhsREREwMzOr8n6zZ8/Gu+++q3icnZ0NLy8v9O7dGzY2mv/aRSqVIiIiAr169YKxsbHGz0fqxz7UfexD3cc+1K6k9DyE77+GvxJL19Rt4GiBRQMD0KGBfZX2Z//pvpruw7Jv5Kui1gZeV9fSgeypqalwc3NTtKempqJFixYV7nP58mWkpaWhVatWijaZTIaTJ0/iq6++QlFREQwNDcvtZ2pqClNT03LtxsbGNfqmq+nzkfqxD3Uf+1D3sQ9rVqFUhnXHb2Dd8RsolslhamSAt57zw8RuDWBqVP7v3Kdh/+m+mupDVc5RawOvr68vXF1dcfToUUXAzc7Oxvnz5/Hmm29WuE/Pnj0RGRmp1DZu3Dj4+/vj/fffrzDsEhERUfWcjk/H/L1RSErPAwB0a+yEJYMD4e1gqeXKiJRpNfDm5uYiISFB8TgpKQlXr16Fvb096tevj3feeQdLly5Fo0aN4Ovri/nz58Pd3R1DhgxR7NOzZ08MHToUU6dOhbW1NYKCgpTOYWlpCQcHh3LtREREVD1pOYVY9ts17L16HwDgbG2KhQMD0S/YFRKJRMvVEZWn1cB76dIl9OjRQ/G4bBztmDFjsGXLFsyaNQt5eXl4/fXXkZmZiS5duuDQoUNK43Nv3LiB9PT0Gq+diIiorpHJBbZduI1Vh2KRU1gCAwkwuqMP3uvdGNZmHIZAtZdWA2/37t0hhKj0eYlEgsWLF2Px4sWVbnPz5s0nnuP48ePVrI6IiIjKRN3Lwtw9UfjnTiYAINjDFsuGBqGZp51W6yKqilo7hpeIiIi0L7eoBJ8euY4tfyVBLgArUyPMDGuCVzp4w9CAwxdINzDwEhERUTlCCByKSsGi/TFIyS4EAAxo5ob5AwLgYlP1pT+JagMGXiIiIlJy51E+Fu6Lxp+xaQAAbwcLLB4chJDGTlqujKh6GHiJiIgIACCVyfHNqSR8fvQ6CqVyGBtKMCmkIab08IOZMZf2JN3FwEtERES4dPMR5uyOxPXUXABAhwb2WDokGH7OVlqujOjZMfASERHVYZn5xfjw91hsv3gHAGBvaYK5/ZpiWCsPrqlLeoOBl4iIqA4SQmDXlXtYdvAaHuUVAwBeauuF9/v4o56liZarI1IvBl4iIqI65saDXMzbHYWziQ8BAI1drLBsaDDa+thruTIizWDgJSIiqiMKpTKsPZaA9ScSUSyTw8zYAG/3bIzXuvjCxMhA2+URaQwDLxERUR1wKv4B5u+Jws2H+QCAHk2csHhwELzsLbRcGZHmMfASERHpsbScQiw9cA37/rkPAHCxMcXCgYHoG+TKSWlUZzDwEhER6SG5XGDbhdtYeSgWOYUlMJAAozv64L3ejWFtZqzt8ohqFAMvERGRnom5n405uyNx9U4mACDYwxbLhwYj2NNWu4URaYlKgTczMxO7d+/GqVOncOvWLeTn58PJyQktW7ZEWFgYOnXqpKk6iYiI6Cnyikqw+o/r2HTmJmRyAStTI8zo3RivdvSBoQGHL1DdVaUpmffv38eECRPg5uaGpUuXoqCgAC1atEDPnj3h6emJY8eOoVevXggICMDPP/+s6ZqJiIjoP45Ep6DXpyew8VQSZHKB/sFu+OPdEIzt7MuwS3Vela7wtmzZEmPGjMHly5cREBBQ4TYFBQXYs2cPVq9ejTt37mDGjBlqLZSIiIjKu5dZgPB90YiISQUAeNYzx5LBQejh76zlyohqjyoF3piYGDg4ODxxG3Nzc4wcORIjR47Ew4cP1VIcERERVaxEJsfmMzfx2R/XkV8sg5GBBK93a4C3nmsEcxNDbZdHVKtUKfA+Lew+6/ZERERUdVduZ2Du7ihcS84GALT1qYdlQ4PR2MVay5UR1U4qr9JgaGiIbt264ddff4W9/f/fgjA1NRXu7u6QyWRqLZCIiIhKZRVIsepQLLZduA0hADsLY8zp2xTPt/aEAcfpElVK5cArhEBRURHatGmD/fv3IzAwUOk5IiIiUi8hBPb9cx9LDlxDem4RAGB4K0/M6ecPBytTLVdHVPupfONsiUSCX3/9FQMHDkTHjh2xd+9epeeIiIhIfW6m52H0pgt4e/tVpOcWoaGTJX6a2AGfjGjOsEtURdW6wmtoaIjPP/8cgYGBePHFFzFv3jxMmDBBE/URERHVSUUlMmw4kYivjiWguEQOEyMDvNXDD6+HNICpESelEanime609vrrr6NRo0Z44YUXcPLkSXXVREREVKf9dSMd8/ZEIfFBHgCgayNHLBkcBB9HSy1XRqSbVA683t7eMDT8/39Z9ujRA+fOncPAgQPVWhgREVFd8zC3CMsOXsOuK/cAAI5WplgwMAADm7lx2CDRM1A58CYlJZVr8/Pzw99//43U1FS1FEVERFSXyOUCv1y6gxW/xyKrQAqJBHilvTdmhDWBrbmxtssj0nlVCrxCiKf+y9LMzAze3t5qKYqIiKiuiEvJwdzdkbh0KwMAEOBmg2VDg9Cyfj0tV0akP6q0SkNgYCC2b9+O4uLiJ24XHx+PN998Ex9++KFaiiMiItJXBcUyfPh7LPp/cQqXbmXAwsQQ8/o3xb6pnRl2idSsSld4v/zyS7z//vuYPHkyevXqhTZt2sDd3R1mZmbIyMhATEwMTp8+jejoaEydOhVvvvmmpusmIiLSWX/GpmLB3mjczSgAAPQOcEH4oEC425lruTIi/VSlwNuzZ09cunQJp0+fxs8//4wff/wRt27dQkFBARwdHdGyZUuMHj0aL7/8MurV479KiYiIKpKSVYhF+6Pxe1QKAMDDzhzhgwLRK8BFy5UR6TeVJq116dIFXbp00VQtREREekkmF/jur5v45Egc8oplMDSQ4LUuvni7ZyNYmj7TCqFEVAUqvcuEEEhISEBxcTGaNGkCIyO+SYmIiJ7k37uZmLM7ElH3sgEALevbYfnQYDR1s9FyZUR1R5UTa1JSEgYNGoSYmBgAgIeHB3799Ve0bdtWY8URERHpquxCKT45HIet525BCMDGzAjv9/XHyLb1YWDANXWJalKVA+/MmTNRUlKCH374AWZmZvj4448xadIkXL58WZP1ERER6RQhBA5GpmDR/mik5RQBAIa0cMfc/gFwsjbVcnVEdVOVA+/p06exc+dOxRjeDh06wNPTE3l5ebC05K0OiYiIbj/Mx4J9UTge9wAA4OtoiSWDg9ClkaOWKyOq26oceNPS0tCoUSPFYzc3N5ibmyMtLQ2+vr4aKY6IiEgXFJfIsfFUIr44Go+iEjlMDA3wZveGeLN7Q5gZG2q7PKI6r8qBVyKRIDc3F+bm/79GoIGBAXJycpCdna1os7HhIHwiIqo7LiQ9wtzdkYhPywUAdGrogCVDgtDQyUrLlRFRmSoHXiEEGjduXK6tZcuWiv+XSCSQyWTqrZCIiKgWysgrxorfr+GXS3cBAA6WJpg3oCmGtPCARMJJaUS1SZUD77FjxzRZBxERkU4QQmDn5btYfvAaMvKlAICR7bzwfh9/2FmYaLk6IqpIlQNvSEiI2k9+8uRJfPTRR7h8+TKSk5Oxe/duDBkyRPG8EAILFy7Exo0bkZmZic6dO2PdunVKY4n/a8WKFdi1axdiY2Nhbm6OTp06YeXKlWjSpIna6yciorolIS0Hc3dH4XzSIwBAExdrLB8WhNbe9lqujIiexECbJ8/Ly0Pz5s2xZs2aCp9ftWoVvvjiC6xfvx7nz5+HpaUlwsLCUFhYWOkxT5w4gSlTpuDcuXOIiIiAVCpF7969kZeXp6mXQUREeq5QKsPHh+PQ9/NTOJ/0CGbGBvigrz8OTOvCsEukA7R6q7S+ffuib9++FT4nhMDq1asxb948DB48GACwdetWuLi4YM+ePXjppZcq3O/QoUNKj7ds2QJnZ2dcvnwZ3bp1U+8LICIivXfy+gPM3xuFWw/zAQA9/Z0RPigQXvYWWq6MiKqq1t4bOCkpCSkpKQgNDVW02draon379jh79mylgfe/srKyAAD29pX/C7yoqAhFRUWKx2WrTkilUkil0uqUr5Kyc9TEuUgz2Ie6j32o+9Tdh2k5RVh+MA6/RaUAAFxsTDG/nz96BzhDIpHwd0XN+B7UfTXdh6qcp9YG3pSU/33AuLgotbu4uCieexq5XI533nkHnTt3RlBQUKXbrVixAosWLSrXfuTIEVhY1Ny/4CMiImrsXKQZ7EPdxz7Ufc/ah3IBnEmV4MBtAxTKJJBAoJurQL/6eZDduozfb6mpUKoQ34O6r6b6MD8/v8rbqhx4N2/ejBdffLFGg2B1TZkyBVFRUTh9+vQTt5s9ezbeffddxePs7Gx4eXmhd+/eNbKusFQqRUREBHr16gVjY2ONn4/Uj32o+9iHuk8dfRh9PxsL9sfg37ul3/QFe9hgyaAABLpzjXlN43tQ99V0Hz5+H4inUTnwfvDBB3j77bfxwgsv4LXXXkOnTp1UPUSVuLq6AgBSU1Ph5uamaE9NTUWLFi2euv/UqVNx4MABnDx5Ep6enk/c1tTUFKam5e9vbmxsXKNvupo+H6kf+1D3sQ91X3X6MLeoBJ9FXMfmM0mQC8DK1Agzw5rglQ7eMDTgmro1ie9B3VdTfajKOVRepeHevXv47rvvkJ6eju7du8Pf3x8rV66s8jCDqvL19YWrqyuOHj2qaMvOzsb58+fRsWPHSvcTQmDq1KnYvXs3/vzzT972mIiIKiWEwKGoFPT69AS+PV0advs3c8PR90IwppMPwy6RnlA58BoZGWHo0KHYu3cv7ty5g4kTJ+LHH39E/fr1MWjQIOzduxdyubxKx8rNzcXVq1dx9epVAKUT1a5evYrbt29DIpHgnXfewdKlS7Fv3z5ERkZi9OjRcHd3V1qrt2fPnvjqq68Uj6dMmYIffvgB27Ztg7W1NVJSUpCSkoKCggJVXyoREemxuxn5mLj1Eib9cBnJWYXwsjfHlnFtsWZUK7jYmGm7PCJSo2eatObi4oIuXbrg+vXruH79OiIjIzFmzBjUq1cPmzdvRvfu3Z+4/6VLl9CjRw/F47JxtGPGjMGWLVswa9Ys5OXl4fXXX0dmZia6dOmCQ4cOwczs/z+Ibty4gfT0dMXjdevWAUC5c2/evBljx459lpdLRER6QCqTY9PpJKz+Ix4FUhmMDSV4vVsDTO3RCOYmhtouj4g0oFqBNzU1Fd9//z02b96MxMREDBkyBAcOHEBoaCjy8vKwePFijBkzBrduPXkqa/fu3SGEqPR5iUSCxYsXY/HixZVuc/PmTaXHTzoeERHVbZdvZWDu7kjEpuQAANr52GPZ0CA0crHWcmVEpEkqB96BAwfi8OHDaNy4MSZOnIjRo0crrXFraWmJ9957Dx999JFaCyUiIqqurHwpVh6OxU8XbkMIoJ6FMWb3a4oXWntCIuE4XSJ9p3LgdXZ2xokTJ544cczJyQlJSUnPVBgREdGzEkJg79X7WPpbDNJziwEAz7f2xJx+TWFvaaLl6oiopqgceENCQtCqVaty7cXFxdi+fTtGjx4NiUQCb29vtRRIRERUHUnpeZi/JwqnE0rnefg5W2HpkCB0aOCg5cqIqKapvErDuHHjFLfrfVxOTg7GjRunlqKIiIiqq6hEjtV/XEfY6pM4nZAOUyMDzOjdGAendWXYJaqjVL7CK4SocLzT3bt3YWtrq5aiiIiIquN6lgSrv/oLSQ9LbznarbETlgwOhLeDpZYrIyJtqnLgbdmyJSQSCSQSCXr27Akjo//fVSaTISkpCX369NFIkURERE+SnluEJfujsTfGEEA+nKxNsWBAAAY0c+OkNCKqeuAtu9nD1atXERYWBisrK8VzJiYm8PHxwfDhw9VeIBERUWXkcoGfL93Bh7/HIqtACgkEXm5fH7P6NoWNGW9PS0Slqhx4Fy5cCADw8fHBiy++qHTzByIiopoWm5KNubujcPlWBgCgqas1+jpm4M0BTWFszLBLRP9P5TG8Y8aM0UQdREREVZJfXILPj8bj21NJKJELWJoY4t3eTTCqjTuOHD6k7fKIqBaqUuC1t7fH9evX4ejoiHr16j1xPNSjR4/UVhwREdHjjl5LxYK90biXWQAA6BPoioWDAuBmaw6pVKrl6oiotqpS4P3ss89gbW2t+H9OACAiopqUnFWARfticCg6BQDgYWeOxYMD0bOpi5YrIyJdUKXA+/gwhrFjx2qqFiIiIiUlMjm+O3sLnx6JQ16xDIYGEkzo6ou3ezaChYnKo/KIqI5S+dPiypUrMDY2RnBwMABg79692Lx5MwICAhAeHg4TE96qkYiInt0/dzIxZ3ckou9nAwBa1bfD8mHB8He10XJlRKRrVL7T2htvvIHr168DABITE/Hiiy/CwsICO3bswKxZs9ReIBER1S3ZhVIs2BuFIWvPIPp+NmzMjLBiWDB2TurEsEtE1aLyFd7r16+jRYsWAIAdO3YgJCQE27Ztw5kzZ/DSSy9h9erVai6RiIjqAiEEfotMxuL9MUjLKQIADG3pgbn9m8LRylTL1RGRLqvWrYXlcjkA4I8//sCAAQMAAF5eXkhPT1dvdUREVCfcfpiP+XujcOL6AwCAr6Mllg4JQmc/Ry1XRkT6QOXA26ZNGyxduhShoaE4ceIE1q1bBwBISkqCiwtnyxIRUdUVl8ix8VQivjgaj6ISOUwMDTC5R0NMCmkIM2NDbZdHRHpC5cC7evVqvPzyy9izZw/mzp0LPz8/AMDOnTvRqVMntRdIRET66XziQ8zdE4WEtFwAQGc/BywZHIQGTlZP2ZOISDUqB95mzZohMjKyXPtHH30EQ0P+a5yIiJ7sUV4xVhy8hh2X7wIAHK1MMK9/AAa3cOc670SkEdVexLC4uBhpaWmK8bxl6tev/8xFERGR/hFCYMflu1hx8Boy8kvvijaqfX28H+YPWwtjLVdHRPqsWqs0vPbaa/jrr7+U2oUQkEgkkMlkaiuOiIj0Q3xqDubuicKFpNLbz/u7WmPZ0GC09q6n5cqIqC5QOfCOGzcORkZGOHDgANzc3Pj1ExERVapQKsOXf8bj65OJkMoEzI0NMb1XI4zr7AtjQ5WXgiciqhaVA+/Vq1dx+fJl+Pv7a6IeIiLSE8fj0rBgbzRuP8oHAIQ2dUb4oEB41rPQcmVEVNeoHHgDAgK43i4REVUqLbsQiw7E4Ld/kwEAbrZmCB8UiN4BLvxWkIi0QuXAu3LlSsyaNQvLly9HcHAwjI2VJxrY2PC2j0REdZFMLvDj+Vv46FAccopKYCABxnX2xfRejWFlWu050kREz0zlT6DQ0FAAQM+ePZXaOWmNiKjuirqXhbm7I/HP3SwAQHMvOywbEoQgD1stV0ZEVI3Ae+zYMU3UQUREOii3qASfHrmOLX8lQS4Aa1MjzOrTBKPae8PQgMMXiKh2UDnwhoSEaKIOIiLSMYejUxC+LxrJWYUAgAHN3LBgQACcbcy0XBkRkbJqDao6deoUNmzYgMTEROzYsQMeHh74/vvv4evriy5duqi7RiIiqkXuZRZg4d5o/HEtFQBQ394CS4YEIaSxk5YrIyKqmMqLIP76668ICwuDubk5rly5gqKiIgBAVlYWli9frvYCiYiodiiRybHxZCJ6fXoCf1xLhbGhBFN6NMSR6d0YdomoVlM58C5duhTr16/Hxo0blVZo6Ny5M65cuaLW4oiIqHa4eicTg746g2UHryG/WIa2PvVwcFpXzAzzh5mxobbLIyJ6IpWHNMTFxaFbt27l2m1tbZGZmamOmoiIqJbILpTio0Nx+OH8LQgB2FkYY3Zff7zQ2gsGnJRGRDpC5cDr6uqKhIQE+Pj4KLWfPn0aDRo0UFddRESkRUII/BaZjEX7Y/Agp3To2rBWHpjbrykcrEy1XB0RkWpUDrwTJ07E22+/jU2bNkEikeD+/fs4e/YsZsyYgfnz52uiRiIiqkG3H+Zj/t4onLj+AADQwNESS4cEoZOfo5YrIyKqHpUD7wcffAC5XI6ePXsiPz8f3bp1g6mpKWbMmIG33npLEzUSEVENkMrk2HgqEZ//EY+iEjlMDA0wuUdDvNm9IUyNOE6XiHSXyoFXIpFg7ty5mDlzJhISEpCbm4uAgABYWVlpoj4iIqoBl24+wpzdkbiemgsA6NTQAUuHBKGBEz/biUj3VWsdXiEEsrOz4eLigoCAAHXXRERENSQzvxgf/h6L7RfvAADsLU0wr39TDG3pAYmEk9KISD+otCxZSkoKRo8ejXr16sHFxQXOzs6oV68exo8fj9TUVE3VSEREaiaEwO6/76LnJycUYfeltl74870QDGvlybBLRHqlyoE3OzsbnTp1wqFDhzBu3DisXbsWa9aswauvvor9+/eja9euyM3NVenkJ0+exMCBA+Hu7g6JRII9e/YoPS+EwIIFC+Dm5gZzc3OEhoYiPj7+qcdds2YNfHx8YGZmhvbt2+PChQsq1UVEpM8SH+TilW/PY/rP/+BhXjEaOVthx6SO+HB4M9hZmGi7PCIitaty4P38889haGiI6OhofPbZZ3jjjTcwadIkfPHFF4iOjoYQAl988YVKJ8/Ly0Pz5s2xZs2aCp9ftWoVvvjiC6xfvx7nz5+HpaUlwsLCUFhYWOkxf/75Z7z77rtYuHAhrly5gubNmyMsLAxpaWkq1UZEpG+KSmRY/cd19Fl9CmcSHsLUyAAzw5rgt2ld0dbHXtvlERFpTJUD72+//YY5c+bAyan87SOdnZ0xe/Zs7N+/X6WT9+3bF0uXLsXQoUPLPSeEwOrVqzFv3jwMHjwYzZo1w9atW3H//v1yV4If9+mnn2LixIkYN24cAgICsH79elhYWGDTpk0q1UZEpE/+upGOvp+fwuo/4lEsk6NbYydETA/BlB5+MDFS+aabREQ6pcqT1q5fv45OnTpV+nynTp0wY8YMtRQFAElJSUhJSUFoaKiizdbWFu3bt8fZs2fx0ksvldunuLgYly9fxuzZsxVtBgYGCA0NxdmzZys9V1FREYqKihSPs7OzAQBSqRRSqVQdL+eJys5RE+cizWAf6j597cOHecVYeSgOu68mAwCcrEwwt58/+gW5QCKR6NXr1dc+rCvYf7qvpvtQlfNUOfBmZ2fDzs6u0uft7OwUQVEdUlJSAAAuLi5K7S4uLorn/is9PR0ymazCfWJjYys914oVK7Bo0aJy7UeOHIGFhYWqpVdbREREjZ2LNIN9qPv0pQ/lAjifJsG+2wbIL5FAAoHOLgL96+dDcucKfr+j7Qo1R1/6sK5i/+m+murD/Pz8Km9b5cArhICBQeVfe0kkEgghqnzi2mT27Nl49913FY+zs7Ph5eWF3r17w8bGRuPnl0qliIiIQK9evWBsbKzx85H6sQ91nz71YXxaLhbsi8GlW5kAAH9XaywZ1BQtvOy0Wpem6VMf1kXsP91X032oyoVWlQJv48aNK12qRt1h19XVFQCQmpoKNzc3RXtqaipatGhR4T6Ojo4wNDQst0Raamqq4ngVMTU1halp+XvDGxsb1+ibrqbPR+rHPtR9utyHhVIZvvwzHl+fTIRUJmBubIh3ezXGuM4+MDKsO+N0dbkPif2nD2qqD1U5R5UD7+bNm6tVTHX5+vrC1dUVR48eVQTc7OxsnD9/Hm+++WaF+5iYmKB169Y4evQohgwZAgCQy+U4evQopk6dWkOVExHVvJPXH2DenijcflT6FV9oU2eEDwqEZ72aG5ZFRFRbVTnwjhkzRu0nz83NRUJCguJxUlISrl69Cnt7e9SvXx/vvPMOli5dikaNGsHX1xfz58+Hu7u7IswCQM+ePTF06FBFoH333XcxZswYtGnTBu3atcPq1auRl5eHcePGqb1+IiJtS8spxJID17D/n/sAAFcbM4QPCkRYoAtvHkFE9D/VurWwuly6dAk9evRQPC4bRztmzBhs2bIFs2bNQl5eHl5//XVkZmaiS5cuOHToEMzMzBT73LhxA+np6YrHL774Ih48eIAFCxYgJSUFLVq0wKFDh8pNZCMi0mVyucC2C7ex8lAscgpLYCABxnTywXu9m8DKVKsf7UREtY5WPxW7d+/+xLG/EokEixcvxuLFiyvd5ubNm+Xapk6dyiEMRKS3riVnY87uSPx9OxMAEOxhixXDghHkYavdwoiIaileBiAi0hH5xSX4/I94fHM6CTK5gJWpEWb0boxXO/rA0IDDF4iIKsPAS0SkA45eS8WCvdG4l1kAAOgb5IqFAwPhamv2lD2JiEjlwHvs2DGlcbdERKQ5KVmFCN8XjUPRpTfc8bAzx5IhgXjOn/MSiIiqSuXA26dPH3h6emLcuHEYM2YMvLy8NFEXEVGdJpMLbD17E58cuY7cohIYGkgwoasv3u7ZCBYm/HKOiEgVKq9Efu/ePUydOhU7d+5EgwYNEBYWhl9++QXFxcWaqI+IqM6JvJuFIWvOYNH+GOQWlaBlfTsceKsLZvdtyrBLRFQNKgdeR0dHTJ8+HVevXsX58+fRuHFjTJ48Ge7u7pg2bRr++ecfTdRJRKT3cotKsGh/NAavOY3Ie1mwNjPCsqFB+HVSJzR10/xtzomI9NUzXSpo1aoVXF1d4eDggA8//BCbNm3C2rVr0bFjR6xfvx6BgYHqqpOISG8JIXA4OgXh+2KQkl0IABjU3B3zBjSFszUnpRERPatq3VxdKpVi586d6NevH7y9vXH48GF89dVXSE1NRUJCAry9vfHCCy+ou1YiIr1zNyMfE767hEk/XEFKdiG8HSywdXw7fDGyJcMuEZGaqHyF96233sJPP/0EIQReffVVrFq1CkFBQYrnLS0t8fHHH8Pd3V2thRIR6ROpTI7NZ5LwWUQ8CqQyGBtK8Ea3hpj6nB/MjA21XR4RkV5ROfDGxMTgyy+/xLBhw2BqalrhNo6Ojjh27NgzF0dEpI+u3M7AnF2RiE3JAQC087HH8mFB8HO21nJlRET6SeXAu3DhQnTq1AlGRsq7lpSU4K+//kK3bt1gZGSEkJAQtRVJRKQPsgqkWHUoFtsu3IYQgJ2FMeb0bYrnW3vCgHdKIyLSGJUDb48ePZCcnAxnZ2el9qysLPTo0QMymUxtxRER6QMhBPb/m4zF+2OQnlsEABjeyhNz+vnDwarib8qIiEh9VA68QghIJOWvRDx8+BCWlpZqKYqISF/cepiHeXuicCo+HQDQwMkSy4YEo2NDBy1XRkRUd1Q58A4bNgwAIJFIMHbsWKXxuzKZDP/++y86deqk/gqJiHRQcYkcG08l4ouj8SgqkcPEyABTuvthUvcGMDXipDQioppU5cBra2sLoPQKr7W1NczNzRXPmZiYoEOHDpg4caL6KyQi0jEXkh5h7u5IxKflAgA6+zlg6ZBg+DryWzAiIm2ocuDdvHkzAMDHxwczZszg8AUiov/IyCvGh7/H4udLdwAADpYmmD8gAINbuFc4FIyIiGpGtVZpICKi/yeEwK4r97Ds4DU8yisGAIxs54X3+/jDzsJEy9UREVGVAm+rVq1w9OhR1KtXDy1btnzilYorV66orTgiotruxoNczNsdhbOJDwEATVyssWxoENr42Gu5MiIiKlOlwDt48GDFJLUhQ4Zosh4iIp1QKJVh3fEbWHf8BoplcpgZG+Dtno0xoasvjA2rddd2IiLSkCoF3seHMXBIAxHVdX8lpGPunigkpecBALo3ccKSwUHwsrfQcmVERFQRlcfwEhHVVem5RVj22zXs/vseAMDJ2hThAwPRL9iVk9KIiGqxKgXeevXqVfnD/NGjR89UEBFRbSOXC/xy6Q5W/B6LrAIpJBLg1Q7emBHWBDZmxtouj4iInqJKgXf16tUaLoOIqHa6npqDubsjcfFmBgAgwM0Gy4cFo4WXnXYLIyKiKqtS4B0zZoym6yAiqlUKimX48s94fH0yESVyAQsTQ7zbqzHGdvKBESelERHplCoF3uzsbNjY2Cj+/0nKtiMi0lXH49Iwf28U7jwqAAD0CnBB+KBAeNiZP2VPIiKqjao8hjc5ORnOzs6ws7OrcDyvEAISiQQymUztRRIR1YSsYuCdn//Fb1EpAAA3WzOEDwpEWKCrlisjIqJnUaXA++eff8LevnQR9WPHjmm0ICKimiaTC/x4/jY+vGqIQlkKDCTAuM6+mN6rMaxMuZgNEZGuq9IneUhISIX/T0Sk66LvZ2HO7ij8cycTgATNPGywfFgzBHnYars0IiJSk2pdusjIyMC3336La9euAQACAgIwbtw4xVVgIqLaLq+oBKv/uI5NZ25CJhewNDVEH7diLBvXHmamJtouj4iI1EjlqcYnT56Ej48PvvjiC2RkZCAjIwNffPEFfH19cfLkSU3USESkVn/EpKLXpyew8VQSZHKB/sFuODytM7q5CRga8AYSRET6RuUrvFOmTMGLL76IdevWwdDQEAAgk8kwefJkTJkyBZGRkWovkohIHZKzChC+LxqHo1MBAJ71zLFkcBB6+DtDKpVquToiItIUlQNvQkICdu7cqQi7AGBoaIh3330XW7duVWtxRETqUCKT47uzt/DpkTjkFctgZCDBhK4N8HbPRjA3MXz6AYiISKepHHhbtWqFa9euoUmTJkrt165dQ/PmzdVWGBGROvx7NxOzd0Ui+n7pGuKtveth2dAg+LtyzXAiorqiSoH333//Vfz/tGnT8PbbbyMhIQEdOnQAAJw7dw5r1qzBhx9+qJkqiYhUlFMoxSdHrmPr2ZuQC8DGzAgf9G2Kl9p6wYDjdImI6pQqBd4WLVpAIpFACKFomzVrVrntRo0ahRdffFF91RERqUgIgd+jUrBofzRSs4sAAENauGNu/wA4WZtquToiItKGKgXepKQkTddBRPTM7jzKx4K9UTgW9wAA4ONggSVDgtC1kZOWKyMiIm2qUuD19vbWdB1ERNUmlcnx7ekkrP7jOgqlchgbSvBmSENM7uEHM2NOSiMiquuqfc/MmJgY3L59G8XFxUrtgwYNeuaiiIiq6vKtDMzdHYnYlBwAQHtfeywbGgQ/Z2stV0ZERLWFyoE3MTERQ4cORWRkpNK4XomkdBKITCZTa4E5OTmYP38+du/ejbS0NLRs2RKff/452rZtW+k+P/74I1atWoX4+HjY2tqib9+++Oijj+Dg4KDW2ohIe7LypVh5OBbbzt8GANSzMMacfk3xfGtPxecRERERUI07rb399tvw9fVFWloaLCwsEB0djZMnT6JNmzY4fvy42gucMGECIiIi8P333yMyMhK9e/dGaGgo7t27V+H2Z86cwejRo/Haa68hOjoaO3bswIULFzBx4kS110ZENU8Igb1X76Hnp8cVYfeF1p44+l53vNDGi2GXiIjKUfkK79mzZ/Hnn3/C0dERBgYGMDAwQJcuXbBixQpMmzYNf//9t9qKKygowK+//oq9e/eiW7duAIDw8HDs378f69atw9KlSyusz8fHB9OmTQMA+Pr64o033sDKlSvVVhcRacfN9DzM3xuFU/HpAICGTpZYNjQYHRrw2xsiIqqcyoFXJpPB2rp0bJyjoyPu37+PJk2awNvbG3FxcWotrqSkBDKZDGZmZkrt5ubmOH36dIX7dOzYEXPmzMHBgwfRt29fpKWlYefOnejXr1+l5ykqKkJRUZHicXZ26QL1Uqm0Rm43WnYO3tpUd7EPNauoRI5vTt/E2hOJKC6Rw8TIAJNDGmBCFx+YGhmo5efOPtR97EPdxv7TfTXdh6qcRyIeX1y3Crp27Yr33nsPQ4YMwahRo5CRkYF58+bh66+/xuXLlxEVFaVywU/SqVMnmJiYYNu2bXBxccFPP/2EMWPGwM/Pr9KAvWPHDowfPx6FhYUoKSnBwIED8euvv8LY2LjC7cPDw7Fo0aJy7du2bYOFhYVaXw8RqSYhG/gl0RCpBaVDFRrbyjHCVw4ncy0XRkREWpWfn49Ro0YhKysLNjZPvnumyoH38OHDyMvLw7Bhw5CQkIABAwbg+vXrcHBwwM8//4znnnvumYr/rxs3bmD8+PE4efIkDA0N0apVKzRu3BiXL1/GtWvXym0fExOD0NBQTJ8+HWFhYUhOTsbMmTPRtm1bfPvttxWeo6IrvF5eXkhPT3/qD1AdpFIpIiIi0KtXr0pDOdVu7EP1e5RXjFVHruPXK/cBAA6WJpjTtwkGNnPVyDhd9qHuYx/qNvaf7qvpPszOzoajo2OVAq/KQxrCwsIU/+/n54fY2Fg8evQI9erV08hfQg0bNsSJEyeQl5eH7OxsuLm54cUXX0SDBg0q3H7FihXo3LkzZs6cCQBo1qwZLC0t0bVrVyxduhRubm7l9jE1NYWpafk7MBkbG9fom66mz0fqxz58dkII/HrlHpb9FoOM/NKvq0a1r4/3w/xha6H5ny37UPexD3Ub+0/31VQfqnKOaq/DCwB37twBAHh5eT3LYarE0tISlpaWyMjIwOHDh7Fq1aoKt8vPz4eRkfLLMjQsXXhexYvZRFTDEtJyMW9PJM4lPgIANHGxxvJhQWjtba/lyoiISJepvCxZSUkJ5s+fD1tbW/j4+MDHxwe2traYN2+eRgYpHz58GIcOHUJSUhIiIiLQo0cP+Pv7Y9y4cQCA2bNnY/To0YrtBw4ciF27dmHdunVITEzEmTNnMG3aNLRr1w7u7u5qr4+Inl2hVIZPj8Sh7+cncS7xEcyMDfBBX38cmNaFYZeIiJ6Zyld433rrLezatQurVq1Cx44dAZQuBRYeHo6HDx9i3bp1ai0wKysLs2fPxt27d2Fvb4/hw4dj2bJlisvYycnJuH37tmL7sWPHIicnB1999RXee+892NnZ4bnnnuOyZES11On4dMzbE4mbD/MBAD2aOGHx4CB42XPCKBERqYfKgXfbtm3Yvn07+vbtq2hr1qwZvLy8MHLkSLUH3hEjRmDEiBGVPr9ly5ZybW+99RbeeusttdZBROqVnluEpQdisOdq6aQ0Z2tThA8KRN8gzUxKIyKiukvlwGtqagofH59y7b6+vjAxMVFHTUSkx+Ryge0X7+DD368hu7AEEgkwuoM33gtrAhszTlQhIiL1UznwTp06FUuWLMHmzZsVKxsUFRVh2bJlmDp1qtoLJCL9EZeSgzm7I3H5VgYAINDdBsuHBqO5l512CyMiIr1WpcA7bNgwpcd//PEHPD090bx5cwDAP//8g+LiYvTs2VP9FRKRzisoluHzo/H45lQiSuQCFiaGeK93E4zp6A0jQ5XnzhIREamkSoHX1tZW6fHw4cOVHtfEsmREpJuOxaZh/t4o3M0oAAD0DnBB+KBAuNvxVmlERFQzqhR4N2/erOk6iEjPpGYXYvH+GPwWmQwAcLc1w6LBQegV4KLlyoiIqK6p9o0nHjx4gLi4OABAkyZN4OTkpLaiiEh3yeQCP5y7hY8PxyGnqASGBhKM6+SD6b0aw9L0me51Q0REVC0q/+2Tl5eHt956C1u3boVcLgdQeiez0aNH48svv4SFBdfOJKqrou5lYe7uSPxzNwsA0NzLDsuHBiHQ3fYpexIREWmOyrNF3n33XZw4cQL79+9HZmYmMjMzsXfvXpw4cQLvvfeeJmokolour6gESw7EYNBXp/HP3SxYmxphyeBA7HqzE8MuERFpncpXeH/99Vfs3LkT3bt3V7T169cP5ubmGDFihNpvPEFEtduR6BQs3BeN5KxCAED/Zm5YMCAALjZmWq6MiIiolMqBNz8/Hy4u5SedODs7Iz8/Xy1FEVHtdz+zAAv3RSMiJhUA4GVvjiWDg9C9ibOWKyMiIlKmcuDt2LEjFi5ciK1bt8LMrPQKTkFBARYtWoSOHTuqvUAiql1KZHJs+esmPo24jvxiGYwMJJjYrQGmPdcI5iaG2i6PiIioHJUD7+rVq9GnT59yN54wMzPD4cOH1V4gEdUe/9zJxOxdkYhJzgYAtPGuh2VDg9HE1VrLlREREVVO5cAbHByM+Ph4/Pjjj4iNjQUAjBw5Ei+//DLMzbmQPJE+yimU4uPDcdh67haEAGzNjTG7rz9GtPGCgYFE2+URERE9kUqBVyqVwt/fHwcOHMDEiRM1VRMR1RJCCByKSkH4/mikZhcBAIa29MDc/k3haGWq5eqIiIiqRqXAa2xsjMLCQk3VQkS1yN2MfCzcG42jsWkAAB8HCywdEowujRy1XBkREZFqVF6Hd8qUKVi5ciVKSko0UQ8RaVmJTI6NJxPR69OTOBqbBmNDCd56zg+H3unGsEtERDpJ5TG8Fy9exNGjR3HkyBEEBwfD0tJS6fldu3aprTgiqll/387AnN1RuPa/SWntfOyxfFgQ/Jw5KY2IiHSXyoHXzs4Ow4cP10QtRKQl2YVSfHQoDj+cL52UZmdhjDl9m+L51p6clEZERDpP5cC7efNmTdRBRFoghMDByBQs2h+NtJzSSWnDWnlgbr+mcOCkNCIi0hNVDrxyuRwfffQR9u3bh+LiYvTs2RMLFy7kUmREOurOo3ws2BuFY3EPAAC+jpZYNiQInfw4TpeIiPRLlQPvsmXLEB4ejtDQUJibm+Pzzz9HWloaNm3apMn6iEjNpDI5vj2dhNV/XEehVA4TQwNM6t4Qk7s3hJkx75RGRET6p8qBd+vWrVi7di3eeOMNAMAff/yB/v3745tvvoGBgcqLPRCRFly+lYG5uyMRm5IDAGjva49lQ4Ph52yl5cqIiIg0p8qB9/bt2+jXr5/icWhoKCQSCe7fvw9PT0+NFEdE6pFVIMVHh2Px4/nbEAKoZ2GMOf1KJ6VJJJyURkRE+q3KgbekpARmZmZKbcbGxpBKpWoviojUQwiBA/8mY/GBGDz436S051t7Yk6/prC3NNFydURERDWjyoFXCIGxY8fC1PT/Z24XFhZi0qRJSmvxch1eotrh9sN8zN8bhRPXSyelNXC0xNKhQejUkJPSiIiobqly4B0zZky5tldeeUWtxRDRs5PK5Nh4KhGf/xGPopLSSWlTevhhUvcGMDXipDQiIqp7qhx4uf4uUe13+dYjzNkVhbjU0klpHRs4YOnQIDR04qQ0IiKqu1S+8QQR1T5Z+VKsPByLbedvAyidlDavfwCGtfLgpDQiIqrzGHiJdJgQAvv+uY8lB2KQnlsMABjRxhOz+zZFPU5KIyIiAsDAS6Szbj3Mw7w9UTgVnw4AaOhkiWVDg9GhgYOWKyMiIqpdGHiJdExxSemktC+O/m9SmpEB3urhh9dDOCmNiIioIgy8RDrk4s1HmLMrEvFpuQCAzn4OWDokGL6Olk/Zk4iIqO5i4CXSAZn5xVh5KBY/XbgDAHCwNMG8AU0xpAUnpRERET0NAy9RLSaEwN6rpZPSHuaVTkp7qa0XPujrDzsLTkojIiKqCgZeolrqZnrppLTTCaWT0vycrbB8aDDa+dpruTIiIiLdwsBLVMsUl8ix4cQNfHksAcX/m5Q27Tk/vN6tIUyMDLRdHhERkc5h4CWqRS4kPcKc3ZFI+N+ktK6NHLFkcBB8OCmNiIio2hh4iWqBjLxifPh7LH6+VDopzdHKBPMHBGBQc3dOSiMiInpGtf770ZycHLzzzjvw9vaGubk5OnXqhIsXLz5xn6KiIsydOxfe3t4wNTWFj48PNm3aVEMVE1WdEAK7rtxFz09PKMLuyHb1cfTd7hjMFRiIiIjUotZf4Z0wYQKioqLw/fffw93dHT/88ANCQ0MRExMDDw+PCvcZMWIEUlNT8e2338LPzw/JycmQy+U1XDnRkyU+yMW8PVH468ZDAEBjl9JJaW18OCmNiIhInWp14C0oKMCvv/6KvXv3olu3bgCA8PBw7N+/H+vWrcPSpUvL7XPo0CGcOHECiYmJsLcvDQ4+Pj41WTbRExWVyLDhRCK++t+kNFMjA0zr2QgTuzbgpDQiIiINqNWBt6SkBDKZDGZmZkrt5ubmOH36dIX77Nu3D23atMGqVavw/fffw9LSEoMGDcKSJUtgbm5e4T5FRUUoKipSPM7OzgYASKVSSKVSNb2aypWdoybORZpR1T48n/QIC/bFIDE9HwDQ1c8B4QObor69BSBkkEplGq+VKsb3oe5jH+o29p/uq+k+VOU8EiGE0GAtz6xTp04wMTHBtm3b4OLigp9++gljxoyBn58f4uLiym3fp08fHD9+HKGhoViwYAHS09MxefJk9OjRA5s3b67wHOHh4Vi0aFG59m3btsHCwkLtr4nqnlwpsPeWAS48KL2Ca20sMMxHjpYOAhymS0REpLr8/HyMGjUKWVlZsLGxeeK2tT7w3rhxA+PHj8fJkydhaGiIVq1aoXHjxrh8+TKuXbtWbvvevXvj1KlTSElJga2tLQBg165deP7555GXl1fhVd6KrvB6eXkhPT39qT9AdZBKpYiIiECvXr1gbGys8fOR+lXWh0II7L56Hx8euo6M/NJ/iY5s64kZvRrBxpx9XZvwfaj72Ie6jf2n+2q6D7Ozs+Ho6FilwFurhzQAQMOGDXHixAnk5eUhOzsbbm5uePHFF9GgQYMKt3dzc4OHh4ci7AJA06ZNIYTA3bt30ahRo3L7mJqawtTUtFy7sbFxjb7pavp8pH6P9+GNB7mYuzsS5xIfAQCauFhj+bAgtPbmpLTajO9D3cc+1G3sP91XU32oyjlqfeAtY2lpCUtLS2RkZODw4cNYtWpVhdt17twZO3bsQG5uLqysrAAA169fh4GBATw9PWuyZKqjCqUyrDt+A+uO30CxTA4zYwO8E9oYr3XxhbEhJ6URERHVtFr/t+/hw4dx6NAhJCUlISIiAj169IC/vz/GjRsHAJg9ezZGjx6t2H7UqFFwcHDAuHHjEBMTg5MnT2LmzJkYP358pZPWiNTlXOIj9Pv8FD4/Go9imRwhjZ0QMT0Ek0IaMuwSERFpSa2/wpuVlYXZs2fj7t27sLe3x/Dhw7Fs2TLFZezk5GTcvn1bsb2VlRUiIiLw1ltvoU2bNnBwcMCIESMqXMKMSF0e5RXjxwQDXDh7CQDgZG2KhQMD0D/YjTePICIi0rJaH3hHjBiBESNGVPr8li1byrX5+/sjIiJCg1URlRJCYMflu1j+2zVkFhhAIgFeae+NGWFNYMtJaURERLVCrQ+8RLVVQlrppLTzSaWT0twtBD5/pT3aNnDScmVERET0OAZeIhUVSmVYeywB607cgFQmYG5siGnPNYRLZgxaeNlpuzwiIiL6DwZeIhWcSUjHvD1RSErPAwD0aOKExYOD4GptjIMHY7RcHREREVWEgZeoCh7mFmHZb9ew6+97AABna1OEDwpE3yBXSCQS3gqTiIioFmPgJXoCuVxgx+U7WH4wFlkFUkgkwOgO3ngvrAlszDgpjYiISBcw8BJVIj41B3N3R+HCzdJJaU3dbLBiWDDH6RIREekYBl6i/yiUyvDVnwnYcPL/J6W917sxxnbygRFvHkFERKRzGHiJHnMq/gHm7YnCrYf5AIDQps4IHxQIz3oWWq6MiIiIqouBlwjAg5wiLP0tBnuv3gcAuNqYIXxQIMICXXinNCIiIh3HwEt1mlwu8POlO1hx8BqyC0sgkQBjOvrgvd6NYc1JaURERHqBgZfqrOupOZizKxKXbmUAAALdbbB8aDCac1IaERGRXmHgpTqnUCrDF0fj8fXJRJTIBSxMDPFuL05KIyIi0lcMvFSnnLj+APP3ROH2o9JJab0CXLBoUCDc7cy1XBkRERFpCgMv1QlpOYVYcuAa9v9TOinNzbZsUpqrlisjIiIiTWPgJb0mlwv8dPE2Pvw9FjmFJTCQAGM7+eLd3o1hZcpffyIiorqAf+OT3opNycacXZG4cjsTABDsYYvlQ4MR7Gmr3cKIiIioRjHwkt4pKJbhiz/jsfF/k9IsTQwxI6wJRnf0gaEB19QlIiKqaxh4Sa8ci0vDgr1RuPOoAAAQFuiC8EGBcLPlpDQiIqK6ioGX9EJadiEWHYjBb/8mAwDcbc2waHAQegW4aLkyIiIi0jYGXtJpcrnAjxduY9XvscgpKp2UNr6zL6b3agxLTkojIiIiMPCSDruWnI05uyPx9/8mpTX3tMWyocEI8uCkNCIiIvp/DLykc/KLS/D5H/H45nQSZHIBK1MjzAxrglc6eHNSGhEREZXDwEs65VhsGubticK9zNJJaX2DXLFwYCBcbc20XBkRERHVVgy8pBNSswuxeH8MfossnZTmYWeOxYMD0bMpJ6URERHRkzHwUq0mkwv8eP4WPjoUh5yiEhgaSPBaF1+8E9oIFib89SUiIqKnY2KgWiv6fhbm7I7CP3cyAQDNveywfGgQAt05KY2IiIiqjoGXap28ohKs/uM6Np25CZlcwNrUCLP6NMGo9pyURkRERKpj4KVa5ei1VCzYG62YlNY/2A0LBgbAxYaT0oiIiKh6GHipVkjJKsSi/dH4PSoFQOmktKVDgtDD31nLlREREZGuY+AlrZLJBb4/exMfH7mO3P9NSpvQ1Rdv9+SkNCIiIlIPJgrSmqh7WZizOxL/3s0CALTwssOKYcFo6maj5cqIiIhInzDwUo3LKyrBpxHXsflMEuQCsDYzwvt9/DGqXX0YcFIaERERqRkDL9WoiJhULNwbhftZhQCAAc3csGBAAJw5KY2IiIg0hIGXasT9zAKE74vGkZhUAICXvTmWDA5C9yaclEZERESaxcBLGiWTC3z31018ciQOecUyGBlIMLFbA0x7rhHMTQy1XR4RERHVAQy8pDGRd0snpUXeK52U1tq7HpYNDYK/KyelERERUc1h4CW1yy0qwSdH4vDdXzchF4CNmRE+6NsUL7X14qQ0IiIiqnEMvKRWh6NTEL4vGsn/m5Q2qLk75g1oCmdrTkojIiIi7TDQdgFPk5OTg3feeQfe3t4wNzdHp06dcPHixSrte+bMGRgZGaFFixaaLZJwL7MAE767hDe+v4zkrELUt7fA1vHt8MXIlgy7REREpFW1/grvhAkTEBUVhe+//x7u7u744YcfEBoaipiYGHh4eFS6X2ZmJkaPHo2ePXsiNTW1BiuuW0pkcmz56yY+jbiO/P9NSnsjpAHeeq4RzIw5KY2IiIi0r1Zf4S0oKMCvv/6KVatWoVu3bvDz80N4eDj8/Pywbt26J+47adIkjBo1Ch07dqyhauuef+9mYvCaM1j62zXkF8vQxrseDr7dFTPD/Bl2iYiIqNao1Vd4S0pKIJPJYGam/JW4ubk5Tp8+Xel+mzdvRmJiIn744QcsXbr0qecpKipCUVGR4nF2djYAQCqVQiqVVrP6qis7R02cSx1yCkvw2dEE/Hj+tmJS2vthjfF8Kw8YGEh05nWok671IZXHPtR97EPdxv7TfTXdh6qcRyKEEBqs5Zl16tQJJiYm2LZtG1xcXPDTTz9hzJgx8PPzQ1xcXLnt4+Pj0aVLF5w6dQqNGzdGeHg49uzZg6tXr1Z6jvDwcCxatKhc+7Zt22BhYaHOl6PThAD+fSTBr0kGyJKWrrbQ2lGOId5y2JhouTgiIiKqU/Lz8zFq1ChkZWXBxubJS57W6iu8APD9999j/Pjx8PDwgKGhIVq1aoWRI0fi8uXL5baVyWQYNWoUFi1ahMaNG1f5HLNnz8a7776reJydnQ0vLy/07t37qT9AdZBKpYiIiECvXr1gbGys8fNVx73MAiw6cA3HrqcDALztLRA+sCm6+DloubLaQRf6kJ6Mfaj72Ie6jf2n+2q6D8u+ka+KWh94GzZsiBMnTiAvLw/Z2dlwc3PDiy++iAYNGpTbNicnB5cuXcLff/+NqVOnAgDkcjmEEDAyMsKRI0fw3HPPldvP1NQUpqam5dqNjY1r9E1X0+erihKZHJvPlE5KK5DKYGwowaSQhpjSw4/jdCtQG/uQVMM+1H3sQ93G/tN9NdWHqpyj1gfeMpaWlrC0tERGRgYOHz6MVatWldvGxsYGkZGRSm1r167Fn3/+iZ07d8LX17emytULV+9kYs6uSMQkl/4Lqp2PPZYPC4Kfs7WWKyMiIiKqulofeA8fPgwhBJo0aYKEhATMnDkT/v7+GDduHIDS4Qj37t3D1q1bYWBggKCgIKX9nZ2dYWZmVq6dKpddKMXHh+Pw/blbEAKwNTfG3H5N8XxrT94pjYiIiHROrQ+8WVlZmD17Nu7evQt7e3sMHz4cy5YtU1zGTk5Oxu3bt7VcpX4QQuD3qNI7paXllK5aMaylB+b0bwpHq/JDPoiIiIh0Qa0PvCNGjMCIESMqfX7Lli1P3D88PBzh4eHqLUoP3XmUj4X7ovFnbBoAwNfREsuGBKGTn6OWKyMiIiJ6NrU+8JJmSWVybDqdhNV/xCsmpb3Z3Q+TuzfkpDQiIiLSCwy8ddiV2xmYsysSsSk5AIB2vvZYPjQYfs5WWq6MiIiISH0YeOugrAIpPjocix/P34YQQD0LY8z536Q0iYST0oiIiEi/MPDWIUII/BaZjEX7Y/Dgf5PSnm/tiTn9msLekrdKIyIiIv3EwFtH3HmUj/l7o3A87gEAoIGjJZYODUKnhpyURkRERPqNgVfPSWVyfHMqCZ8fvY5CqRwmhgaY3KMh3uzeEKZGnJRGRERE+o+BV49dvlU6KS0utXRSWscGDlg6NAgNnTgpjYiIiOoOBl49lJUvxcrDsdh2vvSGHPUsjDGvfwCGtfLgpDQiIiKqcxh49YgQAvv/Tcbi/TFIzy2dlPZCa0/M5qQ0IiIiqsMYePXErYd5mLcnCqfi0wEADZ0ssWxoMDo0cNByZURERETaxcCr44pL5Nh4KhFfHI1HUYkcJkYGmNrDD2+ENOCkNCIiIiIw8Oq0SzcfYc7uSFxPzQUAdPZzwNIhwfB1tNRyZURERES1BwOvDsrKl+LDQ9fw04U7AAB7SxPMH9AUQ1pwUhoRERHRfzHw6hAhBPb9cx9LDsQgPbcYAPBiGy980Ncf9TgpjYiIiKhCDLw64mZ66aS00wmlk9L8nK2wfGgw2vnaa7kyIiIiotqNgbeWKy6R4+uTN/DFnwko/t+ktGnP+eH1bg1hYmSg7fKIiIiIaj0G3lrsQlLppLSEtNJJaV0bOWLJ4CD4cFIaERERUZUx8NZCmfnFWHEwFj9fKp2U5mhlgvkDAjCouTsnpRERERGpiIG3FhFCYPffd7H0wDU8zCudlDaynRfe7+MPOwtOSiMiIiKqDgbeWiKtABi75TL+SnwEAGjsUjoprY0PJ6URERERPQsGXi2TyQW+OnYDa/4xRIl4BFMjA0zr2QgTuzbgpDQiIiIiNWDg1TIDCXD5diZKhARd/BywbGgwvB04KY2IiIhIXXgJUcskEgnCBzbF6EYybBrdimGXiIiISM0YeGsBb3sLtHYUXIGBiIiISAMYeImIiIhIrzHwEhEREZFeY+AlIiIiIr3GwEtEREREeo2Bl4iIiIj0GgMvEREREek1Bl4iIiIi0msMvERERESk1xh4iYiIiEivMfASERERkV5j4CUiIiIivcbAS0RERER6jYGXiIiIiPQaAy8RERER6TUjbRdQGwkhAADZ2dk1cj6pVIr8/HxkZ2fD2Ni4Rs5J6sU+1H3sQ93HPtRt7D/dV9N9WJbTynLbkzDwViAnJwcA4OXlpeVKiIiIiOhJcnJyYGtr+8RtJKIqsbiOkcvluH//PqytrSGRSDR+vuzsbHh5eeHOnTuwsbHR+PlI/diHuo99qPvYh7qN/af7aroPhRDIycmBu7s7DAyePEqXV3grYGBgAE9Pzxo/r42NDd/kOo59qPvYh7qPfajb2H+6ryb78GlXdstw0hoRERER6TUGXiIiIiLSawy8tYCpqSkWLlwIU1NTbZdC1cQ+1H3sQ93HPtRt7D/dV5v7kJPWiIiIiEiv8QovEREREek1Bl4iIiIi0msMvERERESk1xh4iYiIiEivMfBqyJo1a+Dj4wMzMzO0b98eFy5ceOL2O3bsgL+/P8zMzBAcHIyDBw8qPS+EwIIFC+Dm5gZzc3OEhoYiPj5eky+hzlNnH0qlUrz//vsIDg6GpaUl3N3dMXr0aNy/f1/TL6POUvd78HGTJk2CRCLB6tWr1Vw1PU4TfXjt2jUMGjQItra2sLS0RNu2bXH79m1NvYQ6T919mJubi6lTp8LT0xPm5uYICAjA+vXrNfkS6jRV+i86OhrDhw+Hj4/PEz8fVf2dUBtBard9+3ZhYmIiNm3aJKKjo8XEiROFnZ2dSE1NrXD7M2fOCENDQ7Fq1SoRExMj5s2bJ4yNjUVkZKRimw8//FDY2tqKPXv2iH/++UcMGjRI+Pr6ioKCgpp6WXWKuvswMzNThIaGip9//lnExsaKs2fPinbt2onWrVvX5MuqMzTxHiyza9cu0bx5c+Hu7i4+++wzDb+SuksTfZiQkCDs7e3FzJkzxZUrV0RCQoLYu3dvpcekZ6OJPpw4caJo2LChOHbsmEhKShIbNmwQhoaGYu/evTX1suoMVfvvwoULYsaMGeKnn34Srq6uFX4+qnpMdWLg1YB27dqJKVOmKB7LZDLh7u4uVqxYUeH2I0aMEP3791dqa9++vXjjjTeEEELI5XLh6uoqPvroI8XzmZmZwtTUVPz0008aeAWk7j6syIULFwQAcevWLfUUTQqa6r+7d+8KDw8PERUVJby9vRl4NUgTffjiiy+KV155RTMFUzma6MPAwECxePFipW1atWol5s6dq8bKSQjV++9xlX0+PssxnxWHNKhZcXExLl++jNDQUEWbgYEBQkNDcfbs2Qr3OXv2rNL2ABAWFqbYPikpCSkpKUrb2Nraon379pUek6pPE31YkaysLEgkEtjZ2amlbiqlqf6Ty+V49dVXMXPmTAQGBmqmeAKgmT6Uy+X47bff0LhxY4SFhcHZ2Rnt27fHnj17NPY66jJNvQ87deqEffv24d69exBC4NixY7h+/Tp69+6tmRdSR1Wn/7RxTFUw8KpZeno6ZDIZXFxclNpdXFyQkpJS4T4pKSlP3L7sv6ock6pPE334X4WFhXj//fcxcuRI2NjYqKdwAqC5/lu5ciWMjIwwbdo09RdNSjTRh2lpacjNzcWHH36IPn364MiRIxg6dCiGDRuGEydOaOaF1GGaeh9++eWXCAgIgKenJ0xMTNCnTx+sWbMG3bp1U/+LqMOq03/aOKYqjDR+BiJSIpVKMWLECAghsG7dOm2XQ1Vw+fJlfP7557hy5QokEom2y6FqkMvlAIDBgwdj+vTpAIAWLVrgr7/+wvr16xESEqLN8qiKvvzyS5w7dw779u2Dt7c3Tp48iSlTpsDd3b3c1WGix/EKr5o5OjrC0NAQqampSu2pqalwdXWtcB9XV9cnbl/2X1WOSdWniT4sUxZ2b926hYiICF7d1QBN9N+pU6eQlpaG+vXrw8jICEZGRrh16xbee+89+Pj4aOR11GWa6ENHR0cYGRkhICBAaZumTZtylQYN0EQfFhQUYM6cOfj0008xcOBANGvWDFOnTsWLL76Ijz/+WDMvpI6qTv9p45iqYOBVMxMTE7Ru3RpHjx5VtMnlchw9ehQdO3ascJ+OHTsqbQ8AERERiu19fX3h6uqqtE12djbOnz9f6TGp+jTRh8D/h934+Hj88ccfcHBw0MwLqOM00X+vvvoq/v33X1y9elXxx93dHTNnzsThw4c192LqKE30oYmJCdq2bYu4uDilba5fvw5vb281vwLSRB9KpVJIpVIYGChHF0NDQ8UVfFKP6vSfNo6pEo1Pi6uDtm/fLkxNTcWWLVtETEyMeP3114WdnZ1ISUkRQgjx6quvig8++ECx/ZkzZ4SRkZH4+OOPxbVr18TChQsrXJbMzs5O7N27V/z7779i8ODBXJZMg9Tdh8XFxWLQoEHC09NTXL16VSQnJyv+FBUVaeU16jNNvAf/i6s0aJYm+nDXrl3C2NhYfP311yI+Pl58+eWXwtDQUJw6darGX19doIk+DAkJEYGBgeLYsWMiMTFRbN68WZiZmYm1a9fW+OvTd6r2X1FRkfj777/F33//Ldzc3MSMGTPE33//LeLj46t8TE1i4NWQL7/8UtSvX1+YmJiIdu3aiXPnzimeCwkJEWPGjFHa/pdffhGNGzcWJiYmIjAwUPz2229Kz8vlcjF//nzh4uIiTE1NRc+ePUVcXFxNvJQ6S519mJSUJABU+OfYsWM19IrqFnW/B/+LgVfzNNGH3377rfDz8xNmZmaiefPmYs+ePZp+GXWauvswOTlZjB07Vri7uwszMzPRpEkT8cknnwi5XF4TL6fOUaX/Kvt7LiQkpMrH1CSJEEJo/joyEREREZF2cAwvEREREek1Bl4iIiIi0msMvERERESk1xh4iYiIiEivMfASERERkV5j4CUiIiIivcbAS0RERER6jYGXiIiIiPQaAy8RkY44fvw4JBIJMjMz1bJteHg4WrRoUa7NxcUFEokEe/bsUam+o0ePomnTppDJZE/dNj09Hc7Ozrh7965K5yAiqg4GXiIiNRs7diwkEgkkEgmMjY3h6+uLWbNmobCwUNulKZkxYwaOHj2qeHzt2jUsWrQIGzZsQHJyMvr27QsfHx+sXr26SsebNWsW5s2bB0NDw6du6+joiNGjR2PhwoXVLZ+IqMoYeImINKBPnz5ITk5GYmIiPvvsM2zYsKHWhTsrKys4ODgoHt+4cQMAMHjwYLi6usLU1LTKxzp9+jRu3LiB4cOHV3mfcePG4ccff8SjR4+qXjQRUTUw8BIRaYCpqSlcXV3h5eWFIUOGIDQ0FBEREYrn5XI5VqxYAV9fX5ibm6N58+bYuXOn0jEOHjyIxo0bw9zcHD169MDNmzeVnr916xYGDhyIevXqwdLSEoGBgTh48KDSNpcvX0abNm1gYWGBTp06IS4uTvHc40MawsPDMXDgQACAgYEBJBIJunfvjlu3bmH69OmKK9aV2b59O3r16gUzMzNFm4+Pj2K/x/+UCQwMhLu7O3bv3l21HyoRUTUx8BIRaVhUVBT++usvmJiYKNpWrFiBrVu3Yv369YiOjsb06dPxyiuv4MSJEwCAO3fuYNiwYRg4cCCuXr2KCRMm4IMPPlA67pQpU1BUVISTJ08iMjISK1euhJWVldI2c+fOxSeffIJLly7ByMgI48ePr7DGGTNmYPPmzQCA5ORkJCcnY9euXfD09MTixYsVbZU5deoU2rRpo9R28eJFxX53795Fhw4d0LVrV6Vt2rVrh1OnTj3lJ0hE9GyMtF0AEZE+OnDgAKysrFBSUoKioiIYGBjgq6++AgAUFRVh+fLl+OOPP9CxY0cAQIMGDXD69Gls2LABISEhWLduHRo2bIhPPvkEANCkSRNFqC1z+/ZtDB8+HMHBwYpj/NeyZcsQEhICAPjggw/Qv39/FBYWKl2JBUqHN9jZ2QEAXF1dFe2GhoawtrZWaqvIrVu34O7urtTm5OSk+P+3334bycnJuHjxotI27u7u+Pvvv594bCKiZ8XAS0SkAT169MC6deuQl5eHzz77DEZGRorxrQkJCcjPz0evXr2U9ikuLkbLli0BlE4ga9++vdLzZeG4zLRp0/Dmm2/iyJEjCA0NxfDhw9GsWTOlbR5/7ObmBgBIS0tD/fr11fNC/6egoKBciC7z9ddf49tvv8Vff/2lFIIBwNzcHPn5+WqthYjovzikgYhIAywtLeHn54fmzZtj06ZNOH/+PL799lsAQG5uLgDgt99+w9WrVxV/YmJiyo3jfZIJEyYgMTERr776KiIjI9GmTRt8+eWXStsYGxsr/r9s/KxcLn/Wl1eOo6MjMjIyyrUfO3YMb731FrZu3VoujAPAo0ePyoVgIiJ1Y+AlItIwAwMDzJkzB/PmzUNBQQECAgJgamqK27dvw8/PT+mPl5cXAKBp06a4cOGC0nHOnTtX7theXl6YNGkSdu3ahffeew8bN25Ua+0mJiZVWle3ZcuWiImJUWpLSEjA888/jzlz5mDYsGEV7hcVFaW4qk1EpCkMvERENeCFF16AoaEh1qxZA2tra8yYMQPTp0/Hd999hxs3buDKlSv48ssv8d133wEAJk2ahPj4eMycORNxcXHYtm0btmzZonTMd955B4cPH0ZSUhKuXLmCY8eOoWnTpmqt28fHBydPnsS9e/eQnp5e6XZhYWE4ffq04nFBQQEGDhyIli1b4vXXX0dKSoriT5n8/HxcvnwZvXv3VmvNRET/xcBLRFQDjIyMMHXqVKxatQp5eXlYsmQJ5s+fjxUrVqBp06bo06cPfvvtN/j6+gIA6tevj19//RV79uxB8+bNsX79eixfvlzpmDKZDFOmTFHs37hxY6xdu1atdS9evBg3b95Ew4YNnzj04OWXX0Z0dLRi2bPU1FTExsbi6NGjcHd3h5ubm+JPmb1796J+/frlVm4gIlI3iRBCaLsIIiLSfTNnzkR2djY2bNhQpe07dOiAadOmYdSoURqujIjqOl7hJSIitZg7dy68vb2rNCkuPT0dw4YNw8iRI2ugMiKq63iFl4iIiIj0Gq/wEhEREZFeY+AlIiIiIr3GwEtEREREeo2Bl4iIiIj0GgMvEREREek1Bl4iIiIi0msMvERERESk1xh4iYiIiEivMfASERERkV77P5QpdI5Mb9ICAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 800x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Schechter luminosity function parameters\n",
    "M_star = -20.83  # SDSS Mr*\n",
    "alpha = -1.0  # Faint-end slope\n",
    "phi_star = 0.01  # Normalization (Mpc^-3)\n",
    "\n",
    "# Constants\n",
    "H0 = 70  # Hubble constant (km/s/Mpc)\n",
    "c = 3e5  # Speed of light (km/s)\n",
    "\n",
    "# Function to calculate the Schechter luminosity function\n",
    "def schechter_function(M, M_star, alpha, phi_star):\n",
    "    term = 10**(0.4 * (M_star - M))\n",
    "    return 0.4 * np.log(10) * phi_star * term**(1 + alpha) * np.exp(-term)\n",
    "\n",
    "# Function to calculate comoving volume element (dV/dz in Mpc^3 per steradian)\n",
    "def comoving_volume_element(z, H0, Omega_m=0.3, Omega_L=0.7):\n",
    "    \"\"\"Comoving volume element per unit redshift per steradian.\"\"\"\n",
    "    Ez = np.sqrt(Omega_m * (1 + z)**3 + Omega_L)\n",
    "    dV_dz = c / H0 * (1 + z)**2 / Ez\n",
    "    return dV_dz\n",
    "\n",
    "# 1. Number of galaxies with r < 20.5 mag\n",
    "def number_density_r_limit(m_limit, M_star, alpha, phi_star, z_max):\n",
    "    \"\"\"Calculate the number density of galaxies with apparent magnitude limit.\"\"\"\n",
    "    def absolute_magnitude_limit(m, z):\n",
    "        \"Convert apparent magnitude to absolute magnitude at redshift z.\"\n",
    "        d_l = c * z / H0 * (1 + z / 2)  # Luminosity distance approximation (Mpc)\n",
    "        return m - 5 * np.log10(d_l * 1e6) + 5\n",
    "\n",
    "    z_values = np.linspace(0, z_max, 1000)\n",
    "    m_abs_limits = absolute_magnitude_limit(m_limit, z_values)\n",
    "    phi_values = schechter_function(m_abs_limits, M_star, alpha, phi_star)\n",
    "    dV_dz_values = comoving_volume_element(z_values, H0)\n",
    "    \n",
    "    number_density = np.trapz(phi_values * dV_dz_values, z_values)\n",
    "    return number_density\n",
    "\n",
    "# 2. Ratio of galaxy counts at M_r = -16 mag to M_r = -22 mag\n",
    "def galaxy_count_ratio(M1, M2, M_star, alpha):\n",
    "    phi1 = schechter_function(M1, M_star, alpha, phi_star)\n",
    "    phi2 = schechter_function(M2, M_star, alpha, phi_star)\n",
    "    return phi1 / phi2\n",
    "\n",
    "# 3. Redshift distribution for galaxies with M_r = -20 mag\n",
    "def redshift_distribution(M, M_star, alpha, phi_star, z_max):\n",
    "    z_values = np.linspace(0, z_max, 1000)\n",
    "    phi_values = schechter_function(M, M_star, alpha, phi_star)\n",
    "    dV_dz_values = comoving_volume_element(z_values, H0)\n",
    "    \n",
    "    P_z = phi_values * dV_dz_values\n",
    "    P_z /= np.trapz(P_z, z_values)  # Normalize to make it a probability distribution\n",
    "    return z_values, P_z\n",
    "\n",
    "# Inputs\n",
    "m_limit = 20.5\n",
    "z_max = 0.1  # Approximation for shallow surveys\n",
    "M1, M2 = -16, -22\n",
    "M_target = -20\n",
    "\n",
    "# 1. Number density of galaxies\n",
    "num_density = number_density_r_limit(m_limit, M_star, alpha, phi_star, z_max)\n",
    "print(f\"Number density of galaxies with r < {m_limit} mag: {num_density:.2e} per steradian\")\n",
    "\n",
    "# 2. Ratio of galaxy counts\n",
    "ratio = galaxy_count_ratio(M1, M2, M_star, alpha)\n",
    "print(f\"Ratio of galaxy counts (M_r = {M1} / M_r = {M2}): {ratio:.2f}\")\n",
    "\n",
    "# 3. Redshift distribution\n",
    "z_values, P_z = redshift_distribution(M_target, M_star, alpha, phi_star, z_max)\n",
    "\n",
    "# Plot redshift distribution\n",
    "plt.figure(figsize=(8, 5))\n",
    "plt.plot(z_values, P_z, label=f\"M_r = {M_target}\")\n",
    "plt.xlabel(\"Redshift (z)\")\n",
    "plt.ylabel(\"Probability Density P(z)\")\n",
    "plt.title(\"Redshift Distribution for M_r = -20 mag Galaxies\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3645d936",
   "metadata": {},
   "source": [
    "作业2:某旋涡星系和椭圆星系，其观测到的I波段的视星等都是18等；半光度半径内的表面亮度都是20mag/arcsec^2；其中旋涡星系（图像椭率为 0.5 ）观测得到的HI的速度展宽是200km/s，椭圆星系观测到的中心速度弥散度是200km/s，请问这两个星系的距离分别是多少？ • 旋涡星系的 Tully-Fisher关系以及椭圆星系的Fundamental Plane的内秉弥散在0.2个星等左右，请问以上距离计算的误差是多少？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2d1a290b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Spiral Galaxy:\n",
      "Distance: 6.34e+09 pc\n",
      "Uncertainty: 5.84e+08 pc\n",
      "\n",
      "Elliptical Galaxy:\n",
      "Distance: 1.82e-08 pc\n",
      "Uncertainty: 1.68e-09 pc\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def calculate_distance_tully_fisher(m_I, W, a=-10, b=-3, intrinsic_scatter=0.2):\n",
    "    # Tully-Fisher 关系计算绝对星等\n",
    "    M_I = a * np.log10(W) + b\n",
    "\n",
    "    # 距离模数计算\n",
    "    mu = m_I - M_I\n",
    "\n",
    "    # 距离（单位：pc）\n",
    "    D = 10 ** ((mu + 5) / 5)\n",
    "\n",
    "    # 距离误差估计\n",
    "    delta_mu = intrinsic_scatter\n",
    "    delta_D = D * np.log(10) * delta_mu / 5\n",
    "\n",
    "    return D, delta_D\n",
    "\n",
    "def calculate_distance_fundamental_plane(m_I, sigma, surface_brightness, \n",
    "                                         a=1.2, b=-0.8, c=2.5, intrinsic_scatter=0.2):\n",
    "    # 假设 Fundamental Plane 关系估算星系半径\n",
    "    log_R_e = a * np.log10(sigma) + b * surface_brightness + c\n",
    "    R_e = 10 ** log_R_e  # 单位：kpc\n",
    "\n",
    "    # 假设使用半光度半径对应的绝对星等估计\n",
    "    M_I = m_I - 5 * np.log10(R_e * 1e3) + 5  # 半光度半径与距离模数关系\n",
    "\n",
    "    # 距离模数计算\n",
    "    mu = m_I - M_I\n",
    "\n",
    "    # 距离（单位：pc）\n",
    "    D = 10 ** ((mu + 5) / 5)\n",
    "\n",
    "    # 距离误差估计\n",
    "    delta_mu = intrinsic_scatter\n",
    "    delta_D = D * np.log(10) * delta_mu / 5\n",
    "\n",
    "    return D, delta_D\n",
    "\n",
    "# 输入参数\n",
    "m_I_spiral = 18  # 旋涡星系视星等\n",
    "W_spiral = 200   # HI 速度展宽\n",
    "\n",
    "m_I_elliptical = 18  # 椭圆星系视星等\n",
    "sigma_elliptical = 200  # 中心速度弥散\n",
    "surface_brightness = 20  # 表面亮度（mag/arcsec^2）\n",
    "\n",
    "# 计算旋涡星系的距离和误差\n",
    "D_spiral, delta_D_spiral = calculate_distance_tully_fisher(m_I_spiral, W_spiral)\n",
    "\n",
    "# 计算椭圆星系的距离和误差\n",
    "D_elliptical, delta_D_elliptical = calculate_distance_fundamental_plane(\n",
    "    m_I_elliptical, sigma_elliptical, surface_brightness\n",
    ")\n",
    "\n",
    "# 输出结果\n",
    "print(\"Spiral Galaxy:\")\n",
    "print(f\"Distance: {D_spiral:.2e} pc\")\n",
    "print(f\"Uncertainty: {delta_D_spiral:.2e} pc\")\n",
    "\n",
    "print(\"\\nElliptical Galaxy:\")\n",
    "print(f\"Distance: {D_elliptical:.2e} pc\")\n",
    "print(f\"Uncertainty: {delta_D_elliptical:.2e} pc\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09d6c9e6",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}