% !TEX root = ../thesis-example.tex
%
\chapter{Data and Methods}
\label{sec:data}

In this chapter the selection criteria for the data are explained (section \ref{sec:data:data_selection}). Furthermore a detailed description of the algorithm is given in section \ref{sec:data:data_reduction}.

\section{Data selection}
\label{sec:data:data_selection}

This study uses a broad selection of Kepler/K2 and TESS lightcurves (427 for Kepler, 27 for K2, 1579 for TESS), which were downloaded from MAST with the help of the astroquery (\cite{astroquery}) python package. For Kepler and its continuation mission K2 short-cadence data were used. Unlike long-cadence data with a cadence of 30 minutes, this allows the resolution of shorter events too, as the duration of flares can vary between a few tens of seconds to a few hours (\cite{flare_duration1}).
The initial dataset was taken from a list of well known flaring stars from \cite{kepler_411_study}, \cite{kepler_411_210_comparison}, \cite{doyle_2018}, \cite{doyle_2019}, \cite{au_mic_flaring_spi} and \cite{flare_occurance_periodicity}. These stars were selected because previous studies revealed that these stars showed detectable flares. This allowed a test of the algorithm described in the next section and a comparison with previous results.
Additionally the dataset of M to F stars of \cite{althukair_starlist} which can be found at \href{https://github.com/akthukair/AFD}{https://github.com/akthukair/AFD} was parsed to the downloader GUI. As not all data were available as short cadence data from Kepler/K2, only a subset of this large sample was added. The full list of stars, for which datasets were downloaded and parsed can be found in appendix \ref{chap:list_of_stars} table \ref{apA:list_of_all_stars}.

\section{Data reduction algorithms}
\label{sec:data:data_reduction}

This chapter explains the methods used in this study, split into the algorithms for flare detection and folding lightcurves. Everything, like the GUI discussed in chapter \ref{sec:gui}, was written in python 3 (\cite{10.5555/1593511}). It also makes extensive use of the python packages astropy (\cite{astropy:2018}), numpy (\cite{numpy}), scipy (\cite{scipy}), pandas (\cite{pandas}) and lightkurve (\cite{lightkurve}). The focus of this chapter lays in the description of the method used to create the final output, which is less customizable than the GUI, which exposes most parameters offered in the functions of the lightkurve API. Unless stated otherwise, the default parameters are used. Furthermore for the final output only PDCSAP\_FLUX is used, which is set as the default flux in the lightkurve $LightCurve$ objects after reading the fit file. The methods are the same for both $KeplerLightCurve$ and $TessLightCurve$ subclasses.

\subsection{Flare detection}
\label{sec:data:data_reduction:flare_detection}

The first step is to normalize the lightcurve. This is done to apply the same thresholds to all files in later steps. An example of this is shown in figure \ref{fig:full_gui_normal_selection_normalize_options} in chapter \ref{sec:gui:data_display}. The normalization is done via the $normalize()$ function of the $LightCurve$ class of the lightkurve api. Afterwards the lightcurves are flattened by the function $flatten()$ of the $LightCurve$ objects. The resulting object is then used as the base for the detection of flares. This removes all longterm trends like brightness changes due to spot modulation or similar, while retaining short term events like flares or transits. A similar approach was used by \cite{au_mic_flaring_spi}.\\
The next step is then to call $calculateFlareFitsForLightcurve()$ with the flattened lightcurve as well as the normalized lightcurve as parameters. It returns two lists of python dictionaries with the data for the flare peak as well as a fit which is described in the following paragraphs.
It parses the flattened lightcurve with the scipy $find\_peaks$ function. This function returns local maxima, which can be further filtered by their minimum height (optional) as well as the minimum distance of datapoints (optional). Afterwards the found peaks are sorted by height, and the highest 100 are returned. This was found to be a good amount as the most flares per fits file found in this study were around 70 for CD-56 1032A and B.\\
Afterwards each individual peak is checked. The first step is to estimate a flare starting and end point. This is done by checking the datapoints before and after the peak. If it finds that the flux delta is below 0.005 for three consecutive datapoints, it stops, and assumes that the last checked point is the start/end of the flare. In the case it finds an infinite or NaN value (which can happen if there are gaps in the lightcurve data), or it reaches 100 datapoints before/after it will stop. This was found to cover most flares detected and provides a sufficient number of datapoints for the following steps.\\
Afterwards multiple checks are done. It first checks if the datapoints right before and after the peak are above a threshold of 0.003 (0.3\%) above the mean flattened flux, or if the two datapoints right after the peak are above that threshold (which is a similar approach to \cite{kepler_411_study}). Afterwards it is checked if the datapoint at two indices before the peak is larger than the datapoint right before the registered peak. This eliminates possible false positives which were found by manually inspecting previously found events, but could not be guaranteed to be two seperate events. In such cases it will only count one flare. Then the flare lightcurve is fitted. For the impulsive phase, a gaussian function is used and for the gradual phase an exponetial function is used to account for the exponential decay in flare lightcurves (similar approach to \cite{au_mic_flaring_spi} and \cite{doyle_2018}). Then the residual sum of squares (RSS) between the fit and the flux of the flare, as well as the total sum of squares (TSS) are calculated. Afterwards R-squared is calculated, and if it is below 0.8, the flare is rejected as the flare would not have the typical form. In the last step, the locations of the flare in the normalized and flattened lightcurve are compared. This step has been introduced, as in some rare cases the flattening algorithm can produce a large artificial spike (values of 10 or higher when normalized).

\subsection{Lightcurve folding}
\label{sec:data:data_reduction:lightcurve_folding}

This section mainly describes how the $getOptimizedFold()$ function works. This function was written to optimize the lightcurve folding, to improve accuracy for the folding and lightcurve fitting process. It uses the normalized lightcurve as well as the preferred fitting function type (see chapter \ref{sec:gui:data_display}) as parameters. The fit type can either be "sine" for a sine fit, "poly" for a polynomlial fit, or "linear" for a linear fit. The default value is "sine", but can be changed for each individual star in the GUI.
The function at first generates two periodograms with the lightkurve function $to\_periodogram$. The first one uses the lombscargle algorithm, while the second one uses the box-least-squares algorithm. Afterwards the four highest peaks of each are taken and converted into periods (unit in days).
Then each period found by each algorithm is compared with the periods found by the other, and in the case of a match (absolute value of the difference between the values of the two algorithms is smaller than 5\% of the larger of the two periods) this is now used as the rotational period as well as the spot modulation period. If it does not find a match, it uses the period corresponding to the highest peak in the lombscargle periodogram.
Afterwards it calculates the epoch time to the first minimum in the lightcurve using numpys $argrelextrema$ function. In the following loop, which is repeated up to 30 times, the lightcurve is folded (using the lightkurve $fold$ function) with the spot modulation period, and the epoch set to the first minima in the lightcurve. Then a fit is performed using the prefered fitting method. If the fitting process fails for the preferred fitting type, and it was set to "sine", it will fall back to "poly" (polynomial fit), and vice versa. If the resulting fit has two maxima, and they are further away from the edges than 10\% of the used period, the algorithm checks if there is a signal for half the used period in any of the two periodograms. If this is the case, it will use this as the period for spot modulation. It will save the current folded lightcurve as folded by rotational period, and generate a new folded lightcurve with the newly found spot modulation period. Afterwards it checks the location of the minimum of the spot modulation folded lightcurve. If the minimum is within 1\% of half the phase it will stop. If not, it will shift the epoch by the difference of the minimum to zero, which is the center of the folded lightcurve phase, and repeat the folding process.\\
The function then returns the rotational period, the spot modulation period, the periods found by both periodograms as well as the folded lightcurve, the corresponding phase and the fit and fit type ("sine", "poly" or "linear") as well as the epoch. If the spot modulation differs from the found rotational period, it will also return the folded lightcurve, phase, fit and fit type of the rotational period folded lightcurve. The last returned value is "isValid", which is set to false if the algorithm does not find a fit with its minimum at half the phase for the spot modulation period folded lightcurve within 30 iterations.

\section{Spectral type identification}
\label{sec:data:sptype_identification}

After the initial download and cross checking with SIMBAD, not all stars had a spectral type assigned. For these stars magnitudes in different wavelength bands were available. The stars with known magnitudes at least in the B and V bands are given in table \ref{tab:unknown_sptypes_bv}, while the stars with only J, H and K bands are given in table \ref{tab:unknown_sptypes_jhk}.

The B-V color index as well as other color indices (if available) have been compared to the table in \cite{spectral_type_color_table} (\href{https://www.stsci.edu/~inr/intrins.html}{https://www.stsci.edu/~inr/intrins.html}). 

This allowed a rough classification of the stars in table \ref{tab:unknown_sptypes_bv}. Additionally, also the effective temperature was taken into account from the TESS Input Catalogue (\cite{revised_tess_input_catalogue}) from Vizier (\cite{vizier}). The classification table for spectral types based on effective temperature has been taken from \cite{harvard_spectral_types_teff} (\href{https://lweb.cfa.harvard.edu/~pberlind/atlas/htmls/note.html}{https://lweb.cfa.harvard.edu/~pberlind/atlas/htmls/note.html}).\\
The final spectral type for the stars can be found in table \ref{tab:unknown_sptypes_final_output}.

If the spectral type calculated from the effective temperature and the indices is identical, it will use that spectral type as final.

In the case of those differing, like for BD-08 995, which was found to be an early K type star (K0) by the intrinsic color method but being a G type star based on the effective temperature, the latter is prefered. For the four stars from table \ref{tab:unknown_sptypes_jhk} which have no B or V measurements, the spectral type based on the effective temperature is used.

\begin{table}
    \centering
    \caption{List of stars with no spectral type entry in SIMBAD, but available magnitudes in the B and V bands. Main Identifier according to SIMBADs "MAIN\_ID" property, as well as TESS Input Catalogue and Kepler Input Catalogue numbers are given.}
    \label{tab:unknown_sptypes_bv}
    \begin{tabular}{lll}
    \hline
    Main Identifier & TIC & KIC \\
    \hline\hline
    1RXS J064643.6-770027 & TIC 177255827 & - \\
    BD-08   995 & TIC 43472154 & - \\
    CPD-19   878 & TIC 93125144 & - \\
    PM J07058-5848 & TIC 279615427 & - \\
    TYC 1360-957-1 & TIC 247117382 & - \\
    TYC 4595-107-1 & TIC 394030788 & - \\
    \hline
    \end{tabular}
\end{table}


\begin{table}
    \centering
    \caption{List of stars with no spectral type entry in SIMBAD, but available magnitudes in the J, H and K bands. Main Identifier according to SIMBADs "MAIN\_ID" property (2MASS), as well as TESS Input Catalogue and Kepler Input Catalogue numbers are given.}
    \label{tab:unknown_sptypes_jhk}
    \begin{tabular}{lll}
    \hline
    Main Identifier & TIC & KIC \\
    \hline\hline
    2MASS J18524052+4156057 & - & KIC 6500181 \\
    2MASS J18534407+4208274 & - & KIC 6668646 \\
    2MASS J18535462+4135227 & - & KIC 6183672 \\
    2MASS J19335656+4010546 & - & KIC 5016904 \\
    \hline
    \end{tabular}
\end{table}

\begin{landscape}
\begin{table}
    \centering
    \caption{Final determined spectral types. Main Identifier taken from SIMBAD. The column B-V contains the value of B-V calculated from the respective magnitudes taken from SIMBAD. The columns T\textsubscript{eff}, MASS and Radius contain the effective temperature, mass and radius for the star taken from the TESS Input Catalogue from Vizier. The column Spectral Type (B-V) contains the determined spectral type based on \cite{spectral_type_color_table}, while the column Spectral Type (T\textsubscript{eff}) contains the spectral type based on the effective temperature from \cite{harvard_spectral_types_teff}. The last column contains the spectral type which is used in this study.}
    \label{tab:unknown_sptypes_final_output}
    \begin{tabular}{lccccccc}
    \hline
    Main Identifier & B-V & T\textsubscript{eff} & Mass & Radius & Spectral Type (B-V) & Spectral Type (T\textsubscript{eff}) & Final Spectral Type \\
     & & [$K$] & [$M_\odot$] & [$R_\odot$ ]& & & \\
    \hline\hline
    1RXS J064643.6-770027 & 1.320000 & 4082 & 0.640 & 0.654 & K7.0 & K & K \\
    BD-08   995 & 0.820000 & 5316 & 0.920 & 0.905 & K0.0 & G & G \\
    CPD-19   878 & 1.087001 & 4852 & 0.790 & 0.893 & K5.0 & K & K \\
    PM J07058-5848 & 1.320001 & 4250 & - & 0.409 & K7.0 & K & K \\
    TYC 1360-957-1 & 0.970000 & 4870 & 0.790 & 0.636 & K3.0 & K & K \\
    TYC 4595-107-1 & 0.750000 & 5231 & 0.890 & 0.902 & G8.0 & G & G \\
    2MASS J18524052+4156057 & - & 3559 & 0.577 & 0.587 & - & K & K \\
    2MASS J18534407+4208274 & - & 3637 & 0.534 & 0.538 & - & K & K \\
    2MASS J18535462+4135227 & - & 3735 & 0.619 & 0.639 & - & K & K \\
    2MASS J19335656+4010546 & - & 3495 & 0.384 & 0.393 & - & M & M \\
    \hline
    \end{tabular}
\end{table}
\end{landscape}