test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import heyoka as hy\n",
    "import numpy as np\n",
    "import time\n",
    "import os\n",
    "import pickle\n",
    "\n",
    "from scipy.optimize import least_squares, root_scalar\n",
    "from copy import copy, deepcopy\n",
    "\n",
    "from src.hamiltonian_operations import *\n",
    "from src.initial_guesses import *\n",
    "from src.diff_corr import *\n",
    "from src.miscellaneous import *\n",
    "from src.ta_construction import *\n",
    "\n",
    "\n",
    "from matplotlib.pylab import plt\n",
    "\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define input "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get periodic orbit initial guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 0.01215057\n",
    "# We find the postiion of the lagrangian points\n",
    "xl_list = find_lagrange_points(mu)\n",
    "xL1, xL2, xL3, (xL4, yL4), (xL5, yL5) = xl_list[1], xl_list[2], xl_list[0], (xl_list[3][0], xl_list[3][1]), (xl_list[4][0], xl_list[4][1])\n",
    "\n",
    "xL = xL1\n",
    "xL_str = 'xL1'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lyapunov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "use_real = False\n",
    "direction = False\n",
    "po_ic, po_period = get_po_ig(xL, mu, use_real=use_real, direction=direction, conj_mom=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Halo"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# [1.1406936568979191, 0.0, 0.11171730476145951, 0.0, 0.2290282993165571, 0.0]\n",
    "# 3.1544236951257316\n",
    "\n",
    "\n",
    "po_ic = [1.1406936568979191, 0.0, 0.11171730476145951, 0.0, 0.2290282993165571, 0.0]\n",
    "po_ic = convert_vel_to_conj_mom(np.array(po_ic)).tolist()\n",
    "po_period = 3.1544236951257316"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_cond_var = np.eye(6).reshape((36,)).tolist()\n",
    "ta = get_ta_var(po_ic, ic_var=init_cond_var, conj_mom=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsteps = 2000\n",
    "epochs = np.linspace(0, po_period, nsteps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ta.pars[0] = mu\n",
    "ta.time = 0\n",
    "ta.state[:] = po_ic + np.eye(6).reshape((36,)).tolist()\n",
    "out2 = ta.propagate_grid(epochs)\n",
    "out2 = out2[5]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1, projection='3d')\n",
    "\n",
    "#Plot the lagrangian points and primaries\n",
    "#fig_name = f'{xL_str}_ur{use_real}_dir{direction}'\n",
    "#ax.set_title(fig_name)\n",
    "buffer=1.5\n",
    "ax.set_xlim(xL1-buffer, xL1+buffer)\n",
    "ax.set_ylim(-buffer, +buffer)\n",
    "# ax.set_xlim(-buffer, buffer)\n",
    "# ax.set_ylim(-buffer, buffer)\n",
    "ax.set_zlim(-buffer, buffer)\n",
    "ax.scatter(-mu, 0, 0, c='k', s=100)\n",
    "ax.scatter((1-mu), 0, 0, c='k', s=30)\n",
    "ax.scatter(xL1, 0, 0, c='r', s=10)\n",
    "ax.scatter(xL2, 0, 0, c='r', s=10)\n",
    "ax.scatter(xL3, 0, 0, c='r', s=10)\n",
    "ax.scatter(xL4, yL4, 0, c='r', s=10)\n",
    "ax.scatter(xL5, yL5, 0, c='r', s=10)\n",
    "# # zero velocity curve\n",
    "#plot_zero_vel_curves(ax, mu, new_ic, xL, buffer=buffer)\n",
    "ax.plot(out2[:, 0], out2[:, 1], out2[:, 2], linewidth=2, c='b')\n",
    "ax.scatter(out2[-1, 0], out2[-1, 1], out2[-1, 2], s=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ta, ic, err = diff_corr_phasing(ta, po_ic, tol=1e-12, max_iter=100, conj_mom=True)\n",
    "t_final = ic[6]\n",
    "ic = ic[:6]\n",
    "# ic, t_final, ds, out = perform_predictor_corrector_step(ta, np.array(po_ic), po_period, ds=1e-4, verbose=True)\n",
    "# out, _ = find_po_family(ta, np.array(po_ic), po_period, ds=1e-6, max_iter=30, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "timesteps = np.linspace(0, t_final, 2000)\n",
    "ta.time = 0.\n",
    "ta.state[:] = ic.tolist() + init_cond_var\n",
    "out3 = ta.propagate_grid(timesteps)\n",
    "out3 = out3[5]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "out.values()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1, projection='3d')\n",
    "\n",
    "#Plot the lagrangian points and primaries\n",
    "#fig_name = f'{xL_str}_ur{use_real}_dir{direction}'\n",
    "#ax.set_title(fig_name)\n",
    "buffer=1.5\n",
    "ax.set_xlim(xL1-buffer, xL1+buffer)\n",
    "ax.set_ylim(-buffer, +buffer)\n",
    "# ax.set_xlim(-buffer, buffer)\n",
    "# ax.set_ylim(-buffer, buffer)\n",
    "ax.set_zlim(-buffer, buffer)\n",
    "ax.scatter(-mu, 0, 0, c='k', s=100)\n",
    "ax.scatter((1-mu), 0, 0, c='k', s=30)\n",
    "ax.scatter(xL1, 0, 0, c='r', s=10)\n",
    "ax.scatter(xL2, 0, 0, c='r', s=10)\n",
    "ax.scatter(xL3, 0, 0, c='r', s=10)\n",
    "ax.scatter(xL4, yL4, 0, c='r', s=10)\n",
    "ax.scatter(xL5, yL5, 0, c='r', s=10)\n",
    "# # zero velocity curve\n",
    "#plot_zero_vel_curves(ax, mu, new_ic, xL, buffer=buffer)\n",
    "ax.plot(out3[:, 0], out3[:, 1], out3[:, 2], linewidth=2, c='b')\n",
    "ax.scatter(out3[-1, 0], out3[-1, 1], out3[-1, 2], s=5)\n",
    "for i in out.keys():\n",
    "    ax.plot(out[i][:, 0], out[i][:, 1], out[i][:, 2], 'g', alpha=0.3) if out[i] is not None else None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_ic = ic\n",
    "new_period = t_final"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perform PA Continuation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ds: 0.0001 err: [3.11467433e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.00011000000000000002 err: [2.26261497e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.00012100000000000003 err: [8.62584349e-14]\n",
      "Converged - increase ds\n",
      "ds: 0.00013310000000000004 err: [1.24123333e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.00014641000000000006 err: [1.37482178e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.00016105100000000008 err: [8.02700609e-14]\n",
      "Converged - increase ds\n",
      "ds: 0.0001771561000000001 err: [3.09737523e-16]\n",
      "Converged - increase ds\n",
      "ds: 0.00019487171000000012 err: [5.37765391e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.00021435888100000014 err: [1.82241411e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.00023579476910000016 err: [5.94943801e-16]\n",
      "Converged - increase ds\n",
      "ds: 0.00025937424601000017 err: [1.80810246e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0002853116706110002 err: [5.89648407e-14]\n",
      "Converged - increase ds\n",
      "ds: 0.00031384283767210024 err: [1.81638475e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0003452271214393103 err: [1.98228973e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.00037974983358324136 err: [4.04619569e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0004177248169415655 err: [4.17924718e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0004594972986357221 err: [2.62755965e-16]\n",
      "Converged - increase ds\n",
      "ds: 0.0005054470284992944 err: [8.05469136e-16]\n",
      "Converged - increase ds\n",
      "ds: 0.000555991731349224 err: [1.2924179e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0006115909044841464 err: [8.08311939e-14]\n",
      "Converged - increase ds\n",
      "ds: 0.0006727499949325611 err: [2.64473814e-15]\n",
      "Converged - increase ds\n",
      "ds: 0.0007400249944258172 err: [4.64094311e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.000814027493868399 err: [4.14887822e-14]\n",
      "Converged - increase ds\n",
      "ds: 0.000895430243255239 err: [1.81675615e-16]\n",
      "Converged - increase ds\n",
      "ds: 0.000984973267580763 err: [2.97265585e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0010834705943388394 err: [8.25311758e-16]\n",
      "Converged - increase ds\n",
      "ds: 0.0011918176537727234 err: [3.66867611e-16]\n",
      "Converged - increase ds\n",
      "ds: 0.001310999419149996 err: [3.49667825e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0014420993610649957 err: [1.16230463e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0015863092971714955 err: [1.66634633e-14]\n",
      "Converged - increase ds\n",
      "ds: 0.0017449402268886452 err: [3.10523352e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.00191943424957751 err: [3.84503773e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.002111377674535261 err: [3.27421303e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0023225154419887874 err: [3.82761879e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0025547669861876662 err: [3.45514758e-13]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0026764225569585074 err: [0.11506972]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.002548973863770007 err: [4.59765311e-13]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.00254319387768436 err: [3.11321734e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.002797513265452796 err: [0.60265428]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.002664298348050282 err: [0.16470189]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.0025374269981431256 err: [3.77568908e-13]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0025316731954262476 err: [2.34968446e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0027848405149688725 err: [4.21382459e-16]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0024001949542115705 err: [2.30453287e-13]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.00239475233526778 err: [2.68501513e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.002634227568794558 err: [0.04599142]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.002508788160756722 err: [1.0055505e-13]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0026282542636498992 err: [0.00687956]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.0025030992987141896 err: [4.35268846e-16]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0026222945034148656 err: [0.02668149]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.002497423336585586 err: [1.78882338e-13]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.002491760245119406 err: [9.54320144e-14]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.002610415494886997 err: [0.02219074]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.002486109995130473 err: [7.23006501e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0027347209946435207 err: [0.05597475]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.0026044961853747815 err: [0.10408298]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.0024804725574997916 err: [1.7913758e-16]\n",
      "Converged - increase ds\n",
      "ds: 0.002728519813249771 err: [0.74755605]\n",
      "Low Precision - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.002474847903174395 err: [1.26853713e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0027223326934918348 err: [0.0334838]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.002592697803325557 err: [0.06731966]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.002469236003167197 err: [9.01055647e-16]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0024636368285568403 err: [3.61101182e-15]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0024580503504875503 err: [1.00120186e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0027038553855363054 err: [0.12749055]\n",
      "Low Precision - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0024524765401689845 err: [4.71437002e-13]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0024469153688760842 err: [3.38060602e-13]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0024413668079489273 err: [1.35229597e-13]\n",
      "Converged - increase ds\n",
      "Singular matrix - reducing ds\n",
      "ds: 0.0025576223702322098 err: [0.15881714]\n",
      "Low Precision - reducing ds\n",
      "ds: 0.0024358308287925805 err: [2.23333444e-13]\n",
      "Converged - increase ds\n",
      "ds: 0.0026794139116718386 err: [0.16925168]\n",
      "Low Precision - reducing ds\n",
      "Singular matrix - reducing ds\n"
     ]
    }
   ],
   "source": [
    "verbose = True\n",
    "start_time = time.time()\n",
    "max_iter = 100\n",
    "ds = 1e-4# We call the predictor\n",
    "# Main iterations\n",
    "out, _ = find_po_family(ta, new_ic, new_period, ds, max_iter=max_iter, corr_max_iter=100, verbose=verbose)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D Plot for Lyapunov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1563b5289ee94652a5955bbf711776e5",
       "version_major": 2,
       "version_minor": 0
      },
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxf0lEQVR4nO3deXTU1f3/8dckkAkRkrCELBKWhCWKQBAkhrYkSiQsUjhaFVwIFAF3EWSJX2X3i6BFrMUC9QuhfrUuFbdW4QsBDkdMowI5IEUURAKUhM1MCEvQ5P7+4DA/R5KY4CSTzH0+zvkcnDv3zuf9uZ3OvHJnPp9xGGOMAAAAYI0AXxcAAACAukUABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAANZJTU1Vamqqr8u4RPv27TV69Gj37U2bNsnhcGjTpk11Wsfo0aPVvn37Ot0ngLpFAAQsNGvWLDkcDh0/frzKfk8//bR++9vfKjIyUg6HQ7NmzaqbAivRvn17ORyOCrdz5875tLb6IDU1tdL5+fLLL31dHoB6pJGvCwBQfz355JOKiopSz549tXbtWl+XI0lKTEzU5MmTL2kPCgryQTW1q1+/fjp79myNjq1NmzaaP3/+Je0xMTHeLA1AA0cABFCp/fv3q3379jp+/LgiIiJ8UoMxRufOnVOTJk0kSVdeeaXuvvtun9RS1wICAhQcHPyz/c6cOaOQkBBJUlhYmDXzA+Dy8REw4EfOnj2rhIQEJSQk6OzZs+72kydPKjo6Wn379lVZWVm1H89b3wO7+JHzT2VlZcnhcOjbb7/12OfNN9+stWvXqnfv3mrSpImWLVtWrf2sXLlSN954o1q3bi2n06mrr75af/7zn6s19sUXX1TXrl0VEhKi5s2bq3fv3nrttdc8+hw+fFi///3vFRkZKafTqa5du2rFihXVevwfM8Zo3rx5atOmjUJCQnTDDTdo165dl/Sr6DuAqampuuaaa7R161b169dPISEheuKJJ6q13/fee09DhgxRTEyMnE6n4uPjNXfu3Go9J15//XX16tVLzZo1U2hoqLp166YXXnjBo09RUZEmTpyo2NhYOZ1OdezYUQsWLFB5eXm16gNQd1gBBPxIkyZNtGrVKv3qV7/Sf/3Xf2nRokWSpAcffFAul0tZWVkKDAz0cZU/b8+ePRo5cqQmTJigcePGqUuXLu77vv/++0u+uxgSEqKQkBD9+c9/VteuXfXb3/5WjRo10gcffKAHHnhA5eXlevDBByvd31/+8hc98sgj+t3vfqdHH31U586d044dO5Sbm6s777xTklRYWKjrr79eDodDDz30kCIiIvTRRx9p7NixKi4u1sSJE6t9fDNmzNC8efM0ePBgDR48WNu2bdOAAQN0/vz5ao0/ceKEBg0apBEjRujuu+9WZGSk+76ysrJL5ic4OFhNmzZVVlaWmjZtqkmTJqlp06basGGDZsyYoeLiYj377LOV7m/dunUaOXKk+vfvrwULFkiSdu/erS1btujRRx+VdGEVMiUlRYcPH9aECRPUtm1bffLJJ8rMzNSRI0e0ePHias8PgDpgAPidzMxMExAQYDZv3mzeeustI8ksXrzYff/MmTONJHPs2LFqPd6xY8eMJDNz5szLqufi/n5q5cqVRpLZv3+/u61du3ZGklmzZs0l/S/e99PtYl1nzpy5ZEx6erqJi4vzaEtJSTEpKSnu28OGDTNdu3at8hjGjh1roqOjzfHjxz3aR4wYYcLCwircd0WOHj1qgoKCzJAhQ0x5ebm7/YknnjCSTEZGhrtt48aNRpLZuHGjR+2SzNKlSy957Iv3/XS7+JgV1ThhwgQTEhJizp07527LyMgw7dq1c99+9NFHTWhoqPnhhx8qPa65c+eaK664wnz11Vce7dOnTzeBgYEmPz+/0rEA6h4fAQN+aNasWeratasyMjL0wAMPKCUlRY888oivy6q2Dh06KD09vcL7kpKStG7dOo9t1KhRkuT+nqAkuVwuHT9+XCkpKfrmm2/kcrkq3V94eLgOHTqkzz77rML7jTF6++23NXToUBljdPz4cfeWnp4ul8ulbdu2VevY1q9fr/Pnz+vhhx/2+Fi8JiuITqdTY8aMqfC+9u3bXzI/U6dOleQ5P6dOndLx48f1m9/8RmfOnKnyLOHw8HCdPn1a69atq7TPW2+9pd/85jdq3ry5x/ykpaWprKxMmzdvrvbxAah9fAQM+KGgoCCtWLFC1113nYKDg7Vy5coKv4NXX3Xo0KHS+1q1aqW0tLQK79uyZYtmzpypnJwcnTlzxuM+l8ulsLCwCsdNmzZN69evV58+fdSxY0cNGDBAd955p371q19Jko4dO6aioiItX75cy5cvr/Axjh49Wp1D04EDByRJnTp18miPiIhQ8+bNq/UYV155ZaVnBl9xxRWVzs+uXbv05JNPasOGDSouLva4r6qA/MADD+jNN9/UoEGDdOWVV2rAgAG6/fbbNXDgQHefr7/+Wjt27Kj0ZKHqzg+AukEABPzUxcu2nDt3Tl9//XWVoaq2VRY+Kzv54McrVdW1b98+9e/fXwkJCVq0aJFiY2MVFBSkDz/8UM8//3yVJyJcddVV2rNnj/7xj39ozZo1evvtt/XSSy9pxowZmj17tnvs3XffrYyMjAofo3v37jWu+XJdzvwUFRUpJSVFoaGhmjNnjuLj4xUcHKxt27Zp2rRpVc5P69atlZeXp7Vr1+qjjz7SRx99pJUrV2rUqFFatWqVJKm8vFw33XSTe7Xxpzp37lzjmgHUHgIg4Id27NihOXPmaMyYMcrLy9O9996rnTt3VroCVtsurmwVFRUpPDzc3X5xNcwbPvjgA5WWlur9999X27Zt3e0bN26s1vgrrrhCd9xxh+644w6dP39et9xyi55++mllZmYqIiJCzZo1U1lZWaWra9XVrl07SRdWzOLi4tztx44d03ffffeLHrsqmzZt0okTJ7R69Wr169fP3b5///5qjQ8KCtLQoUM1dOhQlZeX64EHHtCyZcv01FNPqWPHjoqPj1dJSckvnh8AdYPvAAJ+5vvvv9fo0aMVExOjF154QVlZWSosLNRjjz3ms5ri4+MlyeN7YKdPn3avHnnDxbObjTHuNpfLpZUrV/7s2BMnTnjcDgoK0tVXXy1jjL7//nsFBgbq1ltv1dtvv60vvvjikvHHjh2rdp1paWlq3LixXnzxRY9aa/ss2Yrm5/z583rppZd+duxP5ycgIMC94llaWipJuv3225WTk1PhBcOLior0ww8/XHbtALyPFUDAz8ybN095eXnKzs5Ws2bN1L17d82YMUNPPvmkfve732nw4MHuvosWLXJfQPiigIAA93XlXnnlFR04cMD9fbrNmzdr3rx5kqR77rnHvZr1cwYMGKC2bdtq7NixmjJligIDA7VixQpFREQoPz/fG4etAQMGuFepJkyYoJKSEv3lL39R69atdeTIkZ8dGxUVpV/96leKjIzU7t279ac//UlDhgxRs2bNJEnPPPOMNm7cqKSkJI0bN05XX321Tp48qW3btmn9+vU6efJkteqMiIjQ448/rvnz5+vmm2/W4MGDtX37dn300Udq1arVL56HyvTt21fNmzdXRkaGHnnkETkcDr3yyisegbAy9957r06ePKkbb7xRbdq00YEDB/Tiiy8qMTFRV111lSRpypQpev/993XzzTdr9OjR6tWrl06fPq2dO3fq73//u7799ttaPT4ANeTDM5ABeNnWrVtNo0aNzMMPP+zR/sMPP5jrrrvOxMTEmO+++859WZaKtsDAQPe4yi4rop9cmqS6tSUlJZmgoCDTtm1bs2jRokovAzNkyJAKH6Oq+4wx5v333zfdu3c3wcHBpn379mbBggVmxYoVl+zjp5eBWbZsmenXr59p2bKlcTqdJj4+3kyZMsW4XC6Pxy8sLDQPPvigiY2NNY0bNzZRUVGmf//+Zvny5TWai7KyMjN79mwTHR1tmjRpYlJTU80XX3xh2rVrV63LwFR2yZqq7jPGmC1btpjrr7/eNGnSxMTExJipU6eatWvXXrKPn14G5u9//7sZMGCAad26tft/vwkTJpgjR454PP6pU6dMZmam6dixowkKCjKtWrUyffv2Nc8995w5f/58jeYIQO1yGFONP/8AAADgN/gOIAAAgGX4DiCAy3b27Nkqrx8nSS1atKj0mnX+5tixY1X+rm5QUJBatGhRhxUBQMX8ZgVw8+bNGjp0qGJiYuRwOPTuu+9W2f/ij6z/dCsoKKibggE/8MYbbyg6OrrK7ZNPPvF1mXXmuuuuq3IubrnlFl+XCACS/GgF8PTp0+rRo4d+//vf1+hFds+ePQoNDXXfbt26dW2UB/il9PT0Kn8eTJJ69OhRR9X43quvvqqzZ89Wen91f+kDAGqb3wTAQYMGadCgQTUe17p1a48L0wKovosrW7jg4k/HAUB95zcfAV+uxMRERUdH66abbtKWLVt8XQ4AAECt85sVwJqKjo7W0qVL1bt3b5WWlurll19WamqqcnNzde2111Y4prS01H3Ve+nCb1+ePHlSLVu2rPS3TgEAQP1ijNGpU6cUExOjgAA718L88jqADodD77zzjoYPH16jcSkpKWrbtq1eeeWVCu+fNWuWZs+e7YUKAQCArx08eFBt2rTxdRk+Ye0KYEX69Omjjz/+uNL7MzMzNWnSJPdtl8ultm3b6uDBgx4nkgAAgPqruLhYsbGx7p96tBEB8Efy8vKq/EK70+mU0+m8pD00NJQACABAA2Pz17f8JgCWlJRo79697tv79+9XXl6eWrRoobZt2yozM1OHDx/WX//6V0nS4sWL1aFDB3Xt2lXnzp3Tyy+/rA0bNuj//u//fHUIAAAAdcJvAuDnn3+uG264wX374ke1GRkZysrK0pEjR5Sfn+++//z585o8ebIOHz6skJAQde/eXevXr/d4DAAAAH/klyeB1JXi4mKFhYXJ5XLxETAAAA0E799cBxAAAMA6BEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAs4zcBcPPmzRo6dKhiYmLkcDj07rvv/uyYTZs26dprr5XT6VTHjh2VlZVV63UCAAD4mt8EwNOnT6tHjx5asmRJtfrv379fQ4YM0Q033KC8vDxNnDhR9957r9auXVvLlQIAAPhWI18X4C2DBg3SoEGDqt1/6dKl6tChg/7whz9Ikq666ip9/PHHev7555Wenl5bZQIAAPic36wA1lROTo7S0tI82tLT05WTk1PpmNLSUhUXF3tsAAAADY21AbCgoECRkZEebZGRkSouLtbZs2crHDN//nyFhYW5t9jY2LooFQAAwKusDYCXIzMzUy6Xy70dPHjQ1yUBAADUmN98B7CmoqKiVFhY6NFWWFio0NBQNWnSpMIxTqdTTqezLsoDAACoNdauACYnJys7O9ujbd26dUpOTvZRRQAAAHXDbwJgSUmJ8vLylJeXJ+nCZV7y8vKUn58v6cLHt6NGjXL3v++++/TNN99o6tSp+vLLL/XSSy/pzTff1GOPPeaL8gEAAOqM3wTAzz//XD179lTPnj0lSZMmTVLPnj01Y8YMSdKRI0fcYVCSOnTooH/+859at26devTooT/84Q96+eWXuQQMAADwew5jjPF1EQ1VcXGxwsLC5HK5FBoa6utyAABANfD+7UcrgAAAAKgeAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAlrH2t4ABALUsN1f66iupc2cpKcnX1QD4EVYAAQDeN22adP310qhRF/6dNs3XFQH4EQIgAMC7cnOlhQs92xYuvNAOoF4gAAIAvOurr2rWDqDOEQABAN7VuXPN2gHUOQIgAMC7kpKkqVM926ZN40QQoB7hLGAAgPctWCDdcgtnAQP1FAEQAFA7kpIIfkA9xUfAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGb8KgEuWLFH79u0VHByspKQkffrpp5X2zcrKksPh8NiCg4PrsFoAAADf8JsA+MYbb2jSpEmaOXOmtm3bph49eig9PV1Hjx6tdExoaKiOHDni3g4cOFCHFQMAAPiG3wTARYsWady4cRozZoyuvvpqLV26VCEhIVqxYkWlYxwOh6KiotxbZGRkHVYMAADgG34RAM+fP6+tW7cqLS3N3RYQEKC0tDTl5ORUOq6kpETt2rVTbGyshg0bpl27dtVFuQAAAD7lFwHw+PHjKisru2QFLzIyUgUFBRWO6dKli1asWKH33ntP//u//6vy8nL17dtXhw4dqnQ/paWlKi4u9tgAAAAaGr8IgJcjOTlZo0aNUmJiolJSUrR69WpFRERo2bJllY6ZP3++wsLC3FtsbGwdVgwAAOAdfhEAW7VqpcDAQBUWFnq0FxYWKioqqlqP0bhxY/Xs2VN79+6ttE9mZqZcLpd7O3jw4C+qGwAAwBf8IgAGBQWpV69eys7OdreVl5crOztbycnJ1XqMsrIy7dy5U9HR0ZX2cTqdCg0N9dgAAAAamka+LsBbJk2apIyMDPXu3Vt9+vTR4sWLdfr0aY0ZM0aSNGrUKF155ZWaP3++JGnOnDm6/vrr1bFjRxUVFenZZ5/VgQMHdO+99/ryMAAAAGqd3wTAO+64Q8eOHdOMGTNUUFCgxMRErVmzxn1iSH5+vgIC/v+C53fffadx48apoKBAzZs3V69evfTJJ5/o6quv9tUhAAAA1AmHMcb4uoiGqri4WGFhYXK5XHwcDABAA8H7t598BxAAAADVRwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyjXxdAADUF8YYnThxQiUlJWratKlatmwph8Ph67IAwOtYAQRgvaKiIr3wwgvq1KmTIiIi1KFDB0VERKhTp0564YUXVFRU5OsSAcCrHMYY4+siGqri4mKFhYXJ5XIpNDTU1+UAuAxr167VrbfeqjNnzki6sAp40cXVv5CQEL399ttKT0/3SY0AvIv3b1YAAVhs7dq1GjJkiM6ePStjjH769/DFtrNnz2rIkCFau3atjyoFAO/yqwC4ZMkStW/fXsHBwUpKStKnn35aZf+33npLCQkJCg4OVrdu3fThhx/WUaUAfK2oqEi33nqrjDEqLy+vsm95ebmMMbr11lv5OBh1pqSkRHPnzlVsbKwCAwMVGxuruXPnqqSkxNelwQ/4TQB84403NGnSJM2cOVPbtm1Tjx49lJ6erqNHj1bY/5NPPtHIkSM1duxYbd++XcOHD9fw4cP1xRdf1HHlUG6u9MorF/4F6siqVat05syZnw1/F5WXl+vMmTP661//WsuVoSG5/rmrFTjboTv+1N+rj1tSUqKUlBTNmjVLhw4dUnl5uQ4dOqRZs2YpJSWldkMgr8l2MH6iT58+5sEHH3TfLisrMzExMWb+/PkV9r/99tvNkCFDPNqSkpLMhAkTqr1Pl8tlJBmXy3V5RcOYqVONkf7/NnWqryuCBcrLy018fLxxOBxGUrU3h8Nh4uPjTXl5ua8PAfXA6qmDjJ6Se/Pm69ecOXNMQEBAhc/DgIAAM2fOHK/ty4Mlr8m8fxvjFyuA58+f19atW5WWluZuCwgIUFpamnJyciock5OT49FfktLT0yvtL0mlpaUqLi722PAL5OZKCxd6ti1cyF+dqHUnTpzQvn37LvnO388xxmjfvn06efJkLVWGBiM3V/ef/sij6clN3nv9Wr58eaWr0+Xl5Vq+fLlX9uOB12Sr+EUAPH78uMrKyhQZGenRHhkZqYKCggrHFBQU1Ki/JM2fP19hYWHuLTY29pcXb7OvvqpZO+Alv/Tjs1OnTnmpEjRYX30lV4hnU160vPb69Z///OcX3X9ZeE22il8EwLqSmZkpl8vl3g4ePOjrkhq2zp1r1g54SdOmTX/R+GbNmnmpEjRYnTsr5icfAv1up7z2+hUTE/OL7r8svCZbxS8CYKtWrRQYGKjCwkKP9sLCQkVFRVU4Jioqqkb9JcnpdCo0NNRjwy+QlCRNnerZNm3ahXagFrVs2VLx8fE1/pUPh8Oh+Ph4tWjRopYqQ4ORlKR/hI32aBp9m/dev8aPH6+AgIrfogMCAjR+/Hiv7McDr8lW8YsAGBQUpF69eik7O9vdVl5eruzsbCUnJ1c4Jjk52aO/JK1bt67S/qglCxZI//qX9Ne/Xvj3mWd8XREs4HA49PDDD1/W2EceeYSfh4Mk6aoFKxXXJFoKlPo3v9arr1+PPfaYEhMTLwmBAQEBSkxM1GOPPea1fXngNdkafvNLIG+88YYyMjK0bNky9enTR4sXL9abb76pL7/8UpGRkRo1apSuvPJKzZ8/X9KFy8CkpKTomWee0ZAhQ/T666/rv//7v7Vt2zZdc8011donVxIHGq6ioiK1adNGZ8+erdalYAICAtSkSRMdOnRI4eHhtV8grFdSUqLnn39ey5cv13/+8x/FxMRo/Pjxeuyxx37x1xhsx/u31MjXBXjLHXfcoWPHjmnGjBkqKChQYmKi1qxZ4z7RIz8/3+Mvqb59++q1117Tk08+qSeeeEKdOnXSu+++W+3wB6BhCw8P19tvv60hQ4YoICCgyhAYEBAgh8Oh1atXE/5QZ5o2baqnnnpKTz31lK9LgR/ymxVAX+AvCKDhq+5vAa9evVoDBgzwSY0AvIv3bz/5DiAAXK709HQdOnRIixcvVlxcnMd9cXFxWrx4sQ4fPkz4A+BXWAH8BfgLAvAvxhidPHlSp06dUrNmzdSiRQtO+AD8EO/ffvQdQAD4pRwOh1q2bKmWLVv6uhQAqFV8BAwAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFjGLwLgyZMndddddyk0NFTh4eEaO3asSkpKqhyTmpoqh8Phsd133311VDEAAIDvNPJ1Ad5w11136ciRI1q3bp2+//57jRkzRuPHj9drr71W5bhx48Zpzpw57tshISG1XSoAAIDPNfgAuHv3bq1Zs0afffaZevfuLUl68cUXNXjwYD333HOKiYmpdGxISIiioqLqqlQAAIB6ocF/BJyTk6Pw8HB3+JOktLQ0BQQEKDc3t8qxr776qlq1aqVrrrlGmZmZOnPmTJX9S0tLVVxc7LEBAAA0NA1+BbCgoECtW7f2aGvUqJFatGihgoKCSsfdeeedateunWJiYrRjxw5NmzZNe/bs0erVqysdM3/+fM2ePdtrtQMAAPhCvQ2A06dP14IFC6rss3v37st+/PHjx7v/u1u3boqOjlb//v21b98+xcfHVzgmMzNTkyZNct8uLi5WbGzsZdcAAADgC/U2AE6ePFmjR4+usk9cXJyioqJ09OhRj/YffvhBJ0+erNH3+5KSkiRJe/furTQAOp1OOZ3Oaj8mAABAfVRvA2BERIQiIiJ+tl9ycrKKioq0detW9erVS5K0YcMGlZeXu0NddeTl5UmSoqOjL6teAACAhqLBnwRy1VVXaeDAgRo3bpw+/fRTbdmyRQ899JBGjBjhPgP48OHDSkhI0KeffipJ2rdvn+bOnautW7fq22+/1fvvv69Ro0apX79+6t69uy8PBwAAoNY1+AAoXTibNyEhQf3799fgwYP161//WsuXL3ff//3332vPnj3us3yDgoK0fv16DRgwQAkJCZo8ebJuvfVWffDBB746BAAAgDrjMMYYXxfRUBUXFyssLEwul0uhoaG+LgcAAFQD799+sgIIAACA6iMAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgmUa+LgAA4Kdyc6WvvpI6d5aSknxdDYAfYQUQAOB906ZJ118vjRp14d9p03xdEYAfIQACALwrN1dauNCzbeHCC+0A6gUCIADAu776qmbtAOocARAA4F2dO9esHUCdIwACALwrKUmaOtWzbdo0TgQB6hHOAgYAeN+CBdItt3AWMFBPEQABALUjKYngB9RTfAQMAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJbxiwD49NNPq2/fvgoJCVF4eHi1xhhjNGPGDEVHR6tJkyZKS0vT119/XbuFAgAA1AN+EQDPnz+v2267Tffff3+1xyxcuFB//OMftXTpUuXm5uqKK65Qenq6zp07V4uVAgAA+J7DGGN8XYS3ZGVlaeLEiSoqKqqynzFGMTExmjx5sh5//HFJksvlUmRkpLKysjRixIhq7a+4uFhhYWFyuVwKDQ39peUDAIA6wPu3n6wA1tT+/ftVUFCgtLQ0d1tYWJiSkpKUk5Pjw8oAAABqn5W/BVxQUCBJioyM9GiPjIx031eR0tJSlZaWum8XFxfXToEAAAC1qN6uAE6fPl0Oh6PK7csvv6zTmubPn6+wsDD3FhsbW6f7BwAA8IZ6uwI4efJkjR49uso+cXFxl/XYUVFRkqTCwkJFR0e72wsLC5WYmFjpuMzMTE2aNMl9u7i4mBAIAAAanHobACMiIhQREVErj92hQwdFRUUpOzvbHfiKi4uVm5tb5ZnETqdTTqezVmoCAACoK/X2I+CayM/PV15envLz81VWVqa8vDzl5eWppKTE3SchIUHvvPOOJMnhcGjixImaN2+e3n//fe3cuVOjRo1STEyMhg8f7qOjAAAAqBv1dgWwJmbMmKFVq1a5b/fs2VOStHHjRqWmpkqS9uzZI5fL5e4zdepUnT59WuPHj1dRUZF+/etfa82aNQoODq7T2gEAAOqaX10HsK5xHSEAABoe3r/95CNgAAAAVB8BEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMv4RQB8+umn1bdvX4WEhCg8PLxaY0aPHi2Hw+GxDRw4sHYLBQAAqAca+boAbzh//rxuu+02JScn63/+53+qPW7gwIFauXKl+7bT6ayN8gAAAOoVvwiAs2fPliRlZWXVaJzT6VRUVFQtVAQAAFB/+UUAvFybNm1S69at1bx5c914442aN2+eWrZsWWn/0tJSlZaWum+7XC5JUnFxca3XCgAAvOPi+7YxxseV+I61AXDgwIG65ZZb1KFDB+3bt09PPPGEBg0apJycHAUGBlY4Zv78+e7Vxh+LjY2t7XIBAICXnThxQmFhYb4uwyccpp7G3+nTp2vBggVV9tm9e7cSEhLct7OysjRx4kQVFRXVeH/ffPON4uPjtX79evXv37/CPj9dASwqKlK7du2Un59v7RPIW4qLixUbG6uDBw8qNDTU1+U0WMyj9zCX3sNcegfz6D0ul0tt27bVd999V+2TR/1NvV0BnDx5skaPHl1ln7i4OK/tLy4uTq1atdLevXsrDYBOp7PCE0XCwsL4P6OXhIaGMpdewDx6D3PpPcyldzCP3hMQ4BcXQ7ks9TYARkREKCIios72d+jQIZ04cULR0dF1tk8AAABf8Ivom5+fr7y8POXn56usrEx5eXnKy8tTSUmJu09CQoLeeecdSVJJSYmmTJmif/3rX/r222+VnZ2tYcOGqWPHjkpPT/fVYQAAANSJersCWBMzZszQqlWr3Ld79uwpSdq4caNSU1MlSXv27HGftRsYGKgdO3Zo1apVKioqUkxMjAYMGKC5c+fW6FqATqdTM2fO5PqBXsBcegfz6D3Mpfcwl97BPHoPc1mPTwIBAABA7fCLj4ABAABQfQRAAAAAyxAAAQAALEMABAAAsAwBsIaefvpp9e3bVyEhIdW+evjo0aPlcDg8toEDB9ZuofXc5cyjMUYzZsxQdHS0mjRporS0NH399de1W2gDcPLkSd11110KDQ1VeHi4xo4d63EJpIqkpqZe8py877776qji+mPJkiVq3769goODlZSUpE8//bTK/m+99ZYSEhIUHBysbt266cMPP6yjSuu/msxlVlbWJc+/4ODgOqy2ftq8ebOGDh2qmJgYORwOvfvuuz87ZtOmTbr22mvldDrVsWNHZWVl1XqdDUFN53LTpk2XPCcdDocKCgrqpmAfIADW0Pnz53Xbbbfp/vvvr9G4gQMH6siRI+7tb3/7Wy1V2DBczjwuXLhQf/zjH7V06VLl5ubqiiuuUHp6us6dO1eLldZ/d911l3bt2qV169bpH//4hzZv3qzx48f/7Lhx48Z5PCcXLlxYB9XWH2+88YYmTZqkmTNnatu2berRo4fS09N19OjRCvt/8sknGjlypMaOHavt27dr+PDhGj58uL744os6rrz+qelcShd+zeLHz78DBw7UYcX10+nTp9WjRw8tWbKkWv3379+vIUOG6IYbblBeXp4mTpyoe++9V2vXrq3lSuu/ms7lRXv27PF4XrZu3bqWKqwHDC7LypUrTVhYWLX6ZmRkmGHDhtVqPQ1VdeexvLzcREVFmWeffdbdVlRUZJxOp/nb3/5WixXWb//+97+NJPPZZ5+52z766CPjcDjM4cOHKx2XkpJiHn300TqosP7q06ePefDBB923y8rKTExMjJk/f36F/W+//XYzZMgQj7akpCQzYcKEWq2zIajpXNbk9dNWksw777xTZZ+pU6earl27erTdcccdJj09vRYra3iqM5cbN240ksx3331XJzXVB6wA1pFNmzapdevW6tKli+6//36dOHHC1yU1KPv371dBQYHS0tLcbWFhYUpKSlJOTo4PK/OtnJwchYeHq3fv3u62tLQ0BQQEKDc3t8qxr776qlq1aqVrrrlGmZmZOnPmTG2XW2+cP39eW7du9Xg+BQQEKC0trdLnU05Ojkd/SUpPT7f6+Sdd3lxKF36RqV27doqNjdWwYcO0a9euuijXr/Cc9L7ExERFR0frpptu0pYtW3xdTq3yi18Cqe8GDhyoW265RR06dNC+ffv0xBNPaNCgQcrJyVFgYKCvy2sQLn4PIzIy0qM9MjLSr7+j8XMKCgou+YiiUaNGatGiRZXzcuedd6pdu3aKiYnRjh07NG3aNO3Zs0erV6+u7ZLrhePHj6usrKzC59OXX35Z4ZiCggKefxW4nLns0qWLVqxYoe7du8vlcum5555T3759tWvXLrVp06YuyvYLlT0ni4uLdfbsWTVp0sRHlTU80dHRWrp0qXr37q3S0lK9/PLLSk1NVW5urq699lpfl1crCICSpk+frgULFlTZZ/fu3UpISLisxx8xYoT7v7t166bu3bsrPj5emzZtUv/+/S/rMeuj2p5Hm1R3Li/Xj78j2K1bN0VHR6t///7at2+f4uPjL/txgepITk5WcnKy+3bfvn111VVXadmyZZo7d64PK4OtunTpoi5durhv9+3bV/v27dPzzz+vV155xYeV1R4CoKTJkydr9OjRVfaJi4vz2v7i4uLUqlUr7d27168CYG3OY1RUlCSpsLBQ0dHR7vbCwkIlJiZe1mPWZ9Wdy6ioqEu+aP/DDz/o5MmT7jmrjqSkJEnS3r17rQiArVq1UmBgoAoLCz3aCwsLK523qKioGvW3xeXM5U81btxYPXv21N69e2ujRL9V2XMyNDSU1T8v6NOnjz7++GNfl1FrCICSIiIiFBERUWf7O3TokE6cOOERZPxBbc5jhw4dFBUVpezsbHfgKy4uVm5ubo3PyG4IqjuXycnJKioq0tatW9WrVy9J0oYNG1ReXu4OddWRl5cnSX73nKxMUFCQevXqpezsbA0fPlySVF5eruzsbD300EMVjklOTlZ2drYmTpzoblu3bp3HSpaNLmcuf6qsrEw7d+7U4MGDa7FS/5OcnHzJpYh4TnpPXl6ef78m+voslIbmwIEDZvv27Wb27NmmadOmZvv27Wb79u3m1KlT7j5dunQxq1evNsYYc+rUKfP444+bnJwcs3//frN+/Xpz7bXXmk6dOplz58756jB8rqbzaIwxzzzzjAkPDzfvvfee2bFjhxk2bJjp0KGDOXv2rC8Ood4YOHCg6dmzp8nNzTUff/yx6dSpkxk5cqT7/kOHDpkuXbqY3NxcY4wxe/fuNXPmzDGff/652b9/v3nvvfdMXFyc6devn68OwSdef/1143Q6TVZWlvn3v/9txo8fb8LDw01BQYExxph77rnHTJ8+3d1/y5YtplGjRua5554zu3fvNjNnzjSNGzc2O3fu9NUh1Bs1ncvZs2ebtWvXmn379pmtW7eaESNGmODgYLNr1y5fHUK9cOrUKfdroSSzaNEis337dnPgwAFjjDHTp08399xzj7v/N998Y0JCQsyUKVPM7t27zZIlS0xgYKBZs2aNrw6h3qjpXD7//PPm3XffNV9//bXZuXOnefTRR01AQIBZv369rw6h1hEAaygjI8NIumTbuHGju48ks3LlSmOMMWfOnDEDBgwwERERpnHjxqZdu3Zm3Lhx7hdGW9V0Ho25cCmYp556ykRGRhqn02n69+9v9uzZU/fF1zMnTpwwI0eONE2bNjWhoaFmzJgxHkF6//79HnObn59v+vXrZ1q0aGGcTqfp2LGjmTJlinG5XD46At958cUXTdu2bU1QUJDp06eP+de//uW+LyUlxWRkZHj0f/PNN03nzp1NUFCQ6dq1q/nnP/9ZxxXXXzWZy4kTJ7r7RkZGmsGDB5tt27b5oOr65eKlSH66XZy7jIwMk5KScsmYxMREExQUZOLi4jxeM21W07lcsGCBiY+PN8HBwaZFixYmNTXVbNiwwTfF1xGHMcbU2XIjAAAAfI7rAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZf4f/ltOg/Jh/ZEAAAAASUVORK5CYII=",
      "text/html": [
       "\n",
       "            <div style=\"display: inline-block;\">\n",
       "                <div class=\"jupyter-widgets widget-label\" style=\"text-align: center;\">\n",
       "                    Figure\n",
       "                </div>\n",
       "                <img src='data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxf0lEQVR4nO3deXTU1f3/8dckkAkRkrCELBKWhCWKQBAkhrYkSiQsUjhaFVwIFAF3EWSJX2X3i6BFrMUC9QuhfrUuFbdW4QsBDkdMowI5IEUURAKUhM1MCEvQ5P7+4DA/R5KY4CSTzH0+zvkcnDv3zuf9uZ3OvHJnPp9xGGOMAAAAYI0AXxcAAACAukUABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAANZJTU1Vamqqr8u4RPv27TV69Gj37U2bNsnhcGjTpk11Wsfo0aPVvn37Ot0ngLpFAAQsNGvWLDkcDh0/frzKfk8//bR++9vfKjIyUg6HQ7NmzaqbAivRvn17ORyOCrdz5875tLb6IDU1tdL5+fLLL31dHoB6pJGvCwBQfz355JOKiopSz549tXbtWl+XI0lKTEzU5MmTL2kPCgryQTW1q1+/fjp79myNjq1NmzaaP3/+Je0xMTHeLA1AA0cABFCp/fv3q3379jp+/LgiIiJ8UoMxRufOnVOTJk0kSVdeeaXuvvtun9RS1wICAhQcHPyz/c6cOaOQkBBJUlhYmDXzA+Dy8REw4EfOnj2rhIQEJSQk6OzZs+72kydPKjo6Wn379lVZWVm1H89b3wO7+JHzT2VlZcnhcOjbb7/12OfNN9+stWvXqnfv3mrSpImWLVtWrf2sXLlSN954o1q3bi2n06mrr75af/7zn6s19sUXX1TXrl0VEhKi5s2bq3fv3nrttdc8+hw+fFi///3vFRkZKafTqa5du2rFihXVevwfM8Zo3rx5atOmjUJCQnTDDTdo165dl/Sr6DuAqampuuaaa7R161b169dPISEheuKJJ6q13/fee09DhgxRTEyMnE6n4uPjNXfu3Go9J15//XX16tVLzZo1U2hoqLp166YXXnjBo09RUZEmTpyo2NhYOZ1OdezYUQsWLFB5eXm16gNQd1gBBPxIkyZNtGrVKv3qV7/Sf/3Xf2nRokWSpAcffFAul0tZWVkKDAz0cZU/b8+ePRo5cqQmTJigcePGqUuXLu77vv/++0u+uxgSEqKQkBD9+c9/VteuXfXb3/5WjRo10gcffKAHHnhA5eXlevDBByvd31/+8hc98sgj+t3vfqdHH31U586d044dO5Sbm6s777xTklRYWKjrr79eDodDDz30kCIiIvTRRx9p7NixKi4u1sSJE6t9fDNmzNC8efM0ePBgDR48WNu2bdOAAQN0/vz5ao0/ceKEBg0apBEjRujuu+9WZGSk+76ysrJL5ic4OFhNmzZVVlaWmjZtqkmTJqlp06basGGDZsyYoeLiYj377LOV7m/dunUaOXKk+vfvrwULFkiSdu/erS1btujRRx+VdGEVMiUlRYcPH9aECRPUtm1bffLJJ8rMzNSRI0e0ePHias8PgDpgAPidzMxMExAQYDZv3mzeeustI8ksXrzYff/MmTONJHPs2LFqPd6xY8eMJDNz5szLqufi/n5q5cqVRpLZv3+/u61du3ZGklmzZs0l/S/e99PtYl1nzpy5ZEx6erqJi4vzaEtJSTEpKSnu28OGDTNdu3at8hjGjh1roqOjzfHjxz3aR4wYYcLCwircd0WOHj1qgoKCzJAhQ0x5ebm7/YknnjCSTEZGhrtt48aNRpLZuHGjR+2SzNKlSy957Iv3/XS7+JgV1ThhwgQTEhJizp07527LyMgw7dq1c99+9NFHTWhoqPnhhx8qPa65c+eaK664wnz11Vce7dOnTzeBgYEmPz+/0rEA6h4fAQN+aNasWeratasyMjL0wAMPKCUlRY888oivy6q2Dh06KD09vcL7kpKStG7dOo9t1KhRkuT+nqAkuVwuHT9+XCkpKfrmm2/kcrkq3V94eLgOHTqkzz77rML7jTF6++23NXToUBljdPz4cfeWnp4ul8ulbdu2VevY1q9fr/Pnz+vhhx/2+Fi8JiuITqdTY8aMqfC+9u3bXzI/U6dOleQ5P6dOndLx48f1m9/8RmfOnKnyLOHw8HCdPn1a69atq7TPW2+9pd/85jdq3ry5x/ykpaWprKxMmzdvrvbxAah9fAQM+KGgoCCtWLFC1113nYKDg7Vy5coKv4NXX3Xo0KHS+1q1aqW0tLQK79uyZYtmzpypnJwcnTlzxuM+l8ulsLCwCsdNmzZN69evV58+fdSxY0cNGDBAd955p371q19Jko4dO6aioiItX75cy5cvr/Axjh49Wp1D04EDByRJnTp18miPiIhQ8+bNq/UYV155ZaVnBl9xxRWVzs+uXbv05JNPasOGDSouLva4r6qA/MADD+jNN9/UoEGDdOWVV2rAgAG6/fbbNXDgQHefr7/+Wjt27Kj0ZKHqzg+AukEABPzUxcu2nDt3Tl9//XWVoaq2VRY+Kzv54McrVdW1b98+9e/fXwkJCVq0aJFiY2MVFBSkDz/8UM8//3yVJyJcddVV2rNnj/7xj39ozZo1evvtt/XSSy9pxowZmj17tnvs3XffrYyMjAofo3v37jWu+XJdzvwUFRUpJSVFoaGhmjNnjuLj4xUcHKxt27Zp2rRpVc5P69atlZeXp7Vr1+qjjz7SRx99pJUrV2rUqFFatWqVJKm8vFw33XSTe7Xxpzp37lzjmgHUHgIg4Id27NihOXPmaMyYMcrLy9O9996rnTt3VroCVtsurmwVFRUpPDzc3X5xNcwbPvjgA5WWlur9999X27Zt3e0bN26s1vgrrrhCd9xxh+644w6dP39et9xyi55++mllZmYqIiJCzZo1U1lZWaWra9XVrl07SRdWzOLi4tztx44d03ffffeLHrsqmzZt0okTJ7R69Wr169fP3b5///5qjQ8KCtLQoUM1dOhQlZeX64EHHtCyZcv01FNPqWPHjoqPj1dJSckvnh8AdYPvAAJ+5vvvv9fo0aMVExOjF154QVlZWSosLNRjjz3ms5ri4+MlyeN7YKdPn3avHnnDxbObjTHuNpfLpZUrV/7s2BMnTnjcDgoK0tVXXy1jjL7//nsFBgbq1ltv1dtvv60vvvjikvHHjh2rdp1paWlq3LixXnzxRY9aa/ss2Yrm5/z583rppZd+duxP5ycgIMC94llaWipJuv3225WTk1PhBcOLior0ww8/XHbtALyPFUDAz8ybN095eXnKzs5Ws2bN1L17d82YMUNPPvmkfve732nw4MHuvosWLXJfQPiigIAA93XlXnnlFR04cMD9fbrNmzdr3rx5kqR77rnHvZr1cwYMGKC2bdtq7NixmjJligIDA7VixQpFREQoPz/fG4etAQMGuFepJkyYoJKSEv3lL39R69atdeTIkZ8dGxUVpV/96leKjIzU7t279ac//UlDhgxRs2bNJEnPPPOMNm7cqKSkJI0bN05XX321Tp48qW3btmn9+vU6efJkteqMiIjQ448/rvnz5+vmm2/W4MGDtX37dn300Udq1arVL56HyvTt21fNmzdXRkaGHnnkETkcDr3yyisegbAy9957r06ePKkbb7xRbdq00YEDB/Tiiy8qMTFRV111lSRpypQpev/993XzzTdr9OjR6tWrl06fPq2dO3fq73//u7799ttaPT4ANeTDM5ABeNnWrVtNo0aNzMMPP+zR/sMPP5jrrrvOxMTEmO+++859WZaKtsDAQPe4yi4rop9cmqS6tSUlJZmgoCDTtm1bs2jRokovAzNkyJAKH6Oq+4wx5v333zfdu3c3wcHBpn379mbBggVmxYoVl+zjp5eBWbZsmenXr59p2bKlcTqdJj4+3kyZMsW4XC6Pxy8sLDQPPvigiY2NNY0bNzZRUVGmf//+Zvny5TWai7KyMjN79mwTHR1tmjRpYlJTU80XX3xh2rVrV63LwFR2yZqq7jPGmC1btpjrr7/eNGnSxMTExJipU6eatWvXXrKPn14G5u9//7sZMGCAad26tft/vwkTJpgjR454PP6pU6dMZmam6dixowkKCjKtWrUyffv2Nc8995w5f/58jeYIQO1yGFONP/8AAADgN/gOIAAAgGX4DiCAy3b27Nkqrx8nSS1atKj0mnX+5tixY1X+rm5QUJBatGhRhxUBQMX8ZgVw8+bNGjp0qGJiYuRwOPTuu+9W2f/ij6z/dCsoKKibggE/8MYbbyg6OrrK7ZNPPvF1mXXmuuuuq3IubrnlFl+XCACS/GgF8PTp0+rRo4d+//vf1+hFds+ePQoNDXXfbt26dW2UB/il9PT0Kn8eTJJ69OhRR9X43quvvqqzZ89Wen91f+kDAGqb3wTAQYMGadCgQTUe17p1a48L0wKovosrW7jg4k/HAUB95zcfAV+uxMRERUdH66abbtKWLVt8XQ4AAECt85sVwJqKjo7W0qVL1bt3b5WWlurll19WamqqcnNzde2111Y4prS01H3Ve+nCb1+ePHlSLVu2rPS3TgEAQP1ijNGpU6cUExOjgAA718L88jqADodD77zzjoYPH16jcSkpKWrbtq1eeeWVCu+fNWuWZs+e7YUKAQCArx08eFBt2rTxdRk+Ye0KYEX69Omjjz/+uNL7MzMzNWnSJPdtl8ultm3b6uDBgx4nkgAAgPqruLhYsbGx7p96tBEB8Efy8vKq/EK70+mU0+m8pD00NJQACABAA2Pz17f8JgCWlJRo79697tv79+9XXl6eWrRoobZt2yozM1OHDx/WX//6V0nS4sWL1aFDB3Xt2lXnzp3Tyy+/rA0bNuj//u//fHUIAAAAdcJvAuDnn3+uG264wX374ke1GRkZysrK0pEjR5Sfn+++//z585o8ebIOHz6skJAQde/eXevXr/d4DAAAAH/klyeB1JXi4mKFhYXJ5XLxETAAAA0E799cBxAAAMA6BEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAs4zcBcPPmzRo6dKhiYmLkcDj07rvv/uyYTZs26dprr5XT6VTHjh2VlZVV63UCAAD4mt8EwNOnT6tHjx5asmRJtfrv379fQ4YM0Q033KC8vDxNnDhR9957r9auXVvLlQIAAPhWI18X4C2DBg3SoEGDqt1/6dKl6tChg/7whz9Ikq666ip9/PHHev7555Wenl5bZQIAAPic36wA1lROTo7S0tI82tLT05WTk1PpmNLSUhUXF3tsAAAADY21AbCgoECRkZEebZGRkSouLtbZs2crHDN//nyFhYW5t9jY2LooFQAAwKusDYCXIzMzUy6Xy70dPHjQ1yUBAADUmN98B7CmoqKiVFhY6NFWWFio0NBQNWnSpMIxTqdTTqezLsoDAACoNdauACYnJys7O9ujbd26dUpOTvZRRQAAAHXDbwJgSUmJ8vLylJeXJ+nCZV7y8vKUn58v6cLHt6NGjXL3v++++/TNN99o6tSp+vLLL/XSSy/pzTff1GOPPeaL8gEAAOqM3wTAzz//XD179lTPnj0lSZMmTVLPnj01Y8YMSdKRI0fcYVCSOnTooH/+859at26devTooT/84Q96+eWXuQQMAADwew5jjPF1EQ1VcXGxwsLC5HK5FBoa6utyAABANfD+7UcrgAAAAKgeAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAlrH2t4ABALUsN1f66iupc2cpKcnX1QD4EVYAAQDeN22adP310qhRF/6dNs3XFQH4EQIgAMC7cnOlhQs92xYuvNAOoF4gAAIAvOurr2rWDqDOEQABAN7VuXPN2gHUOQIgAMC7kpKkqVM926ZN40QQoB7hLGAAgPctWCDdcgtnAQP1FAEQAFA7kpIIfkA9xUfAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGb8KgEuWLFH79u0VHByspKQkffrpp5X2zcrKksPh8NiCg4PrsFoAAADf8JsA+MYbb2jSpEmaOXOmtm3bph49eig9PV1Hjx6tdExoaKiOHDni3g4cOFCHFQMAAPiG3wTARYsWady4cRozZoyuvvpqLV26VCEhIVqxYkWlYxwOh6KiotxbZGRkHVYMAADgG34RAM+fP6+tW7cqLS3N3RYQEKC0tDTl5ORUOq6kpETt2rVTbGyshg0bpl27dtVFuQAAAD7lFwHw+PHjKisru2QFLzIyUgUFBRWO6dKli1asWKH33ntP//u//6vy8nL17dtXhw4dqnQ/paWlKi4u9tgAAAAaGr8IgJcjOTlZo0aNUmJiolJSUrR69WpFRERo2bJllY6ZP3++wsLC3FtsbGwdVgwAAOAdfhEAW7VqpcDAQBUWFnq0FxYWKioqqlqP0bhxY/Xs2VN79+6ttE9mZqZcLpd7O3jw4C+qGwAAwBf8IgAGBQWpV69eys7OdreVl5crOztbycnJ1XqMsrIy7dy5U9HR0ZX2cTqdCg0N9dgAAAAamka+LsBbJk2apIyMDPXu3Vt9+vTR4sWLdfr0aY0ZM0aSNGrUKF155ZWaP3++JGnOnDm6/vrr1bFjRxUVFenZZ5/VgQMHdO+99/ryMAAAAGqd3wTAO+64Q8eOHdOMGTNUUFCgxMRErVmzxn1iSH5+vgIC/v+C53fffadx48apoKBAzZs3V69evfTJJ5/o6quv9tUhAAAA1AmHMcb4uoiGqri4WGFhYXK5XHwcDABAA8H7t598BxAAAADVRwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyjXxdAADUF8YYnThxQiUlJWratKlatmwph8Ph67IAwOtYAQRgvaKiIr3wwgvq1KmTIiIi1KFDB0VERKhTp0564YUXVFRU5OsSAcCrHMYY4+siGqri4mKFhYXJ5XIpNDTU1+UAuAxr167VrbfeqjNnzki6sAp40cXVv5CQEL399ttKT0/3SY0AvIv3b1YAAVhs7dq1GjJkiM6ePStjjH769/DFtrNnz2rIkCFau3atjyoFAO/yqwC4ZMkStW/fXsHBwUpKStKnn35aZf+33npLCQkJCg4OVrdu3fThhx/WUaUAfK2oqEi33nqrjDEqLy+vsm95ebmMMbr11lv5OBh1pqSkRHPnzlVsbKwCAwMVGxuruXPnqqSkxNelwQ/4TQB84403NGnSJM2cOVPbtm1Tjx49lJ6erqNHj1bY/5NPPtHIkSM1duxYbd++XcOHD9fw4cP1xRdf1HHlUG6u9MorF/4F6siqVat05syZnw1/F5WXl+vMmTP661//WsuVoSG5/rmrFTjboTv+1N+rj1tSUqKUlBTNmjVLhw4dUnl5uQ4dOqRZs2YpJSWldkMgr8l2MH6iT58+5sEHH3TfLisrMzExMWb+/PkV9r/99tvNkCFDPNqSkpLMhAkTqr1Pl8tlJBmXy3V5RcOYqVONkf7/NnWqryuCBcrLy018fLxxOBxGUrU3h8Nh4uPjTXl5ua8PAfXA6qmDjJ6Se/Pm69ecOXNMQEBAhc/DgIAAM2fOHK/ty4Mlr8m8fxvjFyuA58+f19atW5WWluZuCwgIUFpamnJyciock5OT49FfktLT0yvtL0mlpaUqLi722PAL5OZKCxd6ti1cyF+dqHUnTpzQvn37LvnO388xxmjfvn06efJkLVWGBiM3V/ef/sij6clN3nv9Wr58eaWr0+Xl5Vq+fLlX9uOB12Sr+EUAPH78uMrKyhQZGenRHhkZqYKCggrHFBQU1Ki/JM2fP19hYWHuLTY29pcXb7OvvqpZO+Alv/Tjs1OnTnmpEjRYX30lV4hnU160vPb69Z///OcX3X9ZeE22il8EwLqSmZkpl8vl3g4ePOjrkhq2zp1r1g54SdOmTX/R+GbNmnmpEjRYnTsr5icfAv1up7z2+hUTE/OL7r8svCZbxS8CYKtWrRQYGKjCwkKP9sLCQkVFRVU4Jioqqkb9JcnpdCo0NNRjwy+QlCRNnerZNm3ahXagFrVs2VLx8fE1/pUPh8Oh+Ph4tWjRopYqQ4ORlKR/hI32aBp9m/dev8aPH6+AgIrfogMCAjR+/Hiv7McDr8lW8YsAGBQUpF69eik7O9vdVl5eruzsbCUnJ1c4Jjk52aO/JK1bt67S/qglCxZI//qX9Ne/Xvj3mWd8XREs4HA49PDDD1/W2EceeYSfh4Mk6aoFKxXXJFoKlPo3v9arr1+PPfaYEhMTLwmBAQEBSkxM1GOPPea1fXngNdkafvNLIG+88YYyMjK0bNky9enTR4sXL9abb76pL7/8UpGRkRo1apSuvPJKzZ8/X9KFy8CkpKTomWee0ZAhQ/T666/rv//7v7Vt2zZdc8011donVxIHGq6ioiK1adNGZ8+erdalYAICAtSkSRMdOnRI4eHhtV8grFdSUqLnn39ey5cv13/+8x/FxMRo/Pjxeuyxx37x1xhsx/u31MjXBXjLHXfcoWPHjmnGjBkqKChQYmKi1qxZ4z7RIz8/3+Mvqb59++q1117Tk08+qSeeeEKdOnXSu+++W+3wB6BhCw8P19tvv60hQ4YoICCgyhAYEBAgh8Oh1atXE/5QZ5o2baqnnnpKTz31lK9LgR/ymxVAX+AvCKDhq+5vAa9evVoDBgzwSY0AvIv3bz/5DiAAXK709HQdOnRIixcvVlxcnMd9cXFxWrx4sQ4fPkz4A+BXWAH8BfgLAvAvxhidPHlSp06dUrNmzdSiRQtO+AD8EO/ffvQdQAD4pRwOh1q2bKmWLVv6uhQAqFV8BAwAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFjGLwLgyZMndddddyk0NFTh4eEaO3asSkpKqhyTmpoqh8Phsd133311VDEAAIDvNPJ1Ad5w11136ciRI1q3bp2+//57jRkzRuPHj9drr71W5bhx48Zpzpw57tshISG1XSoAAIDPNfgAuHv3bq1Zs0afffaZevfuLUl68cUXNXjwYD333HOKiYmpdGxISIiioqLqqlQAAIB6ocF/BJyTk6Pw8HB3+JOktLQ0BQQEKDc3t8qxr776qlq1aqVrrrlGmZmZOnPmTJX9S0tLVVxc7LEBAAA0NA1+BbCgoECtW7f2aGvUqJFatGihgoKCSsfdeeedateunWJiYrRjxw5NmzZNe/bs0erVqysdM3/+fM2ePdtrtQMAAPhCvQ2A06dP14IFC6rss3v37st+/PHjx7v/u1u3boqOjlb//v21b98+xcfHVzgmMzNTkyZNct8uLi5WbGzsZdcAAADgC/U2AE6ePFmjR4+usk9cXJyioqJ09OhRj/YffvhBJ0+erNH3+5KSkiRJe/furTQAOp1OOZ3Oaj8mAABAfVRvA2BERIQiIiJ+tl9ycrKKioq0detW9erVS5K0YcMGlZeXu0NddeTl5UmSoqOjL6teAACAhqLBnwRy1VVXaeDAgRo3bpw+/fRTbdmyRQ899JBGjBjhPgP48OHDSkhI0KeffipJ2rdvn+bOnautW7fq22+/1fvvv69Ro0apX79+6t69uy8PBwAAoNY1+AAoXTibNyEhQf3799fgwYP161//WsuXL3ff//3332vPnj3us3yDgoK0fv16DRgwQAkJCZo8ebJuvfVWffDBB746BAAAgDrjMMYYXxfRUBUXFyssLEwul0uhoaG+LgcAAFQD799+sgIIAACA6iMAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgmUa+LgAA4Kdyc6WvvpI6d5aSknxdDYAfYQUQAOB906ZJ118vjRp14d9p03xdEYAfIQACALwrN1dauNCzbeHCC+0A6gUCIADAu776qmbtAOocARAA4F2dO9esHUCdIwACALwrKUmaOtWzbdo0TgQB6hHOAgYAeN+CBdItt3AWMFBPEQABALUjKYngB9RTfAQMAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJbxiwD49NNPq2/fvgoJCVF4eHi1xhhjNGPGDEVHR6tJkyZKS0vT119/XbuFAgAA1AN+EQDPnz+v2267Tffff3+1xyxcuFB//OMftXTpUuXm5uqKK65Qenq6zp07V4uVAgAA+J7DGGN8XYS3ZGVlaeLEiSoqKqqynzFGMTExmjx5sh5//HFJksvlUmRkpLKysjRixIhq7a+4uFhhYWFyuVwKDQ39peUDAIA6wPu3n6wA1tT+/ftVUFCgtLQ0d1tYWJiSkpKUk5Pjw8oAAABqn5W/BVxQUCBJioyM9GiPjIx031eR0tJSlZaWum8XFxfXToEAAAC1qN6uAE6fPl0Oh6PK7csvv6zTmubPn6+wsDD3FhsbW6f7BwAA8IZ6uwI4efJkjR49uso+cXFxl/XYUVFRkqTCwkJFR0e72wsLC5WYmFjpuMzMTE2aNMl9u7i4mBAIAAAanHobACMiIhQREVErj92hQwdFRUUpOzvbHfiKi4uVm5tb5ZnETqdTTqezVmoCAACoK/X2I+CayM/PV15envLz81VWVqa8vDzl5eWppKTE3SchIUHvvPOOJMnhcGjixImaN2+e3n//fe3cuVOjRo1STEyMhg8f7qOjAAAAqBv1dgWwJmbMmKFVq1a5b/fs2VOStHHjRqWmpkqS9uzZI5fL5e4zdepUnT59WuPHj1dRUZF+/etfa82aNQoODq7T2gEAAOqaX10HsK5xHSEAABoe3r/95CNgAAAAVB8BEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMsQAAEAACxDAAQAALAMARAAAMAyBEAAAADLEAABAAAsQwAEAACwDAEQAADAMgRAAAAAyxAAAQAALEMABAAAsAwBEAAAwDIEQAAAAMv4RQB8+umn1bdvX4WEhCg8PLxaY0aPHi2Hw+GxDRw4sHYLBQAAqAca+boAbzh//rxuu+02JScn63/+53+qPW7gwIFauXKl+7bT6ayN8gAAAOoVvwiAs2fPliRlZWXVaJzT6VRUVFQtVAQAAFB/+UUAvFybNm1S69at1bx5c914442aN2+eWrZsWWn/0tJSlZaWum+7XC5JUnFxca3XCgAAvOPi+7YxxseV+I61AXDgwIG65ZZb1KFDB+3bt09PPPGEBg0apJycHAUGBlY4Zv78+e7Vxh+LjY2t7XIBAICXnThxQmFhYb4uwyccpp7G3+nTp2vBggVV9tm9e7cSEhLct7OysjRx4kQVFRXVeH/ffPON4uPjtX79evXv37/CPj9dASwqKlK7du2Un59v7RPIW4qLixUbG6uDBw8qNDTU1+U0WMyj9zCX3sNcegfz6D0ul0tt27bVd999V+2TR/1NvV0BnDx5skaPHl1ln7i4OK/tLy4uTq1atdLevXsrDYBOp7PCE0XCwsL4P6OXhIaGMpdewDx6D3PpPcyldzCP3hMQ4BcXQ7ks9TYARkREKCIios72d+jQIZ04cULR0dF1tk8AAABf8Ivom5+fr7y8POXn56usrEx5eXnKy8tTSUmJu09CQoLeeecdSVJJSYmmTJmif/3rX/r222+VnZ2tYcOGqWPHjkpPT/fVYQAAANSJersCWBMzZszQqlWr3Ld79uwpSdq4caNSU1MlSXv27HGftRsYGKgdO3Zo1apVKioqUkxMjAYMGKC5c+fW6FqATqdTM2fO5PqBXsBcegfz6D3Mpfcwl97BPHoPc1mPTwIBAABA7fCLj4ABAABQfQRAAAAAyxAAAQAALEMABAAAsAwBsIaefvpp9e3bVyEhIdW+evjo0aPlcDg8toEDB9ZuofXc5cyjMUYzZsxQdHS0mjRporS0NH399de1W2gDcPLkSd11110KDQ1VeHi4xo4d63EJpIqkpqZe8py877776qji+mPJkiVq3769goODlZSUpE8//bTK/m+99ZYSEhIUHBysbt266cMPP6yjSuu/msxlVlbWJc+/4ODgOqy2ftq8ebOGDh2qmJgYORwOvfvuuz87ZtOmTbr22mvldDrVsWNHZWVl1XqdDUFN53LTpk2XPCcdDocKCgrqpmAfIADW0Pnz53Xbbbfp/vvvr9G4gQMH6siRI+7tb3/7Wy1V2DBczjwuXLhQf/zjH7V06VLl5ubqiiuuUHp6us6dO1eLldZ/d911l3bt2qV169bpH//4hzZv3qzx48f/7Lhx48Z5PCcXLlxYB9XWH2+88YYmTZqkmTNnatu2berRo4fS09N19OjRCvt/8sknGjlypMaOHavt27dr+PDhGj58uL744os6rrz+qelcShd+zeLHz78DBw7UYcX10+nTp9WjRw8tWbKkWv3379+vIUOG6IYbblBeXp4mTpyoe++9V2vXrq3lSuu/ms7lRXv27PF4XrZu3bqWKqwHDC7LypUrTVhYWLX6ZmRkmGHDhtVqPQ1VdeexvLzcREVFmWeffdbdVlRUZJxOp/nb3/5WixXWb//+97+NJPPZZ5+52z766CPjcDjM4cOHKx2XkpJiHn300TqosP7q06ePefDBB923y8rKTExMjJk/f36F/W+//XYzZMgQj7akpCQzYcKEWq2zIajpXNbk9dNWksw777xTZZ+pU6earl27erTdcccdJj09vRYra3iqM5cbN240ksx3331XJzXVB6wA1pFNmzapdevW6tKli+6//36dOHHC1yU1KPv371dBQYHS0tLcbWFhYUpKSlJOTo4PK/OtnJwchYeHq3fv3u62tLQ0BQQEKDc3t8qxr776qlq1aqVrrrlGmZmZOnPmTG2XW2+cP39eW7du9Xg+BQQEKC0trdLnU05Ojkd/SUpPT7f6+Sdd3lxKF36RqV27doqNjdWwYcO0a9euuijXr/Cc9L7ExERFR0frpptu0pYtW3xdTq3yi18Cqe8GDhyoW265RR06dNC+ffv0xBNPaNCgQcrJyVFgYKCvy2sQLn4PIzIy0qM9MjLSr7+j8XMKCgou+YiiUaNGatGiRZXzcuedd6pdu3aKiYnRjh07NG3aNO3Zs0erV6+u7ZLrhePHj6usrKzC59OXX35Z4ZiCggKefxW4nLns0qWLVqxYoe7du8vlcum5555T3759tWvXLrVp06YuyvYLlT0ni4uLdfbsWTVp0sRHlTU80dHRWrp0qXr37q3S0lK9/PLLSk1NVW5urq699lpfl1crCICSpk+frgULFlTZZ/fu3UpISLisxx8xYoT7v7t166bu3bsrPj5emzZtUv/+/S/rMeuj2p5Hm1R3Li/Xj78j2K1bN0VHR6t///7at2+f4uPjL/txgepITk5WcnKy+3bfvn111VVXadmyZZo7d64PK4OtunTpoi5durhv9+3bV/v27dPzzz+vV155xYeV1R4CoKTJkydr9OjRVfaJi4vz2v7i4uLUqlUr7d27168CYG3OY1RUlCSpsLBQ0dHR7vbCwkIlJiZe1mPWZ9Wdy6ioqEu+aP/DDz/o5MmT7jmrjqSkJEnS3r17rQiArVq1UmBgoAoLCz3aCwsLK523qKioGvW3xeXM5U81btxYPXv21N69e2ujRL9V2XMyNDSU1T8v6NOnjz7++GNfl1FrCICSIiIiFBERUWf7O3TokE6cOOERZPxBbc5jhw4dFBUVpezsbHfgKy4uVm5ubo3PyG4IqjuXycnJKioq0tatW9WrVy9J0oYNG1ReXu4OddWRl5cnSX73nKxMUFCQevXqpezsbA0fPlySVF5eruzsbD300EMVjklOTlZ2drYmTpzoblu3bp3HSpaNLmcuf6qsrEw7d+7U4MGDa7FS/5OcnHzJpYh4TnpPXl6ef78m+voslIbmwIEDZvv27Wb27NmmadOmZvv27Wb79u3m1KlT7j5dunQxq1evNsYYc+rUKfP444+bnJwcs3//frN+/Xpz7bXXmk6dOplz58756jB8rqbzaIwxzzzzjAkPDzfvvfee2bFjhxk2bJjp0KGDOXv2rC8Ood4YOHCg6dmzp8nNzTUff/yx6dSpkxk5cqT7/kOHDpkuXbqY3NxcY4wxe/fuNXPmzDGff/652b9/v3nvvfdMXFyc6devn68OwSdef/1143Q6TVZWlvn3v/9txo8fb8LDw01BQYExxph77rnHTJ8+3d1/y5YtplGjRua5554zu3fvNjNnzjSNGzc2O3fu9NUh1Bs1ncvZs2ebtWvXmn379pmtW7eaESNGmODgYLNr1y5fHUK9cOrUKfdroSSzaNEis337dnPgwAFjjDHTp08399xzj7v/N998Y0JCQsyUKVPM7t27zZIlS0xgYKBZs2aNrw6h3qjpXD7//PPm3XffNV9//bXZuXOnefTRR01AQIBZv369rw6h1hEAaygjI8NIumTbuHGju48ks3LlSmOMMWfOnDEDBgwwERERpnHjxqZdu3Zm3Lhx7hdGW9V0Ho25cCmYp556ykRGRhqn02n69+9v9uzZU/fF1zMnTpwwI0eONE2bNjWhoaFmzJgxHkF6//79HnObn59v+vXrZ1q0aGGcTqfp2LGjmTJlinG5XD46At958cUXTdu2bU1QUJDp06eP+de//uW+LyUlxWRkZHj0f/PNN03nzp1NUFCQ6dq1q/nnP/9ZxxXXXzWZy4kTJ7r7RkZGmsGDB5tt27b5oOr65eKlSH66XZy7jIwMk5KScsmYxMREExQUZOLi4jxeM21W07lcsGCBiY+PN8HBwaZFixYmNTXVbNiwwTfF1xGHMcbU2XIjAAAAfI7rAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZQiAAAAAliEAAgAAWIYACAAAYBkCIAAAgGUIgAAAAJYhAAIAAFiGAAgAAGAZAiAAAIBlCIAAAACWIQACAABYhgAIAABgGQIgAACAZf4f/ltOg/Jh/ZEAAAAASUVORK5CYII=' width=640.0/>\n",
       "            </div>\n",
       "        "
      ],
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "\n",
    "#Plot the lagrangian points and primaries\n",
    "fig_name = f'{xL_str}_ur{use_real}_dir{direction}'\n",
    "ax.set_title(fig_name)\n",
    "buffer=1.5\n",
    "# ax.set_xlim(xL1-buffer, xL1+buffer)\n",
    "# ax.set_ylim(-buffer, +buffer)\n",
    "ax.set_xlim(-buffer, buffer)\n",
    "ax.set_ylim(-buffer, buffer)\n",
    "ax.scatter(-mu, 0, c='k', s=100)\n",
    "ax.scatter((1-mu), 0, c='k', s=30)\n",
    "ax.scatter(xL1, 0, c='r', s=10)\n",
    "ax.scatter(xL2, 0, c='r', s=10)\n",
    "ax.scatter(xL3, 0, c='r', s=10)\n",
    "ax.scatter(xL4, yL4, c='r', s=10)\n",
    "ax.scatter(xL5, yL5, c='r', s=10)\n",
    "# # zero velocity curve\n",
    "#plot_zero_vel_curves(ax, mu, new_ic, xL, buffer=buffer)\n",
    "# ax.plot(out2[:, 0], out2[:, 1], linewidth=2, c='b')\n",
    "# ax.scatter(out2[-1, 0], out2[-1, 1], s=5)\n",
    "# ax.plot(out3[:, 0], out3[:, 1], linewidth=2, c='r')\n",
    "# ax.scatter(out3[-1, 0], out3[-1, 1], s=5)\n",
    "for i in out.keys():\n",
    "    ax.plot(out[i][:, 0], out[i][:, 1],'g', alpha=0.3) if out[i] is not None else None\n",
    "\n",
    "plt.show()\n",
    "save_fig=False\n",
    "if save_fig:\n",
    "    direc = os.getcwd() + '/figures/'\n",
    "    fig.savefig(direc + f'{fig_name}.png', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3D Plot for Halo etc."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1, projection='3d')\n",
    "\n",
    "#Plot the lagrangian points and primaries\n",
    "fig_name = f'{xL_str}_ur{use_real}_dir{direction}'\n",
    "ax.set_title(fig_name)\n",
    "buffer=1.5\n",
    "# ax.set_xlim(xL1-buffer, xL1+buffer)\n",
    "# ax.set_ylim(-buffer, +buffer)\n",
    "ax.set_xlim(-buffer, buffer)\n",
    "ax.set_ylim(-buffer, buffer)\n",
    "ax.set_zlim(-buffer, buffer)\n",
    "ax.scatter(-mu, 0, 0, c='k', s=100)\n",
    "ax.scatter((1-mu), 0, 0, c='k', s=30)\n",
    "ax.scatter(xL1, 0, 0, c='r', s=10)\n",
    "ax.scatter(xL2, 0, 0, c='r', s=10)\n",
    "ax.scatter(xL3, 0, 0, c='r', s=10)\n",
    "ax.scatter(xL4, yL4, 0, c='r', s=10)\n",
    "ax.scatter(xL5, yL5, 0, c='r', s=10)\n",
    "# # zero velocity curve\n",
    "#plot_zero_vel_curves(ax, mu, new_ic, xL, buffer=buffer)\n",
    "# ax.plot(out2[:, 0], out2[:, 1], linewidth=2, c='b')\n",
    "# ax.scatter(out2[-1, 0], out2[-1, 1], s=5)\n",
    "# ax.plot(out3[:, 0], out3[:, 1], linewidth=2, c='r')\n",
    "# ax.scatter(out3[-1, 0], out3[-1, 1], s=5)\n",
    "plt.show()\n",
    "for i in out.keys():\n",
    "    ax.plot(out[i][:, 0], out[i][:, 1], out[i][:, 2], 'g', alpha=0.3) if out[i] is not None else None\n",
    "save_fig=False\n",
    "if save_fig:\n",
    "    direc = os.getcwd() + '/figures/'\n",
    "    fig.savefig(direc + f'{fig_name}.png', dpi=300)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "%%timeit\n",
    "get_cr3bp_hamiltonian(conj_mom=True)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "%%timeit\n",
    "diff_corr_pa(ta, out[10], new_period, ds, predicted[:,0], taux, tauT, tol=1e-15, conj_mom=True)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "H = get_cr3bp_hamiltonian(conj_mom=True)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "%%timeit\n",
    "get_hamiltonian_state_derivative(H, ta.pars[0], new_ic, conj_mom=True).reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "cf_f = get_compiled_hamiltonian_state_derivative(conj_mom=True)\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "%%timeit\n",
    "cf_f(ta.state[:6], pars=[ta.pars[0]]).reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}