07_NEO_Appearances.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "61e33bd7-a70c-4cd7-814c-5037418d0bbf",
   "metadata": {},
   "source": [
    "# Today's NEO Appearances\n",
    "\n",
    "In this session we compute the overall statistics of the current NEO populations' brightness. We'll see that we have a few weird outliers and compare their \"detectability\" with a generic property of the NEO survey telescope Pan-Starrs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "50dc0a3f-2a87-42e8-ad0e-3209752d49bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard libraries\n",
    "import datetime\n",
    "import math\n",
    "import pathlib\n",
    "import sqlite3\n",
    "import sys\n",
    "\n",
    "# Installed libraries\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import spiceypy\n",
    "\n",
    "# Append to root directory of this repository\n",
    "sys.path.append(\"../\")\n",
    "\n",
    "# Auxiliary module that contains the apparent magnitude\n",
    "from auxiliary import photometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3d760111-fc47-4296-8fa3-a0fc17f76abd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Accessing the NEO database\n",
    "database_dir = pathlib.Path(\"../databases/neos/\")\n",
    "database_file = pathlib.Path(\"neodys.db\")\n",
    "database_filepath = database_dir / database_file\n",
    "\n",
    "# Establish a connection to the database and set a cursor\n",
    "neodys_db_con = sqlite3.connect(database_filepath)\n",
    "neodys_db_cur = neodys_db_con.cursor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "24a1558d-6da8-4873-b556-b9c050e2be55",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all information from the DB. Since the DB is rather small, this won't cause any issues!\n",
    "neo_df = pd.read_sql(\"SELECT * FROM main\", neodys_db_con)\n",
    "\n",
    "# Close the database.\n",
    "neodys_db_con.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bb617cf6-0afe-481e-b78f-f45259823847",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Name</th>\n",
       "      <th>Epoch_MJD</th>\n",
       "      <th>SemMajAxis_AU</th>\n",
       "      <th>Ecc_</th>\n",
       "      <th>Incl_deg</th>\n",
       "      <th>LongAscNode_deg</th>\n",
       "      <th>ArgP_deg</th>\n",
       "      <th>MeanAnom_deg</th>\n",
       "      <th>AbsMag_</th>\n",
       "      <th>SlopeParamG_</th>\n",
       "      <th>Aphel_AU</th>\n",
       "      <th>Perihel_AU</th>\n",
       "      <th>NEOClass</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>433</td>\n",
       "      <td>59600.0</td>\n",
       "      <td>1.458273</td>\n",
       "      <td>0.222727</td>\n",
       "      <td>10.828461</td>\n",
       "      <td>304.296346</td>\n",
       "      <td>178.897177</td>\n",
       "      <td>246.904119</td>\n",
       "      <td>10.87</td>\n",
       "      <td>0.46</td>\n",
       "      <td>1.783070</td>\n",
       "      <td>1.133476</td>\n",
       "      <td>Amor</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>719</td>\n",
       "      <td>59600.0</td>\n",
       "      <td>2.637505</td>\n",
       "      <td>0.546959</td>\n",
       "      <td>11.575272</td>\n",
       "      <td>183.855259</td>\n",
       "      <td>156.227573</td>\n",
       "      <td>278.197185</td>\n",
       "      <td>15.53</td>\n",
       "      <td>0.15</td>\n",
       "      <td>4.080112</td>\n",
       "      <td>1.194899</td>\n",
       "      <td>Amor</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>887</td>\n",
       "      <td>59600.0</td>\n",
       "      <td>2.473158</td>\n",
       "      <td>0.570486</td>\n",
       "      <td>9.394141</td>\n",
       "      <td>110.428769</td>\n",
       "      <td>350.493276</td>\n",
       "      <td>86.607207</td>\n",
       "      <td>13.33</td>\n",
       "      <td>-0.12</td>\n",
       "      <td>3.884060</td>\n",
       "      <td>1.062256</td>\n",
       "      <td>Amor</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1036</td>\n",
       "      <td>59600.0</td>\n",
       "      <td>2.665849</td>\n",
       "      <td>0.533122</td>\n",
       "      <td>26.677908</td>\n",
       "      <td>215.517191</td>\n",
       "      <td>132.428118</td>\n",
       "      <td>140.654843</td>\n",
       "      <td>9.41</td>\n",
       "      <td>0.30</td>\n",
       "      <td>4.087071</td>\n",
       "      <td>1.244627</td>\n",
       "      <td>Amor</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1221</td>\n",
       "      <td>59600.0</td>\n",
       "      <td>1.918731</td>\n",
       "      <td>0.435843</td>\n",
       "      <td>11.883264</td>\n",
       "      <td>171.320251</td>\n",
       "      <td>26.643922</td>\n",
       "      <td>261.044505</td>\n",
       "      <td>17.40</td>\n",
       "      <td>0.15</td>\n",
       "      <td>2.754996</td>\n",
       "      <td>1.082467</td>\n",
       "      <td>Amor</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29044</th>\n",
       "      <td>2022KY2</td>\n",
       "      <td>59600.0</td>\n",
       "      <td>1.439076</td>\n",
       "      <td>0.290852</td>\n",
       "      <td>9.562666</td>\n",
       "      <td>220.831102</td>\n",
       "      <td>344.748373</td>\n",
       "      <td>308.992894</td>\n",
       "      <td>24.82</td>\n",
       "      <td>0.15</td>\n",
       "      <td>1.857634</td>\n",
       "      <td>1.020519</td>\n",
       "      <td>Amor</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29045</th>\n",
       "      <td>2022KZ</td>\n",
       "      <td>59600.0</td>\n",
       "      <td>1.404505</td>\n",
       "      <td>0.389038</td>\n",
       "      <td>10.755270</td>\n",
       "      <td>64.554335</td>\n",
       "      <td>243.245936</td>\n",
       "      <td>257.305079</td>\n",
       "      <td>26.44</td>\n",
       "      <td>0.15</td>\n",
       "      <td>1.950911</td>\n",
       "      <td>0.858099</td>\n",
       "      <td>Apollo</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29046</th>\n",
       "      <td>2022KZ1</td>\n",
       "      <td>59600.0</td>\n",
       "      <td>1.643745</td>\n",
       "      <td>0.562916</td>\n",
       "      <td>31.322592</td>\n",
       "      <td>63.673736</td>\n",
       "      <td>262.310279</td>\n",
       "      <td>276.670518</td>\n",
       "      <td>26.14</td>\n",
       "      <td>0.15</td>\n",
       "      <td>2.569035</td>\n",
       "      <td>0.718455</td>\n",
       "      <td>Apollo</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29047</th>\n",
       "      <td>2022KZ2</td>\n",
       "      <td>59600.0</td>\n",
       "      <td>1.342108</td>\n",
       "      <td>0.412424</td>\n",
       "      <td>13.556790</td>\n",
       "      <td>237.807367</td>\n",
       "      <td>283.141272</td>\n",
       "      <td>319.777066</td>\n",
       "      <td>26.25</td>\n",
       "      <td>0.15</td>\n",
       "      <td>1.895626</td>\n",
       "      <td>0.788590</td>\n",
       "      <td>Apollo</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29048</th>\n",
       "      <td>6344P-L</td>\n",
       "      <td>59600.0</td>\n",
       "      <td>2.820120</td>\n",
       "      <td>0.661492</td>\n",
       "      <td>4.679303</td>\n",
       "      <td>182.813414</td>\n",
       "      <td>234.952874</td>\n",
       "      <td>5.480313</td>\n",
       "      <td>20.40</td>\n",
       "      <td>0.15</td>\n",
       "      <td>4.685608</td>\n",
       "      <td>0.954632</td>\n",
       "      <td>Apollo</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>29049 rows × 13 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          Name  Epoch_MJD  SemMajAxis_AU      Ecc_   Incl_deg  \\\n",
       "0          433    59600.0       1.458273  0.222727  10.828461   \n",
       "1          719    59600.0       2.637505  0.546959  11.575272   \n",
       "2          887    59600.0       2.473158  0.570486   9.394141   \n",
       "3         1036    59600.0       2.665849  0.533122  26.677908   \n",
       "4         1221    59600.0       1.918731  0.435843  11.883264   \n",
       "...        ...        ...            ...       ...        ...   \n",
       "29044  2022KY2    59600.0       1.439076  0.290852   9.562666   \n",
       "29045   2022KZ    59600.0       1.404505  0.389038  10.755270   \n",
       "29046  2022KZ1    59600.0       1.643745  0.562916  31.322592   \n",
       "29047  2022KZ2    59600.0       1.342108  0.412424  13.556790   \n",
       "29048  6344P-L    59600.0       2.820120  0.661492   4.679303   \n",
       "\n",
       "       LongAscNode_deg    ArgP_deg  MeanAnom_deg  AbsMag_  SlopeParamG_  \\\n",
       "0           304.296346  178.897177    246.904119    10.87          0.46   \n",
       "1           183.855259  156.227573    278.197185    15.53          0.15   \n",
       "2           110.428769  350.493276     86.607207    13.33         -0.12   \n",
       "3           215.517191  132.428118    140.654843     9.41          0.30   \n",
       "4           171.320251   26.643922    261.044505    17.40          0.15   \n",
       "...                ...         ...           ...      ...           ...   \n",
       "29044       220.831102  344.748373    308.992894    24.82          0.15   \n",
       "29045        64.554335  243.245936    257.305079    26.44          0.15   \n",
       "29046        63.673736  262.310279    276.670518    26.14          0.15   \n",
       "29047       237.807367  283.141272    319.777066    26.25          0.15   \n",
       "29048       182.813414  234.952874      5.480313    20.40          0.15   \n",
       "\n",
       "       Aphel_AU  Perihel_AU NEOClass  \n",
       "0      1.783070    1.133476     Amor  \n",
       "1      4.080112    1.194899     Amor  \n",
       "2      3.884060    1.062256     Amor  \n",
       "3      4.087071    1.244627     Amor  \n",
       "4      2.754996    1.082467     Amor  \n",
       "...         ...         ...      ...  \n",
       "29044  1.857634    1.020519     Amor  \n",
       "29045  1.950911    0.858099   Apollo  \n",
       "29046  2.569035    0.718455   Apollo  \n",
       "29047  1.895626    0.788590   Apollo  \n",
       "29048  4.685608    0.954632   Apollo  \n",
       "\n",
       "[29049 rows x 13 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "neo_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bf552ed8-0511-4def-aa05-25d5921a0f51",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load SPICE kernels\n",
    "spiceypy.furnsh(\"../kernels/spk/de432s.bsp\")\n",
    "spiceypy.furnsh(\"../kernels/lsk/naif0012.tls\")\n",
    "spiceypy.furnsh(\"../kernels/pck/gm_de431.tpc\")\n",
    "\n",
    "# Get the G*M value of the Sun\n",
    "_, gm_sun_pre = spiceypy.bodvcd(bodyid=10, item='GM', maxn=1)\n",
    "gm_sun = gm_sun_pre[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "01537341-0421-43ca-9dfe-cb6672c3ffc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine today's datetime\n",
    "curr_time_et = spiceypy.utc2et(datetime.datetime.now().strftime(\"%Y-%m-%d\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4285e93a-f2ef-43ca-a23f-aed68d8f1618",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.50093320e+08,  1.64316490e+06, -1.29773492e+03])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sun2earth_position_vec = spiceypy.spkgps(targ=399,\n",
    "                                     et=curr_time_et,\n",
    "                                     ref=\"ECLIPJ2000\",\n",
    "                                     obs=10)[0]\n",
    "sun2earth_position_vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2481b84c-15f1-4171-808d-d73589d4535b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For our computations we need to convert some values from AU to km and from deg to rad\n",
    "neo_df.loc[:, \"Perihel_km\"] = neo_df[\"Perihel_AU\"].apply(lambda x: spiceypy.convrt(x, \"AU\", \"km\"))\n",
    "neo_df.loc[:, \"Incl_rad\"] = neo_df[\"Incl_deg\"].apply(lambda x: math.radians(x))\n",
    "neo_df.loc[:, \"LongAscNode_rad\"] = neo_df[\"LongAscNode_deg\"].apply(lambda x: math.radians(x))\n",
    "neo_df.loc[:, \"ArgP_rad\"] = neo_df[\"ArgP_deg\"].apply(lambda x: math.radians(x))\n",
    "neo_df.loc[:, \"MeanAnom_rad\"] = neo_df[\"MeanAnom_deg\"].apply(lambda x: math.radians(x))\n",
    "neo_df.loc[:, \"Epoch_JD\"] = neo_df[\"Epoch_MJD\"].apply(lambda x: x + 2400000.5)\n",
    "neo_df.loc[:, \"Epoch_et\"] = neo_df[\"Epoch_JD\"].apply(lambda x: spiceypy.utc2et(str(x) + \" JD\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e3819761-8751-46f7-aed9-9726e410510b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the position vector of each NEO as seen from the Sun\n",
    "neo_df.loc[:, \"sun2neo_position_vec\"] = \\\n",
    "    neo_df.apply(lambda x: spiceypy.conics(elts=[x[\"Perihel_km\"],\n",
    "                                                 x[\"Ecc_\"],\n",
    "                                                 x[\"Incl_rad\"],\n",
    "                                                 x[\"LongAscNode_rad\"],\n",
    "                                                 x[\"ArgP_rad\"],\n",
    "                                                 x[\"MeanAnom_rad\"],\n",
    "                                                 x[\"Epoch_et\"],\n",
    "                                                 gm_sun],\n",
    "                                         et=curr_time_et)[:3],\n",
    "                 axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "78208450-dfb0-4597-9217-a8488908d6b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.66166744e+08,  5.62640353e+07, -2.01934353e+07])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's take the position value of (433) Eros and compare it with NASA Horizons\n",
    "neo_df.iloc[0][\"sun2neo_position_vec\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7762bb58-6271-4ce7-948c-fd64d31840df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To compute the apparent magnitude we need to re-compute the positional vectors and convert it to\n",
    "# AU\n",
    "neo_df.loc[:, \"neo2earth_position_vec\"] = \\\n",
    "    neo_df[\"sun2neo_position_vec\"].apply(lambda x: sun2earth_position_vec - x)\n",
    "neo_df.loc[:, \"neo2sun_position_vec\"] = \\\n",
    "    neo_df[\"sun2neo_position_vec\"].apply(lambda x: -1.0 * x)\n",
    "\n",
    "neo_df.loc[:, \"neo2earth_position_vec_AU\"] = \\\n",
    "    neo_df[\"neo2earth_position_vec\"].apply(lambda x: [spiceypy.convrt(k, \"km\", \"AU\") for k in x])\n",
    "neo_df.loc[:, \"neo2sun_position_vec_AU\"] = \\\n",
    "    neo_df[\"neo2sun_position_vec\"].apply(lambda x: [spiceypy.convrt(k, \"km\", \"AU\") for k in x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "342be7b5-0044-48df-97c9-10acb6256207",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the apparent magnitude of each NEO for today!\n",
    "neo_df.loc[:, \"app_mag\"] = \\\n",
    "    neo_df.apply(lambda x: photometry.hg_app_mag(abs_mag=x[\"AbsMag_\"],\n",
    "                                                 vec_obj2obs=x[\"neo2earth_position_vec_AU\"],\n",
    "                                                 vec_obj2ill=x[\"neo2sun_position_vec_AU\"],\n",
    "                                                 slope_g=x[\"SlopeParamG_\"]), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c9538541-7ba3-4d17-8800-a976317ec06f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Name</th>\n",
       "      <th>app_mag</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>433</td>\n",
       "      <td>13.290833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>719</td>\n",
       "      <td>19.527435</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>887</td>\n",
       "      <td>20.360448</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1036</td>\n",
       "      <td>16.235524</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1221</td>\n",
       "      <td>19.483616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29044</th>\n",
       "      <td>2022KY2</td>\n",
       "      <td>27.134045</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29045</th>\n",
       "      <td>2022KZ</td>\n",
       "      <td>27.919159</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29046</th>\n",
       "      <td>2022KZ1</td>\n",
       "      <td>28.624954</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29047</th>\n",
       "      <td>2022KZ2</td>\n",
       "      <td>29.567487</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29048</th>\n",
       "      <td>6344P-L</td>\n",
       "      <td>26.049388</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>29049 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          Name    app_mag\n",
       "0          433  13.290833\n",
       "1          719  19.527435\n",
       "2          887  20.360448\n",
       "3         1036  16.235524\n",
       "4         1221  19.483616\n",
       "...        ...        ...\n",
       "29044  2022KY2  27.134045\n",
       "29045   2022KZ  27.919159\n",
       "29046  2022KZ1  28.624954\n",
       "29047  2022KZ2  29.567487\n",
       "29048  6344P-L  26.049388\n",
       "\n",
       "[29049 rows x 2 columns]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's take a look. We can use e.g., Stellarium to get a \"round about\" feeling whether the\n",
    "# app. mag. of Eros is correct\n",
    "neo_df[[\"Name\", \"app_mag\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ea6a013a-4993-48b7-b096-99736b3689b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Name</th>\n",
       "      <th>AbsMag_</th>\n",
       "      <th>app_mag</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>7444</th>\n",
       "      <td>2010LF64</td>\n",
       "      <td>-9.99</td>\n",
       "      <td>-7.404241</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7437</th>\n",
       "      <td>2010LB67</td>\n",
       "      <td>0.17</td>\n",
       "      <td>4.042399</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>433</td>\n",
       "      <td>10.87</td>\n",
       "      <td>13.290833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>308</th>\n",
       "      <td>65803</td>\n",
       "      <td>18.07</td>\n",
       "      <td>14.523648</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>2100</td>\n",
       "      <td>16.17</td>\n",
       "      <td>15.214294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28673</th>\n",
       "      <td>2022FZ1</td>\n",
       "      <td>25.58</td>\n",
       "      <td>48.914466</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28902</th>\n",
       "      <td>2022JA1</td>\n",
       "      <td>25.38</td>\n",
       "      <td>50.279766</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2407</th>\n",
       "      <td>488698</td>\n",
       "      <td>19.45</td>\n",
       "      <td>72.395760</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7504</th>\n",
       "      <td>2010OJ129</td>\n",
       "      <td>99.90</td>\n",
       "      <td>101.461909</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7360</th>\n",
       "      <td>2010HA71</td>\n",
       "      <td>99.90</td>\n",
       "      <td>104.214032</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>29049 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            Name  AbsMag_     app_mag\n",
       "7444    2010LF64    -9.99   -7.404241\n",
       "7437    2010LB67     0.17    4.042399\n",
       "0            433    10.87   13.290833\n",
       "308        65803    18.07   14.523648\n",
       "25          2100    16.17   15.214294\n",
       "...          ...      ...         ...\n",
       "28673    2022FZ1    25.58   48.914466\n",
       "28902    2022JA1    25.38   50.279766\n",
       "2407      488698    19.45   72.395760\n",
       "7504   2010OJ129    99.90  101.461909\n",
       "7360    2010HA71    99.90  104.214032\n",
       "\n",
       "[29049 rows x 3 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Anyway ... we see that some app. mag. values appear to be outliers! The AbsMag_ does not appear\n",
    "# to be realisitc (a simple placeholder?)\n",
    "# https://newton.spacedys.com/neodys/index.php?pc=1.0\n",
    "neo_df[[\"Name\", \"AbsMag_\", \"app_mag\"]].sort_values(by=\"app_mag\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "38fda1fd-a312-4f4e-af00-0ed246811c55",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We set some NEOs we'd like to remove\n",
    "neos_2_del = [\"2010LF64\", \"2010LB67\", \"2010OJ129\", \"2010HA71\"]\n",
    "neo_df = neo_df.loc[~neo_df[\"Name\"].isin(neos_2_del)].copy()\n",
    "neo_df.reset_index(drop=True, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "fc56d031-2896-4434-ae26-9fb376624f46",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(24.5, 3800, 'Pan-STARRS Lim. Mag.')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvcAAAGBCAYAAADi7J76AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACAsUlEQVR4nO3deVxU9eL/8TcM+74oqwrc1BRUjCVBITXFtO1e019evamVaXuWlpnZNW9lXk0zy7pqXbey7FuZS4tLaaKggQYqLqChsYOyDjCs5/cHMNdJ0fmo8JnDvJ+Px3ncK3OY85nXmTl9HM+csQCggIiIiIiIVM9S9gCIiIiIiOjm4OSeiIiIiKiD4OSeiIiIiKiD4OSeiIiIiKiD4OSeiIiIiKiD4OSeiIiIiKiD4OSeiIiIiKiDkDa5nz17NhRFwfvvv2/w83nz5iEnJwdVVVXYs2cPgoODDW53c3PD+vXrUVpaitLSUqxfvx6urq4G6/Tp0wd79+5FVVUVsrOz8dprr7X54yEiIiIiMgVKey8DBgxQfv/9dyUlJUV5//339T+fNWuWUl5erjzwwANKSEiIsmnTJiUnJ0dxcnLSr/P9998rx48fV6KiopSoqCjl+PHjytatW/W3Ozs7K3l5ecqmTZuUkJAQZcyYMUp5ebkyY8aMdn+cXLhw4cKFCxcuXLi089K+G3RxcVHOnDmjDBkyRNmzZ4/B5D43N1eZM2eO/s92dnZKeXm5Mm3aNAWA0qtXL0VRFGXgwIH6dQYNGqQoiqL07NlTAaA88cQTSllZmWJnZ6df59VXX1Wys7Nlh+bChQsXLly4cOHCpU2Xdj8tZ9WqVfjqq6+wd+9eg58HBQXB19cXO3fu1P9Mp9Nh3759GDhwIAAgOjoaFRUVSEhI0K9z4MABaLVag3Xi4+Oh0+n06+zYsQP+/v4IDAxsuwdGRERERCRZu07uH3vsMXTv3h1z58697DYfHx8AQEFBgcHPCwoK9Lf5+PigqKjost8tLCw0WOdK93HpNoiIiIiIOiKr9tpQz549sWDBAsTExKC+vr69NntNU6dOxbRp0wAAvXv3RlZWFgCgtrYWDQ0NsLe3BwDU19dDp9PByckJAKAoCiorK2Fvbw+NRgMAqKqqgpWVFWxsbAAANTU1UBQFdnZ2V70PBwcHWFo2/T2rsrISNjY2sLa2vuJ91NXVoba2Fo6OjgCAxsZGVFVVXfU+dDodLCwsYGtre9X7cHR0hIWFBQBAq9XCzs4OVlZWV7yP2tpa1NfXw8HBAQDQ0NCA6urqq95HdXU1NBqNvk9r99HSp+U+rtaY+4n7iftJ7n6q1ulgZWUF6+ZxcT+Z5n7i64n7ifvJ9PZTt27d4OXlhbbQLuf/TJ48WVEURamrq9MviqIoDQ0NSl1dndKzZ09FURQlIiLC4Pe2b9+urF27VgGgPPLII0p5efll911RUaE8/PDDCgBl3bp1yvbt2w1uj4iIUBRFUQIDA686xqSkJOnnSbX3EhwcLH0MalrYi83Yy7QW9mIv9jKdhb3Elraad7bbaTnffvst+vTpg/79++uXpKQkfPHFF+jfvz/S09ORl5eHuLg4/e/Y2toiNjZWf459YmIinJ2dER0drV8nOjoaTk5OBuvExsbq/2YIAHFxccjJycG5c+fa58GqiKenp+whqAp7iWMzMWrrddf4ybhr/GRp21dbL9nYSwx7iWEv09Buk/uysjKkpaUZLJWVlSguLkZaWhoAYNmyZXj55ZcxevRohISEYO3atdBqtdi4cSMA4NSpU/jhhx+wcuVKREVFISoqCitXrsS2bduQnp4OANi4cSOqqqqwdu1ahISEYPTo0Zg9ezaWLl3aXg+ViMhsDBx1HwaOuk/2MIiI6BLS/jniz5fCBKDMmzdPyc3NVaqrq5W9e/cqISEhBre7ubkpGzZsUMrKypSysjJlw4YNiqurq8E6ffr0UX755Relurpayc3NVf75z39K/ecRU17+3I4Le7EZe4ks89d/pcxf/xV7qWRhL/ZiL9NZ2nDeKf/BmcpijpP7Ll26SB+Dmhb2YjP2MlxkT+7V1kv20lqv+eu/Uh577S3p4zO1hc8v9mrLpa3mne12tRwyTUFBQcjOzpY9DNVgL3FsJsYcez3z9rsYOnocAKC+rg4X83NxcNcP2PT+YtRUV1/1d2+kl6WlJe5/9EkMHf0gOvt1QV1tDfKzzuOXLV/h+w2f4KOfDsHLv2urv3/81wTMmzQWAPDoq29g5ISHser1l7H7/zYarDd09IN45u1l+j+XXijCmWMp+HTJW8g6k67/+aUdGurrUVxYgCO/7MZn7y5EZXmZfr2AW4Px9+deQo9+t8HRxQVlFy/i7PEUrF34Oopyc6441gefmYnoEffgm8Xzrthr8bOPob6+7trRblBn/y74z0+/oqGhAU8Nj8KFvP+N19HFFav3HYGtnT1mjR2Js8ePtvl4rsUcX483gr1MAyf3REQkXeqBfVj+8rPQWFmjd8QAPPXGO7Czt8eq+a+02TYffHom7powGR+/8SrOHP0Ntg4OCOrdB539ugAAXh47CpbNl74L7BWM1z7+HLPGjsLF/FwATX8RAQAraxvccd9ofLv6AwwbO+GyyT0A6Kqq8PSIaMDCAp7ePpj44lzMWbkBz46M0d/PpR0sNVbo2r0nnn5rCRxdXPHuzKcAAC7uHnh97ZdI2b8XC56YiIrSEnT264LwwcNh7+R83S20ZaXX/bvXo7ggH0MfGIf/W/G/z8Pdcd8DKL94EZ39u7TrWIg6Gk7uzdz58+dlD0FV2Escm4lRW6+Wd65vVF1dDUovNH1J4f7tm9Hn9oG4fdhIrJr/Cu647wHcM+kx+P+lO2p1OqQlJWLNgnkoLszH+fPnEXJ7NP61/mu8/vCDmPDCbAT07I2ss+n4zz9nIfPEsVa3GXHnCOz8Yj0OfL9F/7M/0k/p/395SbH+/1eUljT/7KJ+nC2iRtyNwpxsfPWf5bh74hR07XErsjJOG6yjKIr+90qLCrF93Wq88tE6+AXdYrDNSzsUF+ThwA9bMXT0g/rbe4XdDkcXV6yY84L+LwVFOdk4kXTQiMqtP7/mr/8KWRmn8fEbrwIAPvrpEH7++gt4d+2GqBH3oLK8HOsW/Qsp+/fi8df/jfAhw1FaVIDVb7yK1AO/GLXtS+3d/CWGjn7QYHJ/55jx2PPtl3jw6RkG6z40Yw5ujxuFTr5+KLt4AQk/bMMXyxejrrZGv87oac/g3klTYWvvgEO7vkd+1nnc+cA4PDlsgPDYLqW216Ns7GUa2vUbasn0lJSUyB6CqrCXODYTw15Namt00DR/gY6VtQ02vf8OZv51OBY8MQku7h54YcmHAAx7/WPGK/hs6QK89MBd0JaW4PnFH1x1G6UXChFy+0C4ena6obEOGzse+7Z+jVpdNQ7u/A7Dx0646voOzi6IvXc0AKChrvUvdfTu0g23xQxF/SXrlFwohEajQdRd917XWEWeX/dMegwZR1Pw0gN3IeHHrXh24TK88M4KHNn3E14cHYcTyYcwfdH7sLaxvfad/cnhfT/BxsYWfaNiAABBvfvAp1sAEn7Yetm6uuoqrJgzA9PvGYzV81/BoLv/ijFPPKe/fdDdf8WDT8/AxmX/xktj7kL27xm47+FpwmO6Er4exbCXaeDk3sz1799f9hBUhb3EsZkYtfW6/9EncP+jT9zU++zetz9i7x2NY4n7AQA/f/MFjuz7GQXZf+DMsRSsev0VBEdGwcPb16DXF8sX4fihBORknsGXK95Fl1t6wMPbt9XtrF34Opxc3fBxfAqWbduDJ994BwPiRgmN1cu/K3qHD8D+774FAOz99ivccf8YWFnbGKxn7+iITw9n4LMjZ7Ah6RRi7vkbfv1pB3Iyzxisd1vMUHx6OAMbU87iw90H0bXHrfj24xX62zNSj+Dr/7yHZ99ehnWHTuC1Tz7HA48/i85+/kaNV+T5lXrgF+z4fB3yzmdi0/vvwMbWDnl/nMMvW75C/h/n8H8fvQtXz07o1vNWo++zRWNDPfZu+Qp3jvk7gKa/ICX8sA266qrL1v3qo2U4/VsSinKycWTfz/hm5XLE3PM3/e33TJqCPZu/xE9fbUTeud+xedUHyDj6m/CYrkRtr0fZ2Ms08LQcIiK6buFDhgMAtv73Pzd0Py2TWo2VBhorayT9tAMfv9l0ikhQcF88+PQMBPUKgZObG4Cmr6b/84T23OmT+v9fUpgPAHD17ITigjx8ejhDf9u+bV9j1euzkX02Ay/cNxR/CemH3uG3IzgiCjPfXYnUA79gwROToCjKNcc9bOx4HDu4X38qTdqvCaiprsbtw0cavAutq6rCi6PjoNFYIThyAO5/9AmsnDfrsvs7kXwQ//nnLNjY2WH4//sHfLoG4PsNnxiss3HZv7F17Sr0HTAIPfuHYdiY8Rjz+HQsfOphHDu4/5pjNta50ycMxq+rqsIf6f9rXHbhAgDA1eP6/uXj52++wDvf7IRbp86IuedvWPDEpCuuF3XXPbh30lT4dAuEnYMjLDWWsLTU6G/3D+qO3V8afs4hI/UI/AL/cl3jIlI7Tu7NXGlpqewhqAp7iWMzMebaq2VS21Bfh+LCAjTUN52KYmtvj9c+3ohjifFY/vKzKLt4Ec7uHnhr47ewsrZu6mXvBMDwFBcFTRNzS8umvwi8OPp/335epa3433qKgrPHU3H2eCq2r1uNO+57ANMXf4DgyCik/Zp41TFbWlpi6OgH4e7lgy+P/6H/uYWlJYaPnWAwuVcUBfl/nAMA5GSegVtnb7yw5EPMm/z/DO6zRletX++/b72G+ev+D2OfegFffrDEYD1taQkSd2xH4o7t+Gzp23jnm50Y+9Tz15zcizy/WvbBpY/h0p+1/OXHwvL6TgLIzTyL308cwwtLPkTphSKkpxy+7MO0PULDMGPJR/hyxVL8Fr8HVRXliLhzBB5+ed51bVOUub4erxd7mQZO7s3csWOtf9iMLsde4thMjLn2unRSeyn/v3SHq4cnPlv6NgpzsgAAA7r/79SZY8eOIeT26Gve/5Xu+0qyzja9w2/n4HjNdfvHDoWTmztmjR2F+rpa/c87+/rjlf+sR2f/LijKufJlAbevXYX7Hp6GAXGjcGjXD61u48sVS/Hqqk+x68tPUVJYcMV16uvqkJ91Hu5e3tccs6k9v376+nM8s+BdrFv0ryve3issEsUF+fjqo2X6n7VczahFTuYZdO8bip+/+UL/sx79brsp4zO1XqaOvUwDz7k3czExMbKHoCrsJY7NxLCXoQu5Oait0WHUQ4/Au0s3hA0ehvHP/e90lhvp9eJ7q3Dv5Kno0e82dPbzR8jt0Zj6z7dQUlSI078lX/P3h4+dgN/27UHmiWPIyjitX47s+xm5mWcx7IHxrf5udaUWP321EeOefREWFhatrpf2ayKyz2Zg7BPTATSdBvXcovcRPmQ4fAP/Ar+gW3D/o08g7I478etV/pIAANa2thg9/h8I7BWiX3xv8qkry7/fh1H/eMTo9X/59v/wcFQfbF+3+oq35537HR7ePoi9dzS8u3TDXX+fZHC+PQB8t/4TDBn9IO584O/wDQjCX6c8hR79bjM4rer24SOx/Pt98PDyEXo8fD2KYS/TwHfuzdzV/qNCl2MvcWwmRm29anW6Nr3/8pJivD/7efzjhdkYOeFhnD99Emv//Tpe+/hzADfWK2X/Xgwa9VeMnvqM/sugTh1JwkdzX7zmdd9dPTshfMhwvP/K81e8PfHH7Rj6wDhs+uCdVu/ju/Wf4O6HpmDQPX/D/u2bW11v65qVeGbBUnz78Qpkn0mHrqoSk2f9E54+fmhoqEdhdhbWLfoXvlv/8VXH7BsQhIfmLcZDl/zszPFUvDxW7EPEV+P/l+5wdvcwev3GxkZUlBa3envynl3Y8t+P8Micf8HG1g6pB37BpuWLMe31hfp1Dny/Bd5du+GhmXNgY2ePQ7u+x84vNiBy2F36dRydXeD/l+7QWItNe9T2epSNvUyH9K/fNZWlrb4G2JSXmJgY6WNQ08JebMZeprWwF3tdaZn1/ifKKx+tY692XthLbGmreadF8/8hAElJSYiMjJQ9DCIiIjKSjZ097ho/CSnxe9DQ0ICoEXdj/PSXsfi5x/Dr7h9lD4+oVW017+Q592aub9++soegKuwljs3EqK3X2Cefx9gnn5e2fbX1kq1D9lIUhMXeiTc+/QbvfLMTg0bdj/dmPXtTJvYdslcbYi/TwHPuzZybm5vsIagKe4ljMzFq69U3uukDdJdezaQ9qa2XbB2xV22NDvMfHdcm990Re7Ul9jINfOeeiIiIiKiD4OTezKWkpMgegqqwlzg2E8NeYthLDHuJYS8x7GUaOLk3c+7u7rKHoCrsJY7NxLCXGPYSw15i2EsMe5kGTu7NXEBAgOwhqAp7iWMzMWrrpS0tgba0RNr21dZLNvYSw15i2Ms08AO1RER03RY/N1X2EIiI6BJ8597MZWZmyh6CqrCXODYTw15i2EsMe4lhLzHsZRo4uTdzFRUVsoegKuwljs3EqK3XP2a8gn/MeEXa9tXWSzb2EsNeYtjLNHByb+b69esnewiqwl7i2EyM2nr17B+Onv3DpW1fbb1kYy8x7CWGvUwDJ/dERERERB0EJ/dm7uLFi7KHoCrsJY7NxLCXGPYSw15i2EsMe5kGXi3HzJ04cUL2EFSFvcSxmRj2EmOOvbZuvZHfvv5e999/I9tVJ3N8ft0I9jINfOfezMXGxsoegqqwlzg2E6O2Xhfz83AxP0/a9tXWSzYPD/YSweeXGPYyDXznnoiIrtvyWc/KHgIREV2C79ybufr6etlDUBX2EsdmYthLDHuJURT2EsHnlxj2Mg0WABTZgzAVSUlJiIyMlD0MIiLVeOSV+QCANW/PkzwS83Fj59xfP3M8556oLbXVvJPv3Ju50NBQ2UNQFfYSx2Zi1NYrsHcIAnuHSNu+2nrJ5uLCXiL4/BLDXqah3Sb3Tz31FFJTU1FWVoaysjIkJCTg7rvv1t++Zs0aKIpisCQmJhrch42NDZYvX46ioiJotVps2bIF/v7+But07doVW7duhVarRVFREd577z1YW1u3y2NUIxcXF9lDUBX2EsdmYthLDHuJsbJiLxF8folhL9PQbpP77OxsvPzyywgLC0NERAR+/vlnfPvtt+jbt69+nV27dsHHx0e/XDr5B4Bly5ZhzJgxGD9+PGJjY+Hi4oLt27fD0rLpYVhaWuK7776Ds7MzYmNjMX78eIwdOxZLlixpr4dJRERERCRNu10tZ+ufThKcO3cunnzySURHR+PYsWMAgJqaGhQUFFzx911cXDBlyhQ88sgj2L17NwBg4sSJOH/+PIYPH46dO3dixIgRCAkJQUBAALKzswEAs2bNwscff4xXX30VFRUVbfgI1enw4cOyh6Aq7CWOzcSwlxj2ElNWxl4i+PwSw16mQco595aWlhg3bhycnJyQkJCg/3lMTAwKCgpw+vRprFq1Cp07d9bfFh4eDhsbG+zcuVP/s+zsbJw8eRIDBw4EAERHR+PkyZP6iT0A7NixA3Z2dggPD2+HR6Y+Xl5esoegKuwljs3EqK1X3rnfkXfud2nbV1sv2Wxs2EsEn19i2Ms0tOt17vv06YPExETY2dlBq9Vi9OjROH78OADgxx9/xDfffIPMzEwEBgbizTffxM8//4zw8HDU1tbCx8cH9fX1uHDhgsF9FhQUwMfHBwDg4+Nz2Tv/Fy5cQH19vX6dP5s6dSqmTZum//2WL2DIzMxERUUF+vXrB6DpK5VPnDihv72+vh6JiYkIDQ3Vn2N2+PBheHl5oWvXrgCAM2fOoKamBiEhTR82KyoqQkZGhv4vI7W1tTh06BDCwsLg6OgIoOmT0/7+/vDz8wMApKeno6GhAb179wYA5Ofn49y5c4iKigIAVFdXIzk5GREREbC3twcAHDx4EIGBgfrHfPLkSWg0GvTs2RMAkJubi5ycHERGRiI0NBQJCQk4cuQIBgwYABsbGwBAQkICevToof8LVlpaGmxtbdG9e3cAQFZWFgoLC/V/aSovL0dqaiqio6NhZdX0tIqPj0dwcDA8PT0BAEePHoWzszOCgoIAAOfPn0dJSQn69+8PACgtLcWxY8cQExMDCwsLKIqC/fv3o2/fvnBzcwMApKSkwN3dHQEBAVL208CBA5Gamtru+wkAKisrVbmfQkND8e2335rF6+lm7KfQ0FDEx8er5vW0fuHr8PLy0q/T3vspKioKqampZvN6AgA7u0zU11fAxaVfc9OL0GpP6L+gSlHqUVKSCBeXUP059mVlh2Fj4wVv73tRXp6KysozaGysgbNzSPN9FKGyMgPu7k37qbGxFqWlh+DqGgaNxrF5u+p7Pd3ofrr33nuRmpqqiv8+mcJx76677tJ/S61aXk8y91NbaddLYVpbW6Nbt25wdXXF2LFjMXXqVAwZMgRpaWmXrevr64vz589j3Lhx2Lx5M8aPH4/169df9uHYn376CRkZGXjiiSewcuVKdO/eHcOGDTNYp66uDhMnTsQXX3xx1fGZ46UwY2NjER8fL3sYqsFe4thMDHuJMcdeN3IpTA+PWBQXX18vc7wUpjk+v24Ee4npEJfCrKurw9mzZ3HkyBHMmTMHKSkpeOGFF664bl5eHrKzs9GjRw8ATX/TtLKyQqdOnQzW8/b2Rn5+vn4db29vg9s7deoEKysr/Tpk6MyZM7KHoCrsJY7NxKit1xP/WoQn/rVI2vbV1ku2ykr2EsHnlxj2Mg1Sr3NvaWkJW1vbK97m6ekJf39/5OXlAWj6J5Da2lrExcXp1/H390fv3r315+0nJiaid+/eBpfHjIuLg06n44c8WlFTUyN7CKrCXuLYTIzaevkG/gW+gX+Rtn219ZKtsZG9RPD5JYa9TEO7Te7ffvttxMTEICAgAH369MGCBQswZMgQfPbZZ3B0dMTixYsRFRWFgIAADB48GNu2bUNhYSE2b94MoOlcrE8++QSLFi3CsGHD0L9/f2zYsAFHjx7VXz1n586dSEtLw/r169G/f38MGzYMixcvxurVq3mlnFa0nB9GxmEvcWwmhr3EsJeYlnPsyTh8folhL9PQbh+o9fHxwaeffgofHx+UlZXh6NGjGDVqFHbu3Ak7Ozv07dsXkyZNgpubG/Ly8rBnzx48+OCD0Gq1+vt4/vnnUV9fj02bNsHe3h4//fQTJk2ahMbGRgBAY2Mj7rnnHnz44Yc4cOAAqqur8dlnn+Gll15qr4dJRERERCRNu03uH3nkkVZv0+l0GDly5DXvo7a2Fs899xyee+65VtfJysrCfffdd11jNEdFRUWyh6Aq7CWOzcSwlxj2ElNby14i+PwSw16moV0vhUmmJyMjQ/YQVIW9xLGZGLX1Onfy8qudtSe19ZKtspK9RPD5JYa9TAMn92Zu4MCBvGyVAPYSx2Zi1NZrzdvzpG5fbb1kc3cfeN2XwryRS3DeCJmX4OTzSwx7mQapV8shIiIiIqKbh5N7M1dbWyt7CKrCXuLYTIzaej236H08t+h9adtXWy/ZGhvZSwSfX2LYyzTwtBwzd+jQIdlDUBX2EsdmYtTWy9PHV+r21dZLttJS9hLB55cY9jINfOfezIWFhckegqqwlzg2E8NeYthLjKsre4ng80sMe5kGTu7NnKOjo+whqAp7iWMzMewlhr3EaDTsJYLPLzHsZRp4Wg4REZEgWVeOISK6Fk7uzVxSUpLsIagKe4ljMzFq65Wecljq9tXWS7bSUvYSweeXGPYyDTwtx8z5+/vLHoKqsJc4NhOjtl6fLX0bny19W9r21dZLNjs79hLB55cY9jINnNybOT8/P9lDUBX2EsdmYthLDHuJsbNjLxF8folhL9PAyT0REV23l5avxkvLV8seBhERNeM592YuPT1d9hBUhb3EsZkYtfVycnOXun219ZKtspK9RPD5JYa9TAPfuTdzDQ0NsoegKuwljs3EsJcY9hKjKOwlgs8vMexlGji5N3O9e/eWPQRVYS9xbCaGvcSwlxgnJ/YSweeXGPYyDZzcExERERF1EDzn3szl5+fLHoKqsJc4NhOjtl7HEvdL3b7aeslWU8NeIvj8EsNepoGTezN37tw52UNQFfYSx2Zi1Nbrq4+WSd2+2nrJVlV1TvYQVIXPLzHsZRp4Wo6Zi4qKkj0EVWEvcWwmhr3EsJcYd3f2EsHnlxj2Mg2c3BMR0XV7ddWneHXVp7KHQUREzXhajpmrrq6WPQRVYS9xbCZGbb1s7Oykbl9tvWRraGAvEXx+iWEv08B37s1ccnKy7CGoCnuJYzMx7CWGvcSUlbGXCD6/xLCXaeDk3sxFRETIHoKqsJc4NhPDXmLYS4yrK3uJ4PNLDHuZBk7uzZy9vb3sIagKe4ljMzHsJYa9xGg07CWCzy8x7GUaeM49ERFdt8N7d8seAhERXcICgCJ7EKYiKSkJkZGRsofRrqytrVFXVyd7GKrBXuLYTAx7iZHVa+vWdt/kTWFhYQ1FUdfz6/775W2br0cx7CWmreadPC3HzAUGBsoegqqwlzg2E8NeYthLjINDoOwhqAqfX2LYyzRwcm/mfHx8ZA9BVdhLHJuJUVuv+eu/wvz1X0nbvtp6yWZry14i+PwSw16mgZN7IiIiIqIOgpN7M3fy5EnZQ1AV9hLHZmLYSwx7idFq2UsEn19i2Ms0tNvVcp566ik8/vjj+vOx0tLS8Oabb+L777/XrzNv3jxMmzYN7u7uOHToEJ5++mmcOHFCf7ubmxuWL1+O+5s/XbN161Y8++yzKCsr06/Tp08ffPDBB7j99ttRXFyMlStX4o033mifB6lCGo1G9hBUhb3EsZkY9hLz2Wca1NTIHoV6WFjw+SWCr0cx7GUa2u2d++zsbLz88ssICwtDREQEfv75Z3z77bfo27cvAGDWrFmYOXMmnn32WURGRqKwsBC7du2Ck5OT/j42btyIsLAwjBw5EiNHjkRYWBg2bNigv93Z2Rm7du1CQUEBIiMjMX36dLz00kuYMWNGez1M1enZs6fsIagKe4ljMzHsJcbRkb1EsJcYvh7FsJdpaLd37rf+6bphc+fOxZNPPono6GgcO3YMzz//PBYuXIhvvvkGADB58mQUFhZiwoQJWLVqFXr16oVRo0Zh0KBBOHjwIADg8ccfx/79+9GzZ0+kp6fjH//4BxwcHDB58mTodDqkpaWhV69emDFjBpYuXdpeD5WIyGwk/LBN9hCIiOgSUs65t7S0xLhx4+Dk5ISEhAQEBQXB19cXO3fu1K+j0+mwb98+DBw4EAAQHR2NiooKJCQk6Nc5cOAAtFqtwTrx8fHQ6XT6dXbs2AF/f39enqkVubm5soegKuwljs3EqK3Xjs/XYcfn66RtX6dTVy/Z2EuM2l6PsrGXaWjXb6jt06cPEhMTYWdnB61Wi9GjR+P48eOIjo4GABQUFBisX1BQAH9/fwBNl1cqKiq67D4LCwv1l17y8fFBdnb2ZffRctu5c+cu+/2pU6di2rRp+nViY2MBAJmZmaioqEC/fv0AABcvXsSJEyf0t9fX1yMxMRGhoaFwcXEBABw+fBheXl7o2rUrAODMmTOoqalBSEgIAKCoqAgZGRn6v4zU1tbi0KFDCAsLg6OjI4CmLzTw9/eHn58fACA9PR0NDQ3o3bs3ACA/Px/nzp1DVFQUAKC6uhrJycmIiIjQf+3zwYMHERgYqO9y8uRJaDQa/T+X5ebmIicnB5GRkbCxsYGrqyuOHDmCAQMGwMbGBgCQkJCAHj16oHPnzgCaPiNha2uL7t27AwCysrJQWFiI8PBwAEB5eTlSU1MRHR0NK6ump1V8fDyCg4Ph6ekJADh69CicnZ0RFBQEADh//jxKSkrQv39/AEBpaSmOHTuGmJgYWFhYQFEU7N+/H3379oWbmxsAICUlBe7u7ggICJCyn9zd3eHn59fu+wkAKisrVbmfbGxsUFNTYxavp5uxn2xsbODk5KSa19OxtDR07twZPl5eUvaTg4Mr7Oz8UFJyEA4OgfpLPWq1J2FhodGfhqLT5UKny4GbW9N+amioRFnZEbi5DYClZdN+KilJgKNjD9jYNO2nioo0WFrawtGxe/N2s1BbWwhX1/DmHuUoL0+Fu3s0LCya9lNxcTycnIJhY9O0n8rLj8LKyhkODkHN93EedXUlcHFp2k91daWoqDgGD48YtHyvZHHxfjg794W1tVvzfaTA2tod9vZN+6mqKhP19RVwcenX3PQitNoT8PBo2k+KUo+SkkS4uITCyqppP5WVHYaNjRccHAJhZ+eHysozaGysgbNzSPN9FKGyMgPu7k37qbGxFqWlh+DqGgaNxrH5OZUEOzt/2Nn5NT/X06EoDXByatpPNTX5qKo6B3f3qObG1SgrS4arawQ0GvvmxuL7CZB33AsMDISfn58q/vtkCse9S7erhv8+yd5PbaVdv6HW2toa3bp1g6urK8aOHYupU6diyJAhcHFxQUJCArp164asrCz9+p988gn8/f0xcuRIvPLKK3jsscdwyy23GNzn2bNnsXr1aixcuBA7duxAdnY2pkyZor+9a9eu+OOPPxAdHa0/nac15vgNtbGxsYiPj5c9DNVgL3FsJkZtvVqucT9v0lgp29+/PxbFxerpJZuHh/p6yfyGWrW9HmVjLzEd4htq6+rqcPbsWRw5cgRz5sxBSkoKXnjhBeTn5wMAvL29Ddb39vbW35afn6//29+lvLy8DNa50n203EZERERE1JFJvc69paUlbG1tkZmZiby8PMTFxelvs7W1RWxsrP4c+8TERDg7O+tP4QGazrFvOW+/ZZ3Y2FjY2trq14mLi0NOTs4VT8mhpn9KI+Oxlzg2E8NeYppO2yBjsZcYvh7FsJdpaLfJ/dtvv42YmBgEBASgT58+WLBgAYYMGYLPPvsMALBs2TK8/PLLGD16NEJCQrB27VpotVps3LgRAHDq1Cn88MMPWLlyJaKiohAVFYWVK1di27ZtSE9PB9B0qcyqqiqsXbsWISEhGD16NGbPns0r5VzFkSNHZA9BVdhLHJuJYS8xZWXsJYK9xPD1KIa9TEO7Te59fHzw6aef4vTp0/jpp58QGRmJUaNG4ccffwQALFq0CO+++y5WrFiB5ORk+Pr6YsSIEdBqtfr7mDBhAlJTU7Fjxw7s2LEDqampmDhxov728vJyxMXFwc/PD8nJyVixYgWWLFnCyf1VDBgwQPYQVIW9xLGZGPYS4+bGXiLYSwxfj2LYyzS029VyHnnkkWuuM3/+fMyfP7/V20tLSw0m81dy/PhxDB48WHh85qrl6gNkHPYSx2Zi1NZr7+YvpW6/5Uo3ZBz2EqO216Ns7GUa2vVSmERE1LHskTy5JyIiQ+16KUxTZ46XwtRoNGhoaJA9DNVgL3FsJkZtvZzdPAAAFaXFUra/bZsGiqKeXrJZWKivl8xLYart9Sgbe4npEJfCJNPTo0cP2UNQFfYSx2Zi1NbrxeWr8OLyVdK27+iorl6ysZcYtb0eZWMv08DJvZm70ncHUOvYSxybiWEvMS3fJkvGYS8xfD2KYS/TwMk9EREREVEHwcm9mUtLS5M9BFVhL3FsJoa9xFRUsJcI9hLD16MY9jINvFqOmbv023zp2thLHJuJYS8xlpbsJUKNvbZulbftxx9XXy+ZePwyDXzn3sx1795d9hBUhb3EsZkYtfXa8fl67Ph8vbTtOzqqq5ds7CVGba9H2djLNPCdeyIium4JP0h8W5WIiC7Dd+7NXFZWluwhqAp7iWMzMWrr5enjB08fP2nbr65WVy/Z2EuM2l6PsrGXaeA792ausLBQ9hBUhb3EsZkYtfV6btFyAMC8SWOlbL+2Vl29ZGMvMWp7PcrGXqaB79ybufDwcNlDUBX2EsdmYthLjKsre4lgLzF8PYphL9PAyT0RERERUQfByb2ZKy8vlz0EVWEvcWwmhr3E1Nezlwj2EsPXoxj2Mg2c3Ju51NRU2UNQFfYSx2Zi2EtMeTl7iWAvMXw9imEv08DJvZmLjo6WPQRVYS9xbCZGbb22rVmJbWtWStu+u7u6esnGXmLU9nqUjb1MA6+WY+asrPgUEMFe4thMjNp6Je/ZJXX7Fhbq6iUbe4lR2+tRNvYyDXznnoiIrptf0C3wC7pF9jCIiKiZBQDlWivdcccd0Ol0+PXXXwEAkydPxmOPPYa0tDTMnDkTlZWVbT3OdpGUlITIyEjZwyAiUo35678CIO8691v5BbnUhu6/X/YIqCNrq3mnUe/cL1u2DD4+PgCAnj17YuXKlTh69Ciio6OxePHimz4oaj/BwcGyh6Aq7CWOzcSwlxgnJ/YSwV5i+HoUw16mwajJfffu3XHs2DEAwJgxY7Br1y48/fTTmDp1Ku677742HSC1LU9PT9lDUBX2EsdmYthLjI0Ne4lgLzF8PYphL9Ng1OS+sbERGo0GADBs2DD8+OOPAID8/HzuSCIiIiIiE2HU5D4pKQmvvfYaHnroIcTGxuKHH34AAAQGBiIvL69NB0ht6+jRo7KHoCrsJY7NxLCXmPJy9hLBXmL4ehTDXqbBqMn9888/j9tuuw0ffPAB3nrrLfz+++8AgP/3//4fEhMT23SA1LacnZ1lD0FV2Escm4lRW6+vP3oPX3/0nrTtW1mpq5ds7CVGba9H2djLNBh1QdK0tDT069fvsp+/+OKLaGhouOmDovYTFBSE7Oxs2cNQDfYSx2Zi1NbraGK81O07OARBp1NPL9nYS4zaXo+ysZdpMPrbBvr06YNbb70VAJCeno5jx46hpqamzQZGRESmL7BXCN57D7CoSZM9FCIighGT+7CwMKxZswYhISGwsLAAACiKguPHj+ORRx7Bb7/91uaDpLZz/vx52UNQFfYSx2Zi1NbrkTnz0dAZsMqWc5376mp19ZKNvcSo7fUoG3uZhquec9+zZ0/s2bMHdXV1mDhxIsLCwhAWFoZJkyahoaEBe/bsQc+ePdtrrNQGSkpKZA9BVdhLHJuJYS8xdXXsJYK9xPD1KIa9TMNVJ/evv/469u7di4iICHz++edITU1FamoqNm7ciPDwcPzyyy94/fXX22mo1Bb69+8vewiqwl7i2EwMe4lxcekvewiqwl5i+HoUw16m4aqn5dx5551X/ZKqBQsWYMuWLTd9UEREREREJO6q79y7ubkhNze31duzs7Ph6upq1IZmz56NX3/9FWVlZSgsLMTWrVsREhJisM6aNWugKIrB8udLbdrY2GD58uUoKiqCVqvFli1b4O/vb7BO165dsXXrVmi1WhQVFeG9996DtbW1UeM0N6WlpbKHoCrsJY7NxLCXmLq6UtlDUBX2EsPXoxj2Mg1XndxnZ2df8RKYLfr374+cnByjNjRkyBB8+OGHGDhwIO68807U19dj9+7dcHd3N1hv165d8PHx0S933323we3Lli3DmDFjMH78eMTGxsLFxQXbt2+HpWXTQ7G0tMR3330HZ2dnxMbGYvz48Rg7diyWLFli1DjNzbFjx2QPQVXYSxybiVFbr43vLoTlhYXStl9Roa5esrGXGLW9HmVjL9Nw1cn95s2bsXjxYnh5eV12m7e3NxYuXIjNmzcbtaGRI0di7dq1SEtLw/HjxzFx4kR07twZgwYNMlivpqYGBQUF+uXSD2e4uLhgypQpeOmll7B792789ttvmDhxIvr164fhw4cDAEaMGIGQkBBMnDgRv/32G3bv3o1Zs2Zh6tSp/HKFK4iJiZE9BFVhL3FsJkZtvU7/lgxLXbK07Xt4qKuXbOwlRm2vR9nYyzRcdXL/xhtvwNLSEmfOnMGHH36I5557Ds899xz+85//ICMjA1ZWVnjzzTeva8POzs7QaDSXfbI6JiYGBQUFOH36NFatWoXOnTvrbwsPD4eNjQ127typ/1l2djZOnjyJgQMHAgCio6Nx8uRJgy9R2LFjB+zs7BAeHn5dY+3IWi5vSsZhL3FsJkZtvW69LQKNdhESR6CuXvKxlwi1vR5lYy/TcNUP1JaXlyM6OhoLFizAuHHj9KfQlJaWYsOGDXj11VdRVlZ2XRt+77338NtvvxmcU//jjz/im2++QWZmJgIDA/Hmm2/i559/Rnh4OGpra+Hj44P6+npcuHDB4L4KCgrg4+MDAPDx8UFBQYHB7RcuXEB9fb1+nUtNnToV06ZN0/9ubGwsACAzMxMVFRX605IuXryIEydO6G+vr69HYmIiQkND4eLiAgA4fPgwvLy80LVrVwDAmTNnUFNTo/9sQVFRETIyMvR/EamtrcWhQ4cQFhYGR0dHAEBSUhL8/f3h5+cHoOkLwxoaGtC7d28AQH5+Ps6dO4eoqCgAQHV1NZKTkxEREQF7e3sAwMGDBxEYGKh/vCdPnoRGo9FftjQ3Nxc5OTmIjIxEv379UFlZiSNHjmDAgAGwsbEBACQkJKBHjx76v1ylpaXB1tYW3bt3BwBkZWWhsLBQ/xem8vJypKamIjo6GlZWTU+r+Ph4BAcHw9PTEwBw9OhRODs7IygoCEDT9XBLSkr0n64vLS3FsWPHEBMTAwsLCyiKgv3796Nv375wc3MDAKSkpMDd3R0BAQFS9lPLdtp7PwFQ7X7q16+f2byebsZ+6tevn6peT3c9PgPwsoaH9t/Nj/8MGhtr4Owc0rxfilBZmQF396b91NhYi9LSQ3B1DYNG49g81iTY2fnDzs6v+T7SoSgNcHJq2k81NfmoqjoHd/em/dTQUI2ysmS4ukbAxaVprCUlB+HgEAhb26b9pNWehIWFBo6OTftJp8uFTpcDN7fI5vuoRFnZEbi5DYClpU3zfSTA0bEHbGya9lNFRRosLW3h6Ni9+fmRhdraQri6hjf3KEd5eSrc3aNhYdG0n4qL4+HkFAwbm6b9VF5+FFZWznBwCGq+j/OoqyvRX7Wmrq4UFRXHmt9RtwCgoLh4P5yd+8La2q35PlJgbe0Oe/um/VRVlYn6+gr9Y6+tvQit9gQ8PJr2k6LUo6QkES4uobCyatpPZWWHYWPjpf+d9t5PGo29KvdTy2tBDf99MoXjXss2AHX890n2fmorTUcSI7XsoKKiohva6JIlS/D3v/8dMTExyMzMbHU9X19fnD9/HuPGjcPmzZsxfvx4rF+//rIPx/7000/IyMjAE088gZUrV6J79+4YNmyYwTot1+r/4osvWt1eUlKS/olKRETXNn/9V+jbR96XWBG1pfvvlz0C6sjaat551dNyevToYfDnoqIig4m9lZUVhg4dKrTBpUuXYvz48bjzzjuvOrEHgLy8PGRnZ+vHkZ+fDysrK3Tq1MlgPW9vb+Tn5+vX8fb2Nri9U6dOsLKy0q9D/9O3b1/ZQ1AV9hLHZmLYS4yzM3uJYC8xfD2KYS/TcNXJ/cmTJw3OeT9//jy6deum/7OHhwd27dpl9MaWLVumn9ifPn36mut7enrC398feXl5AJr+GaS2thZxcXH6dfz9/dG7d28kJCQAABITE9G7d2+Dy2PGxcVBp9Ph8OHDRo/VXLT8MxUZh73EsZkY9hLTcuoKGYe9xPD1KIa9TMNVz7n/8wcj3N3dodForrpOaz744ANMnDgRf/vb31BSUqJ/d12r1aKyshKOjo54/fXX8fXXXyMvLw+BgYF4++23UVhYqL8iT3l5OT755BMsWrQIhYWFuHjxIpYuXYqjR49i9+7dAICdO3ciLS0N69evx8yZM+Hp6YnFixdj9erVqKioMK4KEREREZEKXfWde2MoinGn7D/99NNwcXHBzz//jPz8fP3y4osvAgAaGhrQt29fbNmyBenp6Vi3bh1Onz6N6OhoaLVa/f08//zz2Lx5MzZt2oQDBw5Aq9XivvvuQ2NjIwCgsbER99xzD6qqqnDgwAFs2rQJX3/9tX47ZCglJUX2EFSFvcSxmRi19VqzYB40RfOkbb+8PEXattWIvcSo7fUoG3uZhqu+c38zXesdfp1Oh5EjR17zfmpra/WX5GxNVlYW7rvvPuExmiN3d3f+i4YA9hLHZmLU1uvcqTRY1MjbvrW1O+rr1dNLNvYSo7bXo2zsZRqu+s69oihwd3fXL4qiwM3NTf9nDw+P9hontZGWS0GRcdhLHJuJUVuvftGxaHSIlbb9lstDknHYS4zaXo+ysZdpuOY59ydOnDD486XX5Wy5higREZmnMU9OR6MHYFkVL3soRESEa0zuRS9zSepzrcuRkiH2EsdmYthLTFUVe4lgLzF8PYphL9Nw1cn9vn372mscJAnPjRPDXuLYTAx7ieH542LYSwxfj2LYyzTc8NVySN1avm6ZjMNe4thMDHuJcXFhLxHsJYavRzHsZRqu+s59eXn5Nc+pb/mQLRERERERyXXVyf0zzzzT6m0hISF45plnYG1tfdMHRe3n4sWLsoegKuwljs3EqK3Xynkv46MP5W2/tlZdvWRjLzFqez3Kxl6mwQKA0OVufHx88MYbb2Dy5Mn48ccfMXv2bIMr6qhZUlISIiMjZQ+DiEhVtm6VPQKitnH//bJHQB1ZW807jT7n3tnZGW+99RYyMjIQEhKCYcOG4f777+8wE3tzFRsr7/rUasRe4thMjNp6RQyNQ6NjnLTte3ioq5ds7CVGba9H2djLNFzzG2qtrKzw9NNPY+7cubh48SIefvhhfP311+0xNiIiMnH3PfI4Gt0By8pdsodCRES4xuR+woQJeOONN2Bvb4+5c+di9erVaGxsbK+xUTuor6+XPQRVYS9xbCaGvcQoCnuJYC8xfD2KYS/TcNVz7hsaGlBdXY1NmzZBq9W2eifTp09vi7G1O55zT0QkZv76r9C3D2CVPVb2UIhuOp5zT22pread1/wSK0VREBQU1Oo617pUJpm20NBQpKamyh6GarCXODYTw15iXFxCUV7OXsZiLzF8PYphL9Nw1cn90KFD22scJImLi4vsIagKe4ljMzHsJcbKir1EsJcYvh7FsJdp4DfUEhHRdVs+6zlo8p+TPQwiImp2zavlUMd2+PBh2UNQFfYSx2Zi1NbrYn4uLCR+hq6sTF29ZGMvMWp7PcrGXqaB79ybOS8vL9lDUBX2EsdmYtTWa+Co+9HoJO9ThzY26uolG3uJUdvrUTb2Mg2c3Ju5rl27yh6CqrCXODYTo7Zed42fhEa3SdK2b2+vrl6ysZcYtb0eZWMv09Dq5P6TTz6Bk5MTgKZvHNNoNO02KCIiIiIiEtfqOfcPPfQQ5syZA61Wiz179sDX1xdFRUXtOTZqB2fOnJE9BFVhL3FsJoa9xFRWspcI9hIzd+4Z1NS0/3bVen19Hr9MQ6uT+3PnzuHZZ5/Fzp07YWFhgejoaJSUlFxx3fj4+DYbILWtGhlHLRVjL3FsJoa9xDQ2spcI9hLDXmJ4/DINrZ6W89JLL2HatGnYs2cPFEXB5s2bsXfv3suWPXv2tOd46SYLCQmRPQRVYS9xbCaGvcQ4O7OXCPYSw15iePwyDa2+c79161Z4eXnB1dUVxcXFCAkJQWFhYXuOjYiITNw7z03Dp5/KHgUREbW45nXuy8rKMHToUGRkZKChoaE9xkTtiJ+jEMNe4thMjNp6VZQWw6JR3vZra9XVSzb2EsNeYtR2/OqojPoSq3379sHGxgaTJk1CcHAwFEXBiRMnsHHjRtTW1rb1GKkNZWRkyB6CqrCXODYTo7ZeQ0c/iEYXwLL8Synbr6xUVy/Z2EsMe4lR2/GrozLqOve9e/dGeno6li5digEDBiAqKgrvvvsu0tPT0atXr7YeI7WhgQMHyh6CqrCXODYTo7ZeQ0Y/iEaXB6Vt391dXb1kYy8x7CVGbcevjsqoyf17772HlJQUdOvWDXfccQfuuOMOdOvWDampqVi2bFkbD5GIiIiIiIxh1Gk5gwYNQmRkJCoqKvQ/q6iowKuvvoqDBw+22eCo7fG0KjHsJY7NxLCXmMZG9hLBXmLYSwyPX6bBqHfudTod3NzcLvu5q6srdDrdzR4TtaNDhw7JHoKqsJc4NhPDXmJKS9lLBHuJYS8xPH6ZBqMm99u2bcPq1asxcOBAWFpawtLSEoMGDcLKlSuxdetWozY0e/Zs/PrrrygrK0NhYSG2bt16xeuhzps3Dzk5OaiqqsKePXsQHBxscLubmxvWr1+P0tJSlJaWYv369XB1dTVYp0+fPti7dy+qqqqQnZ2N1157zagxmqOwsDDZQ1AV9hLHZmLYS4yrK3uJYC8x7CWGxy/TYNTkfvr06cjIyEB8fDx0Oh10Oh1++eUXpKen4/nnnzdqQ0OGDMGHH36IgQMH4s4770R9fT12794Nd3d3/TqzZs3CzJkz8eyzzyIyMhKFhYXYtWsXnJyc9Ots3LgRYWFhGDlyJEaOHImwsDBs2LBBf7uzszN27dqFgoICREZGYvr06XjppZcwY8YMI5OYF0dHR9lDUBX2EsdmYtTW661pE6HJmSht+xqNunrJxl5i2EuM2o5fHZVR59yXlZXhb3/7G2655Rb07t0bAHDy5EmcPXvW6A2NHDnS4M8TJ05EWVkZBg0ahO3btwMAnn/+eSxcuBDffPMNAGDy5MkoLCzEhAkTsGrVKvTq1QujRo3CoEGD9Of6P/7449i/fz969uyJ9PR0/OMf/4CDgwMmT54MnU6HtLQ09OrVCzNmzMDSpUuNHi8REV1bra4aForsURARUQuj3rlvcfbsWWzfvh3bt28XmthfibOzMzQaDUpKSgAAQUFB8PX1xc6dO/Xr6HQ67Nu3T39ppejoaFRUVCAhIUG/zoEDB6DVag3WafkXhhY7duyAv78/AgMDb2jMHVFSUpLsIagKe4ljMzFq63XX+MlocJ0sbfulperqJRt7iWEvMWo7fnVUQpP7m+m9997Db7/9hsTERACAj48PAKCgoMBgvYKCAv1tPj4+V/z2s8LCQoN1rnQfl26D/sff31/2EFSFvcSxmRi19Ro46j4ozvdJ276dnbp6ycZeYthLjNqOXx2VUafl3GxLlixBTEwMYmJi0Ngo8XvLAUydOhXTpk0D0DT5j42NBQBkZmaioqIC/fr1AwBcvHgRJ06c0N9eX1+PxMREhIaGwsXFBQBw+PBheHl5oWvXrgCAM2fOoKamRv/B4aKiImRkZOj/laG2thaHDh1CWFiY/jy1pKQk+Pv7w8/PDwCQnp6OhoYG/elQ+fn5OHfuHKKiogAA1dXVSE5ORkREBOzt7QEABw8eRGBgoP4vMydPnoRGo0HPnj0BALm5ucjJyUFkZCRCQ0ORkJCAI0eOYMCAAbCxsQEAJCQkoEePHujcuTMAIC0tDba2tujevTsAICsrC4WFhQgPDwcAlJeXIzU1FdHR0bCyanpaxcfHIzg4GJ6engCAo0ePwtnZGUFBQQCA8+fPo6SkBP379wcAlJaW4tixY4iJiYGFhQUURcH+/fvRt29f/dWaUlJS4O7ujoCAACn7aeDAgUhNTW33/QQAlZWVqtxPoaGh+Pbbb83i9XQz9lNoaCji4+NV83rSWGpgaWkND4/Y5sd/Bo2NNXB2DmneL0WorMzQfxlQY2MtSksPwdU1TH8+c2lpEuzs/GFn59d8H+lQlAY4OTXtp5qafFRVnYO7e9N+amioRllZMlxdI+DuHoXy8lSUlByEg0MgbG2b9pNWexIWFho4OjbtJ50uFzpdDtzcIpvvoxJlZUfg5jYAlpZN+6mkJAGOjj1gY9O0nyoq0mBpaQtHx+7Nz48s1NYWwtU1vLlHOcrLU+HuHg0Li6b9VFwcDyenYNjYNO2n8vKjsLJyhoNDUPN9nEddXQlcXJr2U11dKSoqjsHDIwaABQAFxcX74ezcF9bWbs33kQJra3fY2zftp6qqTNTXV8DFpV9z44vQak/o94Gi1KOkJBEuLqGwsmraT2Vlh2Fj4wUvr1EoL09t9/2k0dg3N1bXfmrp1d77KTZWnce94cOH67erhv8+AXLne22l6RnajpYuXYq///3vGDp0KE6fPq3/eVBQEH7//XdERkYiOTlZ//Pt27fjwoULePjhh/HII4/gvffe08dtUVFRgWeffRZr167FunXr4OnpiXvvvVd/e0REBJKSkhAUFIRz5861OrakpCT9E9VcxMbGIj4+XvYwVIO9xLGZGLX1mr/+K/TtA1hlj5WyfQ+PWBQXq6eXbOwlRlav++9v903eFGo7fsnWVvPOa56Wo9Fo8OSTT8LX1/eGN7Zs2TKMHz8ed955p8HEHmj6m1NeXh7i4uL0P7O1tUVsbKz+HPvExEQ4OzsjOjpav050dDScnJwM1omNjYWtra1+nbi4OOTk5Fx1Ym+u0tPTZQ9BVdhLHJuJYS8xlZXsJYK9xLCXGB6/TMM1T8tpaGjA4sWL8d13393Qhj744ANMnDgRf/vb31BSUgJvb28AgFarRWVlJYCmyf+cOXNw6tQppKenY+7cudBqtdi4cSMA4NSpU/jhhx+wcuVK/ak0K1euxLZt2/RPqI0bN2LevHlYu3Yt3nzzTfTs2ROzZ8/G/Pnzb2j8HVVDQ4PsIagKe4ljMzE30svIrx25qeq7tP82L6UofH6JYC8x7CWGx3vTYNQHag8ePHjDX0zw9NNPw8XFBT///DPy8/P1y4svvqhfZ9GiRXj33XexYsUKJCcnw9fXFyNGjIBWq9WvM2HCBKSmpmLHjh3YsWMHUlNTMXHi/66xXF5ejri4OPj5+SE5ORkrVqzAkiVLeBnMVrScg0fGYS9xbCZGbb2sssdKOyUHgP58bzIOe4lhLzFqO351VEZ9oHb16tVYsmQJAgICcPjwYf077S1+++23a96HhYWFUQOaP3/+Vd9lLy0tNZjMX8nx48cxePBgo7ZHRERERNRRGDW5bzkt5krvfiuKov9UM6lPfn6+7CGoCnuJYzMxauvV4P4EAEBT8h8p26+pUVcv2dhLDHuJUdvxq6Myalbecskh6nj4IWMx7CWOzcSorZfiOLzp/0ia3FdVnZOyXbViLzHsJUZtx6+Oyqhz7v/444+rLqReLde5JeOwlzg2E8NeYlquqU7GYS8x7CWGxy/TYPQ31I4cORLbtm1DWloaunRpujzClClTcOedd7bZ4IiIiIiIyHhGTe4nTJiAL7/8EhkZGQgKCoK1tTWApmvgz5o1q00HSG2rurpa9hBUhb3EsZkY9hLT0MBeIthLDHuJ4fHLNBg1uZ81axamTp2KGTNmoL6+Xv/zgwcP6r/yl9Tp0m8DpmtjL3FsJkZ1vRRd0yJJWZnKeknGXmLYS4zqjl8dlFGT+x49eiAxMfGyn2u1Wri4uNz0QVH7iYiIkD0EVWEvcWwmRm29rHIeglXOQ9K27+qqrl6ysZcY9hKjtuNXR2XU5D43Nxc9e/a87Od33HEHzp49e9MHRe3H3t5e9hBUhb3EsZkY9hKj0bCXCPYSw15iePwyDUZN7letWoXly5dj4MCBAICuXbti0qRJWLRoET766KM2HSAREZmuBo/n0eDxvOxhEBFRM6Ouc7948WK4urpi165dsLOzw549e1BTU4N33nkHH374YVuPkdrQwYMHZQ9BVdhLHJuJUVsvxSGm6f8UL5Oy/ZISdfWSjb3EsJcYtR2/OiqjL4U5d+5cdOrUCbfffjuioqLQuXNn/POf/2zLsVE7CAwMlD0EVWEvcWwmhr3EODgEyh6CqrCXGPYSw+OXaTB6cg8AiqJAp9OhqqoKDQ0NbTUmakc+Pj6yh6Aq7CWOzcSwlxhbW/YSwV5i2EsMj1+mwajJvY2NDd59910UFxcjNTUVR48eRXFxMZYtWwZbW9u2HiMRERERERnBqHPuP/roI4wYMQKPPfaY/pKY0dHRePvtt+Hs7IwpU6a06SCp7Zw8eVL2EFSFvcSxmRjV9Wookbp5rVZlvSRjLzHsJUZ1x68OyqjJ/f/7f/8PDzzwAHbv3q3/WWZmJgoLC/H1119zcq9iGo1G9hBUhb3EsZkYtfWyypsqdfsWFurqJRt7iWEvMWo7fnVURp2WU1lZiZycnMt+npOTw68aVrkrfX8BtY69xLGZGPYS4+jIXiLYSwx7ieHxyzQYNbl///33MW/ePNjZ2el/Zmdnh9deew3vv/9+mw2OiIhMW0OnV9DQ6RXZwyAiomatnpazZcsWgz8PGTIEOTk5OHr0KACgb9++sLKygqOjY9uOkNpUbm6u7CGoCnuJYzMxauul2IVL3b5Op65esrGXGPYSo7bjV0fV6uT+4sWLBn/++uuvDf6cmZnZNiOidnWl062odewljs3EsJcYnY69RLCXGPYSw+OXaWh1cv/oo4+25zhIksjISMTHx8sehmqwlzg2E8NeYtzcIlFczF7GYi8x7CWGxy/TIPQlVkREREREZLqMuhSmm5sbXn/9dQwdOhReXl6wtDT8O4G3t3ebDI7aXmVlpewhqAp7iWMzMarrVZ8ndfMNDSrrJRl7iWEvMao7fnVQRk3u169fj5CQEKxbtw4FBQVQFKWtx0Xt5MiRI7KHoCrsJY7NxKitl1X+s1K3X1amrl6ysZcY9hKjtuNXR2XU5H7IkCEYPHgwfvvtt7YeD7WzAQMG4NChQ7KHoRrsJY7NxLCXGDe3ASgtZS9jsZcY9hLD45dpMOqc+7Nnz152Kg51DDY2NrKHoCrsJY7NxKitV0Pn+WjoPF/a9i0t1dVLNvYSw15i1Hb86qiMmrFPnz4db7/9Nvr168dJPhER6Sm2IVBsQ2QPg4iImhl1Ws6ZM2dgb2/f6rlUVlZG3Q2ZoISEBNlDUBX2EsdmYthLTEkJe4lgLzHsJYbHL9Ng1Kz8888/h6urK5577jl+oLaD6dGjB06dOiV7GKrBXuLYTAx7iXF07AGtlr2MxV5i2EsMj1+mwajJfUREBG6//XakpaW19XionXXu3JkvRAHsJY7NxLCXGBubzgDYy1jsJYa9xPD4ZRqMmtyfOHECLi4ubT0WIiJSm9rfZY+AiIguYdSnY+fOnYulS5di2LBh8PLygru7u8FirNjYWGzZsgXZ2dlQFAWTJ082uH3NmjVQFMVgSUxMNFjHxsYGy5cvR1FREbRaLbZs2QJ/f3+Ddbp27YqtW7dCq9WiqKgI7733HqytrY0epznhv8aIYS9xbCZGbb2sCmfBqnCWtO1XVKirl2zsJYa9xKjt+NVRGfXO/ffffw8A2Llzp8H59hYWFlAUxegP1Do5OeH48eNYv3491q9ff8V1du3ahYkTJ+r/XFtba3D7smXL8Ne//hXjx4/HxYsXsXTpUmzfvh3h4eFobGyEpaUlvvvuO1y8eBGxsbHw9PTEunXrYGFhgeeee86ocZoTW1tb2UNQFfYSx2Zi2EuMpSV7iWAvMewlhscv02DUrHzo0KE3ZWM//PADfvjhBwDA2rVrr7hOTU0NCgoKrnibi4sLpkyZgkceeQS7d+8GAEycOBHnz5/H8OHDsXPnTowYMQIhISEICAhAdnY2AGDWrFn4+OOP8eqrr6KiouKmPJaOonv37sjLk/v18WrCXuLYTIzaetV7LQIAae/eOzp2R02NenrJxl5i2EuM2o5fHZVRk/t9+/a19Tj0YmJiUFBQgNLSUvzyyy949dVXUVRUBAAIDw+HjY0Ndu7cqV8/OzsbJ0+exMCBA7Fz505ER0fj5MmT+ok9AOzYsQN2dnYIDw/H3r172+2xEBF1eDZ/kT0CIiK6hFGT+9tuu+2qt//22283ZTA//vgjvvnmG2RmZiIwMBBvvvkmfv75Z4SHh6O2thY+Pj6or6/HhQsXDH6voKAAPj4+AAAfH5/L3vm/cOEC6uvr9etcaurUqZg2bZr+d2NjYwEAmZmZqKioQL9+/QAAFy9exIkTJ/S319fXIzExEaGhofoPGx8+fBheXl7o2rUrgKbvB6ipqUFISNMXvBQVFSEjIwMDBw4E0HTK0aFDhxAWFgZHR0cAQFJSEvz9/eHn5wcASE9PR0NDA3r37g0AyM/Px7lz5xAVFQUAqK6uRnJyMiIiImBvbw8AOHjwIAIDA/WP9+TJk9BoNOjZsycAIDc3Fzk5OYiMjISvry/CwsJw5MgRDBgwQP/tcgkJCejRowc6d+4MoOk8OltbW3Tv3h0AkJWVhcLCQoSHhwMAysvLkZqaiujoaP1pWvHx8QgODoanpycA4OjRo3B2dkZQUBAA4Pz58ygpKUH//v0BAKWlpTh27BhiYmL0p3zt378fffv2hZubGwAgJSUF7u7uCAgIkLKffH19ERsb2+77CQAqKytVuZ98fX3RpUsXs3g93Yz95Ovri9DQ0OvaT1ZWJXBxadpPdXWlqKg4Bg+PGAAWABQUF++Hs3NfWFu7NT8fUmBt7Q57+6bXU1VVJurrK+Di0q+56UVotSfg4dG0nxSlHiUliXBxCYWVVdN+KrTQwNLCWr9OZeUZNDbWwNk5pPk+ilBZmQF396b91NhYi9LSQ3B1DYNG49j8nEqCnZ0/7Oz8mu8jHYrSACenpv1UU5OPqqpzcHdv2k8NDdUoK0uGq2sEbG194eERi5KSg3BwCIStbdN+0mpPwsJCA0fHpv2k0+VCp8uBm1tk831UoqzsCNzcBui/hbSkJAGOjj2ar5DSdL61paUtHB27Nz8/slBbWwhX1/Dm5205ystT4e4eDQuLpv1UXBwPJ6dg2Nh4Njc+CisrZzg4BDXfx3nU1bX/fiorOwwbGy99r/beTxqNfXNjde2nll7tvZ9iY9V53PP09NQf19Xw3ydA7nyvrTQ9Q6+hoaEBiqLAwsJC/7NLz72/ni+xqqiowDPPPIN169a1uo6vry/Onz+PcePGYfPmzRg/fjzWr19/2Ydjf/rpJ2RkZOCJJ57AypUr0b17dwwbNsxgnbq6OkycOBFffPFFq9tLSkrSP1HNhYODA6qqqmQPQzXYSxybibmRXlu33uTBGKG+y1cAAKvsse2/cQAajQMaGvj8MhZ7iZHV6/77232TNwWP92Laat5p1Ky85W9ILaytrXHbbbfh1VdfxSuvvHLTB9UiLy8P2dnZ6NGjB4Cmv21aWVmhU6dOBu/ee3t7Iz4+Xr/OoEGDDO6nU6dOsLKyQn5+fpuNVa3Cw8P17eja2Escm4lhLzGuruEoLmYvY7GXGFm9ZPxFvcWN/MWCxy/TYNTk/o8//rjsZ2fPnkVZWRnmzZuHH3/88aYPDAA8PT3h7++v/3DG4cOHUVtbi7i4OHz++ecAAH9/f/Tu3Vv/lceJiYmYO3cu/P39kZOTAwCIi4uDTqfD4cOH22ScRETmyqKGl74jIjIl4ufTXCIzM1N/npMxHB0d9edaWVpaolu3bggNDUVxcTGKi4vx+uuv4+uvv0ZeXh4CAwPx9ttvo7CwEJs3bwbQdD7WJ598gkWLFqGwsFB/KcyjR4/qr56zc+dOpKWlYf369Zg5cyY8PT2xePFirF69mlfKuYLy8nLZQ1AV9hLHZmLU1ktTNE/q9uvr1dVLNvYSw15i1Hb86qiMOuf+z19UZWFhAV9fX7z++uv4y1/+ov8wxLUMHjz4ilerWbt2LZ588kl8++23uO222+Dm5oa8vDzs2bMHr732msGVb2xsbPDOO+9gwoQJsLe3x08//YSnnnrKYJ2uXbviww8/xJ133onq6mp89tlneOmlly67Zv6fmeM590TUccj8p3wi6hjUer6/GrXVvFPoA7UGv2hhgaysLIwbNw6HDh266QOTwRwn99HR0Zd9CzC1jr3EsZmYG+kl5QO1Pu8DAKzyn23/jQNwd49GSQmfX8ZiLzHm2OtGJvc83ouR+oHaP3+JVWNjI4qKinDmzBk0NDTc9EFR+7meKx2ZM/YSx2ZiVNfLylfq5lsubUjGYS8x7CVGdcevDsrkvsSKiIiIiIiuz1VPy/nzufatKSkpuVnjkcocT8shoo7DHK9zT0Q3F8+5bz9tNe+0vNqNFy5cQFFR0VWXwsLCmz4oaj/BwcGyh6Aq7CWOzcSwlxgnJ/YSwV5i2EsMj1+m4aqn5fz5XPtLjRw5EtOnT0d9ff1NHxS1n5avdCbjsJc4NhOjtl4WOrnfH2Jjo65esrGXGPYSo7bjV0d11cn9lc6179+/PxYvXozY2FisXLkSb7zxRpsNjoiITJvmwtuyh0BERJe46mk5lwoMDMRnn32GX3/9FRcvXkRwcDCmT5+OCxcutOX4qI0dPXpU9hBUhb3EsZkY9hJTXs5eIthLDHuJ4fHLNFxzcu/h4YFly5bh1KlT8PHxwcCBA/H3v/8dv//+e3uMj9qYs7Oz7CGoCnuJYzMxautV77sa9b6rpW3fykpdvWRjLzHsJUZtx6+O6qqT+zlz5uDs2bMYPHgw/vrXv2LYsGFITk5ur7FROwgKCpI9BFVhL3FsJkZ1vTTuTYskDg4q6yUZe4lhLzGqO351UFc95/6NN95AdXU1srOz8dRTT+Gpp5664np//etf22RwRERqc6OXo/TwAIqLb85YiIjI/Fx1cr9+/XooSquXwacO4Pz587KHoCrsJY7NxFRXs5cI9hLDXmLYSwyP96bhqpP7Rx55pL3GQZJ0lC8gay/sJY7NxNTVsZcI9hLDXmLYSwyP96bB6KvlUMfUv39/2UNQFfYSx2ZiXFz6yx6CEIuq/bCo2i9t+2rrJRt7iWEvMTzem4arvnNPRER0NZriZbKHQEREl+A792autLRU9hBUhb3EsZmYurpS2UNQFfYSw15i2EsMj/emgZN7M3fs2DHZQ1AV9hLHZmIqKtTVq97/U9T7fypt+2rrJRt7iWEvMTzemwZO7s1cTEyM7CGoCnuJYzMxHh4q62Vh17RIorpekrGXGPYSw+O9aeDk3sxZWFjIHoKqsJc4NhPFXmLYSwx7iWEvETzemwZO7s0cv8dADHuJYzNR7CWGvcSwlxj2EsHjvWng5N7M7d8v7xJ2asRe4thMTHExe4lgLzHsJYa9xPB4bxo4uTdzffv2lT0EVWEvcWwmxtlZXb0sKnfDonK3tO2rrZds7CWGvcTweG8aeJ17M+fm5iZ7CKrCXuLYTIy1tZvsIQjRlPxH6vbV1ks29hLDXmJ4vDcNfOeeiIiIiKiD4OTezKWkpMgegqqwlzg2E1NeniJ7CELqu3yF+i5fSdu+2nrJxl5i2EsMj/emgZN7M+fu7i57CKrCXuLYTIy1NXuJYC8x7CWGvcTweG8aOLk3cwEBAbKHoCrsJY7NxNjbs5cI9hLDXmLYSwyP96aBk3siIiIiog6Ck3szl5mZKXsIqsJe4thMTFUVe4lgLzHsJYa9xPB4bxp4KUwzV1FRIXsIqsJe4thMTH29unpZVGyTun219ZKNvcSwlxge700D37k3c/369ZM9BFVhL3FsJsbFRV29NGXroClbJ237auslG3uJYS8xPN6bhnad3MfGxmLLli3Izs6GoiiYPHnyZevMmzcPOTk5qKqqwp49exAcHGxwu5ubG9avX4/S0lKUlpZi/fr1cHV1NVinT58+2Lt3L6qqqpCdnY3XXnutTR8XEZG5UizsoVjYyx4GERE1a9fJvZOTE44fP47p06ejqqrqsttnzZqFmTNn4tlnn0VkZCQKCwuxa9cuODk56dfZuHEjwsLCMHLkSIwcORJhYWHYsGGD/nZnZ2fs2rULBQUFiIyMxPTp0/HSSy9hxowZ7fIY1ebixYuyh6Aq7CWOzcTU1qqrV4P/BjT4b7j2im1Ebb1kYy8x7CWGx3vTYAFAkbHhiooKPPPMM1i37n//nJubm4sPPvgACxYsAADY2dmhsLAQL774IlatWoVevXrh5MmTGDRoEBISEgAAgwYNwv79+3HrrbciPT0dTzzxBP7973/D29sbOp0OAPDqq6/iySefRJcuXa46pqSkJERGRrbRIyYic7B1q+wRtK+WL7Cyyh4reSREdDPcf7/sEZiPtpp3msw590FBQfD19cXOnTv1P9PpdNi3bx8GDhwIAIiOjkZFRYV+Yg8ABw4cgFarNVgnPj5eP7EHgB07dsDf3x+BgYHt82BUJDY2VvYQVIW9xLGZGA8P9hLBXmLYSwx7ieHx3jSYzNVyfHx8AAAFBQUGPy8oKIC/v79+naKiost+t7CwUP/7Pj4+yM7Ovuw+Wm47d+6cwW1Tp07FtGnT9Le3PDEzMzNRUVGh/3DIxYsXceLECf3t9fX1SExMRGhoKFxcXAAAhw8fhpeXF7p27QoAOHPmDGpqahASEgIAKCoqQkZGhv4vIrW1tTh06BDCwsLg6OgIoOlvcf7+/vDz8wMApKeno6GhAb179wYA5Ofn49y5c4iKigIAVFdXIzk5GREREbC3bzrv9eDBgwgMDNQ3OXnyJDQaDXr27Amg6V9IcnJyEBkZidDQUFRWVuLIkSMYMGAAbGxsAAAJCQno0aMHOnfuDABIS0uDra0tunfvDgDIyspCYWEhwsPDAQDl5eVITU1FdHQ0rKyanlbx8fEIDg6Gp6cnAODo0aNwdnZGUFAQAOD8+fMoKSlB//79AQClpaU4duwYYmJiYGFhAUVRsH//fvTt2xdubm4Amr7a2t3dXf9FGe29n0JDQ6XsJwCq3U+hoaFm83oCAFfXSpSVHYGb2wBYWjbtp5KSBDg69oCNTdN+qqhIg6WlLRwduzdvNwu1tYVwdQ2Hi0so6uvLUV6eCnf3aFhYNO2n4uJ4ODkFw8bGs3lfHoWVlTMcHIKa7+M86upK4OLStJ/q6kpRUXEMHh4xaPlH2uLi/XB27gtra7fm+0iBtbW7/ot6qqoyUV9fof8QYW3tRWi1J/QTHEWpR0lJIlxcQmFl1bSfCi00sLSw1q9TWXkGjY01cHYOab6PIlRWZsDdvWk/NTbWorT0EFxdw6DRODY/p5JgZ+cPOzu/5vtIh6I0wMmpaT/V1OSjquoc3N2b9lNDQzXKypLh6hoBF5fQ5sYH4eAQCFvbpv2k1Z6EhYUGjo5N+0mny4VOlwM3t8jm+7ix/dT0vFXPfiorOwwbGy99r/beTxqNvUr3U6iU/WRv31XafoqNvf7/Pl36OUk1/PcJkDvfaysmc1pOdHQ0EhIS0K1bN2RlZenX++STT+Dv74+RI0filVdewWOPPYZbbrnF4L7Onj2L1atXY+HChdixYweys7MxZcoU/e1du3bFH3/8gejoaBw8eLDVMZnjaTnR0dFITEyUPQzVYC9x5tbsRk/LcXePRkmJenrJPi1Hbb1kYy8x5tjrRk7LMbfj/Y1qq3mnybxzn5+fDwDw9vY2mNx7e3vrb8vPz9f/DfBSXl5eBut4e3sb3N7y55Z16H/4IhTDXuLYTIzaJhKW5V9K3b7aesnGXmLYSwyP96bBZM65z8zMRF5eHuLi4vQ/s7W1RWxsrP4c+8TERDg7OyM6Olq/TnR0NJycnAzWiY2Nha2trX6duLg45OTkXHZKDkF/mgkZh73EsZmYltMA1MKy/EupE3y19ZKNvcSwlxge701Du07uHR0dERoaitDQUFhaWqJbt24IDQ3Vn7O0bNkyvPzyyxg9ejRCQkKwdu1aaLVabNy4EQBw6tQp/PDDD1i5ciWioqIQFRWFlStXYtu2bUhPTwfQdKnMqqoqrF27FiEhIRg9ejRmz56NpUuXtudDVY2W88fIOOwljs3EtJx7qxaKpQcUSw9p21dbL9nYSwx7ieHx3jS06+Q+IiICKSkpSElJgYODA/71r38hJSUF//rXvwAAixYtwrvvvosVK1YgOTkZvr6+GDFiBLRarf4+JkyYgNTUVOzYsQM7duxAamoqJk6cqL+9vLwccXFx8PPzQ3JyMlasWIElS5Zwck9E1AYa/FahwW+V7GEQEVEzaR+oNUXm+IFaBweHK36hGF0Ze4kzt2Y3+oFajcYBDQ3q6SX7A7Vq6yUbe4kxx1438oFaczve36gOf517ksPLy0v2EFSFvcSxmRgbG/YSwV5i2EsMe4nh8d40cHJv5lo+70DGYS9xbCam5frWZBz2EsNeYthLDI/3psFkLoVJRHQz3ejpMURERGrEyb2ZO3PmjOwhqAp7iWMzMZWV6uplWbpe6vbV1ks29hLDXmJ4vDcNnNybuZqaGtlDUBX2EsdmYhob1dXLUiv3n0jU1ks29hLDXmJ4vDcNPOfezIWEhMgegqqwlzg2E+PsrK5eipUfFCs/adtXWy/Z2EsMe4nh8d408J17IiK6bg0+ywHIuxQmEREZ4uTezBUVFckegqqwlzg2E1Nby14i2EsMe4kxx143cjECJ6ciXPK9o0Ju5Pr6ZIin5Zi5jIwM2UNQFfYSx2ZiKivZSwR7iWEvMewlhr1MAyf3Zm7gwIGyh6Aq7CWOzcS4u7OXCPYSw15i2EsMe5kGTu6JiIiIiDoInnNv5mpra2UPQVXYSxybiWlsVFcvy5KVUrevtl6ysZcY9hLDXqbBAoAiexCmIikpCZGRkbKHQUQ3Ab+hlohIPczxA7VtNe/kaTlmLiwsTPYQVIW9xLGZGFdXdfVSrG+BYn2LtO2rrZds7CWGvcSwl2ngaTlmztHRUfYQVIW9xLGZGI1GXb0avP8NQN517tXWSzb2EsNeYtjLNPCdeyIiIiKiDoKTezOXlJQkewiqwl7i2ExMaSl7iWAvMewlhr3EsJdp4OTezPn7+8segqqwlzg2E2Nnx14i2EsMe4lhLzHsZRo4uTdzfn5+soegKuwljs3E2Nmxlwj2EsNeYthLDHuZBn6gloiIrptl8Xuyh0BERJfg5N7Mpaenyx6CqrCXODYTU1mprl6WVfFSt6+2XrKxlxj2EsNepoGn5Zi5hoYG2UNQFfYSx2ZiFEVdvRTbECi2IfK2r7JesrGXGPYSw16mgZN7M9e7d2/ZQ1AV9hLHZmKcnNTVq6HzfDR0ni9t+2rrJRt7iWEvMexlGji5JyIiIiLqIDi5N3P5+fmyh6Aq7CWOzcTU1LCXCPYSw15i2EsMe5kGTu7N3Llz52QPQVXYSxybiamqOid7CKrCXmLYSwx7iWEv08DJvZmLioqSPQRVYS9xbCbG3Z29RLCXGPYSw15i2Ms08FKYRER03SwvLJQ9BCIiugQn92auurpa9hBUhb3EsZmYhgZ19bLUJUvdvtp6ycZeYthLDHuZBp6WY+aSk+X+h1lt2Escm4kpK1NXr0a7CDTaRUjbvtp6ycZeYthLDHuZBpOa3M+bNw+KohgseXl5l62Tk5ODqqoq7NmzB8HBwQa3u7m5Yf369SgtLUVpaSnWr18PV1fX9nwYqhIRIe8/ymrEXuLYTIyrq7p6NXaajcZOs6VtX229ZGMvMewlhr1Mg0lN7gHg1KlT8PHx0S99+/bV3zZr1izMnDkTzz77LCIjI1FYWIhdu3bByclJv87GjRsRFhaGkSNHYuTIkQgLC8OGDRtkPBRVsLe3lz0EVWEvcWwmRqNhLxHsJYa9xLCXGPYyDSZ3zn19fT0KCgqueNvzzz+PhQsX4ptvvgEATJ48GYWFhZgwYQJWrVqFXr16YdSoURg0aBAOHjwIAHj88cexf/9+9OzZE+np6e32OIiIiIiI2pvJvXP/l7/8BTk5Ofj999/x+eefIygoCAAQFBQEX19f7Ny5U7+uTqfDvn37MHDgQABAdHQ0KioqkJCQoF/nwIED0Gq1+nXIUMtfgsg47CWOzcSUlLCXCPYSw15i2EsMe5kGk3rn/tChQ3j44Ydx6tQpeHl5Ye7cuUhISEBISAh8fHwA4LJ39QsKCuDv7w8A8PHxQVFR0WX3W1hYqP/9P5s6dSqmTZum//3Y2FgAQGZmJioqKtCvXz8AwMWLF3HixAn97fX19UhMTERoaChcXFwAAIcPH4aXlxe6du0KADhz5gxqamoQEhICACgqKkJGRob+Lxq1tbU4dOgQwsLC4OjoCABISkqCv78//Pz8AADp6eloaGhA7969ATR92+e5c+f01w6vrq5GcnIyIiIi9Kc/HDx4EIGBgfrHfPLkSWg0GvTs2RMAkJubi5ycHERGRqJLly44ffo0jhw5ggEDBsDGxgYAkJCQgB49eqBz584AgLS0NNja2qJ79+4AgKysLBQWFiI8PBwAUF5ejtTUVERHR8PKqulpFR8fj+DgYHh6egIAjh49CmdnZ/1f2M6fP4+SkhL0798fAFBaWopjx44hJiYGFhYWUBQF+/fvR9++feHm5gYASElJgbu7OwICAqTsp1tvvRXZ2dntvp8AoLKyUpX7qUuXLoiPj2/315OHB1BamgQ7O3/Y2fk1N0yHojTAyalpP9XU5KOq6pz+2swNDdUoK0uGq2uE/p+XS0oOwsEhELa2TftJqz0JCwsNHB2b9pNOlwudLgdubpHN91GJsrIjcHMbAEtLm+b7SICjYw/Y2DTtp4qKNFha2sLRsXvz8yMLtbWFcHUNh719F1RUnEB5eSrc3aNhYdG0n4qL4+HkFAwbG8/mfXkUVlbOcHAIar6P86irK4GLS9N+qqsrRUXFMXh4xACwAKCguHg/nJ37wtrarfk+UmBt7Q57+6bXU1VVJurrK+Di0q+56UVotSfg4dG0nxSlHiUliXBxCYWVVdN+KrTQwNLCWr9OZeUZNDbWwNk5pPk+ilBZmQF396b91NhYi9LSQ3B1DYNG49j8nLr+/eTk1APV1dntvp+anrflqtlPZWWHYWPjBQ+PaFRXZ7f7fpL1errR/eTtfS+qq7PbfT/Z23dtbqyu/eTj81fU1BRc134KDTWPecSl87220vQMNVGOjo74/fffsXDhQhw8eBAJCQno1q0bsrKy9Ot88skn8Pf3x8iRI/HKK6/gsccewy233GJwP2fPnsXq1auxcOHVr8eclJSkn0iZi9jYWMTHx8sehmqwlzhZzbZubfdN3hQeHrEoLlbPc0yxbfqPmUVNmpTtq62XbOwlhr3E3Eiv+++/yYNRgbaad5rUO/d/VllZibS0NPTo0QPffvstAMDb29tgcu/t7Y38/HwATe+WtryDeSkvLy/9OkTUvt5+Gygulj0KaiuyJvVERHRlJnfO/aVsbW3Rq1cv5OXlITMzE3l5eYiLizO4PTY2Vn+OfWJiIpydnREdHa1fJzo6Gk5OTgbn4dP/nDx5UvYQVIW9xGm1bCZCbb0aHWLR6BArbftq6yUbe4lhLzHsZRpM6p37xYsXY9u2bfjjjz/g5eWF1157DY6Ojli3bh0AYNmyZZgzZw5OnTqF9PR0zJ07F1qtFhs3bgTQdBnNH374AStXrtSfR79y5Ups27aNV8pphUajkT0EVWEvcRYWbCZCbb0aPaYDACyr5Jy6oLZesrGXGPYSw16mwaTeue/SpQs+//xznD59Gt988w1qamoQFRWFP/74AwCwaNEivPvuu1ixYgWSk5Ph6+uLESNGQKvV6u9jwoQJSE1NxY4dO7Bjxw6kpqZi4sSJsh6SyWv58CYZh73EtXwAi4zDXmLYSwx7iWEvMexlGkzqnfvx48dfc5358+dj/vz5rd5eWlrKyTwRERERmSWTeuee2l9ubq7sIagKe4nT6dhMBHuJYS8x7CWGvcSwl2ng5N7M5eTkyB6CqrCXOJ2OzUSwlxj2EsNeYthLDHuZBk7uzZy5Xdf/RrGXuJYvOSHjqK2XpuBlaApelrZ9tfWSjb3EsJcY9jINJnXOPRERqYtF3VnZQyCiDkDmFw92tC/Q4jv3Zq6yslL2EFSFvcQ1NLCZCLX1anSMQ6Nj3LVXbCNq6yUbe4lhLzHsZRr4zr2ZO3LkiOwhqAp7iSsrYzMRauvV6P44AMCycpeU7autl2zsJYa9xLCXaeA792ZuwIABsoegKuwlzs2NzUSwlxj2EsNeYthLDHuZBk7uzZyNjY3sIagKe4mztGQzEewlhr3EsJcY9hLDXqaBk3siIiIiog6Ck3szl5CQIHsIqsJe4kpK2EwEe4lhLzHsJYa9xLCXaeDk3sz16NFD9hBUhb3EOTqymQi19dLkPwdN/nPStq+2XrKxlxj2EsNepoGTezPXuXNn2UNQFfYSZ2PDZiLU1suiPhcW9fK+cl5tvWRjLzHsJYa9TAMn90REdN0ane5Ho1MH+wYYIiIV43XuzVxaWprsIagKe4mrqGAzEWrr1eg2CQBgqZXz9ZJq6yUbe4lhLzHsZRr4zr2Zs7W1lT0EVWEvcZaWbCaCvcSwlxj2EsNeYtjLNHByb+a6d+8uewiqwl7iHB3ZTAR7iWEvMewlhr3EsJdp4Gk5RGZgq5wzJoiIiKid8Z17M5eVlSV7CKrCXuKqq9lMBHuJYS8x7CWGvcSwl2ngO/dmrrCwUPYQVIW9xNXWspkItfXS5E6Tun219ZKNvcSwlxj2Mg18597MhYeHyx6CqrCXOFdXNhOhtl4WjcWwaCyWtn219ZKNvcSwlxj2Mg2c3BMR0XVrdHkQjS4Pyh4GERE142k5Zq68vFz2EFSFvcTV17OZCLX1apnYW5Z/KWX7auslG3uJYS8xau3V0S46wXfuzVxqaqrsIagKe4krL2czEewlhr3EsJcY9hLDXqaB79ybuejoaCQmJsoehmrcaK+O9u6AMdzdo1FSwueYsdhLDHuJYS8x7CWGvUwD37k3c1ZW/PudCPYSZ2HBZiLYSwx7iWEvMewlhr1MAyf3REREREQdhAUARfYgTEVSUhIiIyNlD4M6MHM8LYc6NsXCHgBgoVRLHgkRkbr4+rbNvJPv3Ju54OBg2UNQFfYS5+TEZiLU1stCqZY6sVdbL9nYSwx7iWEv08DJvZnz9PSUPQRVYS9xNjZsJkJtvRpcJ6PBdbK07autl2zsJYa9xLCXaejQk/snn3wSv//+O6qrq5GcnIyYmBjZQyIi6lAU5/ugON8nexhERNSsw36s+cEHH8R7772Hp556Cvv378dTTz2FH374AcHBwcjKypI9PJNx9OhR2UNodzdy3ruV1VG89NLNG4s5KC83v+fYjWAvMewlhr3EsJcY9jINHfad+xkzZmDt2rX4+OOPcerUKTz33HPIy8vDk08+KXtoJsXZ2Vn2EFTFyoq9RLGZGPYSw15i2EsMe4lhL9PQId+5t7a2Rnh4ON555x2Dn+/cuRMDBw6UNCrTFBQUhOzsbCnbVuOVYxwcgqDTyemlVmwmhr3EsJcY9hLDXmLYyzR0yMl9p06dYGVlhYKCAoOfFxQUYPjw4ZJGZZrefhsoLpY9CiIiIiK6GTrk5F7E1KlTMW3aNADArbfeiqSkJMkjal+dOnWCre0F2cNQDfYSx2ZiVNerofl/feUcO1XXSzL2EsNeYthLzK233tpm9610tMXa2lqpq6tTxo4da/DzDz74QNm7d6/08ZnSkpSUJH0MalrYi83Yy7QW9mIv9jKdhb1Mo1eH/EBtXV0dDh8+jLi4OIOfx8XFISEhQdKoiIiIiIjaVoc9LWfp0qXYsGEDfv31Vxw4cABPPPEE/Pz88J///Ef20IiIiIiI2kSHndx/+eWX8PT0xNy5c+Hr64vjx4/j7rvvxh9//CF7aCZl1apVsoegKuwljs3EsJcY9hLDXmLYSwx7iWmrXhZoOj+HiIiIiIhUrkOec09EREREZI44uSciIiIi6iA4ue8gYmNjsWXLFmRnZ0NRFEyePNng9jVr1kBRFIMlMTHxmvd7xx13IDk5GdXV1Th79iwef/zxtnoI7epavf7cqmX54IMPWr3PgICAK/7OXXfd1dYPp83Nnj0bv/76K8rKylBYWIitW7ciJCTksvXmzZuHnJwcVFVVYc+ePQgODr7mfT/wwANIS0uDTqdDWloa/va3v7XBI2hf1+plZWWFhQsXIjU1FVqtFrm5ufjss8/QtWvXq97v4MGDr/gca8trJbcHY55fPIb9jzG9eAz7n6eeegqpqakoKytDWVkZEhIScPfddxusw2PX/1ytF49dV3at55iM45f063xyufFl1KhRyltvvaWMGTNGqaysVCZPnmxw+5o1a5SdO3cq3t7e+sXd3f2q9xkYGKhotVpl+fLlSq9evZTHHntMqa2tVR544AHpj7ete13aydvbW7nnnnsURVGUO+64o9X7DAgIUBRFUUaMGGHwu9bW1tIf740uP/74o/Lwww8rISEhSp8+fZRvvvlGycvLM3gOzZo1SykvL1ceeOABJSQkRNm0aZOSk5OjODk5tXq/UVFRSl1dnTJnzhylV69eypw5c5S6ujrl9ttvl/6Y27KXi4uLsnPnTuXBBx9UevbsqURGRir79u1T0tLSFI1G0+r9Dh48WFEURendu7fBc8zS0lL6Y27r5xePYWK9eAz733L//fcrI0eOVG655RalR48eyptvvqnU1tYqffv2VQAeu0R68dh1fc8xCccv+VG43NyloqLiipP7bdu2Cd3PwoULlfT0dIOfrV69WklISJD+GNu615+XVatWKadOnbrqOi3/YQwPD5f+mNp6cXR0VOrr65V7771X/7Pc3Fxlzpw5+j/b2dkp5eXlyrRp01q9ny+++ELZuXOnwc927dqlbNy4UfpjbOtef1569+6tKIqi9OnTp9V1Wv4D6enpKf0xtXcvHsNu7PnFY5jhcvHiRf2xiccusV5/Xnjsunaz9j5+8bQcMxITE4OCggKcPn0aq1atQufOna+6fnR0NHbu3Gnwsx07diAiIgJWVh32KqqXcXR0xN///nesXr3aqPW/+eYbFBQUYP/+/RgzZkwbj04OZ2dnaDQalJSUAACCgoLg6+tr8HzR6XTYt28fBg4c2Or9tPYcu9rvqNGfe12Ji4sLAFx1nRbJycnIzc3F7t27MWTIkJs1TJPRWi8ew67sWs8vHsP+x9LSEuPGjYOTkxMSEhJ47LqGP/e6Eh67DLXWrL2PX9L/dsPl5i5Xeid63Lhxyn333af06dNHuffee5WUlBTl2LFjio2NTav3c/r0aeW1114z+FlsbKyiKIri4+Mj/XG2Za9Ll6lTpyo6nU7p1KnTVe/H09NTmTFjhjJgwAAlPDxcmT9/vlJfX6/84x//kP4Yb/ayadMm5ciRI/p/Uo2OjlYURVG6du1qsN4nn3yi/Pjjj63eT01NjTJx4kSDn02cOFHR6XTSH2Nb9vrzYm1trezfv1/ZsmXLVe+nZ8+eyuOPP66EhYUpUVFRyooVK5SGhgYlJiZG+mNs6148hl3/84vHMCh9+vRRKioqlLq6OqWkpES5++67FYDHLtFef1547DKumYTjl/wgXG7uYsxpJr6+vkptba0yevToVtcxl/8wXqvXr7/+qmzatOm67nvFihVKamqq9Md4M5clS5YoOTk5SlBQkP5n/A+kWK9LF41Go2zatEk5fvy44uHhIXz/33333TX/w6qm5Vq9WhYew4zvxWNY0yT0lltuUcLCwpQFCxYoRUVFSkhICI9dgr0uXYfHLvFmLUtbH794Wo6ZysvLQ3Z2Nnr06NHqOvn5+fD29jb4mbe3N+rq6nDhwoW2HqJJCA0NRWRkpNH/nP1nhw4dumpjtVm6dCnGjx+PO++8E5mZmfqf5+fnA8AVny8tt11Ja8+xq/2OmrTWq4VGo8Hnn3+Ofv36YdiwYSguLhbeRkd6jl2r16V4DDOuF49hTerq6nD27FkcOXIEc+bMQUpKCl544QUeu1rRWq8WPHZd7lrNLtXWxy9O7s2Up6cn/P39kZeX1+o6iYmJiIuLM/hZXFwckpOTUV9f39ZDNAnTpk3D77//jt27d1/X7/fv3/+qjdVk2bJl+onE6dOnDW7LzMxEXl6ewfPF1tYWsbGxrZ6nCbT+HLva76jF1XoBTZeU27RpE/r164ehQ4eioKDgurbTUZ5j1+r1Z+Z+DDO2F49hV2ZpaQlbW1seu4zU0gvgsctYlzb7s/Y4fkn/pwwuN744OjoqoaGhSmhoqFJZWam89tprSmhoqNK1a1fF0dFRWbx4sRIVFaUEBAQogwcPVhISEpSsrCyDS32tW7dOWbdunf7PLZdhevfdd5VevXopU6ZMUWpqalR/Gblr9WpZx97eXiktLTW4isKly4IFC5Tdu3fr/zxp0iRl/PjxSq9evZSePXsqM2fOVGpqapTnn39e+uO90eWDDz5QysrKlKFDhxpcysvR0VG/zqxZs5TS0lJl9OjRSkhIiPL5559fdjm53bt3KwsWLND/OTo6Wqmrq1Nefvll5dZbb1Vmz56t1NbWqv5yctfqpdFolM2bNyvZ2dnKbbfdZrCOnZ2d/n7+/JqcPn268te//lXp3r27EhwcrCxYsEBRFOWq/7SrhuVavXgME+vVsvAY1rS8/fbbSkxMjBIQEKD06dNHWbBggdLQ0KCMHDlSAXjsEunFY5d4M0nHL/lRuNz40nKZqT9bs2aNYmdnp/z4449KQUGBUlNTo5w7d05Zs2aN0qVLF4P72LNnj7Jnzx6Dn91xxx3K4cOHFZ1Op/z+++/K448/Lv2xtnWvlnUefvhhpa6uTvH19b3ifaxZs0bJzMzU/3nSpElKWlqaotVqlbKyMiUpKalDfBANwBVbKYqizJs3z2C9efPmKbm5uUp1dbWyd+/ey843zMzMNGgMQBkzZoxy8uRJpaamRjlx4kSHONhfq1fLJQev5NLPf/z5NfnSSy8p6enpSlVVlXLx4kVl3759yqhRo6Q/3rbuxWOYWK+Whcew/z3Oc+fOKTqdTikoKFB27dqljBgxwmAdHruM68Vjl3gzGccvi+b/Q0REREREKsdz7omIiIiIOghO7omIiIiIOghO7omIiIiIOghO7omIiIiIOghO7omIiIiIOghO7omIiIiIOghO7omIqEMJCAiAoigIDw9vk/vftm0b1qxZ0yb3TUR0ozi5JyIScNttt6G+vh779++XPZQ2s2fPHrz//vvXXG/NmjVQFAUff/zxZbctXLgQiqJg27ZtbTHEq8rKyoKPjw9SUlIAAIMHD4aiKPD09Gz3sRARtTdO7omIBDz22GP48MMP0adPH/Tq1Uv2cPSsra2lbPePP/7Agw8+CAcHB/3PNBoNJk2ahPPnz0sZU2NjIwoKCtDQ0CBl+0REMnFyT0RkJDs7O0yYMAGrVq3CV199hSlTphjc3nI6yPjx4xEfH4/q6mqcPHkScXFx+nVa3kW+55578Ntvv6G6uhrJyckICwvTr+Ph4YGNGzciKysLVVVVOH78OB5++GGDbe3ZswcffvghFi9ejMLCQhw4cAAA0Lt3b2zfvh3l5eUoKCjAxo0b4e3trf+9NWvWYNu2bXjuueeQnZ2N4uJi/Pe//4W9vb3+9iFDhuCZZ56BoihQFAUBAQGtNjl69CgyMjLw4IMP6n92zz33QKfTYe/evQbrRkREYMeOHSgqKkJZWRni4+MRFRVlsE6PHj2wd+9eVFdX49SpUxg1ahQqKiowefJkg8YPPPAAdu7cicrKSqSlpWH48OGX7Yfw8HAEBATox3HhwgUoiqI/peZK/0LR0qeFvb091qxZg4qKCuTn5+OVV165rIG1tTUWLlyIrKwsVFZW4tdff8WIESNabUZE1JY4uSciMtLYsWNx/vx5HD9+HBs2bMCkSZNgZWV12XqLFi3C8uXL0b9/f+zatQtbtmyBn5+fwTrvvPMOXn75ZUREROD333/H9u3b9RNsOzs7HDlyBPfeey9CQkLw3nvvYeXKlbjzzjsN7uOhhx6ChYUFYmNjMWnSJPj4+GDfvn04fvw4br/9dgwfPhxOTk7YsmULLCws9L8XGxuLPn36YPjw4Rg3bhxGjx6N6dOnAwCmT5+OhIQE/Pe//4WPjw98fHyQlZV11S6ffPIJHn30Uf2fH330Uf0pO5dydnbGhg0bEBsbi9tvvx0pKSn4/vvv4eHhAQCwsLDA5s2bUV9fj6ioKDz88MOYN28ebG1tL9vmW2+9heXLlyM0NBRJSUn44osv4OjoeNl6WVlZeOCBBwAAwcHB8PHx0T9WY7zzzjuIi4vDmDFjMGzYMNx222244447DNZZs2YNBg8ejAkTJqBPnz5Yt24dtm3bhn79+hm9HSKim0nhwoULFy7XXvbs2aPMnDlT/+fMzExlzJgx+j8HBAQoiqIoc+bM0f/MwsJCOX36tPLGG28oAJTBgwcriqIoEyZM0K/j6OiolJSUKFOmTGl1259//rmyevVqg7GkpqYarDN//nxl9+7dBj9zc3NTFEVRIiMjFQDKmjVrlD/++EOxtLTUr7Nq1Spl165dBvf9/vvvX7PHmjVrlG3btilubm5KVVWV0r17d8Xb21vR6XRK165d9bdf7T5yc3OVf/zjHwoAZcSIEUpdXZ3i5+envz06OlpRFEWZPHmyQeNp06bp1/Hz81MURVEGDRpksE54eLhBc09Pz8v2558f56VjdnR0VHQ63RX31Zo1axQAyl/+8heloaFB6dq1q8H9bN68WVmxYoX05ywXLlzMb7n8LSciIrrMLbfcgpiYGEyYMEH/s88++wxTpkzB119/bbBuYmKi/v8rioJDhw4hODi41XUqKytx7Ngx/TqWlpaYPXs2xo0bB39/f9ja2sLGxuay01wOHz5s8Ofw8HDccccdqKiouOL4k5KSAAAnTpxAY2Oj/rbc3FwMGDDAmAxXVFpais2bN+PRRx9FaWkp9u7de8V3+zt37ow33ngDQ4cOhbe3NzQaDezt7dGtWzcAQK9evZCbm4vc3Fz97yQlJV3x3PmjR48ajB8AvLy8rvsxXMktt9wCW1vbK+6rFmFhYbC0tMSJEycMftfW1hY///zzTR0PEZExOLknIjLCY489BisrK/zxxx/6n7Wc6tKlSxdkZ2fftG29+OKLmDlzJqZPn45jx45Bq9ViwYIFl01eKysrDf5saWmJ7777Di+++OJl91lQUKD//3V1dQa3KYoCS8sbO0vzv//9L9atWwetVot//vOfV1xn3bp18Pb2xgsvvIBz586hpqYGP/30E2xsbIS39+fHAED4MTQ2NhqcrgSIfzDZ0tISjY2NiIyMvGxM1dXVQvdFRHQz8Jx7IqJr0Gg0mDx5MmbPno3+/fvrl9DQUBw9ehSPPPKIwfp//pDo7bffjpMnT7a6joODA/r06aNfJyYmBtu2bcOnn36K1NRUnD17Fj179rzmOI8cOYKQkBCcP38eZ8+eNVi0Wq3Rj7e2thYajcbo9QHgp59+Qm1tLTp16oRvv/32iuvExMTg/fffx/fff48TJ06goqICvr6++ttPnToFPz8/g59FREQIj+XPamtrAeCy+ykqKjLYFgCEhobq///Zs2dRW1t7xX3V4rfffoOlpSV8fHwua37pv0AQEbUXTu6JiK7hnnvuQadOnbB69WqkpaUZLF988cVlk/snn3wSY8aMQc+ePbFs2TIEBATgo48+Mlhn7ty5GD58OIKDg/Hf//4XtbW12LhxIwAgPT0dw4YNw6BBg3Drrbfigw8+QFBQ0DXHuWLFCri6umLTpk24/fbbERQUhGHDhmHlypVwcnIy+vGeO3cOt99+OwICAuDp6XnZu9ut6devH4KCgvST6T9LT0/HQw89hN69eyMiIgJffPGFwbq7du3C6dOnsW7dOvTr1w8DBgzA0qVLUVdXd9mHc0WcP38ejY2N+v3Y8sHbn3/+GaNGjcJ9992Hnj17YsmSJejatav+9yorK/HJJ5/g3//+t8G+uvQvCRkZGfj000+xdu1ajBkzBkFBQQgPD8fMmTMxevTo6x4zEdH14uSeiOgapkyZgj179qC4uPiy2/7v//4PQUFBBpe7nD17NmbMmIHU1FSMHDkSo0ePRk5OjsHvzZ49G0uWLMGRI0fQo0cP3HvvvaiqqgIAvPnmm/j111/xww8/YN++faisrMRnn312zXHm5eVh0KBBaGxsxI8//oi0tDSsWLECNTU1qKmpMfrxvvPOO6itrcWJEydw4cIF/Tnx16LVaq94vn+LRx99FE5OTjh8+DC++OIL/Pe//8W5c+f0tyuKgtGjR8PW1ha//vor1q1bh7feeguKokCn0xk9/j/Lzc3FvHnz8NZbb6GgoAAffPABgKZTiVqWAwcOoKKiAps3bzb43RdffBF79uzB5s2bsWfPHhw/fhz79u0zWOeRRx7BmjVrsGjRIpw6dQrbt2/HHXfcIe06/0RE0j/Vy4ULFy4dYfnzVVqutLR25RYuV1769eunKIqihIWFSR8LFy5cuKhh4QdqiYjIZPztb39DZWUlMjIyEBgYiKVLlyIlJQVHjhyRPTQiIlXg5J6IiEyGs7Mz/v3vf6Nr164oKSnB3r178cILL8geFhGRalig6S18IiIiIiJSOX6gloiIiIiog+DknoiIiIiog+DknoiIiIiog+DknoiIiIiog+DknoiIiIiog+DknoiIiIiog/j/cIZVo+CsTB4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now let's plot the distribution of the app. mag.\n",
    "# Let's set a dark background\n",
    "plt.style.use('dark_background')\n",
    "\n",
    "# Set a default font size for better readability\n",
    "plt.rcParams.update({'font.size': 14})\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "counts, bins, _ = plt.hist(neo_df[\"app_mag\"],\n",
    "                           bins=np.arange(13.0, 36.0, 1.0),\n",
    "                           color=\"yellow\",\n",
    "                           alpha=0.7)\n",
    "\n",
    "plt.grid(linestyle=\"dashed\", alpha=0.3)\n",
    "plt.xlim(13, 35)\n",
    "plt.xlabel(\"Apparent Magnitude\")\n",
    "plt.ylabel(\"Number of NEOs\")\n",
    "plt.ylim(0, 4000)\n",
    "\n",
    "# Add the limiting magnitude of Pan-STARSS as a vertical line\n",
    "# https://panstarrs.ifa.hawaii.edu/pswww/?page_id=34\n",
    "plt.vlines(24, 0, 4000, linestyles=\"dashed\", color=\"lightblue\")\n",
    "plt.text(24.5, 3800, \"Pan-STARRS Lim. Mag.\", color=\"lightblue\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "57e6084a-5200-4bfd-88b3-0fe7cb52faa4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of NEOs Pan-STARRS could observe today: 5048\n",
      "Number of NEOs Pan-STARRS could NOT observe today: 23906\n"
     ]
    }
   ],
   "source": [
    "print(f\"Number of NEOs Pan-STARRS could observe today: {int(sum(counts[:11]))}\")\n",
    "print(f\"Number of NEOs Pan-STARRS could NOT observe today: {int(sum(counts[11:]))}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88ddd96b-f496-4c8e-830e-010120e9f206",
   "metadata": {},
   "source": [
    "# What's next?\n",
    "\n",
    "Our next session will continue with cell #18 (a new notebook will be created though)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}