04b. Advanced Phase Curve Modeling#
Phase Curve of Solar System Objects
Contact authors: Yumi Choi and Christina Williams
Last verified to run: 2024-01-31
LSST Science Piplines version: Weekly 2024_04
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.3-related documentation and resources at dp0-3.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.
# general python packages
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
# LSST package for TAP queries
from lsst.rsp import get_tap_service
1.2 Define Functions and Parameters¶
1.2.1 Set up some plotting defaults¶
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.
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.
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
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
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
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.
nobs_thrh = '2000'
min_a = '2.0'
max_a = '3.2'
min_q = '1.666'
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)
uniqueObj = service.search(query).to_table()
uniqueObj
| 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 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 8695259029909021189 | 0.15 | 17.31 | 3.052661942178657 | 2818.7297 | 2112 | 18.032787 | 0.004045596 | 0.4757583 | 0.041904327 | 6.717258e-05 | 17.3829 | 0.0023444758 | 0.42496052 | 0.024744578 | 2.4857925e-05 | 17.182878 | 0.0028847314 | 0.4576209 | 0.02973368 | 3.5280296e-05 | 17.225061 | 0.009263834 | 0.5547573 | 0.097571984 | 0.0005015383 |
| 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.
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']))
indivObs = service.search(query).to_table()
indivObs
| 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 | -48611916884959414 | 21.243 | 0.038 | i | 61051.18583 | 14.526479 | 1.8827857 | 2.685296 |
| 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.
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.
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.
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.
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.
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()
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.
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()
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.
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()
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).
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()
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
Gparameter ofHG_model(average measured value - intrinsic value ofG).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.