DP03_03_Trans-Neptunian_Objects.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6b7f500d-ba87-4a25-89ef-5c3038304aff",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-11T17:31:37.391082Z",
     "iopub.status.busy": "2023-07-11T17:31:37.390105Z",
     "iopub.status.idle": "2023-07-11T17:31:37.396377Z",
     "shell.execute_reply": "2023-07-11T17:31:37.395279Z",
     "shell.execute_reply.started": "2023-07-11T17:31:37.391043Z"
    },
    "tags": []
   },
   "source": [
    "# Trans-Neptunian Objects (TNOs)\n",
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "Contact author(s): Andrés A. Plazas Malagón<br>\n",
    "Last verified to run: 2024-12-17 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_50 <br>\n",
    "Container Size: medium <br>\n",
    "Targeted learning level: intermediate<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7de485fc-1275-4fef-9fd5-ad1374609cc0",
   "metadata": {},
   "source": [
    "**Description:** Explore the trans-Neptunian object populations in DP0.3. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b311418-01c0-49ef-9263-2321448ccf8b",
   "metadata": {},
   "source": [
    "**Skills:** Use of the DP0.3 catalogs to study TNO populations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "121ce17e-38eb-49ca-b95e-2a32d7b14250",
   "metadata": {},
   "source": [
    "**LSST Data Products:**  DP0.3 catalogs `SSObject`, `DiaSource`, and `MPCORB` (10-year catalogs)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c9393c0-0c43-4178-bd58-33af12540613",
   "metadata": {},
   "source": [
    "**Packages:** `lsst.rsp`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f62aabe-9d91-4541-ad95-5a0df25f2688",
   "metadata": {},
   "source": [
    "**Credits:** Developed by Andrés A. Plazas Malagón in collaboration with Melissa Graham, Ryan Lau, Pedro Bernardinelli (University of Washington), and the Rubin Community Science Team for DP0.3.\n",
    "This tutorial is primarily based on a Jupyter Notebook by Pedro Bernardinelli, and a Jupyter Notebook by Meg Schwamb. \n",
    "This notebook has made use of suggestions in the [Accessible Authoring Checklist](https://iota-school.github.io/accessibility_hackathon/hack#checklist), and has made use of NASA's Astrophysics Data System Bibliographic Services."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a231d35-fd48-4caa-8519-c605f05e9330",
   "metadata": {},
   "source": [
    "**Get Support:**\n",
    "Find DP0-related documentation and resources at <a href=\"https://dp0.lsst.io\">dp0.lsst.io</a>. \n",
    "Questions are welcome as new topics in the <a href=\"https://community.lsst.org/c/support/dp0\">Support - Data Preview 0 Category</a> of the Rubin Community Forum. \n",
    "Rubin staff will respond to all questions posted there."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7f445fe-3f60-40b6-b2c2-704029b673ea",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "### 1.1 What are Trans-Neptunian objects (TNOs)\n",
    "\n",
    "A trans-Neptunian object (TNO) is a minor planet within the Solar System which orbits the Sun at an average distance greater than that of Neptune (> 30.1 astronomical units; au). \n",
    "The population of TNOs is scientifically valuable because they are largely intact relics from the era of planetary formation. \n",
    "\n",
    "In the current model explaining the formation of the outer Solar System, the small building blocks of planets, called planetesimals, that eventually became TNOs can be categorized into two broad groups.\n",
    "\n",
    "1. **Outwardly Migrated \"Hot\" TNOs.**\n",
    "Some planetesimals formed within the inner region of the protoplanetary disk, located at a distance less than 30 au from the Sun. \n",
    "These planetesimals were later pushed outward by interactions with a migrating Neptune. \n",
    "As a result, they acquired dynamically \"hot\" orbits, meaning their inclination distribution\n",
    "includes high inclination objects (beyond 5 degrees), with a wide range of eccentricities. \n",
    "\n",
    "2. **In-Situ \"Cold\" TNOs.**\n",
    "Other planetesimals formed directly in the trans-Neptunian region, \n",
    "where they remain with dynamically \"cold\" orbits (inclinations < 5 degrees).\n",
    "\n",
    "\n",
    "#### 1.1.1 Dynamical classes\n",
    "\n",
    "In Section 3, clustering within plots of eccentricity or inclination versus semi-major axis\n",
    "will be used to illustrate the following TNO dynamical classes.\n",
    "\n",
    "**Resonant TNOs** are currently found in stable mean motion resonances with Neptune, meaning they have a specific orbital relationship with the planet. \n",
    "It is presumed that they were transported outward while maintaining their resonant positions as Neptune migrated during the early history of the Solar System.\n",
    "\n",
    "**Scattering TNOs** can have orbits extending to large distances, but their orbits are currently experiencing significant variations due to perturbations caused by Neptune when their closest distance to the Sun (perihelia $q$) is less than approximately 40 au, closest to Neptune. \n",
    "Scattering TNOs are believed to be remnants of a population that was scattered during the early formation of the Solar System and have managed to survive for a long time, and whose orbits are experiencing evolution, primarly due to Neptune. \n",
    "\n",
    "**Detached TNOs**, in contrast, have higher perihelia and are dynamically independent from Neptune, as they currently maintain stable orbits that are not significantly influenced by Neptune's gravitational effects.\n",
    "Detached TNOs might have originated as Scattering TNOs whose perihelia were raised through past resonance interactions with Neptune or through interactions with a distant, massive planet in the outer Solar System (though this hypothetical planet has not been observed yet). \n",
    "This class of TNOs is among the most enigmatic one, and the questions about where they formed and how they got to their orbits are still open.\n",
    "\n",
    "**Extreme TNOs** are characterized by their exceptional distances from the Sun, with semi-major axes ($a$) exceeding 150 au (note that definitions may vary slightly) and $q > 40$ au.\n",
    "The formation mechanism responsible for producing these distant extreme TNOs remains unclear.\n",
    "One prominent hypothesis involves the presence of a massive object often referred to as \"Planet X\" or \"Planet 9\" situated at approximately 400 - 800 au from the Sun.\n",
    "The gravitational influence of this hypothetical Neptune-mass planet could have played a significant role in sculpting the distribution of these extreme TNOs, leading to the anisotropies observed in their orbital characteristics. Note that the \"extreme\" class is composed of a mix of \"scattering\" and \"detached\" objects, and is only meanigful in the presence of \"Planet 9/X\" - if there’s no \"Planet 9/X\", the \"extreme\" class is not required; see [Batygin et al. 2019](https://ui.adsabs.harvard.edu/abs/2019PhR...805....1B/abstract).\n",
    "\n",
    "**Classical TNOs** correspond to objects that are not in any of the previous classes, and are situated between Neptune's 3:2 and 2:1 mean motion resonances, at distances ranging from 39 to 47 au, and characterized by low eccentricity ($e \\lesssim 0.24$) and inclination ($i \\lesssim$ 5 degrees). \n",
    "This group is a mix of two types of objects: some are native and remained dynamically \"cold,\" while others were originally \"hot\" objects that migrated into this region.\n",
    "\n",
    "\n",
    "#### 1.1.2 Comets\n",
    "In addition to the TNO dynamical classes, the DP0.3 simulation also includes other types of solar System objects such as long-period comets, and they will also be visualized in section 4. **Comets** are icy objects in space with irregular shapes that orbit the sun in highly-elliptical paths. They come in two types: short-period comets (less than 200 years to orbit) found in the Kuiper Belt beyond Neptune (30 au - 50 au), and long-period comets (over 200 years) that have more random paths and are thought to come from the Oort cloud, far away from the Sun. The Oort cloud has not been directly observed, but scientists believe it contains many icy objects.\n",
    "\n",
    "#### 1.1.3 References\n",
    "\n",
    "The above information was compiled from the following refrences:\n",
    "\n",
    "*A Search of the Full Six Years of the Dark Energy Survey for Outer Solar System Objects*, [Bernardinelli et al. (2022), ApJS](https://ui.adsabs.harvard.edu/abs/2022ApJS..258...41B/abstract)\n",
    "\n",
    "*Transneptunian Space*, [Gladman and Volk (2021), AR&AA](https://ui.adsabs.harvard.edu/abs/2021ARA%26A..59..203G/abstract)\n",
    "\n",
    "*The Outer Solar System Origins Survey (OSSOS). VII. 800+ Trans-Neptunian Objects — The Complete Data Release*, [Banister et al (2018), ApJS](https://ui.adsabs.harvard.edu/abs/2018ApJS..236...18B/abstract)\n",
    "\n",
    "*Modern celestial mechanics : aspects of solar system dynamics*, [Morbidelli et al. 2002](https://ui.adsabs.harvard.edu/abs/2002mcma.book.....M/abstract)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ee13448-23c9-4c15-9091-dbdb6b52f3f3",
   "metadata": {},
   "source": [
    "### 1.2 Import packages\n",
    "\n",
    "The [matplotlib](https://matplotlib.org/) (and especially sublibrary `matplotlib.pyplot`), [numpy](http://www.numpy.org/), and [astropy](http://www.astropy.org/) libraries are widely used Python libraries for plotting, numerical analysis, and astronomical data analysis. The [seaborn](https://seaborn.pydata.org/index.html) library is a statistical data visualization library based on `matplotlib`. The [Garbage Collector module](https://docs.python.org/3/library/gc.html), `gc`, assists with memory management.\n",
    "\n",
    "The `lsst.rsp` package provides access to the Table Access Protocol (TAP) service for access to the DP0 catalogs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c45fdcfe-be27-46c0-a24e-8173356440b0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "from matplotlib.patches import Polygon\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import gc\n",
    "from lsst.rsp import get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c40b37d5-c522-4d14-8e13-d2cc0b8ad4bb",
   "metadata": {},
   "source": [
    "### 1.3 Define functions and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fa5094e-92a2-4de2-8e10-69e643621cd0",
   "metadata": {},
   "source": [
    "The following functions and parameters are defined once, and used throughout this notebook.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41f068b3-67a2-4d3d-9236-78bf2f7967ed",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-12-12T18:24:07.822868Z",
     "iopub.status.busy": "2023-12-12T18:24:07.821830Z",
     "iopub.status.idle": "2023-12-12T18:24:07.828412Z",
     "shell.execute_reply": "2023-12-12T18:24:07.827220Z",
     "shell.execute_reply.started": "2023-12-12T18:24:07.822834Z"
    }
   },
   "source": [
    "**_Function_**: `remove_figure`:\n",
    "\n",
    "Helper function to clean up plots after creating them to manage memory allocation. See Notebook `DP02_03a`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e4c2e11-9570-45d2-a532-cf92dff66273",
   "metadata": {},
   "outputs": [],
   "source": [
    "def remove_figure(fig):\n",
    "    \"\"\"\n",
    "    Remove a figure to reduce memory footprint.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    fig: matplotlib.figure.Figure\n",
    "        Figure to be removed.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    None\n",
    "    \"\"\"\n",
    "    # get the axes and clear their images\n",
    "    for ax in fig.get_axes():\n",
    "        for im in ax.get_images():\n",
    "            im.remove()\n",
    "    fig.clf()       # clear the figure\n",
    "    plt.close(fig)  # close the figure\n",
    "    gc.collect()    # call the garbage collector"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d0270b0-5f41-45f3-8534-2a718478a725",
   "metadata": {},
   "source": [
    "**_Parameter_**: `a_Neptune`\n",
    "\n",
    "Neptune's mean semimajor axis ($a_{Neptune}$) in Astronomical Units (au)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bf0ecb1-d477-4df7-aa9e-fdb68327b597",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "a_Neptune = 30.1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90bc8f87-06a8-4080-87ab-ca687daa38a2",
   "metadata": {},
   "source": [
    "**_Function_**: `estimateDiameter`\n",
    "\n",
    "Estimate the diameter ($d$), in kilometers, of a Solar System object.\n",
    "The inputs are the object's absolute H magnitude ($H$) and albedo ($A$).\n",
    "The term \"albedo\" refers to the fraction of light reflected, from 0 to 1.\n",
    "\n",
    "$d = 10^{(3.1236\\ -\\ 0.5 \\log(A)\\ -\\ 0.2H)}$\n",
    "\n",
    "This function is used in Section 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c189ce6-ccd1-47e0-8007-c53bdc711167",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def estimateDiameter(H, albedo=0.1):\n",
    "    \"\"\"\n",
    "    Return Solar System object size in kilometers.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    H : float\n",
    "        Solar System absolute magnitude\n",
    "    albedo: float, optional.\n",
    "        Solar system object albedo. Default is 0.1.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    d : float\n",
    "        Estimated diameter of the object, in kilometers.\n",
    "    \"\"\"\n",
    "\n",
    "    d = np.power(10., 3.1236 - 0.5*np.log10(albedo) - 0.2*H)\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a47f0719-8189-4093-a321-a40f116e8111",
   "metadata": {},
   "source": [
    "**_Function_**: `qVector`\n",
    "\n",
    "Compute the three-component coordinate position vector ${\\bf{q}} = (q1, q2, q3)$ of an object on the plane of an ellipse where: <br>\n",
    "$q1 = a (\\cos(E) - e)$ <br>\n",
    "$q2 = a \\sqrt{(1 - e^2)} \\sin(E)$, and <br>\n",
    "$q3 = 0$.\n",
    "\n",
    "The input values are the eccentric anomaly ($E$, describing the position of the object along the orbit), semi-major axis ($a$), and eccentricity ($e$). In this orthogonal reference frame, the origin lies on the focus of the ellipse (where the main body is located), and $q_1$ is oriented towards the pericenter of the orbit (see Chapter 1 of [Morbidelli et al. 2002](https://ui.adsabs.harvard.edu/abs/2002mcma.book.....M/abstract)).\n",
    "\n",
    "This function is used in Section 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76f94b85-04f2-4899-a7c2-5d4e423669d8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def qVector(E, a, e):\n",
    "    \"\"\"\n",
    "    Computes the coordinate vector q =\n",
    "    (a * (cos(E) - e),\n",
    "    a * sqrt(1 - e^2) * sin(E), 0),\n",
    "    on the plane of the ellipse.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    E : `float` or `np.array`\n",
    "        Eccentric anomaly.\n",
    "    a : `float` or `np.array`\n",
    "        Semi-major axis.\n",
    "    e : `float` or `np.array`\n",
    "        Eccentricity.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    coordinates : `numpy.array`\n",
    "        Coordinate vector(s) in the plane\n",
    "        of the ellipse.\n",
    "    \"\"\"\n",
    "    q1 = a * (np.cos(E) - e)\n",
    "    q2 = a * np.sqrt(1. - e ** 2) * np.sin(E)\n",
    "\n",
    "    if isinstance(q1, float):\n",
    "        coordinates = np.array([q1, q2, 0.])\n",
    "    else:\n",
    "        coordinates = np.array([q1, q2, np.zeros_like(q1)])\n",
    "    return coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a3891c1-64cb-4cd1-90a7-89ca97db226b",
   "metadata": {},
   "source": [
    "**_Function_**: `rotationMatrix`\n",
    "\n",
    "Compute the rotation matrix between the position vector $\\bf{q}$ in the plane of the ellipse to a position vector $\\bf{r}$ with respect to an arbitrary orthogonal reference frame $(x, y, z)$: \n",
    "\n",
    "${\\bf{r}} = R_{\\rm {xq}} (\\Omega, \\omega, \\textit{i})\\ {\\bf{q}}$\n",
    "\n",
    "The rotation matrix is a function of the longitude of the node, $\\Omega$, the argument of the perihelion $\\omega$, and the orbit incination,  $\\textit{i}$. \n",
    "\n",
    "Note that this equation provides a one-to-one correspondence between the six orbital elements of an object and the ${\\bf{r}}$ position vector.\n",
    "\n",
    "(see Chapter 1 of [Morbidelli et al. 2002](https://ui.adsabs.harvard.edu/abs/2002mcma.book.....M/abstract)).\n",
    "\n",
    "This function is used in Section 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34ea2c86-562d-492b-8dfb-77dc9ece6b43",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rotationMatrix(peri, node, incl):\n",
    "    \"\"\"\n",
    "    Computes the rotation matrix R = R(peri, node, incl) that\n",
    "    maps from the plane of the ellipse to 3D space aligned\n",
    "    with the ecliptic plane.\n",
    "\n",
    "    node : `float`\n",
    "        Longitude of the node (rad).\n",
    "    peri : `float`\n",
    "        Argument of pericenter (rad).\n",
    "    incl : `float`\n",
    "        Orbit inclination (rad).\n",
    "    \"\"\"\n",
    "    cO = np.cos(np.pi * node / 180)\n",
    "    sO = np.sin(np.pi * node / 180)\n",
    "    co = np.cos(np.pi * peri / 180)\n",
    "    so = np.sin(np.pi * peri / 180)\n",
    "    ci = np.cos(np.pi * incl / 180)\n",
    "    si = np.sin(np.pi * incl / 180)\n",
    "\n",
    "    R = np.array([[cO * co - sO * so * ci, -cO * so - sO * co * ci, sO * si],\n",
    "                  [sO * co + cO * so * ci, -sO * so + cO * co * ci, -cO * si],\n",
    "                  [si * so, si * co, ci]])\n",
    "\n",
    "    return R"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c31845a4-a67b-4257-8d9f-6a38afa6a6f3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:23:23.305455Z",
     "iopub.status.busy": "2023-07-26T05:23:23.304489Z",
     "iopub.status.idle": "2023-07-26T05:23:24.085144Z",
     "shell.execute_reply": "2023-07-26T05:23:24.083942Z",
     "shell.execute_reply.started": "2023-07-26T05:23:23.305420Z"
    },
    "tags": []
   },
   "source": [
    "**_Function_**: `fullEllipse`\n",
    "\n",
    "Compute the coordinates of the each point of the ellipse (\"full\" ellipse) in an arbitrary coordinate frame ${\\bf{r}}=(x,y,z)$, strarting from the coordinates in the plane from the ellipse, ${\\bf{q}}=(q_1, q_2, q_3)$, using the rotation matrix defined above:\n",
    "\n",
    "${\\bf{r}} = R_{\\rm {xq}} (\\Omega, \\omega, \\textit{i})\\ {\\bf{q}}$\n",
    "\n",
    "This function is used in Section 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64f01d19-7d23-481d-8868-6c51d30aa9e3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def fullEllipse(a, e, incl, node, peri):\n",
    "    \"\"\" Computes the coordinates of the\n",
    "    full ellipse in 3D space.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    a : `float`\n",
    "        Semi-major axis (au).\n",
    "    e : `float`\n",
    "        Eccentricity.\n",
    "    incl : `float`\n",
    "        Orbit inclination (rad).\n",
    "    node : `float`\n",
    "        Longitud eof the node (rad).\n",
    "    peri : `float`\n",
    "        Argument of the preihelion (rad).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    r : `numpy.array`\n",
    "        (x, y, z) ellipse coordinates\n",
    "        in 3D space.\n",
    "\n",
    "    Notes:\n",
    "    `E` is the eccentric anomaly.\n",
    "    \"\"\"\n",
    "    E = np.linspace(0, 2 * np.pi, 2000)\n",
    "    q_vec = qVector(E, a, e)\n",
    "    Rot = rotationMatrix(peri, node, incl)\n",
    "    r = np.einsum('ij, j...', Rot, q_vec)\n",
    "\n",
    "    return r"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aa4dff9-38a1-4f1d-860e-f9e37a57531c",
   "metadata": {},
   "source": [
    "**_Function_**: `gradientFill`\n",
    "\n",
    "Plot a line with a linear gradient filled beneath it.\n",
    "\n",
    "This function is used in Section 3.2 when plotting eccentricy (or inclination) versus semi-major axis, and\n",
    "was taken from [Stack Overflow](https://stackoverflow.com/questions/29321835/is-it-possible-to-get-color-gradients-under-curve-in-matplotlib)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "455b69fa-5e11-4cb7-9c67-eced6a9381cb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def gradientFill(x, y, alpha=1, fill_color=None, ax=None,\n",
    "                 invert=False, **kwargs):\n",
    "    \"\"\"\n",
    "    Plot a line with a linear alpha gradient filled beneath it.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x, y : array-like\n",
    "        The data values of the line.\n",
    "    alpha : float\n",
    "        Transparency for plotted element, 0 for transparent, 1 for opaque.\n",
    "    fill_color : a matplotlib color specifier (string, tuple) or None\n",
    "        The color for the fill. If None, the color of the line will be used.\n",
    "    ax : a matplotlib Axes instance\n",
    "        The axes to plot on. If None, the current pyplot axes will be used.\n",
    "    invert : boolean\n",
    "        Whether to invert the gradient.\n",
    "\n",
    "    Additional arguments are passed on to matplotlib's `plot` function.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    line : a Line2D instance\n",
    "        The line plotted.\n",
    "    im : an AxesImage instance\n",
    "        The transparent gradient clipped to just the area beneath the curve.\n",
    "    \"\"\"\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "\n",
    "    line, = ax.plot(x, y, alpha=0, **kwargs)\n",
    "    if fill_color is None:\n",
    "        fill_color = line.get_color()\n",
    "\n",
    "    zorder = line.get_zorder()\n",
    "\n",
    "    z = np.empty((100, 1, 4), dtype=float)\n",
    "    rgb = mcolors.colorConverter.to_rgb(fill_color)\n",
    "    z[:, :, :3] = rgb\n",
    "    z[:, :, -1] = np.linspace(0, alpha, 100)[:, None]\n",
    "\n",
    "    xmin, xmax, ymin, ymax = x.min(), x.max(), y.min(), y.max()\n",
    "    im = ax.imshow(z, aspect='auto', extent=[xmin, xmax, ymin, ymax],\n",
    "                   origin='upper', zorder=zorder)\n",
    "\n",
    "    xy = np.column_stack([x, y])\n",
    "\n",
    "    if invert == 0:\n",
    "        xy = np.vstack([[xmin, ymin], xy, [xmax, ymin], [xmin, ymin]])\n",
    "    else:\n",
    "        xy = np.vstack([[xmax, ymax], xy, [xmin, ymax], [xmin, ymin]])\n",
    "\n",
    "    clip_path = Polygon(xy, facecolor='none', edgecolor='none', closed=True)\n",
    "    ax.add_patch(clip_path)\n",
    "    im.set_clip_path(clip_path)\n",
    "\n",
    "    return line, im"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6df7032a-0860-4f64-af46-f9db6b48c10b",
   "metadata": {},
   "source": [
    "**_Function_**: `plotResonances`\n",
    "\n",
    "Calculate mean motion resonances with Neptune and plot them as vertical dashed lines, using $P^2 \\propto a^3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79632108-f376-4a5d-83fa-3a52ff22029b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plotResonances(ax, k, k_N, y_loc=None, a_shift=0.75, alpha=0.2):\n",
    "    \"\"\"Plot Neptune resonances.\n",
    "\n",
    "    Implements Kepler's Third law\n",
    "    `a = a_N*(k/k_N)**(2/3)`, where `a`\n",
    "    is the semi-major axis of a TNO in\n",
    "    orbital resonance with Neptune in au,\n",
    "    `a_N` is Neptune's mean semi-major\n",
    "    axis (`a_N = 30.1` au) and (k/K_N)\n",
    "    is the orbital period ratio.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    ax : `matplotlib.axes`\n",
    "        Figure axes.\n",
    "    k : `int`\n",
    "        Numerator of the period ratio\n",
    "        with respect to Neptune's\n",
    "        orbital period.\n",
    "    k_N : `int`\n",
    "        Denominator in period ratio,\n",
    "        corresponding to Neptune.\n",
    "    y_loc : `float`\n",
    "        Location of the vertical line\n",
    "        showing the resonance.\n",
    "    alpha : `float`\n",
    "       Line transparency.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    ax : `matplotlib.axes`\n",
    "       Figure axes.\n",
    "    \"\"\"\n",
    "    a_res = a_Neptune * np.cbrt(float(k)**2/float(k_N)**2)\n",
    "    ax.axvline(a_res, linewidth=2, alpha=alpha, linestyle='--')\n",
    "    if y_loc is not None:\n",
    "        ax.annotate(r'${}:{}$'.format(k, k_N),\n",
    "                    (a_res + a_shift, y_loc),\n",
    "                    rotation=90, fontsize=12)\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3f85db0-7688-4928-b011-8587fffa59e3",
   "metadata": {},
   "source": [
    "**_Function_**: `plotHistogram`\n",
    "\n",
    "Plot a histogram of the specified data for a given class label.\n",
    "\n",
    "This function is used in Section 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26a80311-bf85-4c00-926d-1b16ea4fafb9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plotHistogram(data, label, nbins, color, linestyle,\n",
    "                  histtype='step',\n",
    "                  density=False, cumulative=False,\n",
    "                  log=False):\n",
    "    \"\"\"Plot a histogram of the specified data for a given class label.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    data : numpy.ndarray\n",
    "        The data to be plotted in the histogram.\n",
    "    label : str\n",
    "        The class label for which the histogram is plotted.\n",
    "    nbins : int\n",
    "        The number of bins to use for the histogram.\n",
    "    color : str\n",
    "        The color of the histogram plot.\n",
    "    linestyle : str\n",
    "        The line style of the histogram plot.\n",
    "    histtype : str, optional\n",
    "        The type of histogram bars. Default is 'step'.\n",
    "    density : bool, optional\n",
    "        If True, the histogram is normalized to form a probability density.\n",
    "        Default is False.\n",
    "    cumulative : bool, optional\n",
    "        If True, a cumulative histogram is plotted. Default is False.\n",
    "    log : bool, optional\n",
    "        use logarithmic scale. Default is False.\n",
    "    \"\"\"\n",
    "    tx = np.where(tnos['class'] == label)[0]\n",
    "    plt.hist(data[tx], bins=nbins, histtype=histtype,\n",
    "             density=density, cumulative=cumulative,\n",
    "             color=color, label=label, linestyle=linestyle,\n",
    "             log=log)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a9f7ec6-8cb9-4c4e-b274-42cff493320a",
   "metadata": {},
   "source": [
    "**_Function_**: `createSubplot`\n",
    "\n",
    "Create a subplot within a grid layout.\n",
    "\n",
    "This function is used in Section 3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc5f8992-1341-480a-a126-056437fd27d2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def createSubplot(rows, cols, plot_num, title, x_label, y_label):\n",
    "    \"\"\"\n",
    "    Create a subplot within a grid layout.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    rows : int\n",
    "        Number of rows in the grid.\n",
    "    cols : int\n",
    "        Number of columns in the grid.\n",
    "    plot_num : int\n",
    "        Position of the subplot in the grid (starting from 1).\n",
    "    title : str\n",
    "        Title of the subplot.\n",
    "    x_label : str\n",
    "        Label for the x-axis of the subplot.\n",
    "    y_label : str\n",
    "        Label for the y-axis of the subplot.\n",
    "    \"\"\"\n",
    "    plt.subplot(rows, cols, plot_num)\n",
    "    plt.title(title)\n",
    "    plt.xlabel(x_label)\n",
    "    plt.ylabel(y_label)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "506c6836-efe7-4f6a-b5a1-81e0f65a0bf0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-25T00:26:00.464887Z",
     "iopub.status.busy": "2023-07-25T00:26:00.464477Z",
     "iopub.status.idle": "2023-07-25T00:26:01.116431Z",
     "shell.execute_reply": "2023-07-25T00:26:01.115215Z",
     "shell.execute_reply.started": "2023-07-25T00:26:00.464861Z"
    },
    "tags": []
   },
   "source": [
    "**_Parameters_**: Set plotting parameters to make visually accessible plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20b76f7e-8bee-482f-8a59-7e0ddfde813e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "params = {'lines.linewidth': 2,\n",
    "          'lines.linestyle': '-',\n",
    "          'lines.color': 'black',\n",
    "          'patch.linewidth': 2,\n",
    "          'font.family': 'serif',\n",
    "          'font.weight': 'normal',\n",
    "          'font.size': 11.0,\n",
    "          'text.color': 'black',\n",
    "          'axes.edgecolor': 'black',\n",
    "          'axes.linewidth': 2.0,\n",
    "          'axes.grid': False,\n",
    "          'axes.titlesize': 'large',\n",
    "          'axes.labelsize': 'large',\n",
    "          'axes.labelweight': 'bold',\n",
    "          'axes.labelcolor': 'black',\n",
    "          'axes.formatter.limits': (-4, 4),\n",
    "          'xtick.major.size': 10,\n",
    "          'xtick.minor.size': 5,\n",
    "          'xtick.major.pad': 8,\n",
    "          'xtick.minor.pad': 8,\n",
    "          'xtick.labelsize': 'large',\n",
    "          'xtick.minor.width': 2.0,\n",
    "          'xtick.major.width': 2.0,\n",
    "          'xtick.minor.visible': True,\n",
    "          'ytick.major.size': 10,\n",
    "          'ytick.minor.size': 5,\n",
    "          'ytick.major.pad': 8,\n",
    "          'ytick.minor.pad': 8,\n",
    "          'ytick.labelsize': 'large',\n",
    "          'ytick.minor.width': 1.0,\n",
    "          'ytick.major.width': 1.0,\n",
    "          'ytick.minor.visible': True,\n",
    "          'legend.numpoints': 1,\n",
    "          'legend.fontsize': 'large',\n",
    "          'legend.shadow': False,\n",
    "          'legend.frameon': True,\n",
    "          'legend.fancybox': True,\n",
    "          'figure.titlesize': 'large',\n",
    "          'figure.figsize': (8, 7),\n",
    "          'figure.facecolor': 'white'}\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15715491-5943-4d01-a148-75c1ee42c57d",
   "metadata": {},
   "source": [
    "## 2. Query the `MPCORB` year 10 catalog for all TNOs\n",
    "\n",
    "Start the TAP service for DP0.3 Solar System data, ``ssotap``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "061db9be-1333-4a21-a2be-4c6425d0512b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"ssotap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8130937d-480a-4cd7-8d08-832c11916bbf",
   "metadata": {},
   "source": [
    "The contents of the `MPCORB` table are explored and explained in the introductory DP0.3 notebook. In this notebook, we willl use the year 10 catalog, `dp03_catalogs_10yr.MPCORB`. A year 1 catalog is also availabe (`dp03_catalogs_1yr.MPCORB`) \n",
    "\n",
    "**Option:** Uncomment the following cell to view the column descriptions for the `MPCORB` catalog as a `pandas` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07f8002b-d500-4c74-9e71-74dca7d7ef7d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT column_name, datatype, description, \"\n",
    "#                         \"unit from TAP_SCHEMA.columns \"\n",
    "#                         \"WHERE table_name = 'dp03_catalogs_10yr.MPCORB'\")\n",
    "# results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80da6bcf-ffc5-4f6e-b687-a0fb0b017517",
   "metadata": {},
   "source": [
    "### 2.1 Create and execute the query\n",
    "\n",
    "Create a query that will retrive the Solar System object ID (`ssObjectid`);\n",
    "five orbital elements eccentricity (`e`), \n",
    "inclination (`incl`; degrees),\n",
    "perihelion distance (`q`; au),\n",
    "longitude of the ascending node (`node`; degrees), \n",
    "and argument of the perihelion (`peri`; degrees);\n",
    "and the Minor Planet Center's measured absolute $H$ and $G$ magnitudes (`mpcH` and `mpcG`).\n",
    "\n",
    "**Orbital elements:**\n",
    "For a refresher on what each orbital element means, the \n",
    "[Wikipedia page for orbital elements](https://en.wikipedia.org/wiki/Orbital_elements)\n",
    "is a good place to start and has an informative diagram.\n",
    "\n",
    "**Absolute $H$ magnitudes:**\n",
    "For Solar System objects, absolute $H$ magnitudes are defined to be for an object 1 au from the Sun and 1 au \n",
    "from the observer, and at a phase angle (the angle Sun-object-Earth) of 0 degrees.\n",
    "Absolute $H$ magnitudes (`mpcH` in the `MPCORB` table) are derived by correcting for distance, fitting a function to the relationship between\n",
    "absolute magnitude and phase (the slope of which is referred to as the $G$ magnitude; `mpcG` in the `MPCORB` table), \n",
    "and evaluating the function at a phase of 0 deg.\n",
    "The $H$ and $G$ magnitudes in different filters can reveal the surface properties and composition of asteroids.\n",
    "The $H$ magnitude is specifically used for objects that primarily reflect sunlight, such as asteroids, since their brightness depends on their distance from the Sun and the observer, \n",
    "and their \"albedo\" (the fraction of light reflected, from 0 to 1; see Section 3).\n",
    "\n",
    "**Minimum semi-major axis:**\n",
    "TNOs, by definition, have a semi-major axis greater than the mean semi-amjor axis of Neptune, 30.1 au (`a_Neptune`).\n",
    "The semi-major axis ($a$) is not included in the `MPCORB` table because it\n",
    "can be derived from the eccentricity and perihelion distance as $a = \\frac{q}{1-e}$.\n",
    "Include a constraint in the query to only return objects with semi-major axes greater than this limit:<br>\n",
    "`WHERE (mpc.q / (1 - mpc.e) > {a_Neptune})`.\n",
    "\n",
    "**Minimum number of observations:**\n",
    "In order to robustly measure the absolute $H$ magnitudes and the orbital elements for a TNO,\n",
    "multiple observations over time are required. \n",
    "The minimum number of observations needed for a given analysis will depend on the science goals\n",
    "and the precision required, and is ultimately left to the discretion\n",
    "of the scientist.\n",
    "Here, as an example, use a limit of 10 observations (`min_num_obs`).\n",
    "The number of observations is stored in the `numObs` column of the `SSObject` table,\n",
    "so including this constraint requires a `JOIN` with this table.\n",
    "\n",
    "**Maximum eccentricity:**\n",
    "Finally, include the constraint that the eccentricity must be less than 1 to reject any \n",
    "hyperbolic orbits (which belong to different populations of small bodies).\n",
    "\n",
    "**Execute the query:**\n",
    "The following query returns 27239 TNOs, which is about an order of magnitude greater than \n",
    "the current number of confirmed TNOs, and is consistent with expectations and forecast for LSST Y10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4aac9a02-b820-4d64-8291-420ca8c645ed",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "min_num_obs = 10\n",
    "query = \"SELECT mpc.ssObjectId, mpc.e, mpc.incl, mpc.q, \" + \\\n",
    "        \"mpc.node, mpc.peri, mpc.mpcH, mpc.mpcG \" + \\\n",
    "        \"FROM dp03_catalogs_10yr.MPCORB as mpc \" + \\\n",
    "        \"JOIN dp03_catalogs_10yr.SSObject as sso \" + \\\n",
    "        \"ON mpc.ssObjectId = sso.ssObjectId \" + \\\n",
    "        f\"WHERE (mpc.q / (1 - mpc.e) > {a_Neptune}) \" + \\\n",
    "        f\"AND sso.numObs > {min_num_obs} \" + \\\n",
    "        \"AND mpc.e < 1 \"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6316e92d-5922-4e43-b9ea-2c84a9ecef75",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)\n",
    "assert job.phase == 'COMPLETED'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e960d84f-2e56-4de9-b018-3050a009b407",
   "metadata": {},
   "outputs": [],
   "source": [
    "tnos = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d500702-6790-4be2-80b4-669d0598b182",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(tnos))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff6a1146-a2aa-406d-91f3-98ecdf6677af",
   "metadata": {},
   "source": [
    "### 2.1 Calculate semi-major axis\n",
    "\n",
    "Calculate the semi-major axis, `a`, for all TNOs and add it to the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5db1fea4-bacd-49c6-be75-b6d160989418",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tnos['a'] = tnos['q'] / (1.0 - tnos['e'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cecd564-65b6-44f7-9b5a-0cf43f574435",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 3. Explore the dynamical classes of TNOs\n",
    "\n",
    "The dynamical classes of TNOs are defined by their orbital parameters, \n",
    "which show clustering in the parameter space of eccentricity or inclination\n",
    "versus semi-major axis.\n",
    "\n",
    "Plot `eccentricity` and `inclination` as a function of semi-major axis `a`, zooming in the range\n",
    "from 28 au to 100 au to better display clustering at certain values of `a`. \n",
    "\n",
    "Explore how the representation changes by adjusting the ranges define din `set_ylim` and `set_xlim`, and re-executing the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "064d22fb-6dce-475b-a404-fdf495887b5a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)\n",
    "\n",
    "ax1.plot(tnos['a'], tnos['e'], '.', ms=2, alpha=0.2, mew=0, color='black')\n",
    "ax1.set_ylabel('Eccentricity')\n",
    "\n",
    "ax2.plot(tnos['a'], tnos['incl'], '.', ms=2, alpha=0.2, mew=0, color='black')\n",
    "ax2.set_ylabel('Inclination (deg)')\n",
    "ax2.set_xlabel('Semi-major axis (au)')\n",
    "\n",
    "ax1.set_ylim([-0.03, 0.7])\n",
    "ax2.set_xlim([28, 100])\n",
    "\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22a3a0d6-0161-43a3-954c-2717acfa386b",
   "metadata": {},
   "source": [
    "> **Figure 1:** Eccentricity (top) and inclination (bottom) as a function of semi-major axis for TNOs.\n",
    "\n",
    "Now use the `seaborn` statistical visualization library to plot the eccentricity and inclination as functions of semi-major axis in a different manner, with a focus on objects having a semi-major axis (`a`) less than 70 au."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18b387df-9cbe-4705-8d55-f9bf5fa67bd0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "sns.set(style=\"ticks\")\n",
    "tmp = tnos[tnos['a'] <= 70]\n",
    "sns.jointplot(x=tmp['a'], y=tmp['e'], kind=\"hex\", color=\"#0077BB\",\n",
    "              gridsize=50, height=5)\n",
    "plt.xlabel('semi-major axis (au)')\n",
    "plt.ylabel('eccentricity')\n",
    "plt.show()\n",
    "remove_figure(fig)\n",
    "\n",
    "fig = plt.figure(figsize=(5, 4))\n",
    "sns.histplot(x=tmp['a'], y=tmp['e'], bins=100, cmap=\"flare\",\n",
    "             cbar=True, cbar_kws=dict(shrink=.95))\n",
    "plt.xlabel('semi-major axis (au)')\n",
    "plt.ylabel('eccentricity')\n",
    "plt.show()\n",
    "remove_figure(fig)\n",
    "\n",
    "fig = plt.figure(figsize=(5, 4))\n",
    "sns.histplot(x=tmp['a'], y=tmp['incl'], bins=100, cmap=\"flare\",\n",
    "             cbar=True, cbar_kws=dict(shrink=.95))\n",
    "plt.xlabel('semi-major axis (au)')\n",
    "plt.ylabel('inclination (degrees)')\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "744002a6-7ecd-4350-b379-cf95521ce79a",
   "metadata": {},
   "source": [
    "> **Figure 2:** Two-dimensional histograms of the orbital elements eccentricity, inclination, and semi-major axis of TNOs.\n",
    "\n",
    "**Why is there clustering in semi-major axis?**\n",
    "\n",
    "The clustering in semi-major axis is due to *resonances*.\n",
    "Since Neptune is so much more massive than a TNO, it influences the regions\n",
    "of orbital element parameter space that are stable, \n",
    "and orbits that are \"in resonance with Neptune\" are stable orbits for TNOs.\n",
    "\n",
    "TNOs with orbital periods that are in, for example, a 1:1 resonance with Neptune\n",
    "have the same period, $P$, and thus same semi-major axis, $a$, because $P^2 \\propto a^3$.\n",
    "TNOs in a 2:1 resonance with Neptune have double the period, and thus $\\sqrt[3]{4}$ times\n",
    "the semi-major axies, which is $\\sim47$ au.\n",
    "The population of 2:1 resonant TNOs at $\\sim47$ will be labeled in the plot in Section 3.2 below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c94ad184-59ac-40dc-9d85-bddb11c0a353",
   "metadata": {},
   "source": [
    "### 3.1 Use known dynamical class boundaries to identify members\n",
    "\n",
    "Below, define the boundaries to classify the TNOs in the following dynamical classes: \n",
    "scattering, detached, classical, and extreme. \n",
    "\n",
    "> **Warning**: In the literature, class definitions may vary!\n",
    "This notebook uses definitions to match the dynamical classifications in [Bernardinelli et al. (2022), ApJS](https://ui.adsabs.harvard.edu/abs/2022ApJS..258...41B/abstract),\n",
    "which found 814 TNOs in the Dark Energy Survey Year 6 data,\n",
    "to recreate their Figure 8 below. Note that the class boundaries are not \"fixed\" but rather \"fluid\" - the cell below uses approximate definitions. The actual dynamical classification process is performed in terms of orbit solutions (integration); see [Gladman et al. 2008](https://ui.adsabs.harvard.edu/abs/2008ssbn.book...43G/abstract) and [Gladman and Volk (2021)](https://ui.adsabs.harvard.edu/abs/2021ARA%26A..59..203G/abstract)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2aaeaa60-82f6-4188-abcd-9c318d46a259",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "max_a_classical = 47\n",
    "min_a_classical = 40\n",
    "max_e_classical = 0.24\n",
    "mask_classical = ((tnos['a'] >= min_a_classical)\n",
    "                  & (tnos['a'] <= max_a_classical)\n",
    "                  & (tnos['e'] <= max_e_classical))\n",
    "\n",
    "min_q_detached = 38\n",
    "mask_detached = ((tnos['a'] > max_a_classical)\n",
    "                 & (tnos['e'] > max_e_classical)\n",
    "                 & (tnos['q'] >= min_q_detached))\n",
    "\n",
    "mask_scattering = ~mask_classical & (tnos['q'] < min_q_detached)\n",
    "\n",
    "min_a_extreme = 150\n",
    "min_q_extreme = 30\n",
    "mask_extreme = (tnos['a'] > min_a_extreme) & (tnos['q'] > min_q_extreme)\n",
    "\n",
    "max_q_comets = 10\n",
    "mask_comets = tnos['q'] < max_q_comets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a0fb75a-d176-498a-9cd8-3c9dd1e8e9bd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-23T04:43:42.971784Z",
     "iopub.status.busy": "2023-07-23T04:43:42.971466Z",
     "iopub.status.idle": "2023-07-23T04:43:42.976329Z",
     "shell.execute_reply": "2023-07-23T04:43:42.975701Z",
     "shell.execute_reply.started": "2023-07-23T04:43:42.971760Z"
    },
    "tags": []
   },
   "source": [
    "Create a new array called `classes` that contains strings labeling the dynamical class for each object. Initialize this array with the string \"Other\" to cover unclassified TNOs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45ff38ad-e7d8-4f6d-8c78-3ef907385531",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "other_string = 'Other     '\n",
    "classes = np.repeat(other_string, len(tnos))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "978e2449-dd8d-4b1d-bac5-eef066ab2eca",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "classes[mask_classical] = 'Classical'\n",
    "classes[mask_detached] = 'Detached'\n",
    "classes[mask_scattering] = 'Scattering'\n",
    "classes[mask_extreme] = 'Extreme'\n",
    "classes[mask_comets] = 'Comets'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98bd9816-77f3-4e92-b7a8-235e007c0d25",
   "metadata": {},
   "source": [
    "\n",
    "Take the whitespace out of the string where `classes` = 'Other'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bdf4596-3c6f-44ca-8e74-5638799a6e3c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tx = np.where(classes == other_string)[0]\n",
    "classes[tx] = 'Other'\n",
    "del tx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5423e164-971f-4e32-a374-6bf9b44db711",
   "metadata": {},
   "source": [
    "Add the contents of the `classes` array as a new column to the TNO table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abf9a96e-720b-4fbf-90ef-241ded618ee1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tnos['class'] = classes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "033e5437-a231-4414-9ec5-fa7a9b68076c",
   "metadata": {},
   "source": [
    "Confirm that the only unique values of the `class` column are the five classes defined above.\n",
    "\n",
    "Print the number of TNOs identified as each class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b47129f-fd16-469d-ad28-63a86c2c56d3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "class_strings, class_counts = np.unique(tnos['class'], return_counts=True)\n",
    "for s, string in enumerate(class_strings):\n",
    "    print('%-12s %6i' % (class_strings[s], class_counts[s]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30e60010-c6d6-48e1-9339-1b7f4e5a9bb5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-25T04:44:26.373430Z",
     "iopub.status.busy": "2023-07-25T04:44:26.372614Z",
     "iopub.status.idle": "2023-07-25T04:44:26.834260Z",
     "shell.execute_reply": "2023-07-25T04:44:26.832870Z",
     "shell.execute_reply.started": "2023-07-25T04:44:26.373401Z"
    },
    "tags": []
   },
   "source": [
    "### 3.2 Visualize the dynamical classes\n",
    "\n",
    "Define five colors, markers, and line styles to use, one set for each class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04932e52-2f8a-4e95-bb2d-eb027a00915d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "class_to_format = {'Other': ('xkcd:black', 'x', (0, (1, 10))),\n",
    "                   'Classical': ('xkcd:light purple', 's', 'dotted'),\n",
    "                   'Detached': ('xkcd:light orange', '+', 'dashed'),\n",
    "                   'Extreme': ('xkcd:scarlet', 'D', 'dashdot'),\n",
    "                   'Scattering': ('xkcd:cerulean blue', '.', 'solid'),\n",
    "                   'Comets': ('xkcd:teal', '^', (5, (10, 3)))}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ebb1759-e0a2-438e-b866-24ff447eb0d8",
   "metadata": {},
   "source": [
    "Recreate Figure 8 from \n",
    "[Bernardinelli et al. (2022), ApJS](https://ui.adsabs.harvard.edu/abs/2022ApJS..258...41B/abstract),\n",
    "but impose a limit on semi-major axis of 100 au which excludes the \"Extreme\" dynamical class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "418000e0-6ed4-4107-b4ec-dad2c4eeca82",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)\n",
    "fig.set_size_inches(9, 7)\n",
    "plt.subplots_adjust(wspace=0, hspace=0)\n",
    "plt.rcParams[\"image.composite_image\"] = False\n",
    "e = np.linspace(0, 0.99999, 10000)\n",
    "\n",
    "for q in [10, 30, 35, 40, 45, 50, 55]:\n",
    "    ax1.plot(q / (1 - e), e, 'k--', alpha=0.3)\n",
    "\n",
    "for i, j in zip([1, 2, 3, 3, 4, 4, 5, 5, 11],\n",
    "                [1, 1, 1, 2, 1, 3, 2, 1, 2]):\n",
    "    ax1 = plotResonances(ax1, i, j, alpha=0.4)\n",
    "    ax2 = plotResonances(ax2, i, j, y_loc=55, alpha=0.4)\n",
    "\n",
    "gradientFill(30 / (1 - e), e, ax=ax1, alpha=0.2, invert=True)\n",
    "color, marker, _ = class_to_format['Detached']\n",
    "gradientFill(min_q_detached / (1 - e)[e > max_e_classical],\n",
    "             e[e > max_e_classical],\n",
    "             ax=ax1, alpha=0.4, color=color, marker=marker)\n",
    "\n",
    "scat = np.append(30 / (1 - e), 35 / (1 - e)[::-1])\n",
    "scat_e = np.append(e, e[::-1])\n",
    "color, marker, _ = class_to_format['Scattering']\n",
    "gradientFill(scat, scat_e, ax=ax1, alpha=0.4,\n",
    "             color=color, marker=marker)\n",
    "\n",
    "cl = np.append(35 / (1 - e)[e < max_e_classical], 100)\n",
    "e_cl = np.append(e[e < max_e_classical], 0.24)\n",
    "color, marker, _ = class_to_format['Classical']\n",
    "gradientFill(cl, e_cl, ax=ax1, alpha=0.4, color=color, marker=marker)\n",
    "\n",
    "for c in class_to_format.keys():\n",
    "    b = tnos[tnos['class'] == c]\n",
    "    use_alpha = 0.3\n",
    "    if c == 'Detached':\n",
    "        use_alpha = 1\n",
    "    color, marker, _ = class_to_format[c]\n",
    "    ax1.plot(b['a'], b['e'], marker, alpha=use_alpha, color=color,\n",
    "             markersize=3)\n",
    "    ax2.plot(b['a'], b['incl'], marker, alpha=use_alpha, color=color,\n",
    "             markersize=3)\n",
    "\n",
    "fontsize_labels = 12\n",
    "ax1.text(45, 0.81, r'$\\bf{Long-period\\ comets}$', rotation=0,\n",
    "         color=class_to_format['Comets'][0], fontsize=fontsize_labels)\n",
    "ax1.text(40, 0.6, r'$\\bf{Scattering}$', rotation=0,\n",
    "         color=class_to_format['Scattering'][0], fontsize=fontsize_labels)\n",
    "ax1.text(85, 0.3, r'$\\bf{Detached}$', rotation=0,\n",
    "         color=class_to_format['Detached'][0], fontsize=fontsize_labels)\n",
    "ax1.text(65, 0.1, r'$\\bf{Classical}$', rotation=0,\n",
    "         color=class_to_format['Classical'][0], fontsize=fontsize_labels)\n",
    "ax1.text(49, 0.12, r'$\\bf{Other}$', rotation=0,\n",
    "         color=class_to_format['Other'][0], fontsize=fontsize_labels)\n",
    "\n",
    "fontsize_q_text = 9\n",
    "ax1.text(70, 0.23, r'$55$ au', rotation=25, fontsize=fontsize_q_text)\n",
    "ax1.text(70, 0.3, r'$50$ au', rotation=23, fontsize=fontsize_q_text)\n",
    "ax1.text(70, 0.37, r'$45$ au', rotation=24, fontsize=fontsize_q_text)\n",
    "ax1.text(70, 0.58, r'$30$ au', rotation=16, fontsize=fontsize_q_text)\n",
    "ax1.text(70, 0.87, r'$10$ au', rotation=3, fontsize=fontsize_q_text)\n",
    "\n",
    "ax2.set_xlabel('Semi-major axis (au)')\n",
    "ax1.set_ylabel('Eccentricity')\n",
    "ax2.set_ylabel('Inclination (deg)')\n",
    "\n",
    "ax1.set_ylim(0, 1)\n",
    "ax2.set_ylim(0, 65)\n",
    "ax1.set_xlim(20, 100)\n",
    "ax2.set_xlim(20, 100)\n",
    "\n",
    "ax1.set_xticks(np.arange(20, 110, 10))\n",
    "ax1.set_yticks(np.arange(0, 1.1, 0.1))\n",
    "ax2.set_yticks(np.arange(0, 70, 10))\n",
    "\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8323d329-5a28-4b46-ba0d-34e8a2d6f341",
   "metadata": {},
   "source": [
    "> **Figure 3:** Eccentricity (top) and inclination (bottom) versus semi-major axis for DP0.3 TNOs. This plot recreates Figure 8 from \n",
    "[Bernardinelli et al. (2022), ApJS](https://ui.adsabs.harvard.edu/abs/2022ApJS..258...41B/abstract).\n",
    "\n",
    "Extend the previous plot in semi-major axies to include \"extreme\" TNOs and more long-period comets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53ef9bdf-3b85-4183-8526-5b597a16273f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True)\n",
    "fig.set_size_inches(15, 7)\n",
    "plt.subplots_adjust(wspace=0, hspace=0)\n",
    "plt.rcParams[\"image.composite_image\"] = False\n",
    "e = np.linspace(0, 0.9999999, 10000)\n",
    "\n",
    "for q in [10, 30, 35, 40, 45, 50, 55]:\n",
    "    ax1.plot(q / (1 - e), e, 'k--', alpha=0.3)\n",
    "    ax2.plot(q / (1 - e), e, 'k--', alpha=0.3)\n",
    "\n",
    "for i, j in zip([1, 2, 3, 3, 4, 4, 5, 5, 11],\n",
    "                [1, 1, 1, 2, 1, 3, 2, 1, 2]):\n",
    "    ax1 = plotResonances(ax1, i, j, 0.7, alpha=0.4)\n",
    "\n",
    "gradientFill(30 / (1 - e), e, ax=ax1, alpha=0.3, invert=True)\n",
    "\n",
    "a_det = 35 / (1 - e)[e > max_e_classical]\n",
    "e_det = e[e > max_e_classical]\n",
    "color, marker, _ = class_to_format['Detached']\n",
    "gradientFill(a_det, e_det, ax=ax1,\n",
    "             alpha=0.4,\n",
    "             color=color, marker=marker)\n",
    "gradientFill(a_det[a_det < 150], e_det[a_det < 150], ax=ax2,\n",
    "             alpha=0.4,\n",
    "             color=color, marker=marker)\n",
    "\n",
    "scat = np.append(30/(1-e), 35/(1-e)[::-1])\n",
    "scat_e = np.append(e, e[::-1])\n",
    "color, marker, _ = class_to_format['Scattering']\n",
    "gradientFill(scat, scat_e, ax=ax1, alpha=0.4,\n",
    "             color=color, marker=marker)\n",
    "gradientFill(scat, scat_e, ax=ax2, alpha=0.4,\n",
    "             color=color, marker=marker)\n",
    "gradientFill(30 / (1 - e), e, ax=ax2, alpha=0.4, invert=True,\n",
    "             color=color, marker=marker)\n",
    "\n",
    "cl = np.append(35 / (1 - e)[e < 0.24], 150)\n",
    "msk = cl < 151\n",
    "e_cl = np.append(e[e < 0.24], 0.24)\n",
    "color, marker, _ = class_to_format['Classical']\n",
    "gradientFill(cl, e_cl, ax=ax1, alpha=0.3,\n",
    "             color=color, marker=marker)\n",
    "gradientFill(cl[msk], e_cl[msk], ax=ax2, alpha=0.3,\n",
    "             color=color, marker=marker)\n",
    "\n",
    "a_ext = min_q_extreme / (1 - e)\n",
    "color, marker, _ = class_to_format['Extreme']\n",
    "gradientFill(np.append(150, a_ext[a_ext > 150]),\n",
    "             np.append(0, e[a_ext > 150]), ax=ax2,\n",
    "             alpha=0.3, color=color, marker=marker)\n",
    "\n",
    "for c in class_to_format.keys():\n",
    "    b = tnos[tnos['class'] == c]\n",
    "    use_alpha = 0.9\n",
    "    color, marker, _ = class_to_format[c]\n",
    "    ax1.plot(b['a'], b['e'], marker, alpha=use_alpha, color=color)\n",
    "    ax2.plot(b['a'], b['e'], marker, alpha=use_alpha, color=color)\n",
    "\n",
    "fontsize_labels = 12\n",
    "ax1.text(55, 0.80, r'$\\bf{Long-period\\ comets}$', rotation=0,\n",
    "         color=class_to_format['Comets'][0], fontsize=fontsize_labels)\n",
    "ax1.text(40, 0.6, r'$\\bf{Scattering}$', rotation=0,\n",
    "         color=class_to_format['Scattering'][0], fontsize=fontsize_labels)\n",
    "ax1.text(75, 0.3, r'$\\bf{Detached}$', rotation=0,\n",
    "         color=class_to_format['Detached'][0], fontsize=fontsize_labels)\n",
    "ax1.text(65, 0.1, r'$\\bf{Classical}$', rotation=0,\n",
    "         color=class_to_format['Classical'][0], fontsize=fontsize_labels)\n",
    "ax1.text(49, 0.12, r'$\\bf{Other}$', rotation=0,\n",
    "         color=class_to_format['Other'][0], fontsize=fontsize_labels)\n",
    "ax2.text(250, 0.5, r'$\\bf{Extreme}$', rotation=0,\n",
    "         color=class_to_format['Extreme'][0], fontsize=fontsize_labels)\n",
    "\n",
    "fontsize_q_text = 9\n",
    "ax1.text(70, 0.23, r'$55$ au', rotation=25, fontsize=fontsize_q_text)\n",
    "ax1.text(70, 0.3, r'$50$ au', rotation=23, fontsize=fontsize_q_text)\n",
    "ax1.text(70, 0.37, r'$45$ au', rotation=24, fontsize=fontsize_q_text)\n",
    "ax1.text(70, 0.58, r'$30$ au', rotation=16, fontsize=fontsize_q_text)\n",
    "ax1.text(70, 0.87, r'$10$ au', rotation=3, fontsize=fontsize_q_text)\n",
    "\n",
    "fig.supxlabel('Semi-major axis (au)')\n",
    "ax1.set_ylabel('Eccentricity')\n",
    "ax1.set_ylim(0, 1.01)\n",
    "ax1.set_xlim(20, 100)\n",
    "ax2.set_xlim(100, 500)\n",
    "ax1.set_yticks(np.arange(0, 1.1, 0.1))\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3170fc9f-1d3a-4367-b105-fa02ead9ed6f",
   "metadata": {},
   "source": [
    "> **Figure 4:** Similar to the top panel of Figure 3, but with the addition of the \"extreme\" TNOs with semi-major axes > 150 au."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1c9830b-dea5-4d77-b153-49aed0287496",
   "metadata": {},
   "source": [
    "### 3.3 Characterize the dynamical classes\n",
    "\n",
    "For the \"Classical\", \"Detached\", and \"Scattering\" classes of TNOs, compare distributions of their\n",
    "orbital elements, estimated diameters, and absolute $H$ magnitudes from the `MPCORB` catalog.\n",
    "\n",
    "Note that in DP0.3 data, the `SSObject` table does not contain any _LSST-measured_ absolute $H$ magnitudes\n",
    "in the LSST's $ugrizy$ filters.\n",
    "The reason why (noisy phase curves) will be demonstrated in Section 5."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f066850f-cedf-4aef-ba6e-f2e5ab262e16",
   "metadata": {},
   "source": [
    "#### 3.3.1 Orbital elements, diameters, and magnitudes\n",
    "\n",
    "**Orbital Elements** \n",
    "\n",
    "Plot histograms of the two orbital elements that have not yet been explored:\n",
    "the longitude of the ascending node (`node`; degrees) and the argument of the perihelion (`peri`; degrees). These will be \"Histogram 1\" and \"Histogram 2\" below.\n",
    "\n",
    "Use a logarithmic scale in the y-axis to simultaneously show the different populations.\n",
    "\n",
    "Notice that the distributions for both the `node` and the `peri` for \"Classical\", \"Detached\", \"Scattering\",  and \"Other\" are all quite flat, but that the distribution for the \"Extreme\" objects  seem to be clustered at certain values. \n",
    "\n",
    "The apparent clustering in longitude of perihelion `node` and ascending node `peri` of \"extreme\" TNOs motivated the hypothesis that the solar system contains a 5–10 Earth-mass planet (Planet X/Planet 9) at 400–800 times Earth's distance from the Sun ([Trujillo & Sheppard 2014](https://ui.adsabs.harvard.edu/abs/2014Natur.507..471T/abstract); [Batygin & Brown 2016](https://ui.adsabs.harvard.edu/abs/2016AJ....151...22B/abstract); [Batygin et al. 2019](https://ui.adsabs.harvard.edu/abs/2019PhR...805....1B/abstract)). \n",
    "\n",
    "However, this hypothesis is still under debate (e.g., [Napier et al. 2021](https://ui.adsabs.harvard.edu/abs/2021PSJ.....2...59N/abstract)), and a statistical analysis that demonstrates whether this apparent clustering could be attributed to small-numbers of objects or not is beyond the scope of this tutorial. \n",
    "\n",
    "Additionally, caution must be exercised when interpreting these results, as the observed clustering could potentially stem from a selection bias. It may be an artifact of the type of model used as inputs for Trans-Neptunian Objects (TNOs) in the \"MPCORB\" catalog. This catalog is a hybrid one that employs a pre-defined population model known as the \"Synthetic Solar System Model, S3M\" [Grav et al. 2011](https://iopscience.iop.org/article/10.1086/659833/meta), wherein model objects are replaced with real, existing objects currently listed in the Minor Planet Center. Consequently, the statistics are influenced by the S3M properties. While the model adequately produces populations for TNOs, it does not fully capture the properties of certain solar system objects, such as \"extreme\" TNOs, many of which have been discovered since 2011 (after the paper's publication) thanks to surveys like _OSSOS_ or _DES_.\n",
    "\n",
    "**Estimated diameters**\n",
    "\n",
    "Estimate the diameter of all TNOs using an assumed mean albedo of 0.1 (see [Vilenius et al. 2012](https://ui.adsabs.harvard.edu/abs/2012A%26A...541A..94V/abstract))  and plot the estimated size distribution for these TNOs (\"Hisogram 3\" below). Note, however, that some dynamical families have higher albedos than others ([Vilenius et al. 2012](https://ui.adsabs.harvard.edu/abs/2012A%26A...541A..94V/abstract)).\n",
    "\n",
    "\n",
    "Notice that the x-axis limit is set to 400 km.\n",
    "Adjust the limits to explore TNOs with larger estimated radii.\n",
    "\n",
    "Compared to the \"classical\" population, the \"detached\" objects are expected to have larger sizes (smaller H magnitude), whereas the \"scattering\" population is expected to have smaller sizes (larger H magnitudes). \n",
    "\n",
    "**Absolute Magnitudes**\n",
    "\n",
    "Plot the cumulative distribution of MPC absolute $H$ magnitudes for these TNOs (\"Histogram 4\").\n",
    "\n",
    "Notice the tail to bright magnitudes for detached objects, similar as the tail to objects with larger size estimates above.\n",
    "\n",
    "As in section 5.1 of [Bernardinelli et al. (2022), ApJS](https://ui.adsabs.harvard.edu/abs/2022ApJS..258...41B/abstract), cumulative plots such as the one in the next cell can be used to test models of TNO formation, such as the CFEPS-L7 model ([Kavelaars et al. 2009](https://ui.adsabs.harvard.edu/abs/2009AJ....137.4917K/abstract); [Petit et al.\n",
    "2011](https://ui.adsabs.harvard.edu/abs/2011AJ....142..131P/abstract); [Gladman et al. 2012](https://ui.adsabs.harvard.edu/abs/2012AJ....144...23G/abstract)) for \"Classical\" TNOs; see Fig. 10 of [Bernardinelli et al. (2022), ApJS](https://ui.adsabs.harvard.edu/abs/2022ApJS..258...41B/abstract).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02e70ec6-4c0c-4471-98e8-9b4bdf39e658",
   "metadata": {},
   "source": [
    "Exclude the comets, which belong to a different type of solar system objects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e85c8ddb-10f8-42b6-bfef-b6d39f793f5b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "classes_tnos = ['Other', 'Classical', 'Detached', 'Extreme', 'Scattering']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7f164b3-3977-497e-85ae-1cd5daba395b",
   "metadata": {},
   "source": [
    "Estimate the diameter of all TNOs using an assumed mean albedo of 0.1 (see [Vilenius et al. 2012](https://ui.adsabs.harvard.edu/abs/2012A%26A...541A..94V/abstract) for measurements of and discussion about TNO's geometric albedos), which is the default in the `estimateDiameter` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "189496c4-a7c1-4b24-b20a-549c72b37552",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tnos['d'] = estimateDiameter(tnos['mpcH'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0d3ee7d-9d8f-441c-b1e8-248af197f8b8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12, 10))\n",
    "\n",
    "createSubplot(2, 2, 1, 'Histogram 1', 'longitude of ascending node (deg)',\n",
    "              'log number of TNOs')\n",
    "for label in classes_tnos:\n",
    "    nbins = 10 if label == 'Extreme' else 25\n",
    "    color, _, linestyle = class_to_format[label]\n",
    "    plotHistogram(tnos['node'], label, nbins, color, linestyle, log=True)\n",
    "\n",
    "createSubplot(2, 2, 2, 'Histogram 2', 'argument of perihelion (deg)',\n",
    "              'log number of TNOs')\n",
    "for label in classes_tnos:\n",
    "    nbins = 10 if label == 'Extreme' else 25\n",
    "    color, _, linestyle = class_to_format[label]\n",
    "    plotHistogram(tnos['peri'], label, nbins, color, linestyle, log=True)\n",
    "\n",
    "createSubplot(2, 2, 3, 'Histogram 3', 'estimated diameter (km)',\n",
    "              'fraction of TNOs')\n",
    "for label in classes_tnos:\n",
    "    nbins = 30 if label == 'Extreme' else 150\n",
    "    color, _, linestyle = class_to_format[label]\n",
    "    plotHistogram(tnos['d'], label, nbins, color, linestyle, density=True)\n",
    "plt.xlim([0, 400])\n",
    "\n",
    "createSubplot(2, 2, 4, 'Histogram 4', 'MPC absolute H magnitude',\n",
    "              'fraction of TNOs')\n",
    "for label in classes_tnos:\n",
    "    nbins = 5 if label == 'Extreme' else 300\n",
    "    color, _, linestyle = class_to_format[label]\n",
    "    plotHistogram(tnos['mpcH'], label, nbins, color, linestyle, density=True,\n",
    "                  cumulative=True)\n",
    "plt.legend(loc='upper left')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b180c3b-eb3c-467a-b6be-7c5aa85d2429",
   "metadata": {},
   "source": [
    "> **Figure 5:** Histograms of the orbital elements, longitude of the ascending node (top left) and argument of perihelion (top right), and object properties, the estimated diameter (bottom left) and absolute H magnitude (bottom right), for five categores of TNOs (see legend in lower-right panel)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0ddc8b7-6026-4446-a6ac-7467cd622135",
   "metadata": {},
   "source": [
    "#### 3.3.2 Absolute G magnitude\n",
    "\n",
    "Print to screen the minimum and maximum values of \n",
    "the MPC absolute $G$ magnitudes, and notice that the MPC assumes a slope \n",
    "of $G=0.15$ for all DP0.3 TNOs.\n",
    "It is not worth creating a plot of the distribution for absolute $G$ magnitudes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62a34f7f-fbb5-470d-9271-bdbd6bbc90c8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for label in class_to_format.keys():\n",
    "    tx = np.where(tnos['class'] == label)[0]\n",
    "    print('%-12s   min G = %4.2f   max G = %4.2f' %\n",
    "          (label, np.min(tnos['mpcG'][tx]), np.max(tnos['mpcG'][tx])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10f88136-da46-468c-b0c0-4a5740fcc633",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### 3.3.3 Summary\n",
    "\n",
    "The distributions of TNO orbital parameters and properties, such as the absolute magnitude $H$, offer valuable insights into the physical characteristics and evolution of various TNO dynamical classes. Additionally, they contribute to our understanding of the formation and evolution of the Solar System as a whole. While a detailed interpretation of the reasons behind these distribution differences among the various TNO dynamical classes and how they relate to different formation models exceeds the scope of this tutorial, it remains an ongoing area of research in planetary science. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37663b79-daa5-45ad-ba1d-3a1e6023a90f",
   "metadata": {},
   "source": [
    "## 4. Visualize orbits for each dynamical class\n",
    "\n",
    "Use the Minor Planet Center-derived orbital elements from the TNOs to visualize their orbits, \n",
    "and look for patterns in orbits of different dynamical classes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbb6a5d2-1428-4bb6-89fb-2939963c40bb",
   "metadata": {},
   "source": [
    "### 4.1 Calculate the full orbital ellipse per object\n",
    "\n",
    "Calculate the full ellipse of the orbit for each TNO, based on their\n",
    "orbital elements, using the `fullEllipse` function. \n",
    "Add the full ellipse to the `tnos` table as a new column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "963f5c42-0e54-4e22-aa9b-acdc8dcc54b6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "temp = []\n",
    "for tno in tnos:\n",
    "    temp.append(fullEllipse(tno['a'], tno['e'], tno['incl'],\n",
    "                            tno['node'], tno['peri']))\n",
    "tnos['fullEllipse'] = temp\n",
    "del temp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5159696-2cdc-4262-b178-5cb5ea0307b2",
   "metadata": {},
   "source": [
    "Plot the full ellipse for the first 100 objects, as an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b06f6f7-dc26-4e41-83bf-1ebca91ba231",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 5))\n",
    "for tno in tnos[:100]:\n",
    "    plt.plot(tno['fullEllipse'][:, 0], tno['fullEllipse'][:, 1],\n",
    "             alpha=0.3)\n",
    "    plt.xlabel(\"X (au)\")\n",
    "    plt.ylabel(\"Y (au)\")\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f7e295b-7260-4e65-ade3-5387e7438fd4",
   "metadata": {},
   "source": [
    "> **Figure 6:** The projection of the heliocentric orbits of 100 TNOs into the XY plane of the solar system."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b79038d-f78e-42f0-9ccb-b0a252539535",
   "metadata": {},
   "source": [
    "### 4.2 Illustrate the orbits of each dynamical class\n",
    "\n",
    "Plot the first 10 examples of orbits in each dynamical class in a 2x2 grid. Note how the orbits for the \"Classical\" objects are more circular, whereas the orbits for the other classes are more eccentric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85483a2a-c004-402b-b82c-4fc3eb513ff0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, 2, figsize=(8, 7))\n",
    "plt.rcParams[\"image.composite_image\"] = False\n",
    "\n",
    "classes = ['Classical', 'Scattering', 'Detached', 'Extreme']\n",
    "positions = [(0, 0), (0, 1), (1, 0), (1, 1)]\n",
    "\n",
    "for use_class, (row, col) in zip(classes, positions):\n",
    "    color, _, linestyle = class_to_format[use_class]\n",
    "    current_ax = ax[row, col]\n",
    "    current_ax.set_title(f'{use_class}')\n",
    "    current_ax.set_xlabel('X (au)')\n",
    "    current_ax.set_ylabel('Y (au)')\n",
    "    # Plot up to 10 instances of the current class\n",
    "    for tno in tnos[tnos['class'] == use_class][:10]:\n",
    "        current_ax.plot(tno['fullEllipse'][:, 0], tno['fullEllipse'][:, 1],\n",
    "                        color=color, linestyle=linestyle, alpha=0.3,\n",
    "                        label=use_class)\n",
    "    current_ax.set_aspect('equal')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "remove_figure(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8935d43f-651a-4aeb-a96f-a7000f075eb8",
   "metadata": {},
   "source": [
    "> **Figure 7:** The projection of the heliocentric orbits of 10 TNOs per class into the XY plane of the solar system: classical (upper left), scattering (upper right), detached (lower left), and extreme (lower right)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9655a9d3-4ece-49d8-be75-a6b2156e8d60",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 5. Visualize a TNO in the deep drilling field\n",
    "\n",
    "Since TNOs are located at a considerable distance from Earth, they exhibit slow motion across the sky, with shifts as small as a few degrees over the 10-year duration of the LSST. The [LSST Deep Drilling Fields](https://www.lsst.org/scientists/survey-design/ddf), each covering an area equivallent to the field-of-view of the LSST Camera (9.6 square degrees), receive thousands of visits or observations throughout the 10-year survey period. Consequently, a single TNO can accumulate a substantial number of observations due to its extended visibility over this extended timeframe. In this section we will follow Section 2 of DP0.3 Portal Tutorial 3, and select a particular TNO from the `dp03_catalogs_10yr.DiaSource` catalog to visualize its orbit in space as a function of time. \n",
    "\n",
    "First, we'll do a query that will get all the TNOs with the most number of observations, singling out a particular TNO with `ssObjectId = -735085100561880491`. We will also select `midPointMjdTai`, which is the effective mid-exposure time for this source (in days) to illustrate how the object moves across the sky.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75174782-c60f-41d6-a901-0e538ae20498",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT ra, dec, midPointMjdTai \" + \\\n",
    "        \"FROM dp03_catalogs_10yr.DiaSource \" + \\\n",
    "        \"WHERE ssObjectId = -735085100561880491\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d141cef5-93e1-4774-905a-969a37da754e",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)\n",
    "assert job.phase == 'COMPLETED'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0aa867f-6661-4d15-9386-6efe04c717f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_object = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4ef7a5f-d793-44d7-b68d-efacf9b9d636",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.scatter(my_object['ra'], my_object['dec'],\n",
    "                  c=my_object['midPointMjdTai'],\n",
    "                  cmap='viridis', alpha=0.5, edgecolors='none')\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Dec (deg)')\n",
    "plt.gca().invert_xaxis()\n",
    "plt.colorbar(label='midPointMjdTai')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e872f146-77c9-47bb-9ffc-beaa42f88836",
   "metadata": {},
   "source": [
    "> **Figure 8:** The coordinates of declination (Dec) versus right ascension (RA) in decimal degrees for observations of a TNO in the deep drilling field. Points are colored by the time of the observation, `midPointMjdTai` (colorbar at right). The loops in the object's path on the sky is a result of Earth's orbital period and the 10-year LSST duration."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c34c9f68-5e03-4ed0-bf37-a614e914960c",
   "metadata": {},
   "source": [
    "## 6. Exercises for the learner\n",
    "\n",
    "1. In section 3, the \"Resonant\" TNOs were labeled among the \"Others\" class. Define cuts to identify these resonant objects based on the resonances shown in the the `e-incl-a` plot of section 3.2.\n",
    "Assign them a new class label within the TNO table, and visualize them in the `e-incl-a` plot, as in Fig. 8 of [Bernardinelli et al. (2022), ApJS](https://ui.adsabs.harvard.edu/abs/2022ApJS..258...41B/abstract).\n",
    "\n",
    "2. In section 3.3.2 we used the absolute magnitude $H$ - size relation assuming a mean albedo of 0.1 for all dynamical classes.\n",
    "[Vilenius et al. 2012](https://ui.adsabs.harvard.edu/abs/2012A%26A...541A..94V/abstract) provide mode detailed distributions and albedo measurenets. Try estimating TNOS sizes with different albedos.\n",
    "\n",
    "3. Section 3.3.3 calculated the cumulative distribution of $H$.\n",
    "Calculate the cumulative distributions for other orbital parameters and compare to section 5.1 of [Bernardinelli et al. (2022), ApJS](https://ui.adsabs.harvard.edu/abs/2022ApJS..258...41B/abstract)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}