examples/ex1_calculate_events.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Calculate events in Circular Restricted Three Body Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: fastrk in c:\\users\\stasb\\pycharmprojects\\fastrk (0.0.2)\n",
      "Requirement already satisfied: orbipy in c:\\users\\stasb\\anaconda3\\envs\\p38\\lib\\site-packages (0.2.5)\n",
      "Collecting numpy~=1.19.2\n",
      "  Downloading numpy-1.19.5-cp38-cp38-win_amd64.whl (13.3 MB)\n",
      "Collecting numba~=0.51.2\n",
      "  Using cached numba-0.51.2-cp38-cp38-win_amd64.whl (2.2 MB)\n",
      "Requirement already satisfied: sympy~=1.6.2 in c:\\users\\stasb\\appdata\\roaming\\python\\python38\\site-packages\\sympy-1.6.2-py3.8.egg (from fastrk) (1.6.2)\n",
      "Requirement already satisfied: matplotlib>=3.0.0 in c:\\users\\stasb\\appdata\\roaming\\python\\python38\\site-packages (from orbipy) (3.4.1)\n",
      "Requirement already satisfied: scipy>=1.1.0 in c:\\users\\stasb\\anaconda3\\envs\\p38\\lib\\site-packages (from orbipy) (1.6.3)\n",
      "Requirement already satisfied: pandas>=0.23.4 in c:\\users\\stasb\\appdata\\roaming\\python\\python38\\site-packages (from orbipy) (1.1.4)\n",
      "Collecting llvmlite<0.35,>=0.34.0.dev0\n",
      "  Using cached llvmlite-0.34.0-cp38-cp38-win_amd64.whl (15.9 MB)\n",
      "Requirement already satisfied: setuptools in c:\\users\\stasb\\anaconda3\\envs\\p38\\lib\\site-packages (from numba~=0.51.2->fastrk) (49.6.0.post20200814)\n",
      "Requirement already satisfied: mpmath>=0.19 in c:\\users\\stasb\\appdata\\roaming\\python\\python38\\site-packages\\mpmath-1.1.0-py3.8.egg (from sympy~=1.6.2->fastrk) (1.1.0)\n",
      "Requirement already satisfied: python-dateutil>=2.7 in c:\\users\\stasb\\appdata\\roaming\\python\\python38\\site-packages\\python_dateutil-2.8.1-py3.8.egg (from matplotlib>=3.0.0->orbipy) (2.8.1)\n",
      "Requirement already satisfied: pyparsing>=2.2.1 in c:\\users\\stasb\\appdata\\roaming\\python\\python38\\site-packages\\pyparsing-3.0.0b1-py3.8.egg (from matplotlib>=3.0.0->orbipy) (3.0.0b1)\n",
      "Requirement already satisfied: pillow>=6.2.0 in c:\\users\\stasb\\appdata\\roaming\\python\\python38\\site-packages\\pillow-8.0.1-py3.8-win-amd64.egg (from matplotlib>=3.0.0->orbipy) (8.0.1)\n",
      "Requirement already satisfied: cycler>=0.10 in c:\\users\\stasb\\appdata\\roaming\\python\\python38\\site-packages\\cycler-0.10.0-py3.8.egg (from matplotlib>=3.0.0->orbipy) (0.10.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in c:\\users\\stasb\\appdata\\roaming\\python\\python38\\site-packages\\kiwisolver-1.3.1-py3.8-win-amd64.egg (from matplotlib>=3.0.0->orbipy) (1.3.1)\n",
      "Requirement already satisfied: pytz>=2017.2 in c:\\users\\stasb\\anaconda3\\envs\\p38\\lib\\site-packages (from pandas>=0.23.4->orbipy) (2021.1)\n",
      "Requirement already satisfied: six>=1.5 in c:\\users\\stasb\\appdata\\roaming\\python\\python38\\site-packages\\six-1.15.0-py3.8.egg (from python-dateutil>=2.7->matplotlib>=3.0.0->orbipy) (1.15.0)\n",
      "Installing collected packages: numpy, llvmlite, numba\n",
      "  Attempting uninstall: numpy\n",
      "    Found existing installation: numpy 1.20.2\n",
      "    Uninstalling numpy-1.20.2:\n",
      "      Successfully uninstalled numpy-1.20.2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Keyring is skipped due to an exception: entry_points() got an unexpected keyword argument 'group'\n",
      "ERROR: Could not install packages due to an EnvironmentError: [WinError 5] Access is denied: 'C:\\\\Users\\\\stasb\\\\AppData\\\\Roaming\\\\Python\\\\Python38\\\\site-packages\\\\~umpy\\\\.libs\\\\libopenblas.GK7GX5KEQ4F6UYO3P26ULGBQYHGQO7J4.gfortran-win_amd64.dll'\n",
      "Consider using the `--user` option or check the permissions.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!pip install fastrk orbipy"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\stasb\\Anaconda3\\envs\\p38\\lib\\site-packages\\numpy\\_distributor_init.py:30: UserWarning: loaded more than 1 DLL from .libs:\n",
      "C:\\Users\\stasb\\Anaconda3\\envs\\p38\\lib\\site-packages\\numpy\\.libs\\libopenblas.NOIJJG62EMASZI6NYURL6JBKM4EVBGM7.gfortran-win_amd64.dll\n",
      "C:\\Users\\stasb\\Anaconda3\\envs\\p38\\lib\\site-packages\\numpy\\.libs\\libopenblas.WCDJNK7YVMPZQ2ME2ZZHJJRJ3JIKNDB7.gfortran-win_amd64.dll\n",
      "  warnings.warn(\"loaded more than 1 DLL from .libs:\"\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.integrate import solve_ivp\n",
    "from fastrk import BT8713M, RKCodeGen, EventsCodeGen\n",
    "from model_crtbp import crtbp\n",
    "from timeit import timeit\n",
    "from numba import njit\n",
    "import orbipy as op"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [],
   "source": [
    "# Orbipy uses scipy.integrate.ode\n",
    "model = op.crtbp3_model()\n",
    "\n",
    "# Events\n",
    "\n",
    "ev_y = op.eventY(0., terminal=False)\n",
    "ev_z = op.eventZ(0., terminal=False)\n",
    "ev_vx = op.eventVX(0., terminal=False)\n",
    "ev_vy = op.eventVY(0., terminal=False)\n",
    "ev_vz = op.eventVZ(0., terminal=False)\n",
    "\n",
    "events = [ev_y, ev_z, ev_vx, ev_vy, ev_vz]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [],
   "source": [
    "# generate code and import rk_prop function\n",
    "# new module rk_bt8713m.py will be created (at first call) or loaded\n",
    "# CRTBP is autonomous system i.e. equations doesn't have explicit time dependency\n",
    "\n",
    "rk_module = RKCodeGen(BT8713M, autonomous=True).save_and_import()\n",
    "rk_prop = rk_module.rk_prop\n",
    "rk_prop_ev = rk_module.rk_prop_ev\n",
    "\n",
    "# after first run code above can be replaced by\n",
    "# from rk_8713M import rk_prop, rk_prop_ev\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [],
   "source": [
    "# generate code and import call_event function;\n",
    "# new module __evcache__/ev_<hash>.py will be created (at first call) or loaded;\n",
    "# call_event(t, s, values, idx) function calls specific event function by its index in events list\n",
    "\n",
    "call_event = EventsCodeGen(events).save_and_import()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [],
   "source": [
    "# Let's define Cauchy's problem:\n",
    "# initial time and state for halo orbit\n",
    "\n",
    "t0 = 0.\n",
    "t1 = 3*np.pi\n",
    "s0 = np.zeros(6)\n",
    "s0[[0, 2, 4]] = 9.949942666080747733e-01, 4.732924802139452415e-03, -1.973768492871211949e-02\n",
    "mc = np.array([3.001348389698916e-06])"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [],
   "source": [
    "# integration parameters, same as Orbipy default\n",
    "params = {'max_step': np.inf,\n",
    "          'rtol': 1e-12,\n",
    "          'atol': 1e-12,\n",
    "          }"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [],
   "source": [
    "# event detector object\n",
    "dtr = op.event_detector(model, events)\n",
    "\n",
    "df, evdf = dtr.prop(s0, t0, t1)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [],
   "source": [
    "values = np.array([e.value for e in events])\n",
    "terminals = np.array([e.terminal for e in events])\n",
    "directions = np.array([e.direction for e in events])\n",
    "counts = np.array([e.count for e in events])\n",
    "accurates = np.array([e.accurate for e in events])\n",
    "\n",
    "trj, evarr = rk_prop_ev(crtbp, s0, 0, t1, *params.values(),\n",
    "                        values, terminals, directions, counts, accurates,\n",
    "                        call_event, 1e-12, 1e-12, 100, mc)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [],
   "source": [
    "evdf1 = pd.DataFrame(evarr, columns=['e', 'cnt', 't', 'x', 'y', 'z', 'vx', 'vy', 'vz'])"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [
    {
     "data": {
      "text/plain": "0    -4.922285e-12\n1     4.164447e-12\n2     7.951861e-12\n3     3.029155e-11\n4     4.951595e-13\n5     5.851453e-11\n6     1.755840e-11\n7     1.680007e-10\n8    -6.343224e-10\n9     1.618141e-10\n10    4.295067e-10\n11    1.630771e-10\n12   -7.911112e-10\n13    6.959928e-09\n14    2.087988e-09\n15    1.687998e-08\n16    8.396579e-09\n17    5.646641e-08\n18   -2.277712e-07\n19    5.820032e-08\n20    1.568987e-07\n21    5.598841e-08\n22   -2.966316e-07\n23    2.560417e-06\n24    7.688840e-07\n25    6.216301e-06\n26    3.096731e-06\n27    2.080954e-05\n28   -8.392425e-05\n29    2.144755e-05\nName: t, dtype: float64"
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check event times\n",
    "evdf.t - evdf1.t"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [
    {
     "data": {
      "text/plain": "2     2.089539e-17\n7    -1.923244e-18\n12   -3.282024e-16\n17   -1.551711e-13\n22    1.787718e-16\n27   -2.841581e-13\nName: y, dtype: float64"
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evdf[evdf.e == 0].y - evdf1[evdf1.e == 0].y"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 1080x360 with 3 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAFgCAYAAABNIolGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAu3UlEQVR4nO3df7xcdX3v+9ebROC2VkWTkghIUIPKKS3iFtp62lrESj29QKk/oD1HaLE8Tlvs6fVogg9b9Vg9h+TeXq23tpqLVLRVVA7WeMVSRFBPlR8bi/wsJGKooYlEfvjjYUUCn/vHrIRJ2Dt7kr1n1tp7Xs/HYx57re/3u9Z8JjPzmcln1vquVBWSJEmSJEldtl/bAUiSJEmSJM3EAoYkSZIkSeo8CxiSJEmSJKnzLGBIkiRJkqTOs4AhSZIkSZI6zwKGJEmSJEnqvAVbwEhyYZJ7k9wy4PhXJbktya1JPjLs+CRp3JiXJal7zM2S5pNUVdsxDEWSXwS+D3yoqn5qhrErgY8DJ1TVA0l+sqruHUWckjQuzMuS1D3mZknzyYI9AqOqvgjc39+W5FlJ/j7JDUm+lOS5TdfvAu+tqgeabU3EkjTHzMuS1D3mZknzyYItYExjHfC6qnoB8AbgL5v2I4Ejk/xjkmuSnNRahJI0XszLktQ95mZJnbS47QBGJckTgZ8HPpFkR/MBzd/FwErgxcChwBeTHF1VD444TEkaG+ZlSeoec7OkLhubAga9o00erKpjpujbDFxbVQ8D30hyJ73kfP0I45OkcWNelqTuMTdL6qyxOYWkqr5LL9G+EiA9P9N0/x29SjJJltA7PO6uFsKUpLFhXpak7jE3S+qyBVvASPJR4CvAc5JsTnI28FvA2Um+BtwKnNIMvxy4L8ltwFXAG6vqvjbilqSFyrwsSd1jbpY0nyzYy6hKkiRJkqSFY8EegSFJkiRJkhaOBTeJ55IlS2rFihVthyFJU7rhhhu+XVVL245j1MzNkrpqXPMymJslddd0uXnBFTBWrFjB5ORk22FI0pSS3N12DG0wN0vqqnHNy2BultRd0+VmTyGRJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJGkMJTkpyR1JNiY5b4r+dyW5sbndmeTBvr5H+vrWjzRwSVrAzM2StGcLbhJPSdKeJVkEvBd4KbAZuD7J+qq6bceYqvo/+sa/Dnh+3y7+raqOGVG4kjQWzM2SNDOPwJCk8XMcsLGq7qqqHwEXA6fsYfwZwEdHEpkkjS9zsyTNwAKGJI2fQ4Bv9q1vbtoeJ8nhwBHA5/uaD0wymeSaJKdOdydJzmnGTW7btm0OwpakBc3cLEkzsIAhdciKFW1HID3O6cAlVfVIX9vhVTUB/Cbw7iTPmmrDqlpXVRNVNbF06dJRxCppzCVtRzAy5mZ1zhi9/9QiCxhSh9x9d9sRaEzcAxzWt35o0zaV09ntEOWquqf5exdwNbuegy1J2jfmZkmagQUMSRo/1wMrkxyRZH96X4QfN2N9kucCBwFf6Ws7KMkBzfIS4EXAbbtvK0naa+ZmSZqBBQypK1av5g5WwurVbUeiBa6qtgPnApcDtwMfr6pbk7w9ycl9Q08HLq6q6mt7HjCZ5GvAVcD5/TPkS1Jrli/nEQLLl7cdyT4xN2tem+fvP80fXkZV6oLVq2HtWlYCrF3ba1uzps2ItMBV1WXAZbu1vWW39bdNsd2XgaOHGpwk7a3ly2HrVgKwdWtvfcuWtqPaa+ZmzUsL5P2n+cEjMKSWrFjRm+wogTvXXkoBAYre+o6+Lkzs2YUYJEnqt+NzMoFHt27d5XP00a1bd+mXNLd8/6ktFjCklmzaBFW925GrTtuZ9ENvfUffpk2thgk4uagkqXt2fE5WwX7Llu3yObrfsmW79EuaW77/1BYLGFIXrFkDq1axgWfDqlWePiJJ0t7YsgWWLaMAli3z8HVplHz/aYQsYEhdsWYNz2FD94oXTi4qSZoPtmxhEeV/nqQ2+P7TiDiJp6TpObmoJEmSpI7wCAypQw4/vO0I5tfkopIkSZLGhwUMqUO6MGHnfJpcVJKkfk4YKLXH959GwVNIJE2vOV1kw9pLOXLVaZ4+IkmSJKk1HoEhac+6OrmoJEmSpLFiAUOSJEmSJHWeBQxJM+rC5KKSJEmSxlurBYwkFya5N8kt0/QnyXuSbExyU5JjRx2jJCfslCRJktS+to/A+CBw0h76fxVY2dzOAf5qBDFJkiRJkqSOabWAUVVfBO7fw5BTgA9VzzXAU5IsH010kiRJkiSpK9o+AmMmhwDf7Fvf3LTtIsk5SSaTTG7btm1kwUmSJEmSpNHoegFjIFW1rqomqmpi6dKlbYcjSZIkSZLmWNcLGPcAh/WtH9q0SZIkSZKkMdL1AsZ64DXN1Uh+FvhOVW1pOyhJkiRJkjRai9u88yQfBV4MLEmyGXgr8ASAqnofcBnwcmAj8APgt9uJVJIkSZIktanVAkZVnTFDfwF/MKJwJEmSJElSR3X9FBJJkiRJkiQLGJIkSZIkqfssYEiSJEmSpM6zgCFJkiRJkjrPAoYkSZIkSeo8CxiSJEmSJKnzLGBIkiRJkqTOs4AhSWMoyUlJ7kiyMcl5U/SflWRbkhub22v7+s5MsqG5nTnayCVp4TI3S9KeLW47AEnSaCVZBLwXeCmwGbg+yfqqum23oR+rqnN32/apwFuBCaCAG5ptHxhB6JK0YJmbJWlmHoEhSePnOGBjVd1VVT8CLgZOGXDblwFXVNX9zRfjK4CThhSnJI0Tc7MkzcAChiSNn0OAb/atb27advcbSW5KckmSw/ZyW5Kck2QyyeS2bdvmIm5JWsjMzZI0AwsYkqSpfBpYUVU/Te+XvIv2dgdVta6qJqpqYunSpXMeoCSNIXOzpLFmAUOSxs89wGF964c2bTtV1X1V9VCzegHwgkG3lSTtE3OzJM3AAoYkjZ/rgZVJjkiyP3A6sL5/QJLlfasnA7c3y5cDv5LkoCQHAb/StEmSZsfcLEkz8CokkjRmqmp7knPpfbldBFxYVbcmeTswWVXrgT9McjKwHbgfOKvZ9v4kf0rvizbA26vq/pE/CElaYMzNkjSzVFXbMcypiYmJmpycbDsMSZpSkhuqaqLtOEbN3Cypq8Y1L4O5WVJ3TZebPYVEkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkKQxlOSkJHck2ZjkvCn6X5/ktiQ3JbkyyeF9fY8kubG5rR9t5JK0cJmbJWnPFrcdgCRptJIsAt4LvBTYDFyfZH1V3dY37J+Aiar6QZLfA9YCr276/q2qjhllzJK00JmbJWlmHoEhSePnOGBjVd1VVT8CLgZO6R9QVVdV1Q+a1WuAQ0ccoySNG3OzJM3AAoYkjZ9DgG/2rW9u2qZzNvDZvvUDk0wmuSbJqUOIT5LGkblZkmbgKSSSpGkl+Y/ABPBLfc2HV9U9SZ4JfD7JzVX19Sm2PQc4B+AZz3jGSOKVpHFgbpY0rjwCQ5LGzz3AYX3rhzZtu0hyIvBm4OSqemhHe1Xd0/y9C7gaeP5Ud1JV66pqoqomli5dOnfRS9LCZG6WpBlYwJCk8XM9sDLJEUn2B04HdpmxPsnzgffT+4J8b1/7QUkOaJaXAC8C+ieYkyTtG3OzJM3AU0gkacxU1fYk5wKXA4uAC6vq1iRvByaraj3wfwJPBD6RBOBfqupk4HnA+5M8Sq8Ifv5uM+RLkvaBuVmSZtZqASPJScCf00vSF1TV+bv1n0UvUe84fO4vquqCkQYpSQtQVV0GXLZb21v6lk+cZrsvA0cPNzpJGk/mZknas9YKGANe6xrgY1V17sgDlCRJkiRJndHmHBgzXutakiRJkiQJ2i1gDHqt699IclOSS5IcNkU/Sc5prns9uW3btmHEKkmSJEmSWtT1q5B8GlhRVT8NXAFcNNUgLwclSZIkSdLC1mYBY8ZrXVfVfX3Xt74AeMGIYpMkSZIkSR3SZgFjkGtdL+9bPRm4fYTxSZIkSZKkjmjtKiQDXuv6D5OcDGwH7gfOaiteSZIkSZLUntYKGDDQta7fBLxp1HFJkiRJkqRu6foknpIkSZIkSRYwJEmSJElS91nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZpnFi9uOwJpVytWtB2BpLaNIg9YwJAkSZpnHnmk7QikXd19d9sRSGrbKPKABQxJGlNJTkpyR5KNSc6bov+AJB9r+q9NsqKv701N+x1JXjbSwCVpgTIvS9KeWcCQpDGUZBHwXuBXgaOAM5Ictduws4EHqurZwLuANc22RwGnA/8OOAn4y2Z/kkbh+OP5IU+A449vOxLNoXmdl1ev5g5WwurVI7tLSR0zojxgAUOSxtNxwMaququqfgRcDJyy25hTgIua5UuAlyRJ035xVT1UVd8ANjb7kzRsxx8P113H/myH666ziLGwzM+8vHo1rF3LSjbC2rUWMaRxNMI8YAFDksbTIcA3+9Y3N21Tjqmq7cB3gKcNuC1JzkkymWRy27Ztcxi6NF4WL4akd3vouq9SQICit76jz4k9572h52WYm9y8YsVjr8k71166y2vyzrWX7uxzYk9p4WorD1jAkCQNRVWtq6qJqppYunRp2+FI89b27VDVux1w3LE7vyCG3vqOvu3bWw5U88Jc5OZNmx57TR656rRdXpNHrjptZ9+mTXMXt6RuaSsPWMCQpPF0D3BY3/qhTduUY5IsBp4M3DfgtpKG4dpr4bjj+BGL4bjjeutaKOZnXl6zBlatYgPPhlWreuuSxssI84AFDEkaT9cDK5MckWR/epO/rd9tzHrgzGb5FcDnq6qa9tOb2fCPAFYC140obknXXsuBPGzxYuGZv3l5zRqewwaLF9I4G1Ee8GxJSRpDVbU9ybnA5cAi4MKqujXJ24HJqloPfAD4cJKNwP30vkzTjPs4cBuwHfiDqnqklQciSQuEeVmSZmYBQ5LGVFVdBly2W9tb+pZ/CLxymm3fCbxzqAFKmtYiL1y8IM3nvHz44W3ds6SuGEUe8BQSSZKkecYJO9U1TtgpaRR5wAKGJEmSJEnqPAsYkiRJkiSp8yxgSJIkSZKkzrOAIUmSJEmSOs8ChiRJkiRJ6jwLGJIkSZIkqfMsYEiSJEmSpM6zgCFJkiRJkjrPAoYkSZIkSeo8CxiSJEmSJKnzLGBIkiRJkqTOs4AhSZIkSZI6zwKGJEmSJEnqvMUzDUhyBPA6YEX/+Ko6eXhhSZIkSZIkPWbGAgbwd8AHgE8Djw41GknSTkl+DfhT4HB6+TpAVdWTWg1MksaYuVmS2jNIAeOHVfWeoUciSdrdu4HTgJurqlqORZLU827MzZLUikEKGH+e5K3APwAP7Wisqq8OLSpJEsA3gVv8gixJnWJulqSWDFLAOBr4T8AJPHYKSTXrkqThWQVcluQL7FpA/r/bC0mSxp65WZJaMkgB45XAM6vqR8MORpK0i3cC3wcOBPZvORZJUo+5WZJaMkgB4xbgKcC9ww1FkrSbp1fVT7UdhCRpF+ZmSWrJIAWMpwD/nOR6dj1MbkFdRjUBz2SUNKgR5YzLkvxKVf3D0O9JkjQoc7MktWSQAsZbh3XnSU4C/hxYBFxQVefv1n8A8CHgBcB9wKuratOw4pGkjvk94A1JHgIexkv1SVIXmJslqSWDFDB+rKo+29+Q5D8DX5jNHSdZBLwXeCmwGbg+yfqquq1v2NnAA1X17CSnA2uAV8/mfiVpvqiqn9i9LUnaiEWS1GNulqT27DfAmD9JsvOKI0lWAafMwX0fB2ysqruaCUIvnmK/pwAXNcuXAC8ZygfE8uU8QmD58jnftaQFaEQ5I8nbd1vfD/ibod6pJGmPzM2S1J5BChgnA/89yS8keSdwPHNTwDiE3nW0d9jctE05pqq2A98BnjYH9/2Y5cth61YCsHWrRQxJezbanHFYkjfBzlPqPglsGOYdSpJmZG6WpJbMWMCoqm/TK2K8F3g68IquXVI1yTlJJpNMbtu2bYDxj90e3bqVojl5kd56f78ktZgzfgc4uvmi/Gngqqp625zfiyRpb5ibJakl0xYwknwvyXeTfBfYCBwJvBLY0TZb9wCH9a0f2rRNOSbJYuDJ9Cbz3EVVrauqiaqaWLp06Yx3XPXYbb9ly3b+RyT01vv7JWnUOSPJsUmOBZ5Pb6LjV9P7de+LTfts9v3UJFck2dD8PWiKMcck+UqSW5PclOTVfX0fTPKNJDc2t2NmE48kzRfmZklq37STeE41QdEcux5YmeQIeoWK04Hf3G3MeuBM4CvAK4DPV81xWWHLFli+nNq6lSxb1luXpOmMJmf82W7rDwBHNe0FnPC4LQZ3HnBlVZ2f5LxmffVuY34AvKaqNiR5OnBDksur6sGm/41VdcksYpDmlcWLYfv2tqNQB5ibgRUrYNOmYd+LpIVirnPGIFchGYqq2p7kXOByepdRvbCqbm0mRpqsqvXAB4APJ9kI3E+vyDH3tmxhUaCsXUgaxJBzRlX98nD2DPTmMHpxs3wRcDW7fUmuqjv7lv81yb3AUuDBIcYlddYjj7QdgbrA3Nxz992jvDdJ891c54zWChgAVXUZcNlubW/pW/4hvdNWJElz4+CqnaWXrcDBexqc5Dhgf+Drfc3vTPIW4ErgvKp6aJptzwHOAXjGM54x27glaSEzN0vSAAa5CokkaR5J8rkkt0xx2+UKUs0pedOelpdkOfBh4Ler6tGm+U3Ac4EXAk/l8Yc49+9/r+Ynkjrp+OP5IU+A449vOxLNcwsiN69ezR2shNXT7l6SHjOEnNHqERhd4oSdkvZGl3NGVZ04XV+SbyVZXlVbmi/B904z7knAZ4A3V9U1ffve8QvhQ0n+GnjDHIYudcvxx8N117E/wHXX9davvbbtqDRPzfvcvHo1rF3LSoC1a3tta9bM+d1IWiCGlDP2+giMJLc3t3Nnfe+SpIElmWgmbpuNHZMj0/z91BT3sz/wSeBDu08I13yxJkmAU4FbZhmP1CmLFz92WeSHrvvqLpdNfui6r+7sW+xPQGos5Ny8YsVj74c71166y/vhzrWX7uxbsWIu7k3SfDeKnLHXBYyqeh7w74Fv7PvdSpL2weuAzyT52Cz2cT7w0iQbgBOb9R1fwC9oxrwK+EXgrCkuyfe3SW4GbgaWAO+YRSxS52zf/thlkQ847thdLpt8wHHH7uzzqiTqs2Bz86ZNj70fjlx12i7vhyNXnbazz6uSSILR5IzMdFXSJK8D/qaqHtj3uxmdiYmJmpycbDsMSZpSkhuqamKW+/iJqvreXMU0CuZmzVvHH89D132VA4471tNHFqi5yMvNfhZ+bl69mjvXXsqRq07z9BFJM5tFzpguNw9yBMbBwPVJPp7kpObQNEnSkCW5MsnLd2tbN9++IEvz2rXXciAPW7zQTmOdm9es4TlssHghaTBDyBkzFjCq6o+BlcAHgLOADUn+e5JnzVkUkqSpHAGsTvLWvrZZ/0ooSZoVc7MktWSgOTCayzltbW7bgYOAS5KsHWJskjTuHgReAhyc5NNJntxyPNJYWrSo7QjUMQ8yxrn58MPbjkDSfDLXOWPGObST/BfgNcC3gQuAN1bVw0n2AzYAq+Y2JElSI1W1Hfj9JGcB/4teAVnSCDlhp3Yz1rnZCTsl7Y25zhmDXATsqcBpVXV3f2NVPZrk1+Y2HElSn/ftWKiqDzYzzP9Bi/FIkszNktSaGQsYVfXWPfTdPrfhSJJ2qKr377Z+A/A7LYUjScLcLEltGmgODEmSJEmSpDZZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5LGTJKnJrkiyYbm70HTjHskyY3NbX1f+xFJrk2yMcnHkuw/uuglaeExL0vSYCxgSNL4OQ+4sqpWAlc261P5t6o6prmd3Ne+BnhXVT0beAA4e7jhStKCZ16WpAFYwJCk8XMKcFGzfBFw6qAbJglwAnDJvmwvSZqSeVmSBmABQ5LGz8FVtaVZ3gocPM24A5NMJrkmyalN29OAB6tqe7O+GThkqo2TnNNsP7lt27a5il2SFqKR5GUwN0ua3xa3HYAkae4l+RywbIquN/evVFUlqWl2c3hV3ZPkmcDnk9wMfGfQGKpqHbAOYGJiYrr7kKSx0IW83Ozf3Cxp3rKAIUkLUFWdOF1fkm8lWV5VW5IsB+6dZh/3NH/vSnI18HzgfwJPSbK4+bXvUOCeOX8AkrTAmJclafY8hUSSxs964Mxm+UzgU7sPSHJQkgOa5SXAi4DbqqqAq4BX7Gl7SdJeMS9L0gAsYEjS+DkfeGmSDcCJzTpJJpJc0Ix5HjCZ5Gv0vhifX1W3NX2rgdcn2Ujv3OsPjDR6SVp4zMuSNABPIZGkMVNV9wEvmaJ9Enhts/xl4Ohptr8LOG6YMUrSODEvS9JgPAJDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5rRQwkjw1yRVJNjR/D5pm3CNJbmxu60cdpyRJkiRJ6oa2jsA4D7iyqlYCVzbrU/m3qjqmuZ08uvAkSZIkSVKXtFXAOAW4qFm+CDi1pTgkSZIkSdI80FYB4+Cq2tIsbwUOnmbcgUkmk1yT5NTpdpbknGbc5LZt2+Y6VkmSJEmS1LLFw9pxks8By6boenP/SlVVkppmN4dX1T1Jngl8PsnNVfX13QdV1TpgHcDExMR0+5IkSZIkSfPU0AoYVXXidH1JvpVkeVVtSbIcuHeafdzT/L0rydXA84HHFTAkSZIkSdLC1tYpJOuBM5vlM4FP7T4gyUFJDmiWlwAvAm4bWYSSJEmSJKkz2ipgnA+8NMkG4MRmnSQTSS5oxjwPmEzyNeAq4PyqsoAhSZIkSdIYGtopJHtSVfcBL5mifRJ4bbP8ZeDoEYcmSZIkSZI6qK0jMCRJkiRJkgZmAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkKQxkuSpSa5IsqH5e9AUY345yY19tx8mObXp+2CSb/T1HTPqxyBJC425WZIGYwFDksbLecCVVbUSuLJZ30VVXVVVx1TVMcAJwA+Af+gb8sYd/VV14whilqSFztwsSQOwgCFJ4+UU4KJm+SLg1BnGvwL4bFX9YJhBSdKYMzdL0gAsYEjSeDm4qrY0y1uBg2cYfzrw0d3a3pnkpiTvSnLAdBsmOSfJZJLJbdu2zSJkSVrwzM2SNAALGJK0wCT5XJJbprid0j+uqgqoPexnOXA0cHlf85uA5wIvBJ4KrJ5u+6paV1UTVTWxdOnS2TwkSZr3zM2SNHuL2w5AkjS3qurE6fqSfCvJ8qra0nwJvncPu3oV8Mmqerhv3zt+IXwoyV8Db5iToCVpgTM3S9LseQSGJI2X9cCZzfKZwKf2MPYMdjtEufliTZLQO0f7lrkPUZLGjrlZkgZgAUOSxsv5wEuTbABObNZJMpHkgh2DkqwADgO+sNv2f5vkZuBmYAnwjlEELUkLnLlZkgbgKSSSNEaq6j7gJVO0TwKv7VvfBBwyxbgThhmfJI0jc7MkDcYjMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkacwkeWWSW5M8mmRiD+NOSnJHko1JzutrPyLJtU37x5LsP5rIJWlhMi9L0mAsYEjS+LkFOA344nQDkiwC3gv8KnAUcEaSo5ruNcC7qurZwAPA2cMNV5IWPPOyJA3AAoYkjZmqur2q7phh2HHAxqq6q6p+BFwMnJIkwAnAJc24i4BThxasJI0B87IkDcYChiRpKocA3+xb39y0PQ14sKq279b+OEnOSTKZZHLbtm1DDVaSxsCs8zKYmyXNb4vbDkCSNPeSfA5YNkXXm6vqU6OIoarWAesAJiYmahT3KUld1YW8DOZmSfObBQxJWoCq6sRZ7uIe4LC+9UObtvuApyRZ3Pzat6NdkrQH5mVJmr1WTiGZ7UzLkqShux5Y2cxsvz9wOrC+qgq4CnhFM+5MYGS/HErSGDMvSxp7bc2BMduZliVJ+yjJryfZDPwc8JkklzftT09yGUDzK965wOXA7cDHq+rWZhergdcn2Ujv3OsPjPoxSNJCYl6WpMG0cgpJVd0O0Js0eVo7Z1puxl4MnALcNvQAJWkBq6pPAp+cov1fgZf3rV8GXDbFuLvo5WhJ0hwwL0vSYLp8FZLpZlp+HGdTliRJkiRpYRvaERijnGnZ2ZQlSZIkSVrYhlbAGOJMy5IkSZIkacx0+RSSKWdabjkmSZIkSZLUgrYuozrbmZYljdCKFW1HIEmSJGnctXUVklnNtCxptO6+u+0IJEmSJI27Lp9CIkmSJEmSBFjAkDST1au5g5WwenXbkUiSJEkaY62cQiJpnli9GtauZSXA2rW9tjVr2oxIkiRJ0pjyCAypQ7owWeaKFZD0bneuvZQCAhS99R19XYhVC1vSdgSSJGlQfm5rFCxgSB3ShckyN22Cqt7tyFWn7SxehN76jr5Nm1oNU5IkSdKY8RQSSdNrThfZsPZSjlx1mqePSJIkSWqNR2BIXdHVyTLXrOE5bLB4odFavpxHCCxf3nYkkiRpJn5ua0QsYEhdsHOyzI29yTK7VsSQRmn5cti6lQBs3eqXIUmSuszPbY2QBQypJfNpsszDD287Ai10O17vCTy6desu74dHt27dpV+SJLXLz221xQKG1JL5NFlmF2LQwrbj9V4F+y1btsv7Yb9ly3bplyRJ7fJzW22xgCF1wZo1sGoVG3g2rFrlfBMab1u2wLJlFMCyZb11SZLUTX5ua4QsYEhd4WSZ0mO2bGER5ZcgSZLmAz+3NSIWMCRJkiRJUudZwJA6xMkyJUmSJGlqFjCkDnGyTOkxTvwlSdL84ee2RsEChiRJkiRJ6jwLGJIkSZIkqfMsYEiSJEmSpM6zgCFJkiRJkjrPAoYkSZIkSeo8CxiSJEmSJKnzLGBIkiRJkqTOSy2wC/Ym2QbcPce7XQJ8e473ua+MZWrG8nhdiQOMpd/hVbW0xftvxQLPzV2JA4xlOsYyNWPpGcu8DEPLzdPp0uttFMbt8cL4PWYf73BNmZsXXAFjGJJMVtVE23GAsUzHWLobBxiLhqMrz2VX4gBjmY6xTM1YNErj9hyP2+OF8XvMPt52eAqJJEmSJEnqPAsYkiRJkiSp8yxgDGZd2wH0MZapGcvjdSUOMBYNR1eey67EAcYyHWOZmrFolMbtOR63xwvj95h9vC1wDgxJkiRJktR5HoEhSZIkSZI6zwKGJEmSJEnqvLErYCS5MMm9SW6ZYdwLk2xP8oq+trVJbk1ye5L3JEnT/s4k30zy/TZjSfJjST6T5J+bvvPbiqVp//skX2v63pdkUVux9PWvn2m/w4wjydVJ7khyY3P7yRZj2T/JuiR3Nq+Z32gjliQ/0ffvcWOSbyd5dxuxNO1nJLk5yU3Na3jJILFo3w3peTQvL8C8PKxYYm42N2vOJXll81w+mmTaSy8mOal5/21Mct4oY5xLSZ6a5IokG5q/B00z7pG+99X6Ucc5WzM9X0kOSPKxpv/aJCtaCHNODfCYz0qyre95fW0bcc6FmfJ38/nwnubf4qYkx446RqpqrG7ALwLHArfsYcwi4PPAZcArmrafB/6x6VsEfAV4cdP3s8By4PttxgL8GPDLzZj9gS8Bv9riv8uTmr8B/idweluxNP2nAR/Z035H8G9yNTDRkdftfwPe0SzvByxp8/np2/YG4Bdbeg8tBu7d8W8BrAXetrfPl7dOvL7NywswLw/x3+VqzM3mZm9zegOeBzxnT++v5nn+OvBMejnya8BRbce+j493LXBes3wesGaacXv1udSl2yDPF/D7wPua5dOBj7Ud9wge81nAX7Qd6xw93j3mb+DlwGfpfY/4WeDaUcc4dkdgVNUXgftnGPY6el/s7u3fFDiQ3gv3AOAJwLeafV5TVVvajqWqflBVVzX7/hHwVeDQNmJp9vndZszipn+gGWOHEUuSJwKvB94xSAzDimNfDSmW3wH+R7P/R6vq2y3GAkCSI4GfpPefvDZiSXP78eZXvycB/zpILNp35uXRxNLsc17n5WHFsq/MzSOLxdw8D1XV7VV1xwzDjgM2VtVdTY68GDhl+NENxSnARc3yRcCp7YUyNIM8X/3/DpcAL9lxJNU8tZBeozMaIH+fAnyoeq4BnpJk+Wii6xm7AsZMkhwC/DrwV/3tVfUV4CpgS3O7vKpu72osSZ4C/O/AlW3GkuRyel9OvkcvibUVy58Cfwb8YC5imEUcAH/dHF72J3OV0Pc2lub1AfCnSb6a5BNJDm4jlt0231Gpn5PLI+1tLFX1MPB7wM30vhwfBXxgLmLRvjMvz20sCzkvzyIWMDebm9WGQ4Bv9q1vbtrmo4P7CudbgeneuwcmmUxyTZJTRxPanBnk+do5pqq2A98BnjaS6IZj0NfobzSnVFyS5LDRhNaK1t+zFjAe793A6qp6tL8xybPpHQp3KL0n6YQkv9DFWJIsBj4KvKeq7mozlqp6Gb3DuA8ATmgjliTHAM+qqk/O0f3vUxxN929V1dHALzS3/9RSLIubti9X1bH0DtP9v1qKpd/p9F67c2WvYknyBHpfkp8PPB24CXjTHMajffNuzMtzFssCz8t7HUvTbW42N2sfJPlcklumuC3IX6gHfbxNsW+6gt/hVTUB/Cbw7iTPGnbcGrpPAyuq6qeBK3jsCBQNweK2A+igCeDi5seXJcDLk2wHVgLXVNX3AZJ8Fvg5BjyccsSxrAM2VNW7OxALVfXDJJ+id8jRFS3E8j1gIskmeq/5n0xydVW9eMRxfKmq7gGoqu8l+Qi9w9I+NMs49iWW/0XvV89Lm+0/AZw9B3HsSyxfatZ/BlhcVTfMURz7EssPAarq6037x+mdx6p2mZfnNpaFnJf3JRZzs7lZ+6iqTpzlLu4B+n+tPrRp66Q9Pd4k30qyvKq2pHdI/b1TjevLN3cluZpeYe7rw4h3CAZ5vnaM2dwU758M3Dea8IZixsdcVf2P7wJ686EsVK2/Zz0CYzdVdURVraiqFfQOrf39qvo74F+AX0qyuPkl4JeAoR6qvC+xJHkHvUTxR23GkuSJTfLe8cvjfwD+uY1Yquqvqurpzfh/D9w5F1+S9+HfZHGaWdOb9l8DBp55fy5jaX4Z+DS9ydEAXgLc1kYsfZuewdz+wrcvsdwDHJVkabOLlzLk97lmZl6em1jGIS/vSyzmZnOzWnU9sDLJEUn2p3e0z7y7MkdjPXBms3wm8KndByQ5KMkBzfIS4EXM0Xt8RAZ5vvr/HV4BfL7JbfPVjI85u84BcTILOz+tB16Tnp8FvlP7MOfYbIzdERhJPkrvi8GSJJuBt9KbJIqqet8eNr2E3qG2N9M7JOzvq+rTzT7X0jsM7MeafV5QVW8bdSxJDgXeTO8L6VebXzb+oqouaCGWg4H1TZLej945rXvaz9BiGeQ+RxFHkh8HLm++lC0CPgf8v23E0vStBj6c3mXxtgG/3WIsAK+iN7PxwIb0fv5vwBeTPAzcTW9maQ2ReXlkscz7vDyMWMzN5mYNR5JfB/4fYCnwmSQ3VtXLkjydXk5+eVVtT3IucDm999+FVXVri2HPxvnAx5OcTe81+iqA9C4h+5+r6rX0TpF6f5JH6eXh86tq3hQwpnu+krwdmKyq9fTmp/lwko30JoM8vb2IZ2/Ax/yHSU4GttN7zGe1FvAsDZC/L6P3mbCR3hGDA31GzWmM87sgJkmSJEmSxoGnkEiSJEmSpM6zgCFJkiRJkjrPAoYkSZIkSeo8CxiSJEmSJKnzLGBI0gCSXJjk3iQzXl4xyTOSXJXkn5LclGSvZvCXJM3MvCxJ3bM3ubkZ/6oktyW5NclHZhpvAUPagyTfbzsGdcYHgZMGHPvHwMer6vn0Lh/2l8MKSho35mX1+SDmZakTzM3q80EGzM1JVgJvAl5UVf8O+KOZtrGAIUkDqKov0ru2905JnpXk75PckORLSZ67YzjwpGb5ycC/jjBUSRoL5mVJ6p69zM2/C7y3qh5otr13pv1bwNDYSPLC5rDRA5P8eHOY0k8NuO2SJF9J8h+SvDjJF5J8KsldSc5P8ltJrktyc5JnDfuxqDPWAa+rqhcAb+CxX/TeBvzHJJuBy4DXtROe1G3mZQ2BeVmapSRvT/JHfevvTPJfBtzW3KypTJebjwSOTPKPSa5JMuORG4uHGKTUKVV1fZL1wDuA/w34m6oa5LzZg4H1wB9X1RVJXgz8DPA8etXFu4ALquq4Jrm/jgEOf9L8luSJwM8Dn0iyo/mA5u8ZwAer6s+S/Bzw4SQ/VVWPthCq1FnmZc0l87I0Zy4ELgXenWQ/eqddHTfTRuZmTWWG3LwYWAm8GDgU+GKSo6vqwen2ZwFD4+btwPXAD4E/HGD8E4ArgT+oqi/0tV9fVVsAknwd+Iem/Wbgl+cuXHXYfsCDVXXMFH1n05z7V1VfSXIgsASY8bA4aQyZlzVXzMvSHKiqTUnuS/J84GDgn6rqvhk2MzdrOnvKzZuBa6vqYeAbSe6kV9C4fk87k8bJ04AnAj8BHDjA+O3ADcDLdmt/qG/50b71R7EwOBaq6rv0Eu0rAdLzM033vwAvadqfR++1tq2VQKXuMy9rTpiXpTl1AXAW8Nv0jsiYiblZU5ohN/8dvaMvSLKE3ikld+1pfxYwNG7eD/wJ8LfAmh2NSf55mvEF/A7w3CSrhx+euirJR4GvAM9JsjnJ2cBvAWcn+RpwK3BKM/y/Ar/btH8UOKuqqo24pXnAvKx9Yl6WhuqT9I5aeiFw+Y5Gc7Nmspe5+XLgviS3AVcBb5zpaB+rXhobSV4DPFxVH0myCPhykhOAm4BMt11VPZLkDGB9ku8Bt40mYnVJVZ0xTdfjJhuqqtuAFw03Imn+My9rNszL0vBU1Y+SXEXv0P9HYOcv5OZm7dFe5uYCXt/cBhKLzxp3SX4NeGZVvaftWCRJ5mVJalszeedXgVdW1Yamzdys1lnAkCRJkiQBkOQo4P8DPllV/7XteKR+FjAkSZIkSVLnOYmnJEmSJEnqPAsYkiRJkiSp8yxgSJIkSZKkzrOAIUmSJEmSOs8ChiRJkiRJ6rz/H5OAfsRyAUTYAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pltr = op.plotter.from_model(model)\n",
    "\n",
    "ax = pltr.plot_proj(evdf, color='b', ls='', marker='+', ms=10)\n",
    "pltr.plot_proj(evdf1, color='r', ls='', marker='.', ax=ax);"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "orbipy DOP853 time 12.63\n",
      "fastrk DOP8713M time 2.11\n",
      "speedup x 6.00\n"
     ]
    }
   ],
   "source": [
    "loops = 1000\n",
    "str0 = \"dtr.prop(s0, t0, t1)\"\n",
    "str1 = \"\"\"rk_prop_ev(crtbp, s0, 0, t1, *params.values(),\n",
    "                     values, terminals, directions, counts, accurates,\n",
    "                     call_event, 1e-12, 1e-12, 100, mc)\"\"\"\n",
    "\n",
    "r0 = timeit(str0, number=loops, globals=globals())\n",
    "r1 = timeit(str1, number=loops, globals=globals())\n",
    "\n",
    "print(f\"orbipy DOP853 time {r0:.2f}\")\n",
    "print(f\"fastrk DOP8713M time {r1:.2f}\")\n",
    "print(f\"speedup x {r0/r1:.2f}\")\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [],
   "source": [
    "@njit\n",
    "def event_0(t, s, *args):\n",
    "    return s[1]\n",
    "\n",
    "@njit\n",
    "def event_1(t, s, *args):\n",
    "    return s[2]\n",
    "\n",
    "@njit\n",
    "def event_2(t, s, *args):\n",
    "    return s[3]\n",
    "\n",
    "@njit\n",
    "def event_3(t, s, *args):\n",
    "    return s[4]\n",
    "\n",
    "@njit\n",
    "def event_4(t, s, *args):\n",
    "    return s[5]\n",
    "\n",
    "# use solve_ivp to calculate events\n",
    "sol = solve_ivp(crtbp, (t0, t1), s0,\n",
    "                method='DOP853', **params,\n",
    "                args=(mc,),\n",
    "                events=[event_0, event_1, event_2, event_3, event_4])"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solve_ivp DOP853 time 3.80\n",
      "fastrk DOP8713M time 0.22\n",
      "speedup x 17.65\n"
     ]
    }
   ],
   "source": [
    "loops = 100\n",
    "str0 = \"\"\"solve_ivp(crtbp, (t0, t1), s0,\n",
    "                method='DOP853', **params,\n",
    "                args=(mc,),\n",
    "                events=[event_0, event_1, event_2, event_3, event_4])\"\"\"\n",
    "str1 = \"\"\"rk_prop_ev(crtbp, s0, 0, t1, *params.values(),\n",
    "                     values, terminals, directions, counts, accurates,\n",
    "                     call_event, 1e-12, 1e-12, 100, mc)\"\"\"\n",
    "\n",
    "r0 = timeit(str0, number=loops, globals=globals())\n",
    "r1 = timeit(str1, number=loops, globals=globals())\n",
    "\n",
    "print(f\"solve_ivp DOP853 time {r0:.2f}\")\n",
    "print(f\"fastrk DOP8713M time {r1:.2f}\")\n",
    "print(f\"speedup x {r0/r1:.2f}\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "outputs": [],
   "source": [
    "dfs = []\n",
    "for i, evarr in enumerate(sol.y_events):\n",
    "    df = pd.DataFrame(evarr, columns=['x','y','z','vx','vy','vz'])\n",
    "    df['e'] = i\n",
    "    df['t'] = sol.t_events[i]\n",
    "    dfs.append(df)\n",
    "    \n",
    "evdf2 = pd.concat(dfs, axis=0)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "outputs": [
    {
     "data": {
      "text/plain": "0    -7.969181e-13\n1    -4.184764e-12\n2     5.084377e-12\n3     5.291101e-12\n4    -4.860801e-11\n5    -6.075274e-11\n6    -7.313705e-12\n7     1.243031e-09\n8    -4.274163e-10\n9    -1.897766e-10\n10   -4.538290e-10\n11   -1.785923e-10\n12    8.427596e-10\n13   -7.195746e-10\n14   -8.946062e-09\n15   -1.796951e-08\n16   -8.984693e-09\n17    3.321368e-07\n18   -1.491594e-07\n19   -6.215807e-08\n20   -1.675911e-07\n21   -5.982199e-08\n22    3.168344e-07\n23   -2.611800e-07\n24   -3.294949e-06\n25   -6.639791e-06\n26   -3.307739e-06\n27    1.223962e-04\n28   -5.498337e-05\n29   -2.290884e-05\nName: t, dtype: float64"
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# solve ivp calculated additional events at initial time\n",
    "# let's drop it and sort by time\n",
    "evdf3 = evdf2[evdf2.t > 1e-19].sort_values(by='t').reset_index(drop=True)\n",
    "\n",
    "# check event times\n",
    "evdf1.t - evdf3.t"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "outputs": [
    {
     "data": {
      "text/plain": "array([-8.60585474e-19,  2.32595247e-18,  9.75781955e-19, -1.55185442e-13,\n        8.34835673e-18, -2.84128339e-13])"
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evdf[evdf.e == 0].y.to_numpy() - evdf3[evdf3.e == 0].y.to_numpy()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 1080x360 with 3 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAFgCAYAAABNIolGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAu3UlEQVR4nO3df7xcdX3v+9ebROC2VkWTkghIUIPKKS3iFtp62lrESj29QKk/oD1HaLE8Tlvs6fVogg9b9Vg9h+TeXq23tpqLVLRVVA7WeMVSRFBPlR8bi/wsJGKooYlEfvjjYUUCn/vHrIRJ2Dt7kr1n1tp7Xs/HYx57re/3u9Z8JjPzmcln1vquVBWSJEmSJEldtl/bAUiSJEmSJM3EAoYkSZIkSeo8CxiSJEmSJKnzLGBIkiRJkqTOs4AhSZIkSZI6zwKGJEmSJEnqvAVbwEhyYZJ7k9wy4PhXJbktya1JPjLs+CRp3JiXJal7zM2S5pNUVdsxDEWSXwS+D3yoqn5qhrErgY8DJ1TVA0l+sqruHUWckjQuzMuS1D3mZknzyYI9AqOqvgjc39+W5FlJ/j7JDUm+lOS5TdfvAu+tqgeabU3EkjTHzMuS1D3mZknzyYItYExjHfC6qnoB8AbgL5v2I4Ejk/xjkmuSnNRahJI0XszLktQ95mZJnbS47QBGJckTgZ8HPpFkR/MBzd/FwErgxcChwBeTHF1VD444TEkaG+ZlSeoec7OkLhubAga9o00erKpjpujbDFxbVQ8D30hyJ73kfP0I45OkcWNelqTuMTdL6qyxOYWkqr5LL9G+EiA9P9N0/x29SjJJltA7PO6uFsKUpLFhXpak7jE3S+qyBVvASPJR4CvAc5JsTnI28FvA2Um+BtwKnNIMvxy4L8ltwFXAG6vqvjbilqSFyrwsSd1jbpY0nyzYy6hKkiRJkqSFY8EegSFJkiRJkhaOBTeJ55IlS2rFihVthyFJU7rhhhu+XVVL245j1MzNkrpqXPMymJslddd0uXnBFTBWrFjB5ORk22FI0pSS3N12DG0wN0vqqnHNy2BultRd0+VmTyGRJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJGkMJTkpyR1JNiY5b4r+dyW5sbndmeTBvr5H+vrWjzRwSVrAzM2StGcLbhJPSdKeJVkEvBd4KbAZuD7J+qq6bceYqvo/+sa/Dnh+3y7+raqOGVG4kjQWzM2SNDOPwJCk8XMcsLGq7qqqHwEXA6fsYfwZwEdHEpkkjS9zsyTNwAKGJI2fQ4Bv9q1vbtoeJ8nhwBHA5/uaD0wymeSaJKdOdydJzmnGTW7btm0OwpakBc3cLEkzsIAhdciKFW1HID3O6cAlVfVIX9vhVTUB/Cbw7iTPmmrDqlpXVRNVNbF06dJRxCppzCVtRzAy5mZ1zhi9/9QiCxhSh9x9d9sRaEzcAxzWt35o0zaV09ntEOWquqf5exdwNbuegy1J2jfmZkmagQUMSRo/1wMrkxyRZH96X4QfN2N9kucCBwFf6Ws7KMkBzfIS4EXAbbtvK0naa+ZmSZqBBQypK1av5g5WwurVbUeiBa6qtgPnApcDtwMfr6pbk7w9ycl9Q08HLq6q6mt7HjCZ5GvAVcD5/TPkS1Jrli/nEQLLl7cdyT4xN2tem+fvP80fXkZV6oLVq2HtWlYCrF3ba1uzps2ItMBV1WXAZbu1vWW39bdNsd2XgaOHGpwk7a3ly2HrVgKwdWtvfcuWtqPaa+ZmzUsL5P2n+cEjMKSWrFjRm+wogTvXXkoBAYre+o6+Lkzs2YUYJEnqt+NzMoFHt27d5XP00a1bd+mXNLd8/6ktFjCklmzaBFW925GrTtuZ9ENvfUffpk2thgk4uagkqXt2fE5WwX7Llu3yObrfsmW79EuaW77/1BYLGFIXrFkDq1axgWfDqlWePiJJ0t7YsgWWLaMAli3z8HVplHz/aYQsYEhdsWYNz2FD94oXTi4qSZoPtmxhEeV/nqQ2+P7TiDiJp6TpObmoJEmSpI7wCAypQw4/vO0I5tfkopIkSZLGhwUMqUO6MGHnfJpcVJKkfk4YKLXH959GwVNIJE2vOV1kw9pLOXLVaZ4+IkmSJKk1HoEhac+6OrmoJEmSpLFiAUOSJEmSJHWeBQxJM+rC5KKSJEmSxlurBYwkFya5N8kt0/QnyXuSbExyU5JjRx2jJCfslCRJktS+to/A+CBw0h76fxVY2dzOAf5qBDFJkiRJkqSOabWAUVVfBO7fw5BTgA9VzzXAU5IsH010kiRJkiSpK9o+AmMmhwDf7Fvf3LTtIsk5SSaTTG7btm1kwUmSJEmSpNHoegFjIFW1rqomqmpi6dKlbYcjSZIkSZLmWNcLGPcAh/WtH9q0SZIkSZKkMdL1AsZ64DXN1Uh+FvhOVW1pOyhJkiRJkjRai9u88yQfBV4MLEmyGXgr8ASAqnofcBnwcmAj8APgt9uJVJIkSZIktanVAkZVnTFDfwF/MKJwJEmSJElSR3X9FBJJkiRJkiQLGJIkSZIkqfssYEiSJEmSpM6zgCFJkiRJkjrPAoYkSZIkSeo8CxiSJEmSJKnzLGBIkiRJkqTOs4AhSWMoyUlJ7kiyMcl5U/SflWRbkhub22v7+s5MsqG5nTnayCVp4TI3S9KeLW47AEnSaCVZBLwXeCmwGbg+yfqqum23oR+rqnN32/apwFuBCaCAG5ptHxhB6JK0YJmbJWlmHoEhSePnOGBjVd1VVT8CLgZOGXDblwFXVNX9zRfjK4CThhSnJI0Tc7MkzcAChiSNn0OAb/atb27advcbSW5KckmSw/ZyW5Kck2QyyeS2bdvmIm5JWsjMzZI0AwsYkqSpfBpYUVU/Te+XvIv2dgdVta6qJqpqYunSpXMeoCSNIXOzpLFmAUOSxs89wGF964c2bTtV1X1V9VCzegHwgkG3lSTtE3OzJM3AAoYkjZ/rgZVJjkiyP3A6sL5/QJLlfasnA7c3y5cDv5LkoCQHAb/StEmSZsfcLEkz8CokkjRmqmp7knPpfbldBFxYVbcmeTswWVXrgT9McjKwHbgfOKvZ9v4kf0rvizbA26vq/pE/CElaYMzNkjSzVFXbMcypiYmJmpycbDsMSZpSkhuqaqLtOEbN3Cypq8Y1L4O5WVJ3TZebPYVEkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkKQxlOSkJHck2ZjkvCn6X5/ktiQ3JbkyyeF9fY8kubG5rR9t5JK0cJmbJWnPFrcdgCRptJIsAt4LvBTYDFyfZH1V3dY37J+Aiar6QZLfA9YCr276/q2qjhllzJK00JmbJWlmHoEhSePnOGBjVd1VVT8CLgZO6R9QVVdV1Q+a1WuAQ0ccoySNG3OzJM3AAoYkjZ9DgG/2rW9u2qZzNvDZvvUDk0wmuSbJqUOIT5LGkblZkmbgKSSSpGkl+Y/ABPBLfc2HV9U9SZ4JfD7JzVX19Sm2PQc4B+AZz3jGSOKVpHFgbpY0rjwCQ5LGzz3AYX3rhzZtu0hyIvBm4OSqemhHe1Xd0/y9C7gaeP5Ud1JV66pqoqomli5dOnfRS9LCZG6WpBlYwJCk8XM9sDLJEUn2B04HdpmxPsnzgffT+4J8b1/7QUkOaJaXAC8C+ieYkyTtG3OzJM3AU0gkacxU1fYk5wKXA4uAC6vq1iRvByaraj3wfwJPBD6RBOBfqupk4HnA+5M8Sq8Ifv5uM+RLkvaBuVmSZtZqASPJScCf00vSF1TV+bv1n0UvUe84fO4vquqCkQYpSQtQVV0GXLZb21v6lk+cZrsvA0cPNzpJGk/mZknas9YKGANe6xrgY1V17sgDlCRJkiRJndHmHBgzXutakiRJkiQJ2i1gDHqt699IclOSS5IcNkU/Sc5prns9uW3btmHEKkmSJEmSWtT1q5B8GlhRVT8NXAFcNNUgLwclSZIkSdLC1mYBY8ZrXVfVfX3Xt74AeMGIYpMkSZIkSR3SZgFjkGtdL+9bPRm4fYTxSZIkSZKkjmjtKiQDXuv6D5OcDGwH7gfOaiteSZIkSZLUntYKGDDQta7fBLxp1HFJkiRJkqRu6foknpIkSZIkSRYwJEmSJElS91nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZpnFi9uOwJpVytWtB2BpLaNIg9YwJAkSZpnHnmk7QikXd19d9sRSGrbKPKABQxJGlNJTkpyR5KNSc6bov+AJB9r+q9NsqKv701N+x1JXjbSwCVpgTIvS9KeWcCQpDGUZBHwXuBXgaOAM5Ictduws4EHqurZwLuANc22RwGnA/8OOAn4y2Z/kkbh+OP5IU+A449vOxLNoXmdl1ev5g5WwurVI7tLSR0zojxgAUOSxtNxwMaququqfgRcDJyy25hTgIua5UuAlyRJ035xVT1UVd8ANjb7kzRsxx8P113H/myH666ziLGwzM+8vHo1rF3LSjbC2rUWMaRxNMI8YAFDksbTIcA3+9Y3N21Tjqmq7cB3gKcNuC1JzkkymWRy27Ztcxi6NF4WL4akd3vouq9SQICit76jz4k9572h52WYm9y8YsVjr8k71166y2vyzrWX7uxzYk9p4WorD1jAkCQNRVWtq6qJqppYunRp2+FI89b27VDVux1w3LE7vyCG3vqOvu3bWw5U88Jc5OZNmx57TR656rRdXpNHrjptZ9+mTXMXt6RuaSsPWMCQpPF0D3BY3/qhTduUY5IsBp4M3DfgtpKG4dpr4bjj+BGL4bjjeutaKOZnXl6zBlatYgPPhlWreuuSxssI84AFDEkaT9cDK5MckWR/epO/rd9tzHrgzGb5FcDnq6qa9tOb2fCPAFYC140obknXXsuBPGzxYuGZv3l5zRqewwaLF9I4G1Ee8GxJSRpDVbU9ybnA5cAi4MKqujXJ24HJqloPfAD4cJKNwP30vkzTjPs4cBuwHfiDqnqklQciSQuEeVmSZmYBQ5LGVFVdBly2W9tb+pZ/CLxymm3fCbxzqAFKmtYiL1y8IM3nvHz44W3ds6SuGEUe8BQSSZKkecYJO9U1TtgpaRR5wAKGJEmSJEnqPAsYkiRJkiSp8yxgSJIkSZKkzrOAIUmSJEmSOs8ChiRJkiRJ6jwLGJIkSZIkqfMsYEiSJEmSpM6zgCFJkiRJkjrPAoYkSZIkSeo8CxiSJEmSJKnzLGBIkiRJkqTOs4AhSZIkSZI6zwKGJEmSJEnqvMUzDUhyBPA6YEX/+Ko6eXhhSZIkSZIkPWbGAgbwd8AHgE8Djw41GknSTkl+DfhT4HB6+TpAVdWTWg1MksaYuVmS2jNIAeOHVfWeoUciSdrdu4HTgJurqlqORZLU827MzZLUikEKGH+e5K3APwAP7Wisqq8OLSpJEsA3gVv8gixJnWJulqSWDFLAOBr4T8AJPHYKSTXrkqThWQVcluQL7FpA/r/bC0mSxp65WZJaMkgB45XAM6vqR8MORpK0i3cC3wcOBPZvORZJUo+5WZJaMkgB4xbgKcC9ww1FkrSbp1fVT7UdhCRpF+ZmSWrJIAWMpwD/nOR6dj1MbkFdRjUBz2SUNKgR5YzLkvxKVf3D0O9JkjQoc7MktWSQAsZbh3XnSU4C/hxYBFxQVefv1n8A8CHgBcB9wKuratOw4pGkjvk94A1JHgIexkv1SVIXmJslqSWDFDB+rKo+29+Q5D8DX5jNHSdZBLwXeCmwGbg+yfqquq1v2NnAA1X17CSnA2uAV8/mfiVpvqiqn9i9LUnaiEWS1GNulqT27DfAmD9JsvOKI0lWAafMwX0fB2ysqruaCUIvnmK/pwAXNcuXAC8ZygfE8uU8QmD58jnftaQFaEQ5I8nbd1vfD/ibod6pJGmPzM2S1J5BChgnA/89yS8keSdwPHNTwDiE3nW0d9jctE05pqq2A98BnjYH9/2Y5cth61YCsHWrRQxJezbanHFYkjfBzlPqPglsGOYdSpJmZG6WpJbMWMCoqm/TK2K8F3g68IquXVI1yTlJJpNMbtu2bYDxj90e3bqVojl5kd56f78ktZgzfgc4uvmi/Gngqqp625zfiyRpb5ibJakl0xYwknwvyXeTfBfYCBwJvBLY0TZb9wCH9a0f2rRNOSbJYuDJ9Cbz3EVVrauqiaqaWLp06Yx3XPXYbb9ly3b+RyT01vv7JWnUOSPJsUmOBZ5Pb6LjV9P7de+LTfts9v3UJFck2dD8PWiKMcck+UqSW5PclOTVfX0fTPKNJDc2t2NmE48kzRfmZklq37STeE41QdEcux5YmeQIeoWK04Hf3G3MeuBM4CvAK4DPV81xWWHLFli+nNq6lSxb1luXpOmMJmf82W7rDwBHNe0FnPC4LQZ3HnBlVZ2f5LxmffVuY34AvKaqNiR5OnBDksur6sGm/41VdcksYpDmlcWLYfv2tqNQB5ibgRUrYNOmYd+LpIVirnPGIFchGYqq2p7kXOByepdRvbCqbm0mRpqsqvXAB4APJ9kI3E+vyDH3tmxhUaCsXUgaxJBzRlX98nD2DPTmMHpxs3wRcDW7fUmuqjv7lv81yb3AUuDBIcYlddYjj7QdgbrA3Nxz992jvDdJ891c54zWChgAVXUZcNlubW/pW/4hvdNWJElz4+CqnaWXrcDBexqc5Dhgf+Drfc3vTPIW4ErgvKp6aJptzwHOAXjGM54x27glaSEzN0vSAAa5CokkaR5J8rkkt0xx2+UKUs0pedOelpdkOfBh4Ler6tGm+U3Ac4EXAk/l8Yc49+9/r+Ynkjrp+OP5IU+A449vOxLNcwsiN69ezR2shNXT7l6SHjOEnNHqERhd4oSdkvZGl3NGVZ04XV+SbyVZXlVbmi/B904z7knAZ4A3V9U1ffve8QvhQ0n+GnjDHIYudcvxx8N117E/wHXX9davvbbtqDRPzfvcvHo1rF3LSoC1a3tta9bM+d1IWiCGlDP2+giMJLc3t3Nnfe+SpIElmWgmbpuNHZMj0/z91BT3sz/wSeBDu08I13yxJkmAU4FbZhmP1CmLFz92WeSHrvvqLpdNfui6r+7sW+xPQGos5Ny8YsVj74c71166y/vhzrWX7uxbsWIu7k3SfDeKnLHXBYyqeh7w74Fv7PvdSpL2weuAzyT52Cz2cT7w0iQbgBOb9R1fwC9oxrwK+EXgrCkuyfe3SW4GbgaWAO+YRSxS52zf/thlkQ847thdLpt8wHHH7uzzqiTqs2Bz86ZNj70fjlx12i7vhyNXnbazz6uSSILR5IzMdFXSJK8D/qaqHtj3uxmdiYmJmpycbDsMSZpSkhuqamKW+/iJqvreXMU0CuZmzVvHH89D132VA4471tNHFqi5yMvNfhZ+bl69mjvXXsqRq07z9BFJM5tFzpguNw9yBMbBwPVJPp7kpObQNEnSkCW5MsnLd2tbN9++IEvz2rXXciAPW7zQTmOdm9es4TlssHghaTBDyBkzFjCq6o+BlcAHgLOADUn+e5JnzVkUkqSpHAGsTvLWvrZZ/0ooSZoVc7MktWSgOTCayzltbW7bgYOAS5KsHWJskjTuHgReAhyc5NNJntxyPNJYWrSo7QjUMQ8yxrn58MPbjkDSfDLXOWPGObST/BfgNcC3gQuAN1bVw0n2AzYAq+Y2JElSI1W1Hfj9JGcB/4teAVnSCDlhp3Yz1rnZCTsl7Y25zhmDXATsqcBpVXV3f2NVPZrk1+Y2HElSn/ftWKiqDzYzzP9Bi/FIkszNktSaGQsYVfXWPfTdPrfhSJJ2qKr377Z+A/A7LYUjScLcLEltGmgODEmSJEmSpDZZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5LGTJKnJrkiyYbm70HTjHskyY3NbX1f+xFJrk2yMcnHkuw/uuglaeExL0vSYCxgSNL4OQ+4sqpWAlc261P5t6o6prmd3Ne+BnhXVT0beAA4e7jhStKCZ16WpAFYwJCk8XMKcFGzfBFw6qAbJglwAnDJvmwvSZqSeVmSBmABQ5LGz8FVtaVZ3gocPM24A5NMJrkmyalN29OAB6tqe7O+GThkqo2TnNNsP7lt27a5il2SFqKR5GUwN0ua3xa3HYAkae4l+RywbIquN/evVFUlqWl2c3hV3ZPkmcDnk9wMfGfQGKpqHbAOYGJiYrr7kKSx0IW83Ozf3Cxp3rKAIUkLUFWdOF1fkm8lWV5VW5IsB+6dZh/3NH/vSnI18HzgfwJPSbK4+bXvUOCeOX8AkrTAmJclafY8hUSSxs964Mxm+UzgU7sPSHJQkgOa5SXAi4DbqqqAq4BX7Gl7SdJeMS9L0gAsYEjS+DkfeGmSDcCJzTpJJpJc0Ix5HjCZ5Gv0vhifX1W3NX2rgdcn2Ujv3OsPjDR6SVp4zMuSNABPIZGkMVNV9wEvmaJ9Enhts/xl4Ohptr8LOG6YMUrSODEvS9JgPAJDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5rRQwkjw1yRVJNjR/D5pm3CNJbmxu60cdpyRJkiRJ6oa2jsA4D7iyqlYCVzbrU/m3qjqmuZ08uvAkSZIkSVKXtFXAOAW4qFm+CDi1pTgkSZIkSdI80FYB4+Cq2tIsbwUOnmbcgUkmk1yT5NTpdpbknGbc5LZt2+Y6VkmSJEmS1LLFw9pxks8By6boenP/SlVVkppmN4dX1T1Jngl8PsnNVfX13QdV1TpgHcDExMR0+5IkSZIkSfPU0AoYVXXidH1JvpVkeVVtSbIcuHeafdzT/L0rydXA84HHFTAkSZIkSdLC1tYpJOuBM5vlM4FP7T4gyUFJDmiWlwAvAm4bWYSSJEmSJKkz2ipgnA+8NMkG4MRmnSQTSS5oxjwPmEzyNeAq4PyqsoAhSZIkSdIYGtopJHtSVfcBL5mifRJ4bbP8ZeDoEYcmSZIkSZI6qK0jMCRJkiRJkgZmAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkKQxkuSpSa5IsqH5e9AUY345yY19tx8mObXp+2CSb/T1HTPqxyBJC425WZIGYwFDksbLecCVVbUSuLJZ30VVXVVVx1TVMcAJwA+Af+gb8sYd/VV14whilqSFztwsSQOwgCFJ4+UU4KJm+SLg1BnGvwL4bFX9YJhBSdKYMzdL0gAsYEjSeDm4qrY0y1uBg2cYfzrw0d3a3pnkpiTvSnLAdBsmOSfJZJLJbdu2zSJkSVrwzM2SNAALGJK0wCT5XJJbprid0j+uqgqoPexnOXA0cHlf85uA5wIvBJ4KrJ5u+6paV1UTVTWxdOnS2TwkSZr3zM2SNHuL2w5AkjS3qurE6fqSfCvJ8qra0nwJvncPu3oV8Mmqerhv3zt+IXwoyV8Db5iToCVpgTM3S9LseQSGJI2X9cCZzfKZwKf2MPYMdjtEufliTZLQO0f7lrkPUZLGjrlZkgZgAUOSxsv5wEuTbABObNZJMpHkgh2DkqwADgO+sNv2f5vkZuBmYAnwjlEELUkLnLlZkgbgKSSSNEaq6j7gJVO0TwKv7VvfBBwyxbgThhmfJI0jc7MkDcYjMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkSZIkSVLnWcCQJEmSJEmdZwFDkiRJkiR1ngUMSZIkSZLUeRYwJEmSJElS51nAkCRJkiRJnWcBQ5IkSZIkdZ4FDEmSJEmS1HkWMCRJkiRJUudZwJAkSZIkSZ1nAUOSJEmSJHWeBQxJkiRJktR5FjAkacwkeWWSW5M8mmRiD+NOSnJHko1JzutrPyLJtU37x5LsP5rIJWlhMi9L0mAsYEjS+LkFOA344nQDkiwC3gv8KnAUcEaSo5ruNcC7qurZwAPA2cMNV5IWPPOyJA3AAoYkjZmqur2q7phh2HHAxqq6q6p+BFwMnJIkwAnAJc24i4BThxasJI0B87IkDcYChiRpKocA3+xb39y0PQ14sKq279b+OEnOSTKZZHLbtm1DDVaSxsCs8zKYmyXNb4vbDkCSNPeSfA5YNkXXm6vqU6OIoarWAesAJiYmahT3KUld1YW8DOZmSfObBQxJWoCq6sRZ7uIe4LC+9UObtvuApyRZ3Pzat6NdkrQH5mVJmr1WTiGZ7UzLkqShux5Y2cxsvz9wOrC+qgq4CnhFM+5MYGS/HErSGDMvSxp7bc2BMduZliVJ+yjJryfZDPwc8JkklzftT09yGUDzK965wOXA7cDHq+rWZhergdcn2Ujv3OsPjPoxSNJCYl6WpMG0cgpJVd0O0Js0eVo7Z1puxl4MnALcNvQAJWkBq6pPAp+cov1fgZf3rV8GXDbFuLvo5WhJ0hwwL0vSYLp8FZLpZlp+HGdTliRJkiRpYRvaERijnGnZ2ZQlSZIkSVrYhlbAGOJMy5IkSZIkacx0+RSSKWdabjkmSZIkSZLUgrYuozrbmZYljdCKFW1HIEmSJGnctXUVklnNtCxptO6+u+0IJEmSJI27Lp9CIkmSJEmSBFjAkDST1au5g5WwenXbkUiSJEkaY62cQiJpnli9GtauZSXA2rW9tjVr2oxIkiRJ0pjyCAypQ7owWeaKFZD0bneuvZQCAhS99R19XYhVC1vSdgSSJGlQfm5rFCxgSB3ShckyN22Cqt7tyFWn7SxehN76jr5Nm1oNU5IkSdKY8RQSSdNrThfZsPZSjlx1mqePSJIkSWqNR2BIXdHVyTLXrOE5bLB4odFavpxHCCxf3nYkkiRpJn5ua0QsYEhdsHOyzI29yTK7VsSQRmn5cti6lQBs3eqXIUmSuszPbY2QBQypJfNpsszDD287Ai10O17vCTy6desu74dHt27dpV+SJLXLz221xQKG1JL5NFlmF2LQwrbj9V4F+y1btsv7Yb9ly3bplyRJ7fJzW22xgCF1wZo1sGoVG3g2rFrlfBMab1u2wLJlFMCyZb11SZLUTX5ua4QsYEhd4WSZ0mO2bGER5ZcgSZLmAz+3NSIWMCRJkiRJUudZwJA6xMkyJUmSJGlqFjCkDnGyTOkxTvwlSdL84ee2RsEChiRJkiRJ6jwLGJIkSZIkqfMsYEiSJEmSpM6zgCFJkiRJkjrPAoYkSZIkSeo8CxiSJEmSJKnzLGBIkiRJkqTOSy2wC/Ym2QbcPce7XQJ8e473ua+MZWrG8nhdiQOMpd/hVbW0xftvxQLPzV2JA4xlOsYyNWPpGcu8DEPLzdPp0uttFMbt8cL4PWYf73BNmZsXXAFjGJJMVtVE23GAsUzHWLobBxiLhqMrz2VX4gBjmY6xTM1YNErj9hyP2+OF8XvMPt52eAqJJEmSJEnqPAsYkiRJkiSp8yxgDGZd2wH0MZapGcvjdSUOMBYNR1eey67EAcYyHWOZmrFolMbtOR63xwvj95h9vC1wDgxJkiRJktR5HoEhSZIkSZI6zwKGJEmSJEnqvLErYCS5MMm9SW6ZYdwLk2xP8oq+trVJbk1ye5L3JEnT/s4k30zy/TZjSfJjST6T5J+bvvPbiqVp//skX2v63pdkUVux9PWvn2m/w4wjydVJ7khyY3P7yRZj2T/JuiR3Nq+Z32gjliQ/0ffvcWOSbyd5dxuxNO1nJLk5yU3Na3jJILFo3w3peTQvL8C8PKxYYm42N2vOJXll81w+mmTaSy8mOal5/21Mct4oY5xLSZ6a5IokG5q/B00z7pG+99X6Ucc5WzM9X0kOSPKxpv/aJCtaCHNODfCYz0qyre95fW0bcc6FmfJ38/nwnubf4qYkx446RqpqrG7ALwLHArfsYcwi4PPAZcArmrafB/6x6VsEfAV4cdP3s8By4PttxgL8GPDLzZj9gS8Bv9riv8uTmr8B/idweluxNP2nAR/Z035H8G9yNTDRkdftfwPe0SzvByxp8/np2/YG4Bdbeg8tBu7d8W8BrAXetrfPl7dOvL7NywswLw/x3+VqzM3mZm9zegOeBzxnT++v5nn+OvBMejnya8BRbce+j493LXBes3wesGaacXv1udSl2yDPF/D7wPua5dOBj7Ud9wge81nAX7Qd6xw93j3mb+DlwGfpfY/4WeDaUcc4dkdgVNUXgftnGPY6el/s7u3fFDiQ3gv3AOAJwLeafV5TVVvajqWqflBVVzX7/hHwVeDQNmJp9vndZszipn+gGWOHEUuSJwKvB94xSAzDimNfDSmW3wH+R7P/R6vq2y3GAkCSI4GfpPefvDZiSXP78eZXvycB/zpILNp35uXRxNLsc17n5WHFsq/MzSOLxdw8D1XV7VV1xwzDjgM2VtVdTY68GDhl+NENxSnARc3yRcCp7YUyNIM8X/3/DpcAL9lxJNU8tZBeozMaIH+fAnyoeq4BnpJk+Wii6xm7AsZMkhwC/DrwV/3tVfUV4CpgS3O7vKpu72osSZ4C/O/AlW3GkuRyel9OvkcvibUVy58Cfwb8YC5imEUcAH/dHF72J3OV0Pc2lub1AfCnSb6a5BNJDm4jlt0231Gpn5PLI+1tLFX1MPB7wM30vhwfBXxgLmLRvjMvz20sCzkvzyIWMDebm9WGQ4Bv9q1vbtrmo4P7CudbgeneuwcmmUxyTZJTRxPanBnk+do5pqq2A98BnjaS6IZj0NfobzSnVFyS5LDRhNaK1t+zFjAe793A6qp6tL8xybPpHQp3KL0n6YQkv9DFWJIsBj4KvKeq7mozlqp6Gb3DuA8ATmgjliTHAM+qqk/O0f3vUxxN929V1dHALzS3/9RSLIubti9X1bH0DtP9v1qKpd/p9F67c2WvYknyBHpfkp8PPB24CXjTHMajffNuzMtzFssCz8t7HUvTbW42N2sfJPlcklumuC3IX6gHfbxNsW+6gt/hVTUB/Cbw7iTPGnbcGrpPAyuq6qeBK3jsCBQNweK2A+igCeDi5seXJcDLk2wHVgLXVNX3AZJ8Fvg5BjyccsSxrAM2VNW7OxALVfXDJJ+id8jRFS3E8j1gIskmeq/5n0xydVW9eMRxfKmq7gGoqu8l+Qi9w9I+NMs49iWW/0XvV89Lm+0/AZw9B3HsSyxfatZ/BlhcVTfMURz7EssPAarq6037x+mdx6p2mZfnNpaFnJf3JRZzs7lZ+6iqTpzlLu4B+n+tPrRp66Q9Pd4k30qyvKq2pHdI/b1TjevLN3cluZpeYe7rw4h3CAZ5vnaM2dwU758M3Dea8IZixsdcVf2P7wJ686EsVK2/Zz0CYzdVdURVraiqFfQOrf39qvo74F+AX0qyuPkl4JeAoR6qvC+xJHkHvUTxR23GkuSJTfLe8cvjfwD+uY1Yquqvqurpzfh/D9w5F1+S9+HfZHGaWdOb9l8DBp55fy5jaX4Z+DS9ydEAXgLc1kYsfZuewdz+wrcvsdwDHJVkabOLlzLk97lmZl6em1jGIS/vSyzmZnOzWnU9sDLJEUn2p3e0z7y7MkdjPXBms3wm8KndByQ5KMkBzfIS4EXM0Xt8RAZ5vvr/HV4BfL7JbfPVjI85u84BcTILOz+tB16Tnp8FvlP7MOfYbIzdERhJPkrvi8GSJJuBt9KbJIqqet8eNr2E3qG2N9M7JOzvq+rTzT7X0jsM7MeafV5QVW8bdSxJDgXeTO8L6VebXzb+oqouaCGWg4H1TZLej945rXvaz9BiGeQ+RxFHkh8HLm++lC0CPgf8v23E0vStBj6c3mXxtgG/3WIsAK+iN7PxwIb0fv5vwBeTPAzcTW9maQ2ReXlkscz7vDyMWMzN5mYNR5JfB/4fYCnwmSQ3VtXLkjydXk5+eVVtT3IucDm999+FVXVri2HPxvnAx5OcTe81+iqA9C4h+5+r6rX0TpF6f5JH6eXh86tq3hQwpnu+krwdmKyq9fTmp/lwko30JoM8vb2IZ2/Ax/yHSU4GttN7zGe1FvAsDZC/L6P3mbCR3hGDA31GzWmM87sgJkmSJEmSxoGnkEiSJEmSpM6zgCFJkiRJkjrPAoYkSZIkSeo8CxiSJEmSJKnzLGBI0gCSXJjk3iQzXl4xyTOSXJXkn5LclGSvZvCXJM3MvCxJ3bM3ubkZ/6oktyW5NclHZhpvAUPagyTfbzsGdcYHgZMGHPvHwMer6vn0Lh/2l8MKSho35mX1+SDmZakTzM3q80EGzM1JVgJvAl5UVf8O+KOZtrGAIUkDqKov0ru2905JnpXk75PckORLSZ67YzjwpGb5ycC/jjBUSRoL5mVJ6p69zM2/C7y3qh5otr13pv1bwNDYSPLC5rDRA5P8eHOY0k8NuO2SJF9J8h+SvDjJF5J8KsldSc5P8ltJrktyc5JnDfuxqDPWAa+rqhcAb+CxX/TeBvzHJJuBy4DXtROe1G3mZQ2BeVmapSRvT/JHfevvTPJfBtzW3KypTJebjwSOTPKPSa5JMuORG4uHGKTUKVV1fZL1wDuA/w34m6oa5LzZg4H1wB9X1RVJXgz8DPA8etXFu4ALquq4Jrm/jgEOf9L8luSJwM8Dn0iyo/mA5u8ZwAer6s+S/Bzw4SQ/VVWPthCq1FnmZc0l87I0Zy4ELgXenWQ/eqddHTfTRuZmTWWG3LwYWAm8GDgU+GKSo6vqwen2ZwFD4+btwPXAD4E/HGD8E4ArgT+oqi/0tV9fVVsAknwd+Iem/Wbgl+cuXHXYfsCDVXXMFH1n05z7V1VfSXIgsASY8bA4aQyZlzVXzMvSHKiqTUnuS/J84GDgn6rqvhk2MzdrOnvKzZuBa6vqYeAbSe6kV9C4fk87k8bJ04AnAj8BHDjA+O3ADcDLdmt/qG/50b71R7EwOBaq6rv0Eu0rAdLzM033vwAvadqfR++1tq2VQKXuMy9rTpiXpTl1AXAW8Nv0jsiYiblZU5ohN/8dvaMvSLKE3ikld+1pfxYwNG7eD/wJ8LfAmh2NSf55mvEF/A7w3CSrhx+euirJR4GvAM9JsjnJ2cBvAWcn+RpwK3BKM/y/Ar/btH8UOKuqqo24pXnAvKx9Yl6WhuqT9I5aeiFw+Y5Gc7Nmspe5+XLgviS3AVcBb5zpaB+rXhobSV4DPFxVH0myCPhykhOAm4BMt11VPZLkDGB9ku8Bt40mYnVJVZ0xTdfjJhuqqtuAFw03Imn+My9rNszL0vBU1Y+SXEXv0P9HYOcv5OZm7dFe5uYCXt/cBhKLzxp3SX4NeGZVvaftWCRJ5mVJalszeedXgVdW1Yamzdys1lnAkCRJkiQBkOQo4P8DPllV/7XteKR+FjAkSZIkSVLnOYmnJEmSJEnqPAsYkiRJkiSp8yxgSJIkSZKkzrOAIUmSJEmSOs8ChiRJkiRJ6rz/H5OAfsRyAUTYAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = pltr.plot_proj(evdf, color='b', ls='', marker='+', ms=10)\n",
    "pltr.plot_proj(evdf3, color='r', ls='', marker='.', ax=ax);"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}