305.1. Variable star light curves

305.1. Variable star light curves¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 1
Container Size: large
LSST Science Pipelines version: r29.2.0
Last verified to run: 2025-09-02
Repository: github.com/lsst/tutorial-notebooks

Learning objective: Learn how to extract and plot a light curve for a variable star from DP1.

LSST data products: ForcedSourceOnDiaObject, DiaObject, Visit

Packages: matplotlib, numpy, lsst.rsp, lsst.utils.plotting

Credit: Originally developed by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

This notebook demonstrates how to obtain and plot light curves of a variable star captured in Data Preview 1 (DP1). Time series are extracted for forced photometry on difference images and direct ("visit") images.

Related tutorials: The 100-level tutorials on the TAP service, and the 200-level tutorials on the difference image analysis (DIA) catalogs.

1.1. Import packages¶

Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), and matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery). Also import the units module from Astropy, which will be used to convert fluxes to magnitudes.

From the lsst package, import the module for accessing the Table Access Protocol (TAP) service, and some plotting utilities.

In [1]:

import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u

from lsst.rsp import get_tap_service

from lsst.utils.plotting import (get_multiband_plot_colors,
                                 get_multiband_plot_symbols)

1.2. Define parameters and functions¶

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

In [2]:

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

Define filter names, plot markers, and colors for plotting.

In [3]:

filter_names = ['u', 'g', 'r', 'i', 'z', 'y']
filter_colors = get_multiband_plot_colors()
filter_symbols = get_multiband_plot_symbols()

2. Retrieve light curve data¶

Obtain light curve data from the ForcedSourceOnDiaObject table, which contains forced PSF photometry for all difference images as well as non-difference images (visit_images). The relevant table columns for plotting light curves are:

• psfDiffFlux: Forced PSF photometry on difference images at DiaObject position
• psfFlux: Forced PSF photometry on visit_images at DiaObject position
• diaObjectId: Unique DiaObject identifier
• band: Filter associated with flux measurement
• visit: Identifier of the visit where the forced photometry was measured (used for table JOIN on Visit table)

The date information expMidptMJD, which is mid-point time for the visit in MJD, can be obtained from the Visit table.

The units used to plot the light curve are in flux (nJy) and not magnitudes since the negative flux measurements would be omitted when converting to magnitudes.

2.1 Look up diaObjectId from coordinates¶

The first step to obtain light curve data is to define the coordinates (ra and dec) of a variable star captured in DP1: Gaia DR3 2912281258855051520.

In [4]:

ra_var = 94.9226329830
dec_var = -25.2318482104

Define a 0.5 arcsecond search radius, and convert it to degrees.

In [5]:

search_rad = 0.5/3600

Search at this position in the DiaObject table to retrieve the diaObjectId of the corresponding object.

In [6]:

query = "SELECT ra, dec, diaObjectId "\
        "FROM dp1.DiaObject "\
        "WHERE CONTAINS(POINT('ICRS', ra, dec), "\
        "CIRCLE('ICRS', " + str(ra_var) + ", "\
        + str(dec_var) + ", " + str(search_rad) + ")) = 1"
print(query)

SELECT ra, dec, diaObjectId FROM dp1.DiaObject WHERE CONTAINS(POINT('ICRS', ra, dec), CIRCLE('ICRS', 94.922632983, -25.2318482104, 0.0001388888888888889)) = 1

In [7]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()

Job phase is COMPLETED

Fetch the results of the TAP search in table form and print the unique entries for diaObjectId in the output table.

In [8]:

assert job.phase == 'COMPLETED'
DiaObj = job.fetch_result().to_table()
print(np.unique(DiaObj['diaObjectId']))

   diaObjectId
------------------
614435753623027782

2.2. Search based on diaObjectId¶

Provide this diaObjectId as a constraint to a TAP query to obtain the time series from ForcedSourceOnDiaObject. The use of an ID in this manner is always the recommended way to access data from a ForcedSource table.

In [9]:

DiaObjID = DiaObj['diaObjectId'][0]

In [10]:

query = "SELECT fsodo.coord_ra, fsodo.coord_dec, "\
        "fsodo.diaObjectId, fsodo.visit, fsodo.band, "\
        "fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, "\
        "fsodo.psfFlux as psfFlux, fsodo.psfFluxErr, "\
        "vis.expMidptMJD "\
        "FROM dp1.ForcedSourceOnDiaObject as fsodo "\
        "JOIN dp1.Visit as vis ON vis.visit = fsodo.visit "\
        "WHERE fsodo.diaObjectId = "+str(DiaObjID)
print(query)

SELECT fsodo.coord_ra, fsodo.coord_dec, fsodo.diaObjectId, fsodo.visit, fsodo.band, fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, fsodo.psfFlux as psfFlux, fsodo.psfFluxErr, vis.expMidptMJD FROM dp1.ForcedSourceOnDiaObject as fsodo JOIN dp1.Visit as vis ON vis.visit = fsodo.visit WHERE fsodo.diaObjectId = 614435753623027782

Run the TAP query.

In [11]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()

Job phase is COMPLETED

Fetch the results of the TAP search in table form and print the number of rows.

In [12]:

assert job.phase == 'COMPLETED'
forced_source = job.fetch_result().to_table()
print(len(forced_source))

225

Uncomment the following to view the table.

In [13]:

# forced_source

Count the number of visits in each band and print them to the screen.

In [14]:

for band in filter_names:
    print(f"{band}: {np.sum((forced_source['band'] == band))} visits")

u: 20 visits
g: 66 visits
r: 71 visits
i: 17 visits
z: 43 visits
y: 8 visits

3. Plot the light curve¶

3.1. Plot flux vs. time¶

Plot the forced PSF photometry on the difference images and the direct (visit_image) images covering the variable star.

In [15]:

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7, 7), sharex=True)
fig.subplots_adjust(wspace=0, hspace=0)

for band in filter_names:
    pickband = np.where(forced_source['band'] == band)
    ax1.plot(forced_source[pickband]['expMidptMJD'], forced_source[pickband]['psfDiffFlux'],
             marker=filter_symbols[band], color=filter_colors[band],
             label=f"{band}", markersize=7, linestyle='none', fillstyle='none')
    ax2.plot(forced_source[pickband]['expMidptMJD'], forced_source[pickband]['psfFlux'],
             marker=filter_symbols[band], color=filter_colors[band], linestyle='none',
             markersize=7, fillstyle='none')

ax1.set_ylabel('PSF difference image flux (nJy)')
ax2.set_ylabel('PSF direct image flux (nJy)')
ax2.set_xlabel('exposure midpoint, MJD (d)')
ax1.minorticks_on()
ax2.minorticks_on()
ax1.legend(ncols=2)
plt.show()

Figure 1: Forced PSF photometry of a variable star in DP1 from the ForcedSourceonDiaObject table. The two panels represent the flux measured on the difference images (top) and the "direct" images (visit_image dataset type).

3.2. Plot flux vs. phase¶

Use the known period of the variable star to "phase-fold" the light curve.

Define the period in days for the SX Phe star from Carlin+2025.

In [16]:

period = 0.07670835

Place each point of the time series into the correct phase of the variable's pulsation. Use the time of maximum $r$-band flux as zero phase.

In [17]:

pickr = np.where(forced_source['band'] == 'r')
ind_max = np.argmax(forced_source['psfFlux'][pickr])
t0 = forced_source['expMidptMJD'][pickr][ind_max]
mjd_norm = (forced_source['expMidptMJD'] - t0) / period
phase = np.mod(mjd_norm, 1.0)

Plot the phase-folded light curves in difference-image and direct-image flux, then plot the curve in magnitudes (measured from the direct images).

In [18]:

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7, 7), sharex=True)
fig.subplots_adjust(wspace=0, hspace=0)

for band in filter_names:
    pickband = np.where((forced_source['band'] == band)
                        & (forced_source['psfDiffFlux'] > -50000))
    ax1.plot(phase[pickband], forced_source[pickband]['psfDiffFlux'],
             marker=filter_symbols[band], color=filter_colors[band],
             label=f"{band}", markersize=7, linestyle='none', fillstyle='none')
    ax2.plot(phase[pickband], forced_source[pickband]['psfFlux'],
             marker=filter_symbols[band], color=filter_colors[band], linestyle='none',
             markersize=7, fillstyle='none')

ax1.set_ylabel('PSF difference image flux (nJy)')
ax2.set_ylabel('PSF direct image flux (nJy)')
ax2.set_xlabel('phase')
ax1.minorticks_on()
ax2.minorticks_on()
ax1.legend(ncols=2)
plt.show()

Figure 2: Forced PSF photometry of a variable star in DP1 from the ForcedSourceonDiaObject table, with each point placed into its observed phase in the variable's pulsation cycle. The two panels represent the flux measured on the difference images (top) and the "direct" images (visit_image dataset type).

3.3. Plot magnitude vs. phase¶

Convert the direct image forced fluxes to magnitudes, then plot a light curve.

Warning: the following code cell produces a pink RuntimeWarning about invalid values encountered in log10, which represents the danger of converting fluxes to magnitudes: direct-image forced fluxes can be negative.

In [19]:

mag = forced_source['psfFlux'].to(u.ABmag).value

/opt/lsst/software/stack/conda/envs/lsst-scipipe-10.1.0/lib/python3.12/site-packages/astropy/units/function/logarithmic.py:67: RuntimeWarning: invalid value encountered in log10
  return dex.to(self._function_unit, np.log10(x))

In [20]:

fig, ax = plt.subplots(1, 1, figsize=(7, 5))

for band in filter_names:
    pickband = np.where(forced_source['band'] == band)
    ax.plot(phase[pickband], mag[pickband],
            marker=filter_symbols[band], color=filter_colors[band],
            label=f"{band}", markersize=7, linestyle='none', fillstyle='none')

ax.set_ylabel('PSF magnitude')
ax.set_xlabel('phase')
ax.minorticks_on()
ax.legend(ncols=2, loc='lower left')
plt.ylim(19.9, 18.2)
plt.show()

Figure 3: Forced PSF photometry, in AB magnitudes, of a variable star in DP1 from the direct-image measurements in the ForcedSourceonDiaObject table. Each point is placed at its observed phase in the variable's pulsation cycle.

In [ ]: