# Other Fun Features¶

Some other features that are included in the $$\texttt{stella}$$ package include finding rotation periods and fitting flares with a simple analytic model to extract parameters such as the equivalent duration. Here, we will go through each of these additional modules.

[2]:

import os, sys
sys.path.insert(1, '/Users/arcticfox/Documents/GitHub/stella/')
import stella
import numpy as np
from tqdm import tqdm_notebook
import matplotlib.pyplot as plt

plt.rcParams['font.size'] = 20


First thing’s first: we need a light curve! We’re going to use the same one that has been used in the previous demonstrations.

[3]:

from lightkurve.search import search_lightcurvefile
lc = search_lightcurvefile(target='tic62124646', mission='TESS')
lc.plot()
lc = lc.remove_nans()

//anaconda3/lib/python3.7/site-packages/lightkurve/lightcurvefile.py:47: LightkurveWarning: LightCurveFile.header is deprecated, please use LightCurveFile.get_header() instead.
LightkurveWarning)


## 1.1 Measuring Rotation Periods¶

To measure the rotation period of a light curve, you need the following information: a target ID (can be anything you want), a time array, a flux array, and a flux error array. These can all be retrieved from the $$\texttt{lightkurve}$$ object. Then, you initialize the class $$\texttt{stella.MeasureProt}$$.

[4]:

mProt = stella.MeasureProt([lc.targetid], [lc.time], [lc.flux], [lc.flux_err])


The rotation measurement tool used in this class is a Lomb-Scargle Periodogram. You have the option the set the minimum frequency (minf=1/12.5) and maximum frequency (maxf=1/0.1) that the periodogram searches over. You can also set the samples per peak (spp=50). The default values are noted in parentheses.

[5]:

mProt.run_LS()

Finding most likely periods: 100%|██████████| 1/1 [00:00<00:00, 29.80it/s]


Calling this will create an $$\texttt{astropy.table.Table}$$ attribute with some metrics about the rotation period measured in that light curve. The columns in this table include: - Target_ID: the ID for that light curve. - period_days: the period as measured in the first orbit. - secondary_period_days: the period as measured in the second orbit. - gauss_width: the width of a best-fit Gaussian to the most likely period. - max_power: the power of the periodogram at the most likely period. - orbit_flag: a combined flag per orbit (0 = reliable measurement). - oflag1: period flag for the first orbit (0 = reliable measurement). - oflag2: period flag for the first orbit (0 = reliable measurement). - Flags: Defines a flag based on all observations of the given target. This mainly applies to stars observed in multiple sectors, otherwise it is the same value as the orbit_flag (0 = reliable measurement). - avg_period_days: the average most likely period for that target. This is averaged over multiple sectors when available.

Note that by fitting to each orbit, we are limiting the period search space to relatively short rotation periods.

[6]:

mProt.LS_results

[6]:

Table length=1
Target_IDperiod_dayssecondary_period_daysgauss_widthmax_powersecondary_max_powerorbit_flagoflag1oflag2Flagsavg_period_days
int64float64float64float64float64float64float64float64float64int64float64
621246463.22967914764591773.21636224711885930.33188306480480350.948139023632630.19343617519244910.00.00.003.2296791476459177

For this star, we find a period of 3.23 days. We can fold over this period in $$\texttt{lightkurve}$$ to see if this period is found in the light curve:

[8]:

lc.fold(mProt.LS_results['avg_period_days'].data[0]).plot()

[8]:

<matplotlib.axes._subplots.AxesSubplot at 0x14ca0c400>


That looks nice! (And so do those flares!!)

## 1.2 Fitting Flares¶

In $$\texttt{stella}$$, we use a very basic model flare fit of a sharp rise and an exponential decay. This can be done through a class included in the code. The first thing that needs to be done is get predictions for where the flares occur.

[9]:

MODEL_DIR = '/Users/arcticfox/Documents/flares/run01/'
MODEL = [os.path.join(MODEL_DIR,i) for i in
os.listdir(MODEL_DIR) if i.endswith('.h5')][0]

[10]:

cnn = stella.ConvNN(output_dir='.')
cnn.predict(modelname=MODEL,
times=lc.time,
fluxes=lc.flux,
errs=lc.flux_err)

WARNING: No stella.DataSet object passed in.
Can only use stella.ConvNN.predict().

100%|██████████| 1/1 [00:00<00:00,  1.16it/s]


Again, now we have a light curve with predictions from just one of the $$\texttt{stella}$$ models.

[11]:

plt.figure(figsize=(14,4))
plt.scatter(cnn.predict_time[0], cnn.predict_flux[0], c=cnn.predictions[0], vmin=0, vmax=1)
plt.xlabel('Time [BJD - 2457000]')
plt.ylabel('Normalized Flux');


Now we can initiate our flare fitting class, stella.FitFlares, and start fitting those flares!

[12]:

ff = stella.FitFlares(id=[lc.targetid],
time=[lc.time],
flux=[lc.flux],
flux_err=[lc.flux_err],
predictions=[cnn.predictions[0]])


The flares are identified given a certain probability threshold (here, I use a threshold of 0.5, which is what was used in Feinstein et al. (2020)). For the cleanest sample of flares, one could explore a higher probability threshold.

Several parameters about the flare are returned in an astropy Table: - tpeak: the peak time of the flare - amp: the amplitude of the flare - dur_min: the duration of the flare, given in minutes - rise: how quickly the flare goes off (this is very fast and generally the same across all flares) - fall: the exponential decay - prob: the probability at the tpeak time

[13]:

ff.identify_flare_peaks(threshold=0.5)

ff.flare_table

Finding & Fitting Flares: 100%|██████████| 1/1 [00:00<00:00,  3.06it/s]

[13]:

Table length=16
Target_IDtpeakampdur_minrisefallprob
float64float64float64float64float64float64float64
62124646.01656.45766947437780.00517271866554260139.0307026803608960.00010.00688458976865427650.9807868599891663
62124646.01657.31324338198560.00935070228511011346.0427474926023450.00010.005440531357077170.9993744492530823
62124646.01659.18828153893080.0165110534040283139.068943053366020.00010.0048665709253705070.9999998807907104
62124646.01659.93274023421740.00531299853780956839.0434203414553450.00010.0099615092212989240.6145604252815247
62124646.01661.35915489099530.01616045832644144640.07148748342360.00010.0054786545406301720.9998658895492554
62124646.01664.60226137635350.0192907510039818140.1231279321261950.00010.0083774666962270731.0
62124646.01665.10504612819890.01314713709421687539.0896531890076840.00010.0080857918850959530.9999781847000122
62124646.01671.53841797829930.02035861226105141839.1231058612179940.00010.0077604368901831541.0
62124646.01673.18702736909740.0205590406398600839.114627785012330.00010.007054779840450441.0
62124646.01676.2536836973380.03055554684541682339.214232149922520.00010.0090138280819179881.0
62124646.01676.95506884504040.00673877124922315639.025499944224180.00010.0046305339257480510.997788667678833
62124646.01677.32451111334220.2810143935238597640.429627767115780.00010.0065647640076292621.0
62124646.01677.79534147014150.00859558528051674539.04217474905530.00010.00624890500430893250.9562909603118896
62124646.01679.16477600030750.006115484082584130539.035647558860110.00010.00713760243215867150.9747272729873657
62124646.01679.5717171545920.00905349498155829339.06350144752650.00010.0088848388797201030.9999923706054688
62124646.01681.73419719783940.01064266441690119239.0321765605989060.00010.0034022541541514470.9999966621398926

Cool! We can mark all of our flares as well, just to see them better by-eye:

[14]:

plt.figure(figsize=(14,4))
plt.scatter(ff.time[0], ff.flux[0], c=cnn.predictions[0], s=5)

for tpeak in ff.flare_table['tpeak']:
plt.vlines(tpeak, 0,2, color='k', alpha=0.5, linewidth=5, zorder=0)

plt.ylim(0.94,1.3)
plt.xlabel('Time [BJD - 2457000]')
plt.ylabel('Normalized Flux');

[15]:

plt.figure(figsize=(14,4))
plt.scatter(ff.time[0], ff.flux[0], c=cnn.predictions[0], s=5)
for tpeak in ff.flare_table['tpeak']:
plt.vlines(tpeak, 0,2, color='k', alpha=0.5, linewidth=5, zorder=0)
plt.xlim(1660,1666)
plt.ylim(0.96,1.05)

[15]:

(0.96, 1.05)


You’ll also note there is a flare at time$:nbsphinx-math:sim$1665.5 that was not identified. As there was a significant gap in the data and the CNN cannot handle data gap, the cadences around this region were ignored. Hence, that flare was not identified.

Additionally, the light curve looks to have a few flares at time$:nbsphinx-math:sim$1663. Zooming in, these don’t look like flares. This is better seen when the eye is not guided by the prediction assignment of the neural network:

[16]:

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(14,4), sharex=True, sharey=True)
ax1.scatter(ff.time[0], ff.flux[0], c=cnn.predictions[0], s=8)
ax2.plot(ff.time[0], ff.flux[0], 'k.')
for tpeak in ff.flare_table['tpeak']:
plt.vlines(tpeak, 0,2, color='k', alpha=0.5, linewidth=5, zorder=0)
plt.xlim(1662.5,1663.5)
plt.ylim(1.01,1.035)
plt.xticks(np.arange(1662.5, 1664.0, 0.5))
ax1.set_ylabel('Normalized Flux')
ax1.set_xlabel('Time [BJD - 2457000]', x=1.1)

[16]:

Text(1.1, 0, 'Time [BJD - 2457000]')