{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2993ee30",
   "metadata": {},
   "source": [
    "作业1:• 单位平方度内，有多少r<20.5 mag 的星系？\n",
    "• SDSS给出的本地星系的Schechter光度函数的参数为Mr*=-20.83，α=-1.0，请问单位体积范围内Mr=-16 mag的星系数目是Mr=-22mag星系数目的多少倍？SDSS是星等极限样本（r<17.77）,请问在SDSS的观测样本中，Mr=-22 mag的星系数目是Mr=-16mag星系数目的多少倍？\n",
    "• 根据以上光度函数，画出一个Mr=20 mag星系的红移分布概率不考虑星系的K改正和演化效应"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ff98f236",
   "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": [
      "/var/folders/mz/pg3ztpgj1f10l242jx2qkjz80000gn/T/ipykernel_27431/3770505998.py:31: RuntimeWarning: divide by zero encountered in log10\n",
      "  return m - 5 * np.log10(d_l * 1e6) + 5\n",
      "/var/folders/mz/pg3ztpgj1f10l242jx2qkjz80000gn/T/ipykernel_27431/3770505998.py:38: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.\n",
      "  number_density = np.trapz(phi_values * dV_dz_values, z_values)\n",
      "/var/folders/mz/pg3ztpgj1f10l242jx2qkjz80000gn/T/ipykernel_27431/3770505998.py:54: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.\n",
      "  P_z /= np.trapz(P_z, z_values)  # Normalize to make it a probability distribution\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAHUCAYAAADY9fvpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACA/klEQVR4nO3dd1RU1xYG8G8YhqEroFQBERFEUFETYwUb9hqjiTH2FDWxJSaxY8PEGGPyLHkxxhKjKcaWWDH22AuJFAsKYkOs9DLlvD8I84KAMjjDMMP3W4u1nHPPvXdfNgPbO+ecKxFCCBARERERmSgzQwdARERERKRPLHiJiIiIyKSx4CUiIiIik8aCl4iIiIhMGgteIiIiIjJpLHiJiIiIyKSx4CUiIiIik8aCl4iIiIhMGgteIiIiIjJpLHjJZKxZswYSiUTzZW5uDjc3N7z66qu4cuWKXs6VlJRU5n2SkpIgkUiwaNGiZ/aNiIiARCIp0pafn4933nkHbm5ukEqlaNy4MW7fvo2IiAhER0eXKYaDBw8W+R5ZWFigZs2aaNWqFaZNm4br168X26c81woAkZGR2Lp1q1b7lHSusLAwBAUFaXWcZ9m5cyciIiJK3Fa7dm0MGzZMp+crqz/++APNmjWDjY0NJBKJ1t8/bRT+PEokklK/FyNGjND0MXZ37tzB9OnT0aJFC9SoUQP29vZo2rQpvvnmG6hUqmL9MzMzMWHCBLi7u8PS0hKNGzfGjz/+aIDI9S8xMRHjxo1D/fr1YWNjA0tLS9SuXRuDBw/GgQMHUJ4Hshb+fK1Zs0b3Af+j8PfZwYMH9XYOMh0seMnkrF69GsePH8e+ffvw7rvvYvv27WjdujUePXpk6NDKbNSoUTh+/HiRthUrVuC///0vpk2bhqNHj+L777/H7du3MXv27DIXvIUiIyNx/PhxHDhwAKtWrUJYWBi+++471K9fHz/88EORvt27d8fx48fh5uam9Tm0LdjKey5t7dy5E7Nnzy5x25YtWzBjxgy9nr8kQggMGDAAMpkM27dvx/HjxxEaGqr389rZ2WHNmjVQq9VF2jMzM/HLL7/A3t5e7zFUhLNnz2LdunXo0KED1q1bh19//RWhoaEYPXo03nzzzWL9+/Xrh7Vr12LWrFnYtWsXXnjhBbz22mvYsGGDAaLXn+3btyM4OBjbt2/H0KFDsWXLFuzZswczZszAgwcP0L59e+zfv9/QYZaoSZMmOH78OJo0aWLoUMgYCCITsXr1agFAnD59ukj77NmzBQDx3Xff6fxciYmJZd4nMTFRABCfffZZuc45atQoYWVlVaTt9OnTAoBYvXp1mY5x4MABAUD88ssvxbY9ePBAhISECHNzc/H333+XK8Z/s7GxEUOHDi1T3+zsbKFWq0vcFhoaKho0aPDc8fzb2LFjRWX79Xfz5k0BQHz66ac6O+bTvq+FP4+jRo0SAMTevXuLbP/222+FlZWVGDx4sE6+V1lZWc99jOfx8OFDkZ+fX6y98GchOTlZ07Zjxw4BQGzYsKFI306dOgl3d3ehVCr1Hm9FSEhIENbW1uKFF14QaWlpJfY5cOCAiI6O1vrYhT9fZf3dRKRvvMNLJq9Zs2YAgLt37xZpP3PmDHr16gVHR0dYWloiJCQEP//8c7H9T5w4gVatWsHS0hLu7u6YMmUKFApFsX779+9HWFgYnJycYGVlBS8vL7z88svIzs4u1nfx4sXw8fGBra0tWrRogRMnThTZ/uSQBolEgm+//RY5OTmaj5jXrFmDF154AQAwfPjwZ348/SyOjo7473//C6VSiS+++ELTXtIwg/Pnz6NHjx5wdnaGXC6Hu7s7unfvjps3b2rizcrKwtq1azVxhYWFFTne3r17MWLECNSsWRPW1tbIy8t76vCJI0eO4KWXXoKVlRU8PDwwY8aMIh9Fl/bx5pMfrQ4bNgzLli3TxFn4VXjOkoY0JCcnY/DgwZrrrV+/Pj7//PMid0X/PWTlWfl9UkREBGrVqgUA+OijjyCRSFC7dm3N9qNHj6JDhw6ws7ODtbU1WrZsiR07dhQ5xtO+r0/j7++Pli1b4rvvvivS/t1336Ffv36oVq3aU/cv7XokEgnOnTuH/v37w8HBAb6+vlofR5ccHBwgk8mKtb/44osAoPnZBQru8tva2uKVV14p0nf48OG4ffs2Tp48+dRzDRs2DLa2trh48SI6d+4MGxsbuLm54ZNPPgFQ8DuldevWsLGxQb169bB27doi+9+7dw9jxoxBYGAgbG1t4ezsjPbt2+PIkSPFznXz5k30798fdnZ2qF69Ol5//XWcPn26TMMJFi9ejOzsbCxfvrzUO/lhYWFo1KiR5nVCQgKGDx8OPz8/WFtbw8PDAz179sSFCxeeei5t9n3nnXdgaWmJs2fPatrUajU6dOgAFxcX3LlzB0Dp7/my/G7Pzs7GBx98AB8fH1haWsLR0RHNmjXDxo0bn3kdZJxY8JLJS0xMBADUq1dP03bgwAG0atUKjx8/xtdff41t27ahcePGGDhwYJE/EnFxcejQoQMeP36MNWvW4Ouvv8b58+cxb968IudISkpC9+7dYWFhge+++w67d+/GJ598AhsbG+Tn5xfpu2zZMkRFRWHJkiX44YcfkJWVhW7duiEtLa3Uazh+/Di6desGKysrHD9+HMePH0e7du2wevVqAMD06dM17aNGjSr39+qFF16Am5sbDh8+XGqfrKwsdOrUCXfv3i1yLV5eXsjIyNDEa2VlhW7dumniWr58eZHjjBgxAjKZDN9//z02bdpUYjFSKCUlBa+++ipef/11bNu2Df3798e8efMwfvx4ra9xxowZ6N+/vybOwq/ShlHcu3cPLVu2xN69ezF37lxs374dHTt2xAcffIB33323WP/y5HfUqFHYvHkzAOC9997D8ePHsWXLFgDAoUOH0L59e6SlpWHVqlXYuHEj7Ozs0LNnT/z000/FjqXN97XQyJEjsXXrVs2wn0uXLuHYsWMYOXLkM/d9mn79+qFu3br45Zdf8PXXXz+1r1KpLNOXKMd40qfZv38/zM3Ni/x+iImJQf369WFubl6kb8OGDTXbn0WhUKBfv37o3r07tm3bhq5du2LKlCmYOnUqhg4dihEjRmDLli3w9/fHsGHDihR3Dx8+BADMmjULO3bswOrVq1GnTh2EhYUVKe6ysrLQrl07HDhwAJ9++il+/vlnuLi4YODAgWW69qioKLi5uWluCpTF7du34eTkhE8++QS7d+/GsmXLYG5ujubNm+PSpUs62XfJkiWoX78+BgwYgMePHwMAZs+ejYMHD2L9+vVPHfJU1t/tkyZNwooVKzBu3Djs3r0b33//PV555RU8ePCgzN8LMjKGvsVMpCuFwwxOnDghFAqFyMjIELt37xaurq6ibdu2QqFQaPoGBASIkJCQIm1CCNGjRw/h5uYmVCqVEEKIgQMHCisrK5GSkqLpo1QqRUBAQJEhDZs2bRIAnvrRX+FHfMHBwUU+Ej116pQAIDZu3KhpmzVrVrGPkYcOHSpsbGyKtOlySEOh5s2bFxk68eTwjTNnzggAYuvWrU89V2lDGgqPN2TIkFK3/XuoSGhoqAAgtm3bVqTvm2++KczMzMT169eLXNuBAweK9Cvpo9WnDWnw9vYuEvfHH38sAIiTJ08W6Td69GghkUjEpUuXipynLPktSWlDXl566SXh7OwsMjIyNG1KpVIEBQWJWrVqaYYsPO37+qzzZWRkCFtbW7F06VIhhBCTJ08WPj4+Qq1Wl2v4R+HP78yZM7WKpSxfT+b3eezZs0eYmZmJiRMnFmn38/MTnTt3Ltb/9u3bAoCIjIx86nGHDh0qAIhff/1V06ZQKETNmjUFAHHu3DlN+4MHD4RUKhWTJk0q9XhKpVIoFArRoUMH0bdvX037smXLBACxa9euIv3ffvvtMv1esLS0FC+99FKxdpVKJRQKhear8PdhabHl5+cLPz+/It/HsgxpKG1fIYS4cuWKsLe3F3369BH79u0TZmZmYvr06UX6lPSeL+vv9qCgINGnT59SYyPTwzu8ZHJeeuklyGQy2NnZoUuXLnBwcMC2bds0d2sSEhJw8eJFvP766wCK3lnq1q0b7ty5o7nbcODAAc3HaIWkUmmxOyiNGzeGhYUF3nrrLaxduxbXrl0rNb7u3btDKpVqXhfeNSpphQRDEM+4g1a3bl04ODjgo48+wtdff424uLhynefll18uc187Ozv06tWrSNugQYOgVqufejdaF/bv34/AwEDNR9+Fhg0bBiFEsQk9usxvVlYWTp48if79+8PW1lbTLpVK8cYbb+DmzZvF7qpp830tVPjx/XfffQelUol169Zphsk8j7LG4u7ujtOnT5fpq2nTpk89lkqlKvKefnIyXqFz585hwIABeOmll7BgwYJi25927WX5vkgkEnTr1k3z2tzcHHXr1oWbmxtCQkI07Y6OjnB2di728/H111+jSZMmsLS0hLm5OWQyGf744w/Ex8dr+hw6dEjze+7fXnvttWfG9zT9+vWDTCbTfI0bN06zTalUIjIyEoGBgbCwsIC5uTksLCxw5cqVIrGVRJt969ati5UrV2Lr1q3o0aMH2rRp88zhWtr8bn/xxRexa9cufPzxxzh48CBycnLK8Z0iY8KCl0zOunXrcPr0aezfvx9vv/024uPji/wBKBzL+8EHHxT5pS6TyTBmzBgAwP379wEADx48gKura7FzPNnm6+uLffv2wdnZGWPHjoWvry98fX3x5ZdfFtvXycmpyGu5XA4AleYXbnJyMtzd3UvdXq1aNRw6dAiNGzfG1KlT0aBBA7i7u2PWrFkljm0ujTYrMfz7PxyFCnOg748gHzx4UGKshd+jJ8+vy/w+evQIQgitzl/eFS5GjhyJc+fOYf78+bh3755OlmYraywWFhZo3Lhxmb7+XfiXxNfXt8h7es6cOcX6nD9/Hp06dYKfnx927typyVEhJyenEn+uCocaODo6PvOarK2tYWlpWew6S9rXwsICubm5mteLFy/G6NGj0bx5c/z66684ceIETp8+jS5duhT5OXrw4EGJ742S2kri5eVV4n/EPv/8c81/MJ40adIkzJgxA3369MFvv/2GkydP4vTp02jUqNEzf8a13bd79+5wcXFBbm4uJk2aVOQ/kiXR5nf7V199hY8++ghbt25Fu3bt4OjoiD59+uh8CUuqPMyf3YXIuNSvX18zJq1du3ZQqVT49ttvsWnTJvTv3x81atQAAEyZMgX9+vUr8Rj+/v4ACv7wpaSkFNteUlubNm3Qpk0bqFQqnDlzBv/5z38wYcIEuLi44NVXX9XV5enVqVOnkJKS8syxm8HBwfjxxx8hhMDff/+NNWvWYM6cObCyssLHH39cpnNpc/fwyQmHwP9zUFhgFhYXT07SKvwDV15OTk6aSTL/dvv2bQDQ/Dzpg4ODA8zMzLQ6f3nvyrZq1Qr+/v6YM2cOOnXqBE9Pz3IdpzyxJCUlwcfHp0x9Dxw4oJkAWZLffvutyM/Ak/95O3/+PDp27Ahvb2/s3bu3xEl5wcHB2LhxI5RKZZFxvIWTq3S9LvST1q9fj7CwMKxYsaJIe+EY+UJOTk44depUsf1L+v1Ukk6dOmHZsmU4c+ZMkXG8T5tguH79egwZMgSRkZFF2u/fv4/q1as/9Xza7vvOO+8gIyMDDRo0wLhx49CmTRs4ODiUenxtfrfb2Nhg9uzZmD17Nu7evau529uzZ09cvHjxqddBxol3eMnkLVy4EA4ODpg5cybUajX8/f3h5+eHv/76C82aNSvxy87ODkBBwfzHH38UKbhUKlWJk4UKSaVSNG/eXLMSwLlz5/R2bbq8O/zw4UO88847kMlkmDhxYpn2kUgkaNSoEb744gtUr169yLXK5XKd3bXOyMjA9u3bi7Rt2LABZmZmaNu2LQBoVjX4+++/i/R7cr/C2ICyfd86dOiAuLi4Ynlct24dJBIJ2rVrV+br0JaNjQ2aN2+OzZs3F4lVrVZj/fr1qFWrVpHJVs9r+vTp6NmzJ95//32dHbMsdDmkITg4uMh7+d8Fb3R0NDp27IhatWohKiqq1OKpb9++yMzMxK+//lqkfe3atXB3d0fz5s2f/6KfQiKRFLvr/Pfffxdbmzs0NBQZGRnYtWtXkfayPiBj4sSJsLa2xtixY4sV09rEtmPHDty6dUun+3777bdYv349li5diu3bt+Px48cYPnz4U4+vze/2f3NxccGwYcPw2muv4dKlSyWurEPGj3d4yeQ5ODhgypQp+PDDD7FhwwYMHjwY//3vf9G1a1d07twZw4YNg4eHBx4+fIj4+HicO3cOv/zyC4CCAmD79u1o3749Zs6cCWtrayxbtgxZWVlFzvH1119j//796N69O7y8vJCbm6tZ5qljx456uzZfX19YWVnhhx9+QP369WFrawt3d/enDkkAgCtXruDEiRNQq9V48OABTp48iVWrViE9PR3r1q1DgwYNSt33999/x/Lly9GnTx/UqVMHQghs3rwZjx8/RqdOnTT9goODcfDgQfz2229wc3ODnZ2d5u6KtpycnDB69GgkJyejXr162LlzJ1auXInRo0fDy8sLQMEQh44dO2LBggVwcHCAt7c3/vjjD83qB/8WHBwMAPj000/RtWtXSKVSNGzYEBYWFsX6Tpw4EevWrUP37t0xZ84ceHt7Y8eOHVi+fDlGjx6t04KzJAsWLECnTp3Qrl07fPDBB7CwsMDy5csRExODjRs36vQpaIMHD8bgwYN1dryysrCw0GqlgPK4dOmS5r04f/58XLlypcjH176+vqhZsyYAoGvXrujUqRNGjx6N9PR01K1bFxs3bsTu3buxfv36Z360/rx69OiBuXPnYtasWQgNDcWlS5cwZ84c+Pj4QKlUavoNHToUX3zxBQYPHox58+ahbt262LVrF/bs2QMAMDN7+j0tX19fbNy4Ea+99hqCg4MxevRoNGnSBHK5HKmpqdi7dy8AFFmyrEePHlizZg0CAgLQsGFDnD17Fp999plmWb1nXVdZ9r1w4QLGjRuHoUOHaorcVatWoX///liyZAkmTJhQ6jnK+ru9efPm6NGjBxo2bAgHBwfEx8fj+++/R4sWLWBtbf3MayEjZNApc0Q6VNqDJ4QQIicnR3h5eQk/Pz/NDPq//vpLDBgwQDg7OwuZTCZcXV1F+/btxddff11k3z///FO89NJLQi6XC1dXVzF58mTxzTffFFlN4Pjx46Jv377C29tbyOVy4eTkJEJDQ8X27ds1x3nagycAiFmzZmlel3WVBiGE2LhxowgICBAymazYcZ5UOKu58Mvc3Fw4OTmJFi1aiKlTp4qkpKRSv6+F13rx4kXx2muvCV9fX2FlZSWqVasmXnzxRbFmzZoi+0VHR4tWrVoJa2trAUCEhoYWOV5JeSptlYYGDRqIgwcPimbNmgm5XC7c3NzE1KlTi83EvnPnjujfv79wdHQU1apVE4MHD9asKvHv2eJ5eXli1KhRombNmkIikRQ555OrNAghxPXr18WgQYOEk5OTkMlkwt/fX3z22WdFZq9rk9+SPG3/I0eOiPbt2wsbGxthZWUlXnrpJfHbb7+V+L0r6fuq7fn+7XlWabh3755W++lT4fentK8nVxPIyMgQ48aNE66ursLCwkI0bNjwmSttFCrtvVraQ1S8vb1F9+7dNa/z8vLEBx98IDw8PISlpaVo0qSJ2Lp1qxg6dKjw9vYusm9ycrLo16+fsLW1FXZ2duLll18WO3fuLHFlk9JcvXpVvPfee8Lf319YWVkJuVwuvL29xSuvvCK2bNlS5OEljx49EiNHjhTOzs7C2tpatG7dWhw5ckSEhoZq3uNClLxKQ1n2zczMFAEBASIwMLDYw0rGjh0rZDKZZsWU0lZmKcvv9o8//lg0a9ZMODg4CLlcLurUqSMmTpwo7t+/X6bvGRkfiRA6XtSQiIiIDCYyMhLTp09HcnJyme68ElUFHNJARERkpJYuXQoACAgIgEKhwP79+/HVV19h8ODBLHaJ/oUFLxERlYlarS51XdtCTz6djPTL2toaX3zxBZKSkpCXlwcvLy989NFHmD59uqFDI6pUOKSBiIjKJCIiArNnz35qn8TERM2KGURElQULXiIiKpPbt29r1v8tTWmrXRARGRILXiIiIiIyaXzwBBERERGZNM4uKIFarcbt27dhZ2en00XdiYiIiEg3hBDIyMiAu7v7Mx+0woK3BLdv39bJc+SJiIiISL9u3LjxzGX4WPCWoPBZ2zdu3CjySEV9USgU2Lt3L8LDwyGTyfR+PtI95tD4MYfGjzk0bsyf8avoHKanp8PT01NTtz0NC94SFA5jsLe3r7CC19raGvb29nyTGynm0Pgxh8aPOTRuzJ/xM1QOyzL8lJPWiIiIiMikseAlIiIiIpPGgpeIiIiITBrH8JaTEAJKpRIqleq5j6VQKGBubo7c3FydHI+0I5VKYW5uziXoiIiITBQL3nLIz8/HnTt3kJ2drZPjCSHg6uqKGzdusOgyEGtra7i5ufGRqERERCaIBa+W1Go1EhMTIZVK4e7uDgsLi+cuUtVqNTIzM2Fra/vMhZNJt4QQyM/Px71795CYmAg/Pz/mgIiIyMSw4NVSfn4+1Go1PD09YW1trZNjqtVq5Ofnw9LSksWWAVhZWUEmk+H69euaPBAREZHpYHVVTixMTQvzSUREZLr4V56IiIiITBoLXiIiIiIyaSx4iYiIiMikseCtQoYNGwaJRIJ33nmn2LYxY8ZAIpFg2LBhFR+YllauXIk2bdrAwcEBDg4O6NixI06dOlWs3/Lly+Hj4wNLS0s0bdoUR44cMUC0REREZGgseKsYT09P/Pjjj8jJydG05ebmYuPGjfDy8ir3cfPz83URXpkcPHgQr732Gg4cOIDjx4/Dy8sL4eHhuHXrlqbPTz/9hAkTJmDatGk4f/482rRpg65duyI5ObnC4iQiIqpqcivp87O4LJkOCCGQoyh/htVqNXLyVTDPV2q9WoCVTKrVOsBNmjTBtWvXsHnzZrz++usAgM2bN8PT0xN16tQp83HCwsIQFBQECwsLrFu3Dg0aNMChQ4e0ir28fvjhhyKvV65ciU2bNuGPP/7AkCFDAACLFy/GyJEjMWrUKADAkiVLsGfPHqxYsQILFiyokDiJiIiqinsZeZj3eyyOX5Kie1c1ZDJDR1QUC14dyFGoEDhzj0HOHTenM6wttEvj8OHDsXr1ak3B+91332HEiBE4ePCgVsdZu3YtRo8ejT///BNCiBL7/PDDD3j77befepz//ve/mljKIzs7GwqFAo6OjgAK7jafPXsWH3/8cZF+4eHhOHbsWLnPQ0REREWp1QIbTiVj4e6LSM9VQgLgxLUH6NjA3dChFcGCtwp64403MGXKFCQlJUEikeDPP//Ejz/+qHXBW7duXSxcuPCpfXr16oXmzZs/tY+Li4tW533Sxx9/DA8PD3Ts2BEAcP/+fahUqmLHdXFxQUpKynOdi4iIiArE3ErDtK0x+OvGYwBAA3c7dHZ6hNB6NQ0bWAlY8OqAlUyKuDmdy72/Wq1GRnoG7OztyjWkQVs1atRA9+7dsXbtWggh0L17d9SoUUPr4zRr1uyZfezs7GBnZ6f1sQEgOTkZgYGBmtdTp07F1KlTi/RZuHAhNm7ciIMHDxZ7QtqTQz2EEM/9GGgiIqKqLjNPicV7L2PNsUSoBWArN8cH4fXwajMP7Nm9y9DhlYgFrw5IJBKthxX8m1qthtJCCmsL8wp74teIESPw7rvvAgCWLVtWrmPY2Ng8s8/zDGlwd3dHdHS05nXhkIVCixYtQmRkJPbt24eGDRtq2mvUqAGpVFrsbm5qaupz300mIiKqqoQQ2BWTgjm/xSElPRcA0KOhG2b0CISLvSUUCoWBIywdC94qqkuXLpqVFTp3Lv/d6Wd5niEN5ubmqFu3bonbPvvsM8ybNw979uwpdqfZwsICTZs2RVRUFPr27atpj4qKQu/evbW8AiIiIkp+kI2Z22Nw8NI9AIC3kzXm9g5C20o4fKEkLHirKKlUivj4eM2/9eV5hjSUZuHChZgxYwY2bNiA2rVra+7k2trawtbWFgAwadIkvPHGG2jWrBlatGiBb775BsnJySWuQUxEREQly1OqsPLwNfxnfwLylGpYSM3wTpgvxoT5wrIcwyoNhQVvFWZvb2/oEMpl+fLlyM/PR//+/Yu0z5o1CxEREQCAgQMH4sGDB5gzZw7u3LmDoKAg7Ny5E97e3gaImIiIyPgcv/oA07dewNV7WQCAVnWdMLd3EOrUtDVwZNpjwVuFrFmz5qnbt27dWuZjabuigy4lJSWVqd+YMWMwZswY/QZDRERkYu5n5iFyRzw2ny94oFMNWwvM6BGIXo3cjXbyNwteIiIiIoJaLbDxdDI+3fXPmroSYHBzb3zQ2R/VrCrZkyS0xIKXinlyObAnxcXFPddjiImIiKhyib2dhulbY3A++TEAoIG7Peb3DUZjz+oGjUtXWPBSMU8uB1bSdiIiIjJ+mXlKfBF1Gav//P+auu+H18MbL3nDXFoxS6VWBBa8VMzTlgMjIiIi4yeEwO6YFMz+15q63Ru6YeY/a+qaGha85SSEMHQIpEPMJxERVRU3HmZj5rYYHPhnTV0vR2vM6d0AYf7OBo5Mf1jwakkmKxi0nZ2dDSsrKwNHQ7qSnZ0N4P/5JSIiMjX5SjVWHrmGr/64gjylGjKpBKNDfTGmXV2jWlO3PFjwakkqlaJ69epITU0FAFhbWz/3Eh1qtRr5+fnIzc2tsEcLUwEhBLKzs5Gamorq1avr9SEcREREhnLi2gNM3xqDhNRMAECLOk6Y2ycIdZ2Nb03d8mDBWw6urq4AoCl6n5cQAjk5ObCysjLa9e2MXfXq1TV5JSIiMhX3M/MQuTMem8/9f03d6d0D0bux8a6pWx4seMtBIpHAzc0Nzs7OUCgUz308hUKBw4cPo23btvxI3QBkMhnv7BIRkUlRqwV+PH0Dn+6+iLQcBSQSYNCLXviwcwCqWVe9WoMF73OQSqU6KZSkUimUSiUsLS1Z8BIREdFzibudjulbL+DcP2vqBrrZY37fIIR4ORg2MANiwUtERERkAjLzlFgSdRmrjyVBpRawsZDi/XB/DGlhWmvqlgcLXiIiIiIjJoTAnti7mP1bLO6kFayp2y3YFTN7NIBrNdNbU7c8WPASERERGakbD7MRsT0Wf1wsmEjv6WiFOb2D0M6E19QtDxa8REREREZGoVLj2yOJ+PKPy8hVFKyp+1bbOnivvZ/Jr6lbHix4iYiIiIzIqcSHmL71Ai7fLVhT96U6jpjXJwh1ne0MHFnlxYKXiIiIyAg8zMrHJ7vi8fOZmwAARxsLTOtWH/2aeFSpNXXLw6BT9g4fPoyePXvC3b1g8eOtW7cW2S6EQEREBNzd3WFlZYWwsDDExsY+87iPHz/G2LFj4ebmBktLS9SvXx87d+7U01UQERER6Y8QAj+fuYEOnx/UFLuvveiJ/e+H4uWmtVjsloFB7/BmZWWhUaNGGD58OF5++eVi2xcuXIjFixdjzZo1qFevHubNm4dOnTrh0qVLsLMr+bZ9fn4+OnXqBGdnZ2zatAm1atXCjRs3Su1PREREVFldvpuB6VticCrpIQAgwNUO8/sGoam3o4EjMy4GLXi7du2Krl27lrhNCIElS5Zg2rRp6NevHwBg7dq1cHFxwYYNG/D222+XuN93332Hhw8f4tixY5qHOHh7e+vnAoiIiIj0ICdfha/2X8HKw9egVAtYyaSY2MkPw1v5QFbF19Qtj0o7hjcxMREpKSkIDw/XtMnlcoSGhuLYsWOlFrzbt29HixYtMHbsWGzbtg01a9bEoEGD8NFHH5X6VLS8vDzk5eVpXqenpwMoeOSvLh4d/CyF56iIc5F+MIfGjzk0fsyhcWP+/u/ApXuY83s8bj4uWFO3Y0BNzOgeAPfqVoBaBYVaZeAIS1bROdTmPJW24E1JSQEAuLi4FGl3cXHB9evXS93v2rVr2L9/P15//XXs3LkTV65cwdixY6FUKjFz5swS91mwYAFmz55drH3v3r2wtrZ+jqvQTlRUVIWdi/SDOTR+zKHxYw6NW1XO3+M84NckM/z9sOAObnULgf4+agQ73EH0sTuINmx4ZVZROczOzi5z30pb8BZ6ciC2EOKpg7PVajWcnZ3xzTffQCqVomnTprh9+zY+++yzUgveKVOmYNKkSZrX6enp8PT0RHh4OOzt7XVzIU+hUCgQFRWFTp06aYZhkHFhDo0fc2j8mEPjVpXzp1Sp8f3JG/jyjwRk5asgNZNgeEtvvBtWBzbySl+qaVR0Dgs/kS+LSvtddHV1BVBwp9fNzU3TnpqaWuyu77+5ublBJpMVGb5Qv359pKSkID8/HxYWFsX2kcvlkMvlxdplMlmFvukq+nyke8yh8WMOjR9zaNyqWv6ibzzG1M0XEHenoHhr4lUd8/sGo76b/m+46UtF5VCbc1TaUc8+Pj5wdXUtcls8Pz8fhw4dQsuWLUvdr1WrVkhISIBarda0Xb58GW5ubiUWu0REREQVLS1HgelbL6Dv8j8Rdycd1axkWNAvGJveaWnUxW5lZdA7vJmZmUhISNC8TkxMRHR0NBwdHeHl5YUJEyYgMjISfn5+8PPzQ2RkJKytrTFo0CDNPkOGDIGHhwcWLFgAABg9ejT+85//YPz48Xjvvfdw5coVREZGYty4cRV+fURERET/JoTA9r9uY+7v8bifWTBhvl8TD0ztVh81bIt/2ky6YdCC98yZM2jXrp3mdeE42qFDh2LNmjX48MMPkZOTgzFjxuDRo0do3rw59u7dW2RN3eTkZJiZ/f9GtaenJ/bu3YuJEyeiYcOG8PDwwPjx4/HRRx9V3IURERERPeHavUzM2BaDPxMeAADq1LTBvD5BaOlbw8CRmT6DFrxhYWEQQpS6XSKRICIiAhEREaX2OXjwYLG2Fi1a4MSJEzqIkIiIiOj55CpUWHHwKlYcvIp8lRpyczO8174u3mxbB3LzkpdMJd2qtJPWiIiIiIzd0Sv3MWNbDBLvZwEA2taribm9G8DbycbAkVUtLHiJiIiIdCw1Ixfzd8RjW/RtAICznRyzejZAt2DXpy6vSvrBgpeIiIhIR1RqgQ2nkrFw90Vk5CphJgGGtKiN98Prwc6y6iy3Vtmw4CUiIiLSgZhbaZi2NQZ/3XgMAAj2qIb5fYPQsFZ1g8ZFLHiJiIiInktmnhKL917GmmOJUAvAVm6OyZ39Mfglb0jNOHyhMmDBS0RERFQOQgjsjknB7N/ikJKeCwDo0dANM3oEwsXe0sDR0b+x4CUiIiLS0o2H2Zi5LQYHLt0DAHg7WWNO7yCE1qtp4MioJCx4iYiIiMooX6nGyiPX8J/9V5CrUEMmlWB0qC/GtKsLSxnX1K2sWPASERERlcHJaw8wfWsMrqRmAgBa1HHC3D5BqOtsa+DI6FlY8BIRERE9xcOsfCzYGY9fzt4EADjZWGB6j/ro09iDa+oaCRa8RERERCVQqwU2nb2JyF3xeJytAAC89qIXPurij+rWFgaOjrTBgpeIiIjoCZdSMjB96wWcTnoEAAhwtcP8vsFo6u1g4MioPFjwEhEREf0jO1+Jr/5IwLdHrkGpFrC2kGJix3oY1qo2ZFIzQ4dH5cSCl4iIiAjAH/F3MXNbLG49zgEAhAe6YFavBvCobmXgyOh5seAlIiKiKu1OWg4itsdiT+xdAIBHdStE9GqAToEuBo6MdIUFLxEREVVJSpUaa44l4Yuoy8jKV8HcTIKRbXwwvoMfrC1YIpkSZpOIiIiqnHPJjzBtSwzi76QDAJp6O2B+3yAEuNobODLSBxa8REREVGWkZSvw6Z6L2HgqGUIA1a1lmNI1AK809YSZGdfUNVUseImIiMjkCSGwLfo25u2Iw/3MfABA/6a1MKVrAJxs5QaOjvSNBS8RERGZtKv3MjFjawyOXX0AAKjrbIt5fYLwUh0nA0dGFYUFLxEREZmkXIUKyw9exdcHryJfpYbc3AzjOvjhzTZ1YGHONXWrEha8REREZHKOXLmHGVtjkPQgGwAQWq8m5vYOgpeTtYEjI0NgwUtEREQmIzUjF/N+j8f2v24DAFzs5ZjVswG6BrlCIuGktKqKBS8REREZPZVa4IeT1/HZ7kvIyFPCTAIMbVkbkzrVg52lzNDhkYGx4CUiIiKjFnMrDdO2XMBfN9MAAA1rVUNk32AEeVQzcGRUWbDgJSIiIqOUkavA4qjLWHssCWoB2MnN8WEXfwxq7g0p19Slf2HBS0REREZFCIFdMSmY/Vss7qbnAQB6NnLHjO714WxvaeDoqDJiwUtERERGI/lBNmZuj8HBS/cAAN5O1pjbOwht69U0cGRUmbHgJSIiokovX6nGyiPX8NUfV5CnVMNCaoZ3wnwxJswXljKpocOjSo4FLxEREVVqJ649wPStMUhIzQQAtPR1wtw+QfCtaWvgyMhYsOAlIiKiSulBZh4id17Er+duAgBq2FpgevdA9G7szjV1SSsseImIiKhSUasFfj5zAwt2XURajgISCTDoRS982DkA1ay5pi5pjwUvERERVRqXUjIQ8ftFnLn+CABQ380e8/sGoYmXg4EjI2PGgpeIiIgMLjtfie3XzXDo5Ako1QLWFlJM6lQPw1rWhrnUzNDhkZFjwUtEREQGtS/uLmZui8HtNDMAAl0auGJmz0C4V7cydGhkIljwEhERkUHcepyD2dtjsTfuLgDAUS6w4JUm6BzkbuDIyNSw4CUiIqIKpVCpsebPJHyx7zKy81UwN5NgRCtv1M1LQHt/PkCCdI8FLxEREVWYs9cfYdqWC7iYkgEAaObtgPl9g1HHyRI7dyYYODoyVQYdBX748GH07NkT7u4F6+lt3bq1yHYhBCIiIuDu7g4rKyuEhYUhNja2zMf/8ccfIZFI0KdPH90GTkRERFpJy1ZgyuYLeHnFMVxMyUB1axkWvtwQP7/dAv6udoYOj0ycQQverKwsNGrUCEuXLi1x+8KFC7F48WIsXboUp0+fhqurKzp16oSMjIxnHvv69ev44IMP0KZNG12HTURERGUkhMCW8zfR/vOD2HgqGQDwStNa2P9+GAa84AkzMz5AgvTPoEMaunbtiq5du5a4TQiBJUuWYNq0aejXrx8AYO3atXBxccGGDRvw9ttvl3pclUqF119/HbNnz8aRI0fw+PFjfYRPRERET3H1Xiamb4nB8WsPAAB1nW0xv08QmtdxMnBkVNVU2jG8iYmJSElJQXh4uKZNLpcjNDQUx44de2rBO2fOHNSsWRMjR47EkSNHnnmuvLw85OXlaV6np6cDABQKBRQKxXNcRdkUnqMizkX6wRwaP+bQ+DGHlUeuQoWvDyfimyOJUKgELGVmeDfMF8NbesPC3KzEHDF/xq+ic6jNeSptwZuSkgIAcHFxKdLu4uKC69evl7rfn3/+iVWrViE6OrrM51qwYAFmz55drH3v3r2wtrYu83GeV1RUVIWdi/SDOTR+zKHxYw4NK/6xBJuumeF+XsFQhfrV1XjFRwmnzHjs2xv/zP2ZP+NXUTnMzs4uc99KW/AWkkiKju0RQhRrK5SRkYHBgwdj5cqVqFGjRpnPMWXKFEyaNEnzOj09HZ6enggPD4e9vX35AteCQqFAVFQUOnXqBJmMzwg3Rsyh8WMOjR9zaFipGXmI3HkJO+L/uWFlJ8f07gHoHOhc6t/tf2P+jF9F57DwE/myqLQFr6urK4CCO71ubm6a9tTU1GJ3fQtdvXoVSUlJ6Nmzp6ZNrVYDAMzNzXHp0iX4+voW208ul0Mulxdrl8lkFfqmq+jzke4xh8aPOTR+zGHFUqkF1p+4jkV7LiEjTwkzCTCspQ8mhdeDrVz7MoP5M34VlUNtzlFpC14fHx+4uroiKioKISEhAID8/HwcOnQIn376aYn7BAQE4MKFC0Xapk+fjoyMDHz55Zfw9PTUe9xERERVxYWbaZi65QIu3EoDADSqVQ3z+wYjyKOagSMjKsqgBW9mZiYSEv6/yHRiYiKio6Ph6OgILy8vTJgwAZGRkfDz84Ofnx8iIyNhbW2NQYMGafYZMmQIPDw8sGDBAlhaWiIoKKjIOapXrw4AxdqJiIiofNJzFVi89zLWHU+CWgB2lub4sEsABr3oBSmXGaNKyKAF75kzZ9CuXTvN68JxtEOHDsWaNWvw4YcfIicnB2PGjMGjR4/QvHlz7N27F3Z2/1+gOjk5GWZmBl1OmIiIqEoQQmDnhRTM/i0WqRkFqxv1buyOad3rw9nO0sDREZXOoAVvWFgYhBClbpdIJIiIiEBERESpfQ4ePPjUc6xZs6Z8wREREZHG9QdZmLktFocu3wMA+NSwwdzeQWjtV/ZJ4kSGUmnH8BIREZHh5SlVWHn4Gv6zPwF5SjUspGYYHeaL0WG+sJRJDR0eUZmw4CUiIqISHb/6ANO3XsDVe1kAgFZ1nTC3dxDq1LQ1cGRE2mHBS0REREU8yMzD/J3x2HzuFgCghq0FZvQIRK9G7mVaU5eosmHBS0RERAAAtVrgl7M3sGDXRTzOVkAiAQa96IUPuwSgmhXXxiXjxYKXiIiIcPluBqZtuYDTSY8AAPXd7BHZNwghXg4Gjozo+bHgJSIiqsJy8lX4av8VrDx8DUq1gLWFFJM61cOwlrVhLuWyn2QaWPASERFVUQcupmLGthjcfJQDAAgPdMGsXg3gUd3KwJER6RYLXiIioiomJS0Xc36Pxc4LKQAA92qWmN07CJ0CXQwcGZF+sOAlIiKqIlRqgXXHk/D53svIzFNCaibByNY+GN/BDzZylgRkuvjTTUREVAX8ffMxpm2JwYVbaQCAEK/qmN8nGIHu9gaOjEj/WPASERGZsIxcBT7fexnrjidBLQB7S3N81DUAr73gBTMzrqlLVQMLXiIiIhMkhMDOCymY/VssUjPyAAC9G7tjevdA1LSTGzg6oor1XAVvXl4e5HK+aYiIiCqTGw+zMWNbDA5eugcAqO1kjXl9gtHar4aBIyMyDK0K3j179mDjxo04cuQIkpOToVarYW1tjSZNmiA8PBzDhw+Hu7u7vmIlIiKip8hXqvHt0Wv46o8ryFWoYSE1wzthvhgT5gtLmdTQ4REZTJkK3q1bt+Kjjz5CWloaunXrhsmTJ8PDwwNWVlZ4+PAhYmJisG/fPsydOxfDhg3D3LlzUbNmTX3HTkRERP84nfQQ07ZcwOW7mQCAFnWcMK9vEHxr2ho4MiLDK1PBGxkZiUWLFqF79+4wMyv+1JUBAwYAAG7duoUvv/wS69atw/vvv6/bSImIiKiYR1n5+GTXRfx05gYAwMnGAtO610ffEA9IJJyURgSUseA9depUmQ7m4eGBhQsXPldARERE9GxCCPx67hYid8bjYVY+AOC1Fz3xUZcAVLe2MHB0RJWL1pPWDh8+jICAADg7OxdpVygUOH78ONq2bauz4IiIiKi4hNRMTNtyAScTHwIA/F3sML9vEJrVdjRwZESVk9YFb1hYGFxcXLB582a0aNFC0/7w4UO0a9cOKpVKpwESERFRgVyFCssOJODrQ1ehUAlYyswwoWM9jGztA5m0+JBDIipQrnfHq6++ig4dOmDNmjVF2oUQuoiJiIiInnD48j10XnIY/9mfAIVKoH2AM6ImhuKdUF8Wu0TPoPUdXolEgilTpqBNmzYYOnQo/v77b3z++eeabURERKQ7qRm5mPt7PH776zYAwNXeEhG9AtG5gSv/7hKVkdYFb+Fd3H79+sHHxwe9e/dGXFwcvvzyS50HR0REVFWp1AIbTl7Hwj2XkJGrhJkEGNqyNt4P94etnA9KJdLGc71jQkJCcOrUKfTp0wcdOnTQVUxERERVWuztNEzdEoO/bjwGADSsVQ2RfYMR5FHNsIERGSmtC96hQ4fCyspK89rV1RWHDh3CW2+9hcOHD+s0OCIioqokK0+JxVGXsfrPRKgFYCs3x+TO/hj8kjekZhy+QFReWhe8q1evLtYml8uxdu1anQRERERUFe2JTUHE9ljcScsFAHRv6IaZPQLhYm9p4MiIjF+ZCt7k5GR4eXmV+aC3bt2Ch4dHuYMiIiKqKm4+ykbE9jjsi78LAPB0tMLc3kEI83d+xp5EVFZlWsfkhRdewJtvvvnUJ66lpaVh5cqVCAoKwubNm3UWIBERkSlSqNT45vBVdFp8GPvi70ImlWBsO1/snRDKYpdIx8p0hzc+Ph6RkZHo0qULZDIZmjVrBnd3d1haWuLRo0eIi4tDbGwsmjVrhs8++wxdu3bVd9xERERG6+z1R5i25QIupmQAAF6s7Yj5fYPg52Jn4MiITFOZCl5HR0csWrQI8+bNw86dO3HkyBEkJSUhJycHNWrUwOuvv47OnTsjKChI3/ESEREZrbRsBT7dcxEbTyVDCMDBWoYp3erjlaa1uKYukR5pNWnN0tIS/fr1Q79+/fQVDxERkckRQmBb9G3M2xGH+5n5AID+TWtharf6cLSxMHB0RKZPq4L3+vXr2Lt3LxQKBcLCwhAYGKivuIiIiExC4v0szNgag6MJ9wEAdZ1tMa9PEF6q42TgyIiqjjIXvIcPH0a3bt2QnZ1dsKO5OdauXYvXXntNb8EREREZqzylCisOXsXyg1eRr1RDbm6G99rXxVttfWFhXqY540SkI2V+x82YMQPt2rXDzZs38eDBA4wYMQIffvihPmMjIiIySscS7qPrkiNYsu8K8pVqtK1XE3sntsW77f1Y7BIZQJnv8F64cAGHDx+Gu7s7AODzzz/HypUr8ejRIzg4OOgtQCIiImNxPzMP83fEY8v5WwCAmnZyzOwRiB4N3TgpjciAylzwPn78GM7O/18X0MbGBtbW1nj8+DELXiIiqtLUaoGfztzAJ7suIi1HAYkEeOMlb3zQ2R/2ljJDh0dU5Wk1aS0uLg4pKSma10IIxMfHIyMjQ9PWsGFD3UVHRERUyV1MSce0LTE4e/0RACDQzR6R/YLR2LO6YQMjIg2tCt4OHTpACFGkrUePHpBIJBBCQCKRQKVS6TRAIiKiyig7X4kv/7iCVUcSoVQL2FhIMSncH0NbeMNcynG6RJVJmQvexMREfcZBRERkNP6Iv4uZ22Jx63EOAKBLA1fM6hUIt2pWBo6MiEpS5v+Cent7l+lLG4cPH0bPnj3h7u4OiUSCrVu3FtkuhEBERATc3d1hZWWFsLAwxMbGPvWYK1euRJs2beDg4AAHBwd07NgRp06d0iouIiKiktxJy8E735/FyLVncOtxDjyqW2HV0Gb4+o2mLHaJKjGDfuaSlZWFRo0aYenSpSVuX7hwIRYvXoylS5fi9OnTcHV1RadOnYqMGX7SwYMH8dprr+HAgQM4fvw4vLy8EB4ejlu3bunrMoiIyMQpVWqsOpqIjp8fwu7YFEjNJHg7tA6iJrVFh/ouhg6PiJ5BqzG8uta1a1d07dq1xG1CCCxZsgTTpk3TPMp47dq1cHFxwYYNG/D222+XuN8PP/xQ5PXKlSuxadMm/PHHHxgyZIhuL4CIiEzeXzceY+qWC4i9nQ4AaOJVHZH9ghHgam/gyIiorAxa8D5NYmIiUlJSEB4ermmTy+UIDQ3FsWPHSi14n5SdnQ2FQgFHR8dS++Tl5SEvL0/zOj294JeaQqGAQqEo5xWUXeE5KuJcpB/MofFjDo2frnOYkavA4n0J+OHUDQgB2Fua48PO9fBKEw+YmUn4s6JjfA8av4rOoTbnqbQFb+HyZy4uRT8qcnFxwfXr18t8nI8//hgeHh7o2LFjqX0WLFiA2bNnF2vfu3cvrK2ty3yu5xUVFVVh5yL9YA6NH3No/J43h0IA0Q8k2JxkhnRFwcMimtVQo0/tXNil/o3du//WRZhUCr4HjV9F5TA7O7vMfbUueCMiIjB8+HCtJ6iV15NPpilc/qwsFi5ciI0bN+LgwYOwtLQstd+UKVMwadIkzev09HR4enoiPDwc9vb6/8hKoVAgKioKnTp1gkzGBcqNEXNo/JhD46eLHCY/zMbs3+Nx+MoDAEBtJ2vM7lkfLX2ddBkqlYDvQeNX0Tks/ES+LLQueH/77TfMmzcPoaGhGDlyJPr16/fUYrK8XF1dARTc6XVzc9O0p6amFrvrW5JFixYhMjIS+/bte+bDMORyOeRyebF2mUxWoW+6ij4f6R5zaPyYQ+NXnhzmK9VYeeQavvrjCvKUalhIzTA6zBejw3xhKZPqKVIqCd+Dxq+icqjNObRepeHs2bM4d+4cGjZsiIkTJ8LNzQ2jR4/G6dOntT3UU/n4+MDV1bXIbfH8/HwcOnQILVu2fOq+n332GebOnYvdu3ejWbNmOo2LiIhMy8lrD9DtqyP4bM8l5CnVaOnrhN0T2mBip3osdolMRLmWJWvYsCG++OIL3Lp1C9999x1u3bqFVq1aITg4GF9++SXS0tLKdJzMzExER0cjOjoaQMFEtejoaCQnJ0MikWDChAmIjIzEli1bEBMTg2HDhsHa2hqDBg3SHGPIkCGYMmWK5vXChQsxffp0fPfdd6hduzZSUlKQkpKCzMzM8lwqERGZqIdZ+Zj8y18Y+M0JJKRmwsnGAl8MbIQfRjVHnZq2hg6PiHToudbhVavVyM/PR15eHoQQcHR0xIoVK+Dp6YmffvrpmfufOXMGISEhCAkJAQBMmjQJISEhmDlzJgDgww8/xIQJEzBmzBg0a9YMt27dwt69e2FnZ6c5RnJyMu7cuaN5vXz5cuTn56N///5wc3PTfC1atOh5LpWIiEyEEAI/n7mBDp8fxC9nbwIAXnvRC/vfD0PfkFplnidCRMajXKs0nD17FqtXr8bGjRshl8sxZMgQLFu2DHXr1gUAfP755xg3bhwGDhz41OOEhYVBCFHqdolEgoiICERERJTa5+DBg0VeJyUllfUyiIioiklIzcDULTE4lfgQABDgaof5fYPQ1Lv0pSuJyPhpXfA2bNgQ8fHxCA8Px6pVq9CzZ09IpUXHOA0ZMgSTJ0/WWZBERETPI1ehwtL9Cfjv4atQqASsZFJM6OiHEa19IJMa9KGjRFQBtC54X3nlFYwYMQIeHh6l9qlZsybUavVzBUZERKQLhy7fw4ytMUh+WLBmZ4cAZ8zu3QC1HCpunXUiMiyt/1srhICDg0Ox9pycHMyZM0cnQRERET2v1PRcvLvhHIZ+dwrJD7Pham+Jrwc3xbdDm7HYJapitC54Z8+eXeKKB9nZ2SU+rYyIiKgiqdQC644nocPnh/D733dgJgFGtvbBvvdD0SXIlZPSiKogrYc0lPaks7/++guOjhz0T0REhnMzCxjwzUn8favgCUyNalXD/L7BCPKoZuDIiMiQylzwOjg4QCKRQCKRoF69ekWKXpVKhczMTLzzzjt6CZKIiOhpMvOUWLT7Itb+LYVAOuzk5viwiz8GNfeG1Ix3dImqujIXvEuWLIEQAiNGjMDs2bNRrdr//7dsYWGB2rVro0WLFnoJkoiIqCRCCOyJvYvZv8XiTlouAAm6B7liVq8GcLbX/WPvicg4lbngHTp0KICCR/62bNmSz7kmIiKDuvkoG7O2xeKPi6kAAE8HK3R3zcT7AxvybxQRFVGmgjc9PR329vYAgJCQEOTk5CAnJ6fEvoX9iIiI9EGhUuO7o4lYsu8KchQqyKQSvN3WF2+38cb+qD2GDo+IKqEyFbwODg64c+cOnJ2dUb169RInrRVOZlOpVDoPkoiICADOXn+IaVticDElAwDwoo8jIvsGoa6zHRQKhYGjI6LKqkwF7/79+zUrMBw4cECvARERET3pcXY+Pt19CRtPJQMAHKxlmNqtPvo3rcVlxojomcpU8IaGhpb4byIiIn0SQmBr9C3M+z0eD7LyAQADmtXCx13rw9HGwsDREZGx0PrBE7t378bRo0c1r5ctW4bGjRtj0KBBePTokU6DIyKiquvavUy8/u1JTPzpLzzIykddZ1v89NZLWNi/EYtdItKK1gXv5MmTkZ5esKD3hQsXMGnSJHTr1g3Xrl3DpEmTdB4gERFVLbkKFb6IuowuS47g2NUHkJubYXJnf+wc1wbN6zgZOjwiMkJaP2ktMTERgYGBAIBff/0VPXv2RGRkJM6dO4du3brpPEAiIqo6/ky4j+lbY5B4PwsAEFqvJub2DoKXk7WBIyMiY6Z1wWthYYHs7GwAwL59+zBkyBAAgKOjo+bOLxERkTbuZeRh/o44bI2+DQBwtpNjVs8G6BbsyklpRPTctC54W7dujUmTJqFVq1Y4deoUfvrpJwDA5cuXUatWLZ0HSEREpkutFvjx9A18sise6blKSCTAkJe88X5nf9hb8uERRKQbWhe8S5cuxZgxY7Bp0yasWLECHh4eAIBdu3ahS5cuOg+QiIhMU/yddEzbcgHnkh8DAII87BHZNxgNa1U3aFxEZHq0Lni9vLzw+++/F2v/4osvdBIQERGZtux8Jb7cdwXfHk2ESi1gYyHF++H+GNLCG+ZSredSExE9k9YFLwCo1WokJCQgNTUVarW6yLa2bdvqJDAiIjI9++LuYtb2WNx6XPB4+q5BrpjVswFcq1kaODIiMmVaF7wnTpzAoEGDcP36dQghimzjo4WJiKgktx/nYPZvsdgTexcA4FHdCnP7NED7ABcDR0ZEVYHWBe8777yDZs2aYceOHXBzc+PsWSIiKpVSpcaaY0n4IuoysvJVMDeTYFSbOhjXoS6sLcr1ISMRkda0/m1z5coVbNq0CXXr1tVHPEREZCKibzzG1M0XEHenYMnKpt4OmN83CAGu9gaOjIiqGq0L3ubNmyMhIYEFLxERlSg9V4HPdl/C+pPXIQRQzUqGKV0DMKCZJ8zM+KkgEVU8rQve9957D++//z5SUlIQHBwMmazoOokNGzbUWXBERGQ8hBD47e87mPt7HO5l5AEA+oV4YGr3+qhhKzdwdERUlWld8L788ssAgBEjRmjaJBIJhBCctEZEVEVdf5CF6VtjcOTKfQBAnRo2mNc3CC19axg4MiKichS8iYmJ+oiDiIiMUL5SjW8OX8V/9icgT6mGhbkZxobVxTthdSA3lxo6PCIiAOUoeL29vfURBxERGZmT1x5g2tYYJKRmAgBa162BuX2C4FPDxsCREREVVa5H2nz//fdo1aoV3N3dcf36dQDAkiVLsG3bNp0GR0RElc+jrHx8uOkvDPzmBBJSM1HD1gJfvtoY3498kcUuEVVKWhe8K1aswKRJk9CtWzc8fvxYM2a3evXqWLJkia7jIyKiSkIIgU1nb6LD4kP4+cxNAMBrL3rhj0lh6N3Yg+uyE1GlpXXB+5///AcrV67EtGnTIJX+f3xWs2bNcOHCBZ0GR0RElcPVe5l4beUJfPDLX3iYlQ9/Fzv8OroFFvQLRjVr2bMPQERkQOWatBYSElKsXS6XIysrSydBERFR5ZCrUGH5wav4+uBV5KvUsJSZYXyHehjVxgcyablGxRERVTitC14fHx9ER0cXm7y2a9cuBAYG6iwwIiIyrD8T7mP61hgk3i+4mdHOvybm9A6Cp6O1gSMjItKO1gXv5MmTMXbsWOTm5kIIgVOnTmHjxo1YsGABvv32W33ESEREFeh+Zh7m74jHlvO3AADOdnJE9GqArkGuHKdLREZJ64J3+PDhUCqV+PDDD5GdnY1BgwbBw8MDX375JV599VV9xEhERBVArRb4+cwNLNh1EWk5CkgkwJCXvPF+Z3/YW3KcLhEZL60LXgB488038eabb+L+/ftQq9VwdnbWdVxERFSBLt/NwNTNF3Dm+iMAQAN3e0T2DUYjz+qGDYyISAfKVfDev38fSUlJkEgkqF27to5DIiKiipKTr8JX+69g5eFrUKoFrC2kmNSpHoa1rA1zTkojIhOhVcEbGxuL0aNH488//yzSHhoaihUrVsDf31+nwRERkf4cvJSKGdticONhDgAgPNAFEb0awL26lYEjIyLSrTL/9z0lJQWhoaG4d+8eFi9ejJ07d2LHjh347LPPcOfOHbRp0wapqalanfzw4cPo2bMn3N3dIZFIsHXr1iLbhRCIiIiAu7s7rKysEBYWhtjY2Gce99dff0VgYCDkcjkCAwOxZcsWreIiIjJlqem5GLvhHIatPo0bD3PgXs0S37zRFN8MacZil4hMUpkL3i+++ALe3t44f/48xo8fj86dO6NLly6YNGkSzp07B09PT3zxxRdanTwrKwuNGjXC0qVLS9y+cOFCLF68GEuXLsXp06fh6uqKTp06ISMjo9RjHj9+HAMHDsQbb7yBv/76C2+88QYGDBiAkydPahUbEZGpUakF1h1PQofPD2HH33cgNZNgVGsfRE0KRXgDV0OHR0SkN2UueKOiovDRRx/B0tKy2DYrKytMnjwZe/bs0erkXbt2xbx589CvX79i24QQWLJkCaZNm4Z+/fohKCgIa9euRXZ2NjZs2FDqMZcsWYJOnTphypQpCAgIwJQpU9ChQwc+9piIqrTY22not+IYZm6LRUaeEo08q2P7u60wvUcgbOTlms5BRGQ0yvxb7tq1a2jSpEmp25s1a4Zr167pJCig4IluKSkpCA8P17TJ5XKEhobi2LFjePvtt0vc7/jx45g4cWKRts6dOz+14M3Ly0NeXp7mdXp6OgBAoVBAoVA8x1WUTeE5KuJcpB/MofEz1Rxm5Snx1f6rWHsiGSq1gK3cHO93qovXXvCE1ExiUtdrqjmsKpg/41fROdTmPGUueDMyMmBvb1/qdjs7O2RmZpb5xM+SkpICAHBxcSnS7uLiguvXrz91v5L2KTxeSRYsWIDZs2cXa9+7dy+srSvuiUJRUVEVdi7SD+bQ+JlSDi88lGBTohke5xc8LCLESY2+tXNR7UEM9uyOMXB0+mNKOayKmD/jV1E5zM7OLnNfrT7HysjIKHFIA1BwV1QIoc3hyuTJp/oIIZ75pB9t95kyZQomTZqkeZ2eng5PT0+Eh4c/tcjXFYVCgaioKHTq1AkyGRd3N0bMofEzpRzeScvF3B0XEXWpYCJxreqWiOhZH6H1aho4Mv0ypRxWRcyf8avoHBZ+Il8WZS54hRCoV6/eU7fr8pGTrq4FEyhSUlLg5uamaU9NTS12B/fJ/Z68m/usfeRyOeRyebF2mUxWoW+6ij4f6R5zaPyMOYdKlRprjiXhi6jLyMpXwdxMgjfb1sG49n6wspAaOrwKY8w5JObPFFRUDrU5R5kL3gMHDpQrmPLy8fGBq6sroqKiEBISAgDIz8/HoUOH8Omnn5a6X4sWLRAVFVVkHO/evXvRsmVLvcdMRGQof914jKlbLiD2dsEdj6beDojsGwx/VzsDR0ZEZHhlLnhDQ0N1fvLMzEwkJCRoXicmJiI6OhqOjo7w8vLChAkTEBkZCT8/P/j5+SEyMhLW1tYYNGiQZp8hQ4bAw8MDCxYsAACMHz8ebdu2xaefforevXtj27Zt2LdvH44eParz+ImIDC09V4HP91zCuhPXIQRgb2mOKd3qY2AzT5iZ6e5TNyIiY2bQtWjOnDmDdu3aaV4XjqMdOnQo1qxZgw8//BA5OTkYM2YMHj16hObNm2Pv3r2ws/v/HYvk5GSYmf1/dbWWLVvixx9/xPTp0zFjxgz4+vrip59+QvPmzSvuwoiI9EwIgV0xKYjYHovUjIJVZvqGeGBa9/qoYVt8iBYRUVVm0II3LCzsqRPdJBIJIiIiEBERUWqfgwcPFmvr378/+vfvr4MIiYgqnxsPszFzWwwOXLoHAKjtZI15fYLR2q+GgSMjIqqcuNo4EZGRUKjUWHU0EUv2XUauQg0LqRneCfPFmDBfWMqqzqQ0IiJtseAlIjICZ68/xLQtMbiYUvBo9eY+jpjfNxh1nW0NHBkRUeWndcG7Zs0aDBgwoEIfyEBEVFWlZSvwye6L2HgqGQDgYC3DtO6BeLmJh06XgiQiMmVmz+5S1JQpU+Dq6oqRI0fi2LFj+oiJiKjKE0JgW/QtdFh8UFPsDmhWC/vfD0P/prVY7BIRaUHrgvfmzZtYv349Hj16hHbt2iEgIACffvrpUx/dS0REZZd0PwtvrDqF8T9G435mPuo62+Knt17Cwv6N4GBjYejwiIiMjtYFr1QqRa9evbB582bcuHEDb731Fn744Qd4eXmhV69e2LZtG9RqtT5iJSIyaXlKFb764wrClxzG0YT7sDA3wwfh9bBzXBs0r+Nk6PCIiIzWc01ac3Z2RqtWrXDp0iVcvnwZFy5cwLBhw1C9enWsXr0aYWFhOgqTiMi0nbj2ANO2XMDVe1kAgDZ+NTC3dxBq17AxcGRERMZP6zu8AHD37l0sWrQIDRo0QFhYGNLT0/H7778jMTERt2/fRr9+/TB06FBdx0pEZHIeZuXjg1/+wqvfnMDVe1moYSvHl682xroRL7LYJSLSEa3v8Pbs2RN79uxBvXr18Oabb2LIkCFwdHTUbLeyssL777+PL774QqeBEhGZEiEENp29icid8XiUrQAADGruhY86B6CatczA0RERmRatC15nZ2ccOnQILVq0KLWPm5sbEhMTnyswIiJTlZCaiWlbLuBk4kMAQICrHeb3DUZTbwcDR0ZEZJq0LnhDQ0PRpEmTYu35+fn48ccfMWTIEEgkEnh7e+skQCIiU5GrUGHZgQR8fegqFCoBK5kUEzr6YURrH8ik5RphRkREZaD1b9jhw4cjLS2tWHtGRgaGDx+uk6CIiEzN0Sv30WXJYfxnfwIUKoH2Ac7YO7Et3g71ZbFLRKRnWt/hFUKUuOD5zZs3Ua1aNZ0ERURkKu5l5GHejjhsi74NAHCxlyOiZwN0CXLlwyOIiCpImQvekJAQSCQSSCQSdOjQAebm/99VpVIhMTERXbp00UuQRETGRq0W+PH0DXyyKx7puUpIJMDQFrXxfng92FlyUhoRUUUqc8Hbp08fAEB0dDQ6d+4MW1tbzTYLCwvUrl0bL7/8ss4DJCIyNhdT0jFtSwzOXn8EAAjysEdk32A0rFXdsIEREVVRZS54Z82aBQCoXbs2Bg4cCEtLS70FRURkjLLzlfjyjytYdSQRSrWAjYUU74f7Y0gLb5hznC4RkcFoPYaXD5QgIiruwMVUzNgWg5uPcgAAnRu4IKJXA7hVszJwZEREVKaC19HREZcvX0aNGjXg4ODw1IkWDx8+1FlwRESV3d30XMz+LRY7L6QAADyqW2F2rwboGOhi4MiIiKhQmQreL774AnZ2dpp/c2YxEVV1KrXA+hPX8dmeS8jMU0JqJsHI1j4Y38EPNnKtPzwjIiI9KtNv5X8PYxg2bJi+YiEiMgoxt9IwdcsF/H2zYE3yxp7VEdk3GIHu9gaOjIiISlKmgjc9Pb3MB7S35y98IjJNmXlKLN57GWuOJUItADtLc3zYJQCDXvSC1IyffBERVVZlKnirV6/+zGEMhQ+kUKlUOgmMiKgy2RObgojtsbiTlgsA6NHQDTN7BMLZnivWEBFVdmUqeA8cOKDvOIiIKqVbj3Mwa1ss9sXfBQB4OlphXp9ghNaraeDIiIiorMpU8IaGhuo7DiKiSkWpUmPNsSQsjrqM7HwVzM0keKttHbzX3g9WFlJDh0dERFooU8H7999/IygoCGZmZvj777+f2rdhw4Y6CYyIyFCibzzG1M0XEHenYP7CC7UdML9vMOq52Bk4MiIiKo8yFbyNGzdGSkoKnJ2d0bhxY0gkEgghivXjGF4iMmY5SmD27/H44dQNCAFUs5JharcAvNLUE2aclEZEZLTKVPAmJiaiZs2amn8TEZkSIQR2XkhBZLQU6YobAIB+IR6Y2r0+atjKDRwdERE9rzIVvN7e3iX+m4jI2N14mI3pW2Nw6PI9ABL4OFljft9gtKxbw9ChERGRjpTrcUCXLl3Cf/7zH8THx0MikSAgIADvvfce/P39dR0fEZFeKFRqrDxyDV/9cQW5CjVkUgk6uCmxaEQL2FpzqTEiIlNipu0OmzZtQlBQEM6ePYtGjRqhYcOGOHfuHIKCgvDLL7/oI0YiIp06k/QQ3b86goW7LyFXoUaLOk74fWxLdPUUkMu4AgMRkanR+g7vhx9+iClTpmDOnDlF2mfNmoWPPvoIr7zyis6CIyLSpcfZ+fh090VsPFUwTtfRxgLTu9dH3xAPKJVKXDRwfEREpB9aF7wpKSkYMmRIsfbBgwfjs88+00lQRES6JITA1uhbmPd7PB5k5QMABjbzxMddA+BgY2Hg6IiISN+0LnjDwsJw5MgR1K1bt0j70aNH0aZNG50FRkSkC9fuZWLGthj8mfAAAODnbIv5fYPxoo+jgSMjIqKKUqaCd/v27Zp/9+rVCx999BHOnj2Ll156CQBw4sQJ/PLLL5g9e7Z+oiQi0lKeUoWvD17DsoMJyFeqITc3w7gOfnizTR1YmGs9fYGIiIxYmQrePn36FGtbvnw5li9fXqRt7NixeOedd3QSGBFReR2/+gDTtl7AtXtZAIC29Wpibu8G8HayMXBkRERkCGUqeNVqtb7jICJ6bg8y8zB/Zzw2n7sFAKhhK8fMnoHo2dANEgmflEZEVFWVax1eIqLKRAiBX87cROSueDzOVkAiAV5v7oXJnQNQzUpm6PCIiMjAylXwZmVl4dChQ0hOTkZ+fn6RbePGjdNJYEREZXHlbgambYnBqaSHAIAAVztE9gtGEy8HA0dGRESVhdYF7/nz59GtWzdkZ2cjKysLjo6OuH//PqytreHs7KzzgjcjIwMzZszAli1bkJqaipCQEHz55Zd44YUXSt3nhx9+wMKFC3HlyhVUq1YNXbp0waJFi+Dk5KTT2IjIcHIVKizdn4D/Hr4KhUrASibFxE5+GN7KBzIpJ6UREdH/af1XYeLEiejZsycePnwIKysrnDhxAtevX0fTpk2xaNEinQc4atQoREVF4fvvv8eFCxcQHh6Ojh074tatWyX2P3r0KIYMGYKRI0ciNjYWv/zyC06fPo1Ro0bpPDYiMozDl++h85LDWHogAQqVQMf6zoia1BZvtfVlsUtERMVo/ZchOjoa77//PqRSKaRSKfLy8uDp6YmFCxdi6tSpOg0uJycHv/76KxYuXIi2bduibt26iIiIgI+PD1asWFHiPidOnEDt2rUxbtw4+Pj4oHXr1nj77bdx5swZncZGRBUvNSMX4zaex5DvTuH6g2y42lvi68FNsXJIM9RysDZ0eEREVElpPaRBJpNpZju7uLggOTkZ9evXR7Vq1ZCcnKzT4JRKJVQqFSwtLYu0W1lZ4ejRoyXu07JlS0ybNg07d+5E165dkZqaik2bNqF79+6lnicvLw95eXma1+np6QAAhUIBhUKhgyt5usJzVMS5SD+YQ/1SqwV+PHMTi6KuICNXCTMJ8MZLXpjQoS5s5eZQKpXPfQ7m0Pgxh8aN+TN+FZ1Dbc4jEUIIbQ4eHh6OYcOGYdCgQXjnnXdw/vx5jBs3Dt9//z0ePXqEkydPah3w07Rs2RIWFhbYsGEDXFxcsHHjRgwZMgR+fn64dOlSifts2rQJw4cPR25uLpRKJXr16oVNmzZBJit5tnZERESJD83YsGEDrK1514jIkG5lAT9fkyIps+A/2p42AgPrqOBpa+DAiIjIoLKzszFo0CCkpaXB3t7+qX21LnjPnDmDjIwMtGvXDvfu3cPQoUNx9OhR1K1bF6tXr0ajRo2eK/gnXb16FSNGjMDhw4chlUrRpEkT1KtXD+fOnUNcXFyx/nFxcejYsSMmTpyIzp07486dO5g8eTJeeOEFrFq1qsRzlHSH19PTE/fv33/mN1AXFAoFoqKi0KlTp1KLcqrcmEPdy85X4j8HrmH1setQqQVs5FJM7FAXg5t7QWqm+zV1mUPjxxwaN+bP+FV0DtPT01GjRo0yFbxaD2lo1qyZ5t81a9bEzp07tY9QC76+vjh06BCysrKQnp4ONzc3DBw4ED4+PiX2X7BgAVq1aoXJkycDABo2bAgbGxu0adMG8+bNg5ubW7F95HI55HJ5sXaZTFahb7qKPh/pHnOoG3/E38XMbbG49TgHANA1yBWzejaAazXLZ+z5/JhD48ccGjfmz/hVVA61OUe5HzyRmpqKS5cuQSKRwN/fHzVr1izvocrExsYGNjY2ePToEfbs2YOFCxeW2C87Oxvm5kUvSyqVAihYnJ6IKq+UtFzM/i0Wu2JSAAAe1a0wp3cDdKjvYuDIiIjImGld8Kanp2Ps2LH48ccfoVKpABQUlAMHDsSyZctQrVo1nQa4Z88eCCHg7++PhIQETJ48Gf7+/hg+fDgAYMqUKbh16xbWrVsHAOjZsyfefPNNrFixQjOkYcKECXjxxRfh7u6u09iISDdUaoF1x5Pw+d7LyMxTQmomwajWPhjf0Q/WFnwgJBERPR+t/5KMGjUK0dHR+P3339GiRQtIJBIcO3YM48ePx5tvvomff/5ZpwGmpaVhypQpuHnzJhwdHfHyyy9j/vz5mtvYd+7cKbI6xLBhw5CRkYGlS5fi/fffR/Xq1dG+fXt8+umnOo2LiHTjws00TN1yARdupQEAQryqY36fYAS663/8PBERVQ1aF7w7duzAnj170Lp1a01b586dsXLlSnTp0kWnwQHAgAEDMGDAgFK3r1mzpljbe++9h/fee0/nsRCR7mTmKfH53ktYeywJagHYWZrjoy4BGPSiF8z0MCmNiIiqLq0LXicnpxKHLVSrVg0ODnx2PRE9nRACe2JTELE9DinpuQCAXo3cMb1HfTjb6X9SGhERVT1aF7zTp0/HpEmTsG7dOs2KBykpKZg8eTJmzJih8wCJyHTcfJSNiO2x2BefCgDwcrTGvD5BaFtPv5NeiYioaitTwRsSEqJ5uhoAXLlyBd7e3vDy8gIAJCcnQy6X4969e3j77bf1EykRGS2FSo3Vfybii6gryFGoIJNK8HZbX7zbvi4sZVJDh0dERCauTAVvnz599BwGEZmqc8mPMHXzBVxMyQAAvFjbEfP7BsHPxc7AkRERUVVRpoJ31qxZ+o6DiExMWo4Cn+25iB9OJkMIoLq1DFO71kf/prU4KY2IiCpUuRe4PHv2LOLj4yGRSBAYGIiQkBBdxkVERkoIgd/+voO5v8fhXkbBI7tfblILU7sFwMm2+BMNiYiI9E3rgjc1NRWvvvoqDh48iOrVq0MIgbS0NLRr1w4//vij3p+4RkSV1/UHWZi+NQZHrtwHANSpYYN5fYPQ0reGgSMjIqKqzEzbHd577z2kp6cjNjYWDx8+xKNHjxATE4P09HSMGzdOHzESUSWXr1Rj2YEEhH9xGEeu3IeF1AwTO9bDrgltWOwSEZHBaX2Hd/fu3di3bx/q16+vaQsMDMSyZcsQHh6u0+CIqPI7lfgQ07ZcwJXUTABAS18nzOsThDo1bQ0cGRERUQGtC161Wq15rO+/yWQyqNVqnQRFRJXfo6x8fLLrIn46cwMA4GRjgek96qNPY48iyxgSEREZmtYFb/v27TF+/Hhs3LgR7u7uAIBbt25h4sSJ6NChg84DJKLKRQiBLedvYd6OeDzMygcAvPqCJz7uGoDq1hYGjo6IiKg4rQvepUuXonfv3qhduzY8PT0hkUiQnJyM4OBgrF+/Xh8xElElce1eJqZvjcGxqw8AAPVcbDG/bzBeqO1o4MiIiIhKp3XB6+npiXPnziEqKgoXL16EEAKBgYHo2LGjPuIjokogT6nCioNXsfzAVeSr1JCbm2FcBz+82aYOLMy1nvtKRERUobQqeJVKJSwtLREdHY1OnTqhU6dO+oqLiCqJ41cfYNrWC7h2LwsA0LZeTczrHQQvJ2sDR0ZERFQ2WhW85ubm8Pb2hkql0lc8RFRJPMzKR+TOeGw6exMAUMNWjlk9A9GjoRsnpRERkVHRekjD9OnTMWXKFKxfvx6Ojhy3R2RqhBD49dwtzN8Rh0fZCgDA68298GGXAFSzKr5CCxERUWWndcH71VdfISEhAe7u7vD29oaNjU2R7efOndNZcERUsRJSMzF96wWcuPYQABDgaof5fYPR1NvBwJERERGVn9YFb+/evflxJpGJyVWosPzgVXx9sGBSmqXMDBM61sPI1j6QSTkpjYiIjJvWBW9ERIQewiAiQzmWcB/TtsYg8X7BpLQw/5qY2zsIno6clEZERKahzLdusrOzMXbsWHh4eMDZ2RmDBg3C/fv39RkbEenRg8w8TPo5GoO+PYnE+1moaSfHskFNsHrYCyx2iYjIpJT5Du+sWbOwZs0avP7667C0tMTGjRsxevRo/PLLL/qMj4h0TAiBX87cROSueDzOVkAiAQY398bkLv6wt+SkNCIiMj1lLng3b96MVatW4dVXXwUADB48GK1atYJKpYJUKtVbgESkOwmpGZi6JQanEv8/KW1Bv2CEeHFSGhERma4yF7w3btxAmzZtNK9ffPFFmJub4/bt2/D09NRLcESkG7kKFZYfSMCKQ1ehUAlYyaSY2MkPw1txUhoREZm+Mhe8KpUKFhYWRXc2N4dSqdR5UESkO0ev3Mf0rReQ9CAbANA+wBlzejdALQeO0yUioqqhzAWvEALDhg2DXC7XtOXm5uKdd94pshbv5s2bdRshEZXL/cw8zN8Rjy3nbwEAXOzliOjZAF2CXLm0IBERVSllLniHDh1arG3w4ME6DYaInp9aLfDzmRtYsOsi0nIKJqUNbVEb74fXgx0npRERURVU5oJ39erV+oyDiHTgyt0MTN1yAaeTHgEAAt3sEdkvGI09qxs2MCIiIgPS+sETRFT55CpUWLo/Af89/P9Jae+H18OwlrVhzklpRERUxbHgJTJyhy/fw4xtMbj+z6S0jvWdMbt3EDyqWxk4MiIiosqBBS+RkbqXkYd5O+KwLfo2AMDV3hIRvRqgcwMXTkojIiL6Fxa8REZGrRb48fQNfLIrHum5SphJgKEta+P9cH/YyvmWJiIiehL/OhIZkUspBZPSzl4vmJQW5GGPyL7BaFirumEDIyIiqsRY8BIZgZx8Ff6z/wq+OXwNSrWAjYUUk8L9MbSFNyelERERPQMLXqJK7uClVMzYFoMbD3MAAOGBLojo1QDunJRGRERUJix4iSqp1IxczP09Hr/9VTApza1a4aQ0VwNHRkREZFxY8BJVMmq1wIZTyfh090Vk/DMpbVhLH0wKr8dJaUREROXAv55ElcjFlHRM3XwB55IfAwCCPaphQb9gBHlUM2xgRERERowFL1ElkJ2vxJd/XMGqI4maSWmTO/vjjRa1ITXjmrpERETPo9JP787IyMCECRPg7e0NKysrtGzZEqdPn37qPnl5eZg2bRq8vb0hl8vh6+uL7777roIiJtLOgUupCP/iMP57qGAFhi4NXLHv/VAMa+XDYpeIiEgHKv0d3lGjRiEmJgbff/893N3dsX79enTs2BFxcXHw8PAocZ8BAwbg7t27WLVqFerWrYvU1FQolcoKjpzo6VLTczH79zjs+PsOAMC9miVm9w5Cp0AXA0dGRERkWip1wZuTk4Nff/0V27ZtQ9u2bQEAERER2Lp1K1asWIF58+YV22f37t04dOgQrl27BkdHRwBA7dq1KzJsoqdSqwV+OJWMhbsuIiNPCamZBMNb1sbETvVgw0lpREREOlep/7oqlUqoVCpYWloWabeyssLRo0dL3Gf79u1o1qwZFi5ciO+//x42Njbo1asX5s6dCyurktctzcvLQ15enuZ1eno6AEChUEChUOjoakpXeI6KOBfpR1lzGH8nAzO2x+Gvm2kAgIYe9pjTKxAN3O0BCP4MGBDfh8aPOTRuzJ/xq+gcanMeiRBC6DGW59ayZUtYWFhgw4YNcHFxwcaNGzFkyBD4+fnh0qVLxfp36dIFBw8eRMeOHTFz5kzcv38fY8aMQfv27UsdxxsREYHZs2cXa9+wYQOsra11fk1U9eSpgN03zHDwjgRqSCCXCvTwVKO1qwCH6RIREWkvOzsbgwYNQlpaGuzt7Z/at9IXvFevXsWIESNw+PBhSKVSNGnSBPXq1cO5c+cQFxdXrH94eDiOHDmClJQUVKtWsJTT5s2b0b9/f2RlZZV4l7ekO7yenp64f//+M7+BuqBQKBAVFYVOnTpBJpPp/Xyke0/L4YFL9zD793jcepwLAOjSwAXTuvnD1d6ypEORgfB9aPyYQ+PG/Bm/is5heno6atSoUaaCt1IPaQAAX19fHDp0CFlZWUhPT4ebmxsGDhwIHx+fEvu7ubnBw8NDU+wCQP369SGEwM2bN+Hn51dsH7lcDrlcXqxdJpNV6Juuos9HuvfvHN5Nz8Xs32Kx80IKAMCjuhXm9G6ADvU5Ka0y4/vQ+DGHxo35M34VlUNtzlHpC95CNjY2sLGxwaNHj7Bnzx4sXLiwxH6tWrXCL7/8gszMTNja2gIALl++DDMzM9SqVasiQ6YqSqUWWH/iOj7bcwmZ/0xKG9XaB+M7+sHawmjeckRERCaj0q/Du2fPHuzevRuJiYmIiopCu3bt4O/vj+HDhwMApkyZgiFDhmj6Dxo0CE5OThg+fDji4uJw+PBhTJ48GSNGjCh10hqRrsTdSUe/5X9i1vZYZOYp0dizOn57tzWmdKvPYpeIiMhAKv1f4LS0NEyZMgU3b96Eo6MjXn75ZcyfP19zG/vOnTtITk7W9Le1tUVUVBTee+89NGvWDE5OThgwYECJS5gR6UpWnhJbk8xw+ORJqNQCdnJzfNjFH4Oae/PhEURERAZW6QveAQMGYMCAAaVuX7NmTbG2gIAAREVF6TEqov/bF3cXM7bF4E6aGQCB7g3dMLNHIFw4KY2IiKhSqPQFL1FldSctB7O3x2F3bMGkNEe5wMIBTdCxgbuBIyMiIqJ/Y8FLpCWVWmDd8SQs2nMJWfkqSM0kGNnKG355CQitV9PQ4REREdETWPASaSHmVhqmbrmAv/95UlqIV3VE9g1G3RpW2LkzwcDRERERUUlY8BKVQVaeEoujLmP1n4lQC8DO0hwfdQnAoBe9YGYm4aMwiYiIKjEWvETPsDc2BbO2x+JOWsGT0no2cseMHvXhbMdJaURERMaABS9RKW4/zkHE9ljsjbsLAPB0tMLc3kEI83c2cGRERESkDRa8RE9QqtRYe/w6Fu8tmJRmbibBW23r4L32frCykBo6PCIiItISC16if7lwMw1TtvyNmFvpAICm3g6I7BsMf1c7A0dGRERE5cWClwhARq4Cn++9jHXHk6AWgL2lOaZ0q4+BzTxhxielERERGTUWvFSlCSGwJ/YuIrbHIiW9YFJa78bumN49EDXt5AaOjoiIiHSBBS9VWbce52DWtljsiy+YlOblaI15fYLQlg+PICIiMikseKnKUarUWHMsCYujLiM7XwWZVIK32/ri3fZ1YSnjpDQiIiJTw4KXqpS/bjzGlM0XEHenYFLaC7UdML9vMOq5cFIaERGRqWLBS1VCRq4Ci/ZcwroT1yEEUM1KhqndAvBKU05KIyIiMnUseMmkCSGwOyYFEb/F4m56HgCgb4gHpnWvjxq2nJRGRERUFbDgJZN181E2Zm2LxR8XUwEAtZ2sMa9PMFr71TBwZERERFSRWPCSyVGq1Fj9Z8GktBxFwaS00aG+GNOOk9KIiIiqIha8ZFLOJz/C1C0xiP9nUtqLtR0R2S8IdZ05KY2IiKiqYsFLJiE9V4HPdl/C+pMFk9KqW8swtWt99G9ai5PSiIiIqjgWvGTUhBDYeSEFs3+LRWpGwaS0fk08MK1bfThxUhoRERGBBS8ZsRsPszFzWwwOXLoHAPCpYYP5fYLQsi4npREREdH/seAlo6NQqbHqaCKW7LuMXIUaFlIzjA7zxegwX05KIyIiomJY8JJROXv9EaZtuYCLKRkAgOY+jpjfNxh1nW0NHBkRERFVVix4ySik5Sjw2Z6L+OFkMoQAHKxlmNqtYFKaRMJJaURERFQ6FrxUqQkh8PvfdzDn9zjc+2dSWv+mtTC1W3042lgYODoiIiIyBix4qdJKfpCNGdticOhywaS0OjVtML9PMFr4Ohk4MiIiIjImLHip0lGo1Fh55Bq+3HcFecqCSWlj29XFO2F1IDfnpDQiIiLSDgteqlTOXn+IqZtjcOluwaS0FnWcMK9vEHxrclIaERERlQ8LXqoU0rIV+HTPRWw4mQwAcLSxwPTu9dE3xIOT0oiIiOi5sOAlgxJCYPtftzH39zjcz8wHAAxoVgtTutaHAyelERERkQ6w4CWDuf4gC9O3xuDIlfsAAN+aNpjfNxgv1eGkNCIiItIdFrxU4fKVBZPSvvrjn0lp5mZ4r11dvBXKSWlERESkeyx4qUKdTnqIqZsv4EpqJgCgVV0nzOsTDJ8aNgaOjIiIiEwVC16qEI+z8/HJrov48fQNAICTjQVm9AhE78bunJRGREREesWCl/RKCIFt0QWT0h5kFUxKe/UFT3zcNQDVrTkpjYiIiPSPBS/pTdL9gklpRxMKJqX5Odtift9gvOjjaODIiIiIqCphwUs6l69U47+HruI/BxKQr1RDbm6GcR388GabOrAwNzN0eERERFTFsOAlnTp57QGmbY1Bwj+T0tr41cDc3kGozUlpREREZCCV/nZbRkYGJkyYAG9vb1hZWaFly5Y4ffp0mfb9888/YW5ujsaNG+s3SMKjrHx8uOkvDPzmBBJSM1HD1gJfvtoY60a8yGKXiIiIDKrS3+EdNWoUYmJi8P3338Pd3R3r169Hx44dERcXBw8Pj1L3S0tLw5AhQ9ChQwfcvXu3AiOuWoQQ2HL+FubtiMfDfyalvfaiFz7uEoBq1jIDR0dERERUye/w5uTk4Ndff8XChQvRtm1b1K1bFxEREfDx8cGKFSueuu/bb7+NQYMGoUWLFhUUbdVz7V4mXv/2JCb9/BceZuWjnostNr3TAgv6BbPYJSIiokqjUt/hVSqVUKlUsLS0LNJuZWWFo0ePlrrf6tWrcfXqVaxfvx7z5s175nny8vKQl5eneZ2eng4AUCgUUCgU5Yy+7ArPURHn0oU8pRrfHEnE14cTNZPS3mvni+EtvWFhbmY016FLxpZDKo45NH7MoXFj/oxfRedQm/NIhBBCj7E8t5YtW8LCwgIbNmyAi4sLNm7ciCFDhsDPzw+XLl0q1v/KlSto3bo1jhw5gnr16iEiIgJbt25FdHR0qeeIiIjA7Nmzi7Vv2LAB1tbWurwco5eQBvx0TYrU3IKHRQRUU+OVOmrUsHzGjkREREQ6lJ2djUGDBiEtLQ329vZP7Vup7/ACwPfff48RI0bAw8MDUqkUTZo0waBBg3Du3LlifVUqFQYNGoTZs2ejXr16ZT7HlClTMGnSJM3r9PR0eHp6Ijw8/JnfQF1QKBSIiopCp06dIJNVzqEAD7Py8emey9gcdxsAUMPWAtO6+qN7sCuflAbjyCE9HXNo/JhD48b8Gb+KzmHhJ/JlUekLXl9fXxw6dAhZWVlIT0+Hm5sbBg4cCB8fn2J9MzIycObMGZw/fx7vvvsuAECtVkMIAXNzc+zduxft27cvtp9cLodcLi/WLpPJKvRNV9HnKwshBH49dwvzd8ThUXbBRwevN/fCh10CUM2qcsVaGVTGHJJ2mEPjxxwaN+bP+FVUDrU5R6UveAvZ2NjAxsYGjx49wp49e7Bw4cJifezt7XHhwoUibcuXL8f+/fuxadOmEotkKt3Ve5mYtuUCTlx7CAAIcLXD/L7BaOrtYODIiIiIiMqu0he8e/bsgRAC/v7+SEhIwOTJk+Hv74/hw4cDKBiOcOvWLaxbtw5mZmYICgoqsr+zszMsLS2LtVPpchUqrDh4FSsOXkW+Sg1LmRkmdKyHka19IJNW6oU9iIiIiIqp9AVvWloapkyZgps3b8LR0REvv/wy5s+fr7mNfefOHSQnJxs4StNx7Op9TN8Sg2v3swAAYf41Mbd3EDwdOXmPiIiIjFOlL3gHDBiAAQMGlLp9zZo1T90/IiICERERug3KBD3Mysf8HfH49dxNAEBNOzkiejZAN05KIyIiIiNX6Qte0i8hBH45exORO+PxOFsBiQQY3Nwbk7v4w96SkwaIiIjI+LHgrcISUjMwdUsMTiX+f1Lagn7BCPHipDQiIiIyHSx4q6BchQrLDyRgxaGrUKgErGRSTOzkh+GtOCmNiIiITA8L3irmz4T7mL41Bon/TEprH+CMOb0boJYDJ6URERGRaWLBW0U8yMzD/B3x2Hz+FgDA2U6O2b0aoEsQJ6URERGRaWPBa+LUaoFfzt5A5M6LSMspmJQ25CVvvN+Zk9KIiIioamDBa8Ku3M3AtC0xOJVUMCkt0M0ekf2C0dizumEDIyIiIqpALHhNUK5ChaX7E/Dfw/+flPZ+eD0Ma1kb5pyURkRERFUMC14Tc+TKPUzfGoPrD7IBAB3rOyOiFyelERERUdXFgtdE3MvIw7wdcdgWfRsA4GpviYheDdC5gQsnpREREVGVxoLXyKnVAj+duYEFO+ORnquEmQQY0qI23g+vBztOSiMiIiJiwWvMLt/NwNTNF3Dm+iMAQJCHPSL7BqNhreqGDYyIiIioEmHBa4RyFSp89ccVfHP4GpRqAWsLKd4P98fQFt6clEZERET0BBa8RubQ5XuYsTUGyQ8LJqV1CnTB7F4N4F7dysCREREREVVOLHiNRGpGLub+Ho/f/iqYlOZWrXBSmquBIyMiIiKq3FjwVnJqtcDG08n4ZNdFZPwzKW1YSx9MCq8HWznTR0RERPQsrJgqsYsp6Zi6+QLOJT8GAAR7VMOCfsEI8qhm2MCIiIiIjAgL3kooJ1+FL/+4gm+PFExKs7GQ4oPO/hjSojakZlxTl4iIiEgbLHgrmQOXUjFjawxuPsoBAHRu4IKIXg3gVo2T0oiIiIjKgwVvJZGWD4z/6S/sjLkLAHCvZonZvYPQKdDFwJERERERGTcWvJXAhlM3sCBailzVXZhJgBGtfDCxUz3YcFIaERER0XNjRVUJxN3JQK5KgoYe9ojs15CT0oiIiIh0iAVvJTA53A/qB9cxd1hzWMotDB0OERERkUnhc2grgWpWMrR2FVyBgYiIiEgPWPASERERkUljwUtEREREJo0FLxERERGZNBa8RERERGTSWPASERERkUljwUtEREREJo0FLxERERGZNBa8RERERGTSWPASERERkUljwUtEREREJo0FLxERERGZNBa8RERERGTSWPASERERkUljwUtEREREJs3c0AFURkIIAEB6enqFnE+hUCA7Oxvp6emQyWQVck7SLebQ+DGHxo85NG7Mn/Gr6BwW1mmFddvTsOAtQUZGBgDA09PTwJEQERER0dNkZGSgWrVqT+0jEWUpi6sYtVqN27dvw87ODhKJRO/nS09Ph6enJ27cuAF7e3u9n490jzk0fsyh8WMOjRvzZ/wqOodCCGRkZMDd3R1mZk8fpcs7vCUwMzNDrVq1Kvy89vb2fJMbOebQ+DGHxo85NG7Mn/GryBw+685uIU5aIyIiIiKTxoKXiIiIiEwaC95KQC6XY9asWZDL5YYOhcqJOTR+zKHxYw6NG/Nn/CpzDjlpjYiIiIhMGu/wEhEREZFJY8FLRERERCaNBS8RERERmTQWvERERERk0ljw6sny5cvh4+MDS0tLNG3aFEeOHHlq/0OHDqFp06awtLREnTp18PXXXxfr8+uvvyIwMBByuRyBgYHYsmWLvsIn6D6HK1euRJs2beDg4AAHBwd07NgRp06d0uclVGn6eA8W+vHHHyGRSNCnTx8dR03/po8cPn78GGPHjoWbmxssLS1Rv3597Ny5U1+XUOXpI4dLliyBv78/rKys4OnpiYkTJyI3N1dfl1ClaZO/O3fuYNCgQfD394eZmRkmTJhQYj+D1TKCdO7HH38UMplMrFy5UsTFxYnx48cLGxsbcf369RL7X7t2TVhbW4vx48eLuLg4sXLlSiGTycSmTZs0fY4dOyakUqmIjIwU8fHxIjIyUpibm4sTJ05U1GVVKfrI4aBBg8SyZcvE+fPnRXx8vBg+fLioVq2auHnzZkVdVpWhj/wVSkpKEh4eHqJNmzaid+/eer6SqksfOczLyxPNmjUT3bp1E0ePHhVJSUniyJEjIjo6uqIuq0rRRw7Xr18v5HK5+OGHH0RiYqLYs2ePcHNzExMmTKioy6oytM1fYmKiGDdunFi7dq1o3LixGD9+fLE+hqxlWPDqwYsvvijeeeedIm0BAQHi448/LrH/hx9+KAICAoq0vf322+Kll17SvB4wYIDo0qVLkT6dO3cWr776qo6ipn/TRw6fpFQqhZ2dnVi7du3zB0xF6Ct/SqVStGrVSnz77bdi6NChLHj1SB85XLFihahTp47Iz8/XfcBUjD5yOHbsWNG+ffsifSZNmiRat26to6ipkLb5+7fQ0NASC15D1jIc0qBj+fn5OHv2LMLDw4u0h4eH49ixYyXuc/z48WL9O3fujDNnzkChUDy1T2nHpPLTVw6flJ2dDYVCAUdHR90ETgD0m785c+agZs2aGDlypO4DJw195XD79u1o0aIFxo4dCxcXFwQFBSEyMhIqlUo/F1KF6SuHrVu3xtmzZzXDwa5du4adO3eie/fueriKqqs8+SsLQ9Yy5no/QxVz//59qFQquLi4FGl3cXFBSkpKifukpKSU2F+pVOL+/ftwc3MrtU9px6Ty01cOn/Txxx/Dw8MDHTt21F3wpLf8/fnnn1i1ahWio6P1FTr9Q185vHbtGvbv34/XX38dO3fuxJUrVzB27FgolUrMnDlTb9dTFekrh6+++iru3buH1q1bQwgBpVKJ0aNH4+OPP9bbtVRF5clfWRiylmHBqycSiaTIayFEsbZn9X+yXdtj0vPRRw4LLVy4EBs3bsTBgwdhaWmpg2jpSbrMX0ZGBgYPHoyVK1eiRo0aug+WSqTr96BarYazszO++eYbSKVSNG3aFLdv38Znn33GgldPdJ3DgwcPYv78+Vi+fDmaN2+OhIQEjB8/Hm5ubpgxY4aOoyd91B2GqmVY8OpYjRo1IJVKi/1vJTU1tdj/agq5urqW2N/c3BxOTk5P7VPaMan89JXDQosWLUJkZCT27duHhg0b6jZ40kv+YmNjkZSUhJ49e2q2q9VqAIC5uTkuXboEX19fHV9J1aWv96CbmxtkMhmkUqmmT/369ZGSkoL8/HxYWFjo+EqqLn3lcMaMGXjjjTcwatQoAEBwcDCysrLw1ltvYdq0aTAz40hNXShP/srCkLUMfzJ0zMLCAk2bNkVUVFSR9qioKLRs2bLEfVq0aFGs/969e9GsWTPIZLKn9intmFR++sohAHz22WeYO3cudu/ejWbNmuk+eNJL/gICAnDhwgVER0drvnr16oV27dohOjoanp6eerueqkhf78FWrVohISFB858VALh8+TLc3NxY7OqYvnKYnZ1drKiVSqUQBZPwdXgFVVt58lcWBq1l9D4trgoqXMpj1apVIi4uTkyYMEHY2NiIpKQkIYQQH3/8sXjjjTc0/QuXYpk4caKIi4sTq1atKrYUy59//imkUqn45JNPRHx8vPjkk0+4LJke6SOHn376qbCwsBCbNm0Sd+7c0XxlZGRU+PWZOn3k70lcpUG/9JHD5ORkYWtrK959911x6dIl8fvvvwtnZ2cxb968Cr++qkAfOZw1a5aws7MTGzduFNeuXRN79+4Vvr6+YsCAARV+faZO2/wJIcT58+fF+fPnRdOmTcWgQYPE+fPnRWxsrGa7IWsZFrx6smzZMuHt7S0sLCxEkyZNxKFDhzTbhg4dKkJDQ4v0P3jwoAgJCREWFhaidu3aYsWKFcWO+csvvwh/f38hk8lEQECA+PXXX/V9GVWarnPo7e0tABT7mjVrVgVcTdWjj/fgv7Hg1T995PDYsWOiefPmQi6Xizp16oj58+cLpVKp70upsnSdQ4VCISIiIoSvr6+wtLQUnp6eYsyYMeLRo0cVcDVVj7b5K+lvnLe3d5E+hqplJP8ESERERERkkjiGl4iIiIhMGgteIiIiIjJpLHiJiIiIyKSx4CUiIiIik8aCl4iIiIhMGgteIiIiIjJpLHiJiIiIyKSx4CUiIiIik8aCl4jISCQlJUEikSA6OlonfdesWYPq1asXafvmm2/g6ekJMzMzLFmyRKv4Ll26BFdXV2RkZDyzb15eHry8vHD27FmtzkFEVB4seImIdGzYsGGQSCSQSCQwNzeHl5cXRo8ejUePHhk6tCIGDhyIy5cva16np6fj3XffxUcffYRbt27hrbfeQlhYGCZMmFCm402bNg1jx46FnZ3dM/vK5XJ88MEH+Oijj8obPhFRmbHgJSLSgy5duuDOnTtISkrCt99+i99++w1jxowxdFhFWFlZwdnZWfM6OTkZCoUC3bt3h5ubG6ytrct8rJs3b2L79u0YPnx4mfd5/fXXceTIEcTHx2sVNxGRtljwEhHpgVwuh6urK2rVqoXw8HAMHDgQe/fuLdJn9erVqF+/PiwtLREQEIDly5cX2X7q1CmEhITA0tISzZo1w/nz54tsf/ToEV5//XXUrFkTVlZW8PPzw+rVq4v0uXbtGtq1awdra2s0atQIx48f12z795CGNWvWIDg4GABQp04dSCQSDBs2DIcOHcKXX36puWOdlJRU4vX+/PPPaNSoEWrVqqVpCwsL0+z376/CYzg5OaFly5bYuHFjmb+vRETlYW7oAIiITN21a9ewe/duyGQyTdvKlSsxa9YsLF26FCEhITh//jzefPNN2NjYYOjQocjKykKPHj3Qvn17rF+/HomJiRg/fnyR486YMQNxcXHYtWsXatSogYSEBOTk5BTpM23aNCxatAh+fn6YNm0aXnvtNSQkJMDcvOiv/4EDB8LT0xMdO3bEqVOn4OnpCSsrK1y+fBlBQUGYM2cOAKBmzZolXuPhw4fRrFmzIm2bN29Gfn6+5vXYsWMRGxsLFxcXTduLL76II0eOaPHdJCLSHgteIiI9+P3332FrawuVSoXc3FwAwOLFizXb586di88//xz9+vUDAPj4+CAuLg7//e9/MXToUPzwww9QqVT47rvvYG1tjQYNGuDmzZsYPXq05hjJyckICQnRFJq1a9cuFscHH3yA7t27AwBmz56NBg0aICEhAQEBAUX6WVlZwcnJCUBBUevq6goAsLCwgLW1teZ1aZKSktC0adMibY6Ojpp/f/HFF9i/fz9OnjwJKysrTbuHh0epd42JiHSFBS8RkR60a9cOK1asQHZ2Nr799ltcvnwZ7733HgDg3r17uHHjBkaOHIk333xTs49SqUS1atUAAPHx8WjUqFGRcbQtWrQoco7Ro0fj5Zdfxrlz5xAeHo4+ffqgZcuWRfo0bNhQ8283NzcAQGpqarGC93nl5OTA0tKyxG27du3Cxx9/jN9++w316tUrss3KygrZ2dk6jYWI6Ekcw0tEpAc2NjaoW7cuGjZsiK+++gp5eXmYPXs2AECtVgMoGNYQHR2t+YqJicGJEycAAEKIZ56ja9euuH79OiZMmIDbt2+jQ4cO+OCDD4r0+fcwColEUuT8ulSjRo0SV6GIi4vDq6++ik8++QTh4eHFtj98+LDUYRJERLrCgpeIqALMmjULixYtwu3bt+Hi4gIPDw9cu3YNdevWLfLl4+MDAAgMDMRff/1VZExuYTH8bzVr1sSwYcOwfv16LFmyBN98841O47awsIBKpXpmv5CQEMTFxRVpe/DgAXr27Il+/fph4sSJJe4XExODkJAQncRKRFQaFrxERBUgLCwMDRo0QGRkJAAgIiICCxYswJdffonLly/jwoULWL16tWac76BBg2BmZoaRI0ciLi4OO3fuxKJFi4occ+bMmdi2bRsSEhIQGxuL33//HfXr19dp3LVr18bJkyeRlJSE+/fvl3p3uHPnzjh+/HiR4rhfv36wsrJCREQEUlJSNF//7nPkyJES7/wSEekSC14iogoyadIkrFy5Ejdu3MCoUaPw7bffapYDCw0NxZo1azR3eG1tbfHbb78hLi4OISEhmDZtGj799NMix7OwsMCUKVPQsGFDtG3bFlKpFD/++KNOY/7ggw8glUoRGBiImjVrIjk5ucR+3bp1g0wmw759+zRthw8fRmxsLGrXrg03NzfN140bNwAAx48fR1paGvr376/TmImIniQRZRkoRkRE9AzLly/Htm3bsGfPnjL1f+WVVxASEoKpU6fqOTIiquq4SgMREenEW2+9hUePHiEjI+OZjxfOy8tDo0aNSh3bS0SkS7zDS0REREQmjWN4iYiIiMikseAlIiIiIpPGgpeIiIiITBoLXiIiIiIyaSx4iYiIiMikseAlIiIiIpPGgpeIiIiITBoLXiIiIiIyaSx4iYiIiMik/Q+RtVe7ogLa2AAAAABJRU5ErkJggg==",
      "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": "9e9fd968",
   "metadata": {},
   "source": [
    "作业2:某旋涡星系和椭圆星系，其观测到的I波段的视星等都是18等；半光度半径内的表面亮度都是20mag/arcsec^2；其中旋涡星系（图像椭率为\n",
    "0.5 ）观测得到的HI的速度展宽是200km/s，椭圆星系观测到的中心速度弥散度是200km/s，请问这两个星系的距离分别是多少？\n",
    "• 旋涡星系的 Tully-Fisher关系以及椭圆星系的Fundamental Plane的内秉弥散在0.2个星等左右，请问以上距离计算的误差是多少？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8a95ba49",
   "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"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "base"
  },
  "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}