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

In this chapter the selection criteria for the data is explained in 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, 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 was 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 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}. Additionally the dataset of the 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. Due to not all data being available as short cadence data from Kepler/K2, only a subset of this large sample was added.

\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}

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 called $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 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 as well as the minimum distance of datapoints they need to be apart. The minimum distance between points is set to 1 with no minimum required height. 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. For this purpose every datapoint of the normalized flattened lightcurve is subtracted by 1 to move the average from 1 to 0. Additionally the star and end point of the flare are estimated. 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 enough datapoints for the following steps.\\
Afterwards multiple checks are done. The first checking if the 1 datapoint before the peak, and 1 after the peak are above a threshold of 0.003, or if 2 datapoint after the peak after above the same 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. While this eliminates the positive detection of 2 flares in they case of them appearing very shortly after another, it was by visual inspection found to eliminate more false positives. Shortly after another appearing flares are still allowed, if the criteria are met, and theres atleast one more datapoint between the peaks. Then a fit of the flare is generated. The first half of the fit, till the peak, is that of a gaussian function, with the second half being an exponential decay (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 location 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 spike (values of 10 or higher when normalized).

\subsection{Lightcurve folding}

This section will mainly describe how the $getOptimizedFold()$ function works. It takes the normalized lightcurve as well as a fit type 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 boxleastsquares algorithm. Afterwards the 4 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 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 gets 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, it will fold the lightcurve (using the lightkurve $fold$ function) with the period set to the spot modulation, and the epoch set to the first minima in the lightcurve. Then a fit is calculated using the prefered fit method (with a polynomial fit fallback for sine fit prefered and vice versa). If the resulting fit should have 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 now use this as the period for spot modulation. A new fold for the rotational period is then generated. 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 shifts the epoch by the difference of the minimum to 0, which is the center of the folded lightcurves 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 rotational period found, 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 should not break out of the loop in less than 30 tries.

\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 various different bands was known though. The stars which have their magnitudes known atleast in the B and V bands are found in table \ref{tab:unknown_sptypes_bv}, while the stars which had their magnitudes known only in the J, H and C bands are in table \ref{tab:unknown_sptypes_jhc}.
Using $B-V$ (as well as further magnitude differences, if available) have been compared to the table at \cite{spectral_type_color_table} (\href{https://www.stsci.edu/~inr/intrins.html}{https://www.stsci.edu/~inr/intrins.html}). This allowed for a rough classification of the stars in table \ref{tab:unknown_sptypes_bv}. Additionally, and to classify the stars in table \ref{tab:unknown_sptypes_jhc}, 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 base spectral type calculated from the effective temperature and the magnitude differences in multiple different color bands 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 is a G type star based on the effective temperature, the latter is prefered. For the four stars from table \ref{tab:unknown_sptypes_jhc} which did not provide magnitudes in the optical range, the spectral type based on the effective temperature is taken as final.

\begin{table}
    \caption{List of stars with no spectral type entry on 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.}
    \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}
    \caption{List of stars with no spectral type entry on SIMBAD, but available magnitudes in the J, H and C bands. Main Identifier according to SIMBADs "MAIN\_ID" property, as well as TESS Input Catalogue and Kepler Input Catalogue numbers.}
    \label{tab:unknown_sptypes_jhc}
    \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}
    \caption{Final estimated spectral types. Main Identifier shows the MAIN\_ID property for the star on SIMBAD. The column SIMBAD B-V contains the value of B-V with the respective magnitudes taken from SIMBAD. The columns TIC 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 estimated 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 estimated and used in this study.}
    \label{tab:unknown_sptypes_final_output}
    \begin{tabular}{lccccccc}
    \hline
    Main Identifier & SIMBAD B-V & TIC T\textsubscript{eff} & TIC Mass & TIC 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}