{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cdd23276-df4e-4ffc-8ec7-8ea8286dfd52",
   "metadata": {},
   "source": [
    "# Orbital classify\n",
    "\n",
    "Orbit classification based on the resonance angle behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4eaeecde-704d-4757-9ab2-75d440c6c340",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import galport"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1044be4f-7c14-4c4d-9c6e-904bf4fbf0aa",
   "metadata": {},
   "source": [
    "Load coordinates, velocities and actions for the set of orbits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "6b6e128f-d2a5-4651-b009-a0b1be101cbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "xv_act = np.load('data/xv_act_0.npy')\n",
    "\n",
    "t = np.arange(0, 600.125, 0.125)\n",
    "xv = xv_act[:,:,0:6]\n",
    "act = xv_act[:,:,6:9]\n",
    "\n",
    "theta_p = np.load('data/phi_p.npy')\n",
    "Omega_p = np.load('data/omega_p.npy')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18325f2a-ee15-407c-ba80-007ab7248e89",
   "metadata": {},
   "source": [
    "Find actions, angles and frequencies of the every orbit and Orbit classification at time $t=400$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "63c171ca-31d4-4462-985e-ccc2c6715cc0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.33 s, sys: 6.88 ms, total: 2.34 s\n",
      "Wall time: 2.34 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "OT = galport.OrbitTools(t=t, xv=xv, act=act)\n",
    "data = OT.calculate_actions(secular=True)\n",
    "otype_ILR = OT.classify_orbits(t_out=400, family='ILR', theta_p=theta_p)\n",
    "otype_vILR = OT.classify_orbits(t_out=400, family='vILR', theta_p=theta_p)\n",
    "otype_uha = OT.classify_orbits(t_out=400, family='uha', theta_p=theta_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acf0325c-d983-4095-82dd-878c696e99a7",
   "metadata": {},
   "source": [
    "Equivalent code for classification (calling the ``OrbitClassifier`` module directly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "a70d6be3-32ac-48ab-9b5c-8d208d728b0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "ang = OT.angles\n",
    "# or\n",
    "ang = data[:, :, 3:6] # if dJdt=False\n",
    "# ang = OT.data[:, :, 6:9] # if dJdt=True\n",
    "\n",
    "OC = galport.OrbitClassifier(t, angles=ang, theta_p=theta_p, time_resolution=5.)\n",
    "# Not all angles are used in the classification; they are taken with a time step of `time_resolution`\n",
    "result = OC(t_out=400, family='ILR', time_around_res=False, amplitude_res=False)\n",
    "# or, you can also set the resonant angle yourself\n",
    "OC = galport.OrbitClassifier(t, angles=2*(ang[:,:,2] - theta_p) - ang[:,:,0], time_resolution=5.)\n",
    "result = OC(t_out=400, time_around_res=True, amplitude_res=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "051a71a7-8116-43d4-8787-fc81aa0e54c9",
   "metadata": {},
   "source": [
    "**List of types**\n",
    "\n",
    "`0`: Not classify\n",
    "\n",
    "`1`: Increasing angle\n",
    "\n",
    "`2`: Decreasing angle\n",
    "\n",
    "`3`: Resonance around $0$\n",
    "\n",
    "`4`: Resonance around $\\pi$\n",
    "\n",
    "`5`: Passage through $0$ from $\\Omega_{res} > 0$ to $\\Omega_{res} < 0$\n",
    "\n",
    "`6`: Passage through $\\pi$ from $\\Omega_{res} > 0$ to $\\Omega_{res} < 0$\n",
    "\n",
    "`7`: Passage through $0$ from $\\Omega_{res} < 0$ to $\\Omega_{res} > 0$\n",
    "\n",
    "`8`: Passage through $\\pi$ from $\\Omega_{res} < 0$ to $\\Omega_{res} > 0$\n",
    "\n",
    "**Additional variables**\n",
    "\n",
    "`time_around_res`: returns the resonance entry and exit times for resonant orbits\n",
    "\n",
    "`amplitude_res`: returns the maximum libration amplitude of the resonant angle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7972a978-be2c-44b9-9d24-882198db24e9",
   "metadata": {},
   "source": [
    "Secular actions and frequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "a96343a2-8bd9-4f3d-a7ea-1802f4662068",
   "metadata": {},
   "outputs": [],
   "source": [
    "t0 = 400\n",
    "n_t = int(t0)*8\n",
    "\n",
    "JR_sec = data[:,n_t,9]\n",
    "Jz_sec = data[:,n_t,10]\n",
    "Lz_sec = data[:,n_t,11]\n",
    "\n",
    "kappa_sec = data[:,n_t,12]\n",
    "omegaz_sec = data[:,n_t,13]\n",
    "Omega_sec = data[:,n_t,14] - Omega_p[n_t]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69f5a166-ca3a-47f2-abc6-7e1187127650",
   "metadata": {},
   "source": [
    "Plot different orbital types on the plane $2(\\Omega - \\Omega_p)/\\kappa$ vs $\\omega_z\\kappa$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "b07928f5-9b6d-489f-8b24-75eae7ea5a80",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAFJCAYAAAAmKgzMAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQ1RJREFUeJzt3XlcVOX+B/DPAOMAyiIkMCgIrpkarhjumoQbSota+kvUrt3crkZp2oa0edNSb13TWy7ULZfout20vFxTSUW5Lqi5lYmgBmigIqsD8/z+mBgZWZwZzsyZYT7v14vXcM6c5XueOfCd5znPeY5CCCFAREREknCSOwAiIqKGhImViIhIQkysREREEmJiJSIikhATKxERkYSYWImIiCTExEpERCQhJlYiIiIJucgdgL3TarW4c+eO3GEQETkUpVIJZ2dnucOoERNrPdy5cwcZGRnQarVyh0JE5HC8vb0REBAAhUIhdygGmFjNJIRAdnY2nJ2dERQUBCcntqoTEVmDEALFxcW4du0aAECtVssckSEmVjOVl5ejuLgYgYGBcHd3lzscIiKH4ubmBgC4du0a/Pz8bKpZmNUsM1VUVAAAGjVqJHMkRESOqbJSo9FoZI7EEBNrPdla2z4RkaOw1f+/TKxEREQSYmIlkplCocDWrVslX5bIFjji+c3E6mAmTZoEhUKh//H19cXQoUNx8uRJuUNrcKqWtVKphL+/PyIjI7F27VqDW7Sys7MxbNgwGSOl2ixcuBBdunQxarmqf1deXl7o168f9u3bZ/kgJTRp0iTExMQYtZwx5zbgmOc3E6sDGjp0KLKzs5GdnY3du3fDxcUFI0eONHt7HCCjdpVlfenSJXz33XcYNGgQZs+ejZEjR6K8vBwAEBAQAJVKJXOkVF8dO3bU/12lpqaibdu2GDlyJG7dumX2NisqKmz2Pnljzm3AMc9vJlYHpFKpEBAQgICAAHTp0gXz58/H5cuXcf36dQDAK6+8gnbt2sHd3R2tWrXCG2+8YdDrrvJb/OrVqxEaGgpXV1e5DsXmVZZ18+bN0a1bN7z66qvYtm0bvvvuOyQmJgIwbP66c+cOZs6cCbVaDVdXV7Rs2RKLFi2qdfvx8fFQq9VscajBwIED8Ze//AXz5s2Dj48PAgICsHDhQoNlsrKyMHr0aDRp0gSenp4YO3YscnNzAQCJiYlISEjAiRMn9LWzys+sJi4uLvq/q4ceeghvvfUWCgsL8fPPP+uXWbp0KTp37ozGjRsjKCgI06dPR2Fhof79xMREeHt7Y/v27XjooYegUqmQlZVl1PGGhITgvffew5QpU+Dh4YHg4GB8+umnBsucOnUKgwcPhpubG3x9ffH888/r979w4UJ8/vnn2LZtm/549+7dW+v+jDm3Acc8v5lYHVxhYSG+/PJLtGnTBr6+vgAADw8PJCYm4syZM/jb3/6Gzz77DMuWLTNY78KFC/jXv/6FzZs3Iz09XYbI7dfgwYMRFhaGzZs3V3vvo48+wvbt2/H111/j/Pnz+OqrrxASElJtOSEEZs2ahS+++AI//vgjHn74YStEbn8+//xzNG7cGIcPH8bixYvx1ltvITk5GYBuONLRo0cjPz8f+/btQ3JyMi5evIhx48YBAMaNG4eXXnrJoCZa+d79lJWVYd26dfD29kb79u31852cnPDRRx/h9OnT+Pzzz/HDDz9g3rx5BusWFxfj/fffx+rVq3H69Gn4+fkZfbwffvghevTogePHj2P69OmYNm0azp8/DwAoKipCVFQUmjZtiv/9739ISkrCf//7X8ycORMA8PLLL2Ps2LEGLVq9e/c2et9A3ec24EDntyCzlJSUiDNnzoiSkpJ6byvz9yKRdOSyyPy9SILI6hYbGyucnZ1F48aNRePGjQUAoVarxdGjR2tdZ8mSJaJ79+766fj4eKFUKsW1a9csHq/k8i4Kcfwr3auFxcbGitGjR9f43rhx40SHDh2EEEIAEFu2bBFCCDFr1iwxePBgodVqa1wPgEhKShLjx48XHTp0EFeuXLFE6BaTVZAltv6yVWQVZFl8XwMGDBB9+/Y1mNezZ0/xyiuvCCGE+M9//iOcnZ1FVtbdWE6fPi0AiLS0NCGE7lwPCwu7777i4+OFk5OT/u9KoVAIT09P8d1339W5XlJSkvD19dVPr1u3TgAQ6enpxh6mXsuWLcX//d//6ae1Wq3w8/MTK1euFEII8emnn4qmTZuKwsJC/TI7duwQTk5OIicnRwhR9zlblbHnthCWPb+l/D8sJY68JLOsvGJELU9BiaYCbkpn7JrTH8G+lh3JadCgQVi5ciUA4MaNG/jkk08wbNgwpKWloWXLlti0aRM++ugj/PrrrygsLER5eTk8PT0NttGyZUs0a9bMonFKLj8DWNkb0BQDSndg2kHAJ1SWUIQQNd6DN2nSJERGRqJ9+/YYOnQoRo4ciccee8xgmRdffBEqlQqHDh3CAw88YK2Q6+3y7ct4YvsTKC0vhauLKzaP2owgjyCL7vPemo5ardYPg3f27FkEBQUhKOhuDA899BC8vb1x9uxZ9OzZ06R9tW/fHtu3bwcA3L59G5s2bcKYMWOwZ88e9OjRAwDw3//+F4sWLcK5c+dQUFCA8vJylJaWori4WD/YQaNGjcyuoVVdT6FQICAgwOB4w8LC0LhxY/0yffr0gVarxfnz5+Hv72/WPu9V27kNNOzzuyo2Bcss7VI+SjS6UZxKNBVIu5Rv8X02btwYbdq0QZs2bdCzZ0+sXr0aRUVF+Oyzz5CamooJEyZg+PDh+Pbbb3H8+HG89tpr1TooVf3jtBtZqbqkCuhes1JlC+Xs2bMIDa2e1Lt164aMjAy8/fbbKCkpwdixY/HUU08ZLBMZGYmrV69i165d1gpXEsdyj6G0vBQAUFpeimO5xyy+T6VSaTCtUCgs1hmoUaNG+r+rrl274q9//SuaN2+O5cuXAwAuXbqEkSNH4uGHH8a//vUvHD16FCtWrABg2AHQzc3N7IEPrHm8tant3AYa9vldFWusMgsP8YGb0llfYw0P8bF6DAqFAk5OTigpKcHBgwfRsmVLvPbaa/r3MzMzrR6TRQRH6GqqlTXW4AhZwvjhhx9w6tQpvPjiizW+7+npiXHjxmHcuHF46qmnMHToUOTn58PHR3dujBo1CtHR0Rg/fjycnZ3x9NNPWzN8s3Xz7wZXF1d9jbWbfzdZ4+nQoQMuX76My5cv62utZ86cwc2bN/HQQw8B0CXLyuFLzeHs7IySkhIAwNGjR6HVavHhhx/qH9rx9ddf1/MojNehQwckJiaiqKhI/8X4wIEDcHJy0l8Hru/x3u/cBhru+V0VE6vMgn3dsWtOf6Rdykd4iI/Fm4EBXceKnJwcALqm4L///e8oLCxEdHQ0CgoKkJWVhY0bN6Jnz57YsWMHtmzZYvGYrMInVNf8m5WqS6pWaAauLOuKigrk5ubi+++/x6JFizBy5EhMnDix2vJLly6FWq1G165d4eTkhKSkJAQEBMDb29tguccffxz//Oc/8eyzz8LFxaXat35bFOQRhM2jNuNY7jF08+9m8Wbg+xkyZAg6d+6MCRMmYPny5SgvL8f06dMxYMAAfdNtSEgIMjIykJ6ejhYtWsDDw6PWW0fKy8v1f1eVTcFnzpzBK6+8AgBo06YNNBoNPv74Y0RHR+PAgQNYtWqVdQ4WwIQJExAfH4/Y2FgsXLgQ169fx6xZs/Dss8/qm4FDQkKwa9cunD9/Hr6+vvDy8qpWC65k6rkNNOzz24DcF3ntla1eNL+f2NhYAUD/4+HhIXr27Cm++eYb/TJz584Vvr6+okmTJmLcuHFi2bJlwsvLS/++sR06HF3VsnZxcRHNmjUTQ4YMEWvXrhUVFRX65VClc8enn34qunTpIho3biw8PT3Fo48+Ko4dO1bjskIIsWnTJuHq6ir+9a9/Weuw7MaAAQPE7NmzDeaNHj1axMbG6qczMzPFqFGjROPGjYWHh4cYM2aMviOPEEKUlpaKJ598Unh7ewsAYt26dTXuKz4+3uDvyt3dXXTu3FnfcajS0qVLhVqtFm5ubiIqKkp88cUXAoC4ceOGEELXeanq31qlPXv2CAAiIyOj1uNt2bKlWLZsmcG8sLAwER8fr58+efKkGDRokHB1dRU+Pj5i6tSp4vbt2/r3r127JiIjI0WTJk0EALFnz54a92XsuS2EZc9vW/0/rBBCCGsn84agtLQUGRkZvI+TiCxu3bp1eO+993DmzJlaa5COyFb/D7PzEhGRjdu5cyfee+89JlU7wWusREQ2LikpSe4QyASssRIREUmIiZWIiEhCTKz1xL5fRETysNX/v7Im1kWLFqFnz57w8PCAn58fYmJi9ANG1yUpKQkPPvggXF1d0blzZ+zcudMK0RpydnYGwEemERHJpbhYN5KarXXqkrXz0r59+zBjxgz07NkT5eXlePXVV/HYY4/hzJkztQ6Zd/DgQTzzzDP6G5HXr1+PmJgYHDt2DJ06dbJa7C4uLnB3d8f169ehVCr1I6kQEZFlCSFQXFyMa9euwdvbW1/RsRU2dR/r9evX4efnh3379qF///41LjNu3DgUFRXh22+/1c975JFH0KVLF6uOYgLoaqsZGRk2+yBiIqKGzNvbGwEBAWaPrWwpNnW7za1btwBAP2ZkTVJTUxEXF2cwLyoqSv8gXWtq1KgR2rZty+ZgIiIrUyqVNldTrWQziVWr1WLOnDno06dPnU26OTk51R5v5O/vrx+jsyZlZWUoKysz2Fd+fj58fX1t7psOERHdnxQVGiEEbt++jcDAQEkv59lMYp0xYwZ++ukn7N+/X/JtL1q0CAkJCZJvl4iI7N/ly5fRokULybZnE4l15syZ+Pbbb5GSknLfgwsICEBubq7BvNzcXAQEBNS6zoIFCwyaj2/duoXg4GD8/PPPdTY7010ajQZ79uzBoEGDbK4Hni2TpdxObwG+X3B3eugioOPj1tm3BHiumYflZrr8/Hy0a9cOHh4ekm5X1sQqhMCsWbOwZcsW7N27t9aH41YVERGB3bt3Y86cOfp5ycnJiIio/dmaKpWqxkc9+fj4wNfX16zYHY1Go4G7uzt8fX35R2sCWcqt06PAjy53nzvb6VHAx37Oc55r5mG5mU/qS4KyJtYZM2Zg/fr12LZtGzw8PPTXSb28vODm5gYAmDhxIpo3b45FixYBAGbPno0BAwbgww8/xIgRI7Bx40YcOXIEn376qWzHQWRTZHjuLBHdJWtiXblyJQBg4MCBBvPXrVuHSZMmAQCysrIMLir37t0b69evx+uvv45XX30Vbdu2xdatW616DyuRzfMJZUIlkonsTcH3s3fv3mrzxowZgzFjxlggIiIiovrhcEFEREQSYmIlIiKSEBMrERGRhJhYiYiIJMTESkREJCEmViIiIgkxsRLJIT8DSF+veyWiBsUmxgomcij5GcDK3neHHJx2kIM5EDUgrLESWVtWqi6pArrXrFR54yEiSTGxEllbcISupgroXoNrf4AEEdkfNgUTWRsHySdq0JhYieTAQfKJGiw2BRMREUmIiZWIiEhCTKxEREQSYmIlIiKSEBMrERGRhJhYiYiIJMTESkREJCEmViIiIgnJnlhTUlIQHR2NwMBAKBQKbN269b7rfPXVVwgLC4O7uzvUajWmTJmCvLw8ywdLRER0H7In1qKiIoSFhWHFihVGLX/gwAFMnDgRzz33HE6fPo2kpCSkpaVh6tSpFo6UiIjo/mQf0nDYsGEYNmyY0cunpqYiJCQEf/nLXwAAoaGh+POf/4z333/fUiESEREZTfYaq6kiIiJw+fJl7Ny5E0II5Obm4ptvvsHw4cPlDo2IiEj+Gqup+vTpg6+++grjxo1DaWkpysvLER0dXWdTcllZGcrKyvTTBQUFAACNRgONRmPxmBuCynJieZmG5WY6lpl5WG6ms1RZKYQQwiJbNoNCocCWLVsQExNT6zJnzpzBkCFD8OKLLyIqKgrZ2dmYO3cuevbsiTVr1tS4zsKFC5GQkFBt/vr16+Hu7i5V+EREZEeKi4sxfvx43Lp1C56enpJt1+4S67PPPovS0lIkJSXp5+3fvx/9+vXDb7/9BrVaXW2dmmqsQUFByM7Ohq+vr6TH0FBpNBokJycjMjISSqVS7nDsBsvNdCwz87DcTJeXlwe1Wi15YrW7puDi4mK4uBiG7ezsDACo7TuCSqWCSqWqNl+pVPIENBHLzDw2V275GTb/oHWbKzM7wXIznqXKSfbEWlhYiAsXLuinMzIykJ6eDh8fHwQHB2PBggW4evUqvvjiCwBAdHQ0pk6dipUrV+qbgufMmYPw8HAEBgbKdRhE9iM/A1jZG9AUA0p3YNpBm02uRPZI9sR65MgRDBo0SD8dFxcHAIiNjUViYiKys7ORlZWlf3/SpEm4ffs2/v73v+Oll16Ct7c3Bg8ezNttyDLsoGZnsqxUXVIFdK9ZqQ3n2IhsgOyJdeDAgbU24QJAYmJitXmzZs3CrFmzLBgVERpuzS44Qnc8lccVHCF3REQNiuyJlchmNdSanU+o7ktCQ6uJE9kIJlai2jTkmp1PKBMqkYUwsRLVhjU7IjIDEytRXVizIyIT2d1YwURERLaMiZWIiEhCTKxElnQj0/CViBo8JlayL/kZQPp63auty88A1kTqfl8TaR8xm8sSn8uNTPv5rImqYOclsh/2NmBDQ70P9l6W+lzWRAJl+fbxWRNVwRor2Y+aEpUtq7wPFpDnPlhr1e4t9bnY02dNVAVrrGQ/7G3ABp9Q4LlkIPW07tWaNS5r1u4t9bko3YGyUvv4rImqYGIl+2GPAzY0bQng9B+vVmTNZmhLfS7PJQO/pdnPZ030ByZWsi8csME41q7dW+JzadoS8Gsj7TaJrICJlag+jH2s3KkkIMSKNa/KWuTZ7dbZHxHpMbESmcuY65iV969++yLg7GT93q17/6qLb+9f2bOWyErYK5jIXMb0hr2Sdvd3a/dutbde1EQNBBMrkbmMuZ2mRfjd311cpb3Web/baeS+3YfIQbEpmMhccvZSNqYZurb47ndd2NjrxkRUIyZWovq4X2/YK2kAGut+Ly+V7rYXY2+nuTe++yVkexvdisgGsSmYyJKqNgVL2RxrbjPv/a67nt3O67JE9SR7Yk1JSUF0dDQCAwOhUCiwdevW+65TVlaG1157DS1btoRKpUJISAjWrl1r+WCJTFU5MMTIZdLW/iqbeWNWmrbduhJyfgaw572701JfEyZyELI3BRcVFSEsLAxTpkzBE088YdQ6Y8eORW5uLtasWYM2bdogOzsbWq3WwpES1UPnMYBSKe02zRmUoa7rwlmpuubqSoNeZTMwkRlkT6zDhg3DsGHDjF7++++/x759+3Dx4kX4+PgAAEJCQiwUHVEDVFtCvne0pg6jrB8bUQMge1OwqbZv344ePXpg8eLFaN68Odq1a4eXX34ZJSUlcodGJC1rP3vW3OZlIjIge43VVBcvXsT+/fvh6uqKLVu24Pfff8f06dORl5eHdevW1bhOWVkZysrK9NMFBQUAAI1GA41GY5W47V1lObG8TGN2ud3I1D2PtLL2+FyydQby92gBdByj+12mz5rnmnlYbqazVFkphBDCIls2g0KhwJYtWxATE1PrMo899hh+/PFH5OTkwMvLCwCwefNmPPXUUygqKoKbm1u1dRYuXIiEhIRq89evXw93d3fJ4iciIvtRXFyM8ePH49atW/D09JRsu3ZXY1Wr1WjevLk+qQJAhw4dIITAlStX0LZt22rrLFiwAHFxcfrpgoICBAUFYdCgQfD19bVK3PZOo9EgOTkZkZGRUErdCacBu2+51VYztUaNVa5a8X3wXDMPy810eXl5Ftmu3SXWPn36ICkpCYWFhWjSpAkA4Oeff4aTkxNatGhR4zoqlQoqlarafKVSyRPQRCwz89Rabr+lAWX5ut/LSnXTfm10P3/+wbIjINW2bxvBc808LDfjWaqcZO+8VFhYiPT0dKSnpwMAMjIykJ6ejqysLAC62ubEiRP1y48fPx6+vr6YPHkyzpw5g5SUFMydOxdTpkypsRmYyKbVdV+pTyjQZbzlOhFxLGEii5C9xnrkyBEMGjRIP13ZZBsbG4vExERkZ2frkywANGnSBMnJyZg1axZ69OgBX19fjB07Fu+8847VYyeqNznHG5Zz30QNmOyJdeDAgair/1RiYmK1eQ8++CCSk5MtGBWRFZkz0END2DdRAyV7UzAREVFDwsRKDYO1B1MgIqqF7E3BRPXGR50RkQ1hjZXs3/0ehUas0RNZEWusZP/uHTyet40YsrUafX4GeyJTg8bESvaPt43UraYavVxlZGtJnsgCmFipYbCH20bkqqnZUo3elpI8kYUwsRJZw41M4LN+hjU1wDqJVq4afU1fJGwpyRNZCBMrkTVcSTOsqZ3dDuz9q/WaRK1do6+tyZfN9uQAmFiJrKFFuGFNDWjYTaJ1NfnaQ7M9UT0wsRJZQ9OWhjU1wLDGam9Nove7XswmX3JgTKxE1nJvTc1em0SN6dnLJl9yYEysRHKx1yZRY3v22uvxEdUTR14isgX2NDISn+NKVCfWWInkZm+DJrCZl6hOTKxEcrPHQRPYzEtUKzYFk/2wp+ZSU7BplahBYY2V7IO9NZeagk2rRA0KEyvZB3tsLr2fe+8FtffjsRQ+DYfsDBMr2YeqAw44NwK8guSO6P7yM4BLqQAa1/xeQ62BS8lWysmOkvuV/BIcvZKL8BAfBPu6yx2OQ5L9GmtKSgqio6MRGBgIhUKBrVu3Gr3ugQMH4OLigi5dulgsPrIRPqHA+K91SbXiDrB+rG1fa834EVgRDnz7om76Rqbh+2e3N7yHs1viGrjUD7E3J8bK5L51mu7Vls87AI9/cgAvJ51A1PIUZOUVyx2OQ5I9sRYVFSEsLAwrVqwwab2bN29i4sSJePTRRy0UGdmcW5d1SRWw7WSUnwF8+cTdWAHdIPxV39/z3t1pF1f777BkqeQjZccuc2OUOrlbWEl5he5VU4G0S/kyR+OYZG8KHjZsGIYNG2byei+88ALGjx8PZ2dnk2q5ZMfsZfzZrFTDpAroBuGv+n556d3pQa/afPPifVnqGriUHbvMjdFezrs/uLk4o6xCCzelM8JDfOQOxyHJnljNsW7dOly8eBFffvkl3nnnHbnDIWuxl96z914PBnSD8Nf0vtId6DBKnjilZMnkI1XHLnNjtJfz7g9bpvfB0SsFvMYqI7tLrL/88gvmz5+PH3/8ES4uxoVfVlaGsrIy/XRBQQEAQKPRQKPRWCTOhqaynOpVXjcydU2iLcINE40pPFoAHcdUBmV+LJbk0QKY+iNwJQ2agO7A/84ZlluV99EiXDdtq8diLAmPSZJzrSb1ibGm806K81lCleXl7+GC0Q/7G8yjmlmqfOwqsVZUVGD8+PFISEhAu3btjF5v0aJFSEhIqDZ/z549cHfnNzpTJCcn13MLjYHLpwGcliIcG9cYuHwOQG3l1hDLQrpjqv+5Vhspy932PkPLlVvDU1xsmc5dCiGEsMiWzaBQKLBlyxbExMTU+P7NmzfRtGlTODs76+dptVoIIeDs7Iz//Oc/GDx4cLX1aqqxBgUFITs7G76+vpIfh12r5Vu4RqNBcnIyIiMjoVQqTd/uqaS7PWQBYOQyoPMYCQK2bfUuN1t2IxNYE3m3afW5ZElqbnZRZjZ4PttFudmYvLw8qNVq3Lp1C56enpJt165qrJ6enjh16pTBvE8++QQ//PADvvnmG4SG1nztQ6VSQaVSVZuvVCp5AlaVnwF81q/OewbNLrOQCMDZ6e62QyIAByr7Bnmu/ZYGlP3R67SsVDft10ayzdt0mdnw+WzT5WZjLFVOsifWwsJCXLhwQT+dkZGB9PR0+Pj4IDg4GAsWLMDVq1fxxRdfwMnJCZ06dTJY38/PD66urtXmkxksObqRnXUAISPYWW9ZSfF8pjrInliPHDmCQYMG6afj4uIAALGxsUhMTER2djaysrLkCs+xBEfo7qksL7XMvZVS9e60o1FwGjRHTy4chpJqIXtiHThwIOq6zJuYmFjn+gsXLsTChQulDYpsl60McUc6TC5E1cg+8hLZkKoDF5SX2uYIM3Y2Cg4ROR4mVrrLHp4Lag8xEpFDk70pmGyIPVwzs/UYef2XbEhWXjHSLuVzFCYrY2IlQ6ZcM5Mrich1Xe9+x1vT9V+PFtaPkwi6pBq1PAUlmgq4KZ2xa05/JlcrYWIl8zhaJyJjjrem678dG/4gGGSb0i7lo0Rj+KQbJlbr4DVWMo+jdSIy5nh5/ZdsSHiID9yUulHq+KQb6zI6sb755puWjIPsjaMlEWOOt/L6b8zKhl+DJ5tzJb/E4DXY1x275vTHB2PC2AxsZUYn1qSkJHz11VfV5m/evBkPP/ywpEGRHXC0JGLs8fqEAl3GN/zysAf5GUD6eukeum7DsvKK8fgnBwAAj39yAFl5utaVYF93PNW9BZOqlRl9jfXf//43Bg4ciJCQEPTp0wfbt2/HwoULceLECTzxxBOWjJFslaMNDtAQjtdRei07WB+A737KRkn5H9dTy3k9VW5GJ9Y2bdpg48aNGDduHNRqNY4dO4annnoKX3zxBcfpJbIHjpRsLDnutY3JyivGsuSf9dMqZydeT5WZ0U3Bqamp6NatGxYvXoyTJ09i//792LRpE5Mqkb1wpA5nlugDYKNNy2mX8lFartVPzxzchrVVmRldY+3Tpw+cnJzQunVrqFQqvP7664iLi0NYWBhatOC9ekQ2z5GeRiP1QCI2XNuv7P2r1ZYDACI7BMgcERmdWG/cuIGTJ0/ixIkT+p8xY8agtLQUTZs2RV5eniXjJKL6svVRq6Qm5TVxG25aruz9m3bxOvBbOlr4uMkdksMzOrF6eXmhX79+6Nevn36eVqvF+fPncfLkSYsER0QSawgdsORg47X9YF93qD0DsfO3dGxL/w0tfJrgys0SDmUoE7NHXjp27Bi6deuGDh06oEOHDlLGRERkW2y0tl91LOCKcl1T8GtbT6GsQgEAHMpQJmaPvBQeHq5/KHmlnTt31jsgIiKrM6Zjko3do1w5FvDLSScQtTwFyWdzqi1TOZQhWZfZibVz587w9PTE5MmT9fNef/11SYIiIrK4ymSa8aOuY9LWabpXG+v1W5t7xwKuCYcylIfZTcEKhQILFy7E8uXL8dRTT2HDhg0QQkgZGxGRZVTt5evcCKi4o5tvYx2T6lLZG7jy6TWRHQJw8tBZvBvTmddYZWZ2YvX09AQAzJkzB02bNsWoUaNQUlIiWWBEDcKNzLuvfm3kjYXuqtrLt+LO3eRqgx2TaqPvDXzPNdbuwU0R6u8pc3SOzeym4L179+p/j42NxfPPP49r165JERNRw5CfAayJ1P2+JtJumhgdgleQ4fToT+xy3OvKsYAB1DhWMMnD6MS6ePFijB8/HoMHD8aIESMwb948pKbeHbnl8ccfR36+6RfJU1JSEB0djcDAQCgUCmzdurXO5Tdv3ozIyEg0a9YMnp6eiIiIwK5du0zeL5HFOdJIR/bm1mXDaa3GpjommSrtUn61sYJJPkYn1o8//hi///47/Pz8AAAbNmxA3759MXToUNy6dcvsAIqKihAWFoYVK1YYtXxKSgoiIyOxc+dOHD16FIMGDUJ0dDSOHz9udgxEFuFoj9azJw3sswkP8YGbyx/PXnVhhyW5GX2N9fLly9XmHTp0CNOmTcOMGTPw5ZdfmhXAsGHDMGzYMKOXX758ucH0e++9h23btuHf//43unbtalYMRBbhEwo8lwyknta92mltqEGy0ftSzRXs645PJnTD9bOH8MmEbuywJDOzr7ECwCOPPIJ169Zh+/btUsVjMq1Wi9u3b8PHh9/QyAY1bWn4SrbDxu5LrY+svGJM/+oYAGD6V8d4jVVmZvUKXrduHTw8PODq6oqtW7fC19dX6riM9sEHH6CwsBBjx46tdZmysjKUlZXppwsKCgAAGo0GGo3G4jE2BJXlxPIyDcvNdCwz06VdvK4fhF+rLUfaxetQewbKHJXts9Q5ZlZiPXz4MJKSknDz5k2MGDFCthrr+vXrkZCQgG3btumv/dZk0aJFSEhIqDZ/z549cHdnk4kpkpOT5Q7BLrHcTMcyM54SwNs9dL+/3UML/JaOnb+lyxmSXSgutkzNXiHMHNVBCIHvv/8ec+fOxdy5cxEbG1v/YBQKbNmyBTExMfddduPGjZgyZQqSkpIwYsSIOpetqcYaFBSE7OxsWWvb9kSj0SA5ORmRkZFQKpVyh2M3WG6mu2+Z3cgErqQBLcLZxF5F5rXbOH3kR3Ts0Q8t/TzkDscu5OXlQa1W49atW/qxGaRgdI21f//+WLJkCXr16gVAlwSHDRsGtVqN4cOHS5JYjbVhwwZMmTIFGzduvG9SBQCVSgWVSlVtvlKp5D87E7HMzMNyM12NZZafAXzWzyafiyq3ln4eOP3HK88141iqnIxOrB07dkSfPn0QHh6OJ598Ep07d0aTJk2wYcOGeo24VFhYiAsXLuinMzIykJ6eDh8fHwQHB2PBggW4evUqvvjiCwC65t/Y2Fj87W9/Q69evZCToxt42s3NDV5eXmbHQUR2wIafiyq3K/kl+tdQfyZWORndK3jlypU4ceIE2rVrh7feegtDhw5F37598cknn2D+/PlmB3DkyBF07dpVf6tMXFwcunbtijfffBMAkJ2djaysLP3yn376KcrLyzFjxgyo1Wr9z+zZs82OgYjshNz3nxrzFBwZZOUVc+QlG2JS56WOHTsiMTERa9aswa+//oqbN2+iZcuW8Pf3NzuAgQMH1jl4f2JiosF01aEUicjBSHX/aX6G6duoOnC/jTVD1zTyEu9llY9ZvYKdnZ3Rrl07qWMhIro/n9D6JTRzE6QNN0PfHXmpgiMv2YB6DRBBRGR3zB3DWe5m6DoE+7pjy/Q+AIAt0/uwtiozsx8bR0RklyoTZGWN1dgEaePDILbwccPJP15JXkysRORY6pMg69sMTQ6BiZWIHE8DSpBZecVIu5SP7i34cHNbwcRKRGSnsvKKEbU8BSWaCnirnJDQTe6ICGDnJaoPG72nj8hRpF3KR4nm7m02ZBtYYyXz2PA9fdSAmHO/qb2pxzGGh/jATemMEk2F/nabqiqbicNDfNhT2IqYWMk8NnxPH9kAKRLijcyGPy5wPb+gBvu6Y9ec/vprrCcP7dG/V7WZ2E3pjF1z+jO5WgkTK5nH3FsWqOGTqjXjSpr0X95srQYswRfUYF93BPu6Q6PR4GSV+QbNxBqOxmRNTKxkHhu/p49kJFVrRotwab+82eLlCwt+QTVoJlZyNCZrYmIl8zWgWxZIQlIli6Ytpf3yZouXLyz4BbVqMzGvsVoXEysRSUvKZCHllzdbvXxhwS+olc3EZF1MrEQkPVtszeDlC7ISJlYichy2mPCpweEAEURERBJiYiUiIpIQEysREZGEmFiJiIgkxMRKREQkIdkTa0pKCqKjoxEYGAiFQoGtW7fed529e/eiW7duUKlUaNOmDRITEy0eJxERkTFkT6xFRUUICwvDihUrjFo+IyMDI0aMwKBBg5Ceno45c+bgT3/6E3bt2mXhSImMxMfpETk02e9jHTZsGIYNG2b08qtWrUJoaCg+/PBDAECHDh2wf/9+LFu2DFFRUZYKk8g4945HO/VHuSMiB3Elv0T/GuqvlDkaxyZ7jdVUqampGDJkiMG8qKgopKamyhQRURX3jkd7JU3eeMghZOUV4/FPDgAAHv/kALLyimWOyLHJXmM1VU5ODvz9/Q3m+fv7o6CgACUlJXBzc6u2TllZGcrKyvTTBQUFAACNRgONRmPZgBuIynJied1HYDig8tHXWDUB3YHL51huJuC5Zrq0i9eh1ZYDALTacqRdvA61Z6DMUdk+S51jdpdYzbFo0SIkJCRUm79nzx64u3OAalMkJyfLHYLte2j53d//dw4Ay80cLDPjKQG83UP3+9s9tMBv6dj5W7qcIdmF4mLL1OztLrEGBAQgNzfXYF5ubi48PT1rrK0CwIIFCxAXF6efLigoQFBQEAYNGgRfX1+LxttQaDQaJCcnIzIyEkolr98Yi+VmOpaZ6WLXpeGny/l4u4cWbxxxQqcgH3w+OVzusGxeXl6eRbZrd4k1IiICO3fuNJiXnJyMiIjaHwGlUqmgUqmqzVcqlfzDNRHLzDwsN9OxzIz3aAc1jmbeAACUaRV4tIOaZWcES5WR7J2XCgsLkZ6ejvT0dAC622nS09ORlZUFQFfbnDhxon75F154ARcvXsS8efNw7tw5fPLJJ/j666/x4osvyhE+EZHsvNwb1TlN1iV7Yj1y5Ai6du2Krl27AgDi4uLQtWtXvPnmmwCA7OxsfZIFgNDQUOzYsQPJyckICwvDhx9+iNWrV/NWGyJq0LLyivHN0Ss19vgND/GBm4szAMDNxRnhIT7WDo+qkL0peODAgRBC1Pp+TaMqDRw4EMePH7dgVEREtiMrrxhRy1NQoqmA0lmBf07phUda3+0fEuzrji3T++DkoT3YMr0Pgn3ZKVNOstdYiYiobmmX8lGiqQAAaCoEnl17uFrNtYWPrvPm0awb+vfqquWS5cheYyUiorqFh/hA6ayApkLXuqepEEi7lG9QM60ceem1rafg5OSCdZN6YnLi/1CiqYCb0hm75vRnTdZKWGMlItuRnwGcSpI7CpsT7OuOf07pBaWzAgDgpjS8jpqVV4zP9l/UT5doKrAl/aq+lluiqUDapXzrBu3AWGMlItPlZ+iGbwyOAHxCpdvmyt5AhRYI+xS4kQn4tZFm2w3AI619sTtuINIu5SM8xEdf+6y8/qrVlqPXH7euuimd8XiX5tie/pu+xsoOTdbDxEpEprn3QQPTDkqTXCvHWXZy1U1fSWNivUewr3u15tzK668qXadgPNmtBV4Y2A7Bvu7YNad/tURMlsfESkSmufdBA1mp0iTW4Ahdoq7Q6qZbcOQgY4SH+MBN6awfK3hkZ7U+idaUiMnyeI2ViExTmQAB3Wtw7aOemcQnVFf7HblMN920pTTbbeCCfd2xblJPKJ1011+f/+cRHPrVMkP1kXFYYyUi01QmQKmvsVZu26MFcHnn/ZclZOUVI+1SPn69VgiN9o8ew1qB/1tzGD+8NJC1VZkwsRKR6XxCpU2oZLKqg0a4OAF/dBgGAJRrq9+OQ9bDpmAiIjtUddCIcq1hYlW5OLEXsIyYWImI7FBlpyVAl0hdnHT/zpVOCnw+OZy1VRmxKZiIyA5VvZ0mr7AMS/9zFoDuGuuVmyUyR+fYmFiJiGzA5duXcSz3GLr5d0OQR5BR61TeTnNvL+AW3m6WCJGMxMRKRCSzy7cv44ntT6C0vBSuLq7YPGqz0ck1K68YW9KvGsxjjVVeTKxERDI7lnsMpeWlAIDS8lIcyz1mVGKt2jO4cuQlPo9Vfuy8REQks27+3eDqohvK0dXFFd38uxm1XtWewZX4PFb5scZKRCSzII8gbB612eRrrJU9g0s0FXBzcQZQoX8uK8mHiZWIyAYEeQQZnVArVe0Z3L2FJ04e2mOh6MgUTKxERHassmewRqPBSbmDIQC8xkpEZBey8orxzdEryMorljsUug+bSKwrVqxASEgIXF1d0atXL6SlpdW5/PLly9G+fXu4ubkhKCgIL774IkpLS60ULRGRdVX2/n056QSilqcwudo42RPrpk2bEBcXh/j4eBw7dgxhYWGIiorCtWvXalx+/fr1mD9/PuLj43H27FmsWbMGmzZtwquvvmrlyImIrKNq798STQXSLuXLHBHVRfbEunTpUkydOhWTJ0/GQw89hFWrVsHd3R1r166tcfmDBw+iT58+GD9+PEJCQvDYY4/hmWeeuW8tl4jIXlUdF9hNyftUbZ2sifXOnTs4evQohgwZop/n5OSEIUOGIDU1tcZ1evfujaNHj+oT6cWLF7Fz504MHz7cKjETEVna5duXse3CNly+fRnA3d6/H4wJw645/Xmfqo2TtVfw77//joqKCvj7+xvM9/f3x7lz52pcZ/z48fj999/Rt29fCCFQXl6OF154oc6m4LKyMpSVlemnCwoKAAAajQYajUaCI2n4KsuJ5WUalpvpbKXMrhZexYlrJxDmF4bmTZpbdb8Tdk5AWXkZVC4qfDX8KzRv0hxqTyVGP6z7X1lT2dhKudkTS5WVQgghLLJlI/z2229o3rw5Dh48iIiICP38efPmYd++fTh8+HC1dfbu3Yunn34a77zzDnr16oULFy5g9uzZmDp1Kt54440a97Nw4UIkJCRUm79+/Xq4u/ObHxGRIyouLsb48eNx69YteHp6SrZdWRPrnTt34O7ujm+++QYxMTH6+bGxsbh58ya2bdtWbZ1+/frhkUcewZIlS/TzvvzySzz//PMoLCyEk1P11u2aaqxBQUHIzs6Gr6+vtAfVQGk0GiQnJyMyMhJKpVLucOwGy810tlBmOy/uxNuH3tZPv/HIGxjeyjqXm2qrsd6PLZSbvcnLy4NarZY8scraFNyoUSN0794du3fv1idWrVaL3bt3Y+bMmTWuU1xcXC15OjvrLurX9h1BpVJBpVJVm69UKnkCmohlZh6Wm+nkLLNugd2gcFHonzbTLbCb1WIJaRqCjaM2mjy8YSWea8azVDnJPvJSXFwcYmNj0aNHD4SHh2P58uUoKirC5MmTAQATJ05E8+bNsWjRIgBAdHQ0li5diq5du+qbgt944w1ER0frEywRUX2YO3avlPu39j5JOrIn1nHjxuH69et48803kZOTgy5duuD777/Xd2jKysoyqKG+/vrrUCgUeP3113H16lU0a9YM0dHRePfdd+U6BCJqgJjcyFyyJ1YAmDlzZq1Nv3v37jWYdnFxQXx8POLj460QGRERkWlkHyDCoeRnAOnrda9ERNQg2USN1SHkZwArewOaYkDpDkw7CPiEyh0VERFJjDVWa8lK1SVVQPeaVfPIUkREZN+YWK0lOEJXUwV0r8ERdS9PRER2iU3B1uITqmv+zUrVJVU2AxMRNUhMrNbkE8qESkTUwLEpmIiISEJMrERERBJiYiUiIpIQEysREZGEmFiJiIgkxMRKREQkISZWIiIiCTGxEhERSYiJlYiISEJMrERERBJiYiUiIpIQEysREZGEmFiJiIgkxMRKREQkIZtIrCtWrEBISAhcXV3Rq1cvpKWl1bn8zZs3MWPGDKjVaqhUKrRr1w47d+60UrRERES1k/15rJs2bUJcXBxWrVqFXr16Yfny5YiKisL58+fh5+dXbfk7d+4gMjISfn5++Oabb9C8eXNkZmbC29vb+sETERHdQ/bEunTpUkydOhWTJ08GAKxatQo7duzA2rVrMX/+/GrLr127Fvn5+Th48CCUSiUAICQkxJohExER1UrWpuA7d+7g6NGjGDJkiH6ek5MThgwZgtTU1BrX2b59OyIiIjBjxgz4+/ujU6dOeO+991BRUWGtsImIiGola431999/R0VFBfz9/Q3m+/v749y5czWuc/HiRfzwww+YMGECdu7ciQsXLmD69OnQaDSIj4+vcZ2ysjKUlZXppwsKCgAAGo0GGo1GoqNp2CrLieVlGpab6Vhm5mG5mc5SZSV7U7CptFot/Pz88Omnn8LZ2Rndu3fH1atXsWTJkloT66JFi5CQkFBt/p49e+Du7m7pkBuU5ORkuUOwSyw307HMzMNyM15xcbFFtitrYn3ggQfg7OyM3Nxcg/m5ubkICAiocR21Wg2lUglnZ2f9vA4dOiAnJwd37txBo0aNqq2zYMECxMXF6acLCgoQFBSEQYMGwdfXV6Kjadg0Gg2Sk5MRGRmpv7ZN98dyMx3LzDwsN9Pl5eVZZLuyJtZGjRqhe/fu2L17N2JiYgDoaqS7d+/GzJkza1ynT58+WL9+PbRaLZycdJeIf/75Z6jV6hqTKgCoVCqoVKpq85VKJU9AE7HMzMNyMx3LzDwsN+NZqpxkv481Li4On332GT7//HOcPXsW06ZNQ1FRkb6X8MSJE7FgwQL98tOmTUN+fj5mz56Nn3/+GTt27MB7772HGTNmyHUIREREerJfYx03bhyuX7+ON998Ezk5OejSpQu+//57fYemrKwsfc0UAIKCgrBr1y68+OKLePjhh9G8eXPMnj0br7zyilyHQEREpCd7YgWAmTNn1tr0u3fv3mrzIiIicOjQIQtHRUREZDrZm4KJiIgaEiZWIiIiCTGxEhERSYiJlYiISEJMrERERBJiYiUiIpIQEysREZGEmFiJiIgkZBMDRFibEAIAcPv2bY6paSSNRoPi4mIUFBSwzEzAcjMdy8w8LDfT3b59G8DdnCAVh0yslU80CA0NlTkSIiKSW15eHry8vCTbnkMmVh8fHwC6cYilLMyGrPJRe5cvX4anp6fc4dgNlpvpWGbmYbmZ7tatWwgODtbnBKk4ZGKtHNTfy8uLJ6CJPD09WWZmYLmZjmVmHpab6ao+6EWS7Um6NSIiIgfHxEpERCQhh0ysKpUK8fHxUKlUcodiN1hm5mG5mY5lZh6Wm+ksVWYKIXU/YyIiIgfmkDVWIiIiS2FiJSIikhATKxERkYSYWImIiCTkEIn13XffRe/eveHu7g5vb2+j1hFC4M0334RarYabmxuGDBmCX375xbKB2pj8/HxMmDABnp6e8Pb2xnPPPYfCwsI61xk4cCAUCoXBzwsvvGCliOWxYsUKhISEwNXVFb169UJaWlqdyyclJeHBBx+Eq6srOnfujJ07d1opUtthSpklJiZWO6dcXV2tGK38UlJSEB0djcDAQCgUCmzduvW+6+zduxfdunWDSqVCmzZtkJiYaPE4bY2p5bZ3795q55pCoUBOTo5J+3WIxHrnzh2MGTMG06ZNM3qdxYsX46OPPsKqVatw+PBhNG7cGFFRUSgtLbVgpLZlwoQJOH36NJKTk/Htt98iJSUFzz///H3Xmzp1KrKzs/U/ixcvtkK08ti0aRPi4uIQHx+PY8eOISwsDFFRUbh27VqNyx88eBDPPPMMnnvuORw/fhwxMTGIiYnBTz/9ZOXI5WNqmQG60YSqnlOZmZlWjFh+RUVFCAsLw4oVK4xaPiMjAyNGjMCgQYOQnp6OOXPm4E9/+hN27dpl4Uhti6nlVun8+fMG55ufn59pOxYOZN26dcLLy+u+y2m1WhEQECCWLFmin3fz5k2hUqnEhg0bLBih7Thz5owAIP73v//p53333XdCoVCIq1ev1rregAEDxOzZs60QoW0IDw8XM2bM0E9XVFSIwMBAsWjRohqXHzt2rBgxYoTBvF69eok///nPFo3TlphaZsb+3ToKAGLLli11LjNv3jzRsWNHg3njxo0TUVFRFozMthlTbnv27BEAxI0bN+q1L4eosZoqIyMDOTk5GDJkiH6el5cXevXqhdTUVBkjs57U1FR4e3ujR48e+nlDhgyBk5MTDh8+XOe6X331FR544AF06tQJCxYsQHFxsaXDlcWdO3dw9OhRg/PEyckJQ4YMqfU8SU1NNVgeAKKiohzmvDKnzACgsLAQLVu2RFBQEEaPHo3Tp09bI1y75ejnWX116dIFarUakZGROHDggMnrO+Qg/PdT2Z7u7+9vMN/f39/ktnZ7lZOTU635w8XFBT4+PnWWwfjx49GyZUsEBgbi5MmTeOWVV3D+/Hls3rzZ0iFb3e+//46Kiooaz5Nz587VuE5OTo5Dn1fmlFn79u2xdu1aPPzww7h16xY++OAD9O7dG6dPn0aLFi2sEbbdqe08KygoQElJCdzc3GSKzLap1WqsWrUKPXr0QFlZGVavXo2BAwfi8OHD6Natm9HbsdvEOn/+fLz//vt1LnP27Fk8+OCDVorIPhhbbuaqeg22c+fOUKvVePTRR/Hrr7+idevWZm+XHFdERAQiIiL0071790aHDh3wj3/8A2+//baMkVFD0759e7Rv314/3bt3b/z6669YtmwZ/vnPfxq9HbtNrC+99BImTZpU5zKtWrUya9sBAQEAgNzcXKjVav383NxcdOnSxaxt2gpjyy0gIKBaZ5Ly8nLk5+fry8cYvXr1AgBcuHChwSXWBx54AM7OzsjNzTWYn5ubW2sZBQQEmLR8Q2NOmd1LqVSia9euuHDhgiVCbBBqO888PT1ZWzVReHg49u/fb9I6dptYmzVrhmbNmllk26GhoQgICMDu3bv1ibSgoACHDx82qWexLTK23CIiInDz5k0cPXoU3bt3BwD88MMP0Gq1+mRpjPT0dAAw+ILSUDRq1Ajdu3fH7t27ERMTAwDQarXYvXs3Zs6cWeM6ERER2L17N+bMmaOfl5ycbFAja8jMKbN7VVRU4NSpUxg+fLgFI7VvERER1W7jcqTzTErp6emm//+qV9cnO5GZmSmOHz8uEhISRJMmTcTx48fF8ePHxe3bt/XLtG/fXmzevFk//de//lV4e3uLbdu2iZMnT4rRo0eL0NBQUVJSIschyGLo0KGia9eu4vDhw2L//v2ibdu24plnntG/f+XKFdG+fXtx+PBhIYQQFy5cEG+99ZY4cuSIyMjIENu2bROtWrUS/fv3l+sQLG7jxo1CpVKJxMREcebMGfH8888Lb29vkZOTI4QQ4tlnnxXz58/XL3/gwAHh4uIiPvjgA3H27FkRHx8vlEqlOHXqlFyHYHWmlllCQoLYtWuX+PXXX8XRo0fF008/LVxdXcXp06flOgSru337tv7/FgCxdOlScfz4cZGZmSmEEGL+/Pni2Wef1S9/8eJF4e7uLubOnSvOnj0rVqxYIZydncX3338v1yHIwtRyW7Zsmdi6dav45ZdfxKlTp8Ts2bOFk5OT+O9//2vSfh0iscbGxgoA1X727NmjXwaAWLdunX5aq9WKN954Q/j7+wuVSiUeffRRcf78eesHL6O8vDzxzDPPiCZNmghPT08xefJkgy8jGRkZBuWYlZUl+vfvL3x8fIRKpRJt2rQRc+fOFbdu3ZLpCKzj448/FsHBwaJRo0YiPDxcHDp0SP/egAEDRGxsrMHyX3/9tWjXrp1o1KiR6Nixo9ixY4eVI5afKWU2Z84c/bL+/v5i+PDh4tixYzJELZ/K20Du/aksp9jYWDFgwIBq63Tp0kU0atRItGrVyuD/m6Mwtdzef/990bp1a+Hq6ip8fHzEwIEDxQ8//GDyfvnYOCIiIgnxPlYiIiIJMbESERFJiImViIhIQkysREREEmJiJSIikhATKxERkYSYWImIiCTExEpERCQhJlYiIiIJMbES2Yi8vDz4+fnh0qVLcodSL08//TQ+/PBDucMgkg0TK5EVLFq0CD179oSHhwf8/PwQExOD8+fPGyzz7rvvYvTo0QgJCTGYf/LkSTzxxBPw9fWFq6srOnbsiCVLlqC8vFzyOKXY1+uvv453330Xt27dqvbe5MmT8frrr0sZMpHNYWIlsoJ9+/ZhxowZOHToEJKTk6HRaPDYY4+hqKgIAFBcXIw1a9bgueeeM1gvJSUFjzzyCNzc3LBt2zacOHECr7zyCpYuXYonnngCWq1Wshil2lenTp3QunVrfPnllwbzKyoq8O2332LUqFGSxUxkkyR5hAARmeTatWsCgNi3b58QQoikpCTRrFkzg2XKy8tFq1atxIQJE6qtf/bsWaFUKsXq1asliUfqfSUkJIi+ffsazEtJSRFqtVpotVohhBA3btwQAMT+/fuFEEL88ssvon379uK1117TL0Nkj1hjJZJBZTOpj48PAODHH3/UP1C+UlpaGi5evIi5c+dWW//BBx/EiBEjsGnTJknikXpf4eHhSEtLQ1lZmX7e9u3bER0dDYVCAUDX7KxQKBAWFob9+/djwIABePXVV/HOO+/olyGyR0ysRFam1WoxZ84c9OnTB506dQIAZGZmIjAw0GC5jIwMAEDbtm1r3E7btm2RmZkpSUxS7yswMBB37txBTk6Oft62bdsMmoFPnDiB1q1bY9u2bRg3bhw2bNiAiRMnmnkERLaDiZXIymbMmIGffvoJGzdu1M8rKSmBq6urwXKenp4AgPz8/Bq3c+PGDf0yVc2fPx8KhaLOn3Pnzkmyr9q4ubkB0F07BoCzZ8/it99+w6OPPqpf5sSJE8jJycGkSZMQEBCA/v37G719IlvGxEpkRTNnzsS3336LPXv2oEWLFvr5DzzwAG7cuGGwbEREBJRKJf79739X205FRQV27dqFvn37VnvvpZdewtmzZ+v8adWqVb329fjjj+Ppp59Gz5490bp1axw5csRgncoE3axZMwC6ZuDIyEiDLw8nTpxAjx49sHfvXhw/fhxbt26trdiI7IvcF3mJHIFWqxUzZswQgYGB4ueff672/pIlS0RYWFi1+S+99JIIDAwUV69erbZ848aNRWZmpmQxmrKvVq1aiffff18IIcSXX34pnn76aYN1Vq9eLVq0aKGfjoiIEOvWrdNPl5eXC1dXV7Ft2zYhhBBjxowRYWFh7LREDQITK5EVTJs2TXh5eYm9e/eK7Oxs/U9xcbEQQoiTJ08KFxcXkZ+fr1/n9u3bIjMzU/Tu3Vu0a9dOHD16VAghxOLFi4VSqRRr164V2dnZory8vN7xmbKv27dvi8DAQP1+jx49KqKiogy2FxsbK6ZMmSKEECI3N1colUpx/fp1/ftnzpwRAPTJ+uTJk0KhUIikpKR6HwuR3JhYiawAQI0/VWtx4eHhYtWqVfrp+Ph4g2VjY2Nr3FZGRka94zNlXwcPHhSRkZH6dT/77DPx8ssv66dLSkqEl5eXSE1NFULoaq99+vQx2N+GDRuEt7e3wbwnn3xSdOzYUVRUVNT7eIjkpBBCCKu2PRNRjXbs2IG5c+fip59+gpOT7XZ/+Mc//oGlS5fip59+QkFBAQYPHozNmzejdevWAICVK1diy5Yt+M9//gMAGDVqFPr27Yt58+bJGTaR1bjIHQAR6YwYMQK//PILrl69iqCgILnDqdXJkycxfPhwdO/eHUIILF68WJ9UAUCpVOLjjz/WT/ft2xfPPPOMHKESyYI1ViIySb9+/bB+/XqbTv5EcrLd9iYiskm2XqMmkhtrrERERBJijZWIiEhCTKxEREQSYmIlIiKSEBMrERGRhJhYiYiIJMTESkREJCEmViIiIgkxsRIREUmIiZWIiEhCTKxEREQSYmIlIiKSEBMrERGRhP4fuZxrxgvM9voAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(5,5))\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "ax.grid(zorder=-1)\n",
    "\n",
    "test_bar = (otype_ILR > 2)\n",
    "test_notbarnotdisk = (otype_ILR <= 2) & (otype_uha == 1) & (otype_vILR == 2)\n",
    "test_disk = ~test_bar & ~test_notbarnotdisk\n",
    "\n",
    "\n",
    "ax.scatter(\n",
    "    2*Omega_sec[test_bar]/kappa_sec[test_bar],\n",
    "    omegaz_sec[test_bar]/kappa_sec[test_bar],\n",
    "    s=3, label='Bar'\n",
    ")\n",
    "\n",
    "ax.scatter(\n",
    "    2*Omega_sec[test_disk]/kappa_sec[test_disk],\n",
    "    omegaz_sec[test_disk]/kappa_sec[test_disk],\n",
    "    s=3, label='Disk'\n",
    ")\n",
    "\n",
    "ax.scatter(\n",
    "    2*Omega_sec[test_notbarnotdisk]/kappa_sec[test_notbarnotdisk],\n",
    "    omegaz_sec[test_notbarnotdisk]/kappa_sec[test_notbarnotdisk],\n",
    "    s=3, label='not Bar, not Disk'\n",
    ")\n",
    "\n",
    "\n",
    "ax.set_xlabel('$2(\\\\Omega - \\\\Omega_p)/\\\\kappa$')\n",
    "ax.set_ylabel('$\\\\omega_z\\\\kappa$')\n",
    "ax.legend(\n",
    "    loc='upper center',\n",
    "    bbox_to_anchor=(0.5, 1.15),\n",
    "    ncols=3\n",
    ")\n",
    "\n",
    "ax.set_xlim(-1.0,1.5)\n",
    "ax.set_ylim(0.5, 2)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}