{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于作业5中的星系，依据fits图像画出其一维的表面亮度轮廓，并计算其星等<p>\n",
    "注意流量单位\n",
    "• 一维轮廓：\n",
    "✓椭圆测光：椭圆是等面亮度轮廓，半径沿着长轴方向\n",
    "✓圆测光: 圆孔径内的平均表面亮度\n",
    "• 星等：\n",
    "✓积分（求和）到某个等亮度半径处\n",
    "✓根据面亮度轮廓模型，积分到无穷远处"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "星等：14.479804983781277\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArQAAAIqCAYAAADVUtnpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACTSUlEQVR4nOzdd3gUVdvH8e+mkARIQgsl0kFApEPoJfSOFBFF6Spi6HaUKgoiPiiCKIggTXpHkEiVXiNNkI4IiLQEEggp8/4xbxJiAiTLJptkf5/r2uuZOTM7c++e5PHm5Mx9LIZhGIiIiIiIpFNO9g5ARERERORJKKEVERERkXRNCa2IiIiIpGtKaEVEREQkXVNCKyIiIiLpmhJaEREREUnXlNCKiIiISLqmhFZERERE0jUltCIiIiKSrimhFUmjTp48Sd++fSldujRZsmTB3d2d/Pnz4+fnR9++fVmyZEmqx3Tr1i0CAgIoVKgQmTJlwmKx4O/vn+px2FrhwoWxWCwJXlmzZqV8+fJ88MEHXL9+3errz5w5E4vFQvfu3W0S77lz57BYLBQuXNgm13NUifW5h4cHxYoVo0ePHhw6dChV4zl27Bht27Yld+7cODs7Y7FYGDFiBAD+/v5YLBY2b94c7z0jRoyId56Io3KxdwAiktDSpUvp3Lkz4eHh5MyZk1q1auHj48PNmzcJCgpi8uTJzJ8/nw4dOqRqXK+//jqLFi2icOHCtG/fHnd3d0qVKpWqMaSkWrVqUbx4cQCio6O5dOkSO3bsYOzYscyaNYvffvuNokWL2jnKJ1O4cGHOnz/P2bNnlRD/v6ZNm5I3b14Arl69yt69e5k5cyZz585lzpw5vPDCCykeQ2hoKC1btuTcuXNUqVKFpk2b4uzsTIUKFVL83iIZgRJakTTmn3/+oVu3boSHh/PWW28xevRo3N3d452zf/9+Fi9enKpxRUREsGzZMtzd3fn999/x8vJK1funhldffTXBKOqVK1eoV68ef/75J++++65V33u7du2oXr063t7eNopUbOn999+P95eG4OBgOnbsSGBgIK+99hqNGzcme/bsKRrD3r17OXfuHDVr1mT79u0Jjs+aNYuwsDAKFiyYonGIpFeaciCSxqxevZo7d+7g6+vL+PHjEySzAJUrV2bMmDGpGtfly5eJjIwkT548GTKZfZi8efPyzjvvALBhwwarruHt7U2pUqXIly+fLUOTFOLt7c3UqVMBCAkJ4Zdffknxe164cAGAp59+OtHjBQsWpFSpUmTOnDnFYxFJj5TQiqQx//zzDwA+Pj7Jfm/MPMCHedg8vAfbf/vtN1q3bo2Pjw9OTk6x8z8LFSoEwPnz5+PNOYy51r///svEiRNp0aIFRYoUwcPDAy8vL6pUqcJnn33GvXv3HhpXWFgYX375JbVr1yZ79uy4ublRqFAhWrduzbx58xJ9z+LFi2nWrBk+Pj5kypSJp556ildeeYVjx44l70tLgpg/R0dGRiY49rjvDh4/h3bFihXUqVMHT09PvL29qVevHmvWrEnSXFnDMJg6dSqVK1cmS5YseHt706RJE3bu3BnvvJgYzp8/D0CRIkUS7cfNmzfHzo2OiIjgs88+49lnn8XDw4OcOXPSvn17/vjjj4fGc/PmTYYPH06FChXw9PQkc+bMlC1bltGjRxMWFpbg/OjoaKZOnUqtWrXIli0brq6u5M6dm/Lly9OvXz/OnTsX7/zLly8zYMAASpQogbu7O5kzZ6ZAgQI0bNiQ8ePHPzSu5CpcuDA5cuQAiI3hwf6Iiorif//7HxUrViRr1qwJfu9++eUXWrVqRe7cucmUKRO+vr506tSJffv2xTsv5vvu1q0bAD/++GO8fonxsN/dx/nzzz/p3bs3xYoVw93dHW9vb+rWrcucOXOS+Y2IpG2aciCSxsT8SfHIkSNs2LCBhg0bptq9Fy1axLfffkupUqVo1KgRN27cwM3NjW7dunHnzh2WLFlClixZeP7552PfE5Ps/fLLLwwYMICnnnqK4sWLU716df799192797N+++/z4oVK9i0aRNubm7x7vnXX3/RrFkzjh07RubMmalVqxY5c+bk77//5rfffuPw4cN07tw59vzIyEhefvllFi5ciJubG5UrV+app57izz//ZO7cuSxdupSlS5fSrFkzm30ve/bsAeDZZ59N9nf3OOPGjeO9994DoFq1ahQtWpRTp07RqlUr3n333ce+v0ePHsybN486derQqlUrgoKCCAwMZOvWrWzZsoVq1aoBULx4cbp168bixYsJDQ2lQ4cOZM2aNfY6Mf0YIyIighYtWrBjxw7q1q3LM888w549e1i2bBmbNm3i4MGDCRLtY8eO0axZM/766y/y5ctH7dq1cXV1Zc+ePQwdOpQlS5awefPmeFMvXn31VWbMmIG7uzu1a9fGx8eHGzducObMGSZNmkTDhg1j73PlyhWqVKnCpUuXKFiwIM2aNcPd3Z1Lly4RFBTE/v37efvttx/7nSVFdHQ0oaGhAAn60TAM2rdvz7p166hTpw7PPPMMR48ejT0+dOhQRo8ejcVioWbNmhQsWJA//viDhQsXsmTJEqZOnUrPnj1jv/du3bpx6tQptm/fTrFixahdu7ZNPsOiRYvo2rUr9+7do1SpUrRo0YLg4GB2795Nly5d2LhxIz/88INN7iVid4aIpCm3b982nnrqKQMwLBaL4e/vb3z88cfGmjVrjKtXrz7yvYDxqF/revXqGYCxadOmRNsBY/LkyYm+9+zZswZgFCpUKNHjx44dM3bu3Jmg/caNG0aTJk0MwBg3bly8Y1FRUUaVKlUMwGjSpEmCz3f37l1jzZo18dqGDBliAEa1atWMM2fOxDu2aNEiw9nZ2ciePbtx8+bNRONMTKFChQzAmDFjRrzYLl68aHz99deGm5ub4ezsbKxatSrBe5Py3c2YMcMAjG7dusVrP3DggOHs7Gw4OzsbS5cujXds4cKFhpOTU6LfeUxfxBw7ceJE7LHIyEijZ8+esd/pwz7r2bNnE41106ZNsdeuWLGicfny5dhjd+/eNZo2bWoAxuuvvx7vfWFhYUaxYsUMwPjoo4+M8PDw2GOhoaHGSy+9ZABGjx49YtvPnz9vAEb+/Pnj3SfGsWPHjPPnz8fujxw5Mvbe0dHR8c69f/++8euvvyb6mR4m5nP+9/fBMAxj9erVscc3btxoGEb87z1//vzxvvcYa9euNQDD3d3dWL9+fbxj33//vQEYrq6uxpEjR+Ide9jPSIyH/e4OHz7cAIzhw4fHaz906JDh5uZmuLu7G0uWLIl37Ny5c0bZsmUNwPjxxx8TvZ9IeqOEViQNOn78uFGtWrXY/3g++KpQoYIxZcoUIzIyMsH7njShbdCgwUPf+7iE9lFOnDhhAIafn1+89uXLlxuAkS9fPuP27duPvc7169cNDw8Pw93d3bh48WKi57z55psGYHz99ddJji8myXvYy8/Pz9i2bVui703Kd/ewZCUm8XzppZcSfd/zzz//2IR25cqVCd53+fJlAzDc3NyM+/fvJ/pZH5fQWiwWIygoKMHxXbt2GYBRtGjReO1TpkwxAKNVq1aJXvf27dtG7ty5DRcXF+PGjRuGYRjGnj17DMBo06ZNou/5r5i+/W/yb63EEtp///3XmDdvnpE7d+7Y37eoqCjDMOJ/77NmzUr0mg0bNjQAY/DgwYkeb9WqlQEYr732Wrx2Wye0nTp1MgBj/PjxiV4v5ruvXLlyosdF0hvNoRVJg0qWLMmuXbvYvXs3w4YNo2nTprFzaoOCgujTpw/NmjXj/v37Nr3vg1MJrBEVFcWGDRv4+OOPefPNN+nRowfdu3fnk08+AeDEiRPxzl+3bh0AnTt3jvfn74fZtGkTd+/epVatWjz11FOJnhPztPqOHTuSHX+tWrXo1q1b7Ktly5YUKFCAvXv3MmjQIE6ePPnQ91rz3W3ZsgWAl19+OdHjD2uP4eLikujUirx585I9e3bCw8Otrp9bsGBBypcvn6D9mWeeAeDvv/+O175mzRoAOnXqlOj1smbNSpUqVYiMjGTv3r0AlCpVCk9PT37++Wc++eQTzp49+8iYqlatCphVCZYuXcqdO3eS96Eeon79+rFzVn18fOjcuTNXr16lUqVKLF++HCenhP+pTKxkXmRkZGyFgofNl+7Vqxdg/iynlOjoaNauXQs8vD+qVKlC1qxZOXjw4CPnt4ukF5pDK5KGVa1aNfY/4oZhcPDgQT7//HPmz5/Pr7/+yldffRX7BL4tPEld0pMnT9KuXbt4cwn/KyQkJN5+zANKSa1le+bMGcCsNvCoh9/AfEgtuRIr2xUZGcmwYcMYM2YM9erV48SJE3h6eiZ4rzXf3cWLFx/53sddM1++fLi6uiZ6zMvLi5s3b1qdrDysPFRMhYvw8PB47TF906VLF7p06fLIa8f0jaenJzNmzKBHjx589NFHfPTRR+TLl4/q1avTrFmzBP/Q6dKlC4GBgcydO5cOHTrg7OxM6dKlqV27Ns8//zwNGjSw6rM+WIfWzc0NX19f6tSpE5vo/lfu3LkTrTZw/fr12O+7SJEiid6rWLFiQMJ/ENjS9evXY3/XChQokKTzH/YPRJH0QgmtSDphsVioVKkSP/30E2FhYaxcuZLly5cnK6GNjo5+5HEPDw+r43v++ec5evRo7MNMpUuXxsvLC1dXV+7fv5+kB6QeJyb+4sWLU6tWrUeea6sFH1xcXBg9ejTTpk3j8uXLzJo1i4CAgATnPcl397Dk/HFJe2Ijh7aS3GvH9E2zZs3IkyfPI8+NqZgB5khno0aNWLlyJb/99hvbt29n2bJlLFu2jGHDhhEYGEjZsmVjY5ozZw5DhgxhzZo1bN++ne3btzNlyhSmTJlC69atWbZsGc7OzsmK/b91aB/nSfo6NTz4ex5TPeFRbPG7KWJvSmhF0qEmTZqwcuVKrl27Fq/d1dWViIgIbt++negoYsyIqK0dP36cQ4cOkTt3bpYtW4aLS/z/a3nYn+pjRgGPHz+epPvEjDaVLFkytiRWanBycqJw4cJcu3btkSWrkuupp57izJkznDt3jtKlSyc4/t+SVWlZgQIFOH78OL169Ur29Atvb+94I7t//fUX/fr1Y8WKFfTt2zd2akaM0qVLU7p0ad555x0Mw2Djxo107tyZVatWMWvWLHr06GGzz5UcOXPmxM3NjfDwcM6cOUO5cuUSnBMzkp2SI6K5cuXCw8ODu3fvMn78eHLlypVi9xJJKzSHViSNMQzjsefEFGHPnz9/vPaY/0gmlnQdOnSIv/76ywYRJnTjxg0AfH19EySzwENrXsbM//zpp59iSyQ9SsOGDcmUKRObN2/m6tWrTxBx8kRHR8cml0mZ65tUdevWBXhord2HtT+JTJkyAYnX1H0SzZs3B2DhwoVPfK0CBQowcuRIwJwz/igWi4WGDRvGlnZ73PkpycXFJbbk1sP+wRVTJqt+/fopFoezszONGzcGbNMfIumBElqRNOabb76hW7duiT7UZBgGS5cuZdKkSQC8+OKL8Y43atQIgJEjR8ab43ju3Dm6deuWpGTZGiVKlMDZ2ZnDhw8nKPy+atUqJkyYkOj72rRpQ8WKFbl06RIdO3ZM8ADTvXv3Yh9uAciTJw/9+vUjNDSU1q1bc/jw4QTXDA8PZ+XKlUke9X2cyMhIPvroo9jR8DZt2tjkugB9+/bFycmJ+fPns2LFinjHli5dypIlS2x2rxgx/wh61Fxna7z++usUKlSIRYsW8d5773H79u0E51y5coVp06bF7h88eJAFCxZw9+7dBOeuWrUKiD89YdasWezfvz/Bubdv3479uXvwfHt46623AJgyZUqCleVmzpzJypUrcXV1ZcCAASkax/Dhw8mUKRPvvPMOP/74Y6LTjY4cOcLSpUtTNA6R1KIpByJpTEREBLNmzWLWrFn4+PhQsWJFcuXKxa1btzh27FjsSOErr7wS+8R0jCFDhrB48WJ+/vlnSpQogZ+fH//++y979+6lVq1a1KxZ06qn/x8nV65c9O3bl6+++oqGDRtSp04dfH19OXHiBAcOHOCjjz5i9OjRCd7n5OTEsmXLaNq0KWvXrqVgwYLUrl07dmGF33//nWzZssX70/vYsWO5fPky8+bNo0KFCpQvX56iRYvi4uLCxYsXCQoKIjQ0lLVr1yZ7Hu33338fLyG/fv06v//+e+zI9ocffkjNmjWt+o4SU7lyZUaPHs2QIUNo27Yt1atXj11YYc+ePbz11lt88cUXsaOqttChQwc2bdrEK6+8QpMmTciePTsA77zzDiVLlrT6ulmyZGHNmjW0atWKcePGMXXqVMqVK0f+/PkJCwvjzz//5I8//iB37ty89tprgDkF5sUXX8TDw4NKlSpRoEABIiMjOXz4MCdOnCBTpkyMGzcu9h5Lly6lW7du+Pr6UqFCBbJnz87NmzfZvn07wcHBlClTJvba9tK8efPYn/fGjRtTq1YtChYsyPHjxzlw4ADOzs58++23j1ykwxYqVarEnDlz6N69O927d+ejjz6idOnSsQtXHD58mIsXL9KpUyfat2+forGIpAq7Fg0TkQRCQkKM5cuXG/369TOqVq1q5M+f33B1dTU8PDyMYsWKGS+99JKxdu3ah77/2LFjRvv27Y3s2bMbbm5uRsmSJY3Ro0cb9+/ff2wd2sQKzMd4XB3a6OhoY/r06UblypWNrFmzGt7e3kbt2rWN+fPnG4bx6Bq5t2/fNj777DPDz8/P8PT0NNzc3IxChQoZbdq0iX3/f/38889G+/btjaeeespwdXU1smXLZjzzzDPGiy++aMybN88IDQ196Gf5r4fVoc2UKZNRqFAho1OnTg/9bpLy3T2uxujSpUuNWrVqGVmyZDE8PT2N2rVrG8uXLze2bt1qAEaNGjXinZ+UmsAPqzcbFRVljBkzxnj22WcNd3f3BLVYY+rQ1qtX76HXflRfhoSEGOPGjTNq1KhhZMuWzXB1dTXy5ctn+Pn5Ge+8846xY8eO2HMvX75sjB071mjRooVRpEgRI3PmzIaXl5dRunRpIyAgwDh+/Hi8a2/dutUYOHCgUbVqVSNv3rxGpkyZjLx58xo1atQwvv76a+POnTsPjflRn+NRffeg5NRiXrt2rdGiRQsjZ86chouLi5E3b16jY8eOxu7duxM939Z1aB+MedCgQUaZMmWMLFmyGO7u7kahQoUMf39/Y+zYscapU6ce+1lE0gOLYaTQ3yBFROSJjBo1iuHDh9OvXz8mTpxo73BERNIszaEVEbGjkydPcvPmzQTtK1euZMyYMVgsliSVXhIRcWSaQysiYkdz587l008/pWLFihQoUICIiAhOnDgRu6raiBEjqFy5sp2jFBFJ25TQiojYUbNmzTh58iS7du3ijz/+4N69e+TMmZPWrVvz5ptvJrq0rYiIxKc5tCIiIiKSrmkOrYiIiIika0poRURERCRdc8g5tNHR0Vy6dAlPT08sFou9wxERERGR/zAMg9u3b+Pr64uT06PHYB0yob106RIFChSwdxgiIiIi8hh//fVX7LLdD+OQCa2npydgfkFeXl4pfr+IiAjWr19PkyZNcHV1TfH7SdqgfndM6nfHpH53POrzlBcSEkKBAgVi87ZHcciENmaagZeXV6oltJkzZ8bLy0s/9A5E/e6Y1O+OSf3ueNTnqScp00P1UJiIiIiIpGtKaEVEREQkXVNCKyIiIiLpmhJaEREREUnXlNCKiIiISLqmhFZERERE0jUltCIiIiKSrimhFREREZF0TQmtiIiIiKRrSmhFREREJF1TQisiIiIi6ZoSWhERERFJ15TQioiIiEi6poRWRERERNI1JbQiIiIikq4poRURERGRdE0JrZ1ER8N338Hrr8OCBXD7tr0jEhEREUmfXOwdgCO6fBm6dYPAQHN/2jTInBkWLoSWLe0bm4iIiEh6oxHaVHb3LtStG5fMxggLg2HD7BOTiIiISHqmhDaVrV8Pp06Z2/nywfffQ6lS5v6BAxAUZLfQRERERNIlJbSpbN26uO3vvoNevaB//7i2H35I/ZhERERE0jMltKnIMGDtWnM7UyZo0MDcfuklcHc3t+fMgXv37BOfiIiISHqkhDYVHT8O58+b23XrQpYs5na2bNChg7l98yasXGmX8ERERETSJSW0qeiXX+K+7ubN4x/r2TNu+8cfUykgERERkQxACW0q+uUXS+z2fxNaf3/In9/cXr8ebtxIvbhERERE0jMltKnk8uXM/PabmdAWKhRX2SCGkxO88IK5HRkJy5alcoAiIiIi6ZQS2lRw8iR89FFt7t83E9rnngOLJeF5nTrFbS9YkErBiYiIiKRzSmhT2IkT0KiRC9evewDw7LMwdGji5/r5QZEi5vbGjfDvv6kUpIiIiEg6poQ2hQUGwuXL5nBsmTIGmzZBrlyJn2uxxE07iIqCRYtSKUgRERGRdEwJbQrr2xdGjYqiSJFbrF8fiY/Po89/cNrBhx+aI7wiIiIi8nBKaFPB++9HM3bsbw8dmX1QhQrmHFuAW7egdWtVPBARERF5lDSX0G7dupXWrVvj6+uLxWJh+fLl8Y7fuXOHvn37kj9/fjw8PChdujTffvutfYJNBje36CSdZ7HA7NlQtqy5f/IkBASkYGAiIiIi6VyaS2hDQ0MpX748kydPTvT44MGDWbduHXPmzOGPP/5g4MCB9O3bl5UZaHktT09YtQpy5jT358+HvXvtG5OIiIhIWpXmEtrmzZszevRo2rVrl+jxHTt20K1bN/z9/SlcuDCvv/465cuXZ8+ePakcacoqVAiGD4/bf/ddMAz7xSMiIiKSVrnYO4DkqlmzJitXrqRnz574+vqyefNm/vzzTyZMmPDQ94SHhxMeHh67HxISAkBERAQREREpHnPMPZJ7r5494auvXDh92sLmzbBqVSTNmyurTS+s7XdJ39Tvjkn97njU5ykvOd+txTDS7rifxWJh2bJltG3bNrYtPDyc119/nVmzZuHi4oKTkxPTpk2ja9euD73OiBEjGDlyZIL2efPmkTlz5pQI3Wa2bfNl/Hg/AMqW/ZePP95h54hEREREUl5YWBidO3cmODgYLy+vR56b7hLa8ePHM23aNMaPH0+hQoXYunUrH3zwAcuWLaNRo0aJXiexEdoCBQpw7dq1x35BthAREUFgYCCNGzfG1dU1We81DChd2hyldXIyOHcukrx5UyhQsakn6XdJv9Tvjkn97njU5ykvJCSEXLlyJSmhTVdTDu7evcuQIUNYtmwZLVu2BKBcuXIEBQUxfvz4hya0bm5uuLm5JWh3dXVN1R9Ca+/XqRN8+ilER1tYscKVvn1TIDhJMan9cyZpg/rdManfHY/6POUk53tNcw+FPUrMnFcnp/hhOzs7Ex2dtLJY6dGLL8Ztz59vvzhERERE0qI0N0J7584dTp06Fbt/9uxZgoKCyJEjBwULFqRevXq88847eHh4UKhQIbZs2cKsWbP43//+Z8eoU1aZMvDMM/DHH7B9O/z1FxQoYO+oRERERNKGNDdCu2/fPipWrEjFihUBs+5sxYoVGTZsGADz58/Hz8+Pl19+mdKlSzN27Fg++eQT3njjDXuGnaIslvijtIsW2S8WERERkbQmzY3Q+vv786jn1PLmzcuMGTNSMaK0oVOnuLq0CxbA4MH2jUdEREQkrUhzI7SSuJIloXx5c3vPHjhzxr7xiIiIiKQVSmjTkQenHSxcaL84RERERNISJbTpyAsvxG0vWGC/OERERETSEiW06UjRouBnLhpGUBCcOGHXcERERETSBCW06cyD0w40SisiIiKihDbd6dgxbnvNGvvFISIiIpJWKKFNZwoUMCseABw8CHfv2jceEREREXtTQpsO1axp/m9EBOzbZ99YREREROxNCW06VKtW3PaOHfaLQ0RERCQtUEKbDsWM0IISWhEREREltOlQyZKQPbu5vWMHPGKlYBEREZEMTwltOuTkFDdKe+0anDxp33hERERE7EkJbTqlaQciIiIiJiW06dSDD4Zt326/OERERETsTQltOuXnBy4u5vbWrfaNRURERMSelNCmU5kzQ/Xq5vaff8LFi/aNR0RERMRelNCmYw0axG1v2GC/OERERETsSQltOtawYdy2EloRERFxVEpo07Hq1c2pB2AmtKpHKyIiIo5ICW06likT1Kljbl+6ZM6lFREREXE0SmjTOU07EBEREUenhDadU0IrIiIijk4JbTpXvjxky2Zu//ab5tGKiIiI41FCm845O8etGvbvv5pHKyIiIo5HCW0GULt23PZvv9kvDhERERF7UEKbAcRUOgDYts1+cYiIiIjYg4s1b9q0aRMbNmxg+/btXLx4kWvXrpE5c2Z8fHwoW7Ys9erVo1WrVuTNm9fW8UoiqlQBNzcID9cIrYiIiDieJCe0oaGhTJw4kWnTpnH+/HmM/3/6yN3dnRw5cnD37l2OHDnCoUOHmDt3Lq6urrRu3ZpBgwZRK2aSp6QINzeoWtVMZs+cMWvS+vraOyoRERGR1JGkKQfffvstxYsX58MPP8TLy4uPP/6YDRs2EBwcTFhYGBcvXuT69etERERw/PhxfvzxR1544QXWr19P3bp1ad++PWfPnk3pz+LQNO1AREREHFWSEtp+/frRuHFjDh06RFBQEEOGDKF+/fp4enrGO89isVCiRAm6dOnC7Nmz+eeff/juu+/4/fffmT17dop8ADHpwTARERFxVEmacnD06FFKlCiR7It7eHjw6quv0qNHDy5cuJDs90vS1awJTk4QHa0FFkRERMSxJGmE1ppk9kHOzs4UKVLkia4hj+btDdWrm9t//AGnT9s3HhEREZHUorJdGUibNnHbq1bZLw4RERGR1JSshHb69Om8/fbbHD58OF57nz59bBqUWKd167htJbQiIiLiKJKc0I4ePZrevXszbdo0KleuzOTJk2OP7d+/P0WCk+R55hkoWtTc3roVgoPtG4+IiIhIakhyQjtr1iz27t3LjRs3GDp0KAMGDGDs2LEAsTVpxb4slrhR2shIWLfOvvGIiIiIpIYkL6yQNWtWKlasCMDQoUMpU6YML7/8Mt7e3lgslhQLUJKndWv46itze8UK6NTJvvGIiIiIpLQkJ7RRUVHx9tu1a8fSpUt5/vnncXGxagVdSQF160K2bHDrFixbBjdvQvbs9o5KREREJOUkecrBs88+y7r//A27WbNmLFq0iPDwcJsHJtZxdYWuXc3te/dg1iz7xiMiIiKS0pKc0M6bN49mzZolaG/evDl37961aVDyZHr3jtv+7jvQFGcRERHJyFSHNgMqXdqcegDmIgtaCldEREQyMqsS2i+++IJcuXJx6dKlRI9funQJHx8fJk6c+ETBifUeHKX95hv7xSEiIiKS0qxKaBctWkT58uXx9fVN9Livry8VKlRg/vz5TxScWK9DB/DxMbcXL4bz5+0bj4iIiEhKsSqhPXnyJM8+++wjz3n22Wc5efKkVUHJk3Nzg4AAczsqCiZMsG88IiIiIinFqoT27t27ZMmS5ZHnuLu7c+fOHauCEtsICAAPD3N72jS4ccO+8YiIiIikBKsS2oIFC7Jjx45HnrNz507y589vVVBiG7lyQc+e5nZYmObSioiISMZkVULbsmVLtm3bxg8//JDo8e+//55t27bROmYdVrGbwYPB6f97eexYOHPGvvGIiIiI2JpVCe37779Pnjx5eO2112jQoAFjxoxh1qxZjBkzhvr169O7d298fX354IMPkn3trVu30rp1a3x9fbFYLCxfvjze8X/++Yfu3bvj6+tL5syZadasmebqPkLRovD66+Z2aKg5Yhsdbd+YRERERGzJqjVrfXx82LRpE6+88gqbN29m8+bNWCwWjP+v4O/n58fcuXPxiXnMPhlCQ0MpX748PXv2pH379vGOGYZB27ZtcXV1ZcWKFXh5efG///2PRo0acezYscfO63VU48bB2rVmpYMtW2DyZOjXz95RiYiIiNiGVQktQMmSJdm7dy979+5lz549BAcHky1bNqpWrUqVKlWsDqh58+Y0b9480WMnT55k165dHDlyJLbKwpQpU8ibNy8//fQTr776qtX3zcg8PeGHH6BhQ3P/vfegeXMoXty+cYmIiIjYgtUJbQw/Pz/8/PxsEctjhYeHA2YFhRhOTk64ubmxbdu2hya04eHhse8FCAkJASAiIoKIiIgUjJjY+zz4v/ZQpw706ePElCnO3L0L3btHs2FDVOz8WrG9tNDvkvrU745J/e541OcpLznf7RMntMeOHeP48eOEhobSpUuXJ73cI5UqVYqCBQvywQcf8N1335ElSxYmTJjAxYsXuXz58kPfN2bMGEaOHJmgff369WTOnDklQ44nMDAw1e6VmLp1nVm6tD7//JOF7dudePPNo7Rpo6fEUpq9+13sQ/3umNTvjkd9nnLCwsKSfK7FiJn4mkx79+7ltdde4/Dhw7FtUVFRgPlgV7NmzZg/fz5t2rSx5vJmcBYLy5Yto23btrFt+/fvp1evXvz+++84OzvTqFEjnJycMAyDtWvXJnqdxEZoCxQowLVr1/Dy8rI6vqSKiIggMDCQxo0b4+rqmuL3e5TffrPQsKH57xhPT4PjxyOxYqqzJEFa6ndJPep3x6R+dzzq85QXEhJCrly5CA4Ofmy+ZtUI7dGjR2nQoAFOTk4MGjSI48ePx0sm69SpQ65cuVi0aNETJbSJqVy5MkFBQQQHB3P//n18fHyoVq3aI+fturm54ebmlqDd1dU1VX8IU/t+iWnQAF57zVxo4fZtC59+6sqkSXYNKcNLC/0uqU/97pjU745HfZ5ykvO9WjWDcvjw4YA5Wjp+/PgEc2gtFgs1atRg79691lw+Sby9vfHx8eHkyZPs27eP5557LsXuldGMGgUxBSG++w7+/NO+8YiIiIg8CasS2i1bttChQweKP+Ix+YIFCz5yXuvD3Llzh6CgIIKCggA4e/YsQUFBXLhwAYBFixaxefNmzpw5w4oVK2jcuDFt27alSZMm1nwUh5Q3L7z7rrkdGQlDhtg3HhEREZEnYVVCe/v2bXLnzv3Ic+7evRs7pzY59u3bR8WKFalYsSIAgwcPpmLFigwbNgyAy5cv06VLF0qVKkX//v3p0qULP/30U/I/hIN76y0zsQVYuhS0NoWIiIikV1bNoS1QoEC8h8ESc+DAAYoVK5bsa/v7+/Oo59T69+9P//79k31diS9LFhg0yKxJaxgwYQJ88429oxIRERFJPqtGaFu1asX69ev59ddfEz2+cOFCdu3aFa86gaQ9r78OWbOa2zNmwLVr9o1HRERExBpWJbRDhgzB19eXFi1a8Nprr7Fv3z4AvvnmG7p06ULnzp0pXLgwgwcPtmmwYlvZskHMWhT37sGUKXYNR0RERMQqViW0Pj4+bNmyBT8/P6ZPn86aNWswDIO+ffsyd+5c/Pz82LhxI97e3raOV2xswABiVwv76iv4/0XURERERNINq1cKK1q0KNu3bycoKIhdu3Zx48YNvLy8qFatWqothStPrnBh6NwZ5syB69fNpHboUHtHJSIiIpJ0T7z0bYUKFahQoYINQhF7GT4cfvoJoqJg/HgICIAcOewdlYiIiEjSWDXl4GHOnTvHtGnTmD17Nnfu3LHlpSUFFS8OPXqY2yEhZlIrIiIikl5YldB++umnFClShJs3b8a2bd68mTJlyvDGG2/QvXt3KlWqxI0bN2wWqKSsoUMhUyZze/JkuH3bvvGIiIiIJJVVCe3y5cspXLgw2bNnj2177733iI6OZuTIkfTp04dTp07x5Zdf2ipOSWEFC0K3buZ2SAj88IN94xERERFJKqsS2nPnzlG6dOnY/UuXLrF3714CAgL46KOPmDRpEg0bNmTZsmU2C1RS3sCBcdtffWXOqRURERFJ66xKaENCQsiWLVvs/tatW7FYLLRu3Tq2rVKlSly4cOGJA5TUU7o0NG1qbp89CytW2DceERERkaSwKqHNkycP58+fj90PDAzEzc2NatWqxbbdu3cPi8Xy5BFKqnpwlPaLL8xlcUVERETSMqsSWj8/P1asWMHq1av59ddfWbBgAfXr18fNzS32nLNnz+Lr62uzQCV1NG1qjtQC7NgBW7bYNx4RERGRx7F66dvIyEiee+45mjZtyr179xgyZEjs8fDwcLZu3RpvxFbSB4sFPvggbn/0aPvFIiIiIpIUViW0lSpVYteuXQwaNIhBgwaxY8cOatWqFXv84MGD1K9fn86dO9ssUEk9L74IxYqZ2xs2wM6d9o1HRERE5FGsXimsfPnylC9fPtFj1atXV4WDdMzFBYYMgV69zP0hQ8zE1smmy3CIiIiI2IZVKUpUVBQhISFER0cnejw6OpqQkBCiVPcp3erSBQoXNrc3b4YxY+wZjYiIiMjDWZXQjhw5kty5c3P9+vVEj1+/fp08efLwySefPFFwYj+urvD99+acWoBhw2DTJvvGJCIiIpIYqxLa1atX07BhQ3x8fBI97uPjQ6NGjVihQqbpWsOGMGKEuR0dDa++av6viIiISFpiVUJ75swZSpUq9chzSpYsydmzZ60KStKOjz6CmOf9zpwxpx+IiIiIpCVWJbQRERE4PeYJIYvFwr1796wKStIOJyfo1y9uf8YM+8UiIiIikhirEtrixYuzcePGR56zceNGihQpYlVQkrY89xzErHS8ZAkEB9s1HBEREZF4rEpo27dvT1BQEMOGDUtQySAqKoqhQ4cSFBREx44dbRKk2Je7O7z8srl99y4sWGDfeEREREQeZFUd2rfeeov58+fzySefMH/+fOrXr89TTz3F33//zaZNmzh9+jTPPPMMb7/9tq3jFTvp0QMmTza3p0wxa9Q6O9s3JhERERGwMqHNmjUrW7dupU+fPixbtoxTp07FHnNycuL555/nm2++IWvWrDYLVOyrUiWoUAGCgszX+PHw3nt2DkpERESEJ1gpzMfHh8WLF/PPP/+wb98+goODyZYtG1WqVCF37ty2jFHSAIsFvv4a6tYFwzDr0rZoAWXL2jsyERERcXRWJbQNGjSgVq1afPzxx+TJk4eWLVvaOi5Jg2rXhnfegXHj4P596N4d9u7VkrgiIiJiX1alIrt379aytg5q1Ch49llz+8ABmDvXvvGIiIiIWJXQlipVivPnz9s6FkkH3Nzgq6/i9j/8EFRuWEREROzJqoS2X79+rFixgmPHjtk6HkkHGjaE5s3N7b/+gokT7RuPiIiIODar5tAWLVoUf39/qlevTu/evfHz8yNPnjxYLJYE59atW/eJg5S0Z9w4+OUXiI6G0aOhc2fIn9/eUYmIiIgjsiqh9ff3x2KxYBgGX3zxRaKJbAzNtc2YypSBV1+FqVPh9m14801YscKshiAiIiKSmqxKaIcNG/bIJFYcw9ixsHIlXLkCq1bB4sWgxeFEREQktVmV0I4YMcLGYUh6lD27WZs2Jol9802oUUNTD0RERCR1qYKoPJEOHaBtW3P72jV44QWIiLBrSCIiIuJglNDKE7FYYPp0KFTI3N+5E95/374xiYiIiGOxOqH966+/6N27N8WKFcPDwwNnZ+cELxcXq1fWlXQkRw5YtAhcXc39L7+EU6fsGpKIiIg4EKsS2jNnzlCpUiWmT59O1qxZCQ8Pp2DBgpQoUQIXFxcMw6BcuXLUqVPH1vFKGuXnB0OHmtvR0TBmjH3jEREREcdhVUI7cuRIgoOD2bBhA7///jsAPXr04I8//uDcuXO0adOG0NBQFi9ebNNgJW3r3x+yZTO3Z82Cc+fsGY2IiIg4CqsS2l9//ZUWLVpQr1692DbDMADIly8fCxYsAGDIkCE2CFHSC29vM6kFiIzUKK2IiIikDqsS2mvXrlGqVKnYfRcXF8LCwmL33dzcaNy4MatXr37yCCVdGTAAPD3N7WnTzDq1IiIiIinJqoQ2V65chIaGxts/95+/L7u4uHDr1q0niU3SoRw54ubSGga8/DIcPWrfmERERCRjsyqhffrppzl9+nTsftWqVfnll184c+YMAP/++y+LFy+mWLFitolS0pW334ZOncztO3fMWrWqTSsiIiIpxaqEtnnz5mzatCl2BHbgwIHcvn2bcuXK4efnR4kSJbhy5Qr9+vWzZaySTlgs8MMPUKGCuX/iBMyebdeQREREJAOzKqHt06cPmzdvxtnZGQB/f3/mz59PoUKFOHLkCHny5GHixIm89tprNg1W0o/MmWHSpLj9jz+G+/ftF4+IiIhkXFatfODl5UW1atXitXXs2JGOHTvaJCjJGGrVgqZN4ZdfzBJeM2fC66/bOyoRERHJaJI8QvvCCy/w008/ERwcnJLxSAYzcmTc9gcfwObNdgtFREREMqgkJ7SBgYG8/PLL5M6dm0aNGjFp0iQuXLiQkrFJBlCtGjz3nLl94wY0agTff2/fmERERCRjSXJCe+3aNQIDA+nduzenT5+mf//+FClShEqVKjFq1CiCgoJsEtCYMWPw8/PD09OT3Llz07ZtW06cOBHvnHv37hEQEEDOnDnJmjUrHTp04J9//rHJ/cX2Zswwpx4AREXBG2/A8eP2jUlEREQyjiQntM7OzjRs2JCJEydy9uxZDhw4wND/Lzg6YsQIKleuTJEiRRg4cCAbN24kKirKqoC2bNlCQEAAu3btIjAwkIiICJo0aRKv7u2gQYNYtWoVixYtYsuWLVy6dIn27dtbdT9Jedmzw5o1ZiILZlL77rv2jUlEREQyDquqHABUqFCBESNGcODAAc6fP8+ECRMoVqwY33zzDY0bNyZ37tx07dqVJUuWcOfOnSRfd926dXTv3p1nn32W8uXLM3PmTC5cuMD+/fsBCA4OZvr06fzvf/+jQYMGVK5cmRkzZrBjxw527dpl7ceRFObsDF98AU89Ze6vWgUbNtg3JhEREckYrKpy8F8FChSgf//+9O/fn1u3brFmzRqWL1/O8uXLmTNnDu7u7vGWxk2OmIfQcuTIAcD+/fuJiIigUaNGseeUKlWKggULsnPnTqpXr57gGuHh4YSHh8fuh4SEABAREUFEKlT8j7lHatwrLXN1hVGjLPTqZf7YDRpksHNnJJky2TmwFKJ+d0zqd8ekfnc86vOUl5zv1mIYhpFSgdy/f58NGzawcuVKpkyZkuz3R0dH06ZNG27dusW2bdsAmDdvHj169IiXoIK5Wln9+vX57LPPElxnxIgRjHzwcfv/N2/ePDJnzpzsuMR60dHwzjv1OH06GwAdO57g5Zc1oVZERETiCwsLo3PnzgQHB+Pl5fXIc20yQvswmTJlonnz5jRv3tyq9wcEBHDkyJHYZNZaH3zwAYMHD47dDwkJoUCBAjRp0uSxX5AtREREEBgYSOPGjXF1dU3x+6V1vr5Qu7ZBZKSFJUtK0L9/MapVS7F/V9mN+t0xqd8dk/rd8ajPU17MX9STwqqEdtSoUY89x8nJCS8vL0qWLIm/vz9ubm7Jukffvn1ZvXo1W7duJX/+/LHtefPm5f79+9y6dYts2bLFtv/zzz/kzZs30Wu5ubklen9XV9dU/SFM7fulVdWqwbBh5is62kLHji58+21cea+MRv3umNTvjkn97njU5yknOd+rVQntiBEjsFgssfsPzlr4b7vFYiFbtmxMmDCBrl27PvbahmHQr18/li1bxubNmylSpEi845UrV8bV1ZUNGzbQoUMHAE6cOMGFCxeoUaOGNR9H7OCDD+Dnn2HXLrhyBdq2hYCA+MvlioiIiCSFVVUONm3aRKtWrXBzc+O1117jxx9/ZN26dfz444+89tpruLm50bp1axYvXswHH3xAZGQkPXv25Ndff33stQMCApgzZw7z5s3D09OTK1eucOXKFe7evQuAt7c3vXr1YvDgwWzatIn9+/fTo0cPatSokegDYZI2ubjA8uXw4GyUyZNhzx67hSQiIiLplFUJ7cmTJ9myZQsHDhzg22+/pUuXLjRp0oQuXbrw7bffcuDAATZt2sStW7cYPXo0O3fuxMXFhS+++OKx154yZQrBwcH4+/uTL1++2NeCBQtiz5kwYQKtWrWiQ4cO1K1bl7x587J06VJrPorYUZ48Zn3aTz6Ja/vqK/vFIyIiIumTVQntV199RadOnShVqlSix0uVKkWnTp2YMGECAKVLl6Z169bsScLwm2EYib66d+8ee467uzuTJ0/mxo0bhIaGsnTp0ofOn5W0zWKBQYMgVy5zf+FC+Ptv+8YkIiIi6YtVCe2pU6di68I+TM6cOTl9+nTsfrFixZK1wII4Dg+PuFXEIiPNqQciIiIiSWVVQuvj48PatWt5WAlbwzBYu3YtOXPmjG27efMm3t7e1kUpGV6fPubCCwBjxkDt2rB2rX1jEhERkfTBqoT2xRdf5NChQ7Rp04ZDhw7FO3bo0CGee+45Dh8+zEsvvRTbvmfPHp555pkni1YyLF9f6NUrbn/7dmjVCtats19MIiIikj5YVbZr5MiR7Nu3jzVr1vDzzz+TJUsWfHx8+PfffwkNDcUwDOrWrRu7OteVK1coXLgwL7zwgk2Dl4zl66+haFGYOROOHTNXFevUCXbvhodM1xYRERGxboTWw8ODX3/9lWnTplG3bl1cXV25cOECrq6u1KtXj2nTprFx40Y8PDwAczGEZcuWxRuxFfkvFxd45x04fBjatTPbQkLM0l7HtTquiIiIPITVS986OTnRq1cvej34d2IRG3ByglmzzHm0v/8O585B9eqweDE0amTv6ERERCStsWqEViSlZc1qriRWoYK5HxwMzZrBd9/ZNSwRERFJg6weoY0RFRXFtWvXCA8PT/R4wYIFn/QW4qB8feG33+Dll2HlSoiKMst7nT8Pn35q7+hEREQkrbB6hHb//v00bdqUrFmz4uvrS5EiRRK8ihYtastYxQFlzQpLl8Lbb8e1jRkDBw7YLyYRERFJW6xKaIOCgqhTpw47d+6kSZMmGIZBuXLlaNKkCbly5cIwDOrVq0eXLl1sHa84IGdn+Pxz+PjjuDZNPRAREZEYViW0H/9/ZrF7925WrFgBQLt27Vi7di3nzp3jjTfe4MiRIwwfPtx2kYrDGzDAHLEFmDcPbt+2bzwiIiKSNliV0G7bto02bdrEWyghZtUwDw8PJk2ahK+vL0OGDLFNlCKAp6c5nxbgzh2YO9e+8YiIiEjaYFVCGxwcHG9+rKurK3fu3Im7qJMT/v7+bNiw4ckjFHlA795x299+Cw9ZfVlEREQciFUJbe7cubl582bsft68eTl58mS8c+7du0dYWNiTRSfyHxUrQtWq5vbvv5tVEERERMSxWZXQli5dmhMnTsTu16pVi/Xr17Nz504A/vjjDxYuXEgprVcqKaBfv7jtcePsF4eIiIikDVYltC1btmTr1q1cvnwZgPfeew/DMKhduzY+Pj6ULVuWW7duaQ6tpIhOnSCmvPGaNXDkiH3jEREREfuyKqF94403+Pvvv8mZMycA5cuXZ8OGDTRr1oxcuXLRqFEjVq1aRbt27WwarAiAqysMHhy3//nn9otFRERE7M+qlcJcXV3JkydPvLaaNWuyZs0amwQl8jivvgqjRsGNG2YJrzFjzJXFRERExPFYvVKYiD1lyQJ9+pjbkZEwe7Z94xERERH7sWqENsa///7LsWPHuHTpEhEREYme07Vr1ye5hchD9ewJn3xibv/wA7z7Llgs9o1JREREUp9VCe3du3fp168fs2fPJjIyMtFzDMPAYrEooZUUU7Qo+PvD5s3w55+wcyfUrGnvqERERCS1WZXQ9u/fnx9++IFy5crx/PPPky9fPlxcnmiwV8QqPXqYCS2Yo7RKaEVERByPVVnokiVLqFKlCjt37sTZ2dnWMYkkWYcO0Lcv3L4NCxbAhAnmErkiIiLiOKx6KCwqKgp/f38ls2J3WbLAiy+a23fuwKRJ9o1HREREUp9VCa2fn1+CpW5F7OXtt8Hp/3+SP/8cQkLsG4+IiIikLqsS2o8//pj169ezevVqW8cjkmwlSsArr5jbN2/CxIn2jUdERERSl1VzaGvUqMH69etp06YNlSpVonz58nh5eSU4z2KxMHTo0CcOUuRxhg6FuXMhKgq++ALeeANy5bJ3VCIiIpIarEpor1+/zgcffMDNmzfZsGEDGzZsSPQ8JbSSWooXh65dYcYMuHULBg6EOXPsHZWIiIikBqsS2n79+rFt2zZatGjBiy++qLJdkiZ88gksW2YmtHPnwksvQcuW9o5KREREUppVWei6devw9/fXHFpJU/LlM8t29ehh7vfuDSdOmJUQREREJOOy6qEwwzCoUqWKrWMReWLdukHTpub233+bUxBEREQkY7Mqoa1Vqxa///67rWMReWIWC4wbF7c/YYL5oJiIiIhkXFYltOPHj2fv3r1MUhV7SYPKlYPGjc3tM2dgxQr7xiMiIiIpy6o5tOPGjaNcuXIMGDCAiRMnUq5cuYeW7Zo+ffoTBymSXG+9BYGB5vb48dC+vX3jERERkZRjVUI7c+bM2O1Tp05x6tSpRM9TQiv20qQJlCkDR47Azp0weTIEBNg7KhEREUkJViW0Z8+etXUcIjZlscCIEfD88+Z+//5QrBg0a2bXsERERCQFWJXQFipUyNZxiNhchw7w3nvw2WcQHQ2dOsHx42Z5LxEREck4rHooTCS9+PTTuPmzISHw5Zd2DUdERERSQJIS2vnz5z/RTf766y+2b9/+RNcQsYaTE0yaBJkymfvffgvBwfaNSURERGwrSQntK6+8Qrly5fjxxx+5c+dOki++b98+evfuzdNPP82mTZusDlLkSeTLB127mtshIfDdd/aNR0RERGwrSXNo9+7dy+DBg+nRowcBAQE0a9aMatWqUblyZfLkyUO2bNm4d+8eN27c4MSJE+zevZvAwEBOnjyJl5cXo0aNYsCAASn9WUQe6p13YPp0MAxzsYX+/cHd3d5RiYiIiC0kKaGtWLEimzZtYuPGjUyZMoWVK1eydOlSLBZLoucbhkHRokUZO3Ysr776KtmzZ7dp0CLJVaIEtGsHS5fClSswdqxZBUFERETSv2RVOWjQoAENGjTg1q1bbN26le3bt3Px4kWuX7+Oh4cHPj4+lC1blnr16lGuXLmUilnEKiNHwsqVEBkJY8bAiy9CqVL2jkpERESelFVlu7Jly0abNm1o06aNreMRSTFlysDbb5ujs/fvQ+/esHmzWbNWRERE0i+V7RKHMnQoFC1qbm/dClOm2DceEREReXJKaMWhZM4cv8rB22/DyZP2i0dERESenBJacTiNGsGbb5rbd+9Cly4QEWHfmERERMR6aS6hHTNmDH5+fnh6epI7d27atm3LiRMn4p0zdepU/P398fLywmKxcOvWLfsEK+nWuHFQvLi5vXu3WdZLRERE0qc0l9Bu2bKFgIAAdu3aRWBgIBERETRp0oTQ0NDYc8LCwmjWrBlDhgyxY6SSnmXJAnPmgKuruf/VVzB7tn1jEhEREetYVeUgJa1bty7e/syZM8mdOzf79++nbt26AAwcOBCAzZs3p3J0kpFUq2Yui9u7t7nfsyccPQoffQRZs9o3NhEREUm6NJfQ/ldwcDAAOXLksPoa4eHhhIeHx+6HhIQAEBERQUQqTJ6MuUdq3EuSp0cP2LfPiWnTnImMhM8+g0WLDDZvjiRv3ie7tvrdManfHZP63fGoz1Necr5bi2EYhq1ufO7cOQIDA3F3d6ddu3ZkfcJhrujoaNq0acOtW7fYtm1bguObN2+mfv363Lx5k2zZsj30OiNGjGDkyJEJ2ufNm0fmzJmfKEZJ/6KiLCxYUJKlS4sTGekMQKVK/zB06C7VqBUREbGTsLAwOnfuTHBwMF5eXo8816qE9tNPP2XatGkcOHAgdlnbzZs306pVK+7evQtAsWLF2LVr1xONrPbp04e1a9eybds28ufPn+B4UhPaxEZoCxQowLVr1x77BdlCREQEgYGBNG7cGNeYSZuS5pw8CY0auXD5spnFTpgQRUBAtNXXU787JvW7Y1K/Ox71ecoLCQkhV65cSUporZpysHz5cgoXLhybzAK89957REdHM3LkSK5cucI333zDl19+yahRo6y5BX379mX16tVs3bo10WQ2Odzc3HBzc0vQ7urqmqo/hKl9P0me0qVhxgxo1szcf/99Z55+2pmWLZ/suup3x6R+d0zqd8ejPk85yflerapycO7cOUqXLh27f+nSJfbu3UtAQAAfffQRkyZNomHDhixbtizZ1zYMg759+7Js2TI2btxIkSJFrAlRxCpNm0L//uZ2eDi0aRN/IQYRERFJe6xKaENCQuL9iX/r1q1YLBZat24d21apUiUuXLiQ7GsHBAQwZ84c5s2bh6enJ1euXOHKlSuxUxkArly5QlBQEKdOnQLg8OHDBAUFcePGDWs+jkg848bBCy+Y29HR8MYb8Oqr5iIMIiIikvZYldDmyZOH8+fPx+4HBgbi5uZGtWrVYtvu3buHxYonaqZMmUJwcDD+/v7ky5cv9rVgwYLYc7799lsqVqzIa6+9BkDdunWpWLEiK1eutObjiMTj5gY//WQuixtj+nTw8wNVihMREUl7rJpD6+fnx4oVK1i9ejXu7u4sWLCA+vXrx5unevbsWXx9fZN97aQ8ozZixAhGjBiR7GuLJJWTE3z+OZQtC336QFiYWaO2fn1zjm27dmYd21y54N49uHgRvLygXDlwdrZ39CIiIo7FqoR2yJAhrFmzhueeew4Ai8USb9Wu8PBwtm7dSrt27WwTpYiddO0KVapAly5w4IDZtm6d+UpM9uzQpAm8/DI0bJh6cYqIiDgyq6YcVKpUiV27djFo0CAGDRrEjh07qFWrVuzxgwcPUr9+fTp37myzQEXspXRp2LMHfvgBnnrq0efevAkLFpgPkxUv7sKWLfmxXaVnERERSYzVK4WVL1+e8uXLJ3qsevXqVlU4EEmrnJ3NVcW6djVHan/5Bc6dg+vXwdUV8ueHCxdg40YzqQW4fNnChAmV2b8/mtGjoU4du34EERGRDMumS9+GhISwe/du3N3dqV27tlUPhYmkZc7O5sNhfn6JH4+IgPXrYcoUWLPGbNu61Ym6dc0pCD/9BD4+qReviIiII7BqysG0adOoV68eN2OGooDff/+dUqVK0axZM/z9/alTpw5hYWE2C1QkPXB1hZYtYfVqWL48krx578Qe27AB2rY1HyITERER27EqoZ09ezbh4eHxVgp76623uHr1Kj169KBFixbs3LmTKVOm2CxQkfSmRQuDSZM28v33keTJY7bt2GHWtNW8WhEREduxKqH9888/482fvX79Ops2beLVV1/l+++/Z9WqVfj5+TF37lybBSqSHrm4GHTtarBmDWTObLbNnQsffmjfuERERDISqxLaW7du4fPARMDffvsNgPbt28e21a5dm3Pnzj1ZdCIZROXKMGcOxEwrHzMG/vc/+8YkIiKSUViV0ObMmZPLly/H7m/YsAFnZ+d4pbsMwyAiIuLJIxTJINq1g0mT4vbfegtmzbJfPCIiIhmFVQltuXLlWLFiBUeOHOHUqVPMmzePWrVqkSVLlthzzp07R758+WwWqEhG8Oab8OAidz17wqpVdgtHREQkQ7AqoX333Xe5efMm5cuXp2TJkty6dYvBgwfHHo+Ojmbbtm1UrlzZZoGKZBTDhkHfvuZ2VBS88ALMnKkHxURERKxlVUJbv359Vq5cSbt27WjXrh0LFiygdevWsce3b9+Or69vvDm1ImKyWOCrr+Cll8z9e/fMRRvatoW1ayE83K7hiYiIpDtWL6zQsmVLWrZsmeixOnXqcPDgQauDEsnonJzgxx/Bw8NcUhdg5Urz5e0NnTubc24zZzaX2y1c2K7hioiIpGlWjdD+140bN/jrr79scSkRh+HqCtOnw5IlkCtXXHtwsLnSWJMmULs2FCkCL74IZ8/aL1YREZG0zOqENjg4mAEDBpAnTx58fHwoUqRI7LHdu3fTokUL9u/fb5MgRTKy9u3h/HlYuhReeSWuXu2DFiyA4sWhZk349FO4dSvVwxQREUmzrEpob9y4QbVq1fj6668pUKAAzzzzDMYDT7SUK1eO7du3a2EFkSTKnNmcYjB7Nly+bI7cvvMO9O8PMSWfo6Nh505zUYann4apU802ERERR2dVQjtixAj+/PNP5s+fz759++jYsWO84x4eHtSrV4+NGzfaJEgRR+LlZZbzGjfOfHjs5EmzMsIzz8Sdc+0a9O4N3btDZKTdQhUREUkTrEpoV65cSatWrXjhhRceek7hwoW5ePGi1YGJiMnbG0aOhGPH4NQp6NQp7tjs2eb8WlVGEBERR2ZVQnv58mVKly79yHPc3NwIDQ21KigRSVyxYjB/vvkgWaZMZtuSJebc2pMn7RubiIiIvVi99O3jqhocP35cK4WJpJD27c0SXx4e5v6BA1CpEqxYYd+4RERE7MGqhLZu3bqsWLHioVMKjh07xrp162jUqNETBSciD9e0KezYASVKmPt37pgPlk2YoIfFRETEsViV0H744YdERUVRq1Yt5s6dy7Vr1wD4448/mD59Og0aNMDNzY133nnHpsGKSHwVKsD+/eY8WjCXzx08GAoUMCskXLpk1/BERERShVUJbdmyZVmwYAG3bt2ia9eufPPNNxiGQZkyZXjttde4e/cuCxcu5Omnn7Z1vCLyH1mzwrx5MHRoXNulS/D11/Dss+aDYw9U1RMREclwrF76tk2bNpw9e5Yff/yR3bt3c+PGDby8vKhWrRo9evQg14NLH4lIirJYYNQoqFYNvvsO1q83Kx/cugVdu8L48WZd2+efB3d3e0crIiJiW1YntAA5cuRg0KBBtopFRJ5Qy5bm69o1GDgQYtY2OXQIunSBfv3gpZfM0Vw9sykiIhmF1UvfikjalSsXzJkDa9aYo7Yxbt2CKVPMube//mqv6ERERGzriUZo9+zZw969e7l16xZRUVEJjlssFoY+OLFPRFJVixbQvDn89ht8/71ZszYsDK5ehcaNoWxZ8Pc3HyqrUcOcuiAiIpLeWJXQ3rhxg7Zt27J9+3aMRzxtooRWxP4sFqhb13x98YU5p3bdOvPY4cPm6+uvoXhxaNLETHBbtoTMme0atoiISJJZldAOHjyYbdu24e/vT7du3cifPz8uLk802CsiqcDHx5yGMGkSzJoFBw/G1aw9dcp8ffMNZMsG3bubD5L5+tozYhERkcezKgtdvXo1VatWZcOGDVj0N0qRdMXJyaxR278/BAebK47NmAFbtsQlt7duwZdfwtSp8P77EBAAOXLYM2oREZGHs+qhsLt371K3bl0lsyLpnLe3Wf1g40a4eRPWroVu3eJKe4WFwbBhkDcvPPccLFhgtomIiKQlViW0FSpU4Ny5czYORUTsycsLmjWDmTPhr7+gb19wdjaPRUSYI7kvvgh58phJ8M8/w+3bdg1ZREQEsDKhHT58OCtXrmTXrl22jkdE0oBcucwHxY4cMZfSfXAe7Z07Zkmwli3NEd4KFSAw0G6hioiIWDeH9sqVK7Rs2ZJ69erx8ssvU6lSJby8vBI9t2vXrk8UoIjYT6lSZmWEceNg61Zzid3Fi805tmAuqfv779CqFaxaZVZJEBERSW1WJbTdu3fHYrFgGAYzZ85k5syZCebTGoaBxWJRQiuSATg7Q/365mvSJHOu7bp1Zn3bY8fg/n1o2xZGjjTLfp0/DxcvQokSUKUK5M5t708gIiIZmVUJ7Q8//KAHwkQclJubmby2bQuRkdCpEyxdCnfvwrvvJv6eChWgTRvzVamSWRvXMLSQg4iI2IbVI7QiIi4u8NNP8PLL5lSEhwkKMl+jRkHOnOZDZmFhUKsWvPQSvPIKZMmSWlGLiEhGY1VCO2vWLCpUqEC5cuUees6RI0c4cOCAphyIZHCZMsHChXDihDmP9sgRKFoUChaEo0dh82bYvz/u/OvX47a3bDFf48eby/I+4v9SREREHsrqEdoRI0Y8MqFdsWIFw4YNU0Ir4gAsFvMBslKlEj9+6RKsXm2W/tq3z6yOcP8+xFT/O3UKqleHuXOhXbtUC1tERDIIq8p2JUVUVBROTil2eRFJR3x94fXXzaT2yhVzNPfMGdi715xTC+Yc3G7dzOMiIiLJkWIZ58GDB8mhtTJF5CEsFrMCwvbt0L692Xb7Nnz4oX3jEhGR9CfJUw4aNGgQb3/mzJls3rw5wXlRUVFcvHiRc+fO8cILLzxxgCKSsbm7w7ffwoYNEBwMM2ZAnz5msisiIpIUSU5oH0xeLRYL586dS3T5WycnJ3LkyEHHjh358ssvbRCiiGR0Pj4wfLi5KplhmFMP1q+Hp56yd2QiIpIeJHnKQXR0dOzLMAxGjBgRry3mFRkZydWrV5k/fz558uRJydhFJAMJCIh7qOzYMahRw6ySICIi8jhWzaHdtGkT3bp1s3UsIuLAMmUyqyAUKWLu//UX1K5tLrkrIiLyKFYltPXq1aNQoUK2jkVEHNzTT8OOHXGVD27dgsaNYdEiu4YlIiJpnFV1aAHu37/P8uXL2bt3L7du3SIqKirBORaLhenTpyfrumPGjGHp0qUcP34cDw8PatasyWeffUbJkiUBuHHjBsOHD2f9+vVcuHABHx8f2rZty8cff4y3t7e1H0dE0oi8ec3FGDp2hF9+MevVdupk1rIdMMDe0YmISFpkVUJ7/vx5GjduzOnTpzEM46HnWZPQbtmyhYCAAPz8/IiMjGTIkCE0adKEY8eOkSVLFi5dusSlS5cYP348pUuX5vz587zxxhtcunSJxY9ae1NE0g1PT3PVsddegx9/NB8UGzgQTp40VxVzd7d3hCIikpZYldAOGjSIU6dO0aVLF3r27En+/PlxcbF6sDeedevWxdufOXMmuXPnZv/+/dStW5cyZcqwZMmS2OPFihXjk08+4ZVXXiEyMtJmcYiIfbm6miW8ChSA0aPNtsmTzaVy582DsmXtG5+IiKQdVmV/GzdupGHDhvz444+2jieB4OBggEcu0hAcHIyXl9dDk9nw8HDCw8Nj90NCQgCIiIggIiLChtEmLuYeqXEvSTvU77YxbBjkz29hwABnwsMtHDkC1asb/PBDFO3bP/wvRPaifndM6nfHoz5Pecn5bi3Go+YMPISXlxdvvPEG48aNS+5bkyU6Opo2bdpw69Yttm3blug5165do3Llyrzyyit88skniZ4zYsQIRo4cmaB93rx5ZM6c2aYxi0jKOH/ekwkTKnPuXNxc+erVL1G+/L/UrHkJb+/7doxORERsLSwsjM6dO8cOXD6KVQlt48aNcXd3Z9WqVVYHmRR9+vRh7dq1bNu2jfz58yc4HhISQuPGjcmRIwcrV67E1dU10eskNkJboEABrl279tgvyBYiIiIIDAykcePGD41RMh71u+3duwd9+jgzd278Ai3Zsxt88UUUdesaXLpk4dIluHkTGjQwYsuApRb1u2NSvzse9XnKCwkJIVeuXElKaK2acjB27Fjq1q3L4sWLef75560K8nH69u3L6tWr2bp1a6LJ7O3bt2nWrBmenp4sW7bskT9Mbm5uuLm5JWh3dXVN1R/C1L6fpA3qd9txdYXZs6FyZRg50lwqF+DmTQs9eyb8v7OcOWHvXlI9qQX1u6NSvzse9XnKSc73mqSEdtSoUQna6tevT6dOnahXrx6VKlVKNHO2WCwMHTo0ycEAGIZBv379WLZsGZs3b6ZIIv8lCgkJoWnTpri5ubFy5Urc9ciziMOwWGDQIOjXD4KC4H//g59+Svzc69ehfXvYvh00u0hEJONKUkI7YsSIhx7bvHkzmzdvTvSYNQltQEAA8+bNY8WKFXh6enLlyhUAvL298fDwICQkhCZNmhAWFsacOXMICQmJfcjLx8cHZ2fnZN1PRNInFxeoUsWsePDii2ZFBIsFnnoKfH3hhx/g1Ckz6e3UCaZNM2vciohIxpOkhHbTpk0pHUesKVOmAODv7x+vfcaMGXTv3p0DBw6we/duAIoXLx7vnLNnz1K4cOHUCFNE0pA2bczXf9uqVYPQUFi9GkqUgC++MGvbiohIxpKkhLZevXopHUesxz2j5u/v/9hzRESefRYWLoSXXzaX0L19G15/HcLDoW9fe0cnIiK2pFUIRCTDatEC/vwT3n/fnIIA5tzby5ehWTMzyT11CnLlglKl4OpVOH0aatWCChXsGrqIiCSDVQnt1q1bH3uOk5MTXl5eFC9eXLVeRcRufHxg+nTIlw9iSlV/+qn5ehg3N9i8GapXT5UQRUTkCVmV0Pr7+2OxWJJ0rpOTE40bN+bzzz/n2WefteZ2IiJP7OOPISICkrIeTHg4tG0Le/ZAwYIpHpqIiDwhqxLaYcOGsWfPHtatW0fJkiWpWbMmefLk4Z9//mHnzp0cP36c5s2bU6xYMQ4cOMC6devYuXMnu3fvpkSJErb+DCIij2WxwGefQa9esGUL7N8POXLA00+bUw1OnDCnHuzYYZb5+ucfaNrUnIdbtqy9oxcRkUexKqFt2LAhY8eO5YcffqB79+4Jjv/444/06dOHDz74gIkTJzJnzhy6du3K6NGjmTVr1pPGLCJitRIlzNfDqh1cv25ONTh1Co4fN0uDjRgBAwaolq2ISFrl9PhTEho6dCitW7dONJkF6NatGy1btuSjjz4C4JVXXsHf35+NGzdaHaiISGrImRPWroUyZcz9+/dhyBAoVgy++w6io+0bn4iIJGRVQrt//35Kliz5yHNKlizJ/v37Y/crVKjAv//+a83tRERSVfHi5pK5b71lTlUAuHIF3ngDGjQwKyeIiEjaYVVCmylTJoKCgh55zsGDB+OtwRsVFUWWLFmsuZ2ISKpzd4fx4+HwYejQIa59yxazxFeTJrBmjf3iExGROFYltI0aNWLt2rV89tlnRERExDsWERHB559/zrp162jSpEls+7Fjxyiox4VFJJ159llYvBg2bICYhQgNAwIDoVUreOEF86EyERGxH6sS2nHjxpEnTx6GDBlCwYIFad26Nb169aJ169YUKlSI999/n9y5c/PZZ58BcOXKFQ4ePEjr1q1tGryISGpp0MAcrf30UyhaNK590SIoUgRefRUOHbJffCIijsyqhLZQoULs27ePzp07c+vWLdasWcOMGTNYs2YNN2/epHPnzuzdu5dChQoBkDdvXq5du8bHH39s0+BFRFJT1qzwwQdw8iTMmWOW/QIICzMXb6hUCUaMcCIiIml1ukVExDasXvrW19eX2bNnM336dE6cOEFISAheXl6ULFmSTJky2TJGEZE0xckJXn4ZGjUyF2yYNctcRjcqCj791Jn8+evzzz9OdOtmJsEiIpKyrBqhfVCmTJkoW7YstWrVomzZskpmRcRh5MkDkybBpUswciS4/P8QwcWLnvTt60zJknDggH1jFBFxBE+c0IqIOLqsWWHYMNi1C2rUiCtUe+kS1K0LP/9sx+BERBxAkqYcNGjQAIvFwo8//kj+/Plp0KBBki5usVjYsGHDEwUoIpJeVK4MW7ZE8fXXW1m8uC47dzoRGgotW8Lzz8Mnn5irlImIiG0lKaHdvHkzFouFsLCw2P2ksFj0YISIOJ5ixYJZty6Knj2dWLLEbFu8GJYsgeeeg0GDoE6duEUbRETkySRpykF0dDRRUVGU+P+hhejo6CS9oqKiUjR4EZG0ysMDFi6Eb78159qCWb92+XKoVw+qVDErJdy/b9cwRUQyBM2hFRFJIU5O0Ls3nDoFY8aAr2/csQMHoEsXs4btmDFw44b94hQRSe+sSmidnZ15+eWXbR2LiEiGlDUrvP8+nD1rjspWrhx37NIlGDIE8ueHPn3g+HH7xSkikl5ZldB6eXlRoEABW8ciIpKhZcpk1q/duxe2boW2bePm0d69a05PeOYZaNEC1q83pyiIiMjjWZXQVq1ald9//93WsYiIOASLxXwobNkyc9Wx/v3jL8Cwdi00bQrlysH8+eaCDSIi8nBWJbQjRoxg48aNzJo1y9bxiIg4lGLF4Kuv4K+/4Isv4P9XDAfgyBF46SWz1Nd778GhQ/aLU0QkLbNq6dvAwED8/f3p0aMHX3/9NX5+fuTJkydBmS6LxcLQoUNtEqiISEaWLRsMHmyO1q5YAePHmws1AJw5A+PGweefm8lvv352DVVEJM2xKqEdMWJE7Pb+/fvZv39/oucpoRURSR4XF+jQAdq3hw0b4LPPYONGiI4259T27w8XLsCnn4Krq72jFRFJG6xKaDdt2mTrOERE5AEWCzRqZL7+/dccsR03zjw2fjysWweTJ5tL64qIODqrEtp69erZOg4REXkIHx9zpLZoUXjzTXO09sgRc4GGV14xpyLkzWvvKEVE7EcLK4iIpBO9e5vzaqtUiWubMwdKlYJVq+wXl4iIvdkkoY2Ojuann35i8ODBDBw4kKVLl9risiIi8h9+fmZS++23kD272RYcDM89B8OHw82b9o1PRMQekpXQTpkyhcKFC5M1a1bq1avHkSNHCAsLo1atWrzyyit8+eWXTJw4kY4dO9KqVSuio6NTKm4REYfl7GyO1v75p/kAGZgPjI0aBXnywAsvwNWr9o1RRCQ1JTmhXblyJQEBAVy4cAF3d3d+++032rVrx9ixY9mzZw+dO3dm4sSJjBw5koIFC7J27Vq+/fbblIxdRMSh5coFixaZFQ9iqiZGRJhtNWvCqVP2jU9EJLUkOaGdPHkyHh4e7N27l2vXrrF3714uX77Ml19+yYgRI5g9ezZ9+/Zl6NChHDhwgOzZszN37tyUjF1ExOFZLPDBBxAUBIMGmUkuwOnTUKOGOTUhIsKuIYqIpLgkJ7QHDhygTZs2VK5cGYDKlSvTunVrQkND6d69e7xzs2fPTps2bTh69KhNgxURkcSVKwf/+x8cPAhlypht165Bnz7wzDOwd6994xMRSUlJTmhv3rxJsWLF4rUVLVoUgKeeeirB+fnz5+f27dtPGJ6IiCRH/vzw229xc2vBHK2tUwd++MF+cYmIpKQkJ7TR0dFkypQpXlvMvpNTwss4Ozs/YWgiImKNbNlg8WKzGkL16mZbeDj06mVWQ7hwwa7hiYjYnOrQiohkUNWqwZYtEBAQ17ZypTkFYeBAJbYiknEka6WwI0eOsHDhwnj7AIsWLcIwjATnioiIfWXKBJMmmauK9e8PV65AWBh89RVMnGhWQ2jXDtq2hf/MKhMRSTeSldAuWbKEJUuWxO7HJLEvvvhignMNw8ASU0dGRETsqmNHaNwYRoyAqVPh7l2zdu327ebr7behfHkYPBheeglcXe0dsYhI0iU5oR0+fHhKxiEiIiksWzb48kv48EOYMgUWLIBjx+KO//47dOsGQ4fCm29Cz57g42OvaEVEkk4JrYiIg/HxgWHDzNeff8Ly5eZDZDGlvS5cgPffN0dzBw82tz097RmxiMij6aEwEREHVqIEvPsu7NkDW7dCixZxx+7dM1chK1nSPCYiklYpoRUREcCsVbtmDZw8CQMGxM2jvXzZnH87f7594xMReRgltCIiEk/x4uZc22PHoFEjs+3+ffNhsZdeMqcpiIikJUpoRUQkUcWLw9q18OqrcW3z50Pp0uY0hbAw+8UmIvIgJbQiIvJQLi5mma/JkyFXLrMtKgo+/9ws87Vnj33jExEBJbQiIvIYFotZxuvMGbPyQcwq6KdOQe3a8M03Zk1bERF7UUIrIiJJ4ukJw4eb9WqrVzfbIiLMpXUbN4ajR+0bn4g4ridOaI8dO8bSpUuZPXu2LeIREZE0rlQps4zXoEFxbRs2mFMQPv88aaO1f/8Nq1fD9espF6eIOA6rE9q9e/dSoUIFypYtS8eOHenevXvssa1bt5I5c2ZWrlyZ7OuOGTMGPz8/PD09yZ07N23btuXEiRPxzunduzfFihXDw8MDHx8fnnvuOY4fP27tRxERkWRydYX//Q9WroTChc22qCjzYbFu3eDOnYe/999/oUoVaN0a8uWDtm3h7NnUiFpEMiqrEtqjR4/SoEEDzp49y6BBg2jevHm843Xq1CFXrlwsWrQo2dfesmULAQEB7Nq1i8DAQCIiImjSpAmhoaGx51SuXJkZM2bwxx9/8Msvv2AYBk2aNCEqKsqajyMiIlZq3dos7/XBB3Fts2dD0aIwYQKEhyd8z/vvw5Ur5nZEBKxYAW+9lTrxikjGZFVCG7MM7v79+xk/fjx+fn7xjlssFmrUqMHemHUUk2HdunV0796dZ599lvLlyzNz5kwuXLjA/v37Y895/fXXqVu3LoULF6ZSpUqMHj2av/76i3PnzlnzcURE5Al4eJgrii1aBJkzm23//msum1ulCjzwf9/s2AE//GBuZ80ad/6mTRAdnbpxi0jG4WLNm7Zs2UKHDh0oXrz4Q88pWLAg69atszqwGMHBwQDkyJEj0eOhoaHMmDGDIkWKUKBAgUTPCQ8PJ/yBYYKQkBAAIiIiiIiIeOIYHyfmHqlxL0k71O+OyZH7/bnnYO9eGDnSmYULzfGSI0fAz8+gZEnIm9fgwAELYAHg44+j2LLFwvLlTty6BQcPRlCunP3ifxKO3O+OSn2e8pLz3VqV0N6+fZvcuXM/8py7d+8+8RSA6OhoBg4cSK1atShTpky8Y9988w3vvvsuoaGhlCxZksDAQDLF1JL5jzFjxjBy5MgE7evXrydzzPBAKggMDEy1e0naoX53TI7c7507Q7Vq3kyaVIGzZ7NhGBaOH4fjxy2x5xQpcouCBbeSM2cRoCwA3333By1bpu/JtI7c745KfZ5ywpKxeovFMJJfPfDpp5+mePHirF27FoCRI0cyatSoeAlszZo1CQ0N5ffff0/u5WP16dOHtWvXsm3bNvLnzx/vWHBwMFevXuXy5cuMHz+ev//+m+3bt+Pu7p7gOomN0BYoUIBr167h5eVldXxJFRERQWBgII0bN8Y1ZnF0yfDU745J/R4nIgK+/NKJxYudOHIEIiIs5M5tUK2awaefRlGyJBw8CNWqmd9Thw7R/PRT+nwWQv3ueNTnKS8kJIRcuXIRHBz82HzNqhHaVq1aMXHiRH799VcaxSz0/YCFCxeya9cuhg4das3lAejbty+rV69m69atCZJZAG9vb7y9vXn66aepXr062bNnZ9myZbz00ksJznVzc8PNzS1Bu6ura6r+EKb2/SRtUL87JvW7WQlhyBDzFR4OISGQK5cFi8VCzCMclSub9W1v34Zt25xwcXHCYnn0ddMy9bvjUZ+nnOR8r1Y9FDZkyBB8fX1p0aIFr732Gvv27QPMaQBdunShc+fOFC5cmMGDByf72oZh0LdvX5YtW8bGjRspUqRIkt5jGEa8UVgREUk73NzAx4cEyaqzM9SqZW7/84+5+piISHJZNULr4+PDli1b6NKlC9OnT49t79u3LwDVqlXjp59+wtvbO9nXDggIYN68eaxYsQJPT0+u/H9tF29vbzw8PDhz5gwLFiygSZMm+Pj4cPHiRcaOHYuHhwctWrSw5uOIiIgd1a0LMc8Qb90KTz9t33hEJP2xKqEFKFq0KNu3bycoKIhdu3Zx48YNvLy8qFatWoIyXskxZcoUAPz9/eO1z5gxg+7du+Pu7s5vv/3Gl19+yc2bN8mTJw9169Zlx44dj31QTURE0p66deO2162DXr3sF4uIpE9WJ7QxKlSoQIUKFWwQiulxz6j5+vry888/2+x+IiJiX1WqQPbscPMmLF4Me/ZA1ar2jkpE0hOr5tAGBwdz6NChh5ZTCA0N5dChQ7H1XkVERB7GzQ3+f70eAPr21SILIpI8ViW0o0aNolatWg+tMxsVFUWtWrX45JNPnig4ERFxDG++Cc8+a27v3QsPPJ4hIvJYViW069ato3Hjxnh6eiZ63MvLi6ZNm2pqgIiIJImrK3z9ddz+W2/B6dP2i0dE0herEtoLFy7w9GMeQy1WrBgXLlywKigREXE89etD9+7m9u3b8PLL5uIMIiKPY1VCa7FYHlvzNTw8/ImXvhUREccycSIUK2Zu795tVjy4d8++MYlI2mdVQluqVCnWrVv30IoE0dHRrF27lpIlSz5RcCIi4lg8PWHePHD5/xo8s2dDzZpw8qR94xKRtM2qhPall17izz//pGfPngQHB8c7FhwcTM+ePTl16hSvvPKKTYIUERHHUbUqzJkDHh7m/sGDULYsfPop3L9v39hEJG2yKqHt27cvNWvW5Mcff6RIkSI0bdqUnj170rRpU4oUKcKsWbOoU6dO7MphIiIiydGpkznlIOZxjfBw+PBDqFwZdu60b2wikvZYldC6urry66+/MnjwYKKioggMDGTmzJkEBgYSHR3NO++8wy+//IKrq6ut4xUREQdRtqw5Ovv22+D0//+1OnIEatWCF14wE14REbAyoQVwd3dn/Pjx3Lhxg8OHD7Nt2zaOHDnC9evX+eyzz3Bzc7NlnCIi4oCyZIHPPzdr01aubLYZBixaBNWrw8CB5r6IOLYnXvrW2dmZZ2OqYYuIiKSASpVg1y6YPBnGjIF//jHbv/oKcuWCjz6yb3wiYl9Wj9CKiIikJhcXGDAAzp2D//0vrn3oUJg61W5hiUgaYPUI7e3bt5k0aRK//vorly5dSrQurcVi4bSWehERERtyd4dBgyAyEt5912x74w0z4e3Z076xiYh9WJXQ/vvvv9SsWZPTp0/j5eVFSEgI3t7e3L9/n7t37wLg6+urh8JERCTFvP22OfXgiy/MebSvvmo+PBaz2piIOA6rphyMGDGC06dPM2vWLG7evAnAoEGDCA0NZffu3VStWpXChQtz9OhRmwYrIiISw2IxHxgbONDcj0lq1661a1giYgdWJbQ///wzDRs25JVXXsFiscQ75ufnx9q1azl37hwjR460SZAiIiKJsVjM+bQxZc+joqBjR9iyxb5xiUjqsiqhvXz5MhUrVozdd3Z2jp1qAJA9e3aaN2/OwoULnzxCERGRR7BYzGoHzz9v7oeGgr8/tGsHeoxDxDFYldB6e3sTERERu589e3YuXrwY7xwvLy/+iamrIiIikoKcnGD2bDORjbF8Ofj5abRWxBFYldAWLVqUc+fOxe5XrFiRwMBArl+/DsDdu3dZtWoVBQsWtEmQIiIij+PuDuvXw3ffga+v2XbzJjRuDHPm2Dc2EUlZViW0TZo0YcOGDYSFhQHQu3dvrl69Svny5enYsSNlypTh9OnTdNejpiIikopcXeH11+GPP6B5c7MtIgK6dIFvvrFvbCKScqxKaN944w2mTZsWm9C2b9+ezz//nNDQUJYsWcKVK1cYPHgw77zzjk2DFRERSQovL1i5Evr0iWsLCIBevSAoyG5hiUgKSVJCGxISwv3792P38+XLR6dOnciVK1ds21tvvcW1a9e4fPkyd+7c4fPPP8fZ2dn2EYuIiCSBi4u5VO6QIXFtP/wAFSuaI7aJrAckIulUkhLa7Nmz89lnn8Xu9+zZk5UrVyY4z9nZmTx58iQo5SUiImIPFgt88olZ2itr1rj2OXOgZUsICbFfbCJiO0lKaC0WC9HR0bH7M2fOJEh/sxERkXRi0CD4+2/4+mvw8DDbNmyAkiVhwgR4oPKkiKRDSUpofX19OXXqVErHIiIikmK8vMwFGDZuhBw5zLYrV2DwYChTBjZtsm98ImI9l6ScVL9+febOncu1a9fIly8fAMuXL49XuisxFouF6dOnP3GQIiIitlK9OuzZA++/D4sXm21nzkCDBuaDY+PHmyXARCT9SFJCO27cOP755x8CAwOJjo7GYrEQFBT02GkHSmhFRCQtKlYMFi2CQ4fMJHbbNrN98mT47Tfz4bHKle0bo4gkXZIS2jx58rBu3ToiIiK4fPkyhQsXZuDAgQwYMCCl4xMREUkx5cqZK4lNmQLvvGPOpT10CKpUgVq1YPTo+KuPiUjalKSENoarqysFCxakXr16VKhQgUKFCqVUXCIiIqnCyckcpfX3h06d4OhRs337dmjY0KyS8N57ZsUEEUmbrFpYoXDhwrHL3IqIiGQEzz4L+/aZo7WlS5tt0dHwwQdQqhS8+SYcO2bfGEUkcVYltD/99BNXr161dSwiIiJ25e4Ob7wBhw/DyJFxo7J//mkmulWrwtat9o1RRBKyKqEtVqwYly9ftnUsIiIiaYKTEwwbBj//DDVqQMzCl6Gh0Ly5WfpLRNIOqxLanj17smbNGv7++29bxyMiIpJmNGsGO3bAjRvQooXZFhYGjRubUxG0fK5I2mBVQtuhQweqVatGzZo1mTx5Mnv27OH8+fNcuHAhwUtERCS98/KCpUvN5XLBnFs7dqxZDeHAAfvGJiLJrHIQo2jRolgsFgzDoH///g89z2KxEBkZaXVwIiIiaYWbGyxfDuPGwYgREBEBR45AtWrw1ltOVKjgbO8QRRyWVQlt165dsah+iYiIOBgXFxgyBFq1gm7dICgIIiPhs8+cyZ69EX/+6UTLllCxojkPV0RSh1UJ7cyZM20choiISPpRrhzs3m0uvDB2rDlae/OmO8OHw/DhULu2OUXBx8fekYo4Bv37UURExAqZMsGoUeZCDG3aRMc7tm2budLY2bN2Ck7EwSihFREReQJPPw2LF0fx7beBfP11FL6+ZvvJk1CpEnz7rfkQmYikHKsfCksKi8XC6dOnrbmFiIhIupI3bxgtWkTTurUzTZvCiRNw6xb06QMzZ5qJbYUKdg5SJIOyaoQ2OjoawzASvG7dusW5c+c4d+4c4eHhROufpCIi4mAKFTJr13brFte2ezdUrgzvvqvatSIpwaqE9ty5c5w9ezbB68aNG5w5c4a2bdtSuHBhjh49aut4RURE0rwcOcxR2U2boFQpsy06Gj7/3Fx57MQJu4YnkuHYfA5t4cKFWbBgATdv3uTDDz+09eVFRETSDX9/+P13sxpCpkxm28GDZpWEIUPg9m27hieSYaTIQ2Gurq40btyYhQsXpsTlRURE0o1MmeDDD2HXLihZ0my7fx/GjIG8eeH5582RXBGxXopVOQgLC+PGjRspdXkREZF0pWJF2L8f3nsvbrQ2LAyWLIGGDWHKFPvGJ5KepUhC+9tvv/HTTz9RMuafoiIiIkKWLOZCDEePwmuvQa5cZrthwJtvQpcuZjWEM2fsG6dIemNV2a4GDRok2h4ZGcnff//NuXPnABg2bJjVgYmIiGRUxYvD1KnwzTcwdKiZ5ALMmWO+nJ2hVy/o0MF8mKxECUhixUwRh2RVQrt58+ZE2y0WC9mzZ6dJkyYMHjyYxo0bJ/vaY8aMYenSpRw/fhwPDw9q1qzJZ599luhor2EYtGjRgnXr1rFs2TLatm2b7PuJiIjYi4uLOZc2X774Jb2iosyEd+rUuHOfeQYGDjRHdi0Wu4QrkmZZldCmZH3ZLVu2EBAQgJ+fH5GRkQwZMoQmTZpw7NgxsmTJEu/cL7/8Eot+q0VEJJ3r3x+6doVDh2DDBpgwIWEFhD/+gN69wcsLXnzRPnGKpFVWJbQpad26dfH2Z86cSe7cudm/fz9169aNbQ8KCuKLL75g37595MuXL7XDFBERsals2aBuXfPVty/Mng03b5qjtVu2mIs1AAQEmOXA8ua1Z7QiaYvNEtrIyEgOHz4MQJkyZXB1dbXJdYODgwHIkSNHbFtYWBidO3dm8uTJ5E3Cb3R4eDjhDyzNEhISAkBERAQRERE2ifNRYu6RGveStEP97pjU747J1v2eLRv06xe/rXNnZxYvduLGDejaNZpZs6LImdMmtxMr6Hc95SXnu7UYhmEk5cSzZ8+yadMmateuTYkSJeIdW716Nb169eLatWsAZM+enW+++YYXXnghGWEnFB0dTZs2bbh16xbbtm2Lbe/duzdRUVF8//335oewWB45h3bEiBGMHDkyQfu8efPInDnzE8UoIiKSGkJCMtGvX32Cg90B8PCIoH37U7RrdxIXlyT9p1wkXYkZwAwODsbLy+uR5yY5oR0yZAifffYZZ86coVChQrHtp06doly5cty7d49ChQqRJUsW/vjjD5ycnNizZw8VK1a0+oP06dOHtWvXsm3bNvLnzw/AypUreeuttzh48CBZs2Y1P8RjEtrERmgLFCjAtWvXHvsF2UJERASBgYE0btzYZiPXkvap3x2T+t0xpVa/r19v4fnnnbl3L+75kYoVDWbMiKR06RS7rSRCv+spLyQkhFy5ciUpoU3ylINt27ZRoUKFeMkswFdffcW9e/cICAjg66+/BmD58uW0b9+eSZMmMX36dCs+AvTt25fVq1ezdevW2GQWYOPGjZw+fZps2bLFO79Dhw7UqVMn0QoMbm5uuLm5JWh3dXVN1R/C1L6fpA3qd8ekfndMKd3vLVvCn3/CJ5/A99+b82sPHrRQubIrXbuaJcCKFEmx20si9LuecpLzvSZ5YYWzZ89StWrVBO3r1q0jU6ZMfPrpp7Ftbdu2pU6dOvz2229JDiSGYRj07duXZcuWsXHjRor85zfz/fff59ChQwQFBcW+ACZMmMCMGTOSfT8REZH0pEABc/GFXbugVCmzLSoKZsyAsmUhMNC+8YnYQ5JHaP/9919yxSxp8v9u3LjB6dOnqVOnDp6envGOVaxYkX379iU7oICAAObNm8eKFSvw9PTkypUrAHh7e+Ph4UHevHkTfRCsYMGCCZJfERGRjKpKFThwAMaPhy++gOBgCA2FVq1g/nxo187eEYqkniSP0Lq6unL9+vV4bfv37wegSpUqCc7/b83YpJoyZQrBwcH4+/uTL1++2NeCBQusup6IiEhG5eFhTjM4dw5iHiO5fx/atzcT24MH7RmdSOpJ8ghtiRIl2LBhQ7y29evXY7FYqFmzZoLzL126ZFV92CQ+o/bE7xEREckosmWDRYugRw9z6VyANWtg3TpztbGePe0ankiKS/IIbYcOHTh58iRvvPEGhw4dYvHixUydOpWsWbPSrFmzBOdv376d4sWL2zRYERERSZyLC/z4I0ybZs6zBXNuba9e8NFH5sitSEaV5IR24MCBlC1blqlTp1KxYkU6derE7du3GTlyZILpBfv27ePUqVM0btzY5gGLiIhI4pyc4NVX4eRJczndGJ98Yj5A9tVXsH49/P/6QiIZRpKnHGTOnJnt27czYcIEdu3aRc6cOenYsSOtW7dOcO6BAwd47rnnaNOmjU2DFRERkcdzczOT10KF4O23wTDg7FkYONA8niuXOR2hcmW7hiliM8la+jZr1qwMHTr0see9/vrrvP7661YHJSIiIk9u8GDw94f33oNff41rv3bNfIhs715IwgryImlekqcciIiISPpTqZJZm3bPHvjuO7PcF8DFi2ZSe/WqXcMTsQkltCIiIg7Azw9efx1Wr457aGz3bnj2WViwwJyWIJJeKaEVERFxIHnywMqV5jxaMKcfvPgidOgA/7+WkUi6o4RWRETEwVSoAEeOwPPPx7UtWwZlysDWrXYLS8RqSmhFREQcUJ485mIMCxeCj4/Zdv06NG4MP/1k39hEkksJrYiIiAPr2BGOHoWmTc39+/ehc2dz6dyjR+0bm0hSKaEVERFxcD4+sGqVuShDjDVrzCkILVvCf1a+F0lzlNCKiIgIrq4wdSrMmgVPPRXX/vPP0KgRdOkCN2/aLz6RR1FCKyIiIgBYLGbi+uef8PnnULBg3LE5c6BwYahXD4YPh9BQu4UpkoASWhEREYknc2ZzydzTp+GHH8Db22wPCTGrIIwaBdWqmYmvSFqghFZEREQS5eICPXrA4cPmg2K+vnHHjh6FypXhq68gMtJ+MYqAEloRERF5jAIFYO5c+PtvOHYMSpc22+/cgYEDoWxZGDoUNm6ES5e06pikPiW0IiIikmTPPGMumfv663Ftx4/D6NHQsKH5QFmRIrBunf1iFMejhFZERESSJWtW+O472LkTatRIePz8eWjeHPr2NSsnbNmS+jGKY3GxdwAiIiKSPlWvDjt2mNMM1q+H33+HXbvMF8DkyXHnTpwI/frZJ07J+DRCKyIiIk/E1xe6d4cJE8wE9/PPwdk5/jlvvw2LF5s1bUuXhtmzNddWbEcJrYiIiNiMxRJX8mvRInjlFbP9/n1zmd0NG+CPP6BrV6hZEyZNgnPn4t4fGgqXL9sldEnHlNCKiIiIzRUqBM8/D99/D+XLJ37Orl3mNIQiRaBtWzMRzpfPHPF96SU4ezbhe8LC4MaNFA1d0iEltCIiIpJi3Nxg3jzIndvc79sXVq6EUqXin7diBXzxBdy+be7Pnw/Fi5vJcNeuEBAArVpBjhyQMyf4+5tTGG7fhqgo8wG1VavMbYBff4X+/eHdd835u7dumccGDICnn4Zffkmtb0BSgx4KExERkRRVujScOAEREeDjY7a1agWHDsHq1TBlilnjFiBTJrOKwo0bEB1tnnPoUMJrbtlivpydwdPTTFgBXn3VTJqbN4+/4MPXX0PFiuY0CIDeveHUKXPxCEn/NEIrIiIiKS5btrhkFsy5tuXLw4cfmvNtZ82CsWPNaQZnz8Lw4WYC+t+Hy3x9zRHWGFFRccksmFMcWrVKuHrZqVNxySyYpcWWLLHVpxN7U0IrIiIiduXmBl26wHvvmQmrlxeMGAEHDkBwsLk62f79cOQIXLxoLuTwyy/w5pvmQg8+PlC3btz1Ll40/7d0aQgMfPgc3s8/h4ULzWkIO3bEtV+4YK6AVrNmyia9qvJgOxpoFxERkTQrSxYzaX2QxQJNmpivGIYBzz1nzqONOef7782FH7ZtMxPU336Djz82R4IPHjST5E6dzPMnTjTr6kZGQlBQ3Ajv88/DV1+Z83GfVHi4mbwDfPaZuVxwzpzmfGI/PzMpb9w47hxJOo3QioiISLpnscC0aeYIL8Bbb8WtYpY1q5ncnjgBL7xgVlNIzK5dsG9fwukKAwZAr15w9ap1se3cadbf9fAwqzesWAHvv2/OKb5yBTZvNkeLW7eGChXMEegFC8wE/auvEsYjCSmhFRERkQwhT5641crGjXv4eR07Qpky5naZMjBmDBQtau5bLGYZsWHDzKQzxg8/QIkS5uIRERGPjyUyEubONUdda9Y06+8ahlm9oW3buPOyZYv/vuPH4dln4cUXzWoQAweaI8ebNsVVcJCENOVAREREMoxcuczXo7i6mqOmJ06YI6LOzmZ5r8uXzfc++Cf/YsXM0d6QEHM+7+DBMHUq9O9vIXPmTBgGnDxpjqTu2WNe5/nnzTnB8+c/Oo5mzWDNGvj3X3NaxIgR5jzh6Oj45+3fDw0amKXPihc3qzr06BE3XUI0QisiIiIOKGtWqFw5roqCkxM89VTC+auvvgp//mlOObBYzLbjx+HNN13o3r05bm6ulCgBkyfD3r1mzdyffoqfzD7zjDkdYteuuEoPefLAzJnmffPkgQ4dzCT7pZfi3vPtt+aDbTGuXjUfXvvlF3j5ZfPhtUOHoGVLc+TYkR8yU0IrIiIi8gh58phzcPfsMacPPMrdu3HL/YL5sNnRo2ZiXK2aWbnhiy9g+3bzug/KmtVchOLaNfM9vXub58+dC+3amQ/IxYiKMpPYF16An382R47feCNuWsLx42YS/s03CUd8MyJNORARERFJgipVzKkB+/fDwoVRrFlzE0/PHGTK5ESzZmaiGrMgBJhzbvv0iRvZBcif30w+HyVnzrhtNzfo3Nl8gfkQWeHCZsWEL7+M/76pU+HSJXNhiVdeMRPjH34wE97ZsyF7dnNE98cfzXm5HTs+6TeSdiihFREREUkii8VMbMuXj6ZWre20aNECV1fzD95Zs5oVEWIMH277lcjy5jWnG/zwQ/x2FxfzQbTVq83Xg9asMZPgWrVg/fq4UdytW6FOHdvGZy+aciAiIiJiA336xM15LV8+5R7aGjQo/n7nzmb93QdXYouJIWa0NyQE1q6NXylh8OCMMx1BCa2IiIiIDbi6mjVlv//eXKHsv8v22kqZMtCihbnt4QGffmpWTDhxwpxu4OFhLtCwebM5B7dHD3P1NQB397havfv2mXN2MwIltCIiIiI24uNjPoz139FSW5s921z1bNMmKFTIbMueHb7+Gm7fNqcWZMsGBQua0xP++Qd+/RVOnTLn0MZ4/33z/PROCa2IiIhIOpMjB3z0kVk54b8SGxl2d4eGDc3SZI0amaW+AP7+26ydm94poRURERFxMBMnxpUB+/Zbc/Q2MX/8Ya5c9vzzaXulMiW0IiIiIg6maNH4ywP36gX37iU8b8gQOHYMliwxpzekVUpoRURERBzQG29A/frm9oULsHJl/ONXrpjVE2Ls2pV6sSWXEloRERERB+TkBB9+GLc/d2784zNnxp9msHt3qoRlFSW0IiIiIg7K3z+ujNfPP8P16+a2YZjlxx60a5fZnhYpoRURERFxUM7O8NJL5nZkJCxaZG6vXw+nT8c/99o1OHs2deNLKiW0IiIiIg7s5ZfjtqdNg+nToUOHuLYyZeK20+o8WiW0IiIiIg6sQoW4JXsPHIBXX4XQUHO/Zk1zJbIYaXUerRJaEREREQdmsZhL5v7XSy+Z9Wlr145r0whtEo0ZMwY/Pz88PT3JnTs3bdu25cSJE/HO8ff3x2KxxHu98cYbdopYREREJH174w3Yvh0mTIA+fWDqVLPqgYeHuaRuyZLmeUFBEB5u11AT5WLvAP5ry5YtBAQE4OfnR2RkJEOGDKFJkyYcO3aMLDFLWgCvvfYao0aNit3PnDmzPcIVERERSfcsFnN6Qc2aiR+vVg1OnID79+HgQahePXXje5w0l9CuW7cu3v7MmTPJnTs3+/fvp27durHtmTNnJm/evKkdnoiIiIjDqV4dZs0CV1c4c0YJbbIFBwcDkCNHjnjtc+fOZc6cOeTNm5fWrVszdOjQh47ShoeHE/7A+HhISAgAERERREREpFDkcWLukRr3krRD/e6Y1O+OSf3ueBytz1u1grJlLZQvb+DuDqnxsZPz3VoMI62WyIXo6GjatGnDrVu32LZtW2z71KlTKVSoEL6+vhw6dIj33nuPqlWrsnTp0kSvM2LECEaOHJmgfd68eZqqICIiIpIGhYWF0blzZ4KDg/Hy8nrkuWk6oe3Tpw9r165l27Zt5M+f/6Hnbdy4kYYNG3Lq1CmKFSuW4HhiI7QFChTg2rVrj/2CbCEiIoLAwEAaN26Mq6trit9P0gb1u2NSvzsm9bvjUZ+nvJCQEHLlypWkhDbNTjno27cvq1evZuvWrY9MZgGqVasG8NCE1s3NDTc3twTtrq6uqfpDmNr3k7RB/e6Y1O+OSf3ueNTnKSc532uaS2gNw6Bfv34sW7aMzZs3U6RIkce+JygoCIB8+fKlcHQiIiIiktakuYQ2ICCAefPmsWLFCjw9Pbly5QoA3t7eeHh4cPr0aebNm0eLFi3ImTMnhw4dYtCgQdStW5dy5crZOXoRERERSW1pLqGdMmUKYC6e8KAZM2bQvXt3MmXKxK+//sqXX35JaGgoBQoUoEOHDnz00Ud2iFZERERE7C3NJbSPe0atQIECbNmyJZWiEREREZG0Ls0tfSsiIiIikhxKaEVEREQkXVNCKyIiIiLpmhJaEREREUnXlNCKiIiISLqmhFZERERE0jUltCIiIiKSrimhFREREZF0TQmtiIiIiKRrSmhFREREJF1TQisiIiIi6ZqLvQOwB8MwAAgJCUmV+0VERBAWFkZISAiurq6pck+xP/W7Y1K/Oyb1u+NRn6e8mDwtJm97FIdMaG/fvg1AgQIF7ByJiIiIiDzK7du38fb2fuQ5FiMpaW8GEx0dzaVLl/D09MRisaT4/UJCQihQoAB//fUXXl5eKX4/SRvU745J/e6Y1O+OR32e8gzD4Pbt2/j6+uLk9OhZsg45Quvk5ET+/PlT/b5eXl76oXdA6nfHpH53TOp3x6M+T1mPG5mNoYfCRERERCRdU0IrIiIiIumaEtpU4ObmxvDhw3Fzc7N3KJKK1O+OSf3umNTvjkd9nrY45ENhIiIiIpJxaIRWRERERNI1JbQiIiIikq4poRURERGRdE0JrYiIiIika0poU8HkyZMpXLgw7u7uVKtWjT179tg7JLGhrVu30rp1a3x9fbFYLCxfvjzeccMwGDZsGPny5cPDw4NGjRpx8uRJ+wQrNjFmzBj8/Pzw9PQkd+7ctG3blhMnTsQ75969ewQEBJAzZ06yZs1Khw4d+Oeff+wUsdjClClTKFeuXGwh/Ro1arB27drY4+rzjG/s2LFYLBYGDhwY26Z+TxuU0KawBQsWMHjwYIYPH86BAwcoX748TZs25erVq/YOTWwkNDSU8uXLM3ny5ESPjxs3jokTJ/Ltt9+ye/dusmTJQtOmTbl3714qRyq2smXLFgICAti1axeBgYFERETQpEkTQkNDY88ZNGgQq1atYtGiRWzZsoVLly7Rvn17O0YtTyp//vyMHTuW/fv3s2/fPho0aMBzzz3H0aNHAfV5Rrd3716+++47ypUrF69d/Z5GGJKiqlatagQEBMTuR0VFGb6+vsaYMWPsGJWkFMBYtmxZ7H50dLSRN29e4/PPP49tu3XrluHm5mb89NNPdohQUsLVq1cNwNiyZYthGGYfu7q6GosWLYo9548//jAAY+fOnfYKU1JA9uzZje+//159nsHdvn3bePrpp43AwECjXr16xoABAwzD0O96WqIR2hR0//599u/fT6NGjWLbnJycaNSoETt37rRjZJJazp49y5UrV+L9DHh7e1OtWjX9DGQgwcHBAOTIkQOA/fv3ExEREa/fS5UqRcGCBdXvGURUVBTz588nNDSUGjVqqM8zuICAAFq2bBmvf0G/62mJi70DyMiuXbtGVFQUefLkideeJ08ejh8/bqeoJDVduXIFINGfgZhjkr5FR0czcOBAatWqRZkyZQCz3zNlykS2bNninat+T/8OHz5MjRo1uHfvHlmzZmXZsmWULl2aoKAg9XkGNX/+fA4cOMDevXsTHNPvetqhhFZE5AkEBARw5MgRtm3bZu9QJBWULFmSoKAggoODWbx4Md26dWPLli32DktSyF9//cWAAQMIDAzE3d3d3uHII2jKQQrKlSsXzs7OCZ52/Oeff8ibN6+dopLUFNPP+hnImPr27cvq1avZtGkT+fPnj23Pmzcv9+/f59atW/HOV7+nf5kyZaJ48eJUrlyZMWPGUL58eb766iv1eQa1f/9+rl69SqVKlXBxccHFxYUtW7YwceJEXFxcyJMnj/o9jVBCm4IyZcpE5cqV2bBhQ2xbdHQ0GzZsoEaNGnaMTFJLkSJFyJs3b7yfgZCQEHbv3q2fgXTMMAz69u3LsmXL2LhxI0WKFIl3vHLlyri6usbr9xMnTnDhwgX1ewYTHR1NeHi4+jyDatiwIYcPHyYoKCj2VaVKFV5++eXYbfV72qApByls8ODBdOvWjSpVqlC1alW+/PJLQkND6dGjh71DExu5c+cOp06dit0/e/YsQUFB5MiRg4IFCzJw4EBGjx7N008/TZEiRRg6dCi+vr60bdvWfkHLEwkICGDevHmsWLECT0/P2Lly3t7eeHh44O3tTa9evRg8eDA5cuTAy8uLfv36UaNGDapXr27n6MVaH3zwAc2bN6dgwYLcvn2befPmsXnzZn755Rf1eQbl6ekZOzc+RpYsWcj5f+3deUxUZxcG8GeAYRgEWQRHIaggROsG1kqIFBFFsQpFJaA1VbRtUKRRCy4Ri4Jt3etGaq1tXTBNajGE1ra2buOG0iouiSbQAZdGEQWliIg64Pn+MEwc7yiYfopjn19CAud973vPnZuQkzvvnGnXzhTnfX9JtHabhf+C7Oxs6dSpk9jb20twcLAUFha2dkr0f6TX6wWA4icxMVFEHrbuysjIEJ1OJxqNRoYMGSIlJSWtmzT9K5buNwDZvHmzaU59fb1MmzZN3NzcxNHRUUaPHi1Xr15tvaTpX3vvvfekc+fOYm9vL56enjJkyBDZvXu3aZz3/L/h0bZdIrzvLwuViEgr1dJERERERP8a99ASERERkVVjQUtEREREVo0FLRERERFZNRa0RERERGTVWNASERERkVVjQUtEREREVo0FLRERERFZNRa0RERERGTVWNASETXj4sWLUKlUmDRpkll80KBBUKlUrZOUBfn5+VCpVDh69Ghrp/JSMBqN8PPzQ0JCQmunQkTPGQtaInrpNRWUj/6o1Wp4e3sjISEBJ06caO0UW53RaMScOXMQFRWFAQMGtHY6LwW1Wo358+cjNzcXhYWFrZ0OET1Hdq2dABFRS3Xt2hXvvvsuAKCurg5FRUXIzc1Ffn4+9u7di4EDB77QfHJycnDnzp0Xes4n2bZtGwwGAzZs2NDaqbxUEhMTkZ6ejoyMDOzZs6e10yGi54QFLRFZDX9/f2RmZprFli5dinnz5iEjIwMHDx58ofl06tTphZ7vab788kv4+PggIiKitVN5qdjZ2WHcuHHIzs5GaWkp/P39WzslInoOuOWAiKza+++/DwAoKipSjG3atAmxsbHo0qULHBwc4O7ujqioKOj1eotrNTY2YtmyZfD394eDgwP8/f2xZMkSPHjwwOJ8S3toMzMzoVKpcODAAcX8LVu2QKVSYcuWLWZxvV6Pt956C15eXtBoNNDpdAgLC8PGjRtb8AoAZ8+exYkTJxAXF6fI5/79+8jOzkZUVBR8fHyg0WjQvn17jBkzBqdOnXpqjjt37kRoaCicnZ3RpUsXszVXr16N/v37w9nZGU5OTujRowdSU1NRXV1tmmcwGDB58mT4+vpCo9HA3d0dgYGBmDlzJkTE7Ly1tbVYuHAhevbsCa1WC1dXV0RFReHIkSMWr7m2thZZWVno06cPHB0d4eLigr59+yIjIwNGo9FsbkJCAkQEW7dubdHrSUTWh09oieiVYGen/HeWkpKCwMBAREZGwtPTE1euXEF+fj4iIyORl5eH2NhYs/lJSUnYtGkTfH19kZKSgrt372LVqlXP9UNWv/zyC2JiYuDq6orY2Fh07NgRlZWVOHPmDLZt24akpKRm19i3bx8AICQkRDF28+ZNzJw5E2FhYRgxYgTc3Nxw/vx5/PTTT9i1axcOHTqE/v37K47Lzc3F7t27ER0djWnTpuHWrVsAgPr6egwdOhQFBQUICAjA5MmTodFoYDAY8NVXX2HixIlwc3NDeXk5goODUVdXh5EjR2Ls2LGoq6uDwWDA+vXrsXLlStM9u3nzJgYOHIhz584hNDQUU6dOxa1bt/Djjz8iIiICubm5GDVqlCm369evIzw8HMXFxQgKCkJycjIePHiA4uJiLFu2DGlpaXB1dTXN79evH9RqNfbt24dPPvnkWW4PEVkLISJ6yV24cEEASFRUlGJs8eLFAkBGjhypGDt//rwiVl5eLl5eXhIQEGAW1+v1AkACAwPl9u3bpvjly5fFw8NDAEhiYqLZMeHh4fL4v9GFCxcKANHr9Ypzb968WQDI5s2bTbExY8YIADl9+rRiflVVlSJmSXx8vAAQg8GgGLt7965cvnxZET979qw4OTlJZGSkxRxtbGxkz549iuPS0tIEgEyYMEEaGhrMxv755x+pra0VEZF169YJAFmzZo1ijRs3bpj9PX78eAEgX3/9tVn82rVr4uPjI56enlJfX2+Kx8XFCQBJT09XrF1RUSFGo1ER79u3r6jVarl7965ijIisH7ccEJHVKC0tRWZmJjIzMzF79mwMHjwY6enp0Ol0WLFihWK+r6+vItaxY0fExcXBYDDg0qVLpnhOTg4AYMGCBWjTpo0p7u3tjRkzZjyHqzGn1WoVsXbt2rXo2MuXLwMAdDqdYkyj0cDb21sR79mzJyIiInDo0CHFW/QAEBsbi8jISLNYQ0MDNm7cCBcXF6xduxa2trZm4y4uLnBycjKLWboud3d30+9VVVXYvn07Bg8ejA8++MBsXvv27TF79mxUVlZi7969AICKigrk5eWha9euiv3UwMPXwNLTep1OB6PRiOvXryvGiMj6ccsBEVmNsrIyZGVlmcU6dOiAw4cPW/ywz/nz57FkyRLs378fV65cwb1798zGy8vL0blzZwDAmTNnAABhYWGKdSzF/l/GjRuHvLw8hISEYPz48RgyZAjCwsLg4eHR4jVu3LgBW1tbODs7Wxw/ffo0li9fjiNHjqCiokJRwFZVVaFjx45mseDgYMU6xcXFqK2tRWRkJNzc3J6aU0xMDObNm4eUlBTs27cPw4cPR3h4OPz8/MzmHT9+HI2Njbh3757FAtVgMJjOHR0djRMnTkBEEBERAbVa/dQcHtVURFdVVcHHx6fFxxGRdWBBS0RWIyoqCr/99hsAoLKyElu3bsXcuXPx9ttv488//zR7OlhaWorg4GDcunULERERiImJQdu2bWFjY4MDBw7g4MGDZgVuTU0NbGxsLBaSlp58/r/Ex8cjPz8fq1atwoYNG/DFF19ApVIhIiICn3/+OYKCgppdQ6vVorGxEUajUVHkHT16FIMHDwYADBs2DAEBAXBycoJKpUJ+fj7OnDmjKPQBy9dcU1MDABaf+D6uS5cuKCwsRGZmJn799Vf88MMPAIDu3btj0aJFiI+PB/Bw/ywAFBQUoKCg4Inr1dXVPXMOj6qvrwcAODo6PtNxRGQdWNASkVXy9PTErFmzUFNTg08//RQff/wx1qxZYxpfvXo1qqursW3bNlPv2iZTp05VtPhycXHBgwcPUFVVBU9PT7Oxa9eutTgvG5uHO7kaGhoUY03F2ONiY2MRGxuL2tpaFBQUIC8vD99++y2GDx+O4uJisw84WdKU782bNxWF6GeffYZ79+7h8OHDePPNN83GCgsLTU+mH2fpG9Ca8rhy5cpT82nSq1cv7NixA0ajEUVFRdi1axfWrVuHsWPHwsvLC6GhoWjbti0AIC0tDStXrmx2zWfNoUlT4fz4vSWiVwP30BKRVUtPT4eXlxfWr1+PixcvmuJlZWUAoOhkICIWnwQGBgYCAA4fPqwYsxR7kqa34i0VXJbaZD3K2dkZw4cPx8aNGzFp0iRcu3YNf/zxR7Pn7N27NwCgpKREMVZWVgZ3d3dFMXvnzh2cPHmy2bUf1a1bN7Rt2xbHjx83a8/VHLVajZCQEGRlZWHdunUQEfz8888AgP79+0OlUuHYsWMtWuuNN96AjY0N9Hq9xb2/T1JSUgJvb2+z/btE9OpgQUtEVk2r1WLu3LkwGo1mLZma9sY+3sd06dKlOHv2rGKdCRMmAAAWLVpkensbeFiYrl27tsX5NLXAysnJMetfe+zYMXz33XeK+YcOHUJjY6Mi3vThJQcHh2bPGR4eDgAWi9/OnTujuroa586dM8UaGxsxa9YsVFZWNrv2o+zs7DBlyhTU1NRgxowZirxrampw+/ZtAA/7Aje1+npU09Pupuvq0KEDEhIScPToUaxYsULRn7bpupq+kU2n0yEuLs7ifmrg4ev2+NPxv//+GxUVFS/8m+SI6AVq3SYLRETNe1rbLhGR+vp68fLyEjs7OyktLRURkZMnT4parRatViuJiYmSmpoqAwYMEAcHBxk5cqTF1lqTJ08WAOLr6yupqamSkpIiHh4eEh0d3eK2XSIioaGhAkCCg4Nl1qxZEh8fL/b29jJ69GhF267AwEDx9PSU0aNHy0cffSSpqakSHBwsACQkJETRGsuS2tpacXZ2lqFDhyrGdu7cKQDE1dVVkpKSZPr06dKnTx9p166dDBo0SADIhQsXTPMttRZ7/LUOCwsTABIQECDTp0+X2bNnS1xcnLRp00ZOnTolIiIzZswQBwcHGTZsmCQnJ8vcuXMlJiZGbG1txd3dXS5dumRa88aNGxIUFCQApHfv3pKUlCRz5syRd955RwICAgSAXL161TS/srJSXnvtNQEgffv2lbS0NElNTZXo6Gixt7eX6upqs5y/+eYbASDbt29v9rUkIuvEgpaIXnrNFbQiItnZ2ab+qE30er2EhoaKs7OzuLq6yogRI6SoqOiJvWIbGhpkyZIl4ufnJ/b29uLn5yeLFy+W0tLSZypoq6qqZOLEieLu7i5arVZCQkLk999/t1gsfv/995KQkCBdu3YVR0dHcXFxkcDAQFm2bJmpp2tLJCcni62trZSXlyvGduzYIa+//ro4OjqKh4eHJCQkSFlZmSQmJj5zQSvysLftypUrJSgoSLRarTg5OUmPHj0kLS3NVEwWFhbKlClTpFevXuLq6iparVYCAgLkww8/NCtmm9y5c0eWL18u/fr1kzZt2ohWqxVfX18ZNWqU5OTkKHrL1tTUSEZGhnTv3l00Go24uLhIUFCQLFiwQO7fv282d9CgQdK+fXtFnIheHSoRC+/vEBGRVSkpKUGvXr2QmZmJ+fPnt3Y6Lw2DwYBu3bohMzMTCxYsaO10iOg5YUFLRPSKSE5ORm5uLi5cuPDEnrT/NRMmTMD+/fvx119/mX1hBhG9Wti2i4joFZGVlQWdToeLFy+aOh/8lxmNRnTr1g2TJk1iMUv0iuMTWiIiIiKyamzbRURERERWjQUtEREREVk1FrREREREZNVY0BIRERGRVWNBS0RERERWjQUtEREREVk1FrREREREZNVY0BIRERGRVWNBS0RERERW7X+m8k2EvEoIUgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from astropy.io import fits\n",
    "from astropy.wcs import WCS\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "with fits.open('frame-r-004874-3-0692.fits') as hdul:\n",
    "    data = hdul[0].data\n",
    "    header = hdul[0].header\n",
    "\n",
    "data = np.nan_to_num(data, nan=0)\n",
    "wcs = WCS(header)\n",
    "ra, dec = 40.669622, -0.013294\n",
    "\n",
    "center_pixel = wcs.world_to_pixel_values(ra, dec)\n",
    "x_center, y_center = int(center_pixel[0]), int(center_pixel[1])\n",
    "\n",
    "cutout_size_sdss = 200\n",
    "x_min, x_max = x_center - cutout_size_sdss, x_center + cutout_size_sdss\n",
    "y_min, y_max = y_center - cutout_size_sdss, y_center + cutout_size_sdss\n",
    "cutout_data = data[y_min:y_max, x_min:x_max]\n",
    "\n",
    "y, x = np.indices(cutout_data.shape)\n",
    "center = np.unravel_index(np.argmax(cutout_data), cutout_data.shape)\n",
    "r = np.sqrt((x - center[0]) ** 2 + (y - center[1]) ** 2).astype(int)\n",
    "\n",
    "tbin = np.bincount(r.ravel(), cutout_data.ravel())\n",
    "nr = np.bincount(r.ravel())\n",
    "radial_profile = tbin / nr\n",
    "\n",
    "\n",
    "pixel_scale = 0.16\n",
    "r_arcsec = np.arange(len(radial_profile)) * pixel_scale\n",
    "\n",
    "\n",
    "flux_profile = radial_profile * 3.6e-6  # 1 nanomaggy = 3.631e-6 Jy\n",
    "magnitude_profile = -2.5 * np.log10(flux_profile/3631)\n",
    "flux_tot = np.sum(flux_profile)\n",
    "mag_tot = -2.5 * np.log10(flux_tot/3631)\n",
    "print(\"星等：%s\" %mag_tot)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(r_arcsec, magnitude_profile, color='b', lw=2)\n",
    "plt.gca().invert_yaxis() \n",
    "plt.xlabel('Radius (arcsec)', fontsize=14)\n",
    "plt.ylabel('Surface Brightness (mag/arcsec²)', fontsize=14)\n",
    "plt.title('Surface Brightness Profile', fontsize=16)\n",
    "plt.grid(True)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于所有波段K改正为0mag的天体来说，其SED可以由什么函数描述f(λ)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "假设SED为幂律谱，$f(λ) ∝ λ^{−α}$ <p>\n",
    "$K(z)=2.5log_{10}(\\frac{L_{λ}}{(1+z)L_{λ^{'}}})$ <p>\n",
    "要使$K(z)=0$, 则$\\frac{L_{λ}}{(1+z)L_{λ^{'}}}=1$ <p>\n",
    "$λ^{'}=λ(1+z)$ <p>\n",
    "$L_{λ}=λ^{−α}$, $L_{λ^{'}}=((1+z)λ)^{−α}$<p>\n",
    "$\\frac{L_{λ}}{(1+z)L_{λ^{'}}}=(1+z)^{α-1}$<p>\n",
    "因此，$α=1$，即SED为$f(λ) ∝ λ^{−1}$时，K改正为0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 按照AB星等的定义，K改正为0，fv是常数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "某标准烛光B波段的星等测量值为25等（内禀绝对星等为-3.0），测量误差为0.1个星等，该天体方向上银河系尘埃的红化E(B-V)=0.2星等（不确定性为0.05个星等），请问该天体的距离为多少，距离测量的不确定性是多少？<p>\n",
    "假设银河系的消光曲线Rv=3.1<p>\n",
    "银河系尘埃消光轮子<p>\n",
    "• https://ned.ipac.caltech.edu/forms/calculator.html\n",
    "• https://github.com/kbarbary/sfdmap\n",
    "• https://github.com/gregreen/dustmaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "距离为2.73 Mpc\n",
      "距离测量不确定性：0.38 Mpc\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "Rv = 3.1\n",
    "E_bv = 0.2\n",
    "Av = Rv * E_bv\n",
    "Ab = 4.1/3.1 * Av\n",
    "\n",
    "Mag = -3.\n",
    "mag = 25 - Ab\n",
    "distance = 10 ** (0.2 * (mag -Mag)) * 10 / 1e6\n",
    "print(\"距离为%.2f Mpc\" %distance)\n",
    "\n",
    "dE_bv = np.array([-0.025, 0.025])\n",
    "Ab_d = 4.1/3.1 * Rv * (E_bv + dE_bv)\n",
    "dmag = np.array([0.05, -0.05])\n",
    "mag_d = 25 + dmag - Ab_d\n",
    "distance_d = 10 ** (0.2 * (mag_d -Mag)) * 10 / 1e6\n",
    "error_distance = distance_d[0] - distance_d[1]\n",
    "print(\"距离测量不确定性：%.2f Mpc\" %error_distance)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 误差是平方求和"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "某星系（Ra=180，Dec=-10），观测到其Ha谱线的（真空中）波长被红移到7000A处，在r波段观测到的其AB星等为20等，其光谱的能量分布fλ 在静止波长（5000-7500A）范围内近似为常数，请计算其r波段的绝对星等。<p>\n",
    "宇宙学（H0 ~100 h kms-1 Mpc-1 Ω0=0.28，ΩΛ=0.72）<p>\n",
    "r波段的波长范围近似为5500-6700A，滤光片响应曲线近似为常数0.6 <p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r波段上的绝对星等为-17.94\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "import astropy.units as u\n",
    "cosmo = FlatLambdaCDM(H0=70 * u.km / u.s / u.Mpc, Tcmb0=2.725 * u.K, Om0=0.28)\n",
    "\n",
    "wl_0 = 6563\n",
    "wl_obs = 7000\n",
    "z_rs = wl_obs / wl_0 - 1\n",
    "\n",
    "distance = cosmo.luminosity_distance(z_rs).to(u.pc).value\n",
    "\n",
    "#由于该星系在5000-7500A上SED近似为常数，因此r波段上K改正近似为0\n",
    "K = 0\n",
    "\n",
    "mag = 20 + 2.5 * np.log10(0.6)\n",
    "Mag = mag - 5 * np.log10(distance / 10) - K\n",
    "print(\"r波段上的绝对星等为%.2f\" %Mag)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 0.6为什么会其起作用"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
