Intermediate: Finding Flares and Variable Stars in TASOC Light Curves

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.io import fits
from astropy.wcs import WCS
from astropy.timeseries import LombScargle
import astropy.units as u

# For matplotlib plotting
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# For animation display
from matplotlib import rc

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",
    provenance_name="TASOC",
    sequence_number=sector,
)
tasoc_res[
    "dataproduct_type",
    "obs_collection",
    "target_name",
    "t_exptime",
    "filters",
    "provenance_name",
    "project",
    "sequence_number",
    "instrument_name",
]
Table length=2
dataproduct_typeobs_collectiontarget_namet_exptimefiltersprovenance_nameprojectsequence_numberinstrument_name
str10str4str9float64str4str5str4int64str10
timeseriesHLSP141914082120.0TESSTASOCTESS1Photometer
timeseriesHLSP1419140821800.0TESSTASOCTESS1Photometer

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=3
dataproduct_typedescriptiondataURIsize
str10str4str125int64
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/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
tasoc_prod
Table masked=True length=3
obsIDobs_collectiondataproduct_typeobs_iddescriptiontypedataURIproductTypeproductGroupDescriptionproductSubGroupDescriptionproductDocumentationURLprojectprvversionproposal_idproductFilenamesizeparent_obsiddataRightscalib_levelfilters
str8str4str10str65str4str1str125str7str28str1str1str5str1str23str77int64str8str6int64str4
79352767HLSPtimeserieshlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05FITSCmast:HLSP/tasoc/s0001/c0120/0000/0001/4191/4082/hlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05_cbv-lc.fitsSCIENCEMinimum Recommended Products----TASOC5G011160_G011155_G011188hlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05_cbv-lc.fits188064079352767PUBLIC4TESS
79500562HLSPtimeserieshlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05FITSCmast:HLSP/tasoc/s0001/c1800/0000/0001/4191/4082/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_cbv-lc.fitsSCIENCEMinimum Recommended Products----TASOC5G011160_G011155_G011188hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_cbv-lc.fits16416079500562PUBLIC4TESS
79500562HLSPtimeserieshlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05FITSCmast:HLSP/tasoc/s0001/c1800/0000/0001/4191/4082/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_ens-lc.fitsSCIENCEMinimum Recommended Products----TASOC5G011160_G011155_G011188hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_ens-lc.fits16704079500562PUBLIC4TESS

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

INFO: Found cached file ./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 with expected size 1880640. [astroquery.query]
INFO: Found cached file ./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 with expected size 164160. [astroquery.query]

INFO: Found cached file ./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 with expected size 167040. [astroquery.query]
Table length=3
Local PathStatusMessageURL
str163str8objectobject
./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_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

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

short_cad_lc
Table length=19261
TIMETIMECORRCADENCENOFLUX_RAWFLUX_RAW_ERRFLUX_BKGFLUX_CORRFLUX_CORR_ERRQUALITYPIXEL_QUALITYMOM_CENTR1MOM_CENTR2POS_CORR1POS_CORR2
float64float32int32float64float64float64float64float64int32int32float64float64float64float64
1325.29473191055150.000392398777044458018.6914062527.0715332031255857.234375nannan18834.48949134566251609.6587898246nannan
1325.29612078712130.000392386957044558050.1210937527.0800399780273445833.82373046875-15477.966097362183459.2737368157305600834.4306222158931609.59601199911530.0649784728884697-0.03711783513426781
1325.29750966369060.000392375147044657973.3398437527.0668315887451175845.93603515625-16172.59679890426459.333388530095700834.40326775202361609.56349750617460.036016035825014114-0.07080704718828201
1325.2988985402590.000392363327044758019.1210937527.081874847412115858.75341796875-14821.928321840816459.8564875964108300834.4024918155881609.56020550229410.03578070178627968-0.0754164531826973
1325.30028741682870.00039235157044858076.8789062527.0953884124755865857.07421875-13300.352272517579460.3382741139668600834.40166899684981609.54902188860230.03445006534457207-0.08668677508831024
1325.30167629342710.000392339727044958152.54687527.105075836181645852.03125-11505.904315655347460.740051243747900834.39528263329881609.54893874850150.027773790061473846-0.0863100215792656
1325.30306517002460.000392327937045058128.6289062527.1025371551513675861.92041015625-11435.558188753192460.919258179338200834.39424804777581609.54848397795040.026579076424241066-0.08645732700824738
1325.30445404662350.000392316147045158119.804687527.10265541076665862.89794921875-11139.658936588527461.12923498827400834.39877300482511609.54514139779580.030768081545829773-0.09115128219127655
1325.3058429232510.000392304387045258085.26562527.0967693328857425874.02685546875-11311.563370110078461.2230351568018700834.3967953924461609.54730108034280.029517315328121185-0.0876065194606781
..........................................
1353.16393257751320.000159770289051058783.8476562527.283060073852546084.0654296875-2344.7684472455153463.0368494184998600834.3955345686121609.6570044969480.031860228627920150.02544116973876953
1353.16532144717030.000159751689051158778.3164062527.282293319702156090.095703125-2821.0573842525346462.846336336718700834.39918711660771609.65585711531840.035427540540695190.023608732968568802
1353.1667103168270.000159733089051258796.89062527.2799434661865236065.37109375-2920.2006968637174462.6142686650985600834.39861024588971609.6576971215560.034142192453145980.02651386149227619
1353.16809918648320.000159714499051358739.0507812527.2695960998535166071.37109375-4347.702717111712462.231439590482500834.4026523813731609.65414117579580.0381550006568431850.02153876982629299
1353.1694880561250.000159695889051458791.304687527.2804241180419926067.6552734375-3942.961521480437462.191791794548900834.39487756831511609.63782202758920.029868377372622490.005358986556529999
1353.17087692578160.000159677289051558736.6367187527.2757854461669926106.0126953125-5384.693659230999461.8738016466022600834.39719181223641609.65310691793930.033532425761222840.020877502858638763
1353.17226579543850.000159658689051658750.9492187527.2727432250976566067.3798828125-5694.076706754814461.5661618028114500834.40691603675091609.65729299670470.043107997626066210.0256554763764143
1353.17365466509550.000159640089051758763.9023437527.276618957519536087.22607421875-6063.881129877546461.3583458760828600834.40080296929271609.6681374877340.0373148210346698760.037226371467113495
1353.17504353475170.000159621499051858787.7695312527.281454086303716082.94140625-6287.681655247179461.148929515107500834.39268278877651609.63919479225660.027845425531268120.006460509728640318
1353.17643240440840.000159602899051958860.570312527.2950668334960946069.53515625-5724.686801500067461.069455843939500834.40145103540661609.66667314225650.0375200696289539340.034000542014837265
long_cad_lc.columns

<TableColumns names=('TIME','TIMECORR','CADENCENO','FLUX_RAW','FLUX_RAW_ERR','FLUX_BKG','FLUX_CORR','FLUX_CORR_ERR','QUALITY','PIXEL_QUALITY','MOM_CENTR1','MOM_CENTR2','POS_CORR1','POS_CORR2')>

bfig = plotting.figure(
    width=850, height=300, title=f"Detrended Lightcurve (TIC{tic_id})"
)

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

# Long cadence
bfig.scatter(
    long_cad_lc["TIME"],
    long_cad_lc["FLUX_RAW"],
    fill_color="#0384f7",
    size=6,
    line_color=None,
)
bfig.line(long_cad_lc["TIME"], long_cad_lc["FLUX_RAW"], 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. In order to do this, we’ll need a way to get all of the images for this target. For our object of interest, we could use the pre-existing target pixel files (TPFs); however, many interesting targets do not have TPFs from TESS. Instead, we’ll 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.

There is a TESScut API: the astroquery.mast Tesscut class, which we can use to make this cutout. Specifically, we can call the Tesscut.get_cutouts function.

cutout_hdu = Tesscut.get_cutouts(objectname=f"TIC {tic_id}", size=40, sector=1)[0]

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.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/635e3b534614660997fffd38312a8b89ce4e9645903bae83108c0a872a13975d.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: 102

INFO: Found cached file ./mastDownload/HLSP/hlsp_qlp_tess_ffi_s0013-0000000141914317_tess_v01_llc/hlsp_qlp_tess_ffi_s0013-0000000141914317_tess_v01_llc.fits with expected size 83520. [astroquery.query]

INFO: Found cached file ./mastDownload/HLSP/hlsp_qlp_tess_ffi_s0013-0000000141914317_tess_v01_llc/hlsp_qlp_tess_ffi_s0013-0000000141914317_tess_v01_llc.txt with expected size 64163. [astroquery.query]

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
.......................................
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.scatter(
    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.scatter(
    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

Top of Page STScI logo

Contents