data_df_processing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T14:08:32.460402Z",
     "start_time": "2020-02-17T14:08:30.456787Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import datetime as dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T14:08:32.474195Z",
     "start_time": "2020-02-17T14:08:32.470982Z"
    }
   },
   "outputs": [],
   "source": [
    "base_f_loc = '/storage/silver/stfc_cg/hf832176/data/THEMIS/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T16:52:53.392323Z",
     "start_time": "2020-02-17T16:52:53.387842Z"
    }
   },
   "outputs": [],
   "source": [
    "spacecraft = 'the'\n",
    "\n",
    "f_loc = base_f_loc+spacecraft+'/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center> Read in Data </center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T16:52:54.779055Z",
     "start_time": "2020-02-17T16:52:54.508455Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/storage/silver/stfc_cg/hf832176/data/THEMIS/the/the_esa.pkl\r\n",
      "/storage/silver/stfc_cg/hf832176/data/THEMIS/the/the_fgm.pkl\r\n",
      "/storage/silver/stfc_cg/hf832176/data/THEMIS/the/the_sst.pkl\r\n",
      "/storage/silver/stfc_cg/hf832176/data/THEMIS/the/the_state.pkl\r\n"
     ]
    }
   ],
   "source": [
    "!ls {f_loc}*.pkl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T16:53:08.316030Z",
     "start_time": "2020-02-17T16:52:55.478084Z"
    }
   },
   "outputs": [],
   "source": [
    "fgm = pd.read_pickle(f_loc+spacecraft+'_fgm.pkl')\n",
    "sst = pd.read_pickle(f_loc+spacecraft+'_sst.pkl')\n",
    "state = pd.read_pickle(f_loc+spacecraft+'_state.pkl')\n",
    "esa = pd.read_pickle(f_loc+spacecraft+'_esa.pkl')\n",
    "# e_levels = pd.read_pickle(f_loc+spacecraft+'_sst_energy_levels.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center> Resample and Concatenate Data </center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T16:54:29.413384Z",
     "start_time": "2020-02-17T16:53:10.030917Z"
    }
   },
   "outputs": [],
   "source": [
    "sst = sst.resample('3T',loffset=dt.timedelta(minutes=1)).mean()\n",
    "state = state.resample('3T',loffset=dt.timedelta(minutes=1)).mean()\n",
    "fgm = fgm.resample('3T',loffset=dt.timedelta(minutes=1)).mean()\n",
    "esa = esa.resample('3T',loffset=dt.timedelta(minutes=1)).mean()\n",
    "# sst.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T16:54:29.909416Z",
     "start_time": "2020-02-17T16:54:29.904934Z"
    }
   },
   "outputs": [],
   "source": [
    "esa.columns = ['esa_'+i for i in esa.columns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T16:54:30.347979Z",
     "start_time": "2020-02-17T16:54:30.343525Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(445921, 17) (445920, 10) (445701, 3) (445919, 32)\n"
     ]
    }
   ],
   "source": [
    "print(sst.shape,state.shape,fgm.shape,esa.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T16:54:31.762093Z",
     "start_time": "2020-02-17T16:54:30.816163Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "DatetimeIndex: 445922 entries, 2007-09-15 00:01:00 to 2010-04-01 00:04:00\n",
      "Freq: 3T\n",
      "Data columns (total 62 columns):\n",
      " #   Column      Non-Null Count   Dtype  \n",
      "---  ------      --------------   -----  \n",
      " 0   esa_E_1     80907 non-null   float64\n",
      " 1   esa_E_2     253893 non-null  float64\n",
      " 2   esa_E_3     253893 non-null  float64\n",
      " 3   esa_E_4     253893 non-null  float64\n",
      " 4   esa_E_5     253893 non-null  float64\n",
      " 5   esa_E_6     253893 non-null  float64\n",
      " 6   esa_E_7     253893 non-null  float64\n",
      " 7   esa_E_8     253893 non-null  float64\n",
      " 8   esa_E_9     253893 non-null  float64\n",
      " 9   esa_E_10    253893 non-null  float64\n",
      " 10  esa_E_11    253893 non-null  float64\n",
      " 11  esa_E_12    253893 non-null  float64\n",
      " 12  esa_E_13    253893 non-null  float64\n",
      " 13  esa_E_14    253893 non-null  float64\n",
      " 14  esa_E_15    253893 non-null  float64\n",
      " 15  esa_E_16    173111 non-null  float64\n",
      " 16  esa_E_17    173111 non-null  float64\n",
      " 17  esa_E_18    173111 non-null  float64\n",
      " 18  esa_E_19    173111 non-null  float64\n",
      " 19  esa_E_20    173111 non-null  float64\n",
      " 20  esa_E_21    173111 non-null  float64\n",
      " 21  esa_E_22    173111 non-null  float64\n",
      " 22  esa_E_23    173111 non-null  float64\n",
      " 23  esa_E_24    173111 non-null  float64\n",
      " 24  esa_E_25    173111 non-null  float64\n",
      " 25  esa_E_26    173111 non-null  float64\n",
      " 26  esa_E_27    173111 non-null  float64\n",
      " 27  esa_E_28    173111 non-null  float64\n",
      " 28  esa_E_29    173111 non-null  float64\n",
      " 29  esa_E_30    173111 non-null  float64\n",
      " 30  esa_E_31    173111 non-null  float64\n",
      " 31  esa_E_32    173111 non-null  float64\n",
      " 32  E_1         328507 non-null  float64\n",
      " 33  E_2         328507 non-null  float64\n",
      " 34  E_3         328507 non-null  float64\n",
      " 35  E_4         328507 non-null  float64\n",
      " 36  E_5         328507 non-null  float64\n",
      " 37  E_6         328507 non-null  float64\n",
      " 38  E_7         328507 non-null  float64\n",
      " 39  E_8         328507 non-null  float64\n",
      " 40  E_9         328507 non-null  float64\n",
      " 41  E_10        328507 non-null  float64\n",
      " 42  E_11        328507 non-null  float64\n",
      " 43  E_12        0 non-null       float64\n",
      " 44  E_13        0 non-null       float64\n",
      " 45  E_14        0 non-null       float64\n",
      " 46  E_15        0 non-null       float64\n",
      " 47  E_16        0 non-null       float64\n",
      " 48  quality     328507 non-null  float64\n",
      " 49  Bx_gsm      192348 non-null  float32\n",
      " 50  By_gsm      192348 non-null  float32\n",
      " 51  Bz_gsm      192348 non-null  float32\n",
      " 52  pos_x_gsm   350880 non-null  float32\n",
      " 53  pos_y_gsm   350880 non-null  float32\n",
      " 54  pos_z_gsm   350880 non-null  float32\n",
      " 55  spin_ra     350880 non-null  float32\n",
      " 56  spin_dec    350880 non-null  float32\n",
      " 57  spin_per    350880 non-null  float32\n",
      " 58  spin_phase  350880 non-null  float32\n",
      " 59  pos_x_gei   350880 non-null  float32\n",
      " 60  pos_y_gei   350880 non-null  float32\n",
      " 61  pos_z_gei   350880 non-null  float32\n",
      "dtypes: float32(13), float64(49)\n",
      "memory usage: 212.2 MB\n"
     ]
    }
   ],
   "source": [
    "tha = pd.concat([esa,sst,fgm,state],axis=1)\n",
    "tha.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center> Additional Data </center>\n",
    "\n",
    "## GSM Radial Distance\n",
    "\n",
    "R$_E$ = 6378.1 km"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T16:54:36.349909Z",
     "start_time": "2020-02-17T16:54:32.220036Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x2b882fa34730>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEECAYAAAAxqm/oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhU5d3/8fd3MtkJkJ0AgYR9CWFJRBCkLqAgKCrSqq27P621Vqut1as+tY/drPrU2tZqrUutWi1VLEoFVASpKEsAgQAJBMISshMg+zKZ+/fHjC1iQraZOTOZ7+u6cpGczJz5ZAIf7rnnnPuIMQallFKBx2Z1AKWUUt2jBa6UUgFKC1wppQKUFrhSSgUoLXCllApQWuBKKRWg7L58sISEBJOWlubLh1RKqYC3ZcuWSmNM4unbfVrgaWlp5OTk+PIhlVIq4InIoba26xSKUkoFKC1wpZQKUFrgSikVoLTAlVIqQGmBK6VUgNICV0qpAOXTwwjVfxljaHI4aWxppdnhREQIsQk2AZtNsNuEyNAQRMTqqEopP6UF7kHNDidHTzRwpKqew1X1lJxs4FhtM8fqmqlyf5yob6ahpZUmh5OOlmK3CfQJtxMTEUqfcDt9I+0kxUSQ1Dec5L4RJLv/TIuPZkDfCGw2LXulgokWeDc0NLdSUF5LXmk1+aU15JfVsL+8lpLqxi+VcohNiI0KIz46jLjoMMYN7EtsVChRYXYi7DbCQ0OIDA0hzG7DGIPTQKvT4DSGllZDfbODmkbXR21TCycbWthTWs3He5uobXJ8KVNkaAhpCdEMS4gmPSGacQP7kjGwH6lxkTqKV6qX0gLvgDGGA5V1bD10nK2HT7Dt8HH2ltXgdBd1uN3GqOQYpg2LJzUuyvURG0lqXBTJfSMI8dKouLbJQVl1IyUnGjl4rI7CStfH7pJqVu4qpdUdsG+EnfED+5ExqC/ZaXGclRZHXHSYVzIppXxLfHlJtezsbOPvp9IbYzh4rJ5P9lXw732VbDpYxYn6FgBiIuxMHhLLpNT+jB0Qw+gBMQyNj/ZaSXdXY0sre8tqyD1aTW7xSXYVV7OnpJpmhxOAUcl9mJoex9T0eGaOSNBCV8rPicgWY0z2V7ZrgbumRNbtq2BNXjn/3lfJ0RMNAAzqH8mMEfFkDY1lypBYhif2Cdh55iZHKzuLTrKxsIpNhVVsOXSc2iYHIjAptT8XjE7i/DFJjEvpG7A/o1K9lRb4aU42tLAmr5yVuaV8vLeChpZWYsLtTB8ez7kjE5g5MpG0+KheO3/saHWSW1zN2vxy1uRXsKPoBMZAUkw4F48fwPzMFM5Ki/O7VxdKBSMtcFxTC6v3lLN0axEf763A4TQkxYRz0fhk5o5P4exhcYSGBOeh8ZW1TXycX8HqvDI+yiunscVJUkw4l0xIYUFmCllDY3vtf2ZK+btuF7iIvAgsAMqNMRmnfe8+4Akg0RhT2VEIKwrcGEPOoeMs3VrE8h0l1DQ6SO4bzmUTBzI3I4XJqf11yuA09c0OVu8pZ/mOYtbkV9DscDI0PoqvZ6dyVdZgkvtGWB1RqaDSkwKfBdQCfz21wEUkFXgeGANk+VuBn2xo4a0tRby64RAHKuuIDA1hXsYArpwymOnD43VqoJNqmxysyi1lSc4RNhZWYRM4b3QSX89OZfbYJOxB+opFKV9qr8A7PIzQGLNORNLa+NaTwP3Ash6n86BdxSd55bNDLPu8mIaWViYP6c/jV2VyyYQUosP1qMmu6hNuZ1HWYBZlDeZgZR1Lco7w5pYiPsorZ2C/CK4/J42rz0qlf5QeyaKUr3VqDtxd4Mu/GIGLyELgAmPM3SJyEMi2cgRujGFtfgXPrN3PpoNVRITaWDhxENdNH0rGoH4ef7xg52h1sjqvnL+sP8hnB44RGRrCoqxB3HhOOiOS+lgdT6lep0dvYp5a4CISBawBLjLGnOyowEXkNuA2gCFDhmQdOtTmlYG6paXVyfIdxfzp4wPkldYwsF8EN89MZ3FWKv2iQj32OKp9u4ureWl9Icu2F9PscDJ7bDLfu3AEmYP7Wx1NqV7DkwU+AVgN1Lu/PRgoBqYaY0rPtB9PjcCbHU7+nnOEZ9fu5+iJBkYl9+H2WcO5bNLAoD2KxGqVtU288tkhXlpfSHWjg6+NSuR7F44ga2ic1dGUCngeK/A2vncQH02hOFqdLN12lKc+3MfREw1MGdKf75w3ggvGJOmRJH6iprGFVzYc4vl/F1JV18w5w+O5d84ostO0yJXqrp4chfI6cB6QAJQBDxtjXjjl+wfxcoE7nYblO0v47Qd7OVBZR+bgftx30WhmjUzQY5P9VH2zg79tPMyf1h2goqaJOeOSuf/i0YxMjrE6mlIBJ2BP5Pls/zF+tnw3u0uqGTMghnvnjGLOuGQt7gBR3+zgpfUHeXbtfuqaHSzOSuWeOSNJ6RdpdTSlAkbAFfihY3X88r09rNpVxqD+kdw/dzSXZg7UqZIAVVXXzNNrCnjls0OIwG2zhvGd80YQGRZidTSl/F7AFHh1YwtPf1TAS+sPYg8R7jx/BLfMTCciVP+h9wZHqup5fFU+72wvZlD/SH48fyzzMgboKyqlzsDvC9wYw/IdJfzvu7s5VtfEVVMG88OLR5Okp233SpsKq/jJslzySmuYMSKen146XufHlWqHXxf44WP1PLQsl3V7K5gwqB+/uCJDjyMOAo5WJ3/bdJgnVuVT39zKrecO457ZI/XVllKn8csCb2l18ud/H+CpD/dhtwk/uHg0109P03VKgsyx2iZ+vTKPJTlFpCdE8+iVEzh7WLzVsZTyG35X4LuKT3Lfku3kldYwd/wAHr5snB6ZEOTWF1Ty4NKdHK6q55tnD+GBeWOIidAzapXymwLfsHETz6zdz1Or9xEbHcYvLs/govEDfJZB+beG5lZ+80E+L3xSSFJMBI8umsB5o5OsjqWUpfyiwCdMmmKG3fo7thed5NKJA3nksvHE6vUYVRs+P3KCH725g/yyGq6fPpQH543VQw5V0PKLAo9IGWnG3vE0P798AvMzU3z2uCowNba08sSqfJ7/pJDhidE8dfVkXV1SBSW/KPDE9HFm945tJMaE++wxVeBbX1DJfUu2U1nbxL0XjeL2WcP1jW4VVNorcJ8u3Tc0PkrLW3XZjBEJrLznXC4eP4DHVubzrec3Ul7TaHUspSyna6+qgNA/Kow/XDuZx6/KZNuR48z/3Sd8tv+Y1bGUspQWuAoYIsLi7FT+eecMYiLsfPP5DTy9pgCn03fTgEr5Ey1wFXDGDOjLO9+dySUTUnh8VT63vLyZ43XNVsdSyue0wFVA6hNu5/fXTOZnC8ezvuAYlz39CfmlNVbHUsqntMBVwBIRrpuext9vn0ZTi5Mr/7ie93ed8ap+SvUqWuAq4E0eEss7353J8KQ+3P7qFp5eU4AvD49Vyipa4KpXGNAvgiW3T+eyiQN5fFU+d72+jYbmVqtjKeVVdqsDKOUpEaEh/PYbkxgzoC+PrcrjSFU9L9x4Fgl99NwD1TvpCFz1KiLCHecN57nrsskvq2HRM59SWFlndSylvEILXPVKc8Yl87f/N42aRgeLnvmUrYePWx1JKY/rsMBF5EURKReR3FO2PS4ieSKyQ0TeFhG9fI7yO1OGxPLWHecQE2Hn2j9v0CNUVK/TmRH4X4C5p237AMgwxmQCe4EHPZxLKY9IT4jmrTvOYXRyDN9+dQuvbjhkdSSlPKbDAjfGrAOqTtv2vjHG4f5yAzDYC9mU8oiEPuG8fts0zhudxEP/zOVPH++3OpJSHuGJOfCbgRXtfVNEbhORHBHJqaio8MDDKdV1UWF2/nRdFgsyU/jVijz+7/18PVZcBbweHUYoIj8GHMBr7d3GGPMc8By4LqnWk8dTqidCQ2w8dfVkosJC+P1HBdQ1tfI/C8YiomuLK//U7HDywFs72v1+twtcRG4EFgAXGh3KqAARYhMevTKT6HA7L64vpL7ZwS+umKAXiFB+p9nh5DuvbeHDPeXt3qZbUygiMhe4H7jMGFPfzXxKWcJmE36yYBx3XTCCNzYf4d4ln9OqS9IqP3Jqef9s4fh2b9fhCFxEXgfOAxJEpAh4GNdRJ+HAB+6XnxuMMd/2RHClfEFEuO+i0USGhfDYynxCRHh88UQdiSvLnV7e101P4/p2btthgRtjrmlj8ws9i6iUf/jOeSNwOg1PvL8Xm014bFEmNi1xZZG2yvtMdC0UFfS+e8FIHE7Dbz/ch90m/PKKCVriyue6Wt6gBa4UAHdfOJJWp+H3HxUQYhN+fnmGHp2ifMbR6uR7r2/rUnmDFrhSgGtO/N45o2hpNTz78X5CQ2w8fOk4LXHldU6n4YGlO1m5q5SHLx3X6fIGLXCl/kNE+NHc0bS0Onnhk0Jio8K4e/ZIq2OpXswYw8//tYc3txTx/dmjuGlGepfurwWu1ClEhIfmj+VkQwtPfriX2OhQru/CiEiprvjd6gJeXF/IzTPS+d6FI7p8fy1wpU4jIjx65QRO1Lfw8Du76BcZysJJg6yOpXqZl9YX8uSHe7kqazAPze/eGcG6HrhSbbCH2PjDtZM5Ky2O+5ZsZ21++2fDKdVVb28r4n/f3c3F45N59MruH/WkBa5UOyJCQ3j+hmxGJcdwx6tb+fzICasjqV5gfUElP/zHDs4ZHs9TV0/GHtL9GtYCV+oM+kaE8vLNU0mICePWlzdzpEpXjlDdl1dazbdf2cLwxD48e10WEaEhPdqfFrhSHUiMCeelG6fS0mq48aVNnKxvsTqSCkAlJxu48cXNRIWH8NJNZ9E3IrTH+9QCV6oTRiT14bnrsjhS1cBtr+TQ5Gi1OpIKIDWNLdz00mZqmxy8dONUBvaP9Mh+tcCV6qSzh8Xz+OJMNhZW8aM3d+gFIVSntLQ6+c5rWykor+WP35zCuIF9PbZvPYxQqS5YOGkQR6rqeeL9vQyJj+beOaOsjqT83CPv7ubf+yp5bFEms0YlenTfWuBKddGd54/g0LF6frd6H6OTY5ifmWJ1JOWnXt1wiFc2HOK2WcP4+lmpHt+/TqEo1UUiws+vyCBraCw/+Md2dhWftDqS8kOf7T/GT9/ZxfmjE/nR3DFeeQwtcKW6IdwewjPfmkL/qFBu++sWKmubrI6k/MjhY/Xc8doW0hKieeqayV67UIgWuFLdlBQTwXPXZVNZ28R3Xt1Ks8NpdSTlB2oaW7j1r5sxBp6/Ptsjhwu2RwtcqR6YMLgfj12VyaaDVfz03V1Wx1EWM8bwg39sZ39FHX/85hTSEqK9+nha4Er10MJJg/j214bzt42H+fvmw1bHURb607oDrNpVxoPzxjBjRILXH08LXCkP+OHFo5k5IoH/WbaL3KP6pmYw+nR/JY+tzGN+Zgq3zOzaut7d1WGBi8iLIlIuIrmnbIsTkQ9EZJ/7z1jvxlTKv4XYhKeunkR8dBjfeW0rJxv0dPtgUnqyke+9vo30hGh+vSjTZ1dy6swI/C/A3NO2PQCsNsaMBFa7v1YqqMX3CecP106h+EQD9y35HKdTz9QMBl9cjLi+uZU/XZdFn3DfnV7TYYEbY9YBVadtXgi87P78ZeByD+dSKiBlDY3lx/PH8uGecp5dt9/qOMoHfvneHrYePsFjV2UyIinGp4/d3TnwZGNMifvzUiDZQ3mUCng3npPGgswUnliVz6f7K62Oo7xoZW4Jf/n0IDfNSGNB5kCfP36P38Q0rhV92n2tKCK3iUiOiORUVFT09OGU8nsiwq8XZZKeEM09b3xOVV2z1ZGUFxQdr+f+N3cwcXA/Hpw31pIM3S3wMhFJAXD/2e71powxzxljso0x2YmJnl3IRSl/FR1u5/fXTOFEfQs//Md2Xbmwl3G0Ornnjc9xGvjdNZMJs1tzQF93H/Ud4Ab35zcAyzwTR6neY9zAvjx4yRhW55Xz8qcHrY6jPOip1fvIOXScX1yRwdB4756scyadOYzwdeAzYLSIFInILcCjwBwR2QfMdn+tlDrNjeekccGYJH65Io/dxdVWx1Ee8GlBJX9YU8DirMEsnDTI0iziy5d22dnZJicnx2ePp5Q/OFbbxLyn/k1MhJ3ld51LZFjProOorPPF77JPhJ3ld80kKsw3hwyKyBZjTPbp2/VMTKW8LL5POE9+YxIHKut4ZLmulxKojDE8sHQnJ+pb+P01k31W3meiBa6UD8wYkcDts4bz+qYjfLC7zOo4qhve3FLEB7vL+OHFoxk/sJ/VcQAtcKV85t45oxib0pcHl+7UQwsDTNHxev733d2cnR7ns3VOOkMLXCkfCbPb+M3XJ3KyoZmH/rlTDy0MEE6na4lYYwxPLJ6IzUsXZ+gOLXClfGhsSl/umT2K93aW8s72YqvjqE54cX0hGw5U8fCl40mNi7I6zpdogSvlY7fPGsbkIf35ybJdlFU3Wh1HncG+shoeW5XP7LFJLM4ebHWcr9ACV8rH7CE2/m/xRJocrfzorR06leKnHK1O7l2ynT7hdn51pe+WiO0KLXClLDAssQ8PzB3D2vwK/r75iNVxVBv+/O9Cdh49yc8vzyAxJtzqOG3SAlfKItdPT2PasDh+8d4enUrxMwcqannyw73MHT+ASyakWB2nXVrgSlnEZhN+dWUmzQ4nDy/TE3z8hdPpOmEnwm7jkYXjrY5zRlrgSlkoPSGae2aPYuWuUlbmlnR8B+V1f9t0mE2FVTw0fxxJfSOsjnNGWuBKWezWc9MZl9KXnyzbpdfStFjJyQYeXZHHjBHxfnnUyem0wJWyWGiIjV8vyqSytolHV+RZHSdoGWN46O1cWp2GX13hn0ednE4LXCk/MGFwP249dxivbzrMhgPHrI4TlJbvKGF1Xjn3XTSKIfH+dcJOe7TAlfIT3589iiFxUTy4dCeNLa1WxwkqNY0t/Gz5biYM6sdNM/xnrZOOaIEr5Sciw0L4+eUZFFbW8ed1B6yOE1Se/GAfFbVN/PzyDEL8aK2TjmiBK+VHZo1K5JIJA/jDmgKOVNVbHSco7C6u5i+fFnLt1CFMTO1vdZwu0QJXys88NH8cNhEeWb7b6ii9ntNp+J9lufSPCuOHF4+2Ok6XaYEr5WcG9o/kexeO5IPdZazJK7c6Tq/25pYithw6zoPzxtA/KszqOF2mBa6UH7plZjrDEqP56bu79A1NLzle18yvVuwhe2gsi6b4/zHfbdECV8oPhdltPHJZBoeO1fOcvqHpFY+tyqe60cHPLs/wq4s0dEWPClxEvi8iu0QkV0ReFxH/Pu9UqQAyc2QC8zNTeFrf0PS43KMneWPzYW6YnsbYlL5Wx+m2bhe4iAwCvgdkG2MygBDgak8FU0rBQ/PHEmITfvGvPVZH6TWMMfxs+W5io8K4e/ZIq+P0SE+nUOxApIjYgShArxGllAel9Ivkjq8NZ+WuUj1D00NW5paysbCKe+eMol9kqNVxeqTbBW6MOQo8ARwGSoCTxpj3T7+diNwmIjkiklNRUdH9pEoFqf83axgD+0Xw83/txunUq/f0RGNLK79csYfRyTFcfVaq1XF6rCdTKLHAQiAdGAhEi8i3Tr+dMeY5Y0y2MSY7MTGx+0mVClIRoSHcP3cMuUerWbrtqNVxAtpL6w9ypKqB/1kwDntI4B/D0ZOfYDZQaIypMMa0AEuBczwTSyl1qssmDmRian8eX5VHfbPD6jgBqbymkafXFDB7bDIzRyZYHccjelLgh4FpIhIlrnUXLwT0nRalvMBmE36yYCxl1U08+7EeVtgdv3l/L02OVn48f6zVUTymJ3PgG4E3ga3ATve+nvNQLqXUabKGxrEgM4Xn1u2n+ESD1XECyt6yGpbkHOG6aWmkJ0RbHcdjejQJZIx52BgzxhiTYYy5zhjT5KlgSqmv+tHcMTid8OQHe62OElAeW5lPdJiduy4YYXUUjwr8WXylgkhqXBTXTR/KW1uL2FdWY3WcgLD5YBUf7inj2+cNJzY68NY7ORMtcKUCzJ3njyAqzM7jq/KtjuL3jDE8uiKPpJhwbg6gCzV0lha4UgEmLjqM22YN4/3dZWw9fNzqOH7tg91lbDl0nO/PGUVkWIjVcTxOC1ypAHTLzHQS+oTx6xV5GKMn97TF0erksVX5DEuMZnFWYK422BEtcKUCUHS4nbsuGMnGwio+3qtnOLdl6dajFJTXcv/Fo3vFSTtt6Z0/lVJB4JqpQxgSF8WvV+brKfanaXK08tsP9zIptT8Xjx9gdRyv0QJXKkCF2W3cd9Eo9pRUs3xnidVx/MqSzUcoPtnIDy4ajes8w95JC1ypAHZp5kBGJvXhd6v30aqjcMC1YNUf1hRwVlosM0bEWx3Hq7TAlQpgNptw9+yRFJTX8i8dhQPw+qbDlFU38f05o3r16Bu0wJUKeJdkpDAqWUfh4Bp9/3Htfs5Oj+Oc4b1jwaoz0QJXKsDZbMLdF46ioLyW5TuC+5oqr244REWNa/QdDLTAleoF5mUMCPpReH2zg2c/3s+MEfFMG9a7576/oAWuVC/wxSh8f0Vd0I7CX91wiMraZr4/OzhG36AFrlSvEcyj8MaWVp5bV8jMEQlkp8VZHcdntMCV6iVOHYW/F2RHpPxjSxGVtU3ceX7vWi62I1rgSvUi8zIGMCwxmmfW7g+aNVIcrU7+9PF+Jg/pz7RhwTP6Bi1wpXoVm0349qzh7C6pZt2+Sqvj+MS7O4opOt7AneeN6PXHfZ9OC1ypXubyyYMY0DeCZ9YWWB3F65xOwx/X7Gd0cgwXjEmyOo7PaYEr1cuE2W3cem46Gw5U9fr1wj/cU8a+8lq+c/5wbLbgGn2DFrhSvdI1U4fQLzKUZ9futzqK1xhjeHrtfobERTF/QorVcSzRowIXkf4i8qaI5InIHhGZ7qlgSqnuiw63c8M5aby/u4yC8t557czP9h9j+5ET3P61Yb12ve+O9PSnfgpYaYwZA0wE9vQ8klLKE248J43I0BCe/fiA1VG84tl1B0iMCWfRlN55tZ3O6HaBi0g/YBbwAoAxptkYc8JTwZRSPRMXHcbVU1P557ajFJ9osDqOR+0tq2Hd3gpumD6UiNDed63LzurJCDwdqABeEpFtIvK8iER7KJdSygNumZmO0xj++tkhq6N41IufFBJut3Ht2UOtjmKpnhS4HZgCPGOMmQzUAQ+cfiMRuU1EckQkp6JCr92nlC8Njo1ibsYAXt90mPpmh9VxPOJYbRNLtx1lUdZg4qLDrI5jqZ4UeBFQZIzZ6P76TVyF/iXGmOeMMdnGmOzExMQePJxSqjtunpHOyYYWlm49anUUj3h1w2GaHU5unpFudRTLdbvAjTGlwBERGe3edCGw2yOplFIekzU0lszB/XhpfWHAX/y4saWVVzYc5PzRiYxI6mN1HMv19CiUu4DXRGQHMAn4Zc8jKaU8SUS4eUY6+yvqWLcvsKcx39leTGVtM7eeO8zqKH6hRwVujPncPT2SaYy53BjTu0/7UipAXTIhhaSYcF5cf9DqKN1mjOHFTwoZMyCGc4YHxwUbOhKcR78rFWTC7Daunz6UdXsr2FcWmCf2rC84Rl5pDbfMTA+6RavaowWuVJC49uyhhNttATsK/8unB0noE8ZlkwZaHcVvaIErFSTiosO4csoglm4t4kR9s9VxuqToeD0f5ZXxjbNSCbcH74k7p9MCVyqIXD89jSaHkze3FFkdpUte33QYIOhP3DmdFrhSQWRsSl+yhsby2sbDAXPFniZHK3/ffIQLxyYzqH+k1XH8iha4UkHmW9OGUFhZx6f7j1kdpVNW5pZSWdvMddN09H06LXClgsy8jBRio0J5dUNgrI/yymeHSIuPYuaIBKuj+B0tcKWCTERoCIuzU3l/dxll1Y1Wxzmj3cXV5Bw6zremDQ3KK+50RAtcqSB07dQhtDoNb2w6YnWUM3p14yEiQm0szkq1Oopf0gJXKgilJURz7sgE3th8GEer0+o4baprcrBs21EWZA6kX1So1XH8kha4UkHqW9OGUnKykY/yyq2O0qZ/7SihrrmVa6bq6Ls9WuBKBakLxySRFBPOkhz/nEZ5Y/NhhidGM2VIrNVR/JYWuFJByh5iY1HWYNbkV1DuZ29mFpTXsPXwCa4+a4iue3IGWuBKBbHFWYNpdRqWbvOviz38ffMR7DbhiimDrI7i17TAlQpiwxL7cFZaLEtyjvjNmZnNDidvbT3KnHHJJPQJtzqOX9MCVyrILc5O5UBFHVsP+8dy/qv3lFFV18zXz9I3LzuiBa5UkJs/IYWosBCWbPaPBa7e2HyElH4RzBqp19DtiBa4UkEuOtzOgswUlu8opq7J2ivXl5xsYN2+ChZnDSZEz7zskBa4UoqvZ6dS19zKeztLLM3x9rajGANX6ZmXnaIFrpQia2gswxKj+UeOddMoxhje3nqU7KGxDImPsixHINECV0ohIlw5eRCbDlZRdLzekgy7iqvZV17L5ZP10MHO6nGBi0iIiGwTkeWeCKSUssbCSa7iXPZ5sSWP//a2o4SF2FiQmWLJ4wciT4zA7wb2eGA/SikLpcZFkT00ln9uO+rzY8IdrU6WfV7M+WMS6R8V5tPHDmQ9KnARGQzMB573TByllJUunzyIfeW17C6p9unjflJQSWVtE1fo9EmX9HQE/lvgfqDd9ShF5DYRyRGRnIqKih4+nFLKm+ZPSMFuE59Po/xz21H6RYZy/pgknz5uoOt2gYvIAqDcGLPlTLczxjxnjMk2xmQnJuqB+Ur5s9joMM4bncSyz4/S6vTNNEpdk4NVu8qYn5lCuD3EJ4/ZW/RkBD4DuExEDgJvABeIyKseSaWUsszlkwdSVt3ExgO+uejxytxSGlpadfqkG7pd4MaYB40xg40xacDVwEfGmG95LJlSyhKzxybTJ9zO2z5aoXD5jmIG9Y8ke6iu+91Vehy4UupLIkJDmJsxgJW5pTS2tHr1saobW/ikoJJ5GQN03e9u8EiBG2PWGmMWeGJfSinrLZw0kJomB2vzvXvgwUd7ymlpNcybMMCrj9Nb6QhcKfUV04fFExsV6vW1UVbklpDcN5zJqTp90h1a4Eqpr7CH2JibMYDVe8q8No1S5x7hXzx+ADZdebBbtMCVUm26ZEIKdc2tXptGWZtfQZPDybwMPXW+u1xdEDkAAA1USURBVLTAlVJt8vY0yorcEuKjw5iaHueV/QcDLXClVJu8OY3S2NLKmrxyLhqfrBdu6AEtcKVUu7w1jfLvfZXUNbcyV6dPekQLXCnVLm9No6zILaFvhJ3pw+I9ut9gowWulGqXN6ZRmh1OPtxdxuxxyYTZtYJ6Qp89pdQZeXoa5dP9lVQ3OvToEw/QAldKndEX0ygrcj0zjbIyt5TosBDOHZngkf0FMy1wpdQZ2UNszBmXzEd7ymly9GwaxdHq5P3dZVwwNpmIUF06tqe0wJVSHZqXkUJNk4P1BZU92s+mg1VU1TUzL0PXPvEELXClVIfOGRFPTLidFTtLe7SflbmlRITaOG+0XtzFE7TAlVIdCreHcOHYJD7YU0ZLa7tXUDwjp9OwMreUr41KJCrM7uGEwUkLXCnVKXMzUjhR38Kmwqpu3X/r4eOU1zTp0ScepAWulOqUr41KJDI0pNtHo6zILSUsxMYFY/XCxZ6iBa6U6pTIsBDOH5PIql1lOLt4wWNjXNMnM0cm0Dci1EsJg48WuFKq0+ZmpFBR08SWw8e7dL+dR09y9EQDc/XoE4/SAldKddoFY5IIs9u6fDTKitxSQmzCnLHJXkoWnLTAlVKd1ifczqyRCazaVYoxnZtG+WL6ZPqweGKjw7ycMLh0u8BFJFVE1ojIbhHZJSJ3ezKYUso/zc1I4eiJBnYUnezU7fPLaiisrNPpEy/oyQjcAdxnjBkHTAPuFJFxnomllPJXc8YmY7cJK3I7N43y3s5SROCi8Tp94mndLnBjTIkxZqv78xpgDzDIU8GUUv6pX1Qo04fHszK3pFPTKCtzSzgrLY6kmAgfpAsuHpkDF5E0YDKw0RP7U0r5t3kZKRw8Vk9eac0Zb7e/opa9ZbW69omX9LjARaQP8BZwjzGmuo3v3yYiOSKSU1HhnatbK6V8a864ZETocBplhftKPheP1wL3hh4VuIiE4irv14wxS9u6jTHmOWNMtjEmOzFRF7BRqjdIjAnnrLQ4VnZwVubyHSVkDY1lYP9IHyULLj05CkWAF4A9xpjfeC6SUioQzMsYwN6yWvZX1Lb5/b1lNeSV1nDZxIE+ThY8ejICnwFcB1wgIp+7Py7xUC6llJ/74rDAle1Mo7y7vRibwLwJOn3iLd1e09EY8wkgHsyilAogKf0imZTan5W5pdx5/ogvfc8Yw7vbi5k+PF6PPvEiPRNTKdVt8zIGsPPoSQ4dq/vS9q2HT3DwWD2XZur0iTdpgSulum3hpEHYbcIrnx360va/fnaQmHA7l+r8t1dpgSulum1AvwgumZDC3zcfobbJAUB5TSPv7SxhUdZgosP1yjvepAWulOqRW2amU9Pk4M/rDgDwu9X7cBq44Zw0a4MFAf3vUSnVIxNT+3PpxIE8vaaAouMNvLW1iBvPSSM9IdrqaL2eFrhSqsd+eUUGlTVNvLW1iIvGJfPAvDFWRwoKWuBKqR6LiQjl9dum0dDcSmRYiNVxgobOgSulPEbL27e0wJVSKkBpgSulVIDSAldKqQClBa6UUgFKC1wppQKUFrhSSgUoLXCllApQ0pmrSnvswURqgHyfPWDb+gEnLc4A/pHDHzKAf+TwhwygOU6VAFRanAH847kYbYyJOX2jr8/EzDfGZPv4Mb9ERJ4zxtxmZQZ/yeEPGfwlhz9k0BxfyZBjdV+4c/jFc9HW9mCcQnnX6gBu/pDDHzKAf+TwhwygOfyR3z4Xvp5C8Yv/UZVS/k/74r/aey58PQJ/zsePp5QKXNoX/9Xmc+HTEbhSSinPCfg5cBGZKyL5IlIgIg+4t70gIttFZIeIvCkifdq574Pu++WLyMVn2mc3MoiI/EJE9orIHhH5Xjv3vUFE9rk/bjhle5aI7HTv83ciIt3McYGIbBWRXBF5WUTafOPaUzlE5EURKReR3FO2PS4iee7fx9si0r+z+d3b00Vko3v730UkrBsZfioiR0Xkc/fHJd7McIYck0RkgztDjohMbee+nvp9pIrIGhHZLSK7RORu9/bF7q+dItLuNIUnnw9/0M6/EZ/2hUcZYwL2AwgB9gPDgDBgOzAO6HvKbX4DPNDGfce5bx8OpLv3E9LePruR4Sbgr4DNfbukNu4bBxxw/xnr/jzW/b1NwDRAgBXAvG4+F0eAUe7bPALc4uUcs4ApQO4p2y4C7O7Pfw38urP53d9bAlzt/vxZ4I5uZPgp8IPuPIfdyXCGHO9/8RwClwBrvfz7SAGmuD+PAfa6/16MBUYDa4FsXzwfVn+c4d+Iz/rC0x+BPgKfChQYYw4YY5qBN4CFxphqcI2CgUigrXmihcAbxpgmY0whUODeX5v77GoG4A7gEWOME8AYU97GfS8GPjDGVBljjgMfAHNFJAXXX6oNxvW356/A5d14LhYBzcaYve7bfODe5rUcxph1QNVp2943xjjcX24ABncy/0L37/AC4E337V7uToZO8liGM+QwQF/35/2A4jbu6snfR4kxZqv78xpgDzDIGLPHGNPRORkefT7aGf12aiTvodGvP/SFRwV6gQ/CNcL8QpF7GyLyElAKjAF+7952mYg80sF9291nFzMMB77hfpm8QkRGujNki8jznchQ1IUM7e1rAGA/5SXyVUCql3N05GZcI0dEZKCIvNdBhnjgxCn/AfQkw3fdL5NfFJFYizLcAzwuIkeAJ4AH3Tm8/vsQkTRgMrDxDLfxyvMhIiHA08A8XKPZa0RkHK5XZE8aY0YAx4Fb2rjvOOBqYDwwF/ijiIScYZ/t8Ye+8KhAL/B2GWNuAgbiGnF8w73tHWPMT3wUIRxoNK5Df/4MvOjOkGOMudVHGQyuv/hPisgmoAZotSAHACLyY8ABvObOUGyMaXMu2guewfWf6iSgBPg/CzKA65XZ940xqcD3gRfcObz6+3DP674F3PPFiLMtXnw+2hupdmYk7/XRrx/0RbcEeoEfxT2idBvs3gaAMaaV/04ldPa+Z9xnF/ZTBCx1b3sbyOxihsFtbD+TNvdljPnMGHOuMWYqsA7XHKg3c7RJRG4EFgDfdL/872yGY0B/+e+br93KYIwpM8a0uqe0/ozrH79PM7jdwH//Xvyjizm69fsQkVBc5f2aMWZpR7fvRI7uPB/tjVTbHMl7afTrD33hUYFe4JuBke55tDBco813RGQE/GdO6zIgr437vgNcLSLhIpIOjMT1BlGb++xqBuCfwPnu23yNtotzFXCRiMS6X9JfBKwyxpQA1SIyzf0zXA8s6+ZzkeR+LsKBH+F6w8mbOb5CROYC9wOXGWPqu5LfXfZrcE3/gKsAu5Mh5ZQvrwBy27iZVzO4FeP6+wCu0ee+Nm7jsd+H+3YvAHuMMb/pYlZfPB9t8tLo1x/6wrM8/a6orz9wvZO/F9c7wT/G9Z/SemAnrn+kr+F+lxnXL+eRU+77Y/f98jnl3fzT99nVDO5t/YF/uXN8Bkx0b88Gnj/lvjfjeklYANx0yvZsd/79wB9wH7PfjRyP43pZmI/r5TPezAG8jmuKogXXiOgW9z6PAJ+7P55133Yg8F5Hzzuud/g3uffzDyC8Gxlecf8uduD6B5bizQxnyDET2ILraIWNQJaXfx8zcU2l7Tjl+b8E139iRUATUIbrPwhv/k6mf/EY7q8fdH9U8t8jlL50m9Nve8rXq9y3bXOf/t4XnvzQE3mUUl7nnm7ZC1yIa4phM3AtrsM73zLGvCEizwI7jDF/PO2+44G/4ZpuGgisxjUClrb2aYzZ5YufyR8E+hSKUioAGNc893dxjZ73AEvcRfsj4F4RKcB1dMsL8OU5cPftlgC7gZXAncb1fkZ7+wwaOgJXSqkApSNwpZQKUFrgSikVoLTAlVJeIW0s5uXefpe4FjfbJSKPtXPfaeI6xf5zcS0G91OfhA4wvr6kmlIqePwF16GOf/1ig4icj+tsyYnGmKYvzlNow8vA140x292nzI/2dthApCNwpZRXmLYX87oDeNQY0+S+TVuLvAEk4TqGHvcRJ7sBRCTaPbLfJCLbRGShe/uNIrJMRNaKawneh73yQ/kZLXCllC+NAs51T498LCJntXO7J4F8ca0ff7uIRLi3/xj4yLiWhjgf18Jg0e7vTcV1GnwmsFjOsM55b6EFrpTyJTuudc6nAT8ElrhPYf8SY8wjuM46fR/XCT8r3d+6CHhARD7HtZZ5BDDE/b0PjDHHjDENuNabmenFn8Mv6By4UsqXioClxnUCyiYRcQIJ7jczJwP/WQ3RGLMfeEZE/gxUiEg8rrMvF5nT1jIXkbP56jrevf4kFx2BK6V86T+LvInIKFxXsak0xtxkjJn0RXmLyPxTRuYjcS2DfALXWZd3ffE9EZl8yr7niEiciETiWpZ2vU9+IgvpCFwp5RUi8jpwHq4RdhHwMK518V90H1rYDNxg2j4d/Dpc69jX41pD/pvGmFYR+RnwW2CHiNiAQlzLFINrca23cC3p+qoxJsd7P51/0FPplVIBz73efLYx5rtWZ/ElnUJRSqkApSNwpZQKUDoCV0qpAKUFrpRSAUoLXCmlApQWuFJKBSgtcKWUClBa4EopFaD+Pw+mYZPgArZJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "RE = 6378.1\n",
    "\n",
    "tha['r'] = [np.linalg.norm(i)/RE for i in tha[['pos_x_gsm','pos_y_gsm','pos_z_gsm']].values]\n",
    "\n",
    "tha['r'].iloc[:600].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GSM -> MLT\n",
    "\n",
    "If you thought THEMIS GSM$^{TM}$ was the same as usual GSM, guess again!\n",
    "\n",
    "X is positive towards the Sun.\n",
    "Z is positive towards Magnetic North.\n",
    "Y is positive towards... dusk.\n",
    "\n",
    "$\\theta_a = Tan^{-1}(Y/X)$\n",
    "\n",
    "$\\theta_b = Tan^{-1}(X/Y)$\n",
    "\n",
    "(X < 0, Y < 0): $MLT = 6 \\left(\\frac{\\theta_a}{\\pi/2}\\right)$\n",
    "\n",
    "(X > 0, Y < 0): $MLT = 6 \\left(\\frac{\\theta_b}{\\pi/2}\\right) + 6$\n",
    "\n",
    "(X > 0, Y > 0): $MLT = 6 \\left(\\frac{\\theta_a}{\\pi/2}\\right) + 12$\n",
    "\n",
    "(X < 0, Y > 0): $MLT = 6 \\left(\\frac{\\theta_b}{\\pi/2}\\right) + 18$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T16:54:40.296291Z",
     "start_time": "2020-02-17T16:54:36.934527Z"
    }
   },
   "outputs": [],
   "source": [
    "mlt_arr = []\n",
    "for x,y in zip(tha['pos_x_gsm'],tha['pos_y_gsm']):\n",
    "# for x,y in zip([-3,3,-3,3],[2,2,-2,-2]):\n",
    "    if y < 0:\n",
    "        theta_a = np.arctan(np.abs(y)/np.abs(x))\n",
    "        theta_b = np.arctan(np.abs(x)/np.abs(y))\n",
    "        if x <= 0:\n",
    "            mlt = 6*(theta_a/(np.pi/2))\n",
    "            mlt_arr.append(mlt)\n",
    "        elif x > 0:\n",
    "            mlt = 6*(theta_b/(np.pi/2)) + 6\n",
    "            mlt_arr.append(mlt)\n",
    "            \n",
    "    elif y > 0:\n",
    "        theta_a = np.arctan(np.abs(y)/np.abs(x))\n",
    "        theta_b = np.arctan(np.abs(x)/np.abs(y))\n",
    "        if x > 0:\n",
    "            mlt = 6*(theta_a/(np.pi/2)) + 12\n",
    "            mlt_arr.append(mlt)\n",
    "        elif x <= 0:\n",
    "            mlt = 6*(theta_b/(np.pi/2)) + 18\n",
    "            mlt_arr.append(mlt)\n",
    "    else:\n",
    "        mlt_arr.append(np.nan)\n",
    "            \n",
    "tha['mlt'] = mlt_arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-17T16:54:42.495234Z",
     "start_time": "2020-02-17T16:54:40.878174Z"
    }
   },
   "outputs": [],
   "source": [
    "tha.to_pickle(base_f_loc+spacecraft+'.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center> Plotting Data </center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-06T20:03:21.612231Z",
     "start_time": "2019-09-06T20:03:21.586842Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits import mplot3d\n",
    "from matplotlib import cm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-06T20:03:21.824926Z",
     "start_time": "2019-09-06T20:03:21.614796Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/users/hf832176/.conda/envs/tb_env5/lib/python3.6/site-packages/numpy/lib/histograms.py:824: RuntimeWarning: invalid value encountered in greater_equal\n",
      "  keep = (tmp_a >= first_edge)\n",
      "/home/users/hf832176/.conda/envs/tb_env5/lib/python3.6/site-packages/numpy/lib/histograms.py:825: RuntimeWarning: invalid value encountered in less_equal\n",
      "  keep &= (tmp_a <= last_edge)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHZ9JREFUeJzt3X+QVeWd5/H3V8Q0xVig+GPoht5mRzIVTWdwptVUrN1SnNjGbEKcwlmjCSSSJVbpbobkDzs/ZeNYxW4ZKR2zuiQSmskYQphxpSJrS3SyqUyNiqibFh1LQgg2l0gGhGgHXIHv/nHPgcPl3u7b955zzzn3fF5VXX3vc8899zn94/me57e5OyIiUjynpZ0BERFJhwKAiEhBKQCIiBSUAoCISEEpAIiIFJQCgIhIQSkAiIgUlAKAiEhBKQCIiBTU6WlnYCznnHOO9/T0pJ0NEZFc2bp167+6+7njHZfpANDT08Nzzz2XdjZERHLFzH5dz3FqAhIRKSgFABGRglIAEBEpqEz3AYiIpOXdd99lZGSEw4cPp52Vmjo6Opg1axaTJ09u6P0KACIiVYyMjHDmmWfS09ODmaWdnVO4O/v27WNkZIQ5c+Y0dA41AYmIVHH48GFmzJiRycIfwMyYMWNGUzUUBQARkRqyWviHms2fAoCISEGpD0BEpA6Xr3iK3QcOxXa+rulT+KeB+bGdrxEKACIZ1L+hH4ChhUMp50RCuw8cYueKj8Z2vp6Bx2I7V6PUBCSSQaXREqXRUtrZkJR9/etf59577z3+/Ktf/Sr33XdfbOdXABARyaglS5YwODgIwLFjx1i3bh033XRTbOcfNwCYWYeZPWtm/9fMtpnZfw3S55jZM2b2mpn90MzOCNLfEzzfHrzeEznXl4P0V82sP7arEBFpQz09PcyYMYMXXniBJ554gosvvpgZM2bEdv56+gDeAea7+9tmNhn4uZn9b+CLwEp3X2dmDwJLgAeC72+6+wVmdgPw34D/aGYXAjcAFwGdwE/M7L3ufjS2qxERaTOf+9znWLNmDb/5zW+4+eabYz33uDUAL3s7eDo5+HJgPrAhSB8EPhE8XhA8J3j9KisPVl0ArHP3d9z9V8B24NJYrkJEpE1dd911PP7442zZsoX+/ngbTuoaBWRmk4CtwAXAt4FfAgfc/UhwyAjQFTzuAl4HcPcjZnYQmBGkPx05bfQ90c9aCiwF6O7unuDliIgko2v6lFhH7nRNn1LXcWeccQZXXnkl06dPZ9KkSbF9PtQZAIJmmnlmNh14BHhftcOC79WmpvkY6ZWftQpYBdDX13fK6yIiaUhrzP6xY8d4+umn+dGPfhT7uSc0CsjdDwA/BT4ITDezMIDMAsIxayPAbIDg9WnA/mh6lfeIiEiFl19+mQsuuICrrrqKuXPnxn7+cWsAZnYu8K67HzCzKcCfU+7Y/UdgIbAOWAw8GrxlY/D8n4PXn3J3N7ONwMNmdg/lTuC5wLMxX4+ISNu48MIL2bFjR2Lnr6cJaCYwGPQDnAasd/cfm9nLwDoz+2vgBeCh4PiHgL81s+2U7/xvAHD3bWa2HngZOALcqhFAIif0b+inNFqic2pn2lmRghg3ALj7L4CLq6TvoMooHnc/DFxf41x3AXdNPJsi7a80WmJ48TC9g71pZ0UKQjOBRUQKSgFARKSgtBqoiEg9VvbCwV3xnW9aNywbju98DVAAEMmYsBNYq4FmzMFdsPxgfOdbPi2+czVITUAiGdC/of94wT+0cEj7AAgAW7Zs4QMf+ACHDx9mdHSUiy66iJdeeim286sGIJIB4QggkahLLrmEj3/843zta1/j0KFDfOpTn+L9739/bOdXABARybBvfOMbXHLJJXR0dMS6GQyoCUhEJNP279/P22+/zVtvvcXhw4djPbcCgIhIhi1dupQ777yTm266idtvvz3Wc6sJSESkHtO64x25M2385e7Xrl3L6aefzo033sjRo0f50Ic+xFNPPcX8+fGsTKoAICJSjxTG7C9atIhFixYBMGnSJJ555plYz68mIBGRglIAEBEpKDUBiaSkf8OJ/V2rLQHdObWT/g39mhQmiVEAEElJdKmHapPAhhYOaWloSZSagEQS1L+h/6Q7fZEsUQ1AJEFa0E2yTAFARKQO4Zadcemc2pl6/44CgIhIHeJesC8L/TvqAxARyagHH3yQefPmMW/ePObMmcOVV14Z6/lVAxBpsfGGf4qEbrnlFm655Rbeffdd5s+fzxe/+MVYz68AINJi0XbktNuAJR++8IUvMH/+fD72sY/Fel41AYkkLJzQJdKINWvW8Otf/5o77rgj9nMrAIgkbGjhkIaDSkO2bt3K3Xffzfe//31OOy3+4nrcJiAzmw2sBf4QOAascvd7zWw58J+A3waHfsXdNwXv+TKwBDgK/Bd3HwrSrwHuBSYB33X3FfFejohIMjqndsY6cqee/p/777+f/fv3H+/87evr47vf/W5seainD+AI8CV3f97MzgS2mtnm4LWV7n539GAzuxC4AbgI6AR+YmbvDV7+NvBhYATYYmYb3f3lOC5EJMsq1/VR52/+pNFf873vfS/R848bANx9D7AnePyWmb0CdI3xlgXAOnd/B/iVmW0HLg1e2+7uOwDMbF1wrAKAtL3KdX3U+StZMKFGJTPrAS4Gwl0JbjOzX5jZajM7K0jrAl6PvG0kSKuVLiIiKag7AJjZHwB/D/yVu/8OeAD4I2Ae5RrCt8JDq7zdx0iv/JylZvacmT3329/+tspbRERaw/2UIipTms1fXQHAzCZTLvz/zt3/IfjgN9z9qLsfA77DiWaeEWB25O2zgNIY6Sdx91Xu3ufufeeee+5Er0ckszQcNF86OjrYt29fZoOAu7Nv3z46OjoaPkc9o4AMeAh4xd3viaTPDPoHAK4DXgoebwQeNrN7KHcCzwWepVwDmGtmc4DdlDuKb2w45yI5M7RwiP4N/eoAzolZs2YxMjJCllsiOjo6mDVrVsPvr2cU0OXAp4FhM3sxSPsK8Ekzm0e5GWcn8HkAd99mZuspd+4eAW5196MAZnYbMER5GOhqd9/WcM5Fckidv/kxefJk5syZk3Y2ElXPKKCfU739ftMY77kLuKtK+qax3iciIq2jmcAiMUhq569w8pH6DiQJWgxOJAZJLfUQNhllYe14aT+qAYiIFJRqACIxia4Vk4Xt/kTGowAgEpNoga8mG8kDNQGJJECTviQPFABEEqA9ACQPFABERApKfQAiCdGSD5J1CgAiCdEoIMk6NQGJiBSUAoCISEEpAIg0SUs8S16pD0CkSaXREsOLh9POhsiEqQYgIlJQCgAiIgWlACAiUlAKACIiBaUAICJSUAoAIiIFpQAgIlJQCgAiIgWliWBFt7IXDu6Cad2wTJOZJkqzgCXPVAMouoO7YPnB8uPl08pfK6tsZ7iyt/ZrBVYaLbVk1U/tMCZJUACQsmXD5UCw/GA5KFQKA0W11yRx2mFMkjBuADCz2Wb2j2b2ipltM7MvBOlnm9lmM3st+H5WkG5mdp+ZbTezX5jZn0bOtTg4/jUzW5zcZcm4wjv6ad1p50REUlJPH8AR4Evu/ryZnQlsNbPNwGeAJ919hZkNAAPA7cBHgLnB12XAA8BlZnY2cAfQB3hwno3u/mbcFyV1iDb9iEghjVsDcPc97v588Pgt4BWgC1gADAaHDQKfCB4vANZ62dPAdDObCfQDm919f1DobwauifVqJB7TutXeL1IAExoFZGY9wMXAM8D57r4HykHCzM4LDusCXo+8bSRIq5Ve+RlLgaUA3d1qnkhFOBoobCaCE01FYXAIH2vkkEhu1R0AzOwPgL8H/srdf2dmNQ+tkuZjpJ+c4L4KWAXQ19d3yuvShHDIJ9TX9l+tcI+mhYFARHKprgBgZpMpF/5/5+7/ECS/YWYzg7v/mcDeIH0EmB15+yygFKRfUZH+08azLnWLjvVXu7+IBOoZBWTAQ8Ar7n5P5KWNQDiSZzHwaCR9UTAa6IPAwaCpaAi42szOCkYMXR2kSdLCDt+4m2vUVyCSa/XUAC4HPg0Mm9mLQdpXgBXAejNbAuwCrg9e2wRcC2wHfg98FsDd95vZncCW4Lhvuvv+WK5C0hEGFDUFieTSuAHA3X9O9fZ7gKuqHO/ArTXOtRpYPZEMiohIMrQWULuaaIeviBSOAkC70kSvxGkhOMk7BQCRBpVGSwwv1jwIyS8FAGledHJYNE2TxEQyTQFAmletoNfIIJHMUwBoN9FJXyIiY1AAaDdZ6fwNm4XUFCSSWQoAkoxaC8opGIhkhgKAJEuLx7XU5SueYveBQwB0TZ/CPw3MTzlHkmUKANI61ZqFkt6UvmCb3u8+cIidKz4KQM/AYynnRrJOAUBap9raQWGfxVi1g8pZzcuGx08LhSugtntTVHDtOzuAlW14fZIIBQBpvcpNZSrTosctGz65Yzs8JpoW3d+4Vgd4zE1R4SzgzDS5BD+PnoHH2MmN6eRBckcBQFpvvI1mQvUW1Cnc7YazgHsGHlOTi+TWuPsBiKQm2mcwVlqLaQ0gaReqAUh21VtTaDGtASTtQjUAEZGCUgCQ4snpVpadUzvp39CfdjakjagJSIonoa0su6ZPOd4RnMSIoKGFQ/QO5itoSbYpAIjEJFrgJz0iKKwJDC0cSvRzpL0pAIjkSNgMVBotpZ0VaQPqA5DiCvsCctQfEN7xaxiqxEE1ACmunC5UFwYB9QdIs1QDEIHcjgyC8gqgUO536Jo+JeXcSJ6oBiAC444Mig6/zFrzy+4Dh6CD40tSsDzV7EiOjFsDMLPVZrbXzF6KpC03s91m9mLwdW3ktS+b2XYze9XM+iPp1wRp281sIP5LEYlBjX6B0mjp+FdWRt5oXoA0q54awBrgfmBtRfpKd787mmBmFwI3ABcBncBPzOy9wcvfBj4MjABbzGyju7/cRN5F4hftF4iuMnp29u78NS9AmjVuAHD3n5lZT53nWwCsc/d3gF+Z2Xbg0uC17e6+A8DM1gXHKgBIdkWbhc7uzsydv0hcmukDuM3MFgHPAV9y9zeBLuDpyDEjQRrA6xXplzXx2SKpqdwDIHOCZqyfv+cc4KNp50YyrNEA8ABwJ+DB928BNwNW5Vinel+DVzuxmS0FlgJ0d6e35K/IceHS0yt7YdnwSdsu1pL0shBjCmous3I0tFXS0VAAcPc3wsdm9h3gx8HTEWB25NBZQDhlsVZ65blXAasA+vr6qgYJkZZaNgyDFdtNjrPXcCuXhRBpVEMBwMxmuvue4Ol1QDhCaCPwsJndQ7kTeC7wLOWawVwzmwPsptxR3NJ96zKzdZ/kxinr7QRNKzs7ALrH38s4QeHf85nvg50dN6a6QY7k17gBwMx+AFwBnGNmI8AdwBVmNo9yM85O4PMA7r7NzNZT7tw9Atzq7keD89wGDAGTgNXuvi32qxlDtNquOzKpxynr7QR3+j0Dj7FzedAEVLm/cYs2rAn/nnsHB2rvg1yDboYkVM8ooE9WSX5ojOPvAu6qkr4J2DSh3IlkwJhbQNaxnES0P6AyvdnCt/PdI/Rv6J/QCCXdDElIM4Fr0F2SwImx/80MAa31t1NP4Rv9O4wKRx8NjZTonax/Y2mM/nJqqOcuSUGi/YUFf1pDP+sZcVSPevKvv+fkZPVnqwAQUeufZKwqvKrSxRBXQVzLhAqIlZERSXV2/tbKf+VwVf09JyOrzW6FDAC1xmjX+ifJSrSWjIt2CEfTanQMN1z4Htx1csdvE8tB1PrbTnUeg7RMIQOAxmhLIqoV9GMME6238E1DNG+Xr3gqsU5sSVchA0BUXP9sumPKn7H21e3f0I8dOTu+gniciWNRWfvbqacTO6tt3FmUpbKi8AEgrh++ahX5M9a+uqXREm+9tqL5dv+wWWhaDBPHokEkA5rpPyhywMhSWVH4ACCSqEYnhq2sWHoCTgSRCUiyKamZQjurnaJFowAghRRO7hqrFpCIah3Fla8vGz61o7dBSY5cqqXWqLnKYyR9CgBSSKXREsOLh+kd7B2zLyB249UIopvQ5NREawZZahMvGgUAKbxoLSC6yFoqd6lxriWU0jpFE1WrTbxWP0GR+w/ipgDQQvrDzb7dBw7xx333AZ0M5f33U8c6RVkW7SeIDkXNy4S1zG8chAJAS6njK9t6Bh6ja/qU481D0nq1hmXn8WYp6dnjcVAASMBYS0dINnVO7YS++4KN1tPOTQKiw1Ez2hQE6j9oNQWABOiPMH/KBX/jSyqkqXNq5/hLQoeFfrSTOcOBoF61ZiwrGNRHAUAk5yYUvMJCP4d9AuNRMJg4BQAprJ6Bx5j6R9OZdNppaWdFYpal2bZZpgAghVXuoDvRSdc5tZPewd7au3+1k5wMEZVkKQBkgIaHtkb053zm+059vSUTwRJSVz9AVM6HiE5Ekltywvi7tmWZAkAGaHhociqDa/hz7h0cSDNbsctzJ3bSmtmSsx55GO5ZiwJASrKw5nsR5Pmfs2VyMkRU4qcAkBI180hmtPHIoDi1Y1OtAoC0nTxMwZf0TXQSWTs21SoASNtRs0+DCjYyqJmhou3ShDtuADCz1cB/APa6+/uDtLOBHwI9wE7gL939TTMz4F7gWuD3wGfc/fngPYuBrwWn/Wt3b8cJ95IS3fXHoEAjg5rVDs0/UF8NYA1wP7A2kjYAPOnuK8xsIHh+O/ARYG7wdRnwAHBZEDDuAPoAB7aa2UZ3fzOuC5HiqTXCR0TqM24AcPefmVlPRfIC4Irg8SDwU8oBYAGw1t0deNrMppvZzODYze6+H8DMNgPXAD9o+gqksNTUkyCNDCqERvsAznf3PQDuvsfMzgvSu4DXI8eNBGm10kUmRE09tUVnMjc9qa1gI4OKuo1l3J3AViXNx0g/9QRmS4GlAN3d+d0Wr1Fa3nZsuuuvLSz0Y50QVpCO4aL+nzUaAN4ws5nB3f9MYG+QPgLMjhw3CygF6VdUpP+02ondfRWwCqCvr69qkGhnWsRKMkUdw22t0WUQNwKLg8eLgUcj6Yus7IPAwaCpaAi42szOMrOzgKuDNBERSUk9w0B/QPnu/RwzG6E8mmcFsN7MlgC7gOuDwzdRHgK6nfIw0M8CuPt+M7sT2BIc982wQ1jqM9aCU+1SfS3CNYpkST2jgD5Z46WrqhzrwK01zrMaWD2h3Mlxtdq+26mZqNY1Vm7uISnRyKDYpd3np5nAknm6+8+Igo0MaoW0+/wUADKsXaabi0g2KQBkmO58k9G/ob8Yu34lRU1BbUMBQAqnNFpieLEKroapKahtaDdsEZGCUg0g504ZRZByfiZKSzvEb8L7AzeqILOE25kCQM6dMoqgI8XMNEBLO8SvZfsDa5Zw7qkJSESkoNq6BqDmhWzS70UkG9o6ABSteaFr+hQ4XG4KyvLyCUX7vaShZf0AkmtqAmojYYG/c8VHq66pk6ZwOYcwOEmyhhYOURotte4Dww7hlS3oe5DYtHUNoMjSXmOkku7625zmBuSSAkCbSnuNERHJPjUBSaFoGQiRE1QDKICsNQelSctAJEyTw3JFAaAA0moOytpwT939t4Amh+WKAkDBJF0bqCz0s9Txq7t/kZMpABRM0rUBjfYRyQ8FgAJT34BIsSkAFFi0wK/cd1fBQKT9KQAI0FzTUNY6e6VMy0HIeBQA5BQTbRpSu382tWxZaMktBQA5hZqGRIpBAUDGpGAgDdPm8ZmnACB1GysYiJxCC8RlXlMBwMx2Am8BR4Ej7t5nZmcDPwR6gJ3AX7r7m2ZmwL3AtcDvgc+4+/PNfL6kJ093//0b+tPOgkgmxbEY3JXuPs/d+4LnA8CT7j4XeDJ4DvARYG7wtRR4IIbPFhlXabTU2rXxRXIiidVAFwCDweNB4BOR9LVe9jQw3cxmJvD5UnD9G/pPuevvnNpZyHWAwqGgqQr7ArRhTOY02wfgwBNm5sD/dPdVwPnuvgfA3feY2XnBsV3A65H3jgRpe6InNLOllGsIdHd3N5k9KaJqd/tFHQufiaGg0Q7glb2n9gmokzg1zQaAy929FBTym83sX8Y41qqk+SkJ5SCyCqCvr++U10Ukx6oV9NWCQiUFiUQ0FQDcvRR832tmjwCXAm+Y2czg7n8msDc4fASYHXn7LEANsyJFV0/BrpFEiWg4AJjZVOA0d38reHw18E1gI7AYWBF8fzR4y0bgNjNbB1wGHAybikTi1jm1k97B3kK2+7el6EYz9R6vGsO4mqkBnA88Uh7dyenAw+7+uJltAdab2RJgF3B9cPwmykNAt1MeBvrZJj5bZExhm3/q7d8Sj4kW5qox1KXhAODuO4A/qZK+D7iqSroDtzb6eSKNUA3ghHA0UBgcK5+3FW1NWRfNBJa20b+hn9Jo6aRCvy0LtwmKNodFR0i19dyI8UYehQoeHBQApG1oy8fqokGwd7A3/XkBrTZWAV/wYakKACIFEq0dVdYICmmiw1LbLDgoAIgUSGWTWDhrWk1lEWMV8NUCw8peOLir/DhnAUIBQKTAMjFTOE+qDUed1g3LD5Yf52z0kQKAiLT3iKA45ejuvh4KACKivoCCSmI1UJGW69/QrzH/DSrqSqmiGoC0CQ0BbVx0YljYH9A5tVPNQY3I2TaYCgAiAtSeL6BAMAFhoR8OJc14IFAAEJGq1C/QhMpAAJkMBgoAklu6Q01O2CegINCkaktSZCgQKABIbqlwSo5WU01AWOhnaK6AAoCISCtlaKVSBQARkVaKFvgp1wY0D0ByTTt/JS9cL0jaj2oAkmvqAE6WOoMTFp03wIqWf7wCgIjUpM7ghKXcMawmIMklLf3QWp1TO9UMlKRp3ezsuLEcCFa2LtgqAEgulUZLav5poaGFQ2oGStKyYXoOP3zystItCARqAhIRyZIWNgupBiAidQlHXKkpqH2oBiAidVGHcPtRDUByRx3AIvFoeQAws2vM7FUz225mA63+fMk3LQAnEp+WNgGZ2STg28CHgRFgi5ltdPeXW5kPyZ7Kgr3a89JoSRuViMSo1X0AlwLb3X0HgJmtAxYACgAFVznEMCzsoztUaccvkXi1OgB0Aa9Hno8Al7U4D5JRlQW+7vRFkmXu3roPM7se6Hf3zwXPPw1c6u7/OXLMUmBp8PSPgVfrPP05wL/GmN0sK9K1QrGuV9favlp5vf/G3c8d76BW1wBGgNmR57OAk+r+7r4KWDXRE5vZc+7e11z28qFI1wrFul5da/vK4vW2ehTQFmCumc0xszOAG4CNLc6DiIjQ4hqAux8xs9uAIWASsNrdt7UyDyIiUtbymcDuvgnYlMCpJ9xslGNFulYo1vXqWttX5q63pZ3AIiKSHVoKQkSkoHIfAMzsejPbZmbHzKwvkt5jZofM7MXg68E08xmHWtcavPblYHmNV82srZZrNLPlZrY78ru8Nu08xa1oS6SY2U4zGw5+n8+lnZ84mdlqM9trZi9F0s42s81m9lrw/aw08xjKfQAAXgL+AvhZldd+6e7zgq9bWpyvJFS9VjO7kPKIqouAa4D/ESy70U5WRn6XSfQhpSayRMpHgAuBTwa/03Z3ZfD7zNTQyBisofx/GDUAPOnuc4Eng+epy30AcPdX3L3eyWK5Nsa1LgDWufs77v4rYDvlZTckH44vkeLu/w8Il0iRHHL3nwH7K5IXAIPB40HgEy3NVA25DwDjmGNmL5jZ/zGzf5d2ZhJUbYmNrpTykpTbzOwXQfU6E9XnGBXh91fJgSfMbGsw+7/dne/uewCC7+elnB8gJxvCmNlPgD+s8tJX3f3RGm/bA3S7+z4z+zPgf5nZRe7+u8QyGoMGr9WqpOVqeNdY1w08ANxJ+ZruBL4F3Ny63CUu97+/Blzu7iUzOw/YbGb/Etw5SwvlIgC4+5838J53gHeCx1vN7JfAe4FMdzg1cq3UscRG1tV73Wb2HeDHCWen1XL/+5sody8F3/ea2SOUm8HaOQC8YWYz3X2Pmc0E9qadIWjjJiAzOzfsCDWzfwvMBXakm6vEbARuMLP3mNkcytf6bMp5ik3wDxO6jnJneDsp1BIpZjbVzM4MHwNX036/00obgcXB48VArdp8S+WiBjAWM7sO+BvgXOAxM3vR3fuBfw9808yOAEeBW9y9smMmV2pdq7tvM7P1lPdVOALc6u5H08xrzP67mc2j3CyyE/h8utmJVwGXSDkfeMTMoFwGPezuj6ebpfiY2Q+AK4BzzGwEuANYAaw3syXALuD69HJ4gmYCi4gUVNs2AYmIyNgUAERECkoBQESkoBQAREQKSgFARKSgFABERApKAUBEpKAUAERECur/A1SRRw4J1P0UAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(tha['pos_y_gsm']/6378.1,bins=100,histtype='step',label='y')\n",
    "plt.hist(tha['pos_x_gsm']/6378.1,bins=100,histtype='step',label='x')\n",
    "plt.hist(tha['pos_z_gsm']/6378.1,bins=100,histtype='step',label='z')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-06T20:14:19.553702Z",
     "start_time": "2019-09-06T20:03:21.826669Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsvXuUZFld5/vZ+5wTz4x8Z9Yrq6sqq7q6Hl3dTVc13S1ty9LFAlFYAo6gzEUFdO5cvIOjsy6i4gzOFdDBkTU6M0u9zr2zWIreEQVBpxHk+gCBptt+0DQ0lZWPyszKd2ZkZrzP2XvfP06ckyciI1+RUdlV3fFdqyofEbHj5IkT3/jt3+/7+/6EMYY22mijjTYOBvLFPoA22mijjZcT2qTbRhtttHGAaJNuG2200cYBok26bbTRRhsHiDbpttFGG20cIOwdbm9LG9poo4029g6x1Q3tSLeNNtpo4wDRJt022mijjQNEm3TbaKONNg4QbdJto4022jhAtEm3jTbaaOMA0SbdNtpoo40DRJt022ijjTYOEG3SbaONNto4QLRJt4022mjjANEm3TbaaKONA0SbdNtoo402DhBt0m2jjTbaOEC0SbeNNtpo4wCxk8tYG21sCWMMWmvK5TKe52HbNlJKLMtCSomUEiG2NFtqo42XJcQOgynb1o5tbIIxBqUUnufVfB/cFhBtpVKhXC7T29uLZVltMm7j5YQtL/B2pNvGrlFPtkIIpJRorcPvoyiXy8zOzpLJZKhUKjVEG0TEbTJu4+WGNum2sSOMMXieh1Kqhmx3ghACIQSWZW1aD8DzPFzXrbmtTcZtvNTRJt02toQxBtd1mZqa4siRI7sm2wBCCBqlrwICrSfSejKOpioC8g7yxgE5t8m4jdsNbdJtYxO01jVphImJCY4dO7bl/bcivq1Id6d1tiLj+tSGMWbbyLhNyG3cimiTbhshtNZhGgHYc2Rbj72S7nbrRL8GqCfj+sdIKbFtu03GbdxSaJPuyxzGmDCNoLUG2DMxZbNZCoUCmUyGZDJZQ5KtIN2tsBMZB3K2sbEx+vv7yWQyIRnXpyraZNzGQaFNui9TBBpbz/OaIltjDIuLi4yOjhKLxUgkEszPz1MsFhFCkEqlcByHUqlEPp8nmUzuK2reC+rJWCkV5oSDv1spRaVSqXlcNE0RRMdtMm6j1WiT7ssM9WQbkMpeyHZubo6xsTEymQx33303iUQC13VDUtVaUywWWV5eZmVlhbGxMYrFIgDJZJJ0Oh3+Owgyri/INfpbo2QcvT/QMGfcVlS00SzapPsyQSON7V7IVmvNzMwMExMT9PT0cN9995FMJgHCHHAAKSXpdBopJdlslrvvvjtco1gsks/nyefzYWQMLw4ZR9EMGbflbW00gzbpvsSxVUPDbhFsw7/yla/Q39/P5cuXicfju3psfU43ION0Ol1zv63I2BiziYxTqdQtQ8ZBLrzd+NHGXtAm3ZcogoaG6elpUqkUnZ2deyIrz/OYnJzkxo0bGGN45StfieM4ezqG3RbStiPjICecz+dZXFykUCjc8mQMjRs/gg++TCbTJuOXMdqk+xJDQLaBhGp9fX1PTQSu6zIxMcHc3BzHjh3jwQcf5PHHH9+WcFul062HlJJUKkUqlWJgYCD8vTGmJjLejoyTyeRNVVDUYztFxerqKrlcjhMnTtTcFpBvu/Hj5YE26b5EEG1ogNoW3ECdsB3K5TLj4+MsLi5y/PhxHn744X1HjjdLMhaoI3ZLxsVikWKxSFdXV01kXN+efDMREGhAqNFjhnbjx8sJbdK9zdGooaE+v7gd6RaLRcbGxshms5w8eZI777yzIdnWF5F2i4OOMhuR8fPPP8/AwABCCPL5PMvLyxQKBbTWJBKJTWmKm0XGWutN57bd+PHyQ5t0b0PspaFhK9LN5/OMjo6Sy+U4deoU58+f3zFNsNc39K1CAEIIEokEmUyG/v7+8PfGmJqccSMyTqVSdHR0tISM93IOd9v4Uf+YduPHrY826d5GaKahoZ5019fXuXbtGuVymeHhYfr7+3d8QwZr7DXdcLM70vYLIQTJZJJkMtmQjAuFAvl8nqmpKfL5PFpr4vF4TWScTqd3TcbN7hbqjzn6tX79duPHrY826d4G2E9Dg5QS13XJZrOMjo6itWZ4eJje3t5dP/9O5BkcWywW29PjDgp7PYYoGff19dWsUy6Xw8h4enp6SzJOpVLYdu3bK3jtbhb2ojVeW1tDa01fX1+78eOA0SbdWxj7bWgwxlAoFJiZmWFlZYXTp0/T1dW15+PYijy11ty4cYOJiQmklCiliMViNcSzmyLeQaAVJBKkKRKJxI5kXCgUUEqFZJxKpSgWi2FDyUGi0TVTLBbRWrcbP14EtEn3FkRAtjMzM8RiMbq6uva0tTfGsLCwwNjYGAD9/f1cuHCh6eOpT1ForZmammJycpLBwUEuX74cvhld1yWXy5HP55mdnaVYLPL444/va1t+q2M7Mq5UKiEZr66usri4yI0bN2o+nIJ/9ZHxzYRSKsz71qPd+HFz0SbdWwj1ExrW19dJpVJ0d3fv+vGzs7OMj4/T2dnJpUuXyOfzrKys7Ou4gkhXKcXk5CTT09McPnw4bJjQWuO6LkIIYrEYvb29YfpifX2dK1euUKlUQjKORoIHqR44aAghiMfjxONxent7wwLd4OBgDRnPzMyQz+c37RRuJhlvl6NvpvGjTca7R5t0bwHUNzREq9D1vgaNEGzzr1+/Tk9PD694xStIJBLAxjZyv7h+/TqLi4scOXKEBx98cBMRbFckipLPdtvyRuqBjo6OF6Xj7GYgmo+PknGA+sh4dnaWfD6P53k4jhOqKIKve+0QrD+WZgqj0a/R44atyXgredvLFW3SfRGxVUNDgJ0aG5RSTE1NMTU1xeDgIFeuXNlUzNpJp7sdPM9jYmKC5eVl0uk0Dz300JZRaDNvou225YGUK5fLsbi4GH54BB1nARnvxhjnVijmAWHDw1bYjoxd192SjOsj492QsVKqZTuK3WqNGzV+RKVtLxdFRZt0XwTs1NAQwLKsTVpM8Mnw+vXr3LhxgyNHjmzri9AM6bquy/j4OPPz8xw/fpyBgQGOHj3aUsLdDttJuaIdZwsLCxQKBaDWpayjo6PGTP1mHGMzaFYyFqRtYrEYPT09NbdFI+O5ubldk3EzkW4zxx39GqBR48e1a9c4efJkSMIv5caPNukeEJqZ0FBPmJVKhYmJCebn5zl27BgPP/zwjtHKXki3UqkwPj7OwsICd9xxR9gKHMiLXmxs1XG2nWVkKpUK89oBmb9Yb95W6HTrsRsynp+fJ5/P47outm2TTqdZX18P0zb1u6ObjUZkvLq6im3bCCFe8o0fbdK9ydjPhIYgpxv1RYiS4W4QXMTbIRhps7y8zIkTJzhz5kzN+reK3nYrbOdSVigUeOGFF0LiLRaLoZFONApMJBI3/c17s3W6UWxFxkGaYm1tjWw2y/z8fA0ZR/8dJBlHUy/7afwIgprOzs4DOe5m0Cbdm4T9TmgA/w0yPz/P4uLitr4I22G7vHCpVGJsbIyVlRVOnjzJXXfdtadW4lsdUko6OjpIJBIMDQ2FpKyUCrvNVldXuXHjBqVSqYa8g3/xeLxlRLlTTvcg4DgO3d3dJBIJhoeHQ91wNGe8sLDA+Pj4lmTsOM6BR5e7afwA+NznPsczzzzDhz70oQM9vr2gTbotRiCtmp2dpa+vrymyzeVyjI6Osra2RiKRqNHB7hWNCLNYLIbrnzp1inPnzm27/q0e6e6E+mO3LItMJkMmk6n5fUDGuVyOlZUVpqamKJfLWJbVMArc62tyM9ILzaK+kBaQcb080XXd8ANqcXGRiYkJKpVKy87JflH/3lpdXW2qAegg0SbdFqG+e+w73/kOr3rVq/a0xtraGqOjo7iuy6lTpzh58iTj4+P7upCjpBuY3OTzeYaHh7lw4cKu1t5NiuJWx27TOY3I2PO8kIyXlpa4fv06lUplz1vyg0wv7ITdFtIcx6Grq2sTkXmeF0bG0XPSDBm38rysrq7uWtf+YqFNuvtEfUNDkPDfy0W0srLC6OgoAMPDw2EeLujr3w8C74Vnn32WYrG4a5Ob+jW2inTz+TzXrl0LCyGBlCv4+lKAbdt0dnZuyhNGiafRljw4B8GW/FZILwTYr3rBtu1dkfHk5OSm3UKgNQ7IuJXytdXV1U0m8bca2qTbJLZqaKi/z1bkZoxhaWkpHGF+5513bnpT77Y5YisEjmKFQoFz587R29vbVETRKEURkG1A5GfPnkVrHWpr5+fnww60p59+ukbO9VJpAd6KeKL50eA8eJ5HuVxGCEFXV1d4Hg6y9TeKm/UBsB0ZB2mK+tRNPB7HdV2Wlpb2nUdvR7ovQezU0BAgIMxGnVvz8/OMjY2RSqW4cOECHR0dDZ+r2QLW6uoqo6OjKKU4ffo0hUKhpvlgr4jmdKNke/r06XDdYLtd/4b7+te/zoULF0ICvnHjRtjyGm0BDrqsbpVIcD/YKj/6zDPP0Nvbi+d5NQ0O0dbf4Dy8WGR8s7DdbmFxcZFyuVxDxs0WNdfW1tqk+1LBbhsaAtSTrtY69EXo7u7mnnvuIZVKbfuce410s9ks165dA+D06dMtu/iEEBSLxTBFcebMmZqoeaciW70fQ/CYcrkcknHQAhydcxbtOmt0rteKc8xWngAshmIPkkr2bLrPbo7vINHT01MzTbm+9Tf6ofRSNgkKYNs2yWSSjo4Ozpw5E/5eKRWek5WVFaanp0OFSbQVul7ut7a21i6k3c5opqEhQECYWmump6e5fv06/f393H///aEvwk7YbaS7vLzM6OgolmVx5syZll50uVyOmZkZlFJcuHCh6RRFPaItwNGus2ACcC6XI5fLMTc3R7FYRAiBm5pivPPjGLm5S++J/O/ASokuKtgSogGzHoC5XJpz5r0MpV+372NvFo3STdu1/m7l3ftSMwnyPG9TZG9ZVsPIeCu53/r6Or//+7/P4uIiTzzxBLFYjBMnTuxr51QqlXj00Ucpl8t4nscP//AP88EPfrDp9QK0SbcB9tPQEEAIweTkJIuLiwwODvLAAw/sWWy+3fMZY1heXubatWvEYjHuuuuuTVX3/SCXy4UTJnp6ekin0/tKUewWQSTzpdWPM23/T6yuCrILRPW9kxZ5UlIj605N8NG0EX1H1wRDnudzH+L53IdQGu5N/R5Hu87d9L8nir2O69nJl6LeJOjFHkvfLPZSSNtKYVIqlUilUrzvfe/jG9/4Bn/6p39Kb28vf/AHf9D0ccXjcb74xS/S0dGB67o88sgjfP/3fz8PPfRQ02tCm3RrEMi+ggi12YaG69evs7S0xOHDhxs6cu33GBcXFxkdHSWZTG6bE24GuVyOkZERKpVKmLOdmpraV0Fvt2Tz7cW/58nSR7CqPJGqnrYYBTJSIQRstYwEjCmSwcUSoAW4xibrJlgzDhmpSVsF0kIhbZit/CgzC2CAioE14tgc4krmt+hODjf9t26HVkijtvOlCHYI0UnIAIlEokZJsdud1kGhUaS7VyQSCR5++GGMMXz4wx9u2W4seG+5rhval+4XbdKlVmP71FNPcd999+2ZbKO+CENDQxw9epSBgYGWEW5gTD46OkpHRweXLl3aMSe8FwRKB9d1OX36dM1Wdz/NEcFjtzuXi+tzfG71J7Ek2NZGlGrpEj2OW13H/73W4Gk/ehWiTKcpk4yc4ujTJPDI2LlwvY1ETW30ZwNJUwau8631H2Z1vYvLnR9jIHFvU3/zVriZkrEoGW81lj5QlQRj6Z977rld5c5vNhoVnJvBzcjdK6W4fPkyIyMjvOc97+HBBx/c95ova9JtNA6nUqns6c0RbaWN+iKMjIzsKzqMHuPs7CxjY2N0dnZy7733tnTky3ZkG6AVpNsIi9lpPrnyv9IVVyHZehoqGo4669iRl2C54nBDHeKUc4O4LNMvi9hyg4wBtAFVZVZp/JSEEBDQiFW9r9Fgim8gG/8scWGICcJ0RdIyJEyWF9bexTfXz/PqgY839Xc3wovRkdbIJMh1Xb7xjW9w6tQp8vk86+vrzM7OUiqVAMICVUDGN9uXwvO8lkTfwXXWymO1LIunn36abDbLm970Jp577jnuvvvufa35siTd7TS2tm3jed6OOaZCocDY2Bhra2ucPHlyUyvtfjW2xphwosDKykqNMXmz60WPbzdkG2A3Bb2d8s8BlktT/I+F/wNlPNLk6KoW8rWGrLbIiBXuiBufBCNc3euU6bG/Q59wkdYG2SoNU5UkOfqoKOi31zji5BESbFOiU2gsUUu+WGCcP6HPwGw5wZLowqFAWips/PulpcLwHE8uvgKBBDKctD5GT/e9+3pT3wodaUFkGaQbBgcHw9sCk6DAFGdmZuammwS1qjkiyOveDHR3d/PqV7+axx57rE26e0Ejsq2/aALSjcp6ogh8EYrFIqdOndqylbZZ0o1Ogejt7SWdTnPXXXfta1sakKZlWayvrzMyMoLneZw5c2aTC1Uj7DfSBf/v+oPpd+KaAkJAp1zHkRtEGBMFDseD1yVaCCvTJytYgvC2jcX9f6fiRWAKrf2UQ4+1EbkG968/fAE4Eo4nSwyZkh8BNzh+gx99L1NiXL2L0UVg8hdIJ++uyZO+GL4DzWK7brTAJKi+TrBbk6Bop9lu0YqcLviSyVa6iy0sLISa62KxyBe+8AXe97737XvdlwXp7rahATZItx5Bw4HneeEI8510uvWjS3Y6xkBaNjAwEE6B+NrXvrbvlk0pJaurq0xMTIQNE7sh2+jjt4p0Pc9jfHycmZmZcKRMdGsqhGCi+Cx/nf1I9REleqUbFsWMqdBnlWrUCEpDTiuGnQJSREhYbKQQDJoMFVJ2Nb8LEARLQRQMrLo2s5UkyvYjoJISFLXNK5PTfgSMv76U1dRDHfkKIGbBIeNhyJHVNpUTH0FxCbvy4YZeDPXtv7camrmetjMJiuppmzEJalWk2+putJmZGX78x388LKz/yI/8CD/4gz+473Vf0qS714YG8C+uKOkGGlgpJcPDw7t+US3LCnNk2yE6cufQoUObpkBs1dm2W6ytrYWKhDvvvHNPZBugUaQbjPKZnZ1laGiIK1euoLUOq+dBamRM/yMLK/+AAeImT8bRIYlWvDJH4+WNIhlQ0ZJOkedwzAt/71Wfuqhs8m6M88nlMGUAm6NYjB/pSqDX9uiy1rlWktiW4q54PkL4PsG62v8+bm2Q74LrMFIZ5Exsmm4JlvTX7JEe2uRYME8wF/tfuO/OL4RPG23/jU5xiMVidHR0UKlUWFtbe9GbHFrpdbCVnnYrQ5xGJkGtinRb7TB2zz338NRTT7VsvQAvOdLdT0MDbES6gVIgHo83pYHdKmIOEEzWnZqa2nLYIzSfplhbW2NkZAStNel0mrvvvrvpfFfUZSwYFTQzM8PQ0BAPPfQQUsrQYaqnpyck9i8t/r/MF/4BgJTOkYr57KgNrJTL3JWuJdyy1hyz17GrL5XSYKqvm/I0J+M5RCRw1NqXhlWPkvFSgq+UzvB9Hc/RLQzxKlHaEs4mVwGIAwkhkUIgaHRNGDSGtOPSaU3zd8VznI9dJwZkTIGE5RPwIeNS1JM8vfRK7ut7HGjc/hvtOJudnd3U5BCNig9KV3sQo3q286UIHNsCk6C1tTWee+45MpnMnme9RZHNZm/5FmB4CZFuKxoaAnnNzMwMvb293H333U07ZW1FltH5ZseOHeOhhx7a9lN+r/4LAdkaY8JW4KeffnpfbmVSSpRSjI2N1Rx3EC0FUXBQrFsrZfkvM/+CTsv/vaOLG4QLFFQlJFxtACHxtMsd1XRCEN1KKcBUGJAKq5piN8BqBTIx4TdMGFh2JUkR476kxytT30aEl3XwnCZMHjdqnqiFwJKCNJCShh/LfIcVT1MSfgJZCxOmJJIWxMwq11fOIYRFf/zjJJP3164W6TiLxWKcP38+PFf1wzej896iZNxqKVcrI929opFV5Ne//nUuXbq05ay36LnYziTodjC7gZcA6QaRRLlcDh3t90q2WmtmZmaYmJjAtm2OHTvG6dOn93Vc9aTrui4TExPMzc3ter5Zo3W2wurqKteuXash2wDbWTPuhMCQfW5ujuHh4YYTgYPzbYzhy0t/wT+sfZKMNAihwZTpivkRvxFgtOBErOQTrtZYEgRlBu0KUhhMlR9jEuKU6IgoFZYrkLIFmWpjX0xDRtocioMQ1eeo+zMNWxUCDRUt+Oz6IJ7VgTZwd2ISS7gMGOi1NVIIpBD02hJDlbwj0bE2hpyBIgZQLJZ/DMq9HO/+xx3P61ZNDtF5b4GUK6oeiObMmy3eHUSkuxcYY3Y16y1IWSmlNpkEOY5zW5jdwG1MulGNbTabZXp6mosXL+6ZbKemppicnKS/v5/Lly+zvLwcDjXcDwKyjDZNHD9+fE/zzWDnSHd1dZWRkRFga5ObZtzKornm3t5e+vv7OXXq1LaP+cLCJ3gy/zkS5P0CGNBjedWIVoOBHplDCo0whDrbHnzCFWx4JqRNmXj1e1El4v6Y/0PKQEL6l25AyGUNc65NQrr02bUEW1aGSSVRQA/QXyXUuAVv7pqnrOf5q8JxZgopvqdzxdfzEqQujN/tVke4AJYUdApDJ+Bqw7IRCLHMZPYcg8lnNilgdnNtRtUAUSnXdlMtopFgQEDb4VYj3e3QiIzrTYKmp6d5//vfz8jICH19fUxMTHDx4kXe9ra3NV3InJyc5B3veAezs7NIKfnpn/5p3vve97bkb7rtSLdRQ0MsFsPzvF0Trud5TE5OMj09zeHDh2uKVzvlYncLpRSrq6s88cQTex4mGcVWkW6UbHcyudkL6UY/iII25mKxyNjY2LaPe2H1SZ7Mfw4ok7B9MnRMqUqKGluAEBXickM3awxkTG5DfYCfbujEJ9xgGx8EqTbQKezwzgYoaJeSgYSBYzFTUyBbcAUzJsmKF6PfyuNIWAE+nT3BK5PjXIhVcKQkYQne1HEdhamJiivKMO8JErYmLSAhgqgZRirQb0F3NeiPWZJDxlBUsC5gvngvg2wm3max3VSLIEURNVKvV5JE7SJfzPRCK9DIJOhTn/oUv/zLv8zly5cZHBzk+eef39ffaNs2v/mbv8n999/P+vo6ly9f5jWveQ0XLlzY9/HfdqSrlAp7oAMScxxnV/KswBdhdnaWo0ePNsyn7pd0gw615eVlhBBhoalZ1BNmYN8ohNi1o9huSDcqWatXUexmXM+fL/8nMNBlb6gOOiyFYMOcpkuUQr2tNmBUjnhsg1hdAwOi7Ot3g9yuAi2hX1hY1bxq+HcJSAuHVJCTMGC0Ia8U10wcY2LMehJIMuUlOWovkrDguzsm+OvVYYZ7voOkGmEDFgKjDd/2wC+5waxKgYK748sUjKFbCKSAMzEYL0merBzibmeaQcsn7KQliBvDooH54is4Hn9+x9dnP9iqYBVEgrlcbpOHsVKKZDJJLpd70U1xWhl1r62tceedd/Lwww/zutftz03uyJEjHDlyBIBMJsP58+eZnp5+eZJuMO8+ip2IMjrC/Pjx4w3zkgGaVQsE0eDq6ionT/qTdb/61a/u+4IKjqcZsg2wHekGzRgTExMMDg5ukqwFj98uJ/zrL7wbAId8SJb5Uone9AbhLhagP0M1lwseNkecWsIdFOWw9XfdtcgDx2xNp7Q2FcCUAT/TanwVgtk41k4puc8o8rrEmk6x6CqmygPkEgkuiDGGHXhP/3cifg6aNWPolhZSSi44cLXs4lkOA1aBBdXF1XI3Z+NZskaTMpKEhJMJTaG0xHPuAPeKBdIYkpafdhg0hgWtmc6+gWPdnzlwT9+ttuXlcjk0uL9+/fomD+MgMj4oH4ZWj+q5GTnd8fFxnnrqqZb4LsBtSLqNsBUpFItFxsfHWVlZ4cSJE7saYb7XSLdQKDA6Osr6+jrDw8OcP3++pRdruVzmxo0bpFKphiN9doNGpBsUD8fHxxkYGNjWenK7jrRvrzxFyeQBSFe52lVwLOWGUW3Bg+HU2kb0ik2/yNYQbj8bhDtfdkg6mnMWCBEoJGBFwfMqyUWrQKdF2ASRVYoikFWGIQsyVZLusOB7EgVWnAL/pTzEa5LfihTmDBUFywakMcQkLGlFBpuYhDvjMOW6rAmHpMlTkmlmKwkOx0oU0FQ8QacNFxJlzlHeOE/+0SKEYEAaNFfJr57Dsv/yRe9Yi9pFdnZ2hgW8QLUTaKzn5+fDukZ9t1krR9JD67rRwI90m9Ghb4dcLsdb3vIWPvaxj7Ws2+0lQbr1yOfzjI2Nsb6+vqsR41HslnSDsTWFQoHh4eE9F/F2wsrKSuhnOzg4yNmzZ5teK0q6gafD+Pg4fX19u/L53SpSLpVK/PHsbyAlWKoQ/r5gZDWH60e1fzV/iX95x5cB8ABlSsSqV17FQFqUiVUJdyyf4FCqzJBlhfnTvAdfcgcoV+AHuuarsjJDSWmyxjfIKYkYloRJBSXX4UQsR081JdFrC36p/ym8qjrCVYZlfMIuKFgXnfR4edK2YhWFo2y6LcOQA9OuCw4UPY8FuhjQZXotg7ANRgiSOCSr4mFDQOYeRemnwIQxWBKM+gFO3XEM+JumX8dWoX5LHzXFqb9fULyLtv7Wd5sFSopmcCtHuq7r8pa3vIW3v/3tvPnNb27Zurcd6W5HbOvr64yOjlIqlZomwp1IN3iOcrnM8PAwfX19N4VsLcvi7Nmz5HI5yuXNkxL2gkBnOzMzw9jYGL29vVy+fHnXRZ6tIt1fGf1J0tInmo5qEctT0G/lNrrODLx1cENCZYzmmFUKCTkmynRLQPgeB0dTJY5afiXOGEPWVXzD6+Zee4XBpJ+319WtuzGw5Dl4Mo7RFVLSV0TEhMdnlk/yz3rHiQuJU7VytI2gqDRL+MeWM4CIUVIu39IDXJQLpCxNGY/JSorjsTLHHDCeS9xa5ngsKO75x6Z0rRudqJ7rpIyRMD4Rl3FZUy5SQioxTWn9PhKZp5t5GVuG3RLdVj4M9SPpJyYmwinI9W3gO0WxrYx0XddtWeHSGMO73vUuzp8/z8/93M+1ZM0Atx3pNkI2m6VQKPDtb397R8esnbCqFwSiAAAgAElEQVRVqmJtbY1r167hed6+n6MRVlZWGBkZwbbtmg64YCpAszDGhOL7Q4cO7YlsAzSKdD/8/L8nZkBEmgWMMTgC0tZGTjyvNcfigapAkwh9DnxFQyaIiA3YAo5YPlkpY8gbD9uGR5y1kISLSrGGQBuY1zGEjKENzLhdJCQMOav0C5e3908isPD1uAptDAlpk7Qkh43haTdJHA9bQI9lyJoY19x+zst5bAFalphykww5JYZsg8tG/rnsGQrVz9nvlBXnnAqWhLIC29h023a16AcJGSMuHFbdCp7lASVK65dJZJ5s8hXdP1oxfr1R66/rumGKYnZ2llwut2nWW6CkCEi/VZFuq3PmX/7yl/n4xz/OpUuXuO+++wD40Ic+xOtf//p9r33bkW50GGLgixB8wt5zzz0t+6QLEDQdaK33bBQTVP23u8CjZHvu3LlNkqBmC3vB1OFgdtrQ0FDN4L+9IBrpaq35w+f+nGXrGh2WqWoJfI8JC9/D1gpTC5p+m5BwLcCmABj/vgKcag1MAj3S3jhnQtApnNpdhABHSDqBkhFkK1A2kNcCS8boMHlOWF6YSzbAU+tphlPrSAlF7dEtbRwpuT9W5nP5Du6IFbEEHJWL/D+LV5B9LhdiK6SkpkMU0IAUAscI1pRh1QjWPEmvrXAk9EqXFzyH844iJkGJNzLq/al/uNrmdDyBALqcGMo4ZE0RyLM491aKlY/ue3veDG6WTtdxnJo2cNjQ1AZkPDU1FbZBB77QQgjy+TzJZHLfx9WqXecjjzxy04qftx3pRicoJJPJkKieeeaZlm4vgm2+lJLTp083ZaQREGajCymYb+Y4TkOyDRCkBnaL4Pxcu3aNrq4uXvGKV7C8vLyvFEVAhFNTU3xp4mt8Jf0Z0taGjjamNUL44thAvQW6Sn4VYMMArENsRMcJ4yKkT7i2AEtsnjKhtfa38vi3y+ptSQn3xDXKlHm+InjAmakh2yXP4AHHUnnKSOJGg4AF7dEvHSwBr03n+KtCN6eddRJS82B6kpFCH/c6WdKWH50rYxC6utWWMFoBlw4WSi4XkwViEoQHkxXJ8ZjGMn/GmpZ0So2WHn9fFlyxS6StOLaU9JoUK6ZAOvkMrp5laak73J47jrPJpexmjGI/SJ1uVFNbP+utWCwyOTkZyiyD4l0zJuqlUqnlAdfNwm1HuuDnVevH1TiO05KmhuXlZfL5POPj45w9e3ZfFcuAdKMSrN2SbXSN3aQXgtlp165do6Ojg/vuuy+MJJrpSIuuG2wVF5aX+evk5+ixfd/aYLsdszYG4JRM1TDchAsAOszhSsuf6mDwmw2Cx2WEXRNRu8aXg60pTbG6TLGawugS0CstYlJiS7gU12gEFWP46voR8qbM3eml0MdhvuSgHckJ249q55RHv/RVCq9PZfmHYpJBu8L3pCfocwgfp7Rh1UDZCIaEwRKCO214vAKuyLDkavpjJfosl0WdoqRLJC044VisKI0ADpl1CgKkFiRlzG8rJsUyBbo6foIjR74Znuuotjba8ho1xgnkXPuJCG+FjrSgeJdMJuns7Aw1sfUm6o18e4NzEW2DbrXD2M3EbUe6QghOnz69KfS3bXtP/rVRGGNYWloKXcXS6TSXLl3ad5QRTQ0sLS2Fk3vPnz+/62GSO6UXgmO/du0a6XSae+65Z1MVuhnSjaYnuru7SaVS/Cf3d0lXSSkRcZ2tDUIsotPIYsINidXV1NxWkzmIEK5X1d/Oa9/tYKzQQTJeoUP65L1mNBlAG4WDwEJiCUkCOJu4gVt9QqWhCKTjsK7hnwrd3J/K4gjDlNIcMRZJy/DdySIVVNh24WlY0IJENSIvK8PTXoJXxCukJdxBngmRYET10GtmEIBtKox6cc7LMilh6LVTflTubN7yiirxKqNYX3sNsc7P++dqC21t1BhnYWGhJiKMkvFu5Vxb7b5eDASNGgG2M1GPTkCOWkU+++yzjI6OIoRgZWVl37Kxd77znXz2s59lcHCQ5557bl9rNcJtR7rQuJreTCdZsBUfGxsjlUpx8eJF0uk0Tz75ZEuqqlJKlpeX+eY3v0k8Hm9qcu9WhBnktEdGRkgmk9sOqtwL6UYj5kwmE44JevtX3kNCuvjpAhOai+u6tFfJXa8SjR/ZBtIxALf6VQgoeW51uCSktAxfUw/fmnNO+YR7Qzkk4wAxXFXmdExFjM0FBaUoC0WvsJFCcDRmUVCKr5c6iNuQEWWkhLTQzMpOnivZXEoskpGKp8qSh5J+BBuvRsqrHhSrB5xXhrQFXbbh2/k4HUZwZ6LM6YThesklrh0my2mGYnmOWh5H7BiW2GgqiObBNSCqqZPgnxSSTrOAzt9PJf44NLjetjPGicq5pqenN3kxBF/rm12CKSK3Anb7PtvKtzewcL169SozMzO8+c1vZnl5md/+7d/m0UcfbeqYfuInfoKf+Zmf4R3veEdTj98JtyXpNsJuW4GhNorLZDKbCKveyHyvCAhxaWmJcrm8rzHpjSLdgGwTicSu7Cd3S7pLS0uMjIyQSqVqIuZ//uS/JpGARJUkZYRI62fcdFgWvhq3+txBey6GjNTVqNeQibiHxaJRlzHYRnJUSuY9iYVgkCIDsahKwtfXzmj41MxF3n30myxLj5QWpCyblGXzUKrEM26CdSPpNH5n3DF7kWvlAZ7MWzySnuP+pIcrBDEjkUIijKZUtYvMejZlkcQxq8Ql3Jda5enycU6YaeISviu+ypdL/Qw7MU5UOS344NDal50texUKAjzTh8er6eXPSNigFAzYaYQRGAwWhkT5AUre30Nid77N28m5gqg4apEYVRAopTblzl8s7De/7DgOV65cYWVlhVgsxkc/+tHQU7tZPProo4yPjzf9+J1wW5Juo4vFtu0d3cGC/OT4+DhdXV01ec/6tZpVDAQ523g8zsDAAIcOHWqacKGWMAOlg+M4eyLynUg3WDcWi3Hx4sWadf/bC3+CMX4xTAiway5mHTpyBT932jrsyvJdxCQQkG11HQiLXv66/g8SSVxuvAGPWXDUGIyRFFEYYK4i+b+Xz/PGnufolfD2o9+iJCRxrckLQ1m59FgxEhIuOyX+It/HktCcimVJScWV+CxHYi5CyKrWVvNEyeJK0hATgm5lmDagEXgKrplezjvLxIUhrXL8YfYUP9kzRq+I8UPpjbZngJJXYZUKWQUDtiRpSfJKY7FEl3wfq/rPiBv/g6jAu0nxfyGRaKP9bjH1KCX2N6mgkRdDvYKgUqnwT//0TzXtvy/WGPZWTo0IGiP2au160LgtSbcRAj/NRoj65fb09Ow4WXevqYpoXjWZTIaEeO3atX2PYbcsi3K5zBNPPLGlrGwnbEW6a2trXL16FSHEpnWffXqCX/zZ36XrowXf8at6DUfH5PhrB99t+CwEJogi/NmXNEgjwttkg0AkTDFohRFgC4nANzZKGYFnBDMK3tr9PKlIpOxqeLxwhO9Kz/gKBVWh34phS3hjeonPFvroxKPX3nAhK2soGoMl4HjMZV5ZDFqCTluy5AISrntxOgUseDEG7Ar3ppd5nawQl4lIVGtYNQUq1c64ZQ2rxqLL+JMrOgWsGSjpJ9BeN57I4lhQ0n9JQqaxyIP2z5EQglj+ESrpL+3p9d0J9QqCubk5rly5gjEmLFrVe/dGifhmStpapaS4XaZGwG1KultFuvVEGTVzCfxydyMr2S3p1pNt/VZ/v2PYs9ksIyMjlEol7r333j2TbYB60g0mAmut+fD7v0i5rIHP88nP/2tisRj/+SOf5XN/8SzJj/mtvdForoZwqf6+mkKoBrVVB9qNv1vjl9eqN9cgWE/gk1jBeIjqnLKiqlAW0IFFUto4Er4raSgomDOQ9yRlNHHL4UJqmUmV5LhVJCZ94u2TAfEuUsJvAdYGJsswrzoYiq+Rxv8gmXdBYDFgaU468IKnuMNeZ97L4Igkpyw/ZxzkkpVW5KngaUXFwJKuGq8LkCbOt1y4N1YiJgVCG5T+VbzCW8l2/S79BoQYoeQ8S6ryEFJIlnWebiuFJA+lz0Ji/wMQd4IQYkvv3iBFEe04i0raAkLeL2G2KtJdW1trWod+0LgtSbcRojndqCfs4ODgrvwFotiJdHci2wDNkm7glSuE4OzZs+H8qGYRkG4+n2dkZIRKpcKZM2d4xw/9Xs393vKa3+JXf+uH+dxfPIv9G0s4kai2/nPOhLV+F4PEVEnWiCrpGuUb0iDxxVPV26sk3SjjJhGksVmpaCYMxKUhBnhCkUfRiU1MWqQsi+PG5vPlGN02ONovlmWk5vHyEV4VnyEt/RZcaXyyTBiLv8mlOBzPkbJhyM4DFkWlSFow6CiuVpJ0iDIpC05b4ErNq5IFJITpiKLnsoYvl4tXVRPzSjCvHe4QpWq+28UzDhUDcSnowJATswgjQ69fgJL6JFIbkpakR6RwAQdIqA9Q4uaT7lbYqmi1nV1ks5K2VsnXbpdRPfASIt1AMjY+Ps709HTDybq7xVZkGVT2R0dHSaVSOxaxgtTAbhGdb7ZX+8btUCqV+MC/+RsCs5XPfOHnmRxfDG/v7LZZy/ofMr/ys/+Dyt1Z4jGQYoMajXEJBtZANbg11ZZe4VOqZSJbfjQJ49+mjUTiy5Q0EmMUCIGutg37RZ2NXFxvTNADzCtFHigZSVlrkB6W8egRcRwJr0tXeL4syGrBCVGh04KTqTlAhm3Jee2SFDFsKXh1R4HrauMYF13wRDd3yCwxCcftPH9X6udNqTxpa0NRYRBMVCAu8mEKRWsooklZFkcsuOElyOGRsT26hYvSDiOuw4WYS1wK1pXBJD4JQFkbDllxpPp1kNXUghQ41fkUABS+CKnvbcnrH8V+CkzbSdqCfHEjSVtAyI0kba0q6LVyVM+P/uiP8rd/+7csLi4yNDTEBz/4Qd71rne1ZG24TUm3/kXyPI+pqSlWVlbo6+vbcrLublFflKsn2+3kWfXr5PP5He8XbPeVUpw5c6ZlF0+pVGJ0dJTf+D9rna3+95/670yPz4c//+Gnfx6AN3zPr4MokfpJVSXBjcf40V5AWJ5vUGMg2D8YI/HnlPlJBA2I6v9+tFv1rkWghEHjE7VjQAufbGxhhQlh13gcsmK42vAd7SKERcFTDNiQpUyn8RsNLsYNJeOigvxxNYXw/y0d54HuCWwJa16FI46DIwVDCMZdgyscUhasaXi22MFrUmXuiMU5FysihBWqEMYqFVaUxfHIqPiygawGZSyGJcQESKMZVd3cYy3iVBUeFRyUcbGlICkMOHMM24nwwyUgQKUVtrRDyhVA3Pw85X0W1Rqh1RrdqKRtYGAg/H1U0pbNZhuOF9pPgbkerWyO+MQnPtGSdbbCbUm6ATzPY2JigtnZWY4dO0Y6nWZ4eHjf6wbphahmNWiY2MsY853SC7lcjpGREVzX5cyZMy3zAn396/4jAL/0gQc5ffo0i3O10fb42CLC89/g7/ipVwHw3377C+Bp9G8UsWrekx7gR6QBB8swM+vnSY0RgIes0oZ/awxMCS0BIykDMSHRQuMZ0MLXxToCSigs4z86GIvuCP/StC3DRSlYp8JG/xoU8YghcbBICIeSVny75GA7/rN3p2e54fqtuY6Eac9lyHaICcmABVMKhi3oiblYSYkQiZAEZysWeSocsl0GHEOH7Z8/pWHd+CPYbWBZx1jwKgzYcCWZ52+KaVzj53bTFHFJMuU5nHBcjlgJhL1Btkop1kSFhZXz9HQ/z6BJI5Bczf4aZ7rfj8BA/t9A+qPNXwgNcFAa3e0kbUFUPDc3R6FQ4PHHHw8lbdEW6L18OKyurrbcS/dm4bYkXaUUIyMjzM3NMTQ0FE6CuHHjRkvWtyyL9fV1vva1r4VGOnsh2+g6jUg3l8uFXrlnzpxpmWOZ67qMjo6GP3/qU9f5D//hVRvHY4OqpqqD+tc/++ePAPDnf/x1yndmQ6WCMf5xW9TpcjEI4QAuFoR5X2XiGF1A2H7USxgVS4zQeCZOwZRIIqkIG9solICk0VQEpIEKHtIInEhLsBQSB0EvNsYYXKNZM2VcrSlKTQJNh4yRkDblYg9lsUivrThW7QQrGEgJv8CV9yx6bYujtuCoXdtkU1CQVQUqUrOqBD0WGGP53WbAgoIKFgZFukqs18s9xOUy/dVzcVgUuO6lOe3kSVmGdQ3DtqFDxsPncpViTFfotE4QZ4K+7nE/6jcaW1ocSX0yfL1s9XmefOaZTbPObucWYNu26e7upru7G611OH8sKmlbXl7eNNEiIPCtfBjapHsASCQSTQ973ApBh9rVq1dRSnHlypWmyDZAPekGxuelUonTp0/XGIBsh53cyqIR/4kTJ8Lff/MbMwCk0jEK+QpKwS/8yg/wkV/9Sz9cU5qv/MMLfOgDnwZbYr1rc8FMRn6uzwYqseFT4HdYBcfrj9NZ0wkSokRcVqNYYyNROEg8SlhCsm40XQbyxiUtHYSwcI2Hhd+wUHMsUhJH0m+s0EEn0ARj4Lv7VjE4GJzQ9zb659R3ipWVZtpVzANfKpyjW0zzAx0LHHZMlbAVKWERExYVfGnakunmKKskLcN98TmWVJxZD47YcDGe4/FyBwboRnImbhDCDvXAWV1iTgtsAQXdWR2+6aexVk2RXtNBJv40Lhex+SZSGs6dO0culwuVBIWCryiJ5kv3Mo79VhpK6XkeVtWsvpEpznbj6AMitiw/FVQulxtq7m9F3Jaka9s2Q0NDm36/GyvFRog6l3V0dHD+/HkmJib2RbiwQbqFQiGcMhGQ7V6KB4H6oP7vUkrxL//Vf+f6dV+f/Gd/8i9IpTZfeK98+DR/+4VvgYFHvuc8ir+sapscfvUjfwlpG/fnlkhE3ovRC6OWbCWednGqvBcW0wzkkHTo6rTfahKiZOJorXAEaBPDlhWEMORMEssUEMJGaY+yhBR+UckAGI02GktU56MZP/eJAGE2xO8+uUbGCUePO1I0MoCnFXnlsmwsDtsaI2AeiAuXH8s8U9UgbxiUL2jBEQlJS3DcMjyrHEoGRr1uLsgVUlLxvJskITWH8McTHRclhu1E6IamjWFBFYkLUXVJC44tBxFpXVRg5/KLWObHkGxNRo3GsQcWp9tJul7sSDcKpdS2tZftxtEHKornn3+eX/u1X2NycpLXvOY1XLp0iTe84Q187/c2X4R87LHHeO9734tSine/+938wi/8QtNrNcJtSbpbEVYgG9utxVvUBrGzs5N7772XZDJJpVJpiWOZ67pks1meffZZTp8+TX9/f1OV2oC8gwtUa83k5CSTk5Mh4QK8+a2/y2Of+Vn+6rGfC/O6v/KBP+X9v/gDPukCH/mNz0IyCBN9UipSIBarTStslogB+JInpzrpwa/sB5GuAZMgbwAUjigTFxJL+EoGgY0RkpxWdAkPISRr2qZbeqwi6UWzZly6cPwuLQwSgWdUGPVa1QKXJ/w21knXpQAkUHRa4GCRkjYgyGm45hpsyy+sfafSx8PxG9gSlHHJqxgdNlyKGXK6tr1YVyP4JIYpBacl2EJwSHrMViRJy6aoJSlLc8FZYV3HmPHgTsfhXHIjb5tTZda1ZtlYHLE0dlU+5jsKVzD0I8R8+Dq4xiYmPGJ8bOPE52ah43DNa7FVvjRqIh6VdEW36LdK+y9sRLp7RVTSdvToUb7v+76PRx99lE984hN84xvf2FewpJTiPe95D5///OcZGhrigQce4I1vfGNLpgAHuC1JF7Y3vdmJdKPeC52dnZvagfc7hr1YLDI6Osra2hq2bfPggw/u60IPSDc6Jv3w4cM89NBDwBM1933dGz7GY5/5WT9yigm+9swkP/SjvwcOSA1/9/dX/TtqDRW/PVf/UjEk2eCCqD9c14A0HiZMKfiJgoIyWMIQE1QLcAYLC484Ag/LSFxjI0XRD4lFgqyq0GNVQCRY13kyUpJVmm4LykaRsATCiJB4NQatFZbwZVy2sdBojjsOec/jhhGseoaUrShpRa+M0yEFd9gW15XClnA2NsO1cpK7E0Xi0o8nTVXDm0Sy4vnz1oSwwHhkqqN/VtwEz5XgUqLMUQdu6CI5k2HM6+aCXCYlNcctSacVC6/JivJYMRXK2ldHIGBZCQaEIS2hoP3z5Jl5X+mA/6G2ajroN1lsnog0keR2fZ1sZSIebNFzuVyYL83lcpsaHQ7SSB12jnR3i6BxY2BgYF8RLsDjjz/OmTNnwoL82972Nj796U+3SXcr7OSpuxPZBthp5PhWCCRaq6urnD59mvPnz/PVr35135GFEILZ2VlmZmZYL8J//M8RorUFwvOHH2rtv3lf+8aPUZMrAIjb/htZGQ73p1mcWgUhyL9pnbgFwQY3eqi6qrMFjSOFb9gtJEoLkH6RSwv/jVrUgFY4UmELjTBVl3OpEWjyJkXaFKptawmyStFtKSoiTVnnQdoUlQeWxjHaJ9gq8QrjT29wjcKORL3SSDpswSkjWTcuQTy+qst0yTg9liIjNPmqAc+JVAkIPqwNK8qjx3KwpcQoEaZJStpBoMhYmjvsHP9YOsWdZoqkhHOOyz+WXGLS4pCIkbKsMEpWWlEwLgaYcSFj+1mcMlAsnwR7LJIFyROoP6oNfXSySJCsNphqTnsRaL7TKjp0cmBggEwmw9raGidOnAiJeHFxkfHxcVzXJRaLbXIpu1npiFb6LrRqUu/09DTHjx8Pfx4aGuJrX/taS9YOcNuS7laRbiOnMWMMc3NzjI2NbWt00ywC5/tsNsupU6daNoY9OO6lpSWEEDzwwAP84Fv/6+b7Wf5AyF/5xe/n33/kf27slT3tc6mkmhMALMHsch4bMNpgtvgA93O1ojpKx29eCCQPFRMDrZG4WEgcCUJoLCzfM1dKpICitjBUcIT/uHWdIk2BmDAg0mR1iW7psi7SaF2spjg9wKNH2lXnL9+Jy9UaW0iU0XjG/14isaWFZSSOkSzrMvXEa0tJStvk8fyGBg1rWtBj+8S74FYYcBwGHcFKxeBik9MWV90EDyeXsAQcEws8ljvED2XmyEj4/lSBhLTCIpnWmgXtjyxKVf0ibGFhjK93Vhpi8etVOZmgE4MQ2Wr22i/2dYu4/7IFZFu9tKX8IpqH9n0tBQhyuo26zuqNcSYnJ0MVQXSaQ6tGsbeqqNfKbrRGwVYr3stCiDcB/xZuY9JthPpIt55sdzK62SvK5TJjY2MsLy8zPDy8p1Hv26F+5M7g4CBDQ0M1279/+77X8u9+868J+0oFfPDXHyMY+40GIX2FggBQPjngWCAlXspm/U1ZHKp2jdV1qzWr6nveq1EAlAzEVdXkxkg0Dgh/woIUCitYRykcy4+KSzqBESViUiOFJG9SKFMgKQ2GBCtK0m2VyYskqKIf6IXE6/hyNeOP6akYj1jVN9czCikMNpZPzkLQLyTKGGbdItKCrC7TLeM4UmIrixvVPbtlgTIetgQt/RlrHRbcGRM8X4EOS1GsJLha6OJsapXhRI4ez5AQgf+CqEa2mkVdCncHawqKCPodOGbDoq7K7XAwokJWw4AwJGTgSgwpbL8phOpORT5CjH8M9cqWntjkV7EfbEd0O6kIGo1iry/c7SVybWWk26rGiKGhISYnJ8Ofp6amOHr06L7XNcb8OfDncBuT7lamN67rtoxstyo6VCoVxsbGWFpa4uTJk9x1110tI9voFIggIn/hhRf4Vx/4M2bn8+Er9u9+86/9b4I2sar7jPHw3+m66o5gV+eOGYPUAl3wIGXjSQ13gBC1aQW/KKZBqNCAvGz8Fl9wqAj/TsKAEB7SgKxGuYHrTUU4KFUhXp30UDQJyqpExtJIJEU6cFWRTkthRIxlZdEjC+REEleVwRIo7QJuDfE60sYzGmkEspomqGgXS1rV6REWUhiOxfxCijKagq6QkA4dlsUhFHPaT5ssa5s+4Q+xXNEVEsI3xzlpaSa05GJiiX4RY8hKIoWoeubKsLFhxZRDxYHWUAAWtU1KGvqMb95eMpAEkrKEiyRvIOjZMhh62YgWPaPIMUyaCo4xyLDVZH3f11UUzagXoiqCQ4cOhb8PGh1yuVzN9N/dejEopVqSR26lw9gDDzzA1atXGRsb49ixY/zxH/8xf/RHf9SStQPctqTbCLZts7i4yOTkJN3d3fuKbOsVA+CT7fj4OIuLi5w8eZKzZ8/uimx3UzFeWVnh6tWrDY3Jf+YDn2v8oGBJvy3M/4Vd9Zqx/SZcXIMUItSzWpZAlRXr/1ueeEC0wVf8bjI/+pUbSgYcisog8arNEKKqVrD9xxiNEgYpPL8lVhiMiKGMRwLPHzgpkqwqjw7pYgmNJ5KsaEWnKGIJiyWVoUMUwBIs6gp9ElzlYnDpFv5Ic2l8gxljDCXtkZC2P7jTaEw16sUQEpYtLTplMtwy9lmSbgtK2jcrL2tBpzT0CIEwEscI+mxBEOP5H0gyTCHMVFxiliIlLRLCIm8UWeV/xtkSeqWHi0VJ+wKRw5bFqvZTDEZDr9BoI7CkpI8Ne8g8FX9EkejG4atsjPGEwIO4Vaif2bcfRBsdAmw1XijILUcjY8/z9i3LhNY2Rti2ze/8zu/w2te+FqUU73znO7l48WJL1gYQQjgvCdINzMmDUeb333//vtMIgYIhaqQzPz/PyZMnOXPmzK6jhUBju9WWbnV1latXr2JZVkNj8t/6r4+F28vTx7oZm1ypRrcw0JdmYSkfPJH/1eATr66mGALyDVNVAoHBSm8EyCK80S96+QUlSRm1kU5AoLGpgB/lGuPft1poiwuBwaFsDLb2sC2DNDZFYeOIEjEBQtjktI2kRMZSCCxWTQeWcOmUZdZNipzn0mvBonHoECU8ynjSo0dYxKplfokkJi08NGjfdzeIem3ppxtkEJUaFX7gSSmr6W3JVp5tAUH7ul7N1bIkKYqkbINtgxEbnsCuspn3oNP2yFSldK6GWQ0nJSSEIYufE++2FDbgIEnihIS74Gni1e02flQAACAASURBVJpjh3mmgd64tY0MN7sNeKvxQoFePWoXGXy/srISEnEqldrz8a2trbWskAbw+te/nte//vUtW68O//62Jd3gop2dnWVsbIyenh7uuusulpaWWpK3tW2bUqnE1NQUc3NznDhxoqkOuCBirr+Q1tfXuXr1KsaYLacOv/mnfpfV9XI4ZuHa7KrfdwsgBHOrBX8fq42/x60KTKXBf4w0GxaMfkUMZSD7zkJ1i2PCtILGN68Rojq1V4IUDhVjECpoQvBC4xsAK+wYsykqiAsPKQWecHCVR0IaBAZjErhCERcujgRNglVliIsKCelhcFhRDglRIikFS8rBFi7GEuSUjWvKFIRHL4ou6WCEQRjhR70YysYjJqww6tUYLPyCllXtcFvWJWJCkhYOEoEyYAgI1KCNpmw0WkjSlt/UsKo9+pyNuiRAxTMsIbjDhh4blqoSOYM/u00KRcnEcI1HTAq6hUFV0zQ9xMOmCaUVWVxS6S+iSt9HCgcpq4oRjK+gwZDPw/T8GJlMZtfjyLfDi9UcYVkWmUymxqL0W9/6VpiuyOfzTE1Nkf//yXvvOLuqcv//vdbep0xPZtJ7SCGhhRYM6rWCBWyIIspVr1y7WPB7EdvlYkHFgliuYkNQ9IpSBGyAdBQJho4hdZLMTCYzk6nnzCm7rPX7Y629zz5TkklmUMLv4TVk5pR9djuf/ezP83k+z/AwSqlRhbt9bffg4CCLFi36p2zHZEII8RLgzEMWdKOmg+nTp8eZbcQtTTaiOVOPP/44S5YsmVS78chW4KTJzYoVK8blon70i3sZzHsW7eyD48nYpEx6wdjM1uZj0XsqAlt0JgRUbGxjoFkZ83Ack61afauwmV3UxqvRFniNcU2ow3jGmUcaHXpWs+tS1iCVT9oBgUNJOZR1SFr4pKWgrDN4KkVah2Qcj4LOUtKQoky9I+gJ0qREmZAUQwQUdYG89JnluKSEIPovLY2KoKwM+ILhSB1h5GVpkWK2cBlQBQo6oEa4IKCgyombdx3vZ61NG3EdggIapaBfQa0EJJSVbxoZpGCBhE5ttLd11sRc6TI5Bc3CISsdQhT1ouK/UFJlQikRLCadno8oOTEFJITxaJP2gKbqDqOOuqpx5CMLWFE77ETi2dQGHPG/tbW1Vf4jkbY4n8+Ty+Xo7OyMt3uktjiVSjEwMMAxxxzzL9yS/YcQYjrwU+Bthyzo1tbWjpoEMdmmhiAI2LVrF52dnaRSKVasWFFlV3cwkWwF3rp1K8VikRUrVuzX5OZXv3us8ofWrFgyjS1bB8AxoJrJCl568ir+dOdGzCMWZB0rN4rSsxEXi/4zc1Xjd7TWuFG2a6vtUlgqQUXDIE1OaFkN0CEpYTxrEZJSqBCEFlzT+IAXg2+KkgKHAFca4C7hUA7AER5ZCWUhKYcpHBGSFT5KZCgHGaRQ1AqJpzMIQopOhj5Vpl8VmSl9ZropI2ez/6Wkg9Im65dCVjdVIJgmzdjzYe2RES51Mo1QHsZDLNLEwUDoM81JkZIOg75R1Aogo3xcB6YJ6FCwxCoRyh7sUc2sSveSwfC50qpAanDjDE0pRa8q40pBCpBqDzJ/omljtiFIuhaDcF7IrFmzqtpgkwWszs7Oqs6zJBCPlR0+m9qAx1MvJLXFI7c74oq7u7vJ5/N87GMfw/M89uzZw/DwMGvWrJmSRobf/OY3XHzxxWzcuJH169dz4oknTnaR7wdmAd8/ZEE3k8mMOnkOZCJwMoIgoK2tjd27dzN//nzWrVtHa2vrhMeW7yu01mzZsoVyubzPVuAzP/hj+nJG6/njL5wRP65CIAWbOoYgW6EWyiH86f6njZnruLzf6M8xXK5RHyitceNhZQ5lpXEi1YL9nAjMHSyXK+xrNQZoZeR9ICmHAL4FdQO+vgIHD1e6BKGhKFLCsB9KZwjM/TSu8MhICIRx9nJESFoEBMJFKNtgobJkZUAPAR2iRHMwzGFuieluRCbY/6SRjkltANnIyyQOEle6NGoHT3loIaiVaYq+z27byjzTXpCGg4A612VRCrYGUbu0oN4arjthQF64NDiClSnFHHpNL4gw+6mJFK71jdBaU1Q+w4RIKQi1xkXTIDUCD+xrorNZkyi+Oq8YdQzHK2CNlx0mgdj3/WcN6B5o1j3W0M27776bd73rXZx88sns3LmThx9+mK997WuTXrejjjqKG264gfe9732TXhaA1vrLwJfhOaZemOio8SjCMGTXrl1VYBudBJPNmsvlMtu3b6evr49FixZx7LHHjstJvfFDP6I/V/G8ffdnb7DdXBDrtpLUgrK9pRArF5TrmGYIVyJTAhWGlXTW/jt4Wh43Kp4lkmJtwdTc4jr4OmI1tLUgN45ippErJB3bibmUFUBACqN/BTOmBh3iaIXrQECaUIHWASmRQskogw5xLM2hdQY/NMBnpgdr0sJ49Br2JCQXYl+vreWix2MioEZ4rK3p5LBUiJugHRBmOoVAoLRCoWIfh7RMo1D4OqQllSIbBHQDgxqagILQZLXGlZKljqY10CgkBR1SL6HBgVAEcSuxaw9PIyZDTo5jLxFQIDQtx9ohLZ3YBS0eFy5AWNWFpnJ4maDSYCLZ4Z49e+jv749bgOvq6mKu+J89BRimxlBdSkmhUOB1r3tdVSfZZGP16tVTtqyRcciC7lgnyERPmjAMaWtro6Ojg3nz5o05aeJgQTep4V26dGnc+bOvdesfGjHSR4jq2TcAJTNtgZQ2Zq6OrPC22NQ0ZUX2ECGgiagiXxPREEaLYMLBU6bRwCG5D83nSzCFKxH56jp4yhjjpGUkE3YIwQJr9LhDKJz4MVcoHOniayy4GvWCaRnWCUo6NHpbNJ6wgBNqq7kw66VtIQwt8LQgUA539B/BYDmLKPt8cemDrGoYsJSIvbW34BtqhdYaRzj2c0zxrd4FVwkGdRAzNQMqoEWkcIRggaMooGKKRdhdFDVsNFvONgm2PUGZWkeSlpJpIm2zcGEPiSbQAcMqICscsnZYj5moMRfoGf8kO4AYmR2Wy2VWrFgBUKWxLRaLozjT+vr6KWle2FdM1aieQ8VLFw5h0D2YiNy52tvbmTt37j7H+riue0DzzUbKylasWIGUkmKxuM/pEcnn7vr5eZz14SvpGbAsogZV1sTuK1DJXkdE2pV4XoDSWGI2SlfN6wcW5A3AWQCNMqmILtDY0ZI6oh9MNqmEAVMISUkdG5pLaTJirRRaaWpSlccDjC4VDFUhhUOgHbsc87jUAik1oaUTKhGpham63TY0QXTJcPBCQaglpVBSDDLkyikGh9N4w1nO3LwcpyR5y9InOGf1IyxpHrRZvFlSqDVKB0giysG0EgspqMFoj0NtMtQh7dMoUmSlQ6B1rPLIajc24DH7sQK2Q8qnRjq0pDJVGW1sG6k1/dqjTjpkhSRrLy5RJ6ArjgVurxykKYxIe55KpUaN2ElmxV1dXWzbtq2q2SH6+VdkxfuKQqFwUG39p5xyypiF90suuYTXv/71U7FqY8YhC7rjHfSxPHWT04HnzJkzoRlqE53kmyy+LVy4cJTSYeRytNa85D//F4DjV83npJWVFsO712+hZ9CYVCutjTwsVX37JQBdDGjOuuzNhdBgeIKSDqka35ugHxSgj7UmX0LZBNrBC0xHVqTDNftPgxaWm9QxjyuE8VXQOsSV2iobzO0dUlqaQSG0Mp3GFsxDjEwNQrTCZtMOSoCYABMU0R9aC0I0oZYEoYOvBSXfoRymGPZSlIopwkIK8imcsovw4Pon1vCnm46BJs1Lj9jMm054lCPm95B2tJ1abCgHVxs6IIVLqEMUmpR0SEXmM/ZOYZpVRiTPvXjOmVbxbm9001XPRWGMzH08zMUrrSUpu0xtC2+q/nFSw8/b/445yNhXIW0szjQ5eDIC48hIfGQL8IE0XUxmQOZYyzkYmuLPf/7zlKzDgcYhC7pQySCSEdEC6XS6ygrxQKcD749eiCiK9vb2qpFBI8NxnKrlfOKy38a/P/x0B1t3dMV/f+67txqQTM7I0dqg1rCirsahcUEjHXsH2YuGZnv4RtK9aHNkLTYMNReQtmCW9FgAF4VGx36yynoLREAJUYYb877SIdQQhlbtoAMy6Wh5BoBNBgxmbhrWHNy4cUXZNON855IPm00XKPtvqCV+KAlDBy90KQcOJU/ildPokgsFJwZcx4faXQpda9bpzo2ruPPJVTg+1DtlXrz6ad70wsdZPHsA5QS42rQRu1Z+5mEKsiKibhKRVBdU5Hxjb495vcJSzEyXKXscjH0lWMDVmmHxUbLaFteegSwXDly9MN7gychIPJfLxUqCIAjIZrNVXHFtbe2YCdJUN2k8mzLv/cUhDbpjRaRg6O7uZufOncyaNeugRrGPB7rJrHnu3LmsW7dun1nzSJri7091VD0/VDKfIYW9qU5Vg60MNEEKmC7JCcj1RqblApQ2dIIAXI2xTYwlBjbLFYQnatIyxLFysyA0lLD5mKilFyOz0rZ4Fp/DJsNVymS40fvA0MZau1a1AGAyZ1cmQd0I/iMOdLxIXju1JgZapQVKS5QSBFrgBxIvdPFDB99zCDwXShJRdBFlifAN4GY7bYdH2n5uCCIAQk2hnObWB4/m9vuPRiqNKPi4NYJ3ve5+3vrKTThCkBYpy7uGcZv0REPZBo3o0OjIBwNwcSrcLzoGXC00NVyOO3y73XfPHIhMBUCN51JWLpfjrLi7uztuAR7JFWutpwR0D9YIfX9x44038uEPf5ienh5OP/10jj32WG69dZx2/AOM5xToRgf94YcfZvbs2axdu/agDTVGgq5Sit27d7Nz584DyppH0gs1GYfh8gjawqlIw2LALWpqW1LkvaDymM0QDW+rqrrToPIapZRRPdhONmmNZ6IPS2a3pvVXRB458XJCNI6sgK/jOHHjmwHFgJQEN3G+C2FA1lMj+zjshcW8qpIciiQOS3SU+WlQWljQlQQhKG3APwhcgsAh8CXKl1B2EUWJ9BwDuB64PYp0SRA02u1RIJRAKoUMQQYa6Rvae8HcJjrbBlAh/OQXL+DK37yYs055kve/9YF4KKYhIXSiRVijtI5Hx0cTNNCmk6ysYHuomC6VMbwRKeptdlnhdVX8e1771GhpLBfZZC+Egn1epZ6FIYQgm82SzWZHtQAnvRhaW1vj6Szbtm2r4ooPlCaY6hbgKM444wzOOOOM/b/wIOKQBt2k5Kazs5MdO3YghGDFihXMmTNn/wvYR0S0QLLVuKWl5YCB3HEc/vPrtwO3c88VH+SolfN48AljHffdT76a877yR5upEoOr8jQ0CPK+AWetFGhhMq4I54WMka2pNo3WmoGSb/VgUauZZvDkYWuzWJGGxQMa0UamZIs7MqrMC4HWAj8wWVvaragXIGqAc1EaPHv90ASxTjVaTlyss6dZrJ2wKaBW0XXCZn3aNYCLoS+0Lb6FWhCEkjB0CUIDttqXUJaIcgJwfYHTr6jPCbwGe8egjG5WKm0AN/qx2ffgniHzukAhXAe05tfXHcYNt63iqq9cy7zZBaTxNEMJ0y4sELgWHCK5F4ntTUtYI6v535j7VRUpn++B5xoKpiQUjbaRIgyduOD58MMP09DQUMWdPlt0thONsbLioaEhdu3aRVNT0yhjnJFZ8b6Sm6n00v1nxSENuhHYtra20tzczAknnEB7e/uULNtxHEqlEn/729+YNm3aqO63iUflVu5lH/geN331nbzu/KsAcEOPuS11dO4djl+mhIbaigUN0oBS1VjeJLWQkfR7fvxawHjo5kLSgUTUhdbAxiEIo8xVxeCbvJXVKKSoFM9MMU1a2VfkC2Dag12nepWkiNwcEnzweLytrjwdGaQpbdZRY7LpUJvsNgwhDB2UMr+rQIInwXMQnkR6MgZcOaRo6AevlsqYYoVJtO2PDBTS0whH8pK1S7n3r9vMygyWoLEGQoWszaBCwTs+chYnHbeTL154F46UprBotznQQdwLV32Uq0NpZWRqGJmasFNJvDDgVz+bxZvO7TbngqV3tNZ48n+oFRcDRqQf3a63tbUxPGwMjiLHrgiQ/9mjdiYbYRiSyWSYMWPGmMY4uVyuaqJFJpOpAuNoFP1Ueun+s+KQBt2tW7fieV6Vq9hkmxq01uzdu5dt27bh+z4nnXTSpAx00unKVTrU8PTTT8V/r169mm9+Yi5nX3gNYAy1dUIIqu2/Sb2u0sJoeKtmo1vqoaRwCpAOzfN7DivE/goRuCqt7ZcbQCWyX1BKWsA0TQ1CVPSzkb2DthX9wBbrbf5GZF0Q1f+q9qn9XywrjrzNdGXVQ6tQ0FqaLBeHwGa7KgAVSnQowZfgiQrgBhZwc4qGHvAyGGpF65jHlYGhFFxP4wwFUJcCL+De+7YaAroc4E6vJQwU9BVgWq3huB3J3x9fyqtOn8WlX72L447qrngiyBSRi1lRB/Fdl6cVWatIQCncaAJGJCnTCk8F9AvFyW/oBEyBsV6YcywggPQrwL/Ynj9pmpubq9rGo2nAuVwuduxKjtoZCUzPxhivG20sY5zkRIvkKPprr72WLVu2oJTiz3/+M8ccc0xVY8hk4oILLuCWW24hnU6zbNkyfvrTn05ZRi32I914VpNKvu+P6kDbvXs35XKZpUuXHvDyent72bp1K7W1tSxbtoxHH32U5z//+ZNax1KpxBv+60qGPbMr506vpavHZCvvOeMkrv7d3/ECjZIabXtIdcTnJRoVVL+C6bK6Q01BfUpSLKkKd2CRLMwF5N/mxaoEQSQZI16msNUtKcy6CVGF7wa3rBRMokm7o6wcxpMNxzGyQKaJMmYAp5LpquhxhzA04KtCC7ZKmkzVdxC+AE8gPYkMhAHUIUV9DwQSqLdXqlAjFIa79TWyrHnT85dy073bzcr0FJDTjbbzpScu4u6/7QKlkV5o0viij7SjhKUXgOMwu6mXb33jVqZPM0b5igova7YrrJrykKQWDCCHickTmv4BwcwmQdqx9ItW+Ci6CnexpP5VABTqNux7B8f7tgJMuVyOfD4/6na9rq6Obdu2sXbt2gkt85mMPXv24HnepNzBisUiP/7xj3nooYc47LDDePzxx7niiitYtmzZpNfvtttu42Uvexmu63LhhRcCcOmllx7IIsb9ZhzSoBuG4aistru7m8HBwbjrZiIxMDDAli1bSKVSLF++PPa0/etf/3pQoJvLF3jVJ34KwJvXLeTo+XVcfO3T8fMLpjewu2coJj5HAW50f68BT6McTexSY+/FnY6QWg233vgRPv3l33LfU+1QR5xq9p2ax0lXQLFCPeqYL0xgNEIoXKcCvsl/93WKREbnI/OpZMVfCCfxeksrqASlEFMMFbDVShrON5QGTQOB8E2GKzxhADcEkdc0dOlqwFVYwNXIskJ6kHqiDXnEQhQgB0pQmwYBMx1Br2ckcWzdg1xsMiVZDEyHX9FHRrcLJR+RcnjnWzdw1ps2GjWIjvrIqiWMSmsCQnwVskeFsaVjszB8cFY7OLKyX0IVEkpNqQDvPnEN1218GpAU6h4af+dPIJJFrKGhIbq6uqipqfmXNzx0dBgVz/z58ye1nJ/97Gd4nsf5558/Fas1Ztx4441cd911/OIXvziQt427Mw9pemGsOBDTm6GhIbZs2YIQglWrVlXd0kQxkakPI+PJrbvi33/ztzYWvGY5dWmHYVt1mj9nmgFdDb7AGtkCIyoyWgC1olI0CxUyD6nQ4DICbvnjYzz4VDtpgLxGoek5owxSokJj11jdLyHiSntFewtaS1s4A0MvVLrPRBWAVm9rElCrd1y0/ypTipUyr49AXEX0gnbMa5RZDx1asI0y3EAgfIkIRAVwbYbbuNfOVkgCrgYsd5vyBXJwGGZPQwEt02rpH/bNxSlUNM+dzt62fvB95IIZZgP9wNAOWkPfMMw050U6YxQUV/9sDb/59Sq+8Y07WLJ4CIGD0iGhNdYBcKTpdksLl2Vj7CIhKwA9rMtmTpqGs1ev4bNX7DLdgONarU88kkWsGTNmUCqVWLNmzZgND5OdeXYgEQTBQdZIqmNgYGDSwL2/uPLKK3nLW94yZct7zoHuRDjdXC7H1q1bCcOQFStWjEvEjzWyZ19RLpfZtm0bFAerHv/m77Zy29f+g1fZAtqDG9tMk0BtlFLKqnQyUm7FCGfvx0VggFiAMW8tKH7wy/vj9ykh8F0f5UaUr8kWjQpMI6SOtbkgLGdbDb4SC3w64lmrtzGeqZZYvViRUPmHiKuEaroiAmCtrU9EWPk8HQLarnOU3QYCQpPlSt/8LQONuyegPu/gOQKDTaa4KLSAQCMLIanQAT/Ebd1LcNRCUIrBbf3QmDIr01lgc6Te29IFK+aB1hx7zHwef7zTXOSa7YiNIDQXSK1h2KOoM3zwg6fxkudv5WMXPEwmA1JLAh3iaR8hJC4yYfReCY3JbIe0RwaJa/XVj/21HpmaxokvMXdF9/1hOie+edTbDzoiHnW8hoeRlpH5fB6lVGwZGRXtpmoS8FQA+tDQ0EGP05lIG/All1yC67qcc845k1rPZBzSoDvWgd9Xpjs8PMzWrVtj04/9mWQkR/bsK3zfp7W1lZ6eHj57/eOUrBVg8uv2iguuYvWiFjbt6jWAW1eRdeErQkk8ukFoQNu6eJS5gWkJTkE5VDAEuGakDNI0NZRSHgOvMbfW0e19haM1dEUSQyNZmNaKVMoUzqL3RDG6DpOgCpIPi+oHIwohUiZUZbhWO6YtwGpFBWxDaQ0fMIAbVADXCUAEmnRPSE1eVgBXCHthsoDbWybjpgCNu6dIcKRxn3Jae1ALDcjU+CGllhrz3qKHXD4XhKA2JQ3gag25AtTWgNa8+Zzjue4XD5tNjTR4SnHP3Uu4597FfOrTD/DCF+8mJV1c7eBpH0/7qIQqIWqE8FGkkGStNheMZO97n10BQpA1H8mX3t3ADVMIuvvrRtufZeTg4CAdHR2Uy2Vc140z4oaGhgMeszNVTQ1DQ0MHrV7YXxvw1Vdfze9+9zvuuOOOKaVeDmnQHSvGynSLxSLbtm1jeHiY5cuX09zcPKGduL+sOQgCdu7cyZ49e1i8eDH/dW3CeNwVKDunLDrNN+7uNZiUJi6ECaXNC9LEzmFamDHqOoTPvfFFfPn7d+HVCqiz6agjYZqmrARyUJPWGlzof6nNcIUFMJTJjEXE25oCXVRQ04DjGvBTOgK/5BaGcQY83u5KZrEVdcKILDcqyFmgBZPVam28JSPQJcSArLLZbczjYgA3hGy3IjM8BuAqgfA1TnuJdKORT6V6SoSzasxO2ZtDz2mOaYV3nPsifvjTv5gVLATQ6ILWPG/tEu65Z6spoLmpGFxv+tXD0UGvmA/5gdnXWvDlzz+fpsZhvvn9u5g9p0hapEgJl7L2DLVi11MKUelKk0bBEWhjun7ZzY/QtWkaQoDnATLFmS3v4/reH4x7Dh5IHIyV4niWkb7v71PKFv2MRyFMVab7TOl0//SnP3HppZdyzz33TMnwzGQc0qA7FnAmO8BKpRLbt29ncHCQZcuWMXPmzAO6Yu2rFbitrY22trZ9+i5EjtZKCHSgTZKYojortNmZ0XZFqGVbSB34n1vugyxkhSBbElz3rXfztq/8gr0DBXAkarqmFCp6XxiAG/GmEahaUb+uyFbRglAZrlfKiugBiHsqKsBrKQDNCDCuRDVfbN9ugdSsi6wAczLD1dhuCQu2OpHdhgIRJgA3NIBa0x6S9sBzxwbc1OYC7mwrHewvEGL3a6hwhItOm+YHsWEHP1Rmh8g9gzDDivY93wCu1tCTg+Z6m/GW8GutDrbsQzod36GQloZq0ZrBgRrOfetpnHhCB5/64kNksyFZYUDHw0M5IrGPjHysWAYnY667mVo48sQSWsO3PrEw1vTe/IPbeN37RpuZH2hMpd9BKpVi+vTpVXeLkZQtn8/T39/Prl27xpWyTeQOciIxlZOAk3HeeedRLpc59dRTAVi3bh1XXHHFlCz7kFYvRDKZkXH//fczY8YM+vr6OOyww5g9e/ZB3R5s2rSJlpaWWLyttaajo4OdO3cyZ84cFi9eXHXinPyR7wHwwLc/aP7+0P+a5gYhImlqtOLm30Q6qBEJD93k60zWS0HRlBfc+ovzAHjBG76FnispOJLi4QHMVMRGr7qSXkcdwklLx/iQ22zY+tTYN+xfBpaMmEKI24rlKNmZDivgOy7YKmFb0SoZruMLRADC0zT/rZ+wqQ6vIQWNIgG4IHxN01M5vLl1hoIp+Zz1stVcd89msxK7C8hZRh4meoumfbo2a2gFLSBtrlbSUwak/RAZKCMdUwrhBZXffTseOQwRoY75XqHN+aGVQipz0Tz99Vt490f+gelbMDskVIphqi1DlTBTguusJ28hD2cesxZVDhAWJKci2+3r66Ovr4/ly5dPelkHEkk/hnw+T6FQoFAo0NLSQlNT04Q6z8aLV7/61fz+979/NjZI/P9DvRBxq4VCgcbGRg4//PBJcTFR1qy1pquri+3btx9QK/D86fW05/IVqYF1fxFWfaCkjhsdBBqthO1PtcBr0UsD1EoGUiaD3zuYY3g2MFdQFAE0Gy8FpD3SOkphAV3hUKVbAd8ok45u9yNPhZF7aySAjgXKFQUEFSohmeVq6y1h2uss6CbANsSArYoyXJCeQISaVFGT2aMJ508zlEI9xEUzBZQVjTvKBC018WTk7MNtXOcY3lXuGoBZDVaVYFp9dU3KAHbnICxsMSveNQQtDZZu8KDGZLOOClGRYbwfVnZIJMhVpi1YWYc23T8ETQ0I4ObLM/zh+tN55etaef8Fm0m5Aa7j0KjH934NfDjzKDOPKwJcHfhs2bJl0q3AUzGp4WAik8mQyWRoaWmJH3vooYdYuHAhw8PDsR9DpGhIFu32J2XL5/OxxPNQiUMadKODMZJbraurY+7cuZMmvx3Hob+/n9bWVhobG6s6WVNxJQAAIABJREFU3/YVhUKB2tpaOoaHK/f1YQRowty1a20LZAaMdaARjgZHoCWRjMByBXbBGZfnffx7lJtAL0gRuhrmRfgs0MoUa0wx3NIUMSAKy6OCsOCbpBYgQS+MiJh+sH9X3Rxp45OAMguI+jqgmk6oAtsIlJWoFM2iDDeoNDSk+zSZYTNS3qvB8N5gVQpAQVPXXoSURNUagMo81Yl3/CIQgiXzGtk1aEfYaw17i+gWC3h788iFBgSEHyBarDyr6CFr7AcFASra6iBERFcWLzC/+3bKhOMgCgXSXhkvnY6PB5hb7ltvXsrSF7+b/7voer55zYPMmlMp9Ca1vR2tDu95xfOIDILMvjbPtbS0jOJPIwvFiWaKU22nOJnQWh+Ud29kGem6brxvnq1dd+PFIQ26ADt27KCjo4MFCxbEBuKdnZ0EQXBQtytR9Pf3097ejuu6rFmz5oDI9E/+5HaES6yHxYNIB3/miYu4bv1Og2COBVSlkZG3Y6hQAnQ8sjchOdBQmiZQKZPM+jNCQ12EugK8OspYtRlshkZbb10Rga8CZTnNGHOFuTCYMbVUyYZjkI3BFbOwCIwFsVrCyL6wBSaqwVabZQib5aJE/LuhFED6AulrarvB8TQiF+LNcmLHNPMeDYOK+nyILCn8eQZInX4ff9ksQwX4AUNP7oEGe5FsyyFm24yoHCCn1cZFslhr7IdI1608HhtiKKsoMaYQ01OacnsbRZ2FmbY9dzBPOZtBpIxpjQ4CnJSs8HNCMthfx7mnvwx8j5/c9yfqGuxBtUzJe175csZj9C7/jx/z+Zs/Ef+tlIr9bJOZYtT0EIFxciLwvyrTnWjsS8oWNXhE04+feuoprrvuOsrlMrfccgtr1qxh0aJFU6Iy+O///m9uuukmpJTMmjWLq666innz5u3/jROMQ5rTBWhvb6e5ubnqCv7II49w+OGHH1TVMWqYkFIybdo0lFITbit86xd/yY7uARwjMQVgycwGdnUMgRSsnF3Hnq48udDo/qMMLJ6gMOKECV2qZtaU0YQzIHBBNSpIVQAVbQt1kRY3SpA1gBnZkKSLx+s2S9b4xopR79OJn0TRLLYUizJbnaASSGS3UaZrfXDlsKK+T0KgDGUbuYWBAVwFIqdIF3yyrQW85Y3mgpLzcLUkrLPnQZLHHSpDyjUdZn5gWn0zthjmWSDSGumrCuBGnRvFMmeefSw3/N9j5gtd9nA7OyjnPZg5A7d5GmjN3FrN7q4CIpsxEr5cnhe+fCV/ua8VIQTPO20xD/52q7nt8Hy0cfKhacZurrj3KQp5OPfEl6NLJdy6NGEpRJtOFXRgMuMbBq8c56hEx6SSKUatwKVSKZZ3+b5PTU0NS5Ys+ZeD7/r16znppJMO+v1hGPLkk0/y3ve+lzPPPJPHHnuMM844g3e9612TXrekXeS3v/1t/vGPfxxMEe25y+nOnj17lP9CKpU6YNOb4eFhtmzZgu/7rFy5kqamJnp6eujr65vwMi77wGm88XO/NGYw9rv/jXe+jDd/+begNZ9/9+t46+d+aQA3OueLxkJGIdBCG05XAgicUBOmIJqJrtKC0NXoOku+RsUkCZYVriTG2mTDwnabEYLWChENm4jRt3ob9ncJVjaTNZmwjHnjCNsJqaB7BLZRZmsBuCq7VZa/DTR1wwKn3wBg4EioTSwr1IhQkR42u6hu2zCF5Y2GvvEVmb4S3kKrQthbRs60Glw/NOvpSrOvhn2oy4wCXHx7DimFCK30Yssu/uuqd/CNc69BLJhrns8XDOAC6VnTzf4YHiazcDaiu2jkfgCDOTbc/hTU1KGBB2/YXLl4yOjiA4N753H2qjmxnE9rTTDsWZ0ezFw4ne7t3fs6KnGMlylG8q62tjb6+/vp7e0FiD0Zoqx4MneGBxJTMarHcRwWL17MjBkz+MIXvjAFa1WJpAXl8PDwlLdHH/KgO1a4rjvhVuCRGt4k2e+67oTmpEUxt8XqBS0QnXHsHL5z04MxOH7o8hstrQAIwdFLZrPxyU5wDG+59vAFLG1u5Ne3Pwk1gDTZn5fSUCPQtTpuGY44XBEBr4sBOSkqmafldCNqAS0TGSgVuiCZUScz2SQwJ6iEKqCNKYdEGm2vdyLB38aAG1EKAbEDmPA12S5wfLOAIEW1oXsIctAnba9k6Y4yhUV1MZDKxzrxjp0LWvP6lx/Jn+99mmLJzBMSZQVZe5oXA2Sd1Y36qgK4gXVTC0OErxCFEqKjm3ddegbfOPvHsGKJWZdS2VAOjQ186pr38LVP3GKX67FjY9eo24ZiroSsNVMSks9pW3k86vlLefK+reYxZYE+rkiauOwvF/Hvc41ipW1LOwtXLBh13u0vInnXwMAAdXV1zJo1K6Yn8vk8e/fuHUVPRGCcpCemKqaKWx4YGHjGvHQ/85nP8LOf/YympibuuuuuKV32s5fgmWCMdUJMpBW4XC6zceNGHn300XikTxJwJ7qcUSErGPXRs0/jb5vtDDQNfUNFVKIa9cPzz8TVIv6Svf9Nz+PNrzne9EkE5keGkB0G+i3KWcCNC1HR6iVv8+16xEWq5ONQAf40kLW/iwrtapkChGN/EtukQ9CB+SGIyFwLroGlDLRRIRBE62BVCdFPQMzfOmWoa9ekfW1sFrLYYZx2G0PItA4bwBWCdFeJsN6FGtPMIP6xB46bBwhE3uP3v3ucYjEw+3SgDGn75S6HyOh33ygNogxXhgr6c8jAegkXyqz+t5Vcef61UG9nfGlNuKcHrTVuXR0/vPQ2AGobXGL/zIjPDUPmLK1kmmafC4jkhcI4jX3q6vdVnosed5yqQ1VbWxsXOL/wxsvHPOUmGkmwk1LS0NDA3LlzWbFiBccffzxr165lxYoVNDQ0xN2bDz30EBs2bGDz5s3s3r2bXC53QInIWDGVGt2DBd1TTjmFo446atTPTTfdBJj237a2Ns455xy++93vTnpdk/GczHT31QqcHJW+dOlSVq1aNe6V/EBAN5fLsXnzZtKOnaagIJvN4llHb9cRhKGudNGWzD9px/RGIOGjl93MH7/1XvOElZn956uXs27NYbznK7dRaFVwjKwAaEhMRaDtiO9QV0DVLMg8pol5XaAKhMeaAaajDDXKjJMqiuQyIrohWnBUSEtSCdG/tqtMRrKwsqZmqw8NLkGaqukZhn7Q1G/O48+vM/uqLUd9Ns2AVSq4/T7q8NkATGvMkBsoorL2FjnnIeptVusFSCcBuABBiOwv4HohyvcrcrG9/QjpsPEJ25M/13Rh6VIZWSzgNjcSeorBXqO1/crP3sVHXv5N895odlF+mCNOP4o9PzYgbSh3jSRAJ75ySdmhcCQ6CM1+HMG3Hr5uGU//dRt7d02c6hor9tcGvD96Ip/P097ezvDwMFpramtrq9QTEzVSH89L90BjMgbmE50E/La3vY3TTz+dz33ucwf1OWPFcxJ0x6IXwjBk586d7N69m8WLF48alT7ecvYHuoVCIfZzWLlyJX74YAxohUKlceMDr17Dt3/3aPx32iYLKgH4+aLPpsR0YICzXvsSstksZ750J9fd+TR7A0spSIF2dXyrnwRe01ob8b72OYmRpkULHpn5Jh9LgvDILrUIiCH+LAB8Evyr4ZgjXlfYApiZ3CAQITj9irq9iqA5ZQuAorItITj5kNo9Af48ozJweku4BcXAbPPFdnI+qk4SDXcb2JNHRiBb8BER+AYaaXqdK4DbM0S6bJQfGtD1WbO+gU+mVMTLWmXD3EoHY3pgL7/s/TFPP7KTi951FQCplOTbH/+13T+6suPKHnf/4hHktIa4MKe1xi8p3Kwb39m0bupAhyHCcZg+u4niYJ7iUCnet1prLjzlC1x008d5+7wPM9k4WLDbV/fZSCP1iehsp8p34ZlqAd6yZUtsDXvzzTezatWqKV3+IQ+645neFAoFgHh6765du5g/fz4nn3zyhA/4yPHpyYgcxYaGhmIuWAhRhWXv+vq1oKEm7dCbK1T2dqKJzgx9wVr5wfsvvQ6RMgCUSTmxLvimezajpGDaVhiaD6oetGMySB2J9UXlRyhbWLPUwagK2VjXm7GkC1FxLMpmE5xuVeHMbExCZlVRKwgLpNIa1shBQbYgCJrcxHrY7FaB0xdSOxDizUiBFMihEm9Yu5wbnzaz5WQxNNtmu8TIB8gau3PLISJlf4+42qho5imcgZzpPItmnPk+otE6zLTupux5UPQRzdMQDVZiVizyy23fBmDVcYshbVzK/NIwWx+z05kzCfNiQIcaYX0HtO8jSwWora/azT/6n/+LXz/YO8SHLn8r337v1RjvDIVWIVs27KTQVxzjYB14HOj49X1FUjsbRXIacC6XG2UZGQGx7/tT5jD2THSiffKTn2TTpk1IKVm8ePGUtf9GcciD7lgRZbq7d++mtbWV2bNn73dU+lgx1gkaBEHsKLZ06VJWr149GvgtKO3am0MARyyYwQ1/3VJZbuLloQa0plEIBjDj0UVKIENN1hqrFIs+vtKolEBnoKFPMFSvK4oFqNzSg20ntiAZYl6TKIahMXyrYHQmG4KQovr1yc2LCmfJB0eCrU6AbQBS2aJZCKlegaM1OvYQthliCKIM6ZKidlBRbLZD2AKFk8cArhCGu9Eaoq6y7hwyamzwlGlasOspAXIFUiUFNRmkhNk10JXz0DVWv9tiK9UDQ+hy2SSm5TKyzmZoSnHVfZ+KN/UT7/hJnL2SN5OaMzUuZWwGO5hjxdpFbHloF3E3YahYcfJSNj3RS7yTgK3rO+LlqkAzY2VDXNk/59LXcM0Fhl9831EVfe5kYipBd6wYbxpw0jKyo6ODgYEBwjDE87xJzXkbGBiookGmKq6//vopX2YynnOFNK01Q0NDdHZ2MjQ0xNq1a1m+fPmkr6xhGLJjxw4efPBBMpkM69atG9X19ru/mflnTgLcNPCOU4+maPwN41vzAKooED+UvO7fjjBvM438FIvm+X//1M8IHNCRWU6ySp7s8iJ6vnJrbx6y/2kRqweEY1yu4v9scUekRCU7tttg3L+oLpxB3MYrRupuA6rpBF+Y6Q17LdUcjQyOllEMcQqQKitqe0OK06yTmganNyCckTH7LlDIUFeKYyEVwPV1BXADjQxDRPsgs7UPNUYiplo76drZDxZwdV2mcvvfbmgdFRr/hWgyMLkCN/7wbgByuTI7Nhn51owZFe3d2tPWVJ0rna2D8a4TgC55bHqwvdJ9FlFLfvXtx+Y7d8UUkZ+DplkNo87vHTt20NvbS7lc7d8wkZgqLvVAI7KMXLBgAatXr2bp0qUsXLiQxYsXk8lk6Ovr46mnnmL9+vU89thjbNu2ja6uLgqFwj7lZc9UpvtMx3Mq0927dy9bt26lpqaGxsbGKeFitNa0t7ezc+dO5s6dO76jGHDrBpPN1tdk8P2QQtm0iS5fMLMClgHmUpeC//vtbfZDwPNDLnj7S/ndff+IfF8Iy+aE68oV0RksmmvIq4r8SxK3E1edoMqCp8a8YKLSyCiTjeRnZtHVz+tEZpvMckdktyjjnyBLglROmCJ/TE1o6CuTKglEg8lasz0BpUYJjplS/Klzjuaya57EC7SRhgWqYqsYqHgEDr7trFMKiiFyUwfZ6U2EQtAnsqaYNZQ375/WhBDC+l5YedbAMEhpqAYhkCsPs8fKjJW/+cq/8Kq3rOX8t/4EgJraFKX+Qlz0u/8P5mL7ktcew93X/IXhoSJieoPZ2CAA3zf0Q3RoEhmndCsNJT+/+Mb4ON1+5V+49I5P8YFjPmOoEKtFr6mpob+/n7a2NjzPI51Ox4Mc9+dV8ExnuhONqFt0LHrC87y4saOnp4dCoVA10SJqA3Yc5xlzGHum4zkButGMs3Q6zdFHH00mk2HDhokN9BsvtNZ0d3fHWsaJmNy095gMZ25zPe8/fR0f+97vQcD191cmAFvbBRBwx0Y7YSJR4Dr71DX88rZHwRX4Zc3bLrzKNEikLFoNh6RzxO8T0XsjiVjs2YDlW6nQDiN1t9HfiRpQ/LyTeD7J6SZpiGghSd7Wgq5Qhk6Qw4J02a5fBLa+JtsWQp2DarCeCV0efoNEpw0oND/cyfdSmRhwiRzAwGSyscbWrnrRo6anBGXPyLu6emD5QmIXsM4+Y0huDd+psx4MJd8Yy8xo5qR1C1j/WE+8G97xwZdw3ffvojjs8cFXXQ415j3nfeaVfO3DvzIvSllaQQg++Pk38OSdT7J39yCx90WhSKopS1Ae46oXBJz3v2/nOx/5ZfXxAQZ7hpi1wCgnhBDgSM688DXMnj2b2bONYiMJUiM51AiEI5CSUj5r2oDDMByzEUMIEZvjjKQnopbnSLb28Y9/3ChCpCQIAo499th4v0xFfP3rX+eCCy6gp6enal2mIv71R2CSUS6XaW1tZdWqVaxZsya+Ck5GS9jX18f69evp6emhsbGRZcuWTYhvypcNHbBo1nSet3px/PiVf37cfPmjVbJguHVPf+VvG48+1R5LtHQNtPblIRNldIo3HrcCUlZPCmApA1SCIgXzAQmKATBHO2kxGYHvyLMgxKgREvSBiLS4oZF9Cfu6qKuMwErBAqPFlZ7AHYS0Z1dKKSgr3KKmcX03IqVik5p0X0BYL9FZ83fNzmH6V85iaLhs3ucrZMTZBDqhsdWIIIScT7a7aG7NiyUjAZvRFJsNiXwZkU5BbY25TDRk7fEIkb0DhlbJZPj73zuQ2UxMK7zh3Bfx8/UXmc/NZkFrXvCKVXz9IwZwF6xsgQg8CkW623r41l0XGjojlr4pgrJCZzMVCidhcvOiNz7P7nTzd9MMQ5doe9zcjNkna156BG+78A1VhykCqRkzZrB06VKOPvpoTjrpJNasWcPMmTMJw5C2tjY2bNjAQw89RKFQoKOjg/7+/gPXn09hHKhO13VdmpqaWLBgAatWrWLt2rXccccdzJ07l2XLlnHnnXfy9a9/fcrWr62tjdtvv31Sk4r3FYd8ppvNZjnuuOOqHjvYDppIayul5Mgjj6S+vp4NGzaMe2UeGYFv0HPprGnG59fqb833q6IwqBWSgv1WhUBdxmS19z60lR27+2L9rZICakTMO9YMwuFLZ3HD09tNhmdZDqNisAmXTtCuAmLwBSL/2VEKBDFin0VAnFAojM5yo422mB1ltwrwIJUz7c312TTulgL5lko32dqXr+a+vX2gNW6/h3IkqsaK9nOK8pzayq1/iJnGW9U5phDDHoSCFMBAwazPUM7sq8ZaaLZc3+AwEgnN0xChQtWkK9u6aw9aw4zDatnbPQzRYNKyh2rfzWDvIE0tTSZDBlCKv9z8SLwf2tuG7QEwD5x36jc57T9fWAHcouVddUUCJnI5yBhOefW6pfEujyiGyx+6iHOXfzJ+/Ih1K3j8nqf5x183M9FwXXdMidf69etxXZeenh62b99OGIZV888aGhpIp9MH/f2ZaEwFt5xOpykWi7z97W8f1dQ02Tj//PP56le/Gs9Jm+o45EF3Kk6QQqHAli1b8Dwv9l2I4kAaJEILpNObarn/ie2VmpddRWm/WJmUQ6GsTKurC8tmzeDp1h6uv/Uxyr7CAZ539ELu3bwr/gK7fWZhv77/ERDm1j1UZriktm248egf7Acn6YTol2RWG9EIUM35BlR42+i9kR1jvLEYA/Do9t5KwmRRky4BPpyQbeTj738F7/zSbwBw+jwyXsh9PaaK7/aUESmnYlLT70NdqmKHmcxqIzVCex/pvhIssrd85QBXSvA99OAQqrER5tiK9lABsacXPasFISUq68acsNPfx5s+eSq//vLt7N3YDotMG7ExLTeavvcd80mC5KTZUtm8RgpIjKG59pGLeMvqzwDwh6sfsPvQnjOpVNypppRCBwo7TIKN922quiNLZd2qvn+AD3//XN5zxCfwy5PrApNS4jhOlVvWePPPokkPEUVRW1s7pUA8VR1puVxuygtpN998M/Pnz2fNmjX7f/FBxiEPukDV7Voy9jc+faTWdizu5kBAN+I4d3fspqt3KHqQSIQbrUljQw25oo+dX0nzrDpo7WFT6574LQ9sa7eFLA1DipS9vd82bEW+u0LEcifuizAbLEw2G92KJxUIdh1GxQiVAyN3l5WWiUSWK2MfBSo8rq/J7vSQdXY+2UDIxlI/77zkNyCMTWP9zjy5o40VohzyDeDWG6nEm089hut//2jc7EBos1qljWIhV6Lu6T2QzeAvnBFzte6wb7Y3k0bPaIbZ0+1zAXJPL+//2pn84Kt3o1wM/wqQL6LzPkGoodaFFnvctYaSD04KsWgB5ZRrMnfgWze8l4++6lvmj/iuR/PpH7zNFIYaMvjFoLIcDA8r7IA5rTWyXELXGQALggACxZ2/vj8+dy+97cKqSShKKZrnPHOFovHmnyW1tvsrZh1MTJWK4mA9HPY1BfhLX/oSt91226TXbV/xnADdsSICy7FogaTW9rDDDhtbaztiOfuLIAjQFr3KWpKtbwa6q0AsxOzw+bOaaO8eMoCWgruf2kkaKPkaHNBZiR9dRDyN6wp8TxPWVxbU3A99YSxgQAWAa7gF4VvZGcSfHysMkvZiSW43pkBEouPMPm+zXKmxz9tR5xpeOnc+j9y/EycH5ZkGcNO9PstamtjkDccuYM6QjgGXQoAoKsIW8/oV86dz/R8fr5i9RIAbGrtFp6fEcavr2bYzxfCcaRaYFe6guX2XQhBoXQFcpXA6+5l7xAIDuGmnoh7wfdOhNr2JG77/AMywWbFSiHwRVShAQwMynYrBUhDy0dd+D+w0ipg+CAK+dO7PKgfY8sSE1iLSccz7tTZ2jtKNE4RpDS4DQ3DFB69BWK554eoFVclDGIZVmXDkpvdMF8PGmvQwUmsbtQIfjFPZVGS6k3EqG68F+IknnqC1tTXOctvb2zn++ONZv349c+bMOejPGxnPCdAdK9MdC3SjwkJHRweLFi1i3bp1k24FVkrR0dHBrl27SDkOfhjSPxyiw0RaaavbOGZUz0lHLuLBJ9twqNTQvAykyqZ+pdI2U1IKt6zNLXEt6BqzrtmegJSUiM4Q5jkGEB1QgYgdw8TIlDVZYItXnoqiwGp2Ky8WsZlONdia961oquOjRy7hi99/DFeDN9ucSk4+oL7Lo3uBsro3TWZQ40+TsdbWURJtOV48zbYdffF+kpEjUDHgRatm8cRDHXjAUxvz1KxeAHnPrJ/t0vq3Vx3Jhy46nbNe/Z3Kvi4F6JSgc08eVZeJPRGOPWkxs1vque36R8w+nWb9FrQGL0Bn0sgaIzHTgIhsHqOw54qQ8MpzjuFPVz5cvY+LpoX341ecxWXn3WA8cz0PnArYaq3RwwWGIjo4k46lbyMBK51OG6DVoNExAEf/Jgtz+zuPo88+2BhrPHuyFXikkXpSPZHJZKqSmqnIdKNtmUra4+ijj6a7u2KjuWTJEv7+979PuXrhOQG6Y0VkelNTU0NyoOT+tLYjY7xW4EhStm3bNmbMmMFJJ51E/W1bKHgFtnX3cdiM5pFvAGnsXU97wRF859q/AHD4vGY27e4zEjFPoxJTbsUQpALwHU2QMhmWi3lMpWF6l2DvHB1LxaSwGa/CZL0Qy8VGfiWraIko07WyL4RARgW4EJOpWbCVnqJps6Lg5Pnx5u0G42ba06ikaNpVZmBxFvIlUJoF2qF7hgOBMq24pRAdtewGZvQ42GVrAX5IJh9CrsQjf9jMV391Lpd88ib25kvkc2Wzb3qL5sRVintu2sA9D2yvbFjRNzXAhga0K+PMdNnKmQztzfHYX7YbM/PomIShydqdymMCULkcGa9EUUtkg5l5hiO5dsN/EwQBuVyOl51xbKwp/ebbbwRg1uLpzFs8D10wF4XP/vL9HP/iI/jCW7/Dw3duNHRDygU/QLiucRXTmjM/duqoc0xKibQTgcEUjZX1cUj+C8Sz/KSUYwJx9NxURrIVeO7cufHnlEqlWMa2e/duyuVyrMttaGiYEuVEPp+nISp8HmLxnAXdqBW4q6uLbdu20dLSwkknnXTARs1jZbr9/f1s3ryZurq6qrlpC2c20T1YoLtvmCPmjWhP1BhUdDR1dRX52UuPW8Fg6Sn29OUrgKs1FDXpKDMpYvx1tWZpQxN7MIUov97ItBTWScyx4OowWps7XpJji3AGZHXlMUXFElKB9BXNj5eRPlCTwleK3XtyFJql+dBA0bSliLey1qTGSuN0h+yZJcyFwAtx8yGqwZ5yPkZVAKZgpkEM+sYuWCvqfUUZ+MTZVzLzuAWGa9UaCj5O1rUZYxlmVjKva276MLu27uEz5/3SpKOxtCxk+1PVHJ4QoIes6sFxqGmuodA1BEN5jnrBcp64q4MS4EyfZqgBrXnvp1+F67q4rhv7zm7evJnvvOemeLkXXvXvfPwFXwMgXedSM0uwfft2Hv6z1WpbUNQp19D9tsD2wPXr+fdPvoGff/mGcQ5U9HZbCEwkDfsC4uRr/hka3aRTWZIn9jwv5ok9z2PDhg0IIaqoifr6+gknQ4ODg6OKjlMdO3bseEaW+5wA3bFuMYIgYOPGjUybNm3CAyXHCtd145bLfD7P5s1GuhNJypLx2pNXsWFrJ8Mlj9nNjdWOXiEGxSQUCn6cXS6f6/KRU5fwyWuftAUwAb7GsTpZgDCrDYgoTWvPoOnS1RrSEsoaUsJMhAhEtc9CktYdsY9EUlahDEiCsAUzHfO60lfUb1dM6/DRUYGu5FFuTFGcLs1UC61J9ylyh2UrTRVdAbrFNcqKIESGuhpwI0D0FE5/mXTEl4aKm/5wPkIIzjrxC5Sm1bB796B5LlfEkY4xAU85ML3Orr9CBIq3n3a5vaMw3LB0QA/78bYLV7LqmAWcdsbRXHb+dRWiJZenmMvz2nNfwC1X3M2Tf9nKiacfy9/v2xYXY1V3D1d87Occ86LVzFk6i507d9LV1cX/XXQ3QcnsrG/dfgGXv//qeB//5JEv4HmlsRFKAAAgAElEQVQe/7nKeDdIR/CJ376dr591TdwsoZVCF4q0P11k++Ot3PD1WwFI1Uz8q7k/IPZ9n+3bt1NTUxO3nmutcRwn/u4804CcTqdpbm6mubmZrq4u1q5dSxiGMRB3dnaSz+djy8ikjG2sRGkyto7/6nhOgG4yohlnpVKJefPmcdhhh01qea7rUiqVePLJJxkeHmblypXjth6+8oTVXPzzu1DA8csSg+yskkvZbPctX/pF4jmHBctXQ/qpWCYlPF0Z5iAsl6s19VJQQhO2OCjr5SB6FLJOo0Jh/HKxb1IJemHMLFeMaKbAcrwm63QGA+raICUEbq/5ov7q2g9w9lu+j1fnUGoUZngmkOoOCOsEpC2QDmjkDDcGUamEMSbX2ma45qdxZ45yIYCZdZUCmKc445TLmDG3htIMe1HT2uhtlUBnnIrfrNYQKHsBoZKpByHv/djLePVZpvHgrUd/Bs9KrnY94nPZvZvMHnDguse/yNvXfJrCUJnfXfkXFh05l11PdbLh/u1Ie+uv+gdIN6Txch4fPP7TfPzG/2Dh0oXc/I2HaN9iOtje+IGX8rc/Psr2J4yJzb9/9vXU1dXxhbO/Fzc6XL31a4buKl2JtusjHAeZTqH8gPNf8HmkbcK57N5Pj3mOTTQiWmLv3r1s376dxYsXM3v27DizHi8rTtISzwQQJxVFjuOMmggc8cT5fJ7e3l527NgxauCmEIK+vr5nbGrEMx3PCdAVQozS2ubz+QmP7BkvIqeynp4ejjrqKI488sgJE/fHrkzoOyNg84G0Zu+QsZ1MafjtPVu5d9d95nl7S2+yVvNQOS3irHW+k6FVlMxYcKv1PO/Fx/Gd7Y/G9EDFn1ckwAgqCDuCdsAWyjTgKeo6NJmSIJ5s6YWkFHz0/51CXV2as97+PH5829/jbrK67QW8hhQ68okdCJH1TgJwSehuBVLAa09cxvqbnmDI09WAO1iCTAolobunlKAHFLImW8nWbfGrXkIpX7Sj3hPAG4bcd9MjvPwNx5FOp/nenZ/g3S/5GmhNYcjctdQ31XD1A58F4OePfYmzV12AX1a0bdmLqKsB62UhVBkN+GWIKp+XnXEV05bPJ9dnjuOaF61gwcrZfOtDPwdg7mEzeOMHXs5vLv89m/62DYDXf+jl1NfX86lXf4XAAu73H72ED510EToUyFSKeKyPELR17aI3t5dHb316IqfbqCgWizz99NOk02lOOOGEqo7KkWCaBGGtJ1ew21/srxU5yRNHioGRAzd/8pOfcMstt6C15vzzz+e4447jTW9600ENok3GxRdfzI9+9KPYuexLX/oSp5122qSWOVY8J0C3q6uLLVu2VGltS6USxeLB+ZAqpdi1axcdHR3Mnj2blpaWA+7r/t7Nfx2leXUwiiJcDEdbgHt3dMbdXyLQcROZ70DGB1VnFuIOhbT5Ba75xps5++LrYgBavRBqb1QUjxVoISq2kbZSJuIvMqAFrms8f9OuQ3GoTGpQkOlRvOH4FWzYuJOu4RLURiMmNDrtUMyEPHDPZl704sP5ya0bCBvNaZPt9PCbUoRNdnTOYICot6dUBLhgOOEAKAfU9vvccf0jBiOnG8pHCtB7C6jGjOlEixsiQmRES0Z2in4I+TIOhuqOQwjqm9Lk9xpZwKZH23nrcZ+vPgCx1Csk35vnzJWm80s6AqUcQ1lERUOtoFgyl6l0Gu15CCeNlgpZX89QrxlYeMIpR9BQn4kBN1Xj8r/3XcTN3/8zv7zEzFCbuaiZ/7j4zVzx/37OxgfMTLRT/uPfmLtkNjd0/4BLzvkuf//TE2ZXhSFOWrBu3TqKxSIXvdCO6JHwwAMPxMqAhoYGGhsbR80wi87dPXv2sHLlSpqbRxR0x4ioYJeMgy3Y7S/CMDxgudjIiRaXXHIJRx55JLt37+ZFL3oRjz766KSUGck4//zz+a//+q8pWdZ4cciPYAdD0iulqk6+3t5euru7Wb169YSXo7Wms7OT1tZW5syZw5IlS/A8j40bN3L88cdPaBnvuPRXbOm0EqhonE6gjaeszSjDjKgAYTRIMtA4AWhfx1msKmhoMhaHtXtDEIKXr1vO75822ROBYFFzlpOXzOHa+7aSX+lYyZjEFfDpk4/n2z980AxmtOCbysOtv/og6x9v5VOX/xEAt8sjo+DKq9/J2953DbpeWhOdSORr1teXmrDJZLjOHg8nLQgao4aDABkNgFQJ6VdgtLapgiZV8M1yHVHJWiMZ2ohoBF7xmqN5x3mnVD2ulOI7//1bNty/hcAPmT6jljd+4ATqpqfiLOm3P3iQDXdsq87mXXjH/3sVp551Iv9fe18e31SZtn2dLF3SJmnpBm0pbdqmdGErpYIi8ooIigv4juAw6rh9gMoygyKL4lRHAVFccMHiJ+K4wIe+issMbu8MuNGWIkVBaNN9SffSbG325/sjeQ4nbVLSNm3Skuv364+2tCfPOT3nyv1c931f9wOzt6Nb3Yc1ImeH1NfzYe3sBC/koq4vDBbg/5W/iI03PIey4ioAQFBIIA5Uv4Ltd76Owi9+ARgGsSkxeKP42V7He3FFPo4dKHD6Wq+d/jvik+JgMBigVqvZ6oDu7m7WaYzP56OlpQVRUVGQyWQelwdcETEXlIxdvbZOp0N1dTUyMzMHtZY9e/YgLCwMK1asGNRxuMjLy0NoaKinSNfllnhUkC5NFnChVqtRU1ODSZMmuXWMtrY2KBQKhIWFORjcmEwmnDp1Crm5uW4dp7u7G9dufse+MNs/jMVuUGOvi702V4ZvSyov1nHRiQp2o2+LvT3YaidQRmtGsBYAn4EhECAinq1xwE5sExMicL62DYBNi/3ug9XQ6XSora3F6peO0YvEEh2fYWA120if12VGoMp6cSJ8uL1V1mJFQBfBH/40Ax98eQpmPoFZbK9TNRDbudD8hoFcNKSh7boEgNqMYBMQHRWCB/6chtdfOw2Nxv53sjdBcO9MsSQIf99+G2TJA3OLslgsLBmp1WpotVoAQGhoKCQSCRsh0oTT2eJKvLfzCLSdXZBnJ2DlU4sRGBgIvV6Ph/9rOzqbtPalOj4GxGxmNXM6tHLC5HhER0tQdOQ0+3OhYcF4T/ESbo9ZAZPepheFj5PinfMvuTyHV1e/je/e/dHhe4EhQhxq2uvyd3Q6HRQKBbRaLUQiEYxGI+s0RiNi6jTmafQk4p71wD0TdlqtFkqlctC2q8899xwmT56MpUuXDuo4XOTl5WH//v2QSCTIycnBrl27BmMdefmRbldXF0pLS3uZ4fSEWq1GaWkpAgICkJqa2ksXIoTg+PHjuPLKK91ai8ViwexH8u2/bPtg7OTC2C0SLcG0PMvWPMCYAZ6V2KwQ7X62BpG9OgFASLNtj034DPRSm9Yq6LD1FVvFfHA9c4M6LRCAwf7X7sbSR/9hn3xphajVjC4pAwTxOdt34OkHpgPWAOzcdgyGQMAqtjEpT2tltacbb5+CA0W/AwACGQZmAwGhzGoiF31tLQCPEMBgQbCGlp+RizKDHUIhD2/vuwfh4UNfZ8nNkFMitlqtvYhYIBA4bM25UtWFtk68uOZdnP3etsNwIBWrlfVZ4O60QiNF0DR2OjxBqdOT8Py/t15yzVarFX+W/wXqZg1ufPBarNx5l9Ofo7XiNFHGNdWntcT0vHU6m+zi6g3Ik6ARcE+tGABaWlqg1WqRnJw8KJ14y5YtuPnmm3H99df36/f6agGeOXMmIiMjwTAMtm7disbGRuzbt6/fa7NjdJMu9RXl4lIRKk28mUwmyOXyPmv+fv7550uSLiEESqUS1dXVaOhi8NI/S3qRrtUCm3RASc9qz5CZAL7JTrqAjXTFsEecBDyjBcEqwBTMwCyyyREhrbYb++prZfimpNq2CCuBUE3At1phMQGWCBttHtixFF1qI9as/QhdPAC0mcH+gGYmRuDOBSl4Yrdta7vgqvG4flYqNv/tW3TzAWO8jYjDAgMhCQlEfavdV8LUw5DGYoWgwwIhq50Sh6aMoCA+nnn2D0hL81xL5UBgtVodiJjWjppMJkilUkyYMAESiaTPmu6jHxXg1bXv2nYM9g422hrsDE//ayMmX5XmsXPgJspSU1Pdsh519QYUEhLiQMT9rWV3BxaLBeXl5VCr1UhLS2MbPbjoT8Ju9erVePjhh93egfYX1dXVuOmmm3DmzJmBHsIl6Y6KRJozuGrfNRqNqKioQGdnJ1JTUz3S4kelifDwcMyYMQNXBQTg//7vGWi6bRZchG/vVXAYMW7/1z5FwmK6+McQMICBf5HMwOehW2y1G5kDPB2xfRvAtz9XAiL7DcrjwSQlsHQSINyWkJo4IQoxERFABPDlF2tx27p9UGkNIGa7yxnD4Gx1Ozbnt8PWmcBAHCaFSByE+x6cgpf+ZYtwYbBCbdJD3WV32qKWjlbb58IOMwIAlsgX/NcEMLxgTEgYg/nXZyE4eGB10kMBHo8HiUQCiUSCiIgIKBQKCIVCjB8/HkajkbU+NJvNDoTEJeK5t8/E3Ntnssfs7u7GlgXbUf1bPVsiJokKxbNHNmJ8apyzZQwIA0mUUbgq0aIG4bTD0mw2QyQSORBxIMdVrb/o6OhAWVkZ4uLiIJfLeyX+qDxBP4BLV04MxdSIxsZGtrPu008/RVZWlkePTzEqIl0ATmdGcSNUOuOsqakJSUlJveab9QVXka5Go0FpaSmEQqFTaWLWujdsn7AVBQTs4EY7oRI+bIklQhBo35JbQ3gwCQBYrJgji8aPilaANj0ACOq0/RzfSNA9xnYjbrxzLl744Bg76BIWWwPCv/euYtfz8TensefAzwAPkFgAvdYCAw+AlH+x5Mp+TUwASAiVDay2Fl3OGwaPEMBohvCCjasppmRF4InNiwddvjPUsFqtqKurQ2NjI5KTk50OOCSEsIREI2KTyeRASBKJpN8DFQcKlUqF0tJSREREICkpacgaGgghrKcCPXej0chWTtBz71k50RNmsxkKhQLd3d1IT09HsH3yxqXgTsJu0aJF+Pzzzz06mPKuu+5CSUkJGIZBYmIi8vPzWRIeAC6/SJeCa0hDR7AP5GbltlHq9Xr2ZkpLS7t0Z4zJMelCyZdHABjsHWc8BoYQgmSRCBWGbltDhYrg+pmTkT2pC7s//dlOegB9L4xMlaKuQ4sAIR83XJWOG65Kx//Z9hEUte0An8AKgr2f/IgVt82GyWTCmwd+BgPgzpum454lubjv4XdR16wB7B1VC+bKUdHciVJlBwj1bqCEC9gtFwGe1QpBq4UWWSAwkMH1/zUOWZkxkEql7MQNqpX6Gjo6OqBQKBAVFYUZM2a41DVpm2pPbwFKSB0dHaipqYHRaERwcDBLSBKJZFCRYU+YTCZUVFRAp9MhMzMTISEhHju2MzAMg5CQEISEhDjUynIrJxoaGqDX6x1mtEkkEtZ7t7W1FeXl5ZgwYQImTpzYL2MaVx12gO3Ze/HFF1FXV+fRawwA7733nkeP5wq+90QMED2dxgghMJlMKCgoYA1pBqpVCQQCWCwWWK1WVFZWoq2tDSkpKYiKinLrZhLY8mUIFdh8YKiZOc30W7oBiGzE20L0rOFNAMPg2fzvsPvJJRzCBcxjeBB0WFHXbpuU8NR9F8uqEqPDoahuZysjPvzmN3z47RnWrzxUFIB7lth0sM2P3oCH/vYxq7/+71el0IkBq4RTa8sSrm0qhEBrwYGX7kZkpBhNTU2orq5GfHw84uxm3zqdDmq1mq2d5mqG3iZig8GAsrIyWCwWTJ482e3IiwtXhKTX66FWq9HZ2Ym6ujoYDAaHyJAScX/Ip2eiLC0tzaOuWv0Bd7x6T08FSsStra3Q6XQwGAwQCAQYP348xGKxQxfaQMHj8VBSUoJ169bhlltuQVVV1ZBoz8OBUSMvGI1GlnSpIY1Op8MVV1wx6Mjg5MmTkEqlaG5uRkJCAuLi4tyKlqm8IBEIoNWbMHdyPIJ5fHxVXAMhD4gS8ZE4LgyF59phBAC6nWeA5ddMxieHS2xfh/FhtFiROSEKZ6vtgxMJ7ERoxdE3HwJg83S46S+2ibXPPrQQW178p628DLgoDTAMIsNEyE6LxTf/LoMQQLg4EGarCS187rRd7iwyQNBlwfo7r8ZNC6awFR9isRgymazP7TXVDNVqNftw0uoB7lZ1KInYHSnB0+BGhvS89Xo9AgMDHYjY1RZ9IIkyb4IQgubmZlRVVWHChAkICgpi5QmtVguGYRwi4v6Y2xgMBuzcuRNHjx5Ffn4+Jk+ePMRn4xGM7uoFwLYFU6vVrCGNXC5HWVkZMjIyBmx2QyONM2fOYOzYsUhLS3ObHAghuOove0AApEVLoVB2IkISjE+f+TPmPvQmAODdjTegtrYWuw8r0K4xwhjAsLWvkxIiILYIUHC2EQgFwDD4fvcqXNDqsGTLe/aJ3LaSrWheAD7e8wBu/cs+aLqMGCMJxscv/BkAcP1db8DIx8WqiZ4PuK07155Es3/PTFhrx0WZE7DxoRsA2G7+8vJy6PV6yOXyAVvr9UXEPcu4Bgv6BhwZGYnExMQhKZHqD6jtIT1v2txAzzs0NBStra1obm7ud6LMW9Dr9Th//jyEQiHkcrnTCJRWTtDz5pbu0b+3s8qJU6dOYd26dbjtttuwYcOGkRTdjn5NV6FQoK2tzcGQhto7DoR0Ozs7UVpaipCQEERHRyM2NtYtEuBmYnm2HgNkTxwPhbITF9SObclnq5px/exZmDJ1BhY/sh/gE1ZG+K22HWPFAUAIAIZBIGw1jhKJBAkRElQ3q20JuAAGLcSEa1a+Cb79LfLA9uUA7J16ApsHL7qtEJgI9AAggk3zAGBmGNu0YW7pF8Pgz9dPw/032wxjaMa8sbERMpkM0dHRg9ou8ng89iGjsgSXiJuamlhpgkvEEomkX9ERLQmcNGmSzyT26BadG23TLXprayvOnTsHHo8HkUiE1tZWGAwGB63Ul0B9quvr65GamtrngMj+VE5UVlbit99+Y6uC9u/fP2SVBN7AqCHdhIQEyGQyhxtTKBT22zC5q6uL1f0yMjIgFotRWuo4QNAZuGQL2PvFA4TQ6k2obdcAsDWo/fLLL4gQC9GuMSH/X2dw49xcSMQCCAQMjHax98/XTcG7351Gk8ZoI1YrgVELnPytEnyeCbXNavAAPLQwE699ZfNpJUICM4BAAR/KNhUmjIvEopU2qYHPBw6/cT9uXbHP5tFrAAx6Aks4w5IvnRSclRCFNzf8N3tera2tqKioQExMDHJzc4csUhwIEVMy7plwqa+vh1KphEwmc1t39yYYhkFbWxu6urqQm5uLkJAQGI1GNiLmzinjShMikWhYPHKdoaurC+fOnUNISAhycnIGtCvh/s3pwExCCFQqFU6dOsWOXl++fDl2796NuXPnevgsvINRIy+YzeZexKhQKCCVSh2Ef1eg9bsqlarXu3ZFRYVD4oQLZ2RLH/K7dhxEVYsKEZIgXLCPlzn0+BLwBMFYuvVDAMDRN2wlXffuPARF0wXASvDdjntQVKbElne+pS8CaAiCrAAj5sNktiI+UoJdD8+DWq3G/c//r71ui1yMWNkFArZiYftkTIaxOZHx4RDdTkmMwctrb2K3bzqdDmVlZRAIBEhNTR2wRONp9CVNCIVCtLe3IzIyEsnJyV6XEi6FvjrKnMFkMjlIEzqdjiUuSsRD1e7LXTPd9UycONGj9op6vR7bt2/H8ePHkZ+fz/ozUOczX6yE6QOjX15wBneGSlosFtTU1KCxsRFJSUlOy1ucHYe2NtJaQi7ZUsxOT0BVy2/oUOsRKGBgNBNs+38/Y/faJezPvPevYtx1Y46NcAkBuoAl6/+B3VsW20oeqCmOGNBbCWCyQMAweP9pm4TwRP4x2x/RCJi7CRBCLkavwEXTGvBsVhA8ONwOaxZfgWVzL7ZKU8NrlUoFuVzuc56lziJiqilqNBqIxWJ0dnbixIkTfUbE3kZf1ouuIBQKWSNwCm67b01NDbRaLWv8w9XHPUHEWq0W586dY5uAPHk9T5w4gfXr12PZsmU4evSoA8EyDDPSCLdPjJozcRYh0DlpzsBt242Nje1zbhrVhrm/Sy3unLUsEkJsHUMSW8MGAXDLlen4+PvfUVLeDACYlRmP42fr8faXxWju6mLPQQgCg9mK+3d+Ah6AlLgIhEmCUVTaYJMaeARmQvDgsx9hfHQYKuptjmab75uLg1/+gpoGFWAAQAiMFgACAmswgEDHc7s2Kw7r/3s2mwyj+lxdXR0SEhJ6dQ75IgghqK+vR319PVuVQNdMW33VajWUSmUvzwVvEfFgOsqcQSAQIDw83KE7i2v8U19ff0njH3fWXFVVhfb2dqSnp3t0Nll3dze2bduGEydO4P333++XK+BIxaiRFywWS69otKmpCTqdDsnJyQ7fb29vR1lZWS9HMVegx5HJZC6lBIoLFy6wskZSUhKu3WRzHPvxxZWYs9ZWtfD9bpukQKsYzEG2YzywIAenTtfiVHmzvdrg4s/OWZ0PK8/eNsyVEGirsAkX/1pCBqGSAGj1vd9wcuWxeOrOuQ7bc5PJBJPJBLFYjKSkJEilUp+KCp2hs7MTZWVlGDNmDJKSktxaL5eInZnfDDURD1dHmTNcyviHlnH1jChVKhXOnz+P6OhoTJgwwaNrLiwsxKOPPorly5dj3bp1oyqaxeUgL7gT6Wo0GpSVlYHP52PKlCluZ7QFAgG6u7thMpkcesG50Ol0KC8vByHEadfQifM17OcPvfQJ3vjrbUiODUep8gJbP3vPgum4Z8F0zFljI2MYgTcO/oDisgabC5kV+Pjp5Vj81Ic2TRawja6x+zdwwSXcQCEfq27Kxe1zLtY3isVi1ndAIBAgJSUFRqMRjY2NUCgU7M9wH0hvJW24MBqNUCgUMBgM/e7O4nouULiKiLk1pYMl4uHuKHMGV9UD3BllVB8PCQlBaGgoNBoNDAYDsrKyPLrm7u5uPPPMM/jll1/w4YcfIi3Nc0ZAIwGjJtJ1Zu+oUqlQV1eHlJQUlJeXo6urq186JdVs9Xo9SktL0dXVBaFQCIlEAqlUypYwVVVVQaVSISUlpdd2ccHmt6EzmBESIMDi3Ik4cNQ2IYBGsLP/aiNYhjD44eWVWLzlHXRoDQABBAa7JS993i1wcO2anBKJMVGhKKlpg05vhMliAZ/PIDQwAPL4KKxdPAtxUb3PlerYLS0tDhaGPX+GRkU0MuKSFk3aDJcEwZUSPFG21hd6RsQajQaEEAcSdoeI+5so8wVYrVY0NjaygyyplNbTb2Kg9bLHjx/Hhg0bcNddd2Ht2rU+v6MaBEZ/c4Qze0e1Wo3Tp0+Dz+f3q223ryQZralUqVRobm5Gd3c3RCIRoqKi2EiCK1f8dLoSG9+1VSFwJYZJSdEwWK0439AOEGIjWM6+4+XVi7D+5X9e/Av0INzn1t2EGZnx/b5GlARiY2Mxfvz4fkWvNGlDyUin04HP5zsQ8VDUk9JteXh4OGQymVceVGfSBCGk1/acrm2kdZQBrg1qPGH809XVhaeffhq//vor9u7dC7lcPlyn5S1cXqRLazVra2ttnWFXXeU2ufRMkvUkEEIIa/0XHR2NhIQEmM1mqFQq9oE0Go0QiURsNHzLMx8BAFYsnI7y+nYcLamyRbD2KoM4sQiNrTqXEa0kWAg+GIwJFeHvDy5A7Nj+W9pRaSU4OBgpKSkeIwHaCUg/uLsB+hEcHDwgIjYajWwHXFpamle25X3BarU6OHHRiJhhGBgMBiQlJSE2NnZERHNcgxp3InKu8Q/Xkzg4ONihtTssLIyNbu+99148/PDDI+J6eACXB+kaDAa2syUqKgoTJkxAcXGxW1Mf+qq3pVCpVFAoFBCJREhOTnbpckRvSBoRr33nGLrsyseHf7kBy1/418U6We5kBeqnwLnqq/8wC3+4dko/r8ZF0PpjOj6+L7N2T4HuBuhHd3c3AgMDHYi4L/MXbiXFUEsJngRNOtF2XqqXAnAZEXsbRqMRpaWlsFqtmDhx4qCcu7jGP2q1Gtu2bcPx48fR1dWFP/7xj7juuuuwaNGi0ZYwc4XRT7pWqxU///wzG8nRYv5LTX1wh2z1ej3Ky8thNBqRmpo6oJKZ2evznb24Q1QbKGCw/7Hb8d1JBW6dnYUw6cAjOxrtNzQ0ICkpCTExMV4lLq75i1qthl6vR1BQUC8i5koJSUlJI+IB5SbKJk6c2Csi5/oOOJvd5g0i5hrUyGSyfk+7vtSxf/zxR2zcuBH3338/br75ZpSUlOD06dPYunWrTyRkhwGjn3QBm4bb853aFem6Q7ZmsxnV1dVob29HcnLyoKZMHDp6Grs/d5z0OjEyGGVNthrd+69NQtpYqYNGKpVKB7Q1b29vZ/1ifcHkxRnozoTKMiqVChqNBgzDYOzYsYiMjBxWg/CBYDCJMm8SsTsGNQOFVqvF3/72N5SVlWHv3r29yjU9ifvuuw9ffvkloqOj2bE6HR0dWLZsGaqrq5GYmIhDhw55fMKEm7g8SJdr70jx888/Y9asWezD4E4nmdVqhVKpRF1dHcaPH4/Y2FiPvDubzWase/1zyOLG4JE/zHX6M8400p5bc1ctudQ3gsfjITU1dUB+scMNrpSQmJgIiUTikKwzmUy9/Hh9wWlqKBJlroi4Z/naQO/F/hjUDOTYP/zwAzZt2oQVK1Zg1apVQx7Rfv/99wgNDcXdd9/Nku5jjz2GMWPGYNOmTdixYwcuXLiA5557bkjX4QKXB+maTKZeoz0KCwuRk5MDPp9/ySQZYJt3VlFRwRbd+8L2lhsRqtVqGAwGNmEhlUohEolQX1+Pjo4OpKamjgg7QMCmgZaVlUEqlUImkzm91lx9nH5YLBavGaN7uqPsUuiLiPvT5ss1qElJSfHo9dJoNHjyySdRWVmJt956C4mJiR479qXQc4BkWloajh49inHjxqGxsRFz585FaeVo+2YAABjESURBVGnpsK2Hg9HfHOEKAoEARqMRAQEBfUoJWq0WZWVlEAqFA54qMFQIDAxEdHQ0a9xDCEF3dzdUKhVqampw4cIFCIVChIWFsVv0/tggDjdMJhNbN52eno7Q0FCXP8ud1EBH5lzKfUwqlQ7J1pzbUZabmzss2qSzpgYuEXPbfJ0R8VAa1BBCcOzYMWzevBkPPvgg9uzZ43W9trm5mb1Pxo0bh5aWFq+uxxlGNekSQiAQCFBfX4+IiAinD6LBYEBFRQW6urqQmpp66XlnPgCGYWAymVBfXw+JRIJJkyZBIBBAp9NBpVKxRETrSGnpmre7yqjfRW1tLRITE/s9O4vClQ0kJaKGhga2asATXXW+0FHGhSsi7um3YLVaYTQaIRaLkZaW5tHKFY1GgyeeeAK1tbX4/PPPMWHCBI8de7RjVJEuV7elSTKZTIa2tjbU1dVBq9WyiSqxWAytVov29nbIZDKkp6ePiLIkas5tNBp7RYl0iCKXiKg+Wltb26urjEoTw3HedMSPRCLBjBkzPC4HOGvx5UaEdXV10Gg0/bJC9KUZZZcCn89HWFgYwsLCWIOatrY2pKamwmKxQKlUoqysjB22OdA3IkIIjh49ii1btmD16tXIz8/3enTLRUxMDDtKvbGx0S1b1+HGqNN0zWZzn0kyk8mE6upqKJVKCIVCMAyDwMBANnLw1Yy51WpFTU0NmpubB2XObTabHfRRnU4HoVDInntfc7sGAiol0HKqvqSE4YCrrjpKxPSNiGb4R1JHGXBpgxpuezeNjLlWkH0RsVqtxhNPPIGGhgbk5+cjISFhuE7LJXpquhs2bEBERASbSOvo6MDOnTu9sbTLI5G2YcMGhIaGIicnB9OnT4dYLHYgj54OYAEBAU5Ll2jGnBKxN71YuR1wY8eORUJCgscjC25rs7Ma2p6tze6umysljB071mejRK45uEqlgkqlgtlsRmRkJKKjowfVVTdcsFgsqKiogFqtRnp6er8kEFc+G6Ghoejs7GR3TE899RTWrVuHe++91yei2z/+8Y84evQo2traEBMTg6eeegqLFy/G0qVLUVtbi4SEBHz00UfeSixfHqRbWlqKgoICFBYW4pdffoHRaERWVhYSExPxn//8B+vXr8fs2bMveUPSXnNu/Sg1PKFEPBxmLzS5FxgYiJSUlEF1C/UHPTuLaGszt2JAIpG4lAg0Gg3Onz8PiUQCmUzmEyVe7oCbKIuPj3eoGujq6mIHSA7FjmAw6OjoQFlZGeLi4hAfH++RNVEi/vrrr/Hmm2+ioqICCQkJyM3Nxdq1a5GRkeGBlY9qXB6k2xOtra1Ys2YNCgsLkZ2djYqKCohEIkyfPh05OTnIyclx2yOUGw2oVCrodDoIBAIHxzFPPYQ0caPRaCCXy30iuUdLt2g0rNFoYLFYHBJ1QUFBqKqqglarRVpamkfNrocSl+ooo3DWVcetoZZKpcP2xgi4NqjxBAgh+Pbbb/Hkk09i/fr1uPvuu9Hd3Y3Tp08jKSmJrRDwwyUuT9Lt7u7GZ599hqVLl7LlMx0dHThx4gQKCgpw4sQJ1NTUID4+HjNmzGDJODw83C3ydLYt59bP9tcCj2tfOBKsAGnFAK2YoB2BY8aMcaiY8NVzGKz1IpWmuETMraGmH0OhB/fXoKY/6OzsxObNm9HR0YE333yTTcwOBxITE1k5TyAQoLi4eNhe28O4PEnXHVitVlRXV6OwsBCFhYUoLi6GRqNBeno6S8JTpkxxazAj3ZZz9WEaDVIidlXI3tHRAYVC4VNNGe5Ao9GgtLQUoaGhSE5OBo/Hg0ajYa8B1/6RErEv6KNDZb3oSpqhNoj0Y6CSiycNapyt/euvv0ZeXh4effRR3HnnncOu3SYmJqK4uHhQLfc+Aj/p9gcmkwm//fYbS8S//vorBAIBsrOzkZ2djZycHKSmpro9IobWz9JtOW1eoCYv9fX1IIRALpe7Pc3C2+BKIBMnTuxTSuC2NqtUql6uY8O5LR/ujjLAdVcdJWL6ZtzXG+1QGtQAtiTzpk2boFarsWfPHnYk+nDDT7qXKen2BCEEGo0GxcXFKCwsRFFREcrLyxEdHe2gD7vr5GWxWNDZ2Ynq6mqo1WoIhUIHWWK4tcH+gBCCxsZG1NTUDGpr2zMapNtybumapxNwNFEWGRmJxMRErzeKcEfJq9Vql/PahtKghhCCI0eO4KmnnsLGjRuxfPlyr16XpKQkVt5buXIlVqxY4bW1DBJ+0vU0aEkUjYaLiorYYnRaspadnd2r+YAbsdBsM4/HY7VBGhEbDAaHSKivaoHhQk8pwdMPf3d3t8M14HosDKZ0j9t23FeizNtwtisyGo0wm80YN24cxo0b57Fx6oBN0tq4cSO6u7vx+uuv+0RyTKlUIjY2Fi0tLZg/fz5effVVzJkzx9vLGgj8pDscsFgsOHfuHAoLC3HixAn88ssvsFgsmDx5MnJychAcHIyffvoJDz300CWnEHP9FbhbUm4R/3C19ZrNZrYG1NPtpH2B67HgbFbZpTqqRuKMMgpqUCMSiRATE8NeB67PAn0z7qurzhkIIfjnP/+Jv//979iyZQvuuOMOn7wueXl5CA0NxaOPPurtpQwEftL1Bihxfvfdd9i2bRsaGhoQFxcHoVCI6dOnY8aMGcjJyUFcXJxbDw23WoCSkCf8d/taf1NTE6qrq5GQkIDY2FivP5y0tZdeA25rMyWhkdxR5o5BzWCGhra3t2PDhg0wm814/fXXPa4NDwY6nY6dxKzT6TB//nw8+eSTWLhwobeXNhD4Sdeb+OKLL2A0GnHbbbcBsNlHcmWJhoYGJCYmstpwdnY2pFKpWwRH23opCVH/XUpAA+kmA2yNGaWlpRCJREhJSfHpBgdua7NKpUJnZycsFgsiIyMRExPjU40MfUGr1eLcuXPs1Iz+SCk927u7urrYN+TW1lZIpVKUl5djx44dePzxx7Fs2TKfux6VlZVYsmQJANv5LF++HI8//riXVzVg+EnXl2G1WlFeXs6S8MmTJ9HV1YXMzEyWiLOystxOrtEkFSVibjcZd3S8M5jNZlRWVkKlUg2rlOAJUN+BqKgoxMbGOkSD3NZmX/PYoAY17e3tSE9P91hTCa0aOXjwID744APU1tYiKysLV155JR5++GGMHz/eI6/jh1P4SXekwWg0oqSkhCXiM2fOICgoCNOmTWOJWCaTuSVLcLvJerY1c01eWltbUVVVhfHjxyMuLs7nIiFXcCdRNtjW5qHCpQxqBgNCCA4fPowdO3Zg69atuP3229Ha2oqTJ08iJycHUVFRHnstP3rBT7ojHYQQdHZ24sSJE2yirrKyEnFxccjOzmY76iIjI90uW6PaaHt7Oy5cuAA+n4/o6GiEh4dDKpX6/JbcEx1lrsq2aDTsyWoBLgZjUOMOWlpa8Mgjj0AoFOLVV18ddoL96quvsG7dOlgsFjzwwAPYtGnTsL6+D8BPuqMRVqsVdXV1KCgoQFFREU6cOIHOzk6kpaWxibopU6a4TK5RKYH+jkgkcnDaoltyrj7sK9ruUHWU9UxWarVaMAzjsCsYrNnRUBjUUBBC8Mknn2Dnzp3Iy8vDbbfdNuxvnBaLBXK5HN9++y3bYn/gwIHLzSTHT7qXC8xmM86ePct6S5SUlIBhGEydOpVt5EhNTcXhw4cRFxeHhIQEl1ICd0tOSchsNrO2l9wC/uGCNzrKLBaLQzSs1WohEAgcGjncqRoZSoMawDaq5pFHHkFwcDBeeeUVr3V1HT9+HHl5efj6668BANu3bwcAbN682Svr8RJG74y0jz76CHl5eTh37hyKioqQk5PD/t/27dvx9ttvg8/nY/fu3ViwYIEXVzo8EAgEmDJlCqZMmYKVK1eCEAKtVouTJ0+isLAQW7ZsQXFxMZKTkzFnzhx2aKczv1uGYRAcHIzg4GC2tIgeT61WQ6lUsm3NPetGhyK66uzsRGlpKaKiooZtRhlgm8oQHh7uMMqbJqmo2c+lWpu5BjUDHVPkClarFf/zP/+DF154AU8//TQWL17sVVmooaHBIUkXHx+PwsJCr63H1zDiSTcrKwuffPIJVq5c6fD933//HQcPHsTZs2ehVCpx3XXXoayszGeHNQ4VKCHOnTsXKSkp+OKLL3DkyBEkJCSgqKgIBQUFePvtt9HU1ISUlBRWlpg2bZpThzB6PO58Mlo3qlKpUFlZ6TCNgmv7OFBwE2VZWVk+0VEmFAoRERHhMMacm6irq6uD0WhEUFAQ9Ho9BAIBJk+e7PG1NzU1Yf369RCLxTh69KhHx6oPFM52z76cGxhujHjSTU9Pd/r9zz77DHfccQcCAwORlJSElJQUFBUVYdasWcO8Qt9BfHw8fvjhB/YBuPXWW3HrrbcCsBFnWVkZCgoK8NlnnyEvLw9GoxGTJk1iiTgjI8Oppsudz0VhNBpZSaK+vn5A3go9E2WejhA9jaCgIAQFBSE6OpptLKmsrERERAQYhsHvv//OtjYPVp6xWq04dOgQXnrpJTzzzDO45ZZbfObaxMfHo66ujv26vr7eawY6vogRT7qu0NDQgJkzZ7Jfx8fHo6GhwYsr8g24ejD5fD7S09ORnp6Oe++9F4Atcjt16hQKCgrw6quv4uzZswgNDXUw+XE1PiggIABRUVFs1pzrrdDW1obKyspeJujcSgFuomz69Ok+U1PrDrgGNbm5uQ5vLtzWZirPAOiVqOtLOmlqasK6deswZswYHDt2zFvjaFxixowZUCgUrL/IwYMH8eGHH3p7WT6DEUG61113HZqamnp9/9lnn2UjtZ7wb3EGj6CgIMyaNYvdHRBC0N7ezprAHzx4kJ1FRU1+pk+f7tQEnmEYiEQiiEQijB07FoDj2PT6+np2SCJgm9KQnJw8ovwSCCFoaGhAXV0d5HK5062+s/Hx3Lbempoadmo11cnFYjGbqDt48CB2796Nbdu2YdGiRT55bQQCAV577TUsWLAAFosF9913HzIzM729LJ/BiCDd7777rt+/49/ieB4MwyAyMhI33HADbrjhBgAXu6kKCwvx73//Gzt37oRWq0VGRgYbEU+ePNmppsv1C4iPj0dnZyfOnz8PsViMqKgotLa2oqamBgEBAQ76sC/aXlKDmpCQkH6PmHcmz9BhmbR5Yu3atSCEsAYw06ZN80nCpbjxxhtx4403ensZPolRUzI2d+5cvPDCC2z1wtmzZ7F8+XIUFRVBqVRi3rx5UCgUl10izRswGo0OJvC//fYbhEIhpk2bxurDKSkpDlFtZWWly46yntOa6SQGrj7sLdtLdwxqBgOr1YoPP/wQr732Gh577DGEhYWhuLgYPB4PTz75pEdfyw+PYvTW6X766adYs2YNWltbERYWhqlTp7L1gc8++yz27dsHgUCAl19+mY3Ohhp5eXl46623WD1z27Ztl/W7PiEEarXawQS+oqIC0dHRiIiIwJkzZ7B//35kZGT0u62ZVgvQKJA7rXmoS8oGY1DjDpRKJdauXYtx48Zh165dHid0d+G/nweE0Uu6vogR7gM6LGhubsaf/vQnALayv5KSErS3t0Mul7P68LRp03qZwLuC1Wp1mNZMddGhmM02VAY13OO///772LNnD5577jksWLDAq1KC/34eEEZvc4QfIxNisRh5eXmYPXs2+z2LxYLff/8dhYWF+Pjjj/H444+DEMKawOfk5GDixIlOpQQej8dGubQwnzubjTYwDNZpjGtQk5OT4/Four6+HmvXrsX48ePx/fffQyqVevT4vgJCCK6++mo8/vjj7A700KFD2LdvH7766isvr25o4Y90hwB5eXnYv38/JBIJcnJysGvXLoduJj/cA5URTp48iaKiIhQWFqK0tBTh4eFspcSMGTPcdkSjI9O5+rDJZGLrZvsaCTTUBjVWqxX/+Mc/kJ+fj+effx7z58/3mUTZUN3PZ86cwe23345Tp07BYrFg6tSp+Oqrr5CcnOyBVXsdfnnB0+irjG3mzJms29fWrVvR2NiIffv2eWGVow+EELS2tjqYwCuVSiQlJTmYwEskEreJmDuXTK1WA7g4DkcqlcJgMEChUAyJQQ0A1NXVYc2aNZDJZNi5c6dXPIy9dT8/9thjCAkJgU6ng1gsxtatWz1yXB+An3S9herqatx00004c+aMt5cyakFN4Knb2smTJ6HX63uZwLsrJdC62QsXLqChoYH13aWWl9RXYbDka7VasX//frz11lvYtWsX5s2b5zPRrSt4+n7W6XTIzs5GQEAAiouLfbIccIDwa7rDicbGRnay6qeffoqsrCwvr2h0g8fjQS6XQy6X4+677wZgKzOjJvD5+fmsCXx2djZLxElJSU41WT6fD5PJhKamJshkMowbN87B4EapVEKv1yM4ONhBH+6P7WVtbS1Wr14NuVyOn376CaGhoR67Hp7GUN7PISEhWLZsGUJDQ0cT4fYJP+kOAR577DHWUjExMRH5+fnD+vp+A2kgMDAQV1xxBa644goAF03gqTb8ySefsG2qlISnT58Oo9GIY8eOISMjA9nZ2SwRBAQEIDIykrVLpLaX1ASe29bMHRnfk9StVivefvttvPPOO9i1axeuvfZan49uh/p+5vF4w+YY5wvwywujDH4DafdBvXkLCgpQWFiII0eOoK2tDVdffTVmzpzpYALv7vG4+jC1vQwODsb333+PlJQU7N69G5mZmdixY4dPR7fDiVFakuaXFy4XFBUVISUlBTKZDABwxx134LPPPvOTrhPweDwkJiYiMTER7e3tuPLKK7Fjxw4olUoUFBTggw8+wIYNG8Dj8dhuupycHMjlcqcVDlxfBQqLxQKlUoni4mK8+eabYBgGAoEA+fn5eOSRR4bzdP3wEfhJd5TBbyA9MKxatYol0ujoaEydOhWrVq1yMIEvKCjAtm3bUFZWhqioKJaEZ8yYgZiYGKcyQW1tLdasWYNJkybh119/RUhICJqbm1FbWzvcp+izyMvL8/YShhV+0h1l8LurDQyuWni5JvBz584FYLvGjY2NrAn83r170dLSwprA5+TkYMqUKThw4ADee+89vPLKK7j66qvZv0NMTAw7iWMo4Z+q4pvwk+4og99dbejBMAxiY2OxePFiLF68GIBNRigtLUVhYSEOHz6MVatWITc3Fz/99BNEIpFX1umfquKb8JPuKIPfQNo74PP5yMjIQEZGBu69914QQry+w/BPVfFNXD51GpcJuAbS6enpWLp0qd9A2gvwNuH2BWe6v3+qyvDBH+mOQviSgXRiYiLrZyAQCFBcXOztJY0q+KeqjDz4SdePIcd//vMftqnAD8/CP1Vl5MEvL/jhx2WGW265BQcPHoTBYEBVVRUUCgVyc3O9vazLBn7S9WNIwTAMrr/+ekyfPh179+719nIuK3z66aeIj4/H8ePHsWjRIrYsLDMzE0uXLkVGRgYWLlyI119/3V+5MIzwtwGPYNTV1WHOnDk4efIkxowZgwsXLiA7OxtHjx7FhAkTvL08ALaRM7GxsWhpacH8+fPx6quvYs6cOd5elh9+DDVciuT+SHcEY/z48XjwwQdZQ5tNmzZhxYoVPkO4AFitMDo6GkuWLEFRUZGXV+SHH96Fn3RHOP7617+ioKAAL7/8Mn788Uef6ufX6XTQaDTs5998843f5rIf+Oijj5CZmQkej+dQ9VFdXY3g4GBMnTqVbVf2Y+TAX70wwiEUCvH8889j4cKF+Oabb/o982so0dzcjCVLlgAAzGYzli9fjoULF3p5VSMHrjrKACA5ORklJSVeWJUfg4WfdEcBjhw5gnHjxuHMmTOYP3++t5fDQiaT4fTp0157/fvuuw9ffvkloqOj2UkHHR0dWLZsGaqrq5GYmIhDhw757Pw6Vx1lfoxs+OWFEY6SkhJ8++23KCgowEsvvYTGxkZvL8lncM899/SaLLtjxw7MmzcPCoUC8+bNw44dO7y0usGhqqoK06ZNwzXXXIMffvjB28vxoz8ghPT14YcPw2q1kpkzZ5JvvvmGEELI7t27yfLly728Kt9CVVUVyczMZL+Wy+VEqVQSQghRKpVELpd7a2mEEELmzZtHMjMze30cPnyY/ZlrrrmGnDhxgv1ar9eTtrY2QgghxcXFJD4+nqhUqmFfux99wiWvXqpkzA8fBsMwKwDMI4Qss3/NB1AEYD0h5JhXF+cjYBgmEcCXhJAs+9edhJAwzv9fIIT4pr5gB8MwRwE8Sghx2kN9qf/3w7fg13RHMAghewHs5XxtATDdeyvyYzjAMEwUgA5CiIVhGBmAVACVXl6WH27Cr+n6cbmhmWGYcQBg/7fFy+txCYZhljAMUw9gFoB/Mgzztf2/5gD4lWGY0wA+BrCKENLhrXX60T/45QU/RjWcyAvPA2gnhOxgGGYTgDGEkMe8uEQ/LjP4SdePUQuGYQ4AmAsgEkAzgL8BOAzgEIAEALUAbvdHiX4MJ/yk64cffvgxjPBrun744Ycfw4j/DxiuSoZo+jnUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.axes(projection='3d')\n",
    "\n",
    "viridis = cm.get_cmap('viridis', len(tha))\n",
    "for i in range(len(tha)):\n",
    "    colour = viridis(i/len(tha))\n",
    "    ax.plot(tha['pos_x_gsm'][i:i+2]/6378.1,\n",
    "            tha['pos_y_gsm'][i:i+2]/6378.1,\n",
    "            tha['pos_z_gsm'][i:i+2]/6378.1,\n",
    "            c=colour)\n",
    "    \n",
    "ax.set_xlabel('X')\n",
    "ax.set_ylabel('Y')\n",
    "ax.set_zlabel('Z')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-06T20:14:19.597091Z",
     "start_time": "2019-09-06T20:14:19.556053Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "124801 27886 10595 9694\n"
     ]
    }
   ],
   "source": [
    "equa_r = tha[(np.abs(tha['pos_z_gsm']/6378.1)<1) &\n",
    "             (tha['r'] >= 7)]\n",
    "dawn = equa_r[(3 <= equa_r['mlt']) &\n",
    "                  (equa_r['mlt'] <= 9)]\n",
    "dusk = equa_r[(15 <= equa_r['mlt']) &\n",
    "                  (equa_r['mlt'] <= 21)]\n",
    "print(len(tha),len(equa_r),len(dawn),len(dusk))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.1"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}