04a. Introduction to Phase Curves

Phase Curve of Solar System Objects
Contact authors: Christina Williams and Yumi Choi
Last verified to run: 2024-01-31
LSST Science Piplines version: Weekly 2024_04
Container Size: medium
Targeted learning level: intermediate

Description: Explore phase curves for DP0.3 solar system objects.

Skills: Use the TAP service and ADQL to access the DP0.3 tables. Join information from multiple DP0.3 tables. Plot phase curves.

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). 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 explores the properties of solar system bodies in the DP0.3 dataset by using the science example of constructing and exploring their phase curves.

The DP0.3 catalogs contain both real and simulated solar system objects (including asteroids, near-Earth objects, Trojans, trans-Neptunian objects, comets). These objects change position between each Rubin image. In the DP0.3 catalogs, the intrinsic properties and orbital parameters are known, and are used to estimate what the measurements would be in a given image, and how they change between images. From these simulated observations, it is possible to reconstruct their intrinsic properties and orbital parameters in the same way as will be done using the real LSST data. An important way to characterize intrinsic properties of a solar system object is by measuring its "phase curve", which is the object brightness as a function of its "solar phase angle" (the angle made between the line of sight from the object to the Sun, and the line of sight from the object to Earth; see Figure 1.

Figure 1: Diagram showing the definition of phase angle $\alpha$, credit Buchheim 2010

In order to reveal the intrinsic properties of the asteroid (such as its surface properties and albedo, which helps determine its class of solar system body), we first must turn apparent magnitudes as a function of time (what is measured by LSST data) into "reduced magnitude", which takes into account the relative distances between the asteroid and the Sun/Earth (heliocentric/topocentric distances) at each observation. Reduced magnitude is normalized such that it is the brightness of an asteroid as if it is observed at 1 astronomical unit (au) from both the Sun and the Earth (i.e. unit topocentric and heliocentric distance). Note that rotation curves or complex geometry of solar system objects are not included in DP0.3 simulations. Thus, any changes over time in an object’s apparent magnitude are due only to changes in its distance and phase angle. Phase curves can be constructed for each filter.

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. In this tutorial, we will mention three functional forms that are relevant for understanding the DP0.3 data products. These are the HG_model, HG1G2_model, and HG12_model. The HG_model is the simplest model (see Bowell et al. 1989), 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$. The HG_model $G$ and $H$ values are stored in the dp03_catalogs_10yr.MPCORB table.

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 will demonstrate how to plot the phase curve of solar system bodies using the DP0.3 simulated catalogs. Section 2 compares the phase curves for three different types of objects at different orbital radii from the Sun, and compares the data to the LSST pipeline measured phase curve parameters, H, and G12. Section 3 aggregates the phase curve fits for a number of solar system bodies in DP0.3 and shows how the quality of the fit depends on LSST observations (which additionally provides some insight into expectations for real LSST data).

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 seaborn package provides statistical data visualization with aesthetic and informative graphics.

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

In [1]:

import numpy as np
from scipy.interpolate import CubicSpline
import matplotlib.pyplot as plt
import seaborn as sns
from sbpy.photometry import HG12

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¶

Define a function to convert a given perihelion distance (q) and eccentricity (e) to an orbital semi-major axis (a). Their relationship is defined by $q = a(1-e)$.

In [4]:

def calc_semimajor_axis(q, e):
    """
    Given a perihelion distance and orbital eccentricity,
    calculate the semi-major axis of the orbit.

    Parameters
    ----------
    q: ndarray
        Distance at perihelion, in au.
    e: ndarray
        Orbital eccentricity.

    Returns
    -------
    a: ndarray
        Semi-major axis of the orbit, in au.
        q = a(1-e), so a = q/(1-e)
    """

    return q / (1.0 - e)

2. Plot the phase curves of asteroids in DP0.3¶

2.1 Create the Rubin TAP Service Client¶

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

In [5]:

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

2.2 Querying the DP0.3 SSObject and MPCORB catalogs¶

To plot an object's phase curve, we will need its apparent magnitudes and uncertainties across a range of observation times. Additionally, we need the phase angles of the object at those observation times, topocentric ($d_t$) and heliocentric ($d_h$) distances.

To define those mock observations of solar system objects, the DP0.3 model uses the HG_model form of the phase curve to predict the observed parameters (apparent magnitude) 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 purposes of DP0.3, the intrinsic slope, mpcG, for all objects have a constant value of 0.15.

In the SSObject table, the LSST data products contain the fitted phase curve parameters based on the mock observations using the HG12_model (i.e., the fitted absolute magnitude H and slope parameter G12 based on mock observations in $griz$ bands). Note that the value of G12 slope will differ from G owing to the difference in functional form (see Section 1.2). The explanation for the absence of $u$ and $y$ bands in DP0.3 catalogs can be found in the DP0.3 documentation.

Limit the query to select sources with a large number of observations that likely have good phase curve fits (in the SSObject table, this is numObs > 100) and with perihelion distance of less than 50 au (in the MPCORB table, q < 50), and then further narrow down the sample to focus on Main Belt asteroids, Jupiter Trojans, and Trans-Neptunian Obects (TNOs) for this tutorial.

Define the table returned by this query as "unique" since it contains the IDs of unique solar system objects. Each object has multiple individual observations in each filter in LSST.

In [6]:

nobs_thrh = '100'
q_thrh = '50'

In [7]:

query = """
SELECT
    mpc.ssObjectId, mpc.e, mpc.q, mpc.mpcG, mpc.mpcH,
    sso.arc, sso.numObs,
    sso.g_H, sso.g_Herr, sso.g_G12, sso.g_G12err,
    sso.g_H_gG12_Cov, sso.g_Ndata, sso.r_H, sso.r_Herr,
    sso.r_G12, sso.r_G12err, sso.r_H_rG12_Cov, sso.r_Ndata,
    sso.i_H, sso.i_Herr, sso.i_G12, sso.i_G12err, sso.i_H_iG12_Cov,
    sso.i_Ndata, sso.z_H, sso.z_Herr, sso.z_G12, sso.z_G12err,
    sso.z_H_zG12_Cov, sso.z_Ndata
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 < {} ORDER by sso.ssObjectId
""".format(nobs_thrh, q_thrh)

uniqueObj = service.search(query).to_table()
uniqueObj

Out[7]:

Table length=2489134

ssObjectId	e	q	mpcG	mpcH	arc	numObs	g_H	g_Herr	g_G12	g_G12err	g_H_gG12_Cov	g_Ndata	r_H	r_Herr	r_G12	r_G12err	r_H_rG12_Cov	r_Ndata	i_H	i_Herr	i_G12	i_G12err	i_H_iG12_Cov	i_Ndata	z_H	z_Herr	z_G12	z_G12err	z_H_zG12_Cov	z_Ndata
		AU		mag	d		mag	mag	mag	mag	mag2		mag	mag	mag	mag	mag2		mag	mag	mag	mag	mag2		mag	mag	mag	mag	mag2
int64	float64	float64	float32	float32	float32	int32	float32	float32	float32	float32	float32	int32	float32	float32	float32	float32	float32	int32	float32	float32	float32	float32	float32	int32	float32	float32	float32	float32	float32	int32
-9223365099512229794	0.07307953303637367	4.792882289504846	0.15	15.19	3621.013	212	15.7715225	0.06386547	0.97223943	0.38761416	0.024255332	30	15.240818	0.028765978	0.4542997	0.15985987	0.004452128	71	15.203755	0.041300815	0.91806996	0.23844625	0.009588541	69	15.198534	0.076746725	0.8757064	0.4624414	0.034198146	42
-9223354132337021199	0.28382716556600585	2.2297244878397118	0.15	16.73	2832.8926	211	17.455976	0.007145925	0.62499183	0.04307446	0.00028402833	28	16.807323	0.0025084645	0.5700074	0.013169541	2.7217959e-05	63	16.595743	0.0028806308	0.5970628	0.0149014965	3.262696e-05	67	16.653269	0.005202933	0.59098387	0.029086128	0.00011901059	53
-9223352335202182781	0.207805099328402	1.802932279106067	0.15	19.41	3204.884	222	20.145678	0.011630123	0.5147999	0.07490258	0.0004703224	28	19.488098	0.0098226555	0.54235405	0.05792951	0.00042414936	78	19.294615	0.012494921	0.54777527	0.07882866	0.0007578876	78	19.336838	0.023139708	0.4519446	0.1240977	0.002346807	38
-9223343361837155588	0.05210057580115273	2.22249354881098	0.15	17.36	3171.7117	269	18.078428	0.00477421	0.6356637	0.030178472	9.009757e-05	27	17.428768	0.0022698361	0.55891764	0.014113579	2.0452382e-05	84	17.234886	0.0028912583	0.5182386	0.01804207	2.9392451e-05	90	17.288725	0.006028117	0.5423687	0.03401306	0.0001460805	68
-9223341751461022808	0.20525093475763315	1.8029228289455757	0.15	19.3	3159.682	144	20.017632	0.025277533	0.34265444	0.13182785	0.0014116154	15	19.394728	0.011760565	0.5220102	0.08267982	0.00045152276	61	19.211948	0.016445434	0.44210672	0.10456468	0.0011191351	49	19.156347	0.04111035	0.20000058	0.22677861	0.0061788773	19
-9223340293227439626	0.04916731373423095	2.59437247949221	0.15	17.05	1956.854	212	17.794104	0.007732349	0.68999124	0.044719346	0.00029694513	34	17.12589	0.003888508	0.60213214	0.02359108	6.663622e-05	72	16.92747	0.0050112433	0.6473837	0.03259541	0.000108004766	57	16.96853	0.01102775	0.5281418	0.067947574	0.0005537676	49
-9223339508945339754	0.18794223627487194	1.9227222647783788	0.15	19.14	2583.77	161	19.876463	0.010376548	0.6686587	0.058875464	0.00022137184	19	19.213957	0.007350991	0.56991196	0.0425632	0.00021779325	62	19.00801	0.008321819	0.5324443	0.048095364	0.0002580475	55	19.055206	0.015550898	0.58402115	0.09111061	0.000479817	25
-9223333218989556101	0.030064887923693176	2.7960357000662928	0.15	15.17	3255.8918	305	15.878582	0.0017741049	0.54394245	0.010351132	1.4347317e-05	45	15.237398	0.0009780133	0.5801164	0.0059684855	3.9648826e-06	95	15.04045	0.00126829	0.57148296	0.007886324	6.286281e-06	92	15.080025	0.0024144351	0.5277001	0.0138973575	2.6122016e-05	73
-9223315625959504597	0.06772062985106744	2.8658655050926103	0.15	15.97	3586.1482	410	16.704933	0.0048450576	0.588174	0.030300928	0.000113880145	35	16.049337	0.0017468001	0.5926988	0.013372099	1.2332018e-05	131	15.848513	0.0024199192	0.5152747	0.021554343	2.2886024e-05	132	15.906233	0.004723713	0.6144167	0.04187281	8.729006e-05	112
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
9223278184854652110	0.08618164232287248	2.40238377582261	0.15	19.52	2928.6963	105	20.201815	0.12893167	0.45968637	0.6014887	0.07533547	12	19.631618	0.022837553	0.7865213	0.14089221	0.0025065695	46	19.339828	0.035211023	0.20000163	0.18320023	0.005566555	34	19.353107	0.10022238	0.038898293	0.2680323	-0.023706537	13
9223281763650087600	0.05490308401541161	2.2129231702223775	0.15	19.25	2960.2139	149	19.806181	0.018485172	0.61418945	0.12059271	0.0014815063	26	19.31811	0.010659503	0.56265324	0.06728714	0.00044201856	58	19.189806	0.01725591	0.45611915	0.10976626	0.0011951937	43	19.178862	0.034777395	0.74127024	0.23181646	0.0041091023	22
9223287598582163833	0.31296411815497166	2.168844808705169	0.15	16.76	3202.0093	127	17.313372	0.004779948	0.585819	0.046645086	0.00016595656	18	16.837505	0.0018107332	0.5070907	0.019420346	1.5750542e-05	49	16.721825	0.0027882745	0.5825473	0.026204742	3.7817146e-05	33	16.71867	0.004648033	0.5453388	0.051939163	0.00011747823	27
9223287675055310733	0.08703352308888852	2.0969441134957676	0.15	19.32	3169.8274	108	20.019531	0.024438698	0.5568655	0.13325709	0.0026610433	21	19.393578	0.014095212	0.54126483	0.079773195	0.00080782955	38	19.195377	0.01880317	0.56129587	0.10931974	0.0015270964	29	19.29349	0.032630686	0.8370995	0.20273608	0.004804776	20
9223296838638448722	0.14764771534731982	2.1508069571429567	0.15	16.11	3494.9194	241	16.659655	0.0015415418	0.56145924	0.008313641	8.021681e-06	38	16.183912	0.0011536715	0.55148244	0.0063217715	3.9874903e-06	74	16.06012	0.0016653342	0.5469504	0.010476721	8.019832e-06	68	16.05897	0.002486881	0.5056896	0.014769626	1.7281367e-05	61
9223298685954327106	0.2059469925979121	2.1331831755454678	0.15	17.63	3503.6614	380	18.368998	0.0032473137	0.47447705	0.021966053	2.0840052e-06	50	17.710314	0.0016728243	0.48909396	0.011891814	2.567575e-06	123	17.520002	0.002101392	0.5063691	0.015854582	4.4191906e-06	131	17.566658	0.005242885	0.5470947	0.035566498	8.7221975e-05	76
9223311537846468385	0.0765472199248066	2.86547027457465	0.15	18.65	3156.943	111	19.351326	0.06095468	0.5438121	0.3508283	0.019637985	23	18.683033	0.038138554	0.5785368	0.22219639	0.007909983	38	18.547647	0.048371747	0.68295443	0.25766724	0.011507387	38	18.649014	0.14814813	0.68391824	0.771532	0.10888628	12
9223337017087558055	0.26681940940283416	2.298783342116851	0.15	17.7	2414.8093	196	18.246964	0.015251029	0.6186255	0.088942036	0.0006744807	28	17.790382	0.00766924	0.54401815	0.04670926	0.00020148733	87	17.658634	0.0118442	0.54106206	0.0760863	0.00025749207	60	17.702526	0.03984432	0.59096545	0.25900084	0.0049145697	21
9223345276844401162	0.0184994877457838	2.8084919457165305	0.15	17.69	3250.853	192	18.228022	0.014724213	0.49944615	0.08959828	0.0011144761	34	17.747402	0.008122387	0.5772704	0.054273263	0.00031237214	60	17.605227	0.012395228	0.54061943	0.07824657	0.00070312025	56	17.67629	0.025864098	0.7409961	0.15452531	0.0033356848	42
9223352613567217422	0.2242284731621877	1.8231970965726527	0.15	17.88	3111.7158	209	18.599514	0.0050439574	0.53474814	0.02576561	8.931978e-05	36	17.949776	0.0037079619	0.5403619	0.019228322	4.8780264e-05	69	17.758282	0.006202687	0.5732125	0.035943147	0.0001712879	62	17.813639	0.013187968	0.61169356	0.077705696	0.00078590924	42

2.3 Identify Main Belt asteroids, Jupiter Trojans, and TNOs¶

First calculate the semi-major axis (a) of each object's orbit. Second, construct a Main Belt asteroid sample, a Jupiter Trojan sample, and a Trans-Neptunian Object (TNO) sample following the population definitions used by the JPL Horizons small body database query tool (https://ssd.jpl.nasa.gov/tools/sbdb_query.html): Main Belt asteroids (2.0 < a < 3.25 au and q > 1.666 au); Jupiter Trojans (4.6 < a < 5.5 au and e < 0.3), and finally TNOs which are at a > 30.1 au. We further limit our sample to the asteroids with numObs > 2000, to focus on objects with many individual observations to define its phase curve. (Note that the query for unique objects in Section 2.2 used a less-restrictive numObs > 100, which will used to explore the empirical properties that can influence phase curve fitting later in Section 3 and the followup tutorial notebook 04b).

In [8]:

a = calc_semimajor_axis(uniqueObj['q'], uniqueObj['e'])

main_belt = (a > 2.0) & (a < 3.2) & (uniqueObj['q'] > 1.666)
main_belt_2k = (uniqueObj['numObs'] > 2000) & main_belt

JT = (a > 4.6) & (a < 5.5) & (uniqueObj['e'] < 0.3)
JT_2k = (uniqueObj['numObs'] > 2000) & JT

TNO = (a > 30.1)
TNO_2k = (uniqueObj['numObs'] > 2000) & TNO

print('There are %d asteroids in this sample!' % len(a[main_belt_2k]))
print('There are %d Jupiter Trojans in this sample!' % len(a[JT_2k]))
print('There are %d TNOs in this sample!' % len(a[TNO_2k]))

all_2k = main_belt_2k + JT_2k + TNO_2k

There are 296 asteroids in this sample!
There are 124 Jupiter Trojans in this sample!
There are 38 TNOs in this sample!

Visualize the distributions of the semi-major axis a for the all unique objects in the 3 different populations.

In [9]:

h = plt.hist(a[main_belt_2k], bins=25, range=(1, 7),
             label='Main Belt asteroids', alpha=.5)
h = plt.hist(a[JT_2k], bins=25, range=(1, 7),
             label='Jupiter Trojans', alpha=.5)
h = plt.hist(a[TNO_2k], bins=50, range=(1, 60), label='TNOs', alpha=.5)

plt.xlabel('Semi-major axis [au]')
plt.ylabel('Number')
plt.yscale('log')
plt.legend()
plt.show()

2.4 Identify individual mock observations of the solar system bodies in 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 in each filter 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 solar system bodies that have >2000 observations as selected above.

In [10]:

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'][all_2k]))

indivObs = service.search(query).to_table()
indivObs

Out[10]:

Table length=1104253

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 the asteroids with > 2000 observations in uniqueObj, as they should be.

In [11]:

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

2.5 Plot the phase curve for a unique Main Belt asteroid¶

To plot the phase curve, it is necessary to first compute the reduced magnitude $H(\alpha)$ for each observation in each filter ($g,r,i,z$ for DP0.3), and add it as a column to the indivObs table of individual observations. The reduced magnitude $H(\alpha)$ as mentioned in Introduction 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.

Note that the phase curves of solar system objects vary by observed filter, due to a combination of reflectivity of the surface and the intrinsic solar spectrum. For reference, to define the intrinsic absolute magnitude in DP0.3, the wavelength used to input the intrinsic absolute magnitude is in the $V$ band. Thus, one should not expect the absolute magnitude from any one filter to line up with the "truth" absolute magnitude, since none of the LSST filters include $V$ band. The conversion between each of the LSST filters and $V$ are documented for solar system bodies, and can be found in the DP0.3 documentation.

In [12]:

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

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

The cell below will plot example phase curves in all available filters (in DP0.3) $g$,$r$,$i$,$z$ for a single object referenced by its ssObjectId, here called sId. (Explore different objects by changing the iObj index to retrieve different sources). Below, it is clear 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.

It is possible to pick an integer number between 0 and len(uniqueObj[main_belt_2k])-1 for iObj below to explore other objects.

The HG12_model-predicted phase curve is also overplotted to visualize the LSST pipeline phase curve fitting (dashed line). This demonstrates how absolute magnitude H is determined in the LSST data products.

In [13]:

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


for i, ifilt in enumerate(filts):
    idx = tmp['band'] == ifilt

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

    # Plot HG12 model
    HG12_mag = HG12.evaluate(np.deg2rad(phases), 
                             uniqueObj[ifilt+'_H'][main_belt_2k][iObj], 
                             uniqueObj[ifilt+'_G12'][main_belt_2k][iObj])
    plt.plot(phases, HG12_mag, color=filter_colors[ifilt], linestyle='--',
             label='HG12 2-parameter model in %s' % ifilt)    


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', ncol=1)
plt.title('Main Belt Asteroid')
plt.show()

The cell below compares the phase curve of one object from each of the three populations selected in Section 2.3: the Main Belt asteroids, with the more distant Jupiter Trojans and the TNOs. The color scheme is the same as the cell above. This cell must be executed after running the above cell first.

In [14]:

plt.figure(figsize=(13, 4))

plt.subplot(131)
for i, ifilt in enumerate(filts):
    idx = tmp['band'] == ifilt

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

    # Plot HG12 model
    HG12_mag = HG12.evaluate(np.deg2rad(phases), 
                             uniqueObj[ifilt+'_H'][main_belt_2k][iObj],
                             uniqueObj[ifilt+'_G12'][main_belt_2k][iObj])
    plt.plot(phases, HG12_mag, color=filter_colors[ifilt], linestyle='--')


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.title('Main Belt Asteroid')

plt.subplot(132)

iObj = 2
sId = uniqueObj['ssObjectId'][JT_2k][iObj]
tmp = indivObs[indivObs['ssObjectId'] == sId]

for i, ifilt in enumerate(filts):
    idx = tmp['band'] == ifilt

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

    # Plot HG12 model
    HG12_mag = HG12.evaluate(np.deg2rad(phases),
                             uniqueObj[ifilt+'_H'][JT_2k][iObj],
                             uniqueObj[ifilt+'_G12'][JT_2k][iObj])
    plt.plot(phases, HG12_mag, color=filter_colors[ifilt], linestyle='--')


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.title('Jupiter Trojans')

plt.subplot(133)
iObj = 3
sId = uniqueObj['ssObjectId'][TNO_2k][iObj]
tmp = indivObs[indivObs['ssObjectId'] == sId]

for i, ifilt in enumerate(filts):
    idx = tmp['band'] == ifilt

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

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.title('TNO')

plt.tight_layout()
plt.show()

Notice that the the phase coverage of objects at different semi-major axes are dramatically different. The most distant (>30.1 au) TNO will cover at most ~0-3 degrees in phase, and this is expected because they are too distant to have long phase angle coverage. Meanwhile, asteroids in the main belt will typically achieve between 0-30 degrees in terms of phase coverage. The LSST pipeline will only attempt to fit phase curves to things with more than 5 degrees of phase coverage, so the SSObject table does not include H, G12 parameters for TNOs. Thus, TNOs do not have any models over-plotted. In these cases, it could be useful to independently model the phase curve using user-defined functions (e.g. for TNOs, a linear fit can be sufficient).

2.6 Explore phase curve fit uncertainty in each filter¶

For one of the Main Belt asteroids with > 2000 observations, the plot below compares the $g$- and $z$-bands to demonstrate fit uncertainty between two filters that produce different reduced magnitude quality (mostly due to difference in brightness and therefore flux uncertainties between the filters). The cell also computes the model phase curve for a given set of the best-fit phase curve parameters for the HG12_model stored in SSObject table, which is called HG12_mag_sso. The plot shows the uncertainty in the model parameters represented by the shaded regions. Choosing fainter or less well-sampled SSObjects increases the error region.

In [15]:

iObj = 3
sId = uniqueObj['ssObjectId'][main_belt_2k][iObj]
tmp = indivObs[indivObs['ssObjectId'] == sId]

for ifilt in ['g', 'z']:
    idx = tmp['band'] == ifilt
    plt.errorbar(tmp['phaseAngle'][idx], tmp['reducedMag'][idx],
                 yerr=tmp['magErr'][idx], fmt='o',
                 color=filter_colors[ifilt], label=ifilt, zorder=10)

    HG12_mag_sso = HG12.evaluate(np.deg2rad(phases),
                                 uniqueObj[ifilt+'_H'][main_belt_2k][iObj],
                                 uniqueObj[ifilt+'_G12'][main_belt_2k][iObj])
    plt.plot(phases, HG12_mag_sso, '--', alpha=0.3, color=filter_colors[ifilt])

    # Compute min/max values in reduced mag at each phase angle
    p1 = HG12.evaluate(np.deg2rad(phases),
                       uniqueObj[ifilt+'_H'][main_belt_2k][iObj] + 3*uniqueObj[ifilt+'_Herr'][main_belt_2k][iObj],
                       uniqueObj[ifilt+'_G12'][main_belt_2k][iObj] + 3*uniqueObj[ifilt+'_G12err'][main_belt_2k][iObj])

    p2 = HG12.evaluate(np.deg2rad(phases), 
                       uniqueObj[ifilt+'_H'][main_belt_2k][iObj] - 3*uniqueObj[ifilt+'_Herr'][main_belt_2k][iObj],
                       uniqueObj[ifilt+'_G12'][main_belt_2k][iObj] + 3*uniqueObj[ifilt+'_G12err'][main_belt_2k][iObj])

    p3 = HG12.evaluate(np.deg2rad(phases),
                       uniqueObj[ifilt+'_H'][main_belt_2k][iObj] + 3*uniqueObj[ifilt+'_Herr'][main_belt_2k][iObj],
                       uniqueObj[ifilt+'_G12'][main_belt_2k][iObj] - 3*uniqueObj[ifilt+'_G12err'][main_belt_2k][iObj])

    p4 = HG12.evaluate(np.deg2rad(phases),
                       uniqueObj[ifilt+'_H'][main_belt_2k][iObj] - 3*uniqueObj[ifilt+'_Herr'][main_belt_2k][iObj],
                       uniqueObj[ifilt+'_G12'][main_belt_2k][iObj] - 3*uniqueObj[ifilt+'_G12err'][main_belt_2k][iObj])

    HG_magHigh = np.maximum(np.maximum(p1, p2), np.maximum(p3, p4))
    HG_magLow = np.minimum(np.minimum(p1, p2), np.minimum(p3, p4))

    plt.fill_between(phases, HG_magLow, HG_magHigh, alpha=0.3, color=filter_colors[ifilt])

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

3. Exploring phase curve data products from the DP0.3 SSObject and MPCORB catalogs¶

3.1 Distribution of H and G12 parameters in observed filters¶

This section explores the distribution of typical values of the G12 slope parameter as a function of absolute magnitudes for Main Belt asteroids and Jupiter Trojans. Due to small phase angle coverage, nearly all of TNOs do not have measured phase curve parameters. Remember that the input (truth) G value using the HG_model that was used to generate the DP0.3 simulated object's observed properties was fixed across the population to a constant value of G=0.15. The DP0.3 automated phase curve fitting (which uses HG12_model) produces a nearly constant value for G12 with a relatively small spread at bright magnitudes, and the scatter in measured G12 starts to deviate more substantially at fainter magnitudes where its likely harder to recover the intrinsic value.

Warnings: the following cell produces RuntimeWarning messages about invalid values which can be ignored. They are due to taking the log of the occasional negative or NaN value in the catalog. Those points will not show in the plot.

In [16]:

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


for i, ifilt in enumerate(filts):
    sns.histplot(x=uniqueObj[ifilt+'_H'], y=uniqueObj[ifilt+'_G12'],
                 bins=100, log_scale=(False, True), ax=axs[i])
    axs[i].scatter(uniqueObj[ifilt+'_H'][main_belt],
                   uniqueObj[ifilt+'_G12'][main_belt],
                   s=3, c='m', alpha=0.3, label='Main Belt asteroids')
    axs[i].scatter(uniqueObj[ifilt+'_H'][JT],
                   uniqueObj[ifilt+'_G12'][JT],
                   s=3, c='k', alpha=0.1, label='Jupiter Trojans')
    axs[i].text(10, 5, ifilt+'-band')
    axs[i].set(xlabel=None)
    axs[i].set(ylabel=None)

axs[0].legend(loc=3, fontsize=10)
plt.xlim(9, 25)
plt.ylim(1e-2, 1e1)

fig.supxlabel('Absolute magnitude, H [mag]')
fig.supylabel('Slope parameter, G12')
print(len(uniqueObj), len(uniqueObj[main_belt_2k]),
      len(uniqueObj[JT_2k]))
plt.show()

2489134 296 124

3.2 Phase curve fit uncertainty vs. apparent brightness of the Main Belt asteroids¶

This section compares the typical apparent magnitude uncertainties per filter to see how that impacts the fit. First, compute the median apparent magnitude (i.e. as observed) and its median uncertainty in the $r$-band in this example. This enables plotting phase curve parameters vs. median apparent magnitude and its median uncertainty to see how these observational characteristics impact the resulting modeling. One can see in the plot that as objects get fainter and the median apparent magnitude uncertainty increases, so does the uncertainty in the fitted parameters of the phase curve.

In [17]:

mag_med = []
magSigma_med = []
ifilt = 'r'

for iobj in uniqueObj['ssObjectId'][main_belt_2k]:
    idx = indivObs['ssObjectId'] == iobj
    tmp = indivObs[idx]
    idx_filt = tmp['band'] == ifilt

    mag_med.append(np.ma.median(tmp['mag'][idx_filt]))
    magSigma_med.append(np.ma.median(tmp['magErr'][idx_filt]))

plt.figure(figsize=(10, 4))

plt.subplot(121)
plt.plot(mag_med, uniqueObj[ifilt+'_Herr'][main_belt_2k], '.', alpha=.5,
         label=ifilt+'-band H parameter')
plt.plot(mag_med, uniqueObj[ifilt+'_G12err'][main_belt_2k], '.', alpha=.5,
         label=ifilt+'-band G12 parameter')
plt.ylabel('Uncertainty in phase curve fit')
plt.xlabel('Median apparent mag')
plt.legend()

plt.subplot(122)
plt.plot(magSigma_med, uniqueObj[ifilt+'_Herr'][main_belt_2k],
         '.', alpha=.5)
plt.plot(magSigma_med, uniqueObj[ifilt+'_G12err'][main_belt_2k],
         '.', alpha=.5)
plt.ylabel('Uncertainty in phase curve fit')
plt.xlabel('Median uncertainty in apparent mag')

plt.tight_layout()
plt.show()

3.3 Phase curve fit uncertainty vs. total number of LSST observations and perihelion distance¶

This section explores the impact of the total number of observations for a given source (numObs) and the perihelion distance (q) on the quality of phase curve fitting in $i$-band as an example. In the left columns, it is clear that the model uncertainties decrease with number of observations of each source. So as LSST accumulates data over time, precision in the phase curve modeling will improve. The right columns show that uncertainties in the phase curve parameters modestly increase for objects with larger perihelion distances.

In [18]:

fig = plt.figure(figsize=(10, 6))
gs = fig.add_gridspec(2, 2, wspace=0, hspace=0.7)
axs = gs.subplots(sharey=True)
ifilt = 'i'

sns.histplot(x=uniqueObj['numObs'], y=uniqueObj[ifilt+'_Herr'], bins=100,
             log_scale=(False, True), ax=axs[0, 0])
axs[0, 0].set_xlabel('Number of LSST observations')
axs[0, 0].set_ylabel(r'$\Delta$H in %s band' % ifilt)
plt.ylim(1e-4, 3e1)

sns.histplot(x=uniqueObj['q'], y=uniqueObj[ifilt+'_Herr'],
             bins=100, ax=axs[0, 1])
axs[0, 1].set_xlabel('perihelion distance [au]')

sns.histplot(x=uniqueObj['numObs'], y=uniqueObj[ifilt+'_G12err'],
             bins=100, log_scale=(False, True), ax=axs[1, 0])
axs[1, 0].set_xlabel('Number of LSST observations')
axs[1, 0].set_ylabel(r'$\Delta$G12 in %s band' % ifilt)

sns.histplot(x=uniqueObj['q'], y=uniqueObj[ifilt+'_G12err'],
             bins=100, ax=axs[1, 1])
axs[1, 1].set_xlabel('perihelion distance [au]')
plt.show()

3.4 Number of datapoints used for phase curve fitting per band¶

The above plots compare numObs (total in all bands) with model fits (per band) which may not be the ideal metric since the quality of phase curves can vary quite a bit between filters. Instead, one can look at the number of datapoints included in the phase curve modeling on a per filter basis (i.e. rNdata for the r-band in the SSObject table). The next plot looks at the distribution of the number of observations in each filter used to model the phase curve per filter. One can see that generally, $r$- and $i$-bands produce the most data points for recovering phase curves, while $g$-band produces less and $z$-band produces the fewest. Phase curves measured in $r$- and $i$-bands will thus be better sampled.

Thus in the second panel, one can see that poorer sampling drives higher uncertainty in H using $z$-band instead of $r$-band for solar system objects.

In [19]:

fig, axs = plt.subplots(1, 2, figsize=(10, 4))

for i, ifilt in enumerate(filts):
    axs[0].hist(uniqueObj[ifilt+'_Ndata'], bins=100, range=(0, 1300),
                histtype='step', color=filter_colors[ifilt], label=ifilt+'-band')

axs[0].set_yscale('log')
axs[0].set_xscale('log')
axs[0].legend()
axs[0].set_xlabel('Number of data points per filter')
axs[0].set_ylabel('Number of sources')

sns.histplot(x=uniqueObj['r_Herr'], y=uniqueObj['z_Herr'],
             bins=100, log_scale=(True, True), color='b', ax=axs[1], alpha=0.7)

one2one = np.arange(1e-4, 1000, .1)
axs[1].plot(one2one, one2one, '--')
axs[1].set_ylim(1e-4, 1e2)
axs[1].set_xlim(1e-4, 1e2)
axs[1].set_ylabel(r'$\Delta$H using z-band')
axs[1].set_xlabel(r'$\Delta$H using r-band')

plt.tight_layout()
plt.show()

4. Excercises for the learner¶

1. Try out the more advanced companion tutorial notebook 04b on explicitly modeling phase curves and comparing with the HG12 model parameters stored in the DP0.3.