DP03_02_Main_Belt_Asteroids.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "# Properties of Main Belt Asteroids in DP0.3\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): Jeff Carlin <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: beginner <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da1a210-d858-42fe-8591-570965b8be1a",
   "metadata": {},
   "source": [
    "**Description:** An introduction to Main Belt asteroids in the DP0.3 dataset, and illustration of some simple analyses of their orbital families."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0baf5-51ad-40ec-8991-060a7b27c289",
   "metadata": {},
   "source": [
    "**Skills:** Use TAP queries to retrieve Solar System objects. Plot properties and orbits of Main Belt asteroids."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393da88f-7978-4920-aa4a-a9830df6eed9",
   "metadata": {},
   "source": [
    "**LSST Data Products:** TAP DP0.3 MPCORB, SSObject, and SSSource tables."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c67fab9-136a-4adc-bb42-142b91ab69dd",
   "metadata": {},
   "source": [
    "**Packages:** lsst.rsp.get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f72b27f",
   "metadata": {},
   "source": [
    "**Credit:** Originally developed by Jeff Carlin and the Rubin Community Science Team in the context of the Rubin DP0."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28e91cbf-ab7f-4e26-9276-b00299d6065e",
   "metadata": {},
   "source": [
    "**Get Support:**\n",
    "Find DP0-related documentation and resources at <a href=\"https://dp0.lsst.io\">dp0.lsst.io</a>. 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. Rubin staff will respond to all questions posted there."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfc73be0",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This tutorial demonstrates how to use TAP to query the MPCORB and SSObject tables. The tables are queried simultaneously and joined on their SSObjectId, which allows one to retrieve the measured magnitudes from SSObject and orbital properties from MPCORB. We examine the properties of a large sample of Solar System objects from the DP0.3 catalogs, then focus on asteroids in the Main Belt (between about 1.6-4.2 au). We examine the orbital properties of Main Belt asteroids, and plot orbits of some of them to look for differences between the orbital families.\n",
    "\n",
    "A detailed introduction to the DP0.3 data products can be found in Notebook DP03_01_Introduction_to_DP03."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc36f107",
   "metadata": {},
   "source": [
    "### 1.1 Package Imports\n",
    "\n",
    "The `matplotlib` and `numpy` libraries are widely used Python libraries for plotting and scientific computing. We will use these packages below, including the `matplotlib.pyplot` plotting sublibrary.\n",
    "\n",
    "We also use the `lsst.rsp` package to access the TAP service and query the DP0 catalogs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cddc1458",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from lsst.rsp import get_tap_service"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a32b8d50-446a-4458-a0e7-c3d11a79fed6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"ssotap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c217adff-25ed-4fce-95e7-8aa04630f6cc",
   "metadata": {},
   "source": [
    "### 1.2 Define Functions and Parameters\n",
    "\n",
    "Define a function to convert a given perihelion distance (`q`) and eccentricity (`e`) to an orbital semi-major axis (`a`). Their relationship is defined by q = a(1-e)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eff03d7-c3c6-4843-8e6a-2a507947f9f2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def calc_semimajor_axis(q, e):\n",
    "    \"\"\"\n",
    "    Given a perihelion distance and orbital eccentricity,\n",
    "      calculate the semi-major axis of the orbit.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    q: float\n",
    "        Distance at perihelion, in au.\n",
    "    e: float\n",
    "        Orbital eccentricity.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    a: float\n",
    "        Semi-major axis of the orbit, in au.\n",
    "        q = a(1-e), so a = q/(1-e)\n",
    "    \"\"\"\n",
    "\n",
    "    return q / (1.0 - e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee4fe44f-03c7-4cb6-bbbf-d24a8fc430ae",
   "metadata": {},
   "source": [
    "Define a function to calculate the semi-major axis (a) at which a mean-motion resonance with a planet occurs via Kepler's Third Law: (P/P_planet)^2 = (a/a_planet)^3, where `a` is the semi-major axis of a small body in au, `a_planet` is the planet's semi-major axis in au, `P` is the small body's orbital period in yr, and `P_planet` is the planet's orbital period in yr."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2aa99ed-040b-4da1-b956-de33f359cf70",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_mean_motion_res_a(P_over_P_planet, a_planet):\n",
    "    \"\"\"\n",
    "    \n",
    "    Implements Kepler's Third law,\n",
    "    `(P/P_planet)^2 = (a/a_planet)^3`, where\n",
    "    `a` is the semi-major axis in au,\n",
    "    `a_planet` is the planet's semi-major axis in au,\n",
    "    `P` is the small body's orbital period in yr, and\n",
    "    `P_planet` is the planet's orbital period in yr,\n",
    "    to find the semimajor axis at which a mean-motion\n",
    "    resonance with a planet occurs, via\n",
    "    `a = a_planet * (P/P_planet)^(2/3)`.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    P_over_P_planet: ndarray\n",
    "        Integer ratio of a small body's orbital period\n",
    "        to a planet's orbital period.\n",
    "    a_planet: ndarray\n",
    "        Semimajor axis of a planet, in au.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    a_res: ndarray\n",
    "        Semimajor axis of a mean motion resonance, in au.\n",
    "        a_res = (P/P_planet)^(2/3) * a_planet\n",
    "    \"\"\"\n",
    "\n",
    "    return (float(P_over_P_planet))**0.66667 * a_planet"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f4f0362-0ee6-48af-81e1-7863925bdb22",
   "metadata": {},
   "source": [
    "Define a function to take an `ssObjectId` as an input, and return all measurements of that object from the `SSSource` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb494a0a-a9c8-47db-8678-e95ab34576cb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def get_obj_matched(objid):\n",
    "    \"\"\"\n",
    "    Given an objectId, retrieve all measurements from\n",
    "      SSSource for that object.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    objid: int64\n",
    "        ObjectId of the object of interest.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    allsrc: dataframe\n",
    "        Measurements for the object from the DP0.3 SSSource table.\n",
    "    \"\"\"\n",
    "\n",
    "    query = \"\"\"SELECT * FROM dp03_catalogs_10yr.DiaSource as diasrc \\\n",
    "            JOIN dp03_catalogs_10yr.SSSource as sssrc \\\n",
    "            ON diasrc.diaSourceId = sssrc.diaSourceId \\\n",
    "            WHERE diasrc.ssObjectId = \"\"\" + str(objid)\n",
    "\n",
    "    job = service.submit_job(query)\n",
    "    job.run()\n",
    "    job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "    assert job.phase == 'COMPLETED'\n",
    "\n",
    "    allsrc = job.fetch_result().to_table().to_pandas()\n",
    "\n",
    "    return allsrc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23ee4cc8-f4d4-46ef-9213-4e5159611039",
   "metadata": {},
   "source": [
    "Define a function to plot orbits of multiple Solar System objects on the same figure. The orbits are plotted in XYZ Cartesian coordinates, either in a heliocentric (Sun-centered) or topocentric (relative to the Earth's surface) frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21b96758-0aa6-41db-adb5-b11c37b56a94",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def xyz_orbit_plot_manyobjects(tobs_list, x_list, y_list, z_list, title=None):\n",
    "    \"\"\"\n",
    "    Given lists (of lists) of observation times and XYZ coordinates,\n",
    "      plot the orbits of multiple SS objects.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    tobs_list: list of lists\n",
    "        MidpointTAI times of each observation of each object.\n",
    "    {xyz}_list: list of lists\n",
    "        Solar system positions of each observation of each object.\n",
    "        Can be either heliocentric or topocentric.\n",
    "    title: string\n",
    "        (Optional) title to add to the top of the plot.\n",
    "    \"\"\"\n",
    "\n",
    "    params = {\n",
    "       'axes.labelsize': 12,\n",
    "       'font.size': 12,\n",
    "       'legend.fontsize': 8,\n",
    "       'xtick.major.width': 2,\n",
    "       'xtick.minor.width': 1,\n",
    "       'xtick.major.size': 6,\n",
    "       'xtick.minor.size': 4,\n",
    "       'xtick.direction': 'in',\n",
    "       'xtick.top': True,\n",
    "       'lines.linewidth': 2,\n",
    "       'axes.linewidth': 2,\n",
    "       'axes.labelweight': 2,\n",
    "       'axes.titleweight': 3,\n",
    "       'ytick.major.width': 2,\n",
    "       'ytick.minor.width': 1,\n",
    "       'ytick.major.size': 6,\n",
    "       'ytick.minor.size': 4,\n",
    "       'ytick.direction': 'in',\n",
    "       'ytick.right': True,\n",
    "       'text.usetex': False,\n",
    "       'figure.figsize': [6, 6],\n",
    "       'figure.facecolor': 'white'\n",
    "       }\n",
    "\n",
    "    plt.rcParams.update(params)\n",
    "\n",
    "    fig = plt.figure()\n",
    "    fig.subplots_adjust(hspace=0, wspace=0)\n",
    "\n",
    "    for i in range(len(tobs_list)):\n",
    "        col = (np.random.random(), np.random.random(), np.random.random())\n",
    "\n",
    "        ax1 = plt.subplot(221)\n",
    "        ax1.plot(x_list[i], z_list[i], '.', color=col, ms=5)\n",
    "        ax1.set_ylabel('Z (au)')\n",
    "        ax1.set_xlabel('X (au)')\n",
    "        if title:\n",
    "            ax1.set_title(title)\n",
    "        ax1.minorticks_on()\n",
    "\n",
    "        ax2 = plt.subplot(222, sharey=ax1)\n",
    "        ax2.plot(y_list[i], z_list[i], '.', color=col, ms=5)\n",
    "        ax2.set_xlabel('Y (au)')\n",
    "        ax2.minorticks_on()\n",
    "        plt.setp(ax2.get_yticklabels(), visible=False)\n",
    "\n",
    "        ax3 = plt.subplot(223, sharex=ax1)\n",
    "        ax3.plot(x_list[i], y_list[i], '.', color=col, ms=5)\n",
    "        ax3.set_ylabel('Y (au)')\n",
    "        ax3.set_xlabel('X (au)')\n",
    "        ax3.minorticks_on()\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d95cd805-a4e1-4e28-9469-83768b904572",
   "metadata": {},
   "source": [
    "Set some plotting defaults to make plots look nice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05b20297-76ac-4d89-9d59-7ddbc666e281",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.style.use('tableau-colorblind10')\n",
    "%matplotlib inline\n",
    "\n",
    "params = {'axes.labelsize': 20,\n",
    "          'font.size': 20,\n",
    "          'legend.fontsize': 14,\n",
    "          'xtick.major.width': 3,\n",
    "          'xtick.minor.width': 2,\n",
    "          'xtick.major.size': 8,\n",
    "          'xtick.minor.size': 4,\n",
    "          'xtick.direction': 'in',\n",
    "          'xtick.top': True,\n",
    "          'lines.linewidth': 3,\n",
    "          'axes.linewidth': 3,\n",
    "          'axes.labelweight': 3,\n",
    "          'axes.titleweight': 3,\n",
    "          'ytick.major.width': 3,\n",
    "          'ytick.minor.width': 2,\n",
    "          'ytick.major.size': 8,\n",
    "          'ytick.minor.size': 4,\n",
    "          'ytick.direction': 'in',\n",
    "          'ytick.right': True,\n",
    "          'figure.figsize': [6, 6],\n",
    "          'figure.facecolor': 'White'\n",
    "          }\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bc6e346-f760-4752-94d7-43cdec77d710",
   "metadata": {},
   "source": [
    "## 2. Orbital parameters of Solar System objects"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "610b4ab6-72b9-44c6-b79b-20b3dd2780f8",
   "metadata": {},
   "source": [
    "Query the DP0.3 catalogs, joining the MPCORB and SSObject tables on their ssObjectId. One _could_ retrieve the entire catalog (containing millions of rows), but that query may take a long time. To minimize query time but still extract a large sample, we use a range of `ssObjectId` values."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfcbf7c6-47c1-41e6-8688-3ec2349eb239",
   "metadata": {},
   "source": [
    "Define the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3af91287-91ea-4756-a7ed-0e561dc4a4b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"SELECT mpc.ssObjectId, mpc.e, mpc.incl, mpc.q, mpc.peri, \n",
    "        sso.ssObjectId, sso.g_H, sso.r_H, sso.i_H, sso.z_H \n",
    "        FROM dp03_catalogs_10yr.MPCORB as mpc \n",
    "        JOIN dp03_catalogs_10yr.SSObject as sso \n",
    "        ON mpc.ssObjectId = sso.ssObjectId \n",
    "        WHERE mpc.ssObjectId < 9223370875126069107 \n",
    "        AND mpc.ssObjectId > 7331137166374808576 \n",
    "        AND sso.numObs > 50\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71ce24d3-9b97-4712-baab-752cacb75f61",
   "metadata": {},
   "source": [
    "Submit the asynchronous query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92f8e4dd-f794-47f4-81e0-cc59e120bf26",
   "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",
    "results = job.fetch_result()\n",
    "assert job.phase == 'COMPLETED'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6ee9f5a-f73b-4d87-94d2-ee29cab81e14",
   "metadata": {},
   "source": [
    "Store the results as a pandas dataframe called `df_mpc_sso` and display number of results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d975c7b9-8c74-40b9-b516-54153ead7e7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mpc_sso = results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fdab0aa-e789-4f78-a7ec-9c2375a2c4a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The query returned ', len(df_mpc_sso), ' results.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "730055a4-0356-4d34-892f-f72e622edb3f",
   "metadata": {},
   "source": [
    "### 2.1 Overall orbital properties\n",
    "\n",
    "First, calculate the semi-major axis of each object's orbit, using the function we defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0051165a-4721-4adb-b946-6039f454474b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "a = calc_semimajor_axis(df_mpc_sso.q, df_mpc_sso.e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1852f4e0-190a-4210-8803-c96353343ac6",
   "metadata": {},
   "source": [
    "Plot the orbital inclination and the eccentricity vs. semi-major axis. Note that we plot the sine of the inclination so its range is limited to values between 0 and 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7242a621-b35b-4c2d-b22f-3d6bbf62ca6a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2, 1, figsize=(9, 5), sharex=True)\n",
    "fig.subplots_adjust(hspace=0.1)\n",
    "\n",
    "im0 = axs[0].hexbin(a, np.sin(np.deg2rad(df_mpc_sso.incl)),\n",
    "              gridsize=(int(103/0.2), 50), cmap='Spectral_r', bins='log',\n",
    "              extent=(0, 103, 0, 1))\n",
    "axs[0].set_ylabel('sin(incl)')\n",
    "axs[0].minorticks_on()\n",
    "axs[0].set_xlim(0, 103)\n",
    "\n",
    "im1 = axs[1].hexbin(a, df_mpc_sso.e,\n",
    "              gridsize=(int(103/0.2), 50), cmap='Spectral_r', bins='log',\n",
    "              extent=(0, 103, 0, 1))\n",
    "axs[1].set_ylabel('eccentricity')\n",
    "axs[1].set_xlabel('semi-major axis (au)')\n",
    "\n",
    "axs[1].set_xlim(0, 103)\n",
    "axs[1].set_ylim(0, 1.09)\n",
    "axs[1].minorticks_on()\n",
    "plt.colorbar(im0)\n",
    "plt.colorbar(im1)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10338e96-95bd-4642-afc8-28cf3ddfcbbf",
   "metadata": {},
   "source": [
    "> **Figure 1:** Orbital inclination (top) and eccentricity (bottom) as a function of semi-major axis for the retrieved subset of solar system objects. Points represent the number of objects in the two-dimensional bin, as illustrated by the colorbars at right.\n",
    " \n",
    "This plot shows Solar System objects extending out to at least 100 au, but the majority of objects are in the inner Solar System (at less than 6 au). For the rest of this notebook, we will focus on the inner Solar System -- specifically on the so-called \"Main Belt\" asteroids. (The outer Solar System bodies will be the subject of another notebook.)\n",
    "\n",
    "### 2.2 Main Belt asteroids\n",
    "\n",
    "#### 2.2.1 Semi-major axes of Main Belt asteroids\n",
    "\n",
    "This plot shows Solar System objects extending out to at least 100 au, but the majority of objects shown are at less than 10 au. For the rest of this notebook, we will focus on the so-called [\"Main Belt\" asteroids](https://en.wikipedia.org/wiki/Asteroid_belt), which roughly speaking inhabit the space between the orbits of Mars and Jupiter. The exact Main Belt definition is not uniform in the literature; for DP0.3 we adopt the definition used by the [JPL small-body database](https://ssd.jpl.nasa.gov/tools/sbdb_query.htmlhttps://ssd.jpl.nasa.gov/tools/sbdb_query.html), where the definitions of the orbital families are visible as the popup notes associated with the \"orbital class\" info. This definition places the Main Belt between 1.6 au < a < 4.6 au, with the constraint q > 1.6 au to ensure that objects don't cross Mars' orbit. However, in this notebook, we alter the definition a bit and consider objects with 1.6 au < a < 4.2 au, q > 1.6, to allow comparison to figures from other papers.\n",
    "\n",
    "Before we get into the orbital elements of the Main Belt population, let's just see where they are located in the Solar System (via their `semi-major axis (a)`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8db5c14c-9f70-4552-aeef-389296b44bbf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "main_belt = (a > 1.6) & (a < 4.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb910281-8289-4c87-83eb-c24647c99e32",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 4))\n",
    "\n",
    "plt.hist(a[main_belt], bins=np.arange(1.6, 4.3, 0.01), color='Black')\n",
    "plt.xlabel('semi-major axis (au)')\n",
    "plt.ylabel('number of objects')\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f700aede-c2e1-46c2-ade5-5e343cc5a141",
   "metadata": {},
   "source": [
    "> **Figure 2:** A histogram of the number of Main Belt asteroids as a function of semi-major axis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8fe8faa-82bf-46d3-8c6b-bd956c732e4e",
   "metadata": {},
   "source": [
    "Clearly it's not just a uniform \"belt\" between Mars and Jupiter -- it has some \"gaps\" where there are very few (or no) objects. These are known as [Kirkwood gaps](https://en.wikipedia.org/wiki/Kirkwood_gap), which arise due to resonances between the asteroid's orbital period and that of Jupiter. For example, at a=2.5 au, an asteroid orbits 3 times for each single orbit of Jupiter, and is thus in a \"3:1 resonance\". Here's a figure illustrating the positions of some of the resonances, and the names of some of the prominent families ([Chris Mihos's ASTR 221 webpage](http://burro.cwru.edu/Academics/Astr221/Gravity/resonance.html)).\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/rubin-dp0/tutorial-notebooks/main/data/kirkwood_resonances.gif\" width=\"350\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ac9b552-2142-4c23-a59c-64b15bdcc0c7",
   "metadata": {},
   "source": [
    "> **Figure 3:** The number of asteroids as a function of semi-major axis. Jupiter resonances are marked with vertical dotted lines, and four families are labeled. This plot is not based on DP0.3 and represents only real objects."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7568e63a-6ef1-4b6e-834d-e09ed628469f",
   "metadata": {},
   "source": [
    "Compare Figure 3 with Figure 2 and notice that the Hildas, Cybeles, and Hungarias are clearly visible in the DP0.3 dataset.\n",
    "\n",
    "We can also add the locations of the mean-motion resonances present in the main belt to our histogram above. We use Kepler's Third Law, (P/P_planet)^2 = (a/a_planet)^3, with the ratio of the orbital periods of the asteroid (`P`) and planet (`P_planet`) along with the semi-major axis of the planet (`a_planet`) to caclulate the semi-major axis at which the corresponding mean-motion resonance occurs: a_res = (P/P_planet)^(2/3) * a_planet. For example, the 3:1 mean-motion resonance with Jupiter, which corresponds to an asteroid orbital period that is 1/3 as long as Jupiter's orbital period, occurs at a_res = ((1/3)^0.66667 * (5.2 au) = 2.5 au."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dc3040b-a408-41fa-9d0c-5e4070a6a439",
   "metadata": {},
   "outputs": [],
   "source": [
    "a_res_J_3_1 = calc_mean_motion_res_a((1/3), 5.2)\n",
    "a_res_J_5_2 = calc_mean_motion_res_a((2/5), 5.2)\n",
    "a_res_J_7_3 = calc_mean_motion_res_a((3/7), 5.2)\n",
    "a_res_J_2_1 = calc_mean_motion_res_a((1/2), 5.2)\n",
    "a_res_J_4_1 = calc_mean_motion_res_a((1/4), 5.2)\n",
    "a_res_J_3_2 = calc_mean_motion_res_a((2/3), 5.2)\n",
    "\n",
    "fig = plt.figure(figsize=(6, 4))\n",
    "\n",
    "plt.hist(a[main_belt], bins=np.arange(1.6, 4.3, 0.01), color='Black')\n",
    "plt.axvline(x=a_res_J_3_1, color='black', linestyle='dashed', linewidth=0.5)\n",
    "plt.axvline(x=a_res_J_5_2, color='black', linestyle='dashed', linewidth=0.5)\n",
    "plt.axvline(x=a_res_J_7_3, color='black', linestyle='dashed', linewidth=0.5)\n",
    "plt.axvline(x=a_res_J_2_1, color='black', linestyle='dashed', linewidth=0.5)\n",
    "plt.axvline(x=a_res_J_4_1, color='black', linestyle='dashed', linewidth=0.5)\n",
    "plt.axvline(x=a_res_J_3_2, color='black', linestyle='dashed', linewidth=0.5)\n",
    "plt.annotate('J 3:1', xy=(a_res_J_3_1-0.08, 5000.), xytext=(a_res_J_3_1-0.08, 4800.), rotation=90.0, fontsize=10)\n",
    "plt.annotate('J 5:2', xy=(a_res_J_5_2-0.08, 5000.), xytext=(a_res_J_5_2-0.08, 4800.), rotation=90.0, fontsize=10)\n",
    "plt.annotate('J 7:3', xy=(a_res_J_7_3-0.08, 5000.), xytext=(a_res_J_7_3-0.08, 4800.), rotation=90.0, fontsize=10)\n",
    "plt.annotate('J 2:1', xy=(a_res_J_2_1-0.08, 5000.), xytext=(a_res_J_2_1-0.08, 4800.), rotation=90.0, fontsize=10)\n",
    "plt.annotate('J 4:1', xy=(a_res_J_4_1-0.08, 5000.), xytext=(a_res_J_4_1-0.08, 4800.), rotation=90.0, fontsize=10)\n",
    "plt.annotate('J 3:2', xy=(a_res_J_3_2-0.08, 5000.), xytext=(a_res_J_3_2-0.08, 4800.), rotation=90.0, fontsize=10)\n",
    "plt.xlabel('semi-major axis (au)')\n",
    "plt.ylabel('number of objects')\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "307c2b3e-2a7e-463c-a465-ba8f2a55aa71",
   "metadata": {},
   "source": [
    "> **Figure 4:** Same as Figure 2, but with Jupiter resonance lines added like in Figure 3.\n",
    "\n",
    "Now we look at the main-belt asteroid's closest approach to the Sun in their orbits, known as the `perihelion distance (q)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3e495c4-5478-41d1-8f9d-cf8358b352d0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 4))\n",
    "\n",
    "plt.hist(df_mpc_sso[main_belt].q, bins=np.arange(0.0, 4.3, 0.02),\n",
    "         color='Black', log=True, histtype='step')\n",
    "plt.xlabel('perihelion distance (au)')\n",
    "plt.ylabel('number of objects')\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7526ee4-a28a-4b93-b13f-18d24c1f1474",
   "metadata": {
    "tags": []
   },
   "source": [
    "> **Figure 5:** The log number of objects as a function of perihelion distance.\n",
    "\n",
    "Hmmm... Should we be concerned that some of those cross the Earth's orbit at 1 au(!)? That's left as an exercise for the learner to explore. (Note that Main Belt asteroids are by definition not going to approach the Earth, as they are typically defined to have perihelion distances beyond the orbit of Mars.)\n",
    "\n",
    "#### 2.2.2 Orbital properties of Main Belt asteroids\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/rubin-dp0/tutorial-notebooks/main/data/Kirkwood-20060509.png\" width=\"250\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed3e44dc-a91e-414a-a671-ccae6299b854",
   "metadata": {},
   "source": [
    "> **Figure 6:** The above illustration (from the Wikipedia article on [Kirkwood gaps](https://en.wikipedia.org/wiki/Kirkwood_gap)) shows planets as large orange circles, and asteroids as small dots of various colors. The Main Belt is clear as a large ring outside Mars' orbit (and inside Jupiter's), and you can clearly see the Kirkwood gaps in the Main Belt.\n",
    "\n",
    "But we can go even further by classifying objects based on their orbital properties."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3219b071-1e8d-4f68-973e-96b5b4f8cb49",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/rubin-dp0/tutorial-notebooks/main/data/ps1_ss_paper_fig5_mainbelt.jpeg\" width=\"350\">\n",
    "<!-- <img src=\"data/ps1_ss_paper_fig5_mainbelt.jpeg\" width=\"350\"> -->"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ac4d387-6db8-491a-9d1c-426337ec41d6",
   "metadata": {},
   "source": [
    "> **Figure 7:** The plot above is from the [Grav et al. 2011 paper](https://ui.adsabs.harvard.edu/abs/2011PASP..123..423G/abstract) describing the PanSTARRS Synthetic Solar System Model, which (in part) forms the basis of DP0.3's simulated dataset. This figure shows the orbital inclinations and eccentricities of \"Main Belt asteroids\" (according to that paper's definition of the Main Belt) as a function of their semi-major axis.\n",
    "\n",
    "Let's make a similar plot for DP0.3 data and see how it compares to the earlier simulation. (First, we'll make the cut to remove objects with perihelia < 1.6 au.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c4b91ec-5d68-4235-b4d3-6149bc909457",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "not_near_mars_peri = df_mpc_sso.q > 1.6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f843d3af-9728-44d3-a595-771725f627e2",
   "metadata": {},
   "source": [
    "Now plot the properties of these selected objects. Note that because DP0.3 (and, ultimately, LSST) has far more detected asteroids than previous SDSS-based studies, we will need to plot density of points instead of a scatter plot. Thus, in each of the following plots, the color scale represents the number of points, from purple/blue as the fewest to red regions containing the most objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7506a13a-e236-4f51-8052-29a6230b8558",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2, 1, figsize=(9, 5), sharex=True)\n",
    "fig.subplots_adjust(hspace=0.1)\n",
    "\n",
    "im0 = axs[0].hexbin(a[not_near_mars_peri],\n",
    "                    np.sin(np.deg2rad(df_mpc_sso[not_near_mars_peri].incl)),\n",
    "                    gridsize=(int(4.7/0.01), 200),\n",
    "                    cmap='Spectral_r', bins='log',\n",
    "                    extent=(0, 4.7, 0, 1))\n",
    "axs[0].set_ylabel('sin(incl)')\n",
    "axs[0].minorticks_on()\n",
    "axs[0].set_xlim(0, 103)\n",
    "\n",
    "im1 = axs[1].hexbin(a[not_near_mars_peri],\n",
    "                    df_mpc_sso[not_near_mars_peri].e,\n",
    "                    gridsize=(int(4.7/0.01), 200),\n",
    "                    cmap='Spectral_r', bins='log',\n",
    "                    extent=(0, 4.7, 0, 1))\n",
    "axs[1].set_ylabel('eccentricity')\n",
    "axs[1].set_xlabel('semi-major axis (au)')\n",
    "\n",
    "axs[1].set_xlim(1.2, 4.7)\n",
    "axs[1].set_ylim(0, 1.09)\n",
    "axs[1].minorticks_on()\n",
    "plt.colorbar(im0)\n",
    "plt.colorbar(im1)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87334ba6-2b52-42b9-a769-3c48bb74fedf",
   "metadata": {},
   "source": [
    "> **Figure 8:** Similar to Figure 1, but for Main Belt asteroids only."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e7c0142-8773-4ec2-a5a8-0f8e2643ca22",
   "metadata": {},
   "source": [
    "## 3. Orbit families among Main Belt asteroids\n",
    "\n",
    "Besides the Kirkwood gaps, the distribution of asteroids in the above plot clearly shows concentrations (\"blobs\") of objects that share similar orbital properties. Let's examine those in more detail.\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/rubin-dp0/tutorial-notebooks/main/data/parker2008_orbit_families_sini_vs_a.png\" width=300>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7ba41dc-96f8-4669-860e-617d5d4a0aa0",
   "metadata": {},
   "source": [
    "> **Figure 9** Inclination versus semi-major axis of Main Belt asteroids from SDSS, from [Parker et al. 2008](https://ui.adsabs.harvard.edu/abs/2008Icar..198..138P/abstract). (See also [Ivezic et al. 2002](https://ui.adsabs.harvard.edu/abs/2002AJ....124.2943I/abstract) and Chapter 5 of the [LSST Science Book](https://ui.adsabs.harvard.edu/abs/2009arXiv0912.0201L/abstract).)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60856ba3-3247-45fd-8395-38d46eb1552b",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/rubin-dp0/tutorial-notebooks/main/data/parker2008_orbit_families_sini_vs_e.png\" width=600>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b3ffddf-8138-4c30-a8fe-9a2a17c02987",
   "metadata": {},
   "source": [
    "> **Figure 10:** Inclination versus eccentricity for Main Belt asteroids in the inner- (left), mid- (center), and outer-belt (right) from SDSS (Parker et al. 2008). The resolution of this plot is low but the x-axis is \"e\" and y-axis is \"sin(i)\".\n",
    "\n",
    "Below, divide the Main Belt asteroids into the following subsets:\n",
    "- inner belt: 2.0 < a < 2.5 au\n",
    "- mid belt: 2.5 < a < 2.82 au\n",
    "- outer belt: 2.82 < a < 3.25 au"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10b5af75-74c5-440a-8446-5075951fce5f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "inner_belt = (a > 2.0) & (a < 2.5) & not_near_mars_peri\n",
    "mid_belt = (a >= 2.5) & (a < 2.82) & not_near_mars_peri\n",
    "outer_belt = (a > 2.82) & (a < 3.25) & not_near_mars_peri"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec5589b3-ef19-4071-9535-129f896e1e7e",
   "metadata": {
    "tags": []
   },
   "source": [
    "Plot `sin(incl)` (inclination) vs. `eccentricity (e)` to compare to Figure 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "198f74a9-0969-44b1-8112-18950738292c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)\n",
    "\n",
    "axs[0].hexbin(df_mpc_sso[inner_belt].e, np.sin(np.deg2rad(df_mpc_sso[inner_belt].incl)),\n",
    "              gridsize=300, bins='log', cmap='Spectral_r')\n",
    "axs[0].set_xlabel('e')\n",
    "axs[0].set_ylabel('sin(i)')\n",
    "axs[0].set_xlim(0, 1)\n",
    "axs[0].set_ylim(0, 1)\n",
    "axs[0].minorticks_on()\n",
    "axs[0].set_title(r'inner belt (2.0 < a < 2.5)', fontsize=16)\n",
    "\n",
    "axs[1].hexbin(df_mpc_sso[mid_belt].e, np.sin(np.deg2rad(df_mpc_sso[mid_belt].incl)),\n",
    "              gridsize=300, bins='log', cmap='Spectral_r')\n",
    "axs[1].set_xlabel('e')\n",
    "axs[1].set_xlim(0, 1)\n",
    "axs[1].set_ylim(0, 1)\n",
    "axs[1].minorticks_on()\n",
    "axs[1].set_title(r'mid belt (2.5 < a < 2.82)', fontsize=16)\n",
    "\n",
    "axs[2].hexbin(df_mpc_sso[outer_belt].e, np.sin(np.deg2rad(df_mpc_sso[outer_belt].incl)),\n",
    "              gridsize=300, bins='log', cmap='Spectral_r')\n",
    "axs[2].set_xlabel('e')\n",
    "axs[2].set_xlim(0, 1)\n",
    "axs[2].set_ylim(0, 1)\n",
    "axs[2].minorticks_on()\n",
    "axs[2].set_title(r'outer belt (2.82 < a < 3.25)', fontsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccf5efa4-af09-4cdb-9596-86fa360073ba",
   "metadata": {
    "tags": []
   },
   "source": [
    "> **Figure 11:** Inclination versus eccentricity for Main Belt asteroids in the inner- (left), mid- (center), and outer-belt (right) from DP0.3. Compare with Figure 10 from SDSS. As in Figure 1 and 8, objects are binned and the colors purple to red represent lower to higher numbers of objects.\n",
    "\n",
    "As you could see in the plot from the Parker et al. paper above, most of the \"blobs\" are at low inclinations and low eccentricities. Zoom in on the lower-left portion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0c80b79-b9d5-460a-94ac-a012e1870122",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)\n",
    "\n",
    "axs[0].hexbin(df_mpc_sso[inner_belt].e, np.sin(np.deg2rad(df_mpc_sso[inner_belt].incl)),\n",
    "              gridsize=200, bins='log', cmap='Spectral_r')\n",
    "axs[0].set_xlabel('e')\n",
    "axs[0].set_ylabel('sin(i)')\n",
    "axs[0].set_xlim(0, 0.275)\n",
    "axs[0].set_ylim(0, 0.275)\n",
    "axs[0].minorticks_on()\n",
    "\n",
    "axs[1].hexbin(df_mpc_sso[mid_belt].e, np.sin(np.deg2rad(df_mpc_sso[mid_belt].incl)),\n",
    "              gridsize=200, bins='log', cmap='Spectral_r')\n",
    "axs[1].set_xlabel('e')\n",
    "axs[1].set_xlim(0, 0.275)\n",
    "axs[1].set_ylim(0, 0.275)\n",
    "axs[1].minorticks_on()\n",
    "\n",
    "axs[2].hexbin(df_mpc_sso[outer_belt].e, np.sin(np.deg2rad(df_mpc_sso[outer_belt].incl)),\n",
    "              gridsize=200, bins='log', cmap='Spectral_r')\n",
    "axs[2].set_xlabel('e')\n",
    "axs[2].set_xlim(0, 0.275)\n",
    "axs[2].set_ylim(0, 0.275)\n",
    "axs[2].minorticks_on()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d28a589-cfac-4f1e-809f-8324d3a0576c",
   "metadata": {},
   "source": [
    "> **Figure 12:** A zoom in on the lower right corners of Figure 11, to better match the axes limits of Figure 10, for comparison."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c1734a8-9b7c-41fc-882c-50906a9f0754",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/rubin-dp0/tutorial-notebooks/main/data/parker2008_orbit_families_sini_vs_e.png\" width=600>\n",
    "\n",
    "> **Figure 13:** A repeat of Figure 10, to enable direct visual comparison with Figure 12 above.\n",
    "\n",
    "Above, see that many of the same features are present. One could look into those features in more detail (indeed, many of them are _named_ orbital families), but that is beyond the scope of this notebook.\n",
    "\n",
    "Finally, to give some intuition about what the actual _orbits_ of these orbital families look like, we will plot them in heliocentric (Sun-centered) Cartesian (XYZ) coordinates. We'll plot 10 randomly-selected objects from \"blobs\" (i.e., orbital families) in the inner, mid, and outer belts in three separate plots.\n",
    "\n",
    "Create a plot overlaying the boxes to select objects from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d4b7fc0-5344-41f0-8bac-c5c990f929c8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)\n",
    "\n",
    "axs[0].hexbin(df_mpc_sso[inner_belt].e, np.sin(np.deg2rad(df_mpc_sso[inner_belt].incl)),\n",
    "              gridsize=200, bins='log', cmap='Spectral_r')\n",
    "axs[0].plot([0.1, 0.1, 0.15, 0.15, 0.1], [0.09, 0.13, 0.13, 0.09, 0.09], color='Black')\n",
    "axs[0].set_xlabel('e')\n",
    "axs[0].set_ylabel('sin(i)')\n",
    "axs[0].set_xlim(0, 0.275)\n",
    "axs[0].set_ylim(0, 0.275)\n",
    "axs[0].minorticks_on()\n",
    "\n",
    "axs[1].hexbin(df_mpc_sso[mid_belt].e, np.sin(np.deg2rad(df_mpc_sso[mid_belt].incl)),\n",
    "              gridsize=200, bins='log', cmap='Spectral_r')\n",
    "axs[1].plot([0.15, 0.15, 0.2, 0.2, 0.15], [0.21, 0.25, 0.25, 0.21, 0.21], color='Black')\n",
    "axs[1].set_xlabel('e')\n",
    "axs[1].set_xlim(0, 0.275)\n",
    "axs[1].set_ylim(0, 0.275)\n",
    "axs[1].minorticks_on()\n",
    "\n",
    "axs[2].hexbin(df_mpc_sso[outer_belt].e, np.sin(np.deg2rad(df_mpc_sso[outer_belt].incl)),\n",
    "              gridsize=200, bins='log', cmap='Spectral_r')\n",
    "axs[2].plot([0.04, 0.04, 0.1, 0.1, 0.04], [0.15, 0.2, 0.2, 0.15, 0.15], color='Black')\n",
    "axs[2].set_xlabel('e')\n",
    "axs[2].set_xlim(0, 0.275)\n",
    "axs[2].set_ylim(0, 0.275)\n",
    "axs[2].minorticks_on()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97ad7c60-ab2a-4136-b10f-7633b6cb5b51",
   "metadata": {},
   "source": [
    "> **Figure 14:** Same as Figure 12, but with open black squares marking the high-density locations to be used below.\n",
    "\n",
    "The following three cells below pick 10 objects at random, then use the `get_obj_matched` function defined above to select all measurements of each object from the SSSource table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cabfd1d0-970b-4357-ac13-6bd3ddd633d1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "inner_blob = (df_mpc_sso.e > 0.1) & (df_mpc_sso.e < 0.15) & \\\n",
    "             (np.sin(np.deg2rad(df_mpc_sso.incl)) < 0.13) & \\\n",
    "             (np.sin(np.deg2rad(df_mpc_sso.incl)) > 0.09)\n",
    "\n",
    "nobj_to_plot = 10\n",
    "pick_obj_inner = df_mpc_sso[inner_belt & inner_blob].sample(nobj_to_plot)\n",
    "\n",
    "pick_obj_torb_inner = []\n",
    "pick_obj_x_inner = []\n",
    "pick_obj_y_inner = []\n",
    "pick_obj_z_inner = []\n",
    "\n",
    "for i in range(nobj_to_plot):\n",
    "    src_cat_obj = get_obj_matched(pick_obj_inner.ssObjectId.values[i])\n",
    "    pick_obj_torb_inner.append(src_cat_obj['midPointMjdTai'])\n",
    "    pick_obj_x_inner.append(src_cat_obj['heliocentricX'])\n",
    "    pick_obj_y_inner.append(src_cat_obj['heliocentricY'])\n",
    "    pick_obj_z_inner.append(src_cat_obj['heliocentricZ'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77f4c567-3644-4dbc-a347-18c1574e4fdf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mid_blob = (df_mpc_sso.e > 0.15) & (df_mpc_sso.e < 0.2) & \\\n",
    "           (np.sin(np.deg2rad(df_mpc_sso.incl)) < 0.25) & \\\n",
    "           (np.sin(np.deg2rad(df_mpc_sso.incl)) > 0.21)\n",
    "\n",
    "nobj_to_plot = 10\n",
    "pick_obj_mid = df_mpc_sso[mid_belt & mid_blob].sample(nobj_to_plot)\n",
    "\n",
    "pick_obj_torb_mid = []\n",
    "pick_obj_x_mid = []\n",
    "pick_obj_y_mid = []\n",
    "pick_obj_z_mid = []\n",
    "\n",
    "for i in range(nobj_to_plot):\n",
    "    src_cat_obj = get_obj_matched(pick_obj_mid.ssObjectId.values[i])\n",
    "    pick_obj_torb_mid.append(src_cat_obj['midPointMjdTai'])\n",
    "    pick_obj_x_mid.append(src_cat_obj['heliocentricX'])\n",
    "    pick_obj_y_mid.append(src_cat_obj['heliocentricY'])\n",
    "    pick_obj_z_mid.append(src_cat_obj['heliocentricZ'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef8ad893-b0bd-47b5-9709-5f11621d193f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "outer_blob = (df_mpc_sso.e > 0.04) & (df_mpc_sso.e < 0.1) & \\\n",
    "             (np.sin(np.deg2rad(df_mpc_sso.incl)) < 0.2) & \\\n",
    "             (np.sin(np.deg2rad(df_mpc_sso.incl)) > 0.15)\n",
    "\n",
    "nobj_to_plot = 10\n",
    "pick_obj_outer = df_mpc_sso[outer_belt & outer_blob].sample(nobj_to_plot)\n",
    "\n",
    "pick_obj_torb_outer = []\n",
    "pick_obj_x_outer = []\n",
    "pick_obj_y_outer = []\n",
    "pick_obj_z_outer = []\n",
    "\n",
    "for i in range(nobj_to_plot):\n",
    "    src_cat_obj = get_obj_matched(pick_obj_outer.ssObjectId.values[i])\n",
    "    pick_obj_torb_outer.append(src_cat_obj['midPointMjdTai'])\n",
    "    pick_obj_x_outer.append(src_cat_obj['heliocentricX'])\n",
    "    pick_obj_y_outer.append(src_cat_obj['heliocentricY'])\n",
    "    pick_obj_z_outer.append(src_cat_obj['heliocentricZ'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e53780ff-7d7e-4fda-86b3-9ee3bda84cbd",
   "metadata": {},
   "source": [
    "The next three cells use the orbit plotting function we defined above to plot them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a98e0328-a585-4129-8a11-72f85cb251a4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xyz_orbit_plot_manyobjects(pick_obj_torb_inner, pick_obj_x_inner,\n",
    "                           pick_obj_y_inner, pick_obj_z_inner, title='inner')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a13634c0-6303-4727-8efa-8f95cc793cbe",
   "metadata": {},
   "source": [
    "> **Figure 15:** Projections of the detections of 10 inner-belt asteroids in the ZX (top left), ZY (top right), and YX (bottom) planes of the solar system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ecd617e-6085-466d-b2b4-18445efbe7ec",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xyz_orbit_plot_manyobjects(pick_obj_torb_mid, pick_obj_x_mid,\n",
    "                           pick_obj_y_mid, pick_obj_z_mid, title='mid')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed4efe0f-fe5b-4b94-b689-fecefcd7c8bf",
   "metadata": {},
   "source": [
    "> **Figure 16:** Similar to Figure 15, but for mid-belt asteroids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bac41262-7685-4611-80d1-c86284d372a2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xyz_orbit_plot_manyobjects(pick_obj_torb_outer, pick_obj_x_outer,\n",
    "                           pick_obj_y_outer, pick_obj_z_outer, title='outer')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3abae71d-1883-43d5-8971-2577a82a8047",
   "metadata": {},
   "source": [
    "> **Figure 17:** Similar to Figure 15, but for outer-belt asteroids.\n",
    "\n",
    "In the above three plots there are small differences between the orbits of these families, but nothing dramatic. (Though that may not be surprising given how similar the inclinations and eccentricities of the different families we selected were in the first place.) More exploration would be required to really understand and classify the different families.\n",
    "\n",
    "## 4. Exercises for the learner\n",
    "\n",
    "1. Explore the properties of other families in the Main Belt. Do their orbits look different than the ones we've examined so far?\n",
    "\n",
    "2. Figure out whether the objects seen in Section 2.2.1 with perihelion distance less than 1.6 au (i.e., ones that travel inside the orbit of Mars) make a close approach to the Earth or Mars.\n",
    "\n",
    "3. Compare properties of Main Belt asteroids with those elsewhere in the Solar System."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4cb05da-1244-404b-b2a7-0af65db2b22c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}