DP03_04b_Advanced_Phase_Curve_Modeling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "# Phase Curve of Solar System Objects\n",
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "Contact authors: Yumi Choi and Christina Williams<br>\n",
    "Last verified to run: 2024-12-17 <br>\n",
    "LSST Science Piplines version: Weekly 2024_50 <br>\n",
    "Container Size: medium <br>\n",
    "Targeted learning level: advanced <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da1a210-d858-42fe-8591-570965b8be1a",
   "metadata": {},
   "source": [
    "**Description:** Explicitly investigate the derivation of phase curves for Main Belt asteroids."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0baf5-51ad-40ec-8991-060a7b27c289",
   "metadata": {},
   "source": [
    "**Skills:** Use the TAP service and ADQL to access the DP0.3 tables. Join information from multiple DP0.3 tables. Derive phase curves using three different models."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393da88f-7978-4920-aa4a-a9830df6eed9",
   "metadata": {},
   "source": [
    "**LSST Data Products:** TAP tables dp03_catalogs_10yr.SSObject, dp03_catalogs_10yr.MPCORB, dp03_catalogs_10yr.DiaSource, dp03_catalogs_10yr.SSSource"
   ]
  },
  {
   "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:**\n",
    "Inspired by a jupyter notebook developed by Queen's University Belfast Planet Lab (including Brian Rogers, Niall McElroy, and Meg Schwamb). Standalone functions for phase curve fitting were developed by Pedro Bernardinelli. References: <a href=\"https://ui.adsabs.harvard.edu/abs/2010Icar..209..542M/abstract\">Muinonen et al. (2010)</a> and <a href=\"http://astronotes.co.uk/blog/2018/05/28/determining-the-h-g-parameters-of-atlas-asteroid-phase-curves.html\">David Young's webpage.</a> Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a400d1bb-e1ec-4160-bd41-bc7ca953a97c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-27T17:55:39.523317Z",
     "iopub.status.busy": "2023-07-27T17:55:39.522902Z",
     "iopub.status.idle": "2023-07-27T17:55:40.035942Z",
     "shell.execute_reply": "2023-07-27T17:55:40.034503Z",
     "shell.execute_reply.started": "2023-07-27T17:55:39.523290Z"
    },
    "tags": []
   },
   "source": [
    "**Get Support:** Find DP0 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 notebook targets for an advanced user who is interested in understanding phase curve fitting in DP0.3 in detail. This notebook will explicitly explore three different functional forms that are broadly used for modeling a phase curve of an asteroid, and make a comparison with the `H` and `G12` parameters stored in the `SSObject` table. Please refer to the Introduction of the companion tutorial notebook 04a for definition of the phase curve. Only the most relevent parts to this tutorial are repeated here.\n",
    "\n",
    "Modeling the phase curve (reduced magnitude $H(α)$ as a function of phase angle $α$) enables measurement of the absolute magnitude, $H$, defined as the brightness at 0 phase angle. The functional form can depend on the data quality and type of object targeted. This tutorial will explore three functional forms that are relevant for understanding the DP0.3 data products: `HG_model`, `HG1G2_model`, and `HG12_model`. The `HG_model` is the simplest model, and has the form:\n",
    "\n",
    "$$H(α)=H−2.5log_{10}[(1−G)Φ_1(α)+GΦ_2(α)],$$\n",
    "\n",
    "where $Φ_n$ are basis functions (which allow for one to model non-linearity in the data, but maintain linearity in the fitted parameters). $Φ_n$ are normalized to unity at α=0 deg. This model provides a best fit for the slope parameter $G$ (from which surface properties can then be derived) and the absolute magnitude $H$. $H$(α) is the reduced magnitude at a given phase angle α if measured at 1 au away from both Earth and from the Sun (i.e. unit topocentric and heliocentric distance). The `HG_model` $G$ and $H$ values are stored in the `dp03_catalogs_10yr.MPCORB` table. For further info on the `HG_model`, see [Bowell et al. 1989](https://ui.adsabs.harvard.edu/abs/1989aste.conf..524B/abstract). \n",
    "\n",
    "To better accommodate various observational effects (e.g., photometric quality, incomplete phase angle sampling) the more sophisticated `HG1G2_model` (a linear three-parameter function) and its nonlinear two-parameter version `HG12_model` were developed (see [Muinonen et al. 2010](https://ui.adsabs.harvard.edu/abs/2010Icar..209..542M/abstract)). The `HG1G2_model` has the form\n",
    "\n",
    "$$H(α)=H−2.5log_{10}[G_1Φ_1(α)+G_2Φ_2(α) + (1-G_1-G_2)Φ3(α)],$$\n",
    "\n",
    "which now has three free parameters, $H$, $G_1$ and $G_2$. However, a third representation, the `HG12_model`, is generally very effective for deriving reliable values of absolute magnitude when the phase angle sampling is not optimal (e.g., poor phase angle coverage at a range of phase angle). Thus, the LSST data products will compute estimated parameters of the `HG12_model` and this will be the focus of this tutorial. The `HG12_model` expresses the $G_1$ and $G_2$ parameters as a piecewise linear function of a single parameter, $G_{12}$, \n",
    "\n",
    "for $G_{12}$ > 0.2,\n",
    "$$G_1 = 0.9529\\times G_{12} + 0.02162,$$\n",
    "$$G_2 = -0.6125\\times G_{12} + 0.5572,$$\n",
    "for $G_{12}$ < 0.2,\n",
    "$$G_1 = 0.7527\\times G_{12} + 0.06164,$$\n",
    "$$G_2 = -0.9612\\times G_{12} + 0.6270.$$\n",
    "\n",
    "This notebook presents an illustrative example of deriving the phase curve of Main Belt asteroids using DP0.3 simulated catalogs. In Section 2, we will retrieve the necessary data from the DP0.3 database. Section 3 will delve into manual phase curve fitting for a single object using three different parametrizations of the phase curve model. In Section 4, we will then compare this manual fitting to the automated phase curve fitting carried out as part of the LSST data products, which are available in the SSObject table."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bf722b4-8227-441f-9f3e-0f8140040fc1",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1.1 Package Imports\n",
    "\n",
    "The [matplotlib](https://matplotlib.org/) (and especially sublibrary `matplotlib.pyplot`), [numpy](http://www.numpy.org/), and [scipy](https://scipy.org/) libraries are widely used Python libraries for plotting and scientific computing, and model fitting.\n",
    "\n",
    "The `lsst.rsp` package provides access to the Table Access Protocol (TAP) service for queries to the DP0 catalogs.\n",
    "\n",
    "The [pandas](https://pandas.pydata.org/) package enables table and dataset manipulation.\n",
    "\n",
    "The [sbpy](https://sbpy.org) package is an `Astropy` affiliated package for small-body planetary astronomy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cddc1458",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import leastsq\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from sbpy.photometry import HG, HG12, HG1G2\n",
    "from sbpy.data import Obs\n",
    "\n",
    "from lsst.rsp import get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c217adff-25ed-4fce-95e7-8aa04630f6cc",
   "metadata": {},
   "source": [
    "### 1.2 Define Functions and Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a4a8b3b-91cf-4729-a707-99b7b11a5b8d",
   "metadata": {},
   "source": [
    "#### 1.2.1 Set up some plotting defaults"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c780a382-5c5e-4f20-9ffd-7aa517c37146",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.style.use('tableau-colorblind10')\n",
    "params = {'axes.labelsize': 15,\n",
    "          'font.size': 15,\n",
    "          'legend.fontsize': 12}\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c8563f1-2309-4486-a884-de242859b088",
   "metadata": {},
   "source": [
    "Set up colors and plot symbols corresponding to the $g,r,i,z$ bands since there are no $u$ and $y$ band data in the DP0.3 catalogs. These colors are the same as those used for $g,r,i,z$ bands in Dark Energy Survey (DES) publications, and are defined in <a href=\"https://github.com/DarkEnergySurvey/descolors\">this github repository</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9187ddd4-d43b-41a5-922e-f33fee1a4454",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "filts = ['g', 'r', 'i', 'z']\n",
    "filter_colors = {'g': '#008060', 'r': '#ff4000', 'i': '#850000', 'z': '#6600cc'}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4d55d16-f2c0-495c-8f43-0b93bc44ed89",
   "metadata": {},
   "source": [
    "#### 1.2.2 Define functions for phase curve fitting "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf5b442f-9c2e-4717-a1bc-6897a5d40992",
   "metadata": {},
   "source": [
    "The following cells define functions to perform the model fitting using the phase curve evalution functions available from the `sbpy` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "513bfe10-d340-435a-ae0a-d6d669c56bee",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def weighted_dev(params, mag, phase, mag_err, model):\n",
    "    \"\"\"\n",
    "    Compute weighted deviation for a given model.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    params: list\n",
    "        phase curve parameters\n",
    "    mag: ndarray\n",
    "        reduced magnitude\n",
    "    phase: ndarray\n",
    "        phase angle in radians\n",
    "    mag_err: ndarray\n",
    "        uncertainty in magnitude\n",
    "    model: function\n",
    "        phase curve model function\n",
    "    \n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    sol: tuple\n",
    "        best-fit solution\n",
    "    \"\"\"\n",
    "\n",
    "    if model is HG1G2.evaluate:\n",
    "        pred = model(phase, params[0], params[1], params[2])\n",
    "    else:\n",
    "        pred = model(phase, params[0], params[1])\n",
    "\n",
    "    return (mag - pred)/mag_err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37d86162-8557-4f65-9813-ec46818466a9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def fitPhaseCurve(mag, phase, sigma, model=HG12.evaluate, params=[0.1]):\n",
    "    \"\"\"\n",
    "    Fit phase curve for given observations to a designated model.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    mag: ndarray\n",
    "        reduced magnitude\n",
    "    phase: ndarray\n",
    "        phase angle in degrees\n",
    "    sigma: ndarray\n",
    "        uncertainty in magnitude\n",
    "    model: function (default=HG12.evaluate)\n",
    "        phase curve model function\n",
    "    params: list (default=[0.1])\n",
    "        phase curve parameters\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    sol: tuple\n",
    "        best-fit solution\n",
    "    \"\"\"\n",
    "\n",
    "    phase = np.deg2rad(phase)\n",
    "    sol = leastsq(weighted_dev, [mag[0]] + params, (mag, phase, sigma, model),\n",
    "                  full_output=True)\n",
    "\n",
    "    return sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cdf1dd4-10b7-4c1f-8f4f-b4ab252822a6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def fitAllPhaseCurveModels(reducedMag, magSigma, phaseAngle, verbose=False):\n",
    "    \"\"\"\n",
    "    Fit phase curves for given observations to three different models,\n",
    "    the HG (sol_HG), HG12 (sol_HG12) and HG1G2 (sol_HG1G2), and \n",
    "    store the resulting solutions in a dictionary of dictionaries. \n",
    "    Save np.nan values when the fit is not converged.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    reducedMag: ndarray\n",
    "        reduced magnitude\n",
    "    magSigma: ndarray\n",
    "        uncertainty in magnitude\n",
    "    phaseAngle: ndarray\n",
    "        phase angle in degrees\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    solutions: dict\n",
    "        Best-fit solutions for each model\n",
    "    \"\"\"\n",
    "\n",
    "    solutions = {}\n",
    "\n",
    "    sol_HG = fitPhaseCurve(reducedMag, phaseAngle, magSigma, model=HG.evaluate)\n",
    "\n",
    "    solutions['HG'] = {}\n",
    "    try:\n",
    "        solutions['HG']['chi2'] = np.sum(sol_HG[2]['fvec']**2)\n",
    "        solutions['HG']['H'] = sol_HG[0][0]\n",
    "        solutions['HG']['G'] = sol_HG[0][1]\n",
    "        solutions['HG']['H_err'] = np.sqrt(sol_HG[1][0, 0])\n",
    "        solutions['HG']['G_err'] = np.sqrt(sol_HG[1][1, 1])\n",
    "        solutions['HG']['cov'] = sol_HG[1]\n",
    "    except TypeError:\n",
    "        if verbose:\n",
    "            print('HG model is not converging')\n",
    "        solutions['HG']['chi2'] = np.nan\n",
    "        solutions['HG']['H'] = np.nan\n",
    "        solutions['HG']['G'] = np.nan\n",
    "        solutions['HG']['H_err'] = np.nan\n",
    "        solutions['HG']['G_err'] = np.nan\n",
    "        solutions['HG']['cov'] = np.nan\n",
    "\n",
    "    sol_HG12 = fitPhaseCurve(reducedMag, phaseAngle,\n",
    "                             magSigma, model=HG12.evaluate)\n",
    "\n",
    "    solutions['HG12'] = {}\n",
    "    try:\n",
    "        solutions['HG12']['chi2'] = np.sum(sol_HG12[2]['fvec']**2)\n",
    "        solutions['HG12']['H'] = sol_HG12[0][0]\n",
    "        solutions['HG12']['G12'] = sol_HG12[0][1]\n",
    "        solutions['HG12']['H_err'] = np.sqrt(sol_HG12[1][0, 0])\n",
    "        solutions['HG12']['G12_err'] = np.sqrt(sol_HG12[1][1, 1])\n",
    "        solutions['HG12']['cov'] = sol_HG12[1]\n",
    "    except TypeError:\n",
    "        if verbose:\n",
    "            print('HG12 model is not converging')\n",
    "        solutions['HG12']['chi2'] = np.nan\n",
    "        solutions['HG12']['H'] = np.nan\n",
    "        solutions['HG12']['G12'] = np.nan\n",
    "        solutions['HG12']['H_err'] = np.nan\n",
    "        solutions['HG12']['G12_err'] = np.nan\n",
    "        solutions['HG12']['cov'] = np.nan\n",
    "\n",
    "    sol_HG1G2 = fitPhaseCurve(reducedMag, phaseAngle, magSigma,\n",
    "                              model=HG1G2.evaluate, params=[0.1, 0.1])\n",
    "\n",
    "    solutions['HG1G2'] = {}\n",
    "    try:\n",
    "        solutions['HG1G2']['chi2'] = np.sum(sol_HG1G2[2]['fvec']**2)\n",
    "        solutions['HG1G2']['H'] = sol_HG1G2[0][0]\n",
    "        solutions['HG1G2']['G1'] = sol_HG1G2[0][1]\n",
    "        solutions['HG1G2']['G2'] = sol_HG1G2[0][2]\n",
    "        solutions['HG1G2']['H_err'] = np.sqrt(sol_HG1G2[1][0, 0])\n",
    "        solutions['HG1G2']['G1_err'] = np.sqrt(sol_HG1G2[1][1, 1])\n",
    "        solutions['HG1G2']['G2_err'] = np.sqrt(sol_HG1G2[1][2, 2])\n",
    "        solutions['HG1G2']['cov'] = sol_HG1G2[1]\n",
    "    except TypeError:\n",
    "        if verbose:\n",
    "            print('HG1G2 model is not converging')\n",
    "        solutions['HG1G2']['chi2'] = np.nan\n",
    "        solutions['HG1G2']['H'] = np.nan\n",
    "        solutions['HG1G2']['G1'] = np.nan\n",
    "        solutions['HG1G2']['G2'] = np.nan\n",
    "        solutions['HG1G2']['H_err'] = np.nan\n",
    "        solutions['HG1G2']['G1_err'] = np.nan\n",
    "        solutions['HG1G2']['G2_err'] = np.nan\n",
    "        solutions['HG1G2']['cov'] = np.nan\n",
    "\n",
    "    return solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec51ac0b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Query the DP0.3 tables for Main Belt asteroids\n",
    "\n",
    "### 2.1 Create the Rubin TAP Service Client\n",
    "\n",
    "Get an instance of the TAP service, and assert that it exists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81b8cd59-1ba3-4eaa-846f-6478ed0c3cf5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"ssotap\")\n",
    "assert service is not None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf968de5-ce59-41c9-88f0-83cd88bdc083",
   "metadata": {},
   "source": [
    "### 2.2 Query the DP0.3 SSObject and MPCORB catalogs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "171e087c-b314-48f7-8f22-6013115f2cb9",
   "metadata": {},
   "source": [
    "The phase curve modeling requires apparent magnitudes & uncertainties, phase angles, topocentric ($d_t$) and heliocentric ($d_h$) distances.\n",
    "\n",
    "To define the properties of solar system objects, the DP0.3 model uses the `HG_model` form of the phase curve to simulate Rubin observations for each object. These \"truth\" values are defined in the MPCORB table as `mpcH` (intrinsic absolute magnitude in $V$ band) and `mpcG` (intrinsic slope). For the DP0.3 simulations, the intrinsic slope, `mpcG`, for all objects was set to a constant value of 0.15.\n",
    "\n",
    "In the SSObject table, the LSST data products also contain the fitted phase curve parameters based on the mock observations using the `HG12_model` (i.e. contain absolute magnitude `H` and slope parameter `G12` in $g,r,i,z$ bands). Note that the value of `G12` slope will differ from `G` owing to the difference in functional form. The explanation for the absence of $u$ and $y$ bands in DP0.3 catalogs can be found in the <a href=\"https://dp0-3.lsst.io/data-products-dp0-3/data-simulation-dp0-3.html\">DP0.3 documentation</a>.\n",
    "\n",
    "To focus on Main Belt asteroids that likely have good phase curve fits for this tutorial, the query is limited to select sources with a large number of observations (in the SSObject table, this is `numObs` > 2000) and with 2.0 au < semi-major axis (`a`) < 3.2 au and the perihelion distance `q` > 1.666 au following the Main Belt asteroid definition by the JPL Horizons small body database query tool (https://ssd.jpl.nasa.gov/tools/sbdb_query.html).  \n",
    "\n",
    "The table returned by this query is called \"unique\" since it contains the IDs of unique Main Belt asteroids (although each object has many individual observations in LSST). The table should contain 296 unique objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f93c0de6-8730-4adc-a0ab-9b795d3dc4e7",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "nobs_thrh = '2000'\n",
    "min_a = '2.0'\n",
    "max_a = '3.2'\n",
    "min_q = '1.666'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0762fed3-cb68-48a8-b0bd-2354fb0262ea",
   "metadata": {},
   "source": [
    "Define query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1c428cb-7c09-4680-be90-c55d7a73c037",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\n",
    "SELECT\n",
    "    mpc.ssObjectId, mpc.mpcG, mpc.mpcH,\n",
    "    mpc.q/(1-mpc.e) as a,\n",
    "    sso.arc, sso.numObs,\n",
    "    sso.g_H, sso.g_Herr, sso.g_G12, sso.g_G12err, sso.g_H_gG12_Cov,\n",
    "    sso.r_H, sso.r_Herr, sso.r_G12, sso.r_G12err, sso.r_H_rG12_Cov,\n",
    "    sso.i_H, sso.i_Herr, sso.i_G12, sso.i_G12err, sso.i_H_iG12_Cov,\n",
    "    sso.z_H, sso.z_Herr, sso.z_G12, sso.z_G12err, sso.z_H_zG12_Cov\n",
    "FROM\n",
    "    dp03_catalogs_10yr.MPCORB as mpc\n",
    "INNER JOIN dp03_catalogs_10yr.SSObject as sso\n",
    "ON mpc.ssObjectId = sso.ssObjectId\n",
    "WHERE sso.numObs > {} AND mpc.q/(1-mpc.e) > {} AND\n",
    "mpc.q/(1-mpc.e) < {} AND mpc.q > {} ORDER by sso.ssObjectId\n",
    "\"\"\".format(nobs_thrh, min_a, max_a, min_q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a940549d-cd7f-46a5-985c-923272a100b9",
   "metadata": {},
   "source": [
    "Submit asynchronous query, and fetch the results into astropy table `uniqueObj`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5ada866-541d-4f09-ba94-d7ce6ca53949",
   "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",
    "uniqueObj = job.fetch_result().to_table()\n",
    "assert job.phase == 'COMPLETED'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5df6d99-c171-47ce-91ee-d3672a13c7ee",
   "metadata": {},
   "source": [
    "Show the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "311a0a97-e22a-451a-9e89-6dbad9738efc",
   "metadata": {},
   "outputs": [],
   "source": [
    "uniqueObj"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "977dd781-f675-4ed4-89ac-07cb6042515e",
   "metadata": {},
   "source": [
    "### 2.3 Query the DP0.3 DiaSource and SSSource catalogs\n",
    "\n",
    "While there are unique solar system objects in the `SSObject` and `MPCORB` tables, these objects will be observed many times over the full LSST survey. Individual observations of each unique object are recorded in the `SSSource` and `DiaSource` tables. Below, we query these tables to obtain all of the individual observed time series data (we call `indivObs`) for the unique objects (`uniqueObj`) selected above. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6df1db2-9e39-4fdf-9564-7807481cd980",
   "metadata": {},
   "source": [
    "Define query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5727ea1b-0d55-40aa-b548-b1c4ad407cb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\n",
    "SELECT\n",
    "    dia.ssObjectId, dia.diaSourceId, dia.mag,\n",
    "    dia.magErr, dia.band, dia.midPointMjdTai,\n",
    "    sss.phaseAngle, sss.topocentricDist, sss.heliocentricDist\n",
    "FROM\n",
    "    dp03_catalogs_10yr.DiaSource as dia\n",
    "INNER JOIN\n",
    "    dp03_catalogs_10yr.SSSource as sss\n",
    "ON\n",
    "    dia.diaSourceId = sss.diaSourceId\n",
    "WHERE\n",
    "    dia.ssObjectId\n",
    "    IN {}\n",
    "ORDER by dia.ssObjectId\n",
    "\"\"\".format(tuple(uniqueObj['ssObjectId']))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d42927b-e676-451d-b654-0080ce9d6033",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T20:37:20.682035Z",
     "iopub.status.busy": "2024-12-12T20:37:20.681766Z",
     "iopub.status.idle": "2024-12-12T20:37:20.686815Z",
     "shell.execute_reply": "2024-12-12T20:37:20.686074Z",
     "shell.execute_reply.started": "2024-12-12T20:37:20.682015Z"
    }
   },
   "source": [
    "Submit asynchronous query and fetch the results into `indivObs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c588f206-9077-4243-88e0-12f2bd3d6a00",
   "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",
    "indivObs = job.fetch_result().to_table()\n",
    "assert job.phase == 'COMPLETED'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee2adccb-e0e9-4652-b70d-f32fd4a0ab11",
   "metadata": {},
   "source": [
    "Show the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17de01eb-261b-4307-9193-340c6933fc05",
   "metadata": {},
   "outputs": [],
   "source": [
    "indivObs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e366f714-8078-40a5-9eee-65639fa293ca",
   "metadata": {},
   "source": [
    "Confirm that the number of unique objects in `indivObs` is identical to that of `uniqueObj`, as they should be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "077c9400-f3fe-422d-bb2b-7ad0f2587e3d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "assert len(uniqueObj) == len(np.unique(indivObs['ssObjectId']))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f57f57e2-a386-46af-a527-73d397658b72",
   "metadata": {},
   "source": [
    "## 3. Fit phase curves of unique main belt asteroids for three different models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5b36990-98e2-42de-a55f-34d71b4d48bb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-08-23T21:26:56.687471Z",
     "iopub.status.busy": "2023-08-23T21:26:56.687082Z",
     "iopub.status.idle": "2023-08-23T21:26:56.690871Z",
     "shell.execute_reply": "2023-08-23T21:26:56.690216Z",
     "shell.execute_reply.started": "2023-08-23T21:26:56.687445Z"
    },
    "tags": []
   },
   "source": [
    "### 3.1 Compute reduced magnitude"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "392ff14e-eca1-4677-8de8-4dae0e6811bc",
   "metadata": {
    "tags": []
   },
   "source": [
    "To plot the phase curve, we first must compute the reduced magnitude $H(\\alpha)$ for each observation, and add it as a column to the `indivObs` table we produced of individual observations. The reduced magnitude $H(\\alpha)$ is the normalized apparent magnitude of an asteroid as if it is observed at 1 au from both the Sun and the Earth as a function of phase angle $\\alpha$, once accounting for the relative distances between the asteroid, and both Sun and Earth:\n",
    "\n",
    "$$H(α) = m−5log_{10}(d_t\\times\\,d_h),$$\n",
    "\n",
    "where $m$ is the apparent magnitude, and $d_t$ and $d_h$ are the topocentric and heliocentric distances of the object at the time of each observation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a387ad4-6057-4631-84fc-da9f99559b26",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "thdist = indivObs['topocentricDist']*indivObs['heliocentricDist']\n",
    "reduced_mag = indivObs['mag'] - 5.0*np.log10(thdist)\n",
    "\n",
    "indivObs.add_column(reduced_mag, name='reducedMag')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e2dbbf1-bd2b-40b0-8981-282ddcf3848c",
   "metadata": {},
   "source": [
    "### 3.2 Fit phase curve per unique object per filter "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98becb52-7bd3-44d8-9b06-f00b9a419467",
   "metadata": {
    "tags": []
   },
   "source": [
    "The cell below now fits the phase curve for individual unique objects in each LSST filter using the three different fitting functions, `HG_model`, `HG1G2_model` and `HG12_model`. To ensure a good fit, the number of observations `nobs_ifilt` must be greater than the number of fit parameters (3). This cell takes ~3-6 min for 296 unique objects with the medium container size. The output x_fitted contains the parameters that are returned for each model fit.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5771b7ed-a488-47b5-9c19-91beae25843a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fitted_array = []\n",
    "\n",
    "for iobj in uniqueObj['ssObjectId']:\n",
    "    idx = indivObs['ssObjectId'] == iobj\n",
    "    tmp = indivObs[idx]\n",
    "    filts_tmp = np.unique(tmp['band'])\n",
    "    for ifilt in filts_tmp:\n",
    "        idx_filt = tmp['band'] == ifilt\n",
    "        nobs_ifilt = len(tmp[idx_filt])\n",
    "\n",
    "        if nobs_ifilt > 3:\n",
    "            x_fitted = fitAllPhaseCurveModels(tmp['reducedMag'][idx_filt],\n",
    "                                              tmp['magErr'][idx_filt],\n",
    "                                              tmp['phaseAngle'][idx_filt])\n",
    "            fitted_array.append([iobj, ifilt, x_fitted])\n",
    "\n",
    "results = pd.DataFrame(fitted_array)\n",
    "results.columns = ['ssObjectId', 'fname', 'fit_param']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a57c8214-2e6b-4355-af8d-284527705faa",
   "metadata": {},
   "source": [
    "Finally, convert the fit parameter dictionary to individual columns in a pandas dataframe to make it easy to read each parameter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc5cbf77-0299-4063-86f1-533bc42fca60",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "L = ['ssObjectId', 'fname']\n",
    "results = results[L].join(pd.json_normalize(results.fit_param))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a086c0e7-801e-4a65-92a5-400bc6506659",
   "metadata": {},
   "source": [
    "### 3.3 Plot the three different best-fit phase curves in each filter for a single object "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b47093f8-45bf-4bf0-9d17-bdc06eeac679",
   "metadata": {},
   "source": [
    "Now, we will plot example phase curves in all available filters (in DP0.3) $g$,$r$,$i$,$z$ for a single object referenced by its ssObjectId, which we call `sId`. (You can explore different objects by changing the `iObj` index to retrieve different sources). First the observations to be modeled are plotted, and then all three forms of phase curve model defined above. The plot demonstrates that the reduced magnitude and phase curve of the source are offset from each other in each filter, reflecting the variation in brightness of asteroids in different filters. Overall, it is clear that the 3 different fitted models all describe the observations well for well-sampled phase curves in DP0.3."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f9a2b64-6b4e-499d-9733-a465b75b3601",
   "metadata": {},
   "source": [
    "You can pick an integer number between 0 and `len(uniqueObj)-1` for `iObj` below to explore other objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3523f9a4-1773-4731-97c8-8d14c1487846",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "iObj = 100\n",
    "sId = uniqueObj['ssObjectId'][iObj]\n",
    "tmp = indivObs[indivObs['ssObjectId'] == sId]\n",
    "phases = np.linspace(0, 90, 100)\n",
    "\n",
    "for i, ifilt in enumerate(filts):\n",
    "    idx = tmp['band'] == ifilt\n",
    "    if i==0:\n",
    "        labels = ['HG 2-parameter model',\n",
    "                  'HG12 2-parameter model', \n",
    "                  'HG1G2 3-parameter model'\n",
    "                 ]\n",
    "    else:\n",
    "        labels = [None, None, None]\n",
    "\n",
    "    plt.errorbar(tmp['phaseAngle'][idx], tmp['reducedMag'][idx],\n",
    "                 yerr=tmp['magErr'][idx], fmt='.',\n",
    "                 color=filter_colors[ifilt], alpha=0.5, label=ifilt)\n",
    "\n",
    "    HG_mag = HG.evaluate(np.deg2rad(phases),\n",
    "                         results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG.H'].values,\n",
    "                         results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG.G'].values)\n",
    "    plt.plot(phases, HG_mag, color=filter_colors[ifilt],\n",
    "             label=labels[0])\n",
    "\n",
    "    HG12_mag = HG12.evaluate(np.deg2rad(phases),\n",
    "                             results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG12.H'].values,\n",
    "                             results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG12.G12'].values)\n",
    "    plt.plot(phases, HG12_mag, color=filter_colors[ifilt], linestyle='--',\n",
    "             label=labels[1])\n",
    "\n",
    "    HG1G2_mag = HG1G2.evaluate(np.deg2rad(phases),\n",
    "                               results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG1G2.H'].values,\n",
    "                               results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG1G2.G1'].values,\n",
    "                               results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG1G2.G2'].values)\n",
    "    plt.plot(phases, HG1G2_mag, color=filter_colors[ifilt], linestyle='dotted',\n",
    "             label=labels[2])\n",
    "\n",
    "plt.xlim(tmp['phaseAngle'].min()-5, tmp['phaseAngle'].max()+5)\n",
    "plt.ylim(tmp['reducedMag'].max()+0.5, tmp['reducedMag'].min()-0.5)\n",
    "plt.xlabel('Phase Angle [deg]')\n",
    "plt.ylabel('Reduced magnitude [mag]')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')\n",
    "plt.title('ssObjectId = %d' % sId)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "548b4a20-988f-482a-a880-3be192a88748",
   "metadata": {},
   "source": [
    "> **Figure 1:** Three different phase curve models fit to the $griz$ phase curves of a single `SSObject`: the HG (solid) and HG12 (dashed) 2-parameter models, and the HG1G2 3-parameter model (dotted) showing different behaviour in the relatively unconstrained low-angle region."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "192f7650-8beb-4215-b4c5-a04e864f9887",
   "metadata": {},
   "source": [
    "## 4. Compare manually-derived phase curve parameters to those provided in DP0.3\n",
    "\n",
    "As mentioned in Section 2.2, the DP0.3 phase curve simulations were based on the `mpcH` and `mpcG` parameters stored in the `MPCORB` table using the `HG_model`. On the other hand, the simulated LSST observations of phase curves were fitted using the `HG12_model` and the resulting `H` and `G12` parameters are stored in the `SSObject` table. This Section compares those phase curve parameters provided in the DP0.3 tables to the manually-driven parameters using the `sbpy` package in this tutorial. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5d95e14-ac7f-42c4-a99d-e56cef7bb584",
   "metadata": {},
   "source": [
    "### 4.1 Compare to the pipeline-provided phase curve parameters stored in the SSObject table\n",
    "\n",
    "While modeling phase curves manually using these three fitting functions demonstrates the process, the `HG12_model` results (the parametrization that is more stable to limited observations across phase angles) will automatically be tabulated as a data product during the course of the survey in the `SSObject` table. This section will compare the manually derived parameters from the last section with those produced in LSST data products."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d32197ca-5b82-45f8-ab84-0108d5eab863",
   "metadata": {},
   "source": [
    "The figure below compares the manually-derived slope parameter `G12` for all the unique objects to the Rubin pipeline-provided values stored in the `SSObject` table in DP0.3. It is clear that overall the manual measurement recovers the DP0.3 value and uncertainty, which is expected because the `HG12` model was used to automatically fit the DP0.3 simulated data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdf0c5f8-4dbc-4e26-8181-196d7443fa3a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 7))\n",
    "gs = fig.add_gridspec(2, 2, wspace=0, hspace=0)\n",
    "axs = gs.subplots(sharex=True, sharey=True)\n",
    "axs = axs.ravel()\n",
    "\n",
    "one2one = np.arange(0.01, 1, .01)\n",
    "filts = ['g', 'r', 'i', 'z']\n",
    "for i, ifilt in enumerate(filts):\n",
    "\n",
    "    axs[i].errorbar(results[results.fname == ifilt]['HG12.G12'],\n",
    "                    uniqueObj[ifilt+'_G12'],\n",
    "                    xerr=results[results.fname == ifilt]['HG12.G12_err'],\n",
    "                    yerr=uniqueObj[ifilt+'_G12err'], fmt='.', alpha=0.5,\n",
    "                    label='G12 HG12 model')\n",
    "\n",
    "    axs[i].plot(one2one, one2one, '--', label='1:1 correspondence')\n",
    "    axs[i].text(0.75, 0.1, ifilt+'-band')\n",
    "\n",
    "fig.supxlabel('Slope parameter G [HG12_model]')\n",
    "fig.supylabel('Slope parameter G [SSObject Table]')\n",
    "axs[0].legend(loc=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2c10048-b820-422c-b490-9b70b16f3d5f",
   "metadata": {},
   "source": [
    "> **Figure 2:** The slope parameter from the `SSObject` table ($y$-axis) versus the slope parameter fit using the HG12 model in this notebook ($x$-axis), in the four filters $griz$ (upper left to bottom right panels). The 1:1 correspondence line is underplotted with a dashed orange line to demonstrate, visually, how well correlated the two slopes are."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20438935-c6a4-49e4-805f-144c9c3731f3",
   "metadata": {},
   "source": [
    "The below figure compares the manually-derived absolute magnitude `H` for all the unique objects to the Rubin pipeline-provided values stored in the `SSObject` table in DP0.3. Like the `G12` parameter, it is clear that overall the manual measurement recovers the DP0.3 value and uncertainty, which is expected because the `HG12` model was used to automatically fit the DP0.3 simulated data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f88d315c-d039-405d-9a61-5cfcd685481d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 7))\n",
    "gs = fig.add_gridspec(2, 2, wspace=0, hspace=0)\n",
    "axs = gs.subplots(sharex=True, sharey=True)\n",
    "axs = axs.ravel()\n",
    "\n",
    "one2one = np.arange(11, 22, 1)\n",
    "for i, ifilt in enumerate(filts):\n",
    "\n",
    "    axs[i].errorbar(results[results.fname == ifilt]['HG12.H'],\n",
    "                    uniqueObj[ifilt+'_H'],\n",
    "                    xerr=results[results.fname == ifilt]['HG12.H_err'],\n",
    "                    yerr=uniqueObj[ifilt+'_Herr'], fmt='.', alpha=0.5,\n",
    "                    label='H HG12 model')\n",
    "\n",
    "    axs[i].plot(one2one, one2one, '--', label='1:1 correspondence')\n",
    "    axs[i].text(18.5, 12, ifilt+'-band')\n",
    "\n",
    "fig.supxlabel('Absolute magnitude H [HG12_model]')\n",
    "fig.supylabel('Absolute magnitude H [SSObject Table]')\n",
    "axs[0].legend(loc=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2822acb2-906b-4e52-a807-e746057bdff6",
   "metadata": {},
   "source": [
    "> **Figure 3:** Similar to Figure 2, but for the absolute H magnitude instead of the slope."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2589b08-468c-4293-a2f5-e7e344903f60",
   "metadata": {},
   "source": [
    "### 4.2 Compare to the input phase curve parameters stored in the MPCORB table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6b53e88-add4-4fcd-87b6-52e892f7c55d",
   "metadata": {},
   "source": [
    "This section compares the input phase curve parameters, `mpcH` and `mpcG`, that were used to simulate asteroids in DP0.3 with those derived using the `HG_model` in this tutorial to see how well the input parameters are recovered in the simulated LSST results at the end of the 10 year survey.  The left panel shows the difference of $\\sim$0.02 between the input slope `G` and the output slope `G` guided by the red solid lines in all four filters. As shown in the right panel, the distributions of out-in absolute magnitude `H` are more complicated. Recall that `mpcH` is reported in $V$ band; the conversion from Rubin filters to $V$ is summarized in <a href=\"https://dp0-3.lsst.io/data-products-dp0-3/data-simulation-dp0-3.html#synthetic-object-populations\">Table 1 of the DP0.3 documentation</a>. The filter conversion can not fully explain the offsets seen in the out-in distributions (the table of filter conversions between V-band and LSST at the link above are never more than 0.4, but the right-hand panel shows g-bad with an offset of >0.4).   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "693b4d12-4691-4d15-9c6c-07f6e12a4d2a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 4))\n",
    "gs = fig.add_gridspec(1, 2, wspace=0)\n",
    "axs = gs.subplots(sharey=True)\n",
    "axs = axs.ravel()\n",
    "\n",
    "one2one = np.arange(0.01, 1, .01)\n",
    "for i, ifilt in enumerate(filts):\n",
    "\n",
    "    h = axs[0].hist(results[results.fname == ifilt]['HG.G']\n",
    "                    - uniqueObj['mpcG'], bins=50, density=True,\n",
    "                    histtype='step', color=filter_colors[ifilt], label=ifilt)\n",
    "\n",
    "    h = axs[1].hist(results[results.fname == ifilt]['HG.H']\n",
    "                    - uniqueObj['mpcH'], bins=50, density=True,\n",
    "                    histtype='step', color=filter_colors[ifilt])\n",
    "\n",
    "axs[0].set_xlabel('Out - In (Slope parameter G)')\n",
    "axs[0].set_ylabel('N')\n",
    "axs[0].plot([0, 0], [0, 80], 'k--',label='zero bias')\n",
    "axs[0].plot([0.02, 0.02], [0, 80], 'k',label='detected bias')\n",
    "axs[0].set_ylim(0, 70)\n",
    "axs[0].legend(loc=1)\n",
    "axs[1].set_xlabel('Out - In (Absolute magnitude H)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91e0c420-1c0e-4b50-a7d3-34cc6726b667",
   "metadata": {},
   "source": [
    "> **Figure 4:** Histograms of the residual between in the fit output and the simulation input parameters for the slope (left) and absolute H magnitude (right) of the phase curves of a sample of `SSObjects` in four filters, $griz$ (legend). The approximate peak of the residual in slope, $\\sim0.02$, is marked with a solid line. To aid visual comparison, 0.0 is marked with a dashed line."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8525bdd-cfe9-41ca-b599-38ffaf62acf9",
   "metadata": {},
   "source": [
    "## 5. An incomplete list of outstanding questions about DP0.3 phase curves \n",
    "This DP0.3 release revealed some unknown features that can be resolved in the future as people use the simulation. Some things that warrant further exploration:\n",
    "\n",
    "- A small bias (roughly 0.02 in all filters) was identified in the `G` parameter of `HG_model` (average measured value - intrinsic value of `G`).\n",
    "\n",
    "- Offsets between the intrinsic absolute magnitude in $V$ band and recovered absolute magnitude in the LSST filters are larger than listed in the filter-conversion table referenced above.\n",
    "\n",
    "- **Update on Jan, 2024**: The current DP0.3 catalogs reported PSF magnitudes without taking the effect of object trailing into account. When correcting for this effect, the two biases seen above disappear. It will be fixed in subsequent releases."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f880c4d0-6ffc-49e9-86f2-95f1b698e279",
   "metadata": {},
   "source": [
    "## 6. Excercises for the learner\n",
    "\n",
    "1. Compare the 3 fits for objects with less phase coverage (e.g., Jupiter Trojans).\n",
    "2. Fit phase curves for TNOs with a linear function.\n",
    "3. Apply the filter conversion in each Rubin filter and revisit the second plot in Section 4.2."
   ]
  }
 ],
 "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
}