Intermediate: Finding Flares and Variable Stars in TASOC Light Curves#

This notebook demostrates MAST’s programmatic tools for accessing TESS time series data while exploring a flaring star from the literature (Günther et al 2019) and a nearby variable star.

The following topics will be covered:

  • Using the MAST API to get mission pipeline and TASOC light curves

  • Plotting TESS light curves in Python

  • Using the MAST API to make an full frame image (FFI) cutout

  • Creating a movie of TPF frames in Python

  • Using the MAST API to get a list of TESS Input Catalog (TIC) sources

  • Over plotting TIC sources on TESS images

  • Exploring the period of a variable star

Table of contents#

  1. Exploring a known stellar flare

    1. Querying MAST

    2. Downloading the light curves

    3. Plotting the light curves

    4. Making an animation

  2. Exploring a variable star

    1. Overlaying TESS Input Catalog sources

    2. Getting the variable star light curve

    3. Plotting the variable star light curve

    4. Finding the period of the variable star

  3. Additional Resources

  4. About this Notebook


Terminology#

  • TESS: The Transiting Exoplanet Survey Satellite

  • TASOC: The TESS Asteroseismic Science Operations Center

  • Sector: TESS observed the sky in regions of 24x96 degrees along the southern, then northern, ecliptic hemispheres. Each of these regions is referred to as a “sector”, starting with Sector 1.

  • TIC: The TESS input catalog.

  • FFI: TESS periodically reads out the entire frame of all four cameras, nominally every 30 minutes, and stores them as full frame images (FFIs).

  • TPF: Target Pixel File, a fits file containing stacks of small images centered on a target, one image for every timestamp the telescope took data.

  • HDU: Header Data Unit. A FITS file is made up of HDUs that contain data and metadata relating to the file. The first HDU is called the primary HDU, and anything that follows is considered an “extension”, e.g., “the first FITS extension”, “the second FITS extension”, etc.

  • HDUList: A list of HDUs that comprise a fits file.

  • BJD: Barycentric Julian Date, the Julian Date that has been corrected for differences in the Earth’s position with respect to the Solar System center of mass.

  • BTJD: Barycentric TESS Julian Date, the timestamp measured in BJD, but offset by 2457000.0. I.e., BTJD = BJD - 2457000.0

  • WCS: World Coordinate System, the coordinates that locate an astronomical object on the sky.


Imports#

In this notbook we will use the MAST module of Astroquery (astroquery.mast) to query and download data, and various astropy and numpy functions to manipulate the data. We will use both the matplotlib and bokeh plotting packages to visualize our data as they have different strengths and weaknesses.

# For querying for data
from astroquery.mast import Tesscut, Observations, Catalogs

# For manipulating data
import numpy as np

from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy.io import fits
from astropy.wcs import WCS
from astropy.timeseries import LombScargle
from astropy.time import Time
import astropy.units as u

# For matplotlib plotting
import matplotlib
%matplotlib inline

import matplotlib.pyplot as plt
import matplotlib.animation as animation

# For animation display
from matplotlib import rc
from IPython.display import HTML
rc('animation', html='jshtml')

# For bokeh plotting
from bokeh import plotting
from bokeh.models import Span
plotting.output_notebook()
Loading BokehJS ...

Exploring a stellar flare#

Selecting the flare#

We will start with a known flare from the literature, in this case from Günther, M. N., Zhan, Z., Seager, S., et al. 2019, arXiv e-prints, arXiv:1901.00443. We picked a particularly long flare to give us the best chance of finding it in the 30 minute cadence data as well as the 2 minute cadence data.

We’ve made note of the TIC ID and sector for our flare of interest, as well as its peak time in BJD format:

tic_id = 141914082
sector = 1

tpeak = 2458341.89227 # Julian Day

Querying MAST#

Mission light curves#

Using the query_criteria function in the astroquery.mast.Observations class, we will specify that we are looking for TESS mission data using the obs_collection argument, that we want a specific TIC ID using the target_name argument, and that we want observations from a particular sector using the sequence_number argument.

mission_res = Observations.query_criteria(obs_collection="TESS", 
                                          target_name=tic_id, 
                                          sequence_number=sector)
mission_res
Table masked=True length=1
intentTypeobs_collectionprovenance_nameinstrument_nameprojectfilterswavelength_regiontarget_nametarget_classificationobs_ids_ras_decdataproduct_typeproposal_picalib_levelt_mint_maxt_exptimeem_minem_maxobs_titlet_obs_releaseproposal_idproposal_typesequence_numbers_regionjpegURLdataURLdataRightsmtFlagsrcDenobsidobjID
str7str4str4str10str4str4str7str9str1str47float64float64str10str14int64float64float64float64float64float64str1float64str23str1int64str41str1str73str6boolfloat64str8str9
scienceTESSSPOCPhotometerTESSTESSOptical141914082--tess2018206045859-s0001-0000000141914082-0120-s94.6175354047675-72.0448462219924timeseriesRicker, George358324.7932369675958352.67632608796120.0600.01000.0--58458.5833333G011175_G011180_G011176--1CIRCLE 94.6175354 -72.04484622 0.00138889--mast:TESS/product/tess2018206045859-s0001-0000000141914082-0120-s_lc.fitsPUBLICFalsenan60829534110344410

TASOC light curves#

In addition to mission pipeline data, MAST also hosts a variety of community contributed High Level Science Products (HLSPs), all of which are given the mission (obs_collection) “HLSP”. In this case we will specifically search for HLSPs in the TESS project, which will return the light curves provided by the TASOC (note the provenance_name of “TASOC”). All other arguments remain the same.

tasoc_res = Observations.query_criteria(target_name=tic_id, 
                                        obs_collection="HLSP", 
                                        project="TESS",
                                        sequence_number=sector)
tasoc_res['dataproduct_type',"obs_collection","target_name","t_exptime","filters",
          "provenance_name","project","sequence_number","instrument_name"]
Table length=6
dataproduct_typeobs_collectiontarget_namet_exptimefiltersprovenance_nameprojectsequence_numberinstrument_name
str10str4str9float64str4str17str4int64str10
timeseriesHLSP1419140821800.0TESSQLPTESS1Photometer
timeseriesHLSP1419140821800.0TESSTESS-SPOCTESS1Photometer
timeseriesHLSP141914082120.0TESSTASOCTESS1Photometer
timeseriesHLSP1419140821800.0TESSTASOCTESS1Photometer
timeseriesHLSP1419140821800.0TESSGSFC-ELEANOR-LITETESS1Photometer
timeseriesHLSP1419140821800.0TESSTGLCTESS1Photometer

This query returns two light curves. To understand the difference between the two light curves we look at the t_exptime column, and note the different values. These exposure times correspond to 2 minutes (short cadence) and 30 minutes (long cadence). We will explore both light curves.

Downloading the light curves#

For the rest of this notebook we will work with the TASOC light curves only, although we could do the same analysis on the mission light curves.

Querying for the list of associated data products#

Each observation may have one or more associated data products. In the case of the TASOC light curves, there is simply a single light curve file for each observation.

tasoc_prod = Observations.get_product_list(tasoc_res)
tasoc_prod["dataproduct_type", "description", "dataURI", "size"]
Table length=9
dataproduct_typedescriptiondataURIsize
str10str4str129int64
timeseriesFITSmast:HLSP/gsfc-eleanor-lite/s0001/0000/0001/4191/4082/hlsp_gsfc-eleanor-lite_tess_ffi_s0001-0000000141914082_tess_v1.0_lc.fits106560
timeseriesFITSmast:HLSP/qlp/s0001/0000/0001/4191/4082/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc.fits80640
timeseriesTextmast:HLSP/qlp/s0001/0000/0001/4191/4082/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc.txt57710
timeseriesFITSmast:HLSP/tasoc/s0001/c1800/0000/0001/4191/4082/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_cbv-lc.fits164160
timeseriesFITSmast:HLSP/tasoc/s0001/c1800/0000/0001/4191/4082/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_ens-lc.fits167040
timeseriesFITSmast:HLSP/tasoc/s0001/c0120/0000/0001/4191/4082/hlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05_cbv-lc.fits1880640
timeseriesFITSmast:HLSP/tess-spoc/s0001/target/0000/0001/4191/4082/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_lc.fits164160
timeseriesFITSmast:HLSP/tess-spoc/s0001/target/0000/0001/4191/4082/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp.fits3925440
timeseriesFITSmast:HLSP/tglc/s0001/cam4-ccd2/0052/6627/0443/4424/hlsp_tglc_tess_ffi_gaiaid-5266270443442455040-s0001-cam4-ccd2_tess_v1_llc.fits57600

Downloading files#

We can choose to download some or all of the associated data files, in this case since we just have the two light curves, we will download all of the products.

tasoc_manifest = Observations.download_products(tasoc_prod)
tasoc_manifest
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/gsfc-eleanor-lite/s0001/0000/0001/4191/4082/hlsp_gsfc-eleanor-lite_tess_ffi_s0001-0000000141914082_tess_v1.0_lc.fits to ./mastDownload/HLSP/hlsp_gsfc-eleanor-lite_tess_ffi_s0001-0000000141914082_tess_v1.0_lc/hlsp_gsfc-eleanor-lite_tess_ffi_s0001-0000000141914082_tess_v1.0_lc.fits ...
 [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/qlp/s0001/0000/0001/4191/4082/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc.fits to ./mastDownload/HLSP/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc.fits ...
 [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/qlp/s0001/0000/0001/4191/4082/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc.txt to ./mastDownload/HLSP/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc.txt ...
 [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/tasoc/s0001/c1800/0000/0001/4191/4082/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_cbv-lc.fits to ./mastDownload/HLSP/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_cbv-lc.fits ...
 [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/tasoc/s0001/c1800/0000/0001/4191/4082/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_ens-lc.fits to ./mastDownload/HLSP/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_ens-lc.fits ...
 [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/tasoc/s0001/c0120/0000/0001/4191/4082/hlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05_cbv-lc.fits to ./mastDownload/HLSP/hlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05/hlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05_cbv-lc.fits ...
 [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/tess-spoc/s0001/target/0000/0001/4191/4082/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_lc.fits to ./mastDownload/HLSP/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_lc.fits ...
 [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/tess-spoc/s0001/target/0000/0001/4191/4082/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp.fits to ./mastDownload/HLSP/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp.fits ...
 [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/tglc/s0001/cam4-ccd2/0052/6627/0443/4424/hlsp_tglc_tess_ffi_gaiaid-5266270443442455040-s0001-cam4-ccd2_tess_v1_llc.fits to ./mastDownload/HLSP/hlsp_tglc_tess_ffi_gaiaid-5266270443442455040-s0001-cam4-ccd2_tess_v1_llc/hlsp_tglc_tess_ffi_gaiaid-5266270443442455040-s0001-cam4-ccd2_tess_v1_llc.fits ...
 [Done]
Table length=9
Local PathStatusMessageURL
str172str8objectobject
./mastDownload/HLSP/hlsp_gsfc-eleanor-lite_tess_ffi_s0001-0000000141914082_tess_v1.0_lc/hlsp_gsfc-eleanor-lite_tess_ffi_s0001-0000000141914082_tess_v1.0_lc.fitsCOMPLETENoneNone
./mastDownload/HLSP/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc.fitsCOMPLETENoneNone
./mastDownload/HLSP/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc.txtCOMPLETENoneNone
./mastDownload/HLSP/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_cbv-lc.fitsCOMPLETENoneNone
./mastDownload/HLSP/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_ens-lc.fitsCOMPLETENoneNone
./mastDownload/HLSP/hlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05/hlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05_cbv-lc.fitsCOMPLETENoneNone
./mastDownload/HLSP/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_lc.fitsCOMPLETENoneNone
./mastDownload/HLSP/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp.fitsCOMPLETENoneNone
./mastDownload/HLSP/hlsp_tglc_tess_ffi_gaiaid-5266270443442455040-s0001-cam4-ccd2_tess_v1_llc/hlsp_tglc_tess_ffi_gaiaid-5266270443442455040-s0001-cam4-ccd2_tess_v1_llc.fitsCOMPLETENoneNone

Plotting the light curves#

We will use the bokeh plotting library so that we can have interactivity, and will plot both the 2 minute and 30 minute cadence data together.

We can tell which file corresponds to which cadence length by examining the filenames and noting that one contains c0120 (2 minutes) and the other c1800 (30 minutes).

# Loading the short cadence light curve
hdu = fits.open(tasoc_manifest["Local Path"][0])
short_cad_lc = Table(hdu[1].data)
hdu.close()

# Loading the long cadence light curve
hdu = fits.open(tasoc_manifest["Local Path"][1])
long_cad_lc = Table(hdu[1].data)
hdu.close()
long_cad_lc
Table length=1281
TIMECADENCENOSAP_FLUXKSPSAP_FLUXKSPSAP_FLUX_ERRQUALITYORBITIDSAP_XSAP_YSAP_BKGSAP_BKG_ERRKSPSAP_FLUX_SMLKSPSAP_FLUX_LAG
float64int32float32float32float32int32int32float32float32float32float32float32float32
1325.323920580209246970.99075160.999379160.000472778740969833.8761609.046323.96545.070.99941080.99934494
1325.344753729780246980.990659530.999916550.000472778740969833.8771609.046292.76430.710.999935870.99988496
1325.365586880907846990.99092951.00057820.000472778740969833.8771609.046180.94480.81.00057021.0005482
1325.38642003356447000.99066241.00047530.000472778740969833.8771609.04518.43442.451.00047791.0004911
1325.407253187721447010.990464751.00023980.000472778740969833.8781609.04628.41461.861.00025391.0002934
1325.428086343321647020.99092641.00048590.000472778740969833.8781609.044-89.99462.611.00048591.0006546
1325.448919500346847030.991206941.0003830.000472778740969833.8781609.044-118.24472.781.00040151.0002135
1325.469752658739347040.99149771.0001450.000472778740969833.8791609.044-307.69311.051.0001511.0001317
1325.490585818471647050.992134331.00012830.000472778740969833.8791609.044-228.2359.261.00014141.0002065
1325.511418979517247060.9924970.99972690.000472778740969833.8791609.043-80.02372.230.999747930.99954337
.......................................
1352.969511727751360241.01276770.99988230.0004727787010833.8841609.154200.05784.490.999892831.0001168
1352.99034477458160251.01241620.99992710.0004727787010833.8841609.155181.46665.60.99994031.0000422
1353.011177821003560261.01180570.9997620.0004727787010833.8821609.153193.32571.460.99976290.9999271
1353.032010867052160271.0117021.0001450.0004727787010833.8841609.157258.93701.061.00013791.0001836
1353.052843912793860281.01110371.00008610.0004727787010833.8831609.153226.06637.771.00012041.0000738
1353.073676958258460291.01053931.00010720.0004727787010833.8831609.156267.1719.31.00009571.0001739
1353.094510003518560301.0117311.00191320.0004727787010833.8831609.156237.25576.671.001871.0020751
1353.115343048610460311.00921191.00009050.0004727787010833.8821609.156220.36776.461.00009391.0001388
1353.13617609361360321.00855951.00016250.0004727787409610833.8831609.156229.04772.091.00017051.0002179
1353.157009138568460331.00755950.999935030.0004727787409610833.8821609.157417.53742.140.999954340.99971807
long_cad_lc.columns
<TableColumns names=('TIME','CADENCENO','SAP_FLUX','KSPSAP_FLUX','KSPSAP_FLUX_ERR','QUALITY','ORBITID','SAP_X','SAP_Y','SAP_BKG','SAP_BKG_ERR','KSPSAP_FLUX_SML','KSPSAP_FLUX_LAG')>
bfig = plotting.figure(width=850, height=300, title=f"Detrended Lightcurve (TIC{tic_id})")

# Short cadence
bfig.circle(short_cad_lc["TIME"], short_cad_lc["RAW_FLUX"]/np.median(short_cad_lc["RAW_FLUX"]), fill_color="black",size=2, line_color=None)
bfig.line(short_cad_lc["TIME"], short_cad_lc["RAW_FLUX"]/np.median(short_cad_lc["RAW_FLUX"]), line_color='black')

# Long cadence
bfig.circle(long_cad_lc["TIME"],long_cad_lc["SAP_FLUX"], fill_color="#0384f7",size=6, line_color=None)
bfig.line(long_cad_lc["TIME"],long_cad_lc["SAP_FLUX"], line_color='#0384f7')

# Marking the flare (tpeak is in BJD, while the time column in the light curve is BTJD, so we must convert)
vline = Span(location=(tpeak - 2457000), dimension='height', line_color='#bf006e', line_width=3)
bfig.renderers.extend([vline])

# Labeling the axes
bfig.xaxis.axis_label = "Time (BTJD)"
bfig.yaxis.axis_label = "Flux"

plotting.show(bfig)

Experiment with the controls on the right to zoom in and and pan around on this light curve to look at the marked flare and other features (some sort of periodic feature maybe?)

Making an animation#

Looking at the above plot we can see the flare event in both the long and short cadence light curves. Since we can see it even in the 30 minute data, we should be able to make an animation of the area around the flaring star and see the flare happen.

We will use TESScut, the MAST cutout tool for full-frame images to cutout the area around the flaring star across the entire sector, and then make a movie that shows how it changes over time.

We will use the astroquery.mast Tesscut class to make this cutout.
We will use two functions:

  • Find the sky coordinate of our flare star: Observations._resolve_object*

  • Query for cutouts and get the result as a list of HDUList objects: Tesscut.get_cutouts **

* Observations._resolve_object is a private (not documented) function which is being removed in favor of the public function Observations.resolve_object in the next Astroquery release.

** We must start by finding the sky coordinate of our star, however starting with the next Astroquery release, Tesscut functions will be able to take an object name such as a TIC ID as well.

coord = Observations.resolve_object(f"TIC {tic_id}")

**Requesting a cutout target pixel file. **

This query will return a list of HDUList objects, each of which is the cutout target pixel file (TPF) for a single sector. In this case, because we specified a single sector we know that the resulting list will only have one element and can pull it out directly.

cutout_hdu = Tesscut.get_cutouts(coordinates=coord, size=40, sector=1)[0]
cutout_hdu.info()
Filename: <class '_io.BytesIO'>
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      57   ()      
  1  PIXELS        1 BinTableHDU    281   1282R x 12C   [D, E, J, 1600J, 1600E, 1600E, 1600E, 1600E, J, E, E, 38A]   
  2  APERTURE      1 ImageHDU        82   (40, 40)   int32   

A TESScut TPF has three extensions:

  • No. 0 (Primary): This HDU contains meta-data related to the entire file.

  • No. 1 (Pixels): This HDU contains a binary table that holds data like cutout image arrays and times. We will extract information from here to get our set of cutout images.

  • No. 2 (Aperture): This HDU contains the image extension with data collected from the aperture. We will use this extension to get the WCS associated with out cutout.

cutout_table = Table(cutout_hdu[1].data)
cutout_table.columns
<TableColumns names=('TIME','TIMECORR','CADENCENO','RAW_CNTS','FLUX','FLUX_ERR','FLUX_BKG','FLUX_BKG_ERR','QUALITY','POS_CORR1','POS_CORR2','FFI_FILE')>

The Pixels extension contains the cutout data table, which has a number of columns. For our puposes we care about the “TIME” column which has the observation times in BTJD, and the “FLUX” column which contains the cutout images (they are calibrated, but not background subtracted).

Exploring the cutout time series#

We want to explore what is happening within our cutout area over the time that the flare occurs, so we will make an animated plot of the cutout frames.

We can’t make a movie of the whole sector (it would take too long), so we will choose only the time range around the flare.

Use the light curve plot to figure out what time range you want to animate, or use our selections.

start_btjd = 1341.5
end_btjd = 1342.5

start_index = (np.abs(cutout_table['TIME'] - start_btjd)).argmin()
end_index = (np.abs(cutout_table['TIME'] - end_btjd)).argmin()

print(f"Frames {start_index}-{end_index} ({end_index-start_index} frames)")
Frames 721-769 (48 frames)

Looking at the animated cutout#

def make_animation(data_array, start_frame=0, end_frame=None, vmin=None, vmax=None, delay=50):
    """
    Function that takes an array where each frame is a 2D image array and make an animated plot
    that runs through the frames.
    
    Note: This can take a long time to run if you have a lot of frames.    
    Parameters
    ----------
    data_array : array
        Array of 2D images.
    start_frame : int
        The index of the initial frame to show. Default is the first frame.
    end_frame : int
        The index of the final frame to show. Default is the last frame.
    vmin : float
        Data range min for the colormap. Defaults to data minimum value.
    vmax : float
        Data range max for the colormap. Defaults to data maximum value.
    delay: 
        Delay before the next frame is shown in milliseconds.

    Returns
    -------
    response : `animation.FuncAnimation`
    """
    
    if not vmin:
        vmin = np.min(data_array)
    if not vmax:
        vmax = np.max(data_array)
        
    if not end_frame:
        end_frame = len(data_array) - 1 # set to the end of the array
        
    num_frames = end_frame - start_frame + 1 # include the end frame
        
    def animate(i, fig, ax, binarytab, start=0):
        """Function used to update the animation"""
        ax.set_title("Epoch #" + str(i+start))
        im = ax.imshow(binarytab[i+start], cmap=plt.cm.YlGnBu_r, vmin=vmin, vmax=vmax)
        return im,
    
    # Create initial plot.
    fig, ax = plt.subplots(figsize=(10,10))
    ax.imshow(data_array[start_frame], cmap=plt.cm.YlGnBu_r, vmin=vmin, vmax=vmax)

    ani = animation.FuncAnimation(fig, animate, fargs=(fig, ax, data_array, start_frame), frames=num_frames, 
                                  interval=delay, repeat_delay=1000)
    
    plt.close()
    
    return ani
make_animation(cutout_table['FLUX'], start_index, end_index, vmax=700, delay=150)

We can see three things in this plot:

  • The flare that occures in frames 740-743

  • An aberration that appears in frame 754

  • A variable star pulsing in the lower right corner


Exploring a variable star#

Now we will look more closely at the variable star we can see in the animation. We will use the TESS Input Catalog (TIC) to figure out the TESS ID of the variable star, and then using that ID pull down and explore the star’s light curve from MAST.

Overlaying TIC sources#

First we will use the astroquery.mast.Catalog class to query the TESS Input Catalog (TIC) and get a list of sources that appear in our cutout.

Here we do a simple cone search and the apply a magnitude limit. Play around with the magnitude limit to see how it changes the number of sources.

sources = Catalogs.query_object(catalog="TIC", objectname=f"TIC {tic_id}", radius=10*u.arcmin)
sources = sources[sources["Tmag"] < 12]
print(f"Number of sources: {len(sources)}")
print(sources)
Number of sources: 9
    ID           ra               dec        ... wdflag     dstArcSec     
--------- ---------------- ----------------- ... ------ ------------------
141914082 94.6175354047675 -72.0448462219924 ...      0                0.0
141914038 94.7353652417206 -72.0837957851839 ...      0   191.637404585262
141914103 94.7749592275334 -72.0227841314024 ...      0 192.00656069294254
141914130 94.5074642175418  -71.998930582098 ...      0  205.6248622381558
141914317 94.9032150266819 -72.0334742503164 ...      0   319.770063872556
141913994 94.5588022852622 -72.1321979465057 ...      0 321.11922083466345
166975135 94.6067354129356 -71.9255734581109 ...      0 429.55027273236567
141869504 94.2183733431791 -72.0232095513896 ...      0  450.0317520158321
141913929 94.7099836100023 -72.1924762405753 ...      0  541.2030879067767

Plotting sources on an individual cutout#

We will get the WCS infomation associated with our cutout from the Aperture extension, and use it to make a WCS-aware plot of a single cutout image. Then we display the image and sources together, and label the sources with their row number in the catalog table.

cutout_wcs = WCS(cutout_hdu[2].header)
cutout_img = cutout_table["FLUX"][start_index]
fig, ax = plt.subplots(subplot_kw={'projection':cutout_wcs})
fig.set_size_inches(10,10)
plt.grid(color='white', ls='solid')
    
# Setup WCS axes.
xcoords = ax.coords[0]
ycoords = ax.coords[1]
xcoords.set_major_formatter('d.ddd')
ycoords.set_major_formatter('d.ddd')
xcoords.set_axislabel("RA (deg)")
ycoords.set_axislabel("Dec (deg)")
ax.imshow(cutout_img, cmap=plt.cm.YlGnBu_r,vmin=0,vmax=700)
ax.plot(sources['ra'],sources['dec'],'x',transform=ax.get_transform('icrs'),color="red")

# Annotating the sources with their row nnumber in the sources table
for i,star in enumerate(sources):
    ax.text(star['ra']+0.01,star['dec'],i,transform=ax.get_transform('icrs'))

ax.set_xlim(0,cutout_img.shape[1]-1)
ax.set_ylim(cutout_img.shape[0]-1,0)

plt.show()
../../../_images/a60bdd2134c0a0889e6765bc93b7f7a78b4fdb45aa20fab67a11a5243b8bd8f4.png

The variable star is row 4 in the catalog sources table (note that if you changed the magnitude threshold the variale star may be in a different row).

sources["ID","ra","dec"][4]
Row index=4
IDradec
str11float64float64
14191431794.9032150266819-72.0334742503164

Getting the variable star light curve#

Again, we will look specifically for the TASOC light curve(s) associated with this star, rather than the mission pipeline ones. Below we go through the same process as in the Downloading the light curves section to search for the observation, then find the associated data products, and download them.

variable_tic_id = sources["ID"][4]

variable_res = Observations.query_criteria(target_name=str(variable_tic_id), 
                                           obs_collection="HLSP", 
                                           filters="TESS")

print(f"Number of tasoc light curves for {variable_tic_id}: {len(variable_res)}\n")

        
variable_prod = Observations.get_product_list(variable_res[0])
variable_manifest = Observations.download_products(variable_prod)
Number of tasoc light curves for 141914317: 80
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/qlp/s0013/0000/0001/4191/4317/hlsp_qlp_tess_ffi_s0013-0000000141914317_tess_v01_llc.fits to ./mastDownload/HLSP/hlsp_qlp_tess_ffi_s0013-0000000141914317_tess_v01_llc/hlsp_qlp_tess_ffi_s0013-0000000141914317_tess_v01_llc.fits ...
 [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HLSP/qlp/s0013/0000/0001/4191/4317/hlsp_qlp_tess_ffi_s0013-0000000141914317_tess_v01_llc.txt to ./mastDownload/HLSP/hlsp_qlp_tess_ffi_s0013-0000000141914317_tess_v01_llc/hlsp_qlp_tess_ffi_s0013-0000000141914317_tess_v01_llc.txt ...
 [Done]

Note that this time there is only one (30 minute cadence) TASOC light curve. This is because this star was not one of the targets that TESS observed at the shorter cadence.

hdu = fits.open(variable_manifest["Local Path"][0])
variable_lc = Table(hdu[1].data)
hdu.close()

Plotting the variable star light curve#

We will again plot the light curve using bokeh so that we have interactivity.

variable_lc
Table length=1320
TIMECADENCENOSAP_FLUXKSPSAP_FLUXKSPSAP_FLUX_ERRQUALITYORBITIDSAP_XSAP_YSAP_BKGSAP_BKG_ERRKSPSAP_FLUX_SMLKSPSAP_FLUX_LAG
float64int32float32float32float32int32int32float32float32float32float32float32float32
1653.9281264088759204701.02684971.05587150.13665526412833421.720151216.579255.87303.01.04919471.0580579
1653.948959649118204710.93603960.96120990.13665526409633421.724241216.578463.58413.660.966199040.95966774
1653.969792891784204720.82098440.84196810.13665526409633421.731871216.5752-218.77387.510.860418140.83654326
1653.9906261366607204730.76165690.78014080.13665526409633421.7351216.5725-0.87454.050.806358340.7720926
1654.011459383893204740.82817890.84724130.13665526409633421.730221216.5757-157.54447.660.86502090.84183574
1654.0322926332728204750.943153260.963719840.13665526409633421.72361216.579784.09448.030.968248250.96247375
1654.053125884903204761.02805531.04926750.13665526409633421.719151216.5809-48.59350.111.04341131.0512749
1654.0739591385864204771.08233211.10343490.13665526409633421.718631216.5817-22.17435.481.09144061.1074061
1654.0947923943963204781.1055961.12593790.13665526409633421.71661216.58342.58475.551.11168271.1304736
1654.1156256521465204791.09297621.11192580.13665526409633421.716131216.581815.22380.431.09921151.1160737
.......................................
1682.1571617084755218250.87786420.85018370.13665526409634421.727871216.628-97.16516.830.86614540.84522027
1682.1779948309704218260.70115070.67423320.13665526413234421.725431216.6162-2064.341480.940.704956050.66895413
1682.1988279531406218270.78981720.760783430.13665526034421.733251216.6271-102.71397.850.78698680.7528833
1682.2196610751407218280.892728570.85744180.13665526034421.72581216.630920.56437.690.87308710.8527917
1682.2404941969792218290.99515120.952960430.13665526034421.72071216.6326223.22546.280.958030460.95133626
1682.2613273188597218301.05914041.01108620.13665526034421.71721216.633774.16340.841.01001051.0116068
1682.2821604407818218311.10349711.05003270.13665526034421.716641216.6342250.13395.141.04491811.0516453
1682.3029935630043218321.11160961.05421270.13665526034421.716131216.6346306.66411.171.04881791.0560246
1682.3238266855037218331.08300141.02351740.13665526034421.716981216.6335295.34449.061.02172271.0241388
1682.3446598085977218341.02434480.964601460.13665526409634421.7191216.6324146.27408.530.96875610.96340925
bfig = plotting.figure(width=850, height=300, title=f"Detrended Lightcurve (TIC{variable_tic_id})")

bfig.circle(variable_lc["TIME"],variable_lc["SAP_FLUX"], fill_color="black",size=4, line_color=None)
bfig.line(variable_lc["TIME"],variable_lc["SAP_FLUX"], line_color='black')

# Labeling the axes
bfig.xaxis.axis_label = "Time (BTJD)"
bfig.yaxis.axis_label = "SAP Flux"

plotting.show(bfig)

That looks variable all right!

Finding the period#

We’ll run a Lomb Scargle priodogram on this light curve to see if we can quantify the periodic behavior. To do this we will use the astropy.timeseries class LombScargle.

lomb = LombScargle(variable_lc["TIME"], variable_lc["SAP_FLUX"])
frequency, power = lomb.autopower(maximum_frequency=25)

Plotting the periodogram#

bfig = plotting.figure(width=850, height=300, x_range=(0,25),
                       title=f"Periodogram (TIC{variable_tic_id})")

bfig.line(frequency, power, line_color='black')

# Labeling the axes
bfig.xaxis.axis_label = "Frequency (1/day)"
bfig.yaxis.axis_label = "Power"

plotting.show(bfig)

Phasing on the highest power period/frequency#

There is a clear dominant frequency in the above plot, with a small harmonic also visible. We will phase the stellar light curve on the period corresponding to the dominant frequency and plot both it and the corresponding sinusoidal fit.

dominant_freq = frequency[np.argmax(power)].value
print(f"The dominant period: {1/dominant_freq*24:.3} hours")
The dominant period: 5.24 hours
bfig = plotting.figure(width=850, height=300, title=f"Phased Lightcurve (TIC{variable_tic_id})")

# Plotting the phased light curve
bfig.circle(variable_lc["TIME"]%(1/dominant_freq),variable_lc["SAP_FLUX"], fill_color="black",size=4, line_color=None)

# Plotting the periodic fit
t_fit = np.linspace(0,1/dominant_freq,100)
bfig.line(t_fit, lomb.model(t_fit, dominant_freq), color='#1b9f00', line_width=2)

# Labeling the axes
bfig.xaxis.axis_label = "Phase (days)"
bfig.yaxis.axis_label = "Flux"

plotting.show(bfig)

Additional Resources#

About this Notebook#

Author: C. E. Brasseur, STScI Software Engineer
Last Updated: May 2023

Top of Page STScI logo