04b. Advanced Phase Curve Modeling#

DP03_04b_Advanced_Phase_Curve_Modeling

Phase Curve of Solar System Objects¶

Contact authors: Yumi Choi and Christina Williams
Last verified to run: 2024-12-17
LSST Science Piplines version: Weekly 2024_50
Container Size: medium
Targeted learning level: advanced

Description: Explicitly investigate the derivation of phase curves for Main Belt asteroids.

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.

LSST Data Products: TAP tables dp03_catalogs_10yr.SSObject, dp03_catalogs_10yr.MPCORB, dp03_catalogs_10yr.DiaSource, dp03_catalogs_10yr.SSSource

Packages: lsst.rsp.get_tap_service

Credit: 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: Muinonen et al. (2010) and David Young's webpage. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Find DP0 documentation and resources at dp0.lsst.io. Questions are welcome as new topics in the Support - Data Preview 0 Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

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.

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:

$$H(α)=H−2.5log_{10}[(1−G)Φ_1(α)+GΦ_2(α)],$$

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.

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). The HG1G2_model has the form

$$H(α)=H−2.5log_{10}[G_1Φ_1(α)+G_2Φ_2(α) + (1-G_1-G_2)Φ3(α)],$$

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}$,

for $G_{12}$ > 0.2, $$G_1 = 0.9529\times G_{12} + 0.02162,$$ $$G_2 = -0.6125\times G_{12} + 0.5572,$$ for $G_{12}$ < 0.2, $$G_1 = 0.7527\times G_{12} + 0.06164,$$ $$G_2 = -0.9612\times G_{12} + 0.6270.$$

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.

1.1 Package Imports¶

The matplotlib (and especially sublibrary matplotlib.pyplot), numpy, and scipy libraries are widely used Python libraries for plotting and scientific computing, and model fitting.

The lsst.rsp package provides access to the Table Access Protocol (TAP) service for queries to the DP0 catalogs.

The pandas package enables table and dataset manipulation.

The sbpy package is an Astropy affiliated package for small-body planetary astronomy.

In [1]:

import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt
import pandas as pd
from sbpy.photometry import HG, HG12, HG1G2
from sbpy.data import Obs

from lsst.rsp import get_tap_service

1.2 Define Functions and Parameters¶

1.2.1 Set up some plotting defaults¶

In [2]:

plt.style.use('tableau-colorblind10')
params = {'axes.labelsize': 15,
          'font.size': 15,
          'legend.fontsize': 12}
plt.rcParams.update(params)

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 this github repository.

In [3]:

filts = ['g', 'r', 'i', 'z']
filter_colors = {'g': '#008060', 'r': '#ff4000', 'i': '#850000', 'z': '#6600cc'}

1.2.2 Define functions for phase curve fitting¶

The following cells define functions to perform the model fitting using the phase curve evalution functions available from the sbpy package.

In [4]:

def weighted_dev(params, mag, phase, mag_err, model):
    """
    Compute weighted deviation for a given model.
    
    Parameters
    ----------
    params: list
        phase curve parameters
    mag: ndarray
        reduced magnitude
    phase: ndarray
        phase angle in radians
    mag_err: ndarray
        uncertainty in magnitude
    model: function
        phase curve model function
    

    Returns
    -------
    sol: tuple
        best-fit solution
    """

    if model is HG1G2.evaluate:
        pred = model(phase, params[0], params[1], params[2])
    else:
        pred = model(phase, params[0], params[1])

    return (mag - pred)/mag_err

In [5]:

def fitPhaseCurve(mag, phase, sigma, model=HG12.evaluate, params=[0.1]):
    """
    Fit phase curve for given observations to a designated model.

    Parameters
    ----------
    mag: ndarray
        reduced magnitude
    phase: ndarray
        phase angle in degrees
    sigma: ndarray
        uncertainty in magnitude
    model: function (default=HG12.evaluate)
        phase curve model function
    params: list (default=[0.1])
        phase curve parameters

    Returns
    -------
    sol: tuple
        best-fit solution
    """

    phase = np.deg2rad(phase)
    sol = leastsq(weighted_dev, [mag[0]] + params, (mag, phase, sigma, model),
                  full_output=True)

    return sol

In [6]:

def fitAllPhaseCurveModels(reducedMag, magSigma, phaseAngle, verbose=False):
    """
    Fit phase curves for given observations to three different models,
    the HG (sol_HG), HG12 (sol_HG12) and HG1G2 (sol_HG1G2), and 
    store the resulting solutions in a dictionary of dictionaries. 
    Save np.nan values when the fit is not converged.

    Parameters
    ----------
    reducedMag: ndarray
        reduced magnitude
    magSigma: ndarray
        uncertainty in magnitude
    phaseAngle: ndarray
        phase angle in degrees

    Returns
    -------
    solutions: dict
        Best-fit solutions for each model
    """

    solutions = {}

    sol_HG = fitPhaseCurve(reducedMag, phaseAngle, magSigma, model=HG.evaluate)

    solutions['HG'] = {}
    try:
        solutions['HG']['chi2'] = np.sum(sol_HG[2]['fvec']**2)
        solutions['HG']['H'] = sol_HG[0][0]
        solutions['HG']['G'] = sol_HG[0][1]
        solutions['HG']['H_err'] = np.sqrt(sol_HG[1][0, 0])
        solutions['HG']['G_err'] = np.sqrt(sol_HG[1][1, 1])
        solutions['HG']['cov'] = sol_HG[1]
    except TypeError:
        if verbose:
            print('HG model is not converging')
        solutions['HG']['chi2'] = np.nan
        solutions['HG']['H'] = np.nan
        solutions['HG']['G'] = np.nan
        solutions['HG']['H_err'] = np.nan
        solutions['HG']['G_err'] = np.nan
        solutions['HG']['cov'] = np.nan

    sol_HG12 = fitPhaseCurve(reducedMag, phaseAngle,
                             magSigma, model=HG12.evaluate)

    solutions['HG12'] = {}
    try:
        solutions['HG12']['chi2'] = np.sum(sol_HG12[2]['fvec']**2)
        solutions['HG12']['H'] = sol_HG12[0][0]
        solutions['HG12']['G12'] = sol_HG12[0][1]
        solutions['HG12']['H_err'] = np.sqrt(sol_HG12[1][0, 0])
        solutions['HG12']['G12_err'] = np.sqrt(sol_HG12[1][1, 1])
        solutions['HG12']['cov'] = sol_HG12[1]
    except TypeError:
        if verbose:
            print('HG12 model is not converging')
        solutions['HG12']['chi2'] = np.nan
        solutions['HG12']['H'] = np.nan
        solutions['HG12']['G12'] = np.nan
        solutions['HG12']['H_err'] = np.nan
        solutions['HG12']['G12_err'] = np.nan
        solutions['HG12']['cov'] = np.nan

    sol_HG1G2 = fitPhaseCurve(reducedMag, phaseAngle, magSigma,
                              model=HG1G2.evaluate, params=[0.1, 0.1])

    solutions['HG1G2'] = {}
    try:
        solutions['HG1G2']['chi2'] = np.sum(sol_HG1G2[2]['fvec']**2)
        solutions['HG1G2']['H'] = sol_HG1G2[0][0]
        solutions['HG1G2']['G1'] = sol_HG1G2[0][1]
        solutions['HG1G2']['G2'] = sol_HG1G2[0][2]
        solutions['HG1G2']['H_err'] = np.sqrt(sol_HG1G2[1][0, 0])
        solutions['HG1G2']['G1_err'] = np.sqrt(sol_HG1G2[1][1, 1])
        solutions['HG1G2']['G2_err'] = np.sqrt(sol_HG1G2[1][2, 2])
        solutions['HG1G2']['cov'] = sol_HG1G2[1]
    except TypeError:
        if verbose:
            print('HG1G2 model is not converging')
        solutions['HG1G2']['chi2'] = np.nan
        solutions['HG1G2']['H'] = np.nan
        solutions['HG1G2']['G1'] = np.nan
        solutions['HG1G2']['G2'] = np.nan
        solutions['HG1G2']['H_err'] = np.nan
        solutions['HG1G2']['G1_err'] = np.nan
        solutions['HG1G2']['G2_err'] = np.nan
        solutions['HG1G2']['cov'] = np.nan

    return solutions

2. Query the DP0.3 tables for Main Belt asteroids¶

2.1 Create the Rubin TAP Service Client¶

Get an instance of the TAP service, and assert that it exists.

In [7]:

service = get_tap_service("ssotap")
assert service is not None

2.2 Query the DP0.3 SSObject and MPCORB catalogs¶

The phase curve modeling requires apparent magnitudes & uncertainties, phase angles, topocentric ($d_t$) and heliocentric ($d_h$) distances.

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.

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 DP0.3 documentation.

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).

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.

In [8]:

nobs_thrh = '2000'
min_a = '2.0'
max_a = '3.2'
min_q = '1.666'

Define query.

In [9]:

query = """
SELECT
    mpc.ssObjectId, mpc.mpcG, mpc.mpcH,
    mpc.q/(1-mpc.e) as a,
    sso.arc, sso.numObs,
    sso.g_H, sso.g_Herr, sso.g_G12, sso.g_G12err, sso.g_H_gG12_Cov,
    sso.r_H, sso.r_Herr, sso.r_G12, sso.r_G12err, sso.r_H_rG12_Cov,
    sso.i_H, sso.i_Herr, sso.i_G12, sso.i_G12err, sso.i_H_iG12_Cov,
    sso.z_H, sso.z_Herr, sso.z_G12, sso.z_G12err, sso.z_H_zG12_Cov
FROM
    dp03_catalogs_10yr.MPCORB as mpc
INNER JOIN dp03_catalogs_10yr.SSObject as sso
ON mpc.ssObjectId = sso.ssObjectId
WHERE sso.numObs > {} AND mpc.q/(1-mpc.e) > {} AND
mpc.q/(1-mpc.e) < {} AND mpc.q > {} ORDER by sso.ssObjectId
""".format(nobs_thrh, min_a, max_a, min_q)

Submit asynchronous query, and fetch the results into astropy table uniqueObj.

In [10]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
uniqueObj = job.fetch_result().to_table()
assert job.phase == 'COMPLETED'

Job phase is COMPLETED

Show the results.

In [11]:

uniqueObj

Out[11]:

Table length=296

ssObjectId	mpcG	mpcH	a	arc	numObs	g_H	g_Herr	g_G12	g_G12err	g_H_gG12_Cov	r_H	r_Herr	r_G12	r_G12err	r_H_rG12_Cov	i_H	i_Herr	i_G12	i_G12err	i_H_iG12_Cov	z_H	z_Herr	z_G12	z_G12err	z_H_zG12_Cov
		mag		d		mag	mag	mag	mag	mag2	mag	mag	mag	mag	mag2	mag	mag	mag	mag	mag2	mag	mag	mag	mag	mag2
int64	float32	float32	float64	float32	int32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32
-9217466392671047318	0.15	14.48	2.7534679055396287	3099.2915	2866	15.203103	0.00038897144	0.4550666	0.0040478734	7.9422597e-07	14.554703	0.00021039916	0.48015854	0.0020262992	1.7115104e-07	14.354623	0.00025000537	0.4709413	0.0025779733	2.3757198e-07	14.407989	0.0006071515	0.51356506	0.005866336	1.1365114e-06
-9186168674478550259	0.15	15.54	2.735597883342819	2920.819	2064	16.26421	0.0009085277	0.5034483	0.009313577	3.7164054e-06	15.6143055	0.00045167241	0.5149425	0.0049309437	1.8081512e-07	15.414917	0.0005772725	0.50797635	0.00510145	7.0264207e-07	15.471048	0.0011641541	0.5433281	0.009459758	1.9743543e-06
-9023703058184288186	0.15	16.64	2.940021025169084	2817.7273	2331	17.358274	0.0035483497	0.50554883	0.036873408	7.245376e-05	16.716164	0.001927695	0.5128336	0.018558158	2.0433228e-05	16.515278	0.0023644213	0.4923302	0.024016123	3.2025448e-05	16.568754	0.005117289	0.5263797	0.04942868	0.00014278179
-8991757915344909776	0.15	18.32	2.712164954223674	2911.9167	2027	19.041605	0.0039131506	0.49194965	0.038261477	2.652303e-05	18.39314	0.0022384885	0.46744	0.020893808	1.3363821e-05	18.197392	0.0027081787	0.5020755	0.026521428	1.2971898e-05	18.232428	0.00799786	0.31482488	0.08765045	0.0002543192
-8983147311557262481	0.15	16.81	3.1285476112080075	3223.7615	2103	17.358887	0.0022244945	0.5284655	0.028648542	1.4641261e-05	16.87871	0.0015134424	0.5128175	0.018769464	6.977196e-06	16.761572	0.001905264	0.5324352	0.025683051	7.7371615e-06	16.753658	0.0048452453	0.46961442	0.05186777	0.00012877512
-8941112935633645120	0.15	17.46	2.4713281790115262	3111.7393	2357	18.188131	0.0013846699	0.5111664	0.011265048	8.2808214e-07	17.539248	0.0007647462	0.48283827	0.005792018	3.7366075e-07	17.337635	0.0010335406	0.5097	0.007862881	-1.0053393e-06	17.386505	0.0027831802	0.5375992	0.025589064	-2.3151732e-05
-8926263449476792198	0.15	16.98	2.603187433175253	3454.771	2145	17.701952	0.0016780668	0.49290684	0.017006207	-1.5563393e-06	17.05407	0.0009794187	0.46493313	0.00934684	4.142698e-07	16.85306	0.0011744339	0.48502782	0.011868707	-9.0575486e-07	16.907093	0.0032809514	0.50106376	0.031326078	3.8512946e-05
-8898442558843979238	0.15	18.23	2.540751939505157	3007.7913	2022	18.774021	0.0023720781	0.5111807	0.017291997	1.17132095e-05	18.30697	0.0015417456	0.4712598	0.010613613	5.945477e-06	18.186533	0.0020921014	0.49740168	0.015085877	7.081695e-06	18.17325	0.004575927	0.5023939	0.031217633	3.5207522e-05
-8848841739843330269	0.15	16.99	3.043081576855482	3627.043	2617	17.53464	0.0016498883	0.4948764	0.01436597	8.768694e-06	17.061077	0.0011163031	0.48569593	0.0095817065	4.305333e-06	16.940798	0.0014123438	0.48687568	0.012630372	6.103225e-06	16.943336	0.003720296	0.5370097	0.036412716	5.3106625e-05
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
8810278553610239375	0.15	17.35	2.798872785591477	3647.0981	2452	17.901987	0.0013853513	0.46128002	0.014875833	5.8035316e-06	17.427883	0.0009080462	0.4850499	0.009180252	2.6425416e-06	17.306454	0.0012100416	0.4284831	0.012724048	3.7698287e-06	17.300947	0.0033608018	0.43586305	0.02993238	5.3218144e-05
8819108235107391621	0.15	15.01	2.6588872083083026	3468.7954	2858	15.557955	0.00031843685	0.48927864	0.002990388	8.0105494e-08	15.084722	0.00022092048	0.48191044	0.0018312577	1.08507855e-07	14.965879	0.00026620654	0.48024583	0.0022703246	1.0680802e-07	14.96012	0.0006686684	0.47443017	0.0062590926	1.9770528e-06
8880916212552915142	0.15	17.0	3.148844765324078	3610.8494	2238	17.721386	0.0037739042	0.49607134	0.032425635	7.96902e-05	17.069847	0.0021867214	0.47974047	0.019748664	2.867727e-05	16.871984	0.0026700455	0.5010986	0.024540653	4.0883348e-05	16.93313	0.0073610554	0.59678143	0.067584604	0.00029311952
8884270739944543284	0.15	13.65	2.7849987623877737	3209.164	2001	14.197751	0.00020256465	0.50636727	0.001886042	1.7905312e-07	13.724491	0.00012941653	0.46272156	0.0013378037	6.530088e-08	13.605648	0.00016325305	0.46933052	0.0016966178	9.598327e-08	13.596713	0.0003350527	0.5351535	0.0033027532	2.7697655e-07
8885214464318942411	0.15	15.26	3.0410713294788416	3631.0134	2415	15.806544	0.000694903	0.49274746	0.0064901104	2.139879e-06	15.332101	0.00046983914	0.4851357	0.0043359916	1.0186678e-06	15.213519	0.00062158675	0.46541023	0.0061650826	1.7904817e-06	15.204619	0.0014996033	0.488601	0.013780486	9.833656e-06
9074174377699860353	0.15	17.34	2.9450152076891394	3224.8936	2218	17.885767	0.003241317	0.52326035	0.0316753	2.6578367e-05	17.409258	0.0020771963	0.4952344	0.018664544	1.4278274e-05	17.2946	0.002804425	0.5797351	0.026659677	2.3194569e-05	17.285494	0.008001417	0.6097138	0.06698252	0.00029604958
9087912756022919586	0.15	16.92	2.7357773950448507	3276.0242	2486	17.470213	0.0014550819	0.5047754	0.011454952	7.291216e-06	17.000196	0.00090536	0.49189544	0.008134237	1.0244829e-06	16.880655	0.0012810556	0.5054339	0.011883684	1.3454323e-06	16.874708	0.0029589871	0.5324183	0.023864923	1.5999802e-05
9103245968292935281	0.15	17.53	2.841129420564853	3646.1116	2003	18.072693	0.002341692	0.4326196	0.027041448	1.2289271e-05	17.601795	0.0015449624	0.4737852	0.01584573	7.0289516e-06	17.484308	0.0020029848	0.47202682	0.021001507	8.883437e-06	17.47665	0.0057830913	0.56938684	0.05978415	7.257145e-05
9208707402267266207	0.15	17.06	2.734850357707651	2856.8145	2051	17.7885	0.0017436614	0.46866193	0.014496459	1.5021225e-05	17.138794	0.0009869646	0.47262788	0.008634665	4.886307e-06	16.939726	0.00122158	0.46071374	0.01077141	7.193949e-06	16.997112	0.0034783634	0.46373403	0.024120187	5.944238e-05

2.3 Query the DP0.3 DiaSource and SSSource catalogs¶

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.

Define query.

In [12]:

query = """
SELECT
    dia.ssObjectId, dia.diaSourceId, dia.mag,
    dia.magErr, dia.band, dia.midPointMjdTai,
    sss.phaseAngle, sss.topocentricDist, sss.heliocentricDist
FROM
    dp03_catalogs_10yr.DiaSource as dia
INNER JOIN
    dp03_catalogs_10yr.SSSource as sss
ON
    dia.diaSourceId = sss.diaSourceId
WHERE
    dia.ssObjectId
    IN {}
ORDER by dia.ssObjectId
""".format(tuple(uniqueObj['ssObjectId']))

Submit asynchronous query and fetch the results into indivObs.

In [13]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
indivObs = job.fetch_result().to_table()
assert job.phase == 'COMPLETED'

Job phase is COMPLETED

Show the results.

In [14]:

indivObs

Out[14]:

Table length=668765

ssObjectId	diaSourceId	mag	magErr	band	midPointMjdTai	phaseAngle	topocentricDist	heliocentricDist
					d	deg	AU	AU
int64	int64	float32	float32	str1	float64	float32	float32	float32
-9217466392671047318	5074726492115737	20.055	0.012	r	61021.35247	17.301477	2.711414	3.1537251
-9217466392671047318	19481485662802488	19.669	0.007	i	61033.2511	15.801364	2.5499394	3.1503248
-9217466392671047318	39913165364548592	20.21	0.009	r	61010.34128	18.036163	2.8695178	3.1563327
-9217466392671047318	48327202766987632	19.45	0.01	z	61049.34712	12.501292	2.3604758	3.14477
-9217466392671047318	48371439483286632	19.711	0.008	i	61031.33981	16.094881	2.5749133	3.1509116
-9217466392671047318	62341640397910521	20.195	0.008	g	61051.17905	12.03446	2.3417952	3.144068
-9217466392671047318	73744426738833766	20.075	0.016	i	61005.35388	18.184	2.9424236	3.1573431
-9217466392671047318	75023334125204409	20.044	0.014	r	61021.3483	17.301888	2.711473	3.1537263
-9217466392671047318	79046699514151278	20.346	0.009	r	60998.32104	18.215149	3.045487	3.158588
...	...	...	...	...	...	...	...	...
9208707402267266207	-45612237486212510	22.469	0.054	g	61022.32521	20.811378	2.1577377	2.634653
9208707402267266207	-43978633643892617	21.167	0.013	r	61065.28681	9.683634	1.7987144	2.7099721
9208707402267266207	-43961380036767235	21.737	0.03	r	61029.34984	19.76049	2.0822046	2.6469705
9208707402267266207	-38749719734873518	21.199	0.073	z	61051.21659	14.5168705	1.8825521	2.68535
9208707402267266207	-25615344533198694	21.774	0.075	i	61005.35429	22.218847	2.350334	2.6050045
9208707402267266207	-22473258342931197	21.505	0.066	z	62026.9898	18.178625	1.9840554	2.678289
9208707402267266207	-17229010326108879	21.168	0.014	r	61064.29393	10.050997	1.8031344	2.7082386
9208707402267266207	-12082581769521416	21.032	0.018	i	61061.28922	11.144528	1.8179785	2.702988
9208707402267266207	-8420768209572350	21.804	0.031	r	61029.35151	19.760197	2.0821872	2.6469736

Confirm that the number of unique objects in indivObs is identical to that of uniqueObj, as they should be.

In [15]:

assert len(uniqueObj) == len(np.unique(indivObs['ssObjectId']))

3. Fit phase curves of unique main belt asteroids for three different models¶

3.1 Compute reduced magnitude¶

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:

$$H(α) = m−5log_{10}(d_t\times\,d_h),$$

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.

In [16]:

thdist = indivObs['topocentricDist']*indivObs['heliocentricDist']
reduced_mag = indivObs['mag'] - 5.0*np.log10(thdist)

indivObs.add_column(reduced_mag, name='reducedMag')

3.2 Fit phase curve per unique object per filter¶

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.

In [17]:

fitted_array = []

for iobj in uniqueObj['ssObjectId']:
    idx = indivObs['ssObjectId'] == iobj
    tmp = indivObs[idx]
    filts_tmp = np.unique(tmp['band'])
    for ifilt in filts_tmp:
        idx_filt = tmp['band'] == ifilt
        nobs_ifilt = len(tmp[idx_filt])

        if nobs_ifilt > 3:
            x_fitted = fitAllPhaseCurveModels(tmp['reducedMag'][idx_filt],
                                              tmp['magErr'][idx_filt],
                                              tmp['phaseAngle'][idx_filt])
            fitted_array.append([iobj, ifilt, x_fitted])

results = pd.DataFrame(fitted_array)
results.columns = ['ssObjectId', 'fname', 'fit_param']

Finally, convert the fit parameter dictionary to individual columns in a pandas dataframe to make it easy to read each parameter.

In [18]:

L = ['ssObjectId', 'fname']
results = results[L].join(pd.json_normalize(results.fit_param))

3.3 Plot the three different best-fit phase curves in each filter for a single object¶

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.

You can pick an integer number between 0 and len(uniqueObj)-1 for iObj below to explore other objects.

In [19]:

iObj = 100
sId = uniqueObj['ssObjectId'][iObj]
tmp = indivObs[indivObs['ssObjectId'] == sId]
phases = np.linspace(0, 90, 100)

for i, ifilt in enumerate(filts):
    idx = tmp['band'] == ifilt
    if i==0:
        labels = ['HG 2-parameter model',
                  'HG12 2-parameter model', 
                  'HG1G2 3-parameter model'
                 ]
    else:
        labels = [None, None, None]

    plt.errorbar(tmp['phaseAngle'][idx], tmp['reducedMag'][idx],
                 yerr=tmp['magErr'][idx], fmt='.',
                 color=filter_colors[ifilt], alpha=0.5, label=ifilt)

    HG_mag = HG.evaluate(np.deg2rad(phases),
                         results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG.H'].values,
                         results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG.G'].values)
    plt.plot(phases, HG_mag, color=filter_colors[ifilt],
             label=labels[0])

    HG12_mag = HG12.evaluate(np.deg2rad(phases),
                             results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG12.H'].values,
                             results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG12.G12'].values)
    plt.plot(phases, HG12_mag, color=filter_colors[ifilt], linestyle='--',
             label=labels[1])

    HG1G2_mag = HG1G2.evaluate(np.deg2rad(phases),
                               results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG1G2.H'].values,
                               results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG1G2.G1'].values,
                               results[(results.ssObjectId == sId) & (results.fname == ifilt)]['HG1G2.G2'].values)
    plt.plot(phases, HG1G2_mag, color=filter_colors[ifilt], linestyle='dotted',
             label=labels[2])

plt.xlim(tmp['phaseAngle'].min()-5, tmp['phaseAngle'].max()+5)
plt.ylim(tmp['reducedMag'].max()+0.5, tmp['reducedMag'].min()-0.5)
plt.xlabel('Phase Angle [deg]')
plt.ylabel('Reduced magnitude [mag]')
plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')
plt.title('ssObjectId = %d' % sId)
plt.show()

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.

4. Compare manually-derived phase curve parameters to those provided in DP0.3¶

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.

4.1 Compare to the pipeline-provided phase curve parameters stored in the SSObject table¶

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.

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.

In [20]:

fig = plt.figure(figsize=(10, 7))
gs = fig.add_gridspec(2, 2, wspace=0, hspace=0)
axs = gs.subplots(sharex=True, sharey=True)
axs = axs.ravel()

one2one = np.arange(0.01, 1, .01)
filts = ['g', 'r', 'i', 'z']
for i, ifilt in enumerate(filts):

    axs[i].errorbar(results[results.fname == ifilt]['HG12.G12'],
                    uniqueObj[ifilt+'_G12'],
                    xerr=results[results.fname == ifilt]['HG12.G12_err'],
                    yerr=uniqueObj[ifilt+'_G12err'], fmt='.', alpha=0.5,
                    label='G12 HG12 model')

    axs[i].plot(one2one, one2one, '--', label='1:1 correspondence')
    axs[i].text(0.75, 0.1, ifilt+'-band')

fig.supxlabel('Slope parameter G [HG12_model]')
fig.supylabel('Slope parameter G [SSObject Table]')
axs[0].legend(loc=2)
plt.show()

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.

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.

In [21]:

fig = plt.figure(figsize=(10, 7))
gs = fig.add_gridspec(2, 2, wspace=0, hspace=0)
axs = gs.subplots(sharex=True, sharey=True)
axs = axs.ravel()

one2one = np.arange(11, 22, 1)
for i, ifilt in enumerate(filts):

    axs[i].errorbar(results[results.fname == ifilt]['HG12.H'],
                    uniqueObj[ifilt+'_H'],
                    xerr=results[results.fname == ifilt]['HG12.H_err'],
                    yerr=uniqueObj[ifilt+'_Herr'], fmt='.', alpha=0.5,
                    label='H HG12 model')

    axs[i].plot(one2one, one2one, '--', label='1:1 correspondence')
    axs[i].text(18.5, 12, ifilt+'-band')

fig.supxlabel('Absolute magnitude H [HG12_model]')
fig.supylabel('Absolute magnitude H [SSObject Table]')
axs[0].legend(loc=2)
plt.show()

Figure 3: Similar to Figure 2, but for the absolute H magnitude instead of the slope.

4.2 Compare to the input phase curve parameters stored in the MPCORB table¶

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 Table 1 of the DP0.3 documentation. 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).

In [22]:

fig = plt.figure(figsize=(10, 4))
gs = fig.add_gridspec(1, 2, wspace=0)
axs = gs.subplots(sharey=True)
axs = axs.ravel()

one2one = np.arange(0.01, 1, .01)
for i, ifilt in enumerate(filts):

    h = axs[0].hist(results[results.fname == ifilt]['HG.G']
                    - uniqueObj['mpcG'], bins=50, density=True,
                    histtype='step', color=filter_colors[ifilt], label=ifilt)

    h = axs[1].hist(results[results.fname == ifilt]['HG.H']
                    - uniqueObj['mpcH'], bins=50, density=True,
                    histtype='step', color=filter_colors[ifilt])

axs[0].set_xlabel('Out - In (Slope parameter G)')
axs[0].set_ylabel('N')
axs[0].plot([0, 0], [0, 80], 'k--',label='zero bias')
axs[0].plot([0.02, 0.02], [0, 80], 'k',label='detected bias')
axs[0].set_ylim(0, 70)
axs[0].legend(loc=1)
axs[1].set_xlabel('Out - In (Absolute magnitude H)')
plt.show()

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.

5. An incomplete list of outstanding questions about DP0.3 phase curves¶

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:

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).
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.
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.

6. Excercises for the learner¶

Compare the 3 fits for objects with less phase coverage (e.g., Jupiter Trojans).
Fit phase curves for TNOs with a linear function.
Apply the filter conversion in each Rubin filter and revisit the second plot in Section 4.2.