12_NEO_Granvik_et_al_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c48ddd3c-8c20-4951-9bcf-d90e8c35fa40",
   "metadata": {},
   "source": [
    "# Granvik et al. (2018) NEO Model\n",
    "\n",
    "In this second-last NEO tutorial we will take a look at a so-called NEO population model. This model, by Granvik et al (2018), is a numerically generated model with the aim to create an un-biased population of Near-Earth Objects (NEOs). We will download and parse the data automatically and later we'll take a look at how one can create kernel density estimator (KDE) based sky maps!\n",
    "\n",
    "The details of the paper and model generation are beyond our scope. Take a direct look at their paper [here](https://www.sciencedirect.com/science/article/pii/S0019103517307017).\n",
    "\n",
    "Before we dive in, let's download the data first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a5abbe13-05bb-481b-b81c-69f6799326ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard libraries\n",
    "import datetime\n",
    "import math\n",
    "import pathlib\n",
    "import re\n",
    "import requests\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",
    "import tqdm\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": "56076f2e-4477-46d1-9a7a-22a5b68fe8fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "02a7c3ac-5c7a-44c5-817e-2dd62349cffc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Our Granvik et al. (2018) file path: raw_data/Granvik+_2018_Icarus.dat\n"
     ]
    }
   ],
   "source": [
    "# Define the directory and name of the Granvik data. The data are originally stored as a gzip\n",
    "# file and need to be unzipped\n",
    "data_dir = pathlib.Path(\"raw_data/\")\n",
    "data_dir.mkdir(parents=True, exist_ok=True)\n",
    "granvik_neo_model_file_zip = pathlib.Path(\"Granvik+_2018_Icarus.dat.gz\")\n",
    "granvik_neo_model_file = pathlib.Path(\"Granvik+_2018_Icarus.dat\")\n",
    "\n",
    "data_filepath = data_dir / granvik_neo_model_file\n",
    "download_filepath = data_dir / granvik_neo_model_file_zip\n",
    "print(f\"Our Granvik et al. (2018) file path: {data_filepath}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "04060dee-e1bc-4863-9749-f7b0c48f4e62",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download the Granvik Model\n",
    "dl_link = \"https://www.mv.helsinki.fi\" \\\n",
    "          \"/home/mgranvik/data/Granvik+_2018_Icarus/Granvik+_2018_Icarus.dat.gz\"\n",
    "response = requests.get(dl_link)\n",
    "download_file_path = pathlib.Path(download_filepath)\n",
    "with download_file_path.open(mode=\"wb+\") as file_obj:\n",
    "    file_obj.write(response.content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2fe66f22-f67a-4e65-b3a8-b46371649b38",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We use the standard libraries gzip and shutil to unzip the file\n",
    "import gzip\n",
    "import shutil\n",
    "\n",
    "with gzip.open(download_file_path, 'rb') as f_in:\n",
    "    with open(data_filepath, 'wb') as f_out:\n",
    "        shutil.copyfileobj(f_in, f_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ad69afd-ce73-4670-9cfe-935c1cdee8cd",
   "metadata": {},
   "source": [
    "## Parsing the Model\n",
    "\n",
    "Now we can open the model and load it into a Pandas dataframe. Additionally, we need to parse the orbital elements into a format that can later be used by SPICE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ccdafee0-1726-4115-8a7c-20e864165264",
   "metadata": {},
   "outputs": [],
   "source": [
    "granvik_model_df = pd.read_csv(data_filepath, names=[\"SemiMajorAxis_AU\",\n",
    "                                                     \"Ecc_\",\n",
    "                                                     \"Incl_deg\",\n",
    "                                                     \"LongAscNode_deg\",\n",
    "                                                     \"ArgP_deg\",\n",
    "                                                     \"MeanAnom_deg\",\n",
    "                                                     \"H\"],\n",
    "                               delim_whitespace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bd6bbbd-a9ef-48f7-a95e-6a2193a9f9e3",
   "metadata": {},
   "source": [
    "*The dataframe structure has been derived from the corresponding README file of the [repository](https://www.mv.helsinki.fi/home/mgranvik/data/Granvik+_2018_Icarus/Granvik+_2018_Icarus.readme)*:\n",
    "\n",
    "Version 1.0\n",
    "\n",
    "File Granvik+_2018_Icarus.dat contains a realization of the NEO model\n",
    "described in\n",
    "\n",
    "Granvik, Morbidelli, Jedicke, Bolin, Bottke, Beshore, Vokrouhlicky,\n",
    "Nesvorny, and Michel (2018). Debiased orbit and absolute-magnitude\n",
    "distributions for near-Earth objects. Accepted for publication in\n",
    "Icarus.\n",
    "\n",
    "The file Granvik+_2018_Icarus.dat contains 802,000 synthetic NEOs with\n",
    "17<H<25 parameterized by their orbits and absolute magnitudes. The 7\n",
    "columns in the file are\n",
    "\n",
    "1) semimajor axis [au]\n",
    "2) eccentricity\n",
    "3) inclination [deg]\n",
    "4) longitude of ascending node [deg]\n",
    "5) argument of perihelion [deg]\n",
    "6) mean anomaly [deg]\n",
    "7) absolute magnitude\n",
    "\n",
    "Please cite the above article if you use the data in your\n",
    "publications.\n",
    "\n",
    "Mikael Granvik\n",
    "27-04-2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5b21a18d-7e06-44a4-b541-d134cb3fa25e",
   "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>SemiMajorAxis_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>H</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.574981</td>\n",
       "      <td>0.783617</td>\n",
       "      <td>33.520763</td>\n",
       "      <td>278.480591</td>\n",
       "      <td>75.952057</td>\n",
       "      <td>103.833748</td>\n",
       "      <td>21.064367</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.682046</td>\n",
       "      <td>0.275324</td>\n",
       "      <td>25.924536</td>\n",
       "      <td>169.683060</td>\n",
       "      <td>350.926636</td>\n",
       "      <td>46.401562</td>\n",
       "      <td>24.683197</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2.480214</td>\n",
       "      <td>0.546448</td>\n",
       "      <td>14.240792</td>\n",
       "      <td>121.685623</td>\n",
       "      <td>316.509247</td>\n",
       "      <td>279.612305</td>\n",
       "      <td>24.843449</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.941440</td>\n",
       "      <td>0.620780</td>\n",
       "      <td>13.679822</td>\n",
       "      <td>8.245143</td>\n",
       "      <td>175.087372</td>\n",
       "      <td>129.364609</td>\n",
       "      <td>24.954794</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.114843</td>\n",
       "      <td>0.548071</td>\n",
       "      <td>37.906609</td>\n",
       "      <td>209.261627</td>\n",
       "      <td>37.110592</td>\n",
       "      <td>120.817390</td>\n",
       "      <td>24.342936</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",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>801995</th>\n",
       "      <td>2.877656</td>\n",
       "      <td>0.607120</td>\n",
       "      <td>6.684607</td>\n",
       "      <td>118.446182</td>\n",
       "      <td>39.764286</td>\n",
       "      <td>265.229523</td>\n",
       "      <td>23.330488</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>801996</th>\n",
       "      <td>1.563997</td>\n",
       "      <td>0.207188</td>\n",
       "      <td>28.263660</td>\n",
       "      <td>135.846954</td>\n",
       "      <td>82.046822</td>\n",
       "      <td>101.014885</td>\n",
       "      <td>24.990154</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>801997</th>\n",
       "      <td>2.679385</td>\n",
       "      <td>0.535874</td>\n",
       "      <td>7.462100</td>\n",
       "      <td>227.803268</td>\n",
       "      <td>272.786926</td>\n",
       "      <td>153.444550</td>\n",
       "      <td>23.996401</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>801998</th>\n",
       "      <td>1.592931</td>\n",
       "      <td>0.826241</td>\n",
       "      <td>43.307148</td>\n",
       "      <td>337.204224</td>\n",
       "      <td>43.215237</td>\n",
       "      <td>107.381775</td>\n",
       "      <td>24.795914</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>801999</th>\n",
       "      <td>2.591117</td>\n",
       "      <td>0.826896</td>\n",
       "      <td>22.001493</td>\n",
       "      <td>334.001678</td>\n",
       "      <td>286.224030</td>\n",
       "      <td>218.469681</td>\n",
       "      <td>24.542505</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>802000 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        SemiMajorAxis_AU      Ecc_   Incl_deg  LongAscNode_deg    ArgP_deg  \\\n",
       "0               2.574981  0.783617  33.520763       278.480591   75.952057   \n",
       "1               1.682046  0.275324  25.924536       169.683060  350.926636   \n",
       "2               2.480214  0.546448  14.240792       121.685623  316.509247   \n",
       "3               1.941440  0.620780  13.679822         8.245143  175.087372   \n",
       "4               1.114843  0.548071  37.906609       209.261627   37.110592   \n",
       "...                  ...       ...        ...              ...         ...   \n",
       "801995          2.877656  0.607120   6.684607       118.446182   39.764286   \n",
       "801996          1.563997  0.207188  28.263660       135.846954   82.046822   \n",
       "801997          2.679385  0.535874   7.462100       227.803268  272.786926   \n",
       "801998          1.592931  0.826241  43.307148       337.204224   43.215237   \n",
       "801999          2.591117  0.826896  22.001493       334.001678  286.224030   \n",
       "\n",
       "        MeanAnom_deg          H  \n",
       "0         103.833748  21.064367  \n",
       "1          46.401562  24.683197  \n",
       "2         279.612305  24.843449  \n",
       "3         129.364609  24.954794  \n",
       "4         120.817390  24.342936  \n",
       "...              ...        ...  \n",
       "801995    265.229523  23.330488  \n",
       "801996    101.014885  24.990154  \n",
       "801997    153.444550  23.996401  \n",
       "801998    107.381775  24.795914  \n",
       "801999    218.469681  24.542505  \n",
       "\n",
       "[802000 rows x 7 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "granvik_model_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dde92e37-a913-453d-977c-c6e90f59e1fc",
   "metadata": {},
   "source": [
    "## Enrich the dataframe\n",
    "\n",
    "We need to convert some parameters and derive other orbital elements for SPICE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b9b1e757-0843-4885-8d39-c5c969629d19",
   "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": 9,
   "id": "dd00bf74-c055-4e21-bf4d-66a16798ca7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the perihelion and convert the angular values to deg\n",
    "granvik_model_df.loc[:, \"Perihelion_AU\"] = (1.0 - granvik_model_df[\"Ecc_\"]) \\\n",
    "                                               * granvik_model_df[\"SemiMajorAxis_AU\"]\n",
    "\n",
    "granvik_model_df.loc[:, \"Incl_rad\"] = np.radians(granvik_model_df[\"Incl_deg\"])\n",
    "granvik_model_df.loc[:, \"LongAscNode_rad\"] = np.radians(granvik_model_df[\"LongAscNode_deg\"])\n",
    "granvik_model_df.loc[:, \"ArgP_rad\"] = np.radians(granvik_model_df[\"ArgP_deg\"])\n",
    "\n",
    "# Convert AU to km\n",
    "granvik_model_df.loc[:, \"Perihelion_km\"] = \\\n",
    "    granvik_model_df[\"Perihelion_AU\"].apply(lambda x: spiceypy.convrt(x, \"AU\", \"km\"))\n",
    "granvik_model_df.loc[:, \"SemiMajorAxis_km\"] = \\\n",
    "    granvik_model_df[\"SemiMajorAxis_AU\"].apply(lambda x: spiceypy.convrt(x, \"AU\", \"km\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6593656f-f138-4766-8fa1-d97dac1a11d1",
   "metadata": {},
   "source": [
    "## Generate \"generic\" positions\n",
    "\n",
    "To compute a time-invariant, normalized distribution on NEOs in the sky, we need to populate all orbits with (theoretically), an infinite number of objects. Then we can compute their corresponding brightness (apparent magnitude) and consequently their observability along their orbit. However, theoretically, we also need to move Earth along its orbit too and normalize a large number of these brightness distributions.\n",
    "\n",
    "However, we simplify a few things. Feel free though, to extend this project :)!\n",
    "\n",
    "*So are we going to do?*\n",
    "\n",
    "- Of course we cannot create a infinite number of NEOs. Thus, we populate each orbit with 36 NEOs in an angular distance of 10 degrees. On a circular orbit, this makes sense ... but how does it look like on an elliptical orbit? The Mean Anomaly does not correspond to some kind of \"actual\" position of the NEO, but is a projection on a theoretical, circular construct (take a look at Wikipedia: https://en.wikipedia.org/wiki/Mean_anomaly). But there is also something called true anomaly. Like the \"actual\" position of the NEO. Which one shall we take? Let's take a look now. What we need: a positional parameter that somehow considers the \"dwell time\" of an object (e.g., an object at the Perihelion has a shorter dwell time than far away at the Aphelion)\n",
    "\n",
    "- Additionally, we do not want to exceed computational time. So let's put our home planet on a certain position (we do not integrate over its entire orbit). Feel free though to extend this model in this regard!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5722f76-b2b7-49a9-ad29-79ae41bcfedb",
   "metadata": {},
   "source": [
    "### Sample Orbit\n",
    "\n",
    "Let's take a look at some kind of sample orbit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "91b87fc3-5c49-4ec8-8205-87e4d955ab77",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, we need to define a function that converts true anomaly values to mean anomalies (since\n",
    "# SPICE does not handle true anomalies\n",
    "def true2mean(true_anom, ecc):\n",
    "    \n",
    "    atan2_x1 = -1.0 * np.sqrt(1.0 - ecc**2.0) * np.sin(true_anom)\n",
    "    atan2_x2 = -1.0 * ecc - np.cos(true_anom)\n",
    "    \n",
    "    mean_anom = np.arctan2(atan2_x1, atan2_x2) \\\n",
    "                + np.pi \\\n",
    "                - ecc * ((np.sqrt(1.0 - ecc**2.0) * np.sin(true_anom)) \n",
    "                         / (1.0 + ecc * np.cos(true_anom)))\n",
    "    \n",
    "    return mean_anom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2257dc30-6c9d-43c8-b7fc-7a3aee1d9223",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's create an example plot with the true anomaly distribution\n",
    "# Sample orbit:\n",
    "sample_rp = spiceypy.convrt(0.5, \"AU\", \"km\")\n",
    "sample_ecc = 0.7\n",
    "sample_incl = 0.0\n",
    "sample_lnode = 0.0\n",
    "sample_argp = 0.0\n",
    "sample_et = spiceypy.utc2et(\"2000-001T12:00:00\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c6c1190b-9b4f-4d24-a2c6-94c86416ea79",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████| 3600/3600 [00:00<00:00, 57773.00it/s]\n"
     ]
    }
   ],
   "source": [
    "# First we create a \"complete\" orbit that is used to explain our method:\n",
    "mean_anomaly = np.radians(np.arange(0, 360, 0.1))\n",
    "\n",
    "# List that will store the positions in 2D (X-Y plane)\n",
    "sample_trajectory = []\n",
    "\n",
    "for k in tqdm.tqdm(mean_anomaly):\n",
    "\n",
    "    # Compute the state vector and store in positional values in X-Y direction, converted in AU\n",
    "    temp_state = spiceypy.conics([sample_rp, sample_ecc, sample_incl, sample_lnode, sample_argp,\n",
    "                                  k, sample_et, gm_sun], sample_et)\n",
    "    sample_trajectory.append([spiceypy.convrt(temp_state[0], \"km\", \"AU\"),\n",
    "                              spiceypy.convrt(temp_state[1], \"km\", \"AU\")])\n",
    "    \n",
    "sample_trajectory = np.array(sample_trajectory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3c05e39b-60b9-425e-a868-84ca832b48ff",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████████| 36/36 [00:00<00:00, 26905.73it/s]\n"
     ]
    }
   ],
   "source": [
    "# Now we create 36 NEO positions using the mean anomaly. Basically ... we need the same code again!\n",
    "mean_anomaly = np.radians(np.arange(0, 360, 10))\n",
    "\n",
    "# List that will store the positions in 2D (X-Y plane)\n",
    "mean_anom_10deg_trajectory = []\n",
    "\n",
    "for k in tqdm.tqdm(mean_anomaly):\n",
    "\n",
    "    # Compute the state vector and store in positional values in X-Y direction, converted in AU\n",
    "    temp_state = spiceypy.conics([sample_rp, sample_ecc, sample_incl, sample_lnode, sample_argp,\n",
    "                                  k, sample_et, gm_sun], sample_et)\n",
    "    mean_anom_10deg_trajectory.append([spiceypy.convrt(temp_state[0], \"km\", \"AU\"),\n",
    "                                       spiceypy.convrt(temp_state[1], \"km\", \"AU\")])\n",
    "    \n",
    "mean_anom_10deg_trajectory = np.array(mean_anom_10deg_trajectory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "7b257054-e15b-4793-90c2-835a3b4d9cf6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████████| 36/36 [00:00<00:00, 22641.32it/s]\n"
     ]
    }
   ],
   "source": [
    "# And finally a list with the true anomaly values!\n",
    "# Now we create 36 NEO positions using the mean anomaly. Basically ... we need the same code again!\n",
    "true_anomaly = np.radians(np.arange(0, 360, 10))\n",
    "mean_from_true_anom = [true2mean(n, sample_ecc) for n in true_anomaly]\n",
    "\n",
    "# List that will store the positions in 2D (X-Y plane)\n",
    "true_anom_10deg_trajectory = []\n",
    "\n",
    "for k in tqdm.tqdm(mean_from_true_anom):\n",
    "\n",
    "    # Compute the state vector and store in positional values in X-Y direction, converted in AU\n",
    "    temp_state = spiceypy.conics([sample_rp, sample_ecc, sample_incl, sample_lnode, sample_argp,\n",
    "                                  k, sample_et, gm_sun], sample_et)\n",
    "    true_anom_10deg_trajectory.append([spiceypy.convrt(temp_state[0], \"km\", \"AU\"),\n",
    "                                       spiceypy.convrt(temp_state[1], \"km\", \"AU\")])\n",
    "    \n",
    "true_anom_10deg_trajectory = np.array(true_anom_10deg_trajectory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "85524f3b-c384-48ca-97e1-94b66ece93c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.125, 0.5, 'Eclip. y in AU')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAJCCAYAAADpz40lAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAACOqElEQVR4nOzdd1gU19cH8O/usiBVEBAUEOzYW7ArWLAmViQWjBUVE5NoNGoSaxJNNLG3SBRbNCYaY4wxVrCXxZjYxQKi0ovUBZblvH/wMj9XipRdZoHzeZ7z6MxOOTvA2dk7d+5IABAYY4zpLanYCTDGGCsaF2rGGNNzXKgZY0zPcaFmjDE9x4WaMcb0HBdqxhjTc1yoGdMhIsKiRYtE27+zszOICOPGjRMtB1Z2XKgriHHjxoGIQETo2rVrgcs8fPgQRITAwMByzq70du7cCSLCli1bxE6lQsj7HSAiqFQqxMfHIzg4GGvWrEGTJk20tp9Ro0bho48+0tr2WNkRh/7HuHHjiIgoPT2dNm3alO/1Dh06CK8HBgaKnm9xwtjYmJKTk+nJkycUHx9Pcrlc9Jy0HUREixYt0ur2Tp06RWPGjCEfHx96//33aevWrZSYmEhZWVk0c+bMfOsYGRmRVCot0X6OHDlCoaGhoh8/DiFET4CjGJFXqA8cOECxsbFkYGCg8fratWvp7t27dOvWrQpTqEeNGkXZ2dnUrVs3IiIaMmSI6DlpO3RRqDdv3pxvfo0aNejixYtERNS/f/8y74cLtX4FN31UMPv27UONGjXQt29fYZ5UKsW7776LvXv3FrreBx98gJs3b0KpVCI6Oho//vgjrK2tNZZ555138Mcff+DZs2fIyMhAWFgYVqxYASMjI43lAgICoFQqUbt2bRw6dAgpKSmIiYnBypUrIZUW/1fKx8cHQUFBOH/+PIKDgzFmzJh8y7i7u4OIMHLkSMyfPx/Pnj2DUqnEqVOnUL9+/XzLDxs2DAqFAunp6YiLi8PevXvh5ORUYP5OTk44cuQIUlJS8OLFC8yYMQMA0LRpU5w6dQqpqal4+vQpfHx8NNa3srLCihUr8N9//yE5ORkpKSkIDAwstEkqT4MGDUBEmDlzZr7XWrRoASLCtGnT3njcCpKQkICRI0dCpVLh888/F+YX1EZtamqKlStX4smTJ8jIyEBMTAwCAwPRrVs3AEBgYCDefvttuLi4aDS15Jk1axbOnz+P2NhYKJVK3Lx5E5MmTcqXU2hoKI4dO4YuXbrg6tWrUCqVePz4McaOHZtvWQsLC6xYsQKPHz9GRkYGnj9/jp9++gm1a9cWljE0NMTChQsREhIiLLNq1SoYGxtrbKtnz544e/YsEhISkJaWhkePHmH9+vWlOq76RPRPC443R94ZdYcOHSgoKIh++ukn4bU+ffoQEVG9evUKPKPevHkzZWVlkb+/P02ZMoW+/PJLevnyJd28eZOMjIyE5X777Tf69ddfac6cOTRlyhTavn07qVQq2rdvn8b2AgICKDMzk27evEk//vgjTZ06lX799VciIpo2bVqx3o+trS1lZWXRpEmTCADNmjWLlEolWVhYaCzn7u5ORET//PMPKRQK+uijj2jhwoWUmppKV65c0Vh2zJgxRETCcl999RUlJydTeHg4WVtb58v/1q1btGXLFpo2bRqdP3+eiIjGjx9Pz549oxUrVtD7779Pt27douzsbKpfv76wfrt27ejJkyf0zTff0JQpU+jTTz+lBw8eUEZGBrVo0UIjp9fPqC9evEjXr1/Pdzy+/fZbysjIICsrqyKPW2Fn1Hlx8uRJys7OJnNzcwJAzs7OREQ0btw4YZldu3ZRRkYGff/99zRx4kT65JNP6JdffqEPPviAAFDv3r3pn3/+oZiYGBozZowQees/e/aMNm/eTDNmzKAPPviATp48SUREU6dO1cglNDSU7t+/T5GRkfTVV1/R9OnTKTg4mNRqNTVt2lRYzsTEhG7cuEEqlYr8/f1p6tSp9Nlnn9GFCxeoe/fuwnJHjx6l1NRUWr16Nfn6+tKqVasoIyODjh8/LizTpEkTUiqVdO3aNfrwww9p8uTJtHz5cvrnn39E/xsuY4ieAEcx4tVCPWXKFEpNTSUTExMCQDt37qTLly8TgHyFulOnTkRENHbsWI3tdenShYiIfH19hXnGxsb59jt//nxSq9Xk6OgozAsICCAiogULFmgse/36dVIoFMV6PzNmzKDMzEyytLQkAOTg4EBqtVoo3HmRV6jv3r2r0YY9Y8YMIiJq1qwZASADAwOKjIykO3fuULVq1fKtv3Llynz5f/HFF8K86tWrU1paGqnVaho9erQwv1GjRkRE9OWXXwrzDA0NSSKRaORpaWlJUVFR5O/vrzH/9UI9depUIiJq0qSJME8ikVB4eDgdPHjwjcftTYV69erVRETCB0ZBhTohIYHWr19f5H6Kavoo6Pfk+PHj9PDhQ415oaGhRETUrVs3YZ6NjQ0plUqNn8eiRYuIiMjLy6vQfEaNGkVqtVqjcAOg0aNHExGRp6cnAaAPP/yQiEjjg7kyBDd9VEC//vor5HI5hgwZgmrVqmHIkCH46aefClzW29sbKSkp+Pvvv2FtbS3E/fv3ERUVhR49egjLKpVKAIBEIoGFhQWsra1x4cIFSKVStG3bNt+2/f39NabPnz+PevXqFes9+Pj44Pjx43j58iUA4MWLFzh//nyBzR8AsGvXLqhUKo19ARD299Zbb8He3h6bN29GRkaGsNzZs2cRHByMgQMH5tvmjz/+KPw/KSkJDx48QEZGBvbt2yfMDwkJQWJiosb7ysrKEpoCjIyMUKNGDchkMigUCrRr167I971//35kZGRofP338PCAk5MTdu/eXeS6xZGamgoAMDc3L3SZpKQkdOjQQaNZoSTyfk8MDAxgZWUFa2trBAYGokGDBrCwsNBY9sGDB8LPCgDi4uLw4MEDjePp5eWF27dv48CBA4Xu09vbGyEhIbhz547G7/HZs2eRk5Mj/B4nJSUBAIYMGQKJRFKq96ePuFBXQImJiTh+/DjGjBmDQYMGwcTEBPv37y9w2UaNGsHc3BwxMTGIi4vTCHt7e9SsWVNYtlmzZjh69ChSU1ORlJSEuLg4nDt3DgBQvXp1je1mZWUhKioqX141atR4Y/4NGzZE+/btceHCBdSvX1+Is2fPwt3dHQ4ODvnWCQ8Pz7cvILe9GMhtiwVyC8Pr7t27BxcXlzfmn5SUhBcvXmi0x+bNz9sPkPtBNnfuXKE9NT4+HnFxcXj77bfzHafXvXz5En/88QdGjx4tzPPx8UF8fDyOHj1a5LrFYWZmBgBISUkpdJk5c+agWbNmCA8Ph0KhwJdffolGjRoVex+DBg2CQqGAUqlEQkIC4uLisHz5cgD5f09e/7kBuT+7V49n/fr1cfv27SL32ahRI7i6uub7HX7+/DmkUqnwe7x//36cP38eP/74I2JiYrB//36MGjUKMpms2O9PHxmInQArnb1792LXrl2wsLDAqVOnEBsbW+ByUqkUcXFxGDlyZIGv5xU8CwsLBAYGIi0tDZ9//jkePXoEpVIJBwcH7Ny5M99FwpycnFLnnndx7ttvv8W3336b7/UxY8ZgxYoVGvPUanWB2yrtWVNh+RdnP/Pnz8fXX3+NHTt24IsvvkB8fDzUajXmz59f4AXO1+3atQve3t7o3r07rl69iuHDh2Pv3r0a3xhKq3nz5sjOzkZoaGihyxw4cADnz5/H4MGD0adPH3z44Yf49NNPMX78eI1vEwXp0qULDh06hAsXLmDatGmIiIhAVlYWBgwYgFmzZuX7PdHWz00qleLOnTuF9u2OiIgAAGRkZMDd3R3du3fHgAED0LdvX+zduxezZs1Ct27dNL5tVSRcqCuow4cPIzMzE127dsV7771X6HKPHz+Gp6cnrly5grS0tEKX69GjB2xtbeHl5SWcRQNA7969tZo3kFuIg4KCsGHDhnyvzZo1q8BC/SZPnz4FADRu3BgnT57UeM3V1RVhYWGlzvd1I0aMQGBgICZMmKAxf8mSJcVa/++//0Z0dDTGjh0LOzs7VK9eXSvNHk5OTnB3d8fly5eFJpDCREdHY+vWrdi6dSuqV6+OK1euYMmSJUKhfv1bRR4vLy9kZGSgT58+yMzMFOa/2oRWUo8fP0bz5s3fuEy7du1w+vTpN26PiHD27FmcPXsWc+fOxbRp07B582YMGzasyJ5R+oybPioopVIJPz8/LF68GIcOHSp0uf3790Mmk2HhwoX5XpNKpbC0tATwvzOfV890JBIJZs2apdW8O3XqhPr162Pnzp04ePBgvti1axdatmyJFi1alGi7wcHBiIqKwtSpUzW6E3bt2hVubm74888/tfYe1Gp1vjPCTp06oVOnTsVe/6effoKXlxcmTZqEhw8f4vLly2XKycrKCvv27YNMJsPXX39d6HJSqTRfO3JSUhJCQ0OF3wUASEtL05h+NXci0jhztrS0xMSJE0ud+4EDB9C8eXN4eXkVusz+/fthb28PPz+/fK8ZGhoKTT4FNb39888/Qp4VFZ9RV2B79ux54zLnz5/Hhg0b8Omnn6Jly5Y4fvw4MjMz0aBBA3h5eWHhwoXYuXMnLl68iLi4OOzcuRPr16+HSqWCl5eX8AegLT4+PlCr1YUWziNHjmDTpk3w8fHB3Llzi73d7OxszJkzB7t378b58+exZ88e2Nra4sMPP8Tz588LbGIprT/++ANLlizBzp07cf78eTRs2BBTpkzB3bt3i328du3ahVmzZqFv374lHgukYcOGGDNmjHDRt1WrVhgxYgTMzMwwa9YsHD9+vNB1zc3N8eLFCxw8eFDoB96lSxf0799fo69xcHAwRo4ciTVr1uDq1avIycnB/v37ceTIEXzyySc4efIkdu/ejRo1asDX1xdRUVGoVatWid5HnpUrV2L48OHYt28f+vTpg+vXr8PS0hL9+/fHwoULce7cOezZswdeXl7YtGkT3N3dceHCBUgkEjRu3Bje3t4YMWIEzp49iwULFsDDwwNHjx5FWFgYrKysMG3aNKSmpmr1w1oMonc94XhzvNo9r6jlCrszccKECXT16lVKS0ujpKQkunXrFq1cuZKcnJyEZdq3b0/nz5+n1NRUioqKoo0bN1Lz5s3zde8KCAggpVKZbx953awKy83AwIBiY2Pp4sWLRb6Ha9euUXh4OEkkEqF73bvvvquxTEHdzgDQsGHDSKFQkFKppPj4eNq3bx/VqVNHY5nC8g8MDKR79+7lmx8aGkrHjh0TpuVyOX3zzTf07NkzSk9Pp2vXrlHfvn0pICAgX5e217vnvRr//vuv0P+9uL8Hr8rOzqaEhAS6fv06rV69WqPLX2HHSS6X07fffkv//PMPJSYmUmpqKt26dYtmzZpFMplMWM/Y2Jh27NhBcXFxpFarNX6uY8eOpbt375JSqaSHDx/SzJkzafz48URE5OzsXOhxe/U4v/47amlpSWvWrKHw8HDKzMykZ8+e0Z49e6hWrVrCMjKZjD755BO6efMmKZVKSkhIoODgYFq8eLHQ/9zDw4MOHjxIz549o4yMDHrx4gUdOHCAWrZsKfrfcFlC8v//YYyVs2vXriErK+uNdzQyxm3UjImgVatWcHNzw86dO8VOhVUAfEbNWDlq1qwZ2rVrh5kzZ6JWrVqoW7eucAMJY4XhM2rGypGXlxcCAgJgbGyMkSNHcpFmxcJn1Iwxpuf4jJoxxvQc96P+fzExMcLdbYwxVt6cnZ01xt55FRfq//f06VO4ubmJnQZjrIpSKBSFvsZNH4wxpue4UDPGmJ7jQs0YY3qOCzVjjOk5LtSMMabnuFAzxpie40LNGGN6jgs1Y4zpOS7UjDGm57hQM8aYnuNCzRhjeo4LNWOM6Tku1Iwxpue4UDPGmJ7jQs0YY3qOCzVjjOk5LtSMMabnuFAzxpie40LNWBnJzGRwmOgAmalM7FRYJcWFmlUpuiiqVu5WMG9jDisPK61t83X8YVC1caFmVcrrRbVatWqwsvpfge3Tpw969eolTH/00UcYP368ML169Wp88sknwvTW3VvRYHADZCdlw6qrFfx3+2PWrFnC6z/88AMmT54sTH/++efo37+/MD1gwAA0atRImDY1NS1W3qxq4ULNKjVLS0vUrVsXQO5ZqfMAZ1Q3qA6rrlaQmcqwdu1afPvtt8LyXl5eGD58uDBtb28PW1tbYfrp06eIiIgQpn979BuUmUqQigApcDjsMO7fvy+8/uzZM8THxwvTTZs2hb29PQBAIpFg0aJFwgeDRCJBYGAgpk6dKkx//vnnaNulLay6WiE7KRvW3a35rLoKkgAgsZPQBwqFAm5ubmKnwQogM5PB3tseUfujoE5TF7lst27d0KJFC2zatAkA8OWXX6Jly5YYPHgwbAbaoO6gujCRmiAiIQKJ5xLRJLkJpFIpzp49CyD3jDYzMxPZ2dnFyqve5/WQk5UD5ACQAlJDKZ589eSNeeaxtbWFSqXCy5cvYWBgAG9vb9y8eRO3b9+GlZUVfv75Zyw7vQwR1hGQZklRu15t5PyXg7MrzsLc3BzDhg3D6dOn8fz582Ltj+mvomoQn1Ezvffq136JRAJHR0fhtREjRuDYsWOQSnN/lZs3b46BAwdCIpEAAH7++WcsX74cMjMZrLpaISEqAc+fP0dORg6sulrh0j+XhCINAGlpacUq0nl5QYrcIg0IxbokzROxsbF4+fIlACA7Oxt79+7F7du3AQCJiYkYMHwAIqwikJORgxzKQXJ8MmJsYyAzlcHFxQUffPABnJycAAANGzbE5s2b0aBBAwCAgYGBcBxYxcaFmum1ajWqoVbvWlAnq2HV1QqDvAfh0KFDsLOzAwA8f/4cZ8+ehbGxMQBgy5YtGDhwIIhyvyjeuXMHV65c0UpRfZ2xizEkMglk5jIhJDIJjF2My/CONb2ad446By/jXyI7JxtWHla4desWPDw8EBwcDAAwMTGBsbEx0tPTAQA9e/ZEUFCQUMjNzc1RrVo1reXGyo+B2Akw9iq5XA43Nzc8evQIMTExaOLdBJY1LRH/Ih5KUuKW9Ba+/PJLoRhdvnwZly9fFtbPK9Cve7Wovj6/tJ5tfFbqdYvrTXmnpaUJ8/777z+NC5/Pnj3Dn3/+iRcvXgAARo4ciQkTJqBnz57IyMiAjY0NUlNTkZGRobHtkjQ1sfLBhZqJrmHDhsjIyMCzZ89gY2ODtWvXYuXKlTj410GkuaQhOzIbSqUSBIKspQxHvzpa4gJSHkVVF8qS971793Dv3j1h+uLFi0hOThYK88cff4yWLVti0KBBAABra2skJCQITU2Z0ZmIOxpXtjfAtIKbPphOFdT/Vy6XCz0fDAwMsG3bNowePRoAEBkZiSlTpuD333+HlbsVSEJIT0kH5ZBWmiuqsrt372L//v3C9MGDB7FhwwZhes2aNVi+ZrnQwySvZwwTH59RM53KOzvLislC7J+xAIAff/wR6enp8PPzQ3Z2NmbPno0nT54I69y4cQOAbpor2P/kHec8e/bsQXyDeKAOQJmE2nVqw3GGIwK/CRQpQ5aHCzXTmbyeFmYSM9h72+N84Hmo09TYtm0bMjMzheWuXbtW4PoVtbmiojp18RTqedRDTkYOpFIpVGkqJNolQmYqg2U1S8ybNw+bN2/W+FBl5YObPpjWOTs7Y/78+ajZuyYgBVQZKqQp02DbO/fGkXPnzuHq1asiZ8le93oPk5ioGGRkZcDKwwp16tRB8+bNoVbnXhuoU6cOOnTowN3/yoleF+pu3brh8OHDeP78OYgI48aNe+M6zZs3R1BQENLT0/H8+XMsWLCgHDJlNjY2sLS0BJDbDayrZ1dYd7dGTkYOUlNSkRCVAPOO5tzmqceK6m5448YNDBgwAE+fPgUADB8+HKtWrYKZmRkAcMHWMb1u+jAzM8Pt27exa9cu7Nq1643Lm5ub4+TJkzh37hzc3Nzg6uqKgIAApKWlYdWqVeWQcdVkbm6OP/74A7t378bmzZtx+/ZtTFo7CeadzQvst8w9CfTTm5qaXu36uGHDBpw4cQIpKSkAgO+++w5xcXFYvny5TnOsqvS6UB87dgzHjh0DAOzYseONy48ZMwYmJiYYN24cMjIycOfOHbi6umLWrFlcqMugoH61w4cPR61atbBhwwakpKRg+fLl+Oeff4R15I5yvhBYialUKty5c0eYfvjwIZKSkoTpd955B4GBgUhNTRUjvUpHrwt1SXXq1Annz5/X6MB//PhxfPXVV3BxcUFYWJh4yVVgeT03pMlSPDuQe9ZVr1491KlTR1jmyJEjGuvwhcCqZcuWLcL/GzVqhIULF8LAwACHDh0SMavKo1IVant7+3yD00RHRwuvcaEuubyeG0bZRnAY5gCDuwYIvRuKVatWCReWGHtVSEgIfHx8EBoaCgDo1asX+vbtiy+//FJoKmElU6kKdUn5+vpiypQpAHIvhjFNTZo0QY5bDnKkOUhPSUeScRLUzdTAXXCRZkV68OCB8H8zMzNYW1sLzSDm5uZcsEtIr3t9lFRUVJQwWE+evOmoqKh8y/v7+8PNzQ1ubm6Ii+MLXK+Sy+X4atVXQs+NHHUOEqMTYdTWiHtusBI5fPgwJk2aBCKCTCbDnj17MHPmTLHTqlAqVaG+fPkyunXrBiMjI2Gep6cnXrx4wc0exeDi4oLp06cDyL1YtOjXRYhPiNfqiHOsapNKpfj1119x8eJFALknBPXr1xde50eOFUyvC7WpqSlatWqFVq1aQSqVok6dOmjVqpUwbOOyZctw6tQpYfm9e/ciPT0dO3bsQLNmzTB06FDMmzePe3wUU5s2beDt7Q0HBwcAQJJxEkhKOh3Gk1UtKpUKe/bsEe5G9fLywr59++Di4gKAHzlWFNLXcHd3p4IEBAQQAAoICKDQ0FCNdZo3b05nz54lpVJJERERtHDhwmLtS6FQiP5+yyNkZjJymOhAMlMZGRoa0rRp06h79+4EgKRSKVWvXl30HDmqTpibm9OQIUMIyP3dbL6qOTX8uiE1XN6QZKYy0fMrzyiqBun1xcSzZ88WecfThAkT8s27ffs23N3ddZlWhfbqEJYJxxLg7u4OuVyOc+fOIScnR6MvLGO6lpKSgt9//x0AULN3TVjbWiMnPQdRL6P45qhX6HWhZtolM5PB1sMWloaWUHdVIzEoEePHj9cYIIkxMcjMZDDrYIa4iDios9XIyc6BdXdrWD2zwsObD8VOT3R63UbNtMvK3QoGhgawNLeEYTVDWHlYcZFmeiFvQKiM9AyoslRADmBhaYF3v34X1tbWYqcnOj6jrgJat24NaydrhHcOR/rLdDx6+Qg5yH24a2JQIj9uiYmuoLHHU1JTcDT4KOLj4wEAtWrVQmRkpFgpiooLdRUwYcIEPKzxEOER4blDWOb1t+NBkpieeNOQAw0bNsTu3buxdOlS/PXXX+WUlf7gpo9KqnXr1rCwsAAALF68GH8q/tT5E7MZ05Xw8HBs374d58+fB1A1h1UVvVuKPkRl6p5Xo0YNunjxIn300Uei58LBoe2QSCS0fv168vX1FT0XbUaF7Z7HiidvGNKXf7xEWlwaEhIS8Omnn2oMO8pYZWFgYICIiAjExsaKnUq5Ev2TRB+iIp9R2wy0oeabm9Pc3+ZS8+bNRc+Hg6M8o3379tSnTx/R8yhr8Bl1JZY3DGlWQhbOGZxDgjJB7JQYK1cjRoyAk5MTTp8+XWlHdeRCXYFZWlqiy6wuCJGGQJ2hRlxCHLJcs4DHYmfGWPmZN28eLC0toVarIZPJUL16dSQkVK4TFu71UYF17NERGfUzICc5ACAnI7dvNI88xqoStVot9LWePn069u7di+rVq4uclXbxGXUFZGxsDKVSieCsYNhG2SJLmZX7Aj9AllVxf/75J16+fFnpxqzhM+oKZuDAgTh06BBq1aoFYxdjqEnNfaMZ+3+hoaHYvXs3AMDR0RGjRo0SOSPt4DPqCubWrVu4dOkSEhISkLmRx+lgrDBDhgzB0KFDcfz48UrRZi16txR9CH3tniczk5HjZEfq3re76LlwcFSkkEgkVLt2bY1psXMqKoqqQdz0oees3K1g3d4aPT7qgcaNG4udDmMVBhEhIiICADB06FCsW7cO1apVEzmr0uGmDz2W10c6LSYNv6X/hkfPH4mdEmMVkkqlQmZmJlQqldiplAqfUespY2NjDFk4BFKZFKQiZGRm8HPkGCulP//8E7Nnz4ZarYaRkRFMTEzETqlEuFDrqSZtm+Clw0vIJdxHmjFt+u6777B69eoKNQIfN33oqXDLcKRHpkOV+v9f1biPNGNacejQIcjlchCR2KkUGxdqPSKVSrF48WIcOXIEUS5RyEGOxhMvAHAfacbK6MyZM8L/nZ2d8fz5c70fI4QLtR6xsLBAw4YN4eLiAsVGhdjpMFapWVtbY9euXThw4ADWr18vdjpvJHr/QX0IfelHLZfLRc+Bg6OqxLBhw8jW1paA3HsWHCY6kMxUJkou3I9aj8nMZHBb4IaZn82EVCqtsN2HGKuIfvvtN+EBBA2HN4R5G3O97F3FhVpkVu5WUNVRIbVeKmQy7tHBmBiGjh4Khz4OkKZL9bJ3FRdqEeXd0JIanYqg2CDkGOaInRJjVdJV5VUkJCUgMz1T6F2lT7hQi0Qmk2HY4mEwMDQAqUgvfzkYqwpkZjKYtjdFcmxy7gwV9O6smgu1SJwaOSHOLg5GEiMAfEMLY2KxcrfKrYQ5gEQqQV3nurC0ttSrEycu1CJJrZeKiKgIpCSl5M545YYWxlj5MXYxhkQmgcxcBqmpFLGpsVAqlXp1z0KR/aiTk5MLvHsnKSkJDx48wIoVK3Dy5EmdJQcAfn5+mDNnDmrVqoU7d+7g448/xoULFwpc1t3dHUFBQfnmu7q64sGDBzrNs7isra3RqVMn/Of8H0hCfEMLYyJ7tvFZoa9JJBK9uIOxyEL9wQcfFDjf0tIS7dq1w5EjR+Dl5YU///xTJ8l5e3tj7dq1mD59Oi5cuIDp06fj2LFjaNq0KZ49K/zgNm3aVGOg8LzuN/rA29sbY8aMwbBhw/AspvD3wBgT16hRo9C9e3dMnz5dL4p1qTtof/zxx3Tp0iWddQC/cuUKbd26VWNeSEgILVu2rMDl3d3diYjI2tpaq53NtRkSiYQaNmwoSod6Dg6O4kf//v3p66+/JmNj43LZ3xtqUOk33LBhQ0pISNBJ0nK5nFQqFXl5eWnM37BhAwUFBRW4Tl6hDg0NpYiICDp16hR5eHho4yCVKWRmMqrrV5dMbUzL5QfOwcFR8UJndyZWq1YNGRkZZdlEoWxsbGBgYIDo6GiN+dHR0bC3ty9wncjISEybNg3Dhw/HsGHD8ODBA5w+fRpdu3YtcHlfX18oFAooFArY2Nho/T3ksXK3Qg23Gpi4ZmKFGlqRMQbY2tpi2rRpov7tlmlQpsmTJ+Pff//VUiplFxISgpCQEGH6ypUrcHFxwZw5cwq8AOnv7w9/f38AgEKhm0GQ8m5qyYzPxDnDc5CaSKFO0++Ruhhj/9OpUyf4+Pjg9OnTePjwoSg5FFmo165dW+D86tWro23btqhXrx66d++uk8Ti4uKQnZ0NOzs7jfl2dnaIiooq9nauXr2KkSNHaju9Ysvro6lMVUJqIuXxpBmrYP744w9cvHgR8fHxouVQZKFu0aJFgfOTk5Nx7NgxbN68GWFhYbrICyqVCtevX4enpycOHDggzPf09MTBgweLvZ3WrVsjMjJSFym+kcxMBpeBLpDmSBGfHi/c1JIYlMhn1YxVIHlF2s7OLl9zbHkoslD37NmzvPIo0KpVq7B7925cu3YNFy9exLRp01C7dm1s2bIFALBz504AwLhx4wAAH330EcLCwnDnzh0YGhrCx8cHQ4cOxbBhw0TJ38rdCjK5DIY5hrkz+CktjFVYgwYNwvz58zFy5Eg8ffq0XPddpjZqd3d3TJ48GWPHjtVWPhp++eUXWFtb44svvkCtWrVw+/ZtDBgwAOHh4QCAOnXqaCxvaGiIlStXwtHREUqlEnfu3MGAAQNw7NgxneT3JsYuxkhMSgQAjRtb+KYWxiqec+fOiXZGDZSwC4mdnR3NmzePQkJCSKVSUWBgoOjdWrQR2u6eV7t2bWFAcg4ODo43RZm750kkErzzzjs4fPgwnj17hq+++gq7d++Go6MjevToUZxNVDkff/wxdu/ezWNMM1bJtGjRotC7tnWlyEJdv359LF++HM+fP8fq1avxzz//oGnTpsjJycGBAwdE+wqgz2RmMjhMdMCaH9bgq6++0vuHZjLGSqZ169YYOHAgLC0ty3W/hZ5uZ2Zmkr+/P3Xp0kVjflZWFjVp0kT0rwraDG01fdgMtCHX9a5kM9BG9PfEwcGh/ZDL5WRoaKj17RZVg4q8mHjv3j30798fiYmJSE5Oxq1bt4pavMqTmclg62ELS0NLwB3cDY+xSujV55oaGRkhMzNT5/sssumjdevWGDJkCMzNzXHu3Dn8999/mD17NgDoxWhS+sbK3QpyIzkszS15bGnGKjGJRILt27dj1qxZ5bK/N15MDA4Ohp+fH2rXro1Vq1Zh0KBBkMlk2LZtG6ZPn17ouBtVTd6t4mmJaXgY8hCqNBU/sYWxSoqIcP78+XIbQkOC3DaQEmnUqBEmTZqE9957DzY2NpDL5TpIrXwpFAq4ubmVen2bgTao4V5Do6lDaiJF4rlEvrmFMfZGRdWgUo2eFxISgrlz58LR0REjRowoU3KVhbGLMWxq2sDKwQoycxlk5jJIZBK+uYWxSkwul6NLly7lsi/Rr6LqQ5S114dEIqFJkyZR//79RX8vHBwc5RNDhw4lhUJBrq6uZd5WqXt9sOIjImzbtk3sNBhj5ejkyZOIjIzU+TNZ+SnkWiAzk6H5J81hYMafe4xVJampqbhy5YrOe8FxodaCRl6NYNbaDO2ntBc7FcZYOTMzM8PYsWPRsGFDne2DTwHLSGYmg9pVjYy4DCjtlJCZyvgmF8aqECKCn58f1Gq1zp4AU+xC7e3tjV69eqFmzZqQSjVPxAcPHqz1xCoKK3crkISQ8jKFn+DCWBWUlpaGAQMG4OXLlzrbR7GaPlasWIE9e/bAxcUFL1++RHx8vEZUVTIzGay7W8PM0AxSmVR4ggvf5MJY1aLLIg0U84z6vffew6hRo0r0CKyqwMrdCtVMqsHB1gGhGaG5T2TnJ7gwVuUYGxtj9uzZuHDhAgIDA7W+/WKdUUulUr162ri+MHYxhjJTiSdRT6CSq/gmF8aqKKVSiRYtWuR7GLe2FOuMeuvWrfDx8cGSJUt0kkRF9WzjM7FTYIzpCW9vb51tu1iF2tLSEqNHj4anpydu3rypMcwfkPtQ2arIzMwMPj4++Ouvv4TnODLGmLYVq1A3bdpUaPpwdXXVeK0qD3fq7OyMCRMm4N9//+VCzVgV5+DggKVLl2LLli1QKBRa3XaxCnXPnj21utPK4v7T+/Dy90L0PX4kGWNVXWJiIogoX/dlbeAbXsrAyt0Kpq1MYRFpwb08GKvi0tPTMXnyZJ1su9BCffjwYfj4+CAlJQWHDx8uciNV8YYXmZkMLgNdkJ2VDauuVvzYLcaYzhRaqOPj44X256p8U0th8vpQk5SQlpnGfacZYxg5ciS8vLzg5eWl1e0WWqgnTpxY4P/Z/x67Ffk0EshB7k0ufFbNWJUXHR2NmzdvQi6X5+sdVxbcRl0KVu5WubcK5fz/jLxizWfVjFVpgYGBOrkzkQt1KRi7GENuJIdpDVOkpqZCrVYL8xljTBdEf5yNPkRJH8XVvn17Onr0KNWtW1f03Dk4OPQjatasSUePHqU+ffqUeN2iapDePzjAz88PT548gVKpRHBwMLp27Vrk8t27d0dwcDCUSiUeP36MqVOn6iSva9euYeDAgQgNDdXJ9hljFU9KSgquXLmC2NhYrW/7jZVeJpOJ8unk7e1NWVlZNHnyZHJ1daV169ZRSkoKOTk5Fbi8i4sLpaam0rp168jV1ZUmT55MWVlZNGzYsDJ9mnFwcHDoOt5Qg968gZiYGFq5cqVWnrRbkrhy5Qpt3bpVY15ISAgtW7aswOW/+eYbCgkJ0Zjn7+9Ply5dKutByhce/T3onQ3vkMxUnA8xDo6qGFIpaOBA0Bdf5P4rlYqf0+shM5ORw0SHEteGMjd9fPbZZ+jcuTNu376NS5cuYeLEiTA1NS3OqqUml8vRrl07nDhxQmP+iRMn0Llz5wLX6dSpU77ljx8/jrfeegsGBtq9biptLUWUWRSsPKy0ul3GWMGkUuD4cWDfPmDx4tx/jx/Pna9PPtr2EWw72Wq1NhTrLf7444/o0qULmjdvjgsXLuCrr75CZGQktm3bVmjRLCsbGxsYGBggOlpzHI3o6GjY29sXuI69vX2By8vlctjY2GgtN5mZDM8sniHxeSI/0YWxctK/P9ChA2BuDshkuf926JA7X1/IzGS4kHgBlEparQ0l+iy6f/8+Pv30Uzg6OuKzzz7D6NGjce7cOdy7dw9Tp06FRCLRSlLlxdfXFwqFAgqFokSFPK8fNalI6D/NGNOtNm0AExPNeSYmQOvWoqRTICt3KyQmJSL1ZapWa0OJCrVcLse7776LY8eOYfXq1bhy5QrGjx+PgIAALFiwAHv37tVKUgAQFxeH7OzsfE9MsLOzQ1RUVIHrREVFFbi8SqVCXFz+G1H8/f3h5uYGNze3Al8vSN5diSYGJnB0dOTnJDJWTm7cANLTNeelpwP68vCpvNqQk5F7J5w2a0OxCnWbNm2wfv16REZGYs2aNfj333/RtGlT9OjRA3v27MGKFSvQp08frQ7OpFKpcP36dXh6emrM9/T0xKVLlwpc5/LlywUuHxwcjOzsbK3klXc2LZPKIDeUQ0ISPqtmrBwcOwZcvQqkpABqde6/V6/mztcHebXBxtoGxsbGGncsl1WxrrApFAqcOHECU6ZMweHDh4U78V4VFhaGn3/+ucwJvWrVqlXYvXs3rl27hosXL2LatGmoXbs2tmzZAgDYuXMnAGDcuHEAgC1btuCDDz7A6tWr8cMPP6BLly4YP348Ro0apbWcjF2MIZFJkKpORWpsKqRmUmE+Y0x3cnKAvn1z26Rbt849kz52LHe+PjB2MYbUQIpqVtWQJctClkGWMF8b3thtpE6dOqJ1dfHz86PQ0FDKyMig4OBg6tatm/BaYGAgBQYGaizfvXt3un79OmVkZNCTJ09o6tSpZe4aw8HBwVHckEgkZGBgUOL1ytyPuipESQu1i4sLrV27ttz7lnNwcFTOqNC3kOsrtVoNI0sjyN+R84VExhgAoFGjRpg2bRosLS21ul0u1KX07NkzfPHLF1DVUfGFRMYYAKBx48aYMGGC1rsqc6EupbyuONlJ2dw9jzEGADhy5Ag6d+6MxMRErW6XC3UpWblbobpVddjb2HP3PMaYoKBecWVVpkLdrl07dOvWTVu5VBh5Z9OybBlyKIdvemGMAQBmzZqFHj16aH27ZSrUu3fv1sljZ/RdXsf2uNg4REdFa7VjO2OsYpLJZOjatSsaNGig9W2XaUi59957Dyav33xfBeTd9CIzl+WbzxirmtRqNYYNG6aTbZepUAcHB2srjwrl2cZnAABDQ0Ps2rULf/zxh1bHOWGMsVeVqFD36NEDTZs2BQDcu3cPZ86c0UlSFUVWVhbu3buXb2hVxljVM336dADApk2btL7tYhVqFxcXHDx4EC1btkRERAQAoHbt2rh16xaGDx9epZ8buGTJErFTYIzpAUtLS50O9fzGWxtPnz5NQUFBGs8qdHJyojNnztDp06dFv/VSG1GWsT6MLI3IydeJH8vFwcFR6ijzWB/p6enUsmXLfPNbtWpF6enpor/BcjhIhYalpSV9dvgzarW1FdkMtBH9fXBwcFTMKPNYH+Hh4bnjq76mWrVqePbsWXE2UWmlZKfgZMRJqJPU3JeasSpq+fLlWLx4sc62X6xC/cknn2DdunXo0KEDpFIpJBIJOnTogDVr1uCTTz7RWXIVgZW7FV4mv4QyVcl9qRmroh4/fqzTa3US5J5aFyk5ORlGRkaQyWTI+f9RuqVSKdRqNTIyMjSWrV69uk4S1TWFQgE3N7cSrSMzk6He5/WQk5UDA5kBpDIpspGNJ189gTpN+7eRMsYqr6JqULF6fXzwwQdaTaiyyLtDETlA3fp1kZ6ejsjESFh5WCHuaPGewcgYq9hq1aqFyMhIne6jWIV6165dOk2ionr1DsUXiS+QnZ0NiUzCdygyVkUYGRnh559/xoEDB7B+/Xqd7afQQm1lZSUM1WdlVXS7q7aH9Kso8u5QZIxVTTk5Ofjuu+8QEhKi0/0UWqhjY2NRq1YtxMbGIi4uDkT5m7IlEgmICAYGZboTvdJo0qQJOnfujG3btomdCmOsHKhUKhw5ckTn+ym0wvbs2RMJCQkAoJNh+yqjt956C0NHD8UZwzMI3xXOFxQZq8Ts7e3Rpk0bnDp1CiqVSuf7E72jtz6ENp5CbmxsTPaD7cl1vSvf/MLBUcljwoQJdPXqVbK1tdXK9oqqQUW2URdXVW2jfl2WLAvmHc2Fx3MlBiXyWTVjldSOHTtw4cIFxMbG6nxfhRbqwtqlX8Vt1Jqs3K0gkUngZOeElxkvuZseY5UYEeHhw4flsq9CKyy3S5dM3uO51Eo1slRZoCzis2rGKiGZTIatW7di7969OH36dLnss9BCfe7cuXJJoLJ49eaXyIjczu9SEymfVTNWyVhZWUGlUkGpVJbbPovVZuHl5YWsrCz88ccfGvMHDRoEuVyOgwcP6iS5iqSgx3NVM6qGzHqZImbFGNO2uLg4TJs2rdz3+8arkbdv36Y+ffrkm9+rVy+6deuW6FdftRHa6PXxajRr1owUCgUNHTpU9PfGwcGhnWjWrBmZmprqZNtaGY/a2dk533xnZ2dKS0sT/eCVw0EqVXTp0oWkUqno742Dg6PsYWBgQMeOHaPvv/9eJ9svVfe8VyUmJqJhw4Z4+vSpxvxGjRohJSWlOJuoki5evAgAMDA3gN0IO0Ttj+ILi4xVUNnZ2Zg5cyays7PLfd/FGo/68OHDWL16NRo2bCjMa9SoEVatWoXff/9dJ4kZGhpi3bp1iI2NRWpqKg4fPgwHB4ci11m0aBGISCN0ParVmzRp0gTv//A+qrerzmNVM1bB3b9/H48ePSr3/RarUM+dOxdJSUm4e/cuwsPDER4ejjt37iA5ORlz5szRSWJr1qzB8OHDMWrUKHTr1g0WFhb4888/IZUWnfL9+/dhb28vRIsWLXSSX3HFpMTgn8x/IEmT8BNgGKugli5diilTpoi2/2I1faSkpKBr167o3bs3WrduDQC4ceOGzvoQWlhYYNKkSZgwYQJOnToFABg7diyePn2K3r1748SJE4Wum52djejoaJ3kVRqSVhJExUQhJz2Hu+sxVgHJZDJkZ2cjKytLtBxKdEvhqVOnhMKpS+3atYOhoaFGQX7+/Dnu3buHzp07F1mo69WrhxcvXiAzMxNXr17FZ599ptNH5BQl7yaYnIwcQAJUN6kOaXcp3wTDWAWiVquxdOlSUXMosh3h4sWLGo/WWrZsmcYYINbW1vkuMGqDvb09srOzEReneeYZHR0Ne3v7Qte7evUqxo8fj379+sHX1xf29va4dOkSatSoUeDyvr6+UCgUUCgUsLGx0ep7ADRvgjEyNIKdrR2MTY25rZqxCmLYsGGoVauW2GkUXag7duwIQ0NDYfr999+HpaWlMC2Tyd54ge9VX375Zb6Lfa+Hu7t7yd/F//v777/x66+/4tatWzh9+jTefvttSKVSjBs3rsDl/f394ebmBjc3t3wfCtrw6k0w2YbZCI0KRXpGOj8BhrEKwNLSEjNmzIC3t7fYqZSs6UMikZRpZ2vWrMGePXuKXCY8PBwdO3aEgYEBbGxsNAqonZ0dzp8/X+z9paWl4c6dOxq9VcpTYU+Asbe313iCDmNM/7x8+RKjR49GfHy82KmUrFCXVXx8fLHe9PXr15GVlQVPT0/s27cPAODg4IAmTZrg0qVLxd6fkZERXF1dERgYWOqctc3Q0BABAQH477//MG/ePLHTYYwVoGXLlrh586bo3XtfVejdMNnZ2WRj878B8JOTk8nFxUWYrlmzJmVnZ+vkLp1NmzbRs2fPqFevXtS6dWs6c+YM3bhxQ+NOv3v37tH7778vTK9cuZK6d+9OLi4u1L59ezpy5AglJSVRnTp1ynRXkLajR48e5ODgQABIZiYjh4kOJDOVldv+OTg4Co9evXqRQqGgDh06lOt+S31nokQiwZ49e5CZmTuwULVq1eDv74/09HQAuWesuvLxxx8jOzsb+/fvh7GxMU6fPo333nsPOTk5wjKurq4aFwEdHR2xb98+2NjYIDY2FleuXEHHjh0RHh6uszxL49UzfPs+9jBvY47M6EzutseYiGRmMth72+P8gfNYtmwZrl27JnZKAglyK3aBtm/fXqyNTJw4UVv5iEahUMDNza1c9+nj64MX3V/gxfMXIBCefPWEu+0xJhInLyeYeZgh9lisKCdNRdWgIs+oK0MB1mfXMq7BIs0C6kw1pMZ8MwxjYpGZyeA8wBkyiQzZXbP17l6HYt1CzrRPZiYDNSUkRicCBOGJMHyLOWPlz8rdCkkpSUiKTwKk0Lt7HbhQi+TVm2GkMinq1qkLC0sLvfsFYayyq+lcE1ZdrZD+Mh0pKSnIycjRu5MmLtQiefVmGImJBImZichUZfLNMIyVox49emDKhikwqmYE5PVTyIHenVXz48NFUtjNMEDuoFTJycnlmA1jVVNwcDDwCMhSZ2k8Rg+A3p00id5vUR+iPPtRFxWdOnWis2fPUosWLUTPhYOjskbdunVJIpGInserUVQN4qYPPXP79m0cP34cT548AZB70dFhooNetZcxVpHZ2tpix44d8PPzEzuVYuOmDz2TkpKCZcuWAcgd9Kr+0PqQtpHyDTGMaUlsbCzWrFlTonGDxMZn1HpsykdT4NjXETmp+ncVmrGKpnnz5nBycgIAHDp0SCcjZuoKF2o9dvTZUSS8TIA6Q613V6EZq0hkMhm+/vprzJ8/X+xUSoWbPvSUzEwGo7ZGSInPfcq7ARnAoY+D3t0xxVhFoFar8cknn1TYoYX5jFpPvXpDDADY1rBFdavqsO1tK2pejFUk77zzjjDw/6NHj/RibOnS4EKtp169IUZmLkN0SjRiE2Ihd5QDKPtDHBirCrp164YuXbpUir8X0fsP6kPoSz/q4sTYsWPp22+/JblcLnouHBz6FnK5nExNTQkAGRoakkxWMcZ6537UlUx2drYQAPe1ZuxVq1evxnfffQeJRIKsrCyo1RX/mg5fTKyA8h5PBuQ+gNNhqAPUbdTc15oxAEePHkV2djaISOxUtIbPqCu46Z9Mh10vO+SkcF9rVnUNHToUnTp1AgAcO3YMJ0+eFDkj7eJCXcHtvLETsfGxUGfm9rW27mktdkqMlSu5XI4RI0ZgwIABYqeiM9z0UYHJzGQwcTNB+svcZ1hWk1VDba/aMLhrgIjHESJnx5huNWjQAE+fPoVKpcL06dORlJQkdko6w2fUFdjrfa0pm5CVnYXsZtmi5sWYrtWqVQu7du3ChAkTAAAvX76sVG3Sr+NCXYG93tc6yyAL8YnxMHIygoGBAebPn4/atWuLnSZjWiOX595HEBkZieXLl+Pnn38WOaPyI3r/QX2IitSPujjRsGFDCgwMJHd3d9Fz4eDQRnTq1In+/vtvcnBwED0XXQT3o66CHj58iEGDBuHs2bMAgK5du8LJyYn7XLMK69GjR7h161al6BddUlyoK7GUlNwBnWQyGebMmYPZs2fDyt0K5m3MeSQ+ViGMGDECn3/+OYDccaTnzJmDqKgokbMqf9zrowpQq9WYMGECjCyNYOVrhZyUHNj3sueR+JjekZnJYO9tj6j9UVCnqVG9enXY2NhALpdDpVKJnZ5ouFBXEQkJCbDpZANIgRqmNWBsZYzUt1MRtj9M7NQYE9j3sYddFzvUkNbArR9vYdu2bZW6N0dxcdNHFSEzk8GqqxVyMnIQFxeHuMg4yFvLITOVoUuXLqhWrZrYKbIqTmYmg2l7U1gaWMKglQFkpjIu0v+PC3UV8WqfayJCeko6IAUcBzri+++/x6RJk8ROkVVRvXv3xqpVq1DDowZIQnga+hSp6al8HeUVeluofX19cebMGSQmJoKI4OzsXKz1hg0bhjt37iAjIwN37tzBkCFDdJtoBfF6n2uZuQwSmQQ5tjnw9fUVBnqqW7cu+vbtWynG72X6Le93TC6Xw7C6IWp0q4GcjByAgJwMHrvmVXrbRm1iYoITJ07g8OHDWLNmTbHW6dixI/bv349Fixbht99+w7Bhw/Drr7+iS5cuuHbtmm4T1nPPNj4r/DX877UhQ4Zg8ODBuHTp0v96jbx2gYexsrCyssKGDRvw008/4a+//sLff/8NhVQBq+5Wwl22yIHwnFAeEVKPC/XatWsBAO3atSv2Oh9//DECAwOxbNkyAMCyZcvQo0cPfPzxxxg9erRO8qxs1qxZg0OHDglF+tNPP8WZtDNIdU3lYVRZqUkkEtjb2yMyMhKJiYl49uwZUlNTAeQ2xb36je9Vxi7GYqSrd/S2UJdGp06dsH79eo15x48fxwcffCBSRhUPESEsLAwAYGpqiqbtmuJ03GlkJ2bDqqsVUi6mIPNlprhJsgpn/vz56Nq1KwYPHgyVSoV58+ZpvF7UNz6mx23UpWFvb4/o6GiNedHR0bC3ty9weV9fXygUCigUCtjY2JRHihVKWloa5uyeg5S0FJCKUM2kGmb4z0DdunXFTo3pOXNzc0yYMAFmZmYAgD///BPr1q0TnkrESqZcC/WXX34JIioy3N3dyy0ff39/uLm5wc3NDXFx/JX+dTIzGSy7WEKtzG2XzsnIwcWXF/Ei/gUAoG3btmjSpImYKTI9I5PlNl3UqlUL06dPFwbzv3nzJv7++2/ubldK5dr0sWbNGuzZs6fIZcLDw0u9/aioKNjZ2WnMs7Ozq5K3nGrD68OoZqRnIAtZsOhmgbijcXj//fdhbGwstP9LJBL+Q6yiJBIJfvjhB9y7dw+rV69GSEgIBg0ahMjISLFTqxTKtVDHx8cjPj5eZ9u/fPkyPD098d133wnzPD09cenSJZ3tszJ70wWeGTNmCB+MMpkMv/76K3bv3o1Dhw4Jy3KPkcqrfv36aNasGf744w8QEW7evKlxosVFWrtEH96voLCzs6NWrVrRqFGjiIiof//+1KpVK7KyshKWOXXqFC1btkxjGESVSkVz586lxo0b07x58ygrK4vat29fpiEGOd4cFhYWtHDhQurQoQMBIEtLSxo7diw5eTmR63pXshloI3qOHGUPW1tb4f8zZsygc+fOkbGxseh5VYZ4Qw0SP8GCYtGiRVSQcePGCcuEhoZSQECAxnrDhw+ne/fuUWZmJt29e5eGDh2qjYPEUcLo168fnbl8hpp+15QafNWAXFe4km0dW9Hz4nhzyMxk5DDRgWSmMo35Hh4epFAoqHHjxgSArKysyNzcXPR8K0tUyEKtZweJoxRRd2RdavhtQ6q/qD65/eBGn//xORkZGREAMjQ0FD0/joLDZqANua53pToj6pC/vz/16tWLAJC5uTlNnDiRatSoIXqOlTGKqkGVqh810x8yMxkMWhnk3hIMICUhBX9m/Ilsg2wgE1i5ciXUajVmzZolcqbsVWOnjMWtFreQmpQKo7ZGSFGlCBeIU1JSsH37dpEzrJoqVT9qpj9e7zGiylAhPSNdGGjn3LlzuHDhgrB8QEAAxowZU+C2+Kk0utOrVy+MHDlSmDbuYAzDaoYgFQFSYNnRZThz5oyIGTKACzXTkcIGgcrrMXLw4EH89ttvAAAjIyM8ffoUCQkJuesaG+PQoUPo2bMnAPBTabTIzc0Nfn5+wnTXrl0xaNAgALkfiIExgYiPyO2ZxQMj6Q9u+mA6UZJbgjMzM7F48WJh2tzcHCEhIUhMTITMTIaaPWqillUtGPQyQGJQIqRZUkgkEmRlZekg88qlVatWGD58OJYuXYrs7Gw0a9YMw4cPR0BAADIyMrBixQoolUoA+b8F8cBI+oPPqJneiYmJwdy5c3Hjxg2heGSlZUFNalh5WMHd3R3nz58XbmV3dHRE+/btIZfL37jtitqMUlTeVlZWMDIyAgC0adMG+/fvh5OTEwDA2toa7dq1Q61atQAAP/30Ezw9PZGRkQEAQpEG3vwtiImHz6iZ3sp7Kk1mSiYikiJyz+66WuHJnifYsWMHnj9/DiD3pqbp06fD3d0dKpUKvXv3hpubG7799lvk5OTAxMQEmZmZUKvVQjOKNkcCLI+bevLypkRCF4MuCA4ORmhoKJo3b46AgADMnDkTFy5cwMuXL/HixQvhQ+vMmTMabcxFPXeQB0bSX3xGzfRWYV/Fk12SsXnzZqHoHDhwAJMnT0Z6ejqA3DPstm3bIicnd8UZM2bg6NGjQuE3zjGG80Bn4ey0Vq1asLW1LVOeZWlDl8vlMDbOPWuVSqUYPnw4mjdvDiC3GWj/H/vh1M8J2UnZMO9ojumzpsPNzQ0A8PjxY6xatQpPnjwBAISGhmLWrFnCNKscuFAzvVXcr+IpKSn477//hOkdO3ZgxIgRwnRQUBC2b98uFH5TI1NUM6kmFNa5c+di1apVwvJLlizB0qVLhelJkybBx8dHmB44cCB69OgBIPds2rGvI+QquXDhbfDgwejatauw/OTJk9G7d29heu3atRrbCwoKwuTJkwEARITZs2fDw8MDAJCamoqDDw8iOycbpCLkIAfjV4/HL7/8AiC36WLfvn2IiIgo4dFlFQk3fTC9pa2v4levXkXwnWDU+7wecjJyEJkeKTSjJAYlYseOHcIZLZA7MJhU+r9zmMaNG2u05Y4ZMwaRkZEIDAzMPZuubg5JlgSRqkhYeVhh3LhxuHv3rtD9sG/fvvjnn39w6tQpAEB6errGhdB169YhJCQEQG6hHjBgAJKSkgAAUlMpbmTdQE5q7reDnIwcGLYxhOyojMdOqUIkyL3zpcpTKBTC10lW+dgMtIFVdyvkpOcI86QmUiSeSyxxW3Ve+2+OUQ7qfV4PUANQAznIgdRQiojvI5D5MlMrYy9rM2+m34qqQXxGzaoEbT7qKa9t3KaPDSCFcPclAECae9NI2tG0MuX7an78iCrGhZpVCbro0VAeRZR7YjCACzVjpcZFlJUX7vXBGGN6jgs1Y4zpOS7UjDGm57hQM8aYnuNCzRhjeo4LNWOM6Tku1Iwxpue4UDPGmJ7jQs0YY3qOCzVjjOk5LtSMMabnuFAzxpie40LNGGN6jgs1Y4zpOS7UjDGm57hQM8aYnuNCzRhjeo4LNWOM6Tl+Cvn/S05OxoMHD8ROo0KysbFBXBw/Ebs0+NiVTmU8bs7OzqhZs2ahrxMHSKFQiJ5DRQ0+dnzs+LjpNrjpgzHG9BwXasYY03NcqP/f1q1bxU6hwuJjV3p87Eqnqh03vpjIGGN6js+oGWNMz3GhZowxPceFmjHG9BwXasYY03NcqBljTM9xoWaMMT3HhZoxxvQcF2rGGNNzXKgZY0zPcaFmjDE9x4WaMcb0HBdqxhjTc1yoGWNMz3GhZkxLAgMDERgYKGoOoaGhCAgIEDUHpn1cqCsoIipWjBs3TuxUBRMmTAAR4f79+2KnUiEEBgYKP0e1Wo2kpCTcv38fu3btQu/evbW2n06dOmHRokWoXr261rbJtMtA7ARY6fj4+GhMT5kyBR07dsTEiRM15l+6dKk80yqSj48PQkND0bhxY7z11lsIDg4WOyW9FxERgU8//RQAYGpqigYNGmDYsGEYO3Ys9u/fDx8fH2RnZwvLN27cGDk5OSXaR+fOnbF48WLs2LEDSUlJWs2faY/oD27kKHsEBASQUql843ImJiai5Ofg4EDZ2dnk7e1N4eHhtGbNGtGPmbYjMDCQAgMDtbq9e/fu5ZsvlUppw4YNRET0zTfflHk/n3zyCREROTs7i34MOQoObvqoxAICAqBUKuHs7IzDhw8jKSkJR48eBVB4e2pAQABCQ0Pzzf/ggw9w8+ZNKJVKREdH48cff4S1tXWxcxk9ejTS09Pxxx9/YP/+/Xj33Xchleb/9SMibN68GYMHD8atW7eQkZGB27dvo2/fvvmWbdmyJY4ePYqkpCSkpqYiMDAQXbt21Vhm3LhxICK4u7tj7dq1iImJQWJiIvz9/WFoaAgLCwsEBAQgISEBCQkJ+O677yCRSDS2MWvWLJw/fx6xsbFQKpW4efMmJk2aVOT7lUgkePr0KQ4fPpzvNZlMhqioKPz888/FOXT55OTk4MMPP8SdO3fwwQcfwMLCQnitoDZqPz8/3Lx5E6mpqXj58iVu3LiBKVOmAAAWLVqE7777DgAQFhYmNLW4u7sDAN555x388ccfePbsGTIyMhAWFoYVK1bAyMhIYx95v2u1a9fGoUOHkJKSgpiYGKxcubLAn/P06dNx48YNpKenIzY2FidOnMj3sxs1ahSuXbuG9PR0JCQk4JdffoGLi4vGMvXr18f+/fsRERGBjIwMvHjxAgcOHIC9vX2pjq2+4kJdyUmlUpw4cQLJycmYM2cOdu/eXeJtbN68GatWrcLVq1fx0UcfYevWrfDy8kJgYGC+P9jC+Pj44PDhw8jIyMC+fftgb28PT0/PApft1KkTNm7ciJ9//hmffvopqlWrhoMHD6JGjRrCMq6urjh//jzatGmDlStXYvHixahduzZOnTqFbt265dvmmjVr0LhxYyxevBh//PEHJk+ejK+//hqHDx+GkZERPvvsM1y4cAGffPIJxo8fr7HuzJkzcfv2bSxduhRz5swRPqimTp1a6PslIuzZswf9+vXTyBsA+vTpAzs7O+zatatYx64gOTk52LdvH0xNTfMVuFdNnDgRmzZtwv379zFr1ix8/vnnuHz5Mjp37gwA+O2337B3714AwMcffwwfHx/4+Pjg3r17AHKvK2RmZmLdunX48MMPcebMGcycORM7duzIty+pVIq///4b8fHxmD17Ns6ePYvZs2cLHwp5tm7dio0bNyIqKgrz58/H119/jZcvX6J79+7CMnPnzsWePXsQGhqKTz75BN999x26du2KixcvwsbGBgBgYGCA48ePo2vXrti4cSOmT5+OTZs2wc7ODrVr1y71sdVXop/Wc5Q9Cmr6CAgIICKi77//Pt/yhX1NDwgIoNDQUGG6U6dOREQ0duxYjeW6dOlCRES+vr5vzK1FixZERDRw4EBh3oMHD2j37t35liUiyszMpPr16+db//333xfmHTx4kDIzM6lBgwbCPGtra4qNjSWFQiHMGzduHBERnThxQmM/Fy9eJLVaTVu3bhXmSaVSCg8Pp/Pnz2ssa2xsnC/P48eP08OHD4s8po0bNyYiIj8/P43lfvrpJ4qOjiaZTFbkcSus6SMvBg8eTEREM2bMEOaFhoZSQECAMP3bb7/RrVu3itxPUU0fBb33+fPnk1qtJkdHx3y/awsWLNBY9vr16xo/D3d3dyIi2rBhQ6H5ODk5UVZWVr5t1atXj5RKJX399dcEgFq2bElERMOHD9fJ35Q+BZ9RVwGbNm0q9bre3t5ISUnB33//DWtrayHu37+PqKgo9OjR443b8PHxQXx8PI4fPy7M27dvH4YMGQITE5N8ywcGBuLx48fC9K1bt5CUlIR69eoByD1z69u3L44cOYJHjx4Jy8XHx2PHjh146623ULNmTY1tbt++XWP66tWrkEql2LZtmzAvJycHwcHBwn7yKJVKALlncFZWVrC2tkZgYCAaNGig0ezwugcPHuDq1asYO3asMM/U1BSDBw/Gvn37oFarC123OFJTUwEA5ubmhS6TlJQER0dHvPXWW6XaR957l0gksLCwgLW1NS5cuACpVIq2bdvmW97f319j+vz58xrH08vLC0Buk0thhg0bBrlcjv3792v8ziUlJeHWrVvC71xycjIAoG/fvgX+HlUmXKgrObVajbCwsFKv36hRI5ibmyMmJgZxcXEaYW9vn68gvk4ikWDUqFE4e/YsnJ2dUb9+fdSvXx/Xrl2DmZkZhgwZkm+d8PDwfPMSExNhZWUFALC1tYWpqSkePHiQb7m8r+yvt2W+vs283g3Pnj3LNz9vP3kGDRoEhUIBpVKJhIQExMXFYfny5QDwxi5tO3fuRKdOnYRiNWzYMJiampaqCep1ZmZmAICUlJRCl/n222+RkpIChUKBR48eYfPmzfDw8Cj2Ppo1a4ajR48iNTUVSUlJiIuLw7lz5wDkf+9ZWVmIiorSmJeYmKjR9FO/fn1ERUUhPj6+0H02atQIQO4H3eu/c25ubsLvXFhYGL7//nv4+voiLi4OJ0+exIcffpivqaky4O55lZxKpSrwzI2I8l00A3IvdL1KKpUiLi4OI0eOLHD7iYmJRe7fw8MDTk5OcHJywrBhw/K97uPjI7SR5insTLOgfIursG0WNP/V/XTp0gWHDh3ChQsXMG3aNERERCArKwsDBgzArFmzCrxQ9qqff/4Zq1evho+PD5YuXSq0/16/fr3U7yVP8+bNAUDjW8Xr7t+/j8aNG2PAgAHo27cv3n77bUybNg0bN27EBx98UOT2LSwsEBgYiLS0NHz++ed49OgRlEolHBwcsHPnznzvvaTdAguTt91+/fppdD3Mk3eWDwCzZ8/G9u3bMWjQIPTp0wfff/89vvjiC7i7uwsf2pUBF+oqKjExMd9XfABwdnbWmH78+DE8PT1x5coVpKWllXg/Pj4+iI2NhZ+fX77X+vbti/Hjx8PW1haxsbHF3mZsbCzS0tLQuHHjfK+5uroCQJm+RbzKy8sLGRkZ6NOnDzIzM4X5xWnyAXKP859//gkfHx9s3boVPXv2xMKFC8ucl1QqxejRo5GWloYLFy4UuaxSqcTBgwdx8OBByGQy7NixA++//z6WLVuGiIgIEFGB6/Xo0QO2trbw8vISzqIBlOlmm8ePH6N///6wsbFBXFxcocsAud+CilNs7969i7t37+Kbb75BixYtcP36dcycOTPfRcyKjJs+qqjHjx/D1dVVuIIO5HZ369Kli8Zy+/fvh0wmK7C4SKVSWFpaFroPIyMjDB8+HH/99ZdQKF6N77//HnK5vNCz9cLk5OTg77//xjvvvKPxYWNlZYVx48ZBoVAgJiamRNssjFqtBhFpnD1aWlrmu7GoKDt37kTDhg2xZs0aSKVS7Nmzp0w5SaVSrFu3Dk2bNsW6deuKbPp4vRlArVbj1q1bACD87PI+gF9v8sn7tvHqNwyJRIJZs2aVOvcDBw4AABYvXlzoMgcPHkR2dnahH2h53ULNzc3zfQO8d+8elEplkb+XFRGfUVdR27dvx6xZs3D8+HFs27YNNWvWxLRp03Dnzh2NC2Tnz5/Hhg0b8Omnn6Jly5Y4fvw4MjMz0aBBA3h5eWHhwoXYuXNngfsYNGgQqlevjj/++KPA1x88eICQkBD4+Phg/fr1Jcr/iy++QJ8+fXDhwgVs3LgRGRkZ8PX1haWlpXDBShuOHDmCTz75BCdPnsTu3btRo0YN+Pr6IioqCrVq1SrWNo4dO4aYmBi8++67CAwMzNcuXhQLCwuMGTMGAGBiYiLcmdigQQPs27cPCxYsKHL9EydOICYmBhcuXEBUVBQaNGiAGTNm4L///hPOVvPuEF2+fDn27t2LrKwsnDlzBhcvXkRcXBx27tyJ9evXQ6VSwcvLS2gbL42zZ88iICAA77//PurXr49jx44ByO2SefPmTSxfvhyhoaGYN28evvvuOzg7O+P333/Hy5cvUbduXQwePBj79+/HkiVL0LNnT2zcuBEHDhzAgwcPIJFI8O6778LCwgL79+8vdY76SvSuJxxlj8K65xV1t+Lo0aPp0aNHlJGRQf/88w95enrm656XFxMmTKCrV69SWloaJSUl0a1bt2jlypXk5ORU6PYPHz5MGRkZZGZmVugyK1asICKihg0bEpDbPW/z5s35lnu92xmQ2z3r6NGjlJycTKmpqRQUFETdunXTWCave16HDh005i9atIiIiOzs7N54zMaOHUt3794lpVJJDx8+pJkzZ9L48ePzdWkr6s7ENWvWEBHRhAkTiv0zDQwMpFclJyfTgwcPaNeuXdS7d+8C13n9OPn6+lJgYCDFxMRQRkYGPXnyhNauXUu2trYa682dO5eePn1K2dnZRETk7u5OAKh9+/Z0/vx5Sk1NpaioKNq4cSM1b96ciIjGjRv3xt+1vOP86jyJREIff/wx3b59mzIyMig2NpaOHz9OnTt31lhu0KBBFBQUJPx879+/T5s2baKmTZsSAHJxcSF/f38KCQmhtLQ0io+Pp3PnztGgQYNE/3vUQYieAAdHpY8VK1ZQeno6mZubi54LR8ULbqNmTMcMDQ0xduxY4dZqxkqK26gZ0xFbW1v07t0bQ4cORc2aNbF69WqxU2IVmOin9RwclTHybpeOjo6mDz/8UPR8OCpuSP7/P4wxxvQUt1Ezxpie4zbq/xcTE4OnT5+KnQZjrIpydnYudOwcLtT/7+nTp3BzcxM7DcZYFaVQKAp9jZs+GGNMz3GhZowxPceFmjHG9BwXasYY03NcqBljTM9xoWaMMT3HhZoxxvQcF2rGGNNzXKgZY0zPcaFmjDE9x4WaMcb0HBdqxhjTc1yoGWNMz3GhZowxPceFmjHG9BwXasYY03NcqBljTM9xoWaMMT3HhZpVOXM6G8LDRaYxz8NFhjmdDavE/lnFw4WaVTmKCDV+8TIWiqWHiwy/eBlDEaEucHmJRCL8v1atWrC1tRWmmzRpgrp16wrTjRs3Ru3atYVpBwcHWFhYlGn/jAEAcYAUCoXoOXCUT8hkMvJwkVHMbDPa5lOfEufXIA8XGQGgefPm0cyZM4Vld+7cSStXrhSmf/nlF1q+fLkwfejQIVq6dKkwffToUVqwYIEwffz4cZo/f74w/ddff9GMGTOE/V9Y2p9efva//ffv35/q1q0r+jHiKP94Qw0SP0F9CC7UlTO6dOlCI0aMEKa/+eYb2rlzJwGgJR5GRIss6N/V/3t99uzZ9OGHHwrTXl5e1K9fP2G6a9eu1KpVK2G6RYsW5OLiIky3bduW6tevL0z37NmTmjVrJkz7+flR9+7dCQAt7ZG7/wPvtyEAJJfLSaFQ0IQJEwgAGRoa0okTJ2jw4MEE5H7ADB06lGrXri36ceXQfnChLvtB4tDTMDU1pebNmwvTEydOpH379gnTn3/+OR09elSY7t+/P40YMUI4o13RvzrFzjEXzmjLK/L2v8TDiGJmm5GHi4wkEgk5ODiQlZUVASBzc3OaP38+ubm5EQCqU6cOKRQK6t+/PwEge3t72r59O7Vp879Cb2pqqrGfOZ0N8703DxcZzelsKPrPjkMzuFCX/SBx6Ek4OTnRmDFjSCqVEgDy9fWlq1evkrGxMQGgPn360Pz580kikRAAqlatmvD/vMgrknkF7PVpXUdp9y+RSMjOzk4oxnXr1qXNmzdTkyZNCAC5ubmRQqEQCretrS1N69uMYl75ICrv98pR/OBCXfaDxCFSODk50SeffELW1tYEgPr160cKhUJox3VycqKOHTuSXC4v9jbFPsvU1f5r165NEyZMIHNzcwJAw4YNI4VCQcPb1qSY2Wa0YXhtSvysBvVuYCT6z5Ujf3ChLvtB4iinsLGxofnz5wvNGa6urnThwgXh67+xsTFZWlqKnmdFCBMTE+G45bXHP90xhWSy3A+JDh06UI8ePUTPkyM3iqpB3D2P6dSb+gzL5XJMnjwZXbt2BQAolUr07t0bjo6OAIAHDx7Aw8MDCoVCeP3ly5fl9wYqsPT0dCgUCni4yOD3lhxfnc+C6b396OaU+/rIkSMxdepUYXkPDw80a9ZMpGzZm4j+SaIPwWfUuomC2mMT5lvRx0PaCcscOXKEZsyYIUy/3qbMod3jnzctk8moVq1awrKHDx+mb7/9Vphu0aIFGRryRcfyCm76KPtB4ihD9GlkTPHzLIUeDr+v+UToIgegRO3LHCWLkrSHm5mZCYXb3Nycrl69Sn5+fsLrRkbctq3L4EJd9oPEUcJ49ax47ty59HTHFKJFFrTEw4js7Oy4OOt5yGQy6tixIzk6OhIAaty4MZ07d47eeust0S/GVtbgNmpWrtzc3HDs2DHhVuqnZ/fA/P4v+PJcJvzekqOJcRxUKpXIWbKiqNVqXLlyBc+fPweQ29597NgxPHjwAIoINX4bZYEt88bAzMyMb4EvJ6J/khQW3bp1o8OHD9Pz58+JiGjcuHFvXKd58+YUFBRE6enp9Pz5c43beYsKPqMufcjlchoyZIjQn9fOzo5WrFhBdevWFb3PModu4qtJfShjqQMt7VGNYmab0ZjOTmRiYiJ6XhU5KmzTR//+/enrr7+m4cOHU1pa2hsLtbm5OUVGRtL+/fupWbNmNHz4cEpOTqZZs2aV9SBxFBB5N5kYGRlRYGCgxq3XecFfkytvLO1RTWjO2r9/P23ZskX0nCpyVNhC/WqkpKS8sVBPmzaNkpKSqFq1asK8zz//nJ4/f17Wg8TxWixZsoT8/f2FaXt7e9Fz4ii/eP0W+PHu9ah169YE5H7D+v7776lly5ai51mRosq0UXfq1Annz59HRkaGMO/48eNwcHCAi4uLeIlVAtWqVcPbb78tDPl55coVnDlzRng9KipKrNRYOctrk/Y+oMSioEx4H1BihVsMLF/eAgA4OjqiQYMGMDIyAgAYGxvDxMREzJQrvEpVqO3t7REdHa0xL2/a3t4+3/K+vr5QKBRQKBSwsbEplxwrqm7dumHRokVo27YtAODYsWPYt2+fyFkxMbjVlsH7gBJBYbkXD4PC1PA+oIRb7dwbm0JDQzF06FDhJqWRI0fizz//hKWlpVgpV3gGYicgJn9/f/j7+wOA8EvFchkYGGDSpEkIDQ3FiRMncPr0aUycOBG3bt0SOzUmspWXsvLNCwpTC4UbAHJycoT/X7p0CUQk3FHauXNn9Da+g3OPUzXW8XCRwa22rMDtV3WV6ow6KioKdnZ2GvPypvmrefHkNW1kZ2ejc+fOwi3FOTk5XKRZqTx48AA7duwAkNsMsnz5crh0HspPuSmBSnVGffnyZXz77bcwMjJCZmYmAMDT0xMvXrxAWFiYuMlVAH369MGkSZPw3nvvITMzE5MnT+b+zkyrlEol3nvvPWRlZWGXkRK/epviYqYrOhveg/ev6Rpn2Ox/9PqM2tTUFK1atUKrVq0glUpRp04dtGrVCk5OuaPKLFu2DKdOnRKW37t3L9LT07Fjxw40a9YMQ4cOxbx587Bq1Sqx3oLek0gkkMvlAHLb8yMjI2Fubg4AXKSZTjx9+hSRkZEIClPjbFp9DLa8j9135Vyk30D0bimFhbu7OxUkICCAAFBAQACFhoZqrNO8eXM6e/YsKZVKioiIoIULF5a5a0xlDWNjY/rpp59o8uTJoufCUfUir4vf9wNrCDdBTZ48mdq1ayd6bmJEUTVIr5s+zp49q/EE6NdNmDAh37zbt2/D3d1dl2lVeHK5HCqVCkqlEteuXcPjx4/FTolVMa928QsKS8WROzL8MsIEIS3rwcjICNevXxc7Rb0j+ieJPkRVOaPu1q0b/f3332Rrayt6LhxVNwq7Y3V+N2NhlD5XV1eaMGFClRnAq8rc8MLe7NGjR/jvv/8glfKPnoln5aWsfG3SQWFqLD+vFDoCuLu7Y+TIkahWrZoYKeod0T9J9CEqyxl1QWcqqz8aRn8v6C96bhwcJY3q1asL/981tR31aWSs8XplGjeGz6irEEWEOl//1EkW5/FMbVNkez9j+igpKQkA0Lx5c7QfMhW/eJlU2b7Xon+S6ENUljNqANSrviGlLLCljV61eUhRjkoTrq6u1KOuAcXMNqM1g2tSzBzzSvW7zWfUVYwithoiHQdierNUbA5Wcf9UVincv38fgaHZ2HojBx+1zkCIRdcq87vNhboS6dSpEwCgbY00WD06iKVnc5+o8vpTwBmrqDxcZJjSRoodoTXRJPUyPFxkVaZJT/RTfn2Iit700bVrV1IoFDRrmBs/UYWjUkZhTwv67oPB9PXXX5NMVrF/x7npoxLLO5u4cOEC5s2bB4Pom0UOQclYRVXY8KoOiIZSqYRaXbmbQUT/JNGHqIhn1G5ubrRv3z6ysrISPRcODn0IW1tbcnd3Fz2P0gSfUVdSSUlJUCqVwpM0GKvqJk6ciKVLl2JBT4t812Y8XGSY09lQpMzKTvRPEn2IinJGLZPJqH379qLnwcGhjyGXy6lJkyaFtmfr87UaPqOuRMaOHYsNGzagbt26YqfCmN5RqVS4d+8egsLU+OaxK/6aUBPL+5i/MgBUxWzH1uvR81h+e/fuRVhYGEJDQ8VOhTG9divVCneMHTGv0xksPZt/bJGKhM+oKwAHBwcsWrQIhoaGyMrKQlBQkNgpMab3VA/PwDnmVO79BG6GGN62ptgplRoX6gqgSZMm6Nq1KxwcHMROhbEK4dXxrhcFZeKw4VDsHSpH/yZmYqdWKtz0occMDAyQnZ2NU6dO4fLly0hLSxM7JcYqhNf7XC/dFYiwbvXQ3OoqjomcW2mJfrVTH0Lfen00bNiQ/vjjD2rcuLHouXBwVKaoV68eubq6ip7H68G9PiqglJQUPH/+XBjmkTGmHYsWLcK+DzqiR13NBgV972ct+ieJPoS+nFFbWFiIngMHR2WOWrVq0Qg3e73rZ11hH25b1VhZWWHXrl04dOgQtm/fLnY6jFVKkZGR+DUSiI2V4fcxltj2L2Fskyy97mfNTR965OXLlzhx4gQuXLggdiqMVXq3Uizxsv4wzHorW+/HbedCrQdkMhlMTU1BRFi/fj1CQkLETomxSq+F+UuYPfgFX57T/3HbuVDrgTlz5iAgIICftsxYOcnrZ+21PxULAzMx/i8J/vCx0ttizYVaD5w4cQJHjx5FRkaG2KkwViW83s/auEkfhHf4Ev1a6O/di6JfhdWHEKPXh6Fh5XjMPQdHZQhnZ2dR98/9qPWQra0tDh06BE9PT7FTYYwBePr0KeZ0NsS0vs00hmvQh/7VXKhFkpGRgZs3b+Lhw4dip8IY+3//xcmwtmsSvp02EMD/2rIVEeL3CCn0dDs5OZmSkpLyRXh4OJ08eZI8PT1F/7qirdCXG144ODjEjXHd61HMHHNa4mFUrjfBlPqGlw8++KDA+ZaWlmjXrh2OHDkCLy8v/Pnnn0Vtpkz8/PwwZ84c1KpVC3fu3MHHH39caD9jd3f3AocAdXV1xYMHD3SWY0m0adMGo0ePxtKlS5GSkiJ2Ooyx1+w89wT1pEZY6G6Etf8aISgsVeyUih49b9euXUWufOPGDXz22Wc6K9Te3t5Yu3Ytpk+fjgsXLmD69Ok4duwYmjZtimfPnhW6XtOmTZGQkCBMx8bG6iS/0nB0dISTkxOys7PFToUxVgAPFxn83pIjMMcNfm/dx/EQMxy7J36xLvWpesOGDSkhIUFnXwWuXLlCW7du1ZgXEhJCy5YtK3B5d3d3IiKytrbW6tcObYdUKi23fXFwcBQ/Xh3zo169ejR7ePtya/7QWa+PatWq6azvr1wuR7t27XDixAmN+SdOnEDnzp2LXDc4OBgRERE4deoUPDw8dJJfcczpbCh0oG/Tpg3atm0LDxcZPunIQ6wwpo9e7V/95MkTfHfwGrwPKNHeQdy/2TIV6smTJ+Pff//VUiqabGxsYGBggOjoaI350dHRsLe3L3CdyMhITJs2DcOHD8ewYcPw4MEDnD59Gl27di1weV9fXygUCigUCtjY2Gj9PSgi1PjFyxgeLjJMnjwZq2YM1ZsryIyx/FZeyv9sRZVDR7R+PwBmZuI9HabIj4m1a9cWOL969epo27Yt6tWrh+7du+sksdIICQnRGCfjypUrcHFxwZw5cwq8AOnv7w9/f38AgEKh0Ho+QWFqeB9Q4hcvY/gfmoPmbeXw0uMRuhhj+Q2xfQqzuBuwsLBAampuW7WHiwxutWVYeSmrXHIoslC3aNGiwPnJyck4duwYNm/ejLCwMF3khbi4OGRnZ8POzk5jvp2dHaKiooq9natXr2LkyJHaTq/YrsebYMv1LCzoboSlZzO4SDNWwRz9NxLjGmxGI0MlIqD5PMbyUmSh7tmzZ3nlkY9KpcL169fh6emJAwcOCPM9PT1x8ODBYm+ndevWiIyM1EWKxbJ9wQQMVP6KpYEJuVeSw7K5WDNWgQjfjEeYICilHjzMHosydnWpr1K6u7vT7t27dXYV1NvbmzIzM2nSpEnk6upKa9asoZSUFKpTpw4BoJ07d9LOnTuF5T/66CMaPHgwNWjQgJo2bUrLli0jIqKhQ4eW6YpracPDRUaJn9Wgz0d3y3dFWVfHjIODQzexf1pLokUWtNGrtk62/4YaVLKN2dnZ0bx58ygkJIRUKhUFBgbq9OD4+flRaGgoZWRkUHBwMHXr1k14LTAwUGP/c+bMoZCQEEpPT6f4+Hg6d+4c9e/fXxsHqVQxp7NhvqLs4SKjOZ15MCYOjooUuSdZ5rRheC2dnWyVuVBLJBJ655136PDhw5SVlUXZ2dm0YMECsrOzE/0AltNBKnF07tyZ3nnnHZJIJKK/Nw4OjtLH69+E84q2tot1qftR169fH8uXL8fz58+xevVq/PPPP2jatClycnJw4MCBfF3n2P/069cPPj4+YqfBGCuj18eurtl+OG43nY0OjvJyy6HIi4l3797Frl274O3tjYsXL5ZXTpXCwoULYWVlBSISOxXGWBm83gWvu+wGwsJMsek/OYDcG/503V2vyDPqe/fuoX///hg8eHChXfVYfgYGuZ9/iYmJImfCGNO2A5ef4O30X9HOOh1A+QyFWmShbt26NYYMGQJzc3OcO3cO//33H2bPng0AfKZYCDc3Nxw9ehQNGzYUOxXGmA7kddc74G2GdUPthT7Vuuyu98ZbyIODg+Hn54fatWtj1apVGDRoEGQyGbZt24bp06cXejt3VZWamorg4GCd3QjEGBPf2ac5iK83BDNapmNzsKpc+lSX+Opko0aN6Ntvv6XIyEhSqVSiX5XVRvCDAzg4OIobHi4ySphnSSsHWGqtu55W+1G/GjKZjIYMGSL6QdNGaKNQt2jRgiwsLER/LxwcHLqLgrvrlb1Y62yYU7Vajd9//70sm6g0JBIJVqxYgQULFoidCmNMh17trtewYUN0Hj0Xow5lwa22TGf75IGRtYSIMGPGDEgkErFTYYzp0Ktd8D56C+jczg6//uqM05dyH1Sti656/BRyLXr06BE/VZyxKuTglVA4XvkCDqonAHTXVY8LtRYYGhri/fffh5OTk9ipMMbKUWBoNrx/TccvXsZY4mGks656XKi1oFGjRhg7diwcHBzEToUxVs7upNXA81r9sNDdSGdd9YrdRu3t7Y1evXqhZs2akEo16/vgwYO1nlhFcvv2bfTq1Utnz49kjOmvFmYv4ZpyETvD7OH3VoROxpwvVqFesWIFPv74YwQGBiIiIoLvSixAWlqa2CkwxsqZh4sMPw83woAdMQgKi9R4+ou2i/Ub+/dFRUXR8OHDRe+/qMsobT/qunXr0rfffkvOzs6ivwcODo7yjVfHnJdIJCSXy0s95nyZ+1FLpVKdPW28orO1tYWrqysyMzPFToUxVs5WXsrtP923sQmOHz+OkSNHIihMDUWEGnM6G2ptP8Uq1Fu3buWxlQtx7do1DB48uEQP3GWMVR6KCDV2vyPFvb+24v79+zrpolesNmpLS0uMHj0anp6euHnzJlQqlcbrH330kdYSYoyxikR4+K3X33hpqoKfu/bbqItVqJs2bSo0fbi6umq8VpUvLNrZ2WHNmjVYtWoVFAqF2OkwxkQSFKbG5mAVFrobYfkltTi9Pnr27KnVnVYWxsbGiImJQXJystipMMZE5OEiw/sdqiGi4Rj4qX7GiYcZ5X9GzQoWFhbGzT6MVXF5bdLv/pqO2hEpyLyv/TsUCy3Uhw8fho+PD1JSUnD48OEiN1LVb3hhjFVdGg+/fbIbABATmzswk84LdXx8vND+HB8fr5WdVTY//PADHjx4gFWrVomdCmNMJCsvZQld8a7FGKF27doIevQIADCns6FWRtErtFBPnDixwP+z/7l9+zaeP38udhqMMZEpItT4xcsY+7L7of27n+Dz0d3x83AjeB9Qam0fot/dow/Bj+Li4OAoS3i4yCh+bnXaMa4hxcw2L/ETX3T2hBfGGGO5gsLU2HA1A+NcorE5OEurvT64UJdShw4dcPr06Xz9yhljVZOHiwx+bobYeNsM09sbwcNFe4/m4kJdSvHx8fj777/5QitjTOiiN/WkITrMO4w1L1rjFy9jrRbrN7adyGRlfxR6acPPz4+ePHlCSqWSgoODqWvXrkUu3717dwoODialUkmPHz+mqVOnlrl9iIODg6OoyBtFTyaTUffu3cnW1rbEo+i9oQa9eQMxMTG0cuVKcnV1Ldc37+3tTVlZWTR58mRydXWldevWUUpKCjk5ORW4vIuLC6WmptK6devI1dWVJk+eTFlZWTRs2LCyHqR8P5BX55V2WEMODo6ShVQKGjgQ9MUXuf9KpeLnlBdlrQ1lLtSTJ0+mixcvUnZ2Nl26dIkmTpxIpqamOn/jV65coa1bt2rMCwkJoWXLlhW4/DfffEMhISEa8/z9/enSpUtlPUgaBz5mthl9O20gnT9/nnrUNaCY2WYlvsLLwcFRspBKQSdPgpKTQdnZuf+ePKk/xTqvNkzs0ZAcHByE6eLWhjIX6rxwdXWlFStWUEREBCUnJ9O2bduoc+fOOnnTcrmcVCoVeXl5aczfsGEDBQUFFbjO2bNnacOGDRrzvLy8KCsriwwMDMpykPL9QBLmW9GxT7txkebgKKcYODC3OBP9L5KTc+eLnVteeLjIKPMrRzo+r3uJa4PWuufdv38fn376KRwdHfHZZ59h9OjROHfuHO7du4epU6dCIpGUZHNFsrGxgYGBAaKjozXmR0dHw97evsB17O3tC1xeLpfDxsYm3/K+vr5QKBRQKBQFvl6YoDA11l9ORz/j/3T2MEvGmKY2bQATE815JiZA69aipFOgoDA1dt4xQB+jf7VaG0pUqOVyOd59910cO3YMq1evxpUrVzB+/HgEBARgwYIF2Lt3r1aSKi/+/v5wc3ODm5sb4uLiir2eh4sMfm/JsfRsJvzekmv1yi5jrGA3bgDp6Zrz0tMBfXr4lIeLDEMc4nRSG954St6mTRtav349xcXFUWRkJH377bfUsGFDjWWaNm1K6enpWvsKoa9NH3ntTsum9CeFQkG96sm5+YODoxyiorRR59UCbbZRF2uYU4VCgRMnTmDKlCk4fPgw1Or8p/NhYWH4+eefi7O5YlGpVLh+/To8PT1x4MABYb6npycOHjxY4DqXL1/G0KFDNeZ5enoiODgY2dnZWskrb6SsJKsoPJX8iLPhBO8DSq2OlMUYyy8nB+jbF+jfP7e5499/gWPHcufrA7faMoz6LRO+X++CyU8/4a+//tJqbXhjpa9Tp44on1De3t6UmZlJkyZNIldXV1qzZg2lpKQI+ezcuZN27twpLJ/XPW/16tXk6upKkyZNoszMTK12z+Pg4OB4PfK65hkbG9P3339PPXr0KP9+1GKGn58fhYaGUkZGBgUHB1O3bt2E1wIDAykwMFBj+e7du9P169cpIyODnjx5otMbXuRyuag3A3FwcOhHlLXZA6jghbq8oqSFukWLFqRQKKhjx46i587BwSF+5BXnJR5GpbpuxaPn6cDz58+xZcsWvHjxQuxUGGN6IChMjb9ia2GhuxF23zfk0fP0QWJiIrZt24Znz56JnQpjTA94uMgwqFYcrhh2x9gmKh49T18YGBigRo0aYqfBGBNZ3uh5w35ORqfP/4T3r+laHT2vTIW6Xbt26Natm1YSqYjWrVuH7777Tuw0GGMi03jALXKbQfK65mlDsfpRF2b37t1o1KgRDAzKtJkKa9++fZDL5WKnwRgT0ZzOhlBEqBEUpsbvv/+Oo0eP4uHJ7XCrLdPKg22BMhbq9957Dyav33xfhZw/f17sFBhjIst7sO2o3zJx9uxZmMXewC9exlp9sG2ZCnVwcLC28qiwateuDZVKhdjYWLFTYYyJIK+Z4xcvY2wO3gS/t+QazSDaUKJC3aNHDzRt2hQAcO/ePZw5c0ZriVREhoaG+O2337Br1y5s2rRJ7HQYYyIJClNj201gobsRlp7N1MlwEm/siO3i4kLXr18nlUpFT58+padPn5JKpaJ//vmH6tatK3pHc21EaW8h79WrV6FPnOHg4Kga4eEio/TF9hS8clipB2kr8w0v27ZtQ0pKCurVqwdnZ2c4OzujXr16ePnyJX788cfibKLSOn36NPelZqwKy+uat/huXSw+myU0g2h7+OM3Vvr09HRq2bJlvvmtWrXS6tCmYkZpz6gNDAyoU6dOfFbNwVFFQ1vPUS3zGXV4eDiMjY3zza9WrVqVP5s0MjLCmjVr8Pbbb4udCmOsnOV1zXsqqQNzc3MAuWfY2uyal+eNlX7gwIF09epV6tChA0mlUpJIJNShQwe6fPkyvf3226J/omkjyjLMafPmzcnQkJ9CzsFR1SJvIKZjmz6j7du3l2rUvLwoqgZJ/v8/RUpOToaRkRFkMhly/n+UbqlUCrVajYyMDI1lq1ev/qbN6SWFQgE3Nzex02CMVTAeLjIcHGmO38ItMbh2XKm75hVVg4rVPe+DDz4o8U6rEolEAm9vb8TExCAwMFDsdBhj5SgoTI0NV5RY6J6DpWd187DrYhXqXbt2aX3HlQkRYdCgQQgJCeFCzVgV88UYd7zv+C+Wns2A31tyBIZla71YF1qorayskJiYKPy/KHnLVWVTpkxBWlqa2GkwxsrRgKbmWNjkCdZFeWBR0GEEhmULt4+Xy52JsbGxqFWrFmJjYxEXFwei/E3ZEokERFRlB2V6FRdpxqqeZpaZGL4/BefCgwBojppXLoW6Z8+eSEhIAJB76zh7s0GDBmH48OEYP358gR9sjLGK79XR8nK74GXldslrZoiVl7IQFKYuv6aPc+fOFfh/VriUlBRERUXB3NwcycnJYqfDGNOBvNHyPv/XEQ08J+DS3m/h75ml1dHyXldkG3VxcRt1rsDAQL6YyFgll9e88fvoGESY3cKEPip4/6rdNunXFVqoC2uXfhW3UResRo0akEqliIuLEzsVxpgOBIWpsfZSGha6/4GliiydFmmgiELN7dKlY2RkhAMHDuDMmTP46quvxE6HMaYDw9rawu8tJZaezdRZl7zXiX4bpj5EWW4hfz369OnDgzRxcFTSeKeFJWV95UjfTB1IAMp02/irUeZBmby8vDBo0KB88/N6OTBNJ06cqPKDVTFWWTW3zMQ3j5tg4183AWj/QbaFeWOlv337NvXp0yff/F69etGtW7dE/4TTRmjzjBoA2djY0IIFC6hmzZqivzcODg79jzKfUderVw8PHjzIN//Ro0eoV69ecTZR5RgYGKBnz55o3ry52KkwxspgTmdD4SEAY8aMQadOneDhIsOczobllkOxumskJiaiYcOGePr0qcb8Ro0aISUlRSeJVXRRUVEYMGAAlErd9a1kjOme8JTxQ1kYPHgwkv49gsFZt3Tab/p1xTqjPnz4MFavXo2GDRsK8xo1aoRVq1bh999/10lihoaGWLduHWJjY5GamorDhw/DwcGhyHUWLVoEItKIyMhIneRXHB+0UcPDRYa6desK88r7k5gxVjZ5bdD7hhriwZZxGKz6XetjeRTHG9tOzM3N6cKFC6RSqSg8PJzCw8NJpVLRxYsXydzcXCftNZs2baIXL15Q7969qU2bNhQYGEg3btwgqVRa6DqLFi2ie/fukZ2dnRA2NjZlbh8qbXi4yChhvhXdP7aVunXrprWrwxwcHOUbNWvWpKU9qhEtsqAlHkY62ccbalDxN9S7d2+aPXs2zZ49m3r16qWzg2JhYUGZmZk0evRoYZ6joyOp1eoCL2rmxaJFi0p9cVMXhRoA9aonp6QvbGiZpykXaQ6OChgymYzO7viS0hbb0xIPI539HWutUJdX9OjRg4go39nw7du3afHixYWut2jRIkpLS6MXL17QkydPaN++fVS3bl1tHKQyxRIPI51+EnNwcOguPFxklDjfiqb2bSpM66JYl7rXx8WLFzUerbVs2TKNMUCsra3zXWDUBnt7e2RnZ+e7BTs6Ohr29vaFrnf16lWMHz8e/fr1g6+vL+zt7XHp0iXUqFFD6zkWl4eLDH5vybH6uhxzetTA2K7OouXCGCs5t9oyDN2XjB+O3wVQfv2mX1doFVer1WRraytMJyUlaZyh1qxZk7Kzs4v9ifHll1/Sm7i7u9OoUaNIpVLlW//06dO0ZcuWYu/P1NSUoqOjaebMmQW+7uvrSwqFghQKBYWGhurkkzjvk9fc3JxO/7iYEuZbcfMHB4eexpzOhsLfp1QqpbVr19Kn3p1oTmfdP7y6qDPqEo2mJJFISrJ4PmvWrMGePXuKXCY8PBwdO3aEgYEBbGxsNM6q7ezscP78+WLvLy0tDXfu3NHorfIqf39/+Pv7A8h9sKS2udWWvXJ1OAW9fZfA3Vmq9UHFGWPakdcVz/uAEtfjTVBX8hzvNXiAodfE/3sttIq/fkadnJxcpjPq4kbexcRRo0YJ8xwcHN54MfH1MDIyooiICFqwYEGZPs20Hf3796eWLVuW2/44ODiKH3nfhHV54bCgKHUbdV5f5Nfn6VpycjK2bduGFStWoFevXmjdujV2796Nmzdv4tSpU8Jy9+7dw/vvvy9Mr1y5Et27d4eLiwvat2+PAwcOwNTUFDt37tR5zsVlaGiIKVOmYMSIEWKnwhgrwNVoQ1yXtsZCdyNsDtbNU8VLqsimD4lEgj179iAzMxMAUK1aNfj7+yM9PR1A7pCeuvLxxx8jOzsb+/fvh7GxMU6fPo333nsPOTk5wjKurq6wsbERph0dHbFv3z7Y2NggNjYWV65cQceOHREeHq6zPEsqKysL06ZNQ2xsrNipMMYK4NvbFb0snmLLzerweyuxXIYwfRMJck+tC7R9+/ZibWTixInaykc0CoUCbm5u5bpPY2NjtGnTBpcuXSrX/TLGCubhIsMvXsaYesoIh/6NF6bL407EN9Ug0duE9CHKs406L2bOnEmXLl3SuA7AwcEhTrRu3ZpWDnHM1ybt4SKrWL0+mHZt2bIFgYGB3AzCmB6YOXMmDAwMcPawj8Z8XTxVvKSKNSgT0w2lUol///0XANC0aVOYm5uLmxBjVcSrQ5fm+XXFDITu+ahcOkyUFBdqPWBubo7NmzdjxowZYqfCWJWQ11+6V31DDBo0CB4uMmzvn4NT9xLFTq1A3PShB1JSUjB//nzcvn1b7FQYqxLybgP/fXR1pDV1hHFTcwz7OUX0Jo7C8Bm1nrh06RKSk5MhkUjQq1cvsdNhrNILClNj7aU01H64B+uvlP/40iXBhVrP9O7dG9988w06duwodiqMVVojRozAkFbW8HtLjqVnM+H3ljxfm7U+4aYPPXPy5EmkpaXhypUrYqfCWKVka2uLz8d0R73rQRi0O/dMOjAsu9z6S5cGn1HrobwbYGrVqoWRI0eKnA1jlUtsbCyObv0Sg/ckCkVZrKFLi4sLtR4bNmwYpkyZAmtra7FTYazCmzBhAvr16wcA+PyvGASGZmu8HhSmxspLWWKk9kZcqPXYpk2bMHbsWMTHx4udCmMVyuv9pGUyGYa2tsZMt7IN1SwWLtR6jIjw4sULAMCgQYOwZMkSyGT6+dWMMX2S10+6R10DGBoaopsT0PH5D/hmz0mxUysVvphYQVhbW8PKyqrMD29grCrIa3M+Os4WNw3fQv3EIHj/mq6XFwqLg8+oK4iAgABh6FcTExNut2bsDYLC1Pgr1h4ds87pzbjSpcWFugLJG4v7iy++wPbt2yGXy0XOiDH94+LigiZNmsDDRQZ3k0cVop/0m3DTRwUUEBCABg0aQKVSiZ0KY3pnyZIlcFA9Qcv73wv9ovW9n3RxiD4OrD6EGONRayNatWpFvr6+JJFIRM+Fg0OskEgkwt+As7MzLe1bQ7RxpUsbpX5mItN/Hh4e6NevH6pVqyZ2KoyJQi6XY/Xq1Zg6dSoA4OnTp1h4PCHfmbM+95N+Ey7UFdzatWsxfvx4KJVKSCQSrBjskK8tzsNFhjmdDUXKkDHtKGgMaQ8XGT52kyAyMrLSP4BD9FN+fYiK2vTxaowePZru/vUDxc2tLnzt83CRlesj7zk4dBWv/i5LJBJaNqU/xX1qUWl+t/lRXFXE77//jszMTMRfT8UvXsbYekONKW1kFfoCCmN58vpG/+JljF335JhR5xbWPO+GoLCjYqemc1yoK5H09HQcPHgQALDrrhyfd5ViR6gdgsIeipwZY9qRZNUSm4PvYqG7FKuvKDH3z8pfpAFuo66UPFxkGNdcjSuG3fGOfUyF7j/KWJ7+/ftj3zI/zOhojKVnM+HjmlmlfrdFb5vRh6gMbdRA/jbpvOmNs9+lESNGiJ4fB0dJw8rKigBQr/qG9PKzGtSjrgEBle/6C3fPq0Lcamu2SQeFqTHqt0y0tFGjVq1aImfHWMksWbIEW7ZsgVwuR1s7YMjeJGF4Un0fQ1rbRP8k0YeoLGfURYVMlnvm0ahRI5o5cyaZmJiInhMHx+thbGws3LzSqVMnGjJkSJW4oYvPqBkAQK3OPctu164d+vXrB6mUf/xMHIX1iV7SpwZ+++03DBs2DABw+fJl/P777yAiMdLUG/yXWgXt27cPQ4cORWpqKgDg008/RatWrUTOilUleeNF5xXrd1pY4hcvYwSGJOHUqVO4e/euyBnqF+6eV0Wlp6cDAGxsbNCtWzc8evQI//33n8hZsari1T7R16gFPC2fYchPuc8wDPr+e7HT0zt6e0bt6+uLM2fOIDExEUQEZ2fnYq03bNgw3LlzBxkZGbhz5w6GDBmi20QruLi4OAwfPhyHDx8GkDt2yLJly2Bubi5yZqwyc3R0xPV4E2wOVmGg2R38FWuPU48r5jgc5UFvC7WJiQlOnDiBxYsXF3udjh07Yv/+/fjpp5/QunVr/PTTT/j111/Rvn173SVaCWRlZQnt1zVq1ICDg4PQLGJoyGOEMO2qUaMGfv31Vyyd0Bt+b8mx9GwmuhjeRxeHHLFT02uiX+0sKtq1a0dERM7Ozm9c9ueff6YTJ05ozDt58iTt3bu3TFdcq1rkXWGXSqV06NAh8vX1FT0njoodderUof79+wvT897tTLGvjNNR2fpElyaqTK+PTp064cSJExrzjh8/js6dO4uUUcWUd4Xd0NAQp06dwu3btwEA1apV44uOrFRGjx6NOXPmCMPxqp8FY8QvaRr9/atSn+jSEP2TpKgoyRl1ZmYmjR07VmPe2LFjKSMjo8DlfX19SaFQkEKhoNDQUNHfq77HiBEjSKFQUMOGDUXPhUO/o379+rRjxw7h79ba2lq4w5Cj4NCbM+ovv/wSRFRkuLu7l1s+/v7+cHNzg5ubG+Li4sptvxXVkSNH8Nlnn+Hhw9xBnsaOHYspU6aInBUrb4X2gfa0gpOTEwAgISEBcrkcNWrUAADEx8cjMTGx3HOtLMq1e96aNWuwZ8+eIpcJDw8v9fajoqJgZ2enMc/Ozg5RUVGl3ib7n4yMDJw8eVKYdnFxgZmZmTDt6uqKkJAQ4SG8rHLK6wOdN1SBh4sMv3gZI/Stz2DyDJgzZw4SExMxZswYsVOtVEQ/5S8qSnox8fjx4xrzjh8/zhcTdRh5Fx6rV69Oly9fJj8/P43X53Q2rHDPruN4c3i4yOjlZ9Z08csBwkXATp06kYODg+i5VdR4Qw0SP8GCws7Ojlq1akWjRo0iIqL+/ftTq1atNNq5Tp06RcuWLROmO3XqRCqViubOnUuNGzemefPmUVZWFrVv376sB4njDSGTycjd3Z0cHR0JyG2j3LZtG43q6FjgaH5V+ep+RY3q1auTl5cXGRjkjl73q19rokUWtMzTVPTcKkNUyEK9aNEiKsi4ceOEZUJDQykgIEBjveHDh9O9e/coMzOT7t69S0OHDtXGQeIoYbi5udGePXuoevXcx4IlzLOk/dNaUswccy7SFSgsLS3JzMyMAFCXLl1IoVCQm5ub8IG7xMOIP3i1FBWyUOvZQeIoY5z8zINokQUt7WFEAKhFixbC2TeHfkVec5aNjQ1dvXqVRo0aRUDut6b69esXOuY5F+uyhd70+mBVk4eLDK2ygvHdNSmmtZPDw0WGuXPnatx16uzsDJmM+9BqW2E9NAp7Kv0PP/yAOXPmAMgdXmDVqlW4dOkSgNzRFx8/flzgmOfcB1r3RP8k0YfgM2rdRGFnX8Pb1iRXV1cCcs/UAgMDad68ecJ6devWrRJjEIt1/POmZ86cSd9//72w/LRp02jIkCGi510Vg5s+yn6QOEoZxen1IZPJqEePHtSkSRMCQLa2tqRQKGjkyJEEgAwNDal58+Ykl8tFfz8VMXrUNRDak5M+t6Ez2xYLr40ePZpmzpwpeo4cXKi1cZA4yjFMTExowIABVLt2bQJyu2gqFArq1KkTAaCaNWtSz549ydjYWPRc9S0kEgk5OTkJPTP69+9PQUFB9E1fc6JFFrRrgistXLiQqlWrJnquHJrBhbrsB4lDxDA3N6cePXoIjw4bOnQoKRQKoZC7ubnRnDlzhN4Jb2oyqUx9uy0tLWnIkCFUs2bN3Pfh4UEKhYKaNm1KAKhp06a0Ze5oiv3Ugnto6HlwoS77QeLQozAwMKAmTZoIBdnLy4tOnjwpPBNyypQpdOzYMZJKpQTk9jDp2rWrsL7YvRZK8kEhlUqpUaNGZGtrS0BuT4wNGzZQx44dCcjtr65QKKhv374E5D6xe/DgwcL9BmK/V47iBxfqsh8kjgoU3bp1oxkzZgjTS5YsoaNHjwrTc+fOpcNrZgvttomf1aB573YWXq9Vq5ZQGHURA5qaU+wr/clnD29P8XMtycNFRjKZjNavX09vv/02ASBTU1NSKBTCYGMmJia0c+dO6t69OwG5hbxWrVqFfouoTN8eKntwoS77QeKowGFmZkZ169YVpocPH05+fn60xMOIaJEFXft2MG3btk14ffPmzfTjjz8K0xs3bqQVK1YI0wsXLtT4IFiwYIHGmN1ff/01TZ48WZjetm0bvf/++8L06dOnacvc0cIHReZXjrTp01HC6z/88AMNHjxYmHZ3dyd7e3vRjyOHbqOoGsTPTGSVXmpqqvDEGgA4ePBg7mhvXsZYejYTfm+dxoLt/xts6scff9To033hwgVkZf3vMVEqlUpj2sDAQGP57Oxs4Yk5AHDnzh1EREQI02vWrMGzZ88QbanCQncjrLmWiS9O/yG8PnXqVI38z549W9q3zioR0T9J9CH4jLrqhD602/It2ByvB9+ZyNgrxL6zLm9YUO8DSiwKyhSexv36HYSM5eGmD1blrLyU/2nXQWFqoXDrWlEfFOWVA6tYuFAzVs7E/qBgFQ83fTDGmJ7jQs0YY3qOCzVjjOk5LtSMMabnuFAzxpie40LNGGN6jgs1Y4zpOS7UjDGm57hQM8aYnuNCzRhjeo4LNWOM6Tku1Iwxpue4UDPGmJ7jQs0YY3qOCzVjjOk5LtSMMabnuFAzxpie40LNGGN6ToLcp9xWecnJyXjw4IHYaVRINjY2iIuLEzuNComPXelUxuPm7OyMmjVrFvq66I9J14co6lHtHHzs+NjpV1S148ZNH4wxpue4UDPGmJ7jQv3/tm7dKnYKFRYfu9LjY1c6Ve248cVExhjTc3xGzRhjeo4LNWOM6Tku1K/ZunUrHj16hPT0dMTExOD333+Hq6ur2GnpPSsrK6xbtw737t1Deno6wsPDsWnTJtSoUUPs1CoEX19fnDlzBomJiSAiODs7i52S3vLz88OTJ0+gVCoRHByMrl27ip2SznGhfk1wcDDGjx+PJk2aoG/fvpBIJDh16hQMDAzETk2v1a5dGw4ODvj000/RokUL+Pj4oHv37ti3b5/YqVUIJiYmOHHiBBYvXix2KnrN29sba9euxbJly9CmTRtcunQJx44dg5OTk9ip6Zzonbn1OVq0aEFERI0aNRI9l4oW/fv3J7VaTebm5qLnUlGiXbt2RETk7Owsei76GFeuXKGtW7dqzAsJCaFly5aJnpsug8+oi2BiYoIJEybg6dOnCAsLEzudCsfCwgKZmZlIT08XOxVWCcjlcrRr1w4nTpzQmH/ixAl07txZpKzKBxfqAvj5+SElJQVpaWno378/evXqhaysLLHTqlCqV6+OL7/8Ev7+/lCr1WKnwyoBGxsbGBgYIDo6WmN+dHQ07O3tRcqqfFSJQv3ll1+CiIoMd3d3YfmffvoJbdq0Qffu3RESEoJff/0VxsbGIr4D8ZT02AGAqakpjhw5ghcvXuDTTz8VKXPxlebYMVaQKnHDi7W1NWxsbIpcJjw8HEqlMt98uVyOxMRETJs2DXv27NFVinqrpMfO1NQUf/31FyQSCfr374+0tLTySFMvleb3rl27dggODoaLiwuePn2q6xQrFLlcjvT0dIwaNQoHDhwQ5m/YsAHNmzeHh4eHeMnpWJXoyhAfH4/4+PhSrSuRSCCRSGBkZKTlrCqGkhw7MzMzHDt2DBKJBP369avSRRoo2+8dy0+lUuH69evw9PTUKNSenp44ePCgiJnpXpUo1MVVv359DB8+HKdOnUJsbCwcHR0xb948ZGZm4s8//xQ7Pb1mZmaGEydOwMLCAkOGDIGpqSlMTU0BAAkJCVCpVCJnqN/s7Oxgb2+PRo0aAQCaNm0KS0tLhIeHIzExUeTs9MeqVauwe/duXLt2DRcvXsS0adNQu3ZtbNmyRezUdE70rif6Eo6OjvTXX39RdHQ0ZWZmUnh4OO3Zs4caN24sem76Hu7u7lQYd3d30fPT91i0aFGBx27cuHGi56Zv4efnR6GhoZSRkUHBwcHUrVs30XPSdVSJNmrGGKvIqkSvD8YYq8i4UDPGmJ7jQs0YY3qOCzVjjOk5LtSMMabnuFAzxpie40LNKh13d3cQEaytrQucFlNAQACOHDkidhqsAhK9MzcHx6sREBBQ4M0fly9fLtb6eTffWFtbEwCSy+VkZ2cn+vsCQBYWFlS9enWtbe/ff/8llUpFDRs2zPdaYGAgrV+/Pt/8cePGUUpKiujHgqP4wWfUTC+dPHkS9vb2GjFgwIBSbUulUuUbGlMsycnJSEpK0sq23NzcULNmTezatQuTJk3SyjaZfuJCzfRSZmYmoqOjNeLVMS8sLCywadMmREREQKlU4u7du/D29i5wW683fYwbNw4pKSl4++238eDBAyiVSpw5cwZ169YtUY4nTpzAyZMnhWlTU1OEhIRgw4YNha7zetNHYGAgNm7ciK+//hqxsbGIjo7GypUrIZFI3rj/SZMmYe/evQgICMB7770HmUxWovxZxcGFmlVIf/31F9zd3TFhwgQ0bdoUs2bNKtHDHYyMjLBo0SJMmDABnTp1gkwmw2+//VaiHMaNG4fWrVtj9uzZAIB169YhKytLmC6uMWPGIDs7G507d8YHH3yAjz/+GO+++26R65iYmGDkyJHYs2cPLly4AKVSibfffrtE+2UVB4+ex/RSv379kJKSojFv48aNmDdvHnr37o1OnTqhWbNmuH//PgAgNDS0RNuXy+X46KOPcOnSJQDA2LFj8eTJE/Tq1QunT58u1jYiIyMxefJk7N+/HxYWFhgzZgzat2+PjIyMEuVy9+5dLFq0CADw8OFD+Pr6olevXvj5558LXWfEiBF49uwZ/v33XwDAnj17MHnyZBw+fLhE+2YVAxdqppfOnTuHKVOmaMx7+fIlAKBNmzaIjIwUinRpqNVqXLt2TZgODw9HREQEmjZtWuxCDQCHDx/G3r17sWDBAsyZMwc3b94scS6vrxMREYGaNWsWuc7kyZOxe/duYXr37t347LPPUKtWLURGRpY4B6bfuOmD6aX09HQ8fvxYI7Q9CD8RlXkbRkZGcHNzQ3Z2Nho0aFCqbbw+VjcRQSot/E+zcePG6Nq1K5YtWwaVSgWVSoV79+7BwMAAEyZMEJZLTk5G9erV861vaWmptQuarHxwoWYVzo0bN1CrVi24urqWehsymQzt27cXpp2cnFC7dm3cu3evRNtZuXIljIyM4OnpiQkTJuCdd94pdU7FNWnSJFy5cgWtWrVC69athVi0aBEmTpwoLPfgwQO0bds23/pt27bFgwcPdJ4n0y7R+whycLwaAQEBdOLECbKzs9MIGxsbAkASiYQuXbpEd+7coT59+pCLiwv17t2bBg8eTED+ftSvT48bN46ysrLo6tWr1LFjR2rVqhWdOXOG/vvvPyEHNzc3unfvHrm5uRWaZ79+/SgzM5Pat29PAGj+/PkUHR1dZJ/tgIAAOnLkiDBdUF/n15d5NQwMDCgqKopmzJiR7zV7e3vKzs6mHj16EACqW7cupaen0/r166lly5bUqFEj+vjjjykrK4v69esn+s+Zo/jBZ9RML3l6eiIqKkojbty4AQAgIvTv3x8XL17Enj17cO/ePaxduxaGhobF3n5mZia+/vpr7Nq1C1evXoVUKsWwYcOE101MTODq6goTE5MC17exsUFAQAC++uoroa37m2++wb179xAQEFCGd160d955B7a2tgU+IzAqKgoXL17E5MmTAeReYO3evTsaNmyIEydO4Nq1axg5ciRGjBiBv//+W2c5Mu3jJ7ywKmfcuHHYsGEDzM3NxU6FsWLhM2rGGNNzXKgZY0zPcdPH/7VbxzQAAAAAgvq3toYHhHACzDlqgDmhBpgTaoA5oQaYE2qAOaEGmAtlney22dTXRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x576 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# We plot now the results\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(5, 8), sharex=True)\n",
    "\n",
    "ax1.set_title(\"Mean Anomaly Distances\")\n",
    "ax1.scatter(0, 0, color=\"yellow\", s=25)\n",
    "ax1.plot(sample_trajectory[:, 0],\n",
    "        sample_trajectory[:, 1],\n",
    "        alpha=0.8,\n",
    "        color=\"white\",\n",
    "        linestyle=\"dotted\")\n",
    "ax1.plot(mean_anom_10deg_trajectory[:, 0],\n",
    "        mean_anom_10deg_trajectory[:, 1],\n",
    "        alpha=0.8,\n",
    "        color=\"tab:green\",\n",
    "        linestyle=\"None\",\n",
    "        marker=\"^\")\n",
    "ax1.axis('equal')\n",
    "\n",
    "ax2.set_title(\"True Anomaly Distances\")\n",
    "ax2.scatter(0, 0, color=\"yellow\", s=25)\n",
    "ax2.plot(sample_trajectory[:, 0],\n",
    "        sample_trajectory[:, 1],\n",
    "        alpha=0.8,\n",
    "        color=\"white\",\n",
    "        linestyle=\"dotted\")\n",
    "ax2.plot(true_anom_10deg_trajectory[:, 0],\n",
    "        true_anom_10deg_trajectory[:, 1],\n",
    "        alpha=1,\n",
    "        color=\"tab:orange\",\n",
    "        linestyle=\"None\",\n",
    "        marker=\"x\")\n",
    "ax2.axis('equal')\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "ax2.set_xlabel(\"Eclip. x in AU\")\n",
    "ax1.set_ylabel(\"Eclip. y in AU\")\n",
    "ax2.set_ylabel(\"Eclip. y in AU\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e64333bb-dbd1-4699-b461-e871ee9ddf90",
   "metadata": {},
   "source": [
    "Long story short: we take the mean anomaly to \"enrich\" our orbits!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fff0aa9e-19b7-49a6-a4c3-ecc0bdbcfc7e",
   "metadata": {},
   "source": [
    "### Creating more NEOs\n",
    "\n",
    "So let's take our dataframe and drop the Mean Anomaly column. We replace it with an array of meananomalies (in 10 degree steps) and \"explode\" the dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "cc4af1a5-8a44-446f-b18b-a961ce8ab101",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, let's add an epoch:\n",
    "granvik_model_df.loc[:, \"epoch_et\"] = spiceypy.utc2et(\"2000-001T12:00:00\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "eb5b5858-ad52-4856-a9e9-072bb7d24471",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We remove the Mean Anomaly first and ...\n",
    "granvik_model_df.drop(columns=\"MeanAnom_deg\", inplace=True)\n",
    "\n",
    "# ... replace it with an array\n",
    "mean_anomaly_deg_array = np.arange(0, 360, 60)\n",
    "granvik_model_df.loc[:, \"MeanAnom_deg\"] = [mean_anomaly_deg_array for i in granvik_model_df.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ab0866b1-f970-4292-9050-a8a39ba37593",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Explode the list\n",
    "granvik_model_df = granvik_model_df.explode(\"MeanAnom_deg\").copy()\n",
    "granvik_model_df.loc[:, \"MeanAnom_rad\"] = \\\n",
    "    granvik_model_df[\"MeanAnom_deg\"].apply(lambda x: np.radians(x))\n",
    "granvik_model_df.reset_index(drop=True, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "4afee2d2-0aa9-410b-8a2c-7c042f7d084a",
   "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>SemiMajorAxis_AU</th>\n",
       "      <th>Ecc_</th>\n",
       "      <th>Incl_deg</th>\n",
       "      <th>LongAscNode_deg</th>\n",
       "      <th>ArgP_deg</th>\n",
       "      <th>H</th>\n",
       "      <th>Perihelion_AU</th>\n",
       "      <th>Incl_rad</th>\n",
       "      <th>LongAscNode_rad</th>\n",
       "      <th>ArgP_rad</th>\n",
       "      <th>Perihelion_km</th>\n",
       "      <th>SemiMajorAxis_km</th>\n",
       "      <th>epoch_et</th>\n",
       "      <th>MeanAnom_deg</th>\n",
       "      <th>MeanAnom_rad</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.574981</td>\n",
       "      <td>0.783617</td>\n",
       "      <td>33.520763</td>\n",
       "      <td>278.480591</td>\n",
       "      <td>75.952057</td>\n",
       "      <td>21.064367</td>\n",
       "      <td>0.557182</td>\n",
       "      <td>0.585048</td>\n",
       "      <td>4.860403</td>\n",
       "      <td>1.325613</td>\n",
       "      <td>8.335328e+07</td>\n",
       "      <td>3.852117e+08</td>\n",
       "      <td>64.183927</td>\n",
       "      <td>0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.574981</td>\n",
       "      <td>0.783617</td>\n",
       "      <td>33.520763</td>\n",
       "      <td>278.480591</td>\n",
       "      <td>75.952057</td>\n",
       "      <td>21.064367</td>\n",
       "      <td>0.557182</td>\n",
       "      <td>0.585048</td>\n",
       "      <td>4.860403</td>\n",
       "      <td>1.325613</td>\n",
       "      <td>8.335328e+07</td>\n",
       "      <td>3.852117e+08</td>\n",
       "      <td>64.183927</td>\n",
       "      <td>60</td>\n",
       "      <td>1.047198</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2.574981</td>\n",
       "      <td>0.783617</td>\n",
       "      <td>33.520763</td>\n",
       "      <td>278.480591</td>\n",
       "      <td>75.952057</td>\n",
       "      <td>21.064367</td>\n",
       "      <td>0.557182</td>\n",
       "      <td>0.585048</td>\n",
       "      <td>4.860403</td>\n",
       "      <td>1.325613</td>\n",
       "      <td>8.335328e+07</td>\n",
       "      <td>3.852117e+08</td>\n",
       "      <td>64.183927</td>\n",
       "      <td>120</td>\n",
       "      <td>2.094395</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2.574981</td>\n",
       "      <td>0.783617</td>\n",
       "      <td>33.520763</td>\n",
       "      <td>278.480591</td>\n",
       "      <td>75.952057</td>\n",
       "      <td>21.064367</td>\n",
       "      <td>0.557182</td>\n",
       "      <td>0.585048</td>\n",
       "      <td>4.860403</td>\n",
       "      <td>1.325613</td>\n",
       "      <td>8.335328e+07</td>\n",
       "      <td>3.852117e+08</td>\n",
       "      <td>64.183927</td>\n",
       "      <td>180</td>\n",
       "      <td>3.141593</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2.574981</td>\n",
       "      <td>0.783617</td>\n",
       "      <td>33.520763</td>\n",
       "      <td>278.480591</td>\n",
       "      <td>75.952057</td>\n",
       "      <td>21.064367</td>\n",
       "      <td>0.557182</td>\n",
       "      <td>0.585048</td>\n",
       "      <td>4.860403</td>\n",
       "      <td>1.325613</td>\n",
       "      <td>8.335328e+07</td>\n",
       "      <td>3.852117e+08</td>\n",
       "      <td>64.183927</td>\n",
       "      <td>240</td>\n",
       "      <td>4.188790</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",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4811995</th>\n",
       "      <td>2.591117</td>\n",
       "      <td>0.826896</td>\n",
       "      <td>22.001493</td>\n",
       "      <td>334.001678</td>\n",
       "      <td>286.224030</td>\n",
       "      <td>24.542505</td>\n",
       "      <td>0.448532</td>\n",
       "      <td>0.383999</td>\n",
       "      <td>5.829429</td>\n",
       "      <td>4.995552</td>\n",
       "      <td>6.709938e+07</td>\n",
       "      <td>3.876256e+08</td>\n",
       "      <td>64.183927</td>\n",
       "      <td>60</td>\n",
       "      <td>1.047198</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4811996</th>\n",
       "      <td>2.591117</td>\n",
       "      <td>0.826896</td>\n",
       "      <td>22.001493</td>\n",
       "      <td>334.001678</td>\n",
       "      <td>286.224030</td>\n",
       "      <td>24.542505</td>\n",
       "      <td>0.448532</td>\n",
       "      <td>0.383999</td>\n",
       "      <td>5.829429</td>\n",
       "      <td>4.995552</td>\n",
       "      <td>6.709938e+07</td>\n",
       "      <td>3.876256e+08</td>\n",
       "      <td>64.183927</td>\n",
       "      <td>120</td>\n",
       "      <td>2.094395</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4811997</th>\n",
       "      <td>2.591117</td>\n",
       "      <td>0.826896</td>\n",
       "      <td>22.001493</td>\n",
       "      <td>334.001678</td>\n",
       "      <td>286.224030</td>\n",
       "      <td>24.542505</td>\n",
       "      <td>0.448532</td>\n",
       "      <td>0.383999</td>\n",
       "      <td>5.829429</td>\n",
       "      <td>4.995552</td>\n",
       "      <td>6.709938e+07</td>\n",
       "      <td>3.876256e+08</td>\n",
       "      <td>64.183927</td>\n",
       "      <td>180</td>\n",
       "      <td>3.141593</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4811998</th>\n",
       "      <td>2.591117</td>\n",
       "      <td>0.826896</td>\n",
       "      <td>22.001493</td>\n",
       "      <td>334.001678</td>\n",
       "      <td>286.224030</td>\n",
       "      <td>24.542505</td>\n",
       "      <td>0.448532</td>\n",
       "      <td>0.383999</td>\n",
       "      <td>5.829429</td>\n",
       "      <td>4.995552</td>\n",
       "      <td>6.709938e+07</td>\n",
       "      <td>3.876256e+08</td>\n",
       "      <td>64.183927</td>\n",
       "      <td>240</td>\n",
       "      <td>4.188790</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4811999</th>\n",
       "      <td>2.591117</td>\n",
       "      <td>0.826896</td>\n",
       "      <td>22.001493</td>\n",
       "      <td>334.001678</td>\n",
       "      <td>286.224030</td>\n",
       "      <td>24.542505</td>\n",
       "      <td>0.448532</td>\n",
       "      <td>0.383999</td>\n",
       "      <td>5.829429</td>\n",
       "      <td>4.995552</td>\n",
       "      <td>6.709938e+07</td>\n",
       "      <td>3.876256e+08</td>\n",
       "      <td>64.183927</td>\n",
       "      <td>300</td>\n",
       "      <td>5.235988</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>4812000 rows × 15 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         SemiMajorAxis_AU      Ecc_   Incl_deg  LongAscNode_deg    ArgP_deg  \\\n",
       "0                2.574981  0.783617  33.520763       278.480591   75.952057   \n",
       "1                2.574981  0.783617  33.520763       278.480591   75.952057   \n",
       "2                2.574981  0.783617  33.520763       278.480591   75.952057   \n",
       "3                2.574981  0.783617  33.520763       278.480591   75.952057   \n",
       "4                2.574981  0.783617  33.520763       278.480591   75.952057   \n",
       "...                   ...       ...        ...              ...         ...   \n",
       "4811995          2.591117  0.826896  22.001493       334.001678  286.224030   \n",
       "4811996          2.591117  0.826896  22.001493       334.001678  286.224030   \n",
       "4811997          2.591117  0.826896  22.001493       334.001678  286.224030   \n",
       "4811998          2.591117  0.826896  22.001493       334.001678  286.224030   \n",
       "4811999          2.591117  0.826896  22.001493       334.001678  286.224030   \n",
       "\n",
       "                 H  Perihelion_AU  Incl_rad  LongAscNode_rad  ArgP_rad  \\\n",
       "0        21.064367       0.557182  0.585048         4.860403  1.325613   \n",
       "1        21.064367       0.557182  0.585048         4.860403  1.325613   \n",
       "2        21.064367       0.557182  0.585048         4.860403  1.325613   \n",
       "3        21.064367       0.557182  0.585048         4.860403  1.325613   \n",
       "4        21.064367       0.557182  0.585048         4.860403  1.325613   \n",
       "...            ...            ...       ...              ...       ...   \n",
       "4811995  24.542505       0.448532  0.383999         5.829429  4.995552   \n",
       "4811996  24.542505       0.448532  0.383999         5.829429  4.995552   \n",
       "4811997  24.542505       0.448532  0.383999         5.829429  4.995552   \n",
       "4811998  24.542505       0.448532  0.383999         5.829429  4.995552   \n",
       "4811999  24.542505       0.448532  0.383999         5.829429  4.995552   \n",
       "\n",
       "         Perihelion_km  SemiMajorAxis_km   epoch_et MeanAnom_deg  MeanAnom_rad  \n",
       "0         8.335328e+07      3.852117e+08  64.183927            0      0.000000  \n",
       "1         8.335328e+07      3.852117e+08  64.183927           60      1.047198  \n",
       "2         8.335328e+07      3.852117e+08  64.183927          120      2.094395  \n",
       "3         8.335328e+07      3.852117e+08  64.183927          180      3.141593  \n",
       "4         8.335328e+07      3.852117e+08  64.183927          240      4.188790  \n",
       "...                ...               ...        ...          ...           ...  \n",
       "4811995   6.709938e+07      3.876256e+08  64.183927           60      1.047198  \n",
       "4811996   6.709938e+07      3.876256e+08  64.183927          120      2.094395  \n",
       "4811997   6.709938e+07      3.876256e+08  64.183927          180      3.141593  \n",
       "4811998   6.709938e+07      3.876256e+08  64.183927          240      4.188790  \n",
       "4811999   6.709938e+07      3.876256e+08  64.183927          300      5.235988  \n",
       "\n",
       "[4812000 rows x 15 columns]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "granvik_model_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "294f793b-60e3-4737-bd7f-13738b3dd88f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Store the resulting dataframe\n",
    "data_dir = pathlib.Path(\"results_data/Granvik\")\n",
    "data_dir.mkdir(parents=True, exist_ok=True)\n",
    "dataframe_filepath = data_dir / \"enriched_granvik_model.parquet\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "22415a97-62de-46a0-afd2-15de4b577fbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save\n",
    "granvik_model_df.to_parquet(dataframe_filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "401696f7-28e3-4492-a702-d1efc5a3252c",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "\n",
    "We have all building blocks now to create \"normalized\" detection sky maps. We need to compute the positions of the NEOs and their corresponding brightness. Then, wel'll apply the KDE method to get a probability density function!"
   ]
  }
 ],
 "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
}