{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "假设一颗大质量的恒星为40,000 K的O型星，其热光度为太阳的10万倍，假设周围是中性星际气体（全部为氢原子，n~10cm-3）均匀分布，其在10000年内能电离的球形空间的半径为多大？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000年内电离的球星空间的半径为16.52 pc\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "yr = 3.1536e7\n",
    "pc = 3.08568e18\n",
    "Lsun = 3.827e33\n",
    "Lstar = 1e5 * Lsun\n",
    "\n",
    "e_ion = 13.6 * 1.6e-12\n",
    "m_H =1.674e-24\n",
    "n_H = 10\n",
    "\n",
    "E_10000 = Lstar * 1e4 * yr\n",
    "N_H = E_10000 / e_ion\n",
    "R = pow(3/4 * N_H / (n_H * np.pi), 1/3) / pc\n",
    "print(\"10000年内电离的球星空间的半径为%.2f pc\" %R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#量级估计good\n",
    "#更好的是计算大于13.6ev的光子数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "请对银河系中的各种物质成分的质量比进行估算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "分子云：<1% <p>\n",
    "冷中性气体：1-5% <p>\n",
    "温中性气体:10%-20% <p>\n",
    "温电离气体：20%-50% <p>\n",
    "H Ⅱ：<1% <p>\n",
    "热电离气体：30%-70%"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下载一个LAMOST或者SDSS观测到的星系光谱文件，画出相应光谱，测量至少一个Lick指数特征及其误差，并在图中标注"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1sAAAIPCAYAAACFVMRQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB+O0lEQVR4nO3dd3gU1eLG8XdTSQJJgIQkQOhKL4KKQaUoEoqKyFUURFCuXhS8CoqIDRAFBMUuXgtgAdtV0YuCAtJUpGmoGqWjEKJAEpKQfn5/zC9LlgQIkN1Jst/P88zD7MzZ2TObYXfePWfOOIwxRgAAAACAMuVjdwUAAAAAoDIibAEAAACAGxC2AAAAAMANCFsAAAAA4AaELQAAAABwA8IWAAAAALgBYQsAAAAA3ICwBQAAAABuQNgCAAAAADcgbAEAAACAG/jZ+eJTpkzRp59+ql9//VVBQUHq1KmTnn76aTVt2tRZJisrS/fff78++OADZWdnKz4+Xq+++qqioqKcZfbu3au77rpLy5YtU9WqVTVkyBBNmTJFfn7Hd2/58uUaPXq0tm7dqtjYWD366KMaOnRoqepZUFCg/fv3q1q1anI4HGW2/wAAAAAqFmOMjh49qtq1a8vH5zRtV8ZG8fHxZvbs2WbLli0mISHB9O7d29SrV8+kp6c7ywwfPtzExsaapUuXmvXr15tLLrnEdOrUybk+Ly/PtGrVynTv3t38/PPP5quvvjIRERFm3LhxzjI7d+40wcHBZvTo0Wbbtm3mpZdeMr6+vmbRokWlque+ffuMJCYmJiYmJiYmJiYmJiPJ7Nu377Q5wmGMMSon/vrrL9WqVUsrVqxQ586dlZqaqsjISM2bN0//+Mc/JEm//vqrmjdvrtWrV+uSSy7RwoULdfXVV2v//v3O1q7XXntNY8eO1V9//aWAgACNHTtWX375pbZs2eJ8rZtuukkpKSlatGjRaeuVmpqq8PBw7du3T6Ghoe7ZecATcnOl2bOt+dtuk/z97a0PAABABZOWlqbY2FilpKQoLCzslGVt7UZ4otTUVElSjRo1JEkbNmxQbm6uunfv7izTrFkz1atXzxm2Vq9erdatW7t0K4yPj9ddd92lrVu36oILLtDq1atdtlFY5r777iuxHtnZ2crOznY+Pnr0qCQpNDSUsIWKLSNDGjPGmr/rLikkxN76AAAAVFClubyo3AyQUVBQoPvuu0+XXnqpWrVqJUlKSkpSQECAwsPDXcpGRUUpKSnJWaZo0CpcX7juVGXS0tJ07NixYnWZMmWKwsLCnFNsbGyZ7CMAAAAA71FuwtaIESO0ZcsWffDBB3ZXRePGjVNqaqpz2rdvn91VAgAAAFDBlItuhCNHjtSCBQu0cuVK1a1b17k8OjpaOTk5SklJcWndOnjwoKKjo51l1q5d67K9gwcPOtcV/lu4rGiZ0NBQBQUFFatPYGCgAgMDy2TfAAAAAHgnW8OWMUb33HOPPvvsMy1fvlwNGzZ0Wd+hQwf5+/tr6dKl6t+/vyQpMTFRe/fuVVxcnCQpLi5OTz31lJKTk1WrVi1J0uLFixUaGqoWLVo4y3z11Vcu2168eLFzGwAAAMDZMMYoLy9P+fn5dlcFZcjf31++vr7nvB1bw9aIESM0b948ff7556pWrZrzGquwsDAFBQUpLCxMw4YN0+jRo1WjRg2FhobqnnvuUVxcnC655BJJUo8ePdSiRQsNHjxY06ZNU1JSkh599FGNGDHC2To1fPhwvfzyy3rwwQd1++2369tvv9VHH32kL7/80rZ9BwAAQMWWk5OjAwcOKDMz0+6qoIw5HA7VrVtXVatWPbft2Dn0+8lG8Jg9e7bzhsOFNzV+//33XW5qXNhFUJL27Nmju+66S8uXL1dISIiGDBmiqVOnFrup8ahRo7Rt2zbVrVtXjz32WKlvapyWlqawsDClpqYyGiEqtowMqfBDIz2d0QgBADhLBQUF+v333+Xr66vIyEgFBASUanQ6lH/GGP3111/KzMzUeeedV6yF60yyQbm6z1Z5RdhCpZGXJ339tTUfHy/5lYvLNgEAqHCysrK0a9cu1a9fX8HBwXZXB2Xs2LFj2r17txo2bKgqVaq4rDuTbMCZFuBN/PykPn3srgUAAJWGj0+5GdwbZaisWik5OgAAAADADWjZArxJbq40d641P2iQ5O9vb30AAAAqMcIW4E1ycqTbbrPmb7iBsAUAAOBGdCMEAAAAvMxff/2lu+66S/Xq1VNgYKCio6MVHx+v77//3mN1cDgcmj9/vsdezw60bAEAAABepn///srJydHbb7+tRo0a6eDBg1q6dKkOHTpkd9Vc5OTkKCAgwO5qnDVatgAAAIAyYIxRRkaGLdOZ3M0pJSVFq1at0tNPP61u3bqpfv36uvjiizVu3Dhde+21kqxWp5kzZ6pXr14KCgpSo0aN9N///tdlO/v27dONN96o8PBw1ahRQ3379tXu3btdysyaNUstW7ZUYGCgYmJiNHLkSElSgwYNJEn9+vWTw+FwPp4wYYLatWunN99802XY9QYNGuj555932Xa7du00YcIE52OHw6H//Oc/uvrqqxUcHKzmzZtr9erV2r59u7p27aqQkBB16tRJO3bsKPV7da4IWwAAAEAZyMzMVNWqVW2ZMjMzS13PwufMnz9f2dnZJy332GOPqX///tq4caMGDRqkm266Sb/88oskKTc3V/Hx8apWrZpWrVql77//XlWrVlXPnj2Vk5MjSZo5c6ZGjBihO++8U5s3b9YXX3yhJk2aSJLWrVsnSZo9e7YOHDjgfCxJ27dv1yeffKJPP/1UCQkJZ/Q3mDRpkm699VYlJCSoWbNmGjhwoP71r39p3LhxWr9+vYwxzsDnCXQjBAAAALyIn5+f5syZozvuuEOvvfaa2rdvry5duuimm25SmzZtnOVuuOEG/fOf/5RkhZjFixfrpZde0quvvqoPP/xQBQUFevPNN533pJo9e7bCw8O1fPly9ejRQ08++aTuv/9+3Xvvvc5tXnTRRZKkyMhISVJ4eLiio6Nd6peTk6N33nnHWeZM3HbbbbrxxhslSWPHjlVcXJwee+wxxcfHS5Luvfde3VY4WJgHELYAAKjkfvnlF9WoUUNRUVF2VwWo1IKDg5Wenm7ba5+J/v37q0+fPlq1apV+/PFHLVy4UNOmTdObb76poUOHSpLi4uJcnhMXF+dsadq4caO2b9+uatWquZTJysrSjh07lJycrP379+vKK688432pX7/+WQUtSS5hsfAzr3Xr1i7LsrKylJaWptDQ0LN6jTNB2AK8SWCg9NFHx+cBVHp79uxRixYtJOmMrukAcOYcDodCQkLsrkapValSRVdddZWuuuoqPfbYY/rnP/+p8ePHO8PWqaSnp6tDhw6aW3j/ziIiIyPl43P2VyuV9B76+PgU+wzLzc0tVs6/yG1tClvcSlpWUFBw1vU7E1yzBXgTPz/r/lo33GDNA6j0NmzYYHcVAFQQLVq0UEZGhvPxjz/+6LL+xx9/VPPmzSVJ7du31++//65atWqpSZMmLlNYWJiqVaumBg0aaOnSpSd9PX9/f+Xn55eqbpGRkTpw4IDzcVpamnbt2nUmu2cLwhYAAADgRQ4dOqQrrrhC7733njZt2qRdu3bp448/1rRp09S3b19nuY8//lizZs3Sb7/9pvHjx2vt2rXOwSUGDRqkiIgI9e3bV6tWrdKuXbu0fPly/fvf/9Yff/whyRpZ8Nlnn9WLL76o33//XT/99JNeeukl5/YLw1hSUpKOHDlyyjpfccUVevfdd7Vq1Spt3rxZQ4YMka+vrxvenbLFT9uAN8nLkz77zJrv14/WLQAAvFDVqlXVsWNHPffcc9qxY4dyc3MVGxurO+64Qw8//LCz3MSJE/XBBx/o7rvvVkxMjN5//31nt+Tg4GCtXLlSY8eO1fXXX6+jR4+qTp06uvLKK53XQg0ZMkRZWVl67rnn9MADDygiIkL/+Mc/nNt/9tlnNXr0aL3xxhuqU6dOsWHjixo3bpx27dqlq6++WmFhYZo0aVKFaNlyGDpwn1ZaWprCwsKUmprqkQvpALfJyJCqVrXm09OlCtSvHMDZ+eSTT5wnN3zlA2UnKytLu3btcrkXVGXicDj02Wef6brrrrO7KrY41d/3TLIB3QgBAAAAwA0IWwAAVGK0ZgGAfbhgAwAAAIALfqgpG7RsAQBQiXHCBAD2IWwBAAAAgBsQtgAAAADADbhmC/AmAQHS7NnH5wEAAOA2hC3Am/j7S0OH2l0LAB7ENVsAYB+6EQIAAACAGxC2AG+Slyd9+aU15eXZXRsAHkDLFoBzsWjRIrVp00bBwcFq0qSJXn31VburVKEQtgBvkp0tXX21NWVn210bAABgg6FDh+q6664rtnz58uVyOBxKSUmRJB06dEg33nijOnXqpE2bNmn69Om677779OOPP572NT755BN17dpVYWFhqlq1qtq0aaMnnnhChw8fLrP9OLG+5RFhCwAAAEAxa9euVXp6up599lk1adJE/fr1U7t27fTDDz+c8nmPPPKIBgwYoIsuukgLFy7Uli1b9Oyzz2rjxo169913PVT70jPGKM9NPX4IWwAAAEBZysg4+ZSVVfqyx46VrqybJCcnKzQ0VCEhIc5lBw8elJ/fycfYW7t2rSZPnqxnn31W06dPV6dOndSgQQNdddVV+uSTTzRkyBBn2c8//1zt27dXlSpV1KhRI02cONEl9DgcDr355pvq16+fgoODdd555+mLL76QJO3evVvdunWTJFWvXl0Oh0ND/38QsIKCAk2ZMkUNGzZUUFCQ2rZtq//+97/O7Ra2iC1cuFAdOnRQYGCgvvvuuzJ5z07EaIQAAFRiXLMF2KBq1ZOv693buna6UK1aUmZmyWW7dJGWLz/+uEED6e+/i5dz0//zop8f6enpevrpp/Xnn3/qqquuOulz5s6dq6pVq+ruu+8ucX14eLgkadWqVbr11lv14osv6vLLL9eOHTt05513SpLGjx/vLD9x4kRNmzZN06dP10svvaRBgwZpz549io2N1SeffKL+/fsrMTFRoaGhCgoKkiRNmTJF7733nl577TWdd955WrlypW655RZFRkaqS5cuzm0/9NBDeuaZZ9SoUSNVr179rN+nUyFsAQAAAF5mwYIFqnpCKMzPzy+x7LZt29SyZUtJ0vXXX68mTZqcdLu///67GjVqJH9//1O+/sSJE/XQQw85W7oaNWqkSZMm6cEHH3QJW0OHDtXNN98sSZo8ebJefPFFrV27Vj179lSNGjUkSbVq1XKGuOzsbE2ePFlLlixRXFycc9vfffed/vOf/7iErSeeeOKUwbEsELYAAKjEaNkCbJCefvJ1vr6uj5OTT17W54QrfnbvPusqnahbt26aOXOmy7I1a9bolltuKVa2SZMmWr9+vTZv3qxJkyapX79+WrBgQYnbLe1nzsaNG/X999/rqaeeci7Lz89XVlaWMjMzFRwcLElq06aNc31ISIhCQ0OVfIr3bPv27crMzCwWonJycnTBBRe4LLvwwgtLVddzQdgCAKASI2wBNihyjZNtZU+7qZBiLVR//PFHiWUDAgLUoUMHdejQQU2aNNHll1+u3bt3q0GDBsXKnn/++fruu++Um5t7ytat9PR0TZw4Uddff32xdVWqVHHOn7gNh8OhgoKCU25Xkr788kvVqVPHZV1gYKDL45AyfD9PhrAFeJOAAOnll4/PAwAAnEJBQYGMMXI4HJKOh5+igaiogQMH6sUXX9Srr76qe++9t9j6lJQUhYeHq3379kpMTDxll8TTCfj/c5mi3R9btGihwMBA7d2716XLoF0IW4A38feXRoywuxYAAKCCOHr0qG677TaNHTtWubm5evjhh3XBBRcoOjq6xPIdO3bUgw8+qPvvv19//vmn+vXrp9q1a2v79u167bXXdNlll+nee+/V448/rquvvlr16tXTP/7xD/n4+Gjjxo3asmWLnnzyyVLVrX79+nI4HFqwYIF69+6toKAgVatWTQ888IBGjRqlgoICXXbZZUpNTdX333+v0NBQl9EQPYGh3wEAAACUKCgoSCkpKbrkkkvUrVs3hYSE6OOPPz7lc55++mnNmzdPa9asUXx8vFq2bKnRo0erTZs2zrATHx+vBQsW6JtvvtFFF12kSy65RM8995zq169f6rrVqVPHOdBGVFSURo4cKUmaNGmSHnvsMU2ZMkXNmzdXz5499eWXX6phw4Zn/0acJYehM/dppaWlKSwsTKmpqQoNDbW7OsDZy8+XVq2y5i+/vPhFugAqnblz5zoveOcrHyg7WVlZ2rVrlxo2bHjSLnUV3Zw5c3TfffcpJSXF7qp43Kn+vmeSDehGCHiTrCzp/28AqPT0Mr3QFgAAAK7oRggAQCVGaxYA2IewBQAAAKCYoUOHemUXwrJE2AIAAAAANyBsAQAAAGeJrrqVU1n9XQlbAABUYpwIAu5ReHPfzMxMm2sCd8jJyZEk+Z7jyM2MRggAAACcIV9fX4WHhys5OVmSFBwcLIfDYXOtUBYKCgr0119/KTg4WH5+5xaXCFuAN/H3l6ZNOz4PoNKjZQtwn+joaElyBi5UHj4+PqpXr945B2jCFuBNAgKkMWPsrgUAAJWCw+FQTEyMatWqpdzcXLurgzIUEBAgH59zv+LK1rC1cuVKTZ8+XRs2bNCBAwf02Wef6brrrnOuP1mSnDZtmsb8/wljgwYNtGfPHpf1U6ZM0UMPPeR8vGnTJo0YMULr1q1TZGSk7rnnHj344INlv0MAAADwOr6+vud8bQ8qJ1sHyMjIyFDbtm31yiuvlLj+wIEDLtOsWbPkcDjUv39/l3JPPPGES7l77rnHuS4tLU09evRQ/fr1tWHDBk2fPl0TJkzQ66+/7tZ9A8ql/Hxp3Tprys+3uzYAPIBuhABgH1tbtnr16qVevXqddH1hP9hCn3/+ubp166ZGjRq5LK9WrVqxsoXmzp2rnJwczZo1SwEBAWrZsqUSEhI0Y8YM3Xnnnee+E0BFkpUlXXyxNZ+eLoWE2FsfAACASqzCDP1+8OBBffnllxo2bFixdVOnTlXNmjV1wQUXaPr06crLy3OuW716tTp37qyAgADnsvj4eCUmJurIkSMlvlZ2drbS0tJcJgAAKiJatgDAPhVmgIy3335b1apV0/XXX++y/N///rfat2+vGjVq6IcfftC4ceN04MABzZgxQ5KUlJSkhg0bujwnKirKua569erFXmvKlCmaOHGim/YEAAAAgDeoMGFr1qxZGjRokKpUqeKyfPTo0c75Nm3aKCAgQP/61780ZcoUBQYGntVrjRs3zmW7aWlpio2NPbuKAwBgI1q2AMA+FSJsrVq1SomJifrwww9PW7Zjx47Ky8vT7t271bRpU0VHR+vgwYMuZQofn+w6r8DAwLMOagAAAAAgVZBrtt566y116NBBbdu2PW3ZhIQE+fj4qFatWpKkuLg4rVy50uXeB4sXL1bTpk1L7EIIAAAAAGXB1rCVnp6uhIQEJSQkSJJ27dqlhIQE7d2711kmLS1NH3/8sf75z38We/7q1av1/PPPa+PGjdq5c6fmzp2rUaNG6ZZbbnEGqYEDByogIEDDhg3T1q1b9eGHH+qFF15w6SYIAIA3oEshAHiWrd0I169fr27dujkfFwagIUOGaM6cOZKkDz74QMYY3XzzzcWeHxgYqA8++EATJkxQdna2GjZsqFGjRrkEqbCwMH3zzTcaMWKEOnTooIiICD3++OMM+w7v5O8vjR9/fB5ApUfAAgD7OAyfwqeVlpamsLAwpaamKjQ01O7qAABQarNmzXLeNiU/P18+PhXiCgIAKLfOJBvwiQsAQCXGb6oAYJ8KMRohgDJSUCD98os137y5xC/cgFcheAGAZxG2AG9y7JjUqpU1n54uhYTYWx8AAIBKjJ+1AQDwErRsAYBnEbYAAKjECFgAYB/CFgAAlRhhCwDsQ9gCAMBLELwAwLMIWwAAVGIELACwD2ELAAAvQfACAM9i6HfAm/j7Sw88cHweAAAAbkPYArxJQIA0fbrdtQBgE1q2AMCz6EYIAEAlRsACAPvQsgV4k4ICae9ea75ePcmH31sAb0LwAgDPImwB3uTYMalhQ2s+PV0KCbG3PgDcjoAFAPbhZ20AALwEwQsAPIuwBQAAAABuQNgCAMBL0LIFAJ5F2AIAoBIjYAGAfQhbAAB4CYIXAHgWYQsAgEqMgAUA9mHod8Cb+PlJd999fB5ApUfYAgD7cLYFeJPAQOmVV+yuBQCbELwAwLPoRggAAAAAbkDLFuBNjJH+/tuaj4iQHA576wPAo2jZAgDPImwB3iQzU6pVy5pPT5dCQuytDwAAQCVGN0IAALwELVsA4FmELQAAAABwA8IWAABegpYtAPAswhYAAJWYg4FwAMA2hC0AALwELVsA4FmELQDus2iR1KaNFBwsNWkivfrquW1vwgSpXbvSl3c4pPnzz+01z0bXrtJ993n+dQEAQLlC2AK8iZ+fNGSINfmd5Z0fhg6Vrruu+PLly61wk5JiPT50SLrxRqlTJ2nTJmn6dCuA/Pjjybe9e7e1jYSEktc/8IC0dOnZ1bs05syRwsPdt313ysqy/jatW1t/25L+Rifz5ZdSx45SUJBUvXrx5y5dav0dq1WToqOlsWOlvLzj6ydMsP5uJ05Fby2Qmys98YTUuLFUpYrUtq0Vxt1t0ybp8sut14yNlaZNO3X5Q4eknj2l2rWlwEDrOSNHSmlpx8sUHusnTklJx8s0aFBymREjjpfZsUPq10+KjJRCQ63/LwcPluXeF0PLFgB4FmEL8CaBgVagmDPHmnentWute3k9+6zVqtWvn9Uq9cMPZ7/NqlWlmjXLrIqVSn6+FZb+/W+pe/fSP++TT6TBg6XbbpM2bpS+/14aOPD4+o0bpd69rQDy88/Shx9KX3whPfTQ8TIPPCAdOOA6tWgh3XDD8TKPPir95z/SSy9J27ZJw4dbx8TPP5e+rnPmWK2GpZWWJvXoIdWvL23YYAX+CROk118/+XN8fKS+fa19/O036zWXLLHqe6LERNd9LryHnSStW+e6bvFia3nhe5KRYdXN4ZC+/dZ633NypGuukQoKSr+PAIByjbAFwD2Sk61f64u2bhw8ePYtalLJ3QhnzZJatrTCY0yM1QpxMuPHW2U2bTqz13v3XaulIixMuukm6ejR42UyMqRbb7WCYEyMFS5PlJ1tBZI6daz3o2NHq3VEslqkWraU7rzzePkdO6xWpFmzSldPydruzJnSHXdYrU+lkZcn3XuvFUKGD5fOP98KSTfeeLzMhx9aXUEff9wKzV26WK1Dr7xy/H2oWtV6zcLp4EErUA0bdnw7774rPfywFdwaNZLuusuaL+n9Kitz51oBpvAYuekmK4zOmHHy51SvbtXtwgutkHblldLdd0urVhUvW6uW6377FPlKjYx0XbdggdWq16WLtf77762W3DlzrNbI1q2lt9+W1q+3wpeb0LIFAJ5F2AK8iTFWOMjIsObd/VqF0tOlxx6T/vxTuuqqsnuNmTOtbll33ilt3my1RjRpUnJd7rlHeucd66S5TZvSv8aOHdZ1XwsWWNOKFdLUqcfXjxljLfv8c+mbb6wQ9dNPrtsYOVJavVr64AMr6N1wg9VS9PvvVve2uXOtE+3PP7daqG65xXqfbr/9+DYcDuvEvCz99JP1N/HxkS64wAqLvXpJW7YcL5OdbdWxqKAgKyRu2FDydt980wpul19++u18913Z7EtJVq+WOneWAgKOL4uPt1qkjhwp3Tb275c+/fR4SCqqXTvrPbvqKis8nUxOjvTee9bfs3BkwOxsa75oC3OVKtbfoozfE0YjBAD7ELYAb5KZabVCVK1qzZ+tBQuOb6dw6tWr5LLbtlmtNE8+aXXPKikMna0nn5Tuv99qnTn/fOmii4oPTJGXZ4WXpUutk9gzff2CAivktGplhYfBg49fN5aeLr31lvTMM1YLSGHrRNHrmfbulWbPlj7+2Hp+48ZWK9dll1nLJeuk/cknpX/+06r/nj3SG2+41qNpU6tlrSzt3Gn9O2GC1c1vwQKrZadrV+nwYWtdfLzV9fP9960g+Oef1rVXktU97kRZWVZ4LNqqVbidGTOsgFlQYHWr+/TTkrdRaO9e12Ns+HArLBddNnnyyZ+flCRFRbkuK3xc9Pqqktx8szWwS506Vgvtm28eXxcTI732mtUF85NPrOu6unYtHrILzZ9vXcs4dOjxZZdcYrVGjh1r/V/MyLCOi/z8U78nZ4HWLACwD2ELwJnr1s0axKLoVPRktKgmTayuUbNnW+X69SubOiQnW60OV1556nKjRklr1kgrV1onzmeqQQMrLBaKibFeW7JavXJyrG6BhWrUsIJRoc2brRPo8893DQkrVljPL3T//VaZl1+2ur2deG3ar7+W3XtXqPDaoEcekfr3lzp0sP5ODocVDiXruqLCboaBgVYde/e21vmU8BXy2WdW98IhQ1yXv/CCdN55UrNmVkvTyJHWdWIlbaNQ7dqux9gTT1jd+4ouK+laqrLw3HNWePr8c+vvNHr08XVNm0r/+pf1fnXqZP29OnWynlOSt96yfoyoXfv4sshI6z3+3/+s4yEszApk7duf+j05RwQvAPCsc7h4AoDXCgkp3kL0xx8llw0IsE5KO3SwnnP55da1Kg0anFsdgoJKV+6qq6xWma+/lgYNOvPX8fd3fexwnNkABunpkq+v1eXO19d1XdWqx+eTk60BGXx9rdafnj3PvK5nKibG+rdFi+PLAgOta6r27j2+bPRoK7QeOGC1fO3eLY0bZ5U70ZtvSldfXbxFKTLSauHJyrJG/Ktd2xpko6RtFPLzcz3OatWy/u6lbZ0svH6sqMLHp7uurfBaq2bNrAB9+eVWV9jC9+xEF19ccve/PXusATY+/bT4uh49rCD399/WvoaHW695qvfkLNCNEADsQ8sWAPcpKHC9dqswuJx47c7ZqFbNCmynGwr+2mulefOsLnoffHDur1tU48bWPq1Zc3zZkSNWaCp0wQVWy1ZyshUSik5FT/hvv/14N8SxY6VffinbupakQwcrXCUmHl+Wm2uFqfr1Xcs6HFZACgqywmtsrNUKU9SuXdKyZcW7EBZVpYrVwpiXZ3XB69u3zHanmLg4q0UzN/f4ssWLrZap6tVLv53CcJ2dffIyCQklB7HZs62Q2KfPyZ8bEWEFrW+/tY6Ta68tfd3OEC1bAOBZtGwBcJ+jR62uYmPHWie8Dz9shY/TtSoUPfkv1LJl8WUTJljdyGrVsrppHT1qDVRwzz2u5fr1s0bDGzzYakH4xz/OepdcVK1qBYsxY6xuf7VqWV3yinYDO/98q0Xt1lutkfcuuED66y8rJLZpY52Ev/KKNZjDpk1WiPnyS+s5P/54fHCHZs2kKVNO3ZVw2zarW+Phw9Z7UXi/ssIRHNeuteqxdOnxa5GGD7dGaYyNtQLW9OlW2aLDtk+fbrW0+fhYLTRTp0offVS8pW7WrOODbJxozRrreq927ax/J0ywQsyDD558f/LzrfeqUM+e1lT0eqvCbpklGThQmjjR+huNHWsN/PHCC67d/T77zGql+/VX6/FXX1mtXxddZG1361br73vppcdbY59/XmrY0Doms7Ks1rxvv7UGSCmqoMAKWye7r93s2VLz5lar3+rV1rWHo0a5dkMFAFRohC0A7hMUZF2Hcskl1snmpZee+h5HhW66qfiyffuKLxsyxDrZfe45a3CBiIiTB6l//MM6+R082AoN119/RrtyUtOnW10Fr7nGam27/34pNdW1zOzZxwfz+PNPq56XXGJ1t/v1V+tk/q23rMAjSa++agWxxx6Tnn7aWpaYWHy7J+rd2+q2VuiCC6x/C1szMjOt7RRt6Zk+3frbDB4sHTtmXX/27beuLT8LF0pPPWW17LRta13HdGKgKhxIZOjQ4iFMsv5Ojz5qDcpRtapV13ffPfVNpPfts0LNqYwfbwW3koSFWQFoxAirFS8iwhrCvugw+6mpruE+KMganGTUKGt/Y2OtY6XofcVyco7/LYODrb/VkiXWtYxFLVlidccsOqpkUYmJVtA7fNgKco88Yr2uG9GyBQCe5TB88p5WWlqawsLClJqaqtDQULurA5y9jIzjrQDp6a73wCprc+ZYI+ulpLjvNQCc1syZM3X33XdLkpKTkxUZGWlzjQCgYjuTbEDLFuBNfH2Pt/yU1PoAoFLj91UA8CzCFuBNqlQ5PqQ3AAAA3IrRCAG4x9ChdCEEyoGirVm0bAGAZ9katlauXKlrrrlGtWvXlsPh0Pz5813WDx06VA6Hw2XqecK9Zw4fPqxBgwYpNDRU4eHhGjZsmNLT013KbNq0SZdffrmqVKmi2NhYTZs2zd27BgAAAMDL2Rq2MjIy1LZtW73yyisnLdOzZ08dOHDAOb3//vsu6wcNGqStW7dq8eLFWrBggVauXKk7i4w0lZaWph49eqh+/frasGGDpk+frgkTJuj10oyIBlQ2GRnW/ZIcDmsegFehZQsAPMvWa7Z69eqlXiXdj6WIwMBARZ/knjy//PKLFi1apHXr1unCCy+UJL300kvq3bu3nnnmGdWuXVtz585VTk6OZs2apYCAALVs2VIJCQmaMWOGSygrKjs7W9lFbl6ZlpZ2lnsIAIC9CFgAYJ9yf83W8uXLVatWLTVt2lR33XWXDh065Fy3evVqhYeHO4OWJHXv3l0+Pj5as2aNs0znzp0VUHhjUEnx8fFKTEzUkSNHSnzNKVOmKCwszDnFFt77BgCACozgBQCeVa7DVs+ePfXOO+9o6dKlevrpp7VixQr16tVL+fn5kqSkpCTVqlXL5Tl+fn6qUaOGkpKSnGWioqJcyhQ+LixzonHjxik1NdU57SvpZqoAAAAAcArleuj3m266yTnfunVrtWnTRo0bN9by5ct15ZVXuu11AwMDFRgY6LbtAwDgKYxGCAD2KdctWydq1KiRIiIitH37dklSdHS0kpOTXcrk5eXp8OHDzuu8oqOjdfDgQZcyhY9Pdi0YAAAAAJyrChW2/vjjDx06dEgxMTGSpLi4OKWkpGjDhg3OMt9++60KCgrUsWNHZ5mVK1cqNzfXWWbx4sVq2rSpqlev7tkdAAAAAOA1bA1b6enpSkhIUEJCgiRp165dSkhI0N69e5Wenq4xY8boxx9/1O7du7V06VL17dtXTZo0UXx8vCSpefPm6tmzp+644w6tXbtW33//vUaOHKmbbrpJtWvXliQNHDhQAQEBGjZsmLZu3aoPP/xQL7zwgkaPHm3XbgP28fWVeve2Jl9fu2sDwAPoRggA9rH1mq3169erW7duzseFAWjIkCGaOXOmNm3apLffflspKSmqXbu2evTooUmTJrlcTzV37lyNHDlSV155pXx8fNS/f3+9+OKLzvVhYWH65ptvNGLECHXo0EERERF6/PHHTzrsO1CpVakiffml3bUAAADwCraGra5du57yV7avv/76tNuoUaOG5s2bd8oybdq00apVq864fgAAVCa0bAGAZ1Woa7YAAMCZIWABgH0IW4A3yciQQkKsKSPD7toA8DCCFwB4Vrm+zxYAN8jMtLsGAAAAXoGWLQAAKjFGIwQA+xC2AAAAAMANCFsAAHgJWrYAwLMIWwAAVGIELACwD2ELAAAvQfACAM9iNELAm/j4SF26HJ8HAACA2xC2AG8SFCQtX253LQB4EKMRAoB9+GkbAAAAANyAsAUAAAAAbkDYArxJRoYUGWlNGRl21waAB9CNEADswzVbgLf5+2+7awAAAOAVaNkCAMBL0LIFAJ5F2AIAoBIjYAGAfQhbAAB4CYIXAHgWYQsAAAAA3ICwBQBAJcZohABgH0YjBLyJj4904YXH5wEAAOA2hC3AmwQFSevW2V0LADahZQsAPIuftgEAqMQIWABgH8IWAABeguAFAJ5F2AK8SWam1KCBNWVm2l0bAACASo1rtgBvYoy0Z8/xeQCVHqMRAoB9aNkCAAAAADcgbAEAAACAGxC2AACoxOhGCAD2IWwBAAAAgBsQtgAA8BK0bAGAZzEaIeBNHA6pRYvj8wAqPQIWANiHsAV4k+BgaetWu2sBwCYELwDwLLoRAgAAAIAbELYAAKjEGI0QAOxD2AK8SWam1LKlNWVm2l0bAACASo1rtgBvYoy0bdvxeQBehZYtAPAsWrYAAKjECFgAYB/CFgAAXoLgBQCeRdgCAAAAADcgbAEAUIkxGiEA2IewBQAAAABuwGiEgDdxOKT69Y/PAwAAwG0IW4A3CQ6Wdu+2uxYAPIhuhABgH7oRAgAAAIAbELYAAPAStGwBgGfZGrZWrlypa665RrVr15bD4dD8+fOd63JzczV27Fi1bt1aISEhql27tm699Vbt37/fZRsNGjSQw+FwmaZOnepSZtOmTbr88stVpUoVxcbGatq0aZ7YPaD8OXZMuugiazp2zO7aAPAAAhYA2MfWsJWRkaG2bdvqlVdeKbYuMzNTP/30kx577DH99NNP+vTTT5WYmKhrr722WNknnnhCBw4ccE733HOPc11aWpp69Oih+vXra8OGDZo+fbomTJig119/3a37BpRLBQXS+vXWVFBgd20AeBjBCwA8y9YBMnr16qVevXqVuC4sLEyLFy92Wfbyyy/r4osv1t69e1WvXj3n8mrVqik6OrrE7cydO1c5OTmaNWuWAgIC1LJlSyUkJGjGjBm68847S3xOdna2srOznY/T0tLOdNcAAAAAeLkKdc1WamqqHA6HwsPDXZZPnTpVNWvW1AUXXKDp06crLy/PuW716tXq3LmzAgICnMvi4+OVmJioI0eOlPg6U6ZMUVhYmHOKjY11y/4AAOBujEYIAPapMGErKytLY8eO1c0336zQ0FDn8n//+9/64IMPtGzZMv3rX//S5MmT9eCDDzrXJyUlKSoqymVbhY+TkpJKfK1x48YpNTXVOe3bt88NewQAAACgMqsQ99nKzc3VjTfeKGOMZs6c6bJu9OjRzvk2bdooICBA//rXvzRlyhQFBgae1esFBgae9XMBACivaNkCAM8q9y1bhUFrz549Wrx4sUurVkk6duyovLw87f7/G7dGR0fr4MGDLmUKH5/sOi8AACoLAhYA2Kdch63CoPX7779ryZIlqlmz5mmfk5CQIB8fH9WqVUuSFBcXp5UrVyo3N9dZZvHixWratKmqV6/utroD5VZEhDUB8DoELwDwLFu7Eaanp2v79u3Ox7t27VJCQoJq1KihmJgY/eMf/9BPP/2kBQsWKD8/33mNVY0aNRQQEKDVq1drzZo16tatm6pVq6bVq1dr1KhRuuWWW5xBauDAgZo4caKGDRumsWPHasuWLXrhhRf03HPP2bLPgK1CQqS//rK7FgAAAF7B1rC1fv16devWzfm48PqrIUOGaMKECfriiy8kSe3atXN53rJly9S1a1cFBgbqgw8+0IQJE5Sdna2GDRtq1KhRLtdxhYWF6ZtvvtGIESPUoUMHRURE6PHHHz/psO8AAFQmjEYIAPaxNWx17dr1lB/8p/tSaN++vX788cfTvk6bNm20atWqM64fAAAAAJytcn3NFoAyduyY1LWrNR07ZndtAAAAKrUKMfQ7gDJSUCCtWHF8HkClRzdCALAPLVsAAAAA4AaELQAAvAQtWwDgWYQtAAAqMQIWANiHsAUAgJcgeAGAZxG2AAAAAMANGI0Q8DbBwXbXAIAHMRohANiHsAV4k5AQKSPD7loAAAB4BboRAgDgJWjZAgDPImwBAFCJEbAAwD6ELcCbZGVJffpYU1aW3bUB4GEELwDwLK7ZArxJfr701VfH5wEAAOA2tGwBAFCJMRohANiHsAUAAAAAbkDYAgDAS9CyBQCeRdgCAKASI2ABgH0IWwAAAADgBoQtAAC8BK1cAOBZDP0OeJOQEImTLcCrELAAwD60bAEA4CUIXgDgWYQtAAAAAHADwhbgTbKypBtusKasLLtrA8ADuKkxANiHsAV4k/x86b//tab8fLtrAwAAUKkRtgAA8BK0bAGAZxG2AACoxAhYAGAfwhYAAF6C4AUAnkXYAgAAAAA3IGwBAFCJMRohANiHsAUAAAAAbuBndwUAeFBwsJSefnwegFehZQsAPIuwBXgTh0MKCbG7FgA8iIAFAPahGyEAAAAAuAFhC/Am2dnS0KHWlJ1td20AeBitXADgWYQtwJvk5Ulvv21NeXl21waABxCwAMA+hC0AALwEwQsAPIuwBQAAAABuQNgCAKAS46bGAGAfwhYAAAAAuAFhCwAAL0HLFgB4FmELAIBKjIAFAPbxs7sCADwoOFhKTj4+D8CrELwAwLPOqmUrKyvrpOsOHDhw1pUB4GYOhxQZaU0Oh921AQAAqNTOKmy1b99eCQkJxZZ/8sknatOmzbnWCQAAlBFGIwQA+5xV2OratasuueQSPf3005KkjIwMDR06VIMHD9bDDz9cphUEUIays6URI6wpO9vu2gAAAFRqZxW2Xn31VX3yySd6/vnndfnll6tt27ZKSEjQ2rVrNWrUqFJvZ+XKlbrmmmtUu3ZtORwOzZ8/32W9MUaPP/64YmJiFBQUpO7du+v33393KXP48GENGjRIoaGhCg8P17Bhw5Senu5SZtOmTbr88stVpUoVxcbGatq0aWez20DFl5cnvfqqNeXl2V0bAB5GyxYAeNZZj0bYq1cvXX/99fr++++1d+9ePf3002rVqtUZbSMjI0Nt27bVK6+8UuL6adOm6cUXX9Rrr72mNWvWKCQkRPHx8S7XjA0aNEhbt27V4sWLtWDBAq1cuVJ33nmnc31aWpp69Oih+vXra8OGDZo+fbomTJig119//ex2HACACoSABQD2OavRCHfs2KGBAwcqKSlJX3/9tVasWKFrr71W9957r5566in5+/uXaju9evVSr169SlxnjNHzzz+vRx99VH379pUkvfPOO4qKitL8+fN100036ZdfftGiRYu0bt06XXjhhZKkl156Sb1799Yzzzyj2rVra+7cucrJydGsWbMUEBCgli1bKiEhQTNmzHAJZUVlZ2cru0gXq7S0tDN5ewAAAADg7Fq22rVrp4YNG2rjxo266qqr9OSTT2rZsmX69NNPdfHFF5dJxXbt2qWkpCR1797duSwsLEwdO3bU6tWrJUmrV69WeHi4M2hJUvfu3eXj46M1a9Y4y3Tu3FkBAQHOMvHx8UpMTNSRI0dKfO0pU6YoLCzMOcXGxpbJPgEAYCdauQDAs876mq0PPvhA4eHhzmWdOnXSzz//rPbt25dJxZKSkiRJUVFRLsujoqKc65KSklSrVi2X9X5+fqpRo4ZLmZK2UfQ1TjRu3DilpqY6p3379p37DgEAYAMCFgDY56y6EQ4ePLjE5dWqVdNbb711ThUqDwIDAxUYGGh3NQAAKFMELwDwrLMKW++8885J1zkcjpOGsTMRHR0tSTp48KBiYmKcyw8ePKh27do5yyQnJ7s8Ly8vT4cPH3Y+Pzo6WgcPHnQpU/i4sAwAAAAAlLWzClv33nuvy+Pc3FxlZmYqICBAwcHBZRK2GjZsqOjoaC1dutQZrtLS0rRmzRrdddddkqS4uDilpKRow4YN6tChgyTp22+/VUFBgTp27Ogs88gjjyg3N9c5cMfixYvVtGlTVa9e/ZzrCVQoQUHSrl3H5wFUetzUGADsc1bXbB05csRlSk9PV2Jioi677DK9//77pd5Oenq6EhISlJCQIMkaFCMhIUF79+6Vw+HQfffdpyeffFJffPGFNm/erFtvvVW1a9fWddddJ0lq3ry5evbsqTvuuENr167V999/r5EjR+qmm25S7dq1JUkDBw5UQECAhg0bpq1bt+rDDz/UCy+8oNGjR5/NrgMVm4+P1KCBNfmc9Z0fAAAAUApn1bJVkvPOO09Tp07VLbfcol9//bVUz1m/fr26devmfFwYgIYMGaI5c+bowQcfVEZGhu68806lpKTosssu06JFi1SlShXnc+bOnauRI0fqyiuvlI+Pj/r3768XX3zRuT4sLEzffPONRowYoQ4dOigiIkKPP/74SYd9BwCgsqJlCwA8q8zClmSNBLh///5Sl+/atespP/gdDoeeeOIJPfHEEyctU6NGDc2bN++Ur9OmTRutWrWq1PUCKq2cHOmRR6z5p56SitwSAUDlRMACAPucVdj64osvXB4bY3TgwAG9/PLLuvTSS8ukYgDcIDdXeuYZa37CBMIW4GUIXgDgWWcVtgqvmSrkcDgUGRmpK664Qs8++2xZ1AsAAAAAKrSzClsFBQVlXQ8AAOAGjEYIAPZhODIAAAAAcINSt2ydyVDpM2bMOKvKAAAA96FlCwA8q9Rh6+effy5VOYfDcdaVAQAAZYuABQD2KXXYWrZsmXbu3KkGDRrIh5uhAgAAAMApnVFqOu+88/T33387Hw8YMEAHDx4s80oBcJOgIGnLFmsKCrK7NgA8jFYuAPCsMwpbJ35If/XVV8rIyCjTCgFwIx8fqWVLa6KFGvAKBCwAsA9nWwAAeAmCFwB41hndZ8vhcBQbAIMBMYAKJCdHmjzZmn/4YSkgwN76AAAAVGJnFLaMMRo6dKgCAwMlSVlZWRo+fLhCQkJcyn366adlV0MAZSc3V5o40ZofM4awBXiBE29qnJiYqKioKIWHh9tXKQDwEmcUtoYMGeLy+JZbbinTygAAAPfZtm2bBgwYoGrVqiktLc3u6gBApXdGYWv27NnuqgcAAHCzZcuWSZKOHj1qc00AwDswQAYAAJVY0W6EYWFhNtYEALwPYQsAAC9RNGwdO3bMxpoAgHcgbAEA4CWqVq3qnKcrIQC4H2ELAIBK7GT31uKeWwDgfmc0QAaACq5KFWnt2uPzALxKQUGBc56wBQDuR9gCvImvr3TRRXbXAoBNCFgA4Fl0IwQAoBIrGrBo2QIAz6JlC/AmOTnSCy9Y8/feKwUE2FsfAB5VNGARtgDA/QhbgDfJzZUefNCav/tuwhbgZQhbAOBZdCMEAKASO1nAImwBgPsRtgAA8BIELADwLMIWAABeggEyAMCzCFsAAFRidCMEAPsQtgAAqMQY+h0A7EPYAgDASxCwAMCzGPod8CZVqkjLlh2fB1Dp0bIFAPYhbAHexNdX6trV7loAsAnXbAGAZ9GNEACASoyWLQCwDy1bgDfJzZVef92av/NOyd/f3voA8ChatgDAswhbgDfJyZFGjrTmhw4lbAFegIAFAPahGyEAAF6C4AUAnkXYAgCgEuOaLQCwD2ELAAAvQcsWAHgWYQsAgEqMli0AsA9hCwAAL0HAAgDPImwBAFCJnazrIMELANyPod8BbxIYKC1YcHwegFehGyEAeBZhC/Amfn5Snz521wKAB9GyBQD2oRshAABegoAFAJ5FyxbgTXJzpblzrflBgyR/f3vrA8DtaNkCAPsQtgBvkpMj3XabNX/DDYQtwMtwzRYAeFa570bYoEEDORyOYtOIESMkSV27di22bvjw4S7b2Lt3r/r06aPg4GDVqlVLY8aMUV5enh27AwCAR9GyBQD2KfctW+vWrVN+fr7z8ZYtW3TVVVfphhtucC6744479MQTTzgfBwcHO+fz8/PVp08fRUdH64cfftCBAwd06623yt/fX5MnT/bMTgAAUA7QsgUAnlXuw1ZkZKTL46lTp6px48bq0qWLc1lwcLCio6NLfP4333yjbdu2acmSJYqKilK7du00adIkjR07VhMmTFBAQECx52RnZys7O9v5OC0trYz2BgAAz6I1CwDsU+67ERaVk5Oj9957T7fffrscDodz+dy5cxUREaFWrVpp3LhxyszMdK5bvXq1WrduraioKOey+Ph4paWlaevWrSW+zpQpUxQWFuacYmNj3bdTAAC4UdGARcsWAHhWuW/ZKmr+/PlKSUnR0KFDncsGDhyo+vXrq3bt2tq0aZPGjh2rxMREffrpp5KkpKQkl6Alyfk4KSmpxNcZN26cRo8e7XyclpZG4AIAVHi0cgGAZ1WosPXWW2+pV69eql27tnPZnXfe6Zxv3bq1YmJidOWVV2rHjh1q3LjxWb1OYGCgAgMDz7m+AADYjQEyAMA+FSZs7dmzR0uWLHG2WJ1Mx44dJUnbt29X48aNFR0drbVr17qUOXjwoCSd9DovoNIKDJQ++uj4PACvQsACAM+qMNdszZ49W7Vq1VKfPn1OWS4hIUGSFBMTI0mKi4vT5s2blZyc7CyzePFihYaGqkWLFm6rL1Au+flZ99e64QZrHkClxzVbAGCfCnG2VVBQoNmzZ2vIkCHyK3KCuGPHDs2bN0+9e/dWzZo1tWnTJo0aNUqdO3dWmzZtJEk9evRQixYtNHjwYE2bNk1JSUl69NFHNWLECLoKAgC8Ct0IAcCzKkTYWrJkifbu3avbb7/dZXlAQICWLFmi559/XhkZGYqNjVX//v316KOPOsv4+vpqwYIFuuuuuxQXF6eQkBANGTLE5b5cgNfIy5M++8ya79eP1i3AC3DNFgDYp0KcafXo0aPEL4XY2FitWLHitM+vX7++vvrqK3dUDahYsrOlG2+05tPTCVuAl6EbIQB4VoW5ZgsAAJw5WrMAwD6ELQAAvAQtWwDgWYQtAAAqMa7ZAgD7ELYAAPAShC0A8CzCFgAAlRgtWwBgH8IWAABeoug1WwAA92PcZ8CbBARIs2cfnwdQ6dGyBQD2IWwB3sTfXxo61O5aALAJoxECgGfRjRAAgEqMli0AsA8tW4A3ycuTvv7amo+Pl/z4CAC8CddsAYBncaYFeJPsbOnqq6359HTCFuAFaNkCAPvQjRAAAC9B2AIAzyJsAQBQidGyBQD2IWwBAFCJEbYAwD6ELQAAvAQDZACAZxG2AACoxGjZAgD7ELYAAPAS3NQYADyLcZ8BbxIQIL388vF5AJUeLVsAYB/CFuBN/P2lESPsrgUAm3DNFgB4Ft0IAQCoxGjZAgD70LIFeJP8fGnVKmv+8sslX1976wPAowhbAOBZhC3Am2RlSd26WfPp6VJIiL31AeB2tGwBgH3oRggAgJdgNEIA8CzCFgAAlRitWQBgH8IWAABegpYtAPAswhYAAJUY12wBgH0IWwAAeAlatgDAswhbAABUYrRsAYB9GPod8Cb+/tK0acfnAXgVAhYAeBZhC/AmAQHSmDF21wKAB9GyBQD2oRshAABegrAFAJ5FyxbgTfLzpZ9+subbt5d8fe2tDwC3KxqqGCADADyLsAV4k6ws6eKLrfn0dCkkxN76APAoAhYAeBbdCAEAqMRo2QIA+xC2AADwElyzBQCeRdgCAKASo2ULAOxD2AIAoBJj6HcAsA9hCwAAL1G0ZQsA4H6ELQAAKrGTtWDRsgUA7sfQ74A38feXxo8/Pg/AaxG2AMD9CFuANwkIkCZMsLsWADyIli0AsA/dCAEA8EKELQBwP1q2AG9SUCD98os137y55MPvLUBlR6gCAPsQtgBvcuyY1KqVNZ+eLoWE2FsfALYhhAGA+5Xrn7UnTJggh8PhMjVr1sy5PisrSyNGjFDNmjVVtWpV9e/fXwcPHnTZxt69e9WnTx8FBwerVq1aGjNmjPLy8jy9KwAA2IJrtgDAPuW+Zatly5ZasmSJ87Gf3/Eqjxo1Sl9++aU+/vhjhYWFaeTIkbr++uv1/fffS5Ly8/PVp08fRUdH64cfftCBAwd06623yt/fX5MnT/b4vgAAUF4QtgDA/cp92PLz81N0dHSx5ampqXrrrbc0b948XXHFFZKk2bNnq3nz5vrxxx91ySWX6JtvvtG2bdu0ZMkSRUVFqV27dpo0aZLGjh2rCRMmKCAgwNO7AwCARxGqAMA+5boboST9/vvvql27tho1aqRBgwZp7969kqQNGzYoNzdX3bt3d5Zt1qyZ6tWrp9WrV0uSVq9erdatWysqKspZJj4+Xmlpadq6detJXzM7O1tpaWkuEwAAlQkhDADcr1yHrY4dO2rOnDlatGiRZs6cqV27dunyyy/X0aNHlZSUpICAAIWHh7s8JyoqSklJSZKkpKQkl6BVuL5w3clMmTJFYWFhzik2NrZsdwwAAA/hmi0AsE+57kbYq1cv53ybNm3UsWNH1a9fXx999JGCgoLc9rrjxo3T6NGjnY/T0tIIXACASoWwBQDuV67D1onCw8N1/vnna/v27brqqquUk5OjlJQUl9atgwcPOq/xio6O1tq1a122UThaYUnXgRUKDAxUYGBg2e8AYDd/f+mBB47PA6j0aNkCAPuU626EJ0pPT9eOHTsUExOjDh06yN/fX0uXLnWuT0xM1N69exUXFydJiouL0+bNm5WcnOwss3jxYoWGhqpFixYerz9gu4AAafp0a2KAGAAAALcq1y1bDzzwgK655hrVr19f+/fv1/jx4+Xr66ubb75ZYWFhGjZsmEaPHq0aNWooNDRU99xzj+Li4nTJJZdIknr06KEWLVpo8ODBmjZtmpKSkvToo49qxIgRtFwBALwCLVsAYJ9yHbb++OMP3XzzzTp06JAiIyN12WWX6ccff1RkZKQk6bnnnpOPj4/69++v7OxsxcfH69VXX3U+39fXVwsWLNBdd92luLg4hYSEaMiQIXriiSfs2iXAXgUF0v+P6Kl69SSfCtW4DaAMEbYAwP0chk/b00pLS1NYWJhSU1MVGhpqd3WAs5eRIVWtas2np0shIfbWB4DbXXLJJVqzZk2x5bNnz9bQoUM9XyEAqODOJBvwszYAAAAAuAFhCwCASoxrtgDAPoQtAAAqMcIWANiHsAUAgBcibAGA+xG2AACoxGjZAgD7ELYAAAAAwA3K9X22AJQxPz/p7ruPzwOo9GjZAgD7cLYFeJPAQOmVV+yuBYBygLAFAO5HN0IAACoxWrYAwD60bAHexBjp77+t+YgIyeGwtz4AbEPYAgD3I2wB3iQzU6pVy5pPT5dCQuytDwC3I1QBgH3oRgh4qYKCArurAMBGhDAAcD/CFuCl5s6da3cVAHgA12wBgH0IW4CXWr58ud1VAGAjwhYAuB9hC/BS+fn5dlcBgAcQqgDAPoQtwEsRtgDvRggDAPcjbAFeKi8vz+4qAPAArtkCAPsw9DvgTfz8NOf/Z3MYjRDwaoQtAHA/whbgTQIDddv/z/bmRAvwCrRsAYB96EYIeCm6EQIAALgXLVuANzFGwf8/m0/YArwCLVsAYB9atgBvkpmpDEkZkvxycuyuDQA32rJli2JjY7Vly5YS1xO2AMD9CFtlLDMzUx988IFSUlLsrgpwSg6Hw+4qAHCjYcOG6Y8//jjpesIWALgfYauM3X333br55pt18803210VAIAXyzlN6zVhCwDcj7BVxt5++21J0qJFi2yuCQDAm/n48BUPAHbjkxgAgErodGGLli0AcD/CFgAAlRBhCwDsR9gCAKASOt0gOIQtAHA/7rPlJv7+/nZXASjO11cf//9svq0VAeBuXLMFAPYjbAHepEoV3fj/s119fct00ykpKQoICFBwcPDpCwNwO7oRAoD9CFtuwpcYvElWVpYaNWokh8Ohv//+m3t4AeUAYQsA7EcfA8BLleWJ1v79+3XkyBEdPnxY+/btK7PtAjh7hC0AsB9hy034EkO5lJEhI8lIqpJfdldt+Rbpknjs2LEy2y6As+eJATIyMzP1zDPPKDEx8Zy3BQCVEWHLTQhbKO/KsqtfQUFBifMA7OOJATKefvppjRkzRs2aNXP7awFARUTYchOuWUF5V5Y/CBC2gPLHE90IN2zYcM7bAIDKjLDlJoQteBPCFlD+eKIbYfXq1c95GwBQmRG23ISwBW9C2AIqnrIIW2FhYWVQEwCovAhbbsLNJOFNCFuAd/It4/v1AUBlQyJwE1q2UN6V5TGaX2RkQwaHASqGsvi/yv93ADg1bmrsJoQtlEu+vvry/2fLbuB3WraA8sgT12wVlZCQoICAALVo0aJMtwsAFRlhy00IWyiXqlTR1f8/26UMu7oStoCKp6zD1gUXXCDJaummKz0AWPg0dBPCFrwJYQuoeNzVBTC/DG+YDgAVHWHLTfhVD96EsAWgEGELAI4jEbgJLVsolzIylC4pXVKVMjwhImwBFQ8tWwDgflyz5SaELZRXIW7YJmELqHgIWwDgfuW6ZWvKlCm66KKLVK1aNdWqVUvXXXedEhMTXcp07dpVDofDZRo+fLhLmb1796pPnz4KDg5WrVq1NGbMGOXl5bm17oQteBPCFlDxuGvod3d/vwJARVKuW7ZWrFihESNG6KKLLlJeXp4efvhh9ejRQ9u2bVNIyPHf5++44w498cQTzsfBwcHO+fz8fPXp00fR0dH64YcfdODAAd16663y9/fX5MmT3VZ3wha8SdGAxX13gIqBli0AcL9yHbYWLVrk8njOnDmqVauWNmzYoM6dOzuXBwcHKzo6usRtfPPNN9q2bZuWLFmiqKgotWvXTpMmTdLYsWM1YcIEBQQEuKXuhC14E1q2ABQibAHAceW6G+GJUlNTJUk1atRwWT537lxFRESoVatWGjdunDIzM53rVq9erdatWysqKsq5LD4+Xmlpadq6dWuJr5Odna20tDSX6UwxGiG8SdGTK8IWUD54+qbGhQhbAHBcuW7ZKqqgoED33XefLr30UrVq1cq5fODAgapfv75q166tTZs2aezYsUpMTNSnn34qSUpKSnIJWpKcj5OSkkp8rSlTpmjixInnVF9atuBNaNkCKh53hS2u2QKA4ypM2BoxYoS2bNmi7777zmX5nXfe6Zxv3bq1YmJidOWVV2rHjh1q3LjxWb3WuHHjNHr0aOfjtLQ0xcbGntE2CFsol3x8tPz/Z8syEhG2gPLndGGqLMJWSf/fadkCgOMqRF+3kSNHasGCBVq2bJnq1q17yrIdO3aUJG3fvl2SFB0drYMHD7qUKXx8suu8AgMDFRoa6jKdKcIWyqWgIHWT1E1Sjq9vmW2WsAWUP54YrKakYEXYAoDjynXYMsZo5MiR+uyzz/Ttt9+qYcOGp31OQkKCJCkmJkaSFBcXp82bNys5OdlZZvHixQoNDVWLFi3KvL6FCFvwJoQtoPw5Xegpi/+rhC0AOLVy3Y1wxIgRmjdvnj7//HNVq1bNeY1VWFiYgoKCtGPHDs2bN0+9e/dWzZo1tWnTJo0aNUqdO3dWmzZtJEk9evRQixYtNHjwYE2bNk1JSUl69NFHNWLECAUGBpZpfYuGLQbIgDchbAHlj11hi2u2AOC4cp0IZs6cqdTUVHXt2lUxMTHO6cMPP5QkBQQEaMmSJerRo4eaNWum+++/X/3799f//vc/5zZ8fX21YMEC+fr6Ki4uTrfccotuvfVWl/tylZWiXzq0bKFcyshQsqRkSVXK8Ndn7rMFlD+nC1tl0QJFyxYAnFq5btk63UlbbGysVqxYcdrt1K9fX1999VVZVeukip5wErZQHhljFOmG7dKyBZQ/p/u/SDdCAHC/ct2yVdHQsoXyzl2tToQtoPzhmi0AsB9hqwzRsoXyjrAFeA+6EQKA/QhbZYiWLZR37gpbRY99whZQPjBABgDYj7BVhop+cTEaIcojdwUhWraA8oeWLQCwH4mgDNGyhfKOboSA92CADACwX7kejbCi4SQT5V2BpHVF5stsu4QtoNxhgAwAsB8tW2Wo6BcM9xpyv/z8fH300Ufat2+f3VWpMEyVKrpY0sWScnx9y2y73GcLKH/s6kbINVsAcBwtW2WIE07Pmj17tu644w75+fkpNzfX7upUCEWPy7I8RmnZAsofWrYAwH60bJUhRmTzrMWLF0viV9QzwQAZgPc43WcjA2QAgPsRtsoQ3Qg9y7cMu8F5C5ORoV2SdkmqUoahiLAFlD+na/GnZQsA3I9uhGWIboSe5efH4XumTEGBGvz/fFmOl0nYAsqf07VscZ8tAHA/WrbKEN0IPYuWrTNHN0LAe5yuZYtuhADgfoStMkTLlmdx4+gzxwAZgPegGyEA2I+z1TJEyxbKu3MNWJ9//rlWrVpVbDnHPlD+0LIFAPbjopcyRMuWZ/Een7lzCUI7duzQddddJ6n4e8+xD5Q/XLMFAPajZasMMRqhZ/Een7lz6Ua4bdu2kz6XboRA+WKMse2mxrRsAcBxtGyVIU44PYuwdeYKjNHWIvNnIj093Tmfk5OjwMDA49vl2AfKldLc6J1rtgDA/WjZKkO0bHkW7/GZM0FBaiWplaSsMxxgxOE4Plh8dna2y7qiJ3aELcB+penKVxb/V3NycootI2wBwHGErTLEdSuexXt85oq+Z2d6olV09McTw1bREy7CFmC/0rRslUUoOvGzQOKaLQAoirBVhhiRDeXdubS+Fj2BysrKcllH2ALKl6JhKygoSM8++6zzcWhoqKSy+b964meBRMsWABRF2CpDtGx51rm00nir/KNHtUXSFkn+pfjlu6iiJ2+0bAHF7d27t9wc/5mZmZKkwMBApaen67bbbnOui4yMlFS61q/TKalli7AFAMcRtsoQ12x5VtH3mC/30snLzVVLSS0l6QyPUcIWcHJLlixR/fr1NXz4cLurIun4gDZVq1aVj4+Pqlat6lzXsGFDSSUHpTNFN0IAODXCVhnyhhHZdu3apZ9++snuakhyDVt8uZdO0ffpTH8QKBq26EYIuHrttdckSW+88YbNNbEUDVuS5O/vr/Hjx6tNmzb65z//Kencw1ZBQUGJrWNlEeIAoLIgbJWhczmRrSgaNWqkDh06aNeuXXZXxQVhq3TOZdTAU7VsFX1cFl2TgIqmvP3IUBi2QkJCnMsmTJigjRs3qm7dupJKvt7qTJQ0EmFZbBcAKhPCVhkq7CMvVd6wVeiHH36wuwp0IzwL59LV9VRhKyMjwzl/shMwoDIrb9fspqamSpKqVatWbF2VKlUknXsL1MmeT8sWABxH2CpDRU84K+PJf9ETiGPHjtlYk+Jo2Sqdou/Tmf4SX/RvfuIv14UndhJhy5utWrVKAwYM0NatW09fuJIp+v+p6HeBXfbs2SNJqlevXrF1hTckP9cWqLS0tBKX07IFAMcRtspQ0S/Y3NzccvHrZlkqbzeuLRpoK2O4dYdz6epatOX2xF+uU1JSnPOELe/10EMP6aOPPtKtt95qd1U8ruj/iZOFEE/auXOnJKvr94kKw9a5tkD99ddfJS73lrCVmJioVatW2V0NAOUcYasMFQ1bxpgyDQDlIdwU/WIuD/UpGhxo2SqdvPx87Za0W1LOGV5bdaqw9eeffzrnCVveq7B7cXkZRMeTirbulvewVdiN8ExCUX5+vm644QaNGjXKuezgwYMllvWWboRXX321OnfurMWLF9tdFQDlGGGrDBU9GZXK7qTz3//+t3x9ffXRRx8pLS1N8+fPt+WXw6JfoOWhJalowCoP9akIcvz81FBSQ0lpZxhQi/6YUPT4y83NtSVsGWO0cOFCHT582COvh9Lz8fG+r5ZDhw45508Vtt5++20NHz7c7YHkVGErPDxckvV/9+jRo6Xa3oYNG/Tf//5Xzz//vPP/XGJioiSpVatWCg4Odpb1hpatnJwcbd++XZI17L+7rVmzRu3atdNVV12lwYMH68CBA25/TeBMHDp0qFz8EF8eed83ohud2E+/rE46X3rpJUnSgAEDNHjwYPXr108PPvhgqZ5rjNH69evL5BqroicH5eGXS1q2zlzRv9uZHp9HjhwpcTs7duxw+YD11LHx/PPPq3fv3ho2bJhHXg+l5+vra3cVPCorK8tlhNaTha2CggINHTpU//nPfzRr1iy31SctLe2UYatatWoKDQ2VJP3xxx+l2mbRcuvWrZMkfffdd5KkG2+8UTt37tSbb74pyfWzoqLJysrS9OnTtX79+lOWKxp2PPED0y233KKNGzdqyZIleu+99/T6668XK/PGG2+oX79+J21xBNzlf//7nyIiIvT888/bXZVSy8/P1yuvvOL80cidCFtlyB1hq2if+OjoaH3xxReSjgew05k4caIuuugiTZkypcT1CQkJ2r9/v8uylJQU5eXl6a233tJll12mvXv3SnL9tbKk8GaM0Y4dO0psZVq/fn2Zt0CUdcuWMabS/yJb9FfsMw1Fn3zySYnP/eqrr1zKFQ457W7Tp0+XJM2fP99tr1HZrrs8mYKCAiUmJpZqf0vzy6Wfn19ZVKvCSExMdPkMKtrKVdTvv//unN+8ebPb6vPCCy8oOztbDRs2dN7A+ESxsbGSpH379pVqm4XBSrK+NzZs2KD//e9/kqT4+HhFRUWpSZMmZ7RNTzHGlPr7eOrUqXrwwQfVs2fPUx7rRb83k5KSzrmOJSn8Pj106JCzFa3QmjVrXB6np6frzjvv1Pz58/XQQw+d1evt37+/2GdASkqKy7KSbu2RlZXlkRaNw4cPu1wf7I1Od4lKVlaWLT+G33DDDZKk+++//5y3lZubq0GDBqlv377KzMzU0aNHtXDhwjI7xgq3M3PmTI0cOVLNmzd3rjt27JjefPNN9ezZUz///HOxH46ysrK0detW7du3Tz///HOpX9O7vhHdrCzD1vLly/WPf/zD+eujVPwD/eeff1a3bt0UHR2tZcuWKSYmRsYY3XPPPZKkl19+WRMnTpQkTZo0ScOHD1edOnUkSddff73at2+vRx99VFWrVtUjjzyicePGObfdtWtXLV++XJI0a9YsTZgwQZs2bXKuLxq2fvvtN7Vv3965/w0aNNBll12m++67T4GBgWrdunWJ+7hv3z4lJydr1qxZatasmXPZtm3bNGXKFLVq1UqbN29WjRo1nPUuquiHfl5enmbMmKGsrCw98MADOnr0qGrUqKHs7Gzl5eUpMDBQmZmZCgsLc9lGTk6OsrOzVa1aNfXp00cLFy7Uzz//rOjoaNWqVUt5eXmaPn26NmzYoAEDBqhnz57Ov4nD4ZAkbdu2Tb/++qvOP/98xcbG6u+//1ZoaKgiIyP14Ycfatq0abrmmmvUs2dP5efna+3atcrJyVFwcLD++9//6ptvvpGvr69WrFihn376Sf/617+cr3H06FH98MMPuuKKK+Tv76/8/Hzt3btXDRs2lDFG77//vi655BLnr9c//fSTateurejo6BLf88xDh7T2/+evOYMP5Hnz5rk8Lgylxhjnh2vnzp21cuXKYicB7lK021Jubq78/f1P+5z8/HwdO3ZMO3fuVKtWrZSWlqbDhw8rPT1dLVq0kJ+fnzIyMhQYGKgFCxaoX79+6t27tyZOnKjGjRtr//798vX1dR6ve/bsUWRkpHJzc5Wfn69HH31UderUUf369XXLLbe4vO6qVat06aWXyt/fXwUFBZowYYKefPJJLV68WN9//73S0tLUrl07ORwOBQUFqXXr1srNzdXs2bN1xx13KCIiQgkJCSooKFBycrLS0tLUqlUr7d+/X5deeqnWr1+vGTNmaNy4ccrKytK3336rhx56SO+8845WrVqlPn366IorrtBtt92mnTt3Kjk5udj7Ex8fr+7du+uWW25RVFSUnn76aT366KPFvtxr1qypCy+8UOeff76qVq2q8PBwtW/f3rm+aNhas2aNPvzwQ+3bt0/79u3T/fffr6+++kr79+/XTTfdpNtvv91l2zVq1NDLL7+siy66SBEREcrJydHmzZuVn5+vTz/9VP/5z38kSRdffLF8fX314IMP6sMPP1S9evV0/fXX6/zzz9dff/2lGjVqKCIiQn///bdq1KghHx8f5w8qQUFBpz1WjDHO/+OnU/SHCOl4q8e+fftkjFHNmjU1YcIEPfPMM84y33//vQ4dOqTly5crJiZGNWrUUGhoqIKDg53d/CRp9+7dCg8PV9WqVbVlyxa1bdtWhw8fVvXq1ZWfny9/f3+lpqbqq6++UnR0tBYvXuz8ce32228/6T7ExsZq69at6tu3rz788EOlpqZq9+7dioiIUP369bVnzx6tWrVKCxYsKNbVsOjJ/AUXXKALL7xQ0vGRD3fu3KnPPvtMDRs21NGjR9WuXTu9/fbbzmM7JSVFP/zwg/bt26ft27c7f4jr3r27Lr30Un377bd64403FBoaqt9//11Vq1bV7t275e/vr3379ikwMFArVqzQbbfdpu3bt+v6669XjRo15HA4dOzYMd1+++3avn27evfurfj4eA0bNky//vqrQkJCNHbsWP30009q0aKF/vzzT+3cubPEQS4OHTqkBg0a6OOPP1bHjh2Vn58vY4x8fX3lcDhcuk4XfjcnJyerVq1azh8eGzduXOL7v379et1zzz2qUqWKrr32WoWGhqpq1apq2rSptmzZoqeeekq//vpriX+3QgsXLlSLFi0UGRmpO++80+XzZs6cOdqwYYNLoG/atKlatWqldu3a6fHHH3cGqAsvvNClFS88PFy9evXSjz/+6NJaGxcXp9WrVyswMFCXXnqpwsLCtHr1apfzks2bN6tVq1b67bff9NZbb2ngwIGqV6+ekpKSVKdOHZfzmTO1c+dOtW3bVhEREdq2bdtJ/w8bY3To0CHVrFlTDodDxhjl5eUpMzNT1apV05YtWxQZGamoqCh9+OGH+vzzz/Xss8+WeI5xsu3n5+c7P+N27dqlPXv2KDk5WU8//bTS0tL0xRdfaNGiRerWrZuqV6+utWvXqmrVqurZs6fL8WCM0bFjxxQcHKy0tDTdcccduvLKKzVw4EAtXLhQSUlJysnJUatWrfTZZ5+pefPmevXVV5WZmamVK1dq8eLF+uGHH7R06VIdPHhQHTp00I8//ijJ6vbbvHlzJSUlqXr16goPD9fKlSu1ZMkSDR8+XAUFBapWrVqxcyLJCrVTpkzR7bffrvDwcEVHRxc7jjMzM/XWW2/p8ssvV7t27VwCXl5envbv3686deroscce06WXXqo+ffq47PfChQtVs2ZNdezY0WW7GRkZzhuxS9Ls2bM1cuRISdL48eM1ceJEtWvXTjNmzNBff/3l/Ky67LLLVLt2bUnW9/J///tfLVmyRIsWLSq2fz4+Pho2bJjzBvTGGH377bf67bffNHbsWGfPhK+//lpNmzZV+/btlZSUpPj4+LP+IUMGp5WammokmdTU1FOWu/HGG40k57Rjx45iZbKzs829995rbrvtNpOXl2eMMebLL780rVq1MoMHDzY//fSTMca4bKe00x9//GEuvfRS5+NGjRqd1XbK43T//febqVOnmieffNKMGjXK/Prrr6ZDhw7O9d26dSvVdmrUqGHuueceExISYvs+laauZ/vcqlWrmgEDBpgnn3zSTJ482ezfv99s3rzZPDhihDGSMZKpWaWK87j86KOPjCTTpEkTU7t2bSPJfP3112bOnDlm2rRpJb7G7bffbpo0aeJ8/MEHHzjnr732WnPdddeZn376yWzdutXcdtttznX9+/c3devWdT5u166dy3bvu+8+s3LlSnPvvfea3377zYwYMcLUqVPHdOrUyezbt8/07dvXdOnSxeU5Q4YMMY8//rg5evSomTx5sgkMDHTuT7Nmzcz48ePN3XffXar3rTTvr6+vr2nZsqXtx0h5ncLDw40xxkyfPt32uhROUVFRxZYlJiaaBx54wNx6662mZ8+epkePHs5jp3AaOHCgadOmjQkJCTGTJk0yAwcONJ988ol5//33zcSJE83ixYud227QoEGZ1DUiIqLUZa+88krTqlWrYsv9/PxMcnLySb+vTvb/+kymxo0bm4MHDzq3WVBQYOvf2M/Pz7bXPvG4KazPXXfdZQYOHGguuuiic36N66+/3ixYsMDW97is3qfevXubKVOmmN69e5vLL7/cfPTRR2bt2rXmpptuMhERESY+Pt4EBQWdclujR4829erVc1kWFhZ21nW7+uqrzUMPPWT69+9v+/vkyal+/fpm5syZpkWLFubmm28u8Vi+9dZbz+k7b8aMGWb58uUmKirK1KlTx7n8P//5j+nYsaO56aabzHPPPXdO+zF48GBzzTXXePz9O102MMYYwlYplCZs5efnF/sDbNu2zRhjTF5enrn//vvNe++957L+ww8/NF988UWx5y1atMj2/3xMlXMKlpxhK1gyV1xxxRk9/4knnihx+Z133mmOHj1q+/4xnf3Up0+fUpWLi4sz1apVO2WZoKAgk5mZaWrVqmX7fnlqCgsLM5s3by7Vyd5dd93ltnr4+PiYvn37mj179pzyey0/P9+8+OKLpd43SWblypXmgw8+MC1btjRDhw41u3fvLrbdJ5980va/xZlOHTp0MO3bt3c+jo6Otr1OJU2F4fn5558vcf3tt99uHnzwwXN6jXbt2pnatWubyMjIYuvOP/98c8UVV5w2BDF5dnI4HLbXoaJMdevWNQ0bNizTbZYmbDmM8ZKLEs5BWlqawsLClJqaetJm8J9++kkdOnQotnzw4MF6//33z2kAh4ULF+r6668vsxsJ161bt8SLol9//XUlJiaqU6dOql69uo4dO6bVq1frySefPKPtDx48WAEBAXrrrbecy6pXr+7s+zp69GjNmDHD5Tm9e/cudu1PeXLrrbdq586dLtctlEarVq3k4+Pj0gUzMjLSeS1e0f3u1KmTpkyZoi5dupRdxU8QLKmws2uIpMxTlD3R77//rpiYGJcm/kLp6ekKCQnRZ599puuvv74MauperVq10vbt2097jd7MmTP173//W7m5uXrggQfUoEEDZ5cGSbrmmmtUv359vfzyy6V+bX9/f7Vo0UJ33HGHbr75Zn388cdKSkrSNddco5SUFC1btkz169fXgQMHVL9+fW3fvl2TJk2SZHXVvOaaa5Senq5GjRopOztbSUlJevzxx9W/f39lZGQ4u0289dZb+vLLL7Vnzx698847Ou+887RkyRJFRESoXbt2WrFihTp06KBJkybp2muvVdeuXWWM0S+//KJXXnlFr776qiRpy5Ytatmy5Un3Z/PmzXrqqaeUl5fn0p3u008/1fXXX6+YmBjt2rVLgYGBMsZo5cqVql+/vgICAvTDDz+oS5cuCgkJUXJysmJiYlRQUKCPP/5Y//vf/3TRRRdp7NixkqxrAwu76gYEBEiSS1e/o0eP6sUXX9R5552nK664Qh07dlS9evU0ffp0zZ8/X7t379Y333yjvn37OgdyKEsvvfSSRo4cqaNHj6p379767rvvFBERoX/84x/atWuX7rrrLnXu3FkhISEKCAjQL7/8omPHjjk/C+bPn6/Bgwfr0KFD+uqrr3Tw4EGlpaWpV69e+u6777R48WJdffXV6tevnxwOh5KTk7Vu3Tq98847qlq1qrNbWXh4+BmNBllQUKDvvvtOf/31l5o3b66wsDBFRkZKkvN9PhPGGP3555/avXu3YmNjdezYMSUkJKhLly46evSoatasqbCwMP3+++8u10oUPjcnJ0c5OTnasmWLkpOTFRcXp+3bt8vX11exsbHy9/fX0aNH1ahRIx06dEhBQUHO16tevbokqV27dkpISFDz5s21efNm1apVS40aNXJeS1a3bl2XrmYn249Dhw7p66+/1sqVK+Xj46OkpCSXa0Sfe+45xcXFaerUqbr00ksVEhKi//3vf/rhhx+UmpqqZs2a6dixY4qMjJSPj4/279+vP/74Q126dNGsWbOUmpqqzz77TMOGDVPVqlWd/0eqVavmfI2MjAz5+vo6h+svVNgFfvLkyfrzzz81cuRIXXTRRXI4HNq+fbuqVKmi6OhopaWlqUaNGs5tbd++XbVr11ZERIQcDocKCgrk4+Oj/Pz8kw5sk5eXp2PHjrnU68iRIwoODnbes23Lli1au3atfH19ncfR0KFDddddd+mWW27R6tWr9d577+m1115z2bavr+9pr7muVauWWrZsqcjISAUFBenrr79WUlKSqlSp4vIZXq9ePU2ePFlt2rRRZmampk2bpk2bNmnfvn3Kzs7W1VdfrREjRig1NVUrVqxQ//79tWvXLn3xxRf6/vvvS7yu3M/Pr9j5W0xMjPz9/eXn56dx48bpyiuvVGxsrFq0aOG8LvPtt9/Wjh07tHfvXgUGBurrr7/W7t27Xbbp6+urGTNmaPLkyQoODtZbb72lRYsWKSoqSq1atVKbNm1kjNHBgweVm5ur7OxsBQQEKCYmRnPmzFFMTIxuueUW57GRm5urvXv3qkGDBnrxxRe1fPlyNW7cWKmpqbrooot03nnnadasWWrUqNFJz+uCg4NVt25d+fj4KD09/ZQD6AwcOFDz5893jsT95JNPasKECWU2YFnheXdISIgiIiKcN2q/7rrrFBcXpz59+ui7777T8OHDiz23aPffBg0aFFv/3HPP6bXXXtOUKVO0c+dO/e9//9Odd96pCy+8UE2bNnWW++OPP7Rr1y799ttvGjNmjMaNG6cxY8YoMzNTmZmZioyMPGU2cDptHIOzZWvJkiVm4sSJJisryxhjdZlYt26dyc7OdmmhOrFZ+2ynrl27OutwslaDvLw8k5aW5rLs6quvNomJiSWWLygoMMYYc+zYMZOQkGByc3NL9R5s2LChxO3NmDHDbNiwwRw7dsx89913Jj8/3xhjTEZGhlmwYIHL9jMzM53zubm55r777jMff/xxia/3xx9/mOTkZHPo0CEzbNgwM3HiRDNgwAAjyYSGhpZYl++++875/Llz5xpJ5sknn3TZ7vr1680bb7xhJk2aZEaOHGleffVVk52dbbZv3+58n8eMGWO2bt3qfK+K2rZtm1m0aJHzcWJionMfCwoKzN9//12q97Ooffv2mbS0tBLXHTt2zBQUFDjrkpaWZg4cOGA2bdpkJk6caPbu3WsOHDhgdu/ebV5//XWTlpZmkpOTzcGDB8327dvNxx9/bHr16uV8j05s2Sr6/rVo0cLMmTPH3HHHHc4umv369TNTp041f/75p7NOKSkp5o8//jCvvfaamTFjhvNvboz1a/m1115rLrroIvPJJ5+YTz75xLz55ptm0qRJ5uDBg2b37t3m119/NVOnTjVdu3Y1GzduNG+99ZbZuHGjGTVqlLnsssuMZHWXmTlzphk5cqT5/PPPzR9//OFS1xo1apgVK1YYY4y59957TYsWLczbb79t7rjjDvPMM8+YNWvWmEOHDpnvv//e7Ny504wcOdK8++67Jicn56R/hyNHjpTqFypjjNm6dav566+/XJbt3bvX+Xc8cuSIOXbsWKm2VZby8/PNjz/+aLKzsz3+2oWv7+/vbyQ5j7thw4bZUpfTeffdd81tt91mDhw44Py7FXbtPlFBQYFJS0szeXl5Zv/+/Wb37t1m+/bt5tixY87PlKL/Dzzp8OHDxY5FoDzLy8szc+bMMatWrXIuO3DggFmwYIHZuHFjid+9npCenn7S74jt27ebjz766JR1y8nJcZ4flmT9+vXl8v9qVlaW+eabb8zRo0dLXP/UU0+Z6tWrmw0bNpx2W5mZmSYzM9Ps2bPH/PbbbyYvL8989dVXzu/uTz75xMycOdOsXbvWFBQUuJwvz54927mdkj6Lly9fbt58881iy4t2rW3UqNFJ96OslfYSI2PoRlgqhW9oaaagoKCzulYqPz+/2DUFb7zxhks9Pv30U5f1t956q3Pdu+++61z+ww8/GGOM+eabb1zKv/LKK+f0Ppx4/cXZBIuylJiYaBo3blwsSOIU0tNdwlZUVJSzu2t5d/jwYfP++++fMjDBfk2bNnX5nPjss8/srhIAwEY7d+4069atK7Z89erV5uGHHzb79+8/p+2vW7fOTJ482aPnpWcStuhGWAqF3QhLIyYm5rQ3G3zjjTc0btw4/f3335KsEcAWLVqk5ORkRUVFOcu98847Gjx4sMtzhw8f7hyNa9++fapbt64kq8vDjTfeqPT0dH3xxRfOEc/i4uK0dq01/tx3332nSy+9tHQ7XYL8/Hx9/PHHmj17tu6++2717dv3rLdVVg4fPqyff/5ZXbp08brhps9KRoZU2A0wPV0KCbG3Pqh0pkyZoocffliSVKdOHe3bt6/Uo/oBAFARlOYSo0KErVI4k7CVkpKiunXrOu81VDj06IMPPqhp06ZJkv7880+lpaU5+6t//fXX6tGjhyQ5hyeWpL///ls1a9Z02X5+fr4GDBigsLAwvfnmm6c9icnPz9fFF18sY4x+/PHHs+qDj0okI0Mq7L+8ezdhC2Xu0KFDioiIkCS99tpr+te//mVzjQAAKFtnErZoCigDhYM/PPHEEwoLC9OiRYt044036rnnntP555+v6tWrKyYmRnFxcerUqZNq1aql2rVrKz8/XwcPHlRMTIxzW4cPH9bs2bPVpUuXYkFLsi4m/e9//1vquvn6+mr9+vUqKCg46cWv8CIhIVKRG2UDZa1mzZrKzMzU6tWr1a1bN7urAwCArWjZKoWi6XXBggVas2aNnnzySc2fP199+/Y9pxv1AQAAAKg46EZYxs7kDQUAAABQeZ1JNij9zTgAVHzHjkldu1pTGd23DQAAACXjmi3AmxQUSCtWHJ8HAACA29CyBQAAAABu4FVh65VXXlGDBg1UpUoVdezY0Xn/KQAAAAAoa14Ttj788EONHj1a48eP108//aS2bdsqPj5eycnJdlcNAAAAQCXkNWFrxowZuuOOO3TbbbepRYsWeu211xQcHKxZs2bZXTUAAAAAlZBXhK2cnBxt2LBB3bt3dy7z8fFR9+7dtXr16mLls7OzlZaW5jIBAAAAwJnwirD1999/Kz8/X1FRUS7Lo6KilJSUVKz8lClTFBYW5pxiY2M9VVXA/YKDrQkAAABu5RVh60yNGzdOqampzmnfvn12VwkoGyEhUkaGNYWE2F0bAACASs0r7rMVEREhX19fHTx40GX5wYMHFR0dXax8YGCgAgMDPVU9AAAAAJWQV7RsBQQEqEOHDlq6dKlzWUFBgZYuXaq4uDgbawYAAACgsvKKli1JGj16tIYMGaILL7xQF198sZ5//nllZGTotttus7tqgOdkZUn9+1vzn3wiValib30AAAAqMa8JWwMGDNBff/2lxx9/XElJSWrXrp0WLVpUbNAMoFLLz5e++ur4PAAAANzGYYwxdleivEtLS1NYWJhSU1MVGhpqd3WAs5eRIVWtas2npzNIBgAAwBk6k2zgFddsAQAAAICnEbYAAAAAwA0IWwAAAADgBoQtAAAAAHADrxmN8FwUjiGSlpZmc02Ac5SRcXw+LY0RCQEAAM5QYSYozTiDhK1SOHr0qCQpNjbW5poAZah2bbtrAAAAUGEdPXpUYWFhpyzD0O+lUFBQoP3796tatWpyOBweec20tDTFxsZq3759DDePM8Kxg7PFsYOzxbGDs8Wxg3Nh1/FjjNHRo0dVu3Zt+fic+qosWrZKwcfHR3Xr1rXltUNDQ/nwwVnh2MHZ4tjB2eLYwdni2MG5sOP4OV2LViEGyAAAAAAANyBsAQAAAIAbELbKqcDAQI0fP16BgYF2VwUVDMcOzhbHDs4Wxw7OFscOzkVFOH4YIAMAAAAA3ICWLQAAAABwA8IWAAAAALgBYQsAAAAA3ICwBQAAAABuQNjykKlTp8rhcOi+++5zLsvKytKIESNUs2ZNVa1aVf3799fBgwddnrd371716dNHwcHBqlWrlsaMGaO8vDyXMsuXL1f79u0VGBioJk2aaM6cOR7YI7jLhAkT5HA4XKZmzZo513Pc4HT+/PNP3XLLLapZs6aCgoLUunVrrV+/3rneGKPHH39cMTExCgoKUvfu3fX777+7bOPw4cMaNGiQQkNDFR4ermHDhik9Pd2lzKZNm3T55ZerSpUqio2N1bRp0zyyf3CPBg0aFPvscTgcGjFihCQ+e3By+fn5euyxx9SwYUMFBQWpcePGmjRpkoqOwcbnDk7m6NGjuu+++1S/fn0FBQWpU6dOWrdunXN9hT92DNxu7dq1pkGDBqZNmzbm3nvvdS4fPny4iY2NNUuXLjXr1683l1xyienUqZNzfV5enmnVqpXp3r27+fnnn81XX31lIiIizLhx45xldu7caYKDg83o0aPNtm3bzEsvvWR8fX3NokWLPLmLKEPjx483LVu2NAcOHHBOf/31l3M9xw1O5fDhw6Z+/fpm6NChZs2aNWbnzp3m66+/Ntu3b3eWmTp1qgkLCzPz5883GzduNNdee61p2LChOXbsmLNMz549Tdu2bc2PP/5oVq1aZZo0aWJuvvlm5/rU1FQTFRVlBg0aZLZs2WLef/99ExQUZP7zn/94dH9RdpKTk10+dxYvXmwkmWXLlhlj+OzByT311FOmZs2aZsGCBWbXrl3m448/NlWrVjUvvPCCswyfOziZG2+80bRo0cKsWLHC/P7772b8+PEmNDTU/PHHH8aYin/sELbc7OjRo+a8884zixcvNl26dHGGrZSUFOPv728+/vhjZ9lffvnFSDKrV682xhjz1VdfGR8fH5OUlOQsM3PmTBMaGmqys7ONMcY8+OCDpmXLli6vOWDAABMfH+/mPYO7jB8/3rRt27bEdRw3OJ2xY8eayy677KTrCwoKTHR0tJk+fbpzWUpKigkMDDTvv/++McaYbdu2GUlm3bp1zjILFy40DofD/Pnnn8YYY1599VVTvXp15zFV+NpNmzYt612CTe69917TuHFjU1BQwGcPTqlPnz7m9ttvd1l2/fXXm0GDBhlj+NzByWVmZhpfX1+zYMECl+Xt27c3jzzySKU4duhG6GYjRoxQnz591L17d5flGzZsUG5ursvyZs2aqV69elq9erUkafXq1WrdurWioqKcZeLj45WWlqatW7c6y5y47fj4eOc2UDH9/vvvql27tho1aqRBgwZp7969kjhucHpffPGFLrzwQt1www2qVauWLrjgAr3xxhvO9bt27VJSUpLL3z8sLEwdO3Z0OYbCw8N14YUXOst0795dPj4+WrNmjbNM586dFRAQ4CwTHx+vxMREHTlyxN27CTfLycnRe++9p9tvv10Oh4PPHpxSp06dtHTpUv3222+SpI0bN+q7775Tr169JPG5g5PLy8tTfn6+qlSp4rI8KChI3333XaU4dghbbvTBBx/op59+0pQpU4qtS0pKUkBAgMLDw12WR0VFKSkpyVmm6JdW4frCdacqk5aWpmPHjpXVrsCDOnbsqDlz5mjRokWaOXOmdu3apcsvv1xHjx7luMFp7dy5UzNnztR5552nr7/+WnfddZf+/e9/6+2335Z0/Bgo6e9f9PioVauWy3o/Pz/VqFHjjI4zVFzz589XSkqKhg4dKonvLJzaQw89pJtuuknNmjWTv7+/LrjgAt13330aNGiQJD53cHLVqlVTXFycJk2apP379ys/P1/vvfeeVq9erQMHDlSKY8fPrVv3Yvv27dO9996rxYsXF0vrwKkU/hIoSW3atFHHjh1Vv359ffTRRwoKCrKxZqgICgoKdOGFF2ry5MmSpAsuuEBbtmzRa6+9piFDhthcO1QUb731lnr16qXatWvbXRVUAB999JHmzp2refPmqWXLlkpISNB9992n2rVr87mD03r33Xd1++23q06dOvL19VX79u118803a8OGDXZXrUzQsuUmGzZsUHJystq3by8/Pz/5+flpxYoVevHFF+Xn56eoqCjl5OQoJSXF5XkHDx5UdHS0JCk6OrrYSE+Fj09XJjQ0lBPzSiI8PFznn3++tm/frujoaI4bnFJMTIxatGjhsqx58+bOrqiFx0BJf/+ix0dycrLL+ry8PB0+fPiMjjNUTHv27NGSJUv0z3/+07mMzx6cypgxY5ytW61bt9bgwYM1atQoZ88ePndwKo0bN9aKFSuUnp6uffv2ae3atcrNzVWjRo0qxbFD2HKTK6+8Ups3b1ZCQoJzuvDCCzVo0CDnvL+/v5YuXep8TmJiovbu3au4uDhJUlxcnDZv3uxyAC1evFihoaHOk6m4uDiXbRSWKdwGKr709HTt2LFDMTEx6tChA8cNTunSSy9VYmKiy7LffvtN9evXlyQ1bNhQ0dHRLn//tLQ0rVmzxuUYSklJcflV8dtvv1VBQYE6duzoLLNy5Url5uY6yyxevFhNmzZV9erV3bZ/cL/Zs2erVq1a6tOnj3MZnz04lczMTPn4uJ5S+vr6qqCgQBKfOyidkJAQxcTE6MiRI/r666/Vt2/fynHsuH0IDjgVHY3QGGsY3Xr16plvv/3WrF+/3sTFxZm4uDjn+sJhdHv06GESEhLMokWLTGRkZInD6I4ZM8b88ssv5pVXXmEY3Qru/vvvN8uXLze7du0y33//venevbuJiIgwycnJxhiOG5za2rVrjZ+fn3nqqafM77//bubOnWuCg4PNe++95ywzdepUEx4ebj7//HOzadMm07dv3xKH0b3gggvMmjVrzHfffWfOO+88l2F0U1JSTFRUlBk8eLDZsmWL+eCDD0xwcDBDMFdw+fn5pl69embs2LHF1vHZg5MZMmSIqVOnjnPo908//dRERESYBx980FmGzx2czKJFi8zChQvNzp07zTfffGPatm1rOnbsaHJycowxFf/YIWx50Ilh69ixY+buu+821atXN8HBwaZfv37mwIEDLs/ZvXu36dWrlwkKCjIRERHm/vvvN7m5uS5lli1bZtq1a2cCAgJMo0aNzOzZsz2wN3CXAQMGmJiYGBMQEGDq1KljBgwY4HKPJI4bnM7//vc/06pVKxMYGGiaNWtmXn/9dZf1BQUF5rHHHjNRUVEmMDDQXHnllSYxMdGlzKFDh8zNN99sqlatakJDQ81tt91mjh496lJm48aN5rLLLjOBgYGmTp06ZurUqW7fN7jX119/bSQVOx6M4bMHJ5eWlmbuvfdeU69ePVOlShXTqFEj88gjj7gMs83nDk7mww8/NI0aNTIBAQEmOjrajBgxwqSkpDjXV/Rjx2FMkdt7AwAAAADKBNdsAQAAAIAbELYAAAAAwA0IWwAAAADgBoQtAAAAAHADwhYAAAAAuAFhCwAAAADcgLAFAAAAAG5A2AIAAAAANyBsAQBwGhMmTFC7du3sroaTw+HQ/Pnz7a4GAOA0CFsAgHLjtddeU7Vq1ZSXl+dclp6eLn9/f3Xt2tWl7PLly+VwOLRjxw4P19JzylvIAwCcGcIWAKDc6Natm9LT07V+/XrnslWrVik6Olpr1qxRVlaWc/myZctUr149NW7c2I6qAgBwWoQtAEC50bRpU8XExGj58uXOZcuXL1ffvn3VsGFD/fjjjy7Lu3XrpnfffVcXXnihqlWrpujoaA0cOFDJycmSpIKCAtWtW1czZ850eZ2ff/5ZPj4+2rNnjyQpJSVF//znPxUZGanQ0FBdccUV2rhx4ynr+uabb6p58+aqUqWKmjVrpldffdW5bvfu3XI4HPr000/VrVs3BQcHq23btlq9erXLNt544w3FxsYqODhY/fr104wZMxQeHi5JmjNnjiZOnKiNGzfK4XDI4XBozpw5zuf+/fff6tevn4KDg3Xeeefpiy++KPX7DADwDMIWAKBc6datm5YtW+Z8vGzZMnXt2lVdunRxLj927JjWrFmjbt26KTc3V5MmTdLGjRs1f/587d69W0OHDpUk+fj46Oabb9a8efNcXmPu3Lm69NJLVb9+fUnSDTfcoOTkZC1cuFAbNmxQ+/btdeWVV+rw4cMl1nHu3Ll6/PHH9dRTT+mXX37R5MmT9dhjj+ntt992KffII4/ogQceUEJCgs4//3zdfPPNzi6S33//vYYPH657771XCQkJuuqqq/TUU085nztgwADdf//9atmypQ4cOKADBw5owIABzvUTJ07UjTfeqE2bNql3794aNGjQSesLALCJAQCgHHnjjTdMSEiIyc3NNWlpacbPz88kJyebefPmmc6dOxtjjFm6dKmRZPbs2VPs+evWrTOSzNGjR40xxvz888/G4XA4y+bn55s6deqYmTNnGmOMWbVqlQkNDTVZWVku22ncuLH5z3/+Y4wxZvz48aZt27Yu6+bNm+dSftKkSSYuLs4YY8yuXbuMJPPmm28612/dutVIMr/88osxxpgBAwaYPn36uGxj0KBBJiwszPn4xNctJMk8+uijzsfp6elGklm4cGGxsgAA+9CyBQAoV7p27aqMjAytW7dOq1at0vnnn6/IyEh16dLFed3W8uXL1ahRI9WrV08bNmzQNddco3r16qlatWrq0qWLJGnv3r2SpHbt2ql58+bO1q0VK1YoOTlZN9xwgyRp48aNSk9PV82aNVW1alXntGvXrhIH38jIyNCOHTs0bNgwl/JPPvlksfJt2rRxzsfExEiSs4tjYmKiLr74YpfyJz4+laLbDgkJUWhoqHPbAIDywc/uCgAAUFSTJk1Ut25dLVu2TEeOHHGGp9q1ays2NlY//PCDli1bpiuuuEIZGRmKj49XfHy85s6dq8jISO3du1fx8fHKyclxbnPQoEGaN2+eHnroIc2bN089e/ZUzZo1JVmjHZ54nVihwuunikpPT5dkXW/VsWNHl3W+vr4uj/39/Z3zDodDknUdWVkouu3C7ZfVtgEAZYOwBQAod7p166bly5fryJEjGjNmjHN5586dtXDhQq1du1Z33XWXfv31Vx06dEhTp05VbGysJLmMZFho4MCBevTRR7Vhwwb997//1WuvveZc1759eyUlJcnPz08NGjQ4bd2ioqJUu3Zt7dy5U4MGDTrrfWzatKnWrVvnsuzExwEBAcrPzz/r1wAA2IuwBQAod7p166YRI0YoNzfX2bIlSV26dNHIkSOVk5Ojbt26yc/PTwEBAXrppZc0fPhwbdmyRZMmTSq2vQYNGqhTp04aNmyY8vPzde211zrXde/eXXFxcbruuus0bdo0nX/++dq/f7++/PJL9evXTxdeeGGx7U2cOFH//ve/FRYWpp49eyo7O1vr16/XkSNHNHr06FLt4z333KPOnTtrxowZuuaaa/Ttt99q4cKFzhawwnrv2rVLCQkJqlu3rqpVq6bAwMAzeSsBADbimi0AQLnTrVs3HTt2TE2aNFFUVJRzeZcuXXT06FHnEPGRkZGaM2eOPv74Y7Vo0UJTp07VM888U+I2Bw0apI0bN6pfv34KCgpyLnc4HPrqq6/UuXNn3XbbbTr//PN10003ac+ePS6vXdQ///lPvfnmm5o9e7Zat26tLl26aM6cOWrYsGGp9/HSSy/Va6+9phkzZqht27ZatGiRRo0apSpVqjjL9O/fXz179lS3bt0UGRmp999/v9TbBwDYz2GMMXZXAgAASHfccYd+/fVXrVq1yu6qAADKAN0IAQCwyTPPPKOrrrpKISEhWrhwod5++22XmyMDACo2WrYAALDJjTfeqOXLl+vo0aNq1KiR7rnnHg0fPtzuagEAyghhCwAAAADcgAEyAAAAAMANCFsAAAAA4AaELQAAAABwA8IWAAAAALgBYQsAAAAA3ICwBQAAAABuQNgCAAAAADcgbAEAAACAG/wfaWyUfAIdR2kAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from astropy.io import fits\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def read_spectrum(fits_file):\n",
    "    with fits.open(fits_file) as hdul:\n",
    "        coadd_data = hdul[1].data\n",
    "        loglam = coadd_data['loglam'] \n",
    "        wavelength = 10 ** loglam\n",
    "        flux = coadd_data['FLUX']  \n",
    "        return wavelength, flux\n",
    "\n",
    "def calculate_h_delta(wavelength, flux):\n",
    "    l_wavelength = 4861 * (1 + 0.04)\n",
    "    h_delta_range = (l_wavelength - 1, l_wavelength + 1)  \n",
    "    continuum_range = (l_wavelength-50, l_wavelength-10), (l_wavelength+10, l_wavelength+50) \n",
    "\n",
    "    h_delta_mask = (wavelength >= h_delta_range[0]) & (wavelength <= h_delta_range[1])\n",
    "    h_delta_flux = flux[h_delta_mask]\n",
    "\n",
    "    continuum_flux = []\n",
    "    for r in continuum_range:\n",
    "        cont_mask = (wavelength >= r[0]) & (wavelength <= r[1])\n",
    "        continuum_flux.append(np.mean(flux[cont_mask]))    \n",
    "    \n",
    "    continuum_mean = np.mean(continuum_flux)\n",
    "\n",
    "    h_delta_index = np.mean(h_delta_flux) / continuum_mean \n",
    "    h_delta_error = np.std(h_delta_flux) / continuum_mean\n",
    "\n",
    "    return h_delta_index, h_delta_error, l_wavelength\n",
    "\n",
    "\n",
    "fits_file = \"spec-1678-53433-0425.fits\"  \n",
    "wavelength, flux = read_spectrum(fits_file)\n",
    "\n",
    "h_delta_index, h_delta_error, h_wl = calculate_h_delta(wavelength, flux)\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(wavelength, flux, label='Spectrum', color='black')\n",
    "plt.axvline(x=h_wl, color='red', linestyle='--', label='Hβ Center')\n",
    "plt.xlim(np.min(wavelength), np.max(wavelength))\n",
    "\n",
    "plt.text(h_wl+100, np.max(flux) * 0.8, 'Hβ Lick Index: %.4f ± %.4f' %(h_delta_index, h_delta_error),\n",
    "            color='red', fontsize=10)\n",
    "\n",
    "plt.xlabel('Wavelength')\n",
    "plt.ylabel('Flux')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SDSS的光谱分辨率是~2000，某星系其连续谱单位波长的典型信噪比为5，该星系Ha线的等值宽度为50A，请问其Ha线流量的信噪比是多少？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "76.18467164406522\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "SN_cont = 5\n",
    "R = 2000\n",
    "EW_H = 50\n",
    "lam_H = 6563\n",
    "\n",
    "det_lam = lam_H / R\n",
    "SN_H = SN_cont * EW_H / det_lam\n",
    "print(SN_H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#信噪比是按照流量倍数开根号增加的"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用GEHONG生成一个积分视场光谱，并进行简单的可视化（选做）"
   ]
  }
 ],
 "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
}
